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