GeoSwimApplication.cxx

Go to the documentation of this file.

//
// GeoSwimApplication
//
// January 20, 2006  M.Ishitsuka First version for GeoSwimmer
//                   ref. PTSimApplication.cxx by S. Kasahara 05/04
//

#include <iostream>
#include <fstream>
using namespace std;

#include <TGeoManager.h>
#include <TGeoTrack.h>
#include <TVirtualMC.h>
#include <TParticlePDG.h>
#include <TLorentzVector.h>
#include <TVector3.h>
#include <TDatabasePDG.h>
#include <TClonesArray.h>
#include <TArrayI.h>
#include <TMCProcess.h>
#include <TNtuple.h>
#include <TH1D.h>
#include <TMath.h>

#include "GeoSwimmer/GeoSwimApplication.h"
#include "GeoSwimmer/GeoSwimZCondition.h"
#include "MCApplication/MCApplication.h"
#include "MessageService/MsgService.h"
#include "Conventions/Munits.h"

# define g3ekbin  g3ekbin_
# define gcomad gcomad_
# define type_of_call
# define DEFCHARD     const char* 
# define DEFCHARL   , const int 
# define PASSCHARD(string) string 
# define PASSCHARL(string) , strlen(string) 

extern "C"
{
  void type_of_call g3ekbin();
  void type_of_call gcomad(DEFCHARD, Int_t*& DEFCHARL);
}

CVSID("$Id: GeoSwimApplication.cxx,v 1.8 2008/08/01 04:23:12 ishi Exp $");

ClassImp(GeoSwimApplication)

//......................................................................

//__________________________________________________________
GeoSwimApplication::GeoSwimApplication(const char* name, const char* title) 
  : MCAppUser(name,title),fGeoSwimStack(0),
    fExitLiner(false)
{
  // Normal constructor

  MSG("GeoSwim",Msg::kDebug) << "GeoSwimApplication normal ctor @ " 
                               << this << endl;
  fGeoSwimStack = new GeoSwimStack();

  MSG("GeoSwim",Msg::kDebug) << "Get address of common GCTRAK"  << endl;
  gcomad(PASSCHARD("GCTRAK"),(int*&) fGctrak  PASSCHARL("GCTRAK"));
  MSG("GeoSwim",Msg::kDebug) << "Get address of common GCFLAG"  << endl;
  gcomad(PASSCHARD("GCFLAG"),(int*&) fGcflag  PASSCHARL("GCFLAG"));

}

//__________________________________________________________
GeoSwimApplication::~GeoSwimApplication() {
  // delete all the owned sub-objects

  MSG("GeoSwim",Msg::kDebug) << "GeoSwimApplication dtor @ " 
                               << this << endl;
 
  if ( fGeoSwimStack ) delete fGeoSwimStack; fGeoSwimStack = 0;

}

//_____________________________________________________________________________
void GeoSwimApplication::BeginEvent() {    


  bool toswim = true;
  const_cast<MCApplication&>(MCApplication::Instance()).SwitchMedia(toswim);

  MSG("GeoSwim",Msg::kDebug) << "GeoSwimApplication::BeginEvent " 
                           << endl;

}

//_____________________________________________________________________________
void GeoSwimApplication::BeginPrimary() {    
  gGeoManager->SetTminTmax(999999,-999999); // sec
  
}

//_____________________________________________________________________________
void GeoSwimApplication::PreTrack() {    
  //
  // User actions at beginning of each track
  // Called by TGeant3gu::gutrak() before calling g3track.
  gMC->TrackPosition(fInitPosition);
  gMC->TrackMomentum(fInitMomentum);
  fGcflag->nrndm[0] = 0;
  fGcflag->nrndm[1] = 0;
  fStep = 0;
  fS = 0.;
  fRange = 0.;
}

