loon: PTSimApplication Class Reference

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.) {
00040   // Normal constructor
00041 
00042   MSG("PTSim",Msg::kDebug) << "PTSimApplication normal ctor @ " 
00043                            << this << endl;
00044 
00045   fPTSimStack = new PTSimStack();
00046 
00047 }

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

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

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

00550                                    {    
00551   // User actions after finishing of an event
00552 
00553   MSG("PTSim",Msg::kDebug) << "PTSimApplication::FinishEvent " 
00554                            << fSnarl << "/" << fEvent << endl;
00555 
00556   // Fill stdhep array with all primaries and selected secondaries
00557   // Adjust digihit track id's starting with first hit in this event 
00558   // to reflect stdheparray indices
00559   fPTSimStack -> FillStdHepArray(fStdHepArray,fSnarl,fHitArray,fNHitBegEvt);
00560   fPTSimStack -> Reset();
00561 
00562 }

00541                                      {    
00542   // User actions after finishing of a primary track
00543 
00544   MSG("PTSim",Msg::kDebug) << "PTSimApplication::FinishPrimary " 
00545                        << fSnarl << "/" << fEvent << "/" << fPrimary << endl;
00546 
00547 }

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

00068     { return fEvtStdHepMap; }

00049 { return fNEvents; }

00048 { return fNSnarls; }

00051 { return fRkVFactor; }

00050 { return fRkVMinDistInCM; }

00047 { return fPTSimStack; }

