EventListManipulator.cxx
Go to the documentation of this file.
1 //
2 // EventListManipulator.cxx
3 //
4 // Created by Brian Rebel on 3/21/16.
5 //
6 
7 #include <stdexcept>
8 
9 #include "TFile.h"
10 #include "TH2.h"
11 #include "TH3.h"
12 
14 #include "FNEX/core/Event.h"
15 #include "FNEX/core/Correlations.h"
18 #include "FNEX/core/VarVals.h"
19 #include "FNEX/core/VarBinning.h"
20 #include "SummaryData/POTSum.h"
21 
23 
24 namespace fnex{
25 
26  //----------------------------------------------------------------------------
28  fnex::InputPoint const& initialGuess)
29  : fInitialGuess(initialGuess)
30  , fTotalFDDataPOTFHC(0.)
31  , fTotalFDDataPOTRHC(0.)
32  , fTotalNDDataPOTFHC(0.)
33  , fTotalNDDataPOTRHC(0.)
34  {
35  fEventCaps.clear();
36 
37  this->reconfigure(p);
38 
39  return;
40  }
41 
42  //----------------------------------------------------------------------------
44  {
45  return;
46  }
47 
48  //......................................................................
50  {
51  // This call sets the seed to a random value
52  fRandom.SetSeed(pset.get<unsigned int>("RandomSeed", 0));
53 
54  fFNEXEventLabels = pset.get< std::vector<std::string> >("FNEXEventLabels", std::vector<std::string>() );
55  fTreeDirectories = pset.get< std::vector<std::string> >("TreeDirectories", std::vector<std::string>() );
56  fPoissonFakeData = pset.get< bool >("PoissonFakeData", false );
57  fUseCosmicHists = pset.get< bool >("UseCosmicHists", true );
58  fUseEventId = pset.get< bool >("UseEventId", false );
59  fSumSubEpochs = pset.get< bool >("SumSubEpochs", false );
60  fNuMuCosmicHistFile = pset.get<std::string >("NuMuCosmicBackgroundHist", "CAFAna/numu/cosmics_fhc__numu2018.root" );
61  fMaxNuEnergy[fnex::kNuESelection] = pset.get< double >("MaxNuEEnergyForAnalysis", 5.0 );
62  fMaxNuEnergy[fnex::kNuMuSelection] = pset.get< double >("MaxNuMuEnergyForAnalysis", 10.0 );
63  auto fakePars = pset.get< std::vector<fhicl::ParameterSet> >("FakeParameters", std::vector<fhicl::ParameterSet>() );
64  auto epochPOTPars = pset.get< std::vector<fhicl::ParameterSet> >("EpochsToUse" );
65  auto highStatsFake = pset.get< bool >("HighStatsFake", false );
66  auto selections = pset.get< std::vector<fhicl::ParameterSet> >("SelectionsToUse", std::vector<fhicl::ParameterSet>() );
67  fAnalysisYear = pset.get< std::string >("AnalysisYear", "2016" );
68  fMakeNDPeripheral = pset.get< bool >("MakeNDPeripheral", false );
69  fUseExtrapolation = pset.get< bool >("UseExtrapolation", true );
70  fMaxEventsPerTree = pset.get< float >("MaxEventsPerTree", 1.e6);
71  if(fAnalysisYear == "2017"){
72  fNuECosmicHistFile= pset.get<std::string >("NuECosmicBackgroundHist", "CAFAna/nue/cosmic_spectra_prediction_2017.root" );
73  fNuEDataHistFile = pset.get<std::string >("NuEDataHistFile", "CAFAna/nue/data_spectra_fd_2017.root" );
74  fCalibHistFile = pset.get<std::string >("CalibSystFile", "FNEX/data/calibration/CalibrationSystematics2017.root" );
75  }//end of 2017 analysis year condition
76  else if(fAnalysisYear == "2018"){
77  fNuECosmicHistFile= pset.get<std::string >("NuECosmicBackgroundHist", "CAFAna/nue/cosmic_spectra_prediction_2018.root" );
78  fCalibHistFile = pset.get<std::string >("CalibSystFile", "FNEX/data/calibration/CalibrationSystematics2018.root");
79  //fNuEDataHistFile = pset.get<std::string >("NuEDataHistFile", "CAFAna/nue/data_spectra_fd_2017.root" );
80  }//end of 2018 analysis year condition
81  else
82  fNuECosmicHistFile= pset.get<std::string >("NuECosmicBackgroundHist", "CAFAna/nue/NueSA_numi_data.root" );
83 
84  for(size_t s = 0; s < fTreeDirectories.size(); ++s) fTreeDirectories[s].append("/full");
85 
86  // determine which general experiments we have, ie nue or numu or both
87  // allow users to force the use of all event lists if they choose
88  fLoadAllEventLists = pset.get<bool>("LoadAllEventLists", false);
89  this->DetermineSelections(selections);
90 
91  std::string epoch;
93  double exposure;
94  for(auto itr : epochPOTPars){
95  period = itr.get<std::string>("Period");
96  epoch = itr.get<std::string>("Epoch" );
97  exposure = itr.get<double >("Exposure");
98 
99  fPeriodsToUse.insert(std::atoi(period.c_str()) * 1000);
100  fFDPOTSum.emplace(period+epoch, fnex::SpillSummary(exposure,
101  exposure,
102  0.,
103  1,
104  1));
105 
106  LOG_VERBATIM("EventListManipulator")
107  << "Using "
108  << period + epoch
109  << " data for this analysis";
110 
111  // now add the exposure to the total period too
112  fFDPOTSum[period] += fnex::SpillSummary(exposure,
113  exposure,
114  0.,
115  1,
116  1);
117 
118  // TODO: We need to find a way to get set the total FD Data POT automatically
119  // if(period.find("1") != std::string::npos ||
120  // period.find("2") != std::string::npos ||
121  // period.find("3") != std::string::npos ||
122  // period.find("5") != std::string::npos)
123  // fTotalFDDataPOTFHC += exposure;
124 
125  // else if(period.find("4") != std::string::npos ||
126  // period.find("6") != std::string::npos)
127  // fTotalFDDataPOTRHC += exposure;
128 
129  LOG_VERBATIM("EventListManipulator")
130  << "Epoch : "
131  << period+epoch
132  << " Exposure : "
133  << exposure
134  << " Accumulated POT FHC : "
136  << " Accumulated POT RHC ; "
138 
139  }
140 
141  // limits for the fraction of events to use from each MD type
142  // the configuration uses a set fraction for each MD type initially,
143  // but this number can be adjusted when the caps are applied if a given
144  // MD type has too few events that way. Too few is hardcoded to be 2.e4
145  auto caps = pset.get<fhicl::ParameterSet>("EventCap", fhicl::ParameterSet());
146  std::string comboStr;
147  std::string keyStr;
148  for(auto itr : caps.get<std::vector<fhicl::ParameterSet> >("Caps", std::vector<fhicl::ParameterSet>())){
149  // get the strings
150  comboStr = itr.get<std::string>("Detector");
151  comboStr += "P";
152  comboStr += itr.get<std::string>("Period" );
153  comboStr += itr.get<std::string>("FileType");
154 
155  // loop over the selection and interaction types to flesh out the
156  for(int s = 0; s < fnex::kUnknownSelection; ++s){
157 
158  if(itr.get<std::string>("FileType").compare("Data") != 0){
159  for(int i = 1; i < fnex::kCosmicMuon + 1; ++i){
160 
162  LOG_DEBUG("EventListManipulator")
163  << "cap for "
164  << keyStr
165  << " : "
166  << fnex::StringToKey(keyStr)
167  << " is "
168  << itr.get<float>("Cap");
169 
170  fEventCaps[fnex::StringToKey(keyStr)] = itr.get<float>("Cap");
171 
172  } // end loop over interaction types
173  } // end if not data
174  else{
175  keyStr = comboStr + fnex::cSelectionType_Strings[s];
176 
177  LOG_DEBUG("EventListManipulator")
178  << "cap for "
179  << keyStr
180  << " : "
181  << fnex::StringToKey(keyStr)
182  << " is "
183  << itr.get<float>("Cap");
184 
185  fEventCaps[fnex::StringToKey(keyStr)] = itr.get<float>("Cap");
186  }
187 
188  } // end loop over selection types
189  } // end loop over event cap parameter sets
190 
191  // for(auto const& itr : fEventCaps)
192  // LOG_VERBATIM("EventListManipulator")
193  // << "Event caps key: "
194  // << itr.first
195  // << " cap: "
196  // << itr.second;
197 
198  LOG_DEBUG("EventListManipulator")
199  << "There are "
200  << fakePars.size()
201  << " fake parameter sets";
202 
203  fFakeExposure.clear();
204  for(auto fake : fakePars){
205 
206  auto fakeDet = fake.get<std::string >("Detector" );
207  auto fakeEpoch = fake.get<std::string >("Epoch" );
208  auto fakeExposure = fake.get<float >("Exposure" );
209 
210  // if we are summing the sub epochs, remove the alpha bit of the epoch string
211  if(fSumSubEpochs){
212  size_t first = fakeEpoch.find_first_of("abcdefghijklmnopqrstuvwxyz");
213  if(first != std::string::npos)
214  fakeEpoch.erase(first, 1);
215  }
216 
218  if (fakeDet.compare("FarDet") == 0) det = novadaq::cnv::kFARDET;
219  else if(fakeDet.compare("NearDet") == 0) det = novadaq::cnv::kNEARDET;
220 
221  // if we want high stats fakes, multiply the exposure by 1e3
222  // we only ever apply high stats to the far detector
223  if(highStatsFake && det == novadaq::cnv::kFARDET) fakeExposure *= 1.e3;
224 
225  // Use the selections from above to determine which fake selections we want
226  // It is very unlikely that we would use some real and some fake in the same detector
227  for(auto const &sel : fSelectionsToUse[det]){
228  fnex::MetaData md(false,
229  det,
231  sel,
233  fakeEpoch);
234 
235  LOG_DEBUG("EventListManipulator")
236  << "configuring fake parameter set "
237  << std::setw(30)
238  << md.ToString()
239  << " to have "
240  << fakeExposure
241  << " POT";
242 
243  if(fFakeExposure.count(md) < 1)
244  fFakeExposure.emplace(md, fakeExposure);
245  else
246  fFakeExposure[md] += fakeExposure;
247 
248  } // end loop over selections
249 
250  } // end if we have fake detectors and selections to make
251 
252  // Set the fake data point
253  if(fFakeExposure.size() > 0){
254  fnex::InputPoint ip(pset.get< fhicl::ParameterSet>("FakeInputPoint"), true);
255 
256  // Set the fake data correlations
257  fCorrelations = std::make_unique<Correlations>(pset.get<fhicl::ParameterSet>("FakeCorrelations"));
258 
260 
261  // Apply correlations to the fake data point
262  fCorrelations->ApplyCorrelations(fFakeDataPoint);
263  }
264 
265  return;
266  }
267 
268  //......................................................................
269  void EventListManipulator::DetermineSelections(std::vector<fhicl::ParameterSet> const& selections)
270  {
271 
272  fSelectionsToUse.clear();
273 
274  std::string detStr;
275  std::string selStr;
277  std::set<fnex::SelectionType_t> selSet;
278 
279  // the parameter sets also contain the beam type (ie FHC or RHC)
280  // for now we are ignoring that, but it is available here should we want
281  // it.
282  // The parameter sets will attempt to add the same detector and selection
283  // combinations to the map twice, but maps can only have unique keys so
284  // there will be no duplication
285 
286  // std::string beamStr;
287  // fnex::BeamType_t beam;
288 
289  for(auto const& itr : selections){
290 
291  // beamStr = itr.get<std::string>("Beam");
292  // beam = (beamStr.find("FHC") != std::string::npos) ? fnex::kFHC : fnex::kRHC;
293 
294  detStr = itr.get<std::string>("Detector");
295  det = (detStr.find("Far") != std::string::npos) ? novadaq::cnv::kFARDET : novadaq::cnv::kNEARDET;
296 
297  selSet.clear();
298  auto selStrs = itr.get<std::vector<std::string>>("Selections");
299  for(auto& sel: selStrs){
300  if (sel.find("NuESel_Peripheral") != std::string::npos) selSet.insert(fnex::kNuESelectionPeripheral);
301  else if(sel.find("NuESel_LowPID") != std::string::npos) selSet.insert(fnex::kNuESelectionLowPID );
302  else if(sel.find("NuESel_HighPID") != std::string::npos) selSet.insert(fnex::kNuESelectionHighPID );
303  else if(sel.find("NuMuSelQ1") != std::string::npos) selSet.insert(fnex::kNuMuSelectionQ1 );
304  else if(sel.find("NuMuSelQ2") != std::string::npos) selSet.insert(fnex::kNuMuSelectionQ2 );
305  else if(sel.find("NuMuSelQ3") != std::string::npos) selSet.insert(fnex::kNuMuSelectionQ3 );
306  else if(sel.find("NuMuSelQ4") != std::string::npos) selSet.insert(fnex::kNuMuSelectionQ4 );
307  } // end loop over selection strings
308 
309  fSelectionsToUse.emplace(det, selSet);
310 
311  } // end loop over detectors
312 
313  return;
314  }
315 
316  // Skip any event list where the selection and detector combination were
317  // not requested in the configuration
318  //......................................................................
320  fnex::DataMC_t const dataMC)
321  {
322  if(fLoadAllEventLists) return true;
323 
324  // check if this is data and we only want MC events
325  // or vice versa
326  if(dataMC == fnex::kMC && !md.isMC){
327  LOG_DEBUG("EventListManipulator")
328  << md.ToString()
329  << " failed fnex::kMC and !md.isMC";
330  return false;
331  }
332  else if(dataMC == fnex::kData && md.isMC){
333  LOG_DEBUG("EventListManipulator")
334  << md.ToString()
335  << " failed fnex::kData and md.isMC";
336  return false;
337  }
338  else if(dataMC == fnex::kFakeData && !md.isMC){
339  LOG_DEBUG("EventListManipulator")
340  << md.ToString()
341  << " failed fnex::kFakeData and !md.isMC";
342  return false;
343  }
344 
345  //Skip unless the epoch is in the fPeriodsToUse set
346  if(fPeriodsToUse.count(md.Period()) < 1){
347  LOG_VERBATIM("EventListManipulator")
348  << md.ToString()
349  << " failed count of periods to use "
350  << md.Period() * 1e-3;
351  return false;
352  }
353  // check if the detector for this metadata is desired, return false if it
354  // isn't in the map
355  if(fSelectionsToUse.count(md.detector) < 1){
356  LOG_VERBATIM("EventListManipulator")
357  << md.ToString()
358  << " failed detector "
359  << md.detector
360  << " in selections to use";
361  return false;
362  }
363  // check if the selection is one we want
364  auto const& selSet = fSelectionsToUse.find(md.detector)->second;
365  if(selSet.count(md.selectionType) < 1 ){
366  LOG_DEBUG("EventListManipulator")
367  << md.ToString()
368  << " failed want "
369  << md.DetectorString()
370  << " "
371  << md.detector
372  << " "
374 
375  return false;
376  }
377  // Metadata types we are going to skip
378  // - ND flux swap for now
379  //if(md->fileType == fnex::kSwap &&
380  // md->detector == novadaq::cnv::kNEARDET ) return false;
381 
382  return true;
383  }
384 
385  //......................................................................
386  // We want to make sure that the pre-configured cap does not result in
387  // a poor sampling of events, so pass in the number of events in the
388  // tree as well to check how the cap will affect the final number of events
389  // used. If too few events will remain after that cap, up it to ensure
390  //
392  int treeEvents)
393  {
394  float cap = 1.;
395 
396  auto key = md.Key();
397 
398  float testNum = 0.;
399  if(fEventCaps.count(key) > 0){
400  // Capping events from trees with very few events to begin with
401  // doesn't make a lot of sense. Make sure we have at least 20k
402  // events in a tree before capping. The statistical uncertainty on
403  // 20k events is < 1%
404  // if we are going to have more than fMaxEventsPerTree in a tree, we might
405  // want to cap it at fMaxEventsPerTree to keep the memory footprint down
406  testNum = treeEvents * fEventCaps[key];
407  if(testNum < 2e4){
408  if(treeEvents > 2.e4) fEventCaps[key] = 2.e4 / (1. * treeEvents);
409  else fEventCaps[key] = 1.;
410  }
411  else if(testNum > fMaxEventsPerTree) fEventCaps[key] = fMaxEventsPerTree / (1. * treeEvents);
412 
413  cap = fEventCaps[key];
414  }
415 
416  // if no cap is specified, use all the events
417  return cap;
418  }
419 
420  //......................................................................
422  std::string const& dirName,
423  fnex::DataMC_t const& dataMC,
424  std::map<fnex::MetaData, fnex::SpillSummary> & mdToSS)
425  {
426  LOG_DEBUG("EventListManipulator")
427  << "extracting metadata from file";
428 
429  TTree *metadataTree = dynamic_cast<TTree *>(metadataFile->Get((dirName + "/metadata").c_str()));
430 
431  if(!metadataTree){
432  throw cet::exception("EventListManipulator")
433  << "cannot find metadata tree: "
434  << dirName + "/metadata";
435  }
436 
437  // we'll want the cosmic background and FD data metadata later on,
438  // so make sets for those metadata here
439  std::set<fnex::MetaData> fdMD;
440 
441  // loop over the entries in the metadataTree
442  fnex::MetaData *md = nullptr;
443  fnex::SpillSummary *ss = nullptr;
444  metadataTree->SetBranchAddress("metadata", &md);
445  metadataTree->SetBranchAddress("spillsummary", &ss);
446 
447  // We have potentially hadd'ed several files together to get the
448  // eventlisttree, so first loop over the entries in the metadata
449  // tree and sum all the POT values for each metadata type
450  for(int ctr = 0; ctr < metadataTree->GetEntries(); ++ctr){
451  metadataTree->GetEntry(ctr);
452 
453  // there are several different MetaData in the FD for
454  // each of these file types, but only the first one
455  // for each period and beam type will be put into
456  // the set, which should work just fine for us.
457  if(md->detector == novadaq::cnv::kFARDET &&
458  md->fileType == fnex::kDataFile &&
460  fdMD.insert(*md);
461 
462  LOG_DEBUG("EventListManipulator")
463  << md->detector
464  << " "
465  << md->fileType
466  << " "
467  << md->selectionType
468  << " "
469  << ss->goodPOT
470  << " "
471  << ss->liveTime
472  << " "
473  << fdMD.size();
474  }
475 
476  // if summing the epochs in a period together, set the md.epoch value
477  // to be the period
478  if(fSumSubEpochs) (*md).epoch = md->Period();
479 
480  // The MC POT information is in units of POT, while for data it is
481  // in units of 10^12 POT, renormalize the MC
482  //
483  // Note that for FD data we do not add the spill summary
484  // information because we are just using the pot info from
485  // an input configuration, not what is stored in the metadata tree
486  if(md->isMC){
487  ss->goodPOT *= 1.e-12;
488  ss->totalPOT *= 1.e-12;
489  }
490  else{
491  if(md->detector == novadaq::cnv::kFARDET){
492 
493  // ss->goodPOT = fFDPOTSum.find(md->EpochString())->second.goodPOT;
494  // ss->totalPOT = ss->goodPOT;
495 
496  if (md->BeamType() == fnex::kFHC) fTotalFDDataPOTFHC += ss->goodPOT;
497  else if(md->BeamType() == fnex::kRHC) fTotalFDDataPOTRHC += ss->goodPOT;
498  } // end if FD data
499  else if(md->detector == novadaq::cnv::kNEARDET){
501  if(md->BeamType() == fnex::kFHC) fTotalNDDataPOTFHC += ss->goodPOT;
502  else if(md->BeamType() == fnex::kRHC) fTotalNDDataPOTRHC += ss->goodPOT;
503  } // end test for numu selection in ND
504 
505  } // end if ND data
506  } // end looking at data
507 
508  mdToSS[*md] += *ss;
509  } // end loop over the metadata trees
510 
511  /*
512  for(auto const& itr : mdToSS){
513  LOG_VERBATIM("EventListManipulator")
514  << "checking live time "
515  << itr.first.ToString()
516  << " "
517  << itr.second.liveTime;
518  }
519  */
520 
521  LOG_DEBUG("EventListManipulator")
522  << "There are "
523  << fdMD.size()
524  << " fd data metadata types";
525 
526  // loop over the FD metadata and set the POT for the FD cosmic ray background
527  // to be (data POT) x (cosmic ray live time)/(data live time)
528  double crPOTEquiv = 0.;
529  auto mdToSSItr = mdToSS.begin();
530  for(auto fdItr : fdMD){
531 
532  mdToSSItr = mdToSS.find(fdItr);
533 
534  crPOTEquiv = mdToSSItr->second.goodPOT / mdToSSItr->second.liveTime;
535 
536  LOG_DEBUG("EventListManipulator")
537  << "FD data type: "
538  << fdItr.ToString()
539  << " has cosmic ray POT equiv scale: "
540  << crPOTEquiv
541  << " "
542  << mdToSSItr->second.goodPOT
543  << " "
544  << mdToSSItr->second.liveTime;
545 
546  // now loop over the the mdToSS to find the cosmic ray metadata with
547  // the same detector and period and apply the scaling to the POT values
548  for(auto & itr : mdToSS){
549 
550  if(itr.first.detector == fdItr.detector &&
551  itr.first.epoch == fdItr.epoch &&
552  itr.first.fileType == fnex::kCosmicBackgroundFile){
553  itr.second.goodPOT = crPOTEquiv * itr.second.liveTime;
554  itr.second.totalPOT = crPOTEquiv * itr.second.liveTime;
555 
556  LOG_DEBUG("EventListManipulator")
557  << itr.first.ToString()
558  << " live time "
559  << itr.second.liveTime
560  << " pot "
561  << itr.second.goodPOT;
562  }
563 
564  } // end loop over mdToSS to find the cosmic ray metadata entries
565 
566  } // end loop over FD data metadata
567 
568  // loop over the mdToSS map one more time to remove any
569  // metadata that are not desired. We waited until now
570  // to do this to ensure that we had all the necessary
571  // metadata to fix the POT equivalent exposure for the
572  // cosmic ray backgrounds
573  // Check that we need this metadata class based on the desired
574  // selections for each detector, epochs and data vs mc
575  for(auto const& itr : mdToSS){
576  if(!this->UseEventsFromMetaData(itr.first, dataMC)) mdToSS.erase(itr.first);
577  }
578 
579  return;
580  }
581 
582  //......................................................................
584  std::string const& dirName,
585  fnex::DataMC_t const& dataMC,
586  std::pair<fnex::MetaData, fnex::SpillSummary> const& mdToSS,
587  fnex::EventListMap & eventLists)
588  {
589  // find the event tree using the metadata information
590  std::string eventTreeName(dirName + "/");
591  auto const& mdString = mdToSS.first.ToString();
592 
593  // need a non-const spill summary for filling the event list in case we have to
594  // apply a cap
595  fnex::SpillSummary ss = mdToSS.second;
596 
597  TTree *eventTree = dynamic_cast<TTree*>(eventFile->Get((eventTreeName + mdString).c_str()));
598 
599  if( !eventTree ){
600  throw cet::exception("EventListManipulator")
601  << "cannot find event tree: "
602  << eventTreeName + mdString;
603  } // end if we didn't get a tree
604 
605  eventLists.emplace(mdToSS.first, this->FillEventList(eventTree, mdToSS.first, ss));
606 
607  return;
608  }
609 
610  //......................................................................
611  void EventListManipulator::CheckDataMCListConsistency(std::set<fnex::MetaDataLite> const& dataSets,
612  std::set<fnex::MetaDataLite> const& mcSets,
613  fnex::DataMC_t const& dataMC,
614  fnex::EventListMap & eventLists)
615  {
616  // nothing to check if we are using fake data, it is made
617  // from the MC sets, but may not exist at the point this
618  // function is called.
619  if(dataMC == fnex::kFakeData) return;
620 
621  // throw an exception if there are data MetaData not in the MC,
622  // but add empty data sets to the list if there are MC epochs not
623  // in the data. That can happen in the FD
624  if(dataMC != fnex::kData){
625  for(auto itr : dataSets){
626 
627  if(mcSets.count(itr) < 1)
628  throw cet::exception("EventListManipulator")
629  << "data set "
630  << itr.ToString()
631  << " is in the data but not MC";
632  }
633  } // end if we want MC events
634  if(dataMC != fnex::kMC){
635  for(auto itr : mcSets){
636 
637  if(dataSets.count(itr) < 1){
638  LOG_VERBATIM("EventListManipulator")
639  << "data set "
640  << itr.ToString()
641  << " is in the MC but not data, create empty data list";
642 
643  fnex::MetaData md(false,
644  itr.detector,
646  itr.selectionType,
648  itr.EpochString());
649 
650  eventLists.emplace(md,
652  fFDPOTSum.find(itr.EpochString())->second)
653  );
654  }
655  }
656  } // end if we want data events
657 
658  return;
659  }
660 
661  //......................................................................
663  {
664  auto const &selSet = fSelectionsToUse.find(novadaq::cnv::kFARDET)->second;
665  std::set <std::string> periods({"1", "2", "3", "4", "5", "6"});
666  for(auto const &sel : selSet){
667  for(auto const &per : periods){
668  fnex::MetaData md(true,
671  sel,
673  per);
674 
675  // need to get the FD data exposure for this selection
676  auto const &ss = fFDPOTSum[per];
677 
678  auto cosmicList = this->CreateFarDetCosmicBackgroundsFromHist(md, ss);
679 
680  eventLists.emplace(md, cosmicList);
681  } // end loop over periods
682  } // end loop over selections
683 
684  return;
685  }
686 
687  //......................................................................
688  // Use this method to load the TTrees containing the fnex::EventLists
689  // until we are able to store them directly in an art file
690  //
691  // If onlyMC is true, only load the MC lists from the file
693  {
694  // the fFNEXEventLabels are assumed to hold the full path to the
695  // files containing the TTrees from which we load the different sample sets
696  fnex::EventListMap eventLists;
697  std::set<fnex::MetaDataLite> dataSets;
698  std::set<fnex::MetaDataLite> mcSets;
699 
700  size_t totalEvents = 0;
701 
702  std::string mdString;
703 
704  fnex::MetaData sumMD;
705 
706  // each input file is assumed to have a metadata tree and a fnexEvent tree
707  for(auto & fileName : fFNEXEventLabels){
708 
709  TFile *f = TFile::Open(fileName.c_str());
710 
711  if(f->IsOpen())
712  LOG_VERBATIM("EventListManipulator")
713  << "successfully opened "
714  << fileName;
715  else
716  throw cet::exception("EventListManipulator")
717  << "Cannot open file "
718  << fileName;
719 
720  std::map<fnex::MetaData, fnex::SpillSummary> mdToSS;
721 
722  for(auto const& dirName : fTreeDirectories){
723 
724  this->ExtractMetaDataFromFile(f, dirName, dataMC, mdToSS);
725 
726  LOG_DEBUG("EventListManipulator")
727  << "There are "
728  << mdToSS.size()
729  << " metadata categories";
730 
731  // now get the event information by looping over the collection of
732  // metadata and spill summary info. This way we have the correct POT
733  // accounting for each metadata type
734 
735  for(auto const& itr : mdToSS){
736 
737  // if(!itr.first.isMC)
738  // LOG_VERBATIM("EventListManipulator")
739  // << itr.first.ToString()
740  // << " has "
741  // << itr.second.goodPOT
742  // << " x 10^20 POT";
743 
744  // create lists of data and MC metadata to check for consistency later
745  if(itr.first.isMC) mcSets .insert(itr.first.LiteMD());
746  else dataSets.insert(itr.first.LiteMD());
747 
748  this->ExtractEventsFromFile(f, dirName, dataMC, itr, eventLists);
749 
750  } // end loop over metadata and spill summaries
751 
752  //Let's make the Peripheral eventlists for the ND
753  //It is a temporary hack to use existing event list
754  if(fMakeNDPeripheral) this->CreateNDPeripheralLists(eventLists);
755 
756  } // end loop over directory names
757 
758  f->Close();
759 
760  } // end loop over the input file names
761 
762  // check that the data and MC epochs match
763  this->CheckDataMCListConsistency(dataSets, mcSets, dataMC, eventLists);
764 
765  // now add in the FD cosmic ray background for each selection
766  // desired and beam type, unless we only want data
767  //if(dataMC != fnex::kData) this->CreateFarDetCosmicBackgrounds(eventLists);
768 
769  for(auto const& itr : eventLists){
770  LOG_DEBUG("EventListManipulator")
771  << itr.first.ToString();
772  totalEvents += itr.second.size();
773  }
774 
775  LOG_VERBATIM("EventListManipulator")
776  << "There are "
777  << eventLists.size()
778  << " event lists and "
779  << totalEvents
780  << " events";
781 
782 
783  if(dataMC != fnex::kData){
784  //Do this before making histograms so weighters will be initialised
785  if(fFakeExposure.size() > 0) this->FillFakeDataLists(eventLists);
786 
787  // We have already made the fake data so we should be good to set the shifter
788  // and weighter to use the initial guess point to fill the histograms for the
789  // MC lists
792  //this->MakeHistograms(eventLists, "SkimmedLists");
793  }
794 
795  for(auto const& itr : eventLists)
796  LOG_VERBATIM("EventListManipulator")
797  << std::setw(55)
798  << itr.first.ToString()
799  << " of period "
800  << itr.first.Period()
801  << " "
802  << std::setw(15)
803  << " represents "
804  << itr.second.ListSpillSummary().goodPOT*1.e-8
805  << " x 10^20 POT and "
806  << itr.second.size()
807  << " events "
808  << itr.second.size() / itr.second.ListSpillSummary().goodPOT;
809 
810  return eventLists; // we relinquish ownership.
811  }
812 
813  //......................................................................
815  {
816  // to make fake data we need to deserialize both data and MC
817  if(fFakeExposure.size() > 0) return this->Deserialize(fnex::kBoth);
818 
819  return this->Deserialize(fnex::kData);
820  }
821 
822  //......................................................................
824  bool isMC,
825  std::vector<std::string> const& indepVarNames,
826  std::vector<float> const& indepVarVals)
827  {
828  if (indepVarNames.size() != indepVarVals.size())
829  throw cet::exception("EventListManipulator")
830  << "IndepVar names and IndepVar values not of same length for event";
831 
832  std::string varName;
833 
834  //unsigned char sigmaKey;
835  size_t pos = 0;
836  std::string sigmaStr;
837  //bool isMCVar = false;
838 
839  for(size_t varNum = 0; varNum < indepVarNames.size(); ++varNum){
840  varName = indepVarNames[varNum];
841 
842  if(varName.compare("Nue_RecoE") == 0 ||
843  varName.compare("NuMu_RecoE") == 0 ){
844  LOG_WARNING("EventListManipulator")
845  << "Variable name "
846  << varName
847  << " is deprecated and should not be loaded, skip it";
848  continue;
849  }
850 
851  // first check to see if we have a MC related variable
852  pos = varName.find("sigma");
853  sigmaStr = "Unknown";
854  if( pos != std::string::npos){
855  sigmaStr = varName.substr(pos-2, 7);
856  varName.erase(pos-2, 7);
857  }
858 
859  if(fnex::IsMCVar(varName) && isMC){
860  // check if the variable name contains "sigma", if it does, we need
861  // to figure out which one we are looking at and act accoringly
862 
863  LOG_DEBUG("EventListManipulator")
864  << "loading variable "
865  << varName
866  << " with value "
867  << indepVarVals[varNum];
868 
869  if(fnexEvent.mcVals.get() == nullptr) fnexEvent.mcVals = std::make_unique<fnex::MCVarVals>();
870  if(pos == std::string::npos)
871  fnexEvent.mcVals->set_val_at(varName, indepVarVals[varNum]);
872  else
873  fnexEvent.mcVals->set_GENIEWgt(varName, sigmaStr, indepVarVals[varNum]);
874 
875  LOG_DEBUG("EventListManipulator")
876  << "just loaded MC vals are \n"
877  << *(fnexEvent.mcVals);
878  }
879  else if(fnex::IsDataVar(varName)){
880  if(fnexEvent.dataVals.get() == nullptr) fnexEvent.dataVals = std::make_unique<fnex::DataVarVals>();
881  fnexEvent.dataVals->set_val_at(varName, indepVarVals[varNum]);
882  }
883  }
884 
885  return;
886  }
887 
888  //......................................................................
890  fnex::MetaData const& md,
892  {
893  // Make a random number generator and use the key from the MetaData to set
894  // the seed
895  TRandom3 rand(md.Key());
896 
897  fnex::EventId eventId;
898  fnex::DataVars dataVars;
899  fnex::TruthVars truthVars;
900  fnex::WeightVars weightVars;
901  // std::vector<fnex::GENIEVar> genieVars(fnex::kSumSmallXSecJoint2017 + 1 - kMaCCQE);
902  std::vector<fnex::GENIEVar> genieVars(fnex::kCCQEPauliSupViaKF + 1 - kMaCCQE);
903 
904  eventTree->SetBranchAddress("dataVars", &dataVars);
905  eventTree->SetBranchAddress("eventId", &eventId );
906 
907  if(md.isMC){
908  eventTree->SetBranchAddress("truthVars" , &truthVars );
909  eventTree->SetBranchAddress("weightVars", &weightVars);
910 
911  // for(auto key = fnex::kMaCCQE; key < kSumSmallXSecJoint2017 + 1; ++key){
912  //The choice of which GENIE vars to keep changed for A2017 and again for A2018:
913  //check that the branch is there before setting its address
914  for(auto key = fnex::kMaCCQE; key < fnex::kCCQEPauliSupViaKF + 1; ++key){
915  if(eventTree->FindBranch(fnex::KeyToVarName(key).c_str()))
916  eventTree->SetBranchAddress(fnex::KeyToVarName(key).c_str(), &genieVars[key - fnex::kMaCCQE]);
917  else
918  LOG_WARNING("EventListManipulator")
919  << "Could not find branch "
920  << fnex::KeyToVarName(key).c_str()
921  << " in the eventlisttree. If you want to use "
922  << fnex::KeyToVarName(key).c_str()
923  << " as a systematic, you're going to have a bad time.";
924  }
925  } // end if MC
926 
927  // What is the event cap for each event type?
928  size_t treeEvents = eventTree->GetEntries();
929  float eventCap = this->EpochEventCap(md, treeEvents);
930 
931  LOG_VERBATIM("EventListManipulator")
932  << "Filling "
933  << md.ToString()
934  << " "
935  << md.Key()
936  << " with cap "
937  << eventCap
938  << " POT "
939  << ss.goodPOT
940  << " : Events "
941  << treeEvents
942  << " cap POT "
943  << ss.goodPOT * eventCap
944  << " : Events "
945  << treeEvents * eventCap;
946 
947  // Rescale POT based on how many events were actually loaded
948  // The event cap is a fraction of the total events available,
949  // so the rescaling is just the cap value.
950  ss.goodPOT *= eventCap;
951  ss.totalPOT *= eventCap;
952 
953  // Initialize the new EventList
954  fnex::EventList eventList(md, ss);
955 
956  // get the generic identifier for numu or nue selection so
957  // we can then determine the maximum energy of the neutrino
958  // events to use
960 
961  LOG_DEBUG("EventListManipulator")
962  << md.ToString()
963  << " selection is "
964  << selection
965  << " max energy: "
966  << fMaxNuEnergy.find(selection)->second;
967 
968  // check that our event list corresponds to a reasonable exposure
969  if(eventList.ListSpillSummary().goodPOT < 1.e-3)
970  throw cet::exception("EventListManipulator")
971  << "exposure for "
972  << eventList.ListMetaData().ToString()
973  << " is "
974  << eventList.ListSpillSummary().goodPOT
975  << " < 10^9 POT. That cannot be right, bail";
976 
977  for(size_t evNum = 0; evNum < treeEvents; ++evNum){
978 
979  eventTree->GetEntry(evNum);
980 
981  // If the uniformly distributed random number is greater than the
982  // cap forget the current event. This has the correct boundry
983  // conditions - if the cap is 1, then all events are used; if it is
984  // 0 then none are.
985 
986  // At some point we may need to worry about keeping more events in the
987  // tails of the spectrum, but for now it looks like we can keep at least
988  // 80% of all MetaData options, so we should not be biasing the tails
989  if(rand.Rndm() > eventCap) continue;
990 
991  fnex::EventPtr ev = std::make_shared<fnex::Event>();
992 
993  // we have to fill the event data vals first to check the
994  // the energy of this event
995  dataVars.fFake_Weight = 1.;
996  ev->dataVals.reset(new fnex::DataVarVals(dataVars));
997 
998  // don't bother filling events with greater than the maximum
999  // energy for the analysis in the trees - it just wastes time
1000  // in the fitting
1001  if(ev->dataVals->val_at(fnex::kNu_RecoE, md) > fMaxNuEnergy.find(selection)->second) continue;
1002 
1003  // now we can fill the rest of the event info because we are
1004  // keeping this event
1005  if(fUseEventId) ev->eventId.reset(new fnex::EventId(eventId));
1006  if(md.isMC){
1007  ev->mcVals.reset(new fnex::MCVarVals(truthVars,
1008  weightVars,
1009  genieVars));
1010  } // end if mc
1011 
1012  // Print out the event id information for each data event in the FD
1013  if(!md.isMC && md.detector == novadaq::cnv::kFARDET){
1014 
1015  std::string selection("numu");
1016  if(md.IsNuESelected()) selection = "nue";
1017  LOG_VERBATIM("EventListManipulator")
1018  << selection
1019  << " "
1020  << eventId;
1021 
1022  }
1023 
1024  LOG_DEBUG("EventListManipulator")
1025  << *ev;
1026 
1027  eventList.AddEvent(ev);
1028 
1029  } // for (evNum)
1030 
1031  return eventList;
1032  }
1033 
1034  //......................................................................
1036  novadaq::cnv::DetId const& det)
1037  {
1038  // Loop over the fake exposure metadatas and make sure that we
1039  // have a MCList that matches each. If not, we need to create
1040  // copies from the latest epoch lists
1041 
1042  // One could verify that all the possible subsamples are available for
1043  // each epoch too (ie interaction type, selection type, etc) and try to do
1044  // the copy of multiple different epochs, but that seems more trouble than
1045  // it is worth right now.
1046 
1047  int latestMCEpoch = 0;
1048  for(auto const& itr : lists){
1049 
1050  // only worry about lists for the MC
1051  if(!itr.first.isMC || itr.first.detector != det) continue;
1052 
1053  // if we found a later epoch, set that as the latest one and
1054  // clear out the set of MC metadata for the latest epoch
1055  if(itr.first.epoch > latestMCEpoch) latestMCEpoch = itr.first.epoch;
1056 
1057  }
1058 
1059  LOG_DEBUG("EventListManipulator")
1060  << "Latest epoch in MC is "
1061  << latestMCEpoch;
1062 
1063  // loop over the lists again, this time to fill the set of metadata for
1064  // the latest epoch. We are doing a repeat of the loop because the
1065  // fnex::EventList is an unordered map and we can't be sure that the entries
1066  // are in epoch order.
1067  std::set<fnex::MetaData> latestMCEpochMetaData;
1068  for(auto const& itr : lists){
1069 
1070  // only worry about lists for the MC
1071  if(!itr.first.isMC || itr.first.detector != det) continue;
1072 
1073  // put this metadata into the list for the latestMCEpochMetadata
1074  if(itr.first.epoch == latestMCEpoch)
1075  latestMCEpochMetaData.insert(itr.first);
1076  }
1077 
1078  // Remember fFakeExposure maps fnex::MetaData to POT exposures
1080  for(auto fakeItr : fFakeExposure){
1081 
1082 
1083  if(fakeItr.first.epoch > latestMCEpoch &&
1084  fakeItr.first.detector == det){
1085 
1086  // Need to copy the latest epoch MC lists into the current
1087  // fake data epoch
1088  for(auto const& itr : latestMCEpochMetaData){
1089 
1090  md = itr;
1091  md.epoch = fakeItr.first.epoch;
1092 
1093  LOG_DEBUG("EventListManipulator")
1094  << "Copying events for "
1095  << itr.ToString()
1096  << " into "
1097  << md.ToString();
1098 
1099  lists.emplace(md, lists.find(itr)->second);
1100  }
1101 
1102  } // end if the fake exposure epoch is later than the latest MC
1103 
1104  } // end loop over fake exposures
1105 
1106  return;
1107  }
1108 
1109 
1110  //......................................................................
1111  double EventListManipulator::RescalePOTForEventCap(double & totalEvents,
1112  size_t const& eventCap)
1113  {
1114  // recast the total variable
1115  size_t total = (size_t)totalEvents;
1116 
1117  auto rescale = this->RescalePOTForEventCap(total, eventCap);
1118 
1119  // may have had to reset the total events to the cap size
1120  totalEvents = 1. * total;
1121 
1122  return rescale;
1123  }
1124 
1125  //......................................................................
1126  double EventListManipulator::RescalePOTForEventCap(size_t & totalEvents,
1127  size_t const& eventCap)
1128  {
1129  double rescale = 1.;
1130 
1131  // Based on the cap, determine how frequently to sample the event tree to use only
1132  // that many events. Also determine when to 'stop' iterating through the tree,
1133  // because we've already
1134  if(eventCap < totalEvents){
1135  rescale = double(eventCap)/double(totalEvents);
1136  totalEvents = eventCap;
1137  }
1138 
1139  return rescale;
1140  }
1141 
1142 
1143  //......................................................................
1145  fnex::SpillSummary const& ss,
1146  fnex::EventList const& mcList)
1147  {
1148  // which detector are we using?
1149  //to det = mcList.ListMetaData().detector;
1150 
1151  // we will just copy all the events from the MC list into the fake list
1152  // we will need to just have a POT normalization weight factor to account for
1153  // the fact that some lists may have more POT than the requested amount
1154  // The requested amount was determined in FillFakeDataListsForDetector to be
1155  // the minimum used POT for all MC file types in the specified epoch and detector
1156  // of this MD.
1157  auto potNorm = ss.goodPOT / mcList.ListSpillSummary().goodPOT;
1158 
1159  // instantiate the fakeData list
1160  fnex::EventList fakeData(md, ss);
1161 
1162  LOG_DEBUG("EventListManipulator")
1163  << std::setw(50)
1164  << mcList.ListMetaData().ToString()
1165  << std::setw(15)
1166  << "Expect "
1167  << mcList.size()
1168  << " events with POT normalization of "
1169  << potNorm
1170  << " (dataPOT: "
1171  << ss.goodPOT
1172  << ", MCPOT: "
1173  << mcList.ListSpillSummary().goodPOT
1174  << ")";
1175 
1176  double wgt = 1.;
1177  double totalWgt = 0.;
1178 
1179  // make a copy of the MetaData from the original list. Use this copy to
1180  // get the fake weights for the event including any systematic changes or
1181  // oscillations. We need this copy because the ShifterAndWeighter::Weight()
1182  // function returns 1 for all MetaData with isMC = false.
1183  fnex::MetaData mdcp = mcList.ListMetaData();
1184 
1187 
1189 
1190  for(auto const& itr : mcList){
1191 
1192  fnex::EventPtr event(new fnex::Event);
1193  event->dataVals.reset(new fnex::DataVarVals(*(itr->dataVals)));
1194  event->mcVals .reset(new fnex::MCVarVals (*(itr->mcVals )));
1195 
1196  // find the weight of the event given the oscillations and systematic
1197  // uncertainties. The weight is stored in the cvVals for the event
1198  // and are applied when filling the bins in the fnex::Spectrum objects
1199  // later. Multiply the weight from the systematics and oscillations by
1200  // the POT normalizaiton to get the right total weight
1201 
1202  // if this event is in the tails of the spectrum, it has already been
1203  // weighted to account for the fact that we take all of those events
1204  // but potentially cap the events in the peak.
1205  wgt = event->dataVals->val_at(fnex::kFake_Weight, md) * potNorm;
1206  wgt *= fnex::ShifterAndWeighter::Instance()->Weight(event, mdcp, wgtScheme);
1207 
1208  // if(mdcp.selectionType == fnex::kNuESelectionPeripheral)
1209  // LOG_VERBATIM("EventListManipulator") << " peripheral wgt = " << wgt;
1210 
1211  // TODO: If this list is in the FD and is one of the background processes for the
1212  // nue appearance or is the signal in the numu disappearance, strictly speaking
1213  // we should also apply the ND data/MC ratio as a weight for the given energy
1214  // bin as well. However, we don't because we are making the ND fake data
1215  // be exactly the ND MC.
1216 
1217  totalWgt += wgt;
1218 
1219  event->dataVals->set_val_at(fnex::kFake_Weight, wgt);
1220 
1221  // now to set any shifts in the energy
1222  auto lepRecoE = event->dataVals->val_at(fnex::kLep_RecoE, md);
1223  auto hadRecoE = event->dataVals->val_at(fnex::kHad_RecoE, md);
1224 
1225  LOG_DEBUG("EventListManipulator")
1226  << "lepE "
1227  << lepRecoE
1228  << " shifted by "
1229  << lepEShift;
1230 
1231  if(md.IsNuMuSelected()){
1232  event->dataVals->set_val_at(fnex::kLep_RecoE, lepEShift * (lepRecoE - 0.1056583668) + 0.1056583668);
1233  }
1234  else if(md.IsNuESelected()){
1235  event->dataVals->set_val_at(fnex::kLep_RecoE, lepEShift * lepRecoE);
1236  }
1237 
1238  LOG_DEBUG("EventListManipulator")
1239  << " -> "
1240  << event->dataVals->val_at(fnex::kLep_RecoE, md);
1241 
1242  event->dataVals->set_val_at(fnex::kHad_RecoE, hadRecoE * hadEShift);
1243 
1244  // because this is fake data now, we no longer need the mc information,
1245  // so let's save some memory in the job
1246  event->mcVals.reset(nullptr);
1247 
1248  fakeData.AddEvent(event);
1249 
1250  } // end loop to select the fake data events
1251 
1252  LOG_VERBATIM("EventListManipulator")
1253  << "Fake Data "
1254  << std::setw(55)
1255  << md.ToString()
1256  << std::setw(15)
1257  << " represents "
1258  << fakeData.ListSpillSummary().goodPOT*1.e-8
1259  << " x 10^20 POT and "
1260  << totalWgt
1261  << " events "
1262  << totalWgt / fakeData.ListSpillSummary().goodPOT;
1263 
1264  return fakeData;
1265  }
1266 
1267  //......................................................................
1268  // make histograms for all the lists, then concatenate them in different
1269  // ways and make more histograms
1271  std::string const& dirName)
1272  {
1273  this->FillHistograms(lists, dirName);
1274 
1275  // the boolean order is: By Epoch, By Interaction Type, By File Type
1276  this->ConcatenateLists(lists, dirName, true, false, false);
1277  this->ConcatenateLists(lists, dirName, true, false, true);
1278  this->ConcatenateLists(lists, dirName, true, true, true);
1279 
1280  return;
1281  }
1282 
1283 
1284  //............................................................................
1286  fnex::EventList const& evList,
1287  std::string const& listType,
1288  double const& normalization)
1289  {
1290 
1291  // LOG_DEBUG("EventListManipulator") << "FILLING HISTOGRAMS\n";
1292 
1293  // set the name of the histogram
1294  std::string histName("neutrinoE");
1295  histName += listType;
1296 
1297  auto md = evList.ListMetaData();
1298 
1299  if (md.fileType == fnex::kCosmicBackgroundFile) histName += "Cosmic";
1300 
1301  fNeutrinoEnergy[md].push_back(tfdir.make<TH1F>(histName.c_str(),
1302  ";E_{#nu} (GeV);Events",
1303  120, 0., 30.));
1304  auto unweighted = fNeutrinoEnergy[md].back();
1305  auto weighted = fNeutrinoEnergy[md].back();
1306 
1307  // histograms to keep track of reco E --> true E mapping
1308  TH2F *recoToTrue = nullptr;
1309  TH2F *recoToTrueOscillated= nullptr; //Oscillated reco to true 2D histogram
1310  TH2F *recoToTrueCC = nullptr;
1311  TH2F *recoToTrueNC = nullptr;
1312  TH2D *recoToTrueCAF = nullptr;
1313  TH1D *recoDecompCAF = nullptr;
1314  TH1D *recoDecompCAFdata = nullptr;
1315  TH2D *recoToTrueCAFnumu = nullptr;
1316  TH1D *recoDecompCAFnumu = nullptr;
1317  TH2D *recoToTrueCAFanumu = nullptr;
1318  TH1D *recoDecompCAFanumu = nullptr;
1319  TH2D *recoToTrueCAFnue = nullptr;
1320  TH1D *recoDecompCAFnue = nullptr;
1321  TH2D *recoToTrueCAFanue = nullptr;
1322  TH1D *recoDecompCAFanue = nullptr;
1323 
1324  // histograms to evaluate BEN decomposition weights
1325  TH2F *testPion = nullptr;
1326 
1327  TH3F *benDecomp = nullptr;
1328  TH1F *benEnuKnownParents = nullptr;
1329  TH1F *benEnuUnknownParents = nullptr;
1330  TH1F *benEnuData = nullptr;
1331  TH1F *benEnuBkg = nullptr;
1332 
1333  double trueBins[101];
1334  double Emin = .5; // 500 MeV: there's really no events below there
1335  // How many edges to generate. Allow room for 0-Emin bin
1336  double N = 100-1;
1337  double A = N*Emin;
1338  trueBins[0] = 0;
1339  for(int i = 1; i <= N; ++i){
1340  trueBins[100-i] = A/i;
1341  }
1342  trueBins[100] = 120; // Replace the infinity that would be here
1343 
1344  double recoBins[21] = { 0., 0.25, 0.5, 0.75,
1345  1., 1.25, 1.5, 1.75,
1346  2., 2.25, 2.5, 2.75,
1347  3., 3.25, 3.5, 3.75,
1348  4., 4.25, 4.5, 4.75,
1349  5. };
1350 
1351 
1352  if(md.isMC){
1353  histName += "Weighted";
1354  fNeutrinoEnergy[md].push_back(tfdir.make<TH1F>(histName.c_str(),
1355  ";E_{#nu} (GeV);Events",
1356  120, 0., 30.));
1357  weighted = fNeutrinoEnergy[md].back();
1358 
1359  recoToTrue = tfdir.make<TH2F>(("recoToTrue" + listType).c_str(),
1360  ";E^{true}_{#nu};E^{reco}_{#nu}",
1361  120, 0., 30.,
1362  120, 0., 30.);
1363 
1364  recoToTrueOscillated = tfdir.make<TH2F>(("recoToTrueOscillated" + listType).c_str(),
1365  ";E^{true}_{#nu};E^{reco}_{#nu}",
1366  120, 0., 30.,
1367  120, 0., 30.);
1368 
1369  recoToTrueCC = tfdir.make<TH2F>(("recoToTrueCC" + listType).c_str(),
1370  ";E^{true}_{#nu};E^{reco}_{#nu}",
1371  120, 0., 30.,
1372  120, 0., 30.);
1373 
1374  recoToTrueNC = tfdir.make<TH2F>(("recoToTrueNC" + listType).c_str(),
1375  ";E^{true}_{#nu};E^{reco}_{#nu}",
1376  120, 0., 30.,
1377  120, 0., 30.);
1378 
1379  recoDecompCAF = tfdir.make<TH1D>(("recoDecompCAF" + listType).c_str(),
1380  ";E^{reco}_{#nu};Entries",
1381  20, 0., 5.);
1382 
1383  recoToTrueCAF = tfdir.make<TH2D>(("recoToTrueCAF" + listType).c_str(),
1384  ";E^{reco}_{#nu};E^{true}_{#nu}",
1385  20, recoBins,
1386  100, trueBins);
1387 
1388  recoDecompCAFnumu = tfdir.make<TH1D>(("recoDecompCAFnumuOrig" + listType).c_str(),
1389  ";E^{reco}_{#nu};Entries",
1390  20, 0., 5.);
1391 
1392  recoToTrueCAFnumu = tfdir.make<TH2D>(("recoToTrueCAFnumuOrig" + listType).c_str(),
1393  ";E^{reco}_{#nu};E^{true}_{#nu}",
1394  20, recoBins,
1395  100, trueBins);
1396 
1397 
1398  recoDecompCAFanumu = tfdir.make<TH1D>(("recoDecompCAFanumuOrig" + listType).c_str(),
1399  ";E^{reco}_{#nu};Entries",
1400  20, 0., 5.);
1401 
1402  recoToTrueCAFanumu = tfdir.make<TH2D>(("recoToTrueCAFanumuOrig" + listType).c_str(),
1403  ";E^{reco}_{#nu};E^{true}_{#nu}",
1404  20, recoBins,
1405  100, trueBins);
1406 
1407 
1408  recoDecompCAFnue = tfdir.make<TH1D>(("recoDecompCAFnueOrig" + listType).c_str(),
1409  ";E^{reco}_{#nu};Entries",
1410  20, 0., 5.);
1411 
1412  recoToTrueCAFnue = tfdir.make<TH2D>(("recoToTrueCAFnueOrig" + listType).c_str(),
1413  ";E^{reco}_{#nu};E^{true}_{#nu}",
1414  20, recoBins,
1415  100, trueBins);
1416 
1417  recoDecompCAFanue = tfdir.make<TH1D>(("recoDecompCAFanueOrig" + listType).c_str(),
1418  ";E^{reco}_{#nu};Entries",
1419  20, 0., 5.);
1420 
1421  recoToTrueCAFanue = tfdir.make<TH2D>(("recoToTrueCAFanueOrig" + listType).c_str(),
1422  ";E^{reco}_{#nu};E^{true}_{#nu}",
1423  20, recoBins,
1424  100, trueBins);
1425 
1426 
1427 
1428  }
1429 
1430  // BEN weight map is generated from numu selection in ND
1431  if (md.detector == novadaq::cnv::kNEARDET &&
1432  md.selectionType == fnex::kNuMuSelection ){
1433 
1434  if (md.isMC){
1435  benDecomp = tfdir.make<TH3F>(("benDecomp" + listType).c_str(),
1436  ";Ancestor target forward momentum p_{z} (GeV/c);Ancestor target transv. momentum p_{T} (GeV/c);E^{reco}_{#nu}",
1437  600, 0., 60.,
1438  100, 0., 1.,
1439  50, 0., 5.);
1440 
1441 
1442  testPion = tfdir.make<TH2F>(("testPion" + listType).c_str(),
1443  ";Ancestor target forward momentum p_{z} (GeV/c);Ancestor target transv. momentum p_{T} (GeV/c)",
1444  600, 0., 60.,
1445  100, 0., 1.);
1446 
1447  benEnuKnownParents = tfdir.make<TH1F>(("benEnuKnownParents" + listType).c_str(),
1448  ";E_{#nu} (GeV);Events",
1449  50, 0., 5.);
1450  benEnuUnknownParents = tfdir.make<TH1F>(("benEnuUnknownParents" + listType).c_str(),
1451  ";E_{#nu} (GeV);Events",
1452  50, 0., 5.);
1453 
1454  benEnuBkg = tfdir.make<TH1F>(("benEnuBkg" + listType).c_str(),
1455  ";E_{#nu} (GeV);Events",
1456  50, 0., 5.);
1457 
1458 
1459  } else {
1460 
1461  benEnuData = tfdir.make<TH1F>(("benEnuData" + listType).c_str(),
1462  ";E_{#nu} (GeV);Events",
1463  50, 0., 5.);
1464 
1465  recoDecompCAFdata = tfdir.make<TH1D>(("recoDecompCAFdata" + listType).c_str(),
1466  ";E^{reco}_{#nu};Entries",
1467  20, 0., 5.);
1468  }// end if Data instead
1469 
1470  }// end if ND and numu selection
1471 
1472  LOG_DEBUG("EventListManipulator")
1473  << " " << listType
1474  << " " << normalization
1478  << " " << fnex::KeyToVarName(fnex::kTrueE)
1480 
1481  double weight = 0.;
1482  double val = 0.;
1483  double valTrue = 0.;
1484  double trueParentPT = 0.;
1485  double trueParentPZ = 0.;
1486  int trueParentPDG = 0 ;
1487  int parentDecay = 0 ;
1488  bool numuFromPiDecay = false;
1489  bool numuFromKaDecay = false;
1490  int nuPDGorig = 0;
1491 
1492  //auto detector = md.detector;
1493 
1495 
1496  for(auto const& ev : evList){
1497 
1498  val = ev->dataVals->val_at(fnex::kNu_RecoE, md);
1499 
1500  // don't even apply the normalization weight so that we have a histogram
1501  // of exactly what was in the event list
1502  unweighted->Fill(val);
1503 
1504  // shift and weight the event if it is MC
1505  if(md.isMC){
1506 
1507  if(ev->mcVals.get() == nullptr)
1508  throw cet::exception("EventListManipulator::FillHistogram")
1509  << "MC event does not have a valid MCVarVals unique_ptr, bail";
1510 
1511  // get the weight for things like POT and MEC weighting, but don't
1512  // apply any shifteres because we are not trying to fit any systematic
1513  // uncertainties here.
1514  weight = fnex::ShifterAndWeighter::Instance()->Weight(ev, md, wgtScheme);
1515  LOG_DEBUG("EventListManipulator") << "weight = " << weight << " Normalization: " << normalization;
1516 
1517  valTrue = ev->mcVals->val_at(fnex::kTrueE);
1518 
1519  weighted->Fill(val, normalization * weight);
1520 
1521  recoToTrue->Fill(valTrue, val);
1522  recoToTrueOscillated->Fill(valTrue, val, weight);
1523  if(ev->mcVals->val_at(fnex::kTrueCCNC) == simb::kCC) recoToTrueCC->Fill(valTrue, val);
1524  else recoToTrueNC->Fill(valTrue, val);
1525 
1526  recoToTrueCAF->Fill(val, valTrue, weight );
1527  recoDecompCAF->Fill(val, weight );
1528 
1529  nuPDGorig = ev->mcVals->val_at(fnex::kTruePDGOrig);
1530 
1531  switch(nuPDGorig){
1532  case 12:
1533  recoToTrueCAFnue->Fill(val, valTrue, weight );
1534  recoDecompCAFnue->Fill(val, weight );
1535  break;
1536  case 14:
1537  recoToTrueCAFnumu->Fill(val, valTrue, weight );
1538  recoDecompCAFnumu->Fill(val, weight );
1539  break;
1540  case -12:
1541  recoToTrueCAFanue->Fill(val, valTrue, weight );
1542  recoDecompCAFanue->Fill(val, weight );
1543  break;
1544  case -14:
1545  recoToTrueCAFanumu->Fill(val, valTrue, weight );
1546  recoDecompCAFanumu->Fill(val, weight );
1547  break;
1548  }
1549 
1550  // BEND map only for the ND MC
1551  if (md.detector == novadaq::cnv::kNEARDET &&
1552  md.selectionType == fnex::kNuMuSelection ){
1553 
1554 
1555 
1556  benEnuBkg->Fill(val, normalization * weight);
1557 
1558  trueParentPDG = ev->mcVals->val_at(fnex::kTrueParentTargetPDG);
1559  parentDecay = ev->mcVals->val_at(fnex::kTrueParentDecay);
1560 
1561 
1562  numuFromPiDecay = ( (trueParentPDG==211) && (parentDecay==13) );
1563  numuFromKaDecay = ( (trueParentPDG==321) && (parentDecay==5 || parentDecay==7));
1564 
1565  // if parent was pion or K+ we use these events to find the
1566  // BEN decomposition maps
1567  // Otherwise we need to subtract them out
1568 
1569  if (numuFromPiDecay || numuFromKaDecay)
1570  benEnuKnownParents->Fill(val, normalization * weight);
1571  else
1572  benEnuUnknownParents->Fill(val, normalization * weight);
1573 
1574  if ( numuFromPiDecay ){
1575 
1576  trueParentPT = ev->mcVals->val_at(fnex::kTrueParentPT);
1577  trueParentPZ = ev->mcVals->val_at(fnex::kTrueParentPZ);
1578 
1579  benDecomp->Fill(trueParentPZ, trueParentPT, val, normalization * weight);
1580  testPion->Fill(trueParentPZ, trueParentPT, normalization * weight);
1581  }
1582  }
1583 
1584  } else if(md.detector == novadaq::cnv::kNEARDET &&
1585  md.selectionType == fnex::kNuMuSelection) {
1586 
1587  benEnuData->Fill(val, normalization);
1588  recoDecompCAFdata->Fill(val, 1 );
1589 
1590  } // end if/else mc
1591 
1592 
1593  } // end loop over events in the list
1594 
1595  //LOG_DEBUG("EventListManipulator") << "DONE FILLING HISTOGRAMS\n";
1596 
1597  return;
1598  }
1599 
1600  //......................................................................
1602  std::string const& dirName)
1603  {
1604 
1605  //LOG_DEBUG("EventListManipulator") << "FILLING HISTOGRAMS II\n";
1606 
1607  // make a neutrino energy histogram for each of the event lists
1608 
1609  // instantiate the check histograms
1610 
1611  // create a directory in which to store these histograms
1612  auto histDir = fTFS->mkdir(dirName);
1613 
1614  // find the data exposures for the different detectors, if there is no
1615  // data list for a given detector use the fFakeExposure value. If there is
1616  // no fake exposure, set it to be 5e20 POT, ie 5e8 * 1e12 (the unit of POT)
1617 
1618  // we need to know what epochs there are in this list, so make a set of those
1619  std::set<std::string> epochs;
1620  for(auto const& itr : eventLists){
1621  epochs.insert(itr.first.EpochString());
1622  LOG_DEBUG("EventListManipulator")
1623  << "epoch "
1624  << itr.first.EpochString()
1625  << " is present";
1626  }
1627 
1628  std::vector<fnex::MetaData> mds;
1629  for(auto const& epoch : epochs){
1630  LOG_DEBUG("EventListManipulator")
1631  << "making MetaData for data in epoch "
1632  << epoch;
1633 
1653 
1654  }
1655 
1656  ExposureMap mdToPOT;
1657 
1658  for(auto const& md : mds){
1659  if(eventLists.count(md) > 0)
1660  mdToPOT[md] = eventLists.find(md)->second.ListSpillSummary().goodPOT;
1661 
1662  if(fFakeExposure.count(md) > 0)
1663  mdToPOT[md] = fFakeExposure.find(md)->second;
1664 
1665  LOG_DEBUG("EventListManipulator")
1666  << std::setw(20)
1667  << md.ToString()
1668  << " "
1669  << mdToPOT[md]
1670  << " x 10^12 POT";
1671  }
1672 
1673  double normalization = 1.;
1674  double dataPOT = 0.;
1675 
1676  fnex::MetaData dataMD;
1678  std::string listName;
1679  // loop over the event lists and fill the histograms
1680  for(auto const& it : eventLists){
1681 
1682  md = it.first;
1683 
1684  auto const& evList = it.second;
1685 
1686  // already checked that the goodPOT value in each list is at least
1687  // 0.001, ie 10^9 POT
1688  dataMD = md;
1689  dataMD.isMC = false;
1690  dataMD.fileType = fnex::kDataFile;
1692  if( mdToPOT.count(dataMD) > 0 ) dataPOT = mdToPOT[dataMD];
1693  else dataPOT = 0.;
1694 
1695  normalization = dataPOT / evList.ListSpillSummary().goodPOT;
1696 
1697  LOG_DEBUG("EventListManipulator")
1698  << "metadata: "
1699  << std::setw(20) << md.ToString()
1700  << " corresponds to " << evList.size()
1701  << " events in " << evList.ListSpillSummary().goodPOT
1702  << " x 10^12 POT and normalization " << normalization
1703  << " data POT: " << dataPOT;
1704 
1705  // avoid duplicating histograms between data and fake data
1706  listName = md.ToString();
1707  if(dirName.find("Fake") != std::string::npos)
1708  listName.replace(listName.find("Data"), 4, "FakeData");
1709 
1710  FillHistogram(histDir, evList, listName, normalization);
1711 
1712  } // end loop over event lists
1713 
1714  //LOG_DEBUG("EventListManipulator") << "DONE FILLING HISTOGRAMS II\n";
1715 
1716  return;
1717  }
1718 
1719  //......................................................................
1721  std::string const& dirName,
1722  bool const& byEpoch,
1723  bool const& byInteraction,
1724  bool const& byFile)
1725  {
1726 
1727  // make event lists that concatenate the different epochs and interaction types
1729  fnex::EventListMap concatLists;
1731  fnex::FileType_t fType;
1732  std::string epoch;
1733  std::unordered_map<fnex::MetaData, int, fnex::MetaDataHasher> prevEpochMap;
1734  for(auto const& itr : lists){
1735  md = itr.first;
1736 
1737  epoch = byEpoch ? "0z" : md.EpochString();
1738  iType = byInteraction ? fnex::kUnknownInteraction : md.interactionType;
1739  fType = byFile ? fnex::kUnknownFileType : md.fileType;
1740 
1741  fnex::MetaData search(md.isMC,
1742  md.detector,
1743  fType,
1744  md.selectionType,
1745  iType,
1746  epoch);
1747 
1748  if(concatLists.count(search) < 1){
1749  concatLists.emplace(search, fnex::EventList(search));
1750  prevEpochMap[search] = md.epoch;
1751  concatLists.at(search).AddGoodPOT(itr.second.ListSpillSummary().goodPOT);
1752  concatLists.at(search).AddGoodSpills(itr.second.ListSpillSummary().numGoodSpills);
1753  }
1754 
1755  if(md.epoch != prevEpochMap[search]){
1756  prevEpochMap[search] = md.epoch;
1757  concatLists.at(search).AddGoodPOT(itr.second.ListSpillSummary().goodPOT);
1758  concatLists.at(search).AddGoodSpills(itr.second.ListSpillSummary().numGoodSpills);
1759  }
1760 
1761  for(auto const& ev : itr.second) concatLists.at(search).AddEvent(ev);
1762  } // end loop to make combined histograms
1763 
1764  std::string directoryName = dirName + "Concatenated";
1765 
1766  if(byEpoch ) directoryName += "ByEpoch";
1767  if(byInteraction) directoryName += "ByInteractionType";
1768  if(byFile ) directoryName += "ByFileType";
1769 
1770  this->FillHistograms(concatLists, directoryName);
1771 
1772  // release the concatenated lists as we don't need them any more
1773  for(auto itr : concatLists) itr.second.resize(0);
1774  concatLists.clear();
1775 
1776  return;
1777  }
1778 
1779  //......................................................................
1781  {
1782 
1783  LOG_DEBUG("EventListManipulator")
1784  << "FILLING FAKE DATA LISTS - there are "
1785  << fFakeExposure.size()
1786  << " fake exposures to fill";
1787 
1788  // keep track of the minimum exposure for each MC combination of
1789  // file type, epoch and detector to use in the weighting of the
1790  // fake data lists
1791  std::map<int, double> detEpochToMinPOT;
1792  int edfKey = 0;
1793  float curPOT = 0.;
1794  float minPOT = 1.e6;
1795  for(auto const& itr : lists){
1796 
1797  if(itr.first.isMC){
1798  edfKey = itr.first.DetectorKey() + itr.first.EpochKey();
1799  curPOT = itr.second.ListSpillSummary().goodPOT;
1800  if(detEpochToMinPOT.count(edfKey) > 0){
1801  minPOT = detEpochToMinPOT[edfKey];
1802  detEpochToMinPOT[edfKey] = std::min(minPOT, curPOT);
1803  }
1804  else
1805  detEpochToMinPOT[edfKey] = curPOT;
1806  }
1807 
1808  }
1809 
1810  // for(auto itr : detEpochToMinPOT)
1811  // LOG_VERBATIM("EventListManipulator")
1812  // << "Minimum POT "
1813  // << itr.first
1814  // << " "
1815  // << itr.second
1816  // << " x 10^20 POT";
1817 
1818  // check if we are to fit fake data instead of the input data
1819  // Having a fake exposure > 0 means we should fake the data
1820  if(fFakeExposure.size() < 1) return;
1821 
1822  // Check to ensure that we have the MC lists available for any requested FakeData epochs.
1825 
1826  // set up the shifters and weighters for the appropriate detectors
1827  auto shftWgt = fnex::ShifterAndWeighter::Instance();
1828  shftWgt->InitShiftsAndWeightsToUse(fFakeDataPoint, fCalibHistFile);
1829 
1830  // create a directory in which to store these histograms
1831  auto histDir = fTFS->mkdir("FakeData");
1832 
1833  //-----------------------------------
1834  // First fill ND fake data
1835  // Need to do this first to properly
1836  // apply the proportional decomp to
1837  // the FD fake data
1838  this->FillFakeDataListsForDetector(lists,
1840  detEpochToMinPOT);
1841 
1842  //-----------------------------------
1843  // Now add the NuE true energy ratio weighter
1844  if(fUseExtrapolation) shftWgt->CalcNearTrueEnergyRatioWeights(&lists);
1845 
1846  //-----------------------------------
1847  // Now fill FD fake data
1848  this->FillFakeDataListsForDetector(lists,
1850  detEpochToMinPOT);
1851 
1852  // At this point, we have faked the data in the event lists
1853  // to have high stats distributions scaled to the appropriate
1854  // POT exposures.
1855  // This would be the time to use those lists to create Poisson distributed
1856  // fake data
1857  // Use the fPoissonFakeData parameter to turn it on and off in fhicl
1858 
1859  // loop over all the fake data metadata (should just be by epoch)
1860  if(fPoissonFakeData){
1861  for(auto fakeItr : fFakeExposure){
1862  auto const& md = fakeItr.first;
1863 
1864  // get the fake list for this md from the collection of event lists
1865  // and then use it to make a Poisson distributed list
1866  auto const& fakeList = lists.find(md)->second;
1867 
1868  // we made the poisson list, so erase the high stats one and replace it
1869  lists.erase(md);
1870  auto retVal = lists.emplace(md, this->CreatePoissonFakeData(fakeList.ListMetaData(),
1871  fakeList.ListSpillSummary(),
1872  fakeList));
1873 
1874  auto poissonDir = fTFS->mkdir("FakePoisson" + md.ToString());
1875  FillHistogram(poissonDir, retVal.first->second, md.ToString(), 1.);
1876 
1877  } // end loop over fake exposures
1878 
1879  } // end if making Poisson fake data
1880 
1881 
1882 // LOG_VERBATIM("EventListManipulator")
1883 // << "After filling fake data ";
1884 //
1885 // for(auto itr : lists){
1886 //
1887 // if(itr.first.isMC) continue;
1888 //
1889 // LOG_VERBATIM("EventListManipulator")
1890 // << std::setw(55)
1891 // << itr.first.ToString()
1892 // << std::setw(15)
1893 // << " represents "
1894 // << itr.second.ListSpillSummary().goodPOT*1.e-8
1895 // << " x 10^20 POT and "
1896 // << itr.second.size()
1897 // << " events "
1898 // << itr.second.size() / itr.second.ListSpillSummary().goodPOT;
1899 // }
1900 
1901 
1902  return;
1903  }
1904 
1905  //----------------------------------------------------------------------------
1907  novadaq::cnv::DetId const& det,
1908  std::map<int, double> const& detEpochToMinPOT)
1909  {
1910  // we are going to want to keep track of the total
1911  // weight of all events in each detector broken down by
1912  // interaction type and file type
1913  std::map<long, double> keyToFakeWeight;
1914 
1915  long key = 0;
1916 
1917  double fakePOT = 0.;
1918 
1919  fnex::EventListMap fakeLists;
1920 
1921  // Remember fFakeExposure maps fnex::MetaData to POT exposures
1922  for(auto fakeItr : fFakeExposure){
1923 
1924  auto const& md = fakeItr.first;
1925 
1926  if(md.detector != det) continue;
1927 
1928  LOG_DEBUG("EventListManipulator")
1929  << "attempt to make fake data for "
1930  << md.ToString();
1931 
1932  if( !this->UseEventsFromMetaData(md, fnex::kBoth) ){
1933  LOG_DEBUG("EventListManipulator")
1934  << " Fake expsure for "
1935  << md.ToString()
1936  << " failed check to use events from that metadata";
1937  continue;
1938  }
1939 
1940  key = md.DetectorKey() + md.EpochKey();
1941  fakePOT = std::min(fakeItr.second, detEpochToMinPOT.find(key)->second);
1942 
1943  LOG_DEBUG("EventListManipulator")
1944  << "Fake POT for "
1945  << key
1946  << " "
1947  << md.detector
1948  << " "
1949  << fakePOT
1950  << " requested "
1951  << fakeItr.second
1952  << " min "
1953  << detEpochToMinPOT.find(key)->second;
1954 
1955  fnex::SpillSummary ss(fakePOT,
1956  fakePOT,
1957  0.,
1958  1,
1959  1);
1960 
1961 
1962  fnex::EventList fakeData(md);
1963 
1964  for(auto listItr : lists){
1965 
1966  LOG_DEBUG("EventListManipulator")
1967  << listItr.first.ToString();
1968 
1969  // don't test on the interaction type here because data has them all
1970  if(listItr.first.detector == md.detector &&
1971  listItr.first.selectionType == md.selectionType &&
1972  listItr.first.epoch == md.epoch &&
1973  listItr.first.isMC){
1974 
1975  key = (md.DetectorKey() +
1976  listItr.first.FileKey() +
1977  listItr.first.InteractionKey() +
1978  md.EpochKey());
1979 
1980  auto partialList = this->AddFakeDataToEventList(md, ss, listItr.second);
1981 
1982  for(auto itr : partialList){
1983  fakeData.AddEvent(itr);
1984 
1985  // use the fake weight if it exists, otherwise we made Poisson
1986  // distributed fake data and we just add 1 to the total weight
1987  keyToFakeWeight[key] += itr->dataVals->val_at(fnex::kFake_Weight, md);
1988  }
1989 
1990  } // end if the right detector, selection, and epoch
1991  } // done looping over MC lists
1992 
1993  fakeData.AddGoodPOT(ss.goodPOT);
1994  fakeData.AddGoodSpills(1);
1995  fakeLists.emplace(md, fakeData);
1996 
1997  } // end loop over fake data lists
1998 
1999 // for(auto const& itr : keyToFakeWeight)
2000 // LOG_DEBUG("EventListManipulator")
2001 // << itr.first
2002 // << " "
2003 // << itr.second;
2004 
2005  // clear out the map as we no longer need it
2006  keyToFakeWeight.clear();
2007 
2008  // make a histograms of the fake data for future reference
2009  //this->FillHistograms(fakeLists, "FakeData");
2010 
2011  // the boolean order is: By Epoch, By Interaction Type, By File Type
2012  //this->ConcatenateLists(fakeLists, "FakeData", fSumSubEpochs, false, false);
2013 
2014  // now replace the data lists in eventLists with our faked up list
2015  for(auto itr : fakeLists){
2016  lists.erase(itr.first);
2017  lists.emplace(itr.first, itr.second);
2018 
2019  itr.second.resize(0);
2020  LOG_DEBUG("EventListManipulator")
2021  << " REPLACING DATA " << itr.first.ToString();
2022 
2023  }
2024 
2025  // clear out the fakeLists
2026  fakeLists.clear();
2027 
2028  return;
2029  }
2030 
2031  //----------------------------------------------------------------------------
2033  float & scale,
2034  fnex::SpillSummary const& ss)
2035  {
2036  std::unique_ptr<TH1F> hist;
2038  auto epoch = md.epoch;
2039 
2040  // get the file containing the histograms
2042 
2043  // adjust the file name for FHC or RHC
2044  // Only adjust the names for the numu analysis.
2045  // Nue analysis has a single root file
2046  if(md.IsNuMuSelected() == true){
2047  if(md.BeamType() == fnex::kFHC &&
2048  fileName.find("rhc") != std::string::npos)
2049  fileName.replace(fileName.find("rhc"), 3, "fhc");
2050  else if(md.BeamType() == fnex::kRHC &&
2051  fileName.find("fhc") != std::string::npos)
2052  fileName.replace(fileName.find("fhc"), 3, "rhc");
2053  }//end of numu analysis
2054 
2055  std::string filePath;
2056  cet::search_path sp("FW_SEARCH_PATH");
2057  if( !sp.find_file(fileName, filePath) )
2058  throw cet::exception("EventListMaker")
2059  << "Cannot find cosmic background file "
2060  << fileName;
2061 
2062  TFile tf(filePath.c_str(), "READ");
2063 
2064  // set the POT accounting for each epoch. The wrinkle is that each histogram
2065  // is only for a period, not an epoch, which is a subset of the period so we
2066  // want to adjust for that. Use the same relative ratios between the epochs
2067  // in a given period. Take the inverse of those ratios to make the POT
2068  // normalization come out correctly
2069  // the value 1000, 2000, etc is the period * 1000, the 97, 98, 99, etc is a, b, c
2070  LOG_DEBUG("EventListMaker")
2071  << "Inside EventListManipulator::FarDetCosmicBackgroundHistAndScale--"
2072  << " Metadata is : "
2073  << md.ToString();
2075  if (epoch == 1097 ){
2076  scale = 1.;
2077  histName = "hp1_all";
2078  hist = std::make_unique<TH1F>(*(dynamic_cast<TH1F*>(tf.Get(histName.c_str()))));
2079  }
2080  else if(epoch == 2097 ||
2081  epoch == 2098 ||
2082  epoch == 2099 ){
2083  scale = 3.;
2084  histName = "hp2_all";
2085  hist = std::make_unique<TH1F>(*(dynamic_cast<TH1F*>(tf.Get(histName.c_str()))));
2086  }
2087  else if(epoch == 3098 ){
2088  scale = 16.66;
2089  histName = "hp3_all";
2090  hist = std::make_unique<TH1F>(*(dynamic_cast<TH1F*>(tf.Get(histName.c_str()))));
2091  }
2092  else if(epoch == 3099 ){
2093  scale = 1.06;
2094  histName = "hp3_all";
2095  hist = std::make_unique<TH1F>(*(dynamic_cast<TH1F*>(tf.Get(histName.c_str()))));
2096  }
2097  else if(epoch == 1000 ){
2098  scale = 1.;
2099  histName = "hp1_all";
2100  hist = std::make_unique<TH1F>(*(dynamic_cast<TH1F*>(tf.Get(histName.c_str()))));
2101  }
2102  else if(epoch == 2000 ){
2103  scale = 1.;
2104  histName = "hp2_all";
2105  hist = std::make_unique<TH1F>(*(dynamic_cast<TH1F*>(tf.Get(histName.c_str()))));
2106  }
2107  else if(epoch == 3000 ){
2108  scale = 1.0;
2109  histName = "hp3_all";
2110  hist = std::make_unique<TH1F>(*(dynamic_cast<TH1F*>(tf.Get(histName.c_str()))));
2111  }
2112  else if(epoch == 0){
2113  // we are combining all the epochs, so we have to do that for all the
2114  // backgrounds.
2115  scale = 1.;
2116  hist = std::make_unique<TH1F>(*(dynamic_cast<TH1F*>(tf.Get("hp1 _all"))));
2117  hist->Add(dynamic_cast<TH1F*>(tf.Get("hp2_all")));
2118  hist->Add(dynamic_cast<TH1F*>(tf.Get("hp3_all")));
2119  }
2120  else
2121  throw cet::exception("EventListManipulator")
2122  << "cannot determine period for making FD cosmic ray background list";
2123 
2124  if(!hist)
2125  LOG_WARNING("EventListManipulator")
2126  << "Unable to find necessary histograms for numu cosmic backgrounds";
2127  } //////////------------------------------------------------ Numu Quantile 1 -------------------
2128  else if(md.selectionType == fnex::kNuMuSelectionQ1){
2129  // SJB NOTE: 2018 cosmic histograms only seem to have a single histogram for
2130  // each numu quantile (i.e. not split up into periods)
2131  // This may be liable to change or I may be missing something so I'll leave this
2132  // code here for now.
2133 
2134  /*if (epoch >= 1000 && epoch <2000 ){
2135  scale = 1.;
2136  histName = "p1q1";
2137  hist = std::make_unique<TH1F>(*(dynamic_cast<TH1F*>(tf.Get(histName.c_str()))));
2138  }
2139  else if(epoch >= 2000 && epoch < 3000){
2140  scale = 1.;
2141  histName = "p2q1";
2142  hist = std::make_unique<TH1F>(*(dynamic_cast<TH1F*>(tf.Get(histName.c_str()))));
2143  }
2144  else if(epoch >= 3000 && epoch < 4000 ){
2145  scale = 1;
2146  histName = "p3q1";
2147  hist = std::make_unique<TH1F>(*(dynamic_cast<TH1F*>(tf.Get(histName.c_str()))));
2148  }
2149  else if(epoch >= 5000 && epoch < 6000 ){
2150  scale = 1.0;
2151  histName = "p5q1";
2152  hist = std::make_unique<TH1F>(*(dynamic_cast<TH1F*>(tf.Get(histName.c_str()))));
2153  }
2154  else if(epoch == 0){*/
2155  if(epoch == 0 || (epoch > 999 && epoch < 8000 ) ) {
2156  // we are combining all the epochs, so we have to do that for all the
2157  // backgrounds.
2158  if(md.BeamType() == fnex::kFHC) scale = (fTotalFDDataPOTFHC / ss.goodPOT);
2159  else if(md.BeamType() == fnex::kRHC) scale = (fTotalFDDataPOTRHC / ss.goodPOT);
2160  histName = "cosmics_q1";
2161  hist = std::make_unique<TH1F>(*(dynamic_cast<TH1F*>((TH1F*)tf.Get(histName.c_str()))));
2162  // hist = std::make_unique<TH1F>(*(dynamic_cast<TH1F*>(tf.Get("hp1 _all"))));
2163  // hist->Add(dynamic_cast<TH1F*>(tf.Get("hp2_all")));
2164  // hist->Add(dynamic_cast<TH1F*>(tf.Get("hp3_all")));
2165  }
2166  else
2167  throw cet::exception("EventListManipulator")
2168  << "cannot determine period for making FD cosmic ray background list";
2169 
2170  if(!hist)
2171  LOG_WARNING("EventListManipulator")
2172  << "Unable to find necessary histograms for numu cosmic backgrounds";
2173  } //////////------------------------------------------------ Numu Quantile 2 -------------------
2174  else if(md.selectionType == fnex::kNuMuSelectionQ2){
2175  /*if (epoch >= 1000 && epoch <2000 ){
2176  scale = 1.;
2177  histName = "p1q2";
2178  hist = std::make_unique<TH1F>(*(dynamic_cast<TH1F*>(tf.Get(histName.c_str()))));
2179  }
2180  else if(epoch >= 2000 && epoch < 3000){
2181  scale = 1.;
2182  histName = "p2q2";
2183  hist = std::make_unique<TH1F>(*(dynamic_cast<TH1F*>(tf.Get(histName.c_str()))));
2184  }
2185  else if(epoch >= 3000 && epoch < 4000 ){
2186  scale = 1;
2187  histName = "p3q2";
2188  hist = std::make_unique<TH1F>(*(dynamic_cast<TH1F*>(tf.Get(histName.c_str()))));
2189  }
2190  else if(epoch >= 5000 && epoch < 6000 ){
2191  scale = 1.0;
2192  histName = "p5q2";
2193  hist = std::make_unique<TH1F>(*(dynamic_cast<TH1F*>(tf.Get(histName.c_str()))));
2194  }
2195  else if(epoch == 0){*/
2196  if(epoch == 0 || ( epoch > 999 && epoch < 8000 ) ) {
2197  // we are combining all the epochs, so we have to do that for all the
2198  // backgrounds.
2199  if(md.BeamType() == fnex::kFHC) scale = (fTotalFDDataPOTFHC / ss.goodPOT);
2200  else if(md.BeamType() == fnex::kRHC) scale = (fTotalFDDataPOTRHC / ss.goodPOT);
2201  histName = "cosmics_q2";
2202  hist = std::make_unique<TH1F>(*(dynamic_cast<TH1F*>((TH1F*)tf.Get(histName.c_str()))));
2203  // hist = std::make_unique<TH1F>(*(dynamic_cast<TH1F*>(tf.Get("hp1 _all"))));
2204  // hist->Add(dynamic_cast<TH1F*>(tf.Get("hp2_all")));
2205  // hist->Add(dynamic_cast<TH1F*>(tf.Get("hp3_all")));
2206  }
2207  else
2208  throw cet::exception("EventListManipulator")
2209  << "cannot determine period for making FD cosmic ray background list";
2210 
2211  if(!hist)
2212  LOG_WARNING("EventListManipulator")
2213  << "Unable to find necessary histograms for numu cosmic backgrounds";
2214  } //////////------------------------------------------------ Numu Quantile 3 -------------------
2215  else if(md.selectionType == fnex::kNuMuSelectionQ3){
2216  /*if (epoch >= 1000 && epoch <2000 ){
2217  scale = 1.;
2218  histName = "p1q3";
2219  hist = std::make_unique<TH1F>(*(dynamic_cast<TH1F*>(tf.Get(histName.c_str()))));
2220  }
2221  else if(epoch >= 2000 && epoch < 3000){
2222  scale = 1.;
2223  histName = "p2q3";
2224  hist = std::make_unique<TH1F>(*(dynamic_cast<TH1F*>(tf.Get(histName.c_str()))));
2225  }
2226  else if(epoch >= 3000 && epoch < 4000 ){
2227  scale = 1;
2228  histName = "p3q3";
2229  hist = std::make_unique<TH1F>(*(dynamic_cast<TH1F*>(tf.Get(histName.c_str()))));
2230  }
2231  else if(epoch >= 5000 && epoch < 6000 ){
2232  scale = 1.0;
2233  histName = "p5q3";
2234  hist = std::make_unique<TH1F>(*(dynamic_cast<TH1F*>(tf.Get(histName.c_str()))));
2235  }
2236  else if(epoch == 0){*/
2237  if(epoch == 0 || ( epoch > 999 && epoch < 8000 ) ) {
2238  // we are combining all the epochs, so we have to do that for all the
2239  // backgrounds.
2240  if(md.BeamType() == fnex::kFHC) scale = (fTotalFDDataPOTFHC / ss.goodPOT);
2241  else if(md.BeamType() == fnex::kRHC) scale = (fTotalFDDataPOTRHC / ss.goodPOT);
2242  histName = "cosmics_q3";
2243  hist = std::make_unique<TH1F>(*(dynamic_cast<TH1F*>((TH1F*)tf.Get(histName.c_str()))));
2244  // hist = std::make_unique<TH1F>(*(dynamic_cast<TH1F*>(tf.Get("hp1 _all"))));
2245  // hist->Add(dynamic_cast<TH1F*>(tf.Get("hp2_all")));
2246  // hist->Add(dynamic_cast<TH1F*>(tf.Get("hp3_all")));
2247  }
2248  else
2249  throw cet::exception("EventListManipulator")
2250  << "cannot determine period for making FD cosmic ray background list";
2251 
2252  if(!hist)
2253  LOG_WARNING("EventListManipulator")
2254  << "Unable to find necessary histograms for numu cosmic backgrounds";
2255  } //////////------------------------------------------------ Numu Quantile 4 -------------------
2256  else if(md.selectionType == fnex::kNuMuSelectionQ4){
2257  /*if (epoch >= 1000 && epoch <2000 ){
2258  scale = 1.;
2259  histName = "p1q4";
2260  hist = std::make_unique<TH1F>(*(dynamic_cast<TH1F*>(tf.Get(histName.c_str()))));
2261  }
2262  else if(epoch >= 2000 && epoch < 3000){
2263  scale = 1.;
2264  histName = "p2q4";
2265  hist = std::make_unique<TH1F>(*(dynamic_cast<TH1F*>(tf.Get(histName.c_str()))));
2266  }
2267  else if(epoch >= 3000 && epoch < 4000 ){
2268  scale = 1;
2269  histName = "p3q4";
2270  hist = std::make_unique<TH1F>(*(dynamic_cast<TH1F*>(tf.Get(histName.c_str()))));
2271  }
2272  else if(epoch >= 5000 && epoch < 6000 ){
2273  scale = 1.0;
2274  histName = "p5q4";
2275  hist = std::make_unique<TH1F>(*(dynamic_cast<TH1F*>(tf.Get(histName.c_str()))));
2276  }
2277  else if(epoch == 0){*/
2278  if(epoch == 0 || ( epoch > 999 && epoch < 8000 ) ) {
2279  // we are combining all the epochs, so we have to do that for all the
2280  // backgrounds.
2281 
2282  if(md.BeamType() == fnex::kFHC) scale = (fTotalFDDataPOTFHC / ss.goodPOT);
2283  else if(md.BeamType() == fnex::kRHC) scale = (fTotalFDDataPOTRHC / ss.goodPOT);
2284  histName = "cosmics_q4";
2285  hist = std::make_unique<TH1F>(*(dynamic_cast<TH1F*>((TH1F*)tf.Get(histName.c_str()))));
2286 
2287  // hist = std::make_unique<TH1F>(*(dynamic_cast<TH1F*>(tf.Get("hp1 _all"))));
2288  // hist->Add(dynamic_cast<TH1F*>(tf.Get("hp2_all")));
2289  // hist->Add(dynamic_cast<TH1F*>(tf.Get("hp3_all")));
2290  }
2291  else
2292  throw cet::exception("EventListManipulator")
2293  << "cannot determine period for making FD cosmic ray background list";
2294 
2295  if(!hist)
2296  LOG_WARNING("EventListManipulator")
2297  << "Unable to find necessary histograms for numu cosmic backgrounds";
2298  }
2299  //Get the cosmics for the nue selections
2300  else if(md.selectionType > fnex::kNuMuSelection &&
2302 
2303  // scale factor is the inverse of the fraction of each period's POT to total FD Data POT
2304  // different scale factors for the FHC and RHC needed.
2305  if(md.BeamType() == fnex::kFHC) scale = (fTotalFDDataPOTFHC / ss.goodPOT);
2306  else if(md.BeamType() == fnex::kRHC) scale = (fTotalFDDataPOTRHC / ss.goodPOT);
2307 
2308 // if(md.BeamType() == fnex::kFHC){
2309 // LOG_VERBATIM("EventListManipulator")
2310 // << "Inside FarDetCosmicBackgroundHistAndScale---"
2311 // << " Metadata : "
2312 // << md.ToString()
2313 // << " fTotalFDDataPOTFHC : "
2314 // << fTotalFDDataPOTFHC
2315 // << " ss.goodPOT : "
2316 // << ss.goodPOT
2317 // << " Scale : "
2318 // << scale;
2319 // }
2320 // if(md.BeamType() == fnex::kRHC){
2321 // LOG_VERBATIM("EventListManipulator")
2322 // << "Inside FarDetCosmicBackgroundHistAndScale---"
2323 // << " Metadata : "
2324 // << md.ToString()
2325 // << " fTotalFDDataPOTRHC : "
2326 // << fTotalFDDataPOTRHC
2327 // << " ss.goodPOT : "
2328 // << ss.goodPOT
2329 // << " Scale : "
2330 // << scale;
2331 // }
2332 
2333  if(fAnalysisYear=="2018") hist = this->Cosmics2018(md, tf);
2334  else if(fAnalysisYear=="2017") hist = this->Cosmics2017(md);
2335  else{
2336 
2337  int firstbin = 0;
2338 
2339  //TODO: Need to adjust scale factor for each epoch seperately and correctly
2340  if (md.selectionType == fnex::kNuESelectionLowPID ) firstbin = 0;
2341  else if(md.selectionType == fnex::kNuESelectionMidPID ) firstbin = 10;
2342  else if(md.selectionType == fnex::kNuESelectionHighPID) firstbin = 20;
2343 
2344  // for the nue analysis we have to get a couple of histograms and do some
2345  // manipulation before we can make our histogram
2346  auto cosmicLiveTime = dynamic_cast<TH1D*>(tf.Get("spec_cvn_cvn2d_cosmic/livetime"));
2347  auto dataLiveTime = dynamic_cast<TH1D*>(tf.Get("spec_cvn_cvn2d_numi/livetime"));
2348  auto cosmicCounts_org = dynamic_cast<TH1D*>(tf.Get("spec_cvn_cvn2d_cosmic/hist"));
2349 
2350  // check that the histograms all exist, throw an exception otherwise
2351  if(!cosmicLiveTime ||
2352  !dataLiveTime ||
2353  !cosmicCounts_org)
2354  LOG_WARNING("EventListManipulator")
2355  << "Unable to find necessary histograms for nue cosmic backgrounds";
2356  else{
2357 
2358  //Need to remake individual CAF histograms for Low, Mid and High CVN bins because
2359  //in the CAF histogram file all three histograms are present as one histogram
2360  hist = std::make_unique<TH1F>("cosmicCounts", "", 10, 0, 5);
2361  for(int i = firstbin + 1; i < firstbin + 11; ++i){
2362  hist->SetBinContent(i - firstbin, cosmicCounts_org->GetBinContent(i));
2363  }
2364  hist->Scale(dataLiveTime->Integral() / cosmicLiveTime->Integral());
2365  }
2366  }//end of 2016 analysis cosmics
2367  }// end of nue selections
2368 
2369  return hist;
2370  }
2371 
2372  //----------------------------------------------------------------------------
2373  // The cosmic ray background for the Numu and Nue analyses are calculated using histograms
2374  // by CAFAna. Use those histograms here to create a list of fnex::Events reflecting
2375  // the background. Bump up the exposure so that we have even numbers of events
2376  // rather than fractional weights for them.
2378  fnex::SpillSummary const& ss)
2379  {
2380  // scale factor to provide smooth distributions
2381  float scale = 5.e3;
2382 
2383  // this scale factor is determined by the FarDetCosmicBackgroundHistAndScale
2384  // method and is used to account for the different POT weighting of each
2385  // epoch and period.
2386  float extraScale = 1.;
2387 
2388  // get the histogram
2389  auto hist = this->FarDetCosmicBackgroundHistAndScale(md, extraScale, ss);
2390 
2391  // All MC POT is in units of 1e12 when put into the event lists
2392  fnex::SpillSummary outSS(scale * extraScale * ss.totalPOT,
2393  scale * extraScale * ss.goodPOT,
2394  scale * ss.totalNumSpills,
2395  scale * ss.numGoodSpills,
2396  scale * ss.liveTime);
2397 
2398  // create the metadata for this list
2399  fnex::MetaData outMD(true,
2402  md.selectionType,
2404  md.EpochString());
2405 
2406  // make the EventList to return
2407  fnex::EventList eventList(outMD, outSS);
2408 
2409  if(!hist){
2410  LOG_WARNING("EventListManipulator")
2411  << "No cosmic ray background histogram found for selection "
2412  << md.selectionType
2413  << " , metadat: "
2414  << md.ToString()
2415  << " return an empty event list";
2416  return eventList;
2417  }
2418 
2419  LOG_DEBUG("EventListManipulator")
2420  << std::setw(50)
2421  << "Cosmic Background exposure: "
2422  << md.ToString()
2423  << std::setw(20)
2424  << outSS.totalPOT
2425  << std::setw(20)
2426  << ss.totalPOT;
2427 
2428  // make the total number of events be around 1e4
2429  double totalEvents = (int)(hist->Integral() * scale);
2430  double energy = 0.;
2431  for(int e = 0; e < totalEvents; ++e){
2432  energy = hist->GetRandom();
2433 
2434  fnex::EventPtr event(new fnex::Event);
2435  event->dataVals = std::make_unique<fnex::DataVarVals>();
2436  event->mcVals = std::make_unique<fnex::MCVarVals >();
2437 
2438  event->dataVals->set_val_at(fnex::kLep_RecoE, energy );
2439  event->dataVals->set_val_at(fnex::kHad_RecoE, 0. );
2440 
2441  event->mcVals->set_val_at(fnex::kTrueE, 0. );
2442  event->mcVals->set_val_at(fnex::kTruePDG, 13 );
2443  event->mcVals->set_val_at(fnex::kTruePDGOrig, 13 );
2444  event->mcVals->set_val_at(fnex::kTrueCCNC, 0. );
2445  event->mcVals->set_val_at(fnex::kTrueIntType, fnex::kCosmicMuon);
2446  event->mcVals->set_val_at(fnex::kXSecCVPPFX_Weight, 1. );
2447  event->mcVals->set_val_at(fnex::kRPACCQE_Weight, 1. );
2448 
2449  eventList.AddEvent(event);
2450  } // end loop to make event list
2451 
2452  return eventList;
2453 
2454  }
2455 
2456 
2457  //----------------------------------------------------------------------------
2458  //Calculate cosmic background from eventlisttrees
2460  fnex::SpillSummary const& ss,
2461  fnex::EventListMap const& listmap)
2462  {
2463  fnex::MetaData outMD(true,
2466  md.selectionType,
2468  md.EpochString());
2469  fnex::EventList return_cosmicList(outMD);
2470 
2471  for(auto cosmic_itr : listmap){
2472  auto const& cosmicMD = cosmic_itr.first;
2473 
2474  if(md.selectionType == cosmicMD.selectionType &&
2475  md.epoch == cosmicMD.epoch){
2476 
2477  //need to declare this here, so it has the correct spill summary
2478  fnex::EventList temp_cosmicList(outMD, cosmic_itr.second.ListSpillSummary());
2479 
2480  if(cosmicMD.fileType != fnex::kCosmicBackgroundFile)
2481  throw cet::exception("EventListManipulator")
2482  << "Uh oh, cosmicMD filetype is NOT cosmic, something's wrong!";
2483  //if(cosmicMD.OscType != fnex::kCosmicEvent) LOG_DEBUG("EventListManipulator") << "CosmicMD osctype is NOT cosmic";
2484 
2485  LOG_DEBUG("EventListManipulator")
2486  << "Match at "
2487  << md.ToString()
2488  << " and "
2489  << cosmicMD.ToString()
2490  << "\nscale of "
2491  << ss.goodPOT
2492  << " / "
2493  << cosmic_itr.second.ListSpillSummary().goodPOT
2494  << " = "
2495  << ss.goodPOT / cosmic_itr.second.ListSpillSummary().goodPOT;
2496 
2497  float POTScale = ss.goodPOT / cosmic_itr.second.ListSpillSummary().goodPOT;
2498 
2499  // loop over events in cosmic event list, weighting by POTScale
2500  for(auto evt : cosmic_itr.second) {
2501  evt->dataVals->set_val_at(fnex::kFake_Weight, POTScale);
2502  evt->mcVals = std::make_unique<fnex::MCVarVals >();
2503 
2504  evt->mcVals->set_val_at(fnex::kTrueE, 0. );
2505  evt->mcVals->set_val_at(fnex::kTruePDG, 13 );
2506  evt->mcVals->set_val_at(fnex::kTruePDGOrig, 13 );
2507  evt->mcVals->set_val_at(fnex::kTrueCCNC, 0. );
2508  evt->mcVals->set_val_at(fnex::kTrueIntType, fnex::kCosmicMuon);
2509  evt->mcVals->set_val_at(fnex::kXSecCVPPFX_Weight, 1. );
2510  evt->mcVals->set_val_at(fnex::kRPACCQE_Weight, 1. );
2511 
2512  temp_cosmicList.AddEvent(evt);
2513  }//finished adding events to the new cosmicList
2514 
2515  //temp_cosmicList is about to die. Make sure we can return it!
2516  return_cosmicList = temp_cosmicList;
2517 
2518  }//finished if(md matches)
2519  }//finished looping over all entries in initial cosmicList
2520 
2521  return return_cosmicList;
2522 
2523  }
2524 
2525  //......................................................................
2526  // The input list is assumed to have been created using the exact MC
2527  // weighted as appropriate for the POT and any oscillations
2528  // This method fills histograms from that list and then draws the correct
2529  // number of events out of it based on the weighted total
2531  fnex::SpillSummary const& ss,
2532  fnex::EventList const& origList)
2533  {
2534  // make the EventList to return
2535  fnex::EventList outList(md, ss);
2536 
2537  auto hadFixedB = fnex::VarBinning::Instance()->BinInfo("Had_RecoE");
2538  auto lepFixedB = fnex::VarBinning::Instance()->BinInfo("Lep_RecoE");
2539 
2540  // Histogram of hadronic vs. leptonic energies from which to draw energies
2541  // for Poisson events from, using 50 MeV divisions for bins.
2542  TH2F hist("poissonFakeDataDraw", "poissonFakeDataDraw",
2543  2 * hadFixedB->nBins(), &(hadFixedB->BinLowEdges()[0]),
2544  2 * lepFixedB->nBins(), &(lepFixedB->BinLowEdges()[0]));
2545 
2546  double lepE = 0.;
2547  double hadE = 0.;
2548  double totWgt = 0.;
2549 
2550  for(auto const& ev : origList){
2551  hist.Fill(ev->dataVals->val_at(fnex::kHad_RecoE, md),
2552  ev->dataVals->val_at(fnex::kLep_RecoE, md),
2553  ev->dataVals->val_at(fnex::kFake_Weight, md));
2554 
2555  totWgt += ev->dataVals->val_at(fnex::kFake_Weight, md);
2556  }
2557 
2558  int totalEvents = fRandom.Poisson(totWgt);
2559 
2560  LOG_DEBUG("EventListManipulator")
2561  << "total initial events: "
2562  << hist.Integral()
2563  << " total weight: "
2564  << totWgt
2565  << " poisson events: "
2566  << totalEvents;
2567 
2568  for(int e = 0; e < totalEvents; ++e){
2569  hist.GetRandom2(hadE, lepE);
2570 
2571  fnex::EventPtr event(new fnex::Event);
2572  event->dataVals = std::make_unique<fnex::DataVarVals>();
2573 
2574  event->dataVals->set_val_at(fnex::kLep_RecoE, lepE);
2575  event->dataVals->set_val_at(fnex::kHad_RecoE, hadE);
2576  event->dataVals->set_val_at(fnex::kFake_Weight, 1. );
2577 
2578  outList.AddEvent(event);
2579  }
2580 
2581  LOG_DEBUG("EventListManipulator")
2582  << "There are "
2583  << outList.size()
2584  << " Poisson distributed events in the "
2585  << outList.ListMetaData().ToString()
2586  << " list";
2587 
2588  return outList;
2589  }
2590 
2591  //......................................................................
2592  // Complete hackery here - I know the ratio of NuMu CC events / POT x 1e-12
2593  // is 0.003, so find a correction factor based on that number and apply it
2594  // per epoch and interaction type
2595  // TODO: remove this hack
2597  {
2598  std::map<int, double> corrScale;
2599  fnex::MetaData md2;
2600  fnex::SpillSummary ss2;
2601 
2602  //Medbh 11/11/17 hack
2603  //int is PERIOD
2604  std::unordered_map<int, double> ndmcEventsPerPOT;
2605  ndmcEventsPerPOT[1000] = 3.12122e-15;
2606  ndmcEventsPerPOT[2000] = 3.09095e-15;
2607  ndmcEventsPerPOT[3000] = 3.0248e-15;
2608  ndmcEventsPerPOT[5000] = 2.9574e-15;
2609  std::unordered_map<int, int> ndmcEventsPerPeriod;
2610  ndmcEventsPerPeriod[1000] = 0;
2611  ndmcEventsPerPeriod[2000] = 0;
2612  ndmcEventsPerPeriod[3000] = 0;
2613  ndmcEventsPerPeriod[5000] = 0;
2614 
2615  std::unordered_map<int, double> ndmcPOT;
2616  ndmcPOT[1000] = 0;
2617  ndmcPOT[2000] = 0;
2618  ndmcPOT[3000] = 0;
2619  ndmcPOT[5000] = 0;
2620 
2621  for(auto itr : lists){
2622  md2 = itr.first;
2623  if(md2.isMC &&
2624  md2.detector == novadaq::cnv::kNEARDET ){//&&
2625  // md2.interactionType == fnex::kNuMuCC ){//&&
2626  // md2.IsNuMuSelected()){
2627 
2628  for(auto & fileName : fFNEXEventLabels){
2629  TFile *f = TFile::Open(fileName.c_str());
2630  for(auto & dirName : fTreeDirectories){
2631 
2632  // find the event tree using the metadata information
2633  std::string eventTreeName(dirName);
2634  eventTreeName += "/";
2635  auto mdString = md2.ToString();
2636 
2637  TTree *eventTree = dynamic_cast<TTree*>(f->Get((eventTreeName + mdString).c_str()));
2638 
2639  LOG_VERBATIM("EventListManipulator")
2640  << "Adding "
2641  << eventTree->GetEntries()
2642  << " from "
2643  << md2.ToString()
2644  << " to "
2645  << ndmcEventsPerPeriod[md2.Period()];
2646  ndmcEventsPerPeriod[md2.Period()] += eventTree->GetEntries();
2647  }
2648  }
2649  }
2650  }
2651 
2652  for(auto itr : lists){
2653  md2 = itr.first;
2654  if(md2.isMC &&
2655  md2.detector == novadaq::cnv::kNEARDET){ //&&
2656  // md2.interactionType == fnex::kNuMuCC ){//&&
2657  // md2.IsNuMuSelected()){
2658 
2659  // double dataToPOT = (1. * itr.second.size() ) / itr.second.ListSpillSummary().goodPOT;
2660 
2661  // corrScale[md2.epoch] = .003 / dataToPOT;
2662 
2663  // LOG_DEBUG("EventListManipulator")
2664  // << md2.ToString()
2665  // << " dataToPOT: "
2666  // << dataToPOT
2667  // << " correction: "
2668  // << corrScale[md2.epoch];
2669 
2670  //Medbh 11/11/17 hack
2671  LOG_DEBUG("EventListManipulator")
2672  << "Working ND magic: "
2673  << md2.ToString()
2674  << " POT is now "
2675  << ndmcEventsPerPeriod[md2.Period()]
2676  << " / "
2677  << ndmcEventsPerPOT[md2.Period()]
2678  << " = "
2679  << (double)ndmcEventsPerPeriod[md2.Period()] / ndmcEventsPerPOT[md2.Period()];
2680 
2681  ndmcPOT[md2.Period()] = (double)ndmcEventsPerPeriod[md2.Period()] / ndmcEventsPerPOT[md2.Period()];
2682  ndmcPOT[md2.Period()] *= 1.e-12;
2683  }
2684  } // end loop to determine the correction factor for each epoch
2685 
2686  fnex::EventListMap eventLists2;
2687  for(auto itr : lists){
2688 
2689  md2 = itr.first;
2690  ss2 = itr.second.ListSpillSummary();
2691 
2692  if(md2.isMC &&
2693  md2.detector == novadaq::cnv::kNEARDET ){//&&
2694  // md2.IsNuMuSelected() ){
2695  // ss2.totalPOT /= corrScale[md2.epoch];
2696  // ss2.goodPOT /= corrScale[md2.epoch];
2697  ss2.totalPOT = ndmcPOT[md2.Period()];
2698  ss2.goodPOT = ndmcPOT[md2.Period()];
2699  }
2700 
2701  fnex::EventList newlist(md2, ss2);
2702  for(size_t evCtr = 0; evCtr < itr.second.size(); ++evCtr)
2703  newlist.AddEvent(itr.second.EventNum(evCtr));
2704 
2705  eventLists2.emplace(md2, newlist);
2706  }
2707 
2708  // swap out the lists and clear out eventLists2
2709  std::swap(lists, eventLists2);
2710 
2711  for(auto itr : eventLists2) itr.second.resize(0);
2712  eventLists2.clear();
2713 
2714  return;
2715  }
2716 
2717  //----------------------------------------------------------------------------
2718  // method to write out a file with the capped event lists to speed up processing
2720  {
2721 
2722  TFile *file = new TFile("cappedEventListsTrees.root", "RECREATE");
2723 
2724  auto dir = file->mkdir("list");
2725  dir->cd();
2726 
2730  std::string epoch;
2731 
2732  // first the metadata tree
2733  // leave this pointer dangling because the fileService created & manages it
2734  TTree *metadataTree = new TTree( "metadata", "Event list metadata" );
2735  metadataTree->Branch("metadata", &md);
2736  metadataTree->Branch("spillsummary", &ss);
2737 
2738  for(auto const& itr : listMap){
2739 
2740  md = itr.first;
2741  ss = itr.second.ListSpillSummary();
2742 
2743  metadataTree->Fill();
2744 
2745  // same principle as above for the event tree
2746  TTree *eventTree = new TTree( md.ToString().c_str(), "FNEX event records");
2747 
2748  fnex::SerializeEvents(eventTree,
2749  md,
2750  itr.second);
2751 
2752 
2753  eventTree->Write();
2754  delete eventTree;
2755 
2756  } // end loop over the lists
2757 
2758  metadataTree->Write();
2759  file->Write();
2760  delete metadataTree;
2761  delete file;
2762 
2763  return;
2764  }
2765  //----------------------------------------------------------------------------
2766  //2017 nue cosmics
2768 
2769  LOG_DEBUG("EventListManipulator::Cosmics2017")
2770  << "Inside EventListManipulator::Cosmics2017 function---";
2771  LOG_DEBUG("EventListManipulator::Cosmics2017")
2772  << "Metadata: "
2773  << md.ToString();
2774 
2775  std::unique_ptr<TH1F> hist;
2776 
2777  // get the file containing the histograms
2779 
2780  std::string filePath;
2781  cet::search_path sp("FW_SEARCH_PATH");
2782  if( !sp.find_file(fileName, filePath) )
2783  throw cet::exception("EventListMaker")
2784  << "Cannot find cosmic background file "
2785  << fileName;
2786 
2787  TFile tf(filePath.c_str(), "READ");
2788 
2789  //Get the file for FD data live time
2790  std::string FDLiveTimefileName = ( md.IsNuMuSelected() ) ? "" : fNuEDataHistFile;
2791 
2792  std::string FDLiveTimefilePath;
2793  cet::search_path FDLiveTimesp("FW_SEARCH_PATH");
2794  if( !FDLiveTimesp.find_file(FDLiveTimefileName, FDLiveTimefilePath) )
2795  throw cet::exception("EventListMaker")
2796  << "Cannot find data live time file "
2797  << FDLiveTimefileName;
2798 
2799  TFile FDLiveTimetf(FDLiveTimefilePath.c_str(), "READ");
2800 
2801  int firstbin = 0;
2802 
2803  //In the CAF nue 2017 cosmic file, there are only 9 bins for each analysis bin
2804  if (md.selectionType == fnex::kNuESelectionLowPID ) firstbin = 0;
2805  else if(md.selectionType == fnex::kNuESelectionMidPID ) firstbin = 9;
2806  else if(md.selectionType == fnex::kNuESelectionHighPID) firstbin = 18;
2807  else if(md.selectionType == fnex::kNuESelectionPeripheral) firstbin = 27;
2808 
2809  // for the nue analysis we have to get a couple of histograms and do some
2810  // manipulation before we can make our histogram
2811  auto cosmicLiveTime = dynamic_cast<TH1D*>(tf .Get("All/cosm_merged_4bin/livetime"));
2812  auto cosmicCounts_org = dynamic_cast<TH1D*>(tf .Get("All/cosm_merged_4bin/hist"));
2813  auto dataLiveTime = dynamic_cast<TH1D*>(FDLiveTimetf.Get("All/spect_merged_4bin/livetime"));
2814 
2815  // check that the histograms all exist, throw an exception otherwise
2816  if(!cosmicLiveTime ||
2817  !cosmicCounts_org ||
2818  !dataLiveTime)
2819  LOG_WARNING("EventListManipulator")
2820  << "Unable to find necessary histograms for nue cosmic backgrounds";
2821  else{
2822 
2823  //Need to remake individual CAF histograms for Low, Mid, High CVN and Peripheral bins because
2824  //in the CAF histogram file all four histograms are present as one histogram
2825  hist = std::make_unique<TH1F>("cosmicCounts", "", 10, 0, 5);
2826  for(int i = firstbin + 1; i < firstbin + 10; ++i){
2827  hist->SetBinContent(i - firstbin, cosmicCounts_org->GetBinContent(i));
2828  LOG_DEBUG("Cosmics2017")
2829  << " Bin: "
2830  << i
2831  << " Input cosmic bin content: "
2832  << cosmicCounts_org->GetBinContent(i)
2833  << " New bin: "
2834  << i - firstbin
2835  << " Output cosmic bin content: "
2836  << hist->GetBinContent(i - firstbin);
2837  }
2838  hist->Scale(dataLiveTime->Integral() / cosmicLiveTime->Integral());
2839  LOG_DEBUG("Cosmics2017")
2840  << " Data live time: "
2841  << dataLiveTime->Integral()
2842  << " Cosmic live time: "
2843  << cosmicLiveTime->Integral();
2844  }
2845 
2846  return hist;
2847  }//end of Cosmics2017 function
2848 //----------------------------------------------------------------------------
2849 //Make ND Peripheral eventists using HIGH CVN eventlists
2851 
2852  // first check that we want an ND peripheral list
2855  return;
2856  }
2857 
2858  LOG_DEBUG("EventListManipulator::CreateNDPeripheralLists")
2859  << "Inside EventListManipulator::CreateNDPeripheralLists function";
2860 
2861  LOG_DEBUG("EventListManipulator::CreateNDPeripheralLists")
2862  << "Size of Input EventListMap : "
2863  << eventlistmap.size();
2864 
2866 
2867  //loop over input eventlist map and make ND peripheral eventlists
2868  //using nue selected ND HIGH CVN event list
2869  for(const auto it : eventlistmap){
2870  md = it.first;
2871  auto eventlist = it.second;
2872 
2873  //only keep ND High CVN metadata
2874  if(md.detector == novadaq::cnv::kNEARDET &&
2876 
2877  //New metadata
2878  fnex::MetaData new_md(md.isMC,
2879  md.detector,
2880  md.fileType,
2882  md.interactionType,
2883  md.EpochString());
2884 
2885  //Fill new eventlist using new metadata as a key and same POT as input eventlist
2886  fnex::EventList new_EventList(new_md, eventlist.ListSpillSummary());
2887 
2888  //Use events from input eventlist and fill events in ND peripheral lists
2889  for(auto evPtr : eventlist) new_EventList.AddEvent(evPtr);
2890 
2891  LOG_DEBUG("EventListManipulator::CreateNDPeripheralLists")
2892  << " Input metadata: "
2893  << md.ToString()
2894  << " Input eventlist size: "
2895  << eventlist.size()
2896  << " New metadata: "
2897  << new_md.ToString()
2898  << " New eventlist size: "
2899  << new_EventList.size();
2900 
2901  //Add nue selected ND Peripheral metadata to the eventlist map
2902  eventlistmap.emplace(new_md, new_EventList);
2903  }//end of condition
2904  }//end of loop over eventlist map
2905 
2906  LOG_DEBUG("EventListManipulator::CreateNDPeripheralLists")
2907  << "Size of EventList Map: "
2908  << eventlistmap.size();
2909 
2910  return;
2911 
2912 }//end of CreateNDPeripheralLists function
2913 //----------------------------------------------------------------------------
2914 //2018 nue cosmics
2916  TFile & tf){
2917 
2918  LOG_DEBUG("EventListManipulator::Cosmics2018")
2919  << "Inside EventListManipulator::Cosmics2018 function---"
2920  << "\nMetadata: "
2921  << md.ToString()
2922  << " Cosmic file name : "
2923  << *(tf.GetName());
2924 
2925  std::unique_ptr<TH1F> hist;
2926 
2927  // The first and the last bins for each PID
2928  int firstbin = 0;
2929  int lastbin = 0;
2930 
2932  firstbin = 0;
2933  lastbin = 10;
2934  }
2936  firstbin = 9;
2937  lastbin = 19;
2938  }
2940  firstbin = 18;
2941  lastbin = 28;
2942  }
2943 
2944  // for the nue analysis we have to get a couple of histograms and do some
2945  // manipulation before we can make our histogram
2946  std::unique_ptr<TH1D> cosmicLiveTime ;
2947  std::unique_ptr<TH1D> cosmicCounts_org ;
2948 
2949  if(md.BeamType() == fnex::kFHC){
2950  cosmicLiveTime = std::make_unique<TH1D>(*(dynamic_cast<TH1D*>(tf.Get("cosmic_spect_fhc/livetime"))));
2951  cosmicCounts_org = std::make_unique<TH1D>(*(dynamic_cast<TH1D*>(tf.Get("cosmic_spect_fhc/hist"))));
2952  }
2953  if(md.BeamType() == fnex::kRHC){
2954  cosmicLiveTime = std::make_unique<TH1D>(*(dynamic_cast<TH1D*>(tf.Get("cosmic_spect_rhc/livetime"))));
2955  cosmicCounts_org = std::make_unique<TH1D>(*(dynamic_cast<TH1D*>(tf.Get("cosmic_spect_rhc/hist"))));
2956  }
2957 
2958 
2959  // check that the histograms all exist, throw an exception otherwise
2960  if(!cosmicLiveTime ||
2961  !cosmicCounts_org)
2962  LOG_WARNING("EventListManipulator")
2963  << "Unable to find necessary histograms for nue cosmic backgrounds";
2964  else{
2965 
2966  //Need to remake individual CAF histograms for Low, High CVN and Peripheral bins because
2967  //in the CAF histogram file all four histograms are present as one histogram but in different bins
2968  hist = std::make_unique<TH1F>("cosmicCounts", "", 10, 0, 5);
2969  for(int i = firstbin + 1; i < lastbin; ++i){
2970  hist->SetBinContent(i - firstbin, cosmicCounts_org->GetBinContent(i));
2971  LOG_DEBUG("Cosmics2018")
2972  << " Bin: "
2973  << i
2974  << " Input cosmic bin content: "
2975  << cosmicCounts_org->GetBinContent(i)
2976  << " New bin: "
2977  << i - firstbin
2978  << " Output cosmic bin content: "
2979  << hist->GetBinContent(i - firstbin);
2980  }
2981 
2982  //NuMu Live Time from CAFAna/Analysis/Exposures.h
2983  double kAna2018Livetime = 0;
2984  if(md.BeamType() == fnex::kFHC) kAna2018Livetime = 438.1584;
2985  else if(md.BeamType() == fnex::kRHC) kAna2018Livetime = 179.6173;
2986 
2987  hist->Scale(kAna2018Livetime / cosmicLiveTime->Integral());
2988  LOG_VERBATIM("Cosmics2018")
2989  << " Data live time: "
2990  << kAna2018Livetime
2991  << " Cosmic live time: "
2992  << cosmicLiveTime->Integral();
2993  }
2994 
2995  return hist;
2996 }//end of Cosmics2018 function
2997 //----------------------------------------------------------------------------
2998 
2999 } // end namespace
fnex::EventListMap Deserialize(fnex::DataMC_t dataMC=fnex::kBoth)
#define LOG_DEBUG(stream)
Definition: Messenger.h:149
float fMaxEventsPerTree
maximum number of events to use per tree
static const unsigned char kNu_RecoE
Definition: VarVals.h:36
void CheckDataMCListConsistency(std::set< fnex::MetaDataLite > const &dataSets, std::set< fnex::MetaDataLite > const &mcSets, fnex::DataMC_t const &dataMC, fnex::EventListMap &eventLists)
std::string fCalibHistFile
file containing calibration systematic histograms
enum fnex::weight_type WeightType_t
fileName
Definition: plotROC.py:78
bool IsNuESelected() const
Definition: Structs.cxx:265
void FixNearDetMCExposures(fnex::EventListMap &lists)
std::string fNuECosmicHistFile
file containing the cosmic background histograms
bool IsNuMuSelected() const
Definition: Structs.cxx:256
set< int >::iterator it
std::string fNuMuCosmicHistFile
file containing the cosmic background histograms
enum fnex::file_type FileType_t
void CreateNDPeripheralLists(fnex::EventListMap &listmap)
fnex::MetaData const & ListMetaData() const
Definition: Event.h:131
EnergyHistMap fNeutrinoEnergy
neutrino energy
static const unsigned char kTrueParentPT
Definition: VarVals.h:48
void AddGoodSpills(unsigned int nSpills)
Definition: Event.h:128
const Var weight
void AddGoodPOT(double pot)
Definition: Event.h:126
double fTotalNDDataPOTFHC
total ND data POT FHC
Base container for the MC related Vars that constitute an event.
Definition: VarVals.h:930
void InitShiftsAndWeightsToUse(fnex::InputPoint const &curPoint, std::string CalibFilename="FNEX/data/calibration/CalibrationSystematics2018.root")
const char * p
Definition: xmltok.h:285
void SerializeEvents(TTree *eventTree, fnex::MetaData const &md, std::vector< fnex::EventInfo > const &events)
Definition: Event.cxx:82
TFileDirectory mkdir(std::string const &dir, std::string const &descr="")
std::string fAnalysisYear
analysis year, "2016" for SA and "2017" for TA
Float_t ss
Definition: plot.C:24
long const Key() const
Definition: Structs.cxx:234
bool fSumSubEpochs
do we want to sum all sub epochs into a single one
enum fnex::interaction_type InteractionType_t
static const unsigned char kXSecCVPPFX_Weight
Definition: VarVals.h:56
std::string const & eventTreeName()
::xsd::cxx::tree::exception< char > exception
Definition: Database.h:225
std::map< fnex::MetaData, fnex::EventList > EventListMap
Definition: Event.h:186
long const Period() const
Definition: Structs.h:62
TString ip
Definition: loadincs.C:5
static const unsigned char kFake_Weight
Definition: VarVals.h:39
std::string const EpochString() const
Definition: Structs.cxx:101
Create a list of fnex::Events to be used in fits.
std::unique_ptr< TH1F > Cosmics2017(fnex::MetaData const &md)
fnex::EventListMap DeserializeData()
Timing fit.
TRandom3 fRandom
Random number generator to use creating fake data lists.
std::map< fnex::SelectionType_t, double > fMaxNuEnergy
maximum neutrino energy to go into the lists
static ShifterAndWeighter * Instance()
std::string find_file(std::string const &filename) const
std::size_t size() const
Definition: Event.h:110
fnex::EventList CreatePoissonFakeData(fnex::MetaData const &md, fnex::SpillSummary const &ss, fnex::EventList const &origList)
unsigned long int numGoodSpills
Definition: Structs.h:108
void ConcatenateLists(fnex::EventListMap const &lists, std::string const &dirName, bool const &byEpoch, bool const &byInteraction, bool const &byFile)
double fTotalFDDataPOTRHC
total FD data POT RHC
void swap(art::HLTGlobalStatus &lhs, art::HLTGlobalStatus &rhs)
static const unsigned char kTruePDGOrig
Definition: VarVals.h:46
void AddEvent(EventPtr ev)
Definition: Event.h:115
static const unsigned char kTrueCCNC
Definition: VarVals.h:44
const XML_Char * s
Definition: expat.h:262
std::unique_ptr< TH1F > FarDetCosmicBackgroundHistAndScale(fnex::MetaData const &md, float &scale, fnex::SpillSummary const &ss)
Double_t scale
Definition: plot.C:25
std::string const ToString() const
Definition: Structs.cxx:114
Far Detector at Ash River, MN.
enum fnex::sel_type SelectionType_t
art::ServiceHandle< art::TFileService > fTFS
TFileService.
fnex::EventList CreateFarDetCosmicBackgroundsFromHist(fnex::MetaData const &md, fnex::SpillSummary const &ss)
double fTotalFDDataPOTFHC
total FD data POT FHC
std::map< long, float > fEventCaps
bool fPoissonFakeData
make the fake data Poisson distributed
list outList
Definition: runNovaSAM.py:344
fnex::EventList AddFakeDataToEventList(fnex::MetaData const &md, fnex::SpillSummary const &ss, fnex::EventList const &mcList)
enum fnex::dataMC_type DataMC_t
void MakeHistograms(fnex::EventListMap const &lists, std::string const &dirName)
double Weight(fnex::EventPtr const &curEvent, fnex::MetaData const &md, fnex::WeightType_t const wgtType=kAllWgt)
fnex::InputPoint fFakeDataPoint
the fake data input point
void CheckMCListsForFakeDataEpochs(fnex::EventListMap &lists, novadaq::cnv::DetId const &det)
T get(std::string const &key) const
Definition: ParameterSet.h:231
static const unsigned char kMaCCQE
Definition: VarVals.h:70
double energy
Definition: plottest35.C:25
#define pot
fnex::InputPoint fInitialGuess
the fake data input point
float fFake_Weight
Weight for fake data events.
Definition: VarVals.h:735
Near Detector in the NuMI cavern.
float EpochEventCap(fnex::MetaData const &md, int treeEvents)
static const unsigned char kCCQEPauliSupViaKF
Definition: VarVals.h:81
fnex::EventList CreateFarDetCosmicBackgroundsFromEventList(fnex::MetaData const &md, fnex::SpillSummary const &ss, fnex::EventListMap const &listmap)
std::unique_ptr< TH1F > Cosmics2018(fnex::MetaData const &md, TFile &tf)
Get the nue 2018 CAFAna cosmic histogram.
static const unsigned char kTrueIntType
Definition: VarVals.h:45
novadaq::cnv::DetId detector
Definition: Structs.h:50
double ShiftAmount(unsigned char const &param)
void ExtractEventsFromFile(TFile *eventFile, std::string const &dirName, fnex::DataMC_t const &dataMC, std::pair< fnex::MetaData, fnex::SpillSummary > const &mdToSS, fnex::EventListMap &eventLists)
static const unsigned char kTrueParentDecay
Definition: VarVals.h:50
static VarBinning * Instance()
Definition: VarBinning.cxx:15
const ana::Var wgt
fnex::FileType_t fileType
Definition: Structs.h:51
std::string const DetectorString() const
Definition: Structs.h:93
std::set< int > fPeriodsToUse
which periods to use in the analysis
bool fUseEventId
use if you need access to run/subrun/event/etc
#define LOG_WARNING(category)
std::string fNuEDataHistFile
file containing the data live time for cosmic normalization
ExposureMap fFakeExposure
POT in 1e12 of an exposure to fake up.
static const unsigned char kTrueParentTargetPDG
Definition: VarVals.h:51
fnex::SelectionType_t selectionType
Definition: Structs.h:52
void DetermineSelections(std::vector< fhicl::ParameterSet > const &selections)
double fTotalNDDataPOTRHC
total ND data POT RHC
unsigned long int totalNumSpills
Definition: Structs.h:107
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
void FillHistogram(art::TFileDirectory &tfdir, fnex::EventList const &evList, std::string const &listType, double const &normalization)
static const double A
Definition: Units.h:82
T * make(ARGS...args) const
void reconfigure(const fhicl::ParameterSet &p)
std::unique_ptr< Correlations > fCorrelations
const std::string cInteractionType_Strings[9]
Definition: Constants.h:150
bool UseEventsFromMetaData(fnex::MetaData const &md, fnex::DataMC_t const dataMC)
double RescalePOTForEventCap(size_t &totalEvents, size_t const &eventCap)
std::unordered_map< fnex::MetaData, double, fnex::MetaDataHasher > ExposureMap
fnex::InteractionType_t interactionType
Definition: Structs.h:53
bool IsMCVar(unsigned char const &key)
Definition: VarVals.h:604
std::map< std::string, fnex::SpillSummary > fFDPOTSum
data POT per epoch in the FD
bool IsDataVar(unsigned char const &key)
Definition: VarVals.h:617
TDirectory * dir
Definition: macro.C:5
static const unsigned char kRPACCQE_Weight
Definition: VarVals.h:57
void FillFakeDataLists(fnex::EventListMap &lists)
void CalcNearTrueEnergyRatioWeights(const fnex::EventListMap *pEventListMap)
std::shared_ptr< Event > EventPtr
Definition: Event.h:50
std::vector< std::string > fFNEXEventLabels
Labels in input files holding FNEX Events.
static const unsigned char kTrueE
Definition: VarVals.h:42
void WriteCappedListsAsTrees(fnex::EventListMap const &list)
bool fLoadAllEventLists
force all event lists to be loaded
void FillHistograms(fnex::EventListMap const &evLists, std::string const &dirName)
std::vector< std::string > fTreeDirectories
directory holding the input trees
fnex::EventList FillEventList(TTree *eventTree, fnex::MetaData const &md, fnex::SpillSummary &ss)
bool fUseExtrapolation
are we extrapolating or using covariance matrices
TFile * file
Definition: cellShifts.C:17
EventListManipulator(fhicl::ParameterSet const &pset, fnex::InputPoint const &initialGuess)
void ExtractMetaDataFromFile(TFile *metadataFile, std::string const &dirName, fnex::DataMC_t const &dataMC, std::map< fnex::MetaData, fnex::SpillSummary > &mdToSS)
static const unsigned char kHad_RecoE
Definition: VarVals.h:34
std::map< novadaq::cnv::DetId, std::set< fnex::SelectionType_t > > fSelectionsToUse
which selection types to use?
void LoadVarVals(fnex::Event &fnexEvent, bool isMC, std::vector< std::string > const &indepVarNames, std::vector< float > const &indepVarVals)
std::string KeyToVarName(unsigned char const &key)
Definition: VarVals.h:445
std::unique_ptr< MCVarVals > mcVals
Definition: Event.h:36
#define LOG_VERBATIM(category)
fnex::BeamType_t const BeamType() const
Definition: Structs.cxx:165
const std::string cSelectionType_Strings[11]
Definition: Constants.h:101
T min(const caf::Proxy< T > &a, T b)
Float_t e
Definition: plot.C:35
long StringToKey(std::string const &str)
Definition: Structs.cxx:15
static const unsigned char kLep_RecoE
Definition: VarVals.h:35
void FillFakeDataListsForDetector(fnex::EventListMap &lists, novadaq::cnv::DetId const &det, std::map< int, double > const &detEpochToMinPOT)
fnex::Binning const * BinInfo(std::string const &varName)
Definition: VarBinning.cxx:65
fnex::SpillSummary const & ListSpillSummary() const
Definition: Event.h:123
static const unsigned char kTrueParentPZ
Definition: VarVals.h:49
Base container for the MC related Vars that constitute an event.
Definition: VarVals.h:874
void CreateFarDetCosmicBackgrounds(fnex::EventListMap &eventLists)
std::unique_ptr< DataVarVals > dataVals
Definition: Event.h:35
static const unsigned char kTruePDG
Definition: VarVals.h:43