GiBUUWeights.cxx
Go to the documentation of this file.
2 
4 
6 
7 #include "TFile.h"
8 #include "TH1.h"
9 #include "TTree.h"
10 
11 #include <iostream>
12 
13 namespace ana
14 {
15  //----------------------------------------------------------------------
16  const NuTruthCut kIsGibuu_NT([](const caf::SRNeutrinoProxy* nu)
17  {
18  // Positive mode means it's GENIE (COH
19  // or nu-on-e or >30GeV or weird
20  // nucleus).
21  return nu->mode < 0;
22  });
23 
24  //----------------------------------------------------------------------
26  {
27  if(!kIsGibuu_NT(nu)) return 1.f;
28  return float(nu->genweight);
29  });
30 
31  //----------------------------------------------------------------------
33  {
34  for(int i = 0; i < 6; ++i)
35  for(int j = 0; j < 6; ++j)
36  delete fCorr[i][j];
37 
38  if(fFluxNumu){
39  delete fFluxNumu;
40  delete fFluxNue;
41  delete fFluxNumubar;
42  delete fFluxNuebar;
43  }
44  }
45 
46  //----------------------------------------------------------------------
48  {
49  if(!kIsGibuu_NT(nu)) return 1;
50 
51  const int tgtBin = TargetBin(nu->tgtA);
52  const int chanBin = ChannelBin(nu->iscc, nu->pdg);
53 
54  TH1*& h = fCorr[tgtBin][chanBin];
55  if(!h){
56  const std::string path = getenv("GIBUU_LIBS_LIB_PATH");
57  const std::string fname = path+"/"+ChannelString(chanBin)+"/"+TargetString(tgtBin)+"/corr.root";
58  TFile f(fname.c_str());
59  if(f.IsZombie()){
60  std::cout << "GiBUUWeights.cxx: Couldn't find " << fname << std::endl;
61  abort();
62  }
63  h = (TH1*)f.Get("corr");
64  h->SetDirectory(0);
65  }
66 
67  if(!fFluxNumu) LoadFluxes();
68 
69  // The weight that was applied in the file contains a part intrinsic to
70  // gibuu, which we want, and another part correcting for the genie/gibuu
71  // total cross section, which was applied assuming wrong fluxes and we need
72  // to fix up.
73 
74  // It assumed we always used numu FHC flux, so we'll cancel that
75  // contribution
76  const double fold = fFluxNumu->GetBinContent(h->FindBin(nu->E));
77 
78  // Look up the flux we should have actually used.
79  TH1* hnew = 0;
80  if(nu->iscc){
81  if(nu->pdg == +12) hnew = fFluxNue;
82  if(nu->pdg == -12) hnew = fFluxNuebar;
83  if(nu->pdg == +14) hnew = fFluxNumu;
84  if(nu->pdg == -14) hnew = fFluxNumubar;
85  }
86  else{
87  // NCs were simulated with numu(bar) flux
88  if(nu->pdg > 0) hnew = fFluxNumu;
89  if(nu->pdg < 0) hnew = fFluxNumubar;
90  }
91 
92  // Protection against unknown events
93  if(!hnew) return nu->genweight;
94 
95  const double fnew = hnew->GetBinContent(h->FindBin(nu->E));
96 
97  // Avoid divide-by-zero
98  if(fnew == 0) return nu->genweight;
99 
100  // Apply the final correction
101  return nu->genweight * fold / fnew;
102  }
103 
104  //----------------------------------------------------------------------
106  {
107  switch(A){
108  case 1: return 0;
109  case 12: return 1;
110  case 16: return 2;
111  case 35: return 3;
112  case 48: return 4;
113  case 56: return 5;
114  default:
115  std::cout << "GiBUUWeights.cxx: Unknown nuclear target A = "
116  << A << std::endl;
117  abort();
118  }
119  }
120 
121  //----------------------------------------------------------------------
123  {
124  switch(bin){
125  case 0: return "H1";
126  case 1: return "C12";
127  case 2: return "O16";
128  case 3: return "Cl35";
129  case 4: return "Ti48";
130  case 5: return "Fe56";
131  default: abort();
132  }
133  }
134 
135  //----------------------------------------------------------------------
136  int FixGibuuWeight::ChannelBin(bool iscc, int pdg) const
137  {
138  if(iscc){
139  switch(pdg){
140  case +12: return 0;
141  case -12: return 1;
142  case +14: return 2;
143  case -14: return 3;
144  default:
145  std::cout << "GiBUUWeights.cxx: Unknown neutrino flavour: "
146  << pdg << std::endl;
147  abort();
148  }
149  }
150  else{
151  if(pdg > 0) return 4; else return 5;
152  }
153  }
154 
155  //----------------------------------------------------------------------
157  {
158  switch(bin){
159  case 0: return "cc_nue";
160  case 1: return "ccbar_nue";
161  case 2: return "cc_numu";
162  case 3: return "ccbar_numu";
163  case 4: return "nc_numu";
164  case 5: return "ncbar_numu";
165  default: abort();
166  }
167  }
168 
169  //----------------------------------------------------------------------
171  {
172  // For the TTree::Draw() commands to work we apparently need the histograms
173  // to be in the current directory.
174  const bool backup = TH1::AddDirectoryStatus();
175  TH1::AddDirectory(true);
176 
177  const double Emax = 30;
178  const double dE = .05;
179 
180  std::vector<double> edges;
181  for(double x = 0; x < 2.999; x += dE) edges.push_back(x);
182  for(double x = 3; x < Emax-.0001; x += 1) edges.push_back(x);
183  edges.push_back(Emax);
184 
185  const std::string path = FindCAFAnaDir()+"/data/gibuu/2016/";
186 
187  TTree trFluxNumu, trFluxNue, trFluxNumubar, trFluxNuebar;
188  trFluxNumu .ReadFile((path+"NOvA-ND-FHC-numu.dat" ).c_str(), "Enu/F:flux/F");
189  trFluxNue .ReadFile((path+"NOvA-ND-FHC-nue.dat" ).c_str(), "Enu/F:flux/F");
190  trFluxNumubar.ReadFile((path+"NOvA-ND-RHC-anumu.dat").c_str(), "Enu/F:flux/F");
191  trFluxNuebar .ReadFile((path+"NOvA-ND-RHC-anue.dat" ).c_str(), "Enu/F:flux/F");
192 
193  fFluxNumu = new TH1F(UniqueName().c_str(), "", edges.size()-1, &edges[0]);
194  fFluxNue = new TH1F(UniqueName().c_str(), "", edges.size()-1, &edges[0]);
195  fFluxNumubar = new TH1F(UniqueName().c_str(), "", edges.size()-1, &edges[0]);
196  fFluxNuebar = new TH1F(UniqueName().c_str(), "", edges.size()-1, &edges[0]);
197 
198  const std::string E = "Enu>>";
199  trFluxNumu .Draw((E+fFluxNumu ->GetName()).c_str(), "flux", "goff norm");
200  trFluxNue .Draw((E+fFluxNue ->GetName()).c_str(), "flux", "goff norm");
201  trFluxNumubar.Draw((E+fFluxNumubar->GetName()).c_str(), "flux", "goff norm");
202  trFluxNuebar .Draw((E+fFluxNuebar ->GetName()).c_str(), "flux", "goff norm");
203 
204  fFluxNumu ->SetDirectory(0);
205  fFluxNue ->SetDirectory(0);
206  fFluxNumubar->SetDirectory(0);
207  fFluxNuebar ->SetDirectory(0);
208 
209  TH1::AddDirectory(backup);
210  }
211 }
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
TH1 * fCorr[6][6]
Definition: GiBUUWeights.h:38
Proxy for caf::SRNeutrino.
Definition: SRProxy.h:510
caf::Proxy< float > genweight
Definition: SRProxy.h:531
caf::Proxy< int > tgtA
Definition: SRProxy.h:562
caf::Proxy< int > mode
Definition: SRProxy.h:543
caf::Proxy< float > E
Definition: SRProxy.h:520
std::string FindCAFAnaDir()
Definition: Utilities.cxx:204
double dE
_Cut< caf::SRNeutrinoProxy > NuTruthCut
Cut designed to be used over the nuTree, ie all neutrinos, not just those that got slices...
Definition: Cut.h:104
void LoadFluxes() const
std::string TargetString(int bin) const
double operator()(const caf::SRNeutrinoProxy *nu) const
int iscc
std::string GetName(int i)
Float_t E
Definition: plot.C:20
const NuTruthVar kRawGibuuWeight_NT([](const caf::SRNeutrinoProxy *nu){if(!kIsGibuu_NT(nu)) return 1.f;return float(nu->genweight);})
Definition: GiBUUWeights.h:15
int TargetBin(int z) const
std::string getenv(std::string const &name)
caf::Proxy< bool > iscc
Definition: SRProxy.h:538
_Var< caf::SRNeutrinoProxy > NuTruthVar
Var designed to be used over the nuTree, ie all neutrinos, not just those that got slices...
Definition: Var.h:82
std::string ChannelString(int bin) const
const double j
Definition: BetheBloch.cxx:29
float bin[41]
Definition: plottest35.C:14
OStream cout
Definition: OStream.cxx:6
const std::string path
Definition: plot_BEN.C:43
static const double A
Definition: Units.h:82
int ChannelBin(bool iscc, int pdg) const
const NuTruthCut kIsGibuu_NT([](const caf::SRNeutrinoProxy *nu){ return nu->mode< 0;})
Definition: GiBUUWeights.h:10
double Emax
caf::Proxy< short int > pdg
Definition: SRProxy.h:552
std::string UniqueName()
Return a different string each time, for creating histograms.
Definition: Utilities.cxx:29
enum BeamMode string