Public Member Functions | Private Member Functions | Private Attributes | List of all members
fnex::EventListManipulator Class Reference

#include "/cvmfs/nova-development.opensciencegrid.org/novasoft/releases/N20-12-02/FNEX/modules/EventListManipulator.h"

Public Member Functions

 EventListManipulator (fhicl::ParameterSet const &pset, fnex::InputPoint const &initialGuess)
 
virtual ~EventListManipulator ()
 
void reconfigure (const fhicl::ParameterSet &p)
 
fnex::EventListMap Deserialize (fnex::DataMC_t dataMC=fnex::kBoth)
 
fnex::EventListMap DeserializeData ()
 
fnex::EventListMap DeserializeMC ()
 

Private Member Functions

void ExtractMetaDataFromFile (TFile *metadataFile, std::string const &dirName, fnex::DataMC_t const &dataMC, std::map< fnex::MetaData, fnex::SpillSummary > &mdToSS)
 
void ExtractEventsFromFile (TFile *eventFile, std::string const &dirName, fnex::DataMC_t const &dataMC, std::pair< fnex::MetaData, fnex::SpillSummary > const &mdToSS, fnex::EventListMap &eventLists)
 
void CheckDataMCListConsistency (std::set< fnex::MetaDataLite > const &dataSets, std::set< fnex::MetaDataLite > const &mcSets, fnex::DataMC_t const &dataMC, fnex::EventListMap &eventLists)
 
void CreateFarDetCosmicBackgrounds (fnex::EventListMap &eventLists)
 
void LoadVarVals (fnex::Event &fnexEvent, bool isMC, std::vector< std::string > const &indepVarNames, std::vector< float > const &indepVarVals)
 
fnex::EventList FillEventList (TTree *eventTree, fnex::MetaData const &md, fnex::SpillSummary &ss)
 
void CheckMCListsForFakeDataEpochs (fnex::EventListMap &lists, novadaq::cnv::DetId const &det)
 
void FillFakeDataLists (fnex::EventListMap &lists)
 
void FillFakeDataListsForDetector (fnex::EventListMap &lists, novadaq::cnv::DetId const &det, std::map< int, double > const &detEpochToMinPOT)
 
fnex::EventList AddFakeDataToEventList (fnex::MetaData const &md, fnex::SpillSummary const &ss, fnex::EventList const &mcList)
 
fnex::EventList CreateFarDetCosmicBackgroundsFromEventList (fnex::MetaData const &md, fnex::SpillSummary const &ss, fnex::EventListMap const &listmap)
 
fnex::EventList CreateFarDetCosmicBackgroundsFromHist (fnex::MetaData const &md, fnex::SpillSummary const &ss)
 
fnex::EventList CreatePoissonFakeData (fnex::MetaData const &md, fnex::SpillSummary const &ss, fnex::EventList const &origList)
 
void MakeHistograms (fnex::EventListMap const &lists, std::string const &dirName)
 
void FillHistograms (fnex::EventListMap const &evLists, std::string const &dirName)
 
void FillHistogram (art::TFileDirectory &tfdir, fnex::EventList const &evList, std::string const &listType, double const &normalization)
 
void ConcatenateLists (fnex::EventListMap const &lists, std::string const &dirName, bool const &byEpoch, bool const &byInteraction, bool const &byFile)
 
double RescalePOTForEventCap (size_t &totalEvents, size_t const &eventCap)
 
double RescalePOTForEventCap (double &totalEvents, size_t const &eventCap)
 
void FixNearDetMCExposures (fnex::EventListMap &lists)
 
std::unique_ptr< TH1F > FarDetCosmicBackgroundHistAndScale (fnex::MetaData const &md, float &scale, fnex::SpillSummary const &ss)
 
float EpochEventCap (fnex::MetaData const &md, int treeEvents)
 
void WriteCappedListsAsTrees (fnex::EventListMap const &list)
 
void DetermineSelections (std::vector< fhicl::ParameterSet > const &selections)
 
bool UseEventsFromMetaData (fnex::MetaData const &md, fnex::DataMC_t const dataMC)
 
void CreateNDPeripheralLists (fnex::EventListMap &listmap)
 
std::unique_ptr< TH1F > Cosmics2017 (fnex::MetaData const &md)
 
std::unique_ptr< TH1F > Cosmics2018 (fnex::MetaData const &md, TFile &tf)
 Get the nue 2018 CAFAna cosmic histogram. More...
 

Private Attributes

art::ServiceHandle< art::TFileServicefTFS
 TFileService. More...
 
std::vector< std::stringfFNEXEventLabels
 Labels in input files holding FNEX Events. More...
 
std::vector< std::stringfTreeDirectories
 directory holding the input trees More...
 
TRandom3 fRandom
 Random number generator to use creating fake data lists. More...
 
EnergyHistMap fNeutrinoEnergy
 neutrino energy More...
 
ExposureMap fFakeExposure
 POT in 1e12 of an exposure to fake up. More...
 
fnex::InputPoint fFakeDataPoint
 the fake data input point More...
 
fnex::InputPoint fInitialGuess
 the fake data input point More...
 
std::map< long, float > fEventCaps
 
bool fPoissonFakeData
 make the fake data Poisson distributed More...
 
bool fUseCosmicHists
 
bool fUseEventId
 use if you need access to run/subrun/event/etc More...
 
bool fSumSubEpochs
 do we want to sum all sub epochs into a single one More...
 
std::unique_ptr< CorrelationsfCorrelations
 
std::string fNuMuCosmicHistFile
 file containing the cosmic background histograms More...
 
std::string fNuECosmicHistFile
 file containing the cosmic background histograms More...
 
std::string fNuEDataHistFile
 file containing the data live time for cosmic normalization More...
 
std::string fCalibHistFile
 file containing calibration systematic histograms More...
 
double fTotalFDDataPOTFHC
 total FD data POT FHC More...
 
double fTotalFDDataPOTRHC
 total FD data POT RHC More...
 
double fTotalNDDataPOTFHC
 total ND data POT FHC More...
 
double fTotalNDDataPOTRHC
 total ND data POT RHC More...
 
std::map< std::string, fnex::SpillSummaryfFDPOTSum
 data POT per epoch in the FD More...
 
std::set< intfPeriodsToUse
 which periods to use in the analysis More...
 
std::map< novadaq::cnv::DetId, std::set< fnex::SelectionType_t > > fSelectionsToUse
 which selection types to use? More...
 
std::map< fnex::SelectionType_t, double > fMaxNuEnergy
 maximum neutrino energy to go into the lists More...
 
std::string fAnalysisYear
 analysis year, "2016" for SA and "2017" for TA More...
 
bool fMakeNDPeripheral
 
bool fLoadAllEventLists
 force all event lists to be loaded More...
 
float fMaxEventsPerTree
 maximum number of events to use per tree More...
 
bool fUseExtrapolation
 are we extrapolating or using covariance matrices More...
 

Detailed Description

Definition at line 48 of file EventListManipulator.h.

Constructor & Destructor Documentation

fnex::EventListManipulator::EventListManipulator ( fhicl::ParameterSet const &  pset,
fnex::InputPoint const &  initialGuess 
)

Definition at line 27 of file EventListManipulator.cxx.

References fEventCaps, and reconfigure().

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  }
double fTotalNDDataPOTFHC
total ND data POT FHC
const char * p
Definition: xmltok.h:285
double fTotalFDDataPOTRHC
total FD data POT RHC
double fTotalFDDataPOTFHC
total FD data POT FHC
std::map< long, float > fEventCaps
fnex::InputPoint fInitialGuess
the fake data input point
double fTotalNDDataPOTRHC
total ND data POT RHC
void reconfigure(const fhicl::ParameterSet &p)
fnex::EventListManipulator::~EventListManipulator ( )
virtual

