Functions
resolution2018.C File Reference
#include "CAFAna/Core/Loaders.h"
#include "CAFAna/Analysis/Prod4Loaders.h"
#include "CAFAna/Prediction/PredictionGenerator.h"
#include "CAFAna/Decomp/CheatDecomp.h"
#include "CAFAna/Core/Binning.h"
#include "CAFAna/Cuts/Cuts.h"
#include "3FlavorAna/Cuts/QuantileCuts.h"
#include "CAFAna/Analysis/Exposures.h"
#include "3FlavorAna/Systs/EnergySysts.h"
#include "CAFAna/Cuts/NueCutsSecondAna.h"
#include "3FlavorAna/Cuts/NumuCuts.h"
#include "CAFAna/Cuts/SpillCuts.h"
#include "3FlavorAna/Cuts/BeamNueCuts.h"
#include "CAFAna/Core/SystShifts.h"
#include "CAFAna/Core/Spectrum.h"
#include "CAFAna/Core/SpectrumLoader.h"
#include "CAFAna/Vars/Vars.h"
#include "3FlavorAna/Vars/Binnings.h"
#include "3FlavorAna/Vars/HistAxes.h"
#include "CAFAna/Fit/Fit.h"
#include "CAFAna/Analysis/Calcs.h"
#include "CAFAna/Analysis/Style.h"
#include "CAFAna/Systs/Systs.h"
#include "CAFAna/Vars/GenieWeights.h"
#include "CAFAna/Vars/PPFXWeights.h"
#include "CAFAna/Experiment/SingleSampleExperiment.h"
#include "CAFAna/Experiment/MultiExperiment.h"
#include "CAFAna/Prediction/PredictionScaleComp.h"
#include "CAFAna/Prediction/PredictionExtrap.h"
#include "CAFAna/Prediction/PredictionNoExtrap.h"
#include "CAFAna/Extrap/ModularExtrap.h"
#include "3FlavorAna/Decomp/BENDecomp.h"
#include "3FlavorAna/Decomp/MichelDecomp.h"
#include "3FlavorAna/Decomp/NumuDecomp.h"
#include "CAFAna/Decomp/ProportionalDecomp.h"
#include "TCanvas.h"
#include "TH1.h"
#include "TGaxis.h"
#include "TFile.h"
#include "TLegend.h"
#include "TLine.h"
#include "TStyle.h"
#include "TLatex.h"
#include "TROOT.h"
#include "TColor.h"
#include "THStack.h"
#include "TPaveText.h"
#include "TGraphErrors.h"
#include "TMultiGraph.h"
#include "Utilities/rootlogon.C"
#include "OscLib/OscCalcDumb.h"
#include "OscLib/IOscCalc.h"
#include <vector>

Go to the source code of this file.

Functions

void Beam (bool isRHC)
 
TCanvas * ResolutionPanels (TString cname, const std::vector< TString > labels, std::string epoch, bool, double, bool, bool)
 
TGraph * GetBiasErrorLine (TMultiGraph *multig, TString opt, int color)
 
void DrawResBin (PredictionExtrap *pred, osc::IOscCalc *calc, double fdpot, double center, TString sigBkg, int color, bool normalized, bool, TMultiGraph *=NULL, bool draw=true, bool tags=true)
 
void ResolutionBars (std::vector< TH1D * > allResol, const std::vector< int > idc, const std::vector< TString > idl, bool resid=false, bool isRHC=false)
 
void DifferenceOverlay (std::vector< TH1D * > allRes, const std::vector< int > idc, const std::vector< TString > idl, double maxy, bool isResolution=false, bool isRHC=false)
 
void resolution2018 (bool createFile=false, bool isRHC=false, std::string epoch="period3", TString suffix="test", bool unoscillated=true, bool useTrue=true)
 
void FormatResHist (TH1D *hist, TString xaxis, TString yaxis, bool choplabel, double xmin, double xmax)
 

Function Documentation

void Beam ( bool  isRHC)

Definition at line 82 of file resolution2018.C.

References calc, plot_validation_datamc_2018::color, DifferenceOverlay(), DrawResBin(), GetBiasErrorLine(), make_mec_shifts_plots::isRHC, PandAna.Demos.pi0_spectra::labels, maxy, MECModelEnuComparisons::opt, plot_validation_datamc::pred, prelim, ResolutionBars(), ResolutionPanels(), and msf_helper::tags.

Referenced by DifferenceOverlay(), GetNumuCosmicsFuture(), GetNumuPredictionsFuture(), plot_recoE_numu(), Plotting_DataAndPrediction(), ResolutionBars(), and ResolutionPanels().

82  {
83  if(isRHC){
84  TLatex* prelim = new TLatex(.276, .95, "Antineutrino Beam");
85  prelim->SetTextColor(kBlack);
86  prelim->SetNDC();
87  prelim->SetTextSize(1/30.);
88  prelim->SetTextAlign(32);
89  prelim->Draw();
90  }
91  else{
92  TLatex* prelim = new TLatex(.241, .95, "Neutrino Beam");
93  prelim->SetTextColor(kBlack);
94  prelim->SetNDC();
95  prelim->SetTextSize(1/30.);
96  prelim->SetTextAlign(32);
97  prelim->Draw();
98  }
99 }
TLatex * prelim
Definition: Xsec_final.C:133
void DifferenceOverlay ( std::vector< TH1D * >  allRes,
const std::vector< int idc,
const std::vector< TString >  idl,
double  maxy,
bool  isResolution = false,
bool  isRHC = false 
)

Definition at line 1054 of file resolution2018.C.

References Beam(), MECModelEnuComparisons::i, makeTrainCVSamples::int, MECModelEnuComparisons::leg, make_syst_table_plots::line, Simulation(), and ana::UniqueName().

Referenced by Beam(), and resolution2018().

1055  {
1056  gPad->SetFillStyle(0);
1057  TH1D * haux = new TH1D(UniqueName().c_str(),
1058  ";Energy difference (Reco - True) (GeV); A.U. (Area normalized)",100,
1059  -1,0.5);
1060  if(isResolution) haux->GetXaxis()->SetTitle("(Reco - True)/True");
1061  haux->GetYaxis()->SetRangeUser(0,maxy);
1062  haux->GetYaxis()->SetLabelSize(0);
1063  haux->GetXaxis()->CenterTitle();
1064  haux->GetYaxis()->CenterTitle();
1065  haux->DrawClone();
1066  TLegend* leg = new TLegend(0.12,0.2,0.4,0.89);
1067  leg->SetFillStyle(0);
1068  leg->SetTextSize(haux->GetXaxis()->GetTitleSize());
1069  //hack
1070 // if(allRes.size()>7) {
1071 // leg->SetTextAlign(31);
1072 // leg->SetX2(0.4);
1073 // leg->SetTextSize(0.9*haux->GetXaxis()->GetTitleSize());
1074 // }
1075  for (int i=0;i< (int)allRes.size();++i){
1076 // if(allRes[i]->GetMaximum()>maxy) maxy = allRes[i]->GetMaximum();
1077  allRes[i]->SetLineColor(idc[i]);
1078  allRes[i]->SetLineWidth(3);
1079 // allRes[i]->Rebin();
1080  allRes[i]->DrawNormalized("hist same");
1081  leg->AddEntry(allRes[i],idl[i],"l");
1082  }
1083 
1084 // haux->GetYaxis()->SetRangeUser(0,maxy);
1085  leg->SetHeader("FD MC");
1086  leg->Draw();
1087 // gPad->SetGridx();
1088  TLine *line = new TLine(0, 0, 0, maxy);
1089  line->SetVertical();
1090  line->SetLineColor(13);
1091  line->SetLineStyle(7);
1092  line->Draw("same");
1093  gPad->RedrawAxis();
1094  Simulation(); Beam(isRHC);
1095  delete haux;
1096  }
void Simulation()
Definition: tools.h:16
double maxy
void Beam(bool isRHC)
std::string UniqueName()
Return a different string each time, for creating histograms.
Definition: Utilities.cxx:29
void DrawResBin ( PredictionExtrap pred,
osc::IOscCalc calc,
double  fdpot,
double  center,
TString  sigBkg,
int  color,
bool  normalized,
bool  unoscillated,
TMultiGraph *  graph = NULL,
bool  draw = true,
bool  tags = true 
)

Definition at line 972 of file resolution2018.C.

