ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
RawClusterBuilderTemplate.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file RawClusterBuilderTemplate.cc
2 
3 #include "BEmcCluster.h"
4 #include "BEmcRec.h"
5 #include "BEmcRecCEMC.h"
6 #include "BEmcRecEEMC.h"
7 #include "BEmcRecFEMC.h"
8 
11 
12 #include <calobase/RawCluster.h>
13 #include <calobase/RawClusterContainer.h>
14 #include <calobase/RawClusterv1.h>
15 #include <calobase/RawTower.h>
16 #include <calobase/RawTowerContainer.h>
17 #include <calobase/RawTowerDefs.h>
18 #include <calobase/RawTowerGeom.h>
19 #include <calobase/RawTowerGeomContainer.h>
20 
22 #include <fun4all/SubsysReco.h>
23 
24 #include <phool/PHCompositeNode.h>
25 #include <phool/PHIODataNode.h>
26 #include <phool/PHNode.h>
27 #include <phool/PHNodeIterator.h>
28 #include <phool/PHObject.h>
29 #include <phool/getClass.h>
30 #include <phool/phool.h>
31 
32 #include <cmath>
33 #include <exception>
34 #include <fstream>
35 #include <iostream>
36 #include <map>
37 #include <stdexcept>
38 #include <utility>
39 #include <vector>
40 
41 using namespace std;
42 
44  : SubsysReco(name)
45  , _clusters(nullptr)
46  , bemc(nullptr)
47  , fEnergyNorm(1.)
48  , _min_tower_e(0.020)
49  , chkenergyconservation(0)
50  , detector("NONE")
51  , BINX0(0)
52  , NBINX(0)
53  , BINY0(0)
54  , NBINY(0)
55  , bPrintGeom(false)
56  , bProfProb(false)
57  , fProbNoiseParam(0.04)
58 {
59 }
60 
62 {
63  // one can delete null pointers
64  delete bemc;
65 }
66 
67 void RawClusterBuilderTemplate::Detector(const std::string &d)
68 {
69  detector = d;
70 
71  // Create proper BEmcRec object
72 
73  if (detector == "CEMC")
74  {
75  bemc = new BEmcRecCEMC();
76  }
77  else if (detector == "FEMC")
78  {
79  bemc = new BEmcRecFEMC();
80  }
81  else if (detector == "EEMC")
82  {
83  bemc = new BEmcRecEEMC();
84  }
85  else if (detector == "EEMC_crystal")
86  {
87  bemc = new BEmcRecEEMC();
88  }
89  else if (detector == "EEMC_glass")
90  {
91  bemc = new BEmcRecEEMC();
92  }
93  else
94  {
95  cout << "Warning from RawClusterBuilderTemplate::Detector(): no detector specific class "
96  << Name() << " defined for detector " << detector
97  << ". Default BEmcRec will be used" << endl;
98  bemc = new BEmcRec();
99  }
100 
101  // Set vertex
102  float vertex[3] = {0, 0, 0};
103  bemc->SetVertex(vertex);
104  // Set threshold
107 }
108 
110 {
111  // _emcprof = new BEmcProfile(fname);
112  bemc->LoadProfile(fname);
113 }
114 
116 {
117  if (bemc == nullptr)
118  {
119  cout << "Error in RawClusterBuilderTemplate::SetCylindricalGeometry()(): detector is not defined; use RawClusterBuilderTemplate::Detector() to define it" << endl;
120  return;
121  }
122 
124 }
125 
127 {
128  if (bemc == nullptr)
129  {
130  cout << "Error in RawClusterBuilderTemplate::SetPlanarGeometry()(): detector is not defined; use RawClusterBuilderTemplate::Detector() to define it" << endl;
131  return;
132  }
133 
135 }
136 
138 {
139  if (bemc == nullptr)
140  {
141  cout << "Error in RawClusterBuilderTemplate::InitRun(): detector is not defined; use RawClusterBuilderTemplate::Detector() to define it" << endl;
143  }
144 
145  try
146  {
147  CreateNodes(topNode);
148  }
149  catch (std::exception &e)
150  {
151  std::cout << PHWHERE << ": " << e.what() << std::endl;
152  throw;
153  }
154 
155  string towergeomnodename = "TOWERGEOM_" + detector;
156  RawTowerGeomContainer *towergeom = findNode::getClass<RawTowerGeomContainer>(topNode, towergeomnodename);
157  if (!towergeom)
158  {
159  cout << PHWHERE << ": Could not find node " << towergeomnodename << endl;
161  }
162 
163  int Calo_ID = towergeom->get_calorimeter_id();
164  // cout << endl << endl << endl << "Calorimeter ID: " << Calo_ID << endl << endl << endl;
165 
166  int ngeom = 0;
167  int ixmin = 999999;
168  int ixmax = -999999;
169  int iymin = 999999;
170  int iymax = -999999;
171  RawTowerGeomContainer::ConstRange begin_end_geom = towergeom->get_tower_geometries();
172  RawTowerGeomContainer::ConstIterator itr_geom = begin_end_geom.first;
173  for (; itr_geom != begin_end_geom.second; ++itr_geom)
174  {
175  RawTowerGeom *towerg = itr_geom->second;
176  RawTowerDefs::keytype towerid = towerg->get_id();
177  int ix = RawTowerDefs::decode_index2(towerid); // index2 is phi in CYL
178  int iy = RawTowerDefs::decode_index1(towerid); // index1 is eta in CYL
179  if (ixmin > ix) ixmin = ix;
180  if (ixmax < ix) ixmax = ix;
181  if (iymin > iy) iymin = iy;
182  if (iymax < iy) iymax = iy;
183  ngeom++;
184  }
185  if (Verbosity() > 1)
186  {
187  cout << "Info from RawClusterBuilderTemplate::InitRun(): Init geometry for "
188  << detector << ": N of geom towers: " << ngeom << "; ix = "
189  << ixmin << "-" << ixmax << ", iy = "
190  << iymin << "-" << iymax << endl;
191  }
192  if (ixmax < ixmin || iymax < iymin)
193  {
194  cout << "Error in RawClusterBuilderTemplate::InitRun(): wrong geometry data for detector "
195  << detector << endl;
197  }
198 
199  BINX0 = ixmin;
200  NBINX = ixmax - ixmin + 1;
201  BINY0 = iymin;
202  NBINY = iymax - iymin + 1;
203 
204  bemc->SetDim(NBINX, NBINY);
205 
206  itr_geom = begin_end_geom.first;
207  for (; itr_geom != begin_end_geom.second; ++itr_geom)
208  {
209  RawTowerGeom *towerg = itr_geom->second;
210  RawTowerDefs::keytype towerid = towerg->get_id();
211  // int itype = towerg->get_tower_type();
212  // if( itype==2 ) { // PbSc
213  int ix = RawTowerDefs::decode_index2(towerid); // index2 is phi in CYL
214  int iy = RawTowerDefs::decode_index1(towerid); // index1 is eta in CYL
215  ix -= BINX0;
216  iy -= BINY0;
217 
218  bemc->SetTowerGeometry(ix, iy, towerg->get_center_x(), towerg->get_center_y(), towerg->get_center_z());
219  bemc->SetCalotype(Calo_ID);
220  if (Calo_ID == RawTowerDefs::EEMC ||
221  Calo_ID == RawTowerDefs::EEMC_crystal ||
222  Calo_ID == RawTowerDefs::EEMC_glass)
223  {
224  bemc->SetScinSize(towerg->get_size_z());
225  }
226  }
227 
229 
230  if (bPrintGeom)
231  {
232  string fname = "geom_" + detector + ".txt";
233  // bemc->PrintTowerGeometry("geom.txt");
234  bemc->PrintTowerGeometry(fname);
235  // PrintCylGeom(towergeom,"phieta.txt");
236  }
237 
239 }
240 
242 {
243  ofstream outfile(fname);
244  if (!outfile.is_open())
245  {
246  cout << "Error in BEmcRec::RawClusterBuilderTemplate::PrintCylGeom(): Failed to open file "
247  << fname << endl;
248  return;
249  }
250  outfile << NBINX << " " << NBINY << endl;
251  for (int ip = 0; ip < NBINX; ip++)
252  {
253  outfile << ip << " " << towergeom->get_phicenter(ip) << endl;
254  }
255  for (int ip = 0; ip < NBINY; ip++)
256  {
257  outfile << ip << " " << towergeom->get_etacenter(ip) << endl;
258  }
259  outfile.close();
260 }
261 
262 bool RawClusterBuilderTemplate::Cell2Abs(RawTowerGeomContainer */*towergeom*/, float /*phiC*/, float /*etaC*/, float &phi, float &eta)
263 {
264  phi = eta = 0;
265  return false;
266 }
267 
269 {
270  if (bemc == nullptr)
271  {
272  cout << "Error in RawClusterBuilderTemplate::process_event(): detector is not defined; use RawClusterBuilderTemplate::Detector() to define it" << endl;
274  }
275 
276  string towernodename = "TOWER_CALIB_" + detector;
277  // Grab the towers
278  RawTowerContainer *towers = findNode::getClass<RawTowerContainer>(topNode, towernodename);
279  if (!towers)
280  {
281  std::cout << PHWHERE << ": Could not find node " << towernodename << std::endl;
283  }
284  string towergeomnodename = "TOWERGEOM_" + detector;
285  RawTowerGeomContainer *towergeom = findNode::getClass<RawTowerGeomContainer>(topNode, towergeomnodename);
286  if (!towergeom)
287  {
288  cout << PHWHERE << ": Could not find node " << towergeomnodename << endl;
290  }
291 
292  // Get vertex
293  float vx = 0;
294  float vy = 0;
295  float vz = 0;
296  GlobalVertexMap *vertexmap = findNode::getClass<GlobalVertexMap>(topNode, "GlobalVertexMap");
297  if (vertexmap)
298  {
299  if (!vertexmap->empty())
300  {
301  GlobalVertex *vertex = (vertexmap->begin()->second);
302  vx = vertex->get_x();
303  vy = vertex->get_y();
304  vz = vertex->get_z();
305  }
306  }
307 
308  // Set vertex
309  float vertex[3] = {vx, vy, vz};
310  bemc->SetVertex(vertex);
311  // Set threshold
313 
316 
317  // _clusters->Reset(); // !!! Not sure if it is necessarry to do it - ask Chris
318 
319  // make the list of towers above threshold
320  RawTowerContainer::ConstRange begin_end = towers->getTowers();
321  RawTowerContainer::ConstIterator itr = begin_end.first;
322 
323  // Define vector of towers in EmcModule format to input into BEmc
324  EmcModule vhit;
325  std::vector<EmcModule> HitList;
326  HitList.erase(HitList.begin(), HitList.end());
327  int ich;
328 
329  for (; itr != begin_end.second; ++itr)
330  {
331  RawTower *tower = itr->second;
332  // cout << " Tower e = " << tower->get_energy()
333  // << " (" << _min_tower_e << ")" << endl;
334  if (tower->get_energy() > _min_tower_e)
335  {
336  // cout << "(" << tower->get_column() << "," << tower->get_row()
337  // << ") (" << tower->get_binphi() << "," << tower->get_bineta()
338  // << ")" << endl;
339  // ix = tower->get_column();
340  RawTowerDefs::keytype towerid = tower->get_id();
341  int ix = RawTowerDefs::decode_index2(towerid); // index2 is phi in CYL
342  int iy = RawTowerDefs::decode_index1(towerid); // index1 is eta in CYL
343  ix -= BINX0;
344  iy -= BINY0;
345  // ix = tower->get_bineta() - BINX0; // eta: index1
346  // iy = tower->get_binphi() - BINY0; // phi: index2
347  if (ix >= 0 && ix < NBINX && iy >= 0 && iy < NBINY)
348  {
349  ich = iy * NBINX + ix;
350  vhit.ich = ich;
351  vhit.amp = tower->get_energy() * fEnergyNorm; // !!! Global Calibration
352  vhit.tof = 0.;
353  HitList.push_back(vhit);
354  }
355  }
356  }
357 
358  bemc->SetModules(&HitList);
359 
360  // Find clusters (as a set of towers with common edge)
361  int ncl = bemc->FindClusters();
362  if (ncl < 0)
363  {
364  cout << "!!! Error in BEmcRec::FindClusters(): numbers of cluster "
365  << ncl << " ?" << endl;
367  }
368 
369  // Get pointer to clusters
370  std::vector<EmcCluster> *ClusterList = bemc->GetClusters();
371  std::vector<EmcCluster>::iterator pc;
372 
373  std::vector<EmcCluster>::iterator pp;
374  float ecl, ecore, xcg, ycg, xx, xy, yy;
375  // float xcorr, ycorr;
376  EmcModule hmax;
377  RawCluster *cluster;
378 
379  std::vector<EmcCluster> PList;
380  std::vector<EmcModule> Peaks;
381  std::vector<EmcCluster> *pPList = &PList;
382  std::vector<EmcModule> *pPeaks = &Peaks;
383 
384  float prob, chi2;
385  int ndf;
386  float xg, yg, zg;
387 
388  vector<EmcModule>::iterator ph;
389  vector<EmcModule> hlist;
390 
391  // ncl = 0;
392  for (pc = ClusterList->begin(); pc != ClusterList->end(); ++pc)
393  {
394  // ecl = pc->GetTotalEnergy();
395  // pc->GetMoments( &xcg, &ycg, &xx, &xy, &yy );
396 
397  int npk = pc->GetSubClusters(pPList, pPeaks);
398  if (npk < 0) return Fun4AllReturnCodes::ABORTEVENT;
399 
400  // cout << " iCl = " << ncl << " (" << npk << "): E ="
401  // << ecl << " x = " << xcg << " y = " << ycg << endl;
402 
403  for (pp = pPList->begin(); pp != pPList->end(); ++pp)
404  {
405  // Cluster energy
406  ecl = pp->GetTotalEnergy();
407  ecore = pp->GetECoreCorrected();
408  // 3x3 energy around center of gravity
409  //e9 = pp->GetE9();
410  // Ecore (basically near 2x2 energy around center of gravity)
411  //ecore = pp->GetECore();
412  // Center of Gravity etc.
413  pp->GetMoments(xcg, ycg, xx, xy, yy);
414  pp->GetGlobalPos(xg, yg, zg);
415 
416  // Tower with max energy
417  hmax = pp->GetMaxTower();
418 
419  // phi = (xcg-float(NPHI)/2.+0.5)/float(NPHI)*2.*M_PI;
420  // eta = (ycg-float(NETA)/2.+0.5)/float(NETA)*2.2; // -1.1<eta<1.1;
421 
422  // Cell2Abs(towergeom,xcg,ycg,phi,eta);
423 
424  // pp->GetCorrPos(&xcorr, &ycorr);
425  // Cell2Abs(towergeom, xcorr, ycorr, phi, eta);
426  // const double ref_radius = towergeom->get_radius();
427 
428  // phi = 0;
429  // if (phi > M_PI) phi -= 2. * M_PI; // convert to [-pi,pi]]
430 
431  prob = -1;
432  chi2 = 0;
433  ndf = 0;
434  prob = pp->GetProb(chi2, ndf);
435  // cout << "Prob/Chi2/NDF = " << prob << " " << chi2
436  // << " " << ndf << " Ecl = " << ecl << endl;
437 
438  cluster = new RawClusterv1();
439  cluster->set_energy(ecl);
440  cluster->set_ecore(ecore);
441 
442  cluster->set_r(sqrt(xg * xg + yg * yg));
443  cluster->set_phi(atan2(yg, xg));
444  cluster->set_z(zg);
445 
446  cluster->set_prob(prob);
447  if (ndf > 0)
448  cluster->set_chi2(chi2 / ndf);
449  else
450  cluster->set_chi2(0);
451 
452  hlist = pp->GetHitList();
453  ph = hlist.begin();
454  while (ph != hlist.end())
455  {
456  ich = (*ph).ich;
457  int iy = ich / NBINX;
458  int ix = ich % NBINX;
459  // that code needs a closer look - here are the towers
460  // with their energy added to the cluster object where
461  // the id is the tower id
462  // !!!!! Make sure twrkey is correctly extracted
463  // RawTowerDefs::keytype twrkey = RawTowerDefs::encode_towerid(towers->getCalorimeterID(), ix + BINX0, iy + BINY0);
464  RawTowerDefs::keytype twrkey = RawTowerDefs::encode_towerid(towers->getCalorimeterID(), iy + BINY0, ix + BINX0); // Becuase in this part index1 is iy
465  // cout << iphi << " " << ieta << ": "
466  // << twrkey << " e = " << (*ph).amp) << endl;
467  cluster->addTower(twrkey, (*ph).amp / fEnergyNorm);
468  ++ph;
469  }
470 
471  _clusters->AddCluster(cluster);
472  // ncl++;
473 
474  // cout << " ipk = " << ipk << ": E = " << ecl << " E9 = "
475  // << e9 << " x = " << xcg << " y = " << ycg
476  // << " MaxTower: (" << hmax.ich%NPHI << ","
477  // << hmax.ich/NPHI << ") e = " << hmax.amp << endl;
478  }
479  }
480 
482  {
483  double ecluster = _clusters->getTotalEdep();
484  double etower = towers->getTotalEdep();
485  if (ecluster > 0)
486  {
487  if (fabs(etower - ecluster) / ecluster > 1e-9)
488  {
489  cout << "energy conservation violation: ETower: " << etower
490  << " ECluster: " << ecluster
491  << " diff: " << etower - ecluster << endl;
492  }
493  }
494  else
495  {
496  if (etower != 0)
497  {
498  cout << "energy conservation violation: ETower: " << etower
499  << " ECluster: " << ecluster << endl;
500  }
501  }
502  }
504 }
505 
507 {
508  PHNodeIterator iter(topNode);
509 
510  // Grab the cEMC node
511  PHCompositeNode *dstNode = static_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
512  if (!dstNode)
513  {
514  std::cerr << PHWHERE << "DST Node missing, doing nothing." << std::endl;
515  throw std::runtime_error("Failed to find DST node in EmcRawTowerBuilder::CreateNodes");
516  }
517 
518  //Get the _det_name subnode
519  PHCompositeNode *cemcNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", detector));
520 
521  //Check that it is there
522  if (!cemcNode)
523  {
524  cemcNode = new PHCompositeNode(detector);
525  dstNode->addNode(cemcNode);
526  }
527 
529  ClusterNodeName = "CLUSTER_" + detector;
531  cemcNode->addNode(clusterNode);
532 }