+ Modifications of the standard tracking algorithm:
authormartinez <martinez@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 20 Dec 2006 10:58:15 +0000 (10:58 +0000)
committermartinez <martinez@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 20 Dec 2006 10:58:15 +0000 (10:58 +0000)
  - Now we use the track covariance matrix in order to parameterize
    the size of the area where new hits (raw-clusters) are searched;
  - Two extreme options now exist: either we track all possible candidates
    or we keep only the best ones at each step of the  algorithm (using
    the flag AliMUONTrackReconstructor::fgkTrackAllTracks);
+ Removed the class AliMUONSegment which has been replaced by the
  existing class AliMUONObjectPair when needed.
+ Also removed or redefined some hard-coded constants used in the
  tracking algorithms.
             (Ph. Pillot)

20 files changed:
MUON/AliMUONConstants.cxx
MUON/AliMUONConstants.h
MUON/AliMUONHitForRec.cxx
MUON/AliMUONHitForRec.h
MUON/AliMUONTrack.cxx
MUON/AliMUONTrack.h
MUON/AliMUONTrackExtrap.cxx
MUON/AliMUONTrackExtrap.h
MUON/AliMUONTrackK.cxx
MUON/AliMUONTrackK.h
MUON/AliMUONTrackParam.cxx
MUON/AliMUONTrackParam.h
MUON/AliMUONTrackReconstructor.cxx
MUON/AliMUONTrackReconstructor.h
MUON/AliMUONTrackReconstructorK.cxx
MUON/AliMUONTrackReconstructorK.h
MUON/AliMUONVTrackReconstructor.cxx
MUON/AliMUONVTrackReconstructor.h
MUON/MUONrecLinkDef.h
MUON/libMUONrec.pkg

