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