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(fNTrackHits,fNTrackHits);
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);
607 //__________________________________________________________________________
608 Double_t AliMUONTrack::ComputeGlobalChi2(Bool_t accountForMCS)
610 /// Compute the chi2 of the track accounting for multiple scattering or not according to the flag
611 /// - Assume that track parameters at each hit are corrects
612 /// - Assume the hits weights matrices are corrects
613 /// - Return negative value if chi2 computation failed
615 // Check hits (if the first one exist, assume that the other ones exit too!)
616 AliMUONTrackParam* trackParamAtHit = (AliMUONTrackParam*) fTrackParamAtHit->First();
617 if (!trackParamAtHit->GetHitForRecPtr()) {
618 AliWarning("hit is missing");
626 // Check the weight matrices. If weight matrices are not available compute chi2 without MCS
627 if (!fHitWeightsNonBending || !fHitWeightsBending) {
628 AliWarning("hit weights including multiple scattering effects are not available\n\t\t --> compute chi2 WITHOUT multiple scattering");
629 return ComputeGlobalChi2(kFALSE);
631 if (fHitWeightsNonBending->GetNrows() != fNTrackHits || fHitWeightsBending->GetNcols() != fNTrackHits) {
632 AliWarning("hit weights including multiple scattering effects are not available\n\t\t --> compute chi2 WITHOUT multiple scattering");
633 return ComputeGlobalChi2(kFALSE);
637 AliMUONHitForRec *hitForRec;
638 Double_t *dX = new Double_t[fNTrackHits];
639 Double_t *dY = new Double_t[fNTrackHits];
640 Int_t hitNumber1, hitNumber2;
641 for (hitNumber1 = 0; hitNumber1 < fNTrackHits ; hitNumber1++) {
642 trackParamAtHit = (AliMUONTrackParam*) fTrackParamAtHit->UncheckedAt(hitNumber1);
643 hitForRec = trackParamAtHit->GetHitForRecPtr();
644 dX[hitNumber1] = hitForRec->GetNonBendingCoor() - trackParamAtHit->GetNonBendingCoor();
645 dY[hitNumber1] = hitForRec->GetBendingCoor() - trackParamAtHit->GetBendingCoor();
646 for (hitNumber2 = 0; hitNumber2 < hitNumber1; hitNumber2++) {
647 chi2 += ((*fHitWeightsNonBending)(hitNumber1, hitNumber2) + (*fHitWeightsNonBending)(hitNumber2, hitNumber1)) * dX[hitNumber1] * dX[hitNumber2] +
648 ((*fHitWeightsBending)(hitNumber1, hitNumber2) + (*fHitWeightsBending)(hitNumber2, hitNumber1)) * dY[hitNumber1] * dY[hitNumber2];
650 chi2 += ((*fHitWeightsNonBending)(hitNumber1, hitNumber1) * dX[hitNumber1] * dX[hitNumber1]) +
651 ((*fHitWeightsBending)(hitNumber1, hitNumber1) * dY[hitNumber1] * dY[hitNumber1]);
658 AliMUONHitForRec *hitForRec;
660 for (Int_t hitNumber = 0; hitNumber < fNTrackHits ; hitNumber++) {
661 trackParamAtHit = (AliMUONTrackParam*) fTrackParamAtHit->UncheckedAt(hitNumber);
662 hitForRec = trackParamAtHit->GetHitForRecPtr();
663 dX = hitForRec->GetNonBendingCoor() - trackParamAtHit->GetNonBendingCoor();
664 dY = hitForRec->GetBendingCoor() - trackParamAtHit->GetBendingCoor();
665 chi2 += dX * dX / hitForRec->GetNonBendingReso2() + dY * dY / hitForRec->GetBendingReso2();
674 //__________________________________________________________________________
675 Bool_t AliMUONTrack::ComputeHitWeights(TMatrixD* mcsCovariances)
677 /// Compute the weight matrices of the attached hits, in non bending and bending direction,
678 /// accounting for multiple scattering correlations and hits resolution
679 /// - Use the provided MCS covariance matrix if any (otherwise build it temporarily)
680 /// - Assume that track parameters at each hit are corrects
681 /// - Return kFALSE if computation failed
684 if (!fHitWeightsNonBending) fHitWeightsNonBending = new TMatrixD(fNTrackHits,fNTrackHits);
685 if (!fHitWeightsBending) fHitWeightsBending = new TMatrixD(fNTrackHits,fNTrackHits);
687 // Check hits (if the first one exist, assume that the other ones exit too!)
688 if (!((AliMUONTrackParam*) fTrackParamAtHit->First())->GetHitForRecPtr()) {
689 AliWarning("hit is missing");
690 fHitWeightsNonBending->ResizeTo(0,0);
691 fHitWeightsBending->ResizeTo(0,0);
695 // Compute weights matrices
696 if (!ComputeHitWeights(*fHitWeightsNonBending, *fHitWeightsBending, mcsCovariances)) return kFALSE;
702 //__________________________________________________________________________
703 Bool_t AliMUONTrack::ComputeHitWeights(TMatrixD& hitWeightsNB, TMatrixD& hitWeightsB, TMatrixD* mcsCovariances, AliMUONHitForRec* discardedHit) const
705 /// Compute the weight matrices, in non bending and bending direction,
706 /// of the other attached hits assuming the discarded one does not exist
707 /// accounting for multiple scattering correlations and hits resolution
708 /// - Use the provided MCS covariance matrix if any (otherwise build it temporarily)
709 /// - Return kFALSE if computation failed
711 // Check MCS covariance matrix and recompute it if need
712 Bool_t deleteMCSCov = kFALSE;
713 if (!mcsCovariances) {
714 mcsCovariances = new TMatrixD(fNTrackHits,fNTrackHits);
715 deleteMCSCov = kTRUE;
716 ComputeMCSCovariances(*mcsCovariances);
719 // Resize the weights matrices; alocate memory
721 hitWeightsNB.ResizeTo(fNTrackHits-1,fNTrackHits-1);
722 hitWeightsB.ResizeTo(fNTrackHits-1,fNTrackHits-1);
724 hitWeightsNB.ResizeTo(fNTrackHits,fNTrackHits);
725 hitWeightsB.ResizeTo(fNTrackHits,fNTrackHits);
729 AliMUONHitForRec *hitForRec1, *hitForRec2;
730 Int_t currentHitNumber1, currentHitNumber2;
732 // Compute the covariance matrices
733 currentHitNumber1 = 0;
734 for (Int_t hitNumber1 = 0; hitNumber1 < fNTrackHits; hitNumber1++) {
735 hitForRec1 = ((AliMUONTrackParam*) fTrackParamAtHit->UncheckedAt(hitNumber1))->GetHitForRecPtr();
737 if (hitForRec1 == discardedHit) continue;
739 // Loop over next hits
740 currentHitNumber2 = currentHitNumber1;
741 for (Int_t hitNumber2 = hitNumber1; hitNumber2 < fNTrackHits; hitNumber2++) {
742 hitForRec2 = ((AliMUONTrackParam*) fTrackParamAtHit->UncheckedAt(hitNumber2))->GetHitForRecPtr();
744 if (hitForRec2 == discardedHit) continue;
746 // Fill with MCS covariances
747 hitWeightsNB(currentHitNumber1, currentHitNumber2) = (*mcsCovariances)(hitNumber1,hitNumber2);
749 // Equal contribution from multiple scattering in non bending and bending directions
750 hitWeightsB(currentHitNumber1, currentHitNumber2) = hitWeightsNB(currentHitNumber1, currentHitNumber2);
752 // Add contribution from hit resolution to diagonal element and symmetrize the matrix
753 if (currentHitNumber1 == currentHitNumber2) {
755 // In non bending plane
756 hitWeightsNB(currentHitNumber1, currentHitNumber1) += hitForRec1->GetNonBendingReso2();
758 hitWeightsB(currentHitNumber1, currentHitNumber1) += hitForRec1->GetBendingReso2();
762 // In non bending plane
763 hitWeightsNB(currentHitNumber2, currentHitNumber1) = hitWeightsNB(currentHitNumber1, currentHitNumber2);
765 hitWeightsB(currentHitNumber2, currentHitNumber1) = hitWeightsB(currentHitNumber1, currentHitNumber2);
775 // Inversion of covariance matrices to get the weights
776 if (hitWeightsNB.Determinant() != 0 && hitWeightsB.Determinant() != 0) {
777 hitWeightsNB.Invert();
778 hitWeightsB.Invert();
780 AliWarning(" Determinant = 0");
781 hitWeightsNB.ResizeTo(0,0);
782 hitWeightsB.ResizeTo(0,0);
783 if(deleteMCSCov) delete mcsCovariances;
787 if(deleteMCSCov) delete mcsCovariances;
793 //__________________________________________________________________________
794 void AliMUONTrack::ComputeMCSCovariances(TMatrixD& mcsCovariances) const
796 /// Compute the multiple scattering covariance matrix
797 /// (assume that track parameters at each hit are corrects)
799 // Reset the size of the covariance matrix if needed
800 if (mcsCovariances.GetNrows() != fNTrackHits) mcsCovariances.ResizeTo(fNTrackHits,fNTrackHits);
803 Int_t nChambers = AliMUONConstants::NTrackingCh();
804 AliMUONTrackParam* trackParamAtHit;
805 AliMUONHitForRec *hitForRec;
806 AliMUONTrackParam extrapTrackParam;
807 Int_t currentChamber = 0, expectedChamber = 0, size = 0;
808 Double_t *mcsAngle2 = new Double_t[2*nChambers];
809 Double_t *zMCS = new Double_t[2*nChambers];
810 Int_t *indices = new Int_t[2*fNTrackHits];
812 // Compute multiple scattering dispersion angle at each chamber
813 // and save the z position where it is calculated
814 for (Int_t hitNumber = 0; hitNumber < fNTrackHits; hitNumber++) {
815 trackParamAtHit = (AliMUONTrackParam*) fTrackParamAtHit->UncheckedAt(hitNumber);
816 hitForRec = trackParamAtHit->GetHitForRecPtr();
818 // look for missing chambers if any
819 currentChamber = hitForRec->GetChamberNumber();
820 while (currentChamber > expectedChamber) {
822 // Save the z position where MCS dispersion is calculated
823 zMCS[size] = AliMUONConstants::DefaultChamberZ(expectedChamber);
825 // Do not take into account MCS in chambers prior the first hit
828 // Get track parameters at missing chamber z
829 extrapTrackParam = *trackParamAtHit;
830 AliMUONTrackExtrap::ExtrapToZ(&extrapTrackParam, zMCS[expectedChamber]);
832 // Save multiple scattering dispersion angle in missing chamber
833 mcsAngle2[size] = AliMUONTrackExtrap::GetMCSAngle2(extrapTrackParam,AliMUONConstants::ChamberThicknessInX0(),1.);
835 } else mcsAngle2[size] = 0.;
841 // Save z position where MCS dispersion is calculated
842 zMCS[size] = trackParamAtHit->GetZ();
844 // Save multiple scattering dispersion angle in current chamber
845 mcsAngle2[size] = AliMUONTrackExtrap::GetMCSAngle2(*trackParamAtHit,AliMUONConstants::ChamberThicknessInX0(),1.);
847 // Save indice in zMCS array corresponding to the current cluster
848 indices[hitNumber] = size;
850 expectedChamber = currentChamber + 1;
854 // complete array of z if last hit is on the last but one chamber
855 if (currentChamber != nChambers-1) zMCS[size++] = AliMUONConstants::DefaultChamberZ(nChambers-1);
857 // Compute the covariance matrix
858 for (Int_t hitNumber1 = 0; hitNumber1 < fNTrackHits; hitNumber1++) {
860 for (Int_t hitNumber2 = hitNumber1; hitNumber2 < fNTrackHits; hitNumber2++) {
862 // Initialization to 0 (diagonal plus upper triangular part)
863 mcsCovariances(hitNumber1,hitNumber2) = 0.;
865 // Compute contribution from multiple scattering in upstream chambers
866 for (Int_t k = 0; k < indices[hitNumber1]; k++) {
867 mcsCovariances(hitNumber1,hitNumber2) += (zMCS[indices[hitNumber1]] - zMCS[k]) * (zMCS[indices[hitNumber2]] - zMCS[k]) * mcsAngle2[k];
870 // Symetrize the matrix
871 mcsCovariances(hitNumber2,hitNumber1) = mcsCovariances(hitNumber1,hitNumber2);
882 //__________________________________________________________________________
883 Int_t AliMUONTrack::HitsInCommon(AliMUONTrack* track) const
885 /// Returns the number of hits in common between the current track ("this")
886 /// and the track pointed to by "track".
887 Int_t hitsInCommon = 0;
888 AliMUONTrackParam *trackParamAtHit1, *trackParamAtHit2;
889 // Loop over hits of first track
890 trackParamAtHit1 = (AliMUONTrackParam*) this->fTrackParamAtHit->First();
891 while (trackParamAtHit1) {
892 // Loop over hits of second track
893 trackParamAtHit2 = (AliMUONTrackParam*) track->fTrackParamAtHit->First();
894 while (trackParamAtHit2) {
895 // Increment "hitsInCommon" if both TrackParamAtHits point to the same HitForRec
896 if ((trackParamAtHit1->GetHitForRecPtr()) == (trackParamAtHit2->GetHitForRecPtr())) {
900 trackParamAtHit2 = (AliMUONTrackParam*) track->fTrackParamAtHit->After(trackParamAtHit2);
901 } // trackParamAtHit2
902 trackParamAtHit1 = (AliMUONTrackParam*) this->fTrackParamAtHit->After(trackParamAtHit1);
903 } // trackParamAtHit1
907 //__________________________________________________________________________
908 Double_t AliMUONTrack::GetNormalizedChi2() const
910 /// return the chi2 value divided by the number of degrees of freedom (or 1.e10 if ndf < 0)
912 Double_t numberOfDegFree = (2. * fNTrackHits - 5.);
913 if (numberOfDegFree > 0.) return fGlobalChi2 / numberOfDegFree;
917 //__________________________________________________________________________
918 Bool_t* AliMUONTrack::CompatibleTrack(AliMUONTrack * track, Double_t sigma2Cut) const
920 /// Return kTRUE/kFALSE for each chamber if hit is compatible or not
921 TClonesArray *hitArray, *thisHitArray;
922 AliMUONHitForRec *hit, *thisHit;
925 Float_t deltaZMax = 1.; // 1 cm
927 Bool_t *nCompHit = new Bool_t[AliMUONConstants::NTrackingCh()];
929 for ( Int_t ch = 0; ch < AliMUONConstants::NTrackingCh(); ch++) {
930 nCompHit[ch] = kFALSE;
933 thisHitArray = this->GetHitForRecAtHit();
935 hitArray = track->GetHitForRecAtHit();
937 for (Int_t iHthis = 0; iHthis < thisHitArray->GetEntriesFast(); iHthis++) {
938 thisHit = (AliMUONHitForRec*) thisHitArray->At(iHthis);
939 chamberNumber = thisHit->GetChamberNumber();
940 if (chamberNumber < 0 || chamberNumber > AliMUONConstants::NTrackingCh()) continue;
941 nCompHit[chamberNumber] = kFALSE;
942 for (Int_t iH = 0; iH < hitArray->GetEntriesFast(); iH++) {
943 hit = (AliMUONHitForRec*) hitArray->At(iH);
944 deltaZ = TMath::Abs(thisHit->GetZ() - hit->GetZ());
945 chi2 = thisHit->NormalizedChi2WithHitForRec(hit,sigma2Cut); // set cut to 4 sigmas
946 if (chi2 < 3. && deltaZ < deltaZMax) {
947 nCompHit[chamberNumber] = kTRUE;
956 //__________________________________________________________________________
957 void AliMUONTrack::RecursiveDump(void) const
959 /// Recursive dump of AliMUONTrack, i.e. with dump of TrackParamAtHit's and attached HitForRec's
960 AliMUONTrackParam *trackParamAtHit;
961 AliMUONHitForRec *hitForRec;
962 cout << "Recursive dump of Track: " << this << endl;
965 for (Int_t trackHitIndex = 0; trackHitIndex < fNTrackHits; trackHitIndex++) {
966 trackParamAtHit = (AliMUONTrackParam*) ((*fTrackParamAtHit)[trackHitIndex]);
968 cout << "TrackParamAtHit: " << trackParamAtHit << " (index: " << trackHitIndex << ")" << endl;
969 trackParamAtHit->Dump();
970 hitForRec = trackParamAtHit->GetHitForRecPtr();
972 cout << "HitForRec: " << hitForRec << endl;
978 //_____________________________________________-
979 void AliMUONTrack::Print(Option_t*) const
981 /// Printing Track information
983 cout << "<AliMUONTrack> No.Clusters=" << setw(2) << GetNTrackHits() <<
984 ", Match2Trig=" << setw(1) << GetMatchTrigger() <<
985 ", LoTrgNum=" << setw(3) << GetLoTrgNum() <<
986 ", Chi2-tracking-trigger=" << setw(8) << setprecision(5) << GetChi2MatchTrigger();
987 cout << Form(" HitTriggerPattern %x",fHitsPatternInTrigCh) << endl;
988 GetTrackParamAtHit()->First()->Print("FULL");
991 //__________________________________________________________________________
992 void AliMUONTrack::SetLocalTrigger(Int_t loCirc, Int_t loStripX, Int_t loStripY, Int_t loDev, Int_t loLpt, Int_t loHpt)
994 /// pack the local trigger information and store
996 if (loCirc < 0 || loCirc > 233) return;
999 fLocalTrigger += loCirc;
1000 fLocalTrigger += loStripX << 8;
1001 fLocalTrigger += loStripY << 13;
1002 fLocalTrigger += loDev << 17;
1003 fLocalTrigger += loLpt << 22;
1004 fLocalTrigger += loHpt << 24;