makePredTables.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 // makePredTables.C //
9 // Makes tables with predicted number of events for oscillation chanels //
10 // Writes to screen. Save the log to pick up tables in latex format //
11 // //
12 // quantile 0 is always the sum of all 4 quantiles //
13 // //
14 //////////////////////////////////////////////////////////////////////////
15 
16 #include "includes.h"
17 #include "varsandcuts.h"
18 
19 using namespace ana;
20 
21 
22 
23 void makePredTables(std::string calculator = "maxmix"){
24 
25 
27  std::string CalcName = "";
28 
29  if(calculator == "3a"){
30  CalcName = "3A best fit";
31  calc->SetTh23(asin(sqrt(0.559455)));
32  calc->SetDmsq32(0.00244);
33  }
34  if(calculator == "sa"){
35  CalcName = "SA best fit";
36  calc->SetDmsq32(2.6746e-3);
37  calc->SetTh23(0.68696);
38  }
39  if(calculator == "maxmix"){
40  CalcName = "Maimal Mixing";
41  calc->SetTh23(asin(sqrt(0.5))); // sin2theta23 = 0.5
42  calc->SetDmsq32(0.00244); // 3a best fit
43  }
44 
45 
46  // sample
47  std::string SampleName = "9e20";
48  double pot = kAna2017POT;
49  double livetime = kAna2017Livetime;
50 
51 
52  //inputs
53  TString outDir = "";
54  std::string inDirData = "/nova/ana/nu_mu_ana/Ana2017/Data/";
55  std::string inDir = "/nova/ana/nu_mu_ana/Ana2017/Predictions/";
56  std::string inName = "Prediction_StatsOnly_3ACut_3AEnergy_OptBinning_";
57 
58 
59  //cosmics
60  std::string CosmicsSpecFile;
61  std::string CosmicsHistName;
62  std::string cosmicsName = "NoCosmics";
63  std::cout << "\nCosmics" << std::endl;
64  CosmicsSpecFile = "/nova/ana/nu_mu_ana/Ana2017/Cosmics/cosmicsScaled_NewBinning_" + SampleName + ".root";
65  CosmicsHistName = "cosmics3A";
66 
67  std::cout << "\nLoading cosmic distributions" << std::endl;
68  TFile inCosmicsFile(CosmicsSpecFile.c_str());
69 
70  std::vector<Spectrum> sCosmics;
71  std::vector<TH1*> hCosmics;
72  hCosmics.push_back((TH1*)inCosmicsFile.Get("q1all"));
73  hCosmics.push_back((TH1*)inCosmicsFile.Get("q1all"));
74  hCosmics.push_back((TH1*)inCosmicsFile.Get("q2all"));
75  hCosmics.push_back((TH1*)inCosmicsFile.Get("q3all"));
76  hCosmics.push_back((TH1*)inCosmicsFile.Get("q4all"));
77  for(int quantId = 2; quantId < nQuantPlus; quantId++){
78  hCosmics[0]->Add(hCosmics[quantId]);
79  }
80  for(int quantId = 0; quantId<nQuantPlus; quantId++){
81  Spectrum cosmics(hCosmics[quantId], pot, livetime);
82  sCosmics.push_back(cosmics);
83  }
84 
85 
86 
87 
88  std::cout << "\nDefining predictions and systematics" << std::endl;
89  std::vector<PredictionInterp*> predictionVec;
91  ENumu2017ExtrapType extrap = kNuMu;
92 
93  PredictionSystNumu2017 predQ1(extrap, calc, 1);
94  PredictionSystNumu2017 predQ2(extrap, calc, 2);
95  PredictionSystNumu2017 predQ3(extrap, calc, 3);
96  PredictionSystNumu2017 predQ4(extrap, calc, 4);
97  predictionVec.push_back(&predQ1);
98  predictionVec.push_back(&predQ2);
99  predictionVec.push_back(&predQ3);
100  predictionVec.push_back(&predQ4);
101 
102  /*
103  systs = getAllDirectNumuSysts2017();
104  for(int quantId = 0; quantId < nQuant; quantId++){
105  std::cout << "Quantile " << quantId + 1 << std::endl;
106  std::string input = inDir + inName + Form("Quant%d_", quantId+1) + "full_" + SampleName + ".root";
107  TFile* file = new TFile(input.c_str());
108  std::cout << "\nPrediction vectors" << std::endl;
109  predictionVec.push_back(LoadFrom<PredictionInterp>(file, "prediction")).release();
110  file->Close();
111  std::cout << "Read and closed input files" << std::endl;
112  }//loop quant
113  */
114  std::cout << "\nList of systematics:" << std::endl;
115  for (auto & sys:systs) std::cout <<sys->ShortName() << std::endl;
116 
117 
118  TFile* fout = new TFile(outDir + "plotDataPred.root","RECREATE");
119 
120 
121 
122  // nue oscillation channels
123  std::vector<TH1*> hNuEToNuE;
124  std::vector<TH1*> hNuEToNuMu;
125  std::vector<TH1*> hNuEToNuTau;
126  std::vector<TH1*> hBarNuEToNuE;
127  std::vector<TH1*> hBarNuEToNuMu;
128  std::vector<TH1*> hBarNuEToNuTau;
129  // numu oscillation channels
130  std::vector<TH1*> hNuMuToNuE;
131  std::vector<TH1*> hNuMuToNuMu;
132  std::vector<TH1*> hNuMuToNuTau;
133  std::vector<TH1*> hBarNuMuToNuE;
134  std::vector<TH1*> hBarNuMuToNuMu;
135  std::vector<TH1*> hBarNuMuToNuTau;
136  // others
137  std::vector<TH1*> hNC; // nc only
138  std::vector<TH1*> hPred; // total prediction
139 
140 
141 
142  std::cout << "\n\nFilling histograms" << std::endl;
143  hNuEToNuE.push_back(predictionVec[0]->PredictComponent(calc, Flavors::kNuEToNuE, Current::kCC, Sign::kNu).ToTH1(pot, kGray));
144  hNuEToNuMu.push_back(predictionVec[0]->PredictComponent(calc, Flavors::kNuEToNuMu, Current::kCC, Sign::kNu).ToTH1(pot, kGray));
145  hNuEToNuTau.push_back(predictionVec[0]->PredictComponent(calc, Flavors::kNuEToNuTau, Current::kCC, Sign::kNu).ToTH1(pot, kGray));
146  hBarNuEToNuE.push_back(predictionVec[0]->PredictComponent(calc, Flavors::kNuEToNuE, Current::kCC, Sign::kAntiNu).ToTH1(pot, kGray));
147  hBarNuEToNuMu.push_back(predictionVec[0]->PredictComponent(calc, Flavors::kNuEToNuMu, Current::kCC, Sign::kAntiNu).ToTH1(pot, kGray));
148  hBarNuEToNuTau.push_back(predictionVec[0]->PredictComponent(calc, Flavors::kNuEToNuTau, Current::kCC, Sign::kAntiNu).ToTH1(pot, kGray));
149  hNuMuToNuE.push_back(predictionVec[0]->PredictComponent(calc, Flavors::kNuMuToNuE, Current::kCC, Sign::kNu).ToTH1(pot, kGray));
150  hNuMuToNuMu.push_back(predictionVec[0]->PredictComponent(calc, Flavors::kNuMuToNuMu, Current::kCC, Sign::kNu).ToTH1(pot, kGray));
151  hNuMuToNuTau.push_back(predictionVec[0]->PredictComponent(calc, Flavors::kNuMuToNuTau, Current::kCC, Sign::kNu).ToTH1(pot, kGray));
152  hBarNuMuToNuE.push_back(predictionVec[0]->PredictComponent(calc, Flavors::kNuMuToNuE, Current::kCC, Sign::kAntiNu).ToTH1(pot, kGray));
153  hBarNuMuToNuMu.push_back(predictionVec[0]->PredictComponent(calc, Flavors::kNuMuToNuMu, Current::kCC, Sign::kAntiNu).ToTH1(pot, kGray));
154  hBarNuMuToNuTau.push_back(predictionVec[0]->PredictComponent(calc, Flavors::kNuMuToNuTau, Current::kCC, Sign::kAntiNu).ToTH1(pot, kGray));
155  hNC.push_back(predictionVec[0]->PredictComponent(calc, Flavors::kAll, Current::kNC, Sign::kBoth).ToTH1(pot, kGray));
156  hPred.push_back(predictionVec[0]->Predict(calc).ToTH1(pot, kGray));
157 
158  for(int quantId = 0; quantId < nQuant; quantId++){
159  hNuEToNuE.push_back(predictionVec[quantId]->PredictComponent(calc, Flavors::kNuEToNuE, Current::kCC, Sign::kNu).ToTH1(pot, kGray));
160  hNuEToNuMu.push_back(predictionVec[quantId]->PredictComponent(calc, Flavors::kNuEToNuMu, Current::kCC, Sign::kNu).ToTH1(pot, kGray));
161  hNuEToNuTau.push_back(predictionVec[quantId]->PredictComponent(calc, Flavors::kNuEToNuTau, Current::kCC, Sign::kNu).ToTH1(pot, kGray));
162  hBarNuEToNuE.push_back(predictionVec[quantId]->PredictComponent(calc, Flavors::kNuEToNuE, Current::kCC, Sign::kAntiNu).ToTH1(pot, kGray));
163  hBarNuEToNuMu.push_back(predictionVec[quantId]->PredictComponent(calc, Flavors::kNuEToNuMu, Current::kCC, Sign::kAntiNu).ToTH1(pot, kGray));
164  hBarNuEToNuTau.push_back(predictionVec[quantId]->PredictComponent(calc, Flavors::kNuEToNuTau, Current::kCC, Sign::kAntiNu).ToTH1(pot, kGray));
165  hNuMuToNuE.push_back(predictionVec[quantId]->PredictComponent(calc, Flavors::kNuMuToNuE, Current::kCC, Sign::kNu).ToTH1(pot, kGray));
166  hNuMuToNuMu.push_back(predictionVec[quantId]->PredictComponent(calc, Flavors::kNuMuToNuMu, Current::kCC, Sign::kNu).ToTH1(pot, kGray));
167  hNuMuToNuTau.push_back(predictionVec[quantId]->PredictComponent(calc, Flavors::kNuMuToNuTau, Current::kCC, Sign::kNu).ToTH1(pot, kGray));
168  hBarNuMuToNuE.push_back(predictionVec[quantId]->PredictComponent(calc, Flavors::kNuMuToNuE, Current::kCC, Sign::kAntiNu).ToTH1(pot, kGray));
169  hBarNuMuToNuMu.push_back(predictionVec[quantId]->PredictComponent(calc, Flavors::kNuMuToNuMu, Current::kCC, Sign::kAntiNu).ToTH1(pot, kGray));
170  hBarNuMuToNuTau.push_back(predictionVec[quantId]->PredictComponent(calc, Flavors::kNuMuToNuTau, Current::kCC, Sign::kAntiNu).ToTH1(pot, kGray));
171  hNC.push_back(predictionVec[quantId]->PredictComponent(calc, Flavors::kAll, Current::kNC, Sign::kBoth).ToTH1(pot, kGray));
172  hPred.push_back(predictionVec[quantId]->Predict(calc).ToTH1(pot, kGray));
173 
174  if(quantId > 0){ // already added the first quantile
175  hNuEToNuE[0]->Add(predictionVec[quantId]->PredictComponent(calc, Flavors::kNuEToNuE, Current::kCC, Sign::kNu).ToTH1(pot, kGray));
176  hNuEToNuMu[0]->Add(predictionVec[quantId]->PredictComponent(calc, Flavors::kNuEToNuMu, Current::kCC, Sign::kNu).ToTH1(pot, kGray));
177  hNuEToNuTau[0]->Add(predictionVec[quantId]->PredictComponent(calc, Flavors::kNuEToNuTau, Current::kCC, Sign::kNu).ToTH1(pot, kGray));
178  hBarNuEToNuE[0]->Add(predictionVec[quantId]->PredictComponent(calc, Flavors::kNuEToNuE, Current::kCC, Sign::kAntiNu).ToTH1(pot, kGray));
179  hBarNuEToNuMu[0]->Add(predictionVec[quantId]->PredictComponent(calc, Flavors::kNuEToNuMu, Current::kCC, Sign::kAntiNu).ToTH1(pot, kGray));
180  hBarNuEToNuTau[0]->Add(predictionVec[quantId]->PredictComponent(calc, Flavors::kNuEToNuTau, Current::kCC, Sign::kAntiNu).ToTH1(pot, kGray));
181  hNuMuToNuE[0]->Add(predictionVec[quantId]->PredictComponent(calc, Flavors::kNuMuToNuE, Current::kCC, Sign::kNu).ToTH1(pot, kGray));
182  hNuMuToNuMu[0]->Add(predictionVec[quantId]->PredictComponent(calc, Flavors::kNuMuToNuMu, Current::kCC, Sign::kNu).ToTH1(pot, kGray));
183  hNuMuToNuTau[0]->Add(predictionVec[quantId]->PredictComponent(calc, Flavors::kNuMuToNuTau, Current::kCC, Sign::kNu).ToTH1(pot, kGray));
184  hBarNuMuToNuE[0]->Add(predictionVec[quantId]->PredictComponent(calc, Flavors::kNuMuToNuE, Current::kCC, Sign::kAntiNu).ToTH1(pot, kGray));
185  hBarNuMuToNuMu[0]->Add(predictionVec[quantId]->PredictComponent(calc, Flavors::kNuMuToNuMu, Current::kCC, Sign::kAntiNu).ToTH1(pot, kGray));
186  hBarNuMuToNuTau[0]->Add(predictionVec[quantId]->PredictComponent(calc, Flavors::kNuMuToNuTau, Current::kCC, Sign::kAntiNu).ToTH1(pot, kGray));
187  hNC[0]->Add(predictionVec[quantId]->PredictComponent(calc, Flavors::kAll, Current::kNC, Sign::kBoth).ToTH1(pot, kGray));
188  hPred[0]->Add(predictionVec[quantId]->Predict(calc).ToTH1(pot, kGray));
189  } // for quantId > 0
190  } // for quantId
191 
192 
193 
194 
195  std::cout << "\n\nDeclaring integrals for predicted number of events" << std::endl;
196  // nue oscillation channels
197  std::vector<double> iNuEToNuE;
198  std::vector<double> iNuEToNuMu;
199  std::vector<double> iNuEToNuTau;
200  std::vector<double> iBarNuEToNuE;
201  std::vector<double> iBarNuEToNuMu;
202  std::vector<double> iBarNuEToNuTau;
203  // numu oscillation channels
204  std::vector<double> iNuMuToNuE;
205  std::vector<double> iNuMuToNuMu;
206  std::vector<double> iNuMuToNuTau;
207  std::vector<double> iBarNuMuToNuE;
208  std::vector<double> iBarNuMuToNuMu;
209  std::vector<double> iBarNuMuToNuTau;
210  // others
211  std::vector<double> iNC; // nc only
212  std::vector<double> iPred; // total prediction
213  std::vector<double> iBkg; // beam background
214  std::vector<double> iCos; // cosmic background
215  std::vector<double> iNuMu; // numu + numubar
216 
217  std::cout << "\n\nGetting integrals for predicted number of events" << std::endl;
218  for(int quantId = 0; quantId < nQuantPlus; quantId++){
219  std::cout<< "quant" << quantId << std::endl;
220 
221  iNuEToNuE.push_back(hNuEToNuE[quantId]->Integral());
222  iNuEToNuMu.push_back(hNuEToNuMu[quantId]->Integral());
223  iNuEToNuTau.push_back(hNuEToNuTau[quantId]->Integral());
224  iBarNuEToNuE.push_back(hBarNuEToNuE[quantId]->Integral());
225  iBarNuEToNuMu.push_back(hBarNuEToNuMu[quantId]->Integral());
226  iBarNuEToNuTau.push_back(hBarNuEToNuTau[quantId]->Integral());
227  iNuMuToNuE.push_back(hNuMuToNuE[quantId]->Integral());
228  iNuMuToNuMu.push_back(hNuMuToNuMu[quantId]->Integral());
229  iNuMuToNuTau.push_back(hNuMuToNuTau[quantId]->Integral());
230  iBarNuMuToNuE.push_back(hBarNuMuToNuE[quantId]->Integral());
231  iBarNuMuToNuMu.push_back(hBarNuMuToNuMu[quantId]->Integral());
232  iBarNuMuToNuTau.push_back(hBarNuMuToNuTau[quantId]->Integral());
233  iNC.push_back(hNC[quantId]->Integral());
234  iCos.push_back(hCosmics[quantId]->Integral());
235 
236  float tempNuMu = iNuEToNuMu[quantId] + iBarNuEToNuMu[quantId] + iNuMuToNuMu[quantId] + iBarNuMuToNuMu[quantId];
237  float tempPred = hPred[quantId]->Integral() + iCos[quantId];
238  float tempBkg = tempPred - (iCos[quantId] + iNC[quantId] + tempNuMu);
239 
240  iPred.push_back(tempPred);
241  iNuMu.push_back(tempNuMu);
242  iBkg.push_back(tempBkg);
243 
244  } // for quantId
245 
246 
247  std::cout << "\n\nMake tables for " << CalcName << "\n" << std::endl;
248  for(int quantId = 0; quantId < nQuantPlus; quantId++){
249 
250  cout.setf(ios::fixed, ios::floatfield);
251  cout.width(6);
252  cout.precision(3);
253 
254  std::cout << "\n\nslashbegin{frame}" << std::endl;
255  std::cout << "\n\n" << CalcName << " Quantile " << quantId << "\n" << std::endl;
256 
257  std::cout << "slashbegin{table}" << std::endl;
258  std::cout << "slashcentering" << std::endl;
259  std::cout << "slashbegin{tabular}{cccc}" << std::endl;
260  std::cout << "slashhline" << std::endl;
261  std::cout << "slashmulticolumn{1}{|c}{Channel} & slashmulticolumn{1}{c|}{Selected events} & Channel & slashmulticolumn{1}{c|}{Selected events} doubleslash slashhline" << std::endl;
262  std::cout << "slashmulticolumn{1}{|c}{" << std::endl;
263  std::cout << "slashnumuD slasharrow slashnumuD} & slashmulticolumn{1}{c|}{" << iNuMuToNuMu[quantId] << "} & slashnueD slasharrow slashnumuD & slashmulticolumn{1}{c|}{" << iNuEToNuMu[quantId] << "} doubleslash" << std::endl;
264  std::cout << "slashmulticolumn{1}{|c}{slashnumuD slasharrow slashnueD} & slashmulticolumn{1}{c|}{" << iNuMuToNuE[quantId] << "} & slashnueD slasharrow slashnueD & slashmulticolumn{1}{c|}{" << iNuEToNuE[quantId] << "} doubleslash" << std::endl;
265  std::cout << "slashmulticolumn{1}{|c}{ slashnumuD slasharrow slashnutauD}& slashmulticolumn{1}{c|}{" << iNuMuToNuTau[quantId] << "} & slashnueD slasharrow slashnutauD & slashmulticolumn{1}{c|}{" << iNuEToNuTau[quantId] << "} doubleslash" << std::endl;
266  std::cout << "slashmulticolumn{1}{|c}{slashnumubarDslasharrowslashnumubarD} & slashmulticolumn{1}{c|}{" << iBarNuMuToNuMu[quantId] << "} & slashnuebarD slasharrow slashnumubarD & slashmulticolumn{1}{c|}{" << iBarNuEToNuMu[quantId] << "} doubleslash" << std::endl;
267  std::cout << "slashmulticolumn{1}{|c}{slashnumubarD slasharrow slashnuebarD} & slashmulticolumn{1}{c|}{" << iBarNuMuToNuE[quantId] << "} & slashnuebarD slasharrow slashnuebarD & slashmulticolumn{1}{c|}{" << iBarNuEToNuE[quantId] << "} doubleslash" << std::endl;
268  std::cout << "slashmulticolumn{1}{|c}{slashnumubarD slasharrow slashnutaubarD} & slashmulticolumn{1}{c|}{" << iBarNuMuToNuTau[quantId] << "} & slashnuebarD slasharrow slashnutaubarD & slashmulticolumn{1}{c|}{" << iBarNuEToNuTau[quantId] << "} doubleslash" << std::endl;
269  std::cout << "slashhline" << std::endl;
270  std::cout << "slashhline" << std::endl;
271  std::cout << "slashmulticolumn{1}{|c|}{slashnumuD + slashnumubarD signal} & slashmulticolumn{1}{c|}{NC} & slashmulticolumn{1}{c|}{Other beam bkg} & slashmulticolumn{1}{c|}{Cosmics} doubleslash" << std::endl;
272  std::cout << "slashmulticolumn{1}{|c|}{" << iNuMu[quantId] << "} & slashmulticolumn{1}{c|}{" << iNC[quantId] << "} & slashmulticolumn{1}{c|}{" << iBkg[quantId] << "} & slashmulticolumn{1}{c|}{" << iCos[quantId] << "} doubleslash" << std::endl;
273  std::cout << "slashhline" << std::endl;
274  std::cout << "slashhline" << std::endl;
275  std::cout << "Total: " << iPred[quantId] << std::endl;
276  std::cout << "slashend{tabular}" << std::endl;
277  std::cout << "slashend{table}" << std::endl;
278 
279  std::cout << "slashend{frame}\n\n" << std::endl;
280 
281  } // for quantId
282 
283 
284 
285 
286 } // End of function
void makePredTables(std::string calculator="maxmix")
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
Antineutrinos-only.
Definition: IPrediction.h:50
std::vector< SystGroupDef > systs
Definition: syst_header.h:385
(&#39; appearance&#39;)
Definition: IPrediction.h:18
double Integral(const Spectrum &s, const double pot, int cvnregion)
(&#39;beam &#39;)
Definition: IPrediction.h:15
T sqrt(T number)
Definition: d0nt_math.hpp:156
const double kAna2017POT
Definition: Exposures.h:174
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
TString inDir
Charged-current interactions.
Definition: IPrediction.h:39
const double kAna2017Livetime
Definition: Exposures.h:200
Loads shifted spectra from files.
#define pot
const int nQuantPlus
Definition: varsandcuts.h:5
Neutrinos-only.
Definition: IPrediction.h:49
OStream cout
Definition: OStream.cxx:6
(&#39; survival&#39;)
Definition: IPrediction.h:19
double livetime
Definition: saveFDMCHists.C:21
const int nQuant
Definition: varsandcuts.h:4
virtual void SetTh23(const T &th23)=0
Neutral-current interactions.
Definition: IPrediction.h:40
Both neutrinos and antineutrinos.
Definition: IPrediction.h:52
std::vector< const ISyst * > getAllAna2017Systs(const EAnaType2017 ana, const bool smallgenie)
All neutrinos, any flavor.
Definition: IPrediction.h:26
(&#39; appearance&#39;)
Definition: IPrediction.h:16
Float_t e
Definition: plot.C:35
T asin(T number)
Definition: d0nt_math.hpp:60
enum BeamMode string