loon: PTSimStack Class Reference

00034                        : fStdHepSaveByRange(false),
00035                            fStdHepSaveSibling(1),
00036                            fStdHepSaveAncestor(1),
00037                            fStdHepSelectMask(PTSim::kMomentum),
00038                            fStdHepHitThr(0.001),
00039                            fStdHepThr(0.15) {
00040   // Default constructor
00041 
00042   MSG("PTSim",Msg::kVerbose) << "PTSimStack def ctor @ " << this << endl;
00043   
00044 }

00047                         {
00048   // Destructor
00049 
00050   MSG("PTSim",Msg::kVerbose) << "PTSimStack dtor @ " << this << endl;
00051   Reset();
00052   
00053 }

00480                                                         {
00481   // Protected method.  Adjust digiscinthits' track id to that of the
00482   // tracks stored in the output stdhep array.
00483 
00484   for ( Int_t ihit = firsthit; ihit < hitarray->GetEntriesFast(); ihit++ ) {
00485     DigiScintHit* hit = dynamic_cast<DigiScintHit*>(hitarray->At(ihit));
00486     Int_t trackid = hit->TrackId();
00487     Int_t stdhepid = GetStdHepIndex(trackid);
00488     // transition from fParticles Id to stdhep index
00489     hit -> SetTrackId(stdhepid); 
00490   }
00491 
00492 }

