GenieModule Class Reference

A module for generating nu interaction in MINOS using GENIE. More...

#include <GenieModule.h>

Inheritance diagram for GenieModule:
JobCModule

List of all members.

Public Member Functions

 GenieModule ()
virtual ~GenieModule ()
JobCResult Get (MomNavigator *mom)
JobCResult Ana (const MomNavigator *mom)
void BeginRun ()
void EndRun ()
void EndJob ()
void Report ()
const RegistryDefaultConfig () const
void Config (const Registry &r)

Private Member Functions

VldContext GetNextVldContext (bool advance=true)
void FillStdHep (SimSnarlRecord *simrec, genie::EventRecord *gevrec)
void FillNeuKin (SimSnarlRecord *simrec, genie::EventRecord *gevrec, int nooscpdg)
void FillNuEvtKin (SimSnarlRecord *simrec, genie::EventRecord *gevrec, int nooscpdg)
void FillFluxInfo (SimSnarlRecord *simrec, genie::EventRecord *gevrec)
void CreateGeomVolSelector (std::string volsel)
void GetHallExtent (double *xyzmin, double *xyzmax, double &zlo, double &zhi)
void SetTimeStamp (std::string str, VldTimeStamp &ts)

Private Attributes

Detector::Detector_t fDetector
UShort_t fRun
UShort_t fSubRun
UShort_t fRunType
UInt_t fTrigSrc
Int_t fSpillType
Int_t fSnarl
Int_t fTimeFrame
Double_t fSpillInterval
VldTimeStamp fSpillTimeStamp
Int_t fTrigLimit
Int_t fTrigGen
Double_t fPOTsLimit
Double_t fPOTsGen
Int_t fMinSnarlDump
Int_t fMaxSnarlDump
Int_t fStdHepDumpMode
genie::GFluxI * fGFluxI
genie::geometry::ROOTGeomAnalyzer * fROOTGeomAnalyzer
genie::GMCJDriver * fGMCJDriver
TVector3 fMonoRayOrigin
double fMonoDx
double fMonoDy

Detailed Description

A module for generating nu interaction in MINOS using GENIE.

Author:
(last to touch it)
Author
rhatcher
Version:
Revision
1.7
Date:
Date
2010/09/03 21:19:46

Contact: rhatcher@fnal.gov

Id
GenieModule.h,v 1.7 2010/09/03 21:19:46 rhatcher Exp

Definition at line 49 of file GenieModule.h.


Constructor & Destructor Documentation

GenieModule::GenieModule (  ) 

Definition at line 68 of file GenieModule.cxx.

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

00069     : JobCModule(),
00070       fDetector(Detector::kUnknown),
00071       fRun(0), fSubRun(0), fRunType(0), fTrigSrc(0), fSpillType(0),
00072       fSnarl(0), fTimeFrame(0), fSpillInterval(1.9), fSpillTimeStamp(),
00073       fTrigLimit(99999999), fTrigGen(0), fPOTsLimit(1.0e30), fPOTsGen(0),
00074       fMinSnarlDump(-1), fMaxSnarlDump(-1), fStdHepDumpMode(-1),
00075       fGFluxI(0), // fGNuMIFlux(0), fMonoFlux(0), 
00076       fROOTGeomAnalyzer(0), fGMCJDriver(0),
00077       fMonoRayOrigin(1.48,0.10,0.0), fMonoDx(0.25), fMonoDy(0.25)
00078 {
00079     MSG("GenieModule",Msg::kDebug) << " GenieModule constructor" << endl;
00080 
00081     // initialize GENIE version of TDatabasePDG early
00082     //genie::PDGLibrary* pdglib =
00083     genie::PDGLibrary::Instance();
00084     //genie will load TDatabasePDG w/ GENIE table so we should not!
00085     LoadMinosPDG(false,true);   // fakeout
00086 
00087 }

GenieModule::~GenieModule (  )  [virtual]

Definition at line 91 of file GenieModule.cxx.

References fGMCJDriver, Msg::kDebug, and MSG.

00092 {
00093     MSG("GenieModule",Msg::kDebug) << " GenieModule destructor" << endl;
00094     // delete owned objects, assume deleting null pointers is okay
00095     delete fGMCJDriver;
00096     //owned by GMCJDriver//delete fROOTGeomAnalyzer;
00097     //owned by GMCJDriver//delete fGFluxI;
00098     //delete fMonoFlux;
00099     //delete fGNuMIFlux;
00100 }


Member Function Documentation

JobCResult GenieModule::Ana ( const MomNavigator mom  )  [virtual]

Implement this for read only access to the MomNavigator

Reimplemented from JobCModule.

Definition at line 1034 of file GenieModule.cxx.

References RecDataRecord< T >::FindComponent(), NuEvtKin::FormatToOStream(), MomNavigator::GetFragment(), and JobCResult::kPassed.

01035 {
01036 
01037   SimSnarlRecord* simrec =
01038     dynamic_cast<SimSnarlRecord*>(mom->GetFragment("SimSnarlRecord"));
01039 
01040   const TClonesArray* nuevtkinlist =
01041     dynamic_cast<const TClonesArray*>
01042     (simrec->FindComponent("TClonesArray","NuEvtKinList"));
01043 
01044   // test for finding one ...
01045   for (int ik=0; ik <= nuevtkinlist->GetLast(); ++ik) {
01046     const NuEvtKin* nuevtkin =
01047       dynamic_cast<const NuEvtKin*>(nuevtkinlist->At(ik));
01048     // test for getting object
01049     nuevtkin->FormatToOStream(cout,"v");
01050   }
01051   return JobCResult::kPassed;
01052 }

void GenieModule::BeginRun (  )  [virtual]

Implement for notification of begin of run (meaningful for Daq data only). See GetCurrentRun().

Reimplemented from JobCModule.

Definition at line 104 of file GenieModule.cxx.

References bfld::AsString(), CreateGeomVolSelector(), fDetector, fGFluxI, fGMCJDriver, fMonoDx, fMonoDy, fMonoRayOrigin, fROOTGeomAnalyzer, Registry::Get(), Registry::GetCharString(), JobCModule::GetConfig(), GetNextVldContext(), gGeoManager, gSystem(), UgliLoanPool::Instance(), UgliGeomHandle::IsGeo(), Msg::kDebug, Msg::kFatal, Munits::kilogram, Msg::kInfo, Msg::kWarning, Munits::meter, Munits::meter3, MSG, and UgliLoanPool::SetUseGeo().

