AliceO2Group · davidrohr · May 5, 2021 · Feb 25, 2021
@@ -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
+            )
@@ -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
 };

@@ -50,6 +50,8 @@ class TrackletTransformer
 
   CalibratedTracklet transformTracklet(Tracklet64 tracklet);
 
+  double getTimebin(double x);
+
  private:
   o2::trd::Geometry* mGeo;
   const o2::trd::PadPlane* mPadPlane;

@@ -24,8 +24,10 @@
 #include "TRDBase/PadPlane.h"
 #include <TMath.h>
 #include <FairLogger.h>
+#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;
+}
@@ -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;
+}
@@ -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 <boost/test/unit_test.hpp>
+#include <iostream>
+#include <numeric>
+
+#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
@@ -251,7 +251,7 @@ bool Digitizer::convertHits(const int det, const std::vector<Hit>& 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<Hit>& 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;