GiBUURegen_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // \author Christopher Backhouse -- bckhouse@caltech.edu
3 ////////////////////////////////////////////////////////////////////////
4 
6 
7 #include <cassert>
8 #include <string>
9 #include <vector>
10 
11 #include <wordexp.h>
12 
13 // Framework includes
19 #include "fhiclcpp/ParameterSet.h"
23 
24 // NOvA includes
28 
29 #include "SummaryData/POTSum.h"
30 #include "SummaryData/SpillData.h"
31 #include "Utilities/AssociationUtil.h"
32 
33 #include "dk2nu/tree/dk2nu.h"
34 #include "dk2nu/tree/NuChoice.h"
35 
36 #include "TFile.h"
37 #include "TGraph.h"
38 #include "TH1.h"
39 #include "TTree.h"
40 #include "TRandom3.h"
41 
42 namespace
43 {
44  // That we have GiBUU record files for
45  const std::map<int, std::string> nucleus_to_label = {{1000010010, "H1" },
46  {1000060120, "C12" },
47  {1000080160, "O16" },
48  {1000170350, "Cl35"},
49  {1000220480, "Ti48"},
50  {1000260560, "Fe56"}};
51 
52  std::map<int, std::string> nu_to_genie_label = {{-12, "nu_e_bar"},
53  {+12, "nu_e"},
54  {-14, "nu_mu_bar"},
55  {+14, "nu_mu"}};
56 
57  std::map<int, std::string> nu_to_gibuu_label = {{12, "nue"},
58  {14, "numu"}};
59 }
60 
61 namespace gibuu
62 {
63  namespace
64  {
65  struct GiBUURegenParams
66  {
67  template<class T> using Atom = fhicl::Atom<T>;
68  template<class T> using Sequence = fhicl::Sequence<T>;
69  template<class T> using Table = fhicl::Table<T>;
70  using Comment = fhicl::Comment;
71  using Name = fhicl::Name;
72 
73  Atom<std::string> generatorLabel
74  {
75  Name("GeneratorLabel"),
76  Comment("Where to find the input GENIE MCTruths.")
77  };
78 
79  Atom<std::string> libraryPath
80  {
81  Name("LibraryPath"),
82  Comment("Where to find the GiBUU record files.")
83  };
84 
85  Atom<bool> lowMemory
86  {
87  Name("LowMemory"),
88  Comment("Use less memory at the cost of more disk access.")
89  };
90 
91  Sequence<double> boundingBoxLo
92  {
93  Name("BoundingBoxLo"),
94  Comment("Discard vertices outside this box. Empty = keep all.")
95  };
96 
97  Sequence<double> boundingBoxHi
98  {
99  Name("BoundingBoxHi"),
100  Comment("Discard vertices outside this box. Empty = keep all.")
101  };
102 
103  Atom<bool> copyMCFlux
104  {
105  Name("CopyMCFlux"),
106  Comment("Whether to duplicate generatorLabel's MCFlux objects")
107  };
108 
109  Atom<bool> copyGTruth
110  {
111  Name("CopyGTruth"),
112  Comment("Whether to duplicate generatorLabel's GTruth objects")
113  };
114 
115  Atom<bool> copyBSim
116  {
117  Name("CopyBSim"),
118  Comment("Whether to duplicate generatorLabel's bsim::* objects")
119  };
120  };
121  }
122 
124  {
125  public:
126  // Allows 'nova --print-description' to work
128 
129  explicit GiBUURegen(const Parameters& params);
130  virtual ~GiBUURegen();
131 
132  void beginJob();
133  void produce(art::Event& evt);
134  void endSubRun(art::SubRun& sr);
135 
136  protected:
137  std::string ExpandLibraryPath() const;
138 
139  void LoadGiBUURecords();
140  void LoadGenieXSecs();
141  void LoadGiBUUCorrs();
142 
143  simb::MCTruth GetEvent(const simb::MCTruth& genieTruth,
144  int& trackId) const;
145 
146  struct Key
147  {
148  Key(int _nucl_pdg, int _nu_pdg, bool _iscc) : nucl_pdg(_nucl_pdg), nu_pdg(_nu_pdg), iscc(_iscc) {}
149  bool operator<(const Key& k) const
150  {
151  return (std::make_tuple( nucl_pdg, nu_pdg, iscc) <
152  std::make_tuple(k.nucl_pdg, k.nu_pdg, k.iscc));
153  }
154 
155  int nucl_pdg;
156  int nu_pdg;
157  bool iscc;
158  };
159 
160  Key GetKey(const simb::MCTruth& genieTruth) const;
161  const IRecordList* GetRecordList(const simb::MCTruth& genieTruth) const;
162 
163  std::vector<TVector3> Basis(TVector3 z) const;
164 
165  void Kinematics(const TLorentzVector& nu,
166  const TLorentzVector& lep,
167  double& W, double& x, double& y, double& Qsq) const;
168 
169  double GetGenieXSec(const simb::MCNeutrino& mcn) const;
170 
171  simb::MCTruth CopyGenieEvent(const simb::MCTruth& genieTruth,
172  int& trackId) const;
173 
174  GiBUURegenParams fParams;
178 
179  mutable int fNGiBUU, fNGENIE;
180 
181  std::map<Key, const IRecordList*> fRecords;
182  std::map<Key, TGraph> fGenieXSec;
183  std::map<Key, TH1*> fGiBUUCorr;
184  };
185 
186  //___________________________________________________________________________
188  : fParams(params()),
189  fSpillToken(consumes<sumdata::SpillData>(fParams.generatorLabel())),
190  fTruthToken(consumes<std::vector<simb::MCTruth>>(fParams.generatorLabel())),
191  fPOTSumToken(consumes<sumdata::POTSum, art::InSubRun>(fParams.generatorLabel())),
192  fNGiBUU(0), fNGENIE(0)
193  {
194  produces<std::vector<simb::MCTruth>>();
195  // Associate every truth with the flux it came from
196  produces<art::Assns<simb::MCTruth, simb::MCFlux>>();
197 
198  if(fParams.copyMCFlux()){
199  produces<std::vector<simb::MCFlux>>();
200  }
201  if(fParams.copyGTruth()){
202  produces<std::vector<simb::GTruth>>();
203  produces<art::Assns<simb::MCTruth, simb::GTruth>>();
204  }
205  if(fParams.copyBSim()){
206  produces<std::vector<bsim::Dk2Nu>>();
207  produces<std::vector<bsim::NuChoice>>();
208 #ifdef PUT_DK2NU_ASSN
209  produces<art::Assns<simb::MCTruth, bsim::Dk2Nu>>();
210  produces<art::Assns<simb::MCTruth, bsim::NuChoice>>();
211 #endif
212  }
213 
214  produces<sumdata::SpillData>();
215  produces<sumdata::POTSum, art::InSubRun>();
216 
217  assert(fParams.boundingBoxLo().size() == fParams.boundingBoxHi().size() &&
218  (fParams.boundingBoxLo.size() == 0 ||
219  fParams.boundingBoxLo.size() == 3));
220 
221  const bool onGrid = (getenv("_CONDOR_SCRATCH_DIR") != 0);
222  const std::string libPath = ExpandLibraryPath();
223 
224  // On the grid, absolute path, and it's not to CVMFS = likely bluearc
225  if(onGrid && libPath.find("/") == 0 && libPath.find("/cvmfs/") != 0){
226  std::cout << "Refusing to run on grid with GiBUURegen's LibraryPath set to '"
227  << libPath
228  << "' which is presumably bluearc. Use CVMFS version or arrange your job to copy the files locally." << std::endl;
229  abort();
230  }
231  }
232 
233  //___________________________________________________________________________
235  {
236  if(fNGENIE > 0){
237  std::cout << "GiBUURegen: " << (100*fNGENIE) / (fNGiBUU+fNGENIE)
238  << "% of generated events copied from GENIE" << std::endl;
239  }
240 
241  for(auto it: fRecords) delete it.second;
242  }
243 
244  //___________________________________________________________________________
246  {
247  wordexp_t p;
248  const int status = wordexp(fParams.libraryPath().c_str(), &p, WRDE_SHOWERR | WRDE_UNDEF);
249  if(status != 0){
250  std::cerr << "LibraryPath string '" << fParams.libraryPath()
251  << "' returned error " << status << " from wordexp()."
252  << std::endl;
253  abort();
254  }
255 
256  if(p.we_wordc == 0){
257  std::cerr << "LibraryPath string '" << fParams.libraryPath()
258  << "' didn't expand to anything."
259  << std::endl;
260  abort();
261  }
262 
263  if(p.we_wordc > 1){
264  std::cerr << "LibraryPath string '" << fParams.libraryPath()
265  << "' expanded to " << p.we_wordc << " locations."
266  << std::endl;
267  abort();
268  }
269 
270  const std::string ret = p.we_wordv[0];
271 
272  wordfree(&p);
273 
274  return ret;
275  }
276 
277  //___________________________________________________________________________
279  {
280  LoadGenieXSecs();
281  LoadGiBUUCorrs();
283  }
284 
285  //___________________________________________________________________________
287  {
289 
290  for(bool iscc: {true, false}){
291  for(int sign: {+1, -1}){
292  for(int pdg: {12, 14}){
293  if(!iscc && pdg == 12) continue;
294  for(std::pair<int, std::string> it: nucleus_to_label){
295 
296  const TString gibuuStr =
297  TString::Format("%s/%s%s_%s/%s/records.root",
298  dir.c_str(),
299  iscc ? "cc" : "nc",
300  (sign < 0) ? "bar" : "",
301  nu_to_gibuu_label[pdg].c_str(),
302  it.second.c_str());
303 
304  const Key key(it.first, sign*pdg, iscc);
305 
306  if(fParams.lowMemory())
307  fRecords[key] = new OnDemandRecordList(gibuuStr.Data());
308  else
309  fRecords[key] = new SimpleRecordList(gibuuStr.Data());
310  } // end for nucleus
311  } // end for pdg
312  } // end for sign
313  } // end for iscc
314  }
315 
316  //___________________________________________________________________________
318  {
319  const char* xsec_dir = getenv("GENIEXSECPATH");
320  if(!xsec_dir){
321  std::cout << "GiBUURegen: Environment variable $GENIEXSECPATH is not set." << std::endl;
322  abort();
323  }
324  const std::string xsec_fname = xsec_dir+std::string("/xsec_graphs.root");
325  std::cout << "GiBUURegen: Loading GENIE cross-sections from " << xsec_fname << std::endl;
326 
327  TFile f_xsec(xsec_fname.c_str());
328  assert(!f_xsec.IsZombie());
329 
330  for(bool iscc: {false, true}){
331  for(int sign: {+1, -1}){
332  for(int pdg: {12, 14}){
333  for(std::pair<int, std::string> it: nucleus_to_label){
334  const std::string genieStr = nu_to_genie_label[sign*pdg]+"_"+it.second+"/tot_"+(iscc ? "cc" : "nc");
335  std::cout << " " << genieStr << std::endl;
336  TGraph* g = (TGraph*)f_xsec.Get(genieStr.c_str());
337  assert(g);
338 
339  fGenieXSec[Key(it.first, sign*pdg, iscc)] = TGraph(*g);
340  // ROOT6 only?
341  // fGenieXSec[Key(it.first, sign*pdg, iscc)].SetBit(TGraph::kIsSortedX);
342  } // end for nucleus
343  } // end for pdg
344  } // end for sign
345  } // end for iscc
346  }
347 
348  //___________________________________________________________________________
350  {
352 
353  for(bool iscc: {true, false}){
354  for(int sign: {+1, -1}){
355  for(int pdg: {12, 14}){
356  if(!iscc && pdg == 12) continue;
357  for(std::pair<int, std::string> it: nucleus_to_label){
358 
359  const TString gibuuStr =
360  TString::Format("%s/%s%s_%s/%s/corr.root",
361  dir.c_str(),
362  iscc ? "cc" : "nc",
363  (sign < 0) ? "bar" : "",
364  nu_to_gibuu_label[pdg].c_str(),
365  it.second.c_str());
366 
367  const Key key(it.first, sign*pdg, iscc);
368 
369  TFile f_corr(gibuuStr.Data());
370  assert(!f_corr.IsZombie());
371  fGiBUUCorr[key] = (TH1*)f_corr.Get("corr");
372  assert(fGiBUUCorr[key]);
373  fGiBUUCorr[key]->SetDirectory(0);
374  }
375  }
376  }
377  }
378  }
379 
380  //___________________________________________________________________________
382  {
383  // Always filled
384  auto truthcol = std::make_unique<std::vector<simb::MCTruth>>();
385  auto assns = std::make_unique<art::Assns<simb::MCTruth, simb::MCFlux>>();
386  // Filled if CopyMCFlux
387  auto fluxcol = std::make_unique<std::vector<simb::MCFlux>>();
388  // Filled if CopyGTruth
389  auto gtruthcol = std::make_unique<std::vector<simb::GTruth>>();
390  auto assnsGTruth = std::make_unique<art::Assns<simb::MCTruth, simb::GTruth>>();
391  // Filled if CopyBSim
392  auto dk2nucol = std::make_unique<std::vector<bsim::Dk2Nu>>();
393  auto nuchoicecol = std::make_unique<std::vector<bsim::NuChoice>>();
394 #ifdef PUT_DK2NU_ASSN
395  auto assnsdk2nu = std::make_unique<art::Assns<simb::MCTruth, bsim::Dk2Nu>>();
396  auto assnsnuchoice = std::make_unique<art::Assns<simb::MCTruth, bsim::NuChoice>>();
397 #endif
398 
399  // First up, copy GENIE's POT accounting
401  evt.getByToken(fSpillToken, spilldata);
402  evt.put(std::make_unique<sumdata::SpillData>(*spilldata));
403 
404 
406  evt.getByToken(fTruthToken, genieTruths);
407 
408  art::FindManyP<simb::MCFlux> fmp(genieTruths, evt, fParams.generatorLabel());
409 
410  art::FindManyP<simb::GTruth> fmt(genieTruths, evt, fParams.generatorLabel());
411 
412 #ifdef PUT_DK2NU_ASSN
413  art::FindManyP<bsim::Dk2Nu> fmdk(genieTruths, evt, fParams.generatorLabel());
414  art::FindManyP<bsim::NuChoice> fmnc(genieTruths, evt, fParams.generatorLabel());
415 #else
416  // Otherwise assume they match by index
418  evt.getByLabel(fParams.generatorLabel(), dk2nus);
420  evt.getByLabel(fParams.generatorLabel(), nuchoices);
421 #endif
422 
423  int trackId = 0;
424 
425  for(unsigned int i = 0; i < genieTruths->size(); ++i){
426  const simb::MCTruth& genieTruth = (*genieTruths)[i];
427 
428  const TVector3 vtx = genieTruth.GetNeutrino().Nu().Position().Vect();
429  const std::vector<double> r0 = fParams.boundingBoxLo();
430  const std::vector<double> r1 = fParams.boundingBoxHi();
431 
432  if(!r0.empty() &&
433  (vtx.X() < r0[0] || vtx.Y() < r0[1] || vtx.Z() < r0[2] ||
434  vtx.X() > r1[0] || vtx.Y() > r1[1] || vtx.Z() > r1[2])) continue;
435 
436  truthcol->push_back(GetEvent(genieTruth, trackId));
437 
438  std::vector<art::Ptr<simb::MCFlux>> fluxes = fmp.at(i);
439  assert(fluxes.size() == 1);
440 
441  // Associate most recently added truth with the flux
442  if(fParams.copyMCFlux()){
443  fluxcol->push_back(*fluxes[0]);
444  util::CreateAssn(*this, evt, *truthcol, *fluxcol, *assns, i, i+1);
445  }
446  else{
447  util::CreateAssn(*this, evt, *truthcol, fluxes[0], *assns);
448  }
449 
450  if(fParams.copyGTruth()){
451  std::vector<art::Ptr<simb::GTruth>> gtruths = fmt.at(i);
452  assert(gtruths.size() == 1);
453  gtruthcol->push_back(*gtruths[0]);
454  util::CreateAssn(*this, evt, *truthcol, *gtruthcol, *assnsGTruth, i, i+1);
455  }
456 
457  if(fParams.copyBSim()){
458 #ifdef PUT_DK2NU_ASSN
459  std::vector<art::Ptr<bsim::Dk2Nu>> dk2nus = fmdk.at(i);
460  assert(dk2nus.size() == 1);
461  dk2nucol->push_back(*dk2nus[0]);
462 
463  std::vector<art::Ptr<bsim::NuChoice>> nuchoices = fmnc.at(i);
464  assert(nuchoices.size() == 1);
465  nuchoicecol->push_back(*nuchoices[0]);
466 
467  util::CreateAssn(*this, evt, *truthcol, *dk2nucol, *assnsdk2nu, i, i+1);
468  util::CreateAssn(*this, evt, *truthcol, *nuchoicecol, *assnsnuchoice, i, i+1);
469 #else
470  if(dk2nus->size() > i) dk2nucol->push_back((*dk2nus)[i]);
471  if(nuchoices->size() > i) nuchoicecol->push_back((*nuchoices)[i]);
472 #endif
473  }
474 
475  } // end for i
476 
477  // put the collections in the event
478  evt.put(std::move(truthcol));
479  evt.put(std::move(assns));
480  if(fParams.copyMCFlux()){
481  evt.put(std::move(fluxcol));
482  }
483  if(fParams.copyGTruth()){
484  evt.put(std::move(gtruthcol));
485  evt.put(std::move(assnsGTruth));
486  }
487  if(fParams.copyBSim()){
488  evt.put(std::move(dk2nucol));
489  evt.put(std::move(nuchoicecol));
490 #ifdef PUT_DK2NU_ASSN
491  evt.put(std::move(assnsdk2nu));
492  evt.put(std::move(assnsnuchoice));
493 #endif
494  }
495  }
496 
497  //___________________________________________________________________________
499  {
501  sr.getByToken(fPOTSumToken, potsum);
502  sr.put(std::make_unique<sumdata::POTSum>(*potsum));
503  }
504 
505  //___________________________________________________________________________
507  {
508  const simb::MCNeutrino& mcn = genieTruth.GetNeutrino();
509  const simb::MCParticle& nu = genieTruth.GetNeutrino().Nu();
510  const bool iscc = mcn.CCNC() == simb::kCC;
511  // TODO TODO - when fixing weights etc in GiBUU simulation, NC != NCbar
512  const int pdg = iscc ? nu.PdgCode() : 14; // all NCs are the same
513 
514  return Key(mcn.Target(), pdg, iscc);
515  }
516 
517  //___________________________________________________________________________
519  GetRecordList(const simb::MCTruth& genieTruth) const
520  {
521  const Key k = GetKey(genieTruth);
522 
523  const simb::MCNeutrino& mcn = genieTruth.GetNeutrino();
524 
525  if(k.iscc && abs(k.nu_pdg) != 12 && abs(k.nu_pdg) != 14){
526  std::cout << "GiBUURegen: No match for CC PDG " << k.nu_pdg
527  << ". Copying GENIE event." << std::endl;
528  return 0;
529  }
530 
531  if(nucleus_to_label.find(k.nucl_pdg) == nucleus_to_label.end()){
532  std::cout << "GiBUURegen: No match for target " << k.nucl_pdg
533  << ". Copying GENIE event." << std::endl;
534  return 0;
535  }
536 
537  if(mcn.Mode() == simb::kCoh ||
539  std::cout << "GiBUURegen: COH or nu-on-e interaction. Copying GENIE event."
540  << std::endl;
541  return 0;
542  }
543 
544  auto it = fRecords.find(k);
545  assert(it != fRecords.end());
546  return it->second;
547  }
548 
549  //___________________________________________________________________________
550  std::vector<TVector3> GiBUURegen::Basis(TVector3 z) const
551  {
552  const TVector3 up(0, 1, 0);
553  const TVector3 x = up.Cross(z).Unit(); // Perpendicular to neutrino and up
554  const TVector3 y = x.Cross(z).Unit(); // Defines the third axis
555 
556  // Let's rotate around the neutrino axis to make reuse less obvious
557  const double a = gRandom->Uniform(0, 2*M_PI);
558 
559  const TVector3 xp = cos(a) * x + sin(a) * y;
560  const TVector3 yp = -sin(a) * x + cos(a) * y;
561  const TVector3 zp = z.Unit();
562 
563  return {xp, yp, zp};
564  }
565 
566  //___________________________________________________________________________
567  void GiBUURegen::Kinematics(const TLorentzVector& nu,
568  const TLorentzVector& lep,
569  double& W, double& x, double& y, double& Qsq) const
570  {
571  const TLorentzVector q = nu-lep;
572 
573  const TLorentzVector p_p(0, 0, 0, .9315); // parton
574 
575  Qsq = -q.Mag2();
576 
577  W = (q + p_p).Mag();
578 
579  x = Qsq / (2*p_p.Dot(q));
580 
581  y = 1-lep.E()/nu.E(); // Identical results, no dependence on p_nuc
582  }
583 
584  //___________________________________________________________________________
586  {
587  auto it = fGenieXSec.find(Key(mcn.Target(),
588  mcn.Nu().PdgCode(),
589  mcn.CCNC() == simb::kCC));
590  assert(it != fGenieXSec.end());
591  return it->second.Eval(mcn.Nu().E());
592  }
593 
594  //___________________________________________________________________________
596  int& trackId) const
597  {
598  // Have to do this in order to renumber the trackIDs properly
599 
602 
603  const int nuTrackId = trackId; // neutrino is always first
604 
605  for(int i = 0; i < genieTruth.NParticles(); ++i){
606  const simb::MCParticle& p = genieTruth.GetParticle(i);
607 
608  simb::MCParticle p_out(trackId++, p.PdgCode(), p.Process(),
609  (p.Mother() == -1) ? -1 : nuTrackId,
610  p.Mass(), p.StatusCode());
611  p_out.AddTrajectoryPoint(p.Position(), p.Momentum());
612 
613  if(abs(p.PdgCode()) == 12 ||
614  abs(p.PdgCode()) == 14 ||
615  abs(p.PdgCode()) == 16){
616  // Stick the event weight on the neutrino
617  p_out.SetWeight(1);
618  }
619 
620  ret.Add(p_out);
621  } // end for i
622 
623  const simb::MCNeutrino& mcn = genieTruth.GetNeutrino();
624 
625  ret.SetNeutrino(mcn.CCNC(), mcn.Mode(), mcn.InteractionType(),
626  mcn.Target(), mcn.HitNuc(), mcn.HitQuark(),
627  mcn.W(), mcn.X(), mcn.Y(), mcn.QSqr());
628 
629  return ret;
630  }
631 
632  //___________________________________________________________________________
634  int& trackId) const
635  {
636  const simb::MCNeutrino& mcn = genieTruth.GetNeutrino();
637  const simb::MCParticle& nu = mcn.Nu();
638 
639  const IRecordList* recs = GetRecordList(genieTruth);
640  if(!recs){
641  // The reason has already been printed
642  ++fNGENIE;
643  return CopyGenieEvent(genieTruth, trackId);
644  }
645 
646  const Record* rec = recs->GetRecord(nu.E());
647 
648  if(!rec){
649  std::cout << "GiBUURegen: No match for energy " << nu.E() << " GeV? "
650  << "Copying GENIE event." << std::endl;
651  ++fNGENIE;
652  return CopyGenieEvent(genieTruth, trackId);
653  }
654 
655  const double genie_xsec = GetGenieXSec(mcn);
656  const int tgtA = (mcn.Target()%1000)/10;
657  TH1* hcorr = fGiBUUCorr.find(GetKey(genieTruth))->second;
658  const double gibuu_corr = tgtA * hcorr->GetBinContent(hcorr->FindBin(nu.E()));
659 
662 
663  const int nuTrackId = trackId; // neutrino is always first
664  // Copy the incoming neutrino and the nucleus into the output record
665  for(int i = 0; i < genieTruth.NParticles(); ++i){
666  const simb::MCParticle& p = genieTruth.GetParticle(i);
667  if(p.Mother() != -1) continue; // Only want to copy the initial state particles
668 
669  simb::MCParticle p_out(trackId++, p.PdgCode(), p.Process(), -1,
670  p.Mass(), p.StatusCode());
671  p_out.AddTrajectoryPoint(p.Position(), p.Momentum());
672 
673  if(abs(p.PdgCode()) == 12 ||
674  abs(p.PdgCode()) == 14 ||
675  abs(p.PdgCode()) == 16){
676  // Stick the event weight on the neutrino
677  p_out.SetWeight(rec->weight * gibuu_corr / genie_xsec);
678  }
679 
680  ret.Add(p_out);
681  } // end for i
682 
683  const std::vector<TVector3> basis = Basis(nu.Momentum().Vect());
684 
685  TLorentzVector p_lep;
686 
687  // Then all the particles from GiBUU
688  for(const Particle& part: rec->parts){
689  simb::MCParticle mcp(trackId++, part.pdg, "primary", nuTrackId);
690 
691  mcp.AddTrajectoryPoint(nu.Position(),
692  TLorentzVector(part.px*basis[0] + part.py*basis[1] + part.pz*basis[2],
693  part.E));
694 
695  ret.Add(mcp);
696 
697  if(abs(part.pdg) >= 11 && abs(part.pdg) <= 16) p_lep = mcp.Momentum();
698  }
699 
700 
701  double W, x, y, Qsq;
702  Kinematics(nu.Momentum(), p_lep, W, x, y, Qsq);
703 
705  mcn.Target(), -1, -1,
706  W, x, y, Qsq);
707 
708  ++fNGiBUU;
709  return ret;
710  }
711 
713 } // namespace
double E(const int i=0) const
Definition: MCParticle.h:232
Atom< bool > copyBSim
const TLorentzVector & Position(const int i=0) const
Definition: MCParticle.h:218
static bool CreateAssn(art::EDProducer const &prod, art::Event &evt, std::vector< T > &a, art::Ptr< U > b, art::Assns< T, U > &assn, size_t indx=UINT_MAX, std::string const &instance=std::string())
Create a 1 to 1 association between a new product and one already in the event.
std::map< Key, TH1 * > fGiBUUCorr
int PdgCode() const
Definition: MCParticle.h:211
int CCNC() const
Definition: MCNeutrino.h:148
double QSqr() const
Definition: MCNeutrino.h:157
const simb::MCNeutrino & GetNeutrino() const
Definition: MCTruth.h:76
GiBUURegenParams fParams
int status
Definition: fabricate.py:1613
Atom< bool > copyMCFlux
void AddTrajectoryPoint(TLorentzVector const &position, TLorentzVector const &momentum)
Definition: MCParticle.h:256
set< int >::iterator it
Atom< bool > lowMemory
int HitQuark() const
Definition: MCNeutrino.h:153
void SetOrigin(simb::Origin_t origin)
Definition: MCTruth.h:80
Atom< std::string > generatorLabel
int Mother() const
Definition: MCParticle.h:212
const art::ProductToken< std::vector< simb::MCTruth > > fTruthToken
const char * p
Definition: xmltok.h:285
const simb::MCParticle & Nu() const
Definition: MCNeutrino.h:146
double Mass() const
Definition: MCParticle.h:238
const IRecordList * GetRecordList(const simb::MCTruth &genieTruth) const
Atom< bool > copyGTruth
int HitNuc() const
Definition: MCNeutrino.h:152
OStream cerr
Definition: OStream.cxx:7
caf::StandardRecord * rec
Definition: tutCAFMacro.C:20
simb::MCTruth CopyGenieEvent(const simb::MCTruth &genieTruth, int &trackId) const
Key GetKey(const simb::MCTruth &genieTruth) const
void abs(TH1 *hist)
int StatusCode() const
Definition: MCParticle.h:210
void Kinematics(const TLorentzVector &nu, const TLorentzVector &lep, double &W, double &x, double &y, double &Qsq) const
std::vector< Particle > parts
Definition: RecordList.h:33
int NParticles() const
Definition: MCTruth.h:74
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
DEFINE_ART_MODULE(TestTMapFile)
std::string Process() const
Definition: MCParticle.h:214
double GetGenieXSec(const simb::MCNeutrino &mcn) const
void Add(simb::MCParticle &part)
Definition: MCTruth.h:79
#define M_PI
Definition: SbMath.h:34
object containing MC flux information
Key(int _nucl_pdg, int _nu_pdg, bool _iscc)
ProductID put(std::unique_ptr< PROD > &&product)
Definition: Event.h:102
virtual const Record * GetRecord(float E) const =0
void endSubRun(art::SubRun &sr)
int InteractionType() const
Definition: MCNeutrino.h:150
TString part[npart]
Definition: Style.C:32
int iscc
std::map< Key, const IRecordList * > fRecords
double W() const
Definition: MCNeutrino.h:154
void beginJob()
double Y() const
Definition: MCNeutrino.h:156
std::string ExpandLibraryPath() const
const double a
std::string getenv(std::string const &name)
Sequence< double > boundingBoxHi
int evt
Sequence< double > boundingBoxLo
caf::StandardRecord * sr
double X() const
Definition: MCNeutrino.h:155
Atom< std::string > libraryPath
std::map< Key, TGraph > fGenieXSec
z
Definition: test.py:28
This class describes a particle created in the detector Monte Carlo simulation.
GiBUURegen(const Parameters &params)
OStream cout
Definition: OStream.cxx:6
const simb::MCParticle & GetParticle(int i) const
Definition: MCTruth.h:75
n.b.: this group is similar but not quite, entirely unlike GENIE ScatteringType convention ...
Definition: MCNeutrino.h:80
int Target() const
Definition: MCNeutrino.h:151
T sin(T number)
Definition: d0nt_math.hpp:132
const art::ProductToken< sumdata::POTSum > fPOTSumToken
TDirectory * dir
Definition: macro.C:5
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
bool operator<(const Key &k) const
ProductID put(std::unique_ptr< PROD > &&)
const TLorentzVector & Momentum(const int i=0) const
Definition: MCParticle.h:219
float Mag() const
T cos(T number)
Definition: d0nt_math.hpp:78
TODO.
assert(nhit_max >=nhit_nbins)
Service to store calibration data products (CDP) in the SQLite3 metadatabase of a file...
Definition: FillParentInfo.h:8
simb::MCTruth GetEvent(const simb::MCTruth &genieTruth, int &trackId) const
void produce(art::Event &evt)
void Format(TGraph *gr, int lcol, int lsty, int lwid, int mcol, int msty, double msiz)
Definition: Style.cxx:154
Event generator information.
Definition: MCTruth.h:32
bool getByToken(ProductToken< PROD > const &token, Handle< PROD > &result) const
Definition: DataViewImpl.h:387
float weight
Definition: RecordList.h:31
Event generator information.
Definition: MCNeutrino.h:18
int Mode() const
Definition: MCNeutrino.h:149
#define W(x)
std::vector< TVector3 > Basis(TVector3 z) const
const art::ProductToken< sumdata::SpillData > fSpillToken
def sign(x)
Definition: canMan.py:197
Beam neutrinos.
Definition: MCTruth.h:23
enum BeamMode string