X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=MUON%2FAliMUONTrack.cxx;h=3b98ce745a069242671f70ad4541f0ae372939ec;hb=7f7dd416164e9335d2ed96e06048a8705b6b8d9f;hp=54f370a9fff51732b293f73c780874a136b9fbe2;hpb=30178c30974cdd6a3b59f09e4d479925642e175b;p=u%2Fmrichter%2FAliRoot.git diff --git a/MUON/AliMUONTrack.cxx b/MUON/AliMUONTrack.cxx index 54f370a9fff..3b98ce745a0 100644 --- a/MUON/AliMUONTrack.cxx +++ b/MUON/AliMUONTrack.cxx @@ -15,834 +15,1230 @@ /* $Id$ */ -/////////////////////////////////////////////////// -// -// Reconstructed track -// in -// ALICE -// dimuon -// spectrometer -// -/////////////////////////////////////////////////// - -#include // for exit() - -#include // for cout -#include -#include -#include -#include +//----------------------------------------------------------------------------- +// Class AliMUONTrack +//------------------- +// Reconstructed track in ALICE dimuon spectrometer +//----------------------------------------------------------------------------- -#include "AliMUONEventReconstructor.h" -#include "AliMUONHitForRec.h" -#include "AliMUONSegment.h" #include "AliMUONTrack.h" -#include "AliMUONTrackHit.h" -#include "AliMUONTriggerTrack.h" + +#include "AliMUONReconstructor.h" +#include "AliMUONVCluster.h" +#include "AliMUONVClusterStore.h" +#include "AliMUONObjectPair.h" +#include "AliMUONTrackExtrap.h" #include "AliMUONConstants.h" +#include "AliMUONTrackParam.h" -// 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 "AliLog.h" -void mnvertLocal(Double_t* a, Int_t l, Int_t m, Int_t n, Int_t& ifail); +#include -Double_t MultipleScatteringAngle2(AliMUONTrackHit *TrackHit); +#include +/// \cond CLASSIMP ClassImp(AliMUONTrack) // Class implementation in ROOT context +/// \endcond -TVirtualFitter* AliMUONTrack::fgFitter = NULL; - //__________________________________________________________________________ -AliMUONTrack::AliMUONTrack() - : TObject() -{ - // Default constructor - fgFitter = 0; - fEventReconstructor = 0; - fTrackHitsPtr = 0; - fTrackParamAtHit = new TClonesArray("AliMUONTrackParam",10); -} +const Double_t AliMUONTrack::fgkMaxChi2 = 1.e10; ///< maximum chi2 above which the track can be considered as abnormal - //__________________________________________________________________________ -AliMUONTrack::AliMUONTrack(AliMUONSegment* BegSegment, AliMUONSegment* EndSegment, AliMUONEventReconstructor* EventReconstructor) - : TObject() -{ - // 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 - fTrackParamAtHit = new TClonesArray("AliMUONTrackParam",10); - // set fit conditions... - fFitMCS = 0; - fFitNParam = 3; - fFitStart = 1; - fFitFMin = -1.0; - fMatchTrigger = kFALSE; - fChi2MatchTrigger = 0; - return; -} - //__________________________________________________________________________ -AliMUONTrack::AliMUONTrack(AliMUONSegment* Segment, AliMUONHitForRec* HitForRec, AliMUONEventReconstructor* EventReconstructor) - : TObject() +//__________________________________________________________________________ +AliMUONTrack::AliMUONTrack() + : TObject(), + fTrackParamAtCluster(0x0), + fFitWithVertex(kFALSE), + fVertexErrXY2(), + fFitWithMCS(kFALSE), + fClusterWeightsNonBending(0x0), + fClusterWeightsBending(0x0), + fGlobalChi2(-1.), + fImproved(kFALSE), + fMatchTrigger(-1), + fChi2MatchTrigger(0.), + fTrackID(-1), + fTrackParamAtVertex(0x0), + fHitsPatternInTrigCh(0), + fLocalTrigger(0), + fConnected(kFALSE) { - // 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 - fTrackParamAtHit = new TClonesArray("AliMUONTrackParam",10); - // set fit conditions... - fFitMCS = 0; - fFitNParam = 3; - fFitStart = 1; - fFitFMin = -1.0; - fMatchTrigger = kFALSE; - fChi2MatchTrigger = 0; - return; + /// Default constructor + fVertexErrXY2[0] = 0.; + fVertexErrXY2[1] = 0.; } //__________________________________________________________________________ -AliMUONTrack::~AliMUONTrack() +AliMUONTrack::AliMUONTrack(AliMUONObjectPair *segment, Double_t bendingVertexDispersion) + : TObject(), + fTrackParamAtCluster(new TClonesArray("AliMUONTrackParam",10)), + fFitWithVertex(kFALSE), + fVertexErrXY2(), + fFitWithMCS(kFALSE), + fClusterWeightsNonBending(0x0), + fClusterWeightsBending(0x0), + fGlobalChi2(0.), + fImproved(kFALSE), + fMatchTrigger(-1), + fChi2MatchTrigger(0.), + fTrackID(-1), + fTrackParamAtVertex(0x0), + fHitsPatternInTrigCh(0), + fLocalTrigger(0), + fConnected(kFALSE) { - // Destructor - if (fTrackHitsPtr) { - delete fTrackHitsPtr; // delete the TObjArray of pointers to TrackHit's - fTrackHitsPtr = NULL; - } + /// Constructor from two clusters + + fVertexErrXY2[0] = 0.; + fVertexErrXY2[1] = 0.; + + // Pointers to clusters from the segment + AliMUONVCluster* firstCluster = (AliMUONVCluster*) segment->First(); + AliMUONVCluster* lastCluster = (AliMUONVCluster*) segment->Second(); - if (fTrackParamAtHit) { - // delete the TClonesArray of pointers to TrackParam - delete fTrackParamAtHit; - fTrackParamAtHit = NULL; + // 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 + 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 + 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 + TMatrixD paramCov(5,5); + paramCov.Zero(); + // Non bending plane + paramCov(0,0) = firstCluster->GetErrX2(); + paramCov(0,1) = firstCluster->GetErrX2() / dZ; + paramCov(1,0) = paramCov(0,1); + paramCov(1,1) = ( firstCluster->GetErrX2() + lastCluster->GetErrX2() ) / dZ / dZ; + // Bending plane + paramCov(2,2) = firstCluster->GetErrY2(); + paramCov(2,3) = firstCluster->GetErrY2() / dZ; + paramCov(3,2) = paramCov(2,3); + paramCov(3,3) = ( firstCluster->GetErrY2() + lastCluster->GetErrY2() ) / dZ / dZ; + // Inverse bending momentum (vertex resolution + bending slope resolution + 10% error on dipole parameters+field) + if (AliMUONTrackExtrap::IsFieldON()) { + paramCov(4,4) = ( ( bendingVertexDispersion*bendingVertexDispersion + + (z1 * z1 * lastCluster->GetErrY2() + z2 * z2 * firstCluster->GetErrY2()) / dZ / dZ) / + bendingImpact / bendingImpact + 0.1 * 0.1) * inverseBendingMomentum * inverseBendingMomentum ; + paramCov(2,4) = - z2 * firstCluster->GetErrY2() * inverseBendingMomentum / bendingImpact / dZ; + paramCov(4,2) = paramCov(2,4); + paramCov(3,4) = - (z1 * lastCluster->GetErrY2() + z2 * firstCluster->GetErrY2()) * inverseBendingMomentum / bendingImpact / dZ / dZ; + paramCov(4,3) = paramCov(3,4); + } else paramCov(4,4) = inverseBendingMomentum*inverseBendingMomentum; + trackParamAtFirstCluster.SetCovariances(paramCov); + + // Compute and set track parameters covariances at last cluster + // Non bending plane + paramCov(0,0) = lastCluster->GetErrX2(); + paramCov(0,1) = - lastCluster->GetErrX2() / dZ; + paramCov(1,0) = paramCov(0,1); + // Bending plane + paramCov(2,2) = lastCluster->GetErrY2(); + paramCov(2,3) = - lastCluster->GetErrY2() / dZ; + paramCov(3,2) = paramCov(2,3); + // Inverse bending momentum (vertex resolution + bending slope resolution + 10% error on dipole parameters+field) + if (AliMUONTrackExtrap::IsFieldON()) { + paramCov(2,4) = z1 * lastCluster->GetErrY2() * inverseBendingMomentum / bendingImpact / dZ; + paramCov(4,2) = paramCov(2,4); } + trackParamAtLastCluster.SetCovariances(paramCov); + + // Add track parameters at clusters + AddTrackParamAtCluster(trackParamAtFirstCluster,*firstCluster); + AddTrackParamAtCluster(trackParamAtLastCluster,*lastCluster); + } - //__________________________________________________________________________ -AliMUONTrack::AliMUONTrack (const AliMUONTrack& theMUONTrack) - : TObject(theMUONTrack) +//__________________________________________________________________________ +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), + fChi2MatchTrigger(track.fChi2MatchTrigger), + fTrackID(track.fTrackID), + fTrackParamAtVertex(0x0), + fHitsPatternInTrigCh(track.fHitsPatternInTrigCh), + fLocalTrigger(track.fLocalTrigger), + fConnected(track.fConnected) { - //fEventReconstructor = new AliMUONEventReconstructor(*MUONTrack.fEventReconstructor); - // is it right ? - // NO, because it would use dummy copy constructor - // and AliMUONTrack is not the owner of its EventReconstructor - fEventReconstructor = theMUONTrack.fEventReconstructor; - fTrackParamAtVertex = theMUONTrack.fTrackParamAtVertex; - fTrackHitsPtr = new TObjArray(*theMUONTrack.fTrackHitsPtr); // is it right ? - fTrackParamAtHit = new TClonesArray(*theMUONTrack.fTrackParamAtHit); - fNTrackHits = theMUONTrack.fNTrackHits; - fFitMCS = theMUONTrack.fFitMCS; - fFitNParam = theMUONTrack.fFitNParam; - fFitFMin = theMUONTrack.fFitFMin; - fFitStart = theMUONTrack.fFitStart; - fMatchTrigger = theMUONTrack.fMatchTrigger; - fChi2MatchTrigger = theMUONTrack.fChi2MatchTrigger; + ///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::operator=(const AliMUONTrack& theMUONTrack) +AliMUONTrack & AliMUONTrack::operator=(const AliMUONTrack& track) { - + /// Asignment operator // check assignement to self - if (this == &theMUONTrack) + if (this == &track) return *this; // base class assignement - TObject::operator=(theMUONTrack); - - // fEventReconstructor = new AliMUONEventReconstructor(*MUONTrack.fEventReconstructor); // is it right ? - // is it right ? NO because it would use dummy copy constructor - fEventReconstructor = theMUONTrack.fEventReconstructor; - fTrackParamAtVertex = theMUONTrack.fTrackParamAtVertex; - fTrackHitsPtr = new TObjArray(*theMUONTrack.fTrackHitsPtr); // is it right ? - fTrackParamAtHit = new TClonesArray(*theMUONTrack.fTrackParamAtHit); - fNTrackHits = theMUONTrack.fNTrackHits; - fFitMCS = theMUONTrack.fFitMCS; - fFitNParam = theMUONTrack.fFitNParam; - fFitFMin = theMUONTrack.fFitFMin; - fFitStart = theMUONTrack.fFitStart; - fMatchTrigger = theMUONTrack.fMatchTrigger; - fChi2MatchTrigger = theMUONTrack.fChi2MatchTrigger; + 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; + fChi2MatchTrigger = track.fChi2MatchTrigger; + fTrackID = track.fTrackID; + fHitsPatternInTrigCh = track.fHitsPatternInTrigCh; + fLocalTrigger = track.fLocalTrigger; + fConnected = track.fConnected; return *this; } //__________________________________________________________________________ -void AliMUONTrack::Remove() +AliMUONTrack::~AliMUONTrack() { - // 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; - } - // 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(); + /// Destructor + delete fTrackParamAtCluster; + delete fClusterWeightsNonBending; + delete fClusterWeightsBending; + delete fTrackParamAtVertex; } //__________________________________________________________________________ -void AliMUONTrack::SetFitMCS(Int_t FitMCS) +void AliMUONTrack::Clear(Option_t* opt) { - // 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) ???? + /// Clear arrays + if (opt && opt[0] == 'C' && fTrackParamAtCluster) fTrackParamAtCluster->Clear("C"); else { - cout << "ERROR in AliMUONTrack::SetFitMCS(FitMCS)" << endl; - cout << "FitMCS = " << FitMCS << " is neither 0 nor 1" << endl; - exit(0); + delete fTrackParamAtCluster; + fTrackParamAtCluster = 0x0; } - return; + delete fClusterWeightsNonBending; fClusterWeightsNonBending = 0x0; + delete fClusterWeightsBending; fClusterWeightsBending = 0x0; + delete fTrackParamAtVertex; fTrackParamAtVertex = 0x0; } //__________________________________________________________________________ -void AliMUONTrack::SetFitNParam(Int_t FitNParam) +void AliMUONTrack::Reset() { - // 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); - } - return; + /// Reset to default values + SetUniqueID(0); + fFitWithVertex = kFALSE; + fVertexErrXY2[0] = 0.; + fVertexErrXY2[1] = 0.; + fFitWithMCS = kFALSE; + fGlobalChi2 = -1.; + fImproved = kFALSE; + fMatchTrigger = -1; + fChi2MatchTrigger = 0.; + fTrackID = -1; + fHitsPatternInTrigCh = 0; + fLocalTrigger = 0; + fConnected = kFALSE; + delete fTrackParamAtCluster; fTrackParamAtCluster = 0x0; + delete fClusterWeightsNonBending; fClusterWeightsNonBending = 0x0; + delete fClusterWeightsBending; fClusterWeightsBending = 0x0; + delete fTrackParamAtVertex; fTrackParamAtVertex = 0x0; } //__________________________________________________________________________ -void AliMUONTrack::SetFitStart(Int_t FitStart) +TClonesArray* AliMUONTrack::GetTrackParamAtCluster() const { - // 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); - } - return; + /// return array of track parameters at cluster (create it if needed) + if (!fTrackParamAtCluster) fTrackParamAtCluster = new TClonesArray("AliMUONTrackParam",10); + return fTrackParamAtCluster; } //__________________________________________________________________________ -AliMUONTrackParam* AliMUONTrack::GetTrackParamAtFirstHit(void) const { - // Get pointer to TrackParamAtFirstHit - return ((AliMUONTrackHit*) (fTrackHitsPtr->First()))->GetTrackParam();} +void AliMUONTrack::AddTrackParamAtCluster(const AliMUONTrackParam &trackParam, AliMUONVCluster &cluster, Bool_t copy) +{ + /// 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; + } + + // 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::RecursiveDump(void) const +void AliMUONTrack::RemoveTrackParamAtCluster(AliMUONTrackParam *trackParam) { - // 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(); + /// 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(); } //__________________________________________________________________________ -Int_t AliMUONTrack::HitsInCommon(AliMUONTrack* Track) const +Bool_t AliMUONTrack::UpdateTrackParamAtCluster() { - // 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; + /// Update track parameters at each attached cluster + /// Return kFALSE in case of failure (i.e. extrapolation problem) + + if (GetNClusters() == 0) { + AliWarning("no cluster attached to the track"); + return kFALSE; + } + + Bool_t extrapStatus = kTRUE; + 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 + if (!AliMUONTrackExtrap::ExtrapToZ(trackParamAtCluster, trackParamAtCluster->GetClusterPtr()->GetZ())) extrapStatus = kFALSE; + + // prepare next step + startingTrackParam = trackParamAtCluster; + trackParamAtCluster = (AliMUONTrackParam*) (fTrackParamAtCluster->After(trackParamAtCluster)); + } + + // set global chi2 to max value in case of problem during track extrapolation + if (!extrapStatus) SetGlobalChi2(2.*MaxChi2()); + return extrapStatus; + } //__________________________________________________________________________ -void AliMUONTrack::MatchTriggerTrack(TClonesArray *triggerTrackArray) +Bool_t AliMUONTrack::UpdateCovTrackParamAtCluster() { - // Match this track with one trigger track if possible - AliMUONTrackParam *trackParam; - AliMUONTriggerTrack *triggerTrack; - Double_t xTrack, yTrack, ySlopeTrack, dTrigTrackMin2, dTrigTrack2; - Double_t nSigmaCut2; - - Double_t distSigma[3]={1,1,0.02}; // sigma of distributions (trigger-track) X,Y,slopeY - Double_t distTriggerTrack[3] = {0,0,0}; - - fMatchTrigger = kFALSE; - fChi2MatchTrigger = 0; - - trackParam = (AliMUONTrackParam*) fTrackParamAtHit->Last(); - trackParam->ExtrapToZ(AliMUONConstants::DefaultChamberZ(10)); // extrap to 1st trigger chamber - - nSigmaCut2 = fEventReconstructor->GetMaxSigma2Distance(); // nb of sigma**2 for cut - xTrack = trackParam->GetNonBendingCoor(); - yTrack = trackParam->GetBendingCoor(); - ySlopeTrack = trackParam->GetBendingSlope(); - dTrigTrackMin2 = 999; - - triggerTrack = (AliMUONTriggerTrack*) triggerTrackArray->First(); - while(triggerTrack){ - distTriggerTrack[0] = (triggerTrack->GetX11()-xTrack)/distSigma[0]; - distTriggerTrack[1] = (triggerTrack->GetY11()-yTrack)/distSigma[1]; - distTriggerTrack[2] = (TMath::Tan(triggerTrack->GetThetay())-ySlopeTrack)/distSigma[2]; - dTrigTrack2 = 0; - for (Int_t iVar = 0; iVar < 3; iVar++) - dTrigTrack2 += distTriggerTrack[iVar]*distTriggerTrack[iVar]; - if (dTrigTrack2 < dTrigTrackMin2 && dTrigTrack2 < nSigmaCut2) { - dTrigTrackMin2 = dTrigTrack2; - fMatchTrigger = kTRUE; - fChi2MatchTrigger = dTrigTrack2/3.; // Normalized Chi2, 3 variables (X,Y,slopeY) + /// Update track parameters and their covariances at each attached cluster + /// Include effects of multiple scattering in chambers + /// Return kFALSE in case of failure (i.e. extrapolation problem) + + if (GetNClusters() == 0) { + AliWarning("no cluster attached to the track"); + return kFALSE; + } + + Bool_t extrapStatus = kTRUE; + 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(expectedChamber-1),-1.); + + // add MCS in missing chambers if any + currentChamber = trackParamAtCluster->GetClusterPtr()->GetChamberId(); + while (currentChamber > expectedChamber) { + // extrapolation to the missing chamber + if (!AliMUONTrackExtrap::ExtrapToZCov(trackParamAtCluster, AliMUONConstants::DefaultChamberZ(expectedChamber))) extrapStatus = kFALSE; + // add MCS effect + AliMUONTrackExtrap::AddMCSEffect(trackParamAtCluster,AliMUONConstants::ChamberThicknessInX0(expectedChamber),-1.); + expectedChamber++; } - triggerTrack = (AliMUONTriggerTrack*) triggerTrackArray->After(triggerTrack); + + // extrapolation to the z of the current cluster + if (!AliMUONTrackExtrap::ExtrapToZCov(trackParamAtCluster, trackParamAtCluster->GetClusterPtr()->GetZ())) extrapStatus = kFALSE; + + // prepare next step + expectedChamber = currentChamber + 1; + startingTrackParam = trackParamAtCluster; + trackParamAtCluster = (AliMUONTrackParam*) (fTrackParamAtCluster->After(trackParamAtCluster)); } - + + // set global chi2 to max value in case of problem during track extrapolation + if (!extrapStatus) SetGlobalChi2(2.*MaxChi2()); + return extrapStatus; + } + //__________________________________________________________________________ -void AliMUONTrack::Fit() +Bool_t AliMUONTrack::IsValid(UInt_t requestedStationMask, Bool_t request2ChInSameSt45) { - // 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]; + /// check the validity of the current track: + /// at least one cluster per requested station + /// and at least 2 chambers in stations 4 & 5 that contain cluster(s) + /// + if request2ChInSameSt45 = kTRUE: 2 chambers hit in the same station (4 or 5) + + Int_t nClusters = GetNClusters(); 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); - // Switch off printout - arg[0] = -1; - fgFitter->ExecuteCommand("SET PRINT", arg, 1); // More printing !!!! - // No warnings - fgFitter->ExecuteCommand("SET NOW", arg, 0); - // 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 - // mandatory limits in Bending to avoid NaN values of parameters - fgFitter->SetParameter(3, "X", - trackParam->GetNonBendingCoor(), - 0.03, -500.0, 500.0); - // mandatory limits in non Bending to avoid NaN values of parameters - fgFitter->SetParameter(4, "Y", - trackParam->GetBendingCoor(), - 0.10, -500.0, 500.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); - } - // 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); - } - // global result of the fit - Double_t fedm, errdef; - Int_t npari, nparx; - fgFitter->GetStats(fFitFMin, fedm, errdef, npari, nparx); + Int_t currentCh, currentSt, previousCh = -1, nChHitInSt4 = 0, nChHitInSt5 = 0; + UInt_t presentStationMask(0); + + // first loop over clusters + for (Int_t i = 0; i < nClusters; i++) { + trackParam = (AliMUONTrackParam*) fTrackParamAtCluster->UncheckedAt(i); + + currentCh = trackParam->GetClusterPtr()->GetChamberId(); + currentSt = currentCh/2; + + // build present station mask + presentStationMask |= ( 1 << currentSt ); + + // count the number of chambers hit in station 4 that contain cluster(s) + if (currentSt == 3 && currentCh != previousCh) { + nChHitInSt4++; + previousCh = currentCh; + } + + // count the number of chambers hit in station 5 that contain cluster(s) + if (currentSt == 4 && currentCh != previousCh) { + nChHitInSt5++; + previousCh = currentCh; + } + + } + + // at least one cluster per requested station + if ((requestedStationMask & presentStationMask) != requestedStationMask) return kFALSE; + + // 2 chambers hit in the same station (4 or 5) + if (request2ChInSameSt45) return (nChHitInSt4 == 2 || nChHitInSt5 == 2); + // or 2 chambers hit in station 4 & 5 together + else return (nChHitInSt4+nChHitInSt5 >= 2); + } //__________________________________________________________________________ -void AliMUONTrack::AddSegment(AliMUONSegment* Segment) -{ - // Add Segment to the track - AddHitForRec(Segment->GetHitForRec1()); // 1st hit - AddHitForRec(Segment->GetHitForRec2()); // 2nd hit +void AliMUONTrack::TagRemovableClusters(UInt_t requestedStationMask) { + /// Identify clusters that can be removed from the track, + /// with the only requirements to have at least 1 cluster per requested station + /// and at least 2 chambers over 4 in stations 4 & 5 that contain cluster(s) + + Int_t nClusters = GetNClusters(); + AliMUONTrackParam *trackParam, *nextTrackParam; + Int_t currentCh, nextCh, currentSt, nextSt, previousCh = -1, nChHitInSt45 = 0; + + // first loop over clusters + for (Int_t i = 0; i < nClusters; i++) { + trackParam = (AliMUONTrackParam*) fTrackParamAtCluster->UncheckedAt(i); + + currentCh = trackParam->GetClusterPtr()->GetChamberId(); + currentSt = currentCh/2; + + // reset flags to kFALSE for all clusters in required station + if ((1 << currentSt) & requestedStationMask) trackParam->SetRemovable(kFALSE); + else trackParam->SetRemovable(kTRUE); + + // count the number of chambers in station 4 & 5 that contain cluster(s) + if (currentCh > 5 && currentCh != previousCh) { + nChHitInSt45++; + previousCh = currentCh; + } + + } + + // second loop over clusters + for (Int_t i = 0; i < nClusters; i++) { + trackParam = (AliMUONTrackParam*) fTrackParamAtCluster->UncheckedAt(i); + + currentCh = trackParam->GetClusterPtr()->GetChamberId(); + currentSt = currentCh/2; + + // make sure they are more than 2 clusters in 2 different chambers of stations 4 & 5 + // but 2 clusters in he same chamber will still be flagged as removable + if (nChHitInSt45 < 3 && currentSt > 2) { + + if (i == nClusters-1) { + + trackParam->SetRemovable(kFALSE); + + } else { + + nextTrackParam = (AliMUONTrackParam*) fTrackParamAtCluster->UncheckedAt(i+1); + nextCh = nextTrackParam->GetClusterPtr()->GetChamberId(); + + // set clusters in the same chamber as being removable + if (nextCh == currentCh) { + trackParam->SetRemovable(kTRUE); + nextTrackParam->SetRemovable(kTRUE); + i++; // skip cluster already checked + } else { + trackParam->SetRemovable(kFALSE); + } + + } + + } else { + + // skip clusters already flag as removable + if (trackParam->IsRemovable()) continue; + + // loop over next track parameters + for (Int_t j = i+1; j < nClusters; j++) { + nextTrackParam = (AliMUONTrackParam*) fTrackParamAtCluster->UncheckedAt(j); + + nextCh = nextTrackParam->GetClusterPtr()->GetChamberId(); + nextSt = nextCh/2; + + // set clusters in the same station as being removable + if (nextSt == currentSt) { + trackParam->SetRemovable(kTRUE); + nextTrackParam->SetRemovable(kTRUE); + i++; // skip cluster already checked + } + + } + + } + + } + } //__________________________________________________________________________ -void AliMUONTrack::AddHitForRec(AliMUONHitForRec* HitForRec) +Bool_t AliMUONTrack::ComputeLocalChi2(Bool_t accountForMCS) { - // 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 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::SetTrackParamAtHit(Int_t indexHit, AliMUONTrackParam *TrackParam) const +Double_t AliMUONTrack::ComputeGlobalChi2(Bool_t accountForMCS) { - // Set track parameters at TrackHit with index "indexHit" - // from the track parameters pointed to by "TrackParam". - //PH AliMUONTrackHit* trackHit = (AliMUONTrackHit*) ((*fTrackHitsPtr)[indexHit]); - AliMUONTrackHit* trackHit = (AliMUONTrackHit*) (fTrackHitsPtr->At(indexHit)); - trackHit->SetTrackParam(TrackParam); + /// 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 a value of chi2 higher than the maximum allowed if computation failed + AliDebug(1,"Enter ComputeGlobalChi2"); + + if (!fTrackParamAtCluster) { + AliWarning("no cluster attached to this track"); + return 2.*MaxChi2(); + } + + 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(); + } + + } + + return chi2; + } //__________________________________________________________________________ -void AliMUONTrack::SetTrackParamAtVertex() +Bool_t AliMUONTrack::ComputeClusterWeights(TMatrixD* mcsCovariances) { - // 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); + /// 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 TrackChi2(Int_t &NParam, Double_t * /*Gradient*/, Double_t &Chi2, Double_t *Param, Int_t /*Flag*/) +Bool_t AliMUONTrack::ComputeClusterWeights(TMatrixD& clusterWeightsNB, TMatrixD& clusterWeightsB, + TMatrixD* mcsCovariances, AliMUONVCluster* discardedCluster) 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(); + /// 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 TrackChi2MCS(Int_t &NParam, Double_t * /*Gradient*/, Double_t &Chi2, Double_t *Param, Int_t /*Flag*/) +void AliMUONTrack::ComputeMCSCovariances(TMatrixD& mcsCovariances) 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]); - } - - AliMUONTrackHit *hit; - Int_t chCurrent, chPrev = 0, hitNumber, hitNumber1, hitNumber2, hitNumber3; - Double_t z, z1, z2, z3; - AliMUONTrackHit *hit1, *hit2, *hit3; - Double_t hbc1, hbc2, pbc1, pbc2; - Double_t hnbc1, hnbc2, pnbc1, pnbc2; - Int_t numberOfHit = trackBeingFitted->GetNTrackHits(); - TMatrixD *covBending = new TMatrixD(numberOfHit, numberOfHit); - TMatrixD *covNonBending = new TMatrixD(numberOfHit, numberOfHit); - Double_t *msa2 = new Double_t[numberOfHit]; - - // Predicted coordinates and multiple scattering angles are first calculated - for (hitNumber = 0; hitNumber < numberOfHit; hitNumber++) { - hit = (AliMUONTrackHit*) (*(trackBeingFitted->GetTrackHitsPtr()))[hitNumber]; - z = hit->GetHitForRecPtr()->GetZ(); - // do something special if 2 hits with same Z ???? - // security against infinite loop ???? - (¶m1)->ExtrapToZ(z); // extrapolation - hit->SetTrackParam(¶m1); - // square of multiple scattering angle at current hit, with one chamber - msa2[hitNumber] = MultipleScatteringAngle2(hit); - // correction for eventual missing hits or multiple hits in a chamber, - // according to the number of chambers - // between the current hit and the previous one - chCurrent = hit->GetHitForRecPtr()->GetChamberNumber(); - if (hitNumber > 0) msa2[hitNumber] = msa2[hitNumber] * (chCurrent - chPrev); - chPrev = chCurrent; - } - - // Calculates the covariance matrix - for (hitNumber1 = 0; hitNumber1 < numberOfHit; hitNumber1++) { - hit1 = (AliMUONTrackHit*) (*(trackBeingFitted->GetTrackHitsPtr()))[hitNumber1]; - z1 = hit1->GetHitForRecPtr()->GetZ(); - for (hitNumber2 = hitNumber1; hitNumber2 < numberOfHit; hitNumber2++) { - hit2 = (AliMUONTrackHit*) (*(trackBeingFitted->GetTrackHitsPtr()))[hitNumber2]; - z2 = hit2->GetHitForRecPtr()->GetZ(); - // initialization to 0 (diagonal plus upper triangular part) - (*covBending)(hitNumber2, hitNumber1) = 0.0; - // contribution from multiple scattering in bending plane: - // loop over upstream hits - for (hitNumber3 = 0; hitNumber3 < hitNumber1; hitNumber3++) { - hit3 = (AliMUONTrackHit*) - (*(trackBeingFitted->GetTrackHitsPtr()))[hitNumber3]; - z3 = hit3->GetHitForRecPtr()->GetZ(); - (*covBending)(hitNumber2, hitNumber1) = - (*covBending)(hitNumber2, hitNumber1) + - ((z1 - z3) * (z2 - z3) * msa2[hitNumber3]); - } - // equal contribution from multiple scattering in non bending plane - (*covNonBending)(hitNumber2, hitNumber1) = - (*covBending)(hitNumber2, hitNumber1); - if (hitNumber1 == hitNumber2) { - // Diagonal elements: add contribution from position measurements - // in bending plane - (*covBending)(hitNumber2, hitNumber1) = - (*covBending)(hitNumber2, hitNumber1) + - hit1->GetHitForRecPtr()->GetBendingReso2(); - // and in non bending plane - (*covNonBending)(hitNumber2, hitNumber1) = - (*covNonBending)(hitNumber2, hitNumber1) + - hit1->GetHitForRecPtr()->GetNonBendingReso2(); - } - else { - // Non diagonal elements: symmetrization - // for bending plane - (*covBending)(hitNumber1, hitNumber2) = - (*covBending)(hitNumber2, hitNumber1); - // and non bending plane - (*covNonBending)(hitNumber1, hitNumber2) = - (*covNonBending)(hitNumber2, hitNumber1); - } - } // for (hitNumber2 = hitNumber1;... - } // for (hitNumber1 = 0;... - - // Inversion of covariance matrices - // with "mnvertLocal", local "mnvert" function of Minuit. - // One cannot use directly "mnvert" since "TVirtualFitter" does not know it. - // One will have to replace this local function by the right inversion function - // from a specialized Root package for symmetric positive definite matrices, - // when available!!!! - Int_t ifailBending; - mnvertLocal(&((*covBending)(0,0)), numberOfHit, numberOfHit, numberOfHit, - ifailBending); - Int_t ifailNonBending; - mnvertLocal(&((*covNonBending)(0,0)), numberOfHit, numberOfHit, numberOfHit, - ifailNonBending); - - // It would be worth trying to calculate the inverse of the covariance matrix - // only once per fit, since it cannot change much in principle, - // and it would save a lot of computing time !!!! - - // Calculates Chi2 - if ((ifailBending == 0) && (ifailNonBending == 0)) { - // with Multiple Scattering if inversion correct - for (hitNumber1 = 0; hitNumber1 < numberOfHit ; hitNumber1++) { - hit1 = (AliMUONTrackHit*) (*(trackBeingFitted->GetTrackHitsPtr()))[hitNumber1]; - hbc1 = hit1->GetHitForRecPtr()->GetBendingCoor(); - pbc1 = hit1->GetTrackParam()->GetBendingCoor(); - hnbc1 = hit1->GetHitForRecPtr()->GetNonBendingCoor(); - pnbc1 = hit1->GetTrackParam()->GetNonBendingCoor(); - for (hitNumber2 = 0; hitNumber2 < numberOfHit; hitNumber2++) { - hit2 = (AliMUONTrackHit*) - (*(trackBeingFitted->GetTrackHitsPtr()))[hitNumber2]; - hbc2 = hit2->GetHitForRecPtr()->GetBendingCoor(); - pbc2 = hit2->GetTrackParam()->GetBendingCoor(); - hnbc2 = hit2->GetHitForRecPtr()->GetNonBendingCoor(); - pnbc2 = hit2->GetTrackParam()->GetNonBendingCoor(); - Chi2 = Chi2 + - ((*covBending)(hitNumber2, hitNumber1) * - (hbc1 - pbc1) * (hbc2 - pbc2)) + - ((*covNonBending)(hitNumber2, hitNumber1) * - (hnbc1 - pnbc1) * (hnbc2 - pnbc2)); + /// 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(expectedChamber),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(currentChamber),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); } - } else { - // without Multiple Scattering if inversion impossible - for (hitNumber1 = 0; hitNumber1 < numberOfHit ; hitNumber1++) { - hit1 = (AliMUONTrackHit*) (*(trackBeingFitted->GetTrackHitsPtr()))[hitNumber1]; - hbc1 = hit1->GetHitForRecPtr()->GetBendingCoor(); - pbc1 = hit1->GetTrackParam()->GetBendingCoor(); - hnbc1 = hit1->GetHitForRecPtr()->GetNonBendingCoor(); - pnbc1 = hit1->GetTrackParam()->GetNonBendingCoor(); - Chi2 = Chi2 + - ((hbc1 - pbc1) * (hbc1 - pbc1) / - hit1->GetHitForRecPtr()->GetBendingReso2()) + - ((hnbc1 - pnbc1) * (hnbc1 - pnbc1) / - hit1->GetHitForRecPtr()->GetNonBendingReso2()); + + } + + delete [] mcsAngle2; + delete [] zMCS; + delete [] indices; + +} + + //__________________________________________________________________________ +Int_t AliMUONTrack::ClustersInCommon(AliMUONTrack* track, Int_t stMin, Int_t stMax) const +{ + /// Returns the number of clusters in common in stations [stMin, stMax] + /// between the current track ("this") and the track pointed to by "track". + + if (!fTrackParamAtCluster || !this->fTrackParamAtCluster) return 0; + + Int_t chMin = 2 * stMin; + Int_t chMax = 2 * stMax + 1; + Int_t clustersInCommon = 0; + + // Loop over clusters of first track + Int_t nCl1 = this->GetNClusters(); + for(Int_t iCl1 = 0; iCl1 < nCl1; iCl1++) { + AliMUONVCluster* cl1 = ((AliMUONTrackParam*) this->fTrackParamAtCluster->UncheckedAt(iCl1))->GetClusterPtr(); + + Int_t chCl1 = cl1->GetChamberId(); + if (chCl1 < chMin || chCl1 > chMax) continue; + + // Loop over clusters of second track + Int_t nCl2 = track->GetNClusters(); + for(Int_t iCl2 = 0; iCl2 < nCl2; iCl2++) { + AliMUONVCluster* cl2 = ((AliMUONTrackParam*) track->fTrackParamAtCluster->UncheckedAt(iCl2))->GetClusterPtr(); + + Int_t chCl2 = cl2->GetChamberId(); + if (chCl2 < chMin || chCl2 > chMax) continue; + + // Increment "clustersInCommon" if both clusters have the same ID + if (cl1->GetUniqueID() == cl2->GetUniqueID()) { + clustersInCommon++; + break; + } + } + } - delete covBending; - delete covNonBending; - delete [] msa2; + return clustersInCommon; } -Double_t MultipleScatteringAngle2(AliMUONTrackHit *TrackHit) + //__________________________________________________________________________ +Int_t AliMUONTrack::GetNDF() const { - // 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; + /// return the number of degrees of freedom + + Int_t ndf = 2 * GetNClusters() - 5; + return (ndf > 0) ? ndf : 0; } -//______________________________________________________________________________ - void mnvertLocal(Double_t *a, Int_t l, Int_t, Int_t n, Int_t &ifail) + //__________________________________________________________________________ +Double_t AliMUONTrack::GetNormalizedChi2() const { -//*-*-*-*-*-*-*-*-*-*-*-*Inverts a symmetric matrix*-*-*-*-*-*-*-*-*-*-*-*-* -//*-* ========================== -//*-* inverts a symmetric matrix. matrix is first scaled to -//*-* have all ones on the diagonal (equivalent to change of units) -//*-* but no pivoting is done since matrix is positive-definite. -//*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-* - - // taken from TMinuit package of Root (l>=n) - // fVERTs, fVERTq and fVERTpp changed to localVERTs, localVERTq and localVERTpp - // Double_t localVERTs[n], localVERTq[n], localVERTpp[n]; - Double_t * localVERTs = new Double_t[n]; - Double_t * localVERTq = new Double_t[n]; - Double_t * localVERTpp = new Double_t[n]; - // fMaxint changed to localMaxint - Int_t localMaxint = n; - - /* System generated locals */ - Int_t aOffset; - - /* Local variables */ - Double_t si; - Int_t i, j, k, kp1, km1; - - /* Parameter adjustments */ - aOffset = l + 1; - a -= aOffset; - - /* Function Body */ - ifail = 0; - if (n < 1) goto L100; - if (n > localMaxint) goto L100; -//*-*- scale matrix by sqrt of diag elements - for (i = 1; i <= n; ++i) { - si = a[i + i*l]; - if (si <= 0) goto L100; - localVERTs[i-1] = 1 / TMath::Sqrt(si); - } - for (i = 1; i <= n; ++i) { - for (j = 1; j <= n; ++j) { - a[i + j*l] = a[i + j*l]*localVERTs[i-1]*localVERTs[j-1]; - } - } -//*-*- . . . start main loop . . . . - for (i = 1; i <= n; ++i) { - k = i; -//*-*- preparation for elimination step1 - if (a[k + k*l] != 0) localVERTq[k-1] = 1 / a[k + k*l]; - else goto L100; - localVERTpp[k-1] = 1; - a[k + k*l] = 0; - kp1 = k + 1; - km1 = k - 1; - if (km1 < 0) goto L100; - else if (km1 == 0) goto L50; - else goto L40; -L40: - for (j = 1; j <= km1; ++j) { - localVERTpp[j-1] = a[j + k*l]; - localVERTq[j-1] = a[j + k*l]*localVERTq[k-1]; - a[j + k*l] = 0; - } -L50: - if (k - n < 0) goto L51; - else if (k - n == 0) goto L60; - else goto L100; -L51: - for (j = kp1; j <= n; ++j) { - localVERTpp[j-1] = a[k + j*l]; - localVERTq[j-1] = -a[k + j*l]*localVERTq[k-1]; - a[k + j*l] = 0; - } -//*-*- elimination proper -L60: - for (j = 1; j <= n; ++j) { - for (k = j; k <= n; ++k) { a[j + k*l] += localVERTpp[j-1]*localVERTq[k-1]; } - } + /// return the chi2 value divided by the number of degrees of freedom (or FLT_MAX if ndf <= 0) + + Double_t ndf = (Double_t) GetNDF(); + return (ndf > 0.) ? fGlobalChi2 / ndf : 2.*MaxChi2(); +} + + //__________________________________________________________________________ +Int_t AliMUONTrack::FindCompatibleClusters(AliMUONTrack &track, Double_t sigmaCut, Bool_t compatibleCluster[10]) const +{ + /// Try to match clusters from this track with clusters from the given track within the provided sigma cut: + /// - Fill the array compatibleCluster[iCh] with kTRUE if a compatible cluster has been found in chamber iCh. + /// - Return the number of clusters of "this" track matched with one cluster of the given track. + AliMUONTrackParam *trackParamAtCluster1, *trackParamAtCluster2; + AliMUONVCluster *cluster1, *cluster2; + Double_t chi2, dX, dY; + Double_t chi2Max = sigmaCut * sigmaCut; + + // initialization + Int_t nMatchClusters = 0; + for ( Int_t ch = 0; ch < AliMUONConstants::NTrackingCh(); ch++) compatibleCluster[ch] = kFALSE; + + if (!track.fTrackParamAtCluster || !this->fTrackParamAtCluster) return nMatchClusters; + + // 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); + + // check DE Id + if (cluster1->GetDetElemId() != cluster2->GetDetElemId()) continue; + + // check local chi2 + dX = cluster1->GetX() - cluster2->GetX(); + dY = cluster1->GetY() - cluster2->GetY(); + chi2 = dX * dX / (cluster1->GetErrX2() + cluster2->GetErrX2()) + dY * dY / (cluster1->GetErrY2() + cluster2->GetErrY2()); + if (chi2 > 2. * chi2Max) continue; // 2 because 2 quantities in chi2 + + compatibleCluster[cluster1->GetChamberId()] = kTRUE; + nMatchClusters++; + break; } -//*-*- elements of left diagonal and unscaling - for (j = 1; j <= n; ++j) { - for (k = 1; k <= j; ++k) { - a[k + j*l] = a[k + j*l]*localVERTs[k-1]*localVERTs[j-1]; - a[j + k*l] = a[k + j*l]; - } + + trackParamAtCluster1 = (AliMUONTrackParam*) this->fTrackParamAtCluster->After(trackParamAtCluster1); + } + + return nMatchClusters; +} + +//__________________________________________________________________________ +Bool_t AliMUONTrack::Match(AliMUONTrack &track, Double_t sigmaCut, Int_t &nMatchClusters) const +{ + /// Try to match this track with the given track. Matching conditions: + /// - more than 50% of clusters from this track matched with clusters from the given track + /// - at least 1 cluster matched before and 1 cluster matched after the dipole + + Bool_t compTrack[10]; + nMatchClusters = FindCompatibleClusters(track, sigmaCut, compTrack); + + if ((compTrack[0] || compTrack[1] || compTrack[2] || compTrack[3]) && // at least 1 cluster matched in st 1 & 2 + (compTrack[6] || compTrack[7] || compTrack[8] || compTrack[9]) && // at least 1 cluster matched in st 4 & 5 + 2 * nMatchClusters > GetNClusters()) return kTRUE; // more than 50% of clusters matched + else return kFALSE; + +} + +//__________________________________________________________________________ +void AliMUONTrack::SetTrackParamAtVertex(const AliMUONTrackParam* trackParam) +{ + /// set track parameters at vertex + if (trackParam == 0x0) return; + if (fTrackParamAtVertex) *fTrackParamAtVertex = *trackParam; + else fTrackParamAtVertex = new AliMUONTrackParam(*trackParam); +} + +//__________________________________________________________________________ +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; +} + +//_____________________________________________- +void AliMUONTrack::Print(Option_t*) const +{ + /// Printing Track information + + cout << " No.Clusters=" << setw(2) << GetNClusters() << + ", Match2Trig=" << setw(1) << GetMatchTrigger() << + ", LoTrgNum=" << setw(3) << LoCircuit() << + ", Chi2-tracking-trigger=" << setw(8) << setprecision(5) << GetChi2MatchTrigger(); + cout << Form(" HitTriggerPattern %x",fHitsPatternInTrigCh); + cout << Form(" MClabel=%d",fTrackID) << endl; + if (fTrackParamAtCluster) fTrackParamAtCluster->First()->Print("FULL"); +} + +//__________________________________________________________________________ +void AliMUONTrack::SetLocalTrigger(Int_t loCirc, Int_t loStripX, Int_t loStripY, Int_t loDev, Int_t loLpt, Int_t loHpt, UChar_t respWithoutChamber) +{ + /// 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; + fLocalTrigger += respWithoutChamber << 26; + +} + +//__________________________________________________________________________ +void AliMUONTrack::FindMCLabel() +{ + /// Determine the MC label from the label of the attached clusters and fill fMCLabel data member: + /// More than 50% of clusters, including 1 before and 1 after the dipole, must share the same label + + Int_t nClusters = GetNClusters(); + Int_t halfCluster = nClusters/2; + + // reset MC label + fTrackID = -1; + + // loop over first clusters (if nClusters left < (nClusters-halfCluster) the conditions cannot be fulfilled) + for (Int_t iCluster1 = 0; iCluster1 < nClusters-halfCluster; iCluster1++) { + AliMUONVCluster* cluster1 = ((AliMUONTrackParam*) fTrackParamAtCluster->UncheckedAt(iCluster1))->GetClusterPtr(); + + // if the first cluster is not on station 1 or 2 the conditions cannot be fulfilled + if (cluster1->GetChamberId() > 3) return; + + Int_t label1 = cluster1->GetMCLabel(); + if (label1 < 0) continue; + + Int_t nIdenticalLabel = 1; + + // Loop over next clusters + for (Int_t iCluster2 = iCluster1+1; iCluster2 < nClusters; iCluster2++) { + AliMUONVCluster* cluster2 = ((AliMUONTrackParam*) fTrackParamAtCluster->UncheckedAt(iCluster2))->GetClusterPtr(); + + if (cluster2->GetMCLabel() != label1) continue; + + nIdenticalLabel++; + + // stop as soon as conditions are fulfilled + if (nIdenticalLabel > halfCluster && cluster2->GetChamberId() > 5) { + fTrackID = label1; + return; + } + } - delete localVERTs; - delete localVERTq; - delete localVERTpp; - return; -//*-*- failure return -L100: - delete localVERTs; - delete localVERTq; - delete localVERTpp; - ifail = 1; -} /* mnvertLocal */ + + } + +}