PTSimModule Class Reference

#include <PTSimModule.h>

Inheritance diagram for PTSimModule:
JobCModule

List of all members.

Public Types

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

Public Member Functions

 PTSimModule ()
 ~PTSimModule ()
const RegistryDefaultConfig () const
void Config (const Registry &r)
void HandleCommand (JobCommand *cmd)
void BeginJob ()
void EndJob ()
JobCResult Reco (MomNavigator *mom)
std::string GetMCConfig () const
Int_t GetRandomSeed () const

Private Attributes

std::string fMCConfig
std::string fDecayConfig
PTSimApplicationfMCAppUser
Int_t fRandomSeed
TStopwatch fStopwatch
bool fFailEmptyHitArray

Friends

class PTSimValidate

Detailed Description

Definition at line 26 of file PTSimModule.h.


Member Typedef Documentation

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

Definition at line 32 of file PTSimModule.h.


Constructor & Destructor Documentation

PTSimModule::PTSimModule (  ) 

Definition at line 44 of file PTSimModule.cxx.

References Msg::kFatal, LoadMinosPDG(), and MSG.

00044                         : fMCConfig("<null>"),fDecayConfig("<null>"),
00045                             fMCAppUser(0),fRandomSeed(0),
00046                             fFailEmptyHitArray(false) {
00047   // Default constructor
00048 
00049   LoadMinosPDG(); // Load PDG database w/MINOS specific particles
00050   
00051   fStopwatch.Reset();
00052   fStopwatch.Stop();
00053   
00054 
00055   // The application must be instantiated before the mc implementation
00056   fMCAppUser = new PTSimApplication("PTSimApplication","Minos MC Application");
00057   if ( !fMCAppUser ) {
00058     MSG("PTSim",Msg::kFatal) << "No MC application instance." << endl;
00059     abort();
00060   }
00061 
00062 
00063 }

PTSimModule::~PTSimModule (  ) 

Definition at line 66 of file PTSimModule.cxx.

References MCApplication::Clear(), fMCAppUser, MCApplication::Instance(), Msg::kDebug, and MSG.

00066                           {
00067    // Destructor
00068 
00069   MSG("PTSim",Msg::kDebug) << "PTSimModule dtor @ " << this << endl;
00070   
00071   if ( fMCAppUser ) delete fMCAppUser; fMCAppUser = 0;
00072 
00073   // Clear contents of MCApplication because MCApplication singleton
00074   // destruction is out of order with G4 singleton destruction
00075   MCApplication& mcapp = const_cast<MCApplication&>(MCApplication::Instance());
00076   mcapp.Clear();
00077   
00078 }


Member Function Documentation

void PTSimModule::BeginJob ( void   )  [virtual]

Implement for notification of begin of job

Reimplemented from JobCModule.

Definition at line 187 of file PTSimModule.cxx.

References Msg::kDebug, and MSG.

00187                            {
00188   //
00189   //  Purpose:  Called once at the beginning of the job to
00190   //            handle job module initialization.
00191   //
00192   //  Arguments: none.
00193   //  
00194   //  Return: status.
00195   // 
00196 
00197    MSG("PTSim",Msg::kDebug) << "PTSimModule::BeginJob" << endl;
00198 
00199    // mc implementation initialization is saved until first event
00200    // since vldc is required for geometry construction
00201    
00202 }

void PTSimModule::Config ( const Registry r  )  [virtual]

Return the actual configuration. If your module directly pulls its configuration from the fConfig Registry, you don't need to override this. Override if you have local config variables.

Reimplemented from JobCModule.

Definition at line 150 of file PTSimModule.cxx.

References fDecayConfig, fFailEmptyHitArray, fMCAppUser, fMCConfig, fRandomSeed, Registry::Get(), Msg::kDebug, and MSG.

