File indexing completed on 2025-08-06 08:18:59
0001 #include "PHG4FullProjSpacalCellReco.h"
0002
0003 #include "LightCollectionModel.h"
0004 #include "PHG4Cell.h" // for PHG4Cell
0005 #include "PHG4CylinderCellGeom.h"
0006 #include "PHG4CylinderCellGeomContainer.h"
0007 #include "PHG4CylinderCellGeom_Spacalv1.h"
0008 #include "PHG4CylinderGeom.h" // for PHG4CylinderGeom
0009 #include "PHG4CylinderGeomContainer.h"
0010 #include "PHG4CylinderGeom_Spacalv1.h" // for PHG4CylinderGeom_Spaca...
0011 #include "PHG4CylinderGeom_Spacalv3.h"
0012
0013 #include "PHG4CellContainer.h"
0014 #include "PHG4CellDefs.h"
0015 #include "PHG4Cellv1.h"
0016
0017 #include <phparameter/PHParameterInterface.h> // for PHParameterInterface
0018
0019 #include <g4main/PHG4Hit.h>
0020 #include <g4main/PHG4HitContainer.h>
0021
0022 #include <fun4all/Fun4AllBase.h> // for Fun4AllBase::VERBOSITY...
0023 #include <fun4all/Fun4AllReturnCodes.h>
0024 #include <fun4all/Fun4AllServer.h>
0025 #include <fun4all/SubsysReco.h> // for SubsysReco
0026
0027 #include <phool/PHCompositeNode.h>
0028 #include <phool/PHIODataNode.h>
0029 #include <phool/PHNode.h> // for PHNode
0030 #include <phool/PHNodeIterator.h>
0031 #include <phool/PHObject.h> // for PHObject
0032 #include <phool/getClass.h>
0033 #include <phool/phool.h> // for PHWHERE
0034 #include <phool/recoConsts.h>
0035
0036 #include <ffamodules/CDBInterface.h>
0037
0038 #include <TAxis.h> // for TAxis
0039 #include <TFile.h>
0040 #include <TH1.h>
0041 #include <TH2.h>
0042 #include <TObject.h> // for TObject
0043 #include <TSystem.h>
0044
0045 #include <cassert>
0046 #include <cmath>
0047 #include <cstdlib>
0048 #include <exception>
0049 #include <iostream>
0050 #include <sstream>
0051 #include <utility> // for pair
0052
0053 PHG4FullProjSpacalCellReco::PHG4FullProjSpacalCellReco(const std::string &name)
0054 : SubsysReco(name)
0055 , PHParameterInterface(name)
0056 {
0057 InitializeParameters();
0058 }
0059
0060 int PHG4FullProjSpacalCellReco::ResetEvent(PHCompositeNode * )
0061 {
0062 sum_energy_g4hit = 0.;
0063 return Fun4AllReturnCodes::EVENT_OK;
0064 }
0065
0066 int PHG4FullProjSpacalCellReco::InitRun(PHCompositeNode *topNode)
0067 {
0068 PHNodeIterator iter(topNode);
0069
0070
0071 PHCompositeNode *dstNode;
0072 dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
0073 if (!dstNode)
0074 {
0075 std::cout << PHWHERE << "DST Node missing, doing nothing." << std::endl;
0076 exit(1);
0077 }
0078 hitnodename = "G4HIT_" + detector;
0079 PHG4HitContainer *g4hit = findNode::getClass<PHG4HitContainer>(topNode, hitnodename);
0080 if (!g4hit)
0081 {
0082 std::cout << "PHG4FullProjSpacalCellReco::InitRun - Could not locate g4 hit node "
0083 << hitnodename << std::endl;
0084 topNode->print();
0085 exit(1);
0086 }
0087 cellnodename = "G4CELL_" + detector;
0088 PHG4CellContainer *cells = findNode::getClass<PHG4CellContainer>(topNode, cellnodename);
0089 if (!cells)
0090 {
0091 PHNodeIterator dstiter(dstNode);
0092 PHCompositeNode *DetNode = dynamic_cast<PHCompositeNode *>(dstiter.findFirst("PHCompositeNode", detector));
0093 if (!DetNode)
0094 {
0095 DetNode = new PHCompositeNode(detector);
0096 dstNode->addNode(DetNode);
0097 }
0098 cells = new PHG4CellContainer();
0099 PHIODataNode<PHObject> *newNode = new PHIODataNode<PHObject>(cells, cellnodename, "PHObject");
0100 DetNode->addNode(newNode);
0101 }
0102
0103 geonodename = "CYLINDERGEOM_" + detector;
0104 PHG4CylinderGeomContainer *geo = findNode::getClass<PHG4CylinderGeomContainer>(topNode, geonodename);
0105 if (!geo)
0106 {
0107 std::cout << "PHG4FullProjSpacalCellReco::InitRun - Could not locate geometry node "
0108 << geonodename << std::endl;
0109 topNode->print();
0110 exit(1);
0111 }
0112 if (Verbosity() > 0)
0113 {
0114 std::cout << "PHG4FullProjSpacalCellReco::InitRun - incoming geometry:"
0115 << std::endl;
0116 geo->identify();
0117 }
0118 seggeonodename = "CYLINDERCELLGEOM_" + detector;
0119 PHG4CylinderCellGeomContainer *seggeo = findNode::getClass<PHG4CylinderCellGeomContainer>(topNode, seggeonodename);
0120 if (!seggeo)
0121 {
0122 seggeo = new PHG4CylinderCellGeomContainer();
0123 PHCompositeNode *runNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "RUN"));
0124 PHIODataNode<PHObject> *newNode = new PHIODataNode<PHObject>(seggeo, seggeonodename, "PHObject");
0125 runNode->addNode(newNode);
0126 }
0127
0128 const PHG4CylinderGeom *layergeom_raw = geo->GetFirstLayerGeom();
0129 assert(layergeom_raw);
0130
0131
0132 const PHG4CylinderGeom_Spacalv3 *layergeom =
0133 dynamic_cast<const PHG4CylinderGeom_Spacalv3 *>(layergeom_raw);
0134
0135 if (!layergeom)
0136 {
0137 std::cout << "PHG4FullProjSpacalCellReco::InitRun - Fatal Error -"
0138 << " require to work with a version of SPACAL geometry of PHG4CylinderGeom_Spacalv3 or higher. "
0139 << "However the incoming geometry has version "
0140 << layergeom_raw->ClassName() << std::endl;
0141 exit(1);
0142 }
0143 if (Verbosity() > 1)
0144 {
0145 layergeom->identify();
0146 }
0147
0148 layergeom->subtower_consistency_check();
0149
0150
0151
0152
0153 PHG4CylinderCellGeom_Spacalv1 *layerseggeo = new PHG4CylinderCellGeom_Spacalv1();
0154 layerseggeo->set_layer(layergeom->get_layer());
0155 layerseggeo->set_radius(layergeom->get_radius());
0156 layerseggeo->set_thickness(layergeom->get_thickness());
0157 layerseggeo->set_binning(PHG4CellDefs::spacalbinning);
0158
0159
0160
0161 const PHG4CylinderGeom_Spacalv3::tower_map_t &tower_map = layergeom->get_sector_tower_map();
0162 const PHG4CylinderGeom_Spacalv3::sector_map_t §or_map = layergeom->get_sector_map();
0163 const int nphibin = layergeom->get_azimuthal_n_sec()
0164 * layergeom->get_max_phi_bin_in_sec()
0165 * layergeom->get_n_subtower_phi();
0166 const double deltaphi = 2. * M_PI / nphibin;
0167
0168 using map_z_tower_z_ID_t = std::map<double, int>;
0169 map_z_tower_z_ID_t map_z_tower_z_ID;
0170 double phi_min = NAN;
0171
0172 for (const auto &tower_pair : tower_map)
0173 {
0174 const int &tower_ID = tower_pair.first;
0175 const PHG4CylinderGeom_Spacalv3::geom_tower &tower = tower_pair.second;
0176
0177
0178 std::pair<int, int> tower_z_phi_ID = layergeom->get_tower_z_phi_ID(tower_ID, 0);
0179
0180 const int &tower_ID_z = tower_z_phi_ID.first;
0181 const int &tower_ID_phi = tower_z_phi_ID.second;
0182
0183 if (tower_ID_phi == 0)
0184 {
0185
0186
0187 phi_min = M_PI_2 - deltaphi * (layergeom->get_max_phi_bin_in_sec() * layergeom->get_n_subtower_phi() / 2)
0188 + sector_map.begin()->second;
0189 }
0190
0191 if (tower_ID_phi == layergeom->get_max_phi_bin_in_sec() / 2)
0192 {
0193
0194 map_z_tower_z_ID[tower.centralZ] = tower_ID_z;
0195 }
0196
0197 }
0198
0199 assert(!std::isnan(phi_min));
0200 layerseggeo->set_phimin(phi_min);
0201 layerseggeo->set_phistep(deltaphi);
0202 layerseggeo->set_phibins(nphibin);
0203
0204 PHG4CylinderCellGeom_Spacalv1::tower_z_ID_eta_bin_map_t tower_z_ID_eta_bin_map;
0205 int eta_bin = 0;
0206 for (auto &z_tower_z_ID : map_z_tower_z_ID)
0207 {
0208 tower_z_ID_eta_bin_map[z_tower_z_ID.second] = eta_bin;
0209 eta_bin++;
0210 }
0211 layerseggeo->set_tower_z_ID_eta_bin_map(tower_z_ID_eta_bin_map);
0212 layerseggeo->set_etabins(eta_bin * layergeom->get_n_subtower_eta());
0213 layerseggeo->set_etamin(NAN);
0214 layerseggeo->set_etastep(NAN);
0215
0216
0217 for (const auto &tower_pair : tower_map)
0218 {
0219 const int &tower_ID = tower_pair.first;
0220 const PHG4CylinderGeom_Spacalv3::geom_tower &tower = tower_pair.second;
0221
0222
0223 std::pair<int, int> tower_z_phi_ID = layergeom->get_tower_z_phi_ID(tower_ID, 0);
0224 const int &tower_ID_z = tower_z_phi_ID.first;
0225 const int &tower_ID_phi = tower_z_phi_ID.second;
0226 const int &etabin = tower_z_ID_eta_bin_map[tower_ID_z];
0227
0228 if (tower_ID_phi == layergeom->get_max_phi_bin_in_sec() / 2)
0229 {
0230
0231 const double dz = fabs(0.5 * (tower.pDy1 + tower.pDy2) / sin(tower.pRotationAngleX));
0232 const double tower_radial = layergeom->get_tower_radial_position(tower);
0233
0234 auto z_to_eta = [&tower_radial](const double &z)
0235 { return -log(tan(0.5 * atan2(tower_radial, z))); };
0236
0237 const double eta_central = z_to_eta(tower.centralZ);
0238
0239 const double deta = (z_to_eta(tower.centralZ + dz) - z_to_eta(tower.centralZ - dz)) / 2;
0240 assert(deta > 0);
0241
0242 for (int sub_tower_ID_y = 0; sub_tower_ID_y < tower.NSubtowerY;
0243 ++sub_tower_ID_y)
0244 {
0245 assert(tower.NSubtowerY <= layergeom->get_n_subtower_eta());
0246
0247 const int sub_tower_etabin = etabin * layergeom->get_n_subtower_eta() + sub_tower_ID_y;
0248
0249 const std::pair<double, double> etabounds(eta_central - deta + sub_tower_ID_y * 2 * deta / tower.NSubtowerY,
0250 eta_central - deta + (sub_tower_ID_y + 1) * 2 * deta / tower.NSubtowerY);
0251
0252 const std::pair<double, double> zbounds(tower.centralZ - dz + sub_tower_ID_y * 2 * dz / tower.NSubtowerY,
0253 tower.centralZ - dz + (sub_tower_ID_y + 1) * 2 * dz / tower.NSubtowerY);
0254
0255 layerseggeo->set_etabounds(sub_tower_etabin, etabounds);
0256 layerseggeo->set_zbounds(sub_tower_etabin, zbounds);
0257
0258 if (Verbosity() >= VERBOSITY_SOME)
0259 {
0260 std::cout << "PHG4FullProjSpacalCellReco::InitRun::" << Name()
0261 << "\t tower_ID_z = " << tower_ID_z
0262 << "\t tower_ID_phi = " << tower_ID_phi
0263 << "\t sub_tower_ID_y = " << sub_tower_ID_y
0264 << "\t sub_tower_etabin = " << sub_tower_etabin
0265 << "\t dz = " << dz
0266 << "\t tower_radial = " << tower_radial
0267 << "\t eta_central = " << eta_central
0268 << "\t deta = " << deta
0269 << "\t etabounds = [" << etabounds.first << ", " << etabounds.second << "]"
0270 << "\t zbounds = [" << zbounds.first << ", " << zbounds.second << "]"
0271 << std::endl;
0272 }
0273 }
0274 }
0275
0276 }
0277
0278
0279 seggeo->AddLayerCellGeom(layerseggeo);
0280 if (Verbosity() >= VERBOSITY_SOME)
0281 {
0282 std::cout << "PHG4FullProjSpacalCellReco::InitRun::" << Name()
0283 << " - Done layer" << (layergeom->get_layer())
0284 << ". Print out the cell geometry:" << std::endl;
0285 layerseggeo->identify();
0286 }
0287 UpdateParametersWithMacro();
0288
0289 PHCompositeNode *runNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "RUN"));
0290 PHCompositeNode *parNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "PAR"));
0291 PHNodeIterator runIter(runNode);
0292 PHCompositeNode *RunDetNode = dynamic_cast<PHCompositeNode *>(runIter.findFirst("PHCompositeNode", detector));
0293 if (!RunDetNode)
0294 {
0295 RunDetNode = new PHCompositeNode(detector);
0296 runNode->addNode(RunDetNode);
0297 }
0298 std::string paramnodename = "G4CELLPARAM_" + detector;
0299 SaveToNodeTree(RunDetNode, paramnodename);
0300
0301 PHNodeIterator parIter(parNode);
0302 PHCompositeNode *ParDetNode = dynamic_cast<PHCompositeNode *>(parIter.findFirst("PHCompositeNode", detector));
0303 if (!ParDetNode)
0304 {
0305 ParDetNode = new PHCompositeNode(detector);
0306 parNode->addNode(ParDetNode);
0307 }
0308 std::string cellgeonodename = "G4CELLGEO_" + detector;
0309 PutOnParNode(ParDetNode, cellgeonodename);
0310 tmin = get_double_param("tmin");
0311 tmax = get_double_param("tmax");
0312 m_DeltaT = get_double_param("delta_t");
0313 return Fun4AllReturnCodes::EVENT_OK;
0314 }
0315
0316 int PHG4FullProjSpacalCellReco::process_event(PHCompositeNode *topNode)
0317 {
0318 PHG4HitContainer *g4hit = findNode::getClass<PHG4HitContainer>(topNode, hitnodename);
0319 if (!g4hit)
0320 {
0321 std::cout << "PHG4FullProjSpacalCellReco::process_event - Fatal Error - Could not locate g4 hit node "
0322 << hitnodename << std::endl;
0323 exit(1);
0324 }
0325 PHG4CellContainer *cells = findNode::getClass<PHG4CellContainer>(topNode, cellnodename);
0326 if (!cells)
0327 {
0328 std::cout << "PHG4FullProjSpacalCellReco::process_event - Fatal Error - could not locate cell node "
0329 << cellnodename << std::endl;
0330 exit(1);
0331 }
0332
0333 PHG4CylinderCellGeomContainer *seggeo = findNode::getClass<PHG4CylinderCellGeomContainer>(topNode, seggeonodename);
0334 if (!seggeo)
0335 {
0336 std::cout << "PHG4FullProjSpacalCellReco::process_event - Fatal Error - could not locate cell geometry node "
0337 << seggeonodename << std::endl;
0338 exit(1);
0339 }
0340
0341 PHG4CylinderGeomContainer *layergeo = findNode::getClass<PHG4CylinderGeomContainer>(topNode, geonodename);
0342 if (!layergeo)
0343 {
0344 std::cout << "PHG4FullProjSpacalCellReco::process_event - Fatal Error - Could not locate sim geometry node "
0345 << geonodename << std::endl;
0346 exit(1);
0347 }
0348
0349 PHG4HitContainer::ConstIterator hiter;
0350 PHG4HitContainer::ConstRange hit_begin_end = g4hit->getHits();
0351
0352 PHG4CylinderCellGeom *geo_raw = seggeo->GetFirstLayerCellGeom();
0353 PHG4CylinderCellGeom_Spacalv1 *geo = dynamic_cast<PHG4CylinderCellGeom_Spacalv1 *>(geo_raw);
0354 assert(geo);
0355 const PHG4CylinderGeom *layergeom_raw = layergeo->GetFirstLayerGeom();
0356 assert(layergeom_raw);
0357
0358 const PHG4CylinderGeom_Spacalv3 *layergeom =
0359 dynamic_cast<const PHG4CylinderGeom_Spacalv3 *>(layergeom_raw);
0360 assert(layergeom);
0361
0362 for (hiter = hit_begin_end.first; hiter != hit_begin_end.second; ++hiter)
0363 {
0364
0365 if (hiter->second->get_t(0) > tmax)
0366 {
0367 continue;
0368 }
0369 if (hiter->second->get_t(1) < tmin)
0370 {
0371 continue;
0372 }
0373 if (hiter->second->get_t(1) - hiter->second->get_t(0) > m_DeltaT)
0374 {
0375 continue;
0376 }
0377
0378 sum_energy_g4hit += hiter->second->get_edep();
0379
0380
0381 int scint_id = hiter->second->get_scint_id();
0382
0383
0384 PHG4CylinderGeom_Spacalv3::scint_id_coder decoder(scint_id);
0385
0386
0387 std::pair<int, int> tower_z_phi_ID = layergeom->get_tower_z_phi_ID(decoder.tower_ID, decoder.sector_ID);
0388 const int &tower_ID_z = tower_z_phi_ID.first;
0389 const int &tower_ID_phi = tower_z_phi_ID.second;
0390
0391 PHG4CylinderGeom_Spacalv3::tower_map_t::const_iterator it_tower =
0392 layergeom->get_sector_tower_map().find(decoder.tower_ID);
0393 assert(it_tower != layergeom->get_sector_tower_map().end());
0394
0395 unsigned int key = static_cast<unsigned int>(scint_id);
0396 PHG4Cell *cell = nullptr;
0397 std::map<unsigned int, PHG4Cell *>::iterator it = celllist.find(key);
0398 if (it != celllist.end())
0399 {
0400 cell = it->second;
0401 }
0402 else
0403 {
0404
0405 int etabin = -1;
0406 try
0407 {
0408 etabin = geo->get_etabin_block(tower_ID_z);
0409 }
0410 catch (std::exception &e)
0411 {
0412 std::cout << "Print cell geometry:" << std::endl;
0413 geo->identify();
0414 std::cout << "Print scint_id_coder:" << std::endl;
0415 decoder.identify();
0416 std::cout << "Print the hit:" << std::endl;
0417 hiter->second->print();
0418 std::cout << "PHG4FullProjSpacalCellReco::process_event::"
0419 << Name() << " - Fatal Error - " << e.what() << std::endl;
0420 exit(1);
0421 }
0422
0423 const int sub_tower_ID_x = it_tower->second.get_sub_tower_ID_x(decoder.fiber_ID);
0424 const int sub_tower_ID_y = it_tower->second.get_sub_tower_ID_y(decoder.fiber_ID);
0425 unsigned short fiber_ID = decoder.fiber_ID;
0426 unsigned short etabinshort = etabin * layergeom->get_n_subtower_eta() + sub_tower_ID_y;
0427 unsigned short phibin = tower_ID_phi * layergeom->get_n_subtower_phi() + sub_tower_ID_x;
0428 PHG4CellDefs::keytype cellkey = PHG4CellDefs::SpacalBinning::genkey(etabinshort, phibin, fiber_ID);
0429 cell = new PHG4Cellv1(cellkey);
0430 celllist[key] = cell;
0431 }
0432
0433 double light_yield = hiter->second->get_light_yield();
0434
0435
0436
0437 if (light_collection_model.use_fiber_model())
0438 {
0439 const double z = 0.5 * (hiter->second->get_local_z(0) + hiter->second->get_local_z(1));
0440 assert(not std::isnan(z));
0441
0442 light_yield *= light_collection_model.get_fiber_transmission(z);
0443 }
0444
0445
0446 if (light_collection_model.use_fiber_model())
0447 {
0448 const double x = it_tower->second.get_position_fraction_x_in_sub_tower(decoder.fiber_ID);
0449 const double y = it_tower->second.get_position_fraction_y_in_sub_tower(decoder.fiber_ID);
0450
0451 light_yield *= light_collection_model.get_light_guide_efficiency(x, y);
0452 }
0453 cell->add_edep(hiter->first, hiter->second->get_edep());
0454 cell->add_edep(hiter->second->get_edep());
0455 cell->add_light_yield(light_yield);
0456 cell->add_shower_edep(hiter->second->get_shower_id(), hiter->second->get_edep());
0457
0458 }
0459 int numcells = 0;
0460 for (std::map<unsigned int, PHG4Cell *>::const_iterator mapiter =
0461 celllist.begin();
0462 mapiter != celllist.end(); ++mapiter)
0463 {
0464 cells->AddCell(mapiter->second);
0465 numcells++;
0466 if (Verbosity() > 1)
0467 {
0468 std::cout << "PHG4FullProjSpacalCellReco::process_event::" << Name()
0469 << " - "
0470 << "Adding cell in bin eta "
0471 << PHG4CellDefs::SpacalBinning::get_etabin(mapiter->second->get_cellid())
0472 << " phi "
0473 << PHG4CellDefs::SpacalBinning::get_phibin(mapiter->second->get_cellid())
0474 << " fiber "
0475 << PHG4CellDefs::SpacalBinning::get_fiberid(mapiter->second->get_cellid())
0476 << ", energy dep: "
0477 << mapiter->second->get_edep() << ", light yield: "
0478 << mapiter->second->get_light_yield() << std::endl;
0479 }
0480 }
0481 celllist.clear();
0482 if (Verbosity() > 0)
0483 {
0484 std::cout << "PHG4FullProjSpacalCellReco::process_event::" << Name()
0485 << " - "
0486 << " found " << numcells
0487 << " fibers with energy deposition" << std::endl;
0488 }
0489
0490 if (chkenergyconservation || Verbosity() > 4)
0491 {
0492 CheckEnergy(topNode);
0493 }
0494 return Fun4AllReturnCodes::EVENT_OK;
0495 }
0496
0497 int PHG4FullProjSpacalCellReco::CheckEnergy(PHCompositeNode *topNode)
0498 {
0499 PHG4CellContainer *cells = findNode::getClass<PHG4CellContainer>(topNode, cellnodename);
0500 double sum_energy_cells = 0.;
0501 PHG4CellContainer::ConstRange cell_begin_end = cells->getCells();
0502 PHG4CellContainer::ConstIterator citer;
0503 for (citer = cell_begin_end.first; citer != cell_begin_end.second; ++citer)
0504 {
0505 sum_energy_cells += citer->second->get_edep();
0506 }
0507
0508 if (fabs(sum_energy_cells - sum_energy_g4hit) / sum_energy_g4hit > 1e-6)
0509 {
0510 std::cout
0511 << "PHG4FullProjSpacalCellReco::CheckEnergy - energy mismatch between cells: "
0512 << sum_energy_cells << " and hits: " << sum_energy_g4hit
0513 << " diff sum(cells) - sum(hits): "
0514 << sum_energy_cells - sum_energy_g4hit << std::endl;
0515 return -1;
0516 }
0517 else
0518 {
0519 if (Verbosity() > 0)
0520 {
0521 std::cout << "PHG4FullProjSpacalCellReco::CheckEnergy::" << Name()
0522 << " - total energy for this event: " << sum_energy_g4hit
0523 << " GeV. Passed CheckEnergy" << std::endl;
0524 }
0525 }
0526 return 0;
0527 }
0528
0529 void PHG4FullProjSpacalCellReco::SetDefaultParameters()
0530 {
0531 set_default_double_param("tmax", 60.0);
0532 set_default_double_param("tmin", -20.0);
0533 set_default_double_param("delta_t", 100.0);
0534 return;
0535 }
0536
0537 void PHG4FullProjSpacalCellReco::set_timing_window(const double tmi, const double tma)
0538 {
0539 set_double_param("tmin", tmi);
0540 set_double_param("tmax", tma);
0541 }