ROOT_528-00b_version: include/TEveTrackPropagator.h Source File

00001 // @(#)root/eve:$Id: TEveTrackPropagator.h 33627 2010-05-27 19:19:58Z matevz $
00002 // Authors: Matevz Tadel & Alja Mrak-Tadel: 2006, 2007
00003 
00004 /*************************************************************************
00005  * Copyright (C) 1995-2007, Rene Brun and Fons Rademakers.               *
00006  * All rights reserved.                                                  *
00007  *                                                                       *
00008  * For the licensing terms see $ROOTSYS/LICENSE.                         *
00009  * For the list of contributors see $ROOTSYS/README/CREDITS.             *
00010  *************************************************************************/
00011 
00012 #ifndef ROOT_TEveTrackPropagator
00013 #define ROOT_TEveTrackPropagator
00014 
00015 #include "TEveVector.h"
00016 #include "TEvePathMark.h"
00017 #include "TEveUtil.h"
00018 #include "TEveElement.h"
00019 #include "TMarker.h"
00020 
00021 #include <vector>
00022 
00023 class TEvePointSet;
00024 
00025 
00026 //==============================================================================
00027 // TEveMagField
00028 //==============================================================================
00029 
00030 class TEveMagField
00031 {
00032 protected:
00033    Bool_t  fFieldConstant;
00034 
00035 public:
00036    TEveMagField(): fFieldConstant(kFALSE){}
00037    virtual ~TEveMagField(){}
00038 
00039    virtual Bool_t IsConst() const {return fFieldConstant;};
00040 
00041    virtual void  PrintField(Float_t x, Float_t y, Float_t z) const
00042    {
00043       TEveVector b = GetField(x, y, z);
00044       printf("v(%f, %f, %f) B(%f, %f, %f) \n", x, y, z, b.fX, b.fY, b.fZ);
00045    }
00046 
00047    virtual TEveVector GetField(const TEveVector &v) const { return GetField(v.fX, v.fY, v.fZ);}
00048    virtual TEveVector GetField(Float_t x, Float_t y, Float_t z) const = 0;
00049    virtual Float_t GetMaxFieldMag() const { return 4; } // not abstract because of backward compatibility
00050 
00051    ClassDef(TEveMagField, 0); // Abstract interface to magnetic field
00052 };
00053 
00054 
00055 //==============================================================================
00056 // TEveMagFieldConst
00057 //==============================================================================
00058 
00059 class TEveMagFieldConst : public TEveMagField
00060 {
00061 protected:
00062    TEveVector fB;
00063 
00064 public:
00065    TEveMagFieldConst(Float_t x, Float_t y, Float_t z) : TEveMagField(), fB(x, y, z)
00066    { fFieldConstant = kTRUE; }
00067    virtual ~TEveMagFieldConst() {}
00068 
00069    using   TEveMagField::GetField;
00070    virtual TEveVector GetField(Float_t /*x*/, Float_t /*y*/, Float_t /*z*/) const { return fB; }
00071    virtual Float_t GetMaxFieldMag() const { return fB.Mag(); };
00072 
00073    ClassDef(TEveMagFieldConst, 0); // Interface to constant magnetic field.
00074 };
00075 
00076 
00077 //==============================================================================
00078 // TEveMagFieldDuo
00079 //==============================================================================
00080 
00081 class TEveMagFieldDuo : public TEveMagField
00082 {
00083 protected:
00084    TEveVector fBIn;
00085    TEveVector fBOut;
00086    Float_t    fR2;
00087 
00088 public:
00089    TEveMagFieldDuo(Float_t r, Float_t bIn, Float_t bOut) : TEveMagField(),
00090      fBIn(0,0,bIn), fBOut(0,0,bOut), fR2(r*r)
00091    {
00092       fFieldConstant = kFALSE;
00093    }
00094    virtual ~TEveMagFieldDuo() {}
00095 
00096    using   TEveMagField::GetField;
00097    virtual TEveVector GetField(Float_t x, Float_t y, Float_t /*z*/) const
00098    { return  ((x*x+y*y)<fR2) ? fBIn : fBOut; }
00099    virtual Float_t GetMaxFieldMag() const
00100    { Float_t b1 = fBIn.Mag(), b2 = fBOut.Mag(); return b1 > b2 ? b1 : b2; }
00101 
00102    ClassDef(TEveMagFieldDuo, 0); // Interface to magnetic field with two different values depending of radius.
00103 };
00104 
00105 
00106 //==============================================================================
00107 // TEveTrackPropagator
00108 //==============================================================================
00109 
00110 class TEveTrackPropagator : public TEveElementList,
00111                             public TEveRefBackPtr
00112 {
00113    friend class TEveTrackPropagatorSubEditor;
00114 
00115 public:
00116    struct Helix_t
00117    {
00118       Int_t   fCharge;   // Charge of tracked particle.
00119       Float_t fMaxAng;   // Maximum step angle.
00120       Float_t fMaxStep;  // Maximum allowed step size.
00121       Float_t fDelta;    // Maximum error in the middle of the step.
00122 
00123       Float_t fPhi;      // Accumulated angle to check fMaxOrbs by propagator.
00124       Bool_t  fValid;    // Corner case pT~0 or B~0, possible in variable mag field.
00125 
00126       // ----------------------------------------------------------------
00127 
00128       // helix parameters
00129       Float_t fLam;         // Momentum ratio pT/pZ.
00130       Float_t fR;           // Helix radius in cm.
00131       Float_t fPhiStep;     // Caluclated from fMinAng and fDelta.
00132       Float_t fSin, fCos;   // Current sin/cos(phistep).
00133 
00134       // Runge-Kutta parameters
00135       Float_t fRKStep;      // Step for Runge-Kutta.
00136 
00137       // cached
00138       TEveVector fB;        // Current magnetic field, cached.
00139       TEveVector fE1, fE2, fE3; // Base vectors: E1 -> B dir, E2->pT dir, E3 = E1xE2.
00140       TEveVector fPt, fPl;  // Transverse and longitudinal momentum.
00141       Float_t fPtMag;       // Magnitude of pT.
00142       Float_t fPlMag;       // Momentum parallel to mag field.
00143       Float_t fLStep;       // Transverse step arc-length in cm.
00144 
00145       // ----------------------------------------------------------------
00146 
00147       Helix_t();
00148 
00149       void UpdateCommon(const TEveVector & p, const TEveVector& b);
00150       void UpdateHelix(const TEveVector & p, const TEveVector& b, Bool_t full_update, Bool_t enforce_max_step);
00151       void UpdateRK   (const TEveVector & p, const TEveVector& b);
00152 
00153       void Step(const TEveVector4& v, const TEveVector& p, TEveVector4& vOut, TEveVector& pOut);
00154 
00155       Float_t GetStep()  { return fLStep * TMath::Sqrt(1 + fLam*fLam); }
00156       Float_t GetStep2() { return fLStep * fLStep * (1 + fLam*fLam);   }
00157    };
00158 
00159    enum EStepper_e    { kHelix, kRungeKutta };
00160 
00161    enum EProjTrackBreaking_e { kPTB_Break, kPTB_UseFirstPointPos, kPTB_UseLastPointPos };
00162 
00163 private:
00164    TEveTrackPropagator(const TEveTrackPropagator&);            // Not implemented
00165    TEveTrackPropagator& operator=(const TEveTrackPropagator&); // Not implemented
00166 
00167 protected:
00168    EStepper_e               fStepper;
00169 
00170    TEveMagField*            fMagFieldObj;
00171    Bool_t                   fOwnMagFiledObj;
00172 
00173    // Track extrapolation limits
00174    Float_t                  fMaxR;          // Max radius for track extrapolation
00175    Float_t                  fMaxZ;          // Max z-coordinate for track extrapolation.
00176    Int_t                    fNMax;          // Max steps
00177    // Helix limits
00178    Float_t                  fMaxOrbs;       // Maximal angular path of tracks' orbits (1 ~ 2Pi).
00179 
00180    // Path-mark / first-vertex control
00181    Bool_t                   fEditPathMarks; // Show widgets for path-mark control in GUI editor.
00182    Bool_t                   fFitDaughters;  // Pass through daughter creation points when extrapolating a track.
00183    Bool_t                   fFitReferences; // Pass through given track-references when extrapolating a track.
00184    Bool_t                   fFitDecay;      // Pass through decay point when extrapolating a track.
00185    Bool_t                   fFitCluster2Ds; // Pass through 2D-clusters when extrapolating a track.
00186    Bool_t                   fRnrDaughters;  // Render daughter path-marks.
00187    Bool_t                   fRnrReferences; // Render track-reference path-marks.
00188    Bool_t                   fRnrDecay;      // Render decay path-marks.
00189    Bool_t                   fRnrCluster2Ds; // Render 2D-clusters.
00190    Bool_t                   fRnrFV;         // Render first vertex.
00191    TMarker                  fPMAtt;         // Marker attributes for rendering of path-marks.
00192    TMarker                  fFVAtt;         // Marker attributes for fits vertex.
00193 
00194    // Handling of discontinuities in projections
00195    UChar_t                  fProjTrackBreaking; // Handling of projected-track breaking.
00196    Bool_t                   fRnrPTBMarkers;     // Render break-points on tracks.
00197    TMarker                  fPTBAtt;            // Marker attributes for track break-points.
00198 
00199    // ----------------------------------------------------------------
00200 
00201    // Propagation, state of current track
00202    std::vector<TEveVector4> fPoints;        // Calculated point.
00203    TEveVector               fV;             // Start vertex.
00204    Helix_t                  fH;             // Helix.
00205 
00206    void    RebuildTracks();
00207    void    Update(const TEveVector4& v, const TEveVector& p, Bool_t full_update=kFALSE, Bool_t enforce_max_step=kFALSE);
00208    void    Step(const TEveVector4 &v, const TEveVector &p, TEveVector4 &vOut, TEveVector &pOut);
00209 
00210    Bool_t  LoopToVertex(TEveVector& v, TEveVector& p);
00211    void    LoopToBounds(TEveVector& p);
00212 
00213    Bool_t  LineToVertex (TEveVector& v);
00214    void    LineToBounds (TEveVector& p);
00215 
00216    void    StepRungeKutta(Double_t step, Double_t* vect, Double_t* vout);
00217 
00218    Bool_t  HelixIntersectPlane(const TEveVector& p, const TEveVector& point, const TEveVector& normal,
00219                                TEveVector& itsect);
00220    Bool_t  LineIntersectPlane(const TEveVector& p, const TEveVector& point, const TEveVector& normal,
00221                               TEveVector& itsect);
00222 
00223    Bool_t PointOverVertex(const TEveVector4& v0, const TEveVector4& v, Float_t* p=0);
00224 
00225 public:
00226    TEveTrackPropagator(const char* n="TEveTrackPropagator", const char* t="",
00227                        TEveMagField* field=0, Bool_t own_field=kTRUE);
00228    virtual ~TEveTrackPropagator();
00229 
00230    virtual void OnZeroRefCount();
00231 
00232    virtual void CheckReferenceCount(const TEveException& eh="TEveElement::CheckReferenceCount ");
00233 
00234    virtual void ElementChanged(Bool_t update_scenes=kTRUE, Bool_t redraw=kFALSE);
00235 
00236    // propagation
00237    void   InitTrack(TEveVector& v, Int_t charge);
00238    void   ResetTrack();
00239    void   GoToBounds(TEveVector& p);
00240    Bool_t GoToVertex(TEveVector& v, TEveVector& p);
00241 
00242    Bool_t IntersectPlane(const TEveVector& p, const TEveVector& point, const TEveVector& normal,
00243                          TEveVector& itsect);
00244 
00245    void   FillPointSet(TEvePointSet* ps) const;
00246 
00247    void   SetStepper(EStepper_e s) { fStepper = s; }
00248 
00249    void   SetMagField(Float_t bX, Float_t bY, Float_t bZ);
00250    void   SetMagField(Float_t b) { SetMagField(0.f, 0.f, b); }
00251    void   SetMagFieldObj(TEveMagField* field, Bool_t own_field=kTRUE);
00252 
00253    void   SetMaxR(Float_t x);
00254    void   SetMaxZ(Float_t x);
00255    void   SetMaxOrbs(Float_t x);
00256    void   SetMinAng(Float_t x);
00257    void   SetMaxAng(Float_t x);
00258    void   SetMaxStep(Float_t x);
00259    void   SetDelta(Float_t x);
00260 
00261    void   SetEditPathMarks(Bool_t x) { fEditPathMarks = x; }
00262    void   SetRnrDaughters(Bool_t x);
00263    void   SetRnrReferences(Bool_t x);
00264    void   SetRnrDecay(Bool_t x);
00265    void   SetRnrCluster2Ds(Bool_t x);
00266    void   SetFitDaughters(Bool_t x);
00267    void   SetFitReferences(Bool_t x);
00268    void   SetFitDecay(Bool_t x);
00269    void   SetFitCluster2Ds(Bool_t x);
00270    void   SetRnrFV(Bool_t x);
00271    void   SetProjTrackBreaking(UChar_t x);
00272    void   SetRnrPTBMarkers(Bool_t x);
00273 
00274    TEveVector GetMagField(Float_t x, Float_t y, Float_t z) { return fMagFieldObj->GetField(x, y, z); }
00275    void PrintMagField(Float_t x, Float_t y, Float_t z) const;
00276 
00277    EStepper_e   GetStepper()  const { return fStepper;}
00278 
00279    Float_t GetMaxR()     const { return fMaxR;     }
00280    Float_t GetMaxZ()     const { return fMaxZ;     }
00281    Float_t GetMaxOrbs()  const { return fMaxOrbs;  }
00282    Float_t GetMinAng()   const;
00283    Float_t GetMaxAng()   const { return fH.fMaxAng;   }
00284    Float_t GetMaxStep()  const { return fH.fMaxStep;  }
00285    Float_t GetDelta()    const { return fH.fDelta;    }
00286 
00287    Bool_t  GetEditPathMarks() const { return fEditPathMarks; }
00288    Bool_t  GetRnrDaughters()  const { return fRnrDaughters;  }
00289    Bool_t  GetRnrReferences() const { return fRnrReferences; }
00290    Bool_t  GetRnrDecay()      const { return fRnrDecay;      }
00291    Bool_t  GetRnrCluster2Ds() const { return fRnrCluster2Ds; }
00292    Bool_t  GetFitDaughters()  const { return fFitDaughters;  }
00293    Bool_t  GetFitReferences() const { return fFitReferences; }
00294    Bool_t  GetFitDecay()      const { return fFitDecay;      }
00295    Bool_t  GetFitCluster2Ds() const { return fFitCluster2Ds; }
00296    Bool_t  GetRnrFV()         const { return fRnrFV;         }
00297    UChar_t GetProjTrackBreaking() const { return fProjTrackBreaking; }
00298    Bool_t  GetRnrPTBMarkers()     const { return fRnrPTBMarkers; }
00299 
00300    TMarker& RefPMAtt()  { return fPMAtt; }
00301    TMarker& RefFVAtt()  { return fFVAtt; }
00302    TMarker& RefPTBAtt() { return fPTBAtt; }
00303    
00304 
00305    static Bool_t IsOutsideBounds(const TEveVector& point, Float_t maxRsqr, Float_t maxZ);
00306 
00307    static Float_t             fgDefMagField; // Default value for constant solenoid magnetic field.
00308    static const Float_t       fgkB2C;        // Constant for conversion of momentum to curvature.
00309    static TEveTrackPropagator fgDefault;     // Default track propagator.
00310 
00311    static Float_t             fgEditorMaxR;  // Max R that can be set in GUI editor.
00312    static Float_t             fgEditorMaxZ;  // Max Z that can be set in GUI editor.
00313 
00314    ClassDef(TEveTrackPropagator, 0); // Calculates path of a particle taking into account special path-marks and imposed boundaries.
00315 };
00316 
00317 //______________________________________________________________________________
00318 inline Bool_t TEveTrackPropagator::IsOutsideBounds(const TEveVector& point,
00319                                                    Float_t           maxRsqr,
00320                                                    Float_t           maxZ)
00321 {
00322    // Return true if point% is outside of cylindrical bounds detrmined by
00323    // square radius and z.
00324 
00325    return TMath::Abs(point.fZ) > maxZ ||
00326           point.fX*point.fX + point.fY*point.fY > maxRsqr;
00327 }
00328 
00329 //______________________________________________________________________________
00330 inline Bool_t TEveTrackPropagator::PointOverVertex(const TEveVector4 &v0,
00331                                                    const TEveVector4 &v,
00332                                                    Float_t           *p)
00333 {
00334    static const Float_t kMinPl = 1e-5;
00335 
00336    TEveVector dv; dv.Sub(v0, v);
00337 
00338    Float_t dotV;
00339 
00340    if (TMath::Abs(fH.fPlMag) > kMinPl)
00341    {
00342       // Use longitudinal momentum to determine crossing point.
00343       // Works ok for spiraling helices, also for loopers.
00344 
00345       dotV = fH.fE1.Dot(dv);
00346       if (fH.fPlMag < 0)
00347          dotV = -dotV;
00348    }
00349    else
00350    {
00351       // Use full momentum, which is pT, under this conditions.
00352 
00353       dotV = fH.fE2.Dot(dv);
00354    }
00355 
00356    if (p)
00357       *p = dotV;
00358 
00359    return dotV < 0;
00360 }
00361 
00362 #endif