#include <string>
#include <iostream>
#include <TFile.h>
#include <TTree.h>
#include <TH1D.h>
#include <TCanvas.h>
#include <TMath.h>
#include <TStyle.h>

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

 {
   TFile * finp   = new TFile(evtfile.c_str(), "READ");
   TTree * ginuke = (TTree*)finp->Get("ginuke");
   if(!ginuke) return;
 
   // peek first event to get incident hadron and target nucleus PDG codes
   ginuke->Draw("probe:A","iev==0","goff");
   int probe = (int) ginuke->GetV1()[0];
   int A     = (int) ginuke->GetV2()[0];
 
   cout << "Incident hadron    : " << probe << endl;
   cout << "Target mass number : " << A     << endl;
 
   string probe_label = "";
   switch(probe) {
     case (   22) : probe_label = "gamma"; break;
     case (  211) : probe_label = "pi+";   break;
     case ( -211) : probe_label = "pi-";   break;
     case (  111) : probe_label = "pi0";   break;
     case ( 2212) : probe_label = "p";     break;
     case ( 2112) : probe_label = "n";     break;
     default      : probe_label = "X";     break;
   }
   const char * pl = probe_label.c_str();
 
   TH1D * nall        = new TH1D("nall",        "",                                       nke,kemin,kemax);
   TH1D * sigma_tot   = new TH1D("sigma_tot",   Form("%s A(=%d) total",           pl,A),  nke,kemin,kemax);
   TH1D * sigma_reac  = new TH1D("sigma_reac",  Form("%s A(=%d) reaction",        pl,A),  nke,kemin,kemax);
   TH1D * sigma_cex   = new TH1D("sigma_cex",   Form("%s A(=%d) charge exchange", pl,A),  nke,kemin,kemax);
   TH1D * sigma_el    = new TH1D("sigma_el",    Form("%s A(=%d) elastic",         pl,A),  nke,kemin,kemax);
   TH1D * sigma_inel  = new TH1D("sigma_inel",  Form("%s A(=%d) inelastic",       pl,A),  nke,kemin,kemax);
   TH1D * sigma_abs   = new TH1D("sigma_abs",   Form("%s A(=%d) absorption",      pl,A),  nke,kemin,kemax);
   TH1D * sigma_pipro = new TH1D("sigma_pipro", Form("%s A(=%d) pi production",   pl,A),  nke,kemin,kemax);
 
   ginuke->Draw("1000*ke>>nall","","");
   ginuke->Draw("1000*ke>>sigma_tot",  Form("A==%d && probe==%d && probe_fsi>0",A,probe),                 "goff");
   ginuke->Draw("1000*ke>>sigma_reac", Form("A==%d && probe==%d && probe_fsi>0 && probe_fsi!=2",A,probe), "goff");
   ginuke->Draw("1000*ke>>sigma_cex",  Form("A==%d && probe==%d && probe_fsi==1",A,probe),                "goff");
   ginuke->Draw("1000*ke>>sigma_el",   Form("A==%d && probe==%d && probe_fsi==2",A,probe),                "goff");
   ginuke->Draw("1000*ke>>sigma_inel", Form("A==%d && probe==%d && probe_fsi==3",A,probe),                "goff");
   ginuke->Draw("1000*ke>>sigma_abs",  Form("A==%d && probe==%d && probe_fsi>=4 && probe_fsi<=9",A,probe),"goff");
   ginuke->Draw("1000*ke>>sigma_pipro",Form("A==%d && probe==%d && probe_fsi>9",A,probe),                 "goff");
 
   const double fm2tomb = 10.; // fm^2 -> mb
 
   double R  = NR * R0 * TMath::Power((double)A,0.3333);
   double S  = fm2tomb * TMath::Pi() * TMath::Power(R,2);
 
   sigma_tot   -> Divide (nall);  sigma_tot   -> Scale (S);  sigma_tot   -> Sumw2 ();
   sigma_reac  -> Divide (nall);  sigma_reac  -> Scale (S);  sigma_reac  -> Sumw2 ();
   sigma_cex   -> Divide (nall);  sigma_cex   -> Scale (S);  sigma_cex   -> Sumw2 ();
   sigma_el    -> Divide (nall);  sigma_el    -> Scale (S);  sigma_el    -> Sumw2 ();
   sigma_inel  -> Divide (nall);  sigma_inel  -> Scale (S);  sigma_inel  -> Sumw2 ();
   sigma_abs   -> Divide (nall);  sigma_abs   -> Scale (S);  sigma_abs   -> Sumw2 ();
   sigma_pipro -> Divide (nall);  sigma_pipro -> Scale (S);  sigma_pipro -> Sumw2 ();
 
   sigma_tot   -> SetFillColor ( kBlue  );  
   sigma_reac  -> SetFillColor ( kGreen ); 
   sigma_cex   -> SetFillColor ( kRed   );
   sigma_el    -> SetFillColor ( kRed   );   
   sigma_inel  -> SetFillColor ( kRed   ); 
   sigma_abs   -> SetFillColor ( kRed   );
   sigma_pipro -> SetFillColor ( kRed   );
 
   sigma_tot   -> GetXaxis()->SetTitle(Form("%s KE (MeV)", pl));  
   sigma_tot   -> GetYaxis()->SetTitle("#sigma_{tot} (mb)");
   sigma_reac  -> GetXaxis()->SetTitle(Form("%s KE (MeV)", pl));  
   sigma_reac  -> GetYaxis()->SetTitle("#sigma_{reac} (mb)");
   sigma_cex   -> GetXaxis()->SetTitle(Form("%s KE (MeV)", pl));  
   sigma_cex   -> GetYaxis()->SetTitle("#sigma_{1cex} (mb)");
   sigma_el    -> GetXaxis()->SetTitle(Form("%s KE (MeV)", pl)); 
   sigma_el    -> GetYaxis()->SetTitle("#sigma_{elas} (mb)");
   sigma_inel  -> GetXaxis()->SetTitle(Form("%s KE (MeV)", pl));  
   sigma_inel  -> GetYaxis()->SetTitle("#sigma_{inel} (mb)");
   sigma_abs   -> GetXaxis()->SetTitle(Form("%s KE (MeV)", pl)); 
   sigma_abs   -> GetYaxis()->SetTitle("#sigma_{abs} (mb)");
   sigma_pipro -> GetXaxis()->SetTitle(Form("%s KE (MeV)", pl)); 
   sigma_pipro -> GetYaxis()->SetTitle("#sigma_{#pi prod.} (mb)");
 
   gStyle->SetOptStat(0);
 
   TCanvas * ctot = new TCanvas("ctot","",10,10,500,500);
   sigma_tot  -> Draw("E4");
   sigma_reac -> Draw("E4SAME");
   ctot->Update();
   TCanvas * cfates = new TCanvas("cfates","",110,110,600,600);
   cfates->Divide(2,2);
   cfates->cd(1);
   sigma_cex   -> Draw ("E4");
   cfates->cd(2);
   sigma_inel  -> Draw ("E4");
   cfates->cd(3);
   sigma_abs   -> Draw ("E4");
   cfates->cd(4);
   sigma_pipro -> Draw ("E4");
   cfates->Update();
 
   TFile fout("hadxsec.root","RECREATE");
   sigma_tot   -> Write( "sigma_tot"  );
   sigma_reac  -> Write( "sigma_reac" );
   sigma_cex   -> Write( "sigma_cex"  );
   sigma_el    -> Write( "sigma_el"   );
   sigma_inel  -> Write( "sigma_inel" );
   sigma_abs   -> Write( "sigma_abs"  );
   sigma_pipro -> Write( "sigma_pipro");
   fout.Close();
 
 //  finp.Close();
 }