52 TRandom *
RNG() {
return &fRng; }
65 : fKeepProb(keepProb), fLastHash(0)
70 return Keep(
Hash(sr));
88 if (hash != fLastHash)
94 return RNG()->Uniform() < fKeepProb;
113 KeepMCEventFunc(
double keepProb,
double NCwgt, std::map<std::pair<int, int>,
unsigned int> * kept =
nullptr)
114 :
KeepEvtFunc(keepProb), fNCwgt(NCwgt), fEventsKept(kept)
123 if (!Keep(&sr->
spill))
130 if (sr->
mc.
nu[0].iscc)
131 keep = RNG()->Uniform() < fCalc.P(sr->
mc.
nu[0].pdgorig, sr->
mc.
nu[0].pdg, sr->
mc.
nu[0].E);
133 keep = RNG()->Uniform() < fNCwgt;
135 if (keep && fEventsKept)
139 (*fEventsKept)[flav]++;
153 mutable std::map<std::pair<int, int>,
unsigned int> *
fEventsKept;
178 template <
typename T>
196 const Cut * decafCut;
202 fdDataDef =
"prod4_sumrestricteddecaf_fd_numi_rhc_full_goodruns_nue2018";
204 else if (sample ==
"numu")
208 fdDataDef =
"prod4_sumrestricteddecaf_fd_numi_rhc_full_goodruns_numu2018";
213 std::cerr <<
"I only know how to make fake data for 'numu' or 'nue' samples." <<
std::endl;
222 const std::map<Loaders::SwappingConfig, std::string> kSwaps
237 std::map<Loaders::SwappingConfig, Spectrum> mcDummySpecs;
238 std::map<Loaders::SwappingConfig, TH1D> mcDummySpillHists;
239 for (
const auto swapPair : kSwaps)
242 mcDummySpecs.emplace(std::piecewise_construct,
243 std::forward_as_tuple(swapPair.first),
244 std::forward_as_tuple(
loader, dummyAxis, *decafCut));
245 mcDummySpillHists.emplace(std::piecewise_construct,
246 std::forward_as_tuple(swapPair.first),
247 std::forward_as_tuple(Form(
"spills_%d", swapPair.first),
"", 1, 0, 2));
260 const double fakeDataLivetime = fakeDataPOT / dataDummySpec.POT() * dataDummySpec.Livetime();
261 double totalMCPOT = 0;
262 for (
const auto & swapPair : mcDummySpecs)
263 totalMCPOT += swapPair.second.POT();
268 <<
" total RHC POT, we need to draw " << fakeDataPOT <<
" POT of fake data.)" <<
std::endl;
269 for (
const auto & swapPair : mcDummySpecs)
271 <<
" (" << mcDummySpillHists.at(swapPair.first).Integral() <<
" spills)" <<
std::endl;
275 <<
" need " << fakeDataLivetime <<
"s of livetime in fake data." <<
std::endl;
278 std::map<Loaders::SwappingConfig, std::string> intermediateFileNames;
279 const double MCSwapScale = 1. / mcDummySpecs.size();
280 for (
const auto & swapPair : mcDummySpecs)
282 auto swap = swapPair.first;
292 double keepProbCC = fakeDataPOT / swapPair.second.POT();
293 const auto nSpills = mcDummySpillHists.at(swapPair.first).Integral();
295 if (keepProbSpill < 0 || keepProbSpill > 1)
309 <<
" since we will be sampling from multiple swaps.)" <<
std::endl;
312 <<
" spills to get desired livetime" << std::endl;
314 std::string intermediateFileName =
"2019fakedata.intermediate-" + kSwaps.at(
swap) +
".root";
315 intermediateFileNames[
swap] = intermediateFileName;
317 std::map<std::pair<int, int>,
unsigned int> eventsKept;
320 Cut kKeepEventCut{keepMcEventFunc};
325 MCsampler.AddEventCut(*decafCut);
327 double totalLive = 0;
332 bool potKeep =
false;
335 keepSpillFuncObj.SetKeepProb(keepProbCC * MCSwapScale);
336 potKeep = keepSpillFuncObj.Keep(spillProxy);
340 keepSpillFuncObj.SetKeepProb(keepProbSpill);
341 bool livetimeKeep = potKeep || keepSpillFuncObj.Keep(spillProxy);
360 if (!potKeep && !livetimeKeep)
379 TFile fakeEventsFile(intermediateFileName.c_str(),
"update");
380 auto potHist =
dynamic_cast<TH1 *
>(fakeEventsFile.Get(
"TotalPOT"));
381 potHist->SetBinContent(1, totalPOT);
383 auto liveHist =
dynamic_cast<TH1 *
>(fakeEventsFile.Get(
"TotalLivetime"));
384 liveHist->SetBinContent(1, totalLive);
387 auto evtTree =
dynamic_cast<TTree *
>(fakeEventsFile.Get(
"recTree"));
388 auto evtHist =
dynamic_cast<TH1 *
>(fakeEventsFile.Get(
"TotalEvents"));
389 std::size_t evts = evtTree->GetEntries();
390 evtHist->SetBinContent(1, evts);
392 fakeEventsFile.Close();
396 <<
"% of the desired amount." <<
std::endl;
398 for (
const auto & keptPair : eventsKept)
401 std::pair<int, int> flav = keptPair.first;
402 unsigned int kept = keptPair.second;
408 for (
const auto &
v : {flav.first, flav.second})
424 out = out.substr(0, out.size() - 5);
436 std::vector<std::string>
files;
437 std::transform(intermediateFileNames.begin(), intermediateFileNames.end(), std::back_inserter(files),
438 [](
const auto & swapFilePair){
return swapFilePair.second; }
441 files.reserve(files.size() + samquerysource.NFiles());
442 files.insert(files.end(), samquerysource.GetFileNames().begin(), samquerysource.GetFileNames().end());
452 OOTKeepProb *= fakeDataPOT / dataDummySpec.POT();
462 if (!keepCosmicFunc.Keep(srp))
468 std::deque<std::reference_wrapper<float>> toUpdate
475 toUpdate.push_back(vtx.time);
477 for (
auto &
val : toUpdate)
483 for (
const auto & fileNamePair : intermediateFileNames)
484 system( (
"rm " + fileNamePair.second).c_str() );
caf::Proxy< caf::SRSpill > spill
Shared items for fake data studies for 2019 analysis.
Far Detector at Ash River.
void SetSeed(std::size_t seed)
Cuts and Vars for the 2020 FD DiF Study.
void ResetOscCalcToDefault(osc::IOscCalcAdjustable *calc)
Reset calculator to default assumptions for all parameters.
system("rm -rf microbeam.root")
std::map< std::pair< int, int >, unsigned int > * fEventsKept
const Cut kApplySecondAnalysisMask([](const caf::SRProxy *sr){if(sr->hdr.det!=caf::kFARDET) return true; std::bitset< 14 > binary(sr->hdr.dibmask);std::pair< int, int > planesA=calcFirstLastLivePlane(sr->slc.firstplane, binary);std::pair< int, int > planesB=calcFirstLastLivePlane(sr->slc.lastplane, binary);if((planesA.first!=planesB.first)||(planesA.second!=planesB.second)) return false;return((planesA.second-planesA.first+1)/64 >=4);})
const int kBeamWindowMicroSec
How long is the beam window?
const int kBeamWindowMinMicroSec
virtual bool Keep(const caf::SRProxy *sr)
caf::Proxy< caf::SRHeader > hdr
caf::Proxy< float > spillpot
Proxy for caf::StandardRecord.
caf::Proxy< std::vector< caf::SRNeutrino > > nu
virtual bool Keep(const caf::SRSpillProxy *sr)
float starttime
start time [ns]
const Cut kNueDQ2017CVN([](const caf::SRProxy *sr){if(sr->sel.nuecosrej.hitsperplane >=8) return false;if(!sr->vtx.elastic.IsValid) return false;if(sr->vtx.elastic.fuzzyk.npng==0) return false;return true;})
caf::Proxy< short int > nnu
std::pair< Spectrum *, CheatDecomp * > make_pair(SpectrumLoaderBase &loader_data, SpectrumLoaderBase &loader_mc, HistAxis *axis, Cut *cut, const SystShifts &shift, const Var &wei)
KeepMCEventFunc(double keepProb, double NCwgt, std::map< std::pair< int, int >, unsigned int > *kept=nullptr)
const int kMinTimingSidebandAfterMicroSec
const int kBeamWindowMaxMicroSec
Representation of a spectrum in any variable, with associated POT.
void Go()
Call Go() on all the loaders.
const SpillVar kSpillUnweighted
void AddReductionStep(const std::function< ReductionFunc > &f)
Run the specified reduction function over each event.
std::string GetLoaderPath(caf::Det_t det, DataMC datamc, DataSource src=kBeam, SwappingConfig swap=kNonSwap) const
float endtime
end time [ns]
void SetSpillCut(const SpillCut &cut)
const double FULL_RHC_EXPOSURE
KeepEvtFunc(double keepProb)
caf::Proxy< caf::SRCVNResult > cvn
SpectrumLoaderBase & GetLoader(caf::Det_t det, DataMC datamc, DataSource src=kBeam, SwappingConfig swap=kNonSwap)
Retrieve a specific loader.
caf::Proxy< float > nueid
const double kFixedSpillLivetime
bool isgoodspill
Was the pot for a spill good? (only applicable to data, default true)
void AddEventCut(const Cut &cut)
Only copy records to the output file if they pass this cut.
Optimized version of OscCalcPMNS.
printf("%d Experimental points found\n", nlines)
const Cut kNumuDecafCutFD2018
ColorStream(std::ostream &os, unsigned int col=36)
File source based on a SAM query or dataset (definition)
virtual void Go() override
Load all the registered spectra.
const Var kCVNSSe([](const caf::SRProxy *sr){throw std::runtime_error("kCVNSSe is no longer available. Fix your macro so you don't use it.");return-5.;})
2018 nue PID
bool Keep(std::size_t hash)
float livetime
Length of readout [s].
caf::Proxy< caf::SRTruthBranch > mc
bool operator()(const caf::SRProxy *sr)
The StandardRecord is the primary top-level object in the Common Analysis File trees.
std::size_t Hash(const std::vector< std::size_t > &vals)
Generate a unique hash from a sequence of integers.
const std::string OUT_FILE_STUB
SRElastic elastic
Single vertex found by Elastic Arms.
void generate_fd_fake_events(const std::string &sample)
const Cut kInBeamSpill([](const caf::SRProxy *sr){if(sr->spill.run > util::kLastBadTimingRun) return kInBeamSpill_main(sr);else return kInBeamSpill_main(sr)||kInBeamSpill_shifted(sr);}, [](const caf::SRSpillProxy *spill){if(spill->run > util::kLastBadTimingRun) return kInBeamSpill_main.Livetime(spill);else return kInBeamSpill_main.Livetime(spill)+kInBeamSpill_shifted.Livetime(spill);}, [](const caf::SRSpillProxy *spill) -> double{return spill->spillpot;})
Does the event fall inside the window we call the beam spill?
assert(nhit_max >=nhit_nbins)
SRSlice slc
Slice branch: nhit, extents, time, etc.
std::ostream & operator<<(const T &t)
const SpillCut kStandardSpillCuts
Apply this unless you're doing something special.
base class: random number generator with easily settable seed...
caf::Proxy< caf::SRIDBranch > sel
std::string to_string(ModuleType mt)
For nominal spectra and reweighting systs (xsec/flux)
osc::OscCalcPMNSOpt fCalc
void SetLoaderPath(const std::string &path, caf::Det_t det, DataMC datamc, DataSource src=kBeam, SwappingConfig swap=kNonSwap)
Configure loader via wildcard path.
const Cut kNueFD2018DecafCut
static Binning Simple(int n, double lo, double hi, const std::vector< std::string > &labels={})
bool ismc
data or MC? True if MC
const Var kCVN2017e([](const caf::SRProxy *sr){return sr->sel.cvn.nueid;})
const Var kUnweighted
The simplest possible Var, always 1. Used as a default weight.
float meantime
mean time, weighted by charge [ns]
SRVertexBranch vtx
Vertex branch: location, time, etc.
const int kMaxTimingSidebandAfterMicroSec
void SetKeepProb(double prob)