index 6eb28db..c4067e5 100644 (file)
@@ -80,8 +80,31 @@ Float_t  AliMUONConstants::fgDmax[7]  = {  176.6, 229.0, 308.84, 418.2,  522.0,
  
 Int_t    AliMUONConstants::fgMaxZoom = 20;
 
+// Defaults parameters for dipole magnet
+// From ALICE Dimuon - parameters / geometry table,
+// V7-3 (version 7 created 24/03/2004 updated 25/10/2005)
+Double_t AliMUONConstants::fgCoilZ = -994.05;
+Double_t AliMUONConstants::fgCoilL = 502.1;
+Double_t AliMUONConstants::fgYokeZ = -986.6;
+Double_t AliMUONConstants::fgYokeL = 309.4;
+
+// All material constants are taken from AliRoot (version 3.08)  (spectro. (z<0))
+Double_t AliMUONConstants::fgZAbsorberEnd = -503.; // to be coherent with the Geant absorber geometry !!!!
+Int_t   AliMUONConstants::fgNAbsorberElements = 5;
+Double_t AliMUONConstants::fgZAbsorberElement[5] = {-90., -105., -315., -443., -468.};
+Double_t AliMUONConstants::fgX0AbsorberIn[5] = { 30413.000,  // Air
+                                                   24.282,  // C
+                                                   11.274,  // Concrete
+                                                    1.758,  // Fe
+                                                    0.369}; // W (cm)
+Double_t AliMUONConstants::fgX0AbsorberOut[5] = { 24.282,  // C
+                                                 24.282,  // C
+                                                 11.274,  // Concrete
+                                                  1.758,  // Fe 
+                                                  1.758}; // Fe (cm)
+
 // Defaults parameters for track reconstruction
-Double_t  AliMUONConstants::fgDefaultChamberThicknessInX0 = 0.03;
+Double_t  AliMUONConstants::fgChamberThicknessInX0 = 0.03;
 
 //______________________________________________________________________________
 Int_t AliMUONConstants::ChamberNumber(Float_t z) 
index 2710dc2..6725d7d 100644 (file)
@@ -68,8 +68,26 @@ class AliMUONConstants : public TObject {
     static Float_t Pitch()    {return fgPitch;}
     /// Return wire pitch for Station 1 & 2
     static Float_t PitchSt1() {return fgPitchSt1;}
+    /// Return coil z-position
+    static Double_t CoilZ() {return fgCoilZ;}
+    /// Return coil lenght
+    static Double_t CoilL() {return fgCoilL;}
+    /// Return yoke z-position
+    static Double_t YokeZ() {return fgYokeZ;}
+    /// Return yoke lenght
+    static Double_t YokeL() {return fgYokeL;}
+    /// Return the z position of the end of the absorber
+    static Double_t ZAbsorberEnd() {return fgZAbsorberEnd;}
+    /// Return the number of elements making the absorber
+    static Int_t    NAbsorberElements() {return fgNAbsorberElements;}
+    /// Return the z position of the element "iElement" making the absorber
+    static Double_t ZAbsorberElement(Int_t iElement) {return ((iElement >= 0) && (iElement < fgNAbsorberElements)) ? fgZAbsorberElement[iElement] : 0.;}
+    /// Return the radiation lenght of the element "iElement" making the inner part of the absorber
+    static Double_t X0AbsorberIn(Int_t iElement) {return ((iElement >= 0) && (iElement < fgNAbsorberElements)) ? fgX0AbsorberIn[iElement] : -1.;}
+    /// Return the radiation lenght of the element "iElement" making the outer part of the absorber
+    static Double_t X0AbsorberOut(Int_t iElement) {return ((iElement >= 0) && (iElement < fgNAbsorberElements)) ? fgX0AbsorberOut[iElement] : -1.;}
     /// Return chamber thickness in X0
-    static Double_t DefaultChamberThicknessInX0() {return fgDefaultChamberThicknessInX0;}
+    static Double_t ChamberThicknessInX0() {return fgChamberThicknessInX0;}
     /// Return Trigger ToF Limit (75 ns)
     static Float_t TriggerTofLimit() {return fgkTriggerTofLimit;}
  
@@ -88,11 +106,11 @@ class AliMUONConstants : public TObject {
     static Int_t  fgNDetElem;           ///<  Number of Detection Elements.
     static Int_t  fgNGeomModules;       ///<  Number of Geometry modules   
 
-    static Float_t  fgDefaultChamberZ[14];    //!< Z-positions of chambers
-    static Float_t  fgDefaultRatioTriggerChamber[4]; ///< Ratio between trigger chambers
-    static Float_t  fgSt345inclination;       //!< Inclination with respect the vertical axis of stations 345
-    static Float_t  fgDmin[7];                //!< Inner diameter
-    static Float_t  fgDmax[7];                //!< Outer diameter
+    static Float_t  fgDefaultChamberZ[14];             //!< Z-positions of chambers
+    static Float_t  fgDefaultRatioTriggerChamber[4];   ///< Ratio between trigger chambers
+    static Float_t  fgSt345inclination;                        //!< Inclination with respect the vertical axis of stations 345
+    static Float_t  fgDmin[7];                         //!< Inner diameter
+    static Float_t  fgDmax[7];                         //!< Outer diameter
 
     static Float_t  fgDzCh;             ///< Half-distance between two half-chambers 
     static Float_t  fgDzSlat;           ///< Half-distance between two slat on the same chamber
@@ -107,10 +125,21 @@ class AliMUONConstants : public TObject {
     static Float_t  fgPitch;             ///< Wire pitch for St2 & Slats
     static Float_t  fgPitchSt1;          ///< Wire pitch for Station 1
 
-    static Double_t  fgDefaultChamberThicknessInX0; ///< default chamber thickness in X0 for reconstruction
+    static Double_t  fgChamberThicknessInX0; ///< default chamber thickness in X0 for reconstruction
     
-    static Int_t    fgMaxZoom;           ///< Maximum Zoom for event display
-    static Float_t fgkTriggerTofLimit;   ///< Particle above this threshold are discarded in trigger algorithm
+    static Double_t fgCoilZ; ///< Coil z-position
+    static Double_t fgCoilL; ///< Coil lenght
+    static Double_t fgYokeZ; ///< Yoke z-position
+    static Double_t fgYokeL; ///< Yoke lenght
+
+    static Double_t fgZAbsorberEnd;            ///< Z position of the end of the absorber
+    static Int_t    fgNAbsorberElements;       ///< Number of elements making the absorber
+    static Double_t fgZAbsorberElement[5];     ///< Z positions of elements making the absorber
+    static Double_t fgX0AbsorberIn[5];         ///< Radiation lenght of materials in the inner part of the absorber (theta > 3 degrees)
+    static Double_t fgX0AbsorberOut[5];                ///< Radiation lenght of materials in the outer part of the absorber (theta > 3 degrees)
+    
+    static Int_t    fgMaxZoom;          ///< Maximum Zoom for event display
+    static Float_t  fgkTriggerTofLimit; ///< Particle above this threshold are discarded in trigger algorithm
     
     ClassDef(AliMUONConstants, 0) // MUON global constants 
 };
index 4520629..5346b54 100644 (file)
@@ -46,8 +46,6 @@ AliMUONHitForRec::AliMUONHitForRec()
     fHitNumber(0),
     fTTRTrack(0),
     fTrackRefSignal(0),
-    fIndexOfFirstSegment(0),
-    fNSegments(0),
     fNTrackHits(0)
 {
 /// Default Constructor
@@ -67,8 +65,6 @@ AliMUONHitForRec::AliMUONHitForRec(AliTrackReference* theGhit)
     fHitNumber(0),
     fTTRTrack(0),
     fTrackRefSignal(0),
-    fIndexOfFirstSegment(-1),
-    fNSegments(0),
     fNTrackHits(0)
 {
 /// Constructor for AliMUONHitForRec from a track ref. hit.
@@ -115,8 +111,6 @@ AliMUONHitForRec::AliMUONHitForRec(AliMUONRawCluster* theRawCluster)
     fHitNumber(0),
     fTTRTrack(-1),
     fTrackRefSignal(-1),
-    fIndexOfFirstSegment(-1),
-    fNSegments(0),
     fNTrackHits(0)
 {
 /// Constructor for AliMUONHitForRec from a raw cluster.
@@ -140,8 +134,6 @@ AliMUONHitForRec::AliMUONHitForRec (const AliMUONHitForRec& theMUONHitForRec)
     fHitNumber(theMUONHitForRec.fHitNumber),
     fTTRTrack(theMUONHitForRec.fTTRTrack),
     fTrackRefSignal(theMUONHitForRec.fTrackRefSignal),
-    fIndexOfFirstSegment(theMUONHitForRec.fIndexOfFirstSegment),
-    fNSegments(theMUONHitForRec.fNSegments),
     fNTrackHits(theMUONHitForRec.fNTrackHits)
 {
 /// Copy constructor
@@ -163,8 +155,6 @@ AliMUONHitForRec & AliMUONHitForRec::operator=(const AliMUONHitForRec& theMUONHi
   fHitNumber = theMUONHitForRec.fHitNumber;
   fTTRTrack = theMUONHitForRec.fTTRTrack;
   fTrackRefSignal = theMUONHitForRec.fTrackRefSignal;
-  fIndexOfFirstSegment = theMUONHitForRec.fIndexOfFirstSegment;
-  fNSegments = theMUONHitForRec.fNSegments;
   fNTrackHits = theMUONHitForRec.fNTrackHits;
   return *this;
 }
@@ -195,29 +185,23 @@ Int_t AliMUONHitForRec::Compare(const TObject* Hit) const
   //__________________________________________________________________________
 Double_t AliMUONHitForRec::NormalizedChi2WithHitForRec(AliMUONHitForRec* hitForRec, Double_t Sigma2Cut) const
 {
-/// Calculate the normalized Chi2 between the current hitForRec (this)
-/// and the hitForRec pointed to by "hitForRec",
-/// i.e. the square deviations between the coordinates,
-/// in both the bending and the non bending plane,
+/// Calculate the normalized Chi2 between the current hitForRec (this) and the hitForRec pointed to by "hitForRec",
+/// i.e. the square deviations between the coordinates, in both the bending and the non bending plane,
 /// divided by the variance of the same quantities and by "Sigma2Cut".
-/// Returns 3 if none of the 2 quantities is OK,
-/// something smaller than or equal to 2 otherwise.
-/// Would it be more correct to use a real chi square
-/// including the non diagonal term ????
+/// Returns 3 if none of the 2 quantities is OK, something smaller than or equal to 2 otherwise.
+/// Would it be more correct to use a real chi square including the non diagonal term ????
 
   Double_t chi2, chi2Max, diff, normDiff;
   chi2 = 0.0;
   chi2Max = 3.0;
   // coordinate in bending plane
   diff = fBendingCoor - hitForRec->fBendingCoor;
-  normDiff = diff * diff /
-    (fBendingReso2 + hitForRec->fBendingReso2) / Sigma2Cut;
+  normDiff = diff * diff / (fBendingReso2 + hitForRec->fBendingReso2) / Sigma2Cut;
   if (normDiff > 1.0) return chi2Max;
   chi2 = chi2 + normDiff;
   // coordinate in non bending plane
   diff = fNonBendingCoor - hitForRec->fNonBendingCoor;
-  normDiff = diff * diff /
-    (fNonBendingReso2 + hitForRec->fNonBendingReso2) / Sigma2Cut;
+  normDiff = diff * diff / (fNonBendingReso2 + hitForRec->fNonBendingReso2) / Sigma2Cut;
   if (normDiff > 1.0) return chi2Max;
   chi2 = chi2 + normDiff;
   return chi2;
index 02c034b..b9a39e5 100644 (file)
@@ -47,16 +47,10 @@ class AliMUONHitForRec : public TObject {
   void SetTTRTrack(Int_t TTRTrack) { fTTRTrack = TTRTrack;}
   Int_t GetTrackRefSignal(void) const { return fTrackRefSignal;}
   void SetTrackRefSignal(Int_t TrackRefSignal) { fTrackRefSignal = TrackRefSignal;}
-  Int_t GetIndexOfFirstSegment(void) const { return fIndexOfFirstSegment;}
-  void SetIndexOfFirstSegment(Int_t IndexOfFirstSegment) { fIndexOfFirstSegment = IndexOfFirstSegment;}
-  Int_t GetNSegments(void) const { return fNSegments;}
-  void SetNSegments(Int_t NSegments) { fNSegments = NSegments;}
   Int_t GetNTrackHits(void) const { return fNTrackHits;}
   void SetNTrackHits(Int_t NTrackHits) { fNTrackHits = NTrackHits;}
 
-
   Double_t NormalizedChi2WithHitForRec(AliMUONHitForRec* Hit, Double_t Sigma2Cut) const;
-/*   void UpdateFromChamberTrackParam(AliMUONTrackParam *TrackParam, Double_t MCSfactor); */
 
   // What is necessary for sorting TClonesArray's; sufficient too ????
   Bool_t IsSortable() const { return kTRUE; }
@@ -80,10 +74,6 @@ class AliMUONHitForRec : public TObject {
   Int_t fTTRTrack; ///< track number (0...) in TTR
   Int_t fTrackRefSignal; ///< Track ref. signal (1) or background (0)
 
-  // links forward to the segment(s) if HitForRec in first chamber of a station
-  Int_t fIndexOfFirstSegment; //!<  index of first Segment
-  Int_t fNSegments; //!<  number of Segments
-
   Int_t fNTrackHits; //!<  number of TrackHit's made with HitForRec
   
   ClassDef(AliMUONHitForRec, 2) // Hit for reconstruction in ALICE dimuon spectrometer
index 4ee9392..5859b95 100644 (file)
 
 #include <stdlib.h> // for exit()
 #include <Riostream.h>
+#include <TMatrixD.h>
 
 #include "AliMUONTrack.h"
 
 #include "AliMUONTrackParam.h" 
 #include "AliMUONHitForRec.h" 
-#include "AliMUONSegment.h" 
 #include "AliMUONConstants.h"
+#include "AliMUONTrackExtrap.h" 
 
 #include "AliLog.h"
 
@@ -43,6 +44,9 @@
 ClassImp(AliMUONTrack) // Class implementation in ROOT context
 /// \endcond
 
+const Double_t AliMUONTrack::fgkMaxTrackingDistanceBending    = 2.;
+const Double_t AliMUONTrack::fgkMaxTrackingDistanceNonBending = 2.;
+
 //__________________________________________________________________________
 AliMUONTrack::AliMUONTrack()
   : TObject(),
@@ -50,6 +54,9 @@ AliMUONTrack::AliMUONTrack()
     fTrackParamAtHit(0x0),
     fHitForRecAtHit(0x0),
     fNTrackHits(0),
+    fExtrapTrackParam(),
+    fFitWithVertex(kFALSE),
+    fVertex(0x0),
     fFitFMin(-1.),
     fMatchTrigger(kFALSE),
     fChi2MatchTrigger(0.),
@@ -59,52 +66,72 @@ AliMUONTrack::AliMUONTrack()
 }
 
   //__________________________________________________________________________
-AliMUONTrack::AliMUONTrack(AliMUONSegment* BegSegment, AliMUONSegment* EndSegment)
+AliMUONTrack::AliMUONTrack(AliMUONHitForRec* hitForRec1, AliMUONHitForRec* hitForRec2)
   : TObject(),
     fTrackParamAtVertex(),
     fTrackParamAtHit(0x0),
     fHitForRecAtHit(0x0),
     fNTrackHits(0),
+    fExtrapTrackParam(),
+    fFitWithVertex(kFALSE),
+    fVertex(0x0),
     fFitFMin(-1.),
     fMatchTrigger(kFALSE),
     fChi2MatchTrigger(0.),
     fTrackID(0)
 {
-  /// Constructor from two Segment's
+  /// Constructor from thw hitForRec's
 
   fTrackParamAtHit = new TClonesArray("AliMUONTrackParam",10);
   fHitForRecAtHit = new TClonesArray("AliMUONHitForRec",10);
   
-  if (BegSegment) { //AZ
-    AddTrackParamAtHit(0,BegSegment->GetHitForRec1());
-    AddTrackParamAtHit(0,BegSegment->GetHitForRec2());
-    AddTrackParamAtHit(0,EndSegment->GetHitForRec1());
-    AddTrackParamAtHit(0,EndSegment->GetHitForRec2());
-    fTrackParamAtHit->Sort(); // sort TrackParamAtHit according to increasing Z
-  }
-}
-
-  //__________________________________________________________________________
-AliMUONTrack::AliMUONTrack(AliMUONSegment* Segment, AliMUONHitForRec* HitForRec)
-  : TObject(),
-    fTrackParamAtVertex(),
-    fTrackParamAtHit(0x0),
-    fHitForRecAtHit(0x0),
-    fNTrackHits(0),
-    fFitFMin(-1.),
-    fMatchTrigger(kFALSE),
-    fChi2MatchTrigger(0.),
-    fTrackID(0)
-{
-  /// Constructor from one Segment and one HitForRec
-
-  fTrackParamAtHit = new TClonesArray("AliMUONTrackParam",10);
-  fHitForRecAtHit = new TClonesArray("AliMUONHitForRec",10);
+  if (!hitForRec1) return; //AZ
+  
+  // Add hits to the track
+  AddTrackParamAtHit(0,hitForRec1);
+  AddTrackParamAtHit(0,hitForRec2);
+  
+  // sort TrackParamAtHit according to increasing -Z
+  fTrackParamAtHit->Sort();
+  
+  // Set track parameters at first track hit
+  AliMUONTrackParam* trackParamAtFirstHit = (AliMUONTrackParam*) fTrackParamAtHit->First();
+  AliMUONHitForRec* firstHit = trackParamAtFirstHit->GetHitForRecPtr();
+  AliMUONHitForRec* lastHit = ((AliMUONTrackParam*) fTrackParamAtHit->Last())->GetHitForRecPtr();
+  // Z position
+  Double_t z1 = firstHit->GetZ();
+  trackParamAtFirstHit->SetZ(z1);
+  Double_t dZ = z1 - lastHit->GetZ();
+  // Non bending plane
+  Double_t nonBendingCoor = firstHit->GetNonBendingCoor();
+  trackParamAtFirstHit->SetNonBendingCoor(nonBendingCoor);
+  trackParamAtFirstHit->SetNonBendingSlope((nonBendingCoor - lastHit->GetNonBendingCoor()) / dZ);
+  // Bending plane
+  Double_t bendingCoor = firstHit->GetBendingCoor();
+  trackParamAtFirstHit->SetBendingCoor(bendingCoor);
+  Double_t bendingSlope = (bendingCoor - lastHit->GetBendingCoor()) / dZ;
+  trackParamAtFirstHit->SetBendingSlope(bendingSlope);
+  // Inverse bending momentum
+  Double_t bendingImpact = bendingCoor - z1 * bendingSlope;
+  Double_t inverseBendingMomentum = 1. / AliMUONTrackExtrap::GetBendingMomentumFromImpactParam(bendingImpact);
+  trackParamAtFirstHit->SetInverseBendingMomentum(inverseBendingMomentum);
+  
+  // Evaluate covariances
+  TMatrixD *paramCov = trackParamAtFirstHit->GetCovariances();
+  (*paramCov) = 0;
+  // Non bending plane
+  (*paramCov)(0,0) = firstHit->GetNonBendingReso2();
+  (*paramCov)(0,1) = firstHit->GetNonBendingReso2() / dZ;
+  (*paramCov)(1,0) = (*paramCov)(0,1);
+  (*paramCov)(1,1) = ( firstHit->GetNonBendingReso2() + lastHit->GetNonBendingReso2() ) / dZ / dZ;
+  // Bending plane
+  (*paramCov)(2,2) = firstHit->GetBendingReso2();
+  (*paramCov)(2,3) = firstHit->GetBendingReso2() / dZ;
+  (*paramCov)(3,2) = (*paramCov)(2,3);
+  (*paramCov)(3,3) = ( firstHit->GetBendingReso2() + lastHit->GetBendingReso2() ) / dZ / dZ;
+  // Inverse bending momentum (50% error)
+  (*paramCov)(4,4) = 0.5*inverseBendingMomentum * 0.5*inverseBendingMomentum;
   
-  AddTrackParamAtHit(0,Segment->GetHitForRec1());
-  AddTrackParamAtHit(0,Segment->GetHitForRec2());
-  AddTrackParamAtHit(0,HitForRec);
-  fTrackParamAtHit->Sort(); // sort TrackParamAtHit according to increasing Z
 }
 
   //__________________________________________________________________________
@@ -114,13 +141,19 @@ AliMUONTrack::~AliMUONTrack()
   if (fTrackParamAtHit) {
     // delete the TClonesArray of pointers to TrackParam
     delete fTrackParamAtHit;
-    fTrackParamAtHit = NULL;
+    fTrackParamAtHit = 0x0;
   }
 
   if (fHitForRecAtHit) {
     // delete the TClonesArray of pointers to HitForRec
     delete fHitForRecAtHit;
-    fHitForRecAtHit = NULL;
+    fHitForRecAtHit = 0x0;
+  }
+  
+  if (fVertex) {
+    // delete the vertex used during the tracking procedure
+    delete fVertex;
+    fVertex = 0x0;
   }
 }
 
@@ -131,6 +164,9 @@ AliMUONTrack::AliMUONTrack (const AliMUONTrack& theMUONTrack)
     fTrackParamAtHit(0x0),
     fHitForRecAtHit(0x0),
     fNTrackHits(theMUONTrack.fNTrackHits),
+    fExtrapTrackParam(theMUONTrack.fExtrapTrackParam),
+    fFitWithVertex(theMUONTrack.fFitWithVertex),
+    fVertex(0x0),
     fFitFMin(theMUONTrack.fFitFMin),
     fMatchTrigger(theMUONTrack.fMatchTrigger),
     fChi2MatchTrigger(theMUONTrack.fChi2MatchTrigger),
@@ -142,21 +178,24 @@ AliMUONTrack::AliMUONTrack (const AliMUONTrack& theMUONTrack)
   // necessary to make a copy of the objects and not only the pointers in TClonesArray.
   if (theMUONTrack.fTrackParamAtHit) {
     maxIndex = (theMUONTrack.fTrackParamAtHit)->GetEntriesFast();
-    fTrackParamAtHit  =  new TClonesArray("AliMUONTrackParam",maxIndex);
+    fTrackParamAtHit = new TClonesArray("AliMUONTrackParam",maxIndex);
     for (Int_t index = 0; index < maxIndex; index++) {
       new ((*fTrackParamAtHit)[index]) AliMUONTrackParam(*(AliMUONTrackParam*)theMUONTrack.fTrackParamAtHit->At(index));
     }
-  }  
-
+  }
+  
   // necessary to make a copy of the objects and not only the pointers in TClonesArray.
   if (theMUONTrack.fHitForRecAtHit) {
     maxIndex = (theMUONTrack.fHitForRecAtHit)->GetEntriesFast();
-    fHitForRecAtHit  =  new TClonesArray("AliMUONHitForRec",maxIndex);
+    fHitForRecAtHit = new TClonesArray("AliMUONHitForRec",maxIndex);
     for (Int_t index = 0; index < maxIndex; index++) {
       new ((*fHitForRecAtHit)[index]) AliMUONHitForRec(*(AliMUONHitForRec*)theMUONTrack.fHitForRecAtHit->At(index));
     }
-  }  
-
+  }
+  
+  // copy vertex used during the tracking procedure if any
+  if (theMUONTrack.fVertex) fVertex = new AliMUONHitForRec(*(theMUONTrack.fVertex));
+  
 }
 
   //__________________________________________________________________________
@@ -175,28 +214,46 @@ AliMUONTrack & AliMUONTrack::operator=(const AliMUONTrack& theMUONTrack)
   Int_t maxIndex = 0;
   
   // necessary to make a copy of the objects and not only the pointers in TClonesArray.
-  fTrackParamAtHit = 0;
   if (theMUONTrack.fTrackParamAtHit) {
-    fTrackParamAtHit  =  new TClonesArray("AliMUONTrackParam",10);
+    if (fTrackParamAtHit) fTrackParamAtHit->Clear();
+    else fTrackParamAtHit = new TClonesArray("AliMUONTrackParam",10);
     maxIndex = (theMUONTrack.fTrackParamAtHit)->GetEntriesFast();
     for (Int_t index = 0; index < maxIndex; index++) {
       new ((*fTrackParamAtHit)[fTrackParamAtHit->GetEntriesFast()])
        AliMUONTrackParam(*(AliMUONTrackParam*)(theMUONTrack.fTrackParamAtHit)->At(index));
     }
-  }  
+  } else if (fTrackParamAtHit) {
+    delete fTrackParamAtHit;
+    fTrackParamAtHit = 0x0;
+  }
 
   // necessary to make a copy of the objects and not only the pointers in TClonesArray.
-  fHitForRecAtHit = 0;
   if (theMUONTrack.fHitForRecAtHit) {
-    fHitForRecAtHit  =  new TClonesArray("AliMUONHitForRec",10);
+    if (fHitForRecAtHit) fHitForRecAtHit->Clear();
+    else fHitForRecAtHit = new TClonesArray("AliMUONHitForRec",10);
     maxIndex = (theMUONTrack.fHitForRecAtHit)->GetEntriesFast();
     for (Int_t index = 0; index < maxIndex; index++) {
       new ((*fHitForRecAtHit)[fHitForRecAtHit->GetEntriesFast()])
        AliMUONHitForRec(*(AliMUONHitForRec*)(theMUONTrack.fHitForRecAtHit)->At(index));
     }
-  }  
+  } else if (fHitForRecAtHit) {
+    delete fHitForRecAtHit;
+    fHitForRecAtHit = 0x0;
+  }
+  
+  // copy vertex used during the tracking procedure if any.
+  if (theMUONTrack.fVertex) {
+    if (fVertex) *fVertex = *(theMUONTrack.fVertex);
+    else fVertex = new AliMUONHitForRec(*(theMUONTrack.fVertex));
+  } else if (fVertex) {
+    delete fVertex;
+    fVertex = 0x0;
+  }
+  
+  fExtrapTrackParam = theMUONTrack.fExtrapTrackParam;
   
   fNTrackHits         =  theMUONTrack.fNTrackHits;
+  fFitWithVertex      =  theMUONTrack.fFitWithVertex;
   fFitFMin            =  theMUONTrack.fFitFMin;
   fMatchTrigger       =  theMUONTrack.fMatchTrigger;
   fChi2MatchTrigger   =  theMUONTrack.fChi2MatchTrigger;
@@ -235,7 +292,40 @@ void AliMUONTrack::AddHitForRecAtHit(const AliMUONHitForRec *hitForRec)
 }
 
   //__________________________________________________________________________
-Bool_t* AliMUONTrack::CompatibleTrack(AliMUONTrack * Track, Double_t Sigma2Cut) const
+void AliMUONTrack::SetVertex(AliMUONHitForRec* vertex)
+{
+  /// Set the vertex used during the tracking procedure
+  if (!fVertex) fVertex = new AliMUONHitForRec(*vertex);
+  else *fVertex = *vertex;
+}
+
+  //__________________________________________________________________________
+Int_t AliMUONTrack::HitsInCommon(AliMUONTrack* track) const
+{
+  /// Returns the number of hits in common between the current track ("this")
+  /// and the track pointed to by "track".
+  Int_t hitsInCommon = 0;
+  AliMUONTrackParam *trackParamAtHit1, *trackParamAtHit2;
+  // Loop over hits of first track
+  trackParamAtHit1 = (AliMUONTrackParam*) this->fTrackParamAtHit->First();
+  while (trackParamAtHit1) {
+    // Loop over hits of second track
+    trackParamAtHit2 = (AliMUONTrackParam*) track->fTrackParamAtHit->First();
+    while (trackParamAtHit2) {
+      // Increment "hitsInCommon" if both TrackParamAtHits point to the same HitForRec
+      if ((trackParamAtHit1->GetHitForRecPtr()) == (trackParamAtHit2->GetHitForRecPtr())) {
+        hitsInCommon++;
+       break;
+      }
+      trackParamAtHit2 = (AliMUONTrackParam*) track->fTrackParamAtHit->After(trackParamAtHit2);
+    } // trackParamAtHit2
+    trackParamAtHit1 = (AliMUONTrackParam*) this->fTrackParamAtHit->After(trackParamAtHit1);
+  } // trackParamAtHit1
+  return hitsInCommon;
+}
+
+  //__________________________________________________________________________
+Bool_t* AliMUONTrack::CompatibleTrack(AliMUONTrack * track, Double_t sigma2Cut) const
 {
   /// Return kTRUE/kFALSE for each chamber if hit is compatible or not 
   TClonesArray *hitArray, *thisHitArray;
@@ -252,7 +342,7 @@ Bool_t* AliMUONTrack::CompatibleTrack(AliMUONTrack * Track, Double_t Sigma2Cut)
 
   thisHitArray = this->GetHitForRecAtHit();
 
-  hitArray =  Track->GetHitForRecAtHit();
+  hitArray =  track->GetHitForRecAtHit();
 
   for (Int_t iHthis = 0; iHthis < thisHitArray->GetEntriesFast(); iHthis++) {
     thisHit = (AliMUONHitForRec*) thisHitArray->At(iHthis);
@@ -262,7 +352,7 @@ Bool_t* AliMUONTrack::CompatibleTrack(AliMUONTrack * Track, Double_t Sigma2Cut)
     for (Int_t iH = 0; iH < hitArray->GetEntriesFast(); iH++) {
       hit = (AliMUONHitForRec*) hitArray->At(iH);
       deltaZ = TMath::Abs(thisHit->GetZ() - hit->GetZ());
-      chi2 = thisHit->NormalizedChi2WithHitForRec(hit,Sigma2Cut); // set cut to 4 sigmas
+      chi2 = thisHit->NormalizedChi2WithHitForRec(hit,sigma2Cut); // set cut to 4 sigmas
       if (chi2 < 3. && deltaZ < deltaZMax) {
        nCompHit[chamberNumber] = kTRUE;
        break;
@@ -274,25 +364,153 @@ Bool_t* AliMUONTrack::CompatibleTrack(AliMUONTrack * Track, Double_t Sigma2Cut)
 }
 
   //__________________________________________________________________________
-Int_t AliMUONTrack::HitsInCommon(AliMUONTrack* Track) const
+Double_t AliMUONTrack::TryOneHitForRec(AliMUONHitForRec* hitForRec)
 {
-  /// Returns the number of hits in common between the current track ("this")
-  /// and the track pointed to by "Track".
-  Int_t hitsInCommon = 0;
-  AliMUONTrackParam *trackParamAtHit1, *trackParamAtHit2;
-  // Loop over hits of first track
-  trackParamAtHit1 = (AliMUONTrackParam*) this->fTrackParamAtHit->First();
-  while (trackParamAtHit1) {
-    // Loop over hits of second track
-    trackParamAtHit2 = (AliMUONTrackParam*) Track->fTrackParamAtHit->First();
-    while (trackParamAtHit2) {
-      // Increment "hitsInCommon" if both TrackParamAtHits point to the same HitForRec
-      if ((trackParamAtHit1->GetHitForRecPtr()) == (trackParamAtHit2->GetHitForRecPtr())) hitsInCommon++;
-      trackParamAtHit2 = (AliMUONTrackParam*) Track->fTrackParamAtHit->After(trackParamAtHit2);
-    } // trackParamAtHit2
-    trackParamAtHit1 = (AliMUONTrackParam*) this->fTrackParamAtHit->After(trackParamAtHit1);
-  } // trackParamAtHit1
-  return hitsInCommon;
+/// Test the compatibility between the track and the hitForRec:
+/// return the corresponding Chi2
+  
+  // Get track parameters and their covariances at the z position of hitForRec
+  AliMUONTrackParam extrapTrackParam(fExtrapTrackParam);
+  AliMUONTrackExtrap::ExtrapToZCov(&extrapTrackParam, hitForRec->GetZ());
+  
+  // Set differences between trackParam and hitForRec in the bending and non bending directions
+  TMatrixD dPos(2,1);
+  dPos(0,0) = hitForRec->GetNonBendingCoor() - extrapTrackParam.GetNonBendingCoor();
+  dPos(1,0) = hitForRec->GetBendingCoor() - extrapTrackParam.GetBendingCoor();
+  
+  // quick test of hitForRec compatibility within a wide road of x*y = 10*1 cm2 to save computing time
+  if (TMath::Abs(dPos(0,0)) > fgkMaxTrackingDistanceNonBending ||
+      TMath::Abs(dPos(1,0)) > fgkMaxTrackingDistanceBending) return 1.e10;
+  
+  // Set the error matrix from trackParam covariances and hitForRec resolution
+  TMatrixD* paramCov = extrapTrackParam.GetCovariances();
+  TMatrixD error(2,2);
+  error(0,0) = (*paramCov)(0,0) + hitForRec->GetNonBendingReso2();
+  error(0,1) = (*paramCov)(0,2);
+  error(1,0) = (*paramCov)(2,0);
+  error(1,1) = (*paramCov)(2,2) + hitForRec->GetBendingReso2();
+  
+  // Invert the error matrix for Chi2 calculation
+  if (error.Determinant() != 0) {
+    error.Invert();
+  } else {
+    AliWarning(" Determinant error=0");
+    return 1.e10;
+  }
+  
+  // Compute the Chi2 value
+  TMatrixD tmp(error,TMatrixD::kMult,dPos);
+  TMatrixD result(dPos,TMatrixD::kTransposeMult,tmp);
+  
+  return result(0,0);
+  
+}
+
+  //__________________________________________________________________________
+Double_t AliMUONTrack::TryTwoHitForRec(AliMUONHitForRec* hitForRec1, AliMUONHitForRec* hitForRec2)
+{
+/// Test the compatibility between the track and the 2 hitForRec together:
+/// return the corresponding Chi2 accounting for covariances between the 2 hitForRec
+  
+  // Get track parameters and their covariances at the z position of the first hitForRec
+  AliMUONTrackParam extrapTrackParam1(fExtrapTrackParam);
+  AliMUONTrackExtrap::ExtrapToZCov(&extrapTrackParam1, hitForRec1->GetZ());
+  
+  // Get track parameters at second hitForRec
+  AliMUONTrackParam extrapTrackParam2(extrapTrackParam1);
+  AliMUONTrackExtrap::ExtrapToZ(&extrapTrackParam2, hitForRec2->GetZ());
+  
+  // Set differences between track and the 2 hitForRec in the bending and non bending directions
+  TMatrixD dPos(4,1);
+  dPos(0,0) = hitForRec1->GetNonBendingCoor() - extrapTrackParam1.GetNonBendingCoor();
+  dPos(1,0) = hitForRec1->GetBendingCoor() - extrapTrackParam1.GetBendingCoor();
+  dPos(2,0) = hitForRec2->GetNonBendingCoor() - extrapTrackParam2.GetNonBendingCoor();
+  dPos(3,0) = hitForRec2->GetBendingCoor() - extrapTrackParam2.GetBendingCoor();
+  
+  // quick tests of hitForRec compatibility within a wide road of x*y = 1*1 cm2 to save computing time
+  if (TMath::Abs(dPos(0,0)) > fgkMaxTrackingDistanceNonBending ||
+      TMath::Abs(dPos(1,0)) > fgkMaxTrackingDistanceBending    ||
+      TMath::Abs(dPos(2,0)) > fgkMaxTrackingDistanceNonBending ||
+      TMath::Abs(dPos(3,0)) > fgkMaxTrackingDistanceBending) return 1.e10;
+  
+  // Calculate the error matrix from the track parameter covariances at first hitForRec
+  TMatrixD error(4,4);
+  error = 0.;
+  if (extrapTrackParam1.CovariancesExist()) {
+    // Get the pointer to the parameter covariance matrix at first hitForRec
+    TMatrixD* paramCov = extrapTrackParam1.GetCovariances();
+    
+    // Save track parameters at first hitForRec
+    AliMUONTrackParam extrapTrackParam1Save(extrapTrackParam1);
+    Double_t nonBendingCoor1        = extrapTrackParam1Save.GetNonBendingCoor();
+    Double_t nonBendingSlope1       = extrapTrackParam1Save.GetNonBendingSlope();
+    Double_t bendingCoor1           = extrapTrackParam1Save.GetBendingCoor();
+    Double_t bendingSlope1          = extrapTrackParam1Save.GetBendingSlope();
+    Double_t inverseBendingMomentum1 = extrapTrackParam1Save.GetInverseBendingMomentum();
+    Double_t z1                             = extrapTrackParam1Save.GetZ();
+    
+    // Save track coordinates at second hitForRec
+    Double_t nonBendingCoor2        = extrapTrackParam2.GetNonBendingCoor();
+    Double_t bendingCoor2           = extrapTrackParam2.GetBendingCoor();
+    
+    // Calculate the jacobian related to the transformation between track parameters
+    // at first hitForRec and track coordinates at the 2 hitForRec z-position
+    TMatrixD jacob(4,5);
+    jacob = 0.;
+    // first derivative at the first hitForRec:
+    jacob(0,0) = 1.; // dx1/dx
+    jacob(1,2) = 1.; // dy1/dy
+    // first derivative at the second hitForRec:
+    Double_t dParam[5];
+    for (Int_t i=0; i<5; i++) {
+      // Skip jacobian calculation for parameters with no associated error
+      if ((*paramCov)(i,i) == 0.) continue;
+      // Small variation of parameter i only
+      for (Int_t j=0; j<5; j++) {
+        if (j==i) {
+          dParam[j] = TMath::Sqrt((*paramCov)(i,i));
+         if (j == 4) dParam[j] *= TMath::Sign(1.,-inverseBendingMomentum1); // variation always in the same direction
+        } else dParam[j] = 0.;
+      }
+      // Set new track parameters at first hitForRec
+      extrapTrackParam1Save.SetNonBendingCoor       (nonBendingCoor1         + dParam[0]);
+      extrapTrackParam1Save.SetNonBendingSlope      (nonBendingSlope1        + dParam[1]);
+      extrapTrackParam1Save.SetBendingCoor          (bendingCoor1            + dParam[2]);
+      extrapTrackParam1Save.SetBendingSlope         (bendingSlope1           + dParam[3]);
+      extrapTrackParam1Save.SetInverseBendingMomentum(inverseBendingMomentum1 + dParam[4]);
+      extrapTrackParam1Save.SetZ                    (z1);
+      // Extrapolate new track parameters to the z position of the second hitForRec
+      AliMUONTrackExtrap::ExtrapToZ(&extrapTrackParam1Save,hitForRec2->GetZ());
+      // Calculate the jacobian
+      jacob(2,i) = (extrapTrackParam1Save.GetNonBendingCoor()  - nonBendingCoor2) / dParam[i]; // dx2/dParami
+      jacob(3,i) = (extrapTrackParam1Save.GetBendingCoor()     - bendingCoor2   ) / dParam[i]; // dy2/dParami
+    }
+    
+    // Calculate the error matrix
+    TMatrixD tmp((*paramCov),TMatrixD::kMultTranspose,jacob);
+    error = TMatrixD(jacob,TMatrixD::kMult,tmp);
+  }
+  
+  // Add hitForRec resolution to the error matrix
+  error(0,0) += hitForRec1->GetNonBendingReso2();
+  error(1,1) += hitForRec1->GetBendingReso2();
+  error(2,2) += hitForRec2->GetNonBendingReso2();
+  error(3,3) += hitForRec2->GetBendingReso2();
+  
+  // invert the error matrix for Chi2 calculation
+  if (error.Determinant() != 0) {
+    error.Invert();
+  } else {
+    AliWarning(" Determinant error=0");
+    return 1.e10;
+  }
+  
+  // Compute the Chi2 value
+  TMatrixD tmp2(error,TMatrixD::kMult,dPos);
+  TMatrixD result(dPos,TMatrixD::kTransposeMult,tmp2);
+  
+  return result(0,0);
+  
 }
 
   //__________________________________________________________________________
index be7f2c5..6b7dc65 100644 (file)
@@ -19,7 +19,6 @@
 #include "AliMUONTrackParam.h" // object belongs to the class
 
 class AliMUONHitForRec;
-class AliMUONSegment;
 
 class AliMUONTrack : public TObject 
 {
@@ -29,13 +28,12 @@ class AliMUONTrack : public TObject
   AliMUONTrack (const AliMUONTrack& AliMUONTrack); // copy constructor
   AliMUONTrack& operator=(const AliMUONTrack& AliMUONTrack); // assignment operator
 
-  AliMUONTrack(AliMUONSegment* BegSegment, AliMUONSegment* EndSegment); // Constructor from two Segment's
-  AliMUONTrack(AliMUONSegment* Segment, AliMUONHitForRec* HitForRec); // Constructor from one Segment and one HitForRec
+  AliMUONTrack(AliMUONHitForRec* hitForRec1, AliMUONHitForRec* hitForRec2); // Constructor from a segment
 
        /// return pointeur to track parameters at vertex
   AliMUONTrackParam*         GetTrackParamAtVertex(void) {return &fTrackParamAtVertex;}
        /// set track parameters at vertex
-  void                       SetTrackParamAtVertex(AliMUONTrackParam* TrackParam) {fTrackParamAtVertex = *TrackParam;}
+  void                       SetTrackParamAtVertex(AliMUONTrackParam* trackParam) {fTrackParamAtVertex = *trackParam;}
 
        /// return array of track parameters at hit
   TClonesArray*              GetTrackParamAtHit(void) const {return fTrackParamAtHit;}
@@ -55,38 +53,62 @@ class AliMUONTrack : public TObject
        /// set the number of hits attached to the track
   void                       SetNTrackHits(Int_t nTrackHits) {fNTrackHits = nTrackHits;}
 
+       /// return pointeur to track parameters extrapolated to the next station
+  AliMUONTrackParam*         GetExtrapTrackParam(void) {return &fExtrapTrackParam;}
+       /// set track parameters extrapolated to next station
+  void                       SetExtrapTrackParam(AliMUONTrackParam* trackParam) {fExtrapTrackParam = *trackParam;}
+
+       /// return kTrue if the vertex must be used to constrain the fit, kFalse if not
+  Bool_t                     GetFitWithVertex(void) const {return fFitWithVertex;}
+       /// set the flag telling whether the vertex must be used to constrain the fit or not
+  void                       SetFitWithVertex(Bool_t fitWithVertex) { fFitWithVertex = fitWithVertex; }
+       /// return the vertex used during the tracking procedure
+  AliMUONHitForRec*          GetVertex(void) const {return fVertex;}
+  void                       SetVertex(AliMUONHitForRec* vertex);
+
        /// return the minimum value of the function minimized by the fit
   Double_t                   GetFitFMin(void) const {return fFitFMin;}
        /// set the minimum value of the function minimized by the fit
-  void                       SetFitFMin(Double_t chi2) { fFitFMin = chi2; } // set Chi2
+  void                       SetFitFMin(Double_t chi2) { fFitFMin = chi2; }
        /// return kTrue if track matches with trigger track, kFalse if not
   Bool_t                     GetMatchTrigger(void) const {return fMatchTrigger;}
        /// set the flag telling whether track matches with trigger track or not
-  void                      SetMatchTrigger(Bool_t MatchTrigger) {fMatchTrigger = MatchTrigger;}
+  void                      SetMatchTrigger(Bool_t matchTrigger) {fMatchTrigger = matchTrigger;}
        /// return the chi2 of trigger/track matching 
   Double_t                   GetChi2MatchTrigger(void) const {return fChi2MatchTrigger;}
        /// set the chi2 of trigger/track matching 
-  void                       SetChi2MatchTrigger(Double_t Chi2MatchTrigger) {fChi2MatchTrigger = Chi2MatchTrigger;}
+  void                       SetChi2MatchTrigger(Double_t chi2MatchTrigger) {fChi2MatchTrigger = chi2MatchTrigger;}
   
-  Int_t                      HitsInCommon(AliMUONTrack* Track) const;
-  Bool_t*                    CompatibleTrack(AliMUONTrack* Track, Double_t Sigma2Cut) const; // return array of compatible chamber
+  Int_t                      HitsInCommon(AliMUONTrack* track) const;
+  Bool_t*                    CompatibleTrack(AliMUONTrack* track, Double_t sigma2Cut) const; // return array of compatible chamber
   
        /// return track number in TrackRefs
   Int_t                      GetTrackID() const {return fTrackID;}
        /// set track number in TrackRefs
   void                       SetTrackID(Int_t trackID) {fTrackID = trackID;}
 
+  Double_t                   TryOneHitForRec(AliMUONHitForRec* hitForRec);
+  Double_t                   TryTwoHitForRec(AliMUONHitForRec* hitForRec1, AliMUONHitForRec* hitForRec2); 
+  
   void                       RecursiveDump(void) const; // Recursive dump (with track hits)
 
   virtual void               Print(Option_t* opt="") const;
 
 
  private:
+  static const Double_t fgkMaxTrackingDistanceBending;    ///< Maximum distance to the track to search for compatible hitForRec(s) in bending direction
+  static const Double_t fgkMaxTrackingDistanceNonBending; ///< Maximum distance to the track to search for compatible hitForRec(s) in non bending direction
+  
   AliMUONTrackParam fTrackParamAtVertex; ///< Track parameters at vertex
   TClonesArray *fTrackParamAtHit; ///< Track parameters at hit
   TClonesArray *fHitForRecAtHit; ///< Cluster parameters at hit
   Int_t fNTrackHits; ///< Number of hits attached to the track
   
+  AliMUONTrackParam fExtrapTrackParam; //!< Track parameters extrapolated to a given z position
+  
+  Bool_t fFitWithVertex; //!< 1 if using the vertex to constrain the fit, 0 if not
+  AliMUONHitForRec *fVertex; //!< Vertex used during the tracking procedure if required
+  
   Double_t fFitFMin; ///< minimum value of the function minimized by the fit
   Bool_t fMatchTrigger; ///< 1 if track matches with trigger track, 0 if not
   Double_t fChi2MatchTrigger; ///< chi2 of trigger/track matching 
index f4733bf..bb70f71 100644 (file)
@@ -29,6 +29,7 @@
 ///////////////////////////////////////////////////
 
 #include <Riostream.h>
+#include <TMatrixD.h>
 
 #include "AliMUONTrackExtrap.h" 
 #include "AliMUONTrackParam.h"
 ClassImp(AliMUONTrackExtrap) // Class implementation in ROOT context
 
 const AliMagF* AliMUONTrackExtrap::fgkField = 0x0;
+const Int_t    AliMUONTrackExtrap::fgkMaxStepNumber = 5000;
+const Double_t AliMUONTrackExtrap::fgkStepLength = 6.;
+
+  //__________________________________________________________________________
+Double_t AliMUONTrackExtrap::GetImpactParamFromBendingMomentum(Double_t bendingMomentum)
+{
+  /// Returns impact parameter at vertex in bending plane (cm),
+  /// from the signed bending momentum "BendingMomentum" in bending plane (GeV/c),
+  /// using simple values for dipole magnetic field.
+  /// The sign of "BendingMomentum" is the sign of the charge.
+  
+  if (bendingMomentum == 0.) return 1.e10;
+  
+  Double_t simpleBPosition = 0.5 * (AliMUONConstants::CoilZ() + AliMUONConstants::YokeZ());
+  Double_t simpleBLength = 0.5 * (AliMUONConstants::CoilL() + AliMUONConstants::YokeL());
+  Float_t b[3], x[3] = {0.,0.,(Float_t) simpleBPosition};
+  if (fgkField) fgkField->Field(x,b);
+  else {
+    cout<<"F-AliMUONTrackExtrap::GetField: fgkField = 0x0"<<endl;
+    exit(-1);
+  }
+  Double_t simpleBValue = (Double_t) b[0];
+  
+  return (-0.0003 * simpleBValue * simpleBLength * simpleBPosition / bendingMomentum);
+}
+
+  //__________________________________________________________________________
+Double_t AliMUONTrackExtrap::GetBendingMomentumFromImpactParam(Double_t impactParam)
+{
+  /// Returns signed bending momentum in bending plane (GeV/c),
+  /// the sign being the sign of the charge for particles moving forward in Z,
+  /// from the impact parameter "ImpactParam" at vertex in bending plane (cm),
+  /// using simple values for dipole magnetic field.
+  
+  if (impactParam == 0.) return 1.e10;
+  
+  Double_t simpleBPosition = 0.5 * (AliMUONConstants::CoilZ() + AliMUONConstants::YokeZ());
+  Double_t simpleBLength = 0.5 * (AliMUONConstants::CoilL() + AliMUONConstants::YokeL());
+  Float_t b[3], x[3] = {0.,0.,(Float_t) simpleBPosition};
+  if (fgkField) fgkField->Field(x,b);
+  else {
+    cout<<"F-AliMUONTrackExtrap::GetField: fgkField = 0x0"<<endl;
+    exit(-1);
+  }
+  Double_t simpleBValue = (Double_t) b[0];
+  
+  return (-0.0003 * simpleBValue * simpleBLength * simpleBPosition / impactParam);
+}
 
   //__________________________________________________________________________
 void AliMUONTrackExtrap::ExtrapToZ(AliMUONTrackParam* trackParam, Double_t zEnd)
@@ -52,29 +101,24 @@ void AliMUONTrackExtrap::ExtrapToZ(AliMUONTrackParam* trackParam, Double_t zEnd)
   else forwardBackward = -1.0;
   Double_t v3[7], v3New[7]; // 7 in parameter ????
   Int_t i3, stepNumber;
-  Int_t maxStepNumber = 5000; // in parameter ????
   // For safety: return kTRUE or kFALSE ????
   // Parameter vector for calling EXTRAP_ONESTEP
   ConvertTrackParamForExtrap(trackParam, v3, forwardBackward);
   // sign of charge (sign of fInverseBendingMomentum if forward motion)
   // must be changed if backward extrapolation
-  Double_t chargeExtrap = forwardBackward *
-    TMath::Sign(Double_t(1.0), trackParam->GetInverseBendingMomentum());
-  Double_t stepLength = 6.0; // in parameter ????
+  Double_t chargeExtrap = forwardBackward * TMath::Sign(Double_t(1.0), trackParam->GetInverseBendingMomentum());
   // Extrapolation loop
   stepNumber = 0;
-  while (((-forwardBackward * (v3[2] - zEnd)) <= 0.0) &&  // spectro. z<0
-        (stepNumber < maxStepNumber)) {
+  while (((-forwardBackward * (v3[2] - zEnd)) <= 0.0) && (stepNumber < fgkMaxStepNumber)) { // spectro. z<0
     stepNumber++;
     // Option for switching between helix and Runge-Kutta ???? 
-    //ExtrapOneStepRungekutta(chargeExtrap, stepLength, v3, v3New);
-    ExtrapOneStepHelix(chargeExtrap, stepLength, v3, v3New);
+    //ExtrapOneStepRungekutta(chargeExtrap, fgkStepLength, v3, v3New);
+    ExtrapOneStepHelix(chargeExtrap, fgkStepLength, v3, v3New);
     if ((-forwardBackward * (v3New[2] - zEnd)) > 0.0) break; // one is beyond Z spectro. z<0
     // better use TArray ????
-    for (i3 = 0; i3 < 7; i3++)
-      {v3[i3] = v3New[i3];}
+    for (i3 = 0; i3 < 7; i3++) {v3[i3] = v3New[i3];}
   }
-  // check maxStepNumber ????
+  // check fgkMaxStepNumber ????
   // Interpolation back to exact Z (2nd order)
   // should be in function ???? using TArray ????
   Double_t dZ12 = v3New[2] - v3[2]; // 1->2
@@ -91,8 +135,7 @@ void AliMUONTrackExtrap::ExtrapToZ(AliMUONTrackParam* trackParam, Double_t zEnd)
     Double_t xPrimeI = xPrime - 0.5 * xSecond * (dZi2 - dZ1i);
     Double_t yPrimeI = yPrime - 0.5 * ySecond * (dZi2 - dZ1i);
     // (PX, PY, PZ)/PTOT assuming forward motion
-    v3[5] =
-      1.0 / TMath::Sqrt(1.0 + xPrimeI * xPrimeI + yPrimeI * yPrimeI); // PZ/PTOT
+    v3[5] = 1.0 / TMath::Sqrt(1.0 + xPrimeI * xPrimeI + yPrimeI * yPrimeI); // PZ/PTOT
     v3[3] = xPrimeI * v3[5]; // PX/PTOT
     v3[4] = yPrimeI * v3[5]; // PY/PTOT
   } else {
@@ -136,6 +179,76 @@ void AliMUONTrackExtrap::RecoverTrackParam(Double_t *v3, Double_t charge, AliMUO
 }
 
   //__________________________________________________________________________
+void AliMUONTrackExtrap::ExtrapToZCov(AliMUONTrackParam* trackParam, Double_t zEnd)
+{
+  /// Track parameters and their covariances extrapolated to the plane at "zEnd".
+  /// On return, results from the extrapolation are updated in trackParam.
+  
+  if (trackParam->GetZ() == zEnd) return; // nothing to be done if same z
+  
+  // Save the actual track parameters
+  AliMUONTrackParam trackParamSave(*trackParam);
+  Double_t nonBendingCoor        = trackParamSave.GetNonBendingCoor();
+  Double_t nonBendingSlope       = trackParamSave.GetNonBendingSlope();
+  Double_t bendingCoor                   = trackParamSave.GetBendingCoor();
+  Double_t bendingSlope          = trackParamSave.GetBendingSlope();
+  Double_t inverseBendingMomentum = trackParamSave.GetInverseBendingMomentum();
+  Double_t zBegin                = trackParamSave.GetZ();
+  
+  // Extrapolate track parameters to "zEnd"
+  ExtrapToZ(trackParam,zEnd);
+  Double_t extrapNonBendingCoor        = trackParam->GetNonBendingCoor();
+  Double_t extrapNonBendingSlope       = trackParam->GetNonBendingSlope();
+  Double_t extrapBendingCoor           = trackParam->GetBendingCoor();
+  Double_t extrapBendingSlope          = trackParam->GetBendingSlope();
+  Double_t extrapInverseBendingMomentum = trackParam->GetInverseBendingMomentum();
+  
+  // Get the pointer to the parameter covariance matrix
+  if (!trackParam->CovariancesExist()) {
+    //cout<<"W-AliMUONTrackExtrap::ExtrapToZCov: track parameter covariance matrix does not exist"<<endl;
+    //cout<<"                                    -> nothing to extrapolate !!"<<endl;
+    return;
+  }
+  TMatrixD* paramCov = trackParam->GetCovariances();
+  
+  // Calculate the jacobian related to the track parameters extrapolation to "zEnd"
+  TMatrixD jacob(5,5);
+  jacob = 0.;
+  Double_t dParam[5];
+  for (Int_t i=0; i<5; i++) {
+    // Skip jacobian calculation for parameters with no associated error
+    if ((*paramCov)(i,i) == 0.) continue;
+    // Small variation of parameter i only
+    for (Int_t j=0; j<5; j++) {
+      if (j==i) {
+        dParam[j] = TMath::Sqrt((*paramCov)(i,i));
+       if (j == 4) dParam[j] *= TMath::Sign(1.,-inverseBendingMomentum); // variation always in the same direction
+      } else dParam[j] = 0.;
+    }
+    // Set new parameters
+    trackParamSave.SetNonBendingCoor       (nonBendingCoor         + dParam[0]);
+    trackParamSave.SetNonBendingSlope      (nonBendingSlope        + dParam[1]);
+    trackParamSave.SetBendingCoor          (bendingCoor            + dParam[2]);
+    trackParamSave.SetBendingSlope         (bendingSlope           + dParam[3]);
+    trackParamSave.SetInverseBendingMomentum(inverseBendingMomentum + dParam[4]);
+    trackParamSave.SetZ                            (zBegin);
+    // Extrapolate new track parameters to "zEnd"
+    ExtrapToZ(&trackParamSave,zEnd);
+    // Calculate the jacobian
+    jacob(0,i) = (trackParamSave.GetNonBendingCoor()        - extrapNonBendingCoor        ) / dParam[i];
+    jacob(1,i) = (trackParamSave.GetNonBendingSlope()       - extrapNonBendingSlope       ) / dParam[i];
+    jacob(2,i) = (trackParamSave.GetBendingCoor()           - extrapBendingCoor           ) / dParam[i];
+    jacob(3,i) = (trackParamSave.GetBendingSlope()          - extrapBendingSlope          ) / dParam[i];
+    jacob(4,i) = (trackParamSave.GetInverseBendingMomentum() - extrapInverseBendingMomentum) / dParam[i];
+  }
+  
+  // Extrapolate track parameter covariances to "zEnd"
+  TMatrixD tmp((*paramCov),TMatrixD::kMultTranspose,jacob);
+  (*paramCov) = TMatrixD(jacob,TMatrixD::kMult,tmp);
+  
+}
+
+  //__________________________________________________________________________
 void AliMUONTrackExtrap::ExtrapToStation(AliMUONTrackParam* trackParamIn, Int_t station, AliMUONTrackParam *trackParamOut)
 {
   /// Track parameters extrapolated from "trackParamIn" to both chambers of the station(0..) "station"
@@ -150,7 +263,7 @@ void AliMUONTrackExtrap::ExtrapToStation(AliMUONTrackParam* trackParamIn, Int_t
   if ((z1 > trackParamIn->GetZ()) && (z2 > trackParamIn->GetZ())) {i1 = 0; i2 = 1;}
   else if ((z1 < trackParamIn->GetZ()) && (z2 < trackParamIn->GetZ())) {i1 = 1; i2 = 0;}
   else {
-    cout<<"E-AliMUONTrackExtrap::ExtrapToStationAliError: Starting Z ("<<trackParamIn->GetZ()
+    cout<<"E-AliMUONTrackExtrap::ExtrapToStation: Starting Z ("<<trackParamIn->GetZ()
        <<") in between z1 ("<<z1<<") and z2 ("<<z2<<") of station(0..)"<<station<<endl;
     exit(-1);
   }
@@ -167,20 +280,113 @@ void AliMUONTrackExtrap::ExtrapToStation(AliMUONTrackParam* trackParamIn, Int_t
 }
 
   //__________________________________________________________________________
+void AliMUONTrackExtrap::ExtrapToVertexUncorrected(AliMUONTrackParam* trackParam, Double_t zVtx)
+{
+  /// Extrapolation to the vertex (at the z position "zVtx") without Branson and Field correction.
+  /// Returns the track parameters resulting from the extrapolation in the current TrackParam.
+  /// Include multiple Coulomb scattering effects in trackParam covariances.
+  
+  if (trackParam->GetZ() > zVtx) { // spectro. (z<0)
+    cout<<"W-AliMUONTrackExtrap::ExtrapToVertexUncorrected: Starting Z ("<<trackParam->GetZ()
+       <<") upstream the vertex (zVtx = "<<zVtx<<")"<<endl;
+    exit(-1);
+  }
+  
+  if (zVtx < AliMUONConstants::ZAbsorberEnd()) { // spectro. (z<0)
+    cout<<"W-AliMUONTrackExtrap::ExtrapToVertexUncorrected: Ending Z ("<<zVtx
+       <<") downstream the front absorber (zAbsorberEnd = "<<AliMUONConstants::ZAbsorberEnd()<<")"<<endl;
+    
+    ExtrapToZCov(trackParam,zVtx);
+    return;
+  }
+  
+  // First Extrapolates track parameters upstream to the "Z" end of the front absorber
+  if (trackParam->GetZ() < AliMUONConstants::ZAbsorberEnd()) { // spectro. (z<0)
+    ExtrapToZCov(trackParam,AliMUONConstants::ZAbsorberEnd());
+  } else {
+    cout<<"W-AliMUONTrackExtrap::ExtrapToVertexUncorrected: Starting Z ("<<trackParam->GetZ()
+       <<") upstream or inside the front absorber (zAbsorberEnd = "<<AliMUONConstants::ZAbsorberEnd()<<")"<<endl;
+  }
+  
+  // Then go through all the absorber layers
+  Double_t tan3 = TMath::Tan(3./180.*TMath::Pi());
+  Double_t r0Norm, x0, z, zElement, dZ, nonBendingCoor, bendingCoor;
+  for (Int_t iElement=AliMUONConstants::NAbsorberElements()-1; iElement>=0; iElement--) {
+    zElement = AliMUONConstants::ZAbsorberElement(iElement);
+    z = trackParam->GetZ();
+    if (z > zElement) continue; // spectro. (z<0)
+    nonBendingCoor = trackParam->GetNonBendingCoor();
+    bendingCoor = trackParam->GetBendingCoor();
+    r0Norm = nonBendingCoor * nonBendingCoor + bendingCoor * bendingCoor;
+    r0Norm  = TMath::Sqrt(r0Norm) / TMath::Abs(trackParam->GetZ()) / tan3;
+    if (r0Norm > 1.) x0 = AliMUONConstants::X0AbsorberOut(iElement); // outer part of the absorber
+    else x0 = AliMUONConstants::X0AbsorberIn(iElement); // inner part of the absorber
+    
+    if (zVtx > zElement) {
+      ExtrapToZCov(trackParam,zElement); // extrapolate to absorber element "iElement"
+      dZ = zElement - z;
+      AddMCSEffectInTrackParamCov(trackParam,dZ,x0); // include MCS effect in covariances
+    } else {
+      ExtrapToZCov(trackParam,zVtx); // extrapolate to zVtx
+      dZ = zVtx - z;
+      AddMCSEffectInTrackParamCov(trackParam,dZ,x0); // include MCS effect in covariances
+      break;
+    }
+  }
+  
+  // finally go to the vertex
+  ExtrapToZCov(trackParam,zVtx);
+  
+}
+
+  //__________________________________________________________________________
+void AliMUONTrackExtrap::AddMCSEffectInTrackParamCov(AliMUONTrackParam *param, Double_t dZ, Double_t x0)
+{
+  /// Add to the track parameter covariances the effects of multiple Coulomb scattering
+  /// through a material of thickness "dZ" and of radiation length "x0"
+  /// assuming linear propagation and using the small angle approximation.
+  
+  Double_t bendingSlope = param->GetBendingSlope();
+  Double_t nonBendingSlope = param->GetNonBendingSlope();
+  Double_t inverseTotalMomentum2 = param->GetInverseBendingMomentum() * param->GetInverseBendingMomentum() *
+                                  (1.0 + bendingSlope * bendingSlope) /
+                                  (1.0 + bendingSlope *bendingSlope + nonBendingSlope * nonBendingSlope); 
+  // Path length in the material
+  Double_t pathLength = TMath::Abs(dZ) * TMath::Sqrt(1.0 + bendingSlope*bendingSlope + nonBendingSlope*nonBendingSlope);
+  Double_t pathLength2 = pathLength * pathLength;
+  // relativistic velocity
+  Double_t velo = 1.;
+  // Angular dispersion square of the track (variance) in a plane perpendicular to the trajectory
+  Double_t theta02 = 0.0136 / velo * (1 + 0.038 * TMath::Log(pathLength/x0));
+  theta02 *= theta02 * inverseTotalMomentum2 * pathLength / x0;
+  
+  // Add effects of multiple Coulomb scattering in track parameter covariances
+  TMatrixD* paramCov = param->GetCovariances();
+  Double_t varCoor     = pathLength2 * theta02 / 3.;
+  Double_t varSlop     = theta02;
+  Double_t covCorrSlope = pathLength * theta02 / 2.;
+  // Non bending plane
+  (*paramCov)(0,0) += varCoor;         (*paramCov)(0,1) += covCorrSlope;
+  (*paramCov)(1,0) += covCorrSlope;    (*paramCov)(1,1) += varSlop;
+  // Bending plane
+  (*paramCov)(2,2) += varCoor;         (*paramCov)(2,3) += covCorrSlope;
+  (*paramCov)(3,2) += covCorrSlope;    (*paramCov)(3,3) += varSlop;
+  
+}
+
+  //__________________________________________________________________________
 void AliMUONTrackExtrap::ExtrapToVertex(AliMUONTrackParam* trackParam, Double_t xVtx, Double_t yVtx, Double_t zVtx)
 {
   /// Extrapolation to the vertex.
   /// Returns the track parameters resulting from the extrapolation in the current TrackParam.
   /// Changes parameters according to Branson correction through the absorber 
   
-  Double_t zAbsorber = -503.0; // to be coherent with the Geant absorber geometry !!!!
-                               // spectro. (z<0) 
   // Extrapolates track parameters upstream to the "Z" end of the front absorber
-  ExtrapToZ(trackParam,zAbsorber); // !!!
+  ExtrapToZ(trackParam,AliMUONConstants::ZAbsorberEnd()); // !!!
   // Makes Branson correction (multiple scattering + energy loss)
   BransonCorrection(trackParam,xVtx,yVtx,zVtx);
   // Makes a simple magnetic field correction through the absorber
-  FieldCorrection(trackParam,zAbsorber);
+  FieldCorrection(trackParam,AliMUONConstants::ZAbsorberEnd());
 }
 
 
@@ -306,14 +512,13 @@ void AliMUONTrackExtrap::BransonCorrection(AliMUONTrackParam* trackParam, Double
   Double_t  pYZ, pX, pY, pZ, pTotal, xEndAbsorber, yEndAbsorber, radiusEndAbsorber2, pT, theta;
   Int_t sign;
   static Bool_t first = kTRUE;
-  static Double_t zBP1, zBP2, rLimit, thetaLimit, zEndAbsorber;
+  static Double_t zBP1, zBP2, rLimit, thetaLimit;
   // zBP1 for outer part and zBP2 for inner part (only at the first call)
   if (first) {
     first = kFALSE;
   
-    zEndAbsorber = -503;  // spectro (z<0)
     thetaLimit = 3.0 * (TMath::Pi()) / 180.;
-    rLimit = TMath::Abs(zEndAbsorber) * TMath::Tan(thetaLimit);
+    rLimit = TMath::Abs(AliMUONConstants::ZAbsorberEnd()) * TMath::Tan(thetaLimit);
     zBP1 = -450; // values close to those calculated with EvalAbso.C
     zBP2 = -480;
   }
@@ -335,8 +540,8 @@ void AliMUONTrackExtrap::BransonCorrection(AliMUONTrackParam* trackParam, Double
     zBP = zBP2;
   }
 
-  xBP = xEndAbsorber - (pX / pZ) * (zEndAbsorber - zBP);
-  yBP = yEndAbsorber - (pY / pZ) * (zEndAbsorber - zBP);
+  xBP = xEndAbsorber - (pX / pZ) * (AliMUONConstants::ZAbsorberEnd() - zBP);
+  yBP = yEndAbsorber - (pY / pZ) * (AliMUONConstants::ZAbsorberEnd() - zBP);
 
   // new parameters after Branson and energy loss corrections
 //   Float_t zSmear = zBP - gRandom->Gaus(0.,2.);  // !!! possible smearing of Z vertex position
index 635d481..09983d7 100644 (file)
@@ -30,15 +30,24 @@ class AliMUONTrackExtrap : public TObject
        /// set field map
   static void SetField(const AliMagF* magField) {fgkField = magField;}
   
+  static Double_t GetImpactParamFromBendingMomentum(Double_t bendingMomentum);
+  static Double_t GetBendingMomentumFromImpactParam(Double_t impactParam);
+  
   static void ExtrapToZ(AliMUONTrackParam *trackParam, Double_t Z);
+  static void ExtrapToZCov(AliMUONTrackParam* trackParam, Double_t zEnd);
   static void ExtrapToStation(AliMUONTrackParam *trackParamIn, Int_t station, AliMUONTrackParam *trackParamOut);
+  static void ExtrapToVertexUncorrected(AliMUONTrackParam* trackParam, Double_t zVtx);
   static void ExtrapToVertex(AliMUONTrackParam *trackParam, Double_t xVtx, Double_t yVtx, Double_t zVtx);
   
+  static void AddMCSEffectInTrackParamCov(AliMUONTrackParam *param, Double_t dZ, Double_t x0);
+  
   static void ExtrapOneStepRungekutta(Double_t charge, Double_t step, Double_t* vect, Double_t* vout);
   
   
  private:
-  static const AliMagF* fgkField;     //!< field map
+  static const AliMagF* fgkField;        //!< field map
+  static const Int_t    fgkMaxStepNumber; //!< Maximum number of steps for track extrapolation
+  static const Double_t fgkStepLength;   //!< Step lenght for track extrapolation
 
   // Functions
   AliMUONTrackExtrap(const AliMUONTrackExtrap& trackExtrap);
index aed91dc..afeab8f 100644 (file)
@@ -33,8 +33,8 @@
 
 #include "AliMUONTrackReconstructorK.h"
 #include "AliMagF.h"
-#include "AliMUONSegment.h"
 #include "AliMUONHitForRec.h"
+#include "AliMUONObjectPair.h"
 #include "AliMUONRawCluster.h"
 #include "AliMUONTrackParam.h"
 #include "AliMUONTrackExtrap.h"
@@ -53,7 +53,7 @@ const Double_t AliMUONTrackK::fgkChi2max = 100;
 const Int_t AliMUONTrackK::fgkTriesMax = 10000; 
 const Double_t AliMUONTrackK::fgkEpsilon = 0.002; 
 
-void mnvertLocal(Double_t* a, Int_t l, Int_t m, Int_t n, Int_t& ifail); // from AliMUONTrack
+void mnvertLocal(Double_t* a, Int_t l, Int_t m, Int_t n, Int_t& ifail);
 
 Int_t AliMUONTrackK::fgDebug = -1; //-1;
 Int_t AliMUONTrackK::fgNOfPoints = 0; 
@@ -137,10 +137,9 @@ AliMUONTrackK::AliMUONTrackK(AliMUONTrackReconstructorK *TrackReconstructor, TCl
 }
 
   //__________________________________________________________________________
-AliMUONTrackK::AliMUONTrackK(AliMUONSegment *segment)
-  //: AliMUONTrack(segment, segment, fgTrackReconstructor)
-  : AliMUONTrack(NULL, segment),
-    fStartSegment(segment),
+AliMUONTrackK::AliMUONTrackK(AliMUONObjectPair *segment)
+  : AliMUONTrack(NULL,NULL),
+    fStartSegment(new AliMUONObjectPair(*segment)),
     fPosition(0.),
     fPositionNew(0.),
     fChi2(0.),
@@ -166,20 +165,20 @@ AliMUONTrackK::AliMUONTrackK(AliMUONSegment *segment)
     fChi2Smooth(0x0)
 {
 /// Constructor from a segment
-  Double_t dX, dY, dZ;
+  Double_t dX, dY, dZ, bendingSlope, bendingImpact;
   AliMUONHitForRec *hit1, *hit2;
   AliMUONRawCluster *clus;
   TClonesArray *rawclusters;
 
   // Pointers to hits from the segment
-  hit1 = segment->GetHitForRec1();
-  hit2 = segment->GetHitForRec2();
+  hit1 = (AliMUONHitForRec*) segment->First();
+  hit2 = (AliMUONHitForRec*) segment->Second();
   hit1->SetNTrackHits(hit1->GetNTrackHits()+1); // mark hit as being on track
   hit2->SetNTrackHits(hit2->GetNTrackHits()+1); // mark hit as being on track
   // check sorting in Z
   if (TMath::Abs(hit1->GetZ()) > TMath::Abs(hit2->GetZ())) {
     hit1 = hit2;
-    hit2 = segment->GetHitForRec1();
+    hit2 = (AliMUONHitForRec*) segment->First();
   }
 
   // Fill array of track parameters
@@ -205,17 +204,21 @@ AliMUONTrackK::AliMUONTrackK(AliMUONSegment *segment)
   dZ = hit2->GetZ() - hit1->GetZ();
   dY = hit2->GetBendingCoor() - hit1->GetBendingCoor();
   dX = hit2->GetNonBendingCoor() - hit1->GetNonBendingCoor();
+  bendingSlope = (hit2->GetBendingCoor() - hit1->GetBendingCoor()) / dZ;
+  bendingImpact = hit1->GetBendingCoor() - hit1->GetZ() * bendingSlope;
   (*fTrackPar)(2,0) = TMath::ATan2(dY,dZ); // alpha
   if ((*fTrackPar)(2,0) < 0.) (*fTrackPar)(2,0) += 2*TMath::Pi(); // from 0 to 2*pi
   (*fTrackPar)(3,0) = TMath::ATan2(-dX,dZ/TMath::Cos((*fTrackPar)(2,0))); // beta
   (*fTrackPar)(2,0) -= TMath::Pi();
-  (*fTrackPar)(4,0) = 1/fgTrackReconstructor->GetBendingMomentumFromImpactParam(segment->GetBendingImpact()); // 1/Pt
+  (*fTrackPar)(4,0) = 1./AliMUONTrackExtrap::GetBendingMomentumFromImpactParam(bendingImpact); // 1/Pt
   (*fTrackPar)(4,0) *= TMath::Cos((*fTrackPar)(3,0)); // 1/p
   // Evaluate covariance (and weight) matrix
   EvalCovariance(dZ);
 
   if (fgDebug < 0 ) return;
-  cout << fgTrackReconstructor->GetNRecTracks()-1 << " " << fgTrackReconstructor->GetBendingMomentumFromImpactParam(segment->GetBendingImpact()) << " " << 1/(*fTrackPar)(4,0) << " ";
+  cout << fgTrackReconstructor->GetNRecTracks()-1 << " "
+       << AliMUONTrackExtrap::GetBendingMomentumFromImpactParam(bendingImpact)
+       << " " << 1/(*fTrackPar)(4,0) << " ";
     // from raw clusters
   for (Int_t i=0; i<2; i++) {
     hit1 = (AliMUONHitForRec*) ((*fTrackHits)[i]);
@@ -225,7 +228,7 @@ AliMUONTrackK::AliMUONTrackK(AliMUONSegment *segment)
     if (clus->GetTrack(2) != -1) cout << " " << clus->GetTrack(2);
     if (i == 0) cout << " <--> ";
   }
-  cout << " @ " << fStartSegment->GetHitForRec1()->GetChamberNumber() << endl;
+  cout << " @ " << hit1->GetChamberNumber() << endl;
   
 }
 
@@ -234,6 +237,11 @@ AliMUONTrackK::~AliMUONTrackK()
 {
 /// Destructor
 
+  if (fStartSegment) {
+    delete fStartSegment;
+    fStartSegment = 0x0;
+  }
+  
   if (fTrackHits) {
     //cout << fNmbTrackHits << endl;
     for (Int_t i = 0; i < fNmbTrackHits; i++) {
@@ -323,8 +331,11 @@ AliMUONTrackK & AliMUONTrackK::operator=(const AliMUONTrackK& source)
   // base class assignement
   //AZ TObject::operator=(source);
   AliMUONTrack::operator=(source);
-
-  fStartSegment = source.fStartSegment;
+  
+  if (fStartSegment) delete fStartSegment;
+  if (source.fStartSegment) fStartSegment = new AliMUONObjectPair(*(source.fStartSegment));
+  else fStartSegment = 0x0;
+  
   fNmbTrackHits = source.fNmbTrackHits;
   fChi2 = source.fChi2;
   fPosition = source.fPosition;
@@ -442,7 +453,7 @@ Bool_t AliMUONTrackK::KalmanFilter(Int_t ichamBeg, Int_t ichamEnd, Bool_t Back,
   } else if (fRecover) {
     hit = GetHitLastOk();
     currIndx = fTrackHits->IndexOf(hit);
-    if (currIndx < 0) hit = fStartSegment->GetHitForRec1(); // for station 3
+    if (currIndx < 0) hit = (AliMUONHitForRec*) fStartSegment->First(); // for station 3
     Back = kTRUE;
     ichamb = hit->GetChamberNumber();
     if (hit == fSkipHit || fRecover == 2 && currIndx >= 0) {
@@ -1160,7 +1171,7 @@ void AliMUONTrackK::MSThin(Int_t sign)
   momentum = 1/(*fTrackParNew)(4,0); // particle momentum
   //velo = momentum/TMath::Sqrt(momentum*momentum+muonMass*muonMass); // velocity/c for muon hypothesis
   velo = 1; // relativistic
-  path = TMath::Abs(fgTrackReconstructor->GetChamberThicknessInX0()/cosAlph/cosBeta); // path length
+  path = TMath::Abs(AliMUONConstants::ChamberThicknessInX0()/cosAlph/cosBeta); // path length
   theta0 = 0.0136/velo/momentum*TMath::Sqrt(path)*(1+0.038*TMath::Log(path)); // projected scattering angle
 
   (*fWeight)(2,2) += sign*theta0/cosBeta*theta0/cosBeta; // alpha
@@ -1191,7 +1202,7 @@ void AliMUONTrackK::StartBack(void)
   //__________________________________________________________________________
 void AliMUONTrackK::SortHits(Int_t iflag, TObjArray *array)
 {
-/// Sort hits in Z if the seed segment in the last but one station
+/// Sort hits in Z if the seed segment is in the last but one station
 /// (if iflag==0 in descending order in abs(z), if !=0 - unsort)
   
   if (iflag && ((AliMUONHitForRec*)(array->UncheckedAt(0)))->GetChamberNumber() == 6) return;
@@ -1736,7 +1747,7 @@ Bool_t AliMUONTrackK::Recover(void)
   if (chi2max < 15) imax = fNmbTrackHits - 1; // discard the last point
   // Check if the outlier is not from the seed segment
   AliMUONHitForRec *skipHit = (AliMUONHitForRec*) fTrackHits->UncheckedAt(imax);
-  if (skipHit == fStartSegment->GetHitForRec1() || skipHit == fStartSegment->GetHitForRec2()) {
+  if (skipHit == (AliMUONHitForRec*) fStartSegment->First() || skipHit == (AliMUONHitForRec*) fStartSegment->Second()) {
     //DropBranches(fStartSegment); // drop all tracks with the same seed segment
     return kFALSE; // to be changed probably
   }
@@ -1770,8 +1781,8 @@ Bool_t AliMUONTrackK::Recover(void)
   // Remove all saved steps and smoother matrices after the skipped hit 
   RemoveMatrices(skipHit->GetZ());
 
-  //AZ(z->-z) if (skipHit->GetZ() > fStartSegment->GetHitForRec2()->GetZ() || !fNSteps) {
-  if (TMath::Abs(skipHit->GetZ()) > TMath::Abs(fStartSegment->GetHitForRec2()->GetZ()) || !fNSteps) {
+  //AZ(z->-z) if (skipHit->GetZ() > ((AliMUONHitForRec*) fStartSegment->Second())->GetZ() || !fNSteps) {
+  if (TMath::Abs(skipHit->GetZ()) > TMath::Abs( ((AliMUONHitForRec*) fStartSegment->Second())->GetZ()) || !fNSteps) {
     // Propagation toward high Z or skipped hit next to segment - 
     // start track from segment 
     trackK = new ((*trackPtr)[nRecTracks]) AliMUONTrackK(fStartSegment); 
@@ -2089,7 +2100,7 @@ void AliMUONTrackK::Outlier()
   }
   // Check if the outlier is not from the seed segment
   AliMUONHitForRec *hit = (AliMUONHitForRec*) fTrackHits->UncheckedAt(imax);
-  if (hit == fStartSegment->GetHitForRec1() || hit == fStartSegment->GetHitForRec2()) return; // to be changed probably
+  if (hit == (AliMUONHitForRec*) fStartSegment->First() || hit == (AliMUONHitForRec*) fStartSegment->Second()) return; // to be changed probably
 
   // Check for missing station
   Int_t ok = 1;
@@ -2271,12 +2282,13 @@ void AliMUONTrackK::DropBranches(Int_t imax, TObjArray *hits)
 }
 
   //__________________________________________________________________________
-void AliMUONTrackK::DropBranches(AliMUONSegment *segment)
+void AliMUONTrackK::DropBranches(AliMUONObjectPair *segment)
 {
 /// Drop all candidates with the same seed segment
   Int_t nRecTracks;
   TClonesArray *trackPtr;
   AliMUONTrackK *trackK;
+  AliMUONObjectPair *trackKSegment;
 
   trackPtr = fgTrackReconstructor->GetRecTracksPtr();
   nRecTracks = fgTrackReconstructor->GetNRecTracks();
@@ -2284,9 +2296,11 @@ void AliMUONTrackK::DropBranches(AliMUONSegment *segment)
 
   for (Int_t i=icand+1; i<nRecTracks; i++) {
     trackK = (AliMUONTrackK*) ((*trackPtr)[i]);
+    trackKSegment = trackK->fStartSegment;
     if (trackK->fNmbTrackHits == 2 && trackK->GetRecover() == 0) continue;
     if (trackK->GetRecover() < 0) continue;
-    if (trackK->fStartSegment == segment) trackK->SetRecover(-1);
+    if (trackKSegment->First()  == segment->First() &&
+        trackKSegment->Second() == segment->Second()) trackK->SetRecover(-1);
   }
   if (fgDebug >= 0) cout << " Drop segment " << endl; 
 }
@@ -2304,7 +2318,7 @@ AliMUONHitForRec* AliMUONTrackK::GetHitLastOk(void)
 Int_t AliMUONTrackK::GetStation0(void)
 {
 /// Return seed station number
-  return fStartSegment->GetHitForRec1()->GetChamberNumber() / 2;
+  return ((AliMUONHitForRec*) fStartSegment->First())->GetChamberNumber() / 2;
 }
 
   //__________________________________________________________________________
@@ -2398,3 +2412,103 @@ Bool_t AliMUONTrackK::ExistDouble(void)
   delete hitArray;
   return same;
 }
+
+//______________________________________________________________________________
+ void mnvertLocal(Double_t *a, Int_t l, Int_t, Int_t n, Int_t &ifail)
+{
+///*-*-*-*-*-*-*-*-*-*-*-*Inverts a symmetric matrix*-*-*-*-*-*-*-*-*-*-*-*-*
+///*-*                    ==========================
+///*-*        inverts a symmetric matrix.   matrix is first scaled to
+///*-*        have all ones on the diagonal (equivalent to change of units)
+///*-*        but no pivoting is done since matrix is positive-definite.
+///*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
+
+  // taken from TMinuit package of Root (l>=n)
+  // fVERTs, fVERTq and fVERTpp changed to localVERTs, localVERTq and localVERTpp
+  //  Double_t localVERTs[n], localVERTq[n], localVERTpp[n];
+  Double_t * localVERTs = new Double_t[n];
+  Double_t * localVERTq = new Double_t[n];
+  Double_t * localVERTpp = new Double_t[n];
+  // fMaxint changed to localMaxint
+  Int_t localMaxint = n;
+
+    /* System generated locals */
+    Int_t aOffset;
+
+    /* Local variables */
+    Double_t si;
+    Int_t i, j, k, kp1, km1;
+
+    /* Parameter adjustments */
+    aOffset = l + 1;
+    a -= aOffset;
+
+    /* Function Body */
+    ifail = 0;
+    if (n < 1) goto L100;
+    if (n > localMaxint) goto L100;
+//*-*-                  scale matrix by sqrt of diag elements
+    for (i = 1; i <= n; ++i) {
+        si = a[i + i*l];
+        if (si <= 0) goto L100;
+        localVERTs[i-1] = 1 / TMath::Sqrt(si);
+    }
+    for (i = 1; i <= n; ++i) {
+        for (j = 1; j <= n; ++j) {
+            a[i + j*l] = a[i + j*l]*localVERTs[i-1]*localVERTs[j-1];
+        }
+    }
+//*-*-                                       . . . start main loop . . . .
+    for (i = 1; i <= n; ++i) {
+        k = i;
+//*-*-                  preparation for elimination step1
+        if (a[k + k*l] != 0) localVERTq[k-1] = 1 / a[k + k*l];
+        else goto L100;
+        localVERTpp[k-1] = 1;
+        a[k + k*l] = 0;
+        kp1 = k + 1;
+        km1 = k - 1;
+        if (km1 < 0) goto L100;
+        else if (km1 == 0) goto L50;
+        else               goto L40;
+L40:
+        for (j = 1; j <= km1; ++j) {
+            localVERTpp[j-1] = a[j + k*l];
+            localVERTq[j-1]  = a[j + k*l]*localVERTq[k-1];
+            a[j + k*l]   = 0;
+        }
+L50:
+        if (k - n < 0) goto L51;
+        else if (k - n == 0) goto L60;
+        else                goto L100;
+L51:
+        for (j = kp1; j <= n; ++j) {
+            localVERTpp[j-1] = a[k + j*l];
+            localVERTq[j-1]  = -a[k + j*l]*localVERTq[k-1];
+            a[k + j*l]   = 0;
+        }
+//*-*-                  elimination proper
+L60:
+        for (j = 1; j <= n; ++j) {
+            for (k = j; k <= n; ++k) { a[j + k*l] += localVERTpp[j-1]*localVERTq[k-1]; }
+        }
+    }
+//*-*-                  elements of left diagonal and unscaling
+    for (j = 1; j <= n; ++j) {
+        for (k = 1; k <= j; ++k) {
+            a[k + j*l] = a[k + j*l]*localVERTs[k-1]*localVERTs[j-1];
+            a[j + k*l] = a[k + j*l];
+        }
+    }
+    delete [] localVERTs;
+    delete [] localVERTq;
+    delete [] localVERTpp;
+    return;
+//*-*-                  failure return
+L100:
+    delete [] localVERTs;
+    delete [] localVERTq;
+    delete [] localVERTpp;
+    ifail = 1;
+} /* mnvertLocal */
+
index 292fa30..6dde99a 100644 (file)
@@ -20,7 +20,7 @@ class TObjArray;
 
 class AliMUONEventRecoCombi;
 class AliMUONHitForRec;
-class AliMUONSegment;
+class AliMUONObjectPair;
 class AliMUONTrackReconstructorK;
 #include "AliMUONTrack.h" 
 
@@ -32,7 +32,7 @@ class AliMUONTrackK : public AliMUONTrack {
   virtual ~AliMUONTrackK(); // Destructor
 
   AliMUONTrackK(AliMUONTrackReconstructorK *TrackReconstructor, TClonesArray *hitForRec); // Constructor
-  AliMUONTrackK(AliMUONSegment *segment); // Constructor from a segment
+  AliMUONTrackK(AliMUONObjectPair *segment); // Constructor from a segment
 
   // Pointer to hits on track
   TObjArray* GetTrackHits(void) const {return fTrackHits;} // ptr. to hits on track
@@ -47,7 +47,7 @@ class AliMUONTrackK : public AliMUONTrack {
   void SetBPFlag(Bool_t BPFlag) {fBPFlag = BPFlag;} // set backpropagation flag
   Int_t GetRecover(void) const {return fRecover;} // return recover flag 
   void SetRecover(Int_t iRecover) {fRecover = iRecover;} // set recover flag
-  AliMUONSegment* GetStartSegment(void) const {return fStartSegment;} // return seed segment
+  AliMUONObjectPair* GetStartSegment(void) const {return fStartSegment;} // return seed segment
   Bool_t KalmanFilter(Int_t ichamBeg, Int_t ichamEnd, Bool_t Back, Double_t zDipole1, Double_t zDipole2); // Kalman filter
   void StartBack(void); // start backpropagator
   void SetTrackQuality(Int_t iChi2); // compute track quality or Chi2
@@ -86,7 +86,7 @@ class AliMUONTrackK : public AliMUONTrack {
   static TClonesArray *fgHitForRec; ///< pointer to hits
   static AliMUONEventRecoCombi *fgCombi; ///< pointer to combined cluster/track finder
 
-  AliMUONSegment *fStartSegment; ///< seed segment  
+  AliMUONObjectPair *fStartSegment; ///< seed segment  
   Double_t fPosition; ///< Z-coordinate of track
   Double_t fPositionNew; //!< Z-coordinate of track
   Double_t fChi2; ///< Chi2 of track
@@ -136,7 +136,7 @@ class AliMUONTrackK : public AliMUONTrack {
   void Outlier();
   void SortHits(Int_t iflag, TObjArray *array);
   void DropBranches(Int_t imax, TObjArray *hits);
-  void DropBranches(AliMUONSegment *segment);
+  void DropBranches(AliMUONObjectPair *segment);
   Bool_t ExistDouble(AliMUONHitForRec *hit);
   Bool_t ExistDouble(void);
   void CheckBranches(TArrayD &branchChi2, Int_t nBranch);
index d92b72c..cf86ff5 100644 (file)
@@ -26,6 +26,7 @@
 ///////////////////////////////////////////////////
 
 #include <Riostream.h>
+#include <TMatrixD.h>
 
 #include "AliMUONTrackParam.h"
 #include "AliESDMuonTrack.h"
@@ -39,12 +40,13 @@ ClassImp(AliMUONTrackParam) // Class implementation in ROOT context
   //_________________________________________________________________________
 AliMUONTrackParam::AliMUONTrackParam()
   : TObject(),
-    fInverseBendingMomentum(0.),
-    fBendingSlope(0.),
+    fNonBendingCoor(0.),
     fNonBendingSlope(0.),
-    fZ(0.),
     fBendingCoor(0.),
-    fNonBendingCoor(0.),
+    fBendingSlope(0.),
+    fInverseBendingMomentum(0.),
+    fZ(0.),
+    fCovariances(0x0),
     fHitForRecPtr(0x0)
 {
   /// Constructor
@@ -53,15 +55,17 @@ AliMUONTrackParam::AliMUONTrackParam()
   //_________________________________________________________________________
 AliMUONTrackParam::AliMUONTrackParam(const AliMUONTrackParam& theMUONTrackParam)
   : TObject(theMUONTrackParam),
-    fInverseBendingMomentum(theMUONTrackParam.fInverseBendingMomentum), 
-    fBendingSlope(theMUONTrackParam.fBendingSlope),
+    fNonBendingCoor(theMUONTrackParam.fNonBendingCoor),
     fNonBendingSlope(theMUONTrackParam.fNonBendingSlope),
-    fZ(theMUONTrackParam.fZ),
     fBendingCoor(theMUONTrackParam.fBendingCoor),
-    fNonBendingCoor(theMUONTrackParam.fNonBendingCoor),
+    fBendingSlope(theMUONTrackParam.fBendingSlope),
+    fInverseBendingMomentum(theMUONTrackParam.fInverseBendingMomentum), 
+    fZ(theMUONTrackParam.fZ),
+    fCovariances(0x0),
     fHitForRecPtr(theMUONTrackParam.fHitForRecPtr)
 {
   /// Copy constructor
+  if (theMUONTrackParam.fCovariances) fCovariances = new TMatrixD(*(theMUONTrackParam.fCovariances));
 }
 
   //_________________________________________________________________________
@@ -74,14 +78,21 @@ AliMUONTrackParam& AliMUONTrackParam::operator=(const AliMUONTrackParam& theMUON
   // base class assignement
   TObject::operator=(theMUONTrackParam);
 
-  fInverseBendingMomentum =  theMUONTrackParam.fInverseBendingMomentum; 
-  fBendingSlope           =  theMUONTrackParam.fBendingSlope; 
-  fNonBendingSlope        =  theMUONTrackParam.fNonBendingSlope; 
-  fZ                      =  theMUONTrackParam.fZ; 
-  fBendingCoor            =  theMUONTrackParam.fBendingCoor; 
-  fNonBendingCoor         =  theMUONTrackParam.fNonBendingCoor;
-  fHitForRecPtr           =  theMUONTrackParam.fHitForRecPtr;
-
+  fNonBendingCoor              =  theMUONTrackParam.fNonBendingCoor;
+  fNonBendingSlope             =  theMUONTrackParam.fNonBendingSlope; 
+  fBendingCoor                 =  theMUONTrackParam.fBendingCoor; 
+  fBendingSlope                =  theMUONTrackParam.fBendingSlope; 
+  fInverseBendingMomentum      =  theMUONTrackParam.fInverseBendingMomentum; 
+  fZ                           =  theMUONTrackParam.fZ; 
+  
+  if (theMUONTrackParam.fCovariances) {
+    if (fCovariances) *fCovariances = *(theMUONTrackParam.fCovariances);
+    else fCovariances = new TMatrixD(*(theMUONTrackParam.fCovariances));
+  } else if (fCovariances) {
+    delete fCovariances;
+    fCovariances = 0x0;
+  }
+  
   return *this;
 }
 
@@ -91,19 +102,19 @@ AliMUONTrackParam::~AliMUONTrackParam()
 /// Destructor
 /// Update the number of TrackHit's connected to the attached HitForRec if any
   if (fHitForRecPtr) fHitForRecPtr->SetNTrackHits(fHitForRecPtr->GetNTrackHits() - 1); // decrement NTrackHits of hit
+  DeleteCovariances();
 }
 
   //__________________________________________________________________________
 void AliMUONTrackParam::SetTrackParam(AliMUONTrackParam& theMUONTrackParam)
 {
-  /// Set track parameters from "TrackParam" leaving pointer to fHitForRecPtr unchanged
-  fInverseBendingMomentum =  theMUONTrackParam.fInverseBendingMomentum; 
-  fBendingSlope           =  theMUONTrackParam.fBendingSlope; 
-  fNonBendingSlope        =  theMUONTrackParam.fNonBendingSlope; 
-  fZ                      =  theMUONTrackParam.fZ; 
-  fBendingCoor            =  theMUONTrackParam.fBendingCoor; 
-  fNonBendingCoor         =  theMUONTrackParam.fNonBendingCoor;
-  
+  /// Set track parameters from "TrackParam" leaving pointer to fHitForRecPtr and parameter covariances unchanged
+  fNonBendingCoor              =  theMUONTrackParam.fNonBendingCoor;
+  fNonBendingSlope             =  theMUONTrackParam.fNonBendingSlope; 
+  fBendingCoor                 =  theMUONTrackParam.fBendingCoor; 
+  fBendingSlope                =  theMUONTrackParam.fBendingSlope; 
+  fInverseBendingMomentum      =  theMUONTrackParam.fInverseBendingMomentum; 
+  fZ                           =  theMUONTrackParam.fZ; 
 }
 
   //__________________________________________________________________________
@@ -111,31 +122,28 @@ AliMUONHitForRec* AliMUONTrackParam::GetHitForRecPtr(void) const
 {
 /// return pointer to HitForRec attached to the current TrackParam
 /// this method should not be called when fHitForRecPtr == NULL
-  if (!fHitForRecPtr) AliWarning("AliMUONTrackParam::GetHitForRecPtr: fHitForRecPtr == NULL");
+  if (!fHitForRecPtr) AliWarning("fHitForRecPtr == NULL");
   return fHitForRecPtr;
 }
 
   //__________________________________________________________________________
-Int_t AliMUONTrackParam::Compare(const TObject* TrackParam) const
+void AliMUONTrackParam::SetHitForRecPtr(AliMUONHitForRec* hitForRec)
 {
-/// "Compare" function to sort with decreasing Z (spectro. muon Z <0).
-/// Returns 1 (0, -1) if Z of current TrackHit
-/// is smaller than (equal to, larger than) Z of TrackHit
-  if (fHitForRecPtr->GetZ() < ((AliMUONTrackParam*)TrackParam)->fHitForRecPtr->GetZ()) return(1);
-  else if (fHitForRecPtr->GetZ() == ((AliMUONTrackParam*)TrackParam)->fHitForRecPtr->GetZ()) return(0);
-  else return(-1);
+/// set pointeur to associated HitForRec and update the number of TrackHit's connected to it
+  fHitForRecPtr = hitForRec;
+  fHitForRecPtr->SetNTrackHits(fHitForRecPtr->GetNTrackHits() + 1); // increment NTrackHits of hit
 }
 
   //_________________________________________________________________________
 void AliMUONTrackParam::GetParamFrom(const AliESDMuonTrack& esdMuonTrack)
 {
   /// assigned value form ESD track.
-  fInverseBendingMomentum =  esdMuonTrack.GetInverseBendingMomentum();
-  fBendingSlope           =  TMath::Tan(esdMuonTrack.GetThetaY());
-  fNonBendingSlope        =  TMath::Tan(esdMuonTrack.GetThetaX());
-  fZ                      =  esdMuonTrack.GetZ(); 
-  fBendingCoor            =  esdMuonTrack.GetBendingCoor(); 
-  fNonBendingCoor         =  esdMuonTrack.GetNonBendingCoor();
+  fInverseBendingMomentum      =  esdMuonTrack.GetInverseBendingMomentum();
+  fBendingSlope                =  TMath::Tan(esdMuonTrack.GetThetaY());
+  fNonBendingSlope             =  TMath::Tan(esdMuonTrack.GetThetaX());
+  fZ                           =  esdMuonTrack.GetZ(); 
+  fBendingCoor                 =  esdMuonTrack.GetBendingCoor(); 
+  fNonBendingCoor              =  esdMuonTrack.GetNonBendingCoor();
 }
 
   //_________________________________________________________________________
@@ -203,6 +211,97 @@ Double_t AliMUONTrackParam::P() const
   
 }
 
+  //__________________________________________________________________________
+TMatrixD* AliMUONTrackParam::GetCovariances()
+{
+  /// Return the covariance matrix (create it before if needed)
+  if (!fCovariances) {
+    fCovariances = new TMatrixD(5,5);
+    (*fCovariances) = 0;
+  }
+  return fCovariances;
+  }
+
+  //__________________________________________________________________________
+void AliMUONTrackParam::SetCovariances(TMatrixD* covariances)
+{
+  /// Set the covariance matrix
+  if (covariances == fCovariances) return; // nothing to be done
+  if (fCovariances) *fCovariances = *covariances;
+  else fCovariances = new TMatrixD(*covariances);
+}
+
+  //__________________________________________________________________________
+void AliMUONTrackParam::SetCovariances(Double_t matrix[5][5])
+{
+  /// Set the covariance matrix
+  if (fCovariances) fCovariances->SetMatrixArray(&(matrix[0][0]));
+  else fCovariances = new TMatrixD(5,5,&(matrix[0][0]));
+}
+
+  //__________________________________________________________________________
+void AliMUONTrackParam::SetVariances(Double_t matrix[5][5])
+{
+  /// Set the diagonal terms of the covariance matrix (variances)
+  if (!fCovariances) fCovariances = new TMatrixD(5,5);
+  (*fCovariances) = 0;
+  for (Int_t i=0; i<5; i++) (*fCovariances)(i,i) = matrix[i][i];
+}
+
+  //__________________________________________________________________________
+void AliMUONTrackParam::DeleteCovariances()
+{
+  /// Delete the covariance matrix
+  if (fCovariances) delete fCovariances;
+  fCovariances = 0x0;
+}
+
+  //__________________________________________________________________________
+void AliMUONTrackParam::EvalCovariances(AliMUONHitForRec* hit2)
+{
+  /// Evaluate covariances assuming the track is only a straight line
+  /// between the HitForRec attached to the current TrackParam and hit2.
+  /// Nothing can be done on fInverseBendingMomentum (-> 50% err).
+  
+  // Allocate memory if needed
+  if (!fCovariances) fCovariances = new TMatrixD(5,5);
+  
+  // Reset the covariance matrix
+  (*fCovariances) = 0;
+  
+  if (!fHitForRecPtr) {
+    AliWarning("fHitForRecPtr == NULL: cannot calculate TrackParam covariances");
+    return;
+  }
+  
+  Double_t dz = fHitForRecPtr->GetZ() - hit2->GetZ();
+  
+  // Non bending plane
+  (*fCovariances)(0,0) = fHitForRecPtr->GetNonBendingReso2();
+  (*fCovariances)(0,1) = fHitForRecPtr->GetNonBendingReso2() / dz;
+  (*fCovariances)(1,0) = (*fCovariances)(0,1);
+  (*fCovariances)(1,1) = ( fHitForRecPtr->GetNonBendingReso2() + hit2->GetNonBendingReso2() ) / dz / dz;
+  // Bending plane
+  (*fCovariances)(2,2) = fHitForRecPtr->GetBendingReso2();
+  (*fCovariances)(2,3) = fHitForRecPtr->GetBendingReso2() / dz;
+  (*fCovariances)(3,2) = (*fCovariances)(2,3);
+  (*fCovariances)(3,3) = ( fHitForRecPtr->GetBendingReso2() + hit2->GetBendingReso2() ) / dz / dz;
+  // Inverse bending momentum
+  (*fCovariances)(4,4) = 0.5*fInverseBendingMomentum * 0.5*fInverseBendingMomentum; // error 50%
+  
+}
+
+  //__________________________________________________________________________
+Int_t AliMUONTrackParam::Compare(const TObject* trackParam) const
+{
+  /// "Compare" function to sort with decreasing Z (spectro. muon Z <0).
+  /// Returns 1 (0, -1) if Z of current TrackHit
+  /// is smaller than (equal to, larger than) Z of TrackHit
+  if (fHitForRecPtr->GetZ() < ((AliMUONTrackParam*)trackParam)->fHitForRecPtr->GetZ()) return(1);
+  else if (fHitForRecPtr->GetZ() == ((AliMUONTrackParam*)trackParam)->fHitForRecPtr->GetZ()) return(0);
+  else return(-1);
+}
+
 //_____________________________________________-
 void AliMUONTrackParam::Print(Option_t* opt) const
 {
index fb71200..9b8ad7a 100644 (file)
@@ -15,6 +15,7 @@
 ////////////////////////////////////////////////////
 
 #include <TObject.h>
+#include <TMatrixDfwd.h>
 #include "AliMUONHitForRec.h"
 
 class AliESDMuonTrack;
@@ -35,53 +36,72 @@ class AliMUONTrackParam : public TObject
        /// return inverse bending momentum (GeV/c ** -1) times the charge (assumed forward motion)
   Double_t GetInverseBendingMomentum(void) const {return fInverseBendingMomentum;}
        /// set inverse bending momentum (GeV/c ** -1) times the charge (assumed forward motion)
-  void     SetInverseBendingMomentum(Double_t InverseBendingMomentum) {fInverseBendingMomentum = InverseBendingMomentum;}
+  void     SetInverseBendingMomentum(Double_t inverseBendingMomentum) {fInverseBendingMomentum = inverseBendingMomentum;}
        /// return bending slope (cm ** -1)
   Double_t GetBendingSlope(void) const {return fBendingSlope;}
        /// set bending slope (cm ** -1)
-  void     SetBendingSlope(Double_t BendingSlope) {fBendingSlope = BendingSlope;}
+  void     SetBendingSlope(Double_t bendingSlope) {fBendingSlope = bendingSlope;}
        /// return non bending slope (cm ** -1)
   Double_t GetNonBendingSlope(void) const {return fNonBendingSlope;}
        /// set non bending slope (cm ** -1)
-  void     SetNonBendingSlope(Double_t NonBendingSlope) {fNonBendingSlope = NonBendingSlope;}
+  void     SetNonBendingSlope(Double_t nonBendingSlope) {fNonBendingSlope = nonBendingSlope;}
        /// return Z coordinate (cm)
   Double_t GetZ(void) const {return fZ;}
        /// set Z coordinate (cm)
-  void     SetZ(Double_t Z) {fZ = Z;}
+  void     SetZ(Double_t z) {fZ = z;}
        /// return bending coordinate (cm)
   Double_t GetBendingCoor(void) const {return fBendingCoor;}
        /// set bending coordinate (cm)
-  void     SetBendingCoor(Double_t BendingCoor) {fBendingCoor = BendingCoor;}
+  void     SetBendingCoor(Double_t bendingCoor) {fBendingCoor = bendingCoor;}
        /// return non bending coordinate (cm)
   Double_t GetNonBendingCoor(void) const {return fNonBendingCoor;}
        /// set non bending coordinate (cm)
-  void     SetNonBendingCoor(Double_t NonBendingCoor) {fNonBendingCoor = NonBendingCoor;}
-  void              SetTrackParam(AliMUONTrackParam& TrackParam);
+  void     SetNonBendingCoor(Double_t nonBendingCoor) {fNonBendingCoor = nonBendingCoor;}
+  
+  void     SetTrackParam(AliMUONTrackParam& theMUONTrackParam);
+  
   AliMUONHitForRec* GetHitForRecPtr(void) const;
-       /// set pointeur to associated HitForRec if any
-  void              SetHitForRecPtr(AliMUONHitForRec* HitForRec) {fHitForRecPtr = HitForRec;}
+  void     SetHitForRecPtr(AliMUONHitForRec* hitForRec);
   
   Double_t Px() const;  // return px
   Double_t Py() const;  // return py
   Double_t Pz() const;  // return pz
   Double_t P()  const;  // return total momentum
 
+       /// return kTRUE if the covariance matrix exist, kFALSE if not
+  Bool_t    CovariancesExist(void) {return (fCovariances) ? kTRUE : kFALSE;}
+  TMatrixD* GetCovariances(void);
+  void      SetCovariances(TMatrixD* covariances);
+  void      SetCovariances(Double_t matrix[5][5]);
+  void      SetVariances(Double_t matrix[5][5]);
+  void      DeleteCovariances(void);
+  
+  void EvalCovariances(AliMUONHitForRec* hit2);
+
        /// necessary for sorting TClonesArray of TrackHit's
   Bool_t IsSortable () const {return kTRUE;}
-       /// "Compare" function for sorting
-  Int_t Compare(const TObject* TrackParam) const;
+  Int_t Compare(const TObject* trackParam) const;
 
   virtual void Print(Option_t* opt="") const;
  
 
  private:
-  Double_t fInverseBendingMomentum; ///< Inverse bending momentum (GeV/c ** -1) times the charge (assumed forward motion)
-  Double_t fBendingSlope; ///< Bending slope (cm ** -1)
+  // Parameters
+  Double_t fNonBendingCoor; ///< Non bending coordinate (cm)
   Double_t fNonBendingSlope; ///< Non bending slope (cm ** -1)
+  Double_t fBendingCoor; ///< Bending coordinate (cm)
+  Double_t fBendingSlope; ///< Bending slope (cm ** -1)
+  Double_t fInverseBendingMomentum; ///< Inverse bending momentum (GeV/c ** -1) times the charge (assumed forward motion)
   Double_t fZ; ///< Z coordinate (cm)
-  Double_t fBendingCoor; ///< bending coordinate (cm)
-  Double_t fNonBendingCoor; ///< non bending coordinate (cm)
-
+  
+  /// Covariance matrix of track parameters, ordered as follow:
+  ///    <X,X>      <X,SlopeX>        <X,Y>      <X,SlopeY>       <X,InvP_yz>
+  /// <X,SlopeX>  <SlopeX,SlopeX>  <Y,SlopeX>  <SlopeX,SlopeY>  <SlopeX,InvP_yz>
+  ///    <X,Y>      <Y,SlopeX>        <Y,Y>      <Y,SlopeY>       <Y,InvP_yz>
+  /// <X,SlopeY>  <SlopeX,SlopeY>  <Y,SlopeY>  <SlopeY,SlopeY>  <SlopeY,InvP_yz>
+  /// <X,InvP_yz> <SlopeX,InvP_yz> <Y,InvP_yz> <SlopeY,InvP_yz> <InvP_yz,InvP_yz>
+  TMatrixD *fCovariances; //!
+  
   AliMUONHitForRec *fHitForRecPtr; //!< Pointer to associated HitForRec if any
   
   ClassDef(AliMUONTrackParam, 2) // Track parameters in ALICE dimuon spectrometer
index 9a90ed7..e56892d 100644 (file)
 #include "AliMUONConstants.h"
 #include "AliMUONRawCluster.h"
 #include "AliMUONHitForRec.h"
-#include "AliMUONSegment.h"
+#include "AliMUONObjectPair.h"
 #include "AliMUONTrack.h"
 #include "AliMUONTrackParam.h"
 #include "AliMUONTrackExtrap.h"
 #include "AliLog.h"
-#include <TVirtualFitter.h>
+#include <TMinuit.h>
 
 // Functions to be minimized with Minuit
 void TrackChi2(Int_t &NParam, Double_t *Gradient, Double_t &Chi2, Double_t *Param, Int_t Flag);
 void TrackChi2MCS(Int_t &NParam, Double_t *Gradient, Double_t &Chi2, Double_t *Param, Int_t Flag);
 
-void mnvertLocal(Double_t* a, Int_t l, Int_t m, Int_t n, Int_t& ifail);
-
 Double_t MultipleScatteringAngle2(AliMUONTrackParam *param);
 
 /// \cond CLASSIMP
@@ -55,14 +53,12 @@ ClassImp(AliMUONTrackReconstructor) // Class implementation in ROOT context
 /// \endcond
 
 //************* Defaults parameters for reconstruction
-const Double_t AliMUONTrackReconstructor::fgkDefaultMaxChi2 = 100.0;
-
-TVirtualFitter* AliMUONTrackReconstructor::fgFitter = 0x0; 
+const Double_t AliMUONTrackReconstructor::fgkMaxNormChi2 = 100.0;
+const Bool_t AliMUONTrackReconstructor::fgkTrackAllTracks = kFALSE;
 
 //__________________________________________________________________________
 AliMUONTrackReconstructor::AliMUONTrackReconstructor(AliMUONData* data)
-  : AliMUONVTrackReconstructor(data),
-    fMaxChi2(fgkDefaultMaxChi2)
+  : AliMUONVTrackReconstructor(data)
 {
   /// Constructor for class AliMUONTrackReconstructor
   
@@ -84,7 +80,7 @@ void AliMUONTrackReconstructor::AddHitsForRecFromRawClusters()
   TTree *treeR;
   AliMUONHitForRec *hitForRec;
   AliMUONRawCluster *clus;
-  Int_t iclus, nclus, nTRentries;
+  Int_t iclus, nclus;
   TClonesArray *rawclusters;
   AliDebug(1,"Enter AddHitsForRecFromRawClusters");
   
@@ -94,23 +90,11 @@ void AliMUONTrackReconstructor::AddHitsForRecFromRawClusters()
     exit(0);
   }
   
-  nTRentries = Int_t(treeR->GetEntries());
-  
-  if (!(fMUONData->IsRawClusterBranchesInTree())) {
-    AliError(Form("RawCluster information is not avalaible, nTRentries = %d not equal to 1",nTRentries));
-    exit(0);
-  }
-
   fMUONData->SetTreeAddress("RC");
   fMUONData->GetRawClusters(); // only one entry  
   
   // Loop over tracking chambers
   for (Int_t ch = 0; ch < AliMUONConstants::NTrackingCh(); ch++) {
-    // number of HitsForRec to 0 for the chamber
-    fNHitsForRecPerChamber[ch] = 0;
-    // index of first HitForRec for the chamber
-    if (ch == 0) fIndexOfFirstHitForRecPerChamber[ch] = 0;
-    else fIndexOfFirstHitForRecPerChamber[ch] = fNHitsForRec;
     rawclusters =fMUONData->RawClusters(ch);
     nclus = (Int_t) (rawclusters->GetEntries());
     // Loop over (cathode correlated) raw clusters
@@ -120,7 +104,6 @@ void AliMUONTrackReconstructor::AddHitsForRecFromRawClusters()
       // and increment number of AliMUONHitForRec's (total and in chamber)
       hitForRec = new ((*fHitsForRecPtr)[fNHitsForRec]) AliMUONHitForRec(clus);
       fNHitsForRec++;
-      (fNHitsForRecPerChamber[ch])++;
       // more information into HitForRec
       hitForRec->SetBendingReso2(clus->GetErrY() * clus->GetErrY());
       hitForRec->SetNonBendingReso2(clus->GetErrX() * clus->GetErrX());
@@ -129,7 +112,6 @@ void AliMUONTrackReconstructor::AddHitsForRecFromRawClusters()
       hitForRec->SetHitNumber(iclus);
       // Z coordinate of the raw cluster (cm)
       hitForRec->SetZ(clus->GetZ(0));
-      
       StdoutToAliDebug(3,
                        cout << "Chamber " << ch <<
                        " raw cluster  " << iclus << " : " << endl;
@@ -161,31 +143,6 @@ void AliMUONTrackReconstructor::AddHitsForRecFromRawClusters()
 }
 
   //__________________________________________________________________________
-void AliMUONTrackReconstructor::MakeSegments(void)
-{
-  /// To make the list of segments in all stations,
-  /// from the list of hits to be reconstructed
-  AliDebug(1,"Enter MakeSegments");
-  // Loop over stations
-  for (Int_t st = 0; st < AliMUONConstants::NTrackingCh()/2; st++) MakeSegmentsPerStation(st); 
-  
-  StdoutToAliDebug(3,
-    cout << "end of MakeSegments" << endl;
-    for (Int_t st = 0; st < AliMUONConstants::NTrackingCh()/2; st++) 
-    {
-      cout << "station " << st
-           << "  has " << fNSegments[st] << " segments:"
-           << endl;
-      for (Int_t seg = 0; seg < fNSegments[st]; seg++) 
-      {
-             ((*fSegmentsPtr[st])[seg])->Print();
-      }
-    }
-                   );
-  return;
-}
-
-  //__________________________________________________________________________
 void AliMUONTrackReconstructor::MakeTracks(void)
 {
   /// To make the tracks from the list of segments and points in all stations
@@ -196,489 +153,555 @@ void AliMUONTrackReconstructor::MakeTracks(void)
   FollowTracks();
   // Remove double tracks
   RemoveDoubleTracks();
-  UpdateHitForRecAtHit();
+  // Propagate tracks to the vertex through absorber
+  ExtrapTracksToVertex();
+  // Fill out the AliMUONTrack's
+  FillMUONTrack();
 }
 
   //__________________________________________________________________________
 void AliMUONTrackReconstructor::MakeTrackCandidates(void)
 {
-  /// To make track candidates
-  /// with at least 3 aligned points in stations(1..) 4 and 5
-  /// (two Segment's or one Segment and one HitForRec)
-  Int_t begStation, iBegSegment, nbCan1Seg1Hit, nbCan2Seg;
-  AliMUONSegment *begSegment;
+  /// To make track candidates:
+  /// Start with segments station(1..) 4 or 5 then follow track in station 5 or 4.
+  /// Good candidates are made of at least three hitForRec's.
+  /// Keep only best candidates or all of them according to the flag fgkTrackAllTracks.
+  TClonesArray *segments;
+  AliMUONObjectPair *segment;
+  AliMUONHitForRec *hitForRec1, *hitForRec2;
+  AliMUONTrack *track;
+  AliMUONTrackParam *trackParamAtFirstHit;
+  Int_t iCandidate = 0;
+
   AliDebug(1,"Enter MakeTrackCandidates");
-  // Loop over stations(1..) 5 and 4 for the beginning segment
-  for (begStation = 4; begStation > 2; begStation--) {
-    // Loop over segments in the beginning station
-    for (iBegSegment = 0; iBegSegment < fNSegments[begStation]; iBegSegment++) {
-      // pointer to segment
-      begSegment = (AliMUONSegment*) ((*fSegmentsPtr[begStation])[iBegSegment]);
-      AliDebug(2,Form("Look for TrackCandidate's with Segment %d  in Station(0..) %d", iBegSegment, begStation));
-      // Look for track candidates with two segments,
-      // "begSegment" and all compatible segments in other station.
-      // Only for beginning station(1..) 5
-      // because candidates with 2 segments have to looked for only once.
-      if (begStation == 4)
-       nbCan2Seg = MakeTrackCandidatesWithTwoSegments(begSegment);
-      // Look for track candidates with one segment and one point,
-      // "begSegment" and all compatible HitForRec's in other station.
-      // Only if "begSegment" does not belong already to a track candidate.
-      // Is that a too strong condition ????
-      if (!(begSegment->GetInTrack()))
-       nbCan1Seg1Hit = MakeTrackCandidatesWithOneSegmentAndOnePoint(begSegment);
-    } // for (iBegSegment = 0;...
-  } // for (begStation = 4;...
-  return;
-}
 
-  //__________________________________________________________________________
-Int_t AliMUONTrackReconstructor::MakeTrackCandidatesWithTwoSegments(AliMUONSegment *BegSegment)
-{
-  /// To make track candidates with two segments in stations(1..) 4 and 5,
-  /// the first segment being pointed to by "BegSegment".
-  /// Returns the number of such track candidates.
-  Int_t endStation, iEndSegment, nbCan2Seg;
-  AliMUONSegment *endSegment;
-  AliMUONSegment *extrapSegment = NULL;
-  AliMUONTrack *recTrack;
-  Double_t mcsFactor;
-  AliDebug(1,"Enter MakeTrackCandidatesWithTwoSegments");
-  // Station for the end segment
-  endStation = 7 - (BegSegment->GetHitForRec1())->GetChamberNumber() / 2;
-  // multiple scattering factor corresponding to one chamber
-  mcsFactor = 0.0136 /
-    GetBendingMomentumFromImpactParam(BegSegment->GetBendingImpact());
-  mcsFactor    = fChamberThicknessInX0 * mcsFactor * mcsFactor;
-  // linear extrapolation to end station
-  // number of candidates with 2 segments to 0
-  nbCan2Seg = 0;
-  // Loop over segments in the end station
-  for (iEndSegment = 0; iEndSegment < fNSegments[endStation]; iEndSegment++) {
-    // pointer to segment
-    endSegment = (AliMUONSegment*) ((*fSegmentsPtr[endStation])[iEndSegment]);
-    // test compatibility between current segment and "extrapSegment"
-    // 4 because 4 quantities in chi2
-    extrapSegment =
-      BegSegment->CreateSegmentFromLinearExtrapToStation(endSegment->GetZ(), mcsFactor);
-    if ((endSegment->
-        NormalizedChi2WithSegment(extrapSegment,
-                                  fMaxSigma2Distance)) <= 4.0) {
-      // both segments compatible:
-      // make track candidate from "begSegment" and "endSegment"
-      AliDebug(2,Form("TrackCandidate with Segment %d in Station(0..) %d", iEndSegment, endStation));
-      // flag for both segments in one track:
-      // to be done in track constructor ????
-      BegSegment->SetInTrack(kTRUE);
-      endSegment->SetInTrack(kTRUE);
-      recTrack = new ((*fRecTracksPtr)[fNRecTracks]) AliMUONTrack(BegSegment, endSegment);
-      // Set track parameters at vertex from last stations 4 & 5
-      CalcTrackParamAtVertex(recTrack);
+  // Loop over stations(1..) 5 and 4 and make track candidates
+  for (Int_t istat=4; istat>=3; istat--) {
+    // Make segments in the station
+    segments = MakeSegmentsInStation(istat);
+    // Loop over segments
+    for (Int_t iseg=0; iseg<segments->GetEntriesFast(); iseg++) {
+      AliDebug(1,Form("Making primary candidate(1..) %d",++iCandidate));
+      // Transform segments to tracks and put them at the end of fRecTracksPtr
+      segment = (AliMUONObjectPair*) ((*segments)[iseg]);
+      hitForRec1 = (AliMUONHitForRec*) segment->First();
+      hitForRec2 = (AliMUONHitForRec*) segment->Second();
+      track = new ((*fRecTracksPtr)[fRecTracksPtr->GetLast()+1]) AliMUONTrack(hitForRec1, hitForRec2);
       fNRecTracks++;
-      if (AliLog::GetGlobalDebugLevel() > 1) recTrack->RecursiveDump();
-      // increment number of track candidates with 2 segments
-      nbCan2Seg++;
+      // Add MCS effects in parameter covariances
+      trackParamAtFirstHit = (AliMUONTrackParam*) (track->GetTrackParamAtHit()->First());
+      AliMUONTrackExtrap::AddMCSEffectInTrackParamCov(trackParamAtFirstHit,AliMUONConstants::ChamberThicknessInX0(),1.);
+      // Printout for debuging
+      if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 2) || (AliLog::GetGlobalDebugLevel() >= 2)) {
+        cout<<endl<<"Track parameter covariances at first hit with multiple Coulomb scattering effects:"<<endl;
+        trackParamAtFirstHit->GetCovariances()->Print();
+      }
+      // Look for compatible hitForRec(s) in the other station
+      FollowTrackInStation(track,7-istat);
     }
-    delete extrapSegment; // should not delete HitForRec's it points to !!!!
-  } // for (iEndSegment = 0;...
-  return nbCan2Seg;
+    // delete the array of segments
+    delete segments;
+  }
+  fRecTracksPtr->Compress(); // this is essential before checking tracks
+  
+  // Keep all different tracks or only the best ones as required
+  if (fgkTrackAllTracks) RemoveIdenticalTracks();
+  else RemoveDoubleTracks();
+  
+  AliDebug(1,Form("Number of good candidates = %d",fNRecTracks));
+  
 }
 
   //__________________________________________________________________________
-Int_t AliMUONTrackReconstructor::MakeTrackCandidatesWithOneSegmentAndOnePoint(AliMUONSegment *BegSegment)
+void AliMUONTrackReconstructor::RemoveIdenticalTracks(void)
 {
-  /// To make track candidates with one segment and one point
-  /// in stations(1..) 4 and 5,
-  /// the segment being pointed to by "BegSegment".
-  Int_t ch, ch1, ch2, endStation, iHit, iHitMax, iHitMin, nbCan1Seg1Hit;
-  AliMUONHitForRec *extrapHitForRec= NULL;
-  AliMUONHitForRec *hit;
-  AliMUONTrack *recTrack;
-  Double_t mcsFactor;
-  AliDebug(1,"Enter MakeTrackCandidatesWithOneSegmentAndOnePoint");
-  // station for the end point
-  endStation = 7 - (BegSegment->GetHitForRec1())->GetChamberNumber() / 2;
-  // multiple scattering factor corresponding to one chamber
-  mcsFactor = 0.0136 /
-    GetBendingMomentumFromImpactParam(BegSegment->GetBendingImpact());
-  mcsFactor    = fChamberThicknessInX0 * mcsFactor * mcsFactor;
-  // first and second chambers(0..) in the end station
-  ch1 = 2 * endStation;
-  ch2 = ch1 + 1;
-  // number of candidates to 0
-  nbCan1Seg1Hit = 0;
-  // Loop over chambers of the end station
-  for (ch = ch2; ch >= ch1; ch--) {
-    // limits for the hit index in the loop
-    iHitMin = fIndexOfFirstHitForRecPerChamber[ch];
-    iHitMax = iHitMin + fNHitsForRecPerChamber[ch];
-    // Loop over HitForRec's in the chamber
-    for (iHit = iHitMin; iHit < iHitMax; iHit++) {
-      // pointer to HitForRec
-      hit = (AliMUONHitForRec*) ((*fHitsForRecPtr)[iHit]);
-      // test compatibility between current HitForRec and "extrapHitForRec"
-      // 2 because 2 quantities in chi2
-      // linear extrapolation to chamber
-      extrapHitForRec =
-       BegSegment->CreateHitForRecFromLinearExtrapToChamber( hit->GetZ(), mcsFactor);
-      if ((hit->
-          NormalizedChi2WithHitForRec(extrapHitForRec,
-                                      fMaxSigma2Distance)) <= 2.0) {
-       // both HitForRec's compatible:
-       // make track candidate from begSegment and current HitForRec
-       AliDebug(2, Form("TrackCandidate with HitForRec  %d in Chamber(0..) %d", iHit, ch));
-       // flag for beginning segments in one track:
-       // to be done in track constructor ????
-       BegSegment->SetInTrack(kTRUE);
-       recTrack = new ((*fRecTracksPtr)[fNRecTracks]) AliMUONTrack(BegSegment, hit);
-        // Set track parameters at vertex from last stations 4 & 5
-        CalcTrackParamAtVertex(recTrack);
-       // the right place to eliminate "double counting" ???? how ????
-       fNRecTracks++;
-       if (AliLog::GetGlobalDebugLevel() > 1) recTrack->RecursiveDump();
-       // increment number of track candidates
-       nbCan1Seg1Hit++;
-      }
-      delete extrapHitForRec;
-    } // for (iHit = iHitMin;...
-  } // for (ch = ch2;...
-  return nbCan1Seg1Hit;
+  /// To remove identical tracks:
+  /// Tracks are considered identical if they have all their hits in common.
+  /// One keeps the track with the larger number of hits if need be
+  AliMUONTrack *track1, *track2, *trackToRemove;
+  Int_t hitsInCommon, nHits1, nHits2;
+  Bool_t removedTrack1;
+  // Loop over first track of the pair
+  track1 = (AliMUONTrack*) fRecTracksPtr->First();
+  while (track1) {
+    removedTrack1 = kFALSE;
+    nHits1 = track1->GetNTrackHits();
+    // Loop over second track of the pair
+    track2 = (AliMUONTrack*) fRecTracksPtr->After(track1);
+    while (track2) {
+      nHits2 = track2->GetNTrackHits();
+      // number of hits in common between two tracks
+      hitsInCommon = track1->HitsInCommon(track2);
+      // check for identical tracks
+      if ((hitsInCommon == nHits1) || (hitsInCommon == nHits2)) {
+        // decide which track to remove
+        if (nHits2 > nHits1) {
+         // remove track1 and continue the first loop with the track next to track1
+         trackToRemove = track1;
+         track1 = (AliMUONTrack*) fRecTracksPtr->After(track1);
+          fRecTracksPtr->Remove(trackToRemove);
+         fRecTracksPtr->Compress(); // this is essential to retrieve the TClonesArray afterwards
+         fNRecTracks--;
+         removedTrack1 = kTRUE;
+         break;
+       } else {
+         // remove track2 and continue the second loop with the track next to track2
+         trackToRemove = track2;
+         track2 = (AliMUONTrack*) fRecTracksPtr->After(track2);
+         fRecTracksPtr->Remove(trackToRemove);
+         fRecTracksPtr->Compress(); // this is essential to retrieve the TClonesArray afterwards
+         fNRecTracks--;
+        }
+      } else track2 = (AliMUONTrack*) fRecTracksPtr->After(track2);
+    } // track2
+    if (removedTrack1) continue;
+    track1 = (AliMUONTrack*) fRecTracksPtr->After(track1);
+  } // track1
+  return;
 }
 
   //__________________________________________________________________________
-void AliMUONTrackReconstructor::CalcTrackParamAtVertex(AliMUONTrack *Track) const
+void AliMUONTrackReconstructor::RemoveDoubleTracks(void)
 {
-  /// Set track parameters at vertex.
-  /// TrackHit's are assumed to be only in stations(1..) 4 and 5,
-  /// and sorted according to increasing Z.
-  /// Parameters are calculated from information in HitForRec's
-  /// of first and last TrackHit's.
-  AliMUONTrackParam *trackParamVertex = Track->GetTrackParamAtVertex(); // pointer to track parameters at vertex
-  // Pointer to HitForRec attached to first TrackParamAtHit
-  AliMUONHitForRec *firstHit = ((AliMUONTrackParam*) (Track->GetTrackParamAtHit()->First()))->GetHitForRecPtr();
-  // Pointer to HitForRec attached to last TrackParamAtHit
-  AliMUONHitForRec *lastHit = ((AliMUONTrackParam*) (Track->GetTrackParamAtHit()->Last()))->GetHitForRecPtr();
-  // Z difference between first and last hits
-  Double_t deltaZ = firstHit->GetZ() - lastHit->GetZ();
-  // bending slope in stations(1..) 4 and 5
-  Double_t bendingSlope = (firstHit->GetBendingCoor() - lastHit->GetBendingCoor()) / deltaZ;
-  trackParamVertex->SetBendingSlope(bendingSlope);
-  // impact parameter
-  Double_t impactParam = firstHit->GetBendingCoor() - bendingSlope * firstHit->GetZ();
-  // signed bending momentum
-  trackParamVertex->SetInverseBendingMomentum(1.0 / GetBendingMomentumFromImpactParam(impactParam));
-  // bending slope at vertex
-  trackParamVertex->SetBendingSlope(bendingSlope + impactParam / fSimpleBPosition);
-  // non bending slope
-  trackParamVertex->SetNonBendingSlope((firstHit->GetNonBendingCoor() - lastHit->GetNonBendingCoor()) / deltaZ);
-  // vertex coordinates at (0,0,0)
-  trackParamVertex->SetZ(0.0);
-  trackParamVertex->SetBendingCoor(0.0);
-  trackParamVertex->SetNonBendingCoor(0.0);
+  /// To remove double tracks:
+  /// Tracks are considered identical if more than half of the hits of the track
+  /// which has the smaller number of hits are in common with the other track.
+  /// Among two identical tracks, one keeps the track with the larger number of hits
+  /// or, if these numbers are equal, the track with the minimum chi2.
+  AliMUONTrack *track1, *track2, *trackToRemove;
+  Int_t hitsInCommon, nHits1, nHits2;
+  Bool_t removedTrack1;
+  // Loop over first track of the pair
+  track1 = (AliMUONTrack*) fRecTracksPtr->First();
+  while (track1) {
+    removedTrack1 = kFALSE;
+    nHits1 = track1->GetNTrackHits();
+    // Loop over second track of the pair
+    track2 = (AliMUONTrack*) fRecTracksPtr->After(track1);
+    while (track2) {
+      nHits2 = track2->GetNTrackHits();
+      // number of hits in common between two tracks
+      hitsInCommon = track1->HitsInCommon(track2);
+      // check for identical tracks
+      if (((nHits1 < nHits2) && (2 * hitsInCommon > nHits1)) || (2 * hitsInCommon > nHits2)) {
+        // decide which track to remove
+        if ((nHits1 > nHits2) || ((nHits1 == nHits2) && (track1->GetFitFMin() <= track2->GetFitFMin()))) {
+         // remove track2 and continue the second loop with the track next to track2
+         trackToRemove = track2;
+         track2 = (AliMUONTrack*) fRecTracksPtr->After(track2);
+         fRecTracksPtr->Remove(trackToRemove);
+         fRecTracksPtr->Compress(); // this is essential to retrieve the TClonesArray afterwards
+         fNRecTracks--;
+        } else {
+         // else remove track1 and continue the first loop with the track next to track1
+         trackToRemove = track1;
+         track1 = (AliMUONTrack*) fRecTracksPtr->After(track1);
+          fRecTracksPtr->Remove(trackToRemove);
+         fRecTracksPtr->Compress(); // this is essential to retrieve the TClonesArray afterwards
+         fNRecTracks--;
+         removedTrack1 = kTRUE;
+         break;
+        }
+      } else track2 = (AliMUONTrack*) fRecTracksPtr->After(track2);
+    } // track2
+    if (removedTrack1) continue;
+    track1 = (AliMUONTrack*) fRecTracksPtr->After(track1);
+  } // track1
+  return;
 }
 
   //__________________________________________________________________________
 void AliMUONTrackReconstructor::FollowTracks(void)
 {
   /// Follow tracks in stations(1..) 3, 2 and 1
-  // too long: should be made more modular !!!!
-  AliMUONHitForRec *bestHit, *extrapHit, *hit;
-  AliMUONSegment *bestSegment, *extrapSegment, *segment;
+  AliDebug(1,"Enter FollowTracks");
+  
   AliMUONTrack *track, *nextTrack;
-  AliMUONTrackParam *trackParam1, trackParam[2], trackParamVertex;
-  // -1 to avoid compilation warnings
-  Int_t ch = -1, chInStation, chBestHit = -1, iHit, iSegment, station, trackIndex; 
-  Double_t bestChi2, chi2, dZ1, dZ2, dZ3, maxSigma2Distance, mcsFactor;
-  Double_t bendingMomentum, chi2Norm = 0.;
-
-
-  // local maxSigma2Distance, for easy increase in testing
-  maxSigma2Distance = fMaxSigma2Distance;
-  AliDebug(2,"Enter FollowTracks");
-  // Loop over track candidates
+  AliMUONTrackParam *trackParamAtFirstHit;
+  Double_t numberOfDegFree, chi2Norm;
+  Int_t CurrentNRecTracks;
+  
+  for (Int_t station = 2; station >= 0; station--) {
+    // Save the actual number of reconstructed track in case of
+    // tracks are added or suppressed during the tracking procedure
+    // !! Do not compress fRecTracksPtr until the end of the loop over tracks !!
+    CurrentNRecTracks = fNRecTracks;
+    for (Int_t iRecTrack = 0; iRecTrack <CurrentNRecTracks; iRecTrack++) {
+      AliDebug(1,Form("FollowTracks: track candidate(1..) %d", iRecTrack+1));
+      track = (AliMUONTrack*) fRecTracksPtr->UncheckedAt(iRecTrack);
+      // Fit the track:
+      // Do not take into account the multiple scattering to speed up the fit
+      // Calculate the track parameter covariance matrix
+      // If "station" is station(1..) 3 then use the vertex to better constrain the fit
+      if (station==2) {
+        SetVertexForFit(track);
+        track->SetFitWithVertex(kTRUE);
+      } else track->SetFitWithVertex(kFALSE);
+      Fit(track,kFALSE, kTRUE);
+      // Remove the track if the normalized chi2 is too high
+      numberOfDegFree = (2. * track->GetNTrackHits() - 5.);
+      if (numberOfDegFree > 0) chi2Norm = track->GetFitFMin() / numberOfDegFree;
+      else chi2Norm = 1.e10;
+      if (chi2Norm > fgkMaxNormChi2) {
+       fRecTracksPtr->Remove(track);
+       fNRecTracks--;
+       continue;
+      }
+      // Add MCS effects in parameter covariances
+      trackParamAtFirstHit = (AliMUONTrackParam*) (track->GetTrackParamAtHit()->First());
+      AliMUONTrackExtrap::AddMCSEffectInTrackParamCov(trackParamAtFirstHit,AliMUONConstants::ChamberThicknessInX0(),1.);
+      // Printout for debuging
+      if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 2) || (AliLog::GetGlobalDebugLevel() >= 2)) {
+        cout<<endl<<"Track parameter covariances at first hit with multiple Coulomb scattering effects:"<<endl;
+        trackParamAtFirstHit->GetCovariances()->Print();
+      }
+      // Look for compatible hitForRec in station(0..) "station"
+      FollowTrackInStation(track,station);
+    }
+    // Compress fRecTracksPtr for the next step
+    fRecTracksPtr->Compress();
+    // Keep only the best tracks if required
+    if (!fgkTrackAllTracks) RemoveDoubleTracks();
+  }
+  
+  // Last fit of track candidates with all station
+  // Take into account the multiple scattering and remove bad tracks
+  Int_t trackIndex = -1;
   track = (AliMUONTrack*) fRecTracksPtr->First();
-  trackIndex = -1;
   while (track) {
     trackIndex++;
     nextTrack = (AliMUONTrack*) fRecTracksPtr->After(track); // prepare next track
-    AliDebug(2,Form("FollowTracks: track candidate(0..): %d", trackIndex));
-    // Fit track candidate from parameters at vertex
-    // -> with 3 parameters (X_vertex and Y_vertex are fixed)
-    // without multiple Coulomb scattering
-    Fit(track,0,0);
-    if (AliLog::GetGlobalDebugLevel()> 2) {
-      cout << "FollowTracks: track candidate(0..): " << trackIndex
-          << " after fit in stations(0..) 3 and 4" << endl;
+    track->SetFitWithVertex(kFALSE); // just to be sure
+    Fit(track,kTRUE, kFALSE);
+    // Printout for debuging
+    if (AliLog::GetGlobalDebugLevel() >= 3) {
+      cout << "FollowTracks: track candidate(0..) " << trackIndex << " after final fit" << endl;
       track->RecursiveDump();
+    } 
+    // Remove the track if the normalized chi2 is too high
+    numberOfDegFree = (2.0 * track->GetNTrackHits() - 5);
+    if (numberOfDegFree > 0) chi2Norm = track->GetFitFMin() / numberOfDegFree;
+    else chi2Norm = 1.e10;
+    if (chi2Norm > fgkMaxNormChi2) {
+      fRecTracksPtr->Remove(track);
+      fNRecTracks--;
     }
-    // Loop over stations(1..) 3, 2 and 1
-    // something SPECIAL for stations 2 and 1 for majority 3 coincidence ????
-    // otherwise: majority coincidence 2 !!!!
-    for (station = 2; station >= 0; station--) {
-      // Track parameters at first track hit (smallest Z)
-      trackParam1 = (AliMUONTrackParam*) (track->GetTrackParamAtHit()->First());
-      // extrapolation to station
-      AliMUONTrackExtrap::ExtrapToStation(trackParam1, station, trackParam);
-      extrapSegment = new AliMUONSegment(); //  empty segment
-      // multiple scattering factor corresponding to one chamber
-      // and momentum in bending plane (not total)
-      mcsFactor = 0.0136 * trackParam1->GetInverseBendingMomentum();
-      mcsFactor        = fChamberThicknessInX0 * mcsFactor * mcsFactor;
-      // Z difference from previous station
-      dZ1 = AliMUONConstants::DefaultChamberZ(2 * station) -
-           AliMUONConstants::DefaultChamberZ(2 * station + 2);
-      // Z difference between the two previous stations
-      dZ2 = AliMUONConstants::DefaultChamberZ(2 * station + 2) -
-           AliMUONConstants::DefaultChamberZ(2 * station + 4);
-      // Z difference between the two chambers in the previous station
-      dZ3 = AliMUONConstants::DefaultChamberZ(2 * station) -
-           AliMUONConstants::DefaultChamberZ(2 * station + 1);
-      extrapSegment->SetBendingCoorReso2(fBendingResolution * fBendingResolution);
-      extrapSegment->SetNonBendingCoorReso2(fNonBendingResolution * fNonBendingResolution);
-      extrapSegment->UpdateFromStationTrackParam(trackParam, mcsFactor, dZ1, dZ2, dZ3, station,
-                                                trackParam1->GetInverseBendingMomentum());
-      bestChi2 = 5.0;
-      bestSegment = NULL;
-      if (AliLog::GetGlobalDebugLevel() > 2) {
-       cout << "FollowTracks: track candidate(0..): " << trackIndex
-            << " Look for segment in station(0..): " << station << endl;
-      }
+    track = nextTrack;
+  }
+  fRecTracksPtr->Compress();
+  
+}
 
-      // Loop over segments in station
-      for (iSegment = 0; iSegment < fNSegments[station]; iSegment++) {
-       // Look for best compatible Segment in station
-       // should consider all possibilities ????
-       // multiple scattering ????
-       // separation in 2 functions: Segment and HitForRec ????
-       segment = (AliMUONSegment*) ((*fSegmentsPtr[station])[iSegment]);
-       // correction of corrected segment (fBendingCoor and fNonBendingCoor)
-       // according to real Z value of "segment" and slopes of "extrapSegment"
-       AliMUONTrackExtrap::ExtrapToZ(&(trackParam[0]), segment->GetZ());
-       AliMUONTrackExtrap::ExtrapToZ(&(trackParam[1]), segment->GetZ()); // now same as trackParam[0] !?!?!?!?!?!
-       extrapSegment->SetBendingCoor((&(trackParam[0]))->GetBendingCoor());
-       extrapSegment->SetNonBendingCoor((&(trackParam[0]))->GetNonBendingCoor());
-       extrapSegment->SetBendingSlope((&(trackParam[0]))->GetBendingSlope());
-       extrapSegment->SetNonBendingSlope((&(trackParam[0]))->GetNonBendingSlope());
-       chi2 = segment->NormalizedChi2WithSegment(extrapSegment, maxSigma2Distance);
-       if (chi2 < bestChi2) {
-         // update best Chi2 and Segment if better found
-         bestSegment = segment;
-         bestChi2 = chi2;
-       }
+  //__________________________________________________________________________
+void AliMUONTrackReconstructor::FollowTrackInStation(AliMUONTrack* trackCandidate, Int_t nextStation)
+{
+  /// Follow trackCandidate in station(0..) nextStation and search for compatible HitForRec(s)
+  /// Keep all possibilities or only the best one(s) according to the flag fgkTrackAllTracks:
+  /// kTRUE:  duplicate "trackCandidate" if there are several possibilities and add the new tracks at the end of
+  ///         fRecTracksPtr to avoid conficts with other track candidates at this current stage of the tracking procedure.
+  ///         Remove the obsolete "trackCandidate" at the end.
+  /// kFALSE: add only the best hit(s) to the "trackCandidate". Try to add a couple of hits in priority.
+  AliDebug(1,Form("Enter FollowTrackInStation(1..) %d", nextStation+1));
+  
+  Int_t ch1 = 2*nextStation;
+  Int_t ch2 = 2*nextStation+1;
+  Double_t zCh2 = AliMUONConstants::DefaultChamberZ(ch2);
+  Double_t chi2WithOneHitForRec = 1.e10;
+  Double_t chi2WithTwoHitForRec = 1.e10;
+  Double_t maxChi2WithOneHitForRec = 2.*fgkMaxNormChi2; // 2 because 2 quantities in chi2
+  Double_t maxChi2WithTwoHitForRec = 4.*fgkMaxNormChi2; // 4 because 4 quantities in chi2
+  Double_t bestChi2WithOneHitForRec = maxChi2WithOneHitForRec;
+  Double_t bestChi2WithTwoHitForRec = maxChi2WithTwoHitForRec;
+  AliMUONTrack *newTrack = 0x0;
+  AliMUONHitForRec *hitForRecCh1, *hitForRecCh2;
+  AliMUONHitForRec *bestHitForRec1 = 0x0, *bestHitForRec2 = 0x0;
+  //
+  //Extrapolate trackCandidate to chamber "ch2" to save computing time in the next steps
+  AliMUONTrackParam *extrapTrackParamPtr = trackCandidate->GetExtrapTrackParam();
+  *extrapTrackParamPtr = *((AliMUONTrackParam*)(trackCandidate->GetTrackParamAtHit()->First()));
+  AliMUONTrackExtrap::ExtrapToZCov(extrapTrackParamPtr, zCh2);
+  AliMUONTrackParam extrapTrackParamSave(*extrapTrackParamPtr);
+  //
+  // Printout for debuging
+  if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 2) || (AliLog::GetGlobalDebugLevel() >= 2)) {
+    TMatrixD* paramCovForDebug = extrapTrackParamPtr->GetCovariances();
+    cout<<endl<<"Track parameter covariances at first hit extrapolated to z = "<<zCh2<<":"<<endl;
+    paramCovForDebug->Print();
+  }
+  //
+  // look for candidates in chamber 2 
+  // Printout for debuging
+  if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
+    cout << "FollowTrackInStation: look for hits in chamber(1..): " << ch2+1 << endl;
+  }
+  for (Int_t hit2 = 0; hit2 < fNHitsForRecPerChamber[ch2]; hit2++) {
+    hitForRecCh2 = (AliMUONHitForRec*) fHitsForRecPtr->UncheckedAt(fIndexOfFirstHitForRecPerChamber[ch2]+hit2);
+    // extrapolate track parameters and covariances only once for this hit
+    AliMUONTrackExtrap::ExtrapToZCov(extrapTrackParamPtr, hitForRecCh2->GetZ());
+    chi2WithOneHitForRec = trackCandidate->TryOneHitForRec(hitForRecCh2);
+    // if good chi2 then try to attach a hitForRec in the other chamber too
+    if (chi2WithOneHitForRec < maxChi2WithOneHitForRec) {
+      // Printout for debuging
+      if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
+        cout << "FollowTrackInStation: look for second hits in chamber(1..): " << ch1+1 << endl;
       }
-      if (bestSegment) {
-       // best segment found: add it to track candidate
-       AliMUONTrackExtrap::ExtrapToZ(&(trackParam[0]), bestSegment->GetZ());
-       track->AddTrackParamAtHit(&(trackParam[0]),bestSegment->GetHitForRec1());
-       AliMUONTrackExtrap::ExtrapToZ(&(trackParam[1]), bestSegment->GetZ()); // now same as trackParam[0] !?!?!?!?!?!
-       track->AddTrackParamAtHit(&(trackParam[1]),bestSegment->GetHitForRec2());
-       AliDebug(3, Form("FollowTracks: track candidate(0..): %d  Added segment in station(0..): %d", trackIndex, station));
-       if (AliLog::GetGlobalDebugLevel()>2) track->RecursiveDump();
-      } else {
-       // No best segment found:
-       // Look for best compatible HitForRec in station:
-       // should consider all possibilities ????
-       // multiple scattering ???? do about like for extrapSegment !!!!
-       extrapHit = new AliMUONHitForRec(); //  empty hit
-       bestChi2 = 3.0;
-       bestHit = NULL;
-       AliDebug(3, Form("FollowTracks: track candidate(0..): %d Segment not found, look for hit in station(0..): %d ", 
-                        trackIndex, station));
-       
-       // Loop over chambers of the station
-       for (chInStation = 0; chInStation < 2; chInStation++) {
-         ch = 2 * station + chInStation;
-         for (iHit = fIndexOfFirstHitForRecPerChamber[ch]; iHit < fIndexOfFirstHitForRecPerChamber[ch]+fNHitsForRecPerChamber[ch]; iHit++) {
-           hit = (AliMUONHitForRec*) ((*fHitsForRecPtr)[iHit]);
-           // coordinates of extrapolated hit
-           AliMUONTrackExtrap::ExtrapToZ(&(trackParam[chInStation]), hit->GetZ());
-           extrapHit->SetBendingCoor((&(trackParam[chInStation]))->GetBendingCoor());
-           extrapHit->SetNonBendingCoor((&(trackParam[chInStation]))->GetNonBendingCoor());
-           // resolutions from "extrapSegment"
-           extrapHit->SetBendingReso2(extrapSegment->GetBendingCoorReso2());
-           extrapHit->SetNonBendingReso2(extrapSegment->GetNonBendingCoorReso2());
-           // Loop over hits in the chamber
-           // condition for hit not already in segment ????
-           chi2 = hit->NormalizedChi2WithHitForRec(extrapHit, maxSigma2Distance);
-           if (chi2 < bestChi2) {
-             // update best Chi2 and HitForRec if better found
-             bestHit = hit;
-             bestChi2 = chi2;
-             chBestHit = chInStation;
+      Bool_t foundSecondHit = kFALSE;
+      for (Int_t hit1 = 0; hit1 < fNHitsForRecPerChamber[ch1]; hit1++) {
+        hitForRecCh1 = (AliMUONHitForRec*) fHitsForRecPtr->UncheckedAt(fIndexOfFirstHitForRecPerChamber[ch1]+hit1);
+       chi2WithTwoHitForRec = trackCandidate->TryTwoHitForRec(hitForRecCh2, hitForRecCh1); // order hits like that to save computing time
+        // if good chi2 then create a new track by adding the 2 hitForRec to the "trackCandidate"
+       if (chi2WithTwoHitForRec < maxChi2WithTwoHitForRec) {
+         foundSecondHit = kTRUE;
+          if (fgkTrackAllTracks) {
+           // copy trackCandidate into a new track put at the end of fRecTracksPtr and add the new hitForRec's
+            newTrack = new ((*fRecTracksPtr)[fRecTracksPtr->GetLast()+1]) AliMUONTrack(*trackCandidate);
+           fNRecTracks++;
+            AliMUONTrackParam TrackParam1(extrapTrackParamSave);
+            AliMUONTrackExtrap::ExtrapToZ(&TrackParam1, hitForRecCh1->GetZ());
+           newTrack->AddTrackParamAtHit(&TrackParam1,hitForRecCh1);
+            AliMUONTrackParam TrackParam2(extrapTrackParamSave);
+            AliMUONTrackExtrap::ExtrapToZ(&TrackParam2, hitForRecCh2->GetZ());
+           newTrack->AddTrackParamAtHit(&TrackParam2,hitForRecCh2);
+            // Sort TrackParamAtHit according to increasing Z
+            newTrack->GetTrackParamAtHit()->Sort();
+           // Update the chi2 of the new track
+           if (newTrack->GetFitFMin()<0) newTrack->SetFitFMin(chi2WithTwoHitForRec);
+           else newTrack->SetFitFMin(newTrack->GetFitFMin() + chi2WithTwoHitForRec);
+           // Printout for debuging
+           if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
+             cout << "FollowTrackInStation: added two hits in station(1..): " << nextStation+1
+                  << " (Chi2 = " << chi2WithTwoHitForRec << ")" << endl;
+             if (AliLog::GetGlobalDebugLevel() >= 3) newTrack->RecursiveDump();
            }
-         }
-       }
-       if (bestHit) {
-         // best hit found: add it to track candidate
-         AliMUONTrackExtrap::ExtrapToZ(&(trackParam[chBestHit]), bestHit->GetZ());
-         track->AddTrackParamAtHit(&(trackParam[chBestHit]),bestHit);
-         if (AliLog::GetGlobalDebugLevel() > 2) {
-           cout << "FollowTracks: track candidate(0..): " << trackIndex
-                << " Added hit in station(0..): " << station << endl;
-           track->RecursiveDump();
-         }
-       } else {
-         // Remove current track candidate
-         // and corresponding TrackHit's, ...
-         fRecTracksPtr->Remove(track);
-         fNRecTracks--;
-         delete extrapSegment;
-         delete extrapHit;
-         break; // stop the search for this candidate:
-         // exit from the loop over station
+          } else if (chi2WithTwoHitForRec < bestChi2WithTwoHitForRec) {
+           // keep track of the best couple of hits
+           bestChi2WithTwoHitForRec = chi2WithTwoHitForRec;
+           bestHitForRec1 = hitForRecCh1;
+           bestHitForRec2 = hitForRecCh2;
+          }
        }
-       delete extrapHit;
-      }
-      delete extrapSegment;
-      // Sort TrackParamAtHit according to increasing Z
-      track->GetTrackParamAtHit()->Sort();
-      // Update track parameters at first track hit (smallest Z)
-      trackParam1 = (AliMUONTrackParam*) (track->GetTrackParamAtHit()->First());
-      bendingMomentum = 0.;
-      if (TMath::Abs(trackParam1->GetInverseBendingMomentum()) > 0.)
-       bendingMomentum = TMath::Abs(1/(trackParam1->GetInverseBendingMomentum()));
-      // Track removed if bendingMomentum not in window [min, max]
-      if ((bendingMomentum < fMinBendingMomentum) || (bendingMomentum > fMaxBendingMomentum)) {
-       fRecTracksPtr->Remove(track);
-       fNRecTracks--;
-       break; // stop the search for this candidate:
-       // exit from the loop over station 
       }
-      // Track fit from parameters at first hit
-      // -> with 5 parameters (momentum and position)
-      // with multiple Coulomb scattering if all stations
-      if (station == 0) Fit(track,1,1);
-      // without multiple Coulomb scattering if not all stations
-      else Fit(track,1,0);
-      Double_t numberOfDegFree = (2.0 * track->GetNTrackHits() - 5);
-      if (numberOfDegFree > 0) {
-        chi2Norm =  track->GetFitFMin() / numberOfDegFree;
-      } else {
-       chi2Norm = 1.e10;
+      // if no hitForRecCh1 found then consider to add hitForRecCh2 only
+      if (!foundSecondHit) {
+        if (fgkTrackAllTracks) {
+         // copy trackCandidate into a new track put at the end of fRecTracksPtr and add the new hitForRec's
+          newTrack = new ((*fRecTracksPtr)[fRecTracksPtr->GetLast()+1]) AliMUONTrack(*trackCandidate);
+         fNRecTracks++;
+          AliMUONTrackParam TrackParam1(extrapTrackParamSave);
+          AliMUONTrackExtrap::ExtrapToZ(&TrackParam1, hitForRecCh2->GetZ());
+          newTrack->AddTrackParamAtHit(&TrackParam1,hitForRecCh2);
+          // Sort TrackParamAtHit according to increasing Z
+          newTrack->GetTrackParamAtHit()->Sort();
+         // Update the chi2 of the new track
+         if (newTrack->GetFitFMin()<0) newTrack->SetFitFMin(chi2WithOneHitForRec);
+         else newTrack->SetFitFMin(newTrack->GetFitFMin() + chi2WithOneHitForRec);
+         // Printout for debuging
+         if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
+           cout << "FollowTrackInStation: added one hit in chamber(1..): " << ch2+1
+                << " (Chi2 = " << chi2WithOneHitForRec << ")" << endl;
+           if (AliLog::GetGlobalDebugLevel() >= 3) newTrack->RecursiveDump();
+         }
+       } if (!bestHitForRec2 && chi2WithOneHitForRec < bestChi2WithOneHitForRec) {
+         // keep track of the best single hitForRec except if a couple
+          // of hits has already been found (i.e. bestHitForRec2!=0x0)
+         bestChi2WithOneHitForRec = chi2WithOneHitForRec;
+         bestHitForRec1 = hitForRecCh2;
+        }
       }
-      // Track removed if normalized chi2 too high
-      if (chi2Norm > fMaxChi2) {
-       fRecTracksPtr->Remove(track);
-       fNRecTracks--;
-       break; // stop the search for this candidate:
-       // exit from the loop over station 
+    }
+    // reset the extrapolated track parameter for next step
+    trackCandidate->SetExtrapTrackParam(&extrapTrackParamSave);
+  }
+  //
+  // look for candidates in chamber 1 not already attached to a track
+  // if we want to keep all possible tracks or if no good couple of hitForRec has been found
+  if (fgkTrackAllTracks || !bestHitForRec2) {
+    if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
+      cout << "FollowTrackInStation: look for single hits in chamber(1..): " << ch1+1 << endl;
+    }
+    for (Int_t hit1 = 0; hit1 < fNHitsForRecPerChamber[ch1]; hit1++) {
+      hitForRecCh1 = (AliMUONHitForRec*) fHitsForRecPtr->UncheckedAt(fIndexOfFirstHitForRecPerChamber[ch1]+hit1);
+      if (hitForRecCh1->GetNTrackHits() >= 1) continue; // Skip hitForRec already used
+      chi2WithOneHitForRec = trackCandidate->TryOneHitForRec(hitForRecCh1);
+      // if good chi2 then create a new track by adding the good hitForRec in "ch1" to the "trackCandidate"
+      // We do not try to attach a hitForRec in the other chamber too since it has already been done above
+      if (chi2WithOneHitForRec < maxChi2WithOneHitForRec) {
+       if (fgkTrackAllTracks) {
+         // copy trackCandidate into a new track put at the end of fRecTracksPtr and add the new hitForRec's
+         newTrack = new ((*fRecTracksPtr)[fRecTracksPtr->GetLast()+1]) AliMUONTrack(*trackCandidate);
+         fNRecTracks++;
+         AliMUONTrackParam TrackParam1(extrapTrackParamSave);
+         AliMUONTrackExtrap::ExtrapToZ(&TrackParam1, hitForRecCh1->GetZ());
+         newTrack->AddTrackParamAtHit(&TrackParam1,hitForRecCh1);
+         // Sort TrackParamAtHit according to increasing Z
+         newTrack->GetTrackParamAtHit()->Sort();
+         // Update the chi2 of the new track
+         if (newTrack->GetFitFMin()<0) newTrack->SetFitFMin(chi2WithOneHitForRec);
+         else newTrack->SetFitFMin(newTrack->GetFitFMin() + chi2WithOneHitForRec);
+         // Printout for debuging
+         if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
+           cout << "FollowTrackInStation: added one hit in chamber(1..): " << ch1+1
+                << " (Chi2 = " << chi2WithOneHitForRec << ")" << endl;
+           if (AliLog::GetGlobalDebugLevel() >= 3) newTrack->RecursiveDump();
+         }
+       } if (!bestHitForRec2 && chi2WithOneHitForRec < bestChi2WithOneHitForRec) {
+         // keep track of the best single hitForRec except if a couple
+         // of hits has already been found (i.e. bestHitForRec1!=0x0)
+         bestChi2WithOneHitForRec = chi2WithOneHitForRec;
+         bestHitForRec1 = hitForRecCh1;
+       }
       }
-      if (AliLog::GetGlobalDebugLevel() > 2) {
-       cout << "FollowTracks: track candidate(0..): " << trackIndex
-            << " after fit from station(0..): " << station << " to 4" << endl;
-       track->RecursiveDump();
+    }
+  }
+  //
+  // fill out the best track if required else clean up the fRecTracksPtr array
+  if (!fgkTrackAllTracks && bestHitForRec1) {
+    AliMUONTrackParam TrackParam1(extrapTrackParamSave);
+    AliMUONTrackExtrap::ExtrapToZ(&TrackParam1, bestHitForRec1->GetZ());
+    trackCandidate->AddTrackParamAtHit(&TrackParam1,bestHitForRec1);
+    if (bestHitForRec2) {
+      AliMUONTrackParam TrackParam2(extrapTrackParamSave);
+      AliMUONTrackExtrap::ExtrapToZ(&TrackParam2, bestHitForRec2->GetZ());
+      trackCandidate->AddTrackParamAtHit(&TrackParam2,bestHitForRec2);
+      // Sort TrackParamAtHit according to increasing Z
+      trackCandidate->GetTrackParamAtHit()->Sort();
+      // Update the chi2 of the new track
+      if (trackCandidate->GetFitFMin()<0) trackCandidate->SetFitFMin(bestChi2WithTwoHitForRec);
+      else trackCandidate->SetFitFMin(trackCandidate->GetFitFMin() + bestChi2WithTwoHitForRec);
+      // Printout for debuging
+      if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
+        cout << "FollowTrackInStation: added the two best hits in station(1..): " << nextStation+1
+             << " (Chi2 = " << bestChi2WithTwoHitForRec << ")" << endl;
+        if (AliLog::GetGlobalDebugLevel() >= 3) newTrack->RecursiveDump();
       }
-      // Track extrapolation to the vertex through the absorber (Branson)
-      // after going through the first station
-      if (station == 0) {
-       trackParamVertex = *((AliMUONTrackParam*) (track->GetTrackParamAtHit()->First()));
-       AliMUONTrackExtrap::ExtrapToVertex(&trackParamVertex, 0., 0., 0.);
-       track->SetTrackParamAtVertex(&trackParamVertex);
-       if (AliLog::GetGlobalDebugLevel() > 0) {
-         cout << "FollowTracks: track candidate(0..): " << trackIndex
-              << " after extrapolation to vertex" << endl;
-         track->RecursiveDump();
-       }
+    } else {
+      // Sort TrackParamAtHit according to increasing Z
+      trackCandidate->GetTrackParamAtHit()->Sort();
+      // Update the chi2 of the new track
+      if (trackCandidate->GetFitFMin()<0) trackCandidate->SetFitFMin(bestChi2WithOneHitForRec);
+      else trackCandidate->SetFitFMin(trackCandidate->GetFitFMin() + bestChi2WithOneHitForRec);
+      // Printout for debuging
+      if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
+        cout << "FollowTrackInStation: added the best hit in station(1..): " << nextStation+1
+             << " (Chi2 = " << bestChi2WithOneHitForRec << ")" << endl;
+        if (AliLog::GetGlobalDebugLevel() >= 3) newTrack->RecursiveDump();
       }
-    } // for (station = 2;...
-    // go really to next track
-    track = nextTrack;
-  } // while (track)
-  // Compression of track array (necessary after Remove)
-  fRecTracksPtr->Compress();
-  return;
+    }
+  } else {
+    fRecTracksPtr->Remove(trackCandidate); // obsolete track
+    fNRecTracks--;
+  }
+  
 }
 
   //__________________________________________________________________________
-void AliMUONTrackReconstructor::Fit(AliMUONTrack *Track, Int_t FitStart, Int_t FitMCS)
+void AliMUONTrackReconstructor::SetVertexForFit(AliMUONTrack* trackCandidate)
 {
-  /// Fit the track "Track",
-  /// with or without multiple Coulomb scattering according to "FitMCS",
-  /// starting, according to "FitStart",
-  /// with track parameters at vertex or at the first TrackHit.
+  /// Add the vertex as a measured hit to constrain the fit of the "trackCandidate"
+  /// Compute the vertex resolution from natural vertex dispersion and
+  /// multiple scattering effets according to trackCandidate path in absorber
+  /// It is necessary to account for multiple scattering effects here instead of during the fit of
+  /// the "trackCandidate" to do not influence the result by changing track resolution at vertex
+  AliDebug(1,"Enter SetVertexForFit");
   
-  if ((FitStart != 0) && (FitStart != 1)) {
-    cout << "ERROR in AliMUONTrackReconstructor::Fit(...)" << endl;
-    cout << "FitStart = " << FitStart << " is neither 0 nor 1" << endl;
-    exit(0);
-  }
-  if ((FitMCS != 0) && (FitMCS != 1)) {
-    cout << "ERROR in AliMUONTrackReconstructor::Fit(...)" << endl;
-    cout << "FitMCS = " << FitMCS << " is neither 0 nor 1" << endl;
-    exit(0);
-  }
+  Double_t nonBendingReso2 = fNonBendingVertexDispersion * fNonBendingVertexDispersion;
+  Double_t bendingReso2 = fBendingVertexDispersion * fBendingVertexDispersion;
+  // add multiple scattering effets
+  AliMUONTrackParam paramAtVertex(*((AliMUONTrackParam*)(trackCandidate->GetTrackParamAtHit()->First())));
+  paramAtVertex.DeleteCovariances(); // to be sure to account only for multiple scattering
+  AliMUONTrackExtrap::ExtrapToVertexUncorrected(&paramAtVertex,0.);
+  TMatrixD* paramCov = paramAtVertex.GetCovariances();
+  nonBendingReso2 += (*paramCov)(0,0);
+  bendingReso2 += (*paramCov)(2,2);
+  // Set the vertex
+  AliMUONHitForRec vertex; // Coordinates set to (0.,0.,0.) by default
+  vertex.SetNonBendingReso2(nonBendingReso2);
+  vertex.SetBendingReso2(bendingReso2);
+  trackCandidate->SetVertex(&vertex);
+}
+
+  //__________________________________________________________________________
+void AliMUONTrackReconstructor::Fit(AliMUONTrack *track, Bool_t includeMCS, Bool_t calcCov)
+{
+  /// Fit the track "track" with or without multiple Coulomb scattering according to "includeMCS".
   
-  Double_t arg[1], benC, errorParam, invBenP, lower, nonBenC, upper, x, y;
-  char parName[50];
+  Double_t benC, errorParam, invBenP, nonBenC, x, y;
   AliMUONTrackParam *trackParam;
-  // Check if Minuit is initialized...
-  fgFitter = TVirtualFitter::Fitter(Track,5);
-  fgFitter->Clear(); // necessary ???? probably yes
-  // how to reset the printout number at every fit ????
-  // is there any risk to leave it like that ????
-  // how to go faster ???? choice of Minuit parameters like EDM ????
-  // choice of function to be minimized according to fFitMCS
-  if (FitMCS == 0) fgFitter->SetFCN(TrackChi2);
-  else fgFitter->SetFCN(TrackChi2MCS);
-  // Switch off printout
-  arg[0] = -1;
-  fgFitter->ExecuteCommand("SET PRINT", arg, 1); // More printing !!!!
+  Double_t arg[1], fedm, errdef, fitFMin;
+  Int_t npari, nparx;
+  Int_t status, covStatus;
+  
+  // Clear MINUIT parameters
+  gMinuit->mncler();
+  // Give the fitted track to MINUIT
+  gMinuit->SetObjectFit(track);
+  if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 2) || (AliLog::GetGlobalDebugLevel() >= 2)) {
+    // Define print level
+    arg[0] = 1;
+    gMinuit->mnexcm("SET PRI", arg, 1, status);
+    // Print covariance matrix
+    gMinuit->mnexcm("SHO COV", arg, 0, status);
+  } else {
+    arg[0] = -1;
+    gMinuit->mnexcm("SET PRI", arg, 1, status);
+  }
   // No warnings
-  fgFitter->ExecuteCommand("SET NOW", arg, 0);
-  // Parameters according to "fFitStart"
-  // (should be a function to be used at every place where needed ????)
-  if (FitStart == 0) trackParam = Track->GetTrackParamAtVertex();
-  else trackParam = (AliMUONTrackParam*) (Track->GetTrackParamAtHit()->First());
-  // set first 3 Minuit parameters
+  gMinuit->mnexcm("SET NOW", arg, 0, status);
+  // Define strategy
+  //arg[0] = 2;
+  //gMinuit->mnexcm("SET STR", arg, 1, status);
+  
+  // Switch between available FCN according to "includeMCS"
+  if (includeMCS) gMinuit->SetFCN(TrackChi2MCS);
+  else gMinuit->SetFCN(TrackChi2);
+  
+  // Set fitted parameters (!! The order is very important for the covariance matrix !!)
+  trackParam = (AliMUONTrackParam*) (track->GetTrackParamAtHit()->First());
   // could be tried with no limits for the search (min=max=0) ????
-  fgFitter->SetParameter(0, "InvBenP",
-                        trackParam->GetInverseBendingMomentum(),
-                        0.003, -0.4, 0.4);
-  fgFitter->SetParameter(1, "BenS",
-                        trackParam->GetBendingSlope(),
-                        0.001, -0.5, 0.5);
-  fgFitter->SetParameter(2, "NonBenS",
-                        trackParam->GetNonBendingSlope(),
-                        0.001, -0.5, 0.5);
-  if (FitStart == 1) {
-    // set last 2 Minuit parameters when we start from first track hit
-    // mandatory limits in Bending to avoid NaN values of parameters
-    fgFitter->SetParameter(3, "X",
-                          trackParam->GetNonBendingCoor(),
-                          0.03, -500.0, 500.0);
-    // mandatory limits in non Bending to avoid NaN values of parameters
-    fgFitter->SetParameter(4, "Y",
-                          trackParam->GetBendingCoor(),
-                          0.10, -500.0, 500.0);
-  }
-  // search without gradient calculation in the function
-  fgFitter->ExecuteCommand("SET NOGRADIENT", arg, 0);
+  // mandatory limits in non Bending to avoid NaN values of parameters
+  gMinuit->mnparm(0, "X", trackParam->GetNonBendingCoor(), 0.03, -500.0, 500.0, status);
+  gMinuit->mnparm(1, "NonBenS", trackParam->GetNonBendingSlope(), 0.001, -0.5, 0.5, status);
+  // mandatory limits in Bending to avoid NaN values of parameters
+  gMinuit->mnparm(2, "Y", trackParam->GetBendingCoor(), 0.10, -500.0, 500.0, status);
+  gMinuit->mnparm(3, "BenS", trackParam->GetBendingSlope(), 0.001, -0.5, 0.5, status);
+  gMinuit->mnparm(4, "InvBenP", trackParam->GetInverseBendingMomentum(), 0.003, -0.4, 0.4, status);
+  
   // minimization
-  fgFitter->ExecuteCommand("MINIMIZE", arg, 0);
-  // exit from Minuit
-  //  fgFitter->ExecuteCommand("EXIT", arg, 0); // necessary ????
+  gMinuit->mnexcm("MIGRAD", arg, 0, status);
+  
+  // Calculate the covariance matrix more accurately if required
+  if (calcCov) gMinuit->mnexcm("HESSE", arg, 0, status);
+  
   // get results into "invBenP", "benC", "nonBenC" ("x", "y")
-  fgFitter->GetParameter(0, parName, invBenP, errorParam, lower, upper);
-  trackParam->SetInverseBendingMomentum(invBenP);
-  fgFitter->GetParameter(1, parName, benC, errorParam, lower, upper);
-  trackParam->SetBendingSlope(benC);
-  fgFitter->GetParameter(2, parName, nonBenC, errorParam, lower, upper);
+  gMinuit->GetParameter(0, x, errorParam);
+  trackParam->SetNonBendingCoor(x);
+  gMinuit->GetParameter(1, nonBenC, errorParam);
   trackParam->SetNonBendingSlope(nonBenC);
-  if (FitStart == 1) {
-    fgFitter->GetParameter(3, parName, x, errorParam, lower, upper);
-    trackParam->SetNonBendingCoor(x);
-    fgFitter->GetParameter(4, parName, y, errorParam, lower, upper);
-    trackParam->SetBendingCoor(y);
-  }
+  gMinuit->GetParameter(2, y, errorParam);
+  trackParam->SetBendingCoor(y);
+  gMinuit->GetParameter(3, benC, errorParam);
+  trackParam->SetBendingSlope(benC);
+  gMinuit->GetParameter(4, invBenP, errorParam);
+  trackParam->SetInverseBendingMomentum(invBenP);
+  
   // global result of the fit
-  Double_t fedm, errdef, fitFMin;
-  Int_t npari, nparx;
-  fgFitter->GetStats(fitFMin, fedm, errdef, npari, nparx);
-  Track->SetFitFMin(fitFMin);
+  gMinuit->mnstat(fitFMin, fedm, errdef, npari, nparx, covStatus);
+  track->SetFitFMin(fitFMin);
+  
+  // Get the covariance matrix if required
+  if (calcCov) {
+    // Covariance matrix according to HESSE status
+    // If problem then keep only the diagonal terms (variances)
+    Double_t matrix[5][5];
+    gMinuit->mnemat(&matrix[0][0],5);
+    if (covStatus == 3) trackParam->SetCovariances(matrix);
+    else trackParam->SetVariances(matrix);
+  } else *(trackParam->GetCovariances()) = 0;
+  
 }
 
   //__________________________________________________________________________
-void TrackChi2(Int_t &NParam, Double_t * /*Gradient*/, Double_t &Chi2, Double_t *Param, Int_t /*Flag*/)
+void TrackChi2(Int_t & /*NParam*/, Double_t * /*Gradient*/, Double_t &Chi2, Double_t *Param, Int_t /*Flag*/)
 {
   /// Return the "Chi2" to be minimized with Minuit for track fitting,
   /// with "NParam" parameters
@@ -686,24 +709,37 @@ void TrackChi2(Int_t &NParam, Double_t * /*Gradient*/, Double_t &Chi2, Double_t
   /// Assumes that the track hits are sorted according to increasing Z.
   /// Track parameters at each TrackHit are updated accordingly.
   /// Multiple Coulomb scattering is not taken into account
-  AliMUONTrack *trackBeingFitted;
+  
+  AliMUONTrack *trackBeingFitted = (AliMUONTrack*) gMinuit->GetObjectFit();
+//  AliMUONTrack *trackBeingFitted = (AliMUONTrack*) AliMUONTrackReconstructor::Fitter()->GetObjectFit();
   AliMUONTrackParam param1;
   AliMUONTrackParam* trackParamAtHit;
   AliMUONHitForRec* hitForRec;
   Chi2 = 0.0; // initialize Chi2
+  Double_t dX, dY;
+  
   // copy of track parameters to be fitted
-  trackBeingFitted = (AliMUONTrack*) AliMUONTrackReconstructor::Fitter()->GetObjectFit();
-  // 3 parameters means fit track candidate from parameters at vertex (X_vertex and Y_vertex are fixed)
-  if (NParam == 3) param1 = *(trackBeingFitted->GetTrackParamAtVertex());
-  else param1 = *((AliMUONTrackParam*) (trackBeingFitted->GetTrackParamAtHit()->First()));
-  // Minuit parameters to be fitted into this copy
-  param1.SetInverseBendingMomentum(Param[0]);
-  param1.SetBendingSlope(Param[1]);
-  param1.SetNonBendingSlope(Param[2]);
-  if (NParam == 5) {
-    param1.SetNonBendingCoor(Param[3]);
-    param1.SetBendingCoor(Param[4]);
+  param1 = *((AliMUONTrackParam*) (trackBeingFitted->GetTrackParamAtHit()->First()));
+  param1.SetNonBendingCoor(Param[0]);
+  param1.SetNonBendingSlope(Param[1]);
+  param1.SetBendingCoor(Param[2]);
+  param1.SetBendingSlope(Param[3]);
+  param1.SetInverseBendingMomentum(Param[4]);
+  
+  // Take the vertex into account in the fit if required
+  if (trackBeingFitted->GetFitWithVertex()) {
+    AliMUONTrackParam paramAtVertex(param1);
+    AliMUONTrackExtrap::ExtrapToZ(&paramAtVertex, 0.);
+    AliMUONHitForRec *vertex = trackBeingFitted->GetVertex();
+    if (!vertex) {
+      cout<<"Error in TrackChi2: Want to use the vertex in tracking but it has not been created!!"<<endl;
+      exit(-1);
+    }
+    dX = vertex->GetNonBendingCoor() - paramAtVertex.GetNonBendingCoor();
+    dY = vertex->GetBendingCoor() - paramAtVertex.GetBendingCoor();
+    Chi2 += dX * dX / vertex->GetNonBendingReso2() + dY * dY / vertex->GetBendingReso2();
   }
+  
   // Follow track through all planes of track hits
   trackParamAtHit = (AliMUONTrackParam*) (trackBeingFitted->GetTrackParamAtHit()->First());
   while (trackParamAtHit) {
@@ -715,16 +751,15 @@ void TrackChi2(Int_t &NParam, Double_t * /*Gradient*/, Double_t &Chi2, Double_t
     // Increment Chi2
     // done hit per hit, with hit resolution,
     // and not with point and angle like in "reco_muon.F" !!!!
-    // Needs to add multiple scattering contribution ????
-    Double_t dX = hitForRec->GetNonBendingCoor() - param1.GetNonBendingCoor();
-    Double_t dY = hitForRec->GetBendingCoor() - param1.GetBendingCoor();
+    dX = hitForRec->GetNonBendingCoor() - param1.GetNonBendingCoor();
+    dY = hitForRec->GetBendingCoor() - param1.GetBendingCoor();
     Chi2 = Chi2 + dX * dX / hitForRec->GetNonBendingReso2() + dY * dY / hitForRec->GetBendingReso2();
     trackParamAtHit = (AliMUONTrackParam*) (trackBeingFitted->GetTrackParamAtHit()->After(trackParamAtHit));
   }
 }
 
   //__________________________________________________________________________
-void TrackChi2MCS(Int_t &NParam, Double_t * /*Gradient*/, Double_t &Chi2, Double_t *Param, Int_t /*Flag*/)
+void TrackChi2MCS(Int_t & /*NParam*/, Double_t * /*Gradient*/, Double_t &Chi2, Double_t *Param, Int_t /*Flag*/)
 {
   /// Return the "Chi2" to be minimized with Minuit for track fitting,
   /// with "NParam" parameters
@@ -732,25 +767,13 @@ void TrackChi2MCS(Int_t &NParam, Double_t * /*Gradient*/, Double_t &Chi2, Double
   /// Assumes that the track hits are sorted according to increasing Z.
   /// Track parameters at each TrackHit are updated accordingly.
   /// Multiple Coulomb scattering is taken into account with covariance matrix.
-  AliMUONTrack *trackBeingFitted;
+  
+  AliMUONTrack *trackBeingFitted = (AliMUONTrack*) gMinuit->GetObjectFit();
+//  AliMUONTrack *trackBeingFitted = (AliMUONTrack*) AliMUONTrackReconstructor::Fitter()->GetObjectFit();
   AliMUONTrackParam param1;
   AliMUONTrackParam* trackParamAtHit;
   AliMUONHitForRec* hitForRec;
   Chi2 = 0.0; // initialize Chi2
-  // copy of track parameters to be fitted
-  trackBeingFitted = (AliMUONTrack*) AliMUONTrackReconstructor::Fitter()->GetObjectFit();
-  // 3 parameters means fit track candidate from parameters at vertex (X_vertex and Y_vertex are fixed)
-  if (NParam == 3) param1 = *(trackBeingFitted->GetTrackParamAtVertex());
-  else param1 = *((AliMUONTrackParam*) (trackBeingFitted->GetTrackParamAtHit()->First()));
-  // Minuit parameters to be fitted into this copy
-  param1.SetInverseBendingMomentum(Param[0]);
-  param1.SetBendingSlope(Param[1]);
-  param1.SetNonBendingSlope(Param[2]);
-  if (NParam == 5) {
-    param1.SetNonBendingCoor(Param[3]);
-    param1.SetBendingCoor(Param[4]);
-  }
-
   Int_t chCurrent, chPrev = 0, hitNumber, hitNumber1, hitNumber2, hitNumber3;
   Double_t z1, z2, z3;
   AliMUONTrackParam *trackParamAtHit1, *trackParamAtHit2, *trackParamAtHit3;
@@ -761,8 +784,30 @@ void TrackChi2MCS(Int_t &NParam, Double_t * /*Gradient*/, Double_t &Chi2, Double
   TMatrixD *covBending = new TMatrixD(numberOfHit, numberOfHit);
   TMatrixD *covNonBending = new TMatrixD(numberOfHit, numberOfHit);
   Double_t *msa2 = new Double_t[numberOfHit];
+  
+  // copy of track parameters to be fitted
+  param1 = *((AliMUONTrackParam*) (trackBeingFitted->GetTrackParamAtHit()->First()));
+  param1.SetNonBendingCoor(Param[0]);
+  param1.SetNonBendingSlope(Param[1]);
+  param1.SetBendingCoor(Param[2]);
+  param1.SetBendingSlope(Param[3]);
+  param1.SetInverseBendingMomentum(Param[4]);
 
-  // Predicted coordinates and  multiple scattering angles are first calculated
+  // Take the vertex into account in the fit if required
+  if (trackBeingFitted->GetFitWithVertex()) {
+    AliMUONTrackParam paramAtVertex(param1);
+    AliMUONTrackExtrap::ExtrapToZ(&paramAtVertex, 0.);
+    AliMUONHitForRec *vertex = trackBeingFitted->GetVertex();
+    if (!vertex) {
+      cout<<"Error in TrackChi2MCS: Want to use the vertex in tracking but it has not been created!!"<<endl;
+      exit(-1);
+    }
+    Double_t dX = vertex->GetNonBendingCoor() - paramAtVertex.GetNonBendingCoor();
+    Double_t dY = vertex->GetBendingCoor() - paramAtVertex.GetBendingCoor();
+    Chi2 += dX * dX / vertex->GetNonBendingReso2() + dY * dY / vertex->GetBendingReso2();
+  }
+  
+  // Predicted coordinates and multiple scattering angles are first calculated
   for (hitNumber = 0; hitNumber < numberOfHit; hitNumber++) {
     trackParamAtHit = (AliMUONTrackParam*) (trackBeingFitted->GetTrackParamAtHit()->UncheckedAt(hitNumber));
     hitForRec = trackParamAtHit->GetHitForRecPtr();
@@ -816,15 +861,10 @@ void TrackChi2MCS(Int_t &NParam, Double_t * /*Gradient*/, Double_t &Chi2, Double
   } // for (hitNumber1 = 0;...
     
   // Inversion of covariance matrices
-  // with "mnvertLocal", local "mnvert" function of Minuit.
-  // One cannot use directly "mnvert" since "TVirtualFitter" does not know it.
-  // One will have to replace this local function by the right inversion function
-  // from a specialized Root package for symmetric positive definite matrices,
-  // when available!!!!
   Int_t ifailBending;
-  mnvertLocal(&((*covBending)(0,0)), numberOfHit, numberOfHit, numberOfHit, ifailBending);
+  gMinuit->mnvert(&((*covBending)(0,0)), numberOfHit, numberOfHit, numberOfHit, ifailBending);
   Int_t ifailNonBending;
-  mnvertLocal(&((*covNonBending)(0,0)), numberOfHit, numberOfHit, numberOfHit, ifailNonBending);
+  gMinuit->mnvert(&((*covNonBending)(0,0)), numberOfHit, numberOfHit, numberOfHit, ifailNonBending);
 
   // It would be worth trying to calculate the inverse of the covariance matrix
   // only once per fit, since it cannot change much in principle,
@@ -870,18 +910,18 @@ void TrackChi2MCS(Int_t &NParam, Double_t * /*Gradient*/, Double_t &Chi2, Double
   delete [] msa2;
 }
 
+  //__________________________________________________________________________
 Double_t MultipleScatteringAngle2(AliMUONTrackParam *param)
 {
   /// Returns square of multiple Coulomb scattering angle
   /// from TrackParamAtHit pointed to by "param"
   Double_t slopeBending, slopeNonBending, radiationLength, inverseBendingMomentum2, inverseTotalMomentum2;
   Double_t varMultipleScatteringAngle;
-  // Better implementation in AliMUONTrack class ????
   slopeBending = param->GetBendingSlope();
   slopeNonBending = param->GetNonBendingSlope();
   // thickness in radiation length for the current track,
   // taking local angle into account
-  radiationLength = AliMUONConstants::DefaultChamberThicknessInX0() *
+  radiationLength = AliMUONConstants::ChamberThicknessInX0() *
                    TMath::Sqrt(1.0 + slopeBending*slopeBending + slopeNonBending*slopeNonBending);
   inverseBendingMomentum2 =  param->GetInverseBendingMomentum() * param->GetInverseBendingMomentum();
   inverseTotalMomentum2 = inverseBendingMomentum2 * (1.0 + slopeBending * slopeBending) /
@@ -892,159 +932,44 @@ Double_t MultipleScatteringAngle2(AliMUONTrackParam *param)
   return varMultipleScatteringAngle;
 }
 
-//______________________________________________________________________________
- void mnvertLocal(Double_t *a, Int_t l, Int_t, Int_t n, Int_t &ifail)
+  //__________________________________________________________________________
+void AliMUONTrackReconstructor::ExtrapTracksToVertex()
 {
-///*-*-*-*-*-*-*-*-*-*-*-*Inverts a symmetric matrix*-*-*-*-*-*-*-*-*-*-*-*-*
-///*-*                    ==========================
-///*-*        inverts a symmetric matrix.   matrix is first scaled to
-///*-*        have all ones on the diagonal (equivalent to change of units)
-///*-*        but no pivoting is done since matrix is positive-definite.
-///*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
-
-  // taken from TMinuit package of Root (l>=n)
-  // fVERTs, fVERTq and fVERTpp changed to localVERTs, localVERTq and localVERTpp
-  //  Double_t localVERTs[n], localVERTq[n], localVERTpp[n];
-  Double_t * localVERTs = new Double_t[n];
-  Double_t * localVERTq = new Double_t[n];
-  Double_t * localVERTpp = new Double_t[n];
-  // fMaxint changed to localMaxint
-  Int_t localMaxint = n;
-
-    /* System generated locals */
-    Int_t aOffset;
-
-    /* Local variables */
-    Double_t si;
-    Int_t i, j, k, kp1, km1;
-
-    /* Parameter adjustments */
-    aOffset = l + 1;
-    a -= aOffset;
-
-    /* Function Body */
-    ifail = 0;
-    if (n < 1) goto L100;
-    if (n > localMaxint) goto L100;
-//*-*-                  scale matrix by sqrt of diag elements
-    for (i = 1; i <= n; ++i) {
-        si = a[i + i*l];
-        if (si <= 0) goto L100;
-        localVERTs[i-1] = 1 / TMath::Sqrt(si);
-    }
-    for (i = 1; i <= n; ++i) {
-        for (j = 1; j <= n; ++j) {
-            a[i + j*l] = a[i + j*l]*localVERTs[i-1]*localVERTs[j-1];
-        }
-    }
-//*-*-                                       . . . start main loop . . . .
-    for (i = 1; i <= n; ++i) {
-        k = i;
-//*-*-                  preparation for elimination step1
-        if (a[k + k*l] != 0) localVERTq[k-1] = 1 / a[k + k*l];
-        else goto L100;
-        localVERTpp[k-1] = 1;
-        a[k + k*l] = 0;
-        kp1 = k + 1;
-        km1 = k - 1;
-        if (km1 < 0) goto L100;
-        else if (km1 == 0) goto L50;
-        else               goto L40;
-L40:
-        for (j = 1; j <= km1; ++j) {
-            localVERTpp[j-1] = a[j + k*l];
-            localVERTq[j-1]  = a[j + k*l]*localVERTq[k-1];
-            a[j + k*l]   = 0;
-        }
-L50:
-        if (k - n < 0) goto L51;
-        else if (k - n == 0) goto L60;
-        else                goto L100;
-L51:
-        for (j = kp1; j <= n; ++j) {
-            localVERTpp[j-1] = a[k + j*l];
-            localVERTq[j-1]  = -a[k + j*l]*localVERTq[k-1];
-            a[k + j*l]   = 0;
-        }
-//*-*-                  elimination proper
-L60:
-        for (j = 1; j <= n; ++j) {
-            for (k = j; k <= n; ++k) { a[j + k*l] += localVERTpp[j-1]*localVERTq[k-1]; }
-        }
+  /// propagate tracks to the vertex through the absorber (Branson)
+  AliMUONTrack *track;
+  AliMUONTrackParam trackParamVertex;
+  AliMUONHitForRec *vertex;
+  Double_t vX, vY, vZ;
+  
+  for (Int_t iRecTrack = 0; iRecTrack <fNRecTracks; iRecTrack++) {
+    track = (AliMUONTrack*) fRecTracksPtr->UncheckedAt(iRecTrack);
+    trackParamVertex = *((AliMUONTrackParam*)(track->GetTrackParamAtHit()->First()));
+    vertex = track->GetVertex();
+    if (vertex) { // if track as a vertex defined, use it
+      vX = vertex->GetNonBendingCoor();
+      vY = vertex->GetBendingCoor();
+      vZ = vertex->GetZ();
+    } else { // else vertex = (0.,0.,0.)
+      vX = 0.;
+      vY = 0.;
+      vZ = 0.;
     }
-//*-*-                  elements of left diagonal and unscaling
-    for (j = 1; j <= n; ++j) {
-        for (k = 1; k <= j; ++k) {
-            a[k + j*l] = a[k + j*l]*localVERTs[k-1]*localVERTs[j-1];
-            a[j + k*l] = a[k + j*l];
-        }
+    AliMUONTrackExtrap::ExtrapToVertex(&trackParamVertex, vX, vY, vZ);
+    track->SetTrackParamAtVertex(&trackParamVertex);
+    
+    if (AliLog::GetGlobalDebugLevel() > 0) {
+      cout << "FollowTracks: track candidate(0..): " << iRecTrack
+          << " after extrapolation to vertex" << endl;
+      track->RecursiveDump();
     }
-    delete [] localVERTs;
-    delete [] localVERTq;
-    delete [] localVERTpp;
-    return;
-//*-*-                  failure return
-L100:
-    delete [] localVERTs;
-    delete [] localVERTq;
-    delete [] localVERTpp;
-    ifail = 1;
-} /* mnvertLocal */
-
-  //__________________________________________________________________________
-void AliMUONTrackReconstructor::RemoveDoubleTracks(void)
-{
-  /// To remove double tracks.
-  /// Tracks are considered identical
-  /// if they have at least half of their hits in common.
-  /// Among two identical tracks, one keeps the track with the larger number of hits
-  /// or, if these numbers are equal, the track with the minimum Chi2.
-  AliMUONTrack *track1, *track2, *trackToRemove;
-  Int_t hitsInCommon, nHits1, nHits2;
-  Bool_t removedTrack1;
-  // Loop over first track of the pair
-  track1 = (AliMUONTrack*) fRecTracksPtr->First();
-  while (track1) {
-    removedTrack1 = kFALSE;
-    nHits1 = track1->GetNTrackHits();
-    // Loop over second track of the pair
-    track2 = (AliMUONTrack*) fRecTracksPtr->After(track1);
-    while (track2) {
-      nHits2 = track2->GetNTrackHits();
-      // number of hits in common between two tracks
-      hitsInCommon = track1->HitsInCommon(track2);
-      // check for identical tracks
-      if ((4 * hitsInCommon) >= (nHits1 + nHits2)) {
-        // decide which track to remove
-        if ((nHits1 > nHits2) || ((nHits1 == nHits2) && (track1->GetFitFMin() < track2->GetFitFMin()))) {
-         // remove track2 and continue the second loop with the track next to track2
-         trackToRemove = track2;
-         track2 = (AliMUONTrack*) fRecTracksPtr->After(track2);
-         fRecTracksPtr->Remove(trackToRemove);
-         fNRecTracks--;
-         fRecTracksPtr->Compress(); // this is essential to retrieve the TClonesArray afterwards
-        } else {
-         // else remove track1 and continue the first loop with the track next to track1
-         trackToRemove = track1;
-         track1 = (AliMUONTrack*) fRecTracksPtr->After(track1);
-          fRecTracksPtr->Remove(trackToRemove);
-         fNRecTracks--;
-         fRecTracksPtr->Compress(); // this is essential to retrieve the TClonesArray afterwards
-         removedTrack1 = kTRUE;
-         break;
-        }
-      } else track2 = (AliMUONTrack*) fRecTracksPtr->After(track2);
-    } // track2
-    if (removedTrack1) continue;
-    track1 = (AliMUONTrack*) fRecTracksPtr->After(track1);
-  } // track1
-  return;
+  }
+  
 }
 
   //__________________________________________________________________________
-void AliMUONTrackReconstructor::UpdateHitForRecAtHit()
+void AliMUONTrackReconstructor::FillMUONTrack()
 {
-  /// Set cluster parameters after track fitting. Fill fHitForRecAtHit of AliMUONTrack's
+  /// Fill fHitForRecAtHit of AliMUONTrack's
   AliMUONTrack *track;
   AliMUONTrackParam *trackParamAtHit;
   track = (AliMUONTrack*) fRecTracksPtr->First();
index a4742cc..5a36263 100644 (file)
@@ -12,9 +12,7 @@
 #include <TObject.h>
 #include "AliMUONVTrackReconstructor.h"
 
-class AliMUONSegment;
 class AliMUONTrack;
-class TVirtualFitter;
 
 class AliMUONTrackReconstructor : public AliMUONVTrackReconstructor {
 
@@ -22,9 +20,6 @@ class AliMUONTrackReconstructor : public AliMUONVTrackReconstructor {
   AliMUONTrackReconstructor(AliMUONData* data); // default Constructor
   virtual ~AliMUONTrackReconstructor(); // Destructor
 
-           /// Return track fitter
-  static TVirtualFitter* Fitter(void) {return fgFitter;}
-  
   virtual void EventDump(void);  // dump reconstructed event
 
 
@@ -32,32 +27,28 @@ class AliMUONTrackReconstructor : public AliMUONVTrackReconstructor {
 
   // Functions
   virtual void AddHitsForRecFromRawClusters();
-  virtual void MakeSegments(void);
   virtual void MakeTracks(void);
   virtual void MakeTrackCandidates(void);
   virtual void FollowTracks(void);
   virtual void RemoveDoubleTracks(void);
+  virtual void ExtrapTracksToVertex(void);
+  virtual void FillMUONTrack(void);
   
 
  private:
   
-  // Defaults parameters for reconstruction
-  static const Double_t fgkDefaultMaxChi2; ///< default max. track chi2 for reconstruction
-
-  static TVirtualFitter* fgFitter; //!< Pointer to track fitter
+  // Parameters for reconstruction
+  static const Double_t fgkMaxNormChi2; ///< maximum Chi2 per degree of freedom for reconstruction
+  static const Bool_t fgkTrackAllTracks; /// kTRUE to track all the possible candidates; kFALSE to track only the best ones
 
-  // Parameters for track reconstruction
-  Double_t fMaxChi2; ///< maximum Chi2 per degree of Freedom
-  
   // Functions
   AliMUONTrackReconstructor (const AliMUONTrackReconstructor& rhs); ///< copy constructor
   AliMUONTrackReconstructor& operator=(const AliMUONTrackReconstructor& rhs); ///< assignment operator
   
-  Int_t MakeTrackCandidatesWithTwoSegments(AliMUONSegment *BegSegment);
-  Int_t MakeTrackCandidatesWithOneSegmentAndOnePoint(AliMUONSegment *BegSegment);
-  void CalcTrackParamAtVertex(AliMUONTrack *Track) const;
-  void Fit(AliMUONTrack *Track, Int_t FitStart, Int_t FitMCS);
-  void UpdateHitForRecAtHit(void);
+  void RemoveIdenticalTracks(void);
+  void FollowTrackInStation(AliMUONTrack* trackCandidate, Int_t nextStation);
+  void SetVertexForFit(AliMUONTrack* trackCandidate);
+  void Fit(AliMUONTrack *track, Bool_t includeMCS, Bool_t calcCov);
 
 
   ClassDef(AliMUONTrackReconstructor, 0) // MUON track reconstructor in ALICE
index 858c3b0..1905064 100644 (file)
@@ -38,8 +38,8 @@
 #include "AliMUONData.h"
 #include "AliMUONConstants.h"
 #include "AliMUONHitForRec.h"
+#include "AliMUONObjectPair.h"
 #include "AliMUONRawCluster.h"
-#include "AliMUONSegment.h"
 #include "AliMUONTrackK.h" 
 #include "AliLog.h"
 
@@ -79,7 +79,7 @@ void AliMUONTrackReconstructorK::AddHitsForRecFromRawClusters()
   TTree *treeR;
   AliMUONHitForRec *hitForRec;
   AliMUONRawCluster *clus;
-  Int_t iclus, nclus, nTRentries;
+  Int_t iclus, nclus;
   TClonesArray *rawclusters;
   AliDebug(1,"Enter AddHitsForRecFromRawClusters");
 
@@ -90,24 +90,12 @@ void AliMUONTrackReconstructorK::AddHitsForRecFromRawClusters()
   }
   
   if (fTrackMethod != 3) { //AZ
-    nTRentries = Int_t(treeR->GetEntries());
-    
-    if (!(fMUONData->IsRawClusterBranchesInTree())) {
-      AliError(Form("RawCluster information is not avalaible, nTRentries = %d not equal to 1",nTRentries));
-      exit(0);
-    }
-
     fMUONData->SetTreeAddress("RC"); //AZ
     fMUONData->GetRawClusters(); // only one entry  
   }
 
   // Loop over tracking chambers
   for (Int_t ch = 0; ch < AliMUONConstants::NTrackingCh(); ch++) {
-    // number of HitsForRec to 0 for the chamber
-    fNHitsForRecPerChamber[ch] = 0;
-    // index of first HitForRec for the chamber
-    if (ch == 0) fIndexOfFirstHitForRecPerChamber[ch] = 0;
-    else fIndexOfFirstHitForRecPerChamber[ch] = fNHitsForRec;
     rawclusters =fMUONData->RawClusters(ch);
     nclus = (Int_t) (rawclusters->GetEntries());
     // Loop over (cathode correlated) raw clusters
@@ -117,7 +105,6 @@ void AliMUONTrackReconstructorK::AddHitsForRecFromRawClusters()
       // and increment number of AliMUONHitForRec's (total and in chamber)
       hitForRec = new ((*fHitsForRecPtr)[fNHitsForRec]) AliMUONHitForRec(clus);
       fNHitsForRec++;
-      (fNHitsForRecPerChamber[ch])++;
       // more information into HitForRec
       hitForRec->SetBendingReso2(clus->GetErrY() * clus->GetErrY());
       hitForRec->SetNonBendingReso2(clus->GetErrX() * clus->GetErrX());
@@ -126,7 +113,6 @@ void AliMUONTrackReconstructorK::AddHitsForRecFromRawClusters()
       hitForRec->SetHitNumber(iclus);
       // Z coordinate of the raw cluster (cm)
       hitForRec->SetZ(clus->GetZ(0));
-      
       StdoutToAliDebug(3,
                        cout << "Chamber " << ch <<
                        " raw cluster  " << iclus << " : " << endl;
@@ -158,32 +144,6 @@ void AliMUONTrackReconstructorK::AddHitsForRecFromRawClusters()
 }
 
   //__________________________________________________________________________
-void AliMUONTrackReconstructorK::MakeSegments(void)
-{
-  /// To make the list of segments in all stations,
-  /// from the list of hits to be reconstructed
-  AliDebug(1,"Enter MakeSegments");
-  //AZ ResetSegments();
-  // Loop over stations
-  for (Int_t st = 3; st < AliMUONConstants::NTrackingCh()/2; st++) MakeSegmentsPerStation(st); 
-  
-  StdoutToAliDebug(3,
-    cout << "end of MakeSegments" << endl;
-    for (Int_t st = 0; st < AliMUONConstants::NTrackingCh()/2; st++) 
-    {
-      cout << "station " << st
-           << "  has " << fNSegments[st] << " segments:"
-           << endl;
-      for (Int_t seg = 0; seg < fNSegments[st]; seg++) 
-      {
-             ((*fSegmentsPtr[st])[seg])->Print();
-      }
-    }
-                   );
-  return;
-}
-
-  //__________________________________________________________________________
 void AliMUONTrackReconstructorK::MakeTracks(void)
 {
   /// To make the tracks from the list of segments and points in all stations
@@ -196,8 +156,8 @@ void AliMUONTrackReconstructorK::MakeTracks(void)
   FollowTracks();
   // Remove double tracks
   RemoveDoubleTracks();
-  // Propagate tracks to the vertex thru absorber
-  GoToVertex();
+  // Propagate tracks to the vertex through absorber
+  ExtrapTracksToVertex();
   // Fill AliMUONTrack data members
   FillMUONTrack();
 }
@@ -205,9 +165,10 @@ void AliMUONTrackReconstructorK::MakeTracks(void)
   //__________________________________________________________________________
 void AliMUONTrackReconstructorK::MakeTrackCandidates(void)
 {
-  /// To make initial tracks for Kalman filter from the list of segments
+  /// To make initial tracks for Kalman filter from segments in stations(1..)  4 and 5
   Int_t istat, iseg;
-  AliMUONSegment *segment;
+  TClonesArray *segments;
+  AliMUONObjectPair *segment;
   AliMUONTrackK *trackK;
 
   AliDebug(1,"Enter MakeTrackCandidatesK");
@@ -215,10 +176,12 @@ void AliMUONTrackReconstructorK::MakeTrackCandidates(void)
   AliMUONTrackK a(this, fHitsForRecPtr);
   // Loop over stations(1...) 5 and 4
   for (istat=4; istat>=3; istat--) {
+    // Make segments in the station
+    segments = MakeSegmentsInStation(istat);
     // Loop over segments in the station
-    for (iseg=0; iseg<fNSegments[istat]; iseg++) {
-      // Transform segments to tracks and evaluate covariance matrix
-      segment = (AliMUONSegment*) ((*fSegmentsPtr[istat])[iseg]);
+    for (iseg=0; iseg<segments->GetEntriesFast(); iseg++) {
+      // Transform segments to tracks
+      segment = (AliMUONObjectPair*) ((*segments)[iseg]);
       trackK = new ((*fRecTracksPtr)[fNRecTracks++]) AliMUONTrackK(segment);
     } // for (iseg=0;...)
   } // for (istat=4;...)
@@ -231,15 +194,16 @@ void AliMUONTrackReconstructorK::FollowTracks(void)
   /// Follow tracks using Kalman filter
   Bool_t ok;
   Int_t icand, ichamBeg = 0, ichamEnd, chamBits;
-  Double_t zDipole1, zDipole2;
   AliMUONTrackK *trackK;
   AliMUONHitForRec *hit;
   AliMUONRawCluster *clus;
   TClonesArray *rawclusters;
   clus = 0; rawclusters = 0;
 
-  zDipole1 = fSimpleBPosition + fSimpleBLength/2;
-  zDipole2 = zDipole1 - fSimpleBLength;
+  Double_t simpleBPosition = 0.5 * (AliMUONConstants::CoilZ() + AliMUONConstants::YokeZ());
+  Double_t simpleBLength = 0.5 * (AliMUONConstants::CoilL() + AliMUONConstants::YokeL());
+  Double_t zDipole1 = simpleBPosition + 0.5 * simpleBLength;
+  Double_t zDipole2 = zDipole1 - simpleBLength;
 
   // Print hits
   trackK = (AliMUONTrackK*) ((*fRecTracksPtr)[0]);
@@ -366,6 +330,7 @@ Bool_t AliMUONTrackReconstructorK::CheckCandidate(Int_t icand, Int_t nSeeds) con
   /// the same seed segment hits as hits of a good track found before)
   AliMUONTrackK *track1, *track2;
   AliMUONHitForRec *hit1, *hit2, *hit;
+  AliMUONObjectPair *segment1, *segment2;
 
   track1 = (AliMUONTrackK*) ((*fRecTracksPtr)[icand]);
   hit1 = (AliMUONHitForRec*) (*track1->GetTrackHits())[0]; // 1'st hit
@@ -375,8 +340,10 @@ Bool_t AliMUONTrackReconstructorK::CheckCandidate(Int_t icand, Int_t nSeeds) con
     track2 = (AliMUONTrackK*) ((*fRecTracksPtr)[i]);
     //if (track2->GetRecover() < 0) continue;
     if (track2->GetRecover() < 0 && icand >= nSeeds) continue;
-
-    if (track1->GetStartSegment() == track2->GetStartSegment()) {
+    segment1 = track1->GetStartSegment();
+    segment2 = track2->GetStartSegment();
+    if (segment1->First()  == segment2->First() &&
+        segment1->Second() == segment2->Second()) {
       return kFALSE;
     } else {
       Int_t nSame = 0;
@@ -427,7 +394,7 @@ void AliMUONTrackReconstructorK::RemoveDoubleTracks(void)
 }
 
   //__________________________________________________________________________
-void AliMUONTrackReconstructorK::GoToVertex(void)
+void AliMUONTrackReconstructorK::ExtrapTracksToVertex(void)
 {
   /// Propagates track to the vertex thru absorber
   /// (using Branson correction for now)
index c904e79..ae29f58 100644 (file)
@@ -32,11 +32,12 @@ class AliMUONTrackReconstructorK : public AliMUONVTrackReconstructor {
   
   // Functions
   virtual void AddHitsForRecFromRawClusters();
-  virtual void MakeSegments(void);
   virtual void MakeTracks(void);
   virtual void MakeTrackCandidates(void);
   virtual void FollowTracks(void);
   virtual void RemoveDoubleTracks(void);
+  virtual void ExtrapTracksToVertex(void);
+  virtual void FillMUONTrack(void);
   
 
  private:
@@ -50,8 +51,6 @@ class AliMUONTrackReconstructorK : public AliMUONVTrackReconstructor {
   AliMUONTrackReconstructorK& operator=(const AliMUONTrackReconstructorK& rhs); ///< assignment operator
   
   Bool_t CheckCandidate(Int_t icand, Int_t nSeeds) const;
-  void GoToVertex(void);
-  void FillMUONTrack(void); // set track parameters at hits for Kalman track
 
 
   ClassDef(AliMUONTrackReconstructorK, 0) // MUON track reconstructor in ALICE
index cf4adca..05e50bb 100644 (file)
 #include <stdlib.h>
 #include <Riostream.h>
 #include <TTree.h>
+#include <TClonesArray.h>
 
 #include "AliMUONVTrackReconstructor.h"
 #include "AliMUONData.h"
 #include "AliMUONConstants.h"
 #include "AliMUONHitForRec.h"
+#include "AliMUONObjectPair.h"
 #include "AliMUONTriggerTrack.h"
 #include "AliMUONTriggerCircuit.h"
 #include "AliMUONLocalTrigger.h"
 #include "AliMUONGlobalTrigger.h"
-#include "AliMUONSegment.h"
 #include "AliMUONTrack.h"
 #include "AliMUONTrackParam.h"
 #include "AliMUONTrackExtrap.h"
@@ -57,16 +58,9 @@ const Double_t AliMUONVTrackReconstructor::fgkDefaultMinBendingMomentum = 3.0;
 const Double_t AliMUONVTrackReconstructor::fgkDefaultMaxBendingMomentum = 3000.0;
 const Double_t AliMUONVTrackReconstructor::fgkDefaultBendingResolution = 0.01;
 const Double_t AliMUONVTrackReconstructor::fgkDefaultNonBendingResolution = 0.144;
-const Double_t AliMUONVTrackReconstructor::fgkDefaultMaxSigma2Distance = 16.0;
-// Simple magnetic field:
-// Value taken from macro MUONtracking.C: 0.7 T, hence 7 kG
-// Length and Position from reco_muon.F, with opposite sign:
-// Length = ZMAGEND-ZCOIL
-// Position = (ZMAGEND+ZCOIL)/2
-// to be ajusted differently from real magnetic field ????
-const Double_t AliMUONVTrackReconstructor::fgkDefaultSimpleBValue = 7.0;
-const Double_t AliMUONVTrackReconstructor::fgkDefaultSimpleBLength = 428.0;
-const Double_t AliMUONVTrackReconstructor::fgkDefaultSimpleBPosition = 1019.0;
+const Double_t AliMUONVTrackReconstructor::fgkDefaultBendingVertexDispersion = 10.;
+const Double_t AliMUONVTrackReconstructor::fgkDefaultNonBendingVertexDispersion = 10.;
+const Double_t AliMUONVTrackReconstructor::fgkDefaultMaxNormChi2MatchTrigger = 16.0;
 
 //__________________________________________________________________________
 AliMUONVTrackReconstructor::AliMUONVTrackReconstructor(AliMUONData* data)
@@ -75,19 +69,15 @@ AliMUONVTrackReconstructor::AliMUONVTrackReconstructor(AliMUONData* data)
     fMaxBendingMomentum(fgkDefaultMaxBendingMomentum),
     fBendingResolution(fgkDefaultBendingResolution),
     fNonBendingResolution(fgkDefaultNonBendingResolution),
-    fMaxSigma2Distance(fgkDefaultMaxSigma2Distance),
-    fChamberThicknessInX0(AliMUONConstants::DefaultChamberThicknessInX0()),
-    fSimpleBValue(fgkDefaultSimpleBValue),
-    fSimpleBLength(fgkDefaultSimpleBLength),
-    fSimpleBPosition(fgkDefaultSimpleBPosition),
+    fBendingVertexDispersion(fgkDefaultBendingVertexDispersion),
+    fNonBendingVertexDispersion(fgkDefaultNonBendingVertexDispersion),
+    fMaxNormChi2MatchTrigger(fgkDefaultMaxNormChi2MatchTrigger),
     fSegmentMaxDistBending(0x0),
     fSegmentMaxDistNonBending(0x0),
     fHitsForRecPtr(0x0),
     fNHitsForRec(0),
     fNHitsForRecPerChamber(0x0),
     fIndexOfFirstHitForRecPerChamber(0x0),
-    fSegmentsPtr(0x0),
-    fNSegments(0x0),
     fRecTracksPtr(0x0),
     fNRecTracks(0),
     fMUONData(data),
@@ -106,28 +96,9 @@ AliMUONVTrackReconstructor::AliMUONVTrackReconstructor(AliMUONData* data)
   // Is 10000 the right size ????
   fHitsForRecPtr = new TClonesArray("AliMUONHitForRec", 10000);
 
-  // Memory allocation for the TClonesArray's of segments in stations
-  // Is 2000 the right size ????
-  fSegmentsPtr = new TClonesArray*[AliMUONConstants::NTrackingSt()];
-  fNSegments = new Int_t[AliMUONConstants::NTrackingSt()];
-  for (Int_t st = 0; st < AliMUONConstants::NTrackingSt(); st++) {
-    fSegmentsPtr[st] = new TClonesArray("AliMUONSegment", 2000);
-    fNSegments[st] = 0; // really needed or GetEntriesFast sufficient ????
-  }
-
+  // set the magnetic field for track extrapolations
   const AliMagF* kField = AliTracker::GetFieldMap();
   if (!kField) AliFatal("No field available");
- // Sign of fSimpleBValue according to sign of Bx value at (50,50,-950).
-  Float_t b[3], x[3];
-  x[0] = 50.; x[1] = 50.; x[2] = -950.;
-  kField->Field(x,b);
-
-  fSimpleBValue    = TMath::Sign(fSimpleBValue,(Double_t) b[0]);
-  fSimpleBPosition = TMath::Sign(fSimpleBPosition,(Double_t) x[2]);
-  // See how to get fSimple(BValue, BLength, BPosition)
-  // automatically calculated from the actual magnetic field ????
-  
-  // set the magnetic field for track extrapolations
   AliMUONTrackExtrap::SetField(kField);
 }
 
@@ -141,10 +112,6 @@ AliMUONVTrackReconstructor::~AliMUONVTrackReconstructor(void)
   delete [] fIndexOfFirstHitForRecPerChamber;
   delete fTriggerTrack;
   delete fHitsForRecPtr;
-  // Correct destruction of everything ????
-  for (Int_t st = 0; st < AliMUONConstants::NTrackingSt(); st++) delete fSegmentsPtr[st];
-  delete [] fSegmentsPtr;
-  delete [] fNSegments;
 }
 
   //__________________________________________________________________________
@@ -159,7 +126,7 @@ void AliMUONVTrackReconstructor::SetReconstructionParametersToDefaults(void)
   // Maximum distance in non bending plane
   // 5 * 0.22 just to remember the way it was made in TRACKF_STAT
   // SIGCUT*DYMAX(IZ)
-  for (Int_t st = 0; st < AliMUONConstants::NTrackingCh()/2; st++)
+  for (Int_t st = 0; st < AliMUONConstants::NTrackingSt(); st++)
     fSegmentMaxDistNonBending[st] = 5. * 0.22;
   // Maximum distance in bending plane:
   // values from TRACKF_STAT, corresponding to (J psi 20cm),
@@ -184,10 +151,8 @@ void AliMUONVTrackReconstructor::EventReconstruct(void)
   // To reconstruct one event
   AliDebug(1,"Enter EventReconstruct");
   ResetTracks(); //AZ
-  ResetSegments(); //AZ
   ResetHitsForRec(); //AZ
   AddHitsForRecFromRawClusters();
-  MakeSegments();
   MakeTracks();
   if (fMUONData->IsTriggerTrackBranchesInTree()) 
     ValidateTracksWithTrigger(); 
@@ -211,17 +176,6 @@ void AliMUONVTrackReconstructor::ResetTracks(void)
 }
 
   //__________________________________________________________________________
-void AliMUONVTrackReconstructor::ResetSegments(void)
-{
-  /// To reset the TClonesArray of segments and the number of Segments for all stations
-  for (Int_t st = 0; st < AliMUONConstants::NTrackingCh()/2; st++) {
-    if (fSegmentsPtr[st]) fSegmentsPtr[st]->Clear();
-    fNSegments[st] = 0;
-  }
-  return;
-}
-
-  //__________________________________________________________________________
 void AliMUONVTrackReconstructor::ResetHitsForRec(void)
 {
   /// To reset the array and the number of HitsForRec,
@@ -265,49 +219,31 @@ void AliMUONVTrackReconstructor::SortHitsForRecWithIncreasingChamber()
 }
 
   //__________________________________________________________________________
-void AliMUONVTrackReconstructor::MakeSegmentsPerStation(Int_t Station)
+TClonesArray* AliMUONVTrackReconstructor::MakeSegmentsInStation(Int_t station)
 {
-  /// To make the list of segments in station number "Station" (0...)
-  /// from the list of hits to be reconstructed.
-  /// Updates "fNSegments"[Station].
-  /// Segments in stations 4 and 5 are sorted
-  /// according to increasing absolute value of "impact parameter"
+  /// To make the list of segments in station(0..) "Station" from the list of hits to be reconstructed.
+  /// Return a new TClonesArray of segments.
+  /// It is the responsibility of the user to delete it afterward.
+  AliDebug(1,Form("Enter MakeSegmentsPerStation (0...) %d",station));
+  
   AliMUONHitForRec *hit1Ptr, *hit2Ptr;
-  AliMUONSegment *segment;
-  Bool_t last2st;
-  Double_t bendingSlope, distBend, distNonBend, extBendCoor, extNonBendCoor,
-      impactParam = 0., maxImpactParam = 0., minImpactParam = 0.; // =0 to avoid compilation warnings.
-  AliDebug(1,Form("Enter MakeSegmentsPerStation (0...) %d",Station));
+  AliMUONObjectPair *segment;
+  Double_t bendingSlope, distBend, distNonBend, extBendCoor, extNonBendCoor, extrapFact;
+  Double_t impactParam = 0., bendingMomentum = 0.; // to avoid compilation warning
   // first and second chambers (0...) in the station
-  Int_t ch1 = 2 * Station;
+  Int_t ch1 = 2 * station;
   Int_t ch2 = ch1 + 1;
-  // variable true for stations downstream of the dipole:
-  // Station(0..4) equal to 3 or 4
-  if ((Station == 3) || (Station == 4)) {
-    last2st = kTRUE;
-    // maximum impact parameter (cm) according to fMinBendingMomentum...
-    maxImpactParam =
-      TMath::Abs(GetImpactParamFromBendingMomentum(fMinBendingMomentum));
-    // minimum impact parameter (cm) according to fMaxBendingMomentum...
-    minImpactParam =
-      TMath::Abs(GetImpactParamFromBendingMomentum(fMaxBendingMomentum));
-  }
-  else last2st = kFALSE;
-  // extrapolation factor from Z of first chamber to Z of second chamber
-  // dZ to be changed to take into account fine structure of chambers ????
-  Double_t extrapFact;
-  // index for current segment
-  Int_t segmentIndex = 0;
-  // Loop over HitsForRec in the first chamber of the station
+  // list of segments
+  TClonesArray *segments = new TClonesArray("AliMUONObjectPair", fNHitsForRecPerChamber[ch2]);
+  // Loop over HitForRec's in the first chamber of the station
   for (Int_t hit1 = fIndexOfFirstHitForRecPerChamber[ch1];
        hit1 < fIndexOfFirstHitForRecPerChamber[ch1] + fNHitsForRecPerChamber[ch1];
        hit1++) {
     // pointer to the HitForRec
     hit1Ptr = (AliMUONHitForRec*) ((*fHitsForRecPtr)[hit1]);
-    // extrapolation,
-    // on the straight line joining the HitForRec to the vertex (0,0,0),
-    // to the Z of the second chamber of the station
-    // Loop over HitsForRec in the second chamber of the station
+    // extrapolation, on the straight line joining the HitForRec to the vertex (0,0,0),
+    // to the Z of the HitForRec in the second chamber of the station
+    // Loop over HitsForRec's in the second chamber of the station
     for (Int_t hit2 = fIndexOfFirstHitForRecPerChamber[ch2];
         hit2 < fIndexOfFirstHitForRecPerChamber[ch2] + fNHitsForRecPerChamber[ch2];
         hit2++) {
@@ -321,36 +257,26 @@ void AliMUONVTrackReconstructor::MakeSegmentsPerStation(Int_t Station)
       extNonBendCoor = extrapFact * hit1Ptr->GetNonBendingCoor();
       distBend = TMath::Abs(hit2Ptr->GetBendingCoor() - extBendCoor);
       distNonBend = TMath::Abs(hit2Ptr->GetNonBendingCoor() - extNonBendCoor);
-      if (last2st) {
-       // bending slope
-       if ( hit1Ptr->GetZ() - hit2Ptr->GetZ() != 0.0 ) {
-         bendingSlope = (hit1Ptr->GetBendingCoor() - hit2Ptr->GetBendingCoor()) /
-           (hit1Ptr->GetZ() - hit2Ptr->GetZ());
-         // absolute value of impact parameter
-         impactParam =
-           TMath::Abs(hit1Ptr->GetBendingCoor() - hit1Ptr->GetZ() * bendingSlope);
-        } 
-        else {
-          AliWarning("hit1Ptr->GetZ() = hit2Ptr->GetZ(): impactParam set to maxImpactParam");
-          impactParam = maxImpactParam;   
-        }   
-      }
+      // bending slope
+      if ( hit1Ptr->GetZ() - hit2Ptr->GetZ() != 0. ) {
+        bendingSlope = (hit1Ptr->GetBendingCoor() - hit2Ptr->GetBendingCoor()) / (hit1Ptr->GetZ() - hit2Ptr->GetZ());
+        // impact parameter
+        impactParam = hit1Ptr->GetBendingCoor() - hit1Ptr->GetZ() * bendingSlope;
+        // absolute value of bending momentum
+       bendingMomentum = TMath::Abs(AliMUONTrackExtrap::GetBendingMomentumFromImpactParam(impactParam));
+      } else {
+        AliWarning("hit1Ptr->GetZ() = hit2Ptr->GetZ(): no segment created");
+        continue;
+      }   
       // check for distances not too large,
       // and impact parameter not too big if stations downstream of the dipole.
       // Conditions "distBend" and "impactParam" correlated for these stations ????
-      if ((distBend < fSegmentMaxDistBending[Station]) &&
-         (distNonBend < fSegmentMaxDistNonBending[Station]) &&
-         (!last2st || (impactParam < maxImpactParam)) &&
-         (!last2st || (impactParam > minImpactParam))) {
+      if ((distBend < fSegmentMaxDistBending[station]) && (distNonBend < fSegmentMaxDistNonBending[station]) &&
+         (bendingMomentum < fMaxBendingMomentum) && (bendingMomentum > fMinBendingMomentum)) {
        // make new segment
-       segment = new ((*fSegmentsPtr[Station])[segmentIndex])
-         AliMUONSegment(hit1Ptr, hit2Ptr);
-       // update "link" to this segment from the hit in the first chamber
-       if (hit1Ptr->GetNSegments() == 0)
-         hit1Ptr->SetIndexOfFirstSegment(segmentIndex);
-       hit1Ptr->SetNSegments(hit1Ptr->GetNSegments() + 1);
+       segment = new ((*segments)[segments->GetLast()+1]) AliMUONObjectPair(hit1Ptr, hit2Ptr, kFALSE, kFALSE);
        if (AliLog::GetGlobalDebugLevel() > 1) {
-         cout << "segmentIndex(0...): " << segmentIndex
+         cout << "segmentIndex(0...): " << segments->GetLast()
               << "  distBend: " << distBend
               << "  distNonBend: " << distNonBend
               << endl;
@@ -359,50 +285,20 @@ void AliMUONVTrackReconstructor::MakeSegmentsPerStation(Int_t Station)
          hit1Ptr->Dump();
          cout << "HitForRec in second chamber" << endl;
          hit2Ptr->Dump();
-       };
-       // increment index for current segment
-       segmentIndex++;
+       }
       }
     } //for (Int_t hit2
   } // for (Int_t hit1...
-  fNSegments[Station] = segmentIndex;
-  // Sorting according to "impact parameter" if station(1..5) 4 or 5,
-  // i.e. Station(0..4) 3 or 4, using the function "Compare".
-  // After this sorting, it is impossible to use
-  // the "fNSegments" and "fIndexOfFirstSegment"
-  // of the HitForRec in the first chamber to explore all segments formed with it.
-  // Is this sorting really needed ????
-  if ((Station == 3) || (Station == 4)) (fSegmentsPtr[Station])->Sort();
-  AliDebug(1,Form("Station: %d  NSegments:  %d ", Station, fNSegments[Station]));
-  return;
-}
-
-  //__________________________________________________________________________
-Double_t AliMUONVTrackReconstructor::GetImpactParamFromBendingMomentum(Double_t BendingMomentum) const
-{
-  /// Returns impact parameter at vertex in bending plane (cm),
-  /// from the signed bending momentum "BendingMomentum" in bending plane (GeV/c),
-  /// using simple values for dipole magnetic field.
-  /// The sign of "BendingMomentum" is the sign of the charge.
-  return (-0.0003 * fSimpleBValue * fSimpleBLength * fSimpleBPosition /
-         BendingMomentum);
-}
-
-  //__________________________________________________________________________
-Double_t AliMUONVTrackReconstructor::GetBendingMomentumFromImpactParam(Double_t ImpactParam) const
-{
-  /// Returns signed bending momentum in bending plane (GeV/c),
-  /// the sign being the sign of the charge for particles moving forward in Z,
-  /// from the impact parameter "ImpactParam" at vertex in bending plane (cm),
-  /// using simple values for dipole magnetic field.
-  return (-0.0003 * fSimpleBValue * fSimpleBLength * fSimpleBPosition /
-         ImpactParam);
+  AliDebug(1,Form("Station: %d  NSegments:  %d ", station, segments->GetEntriesFast()));
+  return segments;
 }
 
   //__________________________________________________________________________
 void AliMUONVTrackReconstructor::ValidateTracksWithTrigger(void)
 {
   /// Try to match track from tracking system with trigger track
+  static const Double_t distSigma[3]={1,1,0.02}; // sigma of distributions (trigger-track) X,Y,slopeY
+  
   AliMUONTrack *track;
   AliMUONTrackParam trackParam; 
   AliMUONTriggerTrack *triggerTrack;
@@ -412,9 +308,8 @@ void AliMUONVTrackReconstructor::ValidateTracksWithTrigger(void)
   TClonesArray *recTriggerTracks = fMUONData->RecTriggerTracks();
   
   Bool_t matchTrigger;
-  Double_t distSigma[3]={1,1,0.02}; // sigma of distributions (trigger-track) X,Y,slopeY
   Double_t distTriggerTrack[3];
-  Double_t chi2MatchTrigger, xTrack, yTrack, ySlopeTrack, dTrigTrackMin2, dTrigTrack2 = 0.;
+  Double_t xTrack, yTrack, ySlopeTrack, chi2MatchTrigger, minChi2MatchTrigger, chi2;
   
   track = (AliMUONTrack*) fRecTracksPtr->First();
   while (track) {
@@ -427,19 +322,20 @@ void AliMUONVTrackReconstructor::ValidateTracksWithTrigger(void)
     xTrack = trackParam.GetNonBendingCoor();
     yTrack = trackParam.GetBendingCoor();
     ySlopeTrack = trackParam.GetBendingSlope();
-    dTrigTrackMin2 = 999.;
+    minChi2MatchTrigger = 999.;
   
     triggerTrack = (AliMUONTriggerTrack*) recTriggerTracks->First();
     while(triggerTrack){
       distTriggerTrack[0] = (triggerTrack->GetX11()-xTrack)/distSigma[0];
       distTriggerTrack[1] = (triggerTrack->GetY11()-yTrack)/distSigma[1];
       distTriggerTrack[2] = (TMath::Tan(triggerTrack->GetThetay())-ySlopeTrack)/distSigma[2];
-      dTrigTrack2 = 0.;
-      for (Int_t iVar = 0; iVar < 3; iVar++) dTrigTrack2 += distTriggerTrack[iVar]*distTriggerTrack[iVar];
-      if (dTrigTrack2 < dTrigTrackMin2 && dTrigTrack2 < fMaxSigma2Distance) {
-        dTrigTrackMin2 = dTrigTrack2;
+      chi2 = 0.;
+      for (Int_t iVar = 0; iVar < 3; iVar++) chi2 += distTriggerTrack[iVar]*distTriggerTrack[iVar];
+      chi2 /= 3.; // Normalized Chi2: 3 degrees of freedom (X,Y,slopeY)
+      if (chi2 < minChi2MatchTrigger && chi2 < fMaxNormChi2MatchTrigger) {
+        minChi2MatchTrigger = chi2;
         matchTrigger = kTRUE;
-        chi2MatchTrigger =  dTrigTrack2/3.; // Normalized Chi2, 3 variables (X,Y,slopeY)
+        chi2MatchTrigger = chi2;
       }
       triggerTrack = (AliMUONTriggerTrack*) recTriggerTracks->After(triggerTrack);
     }
@@ -469,7 +365,6 @@ Bool_t AliMUONVTrackReconstructor::MakeTriggerTracks(void)
   AliDebug(1, "Enter MakeTriggerTracks");
   
   TTree* treeR;
-  Int_t nTRentries;
   UChar_t gloTrigPat;
   TClonesArray *localTrigger;
   TClonesArray *globalTrigger;
@@ -482,15 +377,6 @@ Bool_t AliMUONVTrackReconstructor::MakeTriggerTracks(void)
     return kFALSE;
   }
   
-  nTRentries = Int_t(treeR->GetEntries());
-   
-  treeR->GetEvent(0); // only one entry  
-
-  if (!(fMUONData->IsTriggerBranchesInTree())) {
-    AliWarning(Form("Trigger information is not avalaible, nTRentries = %d not equal to 1",nTRentries));
-    return kFALSE;
-  }
-
   fMUONData->SetTreeAddress("TC");
   fMUONData->GetTrigger();
 
index 2fd0f78..58ed418 100644 (file)
@@ -33,8 +33,6 @@ class AliMUONVTrackReconstructor : public TObject {
   Double_t GetBendingResolution(void) const {return fBendingResolution;}
            /// Return chamber resolution (cm) in non-bending plane
   Double_t GetNonBendingResolution(void) const {return fNonBendingResolution;}
-           /// Return chamber thickness in number of radiation lengths
-  Double_t GetChamberThicknessInX0(void) const {return fChamberThicknessInX0;}
 
   // Reconstructed tracks
            /// Return number of reconstructed tracks
@@ -50,10 +48,6 @@ class AliMUONVTrackReconstructor : public TObject {
   virtual void EventDump(void) = 0;  // dump reconstructed event
   void EventDumpTrigger(void);  // dump reconstructed trigger event
   
-  // needed in MakeSegmentsPerStation !!!
-  Double_t GetImpactParamFromBendingMomentum(Double_t BendingMomentum) const;
-  Double_t GetBendingMomentumFromImpactParam(Double_t ImpactParam) const;
-  
           /// Return MUON data
   AliMUONData*  GetMUONData() {return fMUONData;}
 
@@ -68,42 +62,27 @@ class AliMUONVTrackReconstructor : public TObject {
   static const Double_t fgkDefaultMaxBendingMomentum; ///< default max. bending momentum for reconstruction
   static const Double_t fgkDefaultBendingResolution; ///< default bending coordinate resolution for reconstruction 
   static const Double_t fgkDefaultNonBendingResolution; ///< default non bending coordinate resolution for reconstruction
-  static const Double_t fgkDefaultMaxSigma2Distance; ///< default square of max. distance for window size 
-  // Simple magnetic field:
-  // Value taken from macro MUONtracking.C: 0.7 T, hence 7 kG
-  // Length and Position from reco_muon.F, with opposite sign:
-  // Length = ZMAGEND-ZCOIL
-  // Position = (ZMAGEND+ZCOIL)/2
-  // to be ajusted differently from real magnetic field ????
-  static const Double_t fgkDefaultSimpleBValue; ///< default value of magnetic field (dipole)
-  static const Double_t fgkDefaultSimpleBLength; ///< default length of magnetic field (dipole)
-  static const Double_t fgkDefaultSimpleBPosition; ///< default position of magnetic field (dipole)
+  static const Double_t fgkDefaultBendingVertexDispersion; ///< default vertex dispersion in bending plane for reconstruction
+  static const Double_t fgkDefaultNonBendingVertexDispersion; ///< default vertex dispersion in non bending plane for reconstruction
+  static const Double_t fgkDefaultMaxNormChi2MatchTrigger; ///< default maximum normalized chi2 of tracking/trigger track matching
   
   // Parameters for track reconstruction
   Double_t fMinBendingMomentum; ///< minimum value (GeV/c) of momentum in bending plane
   Double_t fMaxBendingMomentum; ///< maximum value (GeV/c) of momentum in bending plane
   Double_t fBendingResolution; ///< chamber resolution (cm) in bending plane
   Double_t fNonBendingResolution; ///< chamber resolution (cm) in non bending plane
-  Double_t fMaxSigma2Distance; ///< maximum square distance in units of the variance (maximum chi2)
-  Double_t fChamberThicknessInX0; ///< chamber thickness in number of radiation lengths
-  Double_t fSimpleBValue; ///< simple magnetic field: value (kG)
-  Double_t fSimpleBLength; ///< simple magnetic field: length (cm)
-  Double_t fSimpleBPosition; ///< simple magnetic field: Z central position (cm)
+  Double_t fBendingVertexDispersion; ///< vextex dispersion (cm) in bending plane
+  Double_t fNonBendingVertexDispersion; ///< vextex dispersion (cm) in non bending plane
+  Double_t fMaxNormChi2MatchTrigger; ///< maximum normalized chi2 of tracking/trigger track matching
   
   Double_t* fSegmentMaxDistBending; ///< maximum distance (cm) for segments in bending plane
   Double_t* fSegmentMaxDistNonBending; ///< maximum distance (cm) for segments in non bending plane
   
-  // Hits for reconstruction (should be in AliMUON ????)
   TClonesArray* fHitsForRecPtr; ///< pointer to the array of hits for reconstruction
   Int_t fNHitsForRec; ///< number of hits for reconstruction
-  // Information per chamber (should be in AliMUONChamber ????)
   Int_t* fNHitsForRecPerChamber; ///< number of HitsForRec
   Int_t* fIndexOfFirstHitForRecPerChamber; ///< index (0...) of first HitForRec
 
-  // Segments inside a station
-  TClonesArray** fSegmentsPtr; ///< array of pointers to the segments for each station
-  Int_t* fNSegments; ///< number of segments for each station
-
   // Reconstructed tracks
   TClonesArray *fRecTracksPtr; ///< pointer to array of reconstructed tracks
   Int_t fNRecTracks; ///< number of reconstructed tracks
@@ -116,14 +95,15 @@ class AliMUONVTrackReconstructor : public TObject {
   AliMUONVTrackReconstructor& operator=(const AliMUONVTrackReconstructor& rhs); ///< assignment operator
   
   void SortHitsForRecWithIncreasingChamber();
-  void MakeSegmentsPerStation(Int_t Station);
+  TClonesArray *MakeSegmentsInStation(Int_t station);
 
   virtual void AddHitsForRecFromRawClusters() = 0;
-  virtual void MakeSegments(void) = 0;
   virtual void MakeTracks(void) = 0;
   virtual void MakeTrackCandidates(void) = 0;
   virtual void FollowTracks(void) = 0;
   virtual void RemoveDoubleTracks(void) = 0;
+  virtual void ExtrapTracksToVertex(void) = 0;
+  virtual void FillMUONTrack(void) = 0;
 
  private:
   
@@ -135,7 +115,6 @@ class AliMUONVTrackReconstructor : public TObject {
   void SetReconstructionParametersToDefaults(void);
   
   void ResetTracks(void);
-  void ResetSegments(void);
   void ResetHitsForRec(void);
   
   void ValidateTracksWithTrigger(void);
@@ -144,6 +123,6 @@ class AliMUONVTrackReconstructor : public TObject {
 
 
   ClassDef(AliMUONVTrackReconstructor, 0) // MUON track reconstructor in ALICE
-    };
+};
        
 #endif
index 6e9da24..a7aba03 100644 (file)
@@ -21,7 +21,6 @@
 #pragma link C++ class AliMUONTriggerTrack+; 
 #pragma link C++ class AliMUONRecoTrack+; 
 #pragma link C++ class AliMUONHitForRec+; 
-#pragma link C++ class AliMUONSegment+;
 #pragma link C++ class AliMUONClusterDrawAZ+;
 #pragma link C++ class AliMUONDetElement+; 
 #pragma link C++ class AliMUONEventRecoCombi+; 
index 94d67a1..8f00ccc 100644 (file)
@@ -17,7 +17,6 @@ SRCS:= AliMUONClusterReconstructor.cxx \
        AliMUONTriggerTrack.cxx \
        AliMUONRecoTrack.cxx \
        AliMUONHitForRec.cxx \
-       AliMUONSegment.cxx \
        AliMUONClusterDrawAZ.cxx \
        AliMUONDetElement.cxx \
        AliMUONEventRecoCombi.cxx \