00150                                           {
00151   //
00152   // Purpose: Configure the module given a registry.
00153   //
00154   // Arguments: Registry to use to configure the module.
00155   //
00156   // Return: none.
00157   //
00158 
00159   MSG("PTSim",Msg::kDebug) << "Config PTSimModule with r=" << r << endl;
00160   
00161   Int_t tmpi;
00162   Double_t tmpd;
00163   const Char_t* tmpc;
00164   if ( r.Get("MCConfig",tmpc) ) fMCConfig = tmpc;
00165   if ( r.Get("DecayConfig",tmpc) ) fDecayConfig = tmpc;
00166 
00167   PTSimStack* stack = fMCAppUser -> GetStack();
00168   if ( r.Get("StdHepThr",tmpd) ) stack -> SetStdHepThr(tmpd);
00169   if ( r.Get("StdHepHitThr",tmpd) ) stack -> SetStdHepHitThr(tmpd);
00170   if ( r.Get("StdHepSelectMask",tmpi) ) stack -> SetStdHepSelectMask(tmpi);
00171   if ( r.Get("StdHepSaveSibling",tmpi) ) stack -> SetStdHepSaveSibling(tmpi);
00172   if ( r.Get("StdHepSaveAncestor",tmpi) ) stack -> SetStdHepSaveAncestor(tmpi);
00173   if ( r.Get("RandomSeed",tmpi) ) fRandomSeed = tmpi;
00174   if ( r.Get("RkVLevel",tmpi) ) fMCAppUser -> SetRkVLevel(tmpi);
00175   if ( r.Get("RkVMinDist",tmpd) ) fMCAppUser -> SetRkVMinDistInCM(tmpd*100.);
00176   if ( r.Get("RkVFactor",tmpd) ) fMCAppUser -> SetRkVFactor(tmpd);
00177   if ( r.Get("UseLinearStep",tmpi) ) fMCAppUser -> UseLinearStep(tmpi);
00178   if ( r.Get("RockOnly",tmpi) ) fMCAppUser -> SetRockOnly(tmpi);
00179   if ( r.Get("FailEmptyHitArray",tmpi) ) 
00180     fFailEmptyHitArray = (tmpi!=0) ? true : false; 
00181   if ( r.Get("SnarlSplitScheme",tmpi) ) fMCAppUser -> SetSplitScheme(tmpi);
00182   if ( r.Get("PrintSplit",tmpi) ) fMCAppUser -> SetPrintSplit(tmpi);
00183 
00184 }

const Registry & PTSimModule::DefaultConfig ( void   )  const [virtual]

Get the default configuration registry. This should normally be overridden. One useful idiom is to implement it like:

