Skip to content

Commit

Permalink
Fixing HMP matcher (#12161)
Browse files Browse the repository at this point in the history 
* Fixing HMP matcher

* Minors

* Addign track charge sign to track momentum

* Added assignable extra time tolerance

* Updating MatchHMP class

* Minors

* fix braces

---------

Co-authored-by: Ruben Shahoyan <[email protected]>
  • Loading branch information
@gvolpe79 @shahor02
gvolpe79 and shahor02 authored Oct 30, 2023
1 parent 654fa61 commit 69b0e45
Show file tree
Hide file tree
Showing 6 changed files with 112 additions and 90 deletions.
31 changes: 21 additions & 10 deletions Detectors/GlobalTracking/include/GlobalTracking/MatchHMP.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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"
Expand Down Expand Up @@ -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,
Expand Down Expand Up @@ -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
Expand Down
157 changes: 81 additions & 76 deletions Detectors/GlobalTracking/src/MatchHMP.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -95,7 +97,6 @@ void MatchHMP::run(const o2::globaltracking::RecoContainer& inp)
}

if (mIsITSTPCused || mIsTPCTRDused || mIsITSTPCTRDused || mIsITSTPCTOFused || mIsTPCTOFused || mIsITSTPCTRDTOFused || mIsTPCTRDTOFused) {
// doFastMatching();
doMatching();
}

Expand Down Expand Up @@ -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<float, 3> 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)
Expand Down Expand Up @@ -305,6 +324,8 @@ bool MatchHMP::prepareHMPClusters()

mNumOfTriggers = 0;

mHMPTriggersWork.clear();

int nTriggersInCurrentChunk = mHMPTriggersArray.size();
LOG(debug) << "nTriggersInCurrentChunk = " << nTriggersInCurrentChunk;
mNumOfTriggers += nTriggersInCurrentChunk;
Expand Down Expand Up @@ -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;

Expand All @@ -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;

Expand All @@ -433,15 +452,21 @@ void MatchHMP::doMatching()
oneEventClusters.push_back(cluster);
double qthre = pParam->qCut();

if (cluster.q() < 150.) {
if (cluster.q() < 150. || cluster.size() > 10) {
continue;
}

isOkQcut = kTRUE;

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;
Expand All @@ -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<float, 2> trkPos{0, gz};
std::array<float, 3> trkCov{0.1 * 0.1, 0., 0.1 * 0.1};
Expand All @@ -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
Expand All @@ -537,41 +547,36 @@ 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) {
isMatched = kTRUE;
} // 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...

double nmean = pParam->meanIdxRad();

// 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)
Expand Down
2 changes: 1 addition & 1 deletion Detectors/GlobalTrackingWorkflow/include/GlobalTrackingWorkflow/HMPMatcherSpec.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
2 changes: 1 addition & 1 deletion Detectors/GlobalTrackingWorkflow/src/HMPMatcherSpec.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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<DataRequest>();
Expand Down
Loading

0 comments on commit 69b0e45

Please sign in to comment.