References bin, om::cerr, analysePickle::hist, chisquared::hpred, MECModelEnuComparisons::i, label, ana::UniqueName(), and submit_syst::y.

Referenced by Beam().

982  {
983  gPad->SetFillStyle(0);
984 
985  TH1D* hist = new TH1D (UniqueName().c_str(),"temp",280,-1,6);
986  TH1D* hpred = new TH1D (UniqueName().c_str(),"p",10,0,1);
987 
988  if(sigBkg=="tot") hpred = (*pred).Predict(calc).ToTH1(fdpot,color);
989  else if(sigBkg=="sig") hpred= ((*pred).ComponentCC(14,14).Oscillated(calc,14,14)+(*pred).ComponentCC(-14,-14).Oscillated(calc,-14,-14)).
990  ToTH1(fdpot,color);
991  else if(sigBkg=="bkg") hpred = ((*pred).Predict(calc)-(*pred).ComponentCC(14,14).Oscillated(calc,14,14)
992  -(*pred).ComponentCC(-14,-14).Oscillated(calc,-14,-14)).
993  ToTH1(fdpot,color);
994  else {std::cerr <<"imconfused\n";}
995 
996  hpred->ResetStats();
997 
998  if(hpred->Integral()>0.01)
999  {
1000  if(tags && draw){
1001 
1002  TLatex * ltx = new TLatex();
1003  ltx->SetTextSize(9);
1004  ltx->SetTextFont(43);
1005  ltx->SetTextSize(1.2*ltx->GetTextSize());
1006  double y, lower1,lower2;
1007  if((int)((center-0.125)/0.25)%2==0){
1008  y=0.61;
1009  lower1=0;
1010  lower2=-0.3;
1011  }
1012  else{
1013  y = 0.47;
1014  lower1=0.5;
1015  lower2=-0.3;
1016  }
1017  TString label;
1018  label.Form
1019  ("#color[%d]{#bf{#splitline{#lower[%.2f]{#mu %.2f}}{#lower[%.2f]{#sigma %.2f}}}}",
1020  color,lower1,hpred->GetMean(),lower2,hpred->GetRMS());
1021  ltx->DrawLatex(center-0.125, y,label);
1022  }
1023  if(graph){
1024  TGraphErrors * mygraph = new TGraphErrors();
1025  mygraph->SetPoint(mygraph->GetN(),center, hpred->GetMean()) ;
1026  mygraph->SetPointError(mygraph->GetN()-1, 0, hpred->GetRMS());
1027  mygraph->SetMarkerStyle(kFullCircle);
1028  mygraph->SetLineColor(color);
1029  mygraph->SetFillColor(color);
1030  mygraph->SetMarkerColor(color);
1031  mygraph->SetLineWidth(3);
1032  graph->Add(mygraph, "lp same");
1033  }
1034  }
1035 
1036  for(int i=1;i<=hpred->GetNbinsX();++i){
1037  int bin = hist->FindBin(hpred->GetBinCenter(i)+center);
1038  hist->SetBinContent(bin,hpred->GetBinContent(i));
1039  }
1040 
1041  hist->SetLineColor(color);
1042 // hist->SetFillColorAlpha(color,.10);
1043  hist->Rebin();
1044  if(draw && hpred->Integral()>0.001){
1045  if(normalized) hist->DrawNormalized("hist same");
1046  else hist->DrawCopy("hist same");
1047  }
1048 
1049  delete hist, hpred;
1050 }
OStream cerr
Definition: OStream.cxx:7
const char * label
float bin[41]
Definition: plottest35.C:14
std::string UniqueName()
Return a different string each time, for creating histograms.
Definition: Utilities.cxx:29
void FormatResHist ( TH1D *  hist,
TString  xaxis,
TString  yaxis,
bool  choplabel,
double  xmin,
double  xmax 
)

Definition at line 853 of file resolution2018.C.

Referenced by ResolutionPanels().

855  {
856  hist->SetNdivisions(124);
857  hist->GetXaxis()->SetTitle(xaxis);
858  hist->GetYaxis()->SetTitle(yaxis);
859  hist->GetYaxis()->CenterTitle();
860  hist->GetYaxis()->SetDecimals();
861 
862  if(choplabel) {
863  hist->GetXaxis()->SetLabelSize(0);
864  hist->GetYaxis()->SetLabelSize(0);
865  hist->GetXaxis()->SetTitleOffset(1.9);
866  }
867  else {
868  hist->GetXaxis()->SetLabelSize(15);
869  hist->GetYaxis()->SetLabelSize(15);
870  hist->GetXaxis()->SetTitleOffset(3);
871  }
872  hist->GetYaxis()->SetRangeUser(xmin,xmax);
873 
874 }
std::map< std::string, double > xmax
std::string yaxis
std::string xaxis
TGraph * GetBiasErrorLine ( TMultiGraph *  multig,
TString  opt,
int  color 
)

Definition at line 830 of file resolution2018.C.

References compare_h5_caf::idx, runNovaSAM::ret, submit_syst::x, and submit_syst::y.

Referenced by Beam().

830  {
831  TGraph * ret = new TGraph();
832  auto ngraphs = multig->GetListOfGraphs()->GetSize();
833 
834 
835  for (int idx=0; idx<ngraphs; ++idx){
836  auto gr = (TGraphErrors*) multig->GetListOfGraphs()->At(idx);
837  auto x = gr->GetX()[0];
838  auto y = gr->GetY()[0];
839  auto erry = gr->GetEY()[0];
840  if(opt=="bias") ret->SetPoint(ret->GetN(),x,y);
841  if(opt=="error") ret->SetPoint(ret->GetN(),x,erry);
842  }
843  ret->SetLineWidth(3);
844  ret->SetLineColor(color);
845  ret->SetMarkerStyle(kFullCircle);
846  ret->SetMarkerColor(color);
847 
848  return ret;
849 
850 }
void resolution2018 ( bool  createFile = false,
bool  isRHC = false,
std::string  epoch = "period3",
TString  suffix = "test",
bool  unoscillated = true,
bool  useTrue = true 
)

Definition at line 117 of file resolution2018.C.

References allInOneTrainingPlots::axis, osc::_NoOscillations< T >::Copy(), om::cout, cut, ana::DefaultOscCalc(), DifferenceOverlay(), allTimeWatchdog::endl, file, shutoffs::filename, ana::Loaders::Go(), MECModelEnuComparisons::i, PandAna.Demos.tute_pid_validation::ids, compare_h5_caf::idx, inFile, makeTrainCVSamples::int, make_mec_shifts_plots::isRHC, ana::kAna2018FHCPOT, ana::kAna2018RHCPOT, ana::kCCE, ana::Loaders::kFHC, ana::kHadEFracAxis, ana::kIsCoh, ana::kIsDIS, ana::kIsDytmanMEC, ana::kIsQE, ana::kIsRes, ana::kNeutronVisEScalePrimariesSyst2018, ana::kNoCut, ana::kNoShift, ana::kNumuCCOptimisedAxis, ana::kNumuConcat, ana::kNumuCutFD2018, ana::kNumuEnergyBinning, ana::kPPFXFluxCVWgt, ana::Loaders::kRHC, ana::kTrueE, ana::kXSecCVWgt2018, loaders, ana::PredictionExtrap::LoadFrom(), gen_hdf5record::names, part, ana::pnfs2xrootd(), ana::PredictionExtrap::Predict(), ana::QuantileCutsFromTH2(), QuantNames, ResolutionBars(), ana::ShiDef::shift, ana::Binning::Simple(), art::to_string(), and ana::Spectrum::ToTH1().

