Modifications used for addendum to Dimuon TDR (JP Cussonneau):
authorgosset <gosset@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 8 Jan 2001 11:01:02 +0000 (11:01 +0000)
committergosset <gosset@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 8 Jan 2001 11:01:02 +0000 (11:01 +0000)
*. MaxBendingMomentum to make both a segment and a track (default 500)
*. MaxChi2 per degree of freedom to make a track (default 100)
*. MinBendingMomentum used also to make a track
   and not only a segment (default 3)
*. wider roads for track search in stations 1 to 3
*. extrapolation to actual Z instead of Z(chamber) in FollowTracks
*. in track fit:
   - limits on parameters X and Y (+/-500)
   - covariance matrices in double precision
   - normalization of covariance matrices before inversion
   - suppression of Minuit printouts
*. correction against memory leak (delete extrapHit) in FollowTracks
*. RMax to 10 degrees with Z(chamber) instead of fixed values;
   RMin and Rmax cuts suppressed in NewHitForRecFromGEANT,
   because useless with realistic geometry

MUON/AliMUONEventReconstructor.cxx
MUON/AliMUONEventReconstructor.h
MUON/AliMUONSegment.cxx
MUON/AliMUONSegment.h
MUON/AliMUONTrack.cxx

index fb97ac1924757d4a8edef8e93aaa56673ef24390..67bd12794f0fb96e809fca468d5fc8eec33536b7 100644 (file)
 
 /*
 $Log$
+Revision 1.19  2000/11/20 11:24:10  gosset
+New package for reconstructed tracks (A. Gheata):
+* tree output in the file "tree_reco.root"
+* display events and make histograms from this file
+
 Revision 1.18  2000/10/26 12:47:03  gosset
 Real distance between chambers of each station taken into account
 for the reconstruction parameters "fSegmentMaxDistBending[5]"
@@ -139,6 +144,8 @@ Addition of files for track reconstruction in C++
 
 //************* Defaults parameters for reconstruction
 static const Double_t kDefaultMinBendingMomentum = 3.0;
+static const Double_t kDefaultMaxBendingMomentum = 500.0;
+static const Double_t kDefaultMaxChi2 = 100.0;
 static const Double_t kDefaultMaxSigma2Distance = 16.0;
 static const Double_t kDefaultBendingResolution = 0.01;
 static const Double_t kDefaultNonBendingResolution = 0.144;
@@ -239,6 +246,8 @@ void AliMUONEventReconstructor::SetReconstructionParametersToDefaults(void)
   // Set reconstruction parameters to default values
   // Would be much more convenient with a structure (or class) ????
   fMinBendingMomentum = kDefaultMinBendingMomentum;
+  fMaxBendingMomentum = kDefaultMaxBendingMomentum;
+  fMaxChi2 = kDefaultMaxChi2;
   fMaxSigma2Distance = kDefaultMaxSigma2Distance;
 
   AliMUON *pMUON = (AliMUON*) gAlice->GetModule("MUON");
@@ -251,14 +260,10 @@ void AliMUONEventReconstructor::SetReconstructionParametersToDefaults(void)
     if (ch < 4) fRMin[ch] = TMath::Abs((&(pMUON->Chamber(ch)))->Z()) *
                  2.0 * TMath::Pi() / 180.0;
     else fRMin[ch] = 30.0;
+    // maximum radius at 10 degrees and Z of chamber
+    fRMax[ch] = TMath::Abs((&(pMUON->Chamber(ch)))->Z()) *
+                 10.0 * TMath::Pi() / 180.0;
   }
-  // maximum radius
-  // like in TRACKF_STAT (10 degrees ????)
-  fRMax[0] = fRMax[1] = 91.5;
-  fRMax[2] = fRMax[3] = 122.5;
-  fRMax[4] = fRMax[5] = 158.3;
-  fRMax[6] = fRMax[7] = 260.0;
-  fRMax[8] = fRMax[9] = 260.0;
 
   // ******** Parameters for making segments
   // should be parametrized ????
@@ -601,7 +606,8 @@ AliMUONHitForRec* AliMUONEventReconstructor::NewHitForRecFromGEANT(AliMUONHit* H
   bendCoor = Hit->Y();
   nonBendCoor = Hit->X();
   radius = TMath::Sqrt((bendCoor * bendCoor) + (nonBendCoor * nonBendCoor));
-  if ((radius < fRMin[chamber]) || (radius > fRMax[chamber])) return NULL;
+  // This cut is not needed with a realistic chamber geometry !!!!
+//   if ((radius < fRMin[chamber]) || (radius > fRMax[chamber])) return NULL;
   // new AliMUONHitForRec from GEANT hit and increment number of AliMUONHitForRec's
   hitForRec = new ((*fHitsForRecPtr)[fNHitsForRec]) AliMUONHitForRec(Hit);
   fNHitsForRec++;
@@ -609,6 +615,9 @@ AliMUONHitForRec* AliMUONEventReconstructor::NewHitForRecFromGEANT(AliMUONHit* H
   hitForRec->SetBendingCoor(bendCoor + gRandom->Gaus(0., fBendingResolution));
   hitForRec->SetNonBendingCoor(nonBendCoor
                               + gRandom->Gaus(0., fNonBendingResolution));
+//   // !!!! without smearing
+//   hitForRec->SetBendingCoor(bendCoor);
+//   hitForRec->SetNonBendingCoor(nonBendCoor);
   // more information into HitForRec
   //  resolution: angular effect to be added here ????
   hitForRec->SetBendingReso2(fBendingResolution * fBendingResolution);
@@ -795,7 +804,7 @@ void AliMUONEventReconstructor::MakeSegmentsPerStation(Int_t Station)
   AliMUONSegment *segment;
   Bool_t last2st;
   Double_t bendingSlope, distBend, distNonBend, extBendCoor, extNonBendCoor,
-      impactParam = 0., maxImpactParam = 0.; // =0 to avoid compilation warnings.
+      impactParam = 0., maxImpactParam = 0., minImpactParam = 0.; // =0 to avoid compilation warnings.
   AliMUON *pMUON  = (AliMUON*) gAlice->GetModule("MUON"); // necessary ????
   if (fPrintLevel >= 1)
     cout << "enter MakeSegmentsPerStation (0...) " << Station << endl;
@@ -809,6 +818,9 @@ void AliMUONEventReconstructor::MakeSegmentsPerStation(Int_t Station)
     // 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
@@ -852,7 +864,8 @@ void AliMUONEventReconstructor::MakeSegmentsPerStation(Int_t Station)
       // Conditions "distBend" and "impactParam" correlated for these stations ????
       if ((distBend < fSegmentMaxDistBending[Station]) &&
          (distNonBend < fSegmentMaxDistNonBending[Station]) &&
-         (!last2st || (impactParam < maxImpactParam))) {
+         (!last2st || (impactParam < maxImpactParam)) &&
+         (!last2st || (impactParam > minImpactParam))) {
        // make new segment
        segment = new ((*fSegmentsPtr[Station])[segmentIndex])
          AliMUONSegment(hit1Ptr, hit2Ptr);
@@ -1061,13 +1074,14 @@ void AliMUONEventReconstructor::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;
+  AliMUONHitForRec *bestHit, *extrapHit, *extrapCorrHit, *hit;
+  AliMUONSegment *bestSegment, *extrapSegment, *extrapCorrSegment, *segment;
   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.;
   AliMUON *pMUON = (AliMUON*) gAlice->GetModule("MUON"); // necessary ????
   // local maxSigma2Distance, for easy increase in testing
   maxSigma2Distance = fMaxSigma2Distance;
@@ -1102,6 +1116,7 @@ void AliMUONEventReconstructor::FollowTracks(void)
       // extrapolation to station
       trackParam1->ExtrapToStation(station, trackParam);
       extrapSegment = new AliMUONSegment(); //  empty segment
+      extrapCorrSegment = new AliMUONSegment(); //  empty corrected segment
       // multiple scattering factor corresponding to one chamber
       // and momentum in bending plane (not total)
       mcsFactor = 0.0136 * trackParam1->GetInverseBendingMomentum();
@@ -1121,6 +1136,14 @@ void AliMUONEventReconstructor::FollowTracks(void)
       extrapSegment->UpdateFromStationTrackParam
        (trackParam, mcsFactor, dZ1, dZ2, dZ3, station,
         trackParam1->GetInverseBendingMomentum());
+      // same thing for corrected segment
+      // better to use copy constructor, after checking that it works properly !!!!
+      extrapCorrSegment->SetBendingCoorReso2(fBendingResolution * fBendingResolution);
+      extrapCorrSegment->
+       SetNonBendingCoorReso2(fNonBendingResolution * fNonBendingResolution);
+      extrapCorrSegment->UpdateFromStationTrackParam
+       (trackParam, mcsFactor, dZ1, dZ2, dZ3, station,
+        trackParam1->GetInverseBendingMomentum());
       bestChi2 = 5.0;
       bestSegment = NULL;
       if (fPrintLevel >= 10) {
@@ -1134,7 +1157,20 @@ void AliMUONEventReconstructor::FollowTracks(void)
        // multiple scattering ????
        // separation in 2 functions: Segment and HitForRec ????
        segment = (AliMUONSegment*) ((*fSegmentsPtr[station])[iSegment]);
-       chi2 = segment->NormalizedChi2WithSegment(extrapSegment, maxSigma2Distance);
+       // correction of corrected segment (fBendingCoor and fNonBendingCoor)
+       // according to real Z value of "segment" and slopes of "extrapSegment"
+       extrapCorrSegment->
+         SetBendingCoor(extrapSegment->GetBendingCoor() +
+                        extrapSegment->GetBendingSlope() *
+                        (segment->GetHitForRec1()->GetZ() -
+                         (&(pMUON->Chamber(2 * station)))->Z()));
+       extrapCorrSegment->
+         SetNonBendingCoor(extrapSegment->GetNonBendingCoor() +
+                           extrapSegment->GetNonBendingSlope() *
+                           (segment->GetHitForRec1()->GetZ() -
+                            (&(pMUON->Chamber(2 * station)))->Z()));
+       chi2 = segment->
+         NormalizedChi2WithSegment(extrapCorrSegment, maxSigma2Distance);
        if (chi2 < bestChi2) {
          // update best Chi2 and Segment if better found
          bestSegment = segment;
@@ -1159,6 +1195,7 @@ void AliMUONEventReconstructor::FollowTracks(void)
        // should consider all possibilities ????
        // multiple scattering ???? do about like for extrapSegment !!!!
        extrapHit = new AliMUONHitForRec(); //  empty hit
+       extrapCorrHit = new AliMUONHitForRec(); //  empty corrected hit
        bestChi2 = 3.0;
        bestHit = NULL;
        if (fPrintLevel >= 10) {
@@ -1176,6 +1213,14 @@ void AliMUONEventReconstructor::FollowTracks(void)
          // resolutions from "extrapSegment"
          extrapHit->SetBendingReso2(extrapSegment->GetBendingCoorReso2());
          extrapHit->SetNonBendingReso2(extrapSegment->GetNonBendingCoorReso2());
+         // same things for corrected hit
+         // better to use copy constructor, after checking that it works properly !!!!
+         extrapCorrHit->
+           SetBendingCoor((&(trackParam[chInStation]))->GetBendingCoor());
+         extrapCorrHit->
+           SetNonBendingCoor((&(trackParam[chInStation]))->GetNonBendingCoor());
+         extrapHit->SetBendingReso2(extrapSegment->GetBendingCoorReso2());
+         extrapHit->SetNonBendingReso2(extrapSegment->GetNonBendingCoorReso2());
          // Loop over hits in the chamber
          ch = 2 * station + chInStation;
          for (iHit = fIndexOfFirstHitForRecPerChamber[ch];
@@ -1183,6 +1228,18 @@ void AliMUONEventReconstructor::FollowTracks(void)
                 fNHitsForRecPerChamber[ch];
               iHit++) {
            hit = (AliMUONHitForRec*) ((*fHitsForRecPtr)[iHit]);
+           // correction of corrected hit (fBendingCoor and fNonBendingCoor)
+           // according to real Z value of "hit" and slopes of right "trackParam"
+           extrapCorrHit->
+             SetBendingCoor((&(trackParam[chInStation]))->GetBendingCoor() +
+                            (&(trackParam[chInStation]))->GetBendingSlope() *
+                            (hit->GetZ() -
+                             (&(trackParam[chInStation]))->GetZ()));
+           extrapCorrHit->
+             SetNonBendingCoor((&(trackParam[chInStation]))->GetNonBendingCoor() +
+                               (&(trackParam[chInStation]))->GetNonBendingSlope() *
+                               (hit->GetZ() -
+                                (&(trackParam[chInStation]))->GetZ()));
            // condition for hit not already in segment ????
            chi2 = hit->NormalizedChi2WithHitForRec(extrapHit, maxSigma2Distance);
            if (chi2 < bestChi2) {
@@ -1210,16 +1267,31 @@ void AliMUONEventReconstructor::FollowTracks(void)
          // and corresponding TrackHit's, ...
          track->Remove();
          delete extrapSegment;
+         delete extrapCorrSegment;
+         delete extrapHit;
+         delete extrapCorrHit;
          break; // stop the search for this candidate:
          // exit from the loop over station
        }
+       delete extrapHit;
+       delete extrapCorrHit;
       }
       delete extrapSegment;
+      delete extrapCorrSegment;
       // Sort track hits according to increasing Z
       track->GetTrackHitsPtr()->Sort();
       // Update track parameters at first track hit (smallest Z)
       trackParam1 = ((AliMUONTrackHit*)
                     (track->GetTrackHitsPtr()->First()))->GetTrackParam();
+      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)) {
+       track->Remove();
+       break; // stop the search for this candidate:
+       // exit from the loop over station 
+      }
       // Track fit
       // with multiple Coulomb scattering if all stations
       if (station == 0) track->SetFitMCS(1);
@@ -1228,6 +1300,18 @@ void AliMUONEventReconstructor::FollowTracks(void)
       track->SetFitNParam(5);  // with 5 parameters (momentum and position)
       track->SetFitStart(1);  // from parameters at first hit
       track->Fit();
+      Double_t numberOfDegFree = (2.0 * track->GetNTrackHits() - 5);
+      if (numberOfDegFree > 0) {
+        chi2Norm =  track->GetFitFMin() / numberOfDegFree;
+      } else {
+       chi2Norm = 1.e10;
+      }
+      // Track removed if normalized chi2 too high
+      if (chi2Norm > fMaxChi2) {
+       track->Remove();
+       break; // stop the search for this candidate:
+       // exit from the loop over station 
+      }
       if (fPrintLevel >= 10) {
        cout << "FollowTracks: track candidate(0..): " << trackIndex
             << " after fit from station(0..): " << station << " to 4" << endl;
index a1d8f96ce015ab81355cc0ddd56c927a15917bb0..5b9013fa62e4346ab6841682169d678563140222 100644 (file)
@@ -31,6 +31,10 @@ class AliMUONEventReconstructor : public TObject {
   // Get and Set, Set to defaults
   Double_t GetMinBendingMomentum(void) {return fMinBendingMomentum;}
   void SetMinBendingMomentum(Double_t MinBendingMomentum) {fMinBendingMomentum = MinBendingMomentum;}
+  Double_t GetMaxBendingMomentum(void) {return fMaxBendingMomentum;}
+  void SetMaxBendingMomentum(Double_t MaxBendingMomentum) {fMaxBendingMomentum = MaxBendingMomentum;}
+  Double_t GetMaxChi2(void) {return fMaxChi2;}
+  void SetMaxChi2(Double_t MaxChi2) {fMaxChi2 = MaxChi2;}
   Double_t GetMaxSigma2Distance(void) {return fMaxSigma2Distance;}
   void SetMaxSigma2Distance(Double_t MaxSigma2Distance) {fMaxSigma2Distance = MaxSigma2Distance;}
   Double_t GetBendingResolution(void) {return fBendingResolution;}
@@ -84,6 +88,9 @@ class AliMUONEventReconstructor : public TObject {
 
   // Parameters for event reconstruction
   Double_t fMinBendingMomentum; // minimum value (GeV/c) of momentum in bending plane
+  // Parameters for event reconstruction
+  Double_t fMaxBendingMomentum; // maximum value (GeV/c) of momentum in bending plane
+  Double_t fMaxChi2; // maximum Chi2 per degree of Freedom
   Double_t fMaxSigma2Distance; // maximum square distance in units of the variance (maximum chi2)
   Double_t fRMin[kMaxMuonTrackingChambers]; // minimum radius (cm)
   Double_t fRMax[kMaxMuonTrackingChambers]; // maximum radius (cm)
index 01671570e80cc2e02eb63c71086850806c7e275b..fcf369a5750deab291b85f31e2c0e9f46b9fa261 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:06:39  hristov
 Inline functions moved from *.cxx to *.h files instead of forward declarations
 
@@ -253,18 +258,27 @@ void AliMUONSegment::UpdateFromStationTrackParam(AliMUONTrackParam *TrackParam,
   // The "road" is parametrized from the old reco_muon.F
   // with 8 cm between stations.
   AliMUONTrackParam *param0;
-  Double_t cReso2, sReso2;
+//   Double_t cReso2, sReso2;
   // parameters to define the widths of the searching roads in station 0,1,2
   // width = p0 + p1/ (momentum)^2
   //                  station number:        0         1          2
-  static Double_t p0BendingCoor[3] =     { 6.43e-2, 1.64e-2,   0.034 };   
-  static Double_t p1BendingCoor[3] =     {    986.,    821.,    446. };  
-  static Double_t p0BendingSlope[3] =    { 3.54e-6, 3.63e-6,  3.6e-6 };  
-  static Double_t p1BendingSlope[3] =    { 4.49e-3,  4.8e-3,   0.011 };  
-  static Double_t p0NonBendingCoor[3] =  { 4.66e-2, 4.83e-2,   0.049 };   
-  static Double_t p1NonBendingCoor[3] =  {   1444.,    866.,    354. };  
-  static Double_t p0NonBendingSlope[3] = { 6.14e-4, 6.49e-4, 6.85e-4 };  
-  static Double_t p1NonBendingSlope[3] = {      0.,      0.,      0. };  
+//   static Double_t p0BendingCoor[3] =     { 6.43e-2, 1.64e-2,   0.034 };   
+//   static Double_t p1BendingCoor[3] =     {    986.,    821.,    446. };  
+//   static Double_t p0BendingSlope[3] =    { 3.54e-6, 3.63e-6,  3.6e-6 };  
+//   static Double_t p1BendingSlope[3] =    { 4.49e-3,  4.8e-3,   0.011 };  
+//   static Double_t p0NonBendingCoor[3] =  { 4.66e-2, 4.83e-2,   0.049 };   
+//   static Double_t p1NonBendingCoor[3] =  {   1444.,    866.,    354. };  
+//   static Double_t p0NonBendingSlope[3] = { 6.14e-4, 6.49e-4, 6.85e-4 };  
+//   static Double_t p1NonBendingSlope[3] = {      0.,      0.,      0. };
+  
+  static Double_t p0BendingCoor[3] =     { 6.43e-2, 6.43e-2,   6.43e-2  };   
+  static Double_t p1BendingCoor[3] =     {    986.,    986.,       986. };  
+  static Double_t p0BendingSlope[3] =    {   3.6e-6,   3.6e-6,     3.6e-6  };  
+  static Double_t p1BendingSlope[3] =    {  1.1e-2,  1.1e-2,    1.1e-2  };  
+  static Double_t p0NonBendingCoor[3] =  {   0.049,   0.049,     0.049  };   
+  static Double_t p1NonBendingCoor[3] =  {   1444.,   1444.,      1444. };  
+  static Double_t p0NonBendingSlope[3] = {   6.8e-4,   6.8e-4,     6.8e-4  };  
+  static Double_t p1NonBendingSlope[3] = {      0.,      0.,         0. };  
   param0 = &(TrackParam[0]);
 
 // OLD version
@@ -296,15 +310,15 @@ void AliMUONSegment::UpdateFromStationTrackParam(AliMUONTrackParam *TrackParam,
   // because they are added in the functions "NormalizedChi2WithSegment"
   // and "NormalizedChi2WithHitForRec"
   // Bending plane
-  cReso2 = fBendingCoorReso2;
-  sReso2 = (2. * cReso2 )/ (Dz3*Dz3) ;
-  fBendingCoorReso2 = p0BendingCoor[Station] + p1BendingCoor[Station]*InverseMomentum*InverseMomentum - cReso2;
-  fBendingSlopeReso2 = p0BendingSlope[Station] + p1BendingSlope[Station]*InverseMomentum*InverseMomentum - sReso2;
+//   cReso2 = fBendingCoorReso2;
+//   sReso2 = (2. * cReso2 )/ (Dz3*Dz3) ;
+  fBendingCoorReso2 = p0BendingCoor[Station] + p1BendingCoor[Station]*InverseMomentum*InverseMomentum ;  // - cReso2
+  fBendingSlopeReso2 = p0BendingSlope[Station] + p1BendingSlope[Station]*InverseMomentum*InverseMomentum; //  - sReso2;
   // Non bending plane
-  cReso2 = fNonBendingCoorReso2;
-  sReso2 =  (2. * cReso2 )/ (Dz3*Dz3) ;
-  fNonBendingCoorReso2 = p0NonBendingCoor[Station] + p1NonBendingCoor[Station]*InverseMomentum*InverseMomentum - cReso2;
-  fNonBendingSlopeReso2 = p0NonBendingSlope[Station] + p1NonBendingSlope[Station]*InverseMomentum*InverseMomentum - sReso2;
+//   cReso2 = fNonBendingCoorReso2;
+//   sReso2 =  (2. * cReso2 )/ (Dz3*Dz3) ;
+  fNonBendingCoorReso2 = p0NonBendingCoor[Station] + p1NonBendingCoor[Station]*InverseMomentum*InverseMomentum; // - cReso2;
+  fNonBendingSlopeReso2 = p0NonBendingSlope[Station] + p1NonBendingSlope[Station]*InverseMomentum*InverseMomentum; //  - sReso2;
   return;
 }
 
index fd5f2e6a77a21cc6e8bca60ae0027c76914dbf0d..e7a6fbf537e7376e22eaf6eebca6c8076eea45cc 100644 (file)
@@ -29,6 +29,30 @@ class AliMUONSegment : public TObject {
   AliMUONHitForRec* GetHitForRec2(void) {
     // Get fHitForRecPtr2
     return fHitForRecPtr2;}
+  Double_t GetBendingCoor(void) {
+    // Get fBendingCoor
+    return fBendingCoor;}
+  void SetBendingCoor(Double_t BendingCoor) {
+    // Set fBendingCoor
+    fBendingCoor = BendingCoor;}
+  Double_t GetBendingSlope(void) {
+    // Get fBendingSlope
+    return fBendingSlope;}
+  void SetBendingSlope(Double_t BendingSlope) {
+    // Set fBendingSlope
+    fBendingSlope = BendingSlope;}
+  Double_t GetNonBendingCoor(void) {
+    // Get fNonBendingCoor
+    return fNonBendingCoor;}
+  void SetNonBendingCoor(Double_t NonBendingCoor) {
+    // Set fNonBendingCoor
+    fNonBendingCoor = NonBendingCoor;}
+  Double_t GetNonBendingSlope(void) {
+    // Get fNonBendingSlope
+    return fNonBendingSlope;}
+  void SetNonBendingSlope(Double_t NonBendingSlope) {
+    // Set fNonBendingSlope
+    fNonBendingSlope = NonBendingSlope;}
   Double_t GetBendingCoorReso2(void) {
     // Get fBendingCoorReso2
     return fBendingCoorReso2;}
index 83fc6e228c6f015520e2a483bfb04a9377d5b08d..6750b083ad7f974742ed352a7b3e983f64560cf0 100644 (file)
 
 /*
 $Log$
+Revision 1.7  2000/09/19 15:50:46  gosset
+TrackChi2MCS function: covariance matrix better calculated,
+taking into account missing planes...
+
 Revision 1.6  2000/07/20 12:45:27  gosset
 New "EventReconstructor..." structure,
        hopefully more adapted to tree/streamer.
@@ -66,7 +70,7 @@ Addition of files for track reconstruction in C++
 
 #include <TClonesArray.h>
 #include <TMath.h>
-#include <TMatrix.h>
+#include <TMatrixD.h>
 #include <TObjArray.h>
 #include <TVirtualFitter.h>
 
@@ -301,7 +305,7 @@ void AliMUONTrack::Fit()
   // choice of function to be minimized according to fFitMCS
   if (fFitMCS == 0) fgFitter->SetFCN(TrackChi2);
   else fgFitter->SetFCN(TrackChi2MCS);
-  arg[0] = 1;
+  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 ????)
@@ -319,13 +323,15 @@ void AliMUONTrack::Fit()
                         trackParam->GetNonBendingSlope(),
                         0.001, -0.5, 0.5);
   if (fFitNParam == 5) {
-    // set last 2 Minuit parameters (no limits for the search: min=max=0)
+    // set last 2 Minuit parameters
+    // mandatory limits in Bending to avoid NaN values of parameters
     fgFitter->SetParameter(3, "X",
                           trackParam->GetNonBendingCoor(),
-                          0.03, 0.0, 0.0);
+                          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, 0.0, 0.0);
+                          0.10, -500.0, 500.0);
   }
   // search without gradient calculation in the function
   fgFitter->ExecuteCommand("SET NOGRADIENT", arg, 0);
@@ -517,8 +523,8 @@ void TrackChi2MCS(Int_t &NParam, Double_t *Gradient, Double_t &Chi2, Double_t *P
   Double_t hbc1, hbc2, pbc1, pbc2;
   Double_t hnbc1, hnbc2, pnbc1, pnbc2;
   Int_t numberOfHit = trackBeingFitted->GetNTrackHits();
-  TMatrix *covBending = new TMatrix(numberOfHit, numberOfHit);
-  TMatrix *covNonBending = new TMatrix(numberOfHit, numberOfHit);
+  TMatrixD *covBending = new TMatrixD(numberOfHit, numberOfHit);
+  TMatrixD *covNonBending = new TMatrixD(numberOfHit, numberOfHit);
   Double_t *msa2 = new Double_t[numberOfHit];
 
   // Predicted coordinates and  multiple scattering angles are first calculated
@@ -583,7 +589,23 @@ void TrackChi2MCS(Int_t &NParam, Double_t *Gradient, Double_t &Chi2, Double_t *P
       }
     } // for (hitNumber2 = hitNumber1;...
   } // for (hitNumber1 = 0;...
-
+  // Normalization of covariance matrices
+  Double_t normCovBending2 = covBending->E2Norm();
+  Double_t normCovNonBending2 = covNonBending->E2Norm();
+  (*covBending) *= 1/normCovBending2;
+  (*covNonBending) *= 1/normCovNonBending2;
+//   if (covBending->Determinant() < 1.e-33) {
+//     printf(" *** covBending *** \n");
+//     covBending->Print();
+//     printf(" *** covNonBending *** \n");
+//     covNonBending->Print();
+//     cout << " number of hits " <<  numberOfHit << endl;
+//     cout << "Momentum = " << 1/Param[0] <<endl;
+//     cout << "normCovBending = " << normCovBending2 << endl; 
+//     cout << "normCovNonBending = " << normCovNonBending2 << endl; 
+//     exit(0);
+    
+//   }
   // Inverts covariance matrix 
   goodDeterminant = kTRUE;
   // check whether the Invert method returns flag if matrix cannot be inverted,
@@ -608,6 +630,9 @@ void TrackChi2MCS(Int_t &NParam, Double_t *Gradient, Double_t &Chi2, Double_t *P
   // Calculates Chi2
   if (goodDeterminant) {
     // with Multiple Scattering if inversion correct
+    // Inverse  matrices without normalization
+    (*covBending) *= 1/normCovBending2;
+    (*covNonBending) *= 1/normCovNonBending2;
     for (hitNumber1 = 0; hitNumber1 < numberOfHit ; hitNumber1++) { 
       hit1 = (AliMUONTrackHit*) (*(trackBeingFitted->GetTrackHitsPtr()))[hitNumber1];
       hbc1 = hit1->GetHitForRecPtr()->GetBendingCoor();