PndFtsHoughSpace.h

Go to the documentation of this file.

#ifndef PndFtsHoughSpace_H
#define PndFtsHoughSpace_H

#include "PndFtsHoughTrackerTask.h"
#include "PndFtsHoughSpacePeak.h"
#include "PndFtsHoughSpaceBinning.h"

#include "TH2.h"
#include <cmath>
#include <vector>
#include "Rtypes.h"     // for Double_t, Int_t, etc
#include "FairLogger.h" // for FairLogger, MESSAGE_ORIGIN

#include "PndTrackCandHit.h"
#include "TClonesArray.h"

// magnetic field
#include "FairField.h"
#include "TVector3.h"

#include "PndFtsHoughTracklet.h"
#include "PndFtsHoughTrackCand.h"

// For error throwing
#include "TString.h"
#include <stdexcept>

// for saving the path through the Hough space
typedef std::vector<Int_t> IdxPath;               // helper -- path for one hit
typedef std::map<Int_t, IdxPath> HitIdxPathMap;   // that is the one I need -- map of hit indices to path for the corresponding hits
typedef std::pair<Int_t, IdxPath> HitIdxPathPair; // helper for inserting one pair of hit index and path in map

// cout for the above types
std::ostream &operator<<(std::ostream &os, const IdxPath &outVector);
std::ostream &operator<<(std::ostream &os, const HitIdxPathMap &outMap);

class PndFtsHoughSpace : public TH2S {
 public:
  // Constructors/Destructors ---------
  PndFtsHoughSpace(const char *name = nullptr, const Int_t refIndex = -1,

                   PndFtsHoughSpaceBinning binning = PndFtsHoughSpaceBinning(),

                   Double_t zRefPos = 0., Double_t interceptZx = 0.,

                   PndFtsHoughTrackCand *associatedTrackCand = nullptr,

                   PndFtsHoughTrackerTask *trackerTask = nullptr);
  ~PndFtsHoughSpace();

  //    TH2S ExportTH2S();

  // General info for peak finders
  //    The peaks contain the following information
  //            peakTheta = theta for peak
  //            peakSecond = Q/pzx for peak (parabola HT)
  //            peakSecond = intercept for peak (line HT) (in z-x- or z-y-plane)
  //
  //            peakThetaHw = half width of peak in theta
  //            peakSecondHw = half width of peak in second value (see above for what it stands for)
  //
  //            actualHeight height of the peak in the histogram (in counts)

  std::vector<PndFtsHoughTracklet> FindAllPeaksScanPathsMergeBins(const UInt_t minHeight);

  std::vector<PndFtsHoughTracklet> FindAllPeaksScanPathsMergeBinsCalculatingPaths(const UInt_t minHeight);

  std::vector<PndFtsHoughTracklet> FindAllPeaksWithTSpectrum2(const UInt_t minHeight);

  std::vector<PndFtsHoughTracklet> FindAllPeaksBinsWoMergingWithSearchWindow(const UInt_t minHeight, const Int_t vicinityLength = 0);

  std::vector<PndFtsHoughTracklet> FindAllPeaksBlanko(const UInt_t minHeight);

  void FillHoughSpace();

  // operators
  // PndFtsHoughSpace are the same if they contain the same hits, that means if the PndTrackCand are the same, therefore no need to implement that operator here

  // Accessors -----------------------
  inline void Print() const;

  void setVerbose(Int_t verbose) { fVerbose = verbose; };

  Double_t getInterceptZx() const { return fInterceptZx; };

  Double_t getZRefPos() const { return fZRefPos; };

  inline UInt_t GetNHits() const { return fHitId.size(); };

 private:
  // for PandaRoot input/output
  PndFtsHoughTrackerTask *fTrackerTask; 

  void throwError(const TString s) const { throw std::runtime_error(s.Data()); };

  Bool_t setParametersForHsOption(); // set parameters according to the kind of Hough transform I want to do
  void filterInputHits(); // copies input hits (based on z coordinate and skewed/non-skewed) from fFtsHitArray (all FTS hits) to fHitId (only the hits that qualify for the specific
                          // Hough transform)
  inline void AddHitToHS(UInt_t hitId, Double_t rho);
  inline void AddHitToHS(FairLink link, Double_t rho);
  inline Bool_t IsHitFromTubeIdAlreadyAdded(const Int_t tubeIdToAdd);
  inline const PndFtsHit *getHitFromHS(UInt_t index) const;                                // gets the FTS hit corresponding to index
  inline Int_t getHitIdFromHS(UInt_t index) const { return fHitId.at(index).GetHitId(); }; // gets the FTS hit Id corresponding to index
  inline const TVector3 CalculateHitPosFromIntersectionsWithZxTrackModel(const PndFtsHit *const myHit) const;
  inline const TVector3 GetRawOrCalculatedHitPos(const PndFtsHit *const myHit) const;