//_____________________________________________________________________________
void GeoSwimApplication::Stepping() {    
  //
  // User actions at each step
  // Called by TGeant3gu::gustep 

  TLorentzVector pos;
  gMC->TrackPosition(pos);
  TLorentzVector mom;
  gMC->TrackMomentum(mom);
  Float_t a=0.,z=0.,dens=0.,radl=0.,absl=0.;
  Float_t small = 1e-04;
  Int_t curmat;
  curmat = gMC->CurrentMaterial(a,z,dens,radl,absl);

  // Increasing energy in case of backward swimming: E' = E - dE + 2*dE.
  if(fDir==-1){
    fGctrak->gekin = fGctrak->gekin + 2.*fGctrak->destep;
    fGctrak->getot = fGctrak->gekin + gMC->TrackMass();
    fGctrak->vect[6] = sqrt((fGctrak->getot+gMC->TrackMass())*fGctrak->gekin);
    g3ekbin();
  }

  fS += gMC->TrackStep();
  fRange += dens*gMC->TrackStep();

  if((fInitPosition.Z() < fZCondition && pos.Z() >= fZCondition)
     || (fInitPosition.Z() > fZCondition  && pos.Z() <= fZCondition)) {

    MSG("GeoSwim",Msg::kDebug) << "Condition satisfied" << endl;
    fSatisfied = true;

    // Correction of position, energy and direction.
    float frac = (TMath::Abs(pos.Z()-fPreVect[2]) > small) ?
      (fZCondition-fPreVect[2])/(pos.Z()-fPreVect[2]):
      0.;
    // Position
    fGctrak->vect[0] = fPreVect[0]+(fGctrak->vect[0]-fPreVect[0])*frac;
    fGctrak->vect[1] = fPreVect[1]+(fGctrak->vect[1]-fPreVect[1])*frac;
    fGctrak->vect[2] = fPreVect[2]+(fGctrak->vect[2]-fPreVect[2])*frac;
    // Direction
    TVector3 direction(fGctrak->vect[3]*frac + fPreVect[3]*(1.-frac),
                       fGctrak->vect[4]*frac + fPreVect[4]*(1.-frac),
                       fGctrak->vect[5]*frac + fPreVect[5]*(1.-frac));
    direction = direction.Unit();
    fGctrak->vect[3] = direction.x();
    fGctrak->vect[4] = direction.y();
    fGctrak->vect[5] = direction.z();
    // Momentum and energy
    fGctrak->vect[6] = (fGctrak->vect[6]*frac + fPreVect[6]*(1.-frac));
    fGctrak->getot = pow(fGctrak->vect[6],2) + pow(gMC->TrackMass(),2);

    // Strop Run and return "satisfied"
    gMC->StopTrack();
    return;
  }

  if ( IsExiting() ) {
    MSG("GeoSwim",Msg::kDebug) 
      << "Track is entering MARS from LINR, call StopTrack" << endl;
    fSatisfied = false;
    gMC -> StopTrack();  
    return;
  }

  if(fInitMomentum.Z()*mom.Z()<0.){
    MSG("GeoSwim",Msg::kDebug)
      << "Negative Momentum\n";
    fSatisfied = false;
    gMC->StopTrack();
    return;
  }    

  fStep++;
  if(fStep>10000){
    MSG("GeoSwim",Msg::kWarning)
      << "GeoSwimApplication::Stepping more than 10000 steps"
      << fInitPosition.X() << " " << fInitPosition.Y() << " "
      << fInitPosition.Z() << " " << fInitMomentum.X() << " "
      << fInitMomentum.Y() << " " << fInitMomentum.Z() << " "
      << pos.X() << " " << pos.Y() << " " << pos.Z() << " "
      << fZCondition << endl;
    fSatisfied = false;
    gMC->StopTrack();
    return;
  }    

  if (pos.X() > fHallMax.X()/Munits::cm ||
      pos.X() < fHallMin.X()/Munits::cm ||
      pos.Y() > fHallMax.Y()/Munits::cm ||
      pos.Y() < fHallMin.Y()/Munits::cm ||
      pos.Z() > fHallMax.Z()/Munits::cm ||
      pos.Z() < fHallMin.Z()/Munits::cm) {
    MSG("GeoSwim",Msg::kDebug)
      << "GeoSwimmer::Swim Out of hall"
      << pos.X() << " " << pos.Y() << " " << pos.Z() << "\n";
    fSatisfied = false;
    gMC -> StopTrack();  
    return;
  }

  // Store /GCTRAK/VECT(7) as the status of previous step
  for(int i=0;i<7;i++){
      fPreVect[i]=fGctrak->vect[i];
  }
}

//_____________________________________________________________________________
bool GeoSwimApplication::IsExiting() {    
  if ( gMC -> IsTrackExiting() ) {
    if ( std::string(gMC -> CurrentVolName() ) == "LINR" ) {
      fExitLiner = true;
    }
  }
  else if ( gMC -> IsTrackEntering() && fExitLiner ) {
    if ( std::string(gMC -> CurrentVolName() ) == "MARS" ) {
      return true;
    }
  }

  return false;
  
}

//_____________________________________________________________________________
void GeoSwimApplication::PostTrack() {    
  fExitLiner = false; // reset exit flag
  
  MSG("GeoSwim",Msg::kDebug) << "GeoSwimApplication::PostTrack" << endl;

}

//_____________________________________________________________________________
void GeoSwimApplication::FinishPrimary() {    

  MSG("GeoSim",Msg::kDebug) << "GeoSwimApplication::FinishPrimary " << endl;

} 

//_____________________________________________________________________________
void GeoSwimApplication::FinishEvent() {    

  MSG("GeoSim",Msg::kDebug) << "GeoSwimApplication::FinishEvent " << endl;

  fGeoSwimStack -> Reset();

} 

//......................................................................

void GeoSwimApplication::SetZCondition(GeoSwimZCondition &c) 
{ 
  fSatisfied = false;
  fZCondition = c.GetZFinal()/Munits::cm;
}

//......................................................................

void GeoSwimApplication::SetDirection(Int_t idir) 
{ 
  fDir = idir;
}

//......................................................................

void GeoSwimApplication::SetHallSize(TVector3 min, TVector3 max) 
{ 
  fHallMin = min;
  fHallMax = max;
}

//......................................................................

bool GeoSwimApplication::GetSatisfied() 
{ 
  return fSatisfied;
}

//......................................................................

Double_t GeoSwimApplication::GetS() 
{ 
  return fS;
}

//......................................................................

Double_t GeoSwimApplication::GetRange() 
{ 
  return fRange;
}

//_____________________________________________________________________________

---