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

Go to the source code of this file.

Functions

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

Variables

const int kNSamples = 3
 
const int kNWCur = 1
 
const int kNEnergies = 2
 
const int kNRunsPerCase = 5
 
const int kNEvtPerRun = 100000
 
const char * kLabel [kNSamples]
 
const int kRunNu1PI3 [kNSamples][kNWCur][kNEnergies][kNRunsPerCase]
 

Function Documentation

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

Definition at line 49 of file nuint09_1pi3.C.

50 {
51  cout << " ***** running: 1PI.3" << endl;
52 
53  if(isample<0 || isample >= kNSamples) return;
54  if(single_pion_sources<0 || single_pion_sources>2) return;
55 
56  const char * label = kLabel[isample];
57 
58  // get cross section graphs
59 
60  TFile fsig("./cc1pip_tmp.root","read"); // generated at [1PI.1]
61  TGraph * sig_graph_cc1pip = (TGraph*) fsig.Get("CC1pip");
62 
63  // range & spacing
64 
65  const int nKEpip = 60;
66  const double KEpipmin = 0.00;
67  const double KEpipmax = 1.50;
68 
69  const int ncosth = 30;
70  const double costhmin = -1;
71  const double costhmax = +1;
72 
73  // create output stream
74 
75  ostringstream out_filename;
76  out_filename << label;
77 
78  if (single_pion_sources==0) out_filename << ".1pi_3a.";
79  else if (single_pion_sources==1) out_filename << ".1pi_3b.";
80  else if (single_pion_sources==2) out_filename << ".1pi_3c.";
81 
82  if(stage==0) out_filename << "no_FSI.";
83 
84  out_filename << label << "d2sig1pi_dKEpidOmega.data";
85 
86  ofstream out_stream(out_filename.str().c_str(), ios::out);
87 
88  // write out txt file
89 
90  out_stream << "# [" << label << "]" << endl;
91  out_stream << "# " << endl;
92  out_stream << "# [1PI.3]:" << endl;
93  out_stream << "# d^Sigma / dOmega_pi+ dKE_pi+ at E_nu= 1.0 and 1.5 GeV" << endl;
94  if(stage==0) {
95  out_stream << "# ***** NO FSI: The {X pi+} state is a primary hadronic state" << endl;
96  }
97  if(single_pion_sources==0) {
98  out_stream << "# 1pi sources: All" << endl;
99  }
100  else if(single_pion_sources==1) {
101  out_stream << "# 1pi sources: P33(1232) resonance only" << endl;
102  }
103  else if(single_pion_sources==2) {
104  out_stream << "# 1pi sources: All resonances only" << endl;
105  }
106  out_stream << "# " << endl;
107  out_stream << "# Note:" << endl;
108  out_stream << "# - pi+ kinetic energy KE in GeV, linear spacing between KEmin = " << KEpipmin << " GeV, KEmax = " << KEpipmax << " GeV " << endl;
109  out_stream << "# - cross sections in 1E-38 cm^2 /GeV /sr" << endl;
110  out_stream << "# - quoted cross section is nuclear cross section divided with number of nucleons A" << endl;
111  out_stream << "# Columns:" << endl;
112  out_stream << "# | KE(pi+) | cos(theta_pi+) | dsig(numu A -> mu- 1pi+ X; Enu = 1.0 GeV) | dsig(numu A -> mu- 1pi+ X; Enu = 1.5 GeV) | " << endl;
113 
114  out_stream << std::fixed << setprecision(6);
115 
116  //
117  // load event data
118  //
119 
120  TChain * chain = new TChain("gst");
121 
122  // loop over CC/NC cases
123  for(int iwkcur=0; iwkcur<kNWCur; iwkcur++) {
124  // loop over energies
125  for(int ie=0; ie<kNEnergies; ie++) {
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 = kRunNu1PI3[isample][iwkcur][ie][ir];
131  filename << "../gst/gntp." << run_number << ".gst.root";
132  // add to chain
133  cout << "Adding " << filename.str() << " to event chain" << endl;
134  chain->Add(filename.str().c_str());
135  }
136  }
137  }
138 
139  //
140  // get CC1pi+ cross sections at given energies for normalization purposes
141  //
142  double sig_cc1pip_1000MeV = sig_graph_cc1pip -> Eval (1.0);
143  double sig_cc1pip_1500MeV = sig_graph_cc1pip -> Eval (1.5);
144 
145  //
146  // book histograms
147  //
148  TH2D * hst_d2sig_dKEpipdOmg_1000MeV = new TH2D("hst_d2sig_dKEpipdOmg_1000MeV",
149  "d2sig / dKEpi+ dOmega, numu A -> mu- 1pi+ X, Enu=1.0 GeV", nKEpip, KEpipmin, KEpipmax, ncosth, costhmin, costhmax);
150  TH2D * hst_d2sig_dKEpipdOmg_1500MeV = new TH2D("hst_d2sig_dKEpipdOmg_1500MeV",
151  "d2sig / dKEpi+ dOmega, numu A -> mu- 1pi+ X, Enu=1.5 GeV", nKEpip, KEpipmin, KEpipmax, ncosth, costhmin, costhmax);
152 
153  //
154  // fill histograms
155  //
156  if(stage==1) {
157  if(single_pion_sources==0) {
158  // all sources
159  chain->Draw("pzf/sqrt(pzf*pzf+pyf*pyf+pxf*pxf):(Ef-0.139)>>hst_d2sig_dKEpipdOmg_1000MeV","cc&&Ev>0.99&&Ev<1.01&&pdgf==211&&nfpip==1&&nfpim==0&&nfpi0==0","GOFF");
160  chain->Draw("pzf/sqrt(pzf*pzf+pyf*pyf+pxf*pxf):(Ef-0.139)>>hst_d2sig_dKEpipdOmg_1500MeV","cc&&Ev>1.49&&Ev<1.51&&pdgf==211&&nfpip==1&&nfpim==0&&nfpi0==0","GOFF");
161  }
162  else if(single_pion_sources==1) {
163  // P33(1232) only
164  chain->Draw("pzf/sqrt(pzf*pzf+pyf*pyf+pxf*pxf):(Ef-0.139)>>hst_d2sig_dKEpipdOmg_1000MeV","cc&&resid==0&&res&&Ev>0.99&&Ev<1.01&&pdgf==211&&nfpip==1&&nfpim==0&&nfpi0==0","GOFF");
165  chain->Draw("pzf/sqrt(pzf*pzf+pyf*pyf+pxf*pxf):(Ef-0.139)>>hst_d2sig_dKEpipdOmg_1500MeV","cc&&resid==0&&res&&Ev>1.49&&Ev<1.51&&pdgf==211&&nfpip==1&&nfpim==0&&nfpi0==0","GOFF");
166  }
167  else if(single_pion_sources==2) {
168  // all resonances only
169  chain->Draw("pzf/sqrt(pzf*pzf+pyf*pyf+pxf*pxf):(Ef-0.139)>>hst_d2sig_dKEpipdOmg_1000MeV","cc&&res&&Ev>0.99&&Ev<1.01&&pdgf==211&&nfpip==1&&nfpim==0&&nfpi0==0","GOFF");
170  chain->Draw("pzf/sqrt(pzf*pzf+pyf*pyf+pxf*pxf):(Ef-0.139)>>hst_d2sig_dKEpipdOmg_1500MeV","cc&&res&&Ev>1.49&&Ev<1.51&&pdgf==211&&nfpip==1&&nfpim==0&&nfpi0==0","GOFF");
171  }
172  }
173  else if(stage==0) {
174  if(single_pion_sources==0) {
175  // all sources
176  chain->Draw("pzi/sqrt(pzi*pzi+pyi*pyi+pxi*pxi):(Ei-0.139)>>hst_d2sig_dKEpipdOmg_1000MeV","cc&&Ev>0.99&&Ev<1.01&&pdgi==211&&nipip==1&&nipim==0&&nipi0==0","GOFF");
177  chain->Draw("pzi/sqrt(pzi*pzi+pyi*pyi+pxi*pxi):(Ei-0.139)>>hst_d2sig_dKEpipdOmg_1500MeV","cc&&Ev>1.49&&Ev<1.51&&pdgi==211&&nipip==1&&nipim==0&&nipi0==0","GOFF");
178  }
179  else if(single_pion_sources==1) {
180  // P33(1232) only
181  chain->Draw("pzi/sqrt(pzi*pzi+pyi*pyi+pxi*pxi):(Ei-0.139)>>hst_d2sig_dKEpipdOmg_1000MeV","cc&&resid==0&&res&&Ev>0.99&&Ev<1.01&&pdgi==211&&nipip==1&&nipim==0&&nipi0==0","GOFF");
182  chain->Draw("pzi/sqrt(pzi*pzi+pyi*pyi+pxi*pxi):(Ei-0.139)>>hst_d2sig_dKEpipdOmg_1500MeV","cc&&resid==0&&res&&Ev>1.49&&Ev<1.51&&pdgi==211&&nipip==1&&nipim==0&&nipi0==0","GOFF");
183  }
184  else if(single_pion_sources==2) {
185  // all resonances only
186  chain->Draw("pzi/sqrt(pzi*pzi+pyi*pyi+pxi*pxi):(Ei-0.139)>>hst_d2sig_dKEpipdOmg_1000MeV","cc&&res&&Ev>0.99&&Ev<1.01&&pdgi==211&&nipip==1&&nipim==0&&nipi0==0","GOFF");
187  chain->Draw("pzi/sqrt(pzi*pzi+pyi*pyi+pxi*pxi):(Ei-0.139)>>hst_d2sig_dKEpipdOmg_1500MeV","cc&&res&&Ev>1.49&&Ev<1.51&&pdgi==211&&nipip==1&&nipim==0&&nipi0==0","GOFF");
188  }
189  }
190 
191  //
192  // normalize
193  //
194  double norm_cc1pip_1000MeV = hst_d2sig_dKEpipdOmg_1000MeV -> Integral("width") * 2*TMath::Pi() / sig_cc1pip_1000MeV;
195  double norm_cc1pip_1500MeV = hst_d2sig_dKEpipdOmg_1500MeV -> Integral("width") * 2*TMath::Pi() / sig_cc1pip_1500MeV;
196 
197  if (norm_cc1pip_1000MeV > 0) hst_d2sig_dKEpipdOmg_1000MeV -> Scale(1./norm_cc1pip_1000MeV);
198  if (norm_cc1pip_1500MeV > 0) hst_d2sig_dKEpipdOmg_1500MeV -> Scale(1./norm_cc1pip_1500MeV);
199 
200  for(int i = 1; i <= hst_d2sig_dKEpipdOmg_1500MeV->GetNbinsX(); i++) {
201  for(int j = 1; j <= hst_d2sig_dKEpipdOmg_1500MeV->GetNbinsY(); j++) {
202 
203  double KEpip = hst_d2sig_dKEpipdOmg_1500MeV -> GetXaxis() -> GetBinCenter(i);
204  double costhpip = hst_d2sig_dKEpipdOmg_1500MeV -> GetYaxis() -> GetBinCenter(j);
205 
206  double d2sig_dKEpipdOmg_1000MeV = hst_d2sig_dKEpipdOmg_1000MeV -> GetBinContent(i,j);
207  double d2sig_dKEpipdOmg_1500MeV = hst_d2sig_dKEpipdOmg_1500MeV -> GetBinContent(i,j);
208 
209  d2sig_dKEpipdOmg_1000MeV = TMath::Max(0., d2sig_dKEpipdOmg_1000MeV);
210  d2sig_dKEpipdOmg_1500MeV = TMath::Max(0., d2sig_dKEpipdOmg_1500MeV);
211 
212  out_stream << setw(15) << KEpip
213  << setw(15) << costhpip
214  << setw(15) << d2sig_dKEpipdOmg_1000MeV
215  << setw(15) << d2sig_dKEpipdOmg_1500MeV
216  << endl;
217  }
218  }
219 
220  out_stream.close();
221 
222  // visual inspection
223  //TCanvas * c1 = new TCanvas("c1","",20,20,500,500);
224  //...
225  //c1->Update();
226 
227  fsig.Close();
228  delete chain;
229 }
i
size_t i(0)
kNRunsPerCase
const int kNRunsPerCase
Definition: nuint09_1pi3.C:22
kRunNu1PI3
const int kRunNu1PI3[kNSamples][kNWCur][kNEnergies][kNRunsPerCase]
Definition: nuint09_1pi3.C:32
GetYaxis
a1 GetYaxis() -> SetTitle("F2")
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
kNWCur
const int kNWCur
Definition: nuint09_1pi3.C:20
ie
var ie
Definition: slidy.js:14
kNEnergies
const int kNEnergies
Definition: nuint09_1pi3.C:21
kNSamples
const int kNSamples
Definition: nuint09_1pi3.C:19
GetXaxis
a1 GetXaxis() -> SetTitle("Q2")
Definition: plot_lowE_sig_spread_genie.C:84
kLabel
const char * kLabel[kNSamples]
Definition: nuint09_1pi3.C:25

Variable Documentation

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

Definition at line 25 of file nuint09_1pi3.C.

const int kNEnergies = 2

Definition at line 21 of file nuint09_1pi3.C.

const int kNEvtPerRun = 100000

Definition at line 23 of file nuint09_1pi3.C.

const int kNRunsPerCase = 5

Definition at line 22 of file nuint09_1pi3.C.

const int kNSamples = 3

Definition at line 19 of file nuint09_1pi3.C.

const int kNWCur = 1

Definition at line 20 of file nuint09_1pi3.C.

const int kRunNu1PI3[kNSamples][kNWCur][kNEnergies][kNRunsPerCase]
Initial value:
=
{
{
{ {900200, 900201, 900202, 900203, 900204},
{900300, 900301, 900302, 900303, 900304} }
},
{
{ {910200, 910201, 910202, 910203, 910204},
{910300, 910301, 910302, 910303, 910304} }
},
{
{ {920200, 920201, 920202, 920203, 920204},
{920300, 920301, 920302, 920303, 920304} }
}
}

Definition at line 32 of file nuint09_1pi3.C.

Generated by   doxygen 1.8.11