Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2025-08-05 08:13:30 

0001 #ifndef VERTEX_H
0002 #define VERTEX_H
0003 
0004 #include <Riostream.h>
0005 #include <fstream>
0006 #include <iostream>
0007 #include <map>
0008 #include <numeric>
0009 #include <stdlib.h>
0010 #include <vector>
0011 #include <tuple>
0012 
0013 #include <TCanvas.h>
0014 #include <TFile.h>
0015 #include <TGaxis.h>
0016 #include <TH1F.h>
0017 #include <TH2F.h>
0018 #include <TLatex.h>
0019 #include <TLegend.h>
0020 #include <TLine.h>
0021 #include <TMarker.h>
0022 #include <TRandom3.h>
0023 #include <TText.h>
0024 #include <TTree.h>
0025 #include <TVector3.h>
0026 
0027 #include "Utilities.h"
0028 #include "Hit.h"
0029 
0030 #include "/sphenix/user/hjheng/TrackletAna/analysis/plot/sPHENIXStyle/sPhenixStyle.C"
0031 
0032 #define NZGAP 0
0033 
0034 int Ncluster_toplot = 1300;
0035 
0036 using namespace std;
0037 
0038 struct MvtxClusComb
0039 {
0040     MvtxClusComb(float x1, float y1, float z1, float x2, float y2, float z2, float extz) : x1(x1), y1(y1), z1(z1), x2(x2), y2(y2), z2(z2), extz(extz){};
0041 
0042     float x1, y1, z1, x2, y2, z2, extz;
0043 };
0044 
0045 struct Cluster
0046 {
0047     Cluster(uint32_t index, uint32_t nz) : index(index), nz(nz){};
0048 
0049     uint32_t index;
0050     uint32_t nz;
0051 };
0052 
0053 struct Vertex
0054 {
0055     Vertex(float vz, float sigma2, vector<float> hitx1, vector<float> hity1, vector<float> hitz1, vector<float> hitx2, vector<float> hity2, vector<float> hitz2, vector<float> candz)
0056         : vz(vz), sigma2(sigma2), hitx1(hitx1), hity1(hity1), hitz1(hitz1), hitx2(hitx2), hity2(hity2), hitz2(hitz2), candz(candz){};
0057 
0058     float vz;
0059     float sigma2;
0060     vector<float> hitx1, hity1, hitz1, hitx2, hity2, hitz2, candz;
0061 };
0062 
0063 struct VtxData
0064 {
0065     bool isdata;
0066     int event, NhitsLayer1, NTruthVtx;
0067     float TruthPV_trig_x, TruthPV_trig_y, TruthPV_trig_z, TruthPV_mostNpart_x, TruthPV_mostNpart_y, TruthPV_mostNpart_z;
0068     // PV positions from two MVTX halves
0069     float PV_x_halves1, PV_y_halves1, PV_z_halves1, PV_x_halves2, PV_y_halves2, PV_z_halves2;
0070     float Centrality_bimp, Centrality_impactparam, Centrality_mbd, Centrality_mbdquantity;
0071 };
0072 
0073 // https://github.com/bi-ran/tracklet/blob/master/analysis/include/tracklet.h
0074 // float TrackletPV_cluster(int evt, vector<Hit *> layer1, vector<Hit *> layer2, int dPhiCutbin, int dZCutbin, bool verbose)
0075 // {
0076 //     float dPhi_cut = dPhiCutbin * 0.01;
0077 //     float dZ_cut = dZCutbin * 0.01;
0078 
0079 //     float vertex = -999.;
0080 
0081 //     vector<double> vertices;
0082 //     vertices.clear();
0083 
0084 //     for (size_t i = 0; i < layer1.size(); i++)
0085 //     {
0086 //         for (size_t j = 0; j < layer2.size(); j++)
0087 //         {
0088 //             if (fabs(deltaPhi(layer1[i]->Phi(), layer2[j]->Phi())) > dPhi_cut)
0089 //                 continue;
0090 
0091 //             float z = layer1[i]->posZ() - (layer2[j]->posZ() - layer1[i]->posZ()) / (layer2[j]->rho() - layer1[i]->rho()) * layer1[i]->rho();
0092 //             vertices.push_back(z);
0093 //         }
0094 //     }
0095 
0096 //     /* vertex clustering */
0097 //     if (vertices.size())
0098 //     {
0099 //         sort(vertices.begin(), vertices.end());
0100 
0101 //         std::vector<Cluster> clusters;
0102 //         std::size_t z = 0, y = 0;
0103 //         for (; z < vertices.size(); ++z)
0104 //         {
0105 //             for (; y < vertices.size() && vertices[y] - vertices[z] < dZ_cut; ++y)
0106 //                 ;
0107 //             clusters.emplace_back(z, y - z); // index, nz
0108 //         }
0109 
0110 //         std::sort(clusters.begin(), clusters.end(), [](const Cluster &a, const Cluster &b) -> bool { return a.nz > b.nz; });
0111 
0112 //         std::vector<Vertex> candidates;
0113 //         uint32_t maxnz = clusters[0].nz;
0114 
0115 //         for (auto &cluster : clusters)
0116 //         {
0117 //             uint32_t nz = cluster.nz;
0118 //             if (nz + NZGAP < maxnz)
0119 //             {
0120 //                 break;
0121 //             }
0122 
0123 //             uint32_t index = cluster.index;
0124 
0125 //             float vz = 0;
0126 //             float vz2 = 0;
0127 //             float sigma2 = 0;
0128 //             for (uint32_t y = 0; y < nz; ++y)
0129 //             {
0130 //                 vz += vertices[index + y];
0131 //                 vz2 += vertices[index + y] * vertices[index + y];
0132 //             }
0133 //             vz /= nz;
0134 //             sigma2 = vz2 / nz - (vz * vz);
0135 
0136 //             candidates.emplace_back(vz, sigma2);
0137 //         }
0138 
0139 //         vertex = std::min_element(candidates.begin(), candidates.end(), [](const Vertex &a, const Vertex &b) -> bool { return a.sigma2 < b.sigma2; })->vz;
0140 //     }
0141 
0142 //     return vertex;
0143 // }
0144 
0145 // TCanvas *Canv_RecoVtxZTklZRho(vector<TLine *> v_line, vector<Hit *> hit_l1, vector<Hit *> hit_l2, float PVz)
0146 void Draw_RecoVtxZTklZRho(TH2F *hM_clusZRho_half1, TH2F *hM_clusZRho_half2, vector<TLine *> vtl_half1, vector<TLine *> vtl_half2, float PVz_half1, float PVz_half2, float PVz, bool isdata, TString plotname_tstr)
0147 {
0148     SetsPhenixStyle();
0149     gStyle->SetPalette(kThermometer);
0150 
0151     int nTLines = vtl_half1.size() + vtl_half2.size();
0152     float tlinealpha = (nTLines > 500) ? 0.05 : 0.4;
0153 
0154     // merge histograms
0155     TH2F *hM_clusZRho = (TH2F *)hM_clusZRho_half1->Clone();
0156     hM_clusZRho->Add(hM_clusZRho_half2);
0157     int Nclus = hM_clusZRho->GetEntries();
0158 
0159     TCanvas *c = new TCanvas("c", "c", 800, 700);
0160     c->cd();
0161     gPad->SetRightMargin(0.08);
0162     gPad->SetTopMargin(0.05);
0163     gPad->SetLeftMargin(0.15);
0164     gPad->SetBottomMargin(0.13);
0165     hM_clusZRho->GetXaxis()->SetTitle("Position Z [cm]");
0166     hM_clusZRho->GetYaxis()->SetTitle("Radius [cm]");
0167     hM_clusZRho->GetXaxis()->SetTitleOffset(1.1);
0168     hM_clusZRho->GetYaxis()->SetTitleOffset(1.3);
0169     hM_clusZRho->SetMarkerSize(0.5);
0170     hM_clusZRho->Draw("SCAT");
0171 
0172     for (size_t i = 0; i < vtl_half1.size(); i++)
0173     {
0174         vtl_half1[i]->SetLineColorAlpha(14, tlinealpha);
0175         vtl_half1[i]->Draw("same");
0176     }
0177     for (size_t i = 0; i < vtl_half2.size(); i++)
0178     {
0179         vtl_half2[i]->SetLineColorAlpha(14, tlinealpha);
0180         vtl_half2[i]->Draw("same");
0181     }
0182 
0183     hM_clusZRho->Draw("SCATsame");
0184 
0185     // MVTX dimension
0186     float MVTXLength = 27.1; // cm
0187     float MVTXL1_min = 2.461;
0188     float MVTXL1_mid = 2.523;
0189     float MVTXL1_max = 2.793;
0190     float MVTXL2_min = 3.198;
0191     float MVTXL2_mid = 3.335;
0192     float MVTXL2_max = 3.625;
0193     TLine *beamline = new TLine(-20, 0, 20, 0);
0194     TLine *MVTXl1_in = new TLine(-(MVTXLength / 2), MVTXL1_min, MVTXLength / 2, MVTXL1_min);
0195     TLine *MVTXl1_out = new TLine(-(MVTXLength / 2), MVTXL1_max, MVTXLength / 2, MVTXL1_max);
0196     TLine *MVTXl2_in = new TLine(-(MVTXLength / 2), MVTXL2_min, MVTXLength / 2, MVTXL2_min);
0197     TLine *MVTXl2_out = new TLine(-(MVTXLength / 2), MVTXL2_max, MVTXLength / 2, MVTXL2_max);
0198     MVTXl1_in->SetLineStyle(kDashed);
0199     MVTXl1_out->SetLineStyle(kDashed);
0200     MVTXl2_in->SetLineStyle(kDashed);
0201     MVTXl2_out->SetLineStyle(kDashed);
0202     beamline->Draw("same");
0203     MVTXl1_in->Draw("same");
0204     MVTXl1_out->Draw("same");
0205     MVTXl2_in->Draw("same");
0206     MVTXl2_out->Draw("same");
0207 
0208     TMarker *m_PVz_half1 = new TMarker(PVz_half1, 0.0, 20);
0209     TMarker *m_PVz_half2 = new TMarker(PVz_half2, 0.0, 20);
0210     TMarker *m_PVz = new TMarker(PVz, 0.0, 20);
0211     m_PVz_half1->SetMarkerSize(1);
0212     m_PVz_half1->SetMarkerColor(38);
0213     m_PVz_half2->SetMarkerSize(1);
0214     m_PVz_half2->SetMarkerColor(38);
0215     m_PVz->SetMarkerSize(1);
0216     m_PVz->SetMarkerColor(46);
0217     if (fabs(PVz_half1) <= 20.)
0218         m_PVz_half1->Draw("same");
0219     if (fabs(PVz_half2) <= 20.)
0220         m_PVz_half2->Draw("same");
0221     if (fabs(PVz) <= 20.)
0222         m_PVz->Draw("same");
0223 
0224     // TText *t_l1 = new TText(MVTXLength / 2, MVTXL1_min, "Layer 0");
0225     // TText *t_l2 = new TText(MVTXLength / 2, MVTXL2_min, "Layer 1");
0226     // t_l1->SetTextSize(0.04);
0227     // t_l2->SetTextSize(0.04);
0228     // t_l1->SetTextAlign(11);
0229     // t_l2->SetTextAlign(11);
0230     // t_l1->Draw("same");
0231     // t_l2->Draw("same");
0232 
0233     TLegend *leg = new TLegend(0.6, 0.8, 0.88, 0.92);
0234     // leg->SetNColumns(3);
0235     leg->SetTextSize(0.03);
0236     leg->SetFillStyle(0);
0237     leg->AddEntry(hM_clusZRho, "Clusters", "p");
0238     leg->AddEntry(m_PVz_half1, "Estimate of PV_{Reco}", "p");
0239     leg->AddEntry(m_PVz, "Final PV_{Reco}", "p");
0240     leg->Draw("same");
0241 
0242     TLatex *ltx = new TLatex();
0243     ltx->SetNDC();
0244     ltx->SetTextSize(0.045);
0245     ltx->SetTextAlign(kHAlignLeft + kVAlignBottom);
0246     if (isdata)
0247         ltx->DrawLatex(gPad->GetLeftMargin() + 0.05, 1 - gPad->GetTopMargin() - 0.08, "#it{#bf{sPHENIX}} Preliminary");
0248     else
0249         ltx->DrawLatex(gPad->GetLeftMargin() + 0.05, 1 - gPad->GetTopMargin() - 0.08, "#it{#bf{sPHENIX}} Simulation");
0250     ltx->SetTextSize(0.04);
0251     // ltx->SetTextAlign(kHAlignRight + kVAlignBottom);
0252     ltx->DrawLatex(gPad->GetLeftMargin() + 0.05, 1 - gPad->GetTopMargin() - 0.13, "Au+Au #sqrt{s_{NN}}=200 GeV");
0253     ltx->SetTextSize(0.035);
0254     ltx->SetTextAlign(kHAlignRight + kVAlignBottom);
0255     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, "2023-08-30");
0256 
0257     ltx->SetTextSize(0.025);
0258     ltx->SetTextAlign(kHAlignLeft + kVAlignBottom);
0259     ltx->DrawLatex(gPad->GetLeftMargin() + 0.05, gPad->GetBottomMargin() + 0.05, Form("N_{Clusters}^{1st + 2nd layer} = %d; Cluster positions assume ideal detector geometry", Nclus));
0260 
0261     c->SaveAs(plotname_tstr + ".pdf");
0262     c->SaveAs(plotname_tstr + ".png");
0263     c->Close();
0264 }
0265 
0266 // TCanvas *Canv_RecoVtxTklXY(vector<TLine *> v_line, vector<Hit *> hit_l1, vector<Hit *> hit_l2, float PVx, float PVy)
0267 void Draw_RecoVtxTklXY(TH2F *hM_clusXY_half1, TH2F *hM_clusXY_half2, vector<TLine *> vtl_half1, vector<TLine *> vtl_half2, float PVx_half1, float PVx_half2, float PVx, float PVy_half1, float PVy_half2, float PVy, int dZbin, bool isdata,
0268                        TString plotname_tstr)
0269 {
0270     SetsPhenixStyle();
0271     gStyle->SetPalette(kThermometer);
0272 
0273     int nTLines = vtl_half1.size() + vtl_half2.size();
0274     float tlinealpha = (nTLines > 50) ? 0.3 : 0.4;
0275     // merge histograms
0276     TH2F *hM_clusXY = (TH2F *)hM_clusXY_half1->Clone();
0277     hM_clusXY->Add(hM_clusXY_half2);
0278     int Nclus = hM_clusXY->GetEntries();
0279 
0280     TCanvas *c = new TCanvas("c", "c", 800, 780);
0281     c->cd();
0282     gPad->SetRightMargin(0.08);
0283     gPad->SetTopMargin(0.05);
0284     gPad->SetLeftMargin(0.15);
0285     gPad->SetBottomMargin(0.13);
0286     hM_clusXY->GetXaxis()->SetTitle("Position X [cm]");
0287     hM_clusXY->GetYaxis()->SetTitle("Position Y [cm]");
0288     hM_clusXY->GetXaxis()->SetTitleOffset(1.1);
0289     hM_clusXY->GetYaxis()->SetTitleOffset(1.3);
0290     hM_clusXY->SetMarkerSize(0.5);
0291     hM_clusXY->Draw("SCAT");
0292 
0293     for (size_t i = 0; i < vtl_half1.size(); i++)
0294     {
0295         vtl_half1[i]->SetLineColorAlpha(14, tlinealpha);
0296         vtl_half1[i]->Draw("same");
0297     }
0298     for (size_t i = 0; i < vtl_half2.size(); i++)
0299     {
0300         vtl_half2[i]->SetLineColorAlpha(14, tlinealpha);
0301         vtl_half2[i]->Draw("same");
0302     }
0303 
0304     hM_clusXY->Draw("SCATsame");
0305 
0306     TMarker *m_PV1 = new TMarker(PVx_half1, PVy_half1, 20);
0307     m_PV1->SetMarkerSize(1);
0308     m_PV1->SetMarkerColor(38);
0309     m_PV1->Draw("same");
0310     TMarker *m_PV2 = new TMarker(PVx_half2, PVy_half2, 20);
0311     m_PV2->SetMarkerSize(1);
0312     m_PV2->SetMarkerColor(38);
0313     m_PV2->Draw("same");
0314     TMarker *m_PV = new TMarker(PVx, PVy, 20);
0315     m_PV->SetMarkerSize(1);
0316     m_PV->SetMarkerColor(46);
0317     m_PV->Draw("same");
0318 
0319     TLine *MVTX_geo = new TLine();
0320     MVTX_geo->SetLineWidth(2);
0321     MVTX_geo->SetLineColorAlpha(kBlack, 0.3);
0322     vector<tuple<float, float, float, float>> MVTXstavepointsXY = MVTXStavePositionXY();
0323     for (auto &stave : MVTXstavepointsXY)
0324     {
0325         // exclude the staves which radius is larger than 3.625 cm
0326         float xmiddle = (get<0>(stave) + get<2>(stave)) / 2.;
0327         float ymiddle = (get<1>(stave) + get<3>(stave)) / 2.;
0328         float radius = sqrt(xmiddle * xmiddle + ymiddle * ymiddle);
0329         if (radius > 3.625)
0330             continue;
0331         MVTX_geo->DrawLine(get<0>(stave), get<1>(stave), get<2>(stave), get<3>(stave));
0332     }
0333 
0334     TLegend *leg = new TLegend(0.6, 0.8, 0.88, 0.92);
0335     // leg->SetNColumns(3);
0336     leg->SetTextSize(0.03);
0337     leg->SetFillStyle(0);
0338     leg->AddEntry(hM_clusXY, "Clusters", "p");
0339     leg->AddEntry(m_PV1, "Estimate of PV_{Reco}", "p");
0340     leg->AddEntry(m_PV, "Final PV_{Reco}", "p");
0341     leg->Draw("same");
0342 
0343     TLatex *ltx = new TLatex();
0344     ltx->SetNDC();
0345     ltx->SetTextSize(0.04);
0346     ltx->SetTextAlign(kHAlignLeft + kVAlignBottom);
0347     if (isdata)
0348         ltx->DrawLatex(gPad->GetLeftMargin() + 0.05, 1 - gPad->GetTopMargin() - 0.08, "#it{#bf{sPHENIX}} Preliminary");
0349     else
0350         ltx->DrawLatex(gPad->GetLeftMargin() + 0.05, 1 - gPad->GetTopMargin() - 0.08, "#it{#bf{sPHENIX}} Simulation");
0351     ltx->SetTextSize(0.035);
0352     // ltx->SetTextAlign(kHAlignRight + kVAlignBottom);
0353     ltx->DrawLatex(gPad->GetLeftMargin() + 0.05, 1 - gPad->GetTopMargin() - 0.13, "Au+Au #sqrt{s_{NN}}=200 GeV");
0354     ltx->SetTextSize(0.035);
0355     ltx->SetTextAlign(kHAlignRight + kVAlignBottom);
0356     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, "2023-08-30");
0357 
0358     ltx->SetTextSize(0.025);
0359     ltx->SetTextAlign(kHAlignLeft + kVAlignBottom);
0360     ltx->DrawLatex(gPad->GetLeftMargin() + 0.05, gPad->GetBottomMargin() + 0.09, Form("Tracklets of the chosen vertex cluster (|#Deltaz| < %g cm)", dZbin * 0.01));
0361     ltx->DrawLatex(gPad->GetLeftMargin() + 0.05, gPad->GetBottomMargin() + 0.05, Form("N_{Clusters}^{1st + 2nd layer} = %d; Cluster positions assume ideal detector geometry", Nclus));
0362 
0363     c->SaveAs(plotname_tstr + ".pdf");
0364     c->SaveAs(plotname_tstr + ".png");
0365     c->Close();
0366 }
0367 
0368 void Draw_1DHistVtxDemo(TH1F *hM, const char *Xtitle, const char *Ytitleunit, int ndivision, bool isdata, bool ispreliminary, const char *date, TString plotname_tstr)
0369 {
0370     SetsPhenixStyle();
0371     // Get bin size of the histogram
0372     double binwidth = hM->GetBinWidth(1);
0373 
0374     TCanvas *c = new TCanvas("c", "c", 800, 750);
0375     c->cd();
0376     gPad->SetRightMargin(0.06);
0377     gPad->SetTopMargin(0.05);
0378     gPad->SetLeftMargin(0.15);
0379     gPad->SetBottomMargin(0.13);
0380     hM->GetXaxis()->SetTitle(Xtitle);
0381     hM->GetYaxis()->SetTitle(Form("Counts / %g %s", binwidth, Ytitleunit));
0382     hM->GetXaxis()->SetTitleOffset(1.2);
0383     hM->GetYaxis()->SetTitleOffset(1.55);
0384     hM->GetYaxis()->SetRangeUser(0, hM->GetMaximum() * 1.35);
0385     hM->GetXaxis()->SetNdivisions(ndivision, kTRUE);
0386     hM->SetLineColor(1);
0387     hM->SetLineWidth(2);
0388     hM->Draw("hist");
0389     TLatex *ltx = new TLatex();
0390     ltx->SetNDC();
0391     ltx->SetTextSize(0.045);
0392     ltx->SetTextAlign(kHAlignLeft + kVAlignBottom);
0393     if (isdata)
0394         ltx->DrawLatex(gPad->GetLeftMargin() + 0.05, 1 - gPad->GetTopMargin() - 0.08, Form("#it{#bf{sPHENIX}} %s", (ispreliminary) ? "Preliminary" : "Internal"));
0395     else
0396         ltx->DrawLatex(gPad->GetLeftMargin() + 0.05, 1 - gPad->GetTopMargin() - 0.08, "#it{#bf{sPHENIX}} Simulation");
0397     ltx->SetTextSize(0.04);
0398     // ltx->SetTextAlign(kHAlignRight + kVAlignBottom);
0399     ltx->DrawLatex(gPad->GetLeftMargin() + 0.05, 1 - gPad->GetTopMargin() - 0.13, "Au+Au #sqrt{s_{NN}}=200 GeV");
0400     ltx->SetTextSize(0.035);
0401     ltx->SetTextAlign(kHAlignRight + kVAlignBottom);
0402     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, date);
0403     c->SaveAs(plotname_tstr + ".pdf");
0404     c->SaveAs(plotname_tstr + ".png");
0405     delete c;
0406 }
0407 
0408 void Draw_2DHistVtxDemo(TH2F *hM, const char *Xtitle, const char *Ytitle, int ndivision, bool isdata, bool ispreliminary, const char *date, TString plotname_tstr)
0409 {
0410     SetsPhenixStyle();
0411     gStyle->SetPalette(kThermometer);
0412 
0413     TCanvas *c = new TCanvas("c", "c", 800, 750);
0414     c->cd();
0415     gPad->SetRightMargin(0.13);
0416     gPad->SetTopMargin(0.05);
0417     gPad->SetLeftMargin(0.15);
0418     gPad->SetBottomMargin(0.13);
0419     hM->GetXaxis()->SetTitle(Xtitle);
0420     hM->GetYaxis()->SetTitle(Ytitle);
0421     hM->GetXaxis()->SetTitleOffset(1.2);
0422     hM->GetYaxis()->SetTitleOffset(1.45);
0423     hM->GetXaxis()->SetNdivisions(ndivision, kTRUE);
0424     hM->GetYaxis()->SetNdivisions(ndivision, kTRUE);
0425     hM->Draw("colz");
0426     TLatex *ltx = new TLatex();
0427     ltx->SetNDC();
0428     ltx->SetTextSize(0.045);
0429     ltx->SetTextAlign(kHAlignLeft + kVAlignBottom);
0430     if (isdata)
0431         ltx->DrawLatex(gPad->GetLeftMargin() + 0.05, 1 - gPad->GetTopMargin() - 0.08, Form("#it{#bf{sPHENIX}} %s", (ispreliminary) ? "Preliminary" : "Internal"));
0432     else
0433         ltx->DrawLatex(gPad->GetLeftMargin() + 0.05, 1 - gPad->GetTopMargin() - 0.08, "#it{#bf{sPHENIX}} Simulation");
0434     ltx->SetTextSize(0.04);
0435     // ltx->SetTextAlign(kHAlignRight + kVAlignBottom);
0436     ltx->DrawLatex(gPad->GetLeftMargin() + 0.05, 1 - gPad->GetTopMargin() - 0.13, "Au+Au #sqrt{s_{NN}}=200 GeV");
0437     ltx->SetTextSize(0.035);
0438     ltx->SetTextAlign(kHAlignRight + kVAlignBottom);
0439     ltx->DrawLatex(1 - gPad->GetRightMargin(), 1 - gPad->GetTopMargin() + 0.01, date);
0440     c->SaveAs(plotname_tstr + ".pdf");
0441     c->SaveAs(plotname_tstr + ".png");
0442     delete c;
0443 }
0444 
0445 // Calculate the DCA between the (extrapolated) tracklets and the PV
0446 TVector3 DCA_tklPV(vector<float> pos_clus1, vector<float> pos_clus2, vector<float> pos_PV)
0447 {
0448     TVector3 v_clus1(pos_clus1[0], pos_clus1[1], pos_clus1[2]); // B
0449     TVector3 v_clus2(pos_clus2[0], pos_clus2[1], pos_clus2[2]); // A
0450     TVector3 v_PV(pos_PV[0], pos_PV[1], pos_PV[2]);             // C
0451 
0452     TVector3 u = (v_clus1 - v_clus2).Unit(); // unit vextor from cluster 2 to cluster 1
0453     TVector3 AC = v_PV - v_clus2;
0454     TVector3 v_DCA = AC - (AC.Dot(u)) * u; // vector from cluster 2 to the closest point on the line of cluster 1 and cluster 2
0455 
0456     return v_DCA;
0457 }
0458 
0459 std::tuple<vector<float>, int, TH2F *, TH2F *, vector<TLine *>, vector<TLine *>, TH1F *, TH2F *, TH1F *, TH2F *, TH1F *, TH2F *, TH2F *> Tracklet3DPV(vector<Hit *> layer1, vector<Hit *> layer2, int dPhiCutbin, int dZCutbin)
0460 {
0461     // get a random number of integer
0462     TRandom3 *r = new TRandom3(0);
0463     int random_int = r->Integer(1E9);
0464     // tuples for (1) vector of PV estimates, (2) number of candidates in the vertex cluster, (3) TH2F of cluster positions in Z-Rho, (4) TH2F of cluster positions in X-Y, (5) vector of TLines in Z-Rho, (6) vector of TLines in X-Y, (7) TH1F of
0465     // difference between extrapolated z and PV z, (8) TH2F of difference between extrapolated x/y and PV x/y, (9) TH1F of z position of vertex candidates, (10) TH2F of x/y position of vertex candidates
0466     tuple<vector<float>, int, TH2F *, TH2F *, vector<TLine *>, vector<TLine *>, TH1F *, TH2F *, TH1F *, TH2F *, TH1F *, TH2F *, TH2F *> t;
0467     TH2F *hM_ZRho = new TH2F(Form("hM_ZRho_%d", random_int), Form("hM_ZRho_%d", random_int), 2000, -20, 20, 900, -0.6, 4.7);
0468     TH2F *hM_XY = new TH2F(Form("hM_XY_%d", random_int), Form("hM_XY_%d", random_int), 1100, -5.5, 5.5, 1100, -5.5, 5.5);
0469     TH1F *hM_Zdiff = new TH1F(Form("hM_Zdiff_%d", random_int), Form("hM_Zdiff_%d", random_int), 60, -0.06, 0.06);
0470     TH2F *hM_XYdiff = new TH2F(Form("hM_XYdiff_%d", random_int), Form("hM_XYdiff_%d", random_int), 100, -2, 2, 100, -2, 2);
0471     TH1F *hM_vtxcandz = new TH1F(Form("hM_vtxcandz_%d", random_int), Form("hM_vtxcandz_%d", random_int), 100, -70, 70);
0472     TH2F *hM_vtxcandxy = new TH2F(Form("hM_vtxcandxy_%d", random_int), Form("hM_vtxcandxy_%d", random_int), 100, -5, 5, 100, -5, 5);
0473     TH1F *hM_DCAxy = new TH1F(Form("hM_DCAxy_%d", random_int), Form("hM_DCAxy_%d", random_int), 100, 0, 0.5);
0474     TH2F *hM_DCAxy_DCAz = new TH2F(Form("hM_DCAxy_DCAz_%d", random_int), Form("hM_DCAxy_DCAz_%d", random_int), 100, 0, 0.5, 100, -0.05, 0.05);
0475     TH2F *hM_DCAx_DCAy = new TH2F(Form("hM_DCAx_DCAy_%d", random_int), Form("hM_DCAx_DCAy_%d", random_int), 100, -0.5, 0.5, 100, -0.5, 0.5);
0476 
0477     int Ncandidates_incluster = 0;
0478 
0479     for (size_t i = 0; i < layer1.size(); i++)
0480     {
0481         hM_ZRho->Fill(layer1[i]->posZ(), layer1[i]->rho());
0482     }
0483     for (size_t i = 0; i < layer2.size(); i++)
0484     {
0485         hM_ZRho->Fill(layer2[i]->posZ(), layer2[i]->rho());
0486     }
0487 
0488     for (size_t i = 0; i < layer1.size(); i++)
0489     {
0490         hM_XY->Fill(layer1[i]->posX(), layer1[i]->posY());
0491     }
0492     for (size_t i = 0; i < layer2.size(); i++)
0493     {
0494         hM_XY->Fill(layer2[i]->posX(), layer2[i]->posY());
0495     }
0496 
0497     float dPhi_cut = dPhiCutbin * 0.01;
0498     float dZ_cut = dZCutbin * 0.01;
0499 
0500     vector<float> vertex = {-999., -999., -999.};
0501 
0502     vector<MvtxClusComb> vertices;
0503     vertices.clear();
0504     vector<TLine *> v_tlkprojZ; // for debugging
0505     v_tlkprojZ.clear();
0506     vector<TLine *> v_tlkprojXY; // for debugging
0507     v_tlkprojXY.clear();
0508 
0509     for (size_t i = 0; i < layer1.size(); i++)
0510     {
0511         for (size_t j = 0; j < layer2.size(); j++)
0512         {
0513             if (fabs(deltaPhi(layer1[i]->Phi(), layer2[j]->Phi())) > dPhi_cut)
0514                 continue;
0515 
0516             float z = layer1[i]->posZ() - (layer2[j]->posZ() - layer1[i]->posZ()) / (layer2[j]->rho() - layer1[i]->rho()) * layer1[i]->rho();
0517             vertices.emplace_back(layer1[i]->posX(), layer1[i]->posY(), layer1[i]->posZ(), layer2[j]->posX(), layer2[j]->posY(), layer2[j]->posZ(), z);
0518 
0519             // find the equation connecting two hits in the X-Y plane
0520             // y = mx + b
0521             float m = (layer2[j]->posY() - layer1[i]->posY()) / (layer2[j]->posX() - layer1[i]->posX());
0522             float b = layer1[i]->posY() - m * layer1[i]->posX();
0523             // calculate the distance of closest approach from the line to (0,0)
0524             float d = fabs(b) / sqrt(1 + m * m);
0525 
0526             if (d > 0.3)
0527                 continue;
0528 
0529             // for debugging - print out the z and rho of the clusters and the extrpolated z position
0530             // cout << "layer1[i]->posZ() = " << layer1[i]->posZ() << ", layer1[i]->rho() = " << layer1[i]->rho() << ", layer2[j]->posZ() = " << layer2[j]->posZ() << ", layer2[j]->rho() = " << layer2[j]->rho() << ", z = " << z << endl;
0531 
0532             hM_vtxcandz->Fill(z);
0533 
0534             if (z <= 30. && z >= -30.)
0535             {
0536                 TLine *l = new TLine(layer2[j]->posZ(), layer2[j]->rho(), z, 0.0);
0537                 v_tlkprojZ.push_back(l);
0538             }
0539         }
0540     }
0541 
0542     /* vertex clustering */
0543     if (vertices.size())
0544     {
0545         // Sort the vertices based on the extrapolated z position
0546         sort(vertices.begin(), vertices.end(), [](const MvtxClusComb &a, const MvtxClusComb &b) -> bool { return a.extz < b.extz; });
0547         // sort(vertices.begin(), vertices.end());
0548 
0549         // vertex candidate clusters
0550         std::vector<Cluster> clusters;
0551         std::size_t z = 0, y = 0;
0552         for (; z < vertices.size(); ++z)
0553         {
0554             for (; y < vertices.size() && vertices[y].extz - vertices[z].extz < dZ_cut; ++y)
0555                 ;
0556             clusters.emplace_back(z, y - z); // index, nz
0557         }
0558 
0559         std::sort(clusters.begin(), clusters.end(), [](const Cluster &a, const Cluster &b) -> bool { return a.nz > b.nz; });
0560 
0561         std::vector<Vertex> candidates;
0562         uint32_t maxnz = clusters[0].nz;
0563 
0564         for (auto &cluster : clusters)
0565         {
0566             uint32_t nz = cluster.nz;
0567             if (nz + NZGAP < maxnz)
0568             {
0569                 break;
0570             }
0571 
0572             uint32_t index = cluster.index;
0573 
0574             float vz = 0;
0575             float vz2 = 0;
0576             float sigma2 = 0;
0577             vector<float> hitl1x, hitl1y, hitl1z, hitl2x, hitl2y, hitl2z, candidatez;
0578 
0579             for (uint32_t y = 0; y < nz; ++y)
0580             {
0581                 vz += vertices[index + y].extz;
0582                 vz2 += vertices[index + y].extz * vertices[index + y].extz;
0583                 candidatez.push_back(vertices[index + y].extz);
0584                 hitl1x.push_back(vertices[index + y].x1);
0585                 hitl1y.push_back(vertices[index + y].y1);
0586                 hitl1z.push_back(vertices[index + y].z1);
0587                 hitl2x.push_back(vertices[index + y].x2);
0588                 hitl2y.push_back(vertices[index + y].y2);
0589                 hitl2z.push_back(vertices[index + y].z2);
0590             }
0591 
0592             vz /= nz;
0593             sigma2 = vz2 / nz - (vz * vz);
0594 
0595             candidates.emplace_back(vz, sigma2, hitl1x, hitl1y, hitl1z, hitl2x, hitl2y, hitl2z, candidatez);
0596         }
0597 
0598         // PV z position
0599         vertex[2] = std::min_element(candidates.begin(), candidates.end(), [](const Vertex &a, const Vertex &b) -> bool { return a.sigma2 < b.sigma2; })->vz;
0600         vector<float> vtxcandz = std::min_element(candidates.begin(), candidates.end(), [](const Vertex &a, const Vertex &b) -> bool { return a.sigma2 < b.sigma2; })->candz;
0601         vector<float> hit1x = std::min_element(candidates.begin(), candidates.end(), [](const Vertex &a, const Vertex &b) -> bool { return a.sigma2 < b.sigma2; })->hitx1;
0602         vector<float> hit1y = std::min_element(candidates.begin(), candidates.end(), [](const Vertex &a, const Vertex &b) -> bool { return a.sigma2 < b.sigma2; })->hity1;
0603         vector<float> hit1z = std::min_element(candidates.begin(), candidates.end(), [](const Vertex &a, const Vertex &b) -> bool { return a.sigma2 < b.sigma2; })->hitz1;
0604         vector<float> hit2x = std::min_element(candidates.begin(), candidates.end(), [](const Vertex &a, const Vertex &b) -> bool { return a.sigma2 < b.sigma2; })->hitx2;
0605         vector<float> hit2y = std::min_element(candidates.begin(), candidates.end(), [](const Vertex &a, const Vertex &b) -> bool { return a.sigma2 < b.sigma2; })->hity2;
0606         vector<float> hit2z = std::min_element(candidates.begin(), candidates.end(), [](const Vertex &a, const Vertex &b) -> bool { return a.sigma2 < b.sigma2; })->hitz2;
0607         vector<float> vtxcandx, vtxcandy;
0608 
0609         for (size_t i = 0; i < hit1x.size(); i++)
0610         {
0611             float extpx = (vertex[2] - hit1z[i]) / (hit1z[i] - hit2z[i]) * (hit1x[i] - hit2x[i]) + hit1x[i];
0612             float extpy = (vertex[2] - hit1z[i]) / (hit1z[i] - hit2z[i]) * (hit1y[i] - hit2y[i]) + hit1y[i];
0613             vtxcandx.push_back(extpx);
0614             vtxcandy.push_back(extpy);
0615 
0616             hM_vtxcandxy->Fill(extpx, extpy);
0617 
0618             TLine *l_tlkprojXY = new TLine(hit2x[i], hit2y[i], extpx, extpy);
0619             v_tlkprojXY.push_back(l_tlkprojXY);
0620         }
0621 
0622         // copy vector vtxcandx and vtxcandy to a new vector for sorting
0623         vector<float> vtxcandx_sorted = vtxcandx;
0624         vector<float> vtxcandy_sorted = vtxcandy;
0625         // find the median of the extrapolated x and y positions as the PV_x and PV_y
0626         std::sort(vtxcandx_sorted.begin(), vtxcandx_sorted.end());
0627         std::sort(vtxcandy_sorted.begin(), vtxcandy_sorted.end());
0628         vertex[0] = (vtxcandx_sorted.size() % 2 == 0) ? (vtxcandx_sorted[vtxcandx_sorted.size() / 2 - 1] + vtxcandx_sorted[vtxcandx_sorted.size() / 2]) / 2 : vtxcandx_sorted[vtxcandx_sorted.size() / 2];
0629         vertex[1] = (vtxcandy_sorted.size() % 2 == 0) ? (vtxcandy_sorted[vtxcandy_sorted.size() / 2 - 1] + vtxcandy_sorted[vtxcandy_sorted.size() / 2]) / 2 : vtxcandy_sorted[vtxcandy_sorted.size() / 2];
0630         Ncandidates_incluster = vtxcandx.size();
0631 
0632         // fill the histograms
0633         for (size_t i = 0; i < vtxcandx.size(); i++)
0634         {
0635             hM_Zdiff->Fill(vtxcandz[i] - vertex[2]);
0636             hM_XYdiff->Fill(vtxcandx[i] - vertex[0], vtxcandy[i] - vertex[1]);
0637             // Calculate the DCA between the (extrapolated) tracklets and the PV in the vertex cluster
0638             TVector3 v_DCA = DCA_tklPV({hit1x[i], hit1y[i], hit1z[i]}, {hit2x[i], hit2y[i], hit2z[i]}, vertex);
0639             hM_DCAxy->Fill(v_DCA.XYvector().Mod());
0640             hM_DCAxy_DCAz->Fill(v_DCA.XYvector().Mod(), v_DCA.Z());
0641             hM_DCAx_DCAy->Fill(v_DCA.X(), v_DCA.Y());
0642         }
0643 
0644         // clean up
0645         CleanVec(clusters);
0646         CleanVec(candidates);
0647         CleanVec(vtxcandz);
0648         CleanVec(hit1x);
0649         CleanVec(hit1y);
0650         CleanVec(hit1z);
0651         CleanVec(hit2x);
0652         CleanVec(hit2y);
0653         CleanVec(hit2z);
0654         CleanVec(vtxcandx);
0655         CleanVec(vtxcandy);
0656         CleanVec(vtxcandx_sorted);
0657         CleanVec(vtxcandy_sorted);
0658     }
0659 
0660     CleanVec(vertices);
0661 
0662     // setup the tuple
0663     get<0>(t) = vertex;
0664     get<1>(t) = Ncandidates_incluster;
0665     get<2>(t) = hM_ZRho;
0666     get<3>(t) = hM_XY;
0667     get<4>(t) = v_tlkprojZ;
0668     get<5>(t) = v_tlkprojXY;
0669     get<6>(t) = hM_Zdiff;
0670     get<7>(t) = hM_XYdiff;
0671     get<8>(t) = hM_vtxcandz;
0672     get<9>(t) = hM_vtxcandxy;
0673     get<10>(t) = hM_DCAxy;
0674     get<11>(t) = hM_DCAxy_DCAz;
0675     get<12>(t) = hM_DCAx_DCAy;
0676 
0677     return t;
0678 }
0679 
0680 void SetVtxMinitree(TTree *outTree, VtxData &vtxdata)
0681 {
0682     outTree->Branch("event", &vtxdata.event);
0683     outTree->Branch("NhitsLayer1", &vtxdata.NhitsLayer1);
0684     if (!vtxdata.isdata)
0685     {
0686         outTree->Branch("NTruthVtx", &vtxdata.NTruthVtx);
0687         outTree->Branch("TruthPV_trig_x", &vtxdata.TruthPV_trig_x);
0688         outTree->Branch("TruthPV_trig_y", &vtxdata.TruthPV_trig_y);
0689         outTree->Branch("TruthPV_trig_z", &vtxdata.TruthPV_trig_z);
0690         outTree->Branch("TruthPV_mostNpart_x", &vtxdata.TruthPV_mostNpart_x);
0691         outTree->Branch("TruthPV_mostNpart_y", &vtxdata.TruthPV_mostNpart_y);
0692         outTree->Branch("TruthPV_mostNpart_z", &vtxdata.TruthPV_mostNpart_z);
0693         outTree->Branch("Centrality_bimp", &vtxdata.Centrality_bimp);
0694         outTree->Branch("Centrality_impactparam", &vtxdata.Centrality_impactparam);
0695     }
0696     outTree->Branch("Centrality_mbd", &vtxdata.Centrality_mbd);
0697     outTree->Branch("Centrality_mbdquantity", &vtxdata.Centrality_mbdquantity);
0698     outTree->Branch("PV_x_halves1", &vtxdata.PV_x_halves1);
0699     outTree->Branch("PV_y_halves1", &vtxdata.PV_y_halves1);
0700     outTree->Branch("PV_z_halves1", &vtxdata.PV_z_halves1);
0701     outTree->Branch("PV_x_halves2", &vtxdata.PV_x_halves2);
0702     outTree->Branch("PV_y_halves2", &vtxdata.PV_y_halves2);
0703     outTree->Branch("PV_z_halves2", &vtxdata.PV_z_halves2);
0704 }
0705 
0706 std::map<int, float> EvtVtx_map_tklcluster(const char *vtxfname)
0707 {
0708     std::map<int, float> EvtVtx_map;
0709 
0710     TFile *f = new TFile(vtxfname, "READ");
0711     TTree *t = (TTree *)f->Get("minitree");
0712     int event;
0713     float PV_z_halves1, PV_z_halves2;
0714     t->SetBranchAddress("event", &event);
0715     t->SetBranchAddress("PV_z_halves1", &PV_z_halves1);
0716     t->SetBranchAddress("PV_z_halves2", &PV_z_halves2);
0717     for (int ev = 0; ev < t->GetEntriesFast(); ev++)
0718     {
0719         t->GetEntry(ev);
0720         // cout << ev << " " << event << " " << PV_z << endl;
0721         EvtVtx_map[event] = (PV_z_halves1 + PV_z_halves2) / 2.;
0722     }
0723 
0724     return EvtVtx_map;
0725 }
0726 
0727 #endif
[ 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!