diff --git a/DataFormats/Detectors/TRD/include/DataFormatsTRD/Constants.h b/DataFormats/Detectors/TRD/include/DataFormatsTRD/Constants.h index 1a4cd01626406..79d6afad0f423 100644 --- a/DataFormats/Detectors/TRD/include/DataFormatsTRD/Constants.h +++ b/DataFormats/Detectors/TRD/include/DataFormatsTRD/Constants.h @@ -56,6 +56,7 @@ constexpr float GRANULARITYTRKLSLOPE = 1.f / PADGRANULARITYTRKLSLOPE; // granula // OS: Should this not be flexible for example in case of Kr calib? constexpr int TIMEBINS = 30; // the number of time bins +constexpr int NBINSANGLEDIFF = 26; // the number of bins for the track angle used for the vDrift and ExB calibration based on the tracking (last bin is for under-/overflow entries) // Trigger parameters constexpr double READOUT_TIME = 3000; // the time the readout takes, as 30 TB = 3 micro-s. diff --git a/Detectors/TRD/CMakeLists.txt b/Detectors/TRD/CMakeLists.txt index 97f3ae0f4f3f7..6adbc1c43d277 100644 --- a/Detectors/TRD/CMakeLists.txt +++ b/Detectors/TRD/CMakeLists.txt @@ -9,6 +9,7 @@ # submit itself to any jurisdiction. add_subdirectory(base) +add_subdirectory(calibration) add_subdirectory(simulation) add_subdirectory(reconstruction) add_subdirectory(macros) diff --git a/Detectors/TRD/calibration/CMakeLists.txt b/Detectors/TRD/calibration/CMakeLists.txt new file mode 100644 index 0000000000000..ba46a857e7c48 --- /dev/null +++ b/Detectors/TRD/calibration/CMakeLists.txt @@ -0,0 +1,17 @@ +#Copyright CERN and copyright holders of ALICE O2.This software is distributed +#under the terms of the GNU General Public License v3(GPL Version 3), copied +#verbatim in the file "COPYING". +# +#See http: //alice-o2.web.cern.ch/license for full licensing information. +# +#In applying this license CERN does not waive the privileges and immunities +#granted to it by virtue of its status as an Intergovernmental Organization or +#submit itself to any jurisdiction. + +o2_add_library(TRDCalibration + SOURCES src/CalibVDrift.cxx + PUBLIC_LINK_LIBRARIES O2::TRDBase + O2::DataFormatsTRD) + + o2_target_root_dictionary(TRDCalibration + HEADERS include/TRDCalibration/CalibVDrift.h) diff --git a/Detectors/TRD/calibration/include/TRDCalibration/CalibVDrift.h b/Detectors/TRD/calibration/include/TRDCalibration/CalibVDrift.h new file mode 100644 index 0000000000000..7e24a72b798b8 --- /dev/null +++ b/Detectors/TRD/calibration/include/TRDCalibration/CalibVDrift.h @@ -0,0 +1,68 @@ +// Copyright CERN and copyright holders of ALICE O2. This software is +// distributed under the terms of the GNU General Public License v3 (GPL +// Version 3), copied verbatim in the file "COPYING". +// +// See http://alice-o2.web.cern.ch/license for full licensing information. +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +#ifndef ALICEO2_TRD_CALIBVDRIFT_H_ +#define ALICEO2_TRD_CALIBVDRIFT_H_ + +/// \file CalibVDrift.h +/// \author Ole Schmidt, ole.schmidt@cern.ch + +#include "DataFormatsTRD/Constants.h" +#include + +namespace o2 +{ +namespace trd +{ + +/// \brief VDrift calibration class +/// +/// This class is used to determine chamber-wise vDrift values +/// +/// origin: TRD +/// \author Ole Schmidt, ole.schmidt@cern.ch + +class CalibVDrift +{ + public: + /// default constructor + CalibVDrift() = default; + + /// default destructor + ~CalibVDrift() = default; + + /// set input angular difference sums + void setAngleDiffSums(float* input) + { + for (int i = 0; i < constants::MAXCHAMBER * constants::NBINSANGLEDIFF; ++i) { + mAngleDiffSums[i] = input[i]; + } + } + + /// set input angular difference bin counters + void setAngleDiffCounters(short* input) + { + for (int i = 0; i < constants::MAXCHAMBER * constants::NBINSANGLEDIFF; ++i) { + mAngleDiffCounters[i] = input[i]; + } + } + + /// main processing function + void process(); + + private: + std::array mAngleDiffSums{}; ///< input TRD track to tracklet angular difference sums per bin + std::array mAngleDiffCounters{}; ///< input bin counters +}; + +} // namespace trd + +} // namespace o2 +#endif diff --git a/Detectors/TRD/calibration/src/CalibVDrift.cxx b/Detectors/TRD/calibration/src/CalibVDrift.cxx new file mode 100644 index 0000000000000..65076f82d21df --- /dev/null +++ b/Detectors/TRD/calibration/src/CalibVDrift.cxx @@ -0,0 +1,50 @@ +// Copyright CERN and copyright holders of ALICE O2. This software is +// distributed under the terms of the GNU General Public License v3 (GPL +// Version 3), copied verbatim in the file "COPYING". +// +// See http://alice-o2.web.cern.ch/license for full licensing information. +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +/// \file CalibVDrift.cxx +/// \author Ole Schmidt, ole.schmidt@cern.ch + +#include "TFile.h" +#include "TH2F.h" + +#include + +#include "TRDCalibration/CalibVDrift.h" + +using namespace o2::trd; +using namespace o2::trd::constants; + +void CalibVDrift::process() +{ + LOG(info) << "Started processing for vDrift calibration"; + + for (int iDet = 0; iDet < constants::MAXCHAMBER; ++iDet) { + for (int iBin = 0; iBin < constants::NBINSANGLEDIFF; ++iBin) { // note: iBin = constants::NBINSANGLEDIFF - 1 is under-/overflow bin + short nEntries = mAngleDiffCounters[iDet * constants::NBINSANGLEDIFF + iBin]; + float angleDiffSum = mAngleDiffSums[iDet * constants::NBINSANGLEDIFF + iBin]; + if (nEntries > 0) { + LOGF(INFO, "Found %i entrie(s) in chamber %i, bin %i. Average angular deviation: %f", nEntries, iDet, iBin, angleDiffSum / nEntries); + } + } + } + + /* + // as an example I loop over the input, create a histogram and write it to a file + auto fOut = TFile::Open("trdcalibdummy.root", "recreate"); + auto hXY = std::make_unique("histDummy", "foo", 100, -60, 60, 100, 250, 400); // xy distribution of TRD space points + for (int i = 0; i < mAngulerDeviationProf.size(); ++i) { + hXY->Fill(mAngulerDeviationProf[i].mX[0], mAngulerDeviationProf[i].mR); + } + fOut->cd(); + hXY->Write(); + hXY.reset(); // delete the histogram before closing the output file + fOut->Close(); + */ +} diff --git a/Detectors/TRD/calibration/src/TRDCalibrationLinkDef.h b/Detectors/TRD/calibration/src/TRDCalibrationLinkDef.h new file mode 100644 index 0000000000000..6f10a82a50b64 --- /dev/null +++ b/Detectors/TRD/calibration/src/TRDCalibrationLinkDef.h @@ -0,0 +1,19 @@ +// Copyright CERN and copyright holders of ALICE O2. This software is +// distributed under the terms of the GNU General Public License v3 (GPL +// Version 3), copied verbatim in the file "COPYING". +// +// See http://alice-o2.web.cern.ch/license for full licensing information. +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +#ifdef __CLING__ + +#pragma link off all globals; +#pragma link off all classes; +#pragma link off all functions; + +#pragma link C++ class o2::trd::CalibVDrift + ; + +#endif diff --git a/Detectors/TRD/workflow/CMakeLists.txt b/Detectors/TRD/workflow/CMakeLists.txt index 589932309585b..402102d857b2e 100644 --- a/Detectors/TRD/workflow/CMakeLists.txt +++ b/Detectors/TRD/workflow/CMakeLists.txt @@ -25,8 +25,8 @@ o2_add_library(TRDWorkflow src/TRDTrackWriterSpec.cxx src/TRDTrackingWorkflow.cxx src/EntropyDecoderSpec.cxx - src/EntropyEncoderSpec.cxx - PUBLIC_LINK_LIBRARIES O2::Framework O2::DPLUtils O2::Steer O2::Algorithm O2::DataFormatsTRD O2::TRDSimulation O2::TRDReconstruction O2::DetectorsBase O2::SimulationDataFormat O2::TRDBase O2::GPUTracking O2::GlobalTrackingWorkflow) + src/EntropyEncoderSpec.cxx + PUBLIC_LINK_LIBRARIES O2::Framework O2::DPLUtils O2::Steer O2::Algorithm O2::DataFormatsTRD O2::TRDSimulation O2::TRDReconstruction O2::DetectorsBase O2::SimulationDataFormat O2::TRDBase O2::TRDCalibration O2::GPUTracking O2::GlobalTrackingWorkflow) #o2_target_root_dictionary(TRDWorkflow # HEADERS include/TRDWorkflow/TRDTrapSimulatorSpec.h) diff --git a/Detectors/TRD/workflow/include/TRDWorkflow/TRDGlobalTrackingSpec.h b/Detectors/TRD/workflow/include/TRDWorkflow/TRDGlobalTrackingSpec.h index 3613cf620ca15..ec7bb972b8bec 100644 --- a/Detectors/TRD/workflow/include/TRDWorkflow/TRDGlobalTrackingSpec.h +++ b/Detectors/TRD/workflow/include/TRDWorkflow/TRDGlobalTrackingSpec.h @@ -17,6 +17,7 @@ #include "Framework/Task.h" #include "TStopwatch.h" #include "TRDBase/GeometryFlat.h" +#include "TRDCalibration/CalibVDrift.h" #include "GPUO2Interface.h" #include "GPUTRDTracker.h" @@ -41,6 +42,7 @@ class TRDGlobalTracking : public o2::framework::Task std::unique_ptr mFlatGeo{nullptr}; ///< flat TRD geometry bool mUseMC{false}; ///< MC flag bool mUseTrackletTransform{false}; ///< if true, output from TrackletTransformer is used instead of uncalibrated Tracklet64 directly + CalibVDrift mCalibVDrift{}; ///< steers the vDrift calibration TStopwatch mTimer; }; diff --git a/Detectors/TRD/workflow/src/TRDGlobalTrackingSpec.cxx b/Detectors/TRD/workflow/src/TRDGlobalTrackingSpec.cxx index 57a79011f8bb8..c7a3596fcece8 100644 --- a/Detectors/TRD/workflow/src/TRDGlobalTrackingSpec.cxx +++ b/Detectors/TRD/workflow/src/TRDGlobalTrackingSpec.cxx @@ -22,6 +22,7 @@ #include "DataFormatsTRD/Tracklet64.h" #include "DataFormatsTRD/CalibratedTracklet.h" #include "DataFormatsTRD/TriggerRecord.h" +#include "DataFormatsTRD/Constants.h" // GPU header #include "GPUReconstruction.h" @@ -68,6 +69,7 @@ void TRDGlobalTracking::init(InitContext& ic) mTracker->SetProcessPerTimeFrame(); mTracker->SetNMaxCollisions(mRec->GetProcessingSettings().trdNMaxCollisions); mTracker->SetTrkltTransformNeeded(!mUseTrackletTransform); + //mTracker->SetDoImpactAngleHistograms(true); mRec->RegisterGPUProcessor(mTracker, false); mChainTracking->SetTRDGeometry(std::move(mFlatGeo)); @@ -99,7 +101,7 @@ void TRDGlobalTracking::run(ProcessingContext& pc) int nTrackletsCal = 0; if (mUseTrackletTransform) { - cTrklts.emplace(pc.inputs().get>("trdctracklets")); // MC labels associated to the input digits + cTrklts.emplace(pc.inputs().get>("trdctracklets")); cTrkltsPtr = &cTrklts.value(); nTrackletsCal = cTrkltsPtr->size(); LOGF(INFO, "Got %i calibrated tracklets as input", nTrackletsCal); @@ -171,11 +173,18 @@ void TRDGlobalTracking::run(ProcessingContext& pc) mTracker->SetTriggerRecordIndices(&(trdTriggerIndices[0])); mTracker->SetNCollisions(nCollisions); //mTracker->DumpTracks(); + mTracker->ResetImpactAngleHistograms(); mTracker->DoTracking(mChainTracking); //mTracker->DumpTracks(); std::vector tracksOut(mTracker->NTracks()); std::copy(mTracker->Tracks(), mTracker->Tracks() + mTracker->NTracks(), tracksOut.begin()); + + // Temporary until it is transferred to its own DPL device for calibrations + mCalibVDrift.setAngleDiffSums(mTracker->AngleDiffSums()); + mCalibVDrift.setAngleDiffCounters(mTracker->AngleDiffCounters()); + mCalibVDrift.process(); + pc.outputs().snapshot(Output{o2::header::gDataOriginTRD, "MATCHTRD", 0, Lifetime::Timeframe}, tracksOut); mTimer.Stop(); diff --git a/GPU/GPUTracking/TRDTracking/GPUTRDInterfaces.h b/GPU/GPUTracking/TRDTracking/GPUTRDInterfaces.h index f6d71d93c60f6..e5a234d6b6eca 100644 --- a/GPU/GPUTracking/TRDTracking/GPUTRDInterfaces.h +++ b/GPU/GPUTracking/TRDTracking/GPUTRDInterfaces.h @@ -301,7 +301,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 resetCovariance(float s) { ResetCovariance(); } GPUd() void setAlpha(float alpha) { mAlpha = alpha; } GPUd() void set(float x, float alpha, const float param[5], const float cov[15]) { diff --git a/GPU/GPUTracking/TRDTracking/GPUTRDTracker.cxx b/GPU/GPUTracking/TRDTracking/GPUTRDTracker.cxx index 838fec471e870..fec9f3f3ff921 100644 --- a/GPU/GPUTracking/TRDTracking/GPUTRDTracker.cxx +++ b/GPU/GPUTracking/TRDTracking/GPUTRDTracker.cxx @@ -75,6 +75,8 @@ void* GPUTRDTracker_t::SetPointersBase(void* base) computePointerWithAlignment(base, mTrackletIndexArray, (kNChambers + 1) * mNMaxCollisions); computePointerWithAlignment(base, mHypothesis, mNCandidates * mMaxThreads); computePointerWithAlignment(base, mCandidates, mNCandidates * 2 * mMaxThreads); + computePointerWithAlignment(base, mAngleDiffSums, kNChambers * (mNAngleHistogramBins + 1)); + computePointerWithAlignment(base, mAngleDiffCounters, kNChambers * (mNAngleHistogramBins + 1)); return base; } @@ -103,7 +105,7 @@ void* GPUTRDTracker_t::SetPointersTracks(void* base) } template -GPUTRDTracker_t::GPUTRDTracker_t() : mR(nullptr), mIsInitialized(false), mTrkltTransfNeeded(true), 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), mTrackletIndices(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), mMaxEta(0.84f), mExtraRoadY(2.f), mRoadZ(18.f), mZCorrCoefNRC(1.4f), mMCEvent(nullptr), mDebug(new GPUTRDTrackerDebug()) +GPUTRDTracker_t::GPUTRDTracker_t() : mR(nullptr), mIsInitialized(false), mTrkltTransfNeeded(true), mDoImpactAngleHistograms(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), mTrackletIndices(nullptr), mMaxThreads(100), mNTracklets(0), mTrackletIndexArray(nullptr), mHypothesis(nullptr), mNAngleHistogramBins(25), mAngleHistogramRange(50), mCandidates(nullptr), mSpacePoints(nullptr), mAngleDiffSums(nullptr), mAngleDiffCounters(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), mMaxEta(0.84f), mExtraRoadY(2.f), mRoadZ(18.f), mZCorrCoefNRC(1.4f), mMCEvent(nullptr), mDebug(new GPUTRDTrackerDebug()) { //-------------------------------------------------------------------- // Default constructor @@ -272,6 +274,18 @@ void GPUTRDTracker_t::DoTracking(GPUChainTracking* chainTracking) #endif } + if (mDoImpactAngleHistograms) { + GPUInfo("Start calculating angular differences"); + for (int iTrk = 0; iTrk < mNTracks; ++iTrk) { + if (mTracks[iTrk].GetNtracklets() > 3) { + auto trkCopy = mTracks[iTrk]; + PROP prop(&Param().polynomialField); + prop.setFitInProjections(true); + prop.setTrack(&trkCopy); + FillImpactAngleHistograms(&prop, &trkCopy); + } + } + } auto duration = std::chrono::high_resolution_clock::now() - timeStart; (void)duration; // suppress warning about unused variable /* @@ -605,6 +619,125 @@ GPUd() bool GPUTRDTracker_t::CalculateSpacePoints(int iCollision) } return result; } +template +GPUd() void GPUTRDTracker_t::ResetImpactAngleHistograms() +{ + GPUInfo("Resetting angle histograms"); + for (int i = 0; i < kNChambers * (mNAngleHistogramBins + 1); i++) { + mAngleDiffSums[i] = 0; + mAngleDiffCounters[i] = 0; + } +} + +template +GPUd() int GPUTRDTracker_t::FillImpactAngleHistograms(PROP* prop, TRDTRK* t) +{ + //-------------------------------------------------------------------- + // To calibrate vDrift and ExB based on the online tracklets this function + // collects the differences between TRD-only tracks and their associated + // tracklets + // returns 0 in case of success + //-------------------------------------------------------------------- + float invBinWidth = mNAngleHistogramBins / mAngleHistogramRange; + t->SetChi2(0.f); + t->resetCovariance(100); + + // first inward propagation (TRD track fit) + for (int iLayer = kNLayers - 1; iLayer >= 0; --iLayer) { + if (t->GetTracklet(iLayer) == -1) { + continue; + } + if (PropagateToLayerAndUpdate(prop, t, iLayer)) { + return 1; + } + } + + // outward propagation (smoothing) + for (int iLayer = 1; iLayer < kNLayers; ++iLayer) { + if (t->GetTracklet(iLayer) == -1) { + continue; + } + if (PropagateToLayerAndUpdate(prop, t, iLayer)) { + return 2; + } + } + + // second inward propagation (collect angular differences between tracklets + TRD track) + for (int iLayer = kNLayers - 1; iLayer >= 0; --iLayer) { + if (t->GetTracklet(iLayer) == -1) { + continue; + } + if (PropagateToLayerAndUpdate(prop, t, iLayer, false)) { + return 3; + } + float radToDeg = 180.f / CAMath::Pi(); + float trkAngle = CAMath::ASin(t->getSnp()) * radToDeg; + float trkltAngle = CAMath::ATan(mSpacePoints[t->GetTracklet(iLayer)].mDy / 3.) * radToDeg; + + int idxOffsetDet = mTracklets[t->GetTracklet(iLayer)].GetDetector() * (mNAngleHistogramBins + 1); + int idxOffsetAngle = (trkAngle + .5 * mAngleHistogramRange) * invBinWidth; + + if (CAMath::Abs(idxOffsetAngle) >= .5 * mAngleHistogramRange) { + idxOffsetAngle = mNAngleHistogramBins; + } + + mAngleDiffSums[idxOffsetDet + idxOffsetAngle] += trkltAngle - trkAngle; + mAngleDiffCounters[idxOffsetDet + idxOffsetAngle]++; + + //GPUInfo("trkAngle(%f), idxOffsetAngle(%i), angleDifference(%f)", trkAngle, idxOffsetAngle, trkltAngle - trkAngle); + } + return 0; +} + +template +GPUd() int GPUTRDTracker_t::PropagateToLayerAndUpdate(PROP* prop, TRDTRK* trkWork, int iLayer, bool doUpdate) +{ + //-------------------------------------------------------------------- + // Propagates the track to TRD layer iLayer and updates the track + // parameters (if requested) + // returns 0 in case of success + //-------------------------------------------------------------------- + int trackletID = trkWork->GetTracklet(iLayer); + int trackletDet = mTracklets[trackletID].GetDetector(); + int trackletSector = trackletDet / (kNLayers * kNStacks); + int trackletStack = (trackletDet % (kNLayers * kNStacks)) / kNLayers; + + if (trackletSector != GetSector(prop->getAlpha())) { + if (!prop->rotate(GetAlphaOfSector(trackletSector))) { + GPUInfo("Track could not be rotated in tracklet coordinate system"); + return 1; + } + } + + if (!prop->propagateToX(mSpacePoints[trackletID].mR, .8f, 2.f)) { + GPUInfo("Track propagation failed in layer %i (pt=%f, xTrk=%f, xToGo=%f)", iLayer, trkWork->getPt(), trkWork->getX(), mSpacePoints[trackletID].mR); + return 2; + } + + if (!doUpdate) { + // nothing more to be done + return 0; + } + + const GPUTRDpadPlane* pad = mGeo->GetPadPlane(iLayer, trackletStack); + float tilt = CAMath::Tan(CAMath::Pi() / 180.f * pad->GetTiltingAngle()); // tilt is signed! + float tiltCorrUp = tilt * (mSpacePoints[trackletID].mX[1] - trkWork->getZ()); + float zPosCorrUp = mSpacePoints[trackletID].mX[1] + mZCorrCoefNRC * trkWork->getTgl(); + float padLength = pad->GetRowSize(mTracklets[trackletID].GetZbin()); + if (!((trkWork->getSigmaZ2() < (padLength * padLength / 12.f)) && (CAMath::Abs(mSpacePoints[trackletID].mX[1] - trkWork->getZ()) < padLength))) { + tiltCorrUp = 0.f; + } + + My_Float trkltPosUp[2] = {mSpacePoints[trackletID].mX[0] - tiltCorrUp, zPosCorrUp}; + My_Float trkltCovUp[3] = {0.f}; + RecalcTrkltCov(tilt, trkWork->getSnp(), pad->GetRowSize(mTracklets[trackletID].GetZbin()), trkltCovUp); + + if (!prop->update(trkltPosUp, trkltCovUp)) { + GPUWarning("Failed to update track with space point in layer %i", iLayer); + return 3; + } + return 0; +} template GPUd() bool GPUTRDTracker_t::FollowProlongation(PROP* prop, TRDTRK* t, int threadId, int collisionId) @@ -733,6 +866,7 @@ GPUd() bool GPUTRDTracker_t::FollowProlongation(PROP* prop, TRDTRK if (currDet == -1) { continue; } + pad = mGeo->GetPadPlane(currDet); int currSec = mGeo->GetSector(currDet); if (currSec != GetSector(prop->getAlpha())) { if (!prop->rotate(GetAlphaOfSector(currSec))) { @@ -762,10 +896,10 @@ GPUd() bool GPUTRDTracker_t::FollowProlongation(PROP* prop, TRDTRK } float projY, projZ; prop->getPropagatedYZ(mSpacePoints[trkltIdxGlb].mR, projY, projZ); - // correction for tilted pads (only applied if deltaZ < l_pad && track z err << l_pad) + // correction for tilted pads (only applied if deltaZ < lPad && track z err << lPad) float tiltCorr = tilt * (mSpacePoints[trkltIdxGlb].mX[1] - projZ); - float l_pad = pad->GetRowSize(mTracklets[trkltIdxGlb].GetZbin()); - if (!((CAMath::Abs(mSpacePoints[trkltIdxGlb].mX[1] - projZ) < l_pad) && (trkWork->getSigmaZ2() < (l_pad * l_pad / 12.f)))) { + float lPad = pad->GetRowSize(mTracklets[trkltIdxGlb].GetZbin()); + if (!((CAMath::Abs(mSpacePoints[trkltIdxGlb].mX[1] - projZ) < lPad) && (trkWork->getSigmaZ2() < (lPad * lPad / 12.f)))) { tiltCorr = 0.f; } // correction for mean z position of tracklet (is not the center of the pad if track eta != 0) @@ -800,6 +934,7 @@ GPUd() bool GPUTRDTracker_t::FollowProlongation(PROP* prop, TRDTRK mDebug->SetNmatchAvail(matchAvailableAll[iLayer].size(), iLayer); int realTrkltId = matchAvailableAll[iLayer].at(0); int realTrkltDet = mTracklets[realTrkltId].GetDetector(); + pad = mGeo->GetPadPlane(realTrkltDet); prop->rotate(GetAlphaOfSector(mGeo->GetSector(realTrkltDet))); if (!prop->propagateToX(mSpacePoints[realTrkltId].mR, .8f, 2.f) || GetSector(prop->getAlpha()) != mGeo->GetSector(realTrkltDet)) { if (ENABLE_WARNING) { @@ -811,8 +946,8 @@ GPUd() bool GPUTRDTracker_t::FollowProlongation(PROP* prop, TRDTRK float zPosCorrReal = mSpacePoints[realTrkltId].mX[1] + mZCorrCoefNRC * trkWork->getTgl(); float deltaZReal = zPosCorrReal - trkWork->getZ(); float tiltCorrReal = tilt * (mSpacePoints[realTrkltId].mX[1] - trkWork->getZ()); - float l_padReal = pad->GetRowSize(mTracklets[realTrkltId].GetZbin()); - if ((trkWork->getSigmaZ2() >= (l_padReal * l_padReal / 12.f)) || (CAMath::Abs(mSpacePoints[realTrkltId].mX[1] - trkWork->getZ()) >= l_padReal)) { + float lPadReal = pad->GetRowSize(mTracklets[realTrkltId].GetZbin()); + if ((trkWork->getSigmaZ2() >= (lPadReal * lPadReal / 12.f)) || (CAMath::Abs(mSpacePoints[realTrkltId].mX[1] - trkWork->getZ()) >= lPadReal)) { tiltCorrReal = 0; } My_Float yzPosReal[2] = {mSpacePoints[realTrkltId].mX[0] - tiltCorrReal, zPosCorrReal}; @@ -890,10 +1025,11 @@ GPUd() bool GPUTRDTracker_t::FollowProlongation(PROP* prop, TRDTRK continue; } + pad = mGeo->GetPadPlane(mTracklets[mHypothesis[iUpdate + hypothesisIdxOffset].mTrackletId].GetDetector()); float tiltCorrUp = tilt * (mSpacePoints[mHypothesis[iUpdate + hypothesisIdxOffset].mTrackletId].mX[1] - trkWork->getZ()); float zPosCorrUp = mSpacePoints[mHypothesis[iUpdate + hypothesisIdxOffset].mTrackletId].mX[1] + mZCorrCoefNRC * trkWork->getTgl(); - float l_padTrklt = pad->GetRowSize(mTracklets[mHypothesis[iUpdate + hypothesisIdxOffset].mTrackletId].GetZbin()); - if (!((trkWork->getSigmaZ2() < (l_padTrklt * l_padTrklt / 12.f)) && (CAMath::Abs(mSpacePoints[mHypothesis[iUpdate + hypothesisIdxOffset].mTrackletId].mX[1] - trkWork->getZ()) < l_padTrklt))) { + float padLength = pad->GetRowSize(mTracklets[mHypothesis[iUpdate + hypothesisIdxOffset].mTrackletId].GetZbin()); + if (!((trkWork->getSigmaZ2() < (padLength * padLength / 12.f)) && (CAMath::Abs(mSpacePoints[mHypothesis[iUpdate + hypothesisIdxOffset].mTrackletId].mX[1] - trkWork->getZ()) < padLength))) { tiltCorrUp = 0.f; } My_Float trkltPosUp[2] = {mSpacePoints[mHypothesis[iUpdate + hypothesisIdxOffset].mTrackletId].mX[0] - tiltCorrUp, zPosCorrUp}; diff --git a/GPU/GPUTracking/TRDTracking/GPUTRDTracker.h b/GPU/GPUTracking/TRDTracking/GPUTRDTracker.h index ea3bbbc9cd060..da91005d9370d 100644 --- a/GPU/GPUTracking/TRDTracking/GPUTRDTracker.h +++ b/GPU/GPUTracking/TRDTracking/GPUTRDTracker.h @@ -122,6 +122,9 @@ 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() int FillImpactAngleHistograms(PROP* prop, TRDTRK* t); + GPUd() void ResetImpactAngleHistograms(); + GPUd() int PropagateToLayerAndUpdate(PROP* prop, TRDTRK* trkWork, int iLayer, bool doUpdate = true); GPUd() int GetDetectorNumber(const float zPos, const float alpha, const int layer) const; GPUd() bool AdjustSector(PROP* prop, TRDTRK* t) const; GPUd() int GetSector(float alpha) const; @@ -146,6 +149,7 @@ class GPUTRDTracker_t : public GPUProcessor // settings GPUd() void SetTrkltTransformNeeded(bool flag) { mTrkltTransfNeeded = flag; } GPUd() void SetProcessPerTimeFrame() { mProcessPerTimeFrame = true; } + GPUd() void SetDoImpactAngleHistograms(bool flag) { mDoImpactAngleHistograms = flag; } GPUd() void SetMCEvent(AliMCEvent* mc) { mMCEvent = mc; } GPUd() void EnableDebugOutput() { mDebugOutput = true; } GPUd() void SetMaxEta(float maxEta) { mMaxEta = maxEta; } @@ -164,6 +168,8 @@ class GPUTRDTracker_t : public GPUProcessor GPUd() TRDTRK* Tracks() const { return mTracks; } GPUd() int NTracklets() const { return mNTracklets; } GPUd() GPUTRDSpacePointInternal* SpacePoints() const { return mSpacePoints; } + GPUd() float* AngleDiffSums() const { return mAngleDiffSums; } + GPUd() short* AngleDiffCounters() const { return mAngleDiffCounters; } GPUd() GPUTRDTrackletWord* Tracklets() const { return mTracklets; } GPUd() void DumpTracks(); @@ -172,6 +178,9 @@ class GPUTRDTracker_t : public GPUProcessor bool mIsInitialized; // flag is set upon initialization bool mTrkltTransfNeeded; // if the output of the TRDTrackletTransformer is used we don't need to do the coordinate transformation for the tracklets bool mProcessPerTimeFrame; // if true, tracking is done per time frame instead of on a single events basis + bool mDoImpactAngleHistograms; // if true, impact angle vs angle difference histograms are filled + short mNAngleHistogramBins; // number of bins per chamber for the angular difference histograms + float mAngleHistogramRange; // range of impact angles covered by each histogram 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) @@ -196,6 +205,8 @@ class GPUTRDTracker_t : public GPUProcessor Hypothesis* mHypothesis; // array with multiple track hypothesis TRDTRK* mCandidates; // array of tracks for multiple hypothesis tracking GPUTRDSpacePointInternal* mSpacePoints; // array with tracklet coordinates in global tracking frame + float* mAngleDiffSums; // array with sum of angular differences for a given bin + short* mAngleDiffCounters; // array with number of entries for a given bin int* mTrackletLabels; // array with MC tracklet labels TRD_GEOMETRY_CONST GPUTRDGeometry* mGeo; // TRD geometry /// ---- error parametrization depending on magnetic field ----