packages/trackreco/PHTruthVertexing.cc

0001 #include "PHTruthVertexing.h"
0002
0003 #include <globalvertex/SvtxVertexMap.h>
0004 #include <globalvertex/SvtxVertex_v1.h>
0005 #include <trackbase_historic/SvtxTrackMap.h>
0006 #include <g4main/PHG4TruthInfoContainer.h>
0007 #include <g4main/PHG4VtxPoint.h>
0008
0009 #include <fun4all/Fun4AllReturnCodes.h>
0010
0011 #include <phool/PHRandomSeed.h>
0012 #include <phool/getClass.h>
0013 #include <phool/phool.h>
0014
0015 // gsl
0016 #include <gsl/gsl_randist.h>
0017 #include <gsl/gsl_rng.h>
0018
0019 #include <iostream>                            // for operator<<, basic_ostream
0020 #include <set>                                 // for _Rb_tree_iterator, set
0021 #include <utility>                             // for pair
0022
0023 class PHCompositeNode;
0024
0025 using namespace std;
0026
0027 PHTruthVertexing::PHTruthVertexing(const std::string& name)
0028   : PHInitVertexing(name)
0029   , _g4truth_container(nullptr)
0030   , _vertex_error({0.0005, 0.0005, 0.0005})
0031   , _embed_only(false)
0032 {
0033
0034   m_RandomGenerator = gsl_rng_alloc(gsl_rng_mt19937);
0035
0036   unsigned int seed = PHRandomSeed();  // fixed seed is handled in this funtcion
0037   //  cout << Name() << " random seed: " << seed << endl;
0038   gsl_rng_set(m_RandomGenerator, seed);
0039 }
0040
0041
0042
0043 PHTruthVertexing::~PHTruthVertexing() {
0044   gsl_rng_free(m_RandomGenerator);
0045 }
0046
0047 int PHTruthVertexing::Setup(PHCompositeNode* topNode)
0048 {
0049   int ret = PHInitVertexing::Setup(topNode);
0050   if (ret != Fun4AllReturnCodes::EVENT_OK) return ret;
0051
0052   ret = GetNodes(topNode);
0053   if (ret != Fun4AllReturnCodes::EVENT_OK) return ret;
0054
0055   return Fun4AllReturnCodes::EVENT_OK;
0056 }
0057
0058 int PHTruthVertexing::Process(PHCompositeNode* topNode)
0059 {
0060   /*
0061   if(Verbosity() > 1)
0062     {
0063       std::cout << "embed only: " << std::boolalpha << _embed_only << std::endl;
0064     }
0065   */
0066
0067   // we just copy all vertices from the truth container to the SvtxVertexMap
0068
0069   auto vrange =  _g4truth_container->GetPrimaryVtxRange();
0070   for (auto iter = vrange.first; iter != vrange.second; ++iter)  // process all primary vertexes
0071     {
0072       const int point_id = iter->first;
0073       int gembed =  _g4truth_container->isEmbededVtx(point_id);
0074
0075       std::vector<float> pos;
0076       pos.clear();
0077       pos.assign(3, 0.0);
0078
0079       pos[0] = iter->second->get_x();
0080       pos[1] = iter->second->get_y();
0081       pos[2] = iter->second->get_z();
0082
0083       if(Verbosity() > 1)
0084     {
0085       cout << " PHTruthVertexing: gembed: " << gembed << " nominal position: " << " vx " << pos[0] << " vy " << pos[1] << " vz " << pos[2] << endl;
0086     }
0087
0088       // We take all primary vertices and copy them to the vertex map, regardless of track embedding ID
0089
0090       pos[0] += _vertex_error[0] * gsl_ran_ugaussian(m_RandomGenerator);
0091       pos[1] += _vertex_error[1] * gsl_ran_ugaussian(m_RandomGenerator);
0092       pos[2] += _vertex_error[2] * gsl_ran_ugaussian(m_RandomGenerator);
0093
0094
0095       if (Verbosity() > 0)
0096     {
0097       cout << __LINE__ << " PHTruthVertexing::Process: point_id " << point_id << " gembed " << gembed << "   {" << pos[0]
0098            << ", " << pos[1] << ", " << pos[2] << "} +- {"
0099            << _vertex_error[0] << ", " << _vertex_error[1] << ", "
0100            << _vertex_error[2] << "}" << endl;
0101     }
0102
0103       SvtxVertex* vertex = new SvtxVertex_v1();
0104
0105       vertex->set_x(pos[0]);
0106       vertex->set_y(pos[1]);
0107       vertex->set_z(pos[2]);
0108
0109       for (int j = 0; j < 3; ++j)
0110     {
0111       for (int i = j; i < 3; ++i)
0112         {
0113           vertex->set_error(i, j,
0114                 (i == j ? _vertex_error[i] * _vertex_error[i] : 0));
0115         }
0116     }
0117
0118       vertex->set_id(0);
0119       vertex->set_t(0);
0120       vertex->set_chisq(0);
0121       vertex->set_ndof(1);
0122       _vertex_map->insert(vertex);
0123
0124     }
0125
0126   if (Verbosity() > 0)
0127       _vertex_map->identify();
0128
0129   if(_associate_tracks)
0130     { assignTracksVertices(topNode); }
0131
0132   return Fun4AllReturnCodes::EVENT_OK;
0133 }
0134
0135 void PHTruthVertexing::assignTracksVertices(PHCompositeNode *topNode)
0136 {
0137   SvtxTrackMap * trackMap = findNode::getClass<SvtxTrackMap>(topNode, _track_map_name.c_str());
0138   if(!trackMap)
0139     {
0140       std::cout << PHWHERE
0141         << "Can't find requested track map. Exiting"
0142         << std::endl;
0143       return;
0144     }
0145
0146   for(SvtxTrackMap::Iter iter = trackMap->begin();
0147       iter != trackMap->end();
0148       ++iter)
0149     {
0150       auto trackKey = iter->first;
0151       auto track = iter->second;
0152
0153       if(Verbosity() > 3)
0154     track->identify();
0155
0156       const double trackZ = track->get_z();
0157
0158       double dz = 9999.;
0159       int trackVertexId = 9999;
0160
0161       for(SvtxVertexMap::Iter viter = _vertex_map->begin();
0162       viter != _vertex_map->end();
0163       ++viter)
0164     {
0165       auto vertexKey = viter->first;
0166       auto vertex = viter->second;
0167       if(Verbosity() > 3)
0168         vertex->identify();
0169
0170       const double vertexZ = vertex->get_z();
0171
0172       if( fabs(trackZ - vertexZ) < dz )
0173         {
0174           dz = fabs(trackZ - vertexZ);
0175           trackVertexId = vertexKey;
0176         }
0177
0178     }
0179
0180       track->set_vertex_id(trackVertexId);
0181       auto vertex = _vertex_map->find(trackVertexId)->second;
0182       vertex->insert_track(trackKey);
0183
0184     }
0185
0186
0187 }
0188
0189 int PHTruthVertexing::GetNodes(PHCompositeNode* topNode)
0190 {
0191   _g4truth_container = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
0192   if (!_g4truth_container)
0193   {
0194     cerr << PHWHERE << " ERROR: Can't find node G4TruthInfo" << endl;
0195     return Fun4AllReturnCodes::ABORTEVENT;
0196   }
0197   return Fun4AllReturnCodes::EVENT_OK;
0198 }
0199
0200 int PHTruthVertexing::End(PHCompositeNode * /*topNode*/)
0201 {
0202   return Fun4AllReturnCodes::EVENT_OK;
0203 }