00105 {
00106     MSG("GenieModule",Msg::kDebug) << " GenieModule BeginRun" << endl;
00107 
00108     Registry& r = GetConfig();
00109 
00110     fGMCJDriver = new genie::GMCJDriver();
00111 
00112     const char* fstr = "none";
00113     r.Get("Flux",fstr);
00114     std::string flxstring(fstr);
00115 
00116     if ( flxstring == "GNuMIFlux" ) {
00117       // construct GNuMIFlux object and config to use
00118       genie::flux::GNuMIFlux* fGNuMIFlux  = new genie::flux::GNuMIFlux();
00119       fGFluxI = fGNuMIFlux;  // ptr held by module is generic
00120       std::string cfg = "MINOS-";
00121       cfg += Detector::AsString(fDetector);
00122       cfg += r.GetCharString("FluxOpt");
00123       MSG("GenieModule",Msg::kInfo)
00124         << " GenieModule GNuMIFlux cfg: \"" << cfg << "\"" << endl;
00125 
00126       // configure flux w/ XML param set and flux file pattern
00127       const char* fluxpatt = r.GetCharString("FluxFiles");
00128       std::string fluxfname(gSystem->ExpandPathName(fluxpatt));
00129       fGNuMIFlux->LoadBeamSimData(fluxfname,cfg);
00130       double m_unit = genie::utils::units::UnitFromString("m");
00131       fGNuMIFlux->SetLengthUnits(m_unit);
00132       cout << " before PrintConfig " << endl;
00133       fGNuMIFlux->PrintConfig();
00134       cout << " after PrintConfig " << endl;
00135 
00136       // allow user to control pushing flux back to fixed z
00137       double upstreamz = -3.4e38;
00138       r.Get("FluxUpstreamZ",upstreamz);
00139       if ( TMath::Abs(upstreamz) < 1.0e30 ) fGNuMIFlux->SetUpstreamZ(upstreamz);
00140 
00141     }
00142 #ifndef MISSING_GSIMPLENTPFLUX
00143     else if ( flxstring == "GNuMIFlux" ) {
00144       // construct GNuMIFlux object and config to use
00145       genie::flux::GSimpleNtpFlux* fGSimpleNtpFlux  = new genie::flux::GSimpleNtpFlux();
00146       fGFluxI = fGSimpleNtpFlux;  // ptr held by module is generic
00147       std::string cfg = "MINOS-";
00148       cfg += Detector::AsString(fDetector);
00149       cfg += "Det";
00150       cfg += r.GetCharString("FluxOpt");
00151       MSG("GenieModule",Msg::kInfo)
00152         << " GenieModule GSimpleNtpFlux cfg: \"" << cfg << "\"" << endl;
00153 
00154       // configure flux w/ XML param set and flux file pattern
00155       const char* fluxpatt = r.GetCharString("FluxFiles");
00156       std::string fluxfname(gSystem->ExpandPathName(fluxpatt));
00157       fGSimpleNtpFlux->LoadBeamSimData(fluxfname,cfg);
00158 
00159       fGSimpleNtpFlux->PrintConfig();
00160 
00161     }
00162 #endif
00163     else if ( flxstring == "Mono" ) {
00164       genie::flux::GMonoEnergeticFlux* fMonoFlux = new genie::flux::GMonoEnergeticFlux(10.,14);
00165       fGFluxI = fMonoFlux;  // ptr held by module is generic
00166       fMonoFlux->SetDirectionCos(0.,0.,1.);
00167       TVector3& xyz0 = fMonoRayOrigin;
00168       fMonoFlux->SetRayOrigin(xyz0.X(),xyz0.Y(),xyz0.Z()); // flux units: meters
00169       r.Get("MonoDx",fMonoDx);
00170       r.Get("MonoDy",fMonoDy);
00171 
00172     }
00173     else {
00174       MSG("GenieModule",Msg::kFatal)
00175         << " GenieModule no flux handling for \"" << flxstring << "\"" << endl;
00176     }
00177 
00178     VldContext vldc = GetNextVldContext(false);
00179     // Force GeoGeometry (TGeo) style geometry
00180     UgliLoanPool::Instance()->SetUseGeo(true);
00181     UgliGeomHandle ugh(vldc);
00182 
00183     if ( ! ugh.IsGeo() || ! gGeoManager ) {
00184       MSG("GenieModule",Msg::kFatal) << " GenieModule not GeoGeometry" << endl;
00185     }
00186 
00187 #ifdef RWH_TEST
00188     int useOldROOTGeomAnalyzer = 0;
00189     r.Get("OldROOTGeomAnalyzer",useOldROOTGeomAnalyzer);
00190     if ( useOldROOTGeomAnalyzer ) {
00191       MSG("GenieModule",Msg::kWarning)
00192         << "Using ROOTGeomAnalyzerOld!" << endl;
00193       fROOTGeomAnalyzer = new genie::geometry::ROOTGeomAnalyzerOld(gGeoManager);
00194     } else 
00195 #endif
00196     fROOTGeomAnalyzer = new genie::geometry::ROOTGeomAnalyzer(gGeoManager);
00197     // ! important otherwise mixtures are weighted by 1/Nelements
00198     fROOTGeomAnalyzer->SetMixtureWeightsSum(1.);
00199     fROOTGeomAnalyzer->SetWeightWithDensity(true);
00200     std::string lunits_name = "cm";
00201     double lunits = genie::utils::units::UnitFromString(lunits_name);
00202     MSG("GenieModule",Msg::kInfo) << "set length units to " << lunits_name
00203                                   << " " << lunits << " ratio "
00204                                   << lunits/genie::units::meter
00205                                   << std::endl;
00206     fROOTGeomAnalyzer->SetLengthUnits(lunits);
00207     std::string dunits_name = "g_cm3";
00208     double dunits = genie::utils::units::UnitFromString(dunits_name);
00209     MSG("GenieModule",Msg::kInfo) << "set density units to " << dunits_name
00210                                   << " " << dunits << " ratio "
00211                                   << dunits/(genie::units::kilogram/genie::units::meter3)
00212                                   << std::endl;
00213     fROOTGeomAnalyzer->SetDensityUnits(dunits);
00214 
00215     std::string volsel = r.GetCharString("VolumeSelector");
00216     std::transform(volsel.begin(),volsel.end(),volsel.begin(),::toupper); // UPCASE
00217     std::string topvol = r.GetCharString("TopVol");
00218     if ( topvol != "" ) {
00219       if ( ( volsel == "BOTH" || volsel == "ROCK" ) && topvol == "HALL" ) {
00220         MSG("GenieModule",Msg::kWarning) 
00221           << "Use of GeomVolSelector '" << volsel << "' and TopVol '" << topvol << "' "
00222           << "might be in conflict" << std::endl;
00223       }
00224       MSG("GenieModule",Msg::kInfo) << "SetTopVolName: '" << topvol << "'" << endl;
00225       fROOTGeomAnalyzer->SetTopVolName(topvol);
00226     }
00227     // must come *after* any topvol is defined
00228     CreateGeomVolSelector(volsel);
00229 
00230     // adjust how careful the scan for max path length is
00231     int nscan = 0;
00232     r.Get("NScan",nscan);
00233     if ( nscan != 0 ) {
00234       // use flux method to scan geometry for max path lengths
00235       MSG("GenieModule",Msg::kInfo) << "Use flux to scan geometry for max pathlength " 
00236                                     << nscan << " particles (-1=default)" << std::endl;
00237       if ( nscan > 0 ) fROOTGeomAnalyzer->SetScannerNParticles(nscan);
00238       fROOTGeomAnalyzer->SetScannerFlux(fGFluxI);
00239     } else {
00240       int npts = 0, nrays = 0;
00241       r.Get("NPoints",npts);
00242       r.Get("NRays",nrays);
00243       MSG("GenieModule",Msg::kInfo) << "Use box to scan geometry for max pathlength " 
00244                                     << std::endl;
00245       if ( npts > 0 ) {
00246         MSG("GenieModule",Msg::kInfo) << "set npoints " << npts << std::endl;
00247         fROOTGeomAnalyzer->SetScannerNPoints(npts);
00248       }
00249       if ( nrays > 0 ) {
00250         MSG("GenieModule",Msg::kInfo) << "set nrays " << nrays << std::endl;
00251         fROOTGeomAnalyzer->SetScannerNRays(nrays);
00252       }
00253     }
00254 
00255     // path length fudge factor adjustment
00256     double safety = 0;
00257     if ( r.Get("MaxPlSafety",safety) && safety > 0 )
00258       fROOTGeomAnalyzer->SetMaxPlSafetyFactor(safety);
00259            
00260     fGMCJDriver->UseFluxDriver(fGFluxI);
00261     fGMCJDriver->UseGeomAnalyzer(fROOTGeomAnalyzer);
00262     std::string maxplxml = r.GetCharString("MaxPathLengthXML");
00263     if ( maxplxml != "" ) {
00264       MSG("GenieModule",Msg::kInfo) << "UseMaxPathLengths: " << maxplxml << endl;
00265       fGMCJDriver->UseMaxPathLengths(maxplxml);
00266     }
00267     fGMCJDriver->Configure();
00268     fGMCJDriver->UseSplines();
00269     fGMCJDriver->ForceSingleProbScale();
00270 
00271 }

void GenieModule::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 387 of file GenieModule.cxx.

References bfld::AsString(), VldTimeStamp::AsString(), fDetector, fMaxSnarlDump, fMinSnarlDump, fPOTsLimit, fRun, fRunType, fSnarl, fSpillInterval, fSpillTimeStamp, fSpillType, fStdHepDumpMode, fSubRun, fTimeFrame, fTrigLimit, fTrigSrc, Registry::Get(), Registry::GetCharString(), Msg::kInfo, MSG, SetTimeStamp(), and Detector::StringToEnum().

00388 {
00389 //======================================================================
00390 // Configure the module using the registry r
00391 //======================================================================
00392   int    tmpi;
00393   double tmpd;
00394   const char* tmpcs = 0;
00395 
00396   if (r.Get("Detector",       tmpcs)) 
00397     fDetector = Detector::StringToEnum(tmpcs);
00398   if (r.Get("Run",            tmpi))  fRun            = tmpi;
00399   if (r.Get("SubRun",         tmpi))  fSubRun         = tmpi;
00400   if (r.Get("RunType",        tmpi))  fRunType        = tmpi;
00401   if (r.Get("TrigSrc",        tmpi))  fTrigSrc        = tmpi;
00402   if (r.Get("SpillType",      tmpi))  fSpillType      = tmpi;
00403   if (r.Get("Snarl",          tmpi))  fSnarl          = tmpi;
00404   if (r.Get("TimeFrame",      tmpi))  fTimeFrame      = tmpi;
00405   if (r.Get("StartTimeStamp", tmpcs)) SetTimeStamp(tmpcs,fSpillTimeStamp);
00406   if (r.Get("SpillInterval",  tmpd))  fSpillInterval  = tmpd;
00407 
00408 
00409   if (r.Get("TrigLimit",      tmpi))  fTrigLimit      = tmpi;
00410   if (r.Get("POTsLimit",      tmpd))  fPOTsLimit      = tmpd;
00411 
00412   if (r.Get("MinSnarlDump",   tmpi))  fMinSnarlDump   = tmpi;
00413   if (r.Get("MaxSnarlDump",   tmpi))  fMaxSnarlDump   = tmpi;
00414   if (r.Get("StdHepDumpMode", tmpi))  fStdHepDumpMode = tmpi;
00415 
00416   MSG("GenieModule",Msg::kInfo)
00417       << "GenieModule::Config() " << endl
00418       << "  Detector " << Detector::AsString(fDetector)
00419       << " Run " << fRun << " SubRun " << fSubRun
00420       << " RunType " << fRunType << " TrigSrc " << fTrigSrc
00421       << " SpillType " << fSpillType << endl
00422       << "  Next:  Snarl " << fSnarl << " TimeFrame " << fTimeFrame
00423       << " @ '" << fSpillTimeStamp.AsString("c") << "'"
00424       << " delta-t " << fSpillInterval
00425       << endl;
00426 
00427   MSG("GenieModule",Msg::kInfo)
00428       << "  Limits: " << fTrigLimit << " triggers, " << fPOTsLimit << " POTs"
00429       << endl;
00430 
00431   MSG("GenieModule",Msg::kInfo)
00432     << "  Flux:        " << r.GetCharString("Flux") << endl;
00433 
00434 }

