test_flux.C
Go to the documentation of this file.
3 #include "CAFAna/Vars/Vars.h"
5 #include "CAFAna/XSec/Flux.h"
6 using namespace ana;
7 
8 #include "TCanvas.h"
9 #include "TFile.h"
10 #include "TH1.h"
11 #include "TPad.h"
12 
13 #include <iostream>
14 
15 //----------------------------------------------------------------------
17 {
18  TFile fin("$NOVA_ANA/nu_int_ana/NumiFlux/nu_flux_v08.root");
19 
20  TH1* h = 0;
21  switch(pdg){
22  case +12: h = (TH1*)fin.Get("nue_flux"); break;
23  case +14: h = (TH1*)fin.Get("numu_flux"); break;
24  case -12: h = (TH1*)fin.Get("nuebar_flux"); break;
25  case -14: h = (TH1*)fin.Get("numubar_flux"); break;
26  default:
27  std::cerr << "Unknown flux flavor " << pdg << std::endl;
28  abort();
29  }
30  assert(h);
31 
32  h->SetTitle("");
33 
34  // According to the title in the file at least this was per 1e10 POT
35  return Spectrum(h, 1e10, 0, "True neutrino energy (GeV)");
36 }
37 
38 //----------------------------------------------------------------------
40 {
41  std::unique_ptr<TH1D> h(NominalFluxFromHist(pdg).ToTH1(1));
42 
43  // This is not the most elegant solution. Probably both sites should call a
44  // third function that just gives the weights. But it works.
47  sr.hdr.det = caf::kNEARDET;
48  sr.mc.nnu = 1;
49  sr.mc.nu.resize(1);
50  sr.mc.nu[0].pdgorig = pdg;
51 
52  for(int i = 1; i <= h->GetNbinsX(); ++i){
53  sr.mc.nu[0].E = h->GetBinCenter(i);
54  h->SetBinContent(i, w(&sr) * h->GetBinContent(i));
55  }
56 
57  return Spectrum(std::move(h), 1, 0, "True neutrino energy (GeV)");
58 }
59 
60 //----------------------------------------------------------------------
62 {
63  std::unique_ptr<TH1D> h(NominalFluxFromHist(pdg).ToTH1(1));
64 
65  // This is not the most elegant solution. Probably both sites should call a
66  // third function that just gives the weights. But it works.
69  sr.hdr.det = caf::kNEARDET;
70  sr.mc.nnu = 1;
71  sr.mc.nu.resize(1);
72  sr.mc.nu[0].pdgorig = pdg;
73 
74  for(int i = 1; i <= h->GetNbinsX(); ++i){
75  sr.mc.nu[0].E = h->GetBinCenter(i);
76  h->SetBinContent(i, w(&sr) * h->GetBinContent(i));
77  }
78 
79  return Spectrum(std::move(h), 1, 0, "True neutrino energy (GeV)");
80 }
81 
82 
83 void test_flux()
84 {
85  // Use FA files because they still have mc.gen in them
86  // Remember you can use --stride to speed things up
87  // SpectrumLoader loader("prod_caf_S15-05-22a_nd_genie_fhc_nonswap_ndnewpos");
88  // At Caltech these files are available here:
89  SpectrumLoader loader("/nfs/raid11/novasam/prod_caf_S15-05-22a_nd_genie_fhc_nonswap_ndnewpos/*.root");
90 
91  TargetCount n;
92  std::cout << n.Mass() << " kg total, "
93  << n.NNucleons() << " nucleons" << std::endl;
94 
95  TargetCount nC(6);
96  std::cout << nC.Mass() << " kg of Carbon, "
97  << nC.NNucleons() << " nucleons" << std::endl;
98 
99  TargetCount nH(1);
100  std::cout << nH.Mass() << " kg of Hydrogen, "
101  << nH.NNucleons() << " nucleons" << std::endl;
102 
103  TargetCount nCl(17);
104  std::cout << nCl.Mass() << " kg of Chlorine, "
105  << nCl.NNucleons() << " nucleons" << std::endl;
106 
107  TFile* fout = new TFile("flux_tests.root", "RECREATE");
108 
109  for(bool hist: {true, false}){
110  new TCanvas;
111 
112  std::map<int, std::map<int, Spectrum*>> fluxes;
113 
114  for(int flav: {-14, -12, +12, +14}){
115  if(hist){
116  fluxes[flav][0] = new Spectrum(NominalFluxFromHist(flav));
117  fluxes[flav][1] = 0;//new Spectrum(MinervaFluxFromHist(flav));
118  fluxes[flav][2] = 0;//new Spectrum(MippNA49FluxFromHist(flav));
119  }
120  else{
121  fluxes[flav][0] = /*NominalFlux*/DeriveFlux(loader, Binning::Simple(200, 0, 20), flav);
122  fluxes[flav][1] = 0;//MinervaFlux(loader, Binning::Simple(200, 0, 20), flav);
123  fluxes[flav][2] = 0;//MippNA49Flux(loader, Binning::Simple(200, 0, 20), flav);
124  }
125  } // end for flav
126 
127  if(!hist) loader.Go();
128 
129  TH1* h = fluxes[14][0]->ToTH1(1, kBlue);
130  h->GetYaxis()->SetTitle("Neutrinos / m^{2} / POT / bin");
131  h->GetYaxis()->SetRangeUser(0, 1e-5);
132  h->GetXaxis()->SetRangeUser(0, 5);
133  h->Draw("hist");
134 
135  for(int flav: {-14, -12, +12, +14}){
136  for(int syst = 0; syst < 3; ++syst){
137  int col;
138  int style = (flav > 0) ? kSolid : 7;
139 
140  if(abs(flav) == 14){
141  if(syst == 0) col = kBlue;
142  if(syst == 1) col = kCyan;
143  if(syst == 2) col = kBlue+2;
144  }
145  if(abs(flav) == 12){
146  if(syst == 0) col = kRed;
147  if(syst == 1) col = kMagenta;
148  if(syst == 2) col = kRed+2;
149  }
150 
151  if(!fluxes[flav][syst]) continue;
152 
153  fluxes[flav][syst]->ToTH1(1, col, style)->Draw("hist same");
154 
155  fluxes[flav][syst]->SaveTo(fout, TString::Format("flux_%s_%d_%d",
156  hist ? "hist" : "derive",
157  flav,
158  syst).Data());
159  } // end for syst
160  } // end for flav
161 
162  gPad->Update();
163 
164  gPad->Print(hist ? "flux_hist.eps" : "flux_derive.eps");
165  } // end for hist
166 }
Near Detector underground.
Definition: SREnums.h:10
TString fin
Definition: Style.C:24
enum BeamMode kRed
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
caf::Proxy< caf::SRHeader > hdr
Definition: SRProxy.h:2136
Proxy for caf::StandardRecord.
Definition: SRProxy.h:2125
caf::Proxy< std::vector< caf::SRNeutrino > > nu
Definition: SRProxy.h:617
OStream cerr
Definition: OStream.cxx:7
void abs(TH1 *hist)
caf::Proxy< short int > nnu
Definition: SRProxy.h:616
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
double NNucleons() const
Number of nucleons (mass * avogadro&#39;s number)
Definition: TargetCount.h:31
void test_flux()
Definition: test_flux.C:83
Int_t col[ntarg]
Definition: Style.C:29
caf::StandardRecord * sr
Spectrum NominalFluxFromHist(int pdg)
Definition: test_flux.C:16
virtual void Go() override
Load all the registered spectra.
loader
Definition: demo0.py:10
std::vector< float > Spectrum
Definition: Constants.h:610
OStream cout
Definition: OStream.cxx:6
caf::Proxy< caf::SRTruthBranch > mc
Definition: SRProxy.h:2137
Spectrum * DeriveFlux(SpectrumLoaderBase &loader, const Binning &bins, int pdg, const TVector3 *min, const TVector3 *max, const SystShifts &shift, const NuTruthVar weight)
Definition: Flux.cxx:56
Collaborates with Spectrum and OscillatableSpectrum to fill spectra from CAF files.
Spectrum MinervaFluxFromHist(int pdg)
Definition: test_flux.C:39
assert(nhit_max >=nhit_nbins)
void Format(TGraph *gr, int lcol, int lsty, int lwid, int mcol, int msty, double msiz)
Definition: Style.cxx:154
Float_t e
Definition: plot.C:35
enum BeamMode kBlue
static Binning Simple(int n, double lo, double hi, const std::vector< std::string > &labels={})
Definition: Binning.cxx:107
Float_t w
Definition: plot.C:20
Count number of targets within the main part of the ND (vSuperBlockND)
Definition: TargetCount.h:10
double Mass() const
Mass in kg.
Definition: TargetCount.h:29
Spectrum MippNA49FluxFromHist(int pdg)
Definition: test_flux.C:61
caf::Proxy< caf::Det_t > det
Definition: SRProxy.h:231