analysis_INTT/src/Vertex.h

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
0012 #include <TCanvas.h>
0013 #include <TFile.h>
0014 #include <TGaxis.h>
0015 #include <TH1F.h>
0016 #include <TH2F.h>
0017 #include <TLatex.h>
0018 #include <TLegend.h>
0019 #include <TLine.h>
0020 #include <TMarker.h>
0021 #include <TRandom3.h>
0022 #include <TText.h>
0023 #include <TTree.h>
0024 #include <TVector3.h>
0025
0026 #include "/sphenix/user/hjheng/TrackletAna/analysis/plot/sPHENIXStyle/sPhenixStyle.C"
0027
0028 #define NZGAP 0
0029
0030 using namespace std;
0031
0032 struct VtxData
0033 {
0034     bool isdata, is_min_bias, isGoodVtx;
0035     vector<int> firedTriggers;
0036     int event, NClusLayer1, NTruthVtx;
0037     uint64_t INTT_BCO;
0038     float TruthPV_x, TruthPV_y, TruthPV_z;
0039     float PV_x, PV_y, PV_z;
0040     float sigmaGaus_PVz, sigmaGaus_err_PVz; // from the Gaussian fit
0041     float maxGroup_ratio, maxGroup_width;   // from the max group method
0042     float PV_z_rand; // For acceptance correction
0043     float Centrality_bimp, Centrality_impactparam, Centrality_mbd;
0044     float mbd_south_charge_sum, mbd_north_charge_sum, mbd_charge_sum, mbd_charge_asymm, mbd_z_vtx;
0045 };
0046
0047 void SetVtxMinitree(TTree *outTree, VtxData &vtxdata)
0048 {
0049     outTree->Branch("event", &vtxdata.event);
0050     outTree->Branch("INTT_BCO", &vtxdata.INTT_BCO);
0051     outTree->Branch("is_min_bias", &vtxdata.is_min_bias);
0052     outTree->Branch("isGoodVtx", &vtxdata.isGoodVtx);
0053     outTree->Branch("firedTriggers", &vtxdata.firedTriggers);
0054     outTree->Branch("NClusLayer1", &vtxdata.NClusLayer1);
0055     if (!vtxdata.isdata)
0056     {
0057         outTree->Branch("NTruthVtx", &vtxdata.NTruthVtx);
0058         outTree->Branch("TruthPV_x", &vtxdata.TruthPV_x);
0059         outTree->Branch("TruthPV_y", &vtxdata.TruthPV_y);
0060         outTree->Branch("TruthPV_z", &vtxdata.TruthPV_z);
0061         outTree->Branch("Centrality_bimp", &vtxdata.Centrality_bimp);
0062         outTree->Branch("Centrality_impactparam", &vtxdata.Centrality_impactparam);
0063     }
0064     outTree->Branch("MBD_centrality", &vtxdata.Centrality_mbd);
0065     outTree->Branch("mbd_south_charge_sum", &vtxdata.mbd_south_charge_sum);
0066     outTree->Branch("mbd_north_charge_sum", &vtxdata.mbd_north_charge_sum);
0067     outTree->Branch("mbd_charge_sum", &vtxdata.mbd_charge_sum);
0068     outTree->Branch("mbd_charge_asymm", &vtxdata.mbd_charge_asymm);
0069     outTree->Branch("mbd_z_vtx", &vtxdata.mbd_z_vtx);
0070     outTree->Branch("PV_x", &vtxdata.PV_x);
0071     outTree->Branch("PV_y", &vtxdata.PV_y);
0072     outTree->Branch("PV_z", &vtxdata.PV_z);
0073     outTree->Branch("sigmaGaus_PVz", &vtxdata.sigmaGaus_PVz);
0074     outTree->Branch("sigmaGaus_err_PVz", &vtxdata.sigmaGaus_err_PVz);
0075     outTree->Branch("maxGroup_ratio", &vtxdata.maxGroup_ratio);
0076     outTree->Branch("maxGroup_width", &vtxdata.maxGroup_width);
0077 }
0078
0079 std::map<int, vector<float>> EvtVtx_map_event(const char *vtxfname)
0080 {
0081     std::map<int, vector<float>> EvtVtx_map;
0082
0083     TFile *f = new TFile(vtxfname, "READ");
0084     TTree *t = (TTree *)f->Get("minitree");
0085     int event;
0086     float PV_x, PV_y, PV_z;
0087     t->SetBranchAddress("event", &event);
0088     t->SetBranchAddress("PV_x", &PV_x);
0089     t->SetBranchAddress("PV_y", &PV_y);
0090     t->SetBranchAddress("PV_z", &PV_z);
0091     for (int ev = 0; ev < t->GetEntriesFast(); ev++)
0092     {
0093         t->GetEntry(ev);
0094         EvtVtx_map[event] = {PV_x, PV_y, PV_z};
0095     }
0096
0097     return EvtVtx_map;
0098 }
0099
0100 std::map<uint64_t, vector<float>> EvtVtx_map_inttbco(const char *vtxfname)
0101 {
0102     std::map<uint64_t, vector<float>> EvtVtx_map;
0103
0104     TFile *f = new TFile(vtxfname, "READ");
0105     TTree *t = (TTree *)f->Get("minitree");
0106     uint64_t intt_bco;
0107     float PV_x, PV_y, PV_z;
0108     t->SetBranchAddress("INTT_BCO", &intt_bco);
0109     t->SetBranchAddress("PV_x", &PV_x);
0110     t->SetBranchAddress("PV_y", &PV_y);
0111     t->SetBranchAddress("PV_z", &PV_z);
0112     for (int ev = 0; ev < t->GetEntriesFast(); ev++)
0113     {
0114         t->GetEntry(ev);
0115         EvtVtx_map[intt_bco] = {PV_x, PV_y, PV_z};
0116     }
0117
0118     return EvtVtx_map;
0119 }
0120
0121 // Get the beamspot from file/minitree
0122 std::vector<float> getBeamspot(const char *fpath)
0123 {
0124     TString mergename = Form("%s/minitree_merged.root", fpath);
0125     // TString cmd = Form("hadd -f -j 20 %s %s/minitree_00*.root", mergename.Data(), fpath);
0126     // system(cmd.Data());
0127
0128     double beamspotx, beamspoty;
0129
0130     TFile *f = new TFile(mergename, "READ");
0131     TTree *t = (TTree *)f->Get("reco_beamspot");
0132     double bs_x, bs_y;
0133     t->SetBranchAddress("x", &bs_x);
0134     t->SetBranchAddress("y", &bs_y);
0135     int nentries = t->GetEntries();
0136     for (int i = 0; i < nentries; i++)
0137     {
0138         t->GetEntry(i);
0139         beamspotx += bs_x;
0140         beamspoty += bs_y;
0141     }
0142     f->Close();
0143
0144     beamspotx /= nentries;
0145     beamspoty /= nentries;
0146
0147     return {static_cast<float>(beamspotx), static_cast<float>(beamspoty)};
0148 }
0149
0150 TH1F *VtxZ_ReweiHist(const char *filename = "/sphenix/user/hjheng/sPHENIXRepo/analysis/dNdEta_Run2023/analysis_INTT/plot/RecoPV_ana/VtxZ_reweight_HIJING_ana419_20240910.root", const char *histname = "VtxZ_reweight_VtxZm10to10")
0151 {
0152     TFile *f = new TFile(filename, "READ");
0153     TH1F *h = (TH1F *)f->Get(histname);
0154     h->SetDirectory(0);
0155     f->Close();
0156
0157     return h;
0158 }
0159
0160 #endif