tpc/fillSpaceChargeMaps/Shifter.cc

0001 #include "Shifter.h"
0002
0003 #include <TFile.h>
0004 #include <TH3.h>
0005 #include <TString.h>
0006 #include <TVector3.h>
0007
0008 #include <cmath>
0009 #include <cstdlib>  // for getenv
0010
0011 Shifter::Shifter(const std::string &truthfilename, const std::string &correctionfilename)
0012 {
0013   // load a 'truth' distortion map and, optionally, a map of a measured correction to those distortions
0014   // this code is currently set up to load a particular correction map that doesn't have distortions
0015   //  in X,Y, and Z components, but rather only in R, R*Phi, and Z components.
0016
0017   // single event distortion file
0018   if (!truthfilename.empty())
0019   {
0020     forward = TFile::Open(truthfilename.c_str(), "READ");
0021     if (forward)
0022     {
0023       forward->GetObject("hIntDistortionX", hX);
0024       forward->GetObject("hIntDistortionY", hY);
0025       forward->GetObject("hIntDistortionZ", hZ);
0026
0027       // not strictly needed, but handy:
0028       forward->GetObject("hIntDistortionR", hR);
0029       forward->GetObject("hIntDistortionP", hPhi);
0030     }
0031   }
0032   if (hX && hY && hZ)
0033   {
0034     hasTruth = true;
0035   }
0036
0037   // single event distortion file
0038   if (!correctionfilename.empty())
0039   {
0040     // average=TFile::Open(correctionfilename,"READ");
0041     //  hardcoded????????
0042     std::string correction_filename = std::string(getenv("CALIBRATIONROOT")) + "/distortion_maps/Distortions_full_realistic_micromegas_all-coarse.root";
0043     average = TFile::Open(correction_filename.c_str(), "READ");
0044     if (average)
0045     {
0046       average->GetObject("hIntDistortionX", hXave);
0047       average->GetObject("hIntDistortionY", hYave);
0048       average->GetObject("hIntDistortionZ", hZave);
0049
0050       average->GetObject("hIntDistortionR", hRave);
0051       average->GetObject("hIntDistortionP", hPhiave);
0052     }
0053   }
0054   if (hXave && hYave && hZave)
0055   {
0056     hasCorrection = true;
0057   }
0058 }
0059
0060 TVector3 Shifter::ShiftForward(const TVector3 &position)
0061 {
0062   double x, y, z, xshift, yshift, zshift;
0063   // const double mm = 1.0;
0064   // const double cm = 10.0;
0065   TVector3 shiftposition;
0066
0067   x = position.X();
0068   y = position.Y();
0069   z = position.Z();
0070
0071   double r = position.Perp();
0072   double phi = position.Phi();
0073   if (position.Phi() < 0.0)
0074   {
0075     phi = position.Phi() + 2.0 * M_PI;
0076   }
0077
0078   // distort coordinate of stripe
0079   xshift = 0;
0080   yshift = 0;
0081   zshift = 0;
0082   if (hasTruth)
0083   {
0084     xshift = hX->Interpolate(phi, r, z);
0085     yshift = hY->Interpolate(phi, r, z);
0086     zshift = hZ->Interpolate(phi, r, z);
0087   }
0088
0089   // remove average distortion
0090   if (hasCorrection)
0091   {
0092     double raveshift = hRave->Interpolate(phi, r, z);
0093     double paveshift = hPhiave->Interpolate(phi, r, z);  // hugo confirms the units are cm
0094     double cosphi = cos(phi);
0095     double sinphi = sin(phi);
0096     xshift -= raveshift * cosphi - paveshift * sinphi;
0097     yshift -= raveshift * sinphi + paveshift * cosphi;
0098
0099     zshift -= hZave->Interpolate(phi, r, z);
0100   }
0101
0102   TVector3 forwardshift(x + xshift, y + yshift, z + zshift);
0103
0104   return forwardshift;
0105 }
0106
0107 TVector3 Shifter::ShiftBack(const TVector3 &forwardshift)
0108 {
0109   double x, y, z;
0110   // const double mm = 1.0;
0111   // const double cm = 10.0;
0112   TVector3 shiftposition;
0113
0114   x = forwardshift.X();
0115   y = forwardshift.Y();
0116   z = forwardshift.Z();
0117
0118   double rforward = forwardshift.Perp();
0119   double phiforward = forwardshift.Phi();
0120   if (forwardshift.Phi() < 0.0)
0121   {
0122     phiforward += 2.0 * M_PI;
0123   }
0124
0125   double xshiftback = -1 * hXBack->Interpolate(phiforward, rforward, z);
0126   double yshiftback = -1 * hYBack->Interpolate(phiforward, rforward, z);
0127   double zshiftback = -1 * hZBack->Interpolate(phiforward, rforward, z);
0128
0129   shiftposition.SetXYZ(x + xshiftback, y + yshiftback, z + zshiftback);
0130
0131   return shiftposition;
0132 }
0133
0134 TVector3 Shifter::Shift(const TVector3 &position)
0135 {
0136   return ShiftBack(ShiftForward(position));
0137 }