const Registry& MyModule::DefaultConfig() const { static Registry cfg; // never is destroyed if (cfg.Size()) return cfg; // already filled it // set defaults: cfg.Set("TheAnswer",42); cfg.Set("Units","unknown"); return cfg; }

Reimplemented from JobCModule.

Definition at line 81 of file PTSimModule.cxx.

References fMCAppUser, fRandomSeed, JobCModule::GetName(), Msg::kDebug, Msg::kFatal, PTSim::kMomentum, Registry::LockValues(), MSG, Registry::Set(), and Registry::UnLockValues().

00081                                                  {
00082   //
00083   // Purpose: Method to return default configuration.
00084   // 
00085   // Arguments: none.
00086   //
00087   // Return: Registry containing default configuration
00088   //
00089 
00090   MSG("PTSim",Msg::kDebug) << "Loading default config" << endl;
00091 
00092   static Registry r; 
00093   std::string name = this->JobCModule::GetName();
00094   name += ".config.default";
00095   r.SetName(name.c_str());
00096 
00097   r.UnLockValues();
00098   r.Set("MCConfig","<null>");
00099   r.Set("DecayConfig","<null>");  // or PTSim_DecayConfig.C
00100   r.Set("StdHepSelectMask",PTSim::kMomentum);
00101   // stdhep select mask default is to save
00102   // all tracks above momentum threshold (StdHepThr) 
00103   r.Set("StdHepThr",0.15);// default p(GeV/c) thresh applied to stdhep trks
00104   // if user has selected mode to save secondaries resulting in energy
00105   // deposition hit, StdHepHitThr determines threshold above which to
00106   // store
00107   r.Set("StdHepHitThr",0.001);// default E(GeV) hit threshold
00108   // flag used to determine if selected secondary siblings are to be saved 
00109   // (0 => false, 1 => true, 2 (default) => true and save siblings of direct 
00110   // ancestors as well if StdHepSaveAncestor true)
00111   r.Set("StdHepSaveSibling",2);// default is to save selected secondary sibs
00112                                // and sibs of ancestors
00113   r.Set("StdHepSaveAncestor",1);// default is to save selected secondary ances
00114   r.Set("RandomSeed",fRandomSeed);  // optional random seed supplied by user
00115   r.Set("RkVLevel",0);  // default is to impose no rock veto cuts
00116   r.Set("RkVMinDist",2.); // default minimum distance at which veto cut applies
00117   // safety factor multiplied in front of E/(dE/dX)_min when determining max 
00118   // distance particle will travel 
00119   r.Set("RkVFactor",1.25); 
00120   r.Set("UseLinearStep",false); // bad idea to set true, but useful for validation against gminos
00121   r.Set("RockOnly",0);  // if set true, primaries not in rock will be skipped 
00122   r.Set("FailEmptyHitArray",false); // does Reco() return JobCResult::kFailed when there are no resulting hits
00123   r.Set("SnarlSplitScheme",3); // event splitting: 3=full, 2=simplified snarl (much more lax about rules), 1=leave as one event (no splitting)
00124   r.Set("PrintSplit",false); // don't print split event by default
00125 
00126   r.LockValues();
00127 
00128   // Set defaults for storing secondaries
00129   // Always store secondary muons, regardless of production mechanism or 
00130   // momentum
00131   PTSimStack* stack = fMCAppUser -> GetStack();
00132   if ( !stack ) {
00133     MSG("PTSim",Msg::kFatal) << "No stack!" << endl;
00134     abort();
00135   }
00136 
00137   // Always store muons at production
00138   stack -> SetStdHepThrByType(0.,kPNoProcess,-13);
00139   stack -> SetStdHepThrByType(0.,kPNoProcess,+13);
00140   // Always store secondaries produced by muon decay, regardless of type
00141   stack -> SetStdHepThrByType(0.,kPDecay,kRootino,-13);
00142   stack -> SetStdHepThrByType(0.,kPDecay,kRootino,+13);
00143 
00144 
00145   return r;
00146 }

void PTSimModule::EndJob (  )  [virtual]

Implement for notification of end of job

Reimplemented from JobCModule.

Definition at line 205 of file PTSimModule.cxx.

References fMCAppUser, fStopwatch, Msg::kDebug, Msg::kInfo, and MSG.

00205                          {
00206   //
00207   //  Purpose:  Called once at the end of the job to handle job
00208   //            module validation.
00209   //
00210   //  Arguments: none.
00211   //  
00212   //  Return: status.
00213   // 
00214 
00215    MSG("PTSim",Msg::kDebug) << "PTSimModule::EndJob" << endl;
00216 
00217    //ProfilerStop();
00218    
00219    fStopwatch.Stop();
00220    Int_t nSnarl = 0;
00221    Int_t nEvent = 0;
00222    if ( fMCAppUser ) {
00223      nSnarl = fMCAppUser -> GetNSnarls();
00224      nEvent = fMCAppUser -> GetNEvents();
00225    }
00226    
00227    MSG("PTSim",Msg::kInfo) << "PTSimModule::EndJob, Time(sec), Real "
00228                            << fStopwatch.RealTime() << ", CPU "
00229                            << fStopwatch.CpuTime() << " # snarls/events " 
00230                            << nSnarl << "/" << nEvent <<  endl;
00231    
00232 }

std::string PTSimModule::GetMCConfig (  )  const [inline]

Definition at line 47 of file PTSimModule.h.

References fMCConfig.

00047 { return fMCConfig; }

Int_t PTSimModule::GetRandomSeed (  )  const [inline]

Definition at line 48 of file PTSimModule.h.

References fRandomSeed.

00048 { return fRandomSeed; }

void PTSimModule::HandleCommand ( JobCommand command  )  [virtual]

Implement to handle a JobCommand

Reimplemented from JobCModule.

Definition at line 235 of file PTSimModule.cxx.

References bfld::AsString(), fMCAppUser, UtilIstHEP::GetTMCProcess(), JobCommand::HaveCmd(), JobCommand::HaveOpt(), Msg::kDebug, Msg::kError, Msg::kInfo, UtilIstHEP::kMNoProcess, UtilIstHEP::kNProdMethod, Msg::kWarning, MSG, MSGSTREAM, JobCommand::PopCmd(), JobCommand::PopFloatOpt(), and JobCommand::PopIntOpt().

00235                                                {
00236   // Used to handle commands that don't work Config
00237   //
00238   
00239   if ( !cmd->HaveCmd() ) return;
00240   string s = cmd->PopCmd();
00241   if ( s == "StdHepThrByType" ) {
00242     UtilIstHEP::EProdMethod prodmethod = UtilIstHEP::kMNoProcess;
00243     Int_t pdgId = kRootino;
00244     Int_t parentId = kRootino;
00245     double thresh = 0.;
00246     if ( cmd->HaveOpt() ) {
00247       thresh = cmd->PopFloatOpt();
00248 
00249       if ( cmd->HaveOpt() ) {
00250         int prodmethodint = cmd -> PopIntOpt();
00251         if ( prodmethodint < 0 || prodmethodint >= UtilIstHEP::kNProdMethod ){
00252           MSG("PTSim",Msg::kWarning) 
00253           << "StdHepThrByType called w/invalid UtilIstHEP::EProdMethod "
00254           << prodmethodint << ". Ignored." << endl;
00255           return;
00256         }
00257         else prodmethod = (UtilIstHEP::EProdMethod)(prodmethodint);
00258 
00259         if ( cmd->HaveOpt() ) {
00260           pdgId = cmd -> PopIntOpt();
00261           if ( cmd->HaveOpt() ) parentId = cmd -> PopIntOpt();
00262         }
00263       }
00264       else {
00265         MSG("PTSim",Msg::kWarning)
00266           << "StdHepThrByType called w/no production method. Ignored." 
00267           << endl;
00268         return;
00269       }  
00270     }
00271     else {
00272       MSG("PTSim",Msg::kWarning)
00273         << "StdHepThrByType called w/no momentum threshold. Ignored."
00274         << endl;
00275       return;
00276     }
00277     
00278     // Particle Name
00279     std::string particlename = " for particles of type ";
00280     if ( pdgId != kRootino ) {
00281       particlename += UtilString::ToString<int>(pdgId) + "/";
00282       const TDatabasePDG& dbpdg = *(TDatabasePDG::Instance());
00283       particlename += "???";
00284       if ( dbpdg.GetParticle(pdgId) ) particlename 
00285                        += dbpdg.GetParticle(pdgId)->GetName();
00286     }
00287     else particlename = "";
00288 
00289     std::string parentname = " for particles w/parent type ";
00290     if ( !particlename.empty() ) parentname = " and parent type ";
00291     if ( parentId != kRootino ) {
00292       parentname += UtilString::ToString<int>(parentId) + "/";
00293       const TDatabasePDG& dbpdg = *(TDatabasePDG::Instance());
00294       parentname += "???";
00295       if ( dbpdg.GetParticle(parentId) ) parentname 
00296                        += dbpdg.GetParticle(parentId)->GetName();
00297     }
00298     else parentname = "";
00299     
00300     if ( prodmethod == UtilIstHEP::kMNoProcess ) {
00301       MSG("PTSim",Msg::kInfo) << "PTSimModule StdHepThrByType " 
00302                               << "configured to " << thresh << "(GeV/c)"
00303                               << particlename.c_str()
00304                               << parentname.c_str()
00305                               << endl;
00306     }
00307     else {
00308       std::string prodmethodname = " for particles produced by ";
00309       if ( !parentname.empty() || !particlename.empty() ) 
00310         prodmethodname = " and produced by ";
00311       MSG("PTSim",Msg::kInfo) << "PTSimModule StdHepThrByType " 
00312                               << "configured to " << thresh << "(GeV/c)"
00313                               << particlename.c_str()
00314                               << parentname.c_str()
00315                               << prodmethodname.c_str() << prodmethod 
00316                               << "/" << UtilIstHEP::AsString(prodmethod) 
00317                               << endl;
00318     }
00319     
00320     PTSimStack* stack = fMCAppUser -> GetStack();
00321     TMCProcess tmcprocess = UtilIstHEP::GetTMCProcess(prodmethod);
00322     stack -> SetStdHepThrByType(thresh,tmcprocess,pdgId,parentId);
00323     return;
00324   }
00325   else if ( s == "StdHepSave" ) {
00326 
00327     std::vector<int> snarllist;
00328     MsgStream& msgDebug = MSGSTREAM("PTSim",Msg::kDebug);
00329     msgDebug << "StdHepSave ";
00330 
00331     bool isrange = false;
00332     Int_t nsnarl = 0;
00333     while ( cmd->HaveOpt() ) {
00334       Int_t snarlno = cmd->PopIntOpt();
00335       msgDebug << snarlno << " ";
00336       snarllist.push_back(TMath::Abs(snarlno));
00337       nsnarl++;
00338       if ( snarlno < 0 ) isrange = true;
00339     }
00340     msgDebug << endl;
00341     if ( isrange && nsnarl != 2 ) {
00342       MSG("PTSim",Msg::kWarning)
00343         << "StdHepSave by range requires two arguments. Ignored." << endl;
00344       snarllist.clear();
00345       isrange = false;
00346     }
00347     
00348     PTSimStack* stack = fMCAppUser -> GetStack();
00349     stack -> SetStdHepSave(snarllist,isrange);
00350 
00351     return;
00352   }
00353 
00354   // Errors fall through to here
00355   MSG("PTSim",Msg::kError) << "Illegal command: " << s.c_str() << endl;
00356   
00357 }

JobCResult PTSimModule::Reco ( MomNavigator mom  )  [virtual]

Implement this for read-write access to the MomNavigator

Reimplemented from JobCModule.

Definition at line 360 of file PTSimModule.cxx.

References RecJobHistory::CreateJobRecord(), fDecayConfig, fFailEmptyHitArray, fMCAppUser, fMCConfig, fRandomSeed, fStopwatch, MomNavigator::GetFragment(), RecRecordImp< T >::GetJobHistory(), MCApplication::GetMC(), MCApplication::GetRandom(), DataUtil::GetRun(), SimSnarlRecord::GetSimSnarlHeader(), GetVldContext(), header, MCApplication::InitMC(), MCApplication::InitSnarl(), MCApplication::Instance(), Msg::kDebug, Msg::kError, JobCResult::kFailed, JobCResult::kPassed, RecJobHistory::kPTSim, Msg::kWarning, MAXMSG, MSG, run(), JobCResult::SetFailed(), MCApplication::SetUserApplication(), and JobCResult::SetWarning().

00360                                               {
00361   //
00362   //  Purpose:  Run the requested mc implementation.
00363   //
00364   //  Arguments: mom.
00365   //  
00366   //  Return: status.
00367   // 
00368 
00369   JobCResult result(JobCResult::kPassed);  
00370   MSG("PTSim",Msg::kDebug) << "PTSimModule::Reco" << endl;
00371 
00372   // Check that mom exists.
00373   assert(mom);
00374 
00375   SimSnarlRecord* record 
00376   = dynamic_cast<SimSnarlRecord*>(mom->GetFragment("SimSnarlRecord"));
00377   if ( !record ) {
00378     MSG("PTSim",Msg::kError) << "No SimSnarlRecord in Mom." << endl;
00379     result.SetWarning().SetFailed();
00380     return result;
00381   }
00382   
00383   RecJobHistory& jobhist 
00384               = const_cast<RecJobHistory&>(record->GetJobHistory());
00385   jobhist.CreateJobRecord(RecJobHistory::kPTSim);
00386   
00387   const SimSnarlHeader* header= record->GetSimSnarlHeader();
00388 
00389   MCApplication& mcapp = const_cast<MCApplication&>(MCApplication::Instance());
00390   const VldContext& vldc = *(record -> GetVldContext());
00391   if ( !(mcapp.GetMC() ) ) {
00392     // VldContext will be used to construct appropriate geometry only once
00393     mcapp.InitMC(fMCConfig.c_str(),vldc,fDecayConfig.c_str());
00394     //ProfilerStart("ptsimprofile");
00395   }
00396   mcapp.InitSnarl(vldc);
00397   
00398   // Tell the MCApplication Instance about the current User application
00399   mcapp.SetUserApplication(fMCAppUser);
00400   
00401   // Extract StdHep array from record
00402   // Particle stack will be initialized using StdHep array contents.
00403   // The particles in this stack represent the initial interaction state 
00404   // particles + secondaries generated during gminos tracking if a reroot 
00405   // file is used as input
00406   TClonesArray* stdheparray = dynamic_cast<TClonesArray*>(const_cast<TObject*>
00407     (record -> FindComponent("TClonesArray","StdHep")));
00408  
00409   if ( !stdheparray || stdheparray->GetEntriesFast() == 0 ) {
00410     MAXMSG("PTSim",Msg::kWarning,20) 
00411     << "No StdHep TClonesArray in SimSnarlRecord or is empty." 
00412     << endl
00413     << "### return kFailed rather than kError"
00414     << endl;
00415     return JobCResult::kFailed;
00416     //rwh//return JobCResult::kError;
00417   } 
00418 
00419   // Extract NeuKinList array from record
00420   // If present, this will be used to determine the number of events
00421   // in the snarl, which may be more than 1.
00422   const TClonesArray* neukinarray = dynamic_cast<const TClonesArray*>
00423     (record -> FindComponent("TClonesArray","NeuKinList"));
00424   Int_t nevent = 1;
00425   if ( neukinarray && neukinarray->GetEntriesFast() != 0 ) 
00426                                 nevent = neukinarray->GetEntriesFast();
00427 
00428   // Initialize output DigiScintHit stack, results from tracking will
00429   // be stored here.
00430   TClonesArray* hitarray = dynamic_cast<TClonesArray*>(const_cast<TObject*>
00431     (record -> FindComponent("TClonesArray","DigiScintHits")));
00432   if ( hitarray ) hitarray -> Clear("C");
00433   else {
00434     hitarray = new TClonesArray("DigiScintHit");
00435     hitarray -> SetName("DigiScintHits");
00436     record -> AdoptComponent(hitarray);
00437   }
00438 
00439   // Set random number seeds for this event(snarl) according to the
00440   // run/snarl number with an optional user specified random seed thrown
00441   // in as in N. Tagg's DigiPEtoRawDigitModule.
00442   Int_t run = header -> GetRun();
00443   Int_t snarl = header -> GetSnarl();
00444   Int_t seed = run+snarl+fRandomSeed;
00445   MSG("PTSim",Msg::kDebug) << "Random number seed for run,snarl("
00446        << run << "," << snarl << ") = " << seed << endl;
00447   mcapp.GetRandom() -> SetSeed(seed);
00448   
00449   fMCAppUser -> SetHitArray(hitarray); // pass the place to store the hits
00450   // Prime snarl using snarl number, contents of stdheparray, and nevent
00451   fMCAppUser -> InitSnarl(snarl,stdheparray,nevent); 
00452 
00453   fStopwatch.Start(false);
00454   // Monte carlo generation takes place here
00455   gMC -> ProcessRun(nevent);
00456   fStopwatch.Stop();
00457 
00458   if ( fFailEmptyHitArray && hitarray->GetLast() < 0 ) {
00459     // JobCResult set to kFailed if requested to reject events with 
00460     // no hits and there are no hits.
00461     // This is generally useful for rock interactions
00462     result.SetFailed();
00463   }
00464 
00465   return result;
00466 
00467 }


Friends And Related Function Documentation

friend class PTSimValidate [friend]

Definition at line 28 of file PTSimModule.h.


Member Data Documentation

std::string PTSimModule::fDecayConfig [private]

Definition at line 54 of file PTSimModule.h.

Referenced by Config(), and Reco().

Definition at line 58 of file PTSimModule.h.

Referenced by Config(), and Reco().

Definition at line 55 of file PTSimModule.h.

Referenced by Config(), DefaultConfig(), EndJob(), HandleCommand(), Reco(), and ~PTSimModule().

std::string PTSimModule::fMCConfig [private]

Definition at line 53 of file PTSimModule.h.

Referenced by Config(), GetMCConfig(), and Reco().

Int_t PTSimModule::fRandomSeed [private]

Definition at line 56 of file PTSimModule.h.

Referenced by Config(), DefaultConfig(), GetRandomSeed(), and Reco().

TStopwatch PTSimModule::fStopwatch [private]

Definition at line 57 of file PTSimModule.h.

Referenced by EndJob(), and Reco().


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

Generated on 22 Nov 2017 for loon by  doxygen 1.6.1