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