MakeNus17CosBkgd.C
Go to the documentation of this file.
1 // To get spectra of cosmics from cosmic trigger data and numi stream,
2 // and GENIE prediction, by period
3 // argument is period: 1, 2, 3, or 5
4 // Based on compare_cos_numi.C by Kirk Bays.
5 // Histograms made: NC Visible Energy, CVN NC, CVN Cosmic, NC BDT
6 // suitable to run on the grid (set bool grid to true to keep track of pot/livetime and not scale per job)
7 
8 #include "CAFAna/Cuts/Cuts.h"
10 #include "NuXAna/Cuts/NusCuts17.h"
11 #include "CAFAna/Cuts/SpillCuts.h"
12 #include "CAFAna/Cuts/TruthCuts.h"
13 #include "CAFAna/Core/Loaders.h"
14 #include "CAFAna/Core/Binning.h"
17 #include "CAFAna/Prediction/PredictionNumuFAHadE.h"
19 #include "CAFAna/Systs/Systs.h"
23 #include "CAFAna/Analysis/Calcs.h"
24 #include "CAFAna/Analysis/Plots.h"
30 #include "CAFAna/Vars/FitVars.h"
31 #include "NuXAna/Vars/NusVars.h"
35 
36 #include "OscLib/OscCalcGeneral.h"
37 #include "OscLib/OscCalcPMNS.h"
38 #include "OscLib/OscCalc.h"
39 #include "OscLib/IOscCalc.h"
40 
41 #include "TCanvas.h"
42 #include "TFile.h"
43 #include "TGraph.h"
44 #include "TH1.h"
45 #include "TH2.h"
46 #include "TMath.h"
47 #include "TGaxis.h"
48 #include "TMultiGraph.h"
49 #include "TLegend.h"
50 #include "TLatex.h"
51 #include "TStyle.h"
52 #include "THStack.h"
53 #include "TPaveText.h"
54 #include "TGraphAsymmErrors.h"
55 #include "TLegendEntry.h"
56 #include "TMarker.h"
57 
58 #include <cmath>
59 #include <iostream>
60 #include <vector>
61 #include <sstream>
62 #include <string>
63 #include <fstream>
64 #include <iomanip>
65 
66 using namespace ana;
67 
69 {
70 
71  double livetime, pot;
72  std::string cosfiles, numifiles, geniefiles, outfile, quantpath, epoch;
73 
74  if(period == 1) {
75  livetime = kSecondAnaPeriod1Livetime;
77  epoch = "1";
78  }
79 
80  else if(period == 2) {
81  livetime = kSecondAnaPeriod2Livetime;
83  epoch = "2";
84  }
85 
86  else if(period == 3) {
89  epoch = "3";
90  }
91 
92  else if(period == 5) {
93  livetime = kAna2017Period5Livetime;
94  pot = kAna2017Period5POT;
95  epoch = "5";
96  }
97 
98  else {
99  std::cout << "NOT A KNOWN PERIOD!!!! " << period << std::endl;
100  exit (EXIT_FAILURE);
101  }
102 
103  cosfiles = "prod_decaf_R17-03-01-prod3reco.h_fd_cosmic_period" + epoch + "_nue_or_numu_or_nus_contain_v1_goodruns";
104  numifiles = "prod_caf_R17-03-01-prod3reco.k_fd_numi_fhc_period" + epoch + "_v1_goodruns";
105  geniefiles = "prod_decaf_R17-03-01-prod3reco.l_fd_genie_fluxswap_fhc_nova_v08_period" + epoch + "_nue_or_numu_or_nus_contain_v1";
106  outfile = "p" + epoch + "_cosmic_hists.root";
107 
108 
109  bool grid = true; // set to true if running on grid
110 
111  SpectrumLoader loaderCosmic(cosfiles, kCosmic);
112  SpectrumLoader loaderNumi(numifiles, kCosmic);
113  SpectrumLoader loaderGenie(geniefiles);
114 
115  loaderCosmic.SetSpillCut(kStandardSpillCuts);
116  loaderNumi.SetSpillCut(kStandardSpillCuts);
117  loaderGenie.SetSpillCut(kStandardSpillCuts);
118 
119  const Var kCVN_cosmic([](const caf::SRProxy *sr)
120  {
121  return sr->sel.cvn.output[14];
122  });
123 
124 
125  std::vector<HistAxis*> kHistAxes; // size 4
126 
127  kHistAxes.push_back(new HistAxis("Visible Energy (GeV)", kNus17EnergyBinning, kNus17Energy));
128  kHistAxes.push_back(new HistAxis("CVN Cosmic Classifier", Binning::Simple(100,0,1), kCVN_cosmic));
129  kHistAxes.push_back(new HistAxis("CVN NC Classifier", Binning::Simple(100,0,1), kCVNnc));
130  kHistAxes.push_back(new HistAxis("NC Cosmic Rejection BDT", Binning::Simple(100,0,1), kNCCosRejAlt));
131 
132  std::vector<std::string> kType;
133  kType.push_back("Cos");
134  kType.push_back("Numi");
135  kType.push_back("Genie");
136 
137  std::vector<std::string> kHistName;
138  kHistName.push_back("nusE");
139  kHistName.push_back("cvncos");
140  kHistName.push_back("cvnnc");
141  kHistName.push_back("nccosrej");
142 
143  std::vector<Cut> kCuts;
144 
145  const Cut cutCosmic = kNus17FDAlt && kInCosmicTimingWindow;
146  const Cut cutNumi = kNus17FDAlt && kInTimingSideband;
147  const Cut cutGenie = kNus17FDAlt;
148 
149  kCuts.push_back(cutCosmic);
150  kCuts.push_back(cutNumi);
151  kCuts.push_back(cutGenie);
152 
153  std::vector<Var> kWeights;
154 
155  kWeights.push_back(kUnweighted);
156  kWeights.push_back(kUnweighted);
157  kWeights.push_back(kXSecCVWgt2017*kPPFXFluxCVWgt);
158 
159  TH1F* livetimecos = new TH1F("livetimecos","livetimecos",1,0,1);
160  TH1F* livetimenumi = new TH1F("livetimenumi","livetimenumi",1,0,1);
161  TH1F* potgenie = new TH1F("potgenie","potgenie",1,0,1);
162 
163  TH1F* countcos = new TH1F("countcos","countcos",1,0,1);
164  TH1F* countnumi = new TH1F("countnumi","countnumi",1,0,1);
165  TH1F* countgenie = new TH1F("countgenie","countgenie",1,0,1);
166 
167  std::vector< std::vector < std::vector<Spectrum*> > > kSpectra; // will be 3(cos/numi/genie) x 9(vars) x 1(all)
168 
169  const unsigned int nvars = 4;
170  std::vector < std::vector<Spectrum*> > kDoubleVecCos;
171  for(unsigned int j = 0; j < nvars; j++) {
172  std::vector <Spectrum*> kSingleVec;
173  kSingleVec.push_back(new Spectrum(loaderCosmic, *(kHistAxes[j]), kCuts[0], kNoShift, kWeights[0]));
174 
175  kDoubleVecCos.push_back(kSingleVec);
176  }
177  kSpectra.push_back(kDoubleVecCos);
178 
179  std::vector <std::vector<Spectrum*> > kDoubleVecNumi;
180  for(unsigned int j = 0; j < nvars; j++) {
181  std::vector<Spectrum*> kSingleVec;
182  kSingleVec.push_back(new Spectrum(loaderNumi, *(kHistAxes[j]), kCuts[1], kNoShift, kWeights[1]));
183  kDoubleVecNumi.push_back(kSingleVec);
184  }
185  kSpectra.push_back(kDoubleVecNumi);
186 
187  std::vector < std::vector <Spectrum*> > kDoubleVecGenie;
188  for(unsigned int j = 0; j < nvars; j++) {
189  std::vector <Spectrum*> kSingleVec;
190  kSingleVec.push_back(new Spectrum(loaderGenie, *(kHistAxes[j]), kCuts[2], kNoShift, kWeights[2]));
191  kDoubleVecGenie.push_back(kSingleVec);
192  }
193  kSpectra.push_back(kDoubleVecGenie);
194 
195  loaderCosmic.Go();
196  loaderNumi.Go();
197  loaderGenie.Go();
198 
199  std::vector < std::vector < std::vector <TH1*> > > kHists;
200 
201  std::vector < std::vector<TH1*> > kDoubleHistCos;
202  for(unsigned int j = 0; j < nvars; j++) {
203  std::vector<TH1*> kSingleHists;
204  for(unsigned int k = 0; k < 1; k++) {
205  if(grid) kSingleHists.push_back(kSpectra[0][j][k]->ToTH1(kSpectra[0][j][k]->Livetime(), kLivetime));
206  else kSingleHists.push_back(kSpectra[0][j][k]->ToTH1(livetime, kLivetime));
207  }
208  kDoubleHistCos.push_back(kSingleHists);
209  }
210  kHists.push_back(kDoubleHistCos);
211 
212  std::vector < std::vector<TH1*> > kDoubleHistNumi;
213  for(unsigned int j = 0; j < nvars; j++) {
214  std::vector<TH1*> kSingleHists;
215  for(unsigned int k = 0; k < 1; k++) {
216  if(grid) kSingleHists.push_back(kSpectra[1][j][k]->ToTH1(kSpectra[1][j][k]->Livetime(), kLivetime));
217  else kSingleHists.push_back(kSpectra[1][j][k]->ToTH1(livetime, kLivetime));
218  }
219  kDoubleHistNumi.push_back(kSingleHists);
220  }
221  kHists.push_back(kDoubleHistNumi);
222 
223  std::vector < std::vector<TH1*> > kDoubleHistGenie;
224  for(unsigned int j = 0; j < nvars; j++) {
225  std::vector<TH1*> kSingleHists;
226  for(unsigned int k = 0; k < 1; k++) {
227  if(grid) kSingleHists.push_back(kSpectra[2][j][k]->ToTH1(kSpectra[2][j][k]->POT(), kPOT));
228  else kSingleHists.push_back(kSpectra[2][j][k]->ToTH1(pot, kPOT));
229  }
230  kDoubleHistGenie.push_back(kSingleHists);
231  }
232  kHists.push_back(kDoubleHistGenie);
233 
234  livetimecos->SetBinContent(1,kSpectra[0][0][0]->Livetime());
235  livetimenumi->SetBinContent(1,kSpectra[1][0][0]->Livetime());
236  potgenie->SetBinContent(1,kSpectra[2][0][0]->POT());
237 
238  countcos->SetBinContent(1,(kSpectra[0][0][0]->ToTH1(kSpectra[0][0][0]->Livetime(),kLivetime))->GetEntries());
239  countnumi->SetBinContent(1,(kSpectra[1][0][0]->ToTH1(kSpectra[1][0][0]->Livetime(),kLivetime))->GetEntries());
240  countgenie->SetBinContent(1,(kSpectra[2][0][0]->ToTH1(kSpectra[2][0][0]->POT()))->GetEntries());
241 
242  TFile fout(outfile.c_str(), "recreate");
243 
244  for(unsigned int i = 0; i < 3; i++) {
245  for(unsigned int j = 0; j < nvars; j++) {
246  for(unsigned int k = 0; k < 1; k++) {
247  std::string name = kType[i]+kHistName[j];
248  kHists[i][j][k]->Write(name.c_str());
249  }
250  }
251  }
252 
253  if(grid) {
254  livetimecos->Write("livetimecos");
255  livetimenumi->Write("livetimenumi");
256  potgenie->Write("potgenie");
257  countcos->Write("countcos");
258  countnumi->Write("countnumi");
259  potgenie->Write("gotgenie");
260  }
261 
262 
263  fout.Close();
264 
265 }
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 kCVNnc
PID
Definition: Vars.cxx:44
const Var kPPFXFluxCVWgt
Definition: PPFXWeights.h:16
Proxy for caf::StandardRecord.
Definition: SRProxy.h:2126
const double kSecondAnaEpoch3dPOT
Definition: Exposures.h:77
const double kAna2017Period5POT
Definition: Exposures.h:171
const double kSecondAnaEpoch3bPOT
Definition: Exposures.h:75
void SetSpillCut(const SpillCut &cut)
const double kSecondAnaPeriod1POT
Definition: Exposures.h:73
cout<< t1-> GetEntries()
Definition: plottest35.C:29
const Var kNCCosRejAlt
Definition: NusVarsTemp.h:57
const Var kNus17Energy([](const caf::SRProxy *sr){double cale=sr->slc.calE;double recoE=0.;if(sr->hdr.det==caf::kFARDET) recoE=FDscaleCalE17 *cale;if(sr->hdr.det==caf::kNEARDET) recoE=NDscaleCalE17 *cale;return recoE;})
Definition: NusVars.h:62
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
const double kAna2017Epoch3eLivetime
Definition: Exposures.h:197
caf::Proxy< caf::SRCVNResult > cvn
Definition: SRProxy.h:1253
#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 SystShifts kNoShift
Definition: SystShifts.cxx:21
OStream cout
Definition: OStream.cxx:6
std::vector< double > POT
const double kSecondAnaPeriod1Livetime
Definition: Exposures.h:99
const double kSecondAnaEpoch3dLivetime
Definition: Exposures.h:103
const int nvars
Collaborates with Spectrum and OscillatableSpectrum to fill spectra from CAF files.
double livetime
Definition: saveFDMCHists.C:21
exit(0)
const Cut kNus17FDAlt
Definition: NusCuts17.h:92
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< caf::SRIDBranch > sel
Definition: SRProxy.h:2141
const double kSecondAnaEpoch3bLivetime
Definition: Exposures.h:101
const Var kXSecCVWgt2017
Definition: XsecTunes.h:37
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
const Binning kNus17EnergyBinning
Energy binnings used in Nus17 and planning for Nus18 for nus analysis.
Definition: Binning.cxx:71
FILE * outfile
Definition: dump_event.C:13
void MakeNus17CosBkgd(int period)
enum BeamMode string