syst_table_from_preds.C
Go to the documentation of this file.
5 #include "CAFAna/Prediction/PredictionInterpTemplates.h"
7 #include "CAFAna/Systs/NueSystsSecondAna.h"
8 using namespace ana;
9 
10 #include "OscLib/IOscCalc.h"
11 #include "OscLib/OscCalcDumb.h"
12 
13 #include "TCanvas.h"
14 
15 #include <iomanip>
16 #include <map>
17 
18 
19 
20 
21 // A function to make "key" string in the map unique
22 const std::string makeKeyUnique (const std::string& syst, const std::string& shift)
23 {
25  key = syst+"_"+shift;
26  std::cout << "makeUniqueKey : " << key << std::endl;
27  return key;
28 }
29 
30 struct ColumnDef
31 {
36 };
37 
38 const std::vector<ColumnDef> cols =
40  {"beam NueCC", Flavors::kNuEToNuE, Current::kCC, Sign::kBoth},
41  {"NC", Flavors::kAll, Current::kNC, Sign::kBoth},
42  {"Numu CC", Flavors::kNuMuToNuMu, Current::kCC, Sign::kBoth}};
43 
44 
45 // A function to calculate integral from a spectrum for the specific CVN region you choose
46 double Integral(const Spectrum& s, const double pot, int cvnregion)
47 {
48  std::unique_ptr<TH1D> h(s.ToTH1(pot));
49 
50  if(cvnregion == 1) return h->Integral(1, 9);
51  if(cvnregion == 2) return h->Integral(10, 18);
52  if(cvnregion == 3) return h->Integral(19, 27);
53  if(cvnregion == 4) return h->Integral(28, 36);
54  if(cvnregion == 5) return s.Integral(pot);
55 }
56 
57 
58 // A function to calculate the total bkg. for nue appearance
59 double TotBkg(std::vector <double> integrals)
60 {
61  double totbkg = 0;
62  for(auto it: integrals){
63  totbkg += it;
64  }
65  return totbkg;
66 }
67 
68 
69 void syst_table_from_preds(int cvnregion, std::string fname = "")
70 {
71  // A function to switch between different PID region
72  switch (cvnregion)
73  {
74  case 1 :
75  std::cout << "We are in Region 0.75 < CVN < 0.87" << std::endl;
76  break;
77  case 2 :
78  std::cout << "We are in Region 0.87 < CVN < 0.95" << std::endl;
79  break;
80  case 3 :
81  std::cout << "We are in Region 0.95 < CVN <= 1.0" << std::endl;
82  break;
83  case 4 :
84  std::cout<< "We are in Peripheral Region" <<std::endl;
85  break;
86  case 5 :
87  std::cout<< "We are in CVN Region > 0.75" <<std::endl;
88  break;
89  default: std::cout << "Run me like <cafe 'syst_table_from_preds.C+(1/2/3/4/5)'> for CVN 1,2,3/peripheral/CVN >.75 region you want to study. \n"
90  << "1 is for 0.75 < CVN1 < 0.87; 2 is for 0.87 < CVN < 0.95; 3 is for CVN > 1.0; 4 is for peripheral region; 5 is for CVN > 0.75. \n"
91  << "Exiting(0)" << std::endl;
92  return;
93  }
94 
95 
96  osc::OscCalcDumb dumbosc;
97 
98  const double kPOT= 9.0000e20;
99 
100  const int kNumPreds = 2;
101  const std::string predNames[kNumPreds] = {"pred_propDecomp", "pred_comboDecomp"};
102 
103  const int kNumSysts = 5;
104  const std::string systNames[kNumSysts] = {"Nominal", "Calibration", "Lightlevel", "Cherenkov", "CalibShape"};
105 
106  const int kNumShifts = 3;
107  const std::string shiftNames[kNumShifts] = {"noShift","plusOne", "minusOne"};
108 
109  const std::string varName = {"EnergyCVN2D"};
110 
111  std::unique_ptr <PredictionExtrap> predsPropNom, predsPropShifted;
112 
113  // map to hold predictions from all systs and then play with it
114  std::map <const std::string, std::unique_ptr <PredictionExtrap> > map_to_preds;
115  std::map <const std::string, std::unique_ptr <PredictionExtrap> >::iterator map_it;
116 
117 
118  // I think the preds of 2017 ana will only consider propDecomp but you can loop over "kNumPreds" in future.
119  for(int predIdx = 0; predIdx < 1; ++predIdx)
120  {
121  const char* predName = predNames[predIdx].c_str();
122  std::cout << "predName :" << predName << std::endl;
123 
124  predsPropNom = LoadFromFile<PredictionExtrap>(fname, TString::Format("%s/%s/%s/nue_pred_%s", predName, systNames[0].c_str(), shiftNames[0].c_str(), varName.c_str()).Data());
125 
126  for(int systIdx = 1; systIdx < kNumSysts; ++systIdx)
127  {
128  const char* systName = systNames[systIdx].c_str();
129  std::cout << "systName :" << systName << std::endl;
130 
131  for(int shiftIdx = 1; shiftIdx < kNumShifts; ++shiftIdx)
132  {
133  const char* shiftName = shiftNames[shiftIdx].c_str();
134  std::cout << "shifttName :" << shiftName << std::endl;
135 
136  if ((systIdx == 3 || systIdx == 4) && shiftIdx == 2) continue;
137 
138  predsPropShifted = LoadFromFile<PredictionExtrap>(fname, TString::Format("%s/%s/%s/nue_pred_%s", predName, systName, shiftName, varName.c_str()).Data());
139  std ::cout << TString::Format("%s/%s/%s/nue_pred_%s", predName, systName, shiftName, varName.c_str()) << std::endl;
140 
141  // Insert predictions for all syst. samples in the map.
142  map_to_preds.insert(make_pair(makeKeyUnique(systNames[systIdx],shiftNames[shiftIdx]), std::move(predsPropShifted)));
143 
144  } // end loop of shiftIdx
145  } // end loop of systIdx
146  } // end loop of predIdx
147 
148 
149 
150  std::cout << "*********************** For Nominal ************************ " << std::endl;
151 
152  std::vector <double> nomintegrals;
153  for(ColumnDef col: cols)
154  {
155  auto spec = predsPropNom->PredictComponent(&dumbosc, col.flav, col.curr, col.sign);
156 
157  std::cout << col.name << " : " << Integral (spec,kPOT,cvnregion) << std::endl;
158 
159  if ((col.name == "beam NueCC") || (col.name == "NC") || (col.name == "Numu CC"))
160  nomintegrals.push_back (Integral(spec,kPOT,cvnregion));
161 
162  }
163  std::cout << "Total bkg." << " : " << TotBkg(nomintegrals) << std::endl;
164 
165 
166 
167  std::cout << "*********************** For Systematic Sample ************************ " << std::endl;
168 
169  for (map_it = map_to_preds.begin(); map_it != map_to_preds.end(); ++map_it)
170  {
171  auto key = map_it->first;
172  auto pred = map_it->second.get();
173 
174  std::vector <double> sysintegrals;
175  std::cout << "Key/'Name of the Syst.' = " << key << " :: Value/'Check the valid address' = " << pred <<std::endl;
176 
177  for(ColumnDef col: cols)
178  {
179  auto spec = pred->PredictComponent(&dumbosc, col.flav, col.curr, col.sign);
180 
181  std::cout << col.name << " : " << Integral (spec,kPOT,cvnregion) << std::endl;
182 
183  if ((col.name == "beam NueCC") || (col.name == "NC") || (col.name == "Numu CC"))
184  sysintegrals.push_back (Integral(spec,kPOT,cvnregion));
185  }
186 
187  std::cout << "Total bkg." << " : " << TotBkg(sysintegrals) << std::endl;
188  std::cout << " --------------------------- moving to next systematics -----------------" << std::endl;
189  std::cout << " ----------------------------------------------------------------------------" << std::endl;
190 
191  } // loop for map_it
192 
193 
194 }
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
set< int >::iterator it
(&#39; appearance&#39;)
Definition: IPrediction.h:18
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
double Integral(const Spectrum &s, const double pot, int cvnregion)
(&#39;beam &#39;)
Definition: IPrediction.h:15
double TotBkg(std::vector< double > integrals)
std::pair< Spectrum *, CheatDecomp * > make_pair(SpectrumLoaderBase &loader_data, SpectrumLoaderBase &loader_mc, HistAxis *axis, Cut *cut, const SystShifts &shift, const Var &wei)
Definition: DataMCLoad.C:336
const std::string makeKeyUnique(const std::string &syst, const std::string &shift)
double Integral(double exposure, double *err=0, EExposureType expotype=kPOT) const
Return total number of events scaled to pot.
Definition: Spectrum.cxx:249
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
const XML_Char * s
Definition: expat.h:262
Current::Current_t curr
Charged-current interactions.
Definition: IPrediction.h:39
Int_t col[ntarg]
Definition: Style.C:29
#define pot
Simple oscillation probability calculator that has no solar term or mass hierarchy or delta so it&#39;s s...
Definition: OscCalcDumb.h:16
void syst_table_from_preds(int cvnregion, std::string fname="")
std::string name
OStream cout
Definition: OStream.cxx:6
(&#39; survival&#39;)
Definition: IPrediction.h:19
Sign::Sign_t sign
std::map< std::string, std::string > systNames
Definition: PlotUnfolding.C:32
Neutral-current interactions.
Definition: IPrediction.h:40
const std::vector< ColumnDef > cols
Both neutrinos and antineutrinos.
Definition: IPrediction.h:52
Flavors::Flavors_t flav
void Format(TGraph *gr, int lcol, int lsty, int lwid, int mcol, int msty, double msiz)
Definition: Style.cxx:154
All neutrinos, any flavor.
Definition: IPrediction.h:26
enum BeamMode string