Functions | Variables
getContProf.C File Reference
#include <fstream>
#include <iostream>
#include "includes.h"
#include "varsandcuts.h"

Go to the source code of this file.

Functions

void getContProf (int sampleCut=9, std::string calculator="3a", bool usesysts=true, bool cosmics=true)
 

Variables

const int nquant = 4
 
const int nquantplus = 5
 

Function Documentation

void getContProf ( int  sampleCut = 9,
std::string  calculator = "3a",
bool  usesysts = true,
bool  cosmics = true 
)

Definition at line 28 of file getContProf.C.

References std::asin(), calc, Close(), make_associated_cosmic_defs::cosmics, om::cout, ana::DefaultOscCalc(), e, allTimeWatchdog::endl, submit_syst::fout, ana::getAllAna2017Systs(), ana::getNumuMarginalisedOscParam(), inDir, ana::Flavors::kAll, ana::kAna2017Livetime, ana::kAna2017POT, ana::Sign::kBoth, ana::Current::kCC, ana::kFitDmSq32Scaled, ana::kFitSinSqTheta23, ana::Current::kNC, ana::kNuMu, ana::kNumuAna2017, ana::kSolarConstraintsPDG2017, ana::kStandardSpillCuts, livetime, nquant, outDir, pot, ana::Profile(), runNovaSAM::release, ana::FrequentistSurface::SaveTo(), osc::_IOscCalcAdjustable< T >::SetDmsq32(), osc::_IOscCalcAdjustable< T >::SetTh23(), std::sqrt(), string, systs, ana::Spectrum::ToTH1(), and ana::WorldReactorConstraint2017().

28  {
29 
30 
31 
32  auto systs = getAllAna2017Systs(kNumuAna2017);
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 
49  osc::IOscCalcAdjustable* calc = DefaultOscCalc();
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
ana
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
systs
std::vector< SystGroupDef > systs
Definition: syst_header.h:385
nquant
const int nquant
Definition: getContProf.C:25
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:148
submit_syst.fout
fout
Definition: submit_syst.py:244
std::sqrt
T sqrt(T number)
Definition: d0nt_math.hpp:156
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:40
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
make_associated_cosmic_defs.cosmics
cosmics
Definition: make_associated_cosmic_defs.py:65
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:324
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::kStandardSpillCuts
const SpillCut kStandardSpillCuts
Apply this unless you&#39;re doing something special.
Definition: SpillCuts.h:49
ana::_Cut
Template for Cut and SpillCut.
Definition: Cut.h:15
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
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
ana::getNumuMarginalisedOscParam
std::vector< const IFitVar * > getNumuMarginalisedOscParam()
Definition: FitVarsNumu.cxx:9
string
enum BeamMode string

Variable Documentation

const int nquant = 4

Definition at line 25 of file getContProf.C.

Referenced by genie::flux::GNuMIFlux::AddFile(), estimate_pots(), getContProf(), and make_rockpred().

const int nquantplus = 5

Definition at line 26 of file getContProf.C.

Generated on Mon Jan 18 2021 14:15:58 for NOvA by   doxygen 1.8.11