PTSimApplication Class Reference

#include <PTSimApplication.h>

Inheritance diagram for PTSimApplication:
MCAppUser

List of all members.

Public Types

typedef std::map< Int_t,
std::vector< TParticle * > > 
EvtStdHepMap

Public Member Functions

 PTSimApplication (const char *name, const char *title)
virtual ~PTSimApplication ()
virtual void GeneratePrimaries ()
virtual void BeginEvent ()
virtual void BeginPrimary ()
virtual void PreTrack ()
virtual void Stepping ()
virtual void PostTrack ()
virtual void FinishPrimary ()
virtual void FinishEvent ()
virtual PTSimStackGetStack () const
virtual Int_t GetNSnarls () const
virtual Int_t GetNEvents () const
virtual Double_t GetRkVMinDistInCM () const
virtual Double_t GetRkVFactor () const
virtual Bool_t GetRockOnly () const
virtual Bool_t GetSplitScheme () const
virtual void SetHitArray (TClonesArray *array)
virtual void InitSnarl (Int_t snarl, TClonesArray *stdheparray, Int_t nneukin=1)
virtual void SetRkVLevel (Int_t rkvlevel)
virtual void SetRkVMinDistInCM (Double_t rkvmindistincm)
virtual void SetRkVFactor (Double_t rkvfactor)
virtual void SetRockOnly (Bool_t rkonly)
virtual void SetSplitScheme (Int_t scheme)
virtual void SetPrintSplit (Bool_t flag)
virtual void UseLinearStep (Bool_t uselinearstep=true)
virtual void PrintEvtStdHepMap () const

Protected Member Functions

virtual void ResetEvtStdHepMap ()
virtual const EvtStdHepMapGetEvtStdHepMap () const
virtual void SnarlSplitFull (Int_t nneukin)
virtual void SnarlSplitSimple (Int_t nneukin)

Private Member Functions

bool IsExiting ()
bool IsNeutrino ()
void InitRockdEdXMin ()
bool IsRockVetoed ()

Private Attributes

Int_t fSnarl
Int_t fEvent
Int_t fPrimary
Int_t fNSnarls
Int_t fNEvents
PTSimStackfPTSimStack
TClonesArray * fHitArray
TClonesArray * fStdHepArray
 array of DigiScintHits, not owned
PTSimHit fMCHit
 array of TParticles, not owned
Bool_t fExitLiner
Int_t fLogLevel
Int_t fLastStackSize
Int_t fNHitBegEvt
EvtStdHepMap fEvtStdHepMap
std::vector< Int_t > fSensitiveMedId
Int_t fNSensitiveMedId
std::string fCurrentPath
Int_t fRkVLevel
Double_t fRkVMinDistInCM
Double_t fRkVFactor
Double_t fRockdEdXMin
Bool_t fRockOnly
Int_t fSplitScheme
Bool_t fPrintSplit

Friends

class PTSimValidate

Detailed Description

Definition at line 22 of file PTSimApplication.h.


Member Typedef Documentation

typedef std::map<Int_t,std::vector<TParticle*> > PTSimApplication::EvtStdHepMap

Definition at line 28 of file PTSimApplication.h.


Constructor & Destructor Documentation

PTSimApplication::PTSimApplication ( const char *  name,
const char *  title 
)

Definition at line 34 of file PTSimApplication.cxx.

References Msg::kDebug, and MSG.

00035   : MCAppUser(name,title),fSnarl(-1),fEvent(-1),fPrimary(-1),fNSnarls(0),
00036     fNEvents(0),fPTSimStack(0),fHitArray(0),fStdHepArray(0),fMCHit(),
00037     fExitLiner(false),fLogLevel(kInfo),fLastStackSize(0),fNHitBegEvt(0),
00038     fNSensitiveMedId(0),fCurrentPath(""),fRkVLevel(0),fRkVMinDistInCM(-1.),
00039     fRkVFactor(-1.),fRockdEdXMin(-1.),fRockOnly(false),fSplitScheme(3),
00040     fPrintSplit(false) {
00041   // Normal constructor
00042 
00043   MSG("PTSim",Msg::kDebug) << "PTSimApplication normal ctor @ " 
00044                            << this << endl;
00045 
00046   fPTSimStack = new PTSimStack();
00047 
00048 }

PTSimApplication::~PTSimApplication (  )  [virtual]

Definition at line 51 of file PTSimApplication.cxx.

References fPTSimStack, Msg::kDebug, MSG, and ResetEvtStdHepMap().

00051                                     {
00052   // delete all the owned sub-objects
00053 
00054   MSG("PTSim",Msg::kDebug) << "PTSimApplication dtor @ " 
00055                            << this << endl;
00056  
00057   if ( fPTSimStack ) delete fPTSimStack; fPTSimStack = 0;
00058   ResetEvtStdHepMap();
00059 
00060 }


Member Function Documentation

void PTSimApplication::BeginEvent (  )  [virtual]

Implements MCAppUser.

Definition at line 144 of file PTSimApplication.cxx.

References fEvent, fHitArray, fLastStackSize, fLogLevel, fNEvents, fNHitBegEvt, fNSensitiveMedId, fPrimary, fPTSimStack, fSensitiveMedId, fSnarl, MsgStream::GetLogLevel(), Nav::GetName(), MCAppStack< T >::GetNtrack(), MsgService::GetStream(), gGeoManager, MCApplication::Instance(), MsgService::Instance(), Msg::kDebug, and MSG.

00144                                   {    
00145   // User actions at beginning of event
00146 
00147   fEvent++; // event in snarl
00148   fNEvents++; // total number of events
00149   
00150   MSG("PTSim",Msg::kDebug) << "  *** PTSimApplication::BeginEvent: Snarl " 
00151                            << fSnarl << " Event " << fEvent 
00152                            << " ***  " << endl;
00153   fPrimary = -1;
00154 
00155   bool toswim = false;
00156   const_cast<MCApplication&>(MCApplication::Instance()).SwitchMedia(toswim);
00157 
00158   // Log the sensitive media here.  This is used to improve the
00159   // performance during stepping, since the sensitive media is likely
00160   // to be a small subset of all the media.
00161   fSensitiveMedId.clear();
00162   TIter next(gGeoManager->GetListOfMedia());
00163   TGeoMedium* med = 0;
00164   fNSensitiveMedId = 0;
00165   while ((med = (TGeoMedium*)next())) {
00166     if ( med->GetParam(0) > 0 ) {
00167       // sensitive medium
00168       std::string medName = med -> GetName();
00169       if ( medName.find("_SWIM") == std::string::npos ) {
00170         fSensitiveMedId.push_back(med->GetId());
00171         fNSensitiveMedId++;
00172       }
00173     }
00174   }
00175 
00176   fLogLevel = MsgService::Instance()->GetStream("PTSim")->GetLogLevel();
00177   fLastStackSize = fPTSimStack->GetNtrack();
00178   fNHitBegEvt    = fHitArray->GetEntriesFast();
00179   
00180     
00181 }

void PTSimApplication::BeginPrimary (  )  [virtual]

Implements MCAppUser.

Definition at line 184 of file PTSimApplication.cxx.

References fEvent, fLogLevel, fPrimary, fPTSimStack, fSnarl, MCAppStack< T >::GetCurrentTrack(), MCAppStack< T >::GetCurrentTrackNumber(), Msg::kDebug, and MSG.

00184                                     {    
00185   // User actions at beginning of a primary track
00186 
00187   fPrimary++;
00188   if ( fLogLevel <= Msg::kDebug ) {
00189     const TParticle* primary = fPTSimStack->GetCurrentTrack();
00190     Int_t pdg = primary->GetPdgCode();
00191     std::string pname = primary->GetName();
00192     Int_t trkno = fPTSimStack->GetCurrentTrackNumber();
00193  
00194     MSG("PTSim",Msg::kDebug) << "PTSimApplication::BeginPrimary " 
00195                              << fSnarl << "/" << fEvent << "/"
00196                              << fPrimary <<  " trkno " << trkno << " pdg "
00197                              << pdg << "/" << pname.c_str() <<  endl;
00198   }  
00199   
00200 }

void PTSimApplication::FinishEvent (  )  [virtual]

Implements MCAppUser.

Definition at line 561 of file PTSimApplication.cxx.

References fEvent, fHitArray, fNHitBegEvt, fPTSimStack, fSnarl, fStdHepArray, Msg::kDebug, and MSG.

