Back to home page

sPhenix code displayed by LXR

 
 

    


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();  // 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     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     // float dalpha_num = (gap_north-1.5)/2.;
0046 
0047     SetsPhenixStyle();
0048 
0049     TString infilename = "";
0050     TString outfilename = "";
0051     TString demoplotname = "";
0052     if (IsData)
0053     {
0054         cout << "Running on data" << endl;
0055         // Field Off - Mock Data
0056         // infilename = "/sphenix/user/hjheng/TrackletAna/data/MVTXRecoClusters/HijingMBwoPUB0_private/MockData/MVTXRecoClusters_HijingMBwoPU0T_MockData.root";
0057         // outfilename = Form("/sphenix/user/hjheng/TrackletAna/minitree/HijingAuAuMB_NoPileup_0T_RecoVtx_Optimization_MockData/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());
0058 
0059         // Real collisions - Run 20124
0060         infilename = "/sphenix/user/hjheng/macros/detectors/sPHENIX/OutFile/ClusterTree/ClusterTree_Run20124_0.root";
0061         // infilename = "/sphenix/user/hjheng/macros/detectors/sPHENIX/OutFile/ClusterTree/ClusterTree_Run20124_0_larger_subset.root";
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         // outfilename = Form("/sphenix/user/hjheng/TrackletAna/minitree/Data_Run2023_AuAu_larger_subset/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());
0064         demoplotname = "./plot/RecoPV_demo/RecoPV_Data/Run20124";
0065         // demoplotname = "./plot/RecoPV_demo/RecoPV_Data/Run20124_larger_subset";
0066     }
0067     else
0068     {
0069         cout << "Running on MC" << endl;
0070 
0071         // Field Off - Mock Sim
0072         // infilename = "/sphenix/user/hjheng/TrackletAna/data/MVTXRecoClusters/HijingMBwoPUB0_private/MockSim/MVTXRecoClusters_HijingMBwoPU0T_MockSim.root";
0073         // outfilename = Form("/sphenix/user/hjheng/TrackletAna/minitree/HijingAuAuMB_NoPileup_0T_RecoVtx_Optimization_MockSim/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());
0074 
0075         // Field On - Mock Sim
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         // Simple pion events
0081         // infilename = "/sphenix/user/hjheng/TrackletAna/data/MVTXRecoClusters/Pion/Ntuple_Pions_Npart500_VtxZmean0cm_VtxZwidth0cm_pT0p001to5GeV.root";
0082         // outfilename = Form("/sphenix/user/hjheng/TrackletAna/minitree/Pion_RecoVtx_Optimization/SimplePion_Npart500_VtxZmean0cm_VtxZwidth0cm_pT0p001to5GeV/TrackletAna_RecoClusters_SimplePion_RecoVtx_TklCluster_dPhiCutbin%d_dZCutbin%d.root", dPhiCutbin, dZCutbin);
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"); // Reference: https://root-forum.cern.ch/t/sort-ttree-entries/13138
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         // t->SetBranchAddress("centrality_bimp", &centrality_bimp);
0109         // t->SetBranchAddress("centrality_impactparam", &centrality_impactparam);
0110     }
0111     // t->SetBranchAddress("centrality_mbd", &centrality_mbd);
0112     // t->SetBranchAddress("centrality_mbdquantity", &centrality_mbdquantity);
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     // for (int i = 0; i < 1; i++)
0124     // for (Long64_t ev = 0; ev < t->GetEntriesFast(); ev++)
0125     {
0126         Long64_t local = t->LoadTree(index->GetIndex()[i]);
0127         t->GetEntry(local);
0128         // cout << "event = " << event << " local = " << local << endl;
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             // Misalignment study for MVTX
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                 // Hit *hit = new Hit(ClusX->at(ihit), ClusY->at(ihit), ClusZ->at(ihit), 0., 0., 0., 0);
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                 // Hit *hit = new Hit(ClusX->at(ihit), ClusY->at(ihit), ClusZ->at(ihit), 0., 0., 0., 1);
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         // cout << "# of clusters in Layer 1 = " << MVTXlayer1.size() << ", Layer 2 = " << MVTXlayer2.size() << endl;
0171         cout << "# of clusters in Layer 1 = " << MVTXlayer1_halves1.size() + MVTXlayer1_halves2.size() << ", Layer 2 = " << MVTXlayer2_halves1.size() + MVTXlayer2_halves2.size() << endl;
0172 
0173         // float PV_z = TrackletPV_cluster(event, MVTXlayer1, MVTXlayer2, dPhiCutbin, dZCutbin, false);
0174         // Tracklet3DPV(int evt, vector<Hit *> layer1, vector<Hit *> layer2, int dPhiCutbin, int dZCutbin, const char *plotname)
0175         // vector<float> PV_halves1 = Tracklet3DPV(event, MVTXlayer1_halves1, MVTXlayer2_halves1, dPhiCutbin, dZCutbin, Form("./plot/RecoPV_Data/Run20124_Evt%d_RecoPV_halve1", event));
0176         // vector<float> PV_halves2 = Tracklet3DPV(event, MVTXlayer1_halves2, MVTXlayer2_halves2, dPhiCutbin, dZCutbin, Form("./plot/RecoPV_Data/Run20124_Evt%d_RecoPV_halve2", event));
0177         // 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
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         // weighted average of the PV z position from the 1st and 2nd halves
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             // Draw_diffExptzPVz(hM_diffExptzPVz_halves1, hM_diffExptzPVz_halves2, IsData, Form("%s_Evt%d_diffExptzPVz", demoplotname.Data(), event));
0217             // Draw_diffExptxyPVxy(hM_diffExptxyPVxy_halves1, hM_diffExptxyPVxy_halves2, IsData, Form("%s_Evt%d_diffExptxyPVxy", demoplotname.Data(), event));
0218         }
0219 
0220         // merge the histograms
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 }