void GenieModule::CreateGeomVolSelector ( std::string  volsel  )  [private]

Definition at line 1056 of file GenieModule.cxx.

References MuELoss::e, fROOTGeomAnalyzer, Registry::Get(), JobCModule::GetConfig(), GetHallExtent(), Msg::kFatal, Msg::kInfo, and MSG.

Referenced by BeginRun().

01057 {
01058 
01059   MSG("GenieModule",Msg::kInfo)
01060     << "CreateGeomVolSelector: '" << volsel << "'" << std::endl;
01061 
01062   Registry& r = GetConfig();
01063 
01064   double hallxyzmin[3], hallxyzmax[3], zlo, zhi;
01065   GetHallExtent(hallxyzmin,hallxyzmax,zlo,zhi); 
01066 
01067   genie::geometry::GeomVolSelectorFiducial* vselfid = 0;
01068   if ( volsel == "DETECTOR" ) {
01069     // simple volume selector
01070     MSG("GenieModule",Msg::kInfo) << "CreateGeomVolSelector: Fiducial" << endl;
01071     vselfid = new genie::geometry::GeomVolSelectorFiducial();
01072   }
01073   if ( volsel == "ROCK" || volsel == "BOTH" ) {
01074     // rock selector
01075 #ifdef MISSING_GEOMVOLSELECTORROCKBOX
01076     vselfid = 0;
01077     MSG("GenieModule",Msg::kFatal) 
01078       << " GenieModule no GeomVolSelectorRockBox class" << endl;
01079     assert(0);
01080 #else
01081     MSG("GenieModule",Msg::kInfo) << "CreateGeomVolSelector: RockBox" << endl;
01082     genie::geometry::GeomVolSelectorRockBox* vselrock =
01083       new genie::geometry::GeomVolSelectorRockBox();
01084     vselfid = vselrock; // base class
01085 
01086     // configure "rock" portion of selector
01087     // extents *should* come from DB or user or whatever ...
01088     /*
01089     double nearhmin[3] = { 0, 0, 0 };
01090     double nearhmax[3] = { 0, 0, 0 };
01091     double farhmin[3]  = { 0, 0, 0 };
01092     double farhmax[3]  = { 0, 0, 0 };
01093     double* hallxyzmin = ( fDetector == Detector::kNear) ? nearhmin : farhmin;
01094     double* hallxyzmax = ( fDetector == Detector::kNear) ? nearhmax : farhmax;
01095     */
01096     vselrock->SetRockBoxMinimal(hallxyzmin,hallxyzmax); // base size of cavern
01097     double wallmin = 800.; // assume geometry in cm ( 8 meters) 
01098     double dedx = 2.5*1.7e-3; // assume GeV/cm, rho = 2.5 "rock-like" = 1.7e-3
01099     r.Get("MinimumWall",wallmin);
01100     r.Get("DeDx",dedx);
01101     vselrock->SetMinimumWall(wallmin);
01102     vselrock->SetDeDx(dedx);
01103 
01104     if ( volsel == "BOTH" ) {
01105       // need to set some minimalistic "fiducial volume" that 
01106       // it can "exclude"; just make it a tiny tiny bubble
01107       vselfid->MakeSphere(0,0,0,1.0e-10);
01108     }
01109 #endif
01110   }
01111   if ( vselfid && volsel != "BOTH" ) {
01112     // configure "detector" portion of selector
01113     // extents *should* come from DB or user or whatever ...
01114     /*
01115     double neardmin[3] = { 0, 0, 0 };
01116     double neardmax[3] = { 0, 0, 0 };
01117     double fardmin[3]  = { 0, 0, 0 };
01118     double fardmax[3]  = { 0, 0, 0 };
01119     double* detxyzmin = ( fDetector == Detector::kNear) ? neardmin : fardmin;
01120     double* detxyzmax = ( fDetector == Detector::kNear) ? neardmax : fardmax;
01121     */
01122     for (unsigned int i=0; i<2; ++i) {
01123       const double epsilon = 0.01;
01124       hallxyzmin[i] += epsilon;
01125       hallxyzmax[i] -= epsilon;
01126     }
01127     hallxyzmin[2] = zlo;
01128     hallxyzmax[2] = zhi;
01129     vselfid->MakeBox(hallxyzmin,hallxyzmax);
01130   }
01131   if ( vselfid ) {
01132     MSG("GenieModule",Msg::kInfo) << "CreateGeomVolSelector: ShapeMaster2Top" << endl;
01133     vselfid->ConvertShapeMaster2Top(fROOTGeomAnalyzer);
01134     fROOTGeomAnalyzer->AdoptGeomVolSelector(vselfid);
01135   }
01136 
01137 }

const Registry & GenieModule::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 299 of file GenieModule.cxx.

References MuELoss::e, JobCModule::GetName(), gSystem(), Registry::LockValues(), Registry::Set(), and Registry::UnLockValues().

00300 {
00301 //======================================================================
00302 // The default configuration for this module
00303 //======================================================================
00304   //int itrue  = 1; // Work around for Registry's lack of bool
00305   //int ifalse = 0; // Work around for Registry's lack of bool
00306 
00307   static Registry r;
00308 
00309   std::string name = this->GetName();
00310   name += ".config.default";
00311   r.SetName(name.c_str());
00312 
00313   r.UnLockValues();
00314 
00315   // General Run Parameters
00316   r.Set("Run",               9999);
00317   r.Set("SubRun",               0);
00318   r.Set("RunType",              0); // need to set this
00319   r.Set("TrigSrc",              0); // need to set this
00320   r.Set("SpillType",            0); // need to set this
00321 
00322   r.Set("Snarl",                0); // next snarl #
00323   r.Set("TimeFrame",            0); // next timeframe #
00324 
00325   r.Set("TrigLimit",   2147483647);
00326   r.Set("POTsLimit",       2.1e30);
00327 
00328   const char* dfltStartTimeStamp = "";  // "now" by default
00329   const char* charStartTimeStamp = gSystem->Getenv("STARTTIMESTAMP");
00330   if ( charStartTimeStamp == 0 ) charStartTimeStamp = dfltStartTimeStamp;
00331   r.Set("StartTimeStamp",     charStartTimeStamp);
00332   r.Set("SpillTimeInterval",  1.9); // time between spills
00333 
00334   r.Set("MinSnarlDump",         -1);
00335   r.Set("MaxSnarlDump",         -1);
00336   r.Set("StdHepDumpMode",        2); // -1);
00337 
00338   //r.Set("InterConfig","ConfigName='MODBYRS' ConfigVersion=2 "
00339   //                    " RecalcXsec=0 "
00340   //                    " Models='DISMODEL=2 DISRES_OVERLAP=1' "
00341   //                    " Params='QEL_MA=1.032 RES_MA=1.032'");
00342 
00343   // Flux Parameters
00344   r.Set("Flux",       "GNuMIFlux");
00345   r.Set("Detector",     "Unknown");
00346   // GNuMIFlux config is MINOS-$(DET)$(FLUXOPT}
00347   // ie. MINOS-NearDet, MINOS-FarDet, MINOS-NearRock, MINOS-FarRock
00348   r.Set("FluxOpt",          "Det");   // modifier for GNuMIFlux XML param set
00349   r.Set("FluxFiles",
00350         "$FLUXPATH/v19/fluka05_le010z185i_run1/root/fluka05_le010z185i_*.root");
00351   //"$FLUXPATH/flugg/flugg_le010z185i/root/flugg_le010z185i_*.root"
00352   //"$FLUXPATH/v19/fluka05_le010z185i_lowcut/root/fluka05_le010z185i_*.root";
00353 
00354   r.Set("FluxUpstreamZ",  -3.4e38); // force nu back to fixed z off window
00355                                     // if |z0| > 1.0e30 leave on flux window
00356 
00357   // GeomAnalyzer Parameters
00358   r.Set("TopVol",          "HALL"); //
00359   r.Set("VolumeSelector",  "NONE"); // "DETECTOR" "ROCK" "BOTH" "NONE"
00360 
00361   r.Set("MaxPathLengthXML",       "");   // override geometry scan for max path lengths
00362 
00363   r.Set("NScan",                  -1);   // for initial material scan use flux (set 0 to use box method) -1 = use default [10000]
00364   r.Set("NPoints",                 0);   // for initial material scan - 0 = use default [200]
00365   r.Set("NRays",                   0);   // for initial material scan - 0 = use default [200]
00366   r.Set("MaxPlSafety",             0);   // use default 1.1 in ROOTGeomAnalyzer (path length fudge)
00367 
00368   r.Set("MaxPlSafety",            -1);   // use default 1.1 in ROOTGeomAnalyzer
00369 
00370   r.Set("MinimumWall",          800.);   // assume geometry in cm ( 8 meters) 
00371   r.Set("DeDx",           2.5*1.7e-3);   // assume GeV/cm, rho = 2.5 "rock-like" = 1.7e-3
00372 
00373   r.Set("MonoDx",               0.25);   // +/- 25cm jiggle in mono x pos
00374   r.Set("MonoDy",               0.25);   // +/- 25cm jiggle in mono y pos
00375 
00376 #ifdef RWH_TEST
00377   r.Set("OldROOTGeomAnalyzer",ifalse);
00378 #endif
00379 
00380   r.LockValues();
00381 
00382   return r;
00383 }

