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