make_xsec_tuning_hists_mp5.C
Go to the documentation of this file.
1 #include "CAFAna/Core/Binning.h"
2 #include "CAFAna/Core/Spectrum.h"
3 #include "CAFAna/Core/SpectrumLoader.h"
4 #include "CAFAna/Core/Utilities.h"
5 #include "CAFAna/Core/HistAxis.h"
6 #include "CAFAna/Core/Loaders.h"
7 #include "CAFAna/Cuts/Cuts.h"
8 #include "CAFAna/Cuts/SpillCuts.h"
9 #include "3FlavorAna/Cuts/NumuCuts.h"
10 #include "3FlavorAna/Cuts/NumuCuts2017.h"
11 #include "3FlavorAna/Cuts/NumuCuts2018.h"
12 #include "3FlavorAna/Cuts/QuantileCuts.h"
13 #include "CAFAna/Vars/Vars.h"
14 #include "3FlavorAna/Vars/NumuVars.h"
15 #include "CAFAna/Vars/TruthVars.h"
16 #include "CAFAna/Vars/GenieWeights.h"
17 #include "CAFAna/Vars/PPFXWeights.h"
18 #include "3FlavorAna/Vars/HistAxes.h"
19 //#include "CAFAna/Vars/XsecTunes.h"
20 #include "3FlavorAna/Systs/NumuSysts.h"
21 #include "CAFAna/Systs/XSecSysts.h"
22 //#include "MCReweight/func/HistWrapper.h"
23 #include "CAFAna/Vars/XsecTunes.h"
24 #include "StandardRecord/Proxy/SRProxy.h"
25 
26 #include "3FlavorAna/Cuts/NumuCuts2020.h"
27 
28 #include "TTree.h"
29 #include "TH2.h"
30 
31 #include <fstream>
32 
33 using namespace ana;
34 using namespace caf;
35 
36 const Var kMinervaMECGaussEnh( []( const caf::SRProxy* sr )
37 {
38  if ( !kIsNumuCC( sr ) || sr->mc.nu[0].mode != caf::kMEC ) return 1.0;
39 
40  double q0 = kTrueQ0( sr );
41  double q3 = kTrueQ3( sr );
42 
43  double norm = 10.58;
44  double mean_q0 = 0.254;
45  double mean_q3 = 0.508;
46  double sigma_q0 = 0.057;
47  double sigma_q3 = 0.129;
48  double corr = 0.875;
49 
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 );
52 
53  return 1.0 + norm * exp( -0.5 * z / ( 1 - corr * corr ) );
54 });
55 
56 const Var kDoubleGaussEnhMEC( []( const caf::SRProxy* sr )
57 {
58  if ( !kIsNumuCC( sr ) || sr->mc.nu[0].mode != caf::kMEC ) return 1.0;
59 
60  double q0 = kTrueQ0( sr );
61  double q3 = kTrueQ3( sr );
62 
63  // FHC double gauss MEC fit from DocDB 42265-v2
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;
69  double corr_1 = 0.91;
70 
71  double norm_2 = 38.8;
72  double mean_q0_2 = 0.036;
73  double mean_q3_2 = 0.46;
74  double sigma_q0_2 = 0.042;
75  double sigma_q3_2 = 0.233;
76  double corr_2 = 0.71;
77 
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 );
80 
81  double z_2 = pow( ( q0 - mean_q0_2 ) / sigma_q0_2, 2 ) + pow( ( q3 - mean_q3_2 ) / sigma_q3_2, 2 )
82  - 2 * corr_2 * ( q0 - mean_q0_2 ) * ( q3 - mean_q3_2 ) / ( sigma_q0_2 * sigma_q3_2 );
83 
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 ) );
85 
86 });
87 
88 const Var kMINOSLowQ2RESSupp( []( const caf::SRProxy * sr )
89 {
90  if ( !kIsRes( sr ) ) return 1.0;
91  if ( sr->mc.nu[0].tgtA == 1 ) return 1.0; // Exclude hydrogen
92  double Q2 = kTrueQ2( sr );
93  return 1.010 / ( 1 + exp( 1 - sqrt( Q2 ) / 0.156 ) ); // From PhysRevD.91.012005
94 });
95 
96 const Var kMINERvALowQ2RESSupp( []( const caf::SRProxy * sr )
97 {
98  if ( !kIsRes( sr ) ) return 1.0;
99  if ( sr->mc.nu[0].tgtA == 1 ) return 1.0; // Exclude hydrogen
100 
101  // arXiv:1903.01558
102  // Joint FrAbs Tune
103  double R1 = 0.32;
104  double R2 = 0.50;
105 
106  double x1 = 0.0;
107  double x2 = 0.35;
108  double x3 = 0.7;
109 
110  double Q2 = kTrueQ2( sr );
111  if ( Q2 > x3 ) return 1.0;
112 
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 );
115  if ( weight < 0 )
116  {
117  std::cout << "kMINERvALowQ2RESSupp = " << weight << std::endl;
118  std::cout << "Q2 = " << Q2 << std::endl;
119  weight = 0;
120  }
121  return weight;
122 });
123 
124 const Var kVarFermiP( []( const caf::SRProxy* sr )
125 {
126  if ( sr->mc.nnu != 1 ) return 0.;
127 
128  int A = sr->mc.nu[0].tgtA;
129  int Z = sr->mc.nu[0].tgtZ;
130  int N = A-Z;
131 
132  if ( A <= 1 ) return 0.;
133 
134  // From FermiMomentumForIsoscalarNucleonParametrization in $GENIE/src/Physics/NuclearState/NuclearUtils.cxx
135  // Compute Fermi momentum for isoscalar nucleon (in GeV)
136  double x = 1.0 / A;
137  double P_Fermi = 0.27+x*(-1.12887857491+x*(9.72670908033-39.53095724456*x));
138 
139  // From $GENIE/src/Physics/Resonance/XSection/ReinSehgalRESPXSec.cxx
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); }
143 
144  return P_Fermi;
145 });
146 
147 const Var kVarRESPauliBlock( []( const caf::SRProxy* sr )
148 {
149 
150  if ( !kIsRes( sr ) ) return 1.0;
151 
152  int hitnuc = sr->mc.nu[0].hitnuc;
153  if ( hitnuc != 2212 && hitnuc != 2112 )
154  {
155  std::cout << "hitnuc = " << hitnuc << std::endl;
156  return 1.0;
157  }
158 
159  int A = sr->mc.nu[0].tgtA;
160  if ( A < 6 ) return 1.0; // From $GENIE/src/Physics/Resonance/XSection/ReinSehgalRESPXSec.cxx
161 
162  double P_Fermi = kVarFermiP( sr );
163 
164  // From $GENIE/src/Framework/Conventions/Constants.h
165  static const double kPionMass = 1.3957018e-01; // GeV
166  static const double kPionMass2 = TMath::Power(kPionMass,2); // GeV^2
167  static const double kProtonMass = 9.38272081e-01; // GeV
168  static const double kNeutronMass = 9.39565413e-01; // GeV
169 
170  //bool is_p = hitnuc == 2212;
171  //double Mnuc = is_p ? kProtonMass : kNeutronMass;
172  double Mnuc = hitnuc == 2212 ? kProtonMass : kNeutronMass;
173  double Mnuc2 = TMath::Power(Mnuc, 2);
174 
175  double W = kTrueW( sr );
176  if ( !( W > 0 ) )
177  {
178  std::cout << "W = " << W << std::endl;
179  return 1.0;
180  }
181 
182  double W2 = TMath::Power(W,2);
183  double Q2 = kTrueQ2( sr );
184  double k = 0.5 * (W2 - Mnuc2)/Mnuc;
185 
186  // From $GENIE/src/Physics/Resonance/XSection/ReinSehgalRESPXSec.cxx
187  double FactorPauli_RES = 1.0;
188 
189  double k0 = 0., q = 0., q0 = 0.;
190 
191  if (P_Fermi > 0.)
192  {
193  k0 = (W2-Mnuc2-Q2)/(2*W);
194  k = TMath::Sqrt(k0*k0+Q2);
195  q0 = (W2-Mnuc2+kPionMass2)/(2*W);
196  q = TMath::Sqrt(q0*q0-kPionMass2);
197  }
198 
199  if (2*P_Fermi < k-q)
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);
205 
206  return FactorPauli_RES;
207 
208 });
209 
210 
211 
212 /*
213 const Var kHitNucBin( []( const caf::SRProxy * sr )
214 {
215  if ( !kIsDytmanMEC( sr ) ) return -1.0;
216  else if ( kHitNuc( sr ) == 2000000200 ) return 0.0;
217  else if ( kHitNuc( sr ) == 2000000201 ) return 1.0;
218  else if ( kHitNuc( sr ) == 2000000202 ) return 2.0;
219  else return 3.0;
220 });
221 */
222 
223 const Cut kIsOther = !kIsNumuCC || !( kIsQE || kIsRes || kIsDIS || kIsDytmanMEC );
224 
225 const Cut kIsNumuCCMEC = kIsNumuCC && kIsDytmanMEC && !kIsAntiNu;
226 const Cut kIsNumubarCCMEC = kIsNumuCC && kIsDytmanMEC && kIsAntiNu;
227 
228 const Cut kIsMuonNeutrino( []( const caf::SRProxy* sr )
229 {
230  return sr->hdr.ismc && sr->mc.nnu == 1 && sr->mc.nu[0].pdg == 14;
231 });
232 
233 const Cut kIsMuonAntiNeutrino( []( const caf::SRProxy* sr )
234 {
235  return sr->hdr.ismc && sr->mc.nnu == 1 && sr->mc.nu[0].pdg == -14;
236 });
237 
238 const Cut kWSelection( []( const caf::SRProxy* sr )
239 {
240  double W = kRecoW( sr );
241  return W > 1.2 && W < 1.5;
242 });
243 
244 class Decomp
245 {
246 
247  public :
248 
249  Decomp(
250  std::string name,
251  bool isMC,
252  SpectrumLoaderBase& loader,
253  const HistAxis& axis,
254  const Cut& cut,
255  const Var& wgtMC = kUnweighted,
256  const SystShifts& shiftMC = kNoShift )
257  {
258  fName = name;
259 
260  if ( isMC )
261  {
262  fSpectra[ "TotMC" ] = std::make_unique<Spectrum>( loader, axis, cut , shiftMC, wgtMC );
263  fSpectra[ "QEL" ] = std::make_unique<Spectrum>( loader, axis, cut && kIsNumuCC && kIsQE , shiftMC, wgtMC );
264  fSpectra[ "MEC" ] = std::make_unique<Spectrum>( loader, axis, cut && kIsNumuCC && kIsDytmanMEC, shiftMC, wgtMC );
265  fSpectra[ "RES" ] = std::make_unique<Spectrum>( loader, axis, cut && kIsNumuCC && kIsRes , shiftMC, wgtMC );
266  fSpectra[ "DIS" ] = std::make_unique<Spectrum>( loader, axis, cut && kIsNumuCC && kIsDIS , shiftMC, wgtMC );
267  fSpectra[ "Other" ] = std::make_unique<Spectrum>( loader, axis, cut && kIsOther , shiftMC, wgtMC );
268 
269  fSpectra[ "RES_PauliBlock" ] = std::make_unique<Spectrum>( loader, axis, cut && kIsNumuCC && kIsRes , shiftMC, wgtMC * kVarRESPauliBlock );
270 
271  fSpectra[ "MEC_MINERvA" ] = std::make_unique<Spectrum>( loader, axis, cut && kIsNumuCC && kIsDytmanMEC, shiftMC, wgtMC * kMinervaMECGaussEnh );
272  fSpectra[ "MEC_DoubleGauss" ] = std::make_unique<Spectrum>( loader, axis, cut && kIsNumuCC && kIsDytmanMEC, shiftMC, wgtMC * kDoubleGaussEnhMEC );
273  fSpectra[ "MEC_DoubleGauss_Numu" ] = std::make_unique<Spectrum>( loader, axis, cut && kIsNumuCCMEC , shiftMC, wgtMC * kDoubleGaussEnhMEC );
274  fSpectra[ "MEC_DoubleGauss_Numubar" ] = std::make_unique<Spectrum>( loader, axis, cut && kIsNumubarCCMEC , shiftMC, wgtMC * kDoubleGaussEnhMEC );
275 
276  /*
277  fSpectra[ "TotMC_Nominal" ] = std::make_unique<Spectrum>( loader, axis, cut , shiftMC, kPPFXFluxCVWgt );
278  fSpectra[ "QEL_Nominal" ] = std::make_unique<Spectrum>( loader, axis, cut && kIsNumuCC && kIsQE , shiftMC, kPPFXFluxCVWgt );
279  fSpectra[ "MEC_Nominal" ] = std::make_unique<Spectrum>( loader, axis, cut && kIsNumuCC && kIsDytmanMEC, shiftMC, kPPFXFluxCVWgt );
280  fSpectra[ "RES_Nominal" ] = std::make_unique<Spectrum>( loader, axis, cut && kIsNumuCC && kIsRes , shiftMC, kPPFXFluxCVWgt );
281  fSpectra[ "DIS_Nominal" ] = std::make_unique<Spectrum>( loader, axis, cut && kIsNumuCC && kIsDIS , shiftMC, kPPFXFluxCVWgt );
282  fSpectra[ "Other_Nominal" ] = std::make_unique<Spectrum>( loader, axis, cut && kIsOther , shiftMC, kPPFXFluxCVWgt );
283  */
284 
285  }
286  else
287  {
288  fSpectra[ "Data" ] = std::make_unique<Spectrum>( loader, axis, cut );
289  }
290  }
291 
292  void Write() const
293  {
294  for ( const auto & s : fSpectra )
295  {
296  double pot = s.second.get()->POT();
297  std::string out_name = Form( "%s_%s", fName.c_str(), s.first.c_str() );
298  int ndim = s.second.get()->NDimensions();
299  if ( ndim == 1 )
300  {
301  s.second.get()->ToTH1( pot )->Write( out_name.c_str() );
302  }
303  else if ( ndim == 2 )
304  {
305  s.second.get()->ToTH2( pot )->Write( out_name.c_str() );
306  }
307  else throw std::runtime_error( Form( "Write operation for %d dimensions not defined", ndim ) );
308  //s.second.get()->ToTH1( s.second.get()->POT() )->Write( Form( "%s_%s", fName.c_str(), s.first.c_str() ) );
309  }
310  }
311 
312  double POT() const { return fSpectra.begin()->second.get()->POT(); }
313 
314  private :
315 
316  std::string fName;
317  std::map<std::string,std::unique_ptr<Spectrum>> fSpectra;
318 };
319 
320 void make_xsec_tuning_hists_mp5
321 (
322  const std::string beam,
323  const std::string type
324 )
325 {
326  if ( type != "data" && type != "mc" )
327  {
328  std::cout << "type must be \"data\" or \"mc\"" << std::endl;
329  return;
330  }
331  const bool isMC = ( type == "mc" );
332 
333  std::string defn = "";
334  if ( beam == "fhc" )
335  {
336  defn = isMC ?
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";
339  }
340  else if ( beam == "rhc" )
341  {
342  defn = isMC ?
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";
345  }
346  if ( defn.empty() ) throw std::runtime_error( "SAM Definition is empty" );
347 
348  //defn = Form( "defname: %s with limit 1", defn.c_str() );
349  std::cout << "defn = " << defn << std::endl;
350 
351  SpectrumLoader loader( defn );
352  loader.SetSpillCut( kStandardSpillCuts );
353 
354  std::map< std::string, HistAxis > axes;
355  axes.emplace( "RecoEhad" , HistAxis( "Reco E_{had} (GeV)", Binning::Simple( 50, 0., 2. ), kHadE ) );
356  axes.emplace( "RecoEhadVis" , HistAxis( "Visible E_{had} (GeV)", Binning::Simple( 50, 0., 1. ), kNumuHadVisE ) );
357  axes.emplace( "RecoEmu" , HistAxis( "Reco E_{#mu} (GeV)", Binning::Simple( 60, 0., 3. ), kMuE ) );
358  axes.emplace( "RecoEnu" , HistAxis( "Reco E_{#nu} (GeV)", Binning::Simple( 40, 0., 4. ), kCCE ) );
359  axes.emplace( "RecoQ2" , HistAxis( "Reco Q^{2} (GeV^{2})", Binning::Simple( 50, 0., 1. ), kRecoQ2 ) );
360  axes.emplace( "RecoQ3" , HistAxis( "Reco |#vec{q}| (GeV)", Binning::Simple( 40, 0., 2. ), kRecoQmag ) );
361  axes.emplace( "RecoW" , HistAxis( "Reco W (GeV)", Binning::Simple( 60, 0., 3. ), kRecoW ) );
362  axes.emplace( "RecoEhadFrac" , HistAxis( "Reco E_{had} Fraction", Binning::Simple( 50, 0., 1. ), kHadEFrac ) );
363  axes.emplace( "RecoCosThetaMu" , HistAxis( "Reco cos(#theta_{#mu})", Binning::Simple( 50, 0., 1. ), kCosNumi ) );
364  axes.emplace( "PrimTrkLen" , HistAxis( "Primary Track Length (m)", Binning::Simple( 60, 0., 15.), kTrkLength ) );
365 
366  if ( isMC )
367  {
368  //axes.emplace( "NucleonPair", HistAxis( "Nucleon Pair" , Binning::Simple( 5, -1, 4 ), kHitNucBin ) );
369  axes.emplace( "TrueQ2" , HistAxis( "True Q^{2} (GeV^{2})", Binning::Simple( 50, 0., 2. ), kTrueQ2 ) );
370  //axes.emplace( "TrueQ0" , HistAxis( "True q_{0} (GeV)" , Binning::Simple( 200, 0., 2. ), kTrueQ0 ) );
371  //axes.emplace( "TrueQ3" , HistAxis( "True |#vec{q}| (GeV)", Binning::Simple( 200, 0., 2. ), kTrueQ3 ) );
372  //axes.emplace( "TrueW" , HistAxis( "True W (GeV)" , Binning::Simple( 60, 0., 3. ), kTrueW ) );
373  //axes.emplace( "TrueY" , HistAxis( "True y" , Binning::Simple( 50, 0., 1. ), kTrueY ) );
374  }
375 
376  axes.emplace
377  (
378  "RecoEhadVis_vs_RecoQ3",
379  HistAxis( "Reco |#vec{q}| (GeV)", Binning::Simple( 40, 0., 2. ), kRecoQmag, "Visible E_{had} (GeV)", Binning::Simple( 50, 0., 1. ), kNumuHadVisE )
380  );
381 
382  if ( isMC )
383  {
384 
385  axes.emplace
386  (
387  "TrueQ0_vs_TrueQ3",
388  HistAxis( "True |#vec{q}| (GeV)", Binning::Simple( 80, 0., 2. ), kTrueQ3, "True q_{0} (GeV)", Binning::Simple( 80,0, 2. ), kTrueQ0 )
389  );
390 
391  }
392 
393  // Hadronic energy fraction quantiles
394 
395  std::string InDirQuant = "/pnfs/nova/persistent/analysis/numu/Ana2018/provisional/quantiles";
396  std::string InFileQuant = InDirQuant + "/quantiles__" + beam + "_full__numu2018.root";
397  std::cout << "\n\n --- get quantile cuts from file ---" << std::endl;
398  TFile* inFile = TFile::Open( pnfs2xrootd(InFileQuant).c_str() );
399  TH2* FDSpec2D = (TH2*)inFile->FindObjectAny("FDSpec2D");
400  std::vector<Cut> HadEFracQuantCuts = QuantileCutsFromTH2(FDSpec2D, kNumuCCOptimisedAxis, kHadEFracAxis, 4);
401  inFile->Close();
402  delete inFile;
403 
404  const Cut kSelCuts = kNumuCutNDMiniprod5;
405  const Var kWeightCV = kPPFXFluxCVWgt;
406 
407  std::vector<Decomp> decomps;
408  decomps.reserve( 100 );
409 
410  for ( const auto & a : axes )
411  {
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 );
414 
415  if ( a.first != "RecoEhadFrac" )
416  {
417  for ( unsigned int icut = 0; icut < HadEFracQuantCuts.size(); ++icut )
418  {
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 );
421  }
422  }
423  }
424 
425  loader.Go();
426 
427  double pot = decomps.front().POT();
428 
429  TTree potTree( "pot", "pot" );
430  potTree.Branch( Form( "pot_%s", type.c_str() ), &pot, Form( "pot_%s/D", type.c_str() ) );
431  potTree.Fill();
432 
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" );
435 
436  potTree.Write();
437 
438  for ( const auto & d : decomps ) { d.Write(); }
439 
440  outputFile.Close();
441 
442 }
kIsMuonNeutrino
const Cut kIsMuonNeutrino([](const caf::SRProxy *sr){return sr->hdr.ismc &&sr->mc.nnu==1 &&sr->mc.nu[0].pdg==14;})
ana::kHadE
const Var kHadE
Definition: NumuVars.h:23
ana::kIsQE
const Cut kIsQE
Definition: TruthCuts.cxx:104
name
const XML_Char * name
Definition: expat.h:151
ana::kIsNumuCC
const Cut kIsNumuCC(CCFlavSel(14, 14))
Select CC .
kMINERvALowQ2RESSupp
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;})
ana::kRecoQ2
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)
Definition: NumuVars.h:146
ana::HistAxis
_HistAxis< Var > HistAxis
Definition: HistAxis.h:103
ana
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
genie::utils::kinematics::Q2
double Q2(const Interaction *const i)
Definition: KineUtils.cxx:991
ana::kPPFXFluxCVWgt
const Var kPPFXFluxCVWgt
Definition: PPFXWeights.h:16
ana::weight
const Var weight
Definition: PIDCutTuningDefs.h:96
ana::kNeutronMass
const Var kNeutronMass([](const caf::SRProxy *sr){return NeutronMass();})
ana::SystShifts
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
x1
Float_t x1[n_points_granero]
Definition: compare.C:5
caf::Proxy< caf::StandardRecord >::hdr
caf::Proxy< caf::SRHeader > hdr
Definition: SRProxy.h:2096
ana::kRecoQmag
const Var kRecoQmag([](const caf::SRProxy *sr){return sqrt(kRecoQ2(sr)+util::sqr(kHadE(sr)));})
Reconstructed three-momentum transfer.
Definition: NumuVars.h:149
ana::kIsRes
const Cut kIsRes
Definition: TruthCuts.cxx:111
caf::Proxy< caf::StandardRecord >
Proxy for caf::StandardRecord.
Definition: SRProxy.h:2085
std::sqrt
T sqrt(T number)
Definition: d0nt_math.hpp:156
caf::Proxy< caf::SRTruthBranch >::nu
caf::Proxy< std::vector< caf::SRNeutrino > > nu
Definition: SRProxy.h:574
cet::pow
constexpr T pow(T x)
Definition: pow.h:75
ana::kHadEFracAxis
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 ...
Definition: HistAxes.h:30
Decomp::POT
double POT() const
Definition: make_xsec_tuning_hists_mp5.C:312
submit_syst.x2
x2
Definition: submit_syst.py:250
ana::kIsAntiNu
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?
Definition: TruthCuts.h:53
caf::Proxy< caf::SRTruthBranch >::nnu
caf::Proxy< short int > nnu
Definition: SRProxy.h:573
ana::kIsDIS
const Cut kIsDIS
Definition: TruthCuts.cxx:118
calib::corr
double corr(TGraph *g, double thresh)
Definition: TestParticleCorrections_module.cc:110
kIsMuonAntiNeutrino
const Cut kIsMuonAntiNeutrino([](const caf::SRProxy *sr){return sr->hdr.ismc &&sr->mc.nnu==1 &&sr->mc.nu[0].pdg==-14;})
ana::SpectrumLoaderBase::SetSpillCut
void SetSpillCut(const SpillCut &cut)
Definition: SpectrumLoaderBase.cxx:64
make_xsec_tuning_hists_mp5
void make_xsec_tuning_hists_mp5(const std::string beam, const std::string type)
Definition: make_xsec_tuning_hists_mp5.C:321
ana::pnfs2xrootd
std::string pnfs2xrootd(std::string loc, bool unauth)
Definition: UtilsExt.cxx:237
ana::kTrueQ0
const Var kTrueQ0
Definition: TruthVars.h:32
kDoubleGaussEnhMEC
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));})
kMINOSLowQ2RESSupp
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));})
ana::kTrkLength
const Var kTrkLength([](const caf::SRProxy *sr){if(sr->trk.kalman.ntracks< 1) return-5.f;return sr->trk.kalman.tracks[0].len/100;})
Definition: NumuVars.h:65
ana::kNumuHadVisE
const Var kNumuHadVisE([](const caf::SRProxy *sr){return kNumuHadCalE(sr)+kNumuHadTrkE(sr);})
Definition: NumuVars.h:124
ana::kTrueQ2
const Var kTrueQ2
Definition: TruthVars.h:27
ana::kIsDytmanMEC
const Cut kIsDytmanMEC
Definition: TruthCuts.cxx:187
Z
Float_t Z
Definition: plot.C:38
inFile
ifstream inFile
Definition: AnaPlotMaker.h:34
s
const XML_Char * s
Definition: expat.h:262
submit_syst.x
x
Definition: submit_syst.py:249
ana::kTrueQ3
const Var kTrueQ3
Definition: TruthVars.h:38
kIsNumuCCMEC
const Cut kIsNumuCCMEC
Definition: make_xsec_tuning_hists_mp5.C:225
if
if(dump)
Definition: checkMeasured_Target.C:110
a
const double a
Definition: gUpMuFluxGen.cxx:164
ana::kHadEFrac
const Var kHadEFrac
Definition: NumuVars.h:24
stan::math::exp
fvar< T > exp(const fvar< T > &x)
Definition: exp.hpp:10
ana::kCCE
const Var kCCE
Definition: NumuVars.h:21
pot
#define pot
Definition: template_basic.C:49
d
Float_t d
Definition: plot.C:236
kWSelection
const Cut kWSelection([](const caf::SRProxy *sr){double W=kRecoW(sr);return W > 1.2 &&W< 1.5;})
R
#define R(x)
sr
caf::StandardRecord * sr
Definition: FillTrainingTrees.h:18
ana::SpectrumLoader::Go
virtual void Go() override
Load all the registered spectra.
Definition: SpectrumLoader.cxx:66
kIsNumubarCCMEC
const Cut kIsNumubarCCMEC
Definition: make_xsec_tuning_hists_mp5.C:226
Decomp::Write
void Write() const
Definition: make_xsec_tuning_hists_mp5.C:292
demo0.loader
loader
Definition: demo0.py:10
test.z
z
Definition: test.py:28
ana::kNumuCCOptimisedAxis
const HistAxis kNumuCCOptimisedAxis("Reconstructed Neutrino Energy (GeV)", kNumuCCEOptimisedBinning, kCCE)
HistAxis that implements optimised numuCCE from L. Vinton. See docdb 16332. This was close to &#39;custC&#39;...
Definition: HistAxes.h:24
ana::kNoShift
const SystShifts kNoShift
Definition: SystShifts.cxx:21
om::cout
OStream cout
Definition: OStream.cxx:6
submit_concat_project.out_name
out_name
Definition: submit_concat_project.py:352
ana::SpectrumLoaderBase
Base class for the various types of spectrum loader.
Definition: SpectrumLoaderBase.h:63
caf::Proxy< caf::StandardRecord >::mc
caf::Proxy< caf::SRTruthBranch > mc
Definition: SRProxy.h:2097
ana::kNumuCutNDMiniprod5
const Cut kNumuCutNDMiniprod5
Definition: NumuCuts2020.h:52
ana::kRecoW
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)
Definition: NumuVars.h:152
xml_schema::string
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
cut
const Cut cut
Definition: exporter_fd.C:30
genie::units::A
static const double A
Definition: Units.h:82
genie::constants::kPionMass
static const double kPionMass
Definition: Constants.h:74
norm
Float_t norm
Definition: comparison_ascii.C:71
ana::_HistAxis< Var >
ana::SpectrumLoader
Collaborates with Spectrum and OscillatableSpectrum to fill spectra from CAF files.
Definition: SpectrumLoader.h:21
ana::Sqrt
Var Sqrt(const Var &v)
Use to take sqrt of a var.
Definition: Var.cxx:326
caf::Proxy< caf::SRHeader >::ismc
caf::Proxy< bool > ismc
Definition: SRProxy.h:241
kVarRESPauliBlock
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;})
ana::QuantileCutsFromTH2
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...
Definition: QuantileCuts.cxx:198
kIsOther
const Cut kIsOther
Definition: make_xsec_tuning_hists_mp5.C:223
file_size_ana.axes
axes
Definition: file_size_ana.py:282
ana::kStandardSpillCuts
const SpillCut kStandardSpillCuts
Apply this unless you&#39;re doing something special.
Definition: SpillCuts.h:49
kVarFermiP
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;})
caf
This module creates Common Analysis Files.
Definition: FileReducer.h:10
ana::_Cut< caf::SRProxy >
Decomp::Decomp
Decomp(std::string name, bool isMC, SpectrumLoaderBase &loader, const HistAxis &axis, const Cut &cut, const Var &wgtMC=kUnweighted, const SystShifts &shiftMC=kNoShift)
Definition: make_xsec_tuning_hists_mp5.C:249
ana::_Var< caf::SRProxy >
ana::kMuE
const Var kMuE
Definition: NumuVars.h:22
Decomp::fSpectra
std::map< std::string, std::unique_ptr< Spectrum > > fSpectra
Definition: make_xsec_tuning_hists_mp5.C:317
ana::Binning::Simple
static Binning Simple(int n, double lo, double hi, const std::vector< std::string > &labels={})
Definition: Binning.cxx:114
ana::kTrueW
const Var kTrueW
Definition: TruthVars.h:22
W
#define W(x)
kMinervaMECGaussEnh
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));})
ana::kUnweighted
const Var kUnweighted
The simplest possible Var, always 1. Used as a default weight.
Definition: Var.h:96
genie::constants::kProtonMass
static const double kProtonMass
Definition: Constants.h:76
ana::kCosNumi
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;})
Definition: NumuVars.h:43
genie::constants::kPionMass2
static const double kPionMass2
Definition: Constants.h:87
Generated on Tue Oct 20 2020 05:11:57 for NOvA by   doxygen 1.8.11