MVTXTracklet/src/Tracklet_RecoVtx.cxx

0001 #include <fstream>
0002 #include <iostream>
0003 #include <sstream>
0004 #include <string>
0005 #include <vector>
0006
0007 #include <TCanvas.h>
0008 #include <TF1.h>
0009 #include <TFile.h>
0010 #include <TH1F.h>
0011 #include <TH2F.h>
0012 #include <TLegend.h>
0013 #include <TLine.h>
0014 #include <TMath.h>
0015 #include <TObjString.h>
0016 #include <TProfile.h>
0017 #include <TRandom3.h>
0018 #include <TText.h>
0019 #include <TTree.h>
0020 #include <TTreeIndex.h>
0021
0022 #include "Hit.h"
0023 #include "Utilities.h"
0024 #include "Vertex.h"
0025 // #include "misalignment.h"
0026
0027 float precision(float f, int places)
0028 {
0029     float n = std::pow(10.0f, places);
0030     return std::round(f * n) / n;
0031 }
0032
0033 int main(int argc, char *argv[])
0034 {
0035     bool IsData = (TString(argv[1]).Atoi() == 1) ? true : false;
0036     int dPhiCutbin = TString(argv[2]).Atoi();
0037     int dZCutbin = TString(argv[3]).Atoi();
0038     // float gap_north = TString(argv[4]).Atof();  // Nominal: 2.9mm
0039     // float gap_upper = TString(argv[5]).Atof();  // Nominal: (gap_north / 2.) mm
0040     // float cent_shift = TString(argv[6]).Atof(); // Nominal: 0.0mm
0041
0042     // TString gap_north_str = TString(argv[4]).Replace(TString(argv[4]).Index("."), 1, "p");
0043     // TString gap_upper_str = TString(argv[5]).Replace(TString(argv[5]).Index("."), 1, "p");
0044     // TString cent_shift_str = TString(argv[6]).Replace(TString(argv[6]).Index("."), 1, "p");
0045
0046     // float dalpha_num = (gap_north-1.5)/2.;
0047
0048     float beamspot_x = 0.544991;
0049     float beamspot_y = -0.00600619;
0050
0051     SetsPhenixStyle();
0052
0053     TString infilename = "";
0054     TString outdirname = "";
0055     TString outfilename = "";
0056     TString demoplotname = "";
0057     if (IsData)
0058     {
0059         cout << "Running on data" << endl;
0060
0061         int runnumber = 45199;
0062         infilename = "/sphenix/user/hjheng/sPHENIXRepo/mvtx_standalone_cluster/macros/outputClusters_000"+ TString::Format("%d", runnumber) + ".root";
0063         outdirname = "/sphenix/user/hjheng/sPHENIXRepo/analysis/MVTXTracklet/MVTX_Run2024_000" + TString::Format("%d", runnumber);
0064         system("mkdir -p " + outdirname);
0065         outfilename = outdirname + "/Minitree_TrackletRecoVtx_000" + TString::Format("%d", runnumber) + ".root";
0066         demoplotname = "./plot/RecoPV_demo/RecoPV_Data/Run000" + TString::Format("%d", runnumber);
0067         system("mkdir -p " + demoplotname);
0068     }
0069     else
0070     {
0071         cout << "Running on MC" << endl;
0072
0073         infilename = "/sphenix/user/hjheng/TrackletAna/data/MVTXRecoClusters/NoPileup_Nevt500_ana325private_singleEvtDst/MVTXRecoClusters.root";
0074         outfilename = "";
0075         demoplotname = "./plot/RecoPV_demo/RecoPV_Sim/AuAuHijing_FieldOn";
0076     }
0077
0078     vector<Hit *> MVTXlayer1_halves1, MVTXlayer1_halves2, MVTXlayer2_halves1, MVTXlayer2_halves2;
0079
0080     VtxData vtxdata = {};
0081
0082     TFile *f = new TFile(infilename, "READ");
0083     TTree *t = (TTree *)f->Get("Hits");
0084     t->BuildIndex("event"); // Reference: https://root-forum.cern.ch/t/sort-ttree-entries/13138
0085     TTreeIndex *index = (TTreeIndex *)t->GetTreeIndex();
0086     int event, NTruthVtx;
0087     float TruthPV_trig_x, TruthPV_trig_y, TruthPV_trig_z, TruthPV_mostNpart_x, TruthPV_mostNpart_y, TruthPV_mostNpart_z;
0088     float centrality_bimp, centrality_impactparam, centrality_mbd, centrality_mbdquantity;
0089     vector<int> *ClusLayer = 0;
0090     vector<float> *ClusX = 0, *ClusY = 0, *ClusZ = 0;
0091     t->SetBranchAddress("event", &event);
0092     if (!IsData)
0093     {
0094         t->SetBranchAddress("NTruthVtx", &NTruthVtx);
0095         t->SetBranchAddress("TruthPV_trig_x", &TruthPV_trig_x);
0096         t->SetBranchAddress("TruthPV_trig_y", &TruthPV_trig_y);
0097         t->SetBranchAddress("TruthPV_trig_z", &TruthPV_trig_z);
0098         t->SetBranchAddress("TruthPV_mostNpart_x", &TruthPV_mostNpart_x);
0099         t->SetBranchAddress("TruthPV_mostNpart_y", &TruthPV_mostNpart_y);
0100         t->SetBranchAddress("TruthPV_mostNpart_z", &TruthPV_mostNpart_z);
0101         // t->SetBranchAddress("centrality_bimp", &centrality_bimp);
0102         // t->SetBranchAddress("centrality_impactparam", &centrality_impactparam);
0103     }
0104     // t->SetBranchAddress("centrality_mbd", &centrality_mbd);
0105     // t->SetBranchAddress("centrality_mbdquantity", &centrality_mbdquantity);
0106     t->SetBranchAddress("clusLayer", &ClusLayer);
0107     t->SetBranchAddress("globalX", &ClusX);
0108     t->SetBranchAddress("globalY", &ClusY);
0109     t->SetBranchAddress("globalZ", &ClusZ);
0110
0111     TFile *outfile = new TFile(outfilename, "RECREATE");
0112     TTree *minitree = new TTree("minitree", "Minitree of reconstructed vertices");
0113     SetVtxMinitree(minitree, vtxdata);
0114
0115     for (int i = 0; i < index->GetN(); i++)
0116     // for (int i = 0; i < 1; i++)
0117     // for (Long64_t ev = 0; ev < t->GetEntriesFast(); ev++)
0118     {
0119         Long64_t local = t->LoadTree(index->GetIndex()[i]);
0120         t->GetEntry(local);
0121         // cout << "event = " << event << " local = " << local << endl;
0122         cout << "event=" << event << " has a total of " << ClusLayer->size() << " clusters" << endl;
0123
0124         CleanVec(MVTXlayer1_halves1);
0125         CleanVec(MVTXlayer1_halves2);
0126         CleanVec(MVTXlayer2_halves1);
0127         CleanVec(MVTXlayer2_halves2);
0128
0129         for (size_t ihit = 0; ihit < ClusLayer->size(); ihit++)
0130         {
0131             // Misalignment study for MVTX
0132             float x0 = ClusX->at(ihit);
0133             float y0 = ClusY->at(ihit);
0134             float z0 = ClusZ->at(ihit);
0135             vector<float> cpos = {x0, y0, z0};
0136             // UpdatePos_GapTwoHalves(cpos, gap_north, cent_shift, gap_upper);
0137
0138             if (ClusLayer->at(ihit) == 0)
0139             {
0140                 Hit *hit = new Hit(cpos[0], cpos[1], cpos[2], beamspot_x, beamspot_y, 0., 0);
0141                 // Hit *hit = new Hit(ClusX->at(ihit), ClusY->at(ihit), ClusZ->at(ihit), 0., 0., 0., 0);
0142                 if (x0 > 0)
0143                 {
0144                     MVTXlayer1_halves1.push_back(hit);
0145                 }
0146
0147                 else
0148                     MVTXlayer1_halves2.push_back(hit);
0149             }
0150             else if (ClusLayer->at(ihit) == 1)
0151             {
0152                 Hit *hit = new Hit(cpos[0], cpos[1], cpos[2], beamspot_x, beamspot_y, 0., 1);
0153                 // Hit *hit = new Hit(ClusX->at(ihit), ClusY->at(ihit), ClusZ->at(ihit), 0., 0., 0., 1);
0154                 if (x0 > 0)
0155                     MVTXlayer2_halves1.push_back(hit);
0156                 else
0157                     MVTXlayer2_halves2.push_back(hit);
0158             }
0159             else
0160                 continue;
0161         }
0162
0163         // cout << "# of clusters in Layer 1 = " << MVTXlayer1.size() << ", Layer 2 = " << MVTXlayer2.size() << endl;
0164         cout << "# of clusters in Layer 1 = " << MVTXlayer1_halves1.size() + MVTXlayer1_halves2.size() << ", Layer 2 = " << MVTXlayer2_halves1.size() + MVTXlayer2_halves2.size() << endl;
0165
0166         // float PV_z = TrackletPV_cluster(event, MVTXlayer1, MVTXlayer2, dPhiCutbin, dZCutbin, false);
0167         // Tracklet3DPV(int evt, vector<Hit *> layer1, vector<Hit *> layer2, int dPhiCutbin, int dZCutbin, const char *plotname)
0168         // vector<float> PV_halves1 = Tracklet3DPV(event, MVTXlayer1_halves1, MVTXlayer2_halves1, dPhiCutbin, dZCutbin, Form("./plot/RecoPV_Data/Run20124_Evt%d_RecoPV_halve1", event));
0169         // vector<float> PV_halves2 = Tracklet3DPV(event, MVTXlayer1_halves2, MVTXlayer2_halves2, dPhiCutbin, dZCutbin, Form("./plot/RecoPV_Data/Run20124_Evt%d_RecoPV_halve2", event));
0170         // tuples for (1) vector of PV estimates, (2) TH2F of cluster positions in Z-Rho, (3) TH2F of cluster positions in X-Y, (4) vector of TLines in Z-Rho, (5) vector of TLines in X-Y, (6) TH1F of difference between extrapolated z and PV z, (7) TH2F of difference between extrapolated x/y and PV x/y
0171         tuple<vector<float>, int, TH2F *, TH2F *, vector<TLine *>, vector<TLine *>, TH1F *, TH2F *, TH1F *, TH2F *, TH1F *, TH2F *, TH2F *> tuple_PVest1 = Tracklet3DPV(MVTXlayer1_halves1, MVTXlayer2_halves1, dPhiCutbin, dZCutbin);
0172         tuple<vector<float>, int, TH2F *, TH2F *, vector<TLine *>, vector<TLine *>, TH1F *, TH2F *, TH1F *, TH2F *, TH1F *, TH2F *, TH2F *> tuple_PVest2 = Tracklet3DPV(MVTXlayer1_halves2, MVTXlayer2_halves2, dPhiCutbin, dZCutbin);
0173         vector<float> PV_halves1 = get<0>(tuple_PVest1);
0174         vector<float> PV_halves2 = get<0>(tuple_PVest2);
0175         int Nvtxcandidates_halves1 = get<1>(tuple_PVest1);
0176         int Nvtxcandidates_halves2 = get<1>(tuple_PVest2);
0177         TH2F *hM_clusterzrho_halves1 = get<2>(tuple_PVest1);
0178         TH2F *hM_clusterzrho_halves2 = get<2>(tuple_PVest2);
0179         TH2F *hM_clusterxy_halves1 = get<3>(tuple_PVest1);
0180         TH2F *hM_clusterxy_halves2 = get<3>(tuple_PVest2);
0181         vector<TLine *> tl_zrho_halves1 = get<4>(tuple_PVest1);
0182         vector<TLine *> tl_zrho_halves2 = get<4>(tuple_PVest2);
0183         vector<TLine *> tl_xy_halves1 = get<5>(tuple_PVest1);
0184         vector<TLine *> tl_xy_halves2 = get<5>(tuple_PVest2);
0185         TH1F *hM_diffExptzPVz_halves1 = get<6>(tuple_PVest1);
0186         TH1F *hM_diffExptzPVz_halves2 = get<6>(tuple_PVest2);
0187         TH2F *hM_diffExptxyPVxy_halves1 = get<7>(tuple_PVest1);
0188         TH2F *hM_diffExptxyPVxy_halves2 = get<7>(tuple_PVest2);
0189         TH1F *hM_VtxCandZ_halves1 = get<8>(tuple_PVest1);
0190         TH1F *hM_VtxCandZ_halves2 = get<8>(tuple_PVest2);
0191         TH2F *hM_VtxCandXY_halves1 = get<9>(tuple_PVest1);
0192         TH2F *hM_VtxCandXY_halves2 = get<9>(tuple_PVest2);
0193         TH1F *hM_DCAxy_VtxClus_halves1 = get<10>(tuple_PVest1);
0194         TH1F *hM_DCAxy_VtxClus_halves2 = get<10>(tuple_PVest2);
0195         TH2F *hM_DCAxy_DCAz_VtxClus_halves1 = get<11>(tuple_PVest1);
0196         TH2F *hM_DCAxy_DCAz_VtxClus_halves2 = get<11>(tuple_PVest2);
0197         TH2F *hM_DCAx_DCAy_VtxClus_halves1 = get<12>(tuple_PVest1);
0198         TH2F *hM_DCAx_DCAy_VtxClus_halves2 = get<12>(tuple_PVest2);
0199
0200         // weighted average of the PV z position from the 1st and 2nd halves
0201         float PV_x = (PV_halves1[0] * Nvtxcandidates_halves1 + PV_halves2[0] * Nvtxcandidates_halves2) / (Nvtxcandidates_halves1 + Nvtxcandidates_halves2);
0202         float PV_y = (PV_halves1[1] * Nvtxcandidates_halves1 + PV_halves2[1] * Nvtxcandidates_halves2) / (Nvtxcandidates_halves1 + Nvtxcandidates_halves2);
0203         float PV_z = (PV_halves1[2] * Nvtxcandidates_halves1 + PV_halves2[2] * Nvtxcandidates_halves2) / (Nvtxcandidates_halves1 + Nvtxcandidates_halves2);
0204
0205         if (MVTXlayer1_halves1.size() + MVTXlayer1_halves2.size() + MVTXlayer2_halves1.size() + MVTXlayer2_halves2.size() < Ncluster_toplot && i <= 60)
0206         {
0207             Draw_RecoVtxZTklZRho(hM_clusterzrho_halves1, hM_clusterzrho_halves2, tl_zrho_halves1, tl_zrho_halves2, PV_halves1[2], PV_halves2[2], PV_z, IsData, Form("%s/Evt%d_RecoPV_ZRho", demoplotname.Data(), event));
0208             Draw_RecoVtxTklXY(hM_clusterxy_halves1, hM_clusterxy_halves2, tl_xy_halves1, tl_xy_halves2, PV_halves1[0], PV_halves2[0], PV_x, PV_halves1[1], PV_halves2[1], PV_y, dZCutbin, IsData, Form("%s/Evt%d_RecoPV_XY", demoplotname.Data(), event));
0209             // Draw_diffExptzPVz(hM_diffExptzPVz_halves1, hM_diffExptzPVz_halves2, IsData, Form("%s_Evt%d_diffExptzPVz", demoplotname.Data(), event));
0210             // Draw_diffExptxyPVxy(hM_diffExptxyPVxy_halves1, hM_diffExptxyPVxy_halves2, IsData, Form("%s_Evt%d_diffExptxyPVxy", demoplotname.Data(), event));
0211         }
0212
0213         // merge the histograms
0214         hM_clusterzrho_halves1->Add(hM_clusterzrho_halves2);
0215         Draw_1DHistVtxDemo(hM_diffExptzPVz_halves1, "v_{Z}^{Candidate} - PV_{Z}^{Estimate} [cm]", "cm", 510, IsData, true, "2023-08-30", Form("%s/Evt%d_diffExptzPVz", demoplotname.Data(), event));
0216         hM_diffExptxyPVxy_halves1->Add(hM_diffExptxyPVxy_halves2);
0217         Draw_2DHistVtxDemo(hM_diffExptxyPVxy_halves1, "v_{X}^{Candidate} - PV_{X}^{Estimate} [cm]", "v_{Y}^{Candidate} - PV_{Y}^{Estimate} [cm]", 510, IsData, true, "2023-08-30", Form("%s/vt%d_diffExptxyPVxy", demoplotname.Data(), event));
0218
0219         hM_VtxCandZ_halves1->Add(hM_VtxCandZ_halves2);
0220         Draw_1DHistVtxDemo(hM_VtxCandZ_halves1, "v_{Z}^{Candidate} [cm]", "cm", 510, IsData, true, "2023-08-30", Form("%s/Evt%d_VtxCandZall", demoplotname.Data(), event));
0221         hM_VtxCandXY_halves1->Add(hM_VtxCandXY_halves2);
0222         Draw_2DHistVtxDemo(hM_VtxCandXY_halves1, "v_{X}^{Candidate} [cm]", "v_{Y}^{Candidate} [cm]", 510, IsData, true, "2023-08-30", Form("%s/Evt%d_VtxCandXYcluster", demoplotname.Data(), event));
0223
0224         hM_DCAxy_VtxClus_halves1->Add(hM_DCAxy_VtxClus_halves2);
0225         Draw_1DHistVtxDemo(hM_DCAxy_VtxClus_halves1, "DCA_{xy} [cm]", "cm", 505, IsData, false, "", Form("%s/Evt%d_DCAxy", demoplotname.Data(), event));
0226         hM_DCAxy_DCAz_VtxClus_halves1->Add(hM_DCAxy_DCAz_VtxClus_halves2);
0227         Draw_2DHistVtxDemo(hM_DCAxy_DCAz_VtxClus_halves1, "DCA_{xy} [cm]", "DCA_{z} [cm]", 505, IsData, false, "", Form("%s/Evt%d_DCAxyDCAz", demoplotname.Data(), event));
0228         hM_DCAx_DCAy_VtxClus_halves1->Add(hM_DCAx_DCAy_VtxClus_halves2);
0229         Draw_2DHistVtxDemo(hM_DCAx_DCAy_VtxClus_halves1, "DCA_{x} [cm]", "DCA_{y} [cm]", 505, IsData, false, "", Form("%s/Evt%d_DCAxDCAy", demoplotname.Data(), event));
0230
0231         CleanVec(tl_zrho_halves1);
0232         CleanVec(tl_zrho_halves2);
0233         CleanVec(tl_xy_halves1);
0234         CleanVec(tl_xy_halves2);
0235         delete hM_clusterzrho_halves1;
0236         delete hM_clusterzrho_halves2;
0237         delete hM_clusterxy_halves1;
0238         delete hM_clusterxy_halves2;
0239         delete hM_diffExptzPVz_halves1;
0240         delete hM_diffExptzPVz_halves2;
0241         delete hM_diffExptxyPVxy_halves1;
0242         delete hM_diffExptxyPVxy_halves2;
0243         delete hM_VtxCandZ_halves1;
0244         delete hM_VtxCandZ_halves2;
0245         delete hM_VtxCandXY_halves1;
0246         delete hM_VtxCandXY_halves2;
0247         delete hM_DCAxy_VtxClus_halves1;
0248         delete hM_DCAxy_VtxClus_halves2;
0249         delete hM_DCAxy_DCAz_VtxClus_halves1;
0250         delete hM_DCAxy_DCAz_VtxClus_halves2;
0251         delete hM_DCAx_DCAy_VtxClus_halves1;
0252         delete hM_DCAx_DCAy_VtxClus_halves2;
0253
0254         if (!IsData)
0255         {
0256             cout << "NTruthVtx = " << NTruthVtx << endl
0257                  << "Truth primary (most N particles) Z position = " << TruthPV_trig_z << endl
0258                  << "Recontructed primary vertex Z position from 1st halve= " << PV_halves1[2] << endl
0259                  << "Recontructed primary vertex Z position from 2nd halve= " << PV_halves2[2] << endl
0260                  << "Average of the Z positions from 1st and 2nd halve= " << (PV_halves1[2] + PV_halves2[2]) / 2. << endl
0261                  << "(Truth-Reco) = " << TruthPV_trig_z - (PV_halves1[2] + PV_halves2[2]) / 2. << endl;
0262
0263             if (NTruthVtx == 1)
0264             {
0265                 if (TruthPV_trig_z != TruthPV_mostNpart_z)
0266                     cout << "For events with only 1 truth vertex, (TruthPV_trig_z,TruthPV_mostNpart_z) = (" << TruthPV_trig_z << "," << TruthPV_mostNpart_z << "). Check!" << endl;
0267             }
0268         }
0269         else
0270         {
0271             cout << "Reconstructed primary vertex X position from 1st halve= " << PV_halves1[0] << endl << "Reconstructed primary vertex X position from 2nd halve= " << PV_halves2[0] << endl << "Average of the X positions from 1st and 2nd halve= " << (PV_halves1[0] + PV_halves2[0]) / 2. << endl;
0272             cout << "Reconstructed primary vertex Y position from 1st halve= " << PV_halves1[1] << endl << "Reconstructed primary vertex Y position from 2nd halve= " << PV_halves2[1] << endl << "Average of the Y positions from 1st and 2nd halve= " << (PV_halves1[1] + PV_halves2[1]) / 2. << endl;
0273             cout << "Recontructed primary vertex Z position from 1st halve= " << PV_halves1[2] << endl << "Recontructed primary vertex Z position from 2nd halve= " << PV_halves2[2] << endl << "Average of the Z positions from 1st and 2nd halve= " << (PV_halves1[2] + PV_halves2[2]) / 2. << endl;
0274         }
0275
0276         vtxdata.event = event;
0277         if (!IsData)
0278         {
0279             vtxdata.NTruthVtx = NTruthVtx;
0280             vtxdata.TruthPV_trig_x = TruthPV_trig_x;
0281             vtxdata.TruthPV_trig_y = TruthPV_trig_y;
0282             vtxdata.TruthPV_trig_z = TruthPV_trig_z;
0283             vtxdata.TruthPV_mostNpart_x = TruthPV_mostNpart_x;
0284             vtxdata.TruthPV_mostNpart_y = TruthPV_mostNpart_y;
0285             vtxdata.TruthPV_mostNpart_z = TruthPV_mostNpart_z;
0286             vtxdata.Centrality_bimp = centrality_bimp;
0287             vtxdata.Centrality_impactparam = centrality_impactparam;
0288         }
0289         vtxdata.NhitsLayer1 = (int)MVTXlayer1_halves1.size() + (int)MVTXlayer1_halves2.size();
0290         vtxdata.PV_x_halves1 = PV_halves1[0];
0291         vtxdata.PV_y_halves1 = PV_halves1[1];
0292         vtxdata.PV_z_halves1 = PV_halves1[2];
0293         vtxdata.PV_x_halves2 = PV_halves2[0];
0294         vtxdata.PV_y_halves2 = PV_halves2[1];
0295         vtxdata.PV_z_halves2 = PV_halves2[2];
0296         vtxdata.Centrality_mbd = centrality_mbd;
0297         vtxdata.Centrality_mbdquantity = centrality_mbdquantity;
0298
0299         minitree->Fill();
0300     }
0301
0302     outfile->cd();
0303     minitree->Write();
0304     outfile->Close();
0305
0306     f->Close();
0307
0308     return 0;
0309 }