Back to home page

sPhenix code displayed by LXR
 
 

    

Warning, file /analysis/HF-Jet/FastSimPP/macros/DrawTrackingEval.C was not indexed or was modified since last indexation (in which case cross-reference links may be missing, inaccurate or erroneous).

0001 #include "SaveCanvas.C"
0002 #include "sPhenixStyle.C"
0003 
0004 #include <TChain.h>
0005 #include <TCut.h>
0006 #include <TEfficiency.h>
0007 #include <TF1.h>
0008 #include <TGraphAsymmErrors.h>
0009 #include <TGraphErrors.h>
0010 #include <TH2.h>
0011 #include <TH3.h>
0012 
0013 #include <TFile.h>
0014 
0015 #include <TColor.h>
0016 #include <TLatex.h>
0017 #include <TLegend.h>
0018 #include <TLine.h>
0019 
0020 #include <TMath.h>
0021 #include <TString.h>
0022 #include <TTree.h>
0023 #include <TVirtualFitter.h>
0024 
0025 #include <cmath>
0026 #include <iostream>
0027 
0028 using namespace std;
0029 
0030 // ROOT6 disabled assert. Well....
0031 #ifdef assert
0032 #undef assert
0033 #endif
0034 #define assert(exp)                                                                             \
0035   {                                                                                             \
0036     if (!(exp))                                                                                 \
0037     {                                                                                           \
0038       cout << "Assert (" << #exp << ") failed at " << __FILE__ << " line " << __LINE__ << endl; \
0039       exit(1);                                                                                  \
0040     }                                                                                           \
0041   }
0042 
0043 TTree *ntp_gtrack(nullptr);
0044 TFile *_file0 = NULL;
0045 TString description;
0046 TCut primPartCut("gpt>.1 && abs(geta)<1 && abs(gvt)<1e-10");  // primary particle in acceptance
0047 TCut recoCut("nlmaps>=2 && pt>0");
0048 const double pTMin(.3);
0049 
0050 //#define ASSERT(exp) {{exit(1);} }
0051 
0052 //! utility function to
0053 void useLogBins(TAxis *axis)
0054 {
0055   assert(axis);
0056   assert(axis->GetXmin() > 0);
0057   assert(axis->GetXmax() > 0);
0058 
0059   const int bins = axis->GetNbins();
0060 
0061   Axis_t from = log10(axis->GetXmin());
0062   Axis_t to = log10(axis->GetXmax());
0063   Axis_t width = (to - from) / bins;
0064   vector<Axis_t> new_bins(bins + 1);
0065 
0066   for (int i = 0; i <= bins; i++)
0067   {
0068     new_bins[i] = TMath::Power(10, from + i * width);
0069   }
0070   axis->Set(bins, new_bins.data());
0071 }
0072 
0073 TGraphErrors *
0074 FitProfile(const TH2 *h2)
0075 {
0076   TProfile *p2 = h2->ProfileX();
0077 
0078   int n = 0;
0079   double x[1000];
0080   double ex[1000];
0081   double y[1000];
0082   double ey[1000];
0083 
0084   for (int i = 1; i <= h2->GetNbinsX(); i++)
0085   {
0086     TH1D *h1 = h2->ProjectionY(Form("htmp_%d", rand()), i, i);
0087 
0088     if (h1->GetSum() < 30)
0089     {
0090       cout << "FitProfile - ignore bin " << i << endl;
0091       continue;
0092     }
0093     else
0094     {
0095       //      cout << "FitProfile - fit bin " << i << endl;
0096     }
0097 
0098     TF1 fgaus("fgaus", "gaus", -p2->GetBinError(i) * 4,
0099               p2->GetBinError(i) * 4);
0100 
0101     //    TF1 f2(Form("dgaus"), "gaus + [3]*exp(-0.5*((x-[1])/[4])**2) + [5]",
0102     //           -p2->GetBinError(i) * 4, p2->GetBinError(i) * 4);
0103 
0104     fgaus.SetParameter(1, p2->GetBinContent(i));
0105     fgaus.SetParameter(2, p2->GetBinError(i));
0106 
0107     h1->Fit(&fgaus, "MQ");
0108 
0109     TF1 f2(Form("dgaus"), "gaus  + [3]",
0110            -fgaus.GetParameter(2) * 1.5, fgaus.GetParameter(2) * 1.5);
0111 
0112     f2.SetParameters(fgaus.GetParameter(0), fgaus.GetParameter(1),
0113                      fgaus.GetParameter(2), 0);
0114 
0115     h1->Fit(&f2, "MQR");
0116 
0117     //      new TCanvas;
0118     //      h1->Draw();
0119     //      fgaus.Draw("same");
0120     //      break;
0121 
0122     x[n] = p2->GetBinCenter(i);
0123     ex[n] = (p2->GetBinCenter(2) - p2->GetBinCenter(1)) / 2;
0124     y[n] = fgaus.GetParameter(2);
0125     ey[n] = fgaus.GetParError(2);
0126 
0127     //      p2->SetBinContent(i, fgaus.GetParameter(1));
0128     //      p2->SetBinError(i, fgaus.GetParameter(2));
0129 
0130     n++;
0131     delete h1;
0132   }
0133 
0134   return new TGraphErrors(n, x, y, ex, ey);
0135 }
0136 
0137 void RecoCheck(const TString &name = "pi", const TCut &cut = "abs(gflavor)==211", const int nBins = 100)
0138 {
0139   assert(_file0);
0140 
0141   TCanvas *c1 = new TCanvas("RecoCheck_" + name, "RecoCheck_" + name, 1800, 960);
0142   c1->Divide(3, 2);
0143   int idx = 1;
0144   TPad *p;
0145 
0146   p = (TPad *) c1->cd(idx++);
0147   c1->Update();
0148   p->SetLogx();
0149 
0150   p->DrawFrame(pTMin, 0, 20, 1.1, "Efficiency;Truth p_{T} [GeV/c];Tracking efficiency");
0151   TLine *l = new TLine(pTMin, 1, 20, 1);
0152   l->Draw();
0153 
0154   TH1 *hTruthpT = new TH1F("hTruthpT", "hTruthpT;Truth p_{T} [GeV/c]", nBins, pTMin, 20);
0155   useLogBins(hTruthpT->GetXaxis());
0156   TH1 *hRecopT = new TH1F("hRecopT", "hRecopT;Truth p_{T} [GeV/c]", nBins, pTMin, 20);
0157   useLogBins(hRecopT->GetXaxis());
0158   TH1 *hRecopTMVTX = new TH1F("hRecopTMVTX", "hRecopT;Truth p_{T} [GeV/c]", nBins, pTMin, 20);
0159   useLogBins(hRecopTMVTX->GetXaxis());
0160 
0161   ntp_gtrack->Draw("gpt>>hTruthpT", cut && primPartCut, "goff");
0162   ntp_gtrack->Draw("gpt>>hRecopT", cut && primPartCut && "pt>0", "goff");
0163   ntp_gtrack->Draw("gpt>>hRecopTMVTX", cut && primPartCut && recoCut, "goff");
0164 
0165   TEfficiency *hEffpTMVTX = new TEfficiency(*hRecopTMVTX, *hTruthpT);
0166   hEffpTMVTX->SetName("hEffpTMVTX");
0167   hEffpTMVTX->SetTitle("Efficiency;Truth p_{T} [GeV/c];Tracking efficiency");
0168   hEffpTMVTX->Draw("p same");
0169   hEffpTMVTX->SetLineColor(kBlue + 2);
0170   hEffpTMVTX->SetMarkerColor(kBlue + 2);
0171 
0172   //  TH1 * hDCAxy = new TH1F("hDCAxy","")
0173   TEfficiency *hEffpT = new TEfficiency(*hRecopT, *hTruthpT);
0174   hEffpT->SetName("hEffpT");
0175   hEffpT->SetTitle("Efficiency;Truth p_{T} [GeV/c];Tracking efficiency");
0176   hEffpT->Draw("P same");
0177   hEffpT->SetLineColor(kRed + 2);
0178   hEffpT->SetMarkerColor(kRed + 2);
0179 
0180   //  TH2F * hDCA;
0181   //  ntp_gtrack -> Draw();
0182 
0183   p = (TPad *) c1->cd(idx++);
0184   c1->Update();
0185   p->SetLogx();
0186   p->SetLogz();
0187 
0188   TH2 *hpTResolution = new TH2F("hpTResolution", "hpTResolution;Truth p_{T} [GeV/c];#Deltap_{T}/p_{T}", nBins, pTMin, 20, 1000, -.2, .2);
0189   useLogBins(hpTResolution->GetXaxis());
0190   hpTResolution->GetYaxis()->SetRangeUser(-.1, .1);
0191 
0192   ntp_gtrack->Draw("pt/gpt - 1 : gpt>>hpTResolution", cut && primPartCut && recoCut, "goff");
0193 
0194   TGraphErrors *gepTResolution = FitProfile(hpTResolution);
0195   gepTResolution->SetName("gepTResolution");
0196 
0197   hpTResolution->Draw("colz");
0198   gepTResolution->Draw("p");
0199 
0200   p = (TPad *) c1->cd(idx++);
0201   c1->Update();
0202   p->SetLogx();
0203   p->SetLogz();
0204 
0205   TH2 *hcotThetaResolution = new TH2F("hcotThetaResolution", "hcotThetaResolution;Truth p_{T} [GeV/c];#Delta[p_{z}/p_{T}]", nBins, pTMin, 20, 1000, -.1, .1);
0206   useLogBins(hcotThetaResolution->GetXaxis());
0207   hcotThetaResolution->GetYaxis()->SetRangeUser(-.01, .01);
0208 
0209   ntp_gtrack->Draw("pz/pt - gpz/gpt: gpt>>hcotThetaResolution", cut && primPartCut && recoCut, "goff");
0210 
0211   TGraphErrors *gecotThetaResolution = FitProfile(hcotThetaResolution);
0212   gecotThetaResolution->SetName("gecotThetaResolution");
0213 
0214   hcotThetaResolution->Draw("colz");
0215   gecotThetaResolution->Draw("p");
0216 
0217   p = (TPad *) c1->cd(idx++);
0218   c1->Update();
0219   p->SetLogx();
0220   p->SetLogz();
0221 
0222   TH2 *hsinPhiResolution = new TH2F("hsinPhiResolution", "hsinPhiResolution;Truth p_{T} [GeV/c];sin(#Delta#phi)", nBins, pTMin, 20, 1000, -.1, .1);
0223   useLogBins(hsinPhiResolution->GetXaxis());
0224   hsinPhiResolution->GetYaxis()->SetRangeUser(-.01, .01);
0225 
0226   ntp_gtrack->Draw("sin(atan2(py,px) - atan2(gpy,gpx)) : gpt>>hsinPhiResolution", cut && primPartCut && recoCut, "goff");
0227 
0228   TGraphErrors *gesinPhiResolution = FitProfile(hsinPhiResolution);
0229   gesinPhiResolution->SetName("gesinPhiResolution");
0230 
0231   hsinPhiResolution->Draw("colz");
0232   gesinPhiResolution->Draw("p");
0233 
0234   p = (TPad *) c1->cd(idx++);
0235   c1->Update();
0236   p->SetLogx();
0237   p->SetLogz();
0238 
0239   TH2 *hDCA3DXYResolution = new TH2F("hDCA3DXYResolution", "hDCA3DXYResolution;Truth p_{T} [GeV/c];DCA_{xy} [cm]", nBins, pTMin, 20, 1000, -.2, .2);
0240   useLogBins(hDCA3DXYResolution->GetXaxis());
0241   hDCA3DXYResolution->GetYaxis()->SetRangeUser(-.02, .02);
0242 
0243   ntp_gtrack->Draw("dca3dxy : gpt>>hDCA3DXYResolution", cut && primPartCut && recoCut, "goff");
0244 
0245   TGraphErrors *geDCA3DXYResolution = FitProfile(hDCA3DXYResolution);
0246   geDCA3DXYResolution->SetName("geDCA3DXYResolution");
0247 
0248   hDCA3DXYResolution->Draw("colz");
0249   geDCA3DXYResolution->Draw("p");
0250 
0251   p = (TPad *) c1->cd(idx++);
0252   c1->Update();
0253   p->SetLogx();
0254   p->SetLogz();
0255 
0256   TH2 *hDCA3DZResolution = new TH2F("hDCA3DZResolution", "hDCA3DZResolution;Truth p_{T} [GeV/c];DCA_{z} [cm]", nBins, pTMin, 20, 1000, -.2, .2);
0257   useLogBins(hDCA3DZResolution->GetXaxis());
0258   hDCA3DZResolution->GetYaxis()->SetRangeUser(-.02, .02);
0259 
0260   ntp_gtrack->Draw("dca3dz : gpt>>hDCA3DZResolution", cut && primPartCut && recoCut, "goff");
0261 
0262   TGraphErrors *geDCA3DZResolution = FitProfile(hDCA3DZResolution);
0263   geDCA3DZResolution->SetName("geDCA3DZResolution");
0264 
0265   hDCA3DZResolution->Draw("colz");
0266   geDCA3DZResolution->Draw("p");
0267 
0268   SaveCanvas(c1, TString(_file0->GetName()) + TString(c1->GetName()), false);
0269 }
0270 
0271 void DCACheck(const TString &name = "pi", const TCut &cut = "abs(gflavor)==211", const int nBins = 14)
0272 {
0273   assert(_file0);
0274 
0275   const int nDCAbins = 2 * (100 / 1 * 2 + (1000 - 100) / 10 + (2000 - 1000) / 100);
0276 
0277   TH3F *hDCA3Dpz = new TH3F("hDCA3Dpz", "hDCA3Dpz;DCA_{xy} [cm];DCA_{z} [cm];Truth p_{T} [GeV/c]",  //
0278                             nDCAbins, -.2, .2,                                                      //
0279                             nDCAbins, -.2, .2,                                                      //
0280                             nBins, pTMin, 20);
0281 
0282   vector<Axis_t> new_bins(nDCAbins + 1);
0283   int binedge = -20000;
0284   for (int i = 0; i <= nDCAbins; i++)
0285   {
0286     new_bins[i] = binedge * 1e-5;  // convert .1um to cm
0287 
0288     double dedge = 1000;
0289     if (abs(binedge + .5) <= 10000)
0290       dedge = 100;
0291     if (abs(binedge + .5) <= 1000)
0292       dedge = 5;
0293 
0294     //    cout << "new_bins[" << i << "] = " << binedge << "* 1e-4 = " << new_bins[i] << ", dedge = " << dedge << endl;
0295 
0296     binedge += dedge;
0297   }
0298   assert(binedge == 21000);
0299   hDCA3Dpz->GetXaxis()->Set(nDCAbins, new_bins.data());
0300   hDCA3Dpz->GetYaxis()->Set(nDCAbins, new_bins.data());
0301   //    cout << "hDCA3Dpz->GetXaxis()->GetBinCenter(nDCAbins) = " << hDCA3Dpz->GetXaxis()->GetBinCenter(nDCAbins) << endl;
0302 
0303   //  ntp_gtrack->Draw("dca3dxy : dca3dz : gpt >> hDCA3Dpz", cut && primPartCut && recoCut, "goff");
0304   ntp_gtrack->Draw("gpt : dca3dz : dca3dxy >> hDCA3Dpz", cut && primPartCut && recoCut, "goff");
0305 
0306   // plot
0307   TCanvas *c1 = new TCanvas("DCACheck_" + name, "DCACheck_" + name, 1800, 960);
0308 
0309   c1->Divide(5, 3);
0310   int idx = 1;
0311   TPad *p;
0312 
0313   p = (TPad *) c1->cd(idx++);
0314   c1->Update();
0315   p->SetLogx();
0316   p->SetLogz();
0317 
0318   TH2 *hDCA3DXYCheck = new TH2F("hDCA3DXYCheck", "hDCA3DXYCheck;Truth p_{T} [GeV/c];DCA_{xy}", nBins, pTMin, 20, 1000, -.2, .2);
0319   useLogBins(hDCA3DXYCheck->GetXaxis());
0320   hDCA3DXYCheck->GetYaxis()->SetRangeUser(-.02, .02);
0321 
0322   ntp_gtrack->Draw("dca3dxy : gpt>>hDCA3DXYCheck", cut && primPartCut && recoCut, "goff");
0323 
0324   TGraphErrors *geDCA3DXYCheck = FitProfile(hDCA3DXYCheck);
0325   geDCA3DXYCheck->SetName("geDCA3DXYCheck");
0326 
0327   hDCA3DXYCheck->Draw("colz");
0328   geDCA3DXYCheck->Draw("p");
0329 
0330   for (int bin = 1; bin <= nBins; bin++)
0331   {
0332     p = (TPad *) c1->cd(idx++);
0333     c1->Update();
0334     //    p->SetLogx();
0335     p->SetLogz();
0336     p->SetRightMargin(.15);
0337 
0338     hDCA3Dpz->GetZaxis()->SetRange(bin, bin);
0339     TH1 *h2D = hDCA3Dpz->Project3D("yx");
0340     h2D->SetName(Form("hDCA3Dpz%d", bin));
0341     h2D->Draw("colz");
0342   }
0343 
0344   // ping the 3D histogram for saving
0345   c1->GetListOfPrimitives()->Add(hDCA3Dpz);
0346 
0347   SaveCanvas(c1, TString(_file0->GetName()) + TString(c1->GetName()), false);
0348 }
0349 
0350 void DrawTrackingEval(
0351     const TString infile = "/gpfs/mnt/gpfs02/sphenix/user/jinhuang/HF-jet/HF-production-piKp-pp200-dataset1_333ALL.cfg_ALL_g4svtx_eval.root",
0352     const TString &name = "pi", const TCut &cut = "abs(gflavor)==211"  //
0353 )
0354 {
0355   SetsPhenixStyle();
0356   TVirtualFitter::SetDefaultFitter("Minuit2");
0357 
0358 //  description = disc;
0359   gSystem->Load("libg4eval.so");
0360 
0361   if (!_file0)
0362   {
0363     TString chian_str = infile;
0364     chian_str.ReplaceAll("ALL", "*");
0365 
0366     TChain *t = new TChain("ntp_gtrack");
0367     const int n = t->Add(chian_str);
0368 
0369     cout << "Loaded " << n << " root files with " << chian_str << endl;
0370     assert(n > 0);
0371     ntp_gtrack = t;
0372 
0373     _file0 = new TFile;
0374     _file0->SetName(infile);
0375   }
0376 
0377 //  gDirectory->mkdir("/pi");
0378 //  gDirectory->Cd("/pi");
0379   RecoCheck(name, cut);
0380   DCACheck(name, cut);
0381 
0382   //  gDirectory->mkdir("/kaon");
0383   //  gDirectory->Cd("/kaon");
0384   //  RecoCheck("kaon", "abs(gflavor)==321");
0385   //  DCACheck("kaon", "abs(gflavor)==321");
0386   //
0387   //  gDirectory->mkdir("/proton");
0388   //  gDirectory->Cd("/proton");
0389   //  RecoCheck("proton", "abs(gflavor)==2212");
0390   //  DCACheck("proton", "abs(gflavor)==2212");
0391 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!