Back to home page

sPhenix code displayed by LXR

 
 

    


File indexing completed on 2025-12-18 09:20:45

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