void GenieModule::EndJob (  )  [virtual]

Implement for notification of end of job

Reimplemented from JobCModule.

Definition at line 282 of file GenieModule.cxx.

References fPOTsGen, fTrigGen, Msg::kInfo, and MSG.

00283 {
00284     MSG("GenieModule",Msg::kInfo)
00285       << "GenieModule::EndJob generated "
00286       << fTrigGen << " triggers using " << fPOTsGen << " POTs"
00287       << endl;
00288 }

void GenieModule::EndRun (  )  [virtual]

Implement for notification of end of run (meaningful for Daq data only). See GetLastRun().

Reimplemented from JobCModule.

Definition at line 275 of file GenieModule.cxx.

References Msg::kDebug, and MSG.

00276 {
00277     MSG("GenieModule",Msg::kDebug) << " GenieModule EndRun" << endl;
00278 }

void GenieModule::FillFluxInfo ( SimSnarlRecord simrec,
genie::EventRecord *  gevrec 
) [private]

Definition at line 880 of file GenieModule.cxx.

References RecDataRecord< T >::AdoptComponent(), FLUXWGT_DEF::Beam, FLUXINFO_DEF::beampx, FLUXINFO_DEF::beampy, FLUXINFO_DEF::beampz, FLUXINFO_DEF::beamx, FLUXINFO_DEF::beamy, FLUXINFO_DEF::beamz, fGFluxI, FLUXINFO_DEF::FluxEvtNo, FLUXINFO_DEF::FluxRun, FLUXWGT_DEF::ID, FLUXINFO_DEF::ID, FLUXINFO_DEF::mupare, FLUXINFO_DEF::muparpx, FLUXINFO_DEF::muparpy, FLUXINFO_DEF::muparpz, FLUXINFO_DEF::Ndecay, FLUXINFO_DEF::Ndxdz, FLUXINFO_DEF::NdxdzFar, FLUXINFO_DEF::NdxdzNear, FLUXINFO_DEF::Ndydz, FLUXINFO_DEF::NdydzFar, FLUXINFO_DEF::NdydzNear, FLUXINFO_DEF::Necm, FLUXINFO_DEF::Nenergy, FLUXINFO_DEF::NenergyF, FLUXINFO_DEF::NenergyN, FLUXINFO_DEF::Nimpwt, FLUXINFO_DEF::Norig, FLUXINFO_DEF::Npz, FLUXINFO_DEF::Ntype, FLUXINFO_DEF::NWtFar, FLUXINFO_DEF::NWtNear, FLUXINFO_DEF::pdPx, FLUXINFO_DEF::pdPy, FLUXINFO_DEF::pdPz, FLUXINFO_DEF::ppdxdz, FLUXINFO_DEF::ppdydz, FLUXINFO_DEF::ppenergy, FLUXINFO_DEF::ppmedium, FLUXINFO_DEF::pppz, FLUXINFO_DEF::ppvx, FLUXINFO_DEF::ppvy, FLUXINFO_DEF::ppvz, FLUXINFO_DEF::ptype, FLUXINFO_DEF::tgen, FLUXINFO_DEF::tgppx, FLUXINFO_DEF::tgppy, FLUXINFO_DEF::tgppz, FLUXINFO_DEF::tgptype, FLUXINFO_DEF::tprivx, FLUXINFO_DEF::tprivy, FLUXINFO_DEF::tprivz, FLUXINFO_DEF::tptype, FLUXINFO_DEF::tpx, FLUXINFO_DEF::tpy, FLUXINFO_DEF::tpz, FLUXINFO_DEF::tvx, FLUXINFO_DEF::tvy, FLUXINFO_DEF::tvz, FLUXWGT_DEF::Version, FLUXINFO_DEF::Vx, FLUXINFO_DEF::Vy, FLUXINFO_DEF::Vz, FLUXWGT_DEF::Weight, FLUXWGT_DEF::WeightErr, FLUXINFO_DEF::xpoint, FLUXINFO_DEF::ypoint, and FLUXINFO_DEF::zpoint.

Referenced by Get().

00882 {
00883 
00884    Int_t entry = 0;
00885 
00886 //======================================================================
00887 // Fill the REROOT_FluxInfo part of simrec
00888 //======================================================================
00889 
00890    // New TClonesArray of the right size (1)
00891    TClonesArray* fluxinfos = new TClonesArray("REROOT_FluxInfo",1);
00892    fluxinfos->SetName("FluxInfoList");
00893    simrec->AdoptComponent(fluxinfos);
00894 
00895    static const genie::flux::GNuMIFluxPassThroughInfo gnumi_null;
00896 
00897    genie::flux::GNuMIFluxPassThroughInfo gnumi_pti(gnumi_null);
00898 
00899    genie::flux::GMonoEnergeticFlux* fMonoFlux = 
00900      dynamic_cast<genie::flux::GMonoEnergeticFlux*>(fGFluxI);
00901 
00902    genie::flux::GNuMIFlux* fGNuMIFlux = 
00903      dynamic_cast<genie::flux::GNuMIFlux*>(fGFluxI);
00904    if ( fGNuMIFlux ) gnumi_pti = fGNuMIFlux->PassThroughInfo();
00905 
00906 #ifndef MISSING_GSIMPLENTPFLUX
00907    genie::flux::GSimpleNtpFlux* fGSimpleNtpFlux = 
00908      dynamic_cast<genie::flux::GSimpleNtpFlux*>(fGFluxI);
00909    if ( fGSimpleNtpFlux ) {
00910      const genie::flux::GSimpleNtpEntry* simple_entry = 
00911        fGSimpleNtpFlux->GetCurrentEntry();
00912 
00913      gnumi_pti.ntype  = simple_entry->pdg;
00914      gnumi_pti.nimpwt = simple_entry->wgt;
00915 
00916      const genie::flux::GSimpleNtpNuMI*  simple_numi = 
00917        fGSimpleNtpFlux->GetCurrentNuMI();
00918      if ( simple_numi ) {
00919        gnumi_pti.run     = simple_numi->run;
00920        gnumi_pti.evtno   = simple_numi->evtno;
00921        gnumi_pti.tpx     = simple_numi->tpx;
00922        gnumi_pti.tpy     = simple_numi->tpy;
00923        gnumi_pti.tpz     = simple_numi->tpz;
00924        gnumi_pti.tptype  = simple_numi->tptype;
00925      }
00926    }
00927 #endif
00928 
00929    FLUXINFO_DEF* flxi = new FLUXINFO_DEF;
00930    flxi->ID        = 0;
00931    flxi->FluxRun   = gnumi_pti.run;
00932    flxi->FluxEvtNo = gnumi_pti.evtno;
00933    flxi->Ndxdz     = gnumi_pti.ndxdz;
00934    flxi->Ndydz     = gnumi_pti.ndydz;
00935    flxi->Npz       = gnumi_pti.npz;
00936    flxi->Nenergy   = gnumi_pti.nenergy;
00937    flxi->NdxdzNear = gnumi_pti.ndxdznea;
00938    flxi->NdydzNear = gnumi_pti.ndydznea;
00939    flxi->NenergyN  = gnumi_pti.nenergyn;
00940    flxi->NWtNear   = gnumi_pti.nwtnear;
00941    flxi->NdxdzFar  = gnumi_pti.ndxdzfar;
00942    flxi->NdydzFar  = gnumi_pti.ndydzfar;
00943    flxi->NenergyF  = gnumi_pti.nenergyf;
00944    flxi->NWtFar    = gnumi_pti.nwtfar;
00945    flxi->Norig     = gnumi_pti.norig;
00946    flxi->Ndecay    = gnumi_pti.ndecay;
00947    flxi->Ntype     = gnumi_pti.ntype;
00948    flxi->Vx        = gnumi_pti.vx;
00949    flxi->Vy        = gnumi_pti.vy;
00950    flxi->Vz        = gnumi_pti.vz;
00951    flxi->pdPx      = gnumi_pti.pdpx;
00952    flxi->pdPy      = gnumi_pti.pdpy;
00953    flxi->pdPz      = gnumi_pti.pdpz;
00954    flxi->ppdxdz    = gnumi_pti.ppdxdz;
00955    flxi->ppdydz    = gnumi_pti.ppdydz;
00956    flxi->pppz      = gnumi_pti.pppz;
00957    flxi->ppenergy  = gnumi_pti.ppenergy;
00958    flxi->ppmedium  = gnumi_pti.ppmedium;
00959    flxi->ptype     = gnumi_pti.ptype;
00960    flxi->ppvx      = gnumi_pti.ppvx;
00961    flxi->ppvy      = gnumi_pti.ppvy;
00962    flxi->ppvz      = gnumi_pti.ppvz;
00963    flxi->muparpx   = gnumi_pti.muparpx;
00964    flxi->muparpy   = gnumi_pti.muparpy;
00965    flxi->muparpz   = gnumi_pti.muparpz;
00966    flxi->mupare    = gnumi_pti.mupare;
00967    flxi->Necm      = gnumi_pti.necm;
00968    flxi->Nimpwt    = gnumi_pti.nimpwt;
00969    flxi->xpoint    = gnumi_pti.xpoint;
00970    flxi->ypoint    = gnumi_pti.ypoint;
00971    flxi->zpoint    = gnumi_pti.zpoint;
00972    flxi->tvx       = gnumi_pti.tvx;
00973    flxi->tvy       = gnumi_pti.tvy;
00974    flxi->tvz       = gnumi_pti.tvz;
00975    flxi->tpx       = gnumi_pti.tpx;
00976    flxi->tpy       = gnumi_pti.tpy;
00977    flxi->tpz       = gnumi_pti.tpz;
00978    flxi->tptype    = gnumi_pti.tptype;
00979    flxi->tgen      = gnumi_pti.tgen;
00980    flxi->tgptype   = gnumi_pti.tgptype;
00981    flxi->tgppx     = gnumi_pti.tgppx;
00982    flxi->tgppy     = gnumi_pti.tgppy;
00983    flxi->tgppz     = gnumi_pti.tgppz;
00984    flxi->tprivx    = gnumi_pti.tprivx;
00985    flxi->tprivy    = gnumi_pti.tprivy;
00986    flxi->tprivz    = gnumi_pti.tprivz;
00987    flxi->beamx     = gnumi_pti.beamx;
00988    flxi->beamy     = gnumi_pti.beamy;
00989    flxi->beamz     = gnumi_pti.beamz;
00990    flxi->beampx    = gnumi_pti.beampx;
00991    flxi->beampy    = gnumi_pti.beampy;
00992    flxi->beampz    = gnumi_pti.beampz;
00993 
00994    // new with placement
00995    new((*fluxinfos)[entry])REROOT_FluxInfo(flxi);  // bogus entry
00996    // fix up fluxinfo elements
00997 
00998    delete flxi;
00999 
01000 //======================================================================
01001 // Fill the REROOT_FluxWgt part of simrec
01002 //======================================================================
01003 
01004    // New TClonesArray of the right size (1)
01005    TClonesArray* fluxwgts = new TClonesArray("REROOT_FluxWgt",1);
01006    fluxwgts->SetName("FluxWgtList");
01007    simrec->AdoptComponent(fluxwgts);
01008 
01009    FLUXWGT_DEF* flxw = new FLUXWGT_DEF;
01010    flxw->ID        = 0;
01011    // beam name buffer is 32 chars, trailing \0 added at REROOT_FluxWgt stage
01012    const char * bname = "Unknown";
01013 
01014    if ( fMonoFlux       ) bname = "MonoFlux";
01015    if ( fGNuMIFlux      ) bname = "GNuMIFlux";
01016 #ifndef MISSING_GSIMPLENTPFLUX
01017    if ( fGSimpleNtpFlux ) bname = "GSimpleNtpFlux";
01018 #endif
01019 
01020    strncpy(flxw->Beam,bname,32);
01021    flxw->Version   = 0;
01022    flxw->Weight    = gevrec->Weight();
01023    flxw->WeightErr = 0;
01024 
01025    // new with placement
01026    new((*fluxwgts)[entry])REROOT_FluxWgt(flxw);  // bogus entry
01027    // fix up fluxinfo elements
01028 
01029    delete flxw;
01030 }

