hA_basic_analysis.C

Go to the documentation of this file.

/*
 * Basic hadron+nucleon event sample analysis.
 *
 * Extract 
 *  - sigma(total)
 *  - sigma(reaction)
 *  - sigma(charge_exchange)
 *  - sigma(elastic) 
 *  - sigma(inelastic) 
 *  - sigma(absorption)
 *  - sigma(pi_production) 
 * from a `MC experiment' shooting a hadron beam towards a nuclear target.
 * The sample can be generated using GENIE's ghAevgen utility.
 * The summary ntuple input here can be generated by analyzing ghAevgen's output
 * event file using GENIE's gntpc utility (use -f ghA)
 *
 * C.Andreopoulos <costas.andreopoulos \at stfc.ac.uk>
 *
 */

#ifndef __CINT__

#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;
const double R0      = 1.4;
const double kemin   = 0;
const double kemax   = 1500;
const double nke     = 50;

void hA_basic_analysis(string evtfile, int A)
{
  TFile * finp = new TFile(evtfile.c_str(), "READ");
  TTree * ghA = (TTree*)finp->Get("ghA");
  if(!ghA) return;

  TH1D * pip_nall        = new TH1D("pip_nall",        "",                                   nke,kemin,kemax);
  TH1D * pip_sigma_tot   = new TH1D("pip_sigma_tot",   Form("pi+ A(=%d) total",A),           nke,kemin,kemax);
  TH1D * pip_sigma_reac  = new TH1D("pip_sigma_reac",  Form("pi+ A(=%d) reaction",A),        nke,kemin,kemax);
  TH1D * pip_sigma_cex   = new TH1D("pip_sigma_cex",   Form("pi+ A(=%d) charge exchange",A), nke,kemin,kemax);
  TH1D * pip_sigma_el    = new TH1D("pip_sigma_el",    Form("pi+ A(=%d) elastic",A),         nke,kemin,kemax);
  TH1D * pip_sigma_inel  = new TH1D("pip_sigma_inel",  Form("pi+ A(=%d) inelastic",A),       nke,kemin,kemax);
  TH1D * pip_sigma_abs   = new TH1D("pip_sigma_abs",   Form("pi+ A(=%d) absorption",A),      nke,kemin,kemax);
  TH1D * pip_sigma_pipro = new TH1D("pip_sigma_pipro", Form("pi+ A(=%d) pi production",A),   nke,kemin,kemax);

  ghA->Draw("1000*ke>>pip_nall","","");
  ghA->Draw("1000*ke>>pip_sigma_tot",  Form("A==%d && probe==%d && fsi>0",A,211),           "goff");
  ghA->Draw("1000*ke>>pip_sigma_reac", Form("A==%d && probe==%d && fsi>0 && fsi!=2",A,211), "goff");
  ghA->Draw("1000*ke>>pip_sigma_cex",  Form("A==%d && probe==%d && fsi==1",A,211),          "goff");
  ghA->Draw("1000*ke>>pip_sigma_el",   Form("A==%d && probe==%d && fsi==2",A,211),          "goff");
  ghA->Draw("1000*ke>>pip_sigma_inel", Form("A==%d && probe==%d && fsi==3",A,211),          "goff");
  ghA->Draw("1000*ke>>pip_sigma_abs",  Form("A==%d && probe==%d && fsi>=4 && fsi<=9",A,211),"goff");
  ghA->Draw("1000*ke>>pip_sigma_pipro",Form("A==%d && probe==%d && fsi>9",A,211),           "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);

  pip_sigma_tot   -> Divide (pip_nall);  pip_sigma_tot   -> Scale (S);  pip_sigma_tot   -> Sumw2 ();
  pip_sigma_reac  -> Divide (pip_nall);  pip_sigma_reac  -> Scale (S);  pip_sigma_reac  -> Sumw2 ();
  pip_sigma_cex   -> Divide (pip_nall);  pip_sigma_cex   -> Scale (S);  pip_sigma_cex   -> Sumw2 ();
  pip_sigma_el    -> Divide (pip_nall);  pip_sigma_el    -> Scale (S);  pip_sigma_el    -> Sumw2 ();
  pip_sigma_inel  -> Divide (pip_nall);  pip_sigma_inel  -> Scale (S);  pip_sigma_inel  -> Sumw2 ();
  pip_sigma_abs   -> Divide (pip_nall);  pip_sigma_abs   -> Scale (S);  pip_sigma_abs   -> Sumw2 ();
  pip_sigma_pipro -> Divide (pip_nall);  pip_sigma_pipro -> Scale (S);  pip_sigma_pipro -> Sumw2 ();

  pip_sigma_tot   -> SetFillColor(kBlue);  
  pip_sigma_reac  -> SetFillColor(kGreen); 
  pip_sigma_cex   -> SetFillColor(kRed);
  pip_sigma_el    -> SetFillColor(kRed);   
  pip_sigma_inel  -> SetFillColor(kRed); 
  pip_sigma_abs   -> SetFillColor(kRed);
  pip_sigma_pipro -> SetFillColor(kRed);

  pip_sigma_tot  ->GetXaxis()->SetTitle("#pi^{+} KE (MeV)");  pip_sigma_tot  ->GetYaxis()->SetTitle("#sigma (mb)");
  pip_sigma_reac ->GetXaxis()->SetTitle("#pi^{+} KE (MeV)");  pip_sigma_reac ->GetYaxis()->SetTitle("#sigma (mb)");
  pip_sigma_cex  ->GetXaxis()->SetTitle("#pi^{+} KE (MeV)");  pip_sigma_cex  ->GetYaxis()->SetTitle("#sigma (mb)");
  pip_sigma_el   ->GetXaxis()->SetTitle("#pi^{+} KE (MeV)");  pip_sigma_el   ->GetYaxis()->SetTitle("#sigma (mb)");
  pip_sigma_inel ->GetXaxis()->SetTitle("#pi^{+} KE (MeV)");  pip_sigma_inel ->GetYaxis()->SetTitle("#sigma (mb)");
  pip_sigma_abs  ->GetXaxis()->SetTitle("#pi^{+} KE (MeV)");  pip_sigma_abs  ->GetYaxis()->SetTitle("#sigma (mb)");
  pip_sigma_pipro->GetXaxis()->SetTitle("#pi^{+} KE (MeV)");  pip_sigma_pipro->GetYaxis()->SetTitle("#sigma (mb)");

  gStyle->SetOptStat(0);

  TCanvas * ctot = new TCanvas("ctot","",10,10,500,500);
  pip_sigma_tot  -> Draw("E4");
  pip_sigma_reac -> Draw("E4SAME");
  ctot->Update();
  TCanvas * cfates = new TCanvas("cfates","",110,110,600,600);
  cfates->Divide(2,2);
  cfates->cd(1);
  pip_sigma_cex   -> Draw ("E4");
  cfates->cd(2);
  pip_sigma_inel  -> Draw ("E4");
  cfates->cd(3);
  pip_sigma_abs   -> Draw ("E4");
  cfates->cd(4);
  pip_sigma_pipro -> Draw ("E4");
  cfates->Update();

  TFile fout("hAana.root","RECREATE");
  pip_sigma_tot   -> Write("pip_sigma_tot");
  pip_sigma_reac  -> Write("pip_sigma_reac");
  pip_sigma_cex   -> Write("pip_sigma_cex");
  pip_sigma_el    -> Write("pip_sigma_el");
  pip_sigma_inel  -> Write("pip_sigma_inel");
  pip_sigma_abs   -> Write("pip_sigma_abs");
  pip_sigma_pipro -> Write("pip_sigma_pipro");
  fout.Close();

//  finp.Close();
}