00224                                                                          {
00225   // Fill stdheparray with all primaries and selected secondaries.
00226   // The snarl number is used to determine if secondaries were requested
00227   // for storage for this snarl.
00228   // Adjust hitarray hit track indices starting with hit "firsthit" to 
00229   // match those in newly filled section of stdheparray.
00230   //
00231   // Note: Neither hitarray or stdheparray are adopted by this method.
00232   //       Ownership belongs to caller.
00233 
00234   MSG("PTSim",Msg::kDebug) << "PTSimStack::FillStdHepArray for snarl "
00235                              << snarl << endl;
00236 
00237   if ( !stdheparray ) {
00238     MSG("PTSim",Msg::kFatal) << "PTSimStack::FillStdHepArray: Null "
00239                              << "stdheparray TClonesArray ptr. Abort." << endl;
00240     abort();
00241   }
00242   if ( !hitarray ) {
00243     MSG("PTSim",Msg::kFatal) << "PTSimStack::FillStdHepArray: Null "
00244                              << "hitarray TClonesArray ptr. Abort." << endl;
00245     abort();
00246   }
00247   
00248   Int_t firsttrk = stdheparray -> GetEntriesFast();
00249   Int_t nstdhep = firsttrk;
00250   Int_t npart = GetNtrack();
00251 
00252   MSG("PTSim",Msg::kDebug) << "PTSimStack::FillStdHepArray.  Start with "
00253                       << nstdhep << " particles in input stdhep array and " 
00254                       << npart << " particles in stack." << endl;
00255   
00256   // Determine if secondaries are to be stored for this snarl
00257   bool savesecond = SaveSecondariesForSnarl(snarl);
00258 
00259   Int_t decayparentcode = UtilIstHEP::kMDecay 
00260                         + UtilIstHEP::kProdMethodOffset + 1000;
00261 
00262   std::stack<const PTSimParticle*> ministack;
00263   // First go through and tag all particles that are to be saved to
00264   // secondary array
00265   for ( Int_t ip = 0; ip < npart; ip++ ) {
00266     // Determine if particle is primary or passes user-defined selection cuts
00267     PTSimParticle* particle = const_cast<PTSimParticle*>(GetParticle(ip));
00268     
00269     if ( particle -> GetProcess() == kPPrimary ) {
00270       // primary
00271 
00272       fPartToStdHep.insert(make_pair(particle->GetID(),nstdhep));
00273       TParticle* stdhep = 
00274        new((*stdheparray)[nstdhep++])TParticle(*(particle->GetTParticle()));
00275       // Convert vertex to meters
00276       stdhep -> SetProductionVertex(stdhep->Vx()*Munits::cm,
00277                                     stdhep->Vy()*Munits::cm,
00278                                     stdhep->Vz()*Munits::cm,stdhep->T());
00279 
00280       // Adjust indices of parents/children
00281       if ( particle->GetParentID(0) != -1 ) 
00282         stdhep->SetFirstMother(particle->GetParentID(0)+firsttrk);
00283       if ( particle->GetParentID(1) != -1 ) 
00284         stdhep->SetLastMother(particle->GetParentID(1)+firsttrk);
00285       
00286       if ( particle -> GetNChildren() > 0 ) {
00287         if ( particle->GetStatusCode() == 1 ) {
00288           ministack.push(particle); // prime stack for pushing secondaries
00289         }
00290         else {
00291           Int_t nchild = particle->GetNChildren();
00292           if ( particle->GetChild(nchild-1)->GetStatusCode() == 999 ) {
00293             // last child is rootino
00294             nchild += -1;
00295           }
00296           if (  nchild > 0 ) {
00297             stdhep->SetFirstDaughter(particle->GetChildID(0)+firsttrk);
00298             stdhep->SetLastDaughter(particle->GetChildID(nchild-1)+firsttrk);
00299           }
00300         }
00301       }      
00302     }
00303     else if ( savesecond ) {
00304       // Never save decay parents directly unless selected by decay products
00305       if ( decayparentcode != particle->GetStatusCode() ) {
00306         if ( ParticleIsSelected(particle) ) {
00307           // Save secondary.  Mark ancestors to save, and optionally
00308           // siblings and siblings of direct ancestors to save as well.
00309           particle -> SetToSave(true,fStdHepSaveAncestor,fStdHepSaveSibling);
00310         }
00311       }
00312     }
00313   }
00314     
00315   while ( !ministack.empty() ) {
00316     const PTSimParticle* particle = ministack.top();
00317     ministack.pop();  
00318     
00319     // Loop over children, push to stdhep array so that they are in sequence
00320     // because VMC w/Geant4 doesn't guarantee this.
00321     // Also push children with children to ministack for further processing
00322     
00323     // Secondary set to save should have mother already saved
00324     Int_t motherId = particle->GetID();
00325     std::map<int,int>::iterator idItr = fPartToStdHep.find(motherId);
00326     if (idItr == fPartToStdHep.end() ) {
00327       if (idItr == fPartToStdHep.end() ) {
00328         MSG("PTSim",Msg::kError) << "PTSimStack::FillStdHepArray motherId " 
00329                  << motherId << " not in map. abort." << endl;
00330         abort();
00331       }
00332     }
00333     Int_t motherstdindex = idItr->second;
00334 
00335     Int_t nchildren = particle->GetNChildren();
00336     for ( Int_t ic = 0; ic < nchildren; ic++ ) {
00337       const PTSimParticle* child = dynamic_cast<const PTSimParticle*>
00338                                              (particle->GetChild(ic));
00339 
00340       if ( child -> IsSetToSave() ) {
00341         // Particle is to be filled in stdhep array
00342         TParticle* stdhep 
00343           = new((*stdheparray)[nstdhep])TParticle(*(child->GetTParticle()));
00344           
00345         fPartToStdHep.insert(make_pair(child->GetID(),nstdhep));
00346         stdhep -> SetFirstMother(TMath::Abs(motherstdindex)); 
00347         stdhep -> SetLastMother(TMath::Abs(motherstdindex));
00348         // Convert vertex to meters
00349         stdhep -> SetProductionVertex(stdhep->Vx()*Munits::cm,
00350                                       stdhep->Vy()*Munits::cm,
00351                                       stdhep->Vz()*Munits::cm,stdhep->T());
00352         TParticle* ancestor = dynamic_cast<TParticle*>
00353                               (stdheparray->At(TMath::Abs(motherstdindex)));
00354         if(ancestor->GetFirstDaughter()< 0)
00355                                   ancestor->SetFirstDaughter(nstdhep);
00356         ancestor -> SetLastDaughter(nstdhep);
00357         nstdhep++;
00358       }
00359       else {
00360         // Required because digihit may not have direct parent stored
00361         // in output stdhep array
00362         fPartToStdHep.insert(make_pair(child->GetID(),
00363                                       -TMath::Abs(motherstdindex)));
00364       }
00365 
00366       // push particle to stack if it has children
00367       if ( child->GetNChildren() > 0 ) {
00368         ministack.push(child);
00369       }
00370 
00371     }
00372   }
00373   
00374   // Now adjust hit track indices
00375   AdjustHitTrackId(hitarray,firsthit);
00376   
00377   Int_t nfinal = stdheparray->GetEntriesFast();
00378   MSG("PTSim",Msg::kDebug) << "FillStdHepArray.  Finish with "
00379                            << nfinal << " particles in stdheparray." << endl;
00380 
00381  
00382 }

