estimate_cosmics.C
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////////
2 // estimate_cosmics.C //
3 // //
4 // follows get_cosmic_sample.C //
5 // uses cosmic trigger sample for spectrum shape (more stats) //
6 // uses numi data for the count (real fd conditions from sidebands) //
7 // if only using the cosmic trigger data, run *_triggeronly.C //
8 // //
9 // Script by: Diana Mendez (d.p.mendez@sussex.ac.uk) April 4 2018 //
10 // //
11 ////////////////////////////////////////////////////////////////////////////
12 
14 
15 using namespace ana;
16 
18 {
19 
20  std::cout << "\nEstimating " << horn << " cosmics" << std::endl;
21  std::cout << "========================== " << std::endl;
22  std::cout << " *\n *\n * FD" << std::endl;
23  std::cout << " *\n * " << std::endl;
24  std::cout << " * " << std::endl;
25  std::cout << "==========================\n " << std::endl;
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  std::string file_name_numi = "results/cosmicnumi_" + horn + "_" + period_cosmic + "__numu2018.root";
45  std::string file_name_cosmic = "results/cosmiccosmic_" + horn + "_" + period_cosmic + "__numu2018.root";
46  TFile *fnumi = TFile::Open(file_name_numi.c_str());
47  TFile *fcosmic = TFile::Open(file_name_cosmic.c_str());
48 
49 
50  // we'll write original numi and trigger histograms in the output
51  // so clone hcosmic for final estimate (instead of directly scaling the trigger)
52  std::vector<TH1*> hnumi;
53  std::vector<TH1*> hcosmic;
54  std::vector<TH1*> hestimate;
55  std::cout << "Getting histograms from files" << std::endl;
56  for(int quantId = 0; quantId < 5; quantId++){
57  hnumi.push_back((TH1*)fnumi->Get(Form("numuEOpt_q%d", quantId)));
58  hcosmic.push_back((TH1*)fcosmic->Get(Form("numuEOpt_q%d", quantId)));
59  TH1* htemp = (TH1*)hcosmic[quantId]->Clone(); hestimate.push_back(htemp);
60  }
61 
62 
63  // dont trust the count in the files (to do: check)
64  // get the count directly from the integral
65  std::cout << "Getting pot, livetime, count" << std::endl;
66  TH1* hlivetime_numi = (TH1*)fnumi->Get("livetime"); float livetime_numi = hlivetime_numi->Integral();
67  TH1* hlivetime_cosmic = (TH1*)fcosmic->Get("livetime"); float livetime_cosmic = hlivetime_cosmic->Integral();
68 
69  // note livetime scaling:
70  // livetime sideband is not equal to livetime from beamspill which is what analysis call livetime
71  float scale = livetime / livetime_numi;
72  std::cout << "\nsideband livetime = " << livetime_numi << ", ana livetime = " << livetime << " -> scale = " << scale << std::endl;
73 
74 
75  std::cout << "Getting counts per quantile" << std::endl;
76  std::vector<float> count_numi, count_cosmic;
77  for(int quantId = 0; quantId < 5; quantId++){
78  float counttemp;
79  counttemp = hnumi[quantId]->Integral(); count_numi.push_back(counttemp);
80  counttemp = hcosmic[quantId]->Integral(); count_cosmic.push_back(counttemp);
81  }
82 
83  std::vector<float> events;
84  std::vector<float> errors;
85  for(int quantId=0; quantId<5; quantId++){
86  events.push_back( count_numi[quantId] * scale );
87  errors.push_back( events[quantId] * sqrt(count_numi[quantId]) / count_numi[quantId] );
88  if(quantId==0) std::cout << "\nTotal estimate = " << events[quantId] << " +/- " << errors[quantId] << " from " << count_numi[quantId] << " raw events" << std::endl;
89  else std::cout << "quantile " << quantId << " = " << events[quantId] << " +/- " << errors[quantId] << " from " << count_numi[quantId] << " raw events" << std::endl;
90  }
91 
92  // we get the shape from the trigger and the count from numi
93  // now scale the trigger to number of events
94  std::cout << "\nChecking cosmic integrals are right" << std::endl;
95  std::cout << "(should agree with event estimates)" << std::endl;
96  for(int quantId=0; quantId<5; quantId++){
97  hestimate[quantId]->Scale( events[quantId] / (hestimate[quantId]->Integral()) );
98  if(quantId==0) std::cout << "Total integral = " << hestimate[quantId]->Integral() << std::endl;
99  else std::cout << "quantile " << quantId << " = " << hestimate[quantId]->Integral() << std::endl;
100  }
101 
102  std::cout << "\nWriting to output file" << std::endl;
103  std::string name_out = "cosmics_" + horn + "__numu2018.root";
104  TFile fout(name_out.c_str(),"recreate");
105  for(int quantId=0; quantId<5; quantId++){
106  hnumi[quantId]->Write(Form("numi_q%d", quantId));
107  hcosmic[quantId]->Write(Form("trigger_q%d", quantId));
108  hestimate[quantId]->Write(Form("cosmics_q%d", quantId));
109  }
110 
111  fout.Close();
112 
113  std::cout << "\nDone!" << std::endl;
114 
115 }
void estimate_cosmics(std::string horn="rhc")
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
double Integral(const Spectrum &s, const double pot, int cvnregion)
T sqrt(T number)
Definition: d0nt_math.hpp:156
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 kAna2018FHCLivetime
Definition: Exposures.h:213
const double kAna2018RHCLivetime
Definition: Exposures.h:214