6 #include "TDirectory.h" 28 TH1* histMC( mc.
ToTH1(pot) );
29 TH1* histNotMC( notMC.
ToTH1(pot) );
30 histMC->SetLineColor(
kRed );
31 histNotMC->SetLineColor( notMCColor );
33 histMC->SetMaximum( 1.5 * histMC->GetMaximum() );
34 if (restrictRange) histMC->GetXaxis()->SetRangeUser( 0.5, 3.5);
35 histMC->SetTitle( (latex+
" "+title).c_str() );
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);
44 histMC->Write( (TString)
"MC"+saveAs.c_str());
45 histNotMC->Write( (TString)notMCLabel.c_str() + saveAs.c_str());
47 c.Write( saveAs.c_str() );
56 TDirectory*
tmp = gDirectory;
59 TH1*
hist( fRecoFD.Unoscillated().ToTH1(potFD) );
62 hist->SetMaximum( 1.5 *
hist->GetMaximum() );
63 hist->SetTitle(
"FD Reco Spectrum" );
75 TDirectory*
tmp = gDirectory;
79 double potND( decompresult.
POT() );
81 fRecoToTrueND.Unoscillated(), decompresult, potND,
82 "Data", kBlack, fLatex,
"ND Reco Spectrum",
"NDReco" );
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");
96 fRecoToTrueND.TrueEnergy(), recoToTrueND.
TrueEnergy(), potND,
97 "Reweighted",
kBlue, fLatex,
"ND True Energy Spectrum",
"NDTrue",
102 "Extrapolated",
kBlue, fLatex,
"FD True Energy Spectrum",
"FDTrue",
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");
112 "Extrapolated",
kBlue, fLatex,
"FD Reco Spectrum",
"FDReco" );
114 TH1* histNDDatReco = decompresult.
ToTH1(potND);
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);
124 fTrueToRecoFD.SaveTo(dir,
"trueToRecoMC");
127 histNDDatReco->Write(
"NDDatReco");
128 histNDRwtTrue->Write(
"NDRwtTrue");
129 histFDExtTrue->Write(
"FDExtTrue");
130 histFDExtReco->Write(
"FDExtReco");
132 histNDMCReco->Write(
"NDMCReco");
133 histNDMCTrue->Write(
"NDMCTrue");
134 histFDMCTrue->Write(
"FDMCTrue");
135 histFDMCReco->Write(
"FDMCReco");
137 delete histNDDatReco;
138 delete histNDRwtTrue;
139 delete histFDExtTrue;
140 delete histFDExtReco;
154 TDirectory*
tmp = gDirectory;
157 Spectrum fdMCReco( fTrueToRecoFD.Unoscillated() );
161 double potND( decompresult.
POT() );
164 fRecoND, decompresult, potND,
165 "Data", kBlack, fLatex,
"ND Reco Spectrum",
"ND" );
168 fdMCReco, fdExtReco, potFD,
169 "Extrapolated",
kBlue, fLatex,
"FD Reco Spectrum",
"FD" );
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 );
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);
189 cDvMC.Write(
"DataOverMC");
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 );
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);
211 cFoN.Write(
"FarOverNear");
215 TH1* histNDDat = decompresult.
ToTH1(potND);
216 TH1* histFDExt = fdExtReco.ToTH1(potFD);
217 TH1* histNDMC = fRecoND.ToTH1(potND);
218 TH1* histFDMC = fdMCReco.ToTH1(potFD);
220 histNDDat->Write(
"NDDat");
221 histFDExt->Write(
"FDExt");
222 histNDMC->Write(
"NDMC");
223 histFDMC->Write(
"FDMC");
240 TH1*
fd = fTrueToRecoFD.Unoscillated().ToTH1(potMCFD);
241 TH1*
nd = fRecoND.ToTH1(potMCND);
243 TH1* ob = OptimalBinningHelper( fd, nd );
259 TH1*
fd = fTrueToRecoFD.Unoscillated().ToTH1(potMCFD);
260 TH1*
nd = fRecoND.ToTH1(potMCND);
262 const double lower( fd->GetXaxis()->GetXmin() );
263 const double upper( fd->GetXaxis()->GetXmax() );
265 TF1* fdGauss =
new TF1(
UniqueName().c_str() ,
"gaus(0)", lower, upper );
266 TF1* ndGauss =
new TF1(
UniqueName().c_str() ,
"gaus(0)", lower, upper );
268 fdGauss->SetParameters( 100., 2., .5 );
269 ndGauss->SetParameters( 100., 2., .5 );
271 fd->Fit( fdGauss,
"WW" );
272 nd->Fit( ndGauss,
"WW" );
274 TH1* fdFine =
new TH1D(
UniqueName().c_str(),
"fdFine",100,lower,upper);
275 TH1* ndFine =
new TH1D(
UniqueName().c_str(),
"ndFine",100,lower,upper);
277 fdFine->Eval(fdGauss);
278 ndFine->Eval(ndGauss);
280 const double binScale = fdFine->GetBinWidth(1) / fd->GetBinWidth(1);
282 fdFine->Scale(binScale);
283 ndFine->Scale(binScale);
285 TH1* ob(OptimalBinningHelper(fdFine,ndFine));
298 TH1* ob( (TH1*)fd->Clone() );
300 int nbins = fd->GetNbinsX();
301 assert( (nbins == nd->GetNbinsX()) &&
"Bin Mismatch" );
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.) );
317 else ob->SetBinContent(
bin, 0. );
319 ob->SetBinContent( 0, 0. );
320 ob->SetBinContent( nbins+1, 0. );
Cuts and Vars for the 2020 FD DiF Study.
bool isfinite(const stan::math::var &v)
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.
void SavePlots(TDirectory *dir, double potFD) const override
TH1 * OptimalBinning(double potMCFD, double potMCND) const
Uses MC spectra and target MC stats to estimate optimal binning.
void ReweightToRecoSpectrum(const Spectrum &target)
Recale bins so that Unweighted will return target.
def ratio(spec1, spec2, nbins, pot, binrange=(0, 1))
TH1 * OptimalBinningFit(double potMCFD, double potMCND) const
Fits gaussians to FD and ND MC to estimate optimal binning.
Representation of a spectrum in any variable, with associated POT.
Spectrum Unoscillated() const
Spectrum TrueEnergy() const
static TH1 * OptimalBinningHelper(TH1 *, TH1 *)
void SavePlots(TDirectory *dir, double potFD) const override
void SavePlots(TDirectory *dir, double potFD) const override
assert(nhit_max >=nhit_nbins)
Spectrum with true energy information, allowing it to be oscillated
std::string UniqueName()
Return a different string each time, for creating histograms.
void SaveTo(TDirectory *dir, const std::string &name) const