getContProf.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.C //
9 // Makes (sin2theta_23, deltam2_32) surface and profiles //
10 // Uses real 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(int sampleCut = 9, std::string calculator = "3a", bool usesysts = true, bool cosmics = true){
29 
30 
31 
33  std::string SystName = "";
34  if(usesysts == true){
35  std::cout << "With systematics" << std::endl;
36  SystName = "WithSysts";
37  }
38 
39  if(usesysts == false){
40  std::cout << "Stats only" << std::endl;
41  SystName = "StatsOnly";
42  systs = {};
43  }
44  for (auto & sys:systs) std::cout <<sys->ShortName() << std::endl;
45 
46 
47 
48 
50  std::string CalcName = "";
51  std::cout << "\nCalculator:" << std::endl;
52  if(calculator == "sa"){
53  std::cout << "SA best fit" << std::endl;
54  CalcName = "SABestFit";
55  calc->SetDmsq32(2.6746e-3); // SA value
56  calc->SetTh23(0.68696); // non max mixing (ssqth23 = 0.4022)
57  }
58  if(calculator == "3a"){
59  std::cout << "3A best fit" <<std::endl;
60  CalcName = "3ABestFit";
61  calc->SetTh23(asin(sqrt(0.559455))); //
62  calc->SetDmsq32(0.00244); // seed from 3A best fit
63  }
64 
65 
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 = kAna2017POT;
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(pot);
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 = "/nova/ana/nu_mu_ana/Ana2017/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 
164  // predictions new scheme
165  ENumu2017ExtrapType extrap = kNuMu;
166  std::cout << "\nDefining PredictionInterp" << std::endl;
167  std::vector<PredictionInterp*> predictionVec;
168  PredictionSystNumu2017 predQ1(extrap, calc, 1);
169  PredictionSystNumu2017 predQ2(extrap, calc, 2);
170  PredictionSystNumu2017 predQ3(extrap, calc, 3);
171  PredictionSystNumu2017 predQ4(extrap, calc, 4);
172  predictionVec.push_back(&predQ1);
173  predictionVec.push_back(&predQ2);
174  predictionVec.push_back(&predQ3);
175  predictionVec.push_back(&predQ4);
176 
177 
178  TH1* hBkg = predictionVec[0]->PredictComponent(calc, Flavors::kAll, Current::kNC, Sign::kBoth).ToTH1(kAna2017POT, kGray);
179  TH1* hPred = predictionVec[0]->PredictComponent(calc, Flavors::kAll, Current::kCC, Sign::kBoth).ToTH1(kAna2017POT, kGray);
180 
181  for(int epochId = 1; epochId < nquant; epochId++){
182  hBkg->Add(predictionVec[epochId]->PredictComponent(calc, Flavors::kAll, Current::kNC, Sign::kBoth).ToTH1(kAna2017POT, kGray));
183  hPred->Add(predictionVec[epochId]->PredictComponent(calc, Flavors::kAll, Current::kCC, Sign::kBoth).ToTH1(kAna2017POT, kGray));
184  }
185  hPred->Add(CosQuant0);
186  hPred->Add(hBkg);
187 
188 
189 
190  //experiment
191  SingleSampleExperiment realExpt_quant1(predictionVec[0], *FD_Quant1, cosmics1, 0.1);
192  SingleSampleExperiment realExpt_quant2(predictionVec[1], *FD_Quant2, cosmics2, 0.1);
193  SingleSampleExperiment realExpt_quant3(predictionVec[2], *FD_Quant3, cosmics3, 0.1);
194  SingleSampleExperiment realExpt_quant4(predictionVec[3], *FD_Quant4, cosmics4, 0.1);
195  predictionVec.clear();
196 
197 
198  std::vector <const IExperiment*> vExpts;
199  vExpts.push_back(&realExpt_quant1);
200  vExpts.push_back(&realExpt_quant2);
201  vExpts.push_back(&realExpt_quant3);
202  vExpts.push_back(&realExpt_quant4);
203  vExpts.push_back(&kSolarConstraintsPDG2017);
204  vExpts.push_back(WorldReactorConstraint2017());
205  MultiExperiment multExpt(vExpts);
206 
207  TFile* fout = new TFile(outDir + "ContProf_" + SystName + "_" + CalcName + "_3ACut_3AEnergy_OptBinning_" + SampleName + "_Cosmics.root","RECREATE");
208 
209  std::vector<const IFitVar*> oscpar = getNumuMarginalisedOscParam();
210  FrequentistSurface* surface = new FrequentistSurface( &multExpt, calc, &kFitSinSqTheta23, 40, 0.3, 0.7, &kFitDmSq32Scaled, 50, 2.0, 3.0, oscpar, systs);
211  //FrequentistSurface* surface = new FrequentistSurface( &multExpt, calc, &kFitSinSqTheta23, 40, 0.3, 0.7, &kFitDmSq32Scaled, 50, 2.0, 3.0, {&kFitDeltaInPiUnits,&kFitSinSq2Theta13}, systs);
212  TH1* hProfileTheta23 = Profile(&multExpt, calc, &kFitSinSqTheta23, 40, 0.3, 0.7, -1, {&kFitDmSq32Scaled}, systs);
213  TH1* hProfileDelta32 = Profile(&multExpt, calc, &kFitDmSq32Scaled, 50, 2.0, 3.0, -1, {&kFitSinSqTheta23}, systs);
214 
215  surface->SaveTo(fout, "surface");
216  hProfileTheta23->Write("profiletheta23");
217  hProfileDelta32->Write("profiledelta32");
218  hData->Write("Data");
219  hPred->Write("Pred");
220 
221  vExpts.clear();
222 
223 
224 } // End of function
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
std::vector< SystGroupDef > systs
Definition: syst_header.h:385
const int nquant
Definition: getContProf.C:25
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:148
T sqrt(T number)
Definition: d0nt_math.hpp:156
const double kAna2017POT
Definition: Exposures.h:174
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
osc::OscCalcDumb calc
osc::IOscCalcAdjustable * DefaultOscCalc()
Create a new calculator with default assumptions for all parameters.
Definition: Calcs.cxx:49
virtual void SetDmsq32(const T &dmsq32)=0
std::string outDir
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
Log-likelihood scan across two parameters.
TString inDir
Charged-current interactions.
Definition: IPrediction.h:39
const double kAna2017Livetime
Definition: Exposures.h:200
Loads shifted spectra from files.
#define pot
const ReactorExperiment * WorldReactorConstraint2017()
Reactor constraint from PDG 2017.
const ConstrainedFitVarWithPrior fitDmSq32Scaled_UniformPrior & kFitDmSq32Scaled
c1 Close()
OStream cout
Definition: OStream.cxx:6
Combine multiple component experiments.
void getContProf(int sampleCut=9, std::string calculator="3a", bool usesysts=true, bool cosmics=true)
Definition: getContProf.C:28
const int nquantplus
Definition: getContProf.C:26
const SolarConstraints kSolarConstraintsPDG2017(7.53e-5, 0.18e-5, 0.851, 0.020)
double livetime
Definition: saveFDMCHists.C:21
void SaveTo(TDirectory *dir, const std::string &name) const
const ConstrainedFitVarWithPrior fitSsqTh23_UniformPriorSsqTh23 & kFitSinSqTheta23
virtual void SetTh23(const T &th23)=0
Neutral-current interactions.
Definition: IPrediction.h:40
Both neutrinos and antineutrinos.
Definition: IPrediction.h:52
const SpillCut kStandardSpillCuts
Apply this unless you&#39;re doing something special.
Definition: SpillCuts.h:49
Template for Cut and SpillCut.
Definition: Cut.h:15
std::vector< const ISyst * > getAllAna2017Systs(const EAnaType2017 ana, const bool smallgenie)
All neutrinos, any flavor.
Definition: IPrediction.h:26
Float_t e
Definition: plot.C:35
T asin(T number)
Definition: d0nt_math.hpp:60
Compare a single data spectrum to the MC + cosmics expectation.
std::vector< const IFitVar * > getNumuMarginalisedOscParam()
Definition: FitVarsNumu.cxx:9
enum BeamMode string