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