From 736c321282f0e00626ccf46dadb46e59c077d7a0 Mon Sep 17 00:00:00 2001 From: Jason Barrella Date: Thu, 25 Feb 2021 11:06:21 +0200 Subject: [PATCH] new coordinate transformation functions in PadPlane with CMake tests --- Detectors/TRD/base/CMakeLists.txt | 7 + Detectors/TRD/base/include/TRDBase/PadPlane.h | 42 +++++- .../include/TRDBase/TrackletTransformer.h | 2 + Detectors/TRD/base/src/PadPlane.cxx | 58 +++++++- .../TRD/base/src/TrackletTransformer.cxx | 20 +++ .../base/test/testCoordinateTransforms.cxx | 135 ++++++++++++++++++ Detectors/TRD/simulation/src/Digitizer.cxx | 4 +- 7 files changed, 259 insertions(+), 9 deletions(-) create mode 100644 Detectors/TRD/base/test/testCoordinateTransforms.cxx diff --git a/Detectors/TRD/base/CMakeLists.txt b/Detectors/TRD/base/CMakeLists.txt index 3b7c8eb85bee7..eccd1f494a909 100644 --- a/Detectors/TRD/base/CMakeLists.txt +++ b/Detectors/TRD/base/CMakeLists.txt @@ -86,3 +86,10 @@ o2_add_test(RawData ENVIRONMENT O2_ROOT=${CMAKE_BINARY_DIR}/stage LABELS trd ) +o2_add_test(Transformations + COMPONENT_NAME trd + PUBLIC_LINK_LIBRARIES O2::TRDBase O2::DataFormatsTRD + SOURCES test/testCoordinateTransforms.cxx + ENVIRONMENT O2_ROOT=${CMAKE_BINARY_DIR}/stage + LABELS trd + ) diff --git a/Detectors/TRD/base/include/TRDBase/PadPlane.h b/Detectors/TRD/base/include/TRDBase/PadPlane.h index 7e72d508c0e03..c81053a97cee4 100644 --- a/Detectors/TRD/base/include/TRDBase/PadPlane.h +++ b/Detectors/TRD/base/include/TRDBase/PadPlane.h @@ -13,7 +13,6 @@ //Forwards to standard header with protection for GPU compilation #include "GPUCommonRtypes.h" // for ClassDef - #include "GPUCommonDef.h" //////////////////////////////////////////////////////////////////////////// @@ -61,19 +60,44 @@ class PadPlane }; void setLength(double l) { mLength = l; }; void setWidth(double w) { mWidth = w; }; - void setLengthOPad(double l) { mLengthOPad = l; }; - void setWidthOPad(double w) { mWidthOPad = w; }; - void setLengthIPad(double l) { mLengthIPad = l; }; - void setWidthIPad(double w) { mWidthIPad = w; }; + void setLengthOPad(double l) + { + mLengthOPad = l; + mInverseLengthOPad = 1.0 / l; + }; + void setWidthOPad(double w) + { + mWidthOPad = w; + mInverseWidthOPad = 1.0 / w; + }; + void setLengthIPad(double l) + { + mLengthIPad = l; + mInverseLengthIPad = 1.0 / l; + }; + void setWidthIPad(double w) + { + mWidthIPad = w; + mInverseWidthIPad = 1.0 / w; + }; void setPadRowSMOffset(double o) { mPadRowSMOffset = o; }; void setAnodeWireOffset(float o) { mAnodeWireOffset = o; }; void setTiltingAngle(double t); GPUd() int getPadRowNumber(double z) const; GPUd() int getPadRowNumberROC(double z) const; + GPUd() double getPadRow(double z) const; GPUd() int getPadColNumber(double rphi) const; + GPUd() double getPad(double y, double z) const; - GPUd() double getTiltOffset(double rowOffset) const { return mTiltingTan * (rowOffset - 0.5 * mLengthIPad); }; + GPUd() double getTiltOffset(int row, double rowOffset) const + { + if (row == 0 || row == mNrows - 1) { + return mTiltingTan * (rowOffset - 0.5 * mLengthOPad); + } else { + return mTiltingTan * (rowOffset - 0.5 * mLengthIPad); + } + }; GPUd() double getPadRowOffset(int row, double z) const { if ((row < 0) || (row >= mNrows)) { @@ -174,6 +198,12 @@ class PadPlane double mAnodeWireOffset; // Distance of first anode wire from pad edge + double mInverseLengthIPad; // 1 / mLengthIPad + double mInverseLengthOPad; // 1 / mLengthOPad + + double mInverseWidthIPad; // 1 / mWidthIPad + double mInverseWidthOPad; // 1 / mWidthOPad + private: ClassDefNV(PadPlane, 1); // TRD ROC pad plane }; diff --git a/Detectors/TRD/base/include/TRDBase/TrackletTransformer.h b/Detectors/TRD/base/include/TRDBase/TrackletTransformer.h index 80614fb2df6d8..51588c36146b4 100644 --- a/Detectors/TRD/base/include/TRDBase/TrackletTransformer.h +++ b/Detectors/TRD/base/include/TRDBase/TrackletTransformer.h @@ -50,6 +50,8 @@ class TrackletTransformer CalibratedTracklet transformTracklet(Tracklet64 tracklet); + double getTimebin(double x); + private: o2::trd::Geometry* mGeo; const o2::trd::PadPlane* mPadPlane; diff --git a/Detectors/TRD/base/src/PadPlane.cxx b/Detectors/TRD/base/src/PadPlane.cxx index bc8a65c485869..cbbf67258d5cd 100644 --- a/Detectors/TRD/base/src/PadPlane.cxx +++ b/Detectors/TRD/base/src/PadPlane.cxx @@ -24,8 +24,10 @@ #include "TRDBase/PadPlane.h" #include #include +#include "DataFormatsTRD/Constants.h" using namespace o2::trd; +using namespace o2::trd::constants; //_____________________________________________________________________________ void PadPlane::setTiltingAngle(double t) @@ -35,7 +37,7 @@ void PadPlane::setTiltingAngle(double t) // mTiltingAngle = t; - mTiltingTan = TMath::Tan(TMath::Pi() / 180.0 * mTiltingAngle); + mTiltingTan = TMath::Tan(TMath::DegToRad() * mTiltingAngle); } //_____________________________________________________________________________ @@ -158,3 +160,57 @@ void PadPlane::setNrows(int n) } mNrows = n; }; + +double PadPlane::getPadRow(double z) const +{ + double lengthCorr = mLengthIPad * mInverseLengthOPad; + + // calculate position based on inner pad length + double padrow = -z * mInverseLengthIPad + mNrows * 0.5; + + // correct row for outer pad rows + if (padrow <= 1.0) { + padrow = 1.0 - (1.0 - padrow) * lengthCorr; + } + + if (padrow >= double(mNrows - 1)) { + padrow = double(mNrows - 1) + (padrow - double(mNrows - 1)) * lengthCorr; + } + + // sanity check: is the padrow coordinate reasonable? + assert(!(padrow < 0.0 || padrow > double(mNrows))); + + return padrow; +} + +double PadPlane::getPad(double y, double z) const +{ + int padrow = getPadRow(z); + double padrowOffset = getPadRowOffsetROC(padrow, z); + double tiltOffsetY = getTiltOffset(padrow, padrowOffset); + + double pad = y * mInverseWidthIPad + mNcols * 0.5; + + double lengthCorr = mWidthIPad * mInverseWidthOPad; + // correct row for outer pad rows + if (pad <= 1.0) { + pad = 1.0 - (1.0 - pad) * lengthCorr; + } + + if (pad >= double(mNcols - 1)) { + pad = double(mNcols - 1) + (pad - double(mNcols - 1)) * lengthCorr; + } + + double tiltOffsetPad; + if (pad <= 1.0 || pad >= double(mNcols - 1)) { + tiltOffsetPad = tiltOffsetY * mInverseWidthOPad; + pad += tiltOffsetPad; + } else { + tiltOffsetPad = tiltOffsetY * mInverseWidthIPad; + pad += tiltOffsetPad; + } + + assert(!(pad < 0.0 || pad > double(mNcols))); + + return pad; +} diff --git a/Detectors/TRD/base/src/TrackletTransformer.cxx b/Detectors/TRD/base/src/TrackletTransformer.cxx index 741e6d85ec5fc..a1a546a4faad6 100644 --- a/Detectors/TRD/base/src/TrackletTransformer.cxx +++ b/Detectors/TRD/base/src/TrackletTransformer.cxx @@ -165,3 +165,23 @@ CalibratedTracklet TrackletTransformer::transformTracklet(Tracklet64 tracklet) return CalibratedTracklet(sectorSpacePoint[0], sectorSpacePoint[1], sectorSpacePoint[2], dy); } + +double TrackletTransformer::getTimebin(double x) +{ + // calculate timebin from x position within chamber + // calibration parameters need to be extracted from CCDB in the future + double vDrift = 1.5625; // in cm/us + double t0 = 4.0; // time (in timebins) of start of drift region + + double timebin; + // x = 0 at anode plane and points toward pad plane. + if (x < -mGeo->camHght() / 2) { + // drift region + timebin = t0 - (x + mGeo->camHght() / 2) / (vDrift * 0.1); + } else { + // anode region: very rough guess + timebin = t0 - 1.0 + fabs(x); + } + + return timebin; +} diff --git a/Detectors/TRD/base/test/testCoordinateTransforms.cxx b/Detectors/TRD/base/test/testCoordinateTransforms.cxx new file mode 100644 index 0000000000000..1bcfbd56314ed --- /dev/null +++ b/Detectors/TRD/base/test/testCoordinateTransforms.cxx @@ -0,0 +1,135 @@ +// 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 testCoordinateTransformscxx +/// \brief Test local to row-column (float) coordinate transformations in PadPlane class +/// \author Jason Barrella - jbarrell@cern.ch, Sean Murray - murrays@cern.ch + +#define BOOST_TEST_MODULE Test CoordinateTransforms +#define BOOST_TEST_MAIN +#define BOOST_TEST_DYN_LINK +#include +#include +#include + +#include "DataFormatsTRD/Constants.h" +#include "TRDBase/TrackletTransformer.h" + +#include "TRDBase/PadPlane.h" +#include "DetectorsBase/GeometryManager.h" +#include "TRDBase/Geometry.h" + +namespace o2 +{ +namespace trd +{ + +using namespace o2::trd::constants; +using namespace std; + +void testRCPoint(double calculatedPoint, double predictedPoint) +{ + float e = 0.000001; + BOOST_CHECK(fabs(predictedPoint - calculatedPoint) <= e); +} + +BOOST_AUTO_TEST_CASE(LocaltoRCTest) +{ + auto mGeo = o2::trd::Geometry::instance(); + mGeo->createPadPlaneArray(); + mGeo->createClusterMatrixArray(); + + int hcid = 776; + // This C1 chamber has 16 pad rows with I pad length = 90mm and O pad length = 75mm + int detector = hcid / 2; + int stack = mGeo->getStack(detector); + int layer = mGeo->getLayer(detector); + + auto padPlane = mGeo->getPadPlane(layer, stack); + double lengthIPad = padPlane->getLengthIPad(); + double lengthOPad = padPlane->getLengthOPad(); + + double padIWidth = padPlane->getWidthIPad(); + double padOWidth = padPlane->getWidthOPad(); + + double tiltingAngle = padPlane->getTiltingAngle(); + + // Test padrows + auto p1 = padPlane->getPadRow(0); + // Center of the chamber. This should return the lower edge of padrow 8. + // Since we are using float values, the lower edge of padrow 8 correspond with float value 8.0. + testRCPoint(p1, 8.); + + auto p2 = padPlane->getPadRow(lengthIPad / 2.); + // With an I pad length of 9 cm, 9. / 2. should put us half way into the preceeding pad (pad 7) since padrow number + // decreses in positive z direction. + testRCPoint(p2, 7.5); + + auto p3 = padPlane->getPadRow(-lengthIPad / 2.); + // Same as above but in the other direction. + testRCPoint(p3, 8.5); + + auto p4 = padPlane->getPadRow(lengthIPad * 4.2); + // Arbitrary distance in z. + testRCPoint(p4, 8 - 4.2); + + auto p5 = padPlane->getPadRow(lengthIPad * 7 + lengthOPad); + // Lower border case. Take center and add 7 pads * 9 cm + 1 pad * 7.5 cm. + testRCPoint(p5, 0); + + auto p6 = padPlane->getPadRow(-lengthIPad * 7 - lengthOPad); + // Upper border case. Take center and subtract 7 pads * 9 cm - 1 pad * 7.5 cm. + testRCPoint(p6, 16); + + // Test pads with pad tilting + double p13 = padPlane->getPad(0, 0.01); + // Center of chamber plus epsilon (in z). Note the discontinuity in pad vs. z at z=0. This puts us at a point on the lower + // (in z direction) end of padrow 7. Since we have a pad tilt of -2 deg (pads tilted clockwise for particle coming + // from interaction vertex). After a lot thought and doodles on Miro (https://miro.com/app/board/o9J_lKgybMc=/) we + // find that we expect a small negative offset which would place us in the upper half of pad 71. + // To calculate that offset, we multiply the tangent of the tilting angle by the distance that our point is away + // from the center of its local padrow. Some unit conversions are also neccessary... + testRCPoint(p13, 72 + TMath::Tan(TMath::DegToRad() * tiltingAngle) * (0.5 * lengthIPad - 0.01) / padIWidth); + + double p14 = padPlane->getPad(0, -lengthIPad / 2.); + // Move from center of chamber to to center of padrow 8 = 8.5. + // This should now place us as edge of pad 72 = 72.0 since no offset is applied from pad tilt. + testRCPoint(p14, 72); + + double p15 = padPlane->getPad(0, lengthIPad / 2.); + // Should be the same in the other direction at padrow = 7.5 + testRCPoint(p15, 72); + + double p16 = padPlane->getPad(padIWidth * 42, lengthIPad / 2.); + // Adding an arbitrary number of pads should just increase position by same number + testRCPoint(p16, 72 + 42); + + double p17 = padPlane->getPad(padIWidth * 42, -lengthIPad * 4 - 2.3); + // Moving to arbitrary point in both y and z. In y, we are 42 pads from the center which would be 72 + 42 = 114 + // if we were in the center of the padrow and no pad tilting is considered. However, we are not in the center of + // a padrow, but rather 2.3 cm into padrow 8 + 4 = 12. This puts us below the center which is at 4.5 cm + // into the padrow and therefore, we expect a small postiive offset since the pad tilting angle is -2 deg. + testRCPoint(p17, 72 + 42 + TMath::Tan(TMath::DegToRad() * tiltingAngle) * (2.3 - 0.5 * lengthIPad) / padIWidth); + + double p18 = padPlane->getPad(padIWidth * 71 + padOWidth, lengthIPad / 2.); + // Border case right on the upper edge of the padrow in y direction and at the center of the padrow in z + testRCPoint(p18, 144); + + double p19 = padPlane->getPad(-padIWidth * 71 - padOWidth, lengthIPad / 2.); + // Border case right on the lower edge of the padrow in y direction and at the center of the padrow in z + testRCPoint(p19, 0); + + double p20 = padPlane->getPad(0, lengthIPad * 7 + 1.5); + // Ensure that shorter length of outer padrows is considered correctly + testRCPoint(p20, 72 + TMath::Tan(TMath::DegToRad() * tiltingAngle) * (0.5 * lengthOPad - 1.5) / padIWidth); +} + +} // namespace trd +} // namespace o2 diff --git a/Detectors/TRD/simulation/src/Digitizer.cxx b/Detectors/TRD/simulation/src/Digitizer.cxx index 2053be4b39070..c51de1cf7f2db 100644 --- a/Detectors/TRD/simulation/src/Digitizer.cxx +++ b/Detectors/TRD/simulation/src/Digitizer.cxx @@ -251,7 +251,7 @@ bool Digitizer::convertHits(const int det, const std::vector& hits, SignalC } double rowOffset = padPlane->getPadRowOffsetROC(rowE, locR); - double offsetTilt = padPlane->getTiltOffset(rowOffset); + double offsetTilt = padPlane->getTiltOffset(rowE, rowOffset); int colE = padPlane->getPadColNumber(locC + offsetTilt); if (colE < 0) { continue; @@ -295,7 +295,7 @@ bool Digitizer::convertHits(const int det, const std::vector& hits, SignalC } rowOffset = padPlane->getPadRowOffsetROC(rowE, locRd); // The pad column (rphi-direction) - offsetTilt = padPlane->getTiltOffset(rowOffset); + offsetTilt = padPlane->getTiltOffset(rowE, rowOffset); colE = padPlane->getPadColNumber(locCd + offsetTilt); if (colE < 0) { continue;