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