SpectrumLoader.cxx
Go to the documentation of this file.
2 
3 #include "CAFAna/Core/Progress.h"
6 #include "CAFAna/Core/Spectrum.h"
8 
12 
14 
15 #include <cassert>
16 #include <iostream>
17 #include <cmath>
18 
19 #include "TFile.h"
20 #include "TH2.h"
21 #include "TTree.h"
22 
23 namespace ana
24 {
25  //----------------------------------------------------------------------
27  : SpectrumLoaderBase(wildcard)
28  {
29  }
30 
31  //----------------------------------------------------------------------
32  SpectrumLoader::SpectrumLoader(const std::vector<std::string>& fnames)
33  : SpectrumLoaderBase(fnames)
34  {
35  }
36 
37  //----------------------------------------------------------------------
39  {
40  }
41 
42  //----------------------------------------------------------------------
44  int fileLimit)
45  {
47  ret.fWildcard = "project "+proj;
48  ret.fFileSource = std::unique_ptr<IFileSource>(new SAMProjectSource(proj, fileLimit));
49  return ret;
50  }
51 
52  //----------------------------------------------------------------------
54  {
55  }
56 
57  struct CompareByID
58  {
59  bool operator()(const Cut& a, const Cut& b)
60  {
61  return a.ID() < b.ID();
62  }
63  };
64 
65  //----------------------------------------------------------------------
67  {
68  if(fGone){
69  std::cerr << "Error: can only call Go() once on a SpectrumLoader" << std::endl;
70  abort();
71  }
72  fGone = true;
73 
74  const int Nfiles = NFiles();
75 
76  Progress* prog = 0;
77 
79 
80  int fileIdx = -1;
81  while(TFile* f = GetNextFile()){
82  ++fileIdx;
83 
84  if(Nfiles >= 0 && !prog){
85  std::string sum = TString::Format("Filling %d spectra", fSink.NSinks()).Data();
86  if(Nfiles > 1){
87  if(!fSpillHistDefs.empty()){
88  sum += TString::Format(", %lu spillTree spectra", fSpillHistDefs.size()).Data();
89  }
90  if(!fNuHistDefs.empty()){
91  sum += TString::Format(", %lu nuTree spectra", fNuHistDefs.size()).Data();
92  }
93  }
94 
95  sum += TString::Format(" from %d files matching '%s'", Nfiles, fWildcard.c_str()).Data();
96  prog = new Progress(sum);
97  }
98 
99  HandleFile(f, Nfiles == 1 ? prog : 0);
100 
101  if(Nfiles > 1 && prog) prog->SetProgress((fileIdx+1.)/Nfiles);
102  } // end for fileIdx
103 
105 
106  if(prog){
107  prog->Done();
108  delete prog;
109  }
110 
111  std::cout << "Total POT before spill cuts: " << fRunPOT << std::endl;
113 
114  // Delete the tree of syst/cut/weight appliers
115  for(auto it: fSinksS) delete it.second;
116  for(auto it: fSinksSC) delete it.second;
117  for(auto it: fSinksSCW) delete it.second;
118  fSinksS.clear();
119  fSinksSC.clear();
120  fSinksSCW.clear();
121 
122  // Delete the actual spectrum sinks themselves
123  for(IRecordSink* s: fSpectrumSinks) delete s;
124  fSpectrumSinks.clear();
125  }
126 
127  //----------------------------------------------------------------------
129  {
130  assert(!f->IsZombie());
131 
132  // Test for flat (trees are inside a directory) or nested (tree is at top
133  // level) cases.
134  TDirectory* dir = 0;
135  TTree* tr = 0;
136 
137  TObject* obj = f->Get("recTree");
138  assert(obj); // Must exist in one form or the other
139 
140  // It might seem like you could use GetObject() to do the type-checking
141  // directly, but that method seems to have a memory leak in the case that
142  // the types don't match.
143  if(obj->ClassName() == std::string("TTree")){
144  // nested case
145  tr = (TTree*)obj;
146  }
147  else{
148  // FlatCAF case
149  dir = (TDirectory*)obj;
150  tr = (TTree*)dir->Get("rec");
151  assert(tr);
152  }
153 
154  long n;
155  caf::SRProxy sr(dir, tr, "rec", n, 0);
156 
157  // FloatingExceptionOnNaN fpnan;
158 
159  const long Nentries = tr->GetEntries();
160  for(n = 0; n < Nentries; ++n){
161  if(!dir) tr->LoadTree(n); // nested mode
162 
163  HandleRecord(&sr);
164 
165  if(prog && n%100 == 0) prog->SetProgress(double(n)/Nentries);
166  } // end for n
167  }
168 
169  //----------------------------------------------------------------------
171  {
172  // First up, check this spill passes the global quality cuts
173  if(fSpillCut && !(*fSpillCut)(&sr->spill)) return;
174 
175  // These aren't set in SRNeutrino, but we need them for flux systs, so copy
176  // them over from SRSpill.
177  //
178  // TODO TODO - remove this hack (and the equivalent in
179  // SpectrumLoaderBase.cxx) once they're available in the CAFs
180  if(sr->mc.nnu > 0){
181  sr->mc.nu[0].isFHC = bool(sr->spill.isFHC);
182  sr->mc.nu[0].is0HC = bool(sr->spill.is0HC);
183  sr->mc.nu[0].isRHC = bool(sr->spill.isRHC);
184  sr->mc.nu[0].det = caf::Det_t(sr->hdr.det);
185  }
186 
187  fSink.HandleRecord(sr, 1);
188  }
189 
190  //----------------------------------------------------------------------
192  {
194  }
195 
196 } // namespace
caf::Proxy< caf::SRSpill > spill
Definition: SRProxy.h:2102
Det_t
Which NOvA detector?
Definition: SREnums.h:7
caf::Proxy< bool > is0HC
Definition: SRProxy.h:1330
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
set< int >::iterator it
TFile * GetNextFile()
Forwards to fFileSource but also accumulates POT and livetime.
void HandleSpill(const caf::SRSpillProxy *spill) override
bool fGone
Has Go() been called? Can&#39;t add more histograms after that.
caf::Proxy< caf::SRHeader > hdr
Definition: SRProxy.h:2096
Proxy for caf::StandardRecord.
Definition: SRProxy.h:2085
caf::Proxy< std::vector< caf::SRNeutrino > > nu
Definition: SRProxy.h:574
double fRunPOT
Crude measure, not including spill cuts.
void AccumulateExposures(const caf::SRSpillProxy *spill) override
int NSinks() const override
OStream cerr
Definition: OStream.cxx:7
void HandleRecord(caf::SRProxy *sr, double wei) override
caf::Proxy< short int > nnu
Definition: SRProxy.h:573
std::unique_ptr< SpillCut > fSpillCut
Cut applied to the spill branch for every event.
std::set< IRecordSink * > fSpectrumSinks
caf::Proxy< bool > isFHC
Definition: SRProxy.h:1331
const XML_Char * s
Definition: expat.h:262
const double a
std::vector< std::string > wildcard(const std::string &wildcardString)
Definition: convert.C:9
std::unordered_map< int, SystApplier * > fSinksS
bool operator()(const Cut &a, const Cut &b)
int ID() const
Cuts with the same definition will have the same ID.
Definition: Cut.h:55
static SpectrumLoader FromSAMProject(const std::string &proj, int fileLimit=-1)
Named constructor for SAM projects.
caf::StandardRecord * sr
virtual void Go() override
Load all the registered spectra.
std::list< SpillHistDef > fSpillHistDefs
std::unordered_map< std::pair< int, int >, CutApplier * > fSinksSC
fileLimit
Definition: ProjMan.py:109
OStream cout
Definition: OStream.cxx:6
Base class for the various types of spectrum loader.
caf::Proxy< caf::SRTruthBranch > mc
Definition: SRProxy.h:2097
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
std::unique_ptr< IFileSource > fFileSource
void SetProgress(double frac)
Update the progress fraction between zero and one.
Definition: Progress.cxx:41
Collaborates with Spectrum and OscillatableSpectrum to fill spectra from CAF files.
TDirectory * dir
Definition: macro.C:5
virtual void HandleRecord(caf::SRProxy *sr)
Float_t proj
Definition: plot.C:35
const hit & b
Definition: hits.cxx:21
assert(nhit_max >=nhit_nbins)
A simple ascii-art progress bar.
Definition: Progress.h:9
Fetch files from a pre-existing SAM project.
virtual void HandleFile(TFile *f, Progress *prog=0)
std::unordered_map< std::tuple< int, int, int >, WeightApplier * > fSinksSCW
void Format(TGraph *gr, int lcol, int lsty, int lwid, int mcol, int msty, double msiz)
Definition: Style.cxx:154
caf::Proxy< bool > isRHC
Definition: SRProxy.h:1332
int NFiles() const
Forwards to fFileSource.
Double_t sum
Definition: plot.C:31
std::list< NuHistDef > fNuHistDefs
void Done()
Call this when action is completed.
Definition: Progress.cxx:90
void Finalize()
POT/livetime is not filled into sinks until this is called!
caf::Proxy< caf::Det_t > det
Definition: SRProxy.h:231
void ReportExposures(double denom=0) const