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=0,
                        const Int_t refIndex = -1,

                        PndFtsHoughSpaceBinning binning=PndFtsHoughSpaceBinning(),

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

                        PndFtsHoughTrackCand *associatedTrackCand=0,

                        PndFtsHoughTrackerTask *trackerTask=0
        );
        ~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