File indexing completed on 2025-08-05 08:17:42
0001 #include "PHG4BlockCellReco.h"
0002 #include "PHG4BlockCellGeom.h"
0003 #include "PHG4BlockCellGeomContainer.h"
0004 #include "PHG4BlockGeom.h"
0005 #include "PHG4BlockGeomContainer.h"
0006 #include "PHG4Cell.h" // for PHG4Cell, PHG...
0007 #include "PHG4CellContainer.h"
0008 #include "PHG4Cellv1.h"
0009
0010 #include "PHG4CellDefs.h"
0011
0012 #include <phparameter/PHParameterContainerInterface.h> // for PHParameterCo...
0013 #include <phparameter/PHParametersContainer.h>
0014
0015 #include <g4main/PHG4Hit.h>
0016 #include <g4main/PHG4HitContainer.h>
0017 #include <g4main/PHG4Utils.h>
0018
0019 #include <fun4all/Fun4AllReturnCodes.h>
0020 #include <fun4all/SubsysReco.h> // for SubsysReco
0021
0022 #include <phool/PHCompositeNode.h>
0023 #include <phool/PHIODataNode.h>
0024 #include <phool/PHNode.h> // for PHNode
0025 #include <phool/PHNodeIterator.h>
0026 #include <phool/PHObject.h> // for PHObject
0027 #include <phool/getClass.h>
0028 #include <phool/phool.h> // for PHWHERE
0029
0030 #include <cmath>
0031 #include <cstdlib>
0032 #include <iostream>
0033 #include <sstream>
0034 #include <vector> // for vector
0035
0036 using namespace std;
0037
0038 static vector<PHG4Cell *> cellptarray;
0039
0040 PHG4BlockCellReco::PHG4BlockCellReco(const string &name)
0041 : SubsysReco(name)
0042 , PHParameterContainerInterface(name)
0043 , sum_energy_g4hit(0)
0044 , chkenergyconservation(0)
0045 {
0046 SetDefaultParameters();
0047 }
0048
0049 int PHG4BlockCellReco::ResetEvent(PHCompositeNode * )
0050 {
0051 sum_energy_g4hit = 0.;
0052 return Fun4AllReturnCodes::EVENT_OK;
0053 }
0054
0055 int PHG4BlockCellReco::InitRun(PHCompositeNode *topNode)
0056 {
0057 PHNodeIterator iter(topNode);
0058
0059
0060 PHCompositeNode *dstNode;
0061 dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
0062 if (!dstNode)
0063 {
0064 cout << Name() << " DST Node missing, doing nothing." << std::endl;
0065 exit(1);
0066 }
0067 PHNodeIterator dstiter(dstNode);
0068 PHCompositeNode *DetNode =
0069 dynamic_cast<PHCompositeNode *>(dstiter.findFirst("PHCompositeNode",
0070 detector));
0071 if (!DetNode)
0072 {
0073 DetNode = new PHCompositeNode(detector);
0074 dstNode->addNode(DetNode);
0075 }
0076
0077 PHCompositeNode *runNode;
0078 runNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "RUN"));
0079 if (!runNode)
0080 {
0081 cout << Name() << "DST Node missing, doing nothing." << endl;
0082 exit(1);
0083 }
0084 PHNodeIterator runiter(runNode);
0085 PHCompositeNode *RunDetNode =
0086 dynamic_cast<PHCompositeNode *>(runiter.findFirst("PHCompositeNode",
0087 detector));
0088 if (!RunDetNode)
0089 {
0090 RunDetNode = new PHCompositeNode(detector);
0091 runNode->addNode(RunDetNode);
0092 }
0093
0094 hitnodename = "G4HIT_" + detector;
0095 PHG4HitContainer *g4hit = findNode::getClass<PHG4HitContainer>(topNode, hitnodename);
0096 if (!g4hit)
0097 {
0098 cout << Name() << " Could not locate g4 hit node " << hitnodename << endl;
0099 exit(1);
0100 }
0101
0102 cellnodename = "G4CELL_" + detector;
0103 PHG4CellContainer *cells = findNode::getClass<PHG4CellContainer>(topNode, cellnodename);
0104 if (!cells)
0105 {
0106 cells = new PHG4CellContainer();
0107 PHIODataNode<PHObject> *newNode = new PHIODataNode<PHObject>(cells, cellnodename, "PHObject");
0108 DetNode->addNode(newNode);
0109 }
0110
0111 geonodename = "BLOCKGEOM_" + detector;
0112 PHG4BlockGeomContainer *geo = findNode::getClass<PHG4BlockGeomContainer>(topNode, geonodename);
0113 if (!geo)
0114 {
0115 cout << Name() << " Could not locate geometry node " << geonodename << endl;
0116 exit(1);
0117 }
0118
0119 if (Verbosity() > 0)
0120 {
0121 geo->identify();
0122 }
0123
0124 seggeonodename = "BLOCKCELLGEOM_" + detector;
0125 PHG4BlockCellGeomContainer *seggeo = findNode::getClass<PHG4BlockCellGeomContainer>(topNode, seggeonodename);
0126 if (!seggeo)
0127 {
0128 seggeo = new PHG4BlockCellGeomContainer();
0129 PHIODataNode<PHObject> *newNode = new PHIODataNode<PHObject>(seggeo, seggeonodename, "PHObject");
0130 runNode->addNode(newNode);
0131 }
0132 GetParamsContainerModify()->set_name(detector);
0133
0134 UpdateParametersWithMacro();
0135
0136 map<int, PHG4BlockGeom *>::const_iterator miter;
0137 pair<map<int, PHG4BlockGeom *>::const_iterator, map<int, PHG4BlockGeom *>::const_iterator> begin_end = geo->get_begin_end();
0138 map<int, std::pair<double, double> >::iterator sizeiter;
0139 for (miter = begin_end.first; miter != begin_end.second; ++miter)
0140 {
0141 PHG4BlockGeom *layergeom = miter->second;
0142 int layer = layergeom->get_layer();
0143 if (!ExistDetid(layer))
0144 {
0145 cout << Name() << ": No parameters for detid/layer " << layer
0146 << ", hits from this detid/layer will not be accumulated into cells" << endl;
0147 continue;
0148 }
0149 implemented_detid.insert(layer);
0150 double radius = sqrt(pow(layergeom->get_center_x(), 2) + pow(layergeom->get_center_y(), 2));
0151 double width = layergeom->get_size_x();
0152 double length_in_z = layergeom->get_size_z();
0153 double zmin = layergeom->get_center_z() - length_in_z / 2.;
0154 double zmax = zmin + length_in_z;
0155 set_size(layer, get_double_param(layer, "deltaeta"), get_double_param(layer, "deltax"), PHG4CellDefs::etaphibinning);
0156 tmin_max.insert(std::make_pair(layer, std::make_pair(get_double_param(layer, "tmin"), get_double_param(layer, "tmax"))));
0157 sizeiter = cell_size.find(layer);
0158 if (sizeiter == cell_size.end())
0159 {
0160 cout << Name() << "no cell sizes for layer " << layer << endl;
0161 exit(1);
0162 }
0163
0164
0165 PHG4BlockCellGeom *layerseggeo = new PHG4BlockCellGeom();
0166 layerseggeo->set_layer(layergeom->get_layer());
0167
0168
0169
0170 if (binning[layer] == PHG4CellDefs::etaphibinning)
0171 {
0172
0173
0174 double etamin = PHG4Utils::get_eta(radius + 0.5 * layergeom->get_size_y(), zmin);
0175 double etamax = PHG4Utils::get_eta(radius + 0.5 * layergeom->get_size_y(), zmax);
0176 zmin_max[layer] = make_pair(etamin, etamax);
0177 double etastepsize = (sizeiter->second).first;
0178 double d_etabins;
0179 double fract = modf((etamax - etamin) / etastepsize, &d_etabins);
0180 if (fract != 0)
0181 {
0182 d_etabins++;
0183 }
0184 etastepsize = (etamax - etamin) / d_etabins;
0185 (sizeiter->second).first = etastepsize;
0186 int etabins = d_etabins;
0187 double etahi = etamin + etastepsize;
0188 for (int i = 0; i < etabins; i++)
0189 {
0190 if (etahi > (etamax + 1.e-6))
0191 {
0192 cout << "etahi: " << etahi << ", etamax: " << etamax << endl;
0193 }
0194 etahi += etastepsize;
0195 }
0196
0197 double xmin = -layergeom->get_width() / 2.;
0198
0199 double xstepsize = (sizeiter->second).second;
0200 double d_xbins;
0201 fract = modf(width / xstepsize, &d_xbins);
0202 if (fract != 0)
0203 {
0204 d_xbins++;
0205 }
0206
0207 xstepsize = width / d_xbins;
0208 (sizeiter->second).second = xstepsize;
0209 int xbins = d_xbins;
0210
0211
0212
0213
0214
0215
0216
0217
0218
0219
0220
0221
0222
0223 pair<int, int> x_z_bin = make_pair(xbins, etabins);
0224 n_x_z_bins[layer] = x_z_bin;
0225 layerseggeo->set_binning(PHG4CellDefs::etaphibinning);
0226 layerseggeo->set_etabins(etabins);
0227 layerseggeo->set_etamin(etamin);
0228 layerseggeo->set_etastep(etastepsize);
0229 layerseggeo->set_xmin(xmin);
0230 layerseggeo->set_xbins(xbins);
0231 layerseggeo->set_xstep(xstepsize);
0232 xstep[layer] = xstepsize;
0233 etastep[layer] = etastepsize;
0234 }
0235
0236
0237 seggeo->AddLayerCellGeom(layerseggeo);
0238 if (Verbosity() > 1)
0239 {
0240 layerseggeo->identify();
0241 }
0242 }
0243 string nodename = "G4CELLPARAM_" + GetParamsContainer()->Name();
0244 SaveToNodeTree(RunDetNode, nodename);
0245 return Fun4AllReturnCodes::EVENT_OK;
0246 }
0247
0248 int PHG4BlockCellReco::process_event(PHCompositeNode *topNode)
0249 {
0250 PHG4HitContainer *g4hit = findNode::getClass<PHG4HitContainer>(topNode, hitnodename);
0251 if (!g4hit)
0252 {
0253 cout << "Could not locate g4 hit node " << hitnodename << endl;
0254 exit(1);
0255 }
0256 PHG4CellContainer *cells = findNode::getClass<PHG4CellContainer>(topNode, cellnodename);
0257 if (!cells)
0258 {
0259 cout << "could not locate cell node " << cellnodename << endl;
0260 exit(1);
0261 }
0262
0263 PHG4BlockCellGeomContainer *seggeo = findNode::getClass<PHG4BlockCellGeomContainer>(topNode, seggeonodename);
0264 if (!seggeo)
0265 {
0266 cout << "could not locate geo node " << seggeonodename << endl;
0267 exit(1);
0268 }
0269
0270 PHG4HitContainer::LayerIter layer;
0271 pair<PHG4HitContainer::LayerIter, PHG4HitContainer::LayerIter> layer_begin_end = g4hit->getLayers();
0272
0273
0274
0275
0276
0277
0278
0279 for (layer = layer_begin_end.first; layer != layer_begin_end.second; layer++)
0280 {
0281
0282 if (implemented_detid.find(*layer) == implemented_detid.end())
0283 {
0284 continue;
0285 }
0286 PHG4HitContainer::ConstIterator hiter;
0287 PHG4HitContainer::ConstRange hit_begin_end = g4hit->getHits(*layer);
0288 PHG4BlockCellGeom *geo = seggeo->GetLayerCellGeom(*layer);
0289 int nxbins = n_x_z_bins[*layer].first;
0290 int nzbins = n_x_z_bins[*layer].second;
0291 unsigned int nbins = nxbins * nzbins;
0292
0293 if (cellptarray.size() < nbins)
0294 {
0295 cellptarray.resize(nbins, nullptr);
0296 }
0297
0298
0299 if (binning[*layer] == PHG4CellDefs::etaphibinning)
0300 {
0301 for (hiter = hit_begin_end.first; hiter != hit_begin_end.second; hiter++)
0302 {
0303
0304 if (hiter->second->get_t(0) > tmin_max[*layer].second)
0305 {
0306 continue;
0307 }
0308 if (hiter->second->get_t(1) < tmin_max[*layer].first)
0309 {
0310 continue;
0311 }
0312
0313 pair<double, double> etax[2];
0314 double xbin[2];
0315 double etabin[2];
0316 for (int i = 0; i < 2; i++)
0317 {
0318 etax[i] = PHG4Utils::get_etaphi(hiter->second->get_x(i), hiter->second->get_y(i), hiter->second->get_z(i));
0319 etabin[i] = geo->get_etabin(etax[i].first);
0320 xbin[i] = geo->get_xbin(etax[i].second);
0321 }
0322
0323
0324 if (xbin[0] < 0 || xbin[0] >= nxbins || xbin[1] < 0 || xbin[1] >= nxbins)
0325 {
0326 continue;
0327 }
0328 if (etabin[0] < 0 || etabin[0] >= nzbins || etabin[1] < 0 || etabin[1] >= nzbins)
0329 {
0330 continue;
0331 }
0332
0333 if (etabin[0] < 0)
0334 {
0335 if (Verbosity() > 0)
0336 {
0337 hiter->second->identify();
0338 }
0339 continue;
0340 }
0341
0342 sum_energy_g4hit += hiter->second->get_edep();
0343 int intxbin = xbin[0];
0344 int intetabin = etabin[0];
0345 int intxbinout = xbin[1];
0346 int intetabinout = etabin[1];
0347
0348
0349
0350 double ax = (etax[0]).second;
0351 double ay = (etax[0]).first;
0352 double bx = (etax[1]).second;
0353 double by = (etax[1]).first;
0354 if (intxbin > intxbinout)
0355 {
0356 int tmp = intxbin;
0357 intxbin = intxbinout;
0358 intxbinout = tmp;
0359 }
0360
0361 if (intetabin > intetabinout)
0362 {
0363 int tmp = intetabin;
0364 intetabin = intetabinout;
0365 intetabinout = tmp;
0366 }
0367
0368 double trklen = sqrt((ax - bx) * (ax - bx) + (ay - by) * (ay - by));
0369
0370
0371
0372
0373
0374
0375 if (trklen == 0)
0376 {
0377 trklen = -1.;
0378 }
0379 vector<int> vx;
0380 vector<int> veta;
0381 vector<double> vdedx;
0382
0383 if (intxbin == intxbinout && intetabin == intetabinout)
0384 {
0385 if (Verbosity() > 0)
0386 {
0387 cout << "SINGLE CELL FIRED: " << intxbin << " " << intetabin << endl;
0388 }
0389 vx.push_back(intxbin);
0390 veta.push_back(intetabin);
0391 vdedx.push_back(trklen);
0392 }
0393 else
0394 {
0395 for (int ibp = intxbin; ibp <= intxbinout; ibp++)
0396 {
0397 for (int ibz = intetabin; ibz <= intetabinout; ibz++)
0398 {
0399 double cx = geo->get_xcenter(ibp) - geo->get_xstep() / 2.;
0400 double dx = geo->get_xcenter(ibp) + geo->get_xstep() / 2.;
0401 double cy = geo->get_etacenter(ibz) - geo->get_etastep() / 2.;
0402 double dy = geo->get_etacenter(ibz) + geo->get_etastep() / 2.;
0403
0404
0405 pair<bool, double> intersect = PHG4Utils::line_and_rectangle_intersect(ax, ay, bx, by, cx, cy, dx, dy);
0406 if (intersect.first)
0407 {
0408 if (Verbosity() > 0)
0409 {
0410 cout << "CELL FIRED: " << ibp << " " << ibz << " " << intersect.second << endl;
0411 }
0412 vx.push_back(ibp);
0413 veta.push_back(ibz);
0414 vdedx.push_back(intersect.second);
0415 }
0416 }
0417 }
0418 }
0419 if (Verbosity() > 0)
0420 {
0421 cout << "NUMBER OF FIRED CELLS = " << vx.size() << endl;
0422 }
0423
0424 double tmpsum = 0.;
0425 for (unsigned int ii = 0; ii < vx.size(); ii++)
0426 {
0427 tmpsum += vdedx[ii];
0428 vdedx[ii] = vdedx[ii] / trklen;
0429 if (Verbosity() > 0)
0430 {
0431 cout << " CELL " << ii << " dE/dX = " << vdedx[ii] << endl;
0432 }
0433 }
0434 if (Verbosity() > 0)
0435 {
0436 cout << " TOTAL TRACK LENGTH = " << tmpsum << " " << trklen << endl;
0437 }
0438
0439 for (unsigned int i1 = 0; i1 < vx.size(); i1++)
0440 {
0441 int ixbin = vx[i1];
0442 int ietabin = veta[i1];
0443 int ibin = ixbin * nzbins + ietabin;
0444 if (!cellptarray[ibin])
0445 {
0446 PHG4CellDefs::keytype key = PHG4CellDefs::EtaXsizeBinning::genkey(*layer, ixbin, ietabin);
0447 cellptarray[ibin] = new PHG4Cellv1(key);
0448 }
0449 cellptarray[ibin]->add_edep(hiter->first, hiter->second->get_edep() * vdedx[i1]);
0450 cellptarray[ibin]->add_edep(hiter->second->get_edep() * vdedx[i1]);
0451 if (hiter->second->has_property(PHG4Hit::prop_light_yield))
0452 {
0453 cellptarray[ibin]->add_light_yield(hiter->second->get_light_yield() * vdedx[i1]);
0454 }
0455 cellptarray[ibin]->add_shower_edep(hiter->second->get_shower_id(), hiter->second->get_edep() * vdedx[i1]);
0456
0457
0458 if (!isfinite(hiter->second->get_edep() * vdedx[i1]))
0459 {
0460 cout << PHWHERE << " invalid energy dep " << hiter->second->get_edep()
0461 << " or path length: " << vdedx[i1] << endl;
0462 }
0463 }
0464
0465 vx.clear();
0466 veta.clear();
0467 }
0468
0469 int numcells = 0;
0470 for (int ix = 0; ix < nxbins; ix++)
0471 {
0472 for (int iz = 0; iz < nzbins; iz++)
0473 {
0474 int ibin = ix * nzbins + iz;
0475
0476 if (cellptarray[ibin])
0477 {
0478 cells->AddCell(cellptarray[ibin]);
0479 numcells++;
0480 if (Verbosity() > 1)
0481 {
0482 cout << "Adding cell in bin x: " << ix
0483 << " x: " << geo->get_xcenter(ix) * 180. / M_PI
0484 << ", eta bin: " << iz
0485 << ", eta: " << geo->get_etacenter(iz)
0486 << ", energy dep: " << cellptarray[ibin]->get_edep()
0487 << endl;
0488 }
0489
0490 cellptarray[ibin] = nullptr;
0491 }
0492 }
0493 }
0494
0495 if (Verbosity() > 0)
0496 {
0497 cout << Name() << ": found " << numcells << " eta/x cells with energy deposition" << endl;
0498 }
0499 }
0500 }
0501
0502 if (chkenergyconservation)
0503 {
0504 CheckEnergy(topNode);
0505 }
0506
0507 return Fun4AllReturnCodes::EVENT_OK;
0508 }
0509
0510 void PHG4BlockCellReco::etaxsize(const int detid, const double deltaeta, const double deltax)
0511 {
0512 set_double_param(detid, "deltaeta", deltaeta);
0513 set_double_param(detid, "deltax", deltax);
0514 return;
0515 }
0516
0517 void PHG4BlockCellReco::set_timing_window(const int detid, const double tmin, const double tmax)
0518 {
0519 set_double_param(detid, "tmin", tmin);
0520 set_double_param(detid, "tmax", tmax);
0521 return;
0522 }
0523
0524 void PHG4BlockCellReco::set_size(const int i, const double sizeA, const double sizeB, const int what)
0525 {
0526 if (binning.find(i) != binning.end())
0527 {
0528 cout << "size for layer " << i << " already set" << endl;
0529 return;
0530 }
0531
0532 binning[i] = what;
0533 cell_size[i] = std::make_pair(sizeA, sizeB);
0534
0535 return;
0536 }
0537
0538
0539
0540 int PHG4BlockCellReco::CheckEnergy(PHCompositeNode *topNode)
0541 {
0542 PHG4CellContainer *cells = findNode::getClass<PHG4CellContainer>(topNode, cellnodename);
0543 double sum_energy_cells = 0.;
0544 double sum_energy_stored_hits = 0.;
0545 double sum_energy_stored_showers = 0.;
0546
0547 PHG4CellContainer::ConstRange cell_begin_end = cells->getCells();
0548 PHG4CellContainer::ConstIterator citer;
0549 for (citer = cell_begin_end.first; citer != cell_begin_end.second; ++citer)
0550 {
0551 sum_energy_cells += citer->second->get_edep();
0552 PHG4Cell::EdepConstRange cellrange = citer->second->get_g4hits();
0553 for (PHG4Cell::EdepConstIterator iter = cellrange.first; iter != cellrange.second; ++iter)
0554 {
0555 sum_energy_stored_hits += iter->second;
0556 }
0557 PHG4Cell::ShowerEdepConstRange shwrrange = citer->second->get_g4showers();
0558 for (PHG4Cell::ShowerEdepConstIterator iter = shwrrange.first; iter != shwrrange.second; ++iter)
0559 {
0560 sum_energy_stored_showers += iter->second;
0561 }
0562 }
0563
0564
0565 if (sum_energy_stored_hits > 0)
0566 {
0567 if (fabs(sum_energy_cells - sum_energy_stored_hits) / sum_energy_cells > 1e-6)
0568 {
0569 cout << "energy mismatch between cell energy " << sum_energy_cells
0570 << " and stored hit energies " << sum_energy_stored_hits
0571 << endl;
0572 }
0573 }
0574 if (sum_energy_stored_showers > 0)
0575 {
0576 if (fabs(sum_energy_cells - sum_energy_stored_showers) / sum_energy_cells > 1e-6)
0577 {
0578 cout << "energy mismatch between cell energy " << sum_energy_cells
0579 << " and stored shower energies " << sum_energy_stored_showers
0580 << endl;
0581 }
0582 }
0583
0584 if (fabs(sum_energy_cells - sum_energy_g4hit) / sum_energy_g4hit > 1e-6)
0585 {
0586 cout << "energy mismatch between cells: " << sum_energy_cells
0587 << " and hits: " << sum_energy_g4hit
0588 << " diff sum(cells) - sum(hits): " << sum_energy_cells - sum_energy_g4hit
0589 << endl;
0590 return -1;
0591 }
0592
0593 else
0594 {
0595 if (Verbosity() > 0)
0596 {
0597 cout << Name() << ": sum hit energy: " << sum_energy_g4hit << " GeV" << endl;
0598 cout << Name() << ": sum cell energy: " << sum_energy_cells << " GeV" << endl;
0599 cout << Name() << ": sum shower energy: " << sum_energy_stored_showers << " GeV" << endl;
0600 cout << Name() << ": sum stored hit energy: " << sum_energy_stored_hits << " GeV" << endl;
0601 }
0602 }
0603
0604 return 0;
0605 }
0606
0607 void PHG4BlockCellReco::SetDefaultParameters()
0608 {
0609 set_default_double_param("deltaeta", NAN);
0610 set_default_double_param("deltax", NAN);
0611 set_default_double_param("tmax", 60.0);
0612 set_default_double_param("tmin", 0.0);
0613 return;
0614 }