eventListToTextFile.C
Go to the documentation of this file.
1 #include "TTree.h"
2 #include "TFile.h"
3 #include <iostream>
4 #include <fstream>
5 #include <set>
6 #include <map>
7 
8 class EventId{
9 
10 public:
12  : run (std::numeric_limits<int>::lowest())
13  , subrun(std::numeric_limits<int>::lowest())
14  , event (std::numeric_limits<int>::lowest())
15  , cycle (std::numeric_limits<int>::lowest())
16  , slice (std::numeric_limits<int>::lowest())
17  {}
18 
19  int run;
20  int subrun;
21  int event;
22  int cycle;
23  int slice;
24 };
25 
26 struct WeightVars{
27 
29  : fXSecCVPPFX_Weight (1.)
30  , fRPACCQE_Weight (1.)
31  , fRPARES_Weight (1.)
32  , fDISvnCC1pi_Weight (1.)
33  , fEmpiricalMEC_Weight (1.)
34  , fEmpiricalMECtoGENIEQE_Weight (1.)
35  , fEmpiricalMECtoGENIERES_Weight(1.)
36  , fOtherDIS_Weight (1.)
37  {}
38 
39  float fXSecCVPPFX_Weight; ///< Was Tufts weight for SA
40  float fRPACCQE_Weight; ///<
41  float fRPARES_Weight; ///< To be used for systematic evaluation ONLY
42  float fDISvnCC1pi_Weight; ///< To be used for systematic evaluation ONLY
43  float fEmpiricalMEC_Weight; ///< Ditto
46  float fOtherDIS_Weight; ///< Ditto
47 };
48 
49 struct DataVars{
50 
51  // initialize the Fake_Weight to be 1, it will be reset if Fake_Data are generated
53  : fNuE_CVN (std::numeric_limits<float>::lowest())
54  , fNuE_NumMichel (std::numeric_limits<float>::lowest())
55  , fHad_RecoE (std::numeric_limits<float>::lowest())
56  , fLep_RecoE (std::numeric_limits<float>::lowest())
57  , fLep_RecoE_MCFrac(0.)
58  , fRecoQ2 (std::numeric_limits<float>::lowest())
59  , fFake_Weight (1.)
60  {}
61 
62  float fNuE_CVN; ///<
63  float fNuE_NumMichel; ///<
64  float fHad_RecoE; ///<
65  float fLep_RecoE; ///<
66  float fLep_RecoE_MCFrac; ///< fraction of leptonic energy in muon catcher
67  float fRecoQ2; ///< reconstructed Q^2
68  float fFake_Weight; ///< Weight for fake data events
69 };
70 
71 struct SelectOutput {
72 
74  : selName("noname")
75  , energy (std::numeric_limits<double>::lowest())
76  , cvn (std::numeric_limits<double>::lowest())
77  , wgt (1.)
78  , run (std::numeric_limits<int>::lowest())
79  , subrun (std::numeric_limits<int>::lowest())
80  , event (std::numeric_limits<int>::lowest())
81  , slice (std::numeric_limits<int>::lowest())
82  , cycle (std::numeric_limits<int>::lowest())
83  {}
84 
85 
87  double e,
88  double id,
89  double w,
90  int r,
91  int sr,
92  int ev,
93  int sl,
94  int cy)
95  : selName(sel)
96  , energy (e)
97  , cvn (id)
98  , wgt (w)
99  , run (r)
100  , subrun (sr)
101  , event (ev)
102  , slice (sl)
103  , cycle (cy)
104  {}
105 
107  double energy;
108  double cvn;
109  double wgt;
110  int run;
111  int subrun;
112  int event;
113  int slice;
114  int cycle;
115 };
116 
117 //-----------------------------------------------------------------------------
118 double nuEEnergy(double em,
119  double had,
120  bool fhc)
121 {
122  // p0 and p3 are 0.0 for both FHC and RCH
123  double p0 = 0.0;
124  double p1 = (fhc) ? 1.00756 : 0.980479;
125  double p2 = (fhc) ? 1.07093 : 1.45170;
126  double p3 = 0.0;
127  double p4 = (fhc) ? 1.28608e-2 : -5.82609e-3;
128  double p5 = (fhc) ? 2.27129e-1 : -2.27599e-1;
129  double fr = (fhc) ? 1. / 1.0501206 : 1. / 1.001766;
130 
131  return fr *(( p5 * had * had ) +
132  ( p4 * em * em ) +
133  ( p3 * em * had ) +
134  ( p2 * had ) +
135  ( p1 * em ) +
136  p0 );
137 }
138 
139 //-----------------------------------------------------------------------------
140 void treeContentsToFile(std::vector<SelectOutput> & outputVec,
141  TFile* rootFile,
142  std::string const& treeName,
143  std::string const& selName)
144 {
145  // get the TTree
146  TTree *tree = dynamic_cast<TTree *>(rootFile->Get(treeName.c_str()));
147 
148  // std::cout << treeName
149  // << " "
150  // << tree
151  // << std::endl;
152  // std::cout << " with "
153  // << tree->GetEntriesFast()
154  // << " events"
155  // << std::endl;
156 
157  // set the locations of the variables we want
158  EventId evId;
159  WeightVars wgtVars;
160  DataVars dataVars;
161 
162  tree->SetBranchAddress("eventId", &evId);
163  tree->SetBranchAddress("dataVars", &dataVars);
164 
165  bool isMC = (treeName.find("MC") != std::string::npos);
166  if(isMC) tree->SetBranchAddress("weightVars", &wgtVars);
167 
168  bool isFHC = (treeName.find("FHC") != std::string::npos);
169 
170  double energy;
171 
172  for(int i = 0; i < tree->GetEntriesFast(); ++i){
173  tree->GetEntry(i);
174 
175  energy = dataVars.fLep_RecoE + dataVars.fHad_RecoE;
176  if(selName.find("NuE") != std::string::npos)
177  energy = nuEEnergy(dataVars.fLep_RecoE,
178  dataVars.fHad_RecoE,
179  isFHC);
180  else if(selName.find("NC") != std::string::npos)
181  energy = dataVars.fHad_RecoE;
182 
183  outputVec.emplace_back(selName,
184  energy,
185  dataVars.fNuE_CVN,
186  wgtVars.fXSecCVPPFX_Weight,
187  evId.run,
188  evId.subrun,
189  evId.event,
190  evId.slice,
191  evId.cycle);
192  }
193 }
194 
195 //-----------------------------------------------------------------------------
197 {
198  // The plan here is to open the file with the event list trees
199  // and loop over each combination of detector/selection/horn current/file type
200  // to produce output text files of the event lists that can be used to compare
201  // CMF and CAF event lists
202  std::map<std::string, std::string> detectors({
203  {"Near", "nd"}
204  ,{"Far", "fd"}
205  });
206  std::set<std::string> epochs({ "Epoch1FHC"
207  , "Epoch2FHC"
208  , "Epoch3FHC"
209  , "Epoch4RHC"
210  , "Epoch5FHC"
211  , "Epoch6RHC"});
212  std::map<std::string, std::string> fileTypes({
213  {"Data", "data"}
214  ,{"Beam", "nonswap"}
215  ,{"FluxSwap", "fluxswap"}
216  ,{"TauSwap", "tauswap"}
217  ,{"CosmicBackground", "cosmicbackground"}
218  ,{"RockNonSwap", "rocknonswap"}
219  ,{"RockFluxSwap", "rockfluxswap"}
220  });
221  std::map<std::string, std::set<std::string> > selections({
222  {"nue", {"NuESel_LowPID", "NuESel_HighPID", "NuESel_Peripheral"}}
223  ,{"numu", {"NuMuSelQ1", "NuMuSelQ2", "NuMuSelQ3", "NuMuSelQ4"}}
224  ,{"nus", {"NCSel"}}
225  });
226  std::set<std::string> intTypes({ "NuMuCC"
227  , "NuMuBarCC"
228  , "NuECC"
229  , "NuEBarCC"
230  , "NuTauCC"
231  , "NuTauBarCC"
232  , "NC"
233  , "CosmicMuon"});
234 
235  std::map<std::string, std::vector<SelectOutput> > selectedEvents;
236 
237  // Get the file containing all event list trees
238  TFile *eventListFile = new TFile("/nova/ana/users/brebel/skimmed/cmf_eventlisttrees_2019Sep06.root", "READ");
239 
240  std::string pathInFile("list/full/");
242  std::string treeName;
243  std::string dataMC;
244  std::string hc;
245  std::string hcFN;
246 
247  // loop over the combinations of strings from the sets above to grab the trees from the file
248  for(auto const& det : detectors){
249 
250  for(auto const& ft : fileTypes){
251  // certain detector/file type combinations aren't allowed, skip them here
252  if(det.first == "Near" &&
253  ft.first != "Data" &&
254  ft.first != "Beam") continue;
255 
256  dataMC = (ft.first == "Data") ? ft.first : "MC";
257 
258  for(auto const& selItr : selections){
259  // certain selection/filetype combinations aren't allowed
260  if(selItr.first != "nue" &&
261  ft.first.find("Rock") != std::string::npos) continue;
262 
263  // no cosmic background files for NC
264  if(selItr.first == "nus" &&
265  ft.first.find("Cosmic") != std::string::npos) continue;
266 
267 
268  // the epoch strings also have FHC or RHC in them
269  for(auto const& epoch : epochs){
270 
271  hc = (epoch.find("RHC") == std::string::npos) ? "fhc" : "rhc";
272 
273  // std::cout << epoch
274  // << " "
275  // << ft.first
276  // << " "
277  // << hc
278  // << " "
279  // << selItr.first
280  // << std::endl;
281 
282  outName = "selectedEvents_" + det.second + "_" + selItr.first + "_" + hc + "_" + ft.second + "_2018.txt";
283  if(selectedEvents.count(outName) < 1) selectedEvents[outName] = std::vector<SelectOutput>(0);
284 
285  for(auto const& sel : selItr.second){
286 
287  // no peripheral selection in the ND
288  if(sel.find("Peripheral") != std::string::npos &&
289  det.first == "Near") continue;
290 
291  treeName = dataMC + det.first + epoch;
292 
293  if(ft.first != "Data"){
294  treeName += ft.first + sel;
295  for(auto const& it : intTypes){
296  if((ft.first == "CosmicBackground" &&
297  it != "CosmicMuon") ||
298  (ft.first != "CosmicBackground" &&
299  it == "CosmicMuon")) continue;
300 
301  treeContentsToFile(selectedEvents.find(outName)->second,
302  eventListFile,
303  pathInFile + treeName + it,
304  sel);
305  } // end if we have to worry about the interaction types
306  } // end if MC
307  else{
308  treeName += sel;
309  treeContentsToFile(selectedEvents.find(outName)->second,
310  eventListFile,
311  pathInFile + treeName,
312  sel);
313  } // end if data
314 
315  } // end loop over specific selections
316 
317  } // end loop over epochs
318  } // end loop over generic selections
319  } // end loop over file types
320  } // end loop over detectors
321 
322  // write the information out to files
323  for(auto const& itr : selectedEvents){
324  std::ofstream outFile(itr.first);
325 
326  std::cout << "creating file: "
327  << itr.first
328  << " with "
329  << itr.second.size()
330  << " entries"
331  << std::endl;
332 
333  for(auto const& sitr : itr.second){
334  outFile << sitr.selName
335  << "\t\t"
336  << sitr.run
337  << "\t"
338  << sitr.subrun
339  << "\t"
340  << sitr.event
341  << "\t"
342  << sitr.slice
343  << "\t"
344  << sitr.cycle
345  << "\t"
346  << sitr.energy
347  << "\t"
348  << sitr.cvn
349  << "\t"
350  << sitr.wgt
351  << std::endl;
352  }
353 
354  outFile.close();
355 
356  } // end loop over map to selected information
357 
358  // done!
359 }
float fOtherDIS_Weight
Ditto.
void treeContentsToFile(std::vector< SelectOutput > &outputVec, TFile *rootFile, std::string const &treeName, std::string const &selName)
float fEmpiricalMECtoGENIERES_Weight
Ditto.
float fLep_RecoE_MCFrac
fraction of leptonic energy in muon catcher
set< int >::iterator it
std::string selName
float fRPARES_Weight
To be used for systematic evaluation ONLY.
float fEmpiricalMECtoGENIEQE_Weight
Ditto.
float fEmpiricalMEC_Weight
Ditto.
Defines an enumeration for prong classification.
SelectOutput(std::string const &sel, double e, double id, double w, int r, int sr, int ev, int sl, int cy)
double nuEEnergy(double em, double had, bool fhc)
float fXSecCVPPFX_Weight
Was Tufts weight for SA.
TFile * outFile
Definition: PlotXSec.C:135
float fRecoQ2
reconstructed Q^2
double energy
Definition: plottest35.C:25
void eventListToTextFile()
float fDISvnCC1pi_Weight
To be used for systematic evaluation ONLY.
caf::StandardRecord * sr
const ana::Var wgt
static bool isFHC
Definition: run.py:1
OStream cout
Definition: OStream.cxx:6
TRandom3 r(0)
Float_t e
Definition: plot.C:35
Float_t w
Definition: plot.C:20
float fFake_Weight
Weight for fake data events.
float fNuE_NumMichel
enum BeamMode string