X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=MUON%2FAliMUONVTrackReconstructor.cxx;h=6eedfe7bb1ce7f5a868a59e41c5d82025dd8a260;hb=0cd74160756eda9c680afc5ee205d6b795914e89;hp=50dd89ffec0ddcedf0a010b42c2eb47cb187098e;hpb=7771752ec862f31222329412e915e79197b8cea3;p=u%2Fmrichter%2FAliRoot.git diff --git a/MUON/AliMUONVTrackReconstructor.cxx b/MUON/AliMUONVTrackReconstructor.cxx index 50dd89ffec0..6eedfe7bb1c 100644 --- a/MUON/AliMUONVTrackReconstructor.cxx +++ b/MUON/AliMUONVTrackReconstructor.cxx @@ -15,28 +15,50 @@ /* $Id$ */ -//////////////////////////////////// -/// +//----------------------------------------------------------------------------- /// \class AliMUONVTrackReconstructor /// Virtual MUON track reconstructor in ALICE (class renamed from AliMUONEventReconstructor) /// -/// This class contains as data: -/// * a pointer to the array of hits to be reconstructed (the event) -/// * a pointer to the array of segments made with these hits inside each station -/// * a pointer to the array of reconstructed tracks +/// This class contains as data a pointer to the array of reconstructed tracks /// /// It contains as methods, among others: /// * EventReconstruct to build the muon tracks /// * EventReconstructTrigger to build the trigger tracks +/// * ValidateTracksWithTrigger to match tracker/trigger tracks /// -/// \author Philippe Pillot +/// Several options and adjustable parameters are available for both KALMAN and ORIGINAL +/// tracking algorithms. They can be changed by using: +/// AliMUONRecoParam *muonRecoParam = AliMUONRecoParam::GetLow(High)FluxParam(); +/// muonRecoParam->Set...(); // see methods in AliMUONRecoParam.h for details +/// AliMUONReconstructor::SetRecoParam(muonRecoParam); +/// +/// Main parameters and options are: +/// - *fgkSigmaToCutForTracking* : quality cut used to select new clusters to be +/// attached to the track candidate and to select good tracks. +/// - *fgkMakeTrackCandidatesFast* : if this flag is set to 'true', the track candidates +/// are made assuming linear propagation between stations 4 and 5. +/// - *fgkTrackAllTracks* : according to the value of this flag, in case that several +/// new clusters pass the quality cut, either we consider all the possibilities +/// (duplicating tracks) or we attach only the best cluster. +/// - *fgkRecoverTracks* : if this flag is set to 'true', we try to recover the tracks +/// lost during the tracking by removing the worst of the 2 clusters attached in the +/// previous station. +/// - *fgkImproveTracks* : if this flag is set to 'true', we try to improve the quality +/// of the tracks at the end of the tracking by removing clusters that do not pass +/// new quality cut (the track is removed is it does not contain enough cluster anymore). +/// - *fgkComplementTracks* : if this flag is set to 'true', we try to improve the quality +/// of the tracks at the end of the tracking by adding potentially missing clusters +/// (we may have 2 clusters in the same chamber because of the overlapping of detection +/// elements, which is not handle by the tracking algorithm). +/// - *fgkSigmaToCutForImprovement* : quality cut used when we try to improve the +/// quality of the tracks. /// -//////////////////////////////////// +/// \author Philippe Pillot +//----------------------------------------------------------------------------- #include "AliMUONVTrackReconstructor.h" -#include "AliMUONData.h" + #include "AliMUONConstants.h" -#include "AliMUONHitForRec.h" #include "AliMUONObjectPair.h" #include "AliMUONTriggerTrack.h" #include "AliMUONTriggerCircuit.h" @@ -46,421 +68,977 @@ #include "AliMUONTrackParam.h" #include "AliMUONTrackExtrap.h" #include "AliMUONTrackHitPattern.h" +#include "AliMUONVTrackStore.h" +#include "AliMUONVClusterStore.h" +#include "AliMUONVCluster.h" +#include "AliMUONVClusterServer.h" +#include "AliMUONVTriggerStore.h" +#include "AliMUONVTriggerTrackStore.h" +#include "AliMpDEManager.h" +#include "AliMpArea.h" #include "AliLog.h" +#include "AliCodeTimer.h" #include "AliTracker.h" -#include #include #include +#include +#include + +#include /// \cond CLASSIMP ClassImp(AliMUONVTrackReconstructor) // Class implementation in ROOT context /// \endcond -//************* Defaults parameters for reconstruction -const Double_t AliMUONVTrackReconstructor::fgkDefaultMinBendingMomentum = 3.0; -const Double_t AliMUONVTrackReconstructor::fgkDefaultMaxBendingMomentum = 3000.0; -const Double_t AliMUONVTrackReconstructor::fgkDefaultBendingResolution = 0.01; -const Double_t AliMUONVTrackReconstructor::fgkDefaultNonBendingResolution = 0.144; -const Double_t AliMUONVTrackReconstructor::fgkDefaultBendingVertexDispersion = 10.; -const Double_t AliMUONVTrackReconstructor::fgkDefaultNonBendingVertexDispersion = 10.; -const Double_t AliMUONVTrackReconstructor::fgkDefaultMaxNormChi2MatchTrigger = 16.0; - -//__________________________________________________________________________ -AliMUONVTrackReconstructor::AliMUONVTrackReconstructor(AliMUONData* data) + //__________________________________________________________________________ +AliMUONVTrackReconstructor::AliMUONVTrackReconstructor(AliMUONVClusterServer& clusterServer) : TObject(), - fMinBendingMomentum(fgkDefaultMinBendingMomentum), - fMaxBendingMomentum(fgkDefaultMaxBendingMomentum), - fBendingResolution(fgkDefaultBendingResolution), - fNonBendingResolution(fgkDefaultNonBendingResolution), - fBendingVertexDispersion(fgkDefaultBendingVertexDispersion), - fNonBendingVertexDispersion(fgkDefaultNonBendingVertexDispersion), - fMaxNormChi2MatchTrigger(fgkDefaultMaxNormChi2MatchTrigger), - fHitsForRecPtr(0x0), - fNHitsForRec(0), - fNHitsForRecPerChamber(0x0), - fIndexOfFirstHitForRecPerChamber(0x0), fRecTracksPtr(0x0), fNRecTracks(0), - fMUONData(data), - fTriggerTrack(new AliMUONTriggerTrack()), - fTriggerCircuit(0x0) + fClusterServer(clusterServer) { /// Constructor for class AliMUONVTrackReconstructor - fNHitsForRecPerChamber = new Int_t[AliMUONConstants::NTrackingCh()]; - fIndexOfFirstHitForRecPerChamber = new Int_t[AliMUONConstants::NTrackingCh()]; - - // Memory allocation for the TClonesArray of hits for reconstruction - // Is 10000 the right size ???? - fHitsForRecPtr = new TClonesArray("AliMUONHitForRec", 10000); - + + // Memory allocation for the TClonesArray of reconstructed tracks + fRecTracksPtr = new TClonesArray("AliMUONTrack", 100); + // set the magnetic field for track extrapolations const AliMagF* kField = AliTracker::GetFieldMap(); if (!kField) AliFatal("No field available"); AliMUONTrackExtrap::SetField(kField); - - fTrackHitPattern = new AliMUONTrackHitPattern(fMUONData); } //__________________________________________________________________________ -AliMUONVTrackReconstructor::~AliMUONVTrackReconstructor(void) +AliMUONVTrackReconstructor::~AliMUONVTrackReconstructor() { /// Destructor for class AliMUONVTrackReconstructor - delete [] fNHitsForRecPerChamber; - delete [] fIndexOfFirstHitForRecPerChamber; - delete fTriggerTrack; - delete fHitsForRecPtr; - delete fTrackHitPattern; + delete fRecTracksPtr; } //__________________________________________________________________________ -void AliMUONVTrackReconstructor::EventReconstruct(void) +void AliMUONVTrackReconstructor::ResetTracks() { - /// To reconstruct one event - AliDebug(1,"Enter EventReconstruct"); - - ResetTracks(); //AZ - ResetHitsForRec(); //AZ - AddHitsForRecFromRawClusters(); - MakeTracks(); - if (fMUONData->IsTriggerTrackBranchesInTree()) - ValidateTracksWithTrigger(); - - fTrackHitPattern->GetHitPattern(fRecTracksPtr); + /// To reset the TClonesArray of reconstructed tracks + if (fRecTracksPtr) fRecTracksPtr->Clear("C"); + fNRecTracks = 0; + return; +} + //__________________________________________________________________________ +void AliMUONVTrackReconstructor::EventReconstruct(AliMUONVClusterStore& clusterStore, AliMUONVTrackStore& trackStore) +{ + /// To reconstruct one event + AliDebug(1,""); + AliCodeTimerAuto(""); + + // Reset array of tracks + ResetTracks(); + + // Look for candidates from clusters in stations(1..) 4 and 5 + MakeTrackCandidates(clusterStore); + + // Look for extra candidates from clusters in stations(1..) 4 and 5 + if (AliMUONReconstructor::GetRecoParam()->MakeMoreTrackCandidates()) MakeMoreTrackCandidates(clusterStore); + + // Stop tracking if no candidate found + if (fRecTracksPtr->GetEntriesFast() == 0) return; + + // Follow tracks in stations(1..) 3, 2 and 1 + FollowTracks(clusterStore); + + // Complement the reconstructed tracks + if (AliMUONReconstructor::GetRecoParam()->ComplementTracks()) ComplementTracks(clusterStore); + + // Improve the reconstructed tracks + if (AliMUONReconstructor::GetRecoParam()->ImproveTracks()) ImproveTracks(); + + // Remove double tracks + RemoveDoubleTracks(); + + // Fill AliMUONTrack data members + Finalize(); + // Add tracks to MUON data container - for(Int_t i=0; iAt(i); - fMUONData->AddRecTrack(*track); + track->SetUniqueID(i+1); + trackStore.Add(*track); } } //__________________________________________________________________________ -void AliMUONVTrackReconstructor::ResetTracks(void) +TClonesArray* AliMUONVTrackReconstructor::MakeSegmentsBetweenChambers(const AliMUONVClusterStore& clusterStore, Int_t ch1, Int_t ch2) { - /// To reset the TClonesArray of reconstructed tracks - if (fRecTracksPtr) fRecTracksPtr->Delete(); - // Delete in order that the Track destructors are called, - // hence the space for the TClonesArray of pointers to TrackHit's is freed - fNRecTracks = 0; - return; + /// To make the list of segments from the list of clusters in the 2 given chambers. + /// Return a new TClonesArray of segments. + /// It is the responsibility of the user to delete it afterward. + AliDebug(1,Form("Enter MakeSegmentsBetweenChambers (1..) %d-%d", ch1+1, ch2+1)); + + AliMUONVCluster *cluster1, *cluster2; + AliMUONObjectPair *segment; + Double_t nonBendingSlope = 0, bendingSlope = 0, impactParam = 0., bendingMomentum = 0.; // to avoid compilation warning + + // Create iterators to loop over clusters in both chambers + TIter nextInCh1(clusterStore.CreateChamberIterator(ch1,ch1)); + TIter nextInCh2(clusterStore.CreateChamberIterator(ch2,ch2)); + + // list of segments + TClonesArray *segments = new TClonesArray("AliMUONObjectPair", 100); + + // Loop over clusters in the first chamber of the station + while ( ( cluster1 = static_cast(nextInCh1()) ) ) { + + // reset cluster iterator of chamber 2 + nextInCh2.Reset(); + + // Loop over clusters in the second chamber of the station + while ( ( cluster2 = static_cast(nextInCh2()) ) ) { + + // non bending slope + nonBendingSlope = (cluster1->GetX() - cluster2->GetX()) / (cluster1->GetZ() - cluster2->GetZ()); + + // bending slope + bendingSlope = (cluster1->GetY() - cluster2->GetY()) / (cluster1->GetZ() - cluster2->GetZ()); + + // impact parameter + impactParam = cluster1->GetY() - cluster1->GetZ() * bendingSlope; + + // absolute value of bending momentum + bendingMomentum = TMath::Abs(AliMUONTrackExtrap::GetBendingMomentumFromImpactParam(impactParam)); + + // check for non bending slope and bending momentum within tolerances + if (TMath::Abs(nonBendingSlope) < AliMUONReconstructor::GetRecoParam()->GetMaxNonBendingSlope() && + bendingMomentum < AliMUONReconstructor::GetRecoParam()->GetMaxBendingMomentum() && + bendingMomentum > AliMUONReconstructor::GetRecoParam()->GetMinBendingMomentum()) { + + // make new segment + segment = new ((*segments)[segments->GetLast()+1]) AliMUONObjectPair(cluster1, cluster2, kFALSE, kFALSE); + + // Printout for debug + if (AliLog::GetGlobalDebugLevel() > 1) { + cout << "segmentIndex(0...): " << segments->GetLast() << endl; + segment->Dump(); + cout << "Cluster in first chamber" << endl; + cluster1->Print(); + cout << "Cluster in second chamber" << endl; + cluster2->Print(); + } + + } + + } + + } + + // Printout for debug + AliDebug(1,Form("chambers%d-%d: NSegments = %d ", ch1+1, ch2+1, segments->GetEntriesFast())); + + return segments; } //__________________________________________________________________________ -void AliMUONVTrackReconstructor::ResetHitsForRec(void) +void AliMUONVTrackReconstructor::RemoveIdenticalTracks() { - /// To reset the array and the number of HitsForRec, - /// and also the number of HitsForRec - /// and the index of the first HitForRec per chamber - if (fHitsForRecPtr) fHitsForRecPtr->Delete(); - fNHitsForRec = 0; - for (Int_t ch = 0; ch < AliMUONConstants::NTrackingCh(); ch++) - fNHitsForRecPerChamber[ch] = fIndexOfFirstHitForRecPerChamber[ch] = 0; + /// To remove identical tracks: + /// Tracks are considered identical if they have all their clusters in common. + /// One keeps the track with the larger number of clusters if need be + AliMUONTrack *track1, *track2, *trackToRemove; + Int_t clustersInCommon, nClusters1, nClusters2; + Bool_t removedTrack1; + // Loop over first track of the pair + track1 = (AliMUONTrack*) fRecTracksPtr->First(); + while (track1) { + removedTrack1 = kFALSE; + nClusters1 = track1->GetNClusters(); + // Loop over second track of the pair + track2 = (AliMUONTrack*) fRecTracksPtr->After(track1); + while (track2) { + nClusters2 = track2->GetNClusters(); + // number of clusters in common between two tracks + clustersInCommon = track1->ClustersInCommon(track2); + // check for identical tracks + if ((clustersInCommon == nClusters1) || (clustersInCommon == nClusters2)) { + // decide which track to remove + if (nClusters2 > nClusters1) { + // 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 AliMUONVTrackReconstructor::SortHitsForRecWithIncreasingChamber() +void AliMUONVTrackReconstructor::RemoveDoubleTracks() { - /// Sort HitsForRec's in increasing order with respect to chamber number. - /// Uses the function "Compare". - /// Update the information for HitsForRec per chamber too. - Int_t ch, nhits, prevch; - fHitsForRecPtr->Sort(); - for (ch = 0; ch < AliMUONConstants::NTrackingCh(); ch++) { - fNHitsForRecPerChamber[ch] = 0; - fIndexOfFirstHitForRecPerChamber[ch] = 0; - } - prevch = 0; // previous chamber - nhits = 0; // number of hits in current chamber - // Loop over HitsForRec - for (Int_t hit = 0; hit < fNHitsForRec; hit++) { - // chamber number (0...) - ch = ((AliMUONHitForRec*) ((*fHitsForRecPtr)[hit]))->GetChamberNumber(); - // increment number of hits in current chamber - (fNHitsForRecPerChamber[ch])++; - // update index of first HitForRec in current chamber - // if chamber number different from previous one - if (ch != prevch) { - fIndexOfFirstHitForRecPerChamber[ch] = hit; - prevch = ch; - } - } + /// To remove double tracks: + /// Tracks are considered identical if more than half of the clusters of the track + /// which has the smaller number of clusters are in common with the other track. + /// Among two identical tracks, one keeps the track with the larger number of clusters + /// or, if these numbers are equal, the track with the minimum chi2. + AliMUONTrack *track1, *track2, *trackToRemove; + Int_t clustersInCommon, nClusters1, nClusters2; + Bool_t removedTrack1; + // Loop over first track of the pair + track1 = (AliMUONTrack*) fRecTracksPtr->First(); + while (track1) { + removedTrack1 = kFALSE; + nClusters1 = track1->GetNClusters(); + // Loop over second track of the pair + track2 = (AliMUONTrack*) fRecTracksPtr->After(track1); + while (track2) { + nClusters2 = track2->GetNClusters(); + // number of clusters in common between two tracks + clustersInCommon = track1->ClustersInCommon(track2); + // check for identical tracks + if (((nClusters1 < nClusters2) && (2 * clustersInCommon > nClusters1)) || (2 * clustersInCommon > nClusters2)) { + // decide which track to remove + if ((nClusters1 > nClusters2) || ((nClusters1 == nClusters2) && (track1->GetGlobalChi2() <= track2->GetGlobalChi2()))) { + // 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; } //__________________________________________________________________________ -TClonesArray* AliMUONVTrackReconstructor::MakeSegmentsInStation(Int_t station) +void AliMUONVTrackReconstructor::AskForNewClustersInChamber(const AliMUONTrackParam &trackParam, + AliMUONVClusterStore& clusterStore, Int_t chamber) { - /// To make the list of segments in station(0..) "Station" from the list of hits to be reconstructed. - /// Return a new TClonesArray of segments. - /// It is the responsibility of the user to delete it afterward. - AliDebug(1,Form("Enter MakeSegmentsPerStation (0...) %d",station)); + /// Ask the clustering to reconstruct new clusters around the track candidate position + + // check if the current chamber is useable + if (!AliMUONReconstructor::GetRecoParam()->UseChamber(chamber)) return; + + // maximum distance between the center of the chamber and a detection element + // (accounting for the inclination of the chamber) + static const Double_t gkMaxDZ = 15.; // 15 cm + + // extrapolate track parameters to the chamber + AliMUONTrackParam extrapTrackParam(trackParam); + AliMUONTrackExtrap::ExtrapToZCov(&extrapTrackParam, AliMUONConstants::DefaultChamberZ(chamber)); + + // build the searching area using the track resolution and the maximum-distance-to-track value + const TMatrixD& kParamCov = extrapTrackParam.GetCovariances(); + Double_t errX2 = kParamCov(0,0) + gkMaxDZ * gkMaxDZ * kParamCov(1,1) + 2. * gkMaxDZ * TMath::Abs(kParamCov(0,1)); + Double_t errY2 = kParamCov(2,2) + gkMaxDZ * gkMaxDZ * kParamCov(3,3) + 2. * gkMaxDZ * TMath::Abs(kParamCov(2,3)); + Double_t dX = TMath::Abs(trackParam.GetNonBendingSlope()) * gkMaxDZ + + AliMUONReconstructor::GetRecoParam()->GetMaxNonBendingDistanceToTrack() + + AliMUONReconstructor::GetRecoParam()->GetSigmaCutForTracking() * TMath::Sqrt(errX2); + Double_t dY = TMath::Abs(trackParam.GetBendingSlope()) * gkMaxDZ + + AliMUONReconstructor::GetRecoParam()->GetMaxBendingDistanceToTrack() + + AliMUONReconstructor::GetRecoParam()->GetSigmaCutForTracking() * TMath::Sqrt(errY2); + TVector2 dimensions(dX, dY); + TVector2 position(extrapTrackParam.GetNonBendingCoor(), extrapTrackParam.GetBendingCoor()); + AliMpArea area(position, dimensions); + + // ask to cluterize in the given area of the given chamber + fClusterServer.Clusterize(chamber, clusterStore, area); - AliMUONHitForRec *hit1Ptr, *hit2Ptr; - AliMUONObjectPair *segment; - Double_t bendingSlope = 0, impactParam = 0., bendingMomentum = 0.; // to avoid compilation warning - // first and second chambers (0...) in the station - Int_t ch1 = 2 * station; - Int_t ch2 = ch1 + 1; - // list of segments - TClonesArray *segments = new TClonesArray("AliMUONObjectPair", fNHitsForRecPerChamber[ch2]); - // Loop over HitForRec's in the first chamber of the station - for (Int_t hit1 = fIndexOfFirstHitForRecPerChamber[ch1]; - hit1 < fIndexOfFirstHitForRecPerChamber[ch1] + fNHitsForRecPerChamber[ch1]; - hit1++) { - // pointer to the HitForRec - hit1Ptr = (AliMUONHitForRec*) ((*fHitsForRecPtr)[hit1]); - // Loop over HitsForRec's in the second chamber of the station - for (Int_t hit2 = fIndexOfFirstHitForRecPerChamber[ch2]; - hit2 < fIndexOfFirstHitForRecPerChamber[ch2] + fNHitsForRecPerChamber[ch2]; - hit2++) { - // pointer to the HitForRec - hit2Ptr = (AliMUONHitForRec*) ((*fHitsForRecPtr)[hit2]); - if ( hit1Ptr->GetZ() - hit2Ptr->GetZ() != 0. ) { - // bending slope - bendingSlope = (hit1Ptr->GetBendingCoor() - hit2Ptr->GetBendingCoor()) / (hit1Ptr->GetZ() - hit2Ptr->GetZ()); - // impact parameter - impactParam = hit1Ptr->GetBendingCoor() - hit1Ptr->GetZ() * bendingSlope; - // absolute value of bending momentum - bendingMomentum = TMath::Abs(AliMUONTrackExtrap::GetBendingMomentumFromImpactParam(impactParam)); - } else { - AliWarning("hit1Ptr->GetZ() = hit2Ptr->GetZ(): no segment created"); - continue; - } - // check for bending momentum within tolerances - if ((bendingMomentum < fMaxBendingMomentum) && (bendingMomentum > fMinBendingMomentum)) { - // make new segment - segment = new ((*segments)[segments->GetLast()+1]) AliMUONObjectPair(hit1Ptr, hit2Ptr, kFALSE, kFALSE); - if (AliLog::GetGlobalDebugLevel() > 1) { - cout << "segmentIndex(0...): " << segments->GetLast() << endl; - segment->Dump(); - cout << "HitForRec in first chamber" << endl; - hit1Ptr->Dump(); - cout << "HitForRec in second chamber" << endl; - hit2Ptr->Dump(); - } - } - } //for (Int_t hit2 - } // for (Int_t hit1... - AliDebug(1,Form("Station: %d NSegments: %d ", station, segments->GetEntriesFast())); - return segments; } //__________________________________________________________________________ -void AliMUONVTrackReconstructor::ValidateTracksWithTrigger(void) +void AliMUONVTrackReconstructor::AskForNewClustersInStation(const AliMUONTrackParam &trackParam, + AliMUONVClusterStore& clusterStore, Int_t station) { - /// Try to match track from tracking system with trigger track - static const Double_t kDistSigma[3]={1,1,0.02}; // sigma of distributions (trigger-track) X,Y,slopeY + /// Ask the clustering to reconstruct new clusters around the track candidate position + /// in the 2 chambers of the given station + AskForNewClustersInChamber(trackParam, clusterStore, 2*station+1); + AskForNewClustersInChamber(trackParam, clusterStore, 2*station); +} + + //__________________________________________________________________________ +Double_t AliMUONVTrackReconstructor::TryOneCluster(const AliMUONTrackParam &trackParam, AliMUONVCluster* cluster, + AliMUONTrackParam &trackParamAtCluster, Bool_t updatePropagator) +{ +/// Test the compatibility between the track and the cluster (using trackParam's covariance matrix): +/// return the corresponding Chi2 +/// return trackParamAtCluster - AliMUONTrack *track; - AliMUONTrackParam trackParam; - AliMUONTriggerTrack *triggerTrack; - AliMUONLocalTrigger *locTrg; + // extrapolate track parameters and covariances at the z position of the tested cluster + trackParamAtCluster = trackParam; + AliMUONTrackExtrap::ExtrapToZCov(&trackParamAtCluster, cluster->GetZ(), updatePropagator); + + // set pointer to cluster into trackParamAtCluster + trackParamAtCluster.SetClusterPtr(cluster); + + // Set differences between trackParam and cluster in the bending and non bending directions + Double_t dX = cluster->GetX() - trackParamAtCluster.GetNonBendingCoor(); + Double_t dY = cluster->GetY() - trackParamAtCluster.GetBendingCoor(); + + // Calculate errors and covariances + const TMatrixD& kParamCov = trackParamAtCluster.GetCovariances(); + Double_t sigmaX2 = kParamCov(0,0) + cluster->GetErrX2(); + Double_t sigmaY2 = kParamCov(2,2) + cluster->GetErrY2(); - fMUONData->SetTreeAddress("RL"); - fMUONData->GetRecTriggerTracks(); - TClonesArray *recTriggerTracks = fMUONData->RecTriggerTracks(); + // Compute chi2 + return dX * dX / sigmaX2 + dY * dY / sigmaY2; + +} - fMUONData->SetTreeAddress("TC"); - fMUONData->GetTrigger(); - TClonesArray *localTrigger = fMUONData->LocalTrigger(); + //__________________________________________________________________________ +Bool_t AliMUONVTrackReconstructor::TryOneClusterFast(const AliMUONTrackParam &trackParam, AliMUONVCluster* cluster) +{ +/// Test the compatibility between the track and the cluster +/// given the track resolution + the maximum-distance-to-track value +/// and assuming linear propagation of the track: +/// return kTRUE if they are compatibles + + Double_t dZ = cluster->GetZ() - trackParam.GetZ(); + Double_t dX = cluster->GetX() - (trackParam.GetNonBendingCoor() + trackParam.GetNonBendingSlope() * dZ); + Double_t dY = cluster->GetY() - (trackParam.GetBendingCoor() + trackParam.GetBendingSlope() * dZ); + const TMatrixD& kParamCov = trackParam.GetCovariances(); + Double_t errX2 = kParamCov(0,0) + dZ * dZ * kParamCov(1,1) + 2. * dZ * kParamCov(0,1); + Double_t errY2 = kParamCov(2,2) + dZ * dZ * kParamCov(3,3) + 2. * dZ * kParamCov(2,3); - Int_t matchTrigger; - Int_t loTrgNum; - Double_t distTriggerTrack[3]; - Double_t xTrack, yTrack, ySlopeTrack, chi2MatchTrigger, minChi2MatchTrigger, chi2; + Double_t dXmax = AliMUONReconstructor::GetRecoParam()->GetSigmaCutForTracking() * TMath::Sqrt(errX2) + + AliMUONReconstructor::GetRecoParam()->GetMaxNonBendingDistanceToTrack(); + Double_t dYmax = AliMUONReconstructor::GetRecoParam()->GetSigmaCutForTracking() * TMath::Sqrt(errY2) + + AliMUONReconstructor::GetRecoParam()->GetMaxBendingDistanceToTrack(); - track = (AliMUONTrack*) fRecTracksPtr->First(); - while (track) { - matchTrigger = -1; - chi2MatchTrigger = 0.; - loTrgNum = -1; - Int_t doubleMatch=-1; // Check if track matches 2 trigger tracks - Double_t doubleChi2 = -1.; + if (TMath::Abs(dX) > dXmax || TMath::Abs(dY) > dYmax) return kFALSE; + + return kTRUE; + +} + + //__________________________________________________________________________ +Double_t AliMUONVTrackReconstructor::TryTwoClustersFast(const AliMUONTrackParam &trackParamAtCluster1, AliMUONVCluster* cluster2, + AliMUONTrackParam &trackParamAtCluster2) +{ +/// Test the compatibility between the track and the 2 clusters together (using trackParam's covariance matrix) +/// assuming linear propagation between the two clusters: +/// return the corresponding Chi2 accounting for covariances between the 2 clusters +/// return trackParamAtCluster2 + + // extrapolate linearly track parameters and covariances at the z position of the second cluster + trackParamAtCluster2 = trackParamAtCluster1; + AliMUONTrackExtrap::LinearExtrapToZ(&trackParamAtCluster2, cluster2->GetZ()); + + // set pointer to cluster2 into trackParamAtCluster2 + trackParamAtCluster2.SetClusterPtr(cluster2); + + // Set differences between track and clusters in the bending and non bending directions + AliMUONVCluster* cluster1 = trackParamAtCluster1.GetClusterPtr(); + Double_t dX1 = cluster1->GetX() - trackParamAtCluster1.GetNonBendingCoor(); + Double_t dX2 = cluster2->GetX() - trackParamAtCluster2.GetNonBendingCoor(); + Double_t dY1 = cluster1->GetY() - trackParamAtCluster1.GetBendingCoor(); + Double_t dY2 = cluster2->GetY() - trackParamAtCluster2.GetBendingCoor(); + + // Calculate errors and covariances + const TMatrixD& kParamCov1 = trackParamAtCluster1.GetCovariances(); + const TMatrixD& kParamCov2 = trackParamAtCluster2.GetCovariances(); + Double_t dZ = trackParamAtCluster2.GetZ() - trackParamAtCluster1.GetZ(); + Double_t sigma2X1 = kParamCov1(0,0) + cluster1->GetErrX2(); + Double_t sigma2X2 = kParamCov2(0,0) + cluster2->GetErrX2(); + Double_t covX1X2 = kParamCov1(0,0) + dZ * kParamCov1(0,1); + Double_t sigma2Y1 = kParamCov1(2,2) + cluster1->GetErrY2(); + Double_t sigma2Y2 = kParamCov2(2,2) + cluster2->GetErrY2(); + Double_t covY1Y2 = kParamCov1(2,2) + dZ * kParamCov1(2,3); + + // Compute chi2 + Double_t detX = sigma2X1 * sigma2X2 - covX1X2 * covX1X2; + Double_t detY = sigma2Y1 * sigma2Y2 - covY1Y2 * covY1Y2; + if (detX == 0. || detY == 0.) return 1.e10; + return (dX1 * dX1 * sigma2X2 + dX2 * dX2 * sigma2X1 - 2. * dX1 * dX2 * covX1X2) / detX + + (dY1 * dY1 * sigma2Y2 + dY2 * dY2 * sigma2Y1 - 2. * dY1 * dY2 * covY1Y2) / detY; + +} - trackParam = *((AliMUONTrackParam*) (track->GetTrackParamAtHit()->Last())); - AliMUONTrackExtrap::ExtrapToZ(&trackParam, AliMUONConstants::DefaultChamberZ(10)); // extrap to 1st trigger chamber - - xTrack = trackParam.GetNonBendingCoor(); - yTrack = trackParam.GetBendingCoor(); - ySlopeTrack = trackParam.GetBendingSlope(); - minChi2MatchTrigger = 999.; - - triggerTrack = (AliMUONTriggerTrack*) recTriggerTracks->First(); - while(triggerTrack){ - distTriggerTrack[0] = (triggerTrack->GetX11()-xTrack)/kDistSigma[0]; - distTriggerTrack[1] = (triggerTrack->GetY11()-yTrack)/kDistSigma[1]; - distTriggerTrack[2] = (TMath::Tan(triggerTrack->GetThetay())-ySlopeTrack)/kDistSigma[2]; - chi2 = 0.; - for (Int_t iVar = 0; iVar < 3; iVar++) chi2 += distTriggerTrack[iVar]*distTriggerTrack[iVar]; - chi2 /= 3.; // Normalized Chi2: 3 degrees of freedom (X,Y,slopeY) - if (chi2 < fMaxNormChi2MatchTrigger) { - Bool_t isDoubleTrack = (TMath::Abs(chi2 - minChi2MatchTrigger)<1.); - if (chi2 < minChi2MatchTrigger && chi2 < fMaxNormChi2MatchTrigger) { - if(isDoubleTrack){ - doubleMatch = loTrgNum; - doubleChi2 = chi2MatchTrigger; - } - minChi2MatchTrigger = chi2; - chi2MatchTrigger = chi2; - loTrgNum=triggerTrack->GetLoTrgNum(); - locTrg = (AliMUONLocalTrigger*)localTrigger->UncheckedAt(loTrgNum); - matchTrigger=0; - if(locTrg->LoLpt()>0)matchTrigger=1; - if(locTrg->LoHpt()>0)matchTrigger=2; + //__________________________________________________________________________ +Bool_t AliMUONVTrackReconstructor::FollowLinearTrackInChamber(AliMUONTrack &trackCandidate, const AliMUONVClusterStore& clusterStore, + Int_t nextChamber) +{ + /// Follow trackCandidate in chamber(0..) nextChamber assuming linear propagation, 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. + /// return kTRUE if new cluster(s) have been found (otherwise return kFALSE) + AliDebug(1,Form("Enter FollowLinearTrackInChamber(1..) %d", nextChamber+1)); + + Double_t chi2WithOneCluster = 1.e10; + Double_t maxChi2WithOneCluster = 2. * AliMUONReconstructor::GetRecoParam()->GetSigmaCutForTracking() * + AliMUONReconstructor::GetRecoParam()->GetSigmaCutForTracking(); // 2 because 2 quantities in chi2 + Double_t bestChi2WithOneCluster = maxChi2WithOneCluster; + Bool_t foundOneCluster = kFALSE; + AliMUONTrack *newTrack = 0x0; + AliMUONVCluster *cluster; + AliMUONTrackParam trackParam; + AliMUONTrackParam extrapTrackParamAtCluster; + AliMUONTrackParam bestTrackParamAtCluster; + + // Get track parameters according to the propagation direction + if (nextChamber > 7) trackParam = *(AliMUONTrackParam*)trackCandidate.GetTrackParamAtCluster()->Last(); + else trackParam = *(AliMUONTrackParam*)trackCandidate.GetTrackParamAtCluster()->First(); + + // Printout for debuging + if ((AliLog::GetDebugLevel("MUON","AliMUONVTrackReconstructor") >= 2) || (AliLog::GetGlobalDebugLevel() >= 2)) { + cout<= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) { + cout << "FollowLinearTrackInChamber: look for cluster in chamber(1..): " << nextChamber+1 << endl; + } + + // Create iterators to loop over clusters in chamber + TIter next(clusterStore.CreateChamberIterator(nextChamber,nextChamber)); + + // look for candidates in chamber + while ( ( cluster = static_cast(next()) ) ) { + + // try to add the current cluster fast + if (!TryOneClusterFast(trackParam, cluster)) continue; + + // try to add the current cluster accuratly + extrapTrackParamAtCluster = trackParam; + AliMUONTrackExtrap::LinearExtrapToZ(&extrapTrackParamAtCluster, cluster->GetZ()); + chi2WithOneCluster = TryOneCluster(extrapTrackParamAtCluster, cluster, extrapTrackParamAtCluster); + + // if good chi2 then consider to add cluster + if (chi2WithOneCluster < maxChi2WithOneCluster) { + foundOneCluster = kTRUE; + + // Printout for debuging + if ((AliLog::GetDebugLevel("MUON","AliMUONVTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) { + cout << "FollowLinearTrackInChamber: found one cluster in chamber(1..): " << nextChamber+1 + << " (Chi2 = " << chi2WithOneCluster << ")" << endl; + cluster->Print(); + } + + 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); + if (AliMUONReconstructor::GetRecoParam()->RequestStation(nextChamber/2)) + extrapTrackParamAtCluster.SetRemovable(kFALSE); + else extrapTrackParamAtCluster.SetRemovable(kTRUE); + newTrack->AddTrackParamAtCluster(extrapTrackParamAtCluster,*cluster); + fNRecTracks++; + + // Printout for debuging + if ((AliLog::GetDebugLevel("MUON","AliMUONVTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) { + cout << "FollowLinearTrackInChamber: added one cluster in chamber(1..): " << nextChamber+1 << endl; + if (AliLog::GetGlobalDebugLevel() >= 3) newTrack->RecursiveDump(); + } + + } else if (chi2WithOneCluster < bestChi2WithOneCluster) { + // keep track of the best cluster + bestChi2WithOneCluster = chi2WithOneCluster; + bestTrackParamAtCluster = extrapTrackParamAtCluster; + } + + } + + } + + // fill out the best track if required else clean up the fRecTracksPtr array + if (!AliMUONReconstructor::GetRecoParam()->TrackAllTracks()) { + if (foundOneCluster) { + if (AliMUONReconstructor::GetRecoParam()->RequestStation(nextChamber/2)) + bestTrackParamAtCluster.SetRemovable(kFALSE); + else bestTrackParamAtCluster.SetRemovable(kTRUE); + trackCandidate.AddTrackParamAtCluster(bestTrackParamAtCluster,*(bestTrackParamAtCluster.GetClusterPtr())); + + // Printout for debuging + if ((AliLog::GetDebugLevel("MUON","AliMUONVTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) { + cout << "FollowLinearTrackInChamber: 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 AliMUONVTrackReconstructor::FollowLinearTrackInStation(AliMUONTrack &trackCandidate, const AliMUONVClusterStore& clusterStore, + Int_t nextStation) +{ + /// Follow trackCandidate in station(0..) nextStation assuming linear propagation, 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. + /// return kTRUE if new cluster(s) have been found (otherwise return kFALSE) + AliDebug(1,Form("Enter FollowLinearTrackInStation(1..) %d", nextStation+1)); + + // 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; + AliMUONVCluster *clusterCh1, *clusterCh2; + AliMUONTrackParam trackParam; + 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 according to the propagation direction + if (nextStation==4) trackParam = *(AliMUONTrackParam*)trackCandidate.GetTrackParamAtCluster()->Last(); + else trackParam = *(AliMUONTrackParam*)trackCandidate.GetTrackParamAtCluster()->First(); + + // Printout for debuging + if ((AliLog::GetDebugLevel("MUON","AliMUONVTrackReconstructor") >= 2) || (AliLog::GetGlobalDebugLevel() >= 2)) { + cout<= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) { + cout << "FollowLinearTrackInStation: look for clusters in chamber(1..): " << ch2+1 << endl; + } + + // 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(trackParam, clusterCh2)) continue; + + // try to add the current cluster accuratly + extrapTrackParamAtCluster2 = trackParam; + AliMUONTrackExtrap::LinearExtrapToZ(&extrapTrackParamAtCluster2, clusterCh2->GetZ()); + chi2WithOneCluster = TryOneCluster(extrapTrackParamAtCluster2, 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","AliMUONVTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) { + cout << "FollowLinearTrackInStation: found one cluster in chamber(1..): " << ch2+1 + << " (Chi2 = " << chi2WithOneCluster << ")" << endl; + clusterCh2->Print(); + cout << " look for second clusters in chamber(1..): " << ch1+1 << " ..." << endl; + } + + // add MCS effect + AliMUONTrackExtrap::AddMCSEffect(&extrapTrackParamAtCluster2,AliMUONConstants::ChamberThicknessInX0(),1.); + + // 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(extrapTrackParamAtCluster2, clusterCh1)) continue; + + // try to add the current cluster in addition to the one found in the previous chamber + chi2WithTwoClusters = TryTwoClustersFast(extrapTrackParamAtCluster2, clusterCh1, extrapTrackParamAtCluster1); + + // if good chi2 then consider to add the 2 clusters to the "trackCandidate" + if (chi2WithTwoClusters < maxChi2WithTwoClusters) { + foundSecondCluster = kTRUE; + foundTwoClusters = kTRUE; + + // Printout for debuging + if ((AliLog::GetDebugLevel("MUON","AliMUONVTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) { + cout << "FollowLinearTrackInStation: found one cluster in chamber(1..): " << ch1+1 + << " (Chi2 = " << chi2WithTwoClusters << ")" << endl; + clusterCh1->Print(); } - else if(isDoubleTrack) { - doubleMatch = triggerTrack->GetLoTrgNum(); - doubleChi2 = chi2; + + 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); + extrapTrackParamAtCluster1.SetRemovable(kTRUE); + newTrack->AddTrackParamAtCluster(extrapTrackParamAtCluster1,*clusterCh1); + extrapTrackParamAtCluster2.SetRemovable(kTRUE); + newTrack->AddTrackParamAtCluster(extrapTrackParamAtCluster2,*clusterCh2); + fNRecTracks++; + + // Tag clusterCh1 as used + clusterCh1Used[iCluster1] = kTRUE; + + // Printout for debuging + if ((AliLog::GetDebugLevel("MUON","AliMUONVTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) { + cout << "FollowLinearTrackInStation: added two clusters in station(1..): " << nextStation+1 << endl; + if (AliLog::GetGlobalDebugLevel() >= 3) newTrack->RecursiveDump(); + } + + } else if (chi2WithTwoClusters < bestChi2WithTwoClusters) { + // keep track of the best couple of clusters + bestChi2WithTwoClusters = chi2WithTwoClusters; + bestTrackParamAtCluster1 = extrapTrackParamAtCluster1; + bestTrackParamAtCluster2 = extrapTrackParamAtCluster2; + } + + } + + } + + // if no cluster found in chamber1 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); + if (AliMUONReconstructor::GetRecoParam()->RequestStation(nextStation)) + extrapTrackParamAtCluster2.SetRemovable(kFALSE); + else extrapTrackParamAtCluster2.SetRemovable(kTRUE); + newTrack->AddTrackParamAtCluster(extrapTrackParamAtCluster2,*clusterCh2); + fNRecTracks++; + + // Printout for debuging + if ((AliLog::GetDebugLevel("MUON","AliMUONVTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) { + cout << "FollowLinearTrackInStation: added one cluster in chamber(1..): " << ch2+1 << endl; + if (AliLog::GetGlobalDebugLevel() >= 3) newTrack->RecursiveDump(); } + + } else if (!foundTwoClusters && chi2WithOneCluster < bestChi2WithOneCluster) { + // keep track of the best cluster except if a couple of clusters has already been found + bestChi2WithOneCluster = chi2WithOneCluster; + bestTrackParamAtCluster1 = extrapTrackParamAtCluster2; + } + } - triggerTrack = (AliMUONTriggerTrack*) recTriggerTracks->After(triggerTrack); + } - if(doubleMatch>=0){ // If two trigger tracks match, select the one passing more trigger cuts - AliDebug(1, Form("Two candidates found: %i and %i",loTrgNum,doubleMatch)); - AliMUONLocalTrigger *locTrg1 = (AliMUONLocalTrigger*)localTrigger->UncheckedAt(doubleMatch); - if((locTrg1->LoLpt()>0 && matchTrigger<1) || (locTrg1->LoHpt() && matchTrigger<2)){ - if(locTrg1->LoHpt()>0)matchTrigger=2; - else matchTrigger=1; - loTrgNum = doubleMatch; - chi2MatchTrigger=doubleChi2; - } + + } + + // 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 clusters has been found + if (AliMUONReconstructor::GetRecoParam()->TrackAllTracks() || !foundTwoClusters) { + + // Printout for debuging + if ((AliLog::GetDebugLevel("MUON","AliMUONVTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) { + cout << "FollowLinearTrackInStation: look for single cluster in chamber(1..): " << ch1+1 << endl; } - track->SetMatchTrigger(matchTrigger); - track->SetLoTrgNum(loTrgNum); - track->SetChi2MatchTrigger(chi2MatchTrigger); - track = (AliMUONTrack*) fRecTracksPtr->After(track); + //Extrapolate trackCandidate to chamber "ch2" + AliMUONTrackExtrap::LinearExtrapToZ(&trackParam, AliMUONConstants::DefaultChamberZ(ch2)); + + // add MCS effect for next step + AliMUONTrackExtrap::AddMCSEffect(&trackParam,AliMUONConstants::ChamberThicknessInX0(),1.); + + // 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 clusters already used + + // try to add the current cluster fast + if (!TryOneClusterFast(trackParam, clusterCh1)) continue; + + // try to add the current cluster accuratly + extrapTrackParamAtCluster1 = trackParam; + AliMUONTrackExtrap::LinearExtrapToZ(&extrapTrackParamAtCluster1, clusterCh1->GetZ()); + chi2WithOneCluster = TryOneCluster(extrapTrackParamAtCluster1, 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","AliMUONVTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) { + cout << "FollowLinearTrackInStation: found one cluster in chamber(1..): " << ch1+1 + << " (Chi2 = " << chi2WithOneCluster << ")" << endl; + clusterCh1->Print(); + } + + 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); + if (AliMUONReconstructor::GetRecoParam()->RequestStation(nextStation)) + extrapTrackParamAtCluster1.SetRemovable(kFALSE); + else extrapTrackParamAtCluster1.SetRemovable(kTRUE); + newTrack->AddTrackParamAtCluster(extrapTrackParamAtCluster1,*clusterCh1); + fNRecTracks++; + + // Printout for debuging + if ((AliLog::GetDebugLevel("MUON","AliMUONVTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) { + cout << "FollowLinearTrackInStation: added one cluster in chamber(1..): " << ch1+1 << endl; + if (AliLog::GetGlobalDebugLevel() >= 3) newTrack->RecursiveDump(); + } + + } else if (chi2WithOneCluster < bestChi2WithOneCluster) { + // keep track of the best 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 (!AliMUONReconstructor::GetRecoParam()->TrackAllTracks()) { + if (foundTwoClusters) { + bestTrackParamAtCluster1.SetRemovable(kTRUE); + trackCandidate.AddTrackParamAtCluster(bestTrackParamAtCluster1,*(bestTrackParamAtCluster1.GetClusterPtr())); + bestTrackParamAtCluster2.SetRemovable(kTRUE); + trackCandidate.AddTrackParamAtCluster(bestTrackParamAtCluster2,*(bestTrackParamAtCluster2.GetClusterPtr())); + + // Printout for debuging + if ((AliLog::GetDebugLevel("MUON","AliMUONVTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) { + cout << "FollowLinearTrackInStation: added the two best clusters in station(1..): " << nextStation+1 << endl; + if (AliLog::GetGlobalDebugLevel() >= 3) newTrack->RecursiveDump(); + } + + } else if (foundOneCluster) { + if (AliMUONReconstructor::GetRecoParam()->RequestStation(nextStation)) + bestTrackParamAtCluster1.SetRemovable(kFALSE); + else bestTrackParamAtCluster1.SetRemovable(kTRUE); + trackCandidate.AddTrackParamAtCluster(bestTrackParamAtCluster1,*(bestTrackParamAtCluster1.GetClusterPtr())); + + // Printout for debuging + if ((AliLog::GetDebugLevel("MUON","AliMUONVTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) { + cout << "FollowLinearTrackInStation: added the best cluster in chamber(1..): " << bestTrackParamAtCluster1.GetClusterPtr()->GetChamberId()+1 << endl; + if (AliLog::GetGlobalDebugLevel() >= 3) newTrack->RecursiveDump(); + } + + } else { + delete [] clusterCh1Used; + return kFALSE; + } + + } else if (foundOneCluster || foundTwoClusters) { + + // remove obsolete track + fRecTracksPtr->Remove(&trackCandidate); + fNRecTracks--; + + } else { + delete [] clusterCh1Used; + return kFALSE; } + + delete [] clusterCh1Used; + return kTRUE; + +} - return; +//__________________________________________________________________________ +void AliMUONVTrackReconstructor::ImproveTracks() +{ + /// 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"); + + AliMUONTrack *track, *nextTrack; + + // Remove double track to improve only "good" tracks + RemoveDoubleTracks(); + + track = (AliMUONTrack*) fRecTracksPtr->First(); + while (track) { + + // prepare next track in case the actual track is suppressed + nextTrack = (AliMUONTrack*) fRecTracksPtr->After(track); + + ImproveTrack(*track); + + // remove track if improvement failed + if (!track->IsImproved()) { + fRecTracksPtr->Remove(track); + fNRecTracks--; + } + + track = nextTrack; + } + + // compress the array in case of some tracks have been removed + fRecTracksPtr->Compress(); + } //__________________________________________________________________________ -void AliMUONVTrackReconstructor::EventReconstructTrigger(void) +void AliMUONVTrackReconstructor::Finalize() { - /// To reconstruct trigger for one event - AliDebug(1,"Enter EventReconstructTrigger"); - MakeTriggerTracks(); - return; + /// Recompute track parameters and covariances at each attached cluster from those at the first one + + AliMUONTrack *track; + + track = (AliMUONTrack*) fRecTracksPtr->First(); + while (track) { + + FinalizeTrack(*track); + + track = (AliMUONTrack*) fRecTracksPtr->After(track); + + } + +} + +//__________________________________________________________________________ +void AliMUONVTrackReconstructor::ValidateTracksWithTrigger(AliMUONVTrackStore& trackStore, + const AliMUONVTriggerTrackStore& triggerTrackStore, + const AliMUONVTriggerStore& triggerStore, + const AliMUONTrackHitPattern& trackHitPattern) +{ + /// Try to match track from tracking system with trigger track + AliCodeTimerAuto(""); + + trackHitPattern.ExecuteValidation(trackStore, triggerTrackStore, triggerStore); } //__________________________________________________________________________ -Bool_t AliMUONVTrackReconstructor::MakeTriggerTracks(void) +void AliMUONVTrackReconstructor::EventReconstructTrigger(const AliMUONTriggerCircuit& circuit, + const AliMUONVTriggerStore& triggerStore, + AliMUONVTriggerTrackStore& triggerTrackStore) { /// To make the trigger tracks from Local Trigger - AliDebug(1, "Enter MakeTriggerTracks"); - - TTree* treeR; - UChar_t gloTrigPat; - TClonesArray *localTrigger; - TClonesArray *globalTrigger; - AliMUONLocalTrigger *locTrg; - AliMUONGlobalTrigger *gloTrg; + AliDebug(1, ""); + AliCodeTimerAuto(""); + + AliMUONGlobalTrigger* globalTrigger = triggerStore.Global(); + + UChar_t gloTrigPat = 0; - treeR = fMUONData->TreeR(); - if (!treeR) { - AliWarning("TreeR is not loaded"); - return kFALSE; + if (globalTrigger) + { + gloTrigPat = globalTrigger->GetGlobalResponse(); } - fMUONData->SetTreeAddress("TC"); - fMUONData->GetTrigger(); - - // global trigger for trigger pattern - gloTrigPat = 0; - globalTrigger = fMUONData->GlobalTrigger(); - gloTrg = (AliMUONGlobalTrigger*)globalTrigger->UncheckedAt(0); - - if (gloTrg) - gloTrigPat = gloTrg->GetGlobalResponse(); - - - // local trigger for tracking - localTrigger = fMUONData->LocalTrigger(); - Int_t nlocals = (Int_t) (localTrigger->GetEntries()); + TIter next(triggerStore.CreateIterator()); + AliMUONLocalTrigger* locTrg(0x0); Float_t z11 = AliMUONConstants::DefaultChamberZ(10); Float_t z21 = AliMUONConstants::DefaultChamberZ(12); - - Float_t y11 = 0.; - Int_t stripX21 = 0; - Float_t y21 = 0.; - Float_t x11 = 0.; - - for (Int_t i=0; iUncheckedAt(i); - Bool_t xTrig=kFALSE; - Bool_t yTrig=kFALSE; + AliMUONTriggerTrack triggerTrack; + + while ( ( locTrg = static_cast(next()) ) ) + { + Bool_t xTrig=kFALSE; + Bool_t yTrig=kFALSE; + + Int_t localBoardId = locTrg->LoCircuit(); + if ( locTrg->LoSdev()==1 && locTrg->LoDev()==0 && + locTrg->LoStripX()==0) xTrig=kFALSE; // no trigger in X + else xTrig=kTRUE; // trigger in X + if (locTrg->LoTrigY()==1 && + locTrg->LoStripY()==15 ) yTrig = kFALSE; // no trigger in Y + else yTrig = kTRUE; // trigger in Y + + if (xTrig && yTrig) + { // make Trigger Track if trigger in X and Y - if ( locTrg->LoSdev()==1 && locTrg->LoDev()==0 && - locTrg->LoStripX()==0) xTrig=kFALSE; // no trigger in X - else xTrig=kTRUE; // trigger in X - if (locTrg->LoTrigY()==1 && - locTrg->LoStripY()==15 ) yTrig = kFALSE; // no trigger in Y - else yTrig = kTRUE; // trigger in Y - - if (xTrig && yTrig) { // make Trigger Track if trigger in X and Y - - AliDebug(1, "AliMUONTrackReconstructor::MakeTriggerTrack using NEW trigger \n"); - AliMUONTriggerCircuit* circuit = - (AliMUONTriggerCircuit*)fTriggerCircuit->At(locTrg->LoCircuit()-1); // -1 !!! - - y11 = circuit->GetY11Pos(locTrg->LoStripX()); -// need first to convert deviation to [0-30] -// (see AliMUONLocalTriggerBoard::LocalTrigger) - Int_t deviation = locTrg->LoDev(); - Int_t sign = 0; - if ( !locTrg->LoSdev() && deviation ) sign=-1; - if ( !locTrg->LoSdev() && !deviation ) sign= 0; - if ( locTrg->LoSdev() == 1 ) sign=+1; - deviation *= sign; - deviation += 15; - stripX21 = locTrg->LoStripX()+deviation+1; - y21 = circuit->GetY21Pos(stripX21); - x11 = circuit->GetX11Pos(locTrg->LoStripY()); - - AliDebug(1, Form(" MakeTriggerTrack %d %d %d %d %d %f %f %f \n",i,locTrg->LoCircuit(), - locTrg->LoStripX(),locTrg->LoStripX()+locTrg->LoDev()+1,locTrg->LoStripY(),y11, y21, x11)); - - Float_t thetax = TMath::ATan2( x11 , z11 ); - Float_t thetay = TMath::ATan2( (y21-y11) , (z21-z11) ); - - fTriggerTrack->SetX11(x11); - fTriggerTrack->SetY11(y11); - fTriggerTrack->SetThetax(thetax); - fTriggerTrack->SetThetay(thetay); - fTriggerTrack->SetGTPattern(gloTrigPat); - fTriggerTrack->SetLoTrgNum(i); - - fMUONData->AddRecTriggerTrack(*fTriggerTrack); - } // board is fired + Float_t y11 = circuit.GetY11Pos(localBoardId, locTrg->LoStripX()); + // need first to convert deviation to [0-30] + // (see AliMUONLocalTriggerBoard::LocalTrigger) + Int_t deviation = locTrg->LoDev(); + Int_t sign = 0; + if ( !locTrg->LoSdev() && deviation ) sign=-1; + if ( !locTrg->LoSdev() && !deviation ) sign= 0; + if ( locTrg->LoSdev() == 1 ) sign=+1; + deviation *= sign; + deviation += 15; + Int_t stripX21 = locTrg->LoStripX()+deviation+1; + Float_t y21 = circuit.GetY21Pos(localBoardId, stripX21); + Float_t x11 = circuit.GetX11Pos(localBoardId, locTrg->LoStripY()); + + AliDebug(1, Form(" MakeTriggerTrack %d %d %d %d %f %f %f \n",locTrg->LoCircuit(), + locTrg->LoStripX(),locTrg->LoStripX()+locTrg->LoDev()+1,locTrg->LoStripY(),y11, y21, x11)); + + Float_t thetax = TMath::ATan2( x11 , z11 ); + Float_t thetay = TMath::ATan2( (y21-y11) , (z21-z11) ); + + triggerTrack.SetX11(x11); + triggerTrack.SetY11(y11); + triggerTrack.SetThetax(thetax); + triggerTrack.SetThetay(thetay); + triggerTrack.SetGTPattern(gloTrigPat); + triggerTrack.SetLoTrgNum(localBoardId); + + triggerTrackStore.Add(triggerTrack); + } // board is fired } // end of loop on Local Trigger - - return kTRUE; -} - -//__________________________________________________________________________ -void AliMUONVTrackReconstructor::EventDumpTrigger(void) -{ - /// Dump reconstructed trigger event - /// and the particle parameters - AliMUONTriggerTrack *triggertrack ; - Int_t nTriggerTracks = fMUONData->RecTriggerTracks()->GetEntriesFast(); - - AliDebug(1, "****** enter EventDumpTrigger ******"); - AliDebug(1, Form("Number of Reconstructed tracks : %d ", nTriggerTracks)); - - // Loop over reconstructed tracks - for (Int_t trackIndex = 0; trackIndex < nTriggerTracks; trackIndex++) { - triggertrack = (AliMUONTriggerTrack*)fMUONData->RecTriggerTracks()->At(trackIndex); - printf(" trigger track number %i x11=%f y11=%f thetax=%f thetay=%f \n", - trackIndex, - triggertrack->GetX11(),triggertrack->GetY11(), - triggertrack->GetThetax(),triggertrack->GetThetay()); - } }