Functions | Variables
nuint09_1pi1.C File Reference
#include <iomanip>

Go to the source code of this file.

Functions

void nuint09_1pi1 (int isample, int single_pion_sources=0, int stage=1)
 

Variables

const int kNSamples = 3
 
const int kNWCur = 1
 
const int kNRunsPerCase = 10
 
const int kNEvtPerRun = 100000
 
const char * kLabel [kNSamples]
 
const int kRunNu1PI1 [kNSamples][kNWCur][kNRunsPerCase]
 
int kA [kNSamples]
 

Function Documentation

void nuint09_1pi1 ( int  isample,
int  single_pion_sources = 0,
int  stage = 1 
)

Definition at line 52 of file nuint09_1pi1.C.

53 {
54  cout << " ***** running: 1PI.1" << endl;
55 
56  if(isample<0 || isample >= kNSamples) return;
57  if(single_pion_sources<0 || single_pion_sources>2) return;
58 
59  const char * label = kLabel[isample];
60 
61  int A = kA[isample];
62 
63  // get cross section graphs
64 
65  TFile fsig("../sig/splines.root","read");
66  TDirectory * sig_dir = (TDirectory *) fsig.Get(label);
67 
68  TGraph * sig_graph_totcc = (TGraph*) sig_dir->Get("tot_cc");
69 
70  // range & spacing
71 
72  const int nEnu = 60;
73  const double Enumin = 0.05;
74  const double Enumax = 30.00;
75 
76  // create output stream
77 
78  ostringstream out_filename;
79  out_filename << label;
80 
81  if (single_pion_sources==0) out_filename << ".1pi_1a.";
82  else if (single_pion_sources==1) out_filename << ".1pi_1b.";
83  else if (single_pion_sources==2) out_filename << ".1pi_1c.";
84 
85  if(stage==0) out_filename << "no_FSI.";
86 
87  out_filename << "sig1pi_vs_Enu.data";
88  ofstream out_stream(out_filename.str().c_str(), ios::out);
89 
90  // write out txt file
91 
92  out_stream << "# [" << label << "]" << endl;
93  out_stream << "# " << endl;
94  out_stream << "# [1PI.1]:" << endl;
95  out_stream << "# Total cross section for 1 pion (pi+ only) production as a function of energy" << endl;
96  if(stage==0) {
97  out_stream << "# ***** NO FSI: The {X pi+} state is a primary hadronic state" << endl;
98  }
99  if(single_pion_sources==0) {
100  out_stream << "# 1pi sources: All" << endl;
101  }
102  else if(single_pion_sources==1) {
103  out_stream << "# 1pi sources: P33(1232) resonance only" << endl;
104  }
105  else if(single_pion_sources==2) {
106  out_stream << "# 1pi sources: All resonances only" << endl;
107  }
108  out_stream << "# " << endl;
109  out_stream << "# Note:" << endl;
110  out_stream << "# - neutrino energy E in GeV, linear spacing between Emin = " << Enumin << " GeV, Emax = " << Enumax << " GeV " << endl;
111  out_stream << "# - cross sections in 1E-38 cm^2 " << endl;
112  out_stream << "# - quoted cross section is nuclear cross section divided with number of nucleons A" << endl;
113  out_stream << "# Columns:" << endl;
114  out_stream << "# | Energy | sig(nu_mu + A -> mu- 1pi+ X) | " << endl;
115 
116  out_stream << std::fixed << setprecision(6);
117 
118  //
119  // load event data
120  //
121 
122  TChain * chain = new TChain("gst");
123 
124  // loop over CC/NC cases
125  for(int iwkcur=0; iwkcur<kNWCur; iwkcur++) {
126  // loop over runs for current case
127  for(int ir=0; ir<kNRunsPerCase; ir++) {
128  // build filename
129  ostringstream filename;
130  int run_number = kRunNu1PI1[isample][iwkcur][ir];
131 
132  cout << "isample = " << isample << ", iwkcur = " << iwkcur << ", ir = " << ir << ", run = " << run_number << endl;
133 
134  filename << "../gst/gntp." << run_number << ".gst.root";
135  // add to chain
136  cout << "Adding " << filename.str() << " to event chain" << endl;
137  chain->Add(filename.str().c_str());
138  }
139  }
140 
141  //
142  // book histograms
143  //
144  TH1D * hst_cc1pip = new TH1D("hst_cc1pip","CC1pi+ events, Enu spectrum", nEnu, Enumin, Enumax);
145  TH1D * hst_allcc = new TH1D("hst_allcc", "all CC events, Enu spectrum", nEnu, Enumin, Enumax);
146 
147  //
148  // fill histograms
149  //
150 
151  chain->Draw("Ev>>hst_allcc", "cc","GOFF");
152 
153  if(stage==1) {
154  if(single_pion_sources==0) {
155  //all sources
156  chain->Draw("Ev>>hst_cc1pip","cc&&pdgf==211&&nfpip==1&&nfpim==0&&nfpi0==0","GOFF");
157  }
158  else if(single_pion_sources==1) {
159  //P33(1232) only
160  chain->Draw("Ev>>hst_cc1pip","cc&&resid==0&&res&&pdgf==211&&nfpip==1&&nfpim==0&&nfpi0==0","GOFF");
161  }
162  else if(single_pion_sources==2) {
163  //all resonances only
164  chain->Draw("Ev>>hst_cc1pip","cc&&res&&pdgf==211&&nfpip==1&&nfpim==0&&nfpi0==0","GOFF");
165  }
166  }
167  else if(stage==0) {
168  if(single_pion_sources==0) {
169  //all sources
170  chain->Draw("Ev>>hst_cc1pip","cc&&pdgi==211&&nipip==1&&nipim==0&&nipi0==0","GOFF");
171  }
172  else if(single_pion_sources==1) {
173  //P33(1232) only
174  chain->Draw("Ev>>hst_cc1pip","cc&&resid==0&&res&&pdgi==211&&nipip==1&&nipim==0&&nipi0==0","GOFF");
175  }
176  else if(single_pion_sources==2) {
177  //all resonances only
178  chain->Draw("Ev>>hst_cc1pip","cc&&res&&pdgi==211&&nipip==1&&nipim==0&&nipi0==0","GOFF");
179  }
180  }
181 
182  cout << "CC1pi+ nevents: " << hst_cc1pip->GetEntries() << endl;
183  cout << "ALLCC nevents: " << hst_allcc->GetEntries() << endl;
184 
185  //
186  // compute sig(CC1pi+) and write out in txt file expected by NuINT organizers
187  // Also write out root graph in temp file to be re-used at later benchmarking calc
188  //
189 
190  double e[nEnu], sig[nEnu];
191 
192  for(int i=1; i <= hst_cc1pip->GetNbinsX(); i++) {
193 
194  double Enu = hst_cc1pip->GetBinCenter(i);
195 
196  double Ncc1pip = hst_cc1pip -> GetBinContent(i);
197  double Nallcc = hst_allcc -> GetBinContent(i);
198 
199  double sig_cc1pip=0;
200  if(Nallcc>0) { sig_cc1pip = (Ncc1pip/Nallcc) * sig_graph_totcc->Eval(Enu) / A; }
201 
202  out_stream << setw(15) << Enu << setw(15) << sig_cc1pip << endl;
203 
204  e [i-1] = Enu;
205  sig[i-1] = sig_cc1pip;
206  }
207 
208  out_stream.close();
209 
210  TFile ftmp("./cc1pip_tmp.root","recreate");
211  TGraph * grCC1pip = new TGraph(nEnu,e,sig);
212  grCC1pip->Write("CC1pip");
213  hst_allcc->Write();
214  hst_cc1pip->Write();
215 
216  // visual inspection
217  TCanvas * c1 = new TCanvas("c1","",20,20,500,500);
218  grCC1pip->Draw("alp");
219  c1->Update();
220 
221  ftmp.Close();
222 
223  delete chain;
224 }
i
size_t i(0)
kNRunsPerCase
const int kNRunsPerCase
Definition: nuint09_1pi1.C:21
cout
the ParameterSet object passed in for the configuration of a destination should be the only source that can affect the behavior of that destination This is to eliminate the dependencies of configuring a destination from multiple mostly from the defaults It suppresses possible glitches about changing the configuration file somewhere outside of a destination segament might still affect the behavior of that destination In the previous configuration for a specific the value of a certain e may come from following and have been suppressed It the configuring ParameterSet object for each destination will be required to carry a parameter list as complete as possible If a parameter still cannot be found in the ParameterSet the configuration code will go look for a hardwired default directly The model is a great simplicity comparing with the previous especially when looking for default values Another great advantage is most of the parameters now have very limited places that allows to appear Usually they can only appear at one certain level in a configuration file For in the old configuring model or in a default ParameterSet object inside of a or in a category or in a severity object This layout of multiple sources for a single parameter brings great confusion in both writing a configuration and in processing the configuration file Under the new the only allowed place for the parameter limit to appear is inside of a category which is inside of a destination object Other improvements simplify the meaning of a destination name In the old a destination name has multiple folds of meanings the e cout and cerr have the special meaning of logging messages to standard output or standard error the name also serves as the output filename if the destination is a file these names are also references to look up for detailed configurations in configuring the MessageFacility The multi purpose of the destination name might cause some unwanted behavior in either writing or parsing the configuration file To amend in the new model the destination name is now merely a name for a which might represent the literal purpose of this or just an id All other meanings of the destinations names now go into the destination ParameterSet as individual such as the type parameter and filename parameter Following is the deatiled rule for the new configuring Everything that is related with MessageFacility configuration must be wrapped in a single ParameterSet object with the name MessageFacility The MessageFacility ParameterSet object contains a series of top level parameters These parameters can be chosen a vector of string listing the name of debug enabled models Or use *to enable debug messages in all modules a vector of string a vector of string a vector of string a ParameterSet object containing the list of all destinations The destinations ParameterSet object is a combination of ParameterSet objects for individual destinations There are two types of destinations that you can insert in the destinations ParameterSet ordinary including cout
Definition: MessageFacilityAsStandAlonePackage.txt:82
e
const double e
Definition: gUpMuFluxGen.cxx:165
kNWCur
const int kNWCur
Definition: nuint09_1pi1.C:20
kA
int kA[kNSamples]
Definition: nuint09_1pi1.C:44
kRunNu1PI1
const int kRunNu1PI1[kNSamples][kNWCur][kNRunsPerCase]
Definition: nuint09_1pi1.C:31
genie::units::A
static const double A
Definition: Units.h:82
kNSamples
const int kNSamples
Definition: nuint09_1pi1.C:19
kLabel
const char * kLabel[kNSamples]
Definition: nuint09_1pi1.C:24