00842                                        {
00843   // Private method to initialize fRockdEdXMin with the minimum dE/dX from 
00844   // LOSS tables for MARS rock material for mu-. Units of fRockdEdXMin are 
00845   // GeV/cm.  This requires the use of TGeant3 specific code to do it right, 
00846   // which is to check the LOSS tables.  For non-TGeant3 concrete MC, will use 
00847   // a nominal dE/dX value of 1.2e-3 GeV/(g/cm^2) * density of the rock. 
00848 
00849 
00850   Double_t nominaldEdXMin = 1.2e-3;  // GeV/(g/cm^2)
00851 
00852   TGeant3* geant3 = dynamic_cast<TGeant3*>(gMC);
00853   if ( !geant3 ) {
00854     MSG("PTSim",Msg::kWarning) 
00855       << "No support for non-TGeant3 concrete VMC in InitRockdEdXMin\n"
00856       << "to determine minimum dE/dX from LOSS tables.  Will use nominal "
00857       << "value of " << nominaldEdXMin <<  " GeV/(g/cm^2) * the density "
00858       << " of the rock." << endl;
00859     fRockdEdXMin = nominaldEdXMin 
00860              * (gGeoManager->GetVolume("MARS")->GetMaterial()->GetDensity());
00861   }
00862   else {
00863     // TGeant3 allows access to LOSS tables, so do it the right way
00864     
00865     // Determine kinetic energy range of interest using gcmulo common block to 
00866     // extract range for which LOSS table has been calculated.  Nominal range 
00867     // is 10^-5 GeV to 10^4 GeV.  
00868     Float_t log10KEmin = log10(geant3 -> Gcmulo() -> ekmin);
00869     log10KEmin = TMath::Max(log10KEmin,(Float_t)-5.); // don't go below 10 keV
00870     Float_t log10KEmax = log10(geant3 -> Gcmulo() -> ekmax);
00871 
00872     // Now prepare to use Gftmat to extract dE/dX
00873 
00874     // Set up 1000 input KE bins (tkin) at which dE/dX will be extracted from 
00875     // table (value).
00876     const Int_t nlog10KEbin = 1000;
00877     Float_t* tkin = new Float_t[nlog10KEbin];
00878     Float_t* value = new Float_t[nlog10KEbin];
00879 
00880     // Initialize output array returned by Gftmat, contains five energy 
00881     // thresholds for material
00882     Float_t pcut[5] = {0};
00883 
00884     // Material index and particle Id, const_cast because 
00885     // TGeoMaterial::GetIndex() isn't const 
00886     TString volName = "MARS";
00887     TString matName;  // material Name
00888     Int_t imate = 0;  // material Index
00889     Double_t a, z, dens, radl, inter;
00890     TArrayD par;
00891     if ( !(geant3 
00892            -> GetMaterial(volName,matName,imate,a,z,dens,radl,inter,par)) ) {
00893       MSG("PTSim",Msg::kError) << "Failed to retrieve material for volume " 
00894                                << volName.Data() << ". Abort." << endl;
00895       abort();
00896     }
00897     MSG("PTSim",Msg::kDebug) 
00898       << "InitRockdEdXMin: GetMaterial for volume MARS retrieved:\n"
00899       << " imate " << imate << " a " << a << " z " << z << " dens " << dens 
00900       << " name " << matName.Data() << endl;
00901       
00902     Int_t ipart = 6;  // mu-
00903     
00904     char mechname[] = "LOSS";
00905 
00906     Float_t de = (log10KEmax - log10KEmin)/(nlog10KEbin);
00907     for ( Int_t ie = 0; ie < nlog10KEbin; ie++ ) {
00908       Float_t log10tkin = log10KEmin + Float_t(ie)*de;
00909       tkin[ie] = pow((Float_t)10., log10tkin);
00910       value[ie] = 0;
00911     }
00912     // Gftmat input:
00913     //    imate - material number
00914     //    ipart - particle number
00915     //    chmeca - name of mechanism bank to be retrieved
00916     //    kdim - dimension of the arrays tkin and value
00917     //    tkin - array of kinetic energies in GeV
00918     // output:
00919     //    value - array of energy losses in MeV/cm or macroscopic cross 
00920     //            sections in 1/cm
00921     //    pcut - array containing the 5 energy thresholds for the materials
00922     //           in GeV (CUTGAM, CUTELE, CUTNEU, CUTHAD, CUTMUO)
00923     //    ixst - return code ( = 1 if filled, = 0 otherwise)
00924     //  
00925     Int_t ixst = 0;
00926     geant3 -> Gftmat(imate,ipart,mechname,nlog10KEbin,tkin,value,pcut,ixst);
00927     
00928     if ( ixst != 1 ) {
00929       MSG("PTSim",Msg::kError) << "Gftmat call to retrieve dE/dX for matId " 
00930                                << imate << " over log10(KE) range "
00931                                << log10KEmin << " to " << log10KEmax 
00932                                << "\nfailed when initializing energy "
00933                                << "threshold cuts.  Abort." << endl;
00934       abort();
00935     }    
00936 
00937     // Now loop over value, looking for min value.  This is min dE/dX
00938     Float_t minValue = 1.e10;
00939     for ( Int_t ie = 0; ie < nlog10KEbin; ie++ ) {
00940       minValue = TMath::Min(value[ie],minValue);
00941     }
00942 
00943     fRockdEdXMin = minValue*1.e-3; // convert to GeV/cm
00944 
00945     delete [] tkin;
00946     delete [] value;
00947     
00948   }
00949   
00950   MSG("PTSim",Msg::kInfo) 
00951     << "Initialize rock minimum dE/dX for application of veto cuts to: "
00952     << fRockdEdXMin << " (GeV/cm)" << endl;
00953 
00954   return;
00955   
00956 }

