Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
INTTVtxZ.cxx
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file INTTVtxZ.cxx
1 #include <cmath>
2 #include <fstream>
3 #include <iostream>
4 #include <sstream>
5 #include <string>
6 #include <vector>
7 
8 #include <TCanvas.h>
9 #include <TF1.h>
10 #include <TFile.h>
11 #include <TH1F.h>
12 #include <TH2F.h>
13 #include <TKDE.h>
14 #include <TRandom.h>
15 #ifndef __CINT__
16 #include "Math/DistFunc.h"
17 #endif
18 #include <TLegend.h>
19 #include <TLine.h>
20 #include <TMath.h>
21 #include <TObjString.h>
22 #include <TProfile.h>
23 #include <TRandom3.h>
24 #include <TText.h>
25 #include <TTree.h>
26 #include <TTreeIndex.h>
27 
28 #include "Hit.h"
29 #include "Utilities.h"
30 #include "Vertex.h"
31 
32 int NevtToPlot = 10;
33 float degreetoradian = TMath::Pi() / 180.;
34 float radianstodegree = 180. / TMath::Pi();
35 
36 double gaus_func(double *x, double *par)
37 {
38  // note : par[0] : normalization
39  // note : par[1] : mean
40  // note : par[2] : width
41  // note : par[3] : constant offset
42  return par[0] * TMath::Gaus(x[0], par[1], par[2]) + par[3];
43 }
44 
45 void draw_demoplot(TH1F *h, TF1 *f, float dcacut, TString plotname)
46 {
47  TCanvas *c = new TCanvas("c", "c", 800, 700);
48  c->cd();
49  h->GetXaxis()->SetTitle("v_{z}^{candidate} segments [cm]");
50  h->GetYaxis()->SetTitle("Counts");
51  h->GetYaxis()->SetTitleOffset(1.5);
52  h->SetMaximum(h->GetMaximum() * 1.2);
53  h->SetLineColor(1);
54  h->SetLineWidth(2);
55  h->Draw("hist");
56  f->SetLineColor(kRed);
57  f->SetLineWidth(2);
58  f->SetNpx(1000);
59  f->Draw("L same");
60  TLegend *leg = new TLegend(0.47, 0.8, 0.9, 0.92);
61  leg->SetHeader();
62  leg->SetBorderSize(0);
63  leg->SetFillStyle(0);
64  leg->SetTextSize(0.035);
65  leg->AddEntry((TObject *)0, Form("DCA < %.2f cm", dcacut), "");
66  leg->AddEntry((TObject *)0, Form("#mu_{Gaussian} = %.2f #pm %.2f cm", f->GetParameter(1), f->GetParError(1)), "");
67  leg->Draw();
68  c->SaveAs(Form("%s.pdf", plotname.Data()));
69  c->SaveAs(Form("%s.png", plotname.Data()));
70 
71  delete c;
72 }
73 
74 int main(int argc, char *argv[])
75 {
76  SetsPhenixStyle();
77  gStyle->SetPalette(kThermometer);
78 
79  bool IsData = (TString(argv[1]).Atoi() == 1) ? true : false;
80  int NevtToRun = TString(argv[2]).Atoi();
81  float avgVtxX = TString(argv[3]).Atof(); // float avgVtxX = -0.0015 (ana.382) / -0.04 (ana.398) /
82  float avgVtxY = TString(argv[4]).Atof(); // avgVtxY = 0.0012 (ana.382) / 0.24 (ana.398) /
83  float dPhi_cut = TString(argv[5]).Atof(); // Example: 0.11 radian; 0.001919862 radian = 0.11 degree
84  float dca_cut = TString(argv[6]).Atof(); // Example: 0.05cm
85  TString infilename = TString(argv[7]); // /sphenix/user/hjheng/TrackletAna/data/INTT/ana382_zvtx-20cm_dummyAlignParams/sim/INTTRecoClusters_sim_merged.root
86  TString outfilename = TString(argv[8]); // /sphenix/user/hjheng/TrackletAna/minitree/INTT/VtxEvtMap_ana382_zvtx-20cm_dummyAlignParams
87  TString demoplotpath = TString(argv[9]); // ./plot/RecoPV_demo/RecoPV_sim/INTTVtxZ_ana382_zvtx-20cm_dummyAlignParams
88  bool debug = (TString(argv[10]).Atoi() == 1) ? true : false;
89  bool makedemoplot = (TString(argv[11]).Atoi() == 1) ? true : false;
90 
91  // int initevt = process * NevtToRun;
92 
93  // loop argv and cout
94  for (int i = 0; i < argc; i++)
95  {
96  cout << "argv[" << i << "] = " << argv[i] << endl;
97  }
98 
99  // system(Form("mkdir -p %s", outfilepath.Data()));
100  if (makedemoplot)
101  system(Form("mkdir -p %s", demoplotpath.Data()));
102 
103  vector<Hit *> INTTlayer1, INTTlayer2;
104  VtxData vtxdata = {};
105 
106  TFile *f = new TFile(infilename, "READ");
107  TTree *t = (TTree *)f->Get("EventTree");
108  t->BuildIndex("event"); // Reference: https://root-forum.cern.ch/t/sort-ttree-entries/13138
109  TTreeIndex *index = (TTreeIndex *)t->GetTreeIndex();
110  int event, NTruthVtx;
111  float TruthPV_trig_x, TruthPV_trig_y, TruthPV_trig_z, centrality_bimp, centrality_impactparam, centrality_mbd, centrality_mbdquantity;
112  vector<int> *ClusLayer = 0;
113  vector<float> *ClusX = 0, *ClusY = 0, *ClusZ = 0, *ClusPhiSize = 0, *ClusZSize = 0;
114  vector<uint8_t> *ClusLadderZId = 0, *ClusLadderPhiId = 0;
115  t->SetBranchAddress("event", &event);
116  if (!IsData)
117  {
118  t->SetBranchAddress("NTruthVtx", &NTruthVtx);
119  t->SetBranchAddress("TruthPV_trig_x", &TruthPV_trig_x);
120  t->SetBranchAddress("TruthPV_trig_y", &TruthPV_trig_y);
121  t->SetBranchAddress("TruthPV_trig_z", &TruthPV_trig_z);
122  t->SetBranchAddress("centrality_bimp", &centrality_bimp);
123  t->SetBranchAddress("centrality_impactparam", &centrality_impactparam);
124  // The data DST doesn't have centrality_mbd and centrality_mbdquantity branches
126  t->SetBranchAddress("centrality_mbd", &centrality_mbd);
127  t->SetBranchAddress("centrality_mbdquantity", &centrality_mbdquantity);
128  }
129  // t->SetBranchAddress("centrality_mbd", &centrality_mbd);
130  // t->SetBranchAddress("centrality_mbdquantity", &centrality_mbdquantity);
131  t->SetBranchAddress("ClusLayer", &ClusLayer);
132  t->SetBranchAddress("ClusX", &ClusX);
133  t->SetBranchAddress("ClusY", &ClusY);
134  t->SetBranchAddress("ClusZ", &ClusZ);
135  t->SetBranchAddress("ClusPhiSize", &ClusPhiSize);
136  t->SetBranchAddress("ClusZSize", &ClusZSize);
137  t->SetBranchAddress("ClusLadderZId", &ClusLadderZId);
138  t->SetBranchAddress("ClusLadderPhiId", &ClusLadderPhiId);
139 
140  TFile *outfile = new TFile(outfilename.Data(), "RECREATE");
141  TTree *minitree = new TTree("minitree", "Minitree of reconstructed vertices");
142  SetVtxMinitree(minitree, vtxdata);
143 
144  for (int ev = 0; ev < NevtToRun; ev++)
145  {
146  Long64_t local = t->LoadTree(index->GetIndex()[ev]);
147  t->GetEntry(local);
148  // cout << "event = " << event << " local = " << local << endl;
149  cout << "event=" << event << " has a total of " << ClusLayer->size() << " clusters" << endl;
150 
151  CleanVec(INTTlayer1);
152  CleanVec(INTTlayer2);
153 
154  for (size_t ihit = 0; ihit < ClusLayer->size(); ihit++)
155  {
156  if ((ClusLayer->at(ihit) < 3 || ClusLayer->at(ihit) > 6) || (ClusLadderZId->at(ihit) < 0 || ClusLadderZId->at(ihit) > 3))
157  {
158  cout << "Cluster (layer, ladderZId) = (" << ClusLayer->at(ihit) << "," << ClusLadderZId->at(ihit) << ") is not in the INTT acceptance. Exit and check." << endl;
159  exit(1);
160  }
161 
162  if (ClusPhiSize->at(ihit) >= 5)
163  continue;
164 
165  int layer = (ClusLayer->at(ihit) == 3 || ClusLayer->at(ihit) == 4) ? 0 : 1;
166 
167  Hit *hit = new Hit(ClusX->at(ihit), ClusY->at(ihit), ClusZ->at(ihit), avgVtxX, avgVtxY, 0., layer);
168  float edge = (ClusLadderZId->at(ihit) == 0 || ClusLadderZId->at(ihit) == 2) ? 0.8 : 0.5;
169  hit->SetEdge(ClusZ->at(ihit) - edge, ClusZ->at(ihit) + edge);
170 
171  if (layer == 0)
172  INTTlayer1.push_back(hit);
173  else
174  INTTlayer2.push_back(hit);
175  }
176 
177  cout << "# of clusters in inner layer (layer ID 3+4, after cluster phi size selection) = " << INTTlayer1.size() << ", outer layer (layer ID 5+6) = " << INTTlayer2.size() << endl;
178 
179  TH1F *hM_vtxzprojseg = new TH1F(Form("hM_vtxzprojseg_ev%d", ev), Form("hM_vtxzprojseg_ev%d", ev), 2800, -70, 70);
180  int goodpaircount = 0;
181  for (size_t i = 0; i < INTTlayer1.size(); i++)
182  {
183  for (size_t j = 0; j < INTTlayer2.size(); j++)
184  {
185  if (fabs(deltaPhi(INTTlayer1[i]->Phi(), INTTlayer2[j]->Phi())) > dPhi_cut)
186  continue;
187 
188  // find the equation connecting two hits in the X-Y plane
189  // y = mx + b
190  float m = (INTTlayer2[j]->posY() - INTTlayer1[i]->posY()) / (INTTlayer2[j]->posX() - INTTlayer1[i]->posX());
191  float b = INTTlayer1[i]->posY() - m * INTTlayer1[i]->posX();
192  // calculate the distance of closest approach from the line to (avgVtxX, avgVtxY)
193  float dca = fabs(m * avgVtxX - avgVtxY + b) / sqrt(m * m + 1);
194 
195  if (dca > dca_cut)
196  continue;
197 
198  float rho1 = sqrt(pow((INTTlayer1[i]->posX() - avgVtxX), 2) + pow((INTTlayer1[i]->posY() - avgVtxY), 2));
199  float rho2 = sqrt(pow((INTTlayer2[j]->posX() - avgVtxX), 2) + pow((INTTlayer2[j]->posY() - avgVtxY), 2));
200  float z = INTTlayer1[i]->posZ() - (INTTlayer2[j]->posZ() - INTTlayer1[i]->posZ()) / (rho2 - rho1) * rho1;
201  float edge1 = INTTlayer1[i]->Edge().first - (INTTlayer2[j]->Edge().second - INTTlayer1[i]->Edge().first) / (rho2 - rho1) * rho1;
202  float edge2 = INTTlayer1[i]->Edge().second - (INTTlayer2[j]->Edge().first - INTTlayer1[i]->Edge().second) / (rho2 - rho1) * rho1;
203 
204  if (debug)
205  {
206  cout << "----------" << endl
207  << "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. << ","
208  << INTTlayer2[j]->posX() * 10. << "," << INTTlayer2[j]->posY() * 10. << "," << INTTlayer2[j]->posZ() * 10. << ")" << endl;
209  cout << "Cluster [phi1 (deg),eta1,phi2(deg),eta2]=[" << INTTlayer1[i]->Phi() * radianstodegree << "," << INTTlayer1[i]->Eta() << "," << INTTlayer2[j]->Phi() * radianstodegree << "," << INTTlayer2[j]->Eta() << "]" << endl
210  << "delta phi (in degree) = " << fabs(deltaPhi(INTTlayer1[i]->Phi(), INTTlayer2[j]->Phi())) * radianstodegree << "-> pass dPhi selection (<" << dPhi_cut
211  << " rad)=" << ((fabs(deltaPhi(INTTlayer1[i]->Phi(), INTTlayer2[j]->Phi())) < dPhi_cut) ? 1 : 0) << endl;
212  cout << "DCA cut = " << dca_cut << " [cm]; vertex candidate (center,edge1,edge2) = (" << z << "," << edge1 << "," << edge2 << "), difference = " << edge2 - edge1 << endl;
213  cout << "dca w.r.t (avgVtxX [cm], avgVtxY [cm]) (in mm) = " << dca * 10 << endl;
214  }
215 
216  goodpaircount++;
217  if (fabs(z) < 70)
218  {
219  int bin1 = hM_vtxzprojseg->FindBin(edge1);
220  int bin2 = hM_vtxzprojseg->FindBin(edge2);
221 
222  for (int ibin = bin1; ibin <= bin2; ibin++)
223  {
224  float bincenter = hM_vtxzprojseg->GetBinCenter(ibin);
225  float binwidth = hM_vtxzprojseg->GetBinWidth(ibin);
226  // for bin1 and bin2, find the distance of the edge1 and edge2 to the respective bin edge
227  float w;
228  if (ibin == bin1)
229  w = ((bincenter + 0.5 * binwidth) - edge1) / binwidth;
230  else if (ibin == bin2)
231  w = (edge2 - (bincenter - 0.5 * binwidth)) / binwidth;
232  else
233  w = 1.;
234 
235  hM_vtxzprojseg->Fill(bincenter, w);
236  // cout << "bincenter = " << bincenter << ", bincontent = " << bincontent << endl;
237  }
238  }
239  }
240  }
241 
242  cout << "Number of entries in hM_vtxzprojseg = " << hM_vtxzprojseg->Integral(-1, -1) << endl;
243  cout << "Number of good pairs = " << goodpaircount << endl;
244  if (debug && !IsData)
245  cout << "Event " << ev << " (Truth PVx, Truth PVy, Truth PVz) = (" << TruthPV_trig_x << ", " << TruthPV_trig_y << ", " << TruthPV_trig_z << ")" << endl;
246 
247  // find the maximum bin of the histogram hM_vtxzedge_dca[2] and fit the histogram with a Gaussian function around the maximum bin
248  float maxbincenter = hM_vtxzprojseg->GetBinCenter(hM_vtxzprojseg->GetMaximumBin());
249  // TF1 *f1 = new TF1("f1", "[0]*exp(-0.5*((x-[1])/[2])^2) + [3] + [4]*x + [5]*x*x + [6]*x*x*x", maxbincenter - 9, maxbincenter + 9);
250  // f1->SetParName(0, "Norm");
251  // f1->SetParName(1, "#mu");
252  // f1->SetParName(2, "#sigma");
253  // f1->SetParName(3, "p0");
254  // f1->SetParName(4, "p1");
255  // f1->SetParName(5, "p2");
256  // f1->SetParName(6, "p3");
257  // f1->SetParameter(0, hM_vtxzprojseg->GetMaximum());
258  // f1->SetParameter(1, maxbincenter);
259  // f1->SetParameter(2, 5);
260  // f1->SetParLimits(0, hM_vtxzprojseg->GetMaximum() * 0.01, hM_vtxzprojseg->GetMaximum() * 100);
261  // f1->SetParLimits(1, maxbincenter - 10, maxbincenter + 10);
262  // f1->SetParLimits(2, 0.1, 100);
263  // hM_vtxzprojseg->Fit("f1", "R");
264  // float mean = f1->GetParameter(1);
265  // float meanErr = f1->GetParError(1);
266 
267  TF1 *f1 = new TF1("gaus_fit", gaus_func, maxbincenter - 10, maxbincenter + 10, 4); // Gaussian + const. offset
268  f1->SetParameters(hM_vtxzprojseg->GetBinContent(hM_vtxzprojseg->GetMaximumBin()), hM_vtxzprojseg->GetBinCenter(hM_vtxzprojseg->GetMaximumBin()), 4, 0);
269  f1->SetParLimits(0, 0, 10000);
270  f1->SetParLimits(2, 0, 1000);
271  f1->SetParLimits(3, 0, 1000);
272  hM_vtxzprojseg->Fit(f1, "NQ", "", hM_vtxzprojseg->GetBinCenter(hM_vtxzprojseg->GetMaximumBin()) - 9, hM_vtxzprojseg->GetBinCenter(hM_vtxzprojseg->GetMaximumBin()) + 9);
273  f1->SetParameters(hM_vtxzprojseg->GetBinContent(hM_vtxzprojseg->GetMaximumBin()), hM_vtxzprojseg->GetBinCenter(hM_vtxzprojseg->GetMaximumBin()), 4, 0);
274  f1->SetParLimits(0, 0, 10000);
275  f1->SetParLimits(2, 0, 1000);
276  f1->SetParLimits(3, -20, 1000);
277  hM_vtxzprojseg->Fit(f1, "NQ", "", hM_vtxzprojseg->GetBinCenter(hM_vtxzprojseg->GetMaximumBin()) - 9, hM_vtxzprojseg->GetBinCenter(hM_vtxzprojseg->GetMaximumBin()) + 9);
278  float mean = f1->GetParameter(1);
279  float meanErr = f1->GetParError(1);
280 
281  if (debug)
282  cout << "Event " << ev << " Reco PVz = " << mean << " +/- " << meanErr << " cm" << endl;
283 
284  if (ev < NevtToPlot && makedemoplot)
285  {
286  TString pname = Form("%s/hM_vtxzprojseg_ev%d", demoplotpath.Data(), event);
287  draw_demoplot(hM_vtxzprojseg, f1, dca_cut, pname);
288  }
289 
290  vtxdata.isdata = IsData;
291  vtxdata.event = event;
292  if (!IsData)
293  {
294  vtxdata.NTruthVtx = NTruthVtx;
295  vtxdata.TruthPV_x = TruthPV_trig_x;
296  vtxdata.TruthPV_y = TruthPV_trig_y;
297  vtxdata.TruthPV_z = TruthPV_trig_z;
298  vtxdata.Centrality_bimp = centrality_bimp;
299  vtxdata.Centrality_impactparam = centrality_impactparam;
300  }
301  vtxdata.NClusLayer1 = (int)INTTlayer1.size();
302  vtxdata.PV_x = avgVtxX;
303  vtxdata.PV_y = avgVtxY;
304  vtxdata.PV_z = mean;
305  vtxdata.Centrality_mbd = centrality_mbd;
306  vtxdata.Centrality_mbdquantity = centrality_mbdquantity;
307 
308  minitree->Fill();
309  }
310 
311  outfile->cd();
312  minitree->Write();
313  outfile->Close();
314 
315  f->Close();
316 
317  return 0;
318 }