Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2025-08-05 08:18:18 

0001 /*!
0002  *  \file       PHG4TrackFastSim.cc
0003  *  \brief      Kalman Filter based on smeared truth PHG4Hit
0004  *  \details    Kalman Filter based on smeared truth PHG4Hit
0005  *  \author     Haiwang Yu <yuhw@nmsu.edu>
0006  */
0007 
0008 #include "PHG4TrackFastSim.h"
0009 
0010 #include <phgenfit/Fitter.h>
0011 #include <phgenfit/Measurement.h>  // for Measurement
0012 #include <phgenfit/PlanarMeasurement.h>
0013 #include <phgenfit/SpacepointMeasurement.h>
0014 #include <phgenfit/Track.h>
0015 #include <phgeom/PHGeomUtility.h>
0016 
0017 #include <globalvertex/SvtxVertexMap_v1.h>
0018 #include <globalvertex/SvtxVertex_v1.h>
0019 #include <trackbase_historic/SvtxTrack.h>
0020 #include <trackbase_historic/SvtxTrackMap.h>
0021 #include <trackbase_historic/SvtxTrackMap_v1.h>
0022 #include <trackbase_historic/SvtxTrackState.h>
0023 #include <trackbase_historic/SvtxTrackState_v1.h>
0024 #include <trackbase_historic/SvtxTrack_FastSim_v1.h>
0025 
0026 #include <calobase/RawTowerGeom.h>
0027 #include <calobase/RawTowerGeomContainer.h>
0028 
0029 #include <phparameter/PHParameters.h>
0030 
0031 #include <phfield/PHFieldUtility.h>
0032 
0033 #include <g4main/PHG4Hit.h>
0034 #include <g4main/PHG4HitContainer.h>  // for PHG4HitContainer
0035 #include <g4main/PHG4Particle.h>
0036 #include <g4main/PHG4TruthInfoContainer.h>
0037 #include <g4main/PHG4VtxPoint.h>
0038 
0039 #include <fun4all/Fun4AllReturnCodes.h>
0040 #include <fun4all/SubsysReco.h>  // for SubsysReco
0041 
0042 #include <phool/PHCompositeNode.h>
0043 #include <phool/PHIODataNode.h>
0044 #include <phool/PHNode.h>  // for PHNode
0045 #include <phool/PHNodeIterator.h>
0046 #include <phool/PHObject.h>  // for PHObject
0047 #include <phool/PHRandomSeed.h>
0048 #include <phool/getClass.h>
0049 #include <phool/phool.h>
0050 
0051 #include <GenFit/AbsMeasurement.h>
0052 #include <GenFit/EventDisplay.h>
0053 #include <GenFit/MeasuredStateOnPlane.h>
0054 #include <GenFit/RKTrackRep.h>
0055 
0056 #include <GenFit/FitStatus.h>              // for FitStatus
0057 #include <GenFit/GFRaveTrackParameters.h>  // for GFRaveTrackParameters
0058 #include <GenFit/GFRaveVertex.h>
0059 #include <GenFit/GFRaveVertexFactory.h>
0060 #include <GenFit/Track.h>
0061 
0062 #include <TMatrixDSymfwd.h>  // for TMatrixDSym
0063 #include <TMatrixTSym.h>     // for TMatrixTSym
0064 #include <TMatrixTUtils.h>   // for TMatrixTRow
0065 #include <TSystem.h>
0066 #include <TVector3.h>     // for TVector3, operator*
0067 #include <TVectorDfwd.h>  // for TVectorD
0068 #include <TVectorT.h>     // for TVectorT
0069 
0070 #include <gsl/gsl_randist.h>
0071 #include <gsl/gsl_rng.h>
0072 
0073 #include <cassert>  // for assert
0074 #include <cmath>
0075 #include <iostream>  // for operator<<, basic_...
0076 #include <map>
0077 #include <memory>  // for unique_ptr, alloca...
0078 #include <utility>
0079 
0080 class PHField;
0081 class TGeoManager;
0082 namespace genfit
0083 {
0084   class AbsTrackRep;
0085 }  // namespace genfit
0086 
0087 // NOLINTNEXTLINE(bugprone-macro-parentheses)
0088 #define LogDebug(exp) if (Verbosity()) std::cout << "PHG4TrackFastSim (DEBUG): " << __FILE__ << ": " << __LINE__ << ": " << exp << "\n"
0089 // NOLINTNEXTLINE(bugprone-macro-parentheses)
0090 #define LogError(exp) std::cout << "PHG4TrackFastSim (ERROR): " << __FILE__ << ": " << __LINE__ << ": " << exp << "\n"
0091 // NOLINTNEXTLINE(bugprone-macro-parentheses)
0092 #define LogWarning(exp) std::cout << "PHG4TrackFastSim (WARNING): " << __FILE__ << ": " << __LINE__ << ": " << exp << "\n"
0093 
0094 // names of our implemented calorimeters where the projections are done
0095 // at 1/2 of their depth, not at the surface
0096 // this is used to avoid user added projections with identical names
0097 std::set<std::string> reserved_cylinder_projection_names{"CEMC", "HCALIN", "HCALOUT", "BECAL"};
0098 std::set<std::string> reserved_zplane_projection_names{"FEMC", "FHCAL", "EEMC", "EHCAL", "LFHCAL"};
0099 
0100 PHG4TrackFastSim::PHG4TrackFastSim(const std::string& name)
0101   : SubsysReco(name)
0102   , m_Fitter(nullptr)
0103   , m_RaveVertexFactory(nullptr)
0104   , m_TruthContainer(nullptr)
0105   , m_SvtxTrackMapOut(nullptr)
0106   , m_SvtxVertexMap(nullptr)
0107   , m_SubTopnodeName("SVTX")
0108   , m_TrackmapOutNodeName("SvtxTrackMap")
0109   , m_VertexingMethod("kalman-smoothing:1")
0110   , m_FitAlgoName("DafRef")  // was ("KalmanFitterRefTrack")
0111   , m_VertexMinNdf(10.)
0112   , m_VertexXYResolution(50E-4)
0113   , m_VertexZResolution(50E-4)
0114   , m_EventCnt(-1)
0115   , m_PrimaryAssumptionPid(211)
0116   , m_SmearingFlag(true)
0117   , m_DoEvtDisplayFlag(false)
0118   , m_UseVertexInFittingFlag(true)
0119   , m_PrimaryTrackingFlag(1)
0120   , m_DoVertexingFlag(false)
0121 {
0122   unsigned int seed = PHRandomSeed();  // fixed seed is handled in this funtcion
0123   m_RandomGenerator = gsl_rng_alloc(gsl_rng_mt19937);
0124   gsl_rng_set(m_RandomGenerator, seed);
0125 
0126   m_Parameter = new PHParameters(Name());
0127 }
0128 
0129 PHG4TrackFastSim::~PHG4TrackFastSim()
0130 {
0131   delete m_Fitter;
0132   delete m_RaveVertexFactory;
0133   gsl_rng_free(m_RandomGenerator);
0134 }
0135 
0136 int PHG4TrackFastSim::InitRun(PHCompositeNode* topNode)
0137 {
0138   m_EventCnt = -1;
0139 
0140   int ret = CreateNodes(topNode);
0141   if (ret != Fun4AllReturnCodes::EVENT_OK)
0142   {
0143     return ret;
0144   }
0145   ret = GetNodes(topNode);
0146   if (ret != Fun4AllReturnCodes::EVENT_OK)
0147   {
0148     return ret;
0149   }
0150   //
0151   //
0152   TGeoManager* tgeo_manager = PHGeomUtility::GetTGeoManager(topNode);
0153   PHField* field = PHFieldUtility::GetFieldMapNode(nullptr, topNode);
0154 
0155   m_Fitter = PHGenFit::Fitter::getInstance(tgeo_manager,
0156                                            field, m_FitAlgoName, "RKTrackRep",
0157                                            m_DoEvtDisplayFlag);
0158 
0159   if (!m_Fitter)
0160   {
0161     std::cout << PHWHERE << std::endl;
0162     return Fun4AllReturnCodes::ABORTRUN;
0163   }
0164 
0165   m_Fitter->set_verbosity(Verbosity());
0166 
0167   // tower geometry for track states
0168 
0169   for (auto& iter : m_ProjectionsMap)
0170   {
0171     if (isfinite(iter.second.second))
0172     {
0173       continue;
0174     }
0175     switch (iter.second.first)
0176     {
0177     case DETECTOR_TYPE::Cylinder:
0178     {
0179       std::string nodename = "TOWERGEOM_" + iter.first;
0180       RawTowerGeomContainer* geo = findNode::getClass<RawTowerGeomContainer>(topNode, nodename);
0181       if (geo)
0182       {
0183         iter.second.second = geo->get_radius();
0184       }
0185       break;
0186     }
0187     case DETECTOR_TYPE::Vertical_Plane:
0188     {
0189       std::string towergeonodename = "TOWERGEOM_" + iter.first;
0190       RawTowerGeomContainer* towergeo = findNode::getClass<RawTowerGeomContainer>(topNode, towergeonodename);
0191       if (!towergeo)
0192       {
0193         std::cout << PHWHERE << " ERROR: Can't find node " << towergeonodename << std::endl;
0194         return Fun4AllReturnCodes::ABORTEVENT;
0195       }
0196 
0197       // Grab the first tower, us it to get the
0198       // location along the beamline
0199       RawTowerGeomContainer::ConstRange twr_range = towergeo->get_tower_geometries();
0200       RawTowerGeomContainer::ConstIterator twr_iter = twr_range.first;
0201       RawTowerGeom* temp_geo = twr_iter->second;
0202 
0203       // Changed by Barak on 12/10/19
0204       iter.second.second = temp_geo->get_center_z();
0205       break;
0206     }
0207     default:
0208       std::cout << "invalid state reference: " << iter.second.first << std::endl;
0209       gSystem->Exit(1);
0210     }
0211   }
0212 
0213   if (m_DoVertexingFlag)
0214   {
0215     m_RaveVertexFactory = new genfit::GFRaveVertexFactory(Verbosity(), true);
0216     // m_RaveVertexFactory->setMethod("kalman-smoothing:1"); //! kalman-smoothing:1 is the defaul method
0217     if (!m_RaveVertexFactory)
0218     {
0219       std::cout << PHWHERE << " no Vertex Finder" << std::endl;
0220       return Fun4AllReturnCodes::ABORTRUN;
0221     }
0222     m_RaveVertexFactory->setMethod(m_VertexingMethod);
0223     // m_RaveVertexFactory->setBeamspot();
0224 
0225     // m_RaveVertexFactory = new PHRaveVertexFactory(Verbosity());
0226   }
0227   return Fun4AllReturnCodes::EVENT_OK;
0228 }
0229 
0230 int PHG4TrackFastSim::End(PHCompositeNode* /*topNode*/)
0231 {
0232   if (m_DoEvtDisplayFlag && m_Fitter)
0233   {
0234     m_Fitter->displayEvent();
0235   }
0236 
0237   return Fun4AllReturnCodes::EVENT_OK;
0238 }
0239 
0240 int PHG4TrackFastSim::process_event(PHCompositeNode* /*topNode*/)
0241 {
0242   m_EventCnt++;
0243 
0244   if (Verbosity() >= 2)
0245   {
0246     std::cout << "PHG4TrackFastSim::process_event: " << m_EventCnt << ".\n";
0247   }
0248 
0249   //    if(_clustermap_out)
0250   //        _clustermap_out->empty();
0251   //    else {
0252   //        LogError("_clustermap_out not found!");
0253   //        return Fun4AllReturnCodes::ABORTRUN;
0254   //    }
0255 
0256   if (not m_SvtxTrackMapOut)
0257   {
0258     LogError("m_SvtxTrackMapOut not found!");
0259     return Fun4AllReturnCodes::ABORTRUN;
0260   }
0261 
0262   std::vector<PHGenFit::Track*> rf_tracks;
0263 
0264   PHG4VtxPoint* truthVtx = m_TruthContainer->GetPrimaryVtx(m_TruthContainer->GetPrimaryVertexIndex());
0265   TVector3 vtxPoint(truthVtx->get_x(), truthVtx->get_y(), truthVtx->get_z());
0266   // Smear the vertex ONCE for all particles in the event
0267   vtxPoint.SetX(vtxPoint.x() + gsl_ran_gaussian(m_RandomGenerator, m_VertexXYResolution));
0268   vtxPoint.SetY(vtxPoint.y() + gsl_ran_gaussian(m_RandomGenerator, m_VertexXYResolution));
0269   vtxPoint.SetZ(vtxPoint.z() + gsl_ran_gaussian(m_RandomGenerator, m_VertexZResolution));
0270 
0271   PHG4TruthInfoContainer::ConstRange itr_range;
0272   if (m_PrimaryTrackingFlag)
0273   {
0274     // Tracking for primaries only
0275     itr_range = m_TruthContainer->GetPrimaryParticleRange();
0276   }
0277   else
0278   {
0279     // Check ALL particles
0280     itr_range = m_TruthContainer->GetParticleRange();
0281   }
0282 
0283   GenFitTrackMap gf_track_map;
0284   // Now we can loop over the particles
0285 
0286   for (PHG4TruthInfoContainer::ConstIterator itr = itr_range.first;
0287        itr != itr_range.second; ++itr)
0288   {
0289     PHG4Particle* particle = itr->second;
0290 
0291     TVector3 seed_pos(vtxPoint.x(), vtxPoint.y(), vtxPoint.z());
0292     TVector3 seed_mom(0, 0, 0);
0293     TMatrixDSym seed_cov(6);
0294 
0295     //! Create measurements
0296     std::vector<PHGenFit::Measurement*> measurements;
0297 
0298     PHGenFit::Measurement* vtx_meas = nullptr;
0299 
0300     if (m_UseVertexInFittingFlag)
0301     {
0302       vtx_meas = VertexMeasurement(TVector3(vtxPoint.x(),
0303                                             vtxPoint.y(),
0304                                             vtxPoint.z()),
0305                                    m_VertexXYResolution,
0306                                    m_VertexZResolution);
0307       measurements.push_back(vtx_meas);
0308     }
0309 
0310     std::unique_ptr<SvtxTrack> svtx_track_out(new SvtxTrack_FastSim_v1());
0311     PseudoPatternRecognition(particle,
0312                              measurements,
0313                              svtx_track_out.get(),
0314                              seed_pos, seed_mom,
0315                              seed_cov);
0316 
0317     if (measurements.size() < 3)
0318     {
0319       if (Verbosity() >= 2)
0320       {
0321         // LogWarning("measurements.size() < 3");
0322         std::cout << "event: " << m_EventCnt << " : measurements.size() < 3"
0323                   << "\n";
0324       }
0325       // Delete the measurements
0326       // We need to also delete the underlying genfit::AbsMeasurement object
0327       for (auto& measurement : measurements)
0328       {
0329         delete measurement->getMeasurement();
0330         delete measurement;
0331       }
0332       continue;
0333     }
0334 
0335     //! Build TrackRep from particle assumption
0336     /*!
0337      * mu+: -13
0338      * mu-: 13
0339      * pi+: 211
0340      * pi-: -211
0341      * e-:  11
0342      * e+:  -11
0343      */
0344     // int pid = 13; //
0345     // SMART(genfit::AbsTrackRep) rep = NEW(genfit::RKTrackRep)(pid);
0346     genfit::AbsTrackRep* rep = new genfit::RKTrackRep(m_PrimaryAssumptionPid);
0347 
0348     // rep->setDebugLvl(1); //DEBUG
0349 
0350     //! Initialize track with seed from pattern recognition
0351 
0352     PHGenFit::Track* track = new PHGenFit::Track(rep, seed_pos, seed_mom, seed_cov);
0353 
0354     // NOTE: We need to keep a list of tracks so they can
0355     // all be cleanly deleted at the end
0356     rf_tracks.push_back(track);
0357 
0358     // LogDEBUG;
0359     //! Add measurements to track
0360     track->addMeasurements(measurements);
0361 
0362     // LogDEBUG;
0363     //! Fit the track
0364     int fitting_err = m_Fitter->processTrack(track, false);
0365 
0366     if (fitting_err != 0)
0367     {
0368       if (Verbosity() >= 2)
0369       {
0370         // LogWarning("measurements.size() < 3");
0371         std::cout << "event: " << m_EventCnt
0372                   << " : fitting_err != 0, next track."
0373                   << "\n";
0374       }
0375       continue;
0376     }
0377 
0378     TVector3 vtx(vtxPoint.x(), vtxPoint.y(), vtxPoint.z());
0379     bool track_made = MakeSvtxTrack(svtx_track_out.get(), track,
0380                                     particle->get_track_id(),
0381                                     measurements.size(), vtx);
0382     if (Verbosity() > 1)
0383     {
0384       svtx_track_out->identify();
0385     }
0386 
0387     if (track_made)
0388     {
0389       //      track -> output container
0390 
0391       const unsigned int track_id = m_SvtxTrackMapOut->insert(svtx_track_out.get())->get_id();
0392       gf_track_map.insert({track->getGenFitTrack(), track_id});
0393     }
0394 
0395   }  // Loop all primary particles
0396 
0397   // vertex finding
0398   if (m_DoVertexingFlag)
0399   {
0400     if (!m_RaveVertexFactory)
0401     {
0402       std::cout << __PRETTY_FUNCTION__ << "Failed to init vertex finder" << std::endl;
0403       return Fun4AllReturnCodes::ABORTRUN;
0404     }
0405     if (!m_SvtxVertexMap)
0406     {
0407       std::cout << __PRETTY_FUNCTION__ << "Failed to init vertex map" << std::endl;
0408       return Fun4AllReturnCodes::ABORTRUN;
0409     }
0410 
0411     //    genfit::GFRaveVertexFactory* m_RaveVertexFactory = new genfit::GFRaveVertexFactory(10, true);
0412     //    m_RaveVertexFactory->setMethod("kalman-smoothing:1");
0413     //    m_RaveVertexFactory->setBeamspot();
0414 
0415     std::vector<genfit::GFRaveVertex*> rave_vertices;
0416     if (rf_tracks.size() >= 2)
0417     {
0418       try
0419       {
0420         std::vector<genfit::Track*> rf_gf_tracks;
0421         for (auto& rf_track : rf_tracks)
0422         {
0423           genfit::Track* track = rf_track->getGenFitTrack();
0424 
0425           if (Verbosity())
0426           {
0427             TVector3 pos, mom;
0428             TMatrixDSym cov;
0429 
0430             track->getFittedState().getPosMomCov(pos, mom, cov);
0431 
0432             std::cout << "Track getCharge = " << track->getFitStatus()->getCharge() << " getChi2 = " << track->getFitStatus()->getChi2() << " getNdf = " << track->getFitStatus()->getNdf() << std::endl;
0433             pos.Print();
0434             mom.Print();
0435             cov.Print();
0436           }
0437           if (track->getFitStatus()->getNdf() > m_VertexMinNdf)
0438           {
0439             rf_gf_tracks.push_back(track);
0440           }
0441         }
0442         m_RaveVertexFactory->findVertices(&rave_vertices, rf_gf_tracks);
0443       }
0444       catch (...)
0445       {
0446         if (Verbosity() > 1)
0447         {
0448           std::cout << PHWHERE << "GFRaveVertexFactory::findVertices failed!";
0449         }
0450       }
0451     }
0452 
0453     if (Verbosity())
0454     {
0455       std::cout << __PRETTY_FUNCTION__ << __LINE__ << " rf_tracks = " << rf_tracks.size() << " rave_vertices = " << rave_vertices.size() << std::endl;
0456     }
0457     FillSvtxVertexMap(rave_vertices, gf_track_map);
0458   }
0459 
0460   //! add tracks to event display
0461   if (m_DoEvtDisplayFlag)
0462   {
0463     std::vector<genfit::Track*> rf_gf_tracks;
0464     rf_gf_tracks.reserve(rf_tracks.size());
0465     for (auto& rf_track : rf_tracks)
0466     {
0467       rf_gf_tracks.push_back(rf_track->getGenFitTrack());
0468     }
0469     m_Fitter->getEventDisplay()->addEvent(rf_gf_tracks);
0470   }
0471   else
0472   {
0473     for (auto& rf_track : rf_tracks)
0474     {
0475       delete rf_track;
0476     }
0477     rf_tracks.clear();
0478   }
0479 
0480   //    if(m_SvtxTrackMapOut->get(0)) {
0481   //        m_SvtxTrackMapOut->get(0)->identify();
0482   //        std::cout<<"DEBUG : "<< m_SvtxTrackMapOut->get(0)->get_px() <<"\n";
0483   //        std::cout<<"DEBUG : "<< m_SvtxTrackMapOut->get(0)->get_truth_track_id() <<"\n";
0484   //    }
0485 
0486   return Fun4AllReturnCodes::EVENT_OK;
0487 }
0488 
0489 /*
0490  * Fill SvtxVertexMap from GFRaveVertexes and Tracks
0491  */
0492 bool PHG4TrackFastSim::FillSvtxVertexMap(
0493     const std::vector<genfit::GFRaveVertex*>& rave_vertices,
0494     const GenFitTrackMap& gf_track_map)
0495 {
0496   for (genfit::GFRaveVertex* rave_vtx : rave_vertices)
0497   {
0498     if (!rave_vtx)
0499     {
0500       std::cout << PHWHERE << std::endl;
0501       return false;
0502     }
0503 
0504     std::shared_ptr<SvtxVertex> svtx_vtx(new SvtxVertex_v1());
0505 
0506     svtx_vtx->set_chisq(rave_vtx->getChi2());
0507     svtx_vtx->set_ndof(rave_vtx->getNdf());
0508     svtx_vtx->set_position(0, rave_vtx->getPos().X());
0509     svtx_vtx->set_position(1, rave_vtx->getPos().Y());
0510     svtx_vtx->set_position(2, rave_vtx->getPos().Z());
0511 
0512     for (int i = 0; i < 3; i++)
0513     {
0514       for (int j = 0; j < 3; j++)
0515       {
0516         svtx_vtx->set_error(i, j, rave_vtx->getCov()[i][j]);
0517       }
0518     }
0519 
0520     for (unsigned int i = 0; i < rave_vtx->getNTracks(); i++)
0521     {
0522       // TODO improve speed
0523       const genfit::Track* rave_track =
0524           rave_vtx->getParameters(i)->getTrack();
0525       //      for(auto iter : gf_track_map) {
0526       //        if (iter.second == rave_track)
0527       //          svtx_vtx->insert_track(iter.first);
0528       //      }
0529       auto iter = gf_track_map.find(rave_track);
0530       if (iter != gf_track_map.end())
0531       {
0532         svtx_vtx->insert_track(iter->second);
0533       }
0534     }
0535 
0536     if (m_SvtxVertexMap)
0537     {
0538       m_SvtxVertexMap->insert_clone(svtx_vtx.get());
0539     }
0540     else
0541     {
0542       LogError("!m_SvtxVertexMap");
0543     }
0544 
0545   }  // loop over RAVE vertices
0546 
0547   return true;
0548 }
0549 
0550 int PHG4TrackFastSim::CreateNodes(PHCompositeNode* topNode)
0551 {
0552   // create nodes...
0553   PHNodeIterator iter(topNode);
0554 
0555   PHCompositeNode* dstNode = static_cast<PHCompositeNode*>(iter.findFirst("PHCompositeNode", "DST"));
0556   if (!dstNode)
0557   {
0558     std::cout << PHWHERE << " DST Node missing, doing nothing." << std::endl;
0559     return Fun4AllReturnCodes::ABORTEVENT;
0560   }
0561   PHNodeIterator iter_dst(dstNode);
0562 
0563   // Create the FGEM node
0564   PHCompositeNode* tb_node = dynamic_cast<PHCompositeNode*>(iter_dst.findFirst(
0565       "PHCompositeNode", m_SubTopnodeName));
0566   if (!tb_node)
0567   {
0568     tb_node = new PHCompositeNode(m_SubTopnodeName);
0569     dstNode->addNode(tb_node);
0570     if (Verbosity() > 0)
0571     {
0572       std::cout << m_SubTopnodeName << " node added" << std::endl;
0573     }
0574   }
0575 
0576   //    _clustermap_out = new SvtxClusterMap_v1;
0577   //
0578   //    PHIODataNode<PHObject>* clusters_node = new PHIODataNode<PHObject>(
0579   //            _clustermap_out, _clustermap_out_name, "PHObject");
0580   //    tb_node->addNode(clusters_node);
0581   //    if (Verbosity() > 0)
0582   //        std::cout << _clustermap_out_name <<" node added" << std::endl;
0583 
0584   m_SvtxTrackMapOut = findNode::getClass<SvtxTrackMap>(topNode, m_TrackmapOutNodeName);
0585   if (!m_SvtxTrackMapOut)
0586   {
0587     m_SvtxTrackMapOut = new SvtxTrackMap_v1;
0588 
0589     PHIODataNode<PHObject>* tracks_node = new PHIODataNode<PHObject>(m_SvtxTrackMapOut, m_TrackmapOutNodeName, "PHObject");
0590     tb_node->addNode(tracks_node);
0591     if (Verbosity() > 0)
0592     {
0593       std::cout << m_TrackmapOutNodeName << " node added" << std::endl;
0594     }
0595   }
0596 
0597   m_SvtxVertexMap = findNode::getClass<SvtxVertexMap_v1>(topNode, "SvtxVertexMap");
0598   if (not m_SvtxVertexMap)
0599   {
0600     m_SvtxVertexMap = new SvtxVertexMap_v1;
0601     PHIODataNode<PHObject>* vertexes_node = new PHIODataNode<PHObject>(m_SvtxVertexMap, "SvtxVertexMap", "PHObject");
0602     tb_node->addNode(vertexes_node);
0603     if (Verbosity() > 0)
0604     {
0605       std::cout << "Svtx/SvtxVertexMap node added" << std::endl;
0606     }
0607   }
0608 
0609   return Fun4AllReturnCodes::EVENT_OK;
0610 }
0611 
0612 int PHG4TrackFastSim::GetNodes(PHCompositeNode* topNode)
0613 {
0614   assert(m_Parameter);
0615   PHNodeIterator iter(topNode);
0616 
0617   // DST objects
0618   // Truth container
0619   m_TruthContainer = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
0620   if (!m_TruthContainer)
0621   {
0622     std::cout << PHWHERE << " PHG4TruthInfoContainer node not found on node tree"
0623          << std::endl;
0624     return Fun4AllReturnCodes::ABORTEVENT;
0625   }
0626 
0627   for (auto& m_PHG4HitsName : m_PHG4HitsNames)
0628   {
0629     PHG4HitContainer* phg4hit = findNode::getClass<PHG4HitContainer>(topNode, m_PHG4HitsName);
0630     if (!phg4hit)
0631     {
0632       std::cout << PHWHERE << m_PHG4HitsName
0633            << " node not found on node tree" << std::endl;
0634       return Fun4AllReturnCodes::ABORTEVENT;
0635     }
0636 
0637     if (Verbosity() > 0)
0638     {
0639       std::cout << "PHG4TrackFastSim::GetNodes - node added: " << m_PHG4HitsName << std::endl;
0640     }
0641     m_PHG4HitContainer.push_back(phg4hit);
0642 
0643     m_Parameter->set_int_param(m_PHG4HitsName, phg4hit->GetID());
0644   }
0645 
0646   // Update Kalman filter parameter node
0647   m_Parameter->set_string_param("SubTopnodeName", m_SubTopnodeName);
0648   m_Parameter->set_string_param("TrackmapOutNodeName", m_TrackmapOutNodeName);
0649   m_Parameter->set_string_param("VertexingMethod", m_VertexingMethod);
0650   m_Parameter->set_string_param("FitAlgoName", m_FitAlgoName);
0651 
0652   PHCompositeNode* runNode = static_cast<PHCompositeNode*>(iter.findFirst(
0653       "PHCompositeNode", "RUN"));
0654   if (!runNode)
0655   {
0656     std::cout << Name() << "::"
0657               << "::" << __PRETTY_FUNCTION__
0658               << "Run Node missing!" << std::endl;
0659     throw std::runtime_error(
0660         "Failed to find Run node in PHG4TrackFastSim::CreateNodes");
0661   }
0662   if (Verbosity())
0663   {
0664     std::cout << __PRETTY_FUNCTION__ << " : ";
0665     m_Parameter->identify();
0666   }
0667   m_Parameter->SaveToNodeTree(runNode, Name() + std::string("_Parameter"));
0668 
0669   // checks
0670   assert(m_PHG4HitsNames.size() == m_PHG4HitContainer.size());
0671   assert(m_phg4_detector_type.size() == m_PHG4HitContainer.size());
0672   assert(m_phg4_detector_radres.size() == m_PHG4HitContainer.size());
0673   assert(m_phg4_detector_phires.size() == m_PHG4HitContainer.size());
0674   assert(m_phg4_detector_lonres.size() == m_PHG4HitContainer.size());
0675   assert(m_phg4_detector_hitfindeff.size() == m_PHG4HitContainer.size());
0676   assert(m_phg4_detector_noise.size() == m_PHG4HitContainer.size());
0677 
0678   //    _clustermap_out = findNode::getClass<SvtxClusterMap>(topNode,
0679   //            _clustermap_out_name);
0680   //    if (!_clustermap_out && m_EventCnt < 2) {
0681   //        std::cout << PHWHERE << _clustermap_out_name << " node not found on node tree"
0682   //                << std::endl;
0683   //        return Fun4AllReturnCodes::ABORTEVENT;
0684   //    }
0685 
0686   m_SvtxTrackMapOut = findNode::getClass<SvtxTrackMap>(topNode, m_TrackmapOutNodeName);
0687   if (!m_SvtxTrackMapOut && m_EventCnt < 2)
0688   {
0689     std::cout << PHWHERE << m_TrackmapOutNodeName
0690          << " node not found on node tree" << std::endl;
0691     return Fun4AllReturnCodes::ABORTEVENT;
0692   }
0693 
0694   return Fun4AllReturnCodes::EVENT_OK;
0695 }
0696 
0697 int PHG4TrackFastSim::PseudoPatternRecognition(const PHG4Particle* particle,
0698                                                std::vector<PHGenFit::Measurement*>& meas_out,
0699                                                SvtxTrack* track_out,
0700                                                TVector3& seed_pos,
0701                                                TVector3& seed_mom, TMatrixDSym& seed_cov, const bool do_smearing)
0702 {
0703   assert(track_out);
0704 
0705   seed_cov.ResizeTo(6, 6);
0706 
0707   seed_pos.SetXYZ(0, 0, 0);
0708   // reset the seed resolution to the approximate position resolution of the last detector
0709   seed_cov[0][0] = .1 * .1;
0710   seed_cov[1][1] = .1 * .1;
0711   seed_cov[2][2] = 30 * 30;
0712   //  for (int i = 0; i < 3; i++)
0713   //  {
0714   //    seed_cov[i][i] = _phi_resolution * _phi_resolution;
0715   //  }
0716 
0717   seed_mom.SetXYZ(0, 0, 10);
0718   for (int i = 3; i < 6; i++)
0719   {
0720     seed_cov[i][i] = 10;
0721   }
0722 
0723   if (particle)
0724   {
0725     TVector3 True_mom(particle->get_px(), particle->get_py(),
0726                       particle->get_pz());
0727 
0728     seed_mom.SetXYZ(particle->get_px(), particle->get_py(),
0729                     particle->get_pz());
0730     if (do_smearing)
0731     {
0732       const double momSmear = 3. / 180. * M_PI;  // rad
0733       const double momMagSmear = 0.1;            // relative
0734 
0735       seed_mom.SetMag(
0736           True_mom.Mag() + gsl_ran_gaussian(m_RandomGenerator,
0737                                             momMagSmear * True_mom.Mag()));
0738       seed_mom.SetTheta(True_mom.Theta() + gsl_ran_gaussian(m_RandomGenerator, momSmear));
0739       seed_mom.SetPhi(True_mom.Phi() + gsl_ran_gaussian(m_RandomGenerator, momSmear));
0740     }
0741   }
0742 
0743   if (Verbosity())
0744   {
0745     std::cout << "PHG4TrackFastSim::PseudoPatternRecognition - DEBUG: "
0746               << "searching for hits from  " << m_PHG4HitContainer.size() << " PHG4Hit nodes" << std::endl;
0747   }
0748 
0749   // order measurement with g4hit time via stl multimap
0750   std::multimap<double, PHGenFit::Measurement*> ordered_measurements;
0751 
0752   for (unsigned int ilayer = 0; ilayer < m_PHG4HitContainer.size(); ilayer++)
0753   {
0754     if (!m_PHG4HitContainer[ilayer])
0755     {
0756       LogError("No m_PHG4HitContainer[i] found!");
0757       continue;
0758     }
0759 
0760     int dettype = m_phg4_detector_type[ilayer];
0761     float detradres = m_phg4_detector_radres[ilayer];
0762     float detphires = m_phg4_detector_phires[ilayer];
0763     float detlonres = m_phg4_detector_lonres[ilayer];
0764     float dethiteff = m_phg4_detector_hitfindeff[ilayer];
0765     float detnoise = m_phg4_detector_noise[ilayer];
0766     if (Verbosity())
0767     {
0768       std::cout << "PHG4TrackFastSim::PseudoPatternRecognition - DEBUG: "
0769                 << "ilayer: "
0770                 << ilayer << ",  " << m_PHG4HitsNames[ilayer]
0771                 << " with nsublayers: " << m_PHG4HitContainer[ilayer]->num_layers()
0772                 << ", detradres = " << detradres
0773                 << ", detphires = " << detphires
0774                 << ", detlonres = " << detlonres
0775                 << ", dethiteff = " << dethiteff
0776                 << ", detnoise = " << detnoise
0777                 << " \n";
0778     }
0779     for (PHG4HitContainer::LayerIter layerit =
0780              m_PHG4HitContainer[ilayer]->getLayers().first;
0781          layerit != m_PHG4HitContainer[ilayer]->getLayers().second; layerit++)
0782     {
0783       for (PHG4HitContainer::ConstIterator itr =
0784                m_PHG4HitContainer[ilayer]->getHits(*layerit).first;
0785            itr != m_PHG4HitContainer[ilayer]->getHits(*layerit).second; ++itr)
0786       {
0787         PHG4Hit* hit = itr->second;
0788         if (!hit)
0789         {
0790           LogDebug("No PHG4Hit Found!");
0791           continue;
0792         }
0793 
0794         if (hit->get_trkid() == particle->get_track_id() || gsl_ran_binomial(m_RandomGenerator, detnoise, 1) > 0)
0795         {
0796           if (gsl_ran_binomial(m_RandomGenerator, dethiteff, 1) > 0)
0797           {
0798             PHGenFit::Measurement* meas = nullptr;
0799             if (dettype == Vertical_Plane)
0800             {
0801               if (Verbosity())
0802               {
0803                 std::cout << "PHG4TrackFastSim::PseudoPatternRecognition -adding vertical plane hit ilayer: "
0804                           << ilayer << "; detphires: " << detphires << "; detradres: " << detradres << " \n";
0805                 hit->identify();
0806               }
0807               meas = PHG4HitToMeasurementVerticalPlane(hit,
0808                                                        detphires, detradres);
0809             }
0810             else if (dettype == Cylinder)
0811             {
0812               if (Verbosity())
0813               {
0814                 std::cout << "PHG4TrackFastSim::PseudoPatternRecognition -adding cylinder hit ilayer: "
0815                           << ilayer << "; detphires: " << detphires << "; detlonres : " << detlonres << " \n";
0816                 hit->identify();
0817               }
0818               meas = PHG4HitToMeasurementCylinder(hit,
0819                                                   detphires, detlonres);
0820             }
0821             else
0822             {
0823               LogError("Type not implemented!");
0824               return Fun4AllReturnCodes::ABORTEVENT;
0825             }
0826             //            meas_out.push_back(meas);
0827             ordered_measurements.insert(std::make_pair(hit->get_avg_t(), meas));
0828             track_out->add_g4hit_id(m_PHG4HitContainer[ilayer]->GetID(), hit->get_hit_id());
0829             // meas->getMeasurement()->Print(); //DEBUG
0830           }
0831         }
0832       }
0833     } /*Loop layers within one detector layer*/
0834   }   /*Loop detector layers*/
0835 
0836   for (auto& pair : ordered_measurements)
0837   {
0838     meas_out.push_back(pair.second);
0839 
0840     if (Verbosity())
0841     {
0842       std::cout << "PHG4TrackFastSim::PseudoPatternRecognition - measruement at t =  " << pair.first << " ns: ";
0843       pair.second->getMeasurement()->Print();
0844     }
0845   }
0846 
0847   if (Verbosity())
0848   {
0849     std::cout << "PHG4TrackFastSim::PseudoPatternRecognition - meas_out.size = " << meas_out.size() << " for "
0850               << "particle: "
0851               << std::endl;
0852     particle->identify();
0853 
0854     std::cout << "PHG4TrackFastSim::PseudoPatternRecognition - seed_pos = "
0855               << seed_pos.x() << ", " << seed_pos.y() << ". " << seed_pos.z() << std::endl;
0856     std::cout << "PHG4TrackFastSim::PseudoPatternRecognition - seed_pos = "
0857               << seed_mom.x() << ", " << seed_mom.y() << ". " << seed_mom.z() << std::endl;
0858     std::cout << "PHG4TrackFastSim::PseudoPatternRecognition - seed_cov = "
0859               << sqrt(seed_cov[0][0]) << ", " << sqrt(seed_cov[1][1]) << ". " << sqrt(seed_cov[2][2])
0860               << ","
0861               << sqrt(seed_cov[3][3]) << ", " << sqrt(seed_cov[4][4]) << ". " << sqrt(seed_cov[5][5]) << std::endl;
0862   }
0863 
0864   return Fun4AllReturnCodes::EVENT_OK;
0865 }
0866 
0867 bool PHG4TrackFastSim::MakeSvtxTrack(SvtxTrack* out_track,
0868                                      const PHGenFit::Track* phgf_track,
0869                                      const unsigned int truth_track_id,
0870                                      const unsigned int nmeas,
0871                                      const TVector3& vtx)
0872 {
0873   assert(out_track);
0874 
0875   double chi2 = phgf_track->get_chi2();
0876   double ndf = phgf_track->get_ndf();
0877 
0878   double pathlenth_from_first_meas = -999999;
0879   std::unique_ptr<genfit::MeasuredStateOnPlane> gf_state(new genfit::MeasuredStateOnPlane());
0880 
0881   //  if (_detector_type == Vertical_Plane)
0882   //  {
0883   //    pathlenth_orig_from_first_meas = phgf_track->extrapolateToPlane(*gf_state, vtx,
0884   //                                                                    TVector3(0., 0., 1.), 0);
0885   //  }
0886   //  else if (_detector_type == Cylinder)
0887   //    pathlenth_orig_from_first_meas = phgf_track->extrapolateToLine(*gf_state, vtx,
0888   //                                                                   TVector3(0., 0., 1.));
0889   //  else
0890   //  {
0891   //    LogError("Detector Type NOT implemented!");
0892   //    return nullptr;
0893   //  }
0894 
0895   // always extrapolate to a z-line through the vertex
0896   double pathlenth_orig_from_first_meas = phgf_track->extrapolateToLine(*gf_state, vtx,
0897                                                                         TVector3(0., 0., 1.));
0898 
0899   if (Verbosity() > 1)
0900   {
0901     std::cout << __PRETTY_FUNCTION__ << __LINE__ << " pathlenth_orig_from_first_meas = " << pathlenth_orig_from_first_meas << std::endl;
0902   }
0903   if (pathlenth_orig_from_first_meas < -999990)
0904   {
0905     LogError("Extraction faild!");
0906     return false;
0907   }
0908 
0909   TVector3 mom = gf_state->getMom();
0910   TVector3 pos = gf_state->getPos();
0911   TMatrixDSym cov = gf_state->get6DCov();
0912 
0913   out_track->set_truth_track_id(truth_track_id);
0914   /*!
0915    * TODO: check the definition
0916    *  1/p, u'/z', v'/z', u, v
0917    *  u is defined as mom X beam line at POCA
0918    *  so u is the dca2d direction
0919    */
0920   double dca2d = gf_state->getState()[3];
0921   out_track->set_dca2d(dca2d);
0922   out_track->set_dca2d_error(gf_state->getCov()[3][3]);
0923   double dca3d = sqrt(dca2d * dca2d + gf_state->getState()[4] * gf_state->getState()[4]);
0924   out_track->set_dca(dca3d);
0925 
0926   out_track->set_chisq(chi2);
0927   out_track->set_ndf(ndf);
0928   out_track->set_charge(phgf_track->get_charge());
0929 
0930   out_track->set_num_measurements(nmeas);
0931 
0932   out_track->set_px(mom.Px());
0933   out_track->set_py(mom.Py());
0934   out_track->set_pz(mom.Pz());
0935 
0936   out_track->set_x(pos.X());
0937   out_track->set_y(pos.Y());
0938   out_track->set_z(pos.Z());
0939   for (int i = 0; i < 6; i++)
0940   {
0941     for (int j = i; j < 6; j++)
0942     {
0943       out_track->set_error(i, j, cov[i][j]);
0944     }
0945   }
0946   // the default name is UNKNOWN - let's set this to ORIGIN since it is at pathlength=0
0947   out_track->begin_states()->second->set_name("ORIGIN");
0948 
0949   // make the projections for all detector types
0950   for (auto& iter : m_ProjectionsMap)
0951   {
0952     switch (iter.second.first)
0953     {
0954     case DETECTOR_TYPE::Cylinder:
0955       pathlenth_from_first_meas = phgf_track->extrapolateToCylinder(*gf_state, iter.second.second, TVector3(0., 0., 0.), TVector3(0., 0., 1.), nmeas - 1);
0956       break;
0957     case DETECTOR_TYPE::Vertical_Plane:
0958       pathlenth_from_first_meas = phgf_track->extrapolateToPlane(*gf_state, TVector3(0., 0., iter.second.second), TVector3(0, 0., 1.), nmeas - 1);
0959       break;
0960     default:
0961       std::cout << "how in the world did you get here??????" << std::endl;
0962       gSystem->Exit(1);
0963     }
0964     if (pathlenth_from_first_meas < -999990)
0965     {
0966       continue;
0967     }
0968     SvtxTrackState* state = new SvtxTrackState_v1(pathlenth_from_first_meas - pathlenth_orig_from_first_meas);
0969     state->set_x(gf_state->getPos().x());
0970     state->set_y(gf_state->getPos().y());
0971     state->set_z(gf_state->getPos().z());
0972 
0973     state->set_px(gf_state->getMom().x());
0974     state->set_py(gf_state->getMom().y());
0975     state->set_pz(gf_state->getMom().z());
0976 
0977     state->set_name(iter.first);
0978     for (int i = 0; i < 6; i++)
0979     {
0980       for (int j = i; j < 6; j++)
0981       {
0982         state->set_error(i, j, gf_state->get6DCov()[i][j]);
0983       }
0984     }
0985     out_track->insert_state(state);
0986     // the state is cloned on insert_state, so delete this copy here!
0987     delete state;
0988   }
0989 
0990   return true;
0991 }
0992 
0993 PHGenFit::PlanarMeasurement* PHG4TrackFastSim::PHG4HitToMeasurementVerticalPlane(
0994     const PHG4Hit* g4hit, const double phi_resolution,
0995     const double r_resolution)
0996 {
0997   TVector3 pos(g4hit->get_avg_x(), g4hit->get_avg_y(), g4hit->get_avg_z());
0998 
0999   TVector3 v(pos.X(), pos.Y(), 0);
1000   v = 1 / v.Mag() * v;
1001 
1002   TVector3 u = v.Cross(TVector3(0, 0, 1));
1003   u = 1 / u.Mag() * u;
1004 
1005   double u_smear = 0.;
1006   double v_smear = 0.;
1007   if (m_SmearingFlag)
1008   {
1009     u_smear = gsl_ran_gaussian(m_RandomGenerator, phi_resolution);
1010     v_smear = gsl_ran_gaussian(m_RandomGenerator, r_resolution);
1011   }
1012   pos.SetX(g4hit->get_avg_x() + u_smear * u.X() + v_smear * v.X());
1013   pos.SetY(g4hit->get_avg_y() + u_smear * u.Y() + v_smear * v.Y());
1014 
1015   PHGenFit::PlanarMeasurement* meas = new PHGenFit::PlanarMeasurement(pos, u, v, phi_resolution,
1016                                                                       r_resolution);
1017 
1018   //    std::cout<<"------------\n";
1019   //    pos.Print();
1020   //    std::cout<<"at "<<istation<<" station, "<<ioctant << " octant \n";
1021   //    u.Print();
1022   //    v.Print();
1023 
1024   // dynamic_cast<PHGenFit::PlanarMeasurement*> (meas)->getMeasurement()->Print();
1025 
1026   return meas;
1027 }
1028 
1029 PHGenFit::PlanarMeasurement* PHG4TrackFastSim::PHG4HitToMeasurementCylinder(
1030     const PHG4Hit* g4hit, const double phi_resolution,
1031     const double z_resolution)
1032 {
1033   TVector3 pos(g4hit->get_avg_x(), g4hit->get_avg_y(), g4hit->get_avg_z());
1034 
1035   TVector3 v(0, 0, 1);
1036 
1037   TVector3 u = v.Cross(TVector3(pos.X(), pos.Y(), 0));
1038   u = 1 / u.Mag() * u;
1039 
1040   double u_smear = 0.;
1041   double v_smear = 0.;
1042   if (m_SmearingFlag)
1043   {
1044     u_smear = gsl_ran_gaussian(m_RandomGenerator, phi_resolution);
1045     v_smear = gsl_ran_gaussian(m_RandomGenerator, z_resolution);
1046   }
1047   pos.SetX(g4hit->get_avg_x() + u_smear * u.X());
1048   pos.SetY(g4hit->get_avg_y() + u_smear * u.Y());
1049   pos.SetZ(g4hit->get_avg_z() + v_smear);
1050 
1051   PHGenFit::PlanarMeasurement* meas = new PHGenFit::PlanarMeasurement(pos, u, v, phi_resolution,
1052                                                                       z_resolution);
1053 
1054   //    std::cout<<"------------\n";
1055   //    pos.Print();
1056   //    std::cout<<"at "<<istation<<" station, "<<ioctant << " octant \n";
1057   //    u.Print();
1058   //    v.Print();
1059 
1060   // dynamic_cast<PHGenFit::PlanarMeasurement*> (meas)->getMeasurement()->Print();
1061 
1062   return meas;
1063 }
1064 
1065 PHGenFit::Measurement* PHG4TrackFastSim::VertexMeasurement(const TVector3& vtx, double dxy, double dz)
1066 {
1067   TMatrixDSym cov(3);
1068   cov.Zero();
1069   cov(0, 0) = dxy * dxy;
1070   cov(1, 1) = dxy * dxy;
1071   cov(2, 2) = dz * dz;
1072 
1073   TVector3 pos = vtx;
1074   pos.SetX(vtx.X());
1075   pos.SetY(vtx.Y());
1076   pos.SetZ(vtx.Z());
1077 
1078   PHGenFit::Measurement* meas = new PHGenFit::SpacepointMeasurement(pos, cov);
1079 
1080   return meas;
1081 }
1082 
1083 void PHG4TrackFastSim::DisplayEvent() const
1084 {
1085   if (m_DoEvtDisplayFlag && m_Fitter)
1086   {
1087     m_Fitter->displayEvent();
1088   }
1089   return;
1090 }
1091 
1092 void PHG4TrackFastSim::add_state_name(const std::string& stateName)
1093 {
1094   if (reserved_zplane_projection_names.find(stateName) != reserved_zplane_projection_names.end())
1095   {
1096     m_ProjectionsMap.insert(std::make_pair(stateName, std::make_pair(DETECTOR_TYPE::Vertical_Plane, NAN)));
1097   }
1098   else if (reserved_cylinder_projection_names.find(stateName) != reserved_cylinder_projection_names.end())
1099   {
1100     m_ProjectionsMap.insert(std::make_pair(stateName, std::make_pair(DETECTOR_TYPE::Cylinder, NAN)));
1101   }
1102   else
1103   {
1104     std::cout << PHWHERE << " Invalid stateName " << stateName << std::endl;
1105     std::cout << std::endl
1106          << "These are implemented for cylinders" << std::endl;
1107     for (const auto& iter : reserved_cylinder_projection_names)
1108     {
1109       std::cout << iter << std::endl;
1110     }
1111     std::cout << std::endl
1112          << "These are implemented are for zplanes" << std::endl;
1113     for (const auto& iter : reserved_zplane_projection_names)
1114     {
1115       std::cout << iter << std::endl;
1116     }
1117     gSystem->Exit(1);
1118   }
1119   return;
1120 }
1121 
1122 void PHG4TrackFastSim::add_cylinder_state(const std::string& stateName, const double radius)
1123 {
1124   if (reserved_cylinder_projection_names.find(stateName) != reserved_cylinder_projection_names.end() ||
1125       reserved_zplane_projection_names.find(stateName) != reserved_zplane_projection_names.end())
1126   {
1127     std::cout << PHWHERE << ": " << stateName << " is a reserved name, used a different name for your cylinder projection" << std::endl;
1128     gSystem->Exit(1);
1129   }
1130   if (m_ProjectionsMap.find(stateName) != m_ProjectionsMap.end())
1131   {
1132     std::cout << PHWHERE << ": " << stateName << " is already a projection, please rename" << std::endl;
1133     gSystem->Exit(1);
1134   }
1135   m_ProjectionsMap.insert(std::make_pair(stateName, std::make_pair(DETECTOR_TYPE::Cylinder, radius)));
1136   return;
1137 }
1138 
1139 void PHG4TrackFastSim::add_zplane_state(const std::string& stateName, const double zplane)
1140 {
1141   if (reserved_cylinder_projection_names.find(stateName) != reserved_cylinder_projection_names.end() ||
1142       reserved_zplane_projection_names.find(stateName) != reserved_zplane_projection_names.end())
1143   {
1144     std::cout << PHWHERE << ": " << stateName << " is  a reserved name, used different name for your zplane projection" << std::endl;
1145     gSystem->Exit(1);
1146   }
1147   if (m_ProjectionsMap.find(stateName) != m_ProjectionsMap.end())
1148   {
1149     std::cout << PHWHERE << ": " << stateName << " is already a projection, please rename" << std::endl;
1150     gSystem->Exit(1);
1151   }
1152   m_ProjectionsMap.insert(std::make_pair(stateName, std::make_pair(DETECTOR_TYPE::Vertical_Plane, zplane)));
1153   return;
1154 }
[ 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!