analysis/FrawleyTrackEvalMacros/cluster_efficiency.C

0001 #include "TChain.h"
0002 #include "TString.h"
0003 #include "TH1D.h"
0004 #include "TH2D.h"
0005 #include "TCanvas.h"
0006 #include "TROOT.h"
0007 #include "TStyle.h"
0008 #include "TLine.h"
0009 #include "TProfile.h"
0010 #include "TCut.h"
0011 #include "TLegend.h"
0012 #include "TLatex.h"
0013
0014 #include <iostream>
0015
0016 void cluster_efficiency()
0017 {
0018   gROOT->SetStyle("Plain");
0019   gStyle->SetOptStat(0);
0020   gStyle->SetOptTitle(1);
0021
0022   TCut good_g4cluster_cut = Form("layer >= 0 && gr > 0 && gembed > 0 && gtrackID > 0");
0023   TCut good_cluster_cut = Form("layer >= 0 && r > 0 && gembed > 0 && gtrackID > 0");
0024   TCut layer_cut = Form("layer > 6");
0025   TCut layer_cut_tpc[3];
0026   layer_cut_tpc[0] = Form("layer > 6 && layer <= 22");
0027   layer_cut_tpc[1] = Form("layer > 22 && layer <= 38");
0028   layer_cut_tpc[2] = Form("layer > 38 && layer <= 54");
0029   TCut layer_cut_intt = Form("layer > 2 && layer < 7");
0030   TCut layer_cut_mvtx = Form("layer < 3");
0031   TCut etacut = Form("geta > -1.0");
0032   //TCut etacut = Form("geta < 0.2 && geta > 0.1");
0033   //TCut etacut = Form("geta < 1.0 && geta > -1.0");
0034   //TCut ng4hitscut = Form("ng4hits > 0");
0035   TCut ng4hitscut = Form("");
0036
0037   TChain* ntp_g4cluster = new TChain("ntp_g4cluster","truth clusters");
0038   TChain* ntp_cluster = new TChain("ntp_cluster","reco clusters");
0039
0040   for(int ifile = 0; ifile < 100; ++ifile)
0041     {
0042       char name[500];
0043
0044       //sprintf(name,"/sphenix/user/frawley/tpc_cluster_qa/macros/macros/g4simulations/repo_100_tracks_occ_eval_output/g4svtx_eval_%i.root_g4svtx_eval.root",ifile);
0045       //sprintf(name,"/sphenix/user/frawley/tpc_cluster_qa/macros/macros/g4simulations/2k_tracks_occ_eval_output/g4svtx_eval_%i.root_g4svtx_eval.root",ifile);
0046       //sprintf(name,"/sphenix/user/frawley/tpc_cluster_qa/macros/macros/g4simulations/4k_tracks_occ_eval_output/g4svtx_eval_%i.root_g4svtx_eval.root",ifile);
0047       //sprintf(name,"/sphenix/user/frawley/tpc_cluster_qa/macros/macros/g4simulations/100muon_eval_output/g4svtx_eval_%i.root_g4svtx_eval.root",ifile);
0048
0049       sprintf(name,"/sphenix/user/frawley/tpc_cluster_qa/macros/macros/g4simulations/100_tracks_occ_eval_output/g4svtx_eval_%i.root_g4svtx_eval.root",ifile);
0050       //sprintf(name,"/sphenix/user/frawley/current_repo/macros/macros/g4simulations/genfit_pion_eval_output/g4svtx_eval_%i.root_g4svtx_eval.root",ifile);
0051       //sprintf(name,"/sphenix/user/frawley/tpc_cluster_qa/macros/macros/g4simulations/eval_output/g4svtx_eval_%i.root_g4svtx_eval.root",ifile);
0052
0053       //sprintf(name,"/sphenix/user/frawley/cluster_efficiency/macros/macros/g4simulations/low_occupancy_fiducial0.5_search2_eval_output/g4svtx_eval_%i.root_g4svtx_eval.root",ifile);
0054
0055       std::cout << "Open file " << name << std::endl;
0056
0057       ntp_g4cluster->Add(name);
0058       ntp_cluster->Add(name);
0059     }
0060
0061   // cluster efficiency vs layer
0062
0063   TCanvas * ceff = new TCanvas("ceff","ceff",5,5,1200,800);
0064   TH1D* h3_den = new TH1D("h3_den","; layer; Efficiency",54, 0, 54);
0065   TH1D* h3_num = (TH1D*)h3_den->Clone("h3_num");;
0066   TH1D* h3_eff = (TH1D*)h3_den->Clone("h3_eff");;
0067
0068   ntp_g4cluster->Draw("layer>>h3_den",good_g4cluster_cut);
0069   ntp_g4cluster->Draw("layer>>h3_num",good_cluster_cut);
0070
0071   for(int i=1;i<=h3_den->GetNbinsX();++i){
0072     double pass = h3_num->GetBinContent(i);
0073     double all = h3_den->GetBinContent(i);
0074     double eff = 0;
0075     if(all > pass)
0076       eff = pass/all;
0077     else if(all > 0)
0078       eff = 1.;
0079     //std::cout << " i " << i << " pass " << pass << " all " << all << " eff " << eff << std::endl;
0080     //double err = BinomialError(pass, all);
0081     h3_eff->SetBinContent(i, eff);
0082     //h3_eff->SetBinError(i, err);
0083   }
0084
0085   //h3_eff->Draw("e,text");
0086   h3_eff->Draw("p");
0087   h3_eff->SetStats(0);
0088   h3_eff->SetTitle("truth clusters with matched reco cluster; layer; cluster efficiency");
0089   h3_eff->SetMarkerStyle(20);
0090   h3_eff->SetMarkerColor(kBlack);
0091   h3_eff->SetLineColor(kBlack);
0092   //h3_eff->GetYaxis()->SetRangeUser(0.0, 1.1);
0093   h3_eff->GetYaxis()->SetRangeUser(0.0, 1.04);
0094
0095   TLine *unit = new TLine(0,1.0,54,1.0);
0096   unit->SetLineStyle(2);
0097   unit->Draw();
0098
0099   TCanvas *crestpc = new TCanvas("crestpc","crestpc",5,5,1600,800);
0100   crestpc->Divide(3,1);
0101   TH1D *hreslayer[3];
0102   for(int i=0;i<3;++i)
0103     {
0104       crestpc->cd(i+1);
0105       gPad->SetLeftMargin(0.15);
0106       char name[500];
0107       sprintf(name,"hreslayer%i",i);
0108       if(i==0)
0109     hreslayer[i] = new TH1D(name,"cluster r-phi resolution TPC inner",300,-0.095, 0.095);
0110       else if(i==1)
0111     hreslayer[i] = new TH1D(name,"cluster r-phi resolution TPC mid",300,-0.095, 0.095);
0112       else
0113       hreslayer[i] = new TH1D(name,"Cluster r-phi resolution TPC outer",300,-0.095, 0.095);
0114       char drawit[500];
0115       sprintf(drawit,"gr*(phi-gphi)>>hreslayer%i",i);
0116       ntp_cluster->Draw(drawit, good_cluster_cut && layer_cut_tpc[i]);
0117       hreslayer[i]->GetXaxis()->SetTitle("gr*(phi-gphi) (cm)");
0118       hreslayer[i]->GetXaxis()->SetNdivisions(504);
0119       hreslayer[i]->Draw("l");
0120
0121       cout << " RMS for TPC r-phi for i = " << i << "  is " << hreslayer[i]->GetRMS() << endl;;
0122
0123       char latstr[500];
0124       sprintf(latstr,"RMS = %.0f (#mum)",  1.0e04 * hreslayer[i]->GetRMS());
0125       TLatex *lat = new TLatex(0.6, 0.915, latstr);
0126       lat->SetNDC();
0127       lat->Draw();
0128     }
0129
0130   TCanvas *crestpcz = new TCanvas("crestpcz","crestpcz",5,5,1600,800);
0131   crestpcz->Divide(3,1);
0132   TH1D *hreslayerz[3];
0133   for(int i=0;i<3;++i)
0134     {
0135       crestpcz->cd(i+1);
0136       gPad->SetLeftMargin(0.15);
0137       char name[500];
0138       sprintf(name,"hreslayerz%i",i);
0139       if(i==0)
0140     hreslayerz[i] = new TH1D(name,"Cluster z resolution TPC inner",300,-0.4, 0.4);
0141       else if(i==1)
0142     hreslayerz[i] = new TH1D(name,"Cluster z resolution TPC mid",300,-0.4, 0.4);
0143       else
0144     hreslayerz[i] = new TH1D(name,"Cluster z resolution TPC outer",300,-0.4, 0.4);
0145       char drawit[500];
0146       sprintf(drawit,"z-gz>>hreslayerz%i",i);
0147       ntp_g4cluster->Draw(drawit, good_cluster_cut && layer_cut_tpc[i] && etacut && ng4hitscut);
0148       hreslayerz[i]->GetXaxis()->SetTitle("z-gz (cm)");
0149       hreslayerz[i]->GetXaxis()->SetNdivisions(504);
0150       hreslayerz[i]->Draw("l");
0151
0152       char latstr[500];
0153       sprintf(latstr,"RMS = %.0f (#mum)",  1.0e04 * hreslayerz[i]->GetRMS());
0154       TLatex *lat = new TLatex(0.6, 0.915, latstr);
0155       lat->SetNDC();
0156       lat->Draw();
0157     }
0158
0159   TCanvas *cressil = new TCanvas("cressil","cressil",5,5,1600,800);
0160   cressil->Divide(2,1);
0161   TH1D *hreslayer_sil[2];
0162   for(int i=0;i<2;++i)
0163     {
0164       cressil->cd(i+1);
0165       gPad->SetLeftMargin(0.15);
0166       char name[500];
0167       sprintf(name,"hreslayer_sil%i",i);
0168       if(i==0)
0169       hreslayer_sil[i] = new TH1D(name,"Cluster resolution INTT",300,-0.0049, 0.0049);
0170       else
0171       hreslayer_sil[i] = new TH1D(name,"Cluster resolution MVTX",300,-0.0025, 0.0025);
0172       char drawit[500];
0173       sprintf(drawit,"gr*(phi-gphi)>>hreslayer_sil%i",i);
0174       if(i==0)
0175     ntp_cluster->Draw(drawit, good_cluster_cut && layer_cut_intt);
0176       else
0177     ntp_cluster->Draw(drawit, good_cluster_cut && layer_cut_mvtx);
0178       hreslayer_sil[i]->GetXaxis()->SetTitle("gr*(phi-gphi) (cm)");
0179       hreslayer_sil[i]->GetXaxis()->SetNdivisions(504);
0180       hreslayer_sil[i]->Draw("l");
0181
0182       char latstr[500];
0183       sprintf(latstr,"RMS = %.1f (#mum)",  1.0e04 * hreslayer_sil[i]->GetRMS());
0184       TLatex *lat = new TLatex(0.6, 0.915, latstr);
0185       lat->SetNDC();
0186       lat->Draw();
0187     }
0188
0189   TCanvas *cressilz = new TCanvas("cressilz","cressilz",5,5,1600,800);
0190   cressilz->Divide(2,1);
0191   TH1D *hreslayer_silz[2];
0192   for(int i=0;i<2;++i)
0193     {
0194       cressilz->cd(i+1);
0195       gPad->SetLeftMargin(0.15);
0196       char name[500];
0197       sprintf(name,"hreslayer_silz%i",i);
0198       if(i==0)
0199     hreslayer_silz[i] = new TH1D(name,"Cluster resolution INTT z",300,-1.0, 1.0);
0200       else
0201     hreslayer_silz[i] = new TH1D(name,"Cluster resolution MVTX z",300,-0.0025, 0.0025);
0202       char drawit[500];
0203       sprintf(drawit,"z-gz>>hreslayer_silz%i",i);
0204       if(i==0)
0205     ntp_cluster->Draw(drawit, good_cluster_cut && layer_cut_intt);
0206       else
0207     ntp_cluster->Draw(drawit, good_cluster_cut && layer_cut_mvtx);
0208       hreslayer_silz[i]->GetXaxis()->SetTitle("z-gz (cm)");
0209       hreslayer_silz[i]->GetXaxis()->SetNdivisions(504);
0210       hreslayer_silz[i]->Draw("l");
0211
0212       char latstr[500];
0213       if(i==0)
0214     sprintf(latstr,"RMS = %.0f (#mum)",  1.0e04 * hreslayer_silz[i]->GetRMS());
0215       else
0216     sprintf(latstr,"RMS = %.1f (#mum)",  1.0e04 * hreslayer_silz[i]->GetRMS());
0217       TLatex *lat = new TLatex(0.6, 0.915, latstr);
0218       lat->SetNDC();
0219       lat->Draw();
0220     }
0221
0222   TCanvas *ctruth = new TCanvas("ctruth","ctruth",5,5,1200,800);
0223   ctruth->SetLeftMargin(0.15);
0224   TH1D *htruth = new TH1D("htruth","Truth clusters per layer",540,0,60);
0225   ntp_g4cluster->Draw("layer>>htruth", good_g4cluster_cut);
0226   htruth->SetMarkerStyle(20);
0227   htruth->SetMarkerSize(2);
0228   htruth->GetXaxis()->SetTitle("layer");
0229   htruth->GetYaxis()->SetTitle("truth clusters / layer");
0230   htruth->GetYaxis()->SetTitleOffset(1.7);
0231   htruth->Draw("p");
0232
0233   TCanvas *cphi = new TCanvas("cphi","cphi", 5,5,800,800);
0234   TH1D *hphi = new TH1D("hphi","hphi",10000, -4.0, 4.0);
0235   ntp_g4cluster->Draw("phi>>hphi", good_cluster_cut && layer_cut);
0236   hphi->Draw();
0237
0238   TCanvas *ceta = new TCanvas("ceta","ceta", 5,5,800,800);
0239   TH1D *heta = new TH1D("heta","heta",10000, -1.2, 1.2);
0240   ntp_g4cluster->Draw("eta>>heta", good_cluster_cut && layer_cut);
0241   heta->Draw();
0242
0243   TCanvas * cetaz = new TCanvas("cetaz","cetaz",5,5,800,800);
0244   TH2D *hetaz = new TH2D("hetaz","hetaz",200, -1.2, 1.2, 200, -0.5, 0.5);
0245   ntp_g4cluster->Draw("(z-gz):eta>>hetaz", good_cluster_cut && layer_cut);
0246   hetaz->Draw("colz");
0247
0248   /*
0249   TCanvas * cetazsize = new TCanvas("cetazsize","cetazsize",5,5,800,800);
0250   TH2D *hetazsize = new TH2D("hetazsize","hetazsize",200, -1.2, 1.2, 200, 0, 5.0);
0251   ntp_g4cluster->Draw("zsize:eta>>hetazsize", good_cluster_cut && layer_cut);
0252   hetazsize->Draw("colz");
0253   */
0254
0255
0256   // zsize
0257
0258   TCanvas *c2 = new TCanvas("c2","c2", 5,5,800,800);
0259
0260   //TH2D *h21 = new TH2D("h21","h21",800, -110, 110, 5000, 0, 8);
0261   TH2D *h21 = new TH2D("h21","h21",800, -1.2, 1.2, 5000, 0, 8);
0262   //ntp_g4cluster->Draw("zsize*.98:z>>h21",good_cluster_cut);
0263   ntp_g4cluster->Draw("zsize*.98:geta>>h21",good_cluster_cut);
0264   h21->SetMarkerStyle(20);
0265   h21->SetMarkerSize(0.5);
0266   //h21->GetXaxis()->SetTitle("Z (cm)");
0267   h21->GetXaxis()->SetTitle("#eta");
0268   h21->GetYaxis()->SetTitle("zsize (cm)");
0269
0270   //TH2D *h22 = new TH2D("h22","h22",800, -110, 110, 5000, 0, 8);
0271   TH2D *h22 = new TH2D("h22","h22",800, -1.2, 1.2, 5000, 0, 8);
0272   //ntp_g4cluster->Draw("gzsize:gz>>h22",good_g4cluster_cut);
0273   ntp_g4cluster->Draw("gzsize:geta>>h22",good_g4cluster_cut);
0274   h22->SetMarkerStyle(20);
0275   h22->SetMarkerSize(0.5);
0276   h22->SetMarkerColor(kRed);
0277   //h21->GetXaxis()->SetTitle("Z (cm)");
0278   h21->GetXaxis()->SetTitle("#eta");
0279   h21->GetYaxis()->SetTitle("zsize (cm)");
0280
0281   h21->Draw();
0282   h22->Draw("same");
0283
0284   // phisize
0285
0286   TCanvas *c4 = new TCanvas("c4","c4", 5,5,800,800);
0287
0288   TH2D *h41 = new TH2D("h41","h41",500, 0, 80.0, 5000, 0, 2.);
0289   ntp_g4cluster->Draw("phisize*0.98:r>>h41",good_cluster_cut);
0290
0291   h41->SetMarkerStyle(20);
0292   h41->SetMarkerSize(0.5);
0293   h41->GetXaxis()->SetTitle("radius (cm)");
0294   h41->GetYaxis()->SetTitle("radius*phisize (cm)");
0295
0296   TH2D *h42 = new TH2D("h42","h42",500, 0, 80.0, 5000, 0, 2.);
0297   ntp_g4cluster->Draw("gphisize:gr>>h42",good_g4cluster_cut);
0298   h42->SetMarkerStyle(20);
0299   h42->SetMarkerSize(0.5);
0300   h42->SetMarkerColor(kRed);
0301   h42->GetXaxis()->SetTitle("radius (cm)");
0302   h42->GetYaxis()->SetTitle("radius*phisize (cm)");
0303
0304
0305   h41->Draw();
0306   h42->Draw("same");
0307
0308
0309   //h42->Draw();
0310   //h41->Draw("same");
0311
0312 }