macros/diskinematics_analysis/makePlot_truth_LoI_Momentum.C

0001 /**
0002
0003  *\Brief Macro to reproduce ePHENIX LoI Feb. 2014 (arXiv:1402.1209v1) Fig. 2.1, 2.2, & 2.4
0004  *\Author: Thomas Krahulik <thomas.krahulik@stonybrook.edu>
0005
0006  **/
0007
0008 int makePlot_truth_LoI_Momentum(
0009                 const TString infile = "../../data/Sample_DISReco_ep.root",
0010                 //const TString infile = "../../data/EventGenAna_Pythia6_DIS_10x250_100k.root",
0011                 const bool save_figures = true
0012                 )
0013 {
0014   gStyle->SetOptStat(kFALSE);
0015
0016   /*--------------Get Input File--------------*/
0017
0018   TFile *f_truth = new TFile( infile, "OPEN");
0019   TTree* t_truth = (TTree*)f_truth->Get("event_truth");
0020   t_truth->Print();
0021
0022   /*------------------------------------------*/
0023   /*--------------Define Cuts-----------------*/
0024
0025   //TCut Q2_cut = "evtgen_Q2 > 1";
0026
0027   TCut electron_cut = "em_evtgen_pid == 11";
0028   TCut hadron_cut = "abs(em_evtgen_pid) > 100";
0029   TCut proton_cut = "em_evtgen_pid == 2212";
0030   TCut neutron_cut = "em_evtgen_pid == 2112";
0031   TCut Kaon_cut = "abs(em_evtgen_pid) == 321 || em_evtgen_pid == 311";
0032   TCut Pion_charged_cut = "abs(em_evtgen_pid) == 211";
0033   //TCut Pion_cut = "abs(em_evtgen_pid) == 211 || em_evtgen_pid == 111";
0034   TCut photon_cut = "em_evtgen_pid == 22";
0035
0036   TCut mother_cut = "1"; //"p.fParent == 0";
0037   TCut status_cut = "1"; //"p.fKS < 10";
0038
0039   TCut eta_cut_n3n2 = "(-1 * log( tan( em_evtgen_theta / 2.0 ) )) > -3 && (-1 * log( tan( em_evtgen_theta / 2.0 ) )) < -2";
0040   TCut eta_cut_n2n1 = "(-1 * log( tan( em_evtgen_theta / 2.0 ) )) > -2 && (-1 * log( tan( em_evtgen_theta / 2.0 ) )) < -1";
0041   TCut eta_cut_n1z0 = "(-1 * log( tan( em_evtgen_theta / 2.0 ) )) > -1 && (-1 * log( tan( em_evtgen_theta / 2.0 ) )) < -0";
0042
0043   /*------------------------------------------*/
0044   /*-------Momentum vs. Pseudorapidity--------*/
0045   /*------------------------------------------*/
0046
0047   /*------------------------------------------*/
0048   /*---------Electrons (LoI Fig. 2-1)---------*/
0049   /*------------------------------------------*/
0050   TH2F *h_peta_e = new TH2F("h_peta_e", "", 100, -4, 3, 100, 0, 50 ); //250x015 10M
0051
0052   TCanvas *c_peta_e = new TCanvas( "c_peta_e" );
0053   c_peta_e->SetLogz();
0054
0055   t_truth->Draw("em_evtgen_ptotal:(-1 * log( tan( em_evtgen_theta / 2.0 ) ))>>h_peta_e", electron_cut && mother_cut && "evtgen_Q2 > 1", "colz");
0056
0057   h_peta_e->GetXaxis()->SetTitle("Pseudorapidity #eta");
0058   h_peta_e->GetYaxis()->SetTitle("Electron Momentum p_{e-} [GeV]");
0059
0060   if ( save_figures )
0061     {
0062       c_peta_e->Print("plots/Pythia_peta_e.eps");
0063       c_peta_e->Print("plots/Pythia_peta_e.png");
0064     }
0065
0066   /*------------------------------------------*/
0067   /*----------Hadrons (LoI Fig. 2-4)----------*/
0068   /*------------------------------------------*/
0069   TH2F *h_p_eta_h = new TH2F("h_p_eta_h", "", 60,-6,6, 200,0,225); //250x015 10M
0070
0071   TCanvas *c_p_eta_h = new TCanvas( "c_p_eta_h" );
0072   t_truth->Draw("em_evtgen_ptotal:(-1 * log( tan( em_evtgen_theta / 2.0 ) ))>>h_p_eta_h",hadron_cut && "evtgen_Q2 > 1 && 0.01 < evtgen_y < 0.80", "colz");
0073
0074   h_p_eta_h->GetXaxis()->SetTitle("Pseudorapidity #eta");
0075   h_p_eta_h->GetYaxis()->SetTitle("Hadron Momentum p_{Hadron} [GeV]");
0076   c_p_eta_h->SetLogz();
0077
0078   if ( save_figures )
0079     {
0080       c_p_eta_h->Print("plots/Pythia_peta_h.eps");
0081       c_p_eta_h->Print("plots/Pythia_peta_h.png");
0082     }
0083
0084   /*------------Momentum Spectra (Figure 2.2)--------------*/
0085
0086   /*-------- -3 < Eta < -2 ---------*/
0087
0088   TH1F* hp_e_n3n2 = new TH1F("hp_e_n3n2", "dN/dp vs. p", 60, 0, 30);
0089   TH1F* hp_p_n3n2 = new TH1F("hp_p_n3n2", "dN/dp vs. p", 60, 0, 30);
0090   TH1F* hp_y_n3n2 = new TH1F("hp_y_n3n2", "dN/dp vs. p", 60, 0, 30);
0091
0092   TCanvas *cp_e_n3n2 = new TCanvas("cp_e_n3n2");
0093   t_truth->Draw("em_evtgen_ptotal>>hp_e_n3n2", electron_cut && eta_cut_n3n2 && "evtgen_Q2 > 0.01", "goff");
0094   t_truth->Draw("em_evtgen_ptotal>>hp_p_n3n2", Pion_charged_cut && eta_cut_n3n2 && "evtgen_Q2 > 0.01", "goff");
0095   t_truth->Draw("em_evtgen_ptotal>>hp_y_n3n2", photon_cut && eta_cut_n3n2 && status_cut && "evtgen_Q2 > 0.01", "goff");
0096
0097
0098   TH1F* htmp_n3n2 = hp_e_n3n2->Clone();
0099   htmp_n3n2->SetTitle("");
0100   htmp_n3n2->GetXaxis()->SetTitle("p [GeV]");
0101   htmp_n3n2->GetYaxis()->SetTitle("dN/dp");
0102   htmp_n3n2->SetMaximum( 0.99e7);
0103   htmp_n3n2->Draw();
0104
0105   hp_e_n3n2->SetLineColor(2);
0106   hp_e_n3n2->Draw("same");
0107   hp_p_n3n2->SetLineColor(1);
0108   hp_p_n3n2->Draw("same");
0109   hp_y_n3n2->SetLineColor(4);
0110   hp_y_n3n2->Draw("same");
0111   cp_e_n3n2->SetLogy();
0112
0113   TLegend* leg_n3n2 = new TLegend(0.53,0.67,0.73,0.90);
0114   hp_e_n3n2->SetTitle("DIS electron");
0115   hp_p_n3n2->SetTitle("#pi#pm");
0116   hp_y_n3n2->SetTitle("Photons");
0117   leg_n3n2->AddEntry(hp_e_n3n2, "", "L");
0118   leg_n3n2->AddEntry(hp_p_n3n2, "", "L");
0119   leg_n3n2->AddEntry(hp_y_n3n2, "", "L");
0120   leg_n3n2->SetTextSize(0.04);
0121   leg_n3n2->Draw();
0122
0123   if ( save_figures )
0124     {
0125       cp_e_n3n2->Print("plots/Pythia_MomSpec_n3n2.eps");
0126       cp_e_n3n2->Print("plots/Pythia_MomSpec_n3n2.png");
0127     }
0128
0129   /*Temporary Electron Purity Test*/
0130
0131
0132   /*-------- -2 < Eta < -1 ---------*/
0133
0134   TH1F* hp_e_n2n1 = new TH1F("hp_e_n2n1", "dN/dp vs. p", 60, 0, 30);
0135   TH1F* hp_p_n2n1 = new TH1F("hp_p_n2n1", "dN/dp vs. p", 60, 0, 30);
0136   TH1F* hp_y_n2n1 = new TH1F("hp_y_n2n1", "dN/dp vs. p", 60, 0, 30);
0137
0138   TCanvas *cp_e_n2n1 = new TCanvas("cp_e_n2n1");
0139   t_truth->Draw("em_evtgen_ptotal>>hp_e_n2n1", electron_cut && eta_cut_n2n1 && "evtgen_Q2 > 0.01", "goff");
0140   t_truth->Draw("em_evtgen_ptotal>>hp_p_n2n1", Pion_charged_cut && eta_cut_n2n1 && "evtgen_Q2 > 0.01", "goff");
0141   t_truth->Draw("em_evtgen_ptotal>>hp_y_n2n1", photon_cut && eta_cut_n2n1 && status_cut && "evtgen_Q2 > 0.01", "goff");
0142
0143   TH1F* htmp_n2n1 = hp_e_n2n1->Clone();
0144   htmp_n2n1->SetTitle("");
0145   htmp_n2n1->GetXaxis()->SetTitle("p [GeV]");
0146   htmp_n2n1->GetYaxis()->SetTitle("dN/dp");
0147   htmp_n2n1->SetMaximum( 0.99e7);
0148   htmp_n2n1->Draw();
0149
0150   hp_e_n2n1->SetLineColor(2);
0151   hp_e_n2n1->Draw("same");
0152   hp_p_n2n1->SetLineColor(1);
0153   hp_p_n2n1->Draw("same");
0154   hp_y_n2n1->SetLineColor(4);
0155   hp_y_n2n1->Draw("same");
0156   cp_e_n2n1->SetLogy();
0157
0158   TLegend* leg_n2n1 = new TLegend(0.53,0.67,0.73,0.90);
0159   hp_e_n2n1->SetTitle("DIS electron");
0160   hp_p_n2n1->SetTitle("#pi#pm");
0161   hp_y_n2n1->SetTitle("Photons");
0162   leg_n2n1->AddEntry(hp_e_n2n1, "", "L");
0163   leg_n2n1->AddEntry(hp_p_n2n1, "", "L");
0164   leg_n2n1->AddEntry(hp_y_n2n1, "", "L");
0165   leg_n2n1->SetTextSize(0.04);
0166   leg_n2n1->Draw();
0167
0168   if ( save_figures )
0169     {
0170       cp_e_n2n1->Print("plots/Pythia_MomSpec_n2n1.eps");
0171       cp_e_n2n1->Print("plots/Pythia_MomSpec_n2n1.png");
0172     }
0173
0174   /*-------- -1 < Eta < 0 ---------*/
0175
0176   TH1F* hp_e_n1z0 = new TH1F("hp_e_n1z0", "dN/dp vs. p", 60, 0, 30);
0177   TH1F* hp_p_n1z0 = new TH1F("hp_p_n1z0", "dN/dp vs. p", 60, 0, 30);
0178   TH1F* hp_y_n1z0 = new TH1F("hp_y_n1z0", "dN/dp vs. p", 60, 0, 30);
0179
0180   TCanvas *cp_e_n1z0 = new TCanvas("cp_e_n1z0");
0181   t_truth->Draw("em_evtgen_ptotal>>hp_e_n1z0", electron_cut && eta_cut_n1z0 && "evtgen_Q2 > 0.01", "goff");
0182   t_truth->Draw("em_evtgen_ptotal>>hp_p_n1z0", Pion_charged_cut && eta_cut_n1z0 && "evtgen_Q2 > 0.01", "goff");
0183   t_truth->Draw("em_evtgen_ptotal>>hp_y_n1z0", photon_cut && eta_cut_n1z0 && status_cut && "evtgen_Q2 > 0.01", "goff");
0184
0185   TH1F* htmp_n1z0 = hp_e_n1z0->Clone();
0186   htmp_n1z0->SetTitle("");
0187   htmp_n1z0->GetXaxis()->SetTitle("p [GeV]");
0188   htmp_n1z0->GetYaxis()->SetTitle("dN/dp");
0189   htmp_n1z0->SetMaximum( 0.99e7 );
0190   htmp_n1z0->Draw();
0191
0192   hp_e_n1z0->SetLineColor(2);
0193   hp_e_n1z0->Draw("same");
0194   hp_p_n1z0->SetLineColor(1);
0195   hp_p_n1z0->Draw("same");
0196   hp_y_n1z0->SetLineColor(4);
0197   hp_y_n1z0->Draw("same");
0198   cp_e_n1z0->SetLogy();
0199
0200   TLegend* leg_n1z0 = new TLegend(0.53,0.67,0.73,0.90);
0201   hp_e_n1z0->SetTitle("DIS electron");
0202   hp_p_n1z0->SetTitle("#pi#pm");
0203   hp_y_n1z0->SetTitle("Photons");
0204   leg_n1z0->AddEntry(hp_e_n1z0, "", "L");
0205   leg_n1z0->AddEntry(hp_p_n1z0, "", "L");
0206   leg_n1z0->AddEntry(hp_y_n1z0, "", "L");
0207   leg_n1z0->SetTextSize(0.04);
0208   leg_n1z0->Draw();
0209
0210   if ( save_figures )
0211     {
0212       cp_e_n1z0->Print("plots/Pythia_MomSpec_n1z0.eps");
0213       cp_e_n1z0->Print("plots/Pythia_MomSpec_n1z0.png");
0214     }
0215
0216   return 0 ;
0217 }