Back to home page

sPhenix code displayed by LXR

 
 

    


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 * /*topNode*/)
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   // Looking for the DST node
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   // a special implimentation of PHG4CylinderGeom is required here.
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   //      int layer = layergeom->get_layer();
0151 
0152   // create geo object and fill with variables common to all binning methods
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   // construct a map to convert tower_ID into the older eta bins.
0160 
0161   const PHG4CylinderGeom_Spacalv3::tower_map_t &tower_map = layergeom->get_sector_tower_map();
0162   const PHG4CylinderGeom_Spacalv3::sector_map_t &sector_map = layergeom->get_sector_map();
0163   const int nphibin = layergeom->get_azimuthal_n_sec()       // sector
0164                       * layergeom->get_max_phi_bin_in_sec()  // blocks per sector
0165                       * layergeom->get_n_subtower_phi();     // subtower per block
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     // inspect index in sector 0
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       // assign phi min according phi bin 0
0186 // NOLINTNEXTLINE(bugprone-integer-division)
0187       phi_min = M_PI_2 - deltaphi * (layergeom->get_max_phi_bin_in_sec() * layergeom->get_n_subtower_phi() / 2)  // shift of first tower in sector
0188                 + sector_map.begin()->second;
0189     }
0190 
0191     if (tower_ID_phi == layergeom->get_max_phi_bin_in_sec() / 2)
0192     {
0193       // assign eta min according phi bin 0
0194       map_z_tower_z_ID[tower.centralZ] = tower_ID_z;
0195     }
0196     // ...
0197   }  //       const auto &tower_pair: tower_map
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   // build eta bin maps
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     // inspect index in sector 0
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       // half z-range
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       // half eta-range
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         // do not overlap to the next bin.
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   }  //       const auto &tower_pair: tower_map
0277 
0278   // add geo object filled by different binning methods
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   // save this to the run wise tree to store on DST
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   // save this to the parNode for use
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   // a special implimentation of PHG4CylinderGeom is required here.
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     // checking ADC timing integration window cut
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     // hit loop
0381     int scint_id = hiter->second->get_scint_id();
0382 
0383     // decode scint_id
0384     PHG4CylinderGeom_Spacalv3::scint_id_coder decoder(scint_id);
0385 
0386     // convert to z_ID, phi_ID
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       // convert tower_ID_z to to eta bin number
0405       int etabin = -1;
0406       try
0407       {
0408         etabin = geo->get_etabin_block(tower_ID_z);  // block eta bin
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     // light yield correction from fiber attenuation:
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     // light yield correction from light guide collection efficiency:
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   }  // end loop over g4hits
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   // the fractional eloss for particles traversing eta bins leads to minute rounding errors
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);  // collision has a timing spread around the triggered event. Accepting negative time too.
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 }