00048 { return fStdHepHitThr; }

00495                                                 {
00496   // Protected method.
00497   // Return stdhep array index of particle corresponding to 
00498   // fParticles id.  The map fPartToStdHep is created during FillStdHepArray
00499 
00500   std::map<int,int>::const_iterator idItr = fPartToStdHep.find(id);
00501   if (idItr == fPartToStdHep.end() ) {
00502     MSG("PTSim",Msg::kError) << "PTSimStack::GetStdHepIndex for trk Id " << id
00503                              << " not in map. abort." << endl;
00504     abort();
00505   }
00506   return idItr->second;
00507 
00508 }

00044 { return fStdHepSave; }

00046 { return fStdHepSaveAncestor; }

00043 { return fStdHepSaveByRange; }

00045 { return fStdHepSaveSibling; }

00047 { return fStdHepSelectMask; }

00049 { return fStdHepThr; }

00174                                                               {
00175   // Return momentum threshold for particle produced by production
00176   // mechanism process, pdg code pdgId, and parentId.   
00177   // The momentum threshold is determined by looking at the 
00178   // user configured momentum thresholds and taking the minimum
00179   // of the thresholds of all those that match the user's type.
00180   // A "match" to a user-configured threshold by type is made if:
00181   // 1)The process is an exact match or kPNoProcess has been specified in
00182   //   user configured momentum threshold by type, and:
00183   // 2)The pdgId is an exact match or kRootino has been specified in
00184   //   user configured momentum threshold by type, and:
00185   // 3)The parentId is an exact match (can be kRootino)
00186   // If no match, -1 is returned.
00187   //
00188 
00189   Double_t pthreshbytype = -1.;
00190   
00191   PTSimStdHepType thrtype(process,pdgId,parentId);
00192   std::map<PTSimStdHepType,double>::const_iterator threshMapItr;
00193   threshMapItr = fStdHepThrByTypeMap.find(thrtype);
00194   if ( threshMapItr != fStdHepThrByTypeMap.end() ) {
00195     // exact match
00196     pthreshbytype = threshMapItr -> second;
00197   }
00198   
00199   // Try again with particle undefined (kRootino) & process
00200   thrtype = PTSimStdHepType(process,kRootino,parentId);
00201   threshMapItr = fStdHepThrByTypeMap.find(thrtype);
00202   if ( threshMapItr != fStdHepThrByTypeMap.end() ) {
00203     if ( pthreshbytype < 0. ) pthreshbytype = threshMapItr -> second;
00204     else pthreshbytype = TMath::Min(pthreshbytype,threshMapItr->second);
00205   }
00206 
00207   // And again with particle & process undefined (kPNoProcess)
00208   thrtype = PTSimStdHepType(kPNoProcess,pdgId,parentId);
00209   threshMapItr = fStdHepThrByTypeMap.find(thrtype);
00210   if ( threshMapItr != fStdHepThrByTypeMap.end() ) {
00211     if ( pthreshbytype < 0. ) pthreshbytype = threshMapItr -> second;
00212     // If particle passed two different particle type definitions
00213     // using generic particle type for first and a generic production 
00214     // mechanism for second, take minimum of two thresholds.
00215     else pthreshbytype = TMath::Min(pthreshbytype,threshMapItr->second);
00216   }
00217 
00218   return pthreshbytype;
00219   
00220 }

