g4simulation/g4main/PHG4InEvent.cc

0001 #include "PHG4InEvent.h"
0002
0003 #include "PHG4Particle.h"
0004 #include "PHG4VtxPoint.h"
0005 #include "PHG4VtxPointv1.h"
0006
0007 #include <algorithm>
0008 #include <cmath>
0009 #include <cstdlib>
0010
0011 PHG4InEvent::~PHG4InEvent()
0012 {
0013   vtxlist.clear();
0014   particlelist.clear();
0015   return;
0016 }
0017
0018 int PHG4InEvent::AddVtx(const int /*id*/, const PHG4VtxPoint &vtx)
0019 {
0020   PHG4VtxPoint *newvtx = new PHG4VtxPoint(vtx);
0021   int iret = AddVtxCommon(newvtx);
0022   return iret;
0023 }
0024
0025 int PHG4InEvent::AddVtxHepMC(const int id, const double x, const double y, const double z, const double t0 = std::numeric_limits<double>::quiet_NaN())
0026 {
0027   PHG4VtxPoint *newvtx = new PHG4VtxPointv1(x, y, z, t0);
0028   vtxlist[id] = newvtx;
0029   return 0;
0030 }
0031
0032 int PHG4InEvent::AddVtx(const double x, const double y, const double z, const double t0 = std::numeric_limits<double>::quiet_NaN())
0033 {
0034   PHG4VtxPoint *newvtx = new PHG4VtxPointv1(x, y, z, t0);
0035   int iret = AddVtxCommon(newvtx);
0036   return iret;
0037 }
0038
0039 int PHG4InEvent::AddVtxCommon(PHG4VtxPoint *newvtx)
0040 {
0041   std::pair<std::map<int, PHG4VtxPoint *>::const_iterator, std::map<int, PHG4VtxPoint *>::const_iterator> vtxbegin_end = GetVertices();
0042   for (std::map<int, PHG4VtxPoint *>::const_iterator viter = vtxbegin_end.first; viter != vtxbegin_end.second; ++viter)
0043   {
0044     /// Do a simultaneous relative + absolute floating point comparison for
0045     /// the vertex four vectors
0046     /// make the distance comparison 1 micron, time comparison 0.1 ps
0047     const auto distepsilon = 0.0001;
0048     const auto timeepsilon = 0.0001;
0049     auto xlist = {1.0, std::abs(newvtx->get_x()), std::abs(viter->second->get_x())};
0050     auto ylist = {1.0, std::abs(newvtx->get_y()), std::abs(viter->second->get_y())};
0051     auto zlist = {1.0, std::abs(newvtx->get_z()), std::abs(viter->second->get_z())};
0052     auto tlist = {1.0, std::abs(newvtx->get_t()), std::abs(viter->second->get_t())};
0053
0054     bool xcomp = std::abs(newvtx->get_x() - viter->second->get_x()) <= distepsilon * (*std::max_element(xlist.begin(), xlist.end()));
0055     bool ycomp = std::abs(newvtx->get_y() - viter->second->get_y()) <= distepsilon * (*std::max_element(ylist.begin(), ylist.end()));
0056     bool zcomp = std::abs(newvtx->get_z() - viter->second->get_z()) <= distepsilon * (*std::max_element(zlist.begin(), zlist.end()));
0057     bool tcomp = std::abs(newvtx->get_t() - viter->second->get_t()) <= timeepsilon * (*std::max_element(tlist.begin(), tlist.end()));
0058
0059     if (*newvtx == *(viter->second) ||
0060         (xcomp && ycomp && zcomp && tcomp))
0061     {
0062       delete newvtx;
0063       return viter->first;
0064     }
0065   }
0066   int id = GetNVtx() + 1;
0067   vtxlist[id] = newvtx;
0068   return id;
0069 }
0070
0071 int PHG4InEvent::AddParticle(const int vtxid, PHG4Particle *particle)
0072 {
0073   if (!vtxlist.contains(vtxid))
0074   {
0075     std::cout << "cannot add particle to non existing vertex, id: " << vtxid << std::endl;
0076     exit(1);
0077   }
0078   // checking for duplicate particles - sometimes interesting
0079   //  std::pair< std::multimap<int,PHG4Particle *>::const_iterator, std::multimap<int,PHG4Particle *>::const_iterator > particles = GetParticles(vtxid);
0080
0081   //   for (multimap<int,PHG4Particle *>::const_iterator piter = particles.first; piter != particles.second; piter++)
0082   //     {
0083   //       if (*particle == *(piter->second))
0084   //    {
0085   //      cout << PHWHERE << "Particle already added (same pid and momentum), dropping it since duplication will cause problems down the line particle:" << endl;
0086   //      particle->identify();
0087   //      delete particle;
0088   //      return -1;
0089   //    }
0090   //     }
0091   particlelist.insert(std::pair<int, PHG4Particle *>(vtxid, particle));
0092   return 0;
0093 }
0094
0095 void PHG4InEvent::Reset()
0096 {
0097   embedded_particlelist.clear();  // just pointers - we can clear it without deleting
0098   while (vtxlist.begin() != vtxlist.end())
0099   {
0100     delete vtxlist.begin()->second;
0101     vtxlist.erase(vtxlist.begin());
0102   }
0103   while (particlelist.begin() != particlelist.end())
0104   {
0105     delete particlelist.begin()->second;
0106     particlelist.erase(particlelist.begin());
0107   }
0108   return;
0109 }
0110
0111 std::pair<std::map<int, PHG4VtxPoint *>::const_iterator, std::map<int, PHG4VtxPoint *>::const_iterator>
0112 PHG4InEvent::GetVertices() const
0113 {
0114   std::pair<std::map<int, PHG4VtxPoint *>::const_iterator, std::map<int, PHG4VtxPoint *>::const_iterator> retpair(vtxlist.begin(), vtxlist.end());
0115   return retpair;
0116 }
0117
0118 std::pair<std::multimap<int, PHG4Particle *>::const_iterator, std::multimap<int, PHG4Particle *>::const_iterator>
0119 PHG4InEvent::GetParticles(const int vtxid) const
0120 {
0121   std::pair<std::multimap<int, PHG4Particle *>::const_iterator, std::multimap<int, PHG4Particle *>::const_iterator> retpair(particlelist.lower_bound(vtxid), particlelist.upper_bound(vtxid));
0122   return retpair;
0123 }
0124
0125 std::pair<std::multimap<int, PHG4Particle *>::const_iterator, std::multimap<int, PHG4Particle *>::const_iterator>
0126 PHG4InEvent::GetParticles() const
0127 {
0128   std::pair<std::multimap<int, PHG4Particle *>::const_iterator, std::multimap<int, PHG4Particle *>::const_iterator> retpair(particlelist.begin(), particlelist.end());
0129   return retpair;
0130 }
0131
0132 std::pair<std::multimap<int, PHG4Particle *>::iterator, std::multimap<int, PHG4Particle *>::iterator>
0133 PHG4InEvent::GetParticles_Modify()
0134 {
0135   std::pair<std::multimap<int, PHG4Particle *>::iterator, std::multimap<int, PHG4Particle *>::iterator> retpair(particlelist.begin(), particlelist.end());
0136   return retpair;
0137 }
0138
0139 void PHG4InEvent::identify(std::ostream &os) const
0140 {
0141   os << "vtx: " << std::endl;
0142   std::multimap<int, PHG4Particle *>::const_iterator particle_iter;
0143   for (auto iter : vtxlist)
0144   {
0145     os << "vtx " << iter.first << " , ";
0146     iter.second->identify(os);
0147     std::pair<std::multimap<int, PHG4Particle *>::const_iterator, std::multimap<int, PHG4Particle *>::const_iterator> particlebegin_end = GetParticles(iter.first);
0148     for (particle_iter = particlebegin_end.first; particle_iter != particlebegin_end.second; ++particle_iter)
0149     {
0150       os << "vtx " << particle_iter->first << ", ";
0151       particle_iter->second->identify(os);
0152     }
0153   }
0154   if (!embedded_particlelist.empty())
0155   {
0156     os << "embedded particles:" << std::endl;
0157     for (auto iter : embedded_particlelist)
0158     {
0159       (iter.first)->identify(os);
0160     }
0161   }
0162   else
0163   {
0164     os << "no embedded particles" << std::endl;
0165   }
0166   return;
0167 }
0168
0169 int PHG4InEvent::isEmbeded(PHG4Particle *p) const
0170 {
0171   std::map<PHG4Particle *, int>::const_iterator iter = embedded_particlelist.find(p);
0172   if (iter == embedded_particlelist.end())
0173   {
0174     return 0;
0175   }
0176
0177   return iter->second;
0178 }
0179
0180 void PHG4InEvent::DeleteParticle(std::multimap<int, PHG4Particle *>::iterator &iter)
0181 {
0182   delete iter->second;
0183   particlelist.erase(iter);
0184 }