LArSoft  v09_93_00
Liquid Argon Software toolkit - https://larsoft.org/
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages
evgb Namespace Reference

Physics generators for neutrinos, cosmic rays, and others. More...

Namespaces

 util
 

Classes

class  CRYHelper
 Interface to the CRY cosmic-ray generator. More...
 
class  EvtTimeFlat
 Flat time distribution. More...
 
class  EvtTimeFNALBeam
 configurable FNAL Beam time distribution More...
 
class  EvtTimeNone
 time distribution that is delta of 0 (no shift) More...
 
class  EvtTimeShiftFactory
 
class  EvtTimeShiftI
 interface for event time distribution More...
 
class  GENIEHelper
 
class  GiBUUHelper
 
class  MCTruthAndFriendsItr
 
class  RNGWrapper
 

Typedefs

typedef evgb::EvtTimeShiftI *(* EvtTimeShiftICtorFuncPtr_t) (const std::string &)
 

Enumerations

enum  { kNeutrinoGenerator = -100, kCosmicRayGenerator = -200 }
 
enum  { kNeutrinoGenerator = -100, kCosmicRayGenerator = -200 }
 

Functions

unsigned int GetRandomNumberSeed ()
 
std::string ExpandEnvVar (const std::string &s)
 
void SetEventGeneratorListAndTune (const std::string &evtlistname="", const std::string &tunename="${GENIE_XSEC_TUNE}")
 
void FillMCTruth (const genie::EventRecord *grec, double spillTime, simb::MCTruth &mctruth, const std::string &genieVersion="unknown", const std::string &genieTune="unknown", bool addGenieVtxTime=false, const std::unordered_map< std::string, std::string > genConfig={})
 
void FillMCTruth (const genie::EventRecord *grec, TLorentzVector &vtxOffset, simb::MCTruth &mctruth, const std::string &genieVersion="unknown", const std::string &genieTune="unknown", bool addGenieVtxTime=false, const std::unordered_map< std::string, std::string > genConfig={})
 
void FillGTruth (const genie::EventRecord *grec, simb::GTruth &gtruth)
 
genie::EventRecord * RetrieveGHEP (const simb::MCTruth &truth, const simb::GTruth &gtruth, bool useFirstTrajPosition=true)
 return genie::EventRecord pointer; callee takes possession More...
 
std::vector< std::unique_ptr< genie::EventRecord > > RetrieveGHEPs (const art::Handle< std::vector< simb::MCTruth >> &mcTruthHandle, const art::Handle< std::vector< simb::GTruth >> &gTruthHandle)
 
std::vector< std::unique_ptr< genie::EventRecord > > RetrieveGHEPs (const art::Event &e, std::string mcTruthLabel, std::string gTruthLabel)
 
void FillMCFlux (genie::GFluxI *fdriver, simb::MCFlux &mcflux)
 
void FillMCFlux (genie::flux::GNuMIFlux *gnumi, simb::MCFlux &mcflux)
 
void FillMCFlux (const genie::flux::GNuMIFluxPassThroughInfo *nflux, double dk2gen, simb::MCFlux &flux)
 
void FillMCFlux (genie::flux::GSimpleNtpFlux *gsimple, simb::MCFlux &mcflux)
 
void FillMCFlux (const genie::flux::GSimpleNtpEntry *nflux_entry, const genie::flux::GSimpleNtpNuMI *nflux_numi, const genie::flux::GSimpleNtpAux *nflux_aux, const genie::flux::GSimpleNtpMeta *nflux_meta, simb::MCFlux &flux)
 
void FillMCFlux (genie::flux::GDk2NuFlux *gdk2nu, simb::MCFlux &mcflux)
 
void FillMCFlux (const bsim::Dk2Nu *dk2nu, const bsim::NuChoice *nuchoice, simb::MCFlux &flux)
 

Variables

template<class T >
double(T::* RNGWrapper )(void)
 
static const int kNue = 0
 
static const int kNueBar = 1
 
static const int kNuMu = 2
 
static const int kNuMuBar = 3
 
static const int kNuTau = 4
 
static const int kNuTauBar = 5
 

Detailed Description

Physics generators for neutrinos, cosmic rays, and others.

Typedef Documentation

typedef evgb::EvtTimeShiftI*(* evgb::EvtTimeShiftICtorFuncPtr_t) (const std::string &)

Definition at line 31 of file EvtTimeShiftFactory.h.

Enumeration Type Documentation

anonymous enum

Enumerate mother codes for primary particles.

Normally the mother code for a primary particle would be set to some arbitrary invalid value like -1, however, we can use this to mark the source of the particle as being either, eg., neutrino induced or from cosmic-rays.

Enumerator
kNeutrinoGenerator 
kCosmicRayGenerator 

Definition at line 14 of file evgenbase.h.

anonymous enum

Enumerate mother codes for primary particles.

Normally the mother code for a primary particle would be set to some arbitrary invalid value like -1, however, we can use this to mark the source of the particle as being either, eg., neutrino induced or from cosmic-rays.

Enumerator
kNeutrinoGenerator 
kCosmicRayGenerator 

Definition at line 14 of file evgenbase.h.

Function Documentation

std::string evgb::ExpandEnvVar ( const std::string &  s)

Definition at line 98 of file GENIE2ART.cxx.

Referenced by SetEventGeneratorListAndTune().

99 {
100  // utility function:
101  // if input "s" starts w/ $, return corresponding env var value,
102  // otherwise return as is
103 
104  if ( s.find('$') != 0 ) return s;
105 
106  // need to remove ${}'s
107  std::string sEnvVar = s;
108  char rmchars[] = "$(){} ";
109  for (unsigned int i = 0; i < strlen(rmchars); ++i) {
110  // remove moves matching characters in [first,last) to end and
111  // returns a past-the-end iterator for the new end of the range [funky!]
112  // erase actually trims the string
113  sEnvVar.erase( std::remove(sEnvVar.begin(), sEnvVar.end(),
114  rmchars[i]), sEnvVar.end() );
115  }
116  const char* charEnvValue = std::getenv(sEnvVar.c_str());
117  if ( ! charEnvValue ) {
118  // resolved into an empty string ... not what one would expect
119 
120  throw cet::exception("UnresolvedEnvVariable")
121  << " can't resolve " << s << " via getenv(\"" << sEnvVar << "\")";
122  return s; // return original (though we won't get here due to throw)
123  }
124  return std::string(charEnvValue);
125 
126 }
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
void evgb::FillGTruth ( const genie::EventRecord *  grec,
simb::GTruth gtruth 
)

Definition at line 391 of file GENIE2ART.cxx.

References simb::GTruth::fCharmHadronPdg, simb::GTruth::fDecayMode, simb::GTruth::fDiffXsec, simb::GTruth::fFinalLeptonPdg, simb::GTruth::fFinalQuarkPdg, simb::GTruth::fFShadSystP4, simb::GTruth::fGint, simb::GTruth::fGPhaseSpace, simb::GTruth::fgQ2, simb::GTruth::fgq2, simb::GTruth::fGscatter, simb::GTruth::fgT, simb::GTruth::fgW, simb::GTruth::fgWrun, simb::GTruth::fgX, simb::GTruth::fgY, simb::GTruth::fHitNucP4, simb::GTruth::fHitNucPos, simb::GTruth::fIsCharm, simb::GTruth::fIsSeaQuark, simb::GTruth::fIsStrange, simb::GTruth::fNumNeutron, simb::GTruth::fNumPi0, simb::GTruth::fNumPiMinus, simb::GTruth::fNumPiPlus, simb::GTruth::fNumProton, simb::GTruth::fNumRho0, simb::GTruth::fNumRhoMinus, simb::GTruth::fNumRhoPlus, simb::GTruth::fNumSingleGammas, simb::GTruth::fprobability, simb::GTruth::fProbeP4, simb::GTruth::fProbePDG, simb::GTruth::fResNum, simb::GTruth::fStrangeHadronPdg, simb::GTruth::ftgtA, simb::GTruth::fTgtP4, simb::GTruth::ftgtPDG, simb::GTruth::ftgtZ, simb::GTruth::fVertex, simb::GTruth::fweight, and simb::GTruth::fXsec.

Referenced by evg::AddGenieEventsToArt::produce(), and evgb::GENIEHelper::Sample().

