27 #include "Utilities/func/MathUtil.h" 49 double dCP = 1.2*TMath::Pi();
59 TH2 *
ret =
new TH2F(
UniqueName().c_str(),
"",50,0.5,1,50,0.35,0.85);
60 for (
int i = 1;
i <= 2500;
i++){
61 int xbin = (
i-1)%50 + 1;
62 int ybin = (
i-1)/50 + 1;
63 ret->SetBinContent(xbin,ybin,h->GetBinContent(
i));
97 std::vector<double> sigs(50,0);
98 std::vector<double> tots(50,0);
101 for (
int i = 1;
i <= 50;
i++){
105 for (
int j = 1;
j <= 50;
j++){
107 double tot = htot->Integral(
i,
i,
j,50);
108 if (tot == 0)
continue;
109 double cfom = sig/
sqrt(tot);
112 bestY = 0.35 + 0.01*(
j-1);
118 for (
int j = 0;
j <
i;
j++){
120 sigs[
j] += hsig->Integral(i,i,(bestY-0.35)/0.01+1,50);
121 tots[
j] += htot->Integral(i,i,(bestY-0.35)/0.01+1,50);
125 TH1 *
h =
new TH1F(
UniqueName().c_str(),
"",50,0.5,1);
126 for (
int i = 50;
i >= 1;
i--){
127 h->SetBinContent(
i,sigs[
i-1]/
sqrt(tots[
i-1]));
130 h->GetXaxis()->SetTitle(
"CVN-SS");
131 h->GetYaxis()->SetTitle(
"Total FOM");
132 h->GetXaxis()->CenterTitle();
133 h->GetYaxis()->CenterTitle();
135 TCanvas *
can =
new TCanvas();
175 for (
int i = 1;
i <= 50;
i++){
176 double sig1 = hsig->Integral(
i,50);
177 double tot1 = htot->Integral(
i,50);
178 double fom1 = sig1/
sqrt(tot1);
182 best1 = 0.5 + 0.01*(
i-1);
190 for (
int i = 1;
i <= 49;
i++){
191 for (
int j =
i+1;
j <= 50;
j++){
192 double sig1 = hsig->Integral(
i,
j-1);
193 double tot1 = htot->Integral(
i,
j-1);
194 double fom1 = sig1/
sqrt(tot1);
195 double sig2 = hsig->Integral(
j,50);
196 double tot2 = htot->Integral(
j,50);
197 double fom2 = sig2/
sqrt(tot2);
198 if (
sqrt(fom1*fom1+fom2*fom2) > fom){
199 fom =
sqrt(fom1*fom1+fom2*fom2);
200 best1 = 0.5 + 0.01*(
i-1);
201 best2 = 0.5 + 0.01*(
j-1);
210 for (
int i = 1;
i <= 48;
i++){
211 for (
int j =
i+1;
j <= 49;
j++){
212 for (
int k =
j+1; k <= 50; k++){
213 double sig1 = hsig->Integral(
i,
j-1);
214 double tot1 = htot->Integral(
i,
j-1);
215 double fom1 = sig1/
sqrt(tot1);
217 double sig2 = hsig->Integral(
j,k-1);
218 double tot2 = htot->Integral(
j,k-1);
219 double fom2 = sig2/
sqrt(tot2);
221 double sig3 = hsig->Integral(k,50);
222 double tot3 = htot->Integral(k,50);
223 double fom3 = sig2/
sqrt(tot2);
225 double tfom =
sqrt(fom1*fom1+fom2*fom2+fom3*fom3);
228 best1 = 0.5 + 0.01*(
i-1);
229 best2 = 0.5 + 0.01*(
j-1);
230 best3 = 0.5 + 0.01*(k-1);
235 std::cout << fom <<
": " << best1 <<
"-" << best2 <<
"-" << best3 <<
"-1" <<
std::endl;
247 double cvn = sr->
sel.cvnProd3Train.nueid;
249 if (cvn > 0.99 || (cvn>0.95 && bdt > 0.57))
return true;
258 double cvn = sr->
sel.cvnProd3Train.nueid;
261 if (cvn > 0.75 && cvn < 0.87)
263 if (cvn > 0.87 && cvn < 0.95)
275 int nuEBin = nuE/0.5;
276 assert(nuEBin <= 8 && "An event with nuE > 4.5 should never happen
"); 278 int anaBin = 9*selBin + nuEBin; 287 const Var kCVNSSe([](const caf::SRProxy* sr){ 288 return sr->sel.cvnProd3Train.nueid; 291 const Cut kNue2018FDCore = kCVNSSe > 0.75 && kNue2017CorePresel; 292 const Cut kNue2018FD = kNue2018FDCore || kNue2018FDPeripheral; 296 // Macro for optimizing PID cuts / binning in the FD. Just works on RHC now 297 // TODO: incorperate FHC. Update presel cuts if we decide to change them 298 void make_nue_pidtuning(bool makePlots=false) 303 SpectrumLoader lSwap("dpershey_prod4-cosrej_rhc_mc_fluxswap
"); 304 SpectrumLoader lNonS("dpershey_prod4-cosrej_rhc_mc_nonswap
"); 305 SpectrumLoader lTauS("dpershey_prod4-cosrej_rhc_mc_tauswap
"); 307 SpectrumLoader lCosmic("dpershey_prod4-cosrej_cosmicdata_rhctune
"); 309 lSwap.SetSpillCut(kStandardSpillCuts); 310 lNonS.SetSpillCut(kStandardSpillCuts); 311 lTauS.SetSpillCut(kStandardSpillCuts); 312 lCosmic.SetSpillCut(kStandardSpillCuts); 314 Cut selCore = kNue2017CorePresel; 315 Cut selPeri = kNue2017PeripheralPresel; 316 Var wei = kXSecCVWgt2017*kPPFXFluxCVWgt; 320 // Peripheral cuts are 2D PID cuts, set up 2D vars now 321 Var kBDTLight = SIMPLEVAR(sel.nuecosrej.cospidlight); 322 Var kBDTCont = SIMPLEVAR(sel.nuecosrej.cospidcontain); 324 Binning bdtbins = Binning::Simple(50,0.35,0.85); 325 Binning cvnbins = Binning::Simple(50,0.5,1.); 327 Var kLightCVN = Var2D(kBDTLight,bdtbins, 329 Var kContCVN = Var2D(kBDTCont,bdtbins, 333 // Prediction as a function of CVN, for main PID tuning 334 // Be sure you're using the right core preselection, it may change 335 PredictionNoExtrap predCVN(lNonS,lSwap,lTauS,"", 336 Binning::Simple(50,0.5,1), 337 kCVNSSe, selCore, kNoShift, wei); 338 Spectrum cosCVN(lCosmic,HistAxis("",Binning::Simple(50,0.5,1),kCVNSSe), 339 kInCosmicTimingWindow&&selCore,kNoShift,wei); 341 // 2D predictions for determining peripheral PID cuts 342 PredictionNoExtrap predLightPeri(lNonS,lSwap,lTauS, 343 "",Binning::Simple(2500,0,2500),kLightCVN, 344 selPeri,kNoShift,wei); // BDTLight 345 Spectrum cosLightPeri(lCosmic, 346 HistAxis("",Binning::Simple(2500,0,2500),kLightCVN), 347 selPeri&&kInCosmicTimingWindow,kNoShift,wei); 348 PredictionNoExtrap predContPeri (lNonS,lSwap,lTauS, 349 "",Binning::Simple(2500,0,2500),kContCVN, 350 selPeri,kNoShift,wei); // BDTContain 351 Spectrum cosContPeri (lCosmic, 352 HistAxis("",Binning::Simple(2500,0,2500),kContCVN), 353 selPeri&&kInCosmicTimingWindow,kNoShift,wei); 355 // Prediction with 2017 bin boundaries, just for fun 356 PredictionNoExtrap pred2017(lNonS,lSwap,lTauS,"", 357 Binning::Simple(36,0,36), 359 kNue2018FD,kNoShift,wei); 360 Spectrum cos2017(lCosmic, 361 HistAxis("",Binning::Simple(36,0,36),kNue2018AnaBin), 362 kInCosmicTimingWindow&&kNue2018FD,kNoShift,wei); 377 TFile* outFile = new TFile("out_fdcosrej.root
","RECREATE
"); 380 pred2017.SaveTo(outFile, "pred2017
"); 381 cos2017.SaveTo(outFile, "cos2017
"); 383 predCVN.SaveTo(outFile, "predCVN
"); 384 cosCVN.SaveTo(outFile, "cosCVN
"); 386 predLightPeri.SaveTo(outFile, "predLightPeri
"); 387 cosLightPeri .SaveTo(outFile, "cosLightPeri
"); 389 predContPeri.SaveTo(outFile, "predContPeri
"); 390 cosContPeri .SaveTo(outFile, "cosContPeri
"); 394 } // End if !makePlots 397 PeripheralCuts("LightPeri
");
Cuts and Vars for the 2020 FD DiF Study.
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.
const FitDmSq32 kFitDmSq32
void SetValue(osc::IOscCalcAdjustable *osc, double val) const override
Proxy for caf::StandardRecord.
void SetValue(osc::IOscCalcAdjustable *osc, double val) const override
Forward to wrapped Var's SetValue()
double Integral(double exposure, double *err=0, EExposureType expotype=kPOT) const
Return total number of events scaled to pot.
void SetValue(osc::IOscCalcAdjustable *osc, double val) const override
Defines an enumeration for prong classification.
osc::IOscCalcAdjustable * DefaultOscCalc()
Create a new calculator with default assumptions for all parameters.
Representation of a spectrum in any variable, with associated POT.
const Var kNue2018AnaBin([](const caf::SRProxy *sr){int selBin=kNue2018SelectionBin(sr);float nuE=kNueEnergy2018(sr);int nuEBin=nuE/0.5;assert(nuEBin<=8 &&"An event with nuE > 4.5 should never happen");int anaBin=9 *selBin+nuEBin;return anaBin;})
Use this Analysis Binning for Ana2018, official Binning.
caf::Proxy< caf::SRNueCosRej > nuecosrej
Charged-current interactions.
const Var kNueEnergy2017([](const caf::SRProxy *sr){if(!sr->vtx.elastic.IsValid) return-1.0;if(sr->vtx.elastic.fuzzyk.npng< 1) return-1.0;return NueRecoE_2017FDFit(kCVNemE(sr), kCVNhadE(sr));})
void SetValue(osc::IOscCalcAdjustable *osc, double val) const override
void CVNCuts(std::string name)
const Var kNue2018SelectionBin([](const caf::SRProxy *sr){bool isPeri=kNue2018FDPeripheral(sr);if(isPeri) return float(2.5);std::cout<< "ERROR::kNue2018SelectionBin. Looking for cvnProd3Train. Branch no longer exists."<< std::endl;abort();return float(-5.0);})
const Cut kNue2017PeripheralPresel
const Cut kNue2018FDPeripheral(kNue2018FDPeripheralFunc)
const Cut kNue2018CVNVsCosPID(kNue2018CVNVsCosPIDFunc)
const ConstrainedFitVarWithPrior fitSsqTh23_UniformPriorSsqTh23 & kFitSinSqTheta23
assert(nhit_max >=nhit_nbins)
caf::Proxy< float > cospidcontain
const FitSinSq2Theta13 kFitSinSq2Theta13
caf::Proxy< caf::SRIDBranch > sel
Template for Cut and SpillCut.
const FitVarWithPrior fitDeltaInPiUnits_UniformPriordCP & kFitDeltaInPiUnits
std::string UniqueName()
Return a different string each time, for creating histograms.
void PeripheralCuts(std::string name)