MakeCosmicBackground.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, 4, 5, 6, 7
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/NusCuts18.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, hc;
73 
74  if (period == 1) {
75  livetime = kSecondAnaPeriod1Livetime;
77  epoch = "1";
78  hc = "fhc";
79  }
80 
81  else if (period == 2) {
82  livetime = kSecondAnaPeriod2Livetime;
84  epoch = "2";
85  hc = "fhc";
86  }
87 
88  else if (period == 3) {
91  epoch = "3";
92  hc = "fhc";
93  }
94 
95  else if (period == 4) {
96  livetime = kAna2018Period4Livetime;
97  pot = kAna2018Period4POT;
98  epoch = "4";
99  hc = "rhc";
100  }
101 
102  else if (period == 5) {
103  livetime = kAna2017Period5Livetime;
104  pot = kAna2017Period5POT;
105  epoch = "5";
106  hc = "fhc";
107  }
108 
109  else if (period == 6) {
110  livetime = kAna2018Period6Livetime;
111  pot = kAna2018Period6POT;
112  epoch = "6";
113  hc = "rhc";
114  }
115 
116  else if (period == 7) {
117  livetime = kAna2018Period7Livetime;
118  pot = kAna2018Period7POT;
119  }
120 
121  else {
122  std::cout << "NOT A KNOWN PERIOD!!!! " << period << std::endl;
123  exit (EXIT_FAILURE);
124  }
125 
126  if (period == 1 or period == 2 or period == 3 or period == 5)
127  cosfiles = "prod_caf_R17-11-14-prod4reco.b_fd_cosmic_fhcTune_period" + epoch + "_v1_goodruns";
128  else if (period == 4)
129  cosfiles = "prod_caf_R17-11-14-prod4reco.a_fd_cosmic_rhcTune_period4_v1_goodruns";
130  else
131  cosfiles = "prod_caf_R17-11-14-prod4reco.b_fd_cosmic_fhcTune_period5_v1_goodruns";
132  if (period == 7)
133  numifiles = "prod4_caf_fd_numi_rhc_epoch7c";
134  else
135  numifiles = "prod4_caf_fd_numi_" + hc + "_period" + epoch + "_goodruns";
136  //geniefiles = "prod_caf_R17-11-14-prod4reco.d_fd_genie_fluxswap_fhc_nova_v08_period" + epoch + "_v1";
137  outfile = "p" + epoch + "_cosmic_hists.root";
138 
139 
140  bool grid = true; // set to true if running on grid
141 
142  SpectrumLoader loaderCosmic(cosfiles, kCosmic);
143  SpectrumLoader loaderNumi(numifiles, kCosmic);
144  //SpectrumLoader loaderGenie(geniefiles);
145 
146  loaderCosmic.SetSpillCut(kStandardSpillCuts);
147  loaderNumi.SetSpillCut(kStandardSpillCuts);
148  //loaderGenie.SetSpillCut(kStandardSpillCuts);
149 
150  const Var kCVN_cosmic([](const caf::SRProxy *sr) {
151  return sr->sel.cvnProd3Train.output[14];
152  });
153 
154  std::vector<HistAxis*> kHistAxes; // size 4
155 
156  kHistAxes.push_back(new HistAxis("Visible Energy (GeV)", kNus18EnergyBinning, kNus18Energy));
157  kHistAxes.push_back(new HistAxis("CVN Cosmic Classifier", Binning::Simple(100,0,1), kCVN_cosmic));
158  kHistAxes.push_back(new HistAxis("CVN NC Classifier", Binning::Simple(100,0,1), kCVNnc));
159  kHistAxes.push_back(new HistAxis("NC Cosmic Rejection BDT", Binning::Simple(100,0,1), kNCCosRejAlt));
160 
161  std::vector<std::string> kType;
162  kType.push_back("Cos");
163  kType.push_back("Numi");
164  //kType.push_back("Genie");
165 
166  std::vector<std::string> kHistName;
167  kHistName.push_back("nusE");
168  kHistName.push_back("cvncos");
169  kHistName.push_back("cvnnc");
170  kHistName.push_back("nccosrej");
171 
172  std::vector<Cut> kCuts;
173 
174  const Cut cutCosmic = kNus18FD && kInCosmicTimingWindow;
175  const Cut cutNumi = kNus18FD && kInTimingSideband;
176  //const Cut cutGenie = kNus18FD;
177 
178  kCuts.push_back(cutCosmic);
179  kCuts.push_back(cutNumi);
180  //kCuts.push_back(cutGenie);
181 
182  std::vector<Var> kWeights;
183 
184  kWeights.push_back(kUnweighted);
185  kWeights.push_back(kUnweighted);
186  //kWeights.push_back(kXSecCVWgt2018*kPPFXFluxCVWgt);
187 
188  TH1F* livetimecos = new TH1F("livetimecos","livetimecos",1,0,1);
189  TH1F* livetimenumi = new TH1F("livetimenumi","livetimenumi",1,0,1);
190  //TH1F* potgenie = new TH1F("potgenie","potgenie",1,0,1);
191 
192  TH1F* countcos = new TH1F("countcos","countcos",1,0,1);
193  TH1F* countnumi = new TH1F("countnumi","countnumi",1,0,1);
194  //TH1F* countgenie = new TH1F("countgenie","countgenie",1,0,1);
195 
196  std::vector< std::vector < std::vector<Spectrum*> > > kSpectra; // will be 3(cos/numi/genie) x 9(vars) x 1(all)
197 
198  const unsigned int nvars = 4;
199  std::vector < std::vector<Spectrum*> > kDoubleVecCos;
200  for(unsigned int j = 0; j < nvars; j++) {
201  std::vector <Spectrum*> kSingleVec;
202  kSingleVec.push_back(new Spectrum(loaderCosmic, *(kHistAxes[j]), kCuts[0], kNoShift, kWeights[0]));
203  kDoubleVecCos.push_back(kSingleVec);
204  }
205  kSpectra.push_back(kDoubleVecCos);
206 
207  std::vector <std::vector<Spectrum*> > kDoubleVecNumi;
208  for(unsigned int j = 0; j < nvars; j++) {
209  std::vector<Spectrum*> kSingleVec;
210  kSingleVec.push_back(new Spectrum(loaderNumi, *(kHistAxes[j]), kCuts[1], kNoShift, kWeights[1]));
211  kDoubleVecNumi.push_back(kSingleVec);
212  }
213  kSpectra.push_back(kDoubleVecNumi);
214 
215  // std::vector < std::vector <Spectrum*> > kDoubleVecGenie;
216  // for(unsigned int j = 0; j < nvars; j++) {
217  // std::vector <Spectrum*> kSingleVec;
218  // kSingleVec.push_back(new Spectrum(loaderGenie, *(kHistAxes[j]), kCuts[2], kNoShift, kWeights[2]));
219  // kDoubleVecGenie.push_back(kSingleVec);
220  // }
221  // kSpectra.push_back(kDoubleVecGenie);
222 
223  if (!(period == 6 or period == 7))
224  loaderCosmic.Go();
225  loaderNumi.Go();
226  //loaderGenie.Go();
227 
228  std::vector < std::vector < std::vector <TH1*> > > kHists;
229 
230  std::vector < std::vector<TH1*> > kDoubleHistCos;
231  for(unsigned int j = 0; j < nvars; j++) {
232  std::vector<TH1*> kSingleHists;
233  for(unsigned int k = 0; k < 1; k++) {
234  if(grid) kSingleHists.push_back(kSpectra[0][j][k]->ToTH1(kSpectra[0][j][k]->Livetime(), kLivetime));
235  else kSingleHists.push_back(kSpectra[0][j][k]->ToTH1(livetime, kLivetime));
236  }
237  kDoubleHistCos.push_back(kSingleHists);
238  }
239  kHists.push_back(kDoubleHistCos);
240 
241  std::vector < std::vector<TH1*> > kDoubleHistNumi;
242  for(unsigned int j = 0; j < nvars; j++) {
243  std::vector<TH1*> kSingleHists;
244  for(unsigned int k = 0; k < 1; k++) {
245  if(grid) kSingleHists.push_back(kSpectra[1][j][k]->ToTH1(kSpectra[1][j][k]->Livetime(), kLivetime));
246  else kSingleHists.push_back(kSpectra[1][j][k]->ToTH1(livetime, kLivetime));
247  }
248  kDoubleHistNumi.push_back(kSingleHists);
249  }
250  kHists.push_back(kDoubleHistNumi);
251 
252  // std::vector < std::vector<TH1*> > kDoubleHistGenie;
253  // for(unsigned int j = 0; j < nvars; j++) {
254  // std::vector<TH1*> kSingleHists;
255  // for(unsigned int k = 0; k < 1; k++) {
256  // if(grid) kSingleHists.push_back(kSpectra[2][j][k]->ToTH1(kSpectra[2][j][k]->POT(), kPOT));
257  // else kSingleHists.push_back(kSpectra[2][j][k]->ToTH1(pot, kPOT));
258  // }
259  // kDoubleHistGenie.push_back(kSingleHists);
260  // }
261  // kHists.push_back(kDoubleHistGenie);
262 
263  livetimecos->SetBinContent(1,kSpectra[0][0][0]->Livetime());
264  livetimenumi->SetBinContent(1,kSpectra[1][0][0]->Livetime());
265  //potgenie->SetBinContent(1,kSpectra[2][0][0]->POT());
266 
267  countcos->SetBinContent(1,(kSpectra[0][0][0]->ToTH1(kSpectra[0][0][0]->Livetime(),kLivetime))->GetEntries());
268  countnumi->SetBinContent(1,(kSpectra[1][0][0]->ToTH1(kSpectra[1][0][0]->Livetime(),kLivetime))->GetEntries());
269  //countgenie->SetBinContent(1,(kSpectra[2][0][0]->ToTH1(kSpectra[2][0][0]->POT()))->GetEntries());
270 
271  TFile fout(outfile.c_str(), "recreate");
272 
273  for(unsigned int i = 0; i < 2; i++) {
274  for(unsigned int j = 0; j < nvars; j++) {
275  for(unsigned int k = 0; k < 1; k++) {
276  std::string name = kType[i]+kHistName[j];
277  kHists[i][j][k]->Write(name.c_str());
278  }
279  }
280  }
281 
282  if(grid) {
283  livetimecos->Write("livetimecos");
284  livetimenumi->Write("livetimenumi");
285  //potgenie->Write("potgenie");
286  countcos->Write("countcos");
287  countnumi->Write("countnumi");
288  //potgenie->Write("gotgenie");
289  }
290 
291 
292  fout.Close();
293 
294 }
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 Var kNus18Energy([](const caf::SRProxy *sr){bool h_FHC=sr->spill.isFHC;bool h_RHC=sr->spill.isRHC;double cale=sr->slc.calE;double recoE=0.;if(h_FHC &&sr->hdr.det==caf::kFARDET) recoE=FDscaleCalE18 *cale;else if(h_FHC &&sr->hdr.det==caf::kNEARDET) recoE=NDscaleCalE18 *cale;else if(h_RHC &&sr->hdr.det==caf::kFARDET) recoE=FDscaleCalE18RHC *cale;else if(h_RHC &&sr->hdr.det==caf::kNEARDET) recoE=NDscaleCalE18RHC *cale;return recoE;})
Definition: NusVars.h:63
const double kSecondAnaPeriod2POT
Definition: Exposures.h:74
const Var kCVNnc
PID
Definition: Vars.cxx:44
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 Cut kNus18FD
Definition: NusCuts18.h:100
const Var kNCCosRejAlt
Definition: NusVarsTemp.h:57
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 Binning kNus18EnergyBinning
Energy binnings for Nus18 for nus analysis.
Definition: Binning.cxx:74
const double kAna2017Epoch3eLivetime
Definition: Exposures.h:197
#define pot
caf::StandardRecord * sr
virtual void Go() override
Load all the registered spectra.
const double j
Definition: BetheBloch.cxx:29
const double kSecondAnaEpoch3cLivetime
Definition: Exposures.h:102
const double kAna2017Period5Livetime
Definition: Exposures.h:198
void MakeNus18CosBkgd(int period)
std::vector< float > Spectrum
Definition: Constants.h:610
const SystShifts kNoShift
Definition: SystShifts.cxx:21
OStream cout
Definition: OStream.cxx:6
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 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 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