analysis/FrawleyTrackEvalMacros/plot_ntp_g4hit.C

0001 #include "TCanvas.h"
0002 #include "TH1D.h"
0003 #include "TH2D.h"
0004 #include "TCut.h"
0005 #include "TChain.h"
0006 #include "TLegend.h"
0007 #include "TFile.h"
0008
0009 //#include "sPhenixStyle.C"
0010 #include "/phenix/u/yuhw/RootMacros/sPHENIXStyle/sPhenixStyle.C"
0011
0012 void plot_ntp_g4hit(
0013             //const char* input = "g4svtx_eval.root",
0014             const char* input = "/sphenix/user/frawley/pileup/macros/macros/g4simulations/eval_output/g4svtx_eval_1.root_g4svtx_eval.root",
0015             const int min_nfromtruth = 30
0016         ){
0017
0018   //SetsPhenixStyle();
0019   gROOT->SetStyle("Plain");
0020   gStyle->SetOptStat(0);
0021   gStyle->SetOptFit(0);
0022   gStyle->SetOptTitle(1);
0023
0024   bool acts = false; // is this output from acts tracking?
0025
0026   // OK for embedded particles only
0027   //TCut good_gtrack_cut = Form("gtrackID>=0 && gntpc > 20");
0028   //TCut good_track_cut = Form("gtrackID>=0 && ntpc > 20 && quality < 10");
0029
0030   // For Hijing events with embedded particles
0031   TCut good_gtrack_cut = Form("gtrackID>=0 && gembed >= 2 && gntpc > 20 && gnmaps > 1");
0032   TCut good_track_cut = Form("gtrackID>=0 && gembed >= 2 && ntpc > 20 && quality < 10 && nmaps > 1");
0033
0034   //TCut good_gtrack_cut = Form("gtrackID>=0 && gembed == 2 && gntpc > 20");
0035   //TCut good_track_cut = Form("gtrackID>=0 && gembed == 2 && ntpc > 20 && quality < 10");
0036
0037
0038   // rough 4 sigma cut
0039   //TCut pt_cut = Form("fabs((pt - gpt)/pt) < 4*(0.02+0.0012*gpt)");  // 4 times sigma (where sigma = 2.2% at 2 and 4.4% at 20
0040
0041   TFile *fout = new TFile("root_files/ntp_track_out.root","recreate");
0042
0043   TChain* ntp_track = new TChain("ntp_track","reco tracks");
0044   TChain* ntp_gtrack = new TChain("ntp_gtrack","g4 tracks");
0045
0046   bool use_list = false;
0047   int n_list = 1000;
0048
0049   int ifile = 0;
0050   int nbadvtx = 0;
0051   for(int i=0;i<n_list;i++)
0052     {
0053       char name[500];
0054       ifile = i;
0055
0056       //sprintf(name,"/sphenix/user/frawley/acts_qa/macros/macros/g4simulations/eval_output/g4svtx_eval_%i.root_g4svtx_eval.root",ifile);
0057       sprintf(name,"/sphenix/user/frawley/acts_qa/macros/macros/g4simulations/eval_output_2/g4svtx_eval_%i.root_g4svtx_eval.root",ifile);
0058       //sprintf(name,"/sphenix/user/frawley/acts_qa/macros/macros/g4simulations/eval_output_3/g4svtx_eval_%i.root_g4svtx_eval.root",ifile);
0059
0060       // Skip any files where the event vertex was not reconstructed properly
0061       TChain* ntp_vertex = new TChain("ntp_vertex","events");
0062       ntp_vertex->Add(name);
0063       TH1D *hzvtx = new TH1D("hzvtx","hzvtx",500,-10.0,10.0);
0064       ntp_vertex->Draw("vz - gvz >>hzvtx");
0065       double meanzdiff =  hzvtx->GetMean() ;
0066       cout << "   file " << i << " has vertex mean Z difference = " << meanzdiff << endl;
0067       if( fabs(meanzdiff) >  0.1 )
0068     {
0069       cout << "  --- Bad event vertex, skip file number " << i << " with name " << name << endl;
0070       nbadvtx++;
0071       delete ntp_vertex;
0072       delete hzvtx;
0073       continue;
0074     }
0075       delete ntp_vertex;
0076       delete hzvtx;
0077
0078       cout << "Adding file number " << i << " with name " << name << endl;
0079
0080       ntp_gtrack->Add(name);
0081       ntp_track->Add(name);
0082     }
0083   cout << " Found " << nbadvtx << " events with badly reconstructed vertices, and skipped them" << endl << endl;
0084
0085   cout << " ntp_gtrack chain entries = " << ntp_gtrack->GetEntries() << endl;
0086   cout << " ntp_track chain entries = " << ntp_track->GetEntries() << endl;
0087
0088   /*
0089   TH1D *hch = new TH1D("hch","hch",100,-2,+2);
0090   TH1D *hgch = new TH1D("hgch","hgch",100,-2,+2);
0091   ntp_track->Draw("charge>>hch",good_track_cut);
0092   ntp_gtrack->Draw("charge>>hgch",good_gtrack_cut);
0093   */
0094
0095   TCanvas *cmvtx = new TCanvas("cmvtx","cmvtx");
0096   TH1D *hnmvtx = new TH1D("hnmvtx","hnmvtx", 100, -1, 4);
0097   ntp_track->Draw("nmaps>>hnmvtx", good_track_cut);
0098   hnmvtx->Draw();
0099   hnmvtx->Write();
0100
0101   int nbinptres = 80;  // number of bins in pT
0102   int nbinpteff = 120;  // number of bins in pT
0103   int nbinptdca = 40;  // number of bins in pT
0104   int nbin = 400;  // number of bins in Y slice
0105   double range = 0.4;
0106   double ptmax = 40.0;
0107
0108   TCanvas *cpt = new TCanvas("cpt","cpt");
0109   TH2D *h1 = new TH2D("h1","h1",nbinptres, 0.0 ,ptmax,nbin,-range,range);
0110   ntp_track->Draw("(pt-gpt)/gpt:gpt>>h1",good_track_cut);
0111   //h1->FitSlicesY(0,0,-1,0,"qnrl");
0112   h1->FitSlicesY();
0113   TH1D*h1_1 = (TH1D*)gDirectory->Get("h1_1");
0114   TH1D*h1_2 = (TH1D*)gDirectory->Get("h1_2");
0115   h1_2->Draw("e");
0116   //h1_1->Draw("same");
0117
0118   h1_1->SetMarkerStyle(4);
0119   h1_1->SetMarkerColor(kBlack);
0120   h1_1->SetLineColor(kBlack);
0121
0122   h1_2->SetMarkerStyle(20);
0123   h1_2->SetMarkerColor(kBlack);
0124   h1_2->SetLineColor(kBlack);
0125   h1_2->GetYaxis()->SetRangeUser(0, 0.10);
0126   h1_2->GetYaxis()->SetTitleOffset(1.5);
0127   h1_2->SetStats(0);
0128   h1_2->SetTitle(";p_{T} [GeV/c];#frac{#Delta p_{T}}{p_{T}}");
0129
0130   h1->Write();
0131   h1_2->Write();
0132
0133   TCanvas * c3 = new TCanvas("c3","c3");
0134   TH1D* h3_den = new TH1D("h3_den","; p_{T}; Efficiency",nbinpteff, 0, ptmax);
0135   TH1D* h3_num = (TH1D*)h3_den->Clone("h3_num");;
0136   TH1D* h3_eff = (TH1D*)h3_den->Clone("h3_eff");;
0137
0138   cout<<__LINE__<<": "<< good_track_cut <<endl;
0139   ntp_gtrack->Draw("gpt>>h3_den",good_gtrack_cut);
0140   //ntp_track->Draw("gpt>>h3_num",good_track_cut && pt_cut);
0141   ntp_track->Draw("gpt>>h3_num",good_track_cut);
0142
0143   for(int i=1;i<=h3_den->GetNbinsX();++i){
0144     double pass = h3_num->GetBinContent(i);
0145     double all = h3_den->GetBinContent(i);
0146     double eff = 0;
0147     if(all > pass)
0148       eff = pass/all;
0149     else if(all > 0)
0150       eff = 1.;
0151
0152     //double err = BinomialError(pass, all);
0153     h3_eff->SetBinContent(i, eff);
0154     //h3_eff->SetBinError(i, err);
0155   }
0156
0157   //h3_eff->Draw("e,text");
0158   h3_eff->Draw("p");
0159   h3_eff->SetStats(0);
0160   h3_eff->SetTitle("; p_{T} [GeV/c]; eff.");
0161   h3_eff->SetMarkerStyle(20);
0162   h3_eff->SetMarkerColor(kBlack);
0163   h3_eff->SetLineColor(kBlack);
0164   h3_eff->GetYaxis()->SetRangeUser(0.0, 1.1);
0165
0166   h3_den->Write();
0167   h3_num->Write();
0168   h3_eff->Write();
0169
0170   TCanvas *c4 = new TCanvas("c4","c4",5,5,800,800);
0171
0172
0173   TH2D *h2 = new TH2D("h2","h2",nbinptdca, 0,ptmax,nbin,-0.06,0.06);
0174   ntp_track->Draw("(dca3dxy):gpt>>h2",good_track_cut);
0175   h2->FitSlicesY();
0176   TH1D*h2_1 = (TH1D*)gDirectory->Get("h2_1");
0177   TH1D*h2_2 = (TH1D*)gDirectory->Get("h2_2");
0178   h2_2->Draw("e");
0179
0180   h2_1->SetMarkerStyle(4);
0181   h2_1->SetMarkerColor(kBlack);
0182   h2_1->SetLineColor(kBlack);
0183
0184   h2_2->SetMarkerStyle(20);
0185   h2_2->SetMarkerColor(kBlack);
0186   h2_2->SetLineColor(kBlack);
0187   //h2_2->GetYaxis()->SetRangeUser(0, 0.1);
0188   h2_2->GetYaxis()->SetRangeUser(0.,0.01);
0189   h2_2->GetYaxis()->SetTitleOffset(1.5);
0190   h2_2->SetStats(0);
0191   h2_2->SetTitle(";p_{T} [GeV/c];dca2d (cm)}");
0192
0193   h2->Write();
0194   h2_2->Write();
0195
0196
0197   TCanvas *c5 = new TCanvas("c5","c5",5,5,800,800);
0198
0199   TH2D *h3 = new TH2D("h3","h3",nbinptdca, 0, ptmax, nbin, -0.1, 0.1);
0200   if(acts)
0201     ntp_track->Draw("(dca3dz-gvz):gpt>>h3",good_track_cut);
0202   else
0203     ntp_track->Draw("(dca3dz):gpt>>h3",good_track_cut);
0204   h3->FitSlicesY();
0205   TH1D*h3_1 = (TH1D*)gDirectory->Get("h3_1");
0206   TH1D*h3_2 = (TH1D*)gDirectory->Get("h3_2");
0207   h3_2->Draw("e");
0208   //h3_1->Draw("same");
0209
0210   h3_1->SetMarkerStyle(4);
0211   h3_1->SetMarkerColor(kBlack);
0212   h3_1->SetLineColor(kBlack);
0213
0214   h3_2->SetMarkerStyle(20);
0215   h3_2->SetMarkerColor(kBlack);
0216   h3_2->SetLineColor(kBlack);
0217   //h3_2->GetYaxis()->SetRangeUser(0, 0.1);
0218   h3_2->GetYaxis()->SetRangeUser(0.,0.01);
0219   h3_2->GetYaxis()->SetTitleOffset(1.5);
0220   h3_2->SetStats(0);
0221   h3_2->SetTitle(";p_{T} [GeV/c];dca2d (cm)}");
0222
0223   h3->Write();
0224   h3_2->Write();
0225
0226   TCanvas *c6 = new TCanvas("c6","c6", 5,5,1200, 800);
0227
0228   TH2D *h6 = new TH2D("h6","h6",nbinptdca, 0, ptmax, nbin, -1.0, 4.0);
0229   ntp_track->Draw("ntrumaps:gpt>>h6",good_track_cut);
0230   h6->Draw();
0231
0232   TH1D *h6_1 = new TH1D("h6_1","h6_1",nbin, 0.0, 1.2);
0233   int binlo = h6->GetYaxis()->FindBin(2.5);
0234   int binhi = h6->GetYaxis()->FindBin(3.5);
0235   h6->ProjectionX("h6_1",binlo,binhi);
0236   h6_1->GetXaxis()->SetTitle("p_{T} (GeV/c)");
0237   h6_1->SetLineWidth(2.0);
0238   h6_1->SetLineColor(kBlue);
0239
0240
0241   TH1D *h6_2 = new TH1D("h6_2","h6_2",nbin, 0.0, 1.2);
0242   binlo = h6->GetYaxis()->FindBin(-0.5);
0243   binhi = h6->GetYaxis()->FindBin(0.5);
0244   h6->ProjectionX("h6_2",binlo,binhi);
0245   h6_2->GetXaxis()->SetTitle("p_{T} (GeV/c)");
0246   h6_2->SetLineColor(kRed);
0247   h6_2->SetLineWidth(2.0);
0248
0249
0250   h6_1->Draw();
0251   h6_2->Draw("same");
0252
0253   TLegend *l = new TLegend(0.15,0.75,0.45,0.85,"","NDC");
0254   l->SetBorderSize(1);
0255   l->AddEntry(h6_2,"ntrumaps = 0","l");
0256   l->AddEntry(h6_1,"ntrumaps = 3","l");
0257   l->Draw();
0258
0259   h6->Write();
0260   h6_1->Write();
0261   h6_2->Write();
0262
0263   /*
0264   TCanvas *c7 = new TCanvas("c7","c7", 5,5,1200, 800);
0265
0266   TH2D *h7 = new TH2D("h7","h7",nbin, 0, 1.1, nbin, -1.0, 50.0);
0267   ntp_track->Draw("ntrutpc:gpt>>h7","ntrumaps == 0");
0268   h7->Draw();
0269   */
0270   /*
0271   TH1D *h7_1 = new TH1D("h7_1","h7_1",100, 0.0, 1.2);
0272   h7->ProjectionX("h7_1");
0273   h7_1->GetXaxis()->SetTitle("p_{T} (GeV/c)");
0274   h7_1->SetLineWidth(2.0);
0275   h7_1->SetLineColor(kBlue);
0276   h7_1->Draw();
0277   */
0278   /*
0279   TH2D *h8 = new TH2D("h8","h8",100, 0, 1.1, 100, -1.0, 50.0);
0280   ntp_track->Draw("ntrutpc:gpt>>h7","ntrumaps == 3");
0281   h8->Draw();
0282
0283   TH1D *h8_1 = new TH1D("h8_1","h8_1",100, 0.0, 1.2);
0284   h8->ProjectionX("h8_1");
0285   h8_1->GetXaxis()->SetTitle("p_{T} (GeV/c)");
0286   h8_1->SetLineColor(kRed);
0287   h8_1->SetLineWidth(2.0);
0288   h8_1->Draw("same");
0289
0290   TLegend *l = new TLegend(0.15,0.75,0.45,0.85,"","NDC");
0291   l->SetBorderSize(1);
0292   l->AddEntry(h7_1,"ntrumaps = 0","l");
0293   l->AddEntry(h8_1,"ntrumaps = 3","l");
0294   l->Draw();
0295   */
0296 }