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