GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
EventLibraryInterface.cxx
Go to the documentation of this file.
1 //____________________________________________________________________________
2 /*
3  Copyright (c) 2003-2024, The GENIE Collaboration
4  For the full text of the license visit http://copyright.genie-mc.org
5 
6 
7 */
8 //____________________________________________________________________________
9 
20 #include "Tools/EvtLib/Utils.h"
22 
23 #include "TFile.h"
24 
25 using namespace genie;
26 using namespace genie::evtlib;
27 
28 //____________________________________________________________________________
29 EventLibraryInterface::EventLibraryInterface() :
30  EventRecordVisitorI("genie::evtlib::EventLibraryInterface"),
31  fRecordFile(0)
32 {
33 
34 }
35 
36 //____________________________________________________________________________
38  EventRecordVisitorI("genie::evtlib::EventLibraryInterface", config),
39  fRecordFile(0)
40 {
41 
42 }
43 
44 //____________________________________________________________________________
46 {
47  Cleanup();
48 }
49 
50 //____________________________________________________________________________
52 {
53 // Get event summary constructed by GENIE
54 //
55  Interaction* interaction = event->Summary();
56  const InitialState & init_state = interaction->InitState();
57 
58  const EvtLibRecord* rec = GetRecord(interaction);
59  if(!rec) return; // Reason has already been printed
60 
61 
62  // at this point we have the event to add to our event record
63  // the energy might be not the same
64  // so first we correct the probe:
65  // we maintain the same direction but we set the energy selected
66  // by the record.
67  // We only change the the particle in the event record,
68  // Not the summary so that we keep the history
69  // of the event construction
70 
71  // Neutrino is a parent to the lepton(s)
72  GHepParticle * probe = event->Particle(0) ;
73 
74  TLorentzVector probe_p4 ( * probe->P4() ) ;
75 
76  LOG("ELI", pINFO) << "Difference between requested neutrino E and used neutrino energy: " << probe_p4.E() - rec -> E ;
77 
78  probe_p4 *= rec -> E / probe_p4.E() ;
79  probe -> SetMomentum( probe_p4 ) ;
80 
81  // because of the selection procedure there might not be enough
82  // events in the library
83  // so we add a roation around the z axis which is random for each event
84  double alpha = RandomGen::Instance()->RndEvg().Uniform(0, 2*constants::kPi );
85 
86 
87  // now we simply have to add the particles from the EvtLibRecord
88  // But we need to rotate the particles along the direction of the beam
89  // so we evaluate this direction once and for all
90  TVector3 unit_nudir = probe->P4()->Vect().Unit();
91 
92 
93  int firstLep = -1 ;
94  const auto & parts = rec -> parts ;
95 
96  for ( unsigned int i = 0; i < parts.size() ; ++i ) {
97 
98  const auto & part = parts[i] ;
99 
100  int pdg = part.pdg;
101 
102  if(pdg::IsLepton(part.pdg)) {
103  if ( firstLep == -1 ) {
104  firstLep = i ;
105  interaction -> KinePtr() -> SetFSLeptonP4( part.px, part.py, part.pz, part.E ) ;
106  }
107  }
108 
109  // Fix up the outgoing lepton for NC events (due to lepton universality
110  // it could be any value in the library)
111 
112  if( interaction->ProcInfo().IsWeakNC() ) {
113  if ( int(i) == firstLep ) {
114  // in the case of NC the process is always stored as numu (or numubar)
115  // as the couplings are the same
116  // here we restore the proper pdg on what we think it is the leading neutrino
117  pdg = init_state.ProbePdg();
118  }
119  }
120 
121  TLorentzVector p4(part.px, part.py, part.pz, part.E) ;
122 
123  // The next is the rotation of the event according to the common alpha
124  p4.RotateZ( alpha ) ;
125 
126  // Now we rotate the events so that it matches the direction of the incoming neutrino
127  p4.RotateUz(unit_nudir);
128 
129  event->AddParticle(pdg,
131  0 , 1, -1, -1, // child of the neutrino and nucleus
132  p4,
133  TLorentzVector(0, 0, 0, 0));
134 
135  }
136 
137  // we now have the particle in the right direction inside the record we just need to
138  // decorate the summary kinematics
139 
140  if ( firstLep < 0 ) {
141  firstLep = 2 ;
142  }
143  else {
144  firstLep += 2 ;
145  }
146  // these are events from outside, there are no assumptions on the lepton
147  // If we are simulating BSM physics there might not even be a lepton in the final state
148  // So if there is a lepton we evaluate Q2 based on that lepton
149  // If not we evaluate Q2 using the first particle we have available
150 
151  FillKinematics( *event, *interaction->KinePtr(), firstLep );
152 
153 }
154 
155 //____________________________________________________________________________
157 GetRecord(const Interaction* interaction) const
158 {
159  const InitialState& init_state = interaction->InitState();
160 
161  const double probe_E = init_state.ProbeE(kRfLab);
162 
163  if(!init_state.Tgt().IsNucleus()){
164  LOG("ELI", pINFO) << "Skippping non-nuclear target " << init_state;
165  return 0;
166  }
167 
168  const int tgt_pdgc = init_state.TgtPdg();
169 
170  const ProcessInfo& proc = interaction->ProcInfo();
171 
172  if(!proc.IsWeakCC() && !proc.IsWeakNC()){
173  LOG("ELI", pINFO) << "Skipping unknown process " << proc;
174  return 0;
175  }
176 
177  int probe_pdgc = init_state.ProbePdg();
178 
179  // Use nu_mu for NC as a convention internal to this code to index into the
180  // records map.
181  if(proc.IsWeakNC()){
182  if(probe_pdgc > 0) probe_pdgc = kPdgNuMu;
183  else probe_pdgc = kPdgAntiNuMu;
184  }
185 
186  const Key key(tgt_pdgc, probe_pdgc, proc.IsWeakCC());
187 
188  const auto rec_it = fRecords.find(key);
189 
190  if(rec_it == fRecords.end()){
191  LOG("ELI", pINFO) << "Skippping " << key << " -- not found in library";
192  return 0;
193  }
194 
195  const EvtLibRecord* rec = rec_it->second->GetRecord(probe_E);
196 
197  if(!rec){
198  LOG("ELI", pINFO) << "Skippping " << key << " at " << probe_E << " GeV -- not found in library";
199  return 0;
200  }
201 
202  return rec;
203 }
204 //____________________________________________________________________________
206 {
207  Algorithm::Configure(config);
208  LoadRecords();
209 }
210 
211 //___________________________________________________________________________
213 {
214  Algorithm::Configure(config);
215  LoadRecords();
216 }
217 
218 //___________________________________________________________________________
220 {
221  for(auto it: fRecords) delete it.second;
222  fRecords.clear();
223  delete fRecordFile;
224  fRecordFile = 0;
225 }
226 
227 //___________________________________________________________________________
229 {
230  Cleanup();
231 
232  std::string libPath;
233  GetParam("EventLibraryPath", libPath);
234  Expand(libPath);
235 
236  bool onDemand;
237  GetParam("OnDemand", onDemand);
238 
239  PDGLibrary* pdglib = PDGLibrary::Instance();
240 
241  fRecordFile = new TFile(libPath.c_str());
242  if(fRecordFile->IsZombie()) exit(1);
243 
244  TIter next(fRecordFile->GetListOfKeys());
245  while(TObject* dir = next()){
246  const std::string& tgtName = dir->GetName();
247  const TParticlePDG* tgtPart = pdglib->DBase()->GetParticle(tgtName.c_str());
248  if(!tgtPart){
249  LOG("ELI", pWARN) << "Unknown nucleus " << tgtName
250  << " found in " << libPath
251  << " -- skipping";
252  continue;
253  }
254 
255  for(int pdg: {kPdgNuE, kPdgAntiNuE,
258 
259  for(bool iscc: {true, false}){
260  // NCs should be the same for all flavours. Use nu_mu as a convention
261  // internal to this code to index into the records map.
262  if(!iscc && abs(pdg) != kPdgNuMu) continue;
263 
264  std::string nuName = pdglib->Find(pdg)->GetName();
265  if(!iscc) nuName = pdg::IsAntiNeutrino(pdg) ? "nu_bar" : "nu";
266 
267  const std::string treeName =
268  TString::Format("%s/%s/%s/records",
269  tgtName.c_str(),
270  iscc ? "cc" : "nc",
271  nuName.c_str()).Data();
272 
273  const Key key(tgtPart->PdgCode(), pdg, iscc);
274 
275  TTree* tr = (TTree*)fRecordFile->Get(treeName.c_str());
276 
277  if(!tr){
278  LOG("ELI", pINFO) << treeName << " not found in "
279  << libPath << " -- skipping";
280  continue;
281  }
282 
283  if(onDemand)
284  fRecords[key] = new OnDemandRecordList(tr, treeName);
285  else
286  fRecords[key] = new SimpleRecordList(tr, treeName);
287  } // end for iscc
288  } // end for pdg
289  } // end for dir
290 
291  // Need to keep the record file open for OnDemand, but not Simple
292  if(!onDemand){delete fRecordFile; fRecordFile = 0;}
293 }
294 
295 //___________________________________________________________________________
297  Kinematics& kine,
298  int primary_lep_id ) const {
299 
300 
301  const TLorentzVector & p_probe = * event.Particle( 0 ) -> P4() ;
302 
303  const TLorentzVector & p_lep = * event.Particle( primary_lep_id ) -> P4() ;
304 
305  const TLorentzVector q = p_probe - p_lep;
306 
307  // Initial hadronic state, semi-arbitrary
308  const TLorentzVector & p_tgt = * event.Particle( 1 ) -> P4() ;
309 
310  kine.SetQ2(-q.Mag2(), true);
311  kine.SetW((q + p_tgt).Mag(), true);
312  kine.Setx(-q.Mag2() / (2*p_tgt.Dot(q)), true);
313  kine.Sety( p_tgt.Dot( q ) / p_tgt.Dot( p_probe ), true ) ;
314 
315 }
bool IsWeakCC(void) const
const int kPdgNuE
Definition: PDGCodes.h:28
std::map< Key, const IEvtLibRecordList * > fRecords
static RandomGen * Instance()
Access instance.
Definition: RandomGen.cxx:71
void SetQ2(double Q2, bool selected=false)
Definition: Kinematics.cxx:255
const TLorentzVector * P4(void) const
Definition: GHepParticle.h:78
Kinematics * KinePtr(void) const
Definition: Interaction.h:76
Defines the EventRecordVisitorI interface. Concrete implementations of this interface use the &#39;Visito...
const int kPdgAntiNuE
Definition: PDGCodes.h:29
void Expand(std::string &s)
Expand env vars/homedirs/wildcards in s.
Definition: Utils.cxx:8
bool IsNucleus(void) const
Definition: Target.cxx:272
TDatabasePDG * DBase(void)
Definition: PDGLibrary.cxx:81
void Configure(const Registry &config) override
const int kPdgNuMu
Definition: PDGCodes.h:30
Generated/set kinematical variables for an event.
Definition: Kinematics.h:39
string dir
Summary information for an interaction.
Definition: Interaction.h:56
void ProcessEventRecord(GHepRecord *event) const override
bool IsWeakNC(void) const
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
A class encapsulating an enumeration of interaction types (EM, Weak-CC, Weak-NC) and scattering types...
Definition: ProcessInfo.h:46
bool IsAntiNeutrino(int pdgc)
Definition: PDGUtils.cxx:118
virtual void Configure(const Registry &config)
Definition: Algorithm.cxx:62
int ProbePdg(void) const
Definition: InitialState.h:64
TRandom3 & RndEvg(void) const
rnd number generator used by the event generation drivers
Definition: RandomGen.h:74
#define pINFO
Definition: Messenger.h:62
const int kPdgAntiNuTau
Definition: PDGCodes.h:33
const int kPdgAntiNuMu
Definition: PDGCodes.h:31
#define pWARN
Definition: Messenger.h:60
void Setx(double x, bool selected=false)
Definition: Kinematics.cxx:231
const int kPdgNuTau
Definition: PDGCodes.h:32
void SetW(double W, bool selected=false)
Definition: Kinematics.cxx:279
static PDGLibrary * Instance(void)
Definition: PDGLibrary.cxx:68
Singleton class to load &amp; serve a TDatabasePDG.
Definition: PDGLibrary.h:35
A registry. Provides the container for algorithm configuration parameters.
Definition: Registry.h:65
void Sety(double y, bool selected=false)
Definition: Kinematics.cxx:243
int TgtPdg(void) const
const InitialState & InitState(void) const
Definition: Interaction.h:69
const ProcessInfo & ProcInfo(void) const
Definition: Interaction.h:70
void Format(TGraph *gr, int lcol, int lsty, int lwid, int mcol, int msty, double msiz)
Definition: Style.cxx:146
TParticlePDG * Find(int pdgc, bool must_exist=true)
Definition: PDGLibrary.cxx:86
bool GetParam(const RgKey &name, T &p, bool is_top_call=true) const
const Target & Tgt(void) const
Definition: InitialState.h:66
const EvtLibRecord * GetRecord(const Interaction *interaction) const
double ProbeE(RefFrame_t rf) const
GENIE&#39;s GHEP MC event record.
Definition: GHepRecord.h:45
Most commonly used PDG codes. A set of utility functions to handle PDG codes is provided in PDGUtils...
STDHEP-like event record entry that can fit a particle or a nucleus.
Definition: GHepParticle.h:39
bool IsLepton(int pdgc)
Definition: PDGUtils.cxx:86
void FillKinematics(const GHepRecord &, genie::Kinematics &kine, int primary_lep_id) const
Initial State information.
Definition: InitialState.h:48