]> git.uio.no Git - u/mrichter/AliRoot.git/blame - MUON/AliMUONTrackK.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Updated list of MUON libraries
[u/mrichter/AliRoot.git] / MUON / AliMUONTrackK.cxx
CommitLineData
83dbc640 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
b035a148 16/* $Id$ */
83dbc640 17
d19b6003 18// ---------------------
1ffe5558 19// Class AliMUONTrackK
d19b6003 20// ---------------------
4c08623d 21// Reconstructed track in the muon system based on the extended
1ffe5558 22// Kalman filter approach
1ffe5558 23// Author: Alexander Zinchenko, JINR Dubna
24
30178c30 25#include "AliMUONTrackK.h"
6b092dfc 26#include "AliMUONRecData.h"
2457f726 27#include "AliMUONConstants.h"
cc87ebcd 28
2457f726 29#include "AliMUONTrackReconstructorK.h"
83dbc640 30#include "AliMUONHitForRec.h"
208f139e 31#include "AliMUONObjectPair.h"
83dbc640 32#include "AliMUONRawCluster.h"
33#include "AliMUONTrackParam.h"
37827b29 34#include "AliMUONTrackExtrap.h"
cc87ebcd 35#include "AliMUONEventRecoCombi.h"
36#include "AliMUONDetElement.h"
c0e54947 37
8c343c7c 38#include "AliLog.h"
c0e54947 39
40#include <Riostream.h>
41#include <TClonesArray.h>
42#include <TArrayD.h>
83dbc640 43
7945aae7 44/// \cond CLASSIMP
45ClassImp(AliMUONTrackK) // Class implementation in ROOT context
46/// \endcond
47
343146bf 48const Int_t AliMUONTrackK::fgkSize = 5;
b21aa6eb 49const Int_t AliMUONTrackK::fgkNSigma = 12;
50const Double_t AliMUONTrackK::fgkChi2max = 100;
343146bf 51const Int_t AliMUONTrackK::fgkTriesMax = 10000;
52const Double_t AliMUONTrackK::fgkEpsilon = 0.002;
53
208f139e 54void mnvertLocal(Double_t* a, Int_t l, Int_t m, Int_t n, Int_t& ifail);
83dbc640 55
7945aae7 56Int_t AliMUONTrackK::fgDebug = -1; //-1;
83dbc640 57Int_t AliMUONTrackK::fgNOfPoints = 0;
2457f726 58AliMUONTrackReconstructorK* AliMUONTrackK::fgTrackReconstructor = NULL;
83dbc640 59TClonesArray* AliMUONTrackK::fgHitForRec = NULL;
cc87ebcd 60AliMUONEventRecoCombi *AliMUONTrackK::fgCombi = NULL;
96a96c2b 61
83dbc640 62 //__________________________________________________________________________
63AliMUONTrackK::AliMUONTrackK()
54d7ba50 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)
83dbc640 89{
d19b6003 90/// Default constructor
83dbc640 91
b21aa6eb 92 fgTrackReconstructor = NULL; // pointer to event reconstructor
83dbc640 93 fgHitForRec = NULL; // pointer to points
94 fgNOfPoints = 0; // number of points
95
83dbc640 96 return;
97}
98
99 //__________________________________________________________________________
2457f726 100AliMUONTrackK::AliMUONTrackK(AliMUONTrackReconstructorK *TrackReconstructor, TClonesArray *hitForRec)
54d7ba50 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)
83dbc640 126{
d19b6003 127/// Constructor
83dbc640 128
29f1b13a 129 if (!TrackReconstructor) return;
b21aa6eb 130 fgTrackReconstructor = TrackReconstructor; // pointer to event reconstructor
83dbc640 131 fgHitForRec = hitForRec; // pointer to points
132 fgNOfPoints = fgHitForRec->GetEntriesFast(); // number of points
cc87ebcd 133 fgCombi = NULL;
134 if (fgTrackReconstructor->GetTrackMethod() == 3) fgCombi = AliMUONEventRecoCombi::Instance();
83dbc640 135}
136
137 //__________________________________________________________________________
208f139e 138AliMUONTrackK::AliMUONTrackK(AliMUONObjectPair *segment)
139 : AliMUONTrack(NULL,NULL),
140 fStartSegment(new AliMUONObjectPair(*segment)),
54d7ba50 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)
83dbc640 164{
d19b6003 165/// Constructor from a segment
208f139e 166 Double_t dX, dY, dZ, bendingSlope, bendingImpact;
83dbc640 167 AliMUONHitForRec *hit1, *hit2;
168 AliMUONRawCluster *clus;
169 TClonesArray *rawclusters;
170
83dbc640 171 // Pointers to hits from the segment
208f139e 172 hit1 = (AliMUONHitForRec*) segment->First();
173 hit2 = (AliMUONHitForRec*) segment->Second();
83dbc640 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
b035a148 177 if (TMath::Abs(hit1->GetZ()) > TMath::Abs(hit2->GetZ())) {
83dbc640 178 hit1 = hit2;
208f139e 179 hit2 = (AliMUONHitForRec*) segment->First();
83dbc640 180 }
83dbc640 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
1ffe5558 188 fTrackHits->Add(hit2); // add hit 2
189 fTrackHits->Add(hit1); // add hit 1
b035a148 190 //AZ(Z->-Z) fTrackDir = -1;
191 fTrackDir = 1;
83dbc640 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
1ffe5558 197 fTrackHits->Add(hit1); // add hit 1
198 fTrackHits->Add(hit2); // add hit 2
b035a148 199 //AZ(Z->-Z) fTrackDir = 1;
200 fTrackDir = -1;
83dbc640 201 }
202 dZ = hit2->GetZ() - hit1->GetZ();
203 dY = hit2->GetBendingCoor() - hit1->GetBendingCoor();
204 dX = hit2->GetNonBendingCoor() - hit1->GetNonBendingCoor();
208f139e 205 bendingSlope = (hit2->GetBendingCoor() - hit1->GetBendingCoor()) / dZ;
206 bendingImpact = hit1->GetBendingCoor() - hit1->GetZ() * bendingSlope;
83dbc640 207 (*fTrackPar)(2,0) = TMath::ATan2(dY,dZ); // alpha
b21aa6eb 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();
208f139e 211 (*fTrackPar)(4,0) = 1./AliMUONTrackExtrap::GetBendingMomentumFromImpactParam(bendingImpact); // 1/Pt
83dbc640 212 (*fTrackPar)(4,0) *= TMath::Cos((*fTrackPar)(3,0)); // 1/p
b035a148 213 // Evaluate covariance (and weight) matrix
214 EvalCovariance(dZ);
215
b035a148 216 if (fgDebug < 0 ) return;
208f139e 217 cout << fgTrackReconstructor->GetNRecTracks()-1 << " "
218 << AliMUONTrackExtrap::GetBendingMomentumFromImpactParam(bendingImpact)
219 << " " << 1/(*fTrackPar)(4,0) << " ";
83dbc640 220 // from raw clusters
8f373020 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 << " <--> ";
83dbc640 228 }
208f139e 229 cout << " @ " << hit1->GetChamberNumber() << endl;
8f373020 230
83dbc640 231}
232
233 //__________________________________________________________________________
234AliMUONTrackK::~AliMUONTrackK()
235{
d19b6003 236/// Destructor
83dbc640 237
208f139e 238 if (fStartSegment) {
239 delete fStartSegment;
240 fStartSegment = 0x0;
241 }
242
1ffe5558 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;
83dbc640 251 }
b035a148 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;
83dbc640 319}
320
83dbc640 321 //__________________________________________________________________________
322AliMUONTrackK & AliMUONTrackK::operator=(const AliMUONTrackK& source)
323{
d19b6003 324/// Assignment operator
325
83dbc640 326 // Members
327 if(&source == this) return *this;
30178c30 328
329 // base class assignement
b035a148 330 //AZ TObject::operator=(source);
331 AliMUONTrack::operator=(source);
208f139e 332
333 if (fStartSegment) delete fStartSegment;
334 if (source.fStartSegment) fStartSegment = new AliMUONObjectPair(*(source.fStartSegment));
335 else fStartSegment = 0x0;
336
1ffe5558 337 fNmbTrackHits = source.fNmbTrackHits;
83dbc640 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
1ffe5558 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 }
83dbc640 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
b035a148 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 }
83dbc640 380 return *this;
381}
382
383 //__________________________________________________________________________
384void AliMUONTrackK::EvalCovariance(Double_t dZ)
385{
d19b6003 386/// Evaluate covariance (and weight) matrix for track candidate
83dbc640 387 Double_t sigmaB, sigmaNonB, tanA, tanB, dAdY, rad, dBdX, dBdY;
388
b21aa6eb 389 dZ = -dZ;
29f1b13a 390 sigmaB = fgTrackReconstructor->GetBendingResolution(); // bending resolution
391 sigmaNonB = fgTrackReconstructor->GetNonBendingResolution(); // non-bending resolution
83dbc640 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
83dbc640 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) {
83dbc640 416 // fWeight->Invert();
b21aa6eb 417 Int_t ifail;
b21aa6eb 418 mnvertLocal(&((*fWeight)(0,0)), fgkSize,fgkSize,fgkSize,ifail);
83dbc640 419 } else {
8c343c7c 420 AliWarning(" Determinant fWeight=0:");
83dbc640 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{
d19b6003 428/// Follows track through detector stations
83dbc640 429 Bool_t miss, success;
b035a148 430 Int_t ichamb, iFB, iMin, iMax, dChamb, ichambOK;
cc87ebcd 431 Int_t ihit, firstIndx = -1, lastIndx = -1, currIndx = -1, iz0 = -1, iz = -1;
83dbc640 432 Double_t zEnd, dChi2;
cc87ebcd 433 AliMUONHitForRec *hitAdd, *firstHit = NULL, *lastHit = NULL, *hit;
83dbc640 434 AliMUONRawCluster *clus;
435 TClonesArray *rawclusters;
436 hit = 0; clus = 0; rawclusters = 0;
437
b035a148 438 miss = success = kTRUE;
83dbc640 439 Int_t endOfProp = 0;
ef4d63ce 440 //iFB = (ichamEnd == ichamBeg) ? fTrackDir : TMath::Sign(1,ichamEnd-ichamBeg);
441 iFB = (ichamEnd == ichamBeg) ? -fTrackDir : TMath::Sign(1,ichamEnd-ichamBeg);
83dbc640 442 iMin = TMath::Min(ichamEnd,ichamBeg);
443 iMax = TMath::Max(ichamEnd,ichamBeg);
444 ichamb = ichamBeg;
445 ichambOK = ichamb;
446
83dbc640 447 if (Back) {
448 // backpropagation
b035a148 449 currIndx = 1;
1ffe5558 450 if (((AliMUONHitForRec*)fTrackHits->First())->GetChamberNumber() != ichamb) currIndx = 0;
b035a148 451 } else if (fRecover) {
452 hit = GetHitLastOk();
1ffe5558 453 currIndx = fTrackHits->IndexOf(hit);
208f139e 454 if (currIndx < 0) hit = (AliMUONHitForRec*) fStartSegment->First(); // for station 3
b035a148 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
1ffe5558 460 fTrackHits->Remove(fSkipHit); // remove hit
461 fNmbTrackHits --;
b035a148 462 fSkipHit->SetNTrackHits(fSkipHit->GetNTrackHits()-1); // unmark hit
463 if (fRecover == 1) {
464 // recovery
465 Back = kFALSE;
466 fRecover = 0;
1ffe5558 467 ichambOK = ((AliMUONHitForRec*)((*fTrackHits)[fNmbTrackHits-1]))->GetChamberNumber();
468 //if (ichambOK >= 8 && ((AliMUONHitForRec*)((*fTrackHits)[0]))->GetChamberNumber() == 6) ichambOK = 6;
cc87ebcd 469 if (fgTrackReconstructor->GetTrackMethod() == 3 &&
470 fSkipHit->GetHitNumber() < 0) {
471 iz0 = fgCombi->IZfromHit(fSkipHit);
472 currIndx = -1;
473 }
474 else currIndx = fgHitForRec->IndexOf(fSkipHit);
b035a148 475 } else {
476 // outlier
1ffe5558 477 fTrackHits->Compress();
b035a148 478 }
479 } // if (hit == fSkipHit)
1ffe5558 480 else if (currIndx < 0) currIndx = fTrackHits->IndexOf(hit);
b035a148 481 } // else if (fRecover)
482 else {
483 // Get indices of the 1'st and last hits on the track candidate
1ffe5558 484 firstHit = (AliMUONHitForRec*) fTrackHits->First();
485 lastHit = (AliMUONHitForRec*) fTrackHits->Last();
cc87ebcd 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 }
83dbc640 493 }
494
cc87ebcd 495 if (iz0 < 0) iz0 = iFB;
496 while (ichamb >= iMin && ichamb <= iMax) {
83dbc640 497 // Find the closest hit in Z, not belonging to the current plane
498 if (Back) {
499 // backpropagation
1ffe5558 500 if (currIndx < fNmbTrackHits) {
501 hitAdd = (AliMUONHitForRec*) fTrackHits->UncheckedAt(currIndx);
b035a148 502 zEnd = hitAdd->GetZ();
503 //AZ(z->-z) } else zEnd = -9999;
504 } else zEnd = 9999;
83dbc640 505 } else {
b035a148 506 //AZ(Z->-Z) zEnd = -9999;
507 zEnd = 9999;
83dbc640 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 }
cc87ebcd 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 }
83dbc640 535 }
b035a148 536 if (zEnd>999 && ichamb==ichamEnd) endOfProp = 1; // end-of-propagation
83dbc640 537 else {
b035a148 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;
83dbc640 541 ichamb += iFB;
f93b9bd8 542 zEnd = AliMUONConstants::DefaultChamberZ(ichamb);
83dbc640 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);
b035a148 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) {
1ffe5558 557 //if (dStatMiss == 2 && ichambOK/2 != 3 || dStatMiss > 2) { // AZ - missing st. 3
b035a148 558 // missing station - abandon track
559 //cout << dChamb << " " << ichambOK << " " << fgNOfPoints << " " << 1/(*fTrackPar)(4,0) << endl;
560 if (fgDebug >= 10) {
83dbc640 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() << " " << " ";
29fc2c86 566 cout << ((AliMUONHitForRec*)((*fgHitForRec)[i1]))->GetTTRTrack() << endl;
83dbc640 567 }
b035a148 568 cout << endl;
1ffe5558 569 cout << fNmbTrackHits << endl;
570 for (Int_t i1=0; i1<fNmbTrackHits; i1++) {
571 hit = (AliMUONHitForRec*) ((*fTrackHits)[i1]);
83dbc640 572 printf(" * %d %10.4f %10.4f %10.4f",
b035a148 573 hit->GetChamberNumber(), hit->GetBendingCoor(),
83dbc640 574 hit->GetNonBendingCoor(), hit->GetZ());
8f373020 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
83dbc640 584 }
b035a148 585 } // if (fgDebug >= 10)
1ffe5558 586 if (fNmbTrackHits>2 && fRecover==0) Recover(); // try to recover track later
b035a148 587 return kFALSE;
588 } // if (dStatMiss > 1)
589 } // if (!Back)
83dbc640 590 if (endOfProp != 0) break;
591
592 // propagate to the found Z
593
594 // Check if track steps into dipole
b035a148 595 //AZ(z->-z) if (fPosition>zDipole2 && zEnd<zDipole2) {
596 if (fPosition<zDipole2 && zEnd>zDipole2) {
83dbc640 597 //LinearPropagation(zDipole2-zBeg);
598 ParPropagation(zDipole2);
599 MSThin(1); // multiple scattering in the chamber
b035a148 600 WeightPropagation(zDipole2, kTRUE); // propagate weight matrix
83dbc640 601 fPosition = fPositionNew;
602 *fTrackPar = *fTrackParNew;
603 //MagnetPropagation(zEnd);
604 ParPropagation(zEnd);
b035a148 605 WeightPropagation(zEnd, kTRUE);
83dbc640 606 fPosition = fPositionNew;
607 }
608 // Check if track steps out of dipole
b035a148 609 //AZ(z->-z) else if (fPosition>zDipole1 && zEnd<zDipole1) {
610 else if (fPosition<zDipole1 && zEnd>zDipole1) {
83dbc640 611 //MagnetPropagation(zDipole1-zBeg);
612 ParPropagation(zDipole1);
613 MSThin(1); // multiple scattering in the chamber
b035a148 614 WeightPropagation(zDipole1, kTRUE);
83dbc640 615 fPosition = fPositionNew;
616 *fTrackPar = *fTrackParNew;
617 //LinearPropagation(zEnd-zDipole1);
618 ParPropagation(zEnd);
b035a148 619 WeightPropagation(zEnd, kTRUE);
83dbc640 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
b035a148 625 WeightPropagation(zEnd, kTRUE);
83dbc640 626 fPosition = fPositionNew;
627 }
628
629 // Add measurement
630 if (fRecover != 0 && hitAdd == fSkipHit && !miss) {
631 // recovered track - remove the hit
632 miss = kTRUE;
b035a148 633 ichamb = fSkipHit->GetChamberNumber();
634 // remove the skipped hit
1ffe5558 635 fTrackHits->Remove(fSkipHit);
636 fNmbTrackHits --;
637 //AZ fSkipHit->SetNTrackHits(fSkipHit->GetNTrackHits()-1); // unmark hit
83dbc640 638 Back = kFALSE;
1ffe5558 639 fRecover = 0;
640 ichambOK = ((AliMUONHitForRec*)((*fTrackHits)[fNmbTrackHits-1]))->GetChamberNumber();
cc87ebcd 641 //? if (fgTrackReconstructor->GetTrackMethod() != 3) currIndx = fgHitForRec->IndexOf(fSkipHit);
83dbc640 642 currIndx = fgHitForRec->IndexOf(fSkipHit);
643 }
644
645 if (Back && !miss) {
646 // backward propagator
b035a148 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;
b035a148 658 if (fTrackDir > 0) AddMatrices (this, dChi2, hitAdd);
659 fChi2 += dChi2; // Chi2
660 } else *fTrackPar = *fTrackParNew; // adjust end point to the last hit
83dbc640 661 if (ichamb==ichamEnd) break;
b035a148 662 currIndx ++;
83dbc640 663 } else {
664 // forward propagator
cc87ebcd 665 if (miss || !FindPoint(ichamb, zEnd, currIndx, iFB, hitAdd, iz)) {
83dbc640 666 // missing point
667 *fTrackPar = *fTrackParNew;
668 } else {
669 //add point
1ffe5558 670 fTrackHits->Add(hitAdd); // add hit
671 fNmbTrackHits ++;
83dbc640 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
1ffe5558 678 if (fgDebug > 0) cout << fNmbTrackHits << " " << fChi2 << " " << 1/(*fTrackPar)(4,0) << " " << fPosition << endl;
29f1b13a 679 if (1/TMath::Abs((*fTrackPar)(4,0)) < fgTrackReconstructor->GetMinBendingMomentum()) success = kFALSE; // p < p_min
4fcb5a03 680 if (GetRecover() < 0) success = kFALSE;
83dbc640 681 return success;
682}
683
684 //__________________________________________________________________________
685void AliMUONTrackK::ParPropagation(Double_t zEnd)
686{
d19b6003 687/// Propagation of track parameters to zEnd
83dbc640 688 Int_t iFB, nTries;
689 Double_t dZ, step, distance, charge;
690 Double_t vGeant3[7], vGeant3New[7];
b21aa6eb 691 AliMUONTrackParam trackParam;
692
83dbc640 693 nTries = 0;
694 // First step using linear extrapolation
695 dZ = zEnd - fPosition;
83dbc640 696 fPositionNew = fPosition;
697 *fTrackParNew = *fTrackPar;
b035a148 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));
83dbc640 703 SetGeantParam(vGeant3,iFB);
b21aa6eb 704 //fTrackParNew->Print();
83dbc640 705
706 // Check if overstep
707 do {
708 step = TMath::Abs(step);
709 // Propagate parameters
37827b29 710 AliMUONTrackExtrap::ExtrapOneStepRungekutta(charge,step,vGeant3,vGeant3New);
b21aa6eb 711 //extrap_onestep_rungekutta(charge,step,vGeant3,vGeant3New);
83dbc640 712 distance = zEnd - vGeant3New[2];
713 step *= dZ/(vGeant3New[2]-fPositionNew);
714 nTries ++;
343146bf 715 } while (distance*iFB < 0 && TMath::Abs(distance) > fgkEpsilon);
83dbc640 716
717 GetFromGeantParam(vGeant3New,iFB);
b21aa6eb 718 //fTrackParNew->Print();
83dbc640 719
b035a148 720 // Position adjustment (until within tolerance)
343146bf 721 while (TMath::Abs(distance) > fgkEpsilon) {
83dbc640 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
37827b29 730 AliMUONTrackExtrap::ExtrapOneStepRungekutta(charge,step,vGeant3,vGeant3New);
b21aa6eb 731 //extrap_onestep_rungekutta(charge,step,vGeant3,vGeant3New);
83dbc640 732 distance = zEnd - vGeant3New[2];
733 step /= 2;
734 nTries ++;
94eb555e 735 if (nTries > fgkTriesMax) AliError(Form(" Too many tries: %d", nTries));
83dbc640 736 } while (distance*iFB < 0);
737
738 GetFromGeantParam(vGeant3New,iFB);
739 }
740 //cout << nTries << endl;
b21aa6eb 741 //fTrackParNew->Print();
83dbc640 742 return;
743}
83dbc640 744
83dbc640 745 //__________________________________________________________________________
b035a148 746void AliMUONTrackK::WeightPropagation(Double_t zEnd, Bool_t smooth)
83dbc640 747{
d19b6003 748/// Propagation of the weight matrix
749/// W = DtWD, where D is Jacobian
83dbc640 750 Int_t i, j;
751 Double_t dPar;
752
30178c30 753 TMatrixD jacob(fgkSize,fgkSize);
754 jacob = 0;
83dbc640 755
756 // Save initial and propagated parameters
757 TMatrixD trackPar0 = *fTrackPar;
758 TMatrixD trackParNew0 = *fTrackParNew;
83dbc640 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) {
b21aa6eb 765 Int_t ifail;
766 mnvertLocal(&((*fCovariance)(0,0)), fgkSize,fgkSize,fgkSize,ifail);
767 //fCovariance->Print();
83dbc640 768 } else {
8c343c7c 769 AliWarning(" Determinant fCovariance=0:");
83dbc640 770 }
771
b21aa6eb 772 // Loop over parameters to find change of the propagated vs initial ones
343146bf 773 for (i=0; i<fgkSize; i++) {
83dbc640 774 dPar = TMath::Sqrt((*fCovariance)(i,i));
b21aa6eb 775 *fTrackPar = trackPar0;
776 if (i == 4) dPar *= TMath::Sign(1.,-trackPar0(4,0)); // 1/p
83dbc640 777 (*fTrackPar)(i,0) += dPar;
778 ParPropagation(zEnd);
343146bf 779 for (j=0; j<fgkSize; j++) {
b21aa6eb 780 jacob(j,i) = ((*fTrackParNew)(j,0)-trackParNew0(j,0))/dPar;
83dbc640 781 }
782 }
783
83dbc640 784 //trackParNew0.Print();
30178c30 785 //TMatrixD par1(jacob,TMatrixD::kMult,trackPar0); //
83dbc640 786 //par1.Print();
b21aa6eb 787 TMatrixD jacob0 = jacob;
30178c30 788 if (jacob.Determinant() != 0) {
b21aa6eb 789 jacob.Invert();
83dbc640 790 } else {
b21aa6eb 791 AliWarning(" Determinant jacob=0:");
83dbc640 792 }
30178c30 793 TMatrixD weight1(*fWeight,TMatrixD::kMult,jacob); // WD
794 *fWeight = TMatrixD(jacob,TMatrixD::kTransposeMult,weight1); // DtWD
83dbc640 795
796 // Restore initial and propagated parameters
797 *fTrackPar = trackPar0;
798 *fTrackParNew = trackParNew0;
b035a148 799
800 // Save for smoother
801 if (!smooth) return; // do not use smoother
b035a148 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) {
b21aa6eb 811 Int_t ifail;
812 mnvertLocal(&((*fCovariance)(0,0)), fgkSize,fgkSize,fgkSize,ifail);
b035a148 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);
b21aa6eb 822 tmp = new TMatrixD(jacob0); // Jacobian
b035a148 823 fJacob->Add(tmp);
824 tmp->SetUniqueID(1);
825 if (fSteps->fN <= fNSteps) fSteps->Set(fSteps->fN+10);
826 fSteps->AddAt(fPositionNew,fNSteps++);
c0e54947 827 if (fgDebug > 0) printf(" WeightPropagation %d %.3f %.3f %.3f \n", fNSteps,
828 (*fTrackParNew)(1,0), (*fTrackParNew)(0,0), fPositionNew);
83dbc640 829 return;
830}
831
832 //__________________________________________________________________________
cc87ebcd 833Bool_t AliMUONTrackK::FindPoint(Int_t ichamb, Double_t zEnd, Int_t currIndx, Int_t iFB, AliMUONHitForRec *&hitAdd, Int_t iz)
83dbc640 834{
d19b6003 835/// Picks up point within a window for the chamber No ichamb
836/// Split the track if there are more than 1 hit
4fcb5a03 837 Int_t ihit, nRecTracks;
30178c30 838 Double_t windowB, windowNonB, dChi2Tmp=0, dChi2, y, x, savePosition=0;
83dbc640 839 TClonesArray *trackPtr;
840 AliMUONHitForRec *hit, *hitLoop;
841 AliMUONTrackK *trackK;
cc87ebcd 842 AliMUONDetElement *detElem = NULL;
83dbc640 843
c0e54947 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);
30178c30 857 Bool_t ok = kFALSE;
cc87ebcd 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;
c0e54947 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.);
cc87ebcd 870 for (Int_t i = 0; i < nDetElem; i++) {
871 detElem = fgCombi->DetElem(pDEatZ[i]);
c0e54947 872 if (detElem->Inside((*fTrackParNew)(1,0), (*fTrackParNew)(0,0), fPosition, windowNonB, windowB)) {
873 detElem->ClusterReco((*fTrackParNew)(1,0), (*fTrackParNew)(0,0), windowNonB, windowB);
cc87ebcd 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
29f1b13a 883 //sigmaB = fgTrackReconstructor->GetBendingResolution(); // bending resolution
884 //sigmaNonB = fgTrackReconstructor->GetNonBendingResolution(); // non-bending resolution
83dbc640 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) {
b21aa6eb 889 Int_t ifail;
890 mnvertLocal(&((*fCovariance)(0,0)), fgkSize,fgkSize,fgkSize,ifail);
83dbc640 891 } else {
b21aa6eb 892 AliWarning(" Determinant fCovariance=0:");
83dbc640 893 }
30178c30 894 //windowB = fgkNSigma*TMath::Sqrt((*fCovariance)(0,0)+sigmaB*sigmaB);
895 //windowNonB = fgkNSigma*TMath::Sqrt((*fCovariance)(1,1)+sigmaNonB*sigmaNonB);
83dbc640 896 // Loop over all hits and take hits from the chamber
343146bf 897 TMatrixD pointWeight(fgkSize,fgkSize);
83dbc640 898 TMatrixD saveWeight = pointWeight;
899 TMatrixD pointWeightTmp = pointWeight;
343146bf 900 TMatrixD point(fgkSize,1);
83dbc640 901 TMatrixD trackPar = point;
902 TMatrixD trackParTmp = point;
cc87ebcd 903 Int_t nHitsOK = 0, ihitB = currIndx, ihitE = fgNOfPoints, iDhit = iFB;
b035a148 904 Double_t zLast;
1ffe5558 905 zLast = ((AliMUONHitForRec*)fTrackHits->Last())->GetZ();
cc87ebcd 906 if (fgTrackReconstructor->GetTrackMethod() == 3 && detElem) {
907 ihitB = 0;
908 ihitE = detElem->NHitsForRec();
909 iDhit = 1;
910 }
911
4fcb5a03 912 TArrayD branchChi2(20);
cc87ebcd 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);
83dbc640 917 if (hit->GetChamberNumber() == ichamb) {
918 //if (TMath::Abs(hit->GetZ()-zEnd) < 0.1) {
919 if (TMath::Abs(hit->GetZ()-zEnd) < 0.5) {
83dbc640 920 y = hit->GetBendingCoor();
921 x = hit->GetNonBendingCoor();
cc87ebcd 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 }
30178c30 929 windowB = fgkNSigma*TMath::Sqrt((*fCovariance)(0,0)+hit->GetBendingReso2());
930 windowNonB = fgkNSigma*TMath::Sqrt((*fCovariance)(1,1)+hit->GetNonBendingReso2());
b035a148 931
b21aa6eb 932 // windowB = TMath::Min (windowB,5.);
933 if (fgkNSigma > 6) windowB = TMath::Min (windowB,5.);
b035a148 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;
30178c30 959 if (TMath::Abs((*fTrackParNew)(0,0)-y) <= windowB &&
960 TMath::Abs((*fTrackParNew)(1,0)-x) <= windowNonB) {
b035a148 961 //if (TMath::Abs((*fTrackParNew)(0,0)-y) <= windowB &&
962 // TMath::Abs((*fTrackParNew)(1,0)-x) <= windowNonB &&
29fc2c86 963 // hit->GetTrackRefSignal() == 1) { // just for test
83dbc640 964 // Vector of measurements and covariance matrix
4c08623d 965 //fprintf(lun1,"%3d %3d %10.4f %10.4f \n", AliRunLoader::GetRunLoader()->GetEventNumber(), ichamb, x, y);
96a96c2b 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 }
83dbc640 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);
29f1b13a 990 if (TMath::Abs(1./(trackPar)(4,0)) < fgTrackReconstructor->GetMinBendingMomentum()) continue; // p < p_min - next hit
cc87ebcd 991 // if (TMath::Sign(1.,(trackPar)(4,0)*(*fTrackPar)(4,0)) < 0) continue; // change of sign
30178c30 992 ok = kTRUE;
83dbc640 993 nHitsOK++;
994 //if (nHitsOK > -1) {
4fcb5a03 995 if (nHitsOK == 1) {
83dbc640 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;
c0e54947 1004 if (fgDebug > 0) printf(" Added point (ch, x, y, Chi2): %d %.3f %.3f %.3f\n", ichamb, x, y, dChi2);
4fcb5a03 1005 branchChi2[0] = dChi2;
83dbc640 1006 } else {
1007 // branching: create a new track
29f1b13a 1008 trackPtr = fgTrackReconstructor->GetRecTracksPtr();
1009 nRecTracks = fgTrackReconstructor->GetNRecTracks();
b035a148 1010 trackK = new ((*trackPtr)[nRecTracks]) AliMUONTrackK(NULL, NULL);
83dbc640 1011 *trackK = *this;
29f1b13a 1012 fgTrackReconstructor->SetNRecTracks(nRecTracks+1);
c0e54947 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);
83dbc640 1014 trackK->fRecover = 0;
1015 *(trackK->fTrackPar) = trackPar;
1016 *(trackK->fWeight) += pointWeight;
1017 trackK->fChi2 += dChi2;
b035a148 1018 // Check
1019 /*
1020 *(trackK->fCovariance) = *(trackK->fWeight);
1021 if (trackK->fCovariance->Determinant() != 0) {
b21aa6eb 1022 Int_t ifail;
1023 mnvertLocal(&((*(trackK->fCovariance))(0,0)), fgkSize,fgkSize,fgkSize,ifail);
b035a148 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
b035a148 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);
94eb555e 1034 AliError(Form(" *** Position adjustment 1: %f %f", hit->GetZ(), (*trackK->fSteps)[i]));
b035a148 1035 }
1036 else break;
1037 }
1038 }
1039 AddMatrices (trackK, dChi2, hit);
1040 }
83dbc640 1041 // Mark hits as being on 2 tracks
1ffe5558 1042 for (Int_t i=0; i<fNmbTrackHits; i++) {
1043 hitLoop = (AliMUONHitForRec*) ((*fTrackHits)[i]);
1044 //AZ hitLoop->SetNTrackHits(hitLoop->GetNTrackHits()+1);
b035a148 1045 if (fgDebug >=10) {
1046 cout << " ** ";
1047 cout << hitLoop->GetChamberNumber() << " ";
1048 cout << hitLoop->GetBendingCoor() << " ";
1049 cout << hitLoop->GetNonBendingCoor() << " ";
1050 cout << hitLoop->GetZ() << " " << " ";
29fc2c86 1051 cout << hitLoop->GetTTRTrack() << endl;
8f373020 1052 printf(" ** %d %10.4f %10.4f %10.4f\n",
b035a148 1053 hitLoop->GetChamberNumber(), hitLoop->GetBendingCoor(),
8f373020 1054 hitLoop->GetNonBendingCoor(), hitLoop->GetZ());
b035a148 1055 }
83dbc640 1056 }
1057 //add point
1ffe5558 1058 trackK->fTrackHits->Add(hit); // add hit
1059 trackK->fNmbTrackHits ++;
83dbc640 1060 hit->SetNTrackHits(hit->GetNTrackHits()+1); // mark hit as being on track
1061 if (ichamb == 9) {
1062 // the last chamber
b035a148 1063 trackK->fTrackDir = 1;
83dbc640 1064 trackK->fBPFlag = kTRUE;
1065 }
4fcb5a03 1066 if (nHitsOK > branchChi2.GetSize()) branchChi2.Set(branchChi2.GetSize()+10);
1067 branchChi2[nHitsOK-1] = dChi2;
83dbc640 1068 }
1069 }
1070 }
1071 } else break; // different chamber
1072 } // for (ihit=currIndx;
30178c30 1073 if (ok) {
b035a148 1074 // Restore members
83dbc640 1075 *fTrackPar = trackParTmp;
1076 *fWeight = saveWeight;
1077 *fWeight += pointWeightTmp;
1078 fChi2 += dChi2Tmp; // Chi2
83dbc640 1079 fPosition = savePosition;
b035a148 1080 // Add filtered matrices for the smoother
b035a148 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);
b21aa6eb 1086 if (fgDebug > 0) cout << " *** Position adjustment 2 " << fPosition << " " << (*fSteps)[i] << endl;
b035a148 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)
4fcb5a03 1098 // Check for maximum number of branches - exclude excessive
1099 if (nHitsOK > 1) CheckBranches(branchChi2, nHitsOK);
83dbc640 1100 }
30178c30 1101 return ok;
83dbc640 1102}
1103
4fcb5a03 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
83dbc640 1133 //__________________________________________________________________________
1134void AliMUONTrackK::TryPoint(TMatrixD &point, const TMatrixD &pointWeight, TMatrixD &trackParTmp, Double_t &dChi2)
1135{
d19b6003 1136/// Adds a measurement point (modifies track parameters and computes
1137/// change of Chi2)
83dbc640 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) {
cc87ebcd 1151 Int_t ifail;
1152 mnvertLocal(&((wu)(0,0)), fgkSize,fgkSize,fgkSize,ifail);
83dbc640 1153 } else {
b21aa6eb 1154 AliWarning(" Determinant wu=0:");
83dbc640 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{
d19b6003 1174/// Adds multiple scattering in a thin layer (only angles are affected)
83dbc640 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) {
b21aa6eb 1180 Int_t ifail;
1181 mnvertLocal(&((*fWeight)(0,0)), fgkSize,fgkSize,fgkSize,ifail);
83dbc640 1182 } else {
b21aa6eb 1183 AliWarning(" Determinant fWeight=0:");
83dbc640 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
208f139e 1191 path = TMath::Abs(AliMUONConstants::ChamberThicknessInX0()/cosAlph/cosBeta); // path length
83dbc640 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
b21aa6eb 1196 Int_t ifail;
1197 mnvertLocal(&((*fWeight)(0,0)), fgkSize,fgkSize,fgkSize,ifail);
83dbc640 1198 return;
1199}
b035a148 1200
83dbc640 1201 //__________________________________________________________________________
1202void AliMUONTrackK::StartBack(void)
1203{
d19b6003 1204/// Starts backpropagator
83dbc640 1205
1206 fBPFlag = kTRUE;
1207 fChi2 = 0;
343146bf 1208 for (Int_t i=0; i<fgkSize; i++) {
1209 for (Int_t j=0; j<fgkSize; j++) {
83dbc640 1210 if (j==i) (*fWeight)(i,i) /= 100;
1ffe5558 1211 //if (j==i) (*fWeight)(i,i) /= fNmbTrackHits*fNmbTrackHits;
83dbc640 1212 else (*fWeight)(j,i) = 0;
1213 }
1214 }
b035a148 1215 // Sort hits on track in descending order in abs(z)
1ffe5558 1216 SortHits(0, fTrackHits);
b035a148 1217}
1218
1219 //__________________________________________________________________________
1220void AliMUONTrackK::SortHits(Int_t iflag, TObjArray *array)
1221{
208f139e 1222/// Sort hits in Z if the seed segment is in the last but one station
d19b6003 1223/// (if iflag==0 in descending order in abs(z), if !=0 - unsort)
b035a148 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 }
83dbc640 1249}
1250
1251 //__________________________________________________________________________
1252void AliMUONTrackK::SetGeantParam(Double_t *VGeant3, Int_t iFB)
1253{
d19b6003 1254/// Set vector of Geant3 parameters pointed to by "VGeant3"
1255/// from track parameters
83dbc640 1256
1257 VGeant3[0] = (*fTrackParNew)(1,0); // X
1258 VGeant3[1] = (*fTrackParNew)(0,0); // Y
1259 VGeant3[2] = fPositionNew; // Z
b21aa6eb 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
83dbc640 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{
d19b6003 1269/// Get track parameters from vector of Geant3 parameters pointed
1270/// to by "VGeant3"
83dbc640 1271
1272 fPositionNew = VGeant3[2]; // Z
1273 (*fTrackParNew)(0,0) = VGeant3[1]; // Y
1274 (*fTrackParNew)(1,0) = VGeant3[0]; // X
b21aa6eb 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
83dbc640 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{
d19b6003 1283/// Computes "track quality" from Chi2 (if iChi2==0) or vice versa
83dbc640 1284
b035a148 1285 if (fChi2 > 500) {
94eb555e 1286 AliWarning(Form(" *** Too high Chi2: %f ", fChi2));
b035a148 1287 fChi2 = 500;
83dbc640 1288 }
1ffe5558 1289 if (iChi2 == 0) fChi2 = fNmbTrackHits + (500.-fChi2)/501;
1290 else fChi2 = 500 - (fChi2-fNmbTrackHits)*501;
83dbc640 1291}
1292
1293 //__________________________________________________________________________
1294Int_t AliMUONTrackK::Compare(const TObject* trackK) const
1295{
d19b6003 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
83dbc640 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 //__________________________________________________________________________
30178c30 1306Bool_t AliMUONTrackK::KeepTrack(AliMUONTrackK* track0) const
83dbc640 1307{
d19b6003 1308/// Check whether or not to keep current track
1309/// (keep, if it has less than half of common hits with track0)
83dbc640 1310 Int_t hitsInCommon, nHits0, i, j, nTrackHits2;
1311 AliMUONHitForRec *hit0, *hit1;
1312
1313 hitsInCommon = 0;
1ffe5558 1314 nHits0 = track0->fNmbTrackHits;
1315 nTrackHits2 = fNmbTrackHits/2;
83dbc640 1316
1317 for (i=0; i<nHits0; i++) {
1318 // Check if hit belongs to several tracks
1ffe5558 1319 hit0 = (AliMUONHitForRec*) (*track0->fTrackHits)[i];
83dbc640 1320 if (hit0->GetNTrackHits() == 1) continue;
1ffe5558 1321 for (j=0; j<fNmbTrackHits; j++) {
1322 hit1 = (AliMUONHitForRec*) (*fTrackHits)[j];
83dbc640 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{
d19b6003 1337/// Kill track candidate
29f1b13a 1338 fgTrackReconstructor->GetRecTracksPtr()->Remove(this);
83dbc640 1339}
1340
b035a148 1341 //__________________________________________________________________________
1342void AliMUONTrackK::FillMUONTrack(void)
1343{
d19b6003 1344/// Compute track parameters at hit positions (as for AliMUONTrack)
b035a148 1345
1ffe5558 1346 // Set Chi2
8cde4af5 1347 SetTrackQuality(1); // compute Chi2
1ffe5558 1348 SetFitFMin(fChi2);
b035a148 1349
8cde4af5 1350 // Set track parameters at hits
b035a148 1351 AliMUONTrackParam trackParam;
1352 SetTrackParam(&trackParam, fTrackPar, fPosition);
1ffe5558 1353 for (Int_t i = fNmbTrackHits-1; i>=0; i--) {
b035a148 1354 if ((*fChi2Smooth)[i] < 0) {
1355 // Propagate through last chambers
37827b29 1356 AliMUONTrackExtrap::ExtrapToZ(&trackParam, ((AliMUONHitForRec*)((*fTrackHits)[i]))->GetZ());
b035a148 1357 } else {
1358 // Take saved info
1ffe5558 1359 SetTrackParam(&trackParam, (TMatrixD*)fParSmooth->UncheckedAt(i), ((AliMUONHitForRec*)((*fTrackHits)[i]))->GetZ());
b035a148 1360 }
de2cd600 1361 AddTrackParamAtHit(&trackParam,(AliMUONHitForRec*)fTrackHits->UncheckedAt(i));
b21aa6eb 1362 // Fill array of HitForRec's
1ffe5558 1363 AddHitForRecAtHit((AliMUONHitForRec*)fTrackHits->UncheckedAt(i));
b035a148 1364 }
1365}
1366
1367 //__________________________________________________________________________
1368void AliMUONTrackK::SetTrackParam(AliMUONTrackParam *trackParam, TMatrixD *par, Double_t z)
1369{
d19b6003 1370/// Fill AliMUONTrackParam object
b035a148 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
83dbc640 1380 //__________________________________________________________________________
30178c30 1381void AliMUONTrackK::MSLine(Double_t dZ, Double_t x0)
83dbc640 1382{
d19b6003 1383/// Adds multiple scattering in a thick layer for linear propagation
83dbc640 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
b035a148 1393 Double_t step = TMath::Abs(dZ/cosAlph/cosBeta); // step length
83dbc640 1394
1395 // Projected scattering angle
30178c30 1396 Double_t theta0 = 0.0136/velo/momentum/TMath::Sqrt(x0)*(1+0.038*TMath::Log(step/x0));
83dbc640 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();
b21aa6eb 1414 Int_t ifail;
1415 mnvertLocal(&((*fCovariance)(0,0)), fgkSize,fgkSize,fgkSize,ifail);
83dbc640 1416 } else {
b21aa6eb 1417 AliWarning(" Determinant fCovariance=0:" );
83dbc640 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
b21aa6eb 1434 (*fCovariance)(1,2) += dl2*(-dXdT*dAdT+dXdB*dAdB); // <xa>
83dbc640 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) {
b21aa6eb 1446 Int_t ifail;
1447 mnvertLocal(&((*fWeight)(0,0)), fgkSize,fgkSize,fgkSize,ifail);
83dbc640 1448 } else {
b21aa6eb 1449 AliWarning(" Determinant fWeight=0:");
83dbc640 1450 }
1451}
1452
1453 //__________________________________________________________________________
b035a148 1454Bool_t AliMUONTrackK::Recover(void)
83dbc640 1455{
d19b6003 1456/// Adds new failed track(s) which can be tried to be recovered
b035a148 1457 Int_t nRecTracks;
83dbc640 1458 TClonesArray *trackPtr;
1459 AliMUONTrackK *trackK;
1460
b035a148 1461 if (fgDebug > 0) cout << " ******** Enter Recover " << endl;
29f1b13a 1462 trackPtr = fgTrackReconstructor->GetRecTracksPtr();
83dbc640 1463
b035a148 1464 // Remove hit with the highest chi2
1465 Double_t chi2 = 0;
1466 if (fgDebug > 0) {
1ffe5558 1467 for (Int_t i=0; i<fNmbTrackHits; i++) {
b035a148 1468 chi2 = fChi2Smooth ? (*fChi2Smooth)[i] : (*fChi2Array)[i];
1469 printf("%10.4f", chi2);
1470 }
1471 printf("\n");
1ffe5558 1472 for (Int_t i=0; i<fNmbTrackHits; i++) {
1473 printf("%10d", ((AliMUONHitForRec*)fTrackHits->UncheckedAt(i))->GetChamberNumber());
b035a148 1474 }
1475 printf("\n");
1476 }
f29ba3e1 1477 Double_t chi2max = 0;
b035a148 1478 Int_t imax = 0;
1ffe5558 1479 for (Int_t i=0; i<fNmbTrackHits; i++) {
b035a148 1480 chi2 = fChi2Smooth ? (*fChi2Smooth)[i] : (*fChi2Array)[i];
f29ba3e1 1481 if (chi2 < chi2max) continue;
1482 chi2max = chi2;
b035a148 1483 imax = i;
1484 }
f29ba3e1 1485 //if (chi2max < 10) return kFALSE; // !!!
1ffe5558 1486 //if (chi2max < 25) imax = fNmbTrackHits - 1;
1487 if (chi2max < 15) imax = fNmbTrackHits - 1; // discard the last point
b035a148 1488 // Check if the outlier is not from the seed segment
1ffe5558 1489 AliMUONHitForRec *skipHit = (AliMUONHitForRec*) fTrackHits->UncheckedAt(imax);
208f139e 1490 if (skipHit == (AliMUONHitForRec*) fStartSegment->First() || skipHit == (AliMUONHitForRec*) fStartSegment->Second()) {
b035a148 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
1ffe5558 1496 TObjArray *hits = new TObjArray(*fTrackHits);
b035a148 1497 Int_t imax0;
1498 imax0 = imax;
1499
1500 // Hits after the found one will be removed
dc9a75bb 1501 if (GetStation0() == 3 && skipHit->GetChamberNumber() >= 7) {
1ffe5558 1502 SortHits(1, fTrackHits); // unsort hits
1503 imax = fTrackHits->IndexOf(skipHit);
b035a148 1504 }
1ffe5558 1505 Int_t nTrackHits = fNmbTrackHits;
b035a148 1506 for (Int_t i=imax+1; i<nTrackHits; i++) {
1ffe5558 1507 AliMUONHitForRec *hit = (AliMUONHitForRec*) fTrackHits->UncheckedAt(i);
1508 fTrackHits->Remove(hit);
b035a148 1509 hit->SetNTrackHits(hit->GetNTrackHits()-1); // unmark hit
1ffe5558 1510 fNmbTrackHits--;
b035a148 1511 }
83dbc640 1512
1513 // Check if the track candidate doesn't exist yet
b035a148 1514 if (ExistDouble()) { delete hits; return kFALSE; }
1515
1516 //DropBranches(imax0, hits); // drop branches downstream the discarded hit
1517 delete hits;
83dbc640 1518
29f1b13a 1519 nRecTracks = fgTrackReconstructor->GetNRecTracks();
1ffe5558 1520 skipHit = (AliMUONHitForRec*) ((*fTrackHits)[fNmbTrackHits-1]);
b035a148 1521 // Remove all saved steps and smoother matrices after the skipped hit
1522 RemoveMatrices(skipHit->GetZ());
1523
208f139e 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) {
b035a148 1526 // Propagation toward high Z or skipped hit next to segment -
1527 // start track from segment
1528 trackK = new ((*trackPtr)[nRecTracks]) AliMUONTrackK(fStartSegment);
29f1b13a 1529 fgTrackReconstructor->SetNRecTracks(nRecTracks+1);
b035a148 1530 trackK->fRecover = 1;
1531 trackK->fSkipHit = skipHit;
1ffe5558 1532 trackK->fNmbTrackHits = fNmbTrackHits;
1533 delete trackK->fTrackHits; // not efficient ?
1534 trackK->fTrackHits = new TObjArray(*fTrackHits);
b035a148 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;
29f1b13a 1541 fgTrackReconstructor->SetNRecTracks(nRecTracks+1);
b035a148 1542 //AZ(z->-z) trackK->fTrackDir = -1;
1543 trackK->fTrackDir = 1;
83dbc640 1544 trackK->fRecover = 1;
b035a148 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);
b035a148 1550 }
1551 *((TMatrixD*)trackK->fParFilter->Last()) = *((TMatrixD*)fParExtrap->Last());
0ce28091 1552 trackK->fParFilter->Last()->SetUniqueID(1);
b035a148 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);
b035a148 1558 }
1559 *((TMatrixD*)trackK->fCovFilter->Last()) = *((TMatrixD*)fCovExtrap->Last());
0ce28091 1560 trackK->fCovFilter->Last()->SetUniqueID(1);
b035a148 1561 *(trackK->fWeight) = *((TMatrixD*)trackK->fCovFilter->Last());
1562 if (trackK->fWeight->Determinant() != 0) {
b21aa6eb 1563 Int_t ifail;
1564 mnvertLocal(&((*(trackK->fWeight))(0,0)), fgkSize,fgkSize,fgkSize,ifail);
83dbc640 1565 } else {
b21aa6eb 1566 AliWarning(" Determinant fWeight=0:");
83dbc640 1567 }
b035a148 1568 trackK->fPosition = trackK->fPositionNew = (*fSteps)[fNSteps-1];
1569 trackK->fChi2 = 0;
1ffe5558 1570 for (Int_t i=0; i<fNmbTrackHits-1; i++) trackK->fChi2 += (*fChi2Array)[i];
b035a148 1571 if (fgDebug > 0) cout << nRecTracks << " " << trackK->fRecover << endl;
1572 return kTRUE;
1573}
83dbc640 1574
b035a148 1575 //__________________________________________________________________________
1576void AliMUONTrackK::AddMatrices(AliMUONTrackK *trackK, Double_t dChi2, AliMUONHitForRec *hitAdd)
1577{
d19b6003 1578/// Adds matrices for the smoother and keep Chi2 for the point
1579/// Track parameters
b035a148 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) {
b21aa6eb 1599 Int_t ifail;
1600 mnvertLocal(&((*(trackK->fCovariance))(0,0)), fgkSize,fgkSize,fgkSize,ifail);
b035a148 1601 } else {
b21aa6eb 1602 AliWarning(" Determinant fCovariance=0:");
b035a148 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
1ffe5558 1614 Int_t indx = trackK->fTrackHits->IndexOf(hitAdd);
1615 if (indx < 0) indx = fNmbTrackHits;
b035a148 1616 if (trackK->fChi2Array->fN <= indx) trackK->fChi2Array->Set(indx+10);
1617 trackK->fChi2Array->AddAt(dChi2,indx);
1618}
1619
1620 //__________________________________________________________________________
1ffe5558 1621void AliMUONTrackK::CreateMatrix(TObjArray *objArray) const
b035a148 1622{
d19b6003 1623/// Create new matrix and add it to TObjArray
b035a148 1624
1625 TMatrixD *matrix = (TMatrixD*) objArray->First();
1626 TMatrixD *tmp = new TMatrixD(*matrix);
b21aa6eb 1627 objArray->AddAtAndExpand(tmp,objArray->GetLast());
b035a148 1628 //cout << " !!! " << tmp << " " << tmp->GetUniqueID() << " " << tmp->GetNoElements() << endl;
1629}
1630
1631 //__________________________________________________________________________
1632void AliMUONTrackK::RemoveMatrices(Double_t zEnd)
1633{
d19b6003 1634/// Remove matrices (and saved steps) in the smoother part with abs(z) < abs(zEnd)
b035a148 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;
83dbc640 1641}
1642
b035a148 1643 //__________________________________________________________________________
1644void AliMUONTrackK::RemoveMatrices(AliMUONTrackK* trackK)
1645{
d19b6003 1646/// Remove last saved matrices and steps in the smoother part
b035a148 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
b035a148 1685 //__________________________________________________________________________
1686Bool_t AliMUONTrackK::Smooth(void)
1687{
d19b6003 1688/// Apply smoother
1ffe5558 1689 Int_t ihit = fNmbTrackHits - 1;
1690 AliMUONHitForRec *hit = (AliMUONHitForRec*) (*fTrackHits)[ihit];
1691 fChi2Smooth = new TArrayD(fNmbTrackHits);
b035a148 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;
b21aa6eb 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();}
b035a148 1704 */
1ffe5558 1705 for (Int_t i=ihit; i>=0; i--) cout << ((AliMUONHitForRec*)(*fTrackHits)[i])->GetZ() << " ";
b035a148 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--) {
1ffe5558 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;
b035a148 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;
f29ba3e1 1724 Double_t chi2max = 0;
b035a148 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) {
b21aa6eb 1731 Int_t ifail;
1732 mnvertLocal(&(weight(0,0)), fgkSize,fgkSize,fgkSize,ifail);
b035a148 1733 } else {
b21aa6eb 1734 AliWarning(" Determinant weight=0:");
b035a148 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--) {
1ffe5558 1758 hit = (AliMUONHitForRec*) (*fTrackHits)[ihit];
b035a148 1759 if (TMath::Abs(hit->GetZ()-(*fSteps)[i-1]) < 0.1) { found = kTRUE; break; }
1ffe5558 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; }
b035a148 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) {
b21aa6eb 1785 Int_t ifail;
1786 mnvertLocal(&(tmp(0,0)), fgkSize,fgkSize,fgkSize,ifail);
b035a148 1787 } else {
b21aa6eb 1788 AliWarning(" Determinant tmp=0:");
b035a148 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);
f29ba3e1 1794 if (chi2(0,0) > chi2max) chi2max = chi2(0,0);
b035a148 1795 fChi2 += chi2(0,0);
1796 if (chi2(0,0) < 0) {
94eb555e 1797 //chi2.Print();
1798 AliError(Form(" *** chi2 < 0: %d %d ", i, iLast));
b035a148 1799 }
1800 // Save smoothed parameters
1801 TMatrixD *par = new TMatrixD(parSmooth);
b21aa6eb 1802 fParSmooth->AddAtAndExpand(par, ihit);
b035a148 1803
1804 } // for (Int_t i=iLast+1;
1805
f29ba3e1 1806 //if (chi2max > 16) {
1807 //if (chi2max > 25) {
1808 //if (chi2max > 50) {
1809 //if (chi2max > 100) {
1810 if (chi2max > fgkChi2max) {
b035a148 1811 //if (Recover()) DropBranches();
1812 //Recover();
1813 Outlier();
1814 return kFALSE;
1815 }
1816 return kTRUE;
1817}
1818
1819 //__________________________________________________________________________
1820void AliMUONTrackK::Outlier()
1821{
d19b6003 1822/// Adds new track with removed hit having the highest chi2
b035a148 1823
1824 if (fgDebug > 0) {
1825 cout << " ******** Enter Outlier " << endl;
1ffe5558 1826 for (Int_t i=0; i<fNmbTrackHits; i++) printf("%10.4f", (*fChi2Smooth)[i]);
b035a148 1827 printf("\n");
1ffe5558 1828 for (Int_t i=0; i<fNmbTrackHits; i++) {
1829 printf("%10d", ((AliMUONHitForRec*)fTrackHits->UncheckedAt(i))->GetChamberNumber());
b035a148 1830 }
1831 printf("\n");
1832 }
1833
f29ba3e1 1834 Double_t chi2max = 0;
b035a148 1835 Int_t imax = 0;
1ffe5558 1836 for (Int_t i=0; i<fNmbTrackHits; i++) {
f29ba3e1 1837 if ((*fChi2Smooth)[i] < chi2max) continue;
1838 chi2max = (*fChi2Smooth)[i];
b035a148 1839 imax = i;
1840 }
1841 // Check if the outlier is not from the seed segment
1ffe5558 1842 AliMUONHitForRec *hit = (AliMUONHitForRec*) fTrackHits->UncheckedAt(imax);
208f139e 1843 if (hit == (AliMUONHitForRec*) fStartSegment->First() || hit == (AliMUONHitForRec*) fStartSegment->Second()) return; // to be changed probably
b035a148 1844
1845 // Check for missing station
1846 Int_t ok = 1;
1847 if (imax == 0) {
1ffe5558 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--;
b035a148 1851 }
1ffe5558 1852 else if (((AliMUONHitForRec*)fTrackHits->UncheckedAt(imax-1))->GetChamberNumber()/2 - ((AliMUONHitForRec*)fTrackHits->UncheckedAt(imax+1))->GetChamberNumber()/2 > 1) ok--;
b035a148 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
29f1b13a 1862 TClonesArray *trackPtr = fgTrackReconstructor->GetRecTracksPtr();
1863 Int_t nRecTracks = fgTrackReconstructor->GetNRecTracks();
b035a148 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);
29f1b13a 1869 fgTrackReconstructor->SetNRecTracks(nRecTracks+1);
b035a148 1870 trackK->fRecover = 2;
1871 trackK->fSkipHit = hit;
1ffe5558 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);
b035a148 1886 if (fgDebug > 0) cout << nRecTracks << " Outlier" << endl;
1887 return;
1888 }
1889 trackK = new ((*trackPtr)[nRecTracks]) AliMUONTrackK(NULL, NULL);
1890 *trackK = *this;
29f1b13a 1891 fgTrackReconstructor->SetNRecTracks(nRecTracks+1);
b035a148 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);
b035a148 1899 }
1900 *((TMatrixD*)trackK->fParFilter->Last()) = *((TMatrixD*)fParExtrap->Last());
0ce28091 1901 trackK->fParFilter->Last()->SetUniqueID(1);
b035a148 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);
b035a148 1907 }
1908 *((TMatrixD*)trackK->fCovFilter->Last()) = *((TMatrixD*)fCovExtrap->Last());
0ce28091 1909 trackK->fCovFilter->Last()->SetUniqueID(1);
b035a148 1910 *(trackK->fWeight) = *((TMatrixD*)trackK->fCovFilter->Last());
1911 if (trackK->fWeight->Determinant() != 0) {
b21aa6eb 1912 Int_t ifail;
1913 mnvertLocal(&((*(trackK->fWeight))(0,0)), fgkSize,fgkSize,fgkSize,ifail);
b035a148 1914 } else {
b21aa6eb 1915 AliWarning(" Determinant fWeight=0:");
b035a148 1916 }
1917 trackK->fPosition = trackK->fPositionNew = (*fSteps)[fNSteps-1];
1918 trackK->fChi2 = 0;
1ffe5558 1919 for (Int_t i=0; i<fNmbTrackHits-1; i++) trackK->fChi2 += (*fChi2Array)[i];
b035a148 1920 if (fgDebug > 0) cout << nRecTracks << " Outlier " << trackK->fChi2 << endl;
1921}
1922
1923 //__________________________________________________________________________
1924Double_t AliMUONTrackK::GetChi2PerPoint(Int_t iPoint) const
1925{
d19b6003 1926/// Return Chi2 at point
b035a148 1927 return fChi2Smooth ? (*fChi2Smooth)[iPoint] : (*fChi2Array)[iPoint];
1928 //return 0.;
1929}
1930
1931 //__________________________________________________________________________
1932void AliMUONTrackK::Print(FILE *lun) const
1933{
d19b6003 1934/// Print out track information
b035a148 1935
1936 Int_t flag = 1;
1937 AliMUONHitForRec *hit = 0;
b035a148 1938 // from raw clusters
8f373020 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;
b035a148 1952 }
8f373020 1953 flag = 0;
1954 break;
1955
b035a148 1956 Int_t sig[2]={1,1}, tid[2]={0};
1ffe5558 1957 for (Int_t i1=0; i1<fNmbTrackHits; i1++) {
b035a148 1958 if (GetChi2PerPoint(i1) < -0.1) continue;
1ffe5558 1959 hit = (AliMUONHitForRec*) ((*fTrackHits)[i1]);
29f1b13a 1960 rawclusters = fgTrackReconstructor->GetMUONData()->RawClusters(hit->GetChamberNumber());
b035a148 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 }
4c08623d 1966 //fprintf(lun,"%3d %3d %10.4f", AliRunLoader::GetRunLoader()->GetEventNumber(), hit->GetChamberNumber(), GetChi2PerPoint(i1));
b035a148 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{
d19b6003 1981/// Drop branches downstream of the skipped hit
b035a148 1982 Int_t nRecTracks;
1983 TClonesArray *trackPtr;
1984 AliMUONTrackK *trackK;
1985
29f1b13a 1986 trackPtr = fgTrackReconstructor->GetRecTracksPtr();
1987 nRecTracks = fgTrackReconstructor->GetNRecTracks();
b035a148 1988 Int_t icand = trackPtr->IndexOf(this);
1ffe5558 1989 if (!hits) hits = fTrackHits;
b035a148 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]);
1ffe5558 1994 if (trackK->fNmbTrackHits == 2 && trackK->GetRecover() == 0) continue;
b035a148 1995 if (trackK->GetRecover() < 0) continue;
1996
1ffe5558 1997 if (trackK->fNmbTrackHits >= imax + 1) {
b035a148 1998 for (Int_t j=0; j<=imax; j++) {
1ffe5558 1999 //if (j != fNmbTrackHits-1 && (*trackK->fTrackHits)[j] != (*fTrackHits)[j]) break;
2000 if ((*trackK->fTrackHits)[j] != (*hits)[j]) break;
b035a148 2001 if (j == imax) {
1ffe5558 2002 if (hits != fTrackHits) {
b035a148 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) {
1ffe5558 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--;
b035a148 2015 //}
1ffe5558 2016 else if (((AliMUONHitForRec*)trackK->fTrackHits->UncheckedAt(imax-1))->GetChamberNumber()/2 - ((AliMUONHitForRec*)trackK->fTrackHits->UncheckedAt(imax+1))->GetChamberNumber()/2 > 1) ok--;
b035a148 2017 if (!ok) trackK->SetRecover(-1);
2018 }
2019 } // for (Int_t j=0;
2020 }
2021 } // for (Int_t i=0;
2022}
2023
2024 //__________________________________________________________________________
208f139e 2025void AliMUONTrackK::DropBranches(AliMUONObjectPair *segment)
b035a148 2026{
d19b6003 2027/// Drop all candidates with the same seed segment
b035a148 2028 Int_t nRecTracks;
2029 TClonesArray *trackPtr;
2030 AliMUONTrackK *trackK;
208f139e 2031 AliMUONObjectPair *trackKSegment;
b035a148 2032
29f1b13a 2033 trackPtr = fgTrackReconstructor->GetRecTracksPtr();
2034 nRecTracks = fgTrackReconstructor->GetNRecTracks();
b035a148 2035 Int_t icand = trackPtr->IndexOf(this);
2036
2037 for (Int_t i=icand+1; i<nRecTracks; i++) {
2038 trackK = (AliMUONTrackK*) ((*trackPtr)[i]);
208f139e 2039 trackKSegment = trackK->fStartSegment;
1ffe5558 2040 if (trackK->fNmbTrackHits == 2 && trackK->GetRecover() == 0) continue;
b035a148 2041 if (trackK->GetRecover() < 0) continue;
208f139e 2042 if (trackKSegment->First() == segment->First() &&
2043 trackKSegment->Second() == segment->Second()) trackK->SetRecover(-1);
b035a148 2044 }
2045 if (fgDebug >= 0) cout << " Drop segment " << endl;
2046}
2047
2048 //__________________________________________________________________________
2049AliMUONHitForRec* AliMUONTrackK::GetHitLastOk(void)
2050{
d19b6003 2051/// Return the hit where track stopped
b035a148 2052
1ffe5558 2053 if (!fNSteps) return (AliMUONHitForRec*)((*fTrackHits)[1]);
b035a148 2054 return fSkipHit;
2055}
2056
2057 //__________________________________________________________________________
2058Int_t AliMUONTrackK::GetStation0(void)
2059{
d19b6003 2060/// Return seed station number
208f139e 2061 return ((AliMUONHitForRec*) fStartSegment->First())->GetChamberNumber() / 2;
b035a148 2062}
2063
2064 //__________________________________________________________________________
2065Bool_t AliMUONTrackK::ExistDouble(AliMUONHitForRec *hit)
2066{
d19b6003 2067/// Check if the track will make a double after outlier removal
b035a148 2068
29f1b13a 2069 TClonesArray *trackPtr = fgTrackReconstructor->GetRecTracksPtr();
2070 Int_t nRecTracks = fgTrackReconstructor->GetNRecTracks();
1ffe5558 2071 TObjArray *hitArray = new TObjArray(*fTrackHits);
b035a148 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]);
1ffe5558 2079 if (trackK->fNmbTrackHits == 2 && trackK->GetRecover() == 0) continue;
b035a148 2080 if (trackK == this) continue;
1ffe5558 2081 if (trackK->fNmbTrackHits == fNmbTrackHits || trackK->fNmbTrackHits == fNmbTrackHits-1) {
2082 TObjArray *hits = new TObjArray(*trackK->fTrackHits);
b035a148 2083 same = kTRUE;
1ffe5558 2084 if (trackK->fNmbTrackHits == fNmbTrackHits) {
2085 for (Int_t j=0; j<fNmbTrackHits; j++) {
b035a148 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;
1ffe5558 2094 for (Int_t j=0; j<fNmbTrackHits-1; j++) {
b035a148 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;
1ffe5558 2104 for (Int_t j=0; j<fNmbTrackHits-1; j++) {
b035a148 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{
d19b6003 2120/// Check if the track will make a double after recovery
b035a148 2121
29f1b13a 2122 TClonesArray *trackPtr = fgTrackReconstructor->GetRecTracksPtr();
2123 Int_t nRecTracks = fgTrackReconstructor->GetNRecTracks();
b035a148 2124
1ffe5558 2125 TObjArray *hitArray = new TObjArray(*fTrackHits);
b035a148 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]);
1ffe5558 2132 if (trackK->fNmbTrackHits == 2 && trackK->GetRecover() == 0) continue;
b035a148 2133 if (trackK == this) continue;
2134 //AZ if (trackK->GetRecover() < 0) continue; //
1ffe5558 2135 if (trackK->fNmbTrackHits >= fNmbTrackHits) {
2136 TObjArray *hits = new TObjArray(*trackK->fTrackHits);
b035a148 2137 if (trackK->GetStation0() == 3) SortHits(0, hits); // sort
2138 //if (trackK->GetStation0() == 3) SortHits(1, hits); // unsort
1ffe5558 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) {
b035a148 2143 if (trackK->fSkipHit && TMath::Abs(((AliMUONHitForRec*)((*hitArray)[j]))->GetZ()-trackK->fSkipHit->GetZ()) < 0.5) same = kTRUE;
1ffe5558 2144 //if (trackK->fNmbTrackHits > fNmbTrackHits &&
b035a148 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}
208f139e 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