  void AddHitsToTrackletByCalculating(PndFtsHoughTracklet *currentTracklet, Int_t locmax);

  HitIdxPathMap fHitThetaYIdxPath; 

  // Private Data Members ------------
  Int_t fVerbose;        
  const Int_t fRefIndex; 

  Int_t fFtsBranchId;
  std::vector<PndTrackCandHit> fHitId; 

  PndFtsHoughTrackCand *fAssociatedTrackCand; 

  // only hits with a z value (in the laboratory system, zreal) between onlyusehitsfromz and onlyusehitsuptoz will be used for building the houghspace
  Double_t fOnlyUseHitsFromZ;
  Double_t fOnlyUseHitsUpToZ;
  Bool_t fUseNonSkewedStraws; // if kTRUE, then hits from non-skewed straws are used for Hough transform
  Bool_t fUseSkewedStraws;    // if kTRUE, then hits from skewed straws are used for Hough transform

  Bool_t fKeepBConstant; // set only to kFALSE for testing
  // If kTRUE the y-component of the B-field is not used in the parabola hough transform
  // if kFALSE the parabola's shape will be adjusted based on the magnetic field maps

  // at which z value the values are calculated
  Double_t fZRefPos; // is used to redefine an origin for the coordinate system (so that the angle definition gives meaningful theta values)
  //            Double_t interceptZx; // cannot be constant, because might need to be reset if set incorrectly (has to be 0 for line HT)
  // is used to shift the true x values of hits so that they hit the point (zOffset|0) in z-x-plane (value is determined by line fit on chambers1+2)
  // (zreal=zOffset, xreal=interceptZx) = (zshifted = 0, xshifted = 0)
  // zshifted = zreal - zRefPos
  // xshifted = xreal - interceptZx
  // zreal = zshifted + zRefPos
  // xreal = xshifted + interceptZx
  // For the z-x-plane parabola, a shift in x (hitshiftinx) needs to be set (which should be the result of the straight line hough transform before the dipole field)
  // For the straight line (stations before dipole field) hitshiftinx HAS TO BE ZERO

  // is used only for parabola in zx plane to shift the true x values of hits so that they hit the point (zOffset|0) in z-x-plane (value is determined by line fit on chambers1+2)
  Double_t fInterceptZx;

  // for B field access
  FairField *fField; 

  // for HoughTransform
  // if kTRUE will correct the pzx prediction according to values which should be obtained from a line fit mc truth momentum VS. reco momentum with high statistics
  static const Bool_t fCorrectPzx = kFALSE;

  inline Bool_t FillHoles(const Int_t lastBinX, const Int_t lastBinY, const Int_t currentBinY, IdxPath *ptrThetaYIdxPathVec);

  inline Double_t equationParabola(Double_t thetaRad, Double_t hitZShifted, Double_t hitXShifted, Double_t By)
  {
    // for parabola equation
    const Double_t n = 1.;
    const Double_t e = 1.;
    const Double_t c = n * e / 2.;

    Double_t yVal = 1. / c / By * (-hitZShifted * sin(thetaRad) + hitXShifted * cos(thetaRad)) / pow((hitZShifted * cos(thetaRad) + hitXShifted * sin(thetaRad)), 2);

    // next line is with rotation as in paper (I believe it is incorrect)
    //                  value = 1. / c / By     * (hitZshifted * sin(realtheta) - hitXshifted * cos(realtheta))/ pow((hitZshifted * cos(realtheta) + hitXshifted * sin(realtheta)), 2);
    return yVal;
  };

  inline Double_t equationParabolaPz(Double_t thetaRad, Double_t hitZShifted, Double_t hitXShifted, Double_t By)
  {
    // for parabola equation
    const Double_t n = 1.;
    const Double_t e = 1.;
    const Double_t c = n * e / 2.;

    // Use shifted x and shifted z for parabola
    Double_t yVal = c * By * pow((hitZShifted * cos(thetaRad) + hitXShifted * sin(thetaRad)), 2) / (-hitZShifted * sin(thetaRad) + hitXShifted * cos(thetaRad));

    // next line is with rotation as in paper (I believe it is incorrect)
    //                  value = c*By*pow((hitZshifted * cos(realtheta) + hitXshifted * sin(realtheta)), 2)/(hitZshifted * sin(realtheta) - hitXshifted * cos(realtheta));

    return yVal;
  };

  inline Double_t equationLineZxOrZy(Double_t thetaRad, Double_t hitZShifted, Double_t hitXLabSys)
  { // TODO: Check if that also works in zy plane
    // calculate b which is the distance of point on line at z = zOffset from z axis
    Double_t yVal = -tan(thetaRad) * hitZShifted + hitXLabSys;
    return yVal;
  }