Definition at line 43 of file EventListManipulator.cxx.

44  {
45  return;
46  }

Member Function Documentation

fnex::EventList fnex::EventListManipulator::AddFakeDataToEventList ( fnex::MetaData const &  md,
fnex::SpillSummary const &  ss,
fnex::EventList const &  mcList 
)
private

Definition at line 1144 of file EventListManipulator.cxx.

References fnex::EventList::AddEvent(), fUseExtrapolation, fnex::SpillSummary::goodPOT, fnex::ShifterAndWeighter::Instance(), fnex::MetaData::IsNuESelected(), fnex::MetaData::IsNuMuSelected(), fnex::kAllWgt, fnex::kFake_Weight, fnex::kHad_RecoE, fnex::kLep_RecoE, fnex::kSystAndOscWgt, fnex::EventList::ListMetaData(), fnex::EventList::ListSpillSummary(), LOG_DEBUG, LOG_VERBATIM, fnex::ShifterAndWeighter::ShiftAmount(), fnex::EventList::size(), fnex::MetaData::ToString(), fnex::ShifterAndWeighter::Weight(), and wgt.

Referenced by FillFakeDataListsForDetector().

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  }
#define LOG_DEBUG(stream)
Definition: Messenger.h:149
enum fnex::weight_type WeightType_t
Base container for the MC related Vars that constitute an event.
Definition: VarVals.h:930
Float_t ss
Definition: plot.C:24
static const unsigned char kFake_Weight
Definition: VarVals.h:39
static ShifterAndWeighter * Instance()
double Weight(fnex::EventPtr const &curEvent, fnex::MetaData const &md, fnex::WeightType_t const wgtType=kAllWgt)
double ShiftAmount(unsigned char const &param)
const ana::Var wgt
std::shared_ptr< Event > EventPtr
Definition: Event.h:50
bool fUseExtrapolation
are we extrapolating or using covariance matrices
static const unsigned char kHad_RecoE
Definition: VarVals.h:34
#define LOG_VERBATIM(category)
static const unsigned char kLep_RecoE
Definition: VarVals.h:35
Base container for the MC related Vars that constitute an event.
Definition: VarVals.h:874
void fnex::EventListManipulator::CheckDataMCListConsistency ( std::set< fnex::MetaDataLite > const &  dataSets,
std::set< fnex::MetaDataLite > const &  mcSets,
fnex::DataMC_t const &  dataMC,
fnex::EventListMap eventLists 
)
private

Definition at line 611 of file EventListManipulator.cxx.

References fFDPOTSum, fnex::kData, fnex::kDataFile, fnex::kFakeData, fnex::kMC, fnex::kUnknownInteraction, LOG_VERBATIM, fetch_tb_beamline_files::md, and moon_position_table_new3::second.

Referenced by Deserialize().

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  }
::xsd::cxx::tree::exception< char > exception
Definition: Database.h:225
std::map< std::string, fnex::SpillSummary > fFDPOTSum
data POT per epoch in the FD
#define LOG_VERBATIM(category)
void fnex::EventListManipulator::CheckMCListsForFakeDataEpochs ( fnex::EventListMap lists,
novadaq::cnv::DetId const &  det 
)
private

Definition at line 1035 of file EventListManipulator.cxx.

References fnex::MetaData::epoch, fFakeExposure, LOG_DEBUG, fetch_tb_beamline_files::md, and fnex::MetaData::ToString().

Referenced by FillFakeDataLists().

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  }
#define LOG_DEBUG(stream)
Definition: Messenger.h:149
std::string const ToString() const
Definition: Structs.cxx:114
ExposureMap fFakeExposure
POT in 1e12 of an exposure to fake up.
void fnex::EventListManipulator::ConcatenateLists ( fnex::EventListMap const &  lists,
std::string const &  dirName,
bool const &  byEpoch,
bool const &  byInteraction,
bool const &  byFile 
)
private

Definition at line 1720 of file EventListManipulator.cxx.

References fnex::MetaData::detector, fnex::MetaData::epoch, fnex::MetaData::EpochString(), fnex::MetaData::fileType, FillHistograms(), fnex::MetaData::interactionType, fnex::MetaData::isMC, fnex::kUnknownFileType, fnex::kUnknownInteraction, fetch_tb_beamline_files::md, fnex::MetaData::selectionType, and string.

Referenced by MakeHistograms().

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  }
enum fnex::file_type FileType_t
enum fnex::interaction_type InteractionType_t
std::map< fnex::MetaData, fnex::EventList > EventListMap
Definition: Event.h:186
std::string const EpochString() const
Definition: Structs.cxx:101
novadaq::cnv::DetId detector
Definition: Structs.h:50
fnex::FileType_t fileType
Definition: Structs.h:51
fnex::SelectionType_t selectionType
Definition: Structs.h:52
std::string dirName
Definition: PlotSpectra.h:47
fnex::InteractionType_t interactionType
Definition: Structs.h:53
void FillHistograms(fnex::EventListMap const &evLists, std::string const &dirName)
enum BeamMode string
std::unique_ptr< TH1F > fnex::EventListManipulator::Cosmics2017 ( fnex::MetaData const &  md)
private

Definition at line 2767 of file EventListManipulator.cxx.

References plotROC::fileName, cet::search_path::find_file(), fNuECosmicHistFile, fNuEDataHistFile, fNuMuCosmicHistFile, analysePickle::hist, MECModelEnuComparisons::i, fnex::MetaData::IsNuMuSelected(), fnex::kNuESelectionHighPID, fnex::kNuESelectionLowPID, fnex::kNuESelectionMidPID, fnex::kNuESelectionPeripheral, LOG_DEBUG, LOG_WARNING, fnex::MetaData::selectionType, string, and fnex::MetaData::ToString().

Referenced by FarDetCosmicBackgroundHistAndScale().

