28 std::string period_numi =
"full", period_cosmic =
"full";
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());
52 std::vector<TH1*> hnumi;
53 std::vector<TH1*> hcosmic;
54 std::vector<TH1*> hestimate;
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);
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();
71 float scale = livetime / livetime_numi;
72 std::cout <<
"\nsideband livetime = " << livetime_numi <<
", ana livetime = " << livetime <<
" -> scale = " << scale <<
std::endl;
76 std::vector<float> count_numi, count_cosmic;
77 for(
int quantId = 0; quantId < 5; quantId++){
79 counttemp = hnumi[quantId]->Integral(); count_numi.push_back(counttemp);
80 counttemp = hcosmic[quantId]->Integral(); count_cosmic.push_back(counttemp);
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;
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;
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));
void estimate_cosmics(std::string horn="rhc")
Cuts and Vars for the 2020 FD DiF Study.
double Integral(const Spectrum &s, const double pot, int cvnregion)
const double kAna2018RHCPOT
const double kAna2018FHCPOT
const double kAna2018FHCLivetime
const double kAna2018RHCLivetime