Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
CaloValid.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file CaloValid.cc
1 #include "CaloValid.h"
2 
3 
4 // Calo includes
5 #include <calobase/RawCluster.h>
6 #include <calobase/RawClusterContainer.h>
7 #include <calobase/RawClusterUtility.h>
8 #include <calobase/RawTowerGeomContainer.h>
9 #include <calobase/TowerInfo.h>
10 #include <calobase/TowerInfoContainer.h>
11 
12 #include <fun4all/Fun4AllHistoManager.h>
13 #include <fun4all/Fun4AllReturnCodes.h>
14 
15 #include <phool/getClass.h>
16 #include <phool/phool.h> // for PHWHERE
17 
18 #include <globalvertex/GlobalVertex.h>
19 #include <globalvertex/GlobalVertexMap.h>
20 
21 // MBD
22 #include <mbd/MbdPmtContainer.h>
23 #include <mbd/MbdPmtHit.h>
24 
25 #include <TFile.h>
26 #include <TH1.h>
27 #include <TH2.h>
28 #include <TLorentzVector.h>
29 #include <TNtuple.h>
30 #include <TProfile2D.h>
31 #include <TSystem.h>
32 #include <TTree.h>
33 
34 #include <cmath> // for log10, pow, sqrt, abs, M_PI
35 #include <iostream> // for operator<<, endl, basic_...
36 #include <limits>
37 #include <map> // for operator!=, _Rb_tree_con...
38 #include <string>
39 #include <utility> // for pair
40 
41 CaloValid::CaloValid(const std::string& name, const std::string& filename)
42  : SubsysReco(name)
43  , detector("HCALIN")
44  , outfilename(filename)
45 {
46 }
47 
48 CaloValid::~CaloValid()
49 {
50  delete hm;
51  delete towerntuple;
52  delete clusterntuple;
53 }
54 
55 int CaloValid::Init(PHCompositeNode* /*unused*/)
56 {
57  delete hm; // this is a null pointer - make cppcheck happy
58  hm = new Fun4AllHistoManager(Name());
59  // create and register your histos (all types) here
60 
61  if (m_debug)
62  {
63  std::cout << "In CaloValid::Init" << std::endl;
64  }
65 
66 
67  outfile = new TFile(outfilename.c_str(), "RECREATE");
68 
69  h_emcal_mbd_correlation = new TH2F("h_emcal_mbd_correlation", ";emcal;mbd", 100, 0, 1, 100, 0, 1);
70  h_ohcal_mbd_correlation = new TH2F("h_ohcal_mbd_correlation", ";ohcal;mbd", 100, 0, 1, 100, 0, 1);
71  h_ihcal_mbd_correlation = new TH2F("h_ihcal_mbd_correlation", ";ihcal;mbd", 100, 0, 1, 100, 0, 1);
72  h_emcal_hcal_correlation = new TH2F("h_emcal_hcal_correlation", ";emcal;hcal", 100, 0, 1, 100, 0, 1);
73  h_cemc_etaphi = new TH2F("h_cemc_etaphi", ";eta;phi", 96, 0, 96, 256, 0, 256);
74  h_hcalin_etaphi = new TH2F("h_ihcal_etaphi", ";eta;phi", 24, 0, 24, 64, 0, 64);
75  h_hcalout_etaphi = new TH2F("h_ohcal_etaphi", ";eta;phi", 24, 0, 24, 64, 0, 64);
76 
77  h_cemc_etaphi_wQA = new TH2F("h_cemc_etaphi_wQA", ";eta;phi", 96, 0, 96, 256, 0, 256);
78  h_hcalin_etaphi_wQA = new TH2F("h_ihcal_etaphi_wQA", ";eta;phi", 24, 0, 24, 64, 0, 64);
79  h_hcalout_etaphi_wQA = new TH2F("h_ohcal_etaphi_wQA", ";eta;phi", 24, 0, 24, 64, 0, 64);
80  h_ihcal_status = new TH1F("h_ihcal_status","",256,0,256);
81  h_ohcal_status = new TH1F("h_ohcal_status","",256,0,256);
82  h_cemc_status = new TH1F("h_cemc_status","",256,0,256);
83 
84  h_cemc_e_chi2 = LogYHist2D("h_cemc_e_chi2", "", 500, -2, 30, 1000, 0.5, 5e6);
85  h_ihcal_e_chi2 = LogYHist2D("h_ihcal_e_chi2", "", 500, -2, 30, 1000, 0.5, 5e6);
86  h_ohcal_e_chi2 = LogYHist2D("h_ohcal_e_chi2", "", 500, -2, 30, 1000, 0.5, 5e6);
87 
88  h_cemc_etaphi_time = new TProfile2D("h_cemc_etaphi_time", ";eta;phi", 96, 0, 96, 256, 0, 256, -10, 10);
89  h_hcalin_etaphi_time = new TProfile2D("h_ihcal_etaphi_time", ";eta;phi", 24, 0, 24, 64, 0, 64, -10, 10);
90  h_hcalout_etaphi_time = new TProfile2D("h_ohcal_etaphi_time", ";eta;phi", 24, 0, 24, 64, 0, 64, -10, 10);
91 
92  h_cemc_etaphi_fracHitADC = new TProfile2D("h_cemc_etaphi_fracHitADC", ";eta;phi", 96, 0, 96, 256, 0, 256, -10, 10);
93  h_hcalin_etaphi_fracHitADC = new TProfile2D("h_ihcal_etaphi_fracHitADC", ";eta;phi", 24, 0, 24, 64, 0, 64, -10, 10);
94  h_hcalout_etaphi_fracHitADC = new TProfile2D("h_ohcal_etaphi_fracHitADC", ";eta;phi", 24, 0, 24, 64, 0, 64, -10, 10);
95 
96  h_cemc_etaphi_badChi2 = new TProfile2D("h_cemc_etaphi_badChi2", ";eta;phi", 96, 0, 96, 256, 0, 256, -10, 10);
97  h_hcalin_etaphi_badChi2 = new TProfile2D("h_ihcal_etaphi_badChi2", ";eta;phi", 24, 0, 24, 64, 0, 64, -10, 10);
98  h_hcalout_etaphi_badChi2 = new TProfile2D("h_ohcal_etaphi_badChi2", ";eta;phi", 24, 0, 24, 64, 0, 64, -10, 10);
99 
100  // 1D distributions
101  h_InvMass = new TH1F("h_InvMass", "Invariant Mass", 120, 0, 1.2);
102 
103  // ZDC QA plots
104  hzdcSouthraw = new TH1D("hzdcSouthraw", "hzdcSouthraw", 1500, 0, 15000);
105  hzdcNorthraw = new TH1D("hzdcNorthraw", "hzdcNorthraw", 1500, 0, 15000);
106  hzdcSouthcalib = new TH1D("hzdcSouthcalib", "hzdcSouthcalib", 1500, 0, 15000);
107  hzdcNorthcalib = new TH1D("hzdcNorthcalib", "hzdcNorthcalib", 1500, 0, 15000);
108  h_totalzdc_e = new TH1D("h_totalzdc_e", "", 200, 0, 2e4);
109  h_emcal_zdc_correlation = new TH2F("h_zdc_emcal_correlation", ";emcal;zdc", 100, 0, 1, 100, 0, 1);
110 
111  // vertex distributions
112  hvtx_z_raw = new TH1D("hvtx_z_raw", "hvtx_z_raw", 201, -100.5, 100.5);
113  hvtx_z_cut = new TH1D("hvtx_z_cut", "hvtx_z_cut", 201, -100.5, 100.5);
114 
115  // raw timing information
116  hzdctime_cut = new TH1D("hzdctime_cut", "hzdctime_cut", 50, -17.5, 32.5);
117  hemcaltime_cut = new TH1D("hemcaltime_cut", "hemcaltime_cut", 50, -17.5, 32.5);
118  hihcaltime_cut = new TH1D("hihcaltime_cut", "hihcaltime_cut", 50, -17.5, 32.5);
119  hohcaltime_cut = new TH1D("hohcaltime_cut", "hohcaltime_cut", 50, -17.5, 32.5);
120 
121  // extracted timing information
122  hzdctime = new TH1D("hzdctime", "hzdctime", 50, -17.5, 32.5);
123  hemcaltime = new TH1D("hemcaltime", "hemcaltime", 50, -17.5, 32.5);
124  hihcaltime = new TH1D("hihcaltime", "hihcaltime", 50, -17.5, 32.5);
125  hohcaltime = new TH1D("hohcaltime", "hohcaltime", 50, -17.5, 32.5);
126 
127  // cluster QA
128  h_etaphi_clus = new TH2F("h_etaphi_clus", "", 140, -1.2, 1.2, 64, -1 * M_PI, M_PI);
129  h_clusE = new TH1F("h_clusE", "", 100, 0, 10);
130 
131  if (m_debug)
132  {
133  std::cout << "Leaving CaloValid::Init" << std::endl;
134  }
135  return 0;
136 }
137 
138 int CaloValid::process_event(PHCompositeNode* topNode)
139 {
140  _eventcounter++;
141 
142  process_towers(topNode);
143 
144  return Fun4AllReturnCodes::EVENT_OK;
145 }
146 
147 int CaloValid::process_towers(PHCompositeNode* topNode)
148 {
149  if (m_debug)
150  {
151  std::cout << _eventcounter << std::endl;
152  }
153 
154  float totalcemc = 0.;
155  float totalihcal = 0.;
156  float totalohcal = 0.;
157  float totalmbd = 0.;
158  float totalzdc = 0.;
159  float totalzdcsouthraw = 0.;
160  float totalzdcnorthraw = 0.;
161  float totalzdcsouthcalib = 0.;
162  float totalzdcnorthcalib = 0.;
163 
164  float emcaldownscale = 1000000. / 800.;
165  float ihcaldownscale = 40000. / 300.;
166  float ohcaldownscale = 250000. / 600.;
167  float mbddownscale = 2000.0;
168  float zdcdownscale = 1e4;
169 
170  float adc_threshold = 15.;
171 
172  float emcal_hit_threshold = 0.5; // GeV
173  float ohcal_hit_threshold = 0.5;
174  float ihcal_hit_threshold = 0.25;
175 
176  int max_zdc_t = -1;
177  int max_emcal_t = -1;
178  int max_ihcal_t = -1;
179  int max_ohcal_t = -1;
180 
181  // get time estimate
182  max_zdc_t = Getpeaktime(hzdctime_cut);
183  max_emcal_t = Getpeaktime(hemcaltime_cut);
184  max_ihcal_t = Getpeaktime(hihcaltime_cut);
185  max_ohcal_t = Getpeaktime(hohcaltime_cut);
186 
187  //----------------------------------vertex------------------------------------------------------//
188  GlobalVertexMap* vertexmap = findNode::getClass<GlobalVertexMap>(topNode, "GlobalVertexMap");
189  if (!vertexmap)
190  {
191  std::cout << "CaloValid GlobalVertexMap node is missing" << std::endl;
192  }
193  float vtx_z = std::numeric_limits<float>::quiet_NaN();
194  if (vertexmap && !vertexmap->empty())
195  {
196  GlobalVertex* vtx = vertexmap->begin()->second;
197  if (vtx)
198  {
199  vtx_z = vtx->get_z();
200  }
201  hvtx_z_raw->Fill(vtx_z);
202  }
203 
204  //---------------------------calibrated towers-------------------------------//
205  {
206  TowerInfoContainer* towers = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_CEMC");
207  if (towers)
208  {
209  int size = towers->size(); // online towers should be the same!
210  for (int channel = 0; channel < size; channel++)
211  {
212  TowerInfo* tower = towers->get_tower_at_channel(channel);
213  float offlineenergy = tower->get_energy();
214  unsigned int towerkey = towers->encode_key(channel);
215  int ieta = towers->getTowerEtaBin(towerkey);
216  int iphi = towers->getTowerPhiBin(towerkey);
217  int _time = tower->get_time();
218  h_cemc_e_chi2->Fill(offlineenergy, tower->get_chi2());
219  float _timef = tower->get_time_float();
220  hemcaltime_cut->Fill(_time);
221  bool isGood = tower->get_isGood();
222  uint8_t status = tower->get_status();
223  for (int is=0; is<8; is++)
224  {
225  if(status & 1U) // clang-tidy mark 1 as unsigned
226  {
227  h_cemc_status->Fill(is);
228  }
229  status = status >> 1U; // clang-tidy mark 1 as unsigned
230  }
231  if (_time > (max_emcal_t - _range) && _time < (max_emcal_t + _range))
232  {
233  totalcemc += offlineenergy;
234  hemcaltime->Fill(_time);
235  if (offlineenergy > emcal_hit_threshold)
236  {
237  h_cemc_etaphi_time->Fill(ieta, iphi, _timef);
238  h_cemc_etaphi->Fill(ieta, iphi);
239  if (isGood)
240  {
241  h_cemc_etaphi_wQA->Fill(ieta, iphi, offlineenergy);
242  }
243  if (tower->get_isBadChi2())
244  {
245  h_cemc_etaphi_badChi2->Fill(ieta, iphi, 1);
246  }
247  else
248  {
249  h_cemc_etaphi_badChi2->Fill(ieta, iphi, 0);
250  }
251  }
252  }
253  }
254  }
255  }
256 
257  {
258  TowerInfoContainer* towers = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_HCALIN");
259  if (towers)
260  {
261  int size = towers->size(); // online towers should be the same!
262  for (int channel = 0; channel < size; channel++)
263  {
264  TowerInfo* tower = towers->get_tower_at_channel(channel);
265  float offlineenergy = tower->get_energy();
266  unsigned int towerkey = towers->encode_key(channel);
267  int ieta = towers->getTowerEtaBin(towerkey);
268  int iphi = towers->getTowerPhiBin(towerkey);
269  int _time = tower->get_time();
270  float _timef = tower->get_time_float();
271  hihcaltime_cut->Fill(_time);
272  h_ihcal_e_chi2->Fill(offlineenergy, tower->get_chi2());
273  bool isGood = tower->get_isGood();
274  h_ihcal_status->Fill(tower->get_status());
275 
276  uint8_t status = tower->get_status();
277  for (int is=0; is<8; is++)
278  {
279  if(status & 1U) // clang-tidy mark 1 as unsigned
280  {
281  h_ihcal_status->Fill(is);
282  }
283  status = status >> 1U; // clang-tidy mark 1 as unsigned
284  }
285 
286  if (_time > (max_ihcal_t - _range) && _time < (max_ihcal_t + _range))
287  {
288  totalihcal += offlineenergy;
289  hihcaltime->Fill(_time);
290 
291  if (offlineenergy > ihcal_hit_threshold)
292  {
293  h_hcalin_etaphi->Fill(ieta, iphi);
294  h_hcalin_etaphi_time->Fill(ieta, iphi, _timef);
295  if (isGood)
296  {
297  h_hcalin_etaphi_wQA->Fill(ieta, iphi, offlineenergy);
298  }
299  if (tower->get_isBadChi2())
300  {
301  h_hcalin_etaphi_badChi2->Fill(ieta, iphi, 1);
302  }
303  else
304  {
305  h_hcalin_etaphi_badChi2->Fill(ieta, iphi, 0);
306  }
307  }
308  }
309  }
310  }
311  }
312 
313  {
314  TowerInfoContainer* towers = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_HCALOUT");
315  if (towers)
316  {
317  int size = towers->size(); // online towers should be the same!
318  for (int channel = 0; channel < size; channel++)
319  {
320  TowerInfo* tower = towers->get_tower_at_channel(channel);
321  float offlineenergy = tower->get_energy();
322  unsigned int towerkey = towers->encode_key(channel);
323  int ieta = towers->getTowerEtaBin(towerkey);
324  int iphi = towers->getTowerPhiBin(towerkey);
325  int _time = tower->get_time();
326  float _timef = tower->get_time_float();
327  hihcaltime_cut->Fill(_time);
328  hohcaltime_cut->Fill(_time);
329  h_ohcal_e_chi2->Fill(offlineenergy, tower->get_chi2());
330  bool isGood = tower->get_isGood();
331  h_ohcal_status->Fill(tower->get_status());
332 
333  uint8_t status = tower->get_status();
334  for (int is=0; is<8; is++)
335  {
336  if(status & 1U) // clang-tidy mark 1 as unsigned
337  {
338  h_ohcal_status->Fill(is);
339  }
340  status = status >> 1U; // clang-tidy mark 1 as unsigned
341  }
342 
343  if (_time > (max_ohcal_t - _range) && _time < (max_ohcal_t + _range))
344  {
345  totalohcal += offlineenergy;
346  hohcaltime->Fill(_time);
347 
348  if (offlineenergy > ohcal_hit_threshold)
349  {
350  h_hcalout_etaphi->Fill(ieta, iphi);
351  h_hcalout_etaphi_time->Fill(ieta, iphi, _timef);
352  if (isGood)
353  {
354  h_hcalout_etaphi_wQA->Fill(ieta, iphi, offlineenergy);
355  }
356  if (tower->get_isBadChi2())
357  {
358  h_hcalout_etaphi_badChi2->Fill(ieta, iphi, 1);
359  }
360  else
361  {
362  h_hcalout_etaphi_badChi2->Fill(ieta, iphi, 0);
363  }
364  }
365  }
366  }
367  }
368  }
369 
370  {
371  TowerInfoContainer* towers = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_ZDC");
372  if (towers)
373  {
374  int size = towers->size(); // online towers should be the same!
375  for (int channel = 0; channel < size; channel++)
376  {
377  TowerInfo* tower = towers->get_tower_at_channel(channel);
378  float offlineenergy = tower->get_energy();
379  int _time = towers->get_tower_at_channel(channel)->get_time();
380  hzdctime_cut->Fill(_time);
381  if (channel == 0 || channel == 2 || channel == 4)
382  {
383  totalzdcsouthcalib += offlineenergy;
384  }
385  if (channel == 8 || channel == 10 || channel == 12)
386  {
387  totalzdcnorthcalib += offlineenergy;
388  }
389  if (channel == 0 || channel == 2 || channel == 4 || channel == 8 || channel == 10 || channel == 12)
390  {
391  if (_time > (max_zdc_t - _range) && _time < (max_zdc_t + _range))
392  {
393  totalzdc += offlineenergy;
394  hzdctime->Fill(_time);
395  }
396  }
397  }
398  }
399  }
400 
401  //-------------------------- raw tower ------------------------------//
402  {
403  TowerInfoContainer* towers = findNode::getClass<TowerInfoContainer>(topNode, "TOWERS_ZDC");
404  if (towers)
405  {
406  int size = towers->size(); // online towers should be the same!
407  for (int channel = 0; channel < size; channel++)
408  {
409  TowerInfo* tower = towers->get_tower_at_channel(channel);
410  float offlineenergy = tower->get_energy();
411  if (channel == 0 || channel == 2 || channel == 4)
412  {
413  totalzdcsouthraw += offlineenergy;
414  }
415  if (channel == 8 || channel == 10 || channel == 12)
416  {
417  totalzdcnorthraw += offlineenergy;
418  }
419  }
420  }
421  }
422 
423  {
424  TowerInfoContainer* towers = findNode::getClass<TowerInfoContainer>(topNode, "TOWERS_CEMC");
425  if (towers)
426  {
427  int size = towers->size(); // online towers should be the same!
428  for (int channel = 0; channel < size; channel++)
429  {
430  TowerInfo* tower = towers->get_tower_at_channel(channel);
431  unsigned int towerkey = towers->encode_key(channel);
432  int ieta = towers->getTowerEtaBin(towerkey);
433  int iphi = towers->getTowerPhiBin(towerkey);
434 
435  float raw_energy = tower->get_energy();
436  if (raw_energy > adc_threshold)
437  {
438  h_cemc_etaphi_fracHitADC->Fill(ieta, iphi, 1);
439  }
440  else
441  {
442  h_cemc_etaphi_fracHitADC->Fill(ieta, iphi, 0);
443  }
444  }
445  }
446  }
447  {
448  TowerInfoContainer* towers = findNode::getClass<TowerInfoContainer>(topNode, "TOWERS_HCALOUT");
449  if (towers)
450  {
451  int size = towers->size(); // online towers should be the same!
452  for (int channel = 0; channel < size; channel++)
453  {
454  TowerInfo* tower = towers->get_tower_at_channel(channel);
455  unsigned int towerkey = towers->encode_key(channel);
456  int ieta = towers->getTowerEtaBin(towerkey);
457  int iphi = towers->getTowerPhiBin(towerkey);
458 
459  float raw_energy = tower->get_energy();
460  if (raw_energy > adc_threshold)
461  {
462  h_hcalout_etaphi_fracHitADC->Fill(ieta, iphi, 1);
463  }
464  else
465  {
466  h_hcalout_etaphi_fracHitADC->Fill(ieta, iphi, 0);
467  }
468  }
469  }
470  }
471  {
472  TowerInfoContainer* towers = findNode::getClass<TowerInfoContainer>(topNode, "TOWERS_HCALIN");
473  if (towers)
474  {
475  int size = towers->size(); // online towers should be the same!
476  for (int channel = 0; channel < size; channel++)
477  {
478  TowerInfo* tower = towers->get_tower_at_channel(channel);
479  unsigned int towerkey = towers->encode_key(channel);
480  int ieta = towers->getTowerEtaBin(towerkey);
481  int iphi = towers->getTowerPhiBin(towerkey);
482 
483  float raw_energy = tower->get_energy();
484  if (raw_energy > adc_threshold)
485  {
486  h_hcalin_etaphi_fracHitADC->Fill(ieta, iphi, 1);
487  }
488  else
489  {
490  h_hcalin_etaphi_fracHitADC->Fill(ieta, iphi, 0);
491  }
492  }
493  }
494  }
495 
496  //--------------------------- MBD ----------------------------------------//
497  MbdPmtContainer* bbcpmts = findNode::getClass<MbdPmtContainer>(topNode, "MbdPmtContainer");
498  if (!bbcpmts)
499  {
500  std::cout << "makeMBDTrees::process_event: Could not find MbdPmtContainer, aborting" << std::endl;
501  return Fun4AllReturnCodes::ABORTEVENT;
502  }
503 
504  int nPMTs = bbcpmts->get_npmt();
505  for (int i = 0; i < nPMTs; i++)
506  {
507  MbdPmtHit* mbdpmt = bbcpmts->get_pmt(i);
508  float pmtadc = mbdpmt->get_q();
509  totalmbd += pmtadc;
510  }
511 
512  h_emcal_mbd_correlation->Fill(totalcemc / emcaldownscale, totalmbd / mbddownscale);
513  h_ihcal_mbd_correlation->Fill(totalihcal / ihcaldownscale, totalmbd / mbddownscale);
514  h_ohcal_mbd_correlation->Fill(totalohcal / ohcaldownscale, totalmbd / mbddownscale);
515  h_emcal_hcal_correlation->Fill(totalcemc / emcaldownscale, totalohcal / ohcaldownscale);
516  h_emcal_zdc_correlation->Fill(totalcemc / emcaldownscale, totalzdc / zdcdownscale);
517  h_totalzdc_e->Fill(totalzdc);
518 
519  hzdcSouthraw->Fill(totalzdcsouthraw);
520  hzdcNorthraw->Fill(totalzdcnorthraw);
521  hzdcSouthcalib->Fill(totalzdcsouthcalib);
522  hzdcNorthcalib->Fill(totalzdcnorthcalib);
523 
524  //------------------------------ clusters & pi0 ------------------------------//
525  RawClusterContainer* clusterContainer = findNode::getClass<RawClusterContainer>(topNode, "CLUSTERINFO_POS_COR_CEMC");
526  if (!clusterContainer)
527  {
528  std::cout << PHWHERE << "funkyCaloStuff::process_event - Fatal Error - CLUSTER_CEMC node is missing. " << std::endl;
529  return 0;
530  }
531 
533  // geometry for hot tower/cluster masking
534  std::string towergeomnodename = "TOWERGEOM_CEMC";
535  RawTowerGeomContainer* m_geometry = findNode::getClass<RawTowerGeomContainer>(topNode, towergeomnodename);
536  if (!m_geometry)
537  {
538  std::cout << Name() << "::"
539  << "CreateNodeTree"
540  << ": Could not find node " << towergeomnodename << std::endl;
541  gSystem->Exit(1);
542  }
543 
544  // cuts
545  float emcMinClusE1 = 1.3; // 0.5;
546  float emcMinClusE2 = 0.7; // 0.5;
547  float emcMaxClusE = 100;
548  float maxAlpha = 0.6;
549 
550  if (totalcemc < 0.2 * emcaldownscale)
551  {
552  RawClusterContainer::ConstRange clusterEnd = clusterContainer->getClusters();
553  RawClusterContainer::ConstIterator clusterIter;
554  RawClusterContainer::ConstIterator clusterIter2;
555 
556  for (clusterIter = clusterEnd.first; clusterIter != clusterEnd.second; clusterIter++)
557  {
558  RawCluster* recoCluster = clusterIter->second;
559 
560  CLHEP::Hep3Vector vertex(0, 0, 0);
561  CLHEP::Hep3Vector E_vec_cluster = RawClusterUtility::GetECoreVec(*recoCluster, vertex);
562 
563  float clusE = E_vec_cluster.mag();
564  float clus_eta = E_vec_cluster.pseudoRapidity();
565  float clus_phi = E_vec_cluster.phi();
566  float clus_pt = E_vec_cluster.perp();
567  float clus_chisq = recoCluster->get_chi2();
568 
569  h_clusE->Fill(clusE);
570 
571  if (clusE < emcMinClusE1 || clusE > emcMaxClusE)
572  {
573  continue;
574  }
575  if (clus_chisq > 4)
576  {
577  continue;
578  }
579 
580  h_etaphi_clus->Fill(clus_eta, clus_phi);
581 
582  TLorentzVector photon1;
583  photon1.SetPtEtaPhiE(clus_pt, clus_eta, clus_phi, clusE);
584 
585  for (clusterIter2 = clusterEnd.first; clusterIter2 != clusterEnd.second; clusterIter2++)
586  {
587  if (clusterIter == clusterIter2)
588  {
589  continue;
590  }
591  RawCluster* recoCluster2 = clusterIter2->second;
592 
593  CLHEP::Hep3Vector vertex2(0, 0, 0);
594  CLHEP::Hep3Vector E_vec_cluster2 = RawClusterUtility::GetECoreVec(*recoCluster2, vertex2);
595 
596  float clus2E = E_vec_cluster2.mag();
597  float clus2_eta = E_vec_cluster2.pseudoRapidity();
598  float clus2_phi = E_vec_cluster2.phi();
599  float clus2_pt = E_vec_cluster2.perp();
600  float clus2_chisq = recoCluster2->get_chi2();
601 
602  if (clus2E < emcMinClusE2 || clus2E > emcMaxClusE)
603  {
604  continue;
605  }
606  if (clus2_chisq > 4)
607  {
608  continue;
609  }
610 
611  TLorentzVector photon2;
612  photon2.SetPtEtaPhiE(clus2_pt, clus2_eta, clus2_phi, clus2E);
613 
614  if (sqrt(pow(clusE - clus2E, 2)) / (clusE + clus2E) > maxAlpha)
615  {
616  continue;
617  }
618 
619  TLorentzVector pi0 = photon1 + photon2;
620  h_InvMass->Fill(pi0.M());
621  }
622  }
623  }
624  return Fun4AllReturnCodes::EVENT_OK;
625 }
626 
627 int CaloValid::End(PHCompositeNode* /*topNode*/)
628 {
629  outfile->cd();
630 
631  outfile->Write();
632  outfile->Close();
633  delete outfile;
634  hm->dumpHistos(outfilename, "UPDATE");
635  return 0;
636 }
637 
638 int CaloValid::Getpeaktime(TH1* h)
639 {
640  int getmaxtime, tcut = -1;
641 
642  for (int bin = 1; bin < h->GetNbinsX() + 1; bin++)
643  {
644  double c = h->GetBinContent(bin);
645  double max = h->GetMaximum();
646  int bincenter = h->GetBinCenter(bin);
647  if (max == c)
648  {
649  getmaxtime = bincenter;
650  if (getmaxtime != -1)
651  {
652  tcut = getmaxtime;
653  }
654  }
655  }
656 
657  return tcut;
658 }
659 
660 TH2F* CaloValid::LogYHist2D(const std::string& name, const std::string& title, int xbins_in, double xmin, double xmax, int ybins_in, double ymin, double ymax)
661 {
662  Double_t logymin = std::log10(ymin);
663  Double_t logymax = std::log10(ymax);
664  Double_t binwidth = (logymax - logymin) / ybins_in;
665  Double_t ybins[ybins_in + 1];
666 
667  for (Int_t i = 0; i <= ybins_in + 1; i++)
668  {
669  ybins[i] = pow(10, logymin + i * binwidth);
670  }
671 
672  TH2F* h = new TH2F(name.c_str(), title.c_str(), xbins_in, xmin, xmax, ybins_in, ybins);
673 
674  return h;
675 }