2767  {
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
2778  std::string fileName = ( md.IsNuMuSelected() ) ? fNuMuCosmicHistFile : fNuECosmicHistFile;
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
#define LOG_DEBUG(stream)
Definition: Messenger.h:149
fileName
Definition: plotROC.py:78
std::string fNuECosmicHistFile
file containing the cosmic background histograms
std::string fNuMuCosmicHistFile
file containing the cosmic background histograms
::xsd::cxx::tree::exception< char > exception
Definition: Database.h:225
Timing fit.
#define LOG_WARNING(category)
std::string fNuEDataHistFile
file containing the data live time for cosmic normalization
enum BeamMode string
std::unique_ptr< TH1F > fnex::EventListManipulator::Cosmics2018 ( fnex::MetaData const &  md,
TFile &  tf 
)
private

Get the nue 2018 CAFAna cosmic histogram.

Definition at line 2915 of file EventListManipulator.cxx.

References fnex::MetaData::BeamType(), analysePickle::hist, MECModelEnuComparisons::i, fnex::kFHC, fnex::kNuESelectionHighPID, fnex::kNuESelectionLowPID, fnex::kNuESelectionPeripheral, fnex::kRHC, LOG_DEBUG, LOG_VERBATIM, LOG_WARNING, fnex::MetaData::selectionType, and fnex::MetaData::ToString().

Referenced by FarDetCosmicBackgroundHistAndScale().

2916  {
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 
2931  if(md.selectionType == fnex::kNuESelectionLowPID ){
2932  firstbin = 0;
2933  lastbin = 10;
2934  }
2935  else if(md.selectionType == fnex::kNuESelectionHighPID){
2936  firstbin = 9;
2937  lastbin = 19;
2938  }
2939  else if(md.selectionType == fnex::kNuESelectionPeripheral){
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
#define LOG_DEBUG(stream)
Definition: Messenger.h:149
Timing fit.
#define LOG_WARNING(category)
#define LOG_VERBATIM(category)
void fnex::EventListManipulator::CreateFarDetCosmicBackgrounds ( fnex::EventListMap eventLists)
private

Definition at line 662 of file EventListManipulator.cxx.

References CreateFarDetCosmicBackgroundsFromHist(), fFDPOTSum, fSelectionsToUse, fnex::kCosmicBackgroundFile, fnex::kCosmicMuon, novadaq::cnv::kFARDET, fetch_tb_beamline_files::md, and ss.

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  }
Float_t ss
Definition: plot.C:24
Far Detector at Ash River, MN.
fnex::EventList CreateFarDetCosmicBackgroundsFromHist(fnex::MetaData const &md, fnex::SpillSummary const &ss)
std::map< std::string, fnex::SpillSummary > fFDPOTSum
data POT per epoch in the FD
std::map< novadaq::cnv::DetId, std::set< fnex::SelectionType_t > > fSelectionsToUse
which selection types to use?
fnex::EventList fnex::EventListManipulator::CreateFarDetCosmicBackgroundsFromEventList ( fnex::MetaData const &  md,
fnex::SpillSummary const &  ss,
fnex::EventListMap const &  listmap 
)
private

Definition at line 2459 of file EventListManipulator.cxx.

References fnex::MetaData::epoch, fnex::MetaData::EpochString(), evt, fnex::SpillSummary::goodPOT, fnex::kCosmicBackgroundFile, fnex::kCosmicMuon, fnex::kFake_Weight, novadaq::cnv::kFARDET, fnex::kRPACCQE_Weight, fnex::kTrueCCNC, fnex::kTrueE, fnex::kTrueIntType, fnex::kTruePDG, fnex::kTruePDGOrig, fnex::kXSecCVPPFX_Weight, LOG_DEBUG, fnex::MetaData::selectionType, and fnex::MetaData::ToString().

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  }
#define LOG_DEBUG(stream)
Definition: Messenger.h:149
Float_t ss
Definition: plot.C:24
static const unsigned char kXSecCVPPFX_Weight
Definition: VarVals.h:56
::xsd::cxx::tree::exception< char > exception
Definition: Database.h:225
static const unsigned char kFake_Weight
Definition: VarVals.h:39
static const unsigned char kTruePDGOrig
Definition: VarVals.h:46
static const unsigned char kTrueCCNC
Definition: VarVals.h:44
Far Detector at Ash River, MN.
int evt
static const unsigned char kTrueIntType
Definition: VarVals.h:45
static const unsigned char kRPACCQE_Weight
Definition: VarVals.h:57
static const unsigned char kTrueE
Definition: VarVals.h:42
static const unsigned char kTruePDG
Definition: VarVals.h:43
fnex::EventList fnex::EventListManipulator::CreateFarDetCosmicBackgroundsFromHist ( fnex::MetaData const &  md,
fnex::SpillSummary const &  ss 
)
private

Definition at line 2377 of file EventListManipulator.cxx.

References e, energy, fnex::MetaData::EpochString(), FarDetCosmicBackgroundHistAndScale(), fnex::SpillSummary::goodPOT, analysePickle::hist, makeTrainCVSamples::int, fnex::kCosmicBackgroundFile, fnex::kCosmicMuon, novadaq::cnv::kFARDET, fnex::kHad_RecoE, fnex::kLep_RecoE, fnex::kRPACCQE_Weight, fnex::kTrueCCNC, fnex::kTrueE, fnex::kTrueIntType, fnex::kTruePDG, fnex::kTruePDGOrig, fnex::kXSecCVPPFX_Weight, fnex::SpillSummary::liveTime, LOG_DEBUG, LOG_WARNING, fnex::SpillSummary::numGoodSpills, scale, fnex::MetaData::selectionType, fnex::MetaData::ToString(), fnex::SpillSummary::totalNumSpills, and fnex::SpillSummary::totalPOT.

Referenced by CreateFarDetCosmicBackgrounds().

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  }
#define LOG_DEBUG(stream)
Definition: Messenger.h:149
Float_t ss
Definition: plot.C:24
static const unsigned char kXSecCVPPFX_Weight
Definition: VarVals.h:56
static const unsigned char kTruePDGOrig
Definition: VarVals.h:46
static const unsigned char kTrueCCNC
Definition: VarVals.h:44
std::unique_ptr< TH1F > FarDetCosmicBackgroundHistAndScale(fnex::MetaData const &md, float &scale, fnex::SpillSummary const &ss)
Double_t scale
Definition: plot.C:25
Far Detector at Ash River, MN.
double energy
Definition: plottest35.C:25
static const unsigned char kTrueIntType
Definition: VarVals.h:45
#define LOG_WARNING(category)
static const unsigned char kRPACCQE_Weight
Definition: VarVals.h:57
std::shared_ptr< Event > EventPtr
Definition: Event.h:50
static const unsigned char kTrueE
Definition: VarVals.h:42
static const unsigned char kHad_RecoE
Definition: VarVals.h:34
Float_t e
Definition: plot.C:35
static const unsigned char kLep_RecoE
Definition: VarVals.h:35
static const unsigned char kTruePDG
Definition: VarVals.h:43
void fnex::EventListManipulator::CreateNDPeripheralLists ( fnex::EventListMap listmap)
private

Definition at line 2850 of file EventListManipulator.cxx.

References fnex::EventList::AddEvent(), fnex::MetaData::detector, fnex::MetaData::EpochString(), fnex::MetaData::fileType, fSelectionsToUse, fnex::MetaData::interactionType, fnex::MetaData::isMC, it, novadaq::cnv::kNEARDET, fnex::kNuESelectionHighPID, fnex::kNuESelectionPeripheral, LOG_DEBUG, fetch_tb_beamline_files::md, fnex::MetaData::selectionType, and fnex::MetaData::ToString().

Referenced by Deserialize().