392  {
393 
394  //interactions info
395  genie::Interaction *inter = record->Summary();
396  const genie::ProcessInfo &procInfo = inter->ProcInfo();
397  truth.fGint = (int)procInfo.InteractionTypeId();
398  truth.fGscatter = (int)procInfo.ScatteringTypeId();
399 
400  //Event info
401  truth.fweight = record->Weight();
402  truth.fprobability = record->Probability();
403  truth.fXsec = record->XSec();
404  truth.fDiffXsec = record->DiffXSec();
405  truth.fGPhaseSpace = (int)record->DiffXSecVars();
406 
407  TLorentzVector vtx;
408  TLorentzVector *erVtx = record->Vertex();
409  vtx.SetXYZT(erVtx->X(), erVtx->Y(), erVtx->Z(), erVtx->T() );
410  truth.fVertex = vtx;
411 
412  //true reaction information and byproducts
413  //(PRE FSI)
414  const genie::XclsTag &exclTag = inter->ExclTag();
415  truth.fIsCharm = exclTag.IsCharmEvent();
416  truth.fCharmHadronPdg = exclTag.CharmHadronPdg();
417  truth.fIsStrange = exclTag.IsStrangeEvent();
418  truth.fStrangeHadronPdg = exclTag.StrangeHadronPdg();
419  truth.fResNum = (int)exclTag.Resonance();
420  truth.fDecayMode = exclTag.DecayMode();
421 
422  truth.fNumPiPlus = truth.fNumPiMinus = truth.fNumPi0 = 0;
423  truth.fNumProton = truth.fNumNeutron = 0;
424  truth.fNumSingleGammas = 0;
425  truth.fNumRho0 = truth.fNumRhoPlus = truth.fNumRhoMinus = 0;
426 
427  //#define FILL_XCLS_OURSELVES
428 #ifndef FILL_XCLS_OURSELVES
429  truth.fNumProton = exclTag.NProtons();
430  truth.fNumNeutron = exclTag.NNeutrons();
431  truth.fNumPi0 = exclTag.NPi0();
432  truth.fNumPiPlus = exclTag.NPiPlus();
433  truth.fNumPiMinus = exclTag.NPiMinus();
434 #if __GENIE_RELEASE_CODE__ >= GRELCODE(3,2,0)
435  truth.fNumSingleGammas = exclTag.NSingleGammas();
436  truth.fNumRho0 = exclTag.NRho0();
437  truth.fNumRhoPlus = exclTag.NRhoPlus();
438  truth.fNumRhoMinus = exclTag.NRhoMinus();
439 #endif
440 #else
441  // try to fill the counts ourselves ...
442  // most events don't have any of these set
443  for (int idx = 0; idx < record->GetEntries(); idx++)
444  {
445  // want hadrons that are about to be sent to the FSI model
446  const genie::GHepParticle * particle = record->Particle(idx);
447  if (particle->Status() != genie::kIStHadronInTheNucleus)
448  continue;
449 
450  int pdg = particle->Pdg();
451  switch ( pdg ) {
452  case genie::kPdgPi0: truth.fNumPi0++; break;
453  case genie::kPdgPiP: truth.fNumPiPlus++; break;
454  case genie::kPdgPiM: truth.fNumPiMinus++; break;
455  case genie::kPdgNeutron: truth.fNumNeutron++; break;
456  case genie::kPdgProton: truth.fNumProton++; break;
457  case genie::kPdgGamma: truth.fNumSingleGammas++; break;
458  case genie::kPdgRho0: truth.fNumRho0++; break;
459  case genie::kPdgRhoP: truth.fNumRhoPlus++; break;
460  case genie::kPdgRhoM: truth.fNumRhoMinus++; break;
461  }
462  /*
463  if (pdg == genie::kPdgPi0)
464  truth.fNumPi0++;
465  else if (pdg == genie::kPdgPiP)
466  truth.fNumPiPlus++;
467  else if (pdg == genie::kPdgPiM)
468  truth.fNumPiMinus++;
469  else if (pdg == genie::kPdgNeutron)
470  truth.fNumNeutron++;
471  else if (pdg == genie::kPdgProton)
472  truth.fNumProton++;
473  else if (pdg == genie::kPdgGamma)
474  truth.fNumSingleGammas++;
475  else if (pdg == genie::kPdgRho0)
476  truth.fNumRho0++;
477  else if (pdg == genie::kPdgRhoP)
478  truth.fNumRhoPlus++;
479  else if (pdg == genie::kPdgRhoM)
480  truth.fNumRhoMinus++;
481  */
482 
483  } // for (idx)
484 #endif
485 
486 #if __GENIE_RELEASE_CODE__ >= GRELCODE(3,2,0)
487  truth.fFinalQuarkPdg = exclTag.FinalQuarkPdg();
488  truth.fFinalLeptonPdg = exclTag.FinalLeptonPdg();
489 #endif
490 
491  // Get the GENIE kinematics info
492  const genie::Kinematics &kine = inter->Kine();
493  // RWH: really should be looping of GENIE Conventions/KineVar_t enum
494  // and only recording/resetting those that were originally there ...
495  truth.fgQ2 = kine.Q2(true);
496  truth.fgq2 = kine.q2(true);
497  truth.fgW = kine.W(true);
498  if ( kine.KVSet(genie::kKVSelt) ) {
499  // only get this if it is set in the Kinematics class
500  // to avoid a warning message
501  truth.fgT = kine.t(true);
502  }
503  truth.fgX = kine.x(true);
504  truth.fgY = kine.y(true);
505  if ( kine.KVSet(genie::kKVW) ) {
506  // only get this if it is set in the Kinematics class
507  // to avoid a warning message
508  truth.fgWrun = kine.W(false);
509  }
510 
511  /*
512  truth.fgQ2 = kine.Q2(false);
513  truth.fgW = kine.W(false);
514  truth.fgT = kine.t(false);
515  truth.fgX = kine.x(false);
516  truth.fgY = kine.y(false);
517  */
518  truth.fFShadSystP4 = kine.HadSystP4();
519 
520  //Initial State info
521  const genie::InitialState &initState = inter->InitState();
522  truth.fProbePDG = initState.ProbePdg();
523  truth.fProbeP4 = *initState.GetProbeP4();
524  truth.fTgtP4 = *initState.GetTgtP4();
525 
526  //Target info
527  const genie::Target &tgt = initState.Tgt();
528  truth.fIsSeaQuark = tgt.HitSeaQrk();
529  truth.fHitNucP4 = tgt.HitNucP4();
530  truth.fHitNucPos = tgt.HitNucPosition();
531  truth.ftgtZ = tgt.Z();
532  truth.ftgtA = tgt.A();
533  truth.ftgtPDG = tgt.Pdg();
534 
535  return;
536 
537 }
void evgb::FillMCFlux ( genie::GFluxI *  fdriver,
simb::MCFlux mcflux 
)

Definition at line 834 of file GENIE2ART.cxx.

Referenced by FillMCFlux(), evg::AddGenieEventsToArt::produce(), and evgb::GENIEHelper::Sample().

835 {
836  // is the real driver hidden behind a blender?
837  genie::flux::GFluxBlender* gblender =
838  dynamic_cast<genie::flux::GFluxBlender *>(fdriver);
839  if ( gblender ) {
840  // it is, it is ... proceed with that ...
841  fdriver = gblender->GetFluxGenerator();
842  }
843 
844  genie::flux::GNuMIFlux* gnumi =
845  dynamic_cast<genie::flux::GNuMIFlux *>(fdriver);
846  if ( gnumi ) {
847  FillMCFlux(gnumi,mcflux);
848  return;
849  } else {
850  genie::flux::GSimpleNtpFlux* gsimple =
851  dynamic_cast<genie::flux::GSimpleNtpFlux *>(fdriver);
852  if ( gsimple ) {
853  FillMCFlux(gsimple,mcflux);
854  return;
855  } else {
856  genie::flux::GDk2NuFlux* gdk2nu =
857  dynamic_cast<genie::flux::GDk2NuFlux *>(fdriver);
858  if ( gdk2nu ) {
859  FillMCFlux(gdk2nu,mcflux);
860  return;
861  } else {
862  static bool first = true;
863  if ( first ) {
864  first = false;
865  std::string dname = typeid(*fdriver).name();
866  // can't use fdriver->GetClass()->GetName(); not derived from TObject
867  mf::LogInfo("GENIE2ART")
868  << " " << __FILE__ << ":" << __LINE__ << "\n"
869  << " no FillMCFlux() for this flux driver: "
870  << dname
871  << " (typeid.name, use \"c++filt -t\" to demangle)"
872  << std::endl;
873  // atmospheric fluxes don't have a method for FillMCFLux
874  // don't abort ... just note the problem, once // abort();
875  }
876  }
877  }
878  }
879 }
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
void FillMCFlux(genie::GFluxI *fdriver, simb::MCFlux &mcflux)
Definition: GENIE2ART.cxx:834
void evgb::FillMCFlux ( genie::flux::GNuMIFlux *  gnumi,
simb::MCFlux mcflux 
)

Definition at line 883 of file GENIE2ART.cxx.

References FillMCFlux().

885 {
886  const genie::flux::GNuMIFluxPassThroughInfo& numiflux =
887  gnumi->PassThroughInfo();
888  const genie::flux::GNuMIFluxPassThroughInfo* nflux = &numiflux;
889  double dk2gen = gnumi->GetDecayDist();
890  evgb::FillMCFlux(nflux,dk2gen,flux);
891 }
void FillMCFlux(genie::GFluxI *fdriver, simb::MCFlux &mcflux)
Definition: GENIE2ART.cxx:834
void evgb::FillMCFlux ( const genie::flux::GNuMIFluxPassThroughInfo *  nflux,
double  dk2gen,
simb::MCFlux flux 
)

Definition at line 892 of file GENIE2ART.cxx.

References simb::MCFlux::fbeampx, simb::MCFlux::fbeampy, simb::MCFlux::fbeampz, simb::MCFlux::fbeamx, simb::MCFlux::fbeamy, simb::MCFlux::fbeamz, simb::MCFlux::fdk2gen, simb::MCFlux::fevtno, simb::MCFlux::fFluxType, simb::MCFlux::fmupare, simb::MCFlux::fmuparpx, simb::MCFlux::fmuparpy, simb::MCFlux::fmuparpz, simb::MCFlux::fndecay, simb::MCFlux::fndxdz, simb::MCFlux::fndxdzfar, simb::MCFlux::fndxdznea, simb::MCFlux::fndydz, simb::MCFlux::fndydzfar, simb::MCFlux::fndydznea, simb::MCFlux::fnecm, simb::MCFlux::fnenergy, simb::MCFlux::fnenergyf, simb::MCFlux::fnenergyn, simb::MCFlux::fnimpwt, simb::MCFlux::fnorig, simb::MCFlux::fnpz, simb::MCFlux::fntype, simb::MCFlux::fnwtfar, simb::MCFlux::fnwtnear, simb::MCFlux::fpdpx, simb::MCFlux::fpdpy, simb::MCFlux::fpdpz, simb::MCFlux::fppdxdz, simb::MCFlux::fppdydz, simb::MCFlux::fppenergy, simb::MCFlux::fppmedium, simb::MCFlux::fpppz, simb::MCFlux::fppvx, simb::MCFlux::fppvy, simb::MCFlux::fppvz, simb::MCFlux::fptype, simb::MCFlux::frun, simb::MCFlux::ftgen, simb::MCFlux::ftgppx, simb::MCFlux::ftgppy, simb::MCFlux::ftgppz, simb::MCFlux::ftgptype, simb::MCFlux::ftprivx, simb::MCFlux::ftprivy, simb::MCFlux::ftprivz, simb::MCFlux::ftptype, simb::MCFlux::ftpx, simb::MCFlux::ftpy, simb::MCFlux::ftpz, simb::MCFlux::ftvx, simb::MCFlux::ftvy, simb::MCFlux::ftvz, simb::MCFlux::fvx, simb::MCFlux::fvy, simb::MCFlux::fvz, simb::MCFlux::fxpoint, simb::MCFlux::fypoint, simb::MCFlux::fzpoint, simb::kNtuple, and simb::MCFlux::Reset().

