estimate_cosmics_triggeronly.C
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////////
2 // estimate_cosmics_triggeronly.C //
3 // //
4 // follows get_cosmic_sample.C //
5 // uses cosmic trigger sample for spectrum shape (more stats) //
6 // //
7 // Script by: Diana Mendez (d.p.mendez@sussex.ac.uk) April 4 2018 //
8 // //
9 ////////////////////////////////////////////////////////////////////////////
10 
12 
13 using namespace ana;
14 
15 
17 {
18 
19  std::cout << "\nEstimating " << horn << " cosmics" << std::endl;
20  std::cout << "========================== " << std::endl;
21  std::cout << " *\n *\n * FD" << std::endl;
22  std::cout << " *\n * " << std::endl;
23  std::cout << " * " << std::endl;
24  std::cout << "==========================\n " << std::endl;
25 
26 
27  double pot, livetime;
28  std::string period_numi = "full", period_cosmic = "full";
29 
30  if(horn=="fhc"){
31  livetime = kAna2018FHCLivetime;
32  pot = kAna2018FHCPOT;
33  }
34  else if(horn=="rhc"){
35  livetime = kAna2018RHCLivetime;
36  pot = kAna2018RHCPOT;
37  }
38  else{
39  std::cout << "=== unrecognised horn " << horn << "===" << std::endl;
40  exit (EXIT_FAILURE);
41  }
42 
43 
44  double pot_sample = kAna2017Period3POT;
45  double livetime_sample = kAna2017Period3Livetime;
46  // // in case of using period5
47  // pot_sample = kAna2017Period5POT;
48  // livetime_sample = kAna2017Period5Livetime;
49 
50 
51  std::string file_name_cosmic = "results/cosmiccosmic_" + horn + "_" + period_cosmic + "__numu2018.root";
52  TFile *fcosmic = TFile::Open(file_name_cosmic.c_str());
53 
54  std::vector<TH1*> hcosmic;
55  std::cout << "Getting histograms from files" << std::endl;
56  for(int quantId = 0; quantId < 5; quantId++){
57  hcosmic.push_back((TH1*)fcosmic->Get(Form("numuEOpt_q%d", quantId)));
58  }
59 
60 
61  std::cout << "Getting pot, livetime, count" << std::endl;
62  TH1* hlivetime_cosmic = (TH1*)fcosmic->Get("livetime");
63  float livetime_cosmic = hlivetime_cosmic->Integral();
64  // check counts agree in total
65  // and get counts per quantile
66 
67  std::cout << "Checking counts per quantile" << std::endl;
68  std::vector<TH1*> hcount_cosmic;
69  std::vector<float> count_cosmic;
70  for(int quantId = 0; quantId < 5; quantId++){
71  float counttemp = hcosmic[quantId]->Integral();
72  count_cosmic.push_back(counttemp);
73  }
74 
75  // note livetime scaling:
76  // livetime sideband is not equal to livetime from beamspill
77  float scale = livetime / livetime_cosmic;
78  std::cout << "\ncosmic trigger livetime = " << livetime_cosmic << ", ana livetime = " << livetime << " -> scale = " << scale << std::endl;
79 
80 
81  std::vector<float> events;
82  std::vector<float> errors;
83  for(int quantId=0; quantId<5; quantId++){
84  events.push_back( count_cosmic[quantId] * scale );
85  errors.push_back( events[quantId] * sqrt(count_cosmic[quantId]) / count_cosmic[quantId] );
86  if(quantId==0) std::cout << "\nTotal estimate = " << events[quantId] << " +/- " << errors[quantId] << " from " << count_cosmic[quantId] << " raw events" << std::endl;
87  else std::cout << "quantile " << quantId << " = " << events[quantId] << " +/- " << errors[quantId] << " from " << count_cosmic[quantId] << " raw events" << std::endl;
88  }
89 
90  // now scale the trigger to number of events
91  std::cout << "\nChecking cosmic integrals are right" << std::endl;
92  std::cout << "(should agree with event estimates)" << std::endl;
93  for(int quantId=0; quantId<5; quantId++){
94  hcosmic[quantId]->Scale( events[quantId] / (hcosmic[quantId]->Integral()) );
95  if(quantId==0) std::cout << "Total integral = " << hcosmic[quantId]->Integral() << std::endl;
96  else std::cout << "quantile " << quantId << " = " << hcosmic[quantId]->Integral() << std::endl;
97  }
98 
99  std::cout << "\nWriting to output file" << std::endl;
100  std::string name_out = "cosmicstriggeronly_" + horn + "__numu2018.root";
101  TFile fout(name_out.c_str(),"recreate");
102  for(int quantId=0; quantId<5; quantId++){
103  hcosmic[quantId]->Write(Form("cosmics_q%d", quantId));
104  }
105 
106  fout.Close();
107 
108  std::cout << "\nDone!" << std::endl;
109 
110 }
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
void estimate_cosmics_triggeronly(std::string horn="fhc")
double Integral(const Spectrum &s, const double pot, int cvnregion)
T sqrt(T number)
Definition: d0nt_math.hpp:156
const double kAna2017Period3Livetime
Definition: Exposures.h:199
const double kAna2018RHCPOT
Definition: Exposures.h:208
Double_t scale
Definition: plot.C:25
#define pot
OStream cout
Definition: OStream.cxx:6
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
double livetime
Definition: saveFDMCHists.C:21
exit(0)
const double kAna2018FHCPOT
Definition: Exposures.h:207
void events(int which)
Definition: Cana.C:52
const double kAna2017Period3POT
Definition: Exposures.h:172
const double kAna2018FHCLivetime
Definition: Exposures.h:213
const double kAna2018RHCLivetime
Definition: Exposures.h:214