00561                                    {    
00562   // User actions after finishing of an event
00563 
00564   MSG("PTSim",Msg::kDebug) << "PTSimApplication::FinishEvent " 
00565                            << fSnarl << "/" << fEvent << endl;
00566 
00567   // Fill stdhep array with all primaries and selected secondaries
00568   // Adjust digihit track id's starting with first hit in this event 
00569   // to reflect stdheparray indices
00570   fPTSimStack -> FillStdHepArray(fStdHepArray,fSnarl,fHitArray,fNHitBegEvt);
00571   fPTSimStack -> Reset();
00572 
00573 } 

void PTSimApplication::FinishPrimary (  )  [virtual]

Implements MCAppUser.

Definition at line 552 of file PTSimApplication.cxx.

References fEvent, fPrimary, fSnarl, Msg::kDebug, and MSG.

00552                                      {    
00553   // User actions after finishing of a primary track
00554 
00555   MSG("PTSim",Msg::kDebug) << "PTSimApplication::FinishPrimary " 
00556                        << fSnarl << "/" << fEvent << "/" << fPrimary << endl;
00557 
00558 }

void PTSimApplication::GeneratePrimaries (  )  [virtual]

Implements MCAppUser.

Definition at line 63 of file PTSimApplication.cxx.

References fEvent, fEvtStdHepMap, fLogLevel, fPTSimStack, fSnarl, geant3, MCAppStack< T >::GetNprimary(), MCAppStack< T >::GetNtrack(), MCApplication::GetRandom(), MCApplication::Instance(), Msg::kDebug, Msg::kError, Msg::kWarning, MSG, MCAppStack< T >::Print(), and PTSimStack::Reset().

00063                                          {
00064   // Fill the user stack with primary particles
00065 
00066   MSG("PTSim",Msg::kDebug) << "PTSimApplication::GeneratePrimaries " 
00067                            << fSnarl << "/" << fEvent << endl;
00068 
00069   // Clear stack
00070   fPTSimStack->Reset();
00071   
00072   // Push primaries to stack for this event
00073   Int_t nevent = fEvtStdHepMap.size();
00074   if ( fEvent >= nevent ) {
00075     MSG("PTSim",Msg::kError) << "GeneratePrimaries called with Event "
00076                              << fEvent << " but only " << nevent
00077                              << " particles in stack. Abort." << endl;
00078     abort();
00079   }
00080 
00081   std::vector<TParticle*>& stdheplist = fEvtStdHepMap[fEvent];
00082 
00083   TVector3 pol;
00084   for ( unsigned int is = 0; is < stdheplist.size(); is++ ) {
00085     TParticle* prim = stdheplist[is];
00086 
00087     Int_t toBeDone = 0; // only track final states
00088     if ( prim->GetStatusCode() == 1 ) {
00089       if ( abs(prim->GetPdgCode()) == 311 ) {
00090         // If Monte Carlo type is TGeant3, convert K0/K0bar to K0L/K0S 
00091         // so that Geant3 can handle them
00092         TGeant3* geant3 = dynamic_cast<TGeant3*>(gMC);
00093         if ( geant3 ) {
00094           Double_t randno = MCApplication::Instance().GetRandom() -> Rndm();
00095           if ( randno <= 0.5 ) prim -> SetPdgCode(130); // K0L
00096           else prim -> SetPdgCode(310); // K0S
00097         }
00098       }
00099       Int_t primPdg = prim->GetPdgCode();
00100       if ( primPdg != 28 && primPdg != 0 ) {
00101         // particles 28 (reggeon) and 0 (rootino) unknown to MC
00102         if ( gMC->IdFromPDG(primPdg) <= 0 ) {
00103           // particle unknown to monte carlo and not reggeon or rootino
00104           // print warning
00105           const TDatabasePDG& dbpdg = *(TDatabasePDG::Instance());
00106           std::string pname = "???";
00107           if (dbpdg.GetParticle(primPdg)) 
00108                                 pname = dbpdg.GetParticle(primPdg)->GetName();
00109           MSG("PTSim",Msg::kWarning) << "Particle of type "
00110                                      << primPdg << "/" << pname.c_str()  
00111                                      << " not defined to Monte Carlo."
00112                                      << " Skipped!" << endl;
00113         }
00114         else toBeDone = 1;
00115       }
00116     }
00117 
00118     Int_t parent = prim->GetMother(0); // parent trkid, -1 => no parent track
00119     prim->GetPolarisation(pol);
00120     Int_t trkId = -1; // filled by stack
00121     // kPPrimary is used for creator process VMC code
00122     fPTSimStack -> PushTrack(toBeDone,parent,prim->GetPdgCode(),
00123                              prim->Px(),prim->Py(),prim->Pz(),prim->Energy(),
00124                              prim->Vx(),prim->Vy(),prim->Vz(),prim->T(),
00125                              pol.X(),pol.Y(),pol.Z(),kPPrimary,trkId,
00126                              prim->GetWeight(),prim->GetStatusCode());
00127      
00128   }
00129   
00130   if ( fLogLevel <= Msg::kDebug ) {
00131     MSG("PTSim",Msg::kDebug) << "PTSimApplication::GeneratePrimaries "
00132                              << fSnarl << "/" << fEvent << " pushed "
00133                              << fPTSimStack->GetNtrack() << " particles w/"
00134                              << fPTSimStack->GetNprimary() << " primaries "
00135                              << "to stack." << endl;
00136     fPTSimStack->Print();
00137   }
00138   
00139   return;
00140 
00141 }

virtual const EvtStdHepMap& PTSimApplication::GetEvtStdHepMap (  )  const [inline, protected, virtual]

Definition at line 75 of file PTSimApplication.h.

References fEvtStdHepMap.

Referenced by PTSimValidate::TestInitSnarl().

00076     { return fEvtStdHepMap; }

virtual Int_t PTSimApplication::GetNEvents (  )  const [inline, virtual]

Definition at line 49 of file PTSimApplication.h.

References fNEvents.

00049 { return fNEvents; }

virtual Int_t PTSimApplication::GetNSnarls (  )  const [inline, virtual]

Definition at line 48 of file PTSimApplication.h.

References fNSnarls.

00048 { return fNSnarls; }

virtual Double_t PTSimApplication::GetRkVFactor (  )  const [inline, virtual]

Definition at line 51 of file PTSimApplication.h.

References fRkVFactor.

00051 { return fRkVFactor; }

virtual Double_t PTSimApplication::GetRkVMinDistInCM (  )  const [inline, virtual]

Definition at line 50 of file PTSimApplication.h.

References fRkVMinDistInCM.

00050 { return fRkVMinDistInCM; }

virtual Bool_t PTSimApplication::GetRockOnly (  )  const [inline, virtual]

Definition at line 52 of file PTSimApplication.h.

References fRockOnly.

00052 { return fRockOnly; }

virtual Bool_t PTSimApplication::GetSplitScheme (  )  const [inline, virtual]

Definition at line 53 of file PTSimApplication.h.

References fSplitScheme.

00053 { return fSplitScheme; }

virtual PTSimStack* PTSimApplication::GetStack (  )  const [inline, virtual]

Implements MCAppUser.

Definition at line 47 of file PTSimApplication.h.

References fPTSimStack.

00047 { return fPTSimStack; }

void PTSimApplication::InitRockdEdXMin (  )  [private]

Definition at line 989 of file PTSimApplication.cxx.

References MuELoss::a, fRockdEdXMin, geant3, gGeoManager, Msg::kDebug, Msg::kError, Msg::kInfo, Msg::kWarning, and MSG.

Referenced by InitSnarl().

