Back to home page

sPhenix code displayed by LXR

 
 

    


File indexing completed on 2025-08-03 08:13:11

0001 #include "TH2D.h"
0002 #include "TH1D.h"
0003 #include "TFile.h"
0004 #include "TCanvas.h"
0005 #include "TROOT.h"
0006 #include "TStyle.h"
0007 #include "TLine.h"
0008 #include "TLegend.h"
0009 #include "TF1.h"
0010 
0011 #include <iostream>
0012 
0013 void plot_ntp_track_out()
0014 {
0015   gROOT->SetStyle("Plain");
0016   gStyle->SetOptStat(0);
0017   gStyle->SetOptFit(1);
0018   gStyle->SetOptTitle(0);
0019 
0020   double ptmax = 40.0;
0021   double slice_low = 30.0;
0022   double slice_high = 31.0;
0023 
0024   std::vector<std::string> finvec;
0025   std::vector<std::string> legvec;
0026   std::vector<int> col;
0027 
0028   //============
0029 
0030 
0031   finvec.push_back("root_files/ntp_track_out_acts_MB_50khz_nomapscut.root");
0032   legvec.push_back( "TpcTracker+stub matcher+ActsTrkFit (MB+50, nomapscut)");
0033   col.push_back(kRed);
0034 
0035   finvec.push_back("root_files/ntp_track_out_acts_MB_50khz_withmapscut.root");
0036   legvec.push_back( "TpcTracker+stub matcher+ActsTrkFit (MB+50, nmaps>1)");
0037   col.push_back(kBlue);
0038 
0039 
0040 
0041   /*
0042   finvec.push_back("root_files/ntp_track_out_lo_hough_genfitprop_genfit.root");
0043   legvec.push_back( "Hough+GenfitTrkProp+GenfitTrkFitter (LO)");
0044   col.push_back(kBlue);
0045   */
0046 
0047 
0048 
0049   /*
0050   finvec.push_back("root_files/ntp_track_out_stub_matcher_lo_4kevts.root");
0051   legvec.push_back( "TpcTracker+stub matcher+ActsTrkFit (lo)");
0052   col.push_back(kRed);
0053 
0054   finvec.push_back("root_files/ntp_track_out_truth_assoc_lo_4kevts.root");
0055   legvec.push_back( "TpcTracker+truthassoc+ActsTrkFit (lo)");
0056   col.push_back(kBlack);
0057   */
0058 
0059   /*
0060   //  repeat after merge with repo of finvec.push_back("root_files/ntp_track_out_4k_tune8_fixed_mag24_both_stub_matcher.root");
0061   finvec.push_back("root_files/ntp_track_out.root");
0062   legvec.push_back( "TpcTracker+stub matcher+ActsTrkFit (0.02,0.015) (4K)");
0063   col.push_back(kBlack);
0064   */
0065 
0066   /*
0067   finvec.push_back("root_files/ntp_track_out_4k_tune8_fixed_mag24_both_stub_matcher.root");
0068   legvec.push_back( "TpcTracker+stub matcher+ActsTrkFit (0.02,0.015) (4K)");
0069   col.push_back(kBlack);
0070   */
0071 
0072   /*  
0073   finvec.push_back("root_files/ntp_track_out_4k_tune7_fixed_mag24_both_stub_matcher.root");
0074   legvec.push_back( "TpcTracker+stub matcher+ActsTrkFit (0.02,0.02) (4K)");
0075   col.push_back(kMagenta);
0076   */
0077 
0078   /*
0079   finvec.push_back("root_files/ntp_track_out_4k_tune1_fixed_mag24_both_stub_matcher.root");
0080   legvec.push_back( "TpcTracker+stub matcher+ActsTrkFit (0.03,0.02) (4K)");
0081   col.push_back(kBlack);
0082   */
0083 
0084   /*  
0085   finvec.push_back("root_files/ntp_track_out_4k_truth_assoc.root");
0086   legvec.push_back( "TPCTracker+TruthSiliconAssoc+ActsTrkFit (4K)");
0087   col.push_back(kRed);
0088   */
0089 
0090   /*
0091   finvec.push_back("root_files/ntp_track_out_4k_tune5_fixed_mag24_both_stub_matcher.root");
0092   legvec.push_back( "TpcTracker+stub matcher+ActsTrkFit (0.06,0.03) (4K)");
0093   col.push_back(kBlack);
0094   */
0095 
0096   /*
0097  finvec.push_back("root_files/ntp_track_out_4k_hough_genfitprop_genfit.root");
0098   legvec.push_back( "Hough+GenfitTrkProp+Genfit (4K)");
0099   col.push_back(kBlue);
0100   */
0101 
0102   /*
0103   finvec.push_back("root_files/ntp_track_out_4k_tune2_fixed_mag24_both_stub_matcher.root");
0104   legvec.push_back( "TpcTracker+stub matcher+ActsTrkFit (0.02,0.01) (4K)");
0105   col.push_back(kMagenta);
0106   */
0107 
0108   //==========
0109 
0110 
0111   unsigned int nfiles = finvec.size();
0112 
0113   TCanvas *ctemp0 = new TCanvas("ctemp0","ctemp0",5,5,1600,800);
0114   ctemp0->Divide(nfiles,1);
0115   TCanvas *cslice = new TCanvas("cslice","cslice",5,5,1600,800);
0116   cslice->Divide(nfiles,1);
0117   TCanvas *cpt = new TCanvas("cpt","cpt",5,5,1500,800); 
0118   cpt->Divide(2,1);
0119   cpt->cd(1); gPad->SetLeftMargin(0.2);
0120   cpt->cd(2); gPad->SetLeftMargin(0.2);
0121   TCanvas *ceff = new TCanvas("ceff","ceff",5,5,1200,800);
0122   ceff->SetLeftMargin(0.18);
0123   TCanvas *ctemp1 = new TCanvas("ctemp1","ctemp1",5,5,1200,800);
0124   ctemp1->Divide(nfiles,1);
0125   TCanvas *cslicexy = new TCanvas("cslicexy","cslicexy",5,5,1600,800);
0126   cslicexy->Divide(nfiles,1);
0127   TCanvas *cdcaxy = new TCanvas("cdcaxy","cdcaxy",5,5,1200,800); 
0128   cdcaxy->SetLeftMargin(0.2);
0129   gPad->SetGrid();
0130   TCanvas *ctemp2 = new TCanvas("ctemp2","ctemp2",5,5,1200,800);
0131   ctemp2->Divide(nfiles,1);
0132   TCanvas *csliceZ = new TCanvas("csliceZ","csliceZ",5,5,1600,800);
0133   csliceZ->Divide(nfiles,1);
0134   TCanvas *cdcaZ = new TCanvas("cdcaZ","cdcaZ",5,5,1200,800); 
0135   cdcaZ->SetLeftMargin(0.2);
0136   gPad->SetGrid();
0137   TCanvas *cmvtx= new TCanvas("cmvtx","cmvtx",5,5,1200,800);
0138   cmvtx->SetLeftMargin(0.2);
0139   TCanvas *ccomb= new TCanvas("ccomb","ccomb",5,5,1200,800);
0140   ccomb->SetLeftMargin(0.2);
0141 
0142   TLegend *lpd = new TLegend(0.25, 0.6, 0.85, 0.75, "", "NDC");
0143   lpd->SetBorderSize(1);
0144   lpd->SetFillColor(kWhite);
0145   lpd->SetFillStyle(1001);
0146 
0147   TF1 *fz = new TF1("fz","gaus");
0148   fz->SetRange(-0.02, 0.02);
0149 
0150   for(unsigned int i=0; i<finvec.size();++i)
0151     {
0152       TFile *fin = new TFile(finvec[i].c_str());
0153 
0154       // pT resolution
0155 
0156       TH2D *hpt2d;
0157       fin->GetObject("h1",hpt2d); 
0158       std::cout << hpt2d->GetEntries() << std::endl;
0159       hpt2d->FitSlicesY();
0160       TH1D *hptres =  (TH1D*)gDirectory->Get("h1_2");
0161       TH1D *hptcent = (TH1D*)gDirectory->Get("h1_1");
0162       ctemp0->cd(i+1);
0163       hpt2d->Draw("colz");
0164 
0165       int binlow = hpt2d->GetXaxis()->FindBin(slice_low);
0166       int binhigh = hpt2d->GetXaxis()->FindBin(slice_high);
0167       char hname[500];
0168       sprintf(hname,"hptslice%i",i);
0169       TH1D *hptslice = hpt2d->ProjectionY(hname,binlow, binhigh);
0170       cslice->cd(i+1);
0171       hptslice->DrawCopy();
0172       hptslice->Fit("gaus");
0173 
0174       cpt->cd(1);
0175       hptres->GetYaxis()->SetTitleOffset(2.1);
0176       hptres->GetXaxis()->SetTitleOffset(1.2);
0177       hptres->SetMarkerStyle(20);
0178       hptres->SetMarkerSize(1.2);
0179       hptres->SetMarkerColor(col[i]);
0180       hptres->GetXaxis()->SetRangeUser(0, ptmax);
0181       hptres->GetYaxis()->SetRangeUser(0.0, 0.08);
0182       hptres->SetTitle(";p_{T} [GeV/c];#frac{#Delta p_{T}}{p_{T}} (resolution)");
0183       if(i == 0)
0184     hptres->DrawCopy("p");
0185       else
0186     hptres->DrawCopy("p same");
0187 
0188       lpd->AddEntry(hptres, legvec[i].c_str());
0189 
0190       cpt->cd(2);
0191       hptcent->GetYaxis()->SetTitleOffset(2.1);
0192       hptcent->GetXaxis()->SetTitleOffset(1.2);
0193       hptcent->SetMarkerStyle(20);
0194       hptcent->SetMarkerSize(1.2);
0195       hptcent->SetMarkerColor(col[i]);
0196       hptcent->GetXaxis()->SetRangeUser(0, ptmax);
0197       hptcent->GetYaxis()->SetRangeUser(-0.1, 0.1);
0198       hptcent->SetTitle(";p_{T} [GeV/c];#frac{#Delta p_{T}}{p_{T}} (offset)");
0199       if(i==0)
0200     hptcent->DrawCopy("p");
0201       else
0202     hptcent->DrawCopy("p same");
0203 
0204       // Efficiency
0205       ceff->cd();
0206       TH1D *hnum = 0;
0207       TH1D *hden = 0;
0208       fin->GetObject("h3_num",hnum);
0209       fin->GetObject("h3_den",hden);
0210 
0211       TH1D* heff = (TH1D*)hden->Clone("heff");;
0212 
0213       for(int i=1;i<=hden->GetNbinsX();++i)
0214     {
0215       double pass = hnum->GetBinContent(i);
0216       double all = hden->GetBinContent(i);
0217       double eff = 0;
0218       if(all > pass)
0219         eff = pass/all;
0220       else if(all > 0)
0221         eff = 1.;
0222       heff->SetBinContent(i, eff);
0223     }  
0224 
0225       heff->GetYaxis()->SetTitle("Efficiency");
0226       heff->GetYaxis()->SetTitleOffset(1.5);
0227       heff->GetXaxis()->SetTitle("p_{T} (GeV/c)");
0228       heff->GetXaxis()->SetTitleOffset(1.2);
0229       heff->SetMarkerStyle(20);
0230       heff->SetMarkerSize(1.2);
0231       heff->SetMarkerColor(col[i]);
0232       heff->GetXaxis()->SetRangeUser(0.0, ptmax);
0233       heff->GetYaxis()->SetRangeUser(0.0, 1.05);
0234       if(i==0)
0235     heff->DrawCopy("p");
0236       else
0237     heff->DrawCopy("p same");
0238 
0239       TLine *unit = new TLine(0,1.0,40,1.0);
0240       unit->Draw();
0241 
0242       // dca xy resolution
0243       cdcaxy->cd();
0244       TH2D *hdca2d;
0245       fin->GetObject("h2",hdca2d);
0246       ctemp1->cd(i+1);
0247       hdca2d->Draw("colz");
0248 
0249       binlow = hdca2d->GetXaxis()->FindBin(slice_low);
0250       binhigh = hdca2d->GetXaxis()->FindBin(slice_high);
0251       TH1D* hdcaxyslice = hdca2d->ProjectionY("hdcaxyslice",binlow, binhigh);
0252       hdcaxyslice->SetLineColor(col[i]);      
0253       cslicexy->cd(i+1);
0254       hdcaxyslice->Draw();
0255       hdcaxyslice->Fit("gaus");
0256 
0257       cdcaxy->cd();  
0258       hdca2d->FitSlicesY();
0259       TH1D*hdcares = (TH1D*)gDirectory->Get("h2_2");
0260       hdcares->GetYaxis()->SetTitleOffset(2.1);
0261       hdcares->GetYaxis()->SetTitle("DCA(r#phi) (cm)");
0262       hdcares->GetXaxis()->SetTitleOffset(1.2);
0263       hdcares->GetXaxis()->SetTitle("p_{T} (GeV/c)");
0264       hdcares->SetMarkerStyle(20);
0265       hdcares->SetMarkerSize(1.2);
0266       hdcares->SetMarkerColor(col[i]);
0267       hdcares->GetYaxis()->SetRangeUser(0, 0.008);
0268       if(i==0) 
0269     hdcares->DrawCopy("p");
0270       else
0271     hdcares->DrawCopy("p same");
0272 
0273       // dca z resolution
0274       TH2D *hdcaZ2d;
0275       fin->GetObject("h3",hdcaZ2d);
0276       ctemp2->cd(i+1);
0277       hdcaZ2d->Draw("colz");
0278 
0279       binlow = hdcaZ2d->GetXaxis()->FindBin(slice_low);
0280       binhigh = hdcaZ2d->GetXaxis()->FindBin(slice_high);
0281       TH1D* hdcazslice = hdcaZ2d->ProjectionY("hdcazslice",binlow, binhigh);
0282       hdcazslice->SetLineColor(col[i]);      
0283       csliceZ->cd(i+1);
0284       hdcazslice->Draw();
0285       hdcazslice->Fit("gaus");
0286 
0287       cdcaZ->cd();
0288       hdcaZ2d->FitSlicesY(fz);
0289       TH1D* hdcaZres = (TH1D*)gDirectory->Get("h3_2");
0290 
0291       hdcaZres->GetYaxis()->SetTitleOffset(2.1);
0292       hdcaZres->GetYaxis()->SetTitle("DCA(Z) (cm)");
0293       hdcaZres->GetXaxis()->SetTitleOffset(1.2);
0294       hdcaZres->GetXaxis()->SetTitle("p_{T} (GeV/c");
0295       hdcaZres->SetMarkerStyle(20);
0296       hdcaZres->SetMarkerSize(1.2);
0297       hdcaZres->SetMarkerColor(col[i]);
0298       hdcaZres->GetYaxis()->SetRangeUser(0, 0.015);
0299       if(i==0)      
0300     hdcaZres->DrawCopy("p");
0301       else
0302     hdcaZres->DrawCopy("p same");
0303 
0304       // Plot the fraction of tracks matched to the MVTX
0305       cmvtx->cd();
0306       TH2D *hnmvtx = 0;
0307       fin->GetObject("h6",hnmvtx);
0308       if(!hnmvtx) cout << " Did not get hnmvtx" << endl;
0309 
0310       //hnmvtx->DrawCopy("p");
0311 
0312       binlow = hnmvtx->GetYaxis()->FindBin(0.5);
0313       binhigh = hnmvtx->GetYaxis()->FindBin(3.5);
0314       cout << "binlow " << binlow << " binhigh " << binhigh << endl;
0315       TH1D *hnmaps_hit = hnmvtx->ProjectionX("hnmaps_hit",binlow, binhigh);;
0316 
0317       binlow = hnmvtx->GetYaxis()->FindBin(-0.5);
0318       binhigh = hnmvtx->GetYaxis()->FindBin(3.5);
0319       cout << "binlow " << binlow << " binhigh " << binhigh << endl;
0320       TH1D *hnmaps_all = hnmvtx->ProjectionX("hnmaps_all",binlow, binhigh);;
0321 
0322       TH1D *hnmaps = (TH1D*) hnmaps_all->Clone();
0323       for(int i=1;i<=hnmaps_hit->GetNbinsX();++i)
0324     {
0325       double hit = hnmaps_hit->GetBinContent(i);
0326       double all = hnmaps_all->GetBinContent(i);
0327       double eff = 0;
0328       if(all > hit)
0329         eff = hit/all;
0330       else if(all > 0)
0331         eff = 1.;
0332       hnmaps->SetBinContent(i, eff);
0333     }  
0334 
0335       hnmaps->GetXaxis()->SetTitle("p_{T} (GeV/c)");
0336       hnmaps->GetYaxis()->SetTitle("MVTX association efficiency");
0337       hnmaps->GetYaxis()->SetTitleOffset(1.5);
0338       hnmaps->GetXaxis()->SetTitleOffset(1.2);
0339       hnmaps->SetMinimum(0);
0340       hnmaps->SetMaximum(1.1);
0341       hnmaps->SetMarkerStyle(20);
0342       hnmaps->SetMarkerSize(1.5);
0343       hnmaps->SetMarkerColor(col[i]);
0344 
0345       if(i==0)
0346     {
0347       hnmaps->DrawCopy("p");
0348       unit->Draw();
0349     }
0350       else
0351       hnmaps->DrawCopy("p same");
0352 
0353       // Combine track efficiency and MVTX association efficiency
0354 
0355       ccomb->cd();
0356 
0357       cout << " hnmaps bins " << hnmaps->GetNbinsX() << " heff bins " << heff->GetNbinsX() << endl;
0358       heff->Rebin(3);
0359       cout << heff->GetNbinsX() << endl;
0360 
0361       TH1D *hcomb = (TH1D*) hnmaps->Clone();
0362       for(int icomb=1; icomb< hnmaps->GetNbinsX(); ++icomb)
0363     {
0364       double vmaps = hnmaps->GetBinContent(icomb);
0365       double vtrack = heff->GetBinContent(icomb) / 3;
0366       double eff = 0;
0367       eff = vmaps * vtrack;
0368       //cout << " icomb " << icomb << " vmaps " << vmaps << " vtrack " << vtrack << " comb " << eff << endl;
0369       hcomb->SetBinContent(icomb, eff);
0370     }  
0371 
0372       hcomb->GetXaxis()->SetRangeUser(0,39);      
0373       hcomb->GetYaxis()->SetTitle("Efficincy of MVTX associated tracks");
0374 
0375       TLine *unit2 = new TLine(0,1.0,39,1.0);
0376 
0377       if(i==0)
0378     {
0379       hcomb->DrawCopy("p");
0380       unit2->Draw();
0381     }
0382       else
0383       hcomb->DrawCopy("p same");
0384 
0385     }
0386   
0387   cpt->cd(2);  lpd->Draw();
0388 
0389   cdcaxy->cd();  lpd->Draw();
0390 
0391   cdcaZ->cd();  lpd->Draw();
0392 
0393   ceff->cd();  lpd->Draw();
0394 
0395   cmvtx->cd();  lpd->Draw();
0396 
0397   ccomb->cd();  lpd->Draw();
0398 
0399 
0400 
0401 }
0402