1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
18 // ---------------------
19 // Class AliMUONTrackK
20 // ---------------------
21 // Reconstructed track in the muons system based on the extended
22 // Kalman filter approach
23 // Author: Alexander Zinchenko, JINR Dubna
25 #include <Riostream.h>
26 #include <TClonesArray.h>
30 #include "AliMUONTrackK.h"
33 #include "AliMUONTrackReconstructor.h"
35 #include "AliMUONSegment.h"
36 #include "AliMUONHitForRec.h"
37 #include "AliMUONRawCluster.h"
38 #include "AliMUONTrackParam.h"
39 #include "AliMUONEventRecoCombi.h"
40 #include "AliMUONDetElement.h"
44 const Int_t AliMUONTrackK::fgkSize = 5;
45 const Int_t AliMUONTrackK::fgkNSigma = 12;
46 const Double_t AliMUONTrackK::fgkChi2max = 100;
47 const Int_t AliMUONTrackK::fgkTriesMax = 10000;
48 const Double_t AliMUONTrackK::fgkEpsilon = 0.002;
50 void mnvertLocal(Double_t* a, Int_t l, Int_t m, Int_t n, Int_t& ifail); // from AliMUONTrack
52 Int_t AliMUONTrackK::fgDebug = -1; //-1;
53 Int_t AliMUONTrackK::fgNOfPoints = 0;
54 AliMUONTrackReconstructor* AliMUONTrackK::fgTrackReconstructor = NULL;
55 TClonesArray* AliMUONTrackK::fgHitForRec = NULL;
56 AliMUONEventRecoCombi *AliMUONTrackK::fgCombi = NULL;
58 ClassImp(AliMUONTrackK) // Class implementation in ROOT context
60 //__________________________________________________________________________
61 AliMUONTrackK::AliMUONTrackK()
88 /// Default constructor
90 fgTrackReconstructor = NULL; // pointer to event reconstructor
91 fgHitForRec = NULL; // pointer to points
92 fgNOfPoints = 0; // number of points
97 //__________________________________________________________________________
98 AliMUONTrackK::AliMUONTrackK(AliMUONTrackReconstructor *TrackReconstructor, TClonesArray *hitForRec)
127 if (!TrackReconstructor) return;
128 fgTrackReconstructor = TrackReconstructor; // pointer to event reconstructor
129 fgHitForRec = hitForRec; // pointer to points
130 fgNOfPoints = fgHitForRec->GetEntriesFast(); // number of points
132 if (fgTrackReconstructor->GetTrackMethod() == 3) fgCombi = AliMUONEventRecoCombi::Instance();
135 //__________________________________________________________________________
136 AliMUONTrackK::AliMUONTrackK(AliMUONSegment *segment)
137 //: AliMUONTrack(segment, segment, fgTrackReconstructor)
138 : AliMUONTrack(NULL, segment, fgTrackReconstructor),
139 fStartSegment(segment),
143 fTrackHits(new TObjArray(13)),
149 fTrackPar(new TMatrixD(fgkSize,1)),
150 fTrackParNew(new TMatrixD(fgkSize,1)),
151 fCovariance(new TMatrixD(fgkSize,fgkSize)),
152 fWeight(new TMatrixD(fgkSize,fgkSize)),
153 fParExtrap(new TObjArray(15)),
154 fParFilter(new TObjArray(15)),
156 fCovExtrap(new TObjArray(15)),
157 fCovFilter(new TObjArray(15)),
158 fJacob(new TObjArray(15)),
160 fSteps(new TArrayD(15)),
161 fChi2Array(new TArrayD(13)),
164 /// Constructor from a segment
166 AliMUONHitForRec *hit1, *hit2;
167 AliMUONRawCluster *clus;
168 TClonesArray *rawclusters;
170 // Pointers to hits from the segment
171 hit1 = segment->GetHitForRec1();
172 hit2 = segment->GetHitForRec2();
173 hit1->SetNTrackHits(hit1->GetNTrackHits()+1); // mark hit as being on track
174 hit2->SetNTrackHits(hit2->GetNTrackHits()+1); // mark hit as being on track
175 // check sorting in Z
176 if (TMath::Abs(hit1->GetZ()) > TMath::Abs(hit2->GetZ())) {
178 hit2 = segment->GetHitForRec1();
181 // Fill array of track parameters
182 if (hit1->GetChamberNumber() > 7) {
183 // last tracking station
184 (*fTrackPar)(0,0) = hit1->GetBendingCoor(); // y
185 (*fTrackPar)(1,0) = hit1->GetNonBendingCoor(); // x
186 fPosition = hit1->GetZ(); // z
187 fTrackHits->Add(hit2); // add hit 2
188 fTrackHits->Add(hit1); // add hit 1
189 //AZ(Z->-Z) fTrackDir = -1;
192 // last but one tracking station
193 (*fTrackPar)(0,0) = hit2->GetBendingCoor(); // y
194 (*fTrackPar)(1,0) = hit2->GetNonBendingCoor(); // x
195 fPosition = hit2->GetZ(); // z
196 fTrackHits->Add(hit1); // add hit 1
197 fTrackHits->Add(hit2); // add hit 2
198 //AZ(Z->-Z) fTrackDir = 1;
201 dZ = hit2->GetZ() - hit1->GetZ();
202 dY = hit2->GetBendingCoor() - hit1->GetBendingCoor();
203 dX = hit2->GetNonBendingCoor() - hit1->GetNonBendingCoor();
204 (*fTrackPar)(2,0) = TMath::ATan2(dY,dZ); // alpha
205 if ((*fTrackPar)(2,0) < 0.) (*fTrackPar)(2,0) += 2*TMath::Pi(); // from 0 to 2*pi
206 (*fTrackPar)(3,0) = TMath::ATan2(-dX,dZ/TMath::Cos((*fTrackPar)(2,0))); // beta
207 (*fTrackPar)(2,0) -= TMath::Pi();
208 (*fTrackPar)(4,0) = 1/fgTrackReconstructor->GetBendingMomentumFromImpactParam(segment->GetBendingImpact()); // 1/Pt
209 (*fTrackPar)(4,0) *= TMath::Cos((*fTrackPar)(3,0)); // 1/p
210 // Evaluate covariance (and weight) matrix
213 if (fgDebug < 0 ) return;
214 cout << fgTrackReconstructor->GetNRecTracks()-1 << " " << fgTrackReconstructor->GetBendingMomentumFromImpactParam(segment->GetBendingImpact()) << " " << 1/(*fTrackPar)(4,0) << " ";
216 for (Int_t i=0; i<2; i++) {
217 hit1 = (AliMUONHitForRec*) ((*fTrackHits)[i]);
218 rawclusters = fgTrackReconstructor->GetMUONData()->RawClusters(hit1->GetChamberNumber());
219 clus = (AliMUONRawCluster*) rawclusters->UncheckedAt(hit1->GetHitNumber());
220 cout << clus->GetTrack(1);
221 if (clus->GetTrack(2) != -1) cout << " " << clus->GetTrack(2);
222 if (i == 0) cout << " <--> ";
224 cout << " @ " << fStartSegment->GetHitForRec1()->GetChamberNumber() << endl;
228 //__________________________________________________________________________
229 AliMUONTrackK::~AliMUONTrackK()
234 //cout << fNmbTrackHits << endl;
235 for (Int_t i = 0; i < fNmbTrackHits; i++) {
236 AliMUONHitForRec *hit = (AliMUONHitForRec*) fTrackHits->UncheckedAt(i);
237 hit->SetNTrackHits(hit->GetNTrackHits()-1);
239 delete fTrackHits; // delete the TObjArray of pointers to TrackHit's
243 delete fTrackPar; delete fTrackParNew; delete fCovariance;
247 if (fSteps) delete fSteps;
248 if (fChi2Array) delete fChi2Array;
249 if (fChi2Smooth) delete fChi2Smooth;
250 if (fParSmooth) {fParSmooth->Delete(); delete fParSmooth; }
251 // Delete only matrices not shared by several tracks
252 if (!fParExtrap) return;
255 for (Int_t i=0; i<fNSteps; i++) {
256 //if (fParExtrap->UncheckedAt(i)->GetUniqueID() > 1)
257 // fParExtrap->UncheckedAt(i)->SetUniqueID(fParExtrap->RemoveAt(i)->GetUniqueID()-1);
258 id = fParExtrap->UncheckedAt(i)->GetUniqueID();
260 fParExtrap->UncheckedAt(i)->SetUniqueID(id-1);
261 fParExtrap->RemoveAt(i);
263 //if (fParFilter->UncheckedAt(i)->GetUniqueID() > 1)
264 // fParFilter->UncheckedAt(i)->SetUniqueID(fParFilter->RemoveAt(i)->GetUniqueID()-1);
265 id = fParFilter->UncheckedAt(i)->GetUniqueID();
267 fParFilter->UncheckedAt(i)->SetUniqueID(id-1);
268 fParFilter->RemoveAt(i);
270 //if (fCovExtrap->UncheckedAt(i)->GetUniqueID() > 1)
271 // fCovExtrap->UncheckedAt(i)->SetUniqueID(fCovExtrap->RemoveAt(i)->GetUniqueID()-1);
272 id = fCovExtrap->UncheckedAt(i)->GetUniqueID();
274 fCovExtrap->UncheckedAt(i)->SetUniqueID(id-1);
275 fCovExtrap->RemoveAt(i);
278 //if (fCovFilter->UncheckedAt(i)->GetUniqueID() > 1)
279 // fCovFilter->UncheckedAt(i)->SetUniqueID(fCovFilter->RemoveAt(i)->GetUniqueID()-1);
280 id = fCovFilter->UncheckedAt(i)->GetUniqueID();
282 fCovFilter->UncheckedAt(i)->SetUniqueID(id-1);
283 fCovFilter->RemoveAt(i);
285 //if (fJacob->UncheckedAt(i)->GetUniqueID() > 1)
286 // fJacob->UncheckedAt(i)->SetUniqueID(fJacob->RemoveAt(i)->GetUniqueID()-1);
287 id = fJacob->UncheckedAt(i)->GetUniqueID();
289 fJacob->UncheckedAt(i)->SetUniqueID(id-1);
294 for (Int_t i=0; i<fNSteps; i++) {
295 if (fParExtrap->UncheckedAt(i)) ((TMatrixD*)fParExtrap->UncheckedAt(i))->Delete();
296 if (fParFilter->UncheckedAt(i)) ((TMatrixD*)fParFilter->UncheckedAt(i))->Delete();
297 if (fCovExtrap->UncheckedAt(i)) ((TMatrixD*)fCovExtrap->UncheckedAt(i))->Delete();
298 cout << fCovFilter->UncheckedAt(i) << " " << (*fSteps)[i] << endl;
299 if (fCovFilter->UncheckedAt(i)) ((TMatrixD*)fCovFilter->UncheckedAt(i))->Delete();
300 if (fJacob->UncheckedAt(i)) ((TMatrixD*)fJacob->UncheckedAt(i))->Delete();
303 fParExtrap->Delete(); fParFilter->Delete();
304 fCovExtrap->Delete(); fCovFilter->Delete();
306 delete fParExtrap; delete fParFilter;
307 delete fCovExtrap; delete fCovFilter;
311 //__________________________________________________________________________
312 AliMUONTrackK & AliMUONTrackK::operator=(const AliMUONTrackK& source)
314 /// Assignment operator
317 if(&source == this) return *this;
319 // base class assignement
320 //AZ TObject::operator=(source);
321 AliMUONTrack::operator=(source);
323 fStartSegment = source.fStartSegment;
324 fNmbTrackHits = source.fNmbTrackHits;
325 fChi2 = source.fChi2;
326 fPosition = source.fPosition;
327 fPositionNew = source.fPositionNew;
328 fTrackDir = source.fTrackDir;
329 fBPFlag = source.fBPFlag;
330 fRecover = source.fRecover;
331 fSkipHit = source.fSkipHit;
334 fTrackHits = new TObjArray(*source.fTrackHits);
335 //source.fTrackHits->Dump();
336 //fTrackHits->Dump();
337 for (Int_t i = 0; i < fNmbTrackHits; i++) {
338 AliMUONHitForRec *hit = (AliMUONHitForRec*) fTrackHits->UncheckedAt(i);
339 hit->SetNTrackHits(hit->GetNTrackHits()+1);
342 fTrackPar = new TMatrixD(*source.fTrackPar); // track parameters
343 fTrackParNew = new TMatrixD(*source.fTrackParNew); // track parameters
344 fCovariance = new TMatrixD(*source.fCovariance); // covariance matrix
345 fWeight = new TMatrixD(*source.fWeight); // weight matrix (inverse of covariance)
348 fParExtrap = new TObjArray(*source.fParExtrap);
349 fParFilter = new TObjArray(*source.fParFilter);
350 fCovExtrap = new TObjArray(*source.fCovExtrap);
351 fCovFilter = new TObjArray(*source.fCovFilter);
352 fJacob = new TObjArray(*source.fJacob);
353 fSteps = new TArrayD(*source.fSteps);
354 fNSteps = source.fNSteps;
356 if (source.fChi2Array) fChi2Array = new TArrayD(*source.fChi2Array);
360 for (Int_t i=0; i<fParExtrap->GetEntriesFast(); i++) {
361 fParExtrap->UncheckedAt(i)->SetUniqueID(fParExtrap->UncheckedAt(i)->GetUniqueID()+1);
362 fParFilter->UncheckedAt(i)->SetUniqueID(fParFilter->UncheckedAt(i)->GetUniqueID()+1);
363 fCovExtrap->UncheckedAt(i)->SetUniqueID(fCovExtrap->UncheckedAt(i)->GetUniqueID()+1);
364 fCovFilter->UncheckedAt(i)->SetUniqueID(fCovFilter->UncheckedAt(i)->GetUniqueID()+1);
365 fJacob->UncheckedAt(i)->SetUniqueID(fJacob->UncheckedAt(i)->GetUniqueID()+1);
370 //__________________________________________________________________________
371 void AliMUONTrackK::EvalCovariance(Double_t dZ)
373 /// Evaluate covariance (and weight) matrix for track candidate
374 Double_t sigmaB, sigmaNonB, tanA, tanB, dAdY, rad, dBdX, dBdY;
377 sigmaB = fgTrackReconstructor->GetBendingResolution(); // bending resolution
378 sigmaNonB = fgTrackReconstructor->GetNonBendingResolution(); // non-bending resolution
380 (*fWeight)(0,0) = sigmaB*sigmaB; // <yy>
382 (*fWeight)(1,1) = sigmaNonB*sigmaNonB; // <xx>
384 tanA = TMath::Tan((*fTrackPar)(2,0));
385 dAdY = 1/(1+tanA*tanA)/dZ;
386 (*fWeight)(2,2) = dAdY*dAdY*(*fWeight)(0,0)*2; // <aa>
387 (*fWeight)(0,2) = dAdY*(*fWeight)(0,0); // <ya>
388 (*fWeight)(2,0) = (*fWeight)(0,2);
390 rad = dZ/TMath::Cos((*fTrackPar)(2,0));
391 tanB = TMath::Tan((*fTrackPar)(3,0));
392 dBdX = 1/(1+tanB*tanB)/rad;
394 (*fWeight)(3,3) = dBdX*dBdX*(*fWeight)(1,1)*2; // <bb>
395 (*fWeight)(1,3) = dBdX*(*fWeight)(1,1); // <xb>
396 (*fWeight)(3,1) = (*fWeight)(1,3);
398 (*fWeight)(4,4) = ((*fTrackPar)(4,0)*0.5)*((*fTrackPar)(4,0)*0.5); // error 50%
400 // check whether the Invert method returns flag if matrix cannot be inverted,
401 // and do not calculate the Determinant in that case !!!!
402 if (fWeight->Determinant() != 0) {
403 // fWeight->Invert();
405 mnvertLocal(&((*fWeight)(0,0)), fgkSize,fgkSize,fgkSize,ifail);
407 AliWarning(" Determinant fWeight=0:");
412 //__________________________________________________________________________
413 Bool_t AliMUONTrackK::KalmanFilter(Int_t ichamBeg, Int_t ichamEnd, Bool_t Back, Double_t zDipole1, Double_t zDipole2)
415 /// Follows track through detector stations
416 Bool_t miss, success;
417 Int_t ichamb, iFB, iMin, iMax, dChamb, ichambOK;
418 Int_t ihit, firstIndx = -1, lastIndx = -1, currIndx = -1, iz0 = -1, iz = -1;
419 Double_t zEnd, dChi2;
420 AliMUONHitForRec *hitAdd, *firstHit = NULL, *lastHit = NULL, *hit;
421 AliMUONRawCluster *clus;
422 TClonesArray *rawclusters;
423 hit = 0; clus = 0; rawclusters = 0;
425 miss = success = kTRUE;
427 //iFB = (ichamEnd == ichamBeg) ? fTrackDir : TMath::Sign(1,ichamEnd-ichamBeg);
428 iFB = (ichamEnd == ichamBeg) ? -fTrackDir : TMath::Sign(1,ichamEnd-ichamBeg);
429 iMin = TMath::Min(ichamEnd,ichamBeg);
430 iMax = TMath::Max(ichamEnd,ichamBeg);
437 if (((AliMUONHitForRec*)fTrackHits->First())->GetChamberNumber() != ichamb) currIndx = 0;
438 } else if (fRecover) {
439 hit = GetHitLastOk();
440 currIndx = fTrackHits->IndexOf(hit);
441 if (currIndx < 0) hit = fStartSegment->GetHitForRec1(); // for station 3
443 ichamb = hit->GetChamberNumber();
444 if (hit == fSkipHit || fRecover == 2 && currIndx >= 0) {
445 // start from the last point or outlier
446 // remove the skipped hit
447 fTrackHits->Remove(fSkipHit); // remove hit
449 fSkipHit->SetNTrackHits(fSkipHit->GetNTrackHits()-1); // unmark hit
454 ichambOK = ((AliMUONHitForRec*)((*fTrackHits)[fNmbTrackHits-1]))->GetChamberNumber();
455 //if (ichambOK >= 8 && ((AliMUONHitForRec*)((*fTrackHits)[0]))->GetChamberNumber() == 6) ichambOK = 6;
456 if (fgTrackReconstructor->GetTrackMethod() == 3 &&
457 fSkipHit->GetHitNumber() < 0) {
458 iz0 = fgCombi->IZfromHit(fSkipHit);
461 else currIndx = fgHitForRec->IndexOf(fSkipHit);
464 fTrackHits->Compress();
466 } // if (hit == fSkipHit)
467 else if (currIndx < 0) currIndx = fTrackHits->IndexOf(hit);
468 } // else if (fRecover)
470 // Get indices of the 1'st and last hits on the track candidate
471 firstHit = (AliMUONHitForRec*) fTrackHits->First();
472 lastHit = (AliMUONHitForRec*) fTrackHits->Last();
473 if (fgTrackReconstructor->GetTrackMethod() == 3 &&
474 lastHit->GetHitNumber() < 0) iz0 = fgCombi->IZfromHit(lastHit);
476 firstIndx = fgHitForRec->IndexOf(firstHit);
477 lastIndx = fgHitForRec->IndexOf(lastHit);
478 currIndx = TMath::Abs (TMath::Max(firstIndx*iFB,lastIndx*iFB));
482 if (iz0 < 0) iz0 = iFB;
483 while (ichamb >= iMin && ichamb <= iMax) {
484 // Find the closest hit in Z, not belonging to the current plane
487 if (currIndx < fNmbTrackHits) {
488 hitAdd = (AliMUONHitForRec*) fTrackHits->UncheckedAt(currIndx);
489 zEnd = hitAdd->GetZ();
490 //AZ(z->-z) } else zEnd = -9999;
493 //AZ(Z->-Z) zEnd = -9999;
495 for (ihit=currIndx+iFB; ihit>=0 && ihit<fgNOfPoints; ihit+=iFB) {
496 hitAdd = (AliMUONHitForRec*) ((*fgHitForRec)[ihit]);
497 //if (TMath::Abs(hitAdd->GetZ()-fPosition) > 0.1) {
498 if (TMath::Abs(hitAdd->GetZ()-fPosition) > 0.5) {
499 zEnd = hitAdd->GetZ();
505 // Combined cluster / track finder
506 if (zEnd > 999 && iFB < 0 && fgTrackReconstructor->GetTrackMethod() == 3) {
508 AliMUONHitForRec hitTmp;
509 for (iz = iz0 - iFB; iz < fgCombi->Nz(); iz++) {
510 if (TMath::Abs(fgCombi->Z(iz)-fPosition) < 0.5) continue;
511 Int_t *pDEatZ = fgCombi->DEatZ(iz);
512 //cout << iz << " " << fgCombi->Z(iz) << endl;
513 zEnd = fgCombi->Z(iz);
515 AliMUONDetElement *detElem = fgCombi->DetElem(pDEatZ[0]);
517 hitAdd->SetChamberNumber(detElem->Chamber());
518 //hitAdd = (AliMUONHitForRec*) detElem->HitsForRec()->First();
523 if (zEnd>999 && ichamb==ichamEnd) endOfProp = 1; // end-of-propagation
525 // Check if there is a chamber without hits
526 if (zEnd>999 || TMath::Abs(hitAdd->GetChamberNumber()-ichamb) > 1) {
527 if (!Back && zEnd<999) currIndx -= iFB;
529 zEnd = AliMUONConstants::DefaultChamberZ(ichamb);
532 ichamb = hitAdd->GetChamberNumber();
536 if (ichamb<iMin || ichamb>iMax) break;
537 // Check for missing station
539 dChamb = TMath::Abs(ichamb-ichambOK);
540 //cout << dChamb << " " << ichambOK << " " << fgNOfPoints << endl;
541 Int_t dStatMiss = TMath::Abs (ichamb/2 - ichambOK/2);
542 if (zEnd > 999) dStatMiss++;
544 //if (dStatMiss == 2 && ichambOK/2 != 3 || dStatMiss > 2) { // AZ - missing st. 3
545 // missing station - abandon track
546 //cout << dChamb << " " << ichambOK << " " << fgNOfPoints << " " << 1/(*fTrackPar)(4,0) << endl;
548 for (Int_t i1=0; i1<fgNOfPoints; i1++) {
549 cout << " Hit #" << ((AliMUONHitForRec*)((*fgHitForRec)[i1]))->GetChamberNumber() << " ";
550 cout << ((AliMUONHitForRec*)((*fgHitForRec)[i1]))->GetBendingCoor() << " ";
551 cout << ((AliMUONHitForRec*)((*fgHitForRec)[i1]))->GetNonBendingCoor() << " ";
552 cout << ((AliMUONHitForRec*)((*fgHitForRec)[i1]))->GetZ() << " " << " ";
553 cout << ((AliMUONHitForRec*)((*fgHitForRec)[i1]))->GetTTRTrack() << endl;
556 cout << fNmbTrackHits << endl;
557 for (Int_t i1=0; i1<fNmbTrackHits; i1++) {
558 hit = (AliMUONHitForRec*) ((*fTrackHits)[i1]);
559 printf(" * %d %10.4f %10.4f %10.4f",
560 hit->GetChamberNumber(), hit->GetBendingCoor(),
561 hit->GetNonBendingCoor(), hit->GetZ());
563 rawclusters = fgTrackReconstructor->GetMUONData()->RawClusters(hit->GetChamberNumber());
564 clus = (AliMUONRawCluster*) rawclusters->UncheckedAt(hit->GetHitNumber());
565 printf("%3d", clus->GetTrack(1)-1);
566 if (clus->GetTrack(2) != 0)
567 printf("%3d \n", clus->GetTrack(2)-1);
572 } // if (fgDebug >= 10)
573 if (fNmbTrackHits>2 && fRecover==0) Recover(); // try to recover track later
575 } // if (dStatMiss > 1)
577 if (endOfProp != 0) break;
579 // propagate to the found Z
581 // Check if track steps into dipole
582 //AZ(z->-z) if (fPosition>zDipole2 && zEnd<zDipole2) {
583 if (fPosition<zDipole2 && zEnd>zDipole2) {
584 //LinearPropagation(zDipole2-zBeg);
585 ParPropagation(zDipole2);
586 MSThin(1); // multiple scattering in the chamber
587 WeightPropagation(zDipole2, kTRUE); // propagate weight matrix
588 fPosition = fPositionNew;
589 *fTrackPar = *fTrackParNew;
590 //MagnetPropagation(zEnd);
591 ParPropagation(zEnd);
592 WeightPropagation(zEnd, kTRUE);
593 fPosition = fPositionNew;
595 // Check if track steps out of dipole
596 //AZ(z->-z) else if (fPosition>zDipole1 && zEnd<zDipole1) {
597 else if (fPosition<zDipole1 && zEnd>zDipole1) {
598 //MagnetPropagation(zDipole1-zBeg);
599 ParPropagation(zDipole1);
600 MSThin(1); // multiple scattering in the chamber
601 WeightPropagation(zDipole1, kTRUE);
602 fPosition = fPositionNew;
603 *fTrackPar = *fTrackParNew;
604 //LinearPropagation(zEnd-zDipole1);
605 ParPropagation(zEnd);
606 WeightPropagation(zEnd, kTRUE);
607 fPosition = fPositionNew;
609 ParPropagation(zEnd);
610 //MSThin(1); // multiple scattering in the chamber
611 if (TMath::Abs(zEnd-fPosition) > 5) MSThin(1); // multiple scattering in the chamber
612 WeightPropagation(zEnd, kTRUE);
613 fPosition = fPositionNew;
617 if (fRecover != 0 && hitAdd == fSkipHit && !miss) {
618 // recovered track - remove the hit
620 ichamb = fSkipHit->GetChamberNumber();
621 // remove the skipped hit
622 fTrackHits->Remove(fSkipHit);
624 //AZ fSkipHit->SetNTrackHits(fSkipHit->GetNTrackHits()-1); // unmark hit
627 ichambOK = ((AliMUONHitForRec*)((*fTrackHits)[fNmbTrackHits-1]))->GetChamberNumber();
628 //? if (fgTrackReconstructor->GetTrackMethod() != 3) currIndx = fgHitForRec->IndexOf(fSkipHit);
629 currIndx = fgHitForRec->IndexOf(fSkipHit);
633 // backward propagator
635 TMatrixD pointWeight(fgkSize,fgkSize);
636 TMatrixD point(fgkSize,1);
637 TMatrixD trackParTmp = point;
638 point(0,0) = hitAdd->GetBendingCoor();
639 point(1,0) = hitAdd->GetNonBendingCoor();
640 pointWeight(0,0) = 1/hitAdd->GetBendingReso2();
641 pointWeight(1,1) = 1/hitAdd->GetNonBendingReso2();
642 TryPoint(point,pointWeight,trackParTmp,dChi2);
643 *fTrackPar = trackParTmp;
644 *fWeight += pointWeight;
645 if (fTrackDir > 0) AddMatrices (this, dChi2, hitAdd);
646 fChi2 += dChi2; // Chi2
647 } else *fTrackPar = *fTrackParNew; // adjust end point to the last hit
648 if (ichamb==ichamEnd) break;
651 // forward propagator
652 if (miss || !FindPoint(ichamb, zEnd, currIndx, iFB, hitAdd, iz)) {
654 *fTrackPar = *fTrackParNew;
657 fTrackHits->Add(hitAdd); // add hit
659 hitAdd->SetNTrackHits(hitAdd->GetNTrackHits()+1); // mark hit as being on track
661 currIndx = fgHitForRec->IndexOf(hitAdd); // Check
665 if (fgDebug > 0) cout << fNmbTrackHits << " " << fChi2 << " " << 1/(*fTrackPar)(4,0) << " " << fPosition << endl;
666 if (1/TMath::Abs((*fTrackPar)(4,0)) < fgTrackReconstructor->GetMinBendingMomentum()) success = kFALSE; // p < p_min
670 //__________________________________________________________________________
671 void AliMUONTrackK::ParPropagation(Double_t zEnd)
673 /// Propagation of track parameters to zEnd
675 Double_t dZ, step, distance, charge;
676 Double_t vGeant3[7], vGeant3New[7];
677 AliMUONTrackParam trackParam;
680 // First step using linear extrapolation
681 dZ = zEnd - fPosition;
682 fPositionNew = fPosition;
683 *fTrackParNew = *fTrackPar;
684 if (dZ == 0) return; //AZ ???
685 iFB = (Int_t)TMath::Sign(Double_t(1.0),dZ);
686 step = dZ/TMath::Cos((*fTrackPar)(2,0))/TMath::Cos((*fTrackPar)(3,0)); // linear estimate
687 //AZ(z->-z) charge = iFB*TMath::Sign(Double_t(1.0),(*fTrackPar)(4,0));
688 charge = -iFB*TMath::Sign(Double_t(1.0),(*fTrackPar)(4,0));
689 SetGeantParam(vGeant3,iFB);
690 //fTrackParNew->Print();
694 step = TMath::Abs(step);
695 // Propagate parameters
696 trackParam.ExtrapOneStepRungekutta(charge,step,vGeant3,vGeant3New);
697 //extrap_onestep_rungekutta(charge,step,vGeant3,vGeant3New);
698 distance = zEnd - vGeant3New[2];
699 step *= dZ/(vGeant3New[2]-fPositionNew);
701 } while (distance*iFB < 0 && TMath::Abs(distance) > fgkEpsilon);
703 GetFromGeantParam(vGeant3New,iFB);
704 //fTrackParNew->Print();
706 // Position adjustment (until within tolerance)
707 while (TMath::Abs(distance) > fgkEpsilon) {
708 dZ = zEnd - fPositionNew;
709 iFB = (Int_t)TMath::Sign(Double_t(1.0),dZ);
710 step = dZ/TMath::Cos((*fTrackParNew)(2,0))/TMath::Cos((*fTrackParNew)(3,0));
711 step = TMath::Abs(step);
712 SetGeantParam(vGeant3,iFB);
715 // Propagate parameters
716 trackParam.ExtrapOneStepRungekutta(charge,step,vGeant3,vGeant3New);
717 //extrap_onestep_rungekutta(charge,step,vGeant3,vGeant3New);
718 distance = zEnd - vGeant3New[2];
721 if (nTries > fgkTriesMax) AliError(Form(" Too many tries: %d", nTries));
722 } while (distance*iFB < 0);
724 GetFromGeantParam(vGeant3New,iFB);
726 //cout << nTries << endl;
727 //fTrackParNew->Print();
731 //__________________________________________________________________________
732 void AliMUONTrackK::WeightPropagation(Double_t zEnd, Bool_t smooth)
734 /// Propagation of the weight matrix
735 /// W = DtWD, where D is Jacobian
739 TMatrixD jacob(fgkSize,fgkSize);
742 // Save initial and propagated parameters
743 TMatrixD trackPar0 = *fTrackPar;
744 TMatrixD trackParNew0 = *fTrackParNew;
746 // Get covariance matrix
747 *fCovariance = *fWeight;
748 // check whether the Invert method returns flag if matrix cannot be inverted,
749 // and do not calculate the Determinant in that case !!!!
750 if (fCovariance->Determinant() != 0) {
752 mnvertLocal(&((*fCovariance)(0,0)), fgkSize,fgkSize,fgkSize,ifail);
753 //fCovariance->Print();
755 AliWarning(" Determinant fCovariance=0:");
758 // Loop over parameters to find change of the propagated vs initial ones
759 for (i=0; i<fgkSize; i++) {
760 dPar = TMath::Sqrt((*fCovariance)(i,i));
761 *fTrackPar = trackPar0;
762 if (i == 4) dPar *= TMath::Sign(1.,-trackPar0(4,0)); // 1/p
763 (*fTrackPar)(i,0) += dPar;
764 ParPropagation(zEnd);
765 for (j=0; j<fgkSize; j++) {
766 jacob(j,i) = ((*fTrackParNew)(j,0)-trackParNew0(j,0))/dPar;
770 //trackParNew0.Print();
771 //TMatrixD par1(jacob,TMatrixD::kMult,trackPar0); //
773 TMatrixD jacob0 = jacob;
774 if (jacob.Determinant() != 0) {
777 AliWarning(" Determinant jacob=0:");
779 TMatrixD weight1(*fWeight,TMatrixD::kMult,jacob); // WD
780 *fWeight = TMatrixD(jacob,TMatrixD::kTransposeMult,weight1); // DtWD
782 // Restore initial and propagated parameters
783 *fTrackPar = trackPar0;
784 *fTrackParNew = trackParNew0;
787 if (!smooth) return; // do not use smoother
788 if (fTrackDir < 0) return; // only for propagation towards int. point
789 TMatrixD *tmp = new TMatrixD(*fTrackParNew); // extrapolated parameters
790 fParExtrap->Add(tmp);
792 tmp = new TMatrixD(*fTrackParNew); // filtered parameters (if no measurement will be found)
793 fParFilter->Add(tmp);
795 *fCovariance = *fWeight;
796 if (fCovariance->Determinant() != 0) {
798 mnvertLocal(&((*fCovariance)(0,0)), fgkSize,fgkSize,fgkSize,ifail);
800 AliWarning(" Determinant fCovariance=0:");
802 tmp = new TMatrixD(*fCovariance); // extrapolated covariance
803 fCovExtrap->Add(tmp);
805 tmp = new TMatrixD(*fCovariance); // filtered covariance (if no measurement will be found)
806 fCovFilter->Add(tmp);
808 tmp = new TMatrixD(jacob0); // Jacobian
811 if (fSteps->fN <= fNSteps) fSteps->Set(fSteps->fN+10);
812 fSteps->AddAt(fPositionNew,fNSteps++);
813 if (fgDebug > 0) cout << " WeightPropagation " << fNSteps << " " << fPositionNew << endl;
817 //__________________________________________________________________________
818 Bool_t AliMUONTrackK::FindPoint(Int_t ichamb, Double_t zEnd, Int_t currIndx, Int_t iFB, AliMUONHitForRec *&hitAdd, Int_t iz)
820 /// Picks up point within a window for the chamber No ichamb
821 /// Split the track if there are more than 1 hit
822 Int_t ihit, nRecTracks;
823 Double_t windowB, windowNonB, dChi2Tmp=0, dChi2, y, x, savePosition=0;
824 TClonesArray *trackPtr;
825 AliMUONHitForRec *hit, *hitLoop;
826 AliMUONTrackK *trackK;
827 AliMUONDetElement *detElem = NULL;
831 if (fgTrackReconstructor->GetTrackMethod() == 3 && iz >= 0) {
832 // Combined cluster / track finder
833 // Check if extrapolated track passes thru det. elems. at iz
834 Int_t *pDEatZ = fgCombi->DEatZ(iz);
835 Int_t nDetElem = pDEatZ[-1];
836 //cout << fgCombi->Z(iz) << " " << nDetElem << endl;
837 for (Int_t i = 0; i < nDetElem; i++) {
838 detElem = fgCombi->DetElem(pDEatZ[i]);
839 if (detElem->Inside((*fTrackParNew)(1,0), (*fTrackParNew)(0,0), fPosition)) {
840 detElem->ClusterReco((*fTrackParNew)(1,0), (*fTrackParNew)(0,0));
841 hitAdd = (AliMUONHitForRec*) detElem->HitsForRec()->First();
846 if (!ok) return ok; // outside det. elem.
850 //sigmaB = fgTrackReconstructor->GetBendingResolution(); // bending resolution
851 //sigmaNonB = fgTrackReconstructor->GetNonBendingResolution(); // non-bending resolution
852 *fCovariance = *fWeight;
853 // check whether the Invert method returns flag if matrix cannot be inverted,
854 // and do not calculate the Determinant in that case !!!!
855 if (fCovariance->Determinant() != 0) {
857 mnvertLocal(&((*fCovariance)(0,0)), fgkSize,fgkSize,fgkSize,ifail);
859 AliWarning(" Determinant fCovariance=0:");
861 //windowB = fgkNSigma*TMath::Sqrt((*fCovariance)(0,0)+sigmaB*sigmaB);
862 //windowNonB = fgkNSigma*TMath::Sqrt((*fCovariance)(1,1)+sigmaNonB*sigmaNonB);
863 // Loop over all hits and take hits from the chamber
864 TMatrixD pointWeight(fgkSize,fgkSize);
865 TMatrixD saveWeight = pointWeight;
866 TMatrixD pointWeightTmp = pointWeight;
867 TMatrixD point(fgkSize,1);
868 TMatrixD trackPar = point;
869 TMatrixD trackParTmp = point;
870 Int_t nHitsOK = 0, ihitB = currIndx, ihitE = fgNOfPoints, iDhit = iFB;
872 zLast = ((AliMUONHitForRec*)fTrackHits->Last())->GetZ();
873 if (fgTrackReconstructor->GetTrackMethod() == 3 && detElem) {
875 ihitE = detElem->NHitsForRec();
880 for (ihit = ihitB; ihit >= 0 && ihit < ihitE; ihit+=iDhit) {
881 if (fgTrackReconstructor->GetTrackMethod() != 3 || !detElem)
882 hit = (AliMUONHitForRec*) ((*fgHitForRec)[ihit]);
883 else hit = (AliMUONHitForRec*) detElem->HitsForRec()->UncheckedAt(ihit);
884 if (hit->GetChamberNumber() == ichamb) {
885 //if (TMath::Abs(hit->GetZ()-zEnd) < 0.1) {
886 if (TMath::Abs(hit->GetZ()-zEnd) < 0.5) {
887 y = hit->GetBendingCoor();
888 x = hit->GetNonBendingCoor();
889 if (hit->GetBendingReso2() < 0) {
890 // Combined cluster / track finder
891 hit->SetBendingReso2(fgTrackReconstructor->GetBendingResolution()*
892 fgTrackReconstructor->GetBendingResolution());
893 hit->SetNonBendingReso2(fgTrackReconstructor->GetNonBendingResolution()*
894 fgTrackReconstructor->GetNonBendingResolution());
896 windowB = fgkNSigma*TMath::Sqrt((*fCovariance)(0,0)+hit->GetBendingReso2());
897 windowNonB = fgkNSigma*TMath::Sqrt((*fCovariance)(1,1)+hit->GetNonBendingReso2());
899 // windowB = TMath::Min (windowB,5.);
900 if (fgkNSigma > 6) windowB = TMath::Min (windowB,5.);
902 if (TMath::Abs(hit->GetZ()-zLast) < 50) {
903 windowB = TMath::Min (windowB,0.5);
904 windowNonB = TMath::Min (windowNonB,3.);
905 } else if (TMath::Abs(hit->GetZ()-zLast) < 200) {
906 windowB = TMath::Min (windowB,1.5);
907 windowNonB = TMath::Min (windowNonB,3.);
908 } else if (TMath::Abs(hit->GetZ()-zLast) < 350) {
909 windowB = TMath::Min (windowB,4.);
910 windowNonB = TMath::Min (windowNonB,6.);
916 if (TMath::Abs(hit->GetZ()-zLast) < 50) {
917 windowB = 5*TMath::Sqrt((*fCovariance)(0,0)+hit->GetBendingReso2());
918 windowNonB = 5*TMath::Sqrt((*fCovariance)(1,1)+hit->GetNonBendingReso2());
920 windowB = fgkNSigma*TMath::Sqrt((*fCovariance)(0,0)+hit->GetBendingReso2());
921 windowNonB = fgkNSigma*TMath::Sqrt((*fCovariance)(1,1)+hit->GetNonBendingReso2());
925 //cout << TMath::Abs((*fTrackParNew)(0,0)-y) << " " << windowB << " " << TMath::Abs((*fTrackParNew)(1,0)-x) << " " << windowNonB << " " << zEnd << endl;
926 if (TMath::Abs((*fTrackParNew)(0,0)-y) <= windowB &&
927 TMath::Abs((*fTrackParNew)(1,0)-x) <= windowNonB) {
928 //if (TMath::Abs((*fTrackParNew)(0,0)-y) <= windowB &&
929 // TMath::Abs((*fTrackParNew)(1,0)-x) <= windowNonB &&
930 // hit->GetTrackRefSignal() == 1) { // just for test
931 // Vector of measurements and covariance matrix
932 //fprintf(lun1,"%3d %3d %10.4f %10.4f \n", gAlice->GetEvNumber(), ichamb, x, y);
933 if (TMath::Abs(hit->GetZ()-zEnd) > 0.05) {
934 // Adjust position: for multiple hits in the chamber or misalignment (Z as a function of X or Y)
935 //AliWarning(Form(" *** adjust %f %f ", zEnd, hit->GetZ()));
937 *fTrackPar = *fTrackParNew;
938 ParPropagation(zEnd);
939 WeightPropagation(zEnd, kTRUE);
940 fPosition = fPositionNew;
941 *fTrackPar = *fTrackParNew;
943 *fCovariance = *fWeight;
944 if (fCovariance->Determinant() != 0) {
946 mnvertLocal(&((*fCovariance)(0,0)), fgkSize,fgkSize,fgkSize,ifail);
948 AliWarning(" Determinant fCovariance=0:" );
954 pointWeight(0,0) = 1/hit->GetBendingReso2();
955 pointWeight(1,1) = 1/hit->GetNonBendingReso2();
956 TryPoint(point,pointWeight,trackPar,dChi2);
957 if (TMath::Abs(1./(trackPar)(4,0)) < fgTrackReconstructor->GetMinBendingMomentum()) continue; // p < p_min - next hit
958 // if (TMath::Sign(1.,(trackPar)(4,0)*(*fTrackPar)(4,0)) < 0) continue; // change of sign
961 //if (nHitsOK > -1) {
963 // Save current members
964 saveWeight = *fWeight;
965 savePosition = fPosition;
966 // temporary storage for the current track
968 trackParTmp = trackPar;
969 pointWeightTmp = pointWeight;
971 if (fgDebug > 0) cout << " Added point: " << x << " " << y << " " << dChi2 << endl;
973 // branching: create a new track
974 trackPtr = fgTrackReconstructor->GetRecTracksPtr();
975 nRecTracks = fgTrackReconstructor->GetNRecTracks();
976 trackK = new ((*trackPtr)[nRecTracks]) AliMUONTrackK(NULL, NULL);
978 fgTrackReconstructor->SetNRecTracks(nRecTracks+1);
979 if (fgDebug > 0) cout << " ******** New track: " << ichamb << " " << hit->GetTTRTrack() << " " << 1/(trackPar)(4,0) << " " << hit->GetBendingCoor() << " " << hit->GetNonBendingCoor() << " " << fNmbTrackHits << " " << nRecTracks << endl;
980 trackK->fRecover = 0;
981 *(trackK->fTrackPar) = trackPar;
982 *(trackK->fWeight) += pointWeight;
983 trackK->fChi2 += dChi2;
986 *(trackK->fCovariance) = *(trackK->fWeight);
987 if (trackK->fCovariance->Determinant() != 0) {
989 mnvertLocal(&((*(trackK->fCovariance))(0,0)), fgkSize,fgkSize,fgkSize,ifail);
991 cout << (*(trackK->fCovariance))(0,0) << " " << (*(trackK->fCovariance))(1,1) << " " << (*fCovariance)(0,0) << " " << (*fCovariance)(1,1) << endl;
993 // Add filtered matrices for the smoother
995 if (nHitsOK > 2) { // check for position adjustment
996 for (Int_t i=trackK->fNSteps-1; i>=0; i--) {
997 if (TMath::Abs(hit->GetZ()-(*trackK->fSteps)[i]) > 0.1) {
998 //if (TMath::Abs(hit->GetZ()-(trackK->fSteps)[i]) > 0.1) {
999 RemoveMatrices(trackK);
1000 AliError(Form(" *** Position adjustment 1: %f %f", hit->GetZ(), (*trackK->fSteps)[i]));
1005 AddMatrices (trackK, dChi2, hit);
1007 // Mark hits as being on 2 tracks
1008 for (Int_t i=0; i<fNmbTrackHits; i++) {
1009 hitLoop = (AliMUONHitForRec*) ((*fTrackHits)[i]);
1010 //AZ hitLoop->SetNTrackHits(hitLoop->GetNTrackHits()+1);
1013 cout << hitLoop->GetChamberNumber() << " ";
1014 cout << hitLoop->GetBendingCoor() << " ";
1015 cout << hitLoop->GetNonBendingCoor() << " ";
1016 cout << hitLoop->GetZ() << " " << " ";
1017 cout << hitLoop->GetTTRTrack() << endl;
1018 printf(" ** %d %10.4f %10.4f %10.4f\n",
1019 hitLoop->GetChamberNumber(), hitLoop->GetBendingCoor(),
1020 hitLoop->GetNonBendingCoor(), hitLoop->GetZ());
1024 trackK->fTrackHits->Add(hit); // add hit
1025 trackK->fNmbTrackHits ++;
1026 hit->SetNTrackHits(hit->GetNTrackHits()+1); // mark hit as being on track
1029 trackK->fTrackDir = 1;
1030 trackK->fBPFlag = kTRUE;
1035 } else break; // different chamber
1036 } // for (ihit=currIndx;
1039 *fTrackPar = trackParTmp;
1040 *fWeight = saveWeight;
1041 *fWeight += pointWeightTmp;
1042 fChi2 += dChi2Tmp; // Chi2
1043 fPosition = savePosition;
1044 // Add filtered matrices for the smoother
1045 if (fTrackDir > 0) {
1046 for (Int_t i=fNSteps-1; i>=0; i--) {
1047 if (TMath::Abs(fPosition-(*fSteps)[i]) > 0.1) {
1048 //if (TMath::Abs(fPosition-fSteps[i]) > 0.1) {
1049 RemoveMatrices(this);
1050 if (fgDebug > 0) cout << " *** Position adjustment 2 " << fPosition << " " << (*fSteps)[i] << endl;
1053 } // for (Int_t i=fNSteps-1;
1054 AddMatrices (this, dChi2Tmp, hitAdd);
1057 for (Int_t i=0; i<fNSteps; i++) cout << (*fSteps)[i] << " "; cout << endl;
1058 for (Int_t i=0; i<trackK->fNSteps; i++) cout << (*trackK->fSteps)[i] << " "; cout << endl;
1061 } // if (fTrackDir > 0)
1066 //__________________________________________________________________________
1067 void AliMUONTrackK::TryPoint(TMatrixD &point, const TMatrixD &pointWeight, TMatrixD &trackParTmp, Double_t &dChi2)
1069 /// Adds a measurement point (modifies track parameters and computes
1072 // Solving linear system (W+U)p' = U(m-p) + (W+U)p
1073 TMatrixD wu = *fWeight;
1074 wu += pointWeight; // W+U
1075 trackParTmp = point;
1076 trackParTmp -= *fTrackParNew; // m-p
1077 TMatrixD right(pointWeight,TMatrixD::kMult,trackParTmp); // U(m-p)
1078 TMatrixD right1(wu,TMatrixD::kMult,*fTrackParNew); // (W+U)p
1079 right += right1; // U(m-p) + (W+U)p
1081 // check whether the Invert method returns flag if matrix cannot be inverted,
1082 // and do not calculate the Determinant in that case !!!!
1083 if (wu.Determinant() != 0) {
1085 mnvertLocal(&((wu)(0,0)), fgkSize,fgkSize,fgkSize,ifail);
1087 AliWarning(" Determinant wu=0:");
1089 trackParTmp = TMatrixD(wu,TMatrixD::kMult,right);
1091 right1 = trackParTmp;
1092 right1 -= point; // p'-m
1093 point = trackParTmp;
1094 point -= *fTrackParNew; // p'-p
1095 right = TMatrixD(*fWeight,TMatrixD::kMult,point); // W(p'-p)
1096 TMatrixD value(point,TMatrixD::kTransposeMult,right); // (p'-p)'W(p'-p)
1098 right = TMatrixD(pointWeight,TMatrixD::kMult,right1); // U(p'-m)
1099 value = TMatrixD(right1,TMatrixD::kTransposeMult,right); // (p'-m)'U(p'-m)
1100 dChi2 += value(0,0);
1104 //__________________________________________________________________________
1105 void AliMUONTrackK::MSThin(Int_t sign)
1107 /// Adds multiple scattering in a thin layer (only angles are affected)
1108 Double_t cosAlph, cosBeta, momentum, velo, path, theta0;
1110 // check whether the Invert method returns flag if matrix cannot be inverted,
1111 // and do not calculate the Determinant in that case !!!!
1112 if (fWeight->Determinant() != 0) {
1114 mnvertLocal(&((*fWeight)(0,0)), fgkSize,fgkSize,fgkSize,ifail);
1116 AliWarning(" Determinant fWeight=0:");
1119 cosAlph = TMath::Cos((*fTrackParNew)(2,0));
1120 cosBeta = TMath::Cos((*fTrackParNew)(3,0));
1121 momentum = 1/(*fTrackParNew)(4,0); // particle momentum
1122 //velo = momentum/TMath::Sqrt(momentum*momentum+muonMass*muonMass); // velocity/c for muon hypothesis
1123 velo = 1; // relativistic
1124 path = TMath::Abs(fgTrackReconstructor->GetChamberThicknessInX0()/cosAlph/cosBeta); // path length
1125 theta0 = 0.0136/velo/momentum*TMath::Sqrt(path)*(1+0.038*TMath::Log(path)); // projected scattering angle
1127 (*fWeight)(2,2) += sign*theta0/cosBeta*theta0/cosBeta; // alpha
1128 (*fWeight)(3,3) += sign*theta0*theta0; // beta
1130 mnvertLocal(&((*fWeight)(0,0)), fgkSize,fgkSize,fgkSize,ifail);
1134 //__________________________________________________________________________
1135 void AliMUONTrackK::StartBack(void)
1137 /// Starts backpropagator
1141 for (Int_t i=0; i<fgkSize; i++) {
1142 for (Int_t j=0; j<fgkSize; j++) {
1143 if (j==i) (*fWeight)(i,i) /= 100;
1144 //if (j==i) (*fWeight)(i,i) /= fNmbTrackHits*fNmbTrackHits;
1145 else (*fWeight)(j,i) = 0;
1148 // Sort hits on track in descending order in abs(z)
1149 SortHits(0, fTrackHits);
1152 //__________________________________________________________________________
1153 void AliMUONTrackK::SortHits(Int_t iflag, TObjArray *array)
1155 /// Sort hits in Z if the seed segment in the last but one station
1156 /// (if iflag==0 in descending order in abs(z), if !=0 - unsort)
1158 if (iflag && ((AliMUONHitForRec*)(array->UncheckedAt(0)))->GetChamberNumber() == 6) return;
1159 Double_t z = 0, zmax = TMath::Abs(((AliMUONHitForRec*)(array->UncheckedAt(0)))->GetZ());
1160 Int_t i = 1, entries = array->GetEntriesFast();
1161 for ( ; i<entries; i++) {
1163 if (((AliMUONHitForRec*)(array->UncheckedAt(i)))->GetChamberNumber() == 6) break;
1165 z = TMath::Abs(((AliMUONHitForRec*)(array->UncheckedAt(i)))->GetZ());
1166 if (z < zmax) break;
1168 if (fgDebug >= 10) cout << " " << zmax << " " << z << " " << i << endl;
1172 for (Int_t j=0; j<=(i-1)/2; j++) {
1173 TObject *hit = array->UncheckedAt(j);
1174 array->AddAt(array->UncheckedAt(i-j),j);
1175 array->AddAt(hit,i-j);
1177 if (fgDebug >= 10) {
1178 for (i=0; i<entries; i++)
1179 cout << ((AliMUONHitForRec*)(array->UncheckedAt(i)))->GetZ() << " ";
1180 cout << " - Sort" << endl;
1184 //__________________________________________________________________________
1185 void AliMUONTrackK::SetGeantParam(Double_t *VGeant3, Int_t iFB)
1187 /// Set vector of Geant3 parameters pointed to by "VGeant3"
1188 /// from track parameters
1190 VGeant3[0] = (*fTrackParNew)(1,0); // X
1191 VGeant3[1] = (*fTrackParNew)(0,0); // Y
1192 VGeant3[2] = fPositionNew; // Z
1193 VGeant3[3] = iFB*TMath::Sin((*fTrackParNew)(3,0)); // Px/Ptot
1194 VGeant3[4] = iFB*TMath::Cos((*fTrackParNew)(3,0))*TMath::Sin((*fTrackParNew)(2,0)); // Py/Ptot
1195 VGeant3[5] = iFB*TMath::Sqrt(1.0-VGeant3[3]*VGeant3[3]-VGeant3[4]*VGeant3[4]); // Pz/Ptot
1196 VGeant3[6] = 1/TMath::Abs((*fTrackParNew)(4,0)); // Ptot
1199 //__________________________________________________________________________
1200 void AliMUONTrackK::GetFromGeantParam(Double_t *VGeant3, Int_t iFB)
1202 /// Get track parameters from vector of Geant3 parameters pointed
1205 fPositionNew = VGeant3[2]; // Z
1206 (*fTrackParNew)(0,0) = VGeant3[1]; // Y
1207 (*fTrackParNew)(1,0) = VGeant3[0]; // X
1208 (*fTrackParNew)(3,0) = TMath::ASin(iFB*VGeant3[3]); // beta
1209 (*fTrackParNew)(2,0) = TMath::ASin(iFB*VGeant3[4]/TMath::Cos((*fTrackParNew)(3,0))); // alpha
1210 (*fTrackParNew)(4,0) = 1/VGeant3[6]*TMath::Sign(Double_t(1.0),(*fTrackPar)(4,0)); // 1/Ptot
1213 //__________________________________________________________________________
1214 void AliMUONTrackK::SetTrackQuality(Int_t iChi2)
1216 /// Computes "track quality" from Chi2 (if iChi2==0) or vice versa
1219 AliWarning(Form(" *** Too high Chi2: %f ", fChi2));
1222 if (iChi2 == 0) fChi2 = fNmbTrackHits + (500.-fChi2)/501;
1223 else fChi2 = 500 - (fChi2-fNmbTrackHits)*501;
1226 //__________________________________________________________________________
1227 Int_t AliMUONTrackK::Compare(const TObject* trackK) const
1229 /// "Compare" function to sort with decreasing "track quality".
1230 /// Returns +1 (0, -1) if quality of current track
1231 /// is smaller than (equal to, larger than) quality of trackK
1233 if (fChi2 < ((AliMUONTrackK*)trackK)->fChi2) return(+1);
1234 else if (fChi2 == ((AliMUONTrackK*)trackK)->fChi2) return(0);
1238 //__________________________________________________________________________
1239 Bool_t AliMUONTrackK::KeepTrack(AliMUONTrackK* track0) const
1241 /// Check whether or not to keep current track
1242 /// (keep, if it has less than half of common hits with track0)
1243 Int_t hitsInCommon, nHits0, i, j, nTrackHits2;
1244 AliMUONHitForRec *hit0, *hit1;
1247 nHits0 = track0->fNmbTrackHits;
1248 nTrackHits2 = fNmbTrackHits/2;
1250 for (i=0; i<nHits0; i++) {
1251 // Check if hit belongs to several tracks
1252 hit0 = (AliMUONHitForRec*) (*track0->fTrackHits)[i];
1253 if (hit0->GetNTrackHits() == 1) continue;
1254 for (j=0; j<fNmbTrackHits; j++) {
1255 hit1 = (AliMUONHitForRec*) (*fTrackHits)[j];
1256 if (hit1->GetNTrackHits() == 1) continue;
1259 if (hitsInCommon >= nTrackHits2) return kFALSE;
1267 //__________________________________________________________________________
1268 void AliMUONTrackK::Kill(void)
1270 /// Kill track candidate
1271 fgTrackReconstructor->GetRecTracksPtr()->Remove(this);
1274 //__________________________________________________________________________
1275 void AliMUONTrackK::FillMUONTrack(void)
1277 /// Compute track parameters at hit positions (as for AliMUONTrack)
1279 // Set number of hits per track
1280 SetNTrackHits(fNmbTrackHits);
1284 // Set track parameters at vertex
1285 AliMUONTrackParam trackParam;
1286 SetTrackParam(&trackParam, fTrackPar, fPosition);
1287 SetTrackParamAtVertex(&trackParam);
1289 // Set track parameters at hits
1290 for (Int_t i = fNmbTrackHits-1; i>=0; i--) {
1291 if ((*fChi2Smooth)[i] < 0) {
1292 // Propagate through last chambers
1293 trackParam.ExtrapToZ(((AliMUONHitForRec*)((*fTrackHits)[i]))->GetZ());
1296 SetTrackParam(&trackParam, (TMatrixD*)fParSmooth->UncheckedAt(i), ((AliMUONHitForRec*)((*fTrackHits)[i]))->GetZ());
1298 AddTrackParamAtHit(&trackParam);
1299 // Fill array of HitForRec's
1300 AddHitForRecAtHit((AliMUONHitForRec*)fTrackHits->UncheckedAt(i));
1304 //__________________________________________________________________________
1305 void AliMUONTrackK::SetTrackParam(AliMUONTrackParam *trackParam, TMatrixD *par, Double_t z)
1307 /// Fill AliMUONTrackParam object
1309 trackParam->SetBendingCoor((*par)(0,0));
1310 trackParam->SetNonBendingCoor((*par)(1,0));
1311 trackParam->SetBendingSlope(TMath::Tan((*par)(2,0)));
1312 trackParam->SetNonBendingSlope(TMath::Tan((*par)(3,0))/TMath::Cos((*par)(2,0)));
1313 trackParam->SetInverseBendingMomentum((*par)(4,0)/TMath::Cos((*par)(3,0)));
1314 trackParam->SetZ(z);
1317 //__________________________________________________________________________
1318 void AliMUONTrackK::Branson(void)
1320 /// Propagates track to the vertex thru absorber using Branson correction
1321 /// (makes use of the AliMUONTrackParam class)
1323 //AliMUONTrackParam *trackParam = new AliMUONTrackParam();
1324 AliMUONTrackParam trackParam = AliMUONTrackParam();
1326 trackParam->SetBendingCoor((*fTrackPar)(0,0));
1327 trackParam->SetNonBendingCoor((*fTrackPar)(1,0));
1328 trackParam->SetBendingSlope(TMath::Tan((*fTrackPar)(2,0)));
1329 trackParam->SetNonBendingSlope(TMath::Tan((*fTrackPar)(3,0))/TMath::Cos((*fTrackPar)(2,0)));
1330 trackParam->SetInverseBendingMomentum((*fTrackPar)(4,0)/TMath::Cos((*fTrackPar)(3,0)));
1331 trackParam->SetZ(fPosition);
1333 SetTrackParam(&trackParam, fTrackPar, fPosition);
1335 trackParam.ExtrapToVertex(Double_t(0.), Double_t(0.), Double_t(0.));
1336 //trackParam.ExtrapToVertex();
1338 (*fTrackPar)(0,0) = trackParam.GetBendingCoor();
1339 (*fTrackPar)(1,0) = trackParam.GetNonBendingCoor();
1340 (*fTrackPar)(2,0) = TMath::ATan(trackParam.GetBendingSlope());
1341 (*fTrackPar)(3,0) = TMath::ATan(TMath::Cos((*fTrackPar)(2,0))*trackParam.GetNonBendingSlope());
1342 (*fTrackPar)(4,0) = TMath::Cos((*fTrackPar)(3,0))*trackParam.GetInverseBendingMomentum();
1343 fPosition = trackParam.GetZ();
1344 //delete trackParam;
1345 if (fgDebug > 0) cout << 1/(*fTrackPar)(4,0) << " " << fPosition << " " << (*fTrackPar)(0,0) << endl;
1347 // Get covariance matrix
1348 *fCovariance = *fWeight;
1349 if (fCovariance->Determinant() != 0) {
1351 mnvertLocal(&((*fCovariance)(0,0)), fgkSize,fgkSize,fgkSize,ifail);
1353 AliWarning(" Determinant fCovariance=0:");
1357 //__________________________________________________________________________
1358 void AliMUONTrackK::GoToZ(Double_t zEnd)
1360 /// Propagates track to given Z
1362 ParPropagation(zEnd);
1363 MSThin(1); // multiple scattering in the chamber
1364 WeightPropagation(zEnd, kFALSE);
1365 fPosition = fPositionNew;
1366 *fTrackPar = *fTrackParNew;
1369 //__________________________________________________________________________
1370 void AliMUONTrackK::GoToVertex(Int_t iflag)
1373 /// Propagates track to the vertex
1374 /// All material constants are taken from AliRoot
1376 static Double_t x01[5] = { 24.282, // C
1381 // inner part theta < 3 degrees
1382 static Double_t x02[5] = { 30413, // Air
1387 // z positions of the materials inside the absober outer part theta > 3 degres
1388 static Double_t zPos[10] = {-90, -105, -315, -443, -468};
1390 Double_t dZ, r0Norm, x0, deltaP, dChi2, pTotal, pOld;
1391 AliMUONHitForRec *hit;
1392 AliMUONRawCluster *clus;
1393 TClonesArray *rawclusters;
1395 // First step to the rear end of the absorber
1396 Double_t zRear = -503;
1398 Double_t tan3 = TMath::Tan(3./180*TMath::Pi());
1400 // Go through absorber
1401 pOld = 1/(*fTrackPar)(4,0);
1402 Double_t r0Rear = (*fTrackPar)(0,0)*(*fTrackPar)(0,0) +
1403 (*fTrackPar)(1,0)*(*fTrackPar)(1,0);
1404 r0Rear = TMath::Sqrt(r0Rear)/TMath::Abs(fPosition)/tan3;
1406 Double_t p0, cos25, cos60;
1407 if (!iflag) goto vertex;
1409 for (Int_t i=4; i>=0; i--) {
1410 ParPropagation(zPos[i]);
1411 WeightPropagation(zPos[i], kFALSE);
1412 dZ = TMath::Abs (fPositionNew-fPosition);
1413 if (r0Norm > 1) x0 = x01[i];
1415 MSLine(dZ,x0); // multiple scattering in the medium (linear approximation)
1416 fPosition = fPositionNew;
1417 *fTrackPar = *fTrackParNew;
1418 r0Norm = (*fTrackPar)(0,0)*(*fTrackPar)(0,0) +
1419 (*fTrackPar)(1,0)*(*fTrackPar)(1,0);
1420 r0Norm = TMath::Sqrt(r0Norm)/TMath::Abs(fPosition)/tan3;
1422 // Correct momentum for energy losses
1423 pTotal = 1/TMath::Abs((*fTrackPar)(4,0));
1425 cos25 = TMath::Cos(2.5/180*TMath::Pi());
1426 cos60 = TMath::Cos(6.0/180*TMath::Pi());
1427 for (Int_t j=0; j<1; j++) {
1431 deltaP = 2.164 + 0.145e-1*p0 - 0.417e-3*p0*p0;
1433 deltaP = 2.275 + 0.102e-2*p0 - 0.674e-6*p0*p0;
1437 deltaP = 2.581 + 0.188e-1*p0 - 0.398e-3*p0*p0;
1439 deltaP = 2.727 + 0.356e-2*p0 + 0.242e-5*p0*p0;
1447 deltaP = 2.737 + 0.0494*p0 - 0.001123*p0*p0;
1449 deltaP = 3.0643 + 0.01346*p0;
1455 deltaP = 2.1380 + 0.0351*p0 - 0.000853*p0*p0;
1457 deltaP = 2.407 + 0.00702*p0;
1466 deltaP = 2.439 + 0.806e-1*p0 - 0.500e-2*p0*p0 + 0.106e-3*p0*p0*p0;
1468 deltaP = 2.767 + 0.742e-2*p0 - 0.196e-4*p0*p0 + 0.403e-7*p0*p0*p0;
1475 deltaP = 2.209 + 0.800e-2*p0;
1477 deltaP = 2.285 + 0.141e-2*p0 - 0.446e-6*p0*p0;
1487 deltaP = 2.5938 + 0.0570 * p0 - 0.001151 * p0 * p0;
1489 deltaP = 3.0714 + 0.011767 * p0;
1494 deltaP = 2.1207 + 0.05478 * p0 - 0.00145079 * p0 * p0;
1496 deltaP = 2.6069 + 0.0051705 * p0;
1502 p0 = pTotal + deltaP/TMath::Abs(TMath::Cos((*fTrackPar)(2,0))/TMath::Cos((*fTrackPar)(3,0)));
1504 (*fTrackPar)(4,0) = 1/p0*TMath::Sign((Double_t)1.,(*fTrackPar)(4,0));
1508 ParPropagation((Double_t)0.);
1509 WeightPropagation((Double_t)0., kFALSE);
1510 fPosition = fPositionNew;
1511 //*fTrackPar = *fTrackParNew;
1512 // Add vertex as a hit
1513 TMatrixD pointWeight(fgkSize,fgkSize);
1514 TMatrixD point(fgkSize,1);
1515 TMatrixD trackParTmp = point;
1516 point(0,0) = 0; // vertex coordinate - should be taken somewhere
1517 point(1,0) = 0; // vertex coordinate - should be taken somewhere
1518 pointWeight(0,0) = 1/1.e-3/1.e-3; // 10 um error
1519 pointWeight(1,1) = 1/1.e-3/1.e-3; // 10 um error
1520 TryPoint(point,pointWeight,trackParTmp,dChi2);
1521 *fTrackPar = trackParTmp;
1522 *fWeight += pointWeight;
1523 fChi2 += dChi2; // Chi2
1524 if (fgDebug < 0) return; // no output
1526 cout << pOld << " " << 1/(*fTrackPar)(4,0) << " " << dChi2 << " " << fChi2 << " " << fNmbTrackHits << endl;
1527 for (Int_t i1=0; i1<fNmbTrackHits; i1++) {
1528 hit = (AliMUONHitForRec*) ((*fTrackHits)[i1]);
1529 printf ("%5d", hit->GetChamberNumber());
1533 for (Int_t i1=0; i1<fNmbTrackHits; i1++) {
1534 hit = (AliMUONHitForRec*) ((*fTrackHits)[i1]);
1535 //cout << ((AliMUONHitForRec*)((*fTrackHits)[i1]))->GetHitNumber() << " ";
1536 //cout << ((AliMUONHitForRec*)((*fTrackHits)[i1]))->GetZ() << " ";
1537 printf ("%5d", fgHitForRec->IndexOf(hit));
1542 // from raw clusters
1543 for (Int_t i1=0; i1<fNmbTrackHits; i1++) {
1544 hit = (AliMUONHitForRec*) ((*fTrackHits)[i1]);
1545 if (hit->GetHitNumber() < 0) { // combined cluster / track finder
1546 Int_t index = -hit->GetHitNumber() / 100000;
1547 Int_t iPos = -hit->GetHitNumber() - index * 100000;
1548 clus = (AliMUONRawCluster*) fgCombi->DetElem(index-1)->RawClusters()->UncheckedAt(iPos);
1550 rawclusters = fgTrackReconstructor->GetMUONData()->RawClusters(hit->GetChamberNumber());
1551 clus = (AliMUONRawCluster*) rawclusters->UncheckedAt(hit->GetHitNumber());
1553 printf ("%5d", clus->GetTrack(1)%10000000);
1556 for (Int_t i1=0; i1<fNmbTrackHits; i1++) {
1557 hit = (AliMUONHitForRec*) ((*fTrackHits)[i1]);
1558 if (hit->GetHitNumber() < 0) { // combined cluster / track finder
1559 Int_t index = -hit->GetHitNumber() / 100000;
1560 Int_t iPos = -hit->GetHitNumber() - index * 100000;
1561 clus = (AliMUONRawCluster*) fgCombi->DetElem(index-1)->RawClusters()->UncheckedAt(iPos);
1563 rawclusters = fgTrackReconstructor->GetMUONData()->RawClusters(hit->GetChamberNumber());
1564 clus = (AliMUONRawCluster*) rawclusters->UncheckedAt(hit->GetHitNumber());
1566 if (clus->GetTrack(2) != -1) printf ("%5d", clus->GetTrack(2)%10000000);
1567 else printf ("%5s", " ");
1571 for (Int_t i1=0; i1<fNmbTrackHits; i1++) {
1572 //cout << ((AliMUONHitForRec*)((*fTrackHits)[i1]))->GetHitNumber() << " ";
1573 cout << ((AliMUONHitForRec*)((*fTrackHits)[i1]))->GetZ() << " ";
1574 //cout << fgHitForRec->IndexOf(((AliMUONHitForRec*)((*fTrackHits)[i1]))) << " ";
1577 for (Int_t i1=0; i1<fNmbTrackHits; i1++) printf("%8.4f", (*fChi2Smooth)[i1]);
1579 cout << "---------------------------------------------------" << endl;
1581 // Get covariance matrix
1582 /* Not needed - covariance matrix is not interesting to anybody
1583 *fCovariance = *fWeight;
1584 if (fCovariance->Determinant() != 0) {
1586 mnvertLocal(&((*fCovariance)(0,0)), fgkSize,fgkSize,fgkSize,ifail);
1588 AliWarning(" Determinant fCovariance=0:" );
1593 //__________________________________________________________________________
1594 void AliMUONTrackK::MSLine(Double_t dZ, Double_t x0)
1596 /// Adds multiple scattering in a thick layer for linear propagation
1598 Double_t cosAlph = TMath::Cos((*fTrackPar)(2,0));
1599 Double_t tanAlph = TMath::Tan((*fTrackPar)(2,0));
1600 Double_t cosBeta = TMath::Cos((*fTrackPar)(3,0));
1602 sinBeta = TMath::Sin((*fTrackPar)(3,0));
1603 Double_t tanBeta = TMath::Tan((*fTrackPar)(3,0));
1604 Double_t momentum = 1/(*fTrackPar)(4,0);
1605 Double_t velo = 1; // relativistic velocity
1606 Double_t step = TMath::Abs(dZ/cosAlph/cosBeta); // step length
1608 // Projected scattering angle
1609 Double_t theta0 = 0.0136/velo/momentum/TMath::Sqrt(x0)*(1+0.038*TMath::Log(step/x0));
1610 Double_t theta02 = theta0*theta0;
1611 Double_t dl2 = step*step/2*theta02;
1612 Double_t dl3 = dl2*step*2/3;
1615 Double_t dYdT = 1/cosAlph;
1616 Double_t dYdB = 0; //(*fTrackPar)(2,0)*sinBeta/cosAlph;
1617 Double_t dXdT = tanAlph*tanBeta;
1618 //Double_t dXdB = (1+(*fTrackPar)(2,0)*tanAlph*sinBeta*sinBeta)/cosBeta;
1619 Double_t dXdB = 1/cosBeta;
1620 Double_t dAdT = 1/cosBeta;
1621 Double_t dAdB = 0; //(*fTrackPar)(2,0)*tanBeta;
1623 // Get covariance matrix
1624 *fCovariance = *fWeight;
1625 if (fCovariance->Determinant() != 0) {
1626 // fCovariance->Invert();
1628 mnvertLocal(&((*fCovariance)(0,0)), fgkSize,fgkSize,fgkSize,ifail);
1630 AliWarning(" Determinant fCovariance=0:" );
1633 (*fCovariance)(0,0) += dl3*(dYdT*dYdT+dYdB*dYdB); // <yy>
1634 (*fCovariance)(1,1) += dl3*(dXdT*dXdT+dXdB*dXdB); // <xx>
1635 (*fCovariance)(2,2) += theta02*step*(dAdT*dAdT+dAdB*dAdB); // <aa>
1636 (*fCovariance)(3,3) += theta02*step; // <bb>
1638 (*fCovariance)(0,1) += dl3*(dYdT*dXdT+dYdB*dXdB); // <yx>
1639 (*fCovariance)(1,0) = (*fCovariance)(0,1);
1641 (*fCovariance)(0,2) += dl2*(dYdT*dAdT+dYdB*dAdB); // <ya>
1642 (*fCovariance)(2,0) = (*fCovariance)(0,2);
1644 (*fCovariance)(0,3) += dl2*dYdB; // <yb>
1645 (*fCovariance)(3,0) = (*fCovariance)(0,3);
1647 (*fCovariance)(1,2) += dl2*(-dXdT*dAdT+dXdB*dAdB); // <xa>
1648 (*fCovariance)(2,1) = (*fCovariance)(1,2);
1650 (*fCovariance)(1,3) += dl2*dXdB; // <xb>
1651 (*fCovariance)(3,1) = (*fCovariance)(1,3);
1653 (*fCovariance)(2,3) += theta02*step*dAdB; // <ab>
1654 (*fCovariance)(3,2) = (*fCovariance)(2,3);
1656 // Get weight matrix
1657 *fWeight = *fCovariance;
1658 if (fWeight->Determinant() != 0) {
1660 mnvertLocal(&((*fWeight)(0,0)), fgkSize,fgkSize,fgkSize,ifail);
1662 AliWarning(" Determinant fWeight=0:");
1666 //__________________________________________________________________________
1667 Bool_t AliMUONTrackK::Recover(void)
1669 /// Adds new failed track(s) which can be tried to be recovered
1671 TClonesArray *trackPtr;
1672 AliMUONTrackK *trackK;
1674 if (fgDebug > 0) cout << " ******** Enter Recover " << endl;
1675 trackPtr = fgTrackReconstructor->GetRecTracksPtr();
1677 // Remove hit with the highest chi2
1680 for (Int_t i=0; i<fNmbTrackHits; i++) {
1681 chi2 = fChi2Smooth ? (*fChi2Smooth)[i] : (*fChi2Array)[i];
1682 printf("%10.4f", chi2);
1685 for (Int_t i=0; i<fNmbTrackHits; i++) {
1686 printf("%10d", ((AliMUONHitForRec*)fTrackHits->UncheckedAt(i))->GetChamberNumber());
1690 Double_t chi2max = 0;
1692 for (Int_t i=0; i<fNmbTrackHits; i++) {
1693 chi2 = fChi2Smooth ? (*fChi2Smooth)[i] : (*fChi2Array)[i];
1694 if (chi2 < chi2max) continue;
1698 //if (chi2max < 10) return kFALSE; // !!!
1699 //if (chi2max < 25) imax = fNmbTrackHits - 1;
1700 if (chi2max < 15) imax = fNmbTrackHits - 1; // discard the last point
1701 // Check if the outlier is not from the seed segment
1702 AliMUONHitForRec *skipHit = (AliMUONHitForRec*) fTrackHits->UncheckedAt(imax);
1703 if (skipHit == fStartSegment->GetHitForRec1() || skipHit == fStartSegment->GetHitForRec2()) {
1704 //DropBranches(fStartSegment); // drop all tracks with the same seed segment
1705 return kFALSE; // to be changed probably
1708 // Make a copy of track hit collection
1709 TObjArray *hits = new TObjArray(*fTrackHits);
1713 // Hits after the found one will be removed
1714 if (GetStation0() == 3 && skipHit->GetChamberNumber() >= 7) {
1715 SortHits(1, fTrackHits); // unsort hits
1716 imax = fTrackHits->IndexOf(skipHit);
1718 Int_t nTrackHits = fNmbTrackHits;
1719 for (Int_t i=imax+1; i<nTrackHits; i++) {
1720 AliMUONHitForRec *hit = (AliMUONHitForRec*) fTrackHits->UncheckedAt(i);
1721 fTrackHits->Remove(hit);
1722 hit->SetNTrackHits(hit->GetNTrackHits()-1); // unmark hit
1726 // Check if the track candidate doesn't exist yet
1727 if (ExistDouble()) { delete hits; return kFALSE; }
1729 //DropBranches(imax0, hits); // drop branches downstream the discarded hit
1732 nRecTracks = fgTrackReconstructor->GetNRecTracks();
1733 skipHit = (AliMUONHitForRec*) ((*fTrackHits)[fNmbTrackHits-1]);
1734 // Remove all saved steps and smoother matrices after the skipped hit
1735 RemoveMatrices(skipHit->GetZ());
1737 //AZ(z->-z) if (skipHit->GetZ() > fStartSegment->GetHitForRec2()->GetZ() || !fNSteps) {
1738 if (TMath::Abs(skipHit->GetZ()) > TMath::Abs(fStartSegment->GetHitForRec2()->GetZ()) || !fNSteps) {
1739 // Propagation toward high Z or skipped hit next to segment -
1740 // start track from segment
1741 trackK = new ((*trackPtr)[nRecTracks]) AliMUONTrackK(fStartSegment);
1742 fgTrackReconstructor->SetNRecTracks(nRecTracks+1);
1743 trackK->fRecover = 1;
1744 trackK->fSkipHit = skipHit;
1745 trackK->fNmbTrackHits = fNmbTrackHits;
1746 delete trackK->fTrackHits; // not efficient ?
1747 trackK->fTrackHits = new TObjArray(*fTrackHits);
1748 if (fgDebug > 0) cout << nRecTracks << " " << trackK->fRecover << endl;
1752 trackK = new ((*trackPtr)[nRecTracks]) AliMUONTrackK(NULL, NULL);
1754 fgTrackReconstructor->SetNRecTracks(nRecTracks+1);
1755 //AZ(z->-z) trackK->fTrackDir = -1;
1756 trackK->fTrackDir = 1;
1757 trackK->fRecover = 1;
1758 trackK->fSkipHit = skipHit;
1759 Int_t iD = trackK->fParFilter->Last()->GetUniqueID();
1761 trackK->fParFilter->Last()->SetUniqueID(iD-1);
1762 CreateMatrix(trackK->fParFilter);
1764 *((TMatrixD*)trackK->fParFilter->Last()) = *((TMatrixD*)fParExtrap->Last());
1765 trackK->fParFilter->Last()->SetUniqueID(1);
1766 *(trackK->fTrackPar) = *((TMatrixD*)trackK->fParFilter->Last());
1767 iD = trackK->fCovFilter->Last()->GetUniqueID();
1769 trackK->fCovFilter->Last()->SetUniqueID(iD-1);
1770 CreateMatrix(trackK->fCovFilter);
1772 *((TMatrixD*)trackK->fCovFilter->Last()) = *((TMatrixD*)fCovExtrap->Last());
1773 trackK->fCovFilter->Last()->SetUniqueID(1);
1774 *(trackK->fWeight) = *((TMatrixD*)trackK->fCovFilter->Last());
1775 if (trackK->fWeight->Determinant() != 0) {
1777 mnvertLocal(&((*(trackK->fWeight))(0,0)), fgkSize,fgkSize,fgkSize,ifail);
1779 AliWarning(" Determinant fWeight=0:");
1781 trackK->fPosition = trackK->fPositionNew = (*fSteps)[fNSteps-1];
1783 for (Int_t i=0; i<fNmbTrackHits-1; i++) trackK->fChi2 += (*fChi2Array)[i];
1784 if (fgDebug > 0) cout << nRecTracks << " " << trackK->fRecover << endl;
1788 //__________________________________________________________________________
1789 void AliMUONTrackK::AddMatrices(AliMUONTrackK *trackK, Double_t dChi2, AliMUONHitForRec *hitAdd)
1791 /// Adds matrices for the smoother and keep Chi2 for the point
1792 /// Track parameters
1793 //trackK->fParFilter->Last()->Print();
1794 Int_t iD = trackK->fParFilter->Last()->GetUniqueID();
1796 trackK->fParFilter->Last()->SetUniqueID(iD-1);
1797 CreateMatrix(trackK->fParFilter);
1800 *((TMatrixD*)(trackK->fParFilter->Last())) = *(trackK->fTrackPar);
1801 trackK->fParFilter->Last()->SetUniqueID(iD);
1803 cout << " Add matrices" << " " << fPosition << " " << dChi2 << " " << fgHitForRec->IndexOf(hitAdd) << endl;
1804 //trackK->fTrackPar->Print();
1805 //trackK->fTrackParNew->Print();
1806 trackK->fParFilter->Last()->Print();
1807 cout << " Add matrices" << endl;
1810 *(trackK->fCovariance) = *(trackK->fWeight);
1811 if (trackK->fCovariance->Determinant() != 0) {
1813 mnvertLocal(&((*(trackK->fCovariance))(0,0)), fgkSize,fgkSize,fgkSize,ifail);
1815 AliWarning(" Determinant fCovariance=0:");
1817 iD = trackK->fCovFilter->Last()->GetUniqueID();
1819 trackK->fCovFilter->Last()->SetUniqueID(iD-1);
1820 CreateMatrix(trackK->fCovFilter);
1823 *((TMatrixD*)(trackK->fCovFilter->Last())) = *(trackK->fCovariance);
1824 trackK->fCovFilter->Last()->SetUniqueID(iD);
1826 // Save Chi2-increment for point
1827 Int_t indx = trackK->fTrackHits->IndexOf(hitAdd);
1828 if (indx < 0) indx = fNmbTrackHits;
1829 if (trackK->fChi2Array->fN <= indx) trackK->fChi2Array->Set(indx+10);
1830 trackK->fChi2Array->AddAt(dChi2,indx);
1833 //__________________________________________________________________________
1834 void AliMUONTrackK::CreateMatrix(TObjArray *objArray) const
1836 /// Create new matrix and add it to TObjArray
1838 TMatrixD *matrix = (TMatrixD*) objArray->First();
1839 TMatrixD *tmp = new TMatrixD(*matrix);
1840 objArray->AddAtAndExpand(tmp,objArray->GetLast());
1841 //cout << " !!! " << tmp << " " << tmp->GetUniqueID() << " " << tmp->GetNoElements() << endl;
1844 //__________________________________________________________________________
1845 void AliMUONTrackK::RemoveMatrices(Double_t zEnd)
1847 /// Remove matrices (and saved steps) in the smoother part with abs(z) < abs(zEnd)
1849 for (Int_t i=fNSteps-1; i>=0; i--) {
1850 if (fgDebug > 1) printf("%10.4f %10.4f \n", (*fSteps)[i], zEnd);
1851 if (TMath::Abs((*fSteps)[i]) > TMath::Abs(zEnd) - 0.01) break;
1852 RemoveMatrices(this);
1853 } // for (Int_t i=fNSteps-1;
1856 //__________________________________________________________________________
1857 void AliMUONTrackK::RemoveMatrices(AliMUONTrackK* trackK)
1859 /// Remove last saved matrices and steps in the smoother part
1862 Int_t i = trackK->fNSteps;
1865 // Delete only matrices not shared by several tracks
1866 id = trackK->fParExtrap->Last()->GetUniqueID();
1868 trackK->fParExtrap->Last()->SetUniqueID(id-1);
1869 trackK->fParExtrap->RemoveAt(i);
1871 else ((TMatrixD*)(trackK->fParExtrap->RemoveAt(i)))->Delete();
1872 id = fParFilter->Last()->GetUniqueID();
1874 trackK->fParFilter->Last()->SetUniqueID(id-1);
1875 trackK->fParFilter->RemoveAt(i);
1877 else ((TMatrixD*)(trackK->fParFilter->RemoveAt(i)))->Delete();
1878 id = trackK->fCovExtrap->Last()->GetUniqueID();
1880 trackK->fCovExtrap->Last()->SetUniqueID(id-1);
1881 trackK->fCovExtrap->RemoveAt(i);
1883 else ((TMatrixD*)(trackK->fCovExtrap->RemoveAt(i)))->Delete();
1884 id = trackK->fCovFilter->Last()->GetUniqueID();
1886 trackK->fCovFilter->Last()->SetUniqueID(id-1);
1887 trackK->fCovFilter->RemoveAt(i);
1889 else ((TMatrixD*)(trackK->fCovFilter->RemoveAt(i)))->Delete();
1890 id = trackK->fJacob->Last()->GetUniqueID();
1892 trackK->fJacob->Last()->SetUniqueID(id-1);
1893 trackK->fJacob->RemoveAt(i);
1895 else ((TMatrixD*)(trackK->fJacob->RemoveAt(i)))->Delete();
1898 //__________________________________________________________________________
1899 Bool_t AliMUONTrackK::Smooth(void)
1902 Int_t ihit = fNmbTrackHits - 1;
1903 AliMUONHitForRec *hit = (AliMUONHitForRec*) (*fTrackHits)[ihit];
1904 fChi2Smooth = new TArrayD(fNmbTrackHits);
1905 fChi2Smooth->Reset(-1);
1907 fParSmooth = new TObjArray(15);
1908 fParSmooth->Clear();
1911 cout << " ******** Enter Smooth " << endl;
1912 cout << (*fSteps)[fNSteps-1] << " " << fPosition << " " << fPositionNew << endl;
1914 for (Int_t i=fNSteps-1; i>=0; i--) cout << (*fSteps)[i] << " ";
1916 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();}
1918 for (Int_t i=ihit; i>=0; i--) cout << ((AliMUONHitForRec*)(*fTrackHits)[i])->GetZ() << " ";
1922 // Find last point corresponding to the last hit
1923 Int_t iLast = fNSteps - 1;
1924 for ( ; iLast>=0; iLast--) {
1925 //AZ(z->-z) if ((*fSteps)[iLast] + 0.01 > ((AliMUONHitForRec*)(*fTrackHits)[fNmbTrackHits-1])->GetZ()) break;
1926 if (TMath::Abs((*fSteps)[iLast]) + 0.01 > TMath::Abs(((AliMUONHitForRec*)(*fTrackHits)[fNmbTrackHits-1])->GetZ())) break;
1929 TMatrixD parSmooth = *((TMatrixD*)(fParFilter->UncheckedAt(iLast))); // last filtered = first smoothed
1931 TMatrixD covSmooth = *((TMatrixD*)(fCovFilter->UncheckedAt(iLast))); // last filtered = first smoothed
1932 TMatrixD tmp=covSmooth, weight=covSmooth, gain = covSmooth;
1933 TMatrixD tmpPar = *fTrackPar;
1934 //parSmooth.Print(); ((TMatrixD*)(fParExtrap->Last()))->Print();
1937 Double_t chi2max = 0;
1938 for (Int_t i=iLast+1; i>0; i--) {
1939 if (i == iLast + 1) goto L33; // temporary fix
1941 // Smoother gain matrix
1942 weight = *((TMatrixD*)(fCovExtrap->UncheckedAt(i)));
1943 if (weight.Determinant() != 0) {
1945 mnvertLocal(&(weight(0,0)), fgkSize,fgkSize,fgkSize,ifail);
1947 AliWarning(" Determinant weight=0:");
1950 tmp = TMatrixD(*((TMatrixD*)(fJacob->UncheckedAt(i))),TMatrixD::kTransposeMult,weight);
1951 gain = TMatrixD(*((TMatrixD*)(fCovFilter->UncheckedAt(i-1))),TMatrixD::kMult,tmp);
1952 //cout << (*fSteps)[i-1] << " " << (*fSteps)[i] << endl; gain.Print();
1954 // Smoothed parameter vector
1955 //((TMatrixD*)(fParExtrap->UncheckedAt(i)))->Dump();
1956 parSmooth -= *((TMatrixD*)(fParExtrap->UncheckedAt(i)));
1957 tmpPar = *((TMatrixD*)(fParFilter->UncheckedAt(i-1)));
1958 tmpPar += TMatrixD(gain,TMatrixD::kMult,parSmooth);
1961 // Smoothed covariance
1962 covSmooth -= *((TMatrixD*)(fCovExtrap->UncheckedAt(i)));
1963 weight = TMatrixD(TMatrixD::kTransposed,gain);
1964 tmp = TMatrixD(covSmooth,TMatrixD::kMult,weight);
1965 covSmooth = *((TMatrixD*)(fCovFilter->UncheckedAt(i-1)));
1966 covSmooth += TMatrixD(gain,TMatrixD::kMult,tmp);
1968 // Check if there was a measurement at given z
1970 for ( ; ihit>=0; ihit--) {
1971 hit = (AliMUONHitForRec*) (*fTrackHits)[ihit];
1972 if (TMath::Abs(hit->GetZ()-(*fSteps)[i-1]) < 0.1) { found = kTRUE; break; }
1973 //AZ(z->-z) else if (ihit < fNmbTrackHits-1 && hit->GetZ() > (*fSteps)[i-1]) { ihit++; break; }
1974 else if (ihit < fNmbTrackHits-1 && TMath::Abs(hit->GetZ()) > TMath::Abs((*fSteps)[i-1])) { ihit++; break; }
1976 if (!found) continue; // no measurement - skip the rest
1977 else if (fgDebug > 1) cout << i << " " << ihit << " " << hit->GetZ() << endl;
1978 if (ihit == 0) continue; // the first hit - skip the rest
1981 // Smoothed residuals
1983 tmpPar(0,0) = hit->GetBendingCoor() - parSmooth(0,0); // bending coordinate
1984 tmpPar(1,0) = hit->GetNonBendingCoor() - parSmooth(1,0); // non-bending coordinate
1986 cout << hit->GetBendingCoor() << " " << parSmooth(0,0) << " " << tmpPar(0,0) << endl;
1987 cout << hit->GetNonBendingCoor() << " " << parSmooth(1,0) << " " << tmpPar(1,0) << endl;
1989 // Cov. matrix of smoothed residuals
1991 tmp(0,0) = hit->GetBendingReso2() - covSmooth(0,0);
1992 tmp(1,1) = hit->GetNonBendingReso2() - covSmooth(1,1);
1993 tmp(0,1) = tmp(1,0) = -covSmooth(0,1);
1994 tmp(2,2) = tmp(3,3) = tmp(4,4) = 1;
1996 // Calculate Chi2 of smoothed residuals
1997 if (tmp.Determinant() != 0) {
1999 mnvertLocal(&(tmp(0,0)), fgkSize,fgkSize,fgkSize,ifail);
2001 AliWarning(" Determinant tmp=0:");
2003 TMatrixD vector(tmp,TMatrixD::kMult,tmpPar);
2004 TMatrixD chi2(tmpPar,TMatrixD::kTransposeMult,vector);
2005 if (fgDebug > 1) chi2.Print();
2006 (*fChi2Smooth)[ihit] = chi2(0,0);
2007 if (chi2(0,0) > chi2max) chi2max = chi2(0,0);
2009 if (chi2(0,0) < 0) {
2011 AliError(Form(" *** chi2 < 0: %d %d ", i, iLast));
2013 // Save smoothed parameters
2014 TMatrixD *par = new TMatrixD(parSmooth);
2015 fParSmooth->AddAtAndExpand(par, ihit);
2017 } // for (Int_t i=iLast+1;
2019 //if (chi2max > 16) {
2020 //if (chi2max > 25) {
2021 //if (chi2max > 50) {
2022 //if (chi2max > 100) {
2023 if (chi2max > fgkChi2max) {
2024 //if (Recover()) DropBranches();
2032 //__________________________________________________________________________
2033 void AliMUONTrackK::Outlier()
2035 /// Adds new track with removed hit having the highest chi2
2038 cout << " ******** Enter Outlier " << endl;
2039 for (Int_t i=0; i<fNmbTrackHits; i++) printf("%10.4f", (*fChi2Smooth)[i]);
2041 for (Int_t i=0; i<fNmbTrackHits; i++) {
2042 printf("%10d", ((AliMUONHitForRec*)fTrackHits->UncheckedAt(i))->GetChamberNumber());
2047 Double_t chi2max = 0;
2049 for (Int_t i=0; i<fNmbTrackHits; i++) {
2050 if ((*fChi2Smooth)[i] < chi2max) continue;
2051 chi2max = (*fChi2Smooth)[i];
2054 // Check if the outlier is not from the seed segment
2055 AliMUONHitForRec *hit = (AliMUONHitForRec*) fTrackHits->UncheckedAt(imax);
2056 if (hit == fStartSegment->GetHitForRec1() || hit == fStartSegment->GetHitForRec2()) return; // to be changed probably
2058 // Check for missing station
2061 if (((AliMUONHitForRec*)fTrackHits->UncheckedAt(1))->GetChamberNumber() < 8) ok--;
2062 } else if (imax == fNmbTrackHits-1) {
2063 if (((AliMUONHitForRec*)fTrackHits->UncheckedAt(fNmbTrackHits-2))->GetChamberNumber() > 1) ok--;
2065 else if (((AliMUONHitForRec*)fTrackHits->UncheckedAt(imax-1))->GetChamberNumber()/2 - ((AliMUONHitForRec*)fTrackHits->UncheckedAt(imax+1))->GetChamberNumber()/2 > 1) ok--;
2066 if (!ok) { Recover(); return; } // try to recover track
2067 //AZ if (!ok) { if (fgDebug >= 0) cout << imax << endl; DropBranches(imax, 0); return; }
2069 // Remove saved steps and smoother matrices after the outlier
2070 RemoveMatrices(hit->GetZ());
2072 // Check for possible double track candidates
2073 //if (ExistDouble(hit)) return;
2075 TClonesArray *trackPtr = fgTrackReconstructor->GetRecTracksPtr();
2076 Int_t nRecTracks = fgTrackReconstructor->GetNRecTracks();
2078 AliMUONTrackK *trackK = 0;
2079 if (!fNSteps || GetStation0() == 3 && hit->GetChamberNumber() > 7) {
2080 // start track from segment
2081 trackK = new ((*trackPtr)[nRecTracks]) AliMUONTrackK(fStartSegment);
2082 fgTrackReconstructor->SetNRecTracks(nRecTracks+1);
2083 trackK->fRecover = 2;
2084 trackK->fSkipHit = hit;
2085 trackK->fNmbTrackHits = fNmbTrackHits;
2087 hit = (AliMUONHitForRec*) trackK->fTrackHits->UncheckedAt(0);
2088 hit->SetNTrackHits(hit->GetNTrackHits()-1);
2089 hit = (AliMUONHitForRec*) trackK->fTrackHits->UncheckedAt(1);
2090 hit->SetNTrackHits(hit->GetNTrackHits()-1);
2091 delete trackK->fTrackHits; // not efficient ?
2092 trackK->fTrackHits = new TObjArray(*fTrackHits);
2093 for (Int_t i = 0; i < fNmbTrackHits; i++) {
2094 hit = (AliMUONHitForRec*) fTrackHits->UncheckedAt(i);
2095 hit->SetNTrackHits(hit->GetNTrackHits()+1);
2098 if (GetStation0() == 3) trackK->SortHits(1, trackK->fTrackHits);
2099 if (fgDebug > 0) cout << nRecTracks << " Outlier" << endl;
2102 trackK = new ((*trackPtr)[nRecTracks]) AliMUONTrackK(NULL, NULL);
2104 fgTrackReconstructor->SetNRecTracks(nRecTracks+1);
2105 trackK->fTrackDir = 1;
2106 trackK->fRecover = 2;
2107 trackK->fSkipHit = hit;
2108 Int_t iD = trackK->fParFilter->Last()->GetUniqueID();
2110 trackK->fParFilter->Last()->SetUniqueID(iD-1);
2111 CreateMatrix(trackK->fParFilter);
2113 *((TMatrixD*)trackK->fParFilter->Last()) = *((TMatrixD*)fParExtrap->Last());
2114 trackK->fParFilter->Last()->SetUniqueID(1);
2115 *(trackK->fTrackPar) = *((TMatrixD*)trackK->fParFilter->Last());
2116 iD = trackK->fCovFilter->Last()->GetUniqueID();
2118 trackK->fCovFilter->Last()->SetUniqueID(iD-1);
2119 CreateMatrix(trackK->fCovFilter);
2121 *((TMatrixD*)trackK->fCovFilter->Last()) = *((TMatrixD*)fCovExtrap->Last());
2122 trackK->fCovFilter->Last()->SetUniqueID(1);
2123 *(trackK->fWeight) = *((TMatrixD*)trackK->fCovFilter->Last());
2124 if (trackK->fWeight->Determinant() != 0) {
2126 mnvertLocal(&((*(trackK->fWeight))(0,0)), fgkSize,fgkSize,fgkSize,ifail);
2128 AliWarning(" Determinant fWeight=0:");
2130 trackK->fPosition = trackK->fPositionNew = (*fSteps)[fNSteps-1];
2132 for (Int_t i=0; i<fNmbTrackHits-1; i++) trackK->fChi2 += (*fChi2Array)[i];
2133 if (fgDebug > 0) cout << nRecTracks << " Outlier " << trackK->fChi2 << endl;
2136 //__________________________________________________________________________
2137 Double_t AliMUONTrackK::GetChi2PerPoint(Int_t iPoint) const
2139 /// Return Chi2 at point
2140 return fChi2Smooth ? (*fChi2Smooth)[iPoint] : (*fChi2Array)[iPoint];
2144 //__________________________________________________________________________
2145 void AliMUONTrackK::Print(FILE *lun) const
2147 /// Print out track information
2150 AliMUONHitForRec *hit = 0;
2151 // from raw clusters
2152 AliMUONRawCluster *clus = 0;
2153 TClonesArray *rawclusters = 0;
2154 for (Int_t i1=0; i1<fNmbTrackHits; i1++) {
2155 hit = (AliMUONHitForRec*) ((*fTrackHits)[i1]);
2156 rawclusters = fgTrackReconstructor->GetMUONData()->RawClusters(hit->GetChamberNumber());
2157 clus = (AliMUONRawCluster*) rawclusters->UncheckedAt(hit->GetHitNumber());
2158 if (TMath::Abs(clus->GetTrack(1)-1) < 2) {
2159 if (clus->GetTrack(2)) flag = 2;
2162 if (clus->GetTrack(2) && TMath::Abs(clus->GetTrack(2)-1) < 2) {
2169 Int_t sig[2]={1,1}, tid[2]={0};
2170 for (Int_t i1=0; i1<fNmbTrackHits; i1++) {
2171 if (GetChi2PerPoint(i1) < -0.1) continue;
2172 hit = (AliMUONHitForRec*) ((*fTrackHits)[i1]);
2173 rawclusters = fgTrackReconstructor->GetMUONData()->RawClusters(hit->GetChamberNumber());
2174 clus = (AliMUONRawCluster*) rawclusters->UncheckedAt(hit->GetHitNumber());
2175 for (Int_t j=0; j<2; j++) {
2176 tid[j] = clus->GetTrack(j+1) - 1;
2177 if (clus->GetTrack(j+1) < 0) { sig[j] = 0; tid[j] = 999; }
2179 fprintf(lun,"%3d %3d %10.4f", gAlice->GetEvNumber(), hit->GetChamberNumber(), GetChi2PerPoint(i1));
2180 if (!(clus->GetTrack(2))) fprintf(lun, "%3d %3d", sig[0], tid[0]); // simple cluster
2181 else { // track overlap
2182 fprintf(lun, "%3d %3d", TMath::Max(sig[0],sig[1]), TMath::Min(tid[0],tid[1]));
2183 //if (tid[0] < 2) flag *= 2;
2184 //else if (tid[1] < 2) flag *= 3;
2186 fprintf (lun, "%3d \n", flag);
2191 //__________________________________________________________________________
2192 void AliMUONTrackK::DropBranches(Int_t imax, TObjArray *hits)
2194 /// Drop branches downstream of the skipped hit
2196 TClonesArray *trackPtr;
2197 AliMUONTrackK *trackK;
2199 trackPtr = fgTrackReconstructor->GetRecTracksPtr();
2200 nRecTracks = fgTrackReconstructor->GetNRecTracks();
2201 Int_t icand = trackPtr->IndexOf(this);
2202 if (!hits) hits = fTrackHits;
2204 // Check if the track candidate doesn't exist yet
2205 for (Int_t i=icand+1; i<nRecTracks; i++) {
2206 trackK = (AliMUONTrackK*) ((*trackPtr)[i]);
2207 if (trackK->fNmbTrackHits == 2 && trackK->GetRecover() == 0) continue;
2208 if (trackK->GetRecover() < 0) continue;
2210 if (trackK->fNmbTrackHits >= imax + 1) {
2211 for (Int_t j=0; j<=imax; j++) {
2212 //if (j != fNmbTrackHits-1 && (*trackK->fTrackHits)[j] != (*fTrackHits)[j]) break;
2213 if ((*trackK->fTrackHits)[j] != (*hits)[j]) break;
2215 if (hits != fTrackHits) {
2216 // Drop all branches downstream the hit (for Recover)
2217 trackK->SetRecover(-1);
2218 if (fgDebug >= 0 )cout << " Recover " << i << endl;
2221 // Check if the removal of the hit breaks the track
2224 if (((AliMUONHitForRec*)trackK->fTrackHits->UncheckedAt(1))->GetChamberNumber() < 8) ok--; }
2225 else if (imax == trackK->fNmbTrackHits-1) continue;
2226 // else if (imax == trackK->fNmbTrackHits-1) {
2227 //if (((AliMUONHitForRec*)trackK->fTrackHits->UncheckedAt(trackK->fNmbTrackHits-2))->GetChamberNumber() > 1) ok--;
2229 else if (((AliMUONHitForRec*)trackK->fTrackHits->UncheckedAt(imax-1))->GetChamberNumber()/2 - ((AliMUONHitForRec*)trackK->fTrackHits->UncheckedAt(imax+1))->GetChamberNumber()/2 > 1) ok--;
2230 if (!ok) trackK->SetRecover(-1);
2232 } // for (Int_t j=0;
2234 } // for (Int_t i=0;
2237 //__________________________________________________________________________
2238 void AliMUONTrackK::DropBranches(AliMUONSegment *segment)
2240 /// Drop all candidates with the same seed segment
2242 TClonesArray *trackPtr;
2243 AliMUONTrackK *trackK;
2245 trackPtr = fgTrackReconstructor->GetRecTracksPtr();
2246 nRecTracks = fgTrackReconstructor->GetNRecTracks();
2247 Int_t icand = trackPtr->IndexOf(this);
2249 for (Int_t i=icand+1; i<nRecTracks; i++) {
2250 trackK = (AliMUONTrackK*) ((*trackPtr)[i]);
2251 if (trackK->fNmbTrackHits == 2 && trackK->GetRecover() == 0) continue;
2252 if (trackK->GetRecover() < 0) continue;
2253 if (trackK->fStartSegment == segment) trackK->SetRecover(-1);
2255 if (fgDebug >= 0) cout << " Drop segment " << endl;
2258 //__________________________________________________________________________
2259 AliMUONHitForRec* AliMUONTrackK::GetHitLastOk(void)
2261 /// Return the hit where track stopped
2263 if (!fNSteps) return (AliMUONHitForRec*)((*fTrackHits)[1]);
2267 //__________________________________________________________________________
2268 Int_t AliMUONTrackK::GetStation0(void)
2270 /// Return seed station number
2271 return fStartSegment->GetHitForRec1()->GetChamberNumber() / 2;
2274 //__________________________________________________________________________
2275 Bool_t AliMUONTrackK::ExistDouble(AliMUONHitForRec *hit)
2277 /// Check if the track will make a double after outlier removal
2279 TClonesArray *trackPtr = fgTrackReconstructor->GetRecTracksPtr();
2280 Int_t nRecTracks = fgTrackReconstructor->GetNRecTracks();
2281 TObjArray *hitArray = new TObjArray(*fTrackHits);
2282 TObjArray *hitArray1 = new TObjArray(*hitArray);
2283 hitArray1->Remove(hit);
2284 hitArray1->Compress();
2286 Bool_t same = kFALSE;
2287 for (Int_t i=0; i<nRecTracks; i++) {
2288 AliMUONTrackK *trackK = (AliMUONTrackK*) ((*trackPtr)[i]);
2289 if (trackK->fNmbTrackHits == 2 && trackK->GetRecover() == 0) continue;
2290 if (trackK == this) continue;
2291 if (trackK->fNmbTrackHits == fNmbTrackHits || trackK->fNmbTrackHits == fNmbTrackHits-1) {
2292 TObjArray *hits = new TObjArray(*trackK->fTrackHits);
2294 if (trackK->fNmbTrackHits == fNmbTrackHits) {
2295 for (Int_t j=0; j<fNmbTrackHits; j++) {
2296 if (hits->UncheckedAt(j) != hitArray->UncheckedAt(j)) { same = kFALSE; break; }
2298 if (same) { delete hits; break; }
2299 if (trackK->fSkipHit) {
2300 TObjArray *hits1 = new TObjArray(*hits);
2301 if (hits1->Remove(trackK->fSkipHit) > 0) {
2304 for (Int_t j=0; j<fNmbTrackHits-1; j++) {
2305 if (hits1->UncheckedAt(j) != hitArray1->UncheckedAt(j)) { same = kFALSE; break; }
2307 if (same) { delete hits1; break; }
2312 // Check with removed outlier
2314 for (Int_t j=0; j<fNmbTrackHits-1; j++) {
2315 if (hits->UncheckedAt(j) != hitArray1->UncheckedAt(j)) { same = kFALSE; break; }
2317 if (same) { delete hits; break; }
2321 } // for (Int_t i=0; i<nRecTracks;
2322 delete hitArray; delete hitArray1;
2323 if (same && fgDebug >= 0) cout << " Same" << endl;
2327 //__________________________________________________________________________
2328 Bool_t AliMUONTrackK::ExistDouble(void)
2330 /// Check if the track will make a double after recovery
2332 TClonesArray *trackPtr = fgTrackReconstructor->GetRecTracksPtr();
2333 Int_t nRecTracks = fgTrackReconstructor->GetNRecTracks();
2335 TObjArray *hitArray = new TObjArray(*fTrackHits);
2336 if (GetStation0() == 3) SortHits(0, hitArray); // sort
2337 //if (GetStation0() == 3) SortHits(1, hitArray); // unsort
2339 Bool_t same = kFALSE;
2340 for (Int_t i=0; i<nRecTracks; i++) {
2341 AliMUONTrackK *trackK = (AliMUONTrackK*) ((*trackPtr)[i]);
2342 if (trackK->fNmbTrackHits == 2 && trackK->GetRecover() == 0) continue;
2343 if (trackK == this) continue;
2344 //AZ if (trackK->GetRecover() < 0) continue; //
2345 if (trackK->fNmbTrackHits >= fNmbTrackHits) {
2346 TObjArray *hits = new TObjArray(*trackK->fTrackHits);
2347 if (trackK->GetStation0() == 3) SortHits(0, hits); // sort
2348 //if (trackK->GetStation0() == 3) SortHits(1, hits); // unsort
2349 for (Int_t j=0; j<fNmbTrackHits; j++) {
2350 //cout << fNmbTrackHits << " " << i << " " << j << " " << (*hitArray)[j] << " " << (*hits)[j] << " " << trackK->fSkipHit << endl;
2351 if (j != fNmbTrackHits-1 && (*hitArray)[j] != (*hits)[j]) break;
2352 if (j == fNmbTrackHits-1) {
2353 if (trackK->fSkipHit && TMath::Abs(((AliMUONHitForRec*)((*hitArray)[j]))->GetZ()-trackK->fSkipHit->GetZ()) < 0.5) same = kTRUE;
2354 //if (trackK->fNmbTrackHits > fNmbTrackHits &&
2355 //if (trackK->fSkipHit) cout << ((AliMUONHitForRec*)((*hitArray)[j]))->GetZ() << " " << trackK->fSkipHit->GetZ() << endl;
2357 } // for (Int_t j=0;
2361 } // for (Int_t i=0;