Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2025-01-18 10:10:41 

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

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!