Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2025-08-05 08:11:21 

0001 #include <algorithm>
0002 #include <cmath>
0003 #include <fstream>
0004 #include <iostream>
0005 #include <numeric>
0006 #include <sstream>
0007 #include <string>
0008 #include <tuple>
0009 #include <vector>
0010 
0011 #include <TCanvas.h>
0012 #include <TF1.h>
0013 #include <TFile.h>
0014 #include <TH1F.h>
0015 #include <TH2F.h>
0016 #include <TKDE.h>
0017 #include <TRandom.h>
0018 #ifndef __CINT__
0019 #include "Math/DistFunc.h"
0020 #endif
0021 #include <TDirectory.h>
0022 #include <TLegend.h>
0023 #include <TLine.h>
0024 #include <TMath.h>
0025 #include <TObjString.h>
0026 #include <TProfile.h>
0027 #include <TRandom3.h>
0028 #include <TText.h>
0029 #include <TTree.h>
0030 #include <TTreeIndex.h>
0031 
0032 #include "Hit.h"
0033 #include "Utilities.h"
0034 #include "Vertex.h"
0035 
0036 int NevtToPlot = 10;
0037 float degreetoradian = TMath::Pi() / 180.;
0038 float radianstodegree = 180. / TMath::Pi();
0039 
0040 double gaus_func(double *x, double *par)
0041 {
0042     // note : par[0] : normalization
0043     // note : par[1] : mean
0044     // note : par[2] : width
0045     // note : par[3] : constant offset
0046     return par[0] * TMath::Gaus(x[0], par[1], par[2]) + par[3];
0047 }
0048 
0049 void draw_demoplot(TH1F *h, TF1 *f, float dcacut, TString plotname)
0050 {
0051     TCanvas *c = new TCanvas("c", "c", 800, 700);
0052     c->cd();
0053     h->GetXaxis()->SetTitle("v_{z}^{candidate} segments [cm]");
0054     h->GetYaxis()->SetTitle("Counts");
0055     h->GetYaxis()->SetTitleOffset(1.5);
0056     h->SetMaximum(h->GetMaximum() * 1.2);
0057     h->SetLineColor(1);
0058     h->SetLineWidth(2);
0059     h->Draw("hist");
0060     f->SetLineColor(kRed);
0061     f->SetLineWidth(2);
0062     f->SetNpx(1000);
0063     f->Draw("L same");
0064     TLegend *leg = new TLegend(0.47, 0.8, 0.9, 0.92);
0065     leg->SetHeader();
0066     leg->SetBorderSize(0);
0067     leg->SetFillStyle(0);
0068     leg->SetTextSize(0.035);
0069     leg->AddEntry((TObject *)0, Form("DCA < %.2f cm", dcacut), "");
0070     leg->AddEntry((TObject *)0, Form("#mu_{Gaussian} = %.2f #pm %.2f cm", f->GetParameter(1), f->GetParError(1)), "");
0071     leg->Draw();
0072     c->SaveAs(Form("%s.pdf", plotname.Data()));
0073     c->SaveAs(Form("%s.png", plotname.Data()));
0074 
0075     delete c;
0076 }
0077 
0078 std::tuple<int, double, double> find_maxGroup(TH1F *hist_in)
0079 {
0080     double maxbin_content = hist_in->GetBinContent(hist_in->GetMaximumBin());
0081 
0082     vector<double> hist_bincontent;
0083 
0084     for (int i = 1; i <= hist_in->GetNbinsX(); i++)
0085     {
0086         double c = (hist_in->GetBinContent(i) < maxbin_content / 2.0) ? 0.0 : (hist_in->GetBinContent(i) - maxbin_content / 2.0);
0087         hist_bincontent.push_back(c);
0088     }
0089 
0090     vector<vector<double>> hist_bincontent_groups;
0091     vector<double> group;
0092     for (size_t i = 0; i < hist_bincontent.size(); i++)
0093     {
0094         if (hist_bincontent[i] > 0.0)
0095         {
0096             group.push_back(hist_bincontent[i]);
0097         }
0098         else if (group.size() > 0)
0099         {
0100             hist_bincontent_groups.push_back(group);
0101             group.clear();
0102         }
0103     }
0104     if (group.size() > 0)
0105     {
0106         hist_bincontent_groups.push_back(group);
0107     }
0108 
0109     int peak_group_id = -1;
0110     int nbins = 0;
0111     double peak_group_ratio = std::numeric_limits<double>::quiet_NaN();
0112     double group_width = std::numeric_limits<double>::quiet_NaN();
0113     std::tuple<int, double, double> result;
0114 
0115     for (size_t i = 0; i < hist_bincontent_groups.size(); i++)
0116     {
0117         if (maxbin_content > hist_bincontent_groups[i].front() && maxbin_content < hist_bincontent_groups[i].back())
0118         {
0119             peak_group_id = i;
0120             nbins = hist_bincontent_groups[i].size();
0121             peak_group_ratio = accumulate(hist_bincontent_groups[i].begin(), hist_bincontent_groups[i].end(), 0.0) / accumulate(hist_bincontent.begin(), hist_bincontent.end(), 0.0);
0122             group_width = nbins * hist_in->GetBinWidth(1);
0123         }
0124     }
0125 
0126     return make_tuple(nbins, group_width, peak_group_ratio);
0127 }
0128 
0129 int main(int argc, char *argv[])
0130 {
0131     SetsPhenixStyle();
0132     gStyle->SetPalette(kThermometer);
0133 
0134     bool IsData = (TString(argv[1]).Atoi() == 1) ? true : false;        //
0135     int NevtToRun = TString(argv[2]).Atoi();                            //
0136     float dPhi_cut = TString(argv[3]).Atof();                           // Example: 0.11 radian; 0.001919862 radian = 0.11 degree
0137     float dca_cut = TString(argv[4]).Atof();                            // Example: 0.05cm
0138     TString bsresfilepath = TString(argv[5]);                           //
0139     TString infilename = TString(argv[6]);                              // /sphenix/user/hjheng/TrackletAna/data/INTT/ana382_zvtx-20cm_dummyAlignParams/sim/INTTRecoClusters_sim_merged.root
0140     TString outfilename = TString(argv[7]);                             // /sphenix/user/hjheng/TrackletAna/minitree/INTT/VtxEvtMap_ana382_zvtx-20cm_dummyAlignParams
0141     TString demoplotpath = TString(argv[8]);                            // ./plot/RecoPV_demo/RecoPV_sim/INTTVtxZ_ana382_zvtx-20cm_dummyAlignParams
0142     bool debug = (TString(argv[9]).Atoi() == 1) ? true : false;         //
0143     bool makedemoplot = (TString(argv[10]).Atoi() == 1) ? true : false; //
0144 
0145     TString idxstr = (IsData) ? "INTT_BCO" : "event";
0146 
0147     // loop argv and cout
0148     for (int i = 0; i < argc; i++)
0149     {
0150         cout << "argv[" << i << "] = " << argv[i] << endl;
0151     }
0152 
0153     // get the beamspot
0154     std::vector<float> beamspot = getBeamspot(bsresfilepath.Data());
0155     float beamspotx = beamspot[0];
0156     float beamspoty = beamspot[1];
0157     cout << "Beamspot: X = " << beamspotx << ", Y = " << beamspoty << endl;
0158 
0159     // system(Form("mkdir -p %s", outfilepath.Data()));
0160     if (makedemoplot)
0161         system(Form("mkdir -p %s", demoplotpath.Data()));
0162 
0163     vector<Hit *> INTTlayer1, INTTlayer2;
0164     VtxData vtxdata = {};
0165 
0166     TFile *f = new TFile(infilename, "READ");
0167     TTree *t = (TTree *)f->Get("EventTree");
0168     t->BuildIndex(idxstr.Data()); // Reference: https://root-forum.cern.ch/t/sort-ttree-entries/13138
0169     TTreeIndex *index = (TTreeIndex *)t->GetTreeIndex();
0170     int event, NTruthVtx;
0171     uint64_t INTT_BCO;
0172     bool is_min_bias, InttBco_IsToBeRemoved;
0173     float TruthPV_trig_x, TruthPV_trig_y, TruthPV_trig_z, centrality_bimp, centrality_impactparam, centrality_mbd;
0174     float mbd_south_charge_sum, mbd_north_charge_sum, mbd_charge_sum, mbd_charge_asymm, mbd_z_vtx;
0175     vector<int> *firedTriggers = 0;
0176     vector<int> *ClusLayer = 0;
0177     vector<float> *ClusX = 0, *ClusY = 0, *ClusZ = 0, *ClusPhiSize = 0, *ClusZSize = 0;
0178     vector<unsigned int> *ClusAdc = 0;
0179     vector<uint8_t> *ClusLadderZId = 0, *ClusLadderPhiId = 0;
0180     t->SetBranchAddress("event", &event);
0181     t->SetBranchAddress("INTT_BCO", &INTT_BCO);
0182     if (!IsData)
0183     {
0184         t->SetBranchAddress("NTruthVtx", &NTruthVtx);
0185         t->SetBranchAddress("TruthPV_trig_x", &TruthPV_trig_x);
0186         t->SetBranchAddress("TruthPV_trig_y", &TruthPV_trig_y);
0187         t->SetBranchAddress("TruthPV_trig_z", &TruthPV_trig_z);
0188         t->SetBranchAddress("centrality_bimp", &centrality_bimp);
0189         t->SetBranchAddress("centrality_impactparam", &centrality_impactparam);
0190         InttBco_IsToBeRemoved = false;
0191     }
0192     t->SetBranchAddress("is_min_bias", &is_min_bias);
0193     t->SetBranchAddress("MBD_centrality", &centrality_mbd);
0194     t->SetBranchAddress("MBD_z_vtx", &mbd_z_vtx);
0195     t->SetBranchAddress("MBD_south_charge_sum", &mbd_south_charge_sum);
0196     t->SetBranchAddress("MBD_north_charge_sum", &mbd_north_charge_sum);
0197     t->SetBranchAddress("MBD_charge_sum", &mbd_charge_sum);
0198     t->SetBranchAddress("MBD_charge_asymm", &mbd_charge_asymm);
0199     if (IsData)
0200     {
0201         t->SetBranchAddress("firedTriggers", &firedTriggers);
0202         t->SetBranchAddress("InttBco_IsToBeRemoved", &InttBco_IsToBeRemoved);
0203     }
0204     t->SetBranchAddress("ClusLayer", &ClusLayer);
0205     t->SetBranchAddress("ClusX", &ClusX);
0206     t->SetBranchAddress("ClusY", &ClusY);
0207     t->SetBranchAddress("ClusZ", &ClusZ);
0208     t->SetBranchAddress("ClusPhiSize", &ClusPhiSize);
0209     t->SetBranchAddress("ClusZSize", &ClusZSize);
0210     t->SetBranchAddress("ClusAdc", &ClusAdc);
0211     t->SetBranchAddress("ClusLadderZId", &ClusLadderZId);
0212     t->SetBranchAddress("ClusLadderPhiId", &ClusLadderPhiId);
0213 
0214     TFile *outfile = new TFile(outfilename.Data(), "RECREATE");
0215     TTree *minitree = new TTree("minitree", "Minitree of reconstructed vertices");
0216     TDirectory *dir = outfile->mkdir("hist_vtxzprojseg");
0217     SetVtxMinitree(minitree, vtxdata);
0218 
0219     for (int ev = 0; ev < NevtToRun; ev++)
0220     {
0221         Long64_t local = t->LoadTree(index->GetIndex()[ev]);
0222         t->GetEntry(local);
0223 
0224         cout << "event=" << event << ", INTT BCO=" << INTT_BCO << ", total number of clusters=" << ClusLayer->size() << endl;
0225 
0226         CleanVec(INTTlayer1);
0227         CleanVec(INTTlayer2);
0228 
0229         int NClusLayer1_beforecut = 0;
0230 
0231         for (size_t ihit = 0; ihit < ClusLayer->size(); ihit++)
0232         {
0233             if ((ClusLayer->at(ihit) < 3 || ClusLayer->at(ihit) > 6) || (ClusLadderZId->at(ihit) < 0 || ClusLadderZId->at(ihit) > 3))
0234             {
0235                 cout << "Cluster (layer, ladderZId) = (" << ClusLayer->at(ihit) << "," << ClusLadderZId->at(ihit) << ") is not in the INTT acceptance. Exit and check." << endl;
0236                 exit(1);
0237             }
0238 
0239             if (ClusLayer->at(ihit) == 3 || ClusLayer->at(ihit) == 4)
0240                 NClusLayer1_beforecut++;
0241 
0242             if (ClusAdc->at(ihit) <= 35)
0243                 continue;
0244 
0245             int layer = (ClusLayer->at(ihit) == 3 || ClusLayer->at(ihit) == 4) ? 0 : 1;
0246 
0247             Hit *hit = new Hit(ClusX->at(ihit), ClusY->at(ihit), ClusZ->at(ihit), beamspotx, beamspoty, 0., layer);
0248             float edge = (ClusLadderZId->at(ihit) == 0 || ClusLadderZId->at(ihit) == 2) ? 0.8 : 0.5;
0249             hit->SetEdge(ClusZ->at(ihit) - edge, ClusZ->at(ihit) + edge);
0250 
0251             if (layer == 0)
0252                 INTTlayer1.push_back(hit);
0253             else
0254                 INTTlayer2.push_back(hit);
0255         }
0256 
0257         cout << "# of clusters in inner layer (layer ID 3+4, after cluster phi size and adc selection) = " << INTTlayer1.size() << ", outer layer (layer ID 5+6) = " << INTTlayer2.size() << endl;
0258 
0259         TH1F *hM_vtxzprojseg = new TH1F(Form("hM_vtxzprojseg_ev%d", ev), Form("hM_vtxzprojseg_ev%d", ev), 2800, -70, 70);
0260         int goodpaircount = 0;
0261         for (size_t i = 0; i < INTTlayer1.size(); i++)
0262         {
0263             for (size_t j = 0; j < INTTlayer2.size(); j++)
0264             {
0265                 if (fabs(deltaPhi(INTTlayer1[i]->Phi(), INTTlayer2[j]->Phi())) > dPhi_cut)
0266                     continue;
0267 
0268                 // find the equation connecting two hits in the X-Y plane
0269                 // y = mx + b
0270                 float m = (INTTlayer2[j]->posY() - INTTlayer1[i]->posY()) / (INTTlayer2[j]->posX() - INTTlayer1[i]->posX());
0271                 float b = INTTlayer1[i]->posY() - m * INTTlayer1[i]->posX();
0272                 // calculate the distance of closest approach from the line to (beamspotx, beamspoty)
0273                 float dca = fabs(m * beamspotx - beamspoty + b) / sqrt(m * m + 1);
0274 
0275                 if (dca > dca_cut)
0276                     continue;
0277 
0278                 float rho1 = sqrt(pow((INTTlayer1[i]->posX() - beamspotx), 2) + pow((INTTlayer1[i]->posY() - beamspoty), 2));
0279                 float rho2 = sqrt(pow((INTTlayer2[j]->posX() - beamspotx), 2) + pow((INTTlayer2[j]->posY() - beamspoty), 2));
0280                 float z = INTTlayer1[i]->posZ() - (INTTlayer2[j]->posZ() - INTTlayer1[i]->posZ()) / (rho2 - rho1) * rho1;
0281                 float edge1 = INTTlayer1[i]->Edge().first - (INTTlayer2[j]->Edge().second - INTTlayer1[i]->Edge().first) / (rho2 - rho1) * rho1;
0282                 float edge2 = INTTlayer1[i]->Edge().second - (INTTlayer2[j]->Edge().first - INTTlayer1[i]->Edge().second) / (rho2 - rho1) * rho1;
0283 
0284                 if (debug)
0285                 {
0286                     cout << "----------" << endl << "Cluster pair (i,j) = (" << i << "," << j << "); position (x1,y1,z1,x2,y2,z2) mm =(" << INTTlayer1[i]->posX() * 10. << "," << INTTlayer1[i]->posY() * 10. << "," << INTTlayer1[i]->posZ() * 10. << "," << INTTlayer2[j]->posX() * 10. << "," << INTTlayer2[j]->posY() * 10. << "," << INTTlayer2[j]->posZ() * 10. << ")" << endl;
0287                     cout << "Cluster [phi1 (deg),eta1,phi2(deg),eta2]=[" << INTTlayer1[i]->Phi() * radianstodegree << "," << INTTlayer1[i]->Eta() << "," << INTTlayer2[j]->Phi() * radianstodegree << "," << INTTlayer2[j]->Eta() << "]" << endl
0288                          << "delta phi (in degree) = " << fabs(deltaPhi(INTTlayer1[i]->Phi(), INTTlayer2[j]->Phi())) * radianstodegree << "-> pass dPhi selection (<" << dPhi_cut << " rad)=" << ((fabs(deltaPhi(INTTlayer1[i]->Phi(), INTTlayer2[j]->Phi())) < dPhi_cut) ? 1 : 0) << endl;
0289                     cout << "DCA cut = " << dca_cut << " [cm]; vertex candidate (center,edge1,edge2) = (" << z << "," << edge1 << "," << edge2 << "), difference = " << edge2 - edge1 << endl;
0290                     cout << "dca w.r.t (beamspotx [cm], beamspoty [cm]) (in mm) = " << dca * 10 << endl;
0291                 }
0292 
0293                 goodpaircount++;
0294                 if (fabs(z) < 70)
0295                 {
0296                     int bin1 = hM_vtxzprojseg->FindBin(edge1);
0297                     int bin2 = hM_vtxzprojseg->FindBin(edge2);
0298 
0299                     for (int ibin = bin1; ibin <= bin2; ibin++)
0300                     {
0301                         float bincenter = hM_vtxzprojseg->GetBinCenter(ibin);
0302                         float binwidth = hM_vtxzprojseg->GetBinWidth(ibin);
0303                         // for bin1 and bin2, find the distance of the edge1 and edge2 to the respective bin edge
0304                         float w;
0305                         if (ibin == bin1)
0306                             w = ((bincenter + 0.5 * binwidth) - edge1) / binwidth;
0307                         else if (ibin == bin2)
0308                             w = (edge2 - (bincenter - 0.5 * binwidth)) / binwidth;
0309                         else
0310                             w = 1.;
0311 
0312                         hM_vtxzprojseg->Fill(bincenter, w);
0313                     }
0314                 }
0315             }
0316         }
0317 
0318         cout << "Number of entries in hM_vtxzprojseg = " << hM_vtxzprojseg->Integral(-1, -1) << endl;
0319         cout << "Number of good pairs = " << goodpaircount << endl;
0320         if (debug && !IsData)
0321             cout << "Event " << ev << " (Truth PVx, Truth PVy, Truth PVz) = (" << TruthPV_trig_x << ", " << TruthPV_trig_y << ", " << TruthPV_trig_z << ")" << endl;
0322 
0323         // find the maximum bin of the histogram hM_vtxzedge_dca[2] and fit the histogram with a Gaussian function around the maximum bin
0324         float maxbincenter = hM_vtxzprojseg->GetBinCenter(hM_vtxzprojseg->GetMaximumBin());
0325 
0326         TF1 *f1 = new TF1("gaus_fit", gaus_func, maxbincenter - 10, maxbincenter + 10, 4); // Gaussian + const. offset
0327         f1->SetParameters(hM_vtxzprojseg->GetBinContent(hM_vtxzprojseg->GetMaximumBin()), hM_vtxzprojseg->GetBinCenter(hM_vtxzprojseg->GetMaximumBin()), 4, 0);
0328         f1->SetParLimits(0, 0, 10000);
0329         f1->SetParLimits(2, 0, 1000);
0330         f1->SetParLimits(3, 0, 1000);
0331         hM_vtxzprojseg->Fit(f1, "NQ", "", hM_vtxzprojseg->GetBinCenter(hM_vtxzprojseg->GetMaximumBin()) - 9, hM_vtxzprojseg->GetBinCenter(hM_vtxzprojseg->GetMaximumBin()) + 9);
0332         f1->SetParameters(hM_vtxzprojseg->GetBinContent(hM_vtxzprojseg->GetMaximumBin()), hM_vtxzprojseg->GetBinCenter(hM_vtxzprojseg->GetMaximumBin()), 4, 0);
0333         f1->SetParLimits(0, 0, 10000);
0334         f1->SetParLimits(2, 0, 1000);
0335         f1->SetParLimits(3, -20, 1000);
0336         hM_vtxzprojseg->Fit(f1, "NQ", "", hM_vtxzprojseg->GetBinCenter(hM_vtxzprojseg->GetMaximumBin()) - 9, hM_vtxzprojseg->GetBinCenter(hM_vtxzprojseg->GetMaximumBin()) + 9);
0337         float mean = f1->GetParameter(1);
0338         float meanErr = f1->GetParError(1);
0339         float sigma = f1->GetParameter(2);
0340         float sigmaErr = f1->GetParError(2);
0341 
0342         std::tuple<int, double, double> maxGroup = find_maxGroup(hM_vtxzprojseg);
0343         int maxGroup_nbins = std::get<0>(maxGroup);
0344         double maxGroup_width = std::get<1>(maxGroup);
0345         double maxGroup_ratio = std::get<2>(maxGroup);
0346 
0347         if (debug)
0348             cout << "Event " << ev << " Reco PVz = " << mean << " +/- " << meanErr << " cm" << endl;
0349 
0350         if (ev < NevtToPlot && makedemoplot)
0351         {
0352             TString pname = Form("%s/hM_vtxzprojseg_ev%d", demoplotpath.Data(), event);
0353             draw_demoplot(hM_vtxzprojseg, f1, dca_cut, pname);
0354         }
0355 
0356         // Vertex quality
0357         bool isValidMBDZvtx = (mbd_z_vtx != mbd_z_vtx) ? false : true;
0358         bool isGoodINTTGausWidth = (sigma >= 3.0 && sigma <= 8.0) ? true : false;
0359         bool isMBDINTTZvtxMatch = (IsData && (mean - mbd_z_vtx > -5. && mean - mbd_z_vtx < 3.)) || (!IsData && fabs(mean - mbd_z_vtx) < 4.);
0360         bool isGoodINTTGoodPeak = (maxGroup_width > 4. && maxGroup_width < 11.) ? true : false;
0361 
0362         vtxdata.isdata = IsData;
0363         vtxdata.is_min_bias = is_min_bias;
0364         vtxdata.isGoodVtx = isValidMBDZvtx && isGoodINTTGausWidth && isGoodINTTGoodPeak && isMBDINTTZvtxMatch;
0365         vtxdata.event = event;
0366         vtxdata.INTT_BCO = INTT_BCO;
0367         vtxdata.firedTriggers = (IsData) ? *firedTriggers : vector<int>();
0368         if (!IsData)
0369         {
0370             vtxdata.NTruthVtx = NTruthVtx;
0371             vtxdata.TruthPV_x = TruthPV_trig_x;
0372             vtxdata.TruthPV_y = TruthPV_trig_y;
0373             vtxdata.TruthPV_z = TruthPV_trig_z;
0374             vtxdata.Centrality_bimp = centrality_bimp;
0375             vtxdata.Centrality_impactparam = centrality_impactparam;
0376         }
0377         vtxdata.NClusLayer1 = NClusLayer1_beforecut;
0378         vtxdata.PV_x = beamspotx;
0379         vtxdata.PV_y = beamspoty;
0380         vtxdata.PV_z = mean;
0381         vtxdata.sigmaGaus_PVz = sigma;
0382         vtxdata.sigmaGaus_err_PVz = sigmaErr;
0383         vtxdata.maxGroup_ratio = maxGroup_ratio;
0384         vtxdata.maxGroup_width = maxGroup_width;
0385         vtxdata.Centrality_mbd = centrality_mbd;
0386         vtxdata.mbd_south_charge_sum = mbd_south_charge_sum;
0387         vtxdata.mbd_north_charge_sum = mbd_north_charge_sum;
0388         vtxdata.mbd_charge_sum = mbd_charge_sum;
0389         vtxdata.mbd_charge_asymm = mbd_charge_asymm;
0390         vtxdata.mbd_z_vtx = mbd_z_vtx;
0391 
0392         minitree->Fill();
0393         if (makedemoplot)
0394         {
0395             dir->cd();
0396             hM_vtxzprojseg->Write();
0397         }
0398     }
0399 
0400     outfile->cd();
0401     minitree->Write();
0402     outfile->Close();
0403 
0404     f->Close();
0405 
0406     return 0;
0407 }
[ 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!