include/root/TEveTrackPropagator.h

0001 // @(#)root/eve:$Id$
0002 // Authors: Matevz Tadel & Alja Mrak-Tadel: 2006, 2007
0003
0004 /*************************************************************************
0005  * Copyright (C) 1995-2007, Rene Brun and Fons Rademakers.               *
0006  * All rights reserved.                                                  *
0007  *                                                                       *
0008  * For the licensing terms see $ROOTSYS/LICENSE.                         *
0009  * For the list of contributors see $ROOTSYS/README/CREDITS.             *
0010  *************************************************************************/
0011
0012 #ifndef ROOT_TEveTrackPropagator
0013 #define ROOT_TEveTrackPropagator
0014
0015 #include "TEveVector.h"
0016 #include "TEvePathMark.h"
0017 #include "TEveUtil.h"
0018 #include "TEveElement.h"
0019 #include "TMarker.h"
0020
0021 #include <vector>
0022
0023 class TEvePointSet;
0024
0025
0026 //==============================================================================
0027 // TEveMagField
0028 //==============================================================================
0029
0030 class TEveMagField
0031 {
0032 protected:
0033    Bool_t  fFieldConstant;
0034
0035 public:
0036    TEveMagField() : fFieldConstant(kFALSE) {}
0037    virtual ~TEveMagField() {}
0038
0039    virtual Bool_t IsConst() const { return fFieldConstant; }
0040
0041    virtual void  PrintField(Double_t x, Double_t y, Double_t z) const
0042    {
0043       TEveVector b = GetField(x, y, z);
0044       printf("v(%f, %f, %f) B(%f, %f, %f) \n", x, y, z, b.fX, b.fY, b.fZ);
0045    }
0046
0047    // Track propagator uses only GetFieldD() and GetMaxFieldMagD().
0048    //
0049    // Here, in this base-class, we have to keep float versions GetField() and GetMaxFieldMag() as
0050    // primary sources of field data (for backward compatibility in ALICE and CMS).
0051    //
0052    // One only needs to redefine the double versions in subclasses.
0053
0054    virtual Double_t    GetMaxFieldMagD() const { return GetMaxFieldMag(); }
0055    virtual TEveVectorD GetFieldD(Double_t x, Double_t y, Double_t z) const { return GetField(x, y, z); }
0056
0057    TEveVectorD GetFieldD(const TEveVectorD &v) const { return GetFieldD(v.fX, v.fY, v.fZ); }
0058
0059    virtual Float_t    GetMaxFieldMag() const { return 4; }
0060    virtual TEveVector GetField(Float_t, Float_t, Float_t) const { return TEveVector(); }
0061
0062    ClassDef(TEveMagField, 0); // Abstract interface to magnetic field
0063 };
0064
0065
0066 //==============================================================================
0067 // TEveMagFieldConst
0068 //==============================================================================
0069
0070 class TEveMagFieldConst : public TEveMagField
0071 {
0072 protected:
0073    TEveVectorD fB;
0074
0075 public:
0076    TEveMagFieldConst(Double_t x, Double_t y, Double_t z) :
0077       TEveMagField(), fB(x, y, z)
0078    { fFieldConstant = kTRUE; }
0079    ~TEveMagFieldConst() override {}
0080
0081    Double_t    GetMaxFieldMagD() const override { return fB.Mag(); };
0082    TEveVectorD GetFieldD(Double_t /*x*/, Double_t /*y*/, Double_t /*z*/) const override { return fB; }
0083
0084    ClassDefOverride(TEveMagFieldConst, 0); // Interface to constant magnetic field.
0085 };
0086
0087
0088 //==============================================================================
0089 // TEveMagFieldDuo
0090 //==============================================================================
0091
0092 class TEveMagFieldDuo : public TEveMagField
0093 {
0094 protected:
0095    TEveVectorD fBIn;
0096    TEveVectorD fBOut;
0097    Double_t    fR2;
0098
0099 public:
0100    TEveMagFieldDuo(Double_t r, Double_t bIn, Double_t bOut) :
0101       TEveMagField(),
0102       fBIn(0,0,bIn), fBOut(0,0,bOut), fR2(r*r)
0103    {
0104       fFieldConstant = kFALSE;
0105    }
0106    ~TEveMagFieldDuo() override {}
0107
0108    Double_t GetMaxFieldMagD() const override { return std::max(fBIn.Mag(), fBOut.Mag()); }
0109
0110    TEveVectorD GetFieldD(Double_t x, Double_t y, Double_t /*z*/) const override
0111    { return  ((x*x+y*y)<fR2) ? fBIn : fBOut; }
0112
0113    ClassDefOverride(TEveMagFieldDuo, 0); // Interface to magnetic field with two different values depending on radius.
0114 };
0115
0116
0117 //==============================================================================
0118 // TEveTrackPropagator
0119 //==============================================================================
0120
0121 class TEveTrackPropagator : public TEveElementList,
0122                             public TEveRefBackPtr
0123 {
0124    friend class TEveTrackPropagatorSubEditor;
0125
0126 public:
0127    enum EStepper_e    { kHelix, kRungeKutta };
0128
0129    enum EProjTrackBreaking_e { kPTB_Break, kPTB_UseFirstPointPos, kPTB_UseLastPointPos };
0130
0131 protected:
0132    struct Helix_t
0133    {
0134       Int_t    fCharge;   // Charge of tracked particle.
0135       Double_t fMaxAng;   // Maximum step angle.
0136       Double_t fMaxStep;  // Maximum allowed step size.
0137       Double_t fDelta;    // Maximum error in the middle of the step.
0138
0139       Double_t fPhi;      // Accumulated angle to check fMaxOrbs by propagator.
0140       Bool_t   fValid;    // Corner case pT~0 or B~0, possible in variable mag field.
0141
0142       // ----------------------------------------------------------------
0143
0144       // helix parameters
0145       Double_t fLam;         // Momentum ratio pT/pZ.
0146       Double_t fR;           // Helix radius in cm.
0147       Double_t fPhiStep;     // Caluclated from fMinAng and fDelta.
0148       Double_t fSin, fCos;   // Current sin/cos(phistep).
0149
0150       // Runge-Kutta parameters
0151       Double_t fRKStep;      // Step for Runge-Kutta.
0152
0153       // cached
0154       TEveVectorD fB;        // Current magnetic field, cached.
0155       TEveVectorD fE1, fE2, fE3; // Base vectors: E1 -> B dir, E2->pT dir, E3 = E1xE2.
0156       TEveVectorD fPt, fPl;  // Transverse and longitudinal momentum.
0157       Double_t fPtMag;       // Magnitude of pT.
0158       Double_t fPlMag;       // Momentum parallel to mag field.
0159       Double_t fLStep;       // Transverse step arc-length in cm.
0160
0161       // ----------------------------------------------------------------
0162
0163       Helix_t();
0164
0165       void UpdateCommon(const TEveVectorD & p, const TEveVectorD& b);
0166       void UpdateHelix (const TEveVectorD & p, const TEveVectorD& b, Bool_t full_update, Bool_t enforce_max_step);
0167       void UpdateRK    (const TEveVectorD & p, const TEveVectorD& b);
0168
0169       void Step(const TEveVector4D& v, const TEveVectorD& p, TEveVector4D& vOut, TEveVectorD& pOut);
0170
0171       Double_t GetStep()  { return fLStep * TMath::Sqrt(1 + fLam*fLam); }
0172       Double_t GetStep2() { return fLStep * fLStep * (1 + fLam*fLam);   }
0173    };
0174
0175
0176 private:
0177    TEveTrackPropagator(const TEveTrackPropagator&);            // Not implemented
0178    TEveTrackPropagator& operator=(const TEveTrackPropagator&); // Not implemented
0179
0180    void DistributeOffset(const TEveVectorD& off, Int_t first_point, Int_t np, TEveVectorD& p);
0181
0182 protected:
0183    EStepper_e               fStepper;
0184
0185    TEveMagField*            fMagFieldObj;
0186    Bool_t                   fOwnMagFiledObj;
0187
0188    // Track extrapolation limits
0189    Double_t                  fMaxR;          // Max radius for track extrapolation
0190    Double_t                  fMaxZ;          // Max z-coordinate for track extrapolation.
0191    Int_t                     fNMax;          // Max steps
0192    // Helix limits
0193    Double_t                  fMaxOrbs;       // Maximal angular path of tracks' orbits (1 ~ 2Pi).
0194
0195    // Path-mark / first-vertex control
0196    Bool_t                   fEditPathMarks; // Show widgets for path-mark control in GUI editor.
0197    Bool_t                   fFitDaughters;  // Pass through daughter creation points when extrapolating a track.
0198    Bool_t                   fFitReferences; // Pass through given track-references when extrapolating a track.
0199    Bool_t                   fFitDecay;      // Pass through decay point when extrapolating a track.
0200    Bool_t                   fFitCluster2Ds; // Pass through 2D-clusters when extrapolating a track.
0201    Bool_t                   fFitLineSegments; // Pass through line when extrapolating a track.
0202    Bool_t                   fRnrDaughters;  // Render daughter path-marks.
0203    Bool_t                   fRnrReferences; // Render track-reference path-marks.
0204    Bool_t                   fRnrDecay;      // Render decay path-marks.
0205    Bool_t                   fRnrCluster2Ds; // Render 2D-clusters.
0206    Bool_t                   fRnrFV;         // Render first vertex.
0207    TMarker                  fPMAtt;         // Marker attributes for rendering of path-marks.
0208    TMarker                  fFVAtt;         // Marker attributes for fits vertex.
0209
0210    // Handling of discontinuities in projections
0211    UChar_t                  fProjTrackBreaking; // Handling of projected-track breaking.
0212    Bool_t                   fRnrPTBMarkers;     // Render break-points on tracks.
0213    TMarker                  fPTBAtt;            // Marker attributes for track break-points.
0214
0215    // ----------------------------------------------------------------
0216
0217    // Propagation, state of current track
0218    std::vector<TEveVector4D> fPoints;        // Calculated point.
0219    std::vector<TEveVector4D> fLastPoints;    // Copy of the latest calculated points.
0220    TEveVectorD               fV;             // Start vertex.
0221    Helix_t                   fH;             // Helix.
0222
0223    void    RebuildTracks();
0224    void    Update(const TEveVector4D& v, const TEveVectorD& p, Bool_t full_update=kFALSE, Bool_t enforce_max_step=kFALSE);
0225    void    Step(const TEveVector4D &v, const TEveVectorD &p, TEveVector4D &vOut, TEveVectorD &pOut);
0226
0227    Bool_t  LoopToVertex(TEveVectorD& v, TEveVectorD& p);
0228    Bool_t  LoopToLineSegment(const TEveVectorD& s, const TEveVectorD& r, TEveVectorD& p);
0229    void    LoopToBounds(TEveVectorD& p);
0230
0231    Bool_t  LineToVertex (TEveVectorD& v);
0232    void    LineToBounds (TEveVectorD& p);
0233
0234    void    StepRungeKutta(Double_t step, Double_t* vect, Double_t* vout);
0235
0236    Bool_t  HelixIntersectPlane(const TEveVectorD& p, const TEveVectorD& point, const TEveVectorD& normal,
0237                                TEveVectorD&itsect);
0238    Bool_t  LineIntersectPlane(const TEveVectorD& p, const TEveVectorD& point, const TEveVectorD& normal,
0239                               TEveVectorD& itsect);
0240    Bool_t  PointOverVertex(const TEveVector4D& v0, const TEveVector4D& v, Double_t* p = nullptr);
0241
0242    void    ClosestPointFromVertexToLineSegment(const TEveVectorD& v, const TEveVectorD& s, const TEveVectorD& r, Double_t rMagInv, TEveVectorD& c);
0243    Bool_t  ClosestPointBetweenLines(const TEveVectorD&, const TEveVectorD&, const TEveVectorD&, const TEveVectorD&, TEveVectorD& out);
0244
0245 public:
0246    TEveTrackPropagator(const char* n="TEveTrackPropagator", const char* t="",
0247                        TEveMagField* field=nullptr, Bool_t own_field=kTRUE);
0248    ~TEveTrackPropagator() override;
0249
0250    void OnZeroRefCount() override;
0251
0252    void CheckReferenceCount(const TEveException& eh="TEveElement::CheckReferenceCount ") override;
0253
0254    void ElementChanged(Bool_t update_scenes=kTRUE, Bool_t redraw=kFALSE) override;
0255
0256    // propagation
0257    void   InitTrack(const TEveVectorD& v, Int_t charge);
0258    void   ResetTrack();
0259
0260    Int_t    GetCurrentPoint() const;
0261    Double_t GetTrackLength(Int_t start_point=0, Int_t end_point=-1) const;
0262
0263    virtual void   GoToBounds(TEveVectorD& p);
0264    virtual Bool_t GoToVertex(TEveVectorD& v, TEveVectorD& p);
0265    virtual Bool_t GoToLineSegment(const TEveVectorD& s, const TEveVectorD& r, TEveVectorD& p);
0266
0267    // TEveVectorF wrappers
0268    void   InitTrack(const TEveVectorF& v, Int_t charge);
0269    void   GoToBounds(TEveVectorF& p);
0270    Bool_t GoToVertex(TEveVectorF& v, TEveVectorF&p);
0271    Bool_t GoToLineSegment(const TEveVectorF& s, const TEveVectorF& r, TEveVectorF& p);
0272
0273    Bool_t IntersectPlane(const TEveVectorD& p, const TEveVectorD& point, const TEveVectorD& normal,
0274                          TEveVectorD& itsect);
0275
0276    void   FillPointSet(TEvePointSet* ps) const;
0277
0278    void   SetStepper(EStepper_e s) { fStepper = s; }
0279
0280    void   SetMagField(Double_t bX, Double_t bY, Double_t bZ);
0281    void   SetMagField(Double_t b) { SetMagField(0, 0, b); }
0282    void   SetMagFieldObj(TEveMagField* field, Bool_t own_field=kTRUE);
0283
0284    void   SetMaxR(Double_t x);
0285    void   SetMaxZ(Double_t x);
0286    void   SetMaxOrbs(Double_t x);
0287    void   SetMinAng(Double_t x);
0288    void   SetMaxAng(Double_t x);
0289    void   SetMaxStep(Double_t x);
0290    void   SetDelta(Double_t x);
0291
0292    void   SetEditPathMarks(Bool_t x) { fEditPathMarks = x; }
0293    void   SetRnrDaughters(Bool_t x);
0294    void   SetRnrReferences(Bool_t x);
0295    void   SetRnrDecay(Bool_t x);
0296    void   SetRnrCluster2Ds(Bool_t x);
0297    void   SetFitDaughters(Bool_t x);
0298    void   SetFitReferences(Bool_t x);
0299    void   SetFitDecay(Bool_t x);
0300    void   SetFitCluster2Ds(Bool_t x);
0301    void   SetFitLineSegments(Bool_t x);
0302    void   SetRnrFV(Bool_t x);
0303    void   SetProjTrackBreaking(UChar_t x);
0304    void   SetRnrPTBMarkers(Bool_t x);
0305
0306    TEveVectorD GetMagField(Double_t x, Double_t y, Double_t z) { return fMagFieldObj->GetField(x, y, z); }
0307    void PrintMagField(Double_t x, Double_t y, Double_t z) const;
0308
0309    EStepper_e   GetStepper()  const { return fStepper;}
0310
0311    Double_t GetMaxR()     const { return fMaxR;     }
0312    Double_t GetMaxZ()     const { return fMaxZ;     }
0313    Double_t GetMaxOrbs()  const { return fMaxOrbs;  }
0314    Double_t GetMinAng()   const;
0315    Double_t GetMaxAng()   const { return fH.fMaxAng;   }
0316    Double_t GetMaxStep()  const { return fH.fMaxStep;  }
0317    Double_t GetDelta()    const { return fH.fDelta;    }
0318
0319    Bool_t  GetEditPathMarks() const { return fEditPathMarks; }
0320    Bool_t  GetRnrDaughters()  const { return fRnrDaughters;  }
0321    Bool_t  GetRnrReferences() const { return fRnrReferences; }
0322    Bool_t  GetRnrDecay()      const { return fRnrDecay;      }
0323    Bool_t  GetRnrCluster2Ds() const { return fRnrCluster2Ds; }
0324    Bool_t  GetFitDaughters()  const { return fFitDaughters;  }
0325    Bool_t  GetFitReferences() const { return fFitReferences; }
0326    Bool_t  GetFitDecay()      const { return fFitDecay;      }
0327    Bool_t  GetFitCluster2Ds() const { return fFitCluster2Ds; }
0328    Bool_t  GetFitLineSegments() const { return fFitLineSegments; }
0329    Bool_t  GetRnrFV()         const { return fRnrFV;         }
0330    UChar_t GetProjTrackBreaking() const { return fProjTrackBreaking; }
0331    Bool_t  GetRnrPTBMarkers()     const { return fRnrPTBMarkers; }
0332
0333    TMarker& RefPMAtt()  { return fPMAtt; }
0334    TMarker& RefFVAtt()  { return fFVAtt; }
0335    TMarker& RefPTBAtt() { return fPTBAtt; }
0336
0337    const std::vector<TEveVector4D>& GetLastPoints() const { return fLastPoints; }
0338
0339    static Bool_t IsOutsideBounds(const TEveVectorD& point, Double_t maxRsqr, Double_t maxZ);
0340
0341    static Double_t             fgDefMagField; // Default value for constant solenoid magnetic field.
0342    static const Double_t       fgkB2C;        // Constant for conversion of momentum to curvature.
0343    static TEveTrackPropagator  fgDefault;     // Default track propagator.
0344
0345    static Double_t             fgEditorMaxR;  // Max R that can be set in GUI editor.
0346    static Double_t             fgEditorMaxZ;  // Max Z that can be set in GUI editor.
0347
0348    ClassDefOverride(TEveTrackPropagator, 0); // Calculates path of a particle taking into account special path-marks and imposed boundaries.
0349 };
0350
0351 //______________________________________________________________________________
0352 inline Bool_t TEveTrackPropagator::IsOutsideBounds(const TEveVectorD& point,
0353                                                    Double_t           maxRsqr,
0354                                                    Double_t           maxZ)
0355 {
0356    // Return true if point% is outside of cylindrical bounds detrmined by
0357    // square radius and z.
0358
0359    return TMath::Abs(point.fZ) > maxZ ||
0360           point.fX*point.fX + point.fY*point.fY > maxRsqr;
0361 }
0362
0363 //______________________________________________________________________________
0364 inline Bool_t TEveTrackPropagator::PointOverVertex(const TEveVector4D &v0,
0365                                                    const TEveVector4D &v,
0366                                                    Double_t           *p)
0367 {
0368    static const Double_t kMinPl = 1e-5;
0369
0370    TEveVectorD dv; dv.Sub(v0, v);
0371
0372    Double_t dotV;
0373
0374    if (TMath::Abs(fH.fPlMag) > kMinPl)
0375    {
0376       // Use longitudinal momentum to determine crossing point.
0377       // Works ok for spiraling helices, also for loopers.
0378
0379       dotV = fH.fE1.Dot(dv);
0380       if (fH.fPlMag < 0)
0381          dotV = -dotV;
0382    }
0383    else
0384    {
0385       // Use full momentum, which is pT, under this conditions.
0386
0387       dotV = fH.fE2.Dot(dv);
0388    }
0389
0390    if (p)
0391       *p = dotV;
0392
0393    return dotV < 0;
0394 }
0395
0396 #endif