Back to home page

sPhenix code displayed by LXR

 
 

    


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 * /*topNode*/)
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   // Looking for the DST node
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     // create geo object and fill with variables common to all binning methods
0165     PHG4BlockCellGeom *layerseggeo = new PHG4BlockCellGeom();
0166     layerseggeo->set_layer(layergeom->get_layer());
0167     // layerseggeo->set_radius(layergeom->get_radius());
0168     // layerseggeo->set_thickness(layergeom->get_thickness());
0169 
0170     if (binning[layer] == PHG4CellDefs::etaphibinning)
0171     {
0172       // calculate eta at radius+ thickness (outer radius)
0173       // length via eta coverage is calculated using the outer radius
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))  // etahi is a tiny bit larger due to numerical uncertainties
0191         {
0192           cout << "etahi: " << etahi << ", etamax: " << etamax << endl;
0193         }
0194         etahi += etastepsize;
0195       }
0196 
0197       double xmin = -layergeom->get_width() / 2.;
0198       // double xmax = -xmin;
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       // double xlow = xmin;
0211       // double xhi = xlow + xstepsize;
0212 
0213       // for (int i = 0; i < xbins; i++)
0214       // {
0215       //   if (xhi > xmax)
0216       //   {
0217       //     cout << "xhi: " << xhi << ", xmax: " << xmax << endl;
0218       //   }
0219       //   xlow = xhi;
0220       //   xhi +=  xstepsize;
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     // add geo object filled by different binning methods
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   //   cout << "number of layers: " << g4hit->num_layers() << endl;
0273   //   cout << "number of hits: " << g4hit->size() << endl;
0274   //   for (layer = layer_begin_end.first; layer != layer_begin_end.second; layer++)
0275   //     {
0276   //       cout << "layer number: " << *layer << endl;
0277   //     }
0278 
0279   for (layer = layer_begin_end.first; layer != layer_begin_end.second; layer++)
0280   {
0281     // only handle layers/detector ids which have parameters set
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     // ------- eta/x binning ------------------------------------------------------------------------
0299     if (binning[*layer] == PHG4CellDefs::etaphibinning)
0300     {
0301       for (hiter = hit_begin_end.first; hiter != hit_begin_end.second; hiter++)
0302       {
0303         // checking ADC timing integration window cut
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         // check bin range
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         // Determine all fired cells
0349 
0350         double ax = (etax[0]).second;  // x
0351         double ay = (etax[0]).first;   // eta
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         // if entry and exit hit are the same (seems to happen rarely), trklen = 0
0370         // which leads to a 0/0 and an NaN in edep later on
0371         // this code does for particles in the same cell a trklen/trklen (vdedx[ii]/trklen)
0372         // so setting this to any non zero number will do just fine
0373         // I just pick -1 here to flag those strange hits in case I want t oanalyze them
0374         // later on
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)  // single cell fired
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               // cout << "##### line: " << ax << " " << ay << " " << bx << " " << by << endl;
0404               // cout << "####### cell: " << cx << " " << cy << " " << dx << " " << dy << endl;
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++)  // loop over all fired cells
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           // just a sanity check - we don't want to mess up by having Nan's or Infs in our energy deposition
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       }  // end loop over g4hits
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   // the fractional eloss for particles traversing eta bins leads to minute rounding errors
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 }