packages/globalvertex/GlobalVertexReco.cc

0001 #include "GlobalVertexReco.h"
0002
0003 //#include "GlobalVertex.h"     // for GlobalVertex, GlobalVe...
0004 #include "GlobalVertexMap.h"  // for GlobalVertexMap
0005 #include "GlobalVertexMapv1.h"
0006 #include "GlobalVertexv2.h"
0007 #include "MbdVertex.h"
0008 #include "MbdVertexMap.h"
0009 #include "CaloVertex.h"
0010 #include "CaloVertexMap.h"
0011 #include "SvtxVertex.h"
0012 #include "SvtxVertexMap.h"
0013 #include "TruthVertex.h"
0014 #include "TruthVertexMap.h"
0015 #include "TruthVertexMap_v1.h"
0016 #include "TruthVertex_v1.h"
0017
0018 #include <trackbase_historic/SvtxTrack.h>
0019 #include <trackbase_historic/SvtxTrackMap.h>
0020
0021 #include <fun4all/Fun4AllReturnCodes.h>
0022 #include <fun4all/SubsysReco.h>  // for SubsysReco
0023
0024 #include <phool/PHCompositeNode.h>
0025 #include <phool/PHIODataNode.h>
0026 #include <phool/PHNode.h>  // for PHNode
0027 #include <phool/PHNodeIterator.h>
0028 #include <phool/PHObject.h>  // for PHObject
0029 #include <phool/getClass.h>
0030 #include <phool/phool.h>  // for PHWHERE
0031
0032 // truth info, for truth z vertex
0033 #include <g4main/PHG4TruthInfoContainer.h>
0034 #include <g4main/PHG4VtxPoint.h>
0035
0036 #include <cmath>
0037 #include <cstdlib>  // for exit
0038 #include <iostream>
0039 #include <limits>
0040 #include <set>      // for _Rb_tree_const_iterator
0041 #include <utility>  // for pair
0042
0043 GlobalVertexReco::GlobalVertexReco(const std::string &name)
0044   : SubsysReco(name)
0045 {
0046 }
0047
0048 int GlobalVertexReco::InitRun(PHCompositeNode *topNode)
0049 {
0050   if (Verbosity() > 0)
0051   {
0052     std::cout << "======================= GlobalVertexReco::InitRun() =======================" << std::endl;
0053     std::cout << "===========================================================================" << std::endl;
0054   }
0055
0056   return CreateNodes(topNode);
0057 }
0058
0059 int GlobalVertexReco::process_event(PHCompositeNode *topNode)
0060 {
0061   if (Verbosity() > 1)
0062   {
0063     std::cout << "GlobalVertexReco::process_event -- entered" << std::endl;
0064   }
0065
0066   //---------------------------------
0067   // Get Objects off of the Node Tree
0068   //---------------------------------
0069   GlobalVertexMap *globalmap = findNode::getClass<GlobalVertexMap>(topNode, "GlobalVertexMap");
0070
0071   if (!globalmap)
0072   {
0073     std::cout << PHWHERE << "::ERROR - cannot find GlobalVertexMap" << std::endl;
0074     exit(-1);
0075   }
0076
0077   SvtxVertexMap *svtxmap = findNode::getClass<SvtxVertexMap>(topNode, "SvtxVertexMap");
0078   MbdVertexMap *mbdmap = findNode::getClass<MbdVertexMap>(topNode, "MbdVertexMap");
0079   CaloVertexMap *calomap = findNode::getClass<CaloVertexMap>(topNode, "CaloVertexMap");
0080   SvtxTrackMap *trackmap = findNode::getClass<SvtxTrackMap>(topNode, "SvtxTrackMap");
0081   TruthVertexMap *truthmap = findNode::getClass<TruthVertexMap>(topNode, "TruthVertexMap");
0082   PHG4TruthInfoContainer *truthinfo = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
0083
0084   // we will make 4 different kinds of global vertexes
0085   //  (1) SVTX+MBD vertexes - we match SVTX vertex to the nearest MBD vertex within 3 sigma in zvertex
0086   //      the spatial point comes from the SVTX, the timing from the MBD
0087   //      number of SVTX+MBD vertexes <= number of SVTX vertexes
0088
0089   //  (2) SVTX only vertexes - for those cases where the MBD wasn't simulated,
0090   //      or all the MBD vertexes are outside the 3 sigma matching requirement
0091   //      we pass forward the 3d point from the SVTX with the default timing info
0092
0093   //  (3) MBD only vertexes - use the default x,y positions on this module and
0094   //      pull in the mbd z and mbd t
0095
0096   //  (4) Truth vertexes - if the G4TruthInfo node is present (should be for simulation only), we will
0097   //      create a GlobalVertex from the primary truth vertex, and create+fill a TruthVertexMap node
0098
0099   // there may be some quirks as we get to large luminosity and the MBD becomes
0100   // untrust worthy, I'm guessing analyzers would resort exclusively to (1) or (2)
0101   // in those cases
0102
0103   std::set<unsigned int> used_svtx_vtxids;
0104   std::set<unsigned int> used_mbd_vtxids;
0105   std::set<unsigned int> used_calo_vtxids;
0106
0107   if (svtxmap && mbdmap && useVertexType(GlobalVertex::VTXTYPE::SVTX_MBD))
0108   {
0109     if (Verbosity())
0110     {
0111       std::cout << "GlobalVertexReco::process_event - svtxmap && mbdmap" << std::endl;
0112     }
0113
0114     for (SvtxVertexMap::ConstIter svtxiter = svtxmap->begin();
0115          svtxiter != svtxmap->end();
0116          ++svtxiter)
0117     {
0118       const SvtxVertex *svtx = svtxiter->second;
0119
0120       const MbdVertex *mbd_best = nullptr;
0121       float min_sigma = std::numeric_limits<float>::max();
0122       for (MbdVertexMap::ConstIter mbditer = mbdmap->begin();
0123            mbditer != mbdmap->end();
0124            ++mbditer)
0125       {
0126         const MbdVertex *mbd = mbditer->second;
0127
0128         float combined_error = sqrt(svtx->get_error(2, 2) + pow(mbd->get_z_err(), 2));
0129         float sigma = std::fabs(svtx->get_z() - mbd->get_z()) / combined_error;
0130         if (sigma < min_sigma)
0131         {
0132           min_sigma = sigma;
0133           mbd_best = mbd;
0134         }
0135       }
0136
0137       if (min_sigma > 3.0 || !mbd_best)
0138       {
0139         continue;
0140       }
0141
0142       // we have a matching pair
0143       GlobalVertex *vertex = new GlobalVertexv2();
0144       vertex->set_id(globalmap->size());
0145
0146       vertex->clone_insert_vtx(GlobalVertex::SVTX, svtx);
0147       vertex->clone_insert_vtx(GlobalVertex::MBD, mbd_best);
0148       used_svtx_vtxids.insert(svtx->get_id());
0149       used_mbd_vtxids.insert(mbd_best->get_id());
0150       vertex->set_id(globalmap->size());
0151
0152       globalmap->insert(vertex);
0153
0154       //! Reset track ids to the new vertex object
0155       if (trackmap)
0156       {
0157         for (auto iter = svtx->begin_tracks(); iter != svtx->end_tracks();
0158              ++iter)
0159         {
0160           auto *track = trackmap->find(*iter)->second;
0161           track->set_vertex_id(vertex->get_id());
0162         }
0163       }
0164
0165       if (Verbosity() > 1)
0166       {
0167         vertex->identify();
0168       }
0169     }
0170   }
0171
0172   // okay now loop over all unused SVTX vertexes (2nd class)...
0173   if (svtxmap && useVertexType(GlobalVertex::VTXTYPE::SVTX))
0174   {
0175     if (Verbosity())
0176     {
0177       std::cout << "GlobalVertexReco::process_event - svtxmap " << std::endl;
0178     }
0179
0180     for (SvtxVertexMap::ConstIter svtxiter = svtxmap->begin();
0181          svtxiter != svtxmap->end();
0182          ++svtxiter)
0183     {
0184       const SvtxVertex *svtx = svtxiter->second;
0185
0186       if (used_svtx_vtxids.contains(svtx->get_id()))
0187       {
0188         continue;
0189       }
0190       if (std::isnan(svtx->get_z()))
0191       {
0192         continue;
0193       }
0194
0195       // we have a standalone SVTX vertex
0196       GlobalVertex *vertex = new GlobalVertexv2();
0197
0198       vertex->set_id(globalmap->size());
0199
0200       vertex->clone_insert_vtx(GlobalVertex::SVTX, svtx);
0201       used_svtx_vtxids.insert(svtx->get_id());
0202
0203       //! Reset track ids to the new vertex object
0204       if (trackmap)
0205       {
0206         for (auto iter = svtx->begin_tracks(); iter != svtx->end_tracks();
0207              ++iter)
0208         {
0209           auto *track = trackmap->find(*iter)->second;
0210           track->set_vertex_id(vertex->get_id());
0211         }
0212       }
0213       globalmap->insert(vertex);
0214
0215       if (Verbosity() > 1)
0216       {
0217         vertex->identify();
0218       }
0219     }
0220   }
0221
0222   // okay now loop over all unused MBD vertexes (3rd class)...
0223   if (mbdmap && useVertexType(GlobalVertex::VTXTYPE::MBD))
0224   {
0225     if (Verbosity())
0226     {
0227       std::cout << "GlobalVertexReco::process_event -  mbdmap" << std::endl;
0228     }
0229
0230     for (MbdVertexMap::ConstIter mbditer = mbdmap->begin();
0231          mbditer != mbdmap->end();
0232          ++mbditer)
0233     {
0234       const MbdVertex *mbd = mbditer->second;
0235
0236       if (used_mbd_vtxids.contains(mbd->get_id()))
0237       {
0238         continue;
0239       }
0240
0241       if (std::isnan(mbd->get_z()))
0242       {
0243         continue;
0244       }
0245
0246       GlobalVertex *vertex = new GlobalVertexv2();
0247       vertex->set_id(globalmap->size());
0248
0249       vertex->clone_insert_vtx(GlobalVertex::MBD, mbd);
0250       used_mbd_vtxids.insert(mbd->get_id());
0251
0252       globalmap->insert(vertex);
0253
0254       if (Verbosity() > 1)
0255       {
0256         vertex->identify();
0257       }
0258     }
0259   }
0260
0261   // okay now loop over all unused MBD vertexes (3rd class)...
0262   if (calomap && useVertexType(GlobalVertex::VTXTYPE::CALO))
0263   {
0264     if (Verbosity())
0265     {
0266       std::cout << "GlobalVertexReco::process_event -  calomap" << std::endl;
0267     }
0268
0269     for (CaloVertexMap::ConstIter caloiter = calomap->begin();
0270          caloiter != calomap->end();
0271          ++caloiter)
0272       {
0273     const CaloVertex *calo = caloiter->second;
0274
0275     if (used_calo_vtxids.contains(calo->get_id()))
0276       {
0277         continue;
0278       }
0279
0280     if (std::isnan(calo->get_z()))
0281       {
0282         continue;
0283       }
0284
0285     GlobalVertex *vertex = new GlobalVertexv2();
0286     vertex->set_id(globalmap->size());
0287
0288     vertex->clone_insert_vtx(GlobalVertex::CALO, calo);
0289     used_calo_vtxids.insert(calo->get_id());
0290
0291     globalmap->insert(vertex);
0292
0293     if (Verbosity() > 1)
0294       {
0295         vertex->identify();
0296       }
0297       }
0298   }
0299
0300 // okay now loop over all unused MBD vertexes (3rd class)...
0301   if (mbdmap && calomap && useVertexType(GlobalVertex::VTXTYPE::MBD_CALO))
0302   {
0303     if (Verbosity())
0304     {
0305       std::cout << "GlobalVertexReco::process_event -  calomap + mbdmap" << std::endl;
0306     }
0307
0308     for (MbdVertexMap::ConstIter mbditer = mbdmap->begin();
0309          mbditer != mbdmap->end();
0310          ++mbditer)
0311     {
0312       const MbdVertex *mbd = mbditer->second;
0313
0314       if (used_mbd_vtxids.contains(mbd->get_id()))
0315       {
0316         continue;
0317       }
0318
0319
0320       bool getcalo =  std::isnan(mbd->get_z());
0321
0322       if (getcalo)
0323     {
0324       for (CaloVertexMap::ConstIter caloiter = calomap->begin();
0325            caloiter != calomap->end();
0326            ++caloiter)
0327         {
0328           const CaloVertex *calo = caloiter->second;
0329
0330           if (used_calo_vtxids.contains(calo->get_id()))
0331         {
0332           continue;
0333         }
0334
0335           if (std::isnan(calo->get_z()))
0336         {
0337           continue;
0338         }
0339
0340           GlobalVertex *vertex = new GlobalVertexv2();
0341           vertex->set_id(globalmap->size());
0342
0343           vertex->clone_insert_vtx(GlobalVertex::CALO, calo);
0344           used_calo_vtxids.insert(calo->get_id());
0345
0346           globalmap->insert(vertex);
0347
0348           if (Verbosity() > 1)
0349         {
0350           vertex->identify();
0351         }
0352         }
0353
0354     }
0355       else
0356     {
0357       GlobalVertex *vertex = new GlobalVertexv2();
0358       vertex->set_id(globalmap->size());
0359
0360       vertex->clone_insert_vtx(GlobalVertex::MBD, mbd);
0361       used_mbd_vtxids.insert(mbd->get_id());
0362
0363       globalmap->insert(vertex);
0364
0365       if (Verbosity() > 1)
0366         {
0367           vertex->identify();
0368         }
0369     }
0370     }
0371   }
0372 // okay now add the truth vertex (4th class)...
0373
0374   if (truthinfo)
0375   {
0376     PHG4VtxPoint *vtxp = truthinfo->GetPrimaryVtx(truthinfo->GetPrimaryVertexIndex());
0377     float truth_z_vertex = std::numeric_limits<float>::quiet_NaN();
0378     float truth_x_vertex = std::numeric_limits<float>::quiet_NaN();
0379     float truth_y_vertex = std::numeric_limits<float>::quiet_NaN();
0380     if (vtxp)
0381     {
0382       truth_z_vertex = vtxp->get_z();
0383       truth_x_vertex = vtxp->get_x();
0384       truth_y_vertex = vtxp->get_y();
0385       TruthVertex *tvertex = new TruthVertex_v1();
0386       tvertex->set_id(globalmap->size());
0387       tvertex->set_z(truth_z_vertex);
0388       tvertex->set_z_err(0);
0389       tvertex->set_x(truth_x_vertex);
0390       tvertex->set_x_err(0);
0391       tvertex->set_y(truth_y_vertex);
0392       tvertex->set_y_err(0);
0393
0394       tvertex->set_t(0);
0395       tvertex->set_t_err(0);  // 0.1
0396       GlobalVertex *vertex = new GlobalVertexv2();
0397       vertex->clone_insert_vtx(GlobalVertex::TRUTH, tvertex);
0398       globalmap->insert(vertex);
0399       if (truthmap)
0400       {
0401         truthmap->insert(tvertex);
0402       }
0403       if (Verbosity() > 1)
0404       {
0405         vertex->identify();
0406       }
0407     }
0408     else if (Verbosity())
0409     {
0410       std::cout << "PHG4TruthInfoContainer has no primary vertex" << std::endl;
0411     }
0412   }
0413   else if (Verbosity())
0414   {
0415     std::cout << "PHG4TruthInfoContainer is missing" << std::endl;
0416   }
0417
0418   /// Associate any tracks that were not assigned a track-vertex
0419   if (trackmap)
0420   {
0421     for (const auto &[tkey, track] : *trackmap)
0422     {
0423       //! Check that the vertex hasn't already been assigned
0424       auto trackvtxid = track->get_vertex_id();
0425       if (globalmap->find(trackvtxid) != globalmap->end())
0426       {
0427         continue;
0428       }
0429
0430       float maxdz = std::numeric_limits<float>::max();
0431       unsigned int vtxid = std::numeric_limits<unsigned int>::max();
0432
0433       for (const auto &[vkey, vertex] : *globalmap)
0434       {
0435         float dz = track->get_z() - vertex->get_z();
0436         if (std::fabs(dz) < maxdz)
0437         {
0438           maxdz = dz;
0439           vtxid = vkey;
0440         }
0441       }
0442
0443       track->set_vertex_id(vtxid);
0444       if (Verbosity())
0445       {
0446         std::cout << "Associated track with z " << track->get_z() << " to GlobalVertex id " << track->get_vertex_id() << std::endl;
0447       }
0448     }
0449   }
0450
0451   if (Verbosity())
0452   {
0453     globalmap->identify();
0454   }
0455
0456   return Fun4AllReturnCodes::EVENT_OK;
0457 }
0458
0459 int GlobalVertexReco::CreateNodes(PHCompositeNode *topNode)
0460 {
0461   PHNodeIterator iter(topNode);
0462
0463   // Looking for the DST node
0464   PHCompositeNode *dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
0465   if (!dstNode)
0466   {
0467     std::cout << PHWHERE << "DST Node missing, doing nothing." << std::endl;
0468     return Fun4AllReturnCodes::ABORTRUN;
0469   }
0470
0471   // store the GLOBAL stuff under a sub-node directory
0472   PHCompositeNode *globalNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "GLOBAL"));
0473   if (!globalNode)
0474   {
0475     globalNode = new PHCompositeNode("GLOBAL");
0476     dstNode->addNode(globalNode);
0477   }
0478
0479   // create the GlobalVertexMap
0480   GlobalVertexMap *vertexes = findNode::getClass<GlobalVertexMap>(topNode, "GlobalVertexMap");
0481   if (!vertexes)
0482   {
0483     vertexes = new GlobalVertexMapv1();
0484     PHIODataNode<PHObject> *VertexMapNode = new PHIODataNode<PHObject>(vertexes, "GlobalVertexMap", "PHObject");
0485     globalNode->addNode(VertexMapNode);
0486   }
0487
0488   TruthVertexMap *truthmap = findNode::getClass<TruthVertexMap>(topNode, "TruthVertexMap");
0489   PHG4TruthInfoContainer *truthinfo = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
0490   // Create the TruthVertexMap only if the PHG4TruthInfoContainer is present
0491   if (!truthmap && truthinfo)
0492   {
0493     if (Verbosity())
0494     {
0495       std::cout << "Creating TruthVertexMap node" << std::endl;
0496     }
0497     truthmap = new TruthVertexMap_v1();
0498     PHIODataNode<PHObject> *TruthVertexMapNode = new PHIODataNode<PHObject>(truthmap, "TruthVertexMap", "PHObject");
0499     globalNode->addNode(TruthVertexMapNode);
0500   }
0501
0502   return Fun4AllReturnCodes::EVENT_OK;
0503 }