make_theorist_text_file.C
Go to the documentation of this file.
6 using namespace ana;
7 
9 
10 #include "TFile.h"
11 #include "TH2.h"
12 
13 #include <iomanip>
14 #include <iostream>
15 #include <vector>
16 
18 {
19  std::cout << "With dmsq21 = " << calc->GetDmsq21() << " eV^2, dmsq32 = " << calc->GetDmsq32() << " eV^2" << std::endl;
20  std::cout << "theta12 = " << calc->GetTh12() << ", theta23 = " << calc->GetTh23() << std::endl;
21  std::cout << "theta13 = " << calc->GetTh13() << ", delta = " << calc->GetdCP() << std::endl;
22  std::cout << "rho = " << calc->GetRho() << " g/cm^3, Z/A = 0.5, L = " << calc->GetL() << " km" << std::endl;
23  std::cout << "Total number of expected events: " << pred->Predict(calc).Integral(3.45e20)+nCosmic << std::endl;
24 
25  // For debugging differences to Alex's x-check
26  /*
27  std::cout << "NC " << pred->PredictComponent(calc, Flavors::kAll, Current::kNC, Sign::kBoth).Integral(3.45e20) << std::endl;
28  for(auto sign: {Sign::kNu, Sign::kAntiNu}){
29  std::cout << "mu->mu " << pred->PredictComponent(calc, Flavors::kNuMuToNuMu, Current::kCC, sign).Integral(3.45e20) << std::endl;
30  std::cout << "mu->e " << pred->PredictComponent(calc, Flavors::kNuMuToNuE, Current::kCC, sign).Integral(3.45e20) << std::endl;
31  std::cout << "e->mu " << pred->PredictComponent(calc, Flavors::kNuEToNuMu, Current::kCC, sign).Integral(3.45e20) << std::endl;
32  std::cout << "e->e " << pred->PredictComponent(calc, Flavors::kNuEToNuE, Current::kCC, sign).Integral(3.45e20) << std::endl;
33  std::cout << "mu->tau " << pred->PredictComponent(calc, Flavors::kNuMuToNuTau, Current::kCC, sign).Integral(3.45e20) << std::endl;
34  std::cout << "e->tau " << pred->PredictComponent(calc, Flavors::kNuEToNuTau, Current::kCC, sign).Integral(3.45e20) << std::endl;
35  }
36  */
37 }
38 
40 {
41  assert(pid == "LEM" || pid == "LID");
42 
43  TFile* f = new TFile(predfile.c_str());
44  assert(!f->IsZombie());
45 
46  std::unique_ptr<IPrediction> pred = LoadFromFile<IPrediction>(predfile, "prediction"+pid);
47  std::unique_ptr<IExtrap> extrap = LoadFromFile<IExtrap>(predfile, "prediction"+pid+"/nom/nom/extrap");
48 
49  TH1* ee = extrap->NueSurvComponent().TrueEnergy().ToTH1(3.45e20);
50  TH1* em = extrap->NumuAppComponent().TrueEnergy().ToTH1(3.45e20);
51  TH1* me = extrap->NueAppComponent().TrueEnergy().ToTH1(3.45e20);
52  TH1* mm = extrap->NumuSurvComponent().TrueEnergy().ToTH1(3.45e20);
53  TH1* mt = extrap->TauFromMuComponent().TrueEnergy().ToTH1(3.45e20);
54  TH1* et = extrap->TauFromEComponent().TrueEnergy().ToTH1(3.45e20);
55 
56  TH1* eebar = extrap->AntiNueSurvComponent().TrueEnergy().ToTH1(3.45e20);
57  TH1* embar = extrap->AntiNumuAppComponent().TrueEnergy().ToTH1(3.45e20);
58  TH1* mebar = extrap->AntiNueAppComponent().TrueEnergy().ToTH1(3.45e20);
59  TH1* mmbar = extrap->AntiNumuSurvComponent().TrueEnergy().ToTH1(3.45e20);
60  TH1* mtbar = extrap->AntiTauFromMuComponent().TrueEnergy().ToTH1(3.45e20);
61  TH1* etbar = extrap->AntiTauFromEComponent().TrueEnergy().ToTH1(3.45e20);
62 
63  TH1* nc = extrap->NCTotalComponent().ToTH1(3.45e20);
64 
65  double nNC = 0;
66  for(int i = 1; i <= nc->GetNbinsX(); ++i) nNC += nc->GetBinContent(i);
67 
68  const double nCosmic = ((TH1*)f->Get(("fdCosmic"+pid+"Scale/hist").c_str()))->Integral(0, -1);
69 
70  const double nData = ((pid == "LID") ? 6 : 11);
71 
72  // Numbers extracted from the paper
73  const double bkgSyst = (pid == "LEM") ? 13.4 : 10.8;
74  const double sigSyst = (pid == "LEM") ? 15.0 : 17.6;
75 
76  std::cout << "*** NOvA nue analysis, Jan 2016, " << pid << " ***" << std::endl << std::endl;
77  std::cout << "The data event count is " << nData << std::endl;
78  std::cout << "The NC background is " << nNC << " and the cosmic background is " << nCosmic << std::endl;
79  std::cout << "These are unaffected by oscillations." << std::endl;
80  std::cout << "The true energy spectra for all other components are given below." << std::endl;
81  std::cout << "For each column, alpha->beta, P(a->b) is set to 1. ie, you should multiply by P(a->b) to get the final prediction." << std::endl;
82  std::cout << "'b' is short for 'bar', ie antineutrinos." << std::endl;
83  std::cout << "The trueE value given is the bin center, at which P should be computed. The bin edges equally spaced in 1/trueE." << std::endl;
84  std::cout << "Everything is already scaled to the full detector 2.74e20 POT equivalent used in the analysis." << std::endl;
85  std::cout << "For systematics: all background components have a correlated gaussian " << bkgSyst << "% error." << std::endl;
86  std::cout << "The mu->e signal channel has an uncorrelated " << sigSyst << "% error" << std::endl;
88  std::cout << "The following expected counts are provided to allow you to test your oscillation calculations:" << std::endl << std::endl;
89 
91  calc->SetTh13(0);
92  PrintCount(pred.get(), nCosmic, calc);
94  calc->SetTh13(asin(sqrt(.1))/2);
95  calc->SetdCP(3*M_PI/2);
96  calc->SetDmsq32(+3e-3);
97  calc->SetTh23(asin(sqrt(.6)));
98  PrintCount(pred.get(), nCosmic, calc);
100  calc->SetTh13(asin(sqrt(.086))/2);
101  calc->SetdCP(M_PI/2);
102  calc->SetDmsq32(-2e-3);
103  calc->SetTh23(asin(sqrt(.4)));
104  PrintCount(pred.get(), nCosmic, calc);
105 
106  // Cross-checking with official background table in the paper
107  /*
108  calc->SetTh23(asin(sqrt(.5)));
109  calc->SetDmsq32(+2.37e-3);
110  calc->SetTh12(asin(sqrt(.846))/2);
111  calc->SetDmsq21(+7.53e-5);
112  calc->SetTh13(asin(sqrt(.086))/2);
113  calc->SetdCP(0);
114  PrintCount(pred.get(), nCosmic, calc);
115  */
116 
117  std::cout << std::endl;
118  std::cout << " trueE (GeV) e->e e->mu mu->e mu->mu mu->tau e->tau eb->eb eb->mub mub->eb mub->mub mub->taub eb->taub" << std::endl;
119 
120  for(int i = 1; i <= ee->GetNbinsX(); ++i){
121  if(ee->GetBinCenter(i) < 1) continue; // only zeros below here
122  std::cout << setw(12) << ee->GetBinCenter(i) << " \t";
123  std::cout << setw(12) << ee->GetBinContent(i) << setw(12) << em->GetBinContent(i) << setw(12) << me->GetBinContent(i) << setw(12) << mm->GetBinContent(i) << setw(12) << mt->GetBinContent(i) << setw(12) << et->GetBinContent(i) << setw(12);
124  std::cout << setw(12) << eebar->GetBinContent(i) << setw(12) << embar->GetBinContent(i) << setw(12) << mebar->GetBinContent(i) << setw(12) << mmbar->GetBinContent(i) << setw(12) << mtbar->GetBinContent(i) << setw(12) << etbar->GetBinContent(i) << std::endl;
125  }
126 }
virtual T GetDmsq21() const
osc::OscCalculatorDumb calc
Oscillation analysis framework, runs over CAF files outside of ART.
osc::IOscCalculatorAdjustable * DefaultOscCalc()
Create a new calculator with default assumptions for all parameters.
Definition: Calcs.cxx:49
virtual double GetRho() const
T sqrt(T number)
Definition: d0nt_math.hpp:156
void PrintCount(IPrediction *pred, double nCosmic, osc::IOscCalculatorAdjustable *calc)
void make_theorist_text_file(std::string predfile, std::string pid)
virtual void SetdCP(const T &dCP)=0
double Integral(double exposure, double *err=0, EExposureType expotype=kPOT) const
Return total number of events scaled to pot.
Definition: Spectrum.cxx:742
virtual double GetL() const
#define M_PI
Definition: SbMath.h:34
virtual Spectrum Predict(osc::IOscCalculator *calc) const =0
General interface to any calculator that lets you set the parameters.
Definition: NueSkimmer.h:24
virtual void SetTh23(const T &th23)=0
virtual void SetDmsq32(const T &dmsq32)=0
OStream cout
Definition: OStream.cxx:6
static constexpr Double_t mm
Definition: Munits.h:136
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
virtual T GetDmsq32() const
Standard interface to all prediction techniques.
Definition: IPrediction.h:57
virtual void SetTh13(const T &th13)=0
Float_t e
Definition: plot.C:35
T asin(T number)
Definition: d0nt_math.hpp:60