00989                                        {
00990   // Private method to initialize fRockdEdXMin with the minimum dE/dX from 
00991   // LOSS tables for MARS rock material for mu-. Units of fRockdEdXMin are 
00992   // GeV/cm.  This requires the use of TGeant3 specific code to do it right, 
00993   // which is to check the LOSS tables.  For non-TGeant3 concrete MC, will use 
00994   // a nominal dE/dX value of 1.2e-3 GeV/(g/cm^2) * density of the rock. 
00995 
00996 
00997   Double_t nominaldEdXMin = 1.2e-3;  // GeV/(g/cm^2)
00998 
00999   TGeant3* geant3 = dynamic_cast<TGeant3*>(gMC);
01000   if ( !geant3 ) {
01001     MSG("PTSim",Msg::kWarning) 
01002       << "No support for non-TGeant3 concrete VMC in InitRockdEdXMin\n"
01003       << "to determine minimum dE/dX from LOSS tables.  Will use nominal "
01004       << "value of " << nominaldEdXMin <<  " GeV/(g/cm^2) * the density "
01005       << " of the rock." << endl;
01006     fRockdEdXMin = nominaldEdXMin 
01007              * (gGeoManager->GetVolume("MARS")->GetMaterial()->GetDensity());
01008   }
01009   else {
01010     // TGeant3 allows access to LOSS tables, so do it the right way
01011     
01012     // Determine kinetic energy range of interest using gcmulo common block to 
01013     // extract range for which LOSS table has been calculated.  Nominal range 
01014     // is 10^-5 GeV to 10^4 GeV.  
01015     Float_t log10KEmin = log10(geant3 -> Gcmulo() -> ekmin);
01016     log10KEmin = TMath::Max(log10KEmin,(Float_t)-5.); // don't go below 10 keV
01017     Float_t log10KEmax = log10(geant3 -> Gcmulo() -> ekmax);
01018 
01019     // Now prepare to use Gftmat to extract dE/dX
01020 
01021     // Set up 1000 input KE bins (tkin) at which dE/dX will be extracted from 
01022     // table (value).
01023     const Int_t nlog10KEbin = 1000;
01024     Float_t* tkin = new Float_t[nlog10KEbin];
01025     Float_t* value = new Float_t[nlog10KEbin];
01026 
01027     // Initialize output array returned by Gftmat, contains five energy 
01028     // thresholds for material
01029     Float_t pcut[5] = {0};
01030 
01031     // Material index and particle Id, const_cast because 
01032     // TGeoMaterial::GetIndex() isn't const 
01033     TString volName = "MARS";
01034     TString matName;  // material Name
01035     Int_t imate = 0;  // material Index
01036     Double_t a, z, dens, radl, inter;
01037     TArrayD par;
01038     if ( !(geant3 
01039            -> GetMaterial(volName,matName,imate,a,z,dens,radl,inter,par)) ) {
01040       MSG("PTSim",Msg::kError) << "Failed to retrieve material for volume " 
01041                                << volName.Data() << ". Abort." << endl;
01042       abort();
01043     }
01044     MSG("PTSim",Msg::kDebug) 
01045       << "InitRockdEdXMin: GetMaterial for volume MARS retrieved:\n"
01046       << " imate " << imate << " a " << a << " z " << z << " dens " << dens 
01047       << " name " << matName.Data() << endl;
01048       
01049     Int_t ipart = 6;  // mu-
01050     
01051     char mechname[] = "LOSS";
01052 
01053     Float_t de = (log10KEmax - log10KEmin)/(nlog10KEbin);
01054     for ( Int_t ie = 0; ie < nlog10KEbin; ie++ ) {
01055       Float_t log10tkin = log10KEmin + Float_t(ie)*de;
01056       tkin[ie] = pow((Float_t)10., log10tkin);
01057       value[ie] = 0;
01058     }
01059     // Gftmat input:
01060     //    imate - material number
01061     //    ipart - particle number
01062     //    chmeca - name of mechanism bank to be retrieved
01063     //    kdim - dimension of the arrays tkin and value
01064     //    tkin - array of kinetic energies in GeV
01065     // output:
01066     //    value - array of energy losses in MeV/cm or macroscopic cross 
01067     //            sections in 1/cm
01068     //    pcut - array containing the 5 energy thresholds for the materials
01069     //           in GeV (CUTGAM, CUTELE, CUTNEU, CUTHAD, CUTMUO)
01070     //    ixst - return code ( = 1 if filled, = 0 otherwise)
01071     //  
01072     Int_t ixst = 0;
01073     geant3 -> Gftmat(imate,ipart,mechname,nlog10KEbin,tkin,value,pcut,ixst);
01074     
01075     if ( ixst != 1 ) {
01076       MSG("PTSim",Msg::kError) << "Gftmat call to retrieve dE/dX for matId " 
01077                                << imate << " over log10(KE) range "
01078                                << log10KEmin << " to " << log10KEmax 
01079                                << "\nfailed when initializing energy "
01080                                << "threshold cuts.  Abort." << endl;
01081       abort();
01082     }    
01083 
01084     // Now loop over value, looking for min value.  This is min dE/dX
01085     Float_t minValue = 1.e10;
01086     for ( Int_t ie = 0; ie < nlog10KEbin; ie++ ) {
01087       minValue = TMath::Min(value[ie],minValue);
01088     }
01089 
01090     fRockdEdXMin = minValue*1.e-3; // convert to GeV/cm
01091 
01092     delete [] tkin;
01093     delete [] value;
01094     
01095   }
01096   
01097   MSG("PTSim",Msg::kInfo) 
01098     << "Initialize rock minimum dE/dX for application of veto cuts to: "
01099     << fRockdEdXMin << " (GeV/cm)" << endl;
01100 
01101   return;
01102   
01103 }

void PTSimApplication::InitSnarl ( Int_t  snarl,
TClonesArray *  stdheparray,
Int_t  nneukin = 1 
) [virtual]

Definition at line 576 of file PTSimApplication.cxx.

References fEvent, fNSnarls, fPrintSplit, fRkVLevel, fRockdEdXMin, fSnarl, fSplitScheme, fStdHepArray, MCApplication::GetRandom(), InitRockdEdXMin(), MCApplication::Instance(), Msg::kInfo, Msg::kSynopsis, MAXMSG, MSG, PrintEvtStdHepMap(), ResetEvtStdHepMap(), SnarlSplitFull(), and SnarlSplitSimple().

Referenced by PTSimValidate::TestInitSnarl().

00577                                                {
00578   //
00579   // Purpose:  Private method used to fill fEvtStdHepMap with events
00580   //           and associated primaries based on information in input
00581   //           stdheparray.  At the end of the method, the stdheparray
00582   //           contents will be erased in preparation for the transport.
00583   //           At the end of transport the stdheparray will be filled
00584   //           with both primaries & secondaries. 
00585   //
00586   // Arguments: snarl: current snarl number
00587   //            stdheparray: filled with primaries and secondaries. Ownership
00588   //                         remains that of the caller.
00589   //            nevent: number of events in snarl (default 1)
00590   //  
00591   // Return:  none.
00592   //
00593 
00594 
00595   // Extract random number seeds to print.
00596   const MCApplication& mcapp = MCApplication::Instance();
00597   std::string evtstr = " event.";
00598   if ( nevent > 1 ) evtstr = " events.";
00599   
00600   MSG("PTSim",Msg::kSynopsis) << "***** PTSimApplication::InitSnarl: Snarl " 
00601                           << snarl << " with " << nevent << evtstr.c_str() 
00602                           << " Starting seed " 
00603                           << mcapp.GetRandom()->GetSeed() << " *****" << endl;
00604 
00605   if ( fRkVLevel > 0 && fRockdEdXMin < 0. ) InitRockdEdXMin();
00606 
00607   ResetEvtStdHepMap();
00608   fSnarl = snarl; // current snarl
00609   fNSnarls++;  // total number of snarls
00610   fEvent=-1;
00611   
00612   fStdHepArray = stdheparray; // reference to stdheparray, not owned.
00613   
00614   switch ( TMath::Abs(fSplitScheme) ) {
00615   case  1:  // no splitting: fall through, also handled by Simple splitter 
00616   case  2:  // normal simple splitter
00617     SnarlSplitSimple(nevent);
00618     break;
00619   case  3: 
00620     SnarlSplitFull(nevent);
00621     break;
00622   default:
00623     MAXMSG("PTSim",Msg::kInfo,5) << "InitSnarl unrecognized split scheme "
00624                                  << fSplitScheme << ", use SnarlSplitFull"
00625                                  << endl;
00626     SnarlSplitFull(nevent);
00627     break;
00628   }
00629 
00630   // Erase stdheparray contents, pointers to primary TParticles
00631   // now contained in map
00632   fStdHepArray -> Clear("C");
00633 
00634   if ( fPrintSplit ) PrintEvtStdHepMap();  // generally for split debugging
00635 
00636   return;
00637 }

bool PTSimApplication::IsExiting (  )  [private]

Definition at line 506 of file PTSimApplication.cxx.

References fExitLiner.

Referenced by Stepping().

00506                                  {    
00507   // Private method called on each step to determine if track is exiting. 
00508   // Determined as track has previously exited liner and is now entering
00509   // mars volume.  The algorithm is taken from R. Hatcher's gminos. 
00510   //
00511 
00512   if ( gMC -> IsTrackExiting() ) {
00513     if ( std::string(gMC -> CurrentVolName() ) == "LINR" ) {
00514       fExitLiner = true;
00515     }
00516   }
00517   else if ( gMC -> IsTrackEntering() && fExitLiner ) {
00518     if ( std::string(gMC -> CurrentVolName() ) == "MARS" ) {
00519       return true;
00520     }
00521   }
00522 
00523   return false;
00524   
00525 }

