ModularExtrapComponentPlot.cxx
Go to the documentation of this file.
1 #include "CAFAna/Extrap/ModularExtrapComponent.h"
2 #include "CAFAna/Core/Ratio.h"
3 #include "CAFAna/Core/Utilities.h"
4 
5 #include "TCanvas.h"
6 #include "TDirectory.h"
7 #include "TLegend.h"
8 #include "TF1.h"
9 #include "TH1.h"
10 #include "TH2.h"
11 
12 namespace ana
13 {
14 
15  //---------------------------------------------------------------------------
16 
17  void ModularExtrapComponent::ComparisonPlot(
18  Spectrum mc,
19  Spectrum notMC,
20  double pot,
21  std::string notMCLabel,
22  int notMCColor,
23  std::string latex,
24  std::string title,
25  std::string saveAs,
26  bool restrictRange
27  ){
28  TH1* histMC( mc.ToTH1(pot) );
29  TH1* histNotMC( notMC.ToTH1(pot) );
30  histMC->SetLineColor( kRed );
31  histNotMC->SetLineColor( notMCColor );
32  TCanvas c;
33  histMC->SetMaximum( 1.5 * histMC->GetMaximum() );
34  if (restrictRange) histMC->GetXaxis()->SetRangeUser( 0.5, 3.5);
35  histMC->SetTitle( (latex+" "+title).c_str() );
36  histMC->Draw("hist");
37  histNotMC->Draw("same");
38  TLegend legend( .15, .7, .4, .85 );
39  legend.SetTextSize(0.05);
40  legend.AddEntry( histNotMC, notMCLabel.c_str() );
41  legend.AddEntry( histMC, "MC" );
42  legend.SetFillStyle(0);
43  legend.Draw();
44  histMC->Write( (TString)"MC"+saveAs.c_str());
45  histNotMC->Write( (TString)notMCLabel.c_str() + saveAs.c_str());
46 
47  c.Write( saveAs.c_str() );
48  delete histMC;
49  delete histNotMC;
50  }
51 
52  //---------------------------------------------------------------------------
53 
54  void NoReweight::SavePlots(TDirectory* dir, double potFD) const
55  {
56  TDirectory* tmp = gDirectory;
57  dir->cd();
58 
59  TH1* hist( fRecoFD.Unoscillated().ToTH1(potFD) );
60  hist->SetLineColor( kRed );
61  TCanvas c;
62  hist->SetMaximum( 1.5 * hist->GetMaximum() );
63  hist->SetTitle( "FD Reco Spectrum" );
64  hist->Draw("hist");
65  c.Write("MC");
66  delete hist;
67 
68  tmp->cd();
69  }
70 
71  //---------------------------------------------------------------------------
72 
73  void TruthReweight::SavePlots(TDirectory* dir, double potFD) const
74  {
75  TDirectory* tmp = gDirectory;
76  dir->cd();
77 
78  Spectrum decompresult( GetDecompResult(fDecomp,fDecompRes) );
79  double potND( decompresult.POT() );
80  ComparisonPlot(
81  fRecoToTrueND.Unoscillated(), decompresult, potND,
82  "Data", kBlack, fLatex, "ND Reco Spectrum", "NDReco" );
83 
84  TH2* histNDrtt = fRecoToTrueND.ToTH2(potND);
85  histNDrtt->GetYaxis()->SetRangeUser(0, 5);
86  histNDrtt->GetXaxis()->SetRangeUser(0, 5);
87  histNDrtt->SetTitle( (fLatex+" ND Reco to True Spectrum (MC)").c_str() );
88  histNDrtt->Draw("colz");
89  histNDrtt->Write("NDRecoToTrue");
90 
91  OscillatableSpectrum recoToTrueND( fRecoToTrueND );
92  OscillatableSpectrum trueToRecoFD( fTrueToRecoFD );
93 
94  recoToTrueND.ReweightToRecoSpectrum( decompresult );
95  ComparisonPlot(
96  fRecoToTrueND.TrueEnergy(), recoToTrueND.TrueEnergy(), potND,
97  "Reweighted", kBlue, fLatex, "ND True Energy Spectrum", "NDTrue",
98  true );
99 
100  ComparisonPlot(
101  fTrueToRecoFD.TrueEnergy(), Return().TrueEnergy(), potFD,
102  "Extrapolated", kBlue, fLatex, "FD True Energy Spectrum", "FDTrue",
103  true );
104 
105  TH2* histFDttr = fTrueToRecoFD.ToTH2(potFD);
106  histFDttr->SetTitle( (fLatex+" FD True To Reco Spectrum (MC)").c_str() );
107  histFDttr->Draw("colz");
108  histFDttr->Write("FDTrueToReco");
109 
110  ComparisonPlot(
111  fTrueToRecoFD.Unoscillated(), Return().Unoscillated(), potFD,
112  "Extrapolated", kBlue, fLatex, "FD Reco Spectrum", "FDReco" );
113 
114  TH1* histNDDatReco = decompresult.ToTH1(potND);
115  TH1* histNDRwtTrue = recoToTrueND.TrueEnergy().ToTH1(potND);
116  TH1* histFDExtTrue = Return().TrueEnergy().ToTH1(potFD);
117  TH1* histFDExtReco = Return().Unoscillated().ToTH1(potFD);
118 
119  TH1* histNDMCReco = fRecoToTrueND.Unoscillated().ToTH1(potND);
120  TH1* histNDMCTrue = fRecoToTrueND.TrueEnergy().ToTH1(potND);
121  TH1* histFDMCTrue = fTrueToRecoFD.TrueEnergy().ToTH1(potFD);
122  TH1* histFDMCReco = fTrueToRecoFD.Unoscillated().ToTH1(potFD);
123 
124  fTrueToRecoFD.SaveTo(dir, "trueToRecoMC");
125  Return().SaveTo(dir, "trueToRecoExt");
126 
127  histNDDatReco->Write("NDDatReco");
128  histNDRwtTrue->Write("NDRwtTrue");
129  histFDExtTrue->Write("FDExtTrue");
130  histFDExtReco->Write("FDExtReco");
131 
132  histNDMCReco->Write("NDMCReco");
133  histNDMCTrue->Write("NDMCTrue");
134  histFDMCTrue->Write("FDMCTrue");
135  histFDMCReco->Write("FDMCReco");
136 
137  delete histNDDatReco;
138  delete histNDRwtTrue;
139  delete histFDExtTrue;
140  delete histFDExtReco;
141 
142  delete histNDMCReco;
143  delete histNDMCTrue;
144  delete histFDMCTrue;
145  delete histFDMCReco;
146 
147  tmp->cd();
148  }
149 
150  //---------------------------------------------------------------------------
151 
152  void RecoReweight::SavePlots(TDirectory* dir, double potFD) const
153  {
154  TDirectory* tmp = gDirectory;
155  dir->cd();
156 
157  Spectrum fdMCReco( fTrueToRecoFD.Unoscillated() );
158  Spectrum fdExtReco( Return().Unoscillated() );
159 
160  Spectrum decompresult( GetDecompResult(*fDecomp,fDecompRes) );
161  double potND( decompresult.POT() );
162 
163  ComparisonPlot(
164  fRecoND, decompresult, potND,
165  "Data", kBlack, fLatex, "ND Reco Spectrum", "ND" );
166 
167  ComparisonPlot(
168  fdMCReco, fdExtReco, potFD,
169  "Extrapolated", kBlue, fLatex, "FD Reco Spectrum", "FD" );
170 
171  TH1* histDvMCND( (decompresult/fRecoND).ToTH1() );
172  TH1* histEvMCFD( (fdExtReco / fdMCReco).ToTH1() );
173  histDvMCND->SetLineColor( kRed );
174  histEvMCFD->SetLineColor( kBlue );
175  histEvMCFD->SetLineWidth( 4 );
176  histEvMCFD->SetLineStyle( 2 );
177  TCanvas cDvMC;
178  histDvMCND->SetMaximum( 1.5 * histDvMCND->GetMaximum() );
179  histDvMCND->SetTitle( (fLatex+" Data/MC Ratio").c_str() );
180  histDvMCND->GetYaxis()->SetTitle("Ratio");
181  histDvMCND->Draw("hist");
182  histEvMCFD->Draw("same""hist");
183  TLegend legendDvMC( .15, .7, .4, .85 );
184  legendDvMC.SetTextSize(0.05);
185  legendDvMC.AddEntry( histEvMCFD, "FD (Extrapolated)" );
186  legendDvMC.AddEntry( histDvMCND, "ND" );
187  legendDvMC.SetFillStyle(0);
188  legendDvMC.Draw();
189  cDvMC.Write("DataOverMC");
190  delete histDvMCND;
191  delete histEvMCFD;
192 
193  TH1* histFoNMC( (fdMCReco/fRecoND).ToTH1() );
194  TH1* histFoNDE( (fdExtReco/decompresult).ToTH1() );
195  histFoNMC->SetLineColor( kRed );
196  histFoNDE->SetLineColor( kBlue );
197  histFoNDE->SetLineWidth( 4 );
198  histFoNDE->SetLineStyle( 2 );
199  TCanvas cFoN;
200  histFoNMC->SetMaximum( 1.5 * histFoNMC->GetMaximum() );
201  histFoNMC->SetTitle( (fLatex+" Far/Near Ratio").c_str() );
202  histFoNMC->GetYaxis()->SetTitle("Ratio");
203  histFoNMC->Draw("hist");
204  histFoNDE->Draw("same""hist");
205  TLegend legendFoN( .15, .7, .4, .85 );
206  legendFoN.SetTextSize(0.05);
207  legendFoN.AddEntry( histFoNDE, "Extrapolation/Data" );
208  legendFoN.AddEntry( histFoNMC, "MC" );
209  legendFoN.SetFillStyle(0);
210  legendFoN.Draw();
211  cFoN.Write("FarOverNear");
212  delete histFoNMC;
213  delete histFoNDE;
214 
215  TH1* histNDDat = decompresult.ToTH1(potND);
216  TH1* histFDExt = fdExtReco.ToTH1(potFD);
217  TH1* histNDMC = fRecoND.ToTH1(potND);
218  TH1* histFDMC = fdMCReco.ToTH1(potFD);
219 
220  histNDDat->Write("NDDat");
221  histFDExt->Write("FDExt");
222  histNDMC->Write("NDMC");
223  histFDMC->Write("FDMC");
224 
225  delete histNDDat;
226  delete histFDExt;
227  delete histNDMC;
228  delete histFDMC;
229 
230  tmp->cd();
231  }
232 
233  //---------------------------------------------------------------------------
234 
235  TH1* RecoReweight::OptimalBinning(
236  double potMCFD,
237  double potMCND
238  ) const {
239 
240  TH1* fd = fTrueToRecoFD.Unoscillated().ToTH1(potMCFD);
241  TH1* nd = fRecoND.ToTH1(potMCND);
242 
243  TH1* ob = OptimalBinningHelper( fd, nd );
244 
245  delete fd;
246  delete nd;
247 
248  return ob;
249 
250  }
251 
252  //---------------------------------------------------------------------------
253 
254  TH1* RecoReweight::OptimalBinningFit(
255  double potMCFD,
256  double potMCND
257  ) const {
258 
259  TH1* fd = fTrueToRecoFD.Unoscillated().ToTH1(potMCFD);
260  TH1* nd = fRecoND.ToTH1(potMCND);
261 
262  const double lower( fd->GetXaxis()->GetXmin() );
263  const double upper( fd->GetXaxis()->GetXmax() );
264 
265  TF1* fdGauss = new TF1( UniqueName().c_str() ,"gaus(0)", lower, upper );
266  TF1* ndGauss = new TF1( UniqueName().c_str() ,"gaus(0)", lower, upper );
267 
268  fdGauss->SetParameters( 100., 2., .5 );
269  ndGauss->SetParameters( 100., 2., .5 );
270 
271  fd->Fit( fdGauss, "WW" );
272  nd->Fit( ndGauss, "WW" );
273 
274  TH1* fdFine = new TH1D(UniqueName().c_str(),"fdFine",100,lower,upper);
275  TH1* ndFine = new TH1D(UniqueName().c_str(),"ndFine",100,lower,upper);
276 
277  fdFine->Eval(fdGauss);
278  ndFine->Eval(ndGauss);
279 
280  const double binScale = fdFine->GetBinWidth(1) / fd->GetBinWidth(1);
281 
282  fdFine->Scale(binScale);
283  ndFine->Scale(binScale);
284 
285  TH1* ob(OptimalBinningHelper(fdFine,ndFine));
286 
287  delete fd;
288  delete nd;
289  return ob;
290 
291  }
292 
293  //---------------------------------------------------------------------------
294 
295  TH1* RecoReweight::OptimalBinningHelper( TH1* fd, TH1* nd)
296  {
297 
298  TH1* ob( (TH1*)fd->Clone() );
299 
300  int nbins = fd->GetNbinsX();
301  assert( (nbins == nd->GetNbinsX()) && "Bin Mismatch" );
302 
303  for (int bin(1); bin<=nbins; ++bin)
304  {
305  double invdensFD( fd->GetBinWidth(bin) / fd->GetBinContent(bin) );
306  double invdensND( nd->GetBinWidth(bin) / nd->GetBinContent(bin) );
307  double ratio( fd->GetBinContent(bin) / nd->GetBinContent(bin) );
308  double ratioPrev( fd->GetBinContent(bin-1) / nd->GetBinContent(bin-1) );
309  double ratioNext( fd->GetBinContent(bin+1) / nd->GetBinContent(bin+1) );
310  double deltaRatio( ratioNext - ratioPrev );
311  double deltaVar( fd->GetBinCenter(bin+1) - fd->GetBinCenter(bin-1) );
312  double invslope( deltaVar / deltaRatio );
313  double widthCubed( 6. * ( ratio*ratio * invslope*invslope )
314  * ( invdensFD + invdensND ) );
315  double optbin( pow(widthCubed, 1./3.) );
316  if ( std::isfinite(optbin) ) ob->SetBinContent( bin, optbin );
317  else ob->SetBinContent( bin, 0. );
318  }
319  ob->SetBinContent( 0, 0. );
320  ob->SetBinContent( nbins+1, 0. );
321 
322  return ob;
323 
324  }
325 
326 }
latex
void latex(double x, double y, std::string txt, double ang=0, int align=22, double size=0.042)
Definition: make_extrap_figure.C:67
nd
bool nd
Definition: attenuation_calibration_comparisons.C:31
ana
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
boost::math::isfinite
bool isfinite(const stan::math::var &v)
Definition: boost_isfinite.hpp:21
ana::Spectrum::ToTH1
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
ratio
TH1 * ratio(TH1 *h1, TH1 *h2)
Definition: HistogramAttr.cxx:59
cet::pow
constexpr T pow(T x)
Definition: pow.h:75
ana::TruthReweight::SavePlots
void SavePlots(TDirectory *dir, double potFD) const override
Definition: ModularExtrapComponentPlot.cxx:73
ana::ModularExtrapComponent::Return
const OscillatableSpectrum & Return() const
Interface so that result of Eval() is called only once and cached.
Definition: ModularExtrapComponent.cxx:68
ana::RecoReweight::OptimalBinning
TH1 * OptimalBinning(double potMCFD, double potMCND) const
Uses MC spectra and target MC stats to estimate optimal binning.
Definition: ModularExtrapComponentPlot.cxx:235
tmp
Float_t tmp
Definition: plot.C:36
ana::ReweightableSpectrum::ReweightToRecoSpectrum
void ReweightToRecoSpectrum(const Spectrum &target)
Recale bins so that Unweighted will return target.
Definition: ReweightableSpectrum.cxx:151
make_mec_shifts_plots.legend
legend
Definition: make_mec_shifts_plots.py:148
fd
bool fd
Definition: attenuation_calibration_comparisons.C:32
ana::RecoReweight::OptimalBinningFit
TH1 * OptimalBinningFit(double potMCFD, double potMCND) const
Fits gaussians to FD and ND MC to estimate optimal binning.
Definition: ModularExtrapComponentPlot.cxx:254
ana::Spectrum
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
nbins
const int nbins
Definition: cellShifts.C:15
ana::RecoReweight::OptimalBinningHelper
static TH1 * OptimalBinningHelper(TH1 *, TH1 *)
Definition: ModularExtrapComponentPlot.cxx:295
plot_validation_datamc.c
c
Definition: plot_validation_datamc.py:198
pot
#define pot
Definition: template_basic.C:49
ana::RecoReweight::SavePlots
void SavePlots(TDirectory *dir, double potFD) const override
Definition: ModularExtrapComponentPlot.cxx:152
ana::NoReweight::SavePlots
void SavePlots(TDirectory *dir, double potFD) const override
Definition: ModularExtrapComponentPlot.cxx:54
bin
float bin[41]
Definition: plottest35.C:14
ana::Spectrum::POT
double POT() const
Definition: Spectrum.h:219
xml_schema::string
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
ana::ModularExtrapComponent::GetDecompResult
static Spectrum GetDecompResult(const IDecomp &, DecompResult)
Helper function to pick out single Spectrum from a decomposition.
Definition: ModularExtrapComponent.cxx:80
dir
TDirectory * dir
Definition: macro.C:5
mc
bool mc
Definition: attenuation_calibration_comparisons.C:33
ana::assert
assert(nhit_max >=nhit_nbins)
ana::OscillatableSpectrum
Spectrum with true energy information, allowing it to be oscillated
Definition: OscillatableSpectrum.h:36
ana::UniqueName
std::string UniqueName()
Return a different string each time, for creating histograms.
Definition: Utilities.cxx:29
ana::OscillatableSpectrum::SaveTo
void SaveTo(TDirectory *dir, const std::string &name) const
Definition: OscillatableSpectrum.cxx:208
ana::ModularExtrapComponent::ComparisonPlot
static void ComparisonPlot(Spectrum mc, Spectrum notMC, double pot, std::string notMCLabel, int notMCColor, std::string latex, std::string title, std::string saveAs, bool restrictRange=false)
Definition: ModularExtrapComponentPlot.cxx:17
Generated on Mon Oct 26 2020 05:13:48 for NOvA by   doxygen 1.8.11