GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
EvtLibRecordList.cxx
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // \author Chris Backhouse -- c.backhouse@ucl.ac.uk
3 ////////////////////////////////////////////////////////////////////////
4 
6 
8 
9 #include "TFile.h"
10 #include "TTree.h"
11 
12 #include <algorithm>
13 #include <iostream>
14 
15 namespace genie{
16 namespace evtlib{
17  //---------------------------------------------------------------------------
19  {
20  }
21 
22  //---------------------------------------------------------------------------
23  EvtLibRecord::EvtLibRecord(float _E, int _prod_id,
24  const std::vector<EvtLibParticle>& _ps)
25  : E(_E), prod_id(_prod_id), parts(_ps)
26  {
27  }
28 
29  //---------------------------------------------------------------------------
30  bool EvtLibRecord::operator<(const EvtLibRecord& rhs) const
31  {
32  return E < rhs.E;
33  }
34 
35  //---------------------------------------------------------------------------
37  : fTree(tree)
38  {
39  fTree->SetBranchAddress("Enu", &Enu);
40  fTree->SetBranchAddress("prod_id", &prod_id);
41 
42  fTree->SetBranchAddress("nparts", &nparts);
43 
44  fTree->SetBranchAddress("pdg", &pdgs);
45  fTree->SetBranchAddress("E", &Es);
46  fTree->SetBranchAddress("px", &px);
47  fTree->SetBranchAddress("py", &py);
48  fTree->SetBranchAddress("pz", &pz);
49  }
50 
51  //---------------------------------------------------------------------------
53  {
54  }
55 
56  //---------------------------------------------------------------------------
58  {
59  return fTree->GetEntries();
60  }
61 
62  //---------------------------------------------------------------------------
64  {
65  fTree->GetEntry(i);
66 
67  // The event library ROOT format is as minimalistic as possible. We use
68  // variable-sized arrays for the particle list. Because the address of
69  // `parts` had to be provided to SetBranchAddress up-front we had to pick
70  // a fixed size. If the list was actually longer we're in trouble and
71  // should bail out. This seems extremely unlikely to be a problem in
72  // practice, and we can always increase the constant to something larger.
73  if(nparts > kEvtLibMaxParts){
74  LOG("ELI", pFATAL) << "Too many particles " << nparts << "(limit " << kEvtLibMaxParts << ")";
75  exit(1);
76  }
77 
78  std::vector<EvtLibParticle> parts(nparts);
79  for(int j = 0; j < nparts; ++j){
80  parts[j].pdg = pdgs[j];
81  parts[j].E = Es[j];
82  parts[j].px = px[j];
83  parts[j].py = py[j];
84  parts[j].pz = pz[j];
85  } // end for j
86 
87  return EvtLibRecord(Enu, prod_id, parts);
88  }
89 
90  //---------------------------------------------------------------------------
91  SimpleRecordList::SimpleRecordList(TTree* tree, const std::string& prettyName)
92  {
93  std::cout << "Loading " << prettyName;
94  RecordLoader loader(tree);
95 
96  const int N = loader.NRecords();
97  fRecs.reserve(N);
98  for(int i = 0; i < N; ++i){
99  if(i%(N/8) == 0) std::cout << "." << std::flush;
100 
101  fRecs.push_back(loader.GetRecord(i));
102  } // end for i
103  std::cout << std::endl;
104 
105  std::sort(fRecs.begin(), fRecs.end());
106  }
107 
108  //---------------------------------------------------------------------------
110  {
111  auto it = std::lower_bound(fRecs.begin(), fRecs.end(),
112  EvtLibRecord(E, 0, {}));
113  if(it == fRecs.end()) return 0;
114  return &(*it);
115  }
116 
117  //---------------------------------------------------------------------------
118  OnDemandRecordList::OnDemandRecordList(TTree* tree, const std::string& prettyName)
119  : fTree(tree), fPrettyName(prettyName), fLoader(0)
120  {
121  }
122 
123  //---------------------------------------------------------------------------
125  {
126  std::cout << "Loading index to " << fPrettyName;
127 
128  float Enu;
129  fTree->SetBranchAddress("Enu", &Enu);
130 
131  const int N = fTree->GetEntries();
132  fEnergies.reserve(N);
133 
134  for(int i = 0; i < N; ++i){
135  if(i%(N/8) == 0) std::cout << "." << std::flush;
136  fTree->GetEntry(i);
137 
138  fEnergies.emplace_back(Enu, i);
139  } // end for i
140  std::cout << std::endl;
141 
142  std::sort(fEnergies.begin(), fEnergies.end());
143  }
144 
145  //---------------------------------------------------------------------------
147  {
148  if(fEnergies.empty()){
149  LoadIndex();
150  // The loader must be created after the indices are loaded, because they
151  // both want to set branch addresses in the TTree, and we ultimately need
152  // the loader's to be active.
153  fLoader = new RecordLoader(fTree);
154  }
155 
156  auto it = std::lower_bound(fEnergies.begin(), fEnergies.end(),
157  std::make_pair(E, 0));
158  if(it == fEnergies.end()) return 0;
159 
160  fRecord = fLoader->GetRecord(it->second);
161 
162  return &fRecord;
163  }
164 }} // namespaces
const int kEvtLibMaxParts
Maximum number of particles supported in a single library event record.
EvtLibRecord GetRecord(int i) const
float px[kEvtLibMaxParts]
#define pFATAL
Definition: Messenger.h:56
float Es[kEvtLibMaxParts]
const EvtLibRecord * GetRecord(float E) const override
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
OnDemandRecordList(TTree *tree, const std::string &prettyName)
float py[kEvtLibMaxParts]
std::vector< EvtLibRecord > fRecs
bool operator<(const EvtLibRecord &rhs) const
Order by energy - this allows OnDemandRecordList to work efficiently.
const EvtLibRecord * GetRecord(float E) const override
float pz[kEvtLibMaxParts]
SimpleRecordList(TTree *tree, const std::string &prettyName)
int pdgs[kEvtLibMaxParts]
std::vector< std::pair< float, int > > fEnergies