Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
TrackCaloDist.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file TrackCaloDist.C
1 #include <iostream>
2 #include <TH2D.h>
3 #include <TH1D.h>
4 #include <TChain.h>
5 #include <TMath.h>
6 #include <TTree.h>
7 #include <TFile.h>
8 #include <sstream> //std::ostringstsream
9 #include <fstream> //std::ifstream
10 #include <iostream> //std::cout, std::endl
11 #include <cmath>
12 #include <TGraphErrors.h>
13 #include <TGraph.h>
14 #include <TSpectrum.h>
15 #include <TF1.h>
16 #include <TTreeReader.h>
17 #include <TTreeReaderValue.h>
18 #include <TTreeReaderArray.h>
19 #include <string>
20 #include <set>
21 
22 using namespace std;
23 
24 void TrackCaloDist()
25 {
26  TFile *out = new TFile("TrackCaloDist.root","RECREATE");
27 
28  TH1D* h_dist[6]; // distance between track (eta,phi) and shower vertex (eta,phi)
29  TH1D* h_emcal[6]; // ratio of energy deposition in emcal to track momentum
30  TH1D* h_ihcal[6]; // ratio of energy deposition in ihcal to track momentum
31  TH1D* h_ohcal[6]; // ratio of energy deposition in ohcal to track momentum
32  TH2D* h_totale_dist[6]; // total energy deposition from the shower vs distance between track and shower
33 
34  // only filled in the case of an ohcal shower
35  TH1D* ohcal_top = new TH1D("ohcal_top","",300,0,3); // ratio of energy in leading tower of ohcal to track momentum
36  TH1D* ohcal_towers = new TH1D("ohcal_towers","",300,0,3); // ratio of energy in ohcal tower region to track mom
37  TH1D* ohcal_total = new TH1D("ohcal_total","",300,0,3); // ratio of energy in all calo tower region to track mom
38  TH1D* ohcal_ihcal_mip = new TH1D("ohcal_ihcal_mip","",300,0,1); // energy in ihcal in case of ohcal shower
39  TH1D* ohcal_emcal_mip = new TH1D("ohcal_emcal_mip","",300,0,1); // energy in emcal in case of ohcal shower
40  TH2D* h_4towers = new TH2D("h_4towers","",5,-2.5,2.5,5,-2.5,2.5); // 5x5 tower energy deposition in ohcal
41 
42  // classifications of shower location
43  string detector[6] = { "mip", "ohcal", "magnet", "ihcal", "emcal", "tracker"};
44  for (int i = 0; i < 6; i++) {
45  h_dist[i] = new TH1D(TString::Format("h_dist_%s",detector[i].c_str()),"",200,0,5);
46  h_emcal[i] = new TH1D(TString::Format("h_emcale_%s",detector[i].c_str()),"",100,0,3);
47  h_ihcal[i] = new TH1D(TString::Format("h_ihcale_%s",detector[i].c_str()),"",100,0,3);
48  h_ohcal[i] = new TH1D(TString::Format("h_ohcale_%s",detector[i].c_str()),"",100,0,3);
49  h_totale_dist[i] = new TH2D(TString::Format("h_totale_dist_%s",detector[i].c_str()),"",100,0,3,200,0,5);
50  }
51 
52  TH2D* h_ohcal_dist = new TH2D("h_ohcale_dist","",100,0,3,200,0,5);
53  TH1D* h_class = new TH1D("h_class","",7,-0.5,6.5);
54 
55  // distance between the track (eta,phi) and the shower vertex for each classification
56  // energy distribution for each calorimeter for each classification
57 
58  TFile *file = new TFile("TruthCaloTree.root");
59  TTreeReader reader("T",file);
60 
61  TTreeReaderValue<int> n(reader, "n_truth");
62  TTreeReaderArray<float> E(reader, "truthenergy");
63  TTreeReaderArray<int> pid(reader, "truthpid");
64  TTreeReaderArray<float> eta(reader, "trutheta");
65  TTreeReaderArray<float> phi(reader, "truthphi");
66  TTreeReaderArray<float> p(reader, "truthp");
67  TTreeReaderArray<float> pt(reader, "truthpt");
68  TTreeReaderArray<int> track_id(reader, "truth_track_id");
69 
70  TTreeReaderValue<int> n_child(reader, "n_child");
71  TTreeReaderArray<int> child_pid(reader, "child_pid");
72  TTreeReaderArray<int> child_parent_id(reader, "child_parent_id");
73  TTreeReaderArray<int> child_vertex_id(reader, "child_vertex_id");
74 
75  TTreeReaderValue<int> BH_n(reader, "BH_n");
76  TTreeReaderArray<int> BH_track_id(reader, "BH_track_id");
77  TTreeReaderArray<float> BH_e(reader, "BH_e");
78 
79  TTreeReaderValue<int> n_vertex(reader, "n_vertex");
80  TTreeReaderArray<int> vertex_id(reader, "vertex_id");
81  TTreeReaderArray<float> x(reader, "vertex_x");
82  TTreeReaderArray<float> y(reader, "vertex_y");
83  TTreeReaderArray<float> z(reader, "vertex_z");
84 
85  TTreeReaderArray<float> tower_sim_E(reader, "tower_sim_E");
86  TTreeReaderArray<int> tower_sim_ieta(reader, "tower_sim_ieta");
87  TTreeReaderArray<int> tower_sim_iphi(reader, "tower_sim_iphi");
88  TTreeReaderArray<float> tower_sim_eta(reader, "tower_sim_eta");
89  TTreeReaderArray<float> tower_sim_phi(reader, "tower_sim_phi");
90  TTreeReaderValue<int> tower_sim_n(reader, "tower_sim_n");
91 
92  TTreeReaderArray<float> EMcal_sim_E(reader, "EMcal_sim_E");
93  TTreeReaderArray<int> EMcal_sim_ieta(reader, "EMcal_sim_ieta");
94  TTreeReaderArray<int> EMcal_sim_iphi(reader, "EMcal_sim_iphi");
95  TTreeReaderArray<float> EMcal_sim_eta(reader, "EMcal_sim_eta");
96  TTreeReaderArray<float> EMcal_sim_phi(reader, "EMcal_sim_phi");
97  TTreeReaderValue<int> EMcal_sim_n(reader, "EMcal_sim_n");
98 
99  TTreeReaderArray<float> hcalIN_sim_E(reader, "hcalIN_sim_E");
100  TTreeReaderArray<int> hcalIN_sim_ieta(reader, "hcalIN_sim_ieta");
101  TTreeReaderArray<int> hcalIN_sim_iphi(reader, "hcalIN_sim_iphi");
102  TTreeReaderArray<float> hcalIN_sim_eta(reader, "hcalIN_sim_eta");
103  TTreeReaderArray<float> hcalIN_sim_phi(reader, "hcalIN_sim_phi");
104  TTreeReaderValue<int> hcalIN_sim_n(reader, "hcalIN_sim_n");
105 
106  const int nphi_hcal = 64; // number of hcal towers in phi direction
107  const int neta_hcal = 24;
108  const int nphi_emcal = 256;
109  const int neta_emcal = 96;
110  const float max_eta = 0.8; // maximum eta value for isolated track
111  const float neigh_max_eta = 1.0;
112  const float track_pt_cut = 0.5; // pt cut for neighboring particles
113  const float high_pt_cut = 4.0;
114  const float de = 0.045833; // distance from center of one calorimeter tower to edge in eta
115  const float dp = 0.049087; // distacne from center of one calorimeter tower to edge in phi
116  const float ihcal_sf = 0.162166; // sampling fractions
117  const float ohcal_sf = 0.0338021;
118  const float emcal_sf = 0.02;
119  const float emcal_inner = 92; // inner radius of emcal
120  const float emcal_outer = 116;
121  const float ihcal_inner = 117.27;
122  const float ihcal_outer = 134.42;
123  const float ohcal_inner = 183.3;
124  const float ohcal_outer = 274.317;
125 
126  double dr, delta_r, delta_eta, delta_phi, dr_charge, dphi, deta, tower_dr;
127  float mipEM, mipIH, mipOH;
128  double vr, v_phi, v_eta;
129  float avgtracks[10] = {0.0};
130  int ntracks_centrality[10] = {0};
131  int centrality_bins = 10;
132  set<int> charged_pid = {-3334,-3312,-3222,-3112,-2212,-431,-411,-321,-211,-13,-11,11,
133  13,211,321,411,431,2212,3112,3222,3312,3334};
134  set<int> neutral_pid = {-3322,-3212,-3122,-2112,-421,-14,-12,12,14,22,111,130,310,421,
135  2112,3122,3212,3322};
136  int cent[] = {40,50,60,70,80};
137 
138  float five_by_five[5][5] = {{0.0}};
139  int child;
140  set<int> vertex;
141  float v_radius, v_rad;
142  int classify;
143 
144  float eta_map[24],eta_mapIH[24],eta_mapEM[96];
145  float phi_map[64],phi_mapIH[64],phi_mapEM[256];
146 
147 
148  // first find the mapping of (ieta,iphi) calorimeter towers to the actual (eta,phi) coordinates of the towers
149  while (reader.Next()) {
150  for (int i = 0; i < *tower_sim_n; i++) {
151  if (tower_sim_ieta[i] >= 0 && tower_sim_ieta[i] < neta_hcal
152  && tower_sim_iphi[i] >= 0 && tower_sim_iphi[i] < nphi_hcal) {
153  phi_map[tower_sim_iphi[i]] = tower_sim_phi[i];
154  eta_map[tower_sim_ieta[i]] = tower_sim_eta[i];
155  }
156  }
157 
158  for (int i = 0; i < *hcalIN_sim_n; i++) {
159  if (hcalIN_sim_ieta[i] >= 0 && hcalIN_sim_ieta[i] < neta_hcal
160  && hcalIN_sim_iphi[i] >= 0 && hcalIN_sim_iphi[i] < nphi_hcal) {
161  phi_mapIH[hcalIN_sim_iphi[i]] = hcalIN_sim_phi[i];
162  eta_mapIH[hcalIN_sim_ieta[i]] = hcalIN_sim_eta[i];
163  }
164  }
165 
166  for (int i = 0; i < nphi_hcal; i++) {
167  if (phi_map[i] < 0.0) phi_map[i] += 2*M_PI;
168  if (phi_mapIH[i] < 0.0) phi_mapIH[i] += 2*M_PI;
169  }
170 
171  for (int i = 0; i < *EMcal_sim_n; i++) {
172  if (EMcal_sim_ieta[i] >= 0 && EMcal_sim_ieta[i] < neta_emcal
173  && EMcal_sim_iphi[i] >= 0 && EMcal_sim_iphi[i] < nphi_emcal) {
174  phi_mapEM[EMcal_sim_iphi[i]] = EMcal_sim_phi[i];
175  eta_mapEM[EMcal_sim_ieta[i]] = EMcal_sim_eta[i];
176  }
177  }
178 
179  for (int i = 0; i < nphi_emcal; i++) {
180  if (phi_mapEM[i] < 0.0) phi_mapEM[i] += 2*M_PI;
181  }
182  }
183 
184  reader.Restart();
185  int f = 0;
186  while (reader.Next()) {
187 
188  if ( f % 1000 == 0) cout << f << endl;
189  if (!vertex.empty()) vertex.clear();
190  classify = -1;
191 
192  float e_sim[24][64] = {{0.0}};
193  float e_simIH[24][64] = {{0.0}};
194  float e_simEM[96][256] = {{0.0}};
195  float total_energy = 0.0;
196 
197  // find the energy distribution in the calo towers for the event
198  for (int i = 0; i < *tower_sim_n; i++) {
199  if (tower_sim_ieta[i] >= 0 && tower_sim_ieta[i] < neta_hcal
200  && tower_sim_iphi[i] >= 0 && tower_sim_iphi[i] < nphi_hcal) {
201  e_sim[tower_sim_ieta[i]][tower_sim_iphi[i]] += tower_sim_E[i]/ohcal_sf;
202  }
203  }
204  for (int i = 0; i < *hcalIN_sim_n; i++) {
205  if (hcalIN_sim_ieta[i] >= 0 && hcalIN_sim_ieta[i] < neta_hcal
206  && hcalIN_sim_iphi[i] >= 0 && hcalIN_sim_iphi[i] < nphi_hcal) {
207  e_simIH[hcalIN_sim_ieta[i]][hcalIN_sim_iphi[i]] += hcalIN_sim_E[i]/ihcal_sf;
208  }
209  }
210  for (int i = 0; i < *EMcal_sim_n; i++) {
211  if (EMcal_sim_ieta[i] >= 0 && EMcal_sim_ieta[i] < neta_emcal
212  && EMcal_sim_iphi[i] >= 0 && EMcal_sim_iphi[i] < nphi_emcal) {
213  e_simEM[EMcal_sim_ieta[i]][EMcal_sim_iphi[i]] += EMcal_sim_E[i]/emcal_sf;
214  }
215  }
216 
217  // loop over truth level particles
218  for (int i = 0; i < *n; i++) {
219 
220  // find charged isolated high pt particle and check that it is isolated
221  if (pt[i] >= high_pt_cut && charged_pid.find(pid[i])!=charged_pid.end() && abs(double(eta[i])) < max_eta) {
222  if (phi[i] < 0.0) phi[i] += 2*M_PI;
223  dr = 100.0;
224  for (int j = 0; j < *n; j++) {
225  if (i != j && pt[j] > track_pt_cut && abs(double(eta[j])) < neigh_max_eta) {
226  if (phi[j] < 0.0) phi[j] += 2*M_PI;
227  delta_phi = double(phi[i]) - double(phi[j]);
228  if (delta_phi > M_PI) delta_phi -= M_PI;
229  delta_eta = double(eta[i]) - double(eta[j]);
230  delta_r = sqrt(delta_phi*delta_phi + delta_eta*delta_eta);
231  if (delta_r < dr) dr = delta_r;
232  }
233  }
234 
235  // check where particle begins showering
236  if (dr >= 0.2) {
237  for (int j = 0; j < *n_child; j++) {
238  if (child_parent_id[j] == track_id[i] && abs(child_pid[j]) != 11 && child_pid[j] != 22) vertex.insert(child_vertex_id[j]);
239  }
240  v_radius = ohcal_outer + 1;
241  for (set<int>::iterator it = vertex.begin(); it != vertex.end(); it++) {
242  child = 0;
243  for (int j = 0 ; j < *n_child; j++) {
244  if (child_vertex_id[j] == *it && abs(child_pid[j]) != 11 && child_pid[j] != 22) {
245  child++;
246  }
247  }
248  if (child > 1) {
249  // now find the location of the shower vertex
250  for (int v = 0; v < *n_vertex; v++) {
251  if (vertex_id[v] == *it) {
252  v_rad = sqrt(x[v]*x[v] + y[v]*y[v]);
253  if (v_rad < v_radius) {
254  v_radius = v_rad;
255  v_phi = atan2(y[v],x[v]);
256  v_eta = atanh(z[v] / sqrt(x[v]*x[v] + y[v]*y[v] + z[v]*z[v]));
257  dphi = abs(double(v_eta - eta[i]));
258  deta = abs(double(v_phi - phi[i]));
259  if (dphi >= M_PI) dphi -= M_PI;
260  vr = sqrt(dphi*dphi + deta*deta);
261  }
262  }
263  }
264  }
265  }
266 
267 
268  // find the energy deposition and number of towers in tower region
269  int towers[200][2];
270  float top_energy = 0.0;
271  int ti, tj;
272  int numOH = 0;
273  int tn = 0;
274  mipEM = 0.0;
275  mipIH = 0.0;
276  mipOH = 0.0;
277  for (int j = 0; j < neta_hcal; j++) {
278  for (int k = 0; k < nphi_hcal; k++) {
279  dphi = abs(double(eta_map[j] - eta[i]));
280  deta = abs(double(phi_map[k] - phi[i]));
281  if (dphi >= M_PI) dphi -= M_PI;
282  tower_dr = sqrt(dphi*dphi + deta*deta);
283  // find towers in (eta,phi) range of track
284  if (tower_dr <= 0.2) {
285  towers[tn][0] = j;
286  towers[tn][1] = k;
287  mipIH += e_simIH[j][k];
288  mipOH += e_sim[j][k];
289  tn++;
290  }
291  }
292  }
293 
294  for (int j = 0; j < neta_emcal; j++) {
295  for(int k = 0; k < nphi_emcal; k++) {
296  dphi = abs(double(eta_mapEM[j] - eta[i]));
297  deta = abs(double(phi_mapEM[k] - phi[i]));
298  if (dphi >= M_PI) dphi -= M_PI;
299  tower_dr = sqrt(dphi*dphi + deta*deta);
300  if (tower_dr <= 0.2) {
301  mipEM += e_simEM[j][k];
302  }
303  }
304  }
305 
306  int j;
307  int k;
308  double mipTotal = mipEM + mipIH + mipOH;
309 
310  for (int t = 0; t < tn; t++) {
311  j = towers[t][0];
312  k = towers[t][1];
313  if (e_sim[j][k] > top_energy) {
314  top_energy = e_sim[j][k];
315  ti = j;
316  tj = k;
317  }
318  }
319 
320  // based on shower location classification fill histograms
321  if (v_radius < emcal_inner) {
322  h_emcal[5]->Fill(mipEM/p[i]);
323  h_ihcal[5]->Fill(mipIH/p[i]);
324  h_ohcal[5]->Fill(mipOH/p[i]);
325  h_dist[5]->Fill(vr);
326  h_class->Fill(5);
327  h_totale_dist[5]->Fill(mipTotal/p[i],vr);
328  } else if (v_radius >= emcal_inner && v_radius <= emcal_outer) {
329  h_emcal[4]->Fill(mipEM/p[i]);
330  h_ihcal[4]->Fill(mipIH/p[i]);
331  h_ohcal[4]->Fill(mipOH/p[i]);
332  h_dist[4]->Fill(vr);
333  h_class->Fill(4);
334  h_totale_dist[4]->Fill(mipTotal/p[i],vr);
335  } else if (v_radius >= ihcal_inner && v_radius <= ihcal_outer) {
336  h_emcal[3]->Fill(mipEM/p[i]);
337  h_ihcal[3]->Fill(mipIH/p[i]);
338  h_ohcal[3]->Fill(mipOH/p[i]);
339  h_dist[3]->Fill(vr);
340  h_class->Fill(3);
341  h_totale_dist[3]->Fill(mipTotal/p[i],vr);
342  } else if (v_radius > ihcal_outer && v_radius < ohcal_inner) {
343  h_emcal[2]->Fill(mipEM/p[i]);
344  h_ihcal[2]->Fill(mipIH/p[i]);
345  h_ohcal[2]->Fill(mipOH/p[i]);
346  h_dist[2]->Fill(vr);
347  h_class->Fill(2);
348  h_totale_dist[2]->Fill(mipTotal/p[i],vr);
349  } else if (v_radius >= ohcal_inner && v_radius <= ohcal_outer) {
350  h_emcal[1]->Fill(mipEM/p[i]);
351  h_ihcal[1]->Fill(mipIH/p[i]);
352  h_ohcal[1]->Fill(mipOH/p[i]);
353  h_dist[1]->Fill(vr);
354  h_class->Fill(1);
355  h_ohcal_dist->Fill(mipOH/p[i],vr);
356  h_totale_dist[1]->Fill(mipTotal/p[i],vr);
357  ohcal_top->Fill(top_energy/p[i]);
358  ohcal_towers->Fill(mipOH/p[i]);
359  ohcal_total->Fill(mipTotal/p[i]);
360  ohcal_ihcal_mip->Fill(mipIH);
361  ohcal_emcal_mip->Fill(mipEM);
362 
363  // fill 5x5 tower energy distribution for ohcal showers
364  for (int fi = -2; fi < 3; fi++) {
365  for (int fj = -2; fj < 3; fj++) {
366  five_by_five[fi+2][fj+2] += e_sim[ti+fi][tj+fj];
367  h_4towers->SetBinContent(h_4towers->FindBin(fi,fj),five_by_five4[fi+2][fj+2]);
368  }
369  }
370  } else if (v_radius > ohcal_outer) {
371  // check if the event mips through the entire system
372  for (int j = 0; j < *BH_n; j++) {
373  if (BH_track_id[j] == track_id[i]) {
374  h_emcal[0]->Fill(mipEM/p[i]);
375  h_ihcal[0]->Fill(mipIH/p[i]);
376  h_ohcal[0]->Fill(mipOH/p[i]);
377  h_class->Fill(0);
378  h_totale_dist[0]->Fill(mipTotal/p[i],vr);
379  }
380  }
381  }
382  }
383  }
384  }
385  f++;
386  }
387  out->Write();
388  out->Close();
389 }