X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=MUON%2FAliMUONVTrackReconstructor.cxx;h=114b73d973455b988c66e310fdc2d95da5f50e8b;hb=88ab1a2c7f0cd98feda2e9819974ee793a60dc8f;hp=f80e0162668ff1c92670cdafaa9414f9139317cd;hpb=89c8d66d74ab291ce455cd4747abf6fa84796e11;p=u%2Fmrichter%2FAliRoot.git diff --git a/MUON/AliMUONVTrackReconstructor.cxx b/MUON/AliMUONVTrackReconstructor.cxx index f80e0162668..114b73d9734 100644 --- a/MUON/AliMUONVTrackReconstructor.cxx +++ b/MUON/AliMUONVTrackReconstructor.cxx @@ -74,10 +74,19 @@ #include "AliMUONVTriggerStore.h" #include "AliMUONVTriggerTrackStore.h" #include "AliMUONRecoParam.h" +#include "AliMUONGeometryTransformer.h" +#include "AliMUONVDigit.h" #include "AliMpDEManager.h" #include "AliMpArea.h" +#include "AliMpDDLStore.h" +#include "AliMpVSegmentation.h" +#include "AliMpSegmentation.h" +#include "AliMpPad.h" +#include "AliMpDetElement.h" +#include "AliMpCathodType.h" + #include "AliLog.h" #include "AliCodeTimer.h" #include "AliTracker.h" @@ -89,18 +98,23 @@ #include +using std::cout; +using std::endl; /// \cond CLASSIMP ClassImp(AliMUONVTrackReconstructor) // Class implementation in ROOT context /// \endcond //__________________________________________________________________________ AliMUONVTrackReconstructor::AliMUONVTrackReconstructor(const AliMUONRecoParam* recoParam, - AliMUONVClusterServer* clusterServer) + AliMUONVClusterServer* clusterServer, + const AliMUONGeometryTransformer* transformer) : TObject(), fRecTracksPtr(0x0), fNRecTracks(0), fClusterServer(clusterServer), -fkRecoParam(recoParam) +fkRecoParam(recoParam), +fkTransformer(transformer), +fMaxMCSAngle2(0x0) { /// Constructor for class AliMUONVTrackReconstructor /// WARNING: if clusterServer=0x0, no clusterization will be possible at this level @@ -110,6 +124,15 @@ fkRecoParam(recoParam) // set the magnetic field for track extrapolations AliMUONTrackExtrap::SetField(); + + // set the maximum MCS angle in chamber from the minimum acceptable momentum + AliMUONTrackParam param; + Double_t inverseBendingP = (GetRecoParam()->GetMinBendingMomentum() > 0.) ? 1./GetRecoParam()->GetMinBendingMomentum() : 1.; + param.SetInverseBendingMomentum(inverseBendingP); + fMaxMCSAngle2 = new Double_t [AliMUONConstants::NTrackingCh()]; + for (Int_t iCh=0; iChImproveTracks()) ImproveTracks(); // Remove connected tracks - if (GetRecoParam()->RemoveConnectedTracksInSt12()) RemoveConnectedTracks(kFALSE); - else RemoveConnectedTracks(kTRUE); + RemoveConnectedTracks(3, 4, kFALSE); + RemoveConnectedTracks(2, 2, kFALSE); + if (GetRecoParam()->RemoveConnectedTracksInSt12()) RemoveConnectedTracks(0, 1, kFALSE); // Fill AliMUONTrack data members Finalize(); + if (!GetRecoParam()->RemoveConnectedTracksInSt12()) TagConnectedTracks(0, 1, kTRUE); + + // Make sure there is no bad track left + RemoveBadTracks(); + + // Refit the reconstructed tracks with a different resolution for mono-cathod clusters + if (GetRecoParam()->DiscardMonoCathodClusters()) DiscardMonoCathodClusters(); // Add tracks to MUON data container - for (Int_t i=0; iAt(i); track->SetUniqueID(i+1); trackStore.Add(*track); } + } - //__________________________________________________________________________ +//__________________________________________________________________________ +Bool_t AliMUONVTrackReconstructor::IsAcceptable(AliMUONTrackParam &trackParam) +{ + /// Return kTRUE if the track is within given limits on momentum/angle/origin + + const TMatrixD& kParamCov = trackParam.GetCovariances(); + Int_t chamber = trackParam.GetClusterPtr()->GetChamberId(); + Double_t z = trackParam.GetZ(); + Double_t sigmaCut = GetRecoParam()->GetSigmaCutForTracking(); + + // MCS dipersion + Double_t angleMCS2 = 0.; + Double_t impactMCS2 = 0.; + if (AliMUONTrackExtrap::IsFieldON() && chamber < 6) { + + // track momentum is known + for (Int_t iCh=0; iCh<=chamber; iCh++) { + Double_t localMCS2 = AliMUONTrackExtrap::GetMCSAngle2(trackParam, AliMUONConstants::ChamberThicknessInX0(iCh), 1.); + angleMCS2 += localMCS2; + impactMCS2 += AliMUONConstants::DefaultChamberZ(chamber) * AliMUONConstants::DefaultChamberZ(chamber) * localMCS2; + } + + } else { + + // track momentum is unknown + for (Int_t iCh=0; iCh<=chamber; iCh++) { + angleMCS2 += fMaxMCSAngle2[iCh]; + impactMCS2 += AliMUONConstants::DefaultChamberZ(chamber) * AliMUONConstants::DefaultChamberZ(chamber) * fMaxMCSAngle2[iCh]; + } + + } + + // ------ track selection in non bending direction ------ + if (GetRecoParam()->SelectOnTrackSlope()) { + + // check if non bending slope is within tolerances + Double_t nonBendingSlopeErr = TMath::Sqrt(kParamCov(1,1) + angleMCS2); + if ((TMath::Abs(trackParam.GetNonBendingSlope()) - sigmaCut * nonBendingSlopeErr) > GetRecoParam()->GetMaxNonBendingSlope()) return kFALSE; + + } else { + + // or check if non bending impact parameter is within tolerances + Double_t nonBendingImpactParam = TMath::Abs(trackParam.GetNonBendingCoor() - z * trackParam.GetNonBendingSlope()); + Double_t nonBendingImpactParamErr = TMath::Sqrt(kParamCov(0,0) + z * z * kParamCov(1,1) - 2. * z * kParamCov(0,1) + impactMCS2); + if ((nonBendingImpactParam - sigmaCut * nonBendingImpactParamErr) > (3. * GetRecoParam()->GetNonBendingVertexDispersion())) return kFALSE; + + } + + // ------ track selection in bending direction ------ + if (AliMUONTrackExtrap::IsFieldON()) { // depending whether the field is ON or OFF + + // check if bending momentum is within tolerances + Double_t bendingMomentum = TMath::Abs(1. / trackParam.GetInverseBendingMomentum()); + Double_t bendingMomentumErr = TMath::Sqrt(kParamCov(4,4)) * bendingMomentum * bendingMomentum; + if (chamber < 6 && (bendingMomentum + sigmaCut * bendingMomentumErr) < GetRecoParam()->GetMinBendingMomentum()) return kFALSE; + else if ((bendingMomentum + 3. * bendingMomentumErr) < GetRecoParam()->GetMinBendingMomentum()) return kFALSE; + + } else { + + if (GetRecoParam()->SelectOnTrackSlope()) { + + // check if bending slope is within tolerances + Double_t bendingSlopeErr = TMath::Sqrt(kParamCov(3,3) + angleMCS2); + if ((TMath::Abs(trackParam.GetBendingSlope()) - sigmaCut * bendingSlopeErr) > GetRecoParam()->GetMaxBendingSlope()) return kFALSE; + + } else { + + // or check if bending impact parameter is within tolerances + Double_t bendingImpactParam = TMath::Abs(trackParam.GetBendingCoor() - z * trackParam.GetBendingSlope()); + Double_t bendingImpactParamErr = TMath::Sqrt(kParamCov(2,2) + z * z * kParamCov(3,3) - 2. * z * kParamCov(2,3) + impactMCS2); + if ((bendingImpactParam - sigmaCut * bendingImpactParamErr) > (3. * GetRecoParam()->GetBendingVertexDispersion())) return kFALSE; + + } + + } + + return kTRUE; + +} + +//__________________________________________________________________________ TClonesArray* AliMUONVTrackReconstructor::MakeSegmentsBetweenChambers(const AliMUONVClusterStore& clusterStore, Int_t ch1, Int_t ch2) { /// To make the list of segments from the list of clusters in the 2 given chambers. - /// Return a new TClonesArray of segments. - /// It is the responsibility of the user to delete it afterward. + /// Return a TClonesArray of new segments (segments made in a previous call of this function are removed). AliDebug(1,Form("Enter MakeSegmentsBetweenChambers (1..) %d-%d", ch1+1, ch2+1)); + AliCodeTimerAuto("",0); AliMUONVCluster *cluster1, *cluster2; AliMUONObjectPair *segment; - Double_t nonBendingSlope = 0, bendingSlope = 0, impactParam = 0., bendingMomentum = 0.; // to avoid compilation warning + Double_t z1 = 0., z2 = 0., dZ = 0.; + Double_t nonBendingSlope = 0., nonBendingSlopeErr = 0., nonBendingImpactParam = 0., nonBendingImpactParamErr = 0.; + Double_t bendingSlope = 0., bendingSlopeErr = 0., bendingImpactParam = 0., bendingImpactParamErr = 0., bendingImpactParamErr2 = 0.; + Double_t bendingMomentum = 0., bendingMomentumErr = 0.; + Double_t bendingVertexDispersion2 = GetRecoParam()->GetBendingVertexDispersion() * GetRecoParam()->GetBendingVertexDispersion(); + Double_t angleMCS2 = 0.; // maximum angular dispersion**2 due to MCS in chamber + Double_t impactMCS2 = 0.; // maximum impact parameter dispersion**2 due to MCS in chamber + for (Int_t iCh=0; iCh<=ch1; iCh++) { + angleMCS2 += fMaxMCSAngle2[iCh]; + impactMCS2 += AliMUONConstants::DefaultChamberZ(iCh) * AliMUONConstants::DefaultChamberZ(iCh) * fMaxMCSAngle2[iCh]; + } + Double_t sigmaCut = GetRecoParam()->GetSigmaCutForTracking(); // 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); + static TClonesArray *segments = new TClonesArray("AliMUONObjectPair", 100); + segments->Clear("C"); // Loop over clusters in the first chamber of the station while ( ( cluster1 = static_cast(nextInCh1()) ) ) { + z1 = cluster1->GetZ(); // reset cluster iterator of chamber 2 nextInCh2.Reset(); // Loop over clusters in the second chamber of the station while ( ( cluster2 = static_cast(nextInCh2()) ) ) { + z2 = cluster2->GetZ(); + dZ = z1 - z2; - // non bending slope - nonBendingSlope = (cluster1->GetX() - cluster2->GetX()) / (cluster1->GetZ() - cluster2->GetZ()); - - // check if non bending slope is within tolerances - if (TMath::Abs(nonBendingSlope) > GetRecoParam()->GetMaxNonBendingSlope()) continue; - - // bending slope - bendingSlope = (cluster1->GetY() - cluster2->GetY()) / (cluster1->GetZ() - cluster2->GetZ()); - - // check the bending momentum of the bending slope depending if the field is ON or OFF - if (AliMUONTrackExtrap::IsFieldON()) { + // ------ track selection in non bending direction ------ + nonBendingSlope = (cluster1->GetX() - cluster2->GetX()) / dZ; + if (GetRecoParam()->SelectOnTrackSlope()) { + + // check if non bending slope is within tolerances + nonBendingSlopeErr = TMath::Sqrt((cluster1->GetErrX2() + cluster2->GetErrX2()) / dZ / dZ + angleMCS2); + if ((TMath::Abs(nonBendingSlope) - sigmaCut * nonBendingSlopeErr) > GetRecoParam()->GetMaxNonBendingSlope()) continue; + + } else { - // impact parameter - impactParam = cluster1->GetY() - cluster1->GetZ() * bendingSlope; + // or check if non bending impact parameter is within tolerances + nonBendingImpactParam = TMath::Abs(cluster1->GetX() - cluster1->GetZ() * nonBendingSlope); + nonBendingImpactParamErr = TMath::Sqrt((z1 * z1 * cluster2->GetErrX2() + z2 * z2 * cluster1->GetErrX2()) / dZ / dZ + impactMCS2); + if ((nonBendingImpactParam - sigmaCut * nonBendingImpactParamErr) > (3. * GetRecoParam()->GetNonBendingVertexDispersion())) continue; - // absolute value of bending momentum - bendingMomentum = TMath::Abs(AliMUONTrackExtrap::GetBendingMomentumFromImpactParam(impactParam)); + } + + // ------ track selection in bending direction ------ + bendingSlope = (cluster1->GetY() - cluster2->GetY()) / dZ; + if (AliMUONTrackExtrap::IsFieldON()) { // depending whether the field is ON or OFF // check if bending momentum is within tolerances - if (bendingMomentum < GetRecoParam()->GetMinBendingMomentum() || - bendingMomentum > GetRecoParam()->GetMaxBendingMomentum()) continue; + bendingImpactParam = cluster1->GetY() - cluster1->GetZ() * bendingSlope; + bendingImpactParamErr2 = (z1 * z1 * cluster2->GetErrY2() + z2 * z2 * cluster1->GetErrY2()) / dZ / dZ + impactMCS2; + bendingMomentum = TMath::Abs(AliMUONTrackExtrap::GetBendingMomentumFromImpactParam(bendingImpactParam)); + bendingMomentumErr = TMath::Sqrt((bendingVertexDispersion2 + bendingImpactParamErr2) / + bendingImpactParam / bendingImpactParam + 0.01) * bendingMomentum; + if ((bendingMomentum + 3. * bendingMomentumErr) < GetRecoParam()->GetMinBendingMomentum()) continue; } else { - // check if non bending slope is within tolerances - if (TMath::Abs(bendingSlope) > GetRecoParam()->GetMaxBendingSlope()) continue; - + if (GetRecoParam()->SelectOnTrackSlope()) { + + // check if bending slope is within tolerances + bendingSlopeErr = TMath::Sqrt((cluster1->GetErrY2() + cluster2->GetErrY2()) / dZ / dZ + angleMCS2); + if ((TMath::Abs(bendingSlope) - sigmaCut * bendingSlopeErr) > GetRecoParam()->GetMaxBendingSlope()) continue; + + } else { + + // or check if bending impact parameter is within tolerances + bendingImpactParam = TMath::Abs(cluster1->GetY() - cluster1->GetZ() * bendingSlope); + bendingImpactParamErr = TMath::Sqrt((z1 * z1 * cluster2->GetErrY2() + z2 * z2 * cluster1->GetErrY2()) / dZ / dZ + impactMCS2); + if ((bendingImpactParam - sigmaCut * bendingImpactParamErr) > (3. * GetRecoParam()->GetBendingVertexDispersion())) continue; + + } + } + // make new segment segment = new ((*segments)[segments->GetLast()+1]) AliMUONObjectPair(cluster1, cluster2, kFALSE, kFALSE); @@ -402,54 +550,109 @@ void AliMUONVTrackReconstructor::RemoveDoubleTracks() } //__________________________________________________________________________ -void AliMUONVTrackReconstructor::RemoveConnectedTracks(Bool_t inSt345) +void AliMUONVTrackReconstructor::RemoveBadTracks() { - /// 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; + /// Remove tracks for which a problem occured somewhere during the tracking + + AliMUONTrack *track, *nextTrack; + Bool_t trackRemoved = kFALSE; + + track = (AliMUONTrack*) fRecTracksPtr->First(); + while (track) { + + nextTrack = (AliMUONTrack*) fRecTracksPtr->After(track); + + if (track->GetGlobalChi2() >= AliMUONTrack::MaxChi2()) { + AliWarning("problem occured somewhere during the tracking --> discard track"); + fRecTracksPtr->Remove(track); + fNRecTracks--; + trackRemoved = kTRUE; + } + + track = nextTrack; + + } + + // compress array of tracks if needed + if (trackRemoved) fRecTracksPtr->Compress(); + +} + + //__________________________________________________________________________ +void AliMUONVTrackReconstructor::RemoveConnectedTracks(Int_t stMin, Int_t stMax, Bool_t all) +{ + /// Find and remove tracks sharing 1 cluster or more in station(s) [stMin, stMax]. + /// For each couple of connected tracks, one removes the one with the smallest + /// number of clusters or with the highest chi2 value in case of equality. + /// If all=kTRUE: both tracks are removed. + + // tag the tracks to be removed + TagConnectedTracks(stMin, stMax, all); + + // remove them + Int_t nTracks = fRecTracksPtr->GetEntriesFast(); + for (Int_t i = 0; i < nTracks; i++) { + if (((AliMUONTrack*) fRecTracksPtr->UncheckedAt(i))->IsConnected()) { + fRecTracksPtr->RemoveAt(i); + fNRecTracks--; + } + } + + // remove holes in the array if any + fRecTracksPtr->Compress(); +} + + //__________________________________________________________________________ +void AliMUONVTrackReconstructor::TagConnectedTracks(Int_t stMin, Int_t stMax, Bool_t all) +{ + /// Find and tag tracks sharing 1 cluster or more in station(s) [stMin, stMax]. + /// For each couple of connected tracks, one tags the one with the smallest + /// number of clusters or with the highest chi2 value in case of equality. + /// If all=kTRUE: both tracks are tagged. + + AliMUONTrack *track1, *track2; + Int_t nClusters1, nClusters2; + Int_t nTracks = fRecTracksPtr->GetEntriesFast(); + + // reset the tags + for (Int_t i = 0; i < nTracks; i++) ((AliMUONTrack*) fRecTracksPtr->UncheckedAt(i))->Connected(kFALSE); + // Loop over first track of the pair - track1 = (AliMUONTrack*) fRecTracksPtr->First(); - while (track1) { - removedTrack1 = kFALSE; - nClusters1 = track1->GetNClusters(); + for (Int_t iTrack1 = 0; iTrack1 < nTracks; iTrack1++) { + track1 = (AliMUONTrack*) fRecTracksPtr->UncheckedAt(iTrack1); + // 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 - if (inSt345) clustersInCommon = track1->ClustersInCommonInSt345(track2); - else clustersInCommon = track1->ClustersInCommon(track2); - // 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; + for (Int_t iTrack2 = iTrack1+1; iTrack2 < nTracks; iTrack2++) { + track2 = (AliMUONTrack*) fRecTracksPtr->UncheckedAt(iTrack2); + + // check for connected tracks and tag them + if (track1->ClustersInCommon(track2, stMin, stMax) > 0) { + + if (all) { + + // tag both tracks + track1->Connected(); + track2->Connected(); + + } else { + + // tag only the worst track + nClusters1 = track1->GetNClusters(); + nClusters2 = track2->GetNClusters(); + if ((nClusters1 > nClusters2) || ((nClusters1 == nClusters2) && (track1->GetGlobalChi2() <= track2->GetGlobalChi2()))) { + track2->Connected(); + } else { + track1->Connected(); + } + + } + + } + + } + + } + } //__________________________________________________________________________ @@ -467,18 +670,20 @@ void AliMUONVTrackReconstructor::AskForNewClustersInChamber(const AliMUONTrackPa // extrapolate track parameters to the chamber AliMUONTrackParam extrapTrackParam(trackParam); - AliMUONTrackExtrap::ExtrapToZCov(&extrapTrackParam, AliMUONConstants::DefaultChamberZ(chamber)); + if (!AliMUONTrackExtrap::ExtrapToZCov(&extrapTrackParam, AliMUONConstants::DefaultChamberZ(chamber))) return; - // build the searching area using the track resolution and the maximum-distance-to-track value + // build the searching area using the track and chamber resolutions and the maximum-distance-to-track value const TMatrixD& kParamCov = extrapTrackParam.GetCovariances(); - Double_t errX2 = kParamCov(0,0) + kMaxDZ * kMaxDZ * kParamCov(1,1) + 2. * kMaxDZ * TMath::Abs(kParamCov(0,1)); - Double_t errY2 = kParamCov(2,2) + kMaxDZ * kMaxDZ * kParamCov(3,3) + 2. * kMaxDZ * TMath::Abs(kParamCov(2,3)); + Double_t errX2 = kParamCov(0,0) + kMaxDZ * kMaxDZ * kParamCov(1,1) + 2. * kMaxDZ * TMath::Abs(kParamCov(0,1)) + + GetRecoParam()->GetDefaultNonBendingReso(chamber) * GetRecoParam()->GetDefaultNonBendingReso(chamber); + Double_t errY2 = kParamCov(2,2) + kMaxDZ * kMaxDZ * kParamCov(3,3) + 2. * kMaxDZ * TMath::Abs(kParamCov(2,3)) + + GetRecoParam()->GetDefaultBendingReso(chamber) * GetRecoParam()->GetDefaultBendingReso(chamber); Double_t dX = TMath::Abs(trackParam.GetNonBendingSlope()) * kMaxDZ + GetRecoParam()->GetMaxNonBendingDistanceToTrack() + - GetRecoParam()->GetSigmaCutForTracking() * TMath::Sqrt(errX2); + GetRecoParam()->GetSigmaCutForTracking() * TMath::Sqrt(2. * errX2); Double_t dY = TMath::Abs(trackParam.GetBendingSlope()) * kMaxDZ + GetRecoParam()->GetMaxBendingDistanceToTrack() + - GetRecoParam()->GetSigmaCutForTracking() * TMath::Sqrt(errY2); + GetRecoParam()->GetSigmaCutForTracking() * TMath::Sqrt(2. * errY2); AliMpArea area(extrapTrackParam.GetNonBendingCoor(), extrapTrackParam.GetBendingCoor(), dX, dY); @@ -507,11 +712,11 @@ Double_t AliMUONVTrackReconstructor::TryOneCluster(const AliMUONTrackParam &trac /// return trackParamAtCluster // extrapolate track parameters and covariances at the z position of the tested cluster + // and set pointer to cluster into trackParamAtCluster trackParamAtCluster = trackParam; - AliMUONTrackExtrap::ExtrapToZCov(&trackParamAtCluster, cluster->GetZ(), updatePropagator); - - // set pointer to cluster into trackParamAtCluster trackParamAtCluster.SetClusterPtr(cluster); + if (!AliMUONTrackExtrap::ExtrapToZCov(&trackParamAtCluster, cluster->GetZ(), updatePropagator)) + return 2.*AliMUONTrack::MaxChi2(); // Set differences between trackParam and cluster in the bending and non bending directions Double_t dX = cluster->GetX() - trackParamAtCluster.GetNonBendingCoor(); @@ -521,9 +726,12 @@ Double_t AliMUONVTrackReconstructor::TryOneCluster(const AliMUONTrackParam &trac const TMatrixD& kParamCov = trackParamAtCluster.GetCovariances(); Double_t sigmaX2 = kParamCov(0,0) + cluster->GetErrX2(); Double_t sigmaY2 = kParamCov(2,2) + cluster->GetErrY2(); + Double_t covXY = kParamCov(0,2); + Double_t det = sigmaX2 * sigmaY2 - covXY * covXY; // Compute chi2 - return dX * dX / sigmaX2 + dY * dY / sigmaY2; + if (det == 0.) return 2.*AliMUONTrack::MaxChi2(); + return (dX * dX * sigmaY2 + dY * dY * sigmaX2 - 2. * dX * dY * covXY) / det; } @@ -531,7 +739,7 @@ Double_t AliMUONVTrackReconstructor::TryOneCluster(const AliMUONTrackParam &trac Bool_t AliMUONVTrackReconstructor::TryOneClusterFast(const AliMUONTrackParam &trackParam, const AliMUONVCluster* cluster) { /// Test the compatibility between the track and the cluster -/// given the track resolution + the maximum-distance-to-track value +/// given the track and cluster resolutions + the maximum-distance-to-track value /// and assuming linear propagation of the track: /// return kTRUE if they are compatibles @@ -539,12 +747,12 @@ Bool_t AliMUONVTrackReconstructor::TryOneClusterFast(const AliMUONTrackParam &tr Double_t dX = cluster->GetX() - (trackParam.GetNonBendingCoor() + trackParam.GetNonBendingSlope() * dZ); Double_t dY = cluster->GetY() - (trackParam.GetBendingCoor() + trackParam.GetBendingSlope() * dZ); const TMatrixD& kParamCov = trackParam.GetCovariances(); - Double_t errX2 = kParamCov(0,0) + dZ * dZ * kParamCov(1,1) + 2. * dZ * kParamCov(0,1); - Double_t errY2 = kParamCov(2,2) + dZ * dZ * kParamCov(3,3) + 2. * dZ * kParamCov(2,3); + Double_t errX2 = kParamCov(0,0) + dZ * dZ * kParamCov(1,1) + 2. * dZ * kParamCov(0,1) + cluster->GetErrX2(); + Double_t errY2 = kParamCov(2,2) + dZ * dZ * kParamCov(3,3) + 2. * dZ * kParamCov(2,3) + cluster->GetErrY2(); - Double_t dXmax = GetRecoParam()->GetSigmaCutForTracking() * TMath::Sqrt(errX2) + + Double_t dXmax = GetRecoParam()->GetSigmaCutForTracking() * TMath::Sqrt(2. * errX2) + GetRecoParam()->GetMaxNonBendingDistanceToTrack(); - Double_t dYmax = GetRecoParam()->GetSigmaCutForTracking() * TMath::Sqrt(errY2) + + Double_t dYmax = GetRecoParam()->GetSigmaCutForTracking() * TMath::Sqrt(2. * errY2) + GetRecoParam()->GetMaxBendingDistanceToTrack(); if (TMath::Abs(dX) > dXmax || TMath::Abs(dY) > dYmax) return kFALSE; @@ -564,7 +772,7 @@ Double_t AliMUONVTrackReconstructor::TryTwoClustersFast(const AliMUONTrackParam // extrapolate linearly track parameters and covariances at the z position of the second cluster trackParamAtCluster2 = trackParamAtCluster1; - AliMUONTrackExtrap::LinearExtrapToZ(&trackParamAtCluster2, cluster2->GetZ()); + AliMUONTrackExtrap::LinearExtrapToZCov(&trackParamAtCluster2, cluster2->GetZ()); // set pointer to cluster2 into trackParamAtCluster2 trackParamAtCluster2.SetClusterPtr(cluster2); @@ -590,7 +798,7 @@ Double_t AliMUONVTrackReconstructor::TryTwoClustersFast(const AliMUONTrackParam // Compute chi2 Double_t detX = sigma2X1 * sigma2X2 - covX1X2 * covX1X2; Double_t detY = sigma2Y1 * sigma2Y2 - covY1Y2 * covY1Y2; - if (detX == 0. || detY == 0.) return 1.e10; + if (detX == 0. || detY == 0.) return 2.*AliMUONTrack::MaxChi2(); return (dX1 * dX1 * sigma2X2 + dX2 * dX2 * sigma2X1 - 2. * dX1 * dX2 * covX1X2) / detX + (dY1 * dY1 * sigma2Y2 + dY2 * dY2 * sigma2Y1 - 2. * dY1 * dY2 * covY1Y2) / detY; @@ -609,7 +817,7 @@ Bool_t AliMUONVTrackReconstructor::FollowLinearTrackInChamber(AliMUONTrack &trac /// 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 chi2WithOneCluster = AliMUONTrack::MaxChi2(); Double_t maxChi2WithOneCluster = 2. * GetRecoParam()->GetSigmaCutForTracking() * GetRecoParam()->GetSigmaCutForTracking(); // 2 because 2 quantities in chi2 Double_t bestChi2WithOneCluster = maxChi2WithOneCluster; @@ -632,7 +840,7 @@ Bool_t AliMUONVTrackReconstructor::FollowLinearTrackInChamber(AliMUONTrack &trac } // Add MCS effect - AliMUONTrackExtrap::AddMCSEffect(&trackParam,AliMUONConstants::ChamberThicknessInX0(),1.); + AliMUONTrackExtrap::AddMCSEffect(&trackParam,AliMUONConstants::ChamberThicknessInX0(trackParam.GetClusterPtr()->GetChamberId()),-1.); // Printout for debuging if ((AliLog::GetDebugLevel("MUON","AliMUONVTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) { @@ -650,7 +858,7 @@ Bool_t AliMUONVTrackReconstructor::FollowLinearTrackInChamber(AliMUONTrack &trac // try to add the current cluster accuratly extrapTrackParamAtCluster = trackParam; - AliMUONTrackExtrap::LinearExtrapToZ(&extrapTrackParamAtCluster, cluster->GetZ()); + AliMUONTrackExtrap::LinearExtrapToZCov(&extrapTrackParamAtCluster, cluster->GetZ()); chi2WithOneCluster = TryOneCluster(extrapTrackParamAtCluster, cluster, extrapTrackParamAtCluster); // if good chi2 then consider to add cluster @@ -744,8 +952,8 @@ Bool_t AliMUONVTrackReconstructor::FollowLinearTrackInStation(AliMUONTrack &trac ch2 = 2*nextStation+1; } - Double_t chi2WithOneCluster = 1.e10; - Double_t chi2WithTwoClusters = 1.e10; + Double_t chi2WithOneCluster = AliMUONTrack::MaxChi2(); + Double_t chi2WithTwoClusters = AliMUONTrack::MaxChi2(); Double_t maxChi2WithOneCluster = 2. * GetRecoParam()->GetSigmaCutForTracking() * GetRecoParam()->GetSigmaCutForTracking(); // 2 because 2 quantities in chi2 Double_t maxChi2WithTwoClusters = 4. * GetRecoParam()->GetSigmaCutForTracking() * @@ -779,7 +987,7 @@ Bool_t AliMUONVTrackReconstructor::FollowLinearTrackInStation(AliMUONTrack &trac } // Add MCS effect - AliMUONTrackExtrap::AddMCSEffect(&trackParam,AliMUONConstants::ChamberThicknessInX0(),1.); + AliMUONTrackExtrap::AddMCSEffect(&trackParam,AliMUONConstants::ChamberThicknessInX0(trackParam.GetClusterPtr()->GetChamberId()),-1.); // Printout for debuging if ((AliLog::GetDebugLevel("MUON","AliMUONVTrackReconstructor") >= 1) || (AliLog::GetGlobalDebugLevel() >= 1)) { @@ -798,7 +1006,7 @@ Bool_t AliMUONVTrackReconstructor::FollowLinearTrackInStation(AliMUONTrack &trac // try to add the current cluster accuratly extrapTrackParamAtCluster2 = trackParam; - AliMUONTrackExtrap::LinearExtrapToZ(&extrapTrackParamAtCluster2, clusterCh2->GetZ()); + AliMUONTrackExtrap::LinearExtrapToZCov(&extrapTrackParamAtCluster2, clusterCh2->GetZ()); chi2WithOneCluster = TryOneCluster(extrapTrackParamAtCluster2, clusterCh2, extrapTrackParamAtCluster2); // if good chi2 then try to attach a cluster in the other chamber too @@ -814,7 +1022,7 @@ Bool_t AliMUONVTrackReconstructor::FollowLinearTrackInStation(AliMUONTrack &trac } // add MCS effect - AliMUONTrackExtrap::AddMCSEffect(&extrapTrackParamAtCluster2,AliMUONConstants::ChamberThicknessInX0(),1.); + AliMUONTrackExtrap::AddMCSEffect(&extrapTrackParamAtCluster2,AliMUONConstants::ChamberThicknessInX0(ch2),-1.); // reset cluster iterator of chamber 1 nextInCh1.Reset(); @@ -914,10 +1122,10 @@ Bool_t AliMUONVTrackReconstructor::FollowLinearTrackInStation(AliMUONTrack &trac } //Extrapolate trackCandidate to chamber "ch2" - AliMUONTrackExtrap::LinearExtrapToZ(&trackParam, AliMUONConstants::DefaultChamberZ(ch2)); + AliMUONTrackExtrap::LinearExtrapToZCov(&trackParam, AliMUONConstants::DefaultChamberZ(ch2)); // add MCS effect for next step - AliMUONTrackExtrap::AddMCSEffect(&trackParam,AliMUONConstants::ChamberThicknessInX0(),1.); + AliMUONTrackExtrap::AddMCSEffect(&trackParam,AliMUONConstants::ChamberThicknessInX0(ch2),-1.); // reset cluster iterator of chamber 1 nextInCh1.Reset(); @@ -934,7 +1142,7 @@ Bool_t AliMUONVTrackReconstructor::FollowLinearTrackInStation(AliMUONTrack &trac // try to add the current cluster accuratly extrapTrackParamAtCluster1 = trackParam; - AliMUONTrackExtrap::LinearExtrapToZ(&extrapTrackParamAtCluster1, clusterCh1->GetZ()); + AliMUONTrackExtrap::LinearExtrapToZCov(&extrapTrackParamAtCluster1, clusterCh1->GetZ()); chi2WithOneCluster = TryOneCluster(extrapTrackParamAtCluster1, clusterCh1, extrapTrackParamAtCluster1); // if good chi2 then consider to add clusterCh1 @@ -1060,19 +1268,134 @@ void AliMUONVTrackReconstructor::ImproveTracks() //__________________________________________________________________________ void AliMUONVTrackReconstructor::Finalize() { - /// Recompute track parameters and covariances at each attached cluster from those at the first one + /// Recompute track parameters and covariances at each attached cluster /// Set the label pointing to the corresponding MC track + /// Remove the track if finalization failed - AliMUONTrack *track; + AliMUONTrack *track, *nextTrack; + Bool_t trackRemoved = kFALSE; track = (AliMUONTrack*) fRecTracksPtr->First(); while (track) { - FinalizeTrack(*track); + nextTrack = (AliMUONTrack*) fRecTracksPtr->After(track); - track->FindMCLabel(); + if (FinalizeTrack(*track)) track->FindMCLabel(); + else { + fRecTracksPtr->Remove(track); + fNRecTracks--; + trackRemoved = kTRUE; + } - track = (AliMUONTrack*) fRecTracksPtr->After(track); + track = nextTrack; + + } + + // compress array of tracks if needed + if (trackRemoved) fRecTracksPtr->Compress(); + +} + +//__________________________________________________________________________ +void AliMUONVTrackReconstructor::DiscardMonoCathodClusters() +{ + /// Assign a different resolution to the mono-cathod clusters + /// in the direction of the missing plane and refit the track + /// Remove the track in case of failure + + if (!fkTransformer) AliFatal("missing geometry transformer"); + + AliMUONTrack *track, *nextTrack; + Bool_t trackRemoved = kFALSE; + + track = (AliMUONTrack*) fRecTracksPtr->First(); + while (track) { + + nextTrack = (AliMUONTrack*) fRecTracksPtr->After(track); + + ChangeMonoCathodClusterRes(*track); + + if (!RefitTrack(*track) || (GetRecoParam()->ImproveTracks() && !track->IsImproved())) { + fRecTracksPtr->Remove(track); + fNRecTracks--; + trackRemoved = kTRUE; + } + + track = nextTrack; + + } + + // compress array of tracks if needed + if (trackRemoved) fRecTracksPtr->Compress(); + +} + +//__________________________________________________________________________ +void AliMUONVTrackReconstructor::ChangeMonoCathodClusterRes(AliMUONTrack &track) +{ + /// Assign a different resolution to the mono-cathod clusters + /// in the direction of the missing plane and refit the track + + // Loop over clusters + AliMUONVCluster *cluster; + Int_t nClusters = track.GetNClusters(); + for (Int_t iCluster = 0; iCluster < nClusters; iCluster++) { + cluster = ((AliMUONTrackParam*) track.GetTrackParamAtCluster()->UncheckedAt(iCluster))->GetClusterPtr(); + + // do it only for stations 3, 4 & 5 + if (cluster->GetChamberId() < 4) continue; + + // get the cathod corresponding to the bending/non-bending plane + Int_t deId = cluster->GetDetElemId(); + AliMpDetElement* de = AliMpDDLStore::Instance()->GetDetElement(deId, kFALSE); + if (!de) continue; + AliMp::CathodType cath1 = de->GetCathodType(AliMp::kBendingPlane); + AliMp::CathodType cath2 = de->GetCathodType(AliMp::kNonBendingPlane); + + // get the corresponding segmentation + const AliMpVSegmentation* seg1 = AliMpSegmentation::Instance()->GetMpSegmentation(deId, cath1); + const AliMpVSegmentation* seg2 = AliMpSegmentation::Instance()->GetMpSegmentation(deId, cath2); + if (!seg1 || !seg2) continue; + + // get local coordinate of the cluster + Double_t lX,lY,lZ; + Double_t gX = cluster->GetX(); + Double_t gY = cluster->GetY(); + Double_t gZ = cluster->GetZ(); + fkTransformer->Global2Local(deId,gX,gY,gZ,lX,lY,lZ); + + // find pads below the cluster + AliMpPad pad1 = seg1->PadByPosition(lX, lY, kFALSE); + AliMpPad pad2 = seg2->PadByPosition(lX, lY, kFALSE); + + // build their ID if pads are valid + UInt_t padId1 = (pad1.IsValid()) ? AliMUONVDigit::BuildUniqueID(deId, pad1.GetManuId(), pad1.GetManuChannel(), cath1) : 0; + UInt_t padId2 = (pad2.IsValid()) ? AliMUONVDigit::BuildUniqueID(deId, pad2.GetManuId(), pad2.GetManuChannel(), cath2) : 0; + + // check if the cluster contains these pads + Bool_t hasNonBending = kFALSE; + Bool_t hasBending = kFALSE; + for (Int_t iDigit = 0; iDigit < cluster->GetNDigits(); iDigit++) { + + UInt_t digitId = cluster->GetDigitId(iDigit); + + if (digitId == padId1) { + + hasBending = kTRUE; + if (hasNonBending) break; + + } else if (digitId == padId2) { + + hasNonBending = kTRUE; + if (hasBending) break; + + } + + } + + // modify the cluster resolution if needed + if (!hasNonBending) cluster->SetErrXY(GetRecoParam()->GetMonoCathodClNonBendingRes(), cluster->GetErrY()); + if (!hasBending) cluster->SetErrXY(cluster->GetErrX(), GetRecoParam()->GetMonoCathodClBendingRes()); } @@ -1085,20 +1408,21 @@ void AliMUONVTrackReconstructor::ValidateTracksWithTrigger(AliMUONVTrackStore& t const AliMUONTrackHitPattern& trackHitPattern) { /// Try to match track from tracking system with trigger track - AliCodeTimerAuto(""); + AliCodeTimerAuto("",0); trackHitPattern.ExecuteValidation(trackStore, triggerTrackStore, triggerStore); } - //__________________________________________________________________________ + +//__________________________________________________________________________ void AliMUONVTrackReconstructor::EventReconstructTrigger(const AliMUONTriggerCircuit& circuit, const AliMUONVTriggerStore& triggerStore, AliMUONVTriggerTrackStore& triggerTrackStore) { - /// To make the trigger tracks from Local Trigger + /// Fill trigger track store from local trigger AliDebug(1, ""); - AliCodeTimerAuto(""); - + AliCodeTimerAuto("",0); + AliMUONGlobalTrigger* globalTrigger = triggerStore.Global(); UChar_t gloTrigPat = 0; @@ -1108,47 +1432,80 @@ void AliMUONVTrackReconstructor::EventReconstructTrigger(const AliMUONTriggerCir gloTrigPat = globalTrigger->GetGlobalResponse(); } + AliMUONTriggerTrack triggerTrack; + TIter next(triggerStore.CreateIterator()); AliMUONLocalTrigger* locTrg(0x0); - - 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) + if ( locTrg->IsTrigX() && locTrg->IsTrigY() ) { // 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()); + if (TriggerToTrack(circuit, *locTrg, triggerTrack, gloTrigPat)) + triggerTrackStore.Add(triggerTrack); + + } // board is fired + } // end of loop on Local Trigger +} + +//__________________________________________________________________________ +Bool_t AliMUONVTrackReconstructor::TriggerToTrack(const AliMUONTriggerCircuit& circuit, + const AliMUONLocalTrigger& locTrg, + AliMUONTriggerTrack& triggerTrack, + UChar_t globalTriggerPattern) +{ + /// To make the trigger tracks from Local Trigger + const Double_t kTrigNonBendReso = AliMUONConstants::TriggerNonBendingReso(); + const Double_t kTrigBendReso = AliMUONConstants::TriggerBendingReso(); + const Double_t kSqrt12 = TMath::Sqrt(12.); + + TMatrixD trigCov(3,3); + + Int_t localBoardId = locTrg.LoCircuit(); - 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 y11 = circuit.GetY11Pos(localBoardId, locTrg.LoStripX()); + Float_t z11 = circuit.GetZ11Pos(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 z21 = circuit.GetZ21Pos(localBoardId, stripX21); + Float_t x11 = circuit.GetX11Pos(localBoardId, locTrg.LoStripY()); - Float_t thetax = TMath::ATan2( x11 , z11 ); - Float_t thetay = TMath::ATan2( (y21-y11) , (z21-z11) ); + AliDebug(1, Form(" MakeTriggerTrack %3d %2d %2d %2d (%f %f %f) (%f %f)\n",locTrg.LoCircuit(), + locTrg.LoStripX(),locTrg.LoStripX()+deviation+1,locTrg.LoStripY(),x11, y11, z11, y21, z21)); - triggerTrack.SetX11(x11); - triggerTrack.SetY11(y11); - triggerTrack.SetThetax(thetax); - triggerTrack.SetThetay(thetay); - triggerTrack.SetGTPattern(gloTrigPat); - triggerTrack.SetLoTrgNum(localBoardId); + if (TMath::Abs(z11) < 0.00001) return kFALSE; + + Double_t deltaZ = z11 - z21; - triggerTrackStore.Add(triggerTrack); - } // board is fired - } // end of loop on Local Trigger -} + Float_t slopeX = x11/z11; + Float_t slopeY = (y11-y21) / deltaZ; + + Float_t sigmaX = circuit.GetX11Width(localBoardId, locTrg.LoStripY()) / kSqrt12; + Float_t sigmaY = circuit.GetY11Width(localBoardId, locTrg.LoStripX()) / kSqrt12; + Float_t sigmaY21 = circuit.GetY21Width(localBoardId, locTrg.LoStripX()) / kSqrt12; + + trigCov.Zero(); + trigCov(0,0) = kTrigNonBendReso * kTrigNonBendReso + sigmaX * sigmaX; + trigCov(1,1) = kTrigBendReso * kTrigBendReso + sigmaY * sigmaY; + trigCov(2,2) = + (2. * kTrigBendReso * kTrigBendReso + sigmaY * sigmaY + sigmaY21 * sigmaY21 ) / deltaZ / deltaZ; + trigCov(1,2) = trigCov(2,1) = trigCov(1,1) / deltaZ; + + triggerTrack.SetX11(x11); + triggerTrack.SetY11(y11); + triggerTrack.SetZ11(z11); + triggerTrack.SetZ21(z21); + triggerTrack.SetSlopeX(slopeX); + triggerTrack.SetSlopeY(slopeY); + triggerTrack.SetGTPattern(globalTriggerPattern); + triggerTrack.SetLoTrgNum(localBoardId); + triggerTrack.SetCovariances(trigCov); + triggerTrack.SetUniqueID(locTrg.GetUniqueID()); + return kTRUE; + +}