Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2025-08-05 08:17:45 

0001 #include "PHG4CylinderCellReco.h"
0002 #include "PHG4CylinderCellGeom.h"
0003 #include "PHG4CylinderCellGeomContainer.h"
0004 #include "PHG4CylinderGeom.h"
0005 #include "PHG4CylinderGeomContainer.h"
0006 
0007 #include "PHG4Cell.h"  // for PHG4Cell, PHG...
0008 #include "PHG4CellContainer.h"
0009 #include "PHG4CellDefs.h"
0010 #include "PHG4Cellv1.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 <cstring>  // for memset
0033 #include <iostream>
0034 #include <iterator>  // for reverse_iterator
0035 #include <sstream>
0036 #include <vector>  // for vector
0037 
0038 PHG4CylinderCellReco::PHG4CylinderCellReco(const std::string &name)
0039   : SubsysReco(name)
0040   , PHParameterContainerInterface(name)
0041 {
0042   SetDefaultParameters();
0043 }
0044 
0045 int PHG4CylinderCellReco::ResetEvent(PHCompositeNode * /*topNode*/)
0046 {
0047   sum_energy_before_cuts = 0.;
0048   sum_energy_g4hit = 0.;
0049   return Fun4AllReturnCodes::EVENT_OK;
0050 }
0051 
0052 int PHG4CylinderCellReco::InitRun(PHCompositeNode *topNode)
0053 {
0054   PHNodeIterator nodeiter(topNode);
0055 
0056   // Looking for the DST node
0057   PHCompositeNode *dstNode;
0058   dstNode = dynamic_cast<PHCompositeNode *>(nodeiter.findFirst("PHCompositeNode", "DST"));
0059   if (!dstNode)
0060   {
0061     std::cout << PHWHERE << "DST Node missing, doing nothing." << std::endl;
0062     exit(1);
0063   }
0064   PHCompositeNode *runNode;
0065   runNode = dynamic_cast<PHCompositeNode *>(nodeiter.findFirst("PHCompositeNode", "RUN"));
0066   if (!runNode)
0067   {
0068     std::cout << Name() << "DST Node missing, doing nothing." << std::endl;
0069     exit(1);
0070   }
0071   PHNodeIterator runiter(runNode);
0072   PHCompositeNode *RunDetNode = dynamic_cast<PHCompositeNode *>(runiter.findFirst("PHCompositeNode", detector));
0073   if (!RunDetNode)
0074   {
0075     RunDetNode = new PHCompositeNode(detector);
0076     runNode->addNode(RunDetNode);
0077   }
0078 
0079   hitnodename = "G4HIT_" + detector;
0080   PHG4HitContainer *g4hit = findNode::getClass<PHG4HitContainer>(topNode, hitnodename);
0081   if (!g4hit)
0082   {
0083     std::cout << "Could not locate g4 hit node " << hitnodename << std::endl;
0084     exit(1);
0085   }
0086   cellnodename = "G4CELL_" + outdetector;
0087   PHG4CellContainer *cells = findNode::getClass<PHG4CellContainer>(topNode, cellnodename);
0088   if (!cells)
0089   {
0090     PHNodeIterator dstiter(dstNode);
0091     PHCompositeNode *DetNode = dynamic_cast<PHCompositeNode *>(dstiter.findFirst("PHCompositeNode", detector));
0092     if (!DetNode)
0093     {
0094       DetNode = new PHCompositeNode(detector);
0095       dstNode->addNode(DetNode);
0096     }
0097     cells = new PHG4CellContainer();
0098     PHIODataNode<PHObject> *newNode = new PHIODataNode<PHObject>(cells, cellnodename, "PHObject");
0099     DetNode->addNode(newNode);
0100   }
0101 
0102   geonodename = "CYLINDERGEOM_" + detector;
0103   PHG4CylinderGeomContainer *geo = findNode::getClass<PHG4CylinderGeomContainer>(topNode, geonodename);
0104   if (!geo)
0105   {
0106     std::cout << "Could not locate geometry node " << geonodename << std::endl;
0107     exit(1);
0108   }
0109   if (Verbosity() > 0)
0110   {
0111     geo->identify();
0112   }
0113   seggeonodename = "CYLINDERCELLGEOM_" + outdetector;
0114   PHG4CylinderCellGeomContainer *seggeo = findNode::getClass<PHG4CylinderCellGeomContainer>(topNode, seggeonodename);
0115   if (!seggeo)
0116   {
0117     seggeo = new PHG4CylinderCellGeomContainer();
0118     PHIODataNode<PHObject> *newNode = new PHIODataNode<PHObject>(seggeo, seggeonodename, "PHObject");
0119     runNode->addNode(newNode);
0120   }
0121 
0122   GetParamsContainerModify()->set_name(detector);
0123 
0124   UpdateParametersWithMacro();
0125 
0126   std::map<int, PHG4CylinderGeom *>::const_iterator miter;
0127   std::pair<std::map<int, PHG4CylinderGeom *>::const_iterator, std::map<int, PHG4CylinderGeom *>::const_iterator> begin_end = geo->get_begin_end();
0128   std::map<int, std::pair<double, double> >::iterator sizeiter;
0129   for (miter = begin_end.first; miter != begin_end.second; ++miter)
0130   {
0131     PHG4CylinderGeom *layergeom = miter->second;
0132     int layer = layergeom->get_layer();
0133     if (!ExistDetid(layer))
0134     {
0135       std::cout << Name() << ": No parameters for detid/layer " << layer
0136                 << ", hits from this detid/layer will not be accumulated into cells" << std::endl;
0137       continue;
0138     }
0139     implemented_detid.insert(layer);
0140     set_size(layer, get_double_param(layer, "size_long"), get_double_param(layer, "size_perp"));
0141     tmin_max.insert(std::make_pair(layer, std::make_pair(get_double_param(layer, "tmin"), get_double_param(layer, "tmax"))));
0142     m_DeltaTMap.insert(std::make_pair(layer, get_double_param(layer, "delta_t")));
0143     double circumference = layergeom->get_radius() * 2 * M_PI;
0144     double length_in_z = layergeom->get_zmax() - layergeom->get_zmin();
0145     sizeiter = cell_size.find(layer);
0146     if (sizeiter == cell_size.end())
0147     {
0148       std::cout << "no cell sizes for layer " << layer << std::endl;
0149       exit(1);
0150     }
0151     // create geo object and fill with variables common to all binning methods
0152     PHG4CylinderCellGeom *layerseggeo = new PHG4CylinderCellGeom();
0153     layerseggeo->set_layer(layergeom->get_layer());
0154     layerseggeo->set_radius(layergeom->get_radius());
0155     layerseggeo->set_thickness(layergeom->get_thickness());
0156     if (binning[layer] == PHG4CellDefs::etaphibinning)
0157     {
0158       // calculate eta at radius+ thickness (outer radius)
0159       // length via eta coverage is calculated using the outer radius
0160       double etamin = PHG4Utils::get_eta(layergeom->get_radius() + layergeom->get_thickness(), layergeom->get_zmin());
0161       double etamax = PHG4Utils::get_eta(layergeom->get_radius() + layergeom->get_thickness(), layergeom->get_zmax());
0162       zmin_max[layer] = std::make_pair(etamin, etamax);
0163       double etastepsize = (sizeiter->second).first;
0164       double d_etabins;
0165       // if the eta cell size is larger than the eta range, make one bin
0166       if (etastepsize > etamax - etamin)
0167       {
0168         d_etabins = 1;
0169       }
0170       else
0171       {
0172         // it is unlikely that the eta range is a multiple of the eta cell size
0173         // then fract is 0, if not - add 1 bin which makes the
0174         // cells a tiny bit smaller but makes them fit
0175         double fract = modf((etamax - etamin) / etastepsize, &d_etabins);
0176         if (fract != 0)
0177         {
0178           d_etabins++;
0179         }
0180       }
0181       etastepsize = (etamax - etamin) / d_etabins;
0182       (sizeiter->second).first = etastepsize;
0183       int etabins = d_etabins;
0184       double etalow = etamin;
0185       double etahi = etalow + etastepsize;
0186       for (int i = 0; i < etabins; i++)
0187       {
0188         if (etahi > (etamax + 1.e-6))  // etahi is a tiny bit larger due to numerical uncertainties
0189         {
0190           std::cout << "etahi: " << etahi << ", etamax: " << etamax << std::endl;
0191         }
0192         etahi += etastepsize;
0193       }
0194 
0195       double phimin = -M_PI;
0196       double phimax = M_PI;
0197       double phistepsize = (sizeiter->second).second;
0198       double d_phibins;
0199       if (phistepsize >= phimax - phimin)
0200       {
0201         d_phibins = 1;
0202       }
0203       else
0204       {
0205         // it is unlikely that the phi range is a multiple of the phi cell size
0206         // then fract is 0, if not - add 1 bin which makes the
0207         // cells a tiny bit smaller but makes them fit
0208         double fract = modf((phimax - phimin) / phistepsize, &d_phibins);
0209         if (fract != 0)
0210         {
0211           d_phibins++;
0212         }
0213       }
0214       phistepsize = (phimax - phimin) / d_phibins;
0215       (sizeiter->second).second = phistepsize;
0216       int phibins = d_phibins;
0217       double philow = phimin;
0218       double phihi = philow + phistepsize;
0219       for (int i = 0; i < phibins; i++)
0220       {
0221         if (phihi > phimax)
0222         {
0223           std::cout << "phihi: " << phihi << ", phimax: " << phimax << std::endl;
0224         }
0225         phihi += phistepsize;
0226       }
0227       std::pair<int, int> phi_z_bin = std::make_pair(phibins, etabins);
0228       n_phi_z_bins[layer] = phi_z_bin;
0229       layerseggeo->set_binning(PHG4CellDefs::etaphibinning);
0230       layerseggeo->set_etabins(etabins);
0231       layerseggeo->set_etamin(etamin);
0232       layerseggeo->set_etastep(etastepsize);
0233       layerseggeo->set_phimin(phimin);
0234       layerseggeo->set_phibins(phibins);
0235       layerseggeo->set_phistep(phistepsize);
0236       phistep[layer] = phistepsize;
0237       etastep[layer] = etastepsize;
0238     }
0239     else if (binning[layer] == PHG4CellDefs::sizebinning)
0240     {
0241       zmin_max[layer] = std::make_pair(layergeom->get_zmin(), layergeom->get_zmax());
0242       double size_z = (sizeiter->second).second;
0243       double size_r = (sizeiter->second).first;
0244       double bins_r;
0245       // if the size is larger than circumference, make it one bin
0246       if (size_r >= circumference)
0247       {
0248         bins_r = 1;
0249       }
0250       else
0251       {
0252         // unlikely but if the circumference is a multiple of the cell size
0253         // use result of division, if not - add 1 bin which makes the
0254         // cells a tiny bit smaller but makes them fit
0255         double fract = modf(circumference / size_r, &bins_r);
0256         if (fract != 0)
0257         {
0258           bins_r++;
0259         }
0260       }
0261       nbins[0] = bins_r;
0262       size_r = circumference / bins_r;
0263       (sizeiter->second).first = size_r;
0264       double phistepsize = 2 * M_PI / bins_r;
0265       double phimin = -M_PI;
0266       double phimax = phimin + phistepsize;
0267       phistep[layer] = phistepsize;
0268       for (int i = 0; i < nbins[0]; i++)
0269       {
0270         if (phimax > (M_PI + 1e-9))
0271         {
0272           std::cout << "phimax: " << phimax << ", M_PI: " << M_PI
0273                     << "phimax-M_PI: " << phimax - M_PI << std::endl;
0274         }
0275         phimax += phistepsize;
0276       }
0277       // if the size is larger than length, make it one bin
0278       if (size_z >= length_in_z)
0279       {
0280         bins_r = 1;
0281       }
0282       else
0283       {
0284         // unlikely but if the length is a multiple of the cell size
0285         // use result of division, if not - add 1 bin which makes the
0286         // cells a tiny bit smaller but makes them fit
0287         double fract = modf(length_in_z / size_z, &bins_r);
0288         if (fract != 0)
0289         {
0290           bins_r++;
0291         }
0292       }
0293       nbins[1] = bins_r;
0294       std::pair<int, int> phi_z_bin = std::make_pair(nbins[0], nbins[1]);
0295       n_phi_z_bins[layer] = phi_z_bin;
0296       // update our map with the new sizes
0297       size_z = length_in_z / bins_r;
0298       (sizeiter->second).second = size_z;
0299       double zlow = layergeom->get_zmin();
0300       double zhigh = zlow + size_z;
0301       ;
0302       for (int i = 0; i < nbins[1]; i++)
0303       {
0304         if (zhigh > (layergeom->get_zmax() + 1e-9))
0305         {
0306           std::cout << "zhigh: " << zhigh << ", zmax "
0307                     << layergeom->get_zmax()
0308                     << ", zhigh-zmax: " << zhigh - layergeom->get_zmax()
0309                     << std::endl;
0310         }
0311         zhigh += size_z;
0312       }
0313       layerseggeo->set_binning(PHG4CellDefs::sizebinning);
0314       layerseggeo->set_zbins(nbins[1]);
0315       layerseggeo->set_zmin(layergeom->get_zmin());
0316       layerseggeo->set_zstep(size_z);
0317       layerseggeo->set_phibins(nbins[0]);
0318       layerseggeo->set_phistep(phistepsize);
0319     }
0320     // add geo object filled by different binning methods
0321     seggeo->AddLayerCellGeom(layerseggeo);
0322     if (Verbosity() > 1)
0323     {
0324       layerseggeo->identify();
0325     }
0326   }
0327 
0328   // print out settings
0329   if (Verbosity() > 0)
0330   {
0331     std::cout << "===================== PHG4CylinderCellReco::InitRun() =====================" << std::endl;
0332     std::cout << " " << outdetector << " Segmentation Description: " << std::endl;
0333     for (std::map<int, int>::iterator iter = binning.begin();
0334          iter != binning.end(); ++iter)
0335     {
0336       int layer = iter->first;
0337 
0338       if (binning[layer] == PHG4CellDefs::etaphibinning)
0339       {
0340         // phi & eta bin is usually used to make projective towers
0341         // so just print the first layer
0342         std::cout << " Layer #" << binning.begin()->first << "-" << binning.rbegin()->first << std::endl;
0343         std::cout << "   Nbins (phi,eta): (" << n_phi_z_bins[layer].first << ", " << n_phi_z_bins[layer].second << ")" << std::endl;
0344         std::cout << "   Cell Size (phi,eta): (" << cell_size[layer].first << " rad, " << cell_size[layer].second << " units)" << std::endl;
0345         break;
0346       }
0347       else if (binning[layer] == PHG4CellDefs::sizebinning)
0348       {
0349         std::cout << " Layer #" << layer << std::endl;
0350         std::cout << "   Nbins (phi,z): (" << n_phi_z_bins[layer].first << ", " << n_phi_z_bins[layer].second << ")" << std::endl;
0351         std::cout << "   Cell Size (phi,z): (" << cell_size[layer].first << " cm, " << cell_size[layer].second << " cm)" << std::endl;
0352       }
0353     }
0354     std::cout << "===========================================================================" << std::endl;
0355   }
0356   std::string nodename = "G4CELLPARAM_" + GetParamsContainer()->Name();
0357   SaveToNodeTree(RunDetNode, nodename);
0358   return Fun4AllReturnCodes::EVENT_OK;
0359 }
0360 
0361 int PHG4CylinderCellReco::process_event(PHCompositeNode *topNode)
0362 {
0363   PHG4HitContainer *g4hit = findNode::getClass<PHG4HitContainer>(topNode, hitnodename);
0364   if (!g4hit)
0365   {
0366     std::cout << "Could not locate g4 hit node " << hitnodename << std::endl;
0367     exit(1);
0368   }
0369   PHG4CellContainer *cells = findNode::getClass<PHG4CellContainer>(topNode, cellnodename);
0370   if (!cells)
0371   {
0372     std::cout << "could not locate cell node " << cellnodename << std::endl;
0373     exit(1);
0374   }
0375 
0376   PHG4CylinderCellGeomContainer *seggeo = findNode::getClass<PHG4CylinderCellGeomContainer>(topNode, seggeonodename);
0377   if (!seggeo)
0378   {
0379     std::cout << "could not locate geo node " << seggeonodename << std::endl;
0380     exit(1);
0381   }
0382 
0383   std::map<int, std::pair<double, double> >::iterator sizeiter;
0384   PHG4HitContainer::LayerIter layer;
0385   std::pair<PHG4HitContainer::LayerIter, PHG4HitContainer::LayerIter> layer_begin_end = g4hit->getLayers();
0386   //   std::cout << "number of layers: " << g4hit->num_layers() << std::endl;
0387   //   std::cout << "number of hits: " << g4hit->size() << std::endl;
0388   //   for (layer = layer_begin_end.first; layer != layer_begin_end.second; layer++)
0389   //     {
0390   //       std::cout << "layer number: " << *layer << std::endl;
0391   //     }
0392   for (layer = layer_begin_end.first; layer != layer_begin_end.second; layer++)
0393   {
0394     // only handle layers/detector ids which have parameters set
0395     if (implemented_detid.find(*layer) == implemented_detid.end())
0396     {
0397       continue;
0398     }
0399     PHG4HitContainer::ConstIterator hiter;
0400     PHG4HitContainer::ConstRange hit_begin_end = g4hit->getHits(*layer);
0401     PHG4CylinderCellGeom *geo = seggeo->GetLayerCellGeom(*layer);
0402     int nphibins = n_phi_z_bins[*layer].first;
0403     int nzbins = n_phi_z_bins[*layer].second;
0404 
0405     // ------- eta/phi binning ------------------------------------------------------------------------
0406     if (binning[*layer] == PHG4CellDefs::etaphibinning)
0407     {
0408       for (hiter = hit_begin_end.first; hiter != hit_begin_end.second; hiter++)
0409       {
0410         sum_energy_before_cuts += hiter->second->get_edep();
0411         // checking ADC timing integration window cut
0412         if (hiter->second->get_t(0) > tmin_max[*layer].second)
0413         {
0414           continue;
0415         }
0416         if (hiter->second->get_t(1) < tmin_max[*layer].first)
0417         {
0418           continue;
0419         }
0420         if (hiter->second->get_t(1) - hiter->second->get_t(0) > m_DeltaTMap[*layer])
0421         {
0422           continue;
0423         }
0424 
0425         std::pair<double, double> etaphi[2];
0426         double phibin[2];
0427         double etabin[2];
0428         for (int i = 0; i < 2; i++)
0429         {
0430           etaphi[i] = PHG4Utils::get_etaphi(hiter->second->get_x(i), hiter->second->get_y(i), hiter->second->get_z(i));
0431           etabin[i] = geo->get_etabin(etaphi[i].first);
0432           phibin[i] = geo->get_phibin(etaphi[i].second);
0433         }
0434         // check bin range
0435         if (phibin[0] < 0 || phibin[0] >= nphibins || phibin[1] < 0 || phibin[1] >= nphibins)
0436         {
0437           continue;
0438         }
0439         if (etabin[0] < 0 || etabin[0] >= nzbins || etabin[1] < 0 || etabin[1] >= nzbins)
0440         {
0441           continue;
0442         }
0443 
0444         if (etabin[0] < 0)
0445         {
0446           if (Verbosity() > 0)
0447           {
0448             hiter->second->identify();
0449           }
0450           continue;
0451         }
0452         sum_energy_g4hit += hiter->second->get_edep();
0453         int intphibin = phibin[0];
0454         int intetabin = etabin[0];
0455         int intphibinout = phibin[1];
0456         int intetabinout = etabin[1];
0457 
0458         // Determine all fired cells
0459 
0460         double ax = (etaphi[0]).second;  // phi
0461         double ay = (etaphi[0]).first;   // eta
0462         double bx = (etaphi[1]).second;
0463         double by = (etaphi[1]).first;
0464         if (intphibin > intphibinout)
0465         {
0466           int tmp = intphibin;
0467           intphibin = intphibinout;
0468           intphibinout = tmp;
0469         }
0470         if (intetabin > intetabinout)
0471         {
0472           int tmp = intetabin;
0473           intetabin = intetabinout;
0474           intetabinout = tmp;
0475         }
0476 
0477         double trklen = sqrt((ax - bx) * (ax - bx) + (ay - by) * (ay - by));
0478         // if entry and exit hit are the same (seems to happen rarely), trklen = 0
0479         // which leads to a 0/0 and an NaN in edep later on
0480         // this code does for particles in the same cell a trklen/trklen (vdedx[ii]/trklen)
0481         // so setting this to any non zero number will do just fine
0482         // I just pick -1 here to flag those strange hits in case I want t oanalyze them
0483         // later on
0484         if (trklen == 0)
0485         {
0486           trklen = -1.;
0487         }
0488         std::vector<int> vphi;
0489         std::vector<int> veta;
0490         std::vector<double> vdedx;
0491 
0492         if (intphibin == intphibinout && intetabin == intetabinout)  // single cell fired
0493         {
0494           if (Verbosity() > 0)
0495           {
0496             std::cout << "SINGLE CELL FIRED: " << intphibin << " " << intetabin << std::endl;
0497           }
0498           vphi.push_back(intphibin);
0499           veta.push_back(intetabin);
0500           vdedx.push_back(trklen);
0501         }
0502         else
0503         {
0504           double phistep_half = geo->get_phistep() / 2.;
0505           double etastep_half = geo->get_etastep() / 2.;
0506           for (int ibp = intphibin; ibp <= intphibinout; ibp++)
0507           {
0508             double cx = geo->get_phicenter(ibp) - phistep_half;
0509             double dx = geo->get_phicenter(ibp) + phistep_half;
0510             for (int ibz = intetabin; ibz <= intetabinout; ibz++)
0511             {
0512               double cy = geo->get_etacenter(ibz) - etastep_half;
0513               double dy = geo->get_etacenter(ibz) + etastep_half;
0514 
0515               // std::cout << "##### line: " << ax << " " << ay << " " << bx << " " << by << std::endl;
0516               // std::cout << "####### cell: " << cx << " " << cy << " " << dx << " " << dy << std::endl;
0517               std::pair<bool, double> intersect = PHG4Utils::line_and_rectangle_intersect(ax, ay, bx, by, cx, cy, dx, dy);
0518               if (intersect.first)
0519               {
0520                 if (Verbosity() > 0)
0521                 {
0522                   std::cout << "CELL FIRED: " << ibp << " " << ibz << " " << intersect.second << std::endl;
0523                 }
0524                 vphi.push_back(ibp);
0525                 veta.push_back(ibz);
0526                 vdedx.push_back(intersect.second);
0527               }
0528             }
0529           }
0530         }
0531         if (Verbosity() > 0)
0532         {
0533           std::cout << "NUMBER OF FIRED CELLS = " << vphi.size() << std::endl;
0534         }
0535 
0536         double tmpsum = 0.;
0537         for (unsigned int ii = 0; ii < vphi.size(); ii++)
0538         {
0539           tmpsum += vdedx[ii];
0540           vdedx[ii] = vdedx[ii] / trklen;
0541           if (Verbosity() > 0)
0542           {
0543             std::cout << "  CELL " << ii << "  dE/dX = " << vdedx[ii] << std::endl;
0544           }
0545         }
0546         if (Verbosity() > 0)
0547         {
0548           std::cout << "    TOTAL TRACK LENGTH = " << tmpsum << " " << trklen << std::endl;
0549         }
0550 
0551         for (unsigned int i1 = 0; i1 < vphi.size(); i1++)  // loop over all fired cells
0552         {
0553           int iphibin = vphi[i1];
0554           int ietabin = veta[i1];
0555 
0556           // This is the key for cellptmap
0557           // It is constructed using the phi and z (or eta) bin index values
0558           // It will be unique for a given phi and z (or eta) bin combination
0559           unsigned long long tmp = iphibin;
0560           unsigned long long key = tmp << 32U;  // clang-tidy: U for unsigned int
0561           key += ietabin;
0562           if (Verbosity() > 1)
0563           {
0564             std::cout << " iphibin " << iphibin << " ietabin " << ietabin << " key 0x" << std::hex << key << std::dec << std::endl;
0565           }
0566           PHG4Cell *cell = nullptr;
0567           it = cellptmap.find(key);
0568           if (it != cellptmap.end())
0569           {
0570             cell = it->second;
0571           }
0572           else
0573           {
0574             PHG4CellDefs::keytype cellkey = PHG4CellDefs::EtaPhiBinning::genkey(*layer, ietabin, iphibin);
0575             cell = new PHG4Cellv1(cellkey);
0576             cellptmap[key] = cell;
0577           }
0578           if (!std::isfinite(hiter->second->get_edep() * vdedx[i1]))
0579           {
0580             std::cout << "hit 0x" << std::hex << hiter->first << std::dec << " not finite, edep: "
0581                       << hiter->second->get_edep() << " weight " << vdedx[i1] << std::endl;
0582           }
0583           cell->add_edep(hiter->first, hiter->second->get_edep() * vdedx[i1]);  // add hit with edep to g4hit list
0584           cell->add_edep(hiter->second->get_edep() * vdedx[i1]);                // add edep to cell
0585           if (hiter->second->has_property(PHG4Hit::prop_light_yield))
0586           {
0587             cell->add_light_yield(hiter->second->get_light_yield() * vdedx[i1]);
0588           }
0589           cell->add_shower_edep(hiter->second->get_shower_id(), hiter->second->get_edep() * vdedx[i1]);
0590           // just a sanity check - we don't want to mess up by having Nan's or Infs in our energy deposition
0591           if (!std::isfinite(hiter->second->get_edep() * vdedx[i1]))
0592           {
0593             std::cout << PHWHERE << " invalid energy dep " << hiter->second->get_edep()
0594                       << " or path length: " << vdedx[i1] << std::endl;
0595           }
0596         }
0597         vphi.clear();
0598         veta.clear();
0599 
0600       }  // end loop over g4hits
0601 
0602       int numcells = 0;
0603 
0604       for (it = cellptmap.begin(); it != cellptmap.end(); ++it)
0605       {
0606         cells->AddCell(it->second);
0607         numcells++;
0608         if (Verbosity() > 1)
0609         {
0610           std::cout << "Adding cell in bin phi: " << PHG4CellDefs::EtaPhiBinning::get_phibin(it->second->get_cellid())
0611                     << " phi: " << geo->get_phicenter(PHG4CellDefs::EtaPhiBinning::get_phibin(it->second->get_cellid())) * 180. / M_PI
0612                     << ", eta bin: " << PHG4CellDefs::EtaPhiBinning::get_etabin(it->second->get_cellid())
0613                     << ", eta: " << geo->get_etacenter(PHG4CellDefs::EtaPhiBinning::get_etabin(it->second->get_cellid()))
0614                     << ", energy dep: " << it->second->get_edep()
0615                     << std::endl;
0616         }
0617       }
0618 
0619       if (Verbosity() > 0)
0620       {
0621         std::cout << Name() << ": found " << numcells << " eta/phi cells with energy deposition" << std::endl;
0622       }
0623     }
0624 
0625     else  // ------ size binning ---------------------------------------------------------------
0626     {
0627       sizeiter = cell_size.find(*layer);
0628       if (sizeiter == cell_size.end())
0629       {
0630         std::cout << "logical screwup!!! no sizes for layer " << *layer << std::endl;
0631         exit(1);
0632       }
0633       double zstepsize = (sizeiter->second).second;
0634       double phistepsize = phistep[*layer];
0635 
0636       for (hiter = hit_begin_end.first; hiter != hit_begin_end.second; hiter++)
0637       {
0638         sum_energy_before_cuts += hiter->second->get_edep();
0639         // checking ADC timing integration window cut
0640         if (hiter->second->get_t(0) > tmin_max[*layer].second)
0641         {
0642           continue;
0643         }
0644         if (hiter->second->get_t(1) < tmin_max[*layer].first)
0645         {
0646           continue;
0647         }
0648         if (hiter->second->get_t(1) - hiter->second->get_t(0) > 100)
0649         {
0650           continue;
0651         }
0652 
0653         double xinout[2];
0654         double yinout[2];
0655         double px[2];
0656         double py[2];
0657         double phi[2];
0658         double z[2];
0659         double phibin[2];
0660         double zbin[2];
0661         if (Verbosity() > 0)
0662         {
0663           std::cout << "--------- new hit in layer # " << *layer << std::endl;
0664         }
0665 
0666         for (int i = 0; i < 2; i++)
0667         {
0668           xinout[i] = hiter->second->get_x(i);
0669           yinout[i] = hiter->second->get_y(i);
0670           px[i] = hiter->second->get_px(i);
0671           py[i] = hiter->second->get_py(i);
0672           phi[i] = std::atan2(hiter->second->get_y(i), hiter->second->get_x(i));
0673           z[i] = hiter->second->get_z(i);
0674           phibin[i] = geo->get_phibin(phi[i]);
0675           zbin[i] = geo->get_zbin(hiter->second->get_z(i));
0676 
0677           if (Verbosity() > 0)
0678           {
0679             std::cout << " " << i << "  phibin: " << phibin[i] << ", phi: " << phi[i] << ", stepsize: " << phistepsize << std::endl;
0680           }
0681           if (Verbosity() > 0)
0682           {
0683             std::cout << " " << i << "  zbin: " << zbin[i] << ", z = " << hiter->second->get_z(i) << ", stepsize: " << zstepsize << " offset: " << zmin_max[*layer].first << std::endl;
0684           }
0685         }
0686         // check bin range
0687         if (phibin[0] < 0 || phibin[0] >= nphibins || phibin[1] < 0 || phibin[1] >= nphibins)
0688         {
0689           continue;
0690         }
0691         if (zbin[0] < 0 || zbin[0] >= nzbins || zbin[1] < 0 || zbin[1] >= nzbins)
0692         {
0693           continue;
0694         }
0695 
0696         if (zbin[0] < 0)
0697         {
0698           hiter->second->identify();
0699           continue;
0700         }
0701         sum_energy_g4hit += hiter->second->get_edep();
0702 
0703         int intphibin = phibin[0];
0704         int intzbin = zbin[0];
0705         int intphibinout = phibin[1];
0706         int intzbinout = zbin[1];
0707 
0708         if (Verbosity() > 0)
0709         {
0710           std::cout << "    phi bin range: " << intphibin << " to " << intphibinout << " phi: " << phi[0] << " to " << phi[1] << std::endl;
0711           std::cout << "    Z bin range: " << intzbin << " to " << intzbinout << " Z: " << z[0] << " to " << z[1] << std::endl;
0712           std::cout << "    phi difference: " << (phi[1] - phi[0]) * 1000. << " milliradians." << std::endl;
0713           std::cout << "    phi difference: " << 2.5 * (phi[1] - phi[0]) * 10000. << " microns." << std::endl;
0714           std::cout << "    path length = " << sqrt((xinout[1] - xinout[0]) * (xinout[1] - xinout[0]) + (yinout[1] - yinout[0]) * (yinout[1] - yinout[0])) << std::endl;
0715           std::cout << "       px = " << px[0] << " " << px[1] << std::endl;
0716           std::cout << "       py = " << py[0] << " " << py[1] << std::endl;
0717           std::cout << "       x = " << xinout[0] << " " << xinout[1] << std::endl;
0718           std::cout << "       y = " << yinout[0] << " " << yinout[1] << std::endl;
0719         }
0720 
0721         // Determine all fired cells
0722 
0723         double ax = phi[0];
0724         double ay = z[0];
0725         double bx = phi[1];
0726         double by = z[1];
0727         if (intphibin > intphibinout)
0728         {
0729           int tmp = intphibin;
0730           intphibin = intphibinout;
0731           intphibinout = tmp;
0732         }
0733         if (intzbin > intzbinout)
0734         {
0735           int tmp = intzbin;
0736           intzbin = intzbinout;
0737           intzbinout = tmp;
0738         }
0739 
0740         double trklen = sqrt((ax - bx) * (ax - bx) + (ay - by) * (ay - by));
0741         // if entry and exit hit are the same (seems to happen rarely), trklen = 0
0742         // which leads to a 0/0 and an NaN in edep later on
0743         // this code does for particles in the same cell a trklen/trklen (vdedx[ii]/trklen)
0744         // so setting this to any non zero number will do just fine
0745         // I just pick -1 here to flag those strange hits in case I want to analyze them
0746         // later on
0747         if (trklen == 0)
0748         {
0749           trklen = -1.;
0750         }
0751         std::vector<int> vphi;
0752         std::vector<int> vz;
0753         std::vector<double> vdedx;
0754 
0755         if (intphibin == intphibinout && intzbin == intzbinout)  // single cell fired
0756         {
0757           if (Verbosity() > 0)
0758           {
0759             std::cout << "SINGLE CELL FIRED: " << intphibin << " " << intzbin << std::endl;
0760           }
0761           vphi.push_back(intphibin);
0762           vz.push_back(intzbin);
0763           vdedx.push_back(trklen);
0764         }
0765         else
0766         {
0767           double phistep_half = geo->get_phistep() / 2.;
0768           double zstep_half = geo->get_zstep() / 2.;
0769           for (int ibp = intphibin; ibp <= intphibinout; ibp++)
0770           {
0771             double cx = geo->get_phicenter(ibp) - phistep_half;
0772             double dx = geo->get_phicenter(ibp) + phistep_half;
0773             for (int ibz = intzbin; ibz <= intzbinout; ibz++)
0774             {
0775               double cy = geo->get_zcenter(ibz) - zstep_half;
0776               double dy = geo->get_zcenter(ibz) + zstep_half;
0777               // std::cout << "##### line: " << ax << " " << ay << " " << bx << " " << by << std::endl;
0778               // std::cout << "####### cell: " << cx << " " << cy << " " << dx << " " << dy << std::endl;
0779               std::pair<bool, double> intersect = PHG4Utils::line_and_rectangle_intersect(ax, ay, bx, by, cx, cy, dx, dy);
0780               if (intersect.first)
0781               {
0782                 if (Verbosity() > 0)
0783                 {
0784                   std::cout << "CELL FIRED: " << ibp << " " << ibz << " " << intersect.second << std::endl;
0785                 }
0786                 vphi.push_back(ibp);
0787                 vz.push_back(ibz);
0788                 vdedx.push_back(intersect.second);
0789               }
0790             }
0791           }
0792         }
0793         if (Verbosity() > 0)
0794         {
0795           std::cout << "NUMBER OF FIRED CELLS = " << vz.size() << std::endl;
0796         }
0797 
0798         double tmpsum = 0.;
0799         for (unsigned int ii = 0; ii < vz.size(); ii++)
0800         {
0801           tmpsum += vdedx[ii];
0802           vdedx[ii] = vdedx[ii] / trklen;
0803           if (Verbosity() > 0)
0804           {
0805             std::cout << "  CELL " << ii << "  dE/dX = " << vdedx[ii] << std::endl;
0806           }
0807         }
0808         if (Verbosity() > 0)
0809         {
0810           std::cout << "    TOTAL TRACK LENGTH = " << tmpsum << " " << trklen << std::endl;
0811         }
0812 
0813         for (unsigned int i1 = 0; i1 < vphi.size(); i1++)  // loop over all fired cells
0814         {
0815           int iphibin = vphi[i1];
0816           int izbin = vz[i1];
0817 
0818           unsigned long long tmp = iphibin;
0819           unsigned long long key = tmp << 32U;  // clang-tidy: U for unsigned int
0820           key += izbin;
0821           if (Verbosity() > 1)
0822           {
0823             std::cout << " iphibin " << iphibin << " izbin " << izbin << " key 0x" << std::hex << key << std::dec << std::endl;
0824           }
0825           // check to see if there is already an entry for this cell
0826           PHG4Cell *cell = nullptr;
0827           it = cellptmap.find(key);
0828 
0829           if (it != cellptmap.end())
0830           {
0831             cell = it->second;
0832             if (Verbosity() > 1)
0833             {
0834               std::cout << "  add energy to existing cell for key = " << cellptmap.find(key)->first << std::endl;
0835             }
0836 
0837             if (Verbosity() > 1 && hiter->second->has_property(PHG4Hit::prop_light_yield) && std::isnan(hiter->second->get_light_yield() * vdedx[i1]))
0838             {
0839               std::cout << "    NAN lighy yield with vdedx[i1] = " << vdedx[i1]
0840                         << " and hiter->second->get_light_yield() = " << hiter->second->get_light_yield() << std::endl;
0841             }
0842           }
0843           else
0844           {
0845             if (Verbosity() > 1)
0846             {
0847               std::cout << "    did not find a previous entry for key = 0x" << std::hex << key << std::dec << " create a new one" << std::endl;
0848             }
0849             PHG4CellDefs::keytype cellkey = PHG4CellDefs::SizeBinning::genkey(*layer, izbin, iphibin);
0850             cell = new PHG4Cellv1(cellkey);
0851             cellptmap[key] = cell;
0852           }
0853           if (!std::isfinite(hiter->second->get_edep() * vdedx[i1]))
0854           {
0855             std::cout << "hit 0x" << std::hex << hiter->first << std::dec << " not finite, edep: "
0856                       << hiter->second->get_edep() << " weight " << vdedx[i1] << std::endl;
0857           }
0858           cell->add_edep(hiter->first, hiter->second->get_edep() * vdedx[i1]);
0859           cell->add_edep(hiter->second->get_edep() * vdedx[i1]);  // add edep to cell
0860           if (hiter->second->has_property(PHG4Hit::prop_light_yield))
0861           {
0862             cell->add_light_yield(hiter->second->get_light_yield() * vdedx[i1]);
0863             if (Verbosity() > 1 && !std::isfinite(hiter->second->get_light_yield() * vdedx[i1]))
0864             {
0865               std::cout << "    NAN lighy yield with vdedx[i1] = " << vdedx[i1]
0866                         << " and hiter->second->get_light_yield() = " << hiter->second->get_light_yield() << std::endl;
0867             }
0868           }
0869           cell->add_shower_edep(hiter->second->get_shower_id(), hiter->second->get_edep() * vdedx[i1]);
0870         }
0871         vphi.clear();
0872         vz.clear();
0873 
0874       }  // end loop over hits
0875 
0876       int numcells = 0;
0877 
0878       for (it = cellptmap.begin(); it != cellptmap.end(); ++it)
0879       {
0880         cells->AddCell(it->second);
0881         numcells++;
0882         if (Verbosity() > 1)
0883         {
0884           std::cout << "Adding cell for key " << std::hex << it->first << std::dec << " in bin phi: " << PHG4CellDefs::SizeBinning::get_phibin(it->second->get_cellid())
0885                     << " phi: " << geo->get_phicenter(PHG4CellDefs::SizeBinning::get_phibin(it->second->get_cellid())) * 180. / M_PI
0886                     << ", z bin: " << PHG4CellDefs::SizeBinning::get_zbin(it->second->get_cellid())
0887                     << ", z: " << geo->get_zcenter(PHG4CellDefs::SizeBinning::get_zbin(it->second->get_cellid()))
0888                     << ", energy dep: " << it->second->get_edep()
0889                     << std::endl;
0890         }
0891       }
0892 
0893       if (Verbosity() > 0)
0894       {
0895         std::cout << "found " << numcells << " z/phi cells with energy deposition" << std::endl;
0896       }
0897     }
0898 
0899     //==========================================================
0900     // now reset the cell map before moving on to the next layer
0901     if (Verbosity() > 1)
0902     {
0903       std::cout << "cellptmap for layer " << *layer << " has final length " << cellptmap.size();
0904     }
0905     while (cellptmap.begin() != cellptmap.end())
0906     {
0907       // Assumes that memmory is freed by the cylinder cell container when it is destroyed
0908       cellptmap.erase(cellptmap.begin());
0909     }
0910     if (Verbosity() > 1)
0911     {
0912       std::cout << " reset it to " << cellptmap.size() << std::endl;
0913     }
0914   }
0915   if (chkenergyconservation)
0916   {
0917     CheckEnergy(topNode);
0918   }
0919 
0920   return Fun4AllReturnCodes::EVENT_OK;
0921 }
0922 
0923 void PHG4CylinderCellReco::cellsize(const int detid, const double sr, const double sz)
0924 {
0925   if (binning.find(detid) != binning.end())
0926   {
0927     std::cout << "size for layer " << detid << " already set" << std::endl;
0928     return;
0929   }
0930   binning[detid] = PHG4CellDefs::sizebinning;
0931   set_double_param(detid, "size_long", sz);
0932   set_double_param(detid, "size_perp", sr);
0933 }
0934 
0935 void PHG4CylinderCellReco::etaphisize(const int detid, const double deltaeta, const double deltaphi)
0936 {
0937   if (binning.find(detid) != binning.end())
0938   {
0939     std::cout << "size for layer " << detid << " already set" << std::endl;
0940     return;
0941   }
0942   binning[detid] = PHG4CellDefs::etaphibinning;
0943   set_double_param(detid, "size_long", deltaeta);
0944   set_double_param(detid, "size_perp", deltaphi);
0945   return;
0946 }
0947 
0948 void PHG4CylinderCellReco::set_size(const int i, const double sizeA, const double sizeB)
0949 {
0950   cell_size[i] = std::make_pair(sizeA, sizeB);
0951   return;
0952 }
0953 
0954 void PHG4CylinderCellReco::set_timing_window(const int detid, const double tmin, const double tmax)
0955 {
0956   set_double_param(detid, "tmin", tmin);
0957   set_double_param(detid, "tmax", tmax);
0958   return;
0959 }
0960 
0961 int PHG4CylinderCellReco::CheckEnergy(PHCompositeNode *topNode)
0962 {
0963   PHG4CellContainer *cells = findNode::getClass<PHG4CellContainer>(topNode, cellnodename);
0964   double sum_energy_cells = 0.;
0965   double sum_energy_stored_hits = 0.;
0966   double sum_energy_stored_showers = 0.;
0967   PHG4CellContainer::ConstRange cell_begin_end = cells->getCells();
0968   PHG4CellContainer::ConstIterator citer;
0969   for (citer = cell_begin_end.first; citer != cell_begin_end.second; ++citer)
0970   {
0971     sum_energy_cells += citer->second->get_edep();
0972     PHG4Cell::EdepConstRange cellrange = citer->second->get_g4hits();
0973     for (PHG4Cell::EdepConstIterator iter = cellrange.first; iter != cellrange.second; ++iter)
0974     {
0975       sum_energy_stored_hits += iter->second;
0976     }
0977     PHG4Cell::ShowerEdepConstRange shwrrange = citer->second->get_g4showers();
0978     for (PHG4Cell::ShowerEdepConstIterator iter = shwrrange.first; iter != shwrrange.second; ++iter)
0979     {
0980       sum_energy_stored_showers += iter->second;
0981     }
0982   }
0983   // the fractional eloss for particles traversing eta bins leads to minute rounding errors
0984   if (sum_energy_stored_hits > 0)
0985   {
0986     if (fabs(sum_energy_cells - sum_energy_stored_hits) / sum_energy_cells > 1e-6)
0987     {
0988       std::cout << "energy mismatch between cell energy " << sum_energy_cells
0989                 << " and stored hit energies " << sum_energy_stored_hits
0990                 << std::endl;
0991     }
0992   }
0993   if (sum_energy_stored_showers > 0)
0994   {
0995     if (fabs(sum_energy_cells - sum_energy_stored_showers) / sum_energy_cells > 1e-6)
0996     {
0997       std::cout << "energy mismatch between cell energy " << sum_energy_cells
0998                 << " and stored shower energies " << sum_energy_stored_showers
0999                 << std::endl;
1000     }
1001   }
1002   if (fabs(sum_energy_cells - sum_energy_g4hit) / sum_energy_g4hit > 1e-6)
1003   {
1004     std::cout << "energy mismatch between cells: " << sum_energy_cells
1005               << " and hits: " << sum_energy_g4hit
1006               << " diff sum(cells) - sum(hits): " << sum_energy_cells - sum_energy_g4hit
1007               << " cut val " << fabs(sum_energy_cells - sum_energy_g4hit) / sum_energy_g4hit
1008               << std::endl;
1009     std::cout << Name() << ":total energy for this event: " << sum_energy_g4hit << " GeV" << std::endl;
1010     std::cout << Name() << ": sum cell energy: " << sum_energy_cells << " GeV" << std::endl;
1011     std::cout << Name() << ": sum shower energy: " << sum_energy_stored_showers << " GeV" << std::endl;
1012     std::cout << Name() << ": sum stored hit energy: " << sum_energy_stored_hits << " GeV" << std::endl;
1013     std::cout << Name() << ": hit energy before cuts: " << sum_energy_before_cuts << " GeV" << std::endl;
1014     return -1;
1015   }
1016   else
1017   {
1018     if (Verbosity() > 0)
1019     {
1020       std::cout << Name() << ":total energy for this event: " << sum_energy_g4hit << " GeV" << std::endl;
1021       std::cout << Name() << ": sum cell energy: " << sum_energy_cells << " GeV" << std::endl;
1022       std::cout << Name() << ": sum shower energy: " << sum_energy_stored_showers << " GeV" << std::endl;
1023       std::cout << Name() << ": sum stored hit energy: " << sum_energy_stored_hits << " GeV" << std::endl;
1024       std::cout << Name() << ": hit energy before cuts: " << sum_energy_before_cuts << " GeV" << std::endl;
1025     }
1026   }
1027   return 0;
1028 }
1029 
1030 void PHG4CylinderCellReco::Detector(const std::string &d)
1031 {
1032   detector = d;
1033   // only set the outdetector nodename if it wasn't set already
1034   // in case the order in the macro is outdetector(); detector();
1035   if (outdetector.size() == 0)
1036   {
1037     OutputDetector(d);
1038   }
1039   return;
1040 }
1041 
1042 void PHG4CylinderCellReco::SetDefaultParameters()
1043 {
1044   set_default_double_param("size_long", NAN);
1045   set_default_double_param("size_perp", NAN);
1046   set_default_double_param("tmax", 60.0);
1047   set_default_double_param("tmin", -20.0);  // collision has a timing spread around the triggered event. Accepting negative time too.
1048   set_default_double_param("delta_t", 100.);
1049   return;
1050 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!