Main Page | Modules | Namespace List | Class Hierarchy | Alphabetical List | Class List | Directories | File List | Namespace Members | Class Members | File Members | Related Pages
MCMerge

MCMerge.cxx

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