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