RKTrackRep.h

Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file RKTrackRep.h

/* Copyright 2008-2009, Technische Universitaet Muenchen,
   Authors: Christian Hoeppner & Sebastian Neubert

   This file is part of GENFIT.

   GENFIT is free software: you can redistribute it and/or modify
   it under the terms of the GNU Lesser General Public License as published
   by the Free Software Foundation, either version 3 of the License, or
   (at your option) any later version.

   GENFIT is distributed in the hope that it will be useful,
   but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
   GNU Lesser General Public License for more details.

   You should have received a copy of the GNU Lesser General Public License
   along with GENFIT.  If not, see <http://www.gnu.org/licenses/>.
*/

/*
 */

#ifndef RKTRACKREP_H
#define RKTRACKREP_H


#include "GFAbsTrackRep.h"
#include "GFDetPlane.h"
#include "GFTrackCand.h"
#include "GFPointPath.h"
#include "RKTools.h"
#include <TMatrixD.h>

class RKTrackRep : public GFAbsTrackRep {

 public:

  // Constructors/Destructors ---------
  RKTrackRep();
  RKTrackRep(const TVector3& pos,
       const TVector3& mom,
       const TVector3& poserr,
       const TVector3& momerr,
       const int& PDGCode);

  RKTrackRep(const TVector3& pos,
       const TVector3& mom,
       const int& PDGCode);

  RKTrackRep(const GFTrackCand* aGFTrackCandPtr, int pdgcode = 0);

  virtual ~RKTrackRep();


  virtual GFAbsTrackRep* clone() const {return new RKTrackRep(*this);}
  virtual GFAbsTrackRep* prototype()const{return new RKTrackRep();}


  double extrapolate(const GFDetPlane&, 
         TMatrixD& statePred,
         TMatrixD& covPred);

  double extrapolate(const GFDetPlane&, 
         TMatrixD& statePred);

  void extrapolateToPoint(const TVector3& pos,
         TVector3& poca,
         TVector3& dirInPoca);

  void extrapolateToLine(const TVector3& point1,
         const TVector3& point2,
         TVector3& poca,
         TVector3& dirInPoca,
         TVector3& poca_onwire);


  double stepalong(double h,
         TVector3& point,
         TVector3& dir);


  TVector3 getPos(const GFDetPlane&);

  TVector3 getMom(const GFDetPlane&);

  void getPosMom(const GFDetPlane&,TVector3& pos,TVector3& mom);

  void getPosMomCov(const GFDetPlane& pl,
                    TVector3& pos, TVector3& mom,
                    TMatrixD& cov);

  double getCharge()const {return fCharge;}

  int getPDG() {return fPdg;};

  void setPropDir(int dir);

  void switchDirection(){fDirection = -1*fDirection;}

  void setPDG(int);


  void setData(const TMatrixD& st,
               const GFDetPlane& pl,
               const TMatrixD* cov=NULL,
               const TMatrixD* aux=NULL);


  void setPosMomCov(const TVector3& pos,
                    const TVector3& mom,
                    const TMatrixD& cov);

  void disableMaterialEffects(bool opt = true){fNoMaterial = opt;}

  const TMatrixD* getAuxInfo(const GFDetPlane& pl);
  
  bool hasAuxInfo() { return true; }



 private:

  void initArrays();

  void calcStateCov(const TVector3& pos,
                    const TVector3& mom,
                    const TVector3& poserr,
                    const TVector3& momerr);

  void calcState(const TVector3& pos,
                 const TVector3& mom);

  void getState7(M1x7& state7);
  void getState7(M1x7& state7, const TMatrixD& state5, const GFDetPlane& pl, const double& spu);
  TMatrixD getState5(const M1x7& state7, const GFDetPlane& pl, double& spu);

  void transformPM7(const TMatrixD& in5x5,
                    M7x7& out7x7,
                    const GFDetPlane& pl,
                    const TMatrixD& state5,
                    const double& spu,
                    TMatrixD* Jac = NULL);

  void transformPM6(const TMatrixD& in5x5,
                    M6x6& out6x6,
                    const GFDetPlane& pl,
                    const TMatrixD& state5,
                    const double& spu,
                    TMatrixD* Jac = NULL);

  void transformM7P(const M7x7& in7x7,
                    TMatrixD& out5x5,
                    const GFDetPlane& pl,
                    const M1x7& state7,
                    TMatrixD* Jac = NULL);

  void transformM6P(const M6x6& in6x6,
                    TMatrixD& out5x5,
                    const GFDetPlane& pl,
                    const M1x7& state7,
                    TMatrixD* Jac = NULL);

  RKTrackRep& operator=(const RKTrackRep*){
    return *this;
  };


  bool RKutta (const GFDetPlane& plane,
               M8x7& P,
               double& coveredDistance,
               std::vector<GFPointPath>& points,
               bool& checkJacProj,
               bool calcCov = true,
               bool onlyOneStep = false,
               double maxStep = 1.E99);

  double estimateStep(std::vector<GFPointPath>& points,
                      const TVector3& pos,
                      const TVector3& dir,
                      const M1x4& SU,
                      const TVector3& MagField,
                      const GFDetPlane& plane,
                      const double& mom,
                      double& relMomLoss,
                      double& deltaAngle,
                      bool& momLossExceeded,
                      bool& atPlane,
                      double maxStep = 1.E99) const;

  TVector3 poca2Line(const TVector3& extr1,
                     const TVector3& extr2,
                     const TVector3& point) const;
    

  double Extrap(const GFDetPlane& plane,
                M1x7& state7,
                M7x7* cov=NULL,
                bool onlyOneStep = false,
                double maxStep = 1.E99);
  
  //RKTrackRep(const RKTrackRep& rhs){};
  
  // data members
  
  int fDirection;   // (-1,0,1) -> (backward prop,decide myself,forward)
  bool fNoMaterial; // don't calculate material effects if true
    
  int fPdg;
  double fMass;
  double fCharge;

  GFDetPlane fCachePlane; 
  double fCacheSpu; 
  double fSpu;
  TMatrixD fAuxInfo;

  // vectors for getState, transform, Extrap etc. functions. Saves a lot of TVector3 constructions/destructions
  TVector3 fO, fU, fV, fW; 
  TVector3 fPos, fDir; 
  TVector3 fpTilde; 
  TVector3 fDirectionBefore, fDirectionAfter; 
  TVector3 fH; 

  // auxiliary variables and arrays
  // needed in Extrap()
  M8x7 fStateJac; 
  M7x7 fNoise; 
  M7x7 fOldCov; 
  // needed in transform...
  M5x7 fJ_pM_5x7; 
  M5x6 fJ_pM_5x6; 
  M7x5 fJ_Mp_7x5; 
  M6x5 fJ_Mp_6x5; 

 public:
  ClassDef(RKTrackRep, 7)

};

#endif