joint_fit_future_tmp_loader_tools.h
Go to the documentation of this file.
1 #pragma once
13 #include "CAFAna/Systs/XSecSysts.h"
14 #include "CAFAna/Vars/FitVars.h"
15 #include "OscLib/IOscCalc.h"
16 #include "CAFAna/Analysis/Style.h"
19 
21 
22 
23 #include "TCanvas.h"
24 #include "TBox.h"
25 #include "TLatex.h"
26 #include "TColor.h"
27 #include "TGraph.h"
28 #include "TVectorD.h"
29 #include "TF1.h"
30 #include "TLegend.h"
31 #include "TLine.h"
32 #include "TMarker.h"
33 #include "TStyle.h"
34 #include "TSystem.h"
35 #include "TGaxis.h"
36 #include "TFile.h"
37 
38 using namespace ana;
39 
40 ////////////////////////////////////////////////////////////////////////////////////////////////
43  bool corrSysts,
45  bool isFakeData,
46  double futureExp,
47  std::string cvntype = "looseptp",
48  bool GetFromUPS=false,
49  bool minimizeMemory = false,
50  bool NERSC = false){
51 
52  cout<<"\n\n ============= You're loading Nue predictions ============="<<endl;
53 
54  std::cout<<"These arguments will be not used : "<<std::endl;
55  std::cout<<"decomp (always like in 2019): "<<decomp<<std::endl;
56  std::cout<<"corrSysts (always false): "<<corrSysts<<std::endl;
57  std::cout<<"isFakeData (always ND true data): "<<isFakeData<<std::endl;
58  std::cout<<"GetFromUPS, minimizeMemory, NERSC (always false) "<<std::endl;
59 
60  std::string path = "/cvmfs/nova.osgstorage.org/analysis/nue/Ana2020/tmp/";
61 
62  std::string name = "";
63  if (decomp.Contains("prop")) name = "extrap_prop";
64  if (decomp.Contains("combo")) name = "extrap_combo";
65 
66  IPrediction* extr = LoadFromFile<IPrediction>(path+"fd_regular_pred_"+beam+"_"+cvntype+"_v4.root", name).release();
67  cout<<"extr "<<extr->Predict(calc).Integral(futureExp)<<endl;
68  PredictionNoExtrap* noextr = LoadFromFile<PredictionNoExtrap>(path+"fd_regular_pred_"+beam+"_"+cvntype+"_v4.root", "noextrap_for_extand").release();
69  cout<<"noextr "<<noextr->Predict(calc).Integral(futureExp)<<endl;
70  IPrediction* pred_fid = new PredictionExtendToPeripheral(extr, noextr);
71 
72  SetFakeCalc(calc, 0.56, 2.48e-3, 0);
73  std::cout << "Scale to "<<futureExp<<std::endl;
74  std::cout << "Nue prediction (no rock)"<< pred_fid->Predict(calc).Integral(futureExp) <<std::endl;
75 
76  // Load rock
77  double nonsScale = 1./10.45;
78  double swapScale = 1./12.91;
79  std::string dirName = "/nova/ana/nu_e_ana/Ana2018/Predictions/FDRock/";
80  if(NERSC) dirName = std::string(getenv("FCHELPERANA2018_LIB_PATH"))+"/Predictions/FDRock/";
81  if(GetFromUPS) dirName = "/cvmfs/nova.osgstorage.org/analysis/nue/Ana2018/";
82  std::string fName2 = "fdrock_nue2018.root";
83  std::string rName;
84 
85  if(beam=="rhc") rName = "pred_FDRock_merg_rhc";
86  else rName = "pred_FDRock_merg_fhc";
87  if( decomp == "noextrap"){
88  if(beam=="rhc") rName = "pred_FDRock_rhc";
89  else rName = "pred_FDRock_fhc";
90  }
91 
92  auto rock = LoadFromFile<PredictionNoExtrap>(dirName + fName2, rName).release();
93 
94  std::cout << "Nue prediction (rock only - not scaled) "<< rock->Predict(calc).Integral(futureExp) <<std::endl;
95  std::cout << "Adding Rock" << endl;
96  auto nue_pred = new PredictionAddRock (pred_fid, rock, nonsScale, swapScale);
97  std::cout << "Nue prediction (rock added) "<< nue_pred->Predict(calc).Integral(futureExp) << std::endl;
98  std::cout << "Pure rock events: "<< nue_pred->Predict(calc).Integral(futureExp) - pred_fid->Predict(calc).Integral(futureExp) << std::endl;
99 
100  return pred_fid;
101 }
102 
103 ////////////////////////////////////////////////////////////////////////////////////////////////
104 std::pair <TH1D*, double > GetNueCosmicsFuture (std::string beam, double futureLive, bool GetFromUPS=false, bool NERSC=false){
105 
106  cout<<"\n\n ============= You're loading Nue cosmics ============="<<endl;
107 
108  cout<<"\n\nCall 2019 cosmic predictions "<<endl;
109  Spectrum* outtime;
111  std::string dir = "/cvmfs/nova.osgstorage.org/analysis/nue/Ana2020/tmp/";
112  /*if(beam=="rhc") name = "cosmic_spect_rhc_2019";
113  else name="cosmic_spect_fhc_2018";*/
114  name = "nue_cos";
115  outtime = LoadFromFile<Spectrum>(dir+"fd_cosmic_pred_"+beam+"_v4.root", name ).release();
116 
117  double outN = outtime->Integral(futureLive, 0, kLivetime);
118 
119  std::cout << "\nAdding " << outtime->Integral(outtime->Livetime(), 0, kLivetime)
120  << " cosmics from out-of-time." << std::endl;
121  std::cout << "Which scale to " << outtime->Integral(futureLive, 0, kLivetime)
122  << " in the spill window." << std::endl;
123  std::cout << " livetime " << futureLive <<" internal livetime "<<outtime->Livetime() <<std::endl;
124  std::cout << "uncertainty " << 1/sqrt(outN) << std::endl;
125 
126  TH1 *h = outtime->ToTH1(futureLive, kBlack, kSolid, kLivetime);
127 
128  return {outtime->ToTH1(futureLive, kBlack, kSolid, kLivetime), 1/sqrt(outN)};
129 
130  outtime->Clear();
131 }
132 
133 
134 ////////////////////////////////////////////////////////////////////////////////////////////////
135 std::vector <const IPrediction*> GetNumuPredictionsFuture(const int nq = 4,
136  bool useSysts = false,
137  std::string Beam = "fhc",
138  std::string OptDip = "OptFOMDip",
139  std::string SepQuant = "CombQuant",
140  bool GetFromUPS = false,
141  ENu2018ExtrapType numuExtrap = kNuMu,
142  bool minimizeMemory = false,
143  bool NERSC = false)
144 {
145 
146  std::string dir = "/nova/ana/users/karlwarb/Analysis2020/Sensitivity/200210/";
147  std::string fname = dir+"fd_numu_preds_extrap_"+Beam+"_"+OptDip+"_"+SepQuant+".root";
148 
149  std::cout<<"\n\n ============= You're loading NuMu Predictions ============="
150  << "\n These arguments will be used : "
151  << "\n\t useSysts : " << useSysts
152  << "\n\t Beam : " << Beam
153  << "\n\t OptDip : " << OptDip
154  << "\n\t SepQuant : " << SepQuant
155  << "\n\t GetFromUPS : " << GetFromUPS
156  << "\n\t numuExtrap : " << numuExtrap
157  << "\n\t minimizeMemory: " << minimizeMemory
158  << "\n\t NERSC : " << NERSC
159  << "\n Loading from: " << fname
160  <<std::endl;
161 
162  auto calc = DefaultOscCalc();
163 
164  std::vector <const IPrediction*> numu_preds;
165 
166  for (int quant = 1; quant <= nq; ++quant){
167  numu_preds.push_back(LoadFromFile<IPrediction>(fname, "extrap_q"+std::to_string(quant)).release());
168  std::cout << "Loaded quant " << quant << ", Integral is " << numu_preds.back()->Predict(calc).Integral(14e20) << std::endl;
169  }
170 
171  return numu_preds;
172 }
173 
174 ////////////////////////////////////////////////////////////////////////////////////////////////
175 std::vector <std::pair <TH1D*, double> > GetNumuCosmicsFuture(const int nq = 4,
176  std::string Beam = "fhc",
177  std::string OptDip = "OptFOMDip",
178  std::string SepQuant = "CombQuant",
179  double futureLive = 1600,
180  bool GetFromUPS = false,
181  bool NERSC = false)
182 {
183 
184  std::string dir = "/nova/ana/users/karlwarb/Analysis2020/Sensitivity/200210/";
185  std::string fname = dir+"fd_numu_cosmics_"+Beam+"_"+OptDip+"_"+SepQuant+".root";
186 
187  std::cout<<"\n\n ============= You're loading Numu Cosmics ============="
188  << "\n These arguments will be used : "
189  << "\n\t Beam : " << Beam
190  << "\n\t OptDip : " << OptDip
191  << "\n\t SepQuant : " << SepQuant
192  << "\n\t futureLive : " << futureLive
193  << "\n\t GetFromUPS : " << GetFromUPS
194  << "\n\t NERSC : " << NERSC
195  << "\n Loading from: " << fname
196  <<std::endl;
197 
198  TFile * fcosm = TFile::Open( (fname).c_str() );
199 
200  if(fcosm->IsZombie()) {std::cerr<< "bad cosmics\n"; exit(1);}
201  std::vector <std::pair <TH1D*, double > > numu_cosmics;
202 
203  for (int quant = 1; quant <=nq; ++quant) {
204  Spectrum* spec = Spectrum::LoadFrom(fcosm, ("numu_cos_q"+std::to_string(quant)).c_str()).release();;
205  TH1D *h = spec->ToTH1(futureLive, kBlack, kSolid, kLivetime);
206  numu_cosmics.push_back({h,1/sqrt(h->Integral())});
207  std::cout << "Loaded quantile " << quant << " " << h->Integral() << " scale error " << 1/sqrt(h->Integral()) << std::endl;
208  }//end quantiles
209 
210  return numu_cosmics;
211 }
212 
213 ////////////////////////////////////////////////////////////////////////////////////////////////
215  const IPrediction* pred,
216  TH1D* cosmics,
217  double futurePOT,
218  bool randomSystShifts = false,
219  bool asimov=false)
220 {
221  auto hMC = pred->Predict(calc).ToTH1(futurePOT);
222  if (randomSystShifts) {};
223  auto total = Spectrum(hMC, futurePOT,0); //Spectrum(cosmics, futurePOT, 0);
224  //total+= Spectrum(hMC, futurePOT,0);
225  auto ret = total.MockData(futurePOT);
226  if(asimov) ret = total.FakeData(futurePOT);
227 
228  //cout<<"data POT is "<<ret.POT();
229 
230  return ret;
231 }
232 
233 ////////////////////////////////////////////////////////////////////////////////////////////////
234 std::vector <Spectrum * > GetMockData(TString ana_options){
235 
236  cout<<"\n\n ============= You're loading Mock data from the file ============="<<endl;
237  std::vector <Spectrum *> spec;
238  std::string Dir = "./";
239  bool RHCmode = ana_options.Contains("rhc");
240  bool FHCmode = ana_options.Contains("fhc");
241  bool NueData = ana_options.Contains("nue");
242  bool NumuData = ana_options.Contains("numu");
243  if (NueData){
244  if (FHCmode) spec.push_back(LoadFromFile<Spectrum>(Dir+"data_for_the_fit.root", "nue_fhc").release());
245  if (RHCmode) spec.push_back(LoadFromFile<Spectrum>(Dir+"data_for_the_fit.root", "nue_rhc").release());
246  }
247  if (NumuData){
248  if (FHCmode){
249  for (int i = 1; i < 5; ++i) {
250  auto temp = LoadFromFile<Spectrum>(Dir+"data_for_the_fit.root", "numu_fhc_"+std::to_string(i)).release();
251  spec.push_back(temp);
252  }
253  }
254  if (RHCmode){
255  for (int i = 1; i < 5; ++i) {
256  auto temp = LoadFromFile<Spectrum>(Dir+"data_for_the_fit.root", "numu_rhc_"+std::to_string(i)).release();
257  spec.push_back(temp);
258  }
259  }
260 
261  }
262  return spec;
263 
264 }
265 
266 ////////////////////////////////////////////////////////////////////////////////////////////////
const XML_Char * name
Definition: expat.h:151
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
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
virtual Spectrum Predict(osc::IOscCalc *calc) const =0
T sqrt(T number)
Definition: d0nt_math.hpp:156
General interface to oscillation calculators.
Definition: StanTypedefs.h:22
OStream cerr
Definition: OStream.cxx:7
void Clear()
Definition: Spectrum.cxx:361
double Integral(double exposure, double *err=0, EExposureType expotype=kPOT) const
Return total number of events scaled to pot.
Definition: Spectrum.cxx:249
osc::OscCalcDumb calc
osc::IOscCalcAdjustable * DefaultOscCalc()
Create a new calculator with default assumptions for all parameters.
Definition: Calcs.cxx:49
std::pair< TH1D *, double > GetNueCosmicsFuture(std::string beam, double futureLive, bool GetFromUPS=false, bool NERSC=false)
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
std::vector< Spectrum * > GetMockData(TString ana_options)
std::vector< std::pair< TH1D *, double > > GetNumuCosmicsFuture(const int nq=4, std::string Beam="fhc", std::string OptDip="OptFOMDip", std::string SepQuant="CombQuant", double futureLive=1600, bool GetFromUPS=false, bool NERSC=false)
static std::unique_ptr< Spectrum > LoadFrom(TDirectory *dir, const std::string &name)
Definition: Spectrum.cxx:535
Spectrum Predict(osc::IOscCalc *calc) const override
std::string getenv(std::string const &name)
std::vector< const IPrediction * > GetNumuPredictionsFuture(const int nq=4, bool useSysts=false, std::string Beam="fhc", std::string OptDip="OptFOMDip", std::string SepQuant="CombQuant", bool GetFromUPS=false, ENu2018ExtrapType numuExtrap=kNuMu, bool minimizeMemory=false, bool NERSC=false)
std::vector< float > Spectrum
Definition: Constants.h:610
void SetFakeCalc(osc::IOscCalcAdjustable *calc, double ssth23=-5, double dmsq32=-5, double dCP_Pi=-5)
Spectrum GenerateFutureData(osc::IOscCalc *calc, const IPrediction *pred, TH1D *cosmics, double futurePOT, bool randomSystShifts=false, bool asimov=false)
OStream cout
Definition: OStream.cxx:6
const std::string path
Definition: plot_BEN.C:43
std::string dirName
Definition: PlotSpectra.h:47
TDirectory * dir
Definition: macro.C:5
exit(0)
TCut outtime("(tave>0.0&&tave<200.0) || (tave>250.0&&tave<450.0)")
Standard interface to all prediction techniques.
Definition: IPrediction.h:57
Prediction that just uses FD MC, with no extrapolation.
void Beam(bool isRHC)
std::string to_string(ModuleType mt)
Definition: ModuleType.h:32
Float_t e
Definition: plot.C:35
double Livetime() const
Seconds. For informational purposes only. No calculations use this.
Definition: Spectrum.h:230
const IPrediction * GetNuePredictionFuture(TString decomp, osc::IOscCalcAdjustable *calc, bool corrSysts, std::string beam, bool isFakeData, double futureExp, std::string cvntype="looseptp", bool GetFromUPS=false, bool minimizeMemory=false, bool NERSC=false)
void rock(std::string suffix="full")
Definition: rock.C:28
enum BeamMode string