From 7786bd93465b27651e98da85e7080a7c7f8eeb30 Mon Sep 17 00:00:00 2001 From: Daniel Brandenburg Date: Sun, 17 Dec 2023 20:21:22 -0500 Subject: [PATCH] FST tracking mode for high mult and data with missing sTGCs --- StRoot/StFwdTrackMaker/StFwdTrackMaker.cxx | 163 ++-- StRoot/StFwdTrackMaker/StFwdTrackMaker.h | 28 +- .../StFwdTrackMaker/include/Tracker/FwdHit.h | 8 +- .../include/Tracker/FwdTracker.h | 809 ++++++++++++------ .../include/Tracker/QualityPlotter.h | 12 +- .../include/Tracker/TrackFitter.h | 263 ++++-- 6 files changed, 853 insertions(+), 430 deletions(-) diff --git a/StRoot/StFwdTrackMaker/StFwdTrackMaker.cxx b/StRoot/StFwdTrackMaker/StFwdTrackMaker.cxx index 46d9313eb16..898660c4ad6 100644 --- a/StRoot/StFwdTrackMaker/StFwdTrackMaker.cxx +++ b/StRoot/StFwdTrackMaker/StFwdTrackMaker.cxx @@ -86,7 +86,6 @@ class GenfitUtils{ }; // GenfitUtils - // Basic sanity cuts on genfit tracks template<> bool GenfitUtils::accept( genfit::Track *track ) { @@ -148,9 +147,7 @@ template<> bool GenfitUtils::accept( genfit::Track *track ) }; - //______________________________________________________________________________________ - class SiRasterizer { public: SiRasterizer() {} @@ -159,7 +156,7 @@ class SiRasterizer { void setup(FwdTrackerConfig &_cfg) { cfg = _cfg; mRasterR = cfg.get("SiRasterizer:r", 3.0); - mRasterPhi = cfg.get("SiRasterizer:phi", 0.1); + mRasterPhi = cfg.get("SiRasterizer:phi", 0.004); } bool active() { @@ -226,8 +223,6 @@ class ForwardTracker : public ForwardTrackMaker { } }; - - //________________________________________________________________________ StFwdTrackMaker::StFwdTrackMaker() : StMaker("fwdTrack"), mGenHistograms(false), mGenTree(false), mForwardTracker(nullptr), mForwardData(nullptr){ SetAttr("useFtt",1); // Default Ftt on @@ -303,6 +298,7 @@ int StFwdTrackMaker::Init() { mTree->Branch("fcsX", &mTreeData. fcsX ); mTree->Branch("fcsY", &mTreeData. fcsY ); mTree->Branch("fcsZ", &mTreeData. fcsZ ); + mTree->Branch("fcsE", &mTreeData. fcsE ); mTree->Branch("fcsDet", &mTreeData. fcsDet ); // mc tracks @@ -312,6 +308,8 @@ int StFwdTrackMaker::Init() { mTree->Branch("mcPhi", &mTreeData. mcPhi ); mTree->Branch("mcCharge", &mTreeData. mcCharge ); mTree->Branch("mcVertexId", &mTreeData. mcVertexId ); + mTree->Branch("mcNumFtt", &mTreeData. mcNumFtt ); + mTree->Branch("mcNumFst", &mTreeData. mcNumFst ); // mcverts mTree->Branch("vmcN", &mTreeData. vmcN, "vmcN/I"); @@ -343,6 +341,8 @@ int StFwdTrackMaker::Init() { mTree->Branch("thaX", &mTreeData. thaX ); mTree->Branch("thaY", &mTreeData. thaY ); mTree->Branch("thaZ", &mTreeData. thaZ ); + mTree->Branch("thdP", &mTreeData. thdP ); + mTree->Branch("thdR", &mTreeData. thdR ); // track projections mTree->Branch("tprojN", &mTreeData. tprojN, "tprojN/I"); @@ -581,8 +581,6 @@ void StFwdTrackMaker::loadFttHits( FwdDataSource::McTrackMap_t &mcTrackMap, FwdD } } // loadFttHits - - void StFwdTrackMaker::loadFttHitsFromStEvent( FwdDataSource::McTrackMap_t &mcTrackMap, FwdDataSource::HitMap_t &hitMap, int count ){ LOG_DEBUG << "Loading FTT Hits from Data" << endm; StEvent *event = (StEvent *)GetDataSet("StEvent"); @@ -720,6 +718,7 @@ void StFwdTrackMaker::loadFttHitsFromGEANT( FwdDataSource::McTrackMap_t &mcTrack // Add hit pointer to the track if (mcTrackMap[track_id]){ mcTrackMap[track_id]->addHit(hit); + mTreeData.mcNumFtt[track_id-1] ++; } else { LOG_ERROR << "Cannot find MC track for GEANT hit (FTT), track_id = " << track_id << endm; } @@ -766,9 +765,10 @@ void StFwdTrackMaker::loadFstHits( FwdDataSource::McTrackMap_t &mcTrackMap, FwdD LOG_DEBUG << "FST HIT: d = " << d << ", x=" << x0 << ", y=" << y0 << ", z=" << vZ << endm; mFstHits.push_back( TVector3( x0, y0, vZ) ); - FwdHit *hit = new FwdHit(count++, x0, y0, vZ, d, 0, hitCov3, nullptr); + // we use d+4 so that both FTT and FST start at 4 + FwdHit *hit = new FwdHit(count++, x0, y0, vZ, d+4, 0, hitCov3, nullptr); // Add the hit to the hit map - hitMap[hit->getSector()].push_back(hit); + hitMap[d+4].push_back(hit); mTreeData.fstX.push_back( x0 ); mTreeData.fstY.push_back( y0 ); @@ -919,6 +919,16 @@ void StFwdTrackMaker::loadFstHitsFromGEANT( FwdDataSource::McTrackMap_t &mcTrack // Add the hit to the hit map hitMap[hit->getSector()].push_back(hit); + + // Add hit pointer to the track + if (mcTrackMap[track_id]){ + mcTrackMap[track_id]->addFstHit(hit); + mTreeData.mcNumFst[track_id-1] ++; + LOG_DEBUG << TString::Format("Adding FST Hit to McTrack, mcPt=%f", mTreeData.mcPt[track_id-1] ) << endm; + + } else { + LOG_ERROR << "Cannot find MC track for GEANT hit (FTT), track_id = " << track_id << endm; + } } } // loadFstHitsFromGEANT @@ -968,6 +978,8 @@ size_t StFwdTrackMaker::loadMcTracks( FwdDataSource::McTrackMap_t &mcTrackMap ){ mTreeData.mcPhi.push_back( phi ); mTreeData.mcCharge.push_back( q ); mTreeData.mcVertexId.push_back( track->start_vertex_p ); + mTreeData.mcNumFtt.push_back(0); + mTreeData.mcNumFst.push_back(0); if (track->is_shower) nShowers++; @@ -1043,9 +1055,10 @@ void StFwdTrackMaker::loadFcs( ) { mTreeData.fcsX.push_back( xyz.x() ); mTreeData.fcsY.push_back( xyz.y() ); mTreeData.fcsZ.push_back( xyz.z() - 2 ); + mTreeData.fcsE.push_back( clu->energy() ); mTreeData.fcsDet.push_back( idet ); - } - } + } // i + } // idet // LOAD PRESHOWER HITS (EPD) for ( int det = 4; det < 6; det ++ ) { @@ -1070,11 +1083,12 @@ void StFwdTrackMaker::loadFcs( ) { mTreeData.fcsX.push_back( x0 ); mTreeData.fcsY.push_back( y0 ); mTreeData.fcsZ.push_back( zepd ); + mTreeData.fcsE.push_back( hit->energy() ); mTreeData.fcsDet.push_back( det ); - } - } - } + } // if det + } // for i + } // for det } // loadFcs @@ -1083,6 +1097,8 @@ int StFwdTrackMaker::Make() { // START time for measuring tracking long long itStart = FwdTrackerUtils::nowNanoSecond(); + StEvent *stEvent = static_cast(GetInputDS("StEvent")); + // Access forward Tracker maps FwdDataSource::McTrackMap_t &mcTrackMap = mForwardData->getMcTracks(); FwdDataSource::HitMap_t &hitMap = mForwardData->getFttHits(); @@ -1135,19 +1151,24 @@ int StFwdTrackMaker::Make() { LOG_DEBUG << ">>START Event Forward Tracking" << endm; mForwardTracker->doEvent(); LOG_DEBUG << "< getRecoTracks().size() << " GenFit Tracks" << endm; - + LOG_DEBUG << "< getTrackSeeds().size() << " GenFit Tracks" << endm; + // perform vertex finding/fitting with Fwd tracks FitVertex(); - - StEvent *stEvent = static_cast(GetInputDS("StEvent")); /**********************************************************************/ - // Run Track finding + fitting - const auto &genfitTracks = mForwardTracker -> globalTracks(); - if ( mVisualize /* && genfitTracks.size() > 0 && genfitTracks.size() < 200*/ ) { - const auto &seed_tracks = mForwardTracker -> getRecoTracks(); + + /**********************************************************************/ + // Output track visualization if configured to do so + std::vector genfitTracks; + for ( auto gtr : mForwardTracker->getTrackResults() ) { + if ( gtr.isFitConvergedFully == false ) continue; + genfitTracks.push_back( gtr.track ); + } + + if ( mVisualize && genfitTracks.size() > 0 && genfitTracks.size() < 400 && eventIndex < 50 ) { + const auto &seed_tracks = mForwardTracker -> getTrackSeeds(); ObjExporter woe; woe.output( @@ -1159,16 +1180,16 @@ int StFwdTrackMaker::Make() { LOG_DEBUG << "Done Writing OBJ " << endm; } else if (mVisualize && genfitTracks.size() == 0) { LOG_DEBUG << "Skipping visualization, no FWD tracks" << endm; - } else if (mVisualize && genfitTracks.size() >= 20) { + } else if (mVisualize && genfitTracks.size() >= 400) { LOG_DEBUG << "Skipping visualization, too many FWD tracks" << endm; } // Fill Track Deltas in ttree for helpful alignment info FillTrackDeltas(); - LOG_INFO << "Forward tracking on this event took " << (FwdTrackerUtils::nowNanoSecond() - itStart) * 1e-6 << " ms" << endm; + LOG_DEBUG << "Forward tracking on this event took " << (FwdTrackerUtils::nowNanoSecond() - itStart) * 1e-6 << " ms" << endm; - if ( true && IAttr("fillEvent") ) { + if ( IAttr("fillEvent") ) { if (!stEvent) { LOG_WARN << "No StEvent found. Forward tracks will not be saved" << endm; @@ -1188,46 +1209,38 @@ StFwdTrack * StFwdTrackMaker::makeStFwdTrack( GenfitTrackResult >r, size_t ind LOG_DEBUG << "StFwdTrackMaker::makeStFwdTrack()" << endm; StFwdTrack *fwdTrack = new StFwdTrack( ); - auto track = gtr.track; - // if FST refit is available save that - if ( gtr.nFST > 0 && gtr.fstTrack != nullptr){ - LOG_DEBUG << "\tSave FST refit track since we have FST points" << endm; - track = gtr.fstTrack; - } else if (gtr.nFST > 0 && gtr.fstTrack == nullptr) { - LOG_WARN << "\tFST refit failed even though we have " << gtr.nFST << " FST points" << endm; - } - // Fit failed beyond use - if ( track == nullptr ){ - LOG_DEBUG << "Track is nullptr, not saving StFwdTrack" << endm; + if ( gtr.track == nullptr ){ + LOG_DEBUG << "GenfitTrack is nullptr, not making StFwdTrack" << endm; return nullptr; } - auto fitStatus = track->getFitStatus(); - if ( !fitStatus ) + auto fitStatus = gtr.track->getFitStatus(); + if ( !fitStatus ) { + LOG_DEBUG << "GenfitTrack FitStatus nullptr, not making StFwdTrack" << endm; return nullptr; - + } // Fill charge and quality info fwdTrack->setDidFitConverge( fitStatus->isFitConverged() ); fwdTrack->setDidFitConvergeFully( fitStatus->isFitConvergedFully() ); fwdTrack->setNumberOfFailedPoints( fitStatus->getNFailedPoints() ); - fwdTrack->setNumberOfFitPoints( track->getNumPoints() ); + fwdTrack->setNumberOfFitPoints( gtr.track->getNumPoints() ); fwdTrack->setChi2( fitStatus->getChi2() ); fwdTrack->setNDF( fitStatus->getNdf() ); fwdTrack->setPval( fitStatus->getPVal() ); - auto cr = track->getCardinalRep(); + auto cr = gtr.track->getCardinalRep(); // charge at first point fwdTrack->setCharge( gtr.charge ); - TVector3 p = cr->getMom( track->getFittedState( 0, cr )); + TVector3 p = cr->getMom( gtr.track->getFittedState( 0, cr )); fwdTrack->setPrimaryMomentum( StThreeVectorD( gtr.momentum.X(), gtr.momentum.Y(), gtr.momentum.Z() ) ); LOG_DEBUG << "Making StFwdTrack with " << TString::Format( "p=(%f, %f, %f)", fwdTrack->momentum().x(), fwdTrack->momentum().y(), fwdTrack->momentum().z() ) << endm; int nSeedPoints = 0; // store the seed points from FTT - for ( auto s : gtr.trackSeed ){ + for ( auto s : gtr.fttSeed ){ FwdHit * fh = static_cast( s ); if (!fh) continue; float cov[9]; @@ -1303,11 +1316,11 @@ StFwdTrack * StFwdTrackMaker::makeStFwdTrack( GenfitTrackResult >r, size_t ind float cov[9]; TVector3 tv3(0, 0, 0); - if ( detIndex != kFcsHcalId ){ - tv3 = ObjExporter::trackPosition( track, z, cov, mom ); + if ( detIndex != kFcsHcalId && detIndex != kFcsWcalId ){ + tv3 = ObjExporter::trackPosition( gtr.track, z, cov, mom ); } else { // use a straight line projection to HCAL since GenFit cannot handle long projections - tv3 = ObjExporter::projectAsStraightLine( track, 575.0, 625.0, z, cov, mom ); + tv3 = ObjExporter::projectAsStraightLine( gtr.track, 575.0, 625.0, z, cov, mom ); } fwdTrack->mProjections.push_back( StFwdTrackProjection( detIndex, StThreeVectorF( tv3.X(), tv3.Y(), tv3.Z() ), StThreeVectorF( mom.X(), mom.Y(), mom.Z() ), cov) ); @@ -1358,36 +1371,49 @@ void StFwdTrackMaker::FillEvent() { } void StFwdTrackMaker::FillTrackDeltas(){ - LOG_DEBUG << "Filling Track Deltas for Alignment" << endm; const auto &fittedTracks = mForwardTracker -> getTrackResults(); for ( size_t i = 0; i < fittedTracks.size(); i++ ){ - auto st = fittedTracks[i].trackSeed; auto gt = fittedTracks[i].track; + Seed_t combinedSeed; + combinedSeed.insert( combinedSeed.begin(), fittedTracks[i].fstSeed.begin(), fittedTracks[i].fstSeed.end() ); // this is goofed but will fix + combinedSeed.insert( combinedSeed.end(), fittedTracks[i].fttSeed.begin(), fittedTracks[i].fttSeed.end() ); + if (fittedTracks[i].isFitConvergedFully == false){ LOG_DEBUG << "Skipping track, failed fit" << endm; continue; } - for ( KiTrack::IHit * hit : st ){ + for ( KiTrack::IHit * hit : combinedSeed ){ TVector3 htv3(hit->getX(), hit->getY(), hit->getZ()); - + TLorentzVector lvHit, lvTrack; + lvHit.SetPxPyPzE( hit->getX(), hit->getY(), hit->getZ(), 1 ); + auto ttv3 = ObjExporter::trackPosition( gt, htv3.Z() ); - + lvTrack.SetPxPyPzE( ttv3.X(), ttv3.Y(), ttv3.Z(), 1 ); + mTreeData.thdX.push_back( (ttv3.X() - htv3.X()) ); mTreeData.thdY.push_back( (ttv3.Y() - htv3.Y()) ); + + mTreeData.thdP.push_back( lvTrack.DeltaPhi( lvHit ) ); + mTreeData.thdR.push_back( lvTrack.Pt() - lvHit.Pt() ); + mTreeData.thaX.push_back( htv3.X() ); mTreeData.thaY.push_back( htv3.Y() ); mTreeData.thaZ.push_back( htv3.Z() ); mTreeData.thdN++; - } - } + } // hit + } // i } //FillTrackDeltas void StFwdTrackMaker::FitVertex(){ - const auto &genfitTracks = mForwardTracker -> globalTracks(); + vector genfitTracks; + const auto &trackResults = mForwardTracker -> getTrackResults(); + for ( auto gtr : trackResults ){ + genfitTracks.push_back( gtr.track ); + } if ( genfitTracks.size() >= 2 ){ genfit::GFRaveVertexFactory gfrvf; @@ -1397,13 +1423,11 @@ void StFwdTrackMaker::FitVertex(){ bscm(1, 1) = bssXY*bssXY; bscm(2, 2) = 50.5 * 50.5; gfrvf.setBeamspot( TVector3( 0, 0, 0 ), bscm ); - // std::vector< genfit::GFRaveVertex * > vertices; - const auto &genfitTracks = mForwardTracker -> globalTracks(); + mRaveVertices.clear(); gfrvf.findVertices( &mRaveVertices, genfitTracks, false ); LOG_DEBUG << "mRaveVertices.size() = " << mRaveVertices.size() << endm; - for ( auto vert : mRaveVertices ){ LOG_DEBUG << TString::Format( "RAVE vertex @(%f, %f, %f)\n\n", vert->getPos().X(), vert->getPos().Y(), vert->getPos().Z() ) << endm; } @@ -1564,7 +1588,17 @@ void StFwdTrackMaker::FillTTree(){ int idt = 0; double qual = 0; - idt = MCTruthUtils::dominantContribution(fittedTracks[i].trackSeed, qual); + // TODO: combine these + Seed_t combinedSeed; + combinedSeed.insert( combinedSeed.begin(), fittedTracks[i].fstSeed.begin(), fittedTracks[i].fstSeed.end() ); // this is goofed but will fix + combinedSeed.insert( combinedSeed.end(), fittedTracks[i].fttSeed.begin(), fittedTracks[i].fttSeed.end() ); + idt = MCTruthUtils::dominantContribution(combinedSeed, qual); + if ( idt == -1 ){ + idt = MCTruthUtils::dominantContribution(fittedTracks[i].fttSeed, qual); + if ( idt == -1 ){ + idt = MCTruthUtils::dominantContribution(fittedTracks[i].fstSeed, qual); + } + } if ( fittedTracks[i].track == nullptr || fittedTracks[i].trackRep == nullptr ) { LOG_INFO << "Skip saving null track" << endm; @@ -1586,15 +1620,15 @@ void StFwdTrackMaker::FillTTree(){ mTreeData.rcPhi.push_back( fittedTracks[i].momentum.Phi() ); mTreeData.rcNumPV.push_back( fittedTracks[i].nPV ); - mTreeData.rcNumFTT.push_back( fittedTracks[i].nFTT ); - mTreeData.rcNumFST.push_back( fittedTracks[i].nFST ); + mTreeData.rcNumFTT.push_back( fittedTracks[i].numFTT() ); + mTreeData.rcNumFST.push_back( fittedTracks[i].numFST() ); mTreeData.rcN ++; } LOG_INFO << "Filling TTree" << endm; mTree->Fill(); } // if mGenTree -} +} // FillTree //________________________________________________________________________ @@ -1638,6 +1672,9 @@ void StFwdTrackMaker::Clear(const Option_t *opts) { mTreeData.mcPhi.clear(); mTreeData.mcVertexId.clear(); mTreeData.mcCharge.clear(); + mTreeData.mcNumFst.clear(); + mTreeData.mcNumFtt.clear(); + mTreeData.vmcX.clear(); mTreeData.vmcY.clear(); mTreeData.vmcZ.clear(); @@ -1653,6 +1690,8 @@ void StFwdTrackMaker::Clear(const Option_t *opts) { mTreeData.vrcZ.clear(); mTreeData.thdX.clear(); mTreeData.thdY.clear(); + mTreeData.thdP.clear(); + mTreeData.thdR.clear(); mTreeData.thaX.clear(); mTreeData.thaY.clear(); mTreeData.thaZ.clear(); diff --git a/StRoot/StFwdTrackMaker/StFwdTrackMaker.h b/StRoot/StFwdTrackMaker/StFwdTrackMaker.h index 5be3f4f0f11..2d0d5a22723 100644 --- a/StRoot/StFwdTrackMaker/StFwdTrackMaker.h +++ b/StRoot/StFwdTrackMaker/StFwdTrackMaker.h @@ -49,50 +49,54 @@ class GenfitTrackResult; const size_t MAX_TREE_ELEMENTS = 4000; struct FwdTreeData { - // hits; + /**@brief Ftt hit related info*/ int fttN; vector fttX, fttY, fttZ; vector fttVolumeId; - // Only avalaible for hits if MC + // Note: Below are only avalaible for hits if MC vector fttPt; vector fttTrackId, fttVertexId; - // hits; + /**@brief Fst hit related info*/ int fstN; vector fstX, fstY, fstZ; vector fstTrackId; + /**@brief Fcs hit related info*/ int fcsN; - vector fcsX, fcsY, fcsZ; + vector fcsX, fcsY, fcsZ, fcsE; vector fcsDet; - // RC tracks + /**@brief RC track related info*/ int rcN; vector rcPt, rcEta, rcPhi, rcQuality; vector rcTrackId, rcNumFST, rcCharge, rcNumFTT, rcNumPV; - // MC Tracks + /**@brief MC Track related info*/ int mcN; vector mcPt, mcEta, mcPhi; vector mcVertexId, mcCharge; + vector mcNumFtt, mcNumFst; - // MC Level vertex info - // maybe use also for TPC vertex if available in data + /**@brief MC Vertex related info*/ int vmcN; vector vmcX, vmcY, vmcZ; + /**@brief Track Projection related info*/ int tprojN; vector tprojX, tprojY, tprojZ; vector tprojPx, tprojPy, tprojPz; vector tprojIdD, tprojIdT; - // RAVE reco vertices + /**@brief RAVE RC Vertex related info*/ int vrcN; vector vrcX, vrcY, vrcZ; + /**@brief Track-to-hit delta related info*/ int thdN; - vector thdX, thdY, thaX, thaY, thaZ; + vector thdX, thdY, thdP, thdR, thaX, thaY, thaZ; + /**@brief Seed finding Criteria related info*/ bool saveCrit = false; std::map> Crits; std::map> CritTrackIds; @@ -132,8 +136,8 @@ class StFwdTrackMaker : public StMaker { // for Wavefront OBJ export - size_t eventIndex = 0; - + size_t eventIndex = 0; // counts up for processed events + size_t mEventNum = 0; // global event num (index) bool mGenHistograms = false; bool mGenTree = false; diff --git a/StRoot/StFwdTrackMaker/include/Tracker/FwdHit.h b/StRoot/StFwdTrackMaker/include/Tracker/FwdHit.h index bcf4cb1fa7b..3ebc7f7fd88 100644 --- a/StRoot/StFwdTrackMaker/include/Tracker/FwdHit.h +++ b/StRoot/StFwdTrackMaker/include/Tracker/FwdHit.h @@ -50,14 +50,14 @@ class McTrack { mStartVertex = start_vertex; } - void addHit(KiTrack::IHit *hit) { mHits.push_back(hit); } - // void addFstHit(KiTrack::IHit *hit) { mFstHits.push_back(hit); } + void addHit(KiTrack::IHit *hit) { mFttHits.push_back(hit); } + void addFstHit(KiTrack::IHit *hit) { mFstHits.push_back(hit); } double mPt, mEta, mPhi; int mTid, mQ, mStartVertex; - std::vector mHits; - // std::vector mFstHits; + std::vector mFttHits; + std::vector mFstHits; }; diff --git a/StRoot/StFwdTrackMaker/include/Tracker/FwdTracker.h b/StRoot/StFwdTrackMaker/include/Tracker/FwdTracker.h index 1316f9eca61..c110e198212 100644 --- a/StRoot/StFwdTrackMaker/include/Tracker/FwdTracker.h +++ b/StRoot/StFwdTrackMaker/include/Tracker/FwdTracker.h @@ -10,6 +10,7 @@ #include "TRandom3.h" #include "TTree.h" #include "TVector3.h" +#include "TLorentzVector.h" #include #include @@ -49,6 +50,10 @@ struct MCTruthUtils { truth[ fhit->_tid ]++; } + if ( truth.size() == 0 ){ + return -1; + } + using namespace std; using P = decltype(truth)::value_type; auto dom = max_element(begin(truth), end(truth), [](P a, P b){ return a.second < b.second; }); @@ -66,14 +71,20 @@ struct MCTruthUtils { class GenfitTrackResult { public: - GenfitTrackResult( size_t nFTT, size_t nFST, - Seed_t &seedTrack, genfit::Track *track ) { - this->nFST = nFST; - this->nFTT = nFTT; - this->trackSeed = seedTrack; + GenfitTrackResult(){} + GenfitTrackResult( Seed_t &fttSeed, + Seed_t &fstSeed, + genfit::Track *track ) { + set( fttSeed, fstSeed, track ); + } +void set( Seed_t &fttSeed, + Seed_t &fstSeed, + genfit::Track *track ){ + this->fttSeed = fttSeed; + this->fstSeed = fstSeed; try { - this->track = track; + this->track = new genfit::Track(*track); this->status = *(this->track->getFitStatus()); this->trackRep = this->track->getCardinalRep(); @@ -83,70 +94,45 @@ class GenfitTrackResult { this->nFailedPoints = this->status.getNFailedPoints(); this->charge = this->status.getCharge(); - this->nPV = this->track->getNumPoints() - (nFTT + nFST); + this->nPV = this->track->getNumPoints() - (numFTT() + numFST()); this->momentum = this->trackRep->getMom( this->track->getFittedState(0, this->trackRep) ); LOG_DEBUG << "GenfitTrackResult::set Track successful" << endm; + size_t numTotal = numFTT() + numFST(); + LOG_DEBUG << "GTR has track->numPoints=" << this->track->getNumPoints() << " vs. numFST: " << numFST() << ", numFTT: " << numFTT() << endm; } catch ( genfit::Exception &e ) { - LOG_ERROR << "GenfitTrackResult cannot get track" << endm; + LOG_ERROR << "CANNOT GET TRACK" << endm; this->track = nullptr; this->trackRep = nullptr; this->isFitConverged = false; this->isFitConvergedFully = false; this->isFitConvergedPartially = false; - this->nFailedPoints = nFST + nFTT; + this->nFailedPoints = 99; this->charge = 0; } } - ~GenfitTrackResult() { - // MEMORY LEAK - // LOG_INFO << "~GenfitTrackResult" << endm; - // if (this->track) - // delete this->track; - // this->track = nullptr; + void addFST( Seed_t &nfstSeed, genfit::Track *track ) { + set( fttSeed, nfstSeed, track ); + } + void addFTT( Seed_t &nfttSeed, genfit::Track *track ) { + set( nfttSeed, fstSeed, track ); } - void setFst( Seed_t &seedFst, genfit::Track *track ){ - LOG_DEBUG << "GenfitTrackResult::setFSt" << endm; - nFST = seedFst.size(); - fstSeed = seedFst; - - try { - this->fstTrack = new genfit::Track(*track); - // make sure the McTrackId is set correctly - this->fstTrack->setMcTrackId( this->track->getMcTrackId() ); - this->fstStatus = *(this->fstTrack->getFitStatus()); - this->fstTrackRep = this->fstTrack->getCardinalRep(); - - this->isFitConverged = this->fstStatus.isFitConverged(); - this->isFitConvergedFully = this->fstStatus.isFitConvergedFully(); - this->isFitConvergedPartially = this->fstStatus.isFitConvergedPartially(); - this->nFailedPoints = this->fstStatus.getNFailedPoints(); - this->charge = this->fstStatus.getCharge(); - this->fstMomentum = this->fstTrackRep->getMom( this->fstTrack->getFittedState(0, this->fstTrackRep) ); + size_t numFTT() const { return fttSeed.size();} + size_t numFST() const { return fstSeed.size();} - } catch ( genfit::Exception &e ) { - LOG_ERROR << "CANNOT GET FST TRACK" << endm; - this->fstTrack = nullptr; - this->fstTrackRep = nullptr; - - this->fstIsFitConverged = false; - this->fstIsFitConvergedFully = false; - this->fstIsFitConvergedPartially = false; - this->fstNFailedPoints = nFST + nFTT; - this->fstCharge = 0; - } + ~GenfitTrackResult() { } - Seed_t trackSeed; + + Seed_t fttSeed; Seed_t fstSeed; + TVector3 pv; TVector3 momentum; double charge; - size_t nFST = 0; - size_t nFTT = 0; size_t nPV = 0; genfit::FitStatus status; genfit::AbsTrackRep *trackRep = nullptr; @@ -155,22 +141,13 @@ class GenfitTrackResult { bool isFitConvergedFully = false; bool isFitConvergedPartially = false; size_t nFailedPoints = 0; + size_t theNumFTT = 0; - // Result after FST refit - genfit::Track *fstTrack = nullptr; - genfit::AbsTrackRep *fstTrackRep = nullptr; - genfit::FitStatus fstStatus; - bool fstIsFitConverged = false; - bool fstIsFitConvergedFully = false; - bool fstIsFitConvergedPartially = false; - size_t fstNFailedPoints = 0; - double fstCharge = 0; - TVector3 fstMomentum; - + void summary() { - LOG_INFO << TString::Format( "TrackResult[p=(%f, %f, %f)/(%f, %f, %f), q=%f, nFTT=%lu, nFST=%lu, nPV=%lu, isFitConvergedFully=%d]", momentum.X(), momentum.Y(), momentum.Z(), momentum.Pt(), momentum.Eta(), momentum.Phi(), charge, nFTT, nFST, nPV, isFitConvergedFully ).Data() << endm; + LOG_INFO << TString::Format( "TrackResult[p=(%f, %f, %f)/(%f, %f, %f), q=%f, nFTT=%lu, nFST=%lu, nPV=%lu, isFitConvergedFully=%d]", momentum.X(), momentum.Y(), momentum.Z(), momentum.Pt(), momentum.Eta(), momentum.Phi(), charge, fttSeed.size(), fstSeed.size(), nPV, isFitConvergedFully ).Data() << endm; } -}; +}; // GenfitTrackResult class ForwardTrackMaker { public: @@ -179,17 +156,27 @@ class ForwardTrackMaker { } const std::vector &getTrackResults() const { return mTrackResults; } - const std::vector &getRecoTracks() const { return mRecoTracks; } - const std::vector &getFitMomenta() const { return mFitMoms; } - const std::vector &getNumFstHits() const { return mNumFstHits; } - const std::vector &getFitStatus() const { return mFitStatus; } - const std::vector &globalTrackReps() const { return mGlobalTrackReps; } - const std::vector &globalTracks() const { return mGlobalTracks; } - + const std::vector &getTrackSeeds() const { return mRecoTracks; } + // const std::vector &getFitMomenta() const { return mFitMoms; } + // const std::vector &getNumFstHits() const { return mNumFstHits; } + // const std::vector &getFitStatus() const { return mFitStatus; } + // const std::vector &globalTrackReps() const { return mGlobalTrackReps; } + // const std::vector &globalTracks() const { return mGlobalTracks; } + + /** + * @brief Set the Config File object + * + * @param cf : config filename + */ void setConfigFile(std::string cf) { mConfigFile = cf; } + /** + * @brief Set the Save Criteria Values object + * + * @param save : true to save crit values + */ void setSaveCriteriaValues(bool save) { mSaveCriteriaValues = save; } @@ -199,6 +186,12 @@ class ForwardTrackMaker { // Adopt external hit loader void setData(std::shared_ptrdata) { mDataSource = data; } + /** + * @brief Initialize FwdTracker + * + * @param geoCache : name of cached geometry file + * @param genHistograms : generate histograms + */ virtual void initialize( TString geoCache, bool genHistograms) { mGenHistograms = genHistograms; if (mGenHistograms) setupHistograms(); @@ -208,9 +201,12 @@ class ForwardTrackMaker { if (!mConfig.exists("TrackFitter")) mDoTrackFitting = false; - } - + } //initialize + /** + * @brief Writes QA histograms and config to ROOT file + * + */ void writeEventHistograms() { // no file, dont write anything @@ -229,14 +225,14 @@ class ForwardTrackMaker { mTrackFitter->writeHistograms(); gDirectory->cd(""); mQualityPlotter->writeHistograms(); - } + } //writeEventHistograms - /** Loads Criteria from XML configuration. - * - * Utility function for loading criteria from XML config. - * - * @return vector of ICriterion pointers - */ + /** + * @brief Loads Criteria from XML configuration. + * Utility function for loading criteria from XML config. + * @param path : path in config to load + * @return vector of ICriterion pointers + */ std::vector loadCriteria(string path) { std::vector crits; @@ -270,8 +266,14 @@ class ForwardTrackMaker { } return crits; - } + } // loadCriteria + /** + * @brief Get the Criteria Values object + * + * @param crit_name : Criteria to get + * @return std::vector : list of values + */ std::vector getCriteriaValues(std::string crit_name) { std::vector em; if (mSaveCriteriaValues != true) { @@ -293,8 +295,14 @@ class ForwardTrackMaker { } return em; - }; + } //getCriteriaValues + /** + * @brief Get the Criteria All Values object + * + * @param crit_name : Criteria values to get + * @return std::vector> : map of values + */ std::vector> getCriteriaAllValues(std::string crit_name) { std::vector> em; if (mSaveCriteriaValues != true) { @@ -316,8 +324,14 @@ class ForwardTrackMaker { } return em; - }; + } // getCriteriaAllValues + /** + * @brief Get the Criteria Track Ids object + * + * @param crit_name : Name of criteria to get track ids for + * @return std::vector : list of track ids + */ std::vector getCriteriaTrackIds(std::string crit_name) { std::vector em; if (mSaveCriteriaValues != true) { @@ -339,8 +353,12 @@ class ForwardTrackMaker { } return em; - }; + } //getCriteriaTrackIds + /** + * @brief Clear the saved values for two hit and three hit criteria + * + */ void clearSavedCriteriaValues() { if (mSaveCriteriaValues != true) { return; @@ -355,8 +373,14 @@ class ForwardTrackMaker { auto critKeeper = static_cast(crit); critKeeper->clear(); } - } + } // clearSavedCriteria + /** + * @brief Determine the total num of hits in the hitmap + * + * @param hitmap : hitmap to consider + * @return size_t : total num of hits + */ size_t nHitsInHitMap(FwdDataSource::HitMap_t &hitmap) { size_t n = 0; @@ -367,6 +391,12 @@ class ForwardTrackMaker { return n; } + /** + * @brief Counts the number of tracks with nHits hits + * + * @param nHits : number of hits + * @return size_t : number of tracks of the given length + */ size_t countRecoTracks(size_t nHits) { size_t n = 0; @@ -378,6 +408,10 @@ class ForwardTrackMaker { return n; } + /** + * @brief Create QA histograms + * + */ void setupHistograms() { mHist["input_nhits"] = new TH1I("input_nhits", ";# hits", 1000, 0, 1000); @@ -392,23 +426,30 @@ class ForwardTrackMaker { FwdTrackerUtils::labelAxis(mHist["FitStatus"]->GetXaxis(), {"Seeds", "AttemptFit", "GoodFit", "BadFit", "GoodCardinal", "PossibleReFit", "AttemptReFit", "GoodReFit", "BadReFit", "w3Si","w2Si", "w1Si", "w0Si" }); mHist["FitDuration"] = new TH1I("FitDuration", ";Duration (ms)", 5000, 0, 50000); - mHist["nSiHitsFound"] = new TH2I( "nSiHitsFound", ";Si Disk; n Hits", 5, 0, 5, 10, 0, 10 ); + mHist["nFstHitsFound"] = new TH2I( "nFstHitsFound", ";FST Disk; n Hits", 5, 0, 5, 10, 0, 10 ); mHist["Step1Duration"] = new TH1I("Step1Duration", ";Duration (ms)", 500, 0, 500); mHist["Step2Duration"] = new TH1I("Step2Duration", ";Duration (ms)", 500, 0, 500); mHist["Step3Duration"] = new TH1I("Step3Duration", ";Duration (ms)", 500, 0, 500); mHist["Step4Duration"] = new TH1I("Step4Duration", ";Duration (ms)", 500, 0, 500); - } - void fillHistograms() { + mHist["FttSearchDeltaX"] = new TH1I("FttSearchDeltaX", ";dx", 500, -100, 100); + mHist["FttSearchDeltaY"] = new TH1I("FttSearchDeltaY", ";dy", 500, -100, 100); + mHist["FttSearchDeltaR"] = new TH1I("FttSearchDeltaR", ";dr", 500, -100, 100); + mHist["FttSearchDeltaP"] = new TH1I("FttSearchDeltaP", ";d#phi", 500, -3.14, 3.4); - if (mGenHistograms && mDataSource != nullptr) { - auto hm = mDataSource->getFttHits(); - for (auto hp : hm) - mHist["input_nhits"]->Fill(hp.second.size()); - } - } + mHist["FttSearchMinDeltaX"] = new TH1I("FttSearchMinDeltaX", ";dx", 500, -100, 100); + mHist["FttSearchMinDeltaY"] = new TH1I("FttSearchMinDeltaY", ";dy", 500, -100, 100); + mHist["FttSearchMinDeltaR"] = new TH1I("FttSearchMinDeltaR", ";dr", 500, -100, 100); + mHist["FttSearchMinDeltaP"] = new TH1I("FttSearchMinDeltaP", ";d#phi", 500, -3.14, 3.4); + + mHist["nFttHitsFound"] = new TH2I( "nFttHitsFound", ";Ftt Plane; n Hits", 5, 0, 5, 10, 0, 10 ); + } //setupHistograms + /** + * @brief Writes histograms to file + * + */ void writeHistograms() { if ( !mGenHistograms ){ return; @@ -420,38 +461,12 @@ class ForwardTrackMaker { } } - // this is the main event loop. doEvent processes a single event iEvent... - void make() { - - int single_event = mConfig.get("Input:event", -1); - - if (single_event >= 0) { - doEvent(single_event); - return; - } - - unsigned long long firstEvent = mConfig.get("Input:first-event", 0); - - if (mConfig.exists("Input:max-events")) { - unsigned long long maxEvents = mConfig.get("Input:max-events", 0); - - if (nEvents > maxEvents) - nEvents = maxEvents; - - } - - // loop over events - - for (unsigned long long iEvent = firstEvent; iEvent < firstEvent + nEvents; iEvent++) { - doEvent(iEvent); - } - - if (mGenHistograms){ - mQualityPlotter->finish(); - writeEventHistograms(); - } - } - + /** + * @brief Remove used hits from the hit map + * + * @param hitmap : hitmap with hits used this round + * @param tracks : tracks formed from hits + */ void removeHits(FwdDataSource::HitMap_t &hitmap, std::vector &tracks) { for (auto track : tracks) { @@ -469,7 +484,12 @@ class ForwardTrackMaker { } // loop on track } // removeHits - void doEvent(unsigned long long int iEvent = 0) { + /** + * @brief Perform track finding and fitting on an event + * + * @param iEvent : event number + */ + void doEvent() { /************** Cleanup ****************************************/ // Moved cleanup to the start of doEvent, so that the fit results // persist after the call @@ -507,33 +527,56 @@ class ForwardTrackMaker { // Load and sort the hits /*************************************************************/ long long itStart = FwdTrackerUtils::nowNanoSecond(); - FwdDataSource::HitMap_t &hitmap = mDataSource->getFttHits();; + FwdDataSource::HitMap_t &fttHitmap = mDataSource->getFttHits(); + FwdDataSource::HitMap_t &fstHitmap = mDataSource->getFstHits(); + + string hitmapSource = mConfig.get("TrackFinder:source", "ftt"); + bool useFttAsSource = !(hitmapSource == "fst"); + + if ( useFttAsSource == false ){ + for (auto hp : fstHitmap){ + LOG_DEBUG << "HITMAP [" << hp.first << ", { "; + for ( auto h: hp.second ){ + LOG_DEBUG << "z=" << h->getZ() << " "; + } + LOG_DEBUG << " }" << endm; + } + } FwdDataSource::McTrackMap_t &mcTrackMap = mDataSource->getMcTracks(); - fillHistograms(); + if (mGenHistograms && mDataSource != nullptr) { + auto hm = mDataSource->getFttHits(); + for (auto hp : hm) + mHist["input_nhits"]->Fill(hp.second.size()); + } + long long duration = (FwdTrackerUtils::nowNanoSecond() - itStart) * 1e-6; // milliseconds if (mGenHistograms) mHist["Step1Duration"]->Fill( duration ); - bool mcTrackFinding = true; - - if (mConfig.exists("TrackFinder")) + if (mConfig.exists("TrackFinder") && mConfig.get( "TrackFinder:mc", false ) == false ){ mcTrackFinding = false; + } + if (mConfig.exists("TrackFinder") && mConfig.get( "TrackFinder:active", false ) == false){ + mcTrackFinding = true; + } /***********************************************/ // MC Track Finding if (mcTrackFinding) { LOG_DEBUG << "MC TRACK FINDING " << endm; - doMcTrackFinding(mcTrackMap); + doMcTrackFinding(mcTrackMap, useFttAsSource); /***********************************************/ - // REFIT with Silicon hits - if (mConfig.get("TrackFitter:refitSi", true)) { - addSiHitsMc(); + // REFIT with FST or FTT hits accordingly + if (mConfig.get("TrackFitter:refit", false)) { + if (useFttAsSource) + addFstHitsMc(); + else + addFttHitsMc(); } else { - LOG_DEBUG << "Skipping FST Hits" << endm; - // skip Si refit + LOG_DEBUG << "Skipping Refit (MC tracking)" << endm; } /***********************************************/ @@ -555,16 +598,24 @@ class ForwardTrackMaker { // plus initial fit size_t nIterations = mConfig.get("TrackFinder:nIterations", 0); for (size_t iIteration = 0; iIteration < nIterations; iIteration++) { - doTrackIteration(iIteration, hitmap); - } + if (useFttAsSource){ + doTrackIteration(iIteration, fttHitmap); + } else { + doTrackIteration(iIteration, fstHitmap); + } + } // iIteration /***********************************************/ /***********************************************/ - // REFIT with Silicon hits - if (mConfig.get("TrackFitter:refitSi", true)) { - addSiHits(); + // REFIT with Fst or Ftt hits + if (mConfig.get("TrackFitter:refit", false)) { + if (useFttAsSource){ + addFstHits(); + } else { + addFttHits(); + } } else { - // Skipping Si Refit + LOG_DEBUG << "Skipping track refit" << endm; } /***********************************************/ @@ -573,65 +624,47 @@ class ForwardTrackMaker { } } // doEvent - void fitTrack(Seed_t &track) { + /** + * @brief Perform a single fit from seed points + * + * @param seed : seed points from either FTT or FST + */ + void fitTrack(Seed_t &seed) { + + // make sure this is consistent with above + string hitmapSource = mConfig.get("TrackFinder:source", "ftt"); + bool useFttAsSource = !(hitmapSource == "fst"); if ( mGenHistograms ){ mHist["FitStatus"]->Fill("Seeds", 1); } + if (!mDoTrackFitting) + return; + + /*******************************************************/ + // Only for Simulation // Calculate the MC info first and check filters int idt = 0; double qual = 0; - idt = MCTruthUtils::dominantContribution(track, qual); + // Get the quality and MC truth id + idt = MCTruthUtils::dominantContribution(seed, qual); - - - TVector3 mcSeedMom; - auto mctm = mDataSource->getMcTracks(); - // get the MC track momentum if we can - if (mctm.count(idt)) { - auto mct = mctm[idt]; - mcSeedMom.SetPtEtaPhi(mct->mPt, mct->mEta, mct->mPhi); - } - - - // Mc Filter - bool bailout = false; - if (qual < mConfig.get("TrackFitter.McFilter:quality-min", 0.0)) { - bailout = true; - // LOG_INFO << "BAIL OUT on Fit bc quality = " << qual << endm; - } + // get the MC track momentum if we can (may be used for state seed) if (mctm.count(idt)) { auto mct = mctm[idt]; mcSeedMom.SetPtEtaPhi(mct->mPt, mct->mEta, mct->mPhi); - if (mct->mPt < mConfig.get("TrackFitter.McFilter:pt-min", 0.0) || - mct->mPt > mConfig.get("TrackFitter.McFilter:pt-max", 1e10)) { - bailout = true; - // LOG_INFO << "BAIL OUT on Fit bc Pt = " << mct->mPt << endm; - } - if (mct->mEta < mConfig.get("TrackFitter.McFilter:eta-min", 0) || - mct->mEta > mConfig.get("TrackFitter.McFilter:eta-max", 1e10)) { - bailout = true; - // LOG_INFO << "BAIL OUT on Fit bc eta = " << mct->mEta << endm; - } - - } else { - // cannot find the track } + /*******************************************************/ - bailout = false; - - TVector3 p; - p.SetPtEtaPhi( 0, -999, 0 ); genfit::FitStatus fitStatus; - genfit::AbsTrackRep *trackRep = nullptr;//new genfit::RKTrackRep(211); // pdg for pi+ genfit::Track *genTrack = nullptr;//new genfit::Track( trackRep, TVector3(0, 0, 0), TVector3(0, 0, 0) ); - if (mDoTrackFitting && !bailout) { + if (mDoTrackFitting) { if ( mGenHistograms ){ mHist["FitStatus"]->Fill("AttemptFit", 1); } @@ -639,59 +672,53 @@ class ForwardTrackMaker { double vertex[3] = { mEventVertex.X(), mEventVertex.Y(), mEventVertex.Z() }; double * pVertex = 0; - if ( fabs(mEventVertex.X()) < 100 ){ + if ( fabs(mEventVertex.X()) < 150 ){ pVertex = vertex; // only use it if it has been set from default } - if (true == mConfig.get("TrackFitter:mcSeed", false)) { // use the MC pt, eta, phi as the seed for fitting - p = mTrackFitter->fitTrack(track, pVertex, &mcSeedMom); + mTrackFitter->fitTrack(seed, pVertex, &mcSeedMom); } else { // Normal case, real data - p = mTrackFitter->fitTrack(track, pVertex); + mTrackFitter->fitTrack(seed, pVertex); } if ( mGenHistograms ){ - if (p.Perp() > 1e-3) { + if (mTrackFitter->getStatus().isFitConvergedFully()) { mHist["FitStatus"]->Fill("GoodFit", 1); } else { mHist["FitStatus"]->Fill("BadFit", 1); } } - genTrack = new genfit::Track(*mTrackFitter->getTrack()); genTrack->setMcTrackId(idt); - GenfitTrackResult gtr( track.size(), 0, track, genTrack ); - - // assign the fit results to be saved - fitStatus = mTrackFitter->getStatus(); - trackRep = mTrackFitter->getTrackRep()->clone(); // Clone the track rep - - if ( mGenHistograms && genTrack->getFitStatus(genTrack->getCardinalRep())->isFitConverged() && p.Perp() > 1e-3) { + Seed_t nonSeeds; // none + GenfitTrackResult gtr; + if ( useFttAsSource ) + gtr.set( seed, nonSeeds, genTrack ); + else + gtr.set( nonSeeds, seed, genTrack ); + + if ( mGenHistograms && genTrack->getFitStatus(genTrack->getCardinalRep())->isFitConverged()) { mHist["FitStatus"]->Fill("GoodCardinal", 1); } - - // Save everything (now use GenfitTrackResult) - mFitMoms.push_back(p); - mGlobalTracks.push_back(genTrack); - mGlobalTrackReps.push_back(trackRep); - mFitStatus.push_back(fitStatus); - mRecoTrackQuality.push_back(qual); - mRecoTrackIdTruth.push_back(idt); - mNumFstHits.push_back(0); - mTrackResults.push_back( gtr ); LOG_DEBUG << "FwdTracker::fitTrack complete" << endm; - } // if (mDoTrackFitting && !bailout) - } + } // if (mDoTrackFitting) + } // fitTrack - void doTrackFitting( std::vector &tracks) { + /** + * @brief Loop on track seeds and fit each one + * + * @param trackSeeds : Track seeds + */ + void doTrackFitting( std::vector &trackSeeds) { long long itStart = FwdTrackerUtils::nowNanoSecond(); // Fit each accepted track seed - for (auto t : tracks) { + for (auto t : trackSeeds) { fitTrack(t); } long long itEnd = FwdTrackerUtils::nowNanoSecond(); @@ -703,7 +730,13 @@ class ForwardTrackMaker { } - void doMcTrackFinding(FwdDataSource::McTrackMap_t &mcTrackMap) { + /** + * @brief MC track finding builds track seeds from available hits using MC association + * + * @param mcTrackMap : Mc tracks + * @param useFttAsSource : Use FTT for seeds or (false) use Fst + */ + void doMcTrackFinding(FwdDataSource::McTrackMap_t &mcTrackMap, bool useFttAsSource = true) { mQualityPlotter->startIteration(); @@ -711,18 +744,31 @@ class ForwardTrackMaker { for (auto kv : mcTrackMap) { auto mc_track = kv.second; - if (mc_track->mHits.size() < 4){ // require min 4 hits on track + + + if (useFttAsSource && mc_track->mFttHits.size() < 4){ // require min 4 FTT hits on track + continue; + } + + if (!useFttAsSource && mc_track->mFstHits.size() < 3 ) { // require min 3 FST hits on track continue; } std::set uvid; Seed_t track; - for (auto h : mc_track->mHits) { - track.push_back(h); - uvid.insert(static_cast(h)->_vid); + if ( useFttAsSource ){ // FTT (DEFAULT) + for (auto h : mc_track->mFttHits) { + track.push_back(h); + uvid.insert(static_cast(h)->_vid); + } + } else { // FST + for (auto h : mc_track->mFstHits) { + track.push_back(h); + uvid.insert(static_cast(h)->_vid); + } } - + if (uvid.size() == track.size()) { // only add tracks that have one hit per volume mRecoTracks.push_back(track); int idt = 0; @@ -742,7 +788,7 @@ class ForwardTrackMaker { if ( mGenHistograms ){ mQualityPlotter->afterIteration(0, mRecoTracks); } - } + } //doMcTrackFinding /** sliceHitMapInPhi @@ -857,6 +903,7 @@ class ForwardTrackMaker { if (doAutomation) { automaton.doAutomaton(); } else { + LOG_DEBUG << "Skipping Automation step" << endm; //Not running Automation Step } @@ -937,6 +984,12 @@ class ForwardTrackMaker { return acceptedTracks; } // doTrackingOnHitmapSubset + /** + * @brief Main tracking procedure + * + * @param iIteration : The track iteration + * @param hitmap : the hitmap of available hits per plane + */ void doTrackIteration(size_t iIteration, FwdDataSource::HitMap_t &hitmap) { // empty the list of reco tracks for the iteration @@ -945,7 +998,8 @@ class ForwardTrackMaker { // check to see if we have hits! size_t nHitsThisIteration = nHitsInHitMap(hitmap); - if (nHitsThisIteration < 4) { + const int minHitsToConsider = 3; + if (nHitsThisIteration < minHitsToConsider) { // No hits left in the hitmap! Skipping this iteration return; } @@ -988,7 +1042,7 @@ class ForwardTrackMaker { size_t nHitsThisSlice = 0; if ( phi_slice_count > 1 ){ nHitsThisSlice = sliceHitMapInPhi( hitmap, slicedHitMap, phi_min, phi_max ); - if ( nHitsThisSlice < 4 ) { + if ( nHitsThisSlice < minHitsToConsider ) { continue; } } else { // no need to slice @@ -1017,7 +1071,7 @@ class ForwardTrackMaker { LOG_DEBUG << " FITTING " << mRecoTracksThisItertion.size() << " now" << endm; - if ( mRecoTracksThisItertion.size() < 201 ){ + if ( mRecoTracksThisItertion.size() < 2001 ){ doTrackFitting( mRecoTracksThisItertion ); } else { LOG_ERROR << "BAILING OUT of fit, too many track candidates" << endm; @@ -1032,14 +1086,18 @@ class ForwardTrackMaker { } // doTrackIteration - void addSiHitsMc() { + /** + * @brief Adds compatible FST hits to tracks seeded with FTT + * + */ + void addFstHitsMc() { FwdDataSource::HitMap_t hitmap = mDataSource->getFstHits(); for (size_t i = 0; i < mTrackResults.size(); i++) { GenfitTrackResult >r = mTrackResults[i]; if ( gtr.status.isFitConverged() == false || gtr.momentum.Perp() < 1e-3) { - LOG_DEBUG << "Skipping addSiHitsMc, fit failed" << endm; + LOG_DEBUG << "Skipping addFstHitsMc, fit failed" << endm; return; } @@ -1047,57 +1105,215 @@ class ForwardTrackMaker { mHist["FitStatus"]->Fill("PossibleReFit", 1); } - Seed_t si_hits_for_this_track(3, nullptr); + Seed_t fst_hits_for_this_track(3, nullptr); for (size_t j = 0; j < 3; j++) { for (auto h0 : hitmap[j]) { if (dynamic_cast(h0)->_tid == gtr.track->getMcTrackId()) { - si_hits_for_this_track[j] = h0; + fst_hits_for_this_track[j] = h0; break; } } // loop on hits in this layer of hitmap } // loop on hitmap layers - size_t nSiHitsFound = 0; - if ( si_hits_for_this_track[0] != nullptr ) nSiHitsFound++; - if ( si_hits_for_this_track[1] != nullptr ) nSiHitsFound++; - if ( si_hits_for_this_track[2] != nullptr ) nSiHitsFound++; - LOG_DEBUG << "Found " << nSiHitsFound << " FST Hits on this track (MC lookup)" << endm; + size_t nFstHitsFound = 0; + if ( fst_hits_for_this_track[0] != nullptr ) nFstHitsFound++; + if ( fst_hits_for_this_track[1] != nullptr ) nFstHitsFound++; + if ( fst_hits_for_this_track[2] != nullptr ) nFstHitsFound++; + LOG_DEBUG << "Found " << nFstHitsFound << " FST Hits on this track (MC lookup)" << endm; if ( mGenHistograms ){ - this->mHist[ "nSiHitsFound" ]->Fill( 1, ( si_hits_for_this_track[0] != nullptr ? 1 : 0 ) ); - this->mHist[ "nSiHitsFound" ]->Fill( 2, ( si_hits_for_this_track[1] != nullptr ? 1 : 0 ) ); - this->mHist[ "nSiHitsFound" ]->Fill( 3, ( si_hits_for_this_track[2] != nullptr ? 1 : 0 ) ); + this->mHist[ "nFstHitsFound" ]->Fill( 1, ( fst_hits_for_this_track[0] != nullptr ? 1 : 0 ) ); + this->mHist[ "nFstHitsFound" ]->Fill( 2, ( fst_hits_for_this_track[1] != nullptr ? 1 : 0 ) ); + this->mHist[ "nFstHitsFound" ]->Fill( 3, ( fst_hits_for_this_track[2] != nullptr ? 1 : 0 ) ); } - if (nSiHitsFound >= 1) { + if (nFstHitsFound >= 1) { if ( mGenHistograms ){ mHist["FitStatus"]->Fill("AttemptReFit", 1); } - TVector3 p = mTrackFitter->refitTrackWithSiHits(gtr.track, si_hits_for_this_track); + TVector3 p = mTrackFitter->refitTrackWithFstHits(gtr.track, fst_hits_for_this_track); if ( mGenHistograms ){ if (p.Perp() == mFitMoms[i].Perp()) { mHist["FitStatus"]->Fill("BadReFit", 1); - LOG_DEBUG << "refitTrackWithSiHits failed refit" << endm; + LOG_DEBUG << "refitTrackWithFstHits failed refit" << endm; } else { mHist["FitStatus"]->Fill("GoodReFit", 1); - gtr.setFst( si_hits_for_this_track, mTrackFitter->getTrack() ); + gtr.addFST( fst_hits_for_this_track, mTrackFitter->getTrack() ); } } - mNumFstHits[i] = nSiHitsFound; + mNumFstHits[i] = nFstHitsFound; mFitMoms[i] = p; } // we have 3 Si hits to refit with if ( mGenHistograms ){ - mHist["FitStatus"]->Fill( TString::Format( "w%luSi", nSiHitsFound ).Data(), 1 ); + mHist["FitStatus"]->Fill( TString::Format( "w%luSi", nFstHitsFound ).Data(), 1 ); } } // loop on the global tracks } // ad Si hits via MC associations - void addSiHits() { + /** + * @brief Adds compatible FTT hits to tracks seeded with FST + * + */ + void addFttHits() { + FwdDataSource::HitMap_t hitmap = mDataSource->getFttHits(); + + for ( auto &t : mTrackResults ){ + if (t.isFitConvergedFully != true) continue; + + Seed_t hits_near_plane0; + Seed_t hits_near_plane1; + Seed_t hits_near_plane2; + Seed_t hits_near_plane3; + Seed_t hits_to_add; + try { + auto msp3 = mTrackFitter->projectToFtt(3, t.track); + auto msp2 = mTrackFitter->projectToFtt(2, t.track); + auto msp1 = mTrackFitter->projectToFtt(1, t.track); + auto msp0 = mTrackFitter->projectToFtt(0, t.track); + + // now look for Ftt hits near these + hits_near_plane3 = findFttHitsNearProjectedState(hitmap[3], msp3); + hits_near_plane2 = findFttHitsNearProjectedState(hitmap[2], msp2); + hits_near_plane1 = findFttHitsNearProjectedState(hitmap[1], msp1); + hits_near_plane0 = findFttHitsNearProjectedState(hitmap[0], msp0); + + LOG_INFO << " Found #FTT hits on planes " << TString::Format( "[%d, %d, %d, %d]", hits_near_plane0.size(), hits_near_plane1.size(), hits_near_plane2.size(), hits_near_plane3.size() ) << endm; + } catch (genfit::Exception &e) { + // Failed to project to Si disk: ", e.what() + LOG_WARN << "Unable to get Ftt projections: " << e.what() << endm; + } + + if ( mGenHistograms ){ + this->mHist[ "nFttHitsFound" ]->Fill( 1, hits_near_plane0.size() ); + this->mHist[ "nFttHitsFound" ]->Fill( 2, hits_near_plane1.size() ); + this->mHist[ "nFttHitsFound" ]->Fill( 3, hits_near_plane2.size() ); + this->mHist[ "nFttHitsFound" ]->Fill( 4, hits_near_plane3.size() ); + } + + if ( hits_near_plane0.size() > 0 ) + hits_to_add.push_back( hits_near_plane0[0] ); + if ( hits_near_plane1.size() > 0 ) + hits_to_add.push_back( hits_near_plane1[0] ); + if ( hits_near_plane2.size() > 0 ) + hits_to_add.push_back( hits_near_plane2[0] ); + if ( hits_near_plane3.size() > 0 ) + hits_to_add.push_back( hits_near_plane3[0] ); + + if ( hits_to_add.size() > 0 ){ + + Seed_t combinedSeed; + LOG_INFO << "Found " << t.fstSeed.size() << " existing FST seed points" << endm; + combinedSeed.insert( combinedSeed.begin(), t.fstSeed.begin(), t.fstSeed.end() ); // this is goofed but will fix + LOG_INFO << "Adding " << hits_to_add.size() << " new FTT seed points" << endm; + combinedSeed.insert( combinedSeed.end(), hits_to_add.begin(), hits_to_add.end() ); + + double vertex[3] = { mEventVertex.X(), mEventVertex.Y(), mEventVertex.Z() }; + LOG_INFO << "Using previous fit momentum as the seed: " << TString::Format( "(pt=%f, eta=%f, phi=%f)", t.momentum.Pt(), t.momentum.Eta(), t.momentum.Phi() ) << endm; + + TVector3 p = mTrackFitter->fitTrack(combinedSeed, vertex, &(t.momentum) ); + + auto status = mTrackFitter->getTrack()->getFitStatus(); + if ( status->isFitConvergedFully() ){ + LOG_INFO << "Track Refit with FTT points converged" << endm; + t.addFTT( hits_to_add, mTrackFitter->getTrack() ); + LOG_INFO << "Track Refit with " << t.track->getNumPoints() << " points" << endm; + } else { + LOG_INFO << "Track Refit with FTT points FAILED" << endm; + } + } + + if ( mGenHistograms ){ + mHist["FitStatus"]->Fill( TString::Format( "w%luFtt", hits_to_add.size() ).Data(), 1 ); + } + + + } // loop over tracks + } + + /** + * @brief Adds compatible FTT hits using MC info + * + */ + void addFttHitsMc() { + FwdDataSource::HitMap_t hitmap = mDataSource->getFttHits(); + + + for (size_t i = 0; i < mTrackResults.size(); i++) { + GenfitTrackResult >r = mTrackResults[i]; + + if ( gtr.status.isFitConverged() == false || gtr.momentum.Perp() < 1e-3) { + LOG_DEBUG << "Skipping addFttHitsMc on this track, fit failed" << endm; + return; + } + + if ( mGenHistograms){ + mHist["FitStatus"]->Fill("PossibleReFit", 1); + } + + Seed_t fttHitsForThisTrack; + + size_t nFttHitsFound = 0; + for (size_t j = 0; j < 4; j++) { + for (auto h0 : hitmap[j]) { + if (dynamic_cast(h0)->_tid == gtr.track->getMcTrackId()) { + // fttHitsForThisTrack[j] = h0; + fttHitsForThisTrack.push_back( h0 ); + nFttHitsFound++; + if ( mGenHistograms ){ + this->mHist[ "nFttHitsFound" ]->Fill( j+1, ( fttHitsForThisTrack[j] != nullptr ? 1 : 0 ) ); + } + break; + } + } // loop on hits in this layer of hitmap + } // loop on hitmap layers + + LOG_DEBUG << "Found " << nFttHitsFound << " FTT Hits on this track (MC lookup)" << endm; + + if (nFttHitsFound >= 1) { + if ( mGenHistograms ){ + mHist["FitStatus"]->Fill("AttemptReFit", 1); + } + + Seed_t combinedSeed; + LOG_INFO << "Found " << gtr.fstSeed.size() << " existing FST seed points" << endm; + combinedSeed.insert( combinedSeed.begin(), gtr.fstSeed.begin(), gtr.fstSeed.end() ); // this is goofed but will fix + LOG_INFO << "Adding " << fttHitsForThisTrack.size() << " new FTT seed points" << endm; + combinedSeed.insert( combinedSeed.end(), fttHitsForThisTrack.begin(), fttHitsForThisTrack.end() ); + + double vertex[3] = { mEventVertex.X(), mEventVertex.Y(), mEventVertex.Z() }; + LOG_INFO << "Using previous fit momentum as the seed: " << TString::Format( "(pt=%f, eta=%f, phi=%f)", gtr.momentum.Pt(), gtr.momentum.Eta(), gtr.momentum.Phi() ) << endm; + + LOG_INFO << "STarting the Fit" << endm; + TVector3 p = mTrackFitter->fitTrack(combinedSeed, vertex, &(gtr.momentum) ); + LOG_INFO << "Fit Complete" << endm; + + auto status = mTrackFitter->getTrack()->getFitStatus(); + if ( status && status->isFitConvergedFully() ){ + LOG_INFO << "Track Refit with FTT points converged" << endm; + gtr.addFTT( fttHitsForThisTrack, mTrackFitter->getTrack() ); + LOG_INFO << "Track Refit with " << gtr.track->getNumPoints() << " points" << endm; + } else { + LOG_INFO << "Track Refit with FTT points FAILED" << endm; + } + + } // we have 3 Si hits to refit with + + if ( mGenHistograms ){ + mHist["FitStatus"]->Fill( TString::Format( "w%luFtt", nFttHitsFound ).Data(), 1 ); + } + } // loop on the global tracks + } // add Ftt hits via MC associations + + /** + * @brief Adds compatible FST hits to a track seeded from FTT + * + */ + void addFstHits() { FwdDataSource::HitMap_t hitmap = mDataSource->getFstHits(); // loop on global tracks @@ -1120,52 +1336,53 @@ class ForwardTrackMaker { auto msp0 = mTrackFitter->projectToFst(0, mGlobalTracks[i]); // now look for Si hits near these - hits_near_disk2 = findSiHitsNearMe(hitmap[2], msp2); - hits_near_disk1 = findSiHitsNearMe(hitmap[1], msp1); - hits_near_disk0 = findSiHitsNearMe(hitmap[0], msp0); + hits_near_disk2 = findFstHitsNearProjectedState(hitmap[2], msp2); + hits_near_disk1 = findFstHitsNearProjectedState(hitmap[1], msp1); + hits_near_disk0 = findFstHitsNearProjectedState(hitmap[0], msp0); } catch (genfit::Exception &e) { // Failed to project to Si disk: ", e.what() + LOG_WARN << "Unable to get Ftt projections: " << e.what() << endm; } vector hits_to_add; - size_t nSiHitsFound = 0; // this is really # of disks on which a hit is found + size_t nFstHitsFound = 0; // this is really # of disks on which a hit is found if ( mGenHistograms ){ - this->mHist[ "nSiHitsFound" ]->Fill( 1, hits_near_disk0.size() ); - this->mHist[ "nSiHitsFound" ]->Fill( 2, hits_near_disk1.size() ); - this->mHist[ "nSiHitsFound" ]->Fill( 3, hits_near_disk2.size() ); + this->mHist[ "nFstHitsFound" ]->Fill( 1, hits_near_disk0.size() ); + this->mHist[ "nFstHitsFound" ]->Fill( 2, hits_near_disk1.size() ); + this->mHist[ "nFstHitsFound" ]->Fill( 3, hits_near_disk2.size() ); } // TODO: HANDLE multiple points found? - if ( hits_near_disk0.size() == 1 ) { + if ( hits_near_disk0.size() >= 1 ) { hits_to_add.push_back( hits_near_disk0[0] ); - nSiHitsFound++; + nFstHitsFound++; } else { hits_to_add.push_back( nullptr ); } - if ( hits_near_disk1.size() == 1 ) { + if ( hits_near_disk1.size() >= 1 ) { hits_to_add.push_back( hits_near_disk1[0] ); - nSiHitsFound++; + nFstHitsFound++; } else { hits_to_add.push_back( nullptr ); } - if ( hits_near_disk2.size() == 1 ) { + if ( hits_near_disk2.size() >= 1 ) { hits_to_add.push_back( hits_near_disk2[0] ); - nSiHitsFound++; + nFstHitsFound++; } else { hits_to_add.push_back( nullptr ); } - if (nSiHitsFound >= 1) { + if (nFstHitsFound >= 1) { if ( mGenHistograms ){ mHist["FitStatus"]->Fill("AttemptReFit", 1); } - // LOG_INFO << "Fitting on GlobalTrack : " << mGlobalTracks[i] << " with " << nSiHitsFound << " si hits" << endm; - TVector3 p = mTrackFitter->refitTrackWithSiHits(mGlobalTracks[i], hits_to_add); + // LOG_INFO << "Fitting on GlobalTrack : " << mGlobalTracks[i] << " with " << nFstHitsFound << " si hits" << endm; + TVector3 p = mTrackFitter->refitTrackWithFstHits(mGlobalTracks[i], hits_to_add); size_t lengthGTR = mTrackResults.size(); if ( lengthGTR >= 1 ){ - mTrackResults[ lengthGTR - 1 ].setFst( hits_to_add, mTrackFitter->getTrack() ); + mTrackResults[ lengthGTR - 1 ].addFST( hits_to_add, mTrackFitter->getTrack() ); } else { LOG_ERROR << "Fit Results not found" << endm; } @@ -1178,8 +1395,7 @@ class ForwardTrackMaker { } } - // mGlobalTracks[i] = mTrackFitter->getTrack(); - mNumFstHits[i] = nSiHitsFound; + mNumFstHits[i] = nFstHitsFound; mFitMoms[i] = p; } else { @@ -1187,13 +1403,22 @@ class ForwardTrackMaker { } if ( mGenHistograms ){ - mHist["FitStatus"]->Fill( TString::Format( "w%luSi", nSiHitsFound ).Data(), 1 ); + mHist["FitStatus"]->Fill( TString::Format( "w%luSi", nFstHitsFound ).Data(), 1 ); } } // loop on globals - } // addSiHits - - Seed_t findSiHitsNearMe(Seed_t &available_hits, genfit::MeasuredStateOnPlane &msp, double dphi = 0.004 * 9.5, double dr = 2.75) { + } // addFstHits + + /** + * @brief Finds FST hits near projected state + * + * @param available_hits : FST hits to consider + * @param msp : measured state on plabe from existing track projection + * @param dphi : search distance in phi + * @param dr : search distance in r + * @return Seed_t : compatible FST hits + */ + Seed_t findFstHitsNearProjectedState(Seed_t &available_hits, genfit::MeasuredStateOnPlane &msp, double dphi = 0.004 * 20.5, double dr = 2.75 * 2) { double probe_phi = TMath::ATan2(msp.getPos().Y(), msp.getPos().X()); double probe_r = sqrt(pow(msp.getPos().X(), 2) + pow(msp.getPos().Y(), 2)); @@ -1203,6 +1428,8 @@ class ForwardTrackMaker { double h_phi = TMath::ATan2(h->getY(), h->getX()); double h_r = sqrt(pow(h->getX(), 2) + pow(h->getY(), 2)); double mdphi = fabs(h_phi - probe_phi); + if (mdphi > 2*3.1415926) + mdphi = mdphi - 2*3.1415926; if ( mdphi < dphi && fabs( h_r - probe_r ) < dr) { // handle 2pi edge found_hits.push_back(h); @@ -1212,6 +1439,74 @@ class ForwardTrackMaker { return found_hits; } + /** + * @brief Finds FTT hits near projected state + * + * @param available_hits : FTT hits to consider + * @param msp : measured state on plane from existing track fit projection + * @param dx : search distance in x + * @param dy : search distance in y + * + * @return compatible FTT hits + */ + Seed_t findFttHitsNearProjectedState(Seed_t &available_hits, genfit::MeasuredStateOnPlane &msp, double dx = 1.5, double dy = 1.5) { + + Seed_t found_hits; + double probe_phi = TMath::ATan2(msp.getPos().Y(), msp.getPos().X()); + double probe_r = sqrt(pow(msp.getPos().X(), 2) + pow(msp.getPos().Y(), 2)); + + TLorentzVector lv1, lv2; + lv1.SetPxPyPzE( msp.getPos().X(), msp.getPos().Y(), 0, 1 ); + + double mindx = 99; + double mindy = 99; + double mindr = 99; + double mindp = 99; + KiTrack::IHit *closest = nullptr; + + for (auto h : available_hits) { + + lv2.SetPxPyPzE( h->getX(), h->getY(), 0, 1 ); + double sr = lv1.Pt() - lv2.Pt(); + double sp = lv1.DeltaPhi( lv2 ); + double sx = h->getX() - msp.getPos().X(); + double sy = h->getY() - msp.getPos().Y(); + + if ( fabs(sr) < fabs(mindr) ) + mindr = sr; + if ( fabs(sp) < fabs(mindp) ){ + mindp = sp; + closest = h; + } + if ( fabs(sx) < fabs(mindx) ) + mindx = sx; + if ( fabs(sy) < fabs(mindy) ) + mindy = sy; + + if ( mGenHistograms ){ + mHist[ "FttSearchDeltaX" ]->Fill( sx ); + mHist[ "FttSearchDeltaY" ]->Fill( sy ); + mHist[ "FttSearchDeltaR" ]->Fill( sr ); + mHist[ "FttSearchDeltaP" ]->Fill( sp ); + } + + + } + + if ( mGenHistograms ){ + mHist[ "FttSearchMinDeltaX" ]->Fill( mindx ); + mHist[ "FttSearchMinDeltaY" ]->Fill( mindy ); + mHist[ "FttSearchMinDeltaR" ]->Fill( mindr ); + mHist[ "FttSearchMinDeltaP" ]->Fill( mindp ); + } + + if ( fabs(mindp) < 0.04*5 && fabs(mindr) < 9 ) { + found_hits.push_back(closest); + } + + return found_hits; + } + bool getSaveCriteriaValues() { return mSaveCriteriaValues; } std::vector getTwoHitCriteria() { return mTwoHitCrit; } std::vector getThreeHitCriteria() { return mThreeHitCrit; } diff --git a/StRoot/StFwdTrackMaker/include/Tracker/QualityPlotter.h b/StRoot/StFwdTrackMaker/include/Tracker/QualityPlotter.h index 973945a3994..a81a97ff410 100644 --- a/StRoot/StFwdTrackMaker/include/Tracker/QualityPlotter.h +++ b/StRoot/StFwdTrackMaker/include/Tracker/QualityPlotter.h @@ -217,12 +217,12 @@ class QualityPlotter { if (kv.second == nullptr) continue; - this->get("nHitsOnTrackMc")->Fill(kv.second->mHits.size()); + this->get("nHitsOnTrackMc")->Fill(kv.second->mFttHits.size()); this->get("McPt")->Fill(kv.second->mPt); this->get("McEta")->Fill(kv.second->mEta); this->get("McPhi")->Fill(kv.second->mPhi); - if (kv.second->mHits.size() >= 4) { + if (kv.second->mFttHits.size() >= 4) { this->get("McPt_4hits")->Fill(kv.second->mPt); this->get("McEta_4hits")->Fill(kv.second->mEta); this->get("McPhi_4hits")->Fill(kv.second->mPhi); @@ -231,25 +231,25 @@ class QualityPlotter { ((TH3 *)this->get("McPtEtaPhi_4hits"))->Fill(kv.second->mPt, kv.second->mEta, kv.second->mPhi); } - if (kv.second->mHits.size() >= 5) { + if (kv.second->mFttHits.size() >= 5) { this->get("McPt_5hits")->Fill(kv.second->mPt); this->get("McEta_5hits")->Fill(kv.second->mEta); this->get("McPhi_5hits")->Fill(kv.second->mPhi); } - if (kv.second->mHits.size() >= 6) { + if (kv.second->mFttHits.size() >= 6) { this->get("McPt_6hits")->Fill(kv.second->mPt); this->get("McEta_6hits")->Fill(kv.second->mEta); this->get("McPhi_6hits")->Fill(kv.second->mPhi); } - if (kv.second->mHits.size() >= 7) { + if (kv.second->mFttHits.size() >= 7) { this->get("McPt_7hits")->Fill(kv.second->mPt); this->get("McEta_7hits")->Fill(kv.second->mEta); this->get("McPhi_7hits")->Fill(kv.second->mPhi); } - for (auto h : kv.second->mHits) { + for (auto h : kv.second->mFttHits) { auto fh = static_cast(h); this->get("McHitMap")->Fill(abs(fh->_vid)); std::string n = "McHitMapLayer" + std::to_string(fh->getLayer()); diff --git a/StRoot/StFwdTrackMaker/include/Tracker/TrackFitter.h b/StRoot/StFwdTrackMaker/include/Tracker/TrackFitter.h index de55a2a53e1..5475b8cfb1a 100644 --- a/StRoot/StFwdTrackMaker/include/Tracker/TrackFitter.h +++ b/StRoot/StFwdTrackMaker/include/Tracker/TrackFitter.h @@ -55,14 +55,26 @@ class TrackFitter { public: - // ctor - // provide the main configuration object + + /** + * @brief Construct a new Track Fitter object + * + * @param _mConfig : Config object + * @param geoCache : Geometry cache filename + */ TrackFitter(FwdTrackerConfig _mConfig, TString geoCache) : mConfig(_mConfig), mGeoCache(geoCache) { mTrackRep = 0; mFitTrack = 0; } - + /** + * @brief Setup the tracker object + * Load geometry + * Setup Material Effects + * Setup the magnetic field + * Setup the fitter + * Setup the fit planes + */ void setup() { // the geometry manager that GenFit will use @@ -223,6 +235,10 @@ class TrackFitter { makeHistograms(); } + /** + * @brief Prepare QA histograms + * + */ void makeHistograms() { std::string n = ""; mHist["ECalProjPosXY"] = new TH2F("ECalProjPosXY", ";X;Y", 1000, -500, 500, 1000, -500, 500); @@ -278,7 +294,10 @@ class TrackFitter { mHist[n] = new TH1F(n.c_str(), "; Duraton (ms)", 500, 0, 50000); } - // writes mHistograms stored in map only if mGenHistograms is true + /** + * @brief writes histograms stored in map only if mGenHistograms is true + * + */ void writeHistograms() { if ( !mGenHistograms ) return; @@ -288,9 +307,12 @@ class TrackFitter { } } - /* Convert the 3x3 covmat to 2x2 by dropping z - * - */ + /** + * @brief Convert the 3x3 covmat to 2x2 by dropping z + * + * @param h : hit with cov matrix + * @return TMatrixDSym : cov matrix 2x2 + */ TMatrixDSym CovMatPlane(KiTrack::IHit *h){ TMatrixDSym cm(2); cm(0, 0) = static_cast(h)->_covmat(0, 0); @@ -300,12 +322,16 @@ class TrackFitter { } - /* FitSimpleCircle - * Used to determine a seed transverse momentum based on space points - * Takes a list of space points KiTrack::IHit * - * Takes three indecise used to lookup three of the possible hits within the list - */ - float fitSimpleCircle(Seed_t trackCand, size_t i0, size_t i1, size_t i2) { + /** + * @brief Fit points to a simple circle + * + * @param trackSeed : seed points to fit + * @param i0 : index of the first hit + * @param i1 : index of the second hit + * @param i2 : index of the third hit + * @return float : curvature + */ + float fitSimpleCircle(Seed_t trackSeed, size_t i0, size_t i1, size_t i2) { float curv = 0; // ensure that no index is outside of range for FST or FTT volumes @@ -313,61 +339,80 @@ class TrackFitter { return 0; try { - KiTrack::SimpleCircle sc(trackCand[i0]->getX(), trackCand[i0]->getY(), trackCand[i1]->getX(), trackCand[i1]->getY(), trackCand[i2]->getX(), trackCand[i2]->getY()); + KiTrack::SimpleCircle sc(trackSeed[i0]->getX(), trackSeed[i0]->getY(), trackSeed[i1]->getX(), trackSeed[i1]->getY(), trackSeed[i2]->getX(), trackSeed[i2]->getY()); curv = sc.getRadius(); } catch (KiTrack::InvalidParameter &e) { // if we got here we failed to get a valid seed. We will still try to move forward but the fit will probably fail + LOG_WARN << "Circle fit failed, FWD track fit will likely faile" << endm; } // make sure the curv is valid - if (isinf(curv)) + if (isinf(curv)){ curv = 999999.9; + } return curv; - } - - /* seedState - * Determines the seed position and momentum for a list of space points + } // fitSimpleCircle + + /** + * @brief Determines the seed state to start the fit + * + * @param trackSeed : seed points + * @param seedPos : position seed (return) + * @param seedMom : momenum seed (return) + * @return float : curvature */ - float seedState(Seed_t trackCand, TVector3 &seedPos, TVector3 &seedMom) { - // we require at least 4 hits, so this should be gauranteed - if(trackCand.size() < 3){ + float seedState(Seed_t trackSeed, TVector3 &seedPos, TVector3 &seedMom) { + // we require at least 3 hits, so this should be gauranteed + if(trackSeed.size() < 3){ // failure return 0.0; } - // we want to use the LAST 3 hits, since silicon doesnt have R information TVector3 p0, p1, p2; // use the closest hit to the interaction point for the seed pos - FwdHit *hit_closest_to_IP = static_cast(trackCand[0]); + FwdHit *hit_closest_to_IP = static_cast(trackSeed[0]); - // maps from to - std::map vol_map; + // maps from to + std::map vol_map; // init the map for (size_t i = 0; i < 13; i++) vol_map[i] = -1; - for (size_t i = 0; i < trackCand.size(); i++) { - auto fwdHit = static_cast(trackCand[i]); - vol_map[abs(fwdHit->_vid)] = i; + vector idx; + for (size_t i = 0; i < trackSeed.size(); i++) { + auto fwdHit = static_cast(trackSeed[i]); + if (vol_map[ abs(fwdHit->_vid) ] == -1) + idx.push_back(fwdHit->_vid); + vol_map[abs(fwdHit->_vid)] = (int)i; + // find the hit closest to IP for the initial position seed if (hit_closest_to_IP->getZ() > fwdHit->getZ()) hit_closest_to_IP = fwdHit; } // now get an estimate of the pT from several overlapping simple circle fits - // enumerate the available partitions + // enumerate the available partitions (example for FTT) // 12 11 10 // 12 11 9 // 12 10 9 // 11 10 9 vector curvs; - curvs.push_back(fitSimpleCircle(trackCand, vol_map[12], vol_map[11], vol_map[10])); - curvs.push_back(fitSimpleCircle(trackCand, vol_map[12], vol_map[11], vol_map[9])); - curvs.push_back(fitSimpleCircle(trackCand, vol_map[12], vol_map[10], vol_map[9])); - curvs.push_back(fitSimpleCircle(trackCand, vol_map[11], vol_map[10], vol_map[9])); + + if (idx.size() < 3){ + return 0.0; + } + + if ( idx.size() == 3 ){ + curvs.push_back(fitSimpleCircle(trackSeed, vol_map[idx[0]], vol_map[idx[1]], vol_map[idx[2]])); + } else if ( idx.size() >= 4 ){ + curvs.push_back(fitSimpleCircle(trackSeed, vol_map[idx[0]], vol_map[idx[1]], vol_map[idx[2]])); + curvs.push_back(fitSimpleCircle(trackSeed, vol_map[idx[0]], vol_map[idx[1]], vol_map[idx[3]])); + curvs.push_back(fitSimpleCircle(trackSeed, vol_map[idx[0]], vol_map[idx[2]], vol_map[idx[3]])); + curvs.push_back(fitSimpleCircle(trackSeed, vol_map[idx[1]], vol_map[idx[2]], vol_map[idx[3]])); + } // average them and exclude failed fits float mcurv = 0; @@ -375,7 +420,7 @@ class TrackFitter { for (size_t i = 0; i < curvs.size(); i++) { if (mGenHistograms) - this->mHist["seed_curv"]->Fill(curvs[i]); + mHist["seed_curv"]->Fill(curvs[i]); if (curvs[i] > 10) { mcurv += curvs[i]; nmeas += 1.0; @@ -385,15 +430,16 @@ class TrackFitter { if (nmeas >= 1) mcurv = mcurv / nmeas; else - mcurv = 10; + mcurv = 100; // Now lets get eta information - // simpler, use farthest points from IP - if (vol_map[9] < 13) - p0.SetXYZ(trackCand[vol_map[9]]->getX(), trackCand[vol_map[9]]->getY(), trackCand[vol_map[9]]->getZ()); + p0.SetXYZ(trackSeed[vol_map[idx[0]]]->getX(), trackSeed[vol_map[idx[0]]]->getY(), trackSeed[vol_map[idx[0]]]->getZ()); + p1.SetXYZ(trackSeed[vol_map[idx[1]]]->getX(), trackSeed[vol_map[idx[1]]]->getY(), trackSeed[vol_map[idx[1]]]->getZ()); + if ( abs(p0.X() - p1.X()) < 1e-6 ){ + p1.SetXYZ(trackSeed[vol_map[idx[2]]]->getX(), trackSeed[vol_map[idx[2]]]->getY(), trackSeed[vol_map[idx[2]]]->getZ()); + } - if (vol_map[10] < 13) - p1.SetXYZ(trackCand[vol_map[10]]->getX(), trackCand[vol_map[10]]->getY(), trackCand[vol_map[10]]->getZ()); + LOG_DEBUG << TString::Format( "Fwd SeedState: p0 (%f, %f, %f), p1 (%f, %f, %f)", p0.X(), p0.Y(), p0.Z(), p1.X(), p1.Y(), p1.Z() ) << endm; const double K = 0.00029979; //K depends on the units used for Bfield double pt = mcurv * K * 5; // pT from average measured curv @@ -403,8 +449,15 @@ class TrackFitter { double phi = TMath::ATan2(dy, dx); double Rxy = sqrt(dx * dx + dy * dy); double theta = TMath::ATan2(Rxy, dz); + if (abs(dx) < 1e-6 || abs(dy) < 1e-6){ + phi = TMath::ATan2( p1.Y(), p1.X() ); + } + Rxy = sqrt( p0.X()*p0.X() + p0.Y()*p0.Y() ); + theta = TMath::ATan2(Rxy, p0.Z()); + + LOG_DEBUG << TString::Format( "pt=%f, dx=%f, dy=%f, dz=%f, phi=%f, theta=%f", pt, dx, dy, dz, phi, theta ) << endm; // double eta = -log( tantheta / 2.0 ); - // these starting conditions can probably be improvd, good study for student + // these starting conditions can probably be improvd, good study for students seedMom.SetPtThetaPhi(pt, theta, phi); seedPos.SetXYZ(hit_closest_to_IP->getX(), hit_closest_to_IP->getY(), hit_closest_to_IP->getZ()); @@ -415,34 +468,57 @@ class TrackFitter { } return mcurv; - } + }//seedState - genfit::MeasuredStateOnPlane projectToFst(size_t si_plane, genfit::Track *fitTrack) { - if (si_plane > 2) { + /** + * @brief Get projection to given FST plane + * + * @param fstPlane : plane index + * @param fitTrack : track to project + * @return genfit::MeasuredStateOnPlane + */ + genfit::MeasuredStateOnPlane projectToFst(size_t fstPlane, genfit::Track *fitTrack) { + if (fstPlane > 2) { genfit::MeasuredStateOnPlane nil; return nil; } - auto detSi = mFSTPlanes[si_plane]; - genfit::MeasuredStateOnPlane tst = fitTrack->getFittedState(1); - auto TCM = fitTrack->getCardinalRep()->get6DCov(tst); - - // this returns the track length if needed - fitTrack->getCardinalRep()->extrapolateToPlane(tst, detSi); - - TCM = fitTrack->getCardinalRep()->get6DCov(tst); + auto detFst = mFSTPlanes[fstPlane]; + // TODO: Why use 1 here? + genfit::MeasuredStateOnPlane tst = fitTrack->getFittedState(1); + // NOTE: this returns the track length if needed + fitTrack->getCardinalRep()->extrapolateToPlane(tst, detFst); - // can get the projected positions if needed - float x = tst.getPos().X(); - float y = tst.getPos().Y(); - float z = tst.getPos().Z(); - // and the uncertainties - LOG_DEBUG << "Track Uncertainty at FST (plane=" << si_plane << ") @ x= " << x << ", y= " << y << ", z= " << z << " : " << sqrt(TCM(0, 0)) << ", " << sqrt(TCM(1, 1)) << endm; + return tst; + } + /** + * @brief Get projection to given FTT plane + * + * @param fttPlane : plane index + * @param fitTrack : track to project + * @return genfit::MeasuredStateOnPlane + */ + genfit::MeasuredStateOnPlane projectToFtt(size_t iFttPlane, genfit::Track *fitTrack) { + if (iFttPlane > 3) { + genfit::MeasuredStateOnPlane nil; + return nil; + } + auto fttPlane = mFTTPlanes[iFttPlane]; + // TODO: why use 1 here? + genfit::MeasuredStateOnPlane tst = fitTrack->getFittedState(1); + // NOTE: this returns the track length if needed + fitTrack->getCardinalRep()->extrapolateToPlane(tst, fttPlane); return tst; } + /** + * @brief Get the Fst Plane object for a given hit + * + * @param h : hit + * @return genfit::SharedPlanePtr + */ genfit::SharedPlanePtr getFstPlane( FwdHit * h ){ size_t planeId = h->getSector(); @@ -467,18 +543,18 @@ class TrackFitter { } return planeCorr; + } // GetFST PLANE - } - - /* RefitTracksWithSiHits - * Takes a previously fit track re-fits it with the newly added silicon hits + /** + * @brief Refit a track with additional FST hits * + * Takes a previously fit track re-fits it with the newly added silicon hits + * @param originalTrack : original fit track + * @param fstHits : new FST hits to add + * @return TVector3 : momentum */ - TVector3 refitTrackWithSiHits(genfit::Track *originalTrack, Seed_t si_hits) { - // mem leak, global track is overwritten without delete. + TVector3 refitTrackWithFstHits(genfit::Track *originalTrack, Seed_t fstHits) { TVector3 pOrig = originalTrack->getCardinalRep()->getMom(originalTrack->getFittedState(1, originalTrack->getCardinalRep())); - - // auto cardinalStatus = originalTrack->getFitStatus(originalTrack->getCardinalRep()); if (originalTrack->getFitStatus(originalTrack->getCardinalRep())->isFitConverged() == false) { // in this case the original track did not converge so we should not refit. @@ -528,7 +604,7 @@ class TrackFitter { int hitId(5); // add the hits to the track - for (auto h : si_hits) { + for (auto h : fstHits) { if ( nullptr == h ) continue; // if no Si hit in this plane, skip hitCoords[0] = h->getX(); @@ -542,7 +618,6 @@ class TrackFitter { return pOrig; } - // auto plane = mFSTPlanes[planeId]; auto plane = getFstPlane( static_cast(h) ); measurement->setPlane(plane, planeId); @@ -566,7 +641,6 @@ class TrackFitter { double cdz = (h->getZ() - planeCorr->getO().Z()); if (mGenHistograms){ - ((TH2*)mHist[ "FstDiffZVsR" ])->Fill( r, dz ); if ( r < 16 ) {// inner @@ -577,10 +651,8 @@ class TrackFitter { mHist["CorrFstDiffZVsPhiSliceOuter"]->Fill( phi_slice, cdz ); mHist["FstDiffZVsPhiOuter"]->Fill( phi, dz ); } - } - // mHist[ "FstDiffZVsPhiSliceInner" ]->Fill( sqrt( pow(hitXYZ.x(), 2), pow(hitXYZ.y(), 2) ), dz ); - - } + } // gen histograms + } // for fstHits // start at 0 if PV not included, 1 otherwise for (size_t i = firstFTTIndex; i < trackPoints.size(); i++) { // clone the track points into this track @@ -700,10 +772,13 @@ class TrackFitter { } //refitwith GBL - - /* Generic method for fitting space points with GenFit - * + /** + * @brief Generic method for fitting space points with GenFit * + * @param spoints : spacepoints + * @param seedPos : seed position + * @param seedMom : seed momentum + * @return TVector3 : momentum from fit */ TVector3 fitSpacePoints( vector spoints, TVector3 &seedPos, TVector3 &seedMom ){ @@ -745,11 +820,15 @@ class TrackFitter { return TVector3(0, 0, 0); } - /* Fit a track - * + /** + * @brief Primary track fitting routine * + * @param trackSeed : + * @param Vertex : Primary Vertex + * @param seedMomentum : seed momentum (can be from MC) + * @return TVector3 : fit momentum */ - TVector3 fitTrack(Seed_t trackCand, double *Vertex = 0, TVector3 *McSeedMom = 0) { + TVector3 fitTrack(Seed_t trackSeed, double *Vertex = 0, TVector3 *seedMomentum = 0) { long long itStart = FwdTrackerUtils::nowNanoSecond(); if (mGenHistograms) this->mHist["FitStatus"]->Fill("Total", 1); @@ -770,10 +849,10 @@ class TrackFitter { // get the seed info from our hits TVector3 seedMom, seedPos; // returns track curvature if needed - seedState(trackCand, seedPos, seedMom); + seedState(trackSeed, seedPos, seedMom); - if (McSeedMom != nullptr) { - seedMom = *McSeedMom; + if (seedMomentum != nullptr) { + seedMom = *seedMomentum; } // If we use the PV, use that as the start pos for the track @@ -828,8 +907,9 @@ class TrackFitter { /****************************************************************************************************************** * loop over the hits, add them to the track ******************************************************************************************************************/ - for (auto h : trackCand) { - + for (auto h : trackSeed) { + + const bool isFTT = h->getZ() > 200; hitCoords[0] = h->getX(); hitCoords[1] = h->getY(); @@ -837,19 +917,24 @@ class TrackFitter { planeId = h->getSector(); - if (mFTTPlanes.size() <= planeId) { - LOG_WARN << "invalid VolumId -> out of bounds DetPlane, vid = " << planeId << endm; + genfit::SharedPlanePtr plane; + if (isFTT && mFTTPlanes.size() <= planeId) { + LOG_ERROR << "invalid VolumId -> out of bounds DetPlane, vid = " << planeId << endm; return TVector3(0, 0, 0); } - auto plane = mFTTPlanes[planeId]; + if (isFTT) + plane = mFTTPlanes[planeId]; + else + plane = getFstPlane( static_cast(h) ); + measurement->setPlane(plane, planeId); fitTrack.insertPoint(new genfit::TrackPoint(measurement, &fitTrack)); if (abs(h->getZ() - plane->getO().Z()) > 0.05) { LOG_WARN << "Z Mismatch h->z = " << h->getZ() << ", plane->z = "<< plane->getO().Z() <<", diff = " << abs(h->getZ() - plane->getO().Z()) << endm; } - } // loop on trackCand + } // loop on trackSeed /******************************************************************************************************************