15 #include "CAFAna/Core/Var.h" 38 while (OutName.find(
"_") != std::string::npos) OutName.replace(OutName.find(
"_"),1,
" ");
40 if (OutName.find(
"Qual") != std::string::npos) OutName.insert(OutName.find(
" Qual"),
", Cut level:");
42 if (OutName.find(
"pot") != std::string::npos) OutName.replace(OutName.find(
" pot"),4,
". With POT normalisation.");
44 if (OutName.find(
"area") != std::string::npos) OutName.replace(OutName.find(
" area"),5,
". With Area normalisation.");
46 if (OutName.find(
"fhc") != std::string::npos) OutName.replace(OutName.find(
" fhc"),4,
", in Neutrino mode");
47 if (OutName.find(
"rhc") != std::string::npos) OutName.replace(OutName.find(
" rhc"),4,
", in Antineutrino mode");
49 std::string Cap =
"Plot showing the number of ND NuMi events in #nu_{#mu} Ana18 for variable: " + OutName;
61 void ND_DataMC_Comp_Systs(
bool fhc =
true,
bool Set1Vars =
true,
bool IsTest=
false,
bool IsOnGrid=
true,
bool Weight18 =
true,
bool Prod4Cut =
true,
bool FullCuts =
false ) {
67 std::string sVarSet = (Set1Vars ==
true ?
"Set1" :
"Set2" );
68 std::string sWeight = (Weight18 ==
true ?
"wgt2018" :
"wgt2017" );
69 std::string sProd4 = (Prod4Cut ==
true ?
"CutProd4":
"CutProd3");
70 std::string sAllCut = (FullCuts ==
true ?
"_AllCuts":
"" );
72 std::cout <<
"\n\nLots of input parameters, so I'm going to write them out now....;" 73 <<
"\n\t IsFHC " << (
IsFHC ==
true ?
"true --> Run in FHC " :
"false --> Run in RHC" ) <<
" ==> sFHC_RHC is " << sFHC_RHC
74 <<
"\n\t Set1Vars " << (Set1Vars ==
true ?
"true --> Making first 9 plots" :
"false --> Making second 9 plots" ) <<
" ==> sVarSet is " << sVarSet
75 <<
"\n\t IsTest " << (IsTest ==
true ?
"true --> Short run use CAFs" :
"false --> Long run use Concat" )
76 <<
"\n\t Weight18 " << (Weight18 ==
true ?
"true --> Run using 2018 XSecs" :
"false --> Run using 2017 XSecs" ) <<
" ==> sWeight is " << sWeight
77 <<
"\n\t Prod4Cut " << (Prod4Cut ==
true ?
"true --> Using Prod4 Cuts" :
"false --> Using Prod3 Cuts" ) <<
" ==> sProd4 is " << sProd4
78 <<
"\n\t FullCuts " << (FullCuts ==
true ?
"true --> Running all cuts" :
"false --> Only running full cuts" ) <<
" ==> sAllCut is " << sAllCut
86 std::string fdspecfile =
"/pnfs/nova/persistent/analysis/numu/Ana2018/official/Quantiles/quantiles__"+sFHC_RHC+
"_full__numu2018.root";
91 TH2 *FDSpec2D = (TH2*)inFile->FindObjectAny(
"FDSpec2D" );
93 const int NHadEFracQuantiles = 4;
121 std::vector<Cut> TieredCuts;
130 CutNames.emplace_back(
"Qual"); TieredCuts.emplace_back(
kNumuQuality);
131 CutNames.emplace_back(
"NoCut"); TieredCuts.emplace_back(
kNoCut);
136 gStyle->SetMarkerStyle(kFullCircle);
137 TGaxis::SetMaxDigits(3);
146 fnameNDMC =
"karlwarb_ND_FHC_MC";
147 fnameNDData =
"karlwarb_ND_FHC_Data";
149 fnameNDMC =
"karlwarb_ND_RHC_MC";
150 fnameNDData =
"karlwarb_ND_RHC_Data";
155 fnameNDMC =
"prod_sumdecaf_R17-11-14-prod4reco.d_nd_genie_nonswap_fhc_nova_v08_full_v1_numu2018_stride22";
157 fnameNDData =
"prod_sumdecaf_R17-09-05-prod4recopreview.f_nd_numi_fhc_full_v1_addShortSimpleCVN_goodruns_numu2018";
159 fnameNDMC =
"prod_sumdecaf_R17-11-14-prod4reco.e_nd_genie_nonswap_rhc_nova_v08_full_v1_numu2018_stride14";
161 fnameNDData =
"prod_sumdecaf_R17-09-05-prod4recopreview.f_nd_numi_rhc_full_v1_addShortSimpleCVN_goodruns_numu2018";
165 std::string fnameODir =
"NDComp_"+sFHC_RHC+
"_"+sVarSet+
"_"+sWeight+
"_"+sProd4+sAllCut;
166 if (IsTest) fnameODir +=
"_Test";
168 std::cout <<
"\n\n ODir is " << fnameODir <<
", I loaded cuts from " << fdspecfile
169 <<
"\n\t fnameNDMC is " << fnameNDMC
170 <<
"\n\t fnameNDData is " << fnameNDData
181 if (IsOnGrid) fnameODir =
"";
182 else gSystem -> MakeDirectory( fnameODir.c_str() );
190 std::vector<Selection> selections;
192 for (
size_t quant=0; quant < HadEFracQuantCuts.size(); ++quant) {
195 selections.emplace_back( TieredCuts[
cuts], sFHC_RHC+
"_"+CutNames[cuts]+
"_QuantAll",
196 "Selected events pass "+CutNames[cuts]+
" all quantiles");
198 selections.emplace_back( TieredCuts[
cuts] && HadEFracQuantCuts[quant], sFHC_RHC+
"_"+CutNames[
cuts]+
"_Quant"+
std::to_string(quant+1),
205 std::vector<TangibleAxis> variables;
208 variables.emplace_back( EnergyAxis ,
"MuonEnergy" ,
"Reconstructed Muon Energy (GeV)" );
209 variables.emplace_back( MuonEPerHit ,
"MuEPerHit" ,
"Muon energy per hit in track" );
210 variables.emplace_back( HadEAxis ,
"HadroEnergy",
"Hadronic Energy (GeV)" );
211 variables.emplace_back( HadEFracAxis,
"HadEFrac" ,
"Hadronic Energy Fraction" );
212 variables.emplace_back( HadEPerHit ,
"HadEPerHit" ,
"Hadronic energy per hit in slice");
213 variables.emplace_back( LengthAxis ,
"KalTrackLen",
"Length of Primary Track (m)" );
214 variables.emplace_back( NHitsInSlc ,
"NHitsInSlc" ,
"Number of hits in the slice" );
215 variables.emplace_back( CosNuMiAxis ,
"TrkCosNuMi" ,
"Kalman Track Cos #theta_{NuMI}" );
216 variables.emplace_back( pTpAxis ,
"pTpNuMi" ,
"p_{T}/p" );
218 variables.emplace_back( HadEFracHRes,
"HadEFacHRes",
"Hadronic Energy Fraction" );
219 variables.emplace_back( NHitsKalTr ,
"NHitsKalTr" ,
"Number of hits in Kalman track" );
220 variables.emplace_back( CosNuMiAxis ,
"TrkCosNuMi" ,
"Kalman Track Cos #theta_{NuMI}" );
221 variables.emplace_back( CVN2017Axis ,
"CVNNuMu2017",
"CVN NuMu 2017 score" );
222 variables.emplace_back( CVN2018Axis ,
"CVNNuMu2018",
"CVN NuMu score" );
223 variables.emplace_back( ReMIDAxis ,
"ReMIdScore" ,
"ReMId score" );
224 variables.emplace_back( ReMIDFrac ,
"ReMIdFrac" ,
"RemID fraction of planes" );
225 variables.emplace_back( ReMIDScat ,
"ReMIdScatt" ,
"ReMId Scattering angle" );
226 variables.emplace_back( ReMIDdEdx ,
"ReMIddEdx" ,
"ReMId dE/dx" );
235 std::vector<DataMCPair> pairs;
240 pairs.reserve(selections.size() * variables.size());
241 for(
const auto& sel:selections){
242 for(
const auto& variable:variables){
243 pairs.emplace_back(sel, variable, loaderData, loaderMC, systs, MyWeight, !
kIsNumuCC, WrongSign );
252 std::vector<TCanvas*>
cans;
257 for(
const auto& pair:pairs) {
259 std::string Name_POT = pair.ShortName()+
"_allSysts_pot";
260 std::string Name_Area = pair.ShortName()+
"_allSysts_area";
261 std::cout <<
"\nNow looking pairs, " << pair.ShortName() <<
" - POT " << Name_POT <<
", Area " << Name_Area <<
", it has purity " << pair.Purity() <<
std::endl;
264 cans.push_back(
new TCanvas( Name_POT.c_str(),Name_POT.c_str() ) );
265 pair.OverlayDataMCSyst();
268 "Full systematic band shown. ",
269 "Near detector data/MC comparison for numu Ana2018." );
272 cans.push_back(
new TCanvas( Name_Area.c_str(), Name_Area.c_str() ) );
273 pair.OverlayDataMCSystNorm();
276 std::string(
"Each systematically shifted histogram (both up and down) ")
277 +
std::string(
"is normalized to the area of the MC distribution, ")
279 "Near detector data/MC comparison for numu Ana2018." );
284 std::string OutFileName = fnameODir +
"NearDet_DataMC_Comp_Systs.root";
285 TFile *
OutFile =
new TFile( OutFileName.c_str(),
"RECREATE" );
288 for(
const auto&
can:cans) {
293 if (!IsFHC) Str =
"Antineutrino beam";
294 TLatex* CornLab =
new TLatex(.15, .93, Str.c_str());
295 CornLab->SetTextColor(kGray+1);
297 CornLab->SetTextSize (2/30.);
298 CornLab->SetTextAlign(11);
302 can->SaveAs ( CanNa.c_str() );
303 can->Write (
can->GetName() );
const Cut kIsNumuCC(CCFlavSel(14, 14))
Select CC .
Cuts and Vars for the 2020 FD DiF Study.
const Binning kRemidBinning
Binning for plotting remid attractively.
const Cut kNumuContainND2017([](const caf::SRProxy *sr){if(!sr->vtx.elastic.IsValid) return false;for(unsigned int i=0;i< sr->vtx.elastic.fuzzyk.nshwlid;++i){TVector3 start=sr->vtx.elastic.fuzzyk.png[i].shwlid.start;TVector3 stop=sr->vtx.elastic.fuzzyk.png[i].shwlid.stop;if(std::min(start.X(), stop.X())< -180.0) return false;if(std::max(start.X(), stop.X()) > 180.0) return false;if(std::min(start.Y(), stop.Y())< -180.0) return false;if(std::max(start.Y(), stop.Y()) > 180.0) return false;if(std::min(start.Z(), stop.Z())< 20.0) return false;if(std::max(start.Z(), stop.Z()) > 1525.0) return false;}if(sr->trk.kalman.ntracks< 1) return false;for(unsigned int i=0;i< sr->trk.kalman.ntracks;++i){if(i==sr->trk.kalman.idxremid) continue;else if(sr->trk.kalman.tracks[i].start.Z() > 1275||sr->trk.kalman.tracks[i].stop.Z() > 1275) return false;}return(sr->trk.kalman.ntracks > sr->trk.kalman.idxremid &&sr->slc.firstplane > 1
&&sr->slc.lastplane< 212 &&sr->trk.kalman.tracks[0].start.Z()< 1100
&&(sr->trk.kalman.tracks[0].stop.Z()< 1275 ||sr->sel.contain.kalyposattrans< 55) &&sr->sel.contain.kalfwdcellnd > 5 &&sr->sel.contain.kalbakcellnd > 10);})
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 ...
const Cut kIsAntiNu([](const caf::SRProxy *sr){if(sr->mc.nnu==0) return false;assert(sr->mc.nnu==1);return sr->mc.nu[0].pdg< 0;})
Is this truly an antineutrino?
void SetSpillCut(const SpillCut &cut)
std::string pnfs2xrootd(std::string loc, bool unauth)
void MakeTextFile(std::string OutName, std::string OutDir)
const Var kTrkLength([](const caf::SRProxy *sr){if(sr->trk.kalman.ntracks< 1) return-5.f;return sr->trk.kalman.tracks[0].len/100;})
const Var kHadEPerNHit([](const caf::SRProxy *sr){if(sr->trk.kalman.ntracks< 1) return 0.0f;int nHit=sr->slc.nhit-sr->trk.kalman.tracks[0].nhit;if(nHit<=0) return 0.0f;float hadE=sr->energy.numu.hadcalE;return hadE/nHit;})
const Cut kIsNu([](const caf::SRProxy *sr){if(sr->mc.nnu==0) return false;assert(sr->mc.nnu==1);return sr->mc.nu[0].pdg > 0;})
std::vector< const ISyst * > getAllDirectNumuSysts2018()
const Var kNumuMuonPtP([](const caf::SRProxy *sr){if(sr->trk.kalman.ntracks > 0 &&sr->trk.kalman.idxremid!=999){if(sr->hdr.det==1){double Zbeam=sr->trk.kalman.tracks[0].dir.Dot(beamDirND);double ptp=sqrt(1-Zbeam *Zbeam);return(float) ptp;}if(sr->hdr.det==2){double Zbeam=sr->trk.kalman.tracks[0].dir.Dot(beamDirFD);double ptp=sqrt(1-Zbeam *Zbeam);return(float) ptp;}}return-5.f;})
virtual void Go() override
Load all the registered spectra.
const Cut kNumuPID2017([](const caf::SRProxy *sr){return(sr->sel.remid.pid > 0.5 &&sr->sel.cvn.numuid > 0.5);})
const HistAxis kNumuCCOptimisedAxis("Reconstructed Neutrino Energy (GeV)", kNumuCCEOptimisedBinning, kCCE)
HistAxis that implements optimised numuCCE from L. Vinton. See docdb 16332. This was close to 'custC'...
Collaborates with Spectrum and OscillatableSpectrum to fill spectra from CAF files.
const Cut kNumuPID2018([](const caf::SRProxy *sr){std::cout<< "ERROR::kNumuPID2018, cutting on both cvnProd3Train and cvn2017."<< " Neither branch exists anymore. Returning False."<< std::endl;abort();return false;})
void ND_DataMC_Comp_Systs(bool fhc=true, bool Set1Vars=true, bool IsTest=false, bool IsOnGrid=true, bool Weight18=true, bool Prod4Cut=true, bool FullCuts=false)
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...
const unsigned int NumBins
const SpillCut kStandardSpillCuts
Apply this unless you're doing something special.
std::vector< const ISyst * > getAllDirectNumuSysts2017()
Get a vector of all the numu group systs.
std::string to_string(ModuleType mt)
static Binning Simple(int n, double lo, double hi, const std::vector< std::string > &labels={})
const Cut kNoCut
The simplest possible cut: pass everything, used as a default.
const Var kCosNumi([](const caf::SRProxy *sr){if(sr->trk.kalman.ntracks > 0 &&sr->trk.kalman.idxremid!=999){if(sr->hdr.det==1){return sr->trk.kalman.tracks[0].dir.Dot(beamDirND);}if(sr->hdr.det==2){return sr->trk.kalman.tracks[0].dir.Dot(beamDirFD);}}return-5.f;})
const Var kTrkNhits([](const caf::SRProxy *sr){if(sr->trk.kalman.ntracks< 1) return 65535;return int(sr->trk.kalman.tracks[0].nhit);})
const Var kXSecCVWgt2018RPAFix_noDIS
void WriteBlurb(const DataMCPair &pair, std::string dir, std::string name, std::string extra="", std::string beginning="")
std::vector< std::string > CutNames