X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=MUON%2FAliMUONVTrackReconstructor.cxx;h=21edde3888f07d0e7f03fb250b7ab41d02a18e32;hb=9570f0874a46c860c2c8f797f105aade49537c27;hp=02c956aefffcecee07aab1bc8ade57e01d21f626;hpb=e1a10d4115aae4682293b6dd274a3b765c06f415;p=u%2Fmrichter%2FAliRoot.git diff --git a/MUON/AliMUONVTrackReconstructor.cxx b/MUON/AliMUONVTrackReconstructor.cxx index 02c956aefff..21edde3888f 100644 --- a/MUON/AliMUONVTrackReconstructor.cxx +++ b/MUON/AliMUONVTrackReconstructor.cxx @@ -15,543 +15,1077 @@ /* $Id$ */ -//////////////////////////////////// -// -// 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 -// -// It contains as methods, among others: -// * EventReconstruct to build the muon tracks -// * EventReconstructTrigger to build the trigger tracks -// -//////////////////////////////////// - -#include -#include -#include +//----------------------------------------------------------------------------- +/// \class AliMUONVTrackReconstructor +/// Virtual MUON track reconstructor in ALICE (class renamed from AliMUONEventReconstructor) +/// +/// 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 +/// +/// Several options and adjustable parameters are available for both KALMAN and ORIGINAL +/// tracking algorithms. They can be changed through the AliMUONRecoParam object +/// set in the reconstruction macro or read from the CDB +/// (see methods in AliMUONRecoParam.h file for details) +/// +/// 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" #include "AliMUONLocalTrigger.h" #include "AliMUONGlobalTrigger.h" -#include "AliMUONSegment.h" #include "AliMUONTrack.h" -#include "AliMagF.h" +#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 "AliMUONRecoParam.h" + +#include "AliMpDEManager.h" +#include "AliMpArea.h" + #include "AliLog.h" +#include "AliCodeTimer.h" #include "AliTracker.h" -ClassImp(AliMUONVTrackReconstructor) // Class implementation in ROOT context +#include +#include +#include +#include -//************* 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::fgkDefaultMaxSigma2Distance = 16.0; -// Simple magnetic field: -// Value taken from macro MUONtracking.C: 0.7 T, hence 7 kG -// Length and Position from reco_muon.F, with opposite sign: -// Length = ZMAGEND-ZCOIL -// Position = (ZMAGEND+ZCOIL)/2 -// to be ajusted differently from real magnetic field ???? -const Double_t AliMUONVTrackReconstructor::fgkDefaultSimpleBValue = 7.0; -const Double_t AliMUONVTrackReconstructor::fgkDefaultSimpleBLength = 428.0; -const Double_t AliMUONVTrackReconstructor::fgkDefaultSimpleBPosition = 1019.0; +#include -//__________________________________________________________________________ -AliMUONVTrackReconstructor::AliMUONVTrackReconstructor(AliMUONData* data) - : TObject(), - fMinBendingMomentum(fgkDefaultMinBendingMomentum), - fMaxBendingMomentum(fgkDefaultMaxBendingMomentum), - fBendingResolution(fgkDefaultBendingResolution), - fNonBendingResolution(fgkDefaultNonBendingResolution), - fMaxSigma2Distance(fgkDefaultMaxSigma2Distance), - fChamberThicknessInX0(AliMUONConstants::DefaultChamberThicknessInX0()), - fSimpleBValue(fgkDefaultSimpleBValue), - fSimpleBLength(fgkDefaultSimpleBLength), - fSimpleBPosition(fgkDefaultSimpleBPosition), - fSegmentMaxDistBending(0x0), - fSegmentMaxDistNonBending(0x0), - fHitsForRecPtr(0x0), - fNHitsForRec(0), - fNHitsForRecPerChamber(0x0), - fIndexOfFirstHitForRecPerChamber(0x0), - fSegmentsPtr(0x0), - fNSegments(0x0), - fRecTracksPtr(0x0), - fNRecTracks(0), - fMUONData(data), - fTriggerTrack(new AliMUONTriggerTrack()), - fTriggerCircuit(0x0) +/// \cond CLASSIMP +ClassImp(AliMUONVTrackReconstructor) // Class implementation in ROOT context +/// \endcond + + //__________________________________________________________________________ +AliMUONVTrackReconstructor::AliMUONVTrackReconstructor(const AliMUONRecoParam* recoParam, + AliMUONVClusterServer* clusterServer) +: TObject(), +fRecTracksPtr(0x0), +fNRecTracks(0), +fClusterServer(clusterServer), +fkRecoParam(recoParam) { /// Constructor for class AliMUONVTrackReconstructor - fSegmentMaxDistBending = new Double_t[AliMUONConstants::NTrackingSt()]; - fSegmentMaxDistNonBending = new Double_t[AliMUONConstants::NTrackingSt()]; - fNHitsForRecPerChamber = new Int_t[AliMUONConstants::NTrackingCh()]; - fIndexOfFirstHitForRecPerChamber = new Int_t[AliMUONConstants::NTrackingCh()]; - - SetReconstructionParametersToDefaults(); - - // Memory allocation for the TClonesArray of hits for reconstruction - // Is 10000 the right size ???? - fHitsForRecPtr = new TClonesArray("AliMUONHitForRec", 10000); - - // Memory allocation for the TClonesArray's of segments in stations - // Is 2000 the right size ???? - fSegmentsPtr = new TClonesArray*[AliMUONConstants::NTrackingSt()]; - fNSegments = new Int_t[AliMUONConstants::NTrackingSt()]; - for (Int_t st = 0; st < AliMUONConstants::NTrackingSt(); st++) { - fSegmentsPtr[st] = new TClonesArray("AliMUONSegment", 2000); - fNSegments[st] = 0; // really needed or GetEntriesFast sufficient ???? - } - - const AliMagF* kField = AliTracker::GetFieldMap(); - if (!kField) AliFatal("No field available"); - // Sign of fSimpleBValue according to sign of Bx value at (50,50,-950). - Float_t b[3], x[3]; - x[0] = 50.; x[1] = 50.; x[2] = -950.; - kField->Field(x,b); - - fSimpleBValue = TMath::Sign(fSimpleBValue,(Double_t) b[0]); - fSimpleBPosition = TMath::Sign(fSimpleBPosition,(Double_t) x[2]); - // See how to get fSimple(BValue, BLength, BPosition) - // automatically calculated from the actual magnetic field ???? + /// WARNING: if clusterServer=0x0, no clusterization will be possible at this level + + // Memory allocation for the TClonesArray of reconstructed tracks + fRecTracksPtr = new TClonesArray("AliMUONTrack", 100); + + // set the magnetic field for track extrapolations + AliMUONTrackExtrap::SetField(); } //__________________________________________________________________________ -AliMUONVTrackReconstructor::~AliMUONVTrackReconstructor(void) +AliMUONVTrackReconstructor::~AliMUONVTrackReconstructor() { /// Destructor for class AliMUONVTrackReconstructor - delete [] fSegmentMaxDistBending; - delete [] fSegmentMaxDistNonBending; - delete [] fNHitsForRecPerChamber; - delete [] fIndexOfFirstHitForRecPerChamber; - delete fTriggerTrack; - delete fHitsForRecPtr; - // Correct destruction of everything ???? - for (Int_t st = 0; st < AliMUONConstants::NTrackingSt(); st++) delete fSegmentsPtr[st]; - delete [] fSegmentsPtr; - delete [] fNSegments; + delete fRecTracksPtr; } //__________________________________________________________________________ -void AliMUONVTrackReconstructor::SetReconstructionParametersToDefaults(void) +void AliMUONVTrackReconstructor::ResetTracks() { - /// Set reconstruction parameters for making segments to default values - // Would be much more convenient with a structure (or class) ???? - - // ******** Parameters for making segments - // should be parametrized ???? - // according to interval between chambers in a station ???? - // Maximum distance in non bending plane - // 5 * 0.22 just to remember the way it was made in TRACKF_STAT - // SIGCUT*DYMAX(IZ) - for (Int_t st = 0; st < AliMUONConstants::NTrackingCh()/2; st++) - fSegmentMaxDistNonBending[st] = 5. * 0.22; - // Maximum distance in bending plane: - // values from TRACKF_STAT, corresponding to (J psi 20cm), - // scaled to the real distance between chambers in a station - fSegmentMaxDistBending[0] = TMath::Abs( 1.5 * - (AliMUONConstants::DefaultChamberZ(1) - AliMUONConstants::DefaultChamberZ(0)) / 20.0); - fSegmentMaxDistBending[1] = TMath::Abs( 1.5 * - (AliMUONConstants::DefaultChamberZ(3) - AliMUONConstants::DefaultChamberZ(2)) / 20.0); - fSegmentMaxDistBending[2] = TMath::Abs( 3.0 * - (AliMUONConstants::DefaultChamberZ(5) - AliMUONConstants::DefaultChamberZ(4)) / 20.0); - fSegmentMaxDistBending[3] = TMath::Abs( 6.0 * - (AliMUONConstants::DefaultChamberZ(7) - AliMUONConstants::DefaultChamberZ(6)) / 20.0); - fSegmentMaxDistBending[4] = TMath::Abs( 6.0 * - (AliMUONConstants::DefaultChamberZ(9) - AliMUONConstants::DefaultChamberZ(8)) / 20.0); - + /// To reset the TClonesArray of reconstructed tracks + if (fRecTracksPtr) fRecTracksPtr->Clear("C"); + fNRecTracks = 0; return; } //__________________________________________________________________________ -void AliMUONVTrackReconstructor::EventReconstruct(void) +void AliMUONVTrackReconstructor::EventReconstruct(AliMUONVClusterStore& clusterStore, AliMUONVTrackStore& trackStore) { - // To reconstruct one event - AliDebug(1,"Enter EventReconstruct"); - ResetTracks(); //AZ - ResetSegments(); //AZ - ResetHitsForRec(); //AZ - AddHitsForRecFromRawClusters(); - MakeSegments(); - MakeTracks(); - if (fMUONData->IsTriggerTrackBranchesInTree()) - ValidateTracksWithTrigger(); + /// 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 (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 (GetRecoParam()->ComplementTracks()) ComplementTracks(clusterStore); + + // Improve the reconstructed tracks + if (GetRecoParam()->ImproveTracks()) ImproveTracks(); + + // Remove connected tracks + if (GetRecoParam()->RemoveConnectedTracksInSt12()) RemoveConnectedTracks(kFALSE); + else RemoveConnectedTracks(kTRUE); + + // 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()); + + // check if non bending slope is within tolerances + if (TMath::Abs(nonBendingSlope) > GetRecoParam()->GetMaxNonBendingSlope()) continue; + + // bending slope + bendingSlope = (cluster1->GetY() - cluster2->GetY()) / (cluster1->GetZ() - cluster2->GetZ()); + + // check the bending momentum of the bending slope depending if the field is ON or OFF + if (AliMUONTrackExtrap::IsFieldON()) { + + // impact parameter + impactParam = cluster1->GetY() - cluster1->GetZ() * bendingSlope; + + // absolute value of bending momentum + bendingMomentum = TMath::Abs(AliMUONTrackExtrap::GetBendingMomentumFromImpactParam(impactParam)); + + // check if bending momentum is within tolerances + if (bendingMomentum < GetRecoParam()->GetMinBendingMomentum() || + bendingMomentum > GetRecoParam()->GetMaxBendingMomentum()) continue; + + } else { + + // check if non bending slope is within tolerances + if (TMath::Abs(bendingSlope) > GetRecoParam()->GetMaxBendingSlope()) continue; + + } + // 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::ResetSegments(void) +void AliMUONVTrackReconstructor::RemoveIdenticalTracks() { - /// To reset the TClonesArray of segments and the number of Segments for all stations - for (Int_t st = 0; st < AliMUONConstants::NTrackingCh()/2; st++) { - if (fSegmentsPtr[st]) fSegmentsPtr[st]->Clear(); - fNSegments[st] = 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::ResetHitsForRec(void) +void AliMUONVTrackReconstructor::RemoveDoubleTracks() { - /// 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 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; } //__________________________________________________________________________ -void AliMUONVTrackReconstructor::SortHitsForRecWithIncreasingChamber() +void AliMUONVTrackReconstructor::RemoveConnectedTracks(Bool_t inSt345) { - /// 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 they share 1 cluster or more. + /// If inSt345=kTRUE only stations 3, 4 and 5 are considered. + /// 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, inSt345); + // check for identical tracks + if (clustersInCommon > 0) { + // 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; } //__________________________________________________________________________ -void AliMUONVTrackReconstructor::MakeSegmentsPerStation(Int_t Station) +void AliMUONVTrackReconstructor::AskForNewClustersInChamber(const AliMUONTrackParam &trackParam, + AliMUONVClusterStore& clusterStore, Int_t chamber) { - /// To make the list of segments in station number "Station" (0...) - /// from the list of hits to be reconstructed. - /// Updates "fNSegments"[Station]. - /// Segments in stations 4 and 5 are sorted - /// according to increasing absolute value of "impact parameter" - AliMUONHitForRec *hit1Ptr, *hit2Ptr; - AliMUONSegment *segment; - Bool_t last2st; - Double_t bendingSlope, distBend, distNonBend, extBendCoor, extNonBendCoor, - impactParam = 0., maxImpactParam = 0., minImpactParam = 0.; // =0 to avoid compilation warnings. - AliDebug(1,Form("Enter MakeSegmentsPerStation (0...) %d",Station)); - // first and second chambers (0...) in the station - Int_t ch1 = 2 * Station; - Int_t ch2 = ch1 + 1; - // variable true for stations downstream of the dipole: - // Station(0..4) equal to 3 or 4 - if ((Station == 3) || (Station == 4)) { - last2st = kTRUE; - // maximum impact parameter (cm) according to fMinBendingMomentum... - maxImpactParam = - TMath::Abs(GetImpactParamFromBendingMomentum(fMinBendingMomentum)); - // minimum impact parameter (cm) according to fMaxBendingMomentum... - minImpactParam = - TMath::Abs(GetImpactParamFromBendingMomentum(fMaxBendingMomentum)); - } - else last2st = kFALSE; - // extrapolation factor from Z of first chamber to Z of second chamber - // dZ to be changed to take into account fine structure of chambers ???? - Double_t extrapFact; - // index for current segment - Int_t segmentIndex = 0; - // Loop over HitsForRec 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]); - // extrapolation, - // on the straight line joining the HitForRec to the vertex (0,0,0), - // to the Z of the second chamber of the station - // Loop over HitsForRec 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]); - // absolute values of distances, in bending and non bending planes, - // between the HitForRec in the second chamber - // and the previous extrapolation - extrapFact = hit2Ptr->GetZ()/ hit1Ptr->GetZ(); - extBendCoor = extrapFact * hit1Ptr->GetBendingCoor(); - extNonBendCoor = extrapFact * hit1Ptr->GetNonBendingCoor(); - distBend = TMath::Abs(hit2Ptr->GetBendingCoor() - extBendCoor); - distNonBend = TMath::Abs(hit2Ptr->GetNonBendingCoor() - extNonBendCoor); - if (last2st) { - // bending slope - if ( hit1Ptr->GetZ() - hit2Ptr->GetZ() != 0.0 ) { - bendingSlope = (hit1Ptr->GetBendingCoor() - hit2Ptr->GetBendingCoor()) / - (hit1Ptr->GetZ() - hit2Ptr->GetZ()); - // absolute value of impact parameter - impactParam = - TMath::Abs(hit1Ptr->GetBendingCoor() - hit1Ptr->GetZ() * bendingSlope); - } - else { - AliWarning("hit1Ptr->GetZ() = hit2Ptr->GetZ(): impactParam set to maxImpactParam"); - impactParam = maxImpactParam; - } - } - // check for distances not too large, - // and impact parameter not too big if stations downstream of the dipole. - // Conditions "distBend" and "impactParam" correlated for these stations ???? - if ((distBend < fSegmentMaxDistBending[Station]) && - (distNonBend < fSegmentMaxDistNonBending[Station]) && - (!last2st || (impactParam < maxImpactParam)) && - (!last2st || (impactParam > minImpactParam))) { - // make new segment - segment = new ((*fSegmentsPtr[Station])[segmentIndex]) - AliMUONSegment(hit1Ptr, hit2Ptr); - // update "link" to this segment from the hit in the first chamber - if (hit1Ptr->GetNSegments() == 0) - hit1Ptr->SetIndexOfFirstSegment(segmentIndex); - hit1Ptr->SetNSegments(hit1Ptr->GetNSegments() + 1); - if (AliLog::GetGlobalDebugLevel() > 1) { - cout << "segmentIndex(0...): " << segmentIndex - << " distBend: " << distBend - << " distNonBend: " << distNonBend - << endl; - segment->Dump(); - cout << "HitForRec in first chamber" << endl; - hit1Ptr->Dump(); - cout << "HitForRec in second chamber" << endl; - hit2Ptr->Dump(); - }; - // increment index for current segment - segmentIndex++; - } - } //for (Int_t hit2 - } // for (Int_t hit1... - fNSegments[Station] = segmentIndex; - // Sorting according to "impact parameter" if station(1..5) 4 or 5, - // i.e. Station(0..4) 3 or 4, using the function "Compare". - // After this sorting, it is impossible to use - // the "fNSegments" and "fIndexOfFirstSegment" - // of the HitForRec in the first chamber to explore all segments formed with it. - // Is this sorting really needed ???? - if ((Station == 3) || (Station == 4)) (fSegmentsPtr[Station])->Sort(); - AliDebug(1,Form("Station: %d NSegments: %d ", Station, fNSegments[Station])); - return; + /// Ask the clustering to reconstruct new clusters around the track candidate position + + // check if the current chamber is useable + if (!fClusterServer || !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 kMaxDZ = 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) + kMaxDZ * kMaxDZ * kParamCov(1,1) + 2. * kMaxDZ * TMath::Abs(kParamCov(0,1)); + Double_t errY2 = kParamCov(2,2) + kMaxDZ * kMaxDZ * kParamCov(3,3) + 2. * kMaxDZ * TMath::Abs(kParamCov(2,3)); + Double_t dX = TMath::Abs(trackParam.GetNonBendingSlope()) * kMaxDZ + + GetRecoParam()->GetMaxNonBendingDistanceToTrack() + + GetRecoParam()->GetSigmaCutForTracking() * TMath::Sqrt(errX2); + Double_t dY = TMath::Abs(trackParam.GetBendingSlope()) * kMaxDZ + + GetRecoParam()->GetMaxBendingDistanceToTrack() + + 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, GetRecoParam()); + } //__________________________________________________________________________ -Double_t AliMUONVTrackReconstructor::GetImpactParamFromBendingMomentum(Double_t BendingMomentum) const +void AliMUONVTrackReconstructor::AskForNewClustersInStation(const AliMUONTrackParam &trackParam, + AliMUONVClusterStore& clusterStore, Int_t station) { - /// Returns impact parameter at vertex in bending plane (cm), - /// from the signed bending momentum "BendingMomentum" in bending plane (GeV/c), - /// using simple values for dipole magnetic field. - /// The sign of "BendingMomentum" is the sign of the charge. - return (-0.0003 * fSimpleBValue * fSimpleBLength * fSimpleBPosition / - BendingMomentum); + /// 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::GetBendingMomentumFromImpactParam(Double_t ImpactParam) const +Double_t AliMUONVTrackReconstructor::TryOneCluster(const AliMUONTrackParam &trackParam, AliMUONVCluster* cluster, + AliMUONTrackParam &trackParamAtCluster, Bool_t updatePropagator) { - /// Returns signed bending momentum in bending plane (GeV/c), - /// the sign being the sign of the charge for particles moving forward in Z, - /// from the impact parameter "ImpactParam" at vertex in bending plane (cm), - /// using simple values for dipole magnetic field. - return (-0.0003 * fSimpleBValue * fSimpleBLength * fSimpleBPosition / - ImpactParam); +/// Test the compatibility between the track and the cluster (using trackParam's covariance matrix): +/// return the corresponding Chi2 +/// return trackParamAtCluster + + // 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(); + + // Compute chi2 + return dX * dX / sigmaX2 + dY * dY / sigmaY2; + } //__________________________________________________________________________ -void AliMUONVTrackReconstructor::ValidateTracksWithTrigger(void) +Bool_t AliMUONVTrackReconstructor::TryOneClusterFast(const AliMUONTrackParam &trackParam, const AliMUONVCluster* cluster) { - /// Try to match track from tracking system with trigger track - AliMUONTrack *track; - AliMUONTrackParam trackParam; - AliMUONTriggerTrack *triggerTrack; +/// 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 - fMUONData->SetTreeAddress("RL"); - fMUONData->GetRecTriggerTracks(); - TClonesArray *recTriggerTracks = fMUONData->RecTriggerTracks(); + 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); + + Double_t dXmax = GetRecoParam()->GetSigmaCutForTracking() * TMath::Sqrt(errX2) + + GetRecoParam()->GetMaxNonBendingDistanceToTrack(); + Double_t dYmax = GetRecoParam()->GetSigmaCutForTracking() * TMath::Sqrt(errY2) + + GetRecoParam()->GetMaxBendingDistanceToTrack(); - Bool_t matchTrigger; - Double_t distSigma[3]={1,1,0.02}; // sigma of distributions (trigger-track) X,Y,slopeY - Double_t distTriggerTrack[3]; - Double_t chi2MatchTrigger, xTrack, yTrack, ySlopeTrack, dTrigTrackMin2, dTrigTrack2 = 0.; + if (TMath::Abs(dX) > dXmax || TMath::Abs(dY) > dYmax) return kFALSE; - track = (AliMUONTrack*) fRecTracksPtr->First(); - while (track) { - matchTrigger = kFALSE; - chi2MatchTrigger = 0.; - - trackParam = *((AliMUONTrackParam*) (track->GetTrackParamAtHit()->Last())); - trackParam.ExtrapToZ(AliMUONConstants::DefaultChamberZ(10)); // extrap to 1st trigger chamber - - xTrack = trackParam.GetNonBendingCoor(); - yTrack = trackParam.GetBendingCoor(); - ySlopeTrack = trackParam.GetBendingSlope(); - dTrigTrackMin2 = 999.; - - triggerTrack = (AliMUONTriggerTrack*) recTriggerTracks->First(); - while(triggerTrack){ - distTriggerTrack[0] = (triggerTrack->GetX11()-xTrack)/distSigma[0]; - distTriggerTrack[1] = (triggerTrack->GetY11()-yTrack)/distSigma[1]; - distTriggerTrack[2] = (TMath::Tan(triggerTrack->GetThetay())-ySlopeTrack)/distSigma[2]; - dTrigTrack2 = 0.; - for (Int_t iVar = 0; iVar < 3; iVar++) dTrigTrack2 += distTriggerTrack[iVar]*distTriggerTrack[iVar]; - if (dTrigTrack2 < dTrigTrackMin2 && dTrigTrack2 < fMaxSigma2Distance) { - dTrigTrackMin2 = dTrigTrack2; - matchTrigger = kTRUE; - chi2MatchTrigger = dTrigTrack2/3.; // Normalized Chi2, 3 variables (X,Y,slopeY) + 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; + +} + + //__________________________________________________________________________ +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. * GetRecoParam()->GetSigmaCutForTracking() * + 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 (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 (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; } - triggerTrack = (AliMUONTriggerTrack*) recTriggerTracks->After(triggerTrack); + } - track->SetMatchTrigger(matchTrigger); - track->SetChi2MatchTrigger(chi2MatchTrigger); - - track = (AliMUONTrack*) fRecTracksPtr->After(track); } - - return; + + // fill out the best track if required else clean up the fRecTracksPtr array + if (!GetRecoParam()->TrackAllTracks()) { + if (foundOneCluster) { + if (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; + } //__________________________________________________________________________ -void AliMUONVTrackReconstructor::EventReconstructTrigger(void) +Bool_t AliMUONVTrackReconstructor::FollowLinearTrackInStation(AliMUONTrack &trackCandidate, const AliMUONVClusterStore& clusterStore, + Int_t nextStation) { - /// To reconstruct trigger for one event - AliDebug(1,"Enter EventReconstructTrigger"); - MakeTriggerTracks(); - return; + /// 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. * GetRecoParam()->GetSigmaCutForTracking() * + GetRecoParam()->GetSigmaCutForTracking(); // 2 because 2 quantities in chi2 + Double_t maxChi2WithTwoClusters = 4. * GetRecoParam()->GetSigmaCutForTracking() * + GetRecoParam()->GetSigmaCutForTracking(); // 4 because 4 quantities in chi2 + Double_t bestChi2WithOneCluster = maxChi2WithOneCluster; + Double_t bestChi2WithTwoClusters = maxChi2WithTwoClusters; + Bool_t foundOneCluster = kFALSE; + Bool_t foundTwoClusters = kFALSE; + AliMUONTrack *newTrack = 0x0; + 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(); + } + + if (GetRecoParam()->TrackAllTracks()) { + // copy trackCandidate into a new track put at the end of fRecTracksPtr and add the new clusters + newTrack = new ((*fRecTracksPtr)[fRecTracksPtr->GetLast()+1]) AliMUONTrack(trackCandidate); + 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 (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 (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; + } + + } + + } + + } + + // 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 (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; + } + + //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 (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 (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 (!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 (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; + } - //__________________________________________________________________________ -Bool_t AliMUONVTrackReconstructor::MakeTriggerTracks(void) +//__________________________________________________________________________ +void AliMUONVTrackReconstructor::ImproveTracks() { - // To make the trigger tracks from Local Trigger - AliDebug(1, "Enter MakeTriggerTracks"); - - Int_t nTRentries; - UChar_t gloTrigPat; - TClonesArray *localTrigger; - TClonesArray *globalTrigger; - AliMUONLocalTrigger *locTrg; - AliMUONGlobalTrigger *gloTrg; - - TTree* treeR = fMUONData->TreeR(); - - nTRentries = Int_t(treeR->GetEntries()); - - treeR->GetEvent(0); // only one entry - - if (!(fMUONData->IsTriggerBranchesInTree())) { - AliWarning(Form("Trigger information is not avalaible, nTRentries = %d not equal to 1",nTRentries)); - return kFALSE; + /// 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; + + 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--; } - - 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(); + + track = nextTrack; + } + + // compress the array in case of some tracks have been removed + fRecTracksPtr->Compress(); +} - // local trigger for tracking - localTrigger = fMUONData->LocalTrigger(); - Int_t nlocals = (Int_t) (localTrigger->GetEntries()); +//__________________________________________________________________________ +void AliMUONVTrackReconstructor::Finalize() +{ + /// 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); + + } + +} - Float_t z11 = AliMUONConstants::DefaultChamberZ(10); - Float_t z21 = AliMUONConstants::DefaultChamberZ(12); +//__________________________________________________________________________ +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(""); - Float_t y11 = 0.; - Int_t stripX21 = 0; - Float_t y21 = 0.; - Float_t x11 = 0.; + trackHitPattern.ExecuteValidation(trackStore, triggerTrackStore, triggerStore); +} - for (Int_t i=0; iUncheckedAt(i); + //__________________________________________________________________________ +void AliMUONVTrackReconstructor::EventReconstructTrigger(const AliMUONTriggerCircuit& circuit, + const AliMUONVTriggerStore& triggerStore, + AliMUONVTriggerTrackStore& triggerTrackStore) +{ + /// To make the trigger tracks from Local Trigger + AliDebug(1, ""); + AliCodeTimerAuto(""); + + AliMUONGlobalTrigger* globalTrigger = triggerStore.Global(); + + UChar_t gloTrigPat = 0; - AliDebug(1, "AliMUONTrackReconstructor::MakeTriggerTrack using NEW trigger \n"); - AliMUONTriggerCircuit* circuit = - (AliMUONTriggerCircuit*)fTriggerCircuit->At(locTrg->LoCircuit()-1); // -1 !!! + if (globalTrigger) + { + gloTrigPat = globalTrigger->GetGlobalResponse(); + } + + TIter next(triggerStore.CreateIterator()); + AliMUONLocalTrigger* locTrg(0x0); - y11 = circuit->GetY11Pos(locTrg->LoStripX()); - stripX21 = locTrg->LoStripX()+locTrg->LoDev()+1; - y21 = circuit->GetY21Pos(stripX21); - x11 = circuit->GetX11Pos(locTrg->LoStripY()); + Float_t z11 = AliMUONConstants::DefaultChamberZ(10); + Float_t z21 = AliMUONConstants::DefaultChamberZ(12); + + AliMUONTriggerTrack triggerTrack; + + while ( ( locTrg = static_cast(next()) ) ) + { + Bool_t xTrig=locTrg->IsTrigX(); + Bool_t yTrig=locTrg->IsTrigY(); + + Int_t localBoardId = locTrg->LoCircuit(); + + if (xTrig && yTrig) + { // make Trigger Track if trigger in X and Y + + Float_t y11 = circuit.GetY11Pos(localBoardId, locTrg->LoStripX()); + // need first to convert deviation to [0-30] + // (see AliMUONLocalTriggerBoard::LocalTrigger) + Int_t deviation = locTrg->GetDeviation(); + 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 %d %f %f %f \n",i,locTrg->LoCircuit(), - locTrg->LoStripX(),locTrg->LoStripX()+locTrg->LoDev()+1,locTrg->LoStripY(),y11, y21, x11)); + 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) ); - fTriggerTrack->SetX11(x11); - fTriggerTrack->SetY11(y11); - fTriggerTrack->SetThetax(thetax); - fTriggerTrack->SetThetay(thetay); - fTriggerTrack->SetGTPattern(gloTrigPat); - - fMUONData->AddRecTriggerTrack(*fTriggerTrack); - } // 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()); - } + 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 }