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

Go to the source code of this file.

Functions

void nuint09_1pi4 (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 kRunNu1PI4 [kNSamples][kNWCur][kNEnergies][kNRunsPerCase]
 

Function Documentation

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

Definition at line 49 of file nuint09_1pi4.C.

50 {
51  cout << " ***** running: 1PI.4" << 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 graph
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 nW = 60;
66  const double Wmin = 0.80;
67  const double Wmax = 2.50;
68 
69  // create output stream
70 
71  ostringstream out_filename;
72  out_filename << label;
73 
74  if (single_pion_sources==0) out_filename << ".1pi_4a.";
75  else if (single_pion_sources==1) out_filename << ".1pi_4b.";
76  else if (single_pion_sources==2) out_filename << ".1pi_4c.";
77 
78  if(stage==0) out_filename << "no_FSI.";
79 
80  out_filename << label << "dsig1pi_dW.data";
81 
82  ofstream out_stream(out_filename.str().c_str(), ios::out);
83 
84  // write out txt file
85 
86  out_stream << "# [" << label << "]" << endl;
87  out_stream << "# " << endl;
88  out_stream << "# [1PI.4]:" << endl;
89  out_stream << "# 1pi+ cross section at E_nu= 1.0 and 1.5 GeV as a function of the hadronic invariant mass W" << endl;
90  if(stage==0) {
91  out_stream << "# ***** NO FSI: The {X pi+} state is a primary hadronic state" << endl;
92  }
93  if(single_pion_sources==0) {
94  out_stream << "# 1pi sources: All" << endl;
95  }
96  else if(single_pion_sources==1) {
97  out_stream << "# 1pi sources: P33(1232) resonance only" << endl;
98  }
99  else if(single_pion_sources==2) {
100  out_stream << "# 1pi sources: All resonances only" << endl;
101  }
102  out_stream << "# " << endl;
103  out_stream << "# Note:" << endl;
104  out_stream << "# - W in GeV, linear spacing between Wmin = " << Wmin << " GeV, Wmax = " << Wmax << " GeV " << endl;
105  out_stream << "# - cross sections in 1E-38 cm^2 / GeV " << endl;
106  out_stream << "# - quoted cross section is nuclear cross section divided with number of nucleons A" << endl;
107  out_stream << "# Columns:" << endl;
108  out_stream << "# | W | dsig(numu A -> mu- 1pi+ X; Enu = 1.0 GeV) | dsig(numu A -> mu- 1pi+ X; Enu = 1.5 GeV) | " << endl;
109 
110  out_stream << std::fixed << setprecision(6);
111 
112  //
113  // load event data
114  //
115 
116  TChain * chain = new TChain("gst");
117 
118  // loop over CC/NC cases
119  for(int iwkcur=0; iwkcur<kNWCur; iwkcur++) {
120  // loop over energies
121  for(int ie=0; ie<kNEnergies; ie++) {
122  // loop over runs for current case
123  for(int ir=0; ir<kNRunsPerCase; ir++) {
124  // build filename
125  ostringstream filename;
126  int run_number = kRunNu1PI4[isample][iwkcur][ie][ir];
127  filename << "../gst/gntp." << run_number << ".gst.root";
128  // add to chain
129  cout << "Adding " << filename.str() << " to event chain" << endl;
130  chain->Add(filename.str().c_str());
131  }
132  }
133  }
134 
135  //
136  // get CC1pi+ cross sections at given energies for normalization purposes
137  //
138  double sig_cc1pip_1000MeV = sig_graph_cc1pip -> Eval (1.0);
139  double sig_cc1pip_1500MeV = sig_graph_cc1pip -> Eval (1.5);
140 
141  //
142  // book histograms
143  //
144  TH1D * hst_dsig_dW_1000MeV = new TH1D("hst_dsig_dW_1000MeV","dsig/dW, numu A -> mu- 1pi+ X, Enu=1.0 GeV", nW, Wmin, Wmax);
145  TH1D * hst_dsig_dW_1500MeV = new TH1D("hst_dsig_dW_1500MeV","dsig/dW, numu A -> mu- 1pi+ X, Enu=1.5 GeV", nW, Wmin, Wmax);
146 
147  //
148  // fill histograms
149  //
150  if(stage==1) {
151  if(single_pion_sources==0) {
152  // all sources
153  chain->Draw("Ws>>hst_dsig_dW_1000MeV","cc&&Ev>0.99&&Ev<1.01&&pdgf==211&&nfpip==1&&nfpim==0&&nfpi0==0","GOFF");
154  chain->Draw("Ws>>hst_dsig_dW_1500MeV","cc&&Ev>1.49&&Ev<1.51&&pdgf==211&&nfpip==1&&nfpim==0&&nfpi0==0","GOFF");
155  }
156  else if(single_pion_sources==1) {
157  // P33(1232) only
158  chain->Draw("Ws>>hst_dsig_dW_1000MeV","cc&&resid==0&&res&&Ev>0.99&&Ev<1.01&&pdgf==211&&nfpip==1&&nfpim==0&&nfpi0==0","GOFF");
159  chain->Draw("Ws>>hst_dsig_dW_1500MeV","cc&&resid==0&&res&&Ev>1.49&&Ev<1.51&&pdgf==211&&nfpip==1&&nfpim==0&&nfpi0==0","GOFF");
160  }
161  else if(single_pion_sources==2) {
162  // all resonances only
163  chain->Draw("Ws>>hst_dsig_dW_1000MeV","cc&&res&&Ev>0.99&&Ev<1.01&&pdgf==211&&nfpip==1&&nfpim==0&&nfpi0==0","GOFF");
164  chain->Draw("Ws>>hst_dsig_dW_1500MeV","cc&&res&&Ev>1.49&&Ev<1.51&&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("Ws>>hst_dsig_dW_1000MeV","cc&&Ev>0.99&&Ev<1.01&&pdgi==211&&nipip==1&&nipim==0&&nipi0==0","GOFF");
171  chain->Draw("Ws>>hst_dsig_dW_1500MeV","cc&&Ev>1.49&&Ev<1.51&&pdgi==211&&nipip==1&&nipim==0&&nipi0==0","GOFF");
172  }
173  else if(single_pion_sources==1) {
174  // P33(1232) only
175  chain->Draw("Ws>>hst_dsig_dW_1000MeV","cc&&resid==0&&res&&Ev>0.99&&Ev<1.01&&pdgi==211&&nipip==1&&nipim==0&&nipi0==0","GOFF");
176  chain->Draw("Ws>>hst_dsig_dW_1500MeV","cc&&resid==0&&res&&Ev>1.49&&Ev<1.51&&pdgi==211&&nipip==1&&nipim==0&&nipi0==0","GOFF");
177  }
178  else if(single_pion_sources==2) {
179  // all resonances only
180  chain->Draw("Ws>>hst_dsig_dW_1000MeV","cc&&res&&Ev>0.99&&Ev<1.01&&pdgi==211&&nipip==1&&nipim==0&&nipi0==0","GOFF");
181  chain->Draw("Ws>>hst_dsig_dW_1500MeV","cc&&res&&Ev>1.49&&Ev<1.51&&pdgi==211&&nipip==1&&nipim==0&&nipi0==0","GOFF");
182  }
183  }
184 
185  //
186  // normalize
187  //
188  double norm_cc1pip_1000MeV = hst_dsig_dW_1000MeV -> Integral("width") / sig_cc1pip_1000MeV;
189  double norm_cc1pip_1500MeV = hst_dsig_dW_1500MeV -> Integral("width") / sig_cc1pip_1500MeV;
190 
191  if (norm_cc1pip_1000MeV > 0) hst_dsig_dW_1000MeV -> Scale(1./norm_cc1pip_1000MeV);
192  if (norm_cc1pip_1500MeV > 0) hst_dsig_dW_1500MeV -> Scale(1./norm_cc1pip_1500MeV);
193 
194  for(int i=1; i <= hst_dsig_dW_1000MeV->GetNbinsX(); i++) {
195 
196  double W = hst_dsig_dW_1000MeV -> GetBinCenter(i);
197 
198  double dsig_dW_1000MeV = hst_dsig_dW_1000MeV -> GetBinContent(i);
199  double dsig_dW_1500MeV = hst_dsig_dW_1500MeV -> GetBinContent(i);
200 
201  dsig_dW_1000MeV = TMath::Max(0., dsig_dW_1000MeV);
202  dsig_dW_1500MeV = TMath::Max(0., dsig_dW_1500MeV);
203 
204  out_stream << setw(15) << W
205  << setw(15) << dsig_dW_1000MeV
206  << setw(15) << dsig_dW_1500MeV
207  << endl;
208  }
209 
210  out_stream.close();
211 
212  // visual inspection
213  //TCanvas * c1 = new TCanvas("c1","",20,20,500,500);
214  //...
215  //c1->Update();
216 
217  delete chain;
218  fsig.Close();
219 }
i
size_t i(0)
Wmax
double Wmax[kNWBins]
Definition: piN_asymmetric_decay.C:40
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
kRunNu1PI4
const int kRunNu1PI4[kNSamples][kNWCur][kNEnergies][kNRunsPerCase]
Definition: nuint09_1pi4.C:32
kNRunsPerCase
const int kNRunsPerCase
Definition: nuint09_1pi4.C:22
dump_to_simple_cpp.W
W
Definition: dump_to_simple_cpp.py:43
kNWCur
const int kNWCur
Definition: nuint09_1pi4.C:20
kLabel
const char * kLabel[kNSamples]
Definition: nuint09_1pi4.C:25
ie
var ie
Definition: slidy.js:14
kNSamples
const int kNSamples
Definition: nuint09_1pi4.C:19
Wmin
double Wmin[kNWBins]
Definition: piN_asymmetric_decay.C:39
kNEnergies
const int kNEnergies
Definition: nuint09_1pi4.C:21

Variable Documentation

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

Definition at line 25 of file nuint09_1pi4.C.

const int kNEnergies = 2

Definition at line 21 of file nuint09_1pi4.C.

const int kNEvtPerRun = 100000

Definition at line 23 of file nuint09_1pi4.C.

const int kNRunsPerCase = 5

Definition at line 22 of file nuint09_1pi4.C.

const int kNSamples = 3

Definition at line 19 of file nuint09_1pi4.C.

const int kNWCur = 1

Definition at line 20 of file nuint09_1pi4.C.

const int kRunNu1PI4[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_1pi4.C.

Generated by   doxygen 1.8.11