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