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