From 4cbc1eb8ea107378312784bbe51f54cffb95ba7d Mon Sep 17 00:00:00 2001 From: Ole Schmidt Date: Thu, 11 Jun 2020 11:24:53 +0200 Subject: [PATCH 1/6] Slow reset of TRD tracker not required anymore --- GPU/GPUTracking/TRDTracking/GPUTRDTracker.cxx | 20 +------------------ GPU/GPUTracking/TRDTracking/GPUTRDTracker.h | 2 +- 2 files changed, 2 insertions(+), 20 deletions(-) diff --git a/GPU/GPUTracking/TRDTracking/GPUTRDTracker.cxx b/GPU/GPUTracking/TRDTracking/GPUTRDTracker.cxx index 9b06a3267d662..d7979cea1f8f2 100644 --- a/GPU/GPUTracking/TRDTracking/GPUTRDTracker.cxx +++ b/GPU/GPUTracking/TRDTracking/GPUTRDTracker.cxx @@ -191,31 +191,13 @@ void GPUTRDTracker_t::InitializeProcessor() } template -void GPUTRDTracker_t::Reset(bool fast) +void GPUTRDTracker_t::Reset() { //-------------------------------------------------------------------- // Reset tracker //-------------------------------------------------------------------- mNTracklets = 0; mNTracks = 0; - if (fast) { - return; - } - for (int i = 0; i < mNMaxSpacePoints; ++i) { - mTracklets[i] = 0x0; - mSpacePoints[i].mR = 0.f; - mSpacePoints[i].mX[0] = 0.f; - mSpacePoints[i].mX[1] = 0.f; - mSpacePoints[i].mCov[0] = 0.f; - mSpacePoints[i].mCov[1] = 0.f; - mSpacePoints[i].mCov[2] = 0.f; - mSpacePoints[i].mDy = 0.f; - mSpacePoints[i].mId = 0; - mSpacePoints[i].mLabel[0] = -1; - mSpacePoints[i].mLabel[1] = -1; - mSpacePoints[i].mLabel[2] = -1; - mSpacePoints[i].mVolumeId = 0; - } } template diff --git a/GPU/GPUTracking/TRDTracking/GPUTRDTracker.h b/GPU/GPUTracking/TRDTracking/GPUTRDTracker.h index 51134e321fbbb..9fbbffd870026 100644 --- a/GPU/GPUTracking/TRDTracking/GPUTRDTracker.h +++ b/GPU/GPUTracking/TRDTracking/GPUTRDTracker.h @@ -107,7 +107,7 @@ class GPUTRDTracker_t : public GPUProcessor short MemoryTracks() const { return mMemoryTracks; } GPUhd() void OverrideGPUGeometry(TRD_GEOMETRY_CONST GPUTRDGeometry* geo) { mGeo = geo; } - void Reset(bool fast = false); + void Reset(); GPUd() int LoadTracklet(const GPUTRDTrackletWord& tracklet, const int* labels = nullptr); //template GPUd() int LoadTrack(const TRDTRK& trk, const int label = -1, const int* nTrkltsOffline = nullptr, const int labelOffline = -1) From 9f30127761a820a20231ea69f4f074d5a5c2048e Mon Sep 17 00:00:00 2001 From: Ole Schmidt Date: Wed, 17 Jun 2020 09:05:53 +0200 Subject: [PATCH 2/6] Add TF support to TRD tracker --- .../TRDTracking/GPUTRDInterfaces.h | 10 +- GPU/GPUTracking/TRDTracking/GPUTRDTracker.cxx | 99 ++++++++++++------- GPU/GPUTracking/TRDTracking/GPUTRDTracker.h | 18 +++- .../TRDTracking/macros/run_trd_tracker.C | 55 ++++++++--- 4 files changed, 130 insertions(+), 52 deletions(-) diff --git a/GPU/GPUTracking/TRDTracking/GPUTRDInterfaces.h b/GPU/GPUTracking/TRDTracking/GPUTRDInterfaces.h index e7800219b3e87..42d762b8480bd 100644 --- a/GPU/GPUTracking/TRDTracking/GPUTRDInterfaces.h +++ b/GPU/GPUTracking/TRDTracking/GPUTRDInterfaces.h @@ -80,6 +80,7 @@ class trackInterface : public AliExternalTrackParam const My_Float* getPar() const { return GetParameter(); } const My_Float* getCov() const { return GetCovariance(); } + float getTime() const { return -1.f; } bool CheckNumericalQuality() const { return true; } // parameter manipulation @@ -179,11 +180,16 @@ class trackInterface : public o2::dataformats::Tra } } - const float* getPar() { return getParams(); } + const float* getPar() const { return getParams(); } + float getTime() const { return mTime; } + void setTime(float t) { mTime = t; } bool CheckNumericalQuality() const { return true; } typedef o2::dataformats::TrackTPCITS baseClass; + + private: + float mTime; }; template <> @@ -294,6 +300,7 @@ class trackInterface : public GPUTPCGMTrackParam GPUd() const float* getPar() const { return GetPar(); } GPUd() const float* getCov() const { return GetCov(); } + GPUd() float getTime() const { return -1.f; } GPUd() void setAlpha(float alpha) { mAlpha = alpha; } GPUd() void set(float x, float alpha, const float param[5], const float cov[15]) @@ -338,7 +345,6 @@ class propagatorInterface : public GPUTPCGMPropagator { //bool ok = PropagateToXAlpha(x, GetAlpha(), true) == 0 ? true : false; int retVal = PropagateToXAlpha(x, GetAlpha(), true); - printf("Return value of PropagateToXAlpha: %i\n", retVal); bool ok = (retVal == 0) ? true : false; ok = mTrack->CheckNumericalQuality(); return ok; diff --git a/GPU/GPUTracking/TRDTracking/GPUTRDTracker.cxx b/GPU/GPUTracking/TRDTracking/GPUTRDTracker.cxx index d7979cea1f8f2..5f3c7814bf63b 100644 --- a/GPU/GPUTracking/TRDTracking/GPUTRDTracker.cxx +++ b/GPU/GPUTracking/TRDTracking/GPUTRDTracker.cxx @@ -75,7 +75,7 @@ void* GPUTRDTracker_t::SetPointersBase(void* base) //-------------------------------------------------------------------- mMaxThreads = mRec->GetMaxThreads(); computePointerWithAlignment(base, mR, kNChambers); - computePointerWithAlignment(base, mTrackletIndexArray, kNChambers + 1); + computePointerWithAlignment(base, mTrackletIndexArray, (kNChambers + 1) * mNMaxCollisions); computePointerWithAlignment(base, mHypothesis, mNCandidates * mMaxThreads); computePointerWithAlignment(base, mCandidates, mNCandidates * 2 * mMaxThreads); return base; @@ -88,9 +88,9 @@ void* GPUTRDTracker_t::SetPointersTracklets(void* base) // Allocate memory for tracklets and space points // (size might change for different events) //-------------------------------------------------------------------- - computePointerWithAlignment(base, mTracklets, mNMaxSpacePoints); - computePointerWithAlignment(base, mSpacePoints, mNMaxSpacePoints); - computePointerWithAlignment(base, mTrackletLabels, 3 * mNMaxSpacePoints); + computePointerWithAlignment(base, mTracklets, mNMaxSpacePoints * mNMaxCollisions); + computePointerWithAlignment(base, mSpacePoints, mNMaxSpacePoints * mNMaxCollisions); + computePointerWithAlignment(base, mTrackletLabels, 3 * mNMaxSpacePoints * mNMaxCollisions); return base; } @@ -105,7 +105,7 @@ void* GPUTRDTracker_t::SetPointersTracks(void* base) } template -GPUTRDTracker_t::GPUTRDTracker_t() : mR(nullptr), mIsInitialized(false), mMemoryPermanent(-1), mMemoryTracklets(-1), mMemoryTracks(-1), mNMaxTracks(0), mNMaxSpacePoints(0), mTracks(nullptr), mNCandidates(1), mNTracks(0), mNEvents(0), mTracklets(nullptr), mMaxThreads(100), mNTracklets(0), mTrackletIndexArray(nullptr), mHypothesis(nullptr), mCandidates(nullptr), mSpacePoints(nullptr), mTrackletLabels(nullptr), mGeo(nullptr), mRPhiA2(0), mRPhiB(0), mRPhiC2(0), mDyA2(0), mDyB(0), mDyC2(0), mAngleToDyA(0), mAngleToDyB(0), mAngleToDyC(0), mDebugOutput(false), mRadialOffset(-0.1), mMinPt(2.f), mMaxEta(0.84f), mExtraRoadY(2.f), mRoadZ(18.f), mMaxChi2(15.0f), mMaxMissingLy(6), mChi2Penalty(12.0f), mZCorrCoefNRC(1.4f), mMCEvent(nullptr), mDebug(new GPUTRDTrackerDebug()) +GPUTRDTracker_t::GPUTRDTracker_t() : mR(nullptr), mIsInitialized(false), mProcessPerTimeFrame(false), mMemoryPermanent(-1), mMemoryTracklets(-1), mMemoryTracks(-1), mNMaxCollisions(1), mNMaxTracks(0), mNMaxSpacePoints(0), mTracks(nullptr), mNCandidates(1), mNCollisions(1), mNTracks(0), mNEvents(0), mTriggerRecordIndices(nullptr), mTriggerRecordTimes(nullptr), mTracklets(nullptr), mMaxThreads(100), mNTracklets(0), mTrackletIndexArray(nullptr), mHypothesis(nullptr), mCandidates(nullptr), mSpacePoints(nullptr), mTrackletLabels(nullptr), mGeo(nullptr), mRPhiA2(0), mRPhiB(0), mRPhiC2(0), mDyA2(0), mDyB(0), mDyC2(0), mAngleToDyA(0), mAngleToDyB(0), mAngleToDyC(0), mDebugOutput(false), mTimeWindow(.1f), mRadialOffset(-0.1), mMinPt(2.f), mMaxEta(0.84f), mExtraRoadY(2.f), mRoadZ(18.f), mMaxChi2(15.0f), mMaxMissingLy(6), mChi2Penalty(12.0f), mZCorrCoefNRC(1.4f), mMCEvent(nullptr), mDebug(new GPUTRDTrackerDebug()) { //-------------------------------------------------------------------- // Default constructor @@ -208,28 +208,38 @@ void GPUTRDTracker_t::DoTracking(GPUChainTracking* chainTracking) //-------------------------------------------------------------------- // sort tracklets and fill index array - CAAlgo::sort(mTracklets, mTracklets + mNTracklets); // tracklets are sorted by HCId - int* trkltIndexArray = &mTrackletIndexArray[1]; - trkltIndexArray[-1] = 0; - int currDet = 0; - int nextDet = 0; - int trkltCounter = 0; - for (int iTrklt = 0; iTrklt < mNTracklets; ++iTrklt) { - if (mTracklets[iTrklt].GetDetector() > currDet) { - nextDet = mTracklets[iTrklt].GetDetector(); - for (int iDet = currDet; iDet < nextDet; ++iDet) { - trkltIndexArray[iDet] = trkltCounter; + for (int iColl = 0; iColl < mNCollisions; ++iColl) { + int nTrklts = 0; + if (mProcessPerTimeFrame) { + // FIXME maybe two nested if statements are not so good in terms of performance? + nTrklts = (iColl < mNCollisions - 1) ? mTriggerRecordIndices[iColl + 1] - mTriggerRecordIndices[iColl] : mNTracklets - mTriggerRecordIndices[iColl]; + } else { + nTrklts = mNTracklets; + } + GPUTRDTrackletWord* tracklets = (mProcessPerTimeFrame) ? &(mTracklets[mTriggerRecordIndices[iColl]]) : mTracklets; + CAAlgo::sort(tracklets, tracklets + nTrklts); // tracklets are sorted by HCId + int* trkltIndexArray = &mTrackletIndexArray[iColl * (kNChambers + 1) + 1]; + trkltIndexArray[-1] = 0; + int currDet = 0; + int nextDet = 0; + int trkltCounter = 0; + for (int iTrklt = 0; iTrklt < nTrklts; ++iTrklt) { + if (tracklets[iTrklt].GetDetector() > currDet) { + nextDet = tracklets[iTrklt].GetDetector(); + for (int iDet = currDet; iDet < nextDet; ++iDet) { + trkltIndexArray[iDet] = trkltCounter; + } + currDet = nextDet; } - currDet = nextDet; + ++trkltCounter; + } + for (int iDet = currDet; iDet <= kNChambers; ++iDet) { + trkltIndexArray[iDet] = trkltCounter; } - ++trkltCounter; - } - for (int iDet = currDet; iDet <= kNChambers; ++iDet) { - trkltIndexArray[iDet] = trkltCounter; - } - if (!CalculateSpacePoints()) { - Error("DoTracking", "Space points for at least one chamber could not be calculated"); + if (!CalculateSpacePoints(iColl)) { + GPUError("Space points for at least one chamber could not be calculated"); + } } auto timeStart = std::chrono::high_resolution_clock::now(); @@ -241,7 +251,7 @@ void GPUTRDTracker_t::DoTracking(GPUChainTracking* chainTracking) #pragma omp parallel for for (int iTrk = 0; iTrk < mNTracks; ++iTrk) { if (omp_get_num_threads() > mMaxThreads) { - Error("DoTracking", "number of parallel threads too high, aborting tracking"); + GPUError("Number of parallel threads too high, aborting tracking"); // break statement not possible in OpenMP for loop iTrk = mNTracks; continue; @@ -422,7 +432,7 @@ GPUd() int GPUTRDTracker_t::LoadTracklet(const GPUTRDTrackletWord& //-------------------------------------------------------------------- // Add single tracklet to tracker //-------------------------------------------------------------------- - if (mNTracklets >= mNMaxSpacePoints) { + if (mNTracklets >= mNMaxSpacePoints * mNMaxCollisions) { Error("LoadTracklet", "Running out of memory for tracklets, skipping tracklet(s). This should actually never happen."); return 1; } @@ -448,22 +458,43 @@ GPUd() void GPUTRDTracker_t::DumpTracks() } } +template +GPUd() int GPUTRDTracker_t::GetCollisionID(float trkTime) const +{ + for (int iColl = 0; iColl < mNCollisions; ++iColl) { + if (CAMath::Abs(trkTime - mTriggerRecordTimes[iColl]) < mTimeWindow) { + GPUInfo("TRD info found from interaction %i at %f for track with time %f", iColl, mTriggerRecordTimes[iColl], trkTime); + return iColl; + } + } + return -1; +} + template GPUd() void GPUTRDTracker_t::DoTrackingThread(int iTrk, int threadId) { //-------------------------------------------------------------------- // perform the tracking for one track (must be threadsafe) //-------------------------------------------------------------------- + int collisionId = 0; + if (mProcessPerTimeFrame) { + collisionId = GetCollisionID(mTracks[iTrk].getTime()); + if (collisionId < 0) { + GPUInfo("Did not find TRD data for track with t=%f", mTracks[iTrk].getTime()); + // no TRD data available for the bunch crossing this track originates from + return; + } + } PROP prop(&Param().polynomialField); auto trkCopy = mTracks[iTrk]; prop.setTrack(&trkCopy); prop.setFitInProjections(true); - FollowProlongation(&prop, &trkCopy, threadId); + FollowProlongation(&prop, &trkCopy, threadId, collisionId); mTracks[iTrk] = trkCopy; // copy back the resulting track } template -GPUd() bool GPUTRDTracker_t::CalculateSpacePoints() +GPUd() bool GPUTRDTracker_t::CalculateSpacePoints(int iCollision) { //-------------------------------------------------------------------- // Calculates TRD space points in sector tracking coordinates @@ -471,10 +502,11 @@ GPUd() bool GPUTRDTracker_t::CalculateSpacePoints() //-------------------------------------------------------------------- bool result = true; + int idxOffset = iCollision * (kNChambers + 1); for (int iDet = 0; iDet < kNChambers; ++iDet) { - int nTracklets = mTrackletIndexArray[iDet + 1] - mTrackletIndexArray[iDet]; + int nTracklets = mTrackletIndexArray[idxOffset + iDet + 1] - mTrackletIndexArray[idxOffset + iDet]; if (nTracklets == 0) { continue; } @@ -491,7 +523,7 @@ GPUd() bool GPUTRDTracker_t::CalculateSpacePoints() float c2 = 1.f / (1.f + t2); // cos^2 (tilt) float sy2 = 0.1f * 0.1f; // sigma_rphi^2, currently assume sigma_rphi = 1 mm - for (int trkltIdx = mTrackletIndexArray[iDet]; trkltIdx < mTrackletIndexArray[iDet + 1]; ++trkltIdx) { + for (int trkltIdx = mTrackletIndexArray[idxOffset + iDet]; trkltIdx < mTrackletIndexArray[idxOffset + iDet + 1]; ++trkltIdx) { int trkltZbin = mTracklets[trkltIdx].GetZbin(); float sz2 = pp->GetRowSize(trkltZbin) * pp->GetRowSize(trkltZbin) / 12.f; // sigma_z = l_pad/sqrt(12) TODO try a larger z error My_Float xTrkltDet[3] = {0.f}; // trklt position in chamber coordinates @@ -523,7 +555,7 @@ GPUd() bool GPUTRDTracker_t::CalculateSpacePoints() } template -GPUd() bool GPUTRDTracker_t::FollowProlongation(PROP* prop, TRDTRK* t, int threadId) +GPUd() bool GPUTRDTracker_t::FollowProlongation(PROP* prop, TRDTRK* t, int threadId, int collisionId) { //-------------------------------------------------------------------- // Propagate TPC track layerwise through TRD and pick up closest @@ -556,6 +588,7 @@ GPUd() bool GPUTRDTracker_t::FollowProlongation(PROP* prop, TRDTRK int candidateIdxOffset = threadId * 2 * mNCandidates; int hypothesisIdxOffset = threadId * mNCandidates; + int trkltIdxOffset = collisionId * (kNChambers + 1); auto trkWork = t; if (mNCandidates > 1) { @@ -671,7 +704,7 @@ GPUd() bool GPUTRDTracker_t::FollowProlongation(PROP* prop, TRDTRK } } // first propagate track to x of tracklet - for (int trkltIdx = mTrackletIndexArray[currDet]; trkltIdx < mTrackletIndexArray[currDet + 1]; ++trkltIdx) { + for (int trkltIdx = mTrackletIndexArray[trkltIdxOffset + currDet]; trkltIdx < mTrackletIndexArray[trkltIdxOffset + currDet + 1]; ++trkltIdx) { if (CAMath::Abs(trkWork->getY() - mSpacePoints[trkltIdx].mX[0]) > roadY || CAMath::Abs(trkWork->getZ() - mSpacePoints[trkltIdx].mX[1]) > roadZ) { // skip tracklets which are too far away // although the radii of space points and tracks may differ by ~ few mm the roads are large enough to allow no efficiency loss by this cut @@ -712,7 +745,7 @@ GPUd() bool GPUTRDTracker_t::FollowProlongation(PROP* prop, TRDTRK } // end candidate loop #ifdef ENABLE_GPUMC - // in case matching tracklet exists in this layer -> store position information for debugging + // in case matching tracklet exists in this layer -> store position information for debugging FIXME: does not yet work for time frames in o2, but here we anyway do not yet have MC labels... if (matchAvailableAll[iLayer].size() > 0 && mDebugOutput) { mDebug->SetNmatchAvail(matchAvailableAll[iLayer].size(), iLayer); int realTrkltId = matchAvailableAll[iLayer].at(0); diff --git a/GPU/GPUTracking/TRDTracking/GPUTRDTracker.h b/GPU/GPUTracking/TRDTracking/GPUTRDTracker.h index 9fbbffd870026..0b79fa6bad325 100644 --- a/GPU/GPUTracking/TRDTracking/GPUTRDTracker.h +++ b/GPU/GPUTracking/TRDTracking/GPUTRDTracker.h @@ -145,9 +145,10 @@ class GPUTRDTracker_t : public GPUProcessor mNTracks++; return (0); } + GPUd() int GetCollisionID(float trkTime) const; GPUd() void DoTrackingThread(int iTrk, int threadId = 0); - GPUd() bool CalculateSpacePoints(); - GPUd() bool FollowProlongation(PROP* prop, TRDTRK* t, int threadId); + GPUd() bool CalculateSpacePoints(int iCollision = 0); + GPUd() bool FollowProlongation(PROP* prop, TRDTRK* t, int threadId, int collisionId); GPUd() float GetPredictedChi2(const My_Float* pTRD, const My_Float* covTRD, const My_Float* pTrk, const My_Float* covTrk) const; GPUd() int GetDetectorNumber(const float zPos, const float alpha, const int layer) const; GPUd() bool AdjustSector(PROP* prop, TRDTRK* t, const int layer) const; @@ -166,7 +167,14 @@ class GPUTRDTracker_t : public GPUProcessor GPUd() void Quicksort(const int left, const int right, const int size); GPUd() void InsertHypothesis(Hypothesis hypo, int& nCurrHypothesis, int idxOffset); + // input from TRD trigger record + GPUd() void SetNMaxCollisions(int nColl) { mNMaxCollisions = nColl; } // can this be fixed to a sufficiently large value? + GPUd() void SetNCollisions(int nColl) { mNCollisions = nColl; } + GPUd() void SetTriggerRecordIndices(int* indices) { mTriggerRecordIndices = indices; } + GPUd() void SetTriggerRecordTimes(float* times) { mTriggerRecordTimes = times; } + // settings + GPUd() void SetProcessPerTimeFrame() { mProcessPerTimeFrame = true; } GPUd() void SetMCEvent(AliMCEvent* mc) { mMCEvent = mc; } GPUd() void EnableDebugOutput() { mDebugOutput = true; } GPUd() void SetPtThreshold(float minPt) { mMinPt = minPt; } @@ -199,15 +207,20 @@ class GPUTRDTracker_t : public GPUProcessor protected: float* mR; // radial position of each TRD chamber, alignment taken into account, radial spread within chambers < 7mm bool mIsInitialized; // flag is set upon initialization + bool mProcessPerTimeFrame; // if true, tracking is done per time frame instead of on a single events basis //FIXME is this needed?? short mMemoryPermanent; // size of permanent memory for the tracker short mMemoryTracklets; // size of memory for TRD tracklets short mMemoryTracks; // size of memory for tracks (used for i/o) + int mNMaxCollisions; // max number of collisions to process (per time frame) int mNMaxTracks; // max number of tracks the tracker can handle (per event) int mNMaxSpacePoints; // max number of space points hold by the tracker (per event) TRDTRK* mTracks; // array of trd-updated tracks int mNCandidates; // max. track hypothesis per layer + int mNCollisions; // number of collisions with TRD tracklet data int mNTracks; // number of TPC tracks to be matched int mNEvents; // number of processed events + int* mTriggerRecordIndices; // index of first tracklet for each collision + float* mTriggerRecordTimes; // time in us for each collision GPUTRDTrackletWord* mTracklets; // array of all tracklets, later sorted by HCId int mMaxThreads; // maximum number of supported threads int mNTracklets; // total number of tracklets in event @@ -229,6 +242,7 @@ class GPUTRDTracker_t : public GPUProcessor float mAngleToDyC; // parameterization for conversion track angle -> tracklet deflection /// ---- end error parametrization ---- bool mDebugOutput; // store debug output + float mTimeWindow; // max. deviation of the ITS-TPC track time w.r.t. TRD trigger record time stamp (in us, default is 100 ns) float mRadialOffset; // due to mis-calibration of t0 float mMinPt; // min pt of TPC tracks for tracking float mMaxEta; // TPC tracks with higher eta are ignored diff --git a/GPU/GPUTracking/TRDTracking/macros/run_trd_tracker.C b/GPU/GPUTracking/TRDTracking/macros/run_trd_tracker.C index c7fc292457c82..198c13fcb243c 100644 --- a/GPU/GPUTracking/TRDTracking/macros/run_trd_tracker.C +++ b/GPU/GPUTracking/TRDTracking/macros/run_trd_tracker.C @@ -15,24 +15,29 @@ #include "GPUTRDGeometry.h" // O2 header +#include "DetectorsCommonDataFormats/NameConf.h" +#include "CommonConstants/LHCConstants.h" #include "DetectorsBase/GeometryManager.h" #include "DetectorsBase/Propagator.h" #include "ReconstructionDataFormats/TrackTPCITS.h" #include "TRDBase/Tracklet.h" +#include "DataFormatsTRD/TriggerRecord.h" #endif using namespace GPUCA_NAMESPACE::gpu; void run_trd_tracker(std::string path = "./", - std::string inputGRP = "o2sim_grp.root", std::string inputTracks = "o2match_itstpc.root", std::string inputTracklets = "trdtracklets.root") { + //-------- debug time information from tracks and tracklets + std::vector trdTriggerTimes; + std::vector trdTriggerIndices; //-------- init geometry and field --------// - o2::base::GeometryManager::loadGeometry(path); - o2::base::Propagator::initFieldFromGRP(path + inputGRP); + o2::base::GeometryManager::loadGeometry(); + o2::base::Propagator::initFieldFromGRP(o2::base::NameConf::getGRPFileName()); auto geo = o2::trd::TRDGeometry::instance(); geo->createPadPlaneArray(); @@ -43,9 +48,11 @@ void run_trd_tracker(std::string path = "./", GPUSettingsEvent cfgEvent; // defaults should be ok GPUSettingsRec cfgRec; // don't care for now, NWaysOuter is set in here for instance GPUSettingsProcessing cfgDeviceProcessing; // also keep defaults here, or adjust debug level - cfgDeviceProcessing.debugLevel = 10; + cfgDeviceProcessing.debugLevel = 5; GPURecoStepConfiguration cfgRecoStep; cfgRecoStep.steps = GPUDataTypes::RecoStep::NoRecoStep; + cfgRecoStep.inputs.clear(); + cfgRecoStep.outputs.clear(); auto rec = GPUReconstruction::CreateInstance("CPU", true); rec->SetSettings(&cfgEvent, &cfgRec, &cfgDeviceProcessing, &cfgRecoStep); @@ -53,12 +60,14 @@ void run_trd_tracker(std::string path = "./", auto tracker = new GPUTRDTracker(); tracker->SetNCandidates(1); // must be set before initialization + tracker->SetProcessPerTimeFrame(); + tracker->SetNMaxCollisions(100); rec->RegisterGPUProcessor(tracker, false); chainTracking->SetTRDGeometry(&geoFlat); - tracker->SetTrackingChain(chainTracking); - rec->Init(); - rec->AllocateRegisteredMemory(nullptr); + if (rec->Init()) { + printf("ERROR: GPUReconstruction not initialized\n"); + } // configure the tracker //tracker->EnableDebugOutput(); @@ -85,23 +94,34 @@ void run_trd_tracker(std::string path = "./", TChain trdTracklets("o2sim"); trdTracklets.AddFile((path + inputTracklets).c_str()); + std::vector* triggerRecordsInArrayPtr = nullptr; + trdTracklets.SetBranchAddress("TrackTrg", &triggerRecordsInArrayPtr); std::vector* trackletsInArrayPtr = nullptr; trdTracklets.SetBranchAddress("Tracklet", &trackletsInArrayPtr); trdTracklets.GetEntry(0); + int nCollisions = triggerRecordsInArrayPtr->size(); int nTracklets = trackletsInArrayPtr->size(); - printf("There are %i tracklets in total\n", nTracklets); - - auto pp = geo->getPadPlane(0, 0); - printf("Tilt=%f\n", pp->getTiltingAngle()); + printf("There are %i tracklets in total from %i trigger records\n", nTracklets, nCollisions); + + for (int iEv = 0; iEv < nCollisions; ++iEv) { + o2::trd::TriggerRecord& trg = triggerRecordsInArrayPtr->at(iEv); + int nTrackletsCurrent = trg.getNumberOfObjects(); + int iFirstTracklet = trg.getFirstEntry(); + int64_t evTime = trg.getBCData().toLong() * o2::constants::lhc::LHCBunchSpacingNS; // event time in ns + trdTriggerTimes.push_back(evTime / 1000.); + trdTriggerIndices.push_back(iFirstTracklet); + printf("Event %i: Occured at %li us after SOR, contains %i tracklets, index of first tracklet is %i\n", iEv, evTime / 1000, nTrackletsCurrent, iFirstTracklet); + } - tracker->Reset(true); + tracker->Reset(); chainTracking->mIOPtrs.nMergedTracks = nTracks; chainTracking->mIOPtrs.nTRDTracklets = nTracklets; chainTracking->AllocateIOMemory(); rec->PrepareEvent(); - rec->AllocateRegisteredMemory(tracker->MemoryTracks()); + rec->SetupGPUProcessor(tracker, true); + printf("Start loading input into TRD tracker\n"); // load everything into the tracker for (int iTrk = 0; iTrk < nTracks; ++iTrk) { auto trk = tracksInArrayPtr->at(iTrk); @@ -114,8 +134,11 @@ void run_trd_tracker(std::string path = "./", for (int i = 0; i < 15; ++i) { trkLoad.setCov(trk.getCov()[i], i); } + trkLoad.setTime(trk.getTimeMUS().getTimeStamp()); tracker->LoadTrack(trkLoad); + printf("Loaded track %i with time %f\n", iTrk, trkLoad.getTime()); } + for (int iTrklt = 0; iTrklt < nTracklets; ++iTrklt) { auto trklt = trackletsInArrayPtr->at(iTrklt); unsigned int trkltWord = trklt.getTrackletWord(); @@ -127,9 +150,11 @@ void run_trd_tracker(std::string path = "./", printf("Could not load tracklet %i\n", iTrklt); } } + tracker->SetTriggerRecordTimes(&(trdTriggerTimes[0])); + tracker->SetTriggerRecordIndices(&(trdTriggerIndices[0])); + tracker->SetNCollisions(nCollisions); tracker->DumpTracks(); - tracker->DoTracking(); + tracker->DoTracking(chainTracking); tracker->DumpTracks(); - printf("Done\n"); } From 6a779cf561d505ee995fdd9b299312d876dadb4f Mon Sep 17 00:00:00 2001 From: Ole Schmidt Date: Mon, 7 Sep 2020 14:27:24 +0200 Subject: [PATCH 3/6] Remove erroneous chi2 calculation causing FPEs --- GPU/GPUTracking/TRDTracking/GPUTRDTracker.cxx | 38 +------------------ GPU/GPUTracking/TRDTracking/GPUTRDTracker.h | 4 +- .../TRDTracking/GPUTRDTrackerDebug.h | 11 ------ 3 files changed, 2 insertions(+), 51 deletions(-) diff --git a/GPU/GPUTracking/TRDTracking/GPUTRDTracker.cxx b/GPU/GPUTracking/TRDTracking/GPUTRDTracker.cxx index 5f3c7814bf63b..d9a54dfa2391c 100644 --- a/GPU/GPUTracking/TRDTracking/GPUTRDTracker.cxx +++ b/GPU/GPUTracking/TRDTracking/GPUTRDTracker.cxx @@ -731,7 +731,7 @@ GPUd() bool GPUTRDTracker_t::FollowProlongation(PROP* prop, TRDTRK float chi2 = prop->getPredictedChi2(trkltPosTmpYZ, trkltCovTmp); // GPUInfo("layer %i: chi2 = %f", iLayer, chi2); if (chi2 < mMaxChi2 && CAMath::Abs(GetAngularPull(mSpacePoints[trkltIdx].mDy, trkWork->getSnp())) < 4) { - Hypothesis hypo(trkWork->GetNlayers(), iCandidate, trkltIdx, trkWork->GetChi2() + chi2, GetPredictedChi2(trkltPosTmpYZ, trkltCovTmp, trkWork->getPar(), trkWork->getCov())); + Hypothesis hypo(trkWork->GetNlayers(), iCandidate, trkltIdx, trkWork->GetChi2() + chi2); InsertHypothesis(hypo, nCurrHypothesis, hypothesisIdxOffset); } // end tracklet chi2 < mMaxChi2 } // end tracklet in window @@ -769,7 +769,6 @@ GPUd() bool GPUTRDTracker_t::FollowProlongation(PROP* prop, TRDTRK My_Float covReal[3] = {0.}; RecalcTrkltCov(tilt, trkWork->getSnp(), pad->GetRowSize(mTracklets[realTrkltId].GetZbin()), covReal); mDebug->SetChi2Real(prop->getPredictedChi2(yzPosReal, covReal), iLayer); - mDebug->SetChi2YZPhiReal(GetPredictedChi2(yzPosReal, covReal, trkWork->getPar(), trkWork->getCov()), iLayer); mDebug->SetRawTrackletPositionReal(mSpacePoints[realTrkltId].mR, mSpacePoints[realTrkltId].mX, iLayer); mDebug->SetCorrectedTrackletPositionReal(yzPosReal, iLayer); mDebug->SetTrackletPropertiesReal(mTracklets[realTrkltId].GetDetector(), iLayer); @@ -778,7 +777,6 @@ GPUd() bool GPUTRDTracker_t::FollowProlongation(PROP* prop, TRDTRK #endif // mDebug->SetChi2Update(mHypothesis[0 + hypothesisIdxOffset].mChi2 - t->GetChi2(), iLayer); // only meaningful for ONE candidate!!! - mDebug->SetChi2YZPhiUpdate(mHypothesis[0 + hypothesisIdxOffset].mChi2YZPhi, iLayer); // only meaningful for ONE candidate!!! mDebug->SetRoad(roadY, roadZ, iLayer); // only meaningful for ONE candidate bool wasTrackStored = false; // -------------------------------------------------------------------------------- @@ -989,40 +987,6 @@ GPUd() bool GPUTRDTracker_t::FollowProlongation(PROP* prop, TRDTRK return true; } -template -GPUd() float GPUTRDTracker_t::GetPredictedChi2(const My_Float* pTRD, const My_Float* covTRD, const My_Float* pTrk, const My_Float* covTrk) const -{ - // FIXME: This function cannot yet be used in production (maybe it is not necessary at all?) - // conversion from tracklet deflection to sin(phi) needs to be parametrized - // predict chi2 for update of track with pTrk and covTrk with TRD space point with pTRD and covTRD - // taking into account y, z and azimuthal angle of the track and tracklet - float deltaY = pTrk[0] - pTRD[0]; - float deltaZ = pTrk[1] - pTRD[1]; - float deltaS = pTrk[2] - pTRD[2]; // FIXME: pTRD[2] is not defined, needs conversion dy -> sin(phi) - - // add errors for track and space point, assume no correlation in y-sin(phi) and z-sin(phi) for space point - float sigmaZ2 = covTrk[2] + covTRD[2]; - float sigmaS2 = covTrk[5] + GetAngularResolution(pTrk[2]); // FIXME: convert angular resolution from dy to sin(phi) for the TRD space point - float sigmaY2 = covTrk[0] + covTRD[0]; - float sigmaZS = covTrk[4]; - float sigmaYS = covTrk[3]; - float sigmaYZ = covTrk[1] + covTRD[1]; - // inverse of the covariance matrix - float c11 = sigmaZ2 * sigmaS2 - sigmaZS * sigmaZS; - float c21 = sigmaZS * sigmaYS - sigmaYZ * sigmaS2; - float c22 = sigmaY2 * sigmaS2 - sigmaYS * sigmaYS; - float c31 = sigmaYZ * sigmaZS - sigmaZ2 * sigmaYS; - float c32 = sigmaYZ * sigmaYS - sigmaY2 * sigmaZS; - float c33 = sigmaY2 * sigmaZ2 - sigmaYZ * sigmaYZ; - // determinant - float det = sigmaY2 * sigmaZ2 * sigmaS2 + 2 * sigmaYZ * sigmaZS * sigmaYS - sigmaYS * sigmaYS * sigmaZ2 - sigmaZS * sigmaZS * sigmaY2 - sigmaS2 * sigmaYZ * sigmaYZ; - if (CAMath::Abs(det) < 1.e-10f) { - printf("Determinant too small: %f\n", det); - det = 1.e-10f; - } - det = 1.f / det; - return (c11 * deltaY * deltaY + 2.f * c21 * deltaY * deltaZ + 2.f * c31 * deltaY * deltaS + c22 * deltaZ * deltaZ + 2.f * c32 * deltaZ * deltaS + c33 * deltaS * deltaS) * det; -} template GPUd() void GPUTRDTracker_t::InsertHypothesis(Hypothesis hypo, int& nCurrHypothesis, int idxOffset) diff --git a/GPU/GPUTracking/TRDTracking/GPUTRDTracker.h b/GPU/GPUTracking/TRDTracking/GPUTRDTracker.h index 0b79fa6bad325..62f2193f023c4 100644 --- a/GPU/GPUTracking/TRDTracking/GPUTRDTracker.h +++ b/GPU/GPUTracking/TRDTracking/GPUTRDTracker.h @@ -95,11 +95,10 @@ class GPUTRDTracker_t : public GPUProcessor int mCandidateId; // to which track candidate the hypothesis belongs int mTrackletId; // tracklet index to be used for update float mChi2; // predicted chi2 for given space point - float mChi2YZPhi; // not yet ready (see GetPredictedChi2 method in cxx file) GPUd() float GetReducedChi2() { return mLayers > 0 ? mChi2 / mLayers : mChi2; } GPUd() Hypothesis() : mLayers(0), mCandidateId(-1), mTrackletId(-1), mChi2(9999.f) {} - GPUd() Hypothesis(int layers, int candidateId, int trackletId, float chi2, float chi2YZPhi = -1.f) : mLayers(layers), mCandidateId(candidateId), mTrackletId(trackletId), mChi2(chi2), mChi2YZPhi(chi2YZPhi) {} + GPUd() Hypothesis(int layers, int candidateId, int trackletId, float chi2, float chi2YZPhi = -1.f) : mLayers(layers), mCandidateId(candidateId), mTrackletId(trackletId), mChi2(chi2) {} }; short MemoryPermanent() const { return mMemoryPermanent; } @@ -149,7 +148,6 @@ class GPUTRDTracker_t : public GPUProcessor GPUd() void DoTrackingThread(int iTrk, int threadId = 0); GPUd() bool CalculateSpacePoints(int iCollision = 0); GPUd() bool FollowProlongation(PROP* prop, TRDTRK* t, int threadId, int collisionId); - GPUd() float GetPredictedChi2(const My_Float* pTRD, const My_Float* covTRD, const My_Float* pTrk, const My_Float* covTrk) const; GPUd() int GetDetectorNumber(const float zPos, const float alpha, const int layer) const; GPUd() bool AdjustSector(PROP* prop, TRDTRK* t, const int layer) const; GPUd() int GetSector(float alpha) const; diff --git a/GPU/GPUTracking/TRDTracking/GPUTRDTrackerDebug.h b/GPU/GPUTracking/TRDTracking/GPUTRDTrackerDebug.h index 209dfdfc807f3..08b453d4f2345 100644 --- a/GPU/GPUTracking/TRDTracking/GPUTRDTrackerDebug.h +++ b/GPU/GPUTracking/TRDTracking/GPUTRDTrackerDebug.h @@ -95,9 +95,7 @@ class GPUTRDTrackerDebug fTrackYReal.ResizeTo(6); fTrackZReal.ResizeTo(6); fTrackSecReal.ResizeTo(6); - fChi2YZPhiUpdate.ResizeTo(6); fChi2Update.ResizeTo(6); - fChi2YZPhiReal.ResizeTo(6); fChi2Real.ResizeTo(6); fNmatchesAvail.ResizeTo(6); fFindable.ResizeTo(6); @@ -151,9 +149,7 @@ class GPUTRDTrackerDebug fTrackYReal.Zero(); fTrackZReal.Zero(); fTrackSecReal.Zero(); - fChi2YZPhiUpdate.Zero(); fChi2Update.Zero(); - fChi2YZPhiReal.Zero(); fChi2Real.Zero(); fNmatchesAvail.Zero(); fFindable.Zero(); @@ -290,8 +286,6 @@ class GPUTRDTrackerDebug // update information void SetChi2Update(float chi2, int ly) { fChi2Update(ly) = chi2; } void SetChi2Real(float chi2, int ly) { fChi2Real(ly) = chi2; } - void SetChi2YZPhiUpdate(float chi2, int ly) { fChi2YZPhiUpdate(ly) = chi2; } - void SetChi2YZPhiReal(float chi2, int ly) { fChi2YZPhiReal(ly) = chi2; } // other infos void SetRoad(float roadY, float roadZ, int ly) @@ -370,8 +364,6 @@ class GPUTRDTrackerDebug "trackletDetReal.=" << &fTrackletDetReal << // detector number for matching or related tracklet if available, otherwise -1 "chi2Update.=" << &fChi2Update << // chi2 for update "chi2Real.=" << &fChi2Real << // chi2 for first tracklet w/ matching MC label - "chi2YZPhiUpdate.=" << &fChi2YZPhiUpdate << // chi2 for update taking into account full tracklet information (y, z, dY aka sin(phi)) - "chi2YZPhiReal.=" << &fChi2YZPhiReal << // chi2 for first tracklet w/ matching MC label taking into account full tracklet information (y, z, dY aka sin(phi)) "chi2Total=" << fChi2 << // total chi2 for track "nLayers=" << fNlayers << // number of layers in which track was findable "nTracklets=" << fNtrklts << // number of attached tracklets @@ -456,8 +448,6 @@ class GPUTRDTrackerDebug TVectorF fTrackletDetReal; TVectorF fChi2Update; TVectorF fChi2Real; - TVectorF fChi2YZPhiUpdate; - TVectorF fChi2YZPhiReal; TVectorF fRoadY; TVectorF fRoadZ; TVectorF fFindable; @@ -512,7 +502,6 @@ class GPUTRDTrackerDebug GPUd() void SetChi2Update(float chi2, int ly) {} GPUd() void SetChi2Real(float chi2, int ly) {} GPUd() void SetChi2YZPhiUpdate(float chi2, int ly) {} - GPUd() void SetChi2YZPhiReal(float chi2, int ly) {} // other infos GPUd() void SetRoad(float roadY, float roadZ, int ly) {} From 526ac156ec746ab52836d28778bd623e5bacf922 Mon Sep 17 00:00:00 2001 From: Ole Schmidt Date: Tue, 15 Sep 2020 15:59:22 +0200 Subject: [PATCH 4/6] Adapt TRD tracker steering macro to new data type --- .../TRDTracking/macros/run_trd_tracker.C | 19 +++++++++++++++---- 1 file changed, 15 insertions(+), 4 deletions(-) diff --git a/GPU/GPUTracking/TRDTracking/macros/run_trd_tracker.C b/GPU/GPUTracking/TRDTracking/macros/run_trd_tracker.C index 198c13fcb243c..9d5212704ef4e 100644 --- a/GPU/GPUTracking/TRDTracking/macros/run_trd_tracker.C +++ b/GPU/GPUTracking/TRDTracking/macros/run_trd_tracker.C @@ -19,14 +19,25 @@ #include "CommonConstants/LHCConstants.h" #include "DetectorsBase/GeometryManager.h" #include "DetectorsBase/Propagator.h" +#include "TRDBase/TRDGeometry.h" #include "ReconstructionDataFormats/TrackTPCITS.h" -#include "TRDBase/Tracklet.h" +#include "DataFormatsTRD/Tracklet64.h" #include "DataFormatsTRD/TriggerRecord.h" #endif using namespace GPUCA_NAMESPACE::gpu; +unsigned int convertTrkltWordToRun2Format(uint64_t trkltWordRun3) +{ + // FIXME: this is currently a dummy function + // need proper functionality to convert the new tracklet data format to + // something compatible with the TRD tracker, but this macro is probably + // not the right place for this + unsigned int trkltWord = 0; + return trkltWord; +} + void run_trd_tracker(std::string path = "./", std::string inputTracks = "o2match_itstpc.root", std::string inputTracklets = "trdtracklets.root") @@ -96,7 +107,7 @@ void run_trd_tracker(std::string path = "./", std::vector* triggerRecordsInArrayPtr = nullptr; trdTracklets.SetBranchAddress("TrackTrg", &triggerRecordsInArrayPtr); - std::vector* trackletsInArrayPtr = nullptr; + std::vector* trackletsInArrayPtr = nullptr; trdTracklets.SetBranchAddress("Tracklet", &trackletsInArrayPtr); trdTracklets.GetEntry(0); int nCollisions = triggerRecordsInArrayPtr->size(); @@ -141,10 +152,10 @@ void run_trd_tracker(std::string path = "./", for (int iTrklt = 0; iTrklt < nTracklets; ++iTrklt) { auto trklt = trackletsInArrayPtr->at(iTrklt); - unsigned int trkltWord = trklt.getTrackletWord(); + unsigned int trkltWord = convertTrkltWordToRun2Format(trklt.getTrackletWord()); GPUTRDTrackletWord trkltLoad; trkltLoad.SetId(iTrklt); - trkltLoad.SetHCId(trklt.getHCId()); + trkltLoad.SetHCId(trklt.getHCID()); trkltLoad.SetTrackletWord(trkltWord); if (tracker->LoadTracklet(trkltLoad) > 0) { printf("Could not load tracklet %i\n", iTrklt); From dbc7467b31b3e8f468cf68100cffd9af681b686d Mon Sep 17 00:00:00 2001 From: Ole Schmidt Date: Wed, 16 Sep 2020 09:00:12 +0200 Subject: [PATCH 5/6] Add missing LinkDef for Tracklet64 class --- DataFormats/Detectors/TRD/src/DataFormatsTRDLinkDef.h | 1 + 1 file changed, 1 insertion(+) diff --git a/DataFormats/Detectors/TRD/src/DataFormatsTRDLinkDef.h b/DataFormats/Detectors/TRD/src/DataFormatsTRDLinkDef.h index f4e87501d4f9f..4fe89c45df5ff 100644 --- a/DataFormats/Detectors/TRD/src/DataFormatsTRDLinkDef.h +++ b/DataFormats/Detectors/TRD/src/DataFormatsTRDLinkDef.h @@ -21,6 +21,7 @@ #pragma link C++ struct o2::trd::TrackletMCMHeader + ; #pragma link C++ struct o2::trd::TrackletMCMData + ; #pragma link C++ class o2::trd::Tracklet64 + ; +#pragma link C++ class std::vector < o2::trd::Tracklet64> + ; #pragma link C++ class std::vector < o2::trd::TriggerRecord > +; #pragma link C++ class std::vector < o2::trd::LinkRecord > +; From a176447c299592ed0fe782a259f9f20e62c0932c Mon Sep 17 00:00:00 2001 From: Ole Schmidt Date: Fri, 25 Sep 2020 09:46:26 +0200 Subject: [PATCH 6/6] Disable material correction for O2 tracking --- GPU/GPUTracking/TRDTracking/GPUTRDInterfaces.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/GPU/GPUTracking/TRDTracking/GPUTRDInterfaces.h b/GPU/GPUTracking/TRDTracking/GPUTRDInterfaces.h index 42d762b8480bd..990b9aa318c51 100644 --- a/GPU/GPUTracking/TRDTracking/GPUTRDInterfaces.h +++ b/GPU/GPUTracking/TRDTracking/GPUTRDInterfaces.h @@ -200,7 +200,7 @@ class propagatorInterface propagatorInterface(const propagatorInterface&) = delete; propagatorInterface& operator=(const propagatorInterface&) = delete; - bool propagateToX(float x, float maxSnp, float maxStep) { return mProp->PropagateToXBxByBz(*mParam, x, 0.13957, maxSnp, maxStep); } + bool propagateToX(float x, float maxSnp, float maxStep) { return mProp->PropagateToXBxByBz(*mParam, x, 0.13957, maxSnp, maxStep, o2::base::Propagator::MatCorrType::USEMatCorrNONE); } int getPropagatedYZ(My_Float x, My_Float& projY, My_Float& projZ) { return static_cast(mParam->getYZAt(x, mProp->getNominalBz(), projY, projZ)); } void setTrack(trackInterface* trk) { mParam = trk; }