In reconstruction classes:
authorivana <ivana@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 22 Apr 2009 10:18:35 +0000 (10:18 +0000)
committerivana <ivana@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 22 Apr 2009 10:18:35 +0000 (10:18 +0000)
- Improve track candidate selection by taking into account the resolution
  on the track parameters.
- Add a flag into RecoParam to switch between track selection on their
  slope or on their impact parameter at vertex (default).
- UPdate REcoParam in OCDB
- Take into account the chamber/cluster resolution when defining
  the "fixed" area to quickly look for new cluster to be attached to the track.
- AliMUONDigitCalibrator: Do not use digit saturation information from
  calibration file when using the calibration option NOGAIN.
(Philippe P.)

MUON/AliMUONDigitCalibrator.cxx
MUON/AliMUONRecoParam.cxx
MUON/AliMUONRecoParam.h
MUON/AliMUONTrackExtrap.cxx
MUON/AliMUONTrackParam.cxx
MUON/AliMUONTrackReconstructor.cxx
MUON/AliMUONTrackReconstructorK.cxx
MUON/AliMUONVTrackReconstructor.cxx
MUON/AliMUONVTrackReconstructor.h
MUON/runDataReconstruction.C
OCDB/MUON/Calib/RecoParam/Run0_999999999_v0_s0.root

