diff --git a/Detectors/GlobalTracking/include/GlobalTracking/MatchHMP.h b/Detectors/GlobalTracking/include/GlobalTracking/MatchHMP.h index 0f7e2214220d9..89d04ac4e2e81 100644 --- a/Detectors/GlobalTracking/include/GlobalTracking/MatchHMP.h +++ b/Detectors/GlobalTracking/include/GlobalTracking/MatchHMP.h @@ -37,6 +37,7 @@ #include "DetectorsBase/GeometryManager.h" #include "DataFormatsHMP/Cluster.h" +#include "GlobalTracking/MatchTPCITS.h" #include "DataFormatsTPC/TrackTPC.h" #include "DataFormatsTRD/TrackTRD.h" #include "ReconstructionDataFormats/PID.h" @@ -85,10 +86,20 @@ class MatchHMP void print() const; void printCandidatesHMP() const; - enum DebugFlagTypes : UInt_t { - MatchTreeAll = 0x1 << 1, // ///< produce matching candidates tree for all candidates - }; + ///< set extra time tolerance + void setExtraTimeToleranceTRD(float val) { mExtraTimeToleranceTRD = val; } + ///< get extra tolerance + float getExtraTimeToleranceTRD() const { return mExtraTimeToleranceTRD; } + ///< set extra time tolerance + void setExtraTimeToleranceTOF(float val) { mExtraTimeToleranceTOF = val; } + ///< get extra tolerance + float getExtraTimeToleranceTOF() const { return mExtraTimeToleranceTOF; } + + /* enum DebugFlagTypes : UInt_t { + MatchTreeAll = 0x1 << 1, // ///< produce matching candidates tree for all candidates + }; + */ enum trackType : int8_t { UNCONS = 0, CONSTR, SIZE, @@ -142,13 +153,13 @@ class MatchHMP float mMaxInvPt = 999.; ///< derived from nominal Bz // to be done later - float mTPCTBinMUS = 0.; ///< TPC time bin duration in microseconds - float mTPCTBinMUSInv = 0.; ///< inverse TPC time bin duration in microseconds - float mTPCBin2Z = 0.; ///< conversion coeff from TPC time-bin to Z - float mTimeTolerance = 1e3; ///< tolerance in ns for track-TOF time bracket matching - float mExtraTimeToleranceTRD = 500E3; ///< extra tolerance in ns for track-TOF time bracket matching - float mExtraTimeToleranceTOF = 500E3; ///< extra tolerance in ns for track-TOF time bracket matching - float mSigmaTimeCut = 1.; ///< number of sigmas to cut on time when matching the track to the TOF cluster + float mTPCTBinMUS = 0.; ///< TPC time bin duration in microseconds + float mTPCTBinMUSInv = 0.; ///< inverse TPC time bin duration in microseconds + float mTPCBin2Z = 0.; ///< conversion coeff from TPC time-bin to Z + // float mTimeTolerance = 1e3; ///< tolerance in ns for track-TOF time bracket matching + float mExtraTimeToleranceTRD = 0.; ///< extra tolerance in ns for track-TOF time bracket matching + float mExtraTimeToleranceTOF = 0.; ///< extra tolerance in ns for track-TOF time bracket matching + float mSigmaTimeCut = 3.; ///< number of sigmas to cut on time when matching the track to the TOF cluster static constexpr Double_t BC_TIME = o2::constants::lhc::LHCBunchSpacingNS; // bunch crossing in ns static constexpr Double_t BC_TIME_INV = 1. / BC_TIME; // inv bunch crossing in ns diff --git a/Detectors/GlobalTracking/src/MatchHMP.cxx b/Detectors/GlobalTracking/src/MatchHMP.cxx index a2d6c9b45b76f..a17165e0bfdcb 100644 --- a/Detectors/GlobalTracking/src/MatchHMP.cxx +++ b/Detectors/GlobalTracking/src/MatchHMP.cxx @@ -75,12 +75,14 @@ void MatchHMP::run(const o2::globaltracking::RecoContainer& inp) } for (int it = 0; it < o2::globaltracking::MatchHMP::trackType::SIZE; it++) { - mMatchedTracksIndex[it].clear(); + mTracksIndexCache[it].clear(); if (mMCTruthON) { mTracksLblWork[it].clear(); } } + mHMPTriggersIndexCache.clear(); + bool isPrepareHMPClusters = prepareHMPClusters(); if (!isPrepareHMPClusters) { // check cluster before of tracks to see also if MC is required @@ -95,7 +97,6 @@ void MatchHMP::run(const o2::globaltracking::RecoContainer& inp) } if (mIsITSTPCused || mIsTPCTRDused || mIsITSTPCTRDused || mIsITSTPCTOFused || mIsTPCTOFused || mIsITSTPCTRDTOFused || mIsTPCTRDTOFused) { - // doFastMatching(); doMatching(); } @@ -246,19 +247,37 @@ void MatchHMP::addTPCTOFSeed(const o2::dataformats::TrackTPCTOF& _tr, o2::datafo void MatchHMP::addConstrainedSeed(o2::track::TrackParCov& trc, o2::dataformats::GlobalTrackID srcGID, timeEst timeMUS) { std::array globalPos; + // current track index int it = mTracksWork[o2::globaltracking::MatchHMP::trackType::CONSTR].size(); - // create working copy of track param - mTracksWork[o2::globaltracking::MatchHMP::trackType::CONSTR].emplace_back(std::make_pair(trc, timeMUS)); + auto prop = o2::base::Propagator::Instance(); + float bxyz[3]; + prop->getFieldXYZ(trc.getXYZGlo(), bxyz); + double bz = -bxyz[2]; - mTrackGid[o2::globaltracking::MatchHMP::trackType::CONSTR].emplace_back(srcGID); + float pCut = 0.; - if (mMCTruthON) { - mTracksLblWork[o2::globaltracking::MatchHMP::trackType::CONSTR].emplace_back(mRecoCont->getTPCITSTrackMCLabel(srcGID)); + if (TMath::Abs(bz - 5.0) < 0.5) { + pCut = 0.450; } - mTracksIndexCache[o2::globaltracking::MatchHMP::trackType::CONSTR].push_back(it); + if (TMath::Abs(bz - 2.0) < 0.5) { + pCut = 0.200; + } + + if (trc.getP() > pCut && TMath::Abs(trc.getTgl()) < 0.544 && TMath::Abs(trc.getPhi() - TMath::Pi()) > (TMath::Pi() * 0.5)) { + // create working copy of track param + mTracksWork[o2::globaltracking::MatchHMP::trackType::CONSTR].emplace_back(std::make_pair(trc, timeMUS)); + + mTrackGid[o2::globaltracking::MatchHMP::trackType::CONSTR].emplace_back(srcGID); + + if (mMCTruthON) { + mTracksLblWork[o2::globaltracking::MatchHMP::trackType::CONSTR].emplace_back(mRecoCont->getTPCITSTrackMCLabel(srcGID)); + } + + mTracksIndexCache[o2::globaltracking::MatchHMP::trackType::CONSTR].push_back(it); + } } //______________________________________________ void MatchHMP::addTPCSeed(const o2::tpc::TrackTPC& _tr, o2::dataformats::GlobalTrackID srcGID, float time0, float terr) @@ -305,6 +324,8 @@ bool MatchHMP::prepareHMPClusters() mNumOfTriggers = 0; + mHMPTriggersWork.clear(); + int nTriggersInCurrentChunk = mHMPTriggersArray.size(); LOG(debug) << "nTriggersInCurrentChunk = " << nTriggersInCurrentChunk; mNumOfTriggers += nTriggersInCurrentChunk; @@ -366,7 +387,7 @@ void MatchHMP::doMatching() for (int ievt = 0; ievt < cacheTriggerHMP.size(); ievt++) { // events loop auto& event = mHMPTriggersWork[cacheTriggerHMP[ievt]]; - auto evtTime = o2::InteractionRecord::bc2ns(event.getBc(), event.getOrbit()); // event(trigger) time in ns + auto evtTime = event.getIr().differenceInBCMUS(mStartIR); int evtTracks = 0; @@ -377,41 +398,39 @@ void MatchHMP::doMatching() auto& trefTrk = trackWork.first; prop->getFieldXYZ(trefTrk.getXYZGlo(), bxyz); - Double_t bz = -bxyz[2]; + double bz = -bxyz[2]; - double timeUncert = 100.; //= trackWork.second.getTimeStampError(); + double timeUncert = trackWork.second.getTimeStampError(); - float minTrkTime = (trackWork.second.getTimeStamp() - mSigmaTimeCut * timeUncert) * 1.E3; // minimum track time in ns - float maxTrkTime = (trackWork.second.getTimeStamp() + mSigmaTimeCut * timeUncert) * 1.E3; // maximum track time in ns + float minTrkTime = (trackWork.second.getTimeStamp() - mSigmaTimeCut * timeUncert); + float maxTrkTime = (trackWork.second.getTimeStamp() + mSigmaTimeCut * timeUncert); - if (evtTime < (maxTrkTime + timeFromTF) && evtTime > (minTrkTime + timeFromTF)) { + // if (evtTime < (maxTrkTime + timeFromTF) && evtTime > (minTrkTime + timeFromTF)) { + if (evtTime < maxTrkTime && evtTime > minTrkTime) { evtTracks++; - MatchInfo* matching = new o2::dataformats::MatchInfoHMP(999999, mTrackGid[type][cacheTrk[itrk]]); + MatchInfo matching(999999, mTrackGid[type][cacheTrk[itrk]]); - matching->setHMPIDtrk(0, 0, 0, 0); // no intersection found - matching->setHMPIDmip(0, 0, 0, 0); // store mip info in any case - matching->setIdxHMPClus(99, 99999); // chamber not found, mip not yet considered - matching->setHMPsignal(Recon::kNotPerformed); // ring reconstruction not yet performed - matching->setIdxTrack(trackGid); + matching.setHMPIDtrk(0, 0, 0, 0); // no intersection found + matching.setHMPIDmip(0, 0, 0, 0); // store mip info in any case + matching.setIdxHMPClus(99, 99999); // chamber not found, mip not yet considered + matching.setHMPsignal(Recon::kNotPerformed); // ring reconstruction not yet performed + matching.setIdxTrack(trackGid); - TrackHMP* hmpTrk = new TrackHMP(trefTrk); // create a hmpid track to be used for propagation and matching - TrackHMP* hmpTrkConstrained = nullptr; // create a hmpid track to be used for propagation and matching + TrackHMP hmpTrk(trefTrk); // create a hmpid track to be used for propagation and matching - hmpTrk->set(trefTrk.getX(), trefTrk.getAlpha(), trefTrk.getParams(), trefTrk.getCharge(), trefTrk.getPID()); + hmpTrk.set(trefTrk.getX(), trefTrk.getAlpha(), trefTrk.getParams(), trefTrk.getCharge(), trefTrk.getPID()); double xPc, yPc, xRa, yRa, theta, phi; Int_t iCh = intTrkCha(&trefTrk, xPc, yPc, xRa, yRa, theta, phi, bz); // find the intersected chamber for this track if (iCh < 0) { - delete hmpTrk; - hmpTrk = nullptr; continue; } // no intersection at all, go next track - matching->setHMPIDtrk(xPc, yPc, theta, phi); // store initial infos - matching->setIdxHMPClus(iCh, 9999); // set chamber, index of cluster + cluster size + matching.setHMPIDtrk(xPc, yPc, theta, phi); // store initial infos + matching.setIdxHMPClus(iCh, 9999); // set chamber, index of cluster + cluster size int index = -1; @@ -433,7 +452,7 @@ void MatchHMP::doMatching() oneEventClusters.push_back(cluster); double qthre = pParam->qCut(); - if (cluster.q() < 150.) { + if (cluster.q() < 150. || cluster.size() > 10) { continue; } @@ -441,7 +460,13 @@ void MatchHMP::doMatching() cluLORS[0] = cluster.x(); cluLORS[1] = cluster.y(); // get the LORS coordinates of the cluster - double dist = TMath::Sqrt((xPc - cluLORS[0]) * (xPc - cluLORS[0]) + (yPc - cluLORS[1]) * (yPc - cluLORS[1])); + + double dist = 0.; + + if (TMath::Abs((xPc - cluLORS[0]) * (xPc - cluLORS[0]) + (yPc - cluLORS[1]) * (yPc - cluLORS[1])) > 0.0001) { + + dist = TMath::Sqrt((xPc - cluLORS[0]) * (xPc - cluLORS[0]) + (yPc - cluLORS[1]) * (yPc - cluLORS[1])); + } if (dist < dmin) { dmin = dist; @@ -454,35 +479,27 @@ void MatchHMP::doMatching() // 2. Propagate track to the MIP cluster using the central method if (!bestHmpCluster) { - delete hmpTrk; - hmpTrk = nullptr; - delete hmpTrkConstrained; - hmpTrkConstrained = nullptr; + oneEventClusters.clear(); continue; } - double Dist = TMath::Sqrt((xPc - bestHmpCluster->x()) * (xPc - bestHmpCluster->x()) + (yPc - bestHmpCluster->y()) * (yPc - bestHmpCluster->y())); - TVector3 vG = pParam->lors2Mars(iCh, bestHmpCluster->x(), bestHmpCluster->y()); float gx = vG.X(); float gy = vG.Y(); float gz = vG.Z(); float alpha = TMath::ATan2(gy, gx); float radiusH = TMath::Sqrt(gy * gy + gx * gx); - if (!(hmpTrk->rotate(alpha))) { + if (!(hmpTrk.rotate(alpha))) { continue; } - if (!prop->PropagateToXBxByBz(*hmpTrk, radiusH, o2::base::Propagator::MAX_SIN_PHI, o2::base::Propagator::MAX_STEP, matCorr)) { - delete hmpTrk; - hmpTrk = nullptr; - delete hmpTrkConstrained; - hmpTrkConstrained = nullptr; + if (!prop->PropagateToXBxByBz(hmpTrk, radiusH, o2::base::Propagator::MAX_SIN_PHI, o2::base::Propagator::MAX_STEP, matCorr)) { + oneEventClusters.clear(); continue; } // 3. Update the track with MIP cluster (Improved angular and position resolution - to be used for Cherenkov angle calculation) - o2::track::TrackParCov trackC(*hmpTrk); + o2::track::TrackParCov trackC(hmpTrk); std::array trkPos{0, gz}; std::array trkCov{0.1 * 0.1, 0., 0.1 * 0.1}; @@ -492,39 +509,32 @@ void MatchHMP::doMatching() // 4. Propagate back the constrained track to the radiator radius - hmpTrkConstrained = new TrackHMP(trackC); - hmpTrkConstrained->set(trackC.getX(), trackC.getAlpha(), trackC.getParams(), trackC.getCharge(), trackC.getPID()); - if (!prop->PropagateToXBxByBz(*hmpTrkConstrained, radiusH - kdRadiator, o2::base::Propagator::MAX_SIN_PHI, o2::base::Propagator::MAX_STEP, matCorr)) { - delete hmpTrk; - hmpTrk = nullptr; - delete hmpTrkConstrained; - hmpTrkConstrained = nullptr; + TrackHMP hmpTrkConstrained(trackC); + hmpTrkConstrained.set(trackC.getX(), trackC.getAlpha(), trackC.getParams(), trackC.getCharge(), trackC.getPID()); + if (!prop->PropagateToXBxByBz(hmpTrkConstrained, radiusH - kdRadiator, o2::base::Propagator::MAX_SIN_PHI, o2::base::Propagator::MAX_STEP, matCorr)) { + oneEventClusters.clear(); continue; } - matching->setHmpMom(hmpTrkConstrained->getP()); + float hmpMom = hmpTrkConstrained.getP() * hmpTrkConstrained.getSign(); + + matching.setHmpMom(hmpMom); // 5. Propagation in the last 10 cm with the fast method double xPc0 = 0., yPc0 = 0.; - intTrkCha(iCh, hmpTrkConstrained, xPc0, yPc0, xRa, yRa, theta, phi, bz); + intTrkCha(iCh, &hmpTrkConstrained, xPc0, yPc0, xRa, yRa, theta, phi, bz); // 6. Set match information int cluSize = bestHmpCluster->size(); - matching->setHMPIDmip(bestHmpCluster->x(), bestHmpCluster->y(), bestHmpCluster->q(), 0); // store mip info in any case - matching->setMipClusSize(bestHmpCluster->size()); - matching->setIdxHMPClus(iCh, index + 1000 * cluSize); // set chamber, index of cluster + cluster size - matching->setHMPIDtrk(xPc, yPc, theta, phi); - - matching->setHMPsignal(pParam->kMipQdcCut); + matching.setHMPIDmip(bestHmpCluster->x(), bestHmpCluster->y(), bestHmpCluster->q(), 0); // store mip info in any case + matching.setMipClusSize(bestHmpCluster->size()); + matching.setIdxHMPClus(iCh, index + 1000 * cluSize); // set chamber, index of cluster + cluster size + matching.setHMPIDtrk(xPc, yPc, theta, phi); if (!isOkQcut) { - delete hmpTrk; - hmpTrk = nullptr; - delete hmpTrkConstrained; - hmpTrkConstrained = nullptr; - continue; + matching.setHMPsignal(pParam->kMipQdcCut); } // dmin recalculated @@ -537,7 +547,7 @@ void MatchHMP::doMatching() // isOkDcut = kTRUE; // switch OFF cut if (!isOkDcut) { - matching->setHMPsignal(pParam->kMipDistCut); // closest cluster with enough charge is still too far from intersection + matching.setHMPsignal(pParam->kMipDistCut); // closest cluster with enough charge is still too far from intersection } if (isOkQcut * isOkDcut) { @@ -545,11 +555,8 @@ void MatchHMP::doMatching() } // MIP-Track matched !! if (!isMatched) { - mMatchedTracks[type].push_back(*matching); - delete hmpTrk; - hmpTrk = nullptr; - delete hmpTrkConstrained; - hmpTrkConstrained = nullptr; + mMatchedTracks[type].push_back(matching); + oneEventClusters.clear(); continue; } // If matched continue... @@ -557,21 +564,19 @@ void MatchHMP::doMatching() // 7. Calculate the Cherenkov angle - recon->setImpPC(xPc, yPc); // store track impact to PC - recon->ckovAngle(matching, oneEventClusters, index, nmean, xRa, yRa); // search for Cerenkov angle of this track + recon->setImpPC(xPc, yPc); // store track impact to PC + recon->ckovAngle(&matching, oneEventClusters, index, nmean, xRa, yRa); // search for Cerenkov angle of this track - mMatchedTracks[type].push_back(*matching); + mMatchedTracks[type].push_back(matching); oneEventClusters.clear(); - delete hmpTrk; - hmpTrk = nullptr; - delete hmpTrkConstrained; - hmpTrkConstrained = nullptr; - } // if matching in time } // tracks loop } // events loop + + delete recon; + recon = nullptr; } //================================================================================================================================================== int MatchHMP::intTrkCha(o2::track::TrackParCov* pTrk, double& xPc, double& yPc, double& xRa, double& yRa, double& theta, double& phi, double bz) diff --git a/Detectors/GlobalTrackingWorkflow/include/GlobalTrackingWorkflow/HMPMatcherSpec.h b/Detectors/GlobalTrackingWorkflow/include/GlobalTrackingWorkflow/HMPMatcherSpec.h index ee38399206eb9..ca515b22859a9 100644 --- a/Detectors/GlobalTrackingWorkflow/include/GlobalTrackingWorkflow/HMPMatcherSpec.h +++ b/Detectors/GlobalTrackingWorkflow/include/GlobalTrackingWorkflow/HMPMatcherSpec.h @@ -25,7 +25,7 @@ namespace globaltracking { /// create a processor spec -framework::DataProcessorSpec getHMPMatcherSpec(o2::dataformats::GlobalTrackID::mask_t src, bool useMC); +framework::DataProcessorSpec getHMPMatcherSpec(o2::dataformats::GlobalTrackID::mask_t src, bool useMC, float extratolerancetrd, float extratolerancetof); } // namespace globaltracking } // namespace o2 diff --git a/Detectors/GlobalTrackingWorkflow/src/HMPMatcherSpec.cxx b/Detectors/GlobalTrackingWorkflow/src/HMPMatcherSpec.cxx index 4ea565525bfaa..13095bfb72ec0 100644 --- a/Detectors/GlobalTrackingWorkflow/src/HMPMatcherSpec.cxx +++ b/Detectors/GlobalTrackingWorkflow/src/HMPMatcherSpec.cxx @@ -144,7 +144,7 @@ void HMPMatcherSpec::endOfStream(EndOfStreamContext& ec) mTimer.CpuTime(), mTimer.RealTime(), mTimer.Counter() - 1); } -DataProcessorSpec getHMPMatcherSpec(GID::mask_t src, bool useMC) +DataProcessorSpec getHMPMatcherSpec(GID::mask_t src, bool useMC, float extratolerancetrd, float extratolerancetof) { // uint32_t ss = o2::globaltracking::getSubSpec(strict ? o2::globaltracking::MatchingType::Strict : o2::globaltracking::MatchingType::Standard); auto dataRequest = std::make_shared(); diff --git a/Detectors/GlobalTrackingWorkflow/src/hmp-matcher-workflow.cxx b/Detectors/GlobalTrackingWorkflow/src/hmp-matcher-workflow.cxx index 89a591fb5305d..41a8cee8400e9 100644 --- a/Detectors/GlobalTrackingWorkflow/src/hmp-matcher-workflow.cxx +++ b/Detectors/GlobalTrackingWorkflow/src/hmp-matcher-workflow.cxx @@ -59,6 +59,8 @@ void customize(std::vector& workflowOptions) {"disable-root-input", o2::framework::VariantType::Bool, false, {"disable root-files input reader"}}, {"disable-root-output", o2::framework::VariantType::Bool, false, {"disable root-files output writer"}}, {"track-sources", VariantType::String, std::string{GID::ALL}, {"comma-separated list of sources to use: allowed TPC,ITS-TPC,TPC-TRD,ITS-TPC-TRD (all)"}}, + {"trd-extra-tolerance", o2::framework::VariantType::Float, 0.0f, {"Extra time tolerance for TRD tracks in microsec"}}, + {"tof-extra-tolerance", o2::framework::VariantType::Float, 0.0f, {"Extra time tolerance for TRD tracks in microsec"}}, {"combine-devices", o2::framework::VariantType::Bool, false, {"merge DPL source/writer devices"}}, {"configKeyValues", VariantType::String, "", {"Semicolon separated key=value strings ..."}}}; o2::raw::HBFUtilsInitializer::addConfigOption(options); @@ -77,6 +79,8 @@ WorkflowSpec defineDataProcessing(ConfigContext const& configcontext) // write the configuration used for the workflow auto useMC = !configcontext.options().get("disable-mc"); + auto extratolerancetrd = configcontext.options().get("trd-extra-tolerance"); + auto extratolerancetof = configcontext.options().get("tof-extra-tolerance"); auto disableRootIn = configcontext.options().get("disable-root-input"); auto disableRootOut = configcontext.options().get("disable-root-output"); @@ -114,7 +118,7 @@ WorkflowSpec defineDataProcessing(ConfigContext const& configcontext) } } - specs.emplace_back(o2::globaltracking::getHMPMatcherSpec(src, useMC)); + specs.emplace_back(o2::globaltracking::getHMPMatcherSpec(src, useMC, extratolerancetrd, extratolerancetof)); if (!disableRootOut) { std::vector writers; diff --git a/Detectors/HMPID/reconstruction/src/Recon.cxx b/Detectors/HMPID/reconstruction/src/Recon.cxx index c626d422dd4eb..867e72cc319c8 100644 --- a/Detectors/HMPID/reconstruction/src/Recon.cxx +++ b/Detectors/HMPID/reconstruction/src/Recon.cxx @@ -87,7 +87,9 @@ void Recon::ckovAngle(o2::dataformats::MatchInfoHMP* match, const std::vector 2 * fParam->qCut() || cluster.size() > 4) { + continue; + } double thetaCer, phiCer; if (findPhotCkov(cluster.x(), cluster.y(), thetaCer, phiCer)) { // find ckov angle for this photon candidate fPhotCkov[fPhotCnt] = thetaCer; // actual theta Cerenkov (in TRS)