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