g4simulation/g4tpc/PHG4TpcDistortion.cc

0001 // $Id: $
0002
0003 /*!
0004  * \file PHG4TpcDistortion.cc
0005  * \brief
0006  * \author Jin Huang <jhuang@bnl.gov>, Henry Klest <henry.klest@stonybrook.edu>, Hugo Pereira Da Costa <hugo.pereira-da-costa@cea.fr>
0007  * \version $Revision:   $
0008  * \date $Date: $
0009  */
0010
0011 #include "PHG4TpcDistortion.h"
0012
0013 #include <TFile.h>
0014 #include <TH3.h>
0015 #include <TTree.h>
0016
0017 #include <cmath>    // for sqrt, fabs, NAN
0018 #include <cstdlib>  // for exit
0019 #include <iostream>
0020
0021 namespace
0022 {
0023   template <class T>
0024   constexpr T square(const T& x)
0025   {
0026     return x * x;
0027   }
0028
0029   // check boundaries in axis
0030   /* for the interpolation to work, the value must be within the range of the provided axis, and not into the first and last bin */
0031   inline bool check_boundaries(const TAxis* axis, double value)
0032   {
0033     const auto bin = axis->FindBin(value);
0034     return (bin >= 2 && bin < axis->GetNbins());
0035   }
0036
0037   // check boundaries in histogram, before interpolation
0038   /* for the interpolation to work, the value must be within the range of the provided axis, and not into the first and last bin */
0039   inline bool check_boundaries(const TH3* h, double r, double phi, double z)
0040   {
0041     return check_boundaries(h->GetXaxis(), r) && check_boundaries(h->GetYaxis(), phi) && check_boundaries(h->GetZaxis(), z);
0042   }
0043
0044   // print histogram
0045   [[maybe_unused]] void print_histogram(TH3* h)
0046   {
0047     std::cout << "PHG4TpcDistortion::print_histogram - name: " << h->GetName() << std::endl;
0048     for (const auto& axis : {h->GetXaxis(), h->GetYaxis(), h->GetZaxis()})
0049     {
0050       std::cout
0051           << "  " << axis->GetName()
0052           << " bins: " << axis->GetNbins()
0053           << " min: " << axis->GetXmin()
0054           << " max: " << axis->GetXmax()
0055           << std::endl;
0056     }
0057     std::cout << std::endl;
0058   }
0059
0060 }  // namespace
0061
0062 //__________________________________________________________________________________________________________
0063 void PHG4TpcDistortion::Init()
0064 {
0065   std::cout << "PHG4TpcDistortion::Init - m_phi_hist_in_radians: " << m_phi_hist_in_radians << std::endl;
0066
0067   if (m_do_static_distortions)
0068   {
0069     std::cout << "PHG4TpcDistortion::Init - m_static_distortion_filename: " << m_static_distortion_filename << std::endl;
0070     m_static_tfile.reset(new TFile(m_static_distortion_filename.c_str()));
0071     if (!m_static_tfile->IsOpen())
0072     {
0073       std::cout << "PHG4TpcDistortion::Init - Static distortion file could not be opened!" << std::endl;
0074       exit(1);
0075     }
0076
0077     // Open Static Space Charge Maps
0078     hDRint[0] = dynamic_cast<TH3*>(m_static_tfile->Get("hIntDistortionR_negz"));
0079     hDRint[1] = dynamic_cast<TH3*>(m_static_tfile->Get("hIntDistortionR_posz"));
0080     hDPint[0] = dynamic_cast<TH3*>(m_static_tfile->Get("hIntDistortionP_negz"));
0081     hDPint[1] = dynamic_cast<TH3*>(m_static_tfile->Get("hIntDistortionP_posz"));
0082     hDZint[0] = dynamic_cast<TH3*>(m_static_tfile->Get("hIntDistortionZ_negz"));
0083     hDZint[1] = dynamic_cast<TH3*>(m_static_tfile->Get("hIntDistortionZ_posz"));
0084     if (m_do_ReachesReadout)
0085     {
0086       hReach[0] = dynamic_cast<TH3*>(m_static_tfile->Get("hReachesReadout_negz"));
0087       hReach[1] = dynamic_cast<TH3*>(m_static_tfile->Get("hReachesReadout_posz"));
0088     }
0089   }
0090
0091   if (m_do_time_ordered_distortions)
0092   {
0093     std::cout << "PHG4TpcDistortion::Init - m_time_ordered_distortion_filename: " << m_time_ordered_distortion_filename << std::endl;
0094     m_time_ordered_tfile.reset(new TFile(m_time_ordered_distortion_filename.c_str()));
0095     if (!m_time_ordered_tfile->IsOpen())
0096     {
0097       std::cout << "PHG4TpcDistortion::Init - TimeOrdered distortion file could not be opened!" << std::endl;
0098       exit(1);
0099     }
0100
0101     TimeTree = static_cast<TTree*>(m_time_ordered_tfile->Get("TimeDists"));
0102     if (!TimeTree)
0103     {
0104       std::cout << "PHG4TpcDistortion::Init - TimeOrdered distortion tree could not be found!" << std::endl;
0105       exit(1);
0106     }
0107
0108     // create histograms
0109     TimehDR[0] = new TH3F();
0110     TimehDR[1] = new TH3F();
0111     TimehDP[0] = new TH3F();
0112     TimehDP[1] = new TH3F();
0113     TimehDZ[0] = new TH3F();
0114     TimehDZ[1] = new TH3F();
0115     TimehRR[0] = new TH3F();  // RR stands for ReachesReadout
0116     TimehRR[1] = new TH3F();  // RR stands for ReachesReadout
0117
0118     // assign to tree branches
0119     TimeTree->SetBranchAddress("hIntDistortionR_negz", &(TimehDR[0]));
0120     TimeTree->SetBranchAddress("hIntDistortionR_posz", &(TimehDR[1]));
0121     TimeTree->SetBranchAddress("hIntDistortionP_negz", &(TimehDP[0]));
0122     TimeTree->SetBranchAddress("hIntDistortionP_posz", &(TimehDP[1]));
0123     TimeTree->SetBranchAddress("hIntDistortionZ_negz", &(TimehDZ[0]));
0124     TimeTree->SetBranchAddress("hIntDistortionZ_posz", &(TimehDZ[1]));
0125     if (m_do_ReachesReadout)
0126     {
0127       TimeTree->SetBranchAddress("hReachesReadout_negz", &(TimehRR[0]));
0128       TimeTree->SetBranchAddress("hReachesReadout_posz", &(TimehRR[1]));
0129     }
0130   }
0131 }
0132
0133 //__________________________________________________________________________________________________________
0134 void PHG4TpcDistortion::load_event(int event_num)
0135 {
0136   if (TimeTree)
0137   {
0138     int nentries = TimeTree->GetEntries();
0139     if (event_num > nentries)
0140     {
0141       event_num = event_num % nentries;
0142     }
0143     if (event_num % nentries == 0 && event_num != 0)
0144     {
0145       std::cout << "Distortion map sequence repeating as of event number " << event_num << std::endl;
0146     }
0147     TimeTree->GetEntry(event_num);
0148   }
0149
0150   return;
0151 }
0152
0153 //__________________________________________________________________________________________________________
0154 double PHG4TpcDistortion::get_x_distortion_cartesian(double x, double y, double z) const
0155 {
0156   double r = sqrt(x * x + y * y);
0157   double phi = std::atan2(y, x);
0158
0159   // get components
0160   double dr = get_distortion('r', r, phi, z);
0161   double dphi = get_distortion('p', r, phi, z);
0162
0163   // rotate into cartesian based on local r phi:
0164   double cosphi = cos(phi);
0165   double sinphi = sin(phi);
0166   double dx = dr * cosphi - dphi * sinphi;
0167   return dx;
0168 }
0169
0170 //__________________________________________________________________________________________________________
0171 double PHG4TpcDistortion::get_y_distortion_cartesian(double x, double y, double z) const
0172 {
0173   double r = sqrt(x * x + y * y);
0174   double phi = std::atan2(y, x);
0175
0176   // get components
0177   double dr = get_distortion('r', r, phi, z);
0178   double dphi = get_distortion('p', r, phi, z);
0179
0180   // rotate into cartesian based on local r phi:
0181   double cosphi = cos(phi);
0182   double sinphi = sin(phi);
0183   double dy = dphi * cosphi + dr * sinphi;
0184   return dy;
0185 }
0186
0187 //__________________________________________________________________________________________________________
0188 double PHG4TpcDistortion::get_z_distortion_cartesian(double x, double y, double z) const
0189 {
0190   double r = sqrt(x * x + y * y);
0191   double phi = std::atan2(y, x);
0192
0193   // get components
0194   double dz = get_distortion('z', r, phi, z);
0195
0196   return dz;
0197 }
0198
0199 //__________________________________________________________________________________________________________
0200 double PHG4TpcDistortion::get_r_distortion(double r, double phi, double z) const
0201 {
0202   return get_distortion('r', r, phi, z);
0203 }
0204
0205 //__________________________________________________________________________________________________________
0206 double PHG4TpcDistortion::get_rphi_distortion(double r, double phi, double z) const
0207 {
0208   if (m_phi_hist_in_radians)
0209   {  // if the hist is in radians, multiply by r to get the rphi distortion
0210     return r * get_distortion('p', r, phi, z);
0211   }
0212
0213   return get_distortion('p', r, phi, z);
0214 }
0215
0216 //__________________________________________________________________________________________________________
0217 double PHG4TpcDistortion::get_z_distortion(double r, double phi, double z) const
0218 {
0219   return get_distortion('z', r, phi, z);
0220 }
0221 //__________________________________________________________________________________________________________
0222 double PHG4TpcDistortion::get_reaches_readout(double r, double phi, double z) const
0223 {
0224   if (m_do_ReachesReadout)
0225   {
0226     return get_distortion('R', r, phi, z);
0227   }
0228
0229   return 1;
0230 }
0231
0232 double PHG4TpcDistortion::get_distortion(char axis, double r, double phi, double z) const
0233 {
0234   if (phi < 0)
0235   {
0236     phi += 2 * M_PI;
0237   }
0238   const int zpart = (z > 0 ? 1 : 0);  // z<0 corresponds to the negative side, which is element 0.
0239
0240   TH3* hdistortion = nullptr;
0241
0242   if (axis != 'r' && axis != 'p' && axis != 'z' && axis != 'R')
0243   {
0244     std::cout << "Distortion Requested along axis " << axis << " which is invalid.  Exiting.\n"
0245               << std::endl;
0246     exit(1);
0247   }
0248
0249   double _distortion = 0.;
0250
0251   // select the appropriate histogram:
0252   if (m_do_static_distortions)
0253   {
0254     if (axis == 'r')
0255     {
0256       hdistortion = hDRint[zpart];
0257     }
0258     else if (axis == 'p')
0259     {
0260       hdistortion = hDPint[zpart];
0261     }
0262     else if (axis == 'z')
0263     {
0264       hdistortion = hDZint[zpart];
0265     }
0266     else if (axis == 'R')
0267     {
0268       hdistortion = hReach[zpart];
0269     }
0270     if (hdistortion)
0271     {
0272       if (check_boundaries(hdistortion, phi, r, z))
0273       {
0274         _distortion += hdistortion->Interpolate(phi, r, z);
0275       }
0276     }
0277     else
0278     {
0279       std::cout << "Static Distortion Requested along axis " << axis << ", but distortion map does not exist.  Exiting.\n"
0280                 << std::endl;
0281       exit(1);
0282     }
0283   }
0284
0285   if (m_do_time_ordered_distortions)
0286   {
0287     if (axis == 'r')
0288     {
0289       hdistortion = TimehDR[zpart];
0290     }
0291     else if (axis == 'p')
0292     {
0293       hdistortion = TimehDP[zpart];
0294     }
0295     else if (axis == 'z')
0296     {
0297       hdistortion = TimehDZ[zpart];
0298     }
0299     else if (axis == 'R')
0300     {
0301       hdistortion = TimehRR[zpart];
0302     }
0303     if (hdistortion)
0304     {
0305       if (check_boundaries(hdistortion, phi, r, z))
0306       {
0307         _distortion += hdistortion->Interpolate(phi, r, z);
0308       }
0309     }
0310     else
0311     {
0312       std::cout << "Time Series Distortion Requested along axis " << axis << ", but distortion map does not exist.  Exiting.\n"
0313                 << std::endl;
0314       exit(1);
0315     }
0316   }
0317
0318   return _distortion;
0319 }