bool PTSimApplication::IsNeutrino (  )  [private]

Definition at line 528 of file PTSimApplication.cxx.

Referenced by Stepping().

00528                                   {    
00529   // Private method called on each step to determine if track particle
00530   // type is neutrino.   
00531 
00532   Int_t abspdgId = TMath::Abs(gMC -> TrackPid());
00533 
00534   if ( abspdgId == 12 ) return true;  // nu_e, nu_e_bar
00535   else if ( abspdgId == 14 ) return true;  // nu_mu, nu_mu_bar
00536   else if ( abspdgId == 16 ) return true;  // nu_tau, nu_tau_bar
00537   else if ( abspdgId == 18 ) return true;  // nu'_tau, nu'_tau_bar
00538   
00539   return false;
00540   
00541 }

bool PTSimApplication::IsRockVetoed (  )  [private]

Definition at line 958 of file PTSimApplication.cxx.

References fRkVFactor, fRkVMinDistInCM, fRockdEdXMin, and gGeoManager.

Referenced by Stepping().

00958                                     {
00959   // Private method called on each primary to determine if track is in rock,
00960   // is charged, and fails to satisfy minimum energy threshold required to
00961   // reach detector cavern.  Return true if track is "Vetoed."
00962 
00963   bool isRockVetoed = false;
00964   
00965   // No veto'ing for uncharged particles
00966   if ( !(gMC->TrackCharge()) ) return isRockVetoed;
00967 
00968   Int_t copyNo = 0;
00969   if ( gMC -> CurrentVolID(copyNo) 
00970        == gGeoManager->GetMasterVolume()->GetNumber() ) {
00971     // In rock
00972     TGeoNode* node = gGeoManager->GetMasterVolume()->GetNode("LINR_1");
00973     TLorentzVector pos;
00974     gMC -> TrackPosition(pos);
00975     Double_t gxyz[3] = {pos.X(),pos.Y(),pos.Z()};
00976     Double_t safety = node -> Safety(gxyz,false);
00977     if ( safety > fRkVMinDistInCM ) {
00978       TLorentzVector mom;
00979       gMC->TrackMomentum(mom);
00980       if ( safety > fRkVFactor*mom.E()/fRockdEdXMin ) isRockVetoed = true;
00981     }
00982   }
00983   
00984   return isRockVetoed;
00985   
00986 }

void PTSimApplication::PostTrack (  )  [virtual]

Implements MCAppUser.

Definition at line 544 of file PTSimApplication.cxx.

References Msg::kVerbose, and MSG.

00544                                  {    
00545   // User actions after finishing of each track
00546 
00547   MSG("PTSim",Msg::kVerbose) << "PTSimApplication::PostTrack" << endl;
00548 
00549 }

void PTSimApplication::PreTrack (  )  [virtual]

Implements MCAppUser.

Definition at line 203 of file PTSimApplication.cxx.

References fEvent, fExitLiner, fLogLevel, fPrimary, fPTSimStack, fSnarl, MCAppStack< T >::GetCurrentParentTrackNumber(), MCAppStack< T >::GetCurrentTrackNumber(), PTSim::GetProcess(), Msg::kDebug, Msg::kVerbose, and MSG.

00203                                 {    
00204   // User actions at beginning of each track
00205 
00206   fExitLiner = false; // reset exit flag
00207 
00208   if ( fLogLevel <= Msg::kDebug ) {
00209     Int_t trkno = fPTSimStack->GetCurrentTrackNumber();
00210     TLorentzVector pos;
00211     gMC->TrackPosition(pos);
00212     TLorentzVector mom;
00213     gMC->TrackMomentum(mom);
00214     Int_t pdg = gMC->TrackPid();
00215     const TDatabasePDG& dbpdg = *(TDatabasePDG::Instance());
00216     std::string pname = "???";
00217     if (dbpdg.GetParticle(pdg)) pname = dbpdg.GetParticle(pdg)->GetName();
00218     Int_t parentno = fPTSimStack->GetCurrentParentTrackNumber();
00219 
00220     if ( fLogLevel == Msg::kDebug ) {
00221       if ((fPTSimStack -> GetParticle(trkno) -> GetProcess()) 
00222                                                          == kPPrimary) {
00223         MSG("PTSim",Msg::kDebug) << "PTSimApplication::PreTrack " 
00224           << fSnarl << "/" << fEvent << "/" << fPrimary << "/" << trkno
00225           << " parent " << parentno   
00226           << "\n   " << pdg << "/" << pname.c_str() 
00227           << " " <<  gMC->CurrentVolName() << " v: " << pos.X() 
00228           << " " << pos.Y() << " " << pos.Z() << " " << pos.T() 
00229           << " p: " << mom.X() << " "  << mom.Y() << " " << mom.Z() 
00230           << " " << mom.E() << endl;
00231       }
00232     }
00233     else {
00234       MSG("PTSim",Msg::kVerbose) << "PTSimApplication::PreTrack " 
00235         << fSnarl << "/" << fEvent << "/" << fPrimary << "/" << trkno
00236         << " parent " << parentno   
00237         << "\n   " << pdg << "/" << pname.c_str() 
00238         << " " <<  gMC->CurrentVolName() << " v: " << pos.X() 
00239         << " " << pos.Y() << " " << pos.Z() << " " << pos.T() 
00240         << " p: " << mom.X() << " "  << mom.Y() << " " << mom.Z() 
00241         << " " << mom.E() << endl;
00242     }
00243   }
00244 
00245 }

void PTSimApplication::PrintEvtStdHepMap (  )  const [virtual]

Definition at line 1106 of file PTSimApplication.cxx.

References fEvtStdHepMap, Msg::kInfo, and MSG.

Referenced by InitSnarl(), SnarlSplitFull(), and SnarlSplitSimple().

01107 {
01108   // Print out how InitSnarl has broken up the particle list
01109   // into sub-events.
01110   //    typedef std::map<Int_t,std::vector<TParticle*> > EvtStdHepMap;
01111 
01112   typedef std::map<Int_t,std::vector<TParticle*> >::const_iterator mitr_t;
01113   mitr_t mapitr = fEvtStdHepMap.begin();
01114   mitr_t mapend = fEvtStdHepMap.end();
01115 
01116   MSG("PTSim",Msg::kInfo) << "PTSimApplication::PrintEvtStdHepMap() " 
01117                           << " " << fEvtStdHepMap.size() << " entries"
01118                           << endl;
01119   for ( ; mapitr != mapend ; ++mapitr ) {
01120     Int_t key = mapitr->first;
01121     const std::vector<TParticle*>& partvec = mapitr->second;
01122     MSG("PTSim",Msg::kInfo) << "event " << setw(2) << key << " " << endl;
01123     for (unsigned int i=0; i < partvec.size(); ++i ) {
01124       TParticle* apart = partvec[i];
01125       if ( ! apart ) {
01126         MSG("PTSim",Msg::kInfo) 
01127           << " [" << setw(3) << i << "] -- no entry --" << endl;
01128         continue;
01129       }
01130       // very condensed ...
01131       MSG("PTSim",Msg::kInfo) 
01132         << " ["  << setw(3)  << i << "]" 
01133         << " "   << setw(4)  << apart->GetStatusCode() 
01134         << " "   << setw(16) << apart->GetName()
01135         << " ("  << setw(2)  << apart->GetMother(0)
01136         << ','   << setw(2)  << apart->GetMother(1)
01137         << ")"
01138         << " ("  << setw(2)  << apart->GetDaughter(0)
01139         << ','   << setw(2)  << apart->GetDaughter(1)
01140         << ")"
01141         << " E=" << apart->Energy()
01142         << " {"  << apart->Vx()
01143         << ","   << apart->Vy()
01144         << ","   << apart->Vz()
01145         << ","   << apart->T()
01146         << "}"
01147         << endl;
01148     }
01149   }
01150 
01151 }

void PTSimApplication::ResetEvtStdHepMap (  )  [protected, virtual]

Definition at line 942 of file PTSimApplication.cxx.

References fEvtStdHepMap.

Referenced by InitSnarl(), and ~PTSimApplication().

00942                                          {
00943 
00944   EvtStdHepMap::reverse_iterator ritr;
00945   for ( ritr = fEvtStdHepMap.rbegin(); ritr != fEvtStdHepMap.rend(); ritr++ ) {
00946     std::vector<TParticle*>& particlelist = ritr->second;
00947     for ( unsigned int ip = 0; ip < particlelist.size(); ip++ ) {
00948       TParticle* particle = particlelist[ip];
00949       if ( particle ) delete particle; particle = 0;
00950     }
00951     particlelist.clear();
00952   }
00953   fEvtStdHepMap.clear();
00954   
00955 }