Variable Documentation

int kA[kNSamples]
Initial value:
=
{
12,
16,
56
}

Definition at line 44 of file nuint09_1pi1.C.

const char* kLabel[kNSamples]
Initial value:
=
{
"nu_mu_C12",
"nu_mu_O16",
"nu_mu_Fe56"
}

Definition at line 24 of file nuint09_1pi1.C.

const int kNEvtPerRun = 100000

Definition at line 22 of file nuint09_1pi1.C.

const int kNRunsPerCase = 10

Definition at line 21 of file nuint09_1pi1.C.

const int kNSamples = 3

Definition at line 19 of file nuint09_1pi1.C.

const int kNWCur = 1

Definition at line 20 of file nuint09_1pi1.C.

const int kRunNu1PI1[kNSamples][kNWCur][kNRunsPerCase]
Initial value:
=
{
{
{909900, 909901, 909902, 909903, 909904, 909905, 909906, 909907, 909908, 909909}
},
{
{919900, 919901, 919902, 919903, 919904, 919905, 919906, 919907, 919908, 919909}
},
{
{929900, 929901, 929902, 929903, 929904, 929905, 929906, 929907, 929908, 929909}
}
}

Definition at line 31 of file nuint09_1pi1.C.

Generated by   doxygen 1.8.11