]> git.uio.no Git - u/mrichter/AliRoot.git/blame_incremental - MUON/AliMUONTrackK.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame_incremental
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (non-incremental) | history | HEAD
Optional geometry without CPV
[u/mrichter/AliRoot.git] / MUON / AliMUONTrackK.cxx
... / ...
CommitLineData
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 AliMUONTrackK
20// ---------------------
21// Reconstructed track in the muons system based on the extended
22// Kalman filter approach
23// Author: Alexander Zinchenko, JINR Dubna
24
25#include <Riostream.h>
26#include <TClonesArray.h>
27#include <TArrayD.h>
28#include <TMatrixD.h>
29
30#include "AliMUONTrackK.h"
31#include "AliMUON.h"
32#include "AliMUONConstants.h"
33
34#include "AliMUONTrackReconstructorK.h"
35#include "AliMagF.h"
36#include "AliMUONHitForRec.h"
37#include "AliMUONObjectPair.h"
38#include "AliMUONRawCluster.h"
39#include "AliMUONTrackParam.h"
40#include "AliMUONTrackExtrap.h"
41#include "AliMUONEventRecoCombi.h"
42#include "AliMUONDetElement.h"
43#include "AliRun.h"
44#include "AliLog.h"
45
46/// \cond CLASSIMP
47ClassImp(AliMUONTrackK) // Class implementation in ROOT context
48/// \endcond
49
50const Int_t AliMUONTrackK::fgkSize = 5;
51const Int_t AliMUONTrackK::fgkNSigma = 12;
52const Double_t AliMUONTrackK::fgkChi2max = 100;
53const Int_t AliMUONTrackK::fgkTriesMax = 10000;
54const Double_t AliMUONTrackK::fgkEpsilon = 0.002;
55
56void mnvertLocal(Double_t* a, Int_t l, Int_t m, Int_t n, Int_t& ifail);
57
58Int_t AliMUONTrackK::fgDebug = -1; //-1;
59Int_t AliMUONTrackK::fgNOfPoints = 0;
60AliMUONTrackReconstructorK* AliMUONTrackK::fgTrackReconstructor = NULL;
61TClonesArray* AliMUONTrackK::fgHitForRec = NULL;
62AliMUONEventRecoCombi *AliMUONTrackK::fgCombi = NULL;
63
64 //__________________________________________________________________________
65AliMUONTrackK::AliMUONTrackK()
66 : AliMUONTrack(),
67 fStartSegment(0x0),
68 fPosition(0.),
69 fPositionNew(0.),
70 fChi2(0.),
71 fTrackHits(0x0),
72 fNmbTrackHits(0),
73 fTrackDir(1),
74 fBPFlag(kFALSE),
75 fRecover(0),
76 fSkipHit(0x0),
77 fTrackPar(0x0),
78 fTrackParNew(0x0),
79 fCovariance(0x0),
80 fWeight(0x0),
81 fParExtrap(0x0),
82 fParFilter(0x0),
83 fParSmooth(0x0),
84 fCovExtrap(0x0),
85 fCovFilter(0x0),
86 fJacob(0x0),
87 fNSteps(0),
88 fSteps(0x0),
89 fChi2Array(0x0),
90 fChi2Smooth(0x0)
91{
92/// Default constructor
93
94 fgTrackReconstructor = NULL; // pointer to event reconstructor
95 fgHitForRec = NULL; // pointer to points
96 fgNOfPoints = 0; // number of points
97
98 return;
99}
100
101 //__________________________________________________________________________
102AliMUONTrackK::AliMUONTrackK(AliMUONTrackReconstructorK *TrackReconstructor, TClonesArray *hitForRec)
103 : AliMUONTrack(),
104 fStartSegment(0x0),
105 fPosition(0.),
106 fPositionNew(0.),
107 fChi2(0.),
108 fTrackHits(0x0),
109 fNmbTrackHits(0),
110 fTrackDir(1),
111 fBPFlag(kFALSE),
112 fRecover(0),
113 fSkipHit(0x0),
114 fTrackPar(0x0),
115 fTrackParNew(0x0),
116 fCovariance(0x0),
117 fWeight(0x0),
118 fParExtrap(0x0),
119 fParFilter(0x0),
120 fParSmooth(0x0),
121 fCovExtrap(0x0),
122 fCovFilter(0x0),
123 fJacob(0x0),
124 fNSteps(0),
125 fSteps(0x0),
126 fChi2Array(0x0),
127 fChi2Smooth(0x0)
128{
129/// Constructor
130
131 if (!TrackReconstructor) return;
132 fgTrackReconstructor = TrackReconstructor; // pointer to event reconstructor
133 fgHitForRec = hitForRec; // pointer to points
134 fgNOfPoints = fgHitForRec->GetEntriesFast(); // number of points
135 fgCombi = NULL;
136 if (fgTrackReconstructor->GetTrackMethod() == 3) fgCombi = AliMUONEventRecoCombi::Instance();
137}
138
139 //__________________________________________________________________________
140AliMUONTrackK::AliMUONTrackK(AliMUONObjectPair *segment)
141 : AliMUONTrack(NULL,NULL),
142 fStartSegment(new AliMUONObjectPair(*segment)),
143 fPosition(0.),
144 fPositionNew(0.),
145 fChi2(0.),
146 fTrackHits(new TObjArray(13)),
147 fNmbTrackHits(2),
148 fTrackDir(1),
149 fBPFlag(kFALSE),
150 fRecover(0),
151 fSkipHit(0x0),
152 fTrackPar(new TMatrixD(fgkSize,1)),
153 fTrackParNew(new TMatrixD(fgkSize,1)),
154 fCovariance(new TMatrixD(fgkSize,fgkSize)),
155 fWeight(new TMatrixD(fgkSize,fgkSize)),
156 fParExtrap(new TObjArray(15)),
157 fParFilter(new TObjArray(15)),
158 fParSmooth(0x0),
159 fCovExtrap(new TObjArray(15)),
160 fCovFilter(new TObjArray(15)),
161 fJacob(new TObjArray(15)),
162 fNSteps(0),
163 fSteps(new TArrayD(15)),
164 fChi2Array(new TArrayD(13)),
165 fChi2Smooth(0x0)
166{
167/// Constructor from a segment
168 Double_t dX, dY, dZ, bendingSlope, bendingImpact;
169 AliMUONHitForRec *hit1, *hit2;
170 AliMUONRawCluster *clus;
171 TClonesArray *rawclusters;
172
173 // Pointers to hits from the segment
174 hit1 = (AliMUONHitForRec*) segment->First();
175 hit2 = (AliMUONHitForRec*) segment->Second();
176 hit1->SetNTrackHits(hit1->GetNTrackHits()+1); // mark hit as being on track
177 hit2->SetNTrackHits(hit2->GetNTrackHits()+1); // mark hit as being on track
178 // check sorting in Z
179 if (TMath::Abs(hit1->GetZ()) > TMath::Abs(hit2->GetZ())) {
180 hit1 = hit2;
181 hit2 = (AliMUONHitForRec*) segment->First();
182 }
183
184 // Fill array of track parameters
185 if (hit1->GetChamberNumber() > 7) {
186 // last tracking station
187 (*fTrackPar)(0,0) = hit1->GetBendingCoor(); // y
188 (*fTrackPar)(1,0) = hit1->GetNonBendingCoor(); // x
189 fPosition = hit1->GetZ(); // z
190 fTrackHits->Add(hit2); // add hit 2
191 fTrackHits->Add(hit1); // add hit 1
192 //AZ(Z->-Z) fTrackDir = -1;
193 fTrackDir = 1;
194 } else {
195 // last but one tracking station
196 (*fTrackPar)(0,0) = hit2->GetBendingCoor(); // y
197 (*fTrackPar)(1,0) = hit2->GetNonBendingCoor(); // x
198 fPosition = hit2->GetZ(); // z
199 fTrackHits->Add(hit1); // add hit 1
200 fTrackHits->Add(hit2); // add hit 2
201 //AZ(Z->-Z) fTrackDir = 1;
202 fTrackDir = -1;
203 }
204 dZ = hit2->GetZ() - hit1->GetZ();
205 dY = hit2->GetBendingCoor() - hit1->GetBendingCoor();
206 dX = hit2->GetNonBendingCoor() - hit1->GetNonBendingCoor();
207 bendingSlope = (hit2->GetBendingCoor() - hit1->GetBendingCoor()) / dZ;
208 bendingImpact = hit1->GetBendingCoor() - hit1->GetZ() * bendingSlope;
209 (*fTrackPar)(2,0) = TMath::ATan2(dY,dZ); // alpha
210 if ((*fTrackPar)(2,0) < 0.) (*fTrackPar)(2,0) += 2*TMath::Pi(); // from 0 to 2*pi
211 (*fTrackPar)(3,0) = TMath::ATan2(-dX,dZ/TMath::Cos((*fTrackPar)(2,0))); // beta
212 (*fTrackPar)(2,0) -= TMath::Pi();
213 (*fTrackPar)(4,0) = 1./AliMUONTrackExtrap::GetBendingMomentumFromImpactParam(bendingImpact); // 1/Pt
214 (*fTrackPar)(4,0) *= TMath::Cos((*fTrackPar)(3,0)); // 1/p
215 // Evaluate covariance (and weight) matrix
216 EvalCovariance(dZ);
217
218 if (fgDebug < 0 ) return;
219 cout << fgTrackReconstructor->GetNRecTracks()-1 << " "
220 << AliMUONTrackExtrap::GetBendingMomentumFromImpactParam(bendingImpact)
221 << " " << 1/(*fTrackPar)(4,0) << " ";
222 // from raw clusters
223 for (Int_t i=0; i<2; i++) {
224 hit1 = (AliMUONHitForRec*) ((*fTrackHits)[i]);
225 rawclusters = fgTrackReconstructor->GetMUONData()->RawClusters(hit1->GetChamberNumber());
226 clus = (AliMUONRawCluster*) rawclusters->UncheckedAt(hit1->GetHitNumber());
227 cout << clus->GetTrack(1);
228 if (clus->GetTrack(2) != -1) cout << " " << clus->GetTrack(2);
229 if (i == 0) cout << " <--> ";
230 }
231 cout << " @ " << hit1->GetChamberNumber() << endl;
232
233}
234
235 //__________________________________________________________________________
236AliMUONTrackK::~AliMUONTrackK()
237{
238/// Destructor
239
240 if (fStartSegment) {
241 delete fStartSegment;
242 fStartSegment = 0x0;
243 }
244
245 if (fTrackHits) {
246 //cout << fNmbTrackHits << endl;
247 for (Int_t i = 0; i < fNmbTrackHits; i++) {
248 AliMUONHitForRec *hit = (AliMUONHitForRec*) fTrackHits->UncheckedAt(i);
249 hit->SetNTrackHits(hit->GetNTrackHits()-1);
250 }
251 delete fTrackHits; // delete the TObjArray of pointers to TrackHit's
252 fTrackHits = NULL;
253 }
254 if (fTrackPar) {
255 delete fTrackPar; delete fTrackParNew; delete fCovariance;
256 delete fWeight;
257 }
258
259 if (fSteps) delete fSteps;
260 if (fChi2Array) delete fChi2Array;
261 if (fChi2Smooth) delete fChi2Smooth;
262 if (fParSmooth) {fParSmooth->Delete(); delete fParSmooth; }
263 // Delete only matrices not shared by several tracks
264 if (!fParExtrap) return;
265
266 Int_t id = 0;
267 for (Int_t i=0; i<fNSteps; i++) {
268 //if (fParExtrap->UncheckedAt(i)->GetUniqueID() > 1)
269 // fParExtrap->UncheckedAt(i)->SetUniqueID(fParExtrap->RemoveAt(i)->GetUniqueID()-1);
270 id = fParExtrap->UncheckedAt(i)->GetUniqueID();
271 if (id > 1) {
272 fParExtrap->UncheckedAt(i)->SetUniqueID(id-1);
273 fParExtrap->RemoveAt(i);
274 }
275 //if (fParFilter->UncheckedAt(i)->GetUniqueID() > 1)
276 // fParFilter->UncheckedAt(i)->SetUniqueID(fParFilter->RemoveAt(i)->GetUniqueID()-1);
277 id = fParFilter->UncheckedAt(i)->GetUniqueID();
278 if (id > 1) {
279 fParFilter->UncheckedAt(i)->SetUniqueID(id-1);
280 fParFilter->RemoveAt(i);
281 }
282 //if (fCovExtrap->UncheckedAt(i)->GetUniqueID() > 1)
283 // fCovExtrap->UncheckedAt(i)->SetUniqueID(fCovExtrap->RemoveAt(i)->GetUniqueID()-1);
284 id = fCovExtrap->UncheckedAt(i)->GetUniqueID();
285 if (id > 1) {
286 fCovExtrap->UncheckedAt(i)->SetUniqueID(id-1);
287 fCovExtrap->RemoveAt(i);
288 }
289
290 //if (fCovFilter->UncheckedAt(i)->GetUniqueID() > 1)
291 // fCovFilter->UncheckedAt(i)->SetUniqueID(fCovFilter->RemoveAt(i)->GetUniqueID()-1);
292 id = fCovFilter->UncheckedAt(i)->GetUniqueID();
293 if (id > 1) {
294 fCovFilter->UncheckedAt(i)->SetUniqueID(id-1);
295 fCovFilter->RemoveAt(i);
296 }
297 //if (fJacob->UncheckedAt(i)->GetUniqueID() > 1)
298 // fJacob->UncheckedAt(i)->SetUniqueID(fJacob->RemoveAt(i)->GetUniqueID()-1);
299 id = fJacob->UncheckedAt(i)->GetUniqueID();
300 if (id > 1) {
301 fJacob->UncheckedAt(i)->SetUniqueID(id-1);
302 fJacob->RemoveAt(i);
303 }
304 }
305 /*
306 for (Int_t i=0; i<fNSteps; i++) {
307 if (fParExtrap->UncheckedAt(i)) ((TMatrixD*)fParExtrap->UncheckedAt(i))->Delete();
308 if (fParFilter->UncheckedAt(i)) ((TMatrixD*)fParFilter->UncheckedAt(i))->Delete();
309 if (fCovExtrap->UncheckedAt(i)) ((TMatrixD*)fCovExtrap->UncheckedAt(i))->Delete();
310 cout << fCovFilter->UncheckedAt(i) << " " << (*fSteps)[i] << endl;
311 if (fCovFilter->UncheckedAt(i)) ((TMatrixD*)fCovFilter->UncheckedAt(i))->Delete();
312 if (fJacob->UncheckedAt(i)) ((TMatrixD*)fJacob->UncheckedAt(i))->Delete();
313 }
314 */
315 fParExtrap->Delete(); fParFilter->Delete();
316 fCovExtrap->Delete(); fCovFilter->Delete();
317 fJacob->Delete();
318 delete fParExtrap; delete fParFilter;
319 delete fCovExtrap; delete fCovFilter;
320 delete fJacob;
321}
322
323 //__________________________________________________________________________
324AliMUONTrackK & AliMUONTrackK::operator=(const AliMUONTrackK& source)
325{
326/// Assignment operator
327
328 // Members
329 if(&source == this) return *this;
330
331 // base class assignement
332 //AZ TObject::operator=(source);
333 AliMUONTrack::operator=(source);
334
335 if (fStartSegment) delete fStartSegment;
336 if (source.fStartSegment) fStartSegment = new AliMUONObjectPair(*(source.fStartSegment));
337 else fStartSegment = 0x0;
338
339 fNmbTrackHits = source.fNmbTrackHits;
340 fChi2 = source.fChi2;
341 fPosition = source.fPosition;
342 fPositionNew = source.fPositionNew;
343 fTrackDir = source.fTrackDir;
344 fBPFlag = source.fBPFlag;
345 fRecover = source.fRecover;
346 fSkipHit = source.fSkipHit;
347
348 // Pointers
349 fTrackHits = new TObjArray(*source.fTrackHits);
350 //source.fTrackHits->Dump();
351 //fTrackHits->Dump();
352 for (Int_t i = 0; i < fNmbTrackHits; i++) {
353 AliMUONHitForRec *hit = (AliMUONHitForRec*) fTrackHits->UncheckedAt(i);
354 hit->SetNTrackHits(hit->GetNTrackHits()+1);
355 }
356
357 fTrackPar = new TMatrixD(*source.fTrackPar); // track parameters
358 fTrackParNew = new TMatrixD(*source.fTrackParNew); // track parameters
359 fCovariance = new TMatrixD(*source.fCovariance); // covariance matrix
360 fWeight = new TMatrixD(*source.fWeight); // weight matrix (inverse of covariance)
361
362 // For smoother
363 fParExtrap = new TObjArray(*source.fParExtrap);
364 fParFilter = new TObjArray(*source.fParFilter);
365 fCovExtrap = new TObjArray(*source.fCovExtrap);
366 fCovFilter = new TObjArray(*source.fCovFilter);
367 fJacob = new TObjArray(*source.fJacob);
368 fSteps = new TArrayD(*source.fSteps);
369 fNSteps = source.fNSteps;
370 fChi2Array = NULL;
371 if (source.fChi2Array) fChi2Array = new TArrayD(*source.fChi2Array);
372 fChi2Smooth = NULL;
373 fParSmooth = NULL;
374
375 for (Int_t i=0; i<fParExtrap->GetEntriesFast(); i++) {
376 fParExtrap->UncheckedAt(i)->SetUniqueID(fParExtrap->UncheckedAt(i)->GetUniqueID()+1);
377 fParFilter->UncheckedAt(i)->SetUniqueID(fParFilter->UncheckedAt(i)->GetUniqueID()+1);
378 fCovExtrap->UncheckedAt(i)->SetUniqueID(fCovExtrap->UncheckedAt(i)->GetUniqueID()+1);
379 fCovFilter->UncheckedAt(i)->SetUniqueID(fCovFilter->UncheckedAt(i)->GetUniqueID()+1);
380 fJacob->UncheckedAt(i)->SetUniqueID(fJacob->UncheckedAt(i)->GetUniqueID()+1);
381 }
382 return *this;
383}
384
385 //__________________________________________________________________________
386void AliMUONTrackK::EvalCovariance(Double_t dZ)
387{
388/// Evaluate covariance (and weight) matrix for track candidate
389 Double_t sigmaB, sigmaNonB, tanA, tanB, dAdY, rad, dBdX, dBdY;
390
391 dZ = -dZ;
392 sigmaB = fgTrackReconstructor->GetBendingResolution(); // bending resolution
393 sigmaNonB = fgTrackReconstructor->GetNonBendingResolution(); // non-bending resolution
394
395 (*fWeight)(0,0) = sigmaB*sigmaB; // <yy>
396
397 (*fWeight)(1,1) = sigmaNonB*sigmaNonB; // <xx>
398
399 tanA = TMath::Tan((*fTrackPar)(2,0));
400 dAdY = 1/(1+tanA*tanA)/dZ;
401 (*fWeight)(2,2) = dAdY*dAdY*(*fWeight)(0,0)*2; // <aa>
402 (*fWeight)(0,2) = dAdY*(*fWeight)(0,0); // <ya>
403 (*fWeight)(2,0) = (*fWeight)(0,2);
404
405 rad = dZ/TMath::Cos((*fTrackPar)(2,0));
406 tanB = TMath::Tan((*fTrackPar)(3,0));
407 dBdX = 1/(1+tanB*tanB)/rad;
408 dBdY = 0; // neglect
409 (*fWeight)(3,3) = dBdX*dBdX*(*fWeight)(1,1)*2; // <bb>
410 (*fWeight)(1,3) = dBdX*(*fWeight)(1,1); // <xb>
411 (*fWeight)(3,1) = (*fWeight)(1,3);
412
413 (*fWeight)(4,4) = ((*fTrackPar)(4,0)*0.5)*((*fTrackPar)(4,0)*0.5); // error 50%
414
415 // check whether the Invert method returns flag if matrix cannot be inverted,
416 // and do not calculate the Determinant in that case !!!!
417 if (fWeight->Determinant() != 0) {
418 // fWeight->Invert();
419 Int_t ifail;
420 mnvertLocal(&((*fWeight)(0,0)), fgkSize,fgkSize,fgkSize,ifail);
421 } else {
422 AliWarning(" Determinant fWeight=0:");
423 }
424 return;
425}
426
427 //__________________________________________________________________________
428Bool_t AliMUONTrackK::KalmanFilter(Int_t ichamBeg, Int_t ichamEnd, Bool_t Back, Double_t zDipole1, Double_t zDipole2)
429{
430/// Follows track through detector stations
431 Bool_t miss, success;
432 Int_t ichamb, iFB, iMin, iMax, dChamb, ichambOK;
433 Int_t ihit, firstIndx = -1, lastIndx = -1, currIndx = -1, iz0 = -1, iz = -1;
434 Double_t zEnd, dChi2;
435 AliMUONHitForRec *hitAdd, *firstHit = NULL, *lastHit = NULL, *hit;
436 AliMUONRawCluster *clus;
437 TClonesArray *rawclusters;
438 hit = 0; clus = 0; rawclusters = 0;
439
440 miss = success = kTRUE;
441 Int_t endOfProp = 0;
442 //iFB = (ichamEnd == ichamBeg) ? fTrackDir : TMath::Sign(1,ichamEnd-ichamBeg);
443 iFB = (ichamEnd == ichamBeg) ? -fTrackDir : TMath::Sign(1,ichamEnd-ichamBeg);
444 iMin = TMath::Min(ichamEnd,ichamBeg);
445 iMax = TMath::Max(ichamEnd,ichamBeg);
446 ichamb = ichamBeg;
447 ichambOK = ichamb;
448
449 if (Back) {
450 // backpropagation
451 currIndx = 1;
452 if (((AliMUONHitForRec*)fTrackHits->First())->GetChamberNumber() != ichamb) currIndx = 0;
453 } else if (fRecover) {
454 hit = GetHitLastOk();
455 currIndx = fTrackHits->IndexOf(hit);
456 if (currIndx < 0) hit = (AliMUONHitForRec*) fStartSegment->First(); // for station 3
457 Back = kTRUE;
458 ichamb = hit->GetChamberNumber();
459 if (hit == fSkipHit || fRecover == 2 && currIndx >= 0) {
460 // start from the last point or outlier
461 // remove the skipped hit
462 fTrackHits->Remove(fSkipHit); // remove hit
463 fNmbTrackHits --;
464 fSkipHit->SetNTrackHits(fSkipHit->GetNTrackHits()-1); // unmark hit
465 if (fRecover == 1) {
466 // recovery
467 Back = kFALSE;
468 fRecover = 0;
469 ichambOK = ((AliMUONHitForRec*)((*fTrackHits)[fNmbTrackHits-1]))->GetChamberNumber();
470 //if (ichambOK >= 8 && ((AliMUONHitForRec*)((*fTrackHits)[0]))->GetChamberNumber() == 6) ichambOK = 6;
471 if (fgTrackReconstructor->GetTrackMethod() == 3 &&
472 fSkipHit->GetHitNumber() < 0) {
473 iz0 = fgCombi->IZfromHit(fSkipHit);
474 currIndx = -1;
475 }
476 else currIndx = fgHitForRec->IndexOf(fSkipHit);
477 } else {
478 // outlier
479 fTrackHits->Compress();
480 }
481 } // if (hit == fSkipHit)
482 else if (currIndx < 0) currIndx = fTrackHits->IndexOf(hit);
483 } // else if (fRecover)
484 else {
485 // Get indices of the 1'st and last hits on the track candidate
486 firstHit = (AliMUONHitForRec*) fTrackHits->First();
487 lastHit = (AliMUONHitForRec*) fTrackHits->Last();
488 if (fgTrackReconstructor->GetTrackMethod() == 3 &&
489 lastHit->GetHitNumber() < 0) iz0 = fgCombi->IZfromHit(lastHit);
490 else {
491 firstIndx = fgHitForRec->IndexOf(firstHit);
492 lastIndx = fgHitForRec->IndexOf(lastHit);
493 currIndx = TMath::Abs (TMath::Max(firstIndx*iFB,lastIndx*iFB));
494 }
495 }
496
497 if (iz0 < 0) iz0 = iFB;
498 while (ichamb >= iMin && ichamb <= iMax) {
499 // Find the closest hit in Z, not belonging to the current plane
500 if (Back) {
501 // backpropagation
502 if (currIndx < fNmbTrackHits) {
503 hitAdd = (AliMUONHitForRec*) fTrackHits->UncheckedAt(currIndx);
504 zEnd = hitAdd->GetZ();
505 //AZ(z->-z) } else zEnd = -9999;
506 } else zEnd = 9999;
507 } else {
508 //AZ(Z->-Z) zEnd = -9999;
509 zEnd = 9999;
510 for (ihit=currIndx+iFB; ihit>=0 && ihit<fgNOfPoints; ihit+=iFB) {
511 hitAdd = (AliMUONHitForRec*) ((*fgHitForRec)[ihit]);
512 //if (TMath::Abs(hitAdd->GetZ()-fPosition) > 0.1) {
513 if (TMath::Abs(hitAdd->GetZ()-fPosition) > 0.5) {
514 zEnd = hitAdd->GetZ();
515 currIndx = ihit;
516 break;
517 }
518 }
519
520 // Combined cluster / track finder
521 if (zEnd > 999 && iFB < 0 && fgTrackReconstructor->GetTrackMethod() == 3) {
522 currIndx = -2;
523 AliMUONHitForRec hitTmp;
524 for (iz = iz0 - iFB; iz < fgCombi->Nz(); iz++) {
525 if (TMath::Abs(fgCombi->Z(iz)-fPosition) < 0.5) continue;
526 Int_t *pDEatZ = fgCombi->DEatZ(iz);
527 //cout << iz << " " << fgCombi->Z(iz) << endl;
528 zEnd = fgCombi->Z(iz);
529 iz0 = iz;
530 AliMUONDetElement *detElem = fgCombi->DetElem(pDEatZ[0]);
531 hitAdd = &hitTmp;
532 hitAdd->SetChamberNumber(detElem->Chamber());
533 //hitAdd = (AliMUONHitForRec*) detElem->HitsForRec()->First();
534 if (hitAdd) break;
535 }
536 }
537 }
538 if (zEnd>999 && ichamb==ichamEnd) endOfProp = 1; // end-of-propagation
539 else {
540 // Check if there is a chamber without hits
541 if (zEnd>999 || TMath::Abs(hitAdd->GetChamberNumber()-ichamb) > 1) {
542 if (!Back && zEnd<999) currIndx -= iFB;
543 ichamb += iFB;
544 zEnd = AliMUONConstants::DefaultChamberZ(ichamb);
545 miss = kTRUE;
546 } else {
547 ichamb = hitAdd->GetChamberNumber();
548 miss = kFALSE;
549 }
550 }
551 if (ichamb<iMin || ichamb>iMax) break;
552 // Check for missing station
553 if (!Back) {
554 dChamb = TMath::Abs(ichamb-ichambOK);
555 //cout << dChamb << " " << ichambOK << " " << fgNOfPoints << endl;
556 Int_t dStatMiss = TMath::Abs (ichamb/2 - ichambOK/2);
557 if (zEnd > 999) dStatMiss++;
558 if (dStatMiss > 1) {
559 //if (dStatMiss == 2 && ichambOK/2 != 3 || dStatMiss > 2) { // AZ - missing st. 3
560 // missing station - abandon track
561 //cout << dChamb << " " << ichambOK << " " << fgNOfPoints << " " << 1/(*fTrackPar)(4,0) << endl;
562 if (fgDebug >= 10) {
563 for (Int_t i1=0; i1<fgNOfPoints; i1++) {
564 cout << " Hit #" << ((AliMUONHitForRec*)((*fgHitForRec)[i1]))->GetChamberNumber() << " ";
565 cout << ((AliMUONHitForRec*)((*fgHitForRec)[i1]))->GetBendingCoor() << " ";
566 cout << ((AliMUONHitForRec*)((*fgHitForRec)[i1]))->GetNonBendingCoor() << " ";
567 cout << ((AliMUONHitForRec*)((*fgHitForRec)[i1]))->GetZ() << " " << " ";
568 cout << ((AliMUONHitForRec*)((*fgHitForRec)[i1]))->GetTTRTrack() << endl;
569 }
570 cout << endl;
571 cout << fNmbTrackHits << endl;
572 for (Int_t i1=0; i1<fNmbTrackHits; i1++) {
573 hit = (AliMUONHitForRec*) ((*fTrackHits)[i1]);
574 printf(" * %d %10.4f %10.4f %10.4f",
575 hit->GetChamberNumber(), hit->GetBendingCoor(),
576 hit->GetNonBendingCoor(), hit->GetZ());
577 // from raw clusters
578 rawclusters = fgTrackReconstructor->GetMUONData()->RawClusters(hit->GetChamberNumber());
579 clus = (AliMUONRawCluster*) rawclusters->UncheckedAt(hit->GetHitNumber());
580 printf("%3d", clus->GetTrack(1)-1);
581 if (clus->GetTrack(2) != 0)
582 printf("%3d \n", clus->GetTrack(2)-1);
583 else
584 printf("\n");
585
586 }
587 } // if (fgDebug >= 10)
588 if (fNmbTrackHits>2 && fRecover==0) Recover(); // try to recover track later
589 return kFALSE;
590 } // if (dStatMiss > 1)
591 } // if (!Back)
592 if (endOfProp != 0) break;
593
594 // propagate to the found Z
595
596 // Check if track steps into dipole
597 //AZ(z->-z) if (fPosition>zDipole2 && zEnd<zDipole2) {
598 if (fPosition<zDipole2 && zEnd>zDipole2) {
599 //LinearPropagation(zDipole2-zBeg);
600 ParPropagation(zDipole2);
601 MSThin(1); // multiple scattering in the chamber
602 WeightPropagation(zDipole2, kTRUE); // propagate weight matrix
603 fPosition = fPositionNew;
604 *fTrackPar = *fTrackParNew;
605 //MagnetPropagation(zEnd);
606 ParPropagation(zEnd);
607 WeightPropagation(zEnd, kTRUE);
608 fPosition = fPositionNew;
609 }
610 // Check if track steps out of dipole
611 //AZ(z->-z) else if (fPosition>zDipole1 && zEnd<zDipole1) {
612 else if (fPosition<zDipole1 && zEnd>zDipole1) {
613 //MagnetPropagation(zDipole1-zBeg);
614 ParPropagation(zDipole1);
615 MSThin(1); // multiple scattering in the chamber
616 WeightPropagation(zDipole1, kTRUE);
617 fPosition = fPositionNew;
618 *fTrackPar = *fTrackParNew;
619 //LinearPropagation(zEnd-zDipole1);
620 ParPropagation(zEnd);
621 WeightPropagation(zEnd, kTRUE);
622 fPosition = fPositionNew;
623 } else {
624 ParPropagation(zEnd);
625 //MSThin(1); // multiple scattering in the chamber
626 if (TMath::Abs(zEnd-fPosition) > 5) MSThin(1); // multiple scattering in the chamber
627 WeightPropagation(zEnd, kTRUE);
628 fPosition = fPositionNew;
629 }
630
631 // Add measurement
632 if (fRecover != 0 && hitAdd == fSkipHit && !miss) {
633 // recovered track - remove the hit
634 miss = kTRUE;
635 ichamb = fSkipHit->GetChamberNumber();
636 // remove the skipped hit
637 fTrackHits->Remove(fSkipHit);
638 fNmbTrackHits --;
639 //AZ fSkipHit->SetNTrackHits(fSkipHit->GetNTrackHits()-1); // unmark hit
640 Back = kFALSE;
641 fRecover = 0;
642 ichambOK = ((AliMUONHitForRec*)((*fTrackHits)[fNmbTrackHits-1]))->GetChamberNumber();
643 //? if (fgTrackReconstructor->GetTrackMethod() != 3) currIndx = fgHitForRec->IndexOf(fSkipHit);
644 currIndx = fgHitForRec->IndexOf(fSkipHit);
645 }
646
647 if (Back && !miss) {
648 // backward propagator
649 if (currIndx) {
650 TMatrixD pointWeight(fgkSize,fgkSize);
651 TMatrixD point(fgkSize,1);
652 TMatrixD trackParTmp = point;
653 point(0,0) = hitAdd->GetBendingCoor();
654 point(1,0) = hitAdd->GetNonBendingCoor();
655 pointWeight(0,0) = 1/hitAdd->GetBendingReso2();
656 pointWeight(1,1) = 1/hitAdd->GetNonBendingReso2();
657 TryPoint(point,pointWeight,trackParTmp,dChi2);
658 *fTrackPar = trackParTmp;
659 *fWeight += pointWeight;
660 if (fTrackDir > 0) AddMatrices (this, dChi2, hitAdd);
661 fChi2 += dChi2; // Chi2
662 } else *fTrackPar = *fTrackParNew; // adjust end point to the last hit
663 if (ichamb==ichamEnd) break;
664 currIndx ++;
665 } else {
666 // forward propagator
667 if (miss || !FindPoint(ichamb, zEnd, currIndx, iFB, hitAdd, iz)) {
668 // missing point
669 *fTrackPar = *fTrackParNew;
670 } else {
671 //add point
672 fTrackHits->Add(hitAdd); // add hit
673 fNmbTrackHits ++;
674 hitAdd->SetNTrackHits(hitAdd->GetNTrackHits()+1); // mark hit as being on track
675 ichambOK = ichamb;
676 currIndx = fgHitForRec->IndexOf(hitAdd); // Check
677 }
678 }
679 } // while
680 if (fgDebug > 0) cout << fNmbTrackHits << " " << fChi2 << " " << 1/(*fTrackPar)(4,0) << " " << fPosition << endl;
681 if (1/TMath::Abs((*fTrackPar)(4,0)) < fgTrackReconstructor->GetMinBendingMomentum()) success = kFALSE; // p < p_min
682 if (GetRecover() < 0) success = kFALSE;
683 return success;
684}
685
686 //__________________________________________________________________________
687void AliMUONTrackK::ParPropagation(Double_t zEnd)
688{
689/// Propagation of track parameters to zEnd
690 Int_t iFB, nTries;
691 Double_t dZ, step, distance, charge;
692 Double_t vGeant3[7], vGeant3New[7];
693 AliMUONTrackParam trackParam;
694
695 nTries = 0;
696 // First step using linear extrapolation
697 dZ = zEnd - fPosition;
698 fPositionNew = fPosition;
699 *fTrackParNew = *fTrackPar;
700 if (dZ == 0) return; //AZ ???
701 iFB = (Int_t)TMath::Sign(Double_t(1.0),dZ);
702 step = dZ/TMath::Cos((*fTrackPar)(2,0))/TMath::Cos((*fTrackPar)(3,0)); // linear estimate
703 //AZ(z->-z) charge = iFB*TMath::Sign(Double_t(1.0),(*fTrackPar)(4,0));
704 charge = -iFB*TMath::Sign(Double_t(1.0),(*fTrackPar)(4,0));
705 SetGeantParam(vGeant3,iFB);
706 //fTrackParNew->Print();
707
708 // Check if overstep
709 do {
710 step = TMath::Abs(step);
711 // Propagate parameters
712 AliMUONTrackExtrap::ExtrapOneStepRungekutta(charge,step,vGeant3,vGeant3New);
713 //extrap_onestep_rungekutta(charge,step,vGeant3,vGeant3New);
714 distance = zEnd - vGeant3New[2];
715 step *= dZ/(vGeant3New[2]-fPositionNew);
716 nTries ++;
717 } while (distance*iFB < 0 && TMath::Abs(distance) > fgkEpsilon);
718
719 GetFromGeantParam(vGeant3New,iFB);
720 //fTrackParNew->Print();
721
722 // Position adjustment (until within tolerance)
723 while (TMath::Abs(distance) > fgkEpsilon) {
724 dZ = zEnd - fPositionNew;
725 iFB = (Int_t)TMath::Sign(Double_t(1.0),dZ);
726 step = dZ/TMath::Cos((*fTrackParNew)(2,0))/TMath::Cos((*fTrackParNew)(3,0));
727 step = TMath::Abs(step);
728 SetGeantParam(vGeant3,iFB);
729 do {
730 // binary search
731 // Propagate parameters
732 AliMUONTrackExtrap::ExtrapOneStepRungekutta(charge,step,vGeant3,vGeant3New);
733 //extrap_onestep_rungekutta(charge,step,vGeant3,vGeant3New);
734 distance = zEnd - vGeant3New[2];
735 step /= 2;
736 nTries ++;
737 if (nTries > fgkTriesMax) AliError(Form(" Too many tries: %d", nTries));
738 } while (distance*iFB < 0);
739
740 GetFromGeantParam(vGeant3New,iFB);
741 }
742 //cout << nTries << endl;
743 //fTrackParNew->Print();
744 return;
745}
746
747 //__________________________________________________________________________
748void AliMUONTrackK::WeightPropagation(Double_t zEnd, Bool_t smooth)
749{
750/// Propagation of the weight matrix
751/// W = DtWD, where D is Jacobian
752 Int_t i, j;
753 Double_t dPar;
754
755 TMatrixD jacob(fgkSize,fgkSize);
756 jacob = 0;
757
758 // Save initial and propagated parameters
759 TMatrixD trackPar0 = *fTrackPar;
760 TMatrixD trackParNew0 = *fTrackParNew;
761
762 // Get covariance matrix
763 *fCovariance = *fWeight;
764 // check whether the Invert method returns flag if matrix cannot be inverted,
765 // and do not calculate the Determinant in that case !!!!
766 if (fCovariance->Determinant() != 0) {
767 Int_t ifail;
768 mnvertLocal(&((*fCovariance)(0,0)), fgkSize,fgkSize,fgkSize,ifail);
769 //fCovariance->Print();
770 } else {
771 AliWarning(" Determinant fCovariance=0:");
772 }
773
774 // Loop over parameters to find change of the propagated vs initial ones
775 for (i=0; i<fgkSize; i++) {
776 dPar = TMath::Sqrt((*fCovariance)(i,i));
777 *fTrackPar = trackPar0;
778 if (i == 4) dPar *= TMath::Sign(1.,-trackPar0(4,0)); // 1/p
779 (*fTrackPar)(i,0) += dPar;
780 ParPropagation(zEnd);
781 for (j=0; j<fgkSize; j++) {
782 jacob(j,i) = ((*fTrackParNew)(j,0)-trackParNew0(j,0))/dPar;
783 }
784 }
785
786 //trackParNew0.Print();
787 //TMatrixD par1(jacob,TMatrixD::kMult,trackPar0); //
788 //par1.Print();
789 TMatrixD jacob0 = jacob;
790 if (jacob.Determinant() != 0) {
791 jacob.Invert();
792 } else {
793 AliWarning(" Determinant jacob=0:");
794 }
795 TMatrixD weight1(*fWeight,TMatrixD::kMult,jacob); // WD
796 *fWeight = TMatrixD(jacob,TMatrixD::kTransposeMult,weight1); // DtWD
797
798 // Restore initial and propagated parameters
799 *fTrackPar = trackPar0;
800 *fTrackParNew = trackParNew0;
801
802 // Save for smoother
803 if (!smooth) return; // do not use smoother
804 if (fTrackDir < 0) return; // only for propagation towards int. point
805 TMatrixD *tmp = new TMatrixD(*fTrackParNew); // extrapolated parameters
806 fParExtrap->Add(tmp);
807 tmp->SetUniqueID(1);
808 tmp = new TMatrixD(*fTrackParNew); // filtered parameters (if no measurement will be found)
809 fParFilter->Add(tmp);
810 tmp->SetUniqueID(1);
811 *fCovariance = *fWeight;
812 if (fCovariance->Determinant() != 0) {
813 Int_t ifail;
814 mnvertLocal(&((*fCovariance)(0,0)), fgkSize,fgkSize,fgkSize,ifail);
815 } else {
816 AliWarning(" Determinant fCovariance=0:");
817 }
818 tmp = new TMatrixD(*fCovariance); // extrapolated covariance
819 fCovExtrap->Add(tmp);
820 tmp->SetUniqueID(1);
821 tmp = new TMatrixD(*fCovariance); // filtered covariance (if no measurement will be found)
822 fCovFilter->Add(tmp);
823 tmp->SetUniqueID(1);
824 tmp = new TMatrixD(jacob0); // Jacobian
825 fJacob->Add(tmp);
826 tmp->SetUniqueID(1);
827 if (fSteps->fN <= fNSteps) fSteps->Set(fSteps->fN+10);
828 fSteps->AddAt(fPositionNew,fNSteps++);
829 if (fgDebug > 0) cout << " WeightPropagation " << fNSteps << " " << fPositionNew << endl;
830 return;
831}
832
833 //__________________________________________________________________________
834Bool_t AliMUONTrackK::FindPoint(Int_t ichamb, Double_t zEnd, Int_t currIndx, Int_t iFB, AliMUONHitForRec *&hitAdd, Int_t iz)
835{
836/// Picks up point within a window for the chamber No ichamb
837/// Split the track if there are more than 1 hit
838 Int_t ihit, nRecTracks;
839 Double_t windowB, windowNonB, dChi2Tmp=0, dChi2, y, x, savePosition=0;
840 TClonesArray *trackPtr;
841 AliMUONHitForRec *hit, *hitLoop;
842 AliMUONTrackK *trackK;
843 AliMUONDetElement *detElem = NULL;
844
845 Bool_t ok = kFALSE;
846
847 if (fgTrackReconstructor->GetTrackMethod() == 3 && iz >= 0) {
848 // Combined cluster / track finder
849 // Check if extrapolated track passes thru det. elems. at iz
850 Int_t *pDEatZ = fgCombi->DEatZ(iz);
851 Int_t nDetElem = pDEatZ[-1];
852 //cout << fgCombi->Z(iz) << " " << nDetElem << endl;
853 for (Int_t i = 0; i < nDetElem; i++) {
854 detElem = fgCombi->DetElem(pDEatZ[i]);
855 if (detElem->Inside((*fTrackParNew)(1,0), (*fTrackParNew)(0,0), fPosition)) {
856 detElem->ClusterReco((*fTrackParNew)(1,0), (*fTrackParNew)(0,0));
857 hitAdd = (AliMUONHitForRec*) detElem->HitsForRec()->First();
858 ok = kTRUE;
859 break;
860 }
861 }
862 if (!ok) return ok; // outside det. elem.
863 ok = kFALSE;
864 }
865
866 //sigmaB = fgTrackReconstructor->GetBendingResolution(); // bending resolution
867 //sigmaNonB = fgTrackReconstructor->GetNonBendingResolution(); // non-bending resolution
868 *fCovariance = *fWeight;
869 // check whether the Invert method returns flag if matrix cannot be inverted,
870 // and do not calculate the Determinant in that case !!!!
871 if (fCovariance->Determinant() != 0) {
872 Int_t ifail;
873 mnvertLocal(&((*fCovariance)(0,0)), fgkSize,fgkSize,fgkSize,ifail);
874 } else {
875 AliWarning(" Determinant fCovariance=0:");
876 }
877 //windowB = fgkNSigma*TMath::Sqrt((*fCovariance)(0,0)+sigmaB*sigmaB);
878 //windowNonB = fgkNSigma*TMath::Sqrt((*fCovariance)(1,1)+sigmaNonB*sigmaNonB);
879 // Loop over all hits and take hits from the chamber
880 TMatrixD pointWeight(fgkSize,fgkSize);
881 TMatrixD saveWeight = pointWeight;
882 TMatrixD pointWeightTmp = pointWeight;
883 TMatrixD point(fgkSize,1);
884 TMatrixD trackPar = point;
885 TMatrixD trackParTmp = point;
886 Int_t nHitsOK = 0, ihitB = currIndx, ihitE = fgNOfPoints, iDhit = iFB;
887 Double_t zLast;
888 zLast = ((AliMUONHitForRec*)fTrackHits->Last())->GetZ();
889 if (fgTrackReconstructor->GetTrackMethod() == 3 && detElem) {
890 ihitB = 0;
891 ihitE = detElem->NHitsForRec();
892 iDhit = 1;
893 }
894
895 TArrayD branchChi2(20);
896 for (ihit = ihitB; ihit >= 0 && ihit < ihitE; ihit+=iDhit) {
897 if (fgTrackReconstructor->GetTrackMethod() != 3 || !detElem)
898 hit = (AliMUONHitForRec*) ((*fgHitForRec)[ihit]);
899 else hit = (AliMUONHitForRec*) detElem->HitsForRec()->UncheckedAt(ihit);
900 if (hit->GetChamberNumber() == ichamb) {
901 //if (TMath::Abs(hit->GetZ()-zEnd) < 0.1) {
902 if (TMath::Abs(hit->GetZ()-zEnd) < 0.5) {
903 y = hit->GetBendingCoor();
904 x = hit->GetNonBendingCoor();
905 if (hit->GetBendingReso2() < 0) {
906 // Combined cluster / track finder
907 hit->SetBendingReso2(fgTrackReconstructor->GetBendingResolution()*
908 fgTrackReconstructor->GetBendingResolution());
909 hit->SetNonBendingReso2(fgTrackReconstructor->GetNonBendingResolution()*
910 fgTrackReconstructor->GetNonBendingResolution());
911 }
912 windowB = fgkNSigma*TMath::Sqrt((*fCovariance)(0,0)+hit->GetBendingReso2());
913 windowNonB = fgkNSigma*TMath::Sqrt((*fCovariance)(1,1)+hit->GetNonBendingReso2());
914
915 // windowB = TMath::Min (windowB,5.);
916 if (fgkNSigma > 6) windowB = TMath::Min (windowB,5.);
917 /*
918 if (TMath::Abs(hit->GetZ()-zLast) < 50) {
919 windowB = TMath::Min (windowB,0.5);
920 windowNonB = TMath::Min (windowNonB,3.);
921 } else if (TMath::Abs(hit->GetZ()-zLast) < 200) {
922 windowB = TMath::Min (windowB,1.5);
923 windowNonB = TMath::Min (windowNonB,3.);
924 } else if (TMath::Abs(hit->GetZ()-zLast) < 350) {
925 windowB = TMath::Min (windowB,4.);
926 windowNonB = TMath::Min (windowNonB,6.);
927 }
928 */
929
930 // Test
931 /*
932 if (TMath::Abs(hit->GetZ()-zLast) < 50) {
933 windowB = 5*TMath::Sqrt((*fCovariance)(0,0)+hit->GetBendingReso2());
934 windowNonB = 5*TMath::Sqrt((*fCovariance)(1,1)+hit->GetNonBendingReso2());
935 } else {
936 windowB = fgkNSigma*TMath::Sqrt((*fCovariance)(0,0)+hit->GetBendingReso2());
937 windowNonB = fgkNSigma*TMath::Sqrt((*fCovariance)(1,1)+hit->GetNonBendingReso2());
938 }
939 */
940
941 //cout << TMath::Abs((*fTrackParNew)(0,0)-y) << " " << windowB << " " << TMath::Abs((*fTrackParNew)(1,0)-x) << " " << windowNonB << " " << zEnd << endl;
942 if (TMath::Abs((*fTrackParNew)(0,0)-y) <= windowB &&
943 TMath::Abs((*fTrackParNew)(1,0)-x) <= windowNonB) {
944 //if (TMath::Abs((*fTrackParNew)(0,0)-y) <= windowB &&
945 // TMath::Abs((*fTrackParNew)(1,0)-x) <= windowNonB &&
946 // hit->GetTrackRefSignal() == 1) { // just for test
947 // Vector of measurements and covariance matrix
948 //fprintf(lun1,"%3d %3d %10.4f %10.4f \n", gAlice->GetEvNumber(), ichamb, x, y);
949 if (TMath::Abs(hit->GetZ()-zEnd) > 0.05) {
950 // Adjust position: for multiple hits in the chamber or misalignment (Z as a function of X or Y)
951 //AliWarning(Form(" *** adjust %f %f ", zEnd, hit->GetZ()));
952 zEnd = hit->GetZ();
953 *fTrackPar = *fTrackParNew;
954 ParPropagation(zEnd);
955 WeightPropagation(zEnd, kTRUE);
956 fPosition = fPositionNew;
957 *fTrackPar = *fTrackParNew;
958 // Get covariance
959 *fCovariance = *fWeight;
960 if (fCovariance->Determinant() != 0) {
961 Int_t ifail;
962 mnvertLocal(&((*fCovariance)(0,0)), fgkSize,fgkSize,fgkSize,ifail);
963 } else {
964 AliWarning(" Determinant fCovariance=0:" );
965 }
966 }
967 point.Zero();
968 point(0,0) = y;
969 point(1,0) = x;
970 pointWeight(0,0) = 1/hit->GetBendingReso2();
971 pointWeight(1,1) = 1/hit->GetNonBendingReso2();
972 TryPoint(point,pointWeight,trackPar,dChi2);
973 if (TMath::Abs(1./(trackPar)(4,0)) < fgTrackReconstructor->GetMinBendingMomentum()) continue; // p < p_min - next hit
974 // if (TMath::Sign(1.,(trackPar)(4,0)*(*fTrackPar)(4,0)) < 0) continue; // change of sign
975 ok = kTRUE;
976 nHitsOK++;
977 //if (nHitsOK > -1) {
978 if (nHitsOK == 1) {
979 // Save current members
980 saveWeight = *fWeight;
981 savePosition = fPosition;
982 // temporary storage for the current track
983 dChi2Tmp = dChi2;
984 trackParTmp = trackPar;
985 pointWeightTmp = pointWeight;
986 hitAdd = hit;
987 if (fgDebug > 0) cout << " Added point: " << x << " " << y << " " << dChi2 << endl;
988 branchChi2[0] = dChi2;
989 } else {
990 // branching: create a new track
991 trackPtr = fgTrackReconstructor->GetRecTracksPtr();
992 nRecTracks = fgTrackReconstructor->GetNRecTracks();
993 trackK = new ((*trackPtr)[nRecTracks]) AliMUONTrackK(NULL, NULL);
994 *trackK = *this;
995 fgTrackReconstructor->SetNRecTracks(nRecTracks+1);
996 if (fgDebug > 0) cout << " ******** New track: " << ichamb << " " << hit->GetTTRTrack() << " " << 1/(trackPar)(4,0) << " " << hit->GetBendingCoor() << " " << hit->GetNonBendingCoor() << " " << fNmbTrackHits << " " << nRecTracks << endl;
997 trackK->fRecover = 0;
998 *(trackK->fTrackPar) = trackPar;
999 *(trackK->fWeight) += pointWeight;
1000 trackK->fChi2 += dChi2;
1001 // Check
1002 /*
1003 *(trackK->fCovariance) = *(trackK->fWeight);
1004 if (trackK->fCovariance->Determinant() != 0) {
1005 Int_t ifail;
1006 mnvertLocal(&((*(trackK->fCovariance))(0,0)), fgkSize,fgkSize,fgkSize,ifail);
1007 }
1008 cout << (*(trackK->fCovariance))(0,0) << " " << (*(trackK->fCovariance))(1,1) << " " << (*fCovariance)(0,0) << " " << (*fCovariance)(1,1) << endl;
1009 */
1010 // Add filtered matrices for the smoother
1011 if (fTrackDir > 0) {
1012 if (nHitsOK > 2) { // check for position adjustment
1013 for (Int_t i=trackK->fNSteps-1; i>=0; i--) {
1014 if (TMath::Abs(hit->GetZ()-(*trackK->fSteps)[i]) > 0.1) {
1015 //if (TMath::Abs(hit->GetZ()-(trackK->fSteps)[i]) > 0.1) {
1016 RemoveMatrices(trackK);
1017 AliError(Form(" *** Position adjustment 1: %f %f", hit->GetZ(), (*trackK->fSteps)[i]));
1018 }
1019 else break;
1020 }
1021 }
1022 AddMatrices (trackK, dChi2, hit);
1023 }
1024 // Mark hits as being on 2 tracks
1025 for (Int_t i=0; i<fNmbTrackHits; i++) {
1026 hitLoop = (AliMUONHitForRec*) ((*fTrackHits)[i]);
1027 //AZ hitLoop->SetNTrackHits(hitLoop->GetNTrackHits()+1);
1028 if (fgDebug >=10) {
1029 cout << " ** ";
1030 cout << hitLoop->GetChamberNumber() << " ";
1031 cout << hitLoop->GetBendingCoor() << " ";
1032 cout << hitLoop->GetNonBendingCoor() << " ";
1033 cout << hitLoop->GetZ() << " " << " ";
1034 cout << hitLoop->GetTTRTrack() << endl;
1035 printf(" ** %d %10.4f %10.4f %10.4f\n",
1036 hitLoop->GetChamberNumber(), hitLoop->GetBendingCoor(),
1037 hitLoop->GetNonBendingCoor(), hitLoop->GetZ());
1038 }
1039 }
1040 //add point
1041 trackK->fTrackHits->Add(hit); // add hit
1042 trackK->fNmbTrackHits ++;
1043 hit->SetNTrackHits(hit->GetNTrackHits()+1); // mark hit as being on track
1044 if (ichamb == 9) {
1045 // the last chamber
1046 trackK->fTrackDir = 1;
1047 trackK->fBPFlag = kTRUE;
1048 }
1049 if (nHitsOK > branchChi2.GetSize()) branchChi2.Set(branchChi2.GetSize()+10);
1050 branchChi2[nHitsOK-1] = dChi2;
1051 }
1052 }
1053 }
1054 } else break; // different chamber
1055 } // for (ihit=currIndx;
1056 if (ok) {
1057 // Restore members
1058 *fTrackPar = trackParTmp;
1059 *fWeight = saveWeight;
1060 *fWeight += pointWeightTmp;
1061 fChi2 += dChi2Tmp; // Chi2
1062 fPosition = savePosition;
1063 // Add filtered matrices for the smoother
1064 if (fTrackDir > 0) {
1065 for (Int_t i=fNSteps-1; i>=0; i--) {
1066 if (TMath::Abs(fPosition-(*fSteps)[i]) > 0.1) {
1067 //if (TMath::Abs(fPosition-fSteps[i]) > 0.1) {
1068 RemoveMatrices(this);
1069 if (fgDebug > 0) cout << " *** Position adjustment 2 " << fPosition << " " << (*fSteps)[i] << endl;
1070 }
1071 else break;
1072 } // for (Int_t i=fNSteps-1;
1073 AddMatrices (this, dChi2Tmp, hitAdd);
1074 /*
1075 if (nHitsOK > 1) {
1076 for (Int_t i=0; i<fNSteps; i++) cout << (*fSteps)[i] << " "; cout << endl;
1077 for (Int_t i=0; i<trackK->fNSteps; i++) cout << (*trackK->fSteps)[i] << " "; cout << endl;
1078 }
1079 */
1080 } // if (fTrackDir > 0)
1081 // Check for maximum number of branches - exclude excessive
1082 if (nHitsOK > 1) CheckBranches(branchChi2, nHitsOK);
1083 }
1084 return ok;
1085}
1086
1087 //__________________________________________________________________________
1088void AliMUONTrackK::CheckBranches(TArrayD &branchChi2, Int_t nBranch)
1089{
1090/// Check for maximum number of branches - exclude excessive
1091
1092 Int_t nBranchMax = 5;
1093 if (nBranch <= nBranchMax) return;
1094
1095 Double_t *chi2 = branchChi2.GetArray();
1096 Int_t *indx = new Int_t [nBranch];
1097 TMath::Sort (nBranch, chi2, indx, kFALSE);
1098 TClonesArray *trackPtr = fgTrackReconstructor->GetRecTracksPtr();
1099 Int_t nRecTracks = fgTrackReconstructor->GetNRecTracks();
1100 Int_t ibeg = nRecTracks - nBranch;
1101
1102 // Discard excessive branches with higher Chi2 contribution
1103 for (Int_t i = nBranchMax; i < nBranch; ++i) {
1104 if (indx[i] == 0) {
1105 // Discard current track
1106 SetRecover(-1);
1107 continue;
1108 }
1109 Int_t j = ibeg + indx[i];
1110 AliMUONTrackK *trackK = (AliMUONTrackK*) trackPtr->UncheckedAt(j);
1111 trackK->SetRecover(-1);
1112 }
1113 delete [] indx;
1114}
1115
1116 //__________________________________________________________________________
1117void AliMUONTrackK::TryPoint(TMatrixD &point, const TMatrixD &pointWeight, TMatrixD &trackParTmp, Double_t &dChi2)
1118{
1119/// Adds a measurement point (modifies track parameters and computes
1120/// change of Chi2)
1121
1122 // Solving linear system (W+U)p' = U(m-p) + (W+U)p
1123 TMatrixD wu = *fWeight;
1124 wu += pointWeight; // W+U
1125 trackParTmp = point;
1126 trackParTmp -= *fTrackParNew; // m-p
1127 TMatrixD right(pointWeight,TMatrixD::kMult,trackParTmp); // U(m-p)
1128 TMatrixD right1(wu,TMatrixD::kMult,*fTrackParNew); // (W+U)p
1129 right += right1; // U(m-p) + (W+U)p
1130
1131 // check whether the Invert method returns flag if matrix cannot be inverted,
1132 // and do not calculate the Determinant in that case !!!!
1133 if (wu.Determinant() != 0) {
1134 Int_t ifail;
1135 mnvertLocal(&((wu)(0,0)), fgkSize,fgkSize,fgkSize,ifail);
1136 } else {
1137 AliWarning(" Determinant wu=0:");
1138 }
1139 trackParTmp = TMatrixD(wu,TMatrixD::kMult,right);
1140
1141 right1 = trackParTmp;
1142 right1 -= point; // p'-m
1143 point = trackParTmp;
1144 point -= *fTrackParNew; // p'-p
1145 right = TMatrixD(*fWeight,TMatrixD::kMult,point); // W(p'-p)
1146 TMatrixD value(point,TMatrixD::kTransposeMult,right); // (p'-p)'W(p'-p)
1147 dChi2 = value(0,0);
1148 right = TMatrixD(pointWeight,TMatrixD::kMult,right1); // U(p'-m)
1149 value = TMatrixD(right1,TMatrixD::kTransposeMult,right); // (p'-m)'U(p'-m)
1150 dChi2 += value(0,0);
1151 return;
1152}
1153
1154 //__________________________________________________________________________
1155void AliMUONTrackK::MSThin(Int_t sign)
1156{
1157/// Adds multiple scattering in a thin layer (only angles are affected)
1158 Double_t cosAlph, cosBeta, momentum, velo, path, theta0;
1159
1160 // check whether the Invert method returns flag if matrix cannot be inverted,
1161 // and do not calculate the Determinant in that case !!!!
1162 if (fWeight->Determinant() != 0) {
1163 Int_t ifail;
1164 mnvertLocal(&((*fWeight)(0,0)), fgkSize,fgkSize,fgkSize,ifail);
1165 } else {
1166 AliWarning(" Determinant fWeight=0:");
1167 }
1168
1169 cosAlph = TMath::Cos((*fTrackParNew)(2,0));
1170 cosBeta = TMath::Cos((*fTrackParNew)(3,0));
1171 momentum = 1/(*fTrackParNew)(4,0); // particle momentum
1172 //velo = momentum/TMath::Sqrt(momentum*momentum+muonMass*muonMass); // velocity/c for muon hypothesis
1173 velo = 1; // relativistic
1174 path = TMath::Abs(AliMUONConstants::ChamberThicknessInX0()/cosAlph/cosBeta); // path length
1175 theta0 = 0.0136/velo/momentum*TMath::Sqrt(path)*(1+0.038*TMath::Log(path)); // projected scattering angle
1176
1177 (*fWeight)(2,2) += sign*theta0/cosBeta*theta0/cosBeta; // alpha
1178 (*fWeight)(3,3) += sign*theta0*theta0; // beta
1179 Int_t ifail;
1180 mnvertLocal(&((*fWeight)(0,0)), fgkSize,fgkSize,fgkSize,ifail);
1181 return;
1182}
1183
1184 //__________________________________________________________________________
1185void AliMUONTrackK::StartBack(void)
1186{
1187/// Starts backpropagator
1188
1189 fBPFlag = kTRUE;
1190 fChi2 = 0;
1191 for (Int_t i=0; i<fgkSize; i++) {
1192 for (Int_t j=0; j<fgkSize; j++) {
1193 if (j==i) (*fWeight)(i,i) /= 100;
1194 //if (j==i) (*fWeight)(i,i) /= fNmbTrackHits*fNmbTrackHits;
1195 else (*fWeight)(j,i) = 0;
1196 }
1197 }
1198 // Sort hits on track in descending order in abs(z)
1199 SortHits(0, fTrackHits);
1200}
1201
1202 //__________________________________________________________________________
1203void AliMUONTrackK::SortHits(Int_t iflag, TObjArray *array)
1204{
1205/// Sort hits in Z if the seed segment is in the last but one station
1206/// (if iflag==0 in descending order in abs(z), if !=0 - unsort)
1207
1208 if (iflag && ((AliMUONHitForRec*)(array->UncheckedAt(0)))->GetChamberNumber() == 6) return;
1209 Double_t z = 0, zmax = TMath::Abs(((AliMUONHitForRec*)(array->UncheckedAt(0)))->GetZ());
1210 Int_t i = 1, entries = array->GetEntriesFast();
1211 for ( ; i<entries; i++) {
1212 if (iflag) {
1213 if (((AliMUONHitForRec*)(array->UncheckedAt(i)))->GetChamberNumber() == 6) break;
1214 } else {
1215 z = TMath::Abs(((AliMUONHitForRec*)(array->UncheckedAt(i)))->GetZ());
1216 if (z < zmax) break;
1217 zmax = z;
1218 if (fgDebug >= 10) cout << " " << zmax << " " << z << " " << i << endl;
1219 }
1220 }
1221 if (!iflag) i--;
1222 for (Int_t j=0; j<=(i-1)/2; j++) {
1223 TObject *hit = array->UncheckedAt(j);
1224 array->AddAt(array->UncheckedAt(i-j),j);
1225 array->AddAt(hit,i-j);
1226 }
1227 if (fgDebug >= 10) {
1228 for (i=0; i<entries; i++)
1229 cout << ((AliMUONHitForRec*)(array->UncheckedAt(i)))->GetZ() << " ";
1230 cout << " - Sort" << endl;
1231 }
1232}
1233
1234 //__________________________________________________________________________
1235void AliMUONTrackK::SetGeantParam(Double_t *VGeant3, Int_t iFB)
1236{
1237/// Set vector of Geant3 parameters pointed to by "VGeant3"
1238/// from track parameters
1239
1240 VGeant3[0] = (*fTrackParNew)(1,0); // X
1241 VGeant3[1] = (*fTrackParNew)(0,0); // Y
1242 VGeant3[2] = fPositionNew; // Z
1243 VGeant3[3] = iFB*TMath::Sin((*fTrackParNew)(3,0)); // Px/Ptot
1244 VGeant3[4] = iFB*TMath::Cos((*fTrackParNew)(3,0))*TMath::Sin((*fTrackParNew)(2,0)); // Py/Ptot
1245 VGeant3[5] = iFB*TMath::Sqrt(1.0-VGeant3[3]*VGeant3[3]-VGeant3[4]*VGeant3[4]); // Pz/Ptot
1246 VGeant3[6] = 1/TMath::Abs((*fTrackParNew)(4,0)); // Ptot
1247}
1248
1249 //__________________________________________________________________________
1250void AliMUONTrackK::GetFromGeantParam(Double_t *VGeant3, Int_t iFB)
1251{
1252/// Get track parameters from vector of Geant3 parameters pointed
1253/// to by "VGeant3"
1254
1255 fPositionNew = VGeant3[2]; // Z
1256 (*fTrackParNew)(0,0) = VGeant3[1]; // Y
1257 (*fTrackParNew)(1,0) = VGeant3[0]; // X
1258 (*fTrackParNew)(3,0) = TMath::ASin(iFB*VGeant3[3]); // beta
1259 (*fTrackParNew)(2,0) = TMath::ASin(iFB*VGeant3[4]/TMath::Cos((*fTrackParNew)(3,0))); // alpha
1260 (*fTrackParNew)(4,0) = 1/VGeant3[6]*TMath::Sign(Double_t(1.0),(*fTrackPar)(4,0)); // 1/Ptot
1261}
1262
1263 //__________________________________________________________________________
1264void AliMUONTrackK::SetTrackQuality(Int_t iChi2)
1265{
1266/// Computes "track quality" from Chi2 (if iChi2==0) or vice versa
1267
1268 if (fChi2 > 500) {
1269 AliWarning(Form(" *** Too high Chi2: %f ", fChi2));
1270 fChi2 = 500;
1271 }
1272 if (iChi2 == 0) fChi2 = fNmbTrackHits + (500.-fChi2)/501;
1273 else fChi2 = 500 - (fChi2-fNmbTrackHits)*501;
1274}
1275
1276 //__________________________________________________________________________
1277Int_t AliMUONTrackK::Compare(const TObject* trackK) const
1278{
1279/// "Compare" function to sort with decreasing "track quality".
1280/// Returns +1 (0, -1) if quality of current track
1281/// is smaller than (equal to, larger than) quality of trackK
1282
1283 if (fChi2 < ((AliMUONTrackK*)trackK)->fChi2) return(+1);
1284 else if (fChi2 == ((AliMUONTrackK*)trackK)->fChi2) return(0);
1285 else return(-1);
1286}
1287
1288 //__________________________________________________________________________
1289Bool_t AliMUONTrackK::KeepTrack(AliMUONTrackK* track0) const
1290{
1291/// Check whether or not to keep current track
1292/// (keep, if it has less than half of common hits with track0)
1293 Int_t hitsInCommon, nHits0, i, j, nTrackHits2;
1294 AliMUONHitForRec *hit0, *hit1;
1295
1296 hitsInCommon = 0;
1297 nHits0 = track0->fNmbTrackHits;
1298 nTrackHits2 = fNmbTrackHits/2;
1299
1300 for (i=0; i<nHits0; i++) {
1301 // Check if hit belongs to several tracks
1302 hit0 = (AliMUONHitForRec*) (*track0->fTrackHits)[i];
1303 if (hit0->GetNTrackHits() == 1) continue;
1304 for (j=0; j<fNmbTrackHits; j++) {
1305 hit1 = (AliMUONHitForRec*) (*fTrackHits)[j];
1306 if (hit1->GetNTrackHits() == 1) continue;
1307 if (hit0 == hit1) {
1308 hitsInCommon++;
1309 if (hitsInCommon >= nTrackHits2) return kFALSE;
1310 break;
1311 }
1312 } // for (j=0;
1313 } // for (i=0;
1314 return kTRUE;
1315}
1316
1317 //__________________________________________________________________________
1318void AliMUONTrackK::Kill(void)
1319{
1320/// Kill track candidate
1321 fgTrackReconstructor->GetRecTracksPtr()->Remove(this);
1322}
1323
1324 //__________________________________________________________________________
1325void AliMUONTrackK::FillMUONTrack(void)
1326{
1327/// Compute track parameters at hit positions (as for AliMUONTrack)
1328
1329 // Set Chi2
1330 SetFitFMin(fChi2);
1331
1332 // Set track parameters at vertex
1333 AliMUONTrackParam trackParam;
1334 SetTrackParam(&trackParam, fTrackPar, fPosition);
1335 SetTrackParamAtVertex(&trackParam);
1336
1337 // Set track parameters at hits
1338 for (Int_t i = fNmbTrackHits-1; i>=0; i--) {
1339 if ((*fChi2Smooth)[i] < 0) {
1340 // Propagate through last chambers
1341 AliMUONTrackExtrap::ExtrapToZ(&trackParam, ((AliMUONHitForRec*)((*fTrackHits)[i]))->GetZ());
1342 } else {
1343 // Take saved info
1344 SetTrackParam(&trackParam, (TMatrixD*)fParSmooth->UncheckedAt(i), ((AliMUONHitForRec*)((*fTrackHits)[i]))->GetZ());
1345 }
1346 AddTrackParamAtHit(&trackParam,(AliMUONHitForRec*)fTrackHits->UncheckedAt(i));
1347 // Fill array of HitForRec's
1348 AddHitForRecAtHit((AliMUONHitForRec*)fTrackHits->UncheckedAt(i));
1349 }
1350}
1351
1352 //__________________________________________________________________________
1353void AliMUONTrackK::SetTrackParam(AliMUONTrackParam *trackParam, TMatrixD *par, Double_t z)
1354{
1355/// Fill AliMUONTrackParam object
1356
1357 trackParam->SetBendingCoor((*par)(0,0));
1358 trackParam->SetNonBendingCoor((*par)(1,0));
1359 trackParam->SetBendingSlope(TMath::Tan((*par)(2,0)));
1360 trackParam->SetNonBendingSlope(TMath::Tan((*par)(3,0))/TMath::Cos((*par)(2,0)));
1361 trackParam->SetInverseBendingMomentum((*par)(4,0)/TMath::Cos((*par)(3,0)));
1362 trackParam->SetZ(z);
1363}
1364
1365 //__________________________________________________________________________
1366void AliMUONTrackK::Branson(void)
1367{
1368/// Propagates track to the vertex thru absorber using Branson correction
1369/// (makes use of the AliMUONTrackParam class)
1370
1371 //AliMUONTrackParam *trackParam = new AliMUONTrackParam();
1372 AliMUONTrackParam trackParam = AliMUONTrackParam();
1373 /*
1374 trackParam->SetBendingCoor((*fTrackPar)(0,0));
1375 trackParam->SetNonBendingCoor((*fTrackPar)(1,0));
1376 trackParam->SetBendingSlope(TMath::Tan((*fTrackPar)(2,0)));
1377 trackParam->SetNonBendingSlope(TMath::Tan((*fTrackPar)(3,0))/TMath::Cos((*fTrackPar)(2,0)));
1378 trackParam->SetInverseBendingMomentum((*fTrackPar)(4,0)/TMath::Cos((*fTrackPar)(3,0)));
1379 trackParam->SetZ(fPosition);
1380 */
1381 SetTrackParam(&trackParam, fTrackPar, fPosition);
1382
1383 AliMUONTrackExtrap::ExtrapToVertex(&trackParam, Double_t(0.), Double_t(0.), Double_t(0.));
1384
1385 (*fTrackPar)(0,0) = trackParam.GetBendingCoor();
1386 (*fTrackPar)(1,0) = trackParam.GetNonBendingCoor();
1387 (*fTrackPar)(2,0) = TMath::ATan(trackParam.GetBendingSlope());
1388 (*fTrackPar)(3,0) = TMath::ATan(TMath::Cos((*fTrackPar)(2,0))*trackParam.GetNonBendingSlope());
1389 (*fTrackPar)(4,0) = TMath::Cos((*fTrackPar)(3,0))*trackParam.GetInverseBendingMomentum();
1390 fPosition = trackParam.GetZ();
1391 //delete trackParam;
1392 if (fgDebug > 0) cout << 1/(*fTrackPar)(4,0) << " " << fPosition << " " << (*fTrackPar)(0,0) << endl;
1393
1394 // Get covariance matrix
1395 *fCovariance = *fWeight;
1396 if (fCovariance->Determinant() != 0) {
1397 Int_t ifail;
1398 mnvertLocal(&((*fCovariance)(0,0)), fgkSize,fgkSize,fgkSize,ifail);
1399 } else {
1400 AliWarning(" Determinant fCovariance=0:");
1401 }
1402}
1403
1404 //__________________________________________________________________________
1405void AliMUONTrackK::GoToZ(Double_t zEnd)
1406{
1407/// Propagates track to given Z
1408
1409 ParPropagation(zEnd);
1410 MSThin(1); // multiple scattering in the chamber
1411 WeightPropagation(zEnd, kFALSE);
1412 fPosition = fPositionNew;
1413 *fTrackPar = *fTrackParNew;
1414}
1415
1416 //__________________________________________________________________________
1417void AliMUONTrackK::GoToVertex(Int_t iflag)
1418{
1419/// Version 3.08
1420/// Propagates track to the vertex
1421/// All material constants are taken from AliRoot
1422
1423 static Double_t x01[5] = { 24.282, // C
1424 24.282, // C
1425 11.274, // Concrete
1426 1.758, // Fe
1427 1.758}; // Fe (cm)
1428 // inner part theta < 3 degrees
1429 static Double_t x02[5] = { 30413, // Air
1430 24.282, // C
1431 11.274, // Concrete
1432 1.758, // Fe
1433 0.369}; // W (cm)
1434 // z positions of the materials inside the absober outer part theta > 3 degres
1435 static Double_t zPos[10] = {-90, -105, -315, -443, -468};
1436
1437 Double_t dZ, r0Norm, x0, deltaP, dChi2, pTotal, pOld;
1438 AliMUONHitForRec *hit;
1439 AliMUONRawCluster *clus;
1440 TClonesArray *rawclusters;
1441
1442 // First step to the rear end of the absorber
1443 Double_t zRear = -503;
1444 GoToZ(zRear);
1445 Double_t tan3 = TMath::Tan(3./180*TMath::Pi());
1446
1447 // Go through absorber
1448 pOld = 1/(*fTrackPar)(4,0);
1449 Double_t r0Rear = (*fTrackPar)(0,0)*(*fTrackPar)(0,0) +
1450 (*fTrackPar)(1,0)*(*fTrackPar)(1,0);
1451 r0Rear = TMath::Sqrt(r0Rear)/TMath::Abs(fPosition)/tan3;
1452 r0Norm = r0Rear;
1453 Double_t p0, cos25, cos60;
1454 if (!iflag) goto vertex;
1455
1456 for (Int_t i=4; i>=0; i--) {
1457 ParPropagation(zPos[i]);
1458 WeightPropagation(zPos[i], kFALSE);
1459 dZ = TMath::Abs (fPositionNew-fPosition);
1460 if (r0Norm > 1) x0 = x01[i];
1461 else x0 = x02[i];
1462 MSLine(dZ,x0); // multiple scattering in the medium (linear approximation)
1463 fPosition = fPositionNew;
1464 *fTrackPar = *fTrackParNew;
1465 r0Norm = (*fTrackPar)(0,0)*(*fTrackPar)(0,0) +
1466 (*fTrackPar)(1,0)*(*fTrackPar)(1,0);
1467 r0Norm = TMath::Sqrt(r0Norm)/TMath::Abs(fPosition)/tan3;
1468 }
1469 // Correct momentum for energy losses
1470 pTotal = 1/TMath::Abs((*fTrackPar)(4,0));
1471 p0 = pTotal;
1472 cos25 = TMath::Cos(2.5/180*TMath::Pi());
1473 cos60 = TMath::Cos(6.0/180*TMath::Pi());
1474 for (Int_t j=0; j<1; j++) {
1475 /*
1476 if (r0Rear > 1) {
1477 if (p0 < 20) {
1478 deltaP = 2.164 + 0.145e-1*p0 - 0.417e-3*p0*p0;
1479 } else {
1480 deltaP = 2.275 + 0.102e-2*p0 - 0.674e-6*p0*p0;
1481 }
1482 } else {
1483 if (p0 < 20) {
1484 deltaP = 2.581 + 0.188e-1*p0 - 0.398e-3*p0*p0;
1485 } else {
1486 deltaP = 2.727 + 0.356e-2*p0 + 0.242e-5*p0*p0;
1487 }
1488 }
1489 */
1490 /*
1491 if (r0Rear < 1) {
1492 //W
1493 if (p0<15) {
1494 deltaP = 2.737 + 0.0494*p0 - 0.001123*p0*p0;
1495 } else {
1496 deltaP = 3.0643 + 0.01346*p0;
1497 }
1498 deltaP *= 0.95;
1499 } else {
1500 //Pb
1501 if (p0<15) {
1502 deltaP = 2.1380 + 0.0351*p0 - 0.000853*p0*p0;
1503 } else {
1504 deltaP = 2.407 + 0.00702*p0;
1505 }
1506 deltaP *= 0.95;
1507 }
1508 */
1509 /*
1510 if (r0Rear < 1) {
1511 //W
1512 if (p0<18) {
1513 deltaP = 2.439 + 0.806e-1*p0 - 0.500e-2*p0*p0 + 0.106e-3*p0*p0*p0;
1514 } else {
1515 deltaP = 2.767 + 0.742e-2*p0 - 0.196e-4*p0*p0 + 0.403e-7*p0*p0*p0;
1516 }
1517 //deltaP += 0.2;
1518 deltaP *= cos25;
1519 } else {
1520 //Pb
1521 if (p0<18) {
1522 deltaP = 2.209 + 0.800e-2*p0;
1523 } else {
1524 deltaP = 2.285 + 0.141e-2*p0 - 0.446e-6*p0*p0;
1525 }
1526 //deltaP += 0.2;
1527 deltaP *= cos60;
1528 }
1529 deltaP *= 1.1;
1530 */
1531 //*
1532 if (r0Rear < 1) {
1533 if (p0 < 20) {
1534 deltaP = 2.5938 + 0.0570 * p0 - 0.001151 * p0 * p0;
1535 } else {
1536 deltaP = 3.0714 + 0.011767 * p0;
1537 }
1538 deltaP *= 0.75;
1539 } else {
1540 if (p0 < 20) {
1541 deltaP = 2.1207 + 0.05478 * p0 - 0.00145079 * p0 * p0;
1542 } else {
1543 deltaP = 2.6069 + 0.0051705 * p0;
1544 }
1545 deltaP *= 0.9;
1546 }
1547 //*/
1548
1549 p0 = pTotal + deltaP/TMath::Abs(TMath::Cos((*fTrackPar)(2,0))/TMath::Cos((*fTrackPar)(3,0)));
1550 }
1551 (*fTrackPar)(4,0) = 1/p0*TMath::Sign((Double_t)1.,(*fTrackPar)(4,0));
1552
1553 // Go to the vertex
1554vertex:
1555 ParPropagation((Double_t)0.);
1556 WeightPropagation((Double_t)0., kFALSE);
1557 fPosition = fPositionNew;
1558 //*fTrackPar = *fTrackParNew;
1559 // Add vertex as a hit
1560 TMatrixD pointWeight(fgkSize,fgkSize);
1561 TMatrixD point(fgkSize,1);
1562 TMatrixD trackParTmp = point;
1563 point(0,0) = 0; // vertex coordinate - should be taken somewhere
1564 point(1,0) = 0; // vertex coordinate - should be taken somewhere
1565 pointWeight(0,0) = 1/1.e-3/1.e-3; // 10 um error
1566 pointWeight(1,1) = 1/1.e-3/1.e-3; // 10 um error
1567 TryPoint(point,pointWeight,trackParTmp,dChi2);
1568 *fTrackPar = trackParTmp;
1569 *fWeight += pointWeight;
1570 fChi2 += dChi2; // Chi2
1571 if (fgDebug < 0) return; // no output
1572
1573 cout << pOld << " " << 1/(*fTrackPar)(4,0) << " " << dChi2 << " " << fChi2 << " " << fNmbTrackHits << endl;
1574 for (Int_t i1=0; i1<fNmbTrackHits; i1++) {
1575 hit = (AliMUONHitForRec*) ((*fTrackHits)[i1]);
1576 printf ("%5d", hit->GetChamberNumber());
1577 }
1578 cout << endl;
1579 if (fgDebug > 0) {
1580 for (Int_t i1=0; i1<fNmbTrackHits; i1++) {
1581 hit = (AliMUONHitForRec*) ((*fTrackHits)[i1]);
1582 //cout << ((AliMUONHitForRec*)((*fTrackHits)[i1]))->GetHitNumber() << " ";
1583 //cout << ((AliMUONHitForRec*)((*fTrackHits)[i1]))->GetZ() << " ";
1584 printf ("%5d", fgHitForRec->IndexOf(hit));
1585 }
1586 cout << endl;
1587 }
1588
1589 // from raw clusters
1590 for (Int_t i1=0; i1<fNmbTrackHits; i1++) {
1591 hit = (AliMUONHitForRec*) ((*fTrackHits)[i1]);
1592 if (hit->GetHitNumber() < 0) { // combined cluster / track finder
1593 Int_t index = -hit->GetHitNumber() / 100000;
1594 Int_t iPos = -hit->GetHitNumber() - index * 100000;
1595 clus = (AliMUONRawCluster*) fgCombi->DetElem(index-1)->RawClusters()->UncheckedAt(iPos);
1596 } else {
1597 rawclusters = fgTrackReconstructor->GetMUONData()->RawClusters(hit->GetChamberNumber());
1598 clus = (AliMUONRawCluster*) rawclusters->UncheckedAt(hit->GetHitNumber());
1599 }
1600 printf ("%5d", clus->GetTrack(1)%10000000);
1601
1602 cout << endl;
1603 for (Int_t i1=0; i1<fNmbTrackHits; i1++) {
1604 hit = (AliMUONHitForRec*) ((*fTrackHits)[i1]);
1605 if (hit->GetHitNumber() < 0) { // combined cluster / track finder
1606 Int_t index = -hit->GetHitNumber() / 100000;
1607 Int_t iPos = -hit->GetHitNumber() - index * 100000;
1608 clus = (AliMUONRawCluster*) fgCombi->DetElem(index-1)->RawClusters()->UncheckedAt(iPos);
1609 } else {
1610 rawclusters = fgTrackReconstructor->GetMUONData()->RawClusters(hit->GetChamberNumber());
1611 clus = (AliMUONRawCluster*) rawclusters->UncheckedAt(hit->GetHitNumber());
1612 }
1613 if (clus->GetTrack(2) != -1) printf ("%5d", clus->GetTrack(2)%10000000);
1614 else printf ("%5s", " ");
1615 }
1616 }
1617 cout << endl;
1618 for (Int_t i1=0; i1<fNmbTrackHits; i1++) {
1619 //cout << ((AliMUONHitForRec*)((*fTrackHits)[i1]))->GetHitNumber() << " ";
1620 cout << ((AliMUONHitForRec*)((*fTrackHits)[i1]))->GetZ() << " ";
1621 //cout << fgHitForRec->IndexOf(((AliMUONHitForRec*)((*fTrackHits)[i1]))) << " ";
1622 }
1623 cout << endl;
1624 for (Int_t i1=0; i1<fNmbTrackHits; i1++) printf("%8.4f", (*fChi2Smooth)[i1]);
1625 cout << endl;
1626 cout << "---------------------------------------------------" << endl;
1627
1628 // Get covariance matrix
1629 /* Not needed - covariance matrix is not interesting to anybody
1630 *fCovariance = *fWeight;
1631 if (fCovariance->Determinant() != 0) {
1632 Int_t ifail;
1633 mnvertLocal(&((*fCovariance)(0,0)), fgkSize,fgkSize,fgkSize,ifail);
1634 } else {
1635 AliWarning(" Determinant fCovariance=0:" );
1636 }
1637 */
1638}
1639
1640 //__________________________________________________________________________
1641void AliMUONTrackK::MSLine(Double_t dZ, Double_t x0)
1642{
1643/// Adds multiple scattering in a thick layer for linear propagation
1644
1645 Double_t cosAlph = TMath::Cos((*fTrackPar)(2,0));
1646 Double_t tanAlph = TMath::Tan((*fTrackPar)(2,0));
1647 Double_t cosBeta = TMath::Cos((*fTrackPar)(3,0));
1648 Double_t sinBeta;
1649 sinBeta = TMath::Sin((*fTrackPar)(3,0));
1650 Double_t tanBeta = TMath::Tan((*fTrackPar)(3,0));
1651 Double_t momentum = 1/(*fTrackPar)(4,0);
1652 Double_t velo = 1; // relativistic velocity
1653 Double_t step = TMath::Abs(dZ/cosAlph/cosBeta); // step length
1654
1655 // Projected scattering angle
1656 Double_t theta0 = 0.0136/velo/momentum/TMath::Sqrt(x0)*(1+0.038*TMath::Log(step/x0));
1657 Double_t theta02 = theta0*theta0;
1658 Double_t dl2 = step*step/2*theta02;
1659 Double_t dl3 = dl2*step*2/3;
1660
1661 //Derivatives
1662 Double_t dYdT = 1/cosAlph;
1663 Double_t dYdB = 0; //(*fTrackPar)(2,0)*sinBeta/cosAlph;
1664 Double_t dXdT = tanAlph*tanBeta;
1665 //Double_t dXdB = (1+(*fTrackPar)(2,0)*tanAlph*sinBeta*sinBeta)/cosBeta;
1666 Double_t dXdB = 1/cosBeta;
1667 Double_t dAdT = 1/cosBeta;
1668 Double_t dAdB = 0; //(*fTrackPar)(2,0)*tanBeta;
1669
1670 // Get covariance matrix
1671 *fCovariance = *fWeight;
1672 if (fCovariance->Determinant() != 0) {
1673 // fCovariance->Invert();
1674 Int_t ifail;
1675 mnvertLocal(&((*fCovariance)(0,0)), fgkSize,fgkSize,fgkSize,ifail);
1676 } else {
1677 AliWarning(" Determinant fCovariance=0:" );
1678 }
1679
1680 (*fCovariance)(0,0) += dl3*(dYdT*dYdT+dYdB*dYdB); // <yy>
1681 (*fCovariance)(1,1) += dl3*(dXdT*dXdT+dXdB*dXdB); // <xx>
1682 (*fCovariance)(2,2) += theta02*step*(dAdT*dAdT+dAdB*dAdB); // <aa>
1683 (*fCovariance)(3,3) += theta02*step; // <bb>
1684
1685 (*fCovariance)(0,1) += dl3*(dYdT*dXdT+dYdB*dXdB); // <yx>
1686 (*fCovariance)(1,0) = (*fCovariance)(0,1);
1687
1688 (*fCovariance)(0,2) += dl2*(dYdT*dAdT+dYdB*dAdB); // <ya>
1689 (*fCovariance)(2,0) = (*fCovariance)(0,2);
1690
1691 (*fCovariance)(0,3) += dl2*dYdB; // <yb>
1692 (*fCovariance)(3,0) = (*fCovariance)(0,3);
1693
1694 (*fCovariance)(1,2) += dl2*(-dXdT*dAdT+dXdB*dAdB); // <xa>
1695 (*fCovariance)(2,1) = (*fCovariance)(1,2);
1696
1697 (*fCovariance)(1,3) += dl2*dXdB; // <xb>
1698 (*fCovariance)(3,1) = (*fCovariance)(1,3);
1699
1700 (*fCovariance)(2,3) += theta02*step*dAdB; // <ab>
1701 (*fCovariance)(3,2) = (*fCovariance)(2,3);
1702
1703 // Get weight matrix
1704 *fWeight = *fCovariance;
1705 if (fWeight->Determinant() != 0) {
1706 Int_t ifail;
1707 mnvertLocal(&((*fWeight)(0,0)), fgkSize,fgkSize,fgkSize,ifail);
1708 } else {
1709 AliWarning(" Determinant fWeight=0:");
1710 }
1711}
1712
1713 //__________________________________________________________________________
1714Bool_t AliMUONTrackK::Recover(void)
1715{
1716/// Adds new failed track(s) which can be tried to be recovered
1717 Int_t nRecTracks;
1718 TClonesArray *trackPtr;
1719 AliMUONTrackK *trackK;
1720
1721 if (fgDebug > 0) cout << " ******** Enter Recover " << endl;
1722 trackPtr = fgTrackReconstructor->GetRecTracksPtr();
1723
1724 // Remove hit with the highest chi2
1725 Double_t chi2 = 0;
1726 if (fgDebug > 0) {
1727 for (Int_t i=0; i<fNmbTrackHits; i++) {
1728 chi2 = fChi2Smooth ? (*fChi2Smooth)[i] : (*fChi2Array)[i];
1729 printf("%10.4f", chi2);
1730 }
1731 printf("\n");
1732 for (Int_t i=0; i<fNmbTrackHits; i++) {
1733 printf("%10d", ((AliMUONHitForRec*)fTrackHits->UncheckedAt(i))->GetChamberNumber());
1734 }
1735 printf("\n");
1736 }
1737 Double_t chi2max = 0;
1738 Int_t imax = 0;
1739 for (Int_t i=0; i<fNmbTrackHits; i++) {
1740 chi2 = fChi2Smooth ? (*fChi2Smooth)[i] : (*fChi2Array)[i];
1741 if (chi2 < chi2max) continue;
1742 chi2max = chi2;
1743 imax = i;
1744 }
1745 //if (chi2max < 10) return kFALSE; // !!!
1746 //if (chi2max < 25) imax = fNmbTrackHits - 1;
1747 if (chi2max < 15) imax = fNmbTrackHits - 1; // discard the last point
1748 // Check if the outlier is not from the seed segment
1749 AliMUONHitForRec *skipHit = (AliMUONHitForRec*) fTrackHits->UncheckedAt(imax);
1750 if (skipHit == (AliMUONHitForRec*) fStartSegment->First() || skipHit == (AliMUONHitForRec*) fStartSegment->Second()) {
1751 //DropBranches(fStartSegment); // drop all tracks with the same seed segment
1752 return kFALSE; // to be changed probably
1753 }
1754
1755 // Make a copy of track hit collection
1756 TObjArray *hits = new TObjArray(*fTrackHits);
1757 Int_t imax0;
1758 imax0 = imax;
1759
1760 // Hits after the found one will be removed
1761 if (GetStation0() == 3 && skipHit->GetChamberNumber() >= 7) {
1762 SortHits(1, fTrackHits); // unsort hits
1763 imax = fTrackHits->IndexOf(skipHit);
1764 }
1765 Int_t nTrackHits = fNmbTrackHits;
1766 for (Int_t i=imax+1; i<nTrackHits; i++) {
1767 AliMUONHitForRec *hit = (AliMUONHitForRec*) fTrackHits->UncheckedAt(i);
1768 fTrackHits->Remove(hit);
1769 hit->SetNTrackHits(hit->GetNTrackHits()-1); // unmark hit
1770 fNmbTrackHits--;
1771 }
1772
1773 // Check if the track candidate doesn't exist yet
1774 if (ExistDouble()) { delete hits; return kFALSE; }
1775
1776 //DropBranches(imax0, hits); // drop branches downstream the discarded hit
1777 delete hits;
1778
1779 nRecTracks = fgTrackReconstructor->GetNRecTracks();
1780 skipHit = (AliMUONHitForRec*) ((*fTrackHits)[fNmbTrackHits-1]);
1781 // Remove all saved steps and smoother matrices after the skipped hit
1782 RemoveMatrices(skipHit->GetZ());
1783
1784 //AZ(z->-z) if (skipHit->GetZ() > ((AliMUONHitForRec*) fStartSegment->Second())->GetZ() || !fNSteps) {
1785 if (TMath::Abs(skipHit->GetZ()) > TMath::Abs( ((AliMUONHitForRec*) fStartSegment->Second())->GetZ()) || !fNSteps) {
1786 // Propagation toward high Z or skipped hit next to segment -
1787 // start track from segment
1788 trackK = new ((*trackPtr)[nRecTracks]) AliMUONTrackK(fStartSegment);
1789 fgTrackReconstructor->SetNRecTracks(nRecTracks+1);
1790 trackK->fRecover = 1;
1791 trackK->fSkipHit = skipHit;
1792 trackK->fNmbTrackHits = fNmbTrackHits;
1793 delete trackK->fTrackHits; // not efficient ?
1794 trackK->fTrackHits = new TObjArray(*fTrackHits);
1795 if (fgDebug > 0) cout << nRecTracks << " " << trackK->fRecover << endl;
1796 return kTRUE;
1797 }
1798
1799 trackK = new ((*trackPtr)[nRecTracks]) AliMUONTrackK(NULL, NULL);
1800 *trackK = *this;
1801 fgTrackReconstructor->SetNRecTracks(nRecTracks+1);
1802 //AZ(z->-z) trackK->fTrackDir = -1;
1803 trackK->fTrackDir = 1;
1804 trackK->fRecover = 1;
1805 trackK->fSkipHit = skipHit;
1806 Int_t iD = trackK->fParFilter->Last()->GetUniqueID();
1807 if (iD > 1) {
1808 trackK->fParFilter->Last()->SetUniqueID(iD-1);
1809 CreateMatrix(trackK->fParFilter);
1810 }
1811 *((TMatrixD*)trackK->fParFilter->Last()) = *((TMatrixD*)fParExtrap->Last());
1812 trackK->fParFilter->Last()->SetUniqueID(1);
1813 *(trackK->fTrackPar) = *((TMatrixD*)trackK->fParFilter->Last());
1814 iD = trackK->fCovFilter->Last()->GetUniqueID();
1815 if (iD > 1) {
1816 trackK->fCovFilter->Last()->SetUniqueID(iD-1);
1817 CreateMatrix(trackK->fCovFilter);
1818 }
1819 *((TMatrixD*)trackK->fCovFilter->Last()) = *((TMatrixD*)fCovExtrap->Last());
1820 trackK->fCovFilter->Last()->SetUniqueID(1);
1821 *(trackK->fWeight) = *((TMatrixD*)trackK->fCovFilter->Last());
1822 if (trackK->fWeight->Determinant() != 0) {
1823 Int_t ifail;
1824 mnvertLocal(&((*(trackK->fWeight))(0,0)), fgkSize,fgkSize,fgkSize,ifail);
1825 } else {
1826 AliWarning(" Determinant fWeight=0:");
1827 }
1828 trackK->fPosition = trackK->fPositionNew = (*fSteps)[fNSteps-1];
1829 trackK->fChi2 = 0;
1830 for (Int_t i=0; i<fNmbTrackHits-1; i++) trackK->fChi2 += (*fChi2Array)[i];
1831 if (fgDebug > 0) cout << nRecTracks << " " << trackK->fRecover << endl;
1832 return kTRUE;
1833}
1834
1835 //__________________________________________________________________________
1836void AliMUONTrackK::AddMatrices(AliMUONTrackK *trackK, Double_t dChi2, AliMUONHitForRec *hitAdd)
1837{
1838/// Adds matrices for the smoother and keep Chi2 for the point
1839/// Track parameters
1840 //trackK->fParFilter->Last()->Print();
1841 Int_t iD = trackK->fParFilter->Last()->GetUniqueID();
1842 if (iD > 1) {
1843 trackK->fParFilter->Last()->SetUniqueID(iD-1);
1844 CreateMatrix(trackK->fParFilter);
1845 iD = 1;
1846 }
1847 *((TMatrixD*)(trackK->fParFilter->Last())) = *(trackK->fTrackPar);
1848 trackK->fParFilter->Last()->SetUniqueID(iD);
1849 if (fgDebug > 1) {
1850 cout << " Add matrices" << " " << fPosition << " " << dChi2 << " " << fgHitForRec->IndexOf(hitAdd) << endl;
1851 //trackK->fTrackPar->Print();
1852 //trackK->fTrackParNew->Print();
1853 trackK->fParFilter->Last()->Print();
1854 cout << " Add matrices" << endl;
1855 }
1856 // Covariance
1857 *(trackK->fCovariance) = *(trackK->fWeight);
1858 if (trackK->fCovariance->Determinant() != 0) {
1859 Int_t ifail;
1860 mnvertLocal(&((*(trackK->fCovariance))(0,0)), fgkSize,fgkSize,fgkSize,ifail);
1861 } else {
1862 AliWarning(" Determinant fCovariance=0:");
1863 }
1864 iD = trackK->fCovFilter->Last()->GetUniqueID();
1865 if (iD > 1) {
1866 trackK->fCovFilter->Last()->SetUniqueID(iD-1);
1867 CreateMatrix(trackK->fCovFilter);
1868 iD = 1;
1869 }
1870 *((TMatrixD*)(trackK->fCovFilter->Last())) = *(trackK->fCovariance);
1871 trackK->fCovFilter->Last()->SetUniqueID(iD);
1872
1873 // Save Chi2-increment for point
1874 Int_t indx = trackK->fTrackHits->IndexOf(hitAdd);
1875 if (indx < 0) indx = fNmbTrackHits;
1876 if (trackK->fChi2Array->fN <= indx) trackK->fChi2Array->Set(indx+10);
1877 trackK->fChi2Array->AddAt(dChi2,indx);
1878}
1879
1880 //__________________________________________________________________________
1881void AliMUONTrackK::CreateMatrix(TObjArray *objArray) const
1882{
1883/// Create new matrix and add it to TObjArray
1884
1885 TMatrixD *matrix = (TMatrixD*) objArray->First();
1886 TMatrixD *tmp = new TMatrixD(*matrix);
1887 objArray->AddAtAndExpand(tmp,objArray->GetLast());
1888 //cout << " !!! " << tmp << " " << tmp->GetUniqueID() << " " << tmp->GetNoElements() << endl;
1889}
1890
1891 //__________________________________________________________________________
1892void AliMUONTrackK::RemoveMatrices(Double_t zEnd)
1893{
1894/// Remove matrices (and saved steps) in the smoother part with abs(z) < abs(zEnd)
1895
1896 for (Int_t i=fNSteps-1; i>=0; i--) {
1897 if (fgDebug > 1) printf("%10.4f %10.4f \n", (*fSteps)[i], zEnd);
1898 if (TMath::Abs((*fSteps)[i]) > TMath::Abs(zEnd) - 0.01) break;
1899 RemoveMatrices(this);
1900 } // for (Int_t i=fNSteps-1;
1901}
1902
1903 //__________________________________________________________________________
1904void AliMUONTrackK::RemoveMatrices(AliMUONTrackK* trackK)
1905{
1906/// Remove last saved matrices and steps in the smoother part
1907
1908 trackK->fNSteps--;
1909 Int_t i = trackK->fNSteps;
1910
1911 Int_t id = 0;
1912 // Delete only matrices not shared by several tracks
1913 id = trackK->fParExtrap->Last()->GetUniqueID();
1914 if (id > 1) {
1915 trackK->fParExtrap->Last()->SetUniqueID(id-1);
1916 trackK->fParExtrap->RemoveAt(i);
1917 }
1918 else ((TMatrixD*)(trackK->fParExtrap->RemoveAt(i)))->Delete();
1919 id = fParFilter->Last()->GetUniqueID();
1920 if (id > 1) {
1921 trackK->fParFilter->Last()->SetUniqueID(id-1);
1922 trackK->fParFilter->RemoveAt(i);
1923 }
1924 else ((TMatrixD*)(trackK->fParFilter->RemoveAt(i)))->Delete();
1925 id = trackK->fCovExtrap->Last()->GetUniqueID();
1926 if (id > 1) {
1927 trackK->fCovExtrap->Last()->SetUniqueID(id-1);
1928 trackK->fCovExtrap->RemoveAt(i);
1929 }
1930 else ((TMatrixD*)(trackK->fCovExtrap->RemoveAt(i)))->Delete();
1931 id = trackK->fCovFilter->Last()->GetUniqueID();
1932 if (id > 1) {
1933 trackK->fCovFilter->Last()->SetUniqueID(id-1);
1934 trackK->fCovFilter->RemoveAt(i);
1935 }
1936 else ((TMatrixD*)(trackK->fCovFilter->RemoveAt(i)))->Delete();
1937 id = trackK->fJacob->Last()->GetUniqueID();
1938 if (id > 1) {
1939 trackK->fJacob->Last()->SetUniqueID(id-1);
1940 trackK->fJacob->RemoveAt(i);
1941 }
1942 else ((TMatrixD*)(trackK->fJacob->RemoveAt(i)))->Delete();
1943}
1944
1945 //__________________________________________________________________________
1946Bool_t AliMUONTrackK::Smooth(void)
1947{
1948/// Apply smoother
1949 Int_t ihit = fNmbTrackHits - 1;
1950 AliMUONHitForRec *hit = (AliMUONHitForRec*) (*fTrackHits)[ihit];
1951 fChi2Smooth = new TArrayD(fNmbTrackHits);
1952 fChi2Smooth->Reset(-1);
1953 fChi2 = 0;
1954 fParSmooth = new TObjArray(15);
1955 fParSmooth->Clear();
1956
1957 if (fgDebug > 0) {
1958 cout << " ******** Enter Smooth " << endl;
1959 cout << (*fSteps)[fNSteps-1] << " " << fPosition << " " << fPositionNew << endl;
1960 /*
1961 for (Int_t i=fNSteps-1; i>=0; i--) cout << (*fSteps)[i] << " ";
1962 cout << endl;
1963 for (Int_t i=fNSteps-1; i>=0; i--) {cout << i << " " << (*fSteps)[i]; ((TMatrixD*)fParFilter->UncheckedAt(i))->Print(); ((TMatrixD*)fParExtrap->UncheckedAt(i))->Print(); ((TMatrixD*)fJacob->UncheckedAt(i))->Print();}
1964 */
1965 for (Int_t i=ihit; i>=0; i--) cout << ((AliMUONHitForRec*)(*fTrackHits)[i])->GetZ() << " ";
1966 cout << endl;
1967 }
1968
1969 // Find last point corresponding to the last hit
1970 Int_t iLast = fNSteps - 1;
1971 for ( ; iLast>=0; iLast--) {
1972 //AZ(z->-z) if ((*fSteps)[iLast] + 0.01 > ((AliMUONHitForRec*)(*fTrackHits)[fNmbTrackHits-1])->GetZ()) break;
1973 if (TMath::Abs((*fSteps)[iLast]) + 0.01 > TMath::Abs(((AliMUONHitForRec*)(*fTrackHits)[fNmbTrackHits-1])->GetZ())) break;
1974 }
1975
1976 TMatrixD parSmooth = *((TMatrixD*)(fParFilter->UncheckedAt(iLast))); // last filtered = first smoothed
1977 //parSmooth.Dump();
1978 TMatrixD covSmooth = *((TMatrixD*)(fCovFilter->UncheckedAt(iLast))); // last filtered = first smoothed
1979 TMatrixD tmp=covSmooth, weight=covSmooth, gain = covSmooth;
1980 TMatrixD tmpPar = *fTrackPar;
1981 //parSmooth.Print(); ((TMatrixD*)(fParExtrap->Last()))->Print();
1982
1983 Bool_t found;
1984 Double_t chi2max = 0;
1985 for (Int_t i=iLast+1; i>0; i--) {
1986 if (i == iLast + 1) goto L33; // temporary fix
1987
1988 // Smoother gain matrix
1989 weight = *((TMatrixD*)(fCovExtrap->UncheckedAt(i)));
1990 if (weight.Determinant() != 0) {
1991 Int_t ifail;
1992 mnvertLocal(&(weight(0,0)), fgkSize,fgkSize,fgkSize,ifail);
1993 } else {
1994 AliWarning(" Determinant weight=0:");
1995 }
1996 // Fk'Wkk+1
1997 tmp = TMatrixD(*((TMatrixD*)(fJacob->UncheckedAt(i))),TMatrixD::kTransposeMult,weight);
1998 gain = TMatrixD(*((TMatrixD*)(fCovFilter->UncheckedAt(i-1))),TMatrixD::kMult,tmp);
1999 //cout << (*fSteps)[i-1] << " " << (*fSteps)[i] << endl; gain.Print();
2000
2001 // Smoothed parameter vector
2002 //((TMatrixD*)(fParExtrap->UncheckedAt(i)))->Dump();
2003 parSmooth -= *((TMatrixD*)(fParExtrap->UncheckedAt(i)));
2004 tmpPar = *((TMatrixD*)(fParFilter->UncheckedAt(i-1)));
2005 tmpPar += TMatrixD(gain,TMatrixD::kMult,parSmooth);
2006 parSmooth = tmpPar;
2007
2008 // Smoothed covariance
2009 covSmooth -= *((TMatrixD*)(fCovExtrap->UncheckedAt(i)));
2010 weight = TMatrixD(TMatrixD::kTransposed,gain);
2011 tmp = TMatrixD(covSmooth,TMatrixD::kMult,weight);
2012 covSmooth = *((TMatrixD*)(fCovFilter->UncheckedAt(i-1)));
2013 covSmooth += TMatrixD(gain,TMatrixD::kMult,tmp);
2014
2015 // Check if there was a measurement at given z
2016 found = kFALSE;
2017 for ( ; ihit>=0; ihit--) {
2018 hit = (AliMUONHitForRec*) (*fTrackHits)[ihit];
2019 if (TMath::Abs(hit->GetZ()-(*fSteps)[i-1]) < 0.1) { found = kTRUE; break; }
2020 //AZ(z->-z) else if (ihit < fNmbTrackHits-1 && hit->GetZ() > (*fSteps)[i-1]) { ihit++; break; }
2021 else if (ihit < fNmbTrackHits-1 && TMath::Abs(hit->GetZ()) > TMath::Abs((*fSteps)[i-1])) { ihit++; break; }
2022 }
2023 if (!found) continue; // no measurement - skip the rest
2024 else if (fgDebug > 1) cout << i << " " << ihit << " " << hit->GetZ() << endl;
2025 if (ihit == 0) continue; // the first hit - skip the rest
2026
2027L33:
2028 // Smoothed residuals
2029 tmpPar = 0;
2030 tmpPar(0,0) = hit->GetBendingCoor() - parSmooth(0,0); // bending coordinate
2031 tmpPar(1,0) = hit->GetNonBendingCoor() - parSmooth(1,0); // non-bending coordinate
2032 if (fgDebug > 1) {
2033 cout << hit->GetBendingCoor() << " " << parSmooth(0,0) << " " << tmpPar(0,0) << endl;
2034 cout << hit->GetNonBendingCoor() << " " << parSmooth(1,0) << " " << tmpPar(1,0) << endl;
2035 }
2036 // Cov. matrix of smoothed residuals
2037 tmp = 0;
2038 tmp(0,0) = hit->GetBendingReso2() - covSmooth(0,0);
2039 tmp(1,1) = hit->GetNonBendingReso2() - covSmooth(1,1);
2040 tmp(0,1) = tmp(1,0) = -covSmooth(0,1);
2041 tmp(2,2) = tmp(3,3) = tmp(4,4) = 1;
2042
2043 // Calculate Chi2 of smoothed residuals
2044 if (tmp.Determinant() != 0) {
2045 Int_t ifail;
2046 mnvertLocal(&(tmp(0,0)), fgkSize,fgkSize,fgkSize,ifail);
2047 } else {
2048 AliWarning(" Determinant tmp=0:");
2049 }
2050 TMatrixD vector(tmp,TMatrixD::kMult,tmpPar);
2051 TMatrixD chi2(tmpPar,TMatrixD::kTransposeMult,vector);
2052 if (fgDebug > 1) chi2.Print();
2053 (*fChi2Smooth)[ihit] = chi2(0,0);
2054 if (chi2(0,0) > chi2max) chi2max = chi2(0,0);
2055 fChi2 += chi2(0,0);
2056 if (chi2(0,0) < 0) {
2057 //chi2.Print();
2058 AliError(Form(" *** chi2 < 0: %d %d ", i, iLast));
2059 }
2060 // Save smoothed parameters
2061 TMatrixD *par = new TMatrixD(parSmooth);
2062 fParSmooth->AddAtAndExpand(par, ihit);
2063
2064 } // for (Int_t i=iLast+1;
2065
2066 //if (chi2max > 16) {
2067 //if (chi2max > 25) {
2068 //if (chi2max > 50) {
2069 //if (chi2max > 100) {
2070 if (chi2max > fgkChi2max) {
2071 //if (Recover()) DropBranches();
2072 //Recover();
2073 Outlier();
2074 return kFALSE;
2075 }
2076 return kTRUE;
2077}
2078
2079 //__________________________________________________________________________
2080void AliMUONTrackK::Outlier()
2081{
2082/// Adds new track with removed hit having the highest chi2
2083
2084 if (fgDebug > 0) {
2085 cout << " ******** Enter Outlier " << endl;
2086 for (Int_t i=0; i<fNmbTrackHits; i++) printf("%10.4f", (*fChi2Smooth)[i]);
2087 printf("\n");
2088 for (Int_t i=0; i<fNmbTrackHits; i++) {
2089 printf("%10d", ((AliMUONHitForRec*)fTrackHits->UncheckedAt(i))->GetChamberNumber());
2090 }
2091 printf("\n");
2092 }
2093
2094 Double_t chi2max = 0;
2095 Int_t imax = 0;
2096 for (Int_t i=0; i<fNmbTrackHits; i++) {
2097 if ((*fChi2Smooth)[i] < chi2max) continue;
2098 chi2max = (*fChi2Smooth)[i];
2099 imax = i;
2100 }
2101 // Check if the outlier is not from the seed segment
2102 AliMUONHitForRec *hit = (AliMUONHitForRec*) fTrackHits->UncheckedAt(imax);
2103 if (hit == (AliMUONHitForRec*) fStartSegment->First() || hit == (AliMUONHitForRec*) fStartSegment->Second()) return; // to be changed probably
2104
2105 // Check for missing station
2106 Int_t ok = 1;
2107 if (imax == 0) {
2108 if (((AliMUONHitForRec*)fTrackHits->UncheckedAt(1))->GetChamberNumber() < 8) ok--;
2109 } else if (imax == fNmbTrackHits-1) {
2110 if (((AliMUONHitForRec*)fTrackHits->UncheckedAt(fNmbTrackHits-2))->GetChamberNumber() > 1) ok--;
2111 }
2112 else if (((AliMUONHitForRec*)fTrackHits->UncheckedAt(imax-1))->GetChamberNumber()/2 - ((AliMUONHitForRec*)fTrackHits->UncheckedAt(imax+1))->GetChamberNumber()/2 > 1) ok--;
2113 if (!ok) { Recover(); return; } // try to recover track
2114 //AZ if (!ok) { if (fgDebug >= 0) cout << imax << endl; DropBranches(imax, 0); return; }
2115
2116 // Remove saved steps and smoother matrices after the outlier
2117 RemoveMatrices(hit->GetZ());
2118
2119 // Check for possible double track candidates
2120 //if (ExistDouble(hit)) return;
2121
2122 TClonesArray *trackPtr = fgTrackReconstructor->GetRecTracksPtr();
2123 Int_t nRecTracks = fgTrackReconstructor->GetNRecTracks();
2124
2125 AliMUONTrackK *trackK = 0;
2126 if (!fNSteps || GetStation0() == 3 && hit->GetChamberNumber() > 7) {
2127 // start track from segment
2128 trackK = new ((*trackPtr)[nRecTracks]) AliMUONTrackK(fStartSegment);
2129 fgTrackReconstructor->SetNRecTracks(nRecTracks+1);
2130 trackK->fRecover = 2;
2131 trackK->fSkipHit = hit;
2132 trackK->fNmbTrackHits = fNmbTrackHits;
2133
2134 hit = (AliMUONHitForRec*) trackK->fTrackHits->UncheckedAt(0);
2135 hit->SetNTrackHits(hit->GetNTrackHits()-1);
2136 hit = (AliMUONHitForRec*) trackK->fTrackHits->UncheckedAt(1);
2137 hit->SetNTrackHits(hit->GetNTrackHits()-1);
2138 delete trackK->fTrackHits; // not efficient ?
2139 trackK->fTrackHits = new TObjArray(*fTrackHits);
2140 for (Int_t i = 0; i < fNmbTrackHits; i++) {
2141 hit = (AliMUONHitForRec*) fTrackHits->UncheckedAt(i);
2142 hit->SetNTrackHits(hit->GetNTrackHits()+1);
2143 }
2144
2145 if (GetStation0() == 3) trackK->SortHits(1, trackK->fTrackHits);
2146 if (fgDebug > 0) cout << nRecTracks << " Outlier" << endl;
2147 return;
2148 }
2149 trackK = new ((*trackPtr)[nRecTracks]) AliMUONTrackK(NULL, NULL);
2150 *trackK = *this;
2151 fgTrackReconstructor->SetNRecTracks(nRecTracks+1);
2152 trackK->fTrackDir = 1;
2153 trackK->fRecover = 2;
2154 trackK->fSkipHit = hit;
2155 Int_t iD = trackK->fParFilter->Last()->GetUniqueID();
2156 if (iD > 1) {
2157 trackK->fParFilter->Last()->SetUniqueID(iD-1);
2158 CreateMatrix(trackK->fParFilter);
2159 }
2160 *((TMatrixD*)trackK->fParFilter->Last()) = *((TMatrixD*)fParExtrap->Last());
2161 trackK->fParFilter->Last()->SetUniqueID(1);
2162 *(trackK->fTrackPar) = *((TMatrixD*)trackK->fParFilter->Last());
2163 iD = trackK->fCovFilter->Last()->GetUniqueID();
2164 if (iD > 1) {
2165 trackK->fCovFilter->Last()->SetUniqueID(iD-1);
2166 CreateMatrix(trackK->fCovFilter);
2167 }
2168 *((TMatrixD*)trackK->fCovFilter->Last()) = *((TMatrixD*)fCovExtrap->Last());
2169 trackK->fCovFilter->Last()->SetUniqueID(1);
2170 *(trackK->fWeight) = *((TMatrixD*)trackK->fCovFilter->Last());
2171 if (trackK->fWeight->Determinant() != 0) {
2172 Int_t ifail;
2173 mnvertLocal(&((*(trackK->fWeight))(0,0)), fgkSize,fgkSize,fgkSize,ifail);
2174 } else {
2175 AliWarning(" Determinant fWeight=0:");
2176 }
2177 trackK->fPosition = trackK->fPositionNew = (*fSteps)[fNSteps-1];
2178 trackK->fChi2 = 0;
2179 for (Int_t i=0; i<fNmbTrackHits-1; i++) trackK->fChi2 += (*fChi2Array)[i];
2180 if (fgDebug > 0) cout << nRecTracks << " Outlier " << trackK->fChi2 << endl;
2181}
2182
2183 //__________________________________________________________________________
2184Double_t AliMUONTrackK::GetChi2PerPoint(Int_t iPoint) const
2185{
2186/// Return Chi2 at point
2187 return fChi2Smooth ? (*fChi2Smooth)[iPoint] : (*fChi2Array)[iPoint];
2188 //return 0.;
2189}
2190
2191 //__________________________________________________________________________
2192void AliMUONTrackK::Print(FILE *lun) const
2193{
2194/// Print out track information
2195
2196 Int_t flag = 1;
2197 AliMUONHitForRec *hit = 0;
2198 // from raw clusters
2199 AliMUONRawCluster *clus = 0;
2200 TClonesArray *rawclusters = 0;
2201 for (Int_t i1=0; i1<fNmbTrackHits; i1++) {
2202 hit = (AliMUONHitForRec*) ((*fTrackHits)[i1]);
2203 rawclusters = fgTrackReconstructor->GetMUONData()->RawClusters(hit->GetChamberNumber());
2204 clus = (AliMUONRawCluster*) rawclusters->UncheckedAt(hit->GetHitNumber());
2205 if (TMath::Abs(clus->GetTrack(1)-1) < 2) {
2206 if (clus->GetTrack(2)) flag = 2;
2207 continue;
2208 }
2209 if (clus->GetTrack(2) && TMath::Abs(clus->GetTrack(2)-1) < 2) {
2210 flag = 3;
2211 continue;
2212 }
2213 flag = 0;
2214 break;
2215
2216 Int_t sig[2]={1,1}, tid[2]={0};
2217 for (Int_t i1=0; i1<fNmbTrackHits; i1++) {
2218 if (GetChi2PerPoint(i1) < -0.1) continue;
2219 hit = (AliMUONHitForRec*) ((*fTrackHits)[i1]);
2220 rawclusters = fgTrackReconstructor->GetMUONData()->RawClusters(hit->GetChamberNumber());
2221 clus = (AliMUONRawCluster*) rawclusters->UncheckedAt(hit->GetHitNumber());
2222 for (Int_t j=0; j<2; j++) {
2223 tid[j] = clus->GetTrack(j+1) - 1;
2224 if (clus->GetTrack(j+1) < 0) { sig[j] = 0; tid[j] = 999; }
2225 }
2226 fprintf(lun,"%3d %3d %10.4f", gAlice->GetEvNumber(), hit->GetChamberNumber(), GetChi2PerPoint(i1));
2227 if (!(clus->GetTrack(2))) fprintf(lun, "%3d %3d", sig[0], tid[0]); // simple cluster
2228 else { // track overlap
2229 fprintf(lun, "%3d %3d", TMath::Max(sig[0],sig[1]), TMath::Min(tid[0],tid[1]));
2230 //if (tid[0] < 2) flag *= 2;
2231 //else if (tid[1] < 2) flag *= 3;
2232 }
2233 fprintf (lun, "%3d \n", flag);
2234 }
2235 }
2236}
2237
2238 //__________________________________________________________________________
2239void AliMUONTrackK::DropBranches(Int_t imax, TObjArray *hits)
2240{
2241/// Drop branches downstream of the skipped hit
2242 Int_t nRecTracks;
2243 TClonesArray *trackPtr;
2244 AliMUONTrackK *trackK;
2245
2246 trackPtr = fgTrackReconstructor->GetRecTracksPtr();
2247 nRecTracks = fgTrackReconstructor->GetNRecTracks();
2248 Int_t icand = trackPtr->IndexOf(this);
2249 if (!hits) hits = fTrackHits;
2250
2251 // Check if the track candidate doesn't exist yet
2252 for (Int_t i=icand+1; i<nRecTracks; i++) {
2253 trackK = (AliMUONTrackK*) ((*trackPtr)[i]);
2254 if (trackK->fNmbTrackHits == 2 && trackK->GetRecover() == 0) continue;
2255 if (trackK->GetRecover() < 0) continue;
2256
2257 if (trackK->fNmbTrackHits >= imax + 1) {
2258 for (Int_t j=0; j<=imax; j++) {
2259 //if (j != fNmbTrackHits-1 && (*trackK->fTrackHits)[j] != (*fTrackHits)[j]) break;
2260 if ((*trackK->fTrackHits)[j] != (*hits)[j]) break;
2261 if (j == imax) {
2262 if (hits != fTrackHits) {
2263 // Drop all branches downstream the hit (for Recover)
2264 trackK->SetRecover(-1);
2265 if (fgDebug >= 0 )cout << " Recover " << i << endl;
2266 continue;
2267 }
2268 // Check if the removal of the hit breaks the track
2269 Int_t ok = 1;
2270 if (imax == 0) {
2271 if (((AliMUONHitForRec*)trackK->fTrackHits->UncheckedAt(1))->GetChamberNumber() < 8) ok--; }
2272 else if (imax == trackK->fNmbTrackHits-1) continue;
2273 // else if (imax == trackK->fNmbTrackHits-1) {
2274 //if (((AliMUONHitForRec*)trackK->fTrackHits->UncheckedAt(trackK->fNmbTrackHits-2))->GetChamberNumber() > 1) ok--;
2275 //}
2276 else if (((AliMUONHitForRec*)trackK->fTrackHits->UncheckedAt(imax-1))->GetChamberNumber()/2 - ((AliMUONHitForRec*)trackK->fTrackHits->UncheckedAt(imax+1))->GetChamberNumber()/2 > 1) ok--;
2277 if (!ok) trackK->SetRecover(-1);
2278 }
2279 } // for (Int_t j=0;
2280 }
2281 } // for (Int_t i=0;
2282}
2283
2284 //__________________________________________________________________________
2285void AliMUONTrackK::DropBranches(AliMUONObjectPair *segment)
2286{
2287/// Drop all candidates with the same seed segment
2288 Int_t nRecTracks;
2289 TClonesArray *trackPtr;
2290 AliMUONTrackK *trackK;
2291 AliMUONObjectPair *trackKSegment;
2292
2293 trackPtr = fgTrackReconstructor->GetRecTracksPtr();
2294 nRecTracks = fgTrackReconstructor->GetNRecTracks();
2295 Int_t icand = trackPtr->IndexOf(this);
2296
2297 for (Int_t i=icand+1; i<nRecTracks; i++) {
2298 trackK = (AliMUONTrackK*) ((*trackPtr)[i]);
2299 trackKSegment = trackK->fStartSegment;
2300 if (trackK->fNmbTrackHits == 2 && trackK->GetRecover() == 0) continue;
2301 if (trackK->GetRecover() < 0) continue;
2302 if (trackKSegment->First() == segment->First() &&
2303 trackKSegment->Second() == segment->Second()) trackK->SetRecover(-1);
2304 }
2305 if (fgDebug >= 0) cout << " Drop segment " << endl;
2306}
2307
2308 //__________________________________________________________________________
2309AliMUONHitForRec* AliMUONTrackK::GetHitLastOk(void)
2310{
2311/// Return the hit where track stopped
2312
2313 if (!fNSteps) return (AliMUONHitForRec*)((*fTrackHits)[1]);
2314 return fSkipHit;
2315}
2316
2317 //__________________________________________________________________________
2318Int_t AliMUONTrackK::GetStation0(void)
2319{
2320/// Return seed station number
2321 return ((AliMUONHitForRec*) fStartSegment->First())->GetChamberNumber() / 2;
2322}
2323
2324 //__________________________________________________________________________
2325Bool_t AliMUONTrackK::ExistDouble(AliMUONHitForRec *hit)
2326{
2327/// Check if the track will make a double after outlier removal
2328
2329 TClonesArray *trackPtr = fgTrackReconstructor->GetRecTracksPtr();
2330 Int_t nRecTracks = fgTrackReconstructor->GetNRecTracks();
2331 TObjArray *hitArray = new TObjArray(*fTrackHits);
2332 TObjArray *hitArray1 = new TObjArray(*hitArray);
2333 hitArray1->Remove(hit);
2334 hitArray1->Compress();
2335
2336 Bool_t same = kFALSE;
2337 for (Int_t i=0; i<nRecTracks; i++) {
2338 AliMUONTrackK *trackK = (AliMUONTrackK*) ((*trackPtr)[i]);
2339 if (trackK->fNmbTrackHits == 2 && trackK->GetRecover() == 0) continue;
2340 if (trackK == this) continue;
2341 if (trackK->fNmbTrackHits == fNmbTrackHits || trackK->fNmbTrackHits == fNmbTrackHits-1) {
2342 TObjArray *hits = new TObjArray(*trackK->fTrackHits);
2343 same = kTRUE;
2344 if (trackK->fNmbTrackHits == fNmbTrackHits) {
2345 for (Int_t j=0; j<fNmbTrackHits; j++) {
2346 if (hits->UncheckedAt(j) != hitArray->UncheckedAt(j)) { same = kFALSE; break; }
2347 }
2348 if (same) { delete hits; break; }
2349 if (trackK->fSkipHit) {
2350 TObjArray *hits1 = new TObjArray(*hits);
2351 if (hits1->Remove(trackK->fSkipHit) > 0) {
2352 hits1->Compress();
2353 same = kTRUE;
2354 for (Int_t j=0; j<fNmbTrackHits-1; j++) {
2355 if (hits1->UncheckedAt(j) != hitArray1->UncheckedAt(j)) { same = kFALSE; break; }
2356 }
2357 if (same) { delete hits1; break; }
2358 }
2359 delete hits1;
2360 }
2361 } else {
2362 // Check with removed outlier
2363 same = kTRUE;
2364 for (Int_t j=0; j<fNmbTrackHits-1; j++) {
2365 if (hits->UncheckedAt(j) != hitArray1->UncheckedAt(j)) { same = kFALSE; break; }
2366 }
2367 if (same) { delete hits; break; }
2368 }
2369 delete hits;
2370 }
2371 } // for (Int_t i=0; i<nRecTracks;
2372 delete hitArray; delete hitArray1;
2373 if (same && fgDebug >= 0) cout << " Same" << endl;
2374 return same;
2375}
2376
2377 //__________________________________________________________________________
2378Bool_t AliMUONTrackK::ExistDouble(void)
2379{
2380/// Check if the track will make a double after recovery
2381
2382 TClonesArray *trackPtr = fgTrackReconstructor->GetRecTracksPtr();
2383 Int_t nRecTracks = fgTrackReconstructor->GetNRecTracks();
2384
2385 TObjArray *hitArray = new TObjArray(*fTrackHits);
2386 if (GetStation0() == 3) SortHits(0, hitArray); // sort
2387 //if (GetStation0() == 3) SortHits(1, hitArray); // unsort
2388
2389 Bool_t same = kFALSE;
2390 for (Int_t i=0; i<nRecTracks; i++) {
2391 AliMUONTrackK *trackK = (AliMUONTrackK*) ((*trackPtr)[i]);
2392 if (trackK->fNmbTrackHits == 2 && trackK->GetRecover() == 0) continue;
2393 if (trackK == this) continue;
2394 //AZ if (trackK->GetRecover() < 0) continue; //
2395 if (trackK->fNmbTrackHits >= fNmbTrackHits) {
2396 TObjArray *hits = new TObjArray(*trackK->fTrackHits);
2397 if (trackK->GetStation0() == 3) SortHits(0, hits); // sort
2398 //if (trackK->GetStation0() == 3) SortHits(1, hits); // unsort
2399 for (Int_t j=0; j<fNmbTrackHits; j++) {
2400 //cout << fNmbTrackHits << " " << i << " " << j << " " << (*hitArray)[j] << " " << (*hits)[j] << " " << trackK->fSkipHit << endl;
2401 if (j != fNmbTrackHits-1 && (*hitArray)[j] != (*hits)[j]) break;
2402 if (j == fNmbTrackHits-1) {
2403 if (trackK->fSkipHit && TMath::Abs(((AliMUONHitForRec*)((*hitArray)[j]))->GetZ()-trackK->fSkipHit->GetZ()) < 0.5) same = kTRUE;
2404 //if (trackK->fNmbTrackHits > fNmbTrackHits &&
2405 //if (trackK->fSkipHit) cout << ((AliMUONHitForRec*)((*hitArray)[j]))->GetZ() << " " << trackK->fSkipHit->GetZ() << endl;
2406 }
2407 } // for (Int_t j=0;
2408 delete hits;
2409 if (same) break;
2410 }
2411 } // for (Int_t i=0;
2412 delete hitArray;
2413 return same;
2414}
2415
2416//______________________________________________________________________________
2417 void mnvertLocal(Double_t *a, Int_t l, Int_t, Int_t n, Int_t &ifail)
2418{
2419///*-*-*-*-*-*-*-*-*-*-*-*Inverts a symmetric matrix*-*-*-*-*-*-*-*-*-*-*-*-*
2420///*-* ==========================
2421///*-* inverts a symmetric matrix. matrix is first scaled to
2422///*-* have all ones on the diagonal (equivalent to change of units)
2423///*-* but no pivoting is done since matrix is positive-definite.
2424///*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
2425
2426 // taken from TMinuit package of Root (l>=n)
2427 // fVERTs, fVERTq and fVERTpp changed to localVERTs, localVERTq and localVERTpp
2428 // Double_t localVERTs[n], localVERTq[n], localVERTpp[n];
2429 Double_t * localVERTs = new Double_t[n];
2430 Double_t * localVERTq = new Double_t[n];
2431 Double_t * localVERTpp = new Double_t[n];
2432 // fMaxint changed to localMaxint
2433 Int_t localMaxint = n;
2434
2435 /* System generated locals */
2436 Int_t aOffset;
2437
2438 /* Local variables */
2439 Double_t si;
2440 Int_t i, j, k, kp1, km1;
2441
2442 /* Parameter adjustments */
2443 aOffset = l + 1;
2444 a -= aOffset;
2445
2446 /* Function Body */
2447 ifail = 0;
2448 if (n < 1) goto L100;
2449 if (n > localMaxint) goto L100;
2450//*-*- scale matrix by sqrt of diag elements
2451 for (i = 1; i <= n; ++i) {
2452 si = a[i + i*l];
2453 if (si <= 0) goto L100;
2454 localVERTs[i-1] = 1 / TMath::Sqrt(si);
2455 }
2456 for (i = 1; i <= n; ++i) {
2457 for (j = 1; j <= n; ++j) {
2458 a[i + j*l] = a[i + j*l]*localVERTs[i-1]*localVERTs[j-1];
2459 }
2460 }
2461//*-*- . . . start main loop . . . .
2462 for (i = 1; i <= n; ++i) {
2463 k = i;
2464//*-*- preparation for elimination step1
2465 if (a[k + k*l] != 0) localVERTq[k-1] = 1 / a[k + k*l];
2466 else goto L100;
2467 localVERTpp[k-1] = 1;
2468 a[k + k*l] = 0;
2469 kp1 = k + 1;
2470 km1 = k - 1;
2471 if (km1 < 0) goto L100;
2472 else if (km1 == 0) goto L50;
2473 else goto L40;
2474L40:
2475 for (j = 1; j <= km1; ++j) {
2476 localVERTpp[j-1] = a[j + k*l];
2477 localVERTq[j-1] = a[j + k*l]*localVERTq[k-1];
2478 a[j + k*l] = 0;
2479 }
2480L50:
2481 if (k - n < 0) goto L51;
2482 else if (k - n == 0) goto L60;
2483 else goto L100;
2484L51:
2485 for (j = kp1; j <= n; ++j) {
2486 localVERTpp[j-1] = a[k + j*l];
2487 localVERTq[j-1] = -a[k + j*l]*localVERTq[k-1];
2488 a[k + j*l] = 0;
2489 }
2490//*-*- elimination proper
2491L60:
2492 for (j = 1; j <= n; ++j) {
2493 for (k = j; k <= n; ++k) { a[j + k*l] += localVERTpp[j-1]*localVERTq[k-1]; }
2494 }
2495 }
2496//*-*- elements of left diagonal and unscaling
2497 for (j = 1; j <= n; ++j) {
2498 for (k = 1; k <= j; ++k) {
2499 a[k + j*l] = a[k + j*l]*localVERTs[k-1]*localVERTs[j-1];
2500 a[j + k*l] = a[k + j*l];
2501 }
2502 }
2503 delete [] localVERTs;
2504 delete [] localVERTq;
2505 delete [] localVERTpp;
2506 return;
2507//*-*- failure return
2508L100:
2509 delete [] localVERTs;
2510 delete [] localVERTq;
2511 delete [] localVERTpp;
2512 ifail = 1;
2513} /* mnvertLocal */
2514