get_cosmic_sample.C
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 // get_cosmic_sample.C //
3 // //
4 // get histograms of reco neutrino energy + cosmic bdt inputs //
5 // from cosmic trigger data and mc prediction //
6 // string set = {numi, cosmic, genie} //
7 // bool fhc = {true=fhc,false=rhc} //
8 // bool grid = true to to keep track of pot/livetime and not scale per job //
9 // //
10 // based on code by: Kirk Bays (bays@caltech.edu) Mar 7 2018 //
11 // modified by: Diana Mendez (d.p.mendez@sussex.ac.uk) April 4 2018 //
12 // //
13 ///////////////////////////////////////////////////////////////////////////////
14 
15 // notes for usage with numi:
16 // * concats currently only available for rhc
17 
18 
20 
21 using namespace ana;
22 
23 
24 void get_cosmic_sample(std::string set = "numi",
25  std::string period = "period1",
26  std::string horn = "fhc",
27  bool grid = true)
28 {
29 
30  double livetime, pot;
31  std::string prodcos, prodgen, setcos;
32  std::string loadfiles, numifiles, cosmicfiles, geniefiles;
33  std::string outfile, quantpath;
34 
35 
36  if(horn=="fhc"){
37  livetime = kAna2018FHCLivetime;
38  pot = kAna2018FHCPOT;
39  prodcos = ".b";
40  prodgen = ".d";
41  setcos = "full";
42  numifiles = "prod4_caf_fd_numi_" + horn + "_" + period + "_goodruns";
43  }
44  else if(horn=="rhc"){
45  livetime = kAna2018RHCLivetime;
46  pot = kAna2018RHCPOT;
47  prodcos = ".a";
48  prodgen = ".e";
49  setcos = "HighGain";
50  numifiles = "prod4_sumdecaf_fd_numi_" + horn + "_" + period + "_goodruns_numu2018";
51  }
52  else{
53  std::cout << "=== unrecognised horn " << horn << "===" << std::endl;
54  exit (EXIT_FAILURE);
55  }
56 
57 
58  // this is basically only for numi
59  // have to think about a neater way to do this
60  // horn is an argument but redefine here for safety
61  if(period=="period1"){
62  horn = "fhc";
63  livetime = kSecondAnaPeriod1Livetime;
65  }
66  else if(period=="period2"){
67  horn = "fhc";
68  livetime = kSecondAnaPeriod2Livetime;
70  }
71  else if(period=="period3"){
72  horn = "fhc";
73  livetime = kAna2017Period3Livetime;
74  pot = kAna2017Period3POT;
75  }
76  else if(period=="period4"){
77  horn = "rhc";
78  livetime = kAna2018Period4Livetime;
79  pot = kAna2018Period4POT;
80  }
81  else if(period=="period5"){
82  horn = "fhc";
83  livetime = kAna2017Period5Livetime;
84  pot = kAna2017Period5POT;
85  }
86  else if(period=="period6"){
87  horn = "rhc";
88  livetime = kAna2018Period6Livetime;
89  pot = kAna2018Period6POT;
90  }
91  else if(period=="epoch7a"){
92  horn = "rhc";
93  livetime = kAna2018Period7Livetime;
94  pot = kAna2018Period7POT;
95  }
96  else if(period=="epoch7b"){
97  horn = "rhc";
98  livetime = kAna2018Period7Livetime;
99  pot = kAna2018Period7POT;
100  }
101  /*
102  else{
103  std::cout << "=== unrecognised period " << period << "===" << std::endl;
104  exit (EXIT_FAILURE);
105  }
106  */
107  quantpath = "/pnfs/nova/persistent/analysis/numu/Ana2018/provisional/quantiles/quantiles__" + horn + "_full__numu2018.root";
108  cosmicfiles = "prod_sumdecaf_R17-11-14-prod4reco" + prodcos + "_fd_cosmic_" + horn + "Tune_" + setcos + "_v1_goodruns_numu2018";
109  geniefiles = "prod_sumdecaf_R17-11-14-prod4reco" + prodgen + "_fd_genie_nonswap_" + horn + "_nova_v08_full_v1_numu2018";
110  outfile = "cosmic" + set + "_" + horn + "_" + period + "__numu2018.root";
111 
112 
113  std::cout << "Defining files, cut and weight" << std::endl;
114  Cut kTempCut = kNoCut;
115  Var kTempWeight = kUnweighted;
116  if(set=="numi"){
117  loadfiles = numifiles;
118  kTempCut = kNumuCutFD2018&&kInTimingSideband;
119  }
120  else if(set=="cosmic"){
121  loadfiles = cosmicfiles;
123  }
124  else if(set=="genie"){
125  loadfiles = geniefiles;
126  kTempCut = kNumuCutFD2018;
127  kTempWeight = kXSecCVWgt2018*kPPFXFluxCVWgt;
128  }
129  else {
130  std::cout << "UNRECOGNISED SAMPLE" << set << std::endl;
131  std::cout << "First argument should be:" << std::endl;
132  std::cout << "numi, cosmic, genie" << std::endl;
133  exit (EXIT_FAILURE);
134  }
135  const Cut kCutSet = kTempCut;
136  const Var kWeightSet = kTempWeight;
137 
138 
139  // ==================================================================
140  // Display settings
141  std::cout << "beam " << horn << ", period " << period << std::endl;
142  std::cout << "sample " << set << std::endl;
143  std::cout << "input " << loadfiles << std::endl;
144  // ==================================================================
145 
146 
147  // DataSource no longer needed for loader
148  std::cout << "Load here" << std::endl;
149  SpectrumLoader loader(loadfiles);
151 
152 
153  const Var kMaxY([](const caf::SRProxy *sr) {
154  if(sr->trk.kalman.ntracks < 1)
155  return -10.0;
156  return max(sr->trk.kalman.tracks[0].start.Y()/100.0, sr->trk.kalman.tracks[0].stop.Y()/100.0);
157  });
158 
159 
160  const Var kCVN_cosmic([](const caf::SRProxy *sr)
161  {
162  return sr->sel.cvn.output[390];
163  //return sr->sel.cvn.output[14];
164  });
165 
166  const Var kCellVar([](const caf::SRProxy *sr)
167  {
168  float kalvar = sr->sel.contain.kalfwdcell + sr->sel.contain.kalbakcell;
169  float cosvar = sr->sel.contain.cosfwdcell + sr->sel.contain.cosbakcell;
170  return min(kalvar,cosvar);
171  });
172 
173  const Var kHitRatio([](const caf::SRProxy *sr)
174  {
175  if(sr->trk.kalman.ntracks < 1)
176  return -5.0;
177  return (sr->trk.kalman.tracks[0].nhit/(int(sr->slc.nhit)*1.0));
178  });
179 
180 
181 
182  std::cout << "Histogram axes" << std::endl;
183  std::vector<HistAxis*> kHistAxes; // size 8
184  kHistAxes.push_back(new HistAxis("Reconstructed Neutrino Energy (GeV)", kNumuEnergyBinning, kCCE));
185  kHistAxes.push_back(new HistAxis("Reconstructed Neutrino Energy (GeV)", kNumuCCEOptimisedBinning, kCCE));
186  kHistAxes.push_back(new HistAxis("cosine angle wrt beam", Binning::Simple(100,0,1), kNumuCosRejAngleKal));
187  kHistAxes.push_back(new HistAxis("kalman y-dir", Binning::Simple(100,-1,1), kDirY));
188  kHistAxes.push_back(new HistAxis("kalman track length", Binning::Simple(100,0,25), kTrkLength));
189  kHistAxes.push_back(new HistAxis("max of Kalman start/stop Y pos", Binning::Simple(100,-8,8), kMaxY));
190  kHistAxes.push_back(new HistAxis("min cos/kal fwd+bak", Binning::Simple(100,0,500), kCellVar));
191  kHistAxes.push_back(new HistAxis("kalman hits / slice hits", Binning::Simple(100,0,1), kHitRatio));
192 
193 
194  std::cout << "Histogram names" << std::endl;
195  std::vector<std::string> kHistName;
196  kHistName.push_back("numuE");
197  kHistName.push_back("numuEOpt");
198  kHistName.push_back("cosangle");
199  kHistName.push_back("kalydir");
200  kHistName.push_back("trklen");
201  kHistName.push_back("maxy");
202  kHistName.push_back("cellvar");
203  kHistName.push_back("hitratio");
204 
205 
206  std::cout << "Quantile names" << std::endl;
207  std::vector<std::string> kQuantile;
208  kQuantile.push_back("q0");
209  kQuantile.push_back("q1");
210  kQuantile.push_back("q2");
211  kQuantile.push_back("q3");
212  kQuantile.push_back("q4");
213 
214 
215  std::cout << "\nGet quantile cuts from file" << std::endl;
216  TFile* quantfile = TFile::Open(pnfs2xrootd(quantpath).c_str());
217  TH2* quanthist = (TH2*)quantfile->FindObjectAny("FDSpec2D");
218  std::vector<Cut> QuantCuts = QuantileCutsFromTH2(quanthist, kNumuCCOptimisedAxis, kHadEFracAxis, 4);
219  quantfile->Close();
220  delete quantfile;
221 
222 
223  std::cout << "\nDeclare crazy spectra" << std::endl;
224  std::vector < std::vector <Spectrum*> > kSpectra; // 8(vars) x 5(all+4quantiles)
225  for(unsigned int varId=0; varId<8; varId++) {
226  std::vector <Spectrum*> kSingleVec;
227  kSingleVec.push_back(new Spectrum(loader, *(kHistAxes[varId]), kCutSet, kNoShift, kWeightSet));
228  for(auto thisCut : QuantCuts){
229  kSingleVec.push_back(new Spectrum(loader, *(kHistAxes[varId]), kCutSet&&thisCut, kNoShift, kWeightSet));
230  }
231  kSpectra.push_back(kSingleVec);
232  }
233 
234 
235  std::cout << "Go loader, go!" << std::endl;
236  loader.Go();
237 
238 
239  std::cout << "\nMake crazy histogram vector" << std::endl;
240  std::vector < std::vector <TH1*> > kHists;
241  for(unsigned int varId=0; varId<8; varId++) {
242  std::vector<TH1*> kSingleHists;
243  for(unsigned int quantId=0; quantId<5; quantId++) {
244  if(set=="genie"){
245  if(grid) kSingleHists.push_back(kSpectra[varId][quantId]->ToTH1(kSpectra[varId][quantId]->POT(), kPOT));
246  else kSingleHists.push_back(kSpectra[varId][quantId]->ToTH1(pot, kPOT));
247  }
248  else{
249  if(grid) kSingleHists.push_back(kSpectra[varId][quantId]->ToTH1(kSpectra[varId][quantId]->Livetime(), kLivetime));
250  else kSingleHists.push_back(kSpectra[varId][quantId]->ToTH1(livetime, kLivetime));
251  }
252  }
253  kHists.push_back(kSingleHists);
254  }
255 
256 
257  std::cout << "Histos for pot count livetime" << std::endl;
258  TH1F* hpot = new TH1F("pot","pot",1,0,1);
259  TH1F* hcount = new TH1F("count","count",1,0,1);
260  TH1F* hlivetime = new TH1F("livetime","livetime",1,0,1);
261 
262 
263  std::cout << "Get pot count livetime" << std::endl;
264  if(set=="genie") hcount->SetBinContent(1,(kSpectra[0][0]->ToTH1(kSpectra[0][0]->POT()))->GetEntries());
265  else hcount->SetBinContent(1,(kSpectra[0][0]->ToTH1(kSpectra[0][0]->Livetime(),kLivetime))->GetEntries());
266  hlivetime->SetBinContent(1,kSpectra[0][0]->Livetime());
267  hpot->SetBinContent(1,kSpectra[0][0]->POT());
268 
269  std::cout << "Generate output" << std::endl;
270  TFile fout(outfile.c_str(), "recreate");
271 
272  if(grid) {
273  std::cout << "Write pot count livetime" << std::endl;
274  hlivetime ->Write("livetime");
275  hcount ->Write("count");
276  hpot ->Write("pot");
277  }
278 
279  std::cout << "Write histos" << std::endl;
280  for(unsigned int varId=0; varId<8; varId++) {
281  for(unsigned int quantId=0; quantId<5; quantId++) {
282  std::string name = kHistName[varId]+"_"+kQuantile[quantId];
283  kHists[varId][quantId]->Write(name.c_str());
284  }
285  }
286 
287 
288  std::cout << "\n\n=== Done! Be happy, go get a bubble tea ===\n\n" << std::endl;
289  fout.Close();
290 
291 }
const XML_Char * name
Definition: expat.h:151
_HistAxis< Var > HistAxis
Definition: HistAxis.h:103
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
const double kSecondAnaPeriod2POT
Definition: Exposures.h:74
const Var kPPFXFluxCVWgt
Definition: PPFXWeights.h:16
caf::Proxy< size_t > ntracks
Definition: SRProxy.h:1778
const Var kNumuCosRejAngleKal
Definition: NumuVars.cxx:547
caf::Proxy< caf::SRContain > contain
Definition: SRProxy.h:1251
Proxy for caf::StandardRecord.
Definition: SRProxy.h:2126
const HistAxis kHadEFracAxis("E_{had.} / E_{#nu}", Binning::Simple(200, 0, 1), kHadEFrac)
HistAxis that implements Hadronic Energy fraction binning used by L.Vinton to derive Hadronic Energy ...
Definition: HistAxes.h:30
const double kAna2017Period5POT
Definition: Exposures.h:171
caf::Proxy< int > cosfwdcell
Definition: SRProxy.h:810
void SetSpillCut(const SpillCut &cut)
caf::Proxy< int > cosbakcell
Definition: SRProxy.h:805
std::string pnfs2xrootd(std::string loc, bool unauth)
Definition: UtilsExt.cxx:237
const double kSecondAnaPeriod1POT
Definition: Exposures.h:73
const double kAna2017Period3Livetime
Definition: Exposures.h:199
const Var kTrkLength([](const caf::SRProxy *sr){if(sr->trk.kalman.ntracks< 1) return-5.f;return sr->trk.kalman.tracks[0].len/100;})
Definition: NumuVars.h:65
cout<< t1-> GetEntries()
Definition: plottest35.C:29
const Var kDirY([](const caf::SRProxy *sr){if(sr->trk.kalman.ntracks< 1) return-5.f;return sr->trk.kalman.tracks[0].dir.Y();})
Definition: NumuVars.h:38
const double kAna2018RHCPOT
Definition: Exposures.h:208
const Cut kInTimingSideband([](const caf::SRProxy *sr){if(sr->spill.run > util::kLastBadTimingRun) return(kInTimingSideband_before(sr)|| kInTimingSideband_after(sr));else return(kInTimingSideband_before(sr)|| kInTimingSideband_afterA(sr)|| kInTimingSideband_afterB(sr));}, [](const caf::SRSpillProxy *spill){if(spill->run > util::kLastBadTimingRun) return(kInTimingSideband_before.Livetime(spill)+ kInTimingSideband_after.Livetime(spill));else return(kInTimingSideband_before.Livetime(spill)+ kInTimingSideband_afterA.Livetime(spill)+ kInTimingSideband_afterB.Livetime(spill));}, [](const caf::SRSpillProxy *spill){return 0;})
Definition: TimingCuts.h:12
const Cut kInCosmicTimingWindow
Is the event far from the start and ends of the spill ? For FD cosmic selection.
Definition: TimingCuts.cxx:165
caf::Proxy< unsigned int > nhit
Definition: SRProxy.h:1315
const Var kMaxY([](const caf::SRProxy *sr){if(!sr->vtx.elastic.IsValid) return-1000.0f;if(sr->vtx.elastic.fuzzyk.nshwlid==0) return-1000.0f;float maxyall=-999.0;for(unsigned int i=0;i< sr->vtx.elastic.fuzzyk.nshwlid;i++){maxyall=std::max(std::max(sr->vtx.elastic.fuzzyk.png[i].shwlid.start.Y(), sr->vtx.elastic.fuzzyk.png[i].shwlid.stop.Y()), maxyall);}return maxyall;})
Definition: NueVars.h:60
caf::Proxy< caf::SRTrackBranch > trk
Definition: SRProxy.h:2145
const Binning kNumuEnergyBinning
Definition: Binnings.cxx:13
caf::Proxy< caf::SRCVNResult > cvn
Definition: SRProxy.h:1253
const Var kCCE
Definition: NumuVars.h:21
#define pot
caf::StandardRecord * sr
virtual void Go() override
Load all the registered spectra.
caf::Proxy< std::vector< float > > output
Definition: SRProxy.h:909
loader
Definition: demo0.py:10
const double kAna2017Period5Livetime
Definition: Exposures.h:198
std::vector< float > Spectrum
Definition: Constants.h:610
const HistAxis kNumuCCOptimisedAxis("Reconstructed Neutrino Energy (GeV)", kNumuCCEOptimisedBinning, kCCE)
HistAxis that implements optimised numuCCE from L. Vinton. See docdb 16332. This was close to &#39;custC&#39;...
Definition: HistAxes.h:24
const SystShifts kNoShift
Definition: SystShifts.cxx:21
OStream cout
Definition: OStream.cxx:6
std::vector< double > POT
static float min(const float a, const float b, const float c)
Definition: absgeo.cxx:45
caf::Proxy< int > kalbakcell
Definition: SRProxy.h:817
caf::Proxy< int > kalfwdcell
Definition: SRProxy.h:822
const double kSecondAnaPeriod1Livetime
Definition: Exposures.h:99
caf::Proxy< caf::SRKalman > kalman
Definition: SRProxy.h:1797
caf::Proxy< caf::SRSlice > slc
Definition: SRProxy.h:2142
Collaborates with Spectrum and OscillatableSpectrum to fill spectra from CAF files.
double livetime
Definition: saveFDMCHists.C:21
const Var kXSecCVWgt2018
Definition: XsecTunes.h:49
const Binning kNumuCCEOptimisedBinning
Optimised binning for numuCCE from L. Vinton. See docdb 16332. This was close to &#39;custC&#39; in that talk...
Definition: Binnings.cxx:28
exit(0)
std::vector< Cut > QuantileCutsFromTH2(TH2 *quantileHist, const HistAxis &independentAxis, const HistAxis &quantileAxis, const unsigned int &numQuantiles, const bool verbose)
: Do the same as the QuantileCuts function but taking in the TH2 instead of making it...
const Cut kNumuCutFD2018
Definition: NumuCuts2018.h:39
const SpillCut kStandardSpillCuts
Apply this unless you&#39;re doing something special.
Definition: SpillCuts.h:49
const double kSecondAnaPeriod2Livetime
Definition: Exposures.h:100
caf::Proxy< std::vector< caf::SRKalmanTrack > > tracks
Definition: SRProxy.h:1780
const double kAna2018FHCPOT
Definition: Exposures.h:207
caf::Proxy< caf::SRIDBranch > sel
Definition: SRProxy.h:2141
T max(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
Definition: statistics.h:68
static Binning Simple(int n, double lo, double hi, const std::vector< std::string > &labels={})
Definition: Binning.cxx:107
const Var kUnweighted
The simplest possible Var, always 1. Used as a default weight.
Definition: Var.h:96
const double kAna2017Period3POT
Definition: Exposures.h:172
void get_cosmic_sample(std::string set="numi", std::string period="period1", std::string horn="fhc", bool grid=true)
const Cut kNoCut
The simplest possible cut: pass everything, used as a default.
Definition: Cut.h:109
const double kAna2018FHCLivetime
Definition: Exposures.h:213
FILE * outfile
Definition: dump_event.C:13
const double kAna2018RHCLivetime
Definition: Exposures.h:214
enum BeamMode string