virtual void PTSimApplication::SetHitArray ( TClonesArray *  array  )  [inline, virtual]

Definition at line 56 of file PTSimApplication.h.

References fHitArray.

00056 { fHitArray = array; }

virtual void PTSimApplication::SetPrintSplit ( Bool_t  flag  )  [inline, virtual]

Definition at line 65 of file PTSimApplication.h.

References fPrintSplit.

00065 { fPrintSplit = flag; }

virtual void PTSimApplication::SetRkVFactor ( Double_t  rkvfactor  )  [inline, virtual]

Definition at line 62 of file PTSimApplication.h.

References fRkVFactor.

00062 { fRkVFactor = rkvfactor; }

virtual void PTSimApplication::SetRkVLevel ( Int_t  rkvlevel  )  [inline, virtual]

Definition at line 59 of file PTSimApplication.h.

References fRkVLevel.

00059 { fRkVLevel = rkvlevel; }

virtual void PTSimApplication::SetRkVMinDistInCM ( Double_t  rkvmindistincm  )  [inline, virtual]

Definition at line 60 of file PTSimApplication.h.

References fRkVMinDistInCM.

00061                          { fRkVMinDistInCM = rkvmindistincm; }

virtual void PTSimApplication::SetRockOnly ( Bool_t  rkonly  )  [inline, virtual]

Definition at line 63 of file PTSimApplication.h.

References fRockOnly.

00063 { fRockOnly = rkonly; }

virtual void PTSimApplication::SetSplitScheme ( Int_t  scheme  )  [inline, virtual]

Definition at line 64 of file PTSimApplication.h.

References fSplitScheme.

00064 { fSplitScheme = scheme; }

void PTSimApplication::SnarlSplitFull ( Int_t  nneukin  )  [protected, virtual]

Definition at line 640 of file PTSimApplication.cxx.

References fEvtStdHepMap, fSplitScheme, fStdHepArray, UtilPDG::isIon(), Msg::kError, Msg::kInfo, UtilPDG::kPDG2006, MSG, PrintEvtStdHepMap(), and UtilPDG::stdIonNumber().

Referenced by InitSnarl().

00640                                                   {
00641   //
00642   // Purpose:  Private method used to fill fEvtStdHepMap with events
00643   //           and associated primaries based on information in 
00644   //           fStdHepArray.
00645   //           Skip check on nevent if fSplitScheme is < 0
00646   //
00647   // Arguments: nevent:  (expected) number of events in snarl
00648   //            
00649   // Return:  none.
00650   //
00651 
00652   if ( nevent <= 0 && fSplitScheme > 0 ) return; // nothing to build
00653 
00654   Int_t nstdhep = fStdHepArray->GetEntriesFast();
00655 
00656   Int_t mcindex = -1;
00657   Int_t nprim = 0;
00658   Int_t idx0 = 0;
00659   bool  prevchild = true;
00660   std::map<int,int> indexmap;
00661 
00662   for ( Int_t istd = 0; istd < nstdhep; istd++ ) {
00663     TParticle* simstdhep = dynamic_cast<TParticle*>(fStdHepArray->At(istd));
00664 
00665     if ( simstdhep->GetMother(0) == -1 && simstdhep->GetMother(1) == -1 ) {
00666       // primary, check to make sure its new event
00667       nprim++;
00668       if ( nprim > 2 || prevchild ) {
00669         if ( mcindex != -1 ) {
00670 
00671           // Clean up previous event stdhep indices
00672           std::vector<TParticle*>& particlelist = fEvtStdHepMap[mcindex];
00673           Int_t nsize = particlelist.size();
00674           for ( int ip = 0; ip < nsize; ip++ ) {
00675             TParticle* particle = particlelist[ip];
00676             if ( particle->GetFirstMother() != -1 ) {
00677               // All mothers should be in index map
00678               std::map<int,int>::iterator inditr 
00679                                 = indexmap.find(particle->GetFirstMother());
00680               if (inditr == indexmap.end() ) {
00681                 MSG("PTSim",Msg::kError) << "No mother index in map. Abort" 
00682                                          << endl;
00683                 abort();
00684               }
00685               particle->SetFirstMother(inditr->second);
00686             }
00687             if ( particle->GetSecondMother() != -1 ) {
00688               // All mothers should be in index map
00689               std::map<int,int>::iterator inditr 
00690                                 = indexmap.find(particle->GetSecondMother());
00691               if (inditr == indexmap.end() ) {
00692                 MSG("PTSim",Msg::kError) << "No mother index in map. Abort" 
00693                                          << endl;
00694                 abort();
00695               }
00696               particle->SetLastMother(inditr->second);
00697             }
00698             if ( particle->GetFirstDaughter() != -1 ) {
00699               // check for first daughter in map
00700               std::map<int,int>::iterator inditr 
00701                                 = indexmap.find(particle->GetFirstDaughter());
00702               if (inditr != indexmap.end() )
00703                 particle->SetFirstDaughter(inditr->second);
00704               else
00705                 particle->SetFirstDaughter(-1);
00706             }             
00707             if ( particle->GetLastDaughter() != -1 ) {
00708               // check for last daughter in map
00709               std::map<int,int>::iterator inditr 
00710                                = indexmap.find(particle->GetLastDaughter());
00711               if (inditr != indexmap.end() )
00712                 particle->SetLastDaughter(inditr->second);
00713               else
00714                 particle->SetLastDaughter(-1);
00715             }
00716           }
00717           indexmap.clear();
00718         }
00719       
00720         mcindex++;
00721         nprim = 1;
00722         idx0 = istd;
00723         if ( mcindex >= nevent && fSplitScheme > 0 ) {
00724           MSG("PTSim",Msg::kInfo)  << flush << endl;
00725           MSG("PTSim",Msg::kError)
00726              << "PTSimModule::BuildEvtStdHepMap:\nBreakdown in procedure "
00727              << "to match stdhep to associated mc entry.\n"
00728              << "Found mcindex " << mcindex << ", nevent " << nevent <<".\n"
00729              << "Unable to process snarl.  Abort." << endl;
00730           PrintEvtStdHepMap();
00731           abort();
00732         }
00733       }
00734     
00735       prevchild = false;
00736     }
00737     
00738     else {
00739       prevchild = true;
00740       nprim = 0;
00741     }
00742     
00743 
00744     Int_t statuscode = simstdhep->GetStatusCode();
00745     bool isprimary = false;
00746     if ( statuscode < 200 || statuscode == 999 ) {
00747       // Primary
00748       isprimary = true;
00749     }
00750     else if ( statuscode == 1205 ) {
00751       // May also be primary if status code of children is < 200
00752       TParticle* child = dynamic_cast<TParticle*>
00753                           (fStdHepArray->At(simstdhep->GetFirstDaughter()));
00754       if ( child->GetStatusCode() < 200 ) isprimary = true;
00755     }
00756     if ( isprimary ) {
00757       // Parent/child indices will be cleaned up at end of this event
00758       // Push primaries only
00759       std::vector<TParticle*>& stdheplist = fEvtStdHepMap[mcindex];
00760       TParticle* stdhepcopy = new TParticle(*simstdhep);
00761       // Convert to cm for use in simulation
00762       stdhepcopy -> SetProductionVertex(stdhepcopy->Vx()*100.,
00763                                         stdhepcopy->Vy()*100.,
00764                                         stdhepcopy->Vz()*100.,
00765                                         stdhepcopy->T());
00766       // Adjust pdg code if ion to be pdg2006 standard if not already done
00767       if ( UtilPDG::isIon(stdhepcopy -> GetPdgCode()) ) {
00768         int pdg2006 = UtilPDG::stdIonNumber(stdhepcopy->GetPdgCode(),
00769                                             UtilPDG::kPDG2006);
00770         stdhepcopy -> SetPdgCode(pdg2006);
00771       }
00772  
00773       indexmap.insert(std::make_pair(istd,stdheplist.size()));
00774       stdheplist.push_back(stdhepcopy);
00775     }
00776     
00777   }
00778 
00779   // At end, check grand total
00780   if ( mcindex != nevent - 1 && fSplitScheme > 0 ) {
00781     MSG("PTSim",Msg::kInfo) << flush << endl;
00782     MSG("PTSim",Msg::kError)
00783       << "PTSimModule::InitSnarl:"
00784       << "\nBreakdown in procedure to match stdhep to "
00785       << "associated mc entry.\nFound " << mcindex + 1 << " primaries "
00786       << "in stdhep array, but nevent = " << nevent << ".\n"
00787       << "Unable to process snarl.  Abort." << endl;
00788     PrintEvtStdHepMap();
00789     abort();
00790   }
00791 
00792   // Clean up last event stdhep indices
00793   std::vector<TParticle*>& particlelist = fEvtStdHepMap[mcindex];
00794   Int_t nsize = particlelist.size();
00795   for ( int ip = 0; ip < nsize; ip++ ) {
00796     TParticle* particle = particlelist[ip];
00797     if ( particle->GetFirstMother() != -1 ) {
00798       // All mothers should be in index map
00799       std::map<int,int>::iterator inditr 
00800                  = indexmap.find(particle->GetFirstMother());
00801       if (inditr == indexmap.end() ) {
00802         MSG("PTSim",Msg::kError) << "No mother index in map. Abort" << endl;
00803         abort();
00804       }
00805       particle->SetFirstMother(inditr->second);
00806     }
00807     if ( particle->GetSecondMother() != -1 ) {
00808       // All mothers should be in index map
00809       std::map<int,int>::iterator inditr 
00810                  = indexmap.find(particle->GetSecondMother());
00811       if (inditr == indexmap.end() ) {
00812         MSG("PTSim",Msg::kError) << "No mother index in map. Abort" << endl;
00813         abort();
00814       }
00815       particle->SetLastMother(inditr->second);
00816     }
00817     if ( particle->GetFirstDaughter() != -1 ) {
00818       // check for daughter in map
00819       std::map<int,int>::iterator inditr 
00820                   = indexmap.find(particle->GetFirstDaughter());
00821       if (inditr != indexmap.end() )
00822         particle->SetFirstDaughter(inditr->second);
00823       else
00824         particle->SetFirstDaughter(-1);
00825     }             
00826     if ( particle->GetLastDaughter() != -1 ) {
00827       // check for daughter in map
00828       std::map<int,int>::iterator inditr 
00829                  = indexmap.find(particle->GetLastDaughter());
00830       if (inditr != indexmap.end() )
00831         particle->SetLastDaughter(inditr->second);
00832       else
00833         particle->SetLastDaughter(-1);
00834     }
00835   }
00836    
00837   indexmap.clear();
00838   
00839   return;
00840 }

