File indexing completed on 2025-12-18 09:20:45
0001
0002
0003
0004
0005
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 }
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
0097
0098
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")
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();
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
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
0201
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
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
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
0227
0228
0229 }
0230 return Fun4AllReturnCodes::EVENT_OK;
0231 }
0232
0233 int PHG4TrackFastSim::End(PHCompositeNode* )
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* )
0244 {
0245 m_EventCnt++;
0246
0247 if (Verbosity() >= 2)
0248 {
0249 std::cout << "PHG4TrackFastSim::process_event: " << m_EventCnt << ".\n";
0250 }
0251
0252
0253
0254
0255
0256
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
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
0278 itr_range = m_TruthContainer->GetPrimaryParticleRange();
0279 }
0280 else
0281 {
0282
0283 itr_range = m_TruthContainer->GetParticleRange();
0284 }
0285
0286 GenFitTrackMap gf_track_map;
0287
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
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
0325 std::cout << "event: " << m_EventCnt << " : measurements.size() < 3"
0326 << "\n";
0327 }
0328
0329
0330 for (auto& measurement : measurements)
0331 {
0332 delete measurement->getMeasurement();
0333 delete measurement;
0334 }
0335 continue;
0336 }
0337
0338
0339
0340
0341
0342
0343
0344
0345
0346
0347
0348
0349 genfit::AbsTrackRep* rep = new genfit::RKTrackRep(m_PrimaryAssumptionPid);
0350
0351
0352
0353
0354
0355 PHGenFit::Track* track = new PHGenFit::Track(rep, seed_pos, seed_mom, seed_cov);
0356
0357
0358
0359 rf_tracks.push_back(track);
0360
0361
0362
0363 track->addMeasurements(measurements);
0364
0365
0366
0367 int fitting_err = m_Fitter->processTrack(track, false);
0368
0369 if (fitting_err != 0)
0370 {
0371 if (Verbosity() >= 2)
0372 {
0373
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
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 }
0399
0400
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
0415
0416
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
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
0485
0486
0487
0488
0489
0490 return Fun4AllReturnCodes::EVENT_OK;
0491 }
0492
0493
0494
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
0527 const genfit::Track* rave_track =
0528 rave_vtx->getParameters(i)->getTrack();
0529
0530
0531
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 }
0550
0551 return true;
0552 }
0553
0554 int PHG4TrackFastSim::CreateNodes(PHCompositeNode* topNode)
0555 {
0556
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
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
0581
0582
0583
0584
0585
0586
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
0622
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
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
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
0683
0684
0685
0686
0687
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
0713 seed_cov[0][0] = .1 * .1;
0714 seed_cov[1][1] = .1 * .1;
0715 seed_cov[2][2] = 30 * 30;
0716
0717
0718
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;
0737 const double momMagSmear = 0.1;
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
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
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
0834 }
0835 }
0836 }
0837 }
0838 }
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
0886
0887
0888
0889
0890
0891
0892
0893
0894
0895
0896
0897
0898
0899
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
0920
0921
0922
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
0951 out_track->begin_states()->second->set_name("ORIGIN");
0952
0953
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
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
1023
1024
1025
1026
1027
1028
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
1059
1060
1061
1062
1063
1064
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 }