2850  {
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
#define LOG_DEBUG(stream)
Definition: Messenger.h:149
set< int >::iterator it
std::string const EpochString() const
Definition: Structs.cxx:101
void AddEvent(EventPtr ev)
Definition: Event.h:115
std::string const ToString() const
Definition: Structs.cxx:114
Near Detector in the NuMI cavern.
novadaq::cnv::DetId detector
Definition: Structs.h:50
fnex::FileType_t fileType
Definition: Structs.h:51
fnex::SelectionType_t selectionType
Definition: Structs.h:52
fnex::InteractionType_t interactionType
Definition: Structs.h:53
std::map< novadaq::cnv::DetId, std::set< fnex::SelectionType_t > > fSelectionsToUse
which selection types to use?
fnex::EventList fnex::EventListManipulator::CreatePoissonFakeData ( fnex::MetaData const &  md,
fnex::SpillSummary const &  ss,
fnex::EventList const &  origList 
)
private

Definition at line 2530 of file EventListManipulator.cxx.

References fnex::EventList::AddEvent(), fnex::VarBinning::BinInfo(), e, fRandom, analysePickle::hist, fnex::VarBinning::Instance(), fnex::kFake_Weight, fnex::kHad_RecoE, fnex::kLep_RecoE, fnex::EventList::ListMetaData(), LOG_DEBUG, runNovaSAM::outList, fnex::EventList::size(), and fnex::MetaData::ToString().

Referenced by FillFakeDataLists().

2533  {
2534  // make the EventList to return
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  }
#define LOG_DEBUG(stream)
Definition: Messenger.h:149
Float_t ss
Definition: plot.C:24
static const unsigned char kFake_Weight
Definition: VarVals.h:39
TRandom3 fRandom
Random number generator to use creating fake data lists.
list outList
Definition: runNovaSAM.py:344
static VarBinning * Instance()
Definition: VarBinning.cxx:15
std::shared_ptr< Event > EventPtr
Definition: Event.h:50
static const unsigned char kHad_RecoE
Definition: VarVals.h:34
Float_t e
Definition: plot.C:35
static const unsigned char kLep_RecoE
Definition: VarVals.h:35
fnex::Binning const * BinInfo(std::string const &varName)
Definition: VarBinning.cxx:65
fnex::EventListMap fnex::EventListManipulator::Deserialize ( fnex::DataMC_t  dataMC = fnex::kBoth)

Definition at line 692 of file EventListManipulator.cxx.

References fnex::ShifterAndWeighter::CalcNearTrueEnergyRatioWeights(), CheckDataMCListConsistency(), CreateNDPeripheralLists(), dirName, ExtractEventsFromFile(), ExtractMetaDataFromFile(), MakeMiniprodValidationCuts::f, fCalibHistFile, fFakeExposure, fFNEXEventLabels, plotROC::fileName, FillFakeDataLists(), fInitialGuess, fMakeNDPeripheral, fTreeDirectories, fUseExtrapolation, fnex::ShifterAndWeighter::InitShiftsAndWeightsToUse(), fnex::ShifterAndWeighter::Instance(), fnex::kData, fnex::kMC, LOG_DEBUG, LOG_VERBATIM, and string.

Referenced by DeserializeData().

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  }
#define LOG_DEBUG(stream)
Definition: Messenger.h:149
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
fileName
Definition: plotROC.py:78
void CreateNDPeripheralLists(fnex::EventListMap &listmap)
void InitShiftsAndWeightsToUse(fnex::InputPoint const &curPoint, std::string CalibFilename="FNEX/data/calibration/CalibrationSystematics2018.root")
::xsd::cxx::tree::exception< char > exception
Definition: Database.h:225
std::map< fnex::MetaData, fnex::EventList > EventListMap
Definition: Event.h:186
static ShifterAndWeighter * Instance()
fnex::InputPoint fInitialGuess
the fake data input point
void ExtractEventsFromFile(TFile *eventFile, std::string const &dirName, fnex::DataMC_t const &dataMC, std::pair< fnex::MetaData, fnex::SpillSummary > const &mdToSS, fnex::EventListMap &eventLists)
ExposureMap fFakeExposure
POT in 1e12 of an exposure to fake up.
std::string dirName
Definition: PlotSpectra.h:47
void FillFakeDataLists(fnex::EventListMap &lists)
void CalcNearTrueEnergyRatioWeights(const fnex::EventListMap *pEventListMap)
std::vector< std::string > fFNEXEventLabels
Labels in input files holding FNEX Events.
std::vector< std::string > fTreeDirectories
directory holding the input trees
bool fUseExtrapolation
are we extrapolating or using covariance matrices
void ExtractMetaDataFromFile(TFile *metadataFile, std::string const &dirName, fnex::DataMC_t const &dataMC, std::map< fnex::MetaData, fnex::SpillSummary > &mdToSS)
#define LOG_VERBATIM(category)
enum BeamMode string
fnex::EventListMap fnex::EventListManipulator::DeserializeData ( )

Definition at line 814 of file EventListManipulator.cxx.

References Deserialize(), fFakeExposure, fnex::kBoth, and fnex::kData.

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  }
fnex::EventListMap Deserialize(fnex::DataMC_t dataMC=fnex::kBoth)
ExposureMap fFakeExposure
POT in 1e12 of an exposure to fake up.
fnex::EventListMap fnex::EventListManipulator::DeserializeMC ( )
inline

Definition at line 60 of file EventListManipulator.h.

References fillBadChanDBTables::det, dirName, fnex::kMC, parse_dependency_file_t::list, fetch_tb_beamline_files::md, scale, ss, and string.

60 { return this->Deserialize(fnex::kMC); }
fnex::EventListMap Deserialize(fnex::DataMC_t dataMC=fnex::kBoth)
void fnex::EventListManipulator::DetermineSelections ( std::vector< fhicl::ParameterSet > const &  selections)
private

Definition at line 269 of file EventListManipulator.cxx.

References fillBadChanDBTables::det, fSelectionsToUse, novadaq::cnv::kFARDET, novadaq::cnv::kNEARDET, fnex::kNuESelectionHighPID, fnex::kNuESelectionLowPID, fnex::kNuESelectionPeripheral, fnex::kNuMuSelectionQ1, fnex::kNuMuSelectionQ2, fnex::kNuMuSelectionQ3, fnex::kNuMuSelectionQ4, and string.

Referenced by reconfigure().

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  }
Far Detector at Ash River, MN.
Near Detector in the NuMI cavern.
std::map< novadaq::cnv::DetId, std::set< fnex::SelectionType_t > > fSelectionsToUse
which selection types to use?
enum BeamMode string
float fnex::EventListManipulator::EpochEventCap ( fnex::MetaData const &  md,
int  treeEvents 
)
private

Definition at line 391 of file EventListManipulator.cxx.

References fEventCaps, fMaxEventsPerTree, findDuplicateFiles::key, and fnex::MetaData::Key().

Referenced by FillEventList().

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  }
float fMaxEventsPerTree
maximum number of events to use per tree
std::map< long, float > fEventCaps
void fnex::EventListManipulator::ExtractEventsFromFile ( TFile *  eventFile,
std::string const &  dirName,
fnex::DataMC_t const &  dataMC,
std::pair< fnex::MetaData, fnex::SpillSummary > const &  mdToSS,
fnex::EventListMap eventLists 
)
private

Definition at line 583 of file EventListManipulator.cxx.

References art::rootNames::eventTreeName(), ss, and string.

Referenced by Deserialize().

588  {
589  // find the event tree using the metadata information
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  }
Float_t ss
Definition: plot.C:24
std::string const & eventTreeName()
::xsd::cxx::tree::exception< char > exception
Definition: Database.h:225
std::string dirName
Definition: PlotSpectra.h:47
enum BeamMode string
void fnex::EventListManipulator::ExtractMetaDataFromFile ( TFile *  metadataFile,
std::string const &  dirName,
fnex::DataMC_t const &  dataMC,
std::map< fnex::MetaData, fnex::SpillSummary > &  mdToSS 
)
private

Definition at line 421 of file EventListManipulator.cxx.

References fnex::MetaData::BeamType(), fnex::MetaData::detector, fnex::MetaData::fileType, fSumSubEpochs, fTotalFDDataPOTFHC, fTotalFDDataPOTRHC, fTotalNDDataPOTFHC, fTotalNDDataPOTRHC, fnex::SpillSummary::goodPOT, fnex::MetaData::isMC, fnex::kCosmicBackgroundFile, fnex::kDataFile, novadaq::cnv::kFARDET, fnex::kFHC, novadaq::cnv::kNEARDET, fnex::kNuMuSelectionQ1, fnex::kRHC, fnex::SpillSummary::liveTime, LOG_DEBUG, fetch_tb_beamline_files::md, fnex::MetaData::Period(), fnex::MetaData::selectionType, ss, fnex::SpillSummary::totalPOT, and UseEventsFromMetaData().

