Functions | Variables
plot_xsecs.C File Reference
#include "TCanvas.h"
#include "TFile.h"
#include "TGraph.h"
#include "TH2.h"
#include "TLegend.h"
#include "TPad.h"
#include "TTree.h"
#include "TVectorD.h"
#include <iostream>
#include <map>
#include <string>

Go to the source code of this file.

Functions

void plot_xsecs ()
 

Variables

const int cols [] = {kRed, kBlue, kGreen+2, kMagenta, kViolet, kGray+1}
 
const std::map< std::string, intnnucs
 

Function Documentation

void plot_xsecs ( )

Definition at line 23 of file plot_xsecs.C.

References ana::assert(), file_size_ana::axes, cols, om::cout, cut, dE, Emax, allTimeWatchdog::endl, exit(), fin, flux, genie::utils::style::Format(), submit_syst::fout, MECModelEnuComparisons::g, cet::getenv(), MECModelEnuComparisons::i, kBlue, MECModelEnuComparisons::leg, submit_nova_art::nevts, ratio(), string, make_root_from_grid_output::tr, and submit_syst::x.

24 {
25  const double Emax = 30;
26  const double dE = .05; // Must match (or be a multiple of) flux .dat binning (.05)
27 
28  std::vector<double> edges;
29  for(double x = 0; x < 0.999; x += dE*2 ) edges.push_back(x);
30  for(double x = 1; x < 1.999; x += dE*4 ) edges.push_back(x);
31  for(double x = 2; x < 2.999; x += dE*10) edges.push_back(x);
32  for(double x = 3; x < 4.999; x += dE*20) edges.push_back(x);
33  for(double x = 5; x < Emax-.0001; x += dE*50) edges.push_back(x);
34  edges.push_back(Emax);
35 
36  if(!getenv("GIBUU_DIR")){
37  std::cout << "GIBUU_DIR not set" << std::endl;
38  exit(1);
39  }
40  std::string fluxpath = getenv("GIBUU_DIR");
41  fluxpath += "/Linux64bit+2.6-2.12-e15/GiBUU/buuinput2019/neutrino/";
42 
43  TFile fin("/pnfs/nova/persistent/users/bckhouse/gibuuLib2020_sigmacut/hadd.lib.root");
44 
45  TFile fout("xsecs.root", "RECREATE");
46 
47  for(std::string probe_base: {"nu_e", "nu_mu"}){
48  for(bool anti: {false, true}){
49  std::string fluxname = "NOvA-ND-";
50  if(anti) fluxname += "RHC-"; else fluxname += "FHC-";
51  if(anti) fluxname += "a";
52  if(probe_base == "nu_mu") fluxname += "numu"; else fluxname += "nue";
53  const std::string plotname = fluxname + ".pdf";
54  fluxname += ".dat";
55 
56  TTree* trFlux = new TTree;
57  trFlux->ReadFile((fluxpath+"/"+fluxname).c_str(), "Enu/F:flux/F");
58 
59  TH1* flux = new TH1F("flux", (fluxname+";True neutrino energy (GeV);Normalized flux").c_str(), edges.size()-1, &edges[0]);
60  trFlux->Draw("Enu>>flux", "flux", "goff norm");
61  flux->Scale(1./flux->Integral(0, -1));
62 
63  flux->Draw("hist");
64  gPad->Print(plotname.c_str());
65 
66  for(std::string current: {"nc", "cc"}){
67  if(current == "nc" && probe_base != "nu_mu") continue;
68 
69  std::string probe = probe_base;
70  if(current == "nc") probe = "nu";
71  if(anti) probe += "_bar";
72 
73  TLegend* leg = new TLegend(.6, .15, .85, .4);
74  leg->SetFillStyle(0);
75  leg->SetNColumns(2);
76 
77  new TCanvas;
78  int colIdx = 0;
79  TH2* axes = new TH2F("", ("GiBUU "+probe+" "+current+";Neutrino energy (GeV);Total cross section / nucleus (10^{-38} cm^{2})").c_str(), 100, 0, 10, 100, 0, 8);
80  axes->SetDirectory(0);
81  axes->Draw();
82  for(std::string tgt: {"H1", "C12", "O16", "Cl35", "Ti48", "Fe56"}){
83  const std::string dirname = TString::Format("%s/%s/%s",
84  tgt.c_str(),
85  current.c_str(),
86  probe.c_str()).Data();
87 
88  TTree* tr = (TTree*)fin.Get((dirname+"/records").c_str());
89  if(!tr) continue;
90 
91  std::string cut = "";
92  // There is some issue with NCbars where these events are produced
93  // with extremely large weights.
94  if(current == "nc" && anti) cut = "prod_id != 32 && prod_id != 33";
95 
96  TH1* hN = new TH1F("hN", "", edges.size()-1, &edges[0]);
97  hN->SetLineColor(kBlue);
98  tr->Draw("Enu>>hN", cut.c_str(), "goff");
99 
100  hN->Scale(1./10000);// have to account for files that produced no output
101 
102  // Unfortunately we have to hardcode this number - how many events we
103  // requested in the gibuu card file.
104  int nevts = 0;
105  if(tgt == "H1") nevts = 6452;
106  if(tgt == "C12") nevts = 40000;
107  if(tgt == "O16") nevts = 1340;
108  if(tgt == "Cl35") nevts = 3302;
109  if(tgt == "Ti48") nevts = 483;
110  if(tgt == "Fe56") nevts = 3576;
111 
112  hN->Scale(3000./nevts); // rejection sampling cutoff / nevts per file
113 
114  // No, per nucleus now
115  // hN->Scale(1./nnucs.at(tgt)); // per nucleon
116 
117  TH1* ratio = (TH1*)hN->Clone();
118  ratio->SetDirectory(0);
119  ratio->Divide(flux);
120 
121  TGraph* gxsec = new TGraph;
122  gxsec->SetPoint(0, 0, 0);
123  for(int i = 1; i <= ratio->GetNbinsX(); ++i){
124  gxsec->SetPoint(i,
125  ratio->GetXaxis()->GetBinCenter(i),
126  ratio->GetBinContent(i));
127  }
128 
129  gxsec->SetLineWidth(2);
130  gxsec->SetLineColor(cols[colIdx++]);
131 
132  gxsec->Draw("same");
133 
134  leg->AddEntry(gxsec, tgt.c_str(), "l");
135 
136  TDirectory* dtgt = fout.GetDirectory(tgt.c_str());
137  if(!dtgt) dtgt = fout.mkdir(tgt.c_str());
138  TDirectory* dcurrent = dtgt->GetDirectory(current.c_str());
139  if(!dcurrent) dcurrent = dtgt->mkdir(current.c_str());
140  TDirectory* dprobe = dcurrent->GetDirectory(probe.c_str());
141  if(!dprobe) dprobe = dcurrent->mkdir(probe.c_str());
142 
143  dprobe->cd();
144  gxsec->Write("xsec");
145  } // end for tgt
146 
147  leg->Draw();
148  gPad->Print(("gibuu_xsec_"+probe+"_"+current+".pdf").c_str());
149  } // end for current
150  } // end for anti
151  } // end for probe_base
152 
153  // Now genie
154 
155  const char* xsec_dir = getenv("GENIEXSECPATH");
156  const std::string xsec_fname = xsec_dir+std::string("/xsec_graphs.root");
157 
158  TFile f_xsec(xsec_fname.c_str());
159  assert(!f_xsec.IsZombie());
160 
161  for(std::string current: {"nc", "cc"}){
162  for(bool anti: {false, true}){
163  for(std::string probe_base: {"nu_e", "nu_mu"}){
164  std::string probe = probe_base;
165  if(anti) probe += "_bar";
166 
167  new TCanvas;
168  TH2* axes = new TH2F("", ";Neutrino energy (GeV);Total cross section / nucleus (10^{-38} cm^{-2})", 100, 0, 10, 100, 0, 8);
169  axes->SetDirectory(0);
170  axes->SetTitle(("GENIE "+probe+" "+current).c_str());
171  axes->Draw();
172  int colIdx = 0;
173 
174  for(std::string tgt: {"H1", "C12", "O16", "Cl35", "Ti48", "Fe56"}){
175  const std::string genieStr = probe+"_"+tgt+"/tot_"+current;
176  TGraph* g = (TGraph*)f_xsec.Get(genieStr.c_str());
177  if(!g){
178  std::cout << genieStr << " not found" << std::endl;
179  continue;
180  }
181 
182  // No, per nucleus
183  /*
184  const int nnuc = nnucs.at(tgt);
185  for(int i = 0; i < g->GetN(); ++i){
186  double x, y;
187  g->GetPoint(i, x, y);
188  g->SetPoint(i, x, y/nnuc);
189  }
190  */
191 
192  g->Draw("l same");
193  g->SetLineWidth(2);
194  g->SetLineColor(cols[colIdx++]);
195  } // end for nucleus
196  // leg->Draw();
197  gPad->Print(TString::Format("genie_xsec_%s_%s.pdf", probe.c_str(), current.c_str()).Data());
198  } // end for pdg
199  } // end for sign
200  } // end for iscc
201 }
TString fin
Definition: Style.C:24
nevts
Setup runNovaSAM #.
TH1 * ratio(TH1 *h1, TH1 *h2)
double dE
Loaders::FluxType flux
std::string getenv(std::string const &name)
OStream cout
Definition: OStream.cxx:6
const Cut cut
Definition: exporter_fd.C:30
const int cols[]
Definition: plot_xsecs.C:14
exit(0)
assert(nhit_max >=nhit_nbins)
double Emax
void Format(TGraph *gr, int lcol, int lsty, int lwid, int mcol, int msty, double msiz)
Definition: Style.cxx:154
enum BeamMode kBlue
enum BeamMode string

Variable Documentation

const int cols[] = {kRed, kBlue, kGreen+2, kMagenta, kViolet, kGray+1}

Definition at line 14 of file plot_xsecs.C.

Referenced by plot_xsecs().

const std::map<std::string, int> nnucs
Initial value:
= {{"H1", 1},
{"C12", 12},
{"O16", 16},
{"Cl35", 35},
{"Ti48", 48},
{"Fe56", 56}}

Definition at line 16 of file plot_xsecs.C.