AliMUONEventReconstructor package:
authorgosset <gosset@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 18 Jul 2000 16:04:06 +0000 (16:04 +0000)
committergosset <gosset@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 18 Jul 2000 16:04:06 +0000 (16:04 +0000)
* a few minor modifications and more comments
* a few corrections
  * right sign for Z of raw clusters
  * right loop over chambers inside station
  * symmetrized covariance matrix for measurements (TrackChi2MCS)
  * right sign of charge in extrapolation (ExtrapToZ)
  * right zEndAbsorber for Branson correction below 3 degrees
* use of TVirtualFitter instead of TMinuit for AliMUONTrack::Fit
* no parameter for AliMUONTrack::Fit() but more fit parameters in Track object

MUON/AliMUONEventReconstructor.cxx
MUON/AliMUONHitForRec.h
MUON/AliMUONTrack.cxx
MUON/AliMUONTrack.h
MUON/AliMUONTrackHit.cxx
MUON/AliMUONTrackHit.h
MUON/AliMUONTrackParam.cxx
MUON/AliMUONTrackParam.h

index b6fa9b4..2023f4b 100644 (file)
@@ -15,6 +15,9 @@
 
 /*
 $Log$
+Revision 1.7  2000/07/03 12:28:06  gosset
+Printout at the right place after extrapolation to vertex
+
 Revision 1.6  2000/06/30 12:01:06  gosset
 Correction for hit search in the right chamber (JPC)
 
@@ -673,9 +676,9 @@ void AliMUONEventReconstructor::AddHitsForRecFromRawClusters(TTree* TR)
       //  original raw cluster
       hitForRec->SetChamberNumber(ch);
       hitForRec->SetHitNumber(iclus);
-      // Z coordinate of the chamber (cm) with sign opposite to GEANT sign
+      // Z coordinate of the chamber (cm)
       // could (should) be more exact from chamber geometry ???? 
-      hitForRec->SetZ(-(&(pMUON->Chamber(ch)))->Z());
+      hitForRec->SetZ((&(pMUON->Chamber(ch)))->Z());
       if (fPrintLevel >= 10) {
        cout << "chamber (0...): " << ch <<
          " raw cluster (0...): " << iclus << endl;
@@ -972,7 +975,7 @@ void AliMUONEventReconstructor::MakeTrackCandidates(void)
       // because candidates with 2 segments have to looked for only once.
       if (begStation == 4)
        nbCan2Seg = MakeTrackCandidatesWithTwoSegments(begSegment);
-      // Look for track candidates with one segments and one point,
+      // 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 ????
@@ -1008,12 +1011,14 @@ void AliMUONEventReconstructor::FollowTracks(void)
     nextTrack = (AliMUONTrack*) fRecTracksPtr->After(track); // prepare next track
     if (fPrintLevel >= 2)
       cout << "FollowTracks: track candidate(0..): " << trackIndex << endl;
-    // Fit track candidate from vertex at X = Y = 0 
-    track->SetFitMCS(0);  // fit without Multiple Scattering
-    track->Fit(track->GetTrackParamAtVertex(), 3);
+    // Fit track candidate
+    track->SetFitMCS(0); // without multiple Coulomb scattering
+    track->SetFitNParam(3); // with 3 parameters (X = Y = 0)
+    track->SetFitStart(0); // from parameters at vertex
+    track->Fit();
     if (fPrintLevel >= 10) {
       cout << "FollowTracks: track candidate(0..): " << trackIndex
-          << " after fit in stations(1..) 4 and 5" << endl;
+          << " after fit in stations(0..) 3 and 4" << endl;
       track->RecursiveDump();
     }
     // Loop over stations(1..) 3, 2 and 1
@@ -1101,7 +1106,7 @@ void AliMUONEventReconstructor::FollowTracks(void)
          extrapHit->SetBendingReso2(extrapSegment->GetBendingCoorReso2());
          extrapHit->SetNonBendingReso2(extrapSegment->GetNonBendingCoorReso2());
          // Loop over hits in the chamber
-         ch = 2 * station + ch;
+         ch = 2 * station + chInStation;
          for (iHit = fIndexOfFirstHitForRecPerChamber[ch];
               iHit < fIndexOfFirstHitForRecPerChamber[ch] +
                 fNHitsForRecPerChamber[ch];
@@ -1120,7 +1125,7 @@ void AliMUONEventReconstructor::FollowTracks(void)
        if (bestHit) {
          // best hit found: add it to track candidate
          track->AddHitForRec(bestHit);
-         // set track parameters at these two TrackHit's
+         // set track parameters at this TrackHit
          track->SetTrackParamAtHit(track->GetNTrackHits() - 1,
                                    &(trackParam[chBestHit]));
          if (fPrintLevel >= 10) {
@@ -1148,17 +1153,21 @@ void AliMUONEventReconstructor::FollowTracks(void)
       // Update track parameters at first track hit (smallest Z)
       trackParam1 = ((AliMUONTrackHit*)
                     (track->GetTrackHitsPtr()->First()))->GetTrackParam();
-      // Track fit from first track hit varying X and Y,
+      // Track fit
       // with multiple Coulomb scattering if all stations
       if (station == 0) track->SetFitMCS(1);
+      // without multiple Coulomb scattering if not all stations
       else track->SetFitMCS(0);
-      track->Fit(trackParam1, 5);
+      track->SetFitNParam(5);  // with 5 parameters (momentum and position)
+      track->SetFitStart(1);  // from parameters at first hit
+      track->Fit();
       if (fPrintLevel >= 10) {
        cout << "FollowTracks: track candidate(0..): " << trackIndex
             << " after fit from station(0..): " << station << " to 4" << endl;
        track->RecursiveDump();
       }
       // Track extrapolation to the vertex through the absorber (Branson)
+      // after going through the first station
       if (station == 0) {
        trackParamVertex = *trackParam1;
        (&trackParamVertex)->ExtrapToVertex();
index c34177c..98deca7 100644 (file)
@@ -26,90 +26,34 @@ class AliMUONHitForRec : public TObject {
   AliMUONHitForRec(AliMUONRawCluster* RawCluster); // Constructor from raw cluster
 
   // Inline functions for Get and Set
-  Double_t GetBendingCoor(void) {
-    // Get fBendingCoor
-    return fBendingCoor;}
-  void SetBendingCoor(Double_t BendingCoor) {
-    // Set fBendingCoor
-    fBendingCoor = BendingCoor;}
-  Double_t GetNonBendingCoor(void) {
-    // Get fNonBendingCoor
-    return fNonBendingCoor;}
-  void SetNonBendingCoor(Double_t NonBendingCoor) {
-    // Set fNonBendingCoor
-    fNonBendingCoor = NonBendingCoor;}
-  Double_t GetZ(void) {
-    // Get fZ
-    return fZ;}
-  void SetZ(Double_t Z) {
-    // Set fZ
-    fZ = Z;}
-  Double_t GetBendingReso2(void) {
-    // Get fBendingReso2
-    return fBendingReso2;}
-  void SetBendingReso2(Double_t BendingReso2) {
-    // Set fBendingReso2
-    fBendingReso2 = BendingReso2;}
-  Double_t GetNonBendingReso2(void) {
-    // Get fNonBendingReso2
-    return fNonBendingReso2;}
-  void SetNonBendingReso2(Double_t NonBendingReso2) {
-    // Set fNonBendingReso2
-    fNonBendingReso2 = NonBendingReso2;}
-  Int_t GetChamberNumber(void) {
-    // Get fChamberNumber
-    return fChamberNumber;}
-  void SetChamberNumber(Int_t ChamberNumber) {
-    // Set fChamberNumber
-    fChamberNumber = ChamberNumber;}
-  Int_t GetHitNumber(void) {
-    // Get fHitNumber
-    return fHitNumber;}
-  void SetHitNumber(Int_t HitNumber) {
-    // Set fHitNumber
-    fHitNumber = HitNumber;}
-  Int_t GetTHTrack(void) {
-    // Get fTHTrack
-    return fTHTrack;}
-  void SetTHTrack(Int_t THTrack) {
-    // Set fTHTrack
-    fTHTrack = THTrack;}
-  Int_t GetGeantSignal(void) {
-    // Get fGeantSignal
-    return fGeantSignal;}
-  void SetGeantSignal(Int_t GeantSignal) {
-    // Set fGeantSignal
-    fGeantSignal = GeantSignal;}
-  Int_t GetIndexOfFirstSegment(void) {
-    // Get fIndexOfFirstSegment
-    return fIndexOfFirstSegment;}
-  void SetIndexOfFirstSegment(Int_t IndexOfFirstSegment) {
-    // Set fIndexOfFirstSegment
-    fIndexOfFirstSegment = IndexOfFirstSegment;}
-  Int_t GetNSegments(void) {
-    // Get fNSegments
-    return fNSegments;}
-  void SetNSegments(Int_t NSegments) {
-    // Set fNSegments
-    fNSegments = NSegments;}
-  AliMUONTrackHit* GetFirstTrackHitPtr(void) {
-    // Get fFirstTrackHitPtr
-    return fFirstTrackHitPtr;}
-  void SetFirstTrackHitPtr(AliMUONTrackHit* FirstTrackHitPtr) {
-    // Set fFirstTrackHitPtr
-    fFirstTrackHitPtr = FirstTrackHitPtr;}
-  AliMUONTrackHit* GetLastTrackHitPtr(void) {
-    // Get fLastTrackHitPtr
-    return fLastTrackHitPtr;}
-  void SetLastTrackHitPtr(AliMUONTrackHit* LastTrackHitPtr) {
-    // Set fLastTrackHitPtr
-    fLastTrackHitPtr = LastTrackHitPtr;}
-  Int_t GetNTrackHits(void) {
-    // Get fNTrackHits
-    return fNTrackHits;}
-  void SetNTrackHits(Int_t NTrackHits) {
-    // Set fNTrackHits
-    fNTrackHits = NTrackHits;}
+  Double_t GetBendingCoor(void) { return fBendingCoor;}
+  void SetBendingCoor(Double_t BendingCoor) { fBendingCoor = BendingCoor;}
+  Double_t GetNonBendingCoor(void) { return fNonBendingCoor;}
+  void SetNonBendingCoor(Double_t NonBendingCoor) { fNonBendingCoor = NonBendingCoor;}
+  Double_t GetZ(void) { return fZ;}
+  void SetZ(Double_t Z) { fZ = Z;}
+  Double_t GetBendingReso2(void) { return fBendingReso2;}
+  void SetBendingReso2(Double_t BendingReso2) { fBendingReso2 = BendingReso2;}
+  Double_t GetNonBendingReso2(void) { return fNonBendingReso2;}
+  void SetNonBendingReso2(Double_t NonBendingReso2) { fNonBendingReso2 = NonBendingReso2;}
+  Int_t GetChamberNumber(void) { return fChamberNumber;}
+  void SetChamberNumber(Int_t ChamberNumber) { fChamberNumber = ChamberNumber;}
+  Int_t GetHitNumber(void) { return fHitNumber;}
+  void SetHitNumber(Int_t HitNumber) { fHitNumber = HitNumber;}
+  Int_t GetTHTrack(void) { return fTHTrack;}
+  void SetTHTrack(Int_t THTrack) { fTHTrack = THTrack;}
+  Int_t GetGeantSignal(void) { return fGeantSignal;}
+  void SetGeantSignal(Int_t GeantSignal) { fGeantSignal = GeantSignal;}
+  Int_t GetIndexOfFirstSegment(void) { return fIndexOfFirstSegment;}
+  void SetIndexOfFirstSegment(Int_t IndexOfFirstSegment) { fIndexOfFirstSegment = IndexOfFirstSegment;}
+  Int_t GetNSegments(void) { return fNSegments;}
+  void SetNSegments(Int_t NSegments) { fNSegments = NSegments;}
+  AliMUONTrackHit* GetFirstTrackHitPtr(void) { return fFirstTrackHitPtr;}
+  void SetFirstTrackHitPtr(AliMUONTrackHit* FirstTrackHitPtr) { fFirstTrackHitPtr = FirstTrackHitPtr;}
+  AliMUONTrackHit* GetLastTrackHitPtr(void) { return fLastTrackHitPtr;}
+  void SetLastTrackHitPtr(AliMUONTrackHit* LastTrackHitPtr) { fLastTrackHitPtr = LastTrackHitPtr;}
+  Int_t GetNTrackHits(void) { return fNTrackHits;}
+  void SetNTrackHits(Int_t NTrackHits) { fNTrackHits = NTrackHits;}
 
 
   Double_t NormalizedChi2WithHitForRec(AliMUONHitForRec* Hit, Double_t Sigma2Cut);
index d1f77f1..8e4f615 100644 (file)
 
 /*
 $Log$
+Revision 1.4  2000/06/30 10:15:48  gosset
+Changes to EventReconstructor...:
+precision fit with multiple Coulomb scattering;
+extrapolation to vertex with Branson correction in absorber (JPC)
+
 Revision 1.3  2000/06/25 13:23:28  hristov
 stdlib.h needed for non-Linux compilation
 
@@ -42,9 +47,9 @@ Addition of files for track reconstruction in C++
 #include <iostream.h>
 
 #include <TClonesArray.h>
-#include <TMinuit.h>
 #include <TMath.h>
 #include <TMatrix.h>
+#include <TVirtualFitter.h>
 
 #include "AliMUONEventReconstructor.h" 
 #include "AliMUONHitForRec.h" 
@@ -53,10 +58,6 @@ Addition of files for track reconstruction in C++
 
 #include <stdlib.h>
 
-// variables to be known from minimization functions
-static AliMUONTrack *trackBeingFitted;
-static AliMUONTrackParam *trackParamBeingFitted;
-
 // 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);
@@ -65,6 +66,8 @@ Double_t MultipleScatteringAngle2(AliMUONTrackHit *TrackHit);
 
 ClassImp(AliMUONTrack) // Class implementation in ROOT context
 
+TVirtualFitter* AliMUONTrack::fgFitter = NULL; 
+
   //__________________________________________________________________________
 AliMUONTrack::AliMUONTrack(AliMUONSegment* BegSegment, AliMUONSegment* EndSegment, AliMUONEventReconstructor* EventReconstructor)
 {
@@ -77,7 +80,10 @@ AliMUONTrack::AliMUONTrack(AliMUONSegment* BegSegment, AliMUONSegment* EndSegmen
   AddSegment(EndSegment); // add hits from EndSegment
   fTrackHitsPtr->Sort(); // sort TrackHits according to increasing Z
   SetTrackParamAtVertex(); // set track parameters at vertex
+  // set fit conditions
   fFitMCS = 0;
+  fFitNParam = 3;
+  fFitStart = 1;
   return;
 }
 
@@ -93,15 +99,20 @@ AliMUONTrack::AliMUONTrack(AliMUONSegment* Segment, AliMUONHitForRec* HitForRec,
   AddHitForRec(HitForRec); // add HitForRec
   fTrackHitsPtr->Sort(); // sort TrackHits according to increasing Z
   SetTrackParamAtVertex(); // set track parameters at vertex
+  // set fit conditions
   fFitMCS = 0;
+  fFitNParam = 3;
+  fFitStart = 1;
   return;
 }
 
+  //__________________________________________________________________________
 AliMUONTrack::AliMUONTrack (const AliMUONTrack& MUONTrack)
 {
 // Dummy copy constructor
 }
 
+  //__________________________________________________________________________
 AliMUONTrack & AliMUONTrack::operator=(const AliMUONTrack& MUONTrack)
 {
 // Dummy assignment operator
@@ -111,26 +122,51 @@ AliMUONTrack & AliMUONTrack::operator=(const AliMUONTrack& MUONTrack)
   //__________________________________________________________________________
 void AliMUONTrack::SetFitMCS(Int_t FitMCS)
 {
-  // Set track fit option with or without multiple Coulomb scattering
+  // Set multiple Coulomb scattering option for track fit "fFitMCS"
   // from "FitMCS" argument: 0 without, 1 with
-  fFitMCS = FitMCS;
+  if ((FitMCS == 0) || (FitMCS == 1)) fFitMCS = FitMCS;
   // better implementation with enum(with, without) ????
+  else {
+    cout << "ERROR in AliMUONTrack::SetFitMCS(FitMCS)" << endl;
+    cout << "FitMCS = " << FitMCS << " is neither 0 nor 1" << endl;
+    exit(0);
+  }
+  return;
+}
+
+  //__________________________________________________________________________
+void AliMUONTrack::SetFitNParam(Int_t FitNParam)
+{
+  // Set number of parameters for track fit "fFitNParam" from "FitNParam":
+  // 3 for momentum, 5 for momentum and position
+  if ((FitNParam == 3) || (FitNParam == 5)) fFitNParam = FitNParam;
+  else {
+    cout << "ERROR in AliMUONTrack::SetFitNParam(FitNParam)" << endl;
+    cout << "FitNParam = " << FitNParam << " is neither 3 nor 5" << endl;
+    exit(0);
+  }
+  return;
+}
+
+  //__________________________________________________________________________
+void AliMUONTrack::SetFitStart(Int_t FitStart)
+{
+  // Set multiple Coulomb scattering option for track fit "fFitStart"
+  // from "FitStart" argument: 0 without, 1 with
+  if ((FitStart == 0) || (FitStart == 1)) fFitStart = FitStart;
+  // better implementation with enum(vertex, firstHit) ????
+  else {
+    cout << "ERROR in AliMUONTrack::SetFitStart(FitStart)" << endl;
+    cout << "FitStart = " << FitStart << " is neither 0 nor 1" << endl;
+    exit(0);
+  }
   return;
 }
 
-// Inline functions for Get and Set: inline removed because it does not work !!!!
-AliMUONTrackParam* AliMUONTrack::GetTrackParamAtVertex(void) {
-  // Get pointer to fTrackParamAtVertex
-  return &fTrackParamAtVertex;}
+  //__________________________________________________________________________
 AliMUONTrackParam* AliMUONTrack::GetTrackParamAtFirstHit(void) {
   // Get pointer to TrackParamAtFirstHit
   return ((AliMUONTrackHit*) (fTrackHitsPtr->First()))->GetTrackParam();}
-TClonesArray* AliMUONTrack::GetTrackHitsPtr(void) {
-  // Get fTrackHitsPtr
-  return fTrackHitsPtr;}
-Int_t AliMUONTrack::GetNTrackHits(void) {
-  // Get fNTrackHits
-  return fNTrackHits;}
 
   //__________________________________________________________________________
 void AliMUONTrack::RecursiveDump(void)
@@ -155,75 +191,76 @@ void AliMUONTrack::RecursiveDump(void)
 }
 
   //__________________________________________________________________________
-void AliMUONTrack::Fit(AliMUONTrackParam *TrackParam, Int_t NParam)
+void AliMUONTrack::Fit()
 {
   // Fit the current track ("this"),
-  // starting with track parameters pointed to by "TrackParam",
-  // and with 3 or 5 parameters ("NParam"):
-  // 3 if one keeps X and Y fixed in "TrackParam",
-  // 5 if one lets them vary.
-  if ((NParam != 3) && (NParam != 5)) {
-    cout << "ERROR in AliMUONTrack::Fit, NParam = " << NParam;
-    cout << " , i.e. neither 3 nor 5 ====> EXIT" << endl;
-    exit(0); // right instruction for exit ????
-  }
-  Int_t error = 0;
+  // with or without multiple Coulomb scattering according to "fFitMCS",
+  // with the number of parameters given by "fFitNParam"
+  // (3 if one keeps X and Y fixed in "TrackParam", 5 if one lets them vary),
+  // starting, according to "fFitStart",
+  // with track parameters at vertex or at the first TrackHit.
+  // "fFitMCS", "fFitNParam" and "fFitStart" have to be set before
+  // by calling the corresponding Set methods.
   Double_t arg[1], benC, errorParam, invBenP, lower, nonBenC, upper, x, y;
-  TString parName;
-  TMinuit *minuit = new TMinuit(5);
-  trackBeingFitted = this; // for the track to be known from the function to minimize
-  trackParamBeingFitted = TrackParam; // for the track parameters to be known from the function to minimize; possible to use only Minuit parameters ????
-  // try to use TVirtualFitter to get this feature automatically !!!!
-  minuit->mninit(5, 10, 7); // sysrd, syswr, syssa: useful ????
+  char parName[50];
+  AliMUONTrackParam *trackParam;
+  // Check if Minuit is initialized...
+  fgFitter = TVirtualFitter::Fitter(this); // add 3 or 5 for the maximum number of parameters ???
+  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 (fFitMCS == 0) minuit->SetFCN(TrackChi2);
-  else minuit->SetFCN(TrackChi2MCS);
-  minuit->SetPrintLevel(1); // More printing !!!!
-  // set first 3 parameters
+  if (fFitMCS == 0) fgFitter->SetFCN(TrackChi2);
+  else fgFitter->SetFCN(TrackChi2MCS);
+  arg[0] = 1;
+  fgFitter->ExecuteCommand("SET PRINT", arg, 1); // More printing !!!!
+  // Parameters according to "fFitStart"
+  // (should be a function to be used at every place where needed ????)
+  if (fFitStart == 0) trackParam = &fTrackParamAtVertex;
+  else trackParam = this->GetTrackParamAtFirstHit();
+  // set first 3 Minuit parameters
   // could be tried with no limits for the search (min=max=0) ????
-  minuit->mnparm(0, "InvBenP",
-                TrackParam->GetInverseBendingMomentum(),
-                0.003, -0.4, 0.4, error);
-  minuit->mnparm(1, "BenS",
-                TrackParam->GetBendingSlope(),
-                0.001, -0.5, 0.5, error);
-  minuit->mnparm(2, "NonBenS",
-                TrackParam->GetNonBendingSlope(),
-                0.001, -0.5, 0.5, error);
-  if (NParam == 5) {
-    // set last 2 parameters (no limits for the search: min=max=0)
-    minuit->mnparm(3, "X",
-                  TrackParam->GetNonBendingCoor(),
-                  0.03, 0.0, 0.0, error);
-    minuit->mnparm(4, "Y",
-                  TrackParam->GetBendingCoor(),
-                  0.10, 0.0, 0.0, error);
+  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 (fFitNParam == 5) {
+    // set last 2 Minuit parameters (no limits for the search: min=max=0)
+    fgFitter->SetParameter(3, "X",
+                          trackParam->GetNonBendingCoor(),
+                          0.03, 0.0, 0.0);
+    fgFitter->SetParameter(4, "Y",
+                          trackParam->GetBendingCoor(),
+                          0.10, 0.0, 0.0);
   }
   // search without gradient calculation in the function
-  minuit->mnexcm("SET NOGRADIENT", arg, 0, error);
+  fgFitter->ExecuteCommand("SET NOGRADIENT", arg, 0);
   // minimization
-  minuit->mnexcm("MINIMIZE", arg, 0, error);
+  fgFitter->ExecuteCommand("MINIMIZE", arg, 0);
   // exit from Minuit
-  minuit->mnexcm("EXIT", arg, 0, error); // necessary ????
-  // print results
-  minuit->mnpout(0, parName, invBenP, errorParam, lower, upper, error);
-  minuit->mnpout(1, parName, benC, errorParam, lower, upper, error);
-  minuit->mnpout(2, parName, nonBenC, errorParam, lower, upper, error);
-  if (NParam == 5) {
-    minuit->mnpout(3, parName, x, errorParam, lower, upper, error);
-    minuit->mnpout(4, parName, y, errorParam, lower, upper, error);
+  fgFitter->ExecuteCommand("EXIT", arg, 0); // necessary ????
+  // get results into "invBenP", "benC", "nonBenC" ("x", "y")
+  fgFitter->GetParameter(0, parName, invBenP, errorParam, lower, upper);
+  fgFitter->GetParameter(1, parName, benC, errorParam, lower, upper);
+  fgFitter->GetParameter(2, parName, nonBenC, errorParam, lower, upper);
+  if (fFitNParam == 5) {
+    fgFitter->GetParameter(3, parName, x, errorParam, lower, upper);
+    fgFitter->GetParameter(4, parName, y, errorParam, lower, upper);
   }
   // result of the fit into track parameters
-  TrackParam->SetInverseBendingMomentum(invBenP);
-  TrackParam->SetBendingSlope(benC);
-  TrackParam->SetNonBendingSlope(nonBenC);
-  if (NParam == 5) {
-    TrackParam->SetNonBendingCoor(x);
-    TrackParam->SetBendingCoor(y);
+  trackParam->SetInverseBendingMomentum(invBenP);
+  trackParam->SetBendingSlope(benC);
+  trackParam->SetNonBendingSlope(nonBenC);
+  if (fFitNParam == 5) {
+    trackParam->SetNonBendingCoor(x);
+    trackParam->SetBendingCoor(y);
   }
-  trackBeingFitted = NULL;
-  delete minuit;
 }
 
   //__________________________________________________________________________
@@ -303,13 +340,17 @@ void TrackChi2(Int_t &NParam, Double_t *Gradient, Double_t &Chi2, Double_t *Para
   // 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;
   AliMUONTrackHit* hit;
   AliMUONTrackParam param1;
   Int_t hitNumber;
   Double_t zHit;
   Chi2 = 0.0; // initialize Chi2
   // copy of track parameters to be fitted
-  param1 = *trackParamBeingFitted;
+  trackBeingFitted = (AliMUONTrack*) AliMUONTrack::Fitter()->GetObjectFit();
+  if (trackBeingFitted->GetFitStart() == 0)
+    param1 = *(trackBeingFitted->GetTrackParamAtVertex());
+  else param1 = *(trackBeingFitted->GetTrackParamAtFirstHit());
   // Minuit parameters to be fitted into this copy
   param1.SetInverseBendingMomentum(Param[0]);
   param1.SetBendingSlope(Param[1]);
@@ -350,10 +391,14 @@ void TrackChi2MCS(Int_t &NParam, Double_t *Gradient, Double_t &Chi2, Double_t *P
   // 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;
   AliMUONTrackParam param1;
   Chi2 = 0.0; // initialize Chi2
   // copy of track parameters to be fitted
-  param1 = *trackParamBeingFitted;
+  trackBeingFitted = (AliMUONTrack*) AliMUONTrack::Fitter()->GetObjectFit();
+  if (trackBeingFitted->GetFitStart() == 0)
+    param1 = *(trackBeingFitted->GetTrackParamAtVertex());
+  else param1 = *(trackBeingFitted->GetTrackParamAtFirstHit());
   // Minuit parameters to be fitted into this copy
   param1.SetInverseBendingMomentum(Param[0]);
   param1.SetBendingSlope(Param[1]);
@@ -363,12 +408,13 @@ void TrackChi2MCS(Int_t &NParam, Double_t *Gradient, Double_t &Chi2, Double_t *P
     param1.SetBendingCoor(Param[4]);
   }
 
-  AliMUONTrackHit* hit, hit1, hit2, hit3;
-  Bool_t GoodDeterminant;
+  AliMUONTrackHit *hit;
+  Bool_t goodDeterminant;
   Int_t hitNumber, hitNumber1, hitNumber2, hitNumber3;
   Double_t zHit[10], paramBendingCoor[10], paramNonBendingCoor[10], ap[10];
   Double_t hitBendingCoor[10], hitNonBendingCoor[10];
   Double_t hitBendingReso2[10], hitNonBendingReso2[10];
+  // dimension 10 in parameter ??? related to AliMUONConstants::NTrackingCh() !!!!
   Int_t numberOfHit = TMath::Min(trackBeingFitted->GetNTrackHits(), 10);
   TMatrix *covBending = new TMatrix(numberOfHit, numberOfHit);
   TMatrix *covNonBending = new TMatrix(numberOfHit, numberOfHit);
@@ -381,59 +427,78 @@ void TrackChi2MCS(Int_t &NParam, Double_t *Gradient, Double_t &Chi2, Double_t *P
     // security against infinite loop ????
     (&param1)->ExtrapToZ(zHit[hitNumber]); // extrapolation
     hit->SetTrackParam(&param1);
-    paramBendingCoor[hitNumber]= (&param1)->GetBendingCoor();
-    paramNonBendingCoor[hitNumber]= (&param1)->GetNonBendingCoor();
-    hitBendingCoor[hitNumber]= hit->GetHitForRecPtr()->GetBendingCoor();
-    hitNonBendingCoor[hitNumber]= hit->GetHitForRecPtr()->GetNonBendingCoor();
-    hitBendingReso2[hitNumber]= hit->GetHitForRecPtr()->GetBendingReso2();
-    hitNonBendingReso2[hitNumber]= hit->GetHitForRecPtr()->GetNonBendingReso2();
+    paramBendingCoor[hitNumber] = (&param1)->GetBendingCoor();
+    paramNonBendingCoor[hitNumber] = (&param1)->GetNonBendingCoor();
+    hitBendingCoor[hitNumber] = hit->GetHitForRecPtr()->GetBendingCoor();
+    hitNonBendingCoor[hitNumber] = hit->GetHitForRecPtr()->GetNonBendingCoor();
+    hitBendingReso2[hitNumber] = hit->GetHitForRecPtr()->GetBendingReso2();
+    hitNonBendingReso2[hitNumber] = hit->GetHitForRecPtr()->GetNonBendingReso2();
     ap[hitNumber] = MultipleScatteringAngle2(hit); // multiple scatt. angle ^2  
   }
 
   // Calculates the covariance matrix
+  // One chamber is taken into account between successive hits.
+  // "ap" should be changed for taking into account the eventual missing hits
+  // by defining an "equivalent" chamber thickness !!!!
   for (hitNumber1 = 0; hitNumber1 < numberOfHit; hitNumber1++) {    
     for (hitNumber2 = hitNumber1; hitNumber2 < numberOfHit; hitNumber2++) {
-      (*covBending)(hitNumber1, hitNumber2) = 0;
-      (*covBending)(hitNumber2, hitNumber1) = 0;
-      if (hitNumber1 == hitNumber2){ // diagonal elements
-       (*covBending)(hitNumber2, hitNumber1) =
-         (*covBending)(hitNumber2, hitNumber1) + hitBendingReso2[hitNumber1];
-      }
-      // Multiple Scattering...  loop on upstream chambers ??
-      for (hitNumber3 = 0; hitNumber3 < hitNumber1; hitNumber3++){     
+      // initialization to 0 (diagonal plus upper triangular part)
+      (*covBending)(hitNumber2, hitNumber1) = 0.0;
+      // contribution from multiple scattering in bending plane:
+      // loop over upstream hits
+      for (hitNumber3 = 0; hitNumber3 < hitNumber1; hitNumber3++) {    
        (*covBending)(hitNumber2, hitNumber1) =
          (*covBending)(hitNumber2, hitNumber1) +
          ((zHit[hitNumber1] - zHit[hitNumber3]) *
           (zHit[hitNumber2] - zHit[hitNumber3]) * ap[hitNumber3]); 
-      }  
-      (*covNonBending)(hitNumber1, hitNumber2) = 0;
+      }
+      // equal contribution from multiple scattering in non bending plane
       (*covNonBending)(hitNumber2, hitNumber1) =
        (*covBending)(hitNumber2, hitNumber1);
-      if (hitNumber1 == hitNumber2) {  // diagonal elements
+      if (hitNumber1 == hitNumber2) {
+       // Diagonal elements: add contribution from position measurements
+       // in bending plane
+       (*covBending)(hitNumber2, hitNumber1) =
+         (*covBending)(hitNumber2, hitNumber1) + hitBendingReso2[hitNumber1];
+       // and in non bending plane
        (*covNonBending)(hitNumber2, hitNumber1) =
-         (*covNonBending)(hitNumber2, hitNumber1) -
-         hitBendingReso2[hitNumber1] + hitNonBendingReso2[hitNumber1] ;
-      }      
-    }
-  }
+         (*covNonBending)(hitNumber2, hitNumber1) + hitNonBendingReso2[hitNumber1];
+      }
+      else {
+       // Non diagonal elements: symmetrization
+       // for bending plane
+       (*covBending)(hitNumber1, hitNumber2) =
+         (*covBending)(hitNumber2, hitNumber1);
+       // and non bending plane
+       (*covNonBending)(hitNumber1, hitNumber2) =
+         (*covNonBending)(hitNumber2, hitNumber1);
+      }
+    } // for (hitNumber2 = hitNumber1;...
+  } // for (hitNumber1 = 0;...
 
   // Inverts covariance matrix 
-  GoodDeterminant = kTRUE;
+  goodDeterminant = kTRUE;
+  // check whether the Invert method returns flag if matrix cannot be inverted,
+  // and do not calculate the Determinant in that case !!!!
   if (covBending->Determinant() != 0) {
     covBending->Invert();
   } else {
-    GoodDeterminant = kFALSE;
+    goodDeterminant = kFALSE;
     cout << "Warning in ChiMCS  Determinant Bending=0: " << endl;  
   }
-  if (covNonBending->Determinant() != 0){
+  if (covNonBending->Determinant() != 0) {
     covNonBending->Invert();
   } else {
-    GoodDeterminant = kFALSE;
+    goodDeterminant = kFALSE;
     cout << "Warning in ChiMCS  Determinant non Bending=0: " << endl;  
   }
+
+  // It would be worth trying to calculate the inverse of the covariance matrix
+  // only once per fit, since it cannot change much in principle,
+  // and it would save a lot of computing time !!!!
   
   // Calculates Chi2
-  if (GoodDeterminant) { // with Multiple Scattering if inversion correct
+  if (goodDeterminant) { // with Multiple Scattering if inversion correct
     for (hitNumber1=0; hitNumber1 < numberOfHit ; hitNumber1++){ 
       for (hitNumber2=0; hitNumber2 < numberOfHit; hitNumber2++){
        Chi2 = Chi2 +
@@ -469,7 +534,9 @@ Double_t MultipleScatteringAngle2(AliMUONTrackHit *TrackHit)
   // at TrackHit pointed to by "TrackHit"
   Double_t slopeBending, slopeNonBending, radiationLength, inverseBendingMomentum2, inverseTotalMomentum2;
   Double_t varMultipleScatteringAngle;
+  AliMUONTrack *trackBeingFitted = (AliMUONTrack*) AliMUONTrack::Fitter()->GetObjectFit();
   AliMUONTrackParam *param = TrackHit->GetTrackParam();
+  // Better implementation in AliMUONTrack class ????
   slopeBending = param->GetBendingSlope();
   slopeNonBending = param->GetNonBendingSlope();
   // thickness in radiation length for the current track,
@@ -481,9 +548,10 @@ Double_t MultipleScatteringAngle2(AliMUONTrackHit *TrackHit)
   inverseBendingMomentum2 = 
     param->GetInverseBendingMomentum() * param->GetInverseBendingMomentum();
   inverseTotalMomentum2 =
-    inverseBendingMomentum2 * (1.0 + slopeBending*slopeBending) /
+    inverseBendingMomentum2 * (1.0 + slopeBending * slopeBending) /
     (1.0 + slopeBending *slopeBending + slopeNonBending * slopeNonBending); 
   varMultipleScatteringAngle = 0.0136 * (1.0 + 0.038 * TMath::Log(radiationLength));
+  // The velocity is assumed to be 1 !!!!
   varMultipleScatteringAngle = inverseTotalMomentum2 * radiationLength *
     varMultipleScatteringAngle * varMultipleScatteringAngle;
   return varMultipleScatteringAngle;
index 4ada4a7..acfc47d 100644 (file)
@@ -7,10 +7,10 @@
 
 #include <TROOT.h>
 #include <TClonesArray.h>
-#include "AliMUONTrackHit.h"
 #include "AliMUONTrackParam.h"
 
 class TClonesArray;
+class TVirtualFitter;
 class AliMUONEventReconstructor;
 class AliMUONHitForRec;
 class AliMUONSegment;
@@ -20,40 +20,47 @@ class AliMUONTrack : public TObject {
   AliMUONTrack(){
     // Constructor
     ;} // Constructor
-  virtual ~AliMUONTrack(){
-    // Destructor
-    ;} // Destructor
+  virtual ~AliMUONTrack() {;} // Destructor
   AliMUONTrack (const AliMUONTrack& AliMUONTrack); // copy constructor
   AliMUONTrack& operator=(const AliMUONTrack& AliMUONTrack); // assignment operator
 
   AliMUONTrack(AliMUONSegment* BegSegment, AliMUONSegment* EndSegment, AliMUONEventReconstructor* EventReconstructor); // Constructor from two Segment's
   AliMUONTrack(AliMUONSegment* Segment, AliMUONHitForRec* HitForRec, AliMUONEventReconstructor* EventReconstructor); // Constructor from one Segment and one HitForRec
 
-  AliMUONEventReconstructor* GetEventReconstructor(void) {return fEventReconstructor;};
-  AliMUONTrackParam* GetTrackParamAtVertex(void);
+  AliMUONEventReconstructor* GetEventReconstructor(void) {return fEventReconstructor;}
+  AliMUONTrackParam* GetTrackParamAtVertex(void) { return &fTrackParamAtVertex;}
   void SetTrackParamAtVertex(void); // Set track parameters at vertex from last stations 4 & 5
-  void SetTrackParamAtVertex(AliMUONTrackParam* TrackParam) {fTrackParamAtVertex = *TrackParam;};
+  void SetTrackParamAtVertex(AliMUONTrackParam* TrackParam) {fTrackParamAtVertex = *TrackParam;}
 
-  TClonesArray* GetTrackHitsPtr(void);
-  Int_t GetNTrackHits(void);
-  Int_t GetFitMCS(void) {return fFitMCS;};
+  TClonesArray* GetTrackHitsPtr(void) { return fTrackHitsPtr;}
+  Int_t GetNTrackHits(void) { return fNTrackHits;}
+  Int_t GetFitMCS(void) {return fFitMCS;}
+  Int_t GetFitNParam(void) {return fFitNParam;}
+  Int_t GetFitStart(void) {return fFitStart;}
   void SetFitMCS(Int_t FitMCS);
+  void SetFitNParam(Int_t FitNParam);
+  void SetFitStart(Int_t FitStart);
 
   AliMUONTrackParam* GetTrackParamAtFirstHit(void);
 
   void RecursiveDump(void); // Recursive dump (with track hits)
-  void Fit(AliMUONTrackParam *TrackParam, Int_t NParam); // Fit
+  void Fit(); // Fit
   void AddSegment(AliMUONSegment* Segment); // Add Segment
   void AddHitForRec(AliMUONHitForRec* HitForRec); // Add HitForRec
   void SetTrackParamAtHit(Int_t indexHit, AliMUONTrackParam *TrackParam);
 
+  static TVirtualFitter* AliMUONTrack::Fitter(void) {return fgFitter;}
+
  protected:
  private:
+  static TVirtualFitter* fgFitter; // Pointer to track fitter
   AliMUONEventReconstructor* fEventReconstructor; // Pointer to EventReconstructor
   AliMUONTrackParam fTrackParamAtVertex; // Track parameters at vertex
   TClonesArray *fTrackHitsPtr; // Pointer to array of TrackHit's
   Int_t fNTrackHits; // Number of TrackHit's
   Int_t fFitMCS; // 0(1) for fit without(with) multiple Coulomb scattering
+  Int_t fFitNParam; // 3(5) for fit with 3(5) parameters
+  Int_t fFitStart; // 0 or 1 for fit starting from parameters at vertex (0) or at first TrackHit(1)
   
   ClassDef(AliMUONTrack, 1) // Reconstructed track in ALICE dimuon spectrometer
     };
index 024cf11..fa14a66 100644 (file)
@@ -15,6 +15,9 @@
 
 /*
 $Log$
+Revision 1.3  2000/06/25 13:06:39  hristov
+Inline functions moved from *.cxx to *.h files instead of forward declarations
+
 Revision 1.2  2000/06/15 07:58:49  morsch
 Code from MUON-dev joined
 
@@ -59,12 +62,13 @@ AliMUONTrackHit::AliMUONTrackHit(AliMUONHitForRec* Hit)
   Hit->SetNTrackHits(Hit->GetNTrackHits() + 1);
 }
 
-
+  //__________________________________________________________________________
 AliMUONTrackHit::AliMUONTrackHit (const AliMUONTrackHit& MUONTrackHit)
 {
 // Dummy copy constructor
 }
 
+  //__________________________________________________________________________
 AliMUONTrackHit & AliMUONTrackHit::operator=(const AliMUONTrackHit& MUONTrackHit)
 {
 // Dummy assignment operator
@@ -73,6 +77,30 @@ AliMUONTrackHit & AliMUONTrackHit::operator=(const AliMUONTrackHit& MUONTrackHit
 
 
   //__________________________________________________________________________
+AliMUONTrackHit::~AliMUONTrackHit()
+{
+  // Destructor
+//   AliMUONHitForRec * hit; // pointer to HitForRec
+//   // remove current TrackHit in HitForRec links
+//   if (this == hit->GetFirstTrackHitPtr())
+//     hit->SetFirstTrackHitPtr(fNextTrackHitWithSameHitForRec); // if first
+//   if (this == hit->GetLastTrackHitPtr())
+//     hit->SetLastTrackHitPtr(fPrevTrackHitWithSameHitForRec); // if last
+//   hit->SetNTrackHits(hit->GetNTrackHits() - 1); // decrement NTrackHits
+//   // update link to next TrackHit of previous TrackHit
+//   if (fPrevTrackHitWithSameHitForRec != NULL)
+//     fPrevTrackHitWithSameHitForRec->
+//       SetNextTrackHitWithSameHitForRec(fNextTrackHitWithSameHitForRec);
+//   // update link to previous TrackHit of next TrackHit
+//   if (fNextTrackHitWithSameHitForRec)
+//     fNextTrackHitWithSameHitForRec->
+//       SetPrevTrackHitWithSameHitForRec(fPrevTrackHitWithSameHitForRec);
+  // to be checked thoroughly !!!!
+  // with Root counter of AliMUONTrackHit objects,
+  // with loop over all these links after the update
+}
+
+  //__________________________________________________________________________
 Int_t AliMUONTrackHit::Compare(TObject* TrackHit)
 {
   // "Compare" function to sort with increasing Z.
index df2996a..f373791 100644 (file)
@@ -15,9 +15,7 @@ class AliMUONTrackHit : public TObject {
   AliMUONTrackHit(){
     // Constructor
     ;} // Constructor
-  virtual ~AliMUONTrackHit(){
-    // Destructor
-    ;} // Destructor
+  virtual ~AliMUONTrackHit(); // Destructor
   AliMUONTrackHit (const AliMUONTrackHit& AliMUONTrackHit); // copy constructor
   AliMUONTrackHit& operator=(const AliMUONTrackHit& AliMUONTrackHit); // assignment operator
   AliMUONTrackHit(AliMUONHitForRec* Hit); // Constructor from one HitForRec
@@ -41,6 +39,9 @@ class AliMUONTrackHit : public TObject {
 
  protected:
  private:
+  void SetNextTrackHitWithSameHitForRec(AliMUONTrackHit *Next) {fNextTrackHitWithSameHitForRec = Next;}
+  void SetPrevTrackHitWithSameHitForRec(AliMUONTrackHit *Prev) {fPrevTrackHitWithSameHitForRec = Prev;}
+
   AliMUONTrackParam fTrackParam; // Track parameters
   AliMUONHitForRec *fHitForRecPtr; // Pointer to HitForRec
   AliMUONTrackHit *fNextTrackHitWithSameHitForRec; // Pointer to next track hit with same HitForRec
index 10a3375..64898d6 100644 (file)
@@ -15,6 +15,9 @@
 
 /*
 $Log$
+Revision 1.4  2000/07/03 07:53:31  morsch
+Double declaration problem on HP solved.
+
 Revision 1.3  2000/06/30 10:15:48  gosset
 Changes to EventReconstructor...:
 precision fit with multiple Coulomb scattering;
@@ -120,9 +123,10 @@ void AliMUONTrackParam::ExtrapToZ(Double_t Z)
   // For use of reco_ghelix...: invert X and Y, PX/PTOT and PY/PTOT !!!!
   temp = vGeant3[0]; vGeant3[0] = vGeant3[1]; vGeant3[1] = temp;
   temp = vGeant3[3]; vGeant3[3] = vGeant3[4]; vGeant3[4] = temp;
-  // charge must be changed with momentum for backward motion
-  Double_t charge =
-    forwardBackward * TMath::Sign(Double_t(1.0), this->fInverseBendingMomentum);
+  // sign of charge (sign of fInverseBendingMomentum if forward motion)
+  // must be also changed if backward extrapolation
+  Double_t chargeExtrap = forwardBackward *
+    TMath::Sign(Double_t(1.0), this->fInverseBendingMomentum);
   Double_t stepLength = 6.0; // in parameter ????
   // Extrapolation loop
   stepNumber = 0;
@@ -131,7 +135,7 @@ void AliMUONTrackParam::ExtrapToZ(Double_t Z)
     stepNumber++;
     // call Geant3 "ghelix" subroutine through a copy in "reco_muon.F":
     // the true function should be called, but how ???? and remove prototyping ...
-    reco_ghelix(charge, stepLength, vGeant3, vGeant3New);
+    reco_ghelix(chargeExtrap, stepLength, vGeant3, vGeant3New);
     if ((forwardBackward * (vGeant3New[2] - Z)) > 0.0) break; // one is beyond Z
     // better use TArray ????
     for (iGeant3 = 0; iGeant3 < 7; iGeant3++)
@@ -160,12 +164,14 @@ void AliMUONTrackParam::ExtrapToZ(Double_t Z)
   vGeant3[2] = Z; // Z
   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
   vGeant3[5] =
     1.0 / TMath::Sqrt(1.0 + xPrimeI * xPrimeI + yPrimeI * yPrimeI); // PZ/PTOT
   vGeant3[3] = xPrimeI * vGeant3[5]; // PX/PTOT
   vGeant3[4] = yPrimeI * vGeant3[5]; // PY/PTOT
-  // Track parameters from Geant3 parameters
-  GetFromGeant3Parameters(vGeant3, charge);
+  // Track parameters from Geant3 parameters,
+  // with charge back for forward motion
+  GetFromGeant3Parameters(vGeant3, chargeExtrap * forwardBackward);
 }
 
   //__________________________________________________________________________
@@ -193,7 +199,8 @@ void AliMUONTrackParam::SetGeant3Parameters(Double_t *VGeant3, Double_t ForwardB
 void AliMUONTrackParam::GetFromGeant3Parameters(Double_t *VGeant3, Double_t Charge)
 {
   // Get track parameters in current AliMUONTrackParam
-  // from Geant3 parameters pointed to by "VGeant3".
+  // from Geant3 parameters pointed to by "VGeant3",
+  // assumed to be calculated for forward motion in Z.
   // "InverseBendingMomentum" is signed with "Charge".
   this->fNonBendingCoor = VGeant3[0]; // X
   this->fBendingCoor = VGeant3[1]; // Y
@@ -243,7 +250,7 @@ void AliMUONTrackParam::ExtrapToVertex()
   // 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 
+  // Changes parameters according to Branson correction through the absorber 
   
   Double_t zAbsorber = 503.0; // to be coherent with the Geant absorber geometry !!!!
   // Extrapolates track parameters upstream to the "Z" end of the front absorber
@@ -261,6 +268,8 @@ void AliMUONTrackParam::BransonCorrection()
   Double_t  pYZ, pX, pY, pZ, pTotal, xEndAbsorber, yEndAbsorber, radiusEndAbsorber2, pT, theta;
   Int_t sign;
   // Would it be possible to calculate all that from Geant configuration ????
+  // and to get the Branson parameters from a function in ABSO module ????
+  // with an eventual contribution from other detectors like START ????
   // Radiation lengths outer part theta > 3 degres
   static Double_t x01[9] = { 18.8,    // C (cm)
                             10.397,   // Concrete (cm)
@@ -324,7 +333,7 @@ void AliMUONTrackParam::BransonCorrection()
     zEndAbsorber = z1[9];
     zBP = zBP1;
   } else {
-    zEndAbsorber = z2[2];
+    zEndAbsorber = z2[3];
     zBP = zBP2;
   }
 
@@ -365,6 +374,7 @@ Double_t AliMUONTrackParam::TotalMomentumEnergyLoss(Double_t rLimit, Double_t pT
   Double_t radiusEndAbsorber2 =
     xEndAbsorber *xEndAbsorber + yEndAbsorber * yEndAbsorber;
   // Parametrization to be redone according to change of absorber material ????
+  // See remark in function BransonCorrection !!!!
   // The name is not so good, and there are many arguments !!!!
   if (radiusEndAbsorber2 < rLimit * rLimit) {
     if (pTotal < 15) {
index 00577f7..7d5d012 100644 (file)
@@ -41,7 +41,7 @@ class AliMUONTrackParam : public TObject {
 
  protected:
  private:
-  Double_t fInverseBendingMomentum; // Inverse bending momentum (GeV/c ** -1)
+  Double_t fInverseBendingMomentum; // Inverse bending momentum (GeV/c ** -1) times the charge (assumed forward motion)
   Double_t fBendingSlope; // Bending slope (cm ** -1)
   Double_t fNonBendingSlope; // Non bending slope (cm ** -1)
   Double_t fZ; // Z coordinate (cm)