1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
18 //-----------------------------------------------------------------------------
21 // Reconstructed track in ALICE dimuon spectrometer
22 //-----------------------------------------------------------------------------
24 #include "AliMUONTrack.h"
26 #include "AliMUONTrackParam.h"
27 #include "AliMUONHitForRec.h"
28 #include "AliMUONObjectPair.h"
29 #include "AliMUONConstants.h"
30 #include "AliMUONTrackExtrap.h"
37 #include <Riostream.h>
40 ClassImp(AliMUONTrack) // Class implementation in ROOT context
43 //__________________________________________________________________________
44 AliMUONTrack::AliMUONTrack()
46 fTrackParamAtVertex(),
47 fTrackParamAtHit(0x0),
50 fFitWithVertex(kFALSE),
53 fHitWeightsNonBending(0x0),
54 fHitWeightsBending(0x0),
59 fChi2MatchTrigger(0.),
61 fHitsPatternInTrigCh(0),
64 /// Default constructor
67 //__________________________________________________________________________
68 AliMUONTrack::AliMUONTrack(AliMUONObjectPair *segment)
70 fTrackParamAtVertex(),
71 fTrackParamAtHit(0x0),
74 fFitWithVertex(kFALSE),
77 fHitWeightsNonBending(0x0),
78 fHitWeightsBending(0x0),
83 fChi2MatchTrigger(0.),
85 fHitsPatternInTrigCh(0),
88 /// Constructor from thw hitForRec's
90 fTrackParamAtHit = new TClonesArray("AliMUONTrackParam",10);
91 fTrackParamAtHit->SetOwner(kTRUE);
92 fHitForRecAtHit = new TClonesArray("AliMUONHitForRec",10);
93 fHitForRecAtHit->SetOwner(kTRUE);
95 if (!segment) return; //AZ
97 // Pointers to hits from the segment
98 AliMUONHitForRec* hit1 = (AliMUONHitForRec*) segment->First();
99 AliMUONHitForRec* hit2 = (AliMUONHitForRec*) segment->Second();
101 // check sorting in -Z (spectro z<0)
102 if (hit1->GetZ() < hit2->GetZ()) {
104 hit2 = (AliMUONHitForRec*) segment->First();
107 // order the hits into the track according to the station the segment belong to
108 //(the hit first attached is the one from which we will start the tracking procedure)
109 if (hit1->GetChamberNumber() == 8) {
110 AddTrackParamAtHit(0,hit1);
111 AddTrackParamAtHit(0,hit2);
113 AddTrackParamAtHit(0,hit2);
114 AddTrackParamAtHit(0,hit1);
117 AliMUONTrackParam* trackParamAtFirstHit = (AliMUONTrackParam*) fTrackParamAtHit->First();
118 AliMUONHitForRec* firstHit = trackParamAtFirstHit->GetHitForRecPtr();
119 AliMUONTrackParam* trackParamAtLastHit = (AliMUONTrackParam*) fTrackParamAtHit->Last();
120 AliMUONHitForRec* lastHit = trackParamAtLastHit->GetHitForRecPtr();
123 // Compute track parameters
124 Double_t dZ = firstHit->GetZ() - lastHit->GetZ();
126 Double_t nonBendingCoor1 = firstHit->GetNonBendingCoor();
127 Double_t nonBendingCoor2 = lastHit->GetNonBendingCoor();
128 Double_t nonBendingSlope = (nonBendingCoor1 - nonBendingCoor2) / dZ;
130 Double_t bendingCoor1 = firstHit->GetBendingCoor();
131 Double_t bendingCoor2 = lastHit->GetBendingCoor();
132 Double_t bendingSlope = (bendingCoor1 - bendingCoor2) / dZ;
133 // Inverse bending momentum
134 Double_t bendingImpact = bendingCoor1 - firstHit->GetZ() * bendingSlope;
135 Double_t inverseBendingMomentum = 1. / AliMUONTrackExtrap::GetBendingMomentumFromImpactParam(bendingImpact);
138 // Set track parameters at first hit
139 trackParamAtFirstHit->SetNonBendingCoor(nonBendingCoor1);
140 trackParamAtFirstHit->SetNonBendingSlope(nonBendingSlope);
141 trackParamAtFirstHit->SetBendingCoor(bendingCoor1);
142 trackParamAtFirstHit->SetBendingSlope(bendingSlope);
143 trackParamAtFirstHit->SetInverseBendingMomentum(inverseBendingMomentum);
146 // Set track parameters at last hit
147 trackParamAtLastHit->SetNonBendingCoor(nonBendingCoor2);
148 trackParamAtLastHit->SetNonBendingSlope(nonBendingSlope);
149 trackParamAtLastHit->SetBendingCoor(bendingCoor2);
150 trackParamAtLastHit->SetBendingSlope(bendingSlope);
151 trackParamAtLastHit->SetInverseBendingMomentum(inverseBendingMomentum);
154 // Compute and set track parameters covariances at first hit
155 TMatrixD paramCov1(5,5);
158 paramCov1(0,0) = firstHit->GetNonBendingReso2();
159 paramCov1(0,1) = firstHit->GetNonBendingReso2() / dZ;
160 paramCov1(1,0) = paramCov1(0,1);
161 paramCov1(1,1) = ( firstHit->GetNonBendingReso2() + lastHit->GetNonBendingReso2() ) / dZ / dZ;
163 paramCov1(2,2) = firstHit->GetBendingReso2();
164 paramCov1(2,3) = firstHit->GetBendingReso2() / dZ;
165 paramCov1(3,2) = paramCov1(2,3);
166 paramCov1(3,3) = ( firstHit->GetBendingReso2() + lastHit->GetBendingReso2() ) / dZ / dZ;
167 // Inverse bending momentum (50% error)
168 paramCov1(4,4) = 0.5*inverseBendingMomentum * 0.5*inverseBendingMomentum;
170 trackParamAtFirstHit->SetCovariances(paramCov1);
173 // Compute and set track parameters covariances at last hit (as if the first hit did not exist)
174 TMatrixD paramCov2(5,5);
177 paramCov2(0,0) = paramCov1(0,0);
178 paramCov2(1,1) = 100.*paramCov1(1,1);
180 paramCov2(2,2) = paramCov1(2,2);
181 paramCov2(3,3) = 100.*paramCov1(3,3);
182 // Inverse bending momentum
183 paramCov2(4,4) = paramCov1(4,4);
185 trackParamAtLastHit->SetCovariances(paramCov2);
188 // Flag first hit as being removable
189 trackParamAtFirstHit->SetRemovable(kTRUE);
191 // Flag last hit as being removable
192 trackParamAtLastHit->SetRemovable(kTRUE);
196 //__________________________________________________________________________
197 AliMUONTrack::AliMUONTrack (const AliMUONTrack& track)
199 fTrackParamAtVertex(track.fTrackParamAtVertex),
200 fTrackParamAtHit(0x0),
201 fHitForRecAtHit(0x0),
202 fNTrackHits(track.fNTrackHits),
203 fFitWithVertex(track.fFitWithVertex),
205 fFitWithMCS(track.fFitWithMCS),
206 fHitWeightsNonBending(0x0),
207 fHitWeightsBending(0x0),
208 fGlobalChi2(track.fGlobalChi2),
209 fImproved(track.fImproved),
210 fMatchTrigger(track.fMatchTrigger),
211 floTrgNum(track.floTrgNum),
212 fChi2MatchTrigger(track.fChi2MatchTrigger),
213 fTrackID(track.fTrackID),
214 fHitsPatternInTrigCh(track.fHitsPatternInTrigCh),
215 fLocalTrigger(track.fLocalTrigger)
220 // necessary to make a copy of the objects and not only the pointers in TClonesArray.
221 if (track.fTrackParamAtHit) {
222 maxIndex = (track.fTrackParamAtHit)->GetEntriesFast();
223 fTrackParamAtHit = new TClonesArray("AliMUONTrackParam",maxIndex);
224 for (Int_t index = 0; index < maxIndex; index++) {
225 new ((*fTrackParamAtHit)[index]) AliMUONTrackParam(*(AliMUONTrackParam*)track.fTrackParamAtHit->At(index));
229 // necessary to make a copy of the objects and not only the pointers in TClonesArray.
230 if (track.fHitForRecAtHit) {
231 maxIndex = (track.fHitForRecAtHit)->GetEntriesFast();
232 fHitForRecAtHit = new TClonesArray("AliMUONHitForRec",maxIndex);
233 for (Int_t index = 0; index < maxIndex; index++) {
234 new ((*fHitForRecAtHit)[index]) AliMUONHitForRec(*(AliMUONHitForRec*)track.fHitForRecAtHit->At(index));
238 // copy vertex used during the tracking procedure if any
239 if (track.fVertex) fVertex = new AliMUONHitForRec(*(track.fVertex));
241 // copy hit weights matrices if any
242 if (track.fHitWeightsNonBending) fHitWeightsNonBending = new TMatrixD(*(track.fHitWeightsNonBending));
243 if (track.fHitWeightsBending) fHitWeightsBending = new TMatrixD(*(track.fHitWeightsBending));
247 //__________________________________________________________________________
248 AliMUONTrack & AliMUONTrack::operator=(const AliMUONTrack& track)
250 /// Asignment operator
251 // check assignement to self
255 // base class assignement
256 TObject::operator=(track);
258 fTrackParamAtVertex = track.fTrackParamAtVertex;
262 // necessary to make a copy of the objects and not only the pointers in TClonesArray.
263 if (track.fTrackParamAtHit) {
264 if (fTrackParamAtHit) fTrackParamAtHit->Clear();
265 else fTrackParamAtHit = new TClonesArray("AliMUONTrackParam",10);
266 maxIndex = (track.fTrackParamAtHit)->GetEntriesFast();
267 for (Int_t index = 0; index < maxIndex; index++) {
268 new ((*fTrackParamAtHit)[fTrackParamAtHit->GetEntriesFast()])
269 AliMUONTrackParam(*(AliMUONTrackParam*)(track.fTrackParamAtHit)->At(index));
271 } else if (fTrackParamAtHit) {
272 delete fTrackParamAtHit;
273 fTrackParamAtHit = 0x0;
276 // necessary to make a copy of the objects and not only the pointers in TClonesArray.
277 if (track.fHitForRecAtHit) {
278 if (fHitForRecAtHit) fHitForRecAtHit->Clear();
279 else fHitForRecAtHit = new TClonesArray("AliMUONHitForRec",10);
280 maxIndex = (track.fHitForRecAtHit)->GetEntriesFast();
281 for (Int_t index = 0; index < maxIndex; index++) {
282 new ((*fHitForRecAtHit)[fHitForRecAtHit->GetEntriesFast()])
283 AliMUONHitForRec(*(AliMUONHitForRec*)(track.fHitForRecAtHit)->At(index));
285 } else if (fHitForRecAtHit) {
286 delete fHitForRecAtHit;
287 fHitForRecAtHit = 0x0;
290 // copy vertex used during the tracking procedure if any.
292 if (fVertex) *fVertex = *(track.fVertex);
293 else fVertex = new AliMUONHitForRec(*(track.fVertex));
294 } else if (fVertex) {
299 // copy hit weights matrix if any
300 if (track.fHitWeightsNonBending) {
301 if (fHitWeightsNonBending) {
302 fHitWeightsNonBending->ResizeTo(*(track.fHitWeightsNonBending));
303 *fHitWeightsNonBending = *(track.fHitWeightsNonBending);
304 } else fHitWeightsNonBending = new TMatrixD(*(track.fHitWeightsNonBending));
305 } else if (fHitWeightsNonBending) {
306 delete fHitWeightsNonBending;
307 fHitWeightsNonBending = 0x0;
310 // copy hit weights matrix if any
311 if (track.fHitWeightsBending) {
312 if (fHitWeightsBending) {
313 fHitWeightsBending->ResizeTo(*(track.fHitWeightsBending));
314 *fHitWeightsBending = *(track.fHitWeightsBending);
315 } else fHitWeightsBending = new TMatrixD(*(track.fHitWeightsBending));
316 } else if (fHitWeightsBending) {
317 delete fHitWeightsBending;
318 fHitWeightsBending = 0x0;
321 fNTrackHits = track.fNTrackHits;
322 fFitWithVertex = track.fFitWithVertex;
323 fFitWithMCS = track.fFitWithMCS;
324 fGlobalChi2 = track.fGlobalChi2;
325 fImproved = track.fImproved;
326 fMatchTrigger = track.fMatchTrigger;
327 floTrgNum = track.floTrgNum;
328 fChi2MatchTrigger = track.fChi2MatchTrigger;
329 fTrackID = track.fTrackID;
330 fHitsPatternInTrigCh = track.fHitsPatternInTrigCh;
331 fLocalTrigger = track.fLocalTrigger;
336 //__________________________________________________________________________
337 AliMUONTrack::~AliMUONTrack()
340 delete fTrackParamAtHit;
341 delete fHitForRecAtHit;
343 delete fHitWeightsNonBending;
344 delete fHitWeightsBending;
347 //__________________________________________________________________________
348 void AliMUONTrack::Clear(Option_t* opt)
351 if ( fTrackParamAtHit ) fTrackParamAtHit->Clear(opt);
352 if ( fHitForRecAtHit ) fHitForRecAtHit->Clear(opt);
353 delete fVertex; fVertex = 0x0;
354 delete fHitWeightsNonBending; fHitWeightsNonBending = 0x0;
355 delete fHitWeightsBending; fHitWeightsBending = 0x0;
358 //__________________________________________________________________________
359 void AliMUONTrack::AddTrackParamAtHit(const AliMUONTrackParam *trackParam, AliMUONHitForRec *hitForRec)
361 /// Add TrackParamAtHit if "trackParam" != NULL
362 /// else create empty TrackParamAtHit and set the z position to the one of "hitForRec" if any
363 /// Update link to HitForRec if "hitForRec" != NULL
364 if (!fTrackParamAtHit) {
365 fTrackParamAtHit = new TClonesArray("AliMUONTrackParam",10);
368 AliMUONTrackParam* trackParamAtHit;
370 trackParamAtHit = new ((*fTrackParamAtHit)[fNTrackHits]) AliMUONTrackParam(*trackParam);
372 if (hitForRec->GetZ() != trackParam->GetZ())
373 AliWarning("Added track parameters at a different z position than the one of the attached hit");
376 trackParamAtHit = new ((*fTrackParamAtHit)[fNTrackHits]) AliMUONTrackParam();
377 if (hitForRec) trackParamAtHit->SetZ(hitForRec->GetZ());
379 if (hitForRec) trackParamAtHit->SetHitForRecPtr(hitForRec);
383 //__________________________________________________________________________
384 void AliMUONTrack::RemoveTrackParamAtHit(AliMUONTrackParam *trackParam)
386 /// Remove trackParam from the array of TrackParamAtHit
387 if (!fTrackParamAtHit) {
388 AliWarning("array fTrackParamAtHit does not exist");
392 if (!fTrackParamAtHit->Remove(trackParam)) {
393 AliWarning("object to remove does not exist in array fTrackParamAtHit");
397 fTrackParamAtHit->Compress();
401 //__________________________________________________________________________
402 void AliMUONTrack::AddHitForRecAtHit(const AliMUONHitForRec *hitForRec)
404 /// Add hitForRec to the array of hitForRec at hit
405 if (!fHitForRecAtHit)
406 fHitForRecAtHit = new TClonesArray("AliMUONHitForRec",10);
409 AliFatal("AliMUONTrack::AddHitForRecAtHit: hitForRec == NULL");
411 new ((*fHitForRecAtHit)[fHitForRecAtHit->GetEntriesFast()]) AliMUONHitForRec(*hitForRec);
414 //__________________________________________________________________________
415 void AliMUONTrack::UpdateTrackParamAtHit()
417 /// Update track parameters at each attached hit
419 if (fNTrackHits == 0) {
420 AliWarning("no hit attached to the track");
425 AliMUONTrackParam* startingTrackParam = (AliMUONTrackParam*) fTrackParamAtHit->First();
426 AliMUONTrackParam* trackParamAtHit = (AliMUONTrackParam*) fTrackParamAtHit->After(startingTrackParam);
427 while (trackParamAtHit) {
430 z = trackParamAtHit->GetZ();
432 // reset track parameters and their covariances
433 trackParamAtHit->SetParameters(startingTrackParam->GetParameters());
434 trackParamAtHit->SetZ(startingTrackParam->GetZ());
436 // extrapolation to the given z
437 AliMUONTrackExtrap::ExtrapToZ(trackParamAtHit, z);
440 startingTrackParam = trackParamAtHit;
441 trackParamAtHit = (AliMUONTrackParam*) (fTrackParamAtHit->After(trackParamAtHit));
446 //__________________________________________________________________________
447 void AliMUONTrack::UpdateCovTrackParamAtHit()
449 /// Update track parameters and their covariances at each attached hit
451 if (fNTrackHits == 0) {
452 AliWarning("no hit attached to the track");
457 AliMUONTrackParam* startingTrackParam = (AliMUONTrackParam*) fTrackParamAtHit->First();
458 AliMUONTrackParam* trackParamAtHit = (AliMUONTrackParam*) fTrackParamAtHit->After(startingTrackParam);
459 while (trackParamAtHit) {
462 z = trackParamAtHit->GetZ();
464 // reset track parameters and their covariances
465 trackParamAtHit->SetParameters(startingTrackParam->GetParameters());
466 trackParamAtHit->SetZ(startingTrackParam->GetZ());
467 trackParamAtHit->SetCovariances(startingTrackParam->GetCovariances());
469 // extrapolation to the given z
470 AliMUONTrackExtrap::ExtrapToZCov(trackParamAtHit, z);
473 startingTrackParam = trackParamAtHit;
474 trackParamAtHit = (AliMUONTrackParam*) (fTrackParamAtHit->After(trackParamAtHit));
479 //__________________________________________________________________________
480 void AliMUONTrack::SetVertex(const AliMUONHitForRec* vertex)
482 /// Set the vertex used during the tracking procedure
483 if (!fVertex) fVertex = new AliMUONHitForRec(*vertex);
484 else *fVertex = *vertex;
488 //__________________________________________________________________________
489 Bool_t AliMUONTrack::ComputeLocalChi2(Bool_t accountForMCS)
491 /// Compute the removable hit contribution to the chi2 of the track
492 /// accounting for multiple scattering or not according to the flag
493 /// - Also recompute the weight matrices of the attached hits if accountForMCS=kTRUE
494 /// - Assume that track parameters at each hit are corrects
495 /// - Return kFALSE if computation failed
497 // Check hits (if the first one exist, assume that the other ones exit too!)
498 AliMUONTrackParam* trackParamAtHit = (AliMUONTrackParam*) fTrackParamAtHit->First();
499 if (!trackParamAtHit->GetHitForRecPtr()) {
500 AliWarning("hit is missing");
504 if (accountForMCS) { // Compute local chi2 taking into account multiple scattering effects
506 // Compute MCS covariance matrix only once
507 TMatrixD mcsCovariances(AliMUONConstants::NTrackingCh(),AliMUONConstants::NTrackingCh());
508 ComputeMCSCovariances(mcsCovariances);
510 // Make sure hit weights are consistent with following calculations
511 if (!ComputeHitWeights(&mcsCovariances)) {
512 AliWarning("cannot take into account the multiple scattering effects");
513 return ComputeLocalChi2(kFALSE);
516 // Compute chi2 of the track
517 Double_t globalChi2 = ComputeGlobalChi2(kTRUE);
518 if (globalChi2 < 0.) return kFALSE;
520 // Loop over removable hits and compute their local chi2
521 AliMUONTrackParam* trackParamAtHit1;
522 AliMUONHitForRec *hitForRec, *discardedHit;
523 Int_t hitNumber1, hitNumber2, currentHitNumber1, currentHitNumber2;
524 TMatrixD hitWeightsNB(fNTrackHits-1,fNTrackHits-1);
525 TMatrixD hitWeightsB(fNTrackHits-1,fNTrackHits-1);
526 Double_t *dX = new Double_t[fNTrackHits-1];
527 Double_t *dY = new Double_t[fNTrackHits-1];
528 Double_t globalChi2b;
529 while (trackParamAtHit) {
531 discardedHit = trackParamAtHit->GetHitForRecPtr();
533 // Recompute hit weights without the current hit
534 if (!ComputeHitWeights(hitWeightsNB, hitWeightsB, &mcsCovariances, discardedHit)) {
535 AliWarning("cannot take into account the multiple scattering effects");
536 ComputeLocalChi2(kFALSE);
539 // Compute track chi2 without the current hit
541 currentHitNumber1 = 0;
542 for (hitNumber1 = 0; hitNumber1 < fNTrackHits ; hitNumber1++) {
543 trackParamAtHit1 = (AliMUONTrackParam*) fTrackParamAtHit->UncheckedAt(hitNumber1);
544 hitForRec = trackParamAtHit1->GetHitForRecPtr();
546 if (hitForRec == discardedHit) continue;
548 // Compute and save residuals
549 dX[currentHitNumber1] = hitForRec->GetNonBendingCoor() - trackParamAtHit1->GetNonBendingCoor();
550 dY[currentHitNumber1] = hitForRec->GetBendingCoor() - trackParamAtHit1->GetBendingCoor();
552 currentHitNumber2 = 0;
553 for (hitNumber2 = 0; hitNumber2 < hitNumber1; hitNumber2++) {
554 hitForRec = ((AliMUONTrackParam*) fTrackParamAtHit->UncheckedAt(hitNumber2))->GetHitForRecPtr();
556 if (hitForRec == discardedHit) continue;
558 // Add contribution from covariances
559 globalChi2b += (hitWeightsNB(currentHitNumber1, currentHitNumber2) +
560 hitWeightsNB(currentHitNumber2, currentHitNumber1)) * dX[currentHitNumber1] * dX[currentHitNumber2] +
561 (hitWeightsB(currentHitNumber1, currentHitNumber2) +
562 hitWeightsB(currentHitNumber2, currentHitNumber1)) * dY[currentHitNumber1] * dY[currentHitNumber2];
567 // Add contribution from variances
568 globalChi2b += hitWeightsNB(currentHitNumber1, currentHitNumber1) * dX[currentHitNumber1] * dX[currentHitNumber1] +
569 hitWeightsB(currentHitNumber1, currentHitNumber1) * dY[currentHitNumber1] * dY[currentHitNumber1];
575 trackParamAtHit->SetLocalChi2(globalChi2 - globalChi2b);
577 trackParamAtHit = (AliMUONTrackParam*) fTrackParamAtHit->After(trackParamAtHit);
583 } else { // without multiple scattering effects
585 AliMUONHitForRec *discardedHit;
587 while (trackParamAtHit) {
589 discardedHit = trackParamAtHit->GetHitForRecPtr();
592 dX = discardedHit->GetNonBendingCoor() - trackParamAtHit->GetNonBendingCoor();
593 dY = discardedHit->GetBendingCoor() - trackParamAtHit->GetBendingCoor();
596 trackParamAtHit->SetLocalChi2(dX * dX / discardedHit->GetNonBendingReso2() + dY * dY / discardedHit->GetBendingReso2());
598 trackParamAtHit = (AliMUONTrackParam*) fTrackParamAtHit->After(trackParamAtHit);
608 //__________________________________________________________________________
609 Double_t AliMUONTrack::ComputeGlobalChi2(Bool_t accountForMCS)
611 /// Compute the chi2 of the track accounting for multiple scattering or not according to the flag
612 /// - Assume that track parameters at each hit are corrects
613 /// - Assume the hits weights matrices are corrects
614 /// - Return negative value if chi2 computation failed
616 // Check hits (if the first one exist, assume that the other ones exit too!)
617 AliMUONTrackParam* trackParamAtHit = (AliMUONTrackParam*) fTrackParamAtHit->First();
618 if (!trackParamAtHit->GetHitForRecPtr()) {
619 AliWarning("hit is missing");
627 // Check the weight matrices
628 Bool_t weightsAvailable = kTRUE;
629 if (!fHitWeightsNonBending || !fHitWeightsBending) weightsAvailable = kFALSE;
630 else if (fHitWeightsNonBending->GetNrows() != fNTrackHits || fHitWeightsNonBending->GetNcols() != fNTrackHits ||
631 fHitWeightsBending->GetNrows() != fNTrackHits || fHitWeightsBending->GetNcols() != fNTrackHits) weightsAvailable = kFALSE;
633 // if weight matrices are not available compute chi2 without MCS
634 if (!weightsAvailable) {
635 AliWarning("hit weights including multiple scattering effects are not available\n\t\t --> compute chi2 WITHOUT multiple scattering");
636 return ComputeGlobalChi2(kFALSE);
640 AliMUONHitForRec *hitForRec;
641 Double_t *dX = new Double_t[fNTrackHits];
642 Double_t *dY = new Double_t[fNTrackHits];
643 Int_t hitNumber1, hitNumber2;
644 for (hitNumber1 = 0; hitNumber1 < fNTrackHits ; hitNumber1++) {
645 trackParamAtHit = (AliMUONTrackParam*) fTrackParamAtHit->UncheckedAt(hitNumber1);
646 hitForRec = trackParamAtHit->GetHitForRecPtr();
647 dX[hitNumber1] = hitForRec->GetNonBendingCoor() - trackParamAtHit->GetNonBendingCoor();
648 dY[hitNumber1] = hitForRec->GetBendingCoor() - trackParamAtHit->GetBendingCoor();
649 for (hitNumber2 = 0; hitNumber2 < hitNumber1; hitNumber2++) {
650 chi2 += ((*fHitWeightsNonBending)(hitNumber1, hitNumber2) + (*fHitWeightsNonBending)(hitNumber2, hitNumber1)) * dX[hitNumber1] * dX[hitNumber2] +
651 ((*fHitWeightsBending)(hitNumber1, hitNumber2) + (*fHitWeightsBending)(hitNumber2, hitNumber1)) * dY[hitNumber1] * dY[hitNumber2];
653 chi2 += ((*fHitWeightsNonBending)(hitNumber1, hitNumber1) * dX[hitNumber1] * dX[hitNumber1]) +
654 ((*fHitWeightsBending)(hitNumber1, hitNumber1) * dY[hitNumber1] * dY[hitNumber1]);
661 AliMUONHitForRec *hitForRec;
663 for (Int_t hitNumber = 0; hitNumber < fNTrackHits ; hitNumber++) {
664 trackParamAtHit = (AliMUONTrackParam*) fTrackParamAtHit->UncheckedAt(hitNumber);
665 hitForRec = trackParamAtHit->GetHitForRecPtr();
666 dX = hitForRec->GetNonBendingCoor() - trackParamAtHit->GetNonBendingCoor();
667 dY = hitForRec->GetBendingCoor() - trackParamAtHit->GetBendingCoor();
668 chi2 += dX * dX / hitForRec->GetNonBendingReso2() + dY * dY / hitForRec->GetBendingReso2();
677 //__________________________________________________________________________
678 Bool_t AliMUONTrack::ComputeHitWeights(TMatrixD* mcsCovariances)
680 /// Compute the weight matrices of the attached hits, in non bending and bending direction,
681 /// accounting for multiple scattering correlations and hits resolution
682 /// - Use the provided MCS covariance matrix if any (otherwise build it temporarily)
683 /// - Assume that track parameters at each hit are corrects
684 /// - Return kFALSE if computation failed
687 if (!fHitWeightsNonBending) fHitWeightsNonBending = new TMatrixD(fNTrackHits,fNTrackHits);
688 if (!fHitWeightsBending) fHitWeightsBending = new TMatrixD(fNTrackHits,fNTrackHits);
690 // Check hits (if the first one exist, assume that the other ones exit too!)
691 if (!((AliMUONTrackParam*) fTrackParamAtHit->First())->GetHitForRecPtr()) {
692 AliWarning("hit is missing");
693 fHitWeightsNonBending->ResizeTo(0,0);
694 fHitWeightsBending->ResizeTo(0,0);
698 // Compute weights matrices
699 if (!ComputeHitWeights(*fHitWeightsNonBending, *fHitWeightsBending, mcsCovariances)) return kFALSE;
705 //__________________________________________________________________________
706 Bool_t AliMUONTrack::ComputeHitWeights(TMatrixD& hitWeightsNB, TMatrixD& hitWeightsB, TMatrixD* mcsCovariances, AliMUONHitForRec* discardedHit) const
708 /// Compute the weight matrices, in non bending and bending direction,
709 /// of the other attached hits assuming the discarded one does not exist
710 /// accounting for multiple scattering correlations and hits resolution
711 /// - Use the provided MCS covariance matrix if any (otherwise build it temporarily)
712 /// - Return kFALSE if computation failed
714 // Check MCS covariance matrix and recompute it if need
715 Bool_t deleteMCSCov = kFALSE;
716 if (!mcsCovariances) {
718 // build MCS covariance matrix
719 mcsCovariances = new TMatrixD(AliMUONConstants::NTrackingCh(),AliMUONConstants::NTrackingCh());
720 deleteMCSCov = kTRUE;
721 ComputeMCSCovariances(*mcsCovariances);
725 // check MCS covariance matrix size
726 if (mcsCovariances->GetNrows() != AliMUONConstants::NTrackingCh() || mcsCovariances->GetNcols() != AliMUONConstants::NTrackingCh()) {
727 ComputeMCSCovariances(*mcsCovariances);
732 // Resize the weights matrices; alocate memory
734 hitWeightsNB.ResizeTo(fNTrackHits-1,fNTrackHits-1);
735 hitWeightsB.ResizeTo(fNTrackHits-1,fNTrackHits-1);
737 hitWeightsNB.ResizeTo(fNTrackHits,fNTrackHits);
738 hitWeightsB.ResizeTo(fNTrackHits,fNTrackHits);
742 AliMUONHitForRec *hitForRec1, *hitForRec2;
743 Int_t chamber1, chamber2, currentHitNumber1, currentHitNumber2;
745 // Compute the covariance matrices
746 currentHitNumber1 = 0;
747 for (Int_t hitNumber1 = 0; hitNumber1 < fNTrackHits; hitNumber1++) {
748 hitForRec1 = ((AliMUONTrackParam*) fTrackParamAtHit->UncheckedAt(hitNumber1))->GetHitForRecPtr();
750 if (hitForRec1 == discardedHit) continue;
752 chamber1 = hitForRec1->GetChamberNumber();
754 // Loop over next hits
755 currentHitNumber2 = currentHitNumber1;
756 for (Int_t hitNumber2 = hitNumber1; hitNumber2 < fNTrackHits; hitNumber2++) {
757 hitForRec2 = ((AliMUONTrackParam*) fTrackParamAtHit->UncheckedAt(hitNumber2))->GetHitForRecPtr();
759 if (hitForRec2 == discardedHit) continue;
761 chamber2 = hitForRec2->GetChamberNumber();
763 // Fill with MCS covariances
764 hitWeightsNB(currentHitNumber1, currentHitNumber2) = (*mcsCovariances)(chamber1,chamber2);
766 // Equal contribution from multiple scattering in non bending and bending directions
767 hitWeightsB(currentHitNumber1, currentHitNumber2) = hitWeightsNB(currentHitNumber1, currentHitNumber2);
769 // Add contribution from hit resolution to diagonal element and symmetrize the matrix
770 if (currentHitNumber1 == currentHitNumber2) {
772 // In non bending plane
773 hitWeightsNB(currentHitNumber1, currentHitNumber1) += hitForRec1->GetNonBendingReso2();
775 hitWeightsB(currentHitNumber1, currentHitNumber1) += hitForRec1->GetBendingReso2();
779 // In non bending plane
780 hitWeightsNB(currentHitNumber2, currentHitNumber1) = hitWeightsNB(currentHitNumber1, currentHitNumber2);
782 hitWeightsB(currentHitNumber2, currentHitNumber1) = hitWeightsB(currentHitNumber1, currentHitNumber2);
792 // Inversion of covariance matrices to get the weights
793 if (hitWeightsNB.Determinant() != 0 && hitWeightsB.Determinant() != 0) {
794 hitWeightsNB.Invert();
795 hitWeightsB.Invert();
797 AliWarning(" Determinant = 0");
798 hitWeightsNB.ResizeTo(0,0);
799 hitWeightsB.ResizeTo(0,0);
800 if(deleteMCSCov) delete mcsCovariances;
804 if(deleteMCSCov) delete mcsCovariances;
810 //__________________________________________________________________________
811 void AliMUONTrack::ComputeMCSCovariances(TMatrixD& mcsCovariances) const
813 /// Compute the multiple scattering covariance matrix
814 /// - Assume that track parameters at each hit are corrects
815 /// - Return kFALSE if computation failed
817 // Make sure the size of the covariance matrix is correct
818 Int_t nChambers = AliMUONConstants::NTrackingCh();
819 mcsCovariances.ResizeTo(nChambers,nChambers);
821 // check for too many track hits
822 if (fNTrackHits > nChambers) {
823 AliWarning("more than 1 hit per chamber!!");
824 mcsCovariances.Zero();
829 AliMUONTrackParam* trackParamAtHit;
830 AliMUONHitForRec *hitForRec;
831 AliMUONTrackParam extrapTrackParam;
832 Int_t currentChamber, expectedChamber;
833 Double_t *mcsAngle2 = new Double_t[nChambers];
834 Double_t *zMCS = new Double_t[nChambers];
836 // Compute multiple scattering dispersion angle at each chamber
837 // and save the z position where it is calculated
840 for (Int_t hitNumber = 0; hitNumber < fNTrackHits; hitNumber++) {
841 trackParamAtHit = (AliMUONTrackParam*) fTrackParamAtHit->UncheckedAt(hitNumber);
842 hitForRec = trackParamAtHit->GetHitForRecPtr();
844 // look for missing chambers if any
845 currentChamber = hitForRec->GetChamberNumber();
846 while (currentChamber > expectedChamber) {
848 // Save the z position where MCS dispersion is calculated
849 zMCS[expectedChamber] = AliMUONConstants::DefaultChamberZ(expectedChamber);
851 // Do not take into account MCS in chambers prior the first hit
854 // Get track parameters at missing chamber z
855 extrapTrackParam = *trackParamAtHit;
856 AliMUONTrackExtrap::ExtrapToZ(&extrapTrackParam, zMCS[expectedChamber]);
858 // Save multiple scattering dispersion angle in missing chamber
859 mcsAngle2[expectedChamber] = AliMUONTrackExtrap::GetMCSAngle2(extrapTrackParam,AliMUONConstants::ChamberThicknessInX0(),1.);
861 } else mcsAngle2[expectedChamber] = 0.;
866 // Save z position where MCS dispersion is calculated
867 zMCS[currentChamber] = trackParamAtHit->GetZ();
869 // Save multiple scattering dispersion angle in current chamber
870 mcsAngle2[currentChamber] = AliMUONTrackExtrap::GetMCSAngle2(*trackParamAtHit,AliMUONConstants::ChamberThicknessInX0(),1.);
875 // complete array of z if last hit is on the last but one chamber
876 if (currentChamber != nChambers-1) zMCS[nChambers-1] = AliMUONConstants::DefaultChamberZ(nChambers-1);
879 // Compute the covariance matrix
880 for (Int_t chamber1 = 0; chamber1 < nChambers; chamber1++) {
882 for (Int_t chamber2 = chamber1; chamber2 < nChambers; chamber2++) {
884 // Initialization to 0 (diagonal plus upper triangular part)
885 mcsCovariances(chamber1, chamber2) = 0.;
887 // Compute contribution from multiple scattering in upstream chambers
888 for (currentChamber = 0; currentChamber < chamber1; currentChamber++) {
889 mcsCovariances(chamber1, chamber2) += (zMCS[chamber1] - zMCS[currentChamber]) * (zMCS[chamber2] - zMCS[currentChamber]) * mcsAngle2[currentChamber];
892 // Symetrize the matrix
893 mcsCovariances(chamber2, chamber1) = mcsCovariances(chamber1, chamber2);
903 //__________________________________________________________________________
904 Int_t AliMUONTrack::HitsInCommon(AliMUONTrack* track) const
906 /// Returns the number of hits in common between the current track ("this")
907 /// and the track pointed to by "track".
908 Int_t hitsInCommon = 0;
909 AliMUONTrackParam *trackParamAtHit1, *trackParamAtHit2;
910 // Loop over hits of first track
911 trackParamAtHit1 = (AliMUONTrackParam*) this->fTrackParamAtHit->First();
912 while (trackParamAtHit1) {
913 // Loop over hits of second track
914 trackParamAtHit2 = (AliMUONTrackParam*) track->fTrackParamAtHit->First();
915 while (trackParamAtHit2) {
916 // Increment "hitsInCommon" if both TrackParamAtHits point to the same HitForRec
917 if ((trackParamAtHit1->GetHitForRecPtr()) == (trackParamAtHit2->GetHitForRecPtr())) {
921 trackParamAtHit2 = (AliMUONTrackParam*) track->fTrackParamAtHit->After(trackParamAtHit2);
922 } // trackParamAtHit2
923 trackParamAtHit1 = (AliMUONTrackParam*) this->fTrackParamAtHit->After(trackParamAtHit1);
924 } // trackParamAtHit1
928 //__________________________________________________________________________
929 Double_t AliMUONTrack::GetNormalizedChi2() const
931 /// return the chi2 value divided by the number of degrees of freedom (or 1.e10 if ndf < 0)
933 Double_t numberOfDegFree = (2. * fNTrackHits - 5.);
934 if (numberOfDegFree > 0.) return fGlobalChi2 / numberOfDegFree;
938 //__________________________________________________________________________
939 Bool_t* AliMUONTrack::CompatibleTrack(AliMUONTrack * track, Double_t sigma2Cut) const
941 /// Return kTRUE/kFALSE for each chamber if hit is compatible or not
942 TClonesArray *hitArray, *thisHitArray;
943 AliMUONHitForRec *hit, *thisHit;
946 Float_t deltaZMax = 1.; // 1 cm
948 Bool_t *nCompHit = new Bool_t[AliMUONConstants::NTrackingCh()];
950 for ( Int_t ch = 0; ch < AliMUONConstants::NTrackingCh(); ch++) {
951 nCompHit[ch] = kFALSE;
954 thisHitArray = this->GetHitForRecAtHit();
956 hitArray = track->GetHitForRecAtHit();
958 for (Int_t iHthis = 0; iHthis < thisHitArray->GetEntriesFast(); iHthis++) {
959 thisHit = (AliMUONHitForRec*) thisHitArray->At(iHthis);
960 chamberNumber = thisHit->GetChamberNumber();
961 if (chamberNumber < 0 || chamberNumber > AliMUONConstants::NTrackingCh()) continue;
962 nCompHit[chamberNumber] = kFALSE;
963 for (Int_t iH = 0; iH < hitArray->GetEntriesFast(); iH++) {
964 hit = (AliMUONHitForRec*) hitArray->At(iH);
965 deltaZ = TMath::Abs(thisHit->GetZ() - hit->GetZ());
966 chi2 = thisHit->NormalizedChi2WithHitForRec(hit,sigma2Cut); // set cut to 4 sigmas
967 if (chi2 < 3. && deltaZ < deltaZMax) {
968 nCompHit[chamberNumber] = kTRUE;
977 //__________________________________________________________________________
978 void AliMUONTrack::RecursiveDump(void) const
980 /// Recursive dump of AliMUONTrack, i.e. with dump of TrackParamAtHit's and attached HitForRec's
981 AliMUONTrackParam *trackParamAtHit;
982 AliMUONHitForRec *hitForRec;
983 cout << "Recursive dump of Track: " << this << endl;
986 for (Int_t trackHitIndex = 0; trackHitIndex < fNTrackHits; trackHitIndex++) {
987 trackParamAtHit = (AliMUONTrackParam*) ((*fTrackParamAtHit)[trackHitIndex]);
989 cout << "TrackParamAtHit: " << trackParamAtHit << " (index: " << trackHitIndex << ")" << endl;
990 trackParamAtHit->Dump();
991 hitForRec = trackParamAtHit->GetHitForRecPtr();
993 cout << "HitForRec: " << hitForRec << endl;
999 //_____________________________________________-
1000 void AliMUONTrack::Print(Option_t*) const
1002 /// Printing Track information
1004 cout << "<AliMUONTrack> No.Clusters=" << setw(2) << GetNTrackHits() <<
1005 ", Match2Trig=" << setw(1) << GetMatchTrigger() <<
1006 ", LoTrgNum=" << setw(3) << GetLoTrgNum() <<
1007 ", Chi2-tracking-trigger=" << setw(8) << setprecision(5) << GetChi2MatchTrigger();
1008 cout << Form(" HitTriggerPattern %x",fHitsPatternInTrigCh) << endl;
1009 GetTrackParamAtHit()->First()->Print("FULL");
1012 //__________________________________________________________________________
1013 void AliMUONTrack::SetLocalTrigger(Int_t loCirc, Int_t loStripX, Int_t loStripY, Int_t loDev, Int_t loLpt, Int_t loHpt)
1015 /// pack the local trigger information and store
1017 if (loCirc < 0 || loCirc > 233) return;
1020 fLocalTrigger += loCirc;
1021 fLocalTrigger += loStripX << 8;
1022 fLocalTrigger += loStripY << 13;
1023 fLocalTrigger += loDev << 17;
1024 fLocalTrigger += loLpt << 22;
1025 fLocalTrigger += loHpt << 24;