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) const
0061 {
0062   double x;
0063   double y;
0064   double z;
0065   double xshift;
0066   double yshift;
0067   double zshift;
0068   // const double mm = 1.0;
0069   // const double cm = 10.0;
0070   TVector3 shiftposition;
0071
0072   x = position.X();
0073   y = position.Y();
0074   z = position.Z();
0075
0076   double r = position.Perp();
0077   double phi = position.Phi();
0078   if (position.Phi() < 0.0)
0079   {
0080     phi = position.Phi() + 2.0 * M_PI;
0081   }
0082
0083   // distort coordinate of stripe
0084   xshift = 0;
0085   yshift = 0;
0086   zshift = 0;
0087   if (hasTruth)
0088   {
0089     xshift = hX->Interpolate(phi, r, z);
0090     yshift = hY->Interpolate(phi, r, z);
0091     zshift = hZ->Interpolate(phi, r, z);
0092   }
0093
0094   // remove average distortion
0095   if (hasCorrection)
0096   {
0097     double raveshift = hRave->Interpolate(phi, r, z);
0098     double paveshift = hPhiave->Interpolate(phi, r, z);  // hugo confirms the units are cm
0099     double cosphi = cos(phi);
0100     double sinphi = sin(phi);
0101     xshift -= raveshift * cosphi - paveshift * sinphi;
0102     yshift -= raveshift * sinphi + paveshift * cosphi;
0103
0104     zshift -= hZave->Interpolate(phi, r, z);
0105   }
0106
0107   TVector3 forwardshift(x + xshift, y + yshift, z + zshift);
0108
0109   return forwardshift;
0110 }
0111
0112 TVector3 Shifter::ShiftBack(const TVector3 &forwardshift) const
0113 {
0114   double x;
0115   double y;
0116   double z;
0117   // const double mm = 1.0;
0118   // const double cm = 10.0;
0119   TVector3 shiftposition;
0120
0121   x = forwardshift.X();
0122   y = forwardshift.Y();
0123   z = forwardshift.Z();
0124
0125   double rforward = forwardshift.Perp();
0126   double phiforward = forwardshift.Phi();
0127   if (forwardshift.Phi() < 0.0)
0128   {
0129     phiforward += 2.0 * M_PI;
0130   }
0131
0132   double xshiftback = -1 * hXBack->Interpolate(phiforward, rforward, z);
0133   double yshiftback = -1 * hYBack->Interpolate(phiforward, rforward, z);
0134   double zshiftback = -1 * hZBack->Interpolate(phiforward, rforward, z);
0135
0136   shiftposition.SetXYZ(x + xshiftback, y + yshiftback, z + zshiftback);
0137
0138   return shiftposition;
0139 }
0140
0141 TVector3 Shifter::Shift(const TVector3 &position) const
0142 {
0143   return ShiftBack(ShiftForward(position));
0144 }