X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=MUON%2FAliMUONTrackReconstructor.cxx;h=8a41b4c726e3fd9e9d862f464255703a981f5acf;hb=b6852c8483ae899e1ba6130778e740c7a8c19753;hp=5e6937d507e8812bbb344f54a245a37873358d58;hpb=6b092dfc5b92d57f52172c8c08f1c01fa953c513;p=u%2Fmrichter%2FAliRoot.git diff --git a/MUON/AliMUONTrackReconstructor.cxx b/MUON/AliMUONTrackReconstructor.cxx index 5e6937d507e..8a41b4c726e 100644 --- a/MUON/AliMUONTrackReconstructor.cxx +++ b/MUON/AliMUONTrackReconstructor.cxx @@ -15,6 +15,7 @@ /* $Id$ */ +//----------------------------------------------------------------------------- /// \class AliMUONTrackReconstructor /// MUON track reconstructor using the original method /// @@ -23,611 +24,1158 @@ /// /// It contains as methods, among others: /// - MakeTracks to build the tracks -/// +//----------------------------------------------------------------------------- #include "AliMUONTrackReconstructor.h" -#include "AliMUONRecData.h" + #include "AliMUONConstants.h" -#include "AliMUONRawCluster.h" -#include "AliMUONHitForRec.h" -#include "AliMUONObjectPair.h" +#include "AliMUONVCluster.h" +#include "AliMUONVClusterServer.h" +#include "AliMUONVClusterStore.h" #include "AliMUONTrack.h" #include "AliMUONTrackParam.h" #include "AliMUONTrackExtrap.h" +#include "AliMUONRecoParam.h" +#include "AliMUONGeometryTransformer.h" + +#include "AliMpArea.h" #include "AliLog.h" #include #include +#include #include +#include // 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); - -Double_t MultipleScatteringAngle2(AliMUONTrackParam *param); +void TrackChi2(Int_t &nParam, Double_t *gradient, Double_t &chi2, Double_t *param, Int_t flag); +using std::endl; +using std::cout; /// \cond CLASSIMP ClassImp(AliMUONTrackReconstructor) // Class implementation in ROOT context /// \endcond -//************* Defaults parameters for reconstruction -const Double_t AliMUONTrackReconstructor::fgkMaxNormChi2 = 100.0; -const Bool_t AliMUONTrackReconstructor::fgkTrackAllTracks = kFALSE; - -//__________________________________________________________________________ -AliMUONTrackReconstructor::AliMUONTrackReconstructor(AliMUONRecData* data) - : AliMUONVTrackReconstructor(data) -{ - /// Constructor for class AliMUONTrackReconstructor - - // Memory allocation for the TClonesArray of reconstructed tracks - fRecTracksPtr = new TClonesArray("AliMUONTrack", 10); -} - //__________________________________________________________________________ -AliMUONTrackReconstructor::~AliMUONTrackReconstructor(void) +AliMUONTrackReconstructor::AliMUONTrackReconstructor(const AliMUONRecoParam* recoParam, AliMUONVClusterServer* clusterServer, + const AliMUONGeometryTransformer* transformer) + : AliMUONVTrackReconstructor(recoParam, clusterServer, transformer) { - /// Destructor for class AliMUONTrackReconstructor - delete fRecTracksPtr; + /// Constructor } //__________________________________________________________________________ -void AliMUONTrackReconstructor::AddHitsForRecFromRawClusters() +AliMUONTrackReconstructor::~AliMUONTrackReconstructor() { - /// To add to the list of hits for reconstruction all the raw clusters - TTree *treeR; - AliMUONHitForRec *hitForRec; - AliMUONRawCluster *clus; - Int_t iclus, nclus; - TClonesArray *rawclusters; - AliDebug(1,"Enter AddHitsForRecFromRawClusters"); - - treeR = fMUONData->TreeR(); - if (!treeR) { - AliError("TreeR must be loaded"); - exit(0); - } - - fMUONData->SetTreeAddress("RC"); - fMUONData->GetRawClusters(); // only one entry - - // Loop over tracking chambers - for (Int_t ch = 0; ch < AliMUONConstants::NTrackingCh(); ch++) { - rawclusters =fMUONData->RawClusters(ch); - nclus = (Int_t) (rawclusters->GetEntries()); - // Loop over (cathode correlated) raw clusters - for (iclus = 0; iclus < nclus; iclus++) { - clus = (AliMUONRawCluster*) rawclusters->UncheckedAt(iclus); - // new AliMUONHitForRec from raw cluster - // and increment number of AliMUONHitForRec's (total and in chamber) - hitForRec = new ((*fHitsForRecPtr)[fNHitsForRec]) AliMUONHitForRec(clus); - fNHitsForRec++; - // more information into HitForRec - hitForRec->SetBendingReso2(clus->GetErrY() * clus->GetErrY()); - hitForRec->SetNonBendingReso2(clus->GetErrX() * clus->GetErrX()); - // original raw cluster - hitForRec->SetChamberNumber(ch); - hitForRec->SetHitNumber(iclus); - // Z coordinate of the raw cluster (cm) - hitForRec->SetZ(clus->GetZ(0)); - if (AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 3) { - cout << "Chamber " << ch <<" raw cluster " << iclus << " : " << endl; - clus->Print("full"); - cout << "AliMUONHitForRec number (1...): " << fNHitsForRec << endl; - hitForRec->Print("full"); - } - } // end of cluster loop - } // end of chamber loop - SortHitsForRecWithIncreasingChamber(); - - AliDebug(1,"End of AddHitsForRecFromRawClusters"); - if (AliLog::GetGlobalDebugLevel() > 0) { - AliDebug(1, Form("NHitsForRec: %d",fNHitsForRec)); - for (Int_t ch = 0; ch < AliMUONConstants::NTrackingCh(); ch++) { - AliDebug(1, Form("Chamber(0...): %d",ch)); - AliDebug(1, Form("NHitsForRec: %d", fNHitsForRecPerChamber[ch])); - AliDebug(1, Form("Index(first HitForRec): %d", fIndexOfFirstHitForRecPerChamber[ch])); - for (Int_t hit = fIndexOfFirstHitForRecPerChamber[ch]; - hit < fIndexOfFirstHitForRecPerChamber[ch] + fNHitsForRecPerChamber[ch]; - hit++) { - AliDebug(1, Form("HitForRec index(0...): %d",hit)); - ((*fHitsForRecPtr)[hit])->Dump(); - } - } - } - - return; -} - - //__________________________________________________________________________ -void AliMUONTrackReconstructor::MakeTracks(void) -{ - /// To make the tracks from the list of segments and points in all stations - AliDebug(1,"Enter MakeTracks"); - // Look for candidates from at least 3 aligned points in stations(1..) 4 and 5 - MakeTrackCandidates(); - // Follow tracks in stations(1..) 3, 2 and 1 - FollowTracks(); - // Remove double tracks - RemoveDoubleTracks(); - // Fill out the AliMUONTrack's - FillMUONTrack(); -} +/// Destructor +} //__________________________________________________________________________ -void AliMUONTrackReconstructor::MakeTrackCandidates(void) +Bool_t AliMUONTrackReconstructor::MakeTrackCandidates(AliMUONVClusterStore& clusterStore) { - /// To make track candidates: + /// To make track candidates (assuming linear propagation if the flag fgkMakeTrackCandidatesFast is set to kTRUE): /// Start with segments station(1..) 4 or 5 then follow track in station 5 or 4. - /// Good candidates are made of at least three hitForRec's. + /// Good candidates are made of at least three clusters. /// Keep only best candidates or all of them according to the flag fgkTrackAllTracks. + TClonesArray *segments; - AliMUONObjectPair *segment; - AliMUONHitForRec *hitForRec1, *hitForRec2; AliMUONTrack *track; - AliMUONTrackParam *trackParamAtFirstHit; Int_t iCandidate = 0; + Bool_t clusterFound; AliDebug(1,"Enter MakeTrackCandidates"); + // Unless we're doing combined tracking, we'll clusterize all stations at once + Int_t firstChamber(0); + Int_t lastChamber(9); + + if (GetRecoParam()->CombineClusterTrackReco()) { + // ... Here's the exception : ask the clustering to reconstruct + // clusters *only* in station 4 and 5 for combined tracking + firstChamber = 6; + } + + for (Int_t i = firstChamber; i <= lastChamber; ++i ) + { + if (fClusterServer && GetRecoParam()->UseChamber(i)) fClusterServer->Clusterize(i, clusterStore, AliMpArea(), GetRecoParam()); + } + // Loop over stations(1..) 5 and 4 and make track candidates for (Int_t istat=4; istat>=3; istat--) { + // Make segments in the station - segments = MakeSegmentsInStation(istat); + segments = MakeSegmentsBetweenChambers(clusterStore, 2*istat, 2*istat+1); + // Loop over segments - for (Int_t iseg=0; isegGetEntriesFast(); iseg++) { + for (Int_t iseg=0; isegGetEntriesFast(); iseg++) + { AliDebug(1,Form("Making primary candidate(1..) %d",++iCandidate)); + // Transform segments to tracks and put them at the end of fRecTracksPtr - segment = (AliMUONObjectPair*) ((*segments)[iseg]); - hitForRec1 = (AliMUONHitForRec*) segment->First(); - hitForRec2 = (AliMUONHitForRec*) segment->Second(); - track = new ((*fRecTracksPtr)[fRecTracksPtr->GetLast()+1]) AliMUONTrack(hitForRec1, hitForRec2); + track = new ((*fRecTracksPtr)[fRecTracksPtr->GetLast()+1]) AliMUONTrack((AliMUONObjectPair*)((*segments)[iseg]),GetRecoParam()->GetBendingVertexDispersion()); fNRecTracks++; - // Add MCS effects in parameter covariances - trackParamAtFirstHit = (AliMUONTrackParam*) (track->GetTrackParamAtHit()->First()); - AliMUONTrackExtrap::AddMCSEffect(trackParamAtFirstHit,AliMUONConstants::ChamberThicknessInX0(),1.); - // Printout for debuging - if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 2) || (AliLog::GetGlobalDebugLevel() >= 2)) { - cout<GetCovariances()->Print(); + + // Look for compatible cluster(s) in the other station + if (GetRecoParam()->MakeTrackCandidatesFast()) + clusterFound = FollowLinearTrackInStation(*track, clusterStore, 7-istat); + else clusterFound = FollowTrackInStation(*track, clusterStore, 7-istat); + + // Remove track if no cluster found on a requested station + // or abort tracking if there are too many candidates + if (GetRecoParam()->RequestStation(7-istat)) { + if (!clusterFound) { + fRecTracksPtr->Remove(track); + fNRecTracks--; + } else if (fNRecTracks > GetRecoParam()->GetMaxTrackCandidates()) { + AliError(Form("Too many track candidates (%d tracks). Stop tracking.", fNRecTracks)); + return kFALSE; + } + } else { + if ((fNRecTracks + segments->GetEntriesFast() - iseg - 1) > GetRecoParam()->GetMaxTrackCandidates()) { + AliError(Form("Too many track candidates (%d tracks). Stop tracking.", fNRecTracks + segments->GetEntriesFast() - iseg - 1)); + return kFALSE; + } } - // Look for compatible hitForRec(s) in the other station - FollowTrackInStation(track,7-istat); + } - // delete the array of segments - delete segments; + } - fRecTracksPtr->Compress(); // this is essential before checking tracks // Keep all different tracks or only the best ones as required - if (fgkTrackAllTracks) RemoveIdenticalTracks(); + if (GetRecoParam()->TrackAllTracks()) RemoveIdenticalTracks(); else RemoveDoubleTracks(); AliDebug(1,Form("Number of good candidates = %d",fNRecTracks)); + return kTRUE; + } //__________________________________________________________________________ -void AliMUONTrackReconstructor::RemoveIdenticalTracks(void) -{ - /// To remove identical tracks: - /// Tracks are considered identical if they have all their hits in common. - /// One keeps the track with the larger number of hits if need be - AliMUONTrack *track1, *track2, *trackToRemove; - Int_t hitsInCommon, nHits1, nHits2; - Bool_t removedTrack1; - // Loop over first track of the pair - track1 = (AliMUONTrack*) fRecTracksPtr->First(); - while (track1) { - removedTrack1 = kFALSE; - nHits1 = track1->GetNTrackHits(); - // Loop over second track of the pair - track2 = (AliMUONTrack*) fRecTracksPtr->After(track1); - while (track2) { - nHits2 = track2->GetNTrackHits(); - // number of hits in common between two tracks - hitsInCommon = track1->HitsInCommon(track2); - // check for identical tracks - if ((hitsInCommon == nHits1) || (hitsInCommon == nHits2)) { - // decide which track to remove - if (nHits2 > nHits1) { - // remove track1 and continue the first loop with the track next to track1 - trackToRemove = track1; - track1 = (AliMUONTrack*) fRecTracksPtr->After(track1); - fRecTracksPtr->Remove(trackToRemove); - fRecTracksPtr->Compress(); // this is essential to retrieve the TClonesArray afterwards - fNRecTracks--; - removedTrack1 = kTRUE; - break; - } else { - // remove track2 and continue the second loop with the track next to track2 - trackToRemove = track2; - track2 = (AliMUONTrack*) fRecTracksPtr->After(track2); - fRecTracksPtr->Remove(trackToRemove); - fRecTracksPtr->Compress(); // this is essential to retrieve the TClonesArray afterwards - fNRecTracks--; - } - } else track2 = (AliMUONTrack*) fRecTracksPtr->After(track2); - } // track2 - if (removedTrack1) continue; - track1 = (AliMUONTrack*) fRecTracksPtr->After(track1); - } // track1 - return; -} - - //__________________________________________________________________________ -void AliMUONTrackReconstructor::RemoveDoubleTracks(void) +Bool_t AliMUONTrackReconstructor::MakeMoreTrackCandidates(AliMUONVClusterStore& clusterStore) { - /// To remove double tracks: - /// Tracks are considered identical if more than half of the hits of the track - /// which has the smaller number of hits are in common with the other track. - /// Among two identical tracks, one keeps the track with the larger number of hits - /// or, if these numbers are equal, the track with the minimum chi2. - AliMUONTrack *track1, *track2, *trackToRemove; - Int_t hitsInCommon, nHits1, nHits2; - Bool_t removedTrack1; - // Loop over first track of the pair - track1 = (AliMUONTrack*) fRecTracksPtr->First(); - while (track1) { - removedTrack1 = kFALSE; - nHits1 = track1->GetNTrackHits(); - // Loop over second track of the pair - track2 = (AliMUONTrack*) fRecTracksPtr->After(track1); - while (track2) { - nHits2 = track2->GetNTrackHits(); - // number of hits in common between two tracks - hitsInCommon = track1->HitsInCommon(track2); - // check for identical tracks - if (((nHits1 < nHits2) && (2 * hitsInCommon > nHits1)) || (2 * hitsInCommon > nHits2)) { - // decide which track to remove - if ((nHits1 > nHits2) || ((nHits1 == nHits2) && (track1->GetFitFMin() <= track2->GetFitFMin()))) { - // remove track2 and continue the second loop with the track next to track2 - trackToRemove = track2; - track2 = (AliMUONTrack*) fRecTracksPtr->After(track2); - fRecTracksPtr->Remove(trackToRemove); - fRecTracksPtr->Compress(); // this is essential to retrieve the TClonesArray afterwards - fNRecTracks--; - } else { - // else remove track1 and continue the first loop with the track next to track1 - trackToRemove = track1; - track1 = (AliMUONTrack*) fRecTracksPtr->After(track1); - fRecTracksPtr->Remove(trackToRemove); - fRecTracksPtr->Compress(); // this is essential to retrieve the TClonesArray afterwards - fNRecTracks--; - removedTrack1 = kTRUE; - break; - } - } else track2 = (AliMUONTrack*) fRecTracksPtr->After(track2); - } // track2 - if (removedTrack1) continue; - track1 = (AliMUONTrack*) fRecTracksPtr->After(track1); - } // track1 - return; + /// To make extra track candidates (assuming linear propagation if the flag fgkMakeTrackCandidatesFast is set to kTRUE): + /// clustering is supposed to be already done + /// Start with segments made of 1 cluster in each of the stations 4 and 5 then follow track in remaining chambers. + /// Good candidates are made of at least three clusters if both station are requested (two otherwise). + /// Keep only best candidates or all of them according to the flag fgkTrackAllTracks. + + TClonesArray *segments; + AliMUONObjectPair *segment; + AliMUONTrack *track; + Int_t iCandidate = 0, iCurrentTrack, nCurrentTracks; + Bool_t clusterFound; + + AliDebug(1,"Enter MakeMoreTrackCandidates"); + + // Double loop over chambers in stations(1..) 4 and 5 to make track candidates + for (Int_t ich1 = 6; ich1 <= 7; ich1++) { + for (Int_t ich2 = 8; ich2 <= 9; ich2++) { + + // Make segments in the station + segments = MakeSegmentsBetweenChambers(clusterStore, ich1, ich2); + + /// Remove segments already attached to a track + RemoveUsedSegments(*segments); + + // Loop over segments + for (Int_t iSegment=0; iSegmentGetEntriesFast(); iSegment++) + { + AliDebug(1,Form("Making primary candidate(1..) %d",++iCandidate)); + segment = (AliMUONObjectPair*) segments->UncheckedAt(iSegment); + + // Transform segments to tracks and put them at the end of fRecTracksPtr + iCurrentTrack = fRecTracksPtr->GetLast()+1; + track = new ((*fRecTracksPtr)[iCurrentTrack]) AliMUONTrack(segment,GetRecoParam()->GetBendingVertexDispersion()); + fNRecTracks++; + + // Look for compatible cluster(s) in the second chamber of station 5 + clusterFound = FollowLinearTrackInChamber(*track, clusterStore, 17-ich2); + + // skip the original track in case it has been removed + if (GetRecoParam()->TrackAllTracks() && clusterFound) iCurrentTrack++; + + // loop over every new tracks + nCurrentTracks = fRecTracksPtr->GetLast()+1; + while (iCurrentTrack < nCurrentTracks) { + track = (AliMUONTrack*) fRecTracksPtr->UncheckedAt(iCurrentTrack); + + // Look for compatible cluster(s) in the second chamber of station 4 + FollowLinearTrackInChamber(*track, clusterStore, 13-ich1); + + iCurrentTrack++; + } + + // abort tracking if there are too many candidates + if ((fNRecTracks + segments->GetEntriesFast() - iSegment - 1) > GetRecoParam()->GetMaxTrackCandidates()) { + AliError(Form("Too many track candidates (%d tracks). Stop tracking.", fNRecTracks + segments->GetEntriesFast() - iSegment - 1)); + return kFALSE; + } + + } + + } + } + + // Keep only the best tracks if required + if (!GetRecoParam()->TrackAllTracks()) RemoveDoubleTracks(); + else fRecTracksPtr->Compress(); + + AliDebug(1,Form("Number of good candidates = %d",fNRecTracks)); + + return kTRUE; + } //__________________________________________________________________________ -void AliMUONTrackReconstructor::FollowTracks(void) +Bool_t AliMUONTrackReconstructor::FollowTracks(AliMUONVClusterStore& clusterStore) { /// Follow tracks in stations(1..) 3, 2 and 1 AliDebug(1,"Enter FollowTracks"); AliMUONTrack *track, *nextTrack; - AliMUONTrackParam *trackParamAtFirstHit; - Double_t numberOfDegFree, chi2Norm; + AliMUONTrackParam *trackParam, *nextTrackParam; Int_t currentNRecTracks; + Double_t sigmaCut2 = GetRecoParam()->GetSigmaCutForTracking() * + GetRecoParam()->GetSigmaCutForTracking(); + for (Int_t station = 2; station >= 0; station--) { + // Save the actual number of reconstructed track in case of // tracks are added or suppressed during the tracking procedure // !! Do not compress fRecTracksPtr until the end of the loop over tracks !! currentNRecTracks = fNRecTracks; + for (Int_t iRecTrack = 0; iRecTrack UncheckedAt(iRecTrack); + // Fit the track: // Do not take into account the multiple scattering to speed up the fit // Calculate the track parameter covariance matrix - // If "station" is station(1..) 3 then use the vertex to better constrain the fit - if (station==2) { - SetVertexForFit(track); - track->SetFitWithVertex(kTRUE); - } else track->SetFitWithVertex(kFALSE); - Fit(track,kFALSE, kTRUE); - // Remove the track if the normalized chi2 is too high - numberOfDegFree = (2. * track->GetNTrackHits() - 5.); - if (numberOfDegFree > 0) chi2Norm = track->GetFitFMin() / numberOfDegFree; - else chi2Norm = 1.e10; - if (chi2Norm > fgkMaxNormChi2) { - fRecTracksPtr->Remove(track); - fNRecTracks--; - continue; + // If there is no cluster out of station 4 or 5 then use the vertex to better constrain the fit + if (((AliMUONTrackParam*) track->GetTrackParamAtCluster()->First())->GetClusterPtr()->GetChamberId() > 5) + Fit(*track, kFALSE, kTRUE, kTRUE); + else Fit(*track, kFALSE, kFALSE, kTRUE); + + // remove tracks out of limits + if (!IsAcceptable(*((AliMUONTrackParam*)track->GetTrackParamAtCluster()->First()))) { + fRecTracksPtr->Remove(track); + fNRecTracks--; + continue; } - // Add MCS effects in parameter covariances - trackParamAtFirstHit = (AliMUONTrackParam*) (track->GetTrackParamAtHit()->First()); - AliMUONTrackExtrap::AddMCSEffect(trackParamAtFirstHit,AliMUONConstants::ChamberThicknessInX0(),1.); - // Printout for debuging - if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 2) || (AliLog::GetGlobalDebugLevel() >= 2)) { - cout<GetCovariances()->Print(); + + // remove track if the normalized chi2 is too high + if (track->GetNormalizedChi2() > sigmaCut2) { + fRecTracksPtr->Remove(track); + fNRecTracks--; + continue; + } + + // save parameters from fit into smoothed parameters to complete track afterward + if (GetRecoParam()->ComplementTracks()) { + + if (station==2) { // save track parameters on stations 4 and 5 + + // extrapolate track parameters and covariances at each cluster + // remove the track in case of failure + if (!track->UpdateCovTrackParamAtCluster()) { + fRecTracksPtr->Remove(track); + fNRecTracks--; + continue; + } + + // save them + trackParam = (AliMUONTrackParam*) track->GetTrackParamAtCluster()->First(); + while (trackParam) { + trackParam->SetSmoothParameters(trackParam->GetParameters()); + trackParam->SetSmoothCovariances(trackParam->GetCovariances()); + trackParam = (AliMUONTrackParam*) track->GetTrackParamAtCluster()->After(trackParam); + } + + } else { // or save track parameters on last station only + + trackParam = (AliMUONTrackParam*) track->GetTrackParamAtCluster()->First(); + if (trackParam->GetClusterPtr()->GetChamberId() < 2*(station+2)) { + + // save parameters from fit + trackParam->SetSmoothParameters(trackParam->GetParameters()); + trackParam->SetSmoothCovariances(trackParam->GetCovariances()); + + // save parameters extrapolated to the second chamber of the same station if it has been hit + nextTrackParam = (AliMUONTrackParam*) track->GetTrackParamAtCluster()->After(trackParam); + if (nextTrackParam->GetClusterPtr()->GetChamberId() < 2*(station+2)) { + + // reset parameters and covariances + nextTrackParam->SetParameters(trackParam->GetParameters()); + nextTrackParam->SetZ(trackParam->GetZ()); + nextTrackParam->SetCovariances(trackParam->GetCovariances()); + + // extrapolate them to the z of the corresponding cluster + // remove the track in case of failure + if (!AliMUONTrackExtrap::ExtrapToZCov(nextTrackParam, nextTrackParam->GetClusterPtr()->GetZ())) { + fRecTracksPtr->Remove(track); + fNRecTracks--; + continue; + } + + // save them + nextTrackParam->SetSmoothParameters(nextTrackParam->GetParameters()); + nextTrackParam->SetSmoothCovariances(nextTrackParam->GetCovariances()); + + } + + } + + } + } - // Look for compatible hitForRec in station(0..) "station" - FollowTrackInStation(track,station); + + // Look for compatible cluster(s) in station(0..) "station" + if (!FollowTrackInStation(*track, clusterStore, station)) { + + // Try to recover track if required + if (GetRecoParam()->RecoverTracks()) { + + // work on a copy of the track if this station is not required + // to keep the case where no cluster is reconstructed as a possible candidate + if (!GetRecoParam()->RequestStation(station)) { + track = new ((*fRecTracksPtr)[fRecTracksPtr->GetLast()+1]) AliMUONTrack(*track); + fNRecTracks++; + } + + // try to recover + if (!RecoverTrack(*track, clusterStore, station)) { + // remove track if no cluster found + fRecTracksPtr->Remove(track); + fNRecTracks--; + } + + } else if (GetRecoParam()->RequestStation(station)) { + // remove track if no cluster found + fRecTracksPtr->Remove(track); + fNRecTracks--; + } + + } + + // abort tracking if there are too many candidates + if (fNRecTracks > GetRecoParam()->GetMaxTrackCandidates()) { + AliError(Form("Too many track candidates (%d tracks). Stop tracking.", fNRecTracks)); + return kFALSE; + } + } + // Compress fRecTracksPtr for the next step fRecTracksPtr->Compress(); + // Keep only the best tracks if required - if (!fgkTrackAllTracks) RemoveDoubleTracks(); + if (!GetRecoParam()->TrackAllTracks()) RemoveDoubleTracks(); + } - // Last fit of track candidates with all station + // Last fit of track candidates with all stations // Take into account the multiple scattering and remove bad tracks Int_t trackIndex = -1; track = (AliMUONTrack*) fRecTracksPtr->First(); while (track) { + trackIndex++; nextTrack = (AliMUONTrack*) fRecTracksPtr->After(track); // prepare next track - track->SetFitWithVertex(kFALSE); // just to be sure - Fit(track,kTRUE, kTRUE); + + Fit(*track, kTRUE, kFALSE, kTRUE); + // Printout for debuging if (AliLog::GetGlobalDebugLevel() >= 3) { cout << "FollowTracks: track candidate(0..) " << trackIndex << " after final fit" << endl; track->RecursiveDump(); } + // Remove the track if the normalized chi2 is too high - numberOfDegFree = (2.0 * track->GetNTrackHits() - 5); - if (numberOfDegFree > 0) chi2Norm = track->GetFitFMin() / numberOfDegFree; - else chi2Norm = 1.e10; - if (chi2Norm > fgkMaxNormChi2) { + if (track->GetNormalizedChi2() > sigmaCut2) { fRecTracksPtr->Remove(track); fNRecTracks--; + track = nextTrack; + continue; } + + // save parameters from fit into smoothed parameters to complete track afterward + if (GetRecoParam()->ComplementTracks()) { + + trackParam = (AliMUONTrackParam*) track->GetTrackParamAtCluster()->First(); + if (trackParam->GetClusterPtr()->GetChamberId() < 2) { + + // save parameters from fit + trackParam->SetSmoothParameters(trackParam->GetParameters()); + trackParam->SetSmoothCovariances(trackParam->GetCovariances()); + + // save parameters extrapolated to the second chamber of the same station if it has been hit + nextTrackParam = (AliMUONTrackParam*) track->GetTrackParamAtCluster()->After(trackParam); + if (nextTrackParam->GetClusterPtr()->GetChamberId() < 2) { + + // reset parameters and covariances + nextTrackParam->SetParameters(trackParam->GetParameters()); + nextTrackParam->SetZ(trackParam->GetZ()); + nextTrackParam->SetCovariances(trackParam->GetCovariances()); + + // extrapolate them to the z of the corresponding cluster + // remove the track in case of failure + if (!AliMUONTrackExtrap::ExtrapToZCov(nextTrackParam, nextTrackParam->GetClusterPtr()->GetZ())) { + fRecTracksPtr->Remove(track); + fNRecTracks--; + track = nextTrack; + continue; + } + + // save them + nextTrackParam->SetSmoothParameters(nextTrackParam->GetParameters()); + nextTrackParam->SetSmoothCovariances(nextTrackParam->GetCovariances()); + + } + + } + + } + track = nextTrack; + } + fRecTracksPtr->Compress(); + return kTRUE; + } //__________________________________________________________________________ -void AliMUONTrackReconstructor::FollowTrackInStation(AliMUONTrack* trackCandidate, Int_t nextStation) +Bool_t AliMUONTrackReconstructor::FollowTrackInChamber(AliMUONTrack &trackCandidate, AliMUONVClusterStore& clusterStore, Int_t nextChamber) { - /// Follow trackCandidate in station(0..) nextStation and search for compatible HitForRec(s) + /// Follow trackCandidate in chamber(0..) nextChamber and search for compatible cluster(s) /// Keep all possibilities or only the best one(s) according to the flag fgkTrackAllTracks: /// kTRUE: duplicate "trackCandidate" if there are several possibilities and add the new tracks at the end of /// fRecTracksPtr to avoid conficts with other track candidates at this current stage of the tracking procedure. /// Remove the obsolete "trackCandidate" at the end. - /// kFALSE: add only the best hit(s) to the "trackCandidate". Try to add a couple of hits in priority. + /// kFALSE: add only the best cluster(s) to the "trackCandidate". Try to add a couple of clusters in priority. + AliDebug(1,Form("Enter FollowTrackInChamber(1..) %d", nextChamber+1)); + + Double_t chi2WithOneCluster = AliMUONTrack::MaxChi2(); + Double_t maxChi2WithOneCluster = 2. * GetRecoParam()->GetSigmaCutForTracking() * + GetRecoParam()->GetSigmaCutForTracking(); // 2 because 2 quantities in chi2 + Double_t bestChi2WithOneCluster = maxChi2WithOneCluster; + Bool_t foundOneCluster = kFALSE; + AliMUONTrack *newTrack = 0x0; + AliMUONVCluster *cluster; + AliMUONTrackParam extrapTrackParam; + AliMUONTrackParam extrapTrackParamAtCh; + AliMUONTrackParam extrapTrackParamAtCluster; + AliMUONTrackParam bestTrackParamAtCluster; + + // Get track parameters according to the propagation direction + if (nextChamber > 7) extrapTrackParamAtCh = *(AliMUONTrackParam*)trackCandidate.GetTrackParamAtCluster()->Last(); + else extrapTrackParamAtCh = *(AliMUONTrackParam*)trackCandidate.GetTrackParamAtCluster()->First(); + + // Printout for debuging + if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 2) || (AliLog::GetGlobalDebugLevel() >= 2)) { + cout<GetChamberId(); + AliMUONTrackExtrap::AddMCSEffect(&extrapTrackParamAtCh,AliMUONConstants::ChamberThicknessInX0(currentChamber),-1.); + + // Add MCS in the missing chamber(s) if any + while (currentChamber > nextChamber + 1) { + // extrapolation to the missing chamber + currentChamber--; + if (!AliMUONTrackExtrap::ExtrapToZCov(&extrapTrackParamAtCh, AliMUONConstants::DefaultChamberZ(currentChamber))) return kFALSE; + // add MCS effect + AliMUONTrackExtrap::AddMCSEffect(&extrapTrackParamAtCh,AliMUONConstants::ChamberThicknessInX0(currentChamber),-1.); + } + + //Extrapolate trackCandidate to chamber + if (!AliMUONTrackExtrap::ExtrapToZCov(&extrapTrackParamAtCh, AliMUONConstants::DefaultChamberZ(nextChamber))) return kFALSE; + + // Printout for debuging + if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 2) || (AliLog::GetGlobalDebugLevel() >= 2)) { + cout<= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) { + cout << "FollowTrackInStation: look for clusters in chamber(1..): " << nextChamber+1 << endl; + } + + // Ask the clustering to reconstruct new clusters around the track position in the current chamber + // except for station 4 and 5 that are already entirely clusterized + if (GetRecoParam()->CombineClusterTrackReco()) { + if (nextChamber < 6) AskForNewClustersInChamber(extrapTrackParamAtCh, clusterStore, nextChamber); + } + + // Create iterators to loop over clusters in both chambers + TIter next(clusterStore.CreateChamberIterator(nextChamber,nextChamber)); + + // look for cluster in chamber + while ( ( cluster = static_cast(next()) ) ) { + + // try to add the current cluster fast + if (!TryOneClusterFast(extrapTrackParamAtCh, cluster)) continue; + + // try to add the current cluster accuratly + chi2WithOneCluster = TryOneCluster(extrapTrackParamAtCh, cluster, extrapTrackParamAtCluster); + + // if good chi2 then consider to add cluster + if (chi2WithOneCluster < maxChi2WithOneCluster) { + foundOneCluster = kTRUE; + + // Printout for debuging + if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) { + cout << "FollowTrackInStation: found one cluster in chamber(1..): " << nextChamber+1 + << " (Chi2 = " << chi2WithOneCluster << ")" << endl; + cluster->Print(); + } + + if (GetRecoParam()->TrackAllTracks()) { + // copy trackCandidate into a new track put at the end of fRecTracksPtr and add the new cluster + newTrack = new ((*fRecTracksPtr)[fRecTracksPtr->GetLast()+1]) AliMUONTrack(trackCandidate); + UpdateTrack(*newTrack,extrapTrackParamAtCluster); + fNRecTracks++; + + // Printout for debuging + if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) { + cout << "FollowTrackInStation: added one cluster in chamber(1..): " << nextChamber+1 << endl; + if (AliLog::GetGlobalDebugLevel() >= 3) newTrack->RecursiveDump(); + } + + } else if (chi2WithOneCluster < bestChi2WithOneCluster) { + // keep track of the best single cluster except if a couple of clusters has already been found + bestChi2WithOneCluster = chi2WithOneCluster; + bestTrackParamAtCluster = extrapTrackParamAtCluster; + } + + } + + } + + // fill out the best track if required else clean up the fRecTracksPtr array + if (!GetRecoParam()->TrackAllTracks()) { + if (foundOneCluster) { + UpdateTrack(trackCandidate,bestTrackParamAtCluster); + + // Printout for debuging + if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) { + cout << "FollowTrackInStation: added the best cluster in chamber(1..): " << bestTrackParamAtCluster.GetClusterPtr()->GetChamberId()+1 << endl; + if (AliLog::GetGlobalDebugLevel() >= 3) newTrack->RecursiveDump(); + } + + } else return kFALSE; + + } else if (foundOneCluster) { + + // remove obsolete track + fRecTracksPtr->Remove(&trackCandidate); + fNRecTracks--; + + } else return kFALSE; + + return kTRUE; + +} + + //__________________________________________________________________________ +Bool_t AliMUONTrackReconstructor::FollowTrackInStation(AliMUONTrack &trackCandidate, AliMUONVClusterStore& clusterStore, Int_t nextStation) +{ + /// Follow trackCandidate in station(0..) nextStation and search for compatible cluster(s) + /// Keep all possibilities or only the best one(s) according to the flag fgkTrackAllTracks: + /// kTRUE: duplicate "trackCandidate" if there are several possibilities and add the new tracks at the end of + /// fRecTracksPtr to avoid conficts with other track candidates at this current stage of the tracking procedure. + /// Remove the obsolete "trackCandidate" at the end. + /// kFALSE: add only the best cluster(s) to the "trackCandidate". Try to add a couple of clusters in priority. AliDebug(1,Form("Enter FollowTrackInStation(1..) %d", nextStation+1)); - Int_t ch1 = 2*nextStation; - Int_t ch2 = 2*nextStation+1; - Double_t zCh2 = AliMUONConstants::DefaultChamberZ(ch2); - Double_t chi2WithOneHitForRec = 1.e10; - Double_t chi2WithTwoHitForRec = 1.e10; - Double_t maxChi2WithOneHitForRec = 2.*fgkMaxNormChi2; // 2 because 2 quantities in chi2 - Double_t maxChi2WithTwoHitForRec = 4.*fgkMaxNormChi2; // 4 because 4 quantities in chi2 - Double_t bestChi2WithOneHitForRec = maxChi2WithOneHitForRec; - Double_t bestChi2WithTwoHitForRec = maxChi2WithTwoHitForRec; + // Order the chamber according to the propagation direction (tracking starts with chamber 2): + // - nextStation == station(1...) 5 => forward propagation + // - nextStation < station(1...) 5 => backward propagation + Int_t ch1, ch2; + if (nextStation==4) { + ch1 = 2*nextStation+1; + ch2 = 2*nextStation; + } else { + ch1 = 2*nextStation; + ch2 = 2*nextStation+1; + } + + Double_t chi2WithOneCluster = AliMUONTrack::MaxChi2(); + Double_t chi2WithTwoClusters = AliMUONTrack::MaxChi2(); + Double_t maxChi2WithOneCluster = 2. * GetRecoParam()->GetSigmaCutForTracking() * + GetRecoParam()->GetSigmaCutForTracking(); // 2 because 2 quantities in chi2 + Double_t maxChi2WithTwoClusters = 4. * GetRecoParam()->GetSigmaCutForTracking() * + GetRecoParam()->GetSigmaCutForTracking(); // 4 because 4 quantities in chi2 + Double_t bestChi2WithOneCluster = maxChi2WithOneCluster; + Double_t bestChi2WithTwoClusters = maxChi2WithTwoClusters; + Bool_t foundOneCluster = kFALSE; + Bool_t foundTwoClusters = kFALSE; AliMUONTrack *newTrack = 0x0; - AliMUONHitForRec *hitForRecCh1, *hitForRecCh2; - AliMUONHitForRec *bestHitForRec1 = 0x0, *bestHitForRec2 = 0x0; - Bool_t *hitForRecCh1Used = new Bool_t[fNHitsForRecPerChamber[ch1]]; - for (Int_t hit1 = 0; hit1 < fNHitsForRecPerChamber[ch1]; hit1++) hitForRecCh1Used[hit1] = kFALSE; - // - //Extrapolate trackCandidate to chamber "ch2" to save computing time in the next steps - AliMUONTrackParam *extrapTrackParamPtr = trackCandidate->GetExtrapTrackParam(); - *extrapTrackParamPtr = *((AliMUONTrackParam*)(trackCandidate->GetTrackParamAtHit()->First())); - AliMUONTrackExtrap::ExtrapToZCov(extrapTrackParamPtr, zCh2); - AliMUONTrackParam extrapTrackParamSave(*extrapTrackParamPtr); - // + AliMUONVCluster *clusterCh1, *clusterCh2; + AliMUONTrackParam extrapTrackParam; + AliMUONTrackParam extrapTrackParamAtCh; + AliMUONTrackParam extrapTrackParamAtCluster1; + AliMUONTrackParam extrapTrackParamAtCluster2; + AliMUONTrackParam bestTrackParamAtCluster1; + AliMUONTrackParam bestTrackParamAtCluster2; + + // Get track parameters according to the propagation direction + if (nextStation==4) extrapTrackParamAtCh = *(AliMUONTrackParam*)trackCandidate.GetTrackParamAtCluster()->Last(); + else extrapTrackParamAtCh = *(AliMUONTrackParam*)trackCandidate.GetTrackParamAtCluster()->First(); + + // Printout for debuging + if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 2) || (AliLog::GetGlobalDebugLevel() >= 2)) { + cout<GetChamberId(); + AliMUONTrackExtrap::AddMCSEffect(&extrapTrackParamAtCh,AliMUONConstants::ChamberThicknessInX0(currentChamber),-1.); + + // Add MCS in the missing chamber(s) if any + while (ch1 < ch2 && currentChamber > ch2 + 1) { + // extrapolation to the missing chamber + currentChamber--; + if (!AliMUONTrackExtrap::ExtrapToZCov(&extrapTrackParamAtCh, AliMUONConstants::DefaultChamberZ(currentChamber))) return kFALSE; + // add MCS effect + AliMUONTrackExtrap::AddMCSEffect(&extrapTrackParamAtCh,AliMUONConstants::ChamberThicknessInX0(currentChamber),-1.); + } + + //Extrapolate trackCandidate to chamber "ch2" + if (!AliMUONTrackExtrap::ExtrapToZCov(&extrapTrackParamAtCh, AliMUONConstants::DefaultChamberZ(ch2))) return kFALSE; + // Printout for debuging if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 2) || (AliLog::GetGlobalDebugLevel() >= 2)) { - TMatrixD* paramCovForDebug = extrapTrackParamPtr->GetCovariances(); - cout<Print(); + cout<= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) { - cout << "FollowTrackInStation: look for hits in chamber(1..): " << ch2+1 << endl; + cout << "FollowTrackInStation: look for clusters in chamber(1..): " << ch2+1 << endl; } - for (Int_t hit2 = 0; hit2 < fNHitsForRecPerChamber[ch2]; hit2++) { - hitForRecCh2 = (AliMUONHitForRec*) fHitsForRecPtr->UncheckedAt(fIndexOfFirstHitForRecPerChamber[ch2]+hit2); - // extrapolate track parameters and covariances only once for this hit - AliMUONTrackExtrap::ExtrapToZCov(extrapTrackParamPtr, hitForRecCh2->GetZ()); - chi2WithOneHitForRec = trackCandidate->TryOneHitForRec(hitForRecCh2); - // if good chi2 then try to attach a hitForRec in the other chamber too - if (chi2WithOneHitForRec < maxChi2WithOneHitForRec) { + + // Ask the clustering to reconstruct new clusters around the track position in the current station + // except for station 4 and 5 that are already entirely clusterized + if (GetRecoParam()->CombineClusterTrackReco()) { + if (nextStation < 3) AskForNewClustersInStation(extrapTrackParamAtCh, clusterStore, nextStation); + } + + Int_t nClusters = clusterStore.GetSize(); + Bool_t *clusterCh1Used = new Bool_t[nClusters]; + for (Int_t i = 0; i < nClusters; i++) clusterCh1Used[i] = kFALSE; + Int_t iCluster1; + + // Create iterators to loop over clusters in both chambers + TIter nextInCh1(clusterStore.CreateChamberIterator(ch1,ch1)); + TIter nextInCh2(clusterStore.CreateChamberIterator(ch2,ch2)); + + // look for candidates in chamber 2 + while ( ( clusterCh2 = static_cast(nextInCh2()) ) ) { + + // try to add the current cluster fast + if (!TryOneClusterFast(extrapTrackParamAtCh, clusterCh2)) continue; + + // try to add the current cluster accuratly + chi2WithOneCluster = TryOneCluster(extrapTrackParamAtCh, clusterCh2, extrapTrackParamAtCluster2); + + // if good chi2 then try to attach a cluster in the other chamber too + if (chi2WithOneCluster < maxChi2WithOneCluster) { + Bool_t foundSecondCluster = kFALSE; + // Printout for debuging if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) { - cout << "FollowTrackInStation: look for second hits in chamber(1..): " << ch1+1 << endl; + cout << "FollowTrackInStation: found one cluster in chamber(1..): " << ch2+1 + << " (Chi2 = " << chi2WithOneCluster << ")" << endl; + clusterCh2->Print(); + cout << " look for second clusters in chamber(1..): " << ch1+1 << " ..." << endl; } - Bool_t foundSecondHit = kFALSE; - for (Int_t hit1 = 0; hit1 < fNHitsForRecPerChamber[ch1]; hit1++) { - hitForRecCh1 = (AliMUONHitForRec*) fHitsForRecPtr->UncheckedAt(fIndexOfFirstHitForRecPerChamber[ch1]+hit1); - chi2WithTwoHitForRec = trackCandidate->TryTwoHitForRec(hitForRecCh2, hitForRecCh1); // order hits like that to save computing time - // if good chi2 then create a new track by adding the 2 hitForRec to the "trackCandidate" - if (chi2WithTwoHitForRec < maxChi2WithTwoHitForRec) { - foundSecondHit = kTRUE; - if (fgkTrackAllTracks) { - // copy trackCandidate into a new track put at the end of fRecTracksPtr and add the new hitForRec's - newTrack = new ((*fRecTracksPtr)[fRecTracksPtr->GetLast()+1]) AliMUONTrack(*trackCandidate); + + // add MCS effect for next step + AliMUONTrackExtrap::AddMCSEffect(&extrapTrackParamAtCluster2,AliMUONConstants::ChamberThicknessInX0(ch2),-1.); + + // copy new track parameters for next step + extrapTrackParam = extrapTrackParamAtCluster2; + + //Extrapolate track parameters to chamber "ch1" + Bool_t normalExtrap = AliMUONTrackExtrap::ExtrapToZ(&extrapTrackParam, AliMUONConstants::DefaultChamberZ(ch1)); + + // reset cluster iterator of chamber 1 + nextInCh1.Reset(); + iCluster1 = -1; + + // look for second candidates in chamber 1 + if (normalExtrap) while ( ( clusterCh1 = static_cast(nextInCh1()) ) ) { + iCluster1++; + + // try to add the current cluster fast + if (!TryOneClusterFast(extrapTrackParam, clusterCh1)) continue; + + // try to add the current cluster accurately + chi2WithTwoClusters = TryTwoClusters(extrapTrackParamAtCluster2, clusterCh1, extrapTrackParamAtCluster1); + + // if good chi2 then create a new track by adding the 2 clusters to the "trackCandidate" + if (chi2WithTwoClusters < maxChi2WithTwoClusters) { + foundSecondCluster = kTRUE; + foundTwoClusters = kTRUE; + + // Printout for debuging + if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) { + cout << "FollowTrackInStation: found second cluster in chamber(1..): " << ch1+1 + << " (Global Chi2 = " << chi2WithTwoClusters << ")" << endl; + clusterCh1->Print(); + } + + if (GetRecoParam()->TrackAllTracks()) { + // copy trackCandidate into a new track put at the end of fRecTracksPtr and add the new clusters + newTrack = new ((*fRecTracksPtr)[fRecTracksPtr->GetLast()+1]) AliMUONTrack(trackCandidate); + UpdateTrack(*newTrack,extrapTrackParamAtCluster1,extrapTrackParamAtCluster2); fNRecTracks++; - AliMUONTrackParam trackParam1(extrapTrackParamSave); - AliMUONTrackExtrap::ExtrapToZ(&trackParam1, hitForRecCh1->GetZ()); - newTrack->AddTrackParamAtHit(&trackParam1,hitForRecCh1); - AliMUONTrackParam trackParam2(extrapTrackParamSave); - AliMUONTrackExtrap::ExtrapToZ(&trackParam2, hitForRecCh2->GetZ()); - newTrack->AddTrackParamAtHit(&trackParam2,hitForRecCh2); - // Sort TrackParamAtHit according to increasing Z - newTrack->GetTrackParamAtHit()->Sort(); - // Update the chi2 of the new track - if (newTrack->GetFitFMin()<0) newTrack->SetFitFMin(chi2WithTwoHitForRec); - else newTrack->SetFitFMin(newTrack->GetFitFMin() + chi2WithTwoHitForRec); - // Tag hitForRecCh1 as used - hitForRecCh1Used[hit1] = kTRUE; + + // Tag clusterCh1 as used + clusterCh1Used[iCluster1] = kTRUE; + // Printout for debuging if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) { - cout << "FollowTrackInStation: added two hits in station(1..): " << nextStation+1 - << " (Chi2 = " << chi2WithTwoHitForRec << ")" << endl; + cout << "FollowTrackInStation: added two clusters in station(1..): " << nextStation+1 << endl; if (AliLog::GetGlobalDebugLevel() >= 3) newTrack->RecursiveDump(); } - } else if (chi2WithTwoHitForRec < bestChi2WithTwoHitForRec) { - // keep track of the best couple of hits - bestChi2WithTwoHitForRec = chi2WithTwoHitForRec; - bestHitForRec1 = hitForRecCh1; - bestHitForRec2 = hitForRecCh2; + + } else if (chi2WithTwoClusters < bestChi2WithTwoClusters) { + // keep track of the best couple of clusters + bestChi2WithTwoClusters = chi2WithTwoClusters; + bestTrackParamAtCluster1 = extrapTrackParamAtCluster1; + bestTrackParamAtCluster2 = extrapTrackParamAtCluster2; } + } + } - // if no hitForRecCh1 found then consider to add hitForRecCh2 only - if (!foundSecondHit) { - if (fgkTrackAllTracks) { - // copy trackCandidate into a new track put at the end of fRecTracksPtr and add the new hitForRec's - newTrack = new ((*fRecTracksPtr)[fRecTracksPtr->GetLast()+1]) AliMUONTrack(*trackCandidate); + + // if no clusterCh1 found then consider to add clusterCh2 only + if (!foundSecondCluster) { + foundOneCluster = kTRUE; + + if (GetRecoParam()->TrackAllTracks()) { + // copy trackCandidate into a new track put at the end of fRecTracksPtr and add the new cluster + newTrack = new ((*fRecTracksPtr)[fRecTracksPtr->GetLast()+1]) AliMUONTrack(trackCandidate); + UpdateTrack(*newTrack,extrapTrackParamAtCluster2); fNRecTracks++; - AliMUONTrackParam trackParam1(extrapTrackParamSave); - AliMUONTrackExtrap::ExtrapToZ(&trackParam1, hitForRecCh2->GetZ()); - newTrack->AddTrackParamAtHit(&trackParam1,hitForRecCh2); - // Sort TrackParamAtHit according to increasing Z - newTrack->GetTrackParamAtHit()->Sort(); - // Update the chi2 of the new track - if (newTrack->GetFitFMin()<0) newTrack->SetFitFMin(chi2WithOneHitForRec); - else newTrack->SetFitFMin(newTrack->GetFitFMin() + chi2WithOneHitForRec); + // Printout for debuging if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) { - cout << "FollowTrackInStation: added one hit in chamber(1..): " << ch2+1 - << " (Chi2 = " << chi2WithOneHitForRec << ")" << endl; + cout << "FollowTrackInStation: added one cluster in chamber(1..): " << ch2+1 << endl; if (AliLog::GetGlobalDebugLevel() >= 3) newTrack->RecursiveDump(); } - } if (!bestHitForRec2 && chi2WithOneHitForRec < bestChi2WithOneHitForRec) { - // keep track of the best single hitForRec except if a couple - // of hits has already been found (i.e. bestHitForRec2!=0x0) - bestChi2WithOneHitForRec = chi2WithOneHitForRec; - bestHitForRec1 = hitForRecCh2; + + } else if (!foundTwoClusters && chi2WithOneCluster < bestChi2WithOneCluster) { + // keep track of the best single cluster except if a couple of clusters has already been found + bestChi2WithOneCluster = chi2WithOneCluster; + bestTrackParamAtCluster1 = extrapTrackParamAtCluster2; } + } + } - // reset the extrapolated track parameter for next step - trackCandidate->SetExtrapTrackParam(&extrapTrackParamSave); + } - // + // look for candidates in chamber 1 not already attached to a track - // if we want to keep all possible tracks or if no good couple of hitForRec has been found - if (fgkTrackAllTracks || !bestHitForRec2) { + // if we want to keep all possible tracks or if no good couple of clusters has been found + if (GetRecoParam()->TrackAllTracks() || !foundTwoClusters) { + + // add MCS effect for next step + AliMUONTrackExtrap::AddMCSEffect(&extrapTrackParamAtCh,AliMUONConstants::ChamberThicknessInX0(ch2),-1.); + + //Extrapolate trackCandidate to chamber "ch1" + Bool_t normalExtrap = AliMUONTrackExtrap::ExtrapToZCov(&extrapTrackParamAtCh, AliMUONConstants::DefaultChamberZ(ch1)); + + // Printout for debuging + if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 2) || (AliLog::GetGlobalDebugLevel() >= 2)) { + cout<= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) { - cout << "FollowTrackInStation: look for single hits in chamber(1..): " << ch1+1 << endl; + cout << "FollowTrackInStation: look for single clusters in chamber(1..): " << ch1+1 << endl; } - for (Int_t hit1 = 0; hit1 < fNHitsForRecPerChamber[ch1]; hit1++) { - hitForRecCh1 = (AliMUONHitForRec*) fHitsForRecPtr->UncheckedAt(fIndexOfFirstHitForRecPerChamber[ch1]+hit1); - if (hitForRecCh1Used[hit1]) continue; // Skip hitForRec already used - chi2WithOneHitForRec = trackCandidate->TryOneHitForRec(hitForRecCh1); - // if good chi2 then create a new track by adding the good hitForRec in "ch1" to the "trackCandidate" - // We do not try to attach a hitForRec in the other chamber too since it has already been done above - if (chi2WithOneHitForRec < maxChi2WithOneHitForRec) { - if (fgkTrackAllTracks) { - // copy trackCandidate into a new track put at the end of fRecTracksPtr and add the new hitForRec's - newTrack = new ((*fRecTracksPtr)[fRecTracksPtr->GetLast()+1]) AliMUONTrack(*trackCandidate); + + // reset cluster iterator of chamber 1 + nextInCh1.Reset(); + iCluster1 = -1; + + // look for second candidates in chamber 1 + if (normalExtrap) while ( ( clusterCh1 = static_cast(nextInCh1()) ) ) { + iCluster1++; + + if (clusterCh1Used[iCluster1]) continue; // Skip cluster already used + + // try to add the current cluster fast + if (!TryOneClusterFast(extrapTrackParamAtCh, clusterCh1)) continue; + + // try to add the current cluster accuratly + chi2WithOneCluster = TryOneCluster(extrapTrackParamAtCh, clusterCh1, extrapTrackParamAtCluster1); + + // if good chi2 then consider to add clusterCh1 + // We do not try to attach a cluster in the other chamber too since it has already been done above + if (chi2WithOneCluster < maxChi2WithOneCluster) { + foundOneCluster = kTRUE; + + // Printout for debuging + if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) { + cout << "FollowTrackInStation: found one cluster in chamber(1..): " << ch1+1 + << " (Chi2 = " << chi2WithOneCluster << ")" << endl; + clusterCh1->Print(); + } + + if (GetRecoParam()->TrackAllTracks()) { + // copy trackCandidate into a new track put at the end of fRecTracksPtr and add the new cluster + newTrack = new ((*fRecTracksPtr)[fRecTracksPtr->GetLast()+1]) AliMUONTrack(trackCandidate); + UpdateTrack(*newTrack,extrapTrackParamAtCluster1); fNRecTracks++; - AliMUONTrackParam trackParam1(extrapTrackParamSave); - AliMUONTrackExtrap::ExtrapToZ(&trackParam1, hitForRecCh1->GetZ()); - newTrack->AddTrackParamAtHit(&trackParam1,hitForRecCh1); - // Sort TrackParamAtHit according to increasing Z - newTrack->GetTrackParamAtHit()->Sort(); - // Update the chi2 of the new track - if (newTrack->GetFitFMin()<0) newTrack->SetFitFMin(chi2WithOneHitForRec); - else newTrack->SetFitFMin(newTrack->GetFitFMin() + chi2WithOneHitForRec); - // Printout for debuging + + // Printout for debuging if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) { - cout << "FollowTrackInStation: added one hit in chamber(1..): " << ch1+1 - << " (Chi2 = " << chi2WithOneHitForRec << ")" << endl; + cout << "FollowTrackInStation: added one cluster in chamber(1..): " << ch1+1 << endl; if (AliLog::GetGlobalDebugLevel() >= 3) newTrack->RecursiveDump(); } - } if (!bestHitForRec2 && chi2WithOneHitForRec < bestChi2WithOneHitForRec) { - // keep track of the best single hitForRec except if a couple - // of hits has already been found (i.e. bestHitForRec1!=0x0) - bestChi2WithOneHitForRec = chi2WithOneHitForRec; - bestHitForRec1 = hitForRecCh1; + + } else if (chi2WithOneCluster < bestChi2WithOneCluster) { + // keep track of the best single cluster except if a couple of clusters has already been found + bestChi2WithOneCluster = chi2WithOneCluster; + bestTrackParamAtCluster1 = extrapTrackParamAtCluster1; } + } + } + } - // + // fill out the best track if required else clean up the fRecTracksPtr array - if (!fgkTrackAllTracks && bestHitForRec1) { - AliMUONTrackParam trackParam1(extrapTrackParamSave); - AliMUONTrackExtrap::ExtrapToZ(&trackParam1, bestHitForRec1->GetZ()); - trackCandidate->AddTrackParamAtHit(&trackParam1,bestHitForRec1); - if (bestHitForRec2) { - AliMUONTrackParam trackParam2(extrapTrackParamSave); - AliMUONTrackExtrap::ExtrapToZ(&trackParam2, bestHitForRec2->GetZ()); - trackCandidate->AddTrackParamAtHit(&trackParam2,bestHitForRec2); - // Sort TrackParamAtHit according to increasing Z - trackCandidate->GetTrackParamAtHit()->Sort(); - // Update the chi2 of the new track - if (trackCandidate->GetFitFMin()<0) trackCandidate->SetFitFMin(bestChi2WithTwoHitForRec); - else trackCandidate->SetFitFMin(trackCandidate->GetFitFMin() + bestChi2WithTwoHitForRec); + if (!GetRecoParam()->TrackAllTracks()) { + if (foundTwoClusters) { + UpdateTrack(trackCandidate,bestTrackParamAtCluster1,bestTrackParamAtCluster2); + // Printout for debuging if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) { - cout << "FollowTrackInStation: added the two best hits in station(1..): " << nextStation+1 - << " (Chi2 = " << bestChi2WithTwoHitForRec << ")" << endl; + cout << "FollowTrackInStation: added the two best clusters in station(1..): " << nextStation+1 << endl; if (AliLog::GetGlobalDebugLevel() >= 3) newTrack->RecursiveDump(); } - } else { - // Sort TrackParamAtHit according to increasing Z - trackCandidate->GetTrackParamAtHit()->Sort(); - // Update the chi2 of the new track - if (trackCandidate->GetFitFMin()<0) trackCandidate->SetFitFMin(bestChi2WithOneHitForRec); - else trackCandidate->SetFitFMin(trackCandidate->GetFitFMin() + bestChi2WithOneHitForRec); + + } else if (foundOneCluster) { + UpdateTrack(trackCandidate,bestTrackParamAtCluster1); + // Printout for debuging if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) { - cout << "FollowTrackInStation: added the best hit in station(1..): " << nextStation+1 - << " (Chi2 = " << bestChi2WithOneHitForRec << ")" << endl; + cout << "FollowTrackInStation: added the best cluster in chamber(1..): " << bestTrackParamAtCluster1.GetClusterPtr()->GetChamberId()+1 << endl; if (AliLog::GetGlobalDebugLevel() >= 3) newTrack->RecursiveDump(); } + + } else { + delete [] clusterCh1Used; + return kFALSE; } - } else { - fRecTracksPtr->Remove(trackCandidate); // obsolete track + + } else if (foundOneCluster || foundTwoClusters) { + + // remove obsolete track + fRecTracksPtr->Remove(&trackCandidate); fNRecTracks--; + + } else { + delete [] clusterCh1Used; + return kFALSE; + } + + delete [] clusterCh1Used; + return kTRUE; + +} + + //__________________________________________________________________________ +Double_t AliMUONTrackReconstructor::TryTwoClusters(const AliMUONTrackParam &trackParamAtCluster1, AliMUONVCluster* cluster2, + AliMUONTrackParam &trackParamAtCluster2) +{ +/// Test the compatibility between the track and the 2 clusters together (using trackParam's covariance matrix): +/// return the corresponding Chi2 accounting for covariances between the 2 clusters +/// return trackParamAtCluster1 & 2 + + // extrapolate track parameters at the z position of the second cluster (no need to extrapolate the covariances) + // and set pointer to cluster into trackParamAtCluster + trackParamAtCluster2.SetParameters(trackParamAtCluster1.GetParameters()); + trackParamAtCluster2.SetZ(trackParamAtCluster1.GetZ()); + trackParamAtCluster2.SetClusterPtr(cluster2); + if (!AliMUONTrackExtrap::ExtrapToZ(&trackParamAtCluster2, cluster2->GetZ())) return 2.*AliMUONTrack::MaxChi2(); + + // Set differences between track and the 2 clusters in the bending and non bending directions + AliMUONVCluster* cluster1 = trackParamAtCluster1.GetClusterPtr(); + TMatrixD dPos(4,1); + dPos(0,0) = cluster1->GetX() - trackParamAtCluster1.GetNonBendingCoor(); + dPos(1,0) = cluster1->GetY() - trackParamAtCluster1.GetBendingCoor(); + dPos(2,0) = cluster2->GetX() - trackParamAtCluster2.GetNonBendingCoor(); + dPos(3,0) = cluster2->GetY() - trackParamAtCluster2.GetBendingCoor(); + + // Calculate the error matrix from the track parameter covariances at first cluster + TMatrixD error(4,4); + error.Zero(); + if (trackParamAtCluster1.CovariancesExist()) { + // Save track parameters at first cluster + TMatrixD paramAtCluster1Save(trackParamAtCluster1.GetParameters()); + + // Save track coordinates at second cluster + Double_t nonBendingCoor2 = trackParamAtCluster2.GetNonBendingCoor(); + Double_t bendingCoor2 = trackParamAtCluster2.GetBendingCoor(); + + // copy track parameters at first cluster for jacobian calculation + AliMUONTrackParam trackParam(trackParamAtCluster1); + + // Get the pointer to the parameter covariance matrix at first cluster + const TMatrixD& kParamCov = trackParam.GetCovariances(); + + // Calculate the jacobian related to the transformation between track parameters + // at first cluster and track coordinates at the 2 cluster z-positions + TMatrixD jacob(4,5); + jacob.Zero(); + // first derivative at the first cluster: + jacob(0,0) = 1.; // dx1/dx + 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; + // Small variation of parameter i only + 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.,direction[j]*paramAtCluster1Save(j,0)); // variation always in the same direction + } else dParam(j,0) = 0.; + } + + // Set new track parameters at first cluster + trackParam.SetParameters(paramAtCluster1Save); + trackParam.AddParameters(dParam); + trackParam.SetZ(cluster1->GetZ()); + + // Extrapolate new track parameters to the z position of the second cluster + if (!AliMUONTrackExtrap::ExtrapToZ(&trackParam,cluster2->GetZ())) return 2.*AliMUONTrack::MaxChi2(); + + // Calculate the jacobian + jacob(2,i) = (trackParam.GetNonBendingCoor() - nonBendingCoor2) / dParam(i,0); // dx2/dParami + jacob(3,i) = (trackParam.GetBendingCoor() - bendingCoor2 ) / dParam(i,0); // dy2/dParami + } + + // Calculate the error matrix + TMatrixD tmp(jacob,TMatrixD::kMult,kParamCov); + error = TMatrixD(tmp,TMatrixD::kMultTranspose,jacob); } + // Add cluster resolution to the error matrix + error(0,0) += cluster1->GetErrX2(); + error(1,1) += cluster1->GetErrY2(); + error(2,2) += cluster2->GetErrX2(); + error(3,3) += cluster2->GetErrY2(); + + // invert the error matrix for Chi2 calculation + if (error.Determinant() != 0) { + error.Invert(); + } else { + AliWarning(" Determinant error=0"); + return 2.*AliMUONTrack::MaxChi2(); + } + + // Compute the Chi2 value + TMatrixD tmp2(dPos,TMatrixD::kTransposeMult,error); + TMatrixD result(tmp2,TMatrixD::kMult,dPos); + + return result(0,0); + +} + + //__________________________________________________________________________ +void AliMUONTrackReconstructor::UpdateTrack(AliMUONTrack &track, AliMUONTrackParam &trackParamAtCluster) +{ + /// Add 1 cluster to the track candidate + /// Update chi2 of the track + + // Compute local chi2 + AliMUONVCluster* cluster = trackParamAtCluster.GetClusterPtr(); + Double_t deltaX = trackParamAtCluster.GetNonBendingCoor() - cluster->GetX(); + Double_t deltaY = trackParamAtCluster.GetBendingCoor() - cluster->GetY(); + Double_t localChi2 = deltaX*deltaX / cluster->GetErrX2() + + deltaY*deltaY / cluster->GetErrY2(); + + // Flag cluster as being not removable + if (GetRecoParam()->RequestStation(cluster->GetChamberId()/2)) + trackParamAtCluster.SetRemovable(kFALSE); + else trackParamAtCluster.SetRemovable(kTRUE); + trackParamAtCluster.SetLocalChi2(0.); // --> Local chi2 not used + + // Update the chi2 of the new track + track.SetGlobalChi2(track.GetGlobalChi2() + localChi2); + + // Update TrackParamAtCluster + track.AddTrackParamAtCluster(trackParamAtCluster,*cluster); + +} + + //__________________________________________________________________________ +void AliMUONTrackReconstructor::UpdateTrack(AliMUONTrack &track, AliMUONTrackParam &trackParamAtCluster1, AliMUONTrackParam &trackParamAtCluster2) +{ + /// Add 2 clusters to the track candidate + /// Update track and local chi2 + + // Update local chi2 at first cluster + AliMUONVCluster* cluster1 = trackParamAtCluster1.GetClusterPtr(); + Double_t deltaX = trackParamAtCluster1.GetNonBendingCoor() - cluster1->GetX(); + Double_t deltaY = trackParamAtCluster1.GetBendingCoor() - cluster1->GetY(); + Double_t localChi2AtCluster1 = deltaX*deltaX / cluster1->GetErrX2() + + deltaY*deltaY / cluster1->GetErrY2(); + trackParamAtCluster1.SetLocalChi2(localChi2AtCluster1); + + // Flag first cluster as being removable + trackParamAtCluster1.SetRemovable(kTRUE); + + // Update local chi2 at second cluster + AliMUONVCluster* cluster2 = trackParamAtCluster2.GetClusterPtr(); + deltaX = trackParamAtCluster2.GetNonBendingCoor() - cluster2->GetX(); + deltaY = trackParamAtCluster2.GetBendingCoor() - cluster2->GetY(); + Double_t localChi2AtCluster2 = deltaX*deltaX / cluster2->GetErrX2() + + deltaY*deltaY / cluster2->GetErrY2(); + trackParamAtCluster2.SetLocalChi2(localChi2AtCluster2); + + // Flag first cluster as being removable + trackParamAtCluster2.SetRemovable(kTRUE); + + // Update the chi2 of the new track + track.SetGlobalChi2(track.GetGlobalChi2() + localChi2AtCluster1 + localChi2AtCluster2); + + // Update TrackParamAtCluster + track.AddTrackParamAtCluster(trackParamAtCluster1,*cluster1); + track.AddTrackParamAtCluster(trackParamAtCluster2,*cluster2); + } //__________________________________________________________________________ -void AliMUONTrackReconstructor::SetVertexForFit(AliMUONTrack* trackCandidate) +Bool_t AliMUONTrackReconstructor::RecoverTrack(AliMUONTrack &trackCandidate, AliMUONVClusterStore& clusterStore, Int_t nextStation) { - /// Add the vertex as a measured hit to constrain the fit of the "trackCandidate" - /// Compute the vertex resolution from natural vertex dispersion and + /// Try to recover the track candidate in the next station + /// by removing the worst of the two clusters attached in the current station + /// Return kTRUE if recovering succeeds + AliDebug(1,"Enter RecoverTrack"); + + // Do not try to recover track until we have attached cluster(s) on station(1..) 3 + if (nextStation > 1) return kFALSE; + + Int_t worstClusterNumber = -1; + Double_t localChi2, worstLocalChi2 = -1.; + + // Look for the cluster to remove + for (Int_t clusterNumber = 0; clusterNumber < 2; clusterNumber++) { + AliMUONTrackParam *trackParamAtCluster = (AliMUONTrackParam*)trackCandidate.GetTrackParamAtCluster()->UncheckedAt(clusterNumber); + + // check if current cluster is in the previous station + if (trackParamAtCluster->GetClusterPtr()->GetChamberId()/2 != nextStation+1) break; + + // check if current cluster is removable + if (!trackParamAtCluster->IsRemovable()) return kFALSE; + + // reset the current cluster as beig not removable if it is on a required station + if (GetRecoParam()->RequestStation(nextStation+1)) trackParamAtCluster->SetRemovable(kFALSE); + + // Pick up cluster with the worst chi2 + localChi2 = trackParamAtCluster->GetLocalChi2(); + if (localChi2 > worstLocalChi2) { + worstLocalChi2 = localChi2; + worstClusterNumber = clusterNumber; + } + + } + + // check if worst cluster found + if (worstClusterNumber < 0) return kFALSE; + + // Remove the worst cluster + trackCandidate.RemoveTrackParamAtCluster((AliMUONTrackParam*)trackCandidate.GetTrackParamAtCluster()->UncheckedAt(worstClusterNumber)); + + // Re-fit the track: + // Do not take into account the multiple scattering to speed up the fit + // Calculate the track parameter covariance matrix + Fit(trackCandidate, kFALSE, kFALSE, kTRUE); + + // skip track if the normalized chi2 is too high + if (trackCandidate.GetNormalizedChi2() > GetRecoParam()->GetSigmaCutForTracking() * GetRecoParam()->GetSigmaCutForTracking()) 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); + +} + + //__________________________________________________________________________ +void AliMUONTrackReconstructor::SetVertexErrXY2ForFit(AliMUONTrack &trackCandidate) +{ + /// Compute the vertex resolution square from natural vertex dispersion and /// multiple scattering effets according to trackCandidate path in absorber /// It is necessary to account for multiple scattering effects here instead of during the fit of /// the "trackCandidate" to do not influence the result by changing track resolution at vertex AliDebug(1,"Enter SetVertexForFit"); - Double_t nonBendingReso2 = fNonBendingVertexDispersion * fNonBendingVertexDispersion; - Double_t bendingReso2 = fBendingVertexDispersion * fBendingVertexDispersion; + Double_t nonBendingReso2 = GetRecoParam()->GetNonBendingVertexDispersion() * + GetRecoParam()->GetNonBendingVertexDispersion(); + Double_t bendingReso2 = GetRecoParam()->GetBendingVertexDispersion() * + GetRecoParam()->GetBendingVertexDispersion(); + // add multiple scattering effets - AliMUONTrackParam paramAtVertex(*((AliMUONTrackParam*)(trackCandidate->GetTrackParamAtHit()->First()))); + AliMUONTrackParam paramAtVertex(*((AliMUONTrackParam*)(trackCandidate.GetTrackParamAtCluster()->First()))); paramAtVertex.DeleteCovariances(); // to be sure to account only for multiple scattering - AliMUONTrackExtrap::ExtrapToVertexUncorrected(¶mAtVertex,0.); - TMatrixD* paramCov = paramAtVertex.GetCovariances(); - nonBendingReso2 += (*paramCov)(0,0); - bendingReso2 += (*paramCov)(2,2); - // Set the vertex - AliMUONHitForRec vertex; // Coordinates set to (0.,0.,0.) by default - vertex.SetNonBendingReso2(nonBendingReso2); - vertex.SetBendingReso2(bendingReso2); - trackCandidate->SetVertex(&vertex); + if (!AliMUONTrackExtrap::ExtrapToZ(¶mAtVertex,AliMUONConstants::AbsZEnd())) { + nonBendingReso2 = 0.; + bendingReso2 = 0.; + } else { + AliMUONTrackExtrap::ExtrapToVertexUncorrected(¶mAtVertex,0.); + const TMatrixD& kParamCov = paramAtVertex.GetCovariances(); + nonBendingReso2 += kParamCov(0,0); + bendingReso2 += kParamCov(2,2); + } + + // Set the vertex resolution square + trackCandidate.SetVertexErrXY2(nonBendingReso2,bendingReso2); } //__________________________________________________________________________ -void AliMUONTrackReconstructor::Fit(AliMUONTrack *track, Bool_t includeMCS, Bool_t calcCov) +void AliMUONTrackReconstructor::Fit(AliMUONTrack &track, Bool_t includeMCS, Bool_t fitWithVertex, Bool_t calcCov) { - /// Fit the track "track" with or without multiple Coulomb scattering according to "includeMCS". + /// Fit the track + /// w/wo multiple Coulomb scattering according to "includeMCS". + /// w/wo constraining the vertex according to "fitWithVertex". + /// calculating or not the covariance matrix according to "calcCov". + AliDebug(1,"Enter Fit"); Double_t benC, errorParam, invBenP, nonBenC, x, y; AliMUONTrackParam *trackParam; - Double_t arg[1], fedm, errdef, fitFMin; + Double_t arg[1], fedm, errdef, globalChi2; Int_t npari, nparx; Int_t status, covStatus; @@ -636,7 +1184,7 @@ void AliMUONTrackReconstructor::Fit(AliMUONTrack *track, Bool_t includeMCS, Bool // Clear MINUIT parameters gMinuit->mncler(); // Give the fitted track to MINUIT - gMinuit->SetObjectFit(track); + gMinuit->SetObjectFit(&track); if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 2) || (AliLog::GetGlobalDebugLevel() >= 2)) { // Define print level arg[0] = 1; @@ -653,20 +1201,31 @@ void AliMUONTrackReconstructor::Fit(AliMUONTrack *track, Bool_t includeMCS, Bool //arg[0] = 2; //gMinuit->mnexcm("SET STR", arg, 1, status); - // Switch between available FCN according to "includeMCS" - if (includeMCS) gMinuit->SetFCN(TrackChi2MCS); - else gMinuit->SetFCN(TrackChi2); + // set flag w/wo multiple scattering according to "includeMCS" + track.FitWithMCS(includeMCS); + if (includeMCS) { + // compute cluster weights only once + if (!track.UpdateTrackParamAtCluster() || !track.ComputeClusterWeights()) { + AliWarning("cannot take into account the multiple scattering effects"); + track.FitWithMCS(kFALSE); + } + } + + track.FitWithVertex(fitWithVertex); + if (fitWithVertex) SetVertexErrXY2ForFit(track); + + // Set fitting function + gMinuit->SetFCN(TrackChi2); // Set fitted parameters (!! The order is very important for the covariance matrix !!) - trackParam = (AliMUONTrackParam*) (track->GetTrackParamAtHit()->First()); - // could be tried with no limits for the search (min=max=0) ???? - // mandatory limits in non Bending to avoid NaN values of parameters + // Mandatory limits to avoid NaN values of parameters + trackParam = (AliMUONTrackParam*) (track.GetTrackParamAtCluster()->First()); + Double_t maxIBM = 1. / GetRecoParam()->GetMinBendingMomentum(); gMinuit->mnparm(0, "X", trackParam->GetNonBendingCoor(), 0.03, -500.0, 500.0, status); - gMinuit->mnparm(1, "NonBenS", trackParam->GetNonBendingSlope(), 0.001, -0.5, 0.5, status); - // mandatory limits in Bending to avoid NaN values of parameters + gMinuit->mnparm(1, "NonBenS", trackParam->GetNonBendingSlope(), 0.001, -1., 1., status); gMinuit->mnparm(2, "Y", trackParam->GetBendingCoor(), 0.10, -500.0, 500.0, status); - gMinuit->mnparm(3, "BenS", trackParam->GetBendingSlope(), 0.001, -0.5, 0.5, status); - gMinuit->mnparm(4, "InvBenP", trackParam->GetInverseBendingMomentum(), 0.003, -0.4, 0.4, status); + gMinuit->mnparm(3, "BenS", trackParam->GetBendingSlope(), 0.001, -1.5, 1.5, status); + gMinuit->mnparm(4, "InvBenP", trackParam->GetInverseBendingMomentum(), 0.003, -maxIBM, maxIBM, status); // minimization gMinuit->mnexcm("MIGRAD", arg, 0, status); @@ -687,347 +1246,286 @@ void AliMUONTrackReconstructor::Fit(AliMUONTrack *track, Bool_t includeMCS, Bool trackParam->SetInverseBendingMomentum(invBenP); // global result of the fit - gMinuit->mnstat(fitFMin, fedm, errdef, npari, nparx, covStatus); - track->SetFitFMin(fitFMin); + gMinuit->mnstat(globalChi2, fedm, errdef, npari, nparx, covStatus); + track.SetGlobalChi2(globalChi2); // Get the covariance matrix if required if (calcCov) { // Covariance matrix according to HESSE status // If problem then keep only the diagonal terms (variances) Double_t matrix[5][5]; + for (Int_t i=0; i<5; i++) for (Int_t j=0; j<5; j++) matrix[i][j] = 0.; gMinuit->mnemat(&matrix[0][0],5); if (covStatus == 3) trackParam->SetCovariances(matrix); else trackParam->SetVariances(matrix); - } else *(trackParam->GetCovariances()) = 0; + } else trackParam->DeleteCovariances(); } //__________________________________________________________________________ -void TrackChi2(Int_t & /*NParam*/, Double_t * /*Gradient*/, Double_t &Chi2, Double_t *Param, Int_t /*Flag*/) +void TrackChi2(Int_t & /*nParam*/, Double_t * /*gradient*/, Double_t &chi2, Double_t *param, Int_t /*flag*/) { - /// Return the "Chi2" to be minimized with Minuit for track fitting, - /// with "NParam" parameters - /// and their current values in array pointed to by "Param". - /// 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 + /// Return the "Chi2" to be minimized with Minuit for track fitting. + /// Assumes that the trackParamAtCluster are sorted according to increasing Z. + /// Track parameters at each cluster are updated accordingly. + /// Vertex is used according to the flag "trackBeingFitted->GetFitWithVertex()". + /// Multiple Coulomb scattering is taken into account according to the flag "trackBeingFitted->GetFitWithMCS()". AliMUONTrack *trackBeingFitted = (AliMUONTrack*) gMinuit->GetObjectFit(); -// AliMUONTrack *trackBeingFitted = (AliMUONTrack*) AliMUONTrackReconstructor::Fitter()->GetObjectFit(); - AliMUONTrackParam param1; - AliMUONTrackParam* trackParamAtHit; - AliMUONHitForRec* hitForRec; - Chi2 = 0.0; // initialize Chi2 + AliMUONTrackParam* trackParamAtCluster = (AliMUONTrackParam*) trackBeingFitted->GetTrackParamAtCluster()->First(); Double_t dX, dY; + chi2 = 0.; // initialize chi2 - // copy of track parameters to be fitted - param1 = *((AliMUONTrackParam*) (trackBeingFitted->GetTrackParamAtHit()->First())); - param1.SetNonBendingCoor(Param[0]); - param1.SetNonBendingSlope(Param[1]); - param1.SetBendingCoor(Param[2]); - param1.SetBendingSlope(Param[3]); - param1.SetInverseBendingMomentum(Param[4]); + // update track parameters + trackParamAtCluster->SetNonBendingCoor(param[0]); + trackParamAtCluster->SetNonBendingSlope(param[1]); + trackParamAtCluster->SetBendingCoor(param[2]); + trackParamAtCluster->SetBendingSlope(param[3]); + trackParamAtCluster->SetInverseBendingMomentum(param[4]); + if (!trackBeingFitted->UpdateTrackParamAtCluster()) { + chi2 = 2.*AliMUONTrack::MaxChi2(); + return; + } // Take the vertex into account in the fit if required - if (trackBeingFitted->GetFitWithVertex()) { - AliMUONTrackParam paramAtVertex(param1); - AliMUONTrackExtrap::ExtrapToZ(¶mAtVertex, 0.); - AliMUONHitForRec *vertex = trackBeingFitted->GetVertex(); - if (!vertex) { - cout<<"Error in TrackChi2: Want to use the vertex in tracking but it has not been created!!"<FitWithVertex()) { + Double_t nonBendingReso2,bendingReso2; + trackBeingFitted->GetVertexErrXY2(nonBendingReso2,bendingReso2); + AliMUONTrackParam paramAtVertex(*trackParamAtCluster); + if (nonBendingReso2 != 0. && bendingReso2 != 0. && AliMUONTrackExtrap::ExtrapToZ(¶mAtVertex, 0.)) { // vextex position = (0,0,0) + dX = paramAtVertex.GetNonBendingCoor(); + dY = paramAtVertex.GetBendingCoor(); + chi2 += dX * dX / nonBendingReso2 + dY * dY / bendingReso2; + } else { + chi2 = 2.*AliMUONTrack::MaxChi2(); + return; } - dX = vertex->GetNonBendingCoor() - paramAtVertex.GetNonBendingCoor(); - dY = vertex->GetBendingCoor() - paramAtVertex.GetBendingCoor(); - Chi2 += dX * dX / vertex->GetNonBendingReso2() + dY * dY / vertex->GetBendingReso2(); } - // Follow track through all planes of track hits - trackParamAtHit = (AliMUONTrackParam*) (trackBeingFitted->GetTrackParamAtHit()->First()); - while (trackParamAtHit) { - hitForRec = trackParamAtHit->GetHitForRecPtr(); - // extrapolation to the plane of the hitForRec attached to the current trackParamAtHit - AliMUONTrackExtrap::ExtrapToZ(¶m1, hitForRec->GetZ()); - // update track parameters of the current hit - trackParamAtHit->SetTrackParam(param1); - // Increment Chi2 - // done hit per hit, with hit resolution, - // and not with point and angle like in "reco_muon.F" !!!! - dX = hitForRec->GetNonBendingCoor() - param1.GetNonBendingCoor(); - dY = hitForRec->GetBendingCoor() - param1.GetBendingCoor(); - Chi2 = Chi2 + dX * dX / hitForRec->GetNonBendingReso2() + dY * dY / hitForRec->GetBendingReso2(); - trackParamAtHit = (AliMUONTrackParam*) (trackBeingFitted->GetTrackParamAtHit()->After(trackParamAtHit)); - } + // compute chi2 w/wo multiple scattering + chi2 += trackBeingFitted->ComputeGlobalChi2(trackBeingFitted->FitWithMCS()); + } //__________________________________________________________________________ -void TrackChi2MCS(Int_t & /*NParam*/, Double_t * /*Gradient*/, Double_t &Chi2, Double_t *Param, Int_t /*Flag*/) +Bool_t AliMUONTrackReconstructor::ComplementTracks(const AliMUONVClusterStore& clusterStore) { - /// Return the "Chi2" to be minimized with Minuit for track fitting, - /// with "NParam" parameters - /// and their current values in array pointed to by "Param". - /// 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. + /// Complete tracks by adding missing clusters (if there is an overlap between + /// two detection elements, the track may have two clusters in the same chamber). + /// Re-fit track parameters and covariances. + /// Return kTRUE if one or more tracks have been complemented. + AliDebug(1,"Enter ComplementTracks"); - AliMUONTrack *trackBeingFitted = (AliMUONTrack*) gMinuit->GetObjectFit(); -// AliMUONTrack *trackBeingFitted = (AliMUONTrack*) AliMUONTrackReconstructor::Fitter()->GetObjectFit(); - AliMUONTrackParam param1; - AliMUONTrackParam* trackParamAtHit; - AliMUONHitForRec* hitForRec; - Chi2 = 0.0; // initialize Chi2 - Int_t chCurrent, chPrev = 0, hitNumber, hitNumber1, hitNumber2, hitNumber3; - Double_t z1, z2, z3; - AliMUONTrackParam *trackParamAtHit1, *trackParamAtHit2, *trackParamAtHit3; - AliMUONHitForRec *hitForRec1, *hitForRec2; - Double_t hbc1, hbc2, pbc1, pbc2; - Double_t hnbc1, hnbc2, pnbc1, pnbc2; - Int_t numberOfHit = trackBeingFitted->GetNTrackHits(); - TMatrixD *covBending = new TMatrixD(numberOfHit, numberOfHit); - TMatrixD *covNonBending = new TMatrixD(numberOfHit, numberOfHit); - Double_t *msa2 = new Double_t[numberOfHit]; - - // copy of track parameters to be fitted - param1 = *((AliMUONTrackParam*) (trackBeingFitted->GetTrackParamAtHit()->First())); - param1.SetNonBendingCoor(Param[0]); - param1.SetNonBendingSlope(Param[1]); - param1.SetBendingCoor(Param[2]); - param1.SetBendingSlope(Param[3]); - param1.SetInverseBendingMomentum(Param[4]); - - // Take the vertex into account in the fit if required - if (trackBeingFitted->GetFitWithVertex()) { - AliMUONTrackParam paramAtVertex(param1); - AliMUONTrackExtrap::ExtrapToZ(¶mAtVertex, 0.); - AliMUONHitForRec *vertex = trackBeingFitted->GetVertex(); - if (!vertex) { - cout<<"Error in TrackChi2MCS: Want to use the vertex in tracking but it has not been created!!"<GetNonBendingCoor() - paramAtVertex.GetNonBendingCoor(); - Double_t dY = vertex->GetBendingCoor() - paramAtVertex.GetBendingCoor(); - Chi2 += dX * dX / vertex->GetNonBendingReso2() + dY * dY / vertex->GetBendingReso2(); - } + Int_t chamberId, detElemId; + Double_t chi2OfCluster, bestChi2OfCluster; + Double_t sigmaCut2 = GetRecoParam()->GetSigmaCutForTracking() * + GetRecoParam()->GetSigmaCutForTracking(); + Bool_t foundOneCluster, trackModified, hasChanged = kFALSE; + AliMUONVCluster* cluster; + AliMUONTrackParam *trackParam, *nextTrackParam, copyOfTrackParam, trackParamAtCluster, bestTrackParamAtCluster; - // Predicted coordinates and multiple scattering angles are first calculated - for (hitNumber = 0; hitNumber < numberOfHit; hitNumber++) { - trackParamAtHit = (AliMUONTrackParam*) (trackBeingFitted->GetTrackParamAtHit()->UncheckedAt(hitNumber)); - hitForRec = trackParamAtHit->GetHitForRecPtr(); - // extrapolation to the plane of the hitForRec attached to the current trackParamAtHit - AliMUONTrackExtrap::ExtrapToZ(¶m1, hitForRec->GetZ()); - // update track parameters of the current hit - trackParamAtHit->SetTrackParam(param1); - // square of multiple scattering angle at current hit, with one chamber - msa2[hitNumber] = MultipleScatteringAngle2(¶m1); - // correction for eventual missing hits or multiple hits in a chamber, - // according to the number of chambers - // between the current hit and the previous one - chCurrent = hitForRec->GetChamberNumber(); - if (hitNumber > 0) msa2[hitNumber] = msa2[hitNumber] * (chCurrent - chPrev); - chPrev = chCurrent; - } - - // Calculates the covariance matrix - for (hitNumber1 = 0; hitNumber1 < numberOfHit; hitNumber1++) { - trackParamAtHit1 = (AliMUONTrackParam*) (trackBeingFitted->GetTrackParamAtHit()->UncheckedAt(hitNumber1)); - hitForRec1 = trackParamAtHit1->GetHitForRecPtr(); - z1 = hitForRec1->GetZ(); - for (hitNumber2 = hitNumber1; hitNumber2 < numberOfHit; hitNumber2++) { - trackParamAtHit2 = (AliMUONTrackParam*) (trackBeingFitted->GetTrackParamAtHit()->UncheckedAt(hitNumber2)); - z2 = trackParamAtHit2->GetHitForRecPtr()->GetZ(); - // 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++) { - trackParamAtHit3 = (AliMUONTrackParam*) (trackBeingFitted->GetTrackParamAtHit()->UncheckedAt(hitNumber3)); - z3 = trackParamAtHit3->GetHitForRecPtr()->GetZ(); - (*covBending)(hitNumber2, hitNumber1) = (*covBending)(hitNumber2, hitNumber1) + ((z1 - z3) * (z2 - z3) * msa2[hitNumber3]); + AliMUONTrack *track = (AliMUONTrack*) fRecTracksPtr->First(); + while (track) { + trackModified = kFALSE; + + trackParam = (AliMUONTrackParam*)track->GetTrackParamAtCluster()->First(); + while (trackParam) { + foundOneCluster = kFALSE; + bestChi2OfCluster = 2. * sigmaCut2; // 2 because 2 quantities in chi2 + chamberId = trackParam->GetClusterPtr()->GetChamberId(); + detElemId = trackParam->GetClusterPtr()->GetDetElemId(); + + // prepare nextTrackParam before adding new cluster because of the sorting + nextTrackParam = (AliMUONTrackParam*)track->GetTrackParamAtCluster()->After(trackParam); + + // recover track parameters from local fit and put them into a copy of trackParam + copyOfTrackParam.SetZ(trackParam->GetZ()); + copyOfTrackParam.SetParameters(trackParam->GetSmoothParameters()); + copyOfTrackParam.SetCovariances(trackParam->GetSmoothCovariances()); + + // Create iterators to loop over clusters in current chamber + TIter nextInCh(clusterStore.CreateChamberIterator(chamberId,chamberId)); + + // look for one second candidate in the same chamber + while ( ( cluster = static_cast(nextInCh()) ) ) { + + // look for a cluster in another detection element + if (cluster->GetDetElemId() == detElemId) continue; + + // try to add the current cluster fast + if (!TryOneClusterFast(copyOfTrackParam, cluster)) continue; + + // try to add the current cluster accurately + chi2OfCluster = TryOneCluster(copyOfTrackParam, cluster, trackParamAtCluster); + + // if better chi2 then prepare to add this cluster to the track + if (chi2OfCluster < bestChi2OfCluster) { + bestChi2OfCluster = chi2OfCluster; + bestTrackParamAtCluster = trackParamAtCluster; + foundOneCluster = kTRUE; + } + } - // equal contribution from multiple scattering in non bending plane - (*covNonBending)(hitNumber2, hitNumber1) = (*covBending)(hitNumber2, hitNumber1); - if (hitNumber1 == hitNumber2) { - // Diagonal elements: add contribution from position measurements - // in bending plane - (*covBending)(hitNumber2, hitNumber1) = (*covBending)(hitNumber2, hitNumber1) + hitForRec1->GetBendingReso2(); - // and in non bending plane - (*covNonBending)(hitNumber2, hitNumber1) = (*covNonBending)(hitNumber2, hitNumber1) + hitForRec1->GetNonBendingReso2(); - } else { - // Non diagonal elements: symmetrization - // for bending plane - (*covBending)(hitNumber1, hitNumber2) = (*covBending)(hitNumber2, hitNumber1); - // and non bending plane - (*covNonBending)(hitNumber1, hitNumber2) = (*covNonBending)(hitNumber2, hitNumber1); + + // add new cluster if any + if (foundOneCluster) { + + // Printout for debuging + if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) { + cout << "ComplementTracks: found one cluster in chamber(1..): " << chamberId+1 << endl; + bestTrackParamAtCluster.GetClusterPtr()->Print(); + cout<SetRemovable(kTRUE); + bestTrackParamAtCluster.SetRemovable(kTRUE); + track->AddTrackParamAtCluster(bestTrackParamAtCluster,*(bestTrackParamAtCluster.GetClusterPtr())); + trackModified = kTRUE; + hasChanged = kTRUE; } - } // for (hitNumber2 = hitNumber1;... - } // for (hitNumber1 = 0;... + + trackParam = nextTrackParam; + } - // Inversion of covariance matrices - Int_t ifailBending; - gMinuit->mnvert(&((*covBending)(0,0)), numberOfHit, numberOfHit, numberOfHit, ifailBending); - Int_t ifailNonBending; - gMinuit->mnvert(&((*covNonBending)(0,0)), numberOfHit, numberOfHit, numberOfHit, ifailNonBending); + // re-fit track parameters if needed + if (trackModified) Fit(*track, kTRUE, kFALSE, kTRUE); + + track = (AliMUONTrack*) fRecTracksPtr->After(track); + } + + return hasChanged; + +} - // 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 ((ifailBending == 0) && (ifailNonBending == 0)) { - // with Multiple Scattering if inversion correct - for (hitNumber1 = 0; hitNumber1 < numberOfHit ; hitNumber1++) { - trackParamAtHit1 = (AliMUONTrackParam*) (trackBeingFitted->GetTrackParamAtHit()->UncheckedAt(hitNumber1)); - hitForRec1 = trackParamAtHit1->GetHitForRecPtr(); - hbc1 = hitForRec1->GetBendingCoor(); - pbc1 = trackParamAtHit1->GetBendingCoor(); - hnbc1 = hitForRec1->GetNonBendingCoor(); - pnbc1 = trackParamAtHit1->GetNonBendingCoor(); - for (hitNumber2 = 0; hitNumber2 < numberOfHit; hitNumber2++) { - trackParamAtHit2 = (AliMUONTrackParam*) (trackBeingFitted->GetTrackParamAtHit()->UncheckedAt(hitNumber2)); - hitForRec2 = trackParamAtHit2->GetHitForRecPtr(); - hbc2 = hitForRec2->GetBendingCoor(); - pbc2 = trackParamAtHit2->GetBendingCoor(); - hnbc2 = hitForRec2->GetNonBendingCoor(); - pnbc2 = trackParamAtHit2->GetNonBendingCoor(); - Chi2 += ((*covBending)(hitNumber2, hitNumber1) * (hbc1 - pbc1) * (hbc2 - pbc2)) + - ((*covNonBending)(hitNumber2, hitNumber1) * (hnbc1 - pnbc1) * (hnbc2 - pnbc2)); + //__________________________________________________________________________ +void AliMUONTrackReconstructor::ImproveTrack(AliMUONTrack &track) +{ + /// Improve the given track by removing clusters with local chi2 highter than the defined cut + /// Recompute track parameters and covariances at the remaining clusters + AliDebug(1,"Enter ImproveTrack"); + + Double_t localChi2, worstLocalChi2; + AliMUONTrackParam *trackParamAtCluster, *worstTrackParamAtCluster; + Double_t sigmaCut2 = GetRecoParam()->GetSigmaCutForImprovement() * + GetRecoParam()->GetSigmaCutForImprovement(); + + while (!track.IsImproved()) { + + // identify removable clusters + track.TagRemovableClusters(GetRecoParam()->RequestedStationMask()); + + // Update track parameters and covariances + if (!track.UpdateCovTrackParamAtCluster()) { + AliWarning("unable to update track parameters and covariances --> stop improvement"); + break; + } + + // Compute local chi2 of each clusters + track.ComputeLocalChi2(kTRUE); + + // Look for the cluster to remove + worstTrackParamAtCluster = NULL; + worstLocalChi2 = 0.; + trackParamAtCluster = (AliMUONTrackParam*) track.GetTrackParamAtCluster()->First(); + while (trackParamAtCluster) { + + // Pick up cluster with the worst chi2 + localChi2 = trackParamAtCluster->GetLocalChi2(); + if (localChi2 > worstLocalChi2) { + worstLocalChi2 = localChi2; + worstTrackParamAtCluster = trackParamAtCluster; } + + trackParamAtCluster = (AliMUONTrackParam*) track.GetTrackParamAtCluster()->After(trackParamAtCluster); } - } else { - // without Multiple Scattering if inversion impossible - for (hitNumber1 = 0; hitNumber1 < numberOfHit ; hitNumber1++) { - trackParamAtHit1 = (AliMUONTrackParam*) (trackBeingFitted->GetTrackParamAtHit()->UncheckedAt(hitNumber1)); - hitForRec1 = trackParamAtHit1->GetHitForRecPtr(); - hbc1 = hitForRec1->GetBendingCoor(); - pbc1 = trackParamAtHit1->GetBendingCoor(); - hnbc1 = hitForRec1->GetNonBendingCoor(); - pnbc1 = trackParamAtHit1->GetNonBendingCoor(); - Chi2 += ((hbc1 - pbc1) * (hbc1 - pbc1) / hitForRec1->GetBendingReso2()) + - ((hnbc1 - pnbc1) * (hnbc1 - pnbc1) / hitForRec1->GetNonBendingReso2()); + + // Check if worst cluster found + if (!worstTrackParamAtCluster) { + AliWarning("Bad local chi2 values?"); + break; } + + // Check whether the worst chi2 is under requirement or not + if (worstLocalChi2 < 2. * sigmaCut2) { // 2 because 2 quantities in chi2 + track.SetImproved(kTRUE); + break; + } + + // if the worst cluster is not removable then stop improvement + if (!worstTrackParamAtCluster->IsRemovable()) break; + + // Remove the worst cluster + track.RemoveTrackParamAtCluster(worstTrackParamAtCluster); + + // Re-fit the track: + // Take into account the multiple scattering + // Calculate the track parameter covariance matrix + Fit(track, kTRUE, kFALSE, kTRUE); + + // Printout for debuging + if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) { + cout << "ImproveTracks: track " << fRecTracksPtr->IndexOf(&track)+1 << " improved " << endl; + } + } - delete covBending; - delete covNonBending; - delete [] msa2; } //__________________________________________________________________________ -Double_t MultipleScatteringAngle2(AliMUONTrackParam *param) +Bool_t AliMUONTrackReconstructor::FinalizeTrack(AliMUONTrack &track) { - /// Returns square of multiple Coulomb scattering angle - /// from TrackParamAtHit pointed to by "param" - Double_t slopeBending, slopeNonBending, radiationLength, inverseBendingMomentum2, inverseTotalMomentum2; - Double_t varMultipleScatteringAngle; - slopeBending = param->GetBendingSlope(); - slopeNonBending = param->GetNonBendingSlope(); - // thickness in radiation length for the current track, - // taking local angle into account - radiationLength = AliMUONConstants::ChamberThicknessInX0() * - TMath::Sqrt(1.0 + slopeBending*slopeBending + slopeNonBending*slopeNonBending); - inverseBendingMomentum2 = param->GetInverseBendingMomentum() * param->GetInverseBendingMomentum(); - inverseTotalMomentum2 = 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; + /// Recompute track parameters and covariances at each attached cluster + /// from those at the first one, if not already done + AliDebug(1,"Enter FinalizeTrack"); + if (!track.IsImproved() && !track.UpdateCovTrackParamAtCluster()) { + AliWarning("finalization failed due to extrapolation problem"); + return kFALSE; + } + return kTRUE; } //__________________________________________________________________________ -void AliMUONTrackReconstructor::FillMUONTrack(void) +Bool_t AliMUONTrackReconstructor::RefitTrack(AliMUONTrack &track, Bool_t enableImprovement) { - /// Fill AliMUONTrack's fHitForRecAtHit array - /// Recompute track parameters and covariances at each attached cluster from those at the first one - AliMUONTrack *track; - AliMUONTrackParam trackParam; - AliMUONTrackParam *trackParamAtHit; - AliMUONHitForRec *hitForRecAtHit; + /// re-fit the given track - track = (AliMUONTrack*) fRecTracksPtr->First(); - while (track) { - trackParamAtHit = (AliMUONTrackParam*) (track->GetTrackParamAtHit()->First()); - trackParam = *trackParamAtHit; - while (trackParamAtHit) { - hitForRecAtHit = trackParamAtHit->GetHitForRecPtr(); - // extrapolation to the plane of the hitForRec attached to the current trackParamAtHit - AliMUONTrackExtrap::ExtrapToZCov(&trackParam, hitForRecAtHit->GetZ()); - // update track parameters of the current hit - trackParamAtHit->SetTrackParam(trackParam); - // update covariance matrix of track parameters of the current hit - trackParamAtHit->SetCovariances(trackParam.GetCovariances()); - // update array of track hit - track->AddHitForRecAtHit(hitForRecAtHit); - // prepare next step, add MCS effects in parameter covariances - AliMUONTrackExtrap::AddMCSEffect(&trackParam,AliMUONConstants::ChamberThicknessInX0(),1.); - trackParamAtHit = (AliMUONTrackParam*) (track->GetTrackParamAtHit()->After(trackParamAtHit)); - } - track = (AliMUONTrack*) fRecTracksPtr->After(track); + // check validity of the track + if (track.GetNClusters() < 3) { + AliWarning("the track does not contain enough clusters --> unable to refit"); + return kFALSE; } - return; -} - - //__________________________________________________________________________ -void AliMUONTrackReconstructor::EventDump(void) -{ - /// Dump reconstructed event (track parameters at vertex and at first hit), - /// and the particle parameters - AliMUONTrack *track; - AliMUONTrackParam trackParam, *trackParam1; - Double_t bendingSlope, nonBendingSlope, pYZ; - Double_t pX, pY, pZ, x, y, z, c; - Int_t trackIndex, nTrackHits; - - AliDebug(1,"****** enter EventDump ******"); - AliDebug(1, Form("Number of Reconstructed tracks : %d", fNRecTracks)); - - fRecTracksPtr->Compress(); // for simple loop without "Next" since no hole - // Loop over reconstructed tracks - for (trackIndex = 0; trackIndex < fNRecTracks; trackIndex++) { - AliDebug(1, Form("track number: %d", trackIndex)); - // function for each track for modularity ???? - track = (AliMUONTrack*) ((*fRecTracksPtr)[trackIndex]); - nTrackHits = track->GetNTrackHits(); - AliDebug(1, Form("Number of track hits: %d ", nTrackHits)); - // track parameters at Vertex - trackParam = (*((AliMUONTrackParam*) track->GetTrackParamAtHit()->First())); - AliMUONTrackExtrap::ExtrapToVertex(&trackParam,0.,0.,0.); - x = trackParam.GetNonBendingCoor(); - y = trackParam.GetBendingCoor(); - z = trackParam.GetZ(); - bendingSlope = trackParam.GetBendingSlope(); - nonBendingSlope = trackParam.GetNonBendingSlope(); - pYZ = 1/TMath::Abs(trackParam.GetInverseBendingMomentum()); - pZ = pYZ/TMath::Sqrt(1+bendingSlope*bendingSlope); - pX = pZ * nonBendingSlope; - pY = pZ * bendingSlope; - c = TMath::Sign(1.0, trackParam.GetInverseBendingMomentum()); - AliDebug(1, Form("Track parameters at Vertex z= %f: X= %f Y= %f pX= %f pY= %f pZ= %f c= %f\n", - z, x, y, pX, pY, pZ, c)); - - // track parameters at first hit - trackParam1 = (AliMUONTrackParam*) track->GetTrackParamAtHit()->First(); - x = trackParam1->GetNonBendingCoor(); - y = trackParam1->GetBendingCoor(); - z = trackParam1->GetZ(); - bendingSlope = trackParam1->GetBendingSlope(); - nonBendingSlope = trackParam1->GetNonBendingSlope(); - pYZ = 1/TMath::Abs(trackParam1->GetInverseBendingMomentum()); - pZ = pYZ/TMath::Sqrt(1.0 + bendingSlope * bendingSlope); - pX = pZ * nonBendingSlope; - pY = pZ * bendingSlope; - c = TMath::Sign(1.0, trackParam1->GetInverseBendingMomentum()); - AliDebug(1, Form("track parameters at z= %f: X= %f Y= %f pX= %f pY= %f pZ= %f c= %f\n", - z, x, y, pX, pY, pZ, c)); + + // reset the seed (i.e. parameters at first cluster) before fitting + AliMUONTrackParam* firstTrackParam = (AliMUONTrackParam*) track.GetTrackParamAtCluster()->First(); + if (firstTrackParam->GetInverseBendingMomentum() == 0.) { + AliWarning("track parameters at first chamber are not initialized --> unable to refit"); + return kFALSE; + } + + // compute track parameters at each cluster from parameters at the first one + // necessary to compute multiple scattering effect during refitting + if (!track.UpdateTrackParamAtCluster()) { + AliWarning("bad track refitting due to extrapolation failure"); + return kFALSE; + } + + // Re-fit the track: + // Take into account the multiple scattering + // Calculate the track parameter covariance matrix + Fit(track, kTRUE, kFALSE, kTRUE); + + // Improve the reconstructed tracks if required + track.SetImproved(kFALSE); + if (enableImprovement && GetRecoParam()->ImproveTracks()) ImproveTrack(track); + + // Fill AliMUONTrack data members + if (track.GetGlobalChi2() < AliMUONTrack::MaxChi2()) return FinalizeTrack(track); + else { + AliWarning("track not finalized due to extrapolation failure"); + return kFALSE; } - // informations about generated particles NO !!!!!!!! - -// for (Int_t iPart = 0; iPart < np; iPart++) { -// p = gAlice->Particle(iPart); -// printf(" particle %d: type= %d px= %f py= %f pz= %f pdg= %d\n", -// iPart, p->GetPdgCode(), p->Px(), p->Py(), p->Pz(), p->GetPdgCode()); -// } - return; + } -