void PTSimApplication::SnarlSplitSimple ( Int_t  nneukin  )  [protected, virtual]

Definition at line 843 of file PTSimApplication.cxx.

References fEvtStdHepMap, fSplitScheme, fStdHepArray, UtilPDG::isIon(), Msg::kError, Msg::kInfo, UtilPDG::kPDG2006, MSG, PrintEvtStdHepMap(), and UtilPDG::stdIonNumber().

Referenced by InitSnarl().

00843                                                     {
00844   //
00845   // Purpose:  Private method used to fill fEvtStdHepMap with events
00846   //           and associated primaries based on information in 
00847   //           fStdHepArray.
00848   //
00849   // Arguments: nevent:  (expected) number of events in snarl
00850   //            
00851   // Return:  none.
00852   //
00853   
00854   // Strategy:  Each "event" is a unspecified series of status 0 
00855   // (or 1001, 2001 for Minerva/ArgoNeut) entries followed by some w/ 
00856   // other status values.  Seeing status 0 again triggers new event.
00857   // Need special case for Minerva/ArgoNeut because sometimes they
00858   // don't pass through the particles that initiated the interaction
00859   // but only the interaction (in their detectors) itself.
00860 
00861   // Special cases for fSplitScheme (normally 2 for this method):
00862   //  -2:   split as usual, but no nevent check
00863   //   1:   no splitting, always leave as one event (but need to fill map)
00864   //        (also no nevent check)
00865 
00866   // Make no attempt at weeding out prior tracking results (GMINOS entries
00867   // w/ status 205 or 1205 for instance for decays).  If this is required 
00868   // then it needs to be done upstream of here.  If one is reading in
00869   // events using EventKinematics/Module/HepevtModule then that module
00870   // can be configured to do so.  Otherwise one can use the tools in
00871   // Util/UtilHepevt (as does HepevtModule).
00872 
00873 
00874   const Int_t kStatProgenitor  =    0; // for NEUGEN3/GENIE either nu or tgt
00875   const Int_t kStatMinervaTrk  = 1001; // particles tracked in Minerva/Argo
00876   const Int_t kStatArgoNeutTrk = 2001; //  but not to be tracked in MINOS
00877 
00878   Int_t nstdhep = fStdHepArray->GetEntriesFast();
00879 
00880   int mcindex = 0;
00881   Int_t indx0 = 0;
00882   bool ready_for_trigger = false;
00883 
00884   for ( Int_t istd = 0; istd < nstdhep; istd++ ) {
00885     TParticle* simstdhep = dynamic_cast<TParticle*>(fStdHepArray->At(istd));
00886     Int_t status = simstdhep->GetStatusCode();
00887     if ( status == kStatProgenitor || 
00888          status == kStatMinervaTrk || 
00889          status == kStatArgoNeutTrk ) {
00890       if ( ready_for_trigger && TMath::Abs(fSplitScheme) != 1 ) {
00891         // start new event
00892         ++mcindex;
00893         ready_for_trigger = false;
00894         indx0 = istd;
00895       }
00896     } else {
00897       ready_for_trigger = true;
00898     }
00899     std::vector<TParticle*>& stdheplist = fEvtStdHepMap[mcindex];
00900     TParticle* stdhepcopy = new TParticle(*simstdhep);
00901 
00902     // Convert to cm for use in simulation
00903     stdhepcopy -> SetProductionVertex(stdhepcopy->Vx()*100.,
00904                                       stdhepcopy->Vy()*100.,
00905                                       stdhepcopy->Vz()*100.,
00906                                       stdhepcopy->T());
00907 
00908     // Adjust pdg code if ion to be pdg2006 standard if not already done
00909     if ( UtilPDG::isIon(stdhepcopy->GetPdgCode()) ) {
00910       int pdg2006 = UtilPDG::stdIonNumber(stdhepcopy->GetPdgCode(),
00911                                           UtilPDG::kPDG2006);
00912       stdhepcopy->SetPdgCode(pdg2006);
00913     }
00914 
00915     // adjust indices, split lists are now independently starting at indx=0
00916     for ( int i = 0; i < 2; ++i ) {
00917       if ( stdhepcopy->GetMother(i) != -1 ) 
00918         stdhepcopy->SetMother(i,stdhepcopy->GetMother(i)-indx0);
00919       if ( stdhepcopy->GetDaughter(i) != -1 ) 
00920         stdhepcopy->SetDaughter(i,stdhepcopy->GetDaughter(i)-indx0);
00921     }
00922 
00923     // push it into list we're creating for this event
00924     stdheplist.push_back(stdhepcopy);
00925   }
00926 
00927   if ( fSplitScheme > 0 && mcindex >= nevent ) {
00928     MSG("PTSim",Msg::kInfo)  << flush << endl;
00929     MSG("PTSim",Msg::kError)
00930       << "PTSimModule::SnarlSplitSimple:\nBreakdown in procedure "
00931       << "to match stdhep to associated mc entry.\n"
00932       << "Found mcindex " << mcindex << ", nevent " << nevent <<".\n"
00933       << "Unable to process snarl.  Abort." << endl;
00934     PrintEvtStdHepMap();
00935     abort();
00936   }
00937 
00938   return;
00939 }

void PTSimApplication::Stepping (  )  [virtual]

Implements MCAppUser.

Definition at line 248 of file PTSimApplication.cxx.

References PTSimHit::AddELoss(), PTSimHit::AddStep(), PTSimHit::Clear(), PTSimHit::CreateDigiScintHit(), fCurrentPath, fEvent, fHitArray, fLastStackSize, fLogLevel, fMCHit, fNSensitiveMedId, fPrimary, fPTSimStack, fRkVLevel, fRockOnly, fSensitiveMedId, fSnarl, MCAppStack< T >::GetCurrentParentTrackNumber(), MCAppStack< T >::GetCurrentTrackNumber(), PTSimHit::GetELoss(), MCAppStack< T >::GetNtrack(), MCAppStack< T >::GetParticle(), MCAppParticle::GetProcess(), PTSimStack::GetStdHepHitThr(), PTSimStack::GetStdHepSaveSibling(), PTSimHit::GetTrkIndex(), gGeoManager, PTSimParticle::HasHitAboveThresh(), PTSimHit::Init(), IsExiting(), IsNeutrino(), IsRockVetoed(), Msg::kDebug, Msg::kError, Msg::kFatal, Msg::kVerbose, MSG, and PTSimHit::SetCurrentState().

