extract_hadnucleus_xsec.C
Go to the documentation of this file.
1 /*
2  * ....................................................................................
3  *
4  * Basic hadron+nucleon event sample analysis.
5  *
6  * ....................................................................................
7  *
8  * Summary:
9  *
10  * This macro extracts:
11  * - sigma(total)
12  * - sigma(reaction)
13  * - sigma(charge_exchange)
14  * - sigma(elastic)
15  * - sigma(inelastic)
16  * - sigma(absorption)
17  * - sigma(pi_production)
18  * from an `MC experiment' shooting a hadron beam towards a nuclear target.
19  *
20  * ....................................................................................
21  *
22  * Detailed instructions:
23  *
24  * The event sample can be generated using GENIE's gevgen_hadron utility.
25  * The summary ntuple input here can be generated by analyzing gevgen_hadron's output
26  * event file using GENIE's gntpc utility (use -f ginuke)
27  *
28  * For example, to get the pi+ Fe56 cross sections, generate 100k events with
29  * incident pion kinetic energies uniformly distributed between 0 and 1500 MeV:
30  * % gevgen_hadron -p 211 -t 1000080160 -k 0.0,1.5 -n 100000 -m hA -r 100 -o pipFe56
31  *
32  * The output GHEP event tree will be saved in pipFe56.100.ghep.root
33  *
34  * Then analyze the output pi+ Fe56 event file to obtain an intranuke summary ntuple:
35  * % gntpc -i gntp.0.ghep.root -f ginuke
36  *
37  * The output summary ntuple will be saved in pipFe56.100.ginuke.root
38  *
39  * Then, run this macro using the output summary ntuple generated at the previous step:
40  * % root
41  * root[0] .L extract_hadnucleus_xsec.C
42  * root[0] extract_hadnucleus_xsec("pipFe56.100.ginuke.root")
43  *
44  * ....................................................................................
45  *
46  * C.Andreopoulos <costas.andreopoulos \at stfc.ac.uk>
47  *
48  */
49 
50 #ifndef __CINT__
51 
52 #include <string>
53 #include <iostream>
54 #include <string>
55 
56 #include <TFile.h>
57 #include <TTree.h>
58 #include <TH1D.h>
59 #include <TCanvas.h>
60 #include <TMath.h>
61 #include <TStyle.h>
62 
63 using namespace std;
64 
65 #endif
66 
67 const double NR = 3.0;
68 const double R0 = 1.4;
69 const double kemin = 0; // MeV
70 const double kemax = 1500; // MeV
71 const double nke = 50;
72 
73 void extract_hadnucleus_xsec(string evtfile)
74 {
75  TFile * finp = new TFile(evtfile.c_str(), "READ");
76  TTree * ginuke = (TTree*)finp->Get("ginuke");
77  if(!ginuke) return;
78 
79  // peek first event to get incident hadron and target nucleus PDG codes
80  ginuke->Draw("probe:A","iev==0","goff");
81  int probe = (int) ginuke->GetV1()[0];
82  int A = (int) ginuke->GetV2()[0];
83 
84  cout << "Incident hadron : " << probe << endl;
85  cout << "Target mass number : " << A << endl;
86 
87  string probe_label = "";
88  switch(probe) {
89  case ( 22) : probe_label = "gamma"; break;
90  case ( 211) : probe_label = "pi+"; break;
91  case ( -211) : probe_label = "pi-"; break;
92  case ( 111) : probe_label = "pi0"; break;
93  case ( 2212) : probe_label = "p"; break;
94  case ( 2112) : probe_label = "n"; break;
95  default : probe_label = "X"; break;
96  }
97  const char * pl = probe_label.c_str();
98 
99  TH1D * nall = new TH1D("nall", "", nke,kemin,kemax);
100  TH1D * sigma_tot = new TH1D("sigma_tot", Form("%s A(=%d) total", pl,A), nke,kemin,kemax);
101  TH1D * sigma_reac = new TH1D("sigma_reac", Form("%s A(=%d) reaction", pl,A), nke,kemin,kemax);
102  TH1D * sigma_cex = new TH1D("sigma_cex", Form("%s A(=%d) charge exchange", pl,A), nke,kemin,kemax);
103  TH1D * sigma_el = new TH1D("sigma_el", Form("%s A(=%d) elastic", pl,A), nke,kemin,kemax);
104  TH1D * sigma_inel = new TH1D("sigma_inel", Form("%s A(=%d) inelastic", pl,A), nke,kemin,kemax);
105  TH1D * sigma_abs = new TH1D("sigma_abs", Form("%s A(=%d) absorption", pl,A), nke,kemin,kemax);
106  TH1D * sigma_pipro = new TH1D("sigma_pipro", Form("%s A(=%d) pi production", pl,A), nke,kemin,kemax);
107 
108  ginuke->Draw("1000*ke>>nall","","");
109  ginuke->Draw("1000*ke>>sigma_tot", Form("A==%d && probe==%d && probe_fsi>0",A,probe), "goff");
110  ginuke->Draw("1000*ke>>sigma_reac", Form("A==%d && probe==%d && probe_fsi>0 && probe_fsi!=2",A,probe), "goff");
111  ginuke->Draw("1000*ke>>sigma_cex", Form("A==%d && probe==%d && probe_fsi==1",A,probe), "goff");
112  ginuke->Draw("1000*ke>>sigma_el", Form("A==%d && probe==%d && probe_fsi==2",A,probe), "goff");
113  ginuke->Draw("1000*ke>>sigma_inel", Form("A==%d && probe==%d && probe_fsi==3",A,probe), "goff");
114  ginuke->Draw("1000*ke>>sigma_abs", Form("A==%d && probe==%d && probe_fsi>=4 && probe_fsi<=9",A,probe),"goff");
115  ginuke->Draw("1000*ke>>sigma_pipro",Form("A==%d && probe==%d && probe_fsi>9",A,probe), "goff");
116 
117  const double fm2tomb = 10.; // fm^2 -> mb
118 
119  double R = NR * R0 * TMath::Power((double)A,0.3333);
120  double S = fm2tomb * TMath::Pi() * TMath::Power(R,2);
121 
122  sigma_tot -> Divide (nall); sigma_tot -> Scale (S); sigma_tot -> Sumw2 ();
123  sigma_reac -> Divide (nall); sigma_reac -> Scale (S); sigma_reac -> Sumw2 ();
124  sigma_cex -> Divide (nall); sigma_cex -> Scale (S); sigma_cex -> Sumw2 ();
125  sigma_el -> Divide (nall); sigma_el -> Scale (S); sigma_el -> Sumw2 ();
126  sigma_inel -> Divide (nall); sigma_inel -> Scale (S); sigma_inel -> Sumw2 ();
127  sigma_abs -> Divide (nall); sigma_abs -> Scale (S); sigma_abs -> Sumw2 ();
128  sigma_pipro -> Divide (nall); sigma_pipro -> Scale (S); sigma_pipro -> Sumw2 ();
129 
130  sigma_tot -> SetFillColor ( kBlue );
131  sigma_reac -> SetFillColor ( kGreen );
132  sigma_cex -> SetFillColor ( kRed );
133  sigma_el -> SetFillColor ( kRed );
134  sigma_inel -> SetFillColor ( kRed );
135  sigma_abs -> SetFillColor ( kRed );
136  sigma_pipro -> SetFillColor ( kRed );
137 
138  sigma_tot -> GetXaxis()->SetTitle(Form("%s KE (MeV)", pl));
139  sigma_tot -> GetYaxis()->SetTitle("#sigma_{tot} (mb)");
140  sigma_reac -> GetXaxis()->SetTitle(Form("%s KE (MeV)", pl));
141  sigma_reac -> GetYaxis()->SetTitle("#sigma_{reac} (mb)");
142  sigma_cex -> GetXaxis()->SetTitle(Form("%s KE (MeV)", pl));
143  sigma_cex -> GetYaxis()->SetTitle("#sigma_{1cex} (mb)");
144  sigma_el -> GetXaxis()->SetTitle(Form("%s KE (MeV)", pl));
145  sigma_el -> GetYaxis()->SetTitle("#sigma_{elas} (mb)");
146  sigma_inel -> GetXaxis()->SetTitle(Form("%s KE (MeV)", pl));
147  sigma_inel -> GetYaxis()->SetTitle("#sigma_{inel} (mb)");
148  sigma_abs -> GetXaxis()->SetTitle(Form("%s KE (MeV)", pl));
149  sigma_abs -> GetYaxis()->SetTitle("#sigma_{abs} (mb)");
150  sigma_pipro -> GetXaxis()->SetTitle(Form("%s KE (MeV)", pl));
151  sigma_pipro -> GetYaxis()->SetTitle("#sigma_{#pi prod.} (mb)");
152 
153  gStyle->SetOptStat(0);
154 
155  TCanvas * ctot = new TCanvas("ctot","",10,10,500,500);
156  sigma_tot -> Draw("E4");
157  sigma_reac -> Draw("E4SAME");
158  ctot->Update();
159  TCanvas * cfates = new TCanvas("cfates","",110,110,600,600);
160  cfates->Divide(2,2);
161  cfates->cd(1);
162  sigma_cex -> Draw ("E4");
163  cfates->cd(2);
164  sigma_inel -> Draw ("E4");
165  cfates->cd(3);
166  sigma_abs -> Draw ("E4");
167  cfates->cd(4);
168  sigma_pipro -> Draw ("E4");
169  cfates->Update();
170 
171  TFile fout("hadxsec.root","RECREATE");
172  sigma_tot -> Write( "sigma_tot" );
173  sigma_reac -> Write( "sigma_reac" );
174  sigma_cex -> Write( "sigma_cex" );
175  sigma_el -> Write( "sigma_el" );
176  sigma_inel -> Write( "sigma_inel" );
177  sigma_abs -> Write( "sigma_abs" );
178  sigma_pipro -> Write( "sigma_pipro");
179  fout.Close();
180 
181 // finp.Close();
182 }
tree Draw("slc.nhit")
void extract_hadnucleus_xsec(string evtfile)
ratio_hxv Divide(hxv, goal_hxv)
bin1_2sigma SetFillColor(3)
enum BeamMode kRed
#define S(x, n)
correl_xv GetYaxis() -> SetDecimals()
const double kemax
const double R0
const double NR
correl_xv GetXaxis() -> SetDecimals()
#define R(x)
OStream cout
Definition: OStream.cxx:6
static const double A
Definition: Units.h:82
const double kemin
enum BeamMode kGreen
enum BeamMode kBlue
simulatedPeak Scale(1/simulationNormalisationFactor)
const double nke
gm Write()