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