index aecf707..a307f6c 100644 (file)
@@ -357,7 +357,7 @@ AliMUONDigitCalibrator::CalibrateDigit(Int_t detElemId, Int_t manuId, Int_t manu
   {
     Int_t saturation(3000);
   
-    if ( gain )
+    if ( gain && ( fApplyGains != fgkNoGain ) )
     {
       saturation = gain->ValueAsInt(manuChannel,4);
     }
index 5291eaf..d7840a7 100644 (file)
@@ -71,7 +71,8 @@ AliMUONRecoParam::AliMUONRecoParam()
   fChargeSigmaCut(4.0),
   fRemoveConnectedTracksInSt12(kFALSE),
   fMaxTriggerTracks(0),
-  fMaxTrackCandidates(0)
+  fMaxTrackCandidates(0),
+  fSelectTrackOnSlope(kFALSE)
 {
   /// Constructor
   
@@ -154,8 +155,9 @@ void AliMUONRecoParam::SetLowFluxParam()
   fMaxBendingMomentum = 3000.;
   fMaxNonBendingSlope = 0.3;
   fMaxBendingSlope = 0.4;
-  fNonBendingVertexDispersion = 10.;
-  fBendingVertexDispersion = 10.;
+  fSelectTrackOnSlope = kFALSE;
+  fNonBendingVertexDispersion = 70.;
+  fBendingVertexDispersion = 70.;
   fMaxNonBendingDistanceToTrack = 1.;
   fMaxBendingDistanceToTrack = 1.;
   fSigmaCutForTracking = 6.;
@@ -196,8 +198,9 @@ void AliMUONRecoParam::SetHighFluxParam()
   fMaxBendingMomentum = 3000.;
   fMaxNonBendingSlope = 0.3;
   fMaxBendingSlope = 0.4;
-  fNonBendingVertexDispersion = 10.;
-  fBendingVertexDispersion = 10.;
+  fSelectTrackOnSlope = kFALSE;
+  fNonBendingVertexDispersion = 70.;
+  fBendingVertexDispersion = 70.;
   fMaxNonBendingDistanceToTrack = 1.;
   fMaxBendingDistanceToTrack = 1.;
   fSigmaCutForTracking = 6.;
@@ -236,12 +239,13 @@ void AliMUONRecoParam::SetCosmicParam()
   SetEventSpecie(AliRecoParam::kCosmic);
   fMinBendingMomentum = 1.;
   fMaxBendingMomentum = 10000000.;
-  fMaxNonBendingSlope = 0.3;
-  fMaxBendingSlope = 0.4;
-  fNonBendingVertexDispersion = 10.;
-  fBendingVertexDispersion = 10.;
-  fMaxNonBendingDistanceToTrack = 10.;
-  fMaxBendingDistanceToTrack = 10.;
+  fMaxNonBendingSlope = 0.4;
+  fMaxBendingSlope = 0.5;
+  fSelectTrackOnSlope = kTRUE;
+  fNonBendingVertexDispersion = 200.;
+  fBendingVertexDispersion = 200.;
+  fMaxNonBendingDistanceToTrack = 1.;
+  fMaxBendingDistanceToTrack = 1.;
   fSigmaCutForTracking = 7.;
   fSigmaCutForImprovement = 7.;
   fSigmaCutForTrigger = 8.;
@@ -325,16 +329,14 @@ void AliMUONRecoParam::Print(Option_t *option) const
   if (fSaveFullClusterInESD) cout<<Form("Save all cluster info in ESD for %5.2f %% of events",fPercentOfFullClusterInESD)<<endl;
   else cout<<"Save partial cluster info in ESD"<<endl;
     
-  cout<<Form("Bending momentum range = [%5.2f,%5.2f]",fMinBendingMomentum,fMaxBendingMomentum)<<endl;
-  
-  cout<<Form("Maximum non bending slope = %5.2f",fMaxNonBendingSlope)<<endl;
-  
-  cout<<Form("Maximum bending slope (used only if B=0) = %5.2f",fMaxBendingSlope)<<endl;
-  
-  if (strstr(fTrackingMode,"ORIGINAL"))
-    cout<<Form("Vertex dispertion = (%5.2f,%5.2f)",fNonBendingVertexDispersion,fBendingVertexDispersion)<<endl;
-  else if (strstr(option,"FULL"))
-    cout<<Form("Vertex dispertion (used for original tracking only) = (%5.2f,%5.2f)",fNonBendingVertexDispersion,fBendingVertexDispersion)<<endl;
+  cout<<"Selection of track candidates:"<<endl;
+  if (fSelectTrackOnSlope) cout<<Form("\t- Non-bending slope < %5.2f",fMaxNonBendingSlope)<<endl;
+  else cout<<"\t- Impact parameter < 3 * vertex dispersion in the non-bending direction"<<endl;
+  cout<<Form("\t- if B!=0: Bending momentum > %5.2f",fMinBendingMomentum)<<endl;
+  if (fSelectTrackOnSlope) cout<<Form("\t  if B==0: Bending slope < %5.2f",fMaxBendingSlope)<<endl;
+  else cout<<"\t  if B==0: Impact parameter < 3 * vertex dispersion in the bending direction"<<endl;
+  
+  cout<<Form("Vertex dispersion (used to estimate initial bending momentum resolution) = (%5.2f,%5.2f)",fNonBendingVertexDispersion,fBendingVertexDispersion)<<endl;
   
   cout<<Form("Maximum distance to track = (%5.2f,%5.2f)",fMaxNonBendingDistanceToTrack,fMaxBendingDistanceToTrack)<<endl;
   
index 9fa4c20..798143c 100644 (file)
@@ -61,6 +61,7 @@ class AliMUONRecoParam : public AliDetectorRecoParam
   void     SetMaxBendingMomentum(Double_t val) {fMaxBendingMomentum = val;}
   /// return the maximum value (GeV/c) of momentum in bending plane
   Double_t GetMaxBendingMomentum() const {return fMaxBendingMomentum;}
+  
   /// set the maximum value of the non bending slope
   void     SetMaxNonBendingSlope(Double_t val) {fMaxNonBendingSlope = val;}
   /// return the maximum value of the non bending slope
@@ -70,6 +71,11 @@ class AliMUONRecoParam : public AliDetectorRecoParam
   /// return the maximum value of the bending slope
   Double_t GetMaxBendingSlope() const {return fMaxBendingSlope;}
   
+  /// switch on/off the track selection according to their slope (instead of their impact parameter)
+  void     SelectOnTrackSlope(Bool_t flag) {fSelectTrackOnSlope = flag;}
+  /// return kTRUE/kFALSE if tracks are selected according to their slope/impact parameter
+  Bool_t   SelectOnTrackSlope() const {return fSelectTrackOnSlope;}
+  
   /// set the vertex dispersion (cm) in non bending plane (used for original tracking only)
   void     SetNonBendingVertexDispersion(Double_t val) {fNonBendingVertexDispersion = val;} 
   /// return the vertex dispersion (cm) in non bending plane (used for original tracking only)
@@ -347,13 +353,15 @@ class AliMUONRecoParam : public AliDetectorRecoParam
   Int_t      fMaxTriggerTracks; ///< maximum number of trigger tracks above which the tracking is cancelled
   Int_t      fMaxTrackCandidates; ///< maximum number of track candidates above which the tracking abort
   
+  Bool_t     fSelectTrackOnSlope; ///< select track candidates according to their slope (instead of their impact parameter)
+  
   // functions
   void SetLowFluxParam();
   void SetHighFluxParam();
   void SetCosmicParam();
   
   
-  ClassDef(AliMUONRecoParam,12) // MUON reco parameters
+  ClassDef(AliMUONRecoParam,13) // MUON reco parameters
 };
 
 #endif
index ba3fac9..97b3f82 100644 (file)
@@ -71,7 +71,7 @@ Double_t AliMUONTrackExtrap::GetImpactParamFromBendingMomentum(Double_t bendingM
   
   if (bendingMomentum == 0.) return 1.e10;
   
-  const Double_t kCorrectionFactor = 0.9; // impact parameter is 10% overestimated
+  const Double_t kCorrectionFactor = 1.1; // impact parameter is 10% underestimated
   
   return kCorrectionFactor * (-0.0003 * fgSimpleBValue * fgkSimpleBLength * fgkSimpleBPosition / bendingMomentum);
 }
@@ -315,6 +315,7 @@ void AliMUONTrackExtrap::ExtrapToZCov(AliMUONTrackParam* trackParam, Double_t zE
   TMatrixD jacob(5,5);
   jacob.Zero();
   TMatrixD dParam(5,1);
+  Double_t direction[5] = {-1.,-1.,1.,1.,-1.};
   for (Int_t i=0; i<5; i++) {
     // Skip jacobian calculation for parameters with no associated error
     if (kParamCov(i,i) <= 0.) continue;
@@ -323,7 +324,7 @@ void AliMUONTrackExtrap::ExtrapToZCov(AliMUONTrackParam* trackParam, Double_t zE
     for (Int_t j=0; j<5; j++) {
       if (j==i) {
         dParam(j,0) = TMath::Sqrt(kParamCov(i,i));
-       dParam(j,0) *= TMath::Sign(1.,paramSave(j,0)); // variation always in the same direction
+       dParam(j,0) *= TMath::Sign(1.,direction[j]*paramSave(j,0)); // variation always in the same direction
       } else dParam(j,0) = 0.;
     }
     
index 045920c..7e5796d 100644 (file)
@@ -52,6 +52,7 @@ AliMUONTrackParam::AliMUONTrackParam()
     fLocalChi2(0.)
 {
   /// Constructor
+  fParameters.Zero();
 }
 
   //_________________________________________________________________________
index 3f0799c..1c3936f 100644 (file)
@@ -118,13 +118,21 @@ Bool_t AliMUONTrackReconstructor::MakeTrackCandidates(AliMUONVClusterStore& clus
       
       // Remove track if no cluster found on a requested station
       // or abort tracking if there are too many candidates
-      if (!clusterFound && GetRecoParam()->RequestStation(7-istat)) {
-        fRecTracksPtr->Remove(track);
-       fNRecTracks--;
-      }        else if (fNRecTracks > GetRecoParam()->GetMaxTrackCandidates()) {
-       AliError(Form("Too many track candidates (%d tracks). Abort tracking.", fNRecTracks));
-       delete segments;
-       return kFALSE;
+      if (GetRecoParam()->RequestStation(7-istat)) {
+       if (!clusterFound) {
+         fRecTracksPtr->Remove(track);
+         fNRecTracks--;
+       } else if (fNRecTracks > GetRecoParam()->GetMaxTrackCandidates()) {
+         AliError(Form("Too many track candidates (%d tracks). Abort tracking.", fNRecTracks));
+         delete segments;
+         return kFALSE;
+       }
+      } else {
+       if ((fNRecTracks + segments->GetEntriesFast() - iseg - 1) > GetRecoParam()->GetMaxTrackCandidates()) {
+         AliError(Form("Too many track candidates (%d tracks). Abort tracking.", fNRecTracks + segments->GetEntriesFast() - iseg - 1));
+         delete segments;
+         return kFALSE;
+       }
       }
       
     }
@@ -199,8 +207,8 @@ Bool_t AliMUONTrackReconstructor::MakeMoreTrackCandidates(AliMUONVClusterStore&
        }
        
        // abort tracking if there are too many candidates
-       if (fNRecTracks > GetRecoParam()->GetMaxTrackCandidates()) {
-         AliError(Form("Too many track candidates (%d tracks). Abort tracking.", fNRecTracks));
+       if ((fNRecTracks + segments->GetEntriesFast() - iSegment - 1) > GetRecoParam()->GetMaxTrackCandidates()) {
+         AliError(Form("Too many track candidates (%d tracks). Abort tracking.", fNRecTracks + segments->GetEntriesFast() - iSegment - 1));
          delete segments;
          return kFALSE;
        }
@@ -255,15 +263,11 @@ Bool_t AliMUONTrackReconstructor::FollowTracks(AliMUONVClusterStore& clusterStor
        Fit(*track, kFALSE, kTRUE, kTRUE);
       else Fit(*track, kFALSE, kFALSE, kTRUE);
       
-      // remove track with absolute bending momentum out of limits
-      if (AliMUONTrackExtrap::IsFieldON()) {
-       Double_t bendingMomentum = TMath::Abs(1. / ((AliMUONTrackParam*)track->GetTrackParamAtCluster()->First())->GetInverseBendingMomentum());
-       if (bendingMomentum < GetRecoParam()->GetMinBendingMomentum() ||
-           bendingMomentum > GetRecoParam()->GetMaxBendingMomentum()) {
-         fRecTracksPtr->Remove(track);
-         fNRecTracks--;
-         continue;
-       }
+      // remove tracks out of limits
+      if (!IsAcceptable(*((AliMUONTrackParam*)track->GetTrackParamAtCluster()->First()))) {
+       fRecTracksPtr->Remove(track);
+       fNRecTracks--;
+       continue;
       }
       
       // remove track if the normalized chi2 is too high
@@ -945,6 +949,7 @@ Double_t AliMUONTrackReconstructor::TryTwoClusters(const AliMUONTrackParam &trac
     jacob(1,2) = 1.; // dy1/dy
     // first derivative at the second cluster:
     TMatrixD dParam(5,1);
+    Double_t direction[5] = {-1.,-1.,1.,1.,-1.};
     for (Int_t i=0; i<5; i++) {
       // Skip jacobian calculation for parameters with no associated error
       if (kParamCov(i,i) == 0.) continue;
@@ -952,7 +957,7 @@ Double_t AliMUONTrackReconstructor::TryTwoClusters(const AliMUONTrackParam &trac
       for (Int_t j=0; j<5; j++) {
         if (j==i) {
           dParam(j,0) = TMath::Sqrt(kParamCov(i,i));
-         if (j == 4) dParam(j,0) *= TMath::Sign(1.,-paramAtCluster1Save(4,0)); // variation always in the same direction
+         dParam(j,0) *= TMath::Sign(1.,direction[j]*paramAtCluster1Save(j,0)); // variation always in the same direction
         } else dParam(j,0) = 0.;
       }
       
@@ -1107,12 +1112,8 @@ Bool_t AliMUONTrackReconstructor::RecoverTrack(AliMUONTrack &trackCandidate, Ali
   // Calculate the track parameter covariance matrix
   Fit(trackCandidate, kFALSE, kFALSE, kTRUE);
   
-  // skip track with absolute bending momentum out of limits
-  if (AliMUONTrackExtrap::IsFieldON()) {
-    Double_t bendingMomentum = TMath::Abs(1. / ((AliMUONTrackParam*)trackCandidate.GetTrackParamAtCluster()->First())->GetInverseBendingMomentum());
-    if (bendingMomentum < GetRecoParam()->GetMinBendingMomentum() ||
-       bendingMomentum > GetRecoParam()->GetMaxBendingMomentum()) return kFALSE;
-  }
+  // skip track out of limits
+  if (!IsAcceptable(*((AliMUONTrackParam*)trackCandidate.GetTrackParamAtCluster()->First()))) return kFALSE;
   
   // Look for new cluster(s) in next station
   return FollowTrackInStation(trackCandidate,clusterStore,nextStation);
index af722e5..c5e98cb 100644 (file)
@@ -116,13 +116,21 @@ Bool_t AliMUONTrackReconstructorK::MakeTrackCandidates(AliMUONVClusterStore& clu
       
       // Remove track if no cluster found on a requested station
       // or abort tracking if there are too many candidates
-      if (!clusterFound && GetRecoParam()->RequestStation(7-istat)) {
-       fRecTracksPtr->Remove(track);
-       fNRecTracks--;
-      } else if (fNRecTracks > GetRecoParam()->GetMaxTrackCandidates()) {
-       AliError(Form("Too many track candidates (%d tracks). Abort tracking.", fNRecTracks));
-       delete segments;
-       return kFALSE;
+      if (GetRecoParam()->RequestStation(7-istat)) {
+       if (!clusterFound) {
+         fRecTracksPtr->Remove(track);
+         fNRecTracks--;
+       } else if (fNRecTracks > GetRecoParam()->GetMaxTrackCandidates()) {
+         AliError(Form("Too many track candidates (%d tracks). Abort tracking.", fNRecTracks));
+         delete segments;
+         return kFALSE;
+       }
+      } else {
+       if ((fNRecTracks + segments->GetEntriesFast() - iSegment - 1) > GetRecoParam()->GetMaxTrackCandidates()) {
+         AliError(Form("Too many track candidates (%d tracks). Abort tracking.", fNRecTracks + segments->GetEntriesFast() - iSegment - 1));
+         delete segments;
+         return kFALSE;
+       }
       }
       
     }
@@ -145,21 +153,17 @@ Bool_t AliMUONTrackReconstructorK::MakeTrackCandidates(AliMUONVClusterStore& clu
     // retrace tracks using Kalman filter
     RetraceTrack(*track,kTRUE);
     
-    // remove tracks with absolute bending momentum out of limits
-    if (GetRecoParam()->MakeTrackCandidatesFast() && AliMUONTrackExtrap::IsFieldON()) {
-      AliMUONTrackParam *trackParam = (AliMUONTrackParam*)track->GetTrackParamAtCluster()->First();
-      Double_t bendingMomentum = TMath::Abs(1. / trackParam->GetInverseBendingMomentum());
-      if (bendingMomentum < GetRecoParam()->GetMinBendingMomentum() || bendingMomentum > GetRecoParam()->GetMaxBendingMomentum()) {
-       fRecTracksPtr->Remove(track);
-       fNRecTracks--;
-      }
+    // remove tracks out of limits
+    if (!IsAcceptable(*((AliMUONTrackParam*)track->GetTrackParamAtCluster()->First()))) {
+      fRecTracksPtr->Remove(track);
+      fNRecTracks--;
     }
     
   }
   
   // Keep only the best tracks if required
   if (!GetRecoParam()->TrackAllTracks()) RemoveDoubleTracks();
-  else if (GetRecoParam()->MakeTrackCandidatesFast() && AliMUONTrackExtrap::IsFieldON()) fRecTracksPtr->Compress();
+  else if (AliMUONTrackExtrap::IsFieldON()) fRecTracksPtr->Compress();
   
   AliDebug(1,Form("Number of good candidates = %d",fNRecTracks));
   
@@ -223,8 +227,8 @@ Bool_t AliMUONTrackReconstructorK::MakeMoreTrackCandidates(AliMUONVClusterStore&
        }
        
        // abort tracking if there are too many candidates
-       if (fNRecTracks > GetRecoParam()->GetMaxTrackCandidates()) {
-         AliError(Form("Too many track candidates (%d tracks). Abort tracking.", fNRecTracks));
+       if ((fNRecTracks + segments->GetEntriesFast() - iSegment - 1) > GetRecoParam()->GetMaxTrackCandidates()) {
+         AliError(Form("Too many track candidates (%d tracks). Abort tracking.", fNRecTracks + segments->GetEntriesFast() - iSegment - 1));
          delete segments;
          return kFALSE;
        }
@@ -248,13 +252,9 @@ Bool_t AliMUONTrackReconstructorK::MakeMoreTrackCandidates(AliMUONVClusterStore&
     RetraceTrack(*track,kTRUE);
     
     // remove tracks with absolute bending momentum out of limits
-    if (AliMUONTrackExtrap::IsFieldON()) {
-      AliMUONTrackParam *trackParam = (AliMUONTrackParam*)track->GetTrackParamAtCluster()->First();
-      Double_t bendingMomentum = TMath::Abs(1. / trackParam->GetInverseBendingMomentum());
-      if (bendingMomentum < GetRecoParam()->GetMinBendingMomentum() || bendingMomentum > GetRecoParam()->GetMaxBendingMomentum()) {
-       fRecTracksPtr->Remove(track);
-       fNRecTracks--;
-      }
+    if (!IsAcceptable(*((AliMUONTrackParam*)track->GetTrackParamAtCluster()->First()))) {
+      fRecTracksPtr->Remove(track);
+      fNRecTracks--;
     }
     
   }
@@ -471,9 +471,16 @@ Bool_t AliMUONTrackReconstructorK::FollowTracks(AliMUONVClusterStore& clusterSto
       }
       
       // abort tracking if there are too many candidates
-      if (fNRecTracks > GetRecoParam()->GetMaxTrackCandidates()) {
-       AliError(Form("Too many track candidates (%d tracks). Abort tracking.", fNRecTracks));
-       return kFALSE;
+      if (GetRecoParam()->RequestStation(station)) {
+       if (fNRecTracks > GetRecoParam()->GetMaxTrackCandidates()) {
+         AliError(Form("Too many track candidates (%d tracks). Abort tracking.", fNRecTracks));
+         return kFALSE;
+       }
+      } else {
+       if ((fNRecTracks + currentNRecTracks - iRecTrack - 1) > GetRecoParam()->GetMaxTrackCandidates()) {
+         AliError(Form("Too many track candidates (%d tracks). Abort tracking.", fNRecTracks + currentNRecTracks - iRecTrack - 1));
+         return kFALSE;
+       }
       }
       
     }
@@ -501,7 +508,6 @@ Bool_t AliMUONTrackReconstructorK::FollowTrackInChamber(AliMUONTrack &trackCandi
   /// return kTRUE if new cluster(s) have been found (otherwise return kFALSE)
   AliDebug(1,Form("Enter FollowTrackInChamber(1..) %d", nextChamber+1));
   
-  Double_t bendingMomentum;
   Double_t chi2OfCluster;
   Double_t maxChi2OfCluster = 2. * GetRecoParam()->GetSigmaCutForTracking() *
                                   GetRecoParam()->GetSigmaCutForTracking(); // 2 because 2 quantities in chi2
@@ -591,12 +597,8 @@ Bool_t AliMUONTrackReconstructorK::FollowTrackInChamber(AliMUONTrack &trackCandi
       // Compute new track parameters including new cluster using kalman filter
       addChi2TrackAtCluster = RunKalmanFilter(extrapTrackParamAtCluster);
       
-      // skip track with absolute bending momentum out of limits
-      if (AliMUONTrackExtrap::IsFieldON()) {
-        bendingMomentum = TMath::Abs(1. / extrapTrackParamAtCluster.GetInverseBendingMomentum());
-        if (bendingMomentum < GetRecoParam()->GetMinBendingMomentum() ||
-           bendingMomentum > GetRecoParam()->GetMaxBendingMomentum()) continue;
-      }
+      // skip track out of limits
+      if (!IsAcceptable(extrapTrackParamAtCluster)) continue;
       
       // remember a cluster was found
       foundOneCluster = kTRUE;
@@ -679,7 +681,6 @@ Bool_t AliMUONTrackReconstructorK::FollowTrackInStation(AliMUONTrack &trackCandi
     ch2 = 2*nextStation+1;
   }
   
-  Double_t bendingMomentum;
   Double_t chi2OfCluster;
   Double_t maxChi2OfCluster = 2. * GetRecoParam()->GetSigmaCutForTracking() *
                                   GetRecoParam()->GetSigmaCutForTracking(); // 2 because 2 quantities in chi2
@@ -782,12 +783,8 @@ Bool_t AliMUONTrackReconstructorK::FollowTrackInStation(AliMUONTrack &trackCandi
       // Compute new track parameters including "clusterCh2" using kalman filter
       addChi2TrackAtCluster2 = RunKalmanFilter(extrapTrackParamAtCluster2);
       
-      // skip track with absolute bending momentum out of limits
-      if (AliMUONTrackExtrap::IsFieldON()) {
-       bendingMomentum = TMath::Abs(1. / extrapTrackParamAtCluster2.GetInverseBendingMomentum());
-       if (bendingMomentum < GetRecoParam()->GetMinBendingMomentum() ||
-           bendingMomentum > GetRecoParam()->GetMaxBendingMomentum()) continue;
-      }
+      // skip track out of limits
+      if (!IsAcceptable(extrapTrackParamAtCluster2)) continue;
       
       // Printout for debuging
       if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructorK") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) {
@@ -840,12 +837,8 @@ Bool_t AliMUONTrackReconstructorK::FollowTrackInStation(AliMUONTrack &trackCandi
           // Compute new track parameters including "clusterCh1" using kalman filter
           addChi2TrackAtCluster1 = RunKalmanFilter(extrapTrackParamAtCluster1);
           
-         // skip track with absolute bending momentum out of limits
-         if (AliMUONTrackExtrap::IsFieldON()) {
-           bendingMomentum = TMath::Abs(1. / extrapTrackParamAtCluster1.GetInverseBendingMomentum());
-           if (bendingMomentum < GetRecoParam()->GetMinBendingMomentum() ||
-               bendingMomentum > GetRecoParam()->GetMaxBendingMomentum()) continue;
-         }
+         // skip track out of limits
+         if (!IsAcceptable(extrapTrackParamAtCluster1)) continue;
          
          // remember a second cluster was found
          foundSecondCluster = kTRUE;
@@ -970,12 +963,8 @@ Bool_t AliMUONTrackReconstructorK::FollowTrackInStation(AliMUONTrack &trackCandi
         // Compute new track parameters including "clusterCh1" using kalman filter
         addChi2TrackAtCluster1 = RunKalmanFilter(extrapTrackParamAtCluster1);
         
-       // skip track with absolute bending momentum out of limits
-       if (AliMUONTrackExtrap::IsFieldON()) {
-         bendingMomentum = TMath::Abs(1. / extrapTrackParamAtCluster1.GetInverseBendingMomentum());
-         if (bendingMomentum < GetRecoParam()->GetMinBendingMomentum() ||
-             bendingMomentum > GetRecoParam()->GetMaxBendingMomentum()) continue;
-       }
+       // skip track out of limits
+       if (!IsAcceptable(extrapTrackParamAtCluster1)) continue;
        
        // remember a cluster was found
        foundOneCluster = kTRUE;
@@ -1232,12 +1221,8 @@ Bool_t AliMUONTrackReconstructorK::RecoverTrack(AliMUONTrack &trackCandidate, Al
   // Re-calculate track parameters at the (new) first cluster
   RetracePartialTrack(trackCandidate,(AliMUONTrackParam*)trackCandidate.GetTrackParamAtCluster()->UncheckedAt(1));
   
-  // skip track with absolute bending momentum out of limits
-  if (AliMUONTrackExtrap::IsFieldON()) {
-    Double_t bendingMomentum = TMath::Abs(1. / ((AliMUONTrackParam*)trackCandidate.GetTrackParamAtCluster()->First())->GetInverseBendingMomentum());
-    if (bendingMomentum < GetRecoParam()->GetMinBendingMomentum() ||
-       bendingMomentum > GetRecoParam()->GetMaxBendingMomentum()) return kFALSE;
-  }
+  // skip track out of limits
+  if (!IsAcceptable(*((AliMUONTrackParam*)trackCandidate.GetTrackParamAtCluster()->First()))) return kFALSE;
   
   // Look for new cluster(s) in next station
   return FollowTrackInStation(trackCandidate, clusterStore, nextStation);
index 6c0d15d..797c704 100644 (file)
@@ -100,7 +100,8 @@ AliMUONVTrackReconstructor::AliMUONVTrackReconstructor(const AliMUONRecoParam* r
 fRecTracksPtr(0x0),
 fNRecTracks(0),
 fClusterServer(clusterServer),
-fkRecoParam(recoParam)
+fkRecoParam(recoParam),
+fMaxMCSAngle2(0.)
 {
   /// Constructor for class AliMUONVTrackReconstructor
   /// WARNING: if clusterServer=0x0, no clusterization will be possible at this level
@@ -110,6 +111,13 @@ fkRecoParam(recoParam)
   
   // set the magnetic field for track extrapolations
   AliMUONTrackExtrap::SetField();
+  
+  // set the maximum MCS angle in chamber from the minimum acceptable momentum
+  AliMUONTrackParam param;
+  Double_t inverseBendingP = (GetRecoParam()->GetMinBendingMomentum() > 0.) ? 1./GetRecoParam()->GetMinBendingMomentum() : 1.;
+  param.SetInverseBendingMomentum(inverseBendingP);
+  fMaxMCSAngle2 = AliMUONTrackExtrap::GetMCSAngle2(param, AliMUONConstants::ChamberThicknessInX0(), 1.);
+  
 }
 
   //__________________________________________________________________________
@@ -176,17 +184,93 @@ void AliMUONVTrackReconstructor::EventReconstruct(AliMUONVClusterStore& clusterS
   }
 }
 
-  //__________________________________________________________________________
+//__________________________________________________________________________
+Bool_t AliMUONVTrackReconstructor::IsAcceptable(AliMUONTrackParam &trackParam)
+{
+  /// Return kTRUE if the track is within given limits on momentum/angle/origin
+  
+  const TMatrixD& kParamCov = trackParam.GetCovariances();
+  Int_t chamber = trackParam.GetClusterPtr()->GetChamberId();
+  Double_t z = trackParam.GetZ();
+  Double_t sigmaCut = GetRecoParam()->GetSigmaCutForTracking();
+  
+  // MCS dipersion
+  Double_t angleMCS2 = 0.;
+  if (AliMUONTrackExtrap::IsFieldON() && chamber < 6)
+    angleMCS2 = AliMUONTrackExtrap::GetMCSAngle2(trackParam, AliMUONConstants::ChamberThicknessInX0(), 1.);
+  else angleMCS2 = fMaxMCSAngle2;
+  Double_t impactMCS2 = 0.;
+  if (!GetRecoParam()->SelectOnTrackSlope()) for (Int_t iCh=0; iCh<=chamber; iCh++)
+    impactMCS2 += AliMUONConstants::DefaultChamberZ(chamber) * AliMUONConstants::DefaultChamberZ(chamber) * angleMCS2;
+  
+  // ------ track selection in non bending direction ------
+  if (GetRecoParam()->SelectOnTrackSlope()) {
+    
+    // check if non bending slope is within tolerances
+    Double_t nonBendingSlopeErr = TMath::Sqrt(kParamCov(1,1) + (chamber + 1.) * angleMCS2);
+    if ((TMath::Abs(trackParam.GetNonBendingSlope()) - sigmaCut * nonBendingSlopeErr) > GetRecoParam()->GetMaxNonBendingSlope()) return kFALSE;
+    
+  } else {
+    
+    // or check if non bending impact parameter is within tolerances
+    Double_t nonBendingImpactParam = TMath::Abs(trackParam.GetNonBendingCoor() - z * trackParam.GetNonBendingSlope());
+    Double_t nonBendingImpactParamErr = TMath::Sqrt(kParamCov(0,0) + z * z * kParamCov(1,1) - 2. * z * kParamCov(0,1) + impactMCS2);
+    if ((nonBendingImpactParam - sigmaCut * nonBendingImpactParamErr) > (3. * GetRecoParam()->GetNonBendingVertexDispersion())) return kFALSE;
+    
+  }
+  
+  // ------ track selection in bending direction ------
+  if (AliMUONTrackExtrap::IsFieldON()) { // depending whether the field is ON or OFF
+    
+    // check if bending momentum is within tolerances
+    Double_t bendingMomentum = TMath::Abs(1. / trackParam.GetInverseBendingMomentum());
+    Double_t bendingMomentumErr = TMath::Sqrt(kParamCov(4,4)) * bendingMomentum * bendingMomentum;
+    if (chamber < 6 && (bendingMomentum + sigmaCut * bendingMomentumErr) < GetRecoParam()->GetMinBendingMomentum()) return kFALSE;
+    else if ((bendingMomentum + 3. * bendingMomentumErr) < GetRecoParam()->GetMinBendingMomentum()) return kFALSE;
+    
+  } else {
+    
+    if (GetRecoParam()->SelectOnTrackSlope()) {
+      
+      // check if bending slope is within tolerances
+      Double_t bendingSlopeErr = TMath::Sqrt(kParamCov(3,3) + (chamber + 1.) * angleMCS2);
+      if ((TMath::Abs(trackParam.GetBendingSlope()) - sigmaCut * bendingSlopeErr) > GetRecoParam()->GetMaxBendingSlope()) return kFALSE;
+      
+    } else {
+      
+      // or check if bending impact parameter is within tolerances
+      Double_t bendingImpactParam = TMath::Abs(trackParam.GetBendingCoor() - z * trackParam.GetBendingSlope());
+      Double_t bendingImpactParamErr = TMath::Sqrt(kParamCov(2,2) + z * z * kParamCov(3,3) - 2. * z * kParamCov(2,3) + impactMCS2);
+      if ((bendingImpactParam - sigmaCut * bendingImpactParamErr) > (3. * GetRecoParam()->GetBendingVertexDispersion())) return kFALSE;
+      
+    }
+    
+  }
+  
+  return kTRUE;
+  
+}
+
+//__________________________________________________________________________
 TClonesArray* AliMUONVTrackReconstructor::MakeSegmentsBetweenChambers(const AliMUONVClusterStore& clusterStore, Int_t ch1, Int_t ch2)
 {
   /// To make the list of segments from the list of clusters in the 2 given chambers.
   /// Return a new TClonesArray of segments.
   /// It is the responsibility of the user to delete it afterward.
   AliDebug(1,Form("Enter MakeSegmentsBetweenChambers (1..) %d-%d", ch1+1, ch2+1));
+  AliCodeTimerAuto("");
   
   AliMUONVCluster *cluster1, *cluster2;
   AliMUONObjectPair *segment;
-  Double_t nonBendingSlope = 0, bendingSlope = 0, impactParam = 0., bendingMomentum = 0.; // to avoid compilation warning
+  Double_t z1 = 0., z2 = 0., dZ = 0.;
+  Double_t nonBendingSlope = 0., nonBendingSlopeErr = 0., nonBendingImpactParam = 0., nonBendingImpactParamErr = 0.;
+  Double_t bendingSlope = 0., bendingSlopeErr = 0., bendingImpactParam = 0., bendingImpactParamErr = 0., bendingImpactParamErr2 = 0.;
+  Double_t bendingMomentum = 0., bendingMomentumErr = 0.;
+  Double_t bendingVertexDispersion2 = GetRecoParam()->GetBendingVertexDispersion() * GetRecoParam()->GetBendingVertexDispersion();
+  Double_t impactMCS2 = 0; // maximum impact parameter dispersion**2 due to MCS in chamber
+  if (!GetRecoParam()->SelectOnTrackSlope() || AliMUONTrackExtrap::IsFieldON()) for (Int_t iCh=0; iCh<=ch1; iCh++)
+    impactMCS2 += AliMUONConstants::DefaultChamberZ(iCh) * AliMUONConstants::DefaultChamberZ(iCh) * fMaxMCSAngle2;
+  Double_t sigmaCut = GetRecoParam()->GetSigmaCutForTracking();
   
   // Create iterators to loop over clusters in both chambers
   TIter nextInCh1(clusterStore.CreateChamberIterator(ch1,ch1));
@@ -197,41 +281,64 @@ TClonesArray* AliMUONVTrackReconstructor::MakeSegmentsBetweenChambers(const AliM
   
   // Loop over clusters in the first chamber of the station
   while ( ( cluster1 = static_cast<AliMUONVCluster*>(nextInCh1()) ) ) {
+    z1 = cluster1->GetZ();
     
     // reset cluster iterator of chamber 2
     nextInCh2.Reset();
     
     // Loop over clusters in the second chamber of the station
     while ( ( cluster2 = static_cast<AliMUONVCluster*>(nextInCh2()) ) ) {
+      z2 = cluster2->GetZ();
+      dZ = z1 - z2;
       
-      // non bending slope
-      nonBendingSlope = (cluster1->GetX() - cluster2->GetX()) / (cluster1->GetZ() - cluster2->GetZ());
-      
-      // check if non bending slope is within tolerances
-      if (TMath::Abs(nonBendingSlope) > GetRecoParam()->GetMaxNonBendingSlope()) continue;
-      
-      // bending slope
-      bendingSlope = (cluster1->GetY() - cluster2->GetY()) / (cluster1->GetZ() - cluster2->GetZ());
-      
-      // check the bending momentum of the bending slope depending if the field is ON or OFF
-      if (AliMUONTrackExtrap::IsFieldON()) {
+      // ------ track selection in non bending direction ------
+      nonBendingSlope = (cluster1->GetX() - cluster2->GetX()) / dZ;
+      if (GetRecoParam()->SelectOnTrackSlope()) {
+       
+       // check if non bending slope is within tolerances
+       nonBendingSlopeErr = TMath::Sqrt((cluster1->GetErrX2() + cluster2->GetErrX2()) / dZ / dZ + (ch1 + 1.) * fMaxMCSAngle2);
+       if ((TMath::Abs(nonBendingSlope) - sigmaCut * nonBendingSlopeErr) > GetRecoParam()->GetMaxNonBendingSlope()) continue;
+       
+      } else {
        
-       // impact parameter
-       impactParam = cluster1->GetY() - cluster1->GetZ() * bendingSlope;
+       // or check if non bending impact parameter is within tolerances
+       nonBendingImpactParam = TMath::Abs(cluster1->GetX() - cluster1->GetZ() * nonBendingSlope);
+       nonBendingImpactParamErr = TMath::Sqrt((z1 * z1 * cluster2->GetErrX2() + z2 * z2 * cluster1->GetErrX2()) / dZ / dZ + impactMCS2);
+       if ((nonBendingImpactParam - sigmaCut * nonBendingImpactParamErr) > (3. * GetRecoParam()->GetNonBendingVertexDispersion())) continue;
        
-       // absolute value of bending momentum
-       bendingMomentum = TMath::Abs(AliMUONTrackExtrap::GetBendingMomentumFromImpactParam(impactParam));
+      }
+      
+      // ------ track selection in bending direction ------
+      bendingSlope = (cluster1->GetY() - cluster2->GetY()) / dZ;
+      if (AliMUONTrackExtrap::IsFieldON()) { // depending whether the field is ON or OFF
        
        // check if bending momentum is within tolerances
-       if (bendingMomentum < GetRecoParam()->GetMinBendingMomentum() ||
-           bendingMomentum > GetRecoParam()->GetMaxBendingMomentum()) continue;
+       bendingImpactParam = cluster1->GetY() - cluster1->GetZ() * bendingSlope;
+       bendingImpactParamErr2 = (z1 * z1 * cluster2->GetErrY2() + z2 * z2 * cluster1->GetErrY2()) / dZ / dZ + impactMCS2;
+       bendingMomentum = TMath::Abs(AliMUONTrackExtrap::GetBendingMomentumFromImpactParam(bendingImpactParam));
+       bendingMomentumErr = TMath::Sqrt((bendingVertexDispersion2 + bendingImpactParamErr2) /
+                                        bendingImpactParam / bendingImpactParam + 0.01) * bendingMomentum;
+       if ((bendingMomentum + 3. * bendingMomentumErr) < GetRecoParam()->GetMinBendingMomentum()) continue;
        
       } else {
        
-       // check if non bending slope is within tolerances
-       if (TMath::Abs(bendingSlope) > GetRecoParam()->GetMaxBendingSlope()) continue;
-      
+       if (GetRecoParam()->SelectOnTrackSlope()) {
+         
+         // check if bending slope is within tolerances
+         bendingSlopeErr = TMath::Sqrt((cluster1->GetErrY2() + cluster2->GetErrY2()) / dZ / dZ + (ch1 + 1.) * fMaxMCSAngle2);
+         if ((TMath::Abs(bendingSlope) - sigmaCut * bendingSlopeErr) > GetRecoParam()->GetMaxBendingSlope()) continue;
+         
+       } else {
+         
+         // or check if bending impact parameter is within tolerances
+         bendingImpactParam = TMath::Abs(cluster1->GetY() - cluster1->GetZ() * bendingSlope);
+         bendingImpactParamErr = TMath::Sqrt((z1 * z1 * cluster2->GetErrY2() + z2 * z2 * cluster1->GetErrY2()) / dZ / dZ + impactMCS2);
+         if ((bendingImpactParam - sigmaCut * bendingImpactParamErr) > (3. * GetRecoParam()->GetBendingVertexDispersion())) continue;
+         
+       }
+       
       }
+      
       // make new segment
       segment = new ((*segments)[segments->GetLast()+1]) AliMUONObjectPair(cluster1, cluster2, kFALSE, kFALSE);
       
@@ -469,16 +576,18 @@ void AliMUONVTrackReconstructor::AskForNewClustersInChamber(const AliMUONTrackPa
   AliMUONTrackParam extrapTrackParam(trackParam);
   AliMUONTrackExtrap::ExtrapToZCov(&extrapTrackParam, AliMUONConstants::DefaultChamberZ(chamber));
   
-  // build the searching area using the track resolution and the maximum-distance-to-track value
+  // build the searching area using the track and chamber resolutions and the maximum-distance-to-track value
   const TMatrixD& kParamCov = extrapTrackParam.GetCovariances();
-  Double_t errX2 = kParamCov(0,0) + kMaxDZ * kMaxDZ * kParamCov(1,1) + 2. * kMaxDZ * TMath::Abs(kParamCov(0,1));
-  Double_t errY2 = kParamCov(2,2) + kMaxDZ * kMaxDZ * kParamCov(3,3) + 2. * kMaxDZ * TMath::Abs(kParamCov(2,3));
+  Double_t errX2 = kParamCov(0,0) + kMaxDZ * kMaxDZ * kParamCov(1,1) + 2. * kMaxDZ * TMath::Abs(kParamCov(0,1)) +
+                   GetRecoParam()->GetDefaultNonBendingReso(chamber) * GetRecoParam()->GetDefaultNonBendingReso(chamber);
+  Double_t errY2 = kParamCov(2,2) + kMaxDZ * kMaxDZ * kParamCov(3,3) + 2. * kMaxDZ * TMath::Abs(kParamCov(2,3)) +
+                  GetRecoParam()->GetDefaultBendingReso(chamber) * GetRecoParam()->GetDefaultBendingReso(chamber);
   Double_t dX = TMath::Abs(trackParam.GetNonBendingSlope()) * kMaxDZ +
                GetRecoParam()->GetMaxNonBendingDistanceToTrack() +
-               GetRecoParam()->GetSigmaCutForTracking() * TMath::Sqrt(errX2);
+               GetRecoParam()->GetSigmaCutForTracking() * TMath::Sqrt(2. * errX2);
   Double_t dY = TMath::Abs(trackParam.GetBendingSlope()) * kMaxDZ +
                GetRecoParam()->GetMaxBendingDistanceToTrack() +
-               GetRecoParam()->GetSigmaCutForTracking() * TMath::Sqrt(errY2);
+               GetRecoParam()->GetSigmaCutForTracking() * TMath::Sqrt(2. * errY2);
   AliMpArea area(extrapTrackParam.GetNonBendingCoor(), 
                  extrapTrackParam.GetBendingCoor(),
                  dX, dY);
@@ -521,9 +630,12 @@ Double_t AliMUONVTrackReconstructor::TryOneCluster(const AliMUONTrackParam &trac
   const TMatrixD& kParamCov = trackParamAtCluster.GetCovariances();
   Double_t sigmaX2 = kParamCov(0,0) + cluster->GetErrX2();
   Double_t sigmaY2 = kParamCov(2,2) + cluster->GetErrY2();
+  Double_t covXY   = kParamCov(0,2);
+  Double_t det     = sigmaX2 * sigmaY2 - covXY * covXY;
   
   // Compute chi2
-  return dX * dX / sigmaX2 + dY * dY / sigmaY2;
+  if (det == 0.) return 1.e10;
+  return (dX * dX * sigmaY2 + dY * dY * sigmaX2 - 2. * dX * dY * covXY) / det;
   
 }
 
@@ -531,7 +643,7 @@ Double_t AliMUONVTrackReconstructor::TryOneCluster(const AliMUONTrackParam &trac
 Bool_t AliMUONVTrackReconstructor::TryOneClusterFast(const AliMUONTrackParam &trackParam, const AliMUONVCluster* cluster)
 {
 /// Test the compatibility between the track and the cluster
-/// given the track resolution + the maximum-distance-to-track value
+/// given the track and cluster resolutions + the maximum-distance-to-track value
 /// and assuming linear propagation of the track:
 /// return kTRUE if they are compatibles
   
@@ -539,12 +651,12 @@ Bool_t AliMUONVTrackReconstructor::TryOneClusterFast(const AliMUONTrackParam &tr
   Double_t dX = cluster->GetX() - (trackParam.GetNonBendingCoor() + trackParam.GetNonBendingSlope() * dZ);
   Double_t dY = cluster->GetY() - (trackParam.GetBendingCoor() + trackParam.GetBendingSlope() * dZ);
   const TMatrixD& kParamCov = trackParam.GetCovariances();
-  Double_t errX2 = kParamCov(0,0) + dZ * dZ * kParamCov(1,1) + 2. * dZ * kParamCov(0,1);
-  Double_t errY2 = kParamCov(2,2) + dZ * dZ * kParamCov(3,3) + 2. * dZ * kParamCov(2,3);
+  Double_t errX2 = kParamCov(0,0) + dZ * dZ * kParamCov(1,1) + 2. * dZ * kParamCov(0,1) + cluster->GetErrX2();
+  Double_t errY2 = kParamCov(2,2) + dZ * dZ * kParamCov(3,3) + 2. * dZ * kParamCov(2,3) + cluster->GetErrY2();
 
-  Double_t dXmax = GetRecoParam()->GetSigmaCutForTracking() * TMath::Sqrt(errX2) +
+  Double_t dXmax = GetRecoParam()->GetSigmaCutForTracking() * TMath::Sqrt(2. * errX2) +
                    GetRecoParam()->GetMaxNonBendingDistanceToTrack();
-  Double_t dYmax = GetRecoParam()->GetSigmaCutForTracking() * TMath::Sqrt(errY2) +
+  Double_t dYmax = GetRecoParam()->GetSigmaCutForTracking() * TMath::Sqrt(2. * errY2) +
                   GetRecoParam()->GetMaxBendingDistanceToTrack();
   
   if (TMath::Abs(dX) > dXmax || TMath::Abs(dY) > dYmax) return kFALSE;
index 4b25fe1..09b630d 100644 (file)
@@ -64,6 +64,8 @@ class AliMUONVTrackReconstructor : public TObject {
 
   const AliMUONRecoParam* fkRecoParam; ///< reference to reco parameters
   
+  Double_t fMaxMCSAngle2; ///< maximum angle dispersion due to MCS
+  
   // Functions
   AliMUONVTrackReconstructor (const AliMUONVTrackReconstructor& rhs); ///< copy constructor
   AliMUONVTrackReconstructor& operator=(const AliMUONVTrackReconstructor& rhs); ///< assignment operator
@@ -83,6 +85,8 @@ class AliMUONVTrackReconstructor : public TObject {
   /// Finalize the given track
   virtual void FinalizeTrack(AliMUONTrack &track) = 0;
   
+  Bool_t IsAcceptable(AliMUONTrackParam &trackParam);
+  
   TClonesArray* MakeSegmentsBetweenChambers(const AliMUONVClusterStore& clusterStore, Int_t ch1, Int_t ch2);
 
   void RemoveUsedSegments(TClonesArray& segments);
index 8720798..c8d1d42 100644 (file)
@@ -82,8 +82,8 @@ void runDataReconstruction(const char* input = "/Users/laurent/Alice/Data/Raw/09
   muonRecoParam->SetCalibrationMode(caliboption.Data());
   
   // cut on (non)bending slopes
-  muonRecoParam->SetMaxNonBendingSlope(0.6);
-  muonRecoParam->SetMaxBendingSlope(0.6);
+  muonRecoParam->SetMaxNonBendingSlope(0.4);
+  muonRecoParam->SetMaxBendingSlope(0.5);
   
   // tracking algorithm
   muonRecoParam->MakeMoreTrackCandidates(kTRUE);
index e7ebd5c..15130b1 100644 (file)
Binary files a/OCDB/MUON/Calib/RecoParam/Run0_999999999_v0_s0.root and b/OCDB/MUON/Calib/RecoParam/Run0_999999999_v0_s0.root differ