X-Git-Url: http://git.uio.no/git/?p=u%2Fmrichter%2FAliRoot.git;a=blobdiff_plain;f=MUON%2FAliMUONTrackReconstructor.cxx;h=2ce0754f3b176bac15276cefbc31ade5c28127e1;hp=3412ce42d4cebf453e4d1fc6072fcf4817ae1830;hb=96ebe67e954961b3c45067cd239082c0d43cd009;hpb=7ec3b9cf11e00fea765d34be1a04f9a0278a174c diff --git a/MUON/AliMUONTrackReconstructor.cxx b/MUON/AliMUONTrackReconstructor.cxx index 3412ce42d4c..2ce0754f3b1 100644 --- a/MUON/AliMUONTrackReconstructor.cxx +++ b/MUON/AliMUONTrackReconstructor.cxx @@ -15,6 +15,7 @@ /* $Id$ */ +//----------------------------------------------------------------------------- /// \class AliMUONTrackReconstructor /// MUON track reconstructor using the original method /// @@ -23,119 +24,101 @@ /// /// It contains as methods, among others: /// - MakeTracks to build the tracks -/// +//----------------------------------------------------------------------------- #include "AliMUONTrackReconstructor.h" #include "AliMUONConstants.h" -#include "AliMUONRawCluster.h" -#include "AliMUONHitForRec.h" -#include "AliMUONObjectPair.h" +#include "AliMUONVCluster.h" +#include "AliMUONVClusterStore.h" #include "AliMUONTrack.h" #include "AliMUONTrackParam.h" #include "AliMUONTrackExtrap.h" + #include "AliLog.h" #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); /// \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() : AliMUONVTrackReconstructor() { - /// Constructor for class AliMUONTrackReconstructor - - // Memory allocation for the TClonesArray of reconstructed tracks - fRecTracksPtr = new TClonesArray("AliMUONTrack", 10); -} - -//__________________________________________________________________________ -AliMUONTrackReconstructor::~AliMUONTrackReconstructor(void) -{ - /// Destructor for class AliMUONTrackReconstructor - delete fRecTracksPtr; + /// Constructor } -//__________________________________________________________________________ -void AliMUONTrackReconstructor::MakeTracks(void) + //__________________________________________________________________________ +AliMUONTrackReconstructor::~AliMUONTrackReconstructor() { - /// 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) +void AliMUONTrackReconstructor::MakeTrackCandidates(const 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"); // Loop over stations(1..) 5 and 4 and make track candidates - for (Int_t istat=4; istat>=3; istat--) - { + for (Int_t istat=4; istat>=3; istat--) { + // Make segments in the station - segments = MakeSegmentsInStation(istat); + segments = MakeSegmentsInStation(clusterStore,istat); + // Loop over segments 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])); 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(); + if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 2) || (AliLog::GetGlobalDebugLevel() >= 2)) { + cout<GetTrackParamAtCluster()->First())->GetCovariances().Print(); + } + + // Look for compatible cluster(s) in the other station + if (AliMUONReconstructor::GetRecoParam()->MakeTrackCandidatesFast()) + clusterFound = FollowLinearTrackInStation(*track, clusterStore, 7-istat); + else clusterFound = FollowTrackInStation(*track, clusterStore, 7-istat); + + // Remove track if no cluster found + if (!clusterFound) { + fRecTracksPtr->Remove(track); + fNRecTracks--; } - // 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 (AliMUONReconstructor::GetRecoParam()->TrackAllTracks()) RemoveIdenticalTracks(); else RemoveDoubleTracks(); AliDebug(1,Form("Number of good candidates = %d",fNRecTracks)); @@ -143,425 +126,721 @@ void AliMUONTrackReconstructor::MakeTrackCandidates(void) } //__________________________________________________________________________ -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) -{ - /// To remove double tracks: - /// Tracks are considered identical if more than half of the hits of the track - /// which has the smaller number of hits are in common with the other track. - /// Among two identical tracks, one keeps the track with the larger number of hits - /// or, if these numbers are equal, the track with the minimum chi2. - AliMUONTrack *track1, *track2, *trackToRemove; - Int_t hitsInCommon, nHits1, nHits2; - Bool_t removedTrack1; - // Loop over first track of the pair - track1 = (AliMUONTrack*) fRecTracksPtr->First(); - while (track1) { - removedTrack1 = kFALSE; - nHits1 = track1->GetNTrackHits(); - // Loop over second track of the pair - track2 = (AliMUONTrack*) fRecTracksPtr->After(track1); - while (track2) { - nHits2 = track2->GetNTrackHits(); - // number of hits in common between two tracks - hitsInCommon = track1->HitsInCommon(track2); - // check for identical tracks - if (((nHits1 < nHits2) && (2 * hitsInCommon > nHits1)) || (2 * hitsInCommon > nHits2)) { - // decide which track to remove - if ((nHits1 > nHits2) || ((nHits1 == nHits2) && (track1->GetFitFMin() <= track2->GetFitFMin()))) { - // remove track2 and continue the second loop with the track next to track2 - trackToRemove = track2; - track2 = (AliMUONTrack*) fRecTracksPtr->After(track2); - fRecTracksPtr->Remove(trackToRemove); - fRecTracksPtr->Compress(); // this is essential to retrieve the TClonesArray afterwards - fNRecTracks--; - } else { - // else remove track1 and continue the first loop with the track next to track1 - trackToRemove = track1; - track1 = (AliMUONTrack*) fRecTracksPtr->After(track1); - fRecTracksPtr->Remove(trackToRemove); - fRecTracksPtr->Compress(); // this is essential to retrieve the TClonesArray afterwards - fNRecTracks--; - removedTrack1 = kTRUE; - break; - } - } else track2 = (AliMUONTrack*) fRecTracksPtr->After(track2); - } // track2 - if (removedTrack1) continue; - track1 = (AliMUONTrack*) fRecTracksPtr->After(track1); - } // track1 - return; -} - - //__________________________________________________________________________ -void AliMUONTrackReconstructor::FollowTracks(void) +void AliMUONTrackReconstructor::FollowTracks(const 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; + Bool_t clusterFound; + + Double_t sigmaCut2 = AliMUONReconstructor::GetRecoParam()->GetSigmaCutForTracking() * + AliMUONReconstructor::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); + if (station==2) Fit(*track, kFALSE, kTRUE, kTRUE); + else Fit(*track, kFALSE, 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) { + if (track->GetNormalizedChi2() > sigmaCut2) { fRecTracksPtr->Remove(track); fNRecTracks--; continue; } - // Add MCS effects in parameter covariances - trackParamAtFirstHit = (AliMUONTrackParam*) (track->GetTrackParamAtHit()->First()); - AliMUONTrackExtrap::AddMCSEffect(trackParamAtFirstHit,AliMUONConstants::ChamberThicknessInX0(),1.); + + // save parameters from fit into smoothed parameters to complete track afterward + if (AliMUONReconstructor::GetRecoParam()->ComplementTracks()) { + + if (station==2) { // save track parameters on stations 4 and 5 + + // extrapolate track parameters and covariances at each cluster + track->UpdateCovTrackParamAtCluster(); + + // 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 + + // save parameters from fit + trackParam = (AliMUONTrackParam*) track->GetTrackParamAtCluster()->First(); + 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 + AliMUONTrackExtrap::ExtrapToZCov(nextTrackParam, nextTrackParam->GetClusterPtr()->GetZ()); + + // save them + nextTrackParam->SetSmoothParameters(nextTrackParam->GetParameters()); + nextTrackParam->SetSmoothCovariances(nextTrackParam->GetCovariances()); + + } + + } + + } + // Printout for debuging if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 2) || (AliLog::GetGlobalDebugLevel() >= 2)) { - cout<GetCovariances()->Print(); + cout<GetTrackParamAtCluster()->First())->GetCovariances().Print(); + } + + // Look for compatible cluster(s) in station(0..) "station" + clusterFound = FollowTrackInStation(*track, clusterStore, station); + + // Try to recover track if required + if (!clusterFound && AliMUONReconstructor::GetRecoParam()->RecoverTracks()) + clusterFound = RecoverTrack(*track, clusterStore, station); + + // remove track if no cluster found + if (!clusterFound) { + fRecTracksPtr->Remove(track); + fNRecTracks--; } - // Look for compatible hitForRec in station(0..) "station" - FollowTrackInStation(track,station); + } + // Compress fRecTracksPtr for the next step fRecTracksPtr->Compress(); + // Keep only the best tracks if required - if (!fgkTrackAllTracks) RemoveDoubleTracks(); + if (!AliMUONReconstructor::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--; } + + // save parameters from fit into smoothed parameters to complete track afterward + if (AliMUONReconstructor::GetRecoParam()->ComplementTracks()) { + + // save parameters from fit + trackParam = (AliMUONTrackParam*) track->GetTrackParamAtCluster()->First(); + 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 + AliMUONTrackExtrap::ExtrapToZCov(nextTrackParam, nextTrackParam->GetClusterPtr()->GetZ()); + + // save them + nextTrackParam->SetSmoothParameters(nextTrackParam->GetParameters()); + nextTrackParam->SetSmoothCovariances(nextTrackParam->GetCovariances()); + + } + + } + track = nextTrack; + } + fRecTracksPtr->Compress(); } //__________________________________________________________________________ -void AliMUONTrackReconstructor::FollowTrackInStation(AliMUONTrack* trackCandidate, Int_t nextStation) +Bool_t AliMUONTrackReconstructor::FollowTrackInStation(AliMUONTrack &trackCandidate, const AliMUONVClusterStore& clusterStore, Int_t nextStation) { - /// Follow trackCandidate in station(0..) nextStation and search for compatible HitForRec(s) + /// 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 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 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 = 1.e10; + Double_t chi2WithTwoClusters = 1.e10; + Double_t maxChi2WithOneCluster = 2. * AliMUONReconstructor::GetRecoParam()->GetSigmaCutForTracking() * + AliMUONReconstructor::GetRecoParam()->GetSigmaCutForTracking(); // 2 because 2 quantities in chi2 + Double_t maxChi2WithTwoClusters = 4. * AliMUONReconstructor::GetRecoParam()->GetSigmaCutForTracking() * + AliMUONReconstructor::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 extrapTrackParamAtCluster1; + AliMUONTrackParam extrapTrackParamAtCluster2; + AliMUONTrackParam bestTrackParamAtCluster1; + AliMUONTrackParam bestTrackParamAtCluster2; + + 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; + + // Get track parameters + AliMUONTrackParam extrapTrackParamAtCh(*(AliMUONTrackParam*)trackCandidate.GetTrackParamAtCluster()->First()); + + // Add MCS effect + AliMUONTrackExtrap::AddMCSEffect(&extrapTrackParamAtCh,AliMUONConstants::ChamberThicknessInX0(),1.); + + // Add MCS in the missing chamber if any (only 1 chamber can be missing according to tracking criteria) + if (ch1 < ch2 && extrapTrackParamAtCh.GetClusterPtr()->GetChamberId() > ch2 + 1) { + // extrapolation to the missing chamber + AliMUONTrackExtrap::ExtrapToZCov(&extrapTrackParamAtCh, AliMUONConstants::DefaultChamberZ(ch2 + 1)); + // add MCS effect + AliMUONTrackExtrap::AddMCSEffect(&extrapTrackParamAtCh,AliMUONConstants::ChamberThicknessInX0(),1.); + } + + //Extrapolate trackCandidate to chamber "ch2" + AliMUONTrackExtrap::ExtrapToZCov(&extrapTrackParamAtCh, AliMUONConstants::DefaultChamberZ(ch2)); + // 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) { + + // 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; + 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(),1.); + + // copy new track parameters for next step + extrapTrackParam = extrapTrackParamAtCluster2; + + //Extrapolate track parameters to chamber "ch1" + AliMUONTrackExtrap::ExtrapToZ(&extrapTrackParam, AliMUONConstants::DefaultChamberZ(ch1)); + + // reset cluster iterator of chamber 1 + nextInCh1.Reset(); + iCluster1 = -1; + + // look for second candidates in chamber 1 + while ( ( clusterCh1 = static_cast(nextInCh1()) ) ) { + iCluster1++; + + // try to add the current cluster fast + if (!TryOneClusterFast(extrapTrackParam, clusterCh1)) continue; + + // try to add the current cluster accuratly + 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; + } + + if (AliMUONReconstructor::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 (AliMUONReconstructor::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 (AliMUONReconstructor::GetRecoParam()->TrackAllTracks() || !foundTwoClusters) { + + // Printout for debuging if ((AliLog::GetDebugLevel("MUON","AliMUONTrackReconstructor") >= 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); + + // add MCS effect for next step + AliMUONTrackExtrap::AddMCSEffect(&extrapTrackParamAtCh,AliMUONConstants::ChamberThicknessInX0(),1.); + + //Extrapolate trackCandidate to chamber "ch1" + AliMUONTrackExtrap::ExtrapToZCov(&extrapTrackParamAtCh, AliMUONConstants::DefaultChamberZ(ch1)); + + // reset cluster iterator of chamber 1 + nextInCh1.Reset(); + iCluster1 = -1; + + // look for second candidates in chamber 1 + 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; + } + + if (AliMUONReconstructor::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 (!AliMUONReconstructor::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) + trackParamAtCluster2.SetParameters(trackParamAtCluster1.GetParameters()); + trackParamAtCluster2.SetZ(trackParamAtCluster1.GetZ()); + AliMUONTrackExtrap::ExtrapToZ(&trackParamAtCluster2, cluster2->GetZ()); + + // set pointer to cluster2 into trackParamAtCluster2 + trackParamAtCluster2.SetClusterPtr(cluster2); + + // 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); + + // add MCS effect to the covariance matrix at first cluster + AliMUONTrackExtrap::AddMCSEffect(&trackParam,AliMUONConstants::ChamberThicknessInX0(),1.); + + // 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); + 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)); + if (j == 4) dParam(j,0) *= TMath::Sign(1.,-paramAtCluster1Save(4,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 + AliMUONTrackExtrap::ExtrapToZ(&trackParam,cluster2->GetZ()); + + // 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 1.e10; + } + + // 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 + trackParamAtCluster.SetRemovable(kFALSE); + 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); + track.GetTrackParamAtCluster()->Sort(); + +} + + //__________________________________________________________________________ +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); + track.GetTrackParamAtCluster()->Sort(); + +} + + //__________________________________________________________________________ +Bool_t AliMUONTrackReconstructor::RecoverTrack(AliMUONTrack &trackCandidate, const AliMUONVClusterStore& clusterStore, Int_t nextStation) +{ + /// 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 = 0.; + + // 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 removable + if (!trackParamAtCluster->IsRemovable()) return kFALSE; + + // Pick up cluster with the worst chi2 + localChi2 = trackParamAtCluster->GetLocalChi2(); + if (localChi2 > worstLocalChi2) { + worstLocalChi2 = localChi2; + worstClusterNumber = clusterNumber; + } } + // Reset best cluster as being NOT removable + ((AliMUONTrackParam*)trackCandidate.GetTrackParamAtCluster()->UncheckedAt((worstClusterNumber+1)%2))->SetRemovable(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); + + // Look for new cluster(s) in next station + return FollowTrackInStation(trackCandidate,clusterStore,nextStation); + } //__________________________________________________________________________ -void AliMUONTrackReconstructor::SetVertexForFit(AliMUONTrack* trackCandidate) +void AliMUONTrackReconstructor::SetVertexErrXY2ForFit(AliMUONTrack &trackCandidate) { - /// Add the vertex as a measured hit to constrain the fit of the "trackCandidate" - /// Compute the vertex resolution from natural vertex dispersion and + /// 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 = AliMUONReconstructor::GetRecoParam()->GetNonBendingVertexDispersion() * + AliMUONReconstructor::GetRecoParam()->GetNonBendingVertexDispersion(); + Double_t bendingReso2 = AliMUONReconstructor::GetRecoParam()->GetBendingVertexDispersion() * + AliMUONReconstructor::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); + 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". 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; @@ -570,7 +849,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; @@ -587,12 +866,24 @@ 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.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()); + trackParam = (AliMUONTrackParam*) (track.GetTrackParamAtCluster()->First()); // could be tried with no limits for the search (min=max=0) ???? // mandatory limits in non Bending to avoid NaN values of parameters gMinuit->mnparm(0, "X", trackParam->GetNonBendingCoor(), 0.03, -500.0, 500.0, status); @@ -621,8 +912,8 @@ 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) { @@ -632,336 +923,279 @@ void AliMUONTrackReconstructor::Fit(AliMUONTrack *track, Bool_t includeMCS, Bool 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]); + trackBeingFitted->UpdateTrackParamAtCluster(); // 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); + if (nonBendingReso2 == 0. || bendingReso2 == 0.) chi2 += 1.e10; + else { + AliMUONTrackParam paramAtVertex(*trackParamAtCluster); + AliMUONTrackExtrap::ExtrapToZ(¶mAtVertex, 0.); // vextex position = (0,0,0) + dX = paramAtVertex.GetNonBendingCoor(); + dY = paramAtVertex.GetBendingCoor(); + chi2 += dX * dX / nonBendingReso2 + dY * dY / bendingReso2; } - 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*/) +void 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 + 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 = AliMUONReconstructor::GetRecoParam()->GetSigmaCutForTracking() * + AliMUONReconstructor::GetRecoParam()->GetSigmaCutForTracking(); + Bool_t foundOneCluster, trackModified; + 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]); - } - // 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); + // Remove double track to complete only "good" tracks + RemoveDoubleTracks(); + + 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; + } + } - } // for (hitNumber2 = hitNumber1;... - } // for (hitNumber1 = 0;... - - // 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); - - // 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)); + + // add new cluster if any + if (foundOneCluster) { + UpdateTrack(*track,bestTrackParamAtCluster); + bestTrackParamAtCluster.SetAloneInChamber(kFALSE); + trackParam->SetAloneInChamber(kFALSE); + trackModified = kTRUE; } + + trackParam = nextTrackParam; } - } 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()); - } + + // re-fit track parameters if needed + if (trackModified) Fit(*track, kTRUE, kFALSE, kTRUE); + + track = (AliMUONTrack*) fRecTracksPtr->After(track); } - delete covBending; - delete covNonBending; - delete [] msa2; -} - - //__________________________________________________________________________ -Double_t MultipleScatteringAngle2(AliMUONTrackParam *param) -{ - /// Returns square of multiple Coulomb scattering angle - /// from TrackParamAtHit pointed to by "param" - Double_t slopeBending, slopeNonBending, radiationLength, inverseBendingMomentum2, inverseTotalMomentum2; - Double_t varMultipleScatteringAngle; - 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; } //__________________________________________________________________________ -void AliMUONTrackReconstructor::FillMUONTrack(void) +void AliMUONTrackReconstructor::ImproveTracks() { - /// 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; + /// Improve tracks by removing clusters with local chi2 highter than the defined cut + /// Recompute track parameters and covariances at the remaining clusters + AliDebug(1,"Enter ImproveTracks"); + + Double_t localChi2, worstLocalChi2; + Int_t worstChamber, previousChamber; + AliMUONTrack *track, *nextTrack; + AliMUONTrackParam *trackParamAtCluster, *worstTrackParamAtCluster, *previousTrackParam, *nextTrackParam; + Double_t sigmaCut2 = AliMUONReconstructor::GetRecoParam()->GetSigmaCutForImprovement() * + AliMUONReconstructor::GetRecoParam()->GetSigmaCutForImprovement(); + + // Remove double track to improve only "good" tracks + RemoveDoubleTracks(); 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)); + + // prepare next track in case the actual track is suppressed + nextTrack = (AliMUONTrack*) fRecTracksPtr->After(track); + + while (!track->IsImproved()) { + + // Update track parameters and covariances + track->UpdateCovTrackParamAtCluster(); + + // 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); + } + + // Check if bad cluster found + if (!worstTrackParamAtCluster) { + track->SetImproved(kTRUE); + 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 remove the entire track + if (!worstTrackParamAtCluster->IsRemovable() && worstTrackParamAtCluster->IsAloneInChamber()) { + fRecTracksPtr->Remove(track); + fNRecTracks--; + break; + } + + // Reset the second cluster in the same station as being not removable + // or reset the second cluster in the same chamber as being alone + worstChamber = worstTrackParamAtCluster->GetClusterPtr()->GetChamberId(); + previousTrackParam = (AliMUONTrackParam*) track->GetTrackParamAtCluster()->Before(worstTrackParamAtCluster); + nextTrackParam = (AliMUONTrackParam*) track->GetTrackParamAtCluster()->After(worstTrackParamAtCluster); + if (worstTrackParamAtCluster->IsAloneInChamber()) { // Worst cluster removable and alone in chamber + + if (worstChamber%2 == 0) { // Modify flags in next chamber + + nextTrackParam->SetRemovable(kFALSE); + if (!nextTrackParam->IsAloneInChamber()) // Make sure both clusters in second chamber are not removable anymore + ((AliMUONTrackParam*) track->GetTrackParamAtCluster()->After(nextTrackParam))->SetRemovable(kFALSE); + + } else { // Modify flags in previous chamber + + previousTrackParam->SetRemovable(kFALSE); + if (!previousTrackParam->IsAloneInChamber()) // Make sure both clusters in second chamber are not removable anymore + ((AliMUONTrackParam*) track->GetTrackParamAtCluster()->Before(previousTrackParam))->SetRemovable(kFALSE); + + } + + } else { // Worst cluster not alone in its chamber + + if (previousTrackParam) previousChamber = previousTrackParam->GetClusterPtr()->GetChamberId(); + else previousChamber = -1; + + if (previousChamber == worstChamber) { // the second cluster on the same chamber is the previous one + + previousTrackParam->SetAloneInChamber(kTRUE); + // transfert the removability to the second cluster + if (worstTrackParamAtCluster->IsRemovable()) previousTrackParam->SetRemovable(kTRUE); + + } else { // the second cluster on the same chamber is the next one + + nextTrackParam->SetAloneInChamber(kTRUE); + // transfert the removability to the second cluster + if (worstTrackParamAtCluster->IsRemovable()) nextTrackParam->SetRemovable(kTRUE); + + } + + } + + // 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; + } + } - track = (AliMUONTrack*) fRecTracksPtr->After(track); + + track = nextTrack; } - return; + + // compress the array in case of some tracks have been removed + fRecTracksPtr->Compress(); + } //__________________________________________________________________________ -void AliMUONTrackReconstructor::EventDump(void) +void AliMUONTrackReconstructor::Finalize() { - /// Dump reconstructed event (track parameters at vertex and at first hit), - /// and the particle parameters + /// Recompute track parameters and covariances at each attached cluster from those at the first one + 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)); + + track = (AliMUONTrack*) fRecTracksPtr->First(); + while (track) { + + // update track parameters if not already done + if (!track->IsImproved()) track->UpdateCovTrackParamAtCluster(); + + track = (AliMUONTrack*) fRecTracksPtr->After(track); + } - // 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; + } -