CAFToEventList_module.cc
Go to the documentation of this file.
1 //
2 // Created by Brian Rebel on 3/13/19.
3 //
4 // Module to take CAF formatted files and turn them into EventList formats
5 //
6 
7 #include <memory>
8 #include <cmath>
9 #include <string>
10 #include <vector>
11 #include <unordered_map>
12 #include <set>
13 
14 // Framework includes
24 #include "fhiclcpp/ParameterSet.h"
26 #include "cetlib_except/exception.h"
27 #include "cetlib/search_path.h"
28 
29 // NOvA includes
30 #include "SummaryData/POTSum.h"
32 #include "FNEX/core/Event.h"
33 #include "FNEX/core/VarVals.h"
34 #include "nugen/NuReweight/ReweightLabels.h"
37 
40 #include "StandardRecord/SRSpill.h"
41 #include "CAFAna/Cuts/Cuts.h"
42 #include "CAFAna/Cuts/SpillCuts.h"
43 #include "CAFAna/Cuts/NumuCuts2018.h"
44 #include "CAFAna/Cuts/TruthCuts.h"
45 #include "CAFAna/Cuts/TimingCuts.h"
46 #include "CAFAna/Cuts/NueCuts2018.h"
47 #include "CAFAna/Cuts/QuantileCuts.h"
48 #include "CAFAna/Vars/NueEnergy2018.h"
49 #include "CAFAna/Vars/NumuVars.h"
50 #include "CAFAna/Vars/HistAxes.h"
51 #include "CAFAna/Vars/Vars.h"
54 #include "CAFAna/Vars/TruthVars.h"
55 #include "CAFAna/Vars/XsecTunes.h"
56 #include "CAFAna/Systs/XSecSysts.h"
57 
58 #include "TTree.h"
59 #include "TFile.h"
60 #include "TH1.h"
61 #include "TRandom3.h"
62 
63 namespace fnex {
64 
65  struct POTLiveTime{
66 
68  : pot (-1.)
69  , liveTime(-1.)
70  {}
71 
72  POTLiveTime(double const& p,
73  double const& l)
74  : pot (p)
75  , liveTime(l)
76  {}
77 
78  double pot;
79  double liveTime;
80 
81  void operator+=(POTLiveTime const& o) { pot = o.pot; liveTime = o.liveTime; }
82 
83  };
84 
86  public:
87  explicit CAFToEventList(fhicl::ParameterSet const& pset);
88  virtual ~CAFToEventList();
89 
90  // Plugins should not be copied or assigned.
91  CAFToEventList(CAFToEventList const &) = delete;
92  CAFToEventList(CAFToEventList &&) = delete;
93  CAFToEventList & operator = (CAFToEventList const &) = delete;
94  CAFToEventList & operator = (CAFToEventList &&) = delete;
95 
96  void analyze(art::Event const& e) override;
97  void reconfigure(fhicl::ParameterSet const& p) ;
98  void endJob() override;
99 
100  private:
101 
102  void clear ();
103  void FindPOTPerPeriod (TFile & tf);
104  void InsertPeriod6Cosmics ();
105  int PeriodAdjustCosmics (int const& inputPeriod);
106  void InitializeEventListMaps();
107  std::string RunToEpoch (int run);
108  bool PassesSelections (caf::SRProxy const& proxy) const;
109  bool PassesNueSelections (caf::SRProxy const& proxy) const;
110  bool PassesNumuSelections (caf::SRProxy const& proxy) const;
111  bool PassesSpillAndBeamCuts (caf::SRProxy const& proxy) const;
112  void FillVariables (std::vector<fnex::DataVarVals> & dataVals,
113  std::vector<fnex::MCVarVals> & mcVals,
114  std::vector<fnex::EventId> & evIds,
115  std::vector<fnex::SelectionType_t> & sels);
116  void FillMCVarVals (fnex::MCVarVals & mcVarVals,
117  caf::SRProxy & srProxy) const;
118  void FillDataVarVals (fnex::DataVarVals & dataVals,
119  caf::SRProxy const& srProxy) const;
120  void FillEventId (fnex::EventId & evId,
121  caf::SRProxy const& srProxy) const;
122  fnex::SelectionType_t FindSelectionType (caf::SRProxy const& srProxy) const;
123  void DefineNumuQuantiles ();
124 
125  std::vector<std::string> fCAFFileNames; ///< name of the CAF file we are using
126  std::string fFHCQuantFileName; ///< name of file in numudata ups product with FHC quantiles
127  std::string fRHCQuantFileName; ///< name of file in numudata ups product with FHC quantiles
128  std::string fSpillTreeName; ///< names of the CAF histogram with the total POT
129  std::string fPOTHistName; ///< names of the CAF histogram with the total POT
130  std::string fCAFRecTreeName; ///< names of the CAF tree with the reco info
131  novadaq::cnv::DetId fDetector; ///< which detector are we using
132  fnex::EventListMap fMDToEvents; ///< each entry in the vector is for a single slice
133  fnex::FileType_t fFileType; ///< get from the configuration for now
134  bool fIsRealData; ///< are we looking at real data?
135  bool fIsFHC; ///< are we looking at FHC files?
136  std::vector<ana::Cut> fNumuFHCQuantCuts; ///< quantile boundaries for the numu FHC data
137  std::vector<ana::Cut> fNumuRHCQuantCuts; ///< quantile boundaries for the numu RHC data
138  unsigned int fNumQuantiles; ///< number of quantiles
139  //std::set<fnex::EventId> fMissingEvents; ///< events missing compared to CAF sample
140  std::map<int, POTLiveTime> fPeriodToPOTLiveTime; ///< keep track of the POT in each period
141  std::set<fnex::SelectionType_t> fSelections; ///< Selections to use
142  bool fMakeSelectedList; ///< configure to true to make a list
143  };
144 
145  //......................................................................
147  : EDAnalyzer (pset)
148  , fDetector (novadaq::cnv::kNEARDET)
149  , fFileType (fnex::kUnknownFileType)
150  , fIsRealData (false)
151  , fIsFHC (true)
152  , fNumQuantiles (1)
153  , fMakeSelectedList(false)
154  {
155  // make sure everything is zeroed out.
156  this->clear();
157 
158  this->reconfigure(pset);
159 
160  // initialize the period to POT map
161  // just allow for many periods, unlikely we will
162  // have more than 50, but it doesn't matter if we
163  // plan ahead
164  for(int i = 1; i < 51; ++i) fPeriodToPOTLiveTime[i * 1000] = POTLiveTime(0., 0.);
165 
166  // fMissingEvents.insert(fnex::EventId(16313, 30, 156, 0, 38));
167  // fMissingEvents.insert(fnex::EventId(17014, 40, 671, 0, 40));
168  // fMissingEvents.insert(fnex::EventId(17920, 36, 971, 0, 17));
169  // fMissingEvents.insert(fnex::EventId(17968, 10, 891, 0, 30));
170  // fMissingEvents.insert(fnex::EventId(18209, 27, 137, 0, 30));
171  // fMissingEvents.insert(fnex::EventId(18234, 10, 854, 0, 39));
172  // fMissingEvents.insert(fnex::EventId(18374, 4, 665, 0, 35));
173  // fMissingEvents.insert(fnex::EventId(18478, 33, 737, 0, 39));
174  // fMissingEvents.insert(fnex::EventId(18507, 41, 474, 0, 28));
175  // fMissingEvents.insert(fnex::EventId(18521, 8, 645, 0, 30));
176  // fMissingEvents.insert(fnex::EventId(18529, 48, 342, 0, 29));
177  // fMissingEvents.insert(fnex::EventId(18556, 58, 347, 0, 29));
178  // fMissingEvents.insert(fnex::EventId(18812, 29, 367, 0, 45));
179  // fMissingEvents.insert(fnex::EventId(18910, 43, 356, 0, 36));
180  // fMissingEvents.insert(fnex::EventId(19198, 44, 210, 0, 38));
181  // fMissingEvents.insert(fnex::EventId(19234, 55, 601, 0, 38));
182  // fMissingEvents.insert(fnex::EventId(19270, 15, 392, 0, 38));
183  // fMissingEvents.insert(fnex::EventId(19302, 9, 411, 0, 37));
184  // fMissingEvents.insert(fnex::EventId(19351, 25, 609, 0, 37));
185  // fMissingEvents.insert(fnex::EventId(19542, 13, 13, 0, 30));
186  // fMissingEvents.insert(fnex::EventId(19610, 47, 898, 0, 35));
187  // fMissingEvents.insert(fnex::EventId(19626, 6, 265, 0, 37));
188  // fMissingEvents.insert(fnex::EventId(19708, 37, 312, 0, 30));
189  // fMissingEvents.insert(fnex::EventId(19712, 25, 937, 0, 43));
190  // fMissingEvents.insert(fnex::EventId(19714, 11, 564, 0, 47));
191  // fMissingEvents.insert(fnex::EventId(19732, 27, 153, 0, 36));
192  // fMissingEvents.insert(fnex::EventId(20945, 40, 821, 0, 32));
193  // fMissingEvents.insert(fnex::EventId(21249, 24, 655, 0, 35));
194  // fMissingEvents.insert(fnex::EventId(21253, 31, 933, 0, 23));
195  // fMissingEvents.insert(fnex::EventId(21267, 43, 687, 0, 27));
196  // fMissingEvents.insert(fnex::EventId(21326, 56, 458, 0, 33));
197  // fMissingEvents.insert(fnex::EventId(21330, 16, 490, 0, 27));
198  // fMissingEvents.insert(fnex::EventId(22064, 42, 266, 0, 31));
199  // fMissingEvents.insert(fnex::EventId(22065, 55, 997, 0, 37));
200  // fMissingEvents.insert(fnex::EventId(22251, 54, 694, 0, 38));
201  // fMissingEvents.insert(fnex::EventId(22287, 35, 518, 0, 34));
202  // fMissingEvents.insert(fnex::EventId(22398, 47, 50, 0, 39));
203  // fMissingEvents.insert(fnex::EventId(22425, 59, 97, 0, 26));
204  // fMissingEvents.insert(fnex::EventId(22427, 12, 863, 0, 30));
205  // fMissingEvents.insert(fnex::EventId(22496, 34, 546, 0, 22));
206  // fMissingEvents.insert(fnex::EventId(22523, 34, 208, 0, 33));
207  // fMissingEvents.insert(fnex::EventId(22628, 2, 755, 0, 25));
208  // fMissingEvents.insert(fnex::EventId(22811, 3, 827, 0, 34));
209  // fMissingEvents.insert(fnex::EventId(22999, 58, 545, 0, 29));
210  // fMissingEvents.insert(fnex::EventId(23029, 57, 726, 0, 34));
211  // fMissingEvents.insert(fnex::EventId(23036, 25, 968, 0, 27));
212  // fMissingEvents.insert(fnex::EventId(23087, 10, 109, 0, 27));
213  // fMissingEvents.insert(fnex::EventId(23120, 22, 16, 0, 31));
214  // fMissingEvents.insert(fnex::EventId(23194, 40, 319, 0, 32));
215  // fMissingEvents.insert(fnex::EventId(23325, 22, 479, 0, 40));
216  // fMissingEvents.insert(fnex::EventId(23364, 26, 871, 0, 34));
217  // fMissingEvents.insert(fnex::EventId(24624, 19, 553, 0, 34));
218  // fMissingEvents.insert(fnex::EventId(24644, 39, 954, 0, 32));
219  // fMissingEvents.insert(fnex::EventId(24687, 10, 92, 0, 31));
220  // fMissingEvents.insert(fnex::EventId(24800, 27, 968, 0, 39));
221  // fMissingEvents.insert(fnex::EventId(24843, 63, 929, 0, 46));
222  // fMissingEvents.insert(fnex::EventId(25084, 47, 58, 0, 34));
223  // fMissingEvents.insert(fnex::EventId(25157, 10, 697, 0, 41));
224  // fMissingEvents.insert(fnex::EventId(25252, 60, 501, 0, 22));
225  // fMissingEvents.insert(fnex::EventId(25254, 61, 805, 0, 38));
226  // fMissingEvents.insert(fnex::EventId(25294, 27, 233, 0, 38));
227  // fMissingEvents.insert(fnex::EventId(25344, 8, 701, 0, 28));
228  // fMissingEvents.insert(fnex::EventId(25351, 55, 456, 0, 42));
229 
230  // fMissingEvents.insert(fnex::EventId(19195, 10, 73, 0, 37));
231  // fMissingEvents.insert(fnex::EventId(19546, 12, 222, 0, 30));
232  // fMissingEvents.insert(fnex::EventId(19580, 19, 584, 0, 36));
233  // fMissingEvents.insert(fnex::EventId(19662, 58, 54, 0, 27));
234  // fMissingEvents.insert(fnex::EventId(21650, 54, 267, 0, 28));
235  // fMissingEvents.insert(fnex::EventId(22797, 43, 678, 0, 22));
236  // fMissingEvents.insert(fnex::EventId(23330, 24, 252, 0, 33));
237 
238  }
239 
240  //......................................................................
242  {
243  this->clear();
244  }
245 
246  //......................................................................
247  // Method to clear out the collections of data products after the
248  // writeResults method is called.
250  {
251  fMDToEvents .clear();
252  fNumuFHCQuantCuts.clear();
253  fNumuRHCQuantCuts.clear();
254  fSelections .clear();
255  }
256 
257  //......................................................................
259  {
260  fCAFFileNames = p.get<std::vector<std::string> >("CAFFileNames" );
261 
262  // assume we don't mix FHC and RHC files in the same job, and the
263  // fIsFHC should only be used for cosmics
264  if(fCAFFileNames.begin()->find("rhc") != std::string::npos)
265  fIsFHC = false;
266 
267  fCAFRecTreeName = p.get<std::string >("CAFRecTreeName", "recTree");
268  fSpillTreeName = p.get<std::string >("SpillTreeName", "spillTree");
269  fPOTHistName = p.get<std::string >("CAFPOTHistName", "TotalPOT");
270  fNumQuantiles = p.get<unsigned int>("NumQuantiles", 4);
271  fMakeSelectedList = p.get<bool >("MakeSelectedList", false);
272  fFHCQuantFileName = p.get<std::string >("FHCQuantFileName", "ana2018/Quantiles/quantiles__fhc__full__numu2018.root");
273  fRHCQuantFileName = p.get<std::string >("FHCQuantFileName", "ana2018/Quantiles/quantiles__rhc__full__numu2018.root");
274 
275  auto detName = p.get<std::string>("DetectorName", "Near");
276  fDetector = (detName == "Near") ? novadaq::cnv::kNEARDET : novadaq::cnv::kFARDET;
277 
278  auto selType = p.get<std::string>("Selections");
279 
280  if(selType == "NuESel"){
284  }
285  else if(selType == "NuMuSel"){
290 
291  this->DefineNumuQuantiles();
292  }
293 
294  auto ftString = p.get<std::string>("FileType");
295 
296  if (ftString == "FluxSwap")
297  fFileType = kSwap;
298  else if(ftString == "NonSwap")
299  fFileType = kBeam;
300  else if(ftString == "TauSwap")
302  else if(ftString == "CosmicBackground")
304  else{
306  fIsRealData = true;
307  } // end tests on ftString
308  }
309 
310  //......................................................................
312  {
313  // set up the numu quantile cuts
314 
315  std::string filePath;
316 
317  cet::search_path sp("FW_SEARCH_PATH");
318  if( !sp.find_file(fFHCQuantFileName, filePath) )
319  throw cet::exception("CAFToEventList")
320  << "Cannot find FHC numu quantile file "
322 
323  TFile fhcTF(filePath.c_str(), "READ");
324  TH2F *FDSpec2D = dynamic_cast<TH2F*>(fhcTF.Get("FDSpec2D"));
325  if(FDSpec2D == nullptr)
326  throw cet::exception("CAFToEventList")
327  << "Cannot find histogram FDSpec2D in "
328  << filePath;
329 
331 
332  if( !sp.find_file(fRHCQuantFileName, filePath) )
333  throw cet::exception("EventListMaker")
334  << "Cannot find RHC numu quantile file "
336 
337  TFile rhcTF(filePath.c_str(), "READ");
338  FDSpec2D = dynamic_cast<TH2F*>(rhcTF.Get("FDSpec2D"));
339  if(FDSpec2D == nullptr)
340  throw cet::exception("CAFToEventList")
341  << "Cannot find histogram FDSpec2D in "
342  << filePath;
343 
345  }
346 
347  //......................................................................
348  // see https://cdcvs.fnal.gov/redmine/projects/novaart/wiki/Period_and_Epoch_Naming
350  {
351  std::string epochStr("Unknown");
352 
354  if (run < 10408) epochStr = "1";
355  else if(run < 10974) epochStr = "2";
356  else if(run < 11629) epochStr = "3";
357  else if(run < 11926) epochStr = "4";
358  else if(run < 12087) epochStr = "5";
359  //else if(run < 12290) epochStr = "6a";
360  //else if(run < 12517) epochStr = "6b";
361  //else if(run < 12580) epochStr = "7a";
362  //else epochStr = "7b";
363  // In the 2018 analysis there is no MC for runs after epoch 6a
364  // so make that the last epoch string we worry about too
365  else epochStr = "6";
366  }
367  else if(fDetector == novadaq::cnv::kFARDET){
368  if (run < 17140) epochStr = "1";
369  else if(run < 19747) epochStr = "2";
370  else if(run < 23420) epochStr = "3";
371  else if(run < 24614) epochStr = "4";
372  else if(run < 25413) epochStr = "5";
373  //else if(run < 26686) epochStr = "6a";
374  //else if(run < 28037) epochStr = "6b";
375  //else if(run < 28550) epochStr = "7a";
376  //else epochStr = "7b";
377  // In the 2018 analysis there is no MC for runs after epoch 6a
378  // so make that the last epoch string we worry about too
379  else epochStr = "6";
380  }
381 
382  return epochStr;
383  }
384 
385  //......................................................................
387  {
388  // here we will create an event list for each possible meta data
389  // for the current epoch and data/MC type.
390 
392 
393  std::set<fnex::InteractionType_t> intTypes({fnex::kNuMuCC,
395  fnex::kNuECC,
399  fnex::kNC});
400 
401  // TODO: make sure to update the epoch list - better yet would be to configure this somehow
402  std::set<std::string> epochs({"1", "2", "3", "4", "5", "6"});
403 
404  // set the bits of the metadata that won't change
405  bool isMC = (fFileType != fnex::kDataFile);
406 
407  for(auto const& e : epochs){
408 
409  // loop over the selection and interaction types to flesh out the lists
410  for(auto const& s : fSelections){
411 
412  // there are no nue peripheral events in the ND
414  s == fnex::kNuESelectionPeripheral) continue;
415 
416  if(fFileType == fnex::kDataFile){
417 
418  // data have an unknown interaction type
419  fnex::MetaData md(isMC,
420  fDetector,
421  fFileType,
422  s,
424  e);
425 
426  // make a new EventList for the data
427  fMDToEvents.emplace(md, fnex::EventList(md, ss));
428  } // end if data
430 
431  fnex::MetaData md(true,
432  fDetector,
433  fFileType,
434  s,
436  e);
437 
438  // make a new EventList for the cosmic background
439  fMDToEvents.emplace(md, fnex::EventList(md, ss));
440  } // end if cosmic background
441  else{
442  // loop over possible interaction types for the metadata
443  // don't use cosmic muons, those are special and are handled
444  // in another way - we don't want to inflate the POT for those lists
445  for(auto const& i : intTypes){
446 
447  fnex::MetaData md(isMC,
448  fDetector,
449  fFileType,
450  s,
451  i,
452  e);
453 
454  // make a new EventList for the MC
455  fMDToEvents.emplace(md, fnex::EventList(md, ss));
456  } // end loop over interaction types
457  } // end if not data or MC, ie cosmic background
458 
459  } // end loop over selection types
460  } // end loop over epochs
461 
462  }
463 
464  //......................................................................
465  // since we are using CAF files, which are not art files, we don't need the
466  // analyze function to do anything
468  {
469  }
470 
471  //......................................................................
473  {
474  if(!this->PassesSpillAndBeamCuts(proxy)) return false;
475 
477  if(this->PassesNumuSelections(proxy)) return true;
478  }
479  else if(fSelections.find(fnex::kNuESelectionLowPID) != fSelections.end()){
480  if(this->PassesNueSelections(proxy)) return true;
481  }
482 
483  return false;
484  }
485 
486  //......................................................................
488  {
490  if(!ana::kNue2018FDAllSamples(&proxy)){
491  LOG_DEBUG("CAFToEventList")
492  << "Fails nue FD All samples";
493 
494  return false;
495  }
496  }
497  else if(fDetector == novadaq::cnv::kNEARDET){
498  if(!ana::kNue2018NDEnergy(&proxy) ||
499  !ana::kNue2018NDPresel(&proxy) ||
500  !ana::kNue2018NDCVNSsb(&proxy) ){
501  LOG_DEBUG("CAFToEventList")
502  << "Fails nue ND energy, presel and CVNSsb "
503  << ana::kNue2018NDEnergy(&proxy)
504  << " "
505  << ana::kNue2018NDPresel(&proxy)
506  << " "
507  << ana::kNue2018NDCVNSsb(&proxy);
508  return false;
509  }
510  }
511 
512  return true;
513  }
514 
515  //......................................................................
517  {
518  // check that this event passes the quality cuts and PID selections
519  // for numu
521  if(!ana::kNumuCutFD2018(&proxy)){
522 
523 #ifdef PROD4
524 
525  LOG_DEBUG("CAFToEventList")
526  << "failed numu FD "
527  << ana::kNumuQuality(&proxy)
528  << " contain "
529  << ana::kNumuContainFD2017(&proxy)
530  << " "
532  << " "
534  << " PID 2018 "
535  << ana::kNumuPID2018(&proxy)
536  << " "
537  << proxy.sel.remid.pid
538  << " "
539  << proxy.sel.cvnProd3Train.numuid
540  << " "
541  << proxy.sel.cvn2017.numuid
542  << " cosRej: "
543  << ana::kNumuCosmicRej2018(&proxy)
544  << " "
545  << proxy.sel.cosrej.anglekal
546  << " "
547  << proxy.sel.cosrej.numucontpid
548  << " "
549  << proxy.slc.nhit
550  << " "
551  << proxy.sel.nuecosrej.pngptp;
552 #endif
553  return false;
554  }
555  }
556  else if(fDetector == novadaq::cnv::kNEARDET){
557  if(!ana::kNumuCutND2018(&proxy)){
558  LOG_DEBUG("CAFToEventList")
559  << "failed numu ND "
560  << ana::kNumuQuality(&proxy)
561  << " "
562  << ana::kNumuContainND2017(&proxy)
563  << " "
564  << ana::kNumuPID2018(&proxy);
565  return false;
566  }
567  }
568 
569  return true;
570  }
571 
572  //......................................................................
574  {
575  if(this->PassesNumuSelections(srProxy)){
576  if(srProxy.spill.isFHC){
577  if(fNumuFHCQuantCuts[0](&srProxy)) return fnex::kNuMuSelectionQ1;
578  else if(fNumuFHCQuantCuts[1](&srProxy)) return fnex::kNuMuSelectionQ2;
579  else if(fNumuFHCQuantCuts[2](&srProxy)) return fnex::kNuMuSelectionQ3;
580  else if(fNumuFHCQuantCuts[3](&srProxy)) return fnex::kNuMuSelectionQ4;
581  } // end if fhc
582  else{
583  if(fNumuRHCQuantCuts[0](&srProxy)) return fnex::kNuMuSelectionQ1;
584  else if(fNumuRHCQuantCuts[1](&srProxy)) return fnex::kNuMuSelectionQ2;
585  else if(fNumuRHCQuantCuts[2](&srProxy)) return fnex::kNuMuSelectionQ3;
586  else if(fNumuRHCQuantCuts[3](&srProxy)) return fnex::kNuMuSelectionQ4;
587  }
588  } // end if numu
589  else if(this->PassesNueSelections(srProxy) && ana::kNue2018CVNCut(&srProxy)){
591  else if(ana::kCVNSSe(&srProxy) < ana::kNue2018CVNFHCHighEdge &&
592  srProxy.spill.isFHC) return fnex::kNuESelectionHighPID;
593  else if(ana::kCVNSSe(&srProxy) < ana::kNue2018CVNRHCHighEdge &&
594  !srProxy.spill.isFHC) return fnex::kNuESelectionHighPID;
595  else return fnex::kNuESelectionLowPID;
596 
597  }
598 
600  }
601 
602  //......................................................................
603  // note that srProxy is not const here because some of the CAF/Var functions
604  // do not take const pointers to it...seems silly to not have it const.
606  caf::SRProxy & srProxy) const
607  {
608 
609  mcVarVals.set_val_at (fnex::kTrueE, (float)ana::kTrueE(&srProxy) );
610  mcVarVals.set_val_at (fnex::kTruePDG, (float)ana::kTruePDG(&srProxy) );
611  mcVarVals.set_val_at (fnex::kTrueCCNC, (srProxy.mc.nu[0].iscc) ? simb::kCC : simb::kNC);
612  mcVarVals.set_val_at (fnex::kTrueIntType, srProxy.mc.nu[0].inttype);
613  mcVarVals.set_val_at (fnex::kTrueIntMode, (float)ana::kMode(&srProxy) );
614  mcVarVals.set_val_at (fnex::kTrueHitNuc, (float)ana::kHitNuc(&srProxy) );
615  mcVarVals.set_val_at (fnex::kTruePDGOrig, srProxy.mc.nu[0].pdgorig );
616  mcVarVals.set_val_at (fnex::kTrueParentPDG, srProxy.mc.nu[0].beam.ptype );
617  mcVarVals.set_val_at (fnex::kTrueParentDecay, srProxy.mc.nu[0].beam.ndecay );
618  mcVarVals.set_val_at (fnex::kTrueParentTargetPDG, srProxy.mc.nu[0].beam.tptype );
619  mcVarVals.set_val_at (fnex::kTrueParentPZ, (float)ana::kTruetpz(&srProxy) );
620  mcVarVals.set_val_at (fnex::kTrueParentPT, (float)ana::kTruetpT(&srProxy));
621  mcVarVals.set_val_at (fnex::kXSecCVPPFX_Weight, (float)(ana::kXSecCVWgt2018(&srProxy) * ana::kPPFXFluxCVWgt(&srProxy)));
622 
623  mcVarVals.set_val_at (fnex::kRPACCQE_Weight, (float)ana::kRPAWeightCCQE2018(&srProxy));
624  mcVarVals.set_val_at (fnex::kRPARES_Weight, (float)ana::kRPAWeightRES2017(&srProxy));
625  mcVarVals.set_val_at (fnex::kEmpiricalMEC_Weight, (float)ana::kEmpiricalMECWgt2018(&srProxy) );
628 
629  double wgtPlus = 1.;
630  double wgtMinus = 1.;
631  ana::kRPACCQEEnhSyst2018.Shift( 1., &srProxy, wgtPlus);
632  ana::kRPACCQEEnhSyst2018.Shift(-1., &srProxy, wgtMinus);
633  mcVarVals.set_GENIEWgt(fnex::kRPACCQEshapeEnh, "-1sigma", (float)wgtMinus);
634  mcVarVals.set_GENIEWgt(fnex::kRPACCQEshapeEnh, "+1sigma", (float)wgtPlus);
635 
636  ana::kRPACCQESuppSyst2018.Shift( 1., &srProxy, wgtPlus);
637  ana::kRPACCQESuppSyst2018.Shift(-1., &srProxy, wgtMinus);
638  mcVarVals.set_GENIEWgt(fnex::kRPACCQEshapeSupp, "-1sigma", (float)wgtMinus);
639  mcVarVals.set_GENIEWgt(fnex::kRPACCQEshapeSupp, "+1sigma", (float)wgtPlus);
640 
645 
646 
647  // Now for the full set of possible GENIE knobs
648 
649  static std::map<rwgt::EReweightLabel, std::string> labelMap
650  {
651  // NC elastic parameters
652  {rwgt::fReweightMaNCEL, "MaNCEL"},
653  {rwgt::fReweightEtaNCEL, "EtaNCEL"},
654 
655  // CCQE parameters
656  {rwgt::fReweightNormCCQE, "NormCCQE"},
657  {rwgt::fReweightNormCCQEenu, "NormCCQE_EnuDep"},
658  {rwgt::fReweightMaCCQEshape, "MaCCQE_shape"},
659  {rwgt::fReweightMaCCQE, "MaCCQE"},
660  {rwgt::fReweightVecCCQEshape, "VecCCQE_shape"},
661 
662  // Resonant production: CC
663  {rwgt::fReweightNormCCRES, "NormCCRES"},
664  {rwgt::fReweightMaCCRESshape, "MaCCRES_shape"},
665  {rwgt::fReweightMvCCRESshape, "MvCCRES_shape"},
666  {rwgt::fReweightMaCCRES, "MaCCRES"},
667  {rwgt::fReweightMvCCRES, "MvCCRES"},
668 
669  // Resonant production: NC
670  {rwgt::fReweightNormNCRES, "NormNCRES"},
671  {rwgt::fReweightMaNCRESshape, "MaNCRES_shape"},
672  {rwgt::fReweightMvNCRESshape, "MvNCRES_shape"},
673  {rwgt::fReweightMaNCRES, "MaNCRES"},
674  {rwgt::fReweightMvNCRES, "MvNCRES"},
675 
676  // Coherent pion production
677  {rwgt::fReweightMaCOHpi, "MaCOHpi"},
678  {rwgt::fReweightR0COHpi, "R0COHpi"},
679 
680  // Non-resonant background
681  {rwgt::fReweightRvpCC1pi, "RvpCC1pi"},
682  {rwgt::fReweightRvpCC2pi, "RvpCC2pi"},
683  {rwgt::fReweightRvpNC1pi, "RvpNC1pi"},
684  {rwgt::fReweightRvpNC2pi, "RvpNC2pi"},
685  {rwgt::fReweightRvnCC1pi, "RvnCC1pi"},
686  {rwgt::fReweightRvnCC2pi, "RvnCC2pi"},
687  {rwgt::fReweightRvnNC1pi, "RvnNC1pi"},
688  {rwgt::fReweightRvnNC2pi, "RvnNC2pi"},
689  {rwgt::fReweightRvbarpCC1pi, "RvbarpCC1pi"},
690  {rwgt::fReweightRvbarpCC2pi, "RvbarpCC2pi"},
691  {rwgt::fReweightRvbarpNC1pi, "RvbarpNC1pi"},
692  {rwgt::fReweightRvbarpNC2pi, "RvbarpNC2pi"},
693  {rwgt::fReweightRvbarnCC1pi, "RvbarnCC1pi"},
694  {rwgt::fReweightRvbarnCC2pi, "RvbarnCC2pi"},
695  {rwgt::fReweightRvbarnNC1pi, "RvbarnNC1pi"},
696  {rwgt::fReweightRvbarnNC2pi, "RvbarnNC2pi"},
697 
698  // DIS tweaking parameters
699  {rwgt::fReweightAhtBY, "AhtBY"},
700  {rwgt::fReweightBhtBY, "BhtBY"},
701  {rwgt::fReweightCV1uBY, "CV1uBY"},
702  {rwgt::fReweightCV2uBY, "CV2uBY"},
703  {rwgt::fReweightAhtBYshape, "AhtBYshape"},
704  {rwgt::fReweightBhtBYshape, "BhtBYshape"},
705  {rwgt::fReweightCV1uBYshape, "CV1uBYshape"},
706  {rwgt::fReweightCV2uBYshape, "CV2uBYshape"},
707  {rwgt::fReweightNormDISCC, "NormDISCC"},
708  {rwgt::fReweightRnubarnuCC, "RnubarnuCC"},
709  {rwgt::fReweightDISNuclMod, "DISNuclMod"},
710 
711  {rwgt::fReweightNC, "NC"},
712 
713  // Hadronization (free nucleon target)
714  {rwgt::fReweightAGKY_xF1pi, "AGKY_xF1pi"},
715  {rwgt::fReweightAGKY_pT1pi, "AGKY_pT1pi"},
716 
717  // Medium-effects to hadronization
718  {rwgt::fReweightFormZone, "FormZone"},
719 
720  // Intranuclear rescattering parameters.
721  {rwgt::fReweightMFP_pi, "MFP_pi"},
722  {rwgt::fReweightMFP_N, "MFP_N"},
723  {rwgt::fReweightFrCEx_pi, "FrCEx_pi"},
724  {rwgt::fReweightFrInel_pi, "FrInel_pi"},
725  {rwgt::fReweightFrAbs_pi, "FrAbs_pi"},
726  {rwgt::fReweightFrPiProd_pi, "FrPiProd_pi"},
727  {rwgt::fReweightFrCEx_N, "FrCEx_N"},
728  // the next two appear no longer to be in use with GENIE v3
729  //{rwgt::fReweightFrElas_pi, "FrElas_pi"},
730  //{rwgt::fReweightFrElas_N, "FrElas_N"},
731  {rwgt::fReweightFrInel_N, "FrInel_N"},
732  {rwgt::fReweightFrAbs_N, "FrAbs_N"},
733  {rwgt::fReweightFrPiProd_N, "FrPiProd_N"},
734 
735  // Nuclear model
736  {rwgt::fReweightCCQEPauliSupViaKF, "CCQEPauliSupViaKF"},
737  {rwgt::fReweightCCQEMomDistroFGtoSF, "CCQEMomDistroFGtoSF"},
738 
739  // Resonance decays
740  {rwgt::fReweightBR1gamma, "BR1gamma"},
741  {rwgt::fReweightBR1eta, "BR1eta"},
742  {rwgt::fReweightTheta_Delta2Npi, "Theta_Delta2Npi"}
743  };
744 
745  std::string rwgtLabel;
746  for (auto const& rwgtPair : labelMap){
747 
748  auto wgtFunc = ana::GetGenieRwgtSyst(rwgtPair.first);
749  rwgtLabel = rwgtPair.second;
750 
751  wgtFunc->Shift(-2., &srProxy, wgtMinus);
752  mcVarVals.set_GENIEWgt(rwgtLabel, "-2sigma", (float)wgtMinus);
753 
754  wgtFunc->Shift(-1., &srProxy, wgtMinus);
755  mcVarVals.set_GENIEWgt(rwgtLabel, "-1sigma", (float)wgtMinus);
756 
757  wgtFunc->Shift( 1., &srProxy, wgtPlus);
758  mcVarVals.set_GENIEWgt(rwgtLabel, "+1sigma", (float)wgtPlus);
759 
760 
761  wgtFunc->Shift( 2., &srProxy, wgtPlus);
762  mcVarVals.set_GENIEWgt(rwgtLabel, "+2sigma", (float)wgtPlus);
763  } // end loop over genie reweight variables
764  }
765 
766  //......................................................................
768  caf::SRProxy const& srProxy) const
769  {
770 
771  if(this->PassesNueSelections(srProxy)){
772  dataVals.set_val_at(fnex::kLep_RecoE, (float)ana::kRecoEME (&srProxy));
773  dataVals.set_val_at(fnex::kHad_RecoE, (float)ana::kRecoHADE(&srProxy));
774  dataVals.set_val_at(fnex::kNuE_CVN, (float)ana::kCVNSSe (&srProxy));
775  dataVals.set_val_at(fnex::kNuE_NumMichel, (float)ana::kNMichels(&srProxy));
776  }
777  else if(this->PassesNumuSelections(srProxy)){
778  dataVals.set_val_at(fnex::kLep_RecoE, (float)(ana::kCCE(&srProxy) - ana::kHadE(&srProxy)));
779  dataVals.set_val_at(fnex::kHad_RecoE, (float)ana::kHadE(&srProxy) );
780  dataVals.set_val_at(fnex::kLep_RecoE_MCFrac, 1. );
781  dataVals.set_val_at(fnex::kRecoQ2, (float)ana::kRecoQ2(&srProxy));
782  }
783 
784  }
785 
786  //......................................................................
788  caf::SRProxy const& srProxy) const
789  {
790  evId.run = srProxy.hdr.run;
791  evId.subrun = srProxy.hdr.subrun;
792  evId.event = srProxy.hdr.evt;
793  evId.slice = srProxy.hdr.subevt;
794  evId.cycle = srProxy.hdr.cycle;
795 
796  }
797 
798  //......................................................................
800  {
801  return (ana::kStandardSpillCuts(&(proxy.spill)) && ana::kInBeamSpill(&proxy));
802  }
803 
804  //......................................................................
805  int CAFToEventList::PeriodAdjustCosmics(int const& inputPeriod)
806  {
807  int period = inputPeriod;
808 
809  // set the period value to be 5000 for the fhc tune with
810  // period === 4000 and set the period value to 4000 with
811  // the rhc tune and period != 4000. Later we set the
812  // period == 6000 exposure to be the same as the period 4000
813  if (period == 4000 && fIsFHC) period = 5000;
814  else if(period != 4000 && !fIsFHC) period = 4000;
815 
816  return period;
817  }
818 
819  //......................................................................
821  {
822  auto *spillTree = dynamic_cast<TTree*>(tf.Get(fSpillTreeName.c_str()));
823 
824  caf::SRSpill *sp = nullptr;
825 
826  spillTree->SetBranchAddress("spill", &sp);
827  spillTree->GetEntry(0);
828 
829  static caf::SRSpillProxy spProxy(nullptr, "");
830 
831  LOG_VERBATIM("CAFToEventList")
832  << "There are "
833  << spillTree->GetEntriesFast()
834  << " spills";
835 
836  int period;
837  float totalPOT = 0.;
838 
839  // loop over the spill tree to add up the POT for each period
840  //for(int i = 0; i < 10; ++i){
841  for(int i = 0; i < spillTree->GetEntriesFast(); ++i){
842 
843  spillTree->GetEntry(i);
844 
845  // get a proxy for each entry
846  spProxy = *sp;
847 
848  // ignore spills that don't pass the standard cuts
849  if(!ana::kStandardSpillCuts(&spProxy)) continue;
850 
851  // now figure out which period this spill belongs to
852  // periods always have a factor of 1000 applied to
853  // them in case we want epoch by epoch
854  period = 1000 * std::atoi((this->RunToEpoch(sp->run)).c_str());
855 
856  // TODO: remove this hack for post 2018 analyses
857  if(fFileType == fnex::kCosmicBackgroundFile) period = this->PeriodAdjustCosmics(period);
858 
859  if(std::isnan(sp->spillpot) ||
860  std::isnan(sp->livetime) ){
861  LOG_DEBUG("CAFToEventList")
862  << tf.GetName()
863  << " "
864  << sp->run
865  << " "
866  << sp->subrun
867  << " "
868  << sp->evt
869  << " spill has nan for pot or livetime "
870  << sp->spillpot
871  << " "
872  << sp->livetime
873  << " period "
874  << period;
875 
876  continue;
877  }
878 
879  // spillpot is in units of 10^20 POT
881  fPeriodToPOTLiveTime[period].liveTime += sp->livetime;
882 
883  totalPOT += spProxy.spillpot;
884 
885  } // end loop over the entries in the spill tree
886 
887  // if the spill tree and histogram totals are off by more than 1%, say something
888  if(std::abs(1. - totalPOT / dynamic_cast<TH1F* >(tf.Get(fPOTHistName.c_str()))->Integral()) > 0.01)
889  LOG_VERBATIM("CAFToEventList")
890  << "pot hist has "
891  << dynamic_cast<TH1F* >(tf.Get(fPOTHistName.c_str()))->Integral()
892  << " POT compared to "
893  << totalPOT
894  << " from the tree";
895 
896  // print out the number of pot for all files so far
897  for(auto const& itr : fPeriodToPOTLiveTime){
898  if(itr.second.pot > 0. ||
899  itr.second.liveTime > 0.)
900  LOG_VERBATIM("CAFToEventList")
901  << "period "
902  << itr.first
903  << " has "
904  << itr.second.pot * 1.e-20
905  << " x 10^20 POT and "
906  << itr.second.liveTime
907  << " s live time";
908  }
909 
910  }
911 
912  //......................................................................
913  void CAFToEventList::FillVariables(std::vector<fnex::DataVarVals> & dataVals,
914  std::vector<fnex::MCVarVals> & mcVals,
915  std::vector<fnex::EventId> & evIds,
916  std::vector<fnex::SelectionType_t> & sels)
917  {
918  dataVals.clear();
919  mcVals .clear();
920  evIds .clear();
921 
922  for(auto const& fileName : fCAFFileNames){
923 
924  // open the CAF file and get the TTrees and histograms for reading
925  TFile *tf = TFile::Open(fileName.c_str());
926 
927  this->FindPOTPerPeriod(*tf);
928 
929  auto *cafTree = dynamic_cast<TTree*>(tf->Get(fCAFRecTreeName.c_str()));
930 
931  // hook the TTree up to a standard record and get the first entry to initialize things
932  caf::StandardRecord *sr = nullptr;
933 
934  cafTree->SetBranchAddress("rec", &sr);
935  cafTree->GetEntry(0);
936 
937  static caf::SRProxy srProxy(nullptr, "");
938 
939  LOG_DEBUG("CAFToEventList")
940  << "There are "
941  << cafTree->GetEntriesFast()
942  << " caf entries";
943 
944  // now loop over the CAF entries and fill the event list variables
945  //for(int i = 0; i < 10; ++i){
946  for(int i = 0; i < cafTree->GetEntriesFast(); ++i){
947 
948  cafTree->GetEntry(i);
949 
950  // get a proxy for each entry
951  srProxy = *sr;
952 
953 
954  // Uncomment the following block to debug cases where events
955  // show up in the official CAF selection but not this one, or
956  // vice versa
957  /*
958  fnex::EventId curEvent(srProxy.hdr.run,
959  srProxy.hdr.subrun,
960  srProxy.hdr.evt,
961  srProxy.hdr.cycle,
962  srProxy.hdr.subevt);
963 
964  if( fMissingEvents.find(curEvent) != fMissingEvents.end())
965  LOG_VERBATIM("MissingEventsCheck")
966  << curEvent
967  << " beam "
968  << ana::kInBeamSpill(&srProxy)
969  << " spill "
970  << ana::kStandardSpillCuts(&(srProxy.spill))
971  << " FD all "
972  << ana::kNue2018FDAllSamples(&srProxy);
973  */
974  if( !this->PassesSelections(srProxy) ) continue;
975 
976  fnex::MCVarVals mcVarVals;
977  fnex::DataVarVals dataVarVals;
978  fnex::EventId evId;
979 
980  // fill the event id information
981  this->FillEventId(evId, srProxy);
982 
983  // fill MC information if this is mc
984  if(srProxy.hdr.ismc){
985  if(srProxy.mc.nu.size() == 1) this->FillMCVarVals(mcVarVals, srProxy);
986  else{
987  LOG_VERBATIM("CAFToEventList")
988  << "There are "
989  << srProxy.mc.nu.size()
990  << " neutrinos in this event "
991  << evId;
992 
993  continue;
994  }
995  } // end if the event is MC
996 
997  // fill the data values
998  this->FillDataVarVals(dataVarVals, srProxy);
999 
1000  // add everything to the vectors
1001  dataVals.push_back(dataVarVals);
1002  mcVals .push_back(mcVarVals );
1003  evIds .push_back(evId );
1004  sels .push_back(this->FindSelectionType(srProxy));
1005 
1006  } // end loop over the CAF StandardRecords
1007 
1008  tf->Close();
1009  } // end loop over CAF file names
1010  }
1011 
1012  //......................................................................
1013  // this method is a hack to account for the fact that there were no
1014  // period 6 cosmic ray background files for the 2018 analysis
1016  {
1017  if(fIsFHC) return;
1018 
1019  // loop over the fMDToEvents map and for all period 4 cosmic metadata
1020  // add a corresponding period 6 cosmic metadata
1021  for(auto const& itr : fMDToEvents){
1022 
1023  if(itr.first.Period() == 4000 &&
1024  itr.first.fileType == fnex::kCosmicBackgroundFile){
1025 
1026  fnex::MetaData md = itr.first;
1027  md.epoch = 6000;
1028 
1029  fMDToEvents.find(md)->second.AddSpillSummary(itr.second.ListSpillSummary());
1030 
1031  // loop over the events in the list and add them to the new list
1032  for(auto const& ev : itr.second)
1033  fMDToEvents.find(md)->second.AddEvent(ev);
1034 
1035  }// end if we are hacking cosmic backgrounds for rhc
1036 
1037  } // end loop over metadata and event lists
1038 
1039  }
1040 
1041  //......................................................................
1043  {
1044  this->InitializeEventListMaps();
1045 
1048  float pdg = 0;
1049 
1050  std::vector<fnex::DataVarVals > dataVals;
1051  std::vector<fnex::MCVarVals > mcVals;
1052  std::vector<fnex::EventId > evIds;
1053  std::vector<fnex::SelectionType_t> sels;
1054 
1055  this->FillVariables(dataVals, mcVals, evIds, sels);
1056 
1057  LOG_VERBATIM("CAFToEventList")
1058  << "there are "
1059  << dataVals.size()
1060  << " data vals "
1061  << mcVals.size()
1062  << " mc vals "
1063  << evIds.size()
1064  << " event ids "
1065  << sels.size()
1066  << " selections ";
1067 
1068  auto periodToPOTItr = fPeriodToPOTLiveTime.begin();
1069 
1070  // loop over the event lists and add SpillSummary info
1071  for(auto & itr : fMDToEvents){
1072 
1074 
1075  periodToPOTItr = fPeriodToPOTLiveTime.find(itr.first.Period());
1076 
1077  // multiply the weight for this period by the total POT
1078  // indicated in the list of CAF files used. We only
1079  // look at either FHC or RHC SAM datalists separately,
1080  // so we don't have to worry about separate accounting
1081  // for those
1082  if(periodToPOTItr != fPeriodToPOTLiveTime.end()){
1083  ss += fnex::SpillSummary(periodToPOTItr->second.pot,
1084  periodToPOTItr->second.pot,
1085  periodToPOTItr->second.liveTime,
1086  1,
1087  1);
1088 
1089  LOG_DEBUG("CAFTToEventList")
1090  << " checking spill summary to add "
1091  << ss
1092  << " "
1093  << periodToPOTItr->second.pot
1094  << " "
1095  << periodToPOTItr->second.liveTime;
1096  }
1097  else{
1098  LOG_WARNING("CAFToEventList")
1099  << itr.first.ToString()
1100  << " could not find period "
1101  << itr.first.Period()
1102  << " in POT weight map, just add zeros";
1103  }
1104 
1105  itr.second.AddSpillSummary(ss);
1106 
1107  LOG_DEBUG("CAFToEventList")
1108  << "adding "
1109  << ss
1110  << " to "
1111  << itr.first.ToString()
1112  << " "
1113  << itr.second.ListSpillSummary();
1114 
1115  } // end loop over the metadata to event lists map
1116 
1117  // now make event lists
1118  for(size_t i = 0; i < dataVals.size(); ++i){
1119 
1120  fnex::EventPtr ev = std::make_unique<fnex::Event>();
1121 
1123 
1124  // figure out the interaction type and then get the correct metadata
1125  // for it
1127  it = fnex::kCosmicMuon;
1128  }
1129  else if(!fIsRealData){
1130  if( mcVals[i].val_at(fnex::kTrueCCNC) == 1.*simb::kNC )
1131  it = fnex::kNC;
1132  else{
1133  pdg = mcVals[i].val_at(fnex::kTruePDG);
1134  if (pdg == 12) it = fnex::kNuECC;
1135  else if(pdg == -12) it = fnex::kNuEBarCC;
1136  else if(pdg == 14) it = fnex::kNuMuCC;
1137  else if(pdg == -14) it = fnex::kNuMuBarCC;
1138  else if(pdg == 16) it = fnex::kNuTauCC;
1139  else if(pdg == -16) it = fnex::kNuTauBarCC;
1140  }
1141  } // end if neutrino MC
1142 
1144  fDetector,
1145  fFileType,
1146  sels[i],
1147  it,
1148  this->RunToEpoch(evIds[i].run));
1149 
1150  ev->eventId = std::make_unique<fnex::EventId >(evIds[i]);
1151  ev->dataVals = std::make_unique<fnex::DataVarVals>(dataVals[i]);
1152 
1154  ev->mcVals = std::make_unique<fnex::MCVarVals>(mcVals[i]);
1155  ev->mcVals->set_val_at(fnex::kTrueE, 0. );
1156  ev->mcVals->set_val_at(fnex::kTruePDG, 13 );
1157  ev->mcVals->set_val_at(fnex::kTruePDGOrig, 13 );
1158  ev->mcVals->set_val_at(fnex::kTrueCCNC, 0. );
1159  ev->mcVals->set_val_at(fnex::kTrueIntType, fnex::kCosmicMuon);
1160  ev->mcVals->set_val_at(fnex::kXSecCVPPFX_Weight, 1. );
1161  ev->mcVals->set_val_at(fnex::kRPACCQE_Weight, 1. );
1162 
1163  } // end if cosmic backgrounds
1164  else if(md.isMC)
1165  ev->mcVals = std::make_unique<fnex::MCVarVals>(mcVals[i]);
1166  else
1167  ev->mcVals.reset(nullptr);
1168 
1169  LOG_DEBUG("CAFToEventList")
1170  << "add the event to "
1171  << md.ToString()
1172  << " "
1173  << fMDToEvents.count(md);
1174 
1175  fMDToEvents.find(md)->second.AddEvent(ev);
1176 
1177  LOG_DEBUG("CAFToEventList")
1178  << i
1179  << " metadata is "
1180  << md.ToString()
1181  << " "
1182  << fMDToEvents.count(md);
1183 
1184  // Handle the nue ND peripheral sample. In the ND as of the 2017 analysis
1185  // the peripheral sample is just the high CVN set so we just copy that
1186  // and write it into a tree
1187  /*
1188  if(md.selectionType == fnex::kNuESelectionHighPID &&
1189  md.detector == novadaq::cnv::kNEARDET){
1190 
1191  // set the selection type to be the peripheral
1192  md.selectionType = fnex::kNuESelectionPeripheral;
1193 
1194  // add the event to the peripheral sample.
1195  fMDToEvents.find(md)->second.AddEvent(ev);
1196  } // end if nue high PID and ND
1197  */
1198 
1199  } // end loop to insert events
1200 
1201  // TODO: remove this hack when not doing the 2018 analysis
1203 
1204  for(auto const & itr : fMDToEvents){
1205 
1206  LOG_VERBATIM("CAFToEventList")
1207  << "metadata "
1208  << itr.first.ToString()
1209  << " has "
1210  << itr.second.ListSpillSummary()
1211  << " and "
1212  << itr.second.size()
1213  << " events";
1214 
1215  if(fMakeSelectedList){
1216  for(auto const& ev : itr.second)
1217  LOG_VERBATIM("SelectedEvents")
1218  << fnex::cSelectionType_Strings[itr.first.selectionType]
1219  << "\t\t"
1220  << ev->eventId->run
1221  << "\t"
1222  << ev->eventId->subrun
1223  << "\t"
1224  << ev->eventId->event
1225  << "\t"
1226  << ev->eventId->slice
1227  << "\t"
1228  << ev->eventId->cycle;
1229  } // end if we are making a list of selected events
1230  }
1231 
1232  LOG_DEBUG("CAFToEventList")
1233  << "Serialize the events";
1234 
1236  fMDToEvents);
1237 
1238  LOG_DEBUG("CAFToEventList")
1239  << "endJob() done";
1240 
1241  }
1242 
1244 
1245 } // end fnex namespace
const Var kTruePDG
Definition: TruthVars.h:48
const Var kHadE
Definition: NumuVars.h:23
caf::Proxy< caf::SRSpill > spill
Definition: SRProxy.h:2010
void analyze(art::Event const &e) override
#define LOG_DEBUG(stream)
Definition: Messenger.h:149
const Var kMode([](const caf::SRProxy *sr){return(sr->mc.nnu==0)?-1:int(sr->mc.nu[0].mode);})
Neutrino interaction mode.
Definition: Vars.h:99
static const unsigned char kNuE_CVN
Definition: VarVals.h:32
const ana::Var kRecoQ2([](const caf::SRProxy *sr){const double M_mu_sqrd=util::sqr(0.1056);double E_mu=kMuE(sr);double p_mu=sqrt(util::sqr(E_mu)-M_mu_sqrd);return 2 *kCCE(sr)*(E_mu-p_mu *kCosNumi(sr))-M_mu_sqrd;})
Reconstructed four-momentum transfer invariant (Q^2)
Definition: NumuVars.h:146
fileName
Definition: plotROC.py:78
std::string fCAFRecTreeName
names of the CAF tree with the reco info
POTLiveTime(double const &p, double const &l)
std::string fSpillTreeName
names of the CAF histogram with the total POT
set< int >::iterator it
enum fnex::file_type FileType_t
const double kNue2018CVNFHCHighEdge
Definition: NueCuts2018.h:41
const NOvARwgtSyst kRPACCQEEnhSyst2018("RPAShapeenh2018","RPA shape: higher-Q^{2} enhancement (2018)", novarwgt::kRPACCQEEnhSyst2018)
Definition: RPASysts.h:13
static const unsigned char kEmpiricalMEC_Weight
Definition: VarVals.h:60
const Var kPPFXFluxCVWgt
Definition: PPFXWeights.h:16
unsigned int fNumQuantiles
number of quantiles
caf::Proxy< unsigned int > subrun
Definition: SRProxy.h:249
static const unsigned char kTrueParentPT
Definition: VarVals.h:48
const Cut kNumuContainND2017([](const caf::SRProxy *sr){if(!sr->vtx.elastic.IsValid) return false;for(unsigned int i=0;i< sr->vtx.elastic.fuzzyk.nshwlid;++i){TVector3 start=sr->vtx.elastic.fuzzyk.png[i].shwlid.start;TVector3 stop=sr->vtx.elastic.fuzzyk.png[i].shwlid.stop;if(std::min(start.X(), stop.X())< -180.0) return false;if(std::max(start.X(), stop.X()) > 180.0) return false;if(std::min(start.Y(), stop.Y())< -180.0) return false;if(std::max(start.Y(), stop.Y()) > 180.0) return false;if(std::min(start.Z(), stop.Z())< 20.0) return false;if(std::max(start.Z(), stop.Z()) > 1525.0) return false;}if(sr->trk.kalman.ntracks< 1) return false;for(unsigned int i=0;i< sr->trk.kalman.ntracks;++i){if(i==sr->trk.kalman.idxremid) continue;else if(sr->trk.kalman.tracks[i].start.Z() > 1275||sr->trk.kalman.tracks[i].stop.Z() > 1275) return false;}return(sr->trk.kalman.ntracks > sr->trk.kalman.idxremid &&sr->slc.firstplane > 1 &&sr->slc.lastplane< 212 &&sr->trk.kalman.tracks[0].start.Z()< 1100 &&(sr->trk.kalman.tracks[0].stop.Z()< 1275 ||sr->sel.contain.kalyposattrans< 55) &&sr->sel.contain.kalfwdcellnd > 5 &&sr->sel.contain.kalbakcellnd > 10);})
Definition: NumuCuts2017.h:11
std::string fFHCQuantFileName
name of file in numudata ups product with FHC quantiles
caf::Proxy< caf::SRHeader > hdr
Definition: SRProxy.h:2004
Base container for the MC related Vars that constitute an event.
Definition: VarVals.h:930
void SerializeEventListMap(std::string const &dirName, fnex::EventListMap const &listMap)
Definition: Event.cxx:176
static const unsigned char kEmpiricalMECtoGENIERES_Weight
Definition: VarVals.h:62
const Cut kNumuCutND2018
Definition: NumuCuts2018.h:41
const char * p
Definition: xmltok.h:285
Proxy for caf::StandardRecord.
Definition: SRProxy.h:1993
void set_val_at(std::string const &varkey, float const &val)
Definition: VarVals.cxx:164
caf::Proxy< std::vector< caf::SRNeutrino > > nu
Definition: SRProxy.h:570
const Cut kNue2018NDCVNSsb
Definition: NueCuts2018.h:153
vector< vector< double > > clear
const HistAxis kHadEFracAxis("E_{had.} / E_{#nu}", Binning::Simple(200, 0, 1), kHadEFrac)
HistAxis that implements Hadronic Energy fraction binning used by L.Vinton to derive Hadronic Energy ...
Definition: HistAxes.h:30
std::vector< ana::Cut > fNumuRHCQuantCuts
quantile boundaries for the numu RHC data
caf::Proxy< float > pid
Definition: SRProxy.h:1074
Float_t ss
Definition: plot.C:24
static const unsigned char kRecoQ2
Definition: VarVals.h:38
static const unsigned char kRPACCQEshapeSupp
Definition: VarVals.h:75
enum fnex::interaction_type InteractionType_t
static const unsigned char kXSecCVPPFX_Weight
Definition: VarVals.h:56
void reconfigure(fhicl::ParameterSet const &p)
const Var kRPAWeightRES2017
Definition: GenieWeights.h:115
static const unsigned char kTrueHitNuc
Definition: VarVals.h:52
::xsd::cxx::tree::exception< char > exception
Definition: Database.h:225
std::map< fnex::MetaData, fnex::EventList > EventListMap
Definition: Event.h:186
virtual void Shift(double sigma, caf::SRProxy *sr, double &weight) const
Perform the systematic shift.
Definition: ISyst.cxx:26
int PeriodAdjustCosmics(int const &inputPeriod)
unsigned int run
run number
Definition: SRSpill.h:25
const Cut kNue2018NDEnergy
Definition: NueCuts2018.h:143
Create a list of fnex::Events to be used in fits.
Timing fit.
const Cut kNue2018CVNCut([](const caf::SRProxy *sr){if(kIsRHC(sr)) return kNue2018RHCCVNCut(sr);else return kNue2018FHCCVNCut(sr);})
Definition: NueCuts2018.h:49
static const unsigned char kTrueParentPDG
Definition: VarVals.h:47
std::string find_file(std::string const &filename) const
const Var kEmpiricalMECWgt2018
See kEmpiricalMECWgt2018_NT.
Definition: GenieWeights.h:254
const Cut kNue2018NDPresel
Definition: NueCuts2018.h:145
void FindPOTPerPeriod(TFile &tf)
float abs(float number)
Definition: d0nt_math.hpp:39
const Var kTruetpT([](const caf::SRProxy *sr){return(sr->mc.nnu==0)?0.:util::pythag(sr->mc.nu[0].beam.tp.X(), sr->mc.nu[0].beam.tp.Y());})
True neutrino ancestor (off-target) pT.
Definition: Vars.h:102
novadaq::cnv::DetId fDetector
which detector are we using
static const unsigned char kRPARES_Weight
Definition: VarVals.h:58
DEFINE_ART_MODULE(SpectrumTest)
caf::Proxy< unsigned int > run
Definition: SRProxy.h:243
const Var kHitNuc
Definition: TruthVars.h:19
caf::Proxy< caf::SRCosRej > cosrej
Definition: SRProxy.h:1190
int isnan(const stan::math::var &a)
Definition: std_isnan.hpp:18
const Var kRPAWeightCCQE2018
Definition: GenieWeights.h:96
bool PassesSpillAndBeamCuts(caf::SRProxy const &proxy) const
const Cut kNumuCosmicRej2018([](const caf::SRProxy *sr){return(sr->sel.cosrej.anglekal > 0.5 && sr->sel.cosrej.numucontpid2019 > 0.53 && sr->slc.nhit< 400 && sr->sel.nuecosrej.pngptp< 0.9 );})
Definition: NumuCuts2018.h:19
static const unsigned char kTruePDGOrig
Definition: VarVals.h:46
caf::Proxy< caf::SRNueCosRej > nuecosrej
Definition: SRProxy.h:1199
static const unsigned char kTrueCCNC
Definition: VarVals.h:44
caf::Proxy< bool > isFHC
Definition: SRProxy.h:1307
const XML_Char * s
Definition: expat.h:262
unsigned int subrun
subrun number
Definition: SRSpill.h:26
std::string const ToString() const
Definition: Structs.cxx:114
Far Detector at Ash River, MN.
void FillVariables(std::vector< fnex::DataVarVals > &dataVals, std::vector< fnex::MCVarVals > &mcVals, std::vector< fnex::EventId > &evIds, std::vector< fnex::SelectionType_t > &sels)
caf::Proxy< unsigned int > nhit
Definition: SRProxy.h:1247
std::string fRHCQuantFileName
name of file in numudata ups product with FHC quantiles
const Var kTrueE([](const caf::SRProxy *sr){return(sr->mc.nnu==0)?0.:float(sr->mc.nu[0].E);})
Definition: Vars.cxx:85
enum fnex::sel_type SelectionType_t
const Cut kNumuContainFD2017
Definition: NumuCuts2017.h:21
const Cut sels[kNumSels]
Definition: vars.h:44
caf::Proxy< short unsigned int > subevt
Definition: SRProxy.h:245
const Var kTruetpz([](const caf::SRProxy *sr){return(sr->mc.nnu==0)?0.:float(sr->mc.nu[0].beam.tp.z);})
True neutrino ancestor (off-target) pz.
Definition: Vars.h:105
fnex::SelectionType_t FindSelectionType(caf::SRProxy const &srProxy) const
static const unsigned char kLep_RecoE_MCFrac
Definition: VarVals.h:37
std::string fPOTHistName
names of the CAF histogram with the total POT
const Cut kNumuOptimizedContainFD2017([](const caf::SRProxy *sr){std::pair< int, int > planes=calcFirstLastLivePlane(sr->slc.firstplane, std::bitset< 14 >(sr->hdr.dibmask));int planestofront=sr->slc.firstplane-planes.first;int planestoback=planes.second-sr->slc.lastplane;return(sr->sel.contain.kalfwdcell > 6 &&sr->sel.contain.kalbakcell > 6 &&sr->sel.contain.cosfwdcell > 0 &&sr->sel.contain.cosbakcell > 7 &&planestofront > 1 &&planestoback > 1);})
Definition: NumuCuts2017.h:16
T get(std::string const &key) const
Definition: ParameterSet.h:231
const Var kCCE
Definition: NumuVars.h:21
static const unsigned char kNuE_NumMichel
Definition: VarVals.h:33
const Var kEmpiricalMECtoGENIERESWgt
Definition: GenieWeights.h:184
unsigned int evt
ART event number, indexes trigger windows.
Definition: SRSpill.h:27
const Var kEmpiricalMECtoGENIEQEWgt
Definition: GenieWeights.h:180
Near Detector in the NuMI cavern.
const Var kRecoEME([](const caf::SRProxy *sr){if(!sr->vtx.elastic.IsValid) return-1.0;if(sr->vtx.elastic.fuzzyk.npng< 1) return-1.0;double CVNem_CalE=0.0;const caf::SRFuzzyKProngProxy &prim_png=sr->vtx.elastic.fuzzyk.png[0];for(const caf::SRFuzzyKProngProxy &png:sr->vtx.elastic.fuzzyk.png){double png_CalE=png.shwlid.calE;double emPID=((double) png.cvnpart.photonid+ (double) png.cvnpart.pizeroid+ (double) png.cvnpart.electronid);double haPID=((double) png.cvnpart.protonid+ (double) png.cvnpart.pionid+ (double) png.cvnpart.neutronid+ (double) png.cvnpart.otherid+ (double) png.cvnpart.muonid);if(emPID< 0) continue;if(emPID >=haPID ) CVNem_CalE+=png_CalE;else continue;}if(CVNem_CalE==0.0) CVNem_CalE=prim_png.shwlid.calE;return CVNem_CalE *CalibrationBugCorrectionFactor(sr->hdr);})
Definition: NueEnergy2018.h:17
bool fMakeSelectedList
configure to true to make a list
const Cut kNumuProngsContainFD2017
Definition: NumuCuts2017.h:19
caf::Proxy< int > cycle
Definition: SRProxy.h:225
static const unsigned char kTrueIntType
Definition: VarVals.h:45
bool fIsFHC
are we looking at FHC files?
caf::Proxy< caf::SRRemid > remid
Definition: SRProxy.h:1202
void FillDataVarVals(fnex::DataVarVals &dataVals, caf::SRProxy const &srProxy) const
static const unsigned char kTrueParentDecay
Definition: VarVals.h:50
const Var kCVNSSe([](const caf::SRProxy *sr){throw std::runtime_error("kCVNSSe is no longer available. Fix your macro so you don't use it.");return-5.;})
2018 nue PID
Definition: Vars.h:52
bool PassesNumuSelections(caf::SRProxy const &proxy) const
const HistAxis kNumuCCOptimisedAxis("Reconstructed Neutrino Energy (GeV)", kNumuCCEOptimisedBinning, kCCE)
HistAxis that implements optimised numuCCE from L. Vinton. See docdb 16332. This was close to &#39;custC&#39;...
Definition: HistAxes.h:24
const Cut kNue2018FDPeripheral(kNue2018FDPeripheralFunc)
const Var kNMichels([](const caf::SRProxy *sr){int n_me=0;for(int i=0;i< (int) sr->me.nslc;i++) if(sr->me.slc[i].mid > 1. &&sr->me.slc[i].deltat > 1200.) n_me++;for(int i=0;i< (int) sr->me.nkalman;i++) if(sr->me.trkkalman[i].mid > 0. &&sr->me.trkkalman[i].deltat > 800.) n_me++;if(n_me > 2) n_me=2;return n_me;})
Number of SlcME&#39;s assoc with slice.
Definition: Vars.h:85
std::map< int, POTLiveTime > fPeriodToPOTLiveTime
keep track of the POT in each period
Definition: run.py:1
#define LOG_WARNING(category)
caf::Proxy< unsigned int > evt
Definition: SRProxy.h:232
void FillMCVarVals(fnex::MCVarVals &mcVarVals, caf::SRProxy &srProxy) const
static float min(const float a, const float b, const float c)
Definition: absgeo.cxx:45
static const unsigned char kTrueParentTargetPDG
Definition: VarVals.h:51
float livetime
Length of readout [s].
Definition: SRSpill.h:40
void set_GENIEWgt(unsigned char const &varkey, std::string const &sigma, float const &wgt)
Definition: VarVals.cxx:282
caf::Proxy< caf::SRTruthBranch > mc
Definition: SRProxy.h:2005
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
The StandardRecord is the primary top-level object in the Common Analysis File trees.
std::vector< std::string > fCAFFileNames
name of the CAF file we are using
caf::Proxy< caf::SRSlice > slc
Definition: SRProxy.h:2009
static const unsigned char kEmpiricalMECtoGENIEQE_Weight
Definition: VarVals.h:61
static const unsigned char kRPACCQE_Weight
Definition: VarVals.h:57
void FillEventId(fnex::EventId &evId, caf::SRProxy const &srProxy) const
std::shared_ptr< Event > EventPtr
Definition: Event.h:50
const Cut kNumuPID2018([](const caf::SRProxy *sr){std::cout<< "ERROR::kNumuPID2018, cutting on both cvnProd3Train and cvn2017."<< " Neither branch exists anymore. Returning False."<< std::endl;abort();return false;})
Definition: NumuCuts2018.h:22
const Var kXSecCVWgt2018
Definition: XsecTunes.h:49
static const unsigned char kTrueE
Definition: VarVals.h:42
caf::Proxy< float > anglekal
Definition: SRProxy.h:811
caf::Proxy< bool > ismc
Definition: SRProxy.h:237
const Cut kNue2018FDAllSamples
Definition: NueCuts2018.h:77
caf::Proxy< float > pngptp
Definition: SRProxy.h:999
std::vector< Cut > QuantileCutsFromTH2(TH2 *quantileHist, const HistAxis &independentAxis, const HistAxis &quantileAxis, const unsigned int &numQuantiles, const bool verbose)
: Do the same as the QuantileCuts function but taking in the TH2 instead of making it...
const Cut kInBeamSpill([](const caf::SRProxy *sr){if(sr->spill.run > util::kLastBadTimingRun) return kInBeamSpill_main(sr);else return kInBeamSpill_main(sr)||kInBeamSpill_shifted(sr);}, [](const caf::SRSpillProxy *spill){if(spill->run > util::kLastBadTimingRun) return kInBeamSpill_main.Livetime(spill);else return kInBeamSpill_main.Livetime(spill)+kInBeamSpill_shifted.Livetime(spill);}, [](const caf::SRSpillProxy *spill) -> double{return spill->spillpot;})
Does the event fall inside the window we call the beam spill?
Definition: TimingCuts.h:8
void set_val_at(std::string const &varkey, float const &val)
Definition: VarVals.cxx:525
const Cut kNumuCutFD2018
Definition: NumuCuts2018.h:39
void operator+=(POTLiveTime const &o)
const NOvARwgtSyst kRPACCQESuppSyst2018("RPAShapesupp2018","RPA shape: low-Q^{2} suppression (2018)", novarwgt::kRPACCQESuppSyst2018)
Definition: RPASysts.h:14
const double kNue2018CVNRHCHighEdge
Definition: NueCuts2018.h:44
static const unsigned char kHad_RecoE
Definition: VarVals.h:34
std::vector< ana::Cut > fNumuFHCQuantCuts
quantile boundaries for the numu FHC data
const SpillCut kStandardSpillCuts
Apply this unless you&#39;re doing something special.
Definition: SpillCuts.h:46
bool fIsRealData
are we looking at real data?
static const unsigned char kTrueIntMode
Definition: VarVals.h:53
float spillpot
Definition: SRSpill.h:38
caf::Proxy< caf::SRIDBranch > sel
Definition: SRProxy.h:2008
#define LOG_VERBATIM(category)
const std::string cSelectionType_Strings[11]
Definition: Constants.h:101
static const unsigned char kRPACCQEshapeEnh
Definition: VarVals.h:74
fnex::EventListMap fMDToEvents
each entry in the vector is for a single slice
bool PassesNueSelections(caf::SRProxy const &proxy) const
const Cut kNumuQuality
Definition: NumuCuts.h:18
Float_t e
Definition: plot.C:35
static const unsigned char kLep_RecoE
Definition: VarVals.h:35
bool PassesSelections(caf::SRProxy const &proxy) const
std::set< fnex::SelectionType_t > fSelections
Selections to use.
std::string RunToEpoch(int run)
CAFToEventList(fhicl::ParameterSet const &pset)
static const unsigned char kTrueParentPZ
Definition: VarVals.h:49
Base container for the MC related Vars that constitute an event.
Definition: VarVals.h:874
static constexpr Double_t sr
Definition: Munits.h:164
detName
Definition: mkDefs.py:106
fnex::FileType_t fFileType
get from the configuration for now
const Var kRecoHADE([](const caf::SRProxy *sr){if(!sr->vtx.elastic.IsValid) return-1.0;if(sr->vtx.elastic.fuzzyk.npng< 1) return-1.0;const double CVNha_CalE=sr->slc.calE *CalibrationBugCorrectionFactor(sr->hdr);return std::max(CVNha_CalE-kRecoEME(sr), 0.);})
Definition: NueEnergy2018.h:18
static const unsigned char kTruePDG
Definition: VarVals.h:43