00248                                 {    
00249   // User actions at each step
00250 
00251 
00252   if ( gMC -> IsNewTrack() ) {
00253     Int_t trkno = fPTSimStack->GetCurrentTrackNumber();
00254     bool isPrimary = false;
00255     if (fPTSimStack->GetParticle(trkno)->GetProcess()==kPPrimary) 
00256       isPrimary = true;
00257     if ( fRockOnly && isPrimary ) {
00258       std::string currentVolName = gMC -> CurrentVolName();
00259       if ( currentVolName != "MARS" && currentVolName != "LINR" ) {
00260         MSG("PTSim",Msg::kDebug) << "Cut primary not in Rock" << std::endl;
00261         gMC -> StopTrack();
00262         return;
00263       }
00264     }
00265     else if ( fRkVLevel > 1 || (fRkVLevel == 1  && isPrimary ) ) {
00266       if ( IsRockVetoed() ) {
00267         if ( fLogLevel <= Msg::kDebug ) {
00268           TLorentzVector pos;
00269           gMC->TrackPosition(pos);
00270           TLorentzVector mom;
00271           gMC->TrackMomentum(mom);
00272           Int_t pdg = gMC->TrackPid();
00273           const TDatabasePDG& dbpdg = *(TDatabasePDG::Instance());
00274           std::string pname = "???";
00275           if (dbpdg.GetParticle(pdg)) pname=dbpdg.GetParticle(pdg)->GetName();
00276           Int_t parentno = fPTSimStack->GetCurrentParentTrackNumber();
00277           MSG("PTSim",Msg::kDebug) 
00278             << "Track is in rock and below required energy, call StopTrack\n" 
00279             << fSnarl << "/" << fEvent << "/" << fPrimary << "/" << trkno
00280             << " parent " << parentno   
00281             << "\n   " << pdg << "/" << pname.c_str() 
00282             << " " <<  gMC->CurrentVolName() << " v: " << pos.X() 
00283             << " " << pos.Y() << " " << pos.Z() << " " << pos.T() 
00284             << " p: " << mom.X() << " "  << mom.Y() << " " << mom.Z() 
00285             << " " << mom.E() << endl;
00286         }
00287         gMC -> StopTrack();
00288         return;
00289       }
00290     }
00291   }
00292 
00293   if ( fLogLevel <= Msg::kVerbose ) {
00294 
00295     TLorentzVector msgpos;
00296     gMC->TrackPosition(msgpos);
00297     TLorentzVector msgmom;
00298     gMC->TrackMomentum(msgmom);
00299 
00300     MSG("PTSim",Msg::kVerbose) << "Stepping:\n"
00301     << "particle Id " << gMC->TrackPid() << " in vol " 
00302     <<  gMC->CurrentVolName() 
00303     << "\npos " << msgpos.X() << " " << msgpos.Y() << " " << msgpos.Z() << " " 
00304     << msgpos.T() << "\nmom " << msgmom.X() << " " << msgmom.Y() << " " 
00305     << msgmom.Z() << " " << msgmom.E() << "\neloss " << gMC->Edep() 
00306     << ", step " << gMC->TrackStep() << endl;
00307 
00308     TArrayI processlist;
00309     Int_t nprocess = gMC->StepProcesses(processlist);
00310     std::string processnamelist = "Step Processes:\n";
00311     for (int ip = 0; ip < nprocess; ip++ ) {
00312       processnamelist +=  std::string(TMCProcessName[processlist[ip]]) + "/";
00313     }
00314     MSG("PTSim",Msg::kVerbose) << processnamelist.c_str() << endl;
00315   
00316     MSG("PTSim",Msg::kVerbose) << "Track Is:" 
00317           << (gMC->IsTrackEntering()    ? "Entering/" : "" )  
00318           << (gMC->IsTrackExiting()     ? "Exiting/" : "" )  
00319           << (gMC->IsTrackOut()         ? "Out/" : "" )  
00320           << (gMC->IsTrackStop()        ? "Stop/" : "" )  
00321           << (gMC->IsNewTrack()         ? "New/" : "" )  
00322           << (gMC->IsTrackInside()      ? "Inside/" : "" )  
00323           << (gMC->IsTrackDisappeared() ? "Disappeared/" : "" )  
00324           << (gMC->IsTrackAlive()       ? "Alive/" : "" )  
00325           << endl;
00326   }
00327   
00328   Int_t nsec = gMC->NSecondaries();
00329   if ( nsec != 0 ) {
00330     if ( fLogLevel <= Msg::kVerbose ) {
00331       std::string secondarystr;
00332       if ( nsec == 1 ) 
00333         secondarystr = " secondary generated in this step.";
00334       else 
00335         secondarystr = " secondaries generated in this step.";
00336       MSG("PTSim",Msg::kVerbose) << nsec << secondarystr.c_str() << endl;
00337     }
00338   
00339     if ( fPTSimStack->GetStdHepSaveSibling() ) {
00340       // Record secondary siblings.
00341       // if we have recorded a decay, increment fLastStackSize by 1
00342       // because the decay parent at point of decay will have been
00343       // pushed to stack
00344       TMCProcess prodprocess = gMC->ProdProcess(0);
00345       if ( prodprocess == kPDecay ) fLastStackSize++;
00346     
00347       if ( nsec != (fPTSimStack->GetNtrack() - fLastStackSize) ) {
00348         // Sanity check, difference between total stack size and
00349         // last step stack size should be equal to number of secondaries
00350         // + the parent particle stored at the decay point.
00351         // If not, shut down.
00352         MSG("PTSim",Msg::kError) 
00353         << "Number of secondaries generated in this step " << nsec 
00354         << "\nnot equal to difference in stack size since last addition."
00355         << "Was: " << fLastStackSize << " Now: " << fPTSimStack->GetNtrack()
00356         << " Abort." << endl;
00357         abort();
00358       }
00359       for ( Int_t isec = 0; isec < nsec; isec++ ) {
00360         Int_t secTrkId = fLastStackSize + isec; 
00361         // Tell secondary track about siblings
00362         PTSimParticle* secParticle = const_cast<PTSimParticle*>
00363                                    (fPTSimStack->GetParticle(secTrkId));
00364         for ( Int_t isec2 = 0; isec2 < nsec; isec2++ ) {
00365           if ( isec2 != isec ) {
00366             Int_t secTrkId2 = fLastStackSize + isec2;
00367             secParticle -> AddSibling(fPTSimStack->GetParticle(secTrkId2));
00368           }
00369         }
00370       }
00371     }
00372   }
00373     
00374   if ( fPTSimStack->GetStdHepSaveSibling() ) 
00375                             fLastStackSize = fPTSimStack->GetNtrack();
00376 
00377   // Flag to stop tracking if entering Mars from Liner
00378   if ( IsExiting() ) {
00379     MSG("PTSim",Msg::kDebug) 
00380       << "Track is entering MARS from LINR, call StopTrack " << endl;
00381     gMC -> StopTrack();  
00382     return;
00383   }
00384   if ( IsNeutrino() ) {
00385     MSG("PTSim",Msg::kDebug)
00386       << "Track is neutrino, call StopTrack " << endl;
00387     gMC -> StopTrack();
00388     return;
00389   }
00390   
00391   if ( !(gMC->TrackCharge()) ) return; // no hits recorded for uncharged
00392 
00393   bool isSensitive = false;
00394   Int_t currentMedId = gMC->CurrentMedium();
00395   for ( int imed = 0; imed < fNSensitiveMedId; imed++ ) {
00396     if ( fSensitiveMedId[imed] == currentMedId ) {
00397       isSensitive = true;
00398       MSG("PTSim",Msg::kVerbose) << "Sensitive volume" << endl;
00399       break;    
00400     }
00401   }
00402   
00403   // Only for sensitive volumes record hits
00404   if ( !isSensitive ) return;
00405 
00406   // Check to see if particle IsTrackStop and yet PTSimHit has not
00407   // been initialized.  This occurs in VMC w/G4 which double steps
00408   // some of the time when track has disappeared, e.g. to Decay or to
00409   // Positron annihilation.  The first call to Stepping will have
00410   // track status IsTrackDisappeared, at which point the PTSimHit
00411   // is closed.  The second call will have IsTrackStop.  The behavior
00412   // is inconsistent, and sometimes Stepping is only called once.
00413   // On the other hand, VMC w/G3 always only calls Stepping once,
00414   // with either status IsTrackDisappeared for decays, positron 
00415   // annhilation, etc. or IsTrackStop for energy below threshold. 
00416   if ( gMC -> IsTrackStop() && fMCHit.GetTrkIndex() < 0 ) return;
00417   
00418   if ( gMC -> IsTrackExiting() ) {
00419     // Don't trust gGeoManager->GetCurrentNode() directly for exiting particles
00420     // because TGeant4 doesn't set this properly before Stepping()
00421     fCurrentPath = gGeoManager->GetPath();
00422     gGeoManager->cd(gMC->CurrentVolPath());
00423   }
00424 
00425   // These should all be scint strip pstyrene insert nodes inside coex
00426   // shell, ask for parent GeoStripNode
00427   GeoStripNode* stp = dynamic_cast<GeoStripNode*>(gGeoManager->GetMother());
00428   if ( !stp ) {
00429     MSG("PTSim",Msg::kFatal) << "Stepping claims to be in sensitive volume "
00430                              << " but not strip node! abort." << endl;
00431     abort();
00432   }
00433   PlexStripEndId seid = stp -> GetSEId();
00434 
00435   // convert global position coordinates to coordinates local to strip
00436   TLorentzVector pos;
00437   gMC->TrackPosition(pos);
00438   TLorentzVector mom;
00439   gMC->TrackMomentum(mom);
00440   Double_t gxyz[3] = {pos.X(),pos.Y(),pos.Z()};
00441   Double_t lxyz[3] = {0};
00442   // Use GeoNode::GlobalToLocal and not TGeoManager::MasterToLocal so 
00443   // that coordinates will be stored in strip local coordinates, even 
00444   // when strip is split...
00445   // Modified S. Kasahara 2006/07/11. The following clause about
00446   // GeoNode::GlobalToLocal altering the state is no longer true.
00447   // Comment out the  "gGeoManager->cd(currentpath)" statement 
00448   // since this is cpu-intensive.
00449   // // but need to be careful because GeoNode::GlobalToLocal will
00450   // //  alter the state (cd to strip node path and not psty insert).
00451   // //  Save current global path and cd back after conversion is done.
00452   //    const char* currentpath = gGeoManager->GetPath();
00453   stp -> GlobalToLocal(gxyz,lxyz);
00454   //   gGeoManager -> cd(currentpath);
00455 
00456   TLorentzVector localpos;
00457   localpos.SetXYZT(lxyz[0],lxyz[1],lxyz[2],pos.T());
00458 
00459   Int_t trackid = fPTSimStack->GetCurrentTrackNumber();
00460   if ( gMC -> IsTrackEntering() ) {
00461     fMCHit.Clear();
00462     fMCHit.Init(seid,gMC->TrackPid(),trackid,pos,localpos,mom);
00463   }
00464   else {
00465     // accumulate summed path length & energy deposition
00466     if ( fMCHit.GetTrkIndex() != trackid ) {
00467       MSG("PTSim",Msg::kError) << "Current track no " << trackid
00468          << " not equal to hit trk index " 
00469          << fMCHit.GetTrkIndex() << ". abort." << endl;
00470       abort();
00471     }
00472     fMCHit.AddELoss(gMC->Edep());
00473     fMCHit.AddStep(gMC->TrackStep());
00474     fMCHit.SetCurrentState(pos,localpos,mom);
00475   }
00476  
00477   if ( gMC -> IsTrackExiting() || gMC -> IsTrackStop()
00478                                || gMC -> IsTrackDisappeared() ) {
00479     if ( fMCHit.GetTrkIndex() != trackid ) {
00480       MSG("PTSim",Msg::kError) << "Current track no " << trackid
00481          << " not equal to hit trk index " 
00482          << fMCHit.GetTrkIndex() << ". abort" << endl;
00483       abort();
00484     }
00485     PTSimParticle* track = const_cast<PTSimParticle*>
00486                              (fPTSimStack->GetParticle(trackid));
00487     if ( !(track->HasHitAboveThresh()) && 
00488          fMCHit.GetELoss() >= fPTSimStack->GetStdHepHitThr() ) {
00489       track -> SetHasHitAboveThresh(true);
00490     }
00491     
00492     fMCHit.CreateDigiScintHit(*fHitArray);
00493     fMCHit.Clear(); // clear for next track (overkill but to be safe)
00494   }
00495 
00496   if ( gMC -> IsTrackExiting() ) {
00497     // Don't trust gGeoManager->GetCurrentNode() directly for exiting particles
00498     // because TGeant4 doesn't set this properly before Stepping()
00499     gGeoManager->cd(fCurrentPath.c_str());
00500   }
00501 
00502 }

