Updating MUON code with the new class AliMUONTrackExtrap (Philippe Pillot)
[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
315 //__________________________________________________________________________
83dbc640 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
83dbc640 671 return success;
672}
673
674 //__________________________________________________________________________
675void AliMUONTrackK::ParPropagation(Double_t zEnd)
676{
d19b6003 677/// Propagation of track parameters to zEnd
83dbc640 678 Int_t iFB, nTries;
679 Double_t dZ, step, distance, charge;
680 Double_t vGeant3[7], vGeant3New[7];
b21aa6eb 681 AliMUONTrackParam trackParam;
682
83dbc640 683 nTries = 0;
684 // First step using linear extrapolation
685 dZ = zEnd - fPosition;
83dbc640 686 fPositionNew = fPosition;
687 *fTrackParNew = *fTrackPar;
b035a148 688 if (dZ == 0) return; //AZ ???
689 iFB = (Int_t)TMath::Sign(Double_t(1.0),dZ);
690 step = dZ/TMath::Cos((*fTrackPar)(2,0))/TMath::Cos((*fTrackPar)(3,0)); // linear estimate
691 //AZ(z->-z) charge = iFB*TMath::Sign(Double_t(1.0),(*fTrackPar)(4,0));
692 charge = -iFB*TMath::Sign(Double_t(1.0),(*fTrackPar)(4,0));
83dbc640 693 SetGeantParam(vGeant3,iFB);
b21aa6eb 694 //fTrackParNew->Print();
83dbc640 695
696 // Check if overstep
697 do {
698 step = TMath::Abs(step);
699 // Propagate parameters
37827b29 700 AliMUONTrackExtrap::ExtrapOneStepRungekutta(charge,step,vGeant3,vGeant3New);
b21aa6eb 701 //extrap_onestep_rungekutta(charge,step,vGeant3,vGeant3New);
83dbc640 702 distance = zEnd - vGeant3New[2];
703 step *= dZ/(vGeant3New[2]-fPositionNew);
704 nTries ++;
343146bf 705 } while (distance*iFB < 0 && TMath::Abs(distance) > fgkEpsilon);
83dbc640 706
707 GetFromGeantParam(vGeant3New,iFB);
b21aa6eb 708 //fTrackParNew->Print();
83dbc640 709
b035a148 710 // Position adjustment (until within tolerance)
343146bf 711 while (TMath::Abs(distance) > fgkEpsilon) {
83dbc640 712 dZ = zEnd - fPositionNew;
713 iFB = (Int_t)TMath::Sign(Double_t(1.0),dZ);
714 step = dZ/TMath::Cos((*fTrackParNew)(2,0))/TMath::Cos((*fTrackParNew)(3,0));
715 step = TMath::Abs(step);
716 SetGeantParam(vGeant3,iFB);
717 do {
718 // binary search
719 // Propagate parameters
37827b29 720 AliMUONTrackExtrap::ExtrapOneStepRungekutta(charge,step,vGeant3,vGeant3New);
b21aa6eb 721 //extrap_onestep_rungekutta(charge,step,vGeant3,vGeant3New);
83dbc640 722 distance = zEnd - vGeant3New[2];
723 step /= 2;
724 nTries ++;
94eb555e 725 if (nTries > fgkTriesMax) AliError(Form(" Too many tries: %d", nTries));
83dbc640 726 } while (distance*iFB < 0);
727
728 GetFromGeantParam(vGeant3New,iFB);
729 }
730 //cout << nTries << endl;
b21aa6eb 731 //fTrackParNew->Print();
83dbc640 732 return;
733}
83dbc640 734
83dbc640 735 //__________________________________________________________________________
b035a148 736void AliMUONTrackK::WeightPropagation(Double_t zEnd, Bool_t smooth)
83dbc640 737{
d19b6003 738/// Propagation of the weight matrix
739/// W = DtWD, where D is Jacobian
83dbc640 740 Int_t i, j;
741 Double_t dPar;
742
30178c30 743 TMatrixD jacob(fgkSize,fgkSize);
744 jacob = 0;
83dbc640 745
746 // Save initial and propagated parameters
747 TMatrixD trackPar0 = *fTrackPar;
748 TMatrixD trackParNew0 = *fTrackParNew;
83dbc640 749
750 // Get covariance matrix
751 *fCovariance = *fWeight;
752 // check whether the Invert method returns flag if matrix cannot be inverted,
753 // and do not calculate the Determinant in that case !!!!
754 if (fCovariance->Determinant() != 0) {
b21aa6eb 755 Int_t ifail;
756 mnvertLocal(&((*fCovariance)(0,0)), fgkSize,fgkSize,fgkSize,ifail);
757 //fCovariance->Print();
83dbc640 758 } else {
8c343c7c 759 AliWarning(" Determinant fCovariance=0:");
83dbc640 760 }
761
b21aa6eb 762 // Loop over parameters to find change of the propagated vs initial ones
343146bf 763 for (i=0; i<fgkSize; i++) {
83dbc640 764 dPar = TMath::Sqrt((*fCovariance)(i,i));
b21aa6eb 765 *fTrackPar = trackPar0;
766 if (i == 4) dPar *= TMath::Sign(1.,-trackPar0(4,0)); // 1/p
83dbc640 767 (*fTrackPar)(i,0) += dPar;
768 ParPropagation(zEnd);
343146bf 769 for (j=0; j<fgkSize; j++) {
b21aa6eb 770 jacob(j,i) = ((*fTrackParNew)(j,0)-trackParNew0(j,0))/dPar;
83dbc640 771 }
772 }
773
83dbc640 774 //trackParNew0.Print();
30178c30 775 //TMatrixD par1(jacob,TMatrixD::kMult,trackPar0); //
83dbc640 776 //par1.Print();
b21aa6eb 777 TMatrixD jacob0 = jacob;
30178c30 778 if (jacob.Determinant() != 0) {
b21aa6eb 779 jacob.Invert();
83dbc640 780 } else {
b21aa6eb 781 AliWarning(" Determinant jacob=0:");
83dbc640 782 }
30178c30 783 TMatrixD weight1(*fWeight,TMatrixD::kMult,jacob); // WD
784 *fWeight = TMatrixD(jacob,TMatrixD::kTransposeMult,weight1); // DtWD
83dbc640 785
786 // Restore initial and propagated parameters
787 *fTrackPar = trackPar0;
788 *fTrackParNew = trackParNew0;
b035a148 789
790 // Save for smoother
791 if (!smooth) return; // do not use smoother
b035a148 792 if (fTrackDir < 0) return; // only for propagation towards int. point
793 TMatrixD *tmp = new TMatrixD(*fTrackParNew); // extrapolated parameters
794 fParExtrap->Add(tmp);
795 tmp->SetUniqueID(1);
796 tmp = new TMatrixD(*fTrackParNew); // filtered parameters (if no measurement will be found)
797 fParFilter->Add(tmp);
798 tmp->SetUniqueID(1);
799 *fCovariance = *fWeight;
800 if (fCovariance->Determinant() != 0) {
b21aa6eb 801 Int_t ifail;
802 mnvertLocal(&((*fCovariance)(0,0)), fgkSize,fgkSize,fgkSize,ifail);
b035a148 803 } else {
804 AliWarning(" Determinant fCovariance=0:");
805 }
806 tmp = new TMatrixD(*fCovariance); // extrapolated covariance
807 fCovExtrap->Add(tmp);
808 tmp->SetUniqueID(1);
809 tmp = new TMatrixD(*fCovariance); // filtered covariance (if no measurement will be found)
810 fCovFilter->Add(tmp);
811 tmp->SetUniqueID(1);
b21aa6eb 812 tmp = new TMatrixD(jacob0); // Jacobian
b035a148 813 fJacob->Add(tmp);
814 tmp->SetUniqueID(1);
815 if (fSteps->fN <= fNSteps) fSteps->Set(fSteps->fN+10);
816 fSteps->AddAt(fPositionNew,fNSteps++);
817 if (fgDebug > 0) cout << " WeightPropagation " << fNSteps << " " << fPositionNew << endl;
83dbc640 818 return;
819}
820
821 //__________________________________________________________________________
cc87ebcd 822Bool_t AliMUONTrackK::FindPoint(Int_t ichamb, Double_t zEnd, Int_t currIndx, Int_t iFB, AliMUONHitForRec *&hitAdd, Int_t iz)
83dbc640 823{
d19b6003 824/// Picks up point within a window for the chamber No ichamb
825/// Split the track if there are more than 1 hit
83dbc640 826 Int_t ihit, nRecTracks;
30178c30 827 Double_t windowB, windowNonB, dChi2Tmp=0, dChi2, y, x, savePosition=0;
83dbc640 828 TClonesArray *trackPtr;
829 AliMUONHitForRec *hit, *hitLoop;
830 AliMUONTrackK *trackK;
cc87ebcd 831 AliMUONDetElement *detElem = NULL;
83dbc640 832
30178c30 833 Bool_t ok = kFALSE;
cc87ebcd 834
835 if (fgTrackReconstructor->GetTrackMethod() == 3 && iz >= 0) {
836 // Combined cluster / track finder
837 // Check if extrapolated track passes thru det. elems. at iz
838 Int_t *pDEatZ = fgCombi->DEatZ(iz);
839 Int_t nDetElem = pDEatZ[-1];
840 //cout << fgCombi->Z(iz) << " " << nDetElem << endl;
841 for (Int_t i = 0; i < nDetElem; i++) {
842 detElem = fgCombi->DetElem(pDEatZ[i]);
843 if (detElem->Inside((*fTrackParNew)(1,0), (*fTrackParNew)(0,0), fPosition)) {
844 detElem->ClusterReco((*fTrackParNew)(1,0), (*fTrackParNew)(0,0));
845 hitAdd = (AliMUONHitForRec*) detElem->HitsForRec()->First();
846 ok = kTRUE;
847 break;
848 }
849 }
850 if (!ok) return ok; // outside det. elem.
851 ok = kFALSE;
852 }
853
29f1b13a 854 //sigmaB = fgTrackReconstructor->GetBendingResolution(); // bending resolution
855 //sigmaNonB = fgTrackReconstructor->GetNonBendingResolution(); // non-bending resolution
83dbc640 856 *fCovariance = *fWeight;
857 // check whether the Invert method returns flag if matrix cannot be inverted,
858 // and do not calculate the Determinant in that case !!!!
859 if (fCovariance->Determinant() != 0) {
b21aa6eb 860 Int_t ifail;
861 mnvertLocal(&((*fCovariance)(0,0)), fgkSize,fgkSize,fgkSize,ifail);
83dbc640 862 } else {
b21aa6eb 863 AliWarning(" Determinant fCovariance=0:");
83dbc640 864 }
30178c30 865 //windowB = fgkNSigma*TMath::Sqrt((*fCovariance)(0,0)+sigmaB*sigmaB);
866 //windowNonB = fgkNSigma*TMath::Sqrt((*fCovariance)(1,1)+sigmaNonB*sigmaNonB);
83dbc640 867 // Loop over all hits and take hits from the chamber
343146bf 868 TMatrixD pointWeight(fgkSize,fgkSize);
83dbc640 869 TMatrixD saveWeight = pointWeight;
870 TMatrixD pointWeightTmp = pointWeight;
343146bf 871 TMatrixD point(fgkSize,1);
83dbc640 872 TMatrixD trackPar = point;
873 TMatrixD trackParTmp = point;
cc87ebcd 874 Int_t nHitsOK = 0, ihitB = currIndx, ihitE = fgNOfPoints, iDhit = iFB;
b035a148 875 Double_t zLast;
1ffe5558 876 zLast = ((AliMUONHitForRec*)fTrackHits->Last())->GetZ();
cc87ebcd 877 if (fgTrackReconstructor->GetTrackMethod() == 3 && detElem) {
878 ihitB = 0;
879 ihitE = detElem->NHitsForRec();
880 iDhit = 1;
881 }
882
83dbc640 883
cc87ebcd 884 for (ihit = ihitB; ihit >= 0 && ihit < ihitE; ihit+=iDhit) {
885 if (fgTrackReconstructor->GetTrackMethod() != 3 || !detElem)
886 hit = (AliMUONHitForRec*) ((*fgHitForRec)[ihit]);
887 else hit = (AliMUONHitForRec*) detElem->HitsForRec()->UncheckedAt(ihit);
83dbc640 888 if (hit->GetChamberNumber() == ichamb) {
889 //if (TMath::Abs(hit->GetZ()-zEnd) < 0.1) {
890 if (TMath::Abs(hit->GetZ()-zEnd) < 0.5) {
83dbc640 891 y = hit->GetBendingCoor();
892 x = hit->GetNonBendingCoor();
cc87ebcd 893 if (hit->GetBendingReso2() < 0) {
894 // Combined cluster / track finder
895 hit->SetBendingReso2(fgTrackReconstructor->GetBendingResolution()*
896 fgTrackReconstructor->GetBendingResolution());
897 hit->SetNonBendingReso2(fgTrackReconstructor->GetNonBendingResolution()*
898 fgTrackReconstructor->GetNonBendingResolution());
899 }
30178c30 900 windowB = fgkNSigma*TMath::Sqrt((*fCovariance)(0,0)+hit->GetBendingReso2());
901 windowNonB = fgkNSigma*TMath::Sqrt((*fCovariance)(1,1)+hit->GetNonBendingReso2());
b035a148 902
b21aa6eb 903 // windowB = TMath::Min (windowB,5.);
904 if (fgkNSigma > 6) windowB = TMath::Min (windowB,5.);
b035a148 905 /*
906 if (TMath::Abs(hit->GetZ()-zLast) < 50) {
907 windowB = TMath::Min (windowB,0.5);
908 windowNonB = TMath::Min (windowNonB,3.);
909 } else if (TMath::Abs(hit->GetZ()-zLast) < 200) {
910 windowB = TMath::Min (windowB,1.5);
911 windowNonB = TMath::Min (windowNonB,3.);
912 } else if (TMath::Abs(hit->GetZ()-zLast) < 350) {
913 windowB = TMath::Min (windowB,4.);
914 windowNonB = TMath::Min (windowNonB,6.);
915 }
916 */
917
918 // Test
919 /*
920 if (TMath::Abs(hit->GetZ()-zLast) < 50) {
921 windowB = 5*TMath::Sqrt((*fCovariance)(0,0)+hit->GetBendingReso2());
922 windowNonB = 5*TMath::Sqrt((*fCovariance)(1,1)+hit->GetNonBendingReso2());
923 } else {
924 windowB = fgkNSigma*TMath::Sqrt((*fCovariance)(0,0)+hit->GetBendingReso2());
925 windowNonB = fgkNSigma*TMath::Sqrt((*fCovariance)(1,1)+hit->GetNonBendingReso2());
926 }
927 */
928
929 //cout << TMath::Abs((*fTrackParNew)(0,0)-y) << " " << windowB << " " << TMath::Abs((*fTrackParNew)(1,0)-x) << " " << windowNonB << " " << zEnd << endl;
30178c30 930 if (TMath::Abs((*fTrackParNew)(0,0)-y) <= windowB &&
931 TMath::Abs((*fTrackParNew)(1,0)-x) <= windowNonB) {
b035a148 932 //if (TMath::Abs((*fTrackParNew)(0,0)-y) <= windowB &&
933 // TMath::Abs((*fTrackParNew)(1,0)-x) <= windowNonB &&
29fc2c86 934 // hit->GetTrackRefSignal() == 1) { // just for test
83dbc640 935 // Vector of measurements and covariance matrix
b21aa6eb 936 //fprintf(lun1,"%3d %3d %10.4f %10.4f \n", gAlice->GetEvNumber(), ichamb, x, y);
96a96c2b 937 if (TMath::Abs(hit->GetZ()-zEnd) > 0.05) {
938 // Adjust position: for multiple hits in the chamber or misalignment (Z as a function of X or Y)
939 //AliWarning(Form(" *** adjust %f %f ", zEnd, hit->GetZ()));
940 zEnd = hit->GetZ();
941 *fTrackPar = *fTrackParNew;
942 ParPropagation(zEnd);
943 WeightPropagation(zEnd, kTRUE);
944 fPosition = fPositionNew;
945 *fTrackPar = *fTrackParNew;
946 // Get covariance
947 *fCovariance = *fWeight;
948 if (fCovariance->Determinant() != 0) {
949 Int_t ifail;
950 mnvertLocal(&((*fCovariance)(0,0)), fgkSize,fgkSize,fgkSize,ifail);
951 } else {
952 AliWarning(" Determinant fCovariance=0:" );
953 }
954 }
83dbc640 955 point.Zero();
956 point(0,0) = y;
957 point(1,0) = x;
958 pointWeight(0,0) = 1/hit->GetBendingReso2();
959 pointWeight(1,1) = 1/hit->GetNonBendingReso2();
960 TryPoint(point,pointWeight,trackPar,dChi2);
29f1b13a 961 if (TMath::Abs(1./(trackPar)(4,0)) < fgTrackReconstructor->GetMinBendingMomentum()) continue; // p < p_min - next hit
cc87ebcd 962 // if (TMath::Sign(1.,(trackPar)(4,0)*(*fTrackPar)(4,0)) < 0) continue; // change of sign
30178c30 963 ok = kTRUE;
83dbc640 964 nHitsOK++;
965 //if (nHitsOK > -1) {
966 if (nHitsOK == 1) {
967 // Save current members
968 saveWeight = *fWeight;
969 savePosition = fPosition;
970 // temporary storage for the current track
971 dChi2Tmp = dChi2;
972 trackParTmp = trackPar;
973 pointWeightTmp = pointWeight;
974 hitAdd = hit;
dc9a75bb 975 if (fgDebug > 0) cout << " Added point: " << x << " " << y << " " << dChi2 << endl;
83dbc640 976 } else {
977 // branching: create a new track
29f1b13a 978 trackPtr = fgTrackReconstructor->GetRecTracksPtr();
979 nRecTracks = fgTrackReconstructor->GetNRecTracks();
b035a148 980 trackK = new ((*trackPtr)[nRecTracks]) AliMUONTrackK(NULL, NULL);
83dbc640 981 *trackK = *this;
29f1b13a 982 fgTrackReconstructor->SetNRecTracks(nRecTracks+1);
1ffe5558 983 if (fgDebug > 0) cout << " ******** New track: " << ichamb << " " << hit->GetTTRTrack() << " " << 1/(trackPar)(4,0) << " " << hit->GetBendingCoor() << " " << hit->GetNonBendingCoor() << " " << fNmbTrackHits << " " << nRecTracks << endl;
83dbc640 984 trackK->fRecover = 0;
985 *(trackK->fTrackPar) = trackPar;
986 *(trackK->fWeight) += pointWeight;
987 trackK->fChi2 += dChi2;
b035a148 988 // Check
989 /*
990 *(trackK->fCovariance) = *(trackK->fWeight);
991 if (trackK->fCovariance->Determinant() != 0) {
b21aa6eb 992 Int_t ifail;
993 mnvertLocal(&((*(trackK->fCovariance))(0,0)), fgkSize,fgkSize,fgkSize,ifail);
b035a148 994 }
995 cout << (*(trackK->fCovariance))(0,0) << " " << (*(trackK->fCovariance))(1,1) << " " << (*fCovariance)(0,0) << " " << (*fCovariance)(1,1) << endl;
996 */
997 // Add filtered matrices for the smoother
b035a148 998 if (fTrackDir > 0) {
999 if (nHitsOK > 2) { // check for position adjustment
1000 for (Int_t i=trackK->fNSteps-1; i>=0; i--) {
1001 if (TMath::Abs(hit->GetZ()-(*trackK->fSteps)[i]) > 0.1) {
1002 //if (TMath::Abs(hit->GetZ()-(trackK->fSteps)[i]) > 0.1) {
1003 RemoveMatrices(trackK);
94eb555e 1004 AliError(Form(" *** Position adjustment 1: %f %f", hit->GetZ(), (*trackK->fSteps)[i]));
b035a148 1005 }
1006 else break;
1007 }
1008 }
1009 AddMatrices (trackK, dChi2, hit);
1010 }
83dbc640 1011 // Mark hits as being on 2 tracks
1ffe5558 1012 for (Int_t i=0; i<fNmbTrackHits; i++) {
1013 hitLoop = (AliMUONHitForRec*) ((*fTrackHits)[i]);
1014 //AZ hitLoop->SetNTrackHits(hitLoop->GetNTrackHits()+1);
b035a148 1015 if (fgDebug >=10) {
1016 cout << " ** ";
1017 cout << hitLoop->GetChamberNumber() << " ";
1018 cout << hitLoop->GetBendingCoor() << " ";
1019 cout << hitLoop->GetNonBendingCoor() << " ";
1020 cout << hitLoop->GetZ() << " " << " ";
29fc2c86 1021 cout << hitLoop->GetTTRTrack() << endl;
8f373020 1022 printf(" ** %d %10.4f %10.4f %10.4f\n",
b035a148 1023 hitLoop->GetChamberNumber(), hitLoop->GetBendingCoor(),
8f373020 1024 hitLoop->GetNonBendingCoor(), hitLoop->GetZ());
b035a148 1025 }
83dbc640 1026 }
1027 //add point
1ffe5558 1028 trackK->fTrackHits->Add(hit); // add hit
1029 trackK->fNmbTrackHits ++;
83dbc640 1030 hit->SetNTrackHits(hit->GetNTrackHits()+1); // mark hit as being on track
1031 if (ichamb == 9) {
1032 // the last chamber
b035a148 1033 trackK->fTrackDir = 1;
83dbc640 1034 trackK->fBPFlag = kTRUE;
1035 }
1036 }
1037 }
1038 }
1039 } else break; // different chamber
1040 } // for (ihit=currIndx;
30178c30 1041 if (ok) {
b035a148 1042 // Restore members
83dbc640 1043 *fTrackPar = trackParTmp;
1044 *fWeight = saveWeight;
1045 *fWeight += pointWeightTmp;
1046 fChi2 += dChi2Tmp; // Chi2
83dbc640 1047 fPosition = savePosition;
b035a148 1048 // Add filtered matrices for the smoother
b035a148 1049 if (fTrackDir > 0) {
1050 for (Int_t i=fNSteps-1; i>=0; i--) {
1051 if (TMath::Abs(fPosition-(*fSteps)[i]) > 0.1) {
1052 //if (TMath::Abs(fPosition-fSteps[i]) > 0.1) {
1053 RemoveMatrices(this);
b21aa6eb 1054 if (fgDebug > 0) cout << " *** Position adjustment 2 " << fPosition << " " << (*fSteps)[i] << endl;
b035a148 1055 }
1056 else break;
1057 } // for (Int_t i=fNSteps-1;
1058 AddMatrices (this, dChi2Tmp, hitAdd);
1059 /*
1060 if (nHitsOK > 1) {
1061 for (Int_t i=0; i<fNSteps; i++) cout << (*fSteps)[i] << " "; cout << endl;
1062 for (Int_t i=0; i<trackK->fNSteps; i++) cout << (*trackK->fSteps)[i] << " "; cout << endl;
1063 }
1064 */
1065 } // if (fTrackDir > 0)
83dbc640 1066 }
30178c30 1067 return ok;
83dbc640 1068}
1069
1070 //__________________________________________________________________________
1071void AliMUONTrackK::TryPoint(TMatrixD &point, const TMatrixD &pointWeight, TMatrixD &trackParTmp, Double_t &dChi2)
1072{
d19b6003 1073/// Adds a measurement point (modifies track parameters and computes
1074/// change of Chi2)
83dbc640 1075
1076 // Solving linear system (W+U)p' = U(m-p) + (W+U)p
1077 TMatrixD wu = *fWeight;
1078 wu += pointWeight; // W+U
1079 trackParTmp = point;
1080 trackParTmp -= *fTrackParNew; // m-p
1081 TMatrixD right(pointWeight,TMatrixD::kMult,trackParTmp); // U(m-p)
1082 TMatrixD right1(wu,TMatrixD::kMult,*fTrackParNew); // (W+U)p
1083 right += right1; // U(m-p) + (W+U)p
1084
1085 // check whether the Invert method returns flag if matrix cannot be inverted,
1086 // and do not calculate the Determinant in that case !!!!
1087 if (wu.Determinant() != 0) {
cc87ebcd 1088 Int_t ifail;
1089 mnvertLocal(&((wu)(0,0)), fgkSize,fgkSize,fgkSize,ifail);
83dbc640 1090 } else {
b21aa6eb 1091 AliWarning(" Determinant wu=0:");
83dbc640 1092 }
1093 trackParTmp = TMatrixD(wu,TMatrixD::kMult,right);
1094
1095 right1 = trackParTmp;
1096 right1 -= point; // p'-m
1097 point = trackParTmp;
1098 point -= *fTrackParNew; // p'-p
1099 right = TMatrixD(*fWeight,TMatrixD::kMult,point); // W(p'-p)
1100 TMatrixD value(point,TMatrixD::kTransposeMult,right); // (p'-p)'W(p'-p)
1101 dChi2 = value(0,0);
1102 right = TMatrixD(pointWeight,TMatrixD::kMult,right1); // U(p'-m)
1103 value = TMatrixD(right1,TMatrixD::kTransposeMult,right); // (p'-m)'U(p'-m)
1104 dChi2 += value(0,0);
1105 return;
1106}
1107
1108 //__________________________________________________________________________
1109void AliMUONTrackK::MSThin(Int_t sign)
1110{
d19b6003 1111/// Adds multiple scattering in a thin layer (only angles are affected)
83dbc640 1112 Double_t cosAlph, cosBeta, momentum, velo, path, theta0;
1113
1114 // check whether the Invert method returns flag if matrix cannot be inverted,
1115 // and do not calculate the Determinant in that case !!!!
1116 if (fWeight->Determinant() != 0) {
b21aa6eb 1117 Int_t ifail;
1118 mnvertLocal(&((*fWeight)(0,0)), fgkSize,fgkSize,fgkSize,ifail);
83dbc640 1119 } else {
b21aa6eb 1120 AliWarning(" Determinant fWeight=0:");
83dbc640 1121 }
1122
1123 cosAlph = TMath::Cos((*fTrackParNew)(2,0));
1124 cosBeta = TMath::Cos((*fTrackParNew)(3,0));
1125 momentum = 1/(*fTrackParNew)(4,0); // particle momentum
1126 //velo = momentum/TMath::Sqrt(momentum*momentum+muonMass*muonMass); // velocity/c for muon hypothesis
1127 velo = 1; // relativistic
29f1b13a 1128 path = TMath::Abs(fgTrackReconstructor->GetChamberThicknessInX0()/cosAlph/cosBeta); // path length
83dbc640 1129 theta0 = 0.0136/velo/momentum*TMath::Sqrt(path)*(1+0.038*TMath::Log(path)); // projected scattering angle
1130
1131 (*fWeight)(2,2) += sign*theta0/cosBeta*theta0/cosBeta; // alpha
1132 (*fWeight)(3,3) += sign*theta0*theta0; // beta
b21aa6eb 1133 Int_t ifail;
1134 mnvertLocal(&((*fWeight)(0,0)), fgkSize,fgkSize,fgkSize,ifail);
83dbc640 1135 return;
1136}
b035a148 1137
83dbc640 1138 //__________________________________________________________________________
1139void AliMUONTrackK::StartBack(void)
1140{
d19b6003 1141/// Starts backpropagator
83dbc640 1142
1143 fBPFlag = kTRUE;
1144 fChi2 = 0;
343146bf 1145 for (Int_t i=0; i<fgkSize; i++) {
1146 for (Int_t j=0; j<fgkSize; j++) {
83dbc640 1147 if (j==i) (*fWeight)(i,i) /= 100;
1ffe5558 1148 //if (j==i) (*fWeight)(i,i) /= fNmbTrackHits*fNmbTrackHits;
83dbc640 1149 else (*fWeight)(j,i) = 0;
1150 }
1151 }
b035a148 1152 // Sort hits on track in descending order in abs(z)
1ffe5558 1153 SortHits(0, fTrackHits);
b035a148 1154}
1155
1156 //__________________________________________________________________________
1157void AliMUONTrackK::SortHits(Int_t iflag, TObjArray *array)
1158{
d19b6003 1159/// Sort hits in Z if the seed segment in the last but one station
1160/// (if iflag==0 in descending order in abs(z), if !=0 - unsort)
b035a148 1161
1162 if (iflag && ((AliMUONHitForRec*)(array->UncheckedAt(0)))->GetChamberNumber() == 6) return;
1163 Double_t z = 0, zmax = TMath::Abs(((AliMUONHitForRec*)(array->UncheckedAt(0)))->GetZ());
1164 Int_t i = 1, entries = array->GetEntriesFast();
1165 for ( ; i<entries; i++) {
1166 if (iflag) {
1167 if (((AliMUONHitForRec*)(array->UncheckedAt(i)))->GetChamberNumber() == 6) break;
1168 } else {
1169 z = TMath::Abs(((AliMUONHitForRec*)(array->UncheckedAt(i)))->GetZ());
1170 if (z < zmax) break;
1171 zmax = z;
1172 if (fgDebug >= 10) cout << " " << zmax << " " << z << " " << i << endl;
1173 }
1174 }
1175 if (!iflag) i--;
1176 for (Int_t j=0; j<=(i-1)/2; j++) {
1177 TObject *hit = array->UncheckedAt(j);
1178 array->AddAt(array->UncheckedAt(i-j),j);
1179 array->AddAt(hit,i-j);
1180 }
1181 if (fgDebug >= 10) {
1182 for (i=0; i<entries; i++)
1183 cout << ((AliMUONHitForRec*)(array->UncheckedAt(i)))->GetZ() << " ";
1184 cout << " - Sort" << endl;
1185 }
83dbc640 1186}
1187
1188 //__________________________________________________________________________
1189void AliMUONTrackK::SetGeantParam(Double_t *VGeant3, Int_t iFB)
1190{
d19b6003 1191/// Set vector of Geant3 parameters pointed to by "VGeant3"
1192/// from track parameters
83dbc640 1193
1194 VGeant3[0] = (*fTrackParNew)(1,0); // X
1195 VGeant3[1] = (*fTrackParNew)(0,0); // Y
1196 VGeant3[2] = fPositionNew; // Z
b21aa6eb 1197 VGeant3[3] = iFB*TMath::Sin((*fTrackParNew)(3,0)); // Px/Ptot
1198 VGeant3[4] = iFB*TMath::Cos((*fTrackParNew)(3,0))*TMath::Sin((*fTrackParNew)(2,0)); // Py/Ptot
83dbc640 1199 VGeant3[5] = iFB*TMath::Sqrt(1.0-VGeant3[3]*VGeant3[3]-VGeant3[4]*VGeant3[4]); // Pz/Ptot
1200 VGeant3[6] = 1/TMath::Abs((*fTrackParNew)(4,0)); // Ptot
1201}
1202
1203 //__________________________________________________________________________
1204void AliMUONTrackK::GetFromGeantParam(Double_t *VGeant3, Int_t iFB)
1205{
d19b6003 1206/// Get track parameters from vector of Geant3 parameters pointed
1207/// to by "VGeant3"
83dbc640 1208
1209 fPositionNew = VGeant3[2]; // Z
1210 (*fTrackParNew)(0,0) = VGeant3[1]; // Y
1211 (*fTrackParNew)(1,0) = VGeant3[0]; // X
b21aa6eb 1212 (*fTrackParNew)(3,0) = TMath::ASin(iFB*VGeant3[3]); // beta
1213 (*fTrackParNew)(2,0) = TMath::ASin(iFB*VGeant3[4]/TMath::Cos((*fTrackParNew)(3,0))); // alpha
83dbc640 1214 (*fTrackParNew)(4,0) = 1/VGeant3[6]*TMath::Sign(Double_t(1.0),(*fTrackPar)(4,0)); // 1/Ptot
1215}
1216
1217 //__________________________________________________________________________
1218void AliMUONTrackK::SetTrackQuality(Int_t iChi2)
1219{
d19b6003 1220/// Computes "track quality" from Chi2 (if iChi2==0) or vice versa
83dbc640 1221
b035a148 1222 if (fChi2 > 500) {
94eb555e 1223 AliWarning(Form(" *** Too high Chi2: %f ", fChi2));
b035a148 1224 fChi2 = 500;
83dbc640 1225 }
1ffe5558 1226 if (iChi2 == 0) fChi2 = fNmbTrackHits + (500.-fChi2)/501;
1227 else fChi2 = 500 - (fChi2-fNmbTrackHits)*501;
83dbc640 1228}
1229
1230 //__________________________________________________________________________
1231Int_t AliMUONTrackK::Compare(const TObject* trackK) const
1232{
d19b6003 1233/// "Compare" function to sort with decreasing "track quality".
1234/// Returns +1 (0, -1) if quality of current track
1235/// is smaller than (equal to, larger than) quality of trackK
83dbc640 1236
1237 if (fChi2 < ((AliMUONTrackK*)trackK)->fChi2) return(+1);
1238 else if (fChi2 == ((AliMUONTrackK*)trackK)->fChi2) return(0);
1239 else return(-1);
1240}
1241
1242 //__________________________________________________________________________
30178c30 1243Bool_t AliMUONTrackK::KeepTrack(AliMUONTrackK* track0) const
83dbc640 1244{
d19b6003 1245/// Check whether or not to keep current track
1246/// (keep, if it has less than half of common hits with track0)
83dbc640 1247 Int_t hitsInCommon, nHits0, i, j, nTrackHits2;
1248 AliMUONHitForRec *hit0, *hit1;
1249
1250 hitsInCommon = 0;
1ffe5558 1251 nHits0 = track0->fNmbTrackHits;
1252 nTrackHits2 = fNmbTrackHits/2;
83dbc640 1253
1254 for (i=0; i<nHits0; i++) {
1255 // Check if hit belongs to several tracks
1ffe5558 1256 hit0 = (AliMUONHitForRec*) (*track0->fTrackHits)[i];
83dbc640 1257 if (hit0->GetNTrackHits() == 1) continue;
1ffe5558 1258 for (j=0; j<fNmbTrackHits; j++) {
1259 hit1 = (AliMUONHitForRec*) (*fTrackHits)[j];
83dbc640 1260 if (hit1->GetNTrackHits() == 1) continue;
1261 if (hit0 == hit1) {
1262 hitsInCommon++;
1263 if (hitsInCommon >= nTrackHits2) return kFALSE;
1264 break;
1265 }
1266 } // for (j=0;
1267 } // for (i=0;
1268 return kTRUE;
1269}
1270
1271 //__________________________________________________________________________
1272void AliMUONTrackK::Kill(void)
1273{
d19b6003 1274/// Kill track candidate
29f1b13a 1275 fgTrackReconstructor->GetRecTracksPtr()->Remove(this);
83dbc640 1276}
1277
1278 //__________________________________________________________________________
b035a148 1279void AliMUONTrackK::FillMUONTrack(void)
1280{
d19b6003 1281/// Compute track parameters at hit positions (as for AliMUONTrack)
b035a148 1282
1ffe5558 1283 // Set Chi2
1284 SetFitFMin(fChi2);
b035a148 1285
1286 // Set track parameters at vertex
1287 AliMUONTrackParam trackParam;
1288 SetTrackParam(&trackParam, fTrackPar, fPosition);
1ffe5558 1289 SetTrackParamAtVertex(&trackParam);
b035a148 1290
1291 // Set track parameters at hits
1ffe5558 1292 for (Int_t i = fNmbTrackHits-1; i>=0; i--) {
b035a148 1293 if ((*fChi2Smooth)[i] < 0) {
1294 // Propagate through last chambers
37827b29 1295 AliMUONTrackExtrap::ExtrapToZ(&trackParam, ((AliMUONHitForRec*)((*fTrackHits)[i]))->GetZ());
b035a148 1296 } else {
1297 // Take saved info
1ffe5558 1298 SetTrackParam(&trackParam, (TMatrixD*)fParSmooth->UncheckedAt(i), ((AliMUONHitForRec*)((*fTrackHits)[i]))->GetZ());
b035a148 1299 }
de2cd600 1300 AddTrackParamAtHit(&trackParam,(AliMUONHitForRec*)fTrackHits->UncheckedAt(i));
b21aa6eb 1301 // Fill array of HitForRec's
1ffe5558 1302 AddHitForRecAtHit((AliMUONHitForRec*)fTrackHits->UncheckedAt(i));
b035a148 1303 }
1304}
1305
1306 //__________________________________________________________________________
1307void AliMUONTrackK::SetTrackParam(AliMUONTrackParam *trackParam, TMatrixD *par, Double_t z)
1308{
d19b6003 1309/// Fill AliMUONTrackParam object
b035a148 1310
1311 trackParam->SetBendingCoor((*par)(0,0));
1312 trackParam->SetNonBendingCoor((*par)(1,0));
1313 trackParam->SetBendingSlope(TMath::Tan((*par)(2,0)));
1314 trackParam->SetNonBendingSlope(TMath::Tan((*par)(3,0))/TMath::Cos((*par)(2,0)));
1315 trackParam->SetInverseBendingMomentum((*par)(4,0)/TMath::Cos((*par)(3,0)));
1316 trackParam->SetZ(z);
1317}
1318
1319 //__________________________________________________________________________
83dbc640 1320void AliMUONTrackK::Branson(void)
1321{
d19b6003 1322/// Propagates track to the vertex thru absorber using Branson correction
1323/// (makes use of the AliMUONTrackParam class)
83dbc640 1324
b21aa6eb 1325 //AliMUONTrackParam *trackParam = new AliMUONTrackParam();
1326 AliMUONTrackParam trackParam = AliMUONTrackParam();
b035a148 1327 /*
83dbc640 1328 trackParam->SetBendingCoor((*fTrackPar)(0,0));
1329 trackParam->SetNonBendingCoor((*fTrackPar)(1,0));
1330 trackParam->SetBendingSlope(TMath::Tan((*fTrackPar)(2,0)));
1331 trackParam->SetNonBendingSlope(TMath::Tan((*fTrackPar)(3,0))/TMath::Cos((*fTrackPar)(2,0)));
1332 trackParam->SetInverseBendingMomentum((*fTrackPar)(4,0)/TMath::Cos((*fTrackPar)(3,0)));
1333 trackParam->SetZ(fPosition);
b035a148 1334 */
b21aa6eb 1335 SetTrackParam(&trackParam, fTrackPar, fPosition);
83dbc640 1336
37827b29 1337 AliMUONTrackExtrap::ExtrapToVertex(&trackParam, Double_t(0.), Double_t(0.), Double_t(0.));
83dbc640 1338
b21aa6eb 1339 (*fTrackPar)(0,0) = trackParam.GetBendingCoor();
1340 (*fTrackPar)(1,0) = trackParam.GetNonBendingCoor();
1341 (*fTrackPar)(2,0) = TMath::ATan(trackParam.GetBendingSlope());
1342 (*fTrackPar)(3,0) = TMath::ATan(TMath::Cos((*fTrackPar)(2,0))*trackParam.GetNonBendingSlope());
1343 (*fTrackPar)(4,0) = TMath::Cos((*fTrackPar)(3,0))*trackParam.GetInverseBendingMomentum();
1344 fPosition = trackParam.GetZ();
1345 //delete trackParam;
94eb555e 1346 if (fgDebug > 0) cout << 1/(*fTrackPar)(4,0) << " " << fPosition << " " << (*fTrackPar)(0,0) << endl;
83dbc640 1347
1348 // Get covariance matrix
1349 *fCovariance = *fWeight;
1350 if (fCovariance->Determinant() != 0) {
cc87ebcd 1351 Int_t ifail;
1352 mnvertLocal(&((*fCovariance)(0,0)), fgkSize,fgkSize,fgkSize,ifail);
83dbc640 1353 } else {
b21aa6eb 1354 AliWarning(" Determinant fCovariance=0:");
83dbc640 1355 }
1356}
1357
1358 //__________________________________________________________________________
1359void AliMUONTrackK::GoToZ(Double_t zEnd)
1360{
d19b6003 1361/// Propagates track to given Z
83dbc640 1362
1363 ParPropagation(zEnd);
1364 MSThin(1); // multiple scattering in the chamber
b035a148 1365 WeightPropagation(zEnd, kFALSE);
83dbc640 1366 fPosition = fPositionNew;
1367 *fTrackPar = *fTrackParNew;
1368}
1369
1370 //__________________________________________________________________________
b035a148 1371void AliMUONTrackK::GoToVertex(Int_t iflag)
83dbc640 1372{
d19b6003 1373/// Version 3.08
1374/// Propagates track to the vertex
1375/// All material constants are taken from AliRoot
83dbc640 1376
1377 static Double_t x01[5] = { 24.282, // C
1378 24.282, // C
1379 11.274, // Concrete
1380 1.758, // Fe
1381 1.758}; // Fe (cm)
1382 // inner part theta < 3 degrees
1383 static Double_t x02[5] = { 30413, // Air
1384 24.282, // C
1385 11.274, // Concrete
1386 1.758, // Fe
1387 0.369}; // W (cm)
1388 // z positions of the materials inside the absober outer part theta > 3 degres
b035a148 1389 static Double_t zPos[10] = {-90, -105, -315, -443, -468};
83dbc640 1390
30178c30 1391 Double_t dZ, r0Norm, x0, deltaP, dChi2, pTotal, pOld;
83dbc640 1392 AliMUONHitForRec *hit;
1393 AliMUONRawCluster *clus;
1394 TClonesArray *rawclusters;
1395
1396 // First step to the rear end of the absorber
b035a148 1397 Double_t zRear = -503;
83dbc640 1398 GoToZ(zRear);
1399 Double_t tan3 = TMath::Tan(3./180*TMath::Pi());
1400
1401 // Go through absorber
1402 pOld = 1/(*fTrackPar)(4,0);
1403 Double_t r0Rear = (*fTrackPar)(0,0)*(*fTrackPar)(0,0) +
1404 (*fTrackPar)(1,0)*(*fTrackPar)(1,0);
b035a148 1405 r0Rear = TMath::Sqrt(r0Rear)/TMath::Abs(fPosition)/tan3;
83dbc640 1406 r0Norm = r0Rear;
b035a148 1407 Double_t p0, cos25, cos60;
1408 if (!iflag) goto vertex;
1409
83dbc640 1410 for (Int_t i=4; i>=0; i--) {
1411 ParPropagation(zPos[i]);
b035a148 1412 WeightPropagation(zPos[i], kFALSE);
83dbc640 1413 dZ = TMath::Abs (fPositionNew-fPosition);
30178c30 1414 if (r0Norm > 1) x0 = x01[i];
1415 else x0 = x02[i];
1416 MSLine(dZ,x0); // multiple scattering in the medium (linear approximation)
83dbc640 1417 fPosition = fPositionNew;
1418 *fTrackPar = *fTrackParNew;
1419 r0Norm = (*fTrackPar)(0,0)*(*fTrackPar)(0,0) +
1420 (*fTrackPar)(1,0)*(*fTrackPar)(1,0);
b035a148 1421 r0Norm = TMath::Sqrt(r0Norm)/TMath::Abs(fPosition)/tan3;
83dbc640 1422 }
1423 // Correct momentum for energy losses
1424 pTotal = 1/TMath::Abs((*fTrackPar)(4,0));
b035a148 1425 p0 = pTotal;
1426 cos25 = TMath::Cos(2.5/180*TMath::Pi());
1427 cos60 = TMath::Cos(6.0/180*TMath::Pi());
1428 for (Int_t j=0; j<1; j++) {
83dbc640 1429 /*
1430 if (r0Rear > 1) {
1431 if (p0 < 20) {
1432 deltaP = 2.164 + 0.145e-1*p0 - 0.417e-3*p0*p0;
1433 } else {
1434 deltaP = 2.275 + 0.102e-2*p0 - 0.674e-6*p0*p0;
1435 }
1436 } else {
1437 if (p0 < 20) {
1438 deltaP = 2.581 + 0.188e-1*p0 - 0.398e-3*p0*p0;
1439 } else {
1440 deltaP = 2.727 + 0.356e-2*p0 + 0.242e-5*p0*p0;
1441 }
1442 }
1443 */
b035a148 1444 /*
83dbc640 1445 if (r0Rear < 1) {
1446 //W
1447 if (p0<15) {
1448 deltaP = 2.737 + 0.0494*p0 - 0.001123*p0*p0;
1449 } else {
1450 deltaP = 3.0643 + 0.01346*p0;
1451 }
b035a148 1452 deltaP *= 0.95;
83dbc640 1453 } else {
1454 //Pb
1455 if (p0<15) {
1456 deltaP = 2.1380 + 0.0351*p0 - 0.000853*p0*p0;
1457 } else {
1458 deltaP = 2.407 + 0.00702*p0;
1459 }
b035a148 1460 deltaP *= 0.95;
1461 }
1462 */
1463 /*
1464 if (r0Rear < 1) {
1465 //W
1466 if (p0<18) {
1467 deltaP = 2.439 + 0.806e-1*p0 - 0.500e-2*p0*p0 + 0.106e-3*p0*p0*p0;
1468 } else {
1469 deltaP = 2.767 + 0.742e-2*p0 - 0.196e-4*p0*p0 + 0.403e-7*p0*p0*p0;
1470 }
1471 //deltaP += 0.2;
1472 deltaP *= cos25;
1473 } else {
1474 //Pb
1475 if (p0<18) {
1476 deltaP = 2.209 + 0.800e-2*p0;
1477 } else {
1478 deltaP = 2.285 + 0.141e-2*p0 - 0.446e-6*p0*p0;
1479 }
1480 //deltaP += 0.2;
1481 deltaP *= cos60;
1482 }
1483 deltaP *= 1.1;
1484 */
1485 //*
1486 if (r0Rear < 1) {
1487 if (p0 < 20) {
1488 deltaP = 2.5938 + 0.0570 * p0 - 0.001151 * p0 * p0;
1489 } else {
1490 deltaP = 3.0714 + 0.011767 * p0;
1491 }
96a96c2b 1492 deltaP *= 0.75;
b035a148 1493 } else {
1494 if (p0 < 20) {
1495 deltaP = 2.1207 + 0.05478 * p0 - 0.00145079 * p0 * p0;
1496 } else {
1497 deltaP = 2.6069 + 0.0051705 * p0;
1498 }
96a96c2b 1499 deltaP *= 0.9;
83dbc640 1500 }
b035a148 1501 //*/
83dbc640 1502
b035a148 1503 p0 = pTotal + deltaP/TMath::Abs(TMath::Cos((*fTrackPar)(2,0))/TMath::Cos((*fTrackPar)(3,0)));
83dbc640 1504 }
1505 (*fTrackPar)(4,0) = 1/p0*TMath::Sign((Double_t)1.,(*fTrackPar)(4,0));
1506
1507 // Go to the vertex
b035a148 1508vertex:
83dbc640 1509 ParPropagation((Double_t)0.);
b035a148 1510 WeightPropagation((Double_t)0., kFALSE);
83dbc640 1511 fPosition = fPositionNew;
1512 //*fTrackPar = *fTrackParNew;
1513 // Add vertex as a hit
343146bf 1514 TMatrixD pointWeight(fgkSize,fgkSize);
1515 TMatrixD point(fgkSize,1);
83dbc640 1516 TMatrixD trackParTmp = point;
1517 point(0,0) = 0; // vertex coordinate - should be taken somewhere
1518 point(1,0) = 0; // vertex coordinate - should be taken somewhere
1519 pointWeight(0,0) = 1/1.e-3/1.e-3; // 10 um error
1520 pointWeight(1,1) = 1/1.e-3/1.e-3; // 10 um error
1521 TryPoint(point,pointWeight,trackParTmp,dChi2);
1522 *fTrackPar = trackParTmp;
1523 *fWeight += pointWeight;
1524 fChi2 += dChi2; // Chi2
b035a148 1525 if (fgDebug < 0) return; // no output
1526
1ffe5558 1527 cout << pOld << " " << 1/(*fTrackPar)(4,0) << " " << dChi2 << " " << fChi2 << " " << fNmbTrackHits << endl;
1528 for (Int_t i1=0; i1<fNmbTrackHits; i1++) {
1529 hit = (AliMUONHitForRec*) ((*fTrackHits)[i1]);
1530 printf ("%5d", hit->GetChamberNumber());
83dbc640 1531 }
1532 cout << endl;
b035a148 1533 if (fgDebug > 0) {
1ffe5558 1534 for (Int_t i1=0; i1<fNmbTrackHits; i1++) {
1535 hit = (AliMUONHitForRec*) ((*fTrackHits)[i1]);
1536 //cout << ((AliMUONHitForRec*)((*fTrackHits)[i1]))->GetHitNumber() << " ";
1537 //cout << ((AliMUONHitForRec*)((*fTrackHits)[i1]))->GetZ() << " ";
1538 printf ("%5d", fgHitForRec->IndexOf(hit));
b035a148 1539 }
1540 cout << endl;
83dbc640 1541 }
8f373020 1542
1543 // from raw clusters
1544 for (Int_t i1=0; i1<fNmbTrackHits; i1++) {
1545 hit = (AliMUONHitForRec*) ((*fTrackHits)[i1]);
1546 if (hit->GetHitNumber() < 0) { // combined cluster / track finder
1547 Int_t index = -hit->GetHitNumber() / 100000;
1548 Int_t iPos = -hit->GetHitNumber() - index * 100000;
1549 clus = (AliMUONRawCluster*) fgCombi->DetElem(index-1)->RawClusters()->UncheckedAt(iPos);
1550 } else {
1551 rawclusters = fgTrackReconstructor->GetMUONData()->RawClusters(hit->GetChamberNumber());
1552 clus = (AliMUONRawCluster*) rawclusters->UncheckedAt(hit->GetHitNumber());
83dbc640 1553 }
8f373020 1554 printf ("%5d", clus->GetTrack(1)%10000000);
1555
83dbc640 1556 cout << endl;
1ffe5558 1557 for (Int_t i1=0; i1<fNmbTrackHits; i1++) {
1558 hit = (AliMUONHitForRec*) ((*fTrackHits)[i1]);
cc87ebcd 1559 if (hit->GetHitNumber() < 0) { // combined cluster / track finder
1560 Int_t index = -hit->GetHitNumber() / 100000;
1561 Int_t iPos = -hit->GetHitNumber() - index * 100000;
1562 clus = (AliMUONRawCluster*) fgCombi->DetElem(index-1)->RawClusters()->UncheckedAt(iPos);
1563 } else {
1564 rawclusters = fgTrackReconstructor->GetMUONData()->RawClusters(hit->GetChamberNumber());
1565 clus = (AliMUONRawCluster*) rawclusters->UncheckedAt(hit->GetHitNumber());
1566 }
1ffe5558 1567 if (clus->GetTrack(2) != -1) printf ("%5d", clus->GetTrack(2)%10000000);
1568 else printf ("%5s", " ");
83dbc640 1569 }
1570 }
1571 cout << endl;
1ffe5558 1572 for (Int_t i1=0; i1<fNmbTrackHits; i1++) {
1573 //cout << ((AliMUONHitForRec*)((*fTrackHits)[i1]))->GetHitNumber() << " ";
1574 cout << ((AliMUONHitForRec*)((*fTrackHits)[i1]))->GetZ() << " ";
1575 //cout << fgHitForRec->IndexOf(((AliMUONHitForRec*)((*fTrackHits)[i1]))) << " ";
83dbc640 1576 }
1577 cout << endl;
1ffe5558 1578 for (Int_t i1=0; i1<fNmbTrackHits; i1++) printf("%8.4f", (*fChi2Smooth)[i1]);
b035a148 1579 cout << endl;
83dbc640 1580 cout << "---------------------------------------------------" << endl;
1581
1582 // Get covariance matrix
b035a148 1583 /* Not needed - covariance matrix is not interesting to anybody
83dbc640 1584 *fCovariance = *fWeight;
1585 if (fCovariance->Determinant() != 0) {
b21aa6eb 1586 Int_t ifail;
1587 mnvertLocal(&((*fCovariance)(0,0)), fgkSize,fgkSize,fgkSize,ifail);
83dbc640 1588 } else {
b21aa6eb 1589 AliWarning(" Determinant fCovariance=0:" );
83dbc640 1590 }
b035a148 1591 */
83dbc640 1592}
1593
1594 //__________________________________________________________________________
30178c30 1595void AliMUONTrackK::MSLine(Double_t dZ, Double_t x0)
83dbc640 1596{
d19b6003 1597/// Adds multiple scattering in a thick layer for linear propagation
83dbc640 1598
1599 Double_t cosAlph = TMath::Cos((*fTrackPar)(2,0));
1600 Double_t tanAlph = TMath::Tan((*fTrackPar)(2,0));
1601 Double_t cosBeta = TMath::Cos((*fTrackPar)(3,0));
1602 Double_t sinBeta;
1603 sinBeta = TMath::Sin((*fTrackPar)(3,0));
1604 Double_t tanBeta = TMath::Tan((*fTrackPar)(3,0));
1605 Double_t momentum = 1/(*fTrackPar)(4,0);
1606 Double_t velo = 1; // relativistic velocity
b035a148 1607 Double_t step = TMath::Abs(dZ/cosAlph/cosBeta); // step length
83dbc640 1608
1609 // Projected scattering angle
30178c30 1610 Double_t theta0 = 0.0136/velo/momentum/TMath::Sqrt(x0)*(1+0.038*TMath::Log(step/x0));
83dbc640 1611 Double_t theta02 = theta0*theta0;
1612 Double_t dl2 = step*step/2*theta02;
1613 Double_t dl3 = dl2*step*2/3;
1614
1615 //Derivatives
1616 Double_t dYdT = 1/cosAlph;
1617 Double_t dYdB = 0; //(*fTrackPar)(2,0)*sinBeta/cosAlph;
1618 Double_t dXdT = tanAlph*tanBeta;
1619 //Double_t dXdB = (1+(*fTrackPar)(2,0)*tanAlph*sinBeta*sinBeta)/cosBeta;
1620 Double_t dXdB = 1/cosBeta;
1621 Double_t dAdT = 1/cosBeta;
1622 Double_t dAdB = 0; //(*fTrackPar)(2,0)*tanBeta;
1623
1624 // Get covariance matrix
1625 *fCovariance = *fWeight;
1626 if (fCovariance->Determinant() != 0) {
1627 // fCovariance->Invert();
b21aa6eb 1628 Int_t ifail;
1629 mnvertLocal(&((*fCovariance)(0,0)), fgkSize,fgkSize,fgkSize,ifail);
83dbc640 1630 } else {
b21aa6eb 1631 AliWarning(" Determinant fCovariance=0:" );
83dbc640 1632 }
1633
1634 (*fCovariance)(0,0) += dl3*(dYdT*dYdT+dYdB*dYdB); // <yy>
1635 (*fCovariance)(1,1) += dl3*(dXdT*dXdT+dXdB*dXdB); // <xx>
1636 (*fCovariance)(2,2) += theta02*step*(dAdT*dAdT+dAdB*dAdB); // <aa>
1637 (*fCovariance)(3,3) += theta02*step; // <bb>
1638
1639 (*fCovariance)(0,1) += dl3*(dYdT*dXdT+dYdB*dXdB); // <yx>
1640 (*fCovariance)(1,0) = (*fCovariance)(0,1);
1641
1642 (*fCovariance)(0,2) += dl2*(dYdT*dAdT+dYdB*dAdB); // <ya>
1643 (*fCovariance)(2,0) = (*fCovariance)(0,2);
1644
1645 (*fCovariance)(0,3) += dl2*dYdB; // <yb>
1646 (*fCovariance)(3,0) = (*fCovariance)(0,3);
1647
b21aa6eb 1648 (*fCovariance)(1,2) += dl2*(-dXdT*dAdT+dXdB*dAdB); // <xa>
83dbc640 1649 (*fCovariance)(2,1) = (*fCovariance)(1,2);
1650
1651 (*fCovariance)(1,3) += dl2*dXdB; // <xb>
1652 (*fCovariance)(3,1) = (*fCovariance)(1,3);
1653
1654 (*fCovariance)(2,3) += theta02*step*dAdB; // <ab>
1655 (*fCovariance)(3,2) = (*fCovariance)(2,3);
1656
1657 // Get weight matrix
1658 *fWeight = *fCovariance;
1659 if (fWeight->Determinant() != 0) {
b21aa6eb 1660 Int_t ifail;
1661 mnvertLocal(&((*fWeight)(0,0)), fgkSize,fgkSize,fgkSize,ifail);
83dbc640 1662 } else {
b21aa6eb 1663 AliWarning(" Determinant fWeight=0:");
83dbc640 1664 }
1665}
1666
1667 //__________________________________________________________________________
b035a148 1668Bool_t AliMUONTrackK::Recover(void)
83dbc640 1669{
d19b6003 1670/// Adds new failed track(s) which can be tried to be recovered
b035a148 1671 Int_t nRecTracks;
83dbc640 1672 TClonesArray *trackPtr;
1673 AliMUONTrackK *trackK;
1674
b035a148 1675 if (fgDebug > 0) cout << " ******** Enter Recover " << endl;
29f1b13a 1676 trackPtr = fgTrackReconstructor->GetRecTracksPtr();
83dbc640 1677
b035a148 1678 // Remove hit with the highest chi2
1679 Double_t chi2 = 0;
1680 if (fgDebug > 0) {
1ffe5558 1681 for (Int_t i=0; i<fNmbTrackHits; i++) {
b035a148 1682 chi2 = fChi2Smooth ? (*fChi2Smooth)[i] : (*fChi2Array)[i];
1683 printf("%10.4f", chi2);
1684 }
1685 printf("\n");
1ffe5558 1686 for (Int_t i=0; i<fNmbTrackHits; i++) {
1687 printf("%10d", ((AliMUONHitForRec*)fTrackHits->UncheckedAt(i))->GetChamberNumber());
b035a148 1688 }
1689 printf("\n");
1690 }
f29ba3e1 1691 Double_t chi2max = 0;
b035a148 1692 Int_t imax = 0;
1ffe5558 1693 for (Int_t i=0; i<fNmbTrackHits; i++) {
b035a148 1694 chi2 = fChi2Smooth ? (*fChi2Smooth)[i] : (*fChi2Array)[i];
f29ba3e1 1695 if (chi2 < chi2max) continue;
1696 chi2max = chi2;
b035a148 1697 imax = i;
1698 }
f29ba3e1 1699 //if (chi2max < 10) return kFALSE; // !!!
1ffe5558 1700 //if (chi2max < 25) imax = fNmbTrackHits - 1;
1701 if (chi2max < 15) imax = fNmbTrackHits - 1; // discard the last point
b035a148 1702 // Check if the outlier is not from the seed segment
1ffe5558 1703 AliMUONHitForRec *skipHit = (AliMUONHitForRec*) fTrackHits->UncheckedAt(imax);
b035a148 1704 if (skipHit == fStartSegment->GetHitForRec1() || skipHit == fStartSegment->GetHitForRec2()) {
1705 //DropBranches(fStartSegment); // drop all tracks with the same seed segment
1706 return kFALSE; // to be changed probably
1707 }
1708
1709 // Make a copy of track hit collection
1ffe5558 1710 TObjArray *hits = new TObjArray(*fTrackHits);
b035a148 1711 Int_t imax0;
1712 imax0 = imax;
1713
1714 // Hits after the found one will be removed
dc9a75bb 1715 if (GetStation0() == 3 && skipHit->GetChamberNumber() >= 7) {
1ffe5558 1716 SortHits(1, fTrackHits); // unsort hits
1717 imax = fTrackHits->IndexOf(skipHit);
b035a148 1718 }
1ffe5558 1719 Int_t nTrackHits = fNmbTrackHits;
b035a148 1720 for (Int_t i=imax+1; i<nTrackHits; i++) {
1ffe5558 1721 AliMUONHitForRec *hit = (AliMUONHitForRec*) fTrackHits->UncheckedAt(i);
1722 fTrackHits->Remove(hit);
b035a148 1723 hit->SetNTrackHits(hit->GetNTrackHits()-1); // unmark hit
1ffe5558 1724 fNmbTrackHits--;
b035a148 1725 }
83dbc640 1726
1727 // Check if the track candidate doesn't exist yet
b035a148 1728 if (ExistDouble()) { delete hits; return kFALSE; }
1729
1730 //DropBranches(imax0, hits); // drop branches downstream the discarded hit
1731 delete hits;
83dbc640 1732
29f1b13a 1733 nRecTracks = fgTrackReconstructor->GetNRecTracks();
1ffe5558 1734 skipHit = (AliMUONHitForRec*) ((*fTrackHits)[fNmbTrackHits-1]);
b035a148 1735 // Remove all saved steps and smoother matrices after the skipped hit
1736 RemoveMatrices(skipHit->GetZ());
1737
1738 //AZ(z->-z) if (skipHit->GetZ() > fStartSegment->GetHitForRec2()->GetZ() || !fNSteps) {
1739 if (TMath::Abs(skipHit->GetZ()) > TMath::Abs(fStartSegment->GetHitForRec2()->GetZ()) || !fNSteps) {
1740 // Propagation toward high Z or skipped hit next to segment -
1741 // start track from segment
1742 trackK = new ((*trackPtr)[nRecTracks]) AliMUONTrackK(fStartSegment);
29f1b13a 1743 fgTrackReconstructor->SetNRecTracks(nRecTracks+1);
b035a148 1744 trackK->fRecover = 1;
1745 trackK->fSkipHit = skipHit;
1ffe5558 1746 trackK->fNmbTrackHits = fNmbTrackHits;
1747 delete trackK->fTrackHits; // not efficient ?
1748 trackK->fTrackHits = new TObjArray(*fTrackHits);
b035a148 1749 if (fgDebug > 0) cout << nRecTracks << " " << trackK->fRecover << endl;
1750 return kTRUE;
1751 }
1752
1753 trackK = new ((*trackPtr)[nRecTracks]) AliMUONTrackK(NULL, NULL);
1754 *trackK = *this;
29f1b13a 1755 fgTrackReconstructor->SetNRecTracks(nRecTracks+1);
b035a148 1756 //AZ(z->-z) trackK->fTrackDir = -1;
1757 trackK->fTrackDir = 1;
83dbc640 1758 trackK->fRecover = 1;
b035a148 1759 trackK->fSkipHit = skipHit;
1760 Int_t iD = trackK->fParFilter->Last()->GetUniqueID();
1761 if (iD > 1) {
1762 trackK->fParFilter->Last()->SetUniqueID(iD-1);
1763 CreateMatrix(trackK->fParFilter);
b035a148 1764 }
1765 *((TMatrixD*)trackK->fParFilter->Last()) = *((TMatrixD*)fParExtrap->Last());
0ce28091 1766 trackK->fParFilter->Last()->SetUniqueID(1);
b035a148 1767 *(trackK->fTrackPar) = *((TMatrixD*)trackK->fParFilter->Last());
1768 iD = trackK->fCovFilter->Last()->GetUniqueID();
1769 if (iD > 1) {
1770 trackK->fCovFilter->Last()->SetUniqueID(iD-1);
1771 CreateMatrix(trackK->fCovFilter);
b035a148 1772 }
1773 *((TMatrixD*)trackK->fCovFilter->Last()) = *((TMatrixD*)fCovExtrap->Last());
0ce28091 1774 trackK->fCovFilter->Last()->SetUniqueID(1);
b035a148 1775 *(trackK->fWeight) = *((TMatrixD*)trackK->fCovFilter->Last());
1776 if (trackK->fWeight->Determinant() != 0) {
b21aa6eb 1777 Int_t ifail;
1778 mnvertLocal(&((*(trackK->fWeight))(0,0)), fgkSize,fgkSize,fgkSize,ifail);
83dbc640 1779 } else {
b21aa6eb 1780 AliWarning(" Determinant fWeight=0:");
83dbc640 1781 }
b035a148 1782 trackK->fPosition = trackK->fPositionNew = (*fSteps)[fNSteps-1];
1783 trackK->fChi2 = 0;
1ffe5558 1784 for (Int_t i=0; i<fNmbTrackHits-1; i++) trackK->fChi2 += (*fChi2Array)[i];
b035a148 1785 if (fgDebug > 0) cout << nRecTracks << " " << trackK->fRecover << endl;
1786 return kTRUE;
1787}
83dbc640 1788
b035a148 1789 //__________________________________________________________________________
1790void AliMUONTrackK::AddMatrices(AliMUONTrackK *trackK, Double_t dChi2, AliMUONHitForRec *hitAdd)
1791{
d19b6003 1792/// Adds matrices for the smoother and keep Chi2 for the point
1793/// Track parameters
b035a148 1794 //trackK->fParFilter->Last()->Print();
1795 Int_t iD = trackK->fParFilter->Last()->GetUniqueID();
1796 if (iD > 1) {
1797 trackK->fParFilter->Last()->SetUniqueID(iD-1);
1798 CreateMatrix(trackK->fParFilter);
1799 iD = 1;
1800 }
1801 *((TMatrixD*)(trackK->fParFilter->Last())) = *(trackK->fTrackPar);
1802 trackK->fParFilter->Last()->SetUniqueID(iD);
1803 if (fgDebug > 1) {
1804 cout << " Add matrices" << " " << fPosition << " " << dChi2 << " " << fgHitForRec->IndexOf(hitAdd) << endl;
1805 //trackK->fTrackPar->Print();
1806 //trackK->fTrackParNew->Print();
1807 trackK->fParFilter->Last()->Print();
1808 cout << " Add matrices" << endl;
1809 }
1810 // Covariance
1811 *(trackK->fCovariance) = *(trackK->fWeight);
1812 if (trackK->fCovariance->Determinant() != 0) {
b21aa6eb 1813 Int_t ifail;
1814 mnvertLocal(&((*(trackK->fCovariance))(0,0)), fgkSize,fgkSize,fgkSize,ifail);
b035a148 1815 } else {
b21aa6eb 1816 AliWarning(" Determinant fCovariance=0:");
b035a148 1817 }
1818 iD = trackK->fCovFilter->Last()->GetUniqueID();
1819 if (iD > 1) {
1820 trackK->fCovFilter->Last()->SetUniqueID(iD-1);
1821 CreateMatrix(trackK->fCovFilter);
1822 iD = 1;
1823 }
1824 *((TMatrixD*)(trackK->fCovFilter->Last())) = *(trackK->fCovariance);
1825 trackK->fCovFilter->Last()->SetUniqueID(iD);
1826
1827 // Save Chi2-increment for point
1ffe5558 1828 Int_t indx = trackK->fTrackHits->IndexOf(hitAdd);
1829 if (indx < 0) indx = fNmbTrackHits;
b035a148 1830 if (trackK->fChi2Array->fN <= indx) trackK->fChi2Array->Set(indx+10);
1831 trackK->fChi2Array->AddAt(dChi2,indx);
1832}
1833
1834 //__________________________________________________________________________
1ffe5558 1835void AliMUONTrackK::CreateMatrix(TObjArray *objArray) const
b035a148 1836{
d19b6003 1837/// Create new matrix and add it to TObjArray
b035a148 1838
1839 TMatrixD *matrix = (TMatrixD*) objArray->First();
1840 TMatrixD *tmp = new TMatrixD(*matrix);
b21aa6eb 1841 objArray->AddAtAndExpand(tmp,objArray->GetLast());
b035a148 1842 //cout << " !!! " << tmp << " " << tmp->GetUniqueID() << " " << tmp->GetNoElements() << endl;
1843}
1844
1845 //__________________________________________________________________________
1846void AliMUONTrackK::RemoveMatrices(Double_t zEnd)
1847{
d19b6003 1848/// Remove matrices (and saved steps) in the smoother part with abs(z) < abs(zEnd)
b035a148 1849
1850 for (Int_t i=fNSteps-1; i>=0; i--) {
1851 if (fgDebug > 1) printf("%10.4f %10.4f \n", (*fSteps)[i], zEnd);
1852 if (TMath::Abs((*fSteps)[i]) > TMath::Abs(zEnd) - 0.01) break;
1853 RemoveMatrices(this);
1854 } // for (Int_t i=fNSteps-1;
83dbc640 1855}
1856
b035a148 1857 //__________________________________________________________________________
1858void AliMUONTrackK::RemoveMatrices(AliMUONTrackK* trackK)
1859{
d19b6003 1860/// Remove last saved matrices and steps in the smoother part
b035a148 1861
1862 trackK->fNSteps--;
1863 Int_t i = trackK->fNSteps;
1864
1865 Int_t id = 0;
1866 // Delete only matrices not shared by several tracks
1867 id = trackK->fParExtrap->Last()->GetUniqueID();
1868 if (id > 1) {
1869 trackK->fParExtrap->Last()->SetUniqueID(id-1);
1870 trackK->fParExtrap->RemoveAt(i);
1871 }
1872 else ((TMatrixD*)(trackK->fParExtrap->RemoveAt(i)))->Delete();
1873 id = fParFilter->Last()->GetUniqueID();
1874 if (id > 1) {
1875 trackK->fParFilter->Last()->SetUniqueID(id-1);
1876 trackK->fParFilter->RemoveAt(i);
1877 }
1878 else ((TMatrixD*)(trackK->fParFilter->RemoveAt(i)))->Delete();
1879 id = trackK->fCovExtrap->Last()->GetUniqueID();
1880 if (id > 1) {
1881 trackK->fCovExtrap->Last()->SetUniqueID(id-1);
1882 trackK->fCovExtrap->RemoveAt(i);
1883 }
1884 else ((TMatrixD*)(trackK->fCovExtrap->RemoveAt(i)))->Delete();
1885 id = trackK->fCovFilter->Last()->GetUniqueID();
1886 if (id > 1) {
1887 trackK->fCovFilter->Last()->SetUniqueID(id-1);
1888 trackK->fCovFilter->RemoveAt(i);
1889 }
1890 else ((TMatrixD*)(trackK->fCovFilter->RemoveAt(i)))->Delete();
1891 id = trackK->fJacob->Last()->GetUniqueID();
1892 if (id > 1) {
1893 trackK->fJacob->Last()->SetUniqueID(id-1);
1894 trackK->fJacob->RemoveAt(i);
1895 }
1896 else ((TMatrixD*)(trackK->fJacob->RemoveAt(i)))->Delete();
1897}
1898
1899 //__________________________________________________________________________
b035a148 1900Bool_t AliMUONTrackK::Smooth(void)
1901{
d19b6003 1902/// Apply smoother
1ffe5558 1903 Int_t ihit = fNmbTrackHits - 1;
1904 AliMUONHitForRec *hit = (AliMUONHitForRec*) (*fTrackHits)[ihit];
1905 fChi2Smooth = new TArrayD(fNmbTrackHits);
b035a148 1906 fChi2Smooth->Reset(-1);
1907 fChi2 = 0;
1908 fParSmooth = new TObjArray(15);
1909 fParSmooth->Clear();
1910
1911 if (fgDebug > 0) {
1912 cout << " ******** Enter Smooth " << endl;
1913 cout << (*fSteps)[fNSteps-1] << " " << fPosition << " " << fPositionNew << endl;
1914 /*
1915 for (Int_t i=fNSteps-1; i>=0; i--) cout << (*fSteps)[i] << " ";
1916 cout << endl;
b21aa6eb 1917 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 1918 */
1ffe5558 1919 for (Int_t i=ihit; i>=0; i--) cout << ((AliMUONHitForRec*)(*fTrackHits)[i])->GetZ() << " ";
b035a148 1920 cout << endl;
1921 }
1922
1923 // Find last point corresponding to the last hit
1924 Int_t iLast = fNSteps - 1;
1925 for ( ; iLast>=0; iLast--) {
1ffe5558 1926 //AZ(z->-z) if ((*fSteps)[iLast] + 0.01 > ((AliMUONHitForRec*)(*fTrackHits)[fNmbTrackHits-1])->GetZ()) break;
1927 if (TMath::Abs((*fSteps)[iLast]) + 0.01 > TMath::Abs(((AliMUONHitForRec*)(*fTrackHits)[fNmbTrackHits-1])->GetZ())) break;
b035a148 1928 }
1929
1930 TMatrixD parSmooth = *((TMatrixD*)(fParFilter->UncheckedAt(iLast))); // last filtered = first smoothed
1931 //parSmooth.Dump();
1932 TMatrixD covSmooth = *((TMatrixD*)(fCovFilter->UncheckedAt(iLast))); // last filtered = first smoothed
1933 TMatrixD tmp=covSmooth, weight=covSmooth, gain = covSmooth;
1934 TMatrixD tmpPar = *fTrackPar;
1935 //parSmooth.Print(); ((TMatrixD*)(fParExtrap->Last()))->Print();
1936
1937 Bool_t found;
f29ba3e1 1938 Double_t chi2max = 0;
b035a148 1939 for (Int_t i=iLast+1; i>0; i--) {
1940 if (i == iLast + 1) goto L33; // temporary fix
1941
1942 // Smoother gain matrix
1943 weight = *((TMatrixD*)(fCovExtrap->UncheckedAt(i)));
1944 if (weight.Determinant() != 0) {
b21aa6eb 1945 Int_t ifail;
1946 mnvertLocal(&(weight(0,0)), fgkSize,fgkSize,fgkSize,ifail);
b035a148 1947 } else {
b21aa6eb 1948 AliWarning(" Determinant weight=0:");
b035a148 1949 }
1950 // Fk'Wkk+1
1951 tmp = TMatrixD(*((TMatrixD*)(fJacob->UncheckedAt(i))),TMatrixD::kTransposeMult,weight);
1952 gain = TMatrixD(*((TMatrixD*)(fCovFilter->UncheckedAt(i-1))),TMatrixD::kMult,tmp);
1953 //cout << (*fSteps)[i-1] << " " << (*fSteps)[i] << endl; gain.Print();
1954
1955 // Smoothed parameter vector
1956 //((TMatrixD*)(fParExtrap->UncheckedAt(i)))->Dump();
1957 parSmooth -= *((TMatrixD*)(fParExtrap->UncheckedAt(i)));
1958 tmpPar = *((TMatrixD*)(fParFilter->UncheckedAt(i-1)));
1959 tmpPar += TMatrixD(gain,TMatrixD::kMult,parSmooth);
1960 parSmooth = tmpPar;
1961
1962 // Smoothed covariance
1963 covSmooth -= *((TMatrixD*)(fCovExtrap->UncheckedAt(i)));
1964 weight = TMatrixD(TMatrixD::kTransposed,gain);
1965 tmp = TMatrixD(covSmooth,TMatrixD::kMult,weight);
1966 covSmooth = *((TMatrixD*)(fCovFilter->UncheckedAt(i-1)));
1967 covSmooth += TMatrixD(gain,TMatrixD::kMult,tmp);
1968
1969 // Check if there was a measurement at given z
1970 found = kFALSE;
1971 for ( ; ihit>=0; ihit--) {
1ffe5558 1972 hit = (AliMUONHitForRec*) (*fTrackHits)[ihit];
b035a148 1973 if (TMath::Abs(hit->GetZ()-(*fSteps)[i-1]) < 0.1) { found = kTRUE; break; }
1ffe5558 1974 //AZ(z->-z) else if (ihit < fNmbTrackHits-1 && hit->GetZ() > (*fSteps)[i-1]) { ihit++; break; }
1975 else if (ihit < fNmbTrackHits-1 && TMath::Abs(hit->GetZ()) > TMath::Abs((*fSteps)[i-1])) { ihit++; break; }
b035a148 1976 }
1977 if (!found) continue; // no measurement - skip the rest
1978 else if (fgDebug > 1) cout << i << " " << ihit << " " << hit->GetZ() << endl;
1979 if (ihit == 0) continue; // the first hit - skip the rest
1980
1981L33:
1982 // Smoothed residuals
1983 tmpPar = 0;
1984 tmpPar(0,0) = hit->GetBendingCoor() - parSmooth(0,0); // bending coordinate
1985 tmpPar(1,0) = hit->GetNonBendingCoor() - parSmooth(1,0); // non-bending coordinate
1986 if (fgDebug > 1) {
1987 cout << hit->GetBendingCoor() << " " << parSmooth(0,0) << " " << tmpPar(0,0) << endl;
1988 cout << hit->GetNonBendingCoor() << " " << parSmooth(1,0) << " " << tmpPar(1,0) << endl;
1989 }
1990 // Cov. matrix of smoothed residuals
1991 tmp = 0;
1992 tmp(0,0) = hit->GetBendingReso2() - covSmooth(0,0);
1993 tmp(1,1) = hit->GetNonBendingReso2() - covSmooth(1,1);
1994 tmp(0,1) = tmp(1,0) = -covSmooth(0,1);
1995 tmp(2,2) = tmp(3,3) = tmp(4,4) = 1;
1996
1997 // Calculate Chi2 of smoothed residuals
1998 if (tmp.Determinant() != 0) {
b21aa6eb 1999 Int_t ifail;
2000 mnvertLocal(&(tmp(0,0)), fgkSize,fgkSize,fgkSize,ifail);
b035a148 2001 } else {
b21aa6eb 2002 AliWarning(" Determinant tmp=0:");
b035a148 2003 }
2004 TMatrixD vector(tmp,TMatrixD::kMult,tmpPar);
2005 TMatrixD chi2(tmpPar,TMatrixD::kTransposeMult,vector);
2006 if (fgDebug > 1) chi2.Print();
2007 (*fChi2Smooth)[ihit] = chi2(0,0);
f29ba3e1 2008 if (chi2(0,0) > chi2max) chi2max = chi2(0,0);
b035a148 2009 fChi2 += chi2(0,0);
2010 if (chi2(0,0) < 0) {
94eb555e 2011 //chi2.Print();
2012 AliError(Form(" *** chi2 < 0: %d %d ", i, iLast));
b035a148 2013 }
2014 // Save smoothed parameters
2015 TMatrixD *par = new TMatrixD(parSmooth);
b21aa6eb 2016 fParSmooth->AddAtAndExpand(par, ihit);
b035a148 2017
2018 } // for (Int_t i=iLast+1;
2019
f29ba3e1 2020 //if (chi2max > 16) {
2021 //if (chi2max > 25) {
2022 //if (chi2max > 50) {
2023 //if (chi2max > 100) {
2024 if (chi2max > fgkChi2max) {
b035a148 2025 //if (Recover()) DropBranches();
2026 //Recover();
2027 Outlier();
2028 return kFALSE;
2029 }
2030 return kTRUE;
2031}
2032
2033 //__________________________________________________________________________
2034void AliMUONTrackK::Outlier()
2035{
d19b6003 2036/// Adds new track with removed hit having the highest chi2
b035a148 2037
2038 if (fgDebug > 0) {
2039 cout << " ******** Enter Outlier " << endl;
1ffe5558 2040 for (Int_t i=0; i<fNmbTrackHits; i++) printf("%10.4f", (*fChi2Smooth)[i]);
b035a148 2041 printf("\n");
1ffe5558 2042 for (Int_t i=0; i<fNmbTrackHits; i++) {
2043 printf("%10d", ((AliMUONHitForRec*)fTrackHits->UncheckedAt(i))->GetChamberNumber());
b035a148 2044 }
2045 printf("\n");
2046 }
2047
f29ba3e1 2048 Double_t chi2max = 0;
b035a148 2049 Int_t imax = 0;
1ffe5558 2050 for (Int_t i=0; i<fNmbTrackHits; i++) {
f29ba3e1 2051 if ((*fChi2Smooth)[i] < chi2max) continue;
2052 chi2max = (*fChi2Smooth)[i];
b035a148 2053 imax = i;
2054 }
2055 // Check if the outlier is not from the seed segment
1ffe5558 2056 AliMUONHitForRec *hit = (AliMUONHitForRec*) fTrackHits->UncheckedAt(imax);
b035a148 2057 if (hit == fStartSegment->GetHitForRec1() || hit == fStartSegment->GetHitForRec2()) return; // to be changed probably
2058
2059 // Check for missing station
2060 Int_t ok = 1;
2061 if (imax == 0) {
1ffe5558 2062 if (((AliMUONHitForRec*)fTrackHits->UncheckedAt(1))->GetChamberNumber() < 8) ok--;
2063 } else if (imax == fNmbTrackHits-1) {
2064 if (((AliMUONHitForRec*)fTrackHits->UncheckedAt(fNmbTrackHits-2))->GetChamberNumber() > 1) ok--;
b035a148 2065 }
1ffe5558 2066 else if (((AliMUONHitForRec*)fTrackHits->UncheckedAt(imax-1))->GetChamberNumber()/2 - ((AliMUONHitForRec*)fTrackHits->UncheckedAt(imax+1))->GetChamberNumber()/2 > 1) ok--;
b035a148 2067 if (!ok) { Recover(); return; } // try to recover track
2068 //AZ if (!ok) { if (fgDebug >= 0) cout << imax << endl; DropBranches(imax, 0); return; }
2069
2070 // Remove saved steps and smoother matrices after the outlier
2071 RemoveMatrices(hit->GetZ());
2072
2073 // Check for possible double track candidates
2074 //if (ExistDouble(hit)) return;
2075
29f1b13a 2076 TClonesArray *trackPtr = fgTrackReconstructor->GetRecTracksPtr();
2077 Int_t nRecTracks = fgTrackReconstructor->GetNRecTracks();
b035a148 2078
2079 AliMUONTrackK *trackK = 0;
2080 if (!fNSteps || GetStation0() == 3 && hit->GetChamberNumber() > 7) {
2081 // start track from segment
2082 trackK = new ((*trackPtr)[nRecTracks]) AliMUONTrackK(fStartSegment);
29f1b13a 2083 fgTrackReconstructor->SetNRecTracks(nRecTracks+1);
b035a148 2084 trackK->fRecover = 2;
2085 trackK->fSkipHit = hit;
1ffe5558 2086 trackK->fNmbTrackHits = fNmbTrackHits;
2087
2088 hit = (AliMUONHitForRec*) trackK->fTrackHits->UncheckedAt(0);
2089 hit->SetNTrackHits(hit->GetNTrackHits()-1);
2090 hit = (AliMUONHitForRec*) trackK->fTrackHits->UncheckedAt(1);
2091 hit->SetNTrackHits(hit->GetNTrackHits()-1);
2092 delete trackK->fTrackHits; // not efficient ?
2093 trackK->fTrackHits = new TObjArray(*fTrackHits);
2094 for (Int_t i = 0; i < fNmbTrackHits; i++) {
2095 hit = (AliMUONHitForRec*) fTrackHits->UncheckedAt(i);
2096 hit->SetNTrackHits(hit->GetNTrackHits()+1);
2097 }
2098
2099 if (GetStation0() == 3) trackK->SortHits(1, trackK->fTrackHits);
b035a148 2100 if (fgDebug > 0) cout << nRecTracks << " Outlier" << endl;
2101 return;
2102 }
2103 trackK = new ((*trackPtr)[nRecTracks]) AliMUONTrackK(NULL, NULL);
2104 *trackK = *this;
29f1b13a 2105 fgTrackReconstructor->SetNRecTracks(nRecTracks+1);
b035a148 2106 trackK->fTrackDir = 1;
2107 trackK->fRecover = 2;
2108 trackK->fSkipHit = hit;
2109 Int_t iD = trackK->fParFilter->Last()->GetUniqueID();
2110 if (iD > 1) {
2111 trackK->fParFilter->Last()->SetUniqueID(iD-1);
2112 CreateMatrix(trackK->fParFilter);
b035a148 2113 }
2114 *((TMatrixD*)trackK->fParFilter->Last()) = *((TMatrixD*)fParExtrap->Last());
0ce28091 2115 trackK->fParFilter->Last()->SetUniqueID(1);
b035a148 2116 *(trackK->fTrackPar) = *((TMatrixD*)trackK->fParFilter->Last());
2117 iD = trackK->fCovFilter->Last()->GetUniqueID();
2118 if (iD > 1) {
2119 trackK->fCovFilter->Last()->SetUniqueID(iD-1);
2120 CreateMatrix(trackK->fCovFilter);
b035a148 2121 }
2122 *((TMatrixD*)trackK->fCovFilter->Last()) = *((TMatrixD*)fCovExtrap->Last());
0ce28091 2123 trackK->fCovFilter->Last()->SetUniqueID(1);
b035a148 2124 *(trackK->fWeight) = *((TMatrixD*)trackK->fCovFilter->Last());
2125 if (trackK->fWeight->Determinant() != 0) {
b21aa6eb 2126 Int_t ifail;
2127 mnvertLocal(&((*(trackK->fWeight))(0,0)), fgkSize,fgkSize,fgkSize,ifail);
b035a148 2128 } else {
b21aa6eb 2129 AliWarning(" Determinant fWeight=0:");
b035a148 2130 }
2131 trackK->fPosition = trackK->fPositionNew = (*fSteps)[fNSteps-1];
2132 trackK->fChi2 = 0;
1ffe5558 2133 for (Int_t i=0; i<fNmbTrackHits-1; i++) trackK->fChi2 += (*fChi2Array)[i];
b035a148 2134 if (fgDebug > 0) cout << nRecTracks << " Outlier " << trackK->fChi2 << endl;
2135}
2136
2137 //__________________________________________________________________________
2138Double_t AliMUONTrackK::GetChi2PerPoint(Int_t iPoint) const
2139{
d19b6003 2140/// Return Chi2 at point
b035a148 2141 return fChi2Smooth ? (*fChi2Smooth)[iPoint] : (*fChi2Array)[iPoint];
2142 //return 0.;
2143}
2144
2145 //__________________________________________________________________________
2146void AliMUONTrackK::Print(FILE *lun) const
2147{
d19b6003 2148/// Print out track information
b035a148 2149
2150 Int_t flag = 1;
2151 AliMUONHitForRec *hit = 0;
b035a148 2152 // from raw clusters
8f373020 2153 AliMUONRawCluster *clus = 0;
2154 TClonesArray *rawclusters = 0;
2155 for (Int_t i1=0; i1<fNmbTrackHits; i1++) {
2156 hit = (AliMUONHitForRec*) ((*fTrackHits)[i1]);
2157 rawclusters = fgTrackReconstructor->GetMUONData()->RawClusters(hit->GetChamberNumber());
2158 clus = (AliMUONRawCluster*) rawclusters->UncheckedAt(hit->GetHitNumber());
2159 if (TMath::Abs(clus->GetTrack(1)-1) < 2) {
2160 if (clus->GetTrack(2)) flag = 2;
2161 continue;
2162 }
2163 if (clus->GetTrack(2) && TMath::Abs(clus->GetTrack(2)-1) < 2) {
2164 flag = 3;
2165 continue;
b035a148 2166 }
8f373020 2167 flag = 0;
2168 break;
2169
b035a148 2170 Int_t sig[2]={1,1}, tid[2]={0};
1ffe5558 2171 for (Int_t i1=0; i1<fNmbTrackHits; i1++) {
b035a148 2172 if (GetChi2PerPoint(i1) < -0.1) continue;
1ffe5558 2173 hit = (AliMUONHitForRec*) ((*fTrackHits)[i1]);
29f1b13a 2174 rawclusters = fgTrackReconstructor->GetMUONData()->RawClusters(hit->GetChamberNumber());
b035a148 2175 clus = (AliMUONRawCluster*) rawclusters->UncheckedAt(hit->GetHitNumber());
2176 for (Int_t j=0; j<2; j++) {
2177 tid[j] = clus->GetTrack(j+1) - 1;
2178 if (clus->GetTrack(j+1) < 0) { sig[j] = 0; tid[j] = 999; }
2179 }
2180 fprintf(lun,"%3d %3d %10.4f", gAlice->GetEvNumber(), hit->GetChamberNumber(), GetChi2PerPoint(i1));
2181 if (!(clus->GetTrack(2))) fprintf(lun, "%3d %3d", sig[0], tid[0]); // simple cluster
2182 else { // track overlap
2183 fprintf(lun, "%3d %3d", TMath::Max(sig[0],sig[1]), TMath::Min(tid[0],tid[1]));
2184 //if (tid[0] < 2) flag *= 2;
2185 //else if (tid[1] < 2) flag *= 3;
2186 }
2187 fprintf (lun, "%3d \n", flag);
2188 }
2189 }
2190}
2191
2192 //__________________________________________________________________________
2193void AliMUONTrackK::DropBranches(Int_t imax, TObjArray *hits)
2194{
d19b6003 2195/// Drop branches downstream of the skipped hit
b035a148 2196 Int_t nRecTracks;
2197 TClonesArray *trackPtr;
2198 AliMUONTrackK *trackK;
2199
29f1b13a 2200 trackPtr = fgTrackReconstructor->GetRecTracksPtr();
2201 nRecTracks = fgTrackReconstructor->GetNRecTracks();
b035a148 2202 Int_t icand = trackPtr->IndexOf(this);
1ffe5558 2203 if (!hits) hits = fTrackHits;
b035a148 2204
2205 // Check if the track candidate doesn't exist yet
2206 for (Int_t i=icand+1; i<nRecTracks; i++) {
2207 trackK = (AliMUONTrackK*) ((*trackPtr)[i]);
1ffe5558 2208 if (trackK->fNmbTrackHits == 2 && trackK->GetRecover() == 0) continue;
b035a148 2209 if (trackK->GetRecover() < 0) continue;
2210
1ffe5558 2211 if (trackK->fNmbTrackHits >= imax + 1) {
b035a148 2212 for (Int_t j=0; j<=imax; j++) {
1ffe5558 2213 //if (j != fNmbTrackHits-1 && (*trackK->fTrackHits)[j] != (*fTrackHits)[j]) break;
2214 if ((*trackK->fTrackHits)[j] != (*hits)[j]) break;
b035a148 2215 if (j == imax) {
1ffe5558 2216 if (hits != fTrackHits) {
b035a148 2217 // Drop all branches downstream the hit (for Recover)
2218 trackK->SetRecover(-1);
2219 if (fgDebug >= 0 )cout << " Recover " << i << endl;
2220 continue;
2221 }
2222 // Check if the removal of the hit breaks the track
2223 Int_t ok = 1;
2224 if (imax == 0) {
1ffe5558 2225 if (((AliMUONHitForRec*)trackK->fTrackHits->UncheckedAt(1))->GetChamberNumber() < 8) ok--; }
2226 else if (imax == trackK->fNmbTrackHits-1) continue;
2227 // else if (imax == trackK->fNmbTrackHits-1) {
2228 //if (((AliMUONHitForRec*)trackK->fTrackHits->UncheckedAt(trackK->fNmbTrackHits-2))->GetChamberNumber() > 1) ok--;
b035a148 2229 //}
1ffe5558 2230 else if (((AliMUONHitForRec*)trackK->fTrackHits->UncheckedAt(imax-1))->GetChamberNumber()/2 - ((AliMUONHitForRec*)trackK->fTrackHits->UncheckedAt(imax+1))->GetChamberNumber()/2 > 1) ok--;
b035a148 2231 if (!ok) trackK->SetRecover(-1);
2232 }
2233 } // for (Int_t j=0;
2234 }
2235 } // for (Int_t i=0;
2236}
2237
2238 //__________________________________________________________________________
2239void AliMUONTrackK::DropBranches(AliMUONSegment *segment)
2240{
d19b6003 2241/// Drop all candidates with the same seed segment
b035a148 2242 Int_t nRecTracks;
2243 TClonesArray *trackPtr;
2244 AliMUONTrackK *trackK;
2245
29f1b13a 2246 trackPtr = fgTrackReconstructor->GetRecTracksPtr();
2247 nRecTracks = fgTrackReconstructor->GetNRecTracks();
b035a148 2248 Int_t icand = trackPtr->IndexOf(this);
2249
2250 for (Int_t i=icand+1; i<nRecTracks; i++) {
2251 trackK = (AliMUONTrackK*) ((*trackPtr)[i]);
1ffe5558 2252 if (trackK->fNmbTrackHits == 2 && trackK->GetRecover() == 0) continue;
b035a148 2253 if (trackK->GetRecover() < 0) continue;
2254 if (trackK->fStartSegment == segment) trackK->SetRecover(-1);
2255 }
2256 if (fgDebug >= 0) cout << " Drop segment " << endl;
2257}
2258
2259 //__________________________________________________________________________
2260AliMUONHitForRec* AliMUONTrackK::GetHitLastOk(void)
2261{
d19b6003 2262/// Return the hit where track stopped
b035a148 2263
1ffe5558 2264 if (!fNSteps) return (AliMUONHitForRec*)((*fTrackHits)[1]);
b035a148 2265 return fSkipHit;
2266}
2267
2268 //__________________________________________________________________________
2269Int_t AliMUONTrackK::GetStation0(void)
2270{
d19b6003 2271/// Return seed station number
b035a148 2272 return fStartSegment->GetHitForRec1()->GetChamberNumber() / 2;
2273}
2274
2275 //__________________________________________________________________________
2276Bool_t AliMUONTrackK::ExistDouble(AliMUONHitForRec *hit)
2277{
d19b6003 2278/// Check if the track will make a double after outlier removal
b035a148 2279
29f1b13a 2280 TClonesArray *trackPtr = fgTrackReconstructor->GetRecTracksPtr();
2281 Int_t nRecTracks = fgTrackReconstructor->GetNRecTracks();
1ffe5558 2282 TObjArray *hitArray = new TObjArray(*fTrackHits);
b035a148 2283 TObjArray *hitArray1 = new TObjArray(*hitArray);
2284 hitArray1->Remove(hit);
2285 hitArray1->Compress();
2286
2287 Bool_t same = kFALSE;
2288 for (Int_t i=0; i<nRecTracks; i++) {
2289 AliMUONTrackK *trackK = (AliMUONTrackK*) ((*trackPtr)[i]);
1ffe5558 2290 if (trackK->fNmbTrackHits == 2 && trackK->GetRecover() == 0) continue;
b035a148 2291 if (trackK == this) continue;
1ffe5558 2292 if (trackK->fNmbTrackHits == fNmbTrackHits || trackK->fNmbTrackHits == fNmbTrackHits-1) {
2293 TObjArray *hits = new TObjArray(*trackK->fTrackHits);
b035a148 2294 same = kTRUE;
1ffe5558 2295 if (trackK->fNmbTrackHits == fNmbTrackHits) {
2296 for (Int_t j=0; j<fNmbTrackHits; j++) {
b035a148 2297 if (hits->UncheckedAt(j) != hitArray->UncheckedAt(j)) { same = kFALSE; break; }
2298 }
2299 if (same) { delete hits; break; }
2300 if (trackK->fSkipHit) {
2301 TObjArray *hits1 = new TObjArray(*hits);
2302 if (hits1->Remove(trackK->fSkipHit) > 0) {
2303 hits1->Compress();
2304 same = kTRUE;
1ffe5558 2305 for (Int_t j=0; j<fNmbTrackHits-1; j++) {
b035a148 2306 if (hits1->UncheckedAt(j) != hitArray1->UncheckedAt(j)) { same = kFALSE; break; }
2307 }
2308 if (same) { delete hits1; break; }
2309 }
2310 delete hits1;
2311 }
2312 } else {
2313 // Check with removed outlier
2314 same = kTRUE;
1ffe5558 2315 for (Int_t j=0; j<fNmbTrackHits-1; j++) {
b035a148 2316 if (hits->UncheckedAt(j) != hitArray1->UncheckedAt(j)) { same = kFALSE; break; }
2317 }
2318 if (same) { delete hits; break; }
2319 }
2320 delete hits;
2321 }
2322 } // for (Int_t i=0; i<nRecTracks;
2323 delete hitArray; delete hitArray1;
2324 if (same && fgDebug >= 0) cout << " Same" << endl;
2325 return same;
2326}
2327
2328 //__________________________________________________________________________
2329Bool_t AliMUONTrackK::ExistDouble(void)
2330{
d19b6003 2331/// Check if the track will make a double after recovery
b035a148 2332
29f1b13a 2333 TClonesArray *trackPtr = fgTrackReconstructor->GetRecTracksPtr();
2334 Int_t nRecTracks = fgTrackReconstructor->GetNRecTracks();
b035a148 2335
1ffe5558 2336 TObjArray *hitArray = new TObjArray(*fTrackHits);
b035a148 2337 if (GetStation0() == 3) SortHits(0, hitArray); // sort
2338 //if (GetStation0() == 3) SortHits(1, hitArray); // unsort
2339
2340 Bool_t same = kFALSE;
2341 for (Int_t i=0; i<nRecTracks; i++) {
2342 AliMUONTrackK *trackK = (AliMUONTrackK*) ((*trackPtr)[i]);
1ffe5558 2343 if (trackK->fNmbTrackHits == 2 && trackK->GetRecover() == 0) continue;
b035a148 2344 if (trackK == this) continue;
2345 //AZ if (trackK->GetRecover() < 0) continue; //
1ffe5558 2346 if (trackK->fNmbTrackHits >= fNmbTrackHits) {
2347 TObjArray *hits = new TObjArray(*trackK->fTrackHits);
b035a148 2348 if (trackK->GetStation0() == 3) SortHits(0, hits); // sort
2349 //if (trackK->GetStation0() == 3) SortHits(1, hits); // unsort
1ffe5558 2350 for (Int_t j=0; j<fNmbTrackHits; j++) {
2351 //cout << fNmbTrackHits << " " << i << " " << j << " " << (*hitArray)[j] << " " << (*hits)[j] << " " << trackK->fSkipHit << endl;
2352 if (j != fNmbTrackHits-1 && (*hitArray)[j] != (*hits)[j]) break;
2353 if (j == fNmbTrackHits-1) {
b035a148 2354 if (trackK->fSkipHit && TMath::Abs(((AliMUONHitForRec*)((*hitArray)[j]))->GetZ()-trackK->fSkipHit->GetZ()) < 0.5) same = kTRUE;
1ffe5558 2355 //if (trackK->fNmbTrackHits > fNmbTrackHits &&
b035a148 2356 //if (trackK->fSkipHit) cout << ((AliMUONHitForRec*)((*hitArray)[j]))->GetZ() << " " << trackK->fSkipHit->GetZ() << endl;
2357 }
2358 } // for (Int_t j=0;
2359 delete hits;
2360 if (same) break;
2361 }
2362 } // for (Int_t i=0;
2363 delete hitArray;
2364 return same;
2365}