895 {
896  flux.Reset();
897  flux.fFluxType = simb::kNtuple;
898 
899  // check the particle codes and the units passed through
900  // nflux->pcodes: 0=original GEANT particle codes, 1=converted to PDG
901  // nflux->units: 0=original GEANT cm, 1=meters
902  if (nflux->pcodes != 1 && nflux->units != 0) {
903  mf::LogError("FillMCFlux")
904  << "either wrong particle codes or units "
905  << "from flux object - beware!!";
906  }
907 
908  // maintained variable names from gnumi ntuples
909  // see http://www.hep.utexas.edu/~zarko/wwwgnumi/v19/[/v19/output_gnumi.html]
910 
911  flux.frun = nflux->run;
912  flux.fevtno = nflux->evtno;
913  flux.fndxdz = nflux->ndxdz;
914  flux.fndydz = nflux->ndydz;
915  flux.fnpz = nflux->npz;
916  flux.fnenergy = nflux->nenergy;
917  flux.fndxdznea = nflux->ndxdznea;
918  flux.fndydznea = nflux->ndydznea;
919  flux.fnenergyn = nflux->nenergyn;
920  flux.fnwtnear = nflux->nwtnear;
921  flux.fndxdzfar = nflux->ndxdzfar;
922  flux.fndydzfar = nflux->ndydzfar;
923  flux.fnenergyf = nflux->nenergyf;
924  flux.fnwtfar = nflux->nwtfar;
925  flux.fnorig = nflux->norig;
926  flux.fndecay = nflux->ndecay;
927  flux.fntype = nflux->ntype;
928  flux.fvx = nflux->vx;
929  flux.fvy = nflux->vy;
930  flux.fvz = nflux->vz;
931  flux.fpdpx = nflux->pdpx;
932  flux.fpdpy = nflux->pdpy;
933  flux.fpdpz = nflux->pdpz;
934  flux.fppdxdz = nflux->ppdxdz;
935  flux.fppdydz = nflux->ppdydz;
936  flux.fpppz = nflux->pppz;
937  flux.fppenergy = nflux->ppenergy;
938  flux.fppmedium = nflux->ppmedium;
939  flux.fptype = nflux->ptype; // converted to PDG
940  flux.fppvx = nflux->ppvx;
941  flux.fppvy = nflux->ppvy;
942  flux.fppvz = nflux->ppvz;
943  flux.fmuparpx = nflux->muparpx;
944  flux.fmuparpy = nflux->muparpy;
945  flux.fmuparpz = nflux->muparpz;
946  flux.fmupare = nflux->mupare;
947  flux.fnecm = nflux->necm;
948  flux.fnimpwt = nflux->nimpwt;
949  flux.fxpoint = nflux->xpoint;
950  flux.fypoint = nflux->ypoint;
951  flux.fzpoint = nflux->zpoint;
952  flux.ftvx = nflux->tvx;
953  flux.ftvy = nflux->tvy;
954  flux.ftvz = nflux->tvz;
955  flux.ftpx = nflux->tpx;
956  flux.ftpy = nflux->tpy;
957  flux.ftpz = nflux->tpz;
958  flux.ftptype = nflux->tptype; // converted to PDG
959  flux.ftgen = nflux->tgen;
960  flux.ftgptype = nflux->tgptype; // converted to PDG
961  flux.ftgppx = nflux->tgppx;
962  flux.ftgppy = nflux->tgppy;
963  flux.ftgppz = nflux->tgppz;
964  flux.ftprivx = nflux->tprivx;
965  flux.ftprivy = nflux->tprivy;
966  flux.ftprivz = nflux->tprivz;
967  flux.fbeamx = nflux->beamx;
968  flux.fbeamy = nflux->beamy;
969  flux.fbeamz = nflux->beamz;
970  flux.fbeampx = nflux->beampx;
971  flux.fbeampy = nflux->beampy;
972  flux.fbeampz = nflux->beampz;
973 
974  flux.fdk2gen = dk2gen;
975 
976  return;
977 }
int frun
Definition: MCFlux.h:35
double fnenergy
Definition: MCFlux.h:40
double ftgppy
Definition: MCFlux.h:86
double fpdpx
Definition: MCFlux.h:55
double ftpx
Definition: MCFlux.h:79
double ftprivy
Definition: MCFlux.h:89
Full flux simulation ntuple.
Definition: MCFlux.h:21
int fppmedium
Definition: MCFlux.h:62
double fbeamx
Definition: MCFlux.h:91
double fbeampy
Definition: MCFlux.h:95
int ftptype
Definition: MCFlux.h:82
simb::flux_code_ fFluxType
Definition: MCFlux.h:98
double fvx
Definition: MCFlux.h:52
double fnwtfar
Definition: MCFlux.h:48
double fppdxdz
Definition: MCFlux.h:58
int ftgen
Definition: MCFlux.h:83
MaybeLogger_< ELseverityLevel::ELsev_error, false > LogError
double fdk2gen
distance from decay to ray origin
Definition: MCFlux.h:103
double ftprivx
Definition: MCFlux.h:88
double fndydzfar
Definition: MCFlux.h:46
int ftgptype
Definition: MCFlux.h:84
double fnenergyf
Definition: MCFlux.h:47
void Reset()
Definition: MCFlux.cxx:94
double ftprivz
Definition: MCFlux.h:90
double fbeamz
Definition: MCFlux.h:93
double fnwtnear
Definition: MCFlux.h:44
int fptype
Definition: MCFlux.h:63
double ftvx
Definition: MCFlux.h:76
int fndecay
Definition: MCFlux.h:50
double fndydz
Definition: MCFlux.h:38
int fnorig
Definition: MCFlux.h:49
double fbeamy
Definition: MCFlux.h:92
double ftgppz
Definition: MCFlux.h:87
double fpdpz
Definition: MCFlux.h:57
double fmuparpz
Definition: MCFlux.h:69
double ftgppx
Definition: MCFlux.h:85
int fevtno
Definition: MCFlux.h:36
double fppdydz
Definition: MCFlux.h:59
double ftpz
Definition: MCFlux.h:81
double ftvy
Definition: MCFlux.h:77
double fbeampz
Definition: MCFlux.h:96
double fndxdzfar
Definition: MCFlux.h:45
double fpdpy
Definition: MCFlux.h:56
double fnpz
Definition: MCFlux.h:39
double fzpoint
Definition: MCFlux.h:75
double fvy
Definition: MCFlux.h:53
double fmupare
Definition: MCFlux.h:70
double fxpoint
Definition: MCFlux.h:73
double fndxdznea
Definition: MCFlux.h:41
double fppenergy
Definition: MCFlux.h:61
int fntype
Definition: MCFlux.h:51
double fnecm
Definition: MCFlux.h:71
double fndydznea
Definition: MCFlux.h:42
double fypoint
Definition: MCFlux.h:74
double fmuparpy
Definition: MCFlux.h:68
double fppvx
Definition: MCFlux.h:64
double fpppz
Definition: MCFlux.h:60
double fppvy
Definition: MCFlux.h:65
double fbeampx
Definition: MCFlux.h:94
double fppvz
Definition: MCFlux.h:66
double ftpy
Definition: MCFlux.h:80
double fndxdz
Definition: MCFlux.h:37
double fvz
Definition: MCFlux.h:54
double fmuparpx
Definition: MCFlux.h:67
double fnenergyn
Definition: MCFlux.h:43
double ftvz
Definition: MCFlux.h:78
double fnimpwt
Definition: MCFlux.h:72
void evgb::FillMCFlux ( genie::flux::GSimpleNtpFlux *  gsimple,
simb::MCFlux mcflux 
)

Definition at line 980 of file GENIE2ART.cxx.

References FillMCFlux().

982 {
983  const genie::flux::GSimpleNtpEntry* nflux_entry =
984  gsimple->GetCurrentEntry();
985  const genie::flux::GSimpleNtpNuMI* nflux_numi =
986  gsimple->GetCurrentNuMI();
987  const genie::flux::GSimpleNtpAux* nflux_aux =
988  gsimple->GetCurrentAux();
989  const genie::flux::GSimpleNtpMeta* nflux_meta =
990  gsimple->GetCurrentMeta();
991  evgb::FillMCFlux(nflux_entry, nflux_numi, nflux_aux, nflux_meta, flux);
992 }
void FillMCFlux(genie::GFluxI *fdriver, simb::MCFlux &mcflux)
Definition: GENIE2ART.cxx:834
void evgb::FillMCFlux ( const genie::flux::GSimpleNtpEntry *  nflux_entry,
const genie::flux::GSimpleNtpNuMI *  nflux_numi,
const genie::flux::GSimpleNtpAux *  nflux_aux,
const genie::flux::GSimpleNtpMeta *  nflux_meta,
simb::MCFlux flux 
)

Definition at line 993 of file GENIE2ART.cxx.

References simb::MCFlux::fdk2gen, simb::MCFlux::fevtno, simb::MCFlux::fFluxType, simb::MCFlux::fmupare, simb::MCFlux::fmuparpx, simb::MCFlux::fmuparpy, simb::MCFlux::fmuparpz, simb::MCFlux::fndecay, simb::MCFlux::fnecm, simb::MCFlux::fnenergyf, simb::MCFlux::fnenergyn, simb::MCFlux::fnimpwt, simb::MCFlux::fntype, simb::MCFlux::fnwtfar, simb::MCFlux::fnwtnear, simb::MCFlux::fpdpx, simb::MCFlux::fpdpy, simb::MCFlux::fpdpz, simb::MCFlux::fppdxdz, simb::MCFlux::fppdydz, simb::MCFlux::fppmedium, simb::MCFlux::fpppz, simb::MCFlux::fptype, simb::MCFlux::frun, simb::MCFlux::ftgen, simb::MCFlux::ftgptype, simb::MCFlux::ftptype, simb::MCFlux::ftpx, simb::MCFlux::ftpy, simb::MCFlux::ftpz, simb::MCFlux::fvx, simb::MCFlux::fvy, simb::MCFlux::fvz, simb::kSimple_Flux, and simb::MCFlux::Reset().