00418                                                                        {
00419   // Protected method.
00420   // According to StdHepXXX configuration parameters, pass or fail
00421   // input particle.
00422   // The criteria for passing a particle are:
00423   // 1) Primaries are always passed.
00424   // 2) Secondaries are passed according to the bits in the fStdHepSelectMask.
00425   //    a) If (fStdHepSelectMask & kMomentum), all secondaries with 
00426   //       momentum above threshold are stored. The momentum threshold is 
00427   //       determined using:
00428   //       i) A momentum threshold appropriate for the secondary particle type/
00429   //          production process (TMCProcess) that resulted in the generation 
00430   //          of the secondary, if specified by the user, OR
00431   //      ii) A default momentum threshold (fStdHepThr) applied to all
00432   //          secondaries generated by processes not covered by i).
00433   //    b) If (fStdHepSelectMask & kHit), all secondaries 
00434   //       resulting in energy deposition hit above threshold are stored.
00435   //       The hit threshold is determined using:
00436   //       i) An energy threshold (fStdHepHitThr).
00437   //    The user may specify any combination of the 2 options and an OR
00438   //    will be taken of the results of the individual options to determine
00439   //    the final pass/fail.
00440 
00441   // Primaries always pass
00442   if ( particle -> GetProcess() == kPPrimary ) return true;
00443   
00444   // Secondaries
00445   if ( fStdHepSelectMask & PTSim::kHit ) {
00446     if ( particle -> HasHitAboveThresh() ) return true;
00447   }
00448   
00449   if ( fStdHepSelectMask & PTSim::kMomentum ) {
00450 
00451     // Pass if above momentum threshold
00452     // Look for a match to StdHepThrByType specification. If multiple matches,
00453     // take minimum threshold.  If no match, take default momentum threshold
00454     // specified by StdHepThr.
00455     Int_t pdgId = particle -> GetPdgCode();
00456     TMCProcess process = particle -> GetProcess();
00457     // try once with no defined parent
00458     Double_t pthresh = GetStdHepThrByType(process,pdgId,kRootino);
00459     // Loop over parents 
00460     for ( UInt_t ip = 0; ip < particle->GetNParents(); ip++ ) {
00461       Int_t parentId = particle->GetParent(ip)->GetPdgCode();
00462       Double_t pthreshbytype = GetStdHepThrByType(process,pdgId,parentId);
00463       if ( pthreshbytype >= 0 ) {
00464         if ( pthresh >= 0 ) pthresh = TMath::Min(pthresh,pthreshbytype);
00465         else pthresh = pthreshbytype;
00466       }
00467     }
00468     if ( pthresh < 0 ) pthresh = GetStdHepThr(); // take def thresh if no match
00469     if ( particle -> GetTParticle() -> P() >= pthresh ) return true; 
00470   }
00471     
00472   // Default is no secondaries stored
00473   return false;
00474 
00475 }