00566                                                {
00567   //
00568   // Purpose:  Private method used to fill fEvtStdHepMap with events
00569   //           and associated primaries based on information in input
00570   //           stdheparray.  At the end of the method, the stdheparray
00571   //           contents will be erased in preparation for the transport.
00572   //           At the end of transport the stdheparray will be filled
00573   //           with both primaries & secondaries. 
00574   //
00575   // Arguments: snarl: current snarl number
00576   //            stdheparray: filled with primaries and secondaries. Ownership
00577   //                         remains that of the caller.
00578   //            nevent: number of events in snarl (default 1)
00579   //  
00580   // Return:  none.
00581   //
00582 
00583 
00584   // Extract random number seeds to print.
00585   const MCApplication& mcapp = MCApplication::Instance();
00586   std::string evtstr = " event.";
00587   if ( nevent > 1 ) evtstr = " events.";
00588   
00589   MSG("PTSim",Msg::kInfo) << "***** PTSimApplication::InitSnarl: Snarl " 
00590                           << snarl << " with " << nevent << evtstr.c_str() 
00591                           << " Starting seed " 
00592                           << mcapp.GetRandom()->GetSeed() << " *****" << endl;
00593 
00594   if ( fRkVLevel > 0 && fRockdEdXMin < 0. ) InitRockdEdXMin();
00595 
00596   ResetEvtStdHepMap();
00597   fSnarl = snarl; // current snarl
00598   fNSnarls++;  // total number of snarls
00599   fEvent=-1;
00600   
00601   fStdHepArray = stdheparray; // reference to stdheparray, not owned.
00602   
00603   if ( nevent <= 0 ) return; // nothing to build
00604 
00605   Int_t nstdhep = fStdHepArray->GetEntriesFast();
00606 
00607   Int_t mcindex = -1;
00608   Int_t nprim = 0;
00609   Int_t idx0 = 0;
00610   bool  prevchild = true;
00611   std::map<int,int> indexmap;
00612 
00613   for ( Int_t istd = 0; istd < nstdhep; istd++ ) {
00614     TParticle* simstdhep = dynamic_cast<TParticle*>(fStdHepArray->At(istd));
00615 
00616     if ( simstdhep->GetMother(0) == -1 && simstdhep->GetMother(1) == -1 ) {
00617       // primary, check to make sure its new event
00618       nprim++;
00619       if ( nprim > 2 || prevchild ) {
00620         if ( mcindex != -1 ) {
00621 
00622           // Clean up previous event stdhep indices
00623           std::vector<TParticle*>& particlelist = fEvtStdHepMap[mcindex];
00624           Int_t nsize = particlelist.size();
00625           for ( int ip = 0; ip < nsize; ip++ ) {
00626             TParticle* particle = particlelist[ip];
00627             if ( particle->GetFirstMother() != -1 ) {
00628               // All mothers should be in index map
00629               std::map<int,int>::iterator inditr 
00630                                 = indexmap.find(particle->GetFirstMother());
00631               if (inditr == indexmap.end() ) {
00632                 MSG("PTSim",Msg::kError) << "No mother index in map. Abort" 
00633                                          << endl;
00634                 abort();
00635               }
00636               particle->SetFirstMother(inditr->second);
00637             }
00638             if ( particle->GetSecondMother() != -1 ) {
00639               // All mothers should be in index map
00640               std::map<int,int>::iterator inditr 
00641                                 = indexmap.find(particle->GetSecondMother());
00642               if (inditr == indexmap.end() ) {
00643                 MSG("PTSim",Msg::kError) << "No mother index in map. Abort" 
00644                                          << endl;
00645                 abort();
00646               }
00647               particle->SetLastMother(inditr->second);
00648             }
00649             if ( particle->GetFirstDaughter() != -1 ) {
00650               // check for first daughter in map
00651               std::map<int,int>::iterator inditr 
00652                                 = indexmap.find(particle->GetFirstDaughter());
00653               if (inditr != indexmap.end() )
00654                 particle->SetFirstDaughter(inditr->second);
00655               else
00656                 particle->SetFirstDaughter(-1);
00657             }             
00658             if ( particle->GetLastDaughter() != -1 ) {
00659               // check for last daughter in map
00660               std::map<int,int>::iterator inditr 
00661                                = indexmap.find(particle->GetLastDaughter());
00662               if (inditr != indexmap.end() )
00663                 particle->SetLastDaughter(inditr->second);
00664               else
00665                 particle->SetLastDaughter(-1);
00666             }
00667           }
00668           indexmap.clear();
00669         }
00670       
00671         mcindex++;
00672         nprim = 1;
00673         idx0 = istd;
00674         if ( mcindex >= nevent ) {
00675           MSG("PTSim",Msg::kError)
00676              << "PTSimModule::BuildEvtStdHepMap:\nBreakdown in procedure "
00677              << "to match stdhep to associated mc entry.\n"
00678              << "Unable to process snarl.  Abort." << endl;
00679           abort();
00680         }
00681       }
00682     
00683       prevchild = false;
00684     }
00685     
00686     else {
00687       prevchild = true;
00688       nprim = 0;
00689     }
00690     
00691 
00692     Int_t statuscode = simstdhep->GetStatusCode();
00693     bool isprimary = false;
00694     if ( statuscode < 200 || statuscode == 999 ) {
00695       // Primary
00696       isprimary = true;
00697     }
00698     else if ( statuscode == 1205 ) {
00699       // May also be primary if status code of children is < 200
00700       TParticle* child = dynamic_cast<TParticle*>
00701                           (fStdHepArray->At(simstdhep->GetFirstDaughter()));
00702       if ( child->GetStatusCode() < 200 ) isprimary = true;
00703     }
00704     if ( isprimary ) {
00705       // Parent/child indices will be cleaned up at end of this event
00706       // Push primaries only
00707       std::vector<TParticle*>& stdheplist = fEvtStdHepMap[mcindex];
00708       TParticle* stdhepcopy = new TParticle(*simstdhep);
00709       // Convert to cm for use in simulation
00710       stdhepcopy -> SetProductionVertex(stdhepcopy->Vx()*100.,
00711                                         stdhepcopy->Vy()*100.,
00712                                         stdhepcopy->Vz()*100.,
00713                                         stdhepcopy->T());
00714       // Adjust pdg code if ion to be pdg2006 standard if not already done
00715       if ( UtilPDG::isIon(stdhepcopy -> GetPdgCode()) ) {
00716         int pdg2006 = UtilPDG::stdIonNumber(stdhepcopy->GetPdgCode(),
00717                                             UtilPDG::kPDG2006);
00718         stdhepcopy -> SetPdgCode(pdg2006);
00719       }
00720  
00721       indexmap.insert(std::make_pair(istd,stdheplist.size()));
00722       stdheplist.push_back(stdhepcopy);
00723     }
00724     
00725   }
00726 
00727   // At end, check grand total
00728   if ( mcindex != nevent - 1  ) {
00729     MSG("PTSim",Msg::kError)
00730       << "PTSimModule::InitSnarl:"
00731       << "\nBreakdown in procedure to match stdhep to "
00732       << "associated mc entry.\nFound " << mcindex + 1 << " primaries "
00733       << "in stdhep array, but nevent = " << nevent << ".\n"
00734       << "Unable to process snarl.  Abort." << endl;
00735     abort();
00736   }
00737 
00738 
00739   // Clean up last event stdhep indices
00740   std::vector<TParticle*>& particlelist = fEvtStdHepMap[mcindex];
00741   Int_t nsize = particlelist.size();
00742   for ( int ip = 0; ip < nsize; ip++ ) {
00743     TParticle* particle = particlelist[ip];
00744     if ( particle->GetFirstMother() != -1 ) {
00745       // All mothers should be in index map
00746       std::map<int,int>::iterator inditr 
00747                  = indexmap.find(particle->GetFirstMother());
00748       if (inditr == indexmap.end() ) {
00749         MSG("PTSim",Msg::kError) << "No mother index in map. Abort" << endl;
00750         abort();
00751       }
00752       particle->SetFirstMother(inditr->second);
00753     }
00754     if ( particle->GetSecondMother() != -1 ) {
00755       // All mothers should be in index map
00756       std::map<int,int>::iterator inditr 
00757                  = indexmap.find(particle->GetSecondMother());
00758       if (inditr == indexmap.end() ) {
00759         MSG("PTSim",Msg::kError) << "No mother index in map. Abort" << endl;
00760         abort();
00761       }
00762       particle->SetLastMother(inditr->second);
00763     }
00764     if ( particle->GetFirstDaughter() != -1 ) {
00765       // check for daughter in map
00766       std::map<int,int>::iterator inditr 
00767                   = indexmap.find(particle->GetFirstDaughter());
00768       if (inditr != indexmap.end() )
00769         particle->SetFirstDaughter(inditr->second);
00770       else
00771         particle->SetFirstDaughter(-1);
00772     }             
00773     if ( particle->GetLastDaughter() != -1 ) {
00774       // check for daughter in map
00775       std::map<int,int>::iterator inditr 
00776                  = indexmap.find(particle->GetLastDaughter());
00777       if (inditr != indexmap.end() )
00778         particle->SetLastDaughter(inditr->second);
00779       else
00780         particle->SetLastDaughter(-1);
00781      }
00782    }
00783    
00784    indexmap.clear();
00785 
00786    // Erase stdheparray contents, pointers to primary TParticles
00787    // now contained in map
00788    fStdHepArray -> Clear("C");
00789 
00790    return;
00791 
00792 }

