EICMCAnalysis/DIS-Pythia6/makePlot_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_LoI_Momentum(
0009                           const TString infile = "data/pythia.ep.100.test.root"
0010                           )
0011 {
0012   gStyle->SetOptStat(kFALSE);
0013
0014   /*--------------Get Input File--------------*/
0015
0016   TFile *f_pyth = new TFile( infile, "OPEN");
0017   TTree* T = (TTree*)f_pyth->Get("tree");
0018   T->Print();
0019
0020   /*------------------------------------------*/
0021   /*--------------Define Cuts-----------------*/
0022
0023   //TCut Q2_cut = "Q2 > 1";
0024
0025   TCut electron_cut = "p.fKF == 11";
0026   TCut hadron_cut = "abs(p.fKF) > 100";
0027   TCut proton_cut = "p.fKF == 2212";
0028   TCut neutron_cut = "p.fKF == 2112";
0029   TCut Kaon_cut = "abs(p.fKF) == 321 || p.fKF == 311";
0030   TCut Pion_charged_cut = "abs(p.fKF) == 211";
0031   //TCut Pion_cut = "abs(p.fKF) == 211 || p.fKF == 111";
0032   TCut photon_cut = "p.fKF == 22";
0033
0034   TCut mother_cut = "p.fParent == 0";
0035   TCut status_cut = "p.fKS < 10";
0036
0037   TCut eta_cut_n3n2 = "TMath::ASinH(p.fPz/sqrt((p.fPx)**2 + (p.fPy)**2)) > -3 && TMath::ASinH(p.fPz/sqrt((p.fPx)**2 + (p.fPy)**2)) < -2";
0038   TCut eta_cut_n2n1 = "TMath::ASinH(p.fPz/sqrt((p.fPx)**2 + (p.fPy)**2)) > -2 && TMath::ASinH(p.fPz/sqrt((p.fPx)**2 + (p.fPy)**2)) < -1";
0039   TCut eta_cut_n1z0 = "TMath::ASinH(p.fPz/sqrt((p.fPx)**2 + (p.fPy)**2)) > -1 && TMath::ASinH(p.fPz/sqrt((p.fPx)**2 + (p.fPy)**2)) < -0";
0040
0041   /*------------------------------------------*/
0042   /*-------Momentum vs. Pseudorapidity--------*/
0043   /*------------------------------------------*/
0044
0045   /*------------------------------------------*/
0046   /*---------Electrons (LoI Fig. 2-1)---------*/
0047   /*------------------------------------------*/
0048   TH2F *h_peta_e = new TH2F("h_peta_e", "", 100, -4, 3, 100, 0, 50 ); //250x015 10M
0049
0050   TCanvas *c_peta_e = new TCanvas( "c_peta_e" );
0051   c_peta_e->SetLogz();
0052
0053   T->Draw("sqrt((p.fPx)**2 + (p.fPy)**2 + (p.fPz)**2 ):TMath::ASinH(p.fPz/sqrt((p.fPx)**2 + (p.fPy)**2))>>h_peta_e", electron_cut && mother_cut && "Q2 > 1", "colz");
0054
0055   h_peta_e->GetXaxis()->SetTitle("Pseudorapidity #eta");
0056   h_peta_e->GetYaxis()->SetTitle("Electron Momentum p_{e-} [GeV]");
0057
0058   //  c_p_eta_e->Print("Plots/Pythia_peta_e_10M_250x010.eps");
0059   //  c_p_eta_e->Print("Plots/Pythia_peta_e_10M_250x010.png");
0060
0061   /*------------------------------------------*/
0062   /*----------Hadrons (LoI Fig. 2-4)----------*/
0063   /*------------------------------------------*/
0064   TH2F *h_p_eta_h = new TH2F("h_p_eta_h", "", 60,-6,6, 200,0,225); //250x015 10M
0065
0066   TCanvas *c_p_eta_h = new TCanvas( "c_p_eta_h" );
0067   // T->Draw("sqrt((p.fPx)**2 + (p.fPy)**2 + (p.fPz)**2 ):TMath::ASinH(p.fPz/sqrt((p.fPx)**2 + (p.fPy)**2))>>h_p_eta_h",hadron_cut && "Q2 > 1 && 0.01 < y < 0.80 && W2 > 10", "colz");
0068   T->Draw("sqrt((p.fPx)**2 + (p.fPy)**2 + (p.fPz)**2 ):TMath::ASinH(p.fPz/sqrt((p.fPx)**2 + (p.fPy)**2))>>h_p_eta_h",hadron_cut && "Q2 > 1 && 0.01 < y < 0.80", "colz");
0069
0070   h_p_eta_h->GetXaxis()->SetTitle("Pseudorapidity #eta");
0071   h_p_eta_h->GetYaxis()->SetTitle("Hadron Momentum p_{Hadron} [GeV]");
0072   c_p_eta_h->SetLogz();
0073
0074   //  c_p_eta_h->Print("Plots/Pythia_peta_h_10M_250x010.eps");
0075   //  c_p_eta_h->Print("Plots/Pythia_peta_h_10M_250x010.png");
0076
0077
0078   /*------------Momentum Spectra (Figure 2.2)--------------*/
0079
0080   /*-------- -3 < Eta < -2 ---------*/
0081
0082   TH1F* hp_e_n3n2 = new TH1F("hp_e_n3n2", "dN/dp vs. p", 60, 0, 30);
0083   TH1F* hp_p_n3n2 = new TH1F("hp_p_n3n2", "dN/dp vs. p", 60, 0, 30);
0084   TH1F* hp_y_n3n2 = new TH1F("hp_y_n3n2", "dN/dp vs. p", 60, 0, 30);
0085
0086   TCanvas *cp_e_n3n2 = new TCanvas("cp_e_n3n2");
0087   T->Draw("sqrt(p.fPx**2 + p.fPy**2 + p.fPz**2)>>hp_e_n3n2", electron_cut && eta_cut_n3n2 && "Q2 > 0.01", "goff");
0088   T->Draw("sqrt(p.fPx**2 + p.fPy**2 + p.fPz**2)>>hp_p_n3n2", Pion_charged_cut && eta_cut_n3n2 && "Q2 > 0.01", "goff");
0089   T->Draw("sqrt(p.fPx**2 + p.fPy**2 + p.fPz**2)>>hp_y_n3n2", photon_cut && eta_cut_n3n2 && status_cut && "Q2 > 0.01", "goff");
0090
0091
0092   TH1F* htmp_n3n2 = hp_e_n3n2->Clone();
0093   htmp_n3n2->SetTitle("");
0094   htmp_n3n2->GetXaxis()->SetTitle("p [GeV]");
0095   htmp_n3n2->GetYaxis()->SetTitle("dN/dp");
0096   htmp_n3n2->SetMaximum( 0.99e7);
0097   htmp_n3n2->Draw();
0098
0099   hp_e_n3n2->SetLineColor(2);
0100   hp_e_n3n2->Draw("same");
0101   hp_p_n3n2->SetLineColor(1);
0102   hp_p_n3n2->Draw("same");
0103   hp_y_n3n2->SetLineColor(4);
0104   hp_y_n3n2->Draw("same");
0105   cp_e_n3n2->SetLogy();
0106
0107   TLegend* leg_n3n2 = new TLegend(0.53,0.67,0.73,0.90);
0108   hp_e_n3n2->SetTitle("DIS electron");
0109   hp_p_n3n2->SetTitle("#pi#pm");
0110   hp_y_n3n2->SetTitle("Photons");
0111   leg_n3n2->AddEntry(hp_e_n3n2, "", "L");
0112   leg_n3n2->AddEntry(hp_p_n3n2, "", "L");
0113   leg_n3n2->AddEntry(hp_y_n3n2, "", "L");
0114   leg_n3n2->SetTextSize(0.04);
0115   leg_n3n2->Draw();
0116
0117   //  cp_e_n3n2->Print("Plots/Pythia_MomSpec_n3n2_10M_250x010.eps");
0118   //  cp_e_n3n2->Print("Plots/Pythia_MomSpec_n3n2_10M_250x010.png");
0119
0120
0121   /*Temporary Electron Purity Test*/
0122
0123
0124
0125
0126
0127
0128   /*-------- -2 < Eta < -1 ---------*/
0129
0130   TH1F* hp_e_n2n1 = new TH1F("hp_e_n2n1", "dN/dp vs. p", 60, 0, 30);
0131   TH1F* hp_p_n2n1 = new TH1F("hp_p_n2n1", "dN/dp vs. p", 60, 0, 30);
0132   TH1F* hp_y_n2n1 = new TH1F("hp_y_n2n1", "dN/dp vs. p", 60, 0, 30);
0133
0134   TCanvas *cp_e_n2n1 = new TCanvas("cp_e_n2n1");
0135   T->Draw("sqrt(p.fPx**2 + p.fPy**2 + p.fPz**2)>>hp_e_n2n1", electron_cut && eta_cut_n2n1 && "Q2 > 0.01", "goff");
0136   T->Draw("sqrt(p.fPx**2 + p.fPy**2 + p.fPz**2)>>hp_p_n2n1", Pion_charged_cut && eta_cut_n2n1 && "Q2 > 0.01", "goff");
0137   T->Draw("sqrt(p.fPx**2 + p.fPy**2 + p.fPz**2)>>hp_y_n2n1", photon_cut && eta_cut_n2n1 && status_cut && "Q2 > 0.01", "goff");
0138
0139   TH1F* htmp_n2n1 = hp_e_n2n1->Clone();
0140   htmp_n2n1->SetTitle("");
0141   htmp_n2n1->GetXaxis()->SetTitle("p [GeV]");
0142   htmp_n2n1->GetYaxis()->SetTitle("dN/dp");
0143   htmp_n2n1->SetMaximum( 0.99e7);
0144   htmp_n2n1->Draw();
0145
0146   hp_e_n2n1->SetLineColor(2);
0147   hp_e_n2n1->Draw("same");
0148   hp_p_n2n1->SetLineColor(1);
0149   hp_p_n2n1->Draw("same");
0150   hp_y_n2n1->SetLineColor(4);
0151   hp_y_n2n1->Draw("same");
0152   cp_e_n2n1->SetLogy();
0153
0154   TLegend* leg_n2n1 = new TLegend(0.53,0.67,0.73,0.90);
0155   hp_e_n2n1->SetTitle("DIS electron");
0156   hp_p_n2n1->SetTitle("#pi#pm");
0157   hp_y_n2n1->SetTitle("Photons");
0158   leg_n2n1->AddEntry(hp_e_n2n1, "", "L");
0159   leg_n2n1->AddEntry(hp_p_n2n1, "", "L");
0160   leg_n2n1->AddEntry(hp_y_n2n1, "", "L");
0161   leg_n2n1->SetTextSize(0.04);
0162   leg_n2n1->Draw();
0163
0164   //  cp_e_n2n1->Print("Plots/Pythia_MomSpec_n2n1_10M_250x010.eps");
0165   //  cp_e_n2n1->Print("Plots/Pythia_MomSpec_n2n1_10M_250x010.png");
0166
0167   /*-------- -1 < Eta < 0 ---------*/
0168
0169   TH1F* hp_e_n1z0 = new TH1F("hp_e_n1z0", "dN/dp vs. p", 60, 0, 30);
0170   TH1F* hp_p_n1z0 = new TH1F("hp_p_n1z0", "dN/dp vs. p", 60, 0, 30);
0171   TH1F* hp_y_n1z0 = new TH1F("hp_y_n1z0", "dN/dp vs. p", 60, 0, 30);
0172
0173   TCanvas *cp_e_n1z0 = new TCanvas("cp_e_n1z0");
0174   T->Draw("sqrt(p.fPx**2 + p.fPy**2 + p.fPz**2)>>hp_e_n1z0", electron_cut && eta_cut_n1z0 && "Q2 > 0.01", "goff");
0175   T->Draw("sqrt(p.fPx**2 + p.fPy**2 + p.fPz**2)>>hp_p_n1z0", Pion_charged_cut && eta_cut_n1z0 && "Q2 > 0.01", "goff");
0176   T->Draw("sqrt(p.fPx**2 + p.fPy**2 + p.fPz**2)>>hp_y_n1z0", photon_cut && eta_cut_n1z0 && status_cut && "Q2 > 0.01", "goff");
0177
0178   TH1F* htmp_n1z0 = hp_e_n1z0->Clone();
0179   htmp_n1z0->SetTitle("");
0180   htmp_n1z0->GetXaxis()->SetTitle("p [GeV]");
0181   htmp_n1z0->GetYaxis()->SetTitle("dN/dp");
0182   htmp_n1z0->SetMaximum( 0.99e7 );
0183   htmp_n1z0->Draw();
0184
0185   hp_e_n1z0->SetLineColor(2);
0186   hp_e_n1z0->Draw("same");
0187   hp_p_n1z0->SetLineColor(1);
0188   hp_p_n1z0->Draw("same");
0189   hp_y_n1z0->SetLineColor(4);
0190   hp_y_n1z0->Draw("same");
0191   cp_e_n1z0->SetLogy();
0192
0193   TLegend* leg_n1z0 = new TLegend(0.53,0.67,0.73,0.90);
0194   hp_e_n1z0->SetTitle("DIS electron");
0195   hp_p_n1z0->SetTitle("#pi#pm");
0196   hp_y_n1z0->SetTitle("Photons");
0197   leg_n1z0->AddEntry(hp_e_n1z0, "", "L");
0198   leg_n1z0->AddEntry(hp_p_n1z0, "", "L");
0199   leg_n1z0->AddEntry(hp_y_n1z0, "", "L");
0200   leg_n1z0->SetTextSize(0.04);
0201   leg_n1z0->Draw();
0202
0203   //  cp_e_n1z0->Print("Plots/Pythia_MomSpec_n1z0_10M_250x010.eps");
0204   //  cp_e_n1z0->Print("Plots/Pythia_MomSpec_n1z0_10M_250x010.png");
0205
0206   return 0 ;
0207 }