00118                                          {
00119   // Specialized MCAppStack<T>::PushTrack(...) to handle case of storing
00120   // decay parent at point of decay.  Needs to be pushed before decay
00121   // children to keep things in the right order.
00122   // See MCAppStack<T>::PushTrack(...) for full list of arguments.
00123   // 
00124   // ntr - track number of this particle in fParticles, output to caller
00125   // 
00126 
00127   Int_t mechstatus = status;
00128   if ( mech != kPPrimary) mechstatus = (Int_t)UtilIstHEP::GetProdMethod(mech)
00129                                      + UtilIstHEP::kProdMethodOffset;
00130   
00131   if ( mech == kPDecay ) {
00132     // Check last daughter position of parent.  
00133     const PTSimParticle* parent = GetParticle(parentId);
00134     const PTSimParticle* lastchild = dynamic_cast<const PTSimParticle*>
00135                             (parent->GetChild(parent->GetNChildren()-1));
00136     Int_t decayparentId = -1;
00137     Int_t decayparentcode = UtilIstHEP::kMDecay 
00138                           + UtilIstHEP::kProdMethodOffset + 1000;
00139     if ( lastchild != 0 ) {
00140       // Check to see if parent at decay point has already been pushed
00141       if ( lastchild -> GetStatusCode() == decayparentcode )
00142                                       decayparentId = lastchild->GetID();
00143     }
00144     if ( decayparentId < 0 ) {
00145       // Push parent at point of decay now
00146       TLorentzVector pos;
00147       gMC->TrackPosition(pos);
00148       TLorentzVector mom;
00149       gMC->TrackMomentum(mom);
00150       Int_t parentpdg = gMC->TrackPid();
00151       // First argument 0 means already tracked (don't push to stack)
00152       MCAppStack<PTSimParticle>::PushTrack(0,parentId,parentpdg,
00153           mom.X(),mom.Y(),mom.Z(),mom.E(),pos.X(),pos.Y(),pos.Z(),pos.T(),
00154           0.,0.,0.,kPTransportation,decayparentId,1.,decayparentcode);
00155     }
00156     // Now push the decay product with the decayparentId as parent
00157     MCAppStack<PTSimParticle>::PushTrack(toBeDone,decayparentId,pdg,
00158                                         px,py,pz,e,vx,vy,vz,tof,
00159                                         polx,poly,polz,mech,ntr,weight,
00160                                         mechstatus);
00161   }
00162   else {
00163     // All other cases (not decay)
00164     MCAppStack<PTSimParticle>::PushTrack(toBeDone,parentId,pdg,
00165                                         px,py,pz,e,vx,vy,vz,tof,
00166                                         polx,poly,polz,mech,ntr,weight,
00167                                         mechstatus);
00168   }
00169   
00170 }

00385                        {
00386   MSG("PTSim",Msg::kVerbose) << "PTSimStack::Reset " << endl;
00387   
00388   MCAppStack<PTSimParticle>::Reset();
00389   fPartToStdHep.clear();
00390   
00391 }

00395                                                           {
00396   // Protected method.  
00397   // Used to determine if secondaries should be saved for input argument 
00398   // snarl given secondary save configuration parameters.
00399   
00400   bool savesecond = false;
00401   if ( !fStdHepSaveByRange ) {
00402     for ( unsigned int is = 0; is < fStdHepSave.size(); is++ ) {
00403       if ( snarl == fStdHepSave[is] ) { savesecond = true; break; }
00404     }
00405   }
00406   else {
00407     if ( snarl <= fStdHepSave[1] && snarl >= fStdHepSave[0] ) {
00408       savesecond = true;
00409     }
00410   }
00411 
00412   return savesecond;
00413   
00414 }

00061 { fStdHepHitThr = ethresh; }

00055     { fStdHepSave = snarllist; fStdHepSaveByRange = isrange; }

00059                              { fStdHepSaveAncestor = saveancestor; }

00057                              { fStdHepSaveSibling = savesib; }

00060 { fStdHepSelectMask = mask; }

00062 { fStdHepThr = pthresh; }