void GenieModule::FillNeuKin ( SimSnarlRecord simrec,
genie::EventRecord *  gevrec,
int  nooscpdg 
) [private]

Definition at line 626 of file GenieModule.cxx.

References NEUKIN_DEF::A, RecDataRecord< T >::AdoptComponent(), NEUKIN_DEF::EMFrac, RecDataRecord< T >::FindComponent(), NEUKIN_DEF::IAction, NEUKIN_DEF::IBoson, NEUKIN_DEF::ID, NEUKIN_DEF::IFlags, NEUKIN_DEF::INu, NEUKIN_DEF::INuNoOsc, NEUKIN_DEF::IResonance, UtilPDG::isElectron(), UtilPDG::isIon(), UtilPDG::isMuon(), UtilPDG::isNeutrino(), UtilPDG::isTau(), NEUKIN_DEF::IStruckQ, NEUKIN_DEF::ITg, NEUKIN_DEF::P4El1, NEUKIN_DEF::P4El2, NEUKIN_DEF::P4Mu1, NEUKIN_DEF::P4Mu2, NEUKIN_DEF::P4Neu, NEUKIN_DEF::P4NeuNoOsc, NEUKIN_DEF::P4Shw, NEUKIN_DEF::P4Tau, NEUKIN_DEF::P4Tgt, NEUKIN_DEF::Q2, NEUKIN_DEF::Sigma, NEUKIN_DEF::SigmaDiff, NEUKIN_DEF::W2, NEUKIN_DEF::X, NEUKIN_DEF::Y, and NEUKIN_DEF::Z.

Referenced by Get().