Referenced by Deserialize().

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  }
#define LOG_DEBUG(stream)
Definition: Messenger.h:149
double fTotalNDDataPOTFHC
total ND data POT FHC
Float_t ss
Definition: plot.C:24
bool fSumSubEpochs
do we want to sum all sub epochs into a single one
::xsd::cxx::tree::exception< char > exception
Definition: Database.h:225
long const Period() const
Definition: Structs.h:62
double fTotalFDDataPOTRHC
total FD data POT RHC
Far Detector at Ash River, MN.
double fTotalFDDataPOTFHC
total FD data POT FHC
Near Detector in the NuMI cavern.
novadaq::cnv::DetId detector
Definition: Structs.h:50
fnex::FileType_t fileType
Definition: Structs.h:51
fnex::SelectionType_t selectionType
Definition: Structs.h:52
double fTotalNDDataPOTRHC
total ND data POT RHC
std::string dirName
Definition: PlotSpectra.h:47
bool UseEventsFromMetaData(fnex::MetaData const &md, fnex::DataMC_t const dataMC)
fnex::BeamType_t const BeamType() const
Definition: Structs.cxx:165
std::unique_ptr< TH1F > fnex::EventListManipulator::FarDetCosmicBackgroundHistAndScale ( fnex::MetaData const &  md,
float &  scale,
fnex::SpillSummary const &  ss 
)
private

Definition at line 2032 of file EventListManipulator.cxx.

References fnex::MetaData::BeamType(), Cosmics2017(), Cosmics2018(), fnex::MetaData::epoch, fAnalysisYear, plotROC::fileName, cet::search_path::find_file(), fNuECosmicHistFile, fNuMuCosmicHistFile, fTotalFDDataPOTFHC, fTotalFDDataPOTRHC, fnex::SpillSummary::goodPOT, analysePickle::hist, compareCafs::histName, MECModelEnuComparisons::i, fnex::MetaData::IsNuMuSelected(), fnex::kFHC, fnex::kNuESelectionHighPID, fnex::kNuESelectionLowPID, fnex::kNuESelectionMidPID, fnex::kNuMuSelection, fnex::kNuMuSelectionQ1, fnex::kNuMuSelectionQ2, fnex::kNuMuSelectionQ3, fnex::kNuMuSelectionQ4, fnex::kRHC, LOG_DEBUG, LOG_WARNING, fnex::MetaData::selectionType, string, and fnex::MetaData::ToString().

Referenced by CreateFarDetCosmicBackgroundsFromHist().

2035  {
2036  std::unique_ptr<TH1F> hist;
2038  auto epoch = md.epoch;
2039 
2040  // get the file containing the histograms
2041  std::string fileName = ( md.IsNuMuSelected() ) ? fNuMuCosmicHistFile : fNuECosmicHistFile;
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();
2074  if(md.selectionType == fnex::kNuMuSelection){
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 &&
2301  md.selectionType < fnex::kNuMuSelectionQ1){
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  }
#define LOG_DEBUG(stream)
Definition: Messenger.h:149
fileName
Definition: plotROC.py:78
std::string fNuECosmicHistFile
file containing the cosmic background histograms
std::string fNuMuCosmicHistFile
file containing the cosmic background histograms
std::string fAnalysisYear
analysis year, "2016" for SA and "2017" for TA
Float_t ss
Definition: plot.C:24
::xsd::cxx::tree::exception< char > exception
Definition: Database.h:225
std::unique_ptr< TH1F > Cosmics2017(fnex::MetaData const &md)
Timing fit.
double fTotalFDDataPOTRHC
total FD data POT RHC
Double_t scale
Definition: plot.C:25
double fTotalFDDataPOTFHC
total FD data POT FHC
std::unique_ptr< TH1F > Cosmics2018(fnex::MetaData const &md, TFile &tf)
Get the nue 2018 CAFAna cosmic histogram.
#define LOG_WARNING(category)
enum BeamMode string
fnex::EventList fnex::EventListManipulator::FillEventList ( TTree *  eventTree,
fnex::MetaData const &  md,
fnex::SpillSummary ss 
)
private

Definition at line 889 of file EventListManipulator.cxx.

References fnex::EventList::AddEvent(), fnex::MetaData::detector, EpochEventCap(), fnex::DataVars::fFake_Weight, fMaxNuEnergy, fUseEventId, fnex::SpillSummary::goodPOT, fnex::MetaData::isMC, fnex::MetaData::IsNuESelected(), fnex::kCCQEPauliSupViaKF, findDuplicateFiles::key, fnex::MetaData::Key(), fnex::KeyToVarName(), novadaq::cnv::kFARDET, fnex::kMaCCQE, fnex::kNu_RecoE, fnex::kNuESelection, fnex::kNuMuSelection, fnex::EventList::ListMetaData(), fnex::EventList::ListSpillSummary(), LOG_DEBUG, LOG_VERBATIM, LOG_WARNING, central_limit::rand, string, fnex::MetaData::ToString(), and fnex::SpillSummary::totalPOT.

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
959  fnex::SelectionType_t selection = (md.IsNuESelected()) ? fnex::kNuESelection : fnex::kNuMuSelection;
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  }
#define LOG_DEBUG(stream)
Definition: Messenger.h:149
static const unsigned char kNu_RecoE
Definition: VarVals.h:36
Base container for the MC related Vars that constitute an event.
Definition: VarVals.h:930
::xsd::cxx::tree::exception< char > exception
Definition: Database.h:225
std::map< fnex::SelectionType_t, double > fMaxNuEnergy
maximum neutrino energy to go into the lists
Far Detector at Ash River, MN.
enum fnex::sel_type SelectionType_t
static const unsigned char kMaCCQE
Definition: VarVals.h:70
float fFake_Weight
Weight for fake data events.
Definition: VarVals.h:735
float EpochEventCap(fnex::MetaData const &md, int treeEvents)
static const unsigned char kCCQEPauliSupViaKF
Definition: VarVals.h:81
bool fUseEventId
use if you need access to run/subrun/event/etc
#define LOG_WARNING(category)
std::shared_ptr< Event > EventPtr
Definition: Event.h:50
std::string KeyToVarName(unsigned char const &key)
Definition: VarVals.h:445
#define LOG_VERBATIM(category)
Base container for the MC related Vars that constitute an event.
Definition: VarVals.h:874
enum BeamMode string
void fnex::EventListManipulator::FillFakeDataLists ( fnex::EventListMap lists)
private

Definition at line 1780 of file EventListManipulator.cxx.

References CheckMCListsForFakeDataEpochs(), CreatePoissonFakeData(), fCalibHistFile, fFakeDataPoint, fFakeExposure, FillFakeDataListsForDetector(), FillHistogram(), fPoissonFakeData, fTFS, fUseExtrapolation, drawEvents::histDir, fnex::ShifterAndWeighter::Instance(), novadaq::cnv::kFARDET, novadaq::cnv::kNEARDET, LOG_DEBUG, fetch_tb_beamline_files::md, std::min(), and art::TFileDirectory::mkdir().

Referenced by Deserialize().

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  }
#define LOG_DEBUG(stream)
Definition: Messenger.h:149
std::string fCalibHistFile
file containing calibration systematic histograms
TFileDirectory mkdir(std::string const &dir, std::string const &descr="")
static ShifterAndWeighter * Instance()
fnex::EventList CreatePoissonFakeData(fnex::MetaData const &md, fnex::SpillSummary const &ss, fnex::EventList const &origList)
Far Detector at Ash River, MN.
art::ServiceHandle< art::TFileService > fTFS
TFileService.
bool fPoissonFakeData
make the fake data Poisson distributed
fnex::InputPoint fFakeDataPoint
the fake data input point
void CheckMCListsForFakeDataEpochs(fnex::EventListMap &lists, novadaq::cnv::DetId const &det)
Near Detector in the NuMI cavern.
ExposureMap fFakeExposure
POT in 1e12 of an exposure to fake up.
void FillHistogram(art::TFileDirectory &tfdir, fnex::EventList const &evList, std::string const &listType, double const &normalization)
bool fUseExtrapolation
are we extrapolating or using covariance matrices
T min(const caf::Proxy< T > &a, T b)
void FillFakeDataListsForDetector(fnex::EventListMap &lists, novadaq::cnv::DetId const &det, std::map< int, double > const &detEpochToMinPOT)
void fnex::EventListManipulator::FillFakeDataListsForDetector ( fnex::EventListMap lists,
novadaq::cnv::DetId const &  det,
std::map< int, double > const &  detEpochToMinPOT 
)
private

