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 #include <memory>
9 #include <iomanip>
10 #include <fstream>
11 
12 #include "TFile.h"
13 
15 
24 
26 
27 namespace cmf{
28 
29  //----------------------------------------------------------------------------
31  {
32  fEventCaps.clear();
33 
34  this->reconfigure(p);
35  }
36 
37  //......................................................................
39  {
40  fExposure.clear();
41  }
42 
43  //......................................................................
45  {
46 
47  fRandom.SetSeed(pset.get<unsigned int>("RandomSeed", 0));
48 
49  auto detectors = pset.get<std::vector<std::string>>("DetectorsToUse", std::vector<std::string>({"NearDet", "FarDet"}));
50  for(auto const& itr : detectors) fDetectors.insert(cmf::StringToDetectorType(itr));
51 
52  fCMFEventLabels = pset.get< std::vector<std::string> >("CMFEventLabels" , std::vector<std::string>() );
53  fTreeDirectories = pset.get< std::vector<std::string> >("TreeDirectories" , std::vector<std::string>() );
54  fUseEventId = pset.get< bool >("UseEventId" , false );
55  fMaxNuEnergy[cmf::kNuESelection ] = pset.get< double >("MaxNuEEnergyForAnalysis" , 5.0 );
56  fMaxNuEnergy[cmf::kNuMuSelection] = pset.get< double >("MaxNuMuEnergyForAnalysis" , 10.0 );
57  fMaxNuEnergy[cmf::kNCSelection ] = pset.get< double >("MaxNCEnergyForAnalysis" , 20.0 );
58  fMaxEventsPerTree = pset.get< float >("MaxEventsPerTree" , 1.e6 );
59  fMinEventsPerTree = pset.get< float >("MinEventsPerTree" , 2.e4 );
60  fFillTextFile = pset.get< bool >("FillTextFile" , false );
61  fDeserializeCosmics = pset.get< bool >("DeserializeCosmics" , true );
62 
63  // if we're printing the event lists we want the event ID information!
64  if (fFillTextFile){
65  fUseEventId = true;
66  fMaxEventsPerTree = 1.e6;
67  }
68 
69  // for(auto & dir : fTreeDirectories) dir.append("/full");
70 
71  // determine which general experiments we have, ie nue or numu or both
72  // allow users to force the use of all event lists if they choose
73  fLoadAllEventLists = pset.get<bool>("LoadAllEventLists", false);
74 
75  this->SetExposures(pset);
76 
77  }
78 
79  // Set the exposures which are configured from the fhicl files
80  //...........................................................................
82  {
83 
84  auto psetExposure = pset.get< std::vector<fhicl::ParameterSet> >("ExposureToUse", std::vector<fhicl::ParameterSet>() );
85  MF_LOG_DEBUG("EventListManipulator")
86  << "There are "
87  << psetExposure.size()
88  << " exposures from fhicl";
89 
90  auto fhcExposureMultiplier = pset.get<double>("FHCExposureMultiplier", 1.);
91  auto rhcExposureMultiplier = pset.get<double>("RHCExposureMultiplier", 1.);
92  fExposure.clear();
93  for(auto const& exposure : psetExposure){
94 
95  auto psetDet = exposure.get<std::string>("Detector" );
96  auto psetPeriod = exposure.get<int >("Period" );
97  auto psetExposure = exposure.get<float >("Exposure" );
98  auto psetLivetime = exposure.get<float >("Livetime" );
99 
100  fPeriodsToUse.insert(psetPeriod);
101 
102  cmf::SpillSummary spillSummary(psetExposure,
103  psetExposure,
104  psetLivetime,
105  0,
106  0);
107 
110 
111  MF_LOG_DEBUG("ELM")
113  << " "
115  << " "
116  << exposure.to_indented_string();
117 
118  // if we want high stats fakes, multiply the exposure by exposure
119  // multiplier - only ever applied to FD
120  if(det == cmf::kFARDET)
121  psetExposure *= (beam == cmf::kFHC) ? fhcExposureMultiplier : rhcExposureMultiplier;
122 
123  auto const& dbsSet = cmf::SelectionUtility::Instance()->SelectionsToUse();
124  for(auto const& sel : (*dbsSet.find(cmf::SelectionUtility::DetBeamSels(det, beam))).sels){
125  cmf::MetaData md(false,
126  det,
128  sel,
130  psetPeriod);
131 
132  if(fExposure.count(md) < 1){
133  fExposure.emplace(md, spillSummary);
134  }
135  else {
136  fExposure[md].totalPOT += psetExposure;
137  fExposure[md].goodPOT += psetExposure;
138  fExposure[md].liveTime += psetLivetime;
139  }
140 
141  } // end loop over selections
142 
143  } // end if we have detectors and selections to make
144 
145  this->PrintPOTSummary (fExposure);
147 
148  //for(auto const& itr : fExposure)
149  // MF_LOG_VERBATIM("EventListManipulator")
150  // << "configuring exposure from fhicl"
151  // << std::setw(30)
152  // << itr.first.ToString()
153  // << " to have "
154  // << itr.second
155  // << " POT and this event list will not be deserialized";
156 
157  }
158 
159  // just print out Exposure summary
160  // TODO - do we want this to work for MINOS data as well?
161  //.....................................................................
163  {
164 
165  std::map<int, double> NDFHC;
166  std::map<int, double> FDFHC;
167  std::map<int, double> NDRHC;
168  std::map<int, double> FDRHC;
169  std::vector<int> ps;
170 
171  for (auto const& itr : exposureMap)
172  {
173  MF_LOG_DEBUG("ELM")
174  << itr.first.ToString()
175  << " "
176  << itr.second;
177 
178  double value;
179  if (isLiveTime) value = itr.second.liveTime;
180  else value = itr.second.goodPOT;
181 
182  // ND FHC
183  if (itr.first.detector == cmf::kNEARDET && itr.first.BeamType() == cmf::BeamType_t::kFHC){
184  if (NDFHC.find(itr.first.period) == NDFHC.end())
185  NDFHC[itr.first.period] = value;
186  if (std::find(ps.begin(), ps.end(), itr.first.period) == ps.end()) ps.push_back(itr.first.period);
187  }
188  // FD FHC
189  if (itr.first.detector == cmf::kFARDET && itr.first.BeamType() == cmf::BeamType_t::kFHC){
190  if (FDFHC.find(itr.first.period) == FDFHC.end())
191  FDFHC[itr.first.period] = value;
192  if (std::find(ps.begin(), ps.end(), itr.first.period) == ps.end()) ps.push_back(itr.first.period);
193  }
194 
195  // ND RHC
196  if (itr.first.detector == cmf::kNEARDET && itr.first.BeamType() == cmf::BeamType_t::kRHC){
197  if (NDRHC.find(itr.first.period) == NDRHC.end())
198  NDRHC[itr.first.period] = value;
199  if (std::find(ps.begin(), ps.end(), itr.first.period) == ps.end()) ps.push_back(itr.first.period);
200  }
201 
202  // FD RHC
203  if (itr.first.detector == cmf::kFARDET && itr.first.BeamType() == cmf::BeamType_t::kRHC){
204  if (FDRHC.find(itr.first.period) == FDRHC.end())
205  FDRHC[itr.first.period] = value;
206  if (std::find(ps.begin(), ps.end(), itr.first.period) == ps.end()) ps.push_back(itr.first.period);
207  }
208  }
209 
210  if (isLiveTime){
211  MF_LOG_VERBATIM("EventListManipulator")
212  << "\n-------------------------------------------------------"
213  << "\nLIVETIME SUMMARY"
214  << "\n-------------------------------------------------------";
215  }
216  else{
217  MF_LOG_VERBATIM("EventListManipulator")
218  << "\n-------------------------------------------------------"
219  << "\nPOT SUMMARY"
220  << "\n-------------------------------------------------------";
221  }
222 
223  std::sort(ps.begin(), ps.end());
224  MF_LOG_VERBATIM("EventListManipulator")
225  << std::setw(20)
226  << ""
227  << std::setw(20)
228  << "ND FHC"
229  << std::setw(20)
230  << "FD FHC"
231  << std::setw(20)
232  << "ND RHC"
233  << std::setw(20)
234  << "FD RHC";
235  for (size_t i = 0; i < ps.size(); ++i){
236  MF_LOG_VERBATIM("EventListManipulator")
237  << std::setw(20)
238  << "Period "+std::to_string(ps[i])
239  << std::setw(20)
240  << NDFHC[ps[i]]
241  << std::setw(20)
242  << FDFHC[ps[i]]
243  << std::setw(20)
244  << NDRHC[ps[i]]
245  << std::setw(20)
246  << FDRHC[ps[i]];
247  }
248 
249  }
250 
251  //----------------------------------------------------------------------------
253  this->PrintExposureSummary(exposureMap, false);
254  }
255 
256  //----------------------------------------------------------------------------
258  this->PrintExposureSummary(exposureMap, true);
259  }
260 
261  // Skip any event list where the selection and detector combination were
262  // not requested in the configuration
263  //......................................................................
265  cmf::DataMC_t const dataMC,
266  std::set<cmf::DetType_t> const& detectors)
267  {
268  if(fLoadAllEventLists) return true;
269 
270  // if we don't want cosmic muons return false right away,
271  // if we do, keep doing checks to be sure this is a selection
272  // we want
274  !fDeserializeCosmics) return false;
275 
276  // check if this is data and we only want MC events
277  // or vice versa
278  if(dataMC == cmf::kMC && !md.isMC){
279  MF_LOG_DEBUG("ELMUEFMMC")
280  << md.ToString()
281  << " failed cmf::kMC and !md.isMC";
282  return false;
283  }
284  else if(dataMC == cmf::kData && md.isMC){
285  MF_LOG_DEBUG("ELMUEFMRealData")
286  << md.ToString()
287  << " failed cmf::kData and md.isMC";
288  return false;
289  }
290  else if(dataMC == cmf::kFakeData && !md.isMC){
291  MF_LOG_DEBUG("ELMUEFMFakeData")
292  << md.ToString()
293  << " failed cmf::kFakeData and !md.isMC";
294  return false;
295  }
296 
297  //Skip unless the epoch is in the fPeriodsToUse set
298  if(fPeriodsToUse.count(md.Period()) < 1){
299  MF_LOG_DEBUG("ELMUEFMPeriod")
300  << md.ToString()
301  << " failed count of periods to use "
302  << md.Period();
303  return false;
304  }
305 
306  // Skip if uncategorized interactions
308  MF_LOG_DEBUG("ELMUEFMInteraction")
309  << md.ToString()
310  << " failed interaction type check to use "
311  << md.interactionType;
312  return false;
313  }
314 
315  // check if the detector for this metadata is desired, return false if it
316  // isn't in the fDetectors set. We have to check detector separately from
317  // the SelectionUtility because the selection utility is set for the entire
318  // job, not just the EventListManipulator which may only want to use
319  // event from a subset of detectors for the job.
320  // Moreover, even though ultimately we want to use both detectors, a job
321  // may want to load one at a time for memory management, so we use the
322  // detectors set to check on that
323  if(fDetectors.count(md.detector) < 1 || detectors.count(md.detector) < 1) return false;
324 
325  // now check that the selection is desired for the detector, beam and
326  // selection of this metadata
328  MF_LOG_DEBUG("ELMUEFMDetector")
329  << md.ToString()
330  << " failed to find "
331  << md.DetectorString()
332  << " "
333  << md.PeriodString()
334  << " "
336  << " in selections to use";
337  return false;
338  }
339 
340  // Metadata types we are going to skip
341  // - ND flux swap for now
342  //if(md->fileType == cmf::kSwap &&
343  // md->detector == cmf::kNEARDET ) return false;
344 
345  return true;
346  }
347 
348  // We want to make sure that the pre-configured cap does not result in
349  // a poor sampling of events, so pass in the number of events in the
350  // tree as well to check how the cap will affect the final number of events
351  // used. If too few events will remain after that cap, up it to ensure
352  //......................................................................
354  long treeEvents)
355  {
356 
357  double cap = (treeEvents > fMaxEventsPerTree - 1) ? fMaxEventsPerTree / (1. * treeEvents) : 1.;
358 
359  MF_LOG_DEBUG("EVentListManipulator")
360  << "Check event cap for "
361  << md.ToString()
362  << ": "
363  << treeEvents
364  << " / "
366  << " / "
368  << " final cap: "
369  << cap;
370 
371  // if no cap is specified, use all the events
372  return cap;
373  }
374 
375  //......................................................................
376  void EventListManipulator::ExtractFromFile(TFile* metadataFile,
377  std::string const& dirName,
378  cmf::DataMC_t const& dataMC,
379  std::set<cmf::DetType_t> const& detectors,
380  std::vector<cmf::FileExtracts> & extractVec)
381  {
382  MF_LOG_DEBUG("EventListManipulator")
383  << "extracting metadata from file";
384 
385  std::string mdString;
386  std::string treeName;
387 
388  TTree *metadataTree = dynamic_cast<TTree *>(metadataFile->Get((dirName + "/metadata").c_str()));
389 
390  if(!metadataTree){
391  throw cet::exception("EventListManipulator")
392  << "cannot find metadata tree: "
393  << dirName + "/metadata";
394  }
395 
396  // The metadata tree may have the same metadata information in it over and over again
397  // because we hadd'ed a bunch of files together. We don't want to read in the event
398  // trees multiple times, so we will make a map of the metadata and spill summary
399  // information so that we only have one entry for each metadata
400  cmf::ExposureMap exposureMap;
401 
402  // loop over the entries in the metadataTree
403  cmf::MetaData *md = nullptr;
404  cmf::SpillSummary *ss = nullptr;
405  metadataTree->SetBranchAddress("metadata", &md);
406  metadataTree->SetBranchAddress("spillsummary", &ss);
407 
408  // We have potentially hadd'ed several files together to get the
409  // eventlisttree, so first loop over the entries in the metadata
410  // tree and sum all the POT values for each metadata type
411  for(int ctr = 0; ctr < metadataTree->GetEntries(); ++ctr){
412  metadataTree->GetEntry(ctr);
413 
414  // now check if we want this metadata for our analysis
415  // and find the minimum POT for all metadata having the same
416  // beam/detector/file/interaction types for concatenated numu or nue
417  // selections
418  if(this->UseEventsFromMetaData(*md, dataMC, detectors)){
419  if(exposureMap.count(*md) < 1) exposureMap[*md] = *ss;
420  else exposureMap[*md] += *ss;
421  }
422  } // end loop over the metadata tree
423 
424  for(auto const& itr : exposureMap){
425 
426  // get the pointer to the event TTree
427  mdString = itr.first.ToString();
428 
429  treeName = dirName + "/" + mdString;
430 
431  if(itr.second.goodPOT > 0 ||
432  (itr.first.fileType == kCosmicBackgroundFile && itr.second.liveTime > 0)){
433  extractVec.emplace_back(itr.first,
434  itr.second,
435  treeName,
436  this->PeriodEventCap(itr.first, dynamic_cast<TTree*>(metadataFile->Get(treeName.c_str()))->GetEntriesFast()));
437  }
438 
439  MF_LOG_DEBUG("EventListManipulator")
440  << "deserializing: "
441  << itr.first.ToString()
442  << " "
443  << itr.second.goodPOT
444  << " "
445  << itr.second.liveTime
446  << " "
447  << exposureMap.size();
448  } // end loop over the metadata to spill summary map
449 
450 
451  // loop over the exposureMap map one more time to check which
452  // selections are being used
453  // for(auto const& itr : exposureMap){
454  // MF_LOG_VERBATIM("EventListManipulator")
455  // << "Using "
456  // << itr.first.ToString()
457  // << " in the analysis";
458  // }
459 
460  }
461 
462  // Use this method to load the TTrees containing the cmf::EventLists
463  // until we are able to store them directly in an art file
464  //
465  // If onlyMC is true, only load the MC lists from the file
466  //......................................................................
468  cmf::DataMC_t dataMC,
469  std::set<cmf::DetType_t> const& detectors)
470  {
471  cmf::DataMC_t dataType = (dataMC == cmf::kData) ? dataMC : cmf::kMC;
472 
473  std::vector<cmf::FileExtracts> extractVec;
474  size_t totalEvents = 0;
475 
478  double normalization = 1.;
479 
480  // each input file is assumed to have a metadata tree and a cmfEvent tree
481  for(auto & fileName : fCMFEventLabels){
482 
483  TFile *f = TFile::Open(fileName.c_str());
484 
485  if(f->IsOpen())
486  MF_LOG_VERBATIM("EventListManipulator")
487  << "successfully opened "
488  << fileName;
489  else
490  throw cet::exception("EventListManipulator")
491  << "Cannot open file "
492  << fileName;
493 
494  for(auto const& dirName : fTreeDirectories){
495 
496  this->ExtractFromFile(f, dirName, dataType, detectors, extractVec);
497 
498  MF_LOG_DEBUG("EventListManipulator")
499  << "There are "
500  << extractVec.size()
501  << " metadata categories";
502  } // end loop over directories
503 
504  // map the detector/file/beam/interaction type + general selection type to the minimum exposure for that
505  // combination. Use the ratio between the capped exposure and the found minimum to adjust the cap for
506  // the other lists making up the general selection.
507  // That is, new cap = (capped exposure for sub list) / (capped exposure for minimum exposure of sub lists)
508  std::map<long, double> keyToMaxExposure;
509  long key;
510  double cappedExposure;
511 
512  for(auto & itr : extractVec){
513  key = itr.Key();
514 
515  // don't waste precision on the exponent for the POT
516  cappedExposure = (itr.metadata.IsCosmicMuon()) ? itr.CappedLiveTime() : itr.CappedPOT() * 1.e-20;
517 
518  if(keyToMaxExposure.count(key) < 1)
519  keyToMaxExposure[key] = cappedExposure;
520  else
521  keyToMaxExposure[key] = std::max(keyToMaxExposure.find(key)->second, cappedExposure);
522 
523  MF_LOG_DEBUG("EventListManipulator")
524  << "input key: "
525  << cmf::KeyToString(itr.metadata.Key())
526  << " concat key: "
527  << cmf::KeyToString(key)
528  << " uncapped POT: "
529  << itr.spillSummary.goodPOT
530  << " uncapped LiveTime: "
531  << itr.spillSummary.liveTime
532  << " / capped exposure (x 1e20) "
533  << keyToMaxExposure[key];
534  }
535 
536  // now loop over the extracted information again to fill the event lists
537  // cosmic ray background livetime normalization is done when the spectra are
538  // filled
539  double exposure;
540  for(auto & itr : extractVec){
541 
542  normalization = 1.;
543 
544  // need a non-const spill summary for filling the event list in case we have to
545  // apply a cap. Any adjustment of the exposure due to capping will be done in the
546  // FillEventList method
547  ss = itr.spillSummary;
548 
549  // check if we are concatenating either the numu or nue selections
550  // and make the appropriate metadata for that case. We need to also
551  // adjust the event cap for this particular tree and change the
552  // total POT for this selection to account for normalization differences
553  // between the selections to be concatenated
554  md = itr.metadata;
555 
556  exposure = (md.IsCosmicMuon()) ? itr.CappedLiveTime() : itr.CappedPOT() * 1.e-20;
557 
560  md = cmf::MetaData(md.isMC,
561  md.detector,
562  md.fileType,
564  md.interactionType,
565  md.Period());
566 
567  normalization = keyToMaxExposure[itr.Key()] / exposure;
568  }
569  else if(cmf::SelectionUtility::Instance()->UsesSelection(cmf::kNuESelection) &&
571  md = cmf::MetaData(md.isMC,
572  md.detector,
573  md.fileType,
575  md.interactionType,
576  md.Period());
577 
578  normalization = keyToMaxExposure[itr.Key()] / exposure;
579  }
580 
581  //if(md.IsNuMuSelected())
582  MF_LOG_DEBUG("EventListManipulator")
583  << "Filling "
584  << itr.metadata.ToString()
585  << " / "
586  << md.ToString()
587  << " with cap "
588  << itr.eventCap * normalization
589  << " Exposure "
590  << exposure * normalization
591  //<< " : Events "
592  //<< itr.eventTree->GetEntries()
593  << " uncapped POT/livetime "
594  << ss.goodPOT
595  << " / "
596  << ss.liveTime;
597 
598  ss *= itr.eventCap * normalization;
599 
600  auto listItr = cmf::FindEventList(md, eventLists);
601  if(listItr == eventLists.end()){
602  eventLists.emplace_back(cmf::EventList(md, ss));
603  listItr = cmf::FindEventList(md, eventLists);
604  }
605 
606  this->FillEventList(f, itr.treeName, itr.eventCap * normalization, *listItr);
607 
608  totalEvents += listItr->TotalEvents();
609  } // end loop over extracted trees
610 
611  f->Close();
612  }// end loop over files
613 
614  MF_LOG_VERBATIM("EventListManipulator")
615  << "There are "
616  << eventLists.size()
617  << " "
618  << cmf::cDataTypeStrings[dataMC]
619  << " event lists and "
620  << totalEvents
621  << " weighted events";
622 
623  std::map<cmf::MetaData, double> totalPOTPerSel;
624 
625  for(auto const& itr : eventLists){
626 
627  md = itr.ListMetaData();
628  ss = itr.ListSpillSummary();
629 
630  if (md.fileType == kCosmicBackgroundFile){
631  MF_LOG_DEBUG("EventListManipulator")
632  << std::setw(55)
633  << md.ToString()
634  << " of period "
635  << md.Period()
636  << " "
637  << std::setw(15)
638  << " represents "
639  << ss.liveTime
640  << " livetime and "
641  << itr.TotalEvents()
642  << " events "
643  << itr.TotalEvents() / ss.liveTime;
644  }
645  else{
646  MF_LOG_DEBUG("EventListManipulator")
647  << std::setw(55)
648  << " "
649  << std::setw(15)
650  << " represents "
651  << ss.goodPOT
652  << " POT and "
653  << itr.TotalEvents()
654  << " events "
655  << itr.TotalEvents() / ss.goodPOT;
656  }
657 
661  else if (md.selectionType == cmf::kNCSelection) selType = cmf::kNCSelection;
662  else continue;
663 
664  cmf::MetaData mdtmp (md.MCKey(),
665  md.detector,
666  md.fileType,
667  selType,
669  md.BeamType() == kFHC ? 1 : 4);
670 
671  if (md.interactionType != cmf::kNuMuCC) continue;
672 
673  if (md.fileType != cmf::kBeam &&
674  md.fileType != cmf::kSwap &&
675  md.fileType != cmf::kTauSwap) continue;
676 
677  if (totalPOTPerSel.find(mdtmp) == totalPOTPerSel.end())
678  totalPOTPerSel[mdtmp] = ss.goodPOT;
679  else
680  totalPOTPerSel[mdtmp] += ss.goodPOT;
681 
682  }
683 
684  for (auto const& itr : totalPOTPerSel){
685  if (itr.first.period == 1 || itr.first.period == 4){
686  MF_LOG_VERBATIM("EventListManipulator")
687  << "-- TOTAL "
688  << cmf::cBeamType_Strings[itr.first.BeamType()]
689  << " POT: \n"
690  << cDetType_Strings[itr.first.detector]
691  << " "
692  << cFileTypeStrings[itr.first.fileType]
693  << " "
694  << cSelectionType_Strings[itr.first.selectionType]
695  << " POT: "
696  << itr.second;
697  }
698  }
699 
700  }
701 
702  //--------------------------------------------------------------------------------------------
704  cmf::EventId & ev,
705  cmf::DataVars & dv,
706  cmf::TruthVars & tv,
707  cmf::WeightVars & wv)
708  {
709 
711  fileName += cBeamType_Strings[md.BeamType()];
712  fileName += cSelectionType_Strings[md.selectionType];
713  fileName += cFileTypeStrings[md.fileType];
714  fileName += ".txt";
715 
716  std::ifstream ifile;
717  std::ofstream ofile;
718  ifile.open(fileName);
719  if(!ifile) {
720  ofile.open(fileName, std::ofstream::out);
721  ofile
722  << "\n"
723  //<< md.ToString()
724  //<< " "
725  << ev.run
726  << " "
727  << ev.subrun
728  << " "
729  << ev.event
730  << " "
731  << ev.cycle
732  << " "
733  << ev.slice
734  << " "
735  << dv.fLep_RecoE
736  << " "
737  << dv.fHad_RecoE
738  << " "
739  << dv.fLep_RecoE + dv.fHad_RecoE
740  << " "
741  << dv.fNu_RecoE
742  << " "
743  << cmf::RecoEnergy(dv.fLep_RecoE, dv.fHad_RecoE, md)
744  << " ";
745  if (md.isMC){
746  ofile
747  << tv.fTrueE
748  << " "
749  << tv.fTruePDG
750  << " "
751  << tv.fTrueCCNC
752  << " "
753  << wv.fXSecCV2020_Weight
754  << " "
755  << wv.fPPFXFluxCV_Weight;
756  }
757 
758 
759  //ofile
760  // << "run"
761  // << " subrun"
762  // << " ev"
763  // << " cyc"
764  // << " slc"
765  // << " Ereco";
766  // if (md.isMC){
767  // ofile
768  // << " Etrue"
769  // << " pdg"
770  // << " ccnc"
771  // << " cvwgt";
772  // }
773  }
774  ofile.open(fileName, std::ofstream::out | std::ofstream::app);
775 
776  ofile
777  << "\n"
778  //<< md.ToString()
779  //<< " "
780  << ev.run
781  << " "
782  << ev.subrun
783  << " "
784  << ev.event
785  << " "
786  << ev.cycle
787  << " "
788  << ev.slice
789  << " "
790  << dv.fLep_RecoE
791  << " "
792  << dv.fHad_RecoE
793  << " "
794  << dv.fHad_RecoE + dv.fLep_RecoE
795  << " "
796  << dv.fNu_RecoE
797  << " "
798  << cmf::RecoEnergy(dv.fLep_RecoE, dv.fHad_RecoE, md)
799  << " ";
800  if (md.isMC){
801  ofile
802  << tv.fTrueE
803  << " "
804  << tv.fTruePDG
805  << " "
806  << tv.fTrueCCNC
807  << " "
808  << wv.fXSecCV2020_Weight
809  << " "
810  << wv.fPPFXFluxCV_Weight;
811  }
812 
813  ofile.close();
814  }
815 
816  //......................................................................
817  void EventListManipulator::FillEventList(TFile * eventFile,
818  std::string const& treeName,
819  double const& eventCap,
820  cmf::EventList & eventList)
821  {
822  cmf::EventId eventId;
823  cmf::DataVars dataVars;
824  cmf::TruthVars truthVars;
825  cmf::WeightVars weightVars;
826 
827  auto const& md = eventList.ListMetaData();
828 
829  std::vector<unsigned char>* systVarKeys = nullptr; // PROD5
830  std::vector<float>* systVarMinus2Wgts = nullptr; // PROD5
831  std::vector<float>* systVarMinus1Wgts = nullptr; // PROD5
832  std::vector<float>* systVarPlus1Wgts = nullptr; // PROD5
833  std::vector<float>* systVarPlus2Wgts = nullptr; // PROD5
834 
835  auto const sysKeys = cmf::ParameterUtility::Instance()->SysParKeys();
836  // for(auto const& itr : sysKeys)
837  // MF_LOG_VERBATIM("EventListManipulator")
838  // << itr
839  // << " "
840  // << cmf::KeyToVarName(itr);
841 
842  TTree *eventTree = dynamic_cast<TTree*>(eventFile->Get(treeName.c_str()));
843  eventTree->SetBranchAddress("dataVars", &dataVars);
844  if(fUseEventId)
845  eventTree->SetBranchAddress("eventId", &eventId);
846 
847  if(md.isMC){
848  eventTree->SetBranchAddress("truthVars" , &truthVars );
849  eventTree->SetBranchAddress("weightVars", &weightVars);
850 
851  eventTree->SetBranchAddress("systVarKeys", &systVarKeys);
852  eventTree->SetBranchAddress("systVarMinus2Wgts", &systVarMinus2Wgts);
853  eventTree->SetBranchAddress("systVarMinus1Wgts", &systVarMinus1Wgts);
854  eventTree->SetBranchAddress("systVarPlus1Wgts", &systVarPlus1Wgts);
855  eventTree->SetBranchAddress("systVarPlus2Wgts", &systVarPlus2Wgts);
856  } // end if MC
857 
858  // get the generic identifier for numu or nue selection so
859  // we can then determine the maximum energy of the neutrino
860  // events to use
862  if ( md.IsNuESelected() ) selection = cmf::kNuESelection;
863  else if( md.IsNCSelected() ) selection = cmf::kNCSelection;
864 
865  MF_LOG_DEBUG("EventListManipulator")
866  << md.ToString()
867  << " selection is "
868  << selection
869  << " max energy: "
870  << fMaxNuEnergy.find(selection)->second;
871 
872  // check that our event list corresponds to a reasonable exposure
873  if(eventList.ListSpillSummary().goodPOT < 1.e17 && md.fileType != kCosmicBackgroundFile)
874  throw cet::exception("EventListManipulator")
875  << "exposure for "
876  << eventList.ListMetaData().ToString()
877  << " is "
878  << eventList.ListSpillSummary().goodPOT
879  << " < 10^9 POT. That cannot be right, bail";
880 
881  auto systs = std::make_unique<cmf::SystVarColl>();
882  uint8_t loc;
883 
884  MF_LOG_DEBUG("ELM")
885  << "There are "
886  << eventTree->GetEntries()
887  << " events in tree before capping, "
888  << eventList.size()
889  << " events in list already, expect "
890  << eventTree->GetEntries() * eventCap + eventList.size()
891  << " events in list after capping";
892 
893  // use maxEvents to adjust the number of entries loaded in testing
894  long maxEvents = eventTree->GetEntries();
895  for(long evNum = 0; evNum < maxEvents; ++evNum){
896 
897  // If the uniformly distributed random number is greater than the
898  // cap forget the current event. This has the correct boundry
899  // conditions - if the cap is 1, then all events are used; if it is
900  // 0 then none are.
901 
902  // At some point we may need to worry about keeping more events in the
903  // tails of the spectrum, but for now it looks like we can keep at least
904  // 80% of all MetaData options, so we should not be biasing the tails
905  if(fRandom.Rndm() > eventCap) continue;
906 
907  eventTree->GetEntry(evNum);
908 
909  if (fFillTextFile)
910  this->FillTextFile(md, eventId, dataVars, truthVars, weightVars);
911 
912  // don't bother filling events outside the energy range for
913  // this selection - it just wastes time in the fitting
914  if(dataVars.fNu_RecoE > fMaxNuEnergy.find(selection)->second ||
916 
917  // set the reco energy bin for this event
919  dataVars.fRecoEBin += cmf::CovarianceBinUtility::Instance()->KeyToOffset(md.DetectorBeamSelectionKey());
920 
921  // now set the dataVals
922  cmf::DataVarVals dataVals(dataVars);
923 
924  // now we can fill the rest of the event info because we are
925  // keeping this event
926  if(md.isMC){
927  systs->clear();
928 
929  // if(systVarKeys->size() < 1)
930  // MF_LOG_VERBATIM("EventListManipulator")
931  // << md.ToString();
932 
933  for(size_t v = 0; v < systVarKeys->size(); ++v){
934  loc = (*systVarKeys)[v];
935 
936  // MF_LOG_VERBATIM("EventListManipulator")
937  // << (*systVarKeys)[v]
938  // << " "
939  // << cmf::KeyToVarName(loc);
940 
941  // check if this sytematic is something we want to evaluate in the job,
942  // ignore it if not
943  if(sysKeys.find(loc) == sysKeys.end()) continue;
944 
945  // don't bother saving a bunch of 1s in the SystVarColl, instead just
946  // return all 1s if that systematic is desired for this event from MCVars.
947  if((*systVarMinus2Wgts)[v] == 1 &&
948  (*systVarMinus1Wgts)[v] == 1 &&
949  (*systVarPlus1Wgts )[v] == 1 &&
950  (*systVarPlus2Wgts )[v] == 1 ) continue;
951 
952  systs->emplace_back(loc,
953  (*systVarMinus2Wgts)[v],
954  (*systVarMinus1Wgts)[v],
955  (*systVarPlus1Wgts )[v],
956  (*systVarPlus2Wgts )[v]);
957 
958  MF_LOG_DEBUG("ELMVarCheck")
959  << "check systematic vector sizes for event "
960  << evNum
961  << ": "
962  << systVarKeys->size()
963  << " "
964  << systVarMinus2Wgts->size()
965  << " "
966  << systVarMinus1Wgts->size()
967  << " "
968  << systVarPlus1Wgts->size()
969  << " "
970  << systVarPlus2Wgts->size()
971  << " "
972  << systs->size()
973  << " "
974  << systs->back();
975  }
976 
977  if(fUseEventId)
978  eventList.AddEvent(eventId,
979  dataVals,
980  cmf::MCVarVals(truthVars, weightVars, *systs));
981  else
982  eventList.AddEvent(dataVals,
983  cmf::MCVarVals(truthVars, weightVars, *systs));
984  } // end if md.isMC
985  else{
986 
987  // Print out the event id information for each data event in the NOvA FD
988  if(md.detector == cmf::kFARDET){
989 
990  std::string selectionStr("numu");
991  if (md.IsNuESelected()) selectionStr = "nue";
992  else if(md.IsNCSelected() ) selectionStr = "nc";
993  MF_LOG_DEBUG("EventListManipulator")
994  << selectionStr
995  << " "
996  << eventId;
997 
998  } // end if FD to print the event information
999 
1000  if(fUseEventId)
1001  eventList.AddEvent(eventId,
1002  dataVals);
1003  else
1004  eventList.AddEvent(dataVals);
1005  }
1006 
1007  } // end loop over events
1008 
1009  // loop over the events we just added and check the systematics
1010  // were loaded correctly
1011  // for(auto const& evt : eventList){
1012  // MF_LOG_VERBATIM("EventListManipulator")
1013  // << "event has systematic uncertainties with values ";
1014  // for(auto const& sitr : evt->MCVals().SystematicShifts())
1015  // MF_LOG_VERBATIM("EventListManipulator")
1016  // << "\t"
1017  // << sitr;
1018  // }
1019 
1020  systs.reset();
1021  eventTree = nullptr;
1022  } // done
1023 
1024 } // end namespace
T max(const caf::Proxy< T > &a, T b)
std::map< cmf::MetaData, cmf::SpillSummary > ExposureMap
Definition: Structs.h:215
cmf::SpillSummary const & ListSpillSummary() const
Definition: Event.h:133
fileName
Definition: plotROC.py:78
int KeyToOffset(long const &key, bool allSels=false)
std::vector< SystGroupDef > systs
Definition: syst_header.h:385
DetBeamSelSet const & SelectionsToUse()
const std::vector< std::string > cFileTypeStrings({"Beam","FluxSwap","TauSwap","Data","CosmicBackground","RockFluxSwap","RockNonSwap","num_file_types","bad_file_type","UnknownFileType"})
long Period() const
Definition: Structs.h:128
cmf::DetType_t detector
Definition: Structs.h:114
float fPPFXFluxCV_Weight
For 2020 Ana, store weights seperately.
Definition: VarVals.h:161
float fTrueCCNC
true cc vs nc
Definition: VarVals.h:44
float fTruePDG
true PDG
Definition: VarVals.h:43
double liveTime
Definition: Structs.h:187
const char * p
Definition: xmltok.h:285
const std::vector< std::string > cDetType_Strings({"UnknownDet","NearDet","FarDet","MINOSNear","MINOSFar","AllDetectors"})
static std::vector< std::string > cDataTypeStrings({"Data","MC","FakeData","Data and MC"})
double const & SelectionLowEdge(cmf::MetaData const &md)
bool fFillTextFile
fills text file with selection information
Float_t ss
Definition: plot.C:24
static SelectionUtility * Instance()
std::vector< std::string > fCMFEventLabels
Labels in input files holding CovarianceMatrixFit Events.
std::set< long > fPeriodsToUse
which periods to use in the analysis
cmf::EventListColl::iterator FindEventList(cmf::MetaData const &md, cmf::EventListColl &eventListColl)
Definition: Event.cxx:44
std::set< uint8_t > SysParKeys() const
enum cmf::det_type DetType_t
std::string PeriodString() const
Definition: Structs.h:162
::xsd::cxx::tree::exception< char > exception
Definition: Database.h:225
float fNu_RecoE
Definition: VarVals.h:127
enum cmf::sel_type SelectionType_t
std::map< cmf::SelectionType_t, double > fMaxNuEnergy
maximum neutrino energy to go into the lists
float fMaxEventsPerTree
maximum number of events to use per tree
long MCKey() const
Definition: Structs.h:120
void FillEventList(TFile *eventFile, std::string const &treeName, double const &eventCap, cmf::EventList &eventList)
std::set< cmf::DetType_t > fDetectors
which detector(s) are we loading events from
static bool IsNuESelected(cmf::SelectionType_t const &sel)
Definition: StaticFuncs.h:380
static bool IsNuMuQuantiles(cmf::SelectionType_t const &sel)
Definition: StaticFuncs.h:373
float fTrueE
True nu energy.
Definition: VarVals.h:42
EventListManipulator(fhicl::ParameterSet const &pset)
cmf::ExposureMap fExposure
POT in 1e12 to normalise to.
float fRecoEBin
logical energy bin
Definition: VarVals.h:133
bool fLoadAllEventLists
force all event lists to be loaded
Base container for the MC related Vars that constitute an event.
Definition: VarVals.h:229
std::string DetectorString() const
Definition: Structs.h:151
int EnergyToBin(double const &energy, cmf::MetaData const &md)
void Deserialize(cmf::EventListColl &eventLists, cmf::DataMC_t dataMC=cmf::kBoth, std::set< cmf::DetType_t > const &detectors=std::set< cmf::DetType_t >({cmf::kNEARDET, cmf::kFARDET}))
static ParameterUtility * Instance()
int subrun
Definition: VarVals.h:80
void PrintLivetimeSummary(cmf::ExposureMap &exposureMap)
cmf::SelectionType_t selectionType
Definition: Structs.h:116
if(dump)
const XML_Char int const XML_Char * value
Definition: expat.h:331
const Cut sels[kNumSels]
Definition: vars.h:44
int cycle
Definition: VarVals.h:82
enum cmf::dataMC_type DataMC_t
T get(std::string const &key) const
Definition: ParameterSet.h:231
cmf::MetaData const & ListMetaData() const
Definition: Event.h:138
static cmf::DetType_t StringToDetectorType(std::string const &str)
Definition: StaticFuncs.h:267
float fMinEventsPerTree
minimum number of events to use per tree
const std::vector< std::string > cSelectionType_Strings({"NuESel_AllPID","NuESel_LowPID","NuESel_MidPID","NuESel_HighPID","NuESel_Peripheral","NuMuSel","NuMuSelQ1","NuMuSelQ2","NuMuSelQ3","NuMuSelQ4","NCSel","UnknownSel"})
static cmf::BeamType_t PeriodToBeamType(std::string const &str)
double PeriodEventCap(cmf::MetaData const &md, long treeEvents)
const std::vector< std::string > cBeamType_Strings({"FHC","RHC","0HC","UnknownBeam"})
float fXSecCV2020_Weight
For 2020 Ana, store weights seperately.
Definition: VarVals.h:160
void FillTextFile(cmf::MetaData const &md, cmf::EventId &ev, cmf::DataVars &dv, cmf::TruthVars &tv, cmf::WeightVars &wv)
std::map< std::string, float > fEventCaps
int slice
Definition: VarVals.h:83
Base container for the MC related Vars that constitute an event.
Definition: VarVals.h:270
bool fDeserializeCosmics
whether or not deserialize cosmic muons
associates each detector-beam pair to a set of selections
std::string ToString() const
Definition: Structs.cxx:135
Module to combine a set of results into a single file currently only does one data product type at a ...
Definition: Event.cxx:24
std::vector< cmf::EventList > EventListColl
Definition: Event.h:150
cmf::FileType_t fileType
Definition: Structs.h:115
double goodPOT
Definition: Structs.h:186
std::string dirName
Definition: PlotSpectra.h:47
void PrintExposureSummary(cmf::ExposureMap &exposureMap, bool isLivetime)
static float RecoEnergy(float leptonicEnergy, float hadronicEnergy, cmf::MetaData const &md)
Definition: VarVals.h:323
TRandom3 fRandom
Random number generator to use creating fake data lists.
#define MF_LOG_VERBATIM(category)
void ExtractFromFile(TFile *metadataFile, std::string const &dirName, cmf::DataMC_t const &dataMC, std::set< cmf::DetType_t > const &detectors, std::vector< cmf::FileExtracts > &extractVec)
app
Definition: demo.py:189
#define MF_LOG_DEBUG(id)
enum cmf::beam_type BeamType_t
bool UseEventsFromMetaData(cmf::MetaData const &md, cmf::DataMC_t dataMC, std::set< cmf::DetType_t > const &detectors)
void AddEvent(EventPtr const &ev)
Definition: Event.h:165
void reconfigure(const fhicl::ParameterSet &p)
void SetExposures(fhicl::ParameterSet const &pset)
void PrintPOTSummary(cmf::ExposureMap &exposureMap)
cmf::BeamType_t BeamType() const
Definition: Structs.cxx:184
float fHad_RecoE
Definition: VarVals.h:128
cmf::InteractionType_t interactionType
Definition: Structs.h:117
bool fUseEventId
use if you need access to run/subrun/event/etc
std::vector< std::string > fTreeDirectories
directory holding the input trees
std::string to_string(ModuleType const mt)
Definition: ModuleType.h:34
bool IsCosmicMuon() const
Definition: Structs.h:139
int event
Definition: VarVals.h:81
std::size_t size() const
Definition: Event.h:116
static std::string KeyToString(long const &key)
Definition: StaticFuncs.h:227
static CovarianceBinUtility * Instance()
float fLep_RecoE
Definition: VarVals.h:129
enum BeamMode string
bool UsesDetBeamAndSel(cmf::MetaData const &md) const