Functions | Variables
extract_hadnucleus_xsec.C File Reference
#include <string>
#include <iostream>
#include <TFile.h>
#include <TTree.h>
#include <TH1D.h>
#include <TCanvas.h>
#include <TMath.h>
#include <TStyle.h>

Go to the source code of this file.

Functions

void extract_hadnucleus_xsec (string evtfile)
 

Variables

const double NR = 3.0
 
const double R0 = 1.4
 
const double kemin = 0
 
const double kemax = 1500
 
const double nke = 50
 

Function Documentation

void extract_hadnucleus_xsec ( string  evtfile)

Definition at line 73 of file extract_hadnucleus_xsec.C.

References genie::units::A, om::cout, Divide(), Draw(), allTimeWatchdog::endl, submit_syst::fout, GetXaxis(), GetYaxis(), makeTrainCVSamples::int, kBlue, kemax, kemin, kGreen, kRed, nke, NR, std_candles::pl, R, R0, S, Scale(), SetFillColor(), and Write().

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")
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()

Variable Documentation

const double kemax = 1500
const double kemin = 0
const double nke = 50

Definition at line 71 of file extract_hadnucleus_xsec.C.

Referenced by extract_hadnucleus_xsec().

const double NR = 3.0
const double R0 = 1.4