1 #include "CAFAna/Core/Binning.h" 5 #include "CAFAna/Core/HistAxis.h" 50 double z =
pow( ( q0 - mean_q0 ) / sigma_q0, 2 ) +
pow( ( q3 - mean_q3 ) / sigma_q3, 2 )
51 - 2 * corr * ( q0 -
mean_q0 ) * ( q3 - mean_q3 ) / ( sigma_q0 *
sigma_q3 );
53 return 1.0 + norm *
exp( -0.5 * z / ( 1 - corr * corr ) );
64 double norm_1 = 12.78;
65 double mean_q0_1 = 0.33;
66 double mean_q3_1 = 0.80;
67 double sigma_q0_1 = 0.10;
68 double sigma_q3_1 = 0.29;
78 double z_1 =
pow( ( q0 - mean_q0_1 ) / sigma_q0_1, 2 ) +
pow( ( q3 - mean_q3_1 ) / sigma_q3_1, 2 )
79 - 2 * corr_1 * ( q0 - mean_q0_1 ) * ( q3 - mean_q3_1 ) / ( sigma_q0_1 * sigma_q3_1 );
81 double z_2 =
pow( ( q0 - mean_q0_2 ) / sigma_q0_2, 2 ) +
pow( ( q3 - mean_q3_2 ) / sigma_q3_2, 2 )
84 return 1.0 + norm_1 *
exp( -0.5 * z_1 / ( 1 - corr_1 * corr_1 ) ) + norm_2 *
exp( -0.5 * z_2 / ( 1 - corr_2 * corr_2 ) );
90 if ( !
kIsRes( sr ) )
return 1.0;
91 if ( sr->
mc.
nu[0].tgtA == 1 )
return 1.0;
93 return 1.010 / ( 1 +
exp( 1 -
sqrt(
Q2 ) / 0.156 ) );
98 if ( !
kIsRes( sr ) )
return 1.0;
99 if ( sr->
mc.
nu[0].tgtA == 1 )
return 1.0;
111 if (
Q2 > x3 )
return 1.0;
113 double R = R2 * (
Q2 -
x1 ) * (
Q2 - x3 ) / ( (
x2 -
x1 ) * (
x2 - x3 ) ) + (
Q2 -
x1 ) * (
Q2 -
x2 ) / ( ( x3 -
x1 ) * ( x3 -
x2 ) );
114 double weight = 1 - ( 1 - R1 ) *
pow( 1 - R, 2 );
126 if ( sr->
mc.
nnu != 1 )
return 0.;
128 int A = sr->
mc.
nu[0].tgtA;
129 int Z = sr->
mc.
nu[0].tgtZ;
132 if (
A <= 1 )
return 0.;
137 double P_Fermi = 0.27+x*(-1.12887857491+x*(9.72670908033-39.53095724456*
x));
140 bool is_p = sr->
mc.
nu[0].hitnuc == 2212;
141 if(is_p) { P_Fermi *= TMath::Power( 2.*
Z/
A, 1./3); }
142 else { P_Fermi *= TMath::Power( 2.*N/
A, 1./3); }
150 if ( !
kIsRes( sr ) )
return 1.0;
152 int hitnuc = sr->
mc.
nu[0].hitnuc;
153 if ( hitnuc != 2212 && hitnuc != 2112 )
159 int A = sr->
mc.
nu[0].tgtA;
160 if ( A < 6 )
return 1.0;
165 static const double kPionMass = 1.3957018e-01;
166 static const double kPionMass2 = TMath::Power(kPionMass,2);
172 double Mnuc = hitnuc == 2212 ? kProtonMass :
kNeutronMass;
173 double Mnuc2 = TMath::Power(Mnuc, 2);
182 double W2 = TMath::Power(W,2);
184 double k = 0.5 * (W2 - Mnuc2)/Mnuc;
187 double FactorPauli_RES = 1.0;
189 double k0 = 0., q = 0., q0 = 0.;
193 k0 = (W2-Mnuc2-
Q2)/(2*W);
200 FactorPauli_RES = 1.0;
201 if (2*P_Fermi >= k+q)
202 FactorPauli_RES = ((3*k*k+q*q)/(2*P_Fermi)-(5*TMath::Power(k,4)+TMath::Power(q,4)+10*k*k*q*q)/(40*TMath::Power(P_Fermi,3)))/(2*k);
203 if (2*P_Fermi >= k-q && 2*P_Fermi <= k+q)
204 FactorPauli_RES = ((q+k)*(q+k)-4*P_Fermi*P_Fermi/5-TMath::Power(k-q, 3)/(2*P_Fermi)+TMath::Power(k-q, 5)/(40*TMath::Power(P_Fermi, 3)))/(4*q*k);
206 return FactorPauli_RES;
241 return W > 1.2 && W < 1.5;
262 fSpectra[
"TotMC" ] = std::make_unique<Spectrum>(
loader,
axis,
cut , shiftMC, wgtMC );
267 fSpectra[
"Other" ] = std::make_unique<Spectrum>(
loader,
axis, cut &&
kIsOther , shiftMC, wgtMC );
288 fSpectra[
"Data" ] = std::make_unique<Spectrum>(
loader,
axis,
cut );
294 for (
const auto &
s : fSpectra )
296 double pot =
s.second.get()->POT();
298 int ndim =
s.second.get()->NDimensions();
301 s.second.get()->ToTH1( pot )->Write( out_name.c_str() );
303 else if ( ndim == 2 )
305 s.second.get()->ToTH2( pot )->Write( out_name.c_str() );
307 else throw std::runtime_error( Form(
"Write operation for %d dimensions not defined", ndim ) );
312 double POT()
const {
return fSpectra.begin()->second.get()->POT(); }
317 std::map<std::string,std::unique_ptr<Spectrum>>
fSpectra;
326 if ( type !=
"data" && type !=
"mc" )
331 const bool isMC = ( type ==
"mc" );
337 "prod_caf_R19-02-23-miniprod5.m_nd_genie_G1810b0211a_nonswap_fhc_nova_v08_full_v1" :
338 "prod_caf_R19-02-23-miniprod5.i_nd_numi_fhc_full_v1_goodruns_stride10";
340 else if ( beam ==
"rhc" )
343 "prod_caf_R19-02-23-miniprod5.m_nd_genie_G1810b0211a_nonswap_rhc_nova_v08_full_v1_nogeniereweight" :
344 "prod_caf_R19-02-23-miniprod5.i_nd_numi_rhc_full_v1_goodruns_stride20";
346 if ( defn.empty() )
throw std::runtime_error(
"SAM Definition is empty" );
354 std::map< std::string, HistAxis >
axes;
378 "RecoEhadVis_vs_RecoQ3",
395 std::string InDirQuant =
"/pnfs/nova/persistent/analysis/numu/Ana2018/provisional/quantiles";
396 std::string InFileQuant = InDirQuant +
"/quantiles__" + beam +
"_full__numu2018.root";
399 TH2* FDSpec2D = (TH2*)inFile->FindObjectAny(
"FDSpec2D");
407 std::vector<Decomp> decomps;
408 decomps.reserve( 100 );
410 for (
const auto &
a : axes )
412 decomps.emplace_back(
a.first, isMC, loader,
a.second, kSelCuts, kWeightCV,
kNoShift );
413 decomps.emplace_back(
a.first +
"_WSelected", isMC, loader,
a.second, kSelCuts &&
kWSelection, kWeightCV,
kNoShift );
415 if (
a.first !=
"RecoEhadFrac" )
417 for (
unsigned int icut = 0; icut < HadEFracQuantCuts.size(); ++icut )
419 std::string dist_quant_name =
a.first + Form(
"_EhadQuantile%d", icut + 1 );
420 decomps.emplace_back( dist_quant_name, isMC, loader,
a.second, kSelCuts && HadEFracQuantCuts[ icut ], kWeightCV,
kNoShift );
427 double pot = decomps.front().POT();
429 TTree potTree(
"pot",
"pot" );
430 potTree.Branch( Form(
"pot_%s", type.c_str() ), &pot, Form(
"pot_%s/D", type.c_str() ) );
433 std::string out_file_name = Form(
"hists_xec_tuning_mp5_%s_%s", beam.c_str(), type.c_str() );
434 TFile outputFile( Form(
"%s.root", out_file_name.c_str() ),
"RECREATE" );
438 for (
const auto &
d : decomps ) {
d.Write(); }
const Cut kIsMuonNeutrino([](const caf::SRProxy *sr){return sr->hdr.ismc &&sr->mc.nnu==1 &&sr->mc.nu[0].pdg==14;})
const Cut kIsNumuCC(CCFlavSel(14, 14))
Select CC .
const Var kMINERvALowQ2RESSupp([](const caf::SRProxy *sr){if(!kIsRes(sr)) return 1.0;if(sr->mc.nu[0].tgtA==1) return 1.0; double R1=0.32;double R2=0.50;double x1=0.0;double x2=0.35;double x3=0.7;double Q2=kTrueQ2(sr);if(Q2 > x3) return 1.0;double R=R2 *(Q2-x1)*(Q2-x3)/((x2-x1)*(x2-x3))+(Q2-x1)*(Q2-x2)/((x3-x1)*(x3-x2));double weight=1-(1-R1)*pow(1-R, 2);if(weight< 0){std::cout<< "kMINERvALowQ2RESSupp = "<< weight<< std::endl;std::cout<< "Q2 = "<< Q2<< std::endl;weight=0;}return weight;})
const ana::Var kRecoQ2([](const caf::SRProxy *sr){const double M_mu_sqrd=util::sqr(0.1056);double E_mu=kMuE(sr);double p_mu=sqrt(util::sqr(E_mu)-M_mu_sqrd);return 2 *kCCE(sr)*(E_mu-p_mu *kCosNumi(sr))-M_mu_sqrd;})
Reconstructed four-momentum transfer invariant (Q^2)
_HistAxis< Var > HistAxis
Cuts and Vars for the 2020 FD DiF Study.
double Q2(const Interaction *const i)
const Var kNeutronMass([](const caf::SRProxy *sr){return NeutronMass();})
Simple record of shifts applied to systematic parameters.
Float_t x1[n_points_granero]
caf::Proxy< caf::SRHeader > hdr
const Var kRecoQmag([](const caf::SRProxy *sr){return sqrt(kRecoQ2(sr)+util::sqr(kHadE(sr)));})
Reconstructed three-momentum transfer.
Proxy for caf::StandardRecord.
caf::Proxy< std::vector< caf::SRNeutrino > > nu
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?
caf::Proxy< short int > nnu
const Cut kIsMuonAntiNeutrino([](const caf::SRProxy *sr){return sr->hdr.ismc &&sr->mc.nnu==1 &&sr->mc.nu[0].pdg==-14;})
void SetSpillCut(const SpillCut &cut)
void make_xsec_tuning_hists_mp5(const std::string beam, const std::string type)
std::string pnfs2xrootd(std::string loc, bool unauth)
const Var kDoubleGaussEnhMEC([](const caf::SRProxy *sr){if(!kIsNumuCC(sr)||sr->mc.nu[0].mode!=caf::kMEC) return 1.0;double q0=kTrueQ0(sr);double q3=kTrueQ3(sr);double norm_1=12.78;double mean_q0_1=0.33;double mean_q3_1=0.80;double sigma_q0_1=0.10;double sigma_q3_1=0.29;double corr_1=0.91;double norm_2=38.8;double mean_q0_2=0.036;double mean_q3_2=0.46;double sigma_q0_2=0.042;double sigma_q3_2=0.233;double corr_2=0.71;double z_1=pow((q0-mean_q0_1)/sigma_q0_1, 2)+pow((q3-mean_q3_1)/sigma_q3_1, 2)-2 *corr_1 *(q0-mean_q0_1)*(q3-mean_q3_1)/(sigma_q0_1 *sigma_q3_1);double z_2=pow((q0-mean_q0_2)/sigma_q0_2, 2)+pow((q3-mean_q3_2)/sigma_q3_2, 2)-2 *corr_2 *(q0-mean_q0_2)*(q3-mean_q3_2)/(sigma_q0_2 *sigma_q3_2);return 1.0+norm_1 *exp(-0.5 *z_1/(1-corr_1 *corr_1))+norm_2 *exp(-0.5 *z_2/(1-corr_2 *corr_2));})
const Var kMINOSLowQ2RESSupp([](const caf::SRProxy *sr){if(!kIsRes(sr)) return 1.0;if(sr->mc.nu[0].tgtA==1) return 1.0;double Q2=kTrueQ2(sr);return 1.010/(1+exp(1-sqrt(Q2)/0.156));})
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 kNumuHadVisE([](const caf::SRProxy *sr){return kNumuHadCalE(sr)+kNumuHadTrkE(sr);})
fvar< T > exp(const fvar< T > &x)
const Cut kWSelection([](const caf::SRProxy *sr){double W=kRecoW(sr);return W > 1.2 &&W< 1.5;})
virtual void Go() override
Load all the registered spectra.
const Cut kIsNumubarCCMEC
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'...
const SystShifts kNoShift
Base class for the various types of spectrum loader.
caf::Proxy< caf::SRTruthBranch > mc
const Cut kNumuCutNDMiniprod5
const Var kRecoW([](const caf::SRProxy *sr){const double M_p=0.938272;double WSq=M_p *M_p+2 *M_p *(kCCE(sr)-kMuE(sr))-kRecoQ2(sr);if(WSq > 0) return sqrt(WSq);else return-5.;})
Reconstructed invariant mass (W)
static const double kPionMass
Collaborates with Spectrum and OscillatableSpectrum to fill spectra from CAF files.
Var Sqrt(const Var &v)
Use to take sqrt of a var.
const Var kVarRESPauliBlock([](const caf::SRProxy *sr){if(!kIsRes(sr)) return 1.0;int hitnuc=sr->mc.nu[0].hitnuc;if(hitnuc!=2212 &&hitnuc!=2112){std::cout<< "hitnuc = "<< hitnuc<< std::endl;return 1.0;}int A=sr->mc.nu[0].tgtA;if(A< 6) return 1.0;double P_Fermi=kVarFermiP(sr);static const double kPionMass=1.3957018e-01;static const double kPionMass2=TMath::Power(kPionMass, 2);static const double kProtonMass=9.38272081e-01;static const double kNeutronMass=9.39565413e-01; double Mnuc=hitnuc==2212?kProtonMass:kNeutronMass;double Mnuc2=TMath::Power(Mnuc, 2);double W=kTrueW(sr);if(!(W > 0)){std::cout<< "W = "<< W<< std::endl;return 1.0;}double W2=TMath::Power(W, 2);double Q2=kTrueQ2(sr);double k=0.5 *(W2-Mnuc2)/Mnuc;double FactorPauli_RES=1.0;double k0=0., q=0., q0=0.;if(P_Fermi > 0.){k0=(W2-Mnuc2-Q2)/(2 *W);k=TMath::Sqrt(k0 *k0+Q2);q0=(W2-Mnuc2+kPionMass2)/(2 *W);q=TMath::Sqrt(q0 *q0-kPionMass2);}if(2 *P_Fermi< k-q) FactorPauli_RES=1.0;if(2 *P_Fermi >=k+q) FactorPauli_RES=((3 *k *k+q *q)/(2 *P_Fermi)-(5 *TMath::Power(k, 4)+TMath::Power(q, 4)+10 *k *k *q *q)/(40 *TMath::Power(P_Fermi, 3)))/(2 *k);if(2 *P_Fermi >=k-q &&2 *P_Fermi<=k+q) FactorPauli_RES=((q+k)*(q+k)-4 *P_Fermi *P_Fermi/5-TMath::Power(k-q, 3)/(2 *P_Fermi)+TMath::Power(k-q, 5)/(40 *TMath::Power(P_Fermi, 3)))/(4 *q *k);return FactorPauli_RES;})
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 SpillCut kStandardSpillCuts
Apply this unless you're doing something special.
const Var kVarFermiP([](const caf::SRProxy *sr){if(sr->mc.nnu!=1) return 0.;int A=sr->mc.nu[0].tgtA;int Z=sr->mc.nu[0].tgtZ;int N=A-Z;if(A<=1) return 0.; double x=1.0/A;double P_Fermi=0.27+x *(-1.12887857491+x *(9.72670908033-39.53095724456 *x));bool is_p=sr->mc.nu[0].hitnuc==2212;if(is_p){P_Fermi *=TMath::Power(2.*Z/A, 1./3);}else{P_Fermi *=TMath::Power(2.*N/A, 1./3);}return P_Fermi;})
This module creates Common Analysis Files.
Decomp(std::string name, bool isMC, SpectrumLoaderBase &loader, const HistAxis &axis, const Cut &cut, const Var &wgtMC=kUnweighted, const SystShifts &shiftMC=kNoShift)
std::map< std::string, std::unique_ptr< Spectrum > > fSpectra
static Binning Simple(int n, double lo, double hi, const std::vector< std::string > &labels={})
const Var kMinervaMECGaussEnh([](const caf::SRProxy *sr){if(!kIsNumuCC(sr)||sr->mc.nu[0].mode!=caf::kMEC) return 1.0;double q0=kTrueQ0(sr);double q3=kTrueQ3(sr);double norm=10.58;double mean_q0=0.254;double mean_q3=0.508;double sigma_q0=0.057;double sigma_q3=0.129;double corr=0.875;double z=pow((q0-mean_q0)/sigma_q0, 2)+pow((q3-mean_q3)/sigma_q3, 2)-2 *corr *(q0-mean_q0)*(q3-mean_q3)/(sigma_q0 *sigma_q3);return 1.0+norm *exp(-0.5 *z/(1-corr *corr));})
const Var kUnweighted
The simplest possible Var, always 1. Used as a default weight.
static const double kProtonMass
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;})
static const double kPionMass2