compare_cos_numi.C
Go to the documentation of this file.
1 // to get histograms of cosmics from cosmic trigger data and numi stream, and GENIE prediction, by period and by quantile
2 // argument is period: 1, 2, 3, or 5
3 // author: Kirk Bays (bays@caltech.edu) Sep 7 2017
4 // Histograms made: cosrej BDT inputs, neutrino E, and optimized neutrino E
5 // suitable to run on the grid (set bool grid to true to keep track of pot/livetime and not scale per job)
6 // another important option is 'quantiletype': 0 means use quantiles determined per period, 1 means use quantiles determined from all periods combined
7 
8 #include "CAFAna/Cuts/Cuts.h"
11 #include "CAFAna/Cuts/SpillCuts.h"
12 #include "CAFAna/Cuts/TruthCuts.h"
14 #include "CAFAna/Core/Loaders.h"
15 #include "CAFAna/Core/Binning.h"
18 #include "CAFAna/Prediction/PredictionNumuFAHadE.h"
20 #include "CAFAna/Systs/Systs.h"
24 #include "CAFAna/Analysis/Calcs.h"
25 #include "CAFAna/Fit/Fit.h"
27 #include "CAFAna/Analysis/Plots.h"
34 #include "CAFAna/Vars/FitVars.h"
38 #include "CAFAna/FC/FCCollection.h"
39 #include "CAFAna/FC/FCPoint.h"
40 #include "CAFAna/FC/FCSurface.h"
41 
42 #include "OscLib/OscCalcGeneral.h"
43 #include "OscLib/OscCalcPMNS.h"
44 #include "OscLib/OscCalc.h"
45 #include "OscLib/IOscCalc.h"
46 
47 #include "TCanvas.h"
48 #include "TFile.h"
49 #include "TGraph.h"
50 #include "TH1.h"
51 #include "TH2.h"
52 #include "TMath.h"
53 #include "TGaxis.h"
54 #include "TMultiGraph.h"
55 #include "TLegend.h"
56 #include "TLatex.h"
57 #include "TStyle.h"
58 #include "THStack.h"
59 #include "TPaveText.h"
60 #include "TGraphAsymmErrors.h"
61 #include "TLegendEntry.h"
62 #include "TMarker.h"
63 
64 #include <cmath>
65 #include <iostream>
66 #include <vector>
67 #include <sstream>
68 #include <string>
69 #include <fstream>
70 #include <iomanip>
71 
72 using namespace ana;
73 
75 {
76  bool grid = false; // set to true if running on grid
77  bool quantiletype = true; // set to all periods combined
78 
79  double livetime, pot;
80  std::string cosfiles, numifiles, geniefiles, outfile, quantpath, epoch;
81 
82  if(period == 1) {
83  livetime = kSecondAnaPeriod1Livetime;
85  epoch = "1";
86  }
87 
88  else if(period == 2) {
89  livetime = kSecondAnaPeriod2Livetime;
91  epoch = "2";
92  }
93 
94  else if(period == 3) {
97  epoch = "3";
98  }
99 
100  else if(period == 5) {
101  livetime = kAna2017Period5Livetime;
102  pot = kAna2017Period5POT;
103  epoch = "5";
104  }
105 
106  else {
107  std::cout << "NOT A KNOWN PERIOD!!!! " << period << std::endl;
108  exit (EXIT_FAILURE);
109  }
110 
111  cosfiles = "prod_decaf_R17-03-01-prod3reco.h_fd_cosmic_period" + epoch + "_nue_or_numu_or_nus_contain_v1_goodruns";
112  numifiles = "prod_caf_R17-03-01-prod3reco.k_fd_numi_fhc_period" + epoch + "_v1_goodruns";
113  geniefiles = "prod_decaf_R17-03-01-prod3reco.l_fd_genie_fluxswap_fhc_nova_v08_period" + epoch + "_nue_or_numu_or_nus_contain_v1";
114  outfile = "p" + epoch + "_cosmic_hists.root";
115  if(quantiletype) quantpath = "/nova/app/users/bays/files/all_FD_histo_for_quantile_cuts.root";
116  else quantpath = "/nova/app/users/mbaird42/3A_FD_spectra/final_period" + epoch + "_FD_histo_for_quantile_cuts.root";
117 
118  SpectrumLoader loaderCosmic(cosfiles, kCosmic);
119  SpectrumLoader loaderNumi(numifiles, kCosmic);
120  SpectrumLoader loaderGenie(geniefiles);
121 
122  loaderCosmic.SetSpillCut(kStandardSpillCuts);
123  loaderNumi.SetSpillCut(kStandardSpillCuts);
124  loaderGenie.SetSpillCut(kStandardSpillCuts);
125 
126  const Var kMaxY([](const caf::SRProxy *sr) {
127  if(sr->trk.kalman.ntracks < 1)
128  return -10.0;
129  return max(sr->trk.kalman.tracks[0].start.Y()/100.0, sr->trk.kalman.tracks[0].stop.Y()/100.0);
130  });
131 
132 
133  const Var kCVN_cosmic([](const caf::SRProxy *sr)
134  {
135  return sr->sel.cvn.output[14];
136  });
137 
138  const Var kCellVar([](const caf::SRProxy *sr)
139  {
140  float kalvar = sr->sel.contain.kalfwdcell + sr->sel.contain.kalbakcell;
141  float cosvar = sr->sel.contain.cosfwdcell + sr->sel.contain.cosbakcell;
142  return min(kalvar,cosvar);
143  });
144 
145  const Var kHitRatio([](const caf::SRProxy *sr)
146  {
147  if(sr->trk.kalman.ntracks < 1)
148  return -5.0;
149  return (sr->trk.kalman.tracks[0].nhit/(int(sr->slc.nhit)*1.0));
150  });
151 
152 
153  TFile quantfile(quantpath.c_str(), "read");
154 
155  gDirectory->cd("dir_FDSpec2D");
156 
157  TH2* quanthist = (TH2*)quantfile.FindObjectAny("FDSpec2D");
158 
159  std::vector<Cut> QuantCuts = QuantileCutsFromTH2(quanthist, kNumuCCOptimisedAxis, kHadEFracAxis, 4);
160 
161  std::vector<HistAxis*> kHistAxes; // size 9
162 
163  kHistAxes.push_back(new HistAxis("Reconstructed Neutrino Energy (GeV)", kNumuEnergyBinning, kCCE));
164  kHistAxes.push_back(new HistAxis("Reconstructed Neutrino Energy (GeV)", kNumuCCEOptimisedBinning, kCCE));
165  kHistAxes.push_back(new HistAxis("cosine angle wrt beam", Binning::Simple(100,0,1), kNumuCosRejAngleKal));
166  kHistAxes.push_back(new HistAxis("kalman y-dir", Binning::Simple(100,-1,1), kDirY));
167  kHistAxes.push_back(new HistAxis("kalman track length", Binning::Simple(100,0,25), kTrkLength));
168  kHistAxes.push_back(new HistAxis("max of Kalman start/stop Y pos", Binning::Simple(100,-8,8), kMaxY));
169  kHistAxes.push_back(new HistAxis("CVN_cosmic", Binning::Simple(100,0,1), kCVN_cosmic));
170  kHistAxes.push_back(new HistAxis("min cos/kal fwd+bak", Binning::Simple(100,0,500), kCellVar));
171  kHistAxes.push_back(new HistAxis("kalman hits / slice hits", Binning::Simple(100,0,1), kHitRatio));
172 
173  std::vector<std::string> kType;
174  kType.push_back("Cos");
175  kType.push_back("Numi");
176  kType.push_back("Genie");
177 
178  std::vector<std::string> kHistName;
179  kHistName.push_back("numuE");
180  kHistName.push_back("numuEOpt");
181  kHistName.push_back("cosangle");
182  kHistName.push_back("kalydir");
183  kHistName.push_back("trklen");
184  kHistName.push_back("maxy");
185  kHistName.push_back("cvncos");
186  kHistName.push_back("cellvar");
187  kHistName.push_back("hitratio");
188 
189  std::vector<std::string> kQuantile;
190  kQuantile.push_back("all");
191  kQuantile.push_back("quant1");
192  kQuantile.push_back("quant2");
193  kQuantile.push_back("quant3");
194  kQuantile.push_back("quant4");
195 
196  std::vector<Cut> kCuts;
197 
198  const Cut cutCosmic = kNumuCutFD2017&&kInCosmicTimingWindow;
199  const Cut cutNumi = kNumuCutFD2017&&kInTimingSideband;
200  const Cut cutGenie = kNumuCutFD2017;
201 
202  kCuts.push_back(cutCosmic);
203  kCuts.push_back(cutNumi);
204  kCuts.push_back(cutGenie);
205 
206  std::vector<Var> kWeights;
207 
208  kWeights.push_back(kUnweighted);
209  kWeights.push_back(kUnweighted);
210  kWeights.push_back(kXSecCVWgt2017*kPPFXFluxCVWgt);
211 
212  TH1F* livetimecos = new TH1F("livetimecos","livetimecos",1,0,1);
213  TH1F* livetimenumi = new TH1F("livetimenumi","livetimenumi",1,0,1);
214  TH1F* potgenie = new TH1F("potgenie","potgenie",1,0,1);
215 
216  TH1F* countcos = new TH1F("countcos","countcos",1,0,1);
217  TH1F* countnumi = new TH1F("countnumi","countnumi",1,0,1);
218  TH1F* countgenie = new TH1F("countgenie","countgenie",1,0,1);
219 
220  std::vector< std::vector < std::vector <Spectrum*> > > kSpectra; // will be 3(cos/numi/genie) x 9(vars) x 5(all+4quantiles)
221 
222  std::vector < std::vector <Spectrum*> > kDoubleVecCos;
223  for(unsigned int j=0; j<9; j++) {
224  std::vector <Spectrum*> kSingleVec;
225  kSingleVec.push_back(new Spectrum(loaderCosmic, *(kHistAxes[j]), kCuts[0], kNoShift, kWeights[0]));
226  for(auto thisCut : QuantCuts){
227  kSingleVec.push_back(new Spectrum(loaderCosmic, *(kHistAxes[j]), kCuts[0]&&thisCut, kNoShift, kWeights[0]));
228  }
229  kDoubleVecCos.push_back(kSingleVec);
230  }
231  kSpectra.push_back(kDoubleVecCos);
232 
233  std::vector < std::vector <Spectrum*> > kDoubleVecNumi;
234  for(unsigned int j=0; j<9; j++) {
235  std::vector <Spectrum*> kSingleVec;
236  kSingleVec.push_back(new Spectrum(loaderNumi, *(kHistAxes[j]), kCuts[1], kNoShift, kWeights[1]));
237  for(auto thisCut : QuantCuts){
238  kSingleVec.push_back(new Spectrum(loaderNumi, *(kHistAxes[j]), kCuts[1]&&thisCut, kNoShift, kWeights[1]));
239  }
240  kDoubleVecNumi.push_back(kSingleVec);
241  }
242  kSpectra.push_back(kDoubleVecNumi);
243 
244  std::vector < std::vector <Spectrum*> > kDoubleVecGenie;
245  for(unsigned int j=0; j<9; j++) {
246  std::vector <Spectrum*> kSingleVec;
247  kSingleVec.push_back(new Spectrum(loaderGenie, *(kHistAxes[j]), kCuts[2], kNoShift, kWeights[2]));
248  for(auto thisCut : QuantCuts){
249  kSingleVec.push_back(new Spectrum(loaderGenie, *(kHistAxes[j]), kCuts[2]&&thisCut, kNoShift, kWeights[2]));
250  }
251  kDoubleVecGenie.push_back(kSingleVec);
252  }
253  kSpectra.push_back(kDoubleVecGenie);
254 
255  loaderCosmic.Go();
256  loaderNumi.Go();
257  loaderGenie.Go();
258 
259  std::vector < std::vector < std::vector <TH1*> > > kHists;
260 
261  std::vector < std::vector<TH1*> > kDoubleHistCos;
262  for(unsigned int j = 0; j < 9; j++) {
263  std::vector<TH1*> kSingleHists;
264  for(unsigned int k = 0; k<5; k++) {
265  if(grid) kSingleHists.push_back(kSpectra[0][j][k]->ToTH1(kSpectra[0][j][k]->Livetime(), kLivetime));
266  else kSingleHists.push_back(kSpectra[0][j][k]->ToTH1(livetime, kLivetime));
267  }
268  kDoubleHistCos.push_back(kSingleHists);
269  }
270  kHists.push_back(kDoubleHistCos);
271 
272  std::vector < std::vector<TH1*> > kDoubleHistNumi;
273  for(unsigned int j = 0; j < 9; j++) {
274  std::vector<TH1*> kSingleHists;
275  for(unsigned int k = 0; k<5; k++) {
276  if(grid) kSingleHists.push_back(kSpectra[1][j][k]->ToTH1(kSpectra[1][j][k]->Livetime(), kLivetime));
277  else kSingleHists.push_back(kSpectra[1][j][k]->ToTH1(livetime, kLivetime));
278  }
279  kDoubleHistNumi.push_back(kSingleHists);
280  }
281  kHists.push_back(kDoubleHistNumi);
282 
283  std::vector < std::vector<TH1*> > kDoubleHistGenie;
284  for(unsigned int j = 0; j < 9; j++) {
285  std::vector<TH1*> kSingleHists;
286  for(unsigned int k = 0; k<5; k++) {
287  if(grid) kSingleHists.push_back(kSpectra[2][j][k]->ToTH1(kSpectra[2][j][k]->POT(), kPOT));
288  else kSingleHists.push_back(kSpectra[2][j][k]->ToTH1(pot, kPOT));
289  }
290  kDoubleHistGenie.push_back(kSingleHists);
291  }
292  kHists.push_back(kDoubleHistGenie);
293 
294  livetimecos->SetBinContent(1,kSpectra[0][0][0]->Livetime());
295  livetimenumi->SetBinContent(1,kSpectra[1][0][0]->Livetime());
296  potgenie->SetBinContent(1,kSpectra[2][0][0]->POT());
297 
298  countcos->SetBinContent(1,(kSpectra[0][0][0]->ToTH1(kSpectra[0][0][0]->Livetime(),kLivetime))->GetEntries());
299  countnumi->SetBinContent(1,(kSpectra[1][0][0]->ToTH1(kSpectra[1][0][0]->Livetime(),kLivetime))->GetEntries());
300  countgenie->SetBinContent(1,(kSpectra[2][0][0]->ToTH1(kSpectra[2][0][0]->POT()))->GetEntries());
301 
302  TFile fout(outfile.c_str(), "recreate");
303 
304  for(unsigned int i=0; i<3; i++) {
305  for(unsigned int j=0; j<9; j++) {
306  for(unsigned int k=0; k<5; k++) {
307  std::string name = kType[i]+kHistName[j]+kQuantile[k];
308  kHists[i][j][k]->Write(name.c_str());
309  }
310  }
311  }
312 
313  if(grid) {
314  livetimecos->Write("livetimecos");
315  livetimenumi->Write("livetimenumi");
316  potgenie->Write("potgenie");
317  countcos->Write("countcos");
318  countnumi->Write("countnumi");
319  countgenie->Write("countgenie");
320  }
321 
322 
323  fout.Close();
324 
325 }
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
const Cut kNumuCutFD2017
Definition: NumuCuts2017.h:39
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 double kSecondAnaEpoch3dPOT
Definition: Exposures.h:77
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
const double kSecondAnaEpoch3bPOT
Definition: Exposures.h:75
caf::Proxy< int > cosfwdcell
Definition: SRProxy.h:810
void SetSpillCut(const SpillCut &cut)
caf::Proxy< int > cosbakcell
Definition: SRProxy.h:805
const double kSecondAnaPeriod1POT
Definition: Exposures.h:73
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
void compare_cos_numi(int period)
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
const double kAna2017Epoch3eLivetime
Definition: Exposures.h:197
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.
const double j
Definition: BetheBloch.cxx:29
caf::Proxy< std::vector< float > > output
Definition: SRProxy.h:909
const double kSecondAnaEpoch3cLivetime
Definition: Exposures.h:102
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
const double kSecondAnaEpoch3dLivetime
Definition: Exposures.h:103
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 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 double kSecondAnaEpoch3cPOT
Definition: Exposures.h:76
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
caf::Proxy< caf::SRIDBranch > sel
Definition: SRProxy.h:2141
const double kSecondAnaEpoch3bLivetime
Definition: Exposures.h:101
const Var kXSecCVWgt2017
Definition: XsecTunes.h:37
T max(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
Definition: statistics.h:68
const double kAna2017Epoch3ePOT
Definition: Exposures.h:170
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
FILE * outfile
Definition: dump_event.C:13
enum BeamMode string