extract_hadnucleus_xsec.C

Go to the documentation of this file.

/*
 * ....................................................................................
 *
 * Basic hadron+nucleon event sample analysis.
 *
 * ....................................................................................
 *
 * Summary:
 *
 * This macro extracts: 
 *  - sigma(total)
 *  - sigma(reaction)
 *  - sigma(charge_exchange)
 *  - sigma(elastic) 
 *  - sigma(inelastic) 
 *  - sigma(absorption)
 *  - sigma(pi_production) 
 * from an `MC experiment' shooting a hadron beam towards a nuclear target.
 *
 * ....................................................................................
 *
 * Detailed instructions:
 *
 * The event sample can be generated using GENIE's gevgen_hadron utility.
 * The summary ntuple input here can be generated by analyzing gevgen_hadron's output
 * event file using GENIE's gntpc utility (use -f ginuke)
 *
 * For example, to get the pi+ Fe56 cross sections, generate 100k events with
 * incident pion kinetic energies uniformly distributed between 0 and 1500 MeV:
 * % gevgen_hadron -p 211 -t 1000080160 -k 0.0,1.5 -n 100000 -m hA -r 100 -o pipFe56 
 *
 * The output GHEP event tree will be saved in pipFe56.100.ghep.root
 *
 * Then analyze the output pi+ Fe56 event file to obtain an intranuke summary ntuple:
 * % gntpc -i gntp.0.ghep.root -f ginuke 
 *
 * The output summary ntuple will be saved in pipFe56.100.ginuke.root
 *
 * Then, run this macro using the output summary ntuple generated at the previous step:
 * % root 
 * root[0] .L extract_hadnucleus_xsec.C
 * root[0] extract_hadnucleus_xsec("pipFe56.100.ginuke.root")
 *
 * ....................................................................................
 *
 * C.Andreopoulos <costas.andreopoulos \at stfc.ac.uk>
 *
 */

#ifndef __CINT__

#include <string>
#include <iostream>
#include <string>

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

using namespace std;

#endif

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

void extract_hadnucleus_xsec(string evtfile)
{
  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();
}