predEventCountWithSystError.C
Go to the documentation of this file.
1 
3 
4 using namespace ana;
5 
6 void predEventCountWithSystError(bool useRHC = false, bool oscillations = false)
7 {
8 
9  osc::IOscCalcAdjustable* calcNoOsc = DefaultOscCalc(); // full systs normal h
10  float theta_real = 0.468;
11  float msq = 2.44;
12  kFitSinSqTheta23.SetValue(calcNoOsc,theta_real);
13  kFitDmSq32Scaled.SetValue(calcNoOsc,msq);
14  kFitDeltaInPiUnits.SetValue(calcNoOsc,1.5);
15  if(!oscillations) calcNoOsc -> SetL(0);
16 
17  std::vector<const ISyst*> totSysts = getAllAna2018Systs(kNumuAna2018);
18  // std::vector<const ISyst*> totSysts;
19  // totSysts.push_back(&kMuEScaleSyst2017);
20  // totSysts.push_back(&kRelMuEScaleSyst2017);
21 
22  ENu2018ExtrapType extrap = kNuMu;
23  const std::string predDir = "/nova/ana/nu_mu_ana/Ana2018/Predictions/";
24  const std::string predFHCFileName = predDir + "pred_fakeNDData_numuconcat_fhc__numu2018.root";
25  const std::string predRHCFileName = predDir + "pred_fakeNDData_numuconcat_rhc__numu2018.root";
26 
27  std::map< std::string, PredictionSystJoint2018* > mapPredsFHC;
28  std::map< std::string, PredictionSystJoint2018* > mapPredsRHC;
29  mapPredsFHC["FHC_Q1"] = new PredictionSystJoint2018(extrap, calcNoOsc, "fhc", 1, totSysts, predFHCFileName);
30  mapPredsFHC["FHC_Q2"] = new PredictionSystJoint2018(extrap, calcNoOsc, "fhc", 2, totSysts, predFHCFileName);
31  mapPredsFHC["FHC_Q3"] = new PredictionSystJoint2018(extrap, calcNoOsc, "fhc", 3, totSysts, predFHCFileName);
32  mapPredsFHC["FHC_Q4"] = new PredictionSystJoint2018(extrap, calcNoOsc, "fhc", 4, totSysts, predFHCFileName);
33  mapPredsRHC["RHC_Q1"] = new PredictionSystJoint2018(extrap, calcNoOsc, "rhc", 1, totSysts, predRHCFileName);
34  mapPredsRHC["RHC_Q2"] = new PredictionSystJoint2018(extrap, calcNoOsc, "rhc", 2, totSysts, predRHCFileName);
35  mapPredsRHC["RHC_Q3"] = new PredictionSystJoint2018(extrap, calcNoOsc, "rhc", 3, totSysts, predRHCFileName);
36  mapPredsRHC["RHC_Q4"] = new PredictionSystJoint2018(extrap, calcNoOsc, "rhc", 4, totSysts, predRHCFileName);
37 
38  double nomCount=0;
39  if(!useRHC)
40  for(auto &thisPred : mapPredsFHC)
41  nomCount += thisPred.second->Predict(calcNoOsc).ToTH1(kAna2018FHCPOT) -> Integral();
42  if(useRHC)
43  for(auto &thisPred : mapPredsRHC)
44  nomCount += thisPred.second->Predict(calcNoOsc).ToTH1(kAna2018RHCPOT) -> Integral();
45 
46  std::map< const ISyst*, double > upShifts;
47  std::map< const ISyst*, double > downShifts;
48  std::vector< double > upSystCount;
49  std::vector< double > dnSystCount;
50 
51  for(const ISyst* syst: totSysts) {
52  std::cout<< syst->ShortName();
53  double downErr = 10000;
54  double upErr = 0;
55 
56  for(double sigma=-1; sigma<=1; sigma+=0.2){
57  SystShifts* shift = new SystShifts(syst, sigma);
58  double eventsSyst=0;
59  if(!useRHC)
60  for(auto &thisPred : mapPredsFHC)
61  eventsSyst += thisPred.second->PredictSyst(calcNoOsc,*shift).ToTH1(kAna2018FHCPOT)->Integral();
62  if(useRHC)
63  for(auto &thisPred : mapPredsRHC)
64  eventsSyst += thisPred.second->PredictSyst(calcNoOsc,*shift).ToTH1(kAna2018RHCPOT)->Integral();
65 
66  std::cout<< " & " << eventsSyst - nomCount;
67 
68  if((eventsSyst<nomCount) && (eventsSyst<downErr)) downErr = eventsSyst;
69  if((eventsSyst>nomCount) && (eventsSyst>upErr)) upErr = eventsSyst;
70 
71  }
72  std::cout<< " \\\\" <<std::endl;
73  if(upErr>0){
74  upSystCount.push_back(upErr);
75  std::cout<< "+" << upErr-nomCount <<std::endl;
76  }
77  if(downErr<1000){
78  dnSystCount.push_back(downErr);
79  std::cout<< "-" << nomCount-downErr <<std::endl;
80  }
81 
82  }
83 
84  // add systs in quadrature
85  double totUpSystCount = 0;
86  double totDnSystCount = 0;
87  for( auto & systCount : upSystCount)
88  totUpSystCount += pow(systCount-nomCount,2);
89  for( auto & systCount : dnSystCount)
90  totDnSystCount += pow(nomCount-systCount,2);
91  double sqrtTotUpSystCount = sqrt(totUpSystCount);
92  double sqrtTotDnSystCount = sqrt(totDnSystCount);
93 
94  //////////////////////////////
95  // add cosmic background:
96  DontAddDirectory guard;
97  TFile* fCosFHC = new TFile("Files/cosmics_fhc__numu2018.root");
98  TFile* fCosRHC = new TFile("Files/cosmics_rhc__numu2018.root");
99  TH1D* hCosFHC_Q1 = (TH1D*) fCosFHC->Get("cosmics_q1");
100  TH1D* hCosFHC_Q2 = (TH1D*) fCosFHC->Get("cosmics_q2");
101  TH1D* hCosFHC_Q3 = (TH1D*) fCosFHC->Get("cosmics_q3");
102  TH1D* hCosFHC_Q4 = (TH1D*) fCosFHC->Get("cosmics_q4");
103  TH1D* hCosRHC_Q1 = (TH1D*) fCosRHC->Get("cosmics_q1");
104  TH1D* hCosRHC_Q2 = (TH1D*) fCosRHC->Get("cosmics_q2");
105  TH1D* hCosRHC_Q3 = (TH1D*) fCosRHC->Get("cosmics_q3");
106  TH1D* hCosRHC_Q4 = (TH1D*) fCosRHC->Get("cosmics_q4");
107  fCosFHC->Close();
108  fCosRHC->Close();
109  delete fCosFHC;
110  delete fCosRHC;
111  Spectrum specCosFHC_Q1(hCosFHC_Q1, kAna2018FHCPOT, kAna2018FHCLivetime);
112  Spectrum specCosFHC_Q2(hCosFHC_Q2, kAna2018FHCPOT, kAna2018FHCLivetime);
113  Spectrum specCosFHC_Q3(hCosFHC_Q3, kAna2018FHCPOT, kAna2018FHCLivetime);
114  Spectrum specCosFHC_Q4(hCosFHC_Q4, kAna2018FHCPOT, kAna2018FHCLivetime);
115  Spectrum specCosRHC_Q1(hCosRHC_Q1, kAna2018RHCPOT, kAna2018RHCLivetime);
116  Spectrum specCosRHC_Q2(hCosRHC_Q2, kAna2018RHCPOT, kAna2018RHCLivetime);
117  Spectrum specCosRHC_Q3(hCosRHC_Q3, kAna2018RHCPOT, kAna2018RHCLivetime);
118  Spectrum specCosRHC_Q4(hCosRHC_Q4, kAna2018RHCPOT, kAna2018RHCLivetime);
119  double cosCount = 0;
120  if(!useRHC){
121  cosCount += specCosFHC_Q1.Integral(kAna2018FHCPOT);
122  cosCount += specCosFHC_Q2.Integral(kAna2018FHCPOT);
123  cosCount += specCosFHC_Q3.Integral(kAna2018FHCPOT);
124  cosCount += specCosFHC_Q4.Integral(kAna2018FHCPOT);
125  }
126  if(useRHC){
127  cosCount += specCosRHC_Q1.Integral(kAna2018RHCPOT);
128  cosCount += specCosRHC_Q2.Integral(kAna2018RHCPOT);
129  cosCount += specCosRHC_Q3.Integral(kAna2018RHCPOT);
130  cosCount += specCosRHC_Q4.Integral(kAna2018RHCPOT);
131  }
132 
133 
134  double totCount = nomCount + cosCount;
135  double errup0 = pow(totCount,0.5);
136  double errdn0 = pow(totCount,0.5);
137 
138  std::cout<< std::fixed << std::setprecision(2)
139  << "\n\n\nEvent count: " << totCount << " +"
140  << sqrtTotUpSystCount << "/-" << sqrtTotDnSystCount << "(syst.) +"
141  << errup0 << "/-" << errdn0 << "(stat.)"
142  <<std::endl;
143 
144 }
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
Loads shifted spectra from files.
double Integral(const Spectrum &s, const double pot, int cvnregion)
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
T sqrt(T number)
Definition: d0nt_math.hpp:156
constexpr T pow(T x)
Definition: pow.h:75
void SetValue(osc::IOscCalcAdjustable *osc, double val) const override
Forward to wrapped Var&#39;s SetValue()
double Integral(double exposure, double *err=0, EExposureType expotype=kPOT) const
Return total number of events scaled to pot.
Definition: Spectrum.cxx:249
osc::IOscCalcAdjustable * DefaultOscCalc()
Create a new calculator with default assumptions for all parameters.
Definition: Calcs.cxx:49
Encapsulate code to systematically shift a caf::SRProxy.
Definition: ISyst.h:14
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
const double kAna2018RHCPOT
Definition: Exposures.h:208
std::vector< const ISyst * > getAllAna2018Systs(const EAnaType2018 ana, const bool smallgenie, const BeamType2018 beam, bool isFit)
void SetValue(osc::IOscCalcAdjustable *osc, double val) const override
const ConstrainedFitVarWithPrior fitDmSq32Scaled_UniformPrior & kFitDmSq32Scaled
double sigma(TH1F *hist, double percentile)
OStream cout
Definition: OStream.cxx:6
const ConstrainedFitVarWithPrior fitSsqTh23_UniformPriorSsqTh23 & kFitSinSqTheta23
void predEventCountWithSystError(bool useRHC=false, bool oscillations=false)
const double kAna2018FHCPOT
Definition: Exposures.h:207
Prevent histograms being added to the current directory.
Definition: UtilsExt.h:46
const FitVarWithPrior fitDeltaInPiUnits_UniformPriordCP & kFitDeltaInPiUnits
const double kAna2018FHCLivetime
Definition: Exposures.h:213
const double kAna2018RHCLivetime
Definition: Exposures.h:214
enum BeamMode string