00495                                  {    
00496   // Private method called on each step to determine if track is exiting. 
00497   // Determined as track has previously exited liner and is now entering
00498   // mars volume.  The algorithm is taken from R. Hatcher's gminos. 
00499   //
00500 
00501   if ( gMC -> IsTrackExiting() ) {
00502     if ( std::string(gMC -> CurrentVolName() ) == "LINR" ) {
00503       fExitLiner = true;
00504     }
00505   }
00506   else if ( gMC -> IsTrackEntering() && fExitLiner ) {
00507     if ( std::string(gMC -> CurrentVolName() ) == "MARS" ) {
00508       return true;
00509     }
00510   }
00511 
00512   return false;
00513   
00514 }

00517                                   {    
00518   // Private method called on each step to determine if track particle
00519   // type is neutrino.   
00520 
00521   Int_t abspdgId = TMath::Abs(gMC -> TrackPid());
00522 
00523   if ( abspdgId == 12 ) return true;  // nu_e, nu_e_bar
00524   else if ( abspdgId == 14 ) return true;  // nu_mu, nu_mu_bar
00525   else if ( abspdgId == 16 ) return true;  // nu_tau, nu_tau_bar
00526   else if ( abspdgId == 18 ) return true;  // nu'_tau, nu'_tau_bar
00527   
00528   return false;
00529   
00530 }