Definition at line 1906 of file EventListManipulator.cxx.

References fnex::EventList::AddEvent(), AddFakeDataToEventList(), fnex::EventList::AddGoodPOT(), fnex::EventList::AddGoodSpills(), fFakeExposure, fnex::SpillSummary::goodPOT, fnex::kBoth, findDuplicateFiles::key, fnex::kFake_Weight, LOG_DEBUG, fetch_tb_beamline_files::md, std::min(), ss, and UseEventsFromMetaData().

Referenced by FillFakeDataLists().

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  }
#define LOG_DEBUG(stream)
Definition: Messenger.h:149
Float_t ss
Definition: plot.C:24
std::map< fnex::MetaData, fnex::EventList > EventListMap
Definition: Event.h:186
static const unsigned char kFake_Weight
Definition: VarVals.h:39
fnex::EventList AddFakeDataToEventList(fnex::MetaData const &md, fnex::SpillSummary const &ss, fnex::EventList const &mcList)
ExposureMap fFakeExposure
POT in 1e12 of an exposure to fake up.
bool UseEventsFromMetaData(fnex::MetaData const &md, fnex::DataMC_t const dataMC)
T min(const caf::Proxy< T > &a, T b)
void fnex::EventListManipulator::FillHistogram ( art::TFileDirectory tfdir,
fnex::EventList const &  evList,
std::string const &  listType,
double const &  normalization 
)
private

Definition at line 1285 of file EventListManipulator.cxx.

References genie::units::A, fNeutrinoEnergy, fUseExtrapolation, compareCafs::histName, MECModelEnuComparisons::i, fnex::ShifterAndWeighter::Instance(), fnex::kAllWgt, simb::kCC, fnex::kCosmicBackgroundFile, fnex::KeyToVarName(), novadaq::cnv::kNEARDET, fnex::kNu_RecoE, fnex::kNuMuSelection, fnex::kSystAndOscWgt, fnex::kTrueCCNC, fnex::kTrueE, fnex::kTrueParentDecay, fnex::kTrueParentPT, fnex::kTrueParentPZ, fnex::kTrueParentTargetPDG, fnex::kTruePDG, fnex::kTruePDGOrig, fnex::EventList::ListMetaData(), LOG_DEBUG, art::TFileDirectory::make(), fetch_tb_beamline_files::md, string, febshutoff_auto::val, ana::weight, and fnex::ShifterAndWeighter::Weight().

Referenced by FillFakeDataLists(), and FillHistograms().

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  }
#define LOG_DEBUG(stream)
Definition: Messenger.h:149
static const unsigned char kNu_RecoE
Definition: VarVals.h:36
enum fnex::weight_type WeightType_t
EnergyHistMap fNeutrinoEnergy
neutrino energy
static const unsigned char kTrueParentPT
Definition: VarVals.h:48
const Var weight
::xsd::cxx::tree::exception< char > exception
Definition: Database.h:225
static ShifterAndWeighter * Instance()
static const unsigned char kTruePDGOrig
Definition: VarVals.h:46
static const unsigned char kTrueCCNC
Definition: VarVals.h:44
double Weight(fnex::EventPtr const &curEvent, fnex::MetaData const &md, fnex::WeightType_t const wgtType=kAllWgt)
Near Detector in the NuMI cavern.
static const unsigned char kTrueParentDecay
Definition: VarVals.h:50
static const unsigned char kTrueParentTargetPDG
Definition: VarVals.h:51
static const double A
Definition: Units.h:82
T * make(ARGS...args) const
static const unsigned char kTrueE
Definition: VarVals.h:42
bool fUseExtrapolation
are we extrapolating or using covariance matrices
std::string KeyToVarName(unsigned char const &key)
Definition: VarVals.h:445
static const unsigned char kTrueParentPZ
Definition: VarVals.h:49
static const unsigned char kTruePDG
Definition: VarVals.h:43
enum BeamMode string
void fnex::EventListManipulator::FillHistograms ( fnex::EventListMap const &  evLists,
std::string const &  dirName 
)
private

Definition at line 1601 of file EventListManipulator.cxx.

References fFakeExposure, fnex::MetaData::fileType, FillHistogram(), fTFS, drawEvents::histDir, fnex::MetaData::interactionType, fnex::MetaData::isMC, it, fnex::kDataFile, novadaq::cnv::kFARDET, novadaq::cnv::kNEARDET, fnex::kNuESelection, fnex::kNuESelectionHighPID, fnex::kNuESelectionLowPID, fnex::kNuESelectionMidPID, fnex::kNuESelectionPeripheral, fnex::kNuMuSelection, fnex::kNuMuSelectionQ1, fnex::kNuMuSelectionQ2, fnex::kNuMuSelectionQ3, fnex::kNuMuSelectionQ4, fnex::kUnknownInteraction, LOG_DEBUG, fetch_tb_beamline_files::md, art::TFileDirectory::mkdir(), string, and fnex::MetaData::ToString().

Referenced by ConcatenateLists(), and MakeHistograms().

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  }
#define LOG_DEBUG(stream)
Definition: Messenger.h:149
set< int >::iterator it
TFileDirectory mkdir(std::string const &dir, std::string const &descr="")
std::string const ToString() const
Definition: Structs.cxx:114
Far Detector at Ash River, MN.
art::ServiceHandle< art::TFileService > fTFS
TFileService.
Near Detector in the NuMI cavern.
fnex::FileType_t fileType
Definition: Structs.h:51
ExposureMap fFakeExposure
POT in 1e12 of an exposure to fake up.
void FillHistogram(art::TFileDirectory &tfdir, fnex::EventList const &evList, std::string const &listType, double const &normalization)
std::string dirName
Definition: PlotSpectra.h:47
std::unordered_map< fnex::MetaData, double, fnex::MetaDataHasher > ExposureMap
fnex::InteractionType_t interactionType
Definition: Structs.h:53
enum BeamMode string
void fnex::EventListManipulator::FixNearDetMCExposures ( fnex::EventListMap lists)
private

Definition at line 2596 of file EventListManipulator.cxx.

References fnex::EventList::AddEvent(), fnex::MetaData::detector, dirName, art::rootNames::eventTreeName(), MakeMiniprodValidationCuts::f, fFNEXEventLabels, plotROC::fileName, fTreeDirectories, fnex::SpillSummary::goodPOT, fnex::MetaData::isMC, novadaq::cnv::kNEARDET, LOG_DEBUG, LOG_VERBATIM, fnex::MetaData::Period(), string, std::swap(), fnex::MetaData::ToString(), and fnex::SpillSummary::totalPOT.

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
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  }
#define LOG_DEBUG(stream)
Definition: Messenger.h:149
fileName
Definition: plotROC.py:78
std::string const & eventTreeName()
std::map< fnex::MetaData, fnex::EventList > EventListMap
Definition: Event.h:186
long const Period() const
Definition: Structs.h:62
void swap(art::HLTGlobalStatus &lhs, art::HLTGlobalStatus &rhs)
std::string const ToString() const
Definition: Structs.cxx:114
Near Detector in the NuMI cavern.
novadaq::cnv::DetId detector
Definition: Structs.h:50
std::string dirName
Definition: PlotSpectra.h:47
std::vector< std::string > fFNEXEventLabels
Labels in input files holding FNEX Events.
std::vector< std::string > fTreeDirectories
directory holding the input trees
#define LOG_VERBATIM(category)
enum BeamMode string
void fnex::EventListManipulator::LoadVarVals ( fnex::Event fnexEvent,
bool  isMC,
std::vector< std::string > const &  indepVarNames,
std::vector< float > const &  indepVarVals 
)
private

