X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=MUON%2FAliMUONTrack.cxx;h=b3b729e1bfca974962bca79afa6edb2a656eb65e;hb=e6168affcdb6417ecf865b47dfd283875500a284;hp=998e651fdd142adacbb3f4d19bad239b0af940b1;hpb=2c3dd8cd8cad5356b853d09d9515108c30444e97;p=u%2Fmrichter%2FAliRoot.git diff --git a/MUON/AliMUONTrack.cxx b/MUON/AliMUONTrack.cxx index 998e651fdd1..b3b729e1bfc 100644 --- a/MUON/AliMUONTrack.cxx +++ b/MUON/AliMUONTrack.cxx @@ -15,834 +15,1054 @@ /* $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 "AliMUONTrackParam.h" +#include "AliMUONRawClusterV2.h" +#include "AliMUONObjectPair.h" #include "AliMUONConstants.h" +#include "AliMUONTrackExtrap.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" +#include "AliESDMuonTrack.h" +#include "AliESDMuonCluster.h" -void mnvertLocal(Double_t* a, Int_t l, Int_t m, Int_t n, Int_t& ifail); +#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() + : TObject(), + fTrackParamAtCluster(new TClonesArray("AliMUONTrackParam",10)), + 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 - fgFitter = 0; - fEventReconstructor = 0; - fTrackHitsPtr = 0; - fTrackParamAtHit = new TClonesArray("AliMUONTrackParam",10); + /// Default constructor + fTrackParamAtCluster->SetOwner(kTRUE); + fVertexErrXY2[0] = 0.; + fVertexErrXY2[1] = 0.; } //__________________________________________________________________________ -AliMUONTrack::AliMUONTrack(AliMUONSegment* BegSegment, AliMUONSegment* EndSegment, AliMUONEventReconstructor* EventReconstructor) - : TObject() +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 - 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() -{ - // 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; + /// Constructor from two clusters + fTrackParamAtCluster->SetOwner(kTRUE); + + 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 + //(the cluster first attached is the one from which we will start the tracking procedure) + AliMUONVCluster *firstCluster, *lastCluster; + if (cluster1->GetChamberId() == 8) { + 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 + 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 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) + 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() +//__________________________________________________________________________ +AliMUONTrack::AliMUONTrack(AliESDMuonTrack &esdTrack) + : TObject(), + fTrackParamAtCluster(new TClonesArray("AliMUONTrackParam",10)), + fFitWithVertex(kFALSE), + fVertexErrXY2(), + fFitWithMCS(kFALSE), + fClusterWeightsNonBending(0x0), + fClusterWeightsBending(0x0), + fGlobalChi2(esdTrack.GetChi2()), + fImproved(kFALSE), + fMatchTrigger(esdTrack.GetMatchTrigger()), + floTrgNum(-1), + fChi2MatchTrigger(esdTrack.GetChi2MatchTrigger()), + fTrackID(0), + fTrackParamAtVertex(new AliMUONTrackParam()), + fHitsPatternInTrigCh(esdTrack.GetHitsPatternInTrigCh()), + fLocalTrigger(0) { - // Destructor - if (fTrackHitsPtr) { - delete fTrackHitsPtr; // delete the TObjArray of pointers to TrackHit's - fTrackHitsPtr = NULL; + /// Constructor from ESD muon track + /// Compute track parameters and covariances at each cluster if available + /// or store parameters and covariances at first (fake) cluster only if not + + fTrackParamAtCluster->SetOwner(kTRUE); + + fVertexErrXY2[0] = 0.; + fVertexErrXY2[1] = 0.; + + // global info + SetLocalTrigger(esdTrack.LoCircuit(), esdTrack.LoStripX(), esdTrack.LoStripY(), + esdTrack.LoDev(), esdTrack.LoLpt(), esdTrack.LoHpt()); + + // track parameters at vertex + fTrackParamAtVertex->GetParamFrom(esdTrack); + + // track parameters at first cluster + AliMUONTrackParam param; + param.GetParamFromUncorrected(esdTrack); + param.GetCovFrom(esdTrack); + + // fill fTrackParamAtCluster with track parameters at each cluster if available + if(esdTrack.ClustersStored()) { + + // loop over ESD clusters + AliESDMuonCluster *esdCluster = (AliESDMuonCluster*) esdTrack.GetClusters().First(); + while (esdCluster) { + + // copy cluster information + AliMUONRawClusterV2 cluster(*esdCluster); + + // only set the Z parameter to avoid error in the AddTrackParamAtCluster(...) method + param.SetZ(cluster.GetZ()); + + // add common track parameters at current cluster + AddTrackParamAtCluster(param, cluster, kTRUE); + + esdCluster = (AliESDMuonCluster*) esdTrack.GetClusters().After(esdCluster); + } + + // sort array of track parameters at clusters + fTrackParamAtCluster->Sort(); + + // check that parameters stored in ESD are parameters at the most upstream cluster + // (convert Z position values to Float_t because of Double32_t in ESD track) + AliMUONTrackParam *firstTrackParam = (AliMUONTrackParam*) fTrackParamAtCluster->First(); + if (((Float_t)firstTrackParam->GetZ()) != ((Float_t)esdTrack.GetZUncorrected())) { + + AliError("track parameters are not given at the most upstream stored cluster"); + fTrackParamAtCluster->Clear("C"); + + } else { + + // Compute track parameters and covariances at each cluster from info at the first one + UpdateCovTrackParamAtCluster(); + + } + } - if (fTrackParamAtHit) { - // delete the TClonesArray of pointers to TrackParam - delete fTrackParamAtHit; - fTrackParamAtHit = NULL; + // fill fTrackParamAtCluster with track parameters at first (fake) cluster + // if first cluster not found or clusters not available + if (GetNClusters() == 0) { + + // get number of the first hit chamber (according to the MUONClusterMap if not empty) + Int_t firstCh = 0; + if (esdTrack.GetMuonClusterMap() != 0) while (!esdTrack.IsInMuonClusterMap(firstCh)) firstCh++; + else firstCh = AliMUONConstants::ChamberNumber(param.GetZ()); + + // produce fake cluster at this chamber + AliMUONRawClusterV2 fakeCluster(firstCh, 0, 0); + fakeCluster.SetXYZ(param.GetNonBendingCoor(), param.GetBendingCoor(), param.GetZ()); + fakeCluster.SetErrXY(0., 0.); + + // add track parameters at first (fake) cluster + AddTrackParamAtCluster(param, fakeCluster, kTRUE); + } + } - //__________________________________________________________________________ -AliMUONTrack::AliMUONTrack (const AliMUONTrack& theMUONTrack) - : TObject(theMUONTrack) +//__________________________________________________________________________ +AliMUONTrack::AliMUONTrack(const AliMUONTrack& track) + : TObject(track), + fTrackParamAtCluster(new TClonesArray("AliMUONTrackParam",10)), + 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) { - //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. + 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); + + // necessary to make a copy of the objects and not only the pointers in TClonesArray. + 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)); + } else if (fClusterWeightsNonBending) { + delete fClusterWeightsNonBending; + fClusterWeightsNonBending = 0x0; + } + + // copy cluster weights matrix if any + if (track.fClusterWeightsBending) { + if (fClusterWeightsBending) { + fClusterWeightsBending->ResizeTo(*(track.fClusterWeightsBending)); + *fClusterWeightsBending = *(track.fClusterWeightsBending); + } else fClusterWeightsBending = new TMatrixD(*(track.fClusterWeightsBending)); + } else if (fClusterWeightsBending) { + delete fClusterWeightsBending; + fClusterWeightsBending = 0x0; + } + + // copy track parameters at vertex if any + if (track.fTrackParamAtVertex) { + if (fTrackParamAtVertex) *fTrackParamAtVertex = *(track.fTrackParamAtVertex); + else fTrackParamAtVertex = new AliMUONTrackParam(*(track.fTrackParamAtVertex)); + } else if (fTrackParamAtVertex) { + delete fTrackParamAtVertex; + fTrackParamAtVertex = 0x0; + } + + 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; } //__________________________________________________________________________ -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) ???? - else { - cout << "ERROR in AliMUONTrack::SetFitMCS(FitMCS)" << endl; - cout << "FitMCS = " << FitMCS << " is neither 0 nor 1" << endl; - exit(0); - } - return; + /// Clear arrays + fTrackParamAtCluster->Clear(opt); + delete fClusterWeightsNonBending; fClusterWeightsNonBending = 0x0; + delete fClusterWeightsBending; fClusterWeightsBending = 0x0; + delete fTrackParamAtVertex; fTrackParamAtVertex = 0x0; } //__________________________________________________________________________ -void AliMUONTrack::SetFitNParam(Int_t FitNParam) +void AliMUONTrack::AddTrackParamAtCluster(const AliMUONTrackParam &trackParam, AliMUONVCluster &cluster, Bool_t copy) { - // 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); + /// 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; } - 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 + 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); } //__________________________________________________________________________ -void AliMUONTrack::SetFitStart(Int_t FitStart) +void AliMUONTrack::RemoveTrackParamAtCluster(AliMUONTrackParam *trackParam) { - // 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); + /// Remove trackParam from the array of TrackParamAtCluster + if (!fTrackParamAtCluster->Remove(trackParam)) { + AliWarning("object to remove does not exist in array fTrackParamAtCluster"); + return; } - return; + + fTrackParamAtCluster->Compress(); } //__________________________________________________________________________ -AliMUONTrackParam* AliMUONTrack::GetTrackParamAtFirstHit(void) const { - // Get pointer to TrackParamAtFirstHit - return ((AliMUONTrackHit*) (fTrackHitsPtr->First()))->GetTrackParam();} - - //__________________________________________________________________________ -void AliMUONTrack::RecursiveDump(void) const +void AliMUONTrack::UpdateTrackParamAtCluster() { - // 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(); + /// 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)); + } + } //__________________________________________________________________________ -Int_t AliMUONTrack::HitsInCommon(AliMUONTrack* Track) const +void AliMUONTrack::UpdateCovTrackParamAtCluster() { - // 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 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; + } + + 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)); + } + } //__________________________________________________________________________ -void AliMUONTrack::MatchTriggerTrack(TClonesArray *triggerTrackArray) +Bool_t AliMUONTrack::ComputeLocalChi2(Bool_t accountForMCS) { - // 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) + /// 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 + + 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); } - triggerTrack = (AliMUONTriggerTrack*) triggerTrackArray->After(triggerTrack); - } + + // 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); - // 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); + /// 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 + + 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 + + // 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 + + // 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) const +void AliMUONTrack::ComputeMCSCovariances(TMatrixD& mcsCovariances) const { - // 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 multiple scattering covariance matrix + /// (assume that track parameters at each cluster are corrects) + + // 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". + 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]); - } - - 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; - } + /// 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; - // 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)); - } - } - } 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()); + // 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); } - delete covBending; - delete covNonBending; - delete [] msa2; + return compatibleCluster; } -Double_t MultipleScatteringAngle2(AliMUONTrackHit *TrackHit) +//__________________________________________________________________________ +void AliMUONTrack::SetTrackParamAtVertex(const AliMUONTrackParam* trackParam) { - // 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; + /// set track parameters at vertex + if (trackParam == 0x0) return; + if (fTrackParamAtVertex) *fTrackParamAtVertex = *trackParam; + else fTrackParamAtVertex = new AliMUONTrackParam(*trackParam); } -//______________________________________________________________________________ - void mnvertLocal(Double_t *a, Int_t l, Int_t, Int_t n, Int_t &ifail) +//__________________________________________________________________________ +void AliMUONTrack::RecursiveDump() 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]; } - } - } -//*-*- 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]; - } - } - delete localVERTs; - delete localVERTq; - delete localVERTpp; - return; -//*-*- failure return -L100: - delete localVERTs; - delete localVERTq; - delete localVERTpp; - ifail = 1; -} /* mnvertLocal */ + /// 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) << GetLoTrgNum() << + ", Chi2-tracking-trigger=" << setw(8) << setprecision(5) << GetChi2MatchTrigger(); + cout << Form(" HitTriggerPattern %x",fHitsPatternInTrigCh) << endl; + 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) +{ + /// 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; + +}