00628 {
00629 //======================================================================
00630 // Fill the REROOT_NeuKin part of simrec
00631 //======================================================================
00632 
00633    // New TClonesArray of the right size (1)
00634    TClonesArray* neukins = new TClonesArray("REROOT_NeuKin",1);
00635    neukins->SetName("NeuKinList");
00636    simrec->AdoptComponent(neukins);
00637 
00638    const genie::Interaction*  inter      = gevrec->Summary();
00639    const genie::InitialState& initState  = inter->InitState();
00640    const genie::ProcessInfo&  procInfo   = inter->ProcInfo();
00641    const genie::Kinematics&   kine       = inter->Kine();
00642    const genie::XclsTag&      exclTag    = inter->ExclTag();
00643    //const genie::KPhaseSpace&  phaseSpace = inter->PhaseSpace();
00644 
00645    //int ipdg_tgt = initState.TgtPdg();
00646    //int ia=0, iz=0, ij;
00647    //if ( UtilPDG::isIon(ipdg_tgt) ) UtilPDG::getIonAZJ(ipdg_tgt,ia,iz,ij);
00648 
00649 
00650    // for Daikon (w/ neugen3 v3_5_0) iflags can be interpreted as:
00651    //  0: non-DIS,  1: old KNO hadronization, 2: modified KNO, 3: charm
00652    // 11: JETSET/string frag, 12: JETSET/cluster frag, 13: JETSET/other frag
00653    int iflags = -1;     // generator flags
00654 
00655    double e_em = 0;
00656 
00657    bool saw_mu1, saw_mu2, saw_el1, saw_el2, saw_tau;
00658    TLorentzVector p4shw, p4mu1, p4mu2, p4el1, p4el2, p4tau;
00659    TLorentzVector p4fslepton = kine.FSLeptonP4();
00660    TLorentzVector p4fshadsys = kine.HadSystP4();
00661 
00662    saw_mu1 = saw_mu2 = saw_el1 = saw_el2 = saw_tau = false;
00663 
00664    const TClonesArray* stdhep =
00665      dynamic_cast<const TClonesArray*>
00666         (simrec->FindComponent("TClonesArray","StdHep"));
00667    size_t nhep = stdhep->GetLast()+1;
00668 
00669    // duplicate what was done for GMINOS
00670    for (size_t ihep=0; ihep < nhep; ++ihep) {
00671      TParticle* apart = dynamic_cast<TParticle*>(stdhep->At(ihep));
00672      int ipdg = apart->GetPdgCode();
00673 
00674      // scan for "generator" flags
00675      // PDG codes for PYTHIA/JETSET special particles
00676      if ( genie::kPdgString  == ipdg ) iflags = 11;
00677      if ( genie::kPdgCluster == ipdg ) iflags = 12;
00678      if ( genie::kPdgIndep   == ipdg ) iflags = 13;
00679 
00680      int iapdg = TMath::Abs(ipdg);
00681      // exclude neutrinos from "visible energy sums"
00682      if ( UtilPDG::isNeutrino(iapdg) ) continue;
00683 
00684      TParticlePDG* pdgentry = apart->GetPDG();
00685      double mass   = ( pdgentry ) ? pdgentry->Mass() : 0;
00686      double energy = apart->Energy();
00687      double evis   = energy;
00688      int isbaryon = (ipdg/1000)%10;
00689      // remove mass of baryon or ion
00690      if ( isbaryon || UtilPDG::isIon(iapdg) ) evis -= mass;
00691 
00692      TVector3 partp3(apart->Px(),apart->Py(),apart->Pz());
00693 
00694      // tau "shower" sub-component is a separate vector; it (or its decay
00695      // products) are also included in the shower vector
00696      // taus may have been decayed by the point this routine is called
00697      // so to see see them we'll have to look at non-final state rows as well
00698      if ( UtilPDG::isTau(iapdg) ) {
00699        double signtau = ( ipdg > 0 ) ? -1 : 1;  // tau- or tau+
00700        if ( ! saw_tau || ( energy > TMath::Abs(p4tau.Energy()) ) ) {
00701          saw_tau = true;
00702          p4tau = TLorentzVector(partp3,signtau*energy);                                
00703        }
00704      }
00705 
00706      // nothing further needs non-final state particles
00707      int istatus = apart->GetStatusCode();
00708      if ( istatus != genie::kIStStableFinalState ) continue;
00709 
00710      // always place first muon/electron in first slot (normally primary)
00711      // if second mu or electron then fill that slot as well
00712      // if more than two mu/e, then fill 2nd slot with most energetic 
00713      // muon/electron (beyond 1st), the remainer get added to the 
00714      // shower.  don't worry about small error in shower energy
00715      // by including third muon's mass rather than just ptot
00716 
00717      if ( UtilPDG::isMuon(ipdg) ) {
00718        double signmu =  ( ipdg > 0 ) ? -1 : 1;  // mu- or mu+
00719        if ( ! saw_mu1 ) {
00720          saw_mu1 = true;
00721          p4mu1 = TLorentzVector(partp3,signmu*energy);
00722        } else if ( ! saw_mu2 ) {
00723          saw_mu2 = true;
00724          p4mu2 = TLorentzVector(partp3,signmu*energy);
00725        } else {
00726          if ( energy > TMath::Abs(p4mu2.Energy()) ) {
00727            // replace mu2 entry, add old mu2 to shower (careful about E sign)
00728            p4shw += TLorentzVector(p4mu2.Vect(),TMath::Abs(p4mu2.Energy()));
00729            p4mu2 = TLorentzVector(partp3,signmu*energy);
00730          } else {
00731            // too small, just part of shower
00732            p4shw += TLorentzVector(partp3,evis);
00733          }
00734        }
00735      } else {       
00736        // "shower" -- everything final state that *isn't* a muon
00737        // in GMINOS this all would make more sense because any
00738        // charm particles would have already have been decayed
00739        // in GENIE this isn't true!
00740        p4shw += TLorentzVector(partp3,evis);
00741 
00742        // "EM" component of shower
00743        // 11=electron/positron, 111=pi_0, 22(9)=photon(alt)
00744        if ( iapdg ==   11 || iapdg ==  111 ||
00745             iapdg ==   22 || iapdg ==    9    ) e_em += evis;
00746 
00747        // looking for separate electrons (also included in "shower" vector)
00748        if ( UtilPDG::isElectron(iapdg) ) {
00749          double signe =  ( ipdg > 0 ) ? -1 : 1;  // e- or e+
00750          if ( ! saw_el1 ) {
00751            saw_el1 = true;
00752            p4el1 = TLorentzVector(partp3,signe*energy);
00753          } else if ( ! saw_el2 || energy > TMath::Abs(p4el2.Energy()) ) {
00754            saw_el2 = true;
00755            p4el2 = TLorentzVector(partp3,signe*energy);
00756          }
00757        } // electron vectors
00758 
00759      } // muon
00760    } // loop over all particles
00761 
00762    // other generator flags
00763    if ( procInfo.IsDeepInelastic() ) {
00764      if ( exclTag.IsCharmEvent() ) iflags = 3;
00765      else if ( iflags == -1 )      iflags = 2;
00766    }
00767    else if ( iflags == -1 )        iflags = 0;
00768 
00769    NEUKIN_DEF* nr = new NEUKIN_DEF;
00770    nr->ID          = 0;
00771    nr->INu         = initState.ProbePdg();
00772    nr->INuNoOsc    = nooscpdg;
00773    nr->ITg         = initState.TgtPdg();
00774    nr->IBoson      = -1;
00775    nr->IResonance  = 4000 + procInfo.ScatteringTypeId();
00776    nr->IAction     = 5000 + procInfo.InteractionTypeId();
00777    nr->IStruckQ    = initState.Tgt().HitQrkPdg();
00778    nr->IFlags      = iflags;
00779    nr->A           = initState.Tgt().A();
00780    nr->Z           = initState.Tgt().Z();
00781    nr->Sigma       = gevrec->XSec();
00782    nr->SigmaDiff   = gevrec->DiffXSec();
00783    nr->X           = kine.x(true);
00784    nr->Y           = kine.y(true);
00785    nr->Q2          = kine.Q2(true);
00786    nr->W2          = kine.W(true) * kine.W(true);
00787    nr->EMFrac      = e_em / p4shw.Energy();
00788    for (int j=0; j<4; ++j) {
00789      nr->P4Neu[j]      = (*initState.GetProbeP4(genie::kRfLab))[j];
00790      nr->P4NeuNoOsc[j] = (*initState.GetProbeP4(genie::kRfLab))[j];
00791      nr->P4Tgt[j]      = (*initState.GetTgtP4(genie::kRfLab))[j];
00792      nr->P4Shw[j]      = p4shw[j];
00793      nr->P4Mu1[j]      = p4mu1[j];
00794      nr->P4Mu2[j]      = p4mu2[j];
00795      nr->P4El1[j]      = p4el1[j];
00796      nr->P4El2[j]      = p4el2[j];
00797      nr->P4Tau[j]      = p4tau[j];
00798    }
00799 
00800    // new with placement
00801    Int_t entry = 0;
00802    new((*neukins)[entry])REROOT_NeuKin(nr);  // bogus entry
00803    // fix up neukin elements
00804 
00805    delete nr;
00806 
00807 }

void GenieModule::FillNuEvtKin ( SimSnarlRecord simrec,
genie::EventRecord *  gevrec,
int  nooscpdg 
) [private]

Definition at line 811 of file GenieModule.cxx.

References RecDataRecord< T >::AdoptComponent(), RecDataRecord< T >::FindComponent(), NuEvtKin::kEM, NuEvtKin::kUnknownInteraction, NuEvtKin::kUnknownScatter, Msg::kWarning, NuEvtKin::kWeakCC, NuEvtKin::kWeakMix, NuEvtKin::kWeakNC, NuEvtKin::LinkStdHepIndexToEvent(), and MSG.

Referenced by Get().

00813 {
00814 //======================================================================
00815 // Fill the EventKinematics/NuEvtKin part of simrec
00816 //======================================================================
00817 
00818   // Find the StdHep list
00819   const TClonesArray* stdhep =
00820     dynamic_cast<const TClonesArray*>
00821       (simrec->FindComponent("TClonesArray","StdHep"));
00822   if (!stdhep) {
00823     MSG("EvtKin",Msg::kWarning)
00824       << "FillNuEvtKin found no StdHep" << endl;
00825   }
00826 
00827    // New TClonesArray of the right size (1)
00828    TClonesArray* nuevtkins = new TClonesArray("NuEvtKin",1);
00829    nuevtkins->SetName("NuEvtKinList");
00830    simrec->AdoptComponent(nuevtkins);
00831 
00832    const genie::ProcessInfo&  procInfo = gevrec->Summary()->ProcInfo();
00833 
00834    NuEvtKin::Scatter_t scatter = NuEvtKin::kUnknownScatter;
00835    int genieScatType =  procInfo.ScatteringTypeId();
00836    scatter = (NuEvtKin::Scatter_t)(4000 + genieScatType);
00837 
00838    NuEvtKin::Interaction_t inter = NuEvtKin::kUnknownInteraction;
00839    switch ( procInfo.InteractionTypeId() ) {
00840    case genie::kIntWeakNC:  inter = NuEvtKin::kWeakNC;  break;
00841    case genie::kIntWeakCC:  inter = NuEvtKin::kWeakCC;  break;
00842    case genie::kIntEM:      inter = NuEvtKin::kEM;      break;
00843    case genie::kIntWeakMix: inter = NuEvtKin::kWeakMix; break;
00844    default:                 inter = NuEvtKin::kUnknownInteraction; break;
00845    }
00846 
00847    //const genie::InitialState& initState  = gevrec->Summary()->InitState();
00848    const genie::Kinematics&   kine       = gevrec->Summary()->Kine();
00849 
00850    //Int_t ipdgNu      = initState.ProbePdg();
00851    Int_t ipdgNuNoOsc = nooscpdg;
00852    Float_t x         = kine.x(true);
00853    Float_t y         = kine.y(true);
00854    Float_t Q2        = kine.Q2(true);
00855    Float_t W2        = kine.W(true) * kine.W(true);
00856 
00857 #ifndef BITTER
00858    cerr << "****** GenieModule::FillNuEvtKin skipped new w/placement in TClonesArray " << endl;
00859 #else 
00860    Int_t entry = 0;  // only one entry at this point
00861    // new with placement
00862    new((*nuevtkins)[entry])NuEvtKin(stdhep,ipdgNuNoOsc,
00863                                     scatter,inter,x,y,Q2,W2);
00864 
00865    NuEvtKin* currentNuEvtKin = (NuEvtKin*)((*nuevtkins)[entry]);
00866    // now fix up the entry by telling it what StdHep entries to use
00867    // assume segments of StdHep start with the documentation neutrino entry
00868    Int_t ihep = 0;
00869    Int_t ihep_next_nu = stdhep->GetLast() + 1;
00870    for ( ; ihep < ihep_next_nu ; ++ihep )
00871      currentNuEvtKin->LinkStdHepIndexToEvent(ihep);
00872 
00873    // prepare for next entry
00874    entry++;
00875 #endif
00876 }

