CMFDSTToEventList_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 // History
7 // - The last production 4 commit is 42775
8 
9 #include <memory>
10 #include <cmath>
11 #include <string>
12 #include <vector>
13 #include <set>
14 #include <limits>
15 
16 // Framework includes
22 #include "art_root_io/TFileService.h"
23 #include "art_root_io/TFileDirectory.h"
25 #include "fhiclcpp/ParameterSet.h"
27 #include "cetlib_except/exception.h"
28 #include "cetlib/search_path.h"
29 
30 // NOvA includes
31 #include "SummaryData/POTSum.h"
36 #include "nugen/NuReweight/ReweightLabels.h"
39 
42 
43 // ROOT includes
44 #include "TTree.h"
45 #include "TFile.h"
46 #include "TH1.h"
47 #include "TRandom3.h"
48 
49 namespace cmf {
50 
51  enum minosRP {
66  };
67 
68  struct POTLiveTime{
69 
71  : pot (-1.)
72  , liveTime(-1.)
73  {}
74 
75  POTLiveTime(double const& p,
76  double const& l)
77  : pot (p)
78  , liveTime(l)
79  {}
80 
81  double pot;
82  double liveTime;
83 
84  void operator+=(POTLiveTime const& o) { pot = o.pot; liveTime = o.liveTime; }
85 
86  };
87 
88  struct PeriodRun{
89 
91  int startRun,
92  int endRun)
93  : fPeriod(period)
94  , fStartRun(startRun)
95  , fEndRun(endRun)
96  {}
97 
98  bool RunInPeriod(int run) const { return ((run > fStartRun - 1) && (run < fEndRun + 1)); }
99 
100  int fPeriod;
101  int fStartRun;
102  int fEndRun;
103 
104  };
105 
107  public:
108  explicit DSTToEventList(fhicl::ParameterSet const& pset);
109  virtual ~DSTToEventList();
110 
111  // Plugins should not be copied or assigned.
112  DSTToEventList(DSTToEventList const &) = delete;
113  DSTToEventList(DSTToEventList &&) = delete;
114  DSTToEventList & operator = (DSTToEventList const &) = delete;
115  DSTToEventList & operator = (DSTToEventList &&) = delete;
116 
117  void analyze(art::Event const& e) override;
118  void reconfigure(fhicl::ParameterSet const& p) ;
119  void endJob() override;
120 
121  private:
122 
123  void clear ();
124  void InitializeEventListColls();
125  int RunToPeriod (int run);
126  void FillVariables (std::vector<cmf::DataVarVals> & dataVals,
127  std::vector<cmf::MCVarVals> & mcVals,
128  std::vector<cmf::EventId> & evIds,
129  std::vector<cmf::SelectionType_t> & sels);
130  void FillMCVals (cmf::TruthVars & truthVars,
131  cmf::WeightVars & weightVars,
132  cmf::SystVarColl & systVars,
133  NuMCEvent & srProxy) const;
134  void FillDataVarVals (cmf::DataVarVals & dataVals,
135  NuEvent const& srProxy) const;
136  void FillEventId (cmf::EventId & evId,
137  NuEvent const& srProxy) const;
138  void FindPOTPerPeriod (TFile & tf);
139  cmf::SelectionType_t FindSelectionType (NuEvent const& srProxy) const;
140 
141  std::vector<std::string> fCAFFileNames; ///< name of the CAF file we are using
142  std::string fFHCQuantFileName; ///< name of file in numudata ups product with FHC quantiles
143  std::string fRHCQuantFileName; ///< name of file in numudata ups product with FHC quantiles
144  std::string fSpillTreeName; ///< names of the Spill tree in the file
145  std::string fPOTHistName; ///< names of the CAF histogram with the total POT
146  std::string fCAFRecTreeName; ///< names of the CAF tree with the reco info
147  cmf::DetType_t fDetector; ///< which detector are we using
148  cmf::EventListColl fMDToEvents; ///< each entry in the vector is for a single slice
149  cmf::FileType_t fFileType; ///< get from the configuration for now
150  std::string fSystType; ///< systematic type
151  bool fIsRealData; ///< are we looking at real data?
152  bool fIsFHC; ///< are we looking at FHC files?
153  unsigned int fNumQuantiles; ///< number of quantiles
154  //std::set<cmf::EventId> fMissingEvents; ///< events missing compared to CAF sample
155  std::map<int, POTLiveTime> fPeriodToPOTLiveTime; ///< keep track of the POT in each period
156  std::set<cmf::SelectionType_t> fSelections; ///< Selections to use
157  bool fMakeSelectedList; ///< configure to true to make a list
158  std::vector<cmf::PeriodRun> fPeriodRuns; ///< vector of periods to runs
159  };
160 
161  //......................................................................
163  : EDAnalyzer (pset)
164  , fDetector (cmf::kNEARDET)
165  , fFileType (cmf::kUnknownFileType)
166  , fIsRealData (false)
167  , fIsFHC (true)
168  , fNumQuantiles (1)
169  , fMakeSelectedList(false)
170  {
171  // make sure everything is zeroed out.
172  this->clear();
173 
174  this->reconfigure(pset);
175 
176  // initialize the period to POT map
177  // just allow for many periods, unlikely we will
178  // have more than 50, but it doesn't matter if we
179  // plan ahead
180  for(int i = 1; i < 51; ++i) fPeriodToPOTLiveTime[i] = POTLiveTime(0., 0.);
181 
182  // fMissingEvents.insert(cmf::EventId(19195, 10, 73, 0, 37));
183  // fMissingEvents.insert(cmf::EventId(19546, 12, 222, 0, 30));
184  // fMissingEvents.insert(cmf::EventId(19580, 19, 584, 0, 36));
185  // fMissingEvents.insert(cmf::EventId(19662, 58, 54, 0, 27));
186  // fMissingEvents.insert(cmf::EventId(21650, 54, 267, 0, 28));
187  // fMissingEvents.insert(cmf::EventId(22797, 43, 678, 0, 22));
188  // fMissingEvents.insert(cmf::EventId(23330, 24, 252, 0, 33));
189 
190  }
191 
192  //......................................................................
194  {
195  this->clear();
196  }
197 
198  //......................................................................
199  // Method to clear out the collections of data products after the
200  // writeResults method is called.
202  {
203  fMDToEvents .clear();
204  fSelections .clear();
205  }
206 
207  //......................................................................
209  {
210  fCAFFileNames = p.get<std::vector<std::string> >("CAFFileNames" );
211 
212  // assume we don't mix FHC and RHC files in the same job, and the
213  // fIsFHC should only be used for cosmics
214  if(fCAFFileNames.begin()->find("rhc") != std::string::npos)
215  fIsFHC = false;
216 
217  fCAFRecTreeName = p.get<std::string >("CAFRecTreeName" , "s" );
218  fSpillTreeName = p.get<std::string >("SpillTreeName" , "spillTree" );
219  fPOTHistName = p.get<std::string >("CAFPOTHistName" , "TotalPOT" );
220  fNumQuantiles = p.get<unsigned int>("NumQuantiles" , 4 );
221  fMakeSelectedList = p.get<bool >("MakeSelectedList" , false );
222  fSystType = p.get<std::string >("SystematicType" , "NOSYST" );
223 
224  auto detName = p.get<std::string>("DetectorName", "MINOSNear");
225  fDetector = (detName == "MINOSNear") ? cmf::kMINOSNEAR : cmf::kMINOSFAR;
226 
227  auto selType = p.get<std::string>("Selections");
228 
229 
230  if(selType == "NuMuSel"){
232  }
233  else{
235  }
236 
237  auto ftString = p.get<std::string>("FileType");
238 
239  if (ftString == "FluxSwap") fFileType = kSwap;
240  else if(ftString == "NonSwap" ) fFileType = kBeam;
241  else if(ftString == "TauSwap" ) fFileType = kTauSwap;
242  else if(ftString == "CosmicBackground") fFileType = kCosmicBackgroundFile;
243  else if(ftString == "RockFluxSwap" ) fFileType = kRockFluxSwap;
244  else if(ftString == "RockNonSwap" ) fFileType = kRockNonSwap;
245  else{
247  fIsRealData = true;
248  } // end tests on ftString
249 
250  // load up the epochs to use
251  auto periodPOTPars = p.get< std::vector<fhicl::ParameterSet> >("PeriodsToUse"); //Setting NOvA periods to use based on a fhicl parameter with definitions in CMF_Periods.fcl
252  //TODO: make fhicl file for minos runs
253 
254  for(auto const& itr : periodPOTPars){//TODO: might need to be more careful about the run numbers but this will do for now
255  std::cout << itr.get<std::string>("Detector") << " " << itr.get<int>("Period")
258  if(cmf::StringToDetectorType(itr.get<std::string>("Detector")) != fDetector) continue;
259 
260  fPeriodRuns.emplace_back(itr.get<cmf::minosRP>("Period"),
261  std::numeric_limits<int>::min(),//itr.get<int>("StartRun"),
262  std::numeric_limits<int>::max());//itr.get<int>("EndRun"));
263  }
264 
265  fPeriodRuns.emplace_back(kRun11,
266  std::numeric_limits<int>::min(),//itr.get<int>("StartRun"),
268 
270 
271  }
272 
273  //......................................................................
274  // since we are using CAF files, which are not art files, we don't need the
275  // analyze function to do anything
277  {
278  }
279 
280 
281 
282 
283  //......................................................................
284  // see https://cdcvs.fnal.gov/redmine/projects/novaart/wiki/Period_and_Epoch_Naming
286  {
287  std::cout << "run is " << run << std::endl;
288  for(auto const& itr : fPeriodRuns){
289  std::cout << itr.fPeriod << std::endl;
290  if(itr.RunInPeriod(run)) return itr.fPeriod;
291  }
292 
293  return kRun11;
294 
295  // no period found for this run and detector, got a problem
296  throw cet::exception("RunToPeriod")
297  << "no period found for run "
298  << run;
299 
300  return 0;
301  }
302 
303  //......................................................................
305  {
306  // here we will create an event list for each possible meta data
307  // for the current epoch and data/MC type.
308 
309  // if we are looking at the ND, there are no cosmic or rock events
310  // there, so don't bother
311  if(fDetector == cmf::kNEARDET &&
314  fFileType == cmf::kRockNonSwap) ) return;
315 
317 
318  std::set<cmf::InteractionType_t> intTypes({cmf::kNuMuCC,
320  cmf::kNuECC,
324  cmf::kNC,
326 
327  // set the bits of the metadata that won't change
328  bool isMC = (fFileType != cmf::kDataFile);
329 
330  for(auto const& itr : fPeriodRuns){//This is empty for minos at the moment
331  // loop over the selection and interaction types to flesh out the lists
332  for(auto const& s : fSelections){
333 
334  // there are no nue peripheral events in the ND
335  if(fDetector == cmf::kNEARDET &&
336  s == cmf::kNuESelectionPeripheral) continue;
337 
338  if(fFileType == cmf::kDataFile){
339 
340  // data have an unknown interaction type
341  cmf::MetaData md(isMC,
342  fDetector,
343  fFileType,
344  s,
346  itr.fPeriod);
347 
348  // make a new EventList for the data
349  if(cmf::FindEventList(md, fMDToEvents) == fMDToEvents.end())
350  fMDToEvents.emplace_back(cmf::EventList(md, ss));
351  std::cout<< "added " << md.ToString() << std::endl;
352  } // end if data
354 
355  cmf::MetaData md(true,
356  fDetector,
357  fFileType,
358  s,
360  itr.fPeriod);
361 
362  // make a new EventList for the cosmic background
363  if(cmf::FindEventList(md, fMDToEvents) == fMDToEvents.end())
364  fMDToEvents.emplace_back(cmf::EventList(md, ss));
365  } // end if cosmic background
366  else{
367  // loop over possible interaction types for the metadata
368  // don't use cosmic muons, those are special and are handled
369  // in another way - we don't want to inflate the POT for those lists
370  for(auto const& i : intTypes){
371 
372  cmf::MetaData md(isMC,
373  fDetector,
374  fFileType,
375  s,
376  i,
377  itr.fPeriod);
378 
379  // make a new EventList for the MC
380  if(cmf::FindEventList(md, fMDToEvents) == fMDToEvents.end())
381  fMDToEvents.emplace_back(cmf::EventList(md, ss));
382  } // end loop over interaction types
383  } // end if not data or MC, ie cosmic background
384 
385  } // end loop over selection types
386  } // end loop over epochs
387  }
388 
389 
390 
391  //......................................................................
392  void DSTToEventList::FindPOTPerPeriod(TFile & tf)//Modify for MINOS (minimally)
393  {
394  int period;
395 
396  auto *spillTree = dynamic_cast<TTree*>(tf.Get(fCAFRecTreeName.c_str()));
397 
398  NuEvent *sp = nullptr;
399 
400  spillTree->SetBranchAddress("s", &sp);
401  spillTree->GetEntry(0);
402 
403  MF_LOG_VERBATIM("DSTToEventList")
404  << "There are "
405  << spillTree->GetEntriesFast()
406  << " reconstructed events";
407 
408 
409 
410 
411  // now figure out which period this spill belongs to
412  period = sp->runPeriod; //this might not work in all cases since minos wasn't careful with the pHE RunI period
413 
415  MF_LOG_VERBATIM("DSTToEventList")
416  << "This is a FD data file"
417  << " ... "
418  << "Filling in POT by hand";
419 
420  if(period == 11)
421  fPeriodToPOTLiveTime[period].pot += 2.9852e20;
422  else
424 
425  }
426  else{
427  if(std::isnan(((TH1F*) tf.Get(fPOTHistName.c_str()))->Integral())){
428  MF_LOG_VERBATIM("DSTToEventList")
429  << tf.GetName()
430  << " "
431  << " file has nan for pot or livetime "
432  << ((TH1F*) tf.Get(fPOTHistName.c_str()))->Integral()
433  << " "
434  << " period "
435  << period;
436 
437  }
438  else
439  fPeriodToPOTLiveTime[period].pot += ((TH1F*)tf.Get(fPOTHistName.c_str()))->Integral();
440  }
442 
443 
444 
445 
446  // print out the number of pot for all files so far
447  for(auto const& itr : fPeriodToPOTLiveTime){
448  if(itr.second.pot > 0. ||
449  itr.second.liveTime > 0.)
450  MF_LOG_VERBATIM("DSTToEventList")
451  << "period "
452  << itr.first
453  << " has "
454  << itr.second.pot * 1.e-20
455  << " x 10^20 POT and "
456  << itr.second.liveTime
457  << " s live time";
458  }
459  }
460 
461  //......................................................................
463  {
464 /* if(this->PassesNumuSelections(srProxy)){
465  if(srProxy.spill.isFHC){
466  if (fNumuFHCQuantCuts[0](&srProxy)) return cmf::kNuMuSelectionQ1;
467  else if(fNumuFHCQuantCuts[1](&srProxy)) return cmf::kNuMuSelectionQ2;
468  else if(fNumuFHCQuantCuts[2](&srProxy)) return cmf::kNuMuSelectionQ3;
469  else if(fNumuFHCQuantCuts[3](&srProxy)) return cmf::kNuMuSelectionQ4;
470  } // end if fhc
471  else{
472  if (fNumuRHCQuantCuts[0](&srProxy)) return cmf::kNuMuSelectionQ1;
473  else if(fNumuRHCQuantCuts[1](&srProxy)) return cmf::kNuMuSelectionQ2;
474  else if(fNumuRHCQuantCuts[2](&srProxy)) return cmf::kNuMuSelectionQ3;
475  else if(fNumuRHCQuantCuts[3](&srProxy)) return cmf::kNuMuSelectionQ4;
476  }
477  } // end if numu
478  else if(this->PassesNueSelections(srProxy)){
479  auto cvnHighEdge = (srProxy.spill.isFHC) ? ana::kNue2020CVNFHCHighEdge : ana::kNue2020CVNRHCHighEdge;
480  if(ana::kNue2020FDPeripheral(&srProxy))
481  return cmf::kNuESelectionPeripheral;
482  if(ana::kCVNe_looseptp(&srProxy) > cvnHighEdge)
483  return cmf::kNuESelectionHighPID;
484 
485  // event is NuE and passes the CVN cut, but isn't high or peripheral,
486  // it must be low
487  return cmf::kNuESelectionLowPID;
488  }
489  else if(this->PassesNCSelections(srProxy)) return cmf::kNCSelection;*/
490 
491 
492  return cmf::kUnknownSelection;
493  }
494 
495  //......................................................................
496  // note that srProxy is not const here because some of the CAF/Var functions
497  // do not take const pointers to it...seems silly to not have it const.
499  cmf::WeightVars & weightVars,
500  cmf::SystVarColl & systShifts,
501  NuMCEvent & srProxy) const
502  {
503 /*
504  truthVars.fTrueE = (float)ana::kTrueE(&srProxy) ;
505  truthVars.fTruePDG = (float)ana::kTruePDG(&srProxy) ;
506  truthVars.fTrueCCNC = (srProxy.mc.nu[0].iscc) ? simb::kCC : simb::kNC;
507  truthVars.fTrueIntType = srProxy.mc.nu[0].inttype ;
508  truthVars.fTrueIntMode = (float)ana::kMode(&srProxy) ;
509  truthVars.fTrueHitNuc = (float)ana::kHitNuc(&srProxy) ;
510  truthVars.fTruePDGOrig = srProxy.mc.nu[0].pdgorig ;
511  truthVars.fTrueParentPDG = srProxy.mc.nu[0].beam.ptype ;
512  truthVars.fTrueParentDecay = srProxy.mc.nu[0].beam.ndecay ;
513  truthVars.fTrueParentTargetPDG = srProxy.mc.nu[0].beam.tptype ;
514  truthVars.fTrueParentPZ = (float)ana::kTruetpz(&srProxy) ;
515  truthVars.fTrueParentPT = (float)ana::kTruetpT(&srProxy) ;
516 
517  // we only care about the following systematics if we're not running a
518  // systematically shifted dataset
519  if (fSystType != "NOSYST") return;
520 
521  // these are weights applied to the event
522  weightVars.fXSecCVPPFX_Weight = (float)(ana::kXSecCVWgt2020(&srProxy) * ana::kPPFXFluxCVWgt(&srProxy));
523  weightVars.fXSecCV2020_Weight = (float)(ana::kXSecCVWgt2020(&srProxy));
524  weightVars.fPPFXFluxCV_Weight = (float)(ana::kPPFXFluxCVWgt(&srProxy));
525 
526  // Now for the full set of systematic knobs
527  // These systematics for the 2020 analysis are taken from Analysis/3FlavorAna2020Systs.cxx
528  //
529  // Note that we don't take:
530  // -- file systematics -> just produce systematically shifted event lists
531  // -- beam systematics -> we have our own handling of that. See core/ShifterAndWeighter.cxx
532 
533  static std::map<const ana::ISyst*, std::string> labelMap =
534  {
535  // xsec systematics
536  {ana::GetGenieKnobSyst(rwgt::fReweightMaCCRES) , "MaCCRES" },
537  {ana::GetGenieKnobSyst(rwgt::fReweightMvCCRES) , "MvCCRES" },
538  {ana::GetGenieKnobSyst(rwgt::fReweightMaNCRES) , "MaNCRES" },
539  {ana::GetGenieKnobSyst(rwgt::fReweightMvNCRES) , "MvNCRES" },
540  {ana::GetGenieKnobSyst(rwgt::fReweightZNormCCQE), "ZNormCCQE" },
541  {&ana::kZExpEV1Syst2020 , "ZExpAxialFFSyst2020_eV1" },
542  {&ana::kZExpEV2Syst2020 , "ZExpAxialFFSyst2020_eV2" },
543  {&ana::kZExpEV3Syst2020 , "ZExpAxialFFSyst2020_eV3" },
544  {&ana::kZExpEV4Syst2020 , "ZExpAxialFFSyst2020_eV4" },
545  {&ana::kRPACCQEEnhSyst2020 , "RPACCQEshapeEnh" }, // ?
546  {&ana::kRPACCQESuppSyst2020 , "RPACCQEshapeSupp" }, // ?
547  {&ana::kRESLowQ2SuppressionSyst2020 , "RESLowQ2SuppSyst2020" },
548  {&ana::kDISvnCC1pi_2020 , "DISvnCC1piWgt_2020" }, // ?
549  {&ana::khNFSISyst2020_MFP , "hNFSISyst2020_MFP" },
550  {&ana::khNFSISyst2020_EV1 , "hNFSISyst2020_EV1" },
551  {&ana::kMECEnuShapeSyst2020Nu , "MECEnuShapeSyst2020Nu" },
552  {&ana::kMECEnuShapeSyst2020AntiNu , "MECEnuShapeSyst2020AntiNu" },
553  {&ana::kMECShapeSyst2020Nu , "MECShapeSyst2020Nu" },
554  {&ana::kMECShapeSyst2020AntiNu , "MECShapeSyst2020AntiNu" },
555  {&ana::kMECInitStateNPFracSyst2020Nu , "MECInitStateNPFracSyst2020Nu" },
556  {&ana::kMECInitStateNPFracSyst2020AntiNu , "MECInitStateNPFracSyst2020AntiNu" },
557  {&ana::kRadCorrNue , "RadCorrNue" },
558  {&ana::kRadCorrNuebar , "RadCorrNuebar" },
559  {&ana::k2ndClassCurrs , "2ndClassCurrs" },
560  // muon energy scale systematics
561  {&ana::kUnCorrNDMuEScaleSyst2020 , "UnCorrNDMuEScaleSyst2020" },
562  {&ana::kUnCorrMuCatMuESyst2020 , "UnCorrMuCatMuESyst2020" },
563  {&ana::kPileupMuESyst2020 , "PileupMuESyst2020" },
564  {&ana::kCorrMuEScaleSyst2020 , "CorrMuEScaleSyst2020" },
565  // neutron systematics
566  {&ana::kNeutronVisEScalePrimariesSyst2018 , "NeutronVisEScalePrimariesSyst2018"},
567  // michel systematics
568  {&ana::kMichelTaggingSyst2020 , "MichelTaggingSyst2020" },
569  // norm systemaitics
570  {&ana::kAna2020NormFHC , "Ana2020NormFHC" },
571  {&ana::kAna2020NormRHC , "Ana2020NormRHC" },
572  // other systematics
573  {&ana::kTauScaleSyst , "TauScaleSyst" },
574  // loads of DIS systs, likely in pca systs, but worth having just in case
575  {&ana::kDISvpCC0pi_2020 , "DISvpCC0pi_2020" },
576  {&ana::kDISvpCC1pi_2020 , "DISvpCC1pi_2020" },
577  {&ana::kDISvpCC2pi_2020 , "DISvpCC2pi_2020" },
578  {&ana::kDISvpCC3pi_2020 , "DISvpCC3pi_2020" },
579  {&ana::kDISvpNC0pi_2020 , "DISvpNC0pi_2020" },
580  {&ana::kDISvpNC1pi_2020 , "DISvpNC1pi_2020" },
581  {&ana::kDISvpNC2pi_2020 , "DISvpNC2pi_2020" },
582  {&ana::kDISvpNC3pi_2020 , "DISvpNC3pi_2020" },
583  {&ana::kDISvnCC0pi_2020 , "DISvnCC0pi_2020" },
584  {&ana::kDISvnCC1pi_2020 , "DISvnCC1pi_2020" },
585  {&ana::kDISvnCC2pi_2020 , "DISvnCC2pi_2020" },
586  {&ana::kDISvnCC3pi_2020 , "DISvnCC3pi_2020" },
587  {&ana::kDISvnNC0pi_2020 , "DISvnNC0pi_2020" },
588  {&ana::kDISvnNC1pi_2020 , "DISvnNC1pi_2020" },
589  {&ana::kDISvnNC2pi_2020 , "DISvnNC2pi_2020" },
590  {&ana::kDISvnNC3pi_2020 , "DISvnNC3pi_2020" },
591  {&ana::kDISvbarpCC0pi_2020 , "DISvbarpCC0pi_2020" },
592  {&ana::kDISvbarpCC1pi_2020 , "DISvbarpCC1pi_2020" },
593  {&ana::kDISvbarpCC2pi_2020 , "DISvbarpCC2pi_2020" },
594  {&ana::kDISvbarpCC3pi_2020 , "DISvbarpCC3pi_2020" },
595  {&ana::kDISvbarpNC0pi_2020 , "DISvbarpNC0pi_2020" },
596  {&ana::kDISvbarpNC1pi_2020 , "DISvbarpNC1pi_2020" },
597  {&ana::kDISvbarpNC2pi_2020 , "DISvbarpNC2pi_2020" },
598  {&ana::kDISvbarpNC3pi_2020 , "DISvbarpNC3pi_2020" },
599  {&ana::kDISvbarnCC0pi_2020 , "DISvbarnCC0pi_2020" },
600  {&ana::kDISvbarnCC1pi_2020 , "DISvbarnCC1pi_2020" },
601  {&ana::kDISvbarnCC2pi_2020 , "DISvbarnCC2pi_2020" },
602  {&ana::kDISvbarnCC3pi_2020 , "DISvbarnCC3pi_2020" },
603  {&ana::kDISvbarnNC0pi_2020 , "DISvbarnNC0pi_2020" },
604  {&ana::kDISvbarnNC1pi_2020 , "DISvbarnNC1pi_2020" },
605  {&ana::kDISvbarnNC2pi_2020 , "DISvbarnNC2pi_2020" },
606  {&ana::kDISvbarnNC3pi_2020 , "DISvbarnNC3pi_2020" },
607  // and some additional likely-pca systs
608  {&ana::kFormZone_2020 , "FormZone_2020" } ,
609  {&ana::khNFSISyst2020_EV2 , "hNFSISyst2020_EV2" } ,
610  {&ana::khNFSISyst2020_EV3 , "hNFSISyst2020_EV3" } ,
611  {ana::GetGenieKnobSyst(rwgt::fReweightMaNCEL) , "MaNCEL" } ,
612  {ana::GetGenieKnobSyst(rwgt::fReweightEtaNCEL) , "EtaNCEL" } ,
613  {&ana::kCOHCCScaleSyst2018 , "COHCCScale2018" } ,
614  {&ana::kCOHNCScaleSyst2018 , "COHNCScale2018" } ,
615  {ana::GetGenieKnobSyst(rwgt::fReweightAhtBY) , "AhtBY" } ,
616  {ana::GetGenieKnobSyst(rwgt::fReweightBhtBY) , "BhtBY" } ,
617  {ana::GetGenieKnobSyst(rwgt::fReweightCV1uBY) , "CV1uBY" } ,
618  {ana::GetGenieKnobSyst(rwgt::fReweightCV2uBY) , "CV2uBY" } ,
619  {ana::GetGenieKnobSyst(rwgt::fReweightAGKY_xF1pi) , "AGKY_xF1pi" } ,
620  {ana::GetGenieKnobSyst(rwgt::fReweightAGKY_pT1pi) , "AGKY_pT1pi" } ,
621  {ana::GetGenieKnobSyst(rwgt::fReweightBR1gamma) , "BR1gamma" } ,
622  {ana::GetGenieKnobSyst(rwgt::fReweightBR1eta) , "BR1eta" } ,
623  {ana::GetGenieKnobSyst(rwgt::fReweightTheta_Delta2Npi) , "Theta_Delta2Npi" } ,
624  // mec for NuX
625  {new ana::MECDoubleGaussEnhSystNux( ana::MECDoubleGaussEnhParamNux::kGauss2DNorm_1Nux , "Norm_1Nux" ), "Norm_1Nux" },
626  {new ana::MECDoubleGaussEnhSystNux( ana::MECDoubleGaussEnhParamNux::kGauss2DMeanQ0_1Nux , "MeanQ0_1Nux" ), "MeanQ0_1Nux" },
627  {new ana::MECDoubleGaussEnhSystNux( ana::MECDoubleGaussEnhParamNux::kGauss2DMeanQ3_1Nux , "MeanQ3_1Nux" ), "MeanQ3_1Nux" },
628  {new ana::MECDoubleGaussEnhSystNux( ana::MECDoubleGaussEnhParamNux::kGauss2DSigmaQ0_1Nux , "SigmaQ0_1Nux" ), "SigmaQ0_1Nux" },
629  {new ana::MECDoubleGaussEnhSystNux( ana::MECDoubleGaussEnhParamNux::kGauss2DSigmaQ3_1Nux , "SigmaQ3_1Nux" ), "SigmaQ3_1Nux" },
630  {new ana::MECDoubleGaussEnhSystNux( ana::MECDoubleGaussEnhParamNux::kGauss2DCorr_1Nux , "Corr_1Nux" ), "Corr_1Nux" },
631  {new ana::MECDoubleGaussEnhSystNux( ana::MECDoubleGaussEnhParamNux::kGauss2DNorm_2Nux , "Norm_2Nux" ), "Norm_2Nux" },
632  {new ana::MECDoubleGaussEnhSystNux( ana::MECDoubleGaussEnhParamNux::kGauss2DMeanQ0_2Nux , "MeanQ0_2Nux" ), "MeanQ0_2Nux" },
633  {new ana::MECDoubleGaussEnhSystNux( ana::MECDoubleGaussEnhParamNux::kGauss2DMeanQ3_2Nux , "MeanQ3_2Nux" ), "MeanQ3_2Nux" },
634  {new ana::MECDoubleGaussEnhSystNux( ana::MECDoubleGaussEnhParamNux::kGauss2DSigmaQ0_2Nux , "SigmaQ0_2Nux" ), "SigmaQ0_2Nux" },
635  {new ana::MECDoubleGaussEnhSystNux( ana::MECDoubleGaussEnhParamNux::kGauss2DSigmaQ3_2Nux , "SigmaQ3_2Nux" ), "SigmaQ3_2Nux" },
636  {new ana::MECDoubleGaussEnhSystNux( ana::MECDoubleGaussEnhParamNux::kGauss2DCorr_2Nux , "Corr_2Nux" ), "Corr_2Nux" },
637  {new ana::MECDoubleGaussEnhSystNux( ana::MECDoubleGaussEnhParamNux::kBaselineNux , "BaselineNux" ), "BaselineNux" }
638  };
639 
640  // add small GENIE PCA scale systs
641  for (int ipca = 0; ipca < 12; ++ipca){
642  std::string gpcaname = "GENIEPCA" + std::to_string(ipca);
643  labelMap[ana::GetGeniePrincipals2020Small(ipca)] = gpcaname;
644  }
645 
646  // now that we have the ISysts for all of the systematics we want to use, get
647  // the -2, -1, +1, +2 shifts and store them in the files.
648  std::string rwgtLabel;
649  std::vector<double> wgt(4, 1.);
650  systShifts.clear();
651  for (auto const& rwgtPair : labelMap){
652 
653  auto wgtFunc = rwgtPair.first;
654  rwgtLabel = rwgtPair.second;
655 
656  wgt[0] = 1.;
657  wgtFunc->Shift(-2., &srProxy, wgt[0]);
658 
659  wgt[1] = 1.;
660  wgtFunc->Shift(-1., &srProxy, wgt[1]);
661 
662  wgt[2] = 1.;
663  wgtFunc->Shift( 1., &srProxy, wgt[2]);
664 
665  wgt[3] = 1.;
666  wgtFunc->Shift( 2., &srProxy, wgt[3]);
667 
668  systShifts.emplace_back(cmf::VarNameToKey(rwgtLabel),
669  (float)wgt[0],
670  (float)wgt[1],
671  (float)wgt[2],
672  (float)wgt[3]);
673 
674  } // end loop over genie reweight variables
675 */
676  }
677 
678  //......................................................................
679  // in this function we save the lepton and hadronic energy separately and
680  // construct the final energy in VarVals
682  NuEvent const& srProxy) const
683  {
684 
685 
686  { //assume reconstructing a NuMu
687  // new implementation
688  dataVals.set_val_at(cmf::kLep_RecoE, (float) srProxy.trkEn);
689  dataVals.set_val_at(cmf::kHad_RecoE, (float) srProxy.shwEnCC);
690  //dataVals.set_val_at(cmf::kRecoQ2, (float) ana::kRecoQ2(&srProxy));
691  //dataVals.set_val_at(cmf::kLep_RecoE_MCFrac, (float)ana::kNumuHadEFrac2020(&srProxy));
692  }
693 
694 
695  }
696 
697  //......................................................................
699  NuEvent const& srProxy) const
700  {
701  evId.run = srProxy.run;
702  evId.subrun = srProxy.subRun;
703  evId.event = srProxy.evt;
704  evId.slice = srProxy.slc;
705  evId.cycle = srProxy.simFlag;//srProxy.hdr.cycle; //TODO: what does cycle mean
706  }
707 
708  //......................................................................
709  void DSTToEventList::FillVariables(std::vector<cmf::DataVarVals> & dataVals,
710  std::vector<cmf::MCVarVals> & mcVals,
711  std::vector<cmf::EventId> & evIds,
712  std::vector<cmf::SelectionType_t> & sels) //TODO: update this method for minos
713  {
714  dataVals.clear();
715  mcVals .clear();
716  evIds .clear();
717 
718  for(auto const& fileName : fCAFFileNames){
719 
720  // open the CAF file and get the TTrees and histograms for reading
721  TFile *tf = TFile::Open(fileName.c_str());
722 
723  this->FindPOTPerPeriod(*tf);
724 
725  auto *cafTree = dynamic_cast<TTree*>(tf->Get(fCAFRecTreeName.c_str()));
726 
727  //TODO add a line for the minos mc tree
728 
729  // hook the TTree up to a standard record and get the first entry to initialize things
730  NuEvent *sr = nullptr;
731 
732  cafTree->SetBranchAddress("s", &sr);
733  cafTree->GetEntry(0);
734 
735  static NuEvent srProxy = NuEvent(); //TODO: is this a sensical thing to do?
736 
737  MF_LOG_DEBUG("DSTToEventList")
738  << "There are "
739  << cafTree->GetEntriesFast()
740  << " caf entries";
741 
742  // now loop over the CAF entries and fill the event list variables
743  //for(int i = 0; i < 10; ++i){
744  for(int i = 0; i < cafTree->GetEntriesFast(); ++i){
745  cafTree->GetEntry(i);
746 
747  // get a proxy for each entry
748  srProxy = *sr;
749 
750  // Uncomment the following block to debug cases where events
751  // show up in the official CAF selection but not this one, or
752  // vice versa
753  /*
754  cmf::EventId curEvent(srProxy.hdr.run,
755  srProxy.hdr.subrun,
756  srProxy.hdr.evt,
757  srProxy.hdr.cycle,
758  srProxy.hdr.subevt);
759 
760  if( fMissingEvents.find(curEvent) != fMissingEvents.end())
761  MF_LOG_VERBATIM("MissingEventsCheck")
762  << curEvent
763  << " beam "
764  << ana::kInBeamSpill(&srProxy)
765  << " spill "
766  << ana::kStandardSpillCuts(&(srProxy.spill))
767  << " FD all "
768  << ana::kNue2020FDAllSamples(&srProxy);
769  */
770 
771  // if ( !this->PassesSelections(srProxy) ) continue; //running over a microDST so all events pass selection
772  sels .emplace_back(cmf::kNuMuSelection);//TODO: for now stick with Numu
773 
774 
775  cmf::TruthVars truthVars;
776  cmf::WeightVars weightVars;
777  cmf::SystVarColl systVars;
778  cmf::DataVarVals dataVarVals;
779  cmf::EventId evId;
780 
781  // TODO: right now we must emplace back the information
782  // immediately after we fill it, else ISyst::Shift, called in
783  // FillMCVals will modify them. This shouldn't happen and really
784  // bothers me, figure that out.
785 
786  // fill the event id information
787  this->FillEventId(evId, srProxy);
788  evIds .emplace_back(evId);
789 
790  // fill the data values
791  this->FillDataVarVals(dataVarVals, srProxy);
792  dataVals.emplace_back(dataVarVals);
793 
794  // fill MC information if this is mc //TODO: addrss MC case
795 /* if(srProxy.hdr.ismc){
796  if(srProxy.mc.nu.size() == 1){
797  this->FillMCVals(truthVars,
798  weightVars,
799  systVars,
800  srProxy);
801  }
802  } // end if the event is MC*/
803 
804 
805  // add everything to the vectors
806  mcVals .emplace_back(truthVars, weightVars, systVars);
807 
808  } // end loop over the CAF StandardRecords
809 
810  tf->Close();
811 
812  } // end loop over CAF file names
813  }
814 
815  //......................................................................
817  {
818  this->InitializeEventListColls();
819  std::cout << "On line " << __LINE__ << std::endl;
822  float pdg = 0;
823 
824  std::vector<cmf::DataVarVals > dataVals;
825  std::vector<cmf::MCVarVals > mcVals;
826  std::vector<cmf::EventId > evIds;
827  std::vector<cmf::SelectionType_t> sels;
828 
829  this->FillVariables(dataVals, mcVals, evIds, sels);
830 
831  MF_LOG_VERBATIM("DSTToEventList")
832  << "there are "
833  << dataVals.size()
834  << " data vals "
835  << mcVals.size()
836  << " mc vals "
837  << evIds.size()
838  << " event ids "
839  << sels.size()
840  << " selections ";
841 
842  auto periodToPOTItr = fPeriodToPOTLiveTime.begin();
843 
844  // need these *Vars for the cosmic muon files
845  cmf::TruthVars tv;
846  cmf::WeightVars wv;
847 
848  // loop over the event lists and add SpillSummary info
849  for(auto & itr : fMDToEvents){
850 
851  auto const& md = itr.ListMetaData();
853 
854  std::cout<< "metadata "
855  << md.ToString()
856  << " has "
857  << itr.ListSpillSummary()
858  << " and "
859  << itr.size()
860  << " events" << std::endl;
861 
862  periodToPOTItr = fPeriodToPOTLiveTime.find(md.Period());
863 
864  // multiply the weight for this period by the total POT
865  // indicated in the list of CAF files used. We only
866  // look at either FHC or RHC SAM datalists separately,
867  // so we don't have to worry about separate accounting
868  // for those
869  if(periodToPOTItr != fPeriodToPOTLiveTime.end()){
870  ss += cmf::SpillSummary(periodToPOTItr->second.pot,
871  periodToPOTItr->second.pot,
872  periodToPOTItr->second.liveTime,
873  1,
874  1);
875 
876  MF_LOG_DEBUG("CAFTToEventList")
877  << " checking spill summary to add "
878  << ss
879  << " "
880  << periodToPOTItr->second.pot
881  << " "
882  << periodToPOTItr->second.liveTime;
883  }
884  else{
885  MF_LOG_WARNING("DSTToEventList")
886  << md.ToString()
887  << " could not find period "
888  << md.Period()
889  << " in POT weight map, just add zeros";
890  }
891 
892  itr.AddSpillSummary(ss);
893 
894  MF_LOG_VERBATIM("DSTToEventList")
895  << "adding "
896  << ss
897  << " to "
898  << md.ToString()
899  << " "
900  << itr.ListSpillSummary();
901 
902  } // end loop over the metadata to event lists map
903 
904  // now make event lists
905  for(size_t i = 0; i < dataVals.size(); ++i){
906 
908 
909  // figure out the interaction type and then get the correct metadata
910  // for it
912  it = cmf::kCosmicMuon;
913  }
914  else if(!fIsRealData){
915  if( mcVals[i].val_at(cmf::kTrueCCNC) == 1. * simb::kNC )
916  it = cmf::kNC;
917  else{
918  pdg = mcVals[i].val_at(cmf::kTruePDG);
919  if (pdg == 12) it = cmf::kNuECC;
920  else if(pdg == -12) it = cmf::kNuEBarCC;
921  else if(pdg == 14) it = cmf::kNuMuCC;
922  else if(pdg == -14) it = cmf::kNuMuBarCC;
923  else if(pdg == 16) it = cmf::kNuTauCC;
924  else if(pdg == -16) it = cmf::kNuTauBarCC;
925  else it = cmf::kUncategorised;
926  }
927  } // end if neutrino MC
928 
930  fDetector,
931  fFileType,
932  sels[i],
933  it,
934  this->RunToPeriod(evIds[i].run));
935 
936  // get an iterator to the event list for this md
937  auto listItr = cmf::FindEventList(md, fMDToEvents);
938  if(listItr == fMDToEvents.end())
939  throw cet::exception("DSTToEventList")
940  << "cannot find metadata "
941  << md.ToString()
942  << " in event list";
943 
945  tv = cmf::TruthVars(mcVals[i].TruthVars());
946  tv.fTrueE = 0. ;
947  tv.fTruePDG = 13 ;
948  tv.fTruePDGOrig = 13 ;
949  tv.fTrueCCNC = 0. ;
951 
952  wv = cmf::WeightVars(mcVals[i].WeightVars());
953  wv.fXSecCVPPFX_Weight = 1. ;
954  wv.fRPACCQE_Weight = 1. ;
955 
956  mcVals[i] = cmf::MCVarVals(tv, wv, mcVals[i].SystematicShifts());
957  listItr->AddEvent(evIds[i],
958  dataVals[i],
959  mcVals[i]);
960  } // end if cosmic backgrounds
961  else if(md.isMC){
962  listItr->AddEvent(evIds[i],
963  dataVals[i],
964  mcVals[i]);
965  }
966  else
967  listItr->AddEvent(evIds[i],
968  dataVals[i]);
969 
970  MF_LOG_DEBUG("DSTToEventList")
971  << "add the event to "
972  << md.ToString()
973  << " "
974  << i;
975 
976  } // end loop to insert events
977 
978 
979  for(auto const& itr : fMDToEvents){
980 
981  auto const& md = itr.ListMetaData();
982 
983  MF_LOG_VERBATIM("DSTToEventList")
984  << "metadata "
985  << md.ToString()
986  << " has "
987  << itr.ListSpillSummary()
988  << " and "
989  << itr.size()
990  << " events";
991 
992  if(fMakeSelectedList){
993  for(auto const& ev : itr){
994  std::stringstream selectionOutput;
995  selectionOutput << cmf::cSelectionType_Strings[md.selectionType]
996  << "\t\t" << ev->EventID().run
997  << "\t" << ev->EventID().subrun
998  << "\t" << ev->EventID().event
999  << "\t" << ev->EventID().slice
1000  << "\t" << ev->EventID().cycle
1001  << "\t" << ev->DataVals().val_at(cmf::kNu_RecoE, md)
1002  << "\t" << ev->DataVals().val_at(cmf::kNuE_CVN, md)
1003  << "\t";
1004  if(md.isMC)
1005  selectionOutput << ev->MCVals().val_at(cmf::kXSecCVPPFX_Weight);
1006 
1007  MF_LOG_VERBATIM("SelectedEvents")
1008  << selectionOutput.str();
1009 
1010  } // end loop over events
1011  } // end if we are making a list of selected events
1012  }
1013 
1014  MF_LOG_DEBUG("DSTToEventList")
1015  << "Serialize the events";
1016 
1018  fMDToEvents);
1019 
1020  MF_LOG_DEBUG("DSTToEventList")
1021  << "endJob() done";
1022 
1023  }
1024 
1026 
1027 } // end cmf namespace
float fTruePDGOrig
true PDG of original neutrino
Definition: VarVals.h:46
POTLiveTime(double const &p, double const &l)
bool RunInPeriod(int run) const
std::string fFHCQuantFileName
name of file in numudata ups product with FHC quantiles
void FillDataVarVals(cmf::DataVarVals &dataVals, NuEvent const &srProxy) const
bool fMakeSelectedList
configure to true to make a list
cmf::DetType_t fDetector
which detector are we using
PeriodRun(int period, int startRun, int endRun)
fileName
Definition: plotROC.py:78
enum cmf::interaction_type InteractionType_t
cmf::FileType_t fFileType
get from the configuration for now
set< int >::iterator it
bool fIsFHC
are we looking at FHC files?
long Period() const
Definition: Structs.h:128
void Initialize(fhicl::ParameterSet const &pset)
float fTrueCCNC
true cc vs nc
Definition: VarVals.h:44
float fTruePDG
true PDG
Definition: VarVals.h:43
double liveTime
Definition: Structs.h:187
void FillEventId(cmf::EventId &evId, NuEvent const &srProxy) const
const char * p
Definition: xmltok.h:285
vector< vector< double > > clear
std::string fPOTHistName
names of the CAF histogram with the total POT
Float_t ss
Definition: plot.C:24
std::string fCAFRecTreeName
names of the CAF tree with the reco info
cmf::EventListColl::iterator FindEventList(cmf::MetaData const &md, cmf::EventListColl &eventListColl)
Definition: Event.cxx:44
enum cmf::det_type DetType_t
::xsd::cxx::tree::exception< char > exception
Definition: Database.h:225
enum cmf::sel_type SelectionType_t
DEFINE_ART_MODULE(TestTMapFile)
Int_t run
fHeader.fRun
Definition: NuEvent.h:37
Timing fit.
std::vector< std::string > fCAFFileNames
name of the CAF file we are using
float fTrueE
True nu energy.
Definition: VarVals.h:42
Base container for the MC related Vars that constitute an event.
Definition: VarVals.h:227
const XML_Char * s
Definition: expat.h:262
static ParameterUtility * Instance()
int subrun
Definition: VarVals.h:80
float fRPACCQE_Weight
Deprecated.
Definition: VarVals.h:151
cmf::SelectionType_t selectionType
Definition: Structs.h:116
const Cut sels[kNumSels]
Definition: vars.h:44
int cycle
Definition: VarVals.h:82
T get(std::string const &key) const
Definition: ParameterSet.h:231
void set_val_at(uint8_t const &varkey, float const &val)
Definition: VarVals.cxx:178
cmf::SelectionType_t FindSelectionType(NuEvent const &srProxy) const
static cmf::DetType_t StringToDetectorType(std::string const &str)
Definition: StaticFuncs.h:267
Float_t shwEnCC
Reconstructed shower energy assuming this is a CC event.
Definition: NuEvent.h:175
caf::StandardRecord * sr
const std::vector< std::string > cSelectionType_Strings({"NuESel_AllPID","NuESel_LowPID","NuESel_MidPID","NuESel_HighPID","NuESel_Peripheral","NuMuSel","NuMuSelQ1","NuMuSelQ2","NuMuSelQ3","NuMuSelQ4","NCSel","UnknownSel"})
Int_t evt
reco&#39;d event number, evt.index
Definition: NuEvent.h:210
Float_t trkEn
Best reconstructed track energy.
Definition: NuEvent.h:149
int slice
Definition: VarVals.h:83
void reconfigure(fhicl::ParameterSet const &p)
Base container for the MC related Vars that constitute an event.
Definition: VarVals.h:268
cmf::EventListColl fMDToEvents
each entry in the vector is for a single slice
float fXSecCVPPFX_Weight
Deprecated, Was Tufts weight for SA.
Definition: VarVals.h:150
Definition: run.py:1
void analyze(art::Event const &e) override
std::vector< cmf::SystVar > SystVarColl
Definition: VarVals.h:223
void SerializeEventListColl(std::string const &dirName, cmf::EventListColl const &listColl)
Definition: Event.cxx:117
OStream cout
Definition: OStream.cxx:6
static float min(const float a, const float b, const float c)
Definition: absgeo.cxx:45
void FillVariables(std::vector< cmf::DataVarVals > &dataVals, std::vector< cmf::MCVarVals > &mcVals, std::vector< cmf::EventId > &evIds, std::vector< cmf::SelectionType_t > &sels)
std::string ToString() const
Definition: Structs.cxx:135
void FillMCVals(cmf::TruthVars &truthVars, cmf::WeightVars &weightVars, cmf::SystVarColl &systVars, NuMCEvent &srProxy) const
Module to combine a set of results into a single file currently only does one data product type at a ...
Definition: Event.cxx:24
std::vector< cmf::EventList > EventListColl
Definition: Event.h:150
unsigned int fNumQuantiles
number of quantiles
Int_t runPeriod
Major analysis run period, e.g. "Run IV" == 4.
Definition: NuEvent.h:1444
#define MF_LOG_VERBATIM(category)
Int_t subRun
fHeader.fSubRun
Definition: NuEvent.h:38
#define MF_LOG_DEBUG(id)
DSTToEventList(fhicl::ParameterSet const &pset)
std::set< cmf::SelectionType_t > fSelections
Selections to use.
std::vector< cmf::PeriodRun > fPeriodRuns
vector of periods to runs
Int_t slc
reco&#39;d slice number, evt.slc
Definition: NuEvent.h:211
void operator+=(POTLiveTime const &o)
#define MF_LOG_WARNING(category)
Float_t e
Definition: plot.C:35
float fTrueIntType
true interaction type
Definition: VarVals.h:45
std::string fRHCQuantFileName
name of file in numudata ups product with FHC quantiles
T max(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
Definition: statistics.h:66
int event
Definition: VarVals.h:81
std::map< int, POTLiveTime > fPeriodToPOTLiveTime
keep track of the POT in each period
Int_t simFlag
fHeader.fVldContext.fSimFlag
Definition: NuEvent.h:47
detName
Definition: mkDefs.py:107
bool fIsRealData
are we looking at real data?
std::string fSystType
systematic type
std::string fSpillTreeName
names of the Spill tree in the file
enum cmf::file_type FileType_t
enum BeamMode string