998 {
999  flux.Reset();
1001 
1002  // maintained variable names from gnumi ntuples
1003  // see http://www.hep.utexas.edu/~zarko/wwwgnumi/v19/[/v19/output_gnumi.html]
1004 
1005 
1006  flux.fntype = nflux_entry->pdg;
1007  flux.fnimpwt = nflux_entry->wgt;
1008  flux.fdk2gen = nflux_entry->dist;
1009  flux.fnenergyn = flux.fnenergyf = nflux_entry->E;
1010 
1011  if ( nflux_numi ) {
1012  flux.frun = nflux_numi->run;
1013  flux.fevtno = nflux_numi->evtno;
1014  flux.ftpx = nflux_numi->tpx;
1015  flux.ftpy = nflux_numi->tpy;
1016  flux.ftpz = nflux_numi->tpz;
1017  flux.ftptype = nflux_numi->tptype; // converted to PDG
1018  flux.fvx = nflux_numi->vx;
1019  flux.fvy = nflux_numi->vy;
1020  flux.fvz = nflux_numi->vz;
1021 
1022  flux.fndecay = nflux_numi->ndecay;
1023  flux.fppmedium = nflux_numi->ppmedium;
1024 
1025  flux.fpdpx = nflux_numi->pdpx;
1026  flux.fpdpy = nflux_numi->pdpy;
1027  flux.fpdpz = nflux_numi->pdpz;
1028 
1029  double apppz = nflux_numi->pppz;
1030  if ( TMath::Abs(nflux_numi->pppz) < 1.0e-30 ) apppz = 1.0e-30;
1031  flux.fppdxdz = nflux_numi->pppx / apppz;
1032  flux.fppdydz = nflux_numi->pppy / apppz;
1033  flux.fpppz = nflux_numi->pppz;
1034 
1035  flux.fptype = nflux_numi->ptype;
1036 
1037  }
1038 
1039  // anything useful stuffed into vdbl or vint?
1040  // need to check the metadata auxintname, auxdblname
1041 
1042  if ( nflux_aux && nflux_meta ) {
1043 
1044  // references just for reducing complexity
1045  const std::vector<std::string>& auxdblname = nflux_meta->auxdblname;
1046  const std::vector<std::string>& auxintname = nflux_meta->auxintname;
1047  const std::vector<int>& auxint = nflux_aux->auxint;
1048  const std::vector<double>& auxdbl = nflux_aux->auxdbl;
1049 
1050  for (size_t id=0; id<auxdblname.size(); ++id) {
1051  if ("muparpx" == auxdblname[id]) flux.fmuparpx = auxdbl[id];
1052  if ("muparpy" == auxdblname[id]) flux.fmuparpy = auxdbl[id];
1053  if ("muparpz" == auxdblname[id]) flux.fmuparpz = auxdbl[id];
1054  if ("mupare" == auxdblname[id]) flux.fmupare = auxdbl[id];
1055  if ("necm" == auxdblname[id]) flux.fnecm = auxdbl[id];
1056  if ("nimpwt" == auxdblname[id]) flux.fnimpwt = auxdbl[id];
1057  if ("fgXYWgt" == auxdblname[id]) {
1058  flux.fnwtnear = flux.fnwtfar = auxdbl[id];
1059  }
1060  }
1061  for (size_t ii=0; ii<auxintname.size(); ++ii) {
1062  if ("tgen" == auxintname[ii]) flux.ftgen = auxint[ii];
1063  if ("tgptype" == auxintname[ii]) flux.ftgptype = auxint[ii];
1064  }
1065 
1066  }
1067 
1068 //#define RWH_TEST
1069 #ifdef RWH_TEST
1070  static bool first = true;
1071  if (first) {
1072  first = false;
1073  mf::LogDebug("GENIE2ART")
1074  << __FILE__ << ":" << __LINE__
1075  << " one time dump of GSimple objects\n";
1076  if ( nflux_meta ) {
1077  mf::LogDebug("GENIE2ART")
1078  << "evgb::FillMCFlux() GSimpleNtpMeta:\n"
1079  << *nflux_meta << "\n";
1080  } else {
1081  mf::LogDebug("GENIE2ART")
1082  << "evgb::FillMCFlux() no GSimpleNtpMeta:\n";
1083  }
1084  }
1085  //mf::LogDebug("GENIEHelper")
1086  mf::LogDebug("GENIE2ART")
1087  << "simb::MCFlux:\n"
1088  << flux << "\n"
1089  << "GSimpleNtpFlux:\n";
1090  if ( nflux_entry) mf::LogDebug("GENIE2ART") << *nflux_entry << "\n";
1091  else mf::LogDebug("GENIE2ART") << "no GSimpleNtpEntry\n";
1092  if ( nflux_numi ) mf::LogDebug("GENIE2ART") << *nflux_numi << "\n";
1093  else mf::LogDebug("GENIE2ART") << "no GSimpleNtpNuMI\n";
1094  if ( nflux_aux ) mf::LogDebug("GENIE2ART") << *nflux_aux << "\n";
1095  else mf::LogDebug("GENIE2ART") << "no GSimpleNtpAux\n";
1096 #endif
1097 
1098  // flux.fndxdz = nflux.ndxdz;
1099  // flux.fndydz = nflux.ndydz;
1100  // flux.fnpz = nflux.npz;
1101  // flux.fnenergy = nflux.nenergy;
1102  // flux.fndxdznea = nflux.ndxdznea;
1103  // flux.fndydznea = nflux.ndydznea;
1104  // flux.fnenergyn = nflux.nenergyn;
1105  // flux.fnwtnear = nflux.nwtnear;
1106  // flux.fndxdzfar = nflux.ndxdzfar;
1107  // flux.fndydzfar = nflux.ndydzfar;
1108  // flux.fnenergyf = nflux.nenergyf;
1109  // flux.fnwtfar = nflux.nwtfar;
1110  // flux.fnorig = nflux.norig;
1111  // in numi // flux.fndecay = nflux.ndecay;
1112  // flux.fntype = nflux.ntype;
1113  // in numi // flux.fvx = nflux.vx;
1114  // in numi // flux.fvy = nflux.vy;
1115  // in numi // flux.fvz = nflux.vz;
1116  // flux.fppenergy = nflux.ppenergy;
1117  // in numi // flux.fppmedium = nflux.ppmedium;
1118  // flux.fppvx = nflux.ppvx;
1119  // flux.fppvy = nflux.ppvy;
1120  // flux.fppvz = nflux.ppvz;
1121  // see above // flux.fmuparpx = nflux.muparpx;
1122  // see above // flux.fmuparpy = nflux.muparpy;
1123  // see above // flux.fmuparpz = nflux.muparpz;
1124  // see above // flux.fmupare = nflux.mupare;
1125  // see above // flux.fnecm = nflux.necm;
1126  // see above // flux.fnimpwt = nflux.nimpwt;
1127  // flux.fxpoint = nflux.xpoint;
1128  // flux.fypoint = nflux.ypoint;
1129  // flux.fzpoint = nflux.zpoint;
1130  // flux.ftvx = nflux.tvx;
1131  // flux.ftvy = nflux.tvy;
1132  // flux.ftvz = nflux.tvz;
1133  // see above // flux.ftgen = nflux.tgen;
1134  // see above // flux.ftgptype = nflux.tgptype; // converted to PDG
1135  // flux.ftgppx = nflux.tgppx;
1136  // flux.ftgppy = nflux.tgppy;
1137  // flux.ftgppz = nflux.tgppz;
1138  // flux.ftprivx = nflux.tprivx;
1139  // flux.ftprivy = nflux.tprivy;
1140  // flux.ftprivz = nflux.tprivz;
1141  // flux.fbeamx = nflux.beamx;
1142  // flux.fbeamy = nflux.beamy;
1143  // flux.fbeamz = nflux.beamz;
1144  // flux.fbeampx = nflux.beampx;
1145  // flux.fbeampy = nflux.beampy;
1146  // flux.fbeampz = nflux.beampz;
1147 
1148  return;
1149 }
int frun
Definition: MCFlux.h:35
double fpdpx
Definition: MCFlux.h:55
double ftpx
Definition: MCFlux.h:79
int fppmedium
Definition: MCFlux.h:62
int ftptype
Definition: MCFlux.h:82
simb::flux_code_ fFluxType
Definition: MCFlux.h:98
double fvx
Definition: MCFlux.h:52
double fnwtfar
Definition: MCFlux.h:48
double fppdxdz
Definition: MCFlux.h:58
int ftgen
Definition: MCFlux.h:83
double fdk2gen
distance from decay to ray origin
Definition: MCFlux.h:103
int ftgptype
Definition: MCFlux.h:84
double fnenergyf
Definition: MCFlux.h:47
void Reset()
Definition: MCFlux.cxx:94
double fnwtnear
Definition: MCFlux.h:44
int fptype
Definition: MCFlux.h:63
int fndecay
Definition: MCFlux.h:50
double fpdpz
Definition: MCFlux.h:57
double fmuparpz
Definition: MCFlux.h:69
int fevtno
Definition: MCFlux.h:36
double fppdydz
Definition: MCFlux.h:59
double ftpz
Definition: MCFlux.h:81
double fpdpy
Definition: MCFlux.h:56
double fvy
Definition: MCFlux.h:53
double fmupare
Definition: MCFlux.h:70
int fntype
Definition: MCFlux.h:51
MaybeLogger_< ELseverityLevel::ELsev_success, false > LogDebug
double fnecm
Definition: MCFlux.h:71
double fmuparpy
Definition: MCFlux.h:68
double fpppz
Definition: MCFlux.h:60
double ftpy
Definition: MCFlux.h:80
double fvz
Definition: MCFlux.h:54
double fmuparpx
Definition: MCFlux.h:67
A simplified flux ntuple for quick running.
Definition: MCFlux.h:22
double fnenergyn
Definition: MCFlux.h:43
double fnimpwt
Definition: MCFlux.h:72
void evgb::FillMCFlux ( genie::flux::GDk2NuFlux *  gdk2nu,
simb::MCFlux mcflux 
)

Definition at line 1152 of file GENIE2ART.cxx.

References simb::MCFlux::fdk2gen, and FillMCFlux().

1154 {
1155  const bsim::Dk2Nu& dk2nu = gdk2nu->GetDk2Nu();
1156  const bsim::NuChoice& nuchoice = gdk2nu->GetNuChoice();
1157  evgb::FillMCFlux(&dk2nu,&nuchoice,flux);
1158  // do this after Fill as that does a Reset()
1159  flux.fdk2gen = gdk2nu->GetDecayDist();
1160 }
void FillMCFlux(genie::GFluxI *fdriver, simb::MCFlux &mcflux)
Definition: GENIE2ART.cxx:834
void evgb::FillMCFlux ( const bsim::Dk2Nu *  dk2nu,
const bsim::NuChoice *  nuchoice,
simb::MCFlux flux 
)

Definition at line 1161 of file GENIE2ART.cxx.

References simb::MCFlux::fevtno, simb::MCFlux::fFluxType, simb::MCFlux::fmupare, simb::MCFlux::fmuparpx, simb::MCFlux::fmuparpy, simb::MCFlux::fmuparpz, simb::MCFlux::fndecay, simb::MCFlux::fnecm, simb::MCFlux::fnenergyf, simb::MCFlux::fnenergyn, simb::MCFlux::fnimpwt, simb::MCFlux::fnorig, simb::MCFlux::fntype, simb::MCFlux::fnwtfar, simb::MCFlux::fnwtnear, simb::MCFlux::fpdpx, simb::MCFlux::fpdpy, simb::MCFlux::fpdpz, simb::MCFlux::fppdxdz, simb::MCFlux::fppdydz, simb::MCFlux::fppenergy, simb::MCFlux::fppmedium, simb::MCFlux::fpppz, simb::MCFlux::fppvx, simb::MCFlux::fppvy, simb::MCFlux::fppvz, simb::MCFlux::fptype, simb::MCFlux::frun, simb::MCFlux::ftgen, simb::MCFlux::ftptype, simb::MCFlux::ftpx, simb::MCFlux::ftpy, simb::MCFlux::ftpz, simb::MCFlux::ftvx, simb::MCFlux::ftvy, simb::MCFlux::ftvz, simb::MCFlux::fvx, simb::MCFlux::fvy, simb::MCFlux::fvz, simb::kDk2Nu, and simb::MCFlux::Reset().