Definition at line 823 of file EventListManipulator.cxx.

References fnex::Event::dataVals, fnex::IsDataVar(), fnex::IsMCVar(), LOG_DEBUG, LOG_WARNING, fnex::Event::mcVals, elec2geo::pos, and string.

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  }
#define LOG_DEBUG(stream)
Definition: Messenger.h:149
::xsd::cxx::tree::exception< char > exception
Definition: Database.h:225
#define LOG_WARNING(category)
bool IsMCVar(unsigned char const &key)
Definition: VarVals.h:604
bool IsDataVar(unsigned char const &key)
Definition: VarVals.h:617
std::unique_ptr< MCVarVals > mcVals
Definition: Event.h:36
std::unique_ptr< DataVarVals > dataVals
Definition: Event.h:35
enum BeamMode string
void fnex::EventListManipulator::MakeHistograms ( fnex::EventListMap const &  lists,
std::string const &  dirName 
)
private

Definition at line 1270 of file EventListManipulator.cxx.

References ConcatenateLists(), and FillHistograms().

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  }
void ConcatenateLists(fnex::EventListMap const &lists, std::string const &dirName, bool const &byEpoch, bool const &byInteraction, bool const &byFile)
std::string dirName
Definition: PlotSpectra.h:47
void FillHistograms(fnex::EventListMap const &evLists, std::string const &dirName)
void fnex::EventListManipulator::reconfigure ( const fhicl::ParameterSet p)

Definition at line 49 of file EventListManipulator.cxx.

References append(), fnex::cInteractionType_Strings, fnex::cSelectionType_Strings, fillBadChanDBTables::det, DetermineSelections(), fAnalysisYear, fCalibHistFile, fCorrelations, fEventCaps, fFakeDataPoint, fFakeExposure, fFDPOTSum, fFNEXEventLabels, fLoadAllEventLists, fMakeNDPeripheral, fMaxEventsPerTree, fMaxNuEnergy, fNuECosmicHistFile, fNuEDataHistFile, fNuMuCosmicHistFile, fPeriodsToUse, fPoissonFakeData, fRandom, fSelectionsToUse, fSumSubEpochs, fTotalFDDataPOTFHC, fTotalFDDataPOTRHC, fTreeDirectories, fUseCosmicHists, fUseEventId, fUseExtrapolation, fhicl::ParameterSet::get(), MECModelEnuComparisons::i, ip, fnex::kCosmicMuon, fnex::kDataFile, novadaq::cnv::kFARDET, novadaq::cnv::kNEARDET, fnex::kNuESelection, fnex::kNuMuSelection, fnex::kUnknownInteraction, fnex::kUnknownSelection, LOG_DEBUG, LOG_VERBATIM, fetch_tb_beamline_files::md, DCS_db_parser::period, string, fnex::StringToKey(), std::swap(), and fnex::MetaData::ToString().

Referenced by EventListManipulator().

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  }
#define LOG_DEBUG(stream)
Definition: Messenger.h:149
float fMaxEventsPerTree
maximum number of events to use per tree
std::string fCalibHistFile
file containing calibration systematic histograms
std::string fNuECosmicHistFile
file containing the cosmic background histograms
std::string fNuMuCosmicHistFile
file containing the cosmic background histograms
std::string fAnalysisYear
analysis year, "2016" for SA and "2017" for TA
bool fSumSubEpochs
do we want to sum all sub epochs into a single one
TString ip
Definition: loadincs.C:5
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
void append()
Definition: append.C:24
double fTotalFDDataPOTRHC
total FD data POT RHC
void swap(art::HLTGlobalStatus &lhs, art::HLTGlobalStatus &rhs)
const XML_Char * s
Definition: expat.h:262
Far Detector at Ash River, MN.
double fTotalFDDataPOTFHC
total FD data POT FHC
std::map< long, float > fEventCaps
bool fPoissonFakeData
make the fake data Poisson distributed
fnex::InputPoint fFakeDataPoint
the fake data input point
Near Detector in the NuMI cavern.
std::set< int > fPeriodsToUse
which periods to use in the analysis
bool fUseEventId
use if you need access to run/subrun/event/etc
std::string fNuEDataHistFile
file containing the data live time for cosmic normalization
ExposureMap fFakeExposure
POT in 1e12 of an exposure to fake up.
void DetermineSelections(std::vector< fhicl::ParameterSet > const &selections)
std::unique_ptr< Correlations > fCorrelations
const std::string cInteractionType_Strings[9]
Definition: Constants.h:150
std::map< std::string, fnex::SpillSummary > fFDPOTSum
data POT per epoch in the FD
std::vector< std::string > fFNEXEventLabels
Labels in input files holding FNEX Events.
bool fLoadAllEventLists
force all event lists to be loaded
std::vector< std::string > fTreeDirectories
directory holding the input trees
bool fUseExtrapolation
are we extrapolating or using covariance matrices
std::map< novadaq::cnv::DetId, std::set< fnex::SelectionType_t > > fSelectionsToUse
which selection types to use?
#define LOG_VERBATIM(category)
const std::string cSelectionType_Strings[11]
Definition: Constants.h:101
long StringToKey(std::string const &str)
Definition: Structs.cxx:15
enum BeamMode string
double fnex::EventListManipulator::RescalePOTForEventCap ( size_t &  totalEvents,
size_t const &  eventCap 
)
private

Definition at line 1126 of file EventListManipulator.cxx.

Referenced by RescalePOTForEventCap().

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  }
double fnex::EventListManipulator::RescalePOTForEventCap ( double &  totalEvents,
size_t const &  eventCap 
)
private

Definition at line 1111 of file EventListManipulator.cxx.

References RescalePOTForEventCap(), and produceInfoGainMatrix::total.

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  }
double RescalePOTForEventCap(size_t &totalEvents, size_t const &eventCap)
bool fnex::EventListManipulator::UseEventsFromMetaData ( fnex::MetaData const &  md,
fnex::DataMC_t const  dataMC 
)
private

Definition at line 319 of file EventListManipulator.cxx.

References fnex::cSelectionType_Strings, fnex::MetaData::detector, fnex::MetaData::DetectorString(), e, fLoadAllEventLists, fPeriodsToUse, fSelectionsToUse, fnex::MetaData::isMC, fnex::kData, fnex::kFakeData, fnex::kMC, LOG_DEBUG, LOG_VERBATIM, fnex::MetaData::Period(), fnex::MetaData::selectionType, and fnex::MetaData::ToString().

Referenced by ExtractMetaDataFromFile(), and FillFakeDataListsForDetector().

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  << " "
373  << fnex::cSelectionType_Strings[md.selectionType];
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  }
#define LOG_DEBUG(stream)
Definition: Messenger.h:149
std::set< int > fPeriodsToUse
which periods to use in the analysis
bool fLoadAllEventLists
force all event lists to be loaded
std::map< novadaq::cnv::DetId, std::set< fnex::SelectionType_t > > fSelectionsToUse
which selection types to use?
#define LOG_VERBATIM(category)
const std::string cSelectionType_Strings[11]
Definition: Constants.h:101
Float_t e
Definition: plot.C:35
void fnex::EventListManipulator::WriteCappedListsAsTrees ( fnex::EventListMap const &  list)
private

