packages/tpc/TpcDistortionCorrection.cc

0001 /*!
0002  * \file TpcDistortionCorrection.cc
0003  * \brief applies provided distortion corrections to a cluster and returns corrected position
0004  * \author Hugo Pereira Da Costa <hugo.pereira-da-costa@cea.fr>
0005  */
0006
0007 #include "TpcDistortionCorrection.h"
0008 #include "TpcDistortionCorrectionContainer.h"
0009
0010 #include <TH1.h>
0011 #include <cmath>
0012
0013 #include <iostream>
0014
0015 namespace
0016 {
0017   template <class T>
0018   inline constexpr T square(const T& x)
0019   {
0020     return x * x;
0021   }
0022
0023   // check boundaries in axis
0024   /* for the interpolation to work, the value must be within the range of the provided axis, and not into the first and last bin */
0025   inline bool check_boundaries(const TAxis* axis, double value)
0026   {
0027     const auto bin = axis->FindBin(value);
0028     return (bin >= 2 && bin < axis->GetNbins());
0029   }
0030
0031   // check boundaries in histogram, before interpolation
0032   /* for the interpolation to work, the value must be within the range of the provided axis, and not into the first and last bin */
0033   inline bool check_boundaries(const TH1* h, double r, double phi, double z)
0034   {
0035     return check_boundaries(h->GetXaxis(), r) && check_boundaries(h->GetYaxis(), phi) && check_boundaries(h->GetZaxis(), z);
0036   }
0037
0038   // check boundaries in histogram, before interpolation
0039   /* for the interpolation to work, the value must be within the range of the provided axis, and not into the first and last bin */
0040   inline bool check_boundaries(const TH1* h, double r, double phi)
0041   {
0042     return check_boundaries(h->GetXaxis(), r) && check_boundaries(h->GetYaxis(), phi);
0043   }
0044
0045 }  // namespace
0046
0047 //________________________________________________________
0048 Acts::Vector3 TpcDistortionCorrection::get_corrected_position(const Acts::Vector3& source, const TpcDistortionCorrectionContainer* dcc, unsigned int mask) const
0049 {
0050   // get cluster radius, phi and z
0051   const auto r = std::sqrt(square(source.x()) + square(source.y()));
0052   auto phi = std::atan2(source.y(), source.x());
0053   if (phi < 0)
0054   {
0055     phi += 2 * M_PI;
0056   }
0057
0058   const auto z = source.z();
0059   const int index = z > 0 ? 1 : 0;
0060
0061   // apply corrections
0062   auto phi_new = phi;
0063   auto r_new = r;
0064   auto z_new = z;
0065
0066   // if the phi correction hist units are cm, we must divide by r to get the dPhi in radians
0067   auto divisor = r;
0068
0069   if (dcc->m_phi_hist_in_radians)
0070   {
0071     // if the phi correction hist units are radians, we must not divide by r.
0072     divisor = 1.0;
0073   }
0074
0075   //set our default corrections to be zero:
0076   //first inherit the same type
0077   auto dphi=phi;
0078   auto dr=r;
0079   auto dz=z;
0080   //then set them to zero:
0081   dphi=0;
0082   dr=0;
0083   dz=0;
0084
0085   //get the corrections from the histograms
0086   if (dcc->m_dimensions == 3)
0087   {
0088     if (dcc->m_hDPint[index] && (mask & COORD_PHI) && check_boundaries(dcc->m_hDPint[index], phi, r, z))
0089     {
0090       dphi=dcc->m_hDPint[index]->Interpolate(phi, r, z) / divisor;
0091     }
0092     if (dcc->m_hDRint[index] && (mask & COORD_R) && check_boundaries(dcc->m_hDRint[index], phi, r, z))
0093     {
0094       dr=dcc->m_hDRint[index]->Interpolate(phi, r, z);
0095     }
0096     if (dcc->m_hDZint[index] && (mask & COORD_Z) && check_boundaries(dcc->m_hDZint[index], phi, r, z))
0097     {
0098       dz=dcc->m_hDZint[index]->Interpolate(phi, r, z);
0099     }
0100   }
0101   else if (dcc->m_dimensions == 2)
0102   {
0103     double zterm = 1.0;
0104
0105     if (dcc->m_interpolate_z){
0106       zterm=(1. - std::abs(z) / 102.605);
0107     }
0108     if (dcc->m_hDPint[index] && (mask & COORD_PHI) && check_boundaries(dcc->m_hDPint[index], phi, r))
0109     {
0110       dphi=dcc->m_hDPint[index]->Interpolate(phi, r) * zterm / divisor;
0111     }
0112     if (dcc->m_hDRint[index] && (mask & COORD_R) && check_boundaries(dcc->m_hDRint[index], phi, r))
0113     {
0114       dr=dcc->m_hDRint[index]->Interpolate(phi, r) * zterm;
0115     }
0116     if (dcc->m_hDZint[index] && (mask & COORD_Z) && check_boundaries(dcc->m_hDZint[index], phi, r))
0117     {
0118       dz=dcc->m_hDZint[index]->Interpolate(phi, r) * zterm;
0119     }
0120
0121   }
0122
0123 //if we are scaling, apply the scale factor to each correction
0124 if(dcc->m_use_scalefactor)
0125   {
0126     dphi *= dcc->m_scalefactor;
0127     dr *= dcc->m_scalefactor;
0128     dz *= dcc->m_scalefactor;
0129   }
0130
0131   phi_new=phi-dphi;
0132   r_new=r-dr;
0133   z_new=z-dz;
0134
0135   // update cluster
0136   const auto x_new = r_new * std::cos(phi_new);
0137   const auto y_new = r_new * std::sin(phi_new);
0138
0139   return {x_new, y_new, z_new};
0140 }