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 
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 }
correl_xv GetYaxis() -> SetDecimals()
void nuint09_coh3(int isample)
Definition: nuint09_coh3.C:59
double Integral(const Spectrum &s, const double pot, int cvnregion)
const int kNRunsPerCase
Definition: nuint09_coh3.C:20
const char * kLabel[kNSamples]
Definition: nuint09_coh3.C:22
string filename
Definition: shutoffs.py:106
const char * label
const int kNWCur
Definition: nuint09_coh3.C:18
correl_xv GetXaxis() -> SetDecimals()
chain
Check that an output directory exists.
const int kRunNu[kNSamples][kNWCur][kNEnergies][kNRunsPerCase]
Definition: nuint09_coh3.C:29
const double j
Definition: BetheBloch.cxx:29
const int kNSamples
Definition: nuint09_coh3.C:17
const int kNEnergies
Definition: nuint09_coh3.C:19
OStream cout
Definition: OStream.cxx:6
simulatedPeak Scale(1/simulationNormalisationFactor)