File indexing completed on 2025-08-05 08:13:29
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();
0039 float gap_upper = TString(argv[5]).Atof();
0040 float cent_shift = TString(argv[6]).Atof();
0041 TString gap_north_str = TString(argv[4]).Replace(TString(argv[4]).Index("."), 1, "p");
0042 TString gap_upper_str = TString(argv[5]).Replace(TString(argv[5]).Index("."), 1, "p");
0043 TString cent_shift_str = TString(argv[6]).Replace(TString(argv[6]).Index("."), 1, "p");
0044
0045
0046
0047 SetsPhenixStyle();
0048
0049 TString infilename = "";
0050 TString outfilename = "";
0051 TString demoplotname = "";
0052 if (IsData)
0053 {
0054 cout << "Running on data" << endl;
0055
0056
0057
0058
0059
0060 infilename = "/sphenix/user/hjheng/macros/detectors/sPHENIX/OutFile/ClusterTree/ClusterTree_Run20124_0.root";
0061
0062 outfilename = Form("/sphenix/user/hjheng/TrackletAna/minitree/Data_Run2023_AuAu/TrackletAna_RecoClusters_RecoVtx_TklCluster_dPhiCutbin%d_dZCutbin%d_GapNorth%s_GapUpper%s_CentShift%s_3DVertex_twohalves.root", dPhiCutbin, dZCutbin, gap_north_str.Data(), gap_upper_str.Data(), cent_shift_str.Data());
0063
0064 demoplotname = "./plot/RecoPV_demo/RecoPV_Data/Run20124";
0065
0066 }
0067 else
0068 {
0069 cout << "Running on MC" << endl;
0070
0071
0072
0073
0074
0075
0076 infilename = "/sphenix/user/hjheng/TrackletAna/data/MVTXRecoClusters/NoPileup_Nevt500_ana325private_singleEvtDst/MVTXRecoClusters.root";
0077 outfilename = Form("/sphenix/user/hjheng/TrackletAna/minitree/AuAu_ana325private_NoPileup_RecoVtx_Optimization/TrackletAna_RecoClusters_RecoVtx_TklCluster_dPhiCutbin%d_dZCutbin%d_GapNorth%s_GapUpper%s_CentShift%s_3DVertex_twohalves.root", dPhiCutbin, dZCutbin, gap_north_str.Data(), gap_upper_str.Data(), cent_shift_str.Data());
0078 demoplotname = "./plot/RecoPV_demo/RecoPV_Sim/AuAuHijing_FieldOn";
0079
0080
0081
0082
0083 }
0084
0085 vector<Hit *> MVTXlayer1_halves1, MVTXlayer1_halves2, MVTXlayer2_halves1, MVTXlayer2_halves2;
0086
0087 VtxData vtxdata = {};
0088
0089 TFile *f = new TFile(infilename, "READ");
0090 TTree *t = (TTree *)f->Get("EventTree");
0091 t->BuildIndex("event");
0092 TTreeIndex *index = (TTreeIndex *)t->GetTreeIndex();
0093 int event, NTruthVtx;
0094 float TruthPV_trig_x, TruthPV_trig_y, TruthPV_trig_z, TruthPV_mostNpart_x, TruthPV_mostNpart_y, TruthPV_mostNpart_z;
0095 float centrality_bimp, centrality_impactparam, centrality_mbd, centrality_mbdquantity;
0096 vector<int> *ClusLayer = 0;
0097 vector<float> *ClusX = 0, *ClusY = 0, *ClusZ = 0;
0098 t->SetBranchAddress("event", &event);
0099 if (!IsData)
0100 {
0101 t->SetBranchAddress("NTruthVtx", &NTruthVtx);
0102 t->SetBranchAddress("TruthPV_trig_x", &TruthPV_trig_x);
0103 t->SetBranchAddress("TruthPV_trig_y", &TruthPV_trig_y);
0104 t->SetBranchAddress("TruthPV_trig_z", &TruthPV_trig_z);
0105 t->SetBranchAddress("TruthPV_mostNpart_x", &TruthPV_mostNpart_x);
0106 t->SetBranchAddress("TruthPV_mostNpart_y", &TruthPV_mostNpart_y);
0107 t->SetBranchAddress("TruthPV_mostNpart_z", &TruthPV_mostNpart_z);
0108
0109
0110 }
0111
0112
0113 t->SetBranchAddress("ClusLayer", &ClusLayer);
0114 t->SetBranchAddress("ClusX", &ClusX);
0115 t->SetBranchAddress("ClusY", &ClusY);
0116 t->SetBranchAddress("ClusZ", &ClusZ);
0117
0118 TFile *outfile = new TFile(outfilename, "RECREATE");
0119 TTree *minitree = new TTree("minitree", "Minitree of reconstructed vertices");
0120 SetVtxMinitree(minitree, vtxdata);
0121
0122 for (int i = 0; i < index->GetN(); i++)
0123
0124
0125 {
0126 Long64_t local = t->LoadTree(index->GetIndex()[i]);
0127 t->GetEntry(local);
0128
0129 cout << "event=" << event << " has a total of " << ClusLayer->size() << " clusters" << endl;
0130
0131 CleanVec(MVTXlayer1_halves1);
0132 CleanVec(MVTXlayer1_halves2);
0133 CleanVec(MVTXlayer2_halves1);
0134 CleanVec(MVTXlayer2_halves2);
0135
0136 for (size_t ihit = 0; ihit < ClusLayer->size(); ihit++)
0137 {
0138
0139 float x0 = ClusX->at(ihit);
0140 float y0 = ClusY->at(ihit);
0141 float z0 = ClusZ->at(ihit);
0142 vector<float> cpos = {x0, y0, z0};
0143 UpdatePos_GapTwoHalves(cpos, gap_north, cent_shift, gap_upper);
0144
0145 if (ClusLayer->at(ihit) == 0)
0146 {
0147 Hit *hit = new Hit(cpos[0], cpos[1], cpos[2], 0., 0., 0., 0);
0148
0149 if (x0 > 0)
0150 {
0151 MVTXlayer1_halves1.push_back(hit);
0152 }
0153
0154 else
0155 MVTXlayer1_halves2.push_back(hit);
0156 }
0157 else if (ClusLayer->at(ihit) == 1)
0158 {
0159 Hit *hit = new Hit(cpos[0], cpos[1], cpos[2], 0., 0., 0., 1);
0160
0161 if (x0 > 0)
0162 MVTXlayer2_halves1.push_back(hit);
0163 else
0164 MVTXlayer2_halves2.push_back(hit);
0165 }
0166 else
0167 continue;
0168 }
0169
0170
0171 cout << "# of clusters in Layer 1 = " << MVTXlayer1_halves1.size() + MVTXlayer1_halves2.size() << ", Layer 2 = " << MVTXlayer2_halves1.size() + MVTXlayer2_halves2.size() << endl;
0172
0173
0174
0175
0176
0177
0178 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);
0179 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);
0180 vector<float> PV_halves1 = get<0>(tuple_PVest1);
0181 vector<float> PV_halves2 = get<0>(tuple_PVest2);
0182 int Nvtxcandidates_halves1 = get<1>(tuple_PVest1);
0183 int Nvtxcandidates_halves2 = get<1>(tuple_PVest2);
0184 TH2F *hM_clusterzrho_halves1 = get<2>(tuple_PVest1);
0185 TH2F *hM_clusterzrho_halves2 = get<2>(tuple_PVest2);
0186 TH2F *hM_clusterxy_halves1 = get<3>(tuple_PVest1);
0187 TH2F *hM_clusterxy_halves2 = get<3>(tuple_PVest2);
0188 vector<TLine *> tl_zrho_halves1 = get<4>(tuple_PVest1);
0189 vector<TLine *> tl_zrho_halves2 = get<4>(tuple_PVest2);
0190 vector<TLine *> tl_xy_halves1 = get<5>(tuple_PVest1);
0191 vector<TLine *> tl_xy_halves2 = get<5>(tuple_PVest2);
0192 TH1F *hM_diffExptzPVz_halves1 = get<6>(tuple_PVest1);
0193 TH1F *hM_diffExptzPVz_halves2 = get<6>(tuple_PVest2);
0194 TH2F *hM_diffExptxyPVxy_halves1 = get<7>(tuple_PVest1);
0195 TH2F *hM_diffExptxyPVxy_halves2 = get<7>(tuple_PVest2);
0196 TH1F *hM_VtxCandZ_halves1 = get<8>(tuple_PVest1);
0197 TH1F *hM_VtxCandZ_halves2 = get<8>(tuple_PVest2);
0198 TH2F *hM_VtxCandXY_halves1 = get<9>(tuple_PVest1);
0199 TH2F *hM_VtxCandXY_halves2 = get<9>(tuple_PVest2);
0200 TH1F *hM_DCAxy_VtxClus_halves1 = get<10>(tuple_PVest1);
0201 TH1F *hM_DCAxy_VtxClus_halves2 = get<10>(tuple_PVest2);
0202 TH2F *hM_DCAxy_DCAz_VtxClus_halves1 = get<11>(tuple_PVest1);
0203 TH2F *hM_DCAxy_DCAz_VtxClus_halves2 = get<11>(tuple_PVest2);
0204 TH2F *hM_DCAx_DCAy_VtxClus_halves1 = get<12>(tuple_PVest1);
0205 TH2F *hM_DCAx_DCAy_VtxClus_halves2 = get<12>(tuple_PVest2);
0206
0207
0208 float PV_x = (PV_halves1[0] * Nvtxcandidates_halves1 + PV_halves2[0] * Nvtxcandidates_halves2) / (Nvtxcandidates_halves1 + Nvtxcandidates_halves2);
0209 float PV_y = (PV_halves1[1] * Nvtxcandidates_halves1 + PV_halves2[1] * Nvtxcandidates_halves2) / (Nvtxcandidates_halves1 + Nvtxcandidates_halves2);
0210 float PV_z = (PV_halves1[2] * Nvtxcandidates_halves1 + PV_halves2[2] * Nvtxcandidates_halves2) / (Nvtxcandidates_halves1 + Nvtxcandidates_halves2);
0211
0212 if (MVTXlayer1_halves1.size() + MVTXlayer1_halves2.size() + MVTXlayer2_halves1.size() + MVTXlayer2_halves2.size() < Ncluster_toplot && i <= 60)
0213 {
0214 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));
0215 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));
0216
0217
0218 }
0219
0220
0221 hM_clusterzrho_halves1->Add(hM_clusterzrho_halves2);
0222 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));
0223 hM_diffExptxyPVxy_halves1->Add(hM_diffExptxyPVxy_halves2);
0224 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_Evt%d_diffExptxyPVxy", demoplotname.Data(), event));
0225
0226 hM_VtxCandZ_halves1->Add(hM_VtxCandZ_halves2);
0227 Draw_1DHistVtxDemo(hM_VtxCandZ_halves1, "v_{Z}^{Candidate} [cm]", "cm", 510, IsData, true, "2023-08-30", Form("%s_Evt%d_VtxCandZall", demoplotname.Data(), event));
0228 hM_VtxCandXY_halves1->Add(hM_VtxCandXY_halves2);
0229 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));
0230
0231 hM_DCAxy_VtxClus_halves1->Add(hM_DCAxy_VtxClus_halves2);
0232 Draw_1DHistVtxDemo(hM_DCAxy_VtxClus_halves1, "DCA_{xy} [cm]", "cm", 505, IsData, false, "", Form("%s_Evt%d_DCAxy", demoplotname.Data(), event));
0233 hM_DCAxy_DCAz_VtxClus_halves1->Add(hM_DCAxy_DCAz_VtxClus_halves2);
0234 Draw_2DHistVtxDemo(hM_DCAxy_DCAz_VtxClus_halves1, "DCA_{xy} [cm]", "DCA_{z} [cm]", 505, IsData, false, "", Form("%s_Evt%d_DCAxyDCAz", demoplotname.Data(), event));
0235 hM_DCAx_DCAy_VtxClus_halves1->Add(hM_DCAx_DCAy_VtxClus_halves2);
0236 Draw_2DHistVtxDemo(hM_DCAx_DCAy_VtxClus_halves1, "DCA_{x} [cm]", "DCA_{y} [cm]", 505, IsData, false, "", Form("%s_Evt%d_DCAxDCAy", demoplotname.Data(), event));
0237
0238 CleanVec(tl_zrho_halves1);
0239 CleanVec(tl_zrho_halves2);
0240 CleanVec(tl_xy_halves1);
0241 CleanVec(tl_xy_halves2);
0242 delete hM_clusterzrho_halves1;
0243 delete hM_clusterzrho_halves2;
0244 delete hM_clusterxy_halves1;
0245 delete hM_clusterxy_halves2;
0246 delete hM_diffExptzPVz_halves1;
0247 delete hM_diffExptzPVz_halves2;
0248 delete hM_diffExptxyPVxy_halves1;
0249 delete hM_diffExptxyPVxy_halves2;
0250 delete hM_VtxCandZ_halves1;
0251 delete hM_VtxCandZ_halves2;
0252 delete hM_VtxCandXY_halves1;
0253 delete hM_VtxCandXY_halves2;
0254 delete hM_DCAxy_VtxClus_halves1;
0255 delete hM_DCAxy_VtxClus_halves2;
0256 delete hM_DCAxy_DCAz_VtxClus_halves1;
0257 delete hM_DCAxy_DCAz_VtxClus_halves2;
0258 delete hM_DCAx_DCAy_VtxClus_halves1;
0259 delete hM_DCAx_DCAy_VtxClus_halves2;
0260
0261 if (!IsData)
0262 {
0263 cout << "NTruthVtx = " << NTruthVtx << endl
0264 << "Truth primary (most N particles) Z position = " << TruthPV_trig_z << endl
0265 << "Recontructed primary vertex Z position from 1st halve= " << PV_halves1[2] << endl
0266 << "Recontructed primary vertex Z position from 2nd halve= " << PV_halves2[2] << endl
0267 << "Average of the Z positions from 1st and 2nd halve= " << (PV_halves1[2] + PV_halves2[2]) / 2. << endl
0268 << "(Truth-Reco) = " << TruthPV_trig_z - (PV_halves1[2] + PV_halves2[2]) / 2. << endl;
0269
0270 if (NTruthVtx == 1)
0271 {
0272 if (TruthPV_trig_z != TruthPV_mostNpart_z)
0273 cout << "For events with only 1 truth vertex, (TruthPV_trig_z,TruthPV_mostNpart_z) = (" << TruthPV_trig_z << "," << TruthPV_mostNpart_z << "). Check!" << endl;
0274 }
0275 }
0276 else
0277 {
0278 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;
0279 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;
0280 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;
0281 }
0282
0283 vtxdata.event = event;
0284 if (!IsData)
0285 {
0286 vtxdata.NTruthVtx = NTruthVtx;
0287 vtxdata.TruthPV_trig_x = TruthPV_trig_x;
0288 vtxdata.TruthPV_trig_y = TruthPV_trig_y;
0289 vtxdata.TruthPV_trig_z = TruthPV_trig_z;
0290 vtxdata.TruthPV_mostNpart_x = TruthPV_mostNpart_x;
0291 vtxdata.TruthPV_mostNpart_y = TruthPV_mostNpart_y;
0292 vtxdata.TruthPV_mostNpart_z = TruthPV_mostNpart_z;
0293 vtxdata.Centrality_bimp = centrality_bimp;
0294 vtxdata.Centrality_impactparam = centrality_impactparam;
0295 }
0296 vtxdata.NhitsLayer1 = (int)MVTXlayer1_halves1.size() + (int)MVTXlayer1_halves2.size();
0297 vtxdata.PV_x_halves1 = PV_halves1[0];
0298 vtxdata.PV_y_halves1 = PV_halves1[1];
0299 vtxdata.PV_z_halves1 = PV_halves1[2];
0300 vtxdata.PV_x_halves2 = PV_halves2[0];
0301 vtxdata.PV_y_halves2 = PV_halves2[1];
0302 vtxdata.PV_z_halves2 = PV_halves2[2];
0303 vtxdata.Centrality_mbd = centrality_mbd;
0304 vtxdata.Centrality_mbdquantity = centrality_mbdquantity;
0305
0306 minitree->Fill();
0307 }
0308
0309 outfile->cd();
0310 minitree->Write();
0311 outfile->Close();
0312
0313 f->Close();
0314
0315 return 0;
0316 }