nuint09_coh2.C
Go to the documentation of this file.
1 //
2 // NuINT09 Conference, Benchmark Calculations (GENIE contribution)
3 //
4 // COH.2:
5 // Cross section as a function of pi0 or pi+ energy at E_nu=0.5, 1.0 and 1.5 GeV
6 //
7 // Inputs:
8 // - sample id: 0 (nu_mu+C12), 1 (nu_mu+O16), 2 (nu_mu+Fe56)
9 //
10 // Costas Andreopoulos, STFC / Rutherford Appleton Laboratory
11 //
12 #include <iomanip>
13 
14 //
15 // consts
16 //
17 const int kNSamples = 3;
18 const int kNWCur = 2;
19 const int kNEnergies = 3;
20 const int kNRunsPerCase = 5;
21 const int kNEvtPerRun = 100000;
22 
23 const char * kLabel[kNSamples] =
24 {
25  /* 0 */ "nu_mu_C12",
26  /* 1 */ "nu_mu_O16",
27  /* 2 */ "nu_mu_Fe56"
28 };
29 
30 const int kRunNu[kNSamples][kNWCur][kNEnergies][kNRunsPerCase] =
31 {
32  /* indices : sample ; cc/nc ; energy */
33  {
34  /* 0,0,0 (nu_mu C12, CC, 0.5 GeV) */ { {100100, 100101, 100102, 100103, 100104},
35  /* 0,0,1 (nu_mu C12, CC, 1.0 GeV) */ {100200, 100201, 100202, 100203, 100204},
36  /* 0,0,2 (nu_mu C12, CC, 1.5 GeV) */ {100300, 100301, 100302, 100303, 100304} },
37  /* 0,1,0 (nu_mu C12, NC, 0.5 GeV) */ { {101100, 101101, 101102, 101103, 101104},
38  /* 0,1,1 (nu_mu C12, NC, 1.0 GeV) */ {101200, 101201, 101202, 101203, 101204},
39  /* 0,1,2 (nu_mu C12, NC, 1.5 GeV) */ {101300, 101301, 101302, 101303, 101304} }
40  },
41  {
42  /* 1,0,0 (nu_mu O16, CC, 0.5 GeV) */ { {110100, 110101, 110102, 110103, 110104},
43  /* 1,0,1 (nu_mu O16, CC, 1.0 GeV) */ {110200, 110201, 110202, 110203, 110204},
44  /* 1,0,2 (nu_mu O16, CC, 1.5 GeV) */ {110300, 110301, 110302, 110303, 110304} },
45  /* 1,1,0 (nu_mu O16, NC, 0.5 GeV) */ { {111100, 111101, 111102, 111103, 111104},
46  /* 1,1,1 (nu_mu O16, NC, 1.0 GeV) */ {111200, 111201, 111202, 111203, 111204},
47  /* 1,1,2 (nu_mu O16, NC, 1.5 GeV) */ {111300, 111301, 111302, 111303, 111304} }
48  },
49  {
50  /* 2,0,0 (nu_mu Fe56, CC, 0.5 GeV) */ { {120100, 120101, 120102, 120103, 120104},
51  /* 2,0,1 (nu_mu Fe56, CC, 1.0 GeV) */ {120200, 120201, 120202, 120203, 120204},
52  /* 2,0,2 (nu_mu Fe56, CC, 1.5 GeV) */ {120300, 120301, 120302, 120303, 120304} },
53  /* 2,1,0 (nu_mu Fe56, NC, 0.5 GeV) */ { {121100, 121101, 121102, 121103, 121104},
54  /* 2,1,1 (nu_mu Fe56, NC, 1.0 GeV) */ {121200, 121201, 121202, 121203, 121204},
55  /* 2,1,2 (nu_mu Fe56, NC, 1.5 GeV) */ {121300, 121301, 121302, 121303, 121304} }
56  }
57 };
58 
59 
60 void nuint09_coh2(int isample)
61 {
62  cout << " ***** running: COH.2" << endl;
63 
64  if(isample<0 || isample >= kNSamples) return;
65 
66  const char * label = kLabel[isample];
67 
68  // get cross section graphs
69 
70  TFile fsig("../sig/splines.root","read");
71  TDirectory * sig_dir = (TDirectory *) fsig.Get(label);
72 
73  TGraph * sig_graph_cohcc = (TGraph*) sig_dir->Get("coh_cc");
74  TGraph * sig_graph_cohnc = (TGraph*) sig_dir->Get("coh_nc");
75 
76  // range & spacing
77 
78  const int nKEpi = 60;
79  const double KEpimin = 0.01;
80  const double KEpimax = 1.50;
81 
82  // create output stream
83 
84  ostringstream out_filename;
85  out_filename << label << ".coh_2.dsig_dKEpi.data";
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 << "# [COH.2]:" << endl;
93  out_stream << "# Cross section as a function of pi0 or pi+ energy at E_nu=0.5, 1.0 and 1.5 GeV" << endl;
94  out_stream << "# " << endl;
95  out_stream << "# Note:" << endl;
96  out_stream << "# - pion energies are _kinetic_ energies " << endl;
97  out_stream << "# - pion kinetic energy KE in GeV, between KEmin = " << KEpimin << " GeV, KEmax = " << KEpimax << " GeV " << endl;
98  out_stream << "# - cross sections in 1E-38 cm^2 / GeV" << endl;
99  out_stream << "# - for coherent scattering we quote _nuclear_ cross section " << endl;
100  out_stream << "# Columns:" << endl;
101  out_stream << "# | KE(pi) | sig(coh; CC; Ev=0.5GeV) | sig(coh; NC; Ev=0.5GeV) | sig(coh; CC; Ev=1.0GeV) | sig(coh; NC; Ev=1.0GeV) | sig(coh; CC; Ev=1.5GeV) | sig(coh; NC; Ev=1.5GeV) " << endl;
102 
103  out_stream << std::fixed << setprecision(6);
104 
105  //
106  // load event data
107  //
108 
109  TChain * chain = new TChain("gst");
110 
111  // loop over CC/NC cases
112  for(int iwkcur=0; iwkcur<kNWCur; iwkcur++) {
113  // loop over energies
114  for(int ie=0; ie<kNEnergies; ie++) {
115  // loop over runs for current case
116  for(int ir=0; ir<kNRunsPerCase; ir++) {
117  // build filename
118  ostringstream filename;
119  int run_number = kRunNu[isample][iwkcur][ie][ir];
120  filename << "../gst/gntp." << run_number << ".gst.root";
121  // add to chain
122  cout << "Adding " << filename.str() << " to event chain" << endl;
123  chain->Add(filename.str().c_str());
124  }
125  }
126  }
127 
128 
129  //
130  // get COH pi cross sections at given energies for normalization purposes
131  //
132  double sig_cohcc_0500MeV = sig_graph_cohcc->Eval(0.5);
133  double sig_cohnc_0500MeV = sig_graph_cohnc->Eval(0.5);
134  double sig_cohcc_1000MeV = sig_graph_cohcc->Eval(1.0);
135  double sig_cohnc_1000MeV = sig_graph_cohnc->Eval(1.0);
136  double sig_cohcc_1500MeV = sig_graph_cohcc->Eval(1.5);
137  double sig_cohnc_1500MeV = sig_graph_cohnc->Eval(1.5);
138 
139  //
140  // book histograms
141  //
142  TH1D * hst_dsig_dKEpi_cohcc_0500MeV = new TH1D("hst_dsig_dKEpi_cohcc_0500MeV","dsig/dKEpi, COH CC, Enu=0.5 GeV", nKEpi, KEpimin, KEpimax);
143  TH1D * hst_dsig_dKEpi_cohnc_0500MeV = new TH1D("hst_dsig_dKEpi_cohnc_0500MeV","dsig/dKEpi, COH NC, Enu=0.5 GeV", nKEpi, KEpimin, KEpimax);
144  TH1D * hst_dsig_dKEpi_cohcc_1000MeV = new TH1D("hst_dsig_dKEpi_cohcc_1000MeV","dsig/dKEpi, COH CC, Enu=1.0 GeV", nKEpi, KEpimin, KEpimax);
145  TH1D * hst_dsig_dKEpi_cohnc_1000MeV = new TH1D("hst_dsig_dKEpi_cohnc_1000MeV","dsig/dKEpi, COH NC, Enu=1.0 GeV", nKEpi, KEpimin, KEpimax);
146  TH1D * hst_dsig_dKEpi_cohcc_1500MeV = new TH1D("hst_dsig_dKEpi_cohcc_1500MeV","dsig/dKEpi, COH CC, Enu=1.5 GeV", nKEpi, KEpimin, KEpimax);
147  TH1D * hst_dsig_dKEpi_cohnc_1500MeV = new TH1D("hst_dsig_dKEpi_cohnc_1500MeV","dsig/dKEpi, COH NC, Enu=1.5 GeV", nKEpi, KEpimin, KEpimax);
148 
149  //
150  // fill histograms
151  //
152  chain->Draw("(Ef-0.13957)>>hst_dsig_dKEpi_cohcc_0500MeV","coh&&cc&&Ev>0.49&&Ev<0.51&&pdgf==211","GOFF");
153  chain->Draw("(Ef-0.13498)>>hst_dsig_dKEpi_cohnc_0500MeV","coh&&nc&&Ev>0.49&&Ev<0.51&&pdgf==111","GOFF");
154  chain->Draw("(Ef-0.13957)>>hst_dsig_dKEpi_cohcc_1000MeV","coh&&cc&&Ev>0.99&&Ev<1.01&&pdgf==211","GOFF");
155  chain->Draw("(Ef-0.13498)>>hst_dsig_dKEpi_cohnc_1000MeV","coh&&nc&&Ev>0.99&&Ev<1.01&&pdgf==111","GOFF");
156  chain->Draw("(Ef-0.13957)>>hst_dsig_dKEpi_cohcc_1500MeV","coh&&cc&&Ev>1.49&&Ev<1.51&&pdgf==211","GOFF");
157  chain->Draw("(Ef-0.13498)>>hst_dsig_dKEpi_cohnc_1500MeV","coh&&nc&&Ev>1.49&&Ev<1.51&&pdgf==111","GOFF");
158 
159  //
160  // normalize
161  //
162  double norm_cohcc_0500MeV = hst_dsig_dKEpi_cohcc_0500MeV -> Integral("width") / sig_cohcc_0500MeV;
163  double norm_cohnc_0500MeV = hst_dsig_dKEpi_cohnc_0500MeV -> Integral("width") / sig_cohnc_0500MeV;
164  double norm_cohcc_1000MeV = hst_dsig_dKEpi_cohcc_1000MeV -> Integral("width") / sig_cohcc_1000MeV;
165  double norm_cohnc_1000MeV = hst_dsig_dKEpi_cohnc_1000MeV -> Integral("width") / sig_cohnc_1000MeV;
166  double norm_cohcc_1500MeV = hst_dsig_dKEpi_cohcc_1500MeV -> Integral("width") / sig_cohcc_1500MeV;
167  double norm_cohnc_1500MeV = hst_dsig_dKEpi_cohnc_1500MeV -> Integral("width") / sig_cohnc_1500MeV;
168 
169  if (norm_cohcc_0500MeV > 0) hst_dsig_dKEpi_cohcc_0500MeV -> Scale(1./norm_cohcc_0500MeV);
170  if (norm_cohnc_0500MeV > 0) hst_dsig_dKEpi_cohnc_0500MeV -> Scale(1./norm_cohnc_0500MeV);
171  if (norm_cohcc_1000MeV > 0) hst_dsig_dKEpi_cohcc_1000MeV -> Scale(1./norm_cohcc_1000MeV);
172  if (norm_cohnc_1000MeV > 0) hst_dsig_dKEpi_cohnc_1000MeV -> Scale(1./norm_cohnc_1000MeV);
173  if (norm_cohcc_1500MeV > 0) hst_dsig_dKEpi_cohcc_1500MeV -> Scale(1./norm_cohcc_1500MeV);
174  if (norm_cohnc_1500MeV > 0) hst_dsig_dKEpi_cohnc_1500MeV -> Scale(1./norm_cohnc_1500MeV);
175 
176  for(int i = 1; i <= hst_dsig_dKEpi_cohnc_1500MeV->GetNbinsX(); i++) {
177 
178  double Epi = hst_dsig_dKEpi_cohnc_1500MeV->GetBinCenter(i);
179 
180  double dsig_dKEpi_cohcc_0500MeV = hst_dsig_dKEpi_cohcc_0500MeV -> GetBinContent(i);
181  double dsig_dKEpi_cohnc_0500MeV = hst_dsig_dKEpi_cohnc_0500MeV -> GetBinContent(i);
182  double dsig_dKEpi_cohcc_1000MeV = hst_dsig_dKEpi_cohcc_1000MeV -> GetBinContent(i);
183  double dsig_dKEpi_cohnc_1000MeV = hst_dsig_dKEpi_cohnc_1000MeV -> GetBinContent(i);
184  double dsig_dKEpi_cohcc_1500MeV = hst_dsig_dKEpi_cohcc_1500MeV -> GetBinContent(i);
185  double dsig_dKEpi_cohnc_1500MeV = hst_dsig_dKEpi_cohnc_1500MeV -> GetBinContent(i);
186 
187  dsig_dKEpi_cohcc_0500MeV = TMath::Max(0., dsig_dKEpi_cohcc_0500MeV);
188  dsig_dKEpi_cohnc_0500MeV = TMath::Max(0., dsig_dKEpi_cohnc_0500MeV);
189  dsig_dKEpi_cohcc_1000MeV = TMath::Max(0., dsig_dKEpi_cohcc_1000MeV);
190  dsig_dKEpi_cohnc_1000MeV = TMath::Max(0., dsig_dKEpi_cohnc_1000MeV);
191  dsig_dKEpi_cohcc_1500MeV = TMath::Max(0., dsig_dKEpi_cohcc_1500MeV);
192  dsig_dKEpi_cohnc_1500MeV = TMath::Max(0., dsig_dKEpi_cohnc_1500MeV);
193 
194  out_stream << setw(15) << Epi
195  << setw(15) << dsig_dKEpi_cohcc_0500MeV
196  << setw(15) << dsig_dKEpi_cohnc_0500MeV
197  << setw(15) << dsig_dKEpi_cohcc_1000MeV
198  << setw(15) << dsig_dKEpi_cohnc_1000MeV
199  << setw(15) << dsig_dKEpi_cohcc_1500MeV
200  << setw(15) << dsig_dKEpi_cohnc_1500MeV
201  << endl;
202  }
203 
204  out_stream.close();
205 
206  // visual inspection
207 
208  TCanvas * c1 = new TCanvas("c1","",20,20,500,500);
209  hst_dsig_dKEpi_cohcc_1500MeV->Draw();
210  hst_dsig_dKEpi_cohnc_1500MeV->Draw("same");
211  c1->Update();
212 
213  delete chain;
214 }
kNRunsPerCase
const int kNRunsPerCase
Definition: nuint09_coh2.C:20
kNWCur
const int kNWCur
Definition: nuint09_coh2.C:18
i
size_t i(0)
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
kLabel
const char * kLabel[kNSamples]
Definition: nuint09_coh2.C:23
kRunNu
const int kRunNu[kNSamples][kNWCur][kNEnergies][kNRunsPerCase]
Definition: nuint09_coh2.C:30
kNSamples
const int kNSamples
Definition: nuint09_coh2.C:17
ie
var ie
Definition: slidy.js:14
kNEvtPerRun
const int kNEvtPerRun
Definition: nuint09_coh2.C:21
kNEnergies
const int kNEnergies
Definition: nuint09_coh2.C:19
nuint09_coh2
void nuint09_coh2(int isample)
Definition: nuint09_coh2.C:60
Generated by   doxygen 1.8.11