1164 {
1165  flux.Reset();
1166  flux.fFluxType = simb::kDk2Nu;
1167 
1168  if ( dk2nu ) {
1169  flux.frun = dk2nu->job;
1170  flux.fevtno = dk2nu->potnum;
1171 
1172  // ignore vector<bsim::NuRay> (see nuchoice above)
1173 
1174  // bsim::Decay object
1175  flux.fnorig = dk2nu->decay.norig;
1176  flux.fndecay = dk2nu->decay.ndecay;
1177  flux.fntype = dk2nu->decay.ntype;
1178  flux.fppmedium = dk2nu->decay.ppmedium;
1179  flux.fptype = dk2nu->decay.ptype;
1180 
1181  flux.fvx = dk2nu->decay.vx;
1182  flux.fvy = dk2nu->decay.vy;
1183  flux.fvz = dk2nu->decay.vz;
1184  flux.fpdpx = dk2nu->decay.pdpx;
1185  flux.fpdpy = dk2nu->decay.pdpy;
1186  flux.fpdpz = dk2nu->decay.pdpz;
1187 
1188  flux.fppdxdz = dk2nu->decay.ppdxdz;
1189  flux.fppdydz = dk2nu->decay.ppdydz;
1190  flux.fpppz = dk2nu->decay.pppz;
1191  flux.fppenergy = dk2nu->decay.ppenergy;
1192 
1193  flux.fmuparpx = dk2nu->decay.muparpx;
1194  flux.fmuparpy = dk2nu->decay.muparpy;
1195  flux.fmuparpz = dk2nu->decay.muparpz;
1196  flux.fmupare = dk2nu->decay.mupare;
1197 
1198  flux.fnecm = dk2nu->decay.necm;
1199  flux.fnimpwt = dk2nu->decay.nimpwt;
1200 
1201  // no place for: vector<bsim::Ancestor>
1202 
1203  // production vertex of nu parent
1204  flux.fppvx = dk2nu->ppvx;
1205  flux.fppvy = dk2nu->ppvy;
1206  flux.fppvz = dk2nu->ppvz;
1207 
1208  // bsim::TgtExit object
1209  flux.ftvx = dk2nu->tgtexit.tvx;
1210  flux.ftvy = dk2nu->tgtexit.tvy;
1211  flux.ftvz = dk2nu->tgtexit.tvz;
1212  flux.ftpx = dk2nu->tgtexit.tpx;
1213  flux.ftpy = dk2nu->tgtexit.tpy;
1214  flux.ftpz = dk2nu->tgtexit.tpz;
1215  flux.ftptype = dk2nu->tgtexit.tptype; // converted to PDG
1216  flux.ftgen = dk2nu->tgtexit.tgen;
1217 
1218  // ignore vector<bsim::Traj>
1219 
1220  }
1221 
1222  if ( nuchoice ) {
1223  flux.fntype = nuchoice->pdgNu;
1224  flux.fnimpwt = nuchoice->impWgt;
1225 
1226  flux.fnenergyn = flux.fnenergyf = nuchoice->p4NuUser.E();
1227  flux.fnwtnear = flux.fnwtfar = nuchoice->xyWgt;
1228  }
1229 
1230  /*
1231  // anything useful stuffed into vdbl or vint?
1232  // need to check the metadata auxintname, auxdblname
1233 
1234  if ( nflux_aux && nflux_meta ) {
1235 
1236  // references just for reducing complexity
1237  const std::vector<std::string>& auxdblname = nflux_meta->auxdblname;
1238  const std::vector<std::string>& auxintname = nflux_meta->auxintname;
1239  const std::vector<int>& auxint = nflux_aux->auxint;
1240  const std::vector<double>& auxdbl = nflux_aux->auxdbl;
1241 
1242  for (size_t id=0; id<auxdblname.size(); ++id) {
1243  if ("muparpx" == auxdblname[id]) flux.fmuparpx = auxdbl[id];
1244  if ("muparpy" == auxdblname[id]) flux.fmuparpy = auxdbl[id];
1245  if ("muparpz" == auxdblname[id]) flux.fmuparpz = auxdbl[id];
1246  if ("mupare" == auxdblname[id]) flux.fmupare = auxdbl[id];
1247  if ("necm" == auxdblname[id]) flux.fnecm = auxdbl[id];
1248  if ("nimpwt" == auxdblname[id]) flux.fnimpwt = auxdbl[id];
1249  if ("fgXYWgt" == auxdblname[id]) {
1250  flux.fnwtnear = flux.fnwtfar = auxdbl[id];
1251  }
1252  }
1253  for (size_t ii=0; ii<auxintname.size(); ++ii) {
1254  if ("tgen" == auxintname[ii]) flux.ftgen = auxint[ii];
1255  if ("tgptype" == auxintname[ii]) flux.ftgptype = auxint[ii];
1256  }
1257 
1258  }
1259 
1260  // probably can get this from vx,vy,vz + NuChoice
1261  flux.fdk2gen = gdk2nu->GetDecayDist();
1262 
1263  */
1264 
1265  // flux.fndxdz = nflux.ndxdz;
1266  // flux.fndydz = nflux.ndydz;
1267  // flux.fnpz = nflux.npz;
1268  // flux.fnenergy = nflux.nenergy;
1269  // flux.fndxdznea = nflux.ndxdznea;
1270  // flux.fndydznea = nflux.ndydznea;
1271  // flux.fnenergyn = nflux.nenergyn;
1272  // flux.fnwtnear = nflux.nwtnear;
1273  // flux.fndxdzfar = nflux.ndxdzfar;
1274  // flux.fndydzfar = nflux.ndydzfar;
1275  // flux.fnenergyf = nflux.nenergyf;
1276  // flux.fnwtfar = nflux.nwtfar;
1277  // flux.fnorig = nflux.norig;
1278  // in numi // flux.fndecay = nflux.ndecay;
1279  // flux.fntype = nflux.ntype;
1280  // in numi // flux.fvx = nflux.vx;
1281  // in numi // flux.fvy = nflux.vy;
1282  // in numi // flux.fvz = nflux.vz;
1283  // flux.fppenergy = nflux.ppenergy;
1284  // in numi // flux.fppmedium = nflux.ppmedium;
1285  // flux.fppvx = nflux.ppvx;
1286  // flux.fppvy = nflux.ppvy;
1287  // flux.fppvz = nflux.ppvz;
1288  // see above // flux.fmuparpx = nflux.muparpx;
1289  // see above // flux.fmuparpy = nflux.muparpy;
1290  // see above // flux.fmuparpz = nflux.muparpz;
1291  // see above // flux.fmupare = nflux.mupare;
1292  // see above // flux.fnecm = nflux.necm;
1293  // see above // flux.fnimpwt = nflux.nimpwt;
1294  // flux.fxpoint = nflux.xpoint;
1295  // flux.fypoint = nflux.ypoint;
1296  // flux.fzpoint = nflux.zpoint;
1297  // flux.ftvx = nflux.tvx;
1298  // flux.ftvy = nflux.tvy;
1299  // flux.ftvz = nflux.tvz;
1300  // see above // flux.ftgen = nflux.tgen;
1301  // see above // flux.ftgptype = nflux.tgptype; // converted to PDG
1302  // flux.ftgppx = nflux.tgppx;
1303  // flux.ftgppy = nflux.tgppy;
1304  // flux.ftgppz = nflux.tgppz;
1305  // flux.ftprivx = nflux.tprivx;
1306  // flux.ftprivy = nflux.tprivy;
1307  // flux.ftprivz = nflux.tprivz;
1308  // flux.fbeamx = nflux.beamx;
1309  // flux.fbeamy = nflux.beamy;
1310  // flux.fbeamz = nflux.beamz;
1311  // flux.fbeampx = nflux.beampx;
1312  // flux.fbeampy = nflux.beampy;
1313  // flux.fbeampz = nflux.beampz;
1314 
1315  return;
1316 }
Unified ntuple flux format (replaces 2)
Definition: MCFlux.h:23
int frun
Definition: MCFlux.h:35
double fpdpx
Definition: MCFlux.h:55
double ftpx
Definition: MCFlux.h:79
int fppmedium
Definition: MCFlux.h:62
int ftptype
Definition: MCFlux.h:82
simb::flux_code_ fFluxType
Definition: MCFlux.h:98
double fvx
Definition: MCFlux.h:52
double fnwtfar
Definition: MCFlux.h:48
double fppdxdz
Definition: MCFlux.h:58
int ftgen
Definition: MCFlux.h:83
double fnenergyf
Definition: MCFlux.h:47
void Reset()
Definition: MCFlux.cxx:94
double fnwtnear
Definition: MCFlux.h:44
int fptype
Definition: MCFlux.h:63
double ftvx
Definition: MCFlux.h:76
int fndecay
Definition: MCFlux.h:50
int fnorig
Definition: MCFlux.h:49
double fpdpz
Definition: MCFlux.h:57
double fmuparpz
Definition: MCFlux.h:69
int fevtno
Definition: MCFlux.h:36
double fppdydz
Definition: MCFlux.h:59
double ftpz
Definition: MCFlux.h:81
double ftvy
Definition: MCFlux.h:77
double fpdpy
Definition: MCFlux.h:56
double fvy
Definition: MCFlux.h:53
double fmupare
Definition: MCFlux.h:70
double fppenergy
Definition: MCFlux.h:61
int fntype
Definition: MCFlux.h:51
double fnecm
Definition: MCFlux.h:71
double fmuparpy
Definition: MCFlux.h:68
double fppvx
Definition: MCFlux.h:64
double fpppz
Definition: MCFlux.h:60
double fppvy
Definition: MCFlux.h:65
double fppvz
Definition: MCFlux.h:66
double ftpy
Definition: MCFlux.h:80
double fvz
Definition: MCFlux.h:54
double fmuparpx
Definition: MCFlux.h:67
double fnenergyn
Definition: MCFlux.h:43
double ftvz
Definition: MCFlux.h:78
double fnimpwt
Definition: MCFlux.h:72
void evgb::FillMCTruth ( const genie::EventRecord *  grec,
double  spillTime,
simb::MCTruth mctruth,
const std::string &  genieVersion = "unknown",
const std::string &  genieTune = "unknown",
bool  addGenieVtxTime = false,
const std::unordered_map< std::string, std::string >  genConfig = {} 
)

Definition at line 182 of file GENIE2ART.cxx.

Referenced by evg::AddGenieEventsToArt::produce(), and evgb::GENIEHelper::Sample().

189 {
190  TLorentzVector vtxOffset(0,0,0,spillTime);
191  FillMCTruth(record,vtxOffset,truth,genieVersion,genieTune,addGenieVtxTime, genConfig);
192 }
void FillMCTruth(const genie::EventRecord *grec, double spillTime, simb::MCTruth &mctruth, const std::string &genieVersion="unknown", const std::string &genieTune="unknown", bool addGenieVtxTime=false, const std::unordered_map< std::string, std::string > genConfig={})
Definition: GENIE2ART.cxx:182
void evgb::FillMCTruth ( const genie::EventRecord *  grec,
TLorentzVector &  vtxOffset,
simb::MCTruth mctruth,
const std::string &  genieVersion = "unknown",
const std::string &  genieTune = "unknown",
bool  addGenieVtxTime = false,
const std::unordered_map< std::string, std::string >  genConfig = {} 
)

Definition at line 194 of file GENIE2ART.cxx.

References simb::MCTruth::Add(), simb::kAMNuGamma, simb::kBeamNeutrino, simb::kCC, simb::kCoh, simb::kCohElastic, simb::kDiffractive, simb::kDIS, simb::kElectronScattering, simb::kEM, simb::kGENIE, simb::kGlashowResonance, simb::kIMDAnnihilation, simb::kInverseBetaDecay, simb::kInverseMuDecay, simb::kMEC, simb::kNC, simb::kNuanceOffset, simb::kNuElectronElastic, simb::kQE, simb::kRes, simb::kUnknownInteraction, simb::kWeakMix, part, simb::MCTruth::SetGeneratorInfo(), simb::MCTruth::SetNeutrino(), simb::MCTruth::SetOrigin(), simb::MCParticle::Vx(), x, and y.

