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