Definition at line 2719 of file EventListManipulator.cxx.

References dir, file, fetch_tb_beamline_files::md, pot, fnex::SerializeEvents(), ss, string, and fnex::MetaData::ToString().

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  }
void SerializeEvents(TTree *eventTree, fnex::MetaData const &md, std::vector< fnex::EventInfo > const &events)
Definition: Event.cxx:82
Float_t ss
Definition: plot.C:24
std::string const ToString() const
Definition: Structs.cxx:114
#define pot
TDirectory * dir
Definition: macro.C:5
TFile * file
Definition: cellShifts.C:17
enum BeamMode string

Member Data Documentation

std::string fnex::EventListManipulator::fAnalysisYear
private

analysis year, "2016" for SA and "2017" for TA

Definition at line 176 of file EventListManipulator.h.

Referenced by FarDetCosmicBackgroundHistAndScale(), and reconfigure().

std::string fnex::EventListManipulator::fCalibHistFile
private

file containing calibration systematic histograms

Definition at line 167 of file EventListManipulator.h.

Referenced by Deserialize(), FillFakeDataLists(), and reconfigure().

std::unique_ptr<Correlations> fnex::EventListManipulator::fCorrelations
private

For correlating shifts in fake data. Uses a single set of correlations here, will have to be a bit fancier for the joint fit

Definition at line 159 of file EventListManipulator.h.

Referenced by reconfigure().

std::map<long, float> fnex::EventListManipulator::fEventCaps
private

maximum fraction of events to use from each combination of detector, filetype and epoch

Definition at line 152 of file EventListManipulator.h.

Referenced by EpochEventCap(), EventListManipulator(), and reconfigure().

fnex::InputPoint fnex::EventListManipulator::fFakeDataPoint
private

the fake data input point

Definition at line 150 of file EventListManipulator.h.

Referenced by FillFakeDataLists(), and reconfigure().

ExposureMap fnex::EventListManipulator::fFakeExposure
private
std::map<std::string, fnex::SpillSummary> fnex::EventListManipulator::fFDPOTSum
private

data POT per epoch in the FD

Definition at line 172 of file EventListManipulator.h.

Referenced by CheckDataMCListConsistency(), CreateFarDetCosmicBackgrounds(), and reconfigure().

std::vector<std::string> fnex::EventListManipulator::fFNEXEventLabels
private

Labels in input files holding FNEX Events.

Definition at line 145 of file EventListManipulator.h.

Referenced by Deserialize(), FixNearDetMCExposures(), and reconfigure().

fnex::InputPoint fnex::EventListManipulator::fInitialGuess
private

the fake data input point

Definition at line 151 of file EventListManipulator.h.

Referenced by Deserialize().

bool fnex::EventListManipulator::fLoadAllEventLists
private

force all event lists to be loaded

Definition at line 179 of file EventListManipulator.h.

Referenced by reconfigure(), and UseEventsFromMetaData().

bool fnex::EventListManipulator::fMakeNDPeripheral
private

configuration to decide if we should make the ND peripheral eventlist or not

Definition at line 177 of file EventListManipulator.h.

Referenced by Deserialize(), and reconfigure().

float fnex::EventListManipulator::fMaxEventsPerTree
private

maximum number of events to use per tree

Definition at line 180 of file EventListManipulator.h.

Referenced by EpochEventCap(), and reconfigure().

std::map<fnex::SelectionType_t, double> fnex::EventListManipulator::fMaxNuEnergy
private

maximum neutrino energy to go into the lists

Definition at line 175 of file EventListManipulator.h.

Referenced by FillEventList(), and reconfigure().

EnergyHistMap fnex::EventListManipulator::fNeutrinoEnergy
private

neutrino energy

Definition at line 148 of file EventListManipulator.h.

Referenced by FillHistogram().

std::string fnex::EventListManipulator::fNuECosmicHistFile
private

file containing the cosmic background histograms

Definition at line 163 of file EventListManipulator.h.

Referenced by Cosmics2017(), FarDetCosmicBackgroundHistAndScale(), and reconfigure().

std::string fnex::EventListManipulator::fNuEDataHistFile
private

file containing the data live time for cosmic normalization

Definition at line 164 of file EventListManipulator.h.

Referenced by Cosmics2017(), and reconfigure().

std::string fnex::EventListManipulator::fNuMuCosmicHistFile
private

file containing the cosmic background histograms

Definition at line 162 of file EventListManipulator.h.

Referenced by Cosmics2017(), FarDetCosmicBackgroundHistAndScale(), and reconfigure().

std::set<int> fnex::EventListManipulator::fPeriodsToUse
private

which periods to use in the analysis

Definition at line 173 of file EventListManipulator.h.

Referenced by reconfigure(), and UseEventsFromMetaData().

bool fnex::EventListManipulator::fPoissonFakeData
private

make the fake data Poisson distributed

Definition at line 154 of file EventListManipulator.h.

Referenced by FillFakeDataLists(), and reconfigure().

TRandom3 fnex::EventListManipulator::fRandom
private

Random number generator to use creating fake data lists.

Definition at line 147 of file EventListManipulator.h.

Referenced by CreatePoissonFakeData(), and reconfigure().

std::map<novadaq::cnv::DetId, std::set<fnex::SelectionType_t> > fnex::EventListManipulator::fSelectionsToUse
private
bool fnex::EventListManipulator::fSumSubEpochs
private

do we want to sum all sub epochs into a single one

Definition at line 158 of file EventListManipulator.h.

Referenced by ExtractMetaDataFromFile(), and reconfigure().

art::ServiceHandle<art::TFileService> fnex::EventListManipulator::fTFS
private

TFileService.

Definition at line 144 of file EventListManipulator.h.

Referenced by FillFakeDataLists(), and FillHistograms().

double fnex::EventListManipulator::fTotalFDDataPOTFHC
private

total FD data POT FHC

Definition at line 168 of file EventListManipulator.h.

Referenced by ExtractMetaDataFromFile(), FarDetCosmicBackgroundHistAndScale(), and reconfigure().

double fnex::EventListManipulator::fTotalFDDataPOTRHC
private

total FD data POT RHC

Definition at line 169 of file EventListManipulator.h.

Referenced by ExtractMetaDataFromFile(), FarDetCosmicBackgroundHistAndScale(), and reconfigure().

double fnex::EventListManipulator::fTotalNDDataPOTFHC
private

total ND data POT FHC

Definition at line 170 of file EventListManipulator.h.

Referenced by ExtractMetaDataFromFile().

double fnex::EventListManipulator::fTotalNDDataPOTRHC
private

total ND data POT RHC

Definition at line 171 of file EventListManipulator.h.

Referenced by ExtractMetaDataFromFile().

std::vector<std::string> fnex::EventListManipulator::fTreeDirectories
private

directory holding the input trees

Definition at line 146 of file EventListManipulator.h.

Referenced by Deserialize(), FixNearDetMCExposures(), and reconfigure().

bool fnex::EventListManipulator::fUseCosmicHists
private

if true, use histograms to make cosmic background if false, use eventlisttrees of cosmic data

Definition at line 155 of file EventListManipulator.h.

Referenced by reconfigure().

bool fnex::EventListManipulator::fUseEventId
private

use if you need access to run/subrun/event/etc

Definition at line 157 of file EventListManipulator.h.

Referenced by FillEventList(), and reconfigure().

bool fnex::EventListManipulator::fUseExtrapolation
private

are we extrapolating or using covariance matrices

Definition at line 181 of file EventListManipulator.h.

Referenced by AddFakeDataToEventList(), Deserialize(), FillFakeDataLists(), FillHistogram(), and reconfigure().


The documentation for this class was generated from the following files: