Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
event-display.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file event-display.C
1 #include <calobase/TowerInfoDefs.h>
2 
3 // c++ includes --
4 #include <string>
5 #include <vector>
6 #include <iostream>
7 #include <fstream>
8 
9 // -- root includes --
10 #include <TH2F.h>
11 #include <TFile.h>
12 #include <TTree.h>
13 #include <TLine.h>
14 #include <TCanvas.h>
15 
16 using std::cout;
17 using std::endl;
18 using std::string;
19 using std::vector;
20 using std::min;
21 using std::max;
22 using std::to_string;
23 using std::ofstream;
24 
25 R__LOAD_LIBRARY(libcalo_io.so)
26 
27 namespace myAnalysis {
28 
29  void init(const string& input_CEMC, const string& outputFile = "event-display/test.json", const string& input_HCALIN = "", const string& input_HCALOUT = "");
30  void analyze();
31  void finalize(const string& run, const string& event);
32 
33  ofstream output;
34 
35  TFile* finput_CEMC = 0;
36  TFile* finput_HCALIN = 0;
37  TFile* finput_HCALOUT = 0;
38 
39  TTree* tree_CEMC;
40  TTree* tree_HCALIN;
41  TTree* tree_HCALOUT;
42 
43  TH2F* h2ADC_CEMC;
44  TH2F* h2ADC_HCALIN;
45  TH2F* h2ADC_HCALOUT;
46 
47  TH1F* hADC_CEMC;
48  TH1F* hADC_HCALIN;
49  TH1F* hADC_HCALOUT;
50 
51  // CEMC
52 
53  UInt_t eta_bins = 96;
54  Float_t eta_min = -0.5;
55  Float_t eta_max = 95.5;
56 
57  UInt_t phi_bins = 256;
58  Float_t phi_min = -0.5;
59  Float_t phi_max = 255.5;
60 
61  UInt_t etabin_max = 95;
62  UInt_t phibin_max = 255;
63 
64  UInt_t channels_bins = 24576;
65  UInt_t channels_min = 0;
66  UInt_t channels_max = 24576;
67 
68  // HCAL
69 
70  UInt_t eta_hcal_bins = 24;
71  Float_t eta_hcal_min = -0.5;
72  Float_t eta_hcal_max = 23.5;
73 
74  UInt_t phi_hcal_bins = 64;
75  Float_t phi_hcal_min = -0.5;
76  Float_t phi_hcal_max = 63.5;
77 
78  UInt_t channels_hcal_bins = 1536;
79  UInt_t channels_hcal_min = 0;
80  UInt_t channels_hcal_max = 1536;
81 
82  const float eta_map[] = {-1.05417, -0.9625, -0.870833,-0.779167,-0.6875, -0.595833,-0.504167,-0.4125,
83  -0.320833,-0.229167,-0.1375, -0.0458333,0.0458333,0.1375, 0.229167, 0.320833,
84  0.4125, 0.504167, 0.595833, 0.6875, 0.779167, 0.870833, 0.9625, 1.05417};
85 
86  const float phi_map[] = {0.0490874, 0.147262, 0.245437, 0.343612, 0.441786, 0.539961, 0.638136, 0.736311,
87  0.834486, 0.93266, 1.03084, 1.12901, 1.22718, 1.32536, 1.42353, 1.52171,
88  1.61988, 1.71806, 1.81623, 1.91441, 2.01258, 2.11076, 2.20893, 2.30711,
89  2.40528, 2.50346, 2.60163, 2.69981, 2.79798, 2.89616, 2.99433, 3.09251,
90  3.19068, 3.28885, 3.38703, 3.4852, 3.58338, 3.68155, 3.77973, 3.8779,
91  3.97608, 4.07425, 4.17243, 4.2706, 4.36878, 4.46695, 4.56513, 4.6633,
92  4.76148, 4.85965, 4.95783, 5.056, 5.15418, 5.25235, 5.35052, 5.4487,
93  5.54687, 5.64505, 5.74322, 5.8414, 5.93957, 6.03775, 6.13592, 6.2341};
94 
95  UInt_t ADC_bins = 65;
96  Float_t ADC_min = 0;
97  Float_t ADC_max = 1.3e3;
98 }
99 
100 void myAnalysis::init(const string& input_CEMC, const string& outputFile, const string& input_HCALIN, const string& input_HCALOUT) {
101 
102  finput_CEMC = TFile::Open(input_CEMC.c_str());
103  if(!input_HCALIN.empty()) finput_HCALIN = TFile::Open(input_HCALIN.c_str());
104  if(!input_HCALOUT.empty()) finput_HCALOUT = TFile::Open(input_HCALOUT.c_str());
105 
106  tree_CEMC = (TTree*)finput_CEMC->Get("W");
107  tree_HCALIN = (TTree*)finput_HCALIN->Get("W");
108  tree_HCALOUT = (TTree*)finput_HCALOUT->Get("W");
109 
110  output.open(outputFile.c_str());
111 
112  hADC_CEMC = new TH1F("hADC_CEMC", "ADC CEMC; ADC; Counts", ADC_bins, ADC_min, ADC_max);
113  h2ADC_CEMC = new TH2F("h2ADC_CEMC", "ADC CEMC; towerid #eta; towerid #phi", eta_bins, eta_min, eta_max,
114  phi_bins, phi_min, phi_max);
115 
116  if(finput_HCALIN) {
117  hADC_HCALIN = new TH1F("hADC_HCALIN", "ADC HCALIN; ADC; Counts", ADC_bins, ADC_min, ADC_max);
118  h2ADC_HCALIN = new TH2F("h2ADC_HCALIN", "ADC HCALIN; towerid #eta; towerid #phi", eta_hcal_bins, eta_hcal_min, eta_hcal_max,
119  phi_hcal_bins, phi_hcal_min, phi_hcal_max);
120  }
121 
122  if(finput_HCALOUT) {
123  hADC_HCALOUT = new TH1F("hADC_HCALOUT", "ADC HCALOUT; ADC; Counts", ADC_bins, ADC_min, ADC_max);
124  h2ADC_HCALOUT = new TH2F("h2ADC_HCALOUT", "ADC HCALOUT; towerid #eta; towerid #phi", eta_hcal_bins, eta_hcal_min, eta_hcal_max,
125  phi_hcal_bins, phi_hcal_min, phi_hcal_max);
126  }
127 
128 output << "{\n \
129  \"EVENT\": {\n \
130  \"runid\": 1,\n \
131  \"evtid\": 1,\n \
132  \"time\": 0,\n \
133  \"type\": \"Collision\",\n \
134  \"s_nn\": 0,\n \
135  \"B\": 3.0,\n \
136  \"pv\": [0,0,0]\n \
137  },\n \
138  \"META\": {\n \
139  \"HITS\": {\n \
140  \"INNERTRACKER\": {\n \
141  \"type\": \"3D\",\n \
142  \"options\": {\n \
143  \"size\": 5,\n \
144  \"color\": 16777215\n \
145  }\n \
146  },\n \
147  \"CEMC\": {\n \
148  \"type\": \"PROJECTIVE\",\n \
149  \"options\": {\n \
150  \"rmin\": 90,\n \
151  \"rmax\": 116.1,\n \
152  \"deta\": 0.024,\n \
153  \"dphi\": 0.024,\n \
154  \"color\": 16766464,\n \
155  \"transparent\": 0.6,\n \
156  \"scaleminmax\": false\n \
157  }\n \
158  },\n \
159  \"HCALIN\": {\n \
160  \"type\": \"PROJECTIVE\",\n \
161  \"options\": {\n \
162  \"rmin\": 117.27,\n \
163  \"rmax\": 134.42,\n \
164  \"deta\": 0.025,\n \
165  \"dphi\": 0.025,\n \
166  \"color\": 4290445312,\n \
167  \"transparent\": 0.6,\n \
168  \"scaleminmax\": false\n \
169  }\n \
170  },\n \
171  \"HCALOUT\": {\n \
172  \"type\": \"PROJECTIVE\",\n \
173  \"options\": {\n \
174  \"rmin\": 183.3,\n \
175  \"rmax\": 274.317,\n \
176  \"deta\": 0.025,\n \
177  \"dphi\": 0.025,\n \
178  \"color\": 24773,\n \
179  \"transparent\": 0.6,\n \
180  \"scaleminmax\": true\n \
181  }\n \
182  }\n \
183  }\n \
184  },\n \
185 \"HITS\": {\n";
186 }
187 
188 void myAnalysis::analyze() {
189 
190  bool first_entry = true;
191  vector<Float_t>* ADC = 0;
192  vector<Int_t>* chan = 0;
193  // vector<Float_t>* time = 0;
194  // vector<Float_t>* ped = 0;
195  // vector<vector<Float_t>>* waveforms = 0; // 2D: channel x time sample
196  Float_t max_ADC = ADC_max;
197  if(finput_CEMC && tree_CEMC->GetEntries()) {
198  tree_CEMC->SetBranchAddress("adc",&ADC);
199  tree_CEMC->SetBranchAddress("chan",&chan);
200  tree_CEMC->GetEntry(0);
201  UInt_t nchannels = ADC->size();
202 
203  for(UInt_t j = 0; j < nchannels; ++j) {
204  // convert from high gain to low gain (factor of 16)
205  Float_t ADC_val = ADC->at(j)/16;
206  // Float_t e = 0.002*ADC_val/16;
207  Float_t e = ADC_val;
208  // if(ADC_val < 4e3) continue;
209  Int_t channel = chan->at(j);
210  if(channel >= channels_max) {
211  cout << "CEMC invalid channel number: " << channel << endl;
212  continue;
213  }
214 
215  UInt_t key = TowerInfoDefs::encode_emcal(channel);
216  UInt_t etabin = TowerInfoDefs::getCaloTowerEtaBin(key);
217  UInt_t phibin = TowerInfoDefs::getCaloTowerPhiBin(key);
218 
219  h2ADC_CEMC->SetBinContent(etabin+1, phibin+1, ADC_val);
220  hADC_CEMC->Fill(ADC_val);
221 
222  max_ADC = max(max_ADC, ADC_val);
223 
224  // mapping: [0,95] -> [-1.1, 1.1]
225  Float_t eta = 2.2/etabin_max*etabin-1.1;
226  // mapping: [0,127] -> [0, pi] and [128, 255] -> [-pi, 0]
227  Float_t phi = (phibin <= phibin_max/2) ? M_PI/(phibin_max/2)*phibin :
228  M_PI/(phibin_max/2)*(phibin*1.-phibin_max);
229 
230  if(first_entry) {
231  output << "\"CEMC\": [{ \"eta\": " << eta << ", \"phi\": " << phi << " , \"e\": " << e << "}\n";
232  first_entry = false;
233  }
234  else output << ",{ \"eta\": " << eta << ", \"phi\": " << phi << " , \"e\": " << e << "}\n";
235  // if(j == 10) break;
236  }
237  cout << "CEMC Max ADC: " << max_ADC << endl;
238  output << "],\n";
239  }
240  else output << "\"CEMC\": [{}],\n";
241 
242  first_entry = true;
243  max_ADC = ADC_max;
244  if(finput_HCALIN && tree_HCALIN->GetEntries()) {
245  ADC = 0;
246  chan = 0;
247 
248  tree_HCALIN->SetBranchAddress("adc",&ADC);
249  tree_HCALIN->SetBranchAddress("chan",&chan);
250  tree_HCALIN->GetEntry(0);
251 
252  UInt_t nchannels = ADC->size();
253 
254  for(UInt_t j = 0; j < nchannels; ++j) {
255  Float_t ADC_val = ADC->at(j);
256  // Float_t e = 0.002*ADC_val;
257  Float_t e = ADC_val;
258  Int_t channel = chan->at(j);
259  if(channel >= channels_max) {
260  cout << "HCALIN invalid channel number: " << channel << endl;
261  continue;
262  }
263 
264  UInt_t key = TowerInfoDefs::encode_hcal(channel);
265  UInt_t etabin = TowerInfoDefs::getCaloTowerEtaBin(key);
266  UInt_t phibin = TowerInfoDefs::getCaloTowerPhiBin(key);
267 
268  h2ADC_HCALIN->SetBinContent(etabin+1, phibin+1, ADC_val);
269  hADC_HCALIN->Fill(ADC_val);
270 
271  max_ADC = max(max_ADC, ADC_val);
272 
273  Float_t eta = eta_map[etabin];
274  Float_t phi = phi_map[phibin];
275 
276  if(first_entry){
277  output << "\"HCALIN\": [{ \"eta\": " << eta << ", \"phi\": " << phi << " , \"e\": " << e << "}\n";
278  first_entry = false;
279  }
280  else output << ",{ \"eta\": " << eta << ", \"phi\": " << phi << " , \"e\": " << e << "}\n";
281  // if(j == 10) break;
282  }
283  cout << "HCALIN Max ADC: " << max_ADC << endl;
284  output << "],\n";
285  }
286  else output << "\"HCALIN\": [{}],\n";
287 
288  first_entry = true;
289  max_ADC = ADC_max;
290  if(finput_HCALOUT && tree_HCALOUT->GetEntries()) {
291  ADC = 0;
292  chan = 0;
293 
294  tree_HCALOUT->SetBranchAddress("adc",&ADC);
295  tree_HCALOUT->SetBranchAddress("chan",&chan);
296  tree_HCALOUT->GetEntry(0);
297 
298  UInt_t nchannels = ADC->size();
299 
300  for(UInt_t j = 0; j < nchannels; ++j) {
301  Float_t ADC_val = ADC->at(j);
302  // Float_t e = 0.002*ADC_val;
303  Float_t e = ADC_val;
304  // if(ADC_val < 100) continue;
305  Int_t channel = chan->at(j);
306  if(channel >= channels_max) {
307  cout << "HCALOUT invalid channel number: " << channel << endl;
308  continue;
309  }
310 
311  UInt_t key = TowerInfoDefs::encode_hcal(channel);
312  UInt_t etabin = TowerInfoDefs::getCaloTowerEtaBin(key);
313  UInt_t phibin = TowerInfoDefs::getCaloTowerPhiBin(key);
314 
315  h2ADC_HCALOUT->SetBinContent(etabin+1, phibin+1, ADC_val);
316  hADC_HCALOUT->Fill(ADC_val);
317 
318  max_ADC = max(max_ADC, ADC_val);
319 
320  Float_t eta = eta_map[etabin];
321  Float_t phi = phi_map[phibin];
322 
323  if(first_entry) {
324  output << "\"HCALOUT\": [{ \"eta\": " << eta << ", \"phi\": " << phi << " , \"e\": " << e << "}\n";
325  first_entry = false;
326  }
327  else output << ",{ \"eta\": " << eta << ", \"phi\": " << phi << " , \"e\": " << e << "}\n";
328  // if(j == 10) break;
329  // cout << "etabin: " << etabin << ", phibin: " << phibin << ", ADC: " << ADC_val << endl;
330  // cout << ",{ \"eta\": " << eta << ", \"phi\": " << phi << " , \"e\": " << ADC_val << "}" << endl;
331  }
332  cout << "HCALOUT Max ADC: " << max_ADC << endl;
333  output << "]\n";
334  }
335  else output << "\"HCALOUT\": [{}]\n";
336 }
337 
338 void myAnalysis::finalize(const string& run, const string& event) {
339 
340  output << "},\n \
341  \"TRACKS\": {\n \
342  \"INNERTRACKER\": []\n \
343  }}";
344 
345 
346  auto c1 = new TCanvas();
347 
348  c1->cd();
349  // c1->SetTickx();
350  // c1->SetTicky();
351 
352  c1->SetCanvasSize(1000, 1500);
353 
354  c1->SetLeftMargin(.12);
355  c1->SetRightMargin(.15);
356  // gPad->SetLogz();
357 
358  string outputFile = "event-display/plots-" + run + "-" + event + ".pdf";
359  c1->Print((outputFile + "[").c_str(), "pdf portrait");
360 
361  gPad->SetLogz();
362  h2ADC_CEMC->SetStats(0);
363  h2ADC_CEMC->SetTitle(("ADC CEMC, Run: " + run + ", event: " + event).c_str());
364  // h2ADC_CEMC->SetTickLength(0,"xy");
365  // h2ADC_CEMC->SetTitleOffset(0.5, "z");
366  // h2ADC_CEMC->SetNdivisions(12,"x");
367  // h2ADC_CEMC->SetNdivisions(32,"y");
368  h2ADC_CEMC->Draw("COLZ1");
369 
370  auto tline = new TLine();
371 
372  for (UInt_t i = 0; i < 32; ++i) {
373  tline->DrawLine(eta_min, 8*i-0.5, eta_max, 8*i-0.5);
374  }
375  for (UInt_t i = 0; i < 12; ++i) {
376  tline->DrawLine(8*i-0.5, phi_min, 8*i-0.5, phi_max);
377  }
378  c1->Print((outputFile).c_str(), "pdf portrait");
379 
380  // hADC_CEMC->SetStats(0);
381  c1->SetCanvasSize(1500, 1000);
382  gPad->SetLogy();
383  hADC_CEMC->SetTitle(("ADC CEMC, Run: " + run + ", event: " + event).c_str());
384  hADC_CEMC->Draw();
385  c1->Print((outputFile).c_str(), "pdf portrait");
386  gPad->SetLogy(0);
387 
388  // gPad->SetLogz(0);
389  // h2ADC_CEMC->SetMinimum(4000);
390  // c1->Print((outputFile).c_str(), "pdf portrait");
391 
392  if(finput_HCALIN) {
393  c1->SetCanvasSize(1000, 1500);
394  gPad->SetLogz();
395  h2ADC_HCALIN->SetStats(0);
396  h2ADC_HCALIN->SetTitle(("ADC HCALIN, Run: " + run + ", event: " + event).c_str());
397  h2ADC_HCALIN->Draw("COLZ1");
398  c1->Print((outputFile).c_str(), "pdf portrait");
399 
400  c1->SetCanvasSize(1500, 1000);
401  gPad->SetLogy();
402  hADC_HCALIN->SetTitle(("ADC HCALIN, Run: " + run + ", event: " + event).c_str());
403  hADC_HCALIN->Draw();
404  c1->Print((outputFile).c_str(), "pdf portrait");
405  gPad->SetLogy(0);
406 
407  // gPad->SetLogz(0);
408  // h2ADC_HCALIN->SetMinimum(100);
409  // c1->Print((outputFile).c_str(), "pdf portrait");
410  }
411 
412 
413  if(finput_HCALOUT) {
414  c1->SetCanvasSize(1000, 1500);
415  gPad->SetLogz();
416  h2ADC_HCALOUT->SetStats(0);
417  h2ADC_HCALOUT->SetTitle(("ADC HCALOUT, Run: " + run + ", event: " + event).c_str());
418  h2ADC_HCALOUT->Draw("COLZ1");
419  c1->Print((outputFile).c_str(), "pdf portrait");
420 
421  c1->SetCanvasSize(1500, 1000);
422  gPad->SetLogy();
423  hADC_HCALOUT->SetTitle(("ADC HCALOUT, Run: " + run + ", event: " + event).c_str());
424  hADC_HCALOUT->Draw();
425  c1->Print((outputFile).c_str(), "pdf portrait");
426  gPad->SetLogy(0);
427 
428  // gPad->SetLogz(0);
429  // h2ADC_HCALOUT->SetMinimum(100);
430  // c1->Print((outputFile).c_str(), "pdf portrait");
431  }
432 
433  c1->Print((outputFile + "]").c_str(), "pdf portrait");
434 
435  // Close files
436  if(finput_CEMC) finput_CEMC->Close();
437  if(finput_HCALIN) finput_HCALIN->Close();
438  if(finput_HCALOUT) finput_HCALOUT->Close();
439  output.close();
440 }
441 
442 void test_hcal() {
443  UInt_t max_eta = 0;
444  UInt_t max_phi = 0;
445  for (UInt_t channel = 0; channel < myAnalysis::channels_hcal_bins; ++channel) {
446  UInt_t key = TowerInfoDefs::encode_hcal(channel);
447  UInt_t etabin = TowerInfoDefs::getCaloTowerEtaBin(key);
448  UInt_t phibin = TowerInfoDefs::getCaloTowerPhiBin(key);
449 
450  cout << "channel: " << channel << ", etabin: " << etabin << ", phibin: " << phibin << endl;
451 
452  max_eta = max(max_eta, etabin);
453  max_phi = max(max_phi, phibin);
454  }
455  cout << "max etabin: " << max_eta << ", max phibin: " << max_phi << endl;
456 }
457 
458 void event_display(const string& run,
459  const string& event,
460  const string& input_CEMC,
461  const string& outputFile="event-display/test.json",
462  const string& input_HCALIN = "",
463  const string& input_HCALOUT = "") {
464  myAnalysis::init(input_CEMC, outputFile, input_HCALIN, input_HCALOUT);
465  myAnalysis::analyze();
466  myAnalysis::finalize(run, event);
467 }
468 
469 # ifndef __CINT__
470 int main(int argc, char* argv[]) {
471  if(argc < 1 || argc > 7){
472  cout << "usage: ./bin/event-display run event input_CEMC outputFile input_HCALIN input_HCALOUT" << endl;
473  cout << "run: run number." << endl;
474  cout << "event: event number." << endl;
475  cout << "input_CEMC: Location of CEMC root file." << endl;
476  cout << "outputFile: output json file. Default = event-display/test.json." << endl;
477  cout << "input_HCALIN: Location of HCALIN root file." << endl;
478  cout << "input_HCALOUT: Location of HCALOUT root file." << endl;
479  return 1;
480  }
481 
482  string run;
483  string event;
484  string input_CEMC;
485  string output = "event-display/test.json";
486  string input_HCALIN = "";
487  string input_HCALOUT = "";
488 
489  if(argc >= 2) {
490  run = argv[1];
491  }
492  if(argc >= 3) {
493  event = argv[2];
494  }
495  if(argc >= 4) {
496  input_CEMC = argv[3];
497  }
498  if(argc >= 5) {
499  output = argv[4];
500  }
501  if(argc >= 6) {
502  input_HCALIN = argv[5];
503  }
504  if(argc >= 7) {
505  input_HCALOUT = argv[6];
506  }
507  event_display(run, event, input_CEMC, output, input_HCALIN, input_HCALOUT);
508 
509  cout << "done" << endl;
510  return 0;
511 }
512 # endif