201 {
202  // offset vector is assmed to be in (cm,ns) which is MCTruth's units
203 
204  // GENIE's vertex is in (meters,seconds)
205  TLorentzVector *vertex = record->Vertex();
206 
207  // get the Interaction object from the record - this is the object
208  // that talks to the event information objects and is in m
209  const genie::Interaction *inter = record->Summary();
210 
211  // get the different components making up the interaction
212  const genie::InitialState &initState = inter->InitState();
213  const genie::ProcessInfo &procInfo = inter->ProcInfo();
214 
215  //const genie::Kinematics &kine = inter->Kine();
216  //const genie::XclsTag &exclTag = inter->ExclTag();
217  //const genie::KPhaseSpace &phaseSpace = inter->PhaseSpace();
218 
219  // add the particles from the interaction
220  TIter partitr(record);
221  genie::GHepParticle *part = 0;
222  // GHepParticles return units of GeV/c for p.
223  // The V_i are all in fermis and are relative to the center
224  // of the struck nucleus.
225  // prior to GENIE R-3_02_00 time was always zero
226  // thereafter that it is given in units of yoctoseconds (10^{-24})
227  // add the lab vertex X/Y/Z to the V_i for everything
228  // (store the true fermi distance in GVtx to be retrievable)
229 
230  int trackid = 0;
231  std::string primary("primary");
232 
233  /*
234  // for debugging purposes ...
235  mf::LogWarning("GENIE2ART")
236  << "addGenieVtxTime is "
237  << (addGenieVtxTime?"true":"false") << " if true, added "
238  << vertex->T() * 1.0e9 << " ns GENIE Vtx "
239  << genie::utils::print::X4AsString(vertex);
240  */
241 
242  while( (part = dynamic_cast<genie::GHepParticle *>(partitr.Next())) ){
243 
244  simb::MCParticle tpart(trackid,
245  part->Pdg(),
246  primary,
247  part->FirstMother(),
248  part->Mass(),
249  part->Status());
250  double vtx[4] = {part->Vx(), part->Vy(), part->Vz(), part->Vt()};
251 
252  // save the "relative to the nucleus" (fermimeter) particle offsets
253  tpart.SetGvtx(vtx);
254 
255  tpart.SetRescatter(part->RescatterCode());
256 
257  // set the vertex location for the neutrino, nucleus and everything
258  // that is to be tracked. GENIE interaction vertex is in meters.
259  // GENIE individual particles are in fermi and
260  // times are in yoctoseconds (10^{-24})
261  // MCTruth uses units of (cm, ns)
262  // GVtx stores position relative to struck nucleus, so we don't
263  // need to do anything special to recover that info for rewgt purposes
264  vtx[0] = 100.*( part->Vx()*1.e-15 + vertex->X()) + vtxOffset.X();
265  vtx[1] = 100.*( part->Vy()*1.e-15 + vertex->Y()) + vtxOffset.Y();
266  vtx[2] = 100.*( part->Vz()*1.e-15 + vertex->Z()) + vtxOffset.Z();
267  const double yocto2ns = 1.0e-15; // 1.0e-24 sec/yoctosec / 1.0e-9 sec/ns
268  vtx[3] = yocto2ns*part->Vt() + vtxOffset.T();
269  // GENIE vertex time is in seconds, MCTruth time in ns
270  if (addGenieVtxTime) vtx[3] += vertex->T() * 1.0e9;
271 
272  TLorentzVector pos(vtx[0], vtx[1], vtx[2], vtx[3]);
273  TLorentzVector mom(part->Px(), part->Py(), part->Pz(), part->E());
274  tpart.AddTrajectoryPoint(pos,mom);
275  if (part->PolzIsSet()) {
276  TVector3 polz;
277  part->GetPolarization(polz);
278  tpart.SetPolarization(polz);
279  }
280  truth.Add(tpart);
281 
282  ++trackid;
283  }// end loop to convert GHepParticles to MCParticles
284 
285  // is the interaction NC or CC
286  int CCNC = simb::kCC;
287  if (procInfo.IsWeakNC()) CCNC = simb::kNC;
288 
289  // what is the interaction type
290  int mode = simb::kUnknownInteraction;
291 
292  if (procInfo.IsQuasiElastic() ) mode = simb::kQE;
293  else if (procInfo.IsDeepInelastic() ) mode = simb::kDIS;
294  else if (procInfo.IsResonant() ) mode = simb::kRes;
295 #if __GENIE_RELEASE_CODE__ >= GRELCODE(3,2,0)
296  else if (procInfo.IsCoherentProduction() ) mode = simb::kCoh;
297  else if (procInfo.IsCoherentElastic() ) mode = simb::kCohElastic;
298 #else
299  else if (procInfo.IsCoherent() ) mode = simb::kCoh;
300  else if (procInfo.IsCoherentElas() ) mode = simb::kCohElastic;
301 #endif
302  else if (procInfo.IsElectronScattering() ) mode = simb::kElectronScattering;
303  else if (procInfo.IsNuElectronElastic() ) mode = simb::kNuElectronElastic;
304  else if (procInfo.IsInverseMuDecay() ) mode = simb::kInverseMuDecay;
305  else if (procInfo.IsIMDAnnihilation() ) mode = simb::kIMDAnnihilation;
306  else if (procInfo.IsInverseBetaDecay() ) mode = simb::kInverseBetaDecay;
307  else if (procInfo.IsGlashowResonance() ) mode = simb::kGlashowResonance;
308  else if (procInfo.IsAMNuGamma() ) mode = simb::kAMNuGamma;
309  else if (procInfo.IsMEC() ) mode = simb::kMEC;
310  else if (procInfo.IsDiffractive() ) mode = simb::kDiffractive;
311  else if (procInfo.IsEM() ) mode = simb::kEM;
312  else if (procInfo.IsWeakMix() ) mode = simb::kWeakMix;
313 
314  int itype = simb::kNuanceOffset + genie::utils::ghep::NuanceReactionCode(record);
315 
316  // set the neutrino information in MCTruth
317  truth.SetOrigin(simb::kBeamNeutrino);
318  std::unordered_map<std::string, std::string> genConfigCopy(genConfig);
319  genConfigCopy.emplace("tune", genieTune);
320  truth.SetGeneratorInfo(simb::Generator_t::kGENIE, genieVersion, genConfigCopy);
321 
322  // The genie event kinematics are subtle different from the event
323  // kinematics that a experimentalist would calculate
324  // Instead of retriving the genie values for these kinematic variables
325  // calcuate them from the the final state particles
326  // while ingnoring the fermi momentum and the off-shellness of the bound nucleon.
327  genie::GHepParticle * hitnucl = record->HitNucleon();
328  TLorentzVector pdummy(0, 0, 0, 0);
329  // these don't exist if it came from nucleon decay .. RWH
330  const TLorentzVector v4_null;
331  genie::GHepParticle* probe = record->Probe();
332  genie::GHepParticle* finallepton = record->FinalStatePrimaryLepton();
333  const TLorentzVector & k1 = ( probe ? *(probe->P4()) : v4_null );
334  const TLorentzVector & k2 = ( finallepton ? *(finallepton->P4()) : v4_null );
335 
336 #ifdef OLD_KINE_CALC
337  //const TLorentzVector & p1 = (hitnucl) ? *(hitnucl->P4()) : pdummy;
338 
339  double M = genie::constants::kNucleonMass;
340  TLorentzVector q = k1-k2; // q=k1-k2, 4-p transfer
341  double Q2 = -1 * q.M2(); // momemtum transfer
342  double v = (hitnucl) ? q.Energy() : -1; // v (E transfer to the nucleus)
343  double x = (hitnucl) ? 0.5*Q2/(M*v) : -1; // Bjorken x
344  double y = (hitnucl) ? v/k1.Energy() : -1; // Inelasticity, y = q*P1/k1*P1
345  double W2 = (hitnucl) ? M*M + 2*M*v - Q2 : -1; // Hadronic Invariant mass ^ 2
346  double W = (hitnucl) ? std::sqrt(W2) : -1;
347 #else
348  // (same strategy as in gNtpConv.cxx::ConvertToGST().)
349 
350  // also note that since most of these variables are calculated purely
351  // from the leptonic system, they have meaning in reactions that didn't
352  // strike a nucleon (or even a hadron) as well.
353  TLorentzVector q = k1-k2; // q=k1-k2, 4-p transfer
354 
355  double Q2 = -1 * q.M2(); // momemtum transfer
356  double v = q.Energy(); // v (E transfer to the had system)
357  double y = v/k1.Energy(); // Inelasticity, y = q*P1/k1*P1
358  double x, W2, W;
359  x = W2 = W = -1;
360 
361 #if __GENIE_RELEASE_CODE__ >= GRELCODE(3,2,0)
362  if ( hitnucl || procInfo.IsCoherentProduction() ) {
363 #else
364  if ( hitnucl || procInfo.IsCoherent() ) {
365 #endif
366  const double M = genie::constants::kNucleonMass;
367  // Bjorken x.
368  // Rein & Sehgal use this same formulation of x even for Coherent
369  x = 0.5*Q2/(M*v);
370  // Hadronic Invariant mass ^ 2.
371  // ("wrong" for Coherent, but it's "experimental", so ok?)
372  W2 = M*M + 2*M*v - Q2;
373  W = std::sqrt(W2);
374  }
375 #endif
376 
377  truth.SetNeutrino(CCNC,
378  mode,
379  itype,
380  initState.Tgt().Pdg(),
381  initState.Tgt().HitNucPdg(),
382  initState.Tgt().HitQrkPdg(),
383  W,
384  x,
385  y,
386  Q2);
387  return;
388 }
Float_t x
Definition: compare.C:6
neutrino electron elastic scatter
Definition: MCNeutrino.h:140
Float_t y
Definition: compare.C:6
offset to account for adding in Nuance codes to this enum
Definition: MCNeutrino.h:95
TString part[npart]
Definition: Style.C:32
double Vx(const int i=0) const
Definition: MCParticle.h:222
n.b.: this group is similar but not quite, entirely unlike GENIE ScatteringType convention ...
Definition: MCNeutrino.h:80
inverse muon decay
Definition: MCNeutrino.h:141
Beam neutrinos.
Definition: MCTruth.h:23
vertex reconstruction
unsigned int evgb::GetRandomNumberSeed ( )
inline

Definition at line 22 of file evgenbase.h.

Referenced by evg::AddGenieEventsToArt::AddGenieEventsToArt(), and evgb::GENIEHelper::GENIEHelper().

23 {
24 
25  // the maximum allowed seed for the art::RandomNumberGenerator
26  // is 900000000. Use the system random number generator to get a pseudo-random
27  // number for the seed value, and take the modulus of the maximum allowed
28  // seed to ensure we don't ever go over that maximum
29 
30  // Set gRandom to be a TRandom3 based on this state in case we need to pull
31  // random values from histograms, etc
32  TRandom3 *rand = new TRandom3(0);
33  gRandom = rand;
34  return rand->Integer(900000000);
35 }
genie::EventRecord * evgb::RetrieveGHEP ( const simb::MCTruth truth,
const simb::GTruth gtruth,
bool  useFirstTrajPosition = true 
)

return genie::EventRecord pointer; callee takes possession

Definition at line 540 of file GENIE2ART.cxx.

References simb::MCParticle::E(), simb::GTruth::fCharmHadronPdg, simb::GTruth::fDecayMode, simb::GTruth::fDiffXsec, simb::GTruth::fFinalLeptonPdg, simb::GTruth::fFinalQuarkPdg, simb::GTruth::fFShadSystP4, simb::GTruth::fGint, simb::GTruth::fGPhaseSpace, simb::GTruth::fgQ2, simb::GTruth::fgq2, simb::GTruth::fGscatter, simb::GTruth::fgT, simb::GTruth::fgW, simb::GTruth::fgWrun, simb::GTruth::fgX, simb::GTruth::fgY, simb::GTruth::fHitNucPos, simb::GTruth::fIsCharm, simb::GTruth::fIsSeaQuark, simb::GTruth::fIsStrange, simb::GTruth::fNumNeutron, simb::GTruth::fNumPi0, simb::GTruth::fNumPiMinus, simb::GTruth::fNumPiPlus, simb::GTruth::fNumProton, simb::GTruth::fNumRho0, simb::GTruth::fNumRhoMinus, simb::GTruth::fNumRhoPlus, simb::GTruth::fNumSingleGammas, simb::GTruth::fprobability, simb::GTruth::fProbePDG, simb::GTruth::fResNum, simb::GTruth::fStrangeHadronPdg, simb::GTruth::ftgtA, simb::GTruth::ftgtPDG, simb::GTruth::ftgtZ, simb::GTruth::fVertex, simb::GTruth::fweight, simb::GTruth::fXsec, simb::MCTruth::GetNeutrino(), simb::MCTruth::GetParticle(), simb::MCParticle::Gvt(), simb::MCParticle::Gvx(), simb::MCParticle::Gvy(), simb::MCParticle::Gvz(), simb::MCNeutrino::HitNuc(), simb::MCNeutrino::HitQuark(), simb::MCNeutrino::Lepton(), simb::MCParticle::Mother(), simb::MCTruth::NParticles(), simb::MCParticle::NumberTrajectoryPoints(), simb::MCParticle::PdgCode(), simb::MCParticle::Polarization(), simb::MCParticle::Px(), simb::MCParticle::Py(), simb::MCParticle::Pz(), simb::MCParticle::Rescatter(), simb::MCParticle::StatusCode(), and target.

Referenced by evg::GenieOutput::analyze(), rwgt::NuReweight::CalcWeight(), evwgh::GenieWeightCalc::GetWeight(), and RetrieveGHEPs().

543 {
544  genie::EventRecord* newEvent = new genie::EventRecord;
545 
546  newEvent->SetWeight(gtruth.fweight);
547  newEvent->SetProbability(gtruth.fprobability);
548  newEvent->SetXSec(gtruth.fXsec);
549 
550  genie::KinePhaseSpace_t space = (genie::KinePhaseSpace_t)gtruth.fGPhaseSpace;
551 
552  newEvent->SetDiffXSec(gtruth.fDiffXsec,space);
553 
554  TLorentzVector vtx = gtruth.fVertex;
555  newEvent->SetVertex(vtx);
556 
557  //mf::LogWarning("GENIE2ART")
558  // << "####### mctruth.NParticles() " << mctruth.NParticles();
559 
560  for (int i = 0; i < mctruth.NParticles(); i++) {
561  simb::MCParticle mcpart = mctruth.GetParticle(i);
562 
563  int gmid = mcpart.PdgCode();
564  genie::GHepStatus_t gmst = (genie::GHepStatus_t)mcpart.StatusCode();
565  int gmmo = mcpart.Mother();
566  int gmfd = -1;
567  int gmld = -1;
568 
569  /*
570  // GENIE will update daughter references as particles are added
571  // without a need to jump through these hoops ... which gets
572  // it wrong anyway (always sets 0th particles daughter to
573  // mctruth.NParticles()-1, and leave the others at -1 (which then
574  // GENIE handles correctly ...
575 
576  int ndaughters = mcpart.NumberDaughters();
577  //find the track ID of the first and last daughter particles
578  int fdtrkid = 0;
579  int ldtrkid = 0;
580  if (ndaughters !=0) {
581  fdtrkid = mcpart.Daughter(0);
582  if (ndaughters == 1) {
583  ldtrkid = 1;
584  }
585  else if (ndaughters >1) {
586  fdtrkid = mcpart.Daughter(ndaughters-1);
587  }
588  }
589 
590  // Genie uses the index in the particle array to reference
591  // the daughter particles.
592  // MCTruth keeps the particles in the same order so use the
593  // track ID to find the proper index.
594  for (int j = 0; j < mctruth.NParticles(); j++) {
595  simb::MCParticle temp = mctruth.GetParticle(i);
596  if (temp.TrackId() == fdtrkid) {
597  gmfd = j;
598  }
599  if (temp.TrackId() == ldtrkid) {
600  gmld = j;
601  }
602  }
603  */
604 
605  double gmpx = mcpart.Px(0);
606  double gmpy = mcpart.Py(0);
607  double gmpz = mcpart.Pz(0);
608  double gme = mcpart.E(0);
609 
610  double gmvx = mcpart.Gvx();
611  double gmvy = mcpart.Gvy();
612  double gmvz = mcpart.Gvz();
613  double gmvt = mcpart.Gvt();
614 
615  int gmri = mcpart.Rescatter();
616 
617  genie::GHepParticle gpart(gmid, gmst, gmmo, -1, gmfd, gmld,
618  gmpx, gmpy, gmpz, gme, gmvx, gmvy, gmvz, gmvt);
619  gpart.SetRescatterCode(gmri);
620  TVector3 polz = mcpart.Polarization();
621  if (polz.x() !=0 || polz.y() !=0 || polz.z() !=0) {
622  gpart.SetPolarization(polz);
623  }
624  newEvent->AddParticle(gpart);
625  }
626 
627  genie::ProcessInfo proc_info;
628  genie::ScatteringType_t gscty = (genie::ScatteringType_t)gtruth.fGscatter;
629  genie::InteractionType_t ginty = (genie::InteractionType_t)gtruth.fGint;
630 
631  proc_info.Set(gscty,ginty);
632 
633  genie::XclsTag gxt;
634 
635  //Set Exclusive Final State particle numbers
636  genie::Resonance_t gres = (genie::Resonance_t)gtruth.fResNum;
637  gxt.SetResonance(gres);
638  gxt.SetDecayMode(gtruth.fDecayMode);
639  gxt.SetNPions(gtruth.fNumPiPlus, gtruth.fNumPi0, gtruth.fNumPiMinus);
640  gxt.SetNNucleons(gtruth.fNumProton, gtruth.fNumNeutron);
641 #if __GENIE_RELEASE_CODE__ >= GRELCODE(3,2,0)
642  gxt.SetNSingleGammas(gtruth.fNumSingleGammas);
643  gxt.SetNRhos(gtruth.fNumRhoPlus, gtruth.fNumRho0, gtruth.fNumRhoMinus);
644  if ( gtruth.fFinalQuarkPdg != 0 )
645  gxt.SetFinalQuark(gtruth.fFinalQuarkPdg);
646  if ( gtruth.fFinalLeptonPdg != 0 )
647  gxt.SetFinalLepton(gtruth.fFinalLeptonPdg);
648 #endif
649 
650  if (gtruth.fIsCharm) {
651  gxt.SetCharm(gtruth.fCharmHadronPdg);
652  } else {
653  gxt.UnsetCharm();
654  }
655 
656  if (gtruth.fIsStrange) {
657  gxt.SetStrange(gtruth.fStrangeHadronPdg);
658  } else {
659  gxt.UnsetStrange();
660  }
661 
662  // Set the GENIE kinematics info
663  genie::Kinematics gkin;
664  // RWH: really should be looping of GENIE Conventions/KineVar_t enum
665  // and only recording/resetting those that were originally there ...
666  const double flagVal = -99999;
667  if ( gtruth.fgX != flagVal) gkin.Setx(gtruth.fgX, true);
668  if ( gtruth.fgY != flagVal) gkin.Sety(gtruth.fgY, true);
669  if ( gtruth.fgT != flagVal) gkin.Sett(gtruth.fgT, true);
670  if ( gtruth.fgW != flagVal) gkin.SetW(gtruth.fgW, true);
671  if ( gtruth.fgQ2 != flagVal) gkin.SetQ2(gtruth.fgQ2, true);
672  if ( gtruth.fgq2 != flagVal) gkin.Setq2(gtruth.fgq2, true);
673  if ( gtruth.fgWrun != flagVal) gkin.SetW(gtruth.fgWrun, false);
674 
675  simb::MCNeutrino nu = mctruth.GetNeutrino();
676  simb::MCParticle lep = nu.Lepton();
677  // is this even real?
678  if ( lep.NumberTrajectoryPoints() > 0 ) {
679  gkin.SetFSLeptonP4(lep.Px(), lep.Py(), lep.Pz(), lep.E());
680  }
681  gkin.SetHadSystP4(gtruth.fFShadSystP4.Px(),
682  gtruth.fFShadSystP4.Py(),
683  gtruth.fFShadSystP4.Pz(),
684  gtruth.fFShadSystP4.E());
685 
686  // reordering this to avoid warning (A=0,Z=0)
687  int probe_pdgc = gtruth.fProbePDG;
688  int tgtZ = gtruth.ftgtZ;
689  int tgtA = gtruth.ftgtA;
690 
691  //std::cerr << " tgtZ " << tgtZ << " tgtA " << tgtA << " probe " << probe_pdgc << std::endl;
692 
693  // genie::InitialState::Init() will fail if target_pdgc or probe_pdgc
694  // come back with nothign from PDGLibrary::Instance()->Find()
695  // fake it ... (what does nucleon decay do here??)
696  if ( tgtZ == 0 || tgtA == 0 ) { tgtZ = tgtA = 1; } // H1
697  if ( probe_pdgc == 0 || probe_pdgc == -1 ) { probe_pdgc = 22; } // gamma
698 
699  //std::cerr << " tgtZ " << tgtZ << " tgtA " << tgtA << " probe " << probe_pdgc << std::endl;
700 
701  int target_pdgc = genie::pdg::IonPdgCode(tgtA,tgtZ);
702 
703  /*
704  TParticlePDG * t = genie::PDGLibrary::Instance()->Find(target_pdgc);
705  TParticlePDG * p = genie::PDGLibrary::Instance()->Find(probe_pdgc );
706 
707  std::cerr << " target " << target_pdgc << " t " << t << " p " << p << std::endl;
708  */
709 
710  int targetNucleon = nu.HitNuc();
711  int struckQuark = nu.HitQuark();
712 
713  //genie::Target tmptgt(gtruth.ftgtZ, gtruth.ftgtA, targetNucleon);
714  // this ctor doesn't copy the state of the Target beyond the PDG value!
715  // so don't bother creating a tmptgt ...
716  //genie::InitialState ginitstate(tmptgt,probe_pdgc);
717 
718  genie::InitialState ginitstate(target_pdgc,probe_pdgc);
719 
720  // do this here _after_ creating InitialState
721  genie::Target* tgtptr = ginitstate.TgtPtr();
722  tgtptr->SetHitNucPdg(targetNucleon);
723  tgtptr->SetHitNucPosition(gtruth.fHitNucPos);
724  tgtptr->SetHitQrkPdg(struckQuark);
725  tgtptr->SetHitSeaQrk(gtruth.fIsSeaQuark);
726 
727  if (newEvent->HitNucleonPosition() >= 0) {
728  genie::GHepParticle * hitnucleon = newEvent->HitNucleon();
729  std::unique_ptr<TLorentzVector> p4hitnucleon(hitnucleon->GetP4());
730  tgtptr->SetHitNucP4(*p4hitnucleon);
731  } else {
732  if ( targetNucleon != 0 ) {
733  mf::LogWarning("GENIE2ART")
734  << "evgb::RetrieveGHEP() no hit nucleon position "
735  << " but targetNucleon is " << targetNucleon
736  << " at " << __FILE__ << ":" << __LINE__
737  << std::endl << std::flush;
738  }
739  TLorentzVector dummy(0.,0.,0.,0.);
740  tgtptr->SetHitNucP4(dummy);
741  }
742 
743  if (newEvent->TargetNucleusPosition() >= 0) {
744  genie::GHepParticle * target = newEvent->TargetNucleus();
745  std::unique_ptr<TLorentzVector> p4target(target->GetP4());
746  ginitstate.SetTgtP4(*p4target);
747  } else {
748  double Erest = 0.;
749  if ( gtruth.ftgtPDG != 0 ) {
750  TParticlePDG* ptmp = genie::PDGLibrary::Instance()->Find(gtruth.ftgtPDG);
751  if ( ptmp ) Erest = ptmp->Mass();
752  } else {
753  mf::LogWarning("GENIE2ART")
754  << "evgb::RetrieveGHEP() no target nucleus position "
755  << " but gtruth.ftgtPDG is " << gtruth.ftgtPDG
756  << " at " << __FILE__ << ":" << __LINE__
757  << std::endl << std::flush;
758  }
759  TLorentzVector dummy(0.,0.,0.,Erest);
760  ginitstate.SetTgtP4(dummy);
761  }
762 
763  genie::GHepParticle * probe = newEvent->Probe();
764  if ( probe ) {
765  std::unique_ptr<TLorentzVector> p4probe(probe->GetP4());
766  ginitstate.SetProbeP4(*p4probe);
767  } else {
768  // this can happen ...
769  mf::LogDebug("GENIE2ART")
770  << "evgb::RetrieveGHEP() no probe "
771  << " at " << __FILE__ << ":" << __LINE__
772  << std::endl << std::flush;
773  TLorentzVector dummy(0.,0.,0.,0.);
774  ginitstate.SetProbeP4(dummy);
775  }
776 
777  genie::Interaction * p_gint = new genie::Interaction(ginitstate,proc_info);
778 
779  p_gint->SetProcInfo(proc_info);
780  p_gint->SetKine(gkin);
781  p_gint->SetExclTag(gxt);
782  newEvent->AttachSummary(p_gint);
783 
784  /*
785  //For temporary debugging purposes
786  genie::Interaction *inter = newEvent->Summary();
787  const genie::InitialState &initState = inter->InitState();
788  const genie::Target &tgt = initState.Tgt();
789  std::cout << "TargetPDG as Recorded: " << gtruth.ftgtPDG << std::endl;
790  std::cout << "TargetZ as Recorded: " << gtruth.ftgtZ << std::endl;
791  std::cout << "TargetA as Recorded: " << gtruth.ftgtA << std::endl;
792  std::cout << "TargetPDG as Recreated: " << tgt.Pdg() << std::endl;
793  std::cout << "TargetZ as Recreated: " << tgt.Z() << std::endl;
794  std::cout << "TargetA as Recreated: " << tgt.A() << std::endl;
795  */
796 
797  return newEvent;
798 
799 }
double fgW
Definition: GTruth.h:64
double E(const int i=0) const
Definition: MCParticle.h:234
int fGint
interaction code
Definition: GTruth.h:56
unsigned int NumberTrajectoryPoints() const
Definition: MCParticle.h:219
const TVector3 & Polarization() const
Definition: MCParticle.h:215
int PdgCode() const
Definition: MCParticle.h:213
double fgq2
Definition: GTruth.h:63
double fgX
Definition: GTruth.h:66
int ftgtA
Definition: GTruth.h:46
double Py(const int i=0) const
Definition: MCParticle.h:232
double Gvz() const
Definition: MCParticle.h:249
int HitQuark() const
Definition: MCNeutrino.h:153
int Mother() const
Definition: MCParticle.h:214
int fGPhaseSpace
phase space system of DiffXSec
Definition: GTruth.h:32
double Px(const int i=0) const
Definition: MCParticle.h:231
int fNumNeutron
number of neutrons after reaction, before FSI
Definition: GTruth.h:81
int HitNuc() const
Definition: MCNeutrino.h:152
double fgWrun
Definition: GTruth.h:70
int fNumRhoPlus
number of pi pluses after reaction, before FSI
Definition: GTruth.h:87
int ftgtZ
Definition: GTruth.h:45
double fXsec
cross section of interaction
Definition: GTruth.h:30
int fNumPiPlus
number of pi pluses after reaction, before FSI
Definition: GTruth.h:83
int StatusCode() const
Definition: MCParticle.h:212
int fNumSingleGammas
number of gammas after reaction, before FSI
Definition: GTruth.h:85
int fNumPiMinus
number of pi minuses after reaction, before FSI
Definition: GTruth.h:84
double Gvx() const
Definition: MCParticle.h:247
double Gvy() const
Definition: MCParticle.h:248
bool fIsStrange
strange production // added version 13
Definition: GTruth.h:78
int fResNum
resonance number
Definition: GTruth.h:89
int fNumProton
number of protons after reaction, before FSI
Definition: GTruth.h:80
int fDecayMode
Definition: GTruth.h:90
double fprobability
interaction probability
Definition: GTruth.h:29
int fFinalLeptonPdg
Definition: GTruth.h:93
const simb::MCParticle & Lepton() const
Definition: MCNeutrino.h:147
int fProbePDG
Definition: GTruth.h:39
int fGscatter
neutrino scattering code
Definition: GTruth.h:55
int fNumRhoMinus
number of pi minuses after reaction, before FSI
Definition: GTruth.h:88
int fNumPi0
number of pi0 after reaction, before FSI
Definition: GTruth.h:82
int fCharmHadronPdg
Definition: GTruth.h:77
bool fIsCharm
did the interaction produce a charmed hadron?
Definition: GTruth.h:76
double fweight
event interaction weight (genie internal)
Definition: GTruth.h:28
double fHitNucPos
Definition: GTruth.h:52
int ftgtPDG
PDG of Target Nucleus, nucleon only if free.
Definition: GTruth.h:47
double fgQ2
Definition: GTruth.h:62
MaybeLogger_< ELseverityLevel::ELsev_success, false > LogDebug
double Gvt() const
Definition: MCParticle.h:250
int fNumRho0
number of pi0 after reaction, before FSI
Definition: GTruth.h:86
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
TLorentzVector fFShadSystP4
generated final state hadronic system (LAB frame)
Definition: GTruth.h:73
double Pz(const int i=0) const
Definition: MCParticle.h:233
cout<< "-> Edep in the target
Definition: analysis.C:53
int Rescatter() const
Definition: MCParticle.h:253
int fFinalQuarkPdg
Definition: GTruth.h:92
double fgT
Definition: GTruth.h:65
Event generator information.
Definition: MCNeutrino.h:18
bool fIsSeaQuark
Definition: GTruth.h:50
TLorentzVector fVertex
Definition: GTruth.h:26
double fgY
a common running variable to be recorded
Definition: GTruth.h:67
double fDiffXsec
differential cross section of interaction
Definition: GTruth.h:31
int fStrangeHadronPdg
Definition: GTruth.h:79
std::vector< std::unique_ptr< genie::EventRecord > > evgb::RetrieveGHEPs ( const art::Handle< std::vector< simb::MCTruth >> &  mcTruthHandle,
const art::Handle< std::vector< simb::GTruth >> &  gTruthHandle 
)

Definition at line 802 of file GENIE2ART.cxx.

References RetrieveGHEP().

Referenced by RetrieveGHEPs().

804  {
805 
806  std::vector<std::unique_ptr<genie::EventRecord>> gheps;
807 
808  size_t NEventUnits = mcTruthHandle->size();
809  if (mcTruthHandle->size() != gTruthHandle->size()) {
810  NEventUnits = std::min(mcTruthHandle->size(), gTruthHandle->size());
811  }
812  for (size_t eu_it = 0; eu_it < NEventUnits; ++eu_it) {
813  gheps.emplace_back(evgb::RetrieveGHEP(mcTruthHandle->at(eu_it),
814  gTruthHandle->at(eu_it)));
815  }
816 
817  return gheps;
818 
819 }
genie::EventRecord * RetrieveGHEP(const simb::MCTruth &truth, const simb::GTruth &gtruth, bool useFirstTrajPosition=true)
return genie::EventRecord pointer; callee takes possession
Definition: GENIE2ART.cxx:540
std::vector< std::unique_ptr< genie::EventRecord > > evgb::RetrieveGHEPs ( const art::Event e,
std::string  mcTruthLabel,
std::string  gTruthLabel 
)

Definition at line 821 of file GENIE2ART.cxx.

References art::ProductRetriever::getHandle(), and RetrieveGHEPs().

824  {
825 
826  art::Handle<std::vector<simb::MCTruth>> mcTruthHandle = e.getHandle<std::vector<simb::MCTruth>>(mcTruthLabel);
827  art::Handle<std::vector<simb::GTruth>> gTruthHandle = e.getHandle<std::vector<simb::GTruth>>(gTruthLabel);
828 
829  return evgb::RetrieveGHEPs(mcTruthHandle, gTruthHandle);
830 
831 }
Handle< PROD > getHandle(SelectorBase const &) const
std::vector< std::unique_ptr< genie::EventRecord > > RetrieveGHEPs(const art::Handle< std::vector< simb::MCTruth >> &mcTruthHandle, const art::Handle< std::vector< simb::GTruth >> &gTruthHandle)
Definition: GENIE2ART.cxx:802
void evgb::SetEventGeneratorListAndTune ( const std::string &  evtlistname = "",
const std::string &  tunename = "${GENIE_XSEC_TUNE}" 
)

Definition at line 128 of file GENIE2ART.cxx.

References ExpandEnvVar().

Referenced by evwgh::EventWeight::EventWeight(), evgb::GENIEHelper::Initialize(), evgen::NeutronOsc::NeutronOsc(), and evgen::NucleonDecay::NucleonDecay().

130 {
131  // set EventGeneratorList name (if non-blank)
132  // set Tune name (if >= R-3_XX_YY)
133 
134 #ifdef GENIE_PRE_R3
135  mf::LogInfo("GENIE2ART") << "GENIE_PRE_R3 ignore setting tune name: \""
136  << tunename << "\"";
137 #else
138  // Constructor automatically calls grunopt->Init();
139  genie::RunOpt* grunopt = genie::RunOpt::Instance();
140 
141  // SetEventGeneratorList wasn't introduced until R-3
142  std::string expEvtGenListName = evgb::ExpandEnvVar(evtgenlistname);
143  if ( expEvtGenListName != "" ) {
144  grunopt->SetEventGeneratorList(expEvtGenListName);
145  }
146 
147  std::string expTuneName = evgb::ExpandEnvVar(tunename);
148  if ( expTuneName != tunename ) {
149  mf::LogInfo("GENIE2ART") << "TuneName started as '" << tunename << "' "
150  << " converted to " << expTuneName;
151  }
152 
153  // If the XSecSplineList returns a non-empty string as the current tune name,
154  // then genie::RunOpt::BuildTune() has already been called.
155  std::string current_tune = genie::XSecSplineList::Instance()->CurrentTune();
156  if ( current_tune.empty() ) {
157  // We need to build the GENIE tune config
158  mf::LogInfo("GENIE2ART") << "Configuring GENIE tune \""
159  << expTuneName << '\"';
160  grunopt->SetTuneName( expTuneName );
161  grunopt->BuildTune();
162  mf::LogInfo("GENIEHelper") << *(grunopt->Tune());
163  }
164  else {
165  // It has already been built, so just check consistency
166  if ( expTuneName != current_tune) {
167  throw cet::exception("TuneNameMismatch") << "Requested GENIE tune \""
168  << expTuneName << "\" does not match previously built tune \""
169  << current_tune << '\"';
170  }
171  }
172 #endif
173 
174 }
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
std::string ExpandEnvVar(const std::string &s)
Definition: GENIE2ART.cxx:98
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33

Variable Documentation

const int evgb::kNue = 0
static

Definition at line 206 of file GENIEHelper.cxx.

Referenced by evgb::GENIEHelper::Sample().

const int evgb::kNueBar = 1
static

Definition at line 207 of file GENIEHelper.cxx.

Referenced by evgb::GENIEHelper::Sample().

const int evgb::kNuMu = 2
static

Definition at line 208 of file GENIEHelper.cxx.

Referenced by evgb::GENIEHelper::Sample().

const int evgb::kNuMuBar = 3
static

Definition at line 209 of file GENIEHelper.cxx.

Referenced by evgb::GENIEHelper::Sample().

const int evgb::kNuTau = 4
static

Definition at line 210 of file GENIEHelper.cxx.

Referenced by evgb::GENIEHelper::Sample().

const int evgb::kNuTauBar = 5
static

Definition at line 211 of file GENIEHelper.cxx.

Referenced by evgb::GENIEHelper::Sample().

template<class T >
double(T::* evgb::RNGWrapper) (void)

Definition at line 82 of file CRYHelper.h.