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