00811                                     {
00812   // Private method called on each primary to determine if track is in rock,
00813   // is charged, and fails to satisfy minimum energy threshold required to
00814   // reach detector cavern.  Return true if track is "Vetoed."
00815 
00816   bool isRockVetoed = false;
00817   
00818   // No veto'ing for uncharged particles
00819   if ( !(gMC->TrackCharge()) ) return isRockVetoed;
00820 
00821   Int_t copyNo = 0;
00822   if ( gMC -> CurrentVolID(copyNo) 
00823        == gGeoManager->GetMasterVolume()->GetNumber() ) {
00824     // In rock
00825     TGeoNode* node = gGeoManager->GetMasterVolume()->GetNode("LINR_1");
00826     TLorentzVector pos;
00827     gMC -> TrackPosition(pos);
00828     Double_t gxyz[3] = {pos.X(),pos.Y(),pos.Z()};
00829     Double_t safety = node -> Safety(gxyz,false);
00830     if ( safety > fRkVMinDistInCM ) {
00831       TLorentzVector mom;
00832       gMC->TrackMomentum(mom);
00833       if ( safety > fRkVFactor*mom.E()/fRockdEdXMin ) isRockVetoed = true;
00834     }
00835   }
00836   
00837   return isRockVetoed;
00838   
00839 }

00533                                  {    
00534   // User actions after finishing of each track
00535 
00536   MSG("PTSim",Msg::kVerbose) << "PTSimApplication::PostTrack" << endl;
00537 
00538 }

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

00795                                          {
00796 
00797   EvtStdHepMap::reverse_iterator ritr;
00798   for ( ritr = fEvtStdHepMap.rbegin(); ritr != fEvtStdHepMap.rend(); ritr++ ) {
00799     std::vector<TParticle*>& particlelist = ritr->second;
00800     for ( unsigned int ip = 0; ip < particlelist.size(); ip++ ) {
00801       TParticle* particle = particlelist[ip];
00802       if ( particle ) delete particle; particle = 0;
00803     }
00804     particlelist.clear();
00805   }
00806   fEvtStdHepMap.clear();
00807   
00808 }

00054 { fHitArray = array; }

00060 { fRkVFactor = rkvfactor; }

00057 { fRkVLevel = rkvlevel; }

00059                          { fRkVMinDistInCM = rkvmindistincm; }

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

00062 { fMCHit.UseLinearStep(uselinearstep); }


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 PTSimStack *	GetStack () const
virtual Int_t	GetNSnarls () const
virtual Int_t	GetNEvents () const
virtual Double_t	GetRkVMinDistInCM () const
virtual Double_t	GetRkVFactor () 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	UseLinearStep (Bool_t uselinearstep=true)
Protected Member Functions
virtual void	ResetEvtStdHepMap ()
virtual const EvtStdHepMap &	GetEvtStdHepMap () const
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
PTSimStack *	fPTSimStack
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
Friends
class	PTSimValidate

PTSimApplication Class Reference

Public Types

Public Member Functions

Protected Member Functions

Private Member Functions

Private Attributes

Friends

Member Typedef Documentation

Constructor & Destructor Documentation

Member Function Documentation

Friends And Related Function Documentation

Member Data Documentation