getContProf_Sensitivity.C
Go to the documentation of this file.
1 /////////////////////////////////////////////////////////////////////////
2 // NOvA Muon Neutrino Appearance Analysis 2017 //
3 // //
4 // Author: Diana Patricia Mendez d.p.mendez@sussex.ac.uk //
5 // //
6 //////////////////////////////////////////////////////////////////////////
7 // //
8 // getContProf_Sensitivity.C //
9 // Makes (sin2theta_23, deltam2_32) sensitivity surface and profiles //
10 // Uses fake FD data and includes cosmics to prediction //
11 // Read pre-made predictions made with systematic file //
12 // but can be used with the old standard predictions //
13 // //
14 //////////////////////////////////////////////////////////////////////////
15 
16 #include <fstream>
17 #include <iostream>
18 
19 #include "includes.h"
20 #include "varsandcuts.h"
21 
22 using namespace ana;
23 
24 
25 const int nquant = 4;
26 const int nquantplus = 5; // includes quant0
27 
28 void getContProf_Sensitivity(int sampleCut = 9, std::string calculator = "3a", bool usesysts = true, bool cosmics = true){
29 
30 
31  auto systs = getAllAna2017Systs(kNumuAna2017);
32  std::string SystName = "";
33  if(usesysts == true){
34  std::cout << "With systematics" << std::endl;
35  SystName = "WithSysts";
36  }
37 
38  if(usesysts == false){
39  std::cout << "Stats only" << std::endl;
40  SystName = "StatsOnly";
41  systs = {};
42  }
43  for (auto & sys:systs) std::cout <<sys->ShortName() << std::endl;
44 
45 
46 
47  osc::IOscCalcAdjustable* calc = DefaultOscCalc();
48  std::string CalcName = "";
49  std::cout << "\nCalculator:" << std::endl;
50  if(calculator == "sa"){
51  std::cout << "SA best fit" << std::endl;
52  CalcName = "SABestFit";
53  calc->SetDmsq32(2.6746e-3); // SA value
54  calc->SetTh23(0.68696); // non max mixing (ssqth23 = 0.4022)
55  }
56  if(calculator == "3a"){
57  std::cout << "3A best fit" <<std::endl;
58  CalcName = "3ABestFit";
59  calc->SetTh23(asin(sqrt(0.559455))); //
60  calc->SetDmsq32(0.00244); // seed from 3A best fit
61  }
62 
63 
64 
65  // sample
66  std::string SampleName = "";
67  std::cout << "\nSample:" << std::endl;
68  SpillCut runSpillCut = (const_cast<SpillCut&>(kStandardSpillCuts));
69  double pot;
70  double livetime;
71  std::cout << "9e20 POT" << std::endl;
72  SampleName = "9e20";
73  pot = kPeriod1POT + kPeriod2POT + kPeriod3POT + kPeriod5POT;
74  livetime = kAna2017Livetime;
75 
76 
77  //inputs
78  std::string ana = "quantiles";
79  TString outDir = "";
80  std::string inDirData = "/nova/ana/nu_mu_ana/Ana2017/Data";
81  std::string inDir = "/nova/ana/nu_mu_ana/Ana2017/Predictions/";
82  std::string inName = "Prediction_StatsOnly_3ACut_3AEnergy_OptBinning_";
83 
84  TFile* fDataQ0 = new TFile( (inDirData + "/fd_data_3ACut_3AEnergy_OptBinning_full_" + SampleName + "_Quant0.root").c_str() );
85  TFile* fDataQ1 = new TFile( (inDirData + "/fd_data_3ACut_3AEnergy_OptBinning_full_" + SampleName + "_Quant1.root").c_str() );
86  TFile* fDataQ2 = new TFile( (inDirData + "/fd_data_3ACut_3AEnergy_OptBinning_full_" + SampleName + "_Quant2.root").c_str() );
87  TFile* fDataQ3 = new TFile( (inDirData + "/fd_data_3ACut_3AEnergy_OptBinning_full_" + SampleName + "_Quant3.root").c_str() );
88  TFile* fDataQ4 = new TFile( (inDirData + "/fd_data_3ACut_3AEnergy_OptBinning_full_" + SampleName + "_Quant4.root").c_str() );
89  Spectrum* FD_Quant0 = LoadFrom< Spectrum >(fDataQ0, "fd_data").release();
90  Spectrum* FD_Quant1 = LoadFrom< Spectrum >(fDataQ1, "fd_data").release();
91  Spectrum* FD_Quant2 = LoadFrom< Spectrum >(fDataQ2, "fd_data").release();
92  Spectrum* FD_Quant3 = LoadFrom< Spectrum >(fDataQ3, "fd_data").release();
93  Spectrum* FD_Quant4 = LoadFrom< Spectrum >(fDataQ4, "fd_data").release();
94  TH1* hData = FD_Quant0->ToTH1(kHereSAPOT);
95  double potData = pot;
96  double liveData = livetime;
97 
98 
99  fDataQ1 -> Close();
100  fDataQ2 -> Close();
101  fDataQ3 -> Close();
102  fDataQ4 -> Close();
103  delete fDataQ1;
104  delete fDataQ2;
105  delete fDataQ3;
106  delete fDataQ4;
107 
108 
109 
110 
111  //cosmics
112  std::string CosmicsSpecFile;
113  std::string CosmicsHistName;
114  std::string cosmicsName = "NoCosmics";
115  if(cosmics){
116  std::cout << "\nCosmics" << std::endl;
117  CosmicsSpecFile = "Cosmics/cosmicsScaled_NewBinning_" + SampleName + ".root";
118  CosmicsHistName = "cosmics3A";
119  cosmicsName = "WithCosmics";
120  }
121 
122 
123 
124  std::cout << "\nLoading cosmic distributions" << std::endl;
125  TFile inCosmicsFile(CosmicsSpecFile.c_str());
126  auto *CosQuant1 = (TH1*)inCosmicsFile.Get("q1all");
127  auto *CosQuant2 = (TH1*)inCosmicsFile.Get("q2all");
128  auto *CosQuant3 = (TH1*)inCosmicsFile.Get("q3all");
129  auto *CosQuant4 = (TH1*)inCosmicsFile.Get("q4all");
130  auto *CosQuant0 = (TH1*)CosQuant1->Clone("CosQuant0");
131  CosQuant0->Add(CosQuant2);
132  CosQuant0->Add(CosQuant3);
133  CosQuant0->Add(CosQuant4);
134  //auto hCosmics = (TH1*)inCosmicsFile.Get(CosmicsHistName.c_str());
135  Spectrum cosmics1(CosQuant1, pot, livetime);
136  Spectrum cosmics2(CosQuant2, pot, livetime);
137  Spectrum cosmics3(CosQuant3, pot, livetime);
138  Spectrum cosmics4(CosQuant4, pot, livetime);
139 
140 
141 
142 
143 
144  /*
145  //predictions
146  std::vector <string> epoch = {"full"};
147  const int nepochs = epoch.size();
148  std::cout << "\nNumber o periods:" << nepochs << std::endl;
149  std::cout << "\nDefining PredictionInterp" << std::endl;
150  std::vector<PredictionInterp*> predictionVec;
151  for(int quantId = 0; quantId < nquant; quantId++){
152  std::cout << "Quantile " << quantId + 1 << std::endl;
153  std::string input = inDir + inName + Form("Quant%d_", quantId+1) + "full_" + SampleName + ".root";
154  TFile* file = new TFile(input.c_str());
155  std::cout << "\nPrediction vectors" << std::endl;
156  predictionVec.push_back(LoadFrom<PredictionInterp>(file, "prediction")).release();
157  file->Close();
158  std::cout << "Read and closed input files" << std::endl;
159  }//loop quant
160  */
161 
162 
163  ENumu2017ExtrapType extrap = kNuMu;
164  std::cout << "\nDefining PredictionInterp" << std::endl;
165  std::vector<PredictionInterp*> predictionVec;
166  PredictionSystNumu2017 predQ1(extrap, calc, 1);
167  PredictionSystNumu2017 predQ2(extrap, calc, 2);
168  PredictionSystNumu2017 predQ3(extrap, calc, 3);
169  PredictionSystNumu2017 predQ4(extrap, calc, 4);
170  predictionVec.push_back(&predQ1);
171  predictionVec.push_back(&predQ2);
172  predictionVec.push_back(&predQ3);
173  predictionVec.push_back(&predQ4);
174  TH1* hBkg = predictionVec[0]->PredictComponent(calc, Flavors::kAll, Current::kNC, Sign::kBoth).ToTH1(kAna2017POT, kBlue);
175  TH1* hPred = predictionVec[0]->PredictComponent(calc, Flavors::kAll, Current::kCC, Sign::kBoth).ToTH1(kAna2017POT, kBlue);
176 
177  for(int epochId = 1; epochId < nquant; epochId++){
178  hBkg->Add(predictionVec[epochId]->PredictComponent(calc, Flavors::kAll, Current::kNC, Sign::kBoth).ToTH1(kAna2017POT, kBlue));
179  hPred->Add(predictionVec[epochId]->PredictComponent(calc, Flavors::kAll, Current::kCC, Sign::kBoth).ToTH1(kAna2017POT, kBlue));
180  }
181  hPred->Add(CosQuant0);
182  hPred->Add(hBkg);
183 
184 
185  // add cosmics to fake data
186  Spectrum fakePred_quant1 = predictionVec[0]->Predict(calc);
187  Spectrum fakePred_quant2 = predictionVec[1]->Predict(calc);
188  Spectrum fakePred_quant3 = predictionVec[2]->Predict(calc);
189  Spectrum fakePred_quant4 = predictionVec[3]->Predict(calc);
190 
191  Spectrum fakeData_quant1 = fakePred_quant1.FakeData(kAna2017POT);
192  Spectrum fakeData_quant2 = fakePred_quant2.FakeData(kAna2017POT);
193  Spectrum fakeData_quant3 = fakePred_quant3.FakeData(kAna2017POT);
194  Spectrum fakeData_quant4 = fakePred_quant4.FakeData(kAna2017POT);
195  fakeData_quant1 += cosmics1;
196  fakeData_quant2 += cosmics2;
197  fakeData_quant3 += cosmics3;
198  fakeData_quant4 += cosmics4;
199  SingleSampleExperiment fakeExpt_quant1(predictionVec[0], fakeData_quant1, cosmics1, 0.1);
200  SingleSampleExperiment fakeExpt_quant2(predictionVec[1], fakeData_quant2, cosmics2, 0.1);
201  SingleSampleExperiment fakeExpt_quant3(predictionVec[2], fakeData_quant3, cosmics3, 0.1);
202  SingleSampleExperiment fakeExpt_quant4(predictionVec[3], fakeData_quant4, cosmics4, 0.1);
203 
204  predictionVec.clear();
205 
206  std::vector <const IExperiment*> vExpts;
207  vExpts.push_back(&fakeExpt_quant1);
208  vExpts.push_back(&fakeExpt_quant2);
209  vExpts.push_back(&fakeExpt_quant3);
210  vExpts.push_back(&fakeExpt_quant4);
211  vExpts.push_back(&kSolarConstraintsPDG2017);
212  vExpts.push_back(WorldReactorConstraint2017());
213  MultiExperiment multExpt(vExpts);
214 
215  TFile* fout = new TFile(outDir + "ContProfSens_" + SystName + "_" + CalcName + "_3ACut_3AEnergy_OptBinning_" + SampleName + "_Cosmics.root","RECREATE");
216 
217  FrequentistSurface* surface = new FrequentistSurface( &multExpt, calc, &kFitSinSqTheta23, 40, 0.3, 0.7, &kFitDmSq32Scaled, 50, 2.0, 3.0, {&kFitDeltaInPiUnits,&kFitSinSq2Theta13}, systs);
218  TH1* hProfileTheta23 = Profile(&multExpt, calc, &kFitSinSqTheta23, 40, 0.3, 0.7, -1, {&kFitDmSq32Scaled}, systs);
219  TH1* hProfileDelta32 = Profile(&multExpt, calc, &kFitDmSq32Scaled, 50, 2.0, 3.0, -1, {&kFitSinSqTheta23}, systs);
220 
221  surface->SaveTo(fout, "surface");
222  hProfileTheta23->Write("profiletheta23");
223  hProfileDelta32->Write("profiledelta32");
224  hData->Write("Data");
225  hPred->Write("Pred");
226 
227  vExpts.clear();
228 
229 
230 } // End of function
nquant
const int nquant
Definition: getContProf_Sensitivity.C:25
ana
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
systs
std::vector< SystGroupDef > systs
Definition: syst_header.h:385
ana::Spectrum::ToTH1
TH1D * ToTH1(double exposure, Color_t col=kBlack, Style_t style=kSolid, EExposureType expotype=kPOT, EBinType bintype=kBinContent) const
Histogram made from this Spectrum, scaled to some exposure.
Definition: Spectrum.cxx:176
submit_syst.fout
fout
Definition: submit_syst.py:244
std::sqrt
T sqrt(T number)
Definition: d0nt_math.hpp:156
nquantplus
const int nquantplus
Definition: getContProf_Sensitivity.C:26
ana::kAna2017POT
const double kAna2017POT
Definition: Exposures.h:174
ana::Profile
TGraph * Profile(const IExperiment *expt, osc::IOscCalcAdjustable *calc, const IFitVar *v, int nbinsx, double minx, double maxx, double input_minchi, const std::vector< const IFitVar * > &profVars, const std::vector< const ISyst * > &profSysts, const SeedList &seedPts, const std::vector< SystShifts > &systSeedPts, std::map< const IFitVar *, TGraph * > &profVarsMap, std::map< const ISyst *, TGraph * > &profSystsMap, MinuitFitter::FitOpts opts)
scan in one variable, profiling over all others
Definition: Fit.cxx:48
calc
osc::OscCalcDumb calc
Definition: make_nueFDprediction_kinematics.C:70
ana::DefaultOscCalc
osc::IOscCalcAdjustable * DefaultOscCalc()
Create a new calculator with default assumptions for all parameters.
Definition: Calcs.cxx:49
osc::_IOscCalcAdjustable::SetDmsq32
virtual void SetDmsq32(const T &dmsq32)=0
outDir
std::string outDir
Definition: plot_ND_DataMC.C:35
ana::Spectrum
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:39
ana::FrequentistSurface
Log-likelihood scan across two parameters.
Definition: FrequentistSurface.h:27
inDir
TString inDir
Definition: plot_ND_DataMC.C:33
ana::Current::kCC
Charged-current interactions.
Definition: IPrediction.h:39
ana::kAna2017Livetime
const double kAna2017Livetime
Definition: Exposures.h:200
ana::Spectrum::FakeData
Spectrum FakeData(double pot) const
Synonymous with AsimovData(). Retained for compatibility.
Definition: Spectrum.cxx:377
make_associated_cosmic_defs.cosmics
cosmics
Definition: make_associated_cosmic_defs.py:68
ana::PredictionSystNumu2017
Loads shifted spectra from files.
Definition: PredictionSystNumu2017.h:115
pot
#define pot
Definition: template_basic.C:49
ana::WorldReactorConstraint2017
const ReactorExperiment * WorldReactorConstraint2017()
Reactor constraint from PDG 2017.
Definition: ReactorExperiment.cxx:115
ana::kFitDmSq32Scaled
const ConstrainedFitVarWithPrior fitDmSq32Scaled_UniformPrior & kFitDmSq32Scaled
Definition: VarsWithPriors.cxx:13
Close
c1 Close()
om::cout
OStream cout
Definition: OStream.cxx:6
ana::MultiExperiment
Combine multiple component experiments.
Definition: MultiExperiment.h:11
ana::kSolarConstraintsPDG2017
const SolarConstraints kSolarConstraintsPDG2017(7.53e-5, 0.18e-5, 0.851, 0.020)
Definition: SolarConstraints.h:56
livetime
double livetime
Definition: saveFDMCHists.C:21
ana::FrequentistSurface::SaveTo
void SaveTo(TDirectory *dir, const std::string &name) const
Definition: FrequentistSurface.cxx:328
ana::kFitSinSqTheta23
const ConstrainedFitVarWithPrior fitSsqTh23_UniformPriorSsqTh23 & kFitSinSqTheta23
Definition: VarsWithPriors.cxx:9
osc::_IOscCalcAdjustable::SetTh23
virtual void SetTh23(const T &th23)=0
osc::_IOscCalcAdjustable< double >
ana::Current::kNC
Neutral-current interactions.
Definition: IPrediction.h:40
ana::Sign::kBoth
Both neutrinos and antineutrinos.
Definition: IPrediction.h:52
ana::kStandardSpillCuts
const SpillCut kStandardSpillCuts
Apply this unless you&#39;re doing something special.
Definition: SpillCuts.h:49
ana::kFitSinSq2Theta13
const FitSinSq2Theta13 kFitSinSq2Theta13
Definition: FitVars.cxx:13
ana::_Cut
Template for Cuts applied to any type of object.
Definition: Cut.h:15
getContProf_Sensitivity
void getContProf_Sensitivity(int sampleCut=9, std::string calculator="3a", bool usesysts=true, bool cosmics=true)
Definition: getContProf_Sensitivity.C:28
ana::getAllAna2017Systs
std::vector< const ISyst * > getAllAna2017Systs(const EAnaType2017 ana, const bool smallgenie)
Definition: JointAna2017Systs.cxx:146
ana::Flavors::kAll
All neutrinos, any flavor.
Definition: IPrediction.h:26
e
Float_t e
Definition: plot.C:35
ana::kFitDeltaInPiUnits
const FitVarWithPrior fitDeltaInPiUnits_UniformPriordCP & kFitDeltaInPiUnits
Definition: VarsWithPriors.cxx:17
kBlue
enum BeamMode kBlue
std::asin
T asin(T number)
Definition: d0nt_math.hpp:60
ana::SingleSampleExperiment
Compare a single data spectrum to the MC + cosmics expectation.
Definition: SingleSampleExperiment.h:21
string
enum BeamMode string
Generated on Fri May 7 2021 14:08:19 for NOvA by   doxygen 1.8.11