void GenieModule::FillStdHep ( SimSnarlRecord simrec,
genie::EventRecord *  gevrec 
) [private]

Definition at line 575 of file GenieModule.cxx.

References RecDataRecord< T >::AdoptComponent(), UtilHepevt::getEmptyStdHep(), gpart, Msg::kDebug, Msg::kInfo, and MSG.

Referenced by Get().

00577 {
00578 //======================================================================
00579 // Fill the StdHep part of simrec
00580 //======================================================================
00581 
00582   TLorentzVector* vertex = gevrec->Vertex();
00583 
00584   //TClonesArray* stdhep_list = fNuTransport.CreateStdHepList();
00585   TClonesArray* stdhepptr = UtilHepevt::getEmptyStdHep();
00586   TClonesArray& stdhep = *stdhepptr;
00587   TIter partitr(gevrec);
00588   genie::GHepParticle* gpart = 0;
00589   // GHepParticles return momentun in units of GeV/c.  Positions are
00590   // in fermi (10^-15m) relative to struck nucleus, i.e.
00591   // need to add vertex pos (given in meters).
00592   int ihep = 0;
00593   while ( ( gpart = dynamic_cast<genie::GHepParticle *>(partitr.Next())) ) {
00594     new(stdhep[ihep])
00595       TParticle(gpart->Pdg(),           gpart->Status(),
00596                 gpart->FirstMother(),   gpart->LastMother(),
00597                 gpart->FirstDaughter(), gpart->LastDaughter(),
00598                 gpart->Px(), gpart->Py(), gpart->Pz(), gpart->E(),
00599                 gpart->Vx(), gpart->Vy(), gpart->Vz(), gpart->Vt());
00600     if ( gpart->Status() == 0 || gpart->Status() == 1 ) {
00601       TParticle* rpart = dynamic_cast<TParticle*>(stdhep[ihep]);
00602       // adjust vertex to global world
00603       const double fm2m = 1.0e-15;
00604       rpart->SetProductionVertex(gpart->Vx()*fm2m + vertex->X(),
00605                                  gpart->Vy()*fm2m + vertex->Y(),
00606                                  gpart->Vz()*fm2m + vertex->Z(),
00607                                  gpart->Vt());
00608     }
00609     ihep++;
00610   }
00611   stdhepptr->SetName("StdHep");
00612   simrec->AdoptComponent(stdhepptr);
00613   MSG("GenieModule",Msg::kInfo)
00614     << " GenieModule evtvtx = [ " 
00615     << vertex->X() << " " 
00616     << vertex->Y() << " " 
00617     << vertex->Z() << " ] " << endl;
00618     
00619   MSG("GenieModule",Msg::kDebug)
00620     << " GenieModule stored StdHep ihep=" << ihep << endl;
00621 
00622 }

JobCResult GenieModule::Get ( MomNavigator mom  )  [virtual]

Implement if your module needs to read data from some external source and fill mom

Reimplemented from JobCModule.

Definition at line 438 of file GenieModule.cxx.

References MomNavigator::AdoptFragment(), dump(), fGFluxI, fGMCJDriver, FillFluxInfo(), FillNeuKin(), FillNuEvtKin(), FillStdHep(), RecDataRecord< T >::FindComponent(), fMaxSnarlDump, fMinSnarlDump, fMonoDx, fMonoDy, fMonoRayOrigin, fPOTsGen, fPOTsLimit, fRun, fRunType, fSnarl, fSpillType, fStdHepDumpMode, fSubRun, fTimeFrame, fTrigGen, fTrigLimit, fTrigSrc, GetNextVldContext(), RecRecordImp< T >::GetTempTags(), gSystem(), hostname, JobCResult::kAOK, Msg::kDebug, JobCResult::kFailed, Msg::kInfo, MSG, Registry::Set(), and SimSnarlHeader::TrimCodename().

00439 {
00440 //======================================================================
00441 // Generate an interaction.
00442 //======================================================================
00443 
00444   if (fTrigGen >= fTrigLimit) {
00445     MSG("GenieModule",Msg::kInfo)
00446       << " GenieModule hit TrigLimit " << fTrigLimit << endl;
00447     raise(SIGTERM);  // signal no more passes through jobc path
00448     return JobCResult::kFailed;
00449   }
00450 
00451   // We need a VldContext
00452   Int_t currentSnarl = fSnarl;
00453   VldContext vldc = GetNextVldContext(true);
00454   MSG("GenieModule",Msg::kDebug) << " GenieModule Reco " << vldc << endl;
00455 
00456   // Start a Record
00457   Int_t errcode = 0;
00458   VldTimeStamp now; // when was this actually generated
00459   std::string  codename = SimSnarlHeader::TrimCodename("$Name:  $");
00460   std::string  hostname(gSystem->HostName());
00461 
00462   SimSnarlHeader simheader(vldc,fRun,fSubRun,fRunType,errcode,
00463                            currentSnarl,fTrigSrc,fTimeFrame,fSpillType,
00464                            now,codename,hostname);
00465   SimSnarlRecord* simrec = new SimSnarlRecord(simheader);
00466   // give the it to MOM to hold as a "fragment"
00467   simrec->GetTempTags().Set("stream","SimSnarl");
00468   simrec->GetTempTags().Set("index",fSnarl-1);
00469   mom->AdoptFragment(simrec);
00470 
00471   genie::flux::GMonoEnergeticFlux* fMonoFlux = 
00472     dynamic_cast<genie::flux::GMonoEnergeticFlux*>(fGFluxI);
00473   if ( fMonoFlux ) {
00474     // fuzz out the position
00475     TVector3& xyz0 = fMonoRayOrigin;
00476     double dx = gRandom->Uniform(-1.,+1.) * fMonoDx;
00477     double dy = gRandom->Uniform(-1.,+1.) * fMonoDy;
00478     fMonoFlux->SetRayOrigin(xyz0.X()+dx,xyz0.Y()+dy,xyz0.Z()); // flux units: meters
00479   }
00480 
00481   // Now start to fill the record
00482   genie::EventRecord* gevrec = fGMCJDriver->GenerateEvent();
00483 
00484   if ( ! gevrec ) {
00485     MSG("GenieModule",Msg::kInfo) << " GenieModule no event generated" << endl;
00486     return JobCResult::kFailed;
00487   }
00488 
00489   // calculate POTs used so far
00490   // the GlobProbScale is supposed to be the interaction probability scale
00491   // used by GENIE (accounts for rejection method of max path then)
00492 
00493   genie::flux::GNuMIFlux* fGNuMIFlux = 
00494     dynamic_cast<genie::flux::GNuMIFlux*>(fGFluxI);
00495   if ( fGNuMIFlux )
00496     fPOTsGen = fGNuMIFlux->UsedPOTs() / fGMCJDriver->GlobProbScale();
00497 
00498 #ifndef MISSING_GSIMPLENTPFLUX
00499   genie::flux::GSimpleNtpFlux* fGSimpleNtpFlux = 
00500     dynamic_cast<genie::flux::GSimpleNtpFlux*>(fGFluxI);
00501   if ( fGSimpleNtpFlux )
00502     fPOTsGen = fGSimpleNtpFlux->UsedPOTs() / fGMCJDriver->GlobProbScale();
00503 #endif
00504 
00505   if (fPOTsGen > fPOTsLimit) {
00506     MSG("GenieModule",Msg::kInfo)
00507       << " GenieModule hit POTsLimit " << fPOTsGen << " > "
00508       << fPOTsLimit << endl;
00509     raise(SIGTERM);  // signal no more passes through jobc path
00510     fPOTsGen = fPOTsLimit;
00511     return JobCResult::kFailed;  // over used POTs, don't accept
00512   }
00513 
00514   fTrigGen++;  // count this as generated
00515 
00516   int nooscpdg = 0; // original flux flavor
00517   if ( fGNuMIFlux ) nooscpdg = fGNuMIFlux->PassThroughInfo().ntype; 
00518 #ifndef MISSING_GSIMPLENTPFLUX
00519   if ( fGSimpleNtpFlux) nooscpdg = fGSimpleNtpFlux->GetCurrentEntry()->pdg;
00520 #endif
00521 
00522   FillStdHep(simrec,gevrec);
00523   FillNeuKin(simrec,gevrec,nooscpdg);
00524   FillNuEvtKin(simrec,gevrec,nooscpdg);
00525   FillFluxInfo(simrec,gevrec);
00526 
00527 #ifdef RWH_DEBUG
00528   MSG("GenieModule",Msg::kInfo)
00529     << "dump snarls " << fMinSnarlDump << " to " << fMaxSnarlDump
00530     << ", current " << currentSnarl << endl;
00531   if (currentSnarl >= fMinSnarlDump && currentSnarl <= fMaxSnarlDump) {
00532     const TClonesArray* stdhep =
00533       dynamic_cast<const TClonesArray*>
00534         (simrec->FindComponent("TClonesArray","StdHep"));
00535     UtilHepevt::dump(stdhep,fStdHepDumpMode);
00536     //was//StdHepUtil::dumpStdHepList(stdhep_list,fStdHepDumpMode);
00537   }
00538   const genie::Interaction*  inter      = gevrec->Summary();
00539   const genie::InitialState& initState  = inter->InitState();
00540   const genie::ProcessInfo&  procInfo   = inter->ProcInfo();
00541   const genie::Kinematics&   kine       = inter->Kine();
00542   const genie::XclsTag&      exclTag    = inter->ExclTag();
00543   const genie::KPhaseSpace&  phaseSpace = inter->PhaseSpace();
00544 
00545   //inter->Print(std::cout);
00546   initState.Print(std::cout);
00547   procInfo.Print(std::cout);
00548   kine.Print(std::cout);
00549   exclTag.Print(std::cout);
00550   phaseSpace.Print("");
00551 #endif
00552 
00553   return JobCResult::kAOK;
00554 }

