loon: MCMerge.cxx Source File

00001 
00002 // $Id: MCMerge.cxx,v 1.15 2006/12/13 01:02:41 arms Exp $
00003 //
00004 // Flexibly merge two sets of Monte Carlo events at the hits/truth level.
00005 // Primary design is for rock/detector event merge (C++ replacement for
00006 // some of reco_minos's capability). Output should be sent to
00007 // PhotonTransport & DetSim next.
00008 //
00009 // K. Arms
00010 // arms@physics.umn.edu
00012 #include <MCMerge/MCMerge.h>
00013 #include <MCMerge/ConcatArrays.h>
00014 using namespace ConcatArrays;
00015 
00016 #include <stdlib.h>
00017 #include <iostream>
00018 #include <sstream>
00019 
00020 #include <MessageService/MsgService.h>
00021 #include <MinosObjectMap/MomNavigator.h>
00022 #include <JobControl/JobCModuleRegistry.h> // For JOBMODULE macro
00023 #include <Record/SimSnarlRecord.h>
00024 #include <Record/SimSnarlHeader.h>
00025 #include <Conventions/Detector.h>
00026 #include <Conventions/SimFlag.h>
00027 
00028 #ifdef EVENT_KINEMATICS_PKGS
00029 #include <EventKinematics/NuEvtKin.h>
00030 #include <EventKinematics/StdHepUtil.h>
00031 using namespace StdHepUtil;
00032 #endif
00033 
00034 #include <Digitization/DigiScintHit.h>
00035 #include <REROOT_Classes/REROOT_StdHep.h>
00036 #include <REROOT_Classes/REROOT_StdHepHead.h>
00037 #include <REROOT_Classes/REROOT_NeuKin.h>
00038 #include <REROOT_Classes/REROOT_FluxInfo.h>
00039 #include <REROOT_Classes/REROOT_FluxWgt.h>
00040 #include "TParticle.h"
00041 
00042 #include "TClonesArray.h"
00043 #include "TRandom3.h"
00044 #include "TSystem.h"
00045 
00046 JOBMODULE(MCMerge, "MCMerge",
00047           "Flexibly merge two sets of Monte Carlo at the hits/truth level");
00048 
00049 CVSID("$Id: MCMerge.cxx,v 1.15 2006/12/13 01:02:41 arms Exp $");
00050 
00051 typedef std::vector< const SimSnarlRecord* >::iterator mcmVecItr;
00052 
00053 //______________________________________________________________________
00054 
00055 MCMerge::MCMerge() :
00056   fBucketFreqMHz ( 53.1 ),
00057   fCurrTS        ( VldTimeStamp() ),
00058   fCurrVld       ( VldContext() ),
00059   fEvtVld        ( true ),
00060   fFirstTS       ( VldTimeStamp() ),
00061   fFirstVld      ( VldContext() ),
00062   fFirstVldSet   ( false ),
00063   fNumBatches    ( 5 ),
00064   fNumBuckets    ( 86 ),
00065   fOutStreamName ( "SimSnarl" ),
00066   fRandomSeed    ( 0 ),
00067   fRanGen        ( 0 ),
00068   fRunn          ( 1000 ),
00069   fSnarlNum      ( -1 ),
00070   fSubRunn       ( 0 ),
00071   fUsrVld        ( false )
00072 {
00073   //
00074   // Default constructor
00075   //
00076 
00077 }
00078 //______________________________________________________________________
00079 
00080 MCMerge::~MCMerge()
00081 {
00082   //
00083   // Default destructor
00084   //
00085 
00086 }
00087 
00088 //______________________________________________________________________
00089 
00090 void MCMerge::BeginJob()
00091 {
00092   //
00093   // Sort out the first VldTimeStamp
00094   //
00095 
00096   std::string strFirstDateTime = fFirstDateTime;
00097   if ( !strFirstDateTime.empty() && strFirstDateTime.length() == 15 ) {
00098     std::string d      =       strFirstDateTime.substr( 0, 1).c_str()  ;
00099     Int_t       year   = atoi( strFirstDateTime.substr( 1, 4).c_str() );
00100     Int_t       month  = atoi( strFirstDateTime.substr( 5, 2).c_str() );
00101     Int_t       day    = atoi( strFirstDateTime.substr( 7, 2).c_str() );
00102     Int_t       hour   = atoi( strFirstDateTime.substr( 9, 2).c_str() );
00103     Int_t       minute = atoi( strFirstDateTime.substr(11, 2).c_str() );
00104     Int_t       second = atoi( strFirstDateTime.substr(13, 2).c_str() );
00105 
00106     // Make sure entered yyyymmddhhmmss are valid values
00107     if (year   > 2038 || year   < 1996) year   = 2038;
00108     if (month  >   12 || month  <    1) month  =    1;
00109     if (day    >   31 || day    <    1) day    =    1;
00110     if (hour   >   23 || hour   <    0) hour   =    0;
00111     if (minute >   59 || minute <    0) minute =    0;
00112     if (second >   59 || second <    0) second =    0;
00113 
00114     // set the 1st TimeStamp according to user input
00115     fFirstTS = VldTimeStamp(year, month,  day,
00116                             hour, minute, second);
00117 
00118     fUsrVld  = true;  // we received a valid time/date from the user
00119     fEvtVld  = false; // user date/time has precedence
00120   } // end if : user entered a valid date/time for 1st output record
00121   
00122   fRanGen = new TRandom3( fRandomSeed ); // initialize the RNG
00123 
00124   MSG("MCMerge",Msg::kDebug)
00125     << "MCMerge random seed : " << fRanGen->GetSeed() << endl;
00126 
00127 }
00128 
00129 //______________________________________________________________________
00130 
00131 void MCMerge::EndJob()
00132 {
00133   //
00134   // Give a summary (if appropriate)
00135   //
00136 
00137   /*
00138   MSG("MCMerge", Msg::kInfo)
00139     << std::string(75,'*') << endl
00140     << "MCMerge Summary: " << endl
00141     << std::string(75,'*') << endl
00142     << endl;
00143   */
00144 }
00145 
00146 //______________________________________________________________________
00147 
00148 Float_t MCMerge::CalcOffsetTime()
00149 {
00150 
00151   // Assume fNumBucket buckets (default 86), fNumBatches (default 5), and
00152   // each bucket (1/53.1 MHz) wide, equally filled & equally
00153   // spaced, and each uniform.
00154 
00155   const Float_t bucketWidthns = (1.e3)/fBucketFreqMHz;
00156   const Float_t batchWidthns  = fNumBuckets*bucketWidthns;
00157 
00158   const Int_t   nBucket       = fRanGen->Integer(fNumBuckets);
00159   const Int_t   nBatch        = fRanGen->Integer(fNumBatches);
00160   const Float_t bucketOffsetT = fRanGen->Uniform(bucketWidthns);
00161 
00162   return ( Float_t(nBucket) * bucketWidthns +
00163            Float_t(nBatch)  * batchWidthns  +
00164            bucketOffsetT                    ); //ns
00165 
00166 }
00167 
00168 //______________________________________________________________________
00169 
00170 void MCMerge::ClearSimSnarlList()
00171 {
00172   //
00173   // Clear the fSimSnarlList without creating a massive memory leak!
00174   //
00175 
00176   if ( fSimSnarlList.empty() ) return;
00177 
00178   for (mcmVecItr it = fSimSnarlList.begin(); it != fSimSnarlList.end(); it++) {
00179     if ( !(*it) ) continue;  // already null
00180     delete (*it); (*it) = 0; // delete the object & null the ptr
00181   } // end loop over vector of SimSnarlRecord ptrs
00182 
00183   fSimSnarlList.clear(); // empty the vector
00184 
00185 }
00186 
00187 //______________________________________________________________________
00188 
00189 bool MCMerge::GrabStreamEvents(MomNavigator* mom)
00190 {
00191 
00192   std::vector<TObject*> fragList      = mom->GetFragmentList();
00193   const Int_t           kNumFragments = fragList.size();
00194 
00195   std::string momFragString;
00196   for (std::vector<TObject*>::const_iterator it = fragList.begin();
00197        it != fragList.end(); ++it) {
00198     // Only want SimSnarlRecords
00199     if ( std::string((*(*it)).ClassName()) != "SimSnarlRecord") continue;
00200 
00201     const char* streamName;
00202     SimSnarlRecord* temssr = dynamic_cast<SimSnarlRecord*>( (*it) );
00203     temssr->GetTempTags().Get( "stream", streamName );
00204     momFragString +=
00205       "\n" + std::string(streamName) + " (" +
00206       std::string(temssr->ClassName()) + ")" +
00207       " with timestamp: ";
00208     momFragString +=
00209       std::string( temssr->GetVldContext()->GetTimeStamp().AsString() );
00210     std::ostringstream runnss; 
00211     runnss << temssr->GetSimSnarlHeader()->GetRun();
00212     momFragString += " (run number " + runnss.str() + ")";
00213 
00214     if (it+1 != fragList.end())
00215       momFragString += ", ";
00216 
00217     fSimSnarlList.push_back( dynamic_cast<SimSnarlRecord*>( (*it)->Clone() ) );
00218   }
00219   MSG("MCMerge", Msg::kDebug) << "MomNavigator had "
00220                               << kNumFragments
00221                               << " fragments"
00222                               << ( (kNumFragments>0)? ": " : "" )
00223                               << momFragString
00224                               << endl;
00225 
00226   // Remove mom's input SimSnarlRecords so BadThings(TM) don't happen later
00227   if ( !fSimSnarlList.empty() ) {
00228     mom->Clear();
00229     MSG("MCMerge", Msg::kVerbose)
00230       << "MomNavigator Cleared ; MCMerge owns "
00231       << fSimSnarlList.size() << " SimSnarlRecords" << endl
00232       << "MomNavigator contents: " << (*mom) << endl;
00233   }
00234 
00235   return ( !fSimSnarlList.empty() );
00236 }
00237 
00238 //______________________________________________________________________
00239 
00240 JobCResult MCMerge::Get(MomNavigator* mom)
00241 {
00242   //
00243   // Take all events pushed to mom from each input stream and push the
00244   // merged SimSnarlRecord to MomNavigator
00245   //
00246   // Merged SimSnarlRecord will only contain:
00247   // DigiScintHits, FluxInfo, FluxWgt, (NuEvtKin), NuKin, and StdHep
00248   //
00249 
00250   if ( !GrabStreamEvents(mom) )
00251     return JobCResult::kFailed; // Yes: fail if there are no SimSnarlRecords!!
00252 
00253   if ( fFirstVldSet ) {
00254     // Increment the base VldTimeStamp & VldContext for this output snarl
00255     fCurrTS.Add(fIncSeconds); 
00256     fCurrVld = VldContext( fCurrVld.GetDetector(),
00257                            fCurrVld.GetSimFlag(),
00258                            fCurrTS               );
00259   } // end if : we incremented the current base timestamp
00260   else {
00261     // We need to set the 1st VldContext & VldTimeStamp
00262     VldContext firstEvtVld = 
00263       (fSimSnarlList[0])->GetSimSnarlHeader()->GetVldContext();
00264     
00265     // Get the correct detector from the input MC
00266     Detector::Detector_t detector = firstEvtVld.GetDetector();
00267 
00268     // Get the correct simFlag from the input MC
00269     SimFlag::SimFlag_t   simFlag  = firstEvtVld.GetSimFlag();
00270 
00271     // Did we get a timestamp from the user?
00272     if (fUsrVld)
00273       fFirstVld = VldContext(detector, simFlag, fFirstTS);
00274     // Did the user override grabbing the VldContext from the input events?
00275     // If so, use 'now'
00276     else if (!fEvtVld)
00277       fFirstVld = VldContext(detector, simFlag, VldTimeStamp());
00278     // Otherwise we use the 1st VldContext encountered from the input events
00279     else
00280       fFirstVld = firstEvtVld;
00281 
00282     fCurrVld     = fFirstVld;
00283     fCurrTS      = fFirstVld.GetTimeStamp();
00284     fFirstVldSet = true;
00285 
00286     MSG("MCMerge",Msg::kInfo) << "\nMCMerge: First VldContext set "
00287                               << fFirstVld << endl << endl;
00288   } // end else : we have now set the first output VldContext
00289 
00290   // Create a seed for our output SimSnarlRecord
00291   // (pull the following from the input events)
00292   Int_t   run       = fRunn;
00293   Int_t   subrun    = fSubRunn;
00294   Short_t runtype   = 0;
00295   UInt_t  errcode   = 0;
00296   Int_t   snarl     = ++fSnarlNum;
00297   UInt_t  trigsrc   = 0;
00298   Int_t   timeframe = fCurrTS.GetSec() - fFirstTS.GetSec();
00299   Int_t   spilltype = -1;
00300 
00301   VldTimeStamp mcGenTS;
00302   std::string  codename = SimSnarlHeader::TrimCodename("$Name: R1-28 $");
00303   std::string  hostname(gSystem->HostName());
00304   /*
00305   MAXMSG("MCMerge",Msg::kWarning,10)
00306     << "using 'now' " << mcGenTS.AsString("sql")
00307     << " for generation timestamp, local machine, & tag"
00308     << endl;
00309   */
00310   // the last three should come from the input data not just locally
00311   // generated ... though what to do if the overlay files are different?
00312 
00313   SimSnarlHeader simheader(fCurrVld,run,subrun,runtype,errcode,
00314                            snarl,trigsrc,timeframe,spilltype,
00315                            mcGenTS,codename,hostname);
00316 
00317   SimSnarlRecord* record = new SimSnarlRecord(simheader);
00318 
00319   MSG("MCMerge",Msg::kVerbose)
00320     << "New SimSnarlRecord\n" << simheader << endl;
00321 
00322   // Loop over our list of input snarls and
00323   // add them to the SimSnarlRecord components
00324   TClonesArray* digiScintHitsList = new TClonesArray("DigiScintHit",1);
00325   digiScintHitsList->SetName("DigiScintHits");
00326 #ifdef EVENT_KINEMATICS_PKGS
00327   TClonesArray* nuEvtKinList      = new TClonesArray("NuEvtKin",1);
00328   nuEvtKinList->SetName("NuEvtKinList");
00329 #endif
00330   TClonesArray* nuKinList         = new TClonesArray("REROOT_NeuKin",1);
00331   nuKinList->SetName("NeuKinList");
00332 
00333   TClonesArray* stdHepList        = new TClonesArray("TParticle",1);
00334   stdHepList->SetName("StdHep");
00335 
00336   TClonesArray* fluxInfoList      = new TClonesArray("REROOT_FluxInfo",1);
00337   fluxInfoList->SetName("FluxInfoList");
00338 
00339   TClonesArray* fluxWgtList       = new TClonesArray("REROOT_FluxWgt",1);
00340   fluxWgtList->SetName("FluxWgtList");
00341 
00342   for (mcmVecItr itr = fSimSnarlList.begin();
00343        itr != fSimSnarlList.end(); itr++) {
00344     const SimSnarlRecord* ssr = (*itr);
00345 
00346     TClonesArray* thisDigiScintHits =
00347       dynamic_cast<TClonesArray*>(const_cast<TObject*>(ssr->FindComponent("TClonesArray","DigiScintHits")));
00348     TClonesArray* thisStdHep =
00349       dynamic_cast<TClonesArray*>(const_cast<TObject*>(ssr->FindComponent("TClonesArray","StdHep")));
00350 #ifdef EVENT_KINEMATICS_PKGS
00351     TClonesArray* thisNuEvtKin =
00352       dynamic_cast<TClonesArray*>(const_cast<TObject*>(ssr->FindComponent("TClonesArray","NuEvtKinList")));
00353 #endif
00354     TClonesArray* thisNuKin =
00355       dynamic_cast<TClonesArray*>(const_cast<TObject*>(ssr->FindComponent("TClonesArray","NeuKinList")));
00356     TClonesArray* thisFluxInfo =
00357       dynamic_cast<TClonesArray*>(const_cast<TObject*>(ssr->FindComponent("TClonesArray","FluxInfoList")));
00358     TClonesArray* thisFluxWgt =
00359       dynamic_cast<TClonesArray*>(const_cast<TObject*>(ssr->FindComponent("TClonesArray","FluxWgtList")));
00360 
00361     // randomly place the events within the spill window
00362     const Float_t offsetns = CalcOffsetTime();
00363 
00364     OffsetTime(thisDigiScintHits, offsetns);
00365     OffsetTime(thisStdHep,        offsetns);
00366 #ifdef EVENT_KINEMATICS_PKGS
00367     OffsetTime(thisNuEvtKin,      offsetns);
00368 #endif
00369     OffsetTime(thisNuKin,         offsetns);
00370     OffsetTime(thisFluxInfo,      offsetns);
00371     OffsetTime(thisFluxWgt,       offsetns);
00372 
00373     const Int_t permStdHepSz = stdHepList->GetEntriesFast();
00374 
00375     ConcatDigiScintHits(thisDigiScintHits, digiScintHitsList,
00376                         permStdHepSz);
00377     ConcatStdHep       (thisStdHep,        stdHepList);
00378 #ifdef EVENT_KINEMATICS_PKGS
00379     ConcatNuEvtKin     (thisNuEvtKin,      nuEvtKinList, stdHepList);
00380 #endif
00381     ConcatNuKin        (thisNuKin,         nuKinList);
00382     ConcatFluxInfo     (thisFluxInfo,      fluxInfoList);
00383     ConcatFluxWgt      (thisFluxWgt,       fluxWgtList);
00384   } // end loop over our input SimSnarlRecord list
00385 
00386   // Empty our local list of SimSnarlRecords (don't need them anymore)
00387   ClearSimSnarlList(); 
00388 
00389   // Give our merged components to the record to own
00390   record->AdoptComponent(digiScintHitsList);
00391   record->AdoptComponent(nuKinList);
00392   record->AdoptComponent(stdHepList);
00393   record->AdoptComponent(fluxInfoList);
00394   record->AdoptComponent(fluxWgtList);
00395 #ifdef EVENT_KINEMATICS_PKGS
00396   record->AdoptComponent(nuEvtKinList);
00397 #endif
00398 
00399   // Give the SimSnarlRecord to MomNavigator to hold as a "fragment"
00400   record->GetTempTags().Set("stream",fOutStreamName);
00401   mom->AdoptFragment(record);
00402 
00403   // Reset the run# on future MsgService calls.
00404   (MsgService::Instance())->SetCurrentRunSnarl(fRunn, fSnarlNum);
00405 
00406   MSG("MCMerge", Msg::kVerbose)
00407     << "MomNavigator contents: " << (*mom) << endl;
00408 
00409   return JobCResult::kAOK; // Done
00410 }
00411 
00412 //___________________________________________________________________________
00413 
00414 const Registry& MCMerge::DefaultConfig() const
00415 {
00416   //
00417   // Supply the default configuration for the module
00418   //
00419   MSG("MCMerge",Msg::kVerbose) << "In MCMerge::DefaultConfig" << endl;
00420 
00421   static Registry r; // Default configuration for module
00422 
00423   // Set name of config
00424   std::string name = this->GetName();
00425   name += ".config.default";
00426   r.SetName(name.c_str());
00427 
00428   // Set values in configuration
00429   r.UnLockValues();
00430   r.Set("BucketFreqMHz",53.1); // bucket frequency
00431   r.Set("EvtVld",1);           // 1st vld is pulled from input events by default
00432   r.Set("FirstDateTime","");   // empty
00433   r.Set("DeltaT",1.9);         // default 1.9 seconds between spills
00434   r.Set("OutStreamName",
00435         "SimSnarl");           // generically the name of the class
00436   r.Set("RandomSeed",0);       // random seed
00437   r.Set("Runn",1000);          // run number
00438   r.Set("NumBatches",5);       // # batches of (NumBuckets) buckets
00439   r.Set("NumBuckets",86);      // # (1/.53.1MHz width) buckets in spill window
00440   r.Set("SubRunn",0);          // subrun number
00441   r.LockValues();
00442 
00443   return r;
00444 }
00445 
00446 //______________________________________________________________________
00447 
00448 void MCMerge::Config(const Registry& r)
00449 {
00450   //
00451   // Configure the module given the Registry r
00452   //
00453 
00454   MSG("MCMerge",Msg::kVerbose) << "In MCMerge::Config" << endl;
00455 
00456   fBucketFreqMHz = r.GetDouble("BucketFreqMHz");
00457   fEvtVld        = r.GetInt("EvtVld");
00458   fFirstDateTime = r.GetCharString("FirstDateTime");
00459   fIncSeconds    = r.GetDouble("DeltaT");
00460   fOutStreamName = r.GetCharString("OutStreamName");
00461   fRandomSeed    = r.GetInt("RandomSeed");
00462   fRunn          = r.GetInt("Runn");
00463   fNumBatches    = r.GetInt("NumBatches");
00464   fNumBuckets    = r.GetInt("NumBuckets");
00465   fSubRunn       = r.GetInt("SubRunn");
00466 
00467 }
00468 
00469 //______________________________________________________________________
00470 
00471 void MCMerge::Help()
00472 {
00473   //
00474   // Config() settings
00475   //
00476 
00477   cout << "Adjustable parameters and default settings:" << endl
00478        << "BucketFreqMHz (53.1)     -- Used to distribute events " << endl
00479        << "                            w/in the snarl window" << endl
00480        << "NumBatches (5)           -- Used to distribute events " << endl
00481        << "                            w/in the snarl window" << endl
00482        << "NumBuckets (86)          -- Used to distribute events " << endl
00483        << "                            w/in the snarl window." << endl
00484        << "The snarl window is NumBatches*NumBuckets/BucketFreqMHz uS wide"
00485        << endl
00486        << "Runn (1000)              -- Output run number" << endl
00487        << "SubRunn (0)              -- Output subrun number" << endl
00488        << endl
00489        << "DeltaT (1.9)             -- The time between timestamps of"
00490        << "                            output snarls (s)"
00491        << endl
00492        << "OutStreamName (SimSnarl) -- Name of the SimSnarlRecord stream "
00493        << endl
00494        << "                            given to MomNavigator for further "
00495        << endl
00496        << "                            processing" << endl
00497        << "EvtVld (1)               -- Determine 1st output VldContext "
00498        << endl
00499        << "                            from 1st input MC event" << endl
00500        << "FirstDateTime (\"\")       -- If filled, this will be the first "
00501        << endl
00502        << "                            VldContext timestamp" << endl
00503        << "If EvtVld == 0 and FirstDateTime == \"\" then the first output "
00504        << endl
00505        << "  VldContext will have the 1st date/time encountered in the input" << endl
00506        << endl;
00507 }
00508 
00509 //______________________________________________________________________
00510