00057                                                                   {
00058   // Sets the pthresh(GeV/c) threshold over which secondaries produced by the
00059   // particle type will be stored to the output stdhep array.
00060   // Particle type is defined as: pdgId*1000 + process.
00061   // Special cases:
00062   // 1)If pdgId=kRootino, the threshold is applied to all particles generated
00063   //   by the specified production mechanism process.  
00064   // 2)If the production mechanism is process=kPNoProcess, the threshold 
00065   //   is applied to all particles of the specified pdgId, regardless of 
00066   //   the production mechanism.
00067   // For the specified particle type, the specified pthresh overrides the
00068   // generic momentum threshold defined by fStdHepThr.
00069 
00070   // Particle Name
00071   std::string particlename = "all particle types";
00072   if ( pdgId != kRootino ) {
00073     const TDatabasePDG& dbpdg = *(TDatabasePDG::Instance());
00074     particlename = "???";
00075     if ( dbpdg.GetParticle(pdgId) ) particlename 
00076                    = dbpdg.GetParticle(pdgId)->GetName();
00077   }
00078   
00079   // Production Process Name
00080   std::string processname = "all production processes";
00081   if ( process != kPNoProcess ) {
00082     processname = TMCProcessName[process];
00083   }
00084   
00085   if ( process == kPNoProcess ) {
00086     MSG("PTSim",Msg::kDebug) << "StdHepThrByType "
00087                              << pthresh << "(GeV/c)" 
00088                              << " for particle type " << pdgId << "("
00089                              << particlename.c_str() << ")" << endl;
00090   }
00091   else if ( pdgId == kRootino ) {
00092     MSG("PTSim",Msg::kDebug) << "StdHepThrByType "
00093                              << pthresh << "(GeV/c)" 
00094                              << " for production process "
00095                              << process << "(" << processname.c_str() << ")"
00096                              << endl;
00097   }
00098   else {                   
00099     MSG("PTSim",Msg::kDebug) << "StdHepThrByType "
00100                              << pthresh << "(GeV/c)" 
00101                              << " for particle type " << pdgId << "("
00102                              << particlename.c_str() << ")" 
00103                              << " and production process "
00104                              << process << "(" << processname.c_str() << ")"
00105                              << endl;
00106   }
00107 
00108   fStdHepThrByTypeMap[PTSimStdHepType(process,pdgId,parentId)] = pthresh;
00109 
00110 }


Public Member Functions
	PTSimStack ()
virtual	~PTSimStack ()
virtual Bool_t	GetStdHepSaveByRange () const
	State testing methods to retrieve secondary save settings set by user.
virtual const std::vector< int > &	GetStdHepSave () const
virtual Int_t	GetStdHepSaveSibling () const
virtual Int_t	GetStdHepSaveAncestor () const
virtual Int_t	GetStdHepSelectMask () const
virtual Double_t	GetStdHepHitThr () const
virtual Double_t	GetStdHepThr () const
virtual Double_t	GetStdHepThrByType (TMCProcess process, Int_t pdgId, Int_t parentId=kRootino) const
virtual void	SetStdHepSave (std::vector< Int_t > snarllist, bool isrange)
virtual void	SetStdHepSaveSibling (Int_t savesib=1)
virtual void	SetStdHepSaveAncestor (Int_t saveancestor=1)
virtual void	SetStdHepSelectMask (Int_t mask)
virtual void	SetStdHepHitThr (double ethresh)
virtual void	SetStdHepThr (double pthresh)
virtual void	SetStdHepThrByType (double pthresh, TMCProcess process, Int_t pdgId, Int_t parentId=kRootino)
virtual void	PushTrack (Int_t toBeDone, Int_t parent, Int_t pdg, Double_t px, Double_t py, Double_t pz, Double_t e, Double_t vx, Double_t vy, Double_t vz, Double_t tof, Double_t polx, Double_t poly, Double_t polz, TMCProcess mech, Int_t &ntr, Double_t weight, Int_t status)
virtual void	FillStdHepArray (TClonesArray stdheparray, Int_t snarl, TClonesArray hitarray, Int_t firsthit)
virtual void	Reset ()
Protected Member Functions
bool	ParticleIsSelected (const PTSimParticle *particle) const
virtual Int_t	GetStdHepIndex (UInt_t id) const
virtual void	AdjustHitTrackId (TClonesArray *digihitarray, Int_t beginhit) const
bool	SaveSecondariesForSnarl (Int_t snarl) const
Private Attributes
Bool_t	fStdHepSaveByRange
std::vector< int >	fStdHepSave
Int_t	fStdHepSaveSibling
Int_t	fStdHepSaveAncestor
Int_t	fStdHepSelectMask
Double_t	fStdHepHitThr
Double_t	fStdHepThr
std::map< PTSimStdHepType, double >	fStdHepThrByTypeMap
std::map< int, int >	fPartToStdHep
Friends
class	PTSimValidate

PTSimStack Class Reference

Public Member Functions

Protected Member Functions

Private Attributes

Friends

Constructor & Destructor Documentation

Member Function Documentation

Friends And Related Function Documentation

Member Data Documentation