LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
GENIE2ART.cxx
Go to the documentation of this file.
1 
11 #include "GENIE2ART.h"
12 //#include <math.h>
13 //#include <map>
14 #include <fstream>
15 #include <memory> // for unique_ptr
16 #include <typeinfo>
17 
18 // ROOT includes
19 #include "TVector3.h"
20 #include "TLorentzVector.h"
21 #include "TSystem.h"
22 
23 //GENIE includes
24 #ifdef GENIE_PRE_R3
25  #include "Conventions/GVersion.h"
26  #include "Conventions/Units.h"
27  #include "EVGCore/EventRecord.h"
28  #include "GHEP/GHepUtils.h"
29  #include "PDG/PDGCodes.h"
30  #include "PDG/PDGLibrary.h"
31  #include "GENIE/Utils/RunOpt.h"
32 
33  #include "Interaction/InitialState.h"
34  #include "Interaction/Interaction.h"
35  #include "Interaction/Kinematics.h"
36  #include "Interaction/KPhaseSpace.h"
37  #include "Interaction/ProcessInfo.h"
38  #include "Interaction/XclsTag.h"
39  #include "GHEP/GHepParticle.h"
40  #include "PDG/PDGCodeList.h"
41  #include "Conventions/Constants.h" //for calculating event kinematics
42 
43  #include "EVGDrivers/GFluxI.h"
44  #include "FluxDrivers/GFluxBlender.h"
45  #include "FluxDrivers/GNuMIFlux.h"
46  #include "FluxDrivers/GSimpleNtpFlux.h"
47 #else
48  // GENIE R-3 reorganized headers
49  #include "GENIE/Framework/Conventions/GVersion.h"
50  #include "GENIE/Framework/Conventions/Units.h"
51  #include "GENIE/Framework/Conventions/Constants.h" //for calculating event kinematics
52  #include "GENIE/Framework/Utils/PrintUtils.h"
53  #include "GENIE/Framework/ParticleData/PDGCodes.h"
54  #include "GENIE/Framework/ParticleData/PDGCodeList.h"
55  #include "GENIE/Framework/ParticleData/PDGLibrary.h"
56  #include "GENIE/Framework/GHEP/GHepUtils.h"
57  #include "GENIE/Framework/GHEP/GHepParticle.h"
58  #include "GENIE/Framework/Utils/RunOpt.h"
59  #include "GENIE/Framework/Utils/XSecSplineList.h"
60 
61  #include "GENIE/Framework/Interaction/InitialState.h"
62  #include "GENIE/Framework/Interaction/Interaction.h"
63  #include "GENIE/Framework/Interaction/Kinematics.h"
64  #include "GENIE/Framework/Interaction/KPhaseSpace.h"
65  #include "GENIE/Framework/Interaction/ProcessInfo.h"
66  #include "GENIE/Framework/Interaction/XclsTag.h"
67 
68  #include "GENIE/Framework/EventGen/EventRecord.h"
69  #include "GENIE/Framework/EventGen/GFluxI.h"
70  #include "GENIE/Tools/Flux/GFluxBlender.h"
71  #include "GENIE/Tools/Flux/GNuMIFlux.h"
72  #include "GENIE/Tools/Flux/GSimpleNtpFlux.h"
73 #endif
74 
75 #ifndef ART_V1
81 #else
82  #include "SimulationBase/MCTruth.h"
85  #include "SimulationBase/GTruth.h"
86  #include "SimulationBase/MCFlux.h"
87 #endif
88 
89 // dk2nu
90 #include "dk2nu/tree/dk2nu.h"
91 #include "dk2nu/tree/NuChoice.h"
92 #include "dk2nu/tree/dkmeta.h"
93 #include "dk2nu/genie/GDk2NuFlux.h"
94 
96 #include "cetlib_except/exception.h"
97 
98 std::string evgb::ExpandEnvVar(const std::string& s)
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 }
127 
128 void evgb::SetEventGeneratorListAndTune(const std::string& evtgenlistname,
129  const std::string& tunename)
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 }
175 
176 
177 //---------------------------------------------------------------------------
178 //choose a spill time (ns) to shift the vertex times by:
179 // double spillTime = fGlobalTimeOffset +
180 // fHelperRandom->Uniform()*fRandomTimeOffset;
181 
182 void evgb::FillMCTruth(const genie::EventRecord *record,
183  double spillTime,
184  simb::MCTruth &truth,
185  const std::string & genieVersion,
186  const std::string & genieTune,
187  bool addGenieVtxTime,
188  const std::unordered_map<std::string, std::string> genConfig)
189 {
190  TLorentzVector vtxOffset(0,0,0,spillTime);
191  FillMCTruth(record,vtxOffset,truth,genieVersion,genieTune,addGenieVtxTime, genConfig);
192 }
193 
194 void evgb::FillMCTruth(const genie::EventRecord *record,
195  TLorentzVector &vtxOffset,
196  simb::MCTruth &truth,
197  const std::string & genieVersion,
198  const std::string & genieTune,
199  bool addGenieVtxTime,
200  const std::unordered_map<std::string, std::string> genConfig)
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
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 }
389 
390 //---------------------------------------------------------------------------
391 void evgb::FillGTruth(const genie::EventRecord* record,
392  simb::GTruth& truth) {
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 }
538 
539 //---------------------------------------------------------------------------
540 genie::EventRecord* evgb::RetrieveGHEP(const simb::MCTruth& mctruth,
541  const simb::GTruth& gtruth,
542  bool useFirstTrajPosition)
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 }
800 
801 //---------------------------------------------------------------------------
802 std::vector<std::unique_ptr<genie::EventRecord>> evgb::RetrieveGHEPs(
803  const art::Handle<std::vector<simb::MCTruth>>& mcTruthHandle,
804  const art::Handle<std::vector<simb::GTruth>>& gTruthHandle){
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 }
820 
821 std::vector<std::unique_ptr<genie::EventRecord>> evgb::RetrieveGHEPs(
822  const art::Event& e,
823  std::string mcTruthLabel,
824  std::string gTruthLabel){
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 }
832 
833 //---------------------------------------------------------------------------
834 void evgb::FillMCFlux(genie::GFluxI* fdriver, simb::MCFlux& mcflux)
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 }
880 
881 //---------------------------------------------------------------------------
882 
883 void evgb::FillMCFlux(genie::flux::GNuMIFlux* gnumi,
884  simb::MCFlux& flux)
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 }
892 void evgb::FillMCFlux(const genie::flux::GNuMIFluxPassThroughInfo* nflux,
893  double dk2gen,
894  simb::MCFlux& flux)
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 }
978 
979 //---------------------------------------------------------------------------
980 void evgb::FillMCFlux(genie::flux::GSimpleNtpFlux* gsimple,
981  simb::MCFlux& flux)
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 }
993 void evgb::FillMCFlux(const genie::flux::GSimpleNtpEntry* nflux_entry,
994  const genie::flux::GSimpleNtpNuMI* nflux_numi,
995  const genie::flux::GSimpleNtpAux* nflux_aux,
996  const genie::flux::GSimpleNtpMeta* nflux_meta,
997  simb::MCFlux& flux)
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 }
1150 
1151 //---------------------------------------------------------------------------
1152 void evgb::FillMCFlux(genie::flux::GDk2NuFlux* gdk2nu,
1153  simb::MCFlux& flux)
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 }
1161 void evgb::FillMCFlux(const bsim::Dk2Nu* dk2nu,
1162  const bsim::NuChoice* nuchoice,
1163  simb::MCFlux& flux)
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 }
1317 
1318 //---------------------------------------------------------------------------
1319 
1320 
double fgW
Definition: GTruth.h:64
Float_t x
Definition: compare.C:6
double E(const int i=0) const
Definition: MCParticle.h:234
int fGint
interaction code
Definition: GTruth.h:56
Unified ntuple flux format (replaces 2)
Definition: MCFlux.h:23
int frun
Definition: MCFlux.h:35
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 fnenergy
Definition: MCFlux.h:40
const simb::MCNeutrino & GetNeutrino() const
Definition: MCTruth.h:77
double ftgppy
Definition: MCFlux.h:86
double fgX
Definition: GTruth.h:66
int ftgtA
Definition: GTruth.h:46
double fpdpx
Definition: MCFlux.h:55
double Py(const int i=0) const
Definition: MCParticle.h:232
double Gvz() const
Definition: MCParticle.h:249
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
double ftpx
Definition: MCFlux.h:79
int HitQuark() const
Definition: MCNeutrino.h:153
double ftprivy
Definition: MCFlux.h:89
Full flux simulation ntuple.
Definition: MCFlux.h:21
void SetOrigin(simb::Origin_t origin)
Definition: MCTruth.h:82
neutrino electron elastic scatter
Definition: MCNeutrino.h:140
int fppmedium
Definition: MCFlux.h:62
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
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
double fbeamx
Definition: MCFlux.h:91
void FillMCFlux(genie::GFluxI *fdriver, simb::MCFlux &mcflux)
Definition: GENIE2ART.cxx:834
int Mother() const
Definition: MCParticle.h:214
double fbeampy
Definition: MCFlux.h:95
int ftptype
Definition: MCFlux.h:82
Float_t y
Definition: compare.C:6
int fGPhaseSpace
phase space system of DiffXSec
Definition: GTruth.h:32
simb::flux_code_ fFluxType
Definition: MCFlux.h:98
double Px(const int i=0) const
Definition: MCParticle.h:231
int fNumNeutron
number of neutrons after reaction, before FSI
Definition: GTruth.h:81
double fvx
Definition: MCFlux.h:52
int HitNuc() const
Definition: MCNeutrino.h:152
double fnwtfar
Definition: MCFlux.h:48
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
Functions for transforming GENIE objects into ART objects (and back)
int fNumPiPlus
number of pi pluses after reaction, before FSI
Definition: GTruth.h:83
offset to account for adding in Nuance codes to this enum
Definition: MCNeutrino.h:95
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 fppdxdz
Definition: MCFlux.h:58
double Gvx() const
Definition: MCParticle.h:247
int NParticles() const
Definition: MCTruth.h:75
void SetNeutrino(int CCNC, int mode, int interactionType, int target, int nucleon, int quark, double w, double x, double y, double qsqr)
Definition: MCTruth.cxx:30
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 Gvy() const
Definition: MCParticle.h:248
Particle class.
double ftprivx
Definition: MCFlux.h:88
bool fIsStrange
strange production // added version 13
Definition: GTruth.h:78
object containing MC flux information
TLorentzVector fProbeP4
Definition: GTruth.h:41
double fndydzfar
Definition: MCFlux.h:46
int ftgptype
Definition: MCFlux.h:84
int fResNum
resonance number
Definition: GTruth.h:89
int fNumProton
number of protons after reaction, before FSI
Definition: GTruth.h:80
TString part[npart]
Definition: Style.C:32
int fDecayMode
Definition: GTruth.h:90
double fprobability
interaction probability
Definition: GTruth.h:29
double fnenergyf
Definition: MCFlux.h:47
void Reset()
Definition: MCFlux.cxx:94
double ftprivz
Definition: MCFlux.h:90
int fFinalLeptonPdg
Definition: GTruth.h:93
double fbeamz
Definition: MCFlux.h:93
double fnwtnear
Definition: MCFlux.h:44
const simb::MCParticle & Lepton() const
Definition: MCNeutrino.h:147
int fptype
Definition: MCFlux.h:63
int fProbePDG
Definition: GTruth.h:39
double ftvx
Definition: MCFlux.h:76
int fGscatter
neutrino scattering code
Definition: GTruth.h:55
int fndecay
Definition: MCFlux.h:50
double fndydz
Definition: MCFlux.h:38
int fNumRhoMinus
number of pi minuses after reaction, before FSI
Definition: GTruth.h:88
std::string ExpandEnvVar(const std::string &s)
Definition: GENIE2ART.cxx:98
int fnorig
Definition: MCFlux.h:49
double fbeamy
Definition: MCFlux.h:92
double ftgppz
Definition: MCFlux.h:87
double fpdpz
Definition: MCFlux.h:57
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 fmuparpz
Definition: MCFlux.h:69
double fweight
event interaction weight (genie internal)
Definition: GTruth.h:28
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
TLorentzVector fHitNucP4
Definition: GTruth.h:51
double fbeampz
Definition: MCFlux.h:96
double fndxdzfar
Definition: MCFlux.h:45
double fpdpy
Definition: MCFlux.h:56
double fnpz
Definition: MCFlux.h:39
const simb::MCParticle & GetParticle(int i) const
Definition: MCTruth.h:76
double fzpoint
Definition: MCFlux.h:75
double fvy
Definition: MCFlux.h:53
double fmupare
Definition: MCFlux.h:70
double fxpoint
Definition: MCFlux.h:73
void SetGeneratorInfo(simb::Generator_t generator, const std::string &genVersion, const std::unordered_map< std::string, std::string > &genConfig)
Definition: MCTruth.h:84
double fndxdznea
Definition: MCFlux.h:41
double Vx(const int i=0) const
Definition: MCParticle.h:222
void Add(simb::MCParticle const &part)
Definition: MCTruth.h:80
n.b.: this group is similar but not quite, entirely unlike GENIE ScatteringType convention ...
Definition: MCNeutrino.h:80
Handle< PROD > getHandle(SelectorBase const &) const
void SetEventGeneratorListAndTune(const std::string &evtlistname="", const std::string &tunename="${GENIE_XSEC_TUNE}")
Definition: GENIE2ART.cxx:128
double fppenergy
Definition: MCFlux.h:61
int fntype
Definition: MCFlux.h:51
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 fnecm
Definition: MCFlux.h:71
inverse muon decay
Definition: MCNeutrino.h:141
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
TLorentzVector fTgtP4
Definition: GTruth.h:42
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
int fFinalQuarkPdg
Definition: GTruth.h:92
double fbeampx
Definition: MCFlux.h:94
double fppvz
Definition: MCFlux.h:66
double fgT
Definition: GTruth.h:65
double ftpy
Definition: MCFlux.h:80
void FillGTruth(const genie::EventRecord *grec, simb::GTruth &gtruth)
Definition: GENIE2ART.cxx:391
double fndxdz
Definition: MCFlux.h:37
Event generator information.
Definition: MCTruth.h:32
double fvz
Definition: MCFlux.h:54
Float_t e
Definition: plot.C:35
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
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
double fmuparpx
Definition: MCFlux.h:67
A simplified flux ntuple for quick running.
Definition: MCFlux.h:22
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
Beam neutrinos.
Definition: MCTruth.h:23
double fnenergyn
Definition: MCFlux.h:43
double ftvz
Definition: MCFlux.h:78
double fnimpwt
Definition: MCFlux.h:72
int fStrangeHadronPdg
Definition: GTruth.h:79
vertex reconstruction