Fix histo name
[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
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
04b5ea16 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
311 //__________________________________________________________________________
b1fea02e 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
334 //__________________________________________________________________________
7332f213 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
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
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
956019b6 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
447 AliMUONTrackExtrap::AddMCSEffect(trackParamAtCluster,AliMUONConstants::ChamberThicknessInX0(),1.);
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
455 AliMUONTrackExtrap::AddMCSEffect(trackParamAtCluster,AliMUONConstants::ChamberThicknessInX0(),1.);
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
208f139e 472}
473
ea94c18b 474 //__________________________________________________________________________
a0dc65b4 475Bool_t AliMUONTrack::IsValid(UInt_t requestedStationMask)
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)
7332f213 480
481 Int_t nClusters = GetNClusters();
482 AliMUONTrackParam *trackParam;
3d7acbb7 483 Int_t currentCh, currentSt, previousCh = -1, nChHitInSt45 = 0;
ad3c6eda 484 UInt_t presentStationMask(0);
7332f213 485
3d7acbb7 486 // first loop over clusters
ad3c6eda 487 for (Int_t i = 0; i < nClusters; i++) {
7332f213 488 trackParam = (AliMUONTrackParam*) fTrackParamAtCluster->UncheckedAt(i);
3d7acbb7 489
490 currentCh = trackParam->GetClusterPtr()->GetChamberId();
491 currentSt = currentCh/2;
492
493 // build present station mask
494 presentStationMask |= ( 1 << currentSt );
495
496 // count the number of chambers in station 4 & 5 that contain cluster(s)
497 if (currentCh > 5 && currentCh != previousCh) {
498 nChHitInSt45++;
499 previousCh = currentCh;
500 }
501
7332f213 502 }
503
3d7acbb7 504 return (((requestedStationMask & presentStationMask) == requestedStationMask) && (nChHitInSt45 >= 2));
7332f213 505}
506
507 //__________________________________________________________________________
a0dc65b4 508void AliMUONTrack::TagRemovableClusters(UInt_t requestedStationMask) {
7332f213 509 /// Identify clusters that can be removed from the track,
3d7acbb7 510 /// with the only requirements to have at least 1 cluster per requested station
511 /// and at least 2 chambers over 4 in stations 4 & 5 that contain cluster(s)
7332f213 512
513 Int_t nClusters = GetNClusters();
514 AliMUONTrackParam *trackParam, *nextTrackParam;
3d7acbb7 515 Int_t currentCh, nextCh, currentSt, nextSt, previousCh = -1, nChHitInSt45 = 0;
7332f213 516
3d7acbb7 517 // first loop over clusters
7332f213 518 for (Int_t i = 0; i < nClusters; i++) {
519 trackParam = (AliMUONTrackParam*) fTrackParamAtCluster->UncheckedAt(i);
3d7acbb7 520
521 currentCh = trackParam->GetClusterPtr()->GetChamberId();
522 currentSt = currentCh/2;
523
524 // reset flags to kFALSE for all clusters in required station
525 if ((1 << currentSt) & requestedStationMask) trackParam->SetRemovable(kFALSE);
9bf6860b 526 else trackParam->SetRemovable(kTRUE);
3d7acbb7 527
528 // count the number of chambers in station 4 & 5 that contain cluster(s)
529 if (currentCh > 5 && currentCh != previousCh) {
530 nChHitInSt45++;
531 previousCh = currentCh;
532 }
533
7332f213 534 }
535
3d7acbb7 536 // second loop over clusters
7332f213 537 for (Int_t i = 0; i < nClusters; i++) {
538 trackParam = (AliMUONTrackParam*) fTrackParamAtCluster->UncheckedAt(i);
539
540 currentCh = trackParam->GetClusterPtr()->GetChamberId();
3d7acbb7 541 currentSt = currentCh/2;
7332f213 542
3d7acbb7 543 // make sure they are more than 2 clusters in 2 different chambers of stations 4 & 5
544 // but 2 clusters in he same chamber will still be flagged as removable
545 if (nChHitInSt45 < 3 && currentSt > 2) {
7332f213 546
3d7acbb7 547 if (i == nClusters-1) {
548
549 trackParam->SetRemovable(kFALSE);
7332f213 550
3d7acbb7 551 } else {
552
553 nextTrackParam = (AliMUONTrackParam*) fTrackParamAtCluster->UncheckedAt(i+1);
554 nextCh = nextTrackParam->GetClusterPtr()->GetChamberId();
555
556 // set clusters in the same chamber as being removable
557 if (nextCh == currentCh) {
558 trackParam->SetRemovable(kTRUE);
559 nextTrackParam->SetRemovable(kTRUE);
560 i++; // skip cluster already checked
561 } else {
562 trackParam->SetRemovable(kFALSE);
563 }
564
565 }
7332f213 566
3d7acbb7 567 } else {
568
569 // skip clusters already flag as removable
570 if (trackParam->IsRemovable()) continue;
571
572 // loop over next track parameters
573 for (Int_t j = i+1; j < nClusters; j++) {
574 nextTrackParam = (AliMUONTrackParam*) fTrackParamAtCluster->UncheckedAt(j);
575
576 nextCh = nextTrackParam->GetClusterPtr()->GetChamberId();
577 nextSt = nextCh/2;
578
579 // set clusters in the same station as being removable
580 if (nextSt == currentSt) {
581 trackParam->SetRemovable(kTRUE);
582 nextTrackParam->SetRemovable(kTRUE);
583 i++; // skip cluster already checked
584 }
585
586 }
7332f213 587
588 }
3d7acbb7 589
7332f213 590 }
591
592}
593
594 //__________________________________________________________________________
ea94c18b 595Bool_t AliMUONTrack::ComputeLocalChi2(Bool_t accountForMCS)
596{
96ebe67e 597 /// Compute each cluster contribution to the chi2 of the track
ea94c18b 598 /// accounting for multiple scattering or not according to the flag
96ebe67e 599 /// - Also recompute the weight matrices of the attached clusters if accountForMCS=kTRUE
600 /// - Assume that track parameters at each cluster are corrects
ea94c18b 601 /// - Return kFALSE if computation failed
7332f213 602 AliDebug(1,"Enter ComputeLocalChi2");
603
604 if (!fTrackParamAtCluster) {
605 AliWarning("no cluster attached to this track");
606 return kFALSE;
607 }
ea94c18b 608
ea94c18b 609 if (accountForMCS) { // Compute local chi2 taking into account multiple scattering effects
610
611 // Compute MCS covariance matrix only once
96ebe67e 612 Int_t nClusters = GetNClusters();
613 TMatrixD mcsCovariances(nClusters,nClusters);
ea94c18b 614 ComputeMCSCovariances(mcsCovariances);
615
96ebe67e 616 // Make sure cluster weights are consistent with following calculations
617 if (!ComputeClusterWeights(&mcsCovariances)) {
ea94c18b 618 AliWarning("cannot take into account the multiple scattering effects");
619 return ComputeLocalChi2(kFALSE);
620 }
621
622 // Compute chi2 of the track
623 Double_t globalChi2 = ComputeGlobalChi2(kTRUE);
624 if (globalChi2 < 0.) return kFALSE;
625
96ebe67e 626 // Loop over removable clusters and compute their local chi2
627 AliMUONTrackParam* trackParamAtCluster1;
628 AliMUONVCluster *cluster, *discardedCluster;
629 Int_t iCluster1, iCluster2, iCurrentCluster1, iCurrentCluster2;
c812dd48 630 TMatrixD clusterWeightsNB(nClusters-1,nClusters-1);
631 TMatrixD clusterWeightsB(nClusters-1,nClusters-1);
96ebe67e 632 Double_t *dX = new Double_t[nClusters-1];
633 Double_t *dY = new Double_t[nClusters-1];
ea94c18b 634 Double_t globalChi2b;
96ebe67e 635 AliMUONTrackParam* trackParamAtCluster = (AliMUONTrackParam*) fTrackParamAtCluster->First();
636 while (trackParamAtCluster) {
ea94c18b 637
96ebe67e 638 discardedCluster = trackParamAtCluster->GetClusterPtr();
ea94c18b 639
96ebe67e 640 // Recompute cluster weights without the current cluster
c812dd48 641 if (!ComputeClusterWeights(clusterWeightsNB, clusterWeightsB, &mcsCovariances, discardedCluster)) {
ea94c18b 642 AliWarning("cannot take into account the multiple scattering effects");
96ebe67e 643 delete [] dX;
644 delete [] dY;
645 return ComputeLocalChi2(kFALSE);
ea94c18b 646 }
647
96ebe67e 648 // Compute track chi2 without the current cluster
ea94c18b 649 globalChi2b = 0.;
96ebe67e 650 iCurrentCluster1 = 0;
651 for (iCluster1 = 0; iCluster1 < nClusters ; iCluster1++) {
652 trackParamAtCluster1 = (AliMUONTrackParam*) fTrackParamAtCluster->UncheckedAt(iCluster1);
653 cluster = trackParamAtCluster1->GetClusterPtr();
ea94c18b 654
96ebe67e 655 if (cluster == discardedCluster) continue;
ea94c18b 656
657 // Compute and save residuals
96ebe67e 658 dX[iCurrentCluster1] = cluster->GetX() - trackParamAtCluster1->GetNonBendingCoor();
659 dY[iCurrentCluster1] = cluster->GetY() - trackParamAtCluster1->GetBendingCoor();
ea94c18b 660
96ebe67e 661 iCurrentCluster2 = 0;
662 for (iCluster2 = 0; iCluster2 < iCluster1; iCluster2++) {
663 cluster = ((AliMUONTrackParam*) fTrackParamAtCluster->UncheckedAt(iCluster2))->GetClusterPtr();
ea94c18b 664
96ebe67e 665 if (cluster == discardedCluster) continue;
ea94c18b 666
667 // Add contribution from covariances
c812dd48 668 globalChi2b += (clusterWeightsNB(iCurrentCluster1, iCurrentCluster2) +
669 clusterWeightsNB(iCurrentCluster2, iCurrentCluster1)) * dX[iCurrentCluster1] * dX[iCurrentCluster2] +
670 (clusterWeightsB(iCurrentCluster1, iCurrentCluster2) +
671 clusterWeightsB(iCurrentCluster2, iCurrentCluster1)) * dY[iCurrentCluster1] * dY[iCurrentCluster2];
ea94c18b 672
96ebe67e 673 iCurrentCluster2++;
ea94c18b 674 }
675
676 // Add contribution from variances
c812dd48 677 globalChi2b += clusterWeightsNB(iCurrentCluster1, iCurrentCluster1) * dX[iCurrentCluster1] * dX[iCurrentCluster1] +
678 clusterWeightsB(iCurrentCluster1, iCurrentCluster1) * dY[iCurrentCluster1] * dY[iCurrentCluster1];
ea94c18b 679
96ebe67e 680 iCurrentCluster1++;
ea94c18b 681 }
682
683 // Set local chi2
96ebe67e 684 trackParamAtCluster->SetLocalChi2(globalChi2 - globalChi2b);
ea94c18b 685
96ebe67e 686 trackParamAtCluster = (AliMUONTrackParam*) fTrackParamAtCluster->After(trackParamAtCluster);
ea94c18b 687 }
688
689 delete [] dX;
690 delete [] dY;
691
692 } else { // without multiple scattering effects
693
96ebe67e 694 AliMUONVCluster *discardedCluster;
ea94c18b 695 Double_t dX, dY;
96ebe67e 696 AliMUONTrackParam* trackParamAtCluster = (AliMUONTrackParam*) fTrackParamAtCluster->First();
697 while (trackParamAtCluster) {
ea94c18b 698
96ebe67e 699 discardedCluster = trackParamAtCluster->GetClusterPtr();
ea94c18b 700
701 // Compute residuals
96ebe67e 702 dX = discardedCluster->GetX() - trackParamAtCluster->GetNonBendingCoor();
703 dY = discardedCluster->GetY() - trackParamAtCluster->GetBendingCoor();
ea94c18b 704
705 // Set local chi2
96ebe67e 706 trackParamAtCluster->SetLocalChi2(dX * dX / discardedCluster->GetErrX2() + dY * dY / discardedCluster->GetErrY2());
ea94c18b 707
96ebe67e 708 trackParamAtCluster = (AliMUONTrackParam*) fTrackParamAtCluster->After(trackParamAtCluster);
ea94c18b 709 }
710
711 }
712
ea94c18b 713 return kTRUE;
714
715}
716
717 //__________________________________________________________________________
718Double_t AliMUONTrack::ComputeGlobalChi2(Bool_t accountForMCS)
719{
720 /// Compute the chi2 of the track accounting for multiple scattering or not according to the flag
96ebe67e 721 /// - Assume that track parameters at each cluster are corrects
722 /// - Assume the cluster weights matrices are corrects
4ea3f013 723 /// - Return a value of chi2 higher than the maximum allowed if computation failed
7332f213 724 AliDebug(1,"Enter ComputeGlobalChi2");
725
726 if (!fTrackParamAtCluster) {
727 AliWarning("no cluster attached to this track");
4ea3f013 728 return 2.*MaxChi2();
7332f213 729 }
ea94c18b 730
ea94c18b 731 Double_t chi2 = 0.;
732
733 if (accountForMCS) {
734
b709ac13 735 // Check the weight matrices. If weight matrices are not available compute chi2 without MCS
96ebe67e 736 if (!fClusterWeightsNonBending || !fClusterWeightsBending) {
737 AliWarning("cluster weights including multiple scattering effects are not available\n\t\t --> compute chi2 WITHOUT multiple scattering");
b709ac13 738 return ComputeGlobalChi2(kFALSE);
739 }
96ebe67e 740 Int_t nClusters = GetNClusters();
741 if (fClusterWeightsNonBending->GetNrows() != nClusters || fClusterWeightsBending->GetNcols() != nClusters) {
742 AliWarning("cluster weights including multiple scattering effects are not available\n\t\t --> compute chi2 WITHOUT multiple scattering");
ea94c18b 743 return ComputeGlobalChi2(kFALSE);
744 }
745
746 // Compute chi2
96ebe67e 747 AliMUONVCluster *cluster;
748 Double_t *dX = new Double_t[nClusters];
749 Double_t *dY = new Double_t[nClusters];
750 AliMUONTrackParam* trackParamAtCluster;
751 for (Int_t iCluster1 = 0; iCluster1 < nClusters; iCluster1++) {
752 trackParamAtCluster = (AliMUONTrackParam*) fTrackParamAtCluster->UncheckedAt(iCluster1);
753 cluster = trackParamAtCluster->GetClusterPtr();
754 dX[iCluster1] = cluster->GetX() - trackParamAtCluster->GetNonBendingCoor();
755 dY[iCluster1] = cluster->GetY() - trackParamAtCluster->GetBendingCoor();
756 for (Int_t iCluster2 = 0; iCluster2 < iCluster1; iCluster2++) {
757 chi2 += ((*fClusterWeightsNonBending)(iCluster1, iCluster2) + (*fClusterWeightsNonBending)(iCluster2, iCluster1)) * dX[iCluster1] * dX[iCluster2] +
758 ((*fClusterWeightsBending)(iCluster1, iCluster2) + (*fClusterWeightsBending)(iCluster2, iCluster1)) * dY[iCluster1] * dY[iCluster2];
ea94c18b 759 }
96ebe67e 760 chi2 += ((*fClusterWeightsNonBending)(iCluster1, iCluster1) * dX[iCluster1] * dX[iCluster1]) +
761 ((*fClusterWeightsBending)(iCluster1, iCluster1) * dY[iCluster1] * dY[iCluster1]);
ea94c18b 762 }
763 delete [] dX;
764 delete [] dY;
765
766 } else {
767
96ebe67e 768 AliMUONVCluster *cluster;
ea94c18b 769 Double_t dX, dY;
96ebe67e 770 AliMUONTrackParam* trackParamAtCluster;
771 Int_t nClusters = GetNClusters();
772 for (Int_t iCluster = 0; iCluster < nClusters ; iCluster++) {
773 trackParamAtCluster = (AliMUONTrackParam*) fTrackParamAtCluster->UncheckedAt(iCluster);
774 cluster = trackParamAtCluster->GetClusterPtr();
775 dX = cluster->GetX() - trackParamAtCluster->GetNonBendingCoor();
776 dY = cluster->GetY() - trackParamAtCluster->GetBendingCoor();
777 chi2 += dX * dX / cluster->GetErrX2() + dY * dY / cluster->GetErrY2();
ea94c18b 778 }
779
780 }
781
782 return chi2;
783
784}
785
786 //__________________________________________________________________________
96ebe67e 787Bool_t AliMUONTrack::ComputeClusterWeights(TMatrixD* mcsCovariances)
ea94c18b 788{
96ebe67e 789 /// Compute the weight matrices of the attached clusters, in non bending and bending direction,
790 /// accounting for multiple scattering correlations and cluster resolution
ea94c18b 791 /// - Use the provided MCS covariance matrix if any (otherwise build it temporarily)
96ebe67e 792 /// - Assume that track parameters at each cluster are corrects
ea94c18b 793 /// - Return kFALSE if computation failed
7332f213 794 AliDebug(1,"Enter ComputeClusterWeights1");
795
796 if (!fTrackParamAtCluster) {
797 AliWarning("no cluster attached to this track");
798 return kFALSE;
799 }
ea94c18b 800
801 // Alocate memory
96ebe67e 802 Int_t nClusters = GetNClusters();
803 if (!fClusterWeightsNonBending) fClusterWeightsNonBending = new TMatrixD(nClusters,nClusters);
804 if (!fClusterWeightsBending) fClusterWeightsBending = new TMatrixD(nClusters,nClusters);
ea94c18b 805
806 // Compute weights matrices
96ebe67e 807 if (!ComputeClusterWeights(*fClusterWeightsNonBending, *fClusterWeightsBending, mcsCovariances)) return kFALSE;
ea94c18b 808
809 return kTRUE;
810
811}
812
813 //__________________________________________________________________________
c812dd48 814Bool_t AliMUONTrack::ComputeClusterWeights(TMatrixD& clusterWeightsNB, TMatrixD& clusterWeightsB,
96ebe67e 815 TMatrixD* mcsCovariances, AliMUONVCluster* discardedCluster) const
ea94c18b 816{
817 /// Compute the weight matrices, in non bending and bending direction,
96ebe67e 818 /// of the other attached clusters assuming the discarded one does not exist
819 /// accounting for multiple scattering correlations and cluster resolution
ea94c18b 820 /// - Use the provided MCS covariance matrix if any (otherwise build it temporarily)
821 /// - Return kFALSE if computation failed
7332f213 822 AliDebug(1,"Enter ComputeClusterWeights2");
ea94c18b 823
824 // Check MCS covariance matrix and recompute it if need
96ebe67e 825 Int_t nClusters = GetNClusters();
ea94c18b 826 Bool_t deleteMCSCov = kFALSE;
827 if (!mcsCovariances) {
96ebe67e 828 mcsCovariances = new TMatrixD(nClusters,nClusters);
ea94c18b 829 deleteMCSCov = kTRUE;
830 ComputeMCSCovariances(*mcsCovariances);
ea94c18b 831 }
832
833 // Resize the weights matrices; alocate memory
96ebe67e 834 if (discardedCluster) {
c812dd48 835 clusterWeightsNB.ResizeTo(nClusters-1,nClusters-1);
836 clusterWeightsB.ResizeTo(nClusters-1,nClusters-1);
ea94c18b 837 } else {
c812dd48 838 clusterWeightsNB.ResizeTo(nClusters,nClusters);
839 clusterWeightsB.ResizeTo(nClusters,nClusters);
ea94c18b 840 }
841
842 // Define variables
96ebe67e 843 AliMUONVCluster *cluster1, *cluster2;
844 Int_t iCurrentCluster1, iCurrentCluster2;
ea94c18b 845
846 // Compute the covariance matrices
96ebe67e 847 iCurrentCluster1 = 0;
848 for (Int_t iCluster1 = 0; iCluster1 < nClusters; iCluster1++) {
849 cluster1 = ((AliMUONTrackParam*) fTrackParamAtCluster->UncheckedAt(iCluster1))->GetClusterPtr();
ea94c18b 850
96ebe67e 851 if (cluster1 == discardedCluster) continue;
ea94c18b 852
96ebe67e 853 // Loop over next clusters
854 iCurrentCluster2 = iCurrentCluster1;
855 for (Int_t iCluster2 = iCluster1; iCluster2 < nClusters; iCluster2++) {
856 cluster2 = ((AliMUONTrackParam*) fTrackParamAtCluster->UncheckedAt(iCluster2))->GetClusterPtr();
ea94c18b 857
96ebe67e 858 if (cluster2 == discardedCluster) continue;
ea94c18b 859
ea94c18b 860 // Fill with MCS covariances
c812dd48 861 clusterWeightsNB(iCurrentCluster1, iCurrentCluster2) = (*mcsCovariances)(iCluster1,iCluster2);
ea94c18b 862
863 // Equal contribution from multiple scattering in non bending and bending directions
c812dd48 864 clusterWeightsB(iCurrentCluster1, iCurrentCluster2) = clusterWeightsNB(iCurrentCluster1, iCurrentCluster2);
ea94c18b 865
96ebe67e 866 // Add contribution from cluster resolution to diagonal element and symmetrize the matrix
867 if (iCurrentCluster1 == iCurrentCluster2) {
ea94c18b 868
869 // In non bending plane
c812dd48 870 clusterWeightsNB(iCurrentCluster1, iCurrentCluster1) += cluster1->GetErrX2();
ea94c18b 871 // In bending plane
c812dd48 872 clusterWeightsB(iCurrentCluster1, iCurrentCluster1) += cluster1->GetErrY2();
ea94c18b 873
874 } else {
875
876 // In non bending plane
c812dd48 877 clusterWeightsNB(iCurrentCluster2, iCurrentCluster1) = clusterWeightsNB(iCurrentCluster1, iCurrentCluster2);
ea94c18b 878 // In bending plane
c812dd48 879 clusterWeightsB(iCurrentCluster2, iCurrentCluster1) = clusterWeightsB(iCurrentCluster1, iCurrentCluster2);
ea94c18b 880
881 }
882
96ebe67e 883 iCurrentCluster2++;
ea94c18b 884 }
885
96ebe67e 886 iCurrentCluster1++;
ea94c18b 887 }
888
889 // Inversion of covariance matrices to get the weights
c812dd48 890 if (clusterWeightsNB.Determinant() != 0 && clusterWeightsB.Determinant() != 0) {
891 clusterWeightsNB.Invert();
892 clusterWeightsB.Invert();
ea94c18b 893 } else {
894 AliWarning(" Determinant = 0");
c812dd48 895 clusterWeightsNB.ResizeTo(0,0);
896 clusterWeightsB.ResizeTo(0,0);
ea94c18b 897 if(deleteMCSCov) delete mcsCovariances;
898 return kFALSE;
899 }
900
901 if(deleteMCSCov) delete mcsCovariances;
902
903 return kTRUE;
904
905}
906
907 //__________________________________________________________________________
908void AliMUONTrack::ComputeMCSCovariances(TMatrixD& mcsCovariances) const
909{
910 /// Compute the multiple scattering covariance matrix
96ebe67e 911 /// (assume that track parameters at each cluster are corrects)
7332f213 912 AliDebug(1,"Enter ComputeMCSCovariances");
ea94c18b 913
b709ac13 914 // Reset the size of the covariance matrix if needed
96ebe67e 915 Int_t nClusters = GetNClusters();
916 if (mcsCovariances.GetNrows() != nClusters) mcsCovariances.ResizeTo(nClusters,nClusters);
ea94c18b 917
918 // Define variables
b709ac13 919 Int_t nChambers = AliMUONConstants::NTrackingCh();
96ebe67e 920 AliMUONTrackParam* trackParamAtCluster;
ea94c18b 921 AliMUONTrackParam extrapTrackParam;
b709ac13 922 Int_t currentChamber = 0, expectedChamber = 0, size = 0;
923 Double_t *mcsAngle2 = new Double_t[2*nChambers];
924 Double_t *zMCS = new Double_t[2*nChambers];
96ebe67e 925 Int_t *indices = new Int_t[2*nClusters];
ea94c18b 926
927 // Compute multiple scattering dispersion angle at each chamber
928 // and save the z position where it is calculated
96ebe67e 929 for (Int_t iCluster = 0; iCluster < nClusters; iCluster++) {
930 trackParamAtCluster = (AliMUONTrackParam*) fTrackParamAtCluster->UncheckedAt(iCluster);
ea94c18b 931
932 // look for missing chambers if any
96ebe67e 933 currentChamber = trackParamAtCluster->GetClusterPtr()->GetChamberId();
ea94c18b 934 while (currentChamber > expectedChamber) {
935
936 // Save the z position where MCS dispersion is calculated
b709ac13 937 zMCS[size] = AliMUONConstants::DefaultChamberZ(expectedChamber);
ea94c18b 938
96ebe67e 939 // Do not take into account MCS in chambers prior the first cluster
940 if (iCluster > 0) {
ea94c18b 941
942 // Get track parameters at missing chamber z
96ebe67e 943 extrapTrackParam = *trackParamAtCluster;
944 AliMUONTrackExtrap::ExtrapToZ(&extrapTrackParam, zMCS[size]);
ea94c18b 945
946 // Save multiple scattering dispersion angle in missing chamber
b709ac13 947 mcsAngle2[size] = AliMUONTrackExtrap::GetMCSAngle2(extrapTrackParam,AliMUONConstants::ChamberThicknessInX0(),1.);
ea94c18b 948
b709ac13 949 } else mcsAngle2[size] = 0.;
ea94c18b 950
951 expectedChamber++;
b709ac13 952 size++;
ea94c18b 953 }
954
955 // Save z position where MCS dispersion is calculated
96ebe67e 956 zMCS[size] = trackParamAtCluster->GetZ();
ea94c18b 957
958 // Save multiple scattering dispersion angle in current chamber
96ebe67e 959 mcsAngle2[size] = AliMUONTrackExtrap::GetMCSAngle2(*trackParamAtCluster,AliMUONConstants::ChamberThicknessInX0(),1.);
b709ac13 960
961 // Save indice in zMCS array corresponding to the current cluster
96ebe67e 962 indices[iCluster] = size;
ea94c18b 963
b709ac13 964 expectedChamber = currentChamber + 1;
965 size++;
ea94c18b 966 }
967
96ebe67e 968 // complete array of z if last cluster is on the last but one chamber
b709ac13 969 if (currentChamber != nChambers-1) zMCS[size++] = AliMUONConstants::DefaultChamberZ(nChambers-1);
ea94c18b 970
971 // Compute the covariance matrix
96ebe67e 972 for (Int_t iCluster1 = 0; iCluster1 < nClusters; iCluster1++) {
ea94c18b 973
96ebe67e 974 for (Int_t iCluster2 = iCluster1; iCluster2 < nClusters; iCluster2++) {
ea94c18b 975
976 // Initialization to 0 (diagonal plus upper triangular part)
96ebe67e 977 mcsCovariances(iCluster1,iCluster2) = 0.;
ea94c18b 978
979 // Compute contribution from multiple scattering in upstream chambers
96ebe67e 980 for (Int_t k = 0; k < indices[iCluster1]; k++) {
981 mcsCovariances(iCluster1,iCluster2) += (zMCS[indices[iCluster1]] - zMCS[k]) * (zMCS[indices[iCluster2]] - zMCS[k]) * mcsAngle2[k];
ea94c18b 982 }
983
984 // Symetrize the matrix
96ebe67e 985 mcsCovariances(iCluster2,iCluster1) = mcsCovariances(iCluster1,iCluster2);
ea94c18b 986 }
987
988 }
989
990 delete [] mcsAngle2;
991 delete [] zMCS;
b709ac13 992 delete [] indices;
ea94c18b 993
994}
995
208f139e 996 //__________________________________________________________________________
89c8d66d 997Int_t AliMUONTrack::ClustersInCommon(AliMUONTrack* track) const
208f139e 998{
96ebe67e 999 /// Returns the number of clusters in common between the current track ("this")
208f139e 1000 /// and the track pointed to by "track".
7332f213 1001 if (!fTrackParamAtCluster || !this->fTrackParamAtCluster) return 0;
89c8d66d 1002 Int_t nCluster1 = this->GetNClusters();
1003 Int_t nCluster2 = track->GetNClusters();
96ebe67e 1004 Int_t clustersInCommon = 0;
1005 AliMUONTrackParam *trackParamAtCluster1, *trackParamAtCluster2;
1006 // Loop over clusters of first track
89c8d66d 1007 for(Int_t iCluster1 = 0; iCluster1 < nCluster1; iCluster1++) {
1008 trackParamAtCluster1 = (AliMUONTrackParam*) this->fTrackParamAtCluster->UncheckedAt(iCluster1);
1009 // Loop over clusters of second track
1010 for(Int_t iCluster2 = 0; iCluster2 < nCluster2; iCluster2++) {
1011 trackParamAtCluster2 = (AliMUONTrackParam*) track->fTrackParamAtCluster->UncheckedAt(iCluster2);
1012 // Increment "clustersInCommon" if both trackParamAtCluster1 & 2 point to the same cluster
1013 if ((trackParamAtCluster1->GetClusterPtr()) == (trackParamAtCluster2->GetClusterPtr())) {
1014 clustersInCommon++;
1015 break;
1016 }
c59f70b9 1017 }
89c8d66d 1018 }
1019 return clustersInCommon;
1020}
1021
1022 //__________________________________________________________________________
1023Int_t AliMUONTrack::ClustersInCommonInSt345(AliMUONTrack* track) const
1024{
1025 /// Returns the number of clusters in common on stations 3, 4 and 5
1026 /// between the current track ("this") and the track pointed to by "track".
1027 if (!fTrackParamAtCluster || !this->fTrackParamAtCluster) return 0;
1028 Int_t nCluster1 = this->GetNClusters();
1029 Int_t nCluster2 = track->GetNClusters();
1030 Int_t clustersInCommon = 0;
1031 AliMUONTrackParam *trackParamAtCluster1, *trackParamAtCluster2;
1032 // Loop over clusters of first track
1033 for(Int_t iCluster1 = 0; iCluster1 < nCluster1; iCluster1++) {
1034 trackParamAtCluster1 = (AliMUONTrackParam*) this->fTrackParamAtCluster->UncheckedAt(iCluster1);
1035 if (trackParamAtCluster1->GetClusterPtr()->GetChamberId() < 4) continue;
1036 // Loop over clusters of second track
1037 for(Int_t iCluster2 = 0; iCluster2 < nCluster2; iCluster2++) {
1038 trackParamAtCluster2 = (AliMUONTrackParam*) track->fTrackParamAtCluster->UncheckedAt(iCluster2);
1039 if (trackParamAtCluster2->GetClusterPtr()->GetChamberId() < 4) continue;
1040 // Increment "clustersInCommon" if both trackParamAtCluster1 & 2 point to the same cluster
1041 if ((trackParamAtCluster1->GetClusterPtr()) == (trackParamAtCluster2->GetClusterPtr())) {
1042 clustersInCommon++;
1043 break;
1044 }
1045 }
1046 }
96ebe67e 1047 return clustersInCommon;
208f139e 1048}
1049
1050 //__________________________________________________________________________
5a240757 1051Int_t AliMUONTrack::GetNDF() const
1052{
1053 /// return the number of degrees of freedom
1054
1055 Int_t ndf = 2 * GetNClusters() - 5;
1056 return (ndf > 0) ? ndf : 0;
1057}
1058
1059 //__________________________________________________________________________
ea94c18b 1060Double_t AliMUONTrack::GetNormalizedChi2() const
1061{
5a240757 1062 /// return the chi2 value divided by the number of degrees of freedom (or FLT_MAX if ndf <= 0)
ea94c18b 1063
5a240757 1064 Double_t ndf = (Double_t) GetNDF();
4ea3f013 1065 return (ndf > 0.) ? fGlobalChi2 / ndf : 2.*MaxChi2();
ea94c18b 1066}
1067
1068 //__________________________________________________________________________
c59f70b9 1069Int_t AliMUONTrack::CompatibleTrack(AliMUONTrack* track, Double_t sigmaCut, Bool_t compatibleCluster[10]) const
b8dc484b 1070{
c59f70b9 1071 /// for each chamber: return kTRUE (kFALSE) if clusters are compatible (not compatible).
1072 /// nMatchClusters = number of clusters of "this" track matched with one cluster of track "track"
96ebe67e 1073 AliMUONTrackParam *trackParamAtCluster1, *trackParamAtCluster2;
1074 AliMUONVCluster *cluster1, *cluster2;
c59f70b9 1075 Double_t chi2, dX, dY;
61fed964 1076 Double_t chi2Max = sigmaCut * sigmaCut;
96ebe67e 1077
c59f70b9 1078 // initialization
1079 Int_t nMatchClusters = 0;
96ebe67e 1080 for ( Int_t ch = 0; ch < AliMUONConstants::NTrackingCh(); ch++) compatibleCluster[ch] = kFALSE;
b8dc484b 1081
c59f70b9 1082 if (!fTrackParamAtCluster || !this->fTrackParamAtCluster) return nMatchClusters;
7332f213 1083
96ebe67e 1084 // Loop over clusters of first track
1085 trackParamAtCluster1 = (AliMUONTrackParam*) this->fTrackParamAtCluster->First();
1086 while (trackParamAtCluster1) {
1087
1088 cluster1 = trackParamAtCluster1->GetClusterPtr();
1089
1090 // Loop over clusters of second track
1091 trackParamAtCluster2 = (AliMUONTrackParam*) track->fTrackParamAtCluster->First();
1092 while (trackParamAtCluster2) {
1093
1094 cluster2 = trackParamAtCluster2->GetClusterPtr();
1095
1096 //prepare next step
1097 trackParamAtCluster2 = (AliMUONTrackParam*) track->fTrackParamAtCluster->After(trackParamAtCluster2);
1098
c59f70b9 1099 // check DE Id
1100 if (cluster1->GetDetElemId() != cluster2->GetDetElemId()) continue;
96ebe67e 1101
2e2d0c44 1102 // check local chi2
96ebe67e 1103 dX = cluster1->GetX() - cluster2->GetX();
96ebe67e 1104 dY = cluster1->GetY() - cluster2->GetY();
2e2d0c44 1105 chi2 = dX * dX / (cluster1->GetErrX2() + cluster2->GetErrX2()) + dY * dY / (cluster1->GetErrY2() + cluster2->GetErrY2());
1106 if (chi2 > 2. * chi2Max) continue; // 2 because 2 quantities in chi2
96ebe67e 1107
1108 compatibleCluster[cluster1->GetChamberId()] = kTRUE;
c59f70b9 1109 nMatchClusters++;
96ebe67e 1110 break;
1111 }
1112
1113 trackParamAtCluster1 = (AliMUONTrackParam*) this->fTrackParamAtCluster->After(trackParamAtCluster1);
b8dc484b 1114 }
1115
c59f70b9 1116 return nMatchClusters;
b8dc484b 1117}
a9e2aefa 1118
96ebe67e 1119//__________________________________________________________________________
96ebe67e 1120void AliMUONTrack::SetTrackParamAtVertex(const AliMUONTrackParam* trackParam)
1121{
1122 /// set track parameters at vertex
1123 if (trackParam == 0x0) return;
1124 if (fTrackParamAtVertex) *fTrackParamAtVertex = *trackParam;
1125 else fTrackParamAtVertex = new AliMUONTrackParam(*trackParam);
1126}
1127
1128//__________________________________________________________________________
1129void AliMUONTrack::RecursiveDump() const
a9e2aefa 1130{
96ebe67e 1131 /// Recursive dump of AliMUONTrack, i.e. with dump of trackParamAtCluster and attached clusters
1132 AliMUONTrackParam *trackParamAtCluster;
1133 AliMUONVCluster *cluster;
de2cd600 1134 cout << "Recursive dump of Track: " << this << endl;
1135 // Track
1136 this->Dump();
96ebe67e 1137 for (Int_t iCluster = 0; iCluster < GetNClusters(); iCluster++) {
1138 trackParamAtCluster = (AliMUONTrackParam*) ((*fTrackParamAtCluster)[iCluster]);
1139 // trackParamAtCluster
1140 cout << "trackParamAtCluster: " << trackParamAtCluster << " (index: " << iCluster << ")" << endl;
1141 trackParamAtCluster->Dump();
1142 cluster = trackParamAtCluster->GetClusterPtr();
1143 // cluster
1144 cout << "cluster: " << cluster << endl;
1145 cluster->Print();
a9e2aefa 1146 }
de2cd600 1147 return;
a9e2aefa 1148}
04b5ea16 1149
6464217e 1150//_____________________________________________-
d2b1e7bb 1151void AliMUONTrack::Print(Option_t*) const
6464217e 1152{
2457f726 1153 /// Printing Track information
d2b1e7bb 1154
96ebe67e 1155 cout << "<AliMUONTrack> No.Clusters=" << setw(2) << GetNClusters() <<
6464217e 1156 ", Match2Trig=" << setw(1) << GetMatchTrigger() <<
c6ba19f7 1157 ", LoTrgNum=" << setw(3) << GetLoTrgNum() <<
d2b1e7bb 1158 ", Chi2-tracking-trigger=" << setw(8) << setprecision(5) << GetChi2MatchTrigger();
2e2d0c44 1159 cout << Form(" HitTriggerPattern %x",fHitsPatternInTrigCh);
1160 cout << Form(" MClabel=%d",fTrackID) << endl;
7332f213 1161 if (fTrackParamAtCluster) fTrackParamAtCluster->First()->Print("FULL");
6464217e 1162}
423b32ca 1163
1164//__________________________________________________________________________
1165void AliMUONTrack::SetLocalTrigger(Int_t loCirc, Int_t loStripX, Int_t loStripY, Int_t loDev, Int_t loLpt, Int_t loHpt)
1166{
1167 /// pack the local trigger information and store
1168
01413742 1169 if (loCirc < 0) return;
423b32ca 1170
1171 fLocalTrigger = 0;
1172 fLocalTrigger += loCirc;
1173 fLocalTrigger += loStripX << 8;
1174 fLocalTrigger += loStripY << 13;
1175 fLocalTrigger += loDev << 17;
1176 fLocalTrigger += loLpt << 22;
1177 fLocalTrigger += loHpt << 24;
1178
1179}
1180
2e2d0c44 1181//__________________________________________________________________________
1182void AliMUONTrack::FindMCLabel()
1183{
1184 /// Determine the MC label from the label of the attached clusters and fill fMCLabel data member:
1185 /// More than 50% of clusters, including 1 before and 1 after the dipole, must share the same label
1186
1187 Int_t nClusters = GetNClusters();
1188 Int_t halfCluster = nClusters/2;
1189
1190 // reset MC label
1191 fTrackID = -1;
1192
1193 // loop over first clusters (if nClusters left < (nClusters-halfCluster) the conditions cannot be fulfilled)
1194 for (Int_t iCluster1 = 0; iCluster1 < nClusters-halfCluster; iCluster1++) {
1195 AliMUONVCluster* cluster1 = ((AliMUONTrackParam*) fTrackParamAtCluster->UncheckedAt(iCluster1))->GetClusterPtr();
1196
1197 // if the first cluster is not on station 1 or 2 the conditions cannot be fulfilled
1198 if (cluster1->GetChamberId() > 3) return;
1199
1200 Int_t label1 = cluster1->GetMCLabel();
1201 if (label1 < 0) continue;
1202
1203 Int_t nIdenticalLabel = 1;
1204
1205 // Loop over next clusters
1206 for (Int_t iCluster2 = iCluster1+1; iCluster2 < nClusters; iCluster2++) {
1207 AliMUONVCluster* cluster2 = ((AliMUONTrackParam*) fTrackParamAtCluster->UncheckedAt(iCluster2))->GetClusterPtr();
1208
1209 if (cluster2->GetMCLabel() != label1) continue;
1210
1211 nIdenticalLabel++;
1212
1213 // stop as soon as conditions are fulfilled
1214 if (nIdenticalLabel > halfCluster && cluster2->GetChamberId() > 5) {
1215 fTrackID = label1;
1216 return;
1217 }
1218
1219 }
1220
1221 }
1222
1223}
1224