virtual void PTSimApplication::UseLinearStep ( Bool_t  uselinearstep = true  )  [inline, virtual]

Definition at line 68 of file PTSimApplication.h.

References fMCHit, and PTSimHit::UseLinearStep().

00068 { fMCHit.UseLinearStep(uselinearstep); }


Friends And Related Function Documentation

friend class PTSimValidate [friend]

Definition at line 24 of file PTSimApplication.h.


Member Data Documentation

std::string PTSimApplication::fCurrentPath [private]

Definition at line 107 of file PTSimApplication.h.

Referenced by Stepping().

Int_t PTSimApplication::fEvent [private]
Bool_t PTSimApplication::fExitLiner [private]

Definition at line 100 of file PTSimApplication.h.

Referenced by IsExiting(), and PreTrack().

TClonesArray* PTSimApplication::fHitArray [private]

Definition at line 97 of file PTSimApplication.h.

Referenced by BeginEvent(), FinishEvent(), SetHitArray(), and Stepping().

Definition at line 102 of file PTSimApplication.h.

Referenced by BeginEvent(), and Stepping().

Int_t PTSimApplication::fLogLevel [private]

Definition at line 101 of file PTSimApplication.h.

Referenced by BeginEvent(), BeginPrimary(), GeneratePrimaries(), PreTrack(), and Stepping().

array of TParticles, not owned

Definition at line 99 of file PTSimApplication.h.

Referenced by Stepping(), and UseLinearStep().

Int_t PTSimApplication::fNEvents [private]

Definition at line 95 of file PTSimApplication.h.

Referenced by BeginEvent(), and GetNEvents().

Definition at line 103 of file PTSimApplication.h.

Referenced by BeginEvent(), and FinishEvent().

Definition at line 106 of file PTSimApplication.h.

Referenced by BeginEvent(), and Stepping().

Int_t PTSimApplication::fNSnarls [private]

Definition at line 94 of file PTSimApplication.h.

Referenced by GetNSnarls(), and InitSnarl().

Int_t PTSimApplication::fPrimary [private]

Definition at line 93 of file PTSimApplication.h.

Referenced by BeginEvent(), BeginPrimary(), FinishPrimary(), PreTrack(), and Stepping().

Definition at line 114 of file PTSimApplication.h.

Referenced by InitSnarl(), and SetPrintSplit().

Double_t PTSimApplication::fRkVFactor [private]

Definition at line 110 of file PTSimApplication.h.

Referenced by GetRkVFactor(), IsRockVetoed(), and SetRkVFactor().

Int_t PTSimApplication::fRkVLevel [private]

Definition at line 108 of file PTSimApplication.h.

Referenced by InitSnarl(), SetRkVLevel(), and Stepping().

Definition at line 109 of file PTSimApplication.h.

Referenced by GetRkVMinDistInCM(), IsRockVetoed(), and SetRkVMinDistInCM().

Double_t PTSimApplication::fRockdEdXMin [private]

Definition at line 111 of file PTSimApplication.h.

Referenced by InitRockdEdXMin(), InitSnarl(), and IsRockVetoed().

Bool_t PTSimApplication::fRockOnly [private]

Definition at line 112 of file PTSimApplication.h.

Referenced by GetRockOnly(), SetRockOnly(), and Stepping().

std::vector<Int_t> PTSimApplication::fSensitiveMedId [private]

Definition at line 105 of file PTSimApplication.h.

Referenced by BeginEvent(), and Stepping().

Int_t PTSimApplication::fSnarl [private]
TClonesArray* PTSimApplication::fStdHepArray [private]

array of DigiScintHits, not owned

Definition at line 98 of file PTSimApplication.h.

Referenced by FinishEvent(), InitSnarl(), SnarlSplitFull(), and SnarlSplitSimple().


The documentation for this class was generated from the following files:

Generated on 2 Nov 2017 for loon by  doxygen 1.6.1