void GenieModule::GetHallExtent ( double *  xyzmin,
double *  xyzmax,
double &  zlo,
double &  zhi 
) [private]

Definition at line 1141 of file GenieModule.cxx.

References fDetector, UgliGeomHandle::GetHallExtentMax(), UgliGeomHandle::GetHallExtentMin(), GetNextVldContext(), Geo::GetScale(), UgliLoanPool::Instance(), Detector::kNear, Geo::kVMC, and UgliLoanPool::SetUseGeo().

Referenced by CreateGeomVolSelector().

01143 {
01144 
01145   VldContext vldc = GetNextVldContext(false);
01146   // Force GeoGeometry (TGeo) style geometry
01147   UgliLoanPool::Instance()->SetUseGeo(true);
01148   UgliGeomHandle ugh(vldc);
01149 
01150   // get scale from user units to internal ROOT geometry units
01151   double scale = Geo::GetScale(Geo::kVMC); 
01152 
01153   ugh.GetHallExtentMin().GetXYZ(xyzmin);
01154   ugh.GetHallExtentMax().GetXYZ(xyzmax);
01155 
01156   cout << "++++++++++++++ GetHallExtent " << " easy way" << endl;
01157   const char* cxyz = "xyz";
01158   for ( int ic = 0; ic < 3; ic++ ) {
01159     xyzmin[ic] *= scale;
01160     xyzmax[ic] *= scale;
01161     cout << " " << cxyz[ic] << ": " << xyzmin[ic] << " " << xyzmax[ic] << endl;
01162   }
01163 
01164   // MINOS DocDB 3717 would indicate that the NearDet coil extends
01165   // 1.15 m in front, and 0.13 m (seems short) behind the first/last plane
01166   // Don't be too picky about the FarDet, or the downstream side
01167   if ( fDetector == Detector::kNear ) {
01168     zlo = -1.20 * scale;  // give 5 cm extra
01169     zhi = xyzmax[2];
01170   } else {
01171     zlo = -2.00 * scale ;  // two meters upstream
01172     zhi = xyzmax[2]; 
01173   }
01174   cout << "GetHallExtent zlo=" <<  zlo << " zhi=" << zhi << endl;
01175 }

VldContext GenieModule::GetNextVldContext ( bool  advance = true  )  [private]

Definition at line 558 of file GenieModule.cxx.

References VldTimeStamp::Add(), fDetector, fSnarl, fSpillInterval, fSpillTimeStamp, fTimeFrame, and SimFlag::kMC.

Referenced by BeginRun(), Get(), and GetHallExtent().

00559 {
00560 //======================================================================
00561 // Generate a vldcontext
00562 //======================================================================
00563 
00564   VldContext vldc(fDetector,SimFlag::kMC,fSpillTimeStamp);
00565   if ( advance ) {
00566     fSnarl++;
00567     fTimeFrame++;
00568     fSpillTimeStamp.Add(fSpillInterval);
00569   }
00570   return vldc;
00571 }

void GenieModule::Report (  )  [virtual]

Implement to spew end of running report

Reimplemented from JobCModule.

Definition at line 292 of file GenieModule.cxx.

References Msg::kDebug, and MSG.

00293 {
00294     MSG("GenieModule",Msg::kDebug) << " GenieModule Report" << endl;
00295 }

void GenieModule::SetTimeStamp ( std::string  str,
VldTimeStamp ts 
) [private]

Definition at line 1178 of file GenieModule.cxx.

References VldTimeStamp::AsString(), Munits::day, VldTimeStamp::GetDate(), VldTimeStamp::GetTime(), Munits::hour, Msg::kInfo, Msg::kWarning, min, month, MSG, and Munits::year.

Referenced by Config().

01179 {
01180   VldTimeStamp now;
01181   if ( str == "" || str == "now" ) {
01182     ts = now;
01183   } else {
01184 
01185     unsigned int year,month,day,hour,min,sec;
01186     now.GetDate(true,0,&year,&month,&day);
01187     now.GetTime(true,0,&hour,&min,&sec);
01188     
01189     const char* fmt = "%d-%d-%d %d:%d:%d";
01190     int nval = sscanf(str.c_str(),fmt,&year,&month,&day,&hour,&min,&sec);
01191     if ( nval != 6 ) {
01192       MSG("GenieModule",Msg::kWarning) 
01193         << " GenieModule::SetTimeStamp nval=" << nval 
01194         << " str='" << str << "'" << endl;
01195     }
01196     if ( year <  38 ) year += 2000;
01197     if ( year < 100 ) year += 1900;
01198     ts = VldTimeStamp(year,month,day,hour,min,sec);
01199   }
01200   MSG("GenieModule",Msg::kInfo) 
01201     << "GenieModule::SetTimeStamp " << ts.AsString("c") << endl;
01202 }


Member Data Documentation

Definition at line 82 of file GenieModule.h.

Referenced by BeginRun(), Config(), GetHallExtent(), and GetNextVldContext().

genie::GFluxI* GenieModule::fGFluxI [private]

Definition at line 106 of file GenieModule.h.

Referenced by BeginRun(), FillFluxInfo(), and Get().

genie::GMCJDriver* GenieModule::fGMCJDriver [private]

Definition at line 110 of file GenieModule.h.

Referenced by BeginRun(), Get(), and ~GenieModule().

Int_t GenieModule::fMaxSnarlDump [private]

Definition at line 101 of file GenieModule.h.

Referenced by Config(), and Get().

Int_t GenieModule::fMinSnarlDump [private]

Definition at line 100 of file GenieModule.h.

Referenced by Config(), and Get().

double GenieModule::fMonoDx [private]

Definition at line 114 of file GenieModule.h.

Referenced by BeginRun(), and Get().

double GenieModule::fMonoDy [private]

Definition at line 115 of file GenieModule.h.

Referenced by BeginRun(), and Get().

TVector3 GenieModule::fMonoRayOrigin [private]

Definition at line 113 of file GenieModule.h.

Referenced by BeginRun(), and Get().

Double_t GenieModule::fPOTsGen [private]

Definition at line 98 of file GenieModule.h.

Referenced by EndJob(), and Get().

Double_t GenieModule::fPOTsLimit [private]

Definition at line 97 of file GenieModule.h.

Referenced by Config(), and Get().

genie::geometry::ROOTGeomAnalyzer* GenieModule::fROOTGeomAnalyzer [private]

Definition at line 109 of file GenieModule.h.

Referenced by BeginRun(), and CreateGeomVolSelector().

UShort_t GenieModule::fRun [private]

Definition at line 84 of file GenieModule.h.

Referenced by Config(), and Get().

UShort_t GenieModule::fRunType [private]

Definition at line 86 of file GenieModule.h.

Referenced by Config(), and Get().

Int_t GenieModule::fSnarl [private]

Definition at line 90 of file GenieModule.h.

Referenced by Config(), Get(), and GetNextVldContext().

Double_t GenieModule::fSpillInterval [private]

Definition at line 92 of file GenieModule.h.

Referenced by Config(), and GetNextVldContext().

Definition at line 93 of file GenieModule.h.

Referenced by Config(), and GetNextVldContext().

Int_t GenieModule::fSpillType [private]

Definition at line 88 of file GenieModule.h.

Referenced by Config(), and Get().

Definition at line 102 of file GenieModule.h.

Referenced by Config(), and Get().

UShort_t GenieModule::fSubRun [private]

Definition at line 85 of file GenieModule.h.

Referenced by Config(), and Get().

Int_t GenieModule::fTimeFrame [private]

Definition at line 91 of file GenieModule.h.

Referenced by Config(), Get(), and GetNextVldContext().

Int_t GenieModule::fTrigGen [private]

Definition at line 96 of file GenieModule.h.

Referenced by EndJob(), and Get().

Int_t GenieModule::fTrigLimit [private]

Definition at line 95 of file GenieModule.h.

Referenced by Config(), and Get().

UInt_t GenieModule::fTrigSrc [private]

Definition at line 87 of file GenieModule.h.

Referenced by Config(), and Get().


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

Generated on 2 Nov 2017 for loon by  doxygen 1.6.1