Coverity fix for uninitialized variables and check for returned null value
[u/mrichter/AliRoot.git] / MUON / AliMUONTrack.cxx
CommitLineData
a9e2aefa 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
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 **************************************************************************/
15
88cb7938 16/* $Id$ */
a9e2aefa 17
3d1463c8 18//-----------------------------------------------------------------------------
19// Class AliMUONTrack
20//-------------------
21// Reconstructed track in ALICE dimuon spectrometer
22//-----------------------------------------------------------------------------
a9e2aefa 23
63ed9c6b 24#include "AliMUONTrack.h"
34f1bfa0 25
9bf6860b 26#include "AliMUONReconstructor.h"
103e6575 27#include "AliMUONVCluster.h"
28#include "AliMUONVClusterStore.h"
7332f213 29#include "AliMUONObjectPair.h"
7332f213 30#include "AliMUONTrackExtrap.h"
933daf4c 31#include "AliMUONConstants.h"
32#include "AliMUONTrackParam.h"
a9e2aefa 33
63ed9c6b 34#include "AliLog.h"
35
63ed9c6b 36#include <TMath.h>
63ed9c6b 37
ea94c18b 38#include <Riostream.h>
39
cdd730d0 40using std::setw;
41using std::endl;
42using std::cout;
43using std::streamsize;
44using std::setprecision;
7945aae7 45/// \cond CLASSIMP
63ed9c6b 46ClassImp(AliMUONTrack) // Class implementation in ROOT context
7945aae7 47/// \endcond
63ed9c6b 48
4ea3f013 49
50const Double_t AliMUONTrack::fgkMaxChi2 = 1.e10; ///< maximum chi2 above which the track can be considered as abnormal
51
52
63ed9c6b 53//__________________________________________________________________________
30178c30 54AliMUONTrack::AliMUONTrack()
54d7ba50 55 : TObject(),
7332f213 56 fTrackParamAtCluster(0x0),
208f139e 57 fFitWithVertex(kFALSE),
96ebe67e 58 fVertexErrXY2(),
ea94c18b 59 fFitWithMCS(kFALSE),
96ebe67e 60 fClusterWeightsNonBending(0x0),
61 fClusterWeightsBending(0x0),
ea94c18b 62 fGlobalChi2(-1.),
63 fImproved(kFALSE),
7771752e 64 fMatchTrigger(-1),
54d7ba50 65 fChi2MatchTrigger(0.),
2e2d0c44 66 fTrackID(-1),
96ebe67e 67 fTrackParamAtVertex(0x0),
423b32ca 68 fHitsPatternInTrigCh(0),
0a2dcc83 69 fHitsPatternInTrigChTrk(0),
5c15a68b 70 fLocalTrigger(0),
71 fConnected(kFALSE)
d837040f 72{
2457f726 73 /// Default constructor
96ebe67e 74 fVertexErrXY2[0] = 0.;
75 fVertexErrXY2[1] = 0.;
d837040f 76}
77
78 //__________________________________________________________________________
a0dc65b4 79AliMUONTrack::AliMUONTrack(AliMUONObjectPair *segment, Double_t bendingVertexDispersion)
54d7ba50 80 : TObject(),
ce8cd162 81 fTrackParamAtCluster(new TObjArray(20)),
208f139e 82 fFitWithVertex(kFALSE),
96ebe67e 83 fVertexErrXY2(),
ea94c18b 84 fFitWithMCS(kFALSE),
96ebe67e 85 fClusterWeightsNonBending(0x0),
86 fClusterWeightsBending(0x0),
ea94c18b 87 fGlobalChi2(0.),
88 fImproved(kFALSE),
7771752e 89 fMatchTrigger(-1),
54d7ba50 90 fChi2MatchTrigger(0.),
2e2d0c44 91 fTrackID(-1),
96ebe67e 92 fTrackParamAtVertex(0x0),
423b32ca 93 fHitsPatternInTrigCh(0),
0a2dcc83 94 fHitsPatternInTrigChTrk(0),
5c15a68b 95 fLocalTrigger(0),
96 fConnected(kFALSE)
a9e2aefa 97{
96ebe67e 98 /// Constructor from two clusters
de2cd600 99
ce8cd162 100 fTrackParamAtCluster->SetOwner(kTRUE);
101
96ebe67e 102 fVertexErrXY2[0] = 0.;
103 fVertexErrXY2[1] = 0.;
208f139e 104
96ebe67e 105 // Pointers to clusters from the segment
9bf6860b 106 AliMUONVCluster* firstCluster = (AliMUONVCluster*) segment->First();
107 AliMUONVCluster* lastCluster = (AliMUONVCluster*) segment->Second();
208f139e 108
ea94c18b 109 // Compute track parameters
96ebe67e 110 Double_t z1 = firstCluster->GetZ();
111 Double_t z2 = lastCluster->GetZ();
112 Double_t dZ = z1 - z2;
208f139e 113 // Non bending plane
96ebe67e 114 Double_t nonBendingCoor1 = firstCluster->GetX();
115 Double_t nonBendingCoor2 = lastCluster->GetX();
ea94c18b 116 Double_t nonBendingSlope = (nonBendingCoor1 - nonBendingCoor2) / dZ;
208f139e 117 // Bending plane
96ebe67e 118 Double_t bendingCoor1 = firstCluster->GetY();
119 Double_t bendingCoor2 = lastCluster->GetY();
ea94c18b 120 Double_t bendingSlope = (bendingCoor1 - bendingCoor2) / dZ;
208f139e 121 // Inverse bending momentum
96ebe67e 122 Double_t bendingImpact = bendingCoor1 - z1 * bendingSlope;
208f139e 123 Double_t inverseBendingMomentum = 1. / AliMUONTrackExtrap::GetBendingMomentumFromImpactParam(bendingImpact);
ea94c18b 124
9bf6860b 125 // Set track parameters at first cluster
96ebe67e 126 AliMUONTrackParam trackParamAtFirstCluster;
127 trackParamAtFirstCluster.SetZ(z1);
128 trackParamAtFirstCluster.SetNonBendingCoor(nonBendingCoor1);
129 trackParamAtFirstCluster.SetNonBendingSlope(nonBendingSlope);
130 trackParamAtFirstCluster.SetBendingCoor(bendingCoor1);
131 trackParamAtFirstCluster.SetBendingSlope(bendingSlope);
132 trackParamAtFirstCluster.SetInverseBendingMomentum(inverseBendingMomentum);
208f139e 133
9bf6860b 134 // Set track parameters at last cluster
96ebe67e 135 AliMUONTrackParam trackParamAtLastCluster;
136 trackParamAtLastCluster.SetZ(z2);
137 trackParamAtLastCluster.SetNonBendingCoor(nonBendingCoor2);
138 trackParamAtLastCluster.SetNonBendingSlope(nonBendingSlope);
139 trackParamAtLastCluster.SetBendingCoor(bendingCoor2);
140 trackParamAtLastCluster.SetBendingSlope(bendingSlope);
141 trackParamAtLastCluster.SetInverseBendingMomentum(inverseBendingMomentum);
ea94c18b 142
9bf6860b 143 // Compute and set track parameters covariances at first cluster
144 TMatrixD paramCov(5,5);
145 paramCov.Zero();
208f139e 146 // Non bending plane
9bf6860b 147 paramCov(0,0) = firstCluster->GetErrX2();
148 paramCov(0,1) = firstCluster->GetErrX2() / dZ;
149 paramCov(1,0) = paramCov(0,1);
150 paramCov(1,1) = ( firstCluster->GetErrX2() + lastCluster->GetErrX2() ) / dZ / dZ;
208f139e 151 // Bending plane
9bf6860b 152 paramCov(2,2) = firstCluster->GetErrY2();
153 paramCov(2,3) = firstCluster->GetErrY2() / dZ;
154 paramCov(3,2) = paramCov(2,3);
155 paramCov(3,3) = ( firstCluster->GetErrY2() + lastCluster->GetErrY2() ) / dZ / dZ;
156 // Inverse bending momentum (vertex resolution + bending slope resolution + 10% error on dipole parameters+field)
89c8d66d 157 if (AliMUONTrackExtrap::IsFieldON()) {
158 paramCov(4,4) = ( ( bendingVertexDispersion*bendingVertexDispersion +
159 (z1 * z1 * lastCluster->GetErrY2() + z2 * z2 * firstCluster->GetErrY2()) / dZ / dZ) /
160 bendingImpact / bendingImpact + 0.1 * 0.1) * inverseBendingMomentum * inverseBendingMomentum ;
161 paramCov(2,4) = - z2 * firstCluster->GetErrY2() * inverseBendingMomentum / bendingImpact / dZ;
162 paramCov(4,2) = paramCov(2,4);
163 paramCov(3,4) = - (z1 * lastCluster->GetErrY2() + z2 * firstCluster->GetErrY2()) * inverseBendingMomentum / bendingImpact / dZ / dZ;
164 paramCov(4,3) = paramCov(3,4);
165 } else paramCov(4,4) = inverseBendingMomentum*inverseBendingMomentum;
9bf6860b 166 trackParamAtFirstCluster.SetCovariances(paramCov);
ea94c18b 167
9bf6860b 168 // Compute and set track parameters covariances at last cluster
89c8d66d 169 // Non bending plane
170 paramCov(0,0) = lastCluster->GetErrX2();
171 paramCov(0,1) = - lastCluster->GetErrX2() / dZ;
172 paramCov(1,0) = paramCov(0,1);
173 // Bending plane
174 paramCov(2,2) = lastCluster->GetErrY2();
175 paramCov(2,3) = - lastCluster->GetErrY2() / dZ;
176 paramCov(3,2) = paramCov(2,3);
177 // Inverse bending momentum (vertex resolution + bending slope resolution + 10% error on dipole parameters+field)
178 if (AliMUONTrackExtrap::IsFieldON()) {
179 paramCov(2,4) = z1 * lastCluster->GetErrY2() * inverseBendingMomentum / bendingImpact / dZ;
180 paramCov(4,2) = paramCov(2,4);
181 }
9bf6860b 182 trackParamAtLastCluster.SetCovariances(paramCov);
ea94c18b 183
96ebe67e 184 // Add track parameters at clusters
185 AddTrackParamAtCluster(trackParamAtFirstCluster,*firstCluster);
186 AddTrackParamAtCluster(trackParamAtLastCluster,*lastCluster);
de2cd600 187
a9e2aefa 188}
189
61fed964 190//__________________________________________________________________________
96ebe67e 191AliMUONTrack::AliMUONTrack(const AliMUONTrack& track)
ea94c18b 192 : TObject(track),
7332f213 193 fTrackParamAtCluster(0x0),
ea94c18b 194 fFitWithVertex(track.fFitWithVertex),
96ebe67e 195 fVertexErrXY2(),
ea94c18b 196 fFitWithMCS(track.fFitWithMCS),
96ebe67e 197 fClusterWeightsNonBending(0x0),
198 fClusterWeightsBending(0x0),
ea94c18b 199 fGlobalChi2(track.fGlobalChi2),
200 fImproved(track.fImproved),
201 fMatchTrigger(track.fMatchTrigger),
ea94c18b 202 fChi2MatchTrigger(track.fChi2MatchTrigger),
203 fTrackID(track.fTrackID),
96ebe67e 204 fTrackParamAtVertex(0x0),
ea94c18b 205 fHitsPatternInTrigCh(track.fHitsPatternInTrigCh),
0a2dcc83 206 fHitsPatternInTrigChTrk(track.fHitsPatternInTrigChTrk),
5c15a68b 207 fLocalTrigger(track.fLocalTrigger),
208 fConnected(track.fConnected)
a9e2aefa 209{
2457f726 210 ///copy constructor
de2cd600 211
ce8cd162 212 // necessary to make a copy of the objects and not only the pointers in TObjArray.
7332f213 213 if (track.fTrackParamAtCluster) {
ce8cd162 214 fTrackParamAtCluster = new TObjArray(track.fTrackParamAtCluster->GetSize());
215 fTrackParamAtCluster->SetOwner(kTRUE);
216 for (Int_t i = 0; i < track.GetNClusters(); i++)
217 fTrackParamAtCluster->AddLast(new AliMUONTrackParam(*static_cast<AliMUONTrackParam*>(track.fTrackParamAtCluster->UncheckedAt(i))));
208f139e 218 }
219
96ebe67e 220 // copy vertex resolution square used during the tracking procedure
221 fVertexErrXY2[0] = track.fVertexErrXY2[0];
222 fVertexErrXY2[1] = track.fVertexErrXY2[1];
208f139e 223
96ebe67e 224 // copy cluster weights matrices if any
225 if (track.fClusterWeightsNonBending) fClusterWeightsNonBending = new TMatrixD(*(track.fClusterWeightsNonBending));
226 if (track.fClusterWeightsBending) fClusterWeightsBending = new TMatrixD(*(track.fClusterWeightsBending));
ea94c18b 227
96ebe67e 228 // copy track parameters at vertex if any
229 if (track.fTrackParamAtVertex) fTrackParamAtVertex = new AliMUONTrackParam(*(track.fTrackParamAtVertex));
208f139e 230
a9e2aefa 231}
232
956019b6 233 //__________________________________________________________________________
ea94c18b 234AliMUONTrack & AliMUONTrack::operator=(const AliMUONTrack& track)
a9e2aefa 235{
2457f726 236 /// Asignment operator
30178c30 237 // check assignement to self
ea94c18b 238 if (this == &track)
a9e2aefa 239 return *this;
61adb9bd 240
30178c30 241 // base class assignement
ea94c18b 242 TObject::operator=(track);
7332f213 243
244 // clear memory
245 Clear();
246
ce8cd162 247 // necessary to make a copy of the objects and not only the pointers in TObjArray
7332f213 248 if (track.fTrackParamAtCluster) {
ce8cd162 249 fTrackParamAtCluster = new TObjArray(track.fTrackParamAtCluster->GetSize());
250 fTrackParamAtCluster->SetOwner(kTRUE);
251 for (Int_t i = 0; i < track.GetNClusters(); i++)
252 fTrackParamAtCluster->AddLast(new AliMUONTrackParam(*static_cast<AliMUONTrackParam*>(track.fTrackParamAtCluster->UncheckedAt(i))));
208f139e 253 }
254
96ebe67e 255 // copy cluster weights matrix if any
256 if (track.fClusterWeightsNonBending) {
257 if (fClusterWeightsNonBending) {
258 fClusterWeightsNonBending->ResizeTo(*(track.fClusterWeightsNonBending));
259 *fClusterWeightsNonBending = *(track.fClusterWeightsNonBending);
260 } else fClusterWeightsNonBending = new TMatrixD(*(track.fClusterWeightsNonBending));
208f139e 261 }
262
96ebe67e 263 // copy cluster weights matrix if any
264 if (track.fClusterWeightsBending) {
265 if (fClusterWeightsBending) {
266 fClusterWeightsBending->ResizeTo(*(track.fClusterWeightsBending));
267 *fClusterWeightsBending = *(track.fClusterWeightsBending);
268 } else fClusterWeightsBending = new TMatrixD(*(track.fClusterWeightsBending));
ea94c18b 269 }
270
96ebe67e 271 // copy track parameters at vertex if any
272 if (track.fTrackParamAtVertex) {
273 if (fTrackParamAtVertex) *fTrackParamAtVertex = *(track.fTrackParamAtVertex);
274 else fTrackParamAtVertex = new AliMUONTrackParam(*(track.fTrackParamAtVertex));
ea94c18b 275 }
de2cd600 276
ea94c18b 277 fFitWithVertex = track.fFitWithVertex;
96ebe67e 278 fVertexErrXY2[0] = track.fVertexErrXY2[0];
279 fVertexErrXY2[1] = track.fVertexErrXY2[1];
ea94c18b 280 fFitWithMCS = track.fFitWithMCS;
281 fGlobalChi2 = track.fGlobalChi2;
282 fImproved = track.fImproved;
283 fMatchTrigger = track.fMatchTrigger;
ea94c18b 284 fChi2MatchTrigger = track.fChi2MatchTrigger;
285 fTrackID = track.fTrackID;
286 fHitsPatternInTrigCh = track.fHitsPatternInTrigCh;
0a2dcc83 287 fHitsPatternInTrigChTrk = track.fHitsPatternInTrigChTrk;
ea94c18b 288 fLocalTrigger = track.fLocalTrigger;
5c15a68b 289 fConnected = track.fConnected;
30178c30 290
61adb9bd 291 return *this;
a9e2aefa 292}
293
04b5ea16 294 //__________________________________________________________________________
ea94c18b 295AliMUONTrack::~AliMUONTrack()
296{
297 /// Destructor
96ebe67e 298 delete fTrackParamAtCluster;
299 delete fClusterWeightsNonBending;
300 delete fClusterWeightsBending;
301 delete fTrackParamAtVertex;
ea94c18b 302}
303
304 //__________________________________________________________________________
ce8cd162 305void AliMUONTrack::Clear(Option_t* /*opt*/)
ea94c18b 306{
307 /// Clear arrays
ce8cd162 308 delete fTrackParamAtCluster; fTrackParamAtCluster = 0x0;
96ebe67e 309 delete fClusterWeightsNonBending; fClusterWeightsNonBending = 0x0;
310 delete fClusterWeightsBending; fClusterWeightsBending = 0x0;
311 delete fTrackParamAtVertex; fTrackParamAtVertex = 0x0;
ea94c18b 312}
313
314 //__________________________________________________________________________
b1fea02e 315void AliMUONTrack::Reset()
316{
317 /// Reset to default values
318 SetUniqueID(0);
319 fFitWithVertex = kFALSE;
320 fVertexErrXY2[0] = 0.;
321 fVertexErrXY2[1] = 0.;
322 fFitWithMCS = kFALSE;
323 fGlobalChi2 = -1.;
324 fImproved = kFALSE;
325 fMatchTrigger = -1;
b1fea02e 326 fChi2MatchTrigger = 0.;
2e2d0c44 327 fTrackID = -1;
b1fea02e 328 fHitsPatternInTrigCh = 0;
0a2dcc83 329 fHitsPatternInTrigChTrk = 0;
b1fea02e 330 fLocalTrigger = 0;
5c15a68b 331 fConnected = kFALSE;
b1fea02e 332 delete fTrackParamAtCluster; fTrackParamAtCluster = 0x0;
333 delete fClusterWeightsNonBending; fClusterWeightsNonBending = 0x0;
334 delete fClusterWeightsBending; fClusterWeightsBending = 0x0;
335 delete fTrackParamAtVertex; fTrackParamAtVertex = 0x0;
336}
337
338 //__________________________________________________________________________
ce8cd162 339TObjArray* AliMUONTrack::GetTrackParamAtCluster() const
7332f213 340{
341 /// return array of track parameters at cluster (create it if needed)
ce8cd162 342 if (!fTrackParamAtCluster) {
343 fTrackParamAtCluster = new TObjArray(20);
344 fTrackParamAtCluster->SetOwner(kTRUE);
345 }
7332f213 346 return fTrackParamAtCluster;
347}
348
349 //__________________________________________________________________________
96ebe67e 350void AliMUONTrack::AddTrackParamAtCluster(const AliMUONTrackParam &trackParam, AliMUONVCluster &cluster, Bool_t copy)
8429a5e4 351{
96ebe67e 352 /// Copy given track parameters into a new TrackParamAtCluster
353 /// Link parameters with the associated cluster
354 /// If copy=kTRUE: the cluster is copied then passed the trackParam which become its owner
355 /// otherwise: make sure to do not delete the cluster until it is used by the track
356
357 // check chamber ID of the associated cluster
358 if (cluster.GetChamberId() < 0 || cluster.GetChamberId() > AliMUONConstants::NTrackingCh()) {
359 AliError(Form("Chamber ID of the associated cluster is not valid (ChamberId=%d)",cluster.GetChamberId()));
ea94c18b 360 return;
361 }
362
96ebe67e 363 // check whether track parameters are given at the correct cluster z position
820b4d9e 364 if (TMath::Abs(cluster.GetZ() - trackParam.GetZ())>1.e-5) { // AU
96ebe67e 365 AliError("track parameters are given at a different z position than the one of the associated cluster");
ea94c18b 366 return;
367 }
368
96ebe67e 369 // add parameters to the array of track parameters
ce8cd162 370 if (!fTrackParamAtCluster) {
371 fTrackParamAtCluster = new TObjArray(20);
372 fTrackParamAtCluster->SetOwner(kTRUE);
373 }
374 AliMUONTrackParam* trackParamAtCluster = new AliMUONTrackParam(trackParam);
375 fTrackParamAtCluster->AddLast(trackParamAtCluster);
96ebe67e 376
377 // link parameters with the associated cluster or its copy
378 if (copy) {
1467f4ba 379 AliMUONVCluster *clusterCopy = static_cast<AliMUONVCluster*>(cluster.Clone());
96ebe67e 380 trackParamAtCluster->SetClusterPtr(clusterCopy, kTRUE);
381 } else trackParamAtCluster->SetClusterPtr(&cluster);
7332f213 382
383 // sort the array of track parameters
384 fTrackParamAtCluster->Sort();
ea94c18b 385}
386
387 //__________________________________________________________________________
96ebe67e 388void AliMUONTrack::RemoveTrackParamAtCluster(AliMUONTrackParam *trackParam)
956019b6 389{
ce8cd162 390 /// Remove trackParam from the array of TrackParamAtCluster and delete it since the array is owner
391
392 if (fTrackParamAtCluster) {
393
394 AliMUONTrackParam* trackParamAtCluster = static_cast<AliMUONTrackParam*>(fTrackParamAtCluster->Remove(trackParam));
395
396 if (trackParamAtCluster) {
397
398 // clean memory
399 delete trackParamAtCluster;
400
401 // remove hole
402 fTrackParamAtCluster->Compress();
403
404 } else AliWarning("object to remove does not exist in array fTrackParamAtCluster");
405
406 } else AliWarning("array fTrackParamAtCluster does not exist");
de2cd600 407
04b5ea16 408}
409
ea94c18b 410 //__________________________________________________________________________
4ea3f013 411Bool_t AliMUONTrack::UpdateTrackParamAtCluster()
d2b1e7bb 412{
96ebe67e 413 /// Update track parameters at each attached cluster
4ea3f013 414 /// Return kFALSE in case of failure (i.e. extrapolation problem)
ea94c18b 415
ce8cd162 416 Int_t nClusters = GetNClusters();
417 if (nClusters == 0) {
96ebe67e 418 AliWarning("no cluster attached to the track");
4ea3f013 419 return kFALSE;
ea94c18b 420 }
421
4ea3f013 422 Bool_t extrapStatus = kTRUE;
ce8cd162 423 AliMUONTrackParam* startingTrackParam = static_cast<AliMUONTrackParam*>(fTrackParamAtCluster->UncheckedAt(0));
424
425 for (Int_t i = 1; i < nClusters; i++) {
426 AliMUONTrackParam* trackParamAtCluster = static_cast<AliMUONTrackParam*>(fTrackParamAtCluster->UncheckedAt(i));
ea94c18b 427
428 // reset track parameters and their covariances
96ebe67e 429 trackParamAtCluster->SetParameters(startingTrackParam->GetParameters());
430 trackParamAtCluster->SetZ(startingTrackParam->GetZ());
ea94c18b 431
432 // extrapolation to the given z
4ea3f013 433 if (!AliMUONTrackExtrap::ExtrapToZ(trackParamAtCluster, trackParamAtCluster->GetClusterPtr()->GetZ())) extrapStatus = kFALSE;
ea94c18b 434
435 // prepare next step
96ebe67e 436 startingTrackParam = trackParamAtCluster;
ea94c18b 437 }
438
4ea3f013 439 // set global chi2 to max value in case of problem during track extrapolation
440 if (!extrapStatus) SetGlobalChi2(2.*MaxChi2());
441 return extrapStatus;
442
d2b1e7bb 443}
444
956019b6 445 //__________________________________________________________________________
4ea3f013 446Bool_t AliMUONTrack::UpdateCovTrackParamAtCluster()
ea94c18b 447{
96ebe67e 448 /// Update track parameters and their covariances at each attached cluster
449 /// Include effects of multiple scattering in chambers
4ea3f013 450 /// Return kFALSE in case of failure (i.e. extrapolation problem)
ea94c18b 451
ce8cd162 452 Int_t nClusters = GetNClusters();
453 if (nClusters == 0) {
96ebe67e 454 AliWarning("no cluster attached to the track");
4ea3f013 455 return kFALSE;
ea94c18b 456 }
457
4ea3f013 458 Bool_t extrapStatus = kTRUE;
ce8cd162 459 AliMUONTrackParam* startingTrackParam = static_cast<AliMUONTrackParam*>(fTrackParamAtCluster->UncheckedAt(0));
96ebe67e 460 Int_t expectedChamber = startingTrackParam->GetClusterPtr()->GetChamberId() + 1;
461 Int_t currentChamber;
ce8cd162 462
463 for (Int_t i = 1; i < nClusters; i++) {
464 AliMUONTrackParam* trackParamAtCluster = static_cast<AliMUONTrackParam*>(fTrackParamAtCluster->UncheckedAt(i));
ea94c18b 465
466 // reset track parameters and their covariances
96ebe67e 467 trackParamAtCluster->SetParameters(startingTrackParam->GetParameters());
468 trackParamAtCluster->SetZ(startingTrackParam->GetZ());
469 trackParamAtCluster->SetCovariances(startingTrackParam->GetCovariances());
ea94c18b 470
96ebe67e 471 // add MCS effect
4663da9f 472 AliMUONTrackExtrap::AddMCSEffect(trackParamAtCluster,AliMUONConstants::ChamberThicknessInX0(expectedChamber-1),-1.);
96ebe67e 473
474 // add MCS in missing chambers if any
475 currentChamber = trackParamAtCluster->GetClusterPtr()->GetChamberId();
476 while (currentChamber > expectedChamber) {
477 // extrapolation to the missing chamber
4ea3f013 478 if (!AliMUONTrackExtrap::ExtrapToZCov(trackParamAtCluster, AliMUONConstants::DefaultChamberZ(expectedChamber))) extrapStatus = kFALSE;
96ebe67e 479 // add MCS effect
4663da9f 480 AliMUONTrackExtrap::AddMCSEffect(trackParamAtCluster,AliMUONConstants::ChamberThicknessInX0(expectedChamber),-1.);
96ebe67e 481 expectedChamber++;
482 }
483
484 // extrapolation to the z of the current cluster
4ea3f013 485 if (!AliMUONTrackExtrap::ExtrapToZCov(trackParamAtCluster, trackParamAtCluster->GetClusterPtr()->GetZ())) extrapStatus = kFALSE;
ea94c18b 486
487 // prepare next step
96ebe67e 488 expectedChamber = currentChamber + 1;
489 startingTrackParam = trackParamAtCluster;
ea94c18b 490 }
96ebe67e 491
4ea3f013 492 // set global chi2 to max value in case of problem during track extrapolation
493 if (!extrapStatus) SetGlobalChi2(2.*MaxChi2());
494 return extrapStatus;
495
208f139e 496}
497
ea94c18b 498 //__________________________________________________________________________
f202486b 499Bool_t AliMUONTrack::IsValid(UInt_t requestedStationMask, Bool_t request2ChInSameSt45)
7332f213 500{
3d7acbb7 501 /// check the validity of the current track:
502 /// at least one cluster per requested station
503 /// and at least 2 chambers in stations 4 & 5 that contain cluster(s)
f202486b 504 /// + if request2ChInSameSt45 = kTRUE: 2 chambers hit in the same station (4 or 5)
7332f213 505
506 Int_t nClusters = GetNClusters();
507 AliMUONTrackParam *trackParam;
f202486b 508 Int_t currentCh, currentSt, previousCh = -1, nChHitInSt4 = 0, nChHitInSt5 = 0;
ad3c6eda 509 UInt_t presentStationMask(0);
7332f213 510
3d7acbb7 511 // first loop over clusters
ad3c6eda 512 for (Int_t i = 0; i < nClusters; i++) {
7332f213 513 trackParam = (AliMUONTrackParam*) fTrackParamAtCluster->UncheckedAt(i);
3d7acbb7 514
515 currentCh = trackParam->GetClusterPtr()->GetChamberId();
516 currentSt = currentCh/2;
517
518 // build present station mask
519 presentStationMask |= ( 1 << currentSt );
520
f202486b 521 // count the number of chambers hit in station 4 that contain cluster(s)
522 if (currentSt == 3 && currentCh != previousCh) {
523 nChHitInSt4++;
524 previousCh = currentCh;
525 }
526
527 // count the number of chambers hit in station 5 that contain cluster(s)
528 if (currentSt == 4 && currentCh != previousCh) {
529 nChHitInSt5++;
3d7acbb7 530 previousCh = currentCh;
531 }
532
7332f213 533 }
534
f202486b 535 // at least one cluster per requested station
536 if ((requestedStationMask & presentStationMask) != requestedStationMask) return kFALSE;
537
538 // 2 chambers hit in the same station (4 or 5)
539 if (request2ChInSameSt45) return (nChHitInSt4 == 2 || nChHitInSt5 == 2);
540 // or 2 chambers hit in station 4 & 5 together
541 else return (nChHitInSt4+nChHitInSt5 >= 2);
542
7332f213 543}
544
545 //__________________________________________________________________________
a0dc65b4 546void AliMUONTrack::TagRemovableClusters(UInt_t requestedStationMask) {
7332f213 547 /// Identify clusters that can be removed from the track,
3d7acbb7 548 /// with the only requirements to have at least 1 cluster per requested station
549 /// and at least 2 chambers over 4 in stations 4 & 5 that contain cluster(s)
7332f213 550
551 Int_t nClusters = GetNClusters();
552 AliMUONTrackParam *trackParam, *nextTrackParam;
3d7acbb7 553 Int_t currentCh, nextCh, currentSt, nextSt, previousCh = -1, nChHitInSt45 = 0;
7332f213 554
3d7acbb7 555 // first loop over clusters
7332f213 556 for (Int_t i = 0; i < nClusters; i++) {
557 trackParam = (AliMUONTrackParam*) fTrackParamAtCluster->UncheckedAt(i);
3d7acbb7 558
559 currentCh = trackParam->GetClusterPtr()->GetChamberId();
560 currentSt = currentCh/2;
561
562 // reset flags to kFALSE for all clusters in required station
563 if ((1 << currentSt) & requestedStationMask) trackParam->SetRemovable(kFALSE);
9bf6860b 564 else trackParam->SetRemovable(kTRUE);
3d7acbb7 565
566 // count the number of chambers in station 4 & 5 that contain cluster(s)
567 if (currentCh > 5 && currentCh != previousCh) {
568 nChHitInSt45++;
569 previousCh = currentCh;
570 }
571
7332f213 572 }
573
3d7acbb7 574 // second loop over clusters
7332f213 575 for (Int_t i = 0; i < nClusters; i++) {
576 trackParam = (AliMUONTrackParam*) fTrackParamAtCluster->UncheckedAt(i);
577
578 currentCh = trackParam->GetClusterPtr()->GetChamberId();
3d7acbb7 579 currentSt = currentCh/2;
7332f213 580
3d7acbb7 581 // make sure they are more than 2 clusters in 2 different chambers of stations 4 & 5
582 // but 2 clusters in he same chamber will still be flagged as removable
583 if (nChHitInSt45 < 3 && currentSt > 2) {
7332f213 584
3d7acbb7 585 if (i == nClusters-1) {
586
587 trackParam->SetRemovable(kFALSE);
7332f213 588
3d7acbb7 589 } else {
590
591 nextTrackParam = (AliMUONTrackParam*) fTrackParamAtCluster->UncheckedAt(i+1);
592 nextCh = nextTrackParam->GetClusterPtr()->GetChamberId();
593
594 // set clusters in the same chamber as being removable
595 if (nextCh == currentCh) {
596 trackParam->SetRemovable(kTRUE);
597 nextTrackParam->SetRemovable(kTRUE);
598 i++; // skip cluster already checked
599 } else {
600 trackParam->SetRemovable(kFALSE);
601 }
602
603 }
7332f213 604
3d7acbb7 605 } else {
606
607 // skip clusters already flag as removable
608 if (trackParam->IsRemovable()) continue;
609
610 // loop over next track parameters
611 for (Int_t j = i+1; j < nClusters; j++) {
612 nextTrackParam = (AliMUONTrackParam*) fTrackParamAtCluster->UncheckedAt(j);
613
614 nextCh = nextTrackParam->GetClusterPtr()->GetChamberId();
615 nextSt = nextCh/2;
616
617 // set clusters in the same station as being removable
618 if (nextSt == currentSt) {
619 trackParam->SetRemovable(kTRUE);
620 nextTrackParam->SetRemovable(kTRUE);
621 i++; // skip cluster already checked
622 }
623
624 }
7332f213 625
626 }
3d7acbb7 627
7332f213 628 }
629
630}
631
632 //__________________________________________________________________________
ea94c18b 633Bool_t AliMUONTrack::ComputeLocalChi2(Bool_t accountForMCS)
634{
96ebe67e 635 /// Compute each cluster contribution to the chi2 of the track
ea94c18b 636 /// accounting for multiple scattering or not according to the flag
96ebe67e 637 /// - Also recompute the weight matrices of the attached clusters if accountForMCS=kTRUE
638 /// - Assume that track parameters at each cluster are corrects
ea94c18b 639 /// - Return kFALSE if computation failed
7332f213 640 AliDebug(1,"Enter ComputeLocalChi2");
641
642 if (!fTrackParamAtCluster) {
643 AliWarning("no cluster attached to this track");
644 return kFALSE;
645 }
ea94c18b 646
ea94c18b 647 if (accountForMCS) { // Compute local chi2 taking into account multiple scattering effects
648
649 // Compute MCS covariance matrix only once
96ebe67e 650 Int_t nClusters = GetNClusters();
651 TMatrixD mcsCovariances(nClusters,nClusters);
ea94c18b 652 ComputeMCSCovariances(mcsCovariances);
653
96ebe67e 654 // Make sure cluster weights are consistent with following calculations
655 if (!ComputeClusterWeights(&mcsCovariances)) {
ea94c18b 656 AliWarning("cannot take into account the multiple scattering effects");
657 return ComputeLocalChi2(kFALSE);
658 }
659
660 // Compute chi2 of the track
661 Double_t globalChi2 = ComputeGlobalChi2(kTRUE);
662 if (globalChi2 < 0.) return kFALSE;
663
96ebe67e 664 // Loop over removable clusters and compute their local chi2
ce8cd162 665 AliMUONTrackParam* trackParamAtCluster;
96ebe67e 666 AliMUONTrackParam* trackParamAtCluster1;
667 AliMUONVCluster *cluster, *discardedCluster;
668 Int_t iCluster1, iCluster2, iCurrentCluster1, iCurrentCluster2;
c812dd48 669 TMatrixD clusterWeightsNB(nClusters-1,nClusters-1);
670 TMatrixD clusterWeightsB(nClusters-1,nClusters-1);
96ebe67e 671 Double_t *dX = new Double_t[nClusters-1];
672 Double_t *dY = new Double_t[nClusters-1];
ea94c18b 673 Double_t globalChi2b;
ce8cd162 674 for (Int_t iCluster = 0; iCluster < nClusters ; iCluster++) {
675 trackParamAtCluster = static_cast<AliMUONTrackParam*>(fTrackParamAtCluster->UncheckedAt(iCluster));
ea94c18b 676
96ebe67e 677 discardedCluster = trackParamAtCluster->GetClusterPtr();
ea94c18b 678
96ebe67e 679 // Recompute cluster weights without the current cluster
c812dd48 680 if (!ComputeClusterWeights(clusterWeightsNB, clusterWeightsB, &mcsCovariances, discardedCluster)) {
ea94c18b 681 AliWarning("cannot take into account the multiple scattering effects");
96ebe67e 682 delete [] dX;
683 delete [] dY;
684 return ComputeLocalChi2(kFALSE);
ea94c18b 685 }
686
96ebe67e 687 // Compute track chi2 without the current cluster
ea94c18b 688 globalChi2b = 0.;
96ebe67e 689 iCurrentCluster1 = 0;
690 for (iCluster1 = 0; iCluster1 < nClusters ; iCluster1++) {
691 trackParamAtCluster1 = (AliMUONTrackParam*) fTrackParamAtCluster->UncheckedAt(iCluster1);
692 cluster = trackParamAtCluster1->GetClusterPtr();
ea94c18b 693
96ebe67e 694 if (cluster == discardedCluster) continue;
ea94c18b 695
696 // Compute and save residuals
96ebe67e 697 dX[iCurrentCluster1] = cluster->GetX() - trackParamAtCluster1->GetNonBendingCoor();
698 dY[iCurrentCluster1] = cluster->GetY() - trackParamAtCluster1->GetBendingCoor();
ea94c18b 699
96ebe67e 700 iCurrentCluster2 = 0;
701 for (iCluster2 = 0; iCluster2 < iCluster1; iCluster2++) {
702 cluster = ((AliMUONTrackParam*) fTrackParamAtCluster->UncheckedAt(iCluster2))->GetClusterPtr();
ea94c18b 703
96ebe67e 704 if (cluster == discardedCluster) continue;
ea94c18b 705
706 // Add contribution from covariances
c812dd48 707 globalChi2b += (clusterWeightsNB(iCurrentCluster1, iCurrentCluster2) +
708 clusterWeightsNB(iCurrentCluster2, iCurrentCluster1)) * dX[iCurrentCluster1] * dX[iCurrentCluster2] +
709 (clusterWeightsB(iCurrentCluster1, iCurrentCluster2) +
710 clusterWeightsB(iCurrentCluster2, iCurrentCluster1)) * dY[iCurrentCluster1] * dY[iCurrentCluster2];
ea94c18b 711
96ebe67e 712 iCurrentCluster2++;
ea94c18b 713 }
714
715 // Add contribution from variances
c812dd48 716 globalChi2b += clusterWeightsNB(iCurrentCluster1, iCurrentCluster1) * dX[iCurrentCluster1] * dX[iCurrentCluster1] +
717 clusterWeightsB(iCurrentCluster1, iCurrentCluster1) * dY[iCurrentCluster1] * dY[iCurrentCluster1];
ea94c18b 718
96ebe67e 719 iCurrentCluster1++;
ea94c18b 720 }
721
722 // Set local chi2
96ebe67e 723 trackParamAtCluster->SetLocalChi2(globalChi2 - globalChi2b);
ea94c18b 724 }
725
726 delete [] dX;
727 delete [] dY;
728
729 } else { // without multiple scattering effects
730
ce8cd162 731 Int_t nClusters = GetNClusters();
732 AliMUONTrackParam* trackParamAtCluster;
96ebe67e 733 AliMUONVCluster *discardedCluster;
ea94c18b 734 Double_t dX, dY;
ce8cd162 735 for (Int_t iCluster = 0; iCluster < nClusters ; iCluster++) {
736 trackParamAtCluster = static_cast<AliMUONTrackParam*>(fTrackParamAtCluster->UncheckedAt(iCluster));
ea94c18b 737
96ebe67e 738 discardedCluster = trackParamAtCluster->GetClusterPtr();
ea94c18b 739
740 // Compute residuals
96ebe67e 741 dX = discardedCluster->GetX() - trackParamAtCluster->GetNonBendingCoor();
742 dY = discardedCluster->GetY() - trackParamAtCluster->GetBendingCoor();
ea94c18b 743
744 // Set local chi2
96ebe67e 745 trackParamAtCluster->SetLocalChi2(dX * dX / discardedCluster->GetErrX2() + dY * dY / discardedCluster->GetErrY2());
ea94c18b 746 }
747
748 }
749
ea94c18b 750 return kTRUE;
751
752}
753
754 //__________________________________________________________________________
755Double_t AliMUONTrack::ComputeGlobalChi2(Bool_t accountForMCS)
756{
757 /// Compute the chi2 of the track accounting for multiple scattering or not according to the flag
96ebe67e 758 /// - Assume that track parameters at each cluster are corrects
759 /// - Assume the cluster weights matrices are corrects
4ea3f013 760 /// - Return a value of chi2 higher than the maximum allowed if computation failed
7332f213 761 AliDebug(1,"Enter ComputeGlobalChi2");
762
763 if (!fTrackParamAtCluster) {
764 AliWarning("no cluster attached to this track");
4ea3f013 765 return 2.*MaxChi2();
7332f213 766 }
ea94c18b 767
ea94c18b 768 Double_t chi2 = 0.;
769
770 if (accountForMCS) {
771
b709ac13 772 // Check the weight matrices. If weight matrices are not available compute chi2 without MCS
96ebe67e 773 if (!fClusterWeightsNonBending || !fClusterWeightsBending) {
774 AliWarning("cluster weights including multiple scattering effects are not available\n\t\t --> compute chi2 WITHOUT multiple scattering");
b709ac13 775 return ComputeGlobalChi2(kFALSE);
776 }
96ebe67e 777 Int_t nClusters = GetNClusters();
778 if (fClusterWeightsNonBending->GetNrows() != nClusters || fClusterWeightsBending->GetNcols() != nClusters) {
779 AliWarning("cluster weights including multiple scattering effects are not available\n\t\t --> compute chi2 WITHOUT multiple scattering");
ea94c18b 780 return ComputeGlobalChi2(kFALSE);
781 }
782
783 // Compute chi2
96ebe67e 784 AliMUONVCluster *cluster;
785 Double_t *dX = new Double_t[nClusters];
786 Double_t *dY = new Double_t[nClusters];
787 AliMUONTrackParam* trackParamAtCluster;
788 for (Int_t iCluster1 = 0; iCluster1 < nClusters; iCluster1++) {
789 trackParamAtCluster = (AliMUONTrackParam*) fTrackParamAtCluster->UncheckedAt(iCluster1);
790 cluster = trackParamAtCluster->GetClusterPtr();
791 dX[iCluster1] = cluster->GetX() - trackParamAtCluster->GetNonBendingCoor();
792 dY[iCluster1] = cluster->GetY() - trackParamAtCluster->GetBendingCoor();
793 for (Int_t iCluster2 = 0; iCluster2 < iCluster1; iCluster2++) {
794 chi2 += ((*fClusterWeightsNonBending)(iCluster1, iCluster2) + (*fClusterWeightsNonBending)(iCluster2, iCluster1)) * dX[iCluster1] * dX[iCluster2] +
795 ((*fClusterWeightsBending)(iCluster1, iCluster2) + (*fClusterWeightsBending)(iCluster2, iCluster1)) * dY[iCluster1] * dY[iCluster2];
ea94c18b 796 }
96ebe67e 797 chi2 += ((*fClusterWeightsNonBending)(iCluster1, iCluster1) * dX[iCluster1] * dX[iCluster1]) +
798 ((*fClusterWeightsBending)(iCluster1, iCluster1) * dY[iCluster1] * dY[iCluster1]);
ea94c18b 799 }
800 delete [] dX;
801 delete [] dY;
802
803 } else {
804
96ebe67e 805 AliMUONVCluster *cluster;
ea94c18b 806 Double_t dX, dY;
96ebe67e 807 AliMUONTrackParam* trackParamAtCluster;
808 Int_t nClusters = GetNClusters();
809 for (Int_t iCluster = 0; iCluster < nClusters ; iCluster++) {
810 trackParamAtCluster = (AliMUONTrackParam*) fTrackParamAtCluster->UncheckedAt(iCluster);
811 cluster = trackParamAtCluster->GetClusterPtr();
812 dX = cluster->GetX() - trackParamAtCluster->GetNonBendingCoor();
813 dY = cluster->GetY() - trackParamAtCluster->GetBendingCoor();
814 chi2 += dX * dX / cluster->GetErrX2() + dY * dY / cluster->GetErrY2();
ea94c18b 815 }
816
817 }
818
819 return chi2;
820
821}
822
823 //__________________________________________________________________________
96ebe67e 824Bool_t AliMUONTrack::ComputeClusterWeights(TMatrixD* mcsCovariances)
ea94c18b 825{
96ebe67e 826 /// Compute the weight matrices of the attached clusters, in non bending and bending direction,
827 /// accounting for multiple scattering correlations and cluster resolution
ea94c18b 828 /// - Use the provided MCS covariance matrix if any (otherwise build it temporarily)
96ebe67e 829 /// - Assume that track parameters at each cluster are corrects
ea94c18b 830 /// - Return kFALSE if computation failed
7332f213 831 AliDebug(1,"Enter ComputeClusterWeights1");
832
833 if (!fTrackParamAtCluster) {
834 AliWarning("no cluster attached to this track");
835 return kFALSE;
836 }
ea94c18b 837
838 // Alocate memory
96ebe67e 839 Int_t nClusters = GetNClusters();
840 if (!fClusterWeightsNonBending) fClusterWeightsNonBending = new TMatrixD(nClusters,nClusters);
841 if (!fClusterWeightsBending) fClusterWeightsBending = new TMatrixD(nClusters,nClusters);
ea94c18b 842
843 // Compute weights matrices
96ebe67e 844 if (!ComputeClusterWeights(*fClusterWeightsNonBending, *fClusterWeightsBending, mcsCovariances)) return kFALSE;
ea94c18b 845
846 return kTRUE;
847
848}
849
850 //__________________________________________________________________________
c812dd48 851Bool_t AliMUONTrack::ComputeClusterWeights(TMatrixD& clusterWeightsNB, TMatrixD& clusterWeightsB,
57e2ad1a 852 TMatrixD* mcsCovariances, const AliMUONVCluster* discardedCluster) const
ea94c18b 853{
854 /// Compute the weight matrices, in non bending and bending direction,
96ebe67e 855 /// of the other attached clusters assuming the discarded one does not exist
856 /// accounting for multiple scattering correlations and cluster resolution
ea94c18b 857 /// - Use the provided MCS covariance matrix if any (otherwise build it temporarily)
858 /// - Return kFALSE if computation failed
7332f213 859 AliDebug(1,"Enter ComputeClusterWeights2");
ea94c18b 860
861 // Check MCS covariance matrix and recompute it if need
96ebe67e 862 Int_t nClusters = GetNClusters();
ea94c18b 863 Bool_t deleteMCSCov = kFALSE;
864 if (!mcsCovariances) {
96ebe67e 865 mcsCovariances = new TMatrixD(nClusters,nClusters);
ea94c18b 866 deleteMCSCov = kTRUE;
867 ComputeMCSCovariances(*mcsCovariances);
ea94c18b 868 }
869
870 // Resize the weights matrices; alocate memory
96ebe67e 871 if (discardedCluster) {
c812dd48 872 clusterWeightsNB.ResizeTo(nClusters-1,nClusters-1);
873 clusterWeightsB.ResizeTo(nClusters-1,nClusters-1);
ea94c18b 874 } else {
c812dd48 875 clusterWeightsNB.ResizeTo(nClusters,nClusters);
876 clusterWeightsB.ResizeTo(nClusters,nClusters);
ea94c18b 877 }
878
879 // Define variables
96ebe67e 880 AliMUONVCluster *cluster1, *cluster2;
881 Int_t iCurrentCluster1, iCurrentCluster2;
ea94c18b 882
883 // Compute the covariance matrices
96ebe67e 884 iCurrentCluster1 = 0;
885 for (Int_t iCluster1 = 0; iCluster1 < nClusters; iCluster1++) {
886 cluster1 = ((AliMUONTrackParam*) fTrackParamAtCluster->UncheckedAt(iCluster1))->GetClusterPtr();
ea94c18b 887
96ebe67e 888 if (cluster1 == discardedCluster) continue;
ea94c18b 889
96ebe67e 890 // Loop over next clusters
891 iCurrentCluster2 = iCurrentCluster1;
892 for (Int_t iCluster2 = iCluster1; iCluster2 < nClusters; iCluster2++) {
893 cluster2 = ((AliMUONTrackParam*) fTrackParamAtCluster->UncheckedAt(iCluster2))->GetClusterPtr();
ea94c18b 894
96ebe67e 895 if (cluster2 == discardedCluster) continue;
ea94c18b 896
ea94c18b 897 // Fill with MCS covariances
c812dd48 898 clusterWeightsNB(iCurrentCluster1, iCurrentCluster2) = (*mcsCovariances)(iCluster1,iCluster2);
ea94c18b 899
900 // Equal contribution from multiple scattering in non bending and bending directions
c812dd48 901 clusterWeightsB(iCurrentCluster1, iCurrentCluster2) = clusterWeightsNB(iCurrentCluster1, iCurrentCluster2);
ea94c18b 902
96ebe67e 903 // Add contribution from cluster resolution to diagonal element and symmetrize the matrix
904 if (iCurrentCluster1 == iCurrentCluster2) {
ea94c18b 905
906 // In non bending plane
c812dd48 907 clusterWeightsNB(iCurrentCluster1, iCurrentCluster1) += cluster1->GetErrX2();
ea94c18b 908 // In bending plane
c812dd48 909 clusterWeightsB(iCurrentCluster1, iCurrentCluster1) += cluster1->GetErrY2();
ea94c18b 910
911 } else {
912
913 // In non bending plane
c812dd48 914 clusterWeightsNB(iCurrentCluster2, iCurrentCluster1) = clusterWeightsNB(iCurrentCluster1, iCurrentCluster2);
ea94c18b 915 // In bending plane
c812dd48 916 clusterWeightsB(iCurrentCluster2, iCurrentCluster1) = clusterWeightsB(iCurrentCluster1, iCurrentCluster2);
ea94c18b 917
918 }
919
96ebe67e 920 iCurrentCluster2++;
ea94c18b 921 }
922
96ebe67e 923 iCurrentCluster1++;
ea94c18b 924 }
925
926 // Inversion of covariance matrices to get the weights
c812dd48 927 if (clusterWeightsNB.Determinant() != 0 && clusterWeightsB.Determinant() != 0) {
928 clusterWeightsNB.Invert();
929 clusterWeightsB.Invert();
ea94c18b 930 } else {
931 AliWarning(" Determinant = 0");
c812dd48 932 clusterWeightsNB.ResizeTo(0,0);
933 clusterWeightsB.ResizeTo(0,0);
ea94c18b 934 if(deleteMCSCov) delete mcsCovariances;
935 return kFALSE;
936 }
937
938 if(deleteMCSCov) delete mcsCovariances;
939
940 return kTRUE;
941
942}
943
944 //__________________________________________________________________________
945void AliMUONTrack::ComputeMCSCovariances(TMatrixD& mcsCovariances) const
946{
947 /// Compute the multiple scattering covariance matrix
96ebe67e 948 /// (assume that track parameters at each cluster are corrects)
7332f213 949 AliDebug(1,"Enter ComputeMCSCovariances");
ea94c18b 950
b709ac13 951 // Reset the size of the covariance matrix if needed
96ebe67e 952 Int_t nClusters = GetNClusters();
953 if (mcsCovariances.GetNrows() != nClusters) mcsCovariances.ResizeTo(nClusters,nClusters);
ea94c18b 954
955 // Define variables
b709ac13 956 Int_t nChambers = AliMUONConstants::NTrackingCh();
96ebe67e 957 AliMUONTrackParam* trackParamAtCluster;
ea94c18b 958 AliMUONTrackParam extrapTrackParam;
b709ac13 959 Int_t currentChamber = 0, expectedChamber = 0, size = 0;
960 Double_t *mcsAngle2 = new Double_t[2*nChambers];
961 Double_t *zMCS = new Double_t[2*nChambers];
96ebe67e 962 Int_t *indices = new Int_t[2*nClusters];
ea94c18b 963
964 // Compute multiple scattering dispersion angle at each chamber
965 // and save the z position where it is calculated
96ebe67e 966 for (Int_t iCluster = 0; iCluster < nClusters; iCluster++) {
967 trackParamAtCluster = (AliMUONTrackParam*) fTrackParamAtCluster->UncheckedAt(iCluster);
ea94c18b 968
969 // look for missing chambers if any
96ebe67e 970 currentChamber = trackParamAtCluster->GetClusterPtr()->GetChamberId();
ea94c18b 971 while (currentChamber > expectedChamber) {
972
973 // Save the z position where MCS dispersion is calculated
b709ac13 974 zMCS[size] = AliMUONConstants::DefaultChamberZ(expectedChamber);
ea94c18b 975
96ebe67e 976 // Do not take into account MCS in chambers prior the first cluster
977 if (iCluster > 0) {
ea94c18b 978
979 // Get track parameters at missing chamber z
96ebe67e 980 extrapTrackParam = *trackParamAtCluster;
981 AliMUONTrackExtrap::ExtrapToZ(&extrapTrackParam, zMCS[size]);
ea94c18b 982
983 // Save multiple scattering dispersion angle in missing chamber
4663da9f 984 mcsAngle2[size] = AliMUONTrackExtrap::GetMCSAngle2(extrapTrackParam,AliMUONConstants::ChamberThicknessInX0(expectedChamber),1.);
ea94c18b 985
b709ac13 986 } else mcsAngle2[size] = 0.;
ea94c18b 987
988 expectedChamber++;
b709ac13 989 size++;
ea94c18b 990 }
991
992 // Save z position where MCS dispersion is calculated
96ebe67e 993 zMCS[size] = trackParamAtCluster->GetZ();
ea94c18b 994
995 // Save multiple scattering dispersion angle in current chamber
4663da9f 996 mcsAngle2[size] = AliMUONTrackExtrap::GetMCSAngle2(*trackParamAtCluster,AliMUONConstants::ChamberThicknessInX0(currentChamber),1.);
b709ac13 997
998 // Save indice in zMCS array corresponding to the current cluster
96ebe67e 999 indices[iCluster] = size;
ea94c18b 1000
b709ac13 1001 expectedChamber = currentChamber + 1;
1002 size++;
ea94c18b 1003 }
1004
96ebe67e 1005 // complete array of z if last cluster is on the last but one chamber
b709ac13 1006 if (currentChamber != nChambers-1) zMCS[size++] = AliMUONConstants::DefaultChamberZ(nChambers-1);
ea94c18b 1007
1008 // Compute the covariance matrix
96ebe67e 1009 for (Int_t iCluster1 = 0; iCluster1 < nClusters; iCluster1++) {
ea94c18b 1010
96ebe67e 1011 for (Int_t iCluster2 = iCluster1; iCluster2 < nClusters; iCluster2++) {
ea94c18b 1012
1013 // Initialization to 0 (diagonal plus upper triangular part)
96ebe67e 1014 mcsCovariances(iCluster1,iCluster2) = 0.;
ea94c18b 1015
1016 // Compute contribution from multiple scattering in upstream chambers
96ebe67e 1017 for (Int_t k = 0; k < indices[iCluster1]; k++) {
1018 mcsCovariances(iCluster1,iCluster2) += (zMCS[indices[iCluster1]] - zMCS[k]) * (zMCS[indices[iCluster2]] - zMCS[k]) * mcsAngle2[k];
ea94c18b 1019 }
1020
1021 // Symetrize the matrix
96ebe67e 1022 mcsCovariances(iCluster2,iCluster1) = mcsCovariances(iCluster1,iCluster2);
ea94c18b 1023 }
1024
1025 }
1026
1027 delete [] mcsAngle2;
1028 delete [] zMCS;
b709ac13 1029 delete [] indices;
ea94c18b 1030
1031}
1032
208f139e 1033 //__________________________________________________________________________
5c15a68b 1034Int_t AliMUONTrack::ClustersInCommon(AliMUONTrack* track, Int_t stMin, Int_t stMax) const
208f139e 1035{
5c15a68b 1036 /// Returns the number of clusters in common in stations [stMin, stMax]
89c8d66d 1037 /// between the current track ("this") and the track pointed to by "track".
5c15a68b 1038
ce8cd162 1039 if (!track || !track->fTrackParamAtCluster || !this->fTrackParamAtCluster) return 0;
5c15a68b 1040
1041 Int_t chMin = 2 * stMin;
1042 Int_t chMax = 2 * stMax + 1;
89c8d66d 1043 Int_t clustersInCommon = 0;
5c15a68b 1044
89c8d66d 1045 // Loop over clusters of first track
5c15a68b 1046 Int_t nCl1 = this->GetNClusters();
1047 for(Int_t iCl1 = 0; iCl1 < nCl1; iCl1++) {
1048 AliMUONVCluster* cl1 = ((AliMUONTrackParam*) this->fTrackParamAtCluster->UncheckedAt(iCl1))->GetClusterPtr();
1049
1050 Int_t chCl1 = cl1->GetChamberId();
1051 if (chCl1 < chMin || chCl1 > chMax) continue;
1052
89c8d66d 1053 // Loop over clusters of second track
5c15a68b 1054 Int_t nCl2 = track->GetNClusters();
1055 for(Int_t iCl2 = 0; iCl2 < nCl2; iCl2++) {
1056 AliMUONVCluster* cl2 = ((AliMUONTrackParam*) track->fTrackParamAtCluster->UncheckedAt(iCl2))->GetClusterPtr();
1057
1058 Int_t chCl2 = cl2->GetChamberId();
1059 if (chCl2 < chMin || chCl2 > chMax) continue;
1060
1061 // Increment "clustersInCommon" if both clusters have the same ID
1062 if (cl1->GetUniqueID() == cl2->GetUniqueID()) {
89c8d66d 1063 clustersInCommon++;
1064 break;
1065 }
5c15a68b 1066
89c8d66d 1067 }
5c15a68b 1068
89c8d66d 1069 }
5c15a68b 1070
96ebe67e 1071 return clustersInCommon;
208f139e 1072}
1073
1074 //__________________________________________________________________________
5a240757 1075Int_t AliMUONTrack::GetNDF() const
1076{
1077 /// return the number of degrees of freedom
1078
1079 Int_t ndf = 2 * GetNClusters() - 5;
1080 return (ndf > 0) ? ndf : 0;
1081}
1082
1083 //__________________________________________________________________________
ea94c18b 1084Double_t AliMUONTrack::GetNormalizedChi2() const
1085{
5a240757 1086 /// return the chi2 value divided by the number of degrees of freedom (or FLT_MAX if ndf <= 0)
ea94c18b 1087
5a240757 1088 Double_t ndf = (Double_t) GetNDF();
4ea3f013 1089 return (ndf > 0.) ? fGlobalChi2 / ndf : 2.*MaxChi2();
ea94c18b 1090}
1091
1092 //__________________________________________________________________________
57e2ad1a 1093Int_t AliMUONTrack::FindCompatibleClusters(const AliMUONTrack &track, Double_t sigmaCut, Bool_t compatibleCluster[10]) const
b8dc484b 1094{
f202486b 1095 /// Try to match clusters from this track with clusters from the given track within the provided sigma cut:
1096 /// - Fill the array compatibleCluster[iCh] with kTRUE if a compatible cluster has been found in chamber iCh.
1097 /// - Return the number of clusters of "this" track matched with one cluster of the given track.
96ebe67e 1098 AliMUONVCluster *cluster1, *cluster2;
c59f70b9 1099 Double_t chi2, dX, dY;
61fed964 1100 Double_t chi2Max = sigmaCut * sigmaCut;
96ebe67e 1101
c59f70b9 1102 // initialization
1103 Int_t nMatchClusters = 0;
96ebe67e 1104 for ( Int_t ch = 0; ch < AliMUONConstants::NTrackingCh(); ch++) compatibleCluster[ch] = kFALSE;
b8dc484b 1105
f202486b 1106 if (!track.fTrackParamAtCluster || !this->fTrackParamAtCluster) return nMatchClusters;
7332f213 1107
96ebe67e 1108 // Loop over clusters of first track
ce8cd162 1109 Int_t nCl1 = this->GetNClusters();
1110 for(Int_t iCl1 = 0; iCl1 < nCl1; iCl1++) {
1111 cluster1 = static_cast<AliMUONTrackParam*>(this->fTrackParamAtCluster->UncheckedAt(iCl1))->GetClusterPtr();
96ebe67e 1112
1113 // Loop over clusters of second track
ce8cd162 1114 Int_t nCl2 = track.GetNClusters();
1115 for(Int_t iCl2 = 0; iCl2 < nCl2; iCl2++) {
1116 cluster2 = static_cast<AliMUONTrackParam*>(track.fTrackParamAtCluster->UncheckedAt(iCl2))->GetClusterPtr();
96ebe67e 1117
c59f70b9 1118 // check DE Id
1119 if (cluster1->GetDetElemId() != cluster2->GetDetElemId()) continue;
96ebe67e 1120
2e2d0c44 1121 // check local chi2
96ebe67e 1122 dX = cluster1->GetX() - cluster2->GetX();
96ebe67e 1123 dY = cluster1->GetY() - cluster2->GetY();
2e2d0c44 1124 chi2 = dX * dX / (cluster1->GetErrX2() + cluster2->GetErrX2()) + dY * dY / (cluster1->GetErrY2() + cluster2->GetErrY2());
1125 if (chi2 > 2. * chi2Max) continue; // 2 because 2 quantities in chi2
96ebe67e 1126
1127 compatibleCluster[cluster1->GetChamberId()] = kTRUE;
c59f70b9 1128 nMatchClusters++;
96ebe67e 1129 break;
1130 }
1131
b8dc484b 1132 }
1133
c59f70b9 1134 return nMatchClusters;
b8dc484b 1135}
a9e2aefa 1136
96ebe67e 1137//__________________________________________________________________________
f202486b 1138Bool_t AliMUONTrack::Match(AliMUONTrack &track, Double_t sigmaCut, Int_t &nMatchClusters) const
1139{
1140 /// Try to match this track with the given track. Matching conditions:
1141 /// - more than 50% of clusters from this track matched with clusters from the given track
1142 /// - at least 1 cluster matched before and 1 cluster matched after the dipole
1143
1144 Bool_t compTrack[10];
1145 nMatchClusters = FindCompatibleClusters(track, sigmaCut, compTrack);
1146
1147 if ((compTrack[0] || compTrack[1] || compTrack[2] || compTrack[3]) && // at least 1 cluster matched in st 1 & 2
1148 (compTrack[6] || compTrack[7] || compTrack[8] || compTrack[9]) && // at least 1 cluster matched in st 4 & 5
1149 2 * nMatchClusters > GetNClusters()) return kTRUE; // more than 50% of clusters matched
1150 else return kFALSE;
1151
1152}
1153
1154//__________________________________________________________________________
96ebe67e 1155void AliMUONTrack::SetTrackParamAtVertex(const AliMUONTrackParam* trackParam)
1156{
1157 /// set track parameters at vertex
1158 if (trackParam == 0x0) return;
1159 if (fTrackParamAtVertex) *fTrackParamAtVertex = *trackParam;
1160 else fTrackParamAtVertex = new AliMUONTrackParam(*trackParam);
1161}
1162
1163//__________________________________________________________________________
1164void AliMUONTrack::RecursiveDump() const
a9e2aefa 1165{
96ebe67e 1166 /// Recursive dump of AliMUONTrack, i.e. with dump of trackParamAtCluster and attached clusters
1167 AliMUONTrackParam *trackParamAtCluster;
1168 AliMUONVCluster *cluster;
de2cd600 1169 cout << "Recursive dump of Track: " << this << endl;
1170 // Track
1171 this->Dump();
96ebe67e 1172 for (Int_t iCluster = 0; iCluster < GetNClusters(); iCluster++) {
1173 trackParamAtCluster = (AliMUONTrackParam*) ((*fTrackParamAtCluster)[iCluster]);
1174 // trackParamAtCluster
1175 cout << "trackParamAtCluster: " << trackParamAtCluster << " (index: " << iCluster << ")" << endl;
1176 trackParamAtCluster->Dump();
1177 cluster = trackParamAtCluster->GetClusterPtr();
1178 // cluster
1179 cout << "cluster: " << cluster << endl;
1180 cluster->Print();
a9e2aefa 1181 }
de2cd600 1182 return;
a9e2aefa 1183}
04b5ea16 1184
6464217e 1185//_____________________________________________-
d2b1e7bb 1186void AliMUONTrack::Print(Option_t*) const
6464217e 1187{
2457f726 1188 /// Printing Track information
d2b1e7bb 1189
7b318ae2 1190 streamsize curW = cout.width();
1191 streamsize curPrecision = cout.precision();
96ebe67e 1192 cout << "<AliMUONTrack> No.Clusters=" << setw(2) << GetNClusters() <<
6464217e 1193 ", Match2Trig=" << setw(1) << GetMatchTrigger() <<
bdfb6eef 1194 ", LoTrgNum=" << setw(3) << LoCircuit() <<
d2b1e7bb 1195 ", Chi2-tracking-trigger=" << setw(8) << setprecision(5) << GetChi2MatchTrigger();
0a2dcc83 1196 cout << Form(" HitTriggerPattern trig %x track %x",fHitsPatternInTrigCh, fHitsPatternInTrigChTrk);
2e2d0c44 1197 cout << Form(" MClabel=%d",fTrackID) << endl;
7332f213 1198 if (fTrackParamAtCluster) fTrackParamAtCluster->First()->Print("FULL");
7b318ae2 1199 cout.width(curW);
1200 cout.precision(curPrecision);
6464217e 1201}
423b32ca 1202
1203//__________________________________________________________________________
00d46f24 1204void AliMUONTrack::SetLocalTrigger(Int_t loCirc, Int_t loStripX, Int_t loStripY, Int_t loDev, Int_t loLpt, Int_t loHpt, UChar_t respWithoutChamber)
423b32ca 1205{
1206 /// pack the local trigger information and store
1207
01413742 1208 if (loCirc < 0) return;
423b32ca 1209
1210 fLocalTrigger = 0;
1211 fLocalTrigger += loCirc;
1212 fLocalTrigger += loStripX << 8;
1213 fLocalTrigger += loStripY << 13;
1214 fLocalTrigger += loDev << 17;
1215 fLocalTrigger += loLpt << 22;
1216 fLocalTrigger += loHpt << 24;
00d46f24 1217 fLocalTrigger += respWithoutChamber << 26;
423b32ca 1218
1219}
1220
2e2d0c44 1221//__________________________________________________________________________
1222void AliMUONTrack::FindMCLabel()
1223{
1224 /// Determine the MC label from the label of the attached clusters and fill fMCLabel data member:
1225 /// More than 50% of clusters, including 1 before and 1 after the dipole, must share the same label
1226
1227 Int_t nClusters = GetNClusters();
1228 Int_t halfCluster = nClusters/2;
1229
1230 // reset MC label
1231 fTrackID = -1;
1232
1233 // loop over first clusters (if nClusters left < (nClusters-halfCluster) the conditions cannot be fulfilled)
1234 for (Int_t iCluster1 = 0; iCluster1 < nClusters-halfCluster; iCluster1++) {
1235 AliMUONVCluster* cluster1 = ((AliMUONTrackParam*) fTrackParamAtCluster->UncheckedAt(iCluster1))->GetClusterPtr();
1236
1237 // if the first cluster is not on station 1 or 2 the conditions cannot be fulfilled
1238 if (cluster1->GetChamberId() > 3) return;
1239
1240 Int_t label1 = cluster1->GetMCLabel();
1241 if (label1 < 0) continue;
1242
1243 Int_t nIdenticalLabel = 1;
1244
1245 // Loop over next clusters
1246 for (Int_t iCluster2 = iCluster1+1; iCluster2 < nClusters; iCluster2++) {
1247 AliMUONVCluster* cluster2 = ((AliMUONTrackParam*) fTrackParamAtCluster->UncheckedAt(iCluster2))->GetClusterPtr();
1248
1249 if (cluster2->GetMCLabel() != label1) continue;
1250
1251 nIdenticalLabel++;
1252
1253 // stop as soon as conditions are fulfilled
1254 if (nIdenticalLabel > halfCluster && cluster2->GetChamberId() > 5) {
1255 fTrackID = label1;
1256 return;
1257 }
1258
1259 }
1260
1261 }
1262
1263}
1264