Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
plotRecoCluster_INTT.cxx
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file plotRecoCluster_INTT.cxx
1 #include <TCanvas.h>
2 #include <TCut.h>
3 #include <TFile.h>
4 #include <TH1F.h>
5 #include <TH2F.h>
6 #include <TMath.h>
7 #include <TObjString.h>
8 #include <TRandom3.h>
9 #include <TTree.h>
10 #include <TTreeIndex.h>
11 
12 #include <fstream>
13 #include <iostream>
14 #include <sstream>
15 #include <string>
16 #include <vector>
17 
18 #include "GenHadron.h"
19 #include "Tracklet.h"
20 #include "Vertex.h"
21 
22 int main(int argc, char *argv[])
23 {
24  bool rundata = true;
25  TString infilename = "/sphenix/user/hjheng/TrackletAna/data/INTT/INTTRecoClusters_test.root";
26  TString outfilename = Form("/sphenix/user/hjheng/TrackletAna/analysis/plot/hists/INTTData/Hists_RecoClusters_INTT.root");
27  // TString EvtVtx_map_filename = "/sphenix/user/hjheng/TrackletAna/minitree/AuAu_ana325private_NoPileup_RecoVtx_Optimization/TrackletAna_RecoClusters_RecoVtx_TklCluster_dPhiCutbin3_dZCutbin10.root";
28 
29  // std::map<int, float> EvtVtx_map = EvtVtx_map_tklcluster(EvtVtx_map_filename.Data());
30 
31  TH1F *hM_ClusX_all = new TH1F("hM_ClusX_all", "hM_ClusX_all", 200, -10, 10);
32  TH1F *hM_ClusX_layer1 = new TH1F("hM_ClusX_layer1", "hM_ClusX_layer1", 200, -10, 10);
33  TH1F *hM_ClusX_layer2 = new TH1F("hM_ClusX_layer2", "hM_ClusX_layer2", 200, -10, 10);
34  TH1F *hM_ClusY_all = new TH1F("hM_ClusY_all", "hM_ClusY_all", 200, -10, 10);
35  TH1F *hM_ClusY_layer1 = new TH1F("hM_ClusY_layer1", "hM_ClusY_layer1", 200, -10, 10);
36  TH1F *hM_ClusY_layer2 = new TH1F("hM_ClusY_layer2", "hM_ClusY_layer2", 200, -10, 10);
37  TH1F *hM_ClusZ_all = new TH1F("hM_ClusZ_all", "hM_ClusZ_all", 200, -25, 25);
38  TH1F *hM_ClusZ_layer1 = new TH1F("hM_ClusZ_layer1", "hM_ClusZ_layer1", 200, -25, 25);
39  TH1F *hM_ClusZ_layer2 = new TH1F("hM_ClusZ_layer2", "hM_ClusZ_layer2", 200, -25, 25);
40  TH1F *hM_ClusR_all = new TH1F("hM_ClusR_all", "hM_ClusR_all", 120, 4, 7);
41  TH1F *hM_ClusR_layer1 = new TH1F("hM_ClusR_layer1", "hM_ClusR_layer1", 120, 4, 7);
42  TH1F *hM_ClusR_layer2 = new TH1F("hM_ClusR_layer2", "hM_ClusR_layer2", 120, 4, 7);
43  // TH1F *hM_ClusEtaPV_all = new TH1F("hM_ClusEtaPV_all", "hM_ClusEtaPV_all", 160, -4, 4);
44  // TH1F *hM_ClusEtaPV_layer1 = new TH1F("hM_ClusEtaPV_layer1", "hM_ClusEtaPV_layer1", 160, -4, 4);
45  // TH1F *hM_ClusEtaPV_layer2 = new TH1F("hM_ClusEtaPV_layer2", "hM_ClusEtaPV_layer2", 160, -4, 4);
46  TH1F *hM_ClusEtaOri_all = new TH1F("hM_ClusEtaOri_all", "hM_ClusEtaOri_all", 120, -3, 3);
47  TH1F *hM_ClusEtaOri_layer1 = new TH1F("hM_ClusEtaOri_layer1", "hM_ClusEtaOri_layer1", 120, -3, 3);
48  TH1F *hM_ClusEtaOri_layer2 = new TH1F("hM_ClusEtaOri_layer2", "hM_ClusEtaOri_layer2", 120, -3, 3);
49  TH1F *hM_ClusPhi_all = new TH1F("hM_ClusPhi_all", "hM_ClusPhi_all", 140, -3.5, 3.5);
50  TH1F *hM_ClusPhi_layer1 = new TH1F("hM_ClusPhi_layer1", "hM_ClusPhi_layer1", 140, -3.5, 3.5);
51  TH1F *hM_ClusPhi_layer2 = new TH1F("hM_ClusPhi_layer2", "hM_ClusPhi_layer2", 140, -3.5, 3.5);
52  TH1F *hM_ClusADC_all = new TH1F("hM_ClusADC_all", "hM_ClusADC_all", 100, 0, 100);
53  TH1F *hM_ClusADC_layer1 = new TH1F("hM_ClusADC_layer1", "hM_ClusADC_layer1", 100, 0, 100);
54  TH1F *hM_ClusADC_layer2 = new TH1F("hM_ClusADC_layer2", "hM_ClusADC_layer2", 100, 0, 100);
55  TH1F *hM_ClusZSize_all = new TH1F("hM_ClusZSize_all", "hM_ClusZSize_all", 20, 0, 20);
56  TH1F *hM_ClusZSize_layer1 = new TH1F("hM_ClusZSize_layer1", "hM_ClusZSize_layer1", 20, 0, 20);
57  TH1F *hM_ClusZSize_layer2 = new TH1F("hM_ClusZSize_layer2", "hM_ClusZSize_layer2", 20, 0, 20);
58  TH1F *hM_ClusPhiSize_all = new TH1F("hM_ClusPhiSize_all", "hM_ClusPhiSize_all", 40, 0, 40);
59  TH1F *hM_ClusPhiSize_layer1 = new TH1F("hM_ClusPhiSize_layer1", "hM_ClusPhiSize_layer1", 40, 0, 40);
60  TH1F *hM_ClusPhiSize_layer2 = new TH1F("hM_ClusPhiSize_layer2", "hM_ClusPhiSize_layer2", 40, 0, 40);
61  // TH1F *hM_ClusPairDr_layer1 = new TH1F("hM_ClusPairDr_layer1", "hM_ClusPairDr_layer1", 50, 0, 0.1);
62  // TH1F *hM_ClusPairDr_layer2 = new TH1F("hM_ClusPairDr_layer2", "hM_ClusPairDr_layer2", 50, 0, 0.1);
63  TH2F *hM_ClusX_ClusY_all = new TH2F("hM_ClusX_ClusY_all", "hM_ClusX_ClusY_all", 260, -13, 13, 260, -13, 13);
64  TH2F *hM_ClusZ_ClusPhi_all = new TH2F("hM_ClusZ_ClusPhi_all", "hM_ClusZ_ClusPhi_all", 1000, -25, 25, 350, -3.5, 3.5);
65  TH2F *hM_ClusZ_ClusPhi_layer1 = new TH2F("hM_ClusZ_ClusPhi_layer1", "hM_ClusZ_ClusPhi_layer1", 1000, -25, 25, 350, -3.5, 3.5);
66  TH2F *hM_ClusZ_ClusPhi_layer2 = new TH2F("hM_ClusZ_ClusPhi_layer2", "hM_ClusZ_ClusPhi_layer2", 1000, -25, 25, 350, -3.5, 3.5);
67  TH2F *hM_ClusEta_ClusZSize_all = new TH2F("hM_ClusEta_ClusZSize_all", "hM_ClusEta_ClusZSize_all", 160, -4, 4, 20, 0, 20);
68  TH2F *hM_ClusEta_ClusZSize_layer1 = new TH2F("hM_ClusEta_ClusZSize_layer1", "hM_ClusEta_ClusZSize_layer1", 160, -4, 4, 20, 0, 20);
69  TH2F *hM_ClusEta_ClusZSize_layer2 = new TH2F("hM_ClusEta_ClusZSize_layer2", "hM_ClusEta_ClusZSize_layer2", 160, -4, 4, 20, 0, 20);
70  TH2F *hM_ClusPhi_ClusPhiSize_all = new TH2F("hM_ClusPhi_ClusPhiSize_all", "hM_ClusPhi_ClusPhiSize_all", 140, -3.5, 3.5, 20, 0, 20);
71  TH2F *hM_ClusPhi_ClusPhiSize_layer1 = new TH2F("hM_ClusPhi_ClusPhiSize_layer1", "hM_ClusPhi_ClusPhiSize_layer1", 140, -3.5, 3.5, 20, 0, 20);
72  TH2F *hM_ClusPhi_ClusPhiSize_layer2 = new TH2F("hM_ClusPhi_ClusPhiSize_layer2", "hM_ClusPhi_ClusPhiSize_layer2", 140, -3.5, 3.5, 20, 0, 20);
73  TH2F *hM_ClusZSize_ClusPhiSize_all = new TH2F("hM_ClusZSize_ClusPhiSize_all", "hM_ClusZSize_ClusPhiSize_all", 20, 0, 20, 20, 0, 20);
74  TH2F *hM_ClusZSize_ClusPhiSize_layer1 = new TH2F("hM_ClusZSize_ClusPhiSize_layer1", "hM_ClusZSize_ClusPhiSize_layer1", 20, 0, 20, 20, 0, 20);
75  TH2F *hM_ClusZSize_ClusPhiSize_layer2 = new TH2F("hM_ClusZSize_ClusPhiSize_layer2", "hM_ClusZSize_ClusPhiSize_layer2", 20, 0, 20, 20, 0, 20);
76  TH2F *hM_ClusEta_ClusADC_all = new TH2F("hM_ClusEta_ClusADC_all", "hM_ClusEta_ClusADC_all", 160, -4, 4, 50, 0, 50);
77  TH2F *hM_ClusEta_ClusADC_layer1 = new TH2F("hM_ClusEta_ClusADC_layer1", "hM_ClusEta_ClusADC_layer1", 160, -4, 4, 50, 0, 50);
78  TH2F *hM_ClusEta_ClusADC_layer2 = new TH2F("hM_ClusEta_ClusADC_layer2", "hM_ClusEta_ClusADC_layer2", 160, -4, 4, 50, 0, 50);
79  TH2F *hM_ClusPhiSize_ClusADC_all = new TH2F("hM_ClusPhiSize_ClusADC_all", "hM_ClusPhiSize_ClusADC_all", 40, 0, 40, 100, 0, 100);
80  TH2F *hM_ClusPhiSize_ClusADC_layer1 = new TH2F("hM_ClusPhiSize_ClusADC_layer1", "hM_ClusPhiSize_ClusADC_layer1", 40, 0, 40, 100, 0, 100);
81  TH2F *hM_ClusPhiSize_ClusADC_layer2 = new TH2F("hM_ClusPhiSize_ClusADC_layer2", "hM_ClusPhiSize_ClusADC_layer2", 40, 0, 40, 100, 0, 100);
82 
83  TFile *f = new TFile(infilename, "READ");
84  TTree *t = (TTree *)f->Get("EventTree");
85  t->BuildIndex("event"); // Reference: https://root-forum.cern.ch/t/sort-ttree-entries/13138
86  TTreeIndex *index = (TTreeIndex *)t->GetTreeIndex();
87  int event, NTruthVtx;
88  float TruthPV_trig_z, TruthPV_mostNpart_z;
89  vector<int> *ClusLayer = 0;
90  vector<float> *ClusX = 0, *ClusY = 0, *ClusZ = 0, *ClusR = 0, *ClusPhi = 0, *ClusEta = 0, *ClusPhiSize = 0, *ClusZSize = 0;
91  vector<unsigned int> *ClusAdc = 0;
92  t->SetBranchAddress("event", &event);
93  t->SetBranchAddress("NTruthVtx", &NTruthVtx);
94  t->SetBranchAddress("TruthPV_trig_z", &TruthPV_trig_z);
95  t->SetBranchAddress("TruthPV_mostNpart_z", &TruthPV_mostNpart_z);
96  t->SetBranchAddress("ClusLayer", &ClusLayer);
97  t->SetBranchAddress("ClusX", &ClusX);
98  t->SetBranchAddress("ClusY", &ClusY);
99  t->SetBranchAddress("ClusZ", &ClusZ);
100  t->SetBranchAddress("ClusR", &ClusR);
101  t->SetBranchAddress("ClusPhi", &ClusPhi);
102  t->SetBranchAddress("ClusEta", &ClusEta);
103  t->SetBranchAddress("ClusAdc", &ClusAdc);
104  t->SetBranchAddress("ClusPhiSize", &ClusPhiSize);
105  t->SetBranchAddress("ClusZSize", &ClusZSize);
106  for (int i = 0; i < t->GetEntriesFast(); i++)
107  {
108  Long64_t local = t->LoadTree(index->GetIndex()[i]);
109  t->GetEntry(local);
110 
111  // float PV_z = EvtVtx_map[event];
112  float PV_z = -20.; // FIX!
113 
114  // vector<Hit *> l1hits, l2hits, l3hits;
115  // l1hits.clear();
116  // l2hits.clear();
117  // l3hits.clear();
118 
119  for (size_t i = 0; i < ClusLayer->size(); i++)
120  {
121  // float x0 = ClusX->at(i);
122  // float y0 = ClusY->at(i);
123  // float z0 = ClusZ->at(i);
124  // vector<float> cpos = {x0, y0, z0};
125  // UpdatePos_GapTwoHalves(cpos, gap_north, cent_shift, gap_upper);
126  // Hit *hit = new Hit(cpos[0], cpos[1], cpos[2], 0., 0., PV_z, ClusLayer->at(i));
127 
128  hM_ClusX_all->Fill(ClusX->at(i));
129  hM_ClusY_all->Fill(ClusY->at(i));
130  hM_ClusZ_all->Fill(ClusZ->at(i));
131  hM_ClusR_all->Fill(ClusR->at(i));
132  // hM_ClusEtaPV_all->Fill(ClusEta->at(i));
133  hM_ClusEtaOri_all->Fill(ClusEta->at(i));
134  hM_ClusPhi_all->Fill(ClusPhi->at(i));
135  hM_ClusADC_all->Fill(ClusAdc->at(i));
136  // hM_ClusZSize_all->Fill(ClusZSize->at(i));
137  hM_ClusPhiSize_all->Fill(ClusPhiSize->at(i));
138  hM_ClusX_ClusY_all->Fill(ClusX->at(i), ClusY->at(i));
139  hM_ClusZ_ClusPhi_all->Fill(ClusZ->at(i), ClusPhi->at(i));
140  // hM_ClusEta_ClusZSize_all->Fill(ClusEta->at(i), ClusZSize->at(i));
141  hM_ClusPhi_ClusPhiSize_all->Fill(ClusPhi->at(i), ClusPhiSize->at(i));
142  // hM_ClusZSize_ClusPhiSize_all->Fill(ClusZSize->at(i), ClusPhiSize->at(i));
143  hM_ClusEta_ClusADC_all->Fill(ClusEta->at(i), ClusAdc->at(i));
144  // hM_ClusZSize_ClusADC_all->Fill(ClusZSize->at(i), ClusAdc->at(i));
145  hM_ClusPhiSize_ClusADC_all->Fill(ClusPhiSize->at(i), ClusAdc->at(i));
146  if (ClusLayer->at(i) == 3 || ClusLayer->at(i) == 4)
147  {
148  // l1hits.push_back(hit);
149  hM_ClusX_layer1->Fill(ClusX->at(i));
150  hM_ClusY_layer1->Fill(ClusY->at(i));
151  hM_ClusZ_layer1->Fill(ClusZ->at(i));
152  hM_ClusR_layer1->Fill(ClusR->at(i));
153  // hM_ClusEtaPV_layer1->Fill(ClusEta->at(i));
154  hM_ClusEtaOri_layer1->Fill(ClusEta->at(i));
155  hM_ClusPhi_layer1->Fill(ClusPhi->at(i));
156  hM_ClusADC_layer1->Fill(ClusAdc->at(i));
157  // hM_ClusZSize_layer1->Fill(ClusZSize->at(i));
158  hM_ClusPhiSize_layer1->Fill(ClusPhiSize->at(i));
159  hM_ClusZ_ClusPhi_layer1->Fill(ClusZ->at(i), ClusPhi->at(i));
160  // hM_ClusEta_ClusZSize_layer1->Fill(ClusEta->at(i), ClusZSize->at(i));
161  hM_ClusPhi_ClusPhiSize_layer1->Fill(ClusPhi->at(i), ClusPhiSize->at(i));
162  // hM_ClusZSize_ClusPhiSize_layer1->Fill(ClusZSize->at(i), ClusPhiSize->at(i));
163  hM_ClusEta_ClusADC_layer1->Fill(ClusEta->at(i), ClusAdc->at(i));
164  // hM_ClusZSize_ClusADC_layer1->Fill(ClusZSize->at(i), ClusAdc->at(i));
165  hM_ClusPhiSize_ClusADC_layer1->Fill(ClusPhiSize->at(i), ClusAdc->at(i));
166  }
167  else if (ClusLayer->at(i) == 5 || ClusLayer->at(i) == 6)
168  {
169  // l2hits.push_back(hit);
170  hM_ClusX_layer2->Fill(ClusX->at(i));
171  hM_ClusY_layer2->Fill(ClusY->at(i));
172  hM_ClusZ_layer2->Fill(ClusZ->at(i));
173  hM_ClusR_layer2->Fill(ClusR->at(i));
174  // hM_ClusEtaPV_layer2->Fill(ClusEta->at(i));
175  hM_ClusEtaOri_layer2->Fill(ClusEta->at(i));
176  hM_ClusPhi_layer2->Fill(ClusPhi->at(i));
177  hM_ClusADC_layer2->Fill(ClusAdc->at(i));
178  // hM_ClusZSize_layer2->Fill(ClusZSize->at(i));
179  hM_ClusPhiSize_layer2->Fill(ClusPhiSize->at(i));
180  hM_ClusZ_ClusPhi_layer2->Fill(ClusZ->at(i), ClusPhi->at(i));
181  // hM_ClusEta_ClusZSize_layer2->Fill(ClusEta->at(i), ClusZSize->at(i));
182  hM_ClusPhi_ClusPhiSize_layer2->Fill(ClusPhi->at(i), ClusPhiSize->at(i));
183  // hM_ClusZSize_ClusPhiSize_layer2->Fill(ClusZSize->at(i), ClusPhiSize->at(i));
184  hM_ClusEta_ClusADC_layer2->Fill(ClusEta->at(i), ClusAdc->at(i));
185  // hM_ClusZSize_ClusADC_layer2->Fill(ClusZSize->at(i), ClusAdc->at(i));
186  hM_ClusPhiSize_ClusADC_layer2->Fill(ClusPhiSize->at(i), ClusAdc->at(i));
187  }
188  else
189  {
190  cout << "[WARNING] ClusLayer = " << ClusLayer->at(i) << ", which is not INTT. Check!" << endl;
191  continue;
192  }
193  }
194  }
195 
196  TFile *fout = new TFile(outfilename, "RECREATE");
197  fout->cd();
198  hM_ClusX_all->Write();
199  hM_ClusX_layer1->Write();
200  hM_ClusX_layer2->Write();
201  hM_ClusY_all->Write();
202  hM_ClusY_layer1->Write();
203  hM_ClusY_layer2->Write();
204  hM_ClusZ_all->Write();
205  hM_ClusZ_layer1->Write();
206  hM_ClusZ_layer2->Write();
207  hM_ClusR_all->Write();
208  hM_ClusR_layer1->Write();
209  hM_ClusR_layer2->Write();
210  // hM_ClusEtaPV_all->Write();
211  // hM_ClusEtaPV_layer1->Write();
212  // hM_ClusEtaPV_layer2->Write();
213  hM_ClusEtaOri_all->Write();
214  hM_ClusEtaOri_layer1->Write();
215  hM_ClusEtaOri_layer2->Write();
216  hM_ClusPhi_all->Write();
217  hM_ClusPhi_layer1->Write();
218  hM_ClusPhi_layer2->Write();
219  hM_ClusADC_all->Write();
220  hM_ClusADC_layer1->Write();
221  hM_ClusADC_layer2->Write();
222  // hM_ClusZSize_all->Write();
223  // hM_ClusZSize_layer1->Write();
224  // hM_ClusZSize_layer2->Write();
225  hM_ClusPhiSize_all->Write();
226  hM_ClusPhiSize_layer1->Write();
227  hM_ClusPhiSize_layer2->Write();
228  // hM_ClusPairDr_layer1->Write();
229  // hM_ClusPairDr_layer2- >Write();
230  hM_ClusX_ClusY_all->Write();
231  hM_ClusZ_ClusPhi_all->Write();
232  hM_ClusZ_ClusPhi_layer1->Write();
233  hM_ClusZ_ClusPhi_layer2->Write();
234  // hM_ClusEta_ClusZSize_all->Write();
235  // hM_ClusEta_ClusZSize_layer1->Write();
236  // hM_ClusEta_ClusZSize_layer2->Write();
237  hM_ClusPhi_ClusPhiSize_all->Write();
238  hM_ClusPhi_ClusPhiSize_layer1->Write();
239  hM_ClusPhi_ClusPhiSize_layer2->Write();
240  hM_ClusEta_ClusADC_all->Write();
241  hM_ClusEta_ClusADC_layer1->Write();
242  hM_ClusEta_ClusADC_layer2->Write();
243  // hM_ClusZSize_ClusPhiSize_all->Write();
244  // hM_ClusZSize_ClusPhiSize_layer1->Write();
245  // hM_ClusZSize_ClusPhiSize_layer2->Write();
246  // hM_ClusZSize_ClusADC_all->Write();
247  // hM_ClusZSize_ClusADC_layer1->Write();
248  // hM_ClusZSize_ClusADC_layer2->Write();
249  hM_ClusPhiSize_ClusADC_all->Write();
250  hM_ClusPhiSize_ClusADC_layer1->Write();
251  hM_ClusPhiSize_ClusADC_layer2->Write();
252  fout->Close();
253 
254  system(Form("cd ./plot/; python plotRecoCluster.py -f %s -i %d", outfilename.Data(), rundata));
255 }