  //------
  // DEBUG
  //-----
  inline TString GetDebugOutName(TString title = "", Int_t param = -1) const;
  void WriteHistoOfHoughSpace() const;
  TH2S MakeEmptyHistoOfSameDimensions(TString specifier = "", Int_t index = -1) const;
  void WriteHistoOfAllPeaks(const std::vector<PndFtsHoughSpacePeak> &peaksToPlot) const;
  void WriteHistoOfAllPaths() const;
  void WriteHistoOfAllPathsForEachMcTruthTrack() const;
  inline void PrintFoundTracklets(const std::vector<PndFtsHoughTracklet> &tracklets) const;
  inline Double_t getByFromBField(Double_t hitXLabSys, Double_t hitYLabSys, Double_t hitZLabSys);

 public:
  ClassDef(PndFtsHoughSpace, 1);
};

// inline functions
void PndFtsHoughSpace::Print() const
{
  std::cout << "=========== PndFtsHoughSpace::Print() ==========" << '\n';
  std::cout << "fZRefPos = " << fZRefPos << '\n';
  std::cout << "fInterceptZx = " << fInterceptZx << 'n';
  TH2S::Print();
}

void PndFtsHoughSpace::PrintFoundTracklets(const std::vector<PndFtsHoughTracklet> &tracklets) const
{
  std::cout << tracklets.size() << " peaks found for " << GetName() << '\n';
  if (10 < fVerbose) {
    for (UInt_t i = 0; i < tracklets.size(); ++i) {
      tracklets[i].Print();
    }
  }
}

const PndFtsHit *PndFtsHoughSpace::getHitFromHS(UInt_t index) const
{
  if (GetNHits() < index)
    throwError("index too high");

  Int_t hitIndex = fHitId.at(index).GetHitId();
  const PndFtsHit *const myHit = (PndFtsHit *)fTrackerTask->GetFtsHit(hitIndex);
  return myHit;
}

const TVector3 PndFtsHoughSpace::GetRawOrCalculatedHitPos(const PndFtsHit *const myHit) const
{
  // get hit position
  TVector3 hitPos;
  if (kFALSE == myHit->GetSkewed()) {
    myHit->Position(hitPos);
  } else {
    hitPos = CalculateHitPosFromIntersectionsWithZxTrackModel(myHit);
  }
  return hitPos;
}

const TVector3 PndFtsHoughSpace::CalculateHitPosFromIntersectionsWithZxTrackModel(const PndFtsHit *const myHit) const
{
  if (0 == fAssociatedTrackCand)
    throwError("No track cand associated with Hough space. Cannot calculate possible hit pos. from track model.");

  const PndFtsTube *ftsTube = fTrackerTask->GetFtsTube(myHit);
  const TVector3 wireDirection = ftsTube->GetWireDirection();
  const TVector3 wireCenter = ftsTube->GetPosition();

  const Double_t hitZLabSys = wireCenter.Z();

  // calculate xTM according to track model at hitZLabSys
  const Double_t xTM = fAssociatedTrackCand->getXLabSys(hitZLabSys);
  // calculate param which is needed for xStraw = xTM
  // xTM = wireCenter.X() + param*wireDirection.X()
  const Double_t param = (xTM - wireCenter.X()) / wireDirection.X();
  // calculate corresponding yStraw
  const Double_t yStraw = wireCenter.Y() + param * wireDirection.Y();

  TVector3 crossingPosition(xTM, yStraw, hitZLabSys);
  return crossingPosition;
}

TString PndFtsHoughSpace::GetDebugOutName(TString title, Int_t param) const
{
  TString debugOut = "/home/plots/";
  debugOut += fTrackerTask->GetEventNr();
  debugOut += title;
  if (-1 != fRefIndex) {
    debugOut += "-rId";
    debugOut += fRefIndex;
  }
  if (-1 != param) {
    debugOut += param;
  }
  debugOut += ".rtg"; // root textual graphics ;) -- actually just a macro // png does not work in this way
  return debugOut;
}

Double_t PndFtsHoughSpace::getByFromBField(Double_t hitXLabSys, Double_t hitYLabSys, Double_t hitZLabSys)
{
  Double_t By = 0.;
  // set y component of B field constant
  if (kTRUE == fKeepBConstant) {
    // do not take B field into account
    By = 1.;
  } else {
    // Use B field information
    Double_t po[3], BB[3];
    po[0] = hitXLabSys; // Use magnetic field at real (not shifted) x position
    po[1] = hitYLabSys;
    po[2] = hitZLabSys;
    fField->GetFieldValue(po, BB); // return value in KG (G3)
    By = BB[1] / 10.;              // By is y-component of magnetic field in Tesla
  }
  return By;
}

#endif