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