X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=MUON%2FAliMUONTrack.cxx;h=8a38f42b8f07055e5ccf7a3a5e269db454c06f09;hb=a6d06a4ffa5aa4bad37d5f5cf31e1d1def9dd2ac;hp=2f0474cc3ce0b2f47c4e879f394db32f9ddfbe4c;hpb=8429a5e4b460a88d69cd5915f9b555d4393d2980;p=u%2Fmrichter%2FAliRoot.git diff --git a/MUON/AliMUONTrack.cxx b/MUON/AliMUONTrack.cxx index 2f0474cc3ce..8a38f42b8f0 100644 --- a/MUON/AliMUONTrack.cxx +++ b/MUON/AliMUONTrack.cxx @@ -13,641 +13,1057 @@ * provided "as is" without express or implied warranty. * **************************************************************************/ -/* -$Log$ -Revision 1.5 2000/07/18 16:04:06 gosset -AliMUONEventReconstructor package: -* a few minor modifications and more comments -* a few corrections - * right sign for Z of raw clusters - * right loop over chambers inside station - * symmetrized covariance matrix for measurements (TrackChi2MCS) - * right sign of charge in extrapolation (ExtrapToZ) - * right zEndAbsorber for Branson correction below 3 degrees -* use of TVirtualFitter instead of TMinuit for AliMUONTrack::Fit -* no parameter for AliMUONTrack::Fit() but more fit parameters in Track object - -Revision 1.4 2000/06/30 10:15:48 gosset -Changes to EventReconstructor...: -precision fit with multiple Coulomb scattering; -extrapolation to vertex with Branson correction in absorber (JPC) - -Revision 1.3 2000/06/25 13:23:28 hristov -stdlib.h needed for non-Linux compilation - -Revision 1.2 2000/06/15 07:58:48 morsch -Code from MUON-dev joined - -Revision 1.1.2.3 2000/06/12 10:11:34 morsch -Dummy copy constructor and assignment operator added - -Revision 1.1.2.2 2000/06/09 12:58:05 gosset -Removed comment beginnings in Log sections of .cxx files -Suppressed most violations of coding rules - -Revision 1.1.2.1 2000/06/07 14:44:53 gosset -Addition of files for track reconstruction in C++ -*/ +/* $Id$ */ -//__________________________________________________________________________ -// +//----------------------------------------------------------------------------- +// Class AliMUONTrack +//------------------- // Reconstructed track in ALICE dimuon spectrometer -//__________________________________________________________________________ +//----------------------------------------------------------------------------- #include "AliMUONTrack.h" -#include - -#include -#include -#include -#include -#include +#include "AliMUONVCluster.h" +#include "AliMUONVClusterStore.h" +#include "AliMUONObjectPair.h" +#include "AliMUONConstants.h" +#include "AliMUONTrackExtrap.h" -#include "AliMUONEventReconstructor.h" -#include "AliMUONHitForRec.h" -#include "AliMUONSegment.h" -#include "AliMUONTrackHit.h" +#include "AliLog.h" -#include - -// Functions to be minimized with Minuit -void TrackChi2(Int_t &NParam, Double_t *Gradient, Double_t &Chi2, Double_t *Param, Int_t Flag); -void TrackChi2MCS(Int_t &NParam, Double_t *Gradient, Double_t &Chi2, Double_t *Param, Int_t Flag); +#include +#include -Double_t MultipleScatteringAngle2(AliMUONTrackHit *TrackHit); +#include +/// \cond CLASSIMP ClassImp(AliMUONTrack) // Class implementation in ROOT context +/// \endcond -TVirtualFitter* AliMUONTrack::fgFitter = NULL; +//__________________________________________________________________________ +AliMUONTrack::AliMUONTrack() + : TObject(), + fTrackParamAtCluster(0x0), + fFitWithVertex(kFALSE), + fVertexErrXY2(), + fFitWithMCS(kFALSE), + fClusterWeightsNonBending(0x0), + fClusterWeightsBending(0x0), + fGlobalChi2(-1.), + fImproved(kFALSE), + fMatchTrigger(-1), + floTrgNum(-1), + fChi2MatchTrigger(0.), + fTrackID(0), + fTrackParamAtVertex(0x0), + fHitsPatternInTrigCh(0), + fLocalTrigger(0) +{ + /// Default constructor + fVertexErrXY2[0] = 0.; + fVertexErrXY2[1] = 0.; +} //__________________________________________________________________________ -AliMUONTrack::AliMUONTrack(AliMUONSegment* BegSegment, AliMUONSegment* EndSegment, AliMUONEventReconstructor* EventReconstructor) +AliMUONTrack::AliMUONTrack(AliMUONObjectPair *segment) + : TObject(), + fTrackParamAtCluster(new TClonesArray("AliMUONTrackParam",10)), + fFitWithVertex(kFALSE), + fVertexErrXY2(), + fFitWithMCS(kFALSE), + fClusterWeightsNonBending(0x0), + fClusterWeightsBending(0x0), + fGlobalChi2(0.), + fImproved(kFALSE), + fMatchTrigger(-1), + floTrgNum(-1), + fChi2MatchTrigger(0.), + fTrackID(0), + fTrackParamAtVertex(0x0), + fHitsPatternInTrigCh(0), + fLocalTrigger(0) { - // Constructor from two Segment's - fEventReconstructor = EventReconstructor; // link back to EventReconstructor - // memory allocation for the TObjArray of pointers to reconstructed TrackHit's - fTrackHitsPtr = new TObjArray(10); - fNTrackHits = 0; - AddSegment(BegSegment); // add hits from BegSegment - AddSegment(EndSegment); // add hits from EndSegment - fTrackHitsPtr->Sort(); // sort TrackHits according to increasing Z - SetTrackParamAtVertex(); // set track parameters at vertex - // set fit conditions... - fFitMCS = 0; - fFitNParam = 3; - fFitStart = 1; - fFitFMin = -1.0; - return; + /// Constructor from two clusters + + fVertexErrXY2[0] = 0.; + fVertexErrXY2[1] = 0.; + + // Pointers to clusters from the segment + AliMUONVCluster* cluster1 = (AliMUONVCluster*) segment->First(); + AliMUONVCluster* cluster2 = (AliMUONVCluster*) segment->Second(); + + // check sorting in -Z (spectro z<0) + if (cluster1->GetZ() < cluster2->GetZ()) { + cluster1 = cluster2; + cluster2 = (AliMUONVCluster*) segment->First(); + } + + // order the clusters into the track according to the station the segment belong to + // to anticipate the direction of propagation in the first tracking step + // (go backward if the segment is on the last station / go forward otherwise) + AliMUONVCluster *firstCluster, *lastCluster; + if (cluster1->GetChamberId() == 8) { // last station + firstCluster = cluster1; + lastCluster = cluster2; + } else { + firstCluster = cluster2; + lastCluster = cluster1; + } + + // Compute track parameters + Double_t z1 = firstCluster->GetZ(); + Double_t z2 = lastCluster->GetZ(); + Double_t dZ = z1 - z2; + // Non bending plane + Double_t nonBendingCoor1 = firstCluster->GetX(); + Double_t nonBendingCoor2 = lastCluster->GetX(); + Double_t nonBendingSlope = (nonBendingCoor1 - nonBendingCoor2) / dZ; + // Bending plane + Double_t bendingCoor1 = firstCluster->GetY(); + Double_t bendingCoor2 = lastCluster->GetY(); + Double_t bendingSlope = (bendingCoor1 - bendingCoor2) / dZ; + // Inverse bending momentum + Double_t bendingImpact = bendingCoor1 - z1 * bendingSlope; + Double_t inverseBendingMomentum = 1. / AliMUONTrackExtrap::GetBendingMomentumFromImpactParam(bendingImpact); + + + // Set track parameters at first cluster (needed by any tracking algorithm) + AliMUONTrackParam trackParamAtFirstCluster; + trackParamAtFirstCluster.SetZ(z1); + trackParamAtFirstCluster.SetNonBendingCoor(nonBendingCoor1); + trackParamAtFirstCluster.SetNonBendingSlope(nonBendingSlope); + trackParamAtFirstCluster.SetBendingCoor(bendingCoor1); + trackParamAtFirstCluster.SetBendingSlope(bendingSlope); + trackParamAtFirstCluster.SetInverseBendingMomentum(inverseBendingMomentum); + + + // Set track parameters at last cluster (used by Kalman only) + AliMUONTrackParam trackParamAtLastCluster; + trackParamAtLastCluster.SetZ(z2); + trackParamAtLastCluster.SetNonBendingCoor(nonBendingCoor2); + trackParamAtLastCluster.SetNonBendingSlope(nonBendingSlope); + trackParamAtLastCluster.SetBendingCoor(bendingCoor2); + trackParamAtLastCluster.SetBendingSlope(bendingSlope); + trackParamAtLastCluster.SetInverseBendingMomentum(inverseBendingMomentum); + + + // Compute and set track parameters covariances at first cluster (needed by any tracking algorithm) + TMatrixD paramCov1(5,5); + paramCov1.Zero(); + // Non bending plane + paramCov1(0,0) = firstCluster->GetErrX2(); + paramCov1(0,1) = firstCluster->GetErrX2() / dZ; + paramCov1(1,0) = paramCov1(0,1); + paramCov1(1,1) = ( firstCluster->GetErrX2() + lastCluster->GetErrX2() ) / dZ / dZ; + // Bending plane + paramCov1(2,2) = firstCluster->GetErrY2(); + paramCov1(2,3) = firstCluster->GetErrY2() / dZ; + paramCov1(3,2) = paramCov1(2,3); + paramCov1(3,3) = ( firstCluster->GetErrY2() + lastCluster->GetErrY2() ) / dZ / dZ; + // Inverse bending momentum (50% error) + paramCov1(4,4) = 0.5*inverseBendingMomentum * 0.5*inverseBendingMomentum; + // Set covariances + trackParamAtFirstCluster.SetCovariances(paramCov1); + + + // Compute and set track parameters covariances at last cluster as if the first cluster did not exist (used by Kalman only) + TMatrixD paramCov2(5,5); + paramCov2.Zero(); + // Non bending plane + paramCov2(0,0) = paramCov1(0,0); + paramCov2(1,1) = 100.*paramCov1(1,1); + // Bending plane + paramCov2(2,2) = paramCov1(2,2); + paramCov2(3,3) = 100.*paramCov1(3,3); + // Inverse bending momentum + paramCov2(4,4) = paramCov1(4,4); + // Set covariances + trackParamAtLastCluster.SetCovariances(paramCov2); + + // Flag clusters as being removable + trackParamAtFirstCluster.SetRemovable(kTRUE); + trackParamAtLastCluster.SetRemovable(kTRUE); + + // Add track parameters at clusters + AddTrackParamAtCluster(trackParamAtFirstCluster,*firstCluster); + AddTrackParamAtCluster(trackParamAtLastCluster,*lastCluster); + +} + +//__________________________________________________________________________ +AliMUONTrack::AliMUONTrack(const AliMUONTrack& track) + : TObject(track), + fTrackParamAtCluster(0x0), + fFitWithVertex(track.fFitWithVertex), + fVertexErrXY2(), + fFitWithMCS(track.fFitWithMCS), + fClusterWeightsNonBending(0x0), + fClusterWeightsBending(0x0), + fGlobalChi2(track.fGlobalChi2), + fImproved(track.fImproved), + fMatchTrigger(track.fMatchTrigger), + floTrgNum(track.floTrgNum), + fChi2MatchTrigger(track.fChi2MatchTrigger), + fTrackID(track.fTrackID), + fTrackParamAtVertex(0x0), + fHitsPatternInTrigCh(track.fHitsPatternInTrigCh), + fLocalTrigger(track.fLocalTrigger) +{ + ///copy constructor + + // necessary to make a copy of the objects and not only the pointers in TClonesArray. + if (track.fTrackParamAtCluster) { + fTrackParamAtCluster = new TClonesArray("AliMUONTrackParam",10); + AliMUONTrackParam *trackParamAtCluster = (AliMUONTrackParam*) track.fTrackParamAtCluster->First(); + while (trackParamAtCluster) { + new ((*fTrackParamAtCluster)[GetNClusters()]) AliMUONTrackParam(*trackParamAtCluster); + trackParamAtCluster = (AliMUONTrackParam*) track.fTrackParamAtCluster->After(trackParamAtCluster); + } + } + + // copy vertex resolution square used during the tracking procedure + fVertexErrXY2[0] = track.fVertexErrXY2[0]; + fVertexErrXY2[1] = track.fVertexErrXY2[1]; + + // copy cluster weights matrices if any + if (track.fClusterWeightsNonBending) fClusterWeightsNonBending = new TMatrixD(*(track.fClusterWeightsNonBending)); + if (track.fClusterWeightsBending) fClusterWeightsBending = new TMatrixD(*(track.fClusterWeightsBending)); + + // copy track parameters at vertex if any + if (track.fTrackParamAtVertex) fTrackParamAtVertex = new AliMUONTrackParam(*(track.fTrackParamAtVertex)); + } //__________________________________________________________________________ -AliMUONTrack::AliMUONTrack(AliMUONSegment* Segment, AliMUONHitForRec* HitForRec, AliMUONEventReconstructor* EventReconstructor) +AliMUONTrack & AliMUONTrack::operator=(const AliMUONTrack& track) { - // Constructor from one Segment and one HitForRec - fEventReconstructor = EventReconstructor; // link back to EventReconstructor - // memory allocation for the TObjArray of pointers to reconstructed TrackHit's - fTrackHitsPtr = new TObjArray(10); - fNTrackHits = 0; - AddSegment(Segment); // add hits from Segment - AddHitForRec(HitForRec); // add HitForRec - fTrackHitsPtr->Sort(); // sort TrackHits according to increasing Z - SetTrackParamAtVertex(); // set track parameters at vertex - // set fit conditions... - fFitMCS = 0; - fFitNParam = 3; - fFitStart = 1; - fFitFMin = -1.0; - return; + /// Asignment operator + // check assignement to self + if (this == &track) + return *this; + + // base class assignement + TObject::operator=(track); + + // clear memory + Clear(); + + // necessary to make a copy of the objects and not only the pointers in TClonesArray + if (track.fTrackParamAtCluster) { + fTrackParamAtCluster = new TClonesArray("AliMUONTrackParam",10); + AliMUONTrackParam *trackParamAtCluster = (AliMUONTrackParam*) track.fTrackParamAtCluster->First(); + while (trackParamAtCluster) { + new ((*fTrackParamAtCluster)[GetNClusters()]) AliMUONTrackParam(*trackParamAtCluster); + trackParamAtCluster = (AliMUONTrackParam*) track.fTrackParamAtCluster->After(trackParamAtCluster); + } + } + + // copy cluster weights matrix if any + if (track.fClusterWeightsNonBending) { + if (fClusterWeightsNonBending) { + fClusterWeightsNonBending->ResizeTo(*(track.fClusterWeightsNonBending)); + *fClusterWeightsNonBending = *(track.fClusterWeightsNonBending); + } else fClusterWeightsNonBending = new TMatrixD(*(track.fClusterWeightsNonBending)); + } + + // copy cluster weights matrix if any + if (track.fClusterWeightsBending) { + if (fClusterWeightsBending) { + fClusterWeightsBending->ResizeTo(*(track.fClusterWeightsBending)); + *fClusterWeightsBending = *(track.fClusterWeightsBending); + } else fClusterWeightsBending = new TMatrixD(*(track.fClusterWeightsBending)); + } + + // copy track parameters at vertex if any + if (track.fTrackParamAtVertex) { + if (fTrackParamAtVertex) *fTrackParamAtVertex = *(track.fTrackParamAtVertex); + else fTrackParamAtVertex = new AliMUONTrackParam(*(track.fTrackParamAtVertex)); + } + + fFitWithVertex = track.fFitWithVertex; + fVertexErrXY2[0] = track.fVertexErrXY2[0]; + fVertexErrXY2[1] = track.fVertexErrXY2[1]; + fFitWithMCS = track.fFitWithMCS; + fGlobalChi2 = track.fGlobalChi2; + fImproved = track.fImproved; + fMatchTrigger = track.fMatchTrigger; + floTrgNum = track.floTrgNum; + fChi2MatchTrigger = track.fChi2MatchTrigger; + fTrackID = track.fTrackID; + fHitsPatternInTrigCh = track.fHitsPatternInTrigCh; + fLocalTrigger = track.fLocalTrigger; + + return *this; } //__________________________________________________________________________ AliMUONTrack::~AliMUONTrack() { - // Destructor - if (fTrackHitsPtr) { - delete fTrackHitsPtr; // delete the TObjArray of pointers to TrackHit's - fTrackHitsPtr = NULL; - } + /// Destructor + delete fTrackParamAtCluster; + delete fClusterWeightsNonBending; + delete fClusterWeightsBending; + delete fTrackParamAtVertex; } //__________________________________________________________________________ -AliMUONTrack::AliMUONTrack (const AliMUONTrack& MUONTrack) +void AliMUONTrack::Clear(Option_t* opt) { -// Dummy copy constructor + /// Clear arrays + if (opt && opt[0] == 'C' && fTrackParamAtCluster) fTrackParamAtCluster->Clear("C"); + else { + delete fTrackParamAtCluster; + fTrackParamAtCluster = 0x0; + } + delete fClusterWeightsNonBending; fClusterWeightsNonBending = 0x0; + delete fClusterWeightsBending; fClusterWeightsBending = 0x0; + delete fTrackParamAtVertex; fTrackParamAtVertex = 0x0; } //__________________________________________________________________________ -AliMUONTrack & AliMUONTrack::operator=(const AliMUONTrack& MUONTrack) +TClonesArray* AliMUONTrack::GetTrackParamAtCluster() const { -// Dummy assignment operator - return *this; + /// return array of track parameters at cluster (create it if needed) + if (!fTrackParamAtCluster) fTrackParamAtCluster = new TClonesArray("AliMUONTrackParam",10); + return fTrackParamAtCluster; } //__________________________________________________________________________ -void AliMUONTrack::Remove() +void AliMUONTrack::AddTrackParamAtCluster(const AliMUONTrackParam &trackParam, AliMUONVCluster &cluster, Bool_t copy) { - // Remove current track from array of tracks, - // and corresponding track hits from array of track hits. - // Compress the TClonesArray it belongs to. - AliMUONTrackHit *nextTrackHit; - AliMUONEventReconstructor *eventRec = this->fEventReconstructor; - TClonesArray *trackHitsPtr = eventRec->GetRecTrackHitsPtr(); - // Loop over all track hits of track - AliMUONTrackHit *trackHit = (AliMUONTrackHit*) fTrackHitsPtr->First(); - while (trackHit) { - nextTrackHit = (AliMUONTrackHit*) fTrackHitsPtr->After(trackHit); - // Remove TrackHit from event TClonesArray. - // Destructor is called, - // hence links between HitForRec's and TrackHit's are updated - trackHitsPtr->Remove(trackHit); - trackHit = nextTrackHit; + /// Copy given track parameters into a new TrackParamAtCluster + /// Link parameters with the associated cluster + /// If copy=kTRUE: the cluster is copied then passed the trackParam which become its owner + /// otherwise: make sure to do not delete the cluster until it is used by the track + + // check chamber ID of the associated cluster + if (cluster.GetChamberId() < 0 || cluster.GetChamberId() > AliMUONConstants::NTrackingCh()) { + AliError(Form("Chamber ID of the associated cluster is not valid (ChamberId=%d)",cluster.GetChamberId())); + return; + } + + // check whether track parameters are given at the correct cluster z position + if (cluster.GetZ() != trackParam.GetZ()) { + AliError("track parameters are given at a different z position than the one of the associated cluster"); + return; } - // Remove the track from event TClonesArray - // Destructor is called, - // hence space for TObjArray of pointers to TrackHit's is freed - eventRec->GetRecTracksPtr()->Remove(this); - // Number of tracks decreased by 1 - eventRec->SetNRecTracks(eventRec->GetNRecTracks() - 1); - // Compress event TClonesArray of Track's: - // this is essential to retrieve the TClonesArray afterwards - eventRec->GetRecTracksPtr()->Compress(); - // Compress event TClonesArray of TrackHit's: - // this is probably also essential to retrieve the TClonesArray afterwards - trackHitsPtr->Compress(); + + // add parameters to the array of track parameters + if (!fTrackParamAtCluster) fTrackParamAtCluster = new TClonesArray("AliMUONTrackParam",10); + AliMUONTrackParam* trackParamAtCluster = new ((*fTrackParamAtCluster)[GetNClusters()]) AliMUONTrackParam(trackParam); + + // link parameters with the associated cluster or its copy + if (copy) { + AliMUONVCluster *clusterCopy = static_cast(cluster.Clone()); + trackParamAtCluster->SetClusterPtr(clusterCopy, kTRUE); + } else trackParamAtCluster->SetClusterPtr(&cluster); + + // sort the array of track parameters + fTrackParamAtCluster->Sort(); } //__________________________________________________________________________ -void AliMUONTrack::SetFitMCS(Int_t FitMCS) +void AliMUONTrack::RemoveTrackParamAtCluster(AliMUONTrackParam *trackParam) { - // Set multiple Coulomb scattering option for track fit "fFitMCS" - // from "FitMCS" argument: 0 without, 1 with - if ((FitMCS == 0) || (FitMCS == 1)) fFitMCS = FitMCS; - // better implementation with enum(with, without) ???? - else { - cout << "ERROR in AliMUONTrack::SetFitMCS(FitMCS)" << endl; - cout << "FitMCS = " << FitMCS << " is neither 0 nor 1" << endl; - exit(0); + /// Remove trackParam from the array of TrackParamAtCluster + if (!fTrackParamAtCluster || !fTrackParamAtCluster->Remove(trackParam)) { + AliWarning("object to remove does not exist in array fTrackParamAtCluster"); + return; } - return; + + fTrackParamAtCluster->Compress(); } //__________________________________________________________________________ -void AliMUONTrack::SetFitNParam(Int_t FitNParam) +void AliMUONTrack::UpdateTrackParamAtCluster() { - // Set number of parameters for track fit "fFitNParam" from "FitNParam": - // 3 for momentum, 5 for momentum and position - if ((FitNParam == 3) || (FitNParam == 5)) fFitNParam = FitNParam; - else { - cout << "ERROR in AliMUONTrack::SetFitNParam(FitNParam)" << endl; - cout << "FitNParam = " << FitNParam << " is neither 3 nor 5" << endl; - exit(0); + /// Update track parameters at each attached cluster + + if (GetNClusters() == 0) { + AliWarning("no cluster attached to the track"); + return; } - return; + + AliMUONTrackParam* startingTrackParam = (AliMUONTrackParam*) fTrackParamAtCluster->First(); + AliMUONTrackParam* trackParamAtCluster = (AliMUONTrackParam*) fTrackParamAtCluster->After(startingTrackParam); + while (trackParamAtCluster) { + + // reset track parameters and their covariances + trackParamAtCluster->SetParameters(startingTrackParam->GetParameters()); + trackParamAtCluster->SetZ(startingTrackParam->GetZ()); + + // extrapolation to the given z + AliMUONTrackExtrap::ExtrapToZ(trackParamAtCluster, trackParamAtCluster->GetClusterPtr()->GetZ()); + + // prepare next step + startingTrackParam = trackParamAtCluster; + trackParamAtCluster = (AliMUONTrackParam*) (fTrackParamAtCluster->After(trackParamAtCluster)); + } + } //__________________________________________________________________________ -void AliMUONTrack::SetFitStart(Int_t FitStart) +void AliMUONTrack::UpdateCovTrackParamAtCluster() { - // Set multiple Coulomb scattering option for track fit "fFitStart" - // from "FitStart" argument: 0 without, 1 with - if ((FitStart == 0) || (FitStart == 1)) fFitStart = FitStart; - // better implementation with enum(vertex, firstHit) ???? - else { - cout << "ERROR in AliMUONTrack::SetFitStart(FitStart)" << endl; - cout << "FitStart = " << FitStart << " is neither 0 nor 1" << endl; - exit(0); + /// Update track parameters and their covariances at each attached cluster + /// Include effects of multiple scattering in chambers + + if (GetNClusters() == 0) { + AliWarning("no cluster attached to the track"); + return; } - return; + + AliMUONTrackParam* startingTrackParam = (AliMUONTrackParam*) fTrackParamAtCluster->First(); + AliMUONTrackParam* trackParamAtCluster = (AliMUONTrackParam*) fTrackParamAtCluster->After(startingTrackParam); + Int_t expectedChamber = startingTrackParam->GetClusterPtr()->GetChamberId() + 1; + Int_t currentChamber; + while (trackParamAtCluster) { + + // reset track parameters and their covariances + trackParamAtCluster->SetParameters(startingTrackParam->GetParameters()); + trackParamAtCluster->SetZ(startingTrackParam->GetZ()); + trackParamAtCluster->SetCovariances(startingTrackParam->GetCovariances()); + + // add MCS effect + AliMUONTrackExtrap::AddMCSEffect(trackParamAtCluster,AliMUONConstants::ChamberThicknessInX0(),1.); + + // add MCS in missing chambers if any + currentChamber = trackParamAtCluster->GetClusterPtr()->GetChamberId(); + while (currentChamber > expectedChamber) { + // extrapolation to the missing chamber + AliMUONTrackExtrap::ExtrapToZCov(trackParamAtCluster, AliMUONConstants::DefaultChamberZ(expectedChamber)); + // add MCS effect + AliMUONTrackExtrap::AddMCSEffect(trackParamAtCluster,AliMUONConstants::ChamberThicknessInX0(),1.); + expectedChamber++; + } + + // extrapolation to the z of the current cluster + AliMUONTrackExtrap::ExtrapToZCov(trackParamAtCluster, trackParamAtCluster->GetClusterPtr()->GetZ()); + + // prepare next step + expectedChamber = currentChamber + 1; + startingTrackParam = trackParamAtCluster; + trackParamAtCluster = (AliMUONTrackParam*) (fTrackParamAtCluster->After(trackParamAtCluster)); + } + } //__________________________________________________________________________ -AliMUONTrackParam* AliMUONTrack::GetTrackParamAtFirstHit(void) { - // Get pointer to TrackParamAtFirstHit - return ((AliMUONTrackHit*) (fTrackHitsPtr->First()))->GetTrackParam();} +Bool_t AliMUONTrack::IsValid() +{ + /// check the validity of the current track (at least one cluster per station) + + Int_t nClusters = GetNClusters(); + AliMUONTrackParam *trackParam; + Int_t currentStation = 0, expectedStation = 0; + + for (Int_t i = 0; i < nClusters; i++) { + trackParam = (AliMUONTrackParam*) fTrackParamAtCluster->UncheckedAt(i); + + currentStation = trackParam->GetClusterPtr()->GetChamberId()/2; + + // missing station + if (currentStation > expectedStation) return kFALSE; + + // found station --> look for next one + if (currentStation == expectedStation) expectedStation++; + + } + + return currentStation == AliMUONConstants::NTrackingSt() - 1; + +} //__________________________________________________________________________ -void AliMUONTrack::RecursiveDump(void) -{ - // Recursive dump of AliMUONTrack, i.e. with dump of TrackHit's and HitForRec's - AliMUONTrackHit *trackHit; - AliMUONHitForRec *hitForRec; - cout << "Recursive dump of Track: " << this << endl; - // Track - this->Dump(); - for (Int_t trackHitIndex = 0; trackHitIndex < fNTrackHits; trackHitIndex++) { - trackHit = (AliMUONTrackHit*) ((*fTrackHitsPtr)[trackHitIndex]); - // TrackHit - cout << "TrackHit: " << trackHit << " (index: " << trackHitIndex << ")" << endl; - trackHit->Dump(); - hitForRec = trackHit->GetHitForRecPtr(); - // HitForRec - cout << "HitForRec: " << hitForRec << endl; - hitForRec->Dump(); +void AliMUONTrack::TagRemovableClusters() { + /// Identify clusters that can be removed from the track, + /// with the only requirement to have at least 1 cluster per station + + Int_t nClusters = GetNClusters(); + AliMUONTrackParam *trackParam, *nextTrackParam; + Int_t currentCh, nextCh; + + // reset flags to default + for (Int_t i = 0; i < nClusters; i++) { + trackParam = (AliMUONTrackParam*) fTrackParamAtCluster->UncheckedAt(i); + trackParam->SetRemovable(kFALSE); } - return; + + // loop over track parameters + for (Int_t i = 0; i < nClusters; i++) { + trackParam = (AliMUONTrackParam*) fTrackParamAtCluster->UncheckedAt(i); + + currentCh = trackParam->GetClusterPtr()->GetChamberId(); + + // loop over next track parameters + for (Int_t j = i+1; j < nClusters; j++) { + nextTrackParam = (AliMUONTrackParam*) fTrackParamAtCluster->UncheckedAt(j); + + nextCh = nextTrackParam->GetClusterPtr()->GetChamberId(); + + // check if the 2 clusters are on the same station + if (nextCh/2 != currentCh/2) break; + + // set clusters in the same station as being removable + trackParam->SetRemovable(kTRUE); + nextTrackParam->SetRemovable(kTRUE); + + } + + } + } //__________________________________________________________________________ -Int_t AliMUONTrack::HitsInCommon(AliMUONTrack* Track) +Bool_t AliMUONTrack::ComputeLocalChi2(Bool_t accountForMCS) { - // Returns the number of hits in common - // between the current track ("this") - // and the track pointed to by "Track". - Int_t hitsInCommon = 0; - AliMUONTrackHit *trackHit1, *trackHit2; - // Loop over hits of first track - trackHit1 = (AliMUONTrackHit*) this->GetTrackHitsPtr()->First(); - while (trackHit1) { - // Loop over hits of second track - trackHit2 = (AliMUONTrackHit*) Track->GetTrackHitsPtr()->First(); - while (trackHit2) { - // Increment "hitsInCommon" if both TrackHits point to the same HitForRec - if ( (trackHit1->GetHitForRecPtr()) == - (trackHit2->GetHitForRecPtr()) ) hitsInCommon++; - trackHit2 = (AliMUONTrackHit*) Track->GetTrackHitsPtr()->After(trackHit2); - } // trackHit2 - trackHit1 = (AliMUONTrackHit*) this->GetTrackHitsPtr()->After(trackHit1); - } // trackHit1 - return hitsInCommon; + /// Compute each cluster contribution to the chi2 of the track + /// accounting for multiple scattering or not according to the flag + /// - Also recompute the weight matrices of the attached clusters if accountForMCS=kTRUE + /// - Assume that track parameters at each cluster are corrects + /// - Return kFALSE if computation failed + AliDebug(1,"Enter ComputeLocalChi2"); + + if (!fTrackParamAtCluster) { + AliWarning("no cluster attached to this track"); + return kFALSE; + } + + if (accountForMCS) { // Compute local chi2 taking into account multiple scattering effects + + // Compute MCS covariance matrix only once + Int_t nClusters = GetNClusters(); + TMatrixD mcsCovariances(nClusters,nClusters); + ComputeMCSCovariances(mcsCovariances); + + // Make sure cluster weights are consistent with following calculations + if (!ComputeClusterWeights(&mcsCovariances)) { + AliWarning("cannot take into account the multiple scattering effects"); + return ComputeLocalChi2(kFALSE); + } + + // Compute chi2 of the track + Double_t globalChi2 = ComputeGlobalChi2(kTRUE); + if (globalChi2 < 0.) return kFALSE; + + // Loop over removable clusters and compute their local chi2 + AliMUONTrackParam* trackParamAtCluster1; + AliMUONVCluster *cluster, *discardedCluster; + Int_t iCluster1, iCluster2, iCurrentCluster1, iCurrentCluster2; + TMatrixD clusterWeightsNB(nClusters-1,nClusters-1); + TMatrixD clusterWeightsB(nClusters-1,nClusters-1); + Double_t *dX = new Double_t[nClusters-1]; + Double_t *dY = new Double_t[nClusters-1]; + Double_t globalChi2b; + AliMUONTrackParam* trackParamAtCluster = (AliMUONTrackParam*) fTrackParamAtCluster->First(); + while (trackParamAtCluster) { + + discardedCluster = trackParamAtCluster->GetClusterPtr(); + + // Recompute cluster weights without the current cluster + if (!ComputeClusterWeights(clusterWeightsNB, clusterWeightsB, &mcsCovariances, discardedCluster)) { + AliWarning("cannot take into account the multiple scattering effects"); + delete [] dX; + delete [] dY; + return ComputeLocalChi2(kFALSE); + } + + // Compute track chi2 without the current cluster + globalChi2b = 0.; + iCurrentCluster1 = 0; + for (iCluster1 = 0; iCluster1 < nClusters ; iCluster1++) { + trackParamAtCluster1 = (AliMUONTrackParam*) fTrackParamAtCluster->UncheckedAt(iCluster1); + cluster = trackParamAtCluster1->GetClusterPtr(); + + if (cluster == discardedCluster) continue; + + // Compute and save residuals + dX[iCurrentCluster1] = cluster->GetX() - trackParamAtCluster1->GetNonBendingCoor(); + dY[iCurrentCluster1] = cluster->GetY() - trackParamAtCluster1->GetBendingCoor(); + + iCurrentCluster2 = 0; + for (iCluster2 = 0; iCluster2 < iCluster1; iCluster2++) { + cluster = ((AliMUONTrackParam*) fTrackParamAtCluster->UncheckedAt(iCluster2))->GetClusterPtr(); + + if (cluster == discardedCluster) continue; + + // Add contribution from covariances + globalChi2b += (clusterWeightsNB(iCurrentCluster1, iCurrentCluster2) + + clusterWeightsNB(iCurrentCluster2, iCurrentCluster1)) * dX[iCurrentCluster1] * dX[iCurrentCluster2] + + (clusterWeightsB(iCurrentCluster1, iCurrentCluster2) + + clusterWeightsB(iCurrentCluster2, iCurrentCluster1)) * dY[iCurrentCluster1] * dY[iCurrentCluster2]; + + iCurrentCluster2++; + } + + // Add contribution from variances + globalChi2b += clusterWeightsNB(iCurrentCluster1, iCurrentCluster1) * dX[iCurrentCluster1] * dX[iCurrentCluster1] + + clusterWeightsB(iCurrentCluster1, iCurrentCluster1) * dY[iCurrentCluster1] * dY[iCurrentCluster1]; + + iCurrentCluster1++; + } + + // Set local chi2 + trackParamAtCluster->SetLocalChi2(globalChi2 - globalChi2b); + + trackParamAtCluster = (AliMUONTrackParam*) fTrackParamAtCluster->After(trackParamAtCluster); + } + + delete [] dX; + delete [] dY; + + } else { // without multiple scattering effects + + AliMUONVCluster *discardedCluster; + Double_t dX, dY; + AliMUONTrackParam* trackParamAtCluster = (AliMUONTrackParam*) fTrackParamAtCluster->First(); + while (trackParamAtCluster) { + + discardedCluster = trackParamAtCluster->GetClusterPtr(); + + // Compute residuals + dX = discardedCluster->GetX() - trackParamAtCluster->GetNonBendingCoor(); + dY = discardedCluster->GetY() - trackParamAtCluster->GetBendingCoor(); + + // Set local chi2 + trackParamAtCluster->SetLocalChi2(dX * dX / discardedCluster->GetErrX2() + dY * dY / discardedCluster->GetErrY2()); + + trackParamAtCluster = (AliMUONTrackParam*) fTrackParamAtCluster->After(trackParamAtCluster); + } + + } + + return kTRUE; + } //__________________________________________________________________________ -void AliMUONTrack::Fit() +Double_t AliMUONTrack::ComputeGlobalChi2(Bool_t accountForMCS) { - // Fit the current track ("this"), - // with or without multiple Coulomb scattering according to "fFitMCS", - // with the number of parameters given by "fFitNParam" - // (3 if one keeps X and Y fixed in "TrackParam", 5 if one lets them vary), - // starting, according to "fFitStart", - // with track parameters at vertex or at the first TrackHit. - // "fFitMCS", "fFitNParam" and "fFitStart" have to be set before - // by calling the corresponding Set methods. - Double_t arg[1], benC, errorParam, invBenP, lower, nonBenC, upper, x, y; - char parName[50]; - AliMUONTrackParam *trackParam; - // Check if Minuit is initialized... - fgFitter = TVirtualFitter::Fitter(this); // add 3 or 5 for the maximum number of parameters ??? - fgFitter->Clear(); // necessary ???? probably yes - // how to reset the printout number at every fit ???? - // is there any risk to leave it like that ???? - // how to go faster ???? choice of Minuit parameters like EDM ???? - // choice of function to be minimized according to fFitMCS - if (fFitMCS == 0) fgFitter->SetFCN(TrackChi2); - else fgFitter->SetFCN(TrackChi2MCS); - arg[0] = 1; - fgFitter->ExecuteCommand("SET PRINT", arg, 1); // More printing !!!! - // Parameters according to "fFitStart" - // (should be a function to be used at every place where needed ????) - if (fFitStart == 0) trackParam = &fTrackParamAtVertex; - else trackParam = this->GetTrackParamAtFirstHit(); - // set first 3 Minuit parameters - // could be tried with no limits for the search (min=max=0) ???? - fgFitter->SetParameter(0, "InvBenP", - trackParam->GetInverseBendingMomentum(), - 0.003, -0.4, 0.4); - fgFitter->SetParameter(1, "BenS", - trackParam->GetBendingSlope(), - 0.001, -0.5, 0.5); - fgFitter->SetParameter(2, "NonBenS", - trackParam->GetNonBendingSlope(), - 0.001, -0.5, 0.5); - if (fFitNParam == 5) { - // set last 2 Minuit parameters (no limits for the search: min=max=0) - fgFitter->SetParameter(3, "X", - trackParam->GetNonBendingCoor(), - 0.03, 0.0, 0.0); - fgFitter->SetParameter(4, "Y", - trackParam->GetBendingCoor(), - 0.10, 0.0, 0.0); - } - // search without gradient calculation in the function - fgFitter->ExecuteCommand("SET NOGRADIENT", arg, 0); - // minimization - fgFitter->ExecuteCommand("MINIMIZE", arg, 0); - // exit from Minuit - fgFitter->ExecuteCommand("EXIT", arg, 0); // necessary ???? - // get results into "invBenP", "benC", "nonBenC" ("x", "y") - fgFitter->GetParameter(0, parName, invBenP, errorParam, lower, upper); - fgFitter->GetParameter(1, parName, benC, errorParam, lower, upper); - fgFitter->GetParameter(2, parName, nonBenC, errorParam, lower, upper); - if (fFitNParam == 5) { - fgFitter->GetParameter(3, parName, x, errorParam, lower, upper); - fgFitter->GetParameter(4, parName, y, errorParam, lower, upper); + /// Compute the chi2 of the track accounting for multiple scattering or not according to the flag + /// - Assume that track parameters at each cluster are corrects + /// - Assume the cluster weights matrices are corrects + /// - Return negative value if chi2 computation failed + AliDebug(1,"Enter ComputeGlobalChi2"); + + if (!fTrackParamAtCluster) { + AliWarning("no cluster attached to this track"); + return 1.e10; } - // result of the fit into track parameters - trackParam->SetInverseBendingMomentum(invBenP); - trackParam->SetBendingSlope(benC); - trackParam->SetNonBendingSlope(nonBenC); - if (fFitNParam == 5) { - trackParam->SetNonBendingCoor(x); - trackParam->SetBendingCoor(y); + + Double_t chi2 = 0.; + + if (accountForMCS) { + + // Check the weight matrices. If weight matrices are not available compute chi2 without MCS + if (!fClusterWeightsNonBending || !fClusterWeightsBending) { + AliWarning("cluster weights including multiple scattering effects are not available\n\t\t --> compute chi2 WITHOUT multiple scattering"); + return ComputeGlobalChi2(kFALSE); + } + Int_t nClusters = GetNClusters(); + if (fClusterWeightsNonBending->GetNrows() != nClusters || fClusterWeightsBending->GetNcols() != nClusters) { + AliWarning("cluster weights including multiple scattering effects are not available\n\t\t --> compute chi2 WITHOUT multiple scattering"); + return ComputeGlobalChi2(kFALSE); + } + + // Compute chi2 + AliMUONVCluster *cluster; + Double_t *dX = new Double_t[nClusters]; + Double_t *dY = new Double_t[nClusters]; + AliMUONTrackParam* trackParamAtCluster; + for (Int_t iCluster1 = 0; iCluster1 < nClusters; iCluster1++) { + trackParamAtCluster = (AliMUONTrackParam*) fTrackParamAtCluster->UncheckedAt(iCluster1); + cluster = trackParamAtCluster->GetClusterPtr(); + dX[iCluster1] = cluster->GetX() - trackParamAtCluster->GetNonBendingCoor(); + dY[iCluster1] = cluster->GetY() - trackParamAtCluster->GetBendingCoor(); + for (Int_t iCluster2 = 0; iCluster2 < iCluster1; iCluster2++) { + chi2 += ((*fClusterWeightsNonBending)(iCluster1, iCluster2) + (*fClusterWeightsNonBending)(iCluster2, iCluster1)) * dX[iCluster1] * dX[iCluster2] + + ((*fClusterWeightsBending)(iCluster1, iCluster2) + (*fClusterWeightsBending)(iCluster2, iCluster1)) * dY[iCluster1] * dY[iCluster2]; + } + chi2 += ((*fClusterWeightsNonBending)(iCluster1, iCluster1) * dX[iCluster1] * dX[iCluster1]) + + ((*fClusterWeightsBending)(iCluster1, iCluster1) * dY[iCluster1] * dY[iCluster1]); + } + delete [] dX; + delete [] dY; + + } else { + + AliMUONVCluster *cluster; + Double_t dX, dY; + AliMUONTrackParam* trackParamAtCluster; + Int_t nClusters = GetNClusters(); + for (Int_t iCluster = 0; iCluster < nClusters ; iCluster++) { + trackParamAtCluster = (AliMUONTrackParam*) fTrackParamAtCluster->UncheckedAt(iCluster); + cluster = trackParamAtCluster->GetClusterPtr(); + dX = cluster->GetX() - trackParamAtCluster->GetNonBendingCoor(); + dY = cluster->GetY() - trackParamAtCluster->GetBendingCoor(); + chi2 += dX * dX / cluster->GetErrX2() + dY * dY / cluster->GetErrY2(); + } + } - // global result of the fit - Double_t fedm, errdef; - Int_t npari, nparx; - fgFitter->GetStats(fFitFMin, fedm, errdef, npari, nparx); + + return chi2; + } //__________________________________________________________________________ -void AliMUONTrack::AddSegment(AliMUONSegment* Segment) +Bool_t AliMUONTrack::ComputeClusterWeights(TMatrixD* mcsCovariances) { - // Add Segment to the track - AddHitForRec(Segment->GetHitForRec1()); // 1st hit - AddHitForRec(Segment->GetHitForRec2()); // 2nd hit + /// Compute the weight matrices of the attached clusters, in non bending and bending direction, + /// accounting for multiple scattering correlations and cluster resolution + /// - Use the provided MCS covariance matrix if any (otherwise build it temporarily) + /// - Assume that track parameters at each cluster are corrects + /// - Return kFALSE if computation failed + AliDebug(1,"Enter ComputeClusterWeights1"); + + if (!fTrackParamAtCluster) { + AliWarning("no cluster attached to this track"); + return kFALSE; + } + + // Alocate memory + Int_t nClusters = GetNClusters(); + if (!fClusterWeightsNonBending) fClusterWeightsNonBending = new TMatrixD(nClusters,nClusters); + if (!fClusterWeightsBending) fClusterWeightsBending = new TMatrixD(nClusters,nClusters); + + // Compute weights matrices + if (!ComputeClusterWeights(*fClusterWeightsNonBending, *fClusterWeightsBending, mcsCovariances)) return kFALSE; + + return kTRUE; + } //__________________________________________________________________________ -void AliMUONTrack::AddHitForRec(AliMUONHitForRec* HitForRec) +Bool_t AliMUONTrack::ComputeClusterWeights(TMatrixD& clusterWeightsNB, TMatrixD& clusterWeightsB, + TMatrixD* mcsCovariances, AliMUONVCluster* discardedCluster) const { - // Add HitForRec to the track: - // actual TrackHit into TClonesArray of TrackHit's for the event; - // pointer to actual TrackHit in TObjArray of pointers to TrackHit's for the track - TClonesArray *recTrackHitsPtr = this->fEventReconstructor->GetRecTrackHitsPtr(); - Int_t eventTrackHits = this->fEventReconstructor->GetNRecTrackHits(); - // event - AliMUONTrackHit* trackHit = - new ((*recTrackHitsPtr)[eventTrackHits]) AliMUONTrackHit(HitForRec); - this->fEventReconstructor->SetNRecTrackHits(eventTrackHits + 1); - // track - fTrackHitsPtr->Add(trackHit); - fNTrackHits++; + /// Compute the weight matrices, in non bending and bending direction, + /// of the other attached clusters assuming the discarded one does not exist + /// accounting for multiple scattering correlations and cluster resolution + /// - Use the provided MCS covariance matrix if any (otherwise build it temporarily) + /// - Return kFALSE if computation failed + AliDebug(1,"Enter ComputeClusterWeights2"); + + // Check MCS covariance matrix and recompute it if need + Int_t nClusters = GetNClusters(); + Bool_t deleteMCSCov = kFALSE; + if (!mcsCovariances) { + mcsCovariances = new TMatrixD(nClusters,nClusters); + deleteMCSCov = kTRUE; + ComputeMCSCovariances(*mcsCovariances); + } + + // Resize the weights matrices; alocate memory + if (discardedCluster) { + clusterWeightsNB.ResizeTo(nClusters-1,nClusters-1); + clusterWeightsB.ResizeTo(nClusters-1,nClusters-1); + } else { + clusterWeightsNB.ResizeTo(nClusters,nClusters); + clusterWeightsB.ResizeTo(nClusters,nClusters); + } + + // Define variables + AliMUONVCluster *cluster1, *cluster2; + Int_t iCurrentCluster1, iCurrentCluster2; + + // Compute the covariance matrices + iCurrentCluster1 = 0; + for (Int_t iCluster1 = 0; iCluster1 < nClusters; iCluster1++) { + cluster1 = ((AliMUONTrackParam*) fTrackParamAtCluster->UncheckedAt(iCluster1))->GetClusterPtr(); + + if (cluster1 == discardedCluster) continue; + + // Loop over next clusters + iCurrentCluster2 = iCurrentCluster1; + for (Int_t iCluster2 = iCluster1; iCluster2 < nClusters; iCluster2++) { + cluster2 = ((AliMUONTrackParam*) fTrackParamAtCluster->UncheckedAt(iCluster2))->GetClusterPtr(); + + if (cluster2 == discardedCluster) continue; + + // Fill with MCS covariances + clusterWeightsNB(iCurrentCluster1, iCurrentCluster2) = (*mcsCovariances)(iCluster1,iCluster2); + + // Equal contribution from multiple scattering in non bending and bending directions + clusterWeightsB(iCurrentCluster1, iCurrentCluster2) = clusterWeightsNB(iCurrentCluster1, iCurrentCluster2); + + // Add contribution from cluster resolution to diagonal element and symmetrize the matrix + if (iCurrentCluster1 == iCurrentCluster2) { + + // In non bending plane + clusterWeightsNB(iCurrentCluster1, iCurrentCluster1) += cluster1->GetErrX2(); + // In bending plane + clusterWeightsB(iCurrentCluster1, iCurrentCluster1) += cluster1->GetErrY2(); + + } else { + + // In non bending plane + clusterWeightsNB(iCurrentCluster2, iCurrentCluster1) = clusterWeightsNB(iCurrentCluster1, iCurrentCluster2); + // In bending plane + clusterWeightsB(iCurrentCluster2, iCurrentCluster1) = clusterWeightsB(iCurrentCluster1, iCurrentCluster2); + + } + + iCurrentCluster2++; + } + + iCurrentCluster1++; + } + + // Inversion of covariance matrices to get the weights + if (clusterWeightsNB.Determinant() != 0 && clusterWeightsB.Determinant() != 0) { + clusterWeightsNB.Invert(); + clusterWeightsB.Invert(); + } else { + AliWarning(" Determinant = 0"); + clusterWeightsNB.ResizeTo(0,0); + clusterWeightsB.ResizeTo(0,0); + if(deleteMCSCov) delete mcsCovariances; + return kFALSE; + } + + if(deleteMCSCov) delete mcsCovariances; + + return kTRUE; + } //__________________________________________________________________________ -void AliMUONTrack::SetTrackParamAtHit(Int_t indexHit, AliMUONTrackParam *TrackParam) +void AliMUONTrack::ComputeMCSCovariances(TMatrixD& mcsCovariances) const { - // Set track parameters at TrackHit with index "indexHit" - // from the track parameters pointed to by "TrackParam". - AliMUONTrackHit* trackHit = (AliMUONTrackHit*) ((*fTrackHitsPtr)[indexHit]); - trackHit->SetTrackParam(TrackParam); + /// Compute the multiple scattering covariance matrix + /// (assume that track parameters at each cluster are corrects) + AliDebug(1,"Enter ComputeMCSCovariances"); + + // Reset the size of the covariance matrix if needed + Int_t nClusters = GetNClusters(); + if (mcsCovariances.GetNrows() != nClusters) mcsCovariances.ResizeTo(nClusters,nClusters); + + // Define variables + Int_t nChambers = AliMUONConstants::NTrackingCh(); + AliMUONTrackParam* trackParamAtCluster; + AliMUONTrackParam extrapTrackParam; + Int_t currentChamber = 0, expectedChamber = 0, size = 0; + Double_t *mcsAngle2 = new Double_t[2*nChambers]; + Double_t *zMCS = new Double_t[2*nChambers]; + Int_t *indices = new Int_t[2*nClusters]; + + // Compute multiple scattering dispersion angle at each chamber + // and save the z position where it is calculated + for (Int_t iCluster = 0; iCluster < nClusters; iCluster++) { + trackParamAtCluster = (AliMUONTrackParam*) fTrackParamAtCluster->UncheckedAt(iCluster); + + // look for missing chambers if any + currentChamber = trackParamAtCluster->GetClusterPtr()->GetChamberId(); + while (currentChamber > expectedChamber) { + + // Save the z position where MCS dispersion is calculated + zMCS[size] = AliMUONConstants::DefaultChamberZ(expectedChamber); + + // Do not take into account MCS in chambers prior the first cluster + if (iCluster > 0) { + + // Get track parameters at missing chamber z + extrapTrackParam = *trackParamAtCluster; + AliMUONTrackExtrap::ExtrapToZ(&extrapTrackParam, zMCS[size]); + + // Save multiple scattering dispersion angle in missing chamber + mcsAngle2[size] = AliMUONTrackExtrap::GetMCSAngle2(extrapTrackParam,AliMUONConstants::ChamberThicknessInX0(),1.); + + } else mcsAngle2[size] = 0.; + + expectedChamber++; + size++; + } + + // Save z position where MCS dispersion is calculated + zMCS[size] = trackParamAtCluster->GetZ(); + + // Save multiple scattering dispersion angle in current chamber + mcsAngle2[size] = AliMUONTrackExtrap::GetMCSAngle2(*trackParamAtCluster,AliMUONConstants::ChamberThicknessInX0(),1.); + + // Save indice in zMCS array corresponding to the current cluster + indices[iCluster] = size; + + expectedChamber = currentChamber + 1; + size++; + } + + // complete array of z if last cluster is on the last but one chamber + if (currentChamber != nChambers-1) zMCS[size++] = AliMUONConstants::DefaultChamberZ(nChambers-1); + + // Compute the covariance matrix + for (Int_t iCluster1 = 0; iCluster1 < nClusters; iCluster1++) { + + for (Int_t iCluster2 = iCluster1; iCluster2 < nClusters; iCluster2++) { + + // Initialization to 0 (diagonal plus upper triangular part) + mcsCovariances(iCluster1,iCluster2) = 0.; + + // Compute contribution from multiple scattering in upstream chambers + for (Int_t k = 0; k < indices[iCluster1]; k++) { + mcsCovariances(iCluster1,iCluster2) += (zMCS[indices[iCluster1]] - zMCS[k]) * (zMCS[indices[iCluster2]] - zMCS[k]) * mcsAngle2[k]; + } + + // Symetrize the matrix + mcsCovariances(iCluster2,iCluster1) = mcsCovariances(iCluster1,iCluster2); + } + + } + + delete [] mcsAngle2; + delete [] zMCS; + delete [] indices; + } //__________________________________________________________________________ -void AliMUONTrack::SetTrackParamAtVertex() +Int_t AliMUONTrack::ClustersInCommon(AliMUONTrack* track) const { - // Set track parameters at vertex. - // TrackHit's are assumed to be only in stations(1..) 4 and 5, - // and sorted according to increasing Z.. - // Parameters are calculated from information in HitForRec's - // of first and last TrackHit's. - AliMUONTrackParam *trackParam = - &fTrackParamAtVertex; // pointer to track parameters - // Pointer to HitForRec of first TrackHit - AliMUONHitForRec *firstHit = - ((AliMUONTrackHit*) (fTrackHitsPtr->First()))->GetHitForRecPtr(); - // Pointer to HitForRec of last TrackHit - AliMUONHitForRec *lastHit = - ((AliMUONTrackHit*) (fTrackHitsPtr->Last()))->GetHitForRecPtr(); - // Z difference between first and last hits - Double_t deltaZ = firstHit->GetZ() - lastHit->GetZ(); - // bending slope in stations(1..) 4 and 5 - Double_t bendingSlope = - (firstHit->GetBendingCoor() - lastHit->GetBendingCoor()) / deltaZ; - trackParam->SetBendingSlope(bendingSlope); - // impact parameter - Double_t impactParam = - firstHit->GetBendingCoor() - bendingSlope * firstHit->GetZ(); // same if from firstHit and lastHit ???? - // signed bending momentum - Double_t signedBendingMomentum = - fEventReconstructor->GetBendingMomentumFromImpactParam(impactParam); - trackParam->SetInverseBendingMomentum(1.0 / signedBendingMomentum); - // bending slope at vertex - trackParam-> - SetBendingSlope(bendingSlope + - impactParam / fEventReconstructor->GetSimpleBPosition()); - // non bending slope - Double_t nonBendingSlope = - (firstHit->GetNonBendingCoor() - lastHit->GetNonBendingCoor()) / deltaZ; - trackParam->SetNonBendingSlope(nonBendingSlope); - // vertex coordinates at (0,0,0) - trackParam->SetZ(0.0); - trackParam->SetBendingCoor(0.0); - trackParam->SetNonBendingCoor(0.0); + /// Returns the number of clusters in common between the current track ("this") + /// and the track pointed to by "track". + if (!fTrackParamAtCluster || !this->fTrackParamAtCluster) return 0; + Int_t clustersInCommon = 0; + AliMUONTrackParam *trackParamAtCluster1, *trackParamAtCluster2; + // Loop over clusters of first track + trackParamAtCluster1 = (AliMUONTrackParam*) this->fTrackParamAtCluster->First(); + while (trackParamAtCluster1) { + // Loop over clusters of second track + trackParamAtCluster2 = (AliMUONTrackParam*) track->fTrackParamAtCluster->First(); + while (trackParamAtCluster2) { + // Increment "clustersInCommon" if both trackParamAtCluster1 & 2 point to the same cluster + if ((trackParamAtCluster1->GetClusterPtr()) == (trackParamAtCluster2->GetClusterPtr())) { + clustersInCommon++; + break; + } + trackParamAtCluster2 = (AliMUONTrackParam*) track->fTrackParamAtCluster->After(trackParamAtCluster2); + } // trackParamAtCluster2 + trackParamAtCluster1 = (AliMUONTrackParam*) this->fTrackParamAtCluster->After(trackParamAtCluster1); + } // trackParamAtCluster1 + return clustersInCommon; } //__________________________________________________________________________ -void TrackChi2(Int_t &NParam, Double_t *Gradient, Double_t &Chi2, Double_t *Param, Int_t Flag) +Double_t AliMUONTrack::GetNormalizedChi2() const { - // Return the "Chi2" to be minimized with Minuit for track fitting, - // with "NParam" parameters - // and their current values in array pointed to by "Param". - // Assumes that the track hits are sorted according to increasing Z. - // Track parameters at each TrackHit are updated accordingly. - // Multiple Coulomb scattering is not taken into account - AliMUONTrack *trackBeingFitted; - AliMUONTrackHit* hit; - AliMUONTrackParam param1; - Int_t hitNumber; - Double_t zHit; - Chi2 = 0.0; // initialize Chi2 - // copy of track parameters to be fitted - trackBeingFitted = (AliMUONTrack*) AliMUONTrack::Fitter()->GetObjectFit(); - if (trackBeingFitted->GetFitStart() == 0) - param1 = *(trackBeingFitted->GetTrackParamAtVertex()); - else param1 = *(trackBeingFitted->GetTrackParamAtFirstHit()); - // Minuit parameters to be fitted into this copy - param1.SetInverseBendingMomentum(Param[0]); - param1.SetBendingSlope(Param[1]); - param1.SetNonBendingSlope(Param[2]); - if (NParam == 5) { - param1.SetNonBendingCoor(Param[3]); - param1.SetBendingCoor(Param[4]); - } - // Follow track through all planes of track hits - for (hitNumber = 0; hitNumber < trackBeingFitted->GetNTrackHits(); hitNumber++) { - hit = (AliMUONTrackHit*) (*(trackBeingFitted->GetTrackHitsPtr()))[hitNumber]; - zHit = hit->GetHitForRecPtr()->GetZ(); - // do something special if 2 hits with same Z ???? - // security against infinite loop ???? - (¶m1)->ExtrapToZ(zHit); // extrapolation - hit->SetTrackParam(¶m1); - // Increment Chi2 - // done hit per hit, with hit resolution, - // and not with point and angle like in "reco_muon.F" !!!! - // Needs to add multiple scattering contribution ???? - Double_t dX = - hit->GetHitForRecPtr()->GetNonBendingCoor() - (¶m1)->GetNonBendingCoor(); - Double_t dY = - hit->GetHitForRecPtr()->GetBendingCoor() - (¶m1)->GetBendingCoor(); - Chi2 = - Chi2 + - dX * dX / hit->GetHitForRecPtr()->GetNonBendingReso2() + - dY * dY / hit->GetHitForRecPtr()->GetBendingReso2(); - } + /// return the chi2 value divided by the number of degrees of freedom (or 1.e10 if ndf < 0) + + Double_t numberOfDegFree = (2. * GetNClusters() - 5.); + if (numberOfDegFree > 0.) return fGlobalChi2 / numberOfDegFree; + else return 1.e10; } //__________________________________________________________________________ -void TrackChi2MCS(Int_t &NParam, Double_t *Gradient, Double_t &Chi2, Double_t *Param, Int_t Flag) +Bool_t* AliMUONTrack::CompatibleTrack(AliMUONTrack *track, Double_t sigmaCut) const { - // Return the "Chi2" to be minimized with Minuit for track fitting, - // with "NParam" parameters - // and their current values in array pointed to by "Param". - // Assumes that the track hits are sorted according to increasing Z. - // Track parameters at each TrackHit are updated accordingly. - // Multiple Coulomb scattering is taken into account with covariance matrix. - AliMUONTrack *trackBeingFitted; - AliMUONTrackParam param1; - Chi2 = 0.0; // initialize Chi2 - // copy of track parameters to be fitted - trackBeingFitted = (AliMUONTrack*) AliMUONTrack::Fitter()->GetObjectFit(); - if (trackBeingFitted->GetFitStart() == 0) - param1 = *(trackBeingFitted->GetTrackParamAtVertex()); - else param1 = *(trackBeingFitted->GetTrackParamAtFirstHit()); - // Minuit parameters to be fitted into this copy - param1.SetInverseBendingMomentum(Param[0]); - param1.SetBendingSlope(Param[1]); - param1.SetNonBendingSlope(Param[2]); - if (NParam == 5) { - param1.SetNonBendingCoor(Param[3]); - param1.SetBendingCoor(Param[4]); - } + /// for each chamber: return kTRUE (kFALSE) if clusters are compatible (not compatible) + AliMUONTrackParam *trackParamAtCluster1, *trackParamAtCluster2; + AliMUONVCluster *cluster1, *cluster2; + Double_t chi2, dX, dY, dZ; + Double_t chi2Max = sigmaCut * sigmaCut; + Double_t dZMax = 1.; // 1 cm + + Bool_t *compatibleCluster = new Bool_t[AliMUONConstants::NTrackingCh()]; + for ( Int_t ch = 0; ch < AliMUONConstants::NTrackingCh(); ch++) compatibleCluster[ch] = kFALSE; - AliMUONTrackHit *hit; - Bool_t goodDeterminant; - Int_t hitNumber, hitNumber1, hitNumber2, hitNumber3; - Double_t zHit[10], paramBendingCoor[10], paramNonBendingCoor[10], ap[10]; - Double_t hitBendingCoor[10], hitNonBendingCoor[10]; - Double_t hitBendingReso2[10], hitNonBendingReso2[10]; - // dimension 10 in parameter ??? related to AliMUONConstants::NTrackingCh() !!!! - Int_t numberOfHit = TMath::Min(trackBeingFitted->GetNTrackHits(), 10); - TMatrix *covBending = new TMatrix(numberOfHit, numberOfHit); - TMatrix *covNonBending = new TMatrix(numberOfHit, numberOfHit); - - // Predicted coordinates and multiple scattering angles are first calculated - for (hitNumber = 0; hitNumber < numberOfHit; hitNumber++) { - hit = (AliMUONTrackHit*) (*(trackBeingFitted->GetTrackHitsPtr()))[hitNumber]; - zHit[hitNumber] = hit->GetHitForRecPtr()->GetZ(); - // do something special if 2 hits with same Z ???? - // security against infinite loop ???? - (¶m1)->ExtrapToZ(zHit[hitNumber]); // extrapolation - hit->SetTrackParam(¶m1); - paramBendingCoor[hitNumber] = (¶m1)->GetBendingCoor(); - paramNonBendingCoor[hitNumber] = (¶m1)->GetNonBendingCoor(); - hitBendingCoor[hitNumber] = hit->GetHitForRecPtr()->GetBendingCoor(); - hitNonBendingCoor[hitNumber] = hit->GetHitForRecPtr()->GetNonBendingCoor(); - hitBendingReso2[hitNumber] = hit->GetHitForRecPtr()->GetBendingReso2(); - hitNonBendingReso2[hitNumber] = hit->GetHitForRecPtr()->GetNonBendingReso2(); - ap[hitNumber] = MultipleScatteringAngle2(hit); // multiple scatt. angle ^2 + if (!fTrackParamAtCluster || !this->fTrackParamAtCluster) return compatibleCluster; + + // Loop over clusters of first track + trackParamAtCluster1 = (AliMUONTrackParam*) this->fTrackParamAtCluster->First(); + while (trackParamAtCluster1) { + + cluster1 = trackParamAtCluster1->GetClusterPtr(); + + // Loop over clusters of second track + trackParamAtCluster2 = (AliMUONTrackParam*) track->fTrackParamAtCluster->First(); + while (trackParamAtCluster2) { + + cluster2 = trackParamAtCluster2->GetClusterPtr(); + + //prepare next step + trackParamAtCluster2 = (AliMUONTrackParam*) track->fTrackParamAtCluster->After(trackParamAtCluster2); + + // z direction + dZ = cluster1->GetZ() - cluster2->GetZ(); + if (dZ > dZMax) continue; + + // non bending direction + dX = cluster1->GetX() - cluster2->GetX(); + chi2 = dX * dX / (cluster1->GetErrX2() + cluster2->GetErrX2()); + if (chi2 > chi2Max) continue; + + // bending direction + dY = cluster1->GetY() - cluster2->GetY(); + chi2 = dY * dY / (cluster1->GetErrY2() + cluster2->GetErrY2()); + if (chi2 > chi2Max) continue; + + compatibleCluster[cluster1->GetChamberId()] = kTRUE; + break; + } + + trackParamAtCluster1 = (AliMUONTrackParam*) this->fTrackParamAtCluster->After(trackParamAtCluster1); } + + return compatibleCluster; +} - // Calculates the covariance matrix - // One chamber is taken into account between successive hits. - // "ap" should be changed for taking into account the eventual missing hits - // by defining an "equivalent" chamber thickness !!!! - for (hitNumber1 = 0; hitNumber1 < numberOfHit; hitNumber1++) { - for (hitNumber2 = hitNumber1; hitNumber2 < numberOfHit; hitNumber2++) { - // initialization to 0 (diagonal plus upper triangular part) - (*covBending)(hitNumber2, hitNumber1) = 0.0; - // contribution from multiple scattering in bending plane: - // loop over upstream hits - for (hitNumber3 = 0; hitNumber3 < hitNumber1; hitNumber3++) { - (*covBending)(hitNumber2, hitNumber1) = - (*covBending)(hitNumber2, hitNumber1) + - ((zHit[hitNumber1] - zHit[hitNumber3]) * - (zHit[hitNumber2] - zHit[hitNumber3]) * ap[hitNumber3]); - } - // equal contribution from multiple scattering in non bending plane - (*covNonBending)(hitNumber2, hitNumber1) = - (*covBending)(hitNumber2, hitNumber1); - if (hitNumber1 == hitNumber2) { - // Diagonal elements: add contribution from position measurements - // in bending plane - (*covBending)(hitNumber2, hitNumber1) = - (*covBending)(hitNumber2, hitNumber1) + hitBendingReso2[hitNumber1]; - // and in non bending plane - (*covNonBending)(hitNumber2, hitNumber1) = - (*covNonBending)(hitNumber2, hitNumber1) + hitNonBendingReso2[hitNumber1]; - } - else { - // Non diagonal elements: symmetrization - // for bending plane - (*covBending)(hitNumber1, hitNumber2) = - (*covBending)(hitNumber2, hitNumber1); - // and non bending plane - (*covNonBending)(hitNumber1, hitNumber2) = - (*covNonBending)(hitNumber2, hitNumber1); - } - } // for (hitNumber2 = hitNumber1;... - } // for (hitNumber1 = 0;... - - // Inverts covariance matrix - goodDeterminant = kTRUE; - // check whether the Invert method returns flag if matrix cannot be inverted, - // and do not calculate the Determinant in that case !!!! - if (covBending->Determinant() != 0) { - covBending->Invert(); - } else { - goodDeterminant = kFALSE; - cout << "Warning in ChiMCS Determinant Bending=0: " << endl; - } - if (covNonBending->Determinant() != 0) { - covNonBending->Invert(); - } else { - goodDeterminant = kFALSE; - cout << "Warning in ChiMCS Determinant non Bending=0: " << endl; - } +//__________________________________________________________________________ +void AliMUONTrack::SetTrackParamAtVertex(const AliMUONTrackParam* trackParam) +{ + /// set track parameters at vertex + if (trackParam == 0x0) return; + if (fTrackParamAtVertex) *fTrackParamAtVertex = *trackParam; + else fTrackParamAtVertex = new AliMUONTrackParam(*trackParam); +} - // It would be worth trying to calculate the inverse of the covariance matrix - // only once per fit, since it cannot change much in principle, - // and it would save a lot of computing time !!!! - - // Calculates Chi2 - if (goodDeterminant) { // with Multiple Scattering if inversion correct - for (hitNumber1=0; hitNumber1 < numberOfHit ; hitNumber1++){ - for (hitNumber2=0; hitNumber2 < numberOfHit; hitNumber2++){ - Chi2 = Chi2 + - ((*covBending)(hitNumber2, hitNumber1) * - (hitBendingCoor[hitNumber1] - paramBendingCoor[hitNumber1]) * - (hitBendingCoor[hitNumber2] - paramBendingCoor[hitNumber2])); - Chi2 = Chi2 + - ((*covNonBending)(hitNumber2, hitNumber1) * - (hitNonBendingCoor[hitNumber1] - paramNonBendingCoor[hitNumber1]) * - (hitNonBendingCoor[hitNumber2] - paramNonBendingCoor[hitNumber2])); - } - } - } else { // without Multiple Scattering if inversion impossible - for (hitNumber1=0; hitNumber1 < numberOfHit ; hitNumber1++) { - Chi2 = Chi2 + - ((hitBendingCoor[hitNumber1] - paramBendingCoor[hitNumber1]) * - (hitBendingCoor[hitNumber1] - paramBendingCoor[hitNumber1]) / - hitBendingReso2[hitNumber1]); - Chi2 = Chi2 + - ((hitNonBendingCoor[hitNumber1] - paramNonBendingCoor[hitNumber1]) * - (hitNonBendingCoor[hitNumber1] - paramNonBendingCoor[hitNumber1]) / - hitNonBendingReso2[hitNumber1]); - } +//__________________________________________________________________________ +void AliMUONTrack::RecursiveDump() const +{ + /// Recursive dump of AliMUONTrack, i.e. with dump of trackParamAtCluster and attached clusters + AliMUONTrackParam *trackParamAtCluster; + AliMUONVCluster *cluster; + cout << "Recursive dump of Track: " << this << endl; + // Track + this->Dump(); + for (Int_t iCluster = 0; iCluster < GetNClusters(); iCluster++) { + trackParamAtCluster = (AliMUONTrackParam*) ((*fTrackParamAtCluster)[iCluster]); + // trackParamAtCluster + cout << "trackParamAtCluster: " << trackParamAtCluster << " (index: " << iCluster << ")" << endl; + trackParamAtCluster->Dump(); + cluster = trackParamAtCluster->GetClusterPtr(); + // cluster + cout << "cluster: " << cluster << endl; + cluster->Print(); } + return; +} - delete covBending; - delete covNonBending; +//_____________________________________________- +void AliMUONTrack::Print(Option_t*) const +{ + /// Printing Track information + + cout << " No.Clusters=" << setw(2) << GetNClusters() << + ", Match2Trig=" << setw(1) << GetMatchTrigger() << + ", LoTrgNum=" << setw(3) << GetLoTrgNum() << + ", Chi2-tracking-trigger=" << setw(8) << setprecision(5) << GetChi2MatchTrigger(); + cout << Form(" HitTriggerPattern %x",fHitsPatternInTrigCh) << endl; + if (fTrackParamAtCluster) fTrackParamAtCluster->First()->Print("FULL"); } -Double_t MultipleScatteringAngle2(AliMUONTrackHit *TrackHit) +//__________________________________________________________________________ +void AliMUONTrack::SetLocalTrigger(Int_t loCirc, Int_t loStripX, Int_t loStripY, Int_t loDev, Int_t loLpt, Int_t loHpt) { - // Returns square of multiple Coulomb scattering angle - // at TrackHit pointed to by "TrackHit" - Double_t slopeBending, slopeNonBending, radiationLength, inverseBendingMomentum2, inverseTotalMomentum2; - Double_t varMultipleScatteringAngle; - AliMUONTrack *trackBeingFitted = (AliMUONTrack*) AliMUONTrack::Fitter()->GetObjectFit(); - AliMUONTrackParam *param = TrackHit->GetTrackParam(); - // Better implementation in AliMUONTrack class ???? - slopeBending = param->GetBendingSlope(); - slopeNonBending = param->GetNonBendingSlope(); - // thickness in radiation length for the current track, - // taking local angle into account - radiationLength = - trackBeingFitted->GetEventReconstructor()->GetChamberThicknessInX0() * - TMath::Sqrt(1.0 + - slopeBending * slopeBending + slopeNonBending * slopeNonBending); - inverseBendingMomentum2 = - param->GetInverseBendingMomentum() * param->GetInverseBendingMomentum(); - inverseTotalMomentum2 = - inverseBendingMomentum2 * (1.0 + slopeBending * slopeBending) / - (1.0 + slopeBending *slopeBending + slopeNonBending * slopeNonBending); - varMultipleScatteringAngle = 0.0136 * (1.0 + 0.038 * TMath::Log(radiationLength)); - // The velocity is assumed to be 1 !!!! - varMultipleScatteringAngle = inverseTotalMomentum2 * radiationLength * - varMultipleScatteringAngle * varMultipleScatteringAngle; - return varMultipleScatteringAngle; + /// pack the local trigger information and store + + if (loCirc < 0) return; + + fLocalTrigger = 0; + fLocalTrigger += loCirc; + fLocalTrigger += loStripX << 8; + fLocalTrigger += loStripY << 13; + fLocalTrigger += loDev << 17; + fLocalTrigger += loLpt << 22; + fLocalTrigger += loHpt << 24; + } +