120 {
121 
122  TString filename = "pred_resol_2018_"+suffix+"_"+epoch+".root";
123 
124  const Var resVar = (kCCE - kTrueE)/kTrueE;
125  const HistAxis resAxis ("Energy Resolution (Reco - True) / True",
126  Binning::Simple(200,-2.5,2.5),resVar);
127 
128  const Var residVar = (kCCE - kTrueE);
129  const HistAxis residAxis ("Energy Difference (Reco - True)",
130  Binning::Simple(200,-2.5,2.5),residVar);
131 
132  const HistAxis trueAxis("True Neutrino Energy (GeV)",kNumuEnergyBinning ,kTrueE);
133 
134 // std::vector <Cut> binCuts;
135 // std::vector <Cut> binCuts_truth;
136 //
137 // for(int i=0;i<20;i++){
138 // double low = i*0.25;
139 // double hig = (i+1.)*0.25;
140 // Cut tempCut = ((kCCE >= low) && (kCCE <hig));
141 // binCuts.emplace_back(tempCut);
142 // Cut tempCut_truth = ((kTrueE >= low) && (kTrueE <hig));
143 // binCuts_truth.emplace_back(tempCut_truth);
144 // }
145 
146  std::vector <Cut> intModes = {kIsQE, kIsRes, kIsDIS, kIsCoh, kIsDytmanMEC,kNoCut};
147  std::vector <TString> intNames = {"QE", "Res", "DIS", "Coh", "DytmanMEC","allModes"};
148 
149  // --- Quantiles
150  std::string fdspecfile = "/pnfs/nova/persistent/analysis/numu/Ana2018/provisional/quantiles/quantiles__rhc_full__numu2018.root";
151  TFile* inFile = TFile::Open( pnfs2xrootd(fdspecfile.c_str()).c_str() );
152  TH2 *FDSpec2D = (TH2*)inFile->FindObjectAny( "FDSpec2D" ); //TString("FD_full_") + TString(sFHC_RHC) );
153  const int NHadEFracQuantiles = 4;
154  std::vector<Cut> HadEFracQuantCuts = QuantileCutsFromTH2(FDSpec2D, kNumuCCOptimisedAxis, kHadEFracAxis, NHadEFracQuantiles);
155  HadEFracQuantCuts.emplace_back( kNoCut );
156  std::vector<TString> QuantNames = { "Quant1", "Quant2", "Quant3", "Quant4", "AllQuant" };
157 
158 
159  struct ShiDef {
160  TString name;
162  };
163  std::vector <ShiDef> gshifts;
164 
165  for(auto idx:{+1,-1}){
166  TString sigs = (idx>0? "plus":"minus");
167  // for now, just adding in the neutron systematic
168 
169  gshifts.push_back({"NeutronVisEScalePrimariesSyst2018_"+sigs,
171 
172 
173 // gshifts.emplace_back({"MaCCQEshape_"+sigs,
174 // SystShifts( new GenieRwgtTableSyst(rwgt::fReweightMaCCQEshape), idx)});
175 // gshifts.emplace_back({"NormCCQE_"+sigs,
176 // SystShifts((ISyst*) new GenieRwgtTableSyst(rwgt::fReweightNormCCQE), idx)});
177 // gshifts.emplace_back({"MvNCRES_"+sigs,
178 // SystShifts((ISyst*) new GenieRwgtTableSyst(rwgt::fReweightMvNCRES), idx)});
179 // gshifts.emplace_back({"MaCCRES_"+sigs,
180 // SystShifts((ISyst*)new GenieRwgtTableSyst(rwgt::fReweightMaCCRES), idx)});
181 // gshifts.emplace_back({"MaNCRES_"+sigs,
182 // SystShifts(new GenieRwgtTableSyst(rwgt::fReweightMaNCRES), idx)});
183 // gshifts.emplace_back({"MvCCRES_"+sigs,
184 // SystShifts(new GenieRwgtTableSyst(rwgt::fReweightMvCCRES), idx)});
185 // gshifts.emplace_back({"CCQEPauliSupViaKF_"+sigs,
186 // SystShifts(new GenieRwgtTableSyst(rwgt::fReweightCCQEPauliSupViaKF), idx)});
187 //
188 // gshifts.emplace_back({"RPACCQEEnh_"+sigs,
189 // SystShifts(kRPACCQEEnhSyst2018, idx)});
190 // gshifts.emplace_back({"kRPACCQESupp_"+sigs,
191 // SystShifts(kRPACCQESuppSyst2018, idx)});
192 //
193 // gshifts.emplace_back({"MECq0Shape_"+sigs,
194 // SystShifts(kMECq0ShapeSyst2018, idx)});
195 // gshifts.emplace_back({"MECInitStateNPFrac_"+sigs,
196 // SystShifts(kMECInitStateNPFracSyst2018, idx)});
197 // gshifts.emplace_back({"MECEnuShape_"+sigs,
198 // SystShifts(kMECEnuShapeSyst2018, idx)});
199 //
200  }
201 
202 
203 
204  if(createFile){
205  ECAFType ltype = kNumuConcat;
207  if(isRHC) ftype = Loaders::kRHC;
208  Prod4NomLoaders loaders(ltype, ftype, epoch, epoch);
209 
210  std::vector <PredictionNoExtrap *> predsNXP;
211  std::vector <CheatDecomp *> decomps;
212 
213  std::vector <TString> names;
214 
215  if(!suffix.Contains("ndonly")){
216  for(auto axis:{kNumuCCOptimisedAxis,resAxis,residAxis,trueAxis}){
217  for(auto cut:intModes){
218  for(auto quant:HadEFracQuantCuts){
219  predsNXP.emplace_back(new PredictionNoExtrap (loaders, axis, kNumuCutFD2018 && cut && quant,
220  kNoShift,
222  }
223  }//end modes
224  for(auto shift:gshifts){
225  for(auto quant:HadEFracQuantCuts){
226  predsNXP.emplace_back(new PredictionNoExtrap (loaders, axis, kNumuCutFD2018 && quant,
227  shift.shift,
229  }
230  }//end shifts
231  }
232  }
233  /*
234  if(suffix.Contains("ndonly")){
235  for(auto axis:{kNumuCCOptimisedAxis,resAxis,residAxis,trueAxis}){
236  decomps.emplace_back(new CheatDecomp(loaders.GetLoader(caf::kNEARDET,Loaders::kMC),
237  axis,kNumuCutND2018,kNoShift,kPPFXFluxCVWgt*kXSecCVWgt2018
238  ));
239  }
240  }*/
241  for(TString name:{"numu_recoE_all","numu_resol_all",
242  "numu_resid_all","numu_trueE_all"}){
243  for(auto intN:intNames){
244  for(auto quantN:QuantNames){
245  names.emplace_back(name+"_" + intN + "_" + quantN);
246  }
247  }
248  for(auto shift:gshifts){
249  for(auto quantN:QuantNames){
250  names.emplace_back(name+"_" + shift.name + "_" + quantN);
251  }
252  }
253  }
254 
255 // for (int i=0;i<20;i++){
256 // if(!suffix.Contains("ndonly")){
257 // for(auto axis:{kNumuCCOptimisedAxis,resAxis,residAxis}){
258 // for(auto cut:intModes){
259 // for(auto quant:HadEFracQuantCuts){
260 // predsNXP.emplace_back(new PredictionNoExtrap(loaders, axis,
261 // kNumuCutFD2018 && binCuts[i] && cut && quant,
262 // kNoShift, kPPFXFluxCVWgt*kXSecCVWgt2018));
263 // }
264 // }
265 // for(auto shift:gshifts){
266 // for(auto quant:HadEFracQuantCuts){
267 // predsNXP.emplace_back(new PredictionNoExtrap(loaders, axis,
268 // kNumuCutFD2018 && binCuts[i] && quant,
269 // shift.shift, kPPFXFluxCVWgt*kXSecCVWgt2018));
270 // }
271 // }
272 // }
273 // for(auto axis:{trueAxis,resAxis,residAxis}){
274 // for(auto cut:intModes){
275 // for(auto quant:HadEFracQuantCuts){
276 // predsNXP.emplace_back(new PredictionNoExtrap(loaders, axis,
277 // kNumuCutFD2018 && binCuts_truth[i] && cut,
278 // kNoShift, kPPFXFluxCVWgt*kXSecCVWgt2018));
279 // }
280 // }
281 // for(auto shift:gshifts){
282 // for(auto quant:HadEFracQuantCuts){
283 // predsNXP.emplace_back(new PredictionNoExtrap(loaders, axis,
284 // kNumuCutFD2018 && binCuts_truth[i],
285 // shift.shift, kPPFXFluxCVWgt*kXSecCVWgt2018));
286 // }
287 // }
288 // }
289 // }
290 //
291 // /* if(suffix.Contains("ndonly")){
292 // for(auto axis:{kNumuCCOptimisedAxis,resAxis,residAxis}){
293 // decomps.emplace_back(new CheatDecomp(loaders.GetLoader(caf::kNEARDET,Loaders::kMC),
294 // axis,kNumuCutND2018 && binCuts[i],kNoShift,kPPFXFluxCVWgt*kXSecCVWgt2018
295 // ));
296 // }
297 // for(auto axis:{trueAxis,resAxis,residAxis}){
298 // decomps.emplace_back(new CheatDecomp(loaders.GetLoader(caf::kNEARDET,Loaders::kMC),
299 // axis,kNumuCutND2018 && binCuts_truth[i],kNoShift,kPPFXFluxCVWgt*kXSecCVWgt2018
300 // ));
301 // }
302 // }*/
303 //
304 // for(TString name:{"numu_recoE_bin_","numu_resol_bin_","numu_resid_bin_"}){
305 // for(auto intN:intNames){
306 // for(auto quantN:QuantNames){
307 // names.emplace_back(name+std::to_string(i)+"_" + intN + "_" + quantN);
308 // }
309 // }
310 // for(auto shift:gshifts){
311 // for(auto quantN:QuantNames){
312 // names.emplace_back(name+std::to_string(i)+"_" + shift.name + "_" + quantN);
313 // }
314 // }
315 // }
316 //
317 // for(TString name:{"numu_trueE_tbin_","numu_resol_tbin_","numu_resid_tbin_"}){
318 // for(auto intN:intNames){
319 // for(auto quantN:QuantNames){
320 // names.emplace_back(name+std::to_string(i)+"_" + intN + "_" + quantN);
321 // }
322 // }
323 // for(auto shift:gshifts){
324 // for(auto quantN:QuantNames){
325 // names.emplace_back(name+std::to_string(i)+"_" + shift.name + "_" + quantN);
326 // }
327 // }
328 // }
329 // }//end bins
330 
331  loaders.Go();
332 
333  TFile * file = new TFile(filename,"recreate");
334 
335 
336  for(int i=0;i<(int)names.size();++i){
337  if(!suffix.Contains("ndonly")){
338  predsNXP[i]->SaveTo(file, "pred_nxp_"+names[i]);
339  }
340 // if(suffix.Contains("extrap") || suffix.Contains("ndonly")){
341 // decomps[i]->SaveTo(file, "decomp_"+names[i]);
342  // }
343  }
344  return;
345  }
346 
347 
348  else{
349  //make some plots
350  TFile * file = new TFile(filename,"read");
351  double fdpot = kAna2018FHCPOT;
352  if(isRHC) double fdpot = kAna2018RHCPOT;
353 
355  auto calc = nos.Copy();
356  if(!unoscillated) calc = (osc::IOscCalc*) DefaultOscCalc();
357  TString oscStr ="_";
358  if(unoscillated) oscStr = "_noosc_";
359  if(useTrue) oscStr="_true" + oscStr;
360 
361  TString pstr = "pred_nxp_numu_";
362 
363  for(int qIdx = 0; qIdx < QuantNames.size(); ++qIdx){
364  auto pred = PredictionExtrap::LoadFrom(file,
365  pstr+(useTrue?"trueE":"recoE") + "_all_allModes_"+QuantNames[qIdx]);
366  auto hp = pred->Predict(calc).ToTH1(fdpot);
367 
368  bool normalized = true;
369  std::vector <TString> ids = {"sig","bkg","tot"};
370  std::vector <int> idc = {kRed,kBlue,kBlack};
371  std::vector <TString> idl = {"#nu_{#mu} CC","Bkgd.","Total"};
372 
373  //all bins
374  auto c1_h = new TCanvas();
375  TString ps_all = pstr + "resid_all_allModes_" + QuantNames[qIdx];
376  auto resid = PredictionExtrap::LoadFrom(file, ps_all);
377  std::vector <TH1D*> allRes;
378  allRes.emplace_back((resid->ComponentCC(14,14).Oscillated(calc,14,14) +
379  resid->ComponentCC(-14,-14).Oscillated(calc,-14,-14)).ToTH1(fdpot));
380  allRes.emplace_back((resid->Predict(calc) -
381  resid->ComponentCC(14,14).Oscillated(calc,14,14) - resid->ComponentCC(-14,-14).Oscillated(calc,-14,-14))
382  .ToTH1(fdpot));
383  allRes.emplace_back(resid->Predict(calc).ToTH1(fdpot));
384 
385  DifferenceOverlay(allRes,idc,idl, 0.1, false, isRHC);
386 
387  auto c1_hu = new TCanvas();
388  std::vector <TH1D*> allResu = {(resid->ComponentCC(14,14).Oscillated(calc,14,14)+resid->ComponentCC(-14,-14).Oscillated(calc,-14,-14)).ToTH1(fdpot)};
389  DifferenceOverlay(allResu,{idc[0]},{idl[0]},0.15, false, isRHC);
390 
391 
392  ps_all = pstr + "resol_all_allModes_"+QuantNames[qIdx];
393  auto resol = PredictionExtrap::LoadFrom(file, ps_all);
394  std::vector <TH1D*> allResol;
395  allResol.emplace_back((resol->ComponentCC(14,14).Oscillated(calc,14,14) +
396  resol->ComponentCC(-14,-14).Oscillated(calc,-14,-14)).ToTH1(fdpot));
397  allResol.emplace_back((resol->Predict(calc) -
398  resol->ComponentCC(14,14).Oscillated(calc,14,14) - resol->ComponentCC(-14,-14).Oscillated(calc,-14,-14))
399  .ToTH1(fdpot));
400  allResol.emplace_back(resol->Predict(calc).ToTH1(fdpot));
401 
402  auto c1_h2 = new TCanvas;
403  DifferenceOverlay(allResol,idc,idl, 0.2, true, isRHC);
404 
405  auto c1_h2u = new TCanvas();
406  std::vector <TH1D*> allResolu = {(resol->ComponentCC(14,14).Oscillated(calc,14,14) +
407  resol->ComponentCC(-14,-14).Oscillated(calc,-14,-14)).ToTH1(fdpot)};
408  DifferenceOverlay(allResolu,{idc[0]},{idl[0]},0.1);
409 
410  auto c1_mg = new TCanvas();
411  ResolutionBars(allResol, idc,idl, false, isRHC);
412 
413  auto c1_mg2 = new TCanvas();
414  ResolutionBars(allRes, idc,idl, true, isRHC);
415 
416 
417  //all bins
418  //repeat for numucc mec
419  //WTF Erika
420  // std::vector <TString> oldN = {"DIS","Res","QE","DytmanMEC"};
421  std::vector <TString> oldN = {"All","DytmanMEC","QE","Res","DIS"};
422  std::vector <TString> interN = {"All", "MEC","QE","RES","DIS"};
423  std::vector <int> interC = {kBlack, kOrange-1,kAzure+7, kGreen+2,kGray+2};
424 
425  std::vector <TString> cuteLabel;
426  for (int i =0; i<5; ++i){
427  TString tstr; tstr.Form("#color[%d]{%s}",interC[i],interN[i].Data());
428  cuteLabel.emplace_back(tstr);
429  }
430 
431  std::vector <TH1D*> allResM;
432  std::vector <TH1D*> allResolM;
433 
434  allResM.emplace_back((resid->ComponentCC(14,14).Oscillated(calc,14,14)+resid->ComponentCC(-14,-14).Oscillated(calc,-14,-14)).ToTH1(fdpot));
435  allResolM.emplace_back((resol->ComponentCC(14,14).Oscillated(calc,14,14)+resol->ComponentCC(-14,-14).Oscillated(calc,-14,-14)).ToTH1(fdpot));
436  for(int i=1;i<5;++i){
437  ps_all = pstr + "resid_all_"+oldN[i]+"_"+QuantNames[qIdx];
438  auto residM = PredictionExtrap::LoadFrom(file, ps_all);
439  allResM.emplace_back((residM->ComponentCC(14,14).Oscillated(calc,14,14)+residM->ComponentCC(-14,-14).Oscillated(calc,-14,-14)).ToTH1(fdpot));
440 
441  ps_all = pstr + "resol_all_"+oldN[i]+"_"+QuantNames[qIdx];
442  auto resolM = PredictionExtrap::LoadFrom(file, ps_all);
443  allResolM.emplace_back((resolM->ComponentCC(14,14).Oscillated(calc,14,14)+resolM->ComponentCC(-14,-14).Oscillated(calc,-14,-14)).ToTH1(fdpot));
444  }
445 
446  auto c2_h = new TCanvas();
447  DifferenceOverlay(allResM,interC,{"All","MEC","QE","RES","DIS"},0.15, false, isRHC);
448 
449  auto c2_h2 = new TCanvas;
450  DifferenceOverlay(allResolM,interC,interN, 0.25, true, isRHC);
451 
452  auto c2_mg = new TCanvas();
453  ResolutionBars(allResolM, interC,{"All","MEC","QE","RES","DIS"}, false, isRHC);
454  auto c2_mg2 = new TCanvas();
455  ResolutionBars(allResM, interC,{"All", "MEC","QE","RES","DIS"}, true, isRHC);
456 
457 
458 
459 
460  //////////////////////////////////////////////////////////////////////////////////
461  //repeat for genie
462  //all bins
463  //repeat for numucc mec
464  std::vector<TString> gOld = {"NeutronVisEScalePrimariesSyst2018_"};
465 
466  std::vector<TString> gShort = {"NeutronSyst"};
467 
468  // std::vector <TString> gOld = {"NormCCQE_",
469  // "MaCCQEshape_",
470  // "CCQEPauliSupViaKF_",
471  // "RPACCQE_",
472  // "MECScale_",
473  // "MaCCRES_",
474  // "MaNCRES_",
475  // "MvCCRES_",
476  // "MvNCRES_" };
477  //
478  // std::vector <TString> gShort= {"Norm^{CCQE}",
479  // "M_{A}^{CCQE,sc}",
480  // "Fermi^{CCQE}",
481  // "RPA^{CCQE}",
482  // "MEC^{sc}",
483  // "M_{A}^{CCRES}",
484  // "M_{A}^{NCRES}",
485  // "M_{v}^{CCRES}",
486  // "M_{v}^{NCRES}"
487  // };
488  //
489  //
490  //
491  for (auto part:{1,1}){
492  int ming = 0; int maxg = 1;
493  // if(part==2) {ming=5; maxg=9;}
494 
495  std::vector <TString> gOldSig;
496  std::vector <TString> gShortSig;
497  std::vector <int> gColors;
498 
499  for(int i=ming;i<maxg;i++){
500  gOldSig.emplace_back(gOld[i]+"plus");
501  gShortSig.emplace_back(gShort[i]+" (#plus1#sigma)");
502  gColors.emplace_back(kPink -i);
503 
504  gOldSig.emplace_back("Nominal");
505  gShortSig.emplace_back("Nominal");
506  gColors.emplace_back(kBlack);
507 
508  gOldSig.emplace_back(gOld[i]+"minus");
509  gShortSig.emplace_back(gShort[i]+" (#minus1#sigma)");
510  gColors.emplace_back(kAzure -i);
511  }
512 
513  std::vector <TH1D*> allResG;
514  std::vector <TH1D*> allResolG;
515 
516  // this is a hack for the neutron systematic only
517 // for(int i=0;i<(int)gOldSig.size();i++){
518  for(int i=0;i<1;i++){
519 
520  ps_all = pstr + "resid_all_"+gOldSig[i]+"_"+QuantNames[qIdx];
521  std::cout <<ps_all << std::endl;
522  auto residM = PredictionExtrap::LoadFrom(file, ps_all);
523  allResG.emplace_back((residM->ComponentCC(14,14).Oscillated(calc,14,14)+residM->ComponentCC(-14,-14).Oscillated(calc,-14,-14)).ToTH1(fdpot));
524  ps_all = pstr + "resol_all_"+gOldSig[i]+"_"+QuantNames[qIdx];
525  auto resolM = PredictionExtrap::LoadFrom(file, ps_all);
526  allResolG.emplace_back((resolM->ComponentCC(14,14).Oscillated(calc,14,14)+resolM->ComponentCC(-14,-14).Oscillated(calc,-14,-14)).ToTH1(fdpot));
527  }
528 
529  // Nominal
530  allResG.emplace_back((resid->ComponentCC(14,14).Oscillated(calc,14,14) +
531  resid->ComponentCC(-14,-14).Oscillated(calc,-14,-14)).ToTH1(fdpot));
532  allResolG.emplace_back((resol->ComponentCC(14,14).Oscillated(calc,14,14) +
533  resol->ComponentCC(-14,-14).Oscillated(calc,-14,-14)).ToTH1(fdpot));
534 
535  for(int i=2;i<3;i++){
536 
537  ps_all = pstr + "resid_all_"+gOldSig[i]+"_"+QuantNames[qIdx];
538  std::cout <<ps_all << std::endl;
539  auto residM = PredictionExtrap::LoadFrom(file, ps_all);
540  allResG.emplace_back((residM->ComponentCC(14,14).Oscillated(calc,14,14)+residM->ComponentCC(-14,-14).Oscillated(calc,-14,-14)).ToTH1(fdpot));
541  ps_all = pstr + "resol_all_"+gOldSig[i]+"_"+QuantNames[qIdx];
542  auto resolM = PredictionExtrap::LoadFrom(file, ps_all);
543  allResolG.emplace_back((resolM->ComponentCC(14,14).Oscillated(calc,14,14)+resolM->ComponentCC(-14,-14).Oscillated(calc,-14,-14)).ToTH1(fdpot));
544  }
545 
546 
547 // }
548 
549  auto c3_h = new TCanvas();
550  DifferenceOverlay(allResG,gColors, gShortSig,0.1, false, isRHC);
551 
552  auto c3_h2 = new TCanvas();
553  DifferenceOverlay(allResolG,gColors, gShortSig,0.2,true, isRHC);
554 
555  auto c3_mg = new TCanvas();
556  ResolutionBars(allResolG, gColors, gShortSig, false, isRHC);
557  auto c3_mg2 = new TCanvas();
558  ResolutionBars(allResG, gColors, gShortSig, true, isRHC);
559 
560  c3_mg->Print("plot4_resol_2018_"+suffix+"_"+epoch + oscStr +"bars_xsec"+std::to_string(part)+"_"+QuantNames[qIdx]+".pdf");
561  c3_mg2->Print("plot4_resol_2018_"+suffix+"_"+epoch + oscStr +"bars2_xsec"+std::to_string(part)+"_"+QuantNames[qIdx]+".pdf");
562  c3_h->Print("plot4_resol_2018_"+suffix+"_"+epoch + oscStr +"overlay_xsec"+std::to_string(part)+"_"+QuantNames[qIdx]+".pdf");
563  c3_h2->Print("plot4_resol_2018_"+suffix+"_"+epoch + oscStr +"overlayf_xsec"+std::to_string(part)+"_"+QuantNames[qIdx]+".pdf");
564  }
565  // // return;
566  //
567  // ////Need to do those before gstyle change im an idiot
568  //
569  //
570  // for (auto part:{1,1}){
571  // int ming = 0; int maxg = 1;
572  // if(part==2) {ming=5; maxg=9;}
573  //
574  // std::vector <TString> gOldSig;
575  // std::vector <TString> gShortSig;
576  // std::vector <int> gColors;
577  //
578  // for(int i=ming;i<maxg;i++){
579  // gOldSig.emplace_back(gOld[i]+"plus");
580  // gShortSig.emplace_back(gShort[i]+" (#plus1#sigma)");
581  // gColors.emplace_back(kPink -i);
582  //
583  // gOldSig.emplace_back(gOld[i]+"minus");
584  // gShortSig.emplace_back(gShort[i]+" (#minus1#sigma)");
585  // gColors.emplace_back(kAzure -i);
586  // }
587  //
588  // auto c3_f = ResolutionPanels("genie", gShortSig, epoch, useTrue,
589  // 1.2*hp->GetMaximum(),true);
590  //
591  //
592  // TLegend * legG = new TLegend(0.08,0.58,0.95,0.84);
593  // legG->SetFillStyle(0);
594  // legG->SetNColumns(4);
595  // legG->SetTextSize(0.06);
596  //
597  // for(int idx = 0; idx <gShortSig.size(); ++idx)
598  // {
599  // TMultiGraph* gr_resol = new TMultiGraph();
600  //
601  // for(int i =0; i<20;i++){
602  // TString ps_bin = (TString)pstr+(useTrue?+"resol_tbin_":"resol_bin_")
603  // + std::to_string(i)+"_";
604  // auto mbin = PredictionExtrap::LoadFrom(file, ps_bin+gOldSig[idx]);
605  //
606  // DrawResBin(mbin.release(),calc,fdpot,(double)i*0.25+.125,"sig",kBlack,
607  // normalized, unoscillated,gr_resol,false, false);
608  // }
609  //
610  // //Transfotm pretty multigraph into boring lines
611  // c3_f->cd(2);
612  // auto bias = GetBiasErrorLine(gr_resol,"bias",gColors[idx]);
613  // bias->Draw("Lp same");
614  // legG->AddEntry(bias,gShortSig[idx],"lp");
615  // gPad->RedrawAxis();
616  // c3_f->cd(1);
617  // GetBiasErrorLine(gr_resol,"error",gColors[idx])->Draw("lp same");
618  // gPad->RedrawAxis();
619  // }
620  //
621  // c3_f->cd(1);
622  // legG->Draw();
623  //
624  // c3_f->cd(2);
625  // legG->Draw();
626  //
627  //
628  // c3_f->Print("plot4_resol_2018_"+suffix+"_"+epoch + oscStr +"fracs_xsec"+
629  // std::to_string(part)+"_"+QuantNames[qIdx]+".pdf");
630  //
631  // }//end parts genie
632 
633 
634 
635  // auto c1 = ResolutionPanels("total", {"#color[2]{#nu_{#mu}CC}","#color[4]{Bkgd.}","#color[1]{Total}"}, epoch,useTrue, 1.2*hp->GetMaximum(),false);
636  // auto c1_f = ResolutionPanels("total", {"True #nu_{#mu} CC signal","Background","Total (Sig. + Bkg.)"}, epoch, useTrue, 1.2*hp->GetMaximum(),true);
637  // auto c1_fu = ResolutionPanels("numu", {"True #nu_{#mu} CC signal"},epoch, useTrue, 1.2*hp->GetMaximum(),true);
638  //
639  //
640  // TLegend * leg3 = new TLegend(0.08,0.74,0.68,0.84);
641  // leg3->SetFillStyle(0);
642  // leg3->SetNColumns(3);
643  //
644  // int mycolors[]={kPink-5,kSpring-7,kOrange+7,kAzure+2,kPink+10,/*kTeal-5,*/kViolet-6,kCyan+3};
645  // auto bkgcolor = kGray+3;
646  //
647  // //Draw the spectrum
648  // //Pretty bin colors
649  // c1->cd(1);
650  // for(int i =0; i<20;i++){
651  // auto pbin = PredictionExtrap::LoadFrom(
652  // file->GetDirectory((TString)pstr+(useTrue?+"trueE_tbin_":"recoE_bin_")+ std::to_string(i)+"_allModes"));
653  // int color = mycolors[i%7];
654  // auto htot_bin = pbin->Predict(calc).ToTH1(fdpot,color);
655  // auto hbkg_bin = (pbin->Predict(calc) - pbin->ComponentCC(14,14).Oscillated(calc,14,14)).ToTH1(fdpot,bkgcolor);
656  // htot_bin->SetFillColorAlpha(color,.1);
657  // hbkg_bin->SetFillColor(bkgcolor);
658  // htot_bin->Draw("hist same");
659  // hbkg_bin->Draw("hist f same");
660  // double offset = 0.1 * hp->GetMaximum();
661  // TLatex * ltx = new TLatex();
662  // ltx->SetTextAlign(22);
663  // ltx->SetTextSize(1.2*ltx->GetTextSize());
664  // TString binStr; binStr.Form("#color[%d]{#bf{B%d}}",color,i+1);
665  // ltx->DrawLatex((double)i*.25+0.125, htot_bin->Integral()+offset,binStr);
666  // }
667  //
668  // for(int idx = 0; idx <3; ++idx)
669  // {
670  // TMultiGraph* gr_resid = new TMultiGraph();
671  // TMultiGraph* gr_resol = new TMultiGraph();
672  //
673  // for(int i =0; i<20;i++){
674  // int color = mycolors[i%7];
675  // //True - reco graphs
676  // auto pbin = PredictionExtrap::LoadFrom(
677  // file->GetDirectory((TString)pstr+(useTrue?+"resid_tbin_":"resid_bin_")+ std::to_string(i)+"_allModes")).release();
678  // c1->cd(idx+2);
679  // DrawResBin(pbin,calc,fdpot,(double)i*0.25+.125,ids[idx],color,normalized,unoscillated,gr_resid, true, true);
680  //
681  // //Get the resol info; false since we are not drawing these anymore
682  // pbin = PredictionExtrap::LoadFrom(
683  // file->GetDirectory((TString)pstr+(useTrue?+"resol_tbin_":"resol_bin_")+ std::to_string(i)+"_allModes")).release();
684  // DrawResBin(pbin,calc,fdpot,(double)i*0.25+.125,ids[idx],color,normalized,unoscillated,gr_resol,false, false);
685  // gPad->RedrawAxis();
686  // }
687  //
688  // //Transfotm pretty multigraph into boring lines
689  // c1_f->cd(2);
690  // auto bias = GetBiasErrorLine(gr_resol,"bias",idc[idx]);
691  // bias->Draw("Lp same");
692  // leg3->AddEntry(bias,idl[idx],"lp");
693  // gPad->RedrawAxis();
694  // c1_f->cd(1);
695  // GetBiasErrorLine(gr_resol,"error",idc[idx])->Draw("lp same");
696  // gPad->RedrawAxis();
697  //
698  // if(idx==0){
699  // c1_fu->cd(1);
700  // GetBiasErrorLine(gr_resol,"error",idc[idx])->Draw("lp same");
701  // leg3->Draw();
702  // gPad->RedrawAxis();
703  // c1_fu->cd(2);
704  // GetBiasErrorLine(gr_resol,"bias",idc[idx])->Draw("lp same");
705  // leg3->Draw();
706  // gPad->RedrawAxis();
707  // c1_fu->Print("plot4_resol_2018_"+suffix+"_"+epoch + oscStr +"fracs_numu"+"_"+QuantNames[qIdx]+.pdf");
708  // }
709  // }
710  //
711  // c1_f->cd(1);
712  // leg3->Draw();
713  //
714  // c1_f->cd(2);
715  // leg3->Draw();
716  //
717  // c1->Print("plot4_resol_2018_"+suffix+"_"+epoch + oscStr +"bins"+"_"+QuantNames[qIdx]+".pdf");
718  // c1_f->Print("plot4_resol_2018_"+suffix+"_"+epoch + oscStr +"fracs"+"_"+QuantNames[qIdx]+".pdf");
719  c1_mg->Print("plot4_resol_2018_"+suffix+"_"+epoch + oscStr +"bars"+"_"+QuantNames[qIdx]+".pdf");
720  c1_mg2->Print("plot4_resol_2018_"+suffix+"_"+epoch + oscStr +"bars2"+"_"+QuantNames[qIdx]+".pdf");
721  c1_h->Print("plot4_resol_2018_"+suffix+"_"+epoch + oscStr +"overlay"+"_"+QuantNames[qIdx]+".pdf");
722  c1_h2->Print("plot4_resol_2018_"+suffix+"_"+epoch + oscStr +"overlayf"+"_"+QuantNames[qIdx]+".pdf");
723  c1_hu->Print("plot4_resol_2018_"+suffix+"_"+epoch + oscStr +"overlay_numu"+"_"+QuantNames[qIdx]+".pdf");
724  c1_h2u->Print("plot4_resol_2018_"+suffix+"_"+epoch + oscStr +"overlayf_numu"+"_"+QuantNames[qIdx]+".pdf");
725 
726 
727  ////////////////////////////////////////////////////////////////////////////////
728 
729  // auto c2 = ResolutionPanels("mode", cuteLabel,epoch,useTrue, 1.2*hp->GetMaximum(),false);
730  // auto c2_f = ResolutionPanels("mode", interN,epoch,useTrue, 1.2*hp->GetMaximum(),true);
731 
732 
733  // c2->cd(1);
734  // for(int i =0; i<20;i++){
735  // TString ps_bin = (TString)pstr+(useTrue?+"trueE_tbin_":"recoE_bin_")
736  // + std::to_string(i)+"_";
737  // auto pbin = PredictionExtrap::LoadFrom(file, ps_bin+"allModes");
738  // auto cbin = PredictionExtrap::LoadFrom(file, ps_bin+"Coh");
739  //
740  // THStack *hs = new THStack("hs","");
741  // //Other is bkg + coh
742  // int color = mycolors[i%7];
743  // auto htot_bin = pbin->Predict(calc).ToTH1(fdpot,/*color*/kBlack);
744  // auto hoth_bin = (pbin->Predict(calc) -
745  // pbin->ComponentCC(14,14).Oscillated(calc,14,14) +
746  // cbin->ComponentCC(14,14).Oscillated(calc,14,14))
747  // .ToTH1(fdpot,bkgcolor);
748  // hoth_bin->SetFillColorAlpha(bkgcolor, 1);
749  // hoth_bin->SetLineColorAlpha(bkgcolor, 0);
750  // hs->Add(hoth_bin,"hist f");
751  //
752  // for (int idx=3;idx>=0;--idx){
753  // auto mbin = PredictionExtrap::LoadFrom(file, ps_bin+oldN[idx]);
754  // auto hmod_bin = mbin->ComponentCC(14,14).Oscillated(calc,14,14).ToTH1(fdpot,interC[idx]);
755  // hmod_bin->SetFillColorAlpha(interC[idx], 0.3);
756  // hmod_bin->SetLineColorAlpha(interC[idx], 0);
757  // hs->Add(hmod_bin,"hist f");
758  // }
759  // hs->Draw("f same");
760  // htot_bin->Draw("hist same");
761  //
762  // double offset = 0.1 * hp->GetMaximum();
763  // TLatex * ltx = new TLatex();
764  // ltx->SetTextAlign(22);
765  // ltx->SetTextSize(1.2*ltx->GetTextSize());
766  // TString binStr; binStr.Form("#color[%d]{#bf{B%d}}",color,i+1);
767  // ltx->DrawLatex((double)i*.25+0.125, htot_bin->Integral()+offset,binStr);
768  // }//end bins
769  //
770  // TLegend * legM = new TLegend(0.08,0.74,0.68,0.84);
771  // legM->SetFillStyle(0);
772  // legM->SetNColumns(4);
773  //
774  // for(int idx = 0; idx <4; ++idx)
775  // {
776  // TMultiGraph* gr_resid = new TMultiGraph();
777  // TMultiGraph* gr_resol = new TMultiGraph();
778  //
779  // for(int i =0; i<20;i++){
780  // //True - reco graphs
781  // int color = mycolors[i%7];
782  // TString ps_bin = (TString)pstr+(useTrue?+"resid_tbin_":"resid_bin_")
783  // + std::to_string(i)+"_";
784  // auto mbin = PredictionExtrap::LoadFrom(file, ps_bin+oldN[idx]);
785  // c2->cd(idx+2);
786  // DrawResBin(mbin.release(),calc,fdpot,(double)i*0.25+.125,"sig",color,
787  // normalized,unoscillated,gr_resid, true, true);
788  // gPad->RedrawAxis();
789  //
790  // ps_bin = (TString)pstr+(useTrue?+"resol_tbin_":"resol_bin_")
791  // + std::to_string(i)+"_";
792  // mbin = PredictionExtrap::LoadFrom(file, ps_bin+oldN[idx]);
793  //
794  // DrawResBin(mbin.release(),calc,fdpot,(double)i*0.25+.125,"sig",color,
795  // normalized, unoscillated,gr_resol,false, false);
796  // }
797  //
798  // //Transfotm pretty multigraph into boring lines
799  // c2_f->cd(2);
800  // auto bias = GetBiasErrorLine(gr_resol,"bias",interC[idx]);
801  // bias->Draw("Lp same");
802  // legM->AddEntry(bias,interN[idx],"lp");
803  // gPad->RedrawAxis();
804  // c2_f->cd(1);
805  // GetBiasErrorLine(gr_resol,"error",interC[idx])->Draw("lp same");
806  // gPad->RedrawAxis();
807  // }
808  //
809  // c2_f->cd(1);
810  // legM->Draw();
811  //
812  // c2_f->cd(2);
813  // legM->Draw();
814  //
815  //
816  //
817  // c2->Print("plot4_resol_2018_"+suffix+"_"+epoch + oscStr +"bins_mode"+"_"+QuantNames[qIdx]+".pdf");
818  // c2_f->Print("plot4_resol_2018_"+suffix+"_"+epoch + oscStr +"fracs_mode"+"_"+QuantNames[qIdx]+".pdf");
819  c2_mg->Print("plot4_resol_2018_"+suffix+"_"+epoch + oscStr +"bars_mode"+"_"+QuantNames[qIdx]+".pdf");
820  c2_mg2->Print("plot4_resol_2018_"+suffix+"_"+epoch + oscStr +"bars2_mode"+"_"+QuantNames[qIdx]+".pdf");
821  c2_h->Print("plot4_resol_2018_"+suffix+"_"+epoch + oscStr +"overlay_mode"+"_"+QuantNames[qIdx]+".pdf");
822  c2_h2->Print("plot4_resol_2018_"+suffix+"_"+epoch + oscStr +"overlayf_mode"+"_"+QuantNames[qIdx]+".pdf");
823  } // end quant
824  }//end plots
825 }
const Cut kIsQE
Definition: TruthCuts.cxx:104
const XML_Char * name
Definition: expat.h:151
Pass neutrinos through unchanged.
Definition: IOscCalc.h:40
const Var kPPFXFluxCVWgt
Definition: PPFXWeights.h:16
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
const Cut kIsRes
Definition: TruthCuts.cxx:111
General interface to oscillation calculators.
Definition: StanTypedefs.h:22
const HistAxis kHadEFracAxis("E_{had.} / E_{#nu}", Binning::Simple(200, 0, 1), kHadEFrac)
HistAxis that implements Hadronic Energy fraction binning used by L.Vinton to derive Hadronic Energy ...
Definition: HistAxes.h:30
void ResolutionBars(std::vector< TH1D * > allResol, const std::vector< int > idc, const std::vector< TString > idl, bool resid=false, bool isRHC=false)
const Cut kIsDIS
Definition: TruthCuts.cxx:118
string filename
Definition: shutoffs.py:106
std::string pnfs2xrootd(std::string loc, bool unauth)
Definition: UtilsExt.cxx:237
const NeutronVisEScaleSyst2018 kNeutronVisEScalePrimariesSyst2018(false)
osc::OscCalcDumb calc
osc::IOscCalcAdjustable * DefaultOscCalc()
Create a new calculator with default assumptions for all parameters.
Definition: Calcs.cxx:49
std::unique_ptr< T > LoadFrom(TDirectory *dir, const std::string &label)
Definition: LoadFromFile.h:17
const double kAna2018RHCPOT
Definition: Exposures.h:208
const Cut kIsDytmanMEC
Definition: TruthCuts.cxx:187
ifstream inFile
Definition: AnaPlotMaker.h:34
TString part[npart]
Definition: Style.C:32
const Var kTrueE([](const caf::SRProxy *sr){return(sr->mc.nnu==0)?0.:float(sr->mc.nu[0].E);})
Definition: Vars.cxx:85
const Binning kNumuEnergyBinning
Definition: Binnings.cxx:13
virtual _IOscCalc< T > * Copy() const override
Definition: IOscCalc.h:44
const Var kCCE
Definition: NumuVars.h:21
const std::string QuantNames[NumQuant]
Definition: DoThePlots.C:207
const HistAxis kNumuCCOptimisedAxis("Reconstructed Neutrino Energy (GeV)", kNumuCCEOptimisedBinning, kCCE)
HistAxis that implements optimised numuCCE from L. Vinton. See docdb 16332. This was close to &#39;custC&#39;...
Definition: HistAxes.h:24
const SystShifts kNoShift
Definition: SystShifts.cxx:21
OStream cout
Definition: OStream.cxx:6
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
const Cut cut
Definition: exporter_fd.C:30
void DifferenceOverlay(std::vector< TH1D * > allRes, const std::vector< int > idc, const std::vector< TString > idl, double maxy, bool isResolution=false, bool isRHC=false)
const Var kXSecCVWgt2018
Definition: XsecTunes.h:49
std::vector< Loaders * > loaders
Definition: syst_header.h:386
std::vector< Cut > QuantileCutsFromTH2(TH2 *quantileHist, const HistAxis &independentAxis, const HistAxis &quantileAxis, const unsigned int &numQuantiles, const bool verbose)
: Do the same as the QuantileCuts function but taking in the TH2 instead of making it...
TFile * file
Definition: cellShifts.C:17
const Cut kNumuCutFD2018
Definition: NumuCuts2018.h:39
Prediction that just uses FD MC, with no extrapolation.
const double kAna2018FHCPOT
Definition: Exposures.h:207
std::string to_string(ModuleType mt)
Definition: ModuleType.h:32
For nominal spectra and reweighting systs (xsec/flux)
Definition: Prod4Loaders.h:96
SystShifts shift
const Cut kIsCoh
Definition: TruthCuts.cxx:133
ECAFType
Definition: Loaders.h:19
void ResolutionBars ( std::vector< TH1D * >  allResol,
const std::vector< int idc,
const std::vector< TString >  idl,
bool  resid = false,
bool  isRHC = false 
)

Definition at line 1099 of file resolution2018.C.

References Beam(), bin, plot_validation_datamc_2018::color, om::cout, allTimeWatchdog::endl, genie::utils::style::Format(), MECModelEnuComparisons::i, nbins, Simulation(), ana::UniqueName(), and yaxis.

Referenced by Beam(), and resolution2018().

1100  {
1101  gPad->SetFillStyle(0);
1102  int nbins = allResol.size();
1103  TH2D* haux2D = new TH2D(UniqueName().c_str(),"; (True-Reco)/True;",1000,-0.5,0.5,
1104  nbins,0,nbins);
1105  if(resid) haux2D ->GetXaxis()->SetTitle("(True-Reco) (GeV)");
1106  haux2D->GetXaxis()->CenterTitle();
1107  auto yaxis = haux2D->GetYaxis();
1108  gPad->SetLeftMargin(.25);
1109  for(int i=0;i<nbins;++i){
1110  int bin = nbins - i;
1111  std::cout << "Test labels" << idl[i] << std::endl;
1112  yaxis->SetBinLabel(bin,idl[i]);
1113  }
1114  yaxis->SetLabelSize(0.07);
1115  haux2D->DrawClone();
1116 
1117  for(int i=0;i<nbins;++i){
1118  int bin = nbins - i;
1119  auto mygraph = new TGraphErrors();
1120  int color = idc[i];
1121  allResol[i]->ResetStats();
1122  mygraph->SetPoint(mygraph->GetN(),allResol[i]->GetMean(), yaxis->GetBinCenter(bin)) ;
1123  mygraph->SetPointError(mygraph->GetN()-1, allResol[i]->GetRMS(),0);
1124  mygraph->SetMarkerStyle(kFullCircle);
1125  mygraph->SetLineColor(color);
1126  mygraph->SetFillColor(color);
1127  mygraph->SetMarkerColor(color);
1128  mygraph->SetLineWidth(4);
1129  mygraph->DrawClone("pl");
1130 
1131  TLatex * ltx = new TLatex();
1132  ltx->SetTextAlign(22);
1133  ltx->SetTextSize(0.04);
1134  if(resid){
1135  ltx->DrawLatex(0.3,yaxis->GetBinCenter(bin),
1136  TString::Format("#color[%d]{#splitline{Mean %.1f MeV}{RMS %.1f MeV}}",
1137  color,
1138  allResol[i]->GetMean()*1000, allResol[i]->GetRMS()*1000));
1139 
1140  }else{
1141  ltx->DrawLatex(0.3,yaxis->GetBinCenter(bin),
1142  TString::Format("#color[%d]{#splitline{Mean %.1f \%}{RMS %.1f \%}}",
1143  color,
1144  allResol[i]->GetMean()*100, allResol[i]->GetRMS()*100));
1145 
1146  }
1147  }
1148  gPad->SetGridx();
1149  Simulation(); Beam(isRHC);
1150  delete haux2D;
1151  }
void Simulation()
Definition: tools.h:16
std::string yaxis
const int nbins
Definition: cellShifts.C:15
float bin[41]
Definition: plottest35.C:14
OStream cout
Definition: OStream.cxx:6
void Format(TGraph *gr, int lcol, int lsty, int lwid, int mcol, int msty, double msiz)
Definition: Style.cxx:154
void Beam(bool isRHC)
std::string UniqueName()
Return a different string each time, for creating histograms.
Definition: Utilities.cxx:29
TCanvas * ResolutionPanels ( TString  cname,
const std::vector< TString >  labels,
std::string  epoch,
bool  useTrue,
double  max1,
bool  fracOnly,
bool  isRHC 
)

Definition at line 878 of file resolution2018.C.

References Beam(), plot_validation_datamc::c, FormatResHist(), analysePickle::hist, MECModelEnuComparisons::i, Simulation(), ana::UniqueName(), xaxis, and yaxis.

Referenced by Beam().

878  {
879 
880  gStyle->SetTitleSize(17, "xyz");
881  gStyle->SetTitleFont(43, "xyz");
882  gStyle->SetLabelSize(15, "xyz");
883  gStyle->SetLabelFont(43, "xyz");
884  gStyle->SetTitleOffset(0.8, "x");
885  gStyle->SetTitleOffset(1.4, "yz");
886 
887  TCanvas * c ;
888  TH1D* hist = new TH1D (UniqueName().c_str(),"",230,-0.25,5.5);
889  TString xaxis,yaxis;
890  xaxis = "Reconstructed neutrino energy (GeV)";
891  if(useTrue) xaxis = "True neutrino energy (GeV)";
892  int npads= labels.size();
893 
894  //Divide the canvas and construct the panels
895  if(fracOnly){
896  c = new TCanvas (cname+"f",cname,700,550);
897  c->Divide(1,2);
898  for(int i=1;i<=2;++i){
899  c->cd(i);
900  if(i==1) FormatResHist(hist, xaxis, "Frac. E resolution", false, 0, 0.18);
901  if(i==2) FormatResHist(hist, xaxis, "Frac. shift in mean", false, -0.18, 0.18);
902 
903  gPad->SetTopMargin(0.15);
904  gPad->SetBottomMargin(0.18);
905  gPad->SetRightMargin(0.02*5./7.);
906  gPad->SetLeftMargin(0.1 *5./7.);
907 
908  hist->GetXaxis()->SetTitleOffset(2);
909  hist->GetYaxis()->SetTitleOffset(1.2);
910 
911  hist->DrawCopy();
912  gPad->SetGridx();
913  gPad->SetFillStyle(0);
914  gPad->SetGridy();
915  }
916  }
917  else{
918  c = new TCanvas(cname,cname,700,700);
919  double del = 0.66/npads;
920  c->Divide(1,1+npads);
921  c->GetPad(1)->SetPad(0,0.66,1,1);
922  for(int i=0;i<npads; ++i){
923  c->GetPad(2+i)->SetPad(0,0.66-(i+1)*del,1,0.66-i*del);
924  }
925  for(int i=1; i<= 1+npads ; ++i){
926  c->cd(i);
927  gPad->SetBottomMargin((i>1?0.18:0.16));
928  gPad->SetRightMargin(0.02);
929  gPad->SetLeftMargin(0.08);
930  if(i>1) gPad->SetTopMargin(0.05);
931  else gPad->SetTopMargin(0.15);
932 
933  if(i==1) FormatResHist(hist, xaxis, "Events/bin", false, 0, max1);
934  else FormatResHist(hist, "Energy difference (True-Reco) (GeV)",
935  "A.U.", true, 0, 0.77);
936  hist->DrawCopy();
937  gPad->SetGridx();
938  gPad->SetFillStyle(0);
939  }
940  }
941 
942 
943  TLatex * ltx = new TLatex();
944  ltx->SetTextSize(19);
945  ltx->SetTextFont(43);
946  ltx->SetTextAlign(11);
947 
948  c->cd();
949  gPad->SetFillStyle(0);
950 
951  if(!fracOnly) {
952  ltx->DrawLatexNDC(0.13,0.9,"#splitline{FD MC}{6 #times 10^{20} POT}");
953  for(int i = 0; i < npads; ++i){
954  c->cd(i+2);
955  ltx->DrawLatexNDC(0.13,0.6, "#bf{"+labels[i]+"}");
956  }
957  }
958  c->cd(1); Simulation(); Beam(isRHC);
959  if(fracOnly) c->cd(2); Simulation(); Beam(isRHC);
960 
961  return c;
962 
963 }
void Simulation()
Definition: tools.h:16
std::string yaxis
std::string xaxis
void FormatResHist(TH1D *hist, TString xaxis, TString yaxis, bool choplabel, double xmin, double xmax)
void Beam(bool isRHC)
std::string UniqueName()
Return a different string each time, for creating histograms.
Definition: Utilities.cxx:29