multiverse_macro.C
Go to the documentation of this file.
1 /*
2 #include "CAFAna/Core/Binning.h"
3 #include "CAFAna/Cuts/Cuts.h"
4 
5 #include "CAFAna/Core/Spectrum.h"
6 #include "CAFAna/Core/SpectrumLoader.h"
7 
8 #include "CAFAna/Core/Var.h"
9 #include "CAFAna/Core/MultiVar.h"
10 #include "CAFAna/Core/HistAxis.h"
11 #include "CAFAna/Core/SystShifts.h"
12 
13 #include "CAFAna/Vars/Vars.h"
14 //#include "CAFAna/Vars/NueVarsExtra.h"
15 
16 #include "NDAna/numucc_1Pi/NumuCC1PiVars.h"
17 #include "NDAna/numucc_1Pi/NumuCC1PiCuts.h"
18 
19 #include "G4RwgtPlots/Geant4WeightVars.h"
20 
21 #include "CAFAna/XSec/Multiverse.h"
22 #include "CAFAna/XSec/SystematicDef.h"
23 //#include "CAFAna/XSec/CutOptimization.h"
24 
25 // for POT
26 #include "CAFAna/Analysis/Exposures.h"
27 
28 // for plotting
29 #include "CAFAna/Analysis/Plots.h"
30 
31 // to get datetime in plot names
32 #include "/nova/app/users/csweeney/xsec/code/Utilities/misc.h"
33 
34 #include "StandardRecord/Proxy/SRProxy.h"
35 
36 #include <vector>
37 #include <string>
38 
39 #include "TFile.h"
40 #include "TH1.h"
41 #include "TCanvas.h"
42 #include "TGraphAsymmErrors.h"
43 #include "TString.h"
44 */
45 
46 #include "CAFAna/Core/Binning.h"
47 #include "CAFAna/Cuts/Cuts.h"
48 
49 #include "CAFAna/Core/Spectrum.h"
51 
52 #include "CAFAna/Core/Var.h"
53 
54 #include "CAFAna/Vars/Vars.h"
55 
56 #include "CAFAna/XSec/Multiverse.h"
58 
60 
61 #include <vector>
62 #include <string>
63 
64 #include "TFile.h"
65 #include "TH1.h"
66 #include "TCanvas.h"
67 #include "TGraphAsymmErrors.h"
68 #include "TString.h"
69 
72 
73 #include "G4RwgtPlots/Geant4WeightVars.h"
74 
75 // for POT
77 
78 // for plotting
79 #include "CAFAna/Analysis/Plots.h"
80 
81 // to get datetime in plot names
82 #include "/nova/app/users/csweeney/xsec-TR/Utilities/misc.h"
83 
84 // for kVtxX
86 
87 
88 namespace ana
89 {
90 
91 
93  ToUpDownHist(Multiverse * mv, const TH1 * h_nominal)
94  //Scavenged from CAFAna/XSec/CutOptimization.h
95  {
96  TH1 * h_up = mv->GetPlusOneSigmaShift (h_nominal);
97  TH1 * h_down = mv->GetMinusOneSigmaShift(h_nominal);
98  // transform into 1 sided shift as the most conservative shift
99  TH1 * ret = (TH1*) h_up->Clone();
100  for(auto ibinx = 1; ibinx <= ret->GetNbinsX(); ibinx++) {
101  for(auto ibiny = 1; ibiny <= ret->GetNbinsY(); ibiny++) {
102  for(auto ibinz = 1; ibinz <= ret->GetNbinsZ(); ibinz++) {
103  double up = h_up ->GetBinContent(ibinx, ibiny, ibinz);
104  double down = h_down ->GetBinContent(ibinx, ibiny, ibinz);
105  double nominal = h_nominal ->GetBinContent(ibinx, ibiny, ibinz);
106  double abs_shift = std::max(std::abs(up - nominal), std::abs(down - nominal));
107  ret->SetBinContent(ibinx, ibiny, ibinz, abs_shift + nominal);
108  }
109  }
110  }
111  return UpDownPair<TH1>{ret};
112  }
113 
114 
115  void PlotDebug(const std::vector<TH1 *> universes,
116  TH1 * h_nominal,
119  //Scavenged from CAFAna/XSec/CutOptimization.h
120  {
121  TDirectory * tmp = gDirectory;
122  TCanvas * c = new TCanvas();
123  c->SetLeftMargin(0.15);
124  c->SetRightMargin(0.01);
125  c->cd();
126  double max_val = -1e9;
127  for(auto ihist = 0u; ihist < universes.size(); ihist++) {
128  max_val = std::max(max_val, universes[ihist]->GetMaximum());
129  }
130  h_nominal->GetYaxis()->SetRangeUser(0, max_val * 1.2);
131  h_nominal->SetTitle(title.c_str());
132  h_nominal->Draw("hist");
133  for(auto ihist = 0u; ihist < universes.size(); ihist++) {
134  universes[ihist]->Draw("hist same");
135  }
136  h_nominal->Draw("hist same");
137 
138  c->Print(name.c_str());
139  tmp->cd();
140  }
141 
142 
143 // extern const Var kWeight([](const caf::SRProxy* sr){
144 // //Just prints out weights for debugging purposes
145 // if(sr->mc.nnu != 1 ) return -5.0;
146 //
147 // std::cout << "foo" << std::endl;
148 // auto weight = sr->mc.nu[0].geantWeights[0].piplusweight;
149 // std::cout << weight << std::endl;
150 //
151 // /*
152 // int nUniv = weightTable.size();
153 // std::cout << "nUniv: " << nUniv << std::endl;
154 // */
155 //
156 // /*
157 // for(int i=0; i<nChannels; ++i){
158 // std::cout << i << std::endl;
159 // std::cout << weightTable[i].piplus << std::endl;
160 //
161 // }
162 // */
163 // std::cout << "----------------" << std::endl;
164 //
165 // return -5.0;
166 // });
167 
168 
169 
170 
171 
172 
173 
174  Var kUnity([](const caf::SRProxy *sr){
175  // To be used as "common_weight" for Multiverse constructor
176  return 1.f;
177  });
178 
179 
180 }
181 
182 
183 using namespace ana;
184 
186 {
187 
188  //std::string inName = "/nova/app/users/csweeney/multirw_files/9file.multirw.caf.root";
189  std::string inName = "/pnfs/nova/scratch/users/csweeney/final_grid_prod5_g4rwgt_miniprod/*";
190 
191  SpectrumLoader* loader = new SpectrumLoader(inName);
192 
193  //std::string outName = "test.root";
194  //TFile *outFile = new TFile(outName.c_str(), "RECREATE");
195 
196  const SystShifts * shift = new SystShifts(kNoShift);
197  const Var * common_weight = new Var(kUnity);
198 
199 
200 
201  //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
202  std::vector<std::vector<Var>> weights_vec{
203  GetkGeantPiplusWeights(), GetkGeantPiplusWeights(), GetkGeantPiplusWeights(),
204  GetkGeantPiminusWeights(), GetkGeantPiminusWeights(), GetkGeantPiminusWeights(),
205  GetkGeantAllPionsWeights(), GetkGeantAllPionsWeights(), GetkGeantAllPionsWeights()
206  };
207  std::vector<std::string> myLabel_vec{
208  "Vertex in X Direction (cm)", "Vertex in Y Direction (cm)", "Vertex in Z Direction (cm)",
209  "Vertex in X Direction (cm)", "Vertex in Y Direction (cm)", "Vertex in Z Direction (cm)",
210  "Vertex in X Direction (cm)", "Vertex in Y Direction (cm)", "Vertex in Z Direction (cm)"
211  };
212  std::vector<std::string> myTitle_vec{
213  "Piplus 100 universes", "Piplus 100 universes", "Piplus 100 universes",
214  "Piminus 100 universes", "Piminus 100 universes", "Piminus 100 universes",
215  "Total 100 universes", "Total 100 universes", "Total 100 universes"
216  };
217  std::vector<std::string> myName_vec = {
218  "eff_piplus_vtxX", "eff_piplus_vtxY", "eff_piplus_vtxZ",
219  "eff_piminus_vtxX", "eff_piminus_vtxY", "eff_piminus_vtxZ",
220  "eff_total_vtxX", "eff_total_vtxY", "eff_total_vtxZ"
221  };
222  std::vector<Binning> myBins_vec{
223  Binning::Simple(10, -200, 200), Binning::Simple(10, -200, 200), Binning::Simple(10, 0, 1270),
224  Binning::Simple(10, -200, 200), Binning::Simple(10, -200, 200), Binning::Simple(10, 0, 1270),
225  Binning::Simple(10, -200, 200), Binning::Simple(10, -200, 200), Binning::Simple(10, 0, 1270)
226  };
227  std::vector<const Var *> kMyVar_vec{
228  new Var(kVtxX), new Var(kVtxY), new Var(kVtxZ),
229  new Var(kVtxX), new Var(kVtxY), new Var(kVtxZ),
230  new Var(kVtxX), new Var(kVtxY), new Var(kVtxZ)
231  };
233  //const Cut * kSignal = new Cut(CutFromNuTruthCut(kTrueDetector_NT) );
234 
235 
236  //-----
237 
238  std::vector<Multiverse *> myMV_vec;
239  std::vector<Spectrum *> mySpec_vec;
240 
241 
242 
243  for(uint i=0; i<myName_vec.size(); ++i){
244 
245  std::string myLabel = myLabel_vec[i];
246  Binning myBins = myBins_vec[i];
247  const Var * kMyVar = kMyVar_vec[i];
248  std::vector<Var> weights = weights_vec[i];
249 
250  // std::string myName = myName_vec[i];
251 
252 
253 
254  const HistAxis * axis = new HistAxis(myLabel,
255  myBins,
256  *kMyVar);
257 
258 
259  Multiverse * myMV = new Multiverse(*loader,
260  *axis,
261  *kSignal,
262  *shift,
263  weights,
264  *common_weight);
265 
266 
267  Spectrum * sSig = new Spectrum(myLabel, myBins,
268  *loader, *kMyVar,
269  *kSignal);
270  myMV_vec.push_back(myMV);
271  mySpec_vec.push_back(sSig);
272  }
273 
274  //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
275 
276 
277 
278 
279  loader->Go();
280 
281 
282  for(uint i=0; i<myMV_vec.size(); ++i){
283 
284  Multiverse * myMV = myMV_vec[i];
285  Spectrum * sSig = mySpec_vec[i];
286 
287  std::string myName = myName_vec[i];
288  std::string myLabel = myLabel_vec[i];
289  std::string myTitle = myTitle_vec[i];
290 
291 
292  //double pot = kAna2020FHCPOT;
293  double pot = sSig->POT();
294 
295  //TDirectory * dir = outFile->mkdir("vtx_x");
296  //TDirectory * nom_dir = dir->mkdir("Nominal");
297  //fNominalSignal->SaveTo(nom_dir->mkdir("Signal"));
298  //TDirectory * mv_dir = dir->mkdir("mv_geant4");
299  //myMV->SaveTo( mv_dir->mkdir("fMVSignal") );
300 
301  //------------------------------------------
302 
303  //TH1 * hSig =(TH1*) fNominalSignal->ToTH1(pot);
304  //std::vector<UpDownPair<TH1> > hSystSignal;
305  //hSystSignal.push_back(ToUpDownHist(myMV,hSig) );
306  TH1 * hSig = sSig->ToTH1(pot);
307 
308  /*
309  hSig->SetAxisRange(0., 2.0*hSig->GetMaximum(), "Y" );
310  const std::vector<TH1*> mv_signal_hists = myMV->GetUniversesHist(kCyan);
311  PlotDebug(mv_signal_hists, hSig, "this is a plot", "plot.png");
312  */
313 
314  Spectrum * sUp = myMV->GetPlusOneSigmaShift(sSig);
315  TH1 * hUp = sUp->ToTH1(pot);
316  std::vector<TH1*> hUp_vec{hUp}; // for some reason syst plotting function wants vector of TH1s
317  //hUp->SetLineColor(kRed);
318 
319  Spectrum * sDown = myMV->GetMinusOneSigmaShift(sSig);
320  TH1 * hDown = sDown->ToTH1(pot);
321  std::vector<TH1*> hDown_vec{hDown};
322  //hDown->SetLineColor(kBlue);
323 
324  // Get fractional uncertainty (take larger of +/-)
325  int nBins = hSig->GetNbinsX();
326  double xLow = hSig->GetXaxis()->GetXmin();
327  double xHigh = hSig->GetXaxis()->GetXmax();
328 
329  TH1D * fracUncert = new TH1D("frac", "", nBins, xLow, xHigh);
330  for(int i=1; i <= nBins; ++i){
331  double up = hUp->GetBinContent(i);
332  double down = hDown->GetBinContent(i);
333  double nom = hSig->GetBinContent(i);
334 
335  double abs_shift = std::max( std::abs(up - nom), std::abs(down - nom) );
336  double frac = abs_shift / nom;
337  //std::cout << "{ " << up << ", " << down << " }" << std::endl;
338  //std::cout << nom << " : " << abs_shift << " : " << frac << std::endl;
339  fracUncert->SetBinContent(i, frac);
340  }
341 
342  fracUncert->GetXaxis()->SetTitle(myLabel.c_str());
343 
344 
345  /*
346  TCanvas * canv2 = new TCanvas("canv2", "canv2", 800, 600);
347  fracUncert->Draw();
348  */
349 
350  TGraphAsymmErrors * graph = PlotWithSystErrorBand(hSig, hUp_vec, hDown_vec);
351  graph->SetTitle(myTitle.c_str());
352 
353  //TCanvas * canv = new TCanvas("canv", "canv", 800, 600);
354 
355  hSig->SetLineColor(1);
356  hSig->SetLineWidth(2);
357  /*
358  hSig->Draw("hist ][");
359  graph->Draw("e2 same");
360  hSig->Draw("hist ][ same");
361 
362  canv->SaveAs("error_band.png");
363  */
364 
365  TString rCname = "rC";
366  TCanvas* rC = new TCanvas(rCname, rCname);
367 
368  rC -> SetBottomMargin(0.);
369  double Spl = 0.3;
370  TPad* P1 = new TPad( "Temp_1", "", 0.0, Spl, 1.0, 1.0, 0 );
371  TPad* P2 = new TPad( "Temp_2", "", 0.0, 0.0, 1.0, Spl, 0 );
372  P2 -> SetRightMargin (.03);
373  P2 -> SetTopMargin (.00);
374  P2 -> SetBottomMargin(.3);
375  P2 -> SetLeftMargin (.1);
376  P2 -> Draw();
377  P1 -> SetRightMargin (.03);
378  P1 -> SetLeftMargin (.1);
379  P1 -> SetTopMargin (.10);
380  P1 -> SetBottomMargin(.00);
381  P1 -> Draw();
382  // Set some label sizes.
383  double Lb1 = 0.07;
384  double Lb2 = 0.15;
385  // --- First, draw the fracUncert so cd onto Pad2
386  P2 -> cd();
387 
388  // Set axis ranges etc.
389  fracUncert->GetYaxis()->SetTitleSize( Lb2 );
390  fracUncert->GetYaxis()->SetTitleOffset(0.4);
391  fracUncert->GetYaxis()->SetLabelSize( Lb2 );
392  fracUncert->GetXaxis()->SetTitleSize( Lb2 );
393  fracUncert->GetXaxis()->SetLabelSize( Lb2 );
394  fracUncert->Draw("");
395 
396  //fracUncert->Draw("axis same");
397  P1->cd();
398  // graph goes on P1
399  hSig->SetTitle(myTitle.c_str());
400  hSig->GetYaxis()->SetTitleSize( Lb1 );
401  hSig->GetYaxis()->SetLabelSize( Lb1 );
402  hSig->GetYaxis()->SetTitleOffset( 0.7 );
403  // Remove the x axis labels
404  hSig->GetXaxis()->SetLabelSize (0 );
405  hSig->GetXaxis()->SetLabelOffset(99);
406  hSig->SetLineColor(4);
407  hSig->DrawCopy("hist");
408 
409  hSig->Draw("hist ][");
410  graph->Draw("e2 same");
411  hSig->Draw("hist ][ same");
412 
413  rC->SaveAs(("./plots/multiverse/" + DateTime() + myName + ".png").c_str());
414 
415  }
416 
417 
418 }
419 
T max(const caf::Proxy< T > &a, T b)
const XML_Char * name
Definition: expat.h:151
tree Draw("slc.nhit")
int nBins
Definition: plotROC.py:16
_HistAxis< Var > HistAxis
Definition: HistAxis.h:103
Represent the binning of a Spectrum&#39;s x-axis.
Definition: Binning.h:16
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
Spectrum * GetPlusOneSigmaShift(const Spectrum *)
Definition: Multiverse.cxx:443
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
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
TGraphAsymmErrors * PlotWithSystErrorBand(IPrediction *pred, const std::vector< const ISyst * > &systs, osc::IOscCalc *calc, double pot, int col, int errCol, float headroom, bool newaxis, EBinType bintype, double alpha)
Plot prediction with +/-1sigma error band.
Definition: Plots.cxx:304
Proxy for caf::StandardRecord.
Definition: SRProxy.h:2126
Float_t tmp
Definition: plot.C:36
const Var kVtxX
Spectrum * GetMinusOneSigmaShift(const Spectrum *)
Definition: Multiverse.cxx:450
float abs(float number)
Definition: d0nt_math.hpp:39
_Var< caf::SRProxy > Var
Representation of a variable to be retrieved from a caf::StandardRecord object.
Definition: Var.h:74
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
const NuTruthCut kTrueDetector_NT([](const caf::SRNeutrinoProxy *nu){bool fid=(nu->vtx.X()< 192 && nu->vtx.X() >-191 && nu->vtx.Y()< 194 && nu->vtx.Y() >-187 && nu->vtx.Z()< 1270 && nu->vtx.Z() > 0);return fid; })
#define pot
_Cut< caf::SRProxy > Cut
Representation of a cut (selection) to be applied to a caf::StandardRecord object.
Definition: Cut.h:96
Var weights
caf::StandardRecord * sr
virtual void Go() override
Load all the registered spectra.
void PlotDebug(std::vector< TH1 * > universes, TH1 *h_nominal, std::string title, std::string name)
void multiverse_macro()
double frac(double x)
Fractional part.
UpDownPair< TH1 > ToUpDownHist(Multiverse *mv, const TH1 *h_nominal)
loader
Definition: demo0.py:10
std::vector< float > Spectrum
Definition: Constants.h:759
double POT() const
Definition: Spectrum.h:227
const NuTruthCut kIsNumuCC1Pi_NT([](const caf::SRNeutrinoProxy *nu){int nPions=0;int nMuons=0;int nPrims=nu->prim.size();for(int prim_idx=0;prim_idx< nPrims;prim_idx++){auto &prim=nu->prim[prim_idx];int pdg=prim.pdg;if(abs(pdg)==211)++nPions;else if(pdg==13)++nMuons;}bool is_numucc=nu->iscc &&nu->pdg==14;bool ret=is_numucc &&(nPions==1)&&(nMuons==1);return ret;})
const SystShifts kNoShift
Definition: SystShifts.cxx:22
const Cut kSignal
Definition: SINCpi0_Cuts.h:400
const Var kVtxY
Collaborates with Spectrum and OscillatableSpectrum to fill spectra from CAF files.
const Var kVtxZ
Var kUnity([](const caf::SRProxy *sr){return 1.f;})
static Binning Simple(int n, double lo, double hi, const std::vector< std::string > &labels={})
Definition: Binning.cxx:107
c cd(1)
Cut CutFromNuTruthCut(const NuTruthCut &stc)
Definition: Cut.cxx:7
enum BeamMode string
unsigned int uint