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 **************************************************************************/
17 //-------------------------------------------------------
18 // Implementation of the TPC tracker
20 // Origin: Marian Ivanov Marian.Ivanov@cern.ch
22 // AliTPC parallel tracker
24 // The track fitting is based on Kalaman filtering approach
26 // The track finding steps:
27 // 1. Seeding - with and without vertex constraint
28 // - seeding with vertex constain done at first n^2 proble
29 // - seeding without vertex constraint n^3 problem
30 // 2. Tracking - follow prolongation road - find cluster - update kalman track
32 // The seeding and tracking is repeated several times, in different seeding region.
33 // This approach enables to find the track which cannot be seeded in some region of TPC
34 // This can happen because of low momenta (track do not reach outer radius), or track is currently in the ded region between sectors, or the track is for the moment overlapped with other track (seed quality is poor) ...
36 // With this approach we reach almost 100 % efficiency also for high occupancy events.
37 // (If the seeding efficiency in a region is about 90 % than with logical or of several
38 // regions we will reach 100% (in theory - supposing independence)
40 // Repeating several seeding - tracking procedures some of the tracks can be find
43 // The procedures to remove multi find tacks are impremented:
44 // RemoveUsed2 - fast procedure n problem -
45 // Algorithm - Sorting tracks according quality
46 // remove tracks with some shared fraction
47 // Sharing in respect to all tacks
48 // Signing clusters in gold region
49 // FindSplitted - slower algorithm n^2
50 // Sort the tracks according quality
51 // Loop over pair of tracks
52 // If overlap with other track bigger than threshold - remove track
54 // FindCurling - Finds the pair of tracks which are curling
55 // - About 10% of tracks can be find with this procedure
56 // The combinatorial background is too big to be used in High
57 // multiplicity environment
58 // - n^2 problem - Slow procedure - currently it is disabled because of
61 // The number of splitted tracks can be reduced disabling the sharing of the cluster.
62 // tpcRecoParam-> SetClusterSharing(kFALSE);
63 // IT IS HIGHLY non recomended to use it in high flux enviroonment
64 // Even using this switch some tracks can be found more than once
65 // (because of multiple seeding and low quality tracks which will not cross full chamber)
68 // The tracker itself can be debugged - the information about tracks can be stored in several // phases of the reconstruction
69 // To enable storage of the TPC tracks in the ESD friend track
70 // use AliTPCReconstructor::SetStreamLevel(n);
72 // The debug level - different procedure produce tree for numerical debugging
73 // To enable them set AliTPCReconstructor::SetStreamLevel(n); where nis bigger 1
77 // Adding systematic errors to the covariance:
79 // The systematic errors due to the misalignment and miscalibration are added to the covariance matrix
80 // of the tracks (not to the clusters as they are dependent):
81 // The parameters form AliTPCRecoParam are used AliTPCRecoParam::GetSystematicError
82 // The systematic errors are expressed there in RMS - position (cm), angle (rad), curvature (1/cm)
83 // The default values are 0.
85 // The sytematic errors are added to the covariance matrix in following places:
87 // 1. During fisrt itteration - AliTPCtrackerMI::FillESD
88 // 2. Second iteration -
89 // 2.a ITS->TPC - AliTPCtrackerMI::ReadSeeds
90 // 2.b TPC->TRD - AliTPCtrackerMI::PropagateBack
91 // 3. Third iteration -
92 // 3.a TRD->TPC - AliTPCtrackerMI::ReadSeeds
93 // 3.b TPC->ITS - AliTPCtrackerMI::RefitInward
95 // There are several places in the code which can be numerically debuged
96 // This code is keeped in order to enable code development and to check the calibration implementtion
98 // 1. ErrParam stream - dump information about
100 // 2.a) cluster error estimate
101 // 3.a) cluster shape estimate
104 // Debug streamer levels:
106 //-------------------------------------------------------
111 #include "Riostream.h"
112 #include <TClonesArray.h>
114 #include <TObjArray.h>
117 #include "AliComplexCluster.h"
118 #include "AliESDEvent.h"
119 #include "AliESDtrack.h"
120 #include "AliESDVertex.h"
123 #include "AliHelix.h"
124 #include "AliRunLoader.h"
125 #include "AliTPCClustersRow.h"
126 #include "AliTPCParam.h"
127 #include "AliTPCReconstructor.h"
128 #include "AliTPCpolyTrack.h"
129 #include "AliTPCreco.h"
130 #include "AliTPCseed.h"
132 #include "AliTPCtrackerSector.h"
133 #include "AliTPCtrackerMI.h"
134 #include "TStopwatch.h"
135 #include "AliTPCReconstructor.h"
136 #include "AliAlignObj.h"
137 #include "AliTrackPointArray.h"
139 #include "AliTPCcalibDB.h"
140 #include "AliTPCTransform.h"
141 #include "AliTPCClusterParam.h"
145 ClassImp(AliTPCtrackerMI)
149 class AliTPCFastMath {
152 static Double_t FastAsin(Double_t x);
154 static Double_t fgFastAsin[20000]; //lookup table for fast asin computation
157 Double_t AliTPCFastMath::fgFastAsin[20000];
158 AliTPCFastMath gAliTPCFastMath; // needed to fill the LUT
160 AliTPCFastMath::AliTPCFastMath(){
162 // initialized lookup table;
163 for (Int_t i=0;i<10000;i++){
164 fgFastAsin[2*i] = TMath::ASin(i/10000.);
165 fgFastAsin[2*i+1] = (TMath::ASin((i+1)/10000.)-fgFastAsin[2*i]);
169 Double_t AliTPCFastMath::FastAsin(Double_t x){
171 // return asin using lookup table
173 Int_t index = int(x*10000);
174 return fgFastAsin[2*index]+(x*10000.-index)*fgFastAsin[2*index+1];
177 Int_t index = int(x*10000);
178 return -(fgFastAsin[2*index]+(x*10000.-index)*fgFastAsin[2*index+1]);
180 //__________________________________________________________________
181 AliTPCtrackerMI::AliTPCtrackerMI()
203 // default constructor
205 for (Int_t irow=0; irow<200; irow++){
212 //_____________________________________________________________________
216 Int_t AliTPCtrackerMI::UpdateTrack(AliTPCseed * track, Int_t accept){
218 //update track information using current cluster - track->fCurrentCluster
221 AliTPCclusterMI* c =track->GetCurrentCluster();
222 if (accept>0) track->SetCurrentClusterIndex1(track->GetCurrentClusterIndex1() | 0x8000); //sign not accepted clusters
224 UInt_t i = track->GetCurrentClusterIndex1();
226 Int_t sec=(i&0xff000000)>>24;
227 //Int_t row = (i&0x00ff0000)>>16;
228 track->SetRow((i&0x00ff0000)>>16);
229 track->SetSector(sec);
230 // Int_t index = i&0xFFFF;
231 if (sec>=fkParam->GetNInnerSector()) track->SetRow(track->GetRow()+fkParam->GetNRowLow());
232 track->SetClusterIndex2(track->GetRow(), i);
233 //track->fFirstPoint = row;
234 //if ( track->fLastPoint<row) track->fLastPoint =row;
235 // if (track->fRow<0 || track->fRow>160) {
236 // printf("problem\n");
238 if (track->GetFirstPoint()>track->GetRow())
239 track->SetFirstPoint(track->GetRow());
240 if (track->GetLastPoint()<track->GetRow())
241 track->SetLastPoint(track->GetRow());
244 track->SetClusterPointer(track->GetRow(),c);
247 Double_t angle2 = track->GetSnp()*track->GetSnp();
249 //SET NEW Track Point
251 if (angle2<1) //PH sometimes angle2 is very big. To be investigated...
253 angle2 = TMath::Sqrt(angle2/(1-angle2));
254 AliTPCTrackerPoint &point =*(track->GetTrackPoint(track->GetRow()));
256 point.SetSigmaY(c->GetSigmaY2()/track->GetCurrentSigmaY2());
257 point.SetSigmaZ(c->GetSigmaZ2()/track->GetCurrentSigmaZ2());
258 point.SetErrY(sqrt(track->GetErrorY2()));
259 point.SetErrZ(sqrt(track->GetErrorZ2()));
261 point.SetX(track->GetX());
262 point.SetY(track->GetY());
263 point.SetZ(track->GetZ());
264 point.SetAngleY(angle2);
265 point.SetAngleZ(track->GetTgl());
266 if (point.IsShared()){
267 track->SetErrorY2(track->GetErrorY2()*4);
268 track->SetErrorZ2(track->GetErrorZ2()*4);
272 Double_t chi2 = track->GetPredictedChi2(track->GetCurrentCluster());
274 // track->SetErrorY2(track->GetErrorY2()*1.3);
275 // track->SetErrorY2(track->GetErrorY2()+0.01);
276 // track->SetErrorZ2(track->GetErrorZ2()*1.3);
277 // track->SetErrorZ2(track->GetErrorZ2()+0.005);
279 if (accept>0) return 0;
280 if (track->GetNumberOfClusters()%20==0){
281 // if (track->fHelixIn){
282 // TClonesArray & larr = *(track->fHelixIn);
283 // Int_t ihelix = larr.GetEntriesFast();
284 // new(larr[ihelix]) AliHelix(*track) ;
287 track->SetNoCluster(0);
288 return track->Update(c,chi2,i);
293 Int_t AliTPCtrackerMI::AcceptCluster(AliTPCseed * seed, AliTPCclusterMI * cluster)
296 // decide according desired precision to accept given
297 // cluster for tracking
299 seed->GetProlongation(cluster->GetX(),yt,zt);
300 Double_t sy2=ErrY2(seed,cluster);
301 Double_t sz2=ErrZ2(seed,cluster);
303 Double_t sdistancey2 = sy2+seed->GetSigmaY2();
304 Double_t sdistancez2 = sz2+seed->GetSigmaZ2();
305 Double_t dy=seed->GetCurrentCluster()->GetY()-yt;
306 Double_t dz=seed->GetCurrentCluster()->GetZ()-zt;
307 Double_t rdistancey2 = (seed->GetCurrentCluster()->GetY()-yt)*
308 (seed->GetCurrentCluster()->GetY()-yt)/sdistancey2;
309 Double_t rdistancez2 = (seed->GetCurrentCluster()->GetZ()-zt)*
310 (seed->GetCurrentCluster()->GetZ()-zt)/sdistancez2;
312 Double_t rdistance2 = rdistancey2+rdistancez2;
315 if (AliTPCReconstructor::StreamLevel()>2 && seed->GetNumberOfClusters()>20) {
316 Float_t rmsy2 = seed->GetCurrentSigmaY2();
317 Float_t rmsz2 = seed->GetCurrentSigmaZ2();
318 Float_t rmsy2p30 = seed->GetCMeanSigmaY2p30();
319 Float_t rmsz2p30 = seed->GetCMeanSigmaZ2p30();
320 Float_t rmsy2p30R = seed->GetCMeanSigmaY2p30R();
321 Float_t rmsz2p30R = seed->GetCMeanSigmaZ2p30R();
322 AliExternalTrackParam param(*seed);
323 static TVectorD gcl(3),gtr(3);
325 param.GetXYZ(gcl.GetMatrixArray());
326 cluster->GetGlobalXYZ(gclf);
327 gcl[0]=gclf[0]; gcl[1]=gclf[1]; gcl[2]=gclf[2];
330 if (AliTPCReconstructor::StreamLevel()>2) {
331 (*fDebugStreamer)<<"ErrParam"<<
344 "rmsy2p30="<<rmsy2p30<<
345 "rmsz2p30="<<rmsz2p30<<
346 "rmsy2p30R="<<rmsy2p30R<<
347 "rmsz2p30R="<<rmsz2p30R<<
348 // normalize distance -
349 "rdisty="<<rdistancey2<<
350 "rdistz="<<rdistancez2<<
351 "rdist="<<rdistance2<< //
355 //return 0; // temporary
356 if (rdistance2>32) return 3;
359 if ((rdistancey2>9. || rdistancez2>9.) && cluster->GetType()==0)
360 return 2; //suspisiouce - will be changed
362 if ((rdistancey2>6.25 || rdistancez2>6.25) && cluster->GetType()>0)
363 // strict cut on overlaped cluster
364 return 2; //suspisiouce - will be changed
366 if ( (rdistancey2>1. || rdistancez2>6.25 )
367 && cluster->GetType()<0){
368 seed->SetNFoundable(seed->GetNFoundable()-1);
378 //_____________________________________________________________________________
379 AliTPCtrackerMI::AliTPCtrackerMI(const AliTPCParam *par):
381 fkNIS(par->GetNInnerSector()/2),
383 fkNOS(par->GetNOuterSector()/2),
400 //---------------------------------------------------------------------
401 // The main TPC tracker constructor
402 //---------------------------------------------------------------------
403 fInnerSec=new AliTPCtrackerSector[fkNIS];
404 fOuterSec=new AliTPCtrackerSector[fkNOS];
407 for (i=0; i<fkNIS; i++) fInnerSec[i].Setup(par,0);
408 for (i=0; i<fkNOS; i++) fOuterSec[i].Setup(par,1);
411 Int_t nrowlow = par->GetNRowLow();
412 Int_t nrowup = par->GetNRowUp();
415 for (i=0;i<nrowlow;i++){
416 fXRow[i] = par->GetPadRowRadiiLow(i);
417 fPadLength[i]= par->GetPadPitchLength(0,i);
418 fYMax[i] = fXRow[i]*TMath::Tan(0.5*par->GetInnerAngle());
422 for (i=0;i<nrowup;i++){
423 fXRow[i+nrowlow] = par->GetPadRowRadiiUp(i);
424 fPadLength[i+nrowlow] = par->GetPadPitchLength(60,i);
425 fYMax[i+nrowlow] = fXRow[i+nrowlow]*TMath::Tan(0.5*par->GetOuterAngle());
428 if (AliTPCReconstructor::StreamLevel()>0) {
429 fDebugStreamer = new TTreeSRedirector("TPCdebug.root");
432 //________________________________________________________________________
433 AliTPCtrackerMI::AliTPCtrackerMI(const AliTPCtrackerMI &t):
454 //------------------------------------
455 // dummy copy constructor
456 //------------------------------------------------------------------
459 AliTPCtrackerMI & AliTPCtrackerMI::operator=(const AliTPCtrackerMI& /*r*/)
461 //------------------------------
463 //--------------------------------------------------------------
466 //_____________________________________________________________________________
467 AliTPCtrackerMI::~AliTPCtrackerMI() {
468 //------------------------------------------------------------------
469 // TPC tracker destructor
470 //------------------------------------------------------------------
477 if (fDebugStreamer) delete fDebugStreamer;
481 void AliTPCtrackerMI::FillESD(const TObjArray* arr)
485 //fill esds using updated tracks
488 // write tracks to the event
489 // store index of the track
490 Int_t nseed=arr->GetEntriesFast();
491 //FindKinks(arr,fEvent);
492 for (Int_t i=0; i<nseed; i++) {
493 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
497 if (AliTPCReconstructor::StreamLevel()>1) {
498 (*fDebugStreamer)<<"Track0"<<
502 // pt->PropagateTo(fkParam->GetInnerRadiusLow());
503 if (pt->GetKinkIndex(0)<=0){ //don't propagate daughter tracks
504 pt->PropagateTo(fkParam->GetInnerRadiusLow());
507 if (( pt->GetPoints()[2]- pt->GetPoints()[0])>5 && pt->GetPoints()[3]>0.8){
509 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
510 iotrack.SetTPCPoints(pt->GetPoints());
511 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
512 iotrack.SetV0Indexes(pt->GetV0Indexes());
513 // iotrack.SetTPCpid(pt->fTPCr);
514 //iotrack.SetTPCindex(i);
515 fEvent->AddTrack(&iotrack);
519 if ( (pt->GetNumberOfClusters()>70)&& (Float_t(pt->GetNumberOfClusters())/Float_t(pt->GetNFoundable()))>0.55) {
521 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
522 iotrack.SetTPCPoints(pt->GetPoints());
523 //iotrack.SetTPCindex(i);
524 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
525 iotrack.SetV0Indexes(pt->GetV0Indexes());
526 // iotrack.SetTPCpid(pt->fTPCr);
527 fEvent->AddTrack(&iotrack);
531 // short tracks - maybe decays
533 if ( (pt->GetNumberOfClusters()>30) && (Float_t(pt->GetNumberOfClusters())/Float_t(pt->GetNFoundable()))>0.70) {
534 Int_t found,foundable,shared;
535 pt->GetClusterStatistic(0,60,found, foundable,shared,kFALSE);
536 if ( (found>20) && (pt->GetNShared()/float(pt->GetNumberOfClusters())<0.2)){
538 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
539 //iotrack.SetTPCindex(i);
540 iotrack.SetTPCPoints(pt->GetPoints());
541 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
542 iotrack.SetV0Indexes(pt->GetV0Indexes());
543 //iotrack.SetTPCpid(pt->fTPCr);
544 fEvent->AddTrack(&iotrack);
549 if ( (pt->GetNumberOfClusters()>20) && (Float_t(pt->GetNumberOfClusters())/Float_t(pt->GetNFoundable()))>0.8) {
550 Int_t found,foundable,shared;
551 pt->GetClusterStatistic(0,60,found, foundable,shared,kFALSE);
552 if (found<20) continue;
553 if (pt->GetNShared()/float(pt->GetNumberOfClusters())>0.2) continue;
556 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
557 iotrack.SetTPCPoints(pt->GetPoints());
558 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
559 iotrack.SetV0Indexes(pt->GetV0Indexes());
560 //iotrack.SetTPCpid(pt->fTPCr);
561 //iotrack.SetTPCindex(i);
562 fEvent->AddTrack(&iotrack);
565 // short tracks - secondaties
567 if ( (pt->GetNumberOfClusters()>30) ) {
568 Int_t found,foundable,shared;
569 pt->GetClusterStatistic(128,158,found, foundable,shared,kFALSE);
570 if ( (found>20) && (pt->GetNShared()/float(pt->GetNumberOfClusters())<0.2) &&float(found)/float(foundable)>0.8){
572 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
573 iotrack.SetTPCPoints(pt->GetPoints());
574 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
575 iotrack.SetV0Indexes(pt->GetV0Indexes());
576 //iotrack.SetTPCpid(pt->fTPCr);
577 //iotrack.SetTPCindex(i);
578 fEvent->AddTrack(&iotrack);
583 if ( (pt->GetNumberOfClusters()>15)) {
584 Int_t found,foundable,shared;
585 pt->GetClusterStatistic(138,158,found, foundable,shared,kFALSE);
586 if (found<15) continue;
587 if (foundable<=0) continue;
588 if (pt->GetNShared()/float(pt->GetNumberOfClusters())>0.2) continue;
589 if (float(found)/float(foundable)<0.8) continue;
592 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
593 iotrack.SetTPCPoints(pt->GetPoints());
594 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
595 iotrack.SetV0Indexes(pt->GetV0Indexes());
596 // iotrack.SetTPCpid(pt->fTPCr);
597 //iotrack.SetTPCindex(i);
598 fEvent->AddTrack(&iotrack);
602 // >> account for suppressed tracks in the kink indices (RS)
603 int nESDtracks = fEvent->GetNumberOfTracks();
604 for (int it=nESDtracks;it--;) {
605 AliESDtrack* esdTr = fEvent->GetTrack(it);
606 if (!esdTr || !esdTr->GetKinkIndex(0)) continue;
607 for (int ik=0;ik<3;ik++) {
609 if (!(knkId=esdTr->GetKinkIndex(ik))) break; // no more kinks for this track
610 AliESDkink* kink = fEvent->GetKink(TMath::Abs(knkId)-1);
612 AliError(Form("ESDTrack%d refers to non-existing kink %d",it,TMath::Abs(knkId)-1));
615 kink->SetIndex(it, knkId<0 ? 0:1); // update track index of the kink: mother at 0, daughter at 1
618 // << account for suppressed tracks in the kink indices (RS)
619 AliInfo(Form("Number of filled ESDs-\t%d\n",fEvent->GetNumberOfTracks()));
627 Double_t AliTPCtrackerMI::ErrY2(AliTPCseed* seed, const AliTPCclusterMI * cl){
630 // Use calibrated cluster error from OCDB
632 AliTPCClusterParam * clparam = AliTPCcalibDB::Instance()->GetClusterParam();
634 Float_t z = TMath::Abs(fkParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
635 Int_t ctype = cl->GetType();
636 Int_t type = (cl->GetRow()<63) ? 0: (cl->GetRow()>126) ? 1:2;
637 Double_t angle = seed->GetSnp()*seed->GetSnp();
638 angle = TMath::Sqrt(TMath::Abs(angle/(1.-angle)));
639 Double_t erry2 = clparam->GetError0Par(0,type, z,angle);
641 erry2+=0.5; // edge cluster
644 seed->SetErrorY2(erry2);
648 //calculate look-up table at the beginning
649 // static Bool_t ginit = kFALSE;
650 // static Float_t gnoise1,gnoise2,gnoise3;
651 // static Float_t ggg1[10000];
652 // static Float_t ggg2[10000];
653 // static Float_t ggg3[10000];
654 // static Float_t glandau1[10000];
655 // static Float_t glandau2[10000];
656 // static Float_t glandau3[10000];
658 // static Float_t gcor01[500];
659 // static Float_t gcor02[500];
660 // static Float_t gcorp[500];
664 // if (ginit==kFALSE){
665 // for (Int_t i=1;i<500;i++){
666 // Float_t rsigma = float(i)/100.;
667 // gcor02[i] = TMath::Max(0.78 +TMath::Exp(7.4*(rsigma-1.2)),0.6);
668 // gcor01[i] = TMath::Max(0.72 +TMath::Exp(3.36*(rsigma-1.2)),0.6);
669 // gcorp[i] = TMath::Max(TMath::Power((rsigma+0.5),1.5),1.2);
673 // for (Int_t i=3;i<10000;i++){
677 // Float_t amp = float(i);
678 // Float_t padlength =0.75;
679 // gnoise1 = 0.0004/padlength;
680 // Float_t nel = 0.268*amp;
681 // Float_t nprim = 0.155*amp;
682 // ggg1[i] = fkParam->GetDiffT()*fkParam->GetDiffT()*(2+0.001*nel/(padlength*padlength))/nel;
683 // glandau1[i] = (2.+0.12*nprim)*0.5* (2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
684 // if (glandau1[i]>1) glandau1[i]=1;
685 // glandau1[i]*=padlength*padlength/12.;
689 // gnoise2 = 0.0004/padlength;
691 // nprim = 0.133*amp;
692 // ggg2[i] = fkParam->GetDiffT()*fkParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
693 // glandau2[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
694 // if (glandau2[i]>1) glandau2[i]=1;
695 // glandau2[i]*=padlength*padlength/12.;
700 // gnoise3 = 0.0004/padlength;
702 // nprim = 0.133*amp;
703 // ggg3[i] = fkParam->GetDiffT()*fkParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
704 // glandau3[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
705 // if (glandau3[i]>1) glandau3[i]=1;
706 // glandau3[i]*=padlength*padlength/12.;
714 // Int_t amp = int(TMath::Abs(cl->GetQ()));
716 // seed->SetErrorY2(1.);
720 // Float_t z = TMath::Abs(fkParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
721 // Int_t ctype = cl->GetType();
722 // Float_t padlength= GetPadPitchLength(seed->GetRow());
723 // Double_t angle2 = seed->GetSnp()*seed->GetSnp();
724 // angle2 = angle2/(1-angle2);
726 // //cluster "quality"
727 // Int_t rsigmay = int(100.*cl->GetSigmaY2()/(seed->GetCurrentSigmaY2()));
730 // if (fSectors==fInnerSec){
731 // snoise2 = gnoise1;
732 // res = ggg1[amp]*z+glandau1[amp]*angle2;
733 // if (ctype==0) res *= gcor01[rsigmay];
736 // res*= gcorp[rsigmay];
740 // if (padlength<1.1){
741 // snoise2 = gnoise2;
742 // res = ggg2[amp]*z+glandau2[amp]*angle2;
743 // if (ctype==0) res *= gcor02[rsigmay];
746 // res*= gcorp[rsigmay];
750 // snoise2 = gnoise3;
751 // res = ggg3[amp]*z+glandau3[amp]*angle2;
752 // if (ctype==0) res *= gcor02[rsigmay];
755 // res*= gcorp[rsigmay];
762 // res*=2.4; // overestimate error 2 times
766 // if (res<2*snoise2)
769 // seed->SetErrorY2(res);
777 Double_t AliTPCtrackerMI::ErrZ2(AliTPCseed* seed, const AliTPCclusterMI * cl){
780 // Use calibrated cluster error from OCDB
782 AliTPCClusterParam * clparam = AliTPCcalibDB::Instance()->GetClusterParam();
784 Float_t z = TMath::Abs(fkParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
785 Int_t ctype = cl->GetType();
786 Int_t type = (cl->GetRow()<63) ? 0: (cl->GetRow()>126) ? 1:2;
788 Double_t angle2 = seed->GetSnp()*seed->GetSnp();
789 angle2 = seed->GetTgl()*seed->GetTgl()*(1+angle2/(1-angle2));
790 Double_t angle = TMath::Sqrt(TMath::Abs(angle2));
791 Double_t errz2 = clparam->GetError0Par(1,type, z,angle);
793 errz2+=0.5; // edge cluster
796 seed->SetErrorZ2(errz2);
802 // //seed->SetErrorY2(0.1);
804 // //calculate look-up table at the beginning
805 // static Bool_t ginit = kFALSE;
806 // static Float_t gnoise1,gnoise2,gnoise3;
807 // static Float_t ggg1[10000];
808 // static Float_t ggg2[10000];
809 // static Float_t ggg3[10000];
810 // static Float_t glandau1[10000];
811 // static Float_t glandau2[10000];
812 // static Float_t glandau3[10000];
814 // static Float_t gcor01[1000];
815 // static Float_t gcor02[1000];
816 // static Float_t gcorp[1000];
820 // if (ginit==kFALSE){
821 // for (Int_t i=1;i<1000;i++){
822 // Float_t rsigma = float(i)/100.;
823 // gcor02[i] = TMath::Max(0.81 +TMath::Exp(6.8*(rsigma-1.2)),0.6);
824 // gcor01[i] = TMath::Max(0.72 +TMath::Exp(2.04*(rsigma-1.2)),0.6);
825 // gcorp[i] = TMath::Max(TMath::Power((rsigma+0.5),1.5),1.2);
829 // for (Int_t i=3;i<10000;i++){
833 // Float_t amp = float(i);
834 // Float_t padlength =0.75;
835 // gnoise1 = 0.0004/padlength;
836 // Float_t nel = 0.268*amp;
837 // Float_t nprim = 0.155*amp;
838 // ggg1[i] = fkParam->GetDiffT()*fkParam->GetDiffT()*(2+0.001*nel/(padlength*padlength))/nel;
839 // glandau1[i] = (2.+0.12*nprim)*0.5* (2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
840 // if (glandau1[i]>1) glandau1[i]=1;
841 // glandau1[i]*=padlength*padlength/12.;
845 // gnoise2 = 0.0004/padlength;
847 // nprim = 0.133*amp;
848 // ggg2[i] = fkParam->GetDiffT()*fkParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
849 // glandau2[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
850 // if (glandau2[i]>1) glandau2[i]=1;
851 // glandau2[i]*=padlength*padlength/12.;
856 // gnoise3 = 0.0004/padlength;
858 // nprim = 0.133*amp;
859 // ggg3[i] = fkParam->GetDiffT()*fkParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
860 // glandau3[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
861 // if (glandau3[i]>1) glandau3[i]=1;
862 // glandau3[i]*=padlength*padlength/12.;
870 // Int_t amp = int(TMath::Abs(cl->GetQ()));
872 // seed->SetErrorY2(1.);
876 // Float_t z = TMath::Abs(fkParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
877 // Int_t ctype = cl->GetType();
878 // Float_t padlength= GetPadPitchLength(seed->GetRow());
880 // Double_t angle2 = seed->GetSnp()*seed->GetSnp();
881 // // if (angle2<0.6) angle2 = 0.6;
882 // angle2 = seed->GetTgl()*seed->GetTgl()*(1+angle2/(1-angle2));
884 // //cluster "quality"
885 // Int_t rsigmaz = int(100.*cl->GetSigmaZ2()/(seed->GetCurrentSigmaZ2()));
888 // if (fSectors==fInnerSec){
889 // snoise2 = gnoise1;
890 // res = ggg1[amp]*z+glandau1[amp]*angle2;
891 // if (ctype==0) res *= gcor01[rsigmaz];
894 // res*= gcorp[rsigmaz];
898 // if (padlength<1.1){
899 // snoise2 = gnoise2;
900 // res = ggg2[amp]*z+glandau2[amp]*angle2;
901 // if (ctype==0) res *= gcor02[rsigmaz];
904 // res*= gcorp[rsigmaz];
908 // snoise2 = gnoise3;
909 // res = ggg3[amp]*z+glandau3[amp]*angle2;
910 // if (ctype==0) res *= gcor02[rsigmaz];
913 // res*= gcorp[rsigmaz];
922 // if ((ctype<0) &&<70){
927 // if (res<2*snoise2)
929 // if (res>3) res =3;
930 // seed->SetErrorZ2(res);
938 void AliTPCtrackerMI::RotateToLocal(AliTPCseed *seed)
940 //rotate to track "local coordinata
941 Float_t x = seed->GetX();
942 Float_t y = seed->GetY();
943 Float_t ymax = x*TMath::Tan(0.5*fSectors->GetAlpha());
946 seed->SetRelativeSector((seed->GetRelativeSector()+1) % fN);
947 if (!seed->Rotate(fSectors->GetAlpha()))
949 } else if (y <-ymax) {
950 seed->SetRelativeSector((seed->GetRelativeSector()-1+fN) % fN);
951 if (!seed->Rotate(-fSectors->GetAlpha()))
959 //_____________________________________________________________________________
960 Double_t AliTPCtrackerMI::F1old(Double_t x1,Double_t y1,
961 Double_t x2,Double_t y2,
962 Double_t x3,Double_t y3) const
964 //-----------------------------------------------------------------
965 // Initial approximation of the track curvature
966 //-----------------------------------------------------------------
967 Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
968 Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
969 (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
970 Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
971 (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
973 Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
974 if ( xr*xr+yr*yr<=0.00000000000001) return 100;
975 return -xr*yr/sqrt(xr*xr+yr*yr);
980 //_____________________________________________________________________________
981 Double_t AliTPCtrackerMI::F1(Double_t x1,Double_t y1,
982 Double_t x2,Double_t y2,
983 Double_t x3,Double_t y3) const
985 //-----------------------------------------------------------------
986 // Initial approximation of the track curvature
987 //-----------------------------------------------------------------
993 Double_t det = x3*y2-x2*y3;
994 if (TMath::Abs(det)<1e-10){
998 Double_t u = 0.5* (x2*(x2-x3)+y2*(y2-y3))/det;
999 Double_t x0 = x3*0.5-y3*u;
1000 Double_t y0 = y3*0.5+x3*u;
1001 Double_t c2 = 1/TMath::Sqrt(x0*x0+y0*y0);
1007 Double_t AliTPCtrackerMI::F2(Double_t x1,Double_t y1,
1008 Double_t x2,Double_t y2,
1009 Double_t x3,Double_t y3) const
1011 //-----------------------------------------------------------------
1012 // Initial approximation of the track curvature
1013 //-----------------------------------------------------------------
1019 Double_t det = x3*y2-x2*y3;
1020 if (TMath::Abs(det)<1e-10) {
1024 Double_t u = 0.5* (x2*(x2-x3)+y2*(y2-y3))/det;
1025 Double_t x0 = x3*0.5-y3*u;
1026 Double_t y0 = y3*0.5+x3*u;
1027 Double_t c2 = 1/TMath::Sqrt(x0*x0+y0*y0);
1036 //_____________________________________________________________________________
1037 Double_t AliTPCtrackerMI::F2old(Double_t x1,Double_t y1,
1038 Double_t x2,Double_t y2,
1039 Double_t x3,Double_t y3) const
1041 //-----------------------------------------------------------------
1042 // Initial approximation of the track curvature times center of curvature
1043 //-----------------------------------------------------------------
1044 Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
1045 Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
1046 (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
1047 Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
1048 (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
1050 Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
1052 return -a/(d*y1-b)*xr/sqrt(xr*xr+yr*yr);
1055 //_____________________________________________________________________________
1056 Double_t AliTPCtrackerMI::F3(Double_t x1,Double_t y1,
1057 Double_t x2,Double_t y2,
1058 Double_t z1,Double_t z2) const
1060 //-----------------------------------------------------------------
1061 // Initial approximation of the tangent of the track dip angle
1062 //-----------------------------------------------------------------
1063 return (z1 - z2)/sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
1067 Double_t AliTPCtrackerMI::F3n(Double_t x1,Double_t y1,
1068 Double_t x2,Double_t y2,
1069 Double_t z1,Double_t z2, Double_t c) const
1071 //-----------------------------------------------------------------
1072 // Initial approximation of the tangent of the track dip angle
1073 //-----------------------------------------------------------------
1077 //angle1 = (z1-z2)*c/(TMath::ASin(c*x1-ni)-TMath::ASin(c*x2-ni));
1079 Double_t d = TMath::Sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
1080 if (TMath::Abs(d*c*0.5)>1) return 0;
1081 // Double_t angle2 = TMath::ASin(d*c*0.5);
1082 // Double_t angle2 = AliTPCFastMath::FastAsin(d*c*0.5);
1083 Double_t angle2 = (d*c*0.5>0.1)? TMath::ASin(d*c*0.5): AliTPCFastMath::FastAsin(d*c*0.5);
1085 angle2 = (z1-z2)*c/(angle2*2.);
1089 Bool_t AliTPCtrackerMI::GetProlongation(Double_t x1, Double_t x2, Double_t x[5], Double_t &y, Double_t &z) const
1090 {//-----------------------------------------------------------------
1091 // This function find proloncation of a track to a reference plane x=x2.
1092 //-----------------------------------------------------------------
1096 if (TMath::Abs(x[4]*x1 - x[2]) >= 0.999) {
1100 Double_t c1=x[4]*x1 - x[2], r1=TMath::Sqrt((1.-c1)*(1.+c1));
1101 Double_t c2=x[4]*x2 - x[2], r2=TMath::Sqrt((1.-c2)*(1.+c2));
1105 Double_t dy = dx*(c1+c2)/(r1+r2);
1108 Double_t delta = x[4]*dx*(c1+c2)/(c1*r2 + c2*r1);
1110 if (TMath::Abs(delta)>0.01){
1111 dz = x[3]*TMath::ASin(delta)/x[4];
1113 dz = x[3]*AliTPCFastMath::FastAsin(delta)/x[4];
1116 //dz = x[3]*AliTPCFastMath::FastAsin(delta)/x[4];
1124 Int_t AliTPCtrackerMI::LoadClusters (TTree *const tree)
1129 return LoadClusters();
1133 Int_t AliTPCtrackerMI::LoadClusters(const TObjArray *arr)
1136 // load clusters to the memory
1137 AliTPCClustersRow *clrow = new AliTPCClustersRow("AliTPCclusterMI");
1138 Int_t lower = arr->LowerBound();
1139 Int_t entries = arr->GetEntriesFast();
1141 for (Int_t i=lower; i<entries; i++) {
1142 clrow = (AliTPCClustersRow*) arr->At(i);
1143 if(!clrow) continue;
1144 if(!clrow->GetArray()) continue;
1148 fkParam->AdjustSectorRow(clrow->GetID(),sec,row);
1150 for (Int_t icl=0; icl<clrow->GetArray()->GetEntriesFast(); icl++){
1151 Transform((AliTPCclusterMI*)(clrow->GetArray()->At(icl)));
1154 if (clrow->GetArray()->GetEntriesFast()<=0) continue;
1155 AliTPCtrackerRow * tpcrow=0;
1158 tpcrow = &(fInnerSec[sec%fkNIS][row]);
1162 tpcrow = &(fOuterSec[(sec-fkNIS*2)%fkNOS][row]);
1163 left = (sec-fkNIS*2)/fkNOS;
1166 tpcrow->SetN1(clrow->GetArray()->GetEntriesFast());
1167 tpcrow->SetClusters1(new AliTPCclusterMI[tpcrow->GetN1()]);
1168 for (Int_t j=0;j<tpcrow->GetN1();++j)
1169 tpcrow->SetCluster1(j, *(AliTPCclusterMI*)(clrow->GetArray()->At(j)));
1172 tpcrow->SetN2(clrow->GetArray()->GetEntriesFast());
1173 tpcrow->SetClusters2(new AliTPCclusterMI[tpcrow->GetN2()]);
1174 for (Int_t j=0;j<tpcrow->GetN2();++j)
1175 tpcrow->SetCluster2(j,*(AliTPCclusterMI*)(clrow->GetArray()->At(j)));
1177 clrow->GetArray()->Clear("C");
1186 Int_t AliTPCtrackerMI::LoadClusters(const TClonesArray *arr)
1189 // load clusters to the memory from one
1192 AliTPCclusterMI *clust=0;
1193 Int_t count[72][96] = { {0} , {0} };
1195 // loop over clusters
1196 for (Int_t icl=0; icl<arr->GetEntriesFast(); icl++) {
1197 clust = (AliTPCclusterMI*)arr->At(icl);
1198 if(!clust) continue;
1199 //printf("cluster: det %d, row %d \n", clust->GetDetector(),clust->GetRow());
1201 // transform clusters
1204 // count clusters per pad row
1205 count[clust->GetDetector()][clust->GetRow()]++;
1208 // insert clusters to sectors
1209 for (Int_t icl=0; icl<arr->GetEntriesFast(); icl++) {
1210 clust = (AliTPCclusterMI*)arr->At(icl);
1211 if(!clust) continue;
1213 Int_t sec = clust->GetDetector();
1214 Int_t row = clust->GetRow();
1216 // filter overlapping pad rows needed by HLT
1217 if(sec<fkNIS*2) { //IROCs
1218 if(row == 30) continue;
1221 if(row == 27 || row == 76) continue;
1227 fInnerSec[sec%fkNIS].InsertCluster(clust, count[sec][row], fkParam);
1230 left = (sec-fkNIS*2)/fkNOS;
1231 fOuterSec[(sec-fkNIS*2)%fkNOS].InsertCluster(clust, count[sec][row], fkParam);
1235 // Load functions must be called behind LoadCluster(TClonesArray*)
1237 //LoadOuterSectors();
1238 //LoadInnerSectors();
1244 Int_t AliTPCtrackerMI::LoadClusters()
1247 // load clusters to the memory
1248 AliTPCClustersRow *clrow= new AliTPCClustersRow("AliTPCclusterMI");
1250 // TTree * tree = fClustersArray.GetTree();
1252 TTree * tree = fInput;
1253 TBranch * br = tree->GetBranch("Segment");
1254 br->SetAddress(&clrow);
1256 Int_t j=Int_t(tree->GetEntries());
1257 for (Int_t i=0; i<j; i++) {
1261 fkParam->AdjustSectorRow(clrow->GetID(),sec,row);
1262 for (Int_t icl=0; icl<clrow->GetArray()->GetEntriesFast(); icl++){
1263 Transform((AliTPCclusterMI*)(clrow->GetArray()->At(icl)));
1266 AliTPCtrackerRow * tpcrow=0;
1269 tpcrow = &(fInnerSec[sec%fkNIS][row]);
1273 tpcrow = &(fOuterSec[(sec-fkNIS*2)%fkNOS][row]);
1274 left = (sec-fkNIS*2)/fkNOS;
1277 tpcrow->SetN1(clrow->GetArray()->GetEntriesFast());
1278 tpcrow->SetClusters1(new AliTPCclusterMI[tpcrow->GetN1()]);
1279 for (Int_t k=0;k<tpcrow->GetN1();++k)
1280 tpcrow->SetCluster1(k, *(AliTPCclusterMI*)(clrow->GetArray()->At(k)));
1283 tpcrow->SetN2(clrow->GetArray()->GetEntriesFast());
1284 tpcrow->SetClusters2(new AliTPCclusterMI[tpcrow->GetN2()]);
1285 for (Int_t k=0;k<tpcrow->GetN2();++k)
1286 tpcrow->SetCluster2(k,*(AliTPCclusterMI*)(clrow->GetArray()->At(k)));
1297 void AliTPCtrackerMI::UnloadClusters()
1300 // unload clusters from the memory
1302 Int_t nrows = fOuterSec->GetNRows();
1303 for (Int_t sec = 0;sec<fkNOS;sec++)
1304 for (Int_t row = 0;row<nrows;row++){
1305 AliTPCtrackerRow* tpcrow = &(fOuterSec[sec%fkNOS][row]);
1307 // if (tpcrow->fClusters1) delete []tpcrow->fClusters1;
1308 // if (tpcrow->fClusters2) delete []tpcrow->fClusters2;
1310 tpcrow->ResetClusters();
1313 nrows = fInnerSec->GetNRows();
1314 for (Int_t sec = 0;sec<fkNIS;sec++)
1315 for (Int_t row = 0;row<nrows;row++){
1316 AliTPCtrackerRow* tpcrow = &(fInnerSec[sec%fkNIS][row]);
1318 // if (tpcrow->fClusters1) delete []tpcrow->fClusters1;
1319 //if (tpcrow->fClusters2) delete []tpcrow->fClusters2;
1321 tpcrow->ResetClusters();
1327 void AliTPCtrackerMI::FillClusterArray(TObjArray* array) const{
1329 // Filling cluster to the array - For visualization purposes
1332 nrows = fOuterSec->GetNRows();
1333 for (Int_t sec = 0;sec<fkNOS;sec++)
1334 for (Int_t row = 0;row<nrows;row++){
1335 AliTPCtrackerRow* tpcrow = &(fOuterSec[sec%fkNOS][row]);
1336 if (!tpcrow) continue;
1337 for (Int_t icl = 0;icl<tpcrow->GetN();icl++){
1338 array->AddLast((TObject*)((*tpcrow)[icl]));
1341 nrows = fInnerSec->GetNRows();
1342 for (Int_t sec = 0;sec<fkNIS;sec++)
1343 for (Int_t row = 0;row<nrows;row++){
1344 AliTPCtrackerRow* tpcrow = &(fInnerSec[sec%fkNIS][row]);
1345 if (!tpcrow) continue;
1346 for (Int_t icl = 0;icl<tpcrow->GetN();icl++){
1347 array->AddLast((TObject*)(*tpcrow)[icl]);
1353 void AliTPCtrackerMI::Transform(AliTPCclusterMI * cluster){
1357 AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ;
1359 AliFatal("Tranformations not in calibDB");
1362 transform->SetCurrentRecoParam((AliTPCRecoParam*)AliTPCReconstructor::GetRecoParam());
1363 Double_t x[3]={cluster->GetRow(),cluster->GetPad(),cluster->GetTimeBin()};
1364 Int_t i[1]={cluster->GetDetector()};
1365 transform->Transform(x,i,0,1);
1366 // if (cluster->GetDetector()%36>17){
1371 // in debug mode check the transformation
1373 if (AliTPCReconstructor::StreamLevel()>2) {
1375 cluster->GetGlobalXYZ(gx);
1376 Int_t event = (fEvent==NULL)? 0: fEvent->GetEventNumberInFile();
1377 TTreeSRedirector &cstream = *fDebugStreamer;
1378 cstream<<"Transform"<<
1389 cluster->SetX(x[0]);
1390 cluster->SetY(x[1]);
1391 cluster->SetZ(x[2]);
1397 //if (!fkParam->IsGeoRead()) fkParam->ReadGeoMatrices();
1398 TGeoHMatrix *mat = fkParam->GetClusterMatrix(cluster->GetDetector());
1400 Double_t pos[3]= {cluster->GetX(),cluster->GetY(),cluster->GetZ()};
1401 Double_t posC[3]={cluster->GetX(),cluster->GetY(),cluster->GetZ()};
1402 if (mat) mat->LocalToMaster(pos,posC);
1404 // chack Loading of Geo matrices from GeoManager - TEMPORARY FIX
1406 cluster->SetX(posC[0]);
1407 cluster->SetY(posC[1]);
1408 cluster->SetZ(posC[2]);
1411 //_____________________________________________________________________________
1412 Int_t AliTPCtrackerMI::LoadOuterSectors() {
1413 //-----------------------------------------------------------------
1414 // This function fills outer TPC sectors with clusters.
1415 //-----------------------------------------------------------------
1416 Int_t nrows = fOuterSec->GetNRows();
1418 for (Int_t sec = 0;sec<fkNOS;sec++)
1419 for (Int_t row = 0;row<nrows;row++){
1420 AliTPCtrackerRow* tpcrow = &(fOuterSec[sec%fkNOS][row]);
1421 Int_t sec2 = sec+2*fkNIS;
1423 Int_t ncl = tpcrow->GetN1();
1425 AliTPCclusterMI *c= (tpcrow->GetCluster1(ncl));
1426 index=(((sec2<<8)+row)<<16)+ncl;
1427 tpcrow->InsertCluster(c,index);
1430 ncl = tpcrow->GetN2();
1432 AliTPCclusterMI *c= (tpcrow->GetCluster2(ncl));
1433 index=((((sec2+fkNOS)<<8)+row)<<16)+ncl;
1434 tpcrow->InsertCluster(c,index);
1437 // write indexes for fast acces
1439 for (Int_t i=0;i<510;i++)
1440 tpcrow->SetFastCluster(i,-1);
1441 for (Int_t i=0;i<tpcrow->GetN();i++){
1442 Int_t zi = Int_t((*tpcrow)[i]->GetZ()+255.);
1443 tpcrow->SetFastCluster(zi,i); // write index
1446 for (Int_t i=0;i<510;i++){
1447 if (tpcrow->GetFastCluster(i)<0)
1448 tpcrow->SetFastCluster(i,last);
1450 last = tpcrow->GetFastCluster(i);
1459 //_____________________________________________________________________________
1460 Int_t AliTPCtrackerMI::LoadInnerSectors() {
1461 //-----------------------------------------------------------------
1462 // This function fills inner TPC sectors with clusters.
1463 //-----------------------------------------------------------------
1464 Int_t nrows = fInnerSec->GetNRows();
1466 for (Int_t sec = 0;sec<fkNIS;sec++)
1467 for (Int_t row = 0;row<nrows;row++){
1468 AliTPCtrackerRow* tpcrow = &(fInnerSec[sec%fkNIS][row]);
1471 Int_t ncl = tpcrow->GetN1();
1473 AliTPCclusterMI *c= (tpcrow->GetCluster1(ncl));
1474 index=(((sec<<8)+row)<<16)+ncl;
1475 tpcrow->InsertCluster(c,index);
1478 ncl = tpcrow->GetN2();
1480 AliTPCclusterMI *c= (tpcrow->GetCluster2(ncl));
1481 index=((((sec+fkNIS)<<8)+row)<<16)+ncl;
1482 tpcrow->InsertCluster(c,index);
1485 // write indexes for fast acces
1487 for (Int_t i=0;i<510;i++)
1488 tpcrow->SetFastCluster(i,-1);
1489 for (Int_t i=0;i<tpcrow->GetN();i++){
1490 Int_t zi = Int_t((*tpcrow)[i]->GetZ()+255.);
1491 tpcrow->SetFastCluster(zi,i); // write index
1494 for (Int_t i=0;i<510;i++){
1495 if (tpcrow->GetFastCluster(i)<0)
1496 tpcrow->SetFastCluster(i,last);
1498 last = tpcrow->GetFastCluster(i);
1510 //_________________________________________________________________________
1511 AliTPCclusterMI *AliTPCtrackerMI::GetClusterMI(Int_t index) const {
1512 //--------------------------------------------------------------------
1513 // Return pointer to a given cluster
1514 //--------------------------------------------------------------------
1515 if (index<0) return 0; // no cluster
1516 Int_t sec=(index&0xff000000)>>24;
1517 Int_t row=(index&0x00ff0000)>>16;
1518 Int_t ncl=(index&0x00007fff)>>00;
1520 const AliTPCtrackerRow * tpcrow=0;
1521 AliTPCclusterMI * clrow =0;
1523 if (sec<0 || sec>=fkNIS*4) {
1524 AliWarning(Form("Wrong sector %d",sec));
1529 AliTPCtrackerSector& tracksec = fInnerSec[sec%fkNIS];
1530 if (tracksec.GetNRows()<=row) return 0;
1531 tpcrow = &(tracksec[row]);
1532 if (tpcrow==0) return 0;
1535 if (tpcrow->GetN1()<=ncl) return 0;
1536 clrow = tpcrow->GetClusters1();
1539 if (tpcrow->GetN2()<=ncl) return 0;
1540 clrow = tpcrow->GetClusters2();
1544 AliTPCtrackerSector& tracksec = fOuterSec[(sec-fkNIS*2)%fkNOS];
1545 if (tracksec.GetNRows()<=row) return 0;
1546 tpcrow = &(tracksec[row]);
1547 if (tpcrow==0) return 0;
1549 if (sec-2*fkNIS<fkNOS) {
1550 if (tpcrow->GetN1()<=ncl) return 0;
1551 clrow = tpcrow->GetClusters1();
1554 if (tpcrow->GetN2()<=ncl) return 0;
1555 clrow = tpcrow->GetClusters2();
1559 return &(clrow[ncl]);
1565 Int_t AliTPCtrackerMI::FollowToNext(AliTPCseed& t, Int_t nr) {
1566 //-----------------------------------------------------------------
1567 // This function tries to find a track prolongation to next pad row
1568 //-----------------------------------------------------------------
1570 Double_t x= GetXrow(nr), ymax=GetMaxY(nr);
1573 AliTPCclusterMI *cl=0;
1574 Int_t tpcindex= t.GetClusterIndex2(nr);
1576 // update current shape info every 5 pad-row
1577 // if ( (nr%5==0) || t.GetNumberOfClusters()<2 || (t.fCurrentSigmaY2<0.0001) ){
1581 if (fIteration>0 && tpcindex>=-1){ //if we have already clusters
1583 if (tpcindex==-1) return 0; //track in dead zone
1585 cl = t.GetClusterPointer(nr);
1586 if ( (cl==0) ) cl = GetClusterMI(tpcindex);
1587 t.SetCurrentClusterIndex1(tpcindex);
1590 Int_t relativesector = ((tpcindex&0xff000000)>>24)%18; // if previously accepted cluster in different sector
1591 Float_t angle = relativesector*fSectors->GetAlpha()+fSectors->GetAlphaShift();
1593 if (angle<-TMath::Pi()) angle += 2*TMath::Pi();
1594 if (angle>=TMath::Pi()) angle -= 2*TMath::Pi();
1596 if (TMath::Abs(angle-t.GetAlpha())>0.001){
1597 Double_t rotation = angle-t.GetAlpha();
1598 t.SetRelativeSector(relativesector);
1599 if (!t.Rotate(rotation)) return 0;
1601 if (!t.PropagateTo(x)) return 0;
1603 t.SetCurrentCluster(cl);
1605 Int_t accept = AcceptCluster(&t,t.GetCurrentCluster());
1606 if ((tpcindex&0x8000)==0) accept =0;
1608 //if founded cluster is acceptible
1609 if (cl->IsUsed(11)) { // id cluster is shared inrease uncertainty
1610 t.SetErrorY2(t.GetErrorY2()+0.03);
1611 t.SetErrorZ2(t.GetErrorZ2()+0.03);
1612 t.SetErrorY2(t.GetErrorY2()*3);
1613 t.SetErrorZ2(t.GetErrorZ2()*3);
1615 t.SetNFoundable(t.GetNFoundable()+1);
1616 UpdateTrack(&t,accept);
1621 if (TMath::Abs(t.GetSnp())>AliTPCReconstructor::GetMaxSnpTracker()) return 0; // cut on angle
1622 if (fIteration>1 && IsFindable(t)){
1623 // not look for new cluster during refitting
1624 t.SetNFoundable(t.GetNFoundable()+1);
1629 // if (TMath::Abs(t.GetSnp())>0.95 || TMath::Abs(x*t.GetC()-t.GetEta())>0.95) return 0;// patch 28 fev 06
1630 Double_t y=t.GetYat(x);
1631 if (TMath::Abs(y)>ymax){
1633 t.SetRelativeSector((t.GetRelativeSector()+1) % fN);
1634 if (!t.Rotate(fSectors->GetAlpha()))
1636 } else if (y <-ymax) {
1637 t.SetRelativeSector((t.GetRelativeSector()-1+fN) % fN);
1638 if (!t.Rotate(-fSectors->GetAlpha()))
1644 if (!t.PropagateTo(x)) {
1645 if (fIteration==0) t.SetRemoval(10);
1649 Double_t z=t.GetZ();
1652 if (!IsActive(t.GetRelativeSector(),nr)) {
1654 t.SetClusterIndex2(nr,-1);
1657 //AliInfo(Form("A - Sector%d phi %f - alpha %f", t.fRelativeSector,y/x, t.GetAlpha()));
1658 Bool_t isActive = IsActive(t.GetRelativeSector(),nr);
1659 Bool_t isActive2 = (nr>=fInnerSec->GetNRows()) ? fOuterSec[t.GetRelativeSector()][nr-fInnerSec->GetNRows()].GetN()>0:fInnerSec[t.GetRelativeSector()][nr].GetN()>0;
1661 if (!isActive || !isActive2) return 0;
1663 const AliTPCtrackerRow &krow=GetRow(t.GetRelativeSector(),nr);
1664 if ( (t.GetSigmaY2()<0) || t.GetSigmaZ2()<0) return 0;
1666 Double_t roadz = 1.;
1668 if (TMath::Abs(TMath::Abs(y)-ymax)<krow.GetDeadZone()){
1670 t.SetClusterIndex2(nr,-1);
1676 // if (TMath::Abs(z)<(AliTPCReconstructor::GetCtgRange()*x+10) && TMath::Abs(z)<fkParam->GetZLength(0) && (TMath::Abs(t.GetSnp())<AliTPCReconstructor::GetMaxSnpTracker()))
1677 t.SetNFoundable(t.GetNFoundable()+1);
1683 // cl = krow.FindNearest2(y+10.,z,roady,roadz,index);
1684 cl = krow.FindNearest2(y,z,roady,roadz,index);
1685 if (cl) t.SetCurrentClusterIndex1(krow.GetIndex(index));
1688 t.SetCurrentCluster(cl);
1690 if (fIteration==2&&cl->IsUsed(10)) return 0;
1691 Int_t accept = AcceptCluster(&t,t.GetCurrentCluster());
1692 if (fIteration==2&&cl->IsUsed(11)) {
1693 t.SetErrorY2(t.GetErrorY2()+0.03);
1694 t.SetErrorZ2(t.GetErrorZ2()+0.03);
1695 t.SetErrorY2(t.GetErrorY2()*3);
1696 t.SetErrorZ2(t.GetErrorZ2()*3);
1699 if (t.fCurrentCluster->IsUsed(10)){
1704 if (t.fNShared>0.7*t.GetNumberOfClusters()) {
1710 if (accept<3) UpdateTrack(&t,accept);
1713 if ( fIteration==0 && t.GetNFoundable()*0.5 > t.GetNumberOfClusters()) t.SetRemoval(10);
1721 //_________________________________________________________________________
1722 Bool_t AliTPCtrackerMI::GetTrackPoint(Int_t index, AliTrackPoint &p ) const
1724 // Get track space point by index
1725 // return false in case the cluster doesn't exist
1726 AliTPCclusterMI *cl = GetClusterMI(index);
1727 if (!cl) return kFALSE;
1728 Int_t sector = (index&0xff000000)>>24;
1729 // Int_t row = (index&0x00ff0000)>>16;
1731 // xyz[0] = fkParam->GetPadRowRadii(sector,row);
1732 xyz[0] = cl->GetX();
1733 xyz[1] = cl->GetY();
1734 xyz[2] = cl->GetZ();
1736 fkParam->AdjustCosSin(sector,cos,sin);
1737 Float_t x = cos*xyz[0]-sin*xyz[1];
1738 Float_t y = cos*xyz[1]+sin*xyz[0];
1740 Float_t sigmaY2 = 0.027*cl->GetSigmaY2();
1741 if (sector < fkParam->GetNInnerSector()) sigmaY2 *= 2.07;
1742 Float_t sigmaZ2 = 0.066*cl->GetSigmaZ2();
1743 if (sector < fkParam->GetNInnerSector()) sigmaZ2 *= 1.77;
1744 cov[0] = sin*sin*sigmaY2;
1745 cov[1] = -sin*cos*sigmaY2;
1747 cov[3] = cos*cos*sigmaY2;
1750 p.SetXYZ(x,y,xyz[2],cov);
1751 AliGeomManager::ELayerID iLayer;
1753 if (sector < fkParam->GetNInnerSector()) {
1754 iLayer = AliGeomManager::kTPC1;
1758 iLayer = AliGeomManager::kTPC2;
1759 idet = sector - fkParam->GetNInnerSector();
1761 UShort_t volid = AliGeomManager::LayerToVolUID(iLayer,idet);
1762 p.SetVolumeID(volid);
1768 Int_t AliTPCtrackerMI::UpdateClusters(AliTPCseed& t, Int_t nr) {
1769 //-----------------------------------------------------------------
1770 // This function tries to find a track prolongation to next pad row
1771 //-----------------------------------------------------------------
1772 t.SetCurrentCluster(0);
1773 t.SetCurrentClusterIndex1(0);
1775 Double_t xt=t.GetX();
1776 Int_t row = GetRowNumber(xt)-1;
1777 Double_t ymax= GetMaxY(nr);
1779 if (row < nr) return 1; // don't prolongate if not information until now -
1780 // if (TMath::Abs(t.GetSnp())>0.9 && t.GetNumberOfClusters()>40. && fIteration!=2) {
1782 // return 0; // not prolongate strongly inclined tracks
1784 // if (TMath::Abs(t.GetSnp())>0.95) {
1786 // return 0; // not prolongate strongly inclined tracks
1787 // }// patch 28 fev 06
1789 Double_t x= GetXrow(nr);
1791 //t.PropagateTo(x+0.02);
1792 //t.PropagateTo(x+0.01);
1793 if (!t.PropagateTo(x)){
1800 if (TMath::Abs(y)>ymax){
1802 t.SetRelativeSector((t.GetRelativeSector()+1) % fN);
1803 if (!t.Rotate(fSectors->GetAlpha()))
1805 } else if (y <-ymax) {
1806 t.SetRelativeSector((t.GetRelativeSector()-1+fN) % fN);
1807 if (!t.Rotate(-fSectors->GetAlpha()))
1810 // if (!t.PropagateTo(x)){
1817 if (TMath::Abs(t.GetSnp())>AliTPCReconstructor::GetMaxSnpTracker()) return 0;
1819 if (!IsActive(t.GetRelativeSector(),nr)) {
1821 t.SetClusterIndex2(nr,-1);
1824 //AliInfo(Form("A - Sector%d phi %f - alpha %f", t.fRelativeSector,y/x, t.GetAlpha()));
1826 AliTPCtrackerRow &krow=GetRow(t.GetRelativeSector(),nr);
1828 if (TMath::Abs(TMath::Abs(y)-ymax)<krow.GetDeadZone()){
1830 t.SetClusterIndex2(nr,-1);
1836 // if (TMath::Abs(t.GetZ())<(AliTPCReconstructor::GetCtgRange()*t.GetX()+10) && (TMath::Abs(t.GetSnp())<AliTPCReconstructor::GetMaxSnpTracker()))
1837 if (IsFindable(t)) t.SetNFoundable(t.GetNFoundable()+1);
1843 if ( (nr%2==0) || t.GetNumberOfClusters()<2){
1844 // t.fCurrentSigmaY = GetSigmaY(&t);
1845 //t.fCurrentSigmaZ = GetSigmaZ(&t);
1849 AliTPCclusterMI *cl=0;
1852 Double_t roady = 1.;
1853 Double_t roadz = 1.;
1857 index = t.GetClusterIndex2(nr);
1858 if ( (index>0) && (index&0x8000)==0){
1859 cl = t.GetClusterPointer(nr);
1860 if ( (cl==0) && (index>0)) cl = GetClusterMI(index);
1861 t.SetCurrentClusterIndex1(index);
1863 t.SetCurrentCluster(cl);
1869 // if (index<0) return 0;
1870 UInt_t uindex = TMath::Abs(index);
1873 //cl = krow.FindNearest2(y+10,z,roady,roadz,uindex);
1874 cl = krow.FindNearest2(y,z,roady,roadz,uindex);
1877 if (cl) t.SetCurrentClusterIndex1(krow.GetIndex(uindex));
1878 t.SetCurrentCluster(cl);
1884 Int_t AliTPCtrackerMI::FollowToNextCluster(AliTPCseed & t, Int_t nr) {
1885 //-----------------------------------------------------------------
1886 // This function tries to find a track prolongation to next pad row
1887 //-----------------------------------------------------------------
1889 //update error according neighborhoud
1891 if (t.GetCurrentCluster()) {
1893 Int_t accept = AcceptCluster(&t,t.GetCurrentCluster());
1895 if (t.GetCurrentCluster()->IsUsed(10)){
1900 t.SetNShared(t.GetNShared()+1);
1901 if (t.GetNShared()>0.7*t.GetNumberOfClusters()) {
1906 if (fIteration>0) accept = 0;
1907 if (accept<3) UpdateTrack(&t,accept);
1911 if ( ( (t.GetSigmaY2()+t.GetSigmaZ2())>0.16)&& t.GetNumberOfClusters()>18) t.SetRemoval(10);
1912 if ( t.GetChi2()/t.GetNumberOfClusters()>6 &&t.GetNumberOfClusters()>18) t.SetRemoval(10);
1914 if (( (t.GetNFoundable()*0.5 > t.GetNumberOfClusters()) || t.GetNoCluster()>15)) t.SetRemoval(10);
1922 //_____________________________________________________________________________
1923 Int_t AliTPCtrackerMI::FollowProlongation(AliTPCseed& t, Int_t rf, Int_t step) {
1924 //-----------------------------------------------------------------
1925 // This function tries to find a track prolongation.
1926 //-----------------------------------------------------------------
1927 Double_t xt=t.GetX();
1929 Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
1930 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
1931 if (alpha < 0. ) alpha += 2.*TMath::Pi();
1933 t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN);
1935 Int_t first = GetRowNumber(xt)-1;
1936 for (Int_t nr= first; nr>=rf; nr-=step) {
1938 if (t.GetKinkIndexes()[0]>0){
1939 for (Int_t i=0;i<3;i++){
1940 Int_t index = t.GetKinkIndexes()[i];
1941 if (index==0) break;
1942 if (index<0) continue;
1944 AliKink * kink = (AliKink*)fEvent->GetKink(index-1);
1946 printf("PROBLEM\n");
1949 Int_t kinkrow = kink->GetTPCRow0()+2+Int_t(0.5/(0.05+kink->GetAngle(2)));
1951 AliExternalTrackParam paramd(t);
1952 kink->SetDaughter(paramd);
1953 kink->SetStatus(2,5);
1960 if (nr==80) t.UpdateReference();
1961 if (nr<fInnerSec->GetNRows())
1962 fSectors = fInnerSec;
1964 fSectors = fOuterSec;
1965 if (FollowToNext(t,nr)==0)
1978 Int_t AliTPCtrackerMI::FollowBackProlongation(AliTPCseed& t, Int_t rf) {
1979 //-----------------------------------------------------------------
1980 // This function tries to find a track prolongation.
1981 //-----------------------------------------------------------------
1983 Double_t xt=t.GetX();
1984 Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
1985 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
1986 if (alpha < 0. ) alpha += 2.*TMath::Pi();
1987 t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN);
1989 Int_t first = t.GetFirstPoint();
1990 if (first<GetRowNumber(xt)+1) first = GetRowNumber(xt)+1;
1992 if (first<0) first=0;
1993 for (Int_t nr=first; nr<=rf; nr++) {
1994 // if ( (TMath::Abs(t.GetSnp())>0.95)) break;//patch 28 fev 06
1995 if (t.GetKinkIndexes()[0]<0){
1996 for (Int_t i=0;i<3;i++){
1997 Int_t index = t.GetKinkIndexes()[i];
1998 if (index==0) break;
1999 if (index>0) continue;
2000 index = TMath::Abs(index);
2001 AliKink * kink = (AliKink*)fEvent->GetKink(index-1);
2003 printf("PROBLEM\n");
2006 Int_t kinkrow = kink->GetTPCRow0()-2-Int_t(0.5/(0.05+kink->GetAngle(2)));
2008 AliExternalTrackParam paramm(t);
2009 kink->SetMother(paramm);
2010 kink->SetStatus(2,1);
2017 if (nr<fInnerSec->GetNRows())
2018 fSectors = fInnerSec;
2020 fSectors = fOuterSec;
2031 Float_t AliTPCtrackerMI::OverlapFactor(AliTPCseed * s1, AliTPCseed * s2, Int_t &sum1, Int_t & sum2)
2039 Float_t dz2 =(s1->GetZ() - s2->GetZ());
2042 Float_t dy2 =TMath::Abs((s1->GetY() - s2->GetY()));
2044 Float_t distance = TMath::Sqrt(dz2+dy2);
2045 if (distance>4.) return 0; // if there are far away - not overlap - to reduce combinatorics
2048 Int_t firstpoint = TMath::Min(s1->GetFirstPoint(),s2->GetFirstPoint());
2049 Int_t lastpoint = TMath::Max(s1->GetLastPoint(),s2->GetLastPoint());
2054 if (firstpoint>lastpoint) {
2055 firstpoint =lastpoint;
2060 for (Int_t i=firstpoint-1;i<lastpoint+1;i++){
2061 if (s1->GetClusterIndex2(i)>0) sum1++;
2062 if (s2->GetClusterIndex2(i)>0) sum2++;
2063 if (s1->GetClusterIndex2(i)==s2->GetClusterIndex2(i) && s1->GetClusterIndex2(i)>0) {
2067 if (sum<5) return 0;
2069 Float_t summin = TMath::Min(sum1+1,sum2+1);
2070 Float_t ratio = (sum+1)/Float_t(summin);
2074 void AliTPCtrackerMI::SignShared(AliTPCseed * s1, AliTPCseed * s2)
2078 Float_t thetaCut = 0.2;//+10.*TMath::Sqrt(s1->GetSigmaTglZ()+ s2->GetSigmaTglZ());
2079 if (TMath::Abs(s1->GetTgl()-s2->GetTgl())>thetaCut) return;
2080 Float_t minCl = TMath::Min(s1->GetNumberOfClusters(),s2->GetNumberOfClusters());
2081 Int_t cutN0 = TMath::Max(5,TMath::Nint(0.1*minCl));
2086 //Int_t firstpoint = TMath::Max(s1->GetFirstPoint(),s2->GetFirstPoint());
2087 //Int_t lastpoint = TMath::Min(s1->GetLastPoint(),s2->GetLastPoint());
2088 Int_t firstpoint = 0;
2089 Int_t lastpoint = 160;
2091 // if (firstpoint>=lastpoint-5) return;;
2093 for (Int_t i=firstpoint;i<lastpoint;i++){
2094 // if ( (s1->GetClusterIndex2(i)&0xFFFF8FFF)==(s2->GetClusterIndex2(i)&0xFFFF8FFF) && s1->GetClusterIndex2(i)>0) {
2095 if ( (s1->GetClusterIndex2(i))==(s2->GetClusterIndex2(i)) && s1->GetClusterIndex2(i)>0) {
2097 s1->SetSharedMapBit(i, kTRUE);
2098 s2->SetSharedMapBit(i, kTRUE);
2100 if (s1->GetClusterIndex2(i)>0)
2101 s1->SetClusterMapBit(i, kTRUE);
2103 if (sumshared>cutN0){
2106 for (Int_t i=firstpoint;i<lastpoint;i++){
2107 // if ( (s1->GetClusterIndex2(i)&0xFFFF8FFF)==(s2->GetClusterIndex2(i)&0xFFFF8FFF) && s1->GetClusterIndex2(i)>0) {
2108 if ( (s1->GetClusterIndex2(i))==(s2->GetClusterIndex2(i)) && s1->GetClusterIndex2(i)>0) {
2109 AliTPCTrackerPoint *p1 = s1->GetTrackPoint(i);
2110 AliTPCTrackerPoint *p2 = s2->GetTrackPoint(i);;
2111 if (s1->IsActive()&&s2->IsActive()){
2112 p1->SetShared(kTRUE);
2113 p2->SetShared(kTRUE);
2119 if (sumshared>cutN0){
2120 for (Int_t i=0;i<4;i++){
2121 if (s1->GetOverlapLabel(3*i)==0){
2122 s1->SetOverlapLabel(3*i, s2->GetLabel());
2123 s1->SetOverlapLabel(3*i+1,sumshared);
2124 s1->SetOverlapLabel(3*i+2,s2->GetUniqueID());
2128 for (Int_t i=0;i<4;i++){
2129 if (s2->GetOverlapLabel(3*i)==0){
2130 s2->SetOverlapLabel(3*i, s1->GetLabel());
2131 s2->SetOverlapLabel(3*i+1,sumshared);
2132 s2->SetOverlapLabel(3*i+2,s1->GetUniqueID());
2139 void AliTPCtrackerMI::SignShared(TObjArray * arr)
2142 //sort trackss according sectors
2144 for (Int_t i=0; i<arr->GetEntriesFast(); i++) {
2145 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2147 //if (pt) RotateToLocal(pt);
2151 arr->Sort(); // sorting according relative sectors
2152 arr->Expand(arr->GetEntries());
2155 Int_t nseed=arr->GetEntriesFast();
2156 for (Int_t i=0; i<nseed; i++) {
2157 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2159 for (Int_t j=0;j<12;j++){
2160 pt->SetOverlapLabel(j,0);
2163 for (Int_t i=0; i<nseed; i++) {
2164 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2166 if (pt->GetRemoval()>10) continue;
2167 for (Int_t j=i+1; j<nseed; j++){
2168 AliTPCseed *pt2=(AliTPCseed*)arr->UncheckedAt(j);
2169 if (TMath::Abs(pt->GetRelativeSector()-pt2->GetRelativeSector())>1) continue;
2171 if (pt2->GetRemoval()<=10) {
2172 //if ( TMath::Abs(pt->GetRelativeSector()-pt2->GetRelativeSector())>0) break;
2180 void AliTPCtrackerMI::SortTracks(TObjArray * arr, Int_t mode) const
2183 //sort tracks in array according mode criteria
2184 Int_t nseed = arr->GetEntriesFast();
2185 for (Int_t i=0; i<nseed; i++) {
2186 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2197 void AliTPCtrackerMI::RemoveUsed2(TObjArray * arr, Float_t factor1, Float_t factor2, Int_t minimal)
2200 // Loop over all tracks and remove overlaped tracks (with lower quality)
2202 // 1. Unsign clusters
2203 // 2. Sort tracks according quality
2204 // Quality is defined by the number of cluster between first and last points
2206 // 3. Loop over tracks - decreasing quality order
2207 // a.) remove - If the fraction of shared cluster less than factor (1- n or 2)
2208 // b.) remove - If the minimal number of clusters less than minimal and not ITS
2209 // c.) if track accepted - sign clusters
2211 //Called in - AliTPCtrackerMI::Clusters2Tracks()
2212 // - AliTPCtrackerMI::PropagateBack()
2213 // - AliTPCtrackerMI::RefitInward()
2216 // factor1 - factor for constrained
2217 // factor2 - for non constrained tracks
2218 // if (Float_t(shared+1)/Float_t(found+1)>factor) - DELETE
2222 Int_t nseed = arr->GetEntriesFast();
2223 Float_t * quality = new Float_t[nseed];
2224 Int_t * indexes = new Int_t[nseed];
2228 for (Int_t i=0; i<nseed; i++) {
2229 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2234 pt->UpdatePoints(); //select first last max dens points
2235 Float_t * points = pt->GetPoints();
2236 if (points[3]<0.8) quality[i] =-1;
2237 quality[i] = (points[2]-points[0])+pt->GetNumberOfClusters();
2238 //prefer high momenta tracks if overlaps
2239 quality[i] *= TMath::Sqrt(TMath::Abs(pt->Pt())+0.5);
2241 TMath::Sort(nseed,quality,indexes);
2244 for (Int_t itrack=0; itrack<nseed; itrack++) {
2245 Int_t trackindex = indexes[itrack];
2246 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(trackindex);
2249 if (quality[trackindex]<0){
2250 delete arr->RemoveAt(trackindex);
2255 Int_t first = Int_t(pt->GetPoints()[0]);
2256 Int_t last = Int_t(pt->GetPoints()[2]);
2257 Double_t factor = (pt->GetBConstrain()) ? factor1: factor2;
2259 Int_t found,foundable,shared;
2260 pt->GetClusterStatistic(first,last, found, foundable,shared,kFALSE); // better to get statistic in "high-dens" region do't use full track as in line bellow
2261 // pt->GetClusterStatistic(0,160, found, foundable,shared,kFALSE);
2262 Bool_t itsgold =kFALSE;
2265 if (pt->GetESD()->GetITSclusters(dummy)>4) itsgold= kTRUE;
2269 if (Float_t(shared+1)/Float_t(found+1)>factor){
2270 if (pt->GetKinkIndexes()[0]!=0) continue; //don't remove tracks - part of the kinks
2271 if( AliTPCReconstructor::StreamLevel()>3){
2272 TTreeSRedirector &cstream = *fDebugStreamer;
2273 cstream<<"RemoveUsed"<<
2274 "iter="<<fIteration<<
2278 delete arr->RemoveAt(trackindex);
2281 if (pt->GetNumberOfClusters()<50&&(found-0.5*shared)<minimal){ //remove short tracks
2282 if (pt->GetKinkIndexes()[0]!=0) continue; //don't remove tracks - part of the kinks
2283 if( AliTPCReconstructor::StreamLevel()>3){
2284 TTreeSRedirector &cstream = *fDebugStreamer;
2285 cstream<<"RemoveShort"<<
2286 "iter="<<fIteration<<
2290 delete arr->RemoveAt(trackindex);
2296 //if (sharedfactor>0.4) continue;
2297 if (pt->GetKinkIndexes()[0]>0) continue;
2298 //Remove tracks with undefined properties - seems
2299 if (pt->GetSigmaY2()<kAlmost0) continue; // ? what is the origin ?
2301 for (Int_t i=first; i<last; i++) {
2302 Int_t index=pt->GetClusterIndex2(i);
2303 // if (index<0 || index&0x8000 ) continue;
2304 if (index<0 || index&0x8000 ) continue;
2305 AliTPCclusterMI *c= pt->GetClusterPointer(i);
2312 Info("RemoveUsed","\n*****\nNumber of good tracks after shared removal\t%d\n",fNtracks);
2318 void AliTPCtrackerMI::DumpClusters(Int_t iter, TObjArray *trackArray)
2321 // Dump clusters after reco
2322 // signed and unsigned cluster can be visualized
2323 // 1. Unsign all cluster
2324 // 2. Sign all used clusters
2327 Int_t nseed = trackArray->GetEntries();
2328 for (Int_t i=0; i<nseed; i++){
2329 AliTPCseed *pt=(AliTPCseed*)trackArray->UncheckedAt(i);
2333 Bool_t isKink=pt->GetKinkIndex(0)!=0;
2334 for (Int_t j=0; j<160; ++j) {
2335 Int_t index=pt->GetClusterIndex2(j);
2336 if (index<0) continue;
2337 AliTPCclusterMI *c= pt->GetClusterPointer(j);
2339 if (isKink) c->Use(100); // kink
2340 c->Use(10); // by default usage 10
2345 for (Int_t sec=0;sec<fkNIS;sec++){
2346 for (Int_t row=0;row<fInnerSec->GetNRows();row++){
2347 AliTPCclusterMI *cl = fInnerSec[sec][row].GetClusters1();
2348 for (Int_t icl =0;icl< fInnerSec[sec][row].GetN1();icl++){
2349 Float_t gx[3]; cl[icl].GetGlobalXYZ(gx);
2350 (*fDebugStreamer)<<"clDump"<<
2358 cl = fInnerSec[sec][row].GetClusters2();
2359 for (Int_t icl =0;icl< fInnerSec[sec][row].GetN2();icl++){
2360 Float_t gx[3]; cl[icl].GetGlobalXYZ(gx);
2361 (*fDebugStreamer)<<"clDump"<<
2372 for (Int_t sec=0;sec<fkNOS;sec++){
2373 for (Int_t row=0;row<fOuterSec->GetNRows();row++){
2374 AliTPCclusterMI *cl = fOuterSec[sec][row].GetClusters1();
2375 for (Int_t icl =0;icl< fOuterSec[sec][row].GetN1();icl++){
2376 Float_t gx[3]; cl[icl].GetGlobalXYZ(gx);
2377 (*fDebugStreamer)<<"clDump"<<
2385 cl = fOuterSec[sec][row].GetClusters2();
2386 for (Int_t icl =0;icl< fOuterSec[sec][row].GetN2();icl++){
2387 Float_t gx[3]; cl[icl].GetGlobalXYZ(gx);
2388 (*fDebugStreamer)<<"clDump"<<
2400 void AliTPCtrackerMI::UnsignClusters()
2403 // loop over all clusters and unsign them
2406 for (Int_t sec=0;sec<fkNIS;sec++){
2407 for (Int_t row=0;row<fInnerSec->GetNRows();row++){
2408 AliTPCclusterMI *cl = fInnerSec[sec][row].GetClusters1();
2409 for (Int_t icl =0;icl< fInnerSec[sec][row].GetN1();icl++)
2410 // if (cl[icl].IsUsed(10))
2412 cl = fInnerSec[sec][row].GetClusters2();
2413 for (Int_t icl =0;icl< fInnerSec[sec][row].GetN2();icl++)
2414 //if (cl[icl].IsUsed(10))
2419 for (Int_t sec=0;sec<fkNOS;sec++){
2420 for (Int_t row=0;row<fOuterSec->GetNRows();row++){
2421 AliTPCclusterMI *cl = fOuterSec[sec][row].GetClusters1();
2422 for (Int_t icl =0;icl< fOuterSec[sec][row].GetN1();icl++)
2423 //if (cl[icl].IsUsed(10))
2425 cl = fOuterSec[sec][row].GetClusters2();
2426 for (Int_t icl =0;icl< fOuterSec[sec][row].GetN2();icl++)
2427 //if (cl[icl].IsUsed(10))
2436 void AliTPCtrackerMI::SignClusters(const TObjArray * arr, Float_t fnumber, Float_t fdensity)
2439 //sign clusters to be "used"
2441 // snumber and sdensity sign number of sigmas - bellow mean value to be accepted
2442 // loop over "primaries"
2456 Int_t nseed = arr->GetEntriesFast();
2457 for (Int_t i=0; i<nseed; i++) {
2458 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2462 if (!(pt->IsActive())) continue;
2463 Float_t dens = pt->GetNumberOfClusters()/Float_t(pt->GetNFoundable());
2464 if ( (dens>0.7) && (pt->GetNumberOfClusters()>70)){
2466 sumdens2+= dens*dens;
2467 sumn += pt->GetNumberOfClusters();
2468 sumn2 += pt->GetNumberOfClusters()*pt->GetNumberOfClusters();
2469 Float_t chi2 = pt->GetChi2()/pt->GetNumberOfClusters();
2472 sumchi2 +=chi2*chi2;
2477 Float_t mdensity = 0.9;
2478 Float_t meann = 130;
2479 Float_t meanchi = 1;
2480 Float_t sdensity = 0.1;
2481 Float_t smeann = 10;
2482 Float_t smeanchi =0.4;
2486 mdensity = sumdens/sum;
2488 meanchi = sumchi/sum;
2490 sdensity = sumdens2/sum-mdensity*mdensity;
2492 sdensity = TMath::Sqrt(sdensity);
2496 smeann = sumn2/sum-meann*meann;
2498 smeann = TMath::Sqrt(smeann);
2502 smeanchi = sumchi2/sum - meanchi*meanchi;
2504 smeanchi = TMath::Sqrt(smeanchi);
2510 //REMOVE SHORT DELTAS or tracks going out of sensitive volume of TPC
2512 for (Int_t i=0; i<nseed; i++) {
2513 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2517 if (pt->GetBSigned()) continue;
2518 if (pt->GetBConstrain()) continue;
2519 //if (!(pt->IsActive())) continue;
2521 Int_t found,foundable,shared;
2522 pt->GetClusterStatistic(0,160,found, foundable,shared);
2523 if (shared/float(found)>0.3) {
2524 if (shared/float(found)>0.9 ){
2525 //delete arr->RemoveAt(i);
2530 Bool_t isok =kFALSE;
2531 if ( (pt->GetNShared()/pt->GetNumberOfClusters()<0.5) &&pt->GetNumberOfClusters()>60)
2533 if ((TMath::Abs(1/pt->GetC())<100.) && (pt->GetNShared()/pt->GetNumberOfClusters()<0.7))
2535 if (TMath::Abs(pt->GetZ()/pt->GetX())>1.1)
2537 if ( (TMath::Abs(pt->GetSnp()>0.7) && pt->GetD(0,0)>60.))
2541 for (Int_t j=0; j<160; ++j) {
2542 Int_t index=pt->GetClusterIndex2(j);
2543 if (index<0) continue;
2544 AliTPCclusterMI *c= pt->GetClusterPointer(j);
2546 //if (!(c->IsUsed(10))) c->Use();
2553 Double_t maxchi = meanchi+2.*smeanchi;
2555 for (Int_t i=0; i<nseed; i++) {
2556 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2560 //if (!(pt->IsActive())) continue;
2561 if (pt->GetBSigned()) continue;
2562 Double_t chi = pt->GetChi2()/pt->GetNumberOfClusters();
2563 if (chi>maxchi) continue;
2566 Float_t dens = pt->GetNumberOfClusters()/Float_t(pt->GetNFoundable());
2568 //sign only tracks with enoug big density at the beginning
2570 if ((pt->GetDensityFirst(40)<0.75) && pt->GetNumberOfClusters()<meann) continue;
2573 Double_t mindens = TMath::Max(double(mdensity-sdensity*fdensity*bfactor),0.65);
2574 Double_t minn = TMath::Max(Int_t(meann-fnumber*smeann*bfactor),50);
2576 // if (pt->fBConstrain) mindens = TMath::Max(mdensity-sdensity*fdensity*bfactor,0.65);
2577 if ( (pt->GetRemoval()==10) && (pt->GetSnp()>0.8)&&(dens>mindens))
2580 if ((dens>mindens && pt->GetNumberOfClusters()>minn) && chi<maxchi ){
2581 //Int_t noc=pt->GetNumberOfClusters();
2582 pt->SetBSigned(kTRUE);
2583 for (Int_t j=0; j<160; ++j) {
2585 Int_t index=pt->GetClusterIndex2(j);
2586 if (index<0) continue;
2587 AliTPCclusterMI *c= pt->GetClusterPointer(j);
2589 // if (!(c->IsUsed(10))) c->Use();
2594 // gLastCheck = nseed;
2602 void AliTPCtrackerMI::StopNotActive(const TObjArray * arr, Int_t row0, Float_t th0, Float_t th1, Float_t th2) const
2604 // stop not active tracks
2605 // take th1 as threshold for number of founded to number of foundable on last 10 active rows
2606 // take th2 as threshold for number of founded to number of foundable on last 20 active rows
2607 Int_t nseed = arr->GetEntriesFast();
2609 for (Int_t i=0; i<nseed; i++) {
2610 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2614 if (!(pt->IsActive())) continue;
2615 StopNotActive(pt,row0,th0, th1,th2);
2621 void AliTPCtrackerMI::StopNotActive(AliTPCseed * seed, Int_t row0, Float_t th0, Float_t th1,
2624 // stop not active tracks
2625 // take th1 as threshold for number of founded to number of foundable on last 10 active rows
2626 // take th2 as threshold for number of founded to number of foundable on last 20 active rows
2629 Int_t foundable = 0;
2630 Int_t maxindex = seed->GetLastPoint(); //last foundable row
2631 if (seed->GetNFoundable()*th0 > seed->GetNumberOfClusters()) {
2632 seed->Desactivate(10) ;
2636 for (Int_t i=row0; i<maxindex; i++){
2637 Int_t index = seed->GetClusterIndex2(i);
2638 if (index!=-1) foundable++;
2640 if (foundable<=30) sumgood1++;
2641 if (foundable<=50) {
2648 if (foundable>=30.){
2649 if (sumgood1<(th1*30.)) seed->Desactivate(10);
2652 if (sumgood2<(th2*50.)) seed->Desactivate(10);
2656 Int_t AliTPCtrackerMI::RefitInward(AliESDEvent *event)
2659 // back propagation of ESD tracks
2662 if (!event) return 0;
2663 const Int_t kMaxFriendTracks=2000;
2667 //PrepareForProlongation(fSeeds,1);
2668 PropagateForward2(fSeeds);
2669 RemoveUsed2(fSeeds,0.4,0.4,20);
2671 TObjArray arraySeed(fSeeds->GetEntries());
2672 for (Int_t i=0;i<fSeeds->GetEntries();i++) {
2673 arraySeed.AddAt(fSeeds->At(i),i);
2675 SignShared(&arraySeed);
2676 // FindCurling(fSeeds, event,2); // find multi found tracks
2677 FindSplitted(fSeeds, event,2); // find multi found tracks
2678 if (AliTPCReconstructor::StreamLevel()>5) FindMultiMC(fSeeds, fEvent,2); // find multi found tracks
2681 Int_t nseed = fSeeds->GetEntriesFast();
2682 for (Int_t i=0;i<nseed;i++){
2683 AliTPCseed * seed = (AliTPCseed*) fSeeds->UncheckedAt(i);
2684 if (!seed) continue;
2685 if (seed->GetKinkIndex(0)>0) UpdateKinkQualityD(seed); // update quality informations for kinks
2686 AliESDtrack *esd=event->GetTrack(i);
2688 if (seed->GetNumberOfClusters()<60 && seed->GetNumberOfClusters()<(esd->GetTPCclusters(0) -5)*0.8){
2689 AliExternalTrackParam paramIn;
2690 AliExternalTrackParam paramOut;
2691 Int_t ncl = seed->RefitTrack(seed,¶mIn,¶mOut);
2692 if (AliTPCReconstructor::StreamLevel()>2) {
2693 (*fDebugStreamer)<<"RecoverIn"<<
2697 "pout.="<<¶mOut<<
2702 seed->Set(paramIn.GetX(),paramIn.GetAlpha(),paramIn.GetParameter(),paramIn.GetCovariance());
2703 seed->SetNumberOfClusters(ncl);
2707 seed->PropagateTo(fkParam->GetInnerRadiusLow());
2708 seed->UpdatePoints();
2709 AddCovariance(seed);
2711 seed->CookdEdx(0.02,0.6);
2712 CookLabel(seed,0.1); //For comparison only
2713 esd->SetTPCClusterMap(seed->GetClusterMap());
2714 esd->SetTPCSharedMap(seed->GetSharedMap());
2716 if (AliTPCReconstructor::StreamLevel()>1 && seed!=0) {
2717 TTreeSRedirector &cstream = *fDebugStreamer;
2724 if (seed->GetNumberOfClusters()>15){
2725 esd->UpdateTrackParams(seed,AliESDtrack::kTPCrefit);
2726 esd->SetTPCPoints(seed->GetPoints());
2727 esd->SetTPCPointsF(seed->GetNFoundable());
2728 Int_t ndedx = seed->GetNCDEDX(0);
2729 Float_t sdedx = seed->GetSDEDX(0);
2730 Float_t dedx = seed->GetdEdx();
2731 esd->SetTPCsignal(dedx, sdedx, ndedx);
2733 // add seed to the esd track in Calib level
2735 Bool_t storeFriend = gRandom->Rndm()<(kMaxFriendTracks)/Float_t(nseed);
2736 if (AliTPCReconstructor::StreamLevel()>0 &&storeFriend){
2737 AliTPCseed * seedCopy = new AliTPCseed(*seed, kTRUE);
2738 esd->AddCalibObject(seedCopy);
2743 //printf("problem\n");
2746 //FindKinks(fSeeds,event);
2747 if (AliTPCReconstructor::StreamLevel()>3) DumpClusters(2,fSeeds);
2748 Info("RefitInward","Number of refitted tracks %d",ntracks);
2753 Int_t AliTPCtrackerMI::PropagateBack(AliESDEvent *event)
2756 // back propagation of ESD tracks
2758 if (!event) return 0;
2762 PropagateBack(fSeeds);
2763 RemoveUsed2(fSeeds,0.4,0.4,20);
2764 //FindCurling(fSeeds, fEvent,1);
2765 FindSplitted(fSeeds, event,1); // find multi found tracks
2766 if (AliTPCReconstructor::StreamLevel()>5) FindMultiMC(fSeeds, fEvent,1); // find multi found tracks
2769 Int_t nseed = fSeeds->GetEntriesFast();
2771 for (Int_t i=0;i<nseed;i++){
2772 AliTPCseed * seed = (AliTPCseed*) fSeeds->UncheckedAt(i);
2773 if (!seed) continue;
2774 if (seed->GetKinkIndex(0)<0) UpdateKinkQualityM(seed); // update quality informations for kinks
2775 seed->UpdatePoints();
2776 AddCovariance(seed);
2777 AliESDtrack *esd=event->GetTrack(i);
2778 if (seed->GetNumberOfClusters()<60 && seed->GetNumberOfClusters()<(esd->GetTPCclusters(0) -5)*0.8){
2779 AliExternalTrackParam paramIn;
2780 AliExternalTrackParam paramOut;
2781 Int_t ncl = seed->RefitTrack(seed,¶mIn,¶mOut);
2782 if (AliTPCReconstructor::StreamLevel()>2) {
2783 (*fDebugStreamer)<<"RecoverBack"<<
2787 "pout.="<<¶mOut<<
2792 seed->Set(paramOut.GetX(),paramOut.GetAlpha(),paramOut.GetParameter(),paramOut.GetCovariance());
2793 seed->SetNumberOfClusters(ncl);
2796 seed->CookdEdx(0.02,0.6);
2797 CookLabel(seed,0.1); //For comparison only
2798 if (seed->GetNumberOfClusters()>15){
2799 esd->UpdateTrackParams(seed,AliESDtrack::kTPCout);
2800 esd->SetTPCPoints(seed->GetPoints());
2801 esd->SetTPCPointsF(seed->GetNFoundable());
2802 Int_t ndedx = seed->GetNCDEDX(0);
2803 Float_t sdedx = seed->GetSDEDX(0);
2804 Float_t dedx = seed->GetdEdx();
2805 esd->SetTPCsignal(dedx, sdedx, ndedx);
2807 Int_t eventnumber = event->GetEventNumberInFile();// patch 28 fev 06
2808 // This is most likely NOT the event number you'd like to use. It has nothing to do with the 'real' event number
2809 if (AliTPCReconstructor::StreamLevel()>1 && esd) {
2810 (*fDebugStreamer)<<"Cback"<<
2813 "EventNrInFile="<<eventnumber<<
2818 if (AliTPCReconstructor::StreamLevel()>3) DumpClusters(1,fSeeds);
2819 //FindKinks(fSeeds,event);
2820 Info("PropagateBack","Number of back propagated tracks %d",ntracks);
2827 void AliTPCtrackerMI::DeleteSeeds()
2832 Int_t nseed = fSeeds->GetEntriesFast();
2833 for (Int_t i=0;i<nseed;i++){
2834 AliTPCseed * seed = (AliTPCseed*)fSeeds->At(i);
2835 if (seed) delete fSeeds->RemoveAt(i);
2842 void AliTPCtrackerMI::ReadSeeds(AliESDEvent *const event, Int_t direction)
2845 //read seeds from the event
2847 Int_t nentr=event->GetNumberOfTracks();
2849 Info("ReadSeeds", "Number of ESD tracks: %d\n", nentr);
2854 fSeeds = new TObjArray(nentr);
2858 for (Int_t i=0; i<nentr; i++) {
2859 AliESDtrack *esd=event->GetTrack(i);
2860 ULong_t status=esd->GetStatus();
2861 if (!(status&AliESDtrack::kTPCin)) continue;
2862 AliTPCtrack t(*esd);
2863 t.SetNumberOfClusters(0);
2864 // AliTPCseed *seed = new AliTPCseed(t,t.GetAlpha());
2865 AliTPCseed *seed = new AliTPCseed(t/*,t.GetAlpha()*/);
2866 seed->SetUniqueID(esd->GetID());
2867 AddCovariance(seed); //add systematic ucertainty
2868 for (Int_t ikink=0;ikink<3;ikink++) {
2869 Int_t index = esd->GetKinkIndex(ikink);
2870 seed->GetKinkIndexes()[ikink] = index;
2871 if (index==0) continue;
2872 index = TMath::Abs(index);
2873 AliESDkink * kink = fEvent->GetKink(index-1);
2874 if (kink&&esd->GetKinkIndex(ikink)<0){
2875 if ((status & AliESDtrack::kTRDrefit) != 0) kink->SetStatus(1,2);
2876 if ((status & AliESDtrack::kITSout) != 0) kink->SetStatus(1,0);
2878 if (kink&&esd->GetKinkIndex(ikink)>0){
2879 if ((status & AliESDtrack::kTRDrefit) != 0) kink->SetStatus(1,6);
2880 if ((status & AliESDtrack::kITSout) != 0) kink->SetStatus(1,4);
2884 if (((status&AliESDtrack::kITSout)==0)&&(direction==1)) seed->ResetCovariance(10.);
2885 if ( direction ==2 &&(status & AliESDtrack::kTRDrefit) == 0 ) seed->ResetCovariance(10.);
2886 //if ( direction ==2 && ((status & AliESDtrack::kTPCout) == 0) ) {
2887 // fSeeds->AddAt(0,i);
2891 if ( direction ==2 &&(status & AliESDtrack::kTRDrefit) > 0 ) {
2892 Double_t par0[5],par1[5],alpha,x;
2893 esd->GetInnerExternalParameters(alpha,x,par0);
2894 esd->GetExternalParameters(x,par1);
2895 Double_t delta1 = TMath::Abs(par0[4]-par1[4])/(0.000000001+TMath::Abs(par0[4]+par1[4]));
2896 Double_t delta2 = TMath::Abs(par0[3]-par1[3]);
2898 if (esd->GetTRDncls()>0) trdchi2 = esd->GetTRDchi2()/esd->GetTRDncls();
2899 //reset covariance if suspicious
2900 if ( (delta1>0.1) || (delta2>0.006) ||trdchi2>7.)
2901 seed->ResetCovariance(10.);
2906 // rotate to the local coordinate system
2908 fSectors=fInnerSec; fN=fkNIS;
2909 Double_t alpha=seed->GetAlpha() - fSectors->GetAlphaShift();
2910 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
2911 if (alpha < 0. ) alpha += 2.*TMath::Pi();
2912 Int_t ns=Int_t(alpha/fSectors->GetAlpha())%fN;
2913 alpha =ns*fSectors->GetAlpha() + fSectors->GetAlphaShift();
2914 if (alpha<-TMath::Pi()) alpha += 2*TMath::Pi();
2915 if (alpha>=TMath::Pi()) alpha -= 2*TMath::Pi();
2916 alpha-=seed->GetAlpha();
2917 if (!seed->Rotate(alpha)) {
2923 if (esd->GetKinkIndex(0)<=0){
2924 for (Int_t irow=0;irow<160;irow++){
2925 Int_t index = seed->GetClusterIndex2(irow);
2928 AliTPCclusterMI * cl = GetClusterMI(index);
2929 seed->SetClusterPointer(irow,cl);
2931 if ((index & 0x8000)==0){
2932 cl->Use(10); // accepted cluster
2934 cl->Use(6); // close cluster not accepted
2937 Info("ReadSeeds","Not found cluster");
2942 fSeeds->AddAt(seed,i);
2948 //_____________________________________________________________________________
2949 void AliTPCtrackerMI::MakeSeeds3(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, Float_t cuts[4],
2950 Float_t deltay, Int_t ddsec) {
2951 //-----------------------------------------------------------------
2952 // This function creates track seeds.
2953 // SEEDING WITH VERTEX CONSTRAIN
2954 //-----------------------------------------------------------------
2955 // cuts[0] - fP4 cut
2956 // cuts[1] - tan(phi) cut
2957 // cuts[2] - zvertex cut
2958 // cuts[3] - fP3 cut
2966 Double_t x[5], c[15];
2967 // Int_t di = i1-i2;
2969 AliTPCseed * seed = new AliTPCseed();
2970 Double_t alpha=fSectors->GetAlpha(), shift=fSectors->GetAlphaShift();
2971 Double_t cs=cos(alpha), sn=sin(alpha);
2973 // Double_t x1 =fOuterSec->GetX(i1);
2974 //Double_t xx2=fOuterSec->GetX(i2);
2976 Double_t x1 =GetXrow(i1);
2977 Double_t xx2=GetXrow(i2);
2979 Double_t x3=GetX(), y3=GetY(), z3=GetZ();
2981 Int_t imiddle = (i2+i1)/2; //middle pad row index
2982 Double_t xm = GetXrow(imiddle); // radius of middle pad-row
2983 const AliTPCtrackerRow& krm=GetRow(sec,imiddle); //middle pad -row
2987 const AliTPCtrackerRow& kr1=GetRow(ns,i1);
2988 Double_t ymax = GetMaxY(i1)-kr1.GetDeadZone()-1.5;
2989 Double_t ymaxm = GetMaxY(imiddle)-kr1.GetDeadZone()-1.5;
2992 // change cut on curvature if it can't reach this layer
2993 // maximal curvature set to reach it
2994 Double_t dvertexmax = TMath::Sqrt((x1-x3)*(x1-x3)+(ymax+5-y3)*(ymax+5-y3));
2995 if (dvertexmax*0.5*cuts[0]>0.85){
2996 cuts[0] = 0.85/(dvertexmax*0.5+1.);
2998 Double_t r2min = 1/(cuts[0]*cuts[0]); //minimal square of radius given by cut
3001 if (deltay>0) ddsec = 0;
3002 // loop over clusters
3003 for (Int_t is=0; is < kr1; is++) {
3005 if (kr1[is]->IsUsed(10)) continue;
3006 Double_t y1=kr1[is]->GetY(), z1=kr1[is]->GetZ();
3007 //if (TMath::Abs(y1)>ymax) continue;
3009 if (deltay>0 && TMath::Abs(ymax-TMath::Abs(y1))> deltay ) continue; // seed only at the edge
3011 // find possible directions
3012 Float_t anglez = (z1-z3)/(x1-x3);
3013 Float_t extraz = z1 - anglez*(x1-xx2); // extrapolated z
3016 //find rotation angles relative to line given by vertex and point 1
3017 Double_t dvertex2 = (x1-x3)*(x1-x3)+(y1-y3)*(y1-y3);
3018 Double_t dvertex = TMath::Sqrt(dvertex2);
3019 Double_t angle13 = TMath::ATan((y1-y3)/(x1-x3));
3020 Double_t cs13 = cos(-angle13), sn13 = sin(-angle13);
3023 // loop over 2 sectors
3029 Double_t dddz1=0; // direction of delta inclination in z axis
3036 for (Int_t dsec = dsec1; dsec<=dsec2;dsec++){
3037 Int_t sec2 = sec + dsec;
3039 // AliTPCtrackerRow& kr2 = fOuterSec[(sec2+fkNOS)%fkNOS][i2];
3040 //AliTPCtrackerRow& kr2m = fOuterSec[(sec2+fkNOS)%fkNOS][imiddle];
3041 AliTPCtrackerRow& kr2 = GetRow((sec2+fkNOS)%fkNOS,i2);
3042 AliTPCtrackerRow& kr2m = GetRow((sec2+fkNOS)%fkNOS,imiddle);
3043 Int_t index1 = TMath::Max(kr2.Find(extraz-0.6-dddz1*TMath::Abs(z1)*0.05)-1,0);
3044 Int_t index2 = TMath::Min(kr2.Find(extraz+0.6+dddz2*TMath::Abs(z1)*0.05)+1,kr2);
3046 // rotation angles to p1-p3
3047 Double_t cs13r = cos(-angle13+dsec*alpha)/dvertex, sn13r = sin(-angle13+dsec*alpha)/dvertex;
3048 Double_t x2, y2, z2;
3050 // Double_t dymax = maxangle*TMath::Abs(x1-xx2);
3053 Double_t dxx0 = (xx2-x3)*cs13r;
3054 Double_t dyy0 = (xx2-x3)*sn13r;
3055 for (Int_t js=index1; js < index2; js++) {
3056 const AliTPCclusterMI *kcl = kr2[js];
3057 if (kcl->IsUsed(10)) continue;
3059 //calcutate parameters
3061 Double_t yy0 = dyy0 +(kcl->GetY()-y3)*cs13r;
3063 if (TMath::Abs(yy0)<0.000001) continue;
3064 Double_t xx0 = dxx0 -(kcl->GetY()-y3)*sn13r;
3065 Double_t y0 = 0.5*(xx0*xx0+yy0*yy0-xx0)/yy0;
3066 Double_t r02 = (0.25+y0*y0)*dvertex2;
3067 //curvature (radius) cut
3068 if (r02<r2min) continue;
3072 Double_t c0 = 1/TMath::Sqrt(r02);
3076 //Double_t dfi0 = 2.*TMath::ASin(dvertex*c0*0.5);
3077 //Double_t dfi1 = 2.*TMath::ASin(TMath::Sqrt(yy0*yy0+(1-xx0)*(1-xx0))*dvertex*c0*0.5);
3078 Double_t dfi0 = 2.*AliTPCFastMath::FastAsin(dvertex*c0*0.5);
3079 Double_t dfi1 = 2.*AliTPCFastMath::FastAsin(TMath::Sqrt(yy0*yy0+(1-xx0)*(1-xx0))*dvertex*c0*0.5);
3082 Double_t z0 = kcl->GetZ();
3083 Double_t zzzz2 = z1-(z1-z3)*dfi1/dfi0;
3084 if (TMath::Abs(zzzz2-z0)>0.5) continue;
3087 Double_t dip = (z1-z0)*c0/dfi1;
3088 Double_t x0 = (0.5*cs13+y0*sn13)*dvertex*c0;
3099 x2= xx2*cs-y2*sn*dsec;
3100 y2=+xx2*sn*dsec+y2*cs;
3110 // do we have cluster at the middle ?
3112 GetProlongation(x1,xm,x,ym,zm);
3114 AliTPCclusterMI * cm=0;
3115 if (TMath::Abs(ym)-ymaxm<0){
3116 cm = krm.FindNearest2(ym,zm,1.0,0.6,dummy);
3117 if ((!cm) || (cm->IsUsed(10))) {
3122 // rotate y1 to system 0
3123 // get state vector in rotated system
3124 Double_t yr1 = (-0.5*sn13+y0*cs13)*dvertex*c0;
3125 Double_t xr2 = x0*cs+yr1*sn*dsec;
3126 Double_t xr[5]={kcl->GetY(),kcl->GetZ(), xr2, dip, c0};
3128 GetProlongation(xx2,xm,xr,ym,zm);
3129 if (TMath::Abs(ym)-ymaxm<0){
3130 cm = kr2m.FindNearest2(ym,zm,1.0,0.6,dummy);
3131 if ((!cm) || (cm->IsUsed(10))) {
3141 dym = ym - cm->GetY();
3142 dzm = zm - cm->GetZ();
3149 Double_t sy1=kr1[is]->GetSigmaY2()*2., sz1=kr1[is]->GetSigmaZ2()*2.;
3150 Double_t sy2=kcl->GetSigmaY2()*2., sz2=kcl->GetSigmaZ2()*2.;
3151 //Double_t sy3=400*3./12., sy=0.1, sz=0.1;
3152 Double_t sy3=25000*x[4]*x[4]+0.1, sy=0.1, sz=0.1;
3153 //Double_t sy3=25000*x[4]*x[4]*60+0.5, sy=0.1, sz=0.1;
3155 Double_t f40=(F1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
3156 Double_t f42=(F1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
3157 Double_t f43=(F1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
3158 Double_t f20=(F2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
3159 Double_t f22=(F2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
3160 Double_t f23=(F2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
3162 Double_t f30=(F3(x1,y1+sy,x2,y2,z1,z2)-x[3])/sy;
3163 Double_t f31=(F3(x1,y1,x2,y2,z1+sz,z2)-x[3])/sz;
3164 Double_t f32=(F3(x1,y1,x2,y2+sy,z1,z2)-x[3])/sy;
3165 Double_t f34=(F3(x1,y1,x2,y2,z1,z2+sz)-x[3])/sz;
3169 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3170 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3171 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3172 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3173 c[13]=f30*sy1*f40+f32*sy2*f42;
3174 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3176 // if (!BuildSeed(kr1[is],kcl,0,x1,x2,x3,x,c)) continue;
3178 UInt_t index=kr1.GetIndex(is);
3179 seed->~AliTPCseed(); // this does not set the pointer to 0...
3180 AliTPCseed *track=new(seed) AliTPCseed(x1, ns*alpha+shift, x, c, index);
3182 track->SetIsSeeding(kTRUE);
3183 track->SetSeed1(i1);
3184 track->SetSeed2(i2);
3185 track->SetSeedType(3);
3189 FollowProlongation(*track, (i1+i2)/2,1);
3190 Int_t foundable,found,shared;
3191 track->GetClusterStatistic((i1+i2)/2,i1, found, foundable, shared, kTRUE);
3192 if ((found<0.55*foundable) || shared>0.5*found || (track->GetSigmaY2()+track->GetSigmaZ2())>0.5){
3194 seed->~AliTPCseed();
3200 FollowProlongation(*track, i2,1);
3204 track->SetBConstrain(1);
3205 // track->fLastPoint = i1+fInnerSec->GetNRows(); // first cluster in track position
3206 track->SetLastPoint(i1); // first cluster in track position
3207 track->SetFirstPoint(track->GetLastPoint());
3209 if (track->GetNumberOfClusters()<(i1-i2)*0.5 ||
3210 track->GetNumberOfClusters() < track->GetNFoundable()*0.6 ||
3211 track->GetNShared()>0.4*track->GetNumberOfClusters() ) {
3213 seed->~AliTPCseed();
3217 // Z VERTEX CONDITION
3218 Double_t zv, bz=GetBz();
3219 if ( !track->GetZAt(0.,bz,zv) ) continue;
3220 if (TMath::Abs(zv-z3)>cuts[2]) {
3221 FollowProlongation(*track, TMath::Max(i2-20,0));
3222 if ( !track->GetZAt(0.,bz,zv) ) continue;
3223 if (TMath::Abs(zv-z3)>cuts[2]){
3224 FollowProlongation(*track, TMath::Max(i2-40,0));
3225 if ( !track->GetZAt(0.,bz,zv) ) continue;
3226 if (TMath::Abs(zv-z3)>cuts[2] &&(track->GetNumberOfClusters() > track->GetNFoundable()*0.7)){
3227 // make seed without constrain
3228 AliTPCseed * track2 = MakeSeed(track,0.2,0.5,1.);
3229 FollowProlongation(*track2, i2,1);
3230 track2->SetBConstrain(kFALSE);
3231 track2->SetSeedType(1);
3232 arr->AddLast(track2);
3234 seed->~AliTPCseed();
3239 seed->~AliTPCseed();
3246 track->SetSeedType(0);
3247 arr->AddLast(track);
3248 seed = new AliTPCseed;
3250 // don't consider other combinations
3251 if (track->GetNumberOfClusters() > track->GetNFoundable()*0.8)
3257 Info("MakeSeeds3","\nSeeding statistic:\t%d\t%d\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin,nout1,nout2);
3263 void AliTPCtrackerMI::MakeSeeds5(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, Float_t cuts[4],
3268 //-----------------------------------------------------------------
3269 // This function creates track seeds.
3270 //-----------------------------------------------------------------
3271 // cuts[0] - fP4 cut
3272 // cuts[1] - tan(phi) cut
3273 // cuts[2] - zvertex cut
3274 // cuts[3] - fP3 cut
3284 Double_t x[5], c[15];
3286 // make temporary seed
3287 AliTPCseed * seed = new AliTPCseed;
3288 Double_t alpha=fOuterSec->GetAlpha(), shift=fOuterSec->GetAlphaShift();
3289 // Double_t cs=cos(alpha), sn=sin(alpha);
3294 Double_t x1 = GetXrow(i1-1);
3295 const AliTPCtrackerRow& kr1=GetRow(sec,i1-1);
3296 Double_t y1max = GetMaxY(i1-1)-kr1.GetDeadZone()-1.5;
3298 Double_t x1p = GetXrow(i1);
3299 const AliTPCtrackerRow& kr1p=GetRow(sec,i1);
3301 Double_t x1m = GetXrow(i1-2);
3302 const AliTPCtrackerRow& kr1m=GetRow(sec,i1-2);
3305 //last 3 padrow for seeding
3306 AliTPCtrackerRow& kr3 = GetRow((sec+fkNOS)%fkNOS,i1-7);
3307 Double_t x3 = GetXrow(i1-7);
3308 // Double_t y3max= GetMaxY(i1-7)-kr3.fDeadZone-1.5;
3310 AliTPCtrackerRow& kr3p = GetRow((sec+fkNOS)%fkNOS,i1-6);
3311 Double_t x3p = GetXrow(i1-6);
3313 AliTPCtrackerRow& kr3m = GetRow((sec+fkNOS)%fkNOS,i1-8);
3314 Double_t x3m = GetXrow(i1-8);
3319 Int_t im = i1-4; //middle pad row index
3320 Double_t xm = GetXrow(im); // radius of middle pad-row
3321 const AliTPCtrackerRow& krm=GetRow(sec,im); //middle pad -row
3322 // Double_t ymmax = GetMaxY(im)-kr1.fDeadZone-1.5;
3325 Double_t deltax = x1-x3;
3326 Double_t dymax = deltax*cuts[1];
3327 Double_t dzmax = deltax*cuts[3];
3329 // loop over clusters
3330 for (Int_t is=0; is < kr1; is++) {
3332 if (kr1[is]->IsUsed(10)) continue;
3333 Double_t y1=kr1[is]->GetY(), z1=kr1[is]->GetZ();
3335 if (deltay>0 && TMath::Abs(y1max-TMath::Abs(y1))> deltay ) continue; // seed only at the edge
3337 Int_t index1 = TMath::Max(kr3.Find(z1-dzmax)-1,0);
3338 Int_t index2 = TMath::Min(kr3.Find(z1+dzmax)+1,kr3);
3344 for (Int_t js=index1; js < index2; js++) {
3345 const AliTPCclusterMI *kcl = kr3[js];
3346 if (kcl->IsUsed(10)) continue;
3348 // apply angular cuts
3349 if (TMath::Abs(y1-y3)>dymax) continue;
3352 if (TMath::Abs(z1-z3)>dzmax) continue;
3354 Double_t angley = (y1-y3)/(x1-x3);
3355 Double_t anglez = (z1-z3)/(x1-x3);
3357 Double_t erry = TMath::Abs(angley)*(x1-x1m)*0.5+0.5;
3358 Double_t errz = TMath::Abs(anglez)*(x1-x1m)*0.5+0.5;
3360 Double_t yyym = angley*(xm-x1)+y1;
3361 Double_t zzzm = anglez*(xm-x1)+z1;
3363 const AliTPCclusterMI *kcm = krm.FindNearest2(yyym,zzzm,erry,errz,index);
3365 if (kcm->IsUsed(10)) continue;
3367 erry = TMath::Abs(angley)*(x1-x1m)*0.4+0.5;
3368 errz = TMath::Abs(anglez)*(x1-x1m)*0.4+0.5;
3375 // look around first
3376 const AliTPCclusterMI *kc1m = kr1m.FindNearest2(angley*(x1m-x1)+y1,
3382 if (kc1m->IsUsed(10)) used++;
3384 const AliTPCclusterMI *kc1p = kr1p.FindNearest2(angley*(x1p-x1)+y1,
3390 if (kc1p->IsUsed(10)) used++;
3392 if (used>1) continue;
3393 if (found<1) continue;
3397 const AliTPCclusterMI *kc3m = kr3m.FindNearest2(angley*(x3m-x3)+y3,
3403 if (kc3m->IsUsed(10)) used++;
3407 const AliTPCclusterMI *kc3p = kr3p.FindNearest2(angley*(x3p-x3)+y3,
3413 if (kc3p->IsUsed(10)) used++;
3417 if (used>1) continue;
3418 if (found<3) continue;
3428 x[4]=F1(x1,y1,x2,y2,x3,y3);
3429 //if (TMath::Abs(x[4]) >= cuts[0]) continue;
3432 x[2]=F2(x1,y1,x2,y2,x3,y3);
3435 x[3]=F3n(x1,y1,x2,y2,z1,z2,x[4]);
3436 //if (TMath::Abs(x[3]) > cuts[3]) continue;
3440 Double_t sy1=0.1, sz1=0.1;
3441 Double_t sy2=0.1, sz2=0.1;
3442 Double_t sy3=0.1, sy=0.1, sz=0.1;
3444 Double_t f40=(F1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
3445 Double_t f42=(F1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
3446 Double_t f43=(F1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
3447 Double_t f20=(F2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
3448 Double_t f22=(F2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
3449 Double_t f23=(F2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
3451 Double_t f30=(F3(x1,y1+sy,x2,y2,z1,z2)-x[3])/sy;
3452 Double_t f31=(F3(x1,y1,x2,y2,z1+sz,z2)-x[3])/sz;
3453 Double_t f32=(F3(x1,y1,x2,y2+sy,z1,z2)-x[3])/sy;
3454 Double_t f34=(F3(x1,y1,x2,y2,z1,z2+sz)-x[3])/sz;
3458 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3459 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3460 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3461 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3462 c[13]=f30*sy1*f40+f32*sy2*f42;
3463 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3465 // if (!BuildSeed(kr1[is],kcl,0,x1,x2,x3,x,c)) continue;
3467 index=kr1.GetIndex(is);
3468 seed->~AliTPCseed();
3469 AliTPCseed *track=new(seed) AliTPCseed(x1, sec*alpha+shift, x, c, index);
3471 track->SetIsSeeding(kTRUE);
3474 FollowProlongation(*track, i1-7,1);
3475 if (track->GetNumberOfClusters() < track->GetNFoundable()*0.75 ||
3476 track->GetNShared()>0.6*track->GetNumberOfClusters() || ( track->GetSigmaY2()+ track->GetSigmaZ2())>0.6){
3478 seed->~AliTPCseed();
3484 FollowProlongation(*track, i2,1);
3485 track->SetBConstrain(0);
3486 track->SetLastPoint(i1+fInnerSec->GetNRows()); // first cluster in track position
3487 track->SetFirstPoint(track->GetLastPoint());
3489 if (track->GetNumberOfClusters()<(i1-i2)*0.5 ||
3490 track->GetNumberOfClusters()<track->GetNFoundable()*0.7 ||
3491 track->GetNShared()>2. || track->GetChi2()/track->GetNumberOfClusters()>6 || ( track->GetSigmaY2()+ track->GetSigmaZ2())>0.5 ) {
3493 seed->~AliTPCseed();
3498 FollowProlongation(*track, TMath::Max(i2-10,0),1);
3499 AliTPCseed * track2 = MakeSeed(track,0.2,0.5,0.9);
3500 FollowProlongation(*track2, i2,1);
3501 track2->SetBConstrain(kFALSE);
3502 track2->SetSeedType(4);
3503 arr->AddLast(track2);
3505 seed->~AliTPCseed();
3509 //arr->AddLast(track);
3510 //seed = new AliTPCseed;
3516 Info("MakeSeeds5","\nSeeding statiistic:\t%d\t%d\t%d\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin,nout1,nout2,nout3);
3522 //_____________________________________________________________________________
3523 void AliTPCtrackerMI::MakeSeeds2(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, Float_t */*cuts[4]*/,
3524 Float_t deltay, Bool_t /*bconstrain*/) {
3525 //-----------------------------------------------------------------
3526 // This function creates track seeds - without vertex constraint
3527 //-----------------------------------------------------------------
3528 // cuts[0] - fP4 cut - not applied
3529 // cuts[1] - tan(phi) cut
3530 // cuts[2] - zvertex cut - not applied
3531 // cuts[3] - fP3 cut
3541 Double_t alpha=fOuterSec->GetAlpha(), shift=fOuterSec->GetAlphaShift();
3542 // Double_t cs=cos(alpha), sn=sin(alpha);
3543 Int_t row0 = (i1+i2)/2;
3544 Int_t drow = (i1-i2)/2;
3545 const AliTPCtrackerRow& kr0=fSectors[sec][row0];
3546 AliTPCtrackerRow * kr=0;
3548 AliTPCpolyTrack polytrack;
3549 Int_t nclusters=fSectors[sec][row0];
3550 AliTPCseed * seed = new AliTPCseed;
3555 for (Int_t is=0; is < nclusters; is++) { //LOOP over clusters
3557 Int_t nfoundable =0;
3558 for (Int_t iter =1; iter<2; iter++){ //iterations
3559 const AliTPCtrackerRow& krm=fSectors[sec][row0-iter];
3560 const AliTPCtrackerRow& krp=fSectors[sec][row0+iter];
3561 const AliTPCclusterMI * cl= kr0[is];
3563 if (cl->IsUsed(10)) {
3569 Double_t x = kr0.GetX();
3570 // Initialization of the polytrack
3575 Double_t y0= cl->GetY();
3576 Double_t z0= cl->GetZ();
3580 Double_t ymax = fSectors->GetMaxY(row0)-kr0.GetDeadZone()-1.5;
3581 if (deltay>0 && TMath::Abs(ymax-TMath::Abs(y0))> deltay ) continue; // seed only at the edge
3583 erry = (0.5)*cl->GetSigmaY2()/TMath::Sqrt(cl->GetQ())*6;
3584 errz = (0.5)*cl->GetSigmaZ2()/TMath::Sqrt(cl->GetQ())*6;
3585 polytrack.AddPoint(x,y0,z0,erry, errz);
3588 if (cl->IsUsed(10)) sumused++;
3591 Float_t roady = (5*TMath::Sqrt(cl->GetSigmaY2()+0.2)+1.)*iter;
3592 Float_t roadz = (5*TMath::Sqrt(cl->GetSigmaZ2()+0.2)+1.)*iter;
3595 AliTPCclusterMI * cl1 = krm.FindNearest(y0,z0,roady,roadz);
3596 if (cl1 && TMath::Abs(ymax-TMath::Abs(y0))) {
3597 erry = (0.5)*cl1->GetSigmaY2()/TMath::Sqrt(cl1->GetQ())*3;
3598 errz = (0.5)*cl1->GetSigmaZ2()/TMath::Sqrt(cl1->GetQ())*3;
3599 if (cl1->IsUsed(10)) sumused++;
3600 polytrack.AddPoint(x,cl1->GetY(),cl1->GetZ(),erry,errz);
3604 AliTPCclusterMI * cl2 = krp.FindNearest(y0,z0,roady,roadz);
3606 erry = (0.5)*cl2->GetSigmaY2()/TMath::Sqrt(cl2->GetQ())*3;
3607 errz = (0.5)*cl2->GetSigmaZ2()/TMath::Sqrt(cl2->GetQ())*3;
3608 if (cl2->IsUsed(10)) sumused++;
3609 polytrack.AddPoint(x,cl2->GetY(),cl2->GetZ(),erry,errz);
3612 if (sumused>0) continue;
3614 polytrack.UpdateParameters();
3620 nfoundable = polytrack.GetN();
3621 nfound = nfoundable;
3623 for (Int_t ddrow = iter+1; ddrow<drow;ddrow++){
3624 Float_t maxdist = 0.8*(1.+3./(ddrow));
3625 for (Int_t delta = -1;delta<=1;delta+=2){
3626 Int_t row = row0+ddrow*delta;
3627 kr = &(fSectors[sec][row]);
3628 Double_t xn = kr->GetX();
3629 Double_t ymax1 = fSectors->GetMaxY(row)-kr->GetDeadZone()-1.5;
3630 polytrack.GetFitPoint(xn,yn,zn);
3631 if (TMath::Abs(yn)>ymax1) continue;
3633 AliTPCclusterMI * cln = kr->FindNearest(yn,zn,roady,roadz);
3635 Float_t dist = TMath::Sqrt( (yn-cln->GetY())*(yn-cln->GetY())+(zn-cln->GetZ())*(zn-cln->GetZ()));
3638 erry = (dist+0.3)*cln->GetSigmaY2()/TMath::Sqrt(cln->GetQ())*(1.+1./(ddrow));
3639 errz = (dist+0.3)*cln->GetSigmaZ2()/TMath::Sqrt(cln->GetQ())*(1.+1./(ddrow));
3640 if (cln->IsUsed(10)) {
3641 // printf("used\n");
3649 polytrack.AddPoint(xn,cln->GetY(),cln->GetZ(),erry, errz);
3654 if ( (sumused>3) || (sumused>0.5*nfound) || (nfound<0.6*nfoundable)) break;
3655 polytrack.UpdateParameters();
3658 if ( (sumused>3) || (sumused>0.5*nfound)) {
3659 //printf("sumused %d\n",sumused);
3664 polytrack.GetFitDerivation(kr0.GetX(),dy,dz);
3665 AliTPCpolyTrack track2;
3667 polytrack.Refit(track2,0.5+TMath::Abs(dy)*0.3,0.4+TMath::Abs(dz)*0.3);
3668 if (track2.GetN()<0.5*nfoundable) continue;
3671 if ((nfound>0.6*nfoundable) &&( nfoundable>0.4*(i1-i2))) {
3673 // test seed with and without constrain
3674 for (Int_t constrain=0; constrain<=0;constrain++){
3675 // add polytrack candidate
3677 Double_t x[5], c[15];
3678 Double_t x1,x2,x3,y1,y2,y3,z1,z2,z3;
3679 track2.GetBoundaries(x3,x1);
3681 track2.GetFitPoint(x1,y1,z1);
3682 track2.GetFitPoint(x2,y2,z2);
3683 track2.GetFitPoint(x3,y3,z3);
3685 //is track pointing to the vertex ?
3688 polytrack.GetFitPoint(x0,y0,z0);
3701 x[4]=F1(x1,y1,x2,y2,x3,y3);
3703 // if (TMath::Abs(x[4]) >= cuts[0]) continue; //
3704 x[2]=F2(x1,y1,x2,y2,x3,y3);
3706 //if (TMath::Abs(x[4]*x1-x[2]) >= cuts[1]) continue;
3707 //x[3]=F3(x1,y1,x2,y2,z1,z2);
3708 x[3]=F3n(x1,y1,x3,y3,z1,z3,x[4]);
3709 //if (TMath::Abs(x[3]) > cuts[3]) continue;
3712 Double_t sy =0.1, sz =0.1;
3713 Double_t sy1=0.02, sz1=0.02;
3714 Double_t sy2=0.02, sz2=0.02;
3718 sy3=25000*x[4]*x[4]+0.1, sy=0.1, sz=0.1;
3721 Double_t f40=(F1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
3722 Double_t f42=(F1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
3723 Double_t f43=(F1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
3724 Double_t f20=(F2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
3725 Double_t f22=(F2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
3726 Double_t f23=(F2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
3728 Double_t f30=(F3(x1,y1+sy,x3,y3,z1,z3)-x[3])/sy;
3729 Double_t f31=(F3(x1,y1,x3,y3,z1+sz,z3)-x[3])/sz;
3730 Double_t f32=(F3(x1,y1,x3,y3+sy,z1,z3)-x[3])/sy;
3731 Double_t f34=(F3(x1,y1,x3,y3,z1,z3+sz)-x[3])/sz;
3736 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3737 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3738 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3739 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3740 c[13]=f30*sy1*f40+f32*sy2*f42;
3741 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3743 //Int_t row1 = fSectors->GetRowNumber(x1);
3744 Int_t row1 = GetRowNumber(x1);
3748 seed->~AliTPCseed();
3749 AliTPCseed *track=new(seed) AliTPCseed(x1,sec*alpha+shift,x,c,index);
3750 track->SetIsSeeding(kTRUE);
3751 Int_t rc=FollowProlongation(*track, i2);
3752 if (constrain) track->SetBConstrain(1);
3754 track->SetBConstrain(0);
3755 track->SetLastPoint(row1+fInnerSec->GetNRows()); // first cluster in track position
3756 track->SetFirstPoint(track->GetLastPoint());
3758 if (rc==0 || track->GetNumberOfClusters()<(i1-i2)*0.5 ||
3759 track->GetNumberOfClusters() < track->GetNFoundable()*0.6 ||
3760 track->GetNShared()>0.4*track->GetNumberOfClusters()) {
3763 seed->~AliTPCseed();
3766 arr->AddLast(track);
3767 seed = new AliTPCseed;
3771 } // if accepted seed
3774 Info("MakeSeeds2","\nSeeding statiistic:\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin3);
3780 AliTPCseed *AliTPCtrackerMI::MakeSeed(AliTPCseed *const track, Float_t r0, Float_t r1, Float_t r2)
3784 //reseed using track points
3785 Int_t p0 = int(r0*track->GetNumberOfClusters()); // point 0
3786 Int_t p1 = int(r1*track->GetNumberOfClusters());
3787 Int_t p2 = int(r2*track->GetNumberOfClusters()); // last point
3789 Double_t x0[3],x1[3],x2[3];
3790 for (Int_t i=0;i<3;i++){
3796 // find track position at given ratio of the length
3797 Int_t sec0=0, sec1=0, sec2=0;
3800 for (Int_t i=0;i<160;i++){
3801 if (track->GetClusterPointer(i)){
3803 AliTPCTrackerPoint *trpoint =track->GetTrackPoint(i);
3804 if ( (index<p0) || x0[0]<0 ){
3805 if (trpoint->GetX()>1){
3806 clindex = track->GetClusterIndex2(i);
3808 x0[0] = trpoint->GetX();
3809 x0[1] = trpoint->GetY();
3810 x0[2] = trpoint->GetZ();
3811 sec0 = ((clindex&0xff000000)>>24)%18;
3816 if ( (index<p1) &&(trpoint->GetX()>1)){
3817 clindex = track->GetClusterIndex2(i);
3819 x1[0] = trpoint->GetX();
3820 x1[1] = trpoint->GetY();
3821 x1[2] = trpoint->GetZ();
3822 sec1 = ((clindex&0xff000000)>>24)%18;
3825 if ( (index<p2) &&(trpoint->GetX()>1)){
3826 clindex = track->GetClusterIndex2(i);
3828 x2[0] = trpoint->GetX();
3829 x2[1] = trpoint->GetY();
3830 x2[2] = trpoint->GetZ();
3831 sec2 = ((clindex&0xff000000)>>24)%18;
3838 Double_t alpha, cs,sn, xx2,yy2;
3840 alpha = (sec1-sec2)*fSectors->GetAlpha();
3841 cs = TMath::Cos(alpha);
3842 sn = TMath::Sin(alpha);
3843 xx2= x1[0]*cs-x1[1]*sn;
3844 yy2= x1[0]*sn+x1[1]*cs;
3848 alpha = (sec0-sec2)*fSectors->GetAlpha();
3849 cs = TMath::Cos(alpha);
3850 sn = TMath::Sin(alpha);
3851 xx2= x0[0]*cs-x0[1]*sn;
3852 yy2= x0[0]*sn+x0[1]*cs;
3858 Double_t x[5],c[15];
3862 x[4]=F1(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]);
3863 // if (x[4]>1) return 0;
3864 x[2]=F2(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]);
3865 x[3]=F3n(x2[0],x2[1],x0[0],x0[1],x2[2],x0[2],x[4]);
3866 //if (TMath::Abs(x[3]) > 2.2) return 0;
3867 //if (TMath::Abs(x[2]) > 1.99) return 0;
3869 Double_t sy =0.1, sz =0.1;
3871 Double_t sy1=0.02+track->GetSigmaY2(), sz1=0.02+track->GetSigmaZ2();
3872 Double_t sy2=0.01+track->GetSigmaY2(), sz2=0.01+track->GetSigmaZ2();
3873 Double_t sy3=0.01+track->GetSigmaY2();
3875 Double_t f40=(F1(x2[0],x2[1]+sy,x1[0],x1[1],x0[0],x0[1])-x[4])/sy;
3876 Double_t f42=(F1(x2[0],x2[1],x1[0],x1[1]+sy,x0[0],x0[1])-x[4])/sy;
3877 Double_t f43=(F1(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]+sy)-x[4])/sy;
3878 Double_t f20=(F2(x2[0],x2[1]+sy,x1[0],x1[1],x0[0],x0[1])-x[2])/sy;
3879 Double_t f22=(F2(x2[0],x2[1],x1[0],x1[1]+sy,x0[0],x0[1])-x[2])/sy;
3880 Double_t f23=(F2(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]+sy)-x[2])/sy;
3882 Double_t f30=(F3(x2[0],x2[1]+sy,x0[0],x0[1],x2[2],x0[2])-x[3])/sy;
3883 Double_t f31=(F3(x2[0],x2[1],x0[0],x0[1],x2[2]+sz,x0[2])-x[3])/sz;
3884 Double_t f32=(F3(x2[0],x2[1],x0[0],x0[1]+sy,x2[2],x0[2])-x[3])/sy;
3885 Double_t f34=(F3(x2[0],x2[1],x0[0],x0[1],x2[2],x0[2]+sz)-x[3])/sz;
3890 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3891 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3892 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3893 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3894 c[13]=f30*sy1*f40+f32*sy2*f42;
3895 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3897 // Int_t row1 = fSectors->GetRowNumber(x2[0]);
3898 AliTPCseed *seed=new AliTPCseed(x2[0], sec2*fSectors->GetAlpha()+fSectors->GetAlphaShift(), x, c, 0);
3899 // Double_t y0,z0,y1,z1, y2,z2;
3900 //seed->GetProlongation(x0[0],y0,z0);
3901 // seed->GetProlongation(x1[0],y1,z1);
3902 //seed->GetProlongation(x2[0],y2,z2);
3904 seed->SetLastPoint(pp2);
3905 seed->SetFirstPoint(pp2);
3912 AliTPCseed *AliTPCtrackerMI::ReSeed(const AliTPCseed *track, Float_t r0, Float_t r1, Float_t r2)
3916 //reseed using founded clusters
3918 // Find the number of clusters
3919 Int_t nclusters = 0;
3920 for (Int_t irow=0;irow<160;irow++){
3921 if (track->GetClusterIndex(irow)>0) nclusters++;
3925 ipos[0] = TMath::Max(int(r0*nclusters),0); // point 0 cluster
3926 ipos[1] = TMath::Min(int(r1*nclusters),nclusters-1); //
3927 ipos[2] = TMath::Min(int(r2*nclusters),nclusters-1); // last point
3930 Double_t xyz[3][3]={{0}};
3931 Int_t row[3]={0},sec[3]={0,0,0};
3933 // find track row position at given ratio of the length
3935 for (Int_t irow=0;irow<160;irow++){
3936 if (track->GetClusterIndex2(irow)<0) continue;
3938 for (Int_t ipoint=0;ipoint<3;ipoint++){
3939 if (index<=ipos[ipoint]) row[ipoint] = irow;
3943 //Get cluster and sector position
3944 for (Int_t ipoint=0;ipoint<3;ipoint++){
3945 Int_t clindex = track->GetClusterIndex2(row[ipoint]);
3946 AliTPCclusterMI * cl = GetClusterMI(clindex);
3949 // AliTPCclusterMI * cl = GetClusterMI(clindex);
3952 sec[ipoint] = ((clindex&0xff000000)>>24)%18;
3953 xyz[ipoint][0] = GetXrow(row[ipoint]);
3954 xyz[ipoint][1] = cl->GetY();
3955 xyz[ipoint][2] = cl->GetZ();
3959 // Calculate seed state vector and covariance matrix
3961 Double_t alpha, cs,sn, xx2,yy2;
3963 alpha = (sec[1]-sec[2])*fSectors->GetAlpha();
3964 cs = TMath::Cos(alpha);
3965 sn = TMath::Sin(alpha);
3966 xx2= xyz[1][0]*cs-xyz[1][1]*sn;
3967 yy2= xyz[1][0]*sn+xyz[1][1]*cs;
3971 alpha = (sec[0]-sec[2])*fSectors->GetAlpha();
3972 cs = TMath::Cos(alpha);
3973 sn = TMath::Sin(alpha);
3974 xx2= xyz[0][0]*cs-xyz[0][1]*sn;
3975 yy2= xyz[0][0]*sn+xyz[0][1]*cs;
3981 Double_t x[5],c[15];
3985 x[4]=F1(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
3986 x[2]=F2(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
3987 x[3]=F3n(xyz[2][0],xyz[2][1],xyz[0][0],xyz[0][1],xyz[2][2],xyz[0][2],x[4]);
3989 Double_t sy =0.1, sz =0.1;
3991 Double_t sy1=0.2, sz1=0.2;
3992 Double_t sy2=0.2, sz2=0.2;
3995 Double_t f40=(F1(xyz[2][0],xyz[2][1]+sy,xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1])-x[4])/sy;
3996 Double_t f42=(F1(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1]+sy,xyz[0][0],xyz[0][1])-x[4])/sy;
3997 Double_t f43=(F1(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]+sy)-x[4])/sy;
3998 Double_t f20=(F2(xyz[2][0],xyz[2][1]+sy,xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1])-x[2])/sy;
3999 Double_t f22=(F2(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1]+sy,xyz[0][0],xyz[0][1])-x[2])/sy;
4000 Double_t f23=(F2(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]+sy)-x[2])/sy;
4002 Double_t f30=(F3(xyz[2][0],xyz[2][1]+sy,xyz[0][0],xyz[0][1],xyz[2][2],xyz[0][2])-x[3])/sy;
4003 Double_t f31=(F3(xyz[2][0],xyz[2][1],xyz[0][0],xyz[0][1],xyz[2][2]+sz,xyz[0][2])-x[3])/sz;
4004 Double_t f32=(F3(xyz[2][0],xyz[2][1],xyz[0][0],xyz[0][1]+sy,xyz[2][2],xyz[0][2])-x[3])/sy;
4005 Double_t f34=(F3(xyz[2][0],xyz[2][1],xyz[0][0],xyz[0][1],xyz[2][2],xyz[0][2]+sz)-x[3])/sz;
4010 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
4011 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
4012 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
4013 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
4014 c[13]=f30*sy1*f40+f32*sy2*f42;
4015 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
4017 // Int_t row1 = fSectors->GetRowNumber(xyz[2][0]);
4018 AliTPCseed *seed=new AliTPCseed(xyz[2][0], sec[2]*fSectors->GetAlpha()+fSectors->GetAlphaShift(), x, c, 0);
4019 seed->SetLastPoint(row[2]);
4020 seed->SetFirstPoint(row[2]);
4025 AliTPCseed *AliTPCtrackerMI::ReSeed(AliTPCseed *track,Int_t r0, Bool_t forward)
4029 //reseed using founded clusters
4032 Int_t row[3]={0,0,0};
4033 Int_t sec[3]={0,0,0};
4035 // forward direction
4037 for (Int_t irow=r0;irow<160;irow++){
4038 if (track->GetClusterIndex(irow)>0){
4043 for (Int_t irow=160;irow>r0;irow--){
4044 if (track->GetClusterIndex(irow)>0){
4049 for (Int_t irow=row[2]-15;irow>row[0];irow--){
4050 if (track->GetClusterIndex(irow)>0){
4058 for (Int_t irow=0;irow<r0;irow++){
4059 if (track->GetClusterIndex(irow)>0){
4064 for (Int_t irow=r0;irow>0;irow--){
4065 if (track->GetClusterIndex(irow)>0){
4070 for (Int_t irow=row[2]-15;irow>row[0];irow--){
4071 if (track->GetClusterIndex(irow)>0){
4078 if ((row[2]-row[0])<20) return 0;
4079 if (row[1]==0) return 0;
4082 //Get cluster and sector position
4083 for (Int_t ipoint=0;ipoint<3;ipoint++){
4084 Int_t clindex = track->GetClusterIndex2(row[ipoint]);
4085 AliTPCclusterMI * cl = GetClusterMI(clindex);
4088 // AliTPCclusterMI * cl = GetClusterMI(clindex);
4091 sec[ipoint] = ((clindex&0xff000000)>>24)%18;
4092 xyz[ipoint][0] = GetXrow(row[ipoint]);
4093 AliTPCTrackerPoint * point = track->GetTrackPoint(row[ipoint]);
4094 if (point&&ipoint<2){
4096 xyz[ipoint][1] = point->GetY();
4097 xyz[ipoint][2] = point->GetZ();
4100 xyz[ipoint][1] = cl->GetY();
4101 xyz[ipoint][2] = cl->GetZ();
4108 // Calculate seed state vector and covariance matrix
4110 Double_t alpha, cs,sn, xx2,yy2;
4112 alpha = (sec[1]-sec[2])*fSectors->GetAlpha();
4113 cs = TMath::Cos(alpha);
4114 sn = TMath::Sin(alpha);
4115 xx2= xyz[1][0]*cs-xyz[1][1]*sn;
4116 yy2= xyz[1][0]*sn+xyz[1][1]*cs;
4120 alpha = (sec[0]-sec[2])*fSectors->GetAlpha();
4121 cs = TMath::Cos(alpha);
4122 sn = TMath::Sin(alpha);
4123 xx2= xyz[0][0]*cs-xyz[0][1]*sn;
4124 yy2= xyz[0][0]*sn+xyz[0][1]*cs;
4130 Double_t x[5],c[15];
4134 x[4]=F1(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
4135 x[2]=F2(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
4136 x[3]=F3n(xyz[2][0],xyz[2][1],xyz[0][0],xyz[0][1],xyz[2][2],xyz[0][2],x[4]);
4138 Double_t sy =0.1, sz =0.1;
4140 Double_t sy1=0.2, sz1=0.2;
4141 Double_t sy2=0.2, sz2=0.2;
4144 Double_t f40=(F1(xyz[2][0],xyz[2][1]+sy,xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1])-x[4])/sy;
4145 Double_t f42=(F1(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1]+sy,xyz[0][0],xyz[0][1])-x[4])/sy;
4146 Double_t f43=(F1(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]+sy)-x[4])/sy;
4147 Double_t f20=(F2(xyz[2][0],xyz[2][1]+sy,xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1])-x[2])/sy;
4148 Double_t f22=(F2(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1]+sy,xyz[0][0],xyz[0][1])-x[2])/sy;
4149 Double_t f23=(F2(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]+sy)-x[2])/sy;
4151 Double_t f30=(F3(xyz[2][0],xyz[2][1]+sy,xyz[0][0],xyz[0][1],xyz[2][2],xyz[0][2])-x[3])/sy;
4152 Double_t f31=(F3(xyz[2][0],xyz[2][1],xyz[0][0],xyz[0][1],xyz[2][2]+sz,xyz[0][2])-x[3])/sz;
4153 Double_t f32=(F3(xyz[2][0],xyz[2][1],xyz[0][0],xyz[0][1]+sy,xyz[2][2],xyz[0][2])-x[3])/sy;
4154 Double_t f34=(F3(xyz[2][0],xyz[2][1],xyz[0][0],xyz[0][1],xyz[2][2],xyz[0][2]+sz)-x[3])/sz;
4159 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
4160 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
4161 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
4162 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
4163 c[13]=f30*sy1*f40+f32*sy2*f42;
4164 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
4166 // Int_t row1 = fSectors->GetRowNumber(xyz[2][0]);
4167 AliTPCseed *seed=new AliTPCseed(xyz[2][0], sec[2]*fSectors->GetAlpha()+fSectors->GetAlphaShift(), x, c, 0);
4168 seed->SetLastPoint(row[2]);
4169 seed->SetFirstPoint(row[2]);
4170 for (Int_t i=row[0];i<row[2];i++){
4171 seed->SetClusterIndex(i, track->GetClusterIndex(i));
4179 void AliTPCtrackerMI::FindMultiMC(const TObjArray * array, AliESDEvent */*esd*/, Int_t iter)
4182 // find multi tracks - THIS FUNCTION IS ONLY FOR DEBUG PURPOSES
4184 // Use AliTPCReconstructor::StreamLevel()>2 if you want to tune parameters - cuts
4186 // Two reasons to have multiple find tracks
4187 // 1. Curling tracks can be find more than once
4188 // 2. Splitted tracks
4189 // a.) Multiple seeding to increase tracking efficiency - (~ 100% reached)
4190 // b.) Edge effect on the sector boundaries
4193 // Algorithm done in 2 phases - because of CPU consumption
4194 // it is n^2 algorithm - for lead-lead 20000x20000 combination are investigated
4196 // Algorihm for curling tracks sign:
4197 // 1 phase -makes a very rough fast cuts to minimize combinatorics
4198 // a.) opposite sign
4199 // b.) one of the tracks - not pointing to the primary vertex -
4200 // c.) delta tan(theta)
4202 // 2 phase - calculates DCA between tracks - time consument
4207 // General cuts - for splitted tracks and for curling tracks
4209 const Float_t kMaxdPhi = 0.2; // maximal distance in phi
4211 // Curling tracks cuts
4216 Int_t nentries = array->GetEntriesFast();
4217 AliHelix *helixes = new AliHelix[nentries];
4218 Float_t *xm = new Float_t[nentries];
4219 Float_t *dz0 = new Float_t[nentries];
4220 Float_t *dz1 = new Float_t[nentries];
4226 // Find track COG in x direction - point with best defined parameters
4228 for (Int_t i=0;i<nentries;i++){
4229 AliTPCseed* track = (AliTPCseed*)array->At(i);
4230 if (!track) continue;
4231 track->SetCircular(0);
4232 new (&helixes[i]) AliHelix(*track);
4236 track->GetDZ(GetX(),GetY(),GetZ(),GetBz(),dz);
4239 for (Int_t icl=0; icl<160; icl++){
4240 AliTPCclusterMI * cl = track->GetClusterPointer(icl);
4246 if (ncl>0) xm[i]/=Float_t(ncl);
4249 for (Int_t i0=0;i0<nentries;i0++){
4250 AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
4251 if (!track0) continue;
4252 Float_t xc0 = helixes[i0].GetHelix(6);
4253 Float_t yc0 = helixes[i0].GetHelix(7);
4254 Float_t r0 = helixes[i0].GetHelix(8);
4255 Float_t rc0 = TMath::Sqrt(xc0*xc0+yc0*yc0);
4256 Float_t fi0 = TMath::ATan2(yc0,xc0);
4258 for (Int_t i1=i0+1;i1<nentries;i1++){
4259 AliTPCseed * track1 = (AliTPCseed*)array->At(i1);
4260 if (!track1) continue;
4261 Int_t lab0=track0->GetLabel();
4262 Int_t lab1=track1->GetLabel();
4263 if (TMath::Abs(lab0)!=TMath::Abs(lab1)) continue;
4265 Float_t xc1 = helixes[i1].GetHelix(6);
4266 Float_t yc1 = helixes[i1].GetHelix(7);
4267 Float_t r1 = helixes[i1].GetHelix(8);
4268 Float_t rc1 = TMath::Sqrt(xc1*xc1+yc1*yc1);
4269 Float_t fi1 = TMath::ATan2(yc1,xc1);
4271 Float_t dfi = fi0-fi1;
4274 if (dfi>1.5*TMath::Pi()) dfi-=TMath::Pi(); // take care about edge effect
4275 if (dfi<-1.5*TMath::Pi()) dfi+=TMath::Pi(); //
4276 if (TMath::Abs(dfi)>kMaxdPhi&&helixes[i0].GetHelix(4)*helixes[i1].GetHelix(4)<0){
4278 // if short tracks with undefined sign
4279 fi1 = -TMath::ATan2(yc1,-xc1);
4282 Float_t dtheta = TMath::Abs(track0->GetTgl()-track1->GetTgl())<TMath::Abs(track0->GetTgl()+track1->GetTgl())? track0->GetTgl()-track1->GetTgl():track0->GetTgl()+track1->GetTgl();
4285 // debug stream to tune "fast cuts"
4287 Double_t dist[3]; // distance at X
4288 Double_t mdist[3]={0,0,0}; // mean distance X+-40cm
4289 track0->GetDistance(track1,0.5*(xm[i0]+xm[i1])-40.,dist,AliTracker::GetBz());
4290 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[i]);
4291 track0->GetDistance(track1,0.5*(xm[i0]+xm[i1])+40.,dist,AliTracker::GetBz());
4292 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[i]);
4293 track0->GetDistance(track1,0.5*(xm[i0]+xm[i1]),dist,AliTracker::GetBz());
4294 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[i]);
4295 for (Int_t i=0;i<3;i++) mdist[i]*=0.33333;
4299 for (Int_t icl=0; icl<160; icl++){
4300 AliTPCclusterMI * cl0 = track0->GetClusterPointer(icl);
4301 AliTPCclusterMI * cl1 = track1->GetClusterPointer(icl);
4304 if (cl0==cl1) sums++;
4308 if (AliTPCReconstructor::StreamLevel()>5) {
4309 TTreeSRedirector &cstream = *fDebugStreamer;
4314 "Tr0.="<<track0<< // seed0
4315 "Tr1.="<<track1<< // seed1
4316 "h0.="<<&helixes[i0]<<
4317 "h1.="<<&helixes[i1]<<
4319 "sum="<<sum<< //the sum of rows with cl in both
4320 "sums="<<sums<< //the sum of shared clusters
4321 "xm0="<<xm[i0]<< // the center of track
4322 "xm1="<<xm[i1]<< // the x center of track
4323 // General cut variables
4324 "dfi="<<dfi<< // distance in fi angle
4325 "dtheta="<<dtheta<< // distance int theta angle
4331 "dist0="<<dist[0]<< //distance x
4332 "dist1="<<dist[1]<< //distance y
4333 "dist2="<<dist[2]<< //distance z
4334 "mdist0="<<mdist[0]<< //distance x
4335 "mdist1="<<mdist[1]<< //distance y
4336 "mdist2="<<mdist[2]<< //distance z
4352 if (AliTPCReconstructor::StreamLevel()>1) {
4353 AliInfo("Time for curling tracks removal DEBUGGING MC");
4360 void AliTPCtrackerMI::FindSplitted(TObjArray * array, AliESDEvent */*esd*/, Int_t /*iter*/){
4362 // Find Splitted tracks and remove the one with worst quality
4363 // Corresponding debug streamer to tune selections - "Splitted2"
4365 // 0. Sort tracks according quility
4366 // 1. Propagate the tracks to the reference radius
4367 // 2. Double_t loop to select close tracks (only to speed up process)
4368 // 3. Calculate cluster overlap ratio - and remove the track if bigger than a threshold
4369 // 4. Delete temporary parameters
4371 const Double_t xref=GetXrow(63); // reference radius -IROC/OROC boundary
4373 const Double_t kCutP1=10; // delta Z cut 10 cm
4374 const Double_t kCutP2=0.15; // delta snp(fi) cut 0.15
4375 const Double_t kCutP3=0.15; // delta tgl(theta) cut 0.15
4376 const Double_t kCutAlpha=0.15; // delta alpha cut
4377 Int_t firstpoint = 0;
4378 Int_t lastpoint = 160;
4380 Int_t nentries = array->GetEntriesFast();
4381 AliExternalTrackParam *params = new AliExternalTrackParam[nentries];
4387 //0. Sort tracks according quality
4388 //1. Propagate the ext. param to reference radius
4389 Int_t nseed = array->GetEntriesFast();
4390 if (nseed<=0) return;
4391 Float_t * quality = new Float_t[nseed];
4392 Int_t * indexes = new Int_t[nseed];
4393 for (Int_t i=0; i<nseed; i++) {
4394 AliTPCseed *pt=(AliTPCseed*)array->UncheckedAt(i);
4399 pt->UpdatePoints(); //select first last max dens points
4400 Float_t * points = pt->GetPoints();
4401 if (points[3]<0.8) quality[i] =-1;
4402 quality[i] = (points[2]-points[0])+pt->GetNumberOfClusters();
4403 //prefer high momenta tracks if overlaps
4404 quality[i] *= TMath::Sqrt(TMath::Abs(pt->Pt())+0.5);
4406 AliTracker::PropagateTrackToBxByBz(&(params[i]),xref,pt->GetMass(),5.,kTRUE);
4407 AliTracker::PropagateTrackToBxByBz(&(params[i]),xref,pt->GetMass(),1.,kTRUE);
4409 TMath::Sort(nseed,quality,indexes);
4411 // 3. Loop over pair of tracks
4413 for (Int_t i0=0; i0<nseed; i0++) {
4414 Int_t index0=indexes[i0];
4415 if (!(array->UncheckedAt(index0))) continue;
4416 AliTPCseed *s1 = (AliTPCseed*)array->UncheckedAt(index0);
4417 if (!s1->IsActive()) continue;
4418 AliExternalTrackParam &par0=params[index0];
4419 for (Int_t i1=i0+1; i1<nseed; i1++) {
4420 Int_t index1=indexes[i1];
4421 if (!(array->UncheckedAt(index1))) continue;
4422 AliTPCseed *s2 = (AliTPCseed*)array->UncheckedAt(index1);
4423 if (!s2->IsActive()) continue;
4424 if (s2->GetKinkIndexes()[0]!=0)
4425 if (s2->GetKinkIndexes()[0] == -s1->GetKinkIndexes()[0]) continue;
4426 AliExternalTrackParam &par1=params[index1];
4427 if (TMath::Abs(par0.GetParameter()[3]-par1.GetParameter()[3])>kCutP3) continue;
4428 if (TMath::Abs(par0.GetParameter()[1]-par1.GetParameter()[1])>kCutP1) continue;
4429 if (TMath::Abs(par0.GetParameter()[2]-par1.GetParameter()[2])>kCutP2) continue;
4430 Double_t dAlpha= TMath::Abs(par0.GetAlpha()-par1.GetAlpha());
4431 if (dAlpha>TMath::Pi()) dAlpha-=TMath::Pi();
4432 if (TMath::Abs(dAlpha)>kCutAlpha) continue;
4437 Int_t firstShared=lastpoint, lastShared=firstpoint;
4438 Int_t firstRow=lastpoint, lastRow=firstpoint;
4440 for (Int_t i=firstpoint;i<lastpoint;i++){
4441 if (s1->GetClusterIndex2(i)>0) nall0++;
4442 if (s2->GetClusterIndex2(i)>0) nall1++;
4443 if (s1->GetClusterIndex2(i)>0 && s2->GetClusterIndex2(i)>0) {
4444 if (i<firstRow) firstRow=i;
4445 if (i>lastRow) lastRow=i;
4447 if ( (s1->GetClusterIndex2(i))==(s2->GetClusterIndex2(i)) && s1->GetClusterIndex2(i)>0) {
4448 if (i<firstShared) firstShared=i;
4449 if (i>lastShared) lastShared=i;
4453 Double_t ratio0 = Float_t(sumShared)/Float_t(TMath::Min(nall0+1,nall1+1));
4454 Double_t ratio1 = Float_t(sumShared)/Float_t(TMath::Max(nall0+1,nall1+1));
4456 if( AliTPCReconstructor::StreamLevel()>1){
4457 TTreeSRedirector &cstream = *fDebugStreamer;
4458 Int_t n0=s1->GetNumberOfClusters();
4459 Int_t n1=s2->GetNumberOfClusters();
4460 Int_t n0F=s1->GetNFoundable();
4461 Int_t n1F=s2->GetNFoundable();
4462 Int_t lab0=s1->GetLabel();
4463 Int_t lab1=s2->GetLabel();
4465 cstream<<"Splitted2"<<
4466 "iter="<<fIteration<<
4467 "lab0="<<lab0<< // MC label if exist
4468 "lab1="<<lab1<< // MC label if exist
4471 "ratio0="<<ratio0<< // shared ratio
4472 "ratio1="<<ratio1<< // shared ratio
4473 "p0.="<<&par0<< // track parameters
4475 "s0.="<<s1<< // full seed
4477 "n0="<<n0<< // number of clusters track 0
4478 "n1="<<n1<< // number of clusters track 1
4479 "nall0="<<nall0<< // number of clusters track 0
4480 "nall1="<<nall1<< // number of clusters track 1
4481 "n0F="<<n0F<< // number of findable
4482 "n1F="<<n1F<< // number of findable
4483 "shared="<<sumShared<< // number of shared clusters
4484 "firstS="<<firstShared<< // first and the last shared row
4485 "lastS="<<lastShared<<
4486 "firstRow="<<firstRow<< // first and the last row with cluster
4487 "lastRow="<<lastRow<< //
4491 // remove track with lower quality
4493 if (ratio0>AliTPCReconstructor::GetRecoParam()->GetCutSharedClusters(0) ||
4494 ratio1>AliTPCReconstructor::GetRecoParam()->GetCutSharedClusters(1)){
4498 delete array->RemoveAt(index1);
4503 // 4. Delete temporary array
4513 void AliTPCtrackerMI::FindCurling(const TObjArray * array, AliESDEvent */*esd*/, Int_t iter)
4516 // find Curling tracks
4517 // Use AliTPCReconstructor::StreamLevel()>1 if you want to tune parameters - cuts
4520 // Algorithm done in 2 phases - because of CPU consumption
4521 // it is n^2 algorithm - for lead-lead 20000x20000 combination are investigated
4522 // see detal in MC part what can be used to cut
4526 const Float_t kMaxC = 400; // maximal curvature to of the track
4527 const Float_t kMaxdTheta = 0.15; // maximal distance in theta
4528 const Float_t kMaxdPhi = 0.15; // maximal distance in phi
4529 const Float_t kPtRatio = 0.3; // ratio between pt
4530 const Float_t kMinDCAR = 2.; // distance to the primary vertex in r - see cpipe cut
4533 // Curling tracks cuts
4536 const Float_t kMaxDeltaRMax = 40; // distance in outer radius
4537 const Float_t kMaxDeltaRMin = 5.; // distance in lower radius - see cpipe cut
4538 const Float_t kMinAngle = 2.9; // angle between tracks
4539 const Float_t kMaxDist = 5; // biggest distance
4541 // The cuts can be tuned using the "MC information stored in Multi tree ==> see FindMultiMC
4544 TCut csign("csign","Tr0.fP[4]*Tr1.fP[4]<0"); //opposite sign
4545 TCut cmax("cmax","abs(Tr0.GetC())>1/400");
4546 TCut cda("cda","sqrt(dtheta^2+dfi^2)<0.15");
4547 TCut ccratio("ccratio","abs((Tr0.fP[4]+Tr1.fP[4])/(abs(Tr0.fP[4])+abs(Tr1.fP[4])))<0.3");
4548 TCut cpipe("cpipe", "min(abs(r0-rc0),abs(r1-rc1))>5");
4550 TCut cdrmax("cdrmax","abs(abs(rc0+r0)-abs(rc1+r1))<40")
4551 TCut cdrmin("cdrmin","abs(abs(rc0+r0)-abs(rc1+r1))<10")
4553 Multi->Draw("dfi","iter==0"+csign+cmax+cda+ccratio); ~94% of curling tracks fulfill
4554 Multi->Draw("min(abs(r0-rc0),abs(r1-rc1))","iter==0&&abs(lab1)==abs(lab0)"+csign+cmax+cda+ccratio+cpipe+cdrmin+cdrmax); //80%
4556 Curling2->Draw("dfi","iter==0&&abs(lab0)==abs(lab1)"+csign+cmax+cdtheta+cdfi+ccratio)
4562 Int_t nentries = array->GetEntriesFast();
4563 AliHelix *helixes = new AliHelix[nentries];
4564 for (Int_t i=0;i<nentries;i++){
4565 AliTPCseed* track = (AliTPCseed*)array->At(i);
4566 if (!track) continue;
4567 track->SetCircular(0);
4568 new (&helixes[i]) AliHelix(*track);
4574 Double_t phase[2][2]={{0,0},{0,0}},radius[2]={0,0};
4580 for (Int_t i0=0;i0<nentries;i0++){
4581 AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
4582 if (!track0) continue;
4583 if (TMath::Abs(track0->GetC())<1/kMaxC) continue;
4584 Float_t xc0 = helixes[i0].GetHelix(6);
4585 Float_t yc0 = helixes[i0].GetHelix(7);
4586 Float_t r0 = helixes[i0].GetHelix(8);
4587 Float_t rc0 = TMath::Sqrt(xc0*xc0+yc0*yc0);
4588 Float_t fi0 = TMath::ATan2(yc0,xc0);
4590 for (Int_t i1=i0+1;i1<nentries;i1++){
4591 AliTPCseed * track1 = (AliTPCseed*)array->At(i1);
4592 if (!track1) continue;
4593 if (TMath::Abs(track1->GetC())<1/kMaxC) continue;
4594 Float_t xc1 = helixes[i1].GetHelix(6);
4595 Float_t yc1 = helixes[i1].GetHelix(7);
4596 Float_t r1 = helixes[i1].GetHelix(8);
4597 Float_t rc1 = TMath::Sqrt(xc1*xc1+yc1*yc1);
4598 Float_t fi1 = TMath::ATan2(yc1,xc1);
4600 Float_t dfi = fi0-fi1;
4603 if (dfi>1.5*TMath::Pi()) dfi-=TMath::Pi(); // take care about edge effect
4604 if (dfi<-1.5*TMath::Pi()) dfi+=TMath::Pi(); //
4605 Float_t dtheta = TMath::Abs(track0->GetTgl()-track1->GetTgl())<TMath::Abs(track0->GetTgl()+track1->GetTgl())? track0->GetTgl()-track1->GetTgl():track0->GetTgl()+track1->GetTgl();
4609 if (track0->GetBConstrain()&&track1->GetBConstrain()) continue; // not constrained
4610 if (track1->GetSigned1Pt()*track0->GetSigned1Pt()>0) continue; // not the same sign
4611 if ( TMath::Abs(track1->GetTgl()+track0->GetTgl())>kMaxdTheta) continue; //distance in the Theta
4612 if ( TMath::Abs(dfi)>kMaxdPhi) continue; //distance in phi
4613 if ( TMath::Sqrt(dfi*dfi+dtheta*dtheta)>kMaxdPhi) continue; //common angular offset
4615 Float_t pt0 = track0->GetSignedPt();
4616 Float_t pt1 = track1->GetSignedPt();
4617 if ((TMath::Abs(pt0+pt1)/(TMath::Abs(pt0)+TMath::Abs(pt1)))>kPtRatio) continue;
4618 if ((iter==1) && TMath::Abs(TMath::Abs(rc0+r0)-TMath::Abs(rc1+r1))>kMaxDeltaRMax) continue;
4619 if ((iter!=1) &&TMath::Abs(TMath::Abs(rc0-r0)-TMath::Abs(rc1-r1))>kMaxDeltaRMin) continue;
4620 if (TMath::Min(TMath::Abs(rc0-r0),TMath::Abs(rc1-r1))<kMinDCAR) continue;
4623 // Now find closest approach
4627 Int_t npoints = helixes[i0].GetRPHIintersections(helixes[i1], phase, radius,10);
4628 if (npoints==0) continue;
4629 helixes[i0].GetClosestPhases(helixes[i1], phase);
4633 Double_t hangles[3];
4634 helixes[i0].Evaluate(phase[0][0],xyz0);
4635 helixes[i1].Evaluate(phase[0][1],xyz1);
4637 helixes[i0].GetAngle(phase[0][0],helixes[i1],phase[0][1],hangles);
4638 Double_t deltah[2],deltabest;
4639 if (TMath::Abs(hangles[2])<kMinAngle) continue;
4643 helixes[i0].ParabolicDCA(helixes[i1],phase[0][0],phase[0][1],radius[0],deltah[0],2);
4645 helixes[i0].ParabolicDCA(helixes[i1],phase[1][0],phase[1][1],radius[1],deltah[1],2);
4646 if (deltah[1]<deltah[0]) ibest=1;
4648 deltabest = TMath::Sqrt(deltah[ibest]);
4649 helixes[i0].Evaluate(phase[ibest][0],xyz0);
4650 helixes[i1].Evaluate(phase[ibest][1],xyz1);
4651 helixes[i0].GetAngle(phase[ibest][0],helixes[i1],phase[ibest][1],hangles);
4652 Double_t radiusbest = TMath::Sqrt(radius[ibest]);
4654 if (deltabest>kMaxDist) continue;
4655 // if (mindcar+mindcaz<40 && (TMath::Abs(hangles[2])<kMinAngle ||deltabest>3)) continue;
4656 Bool_t sign =kFALSE;
4657 if (hangles[2]>kMinAngle) sign =kTRUE;
4660 // circular[i0] = kTRUE;
4661 // circular[i1] = kTRUE;
4662 if (track0->OneOverPt()<track1->OneOverPt()){
4663 track0->SetCircular(track0->GetCircular()+1);
4664 track1->SetCircular(track1->GetCircular()+2);
4667 track1->SetCircular(track1->GetCircular()+1);
4668 track0->SetCircular(track0->GetCircular()+2);
4671 if (AliTPCReconstructor::StreamLevel()>2){
4673 //debug stream to tune "fine" cuts
4674 Int_t lab0=track0->GetLabel();
4675 Int_t lab1=track1->GetLabel();
4676 TTreeSRedirector &cstream = *fDebugStreamer;
4677 cstream<<"Curling2"<<
4693 "npoints="<<npoints<<
4694 "hangles0="<<hangles[0]<<
4695 "hangles1="<<hangles[1]<<
4696 "hangles2="<<hangles[2]<<
4699 "radius="<<radiusbest<<
4700 "deltabest="<<deltabest<<
4701 "phase0="<<phase[ibest][0]<<
4702 "phase1="<<phase[ibest][1]<<
4710 if (AliTPCReconstructor::StreamLevel()>1) {
4711 AliInfo("Time for curling tracks removal");
4720 void AliTPCtrackerMI::FindKinks(TObjArray * array, AliESDEvent *esd)
4727 TObjArray *kinks= new TObjArray(10000);
4728 // TObjArray *v0s= new TObjArray(10000);
4729 Int_t nentries = array->GetEntriesFast();
4730 AliHelix *helixes = new AliHelix[nentries];
4731 Int_t *sign = new Int_t[nentries];
4732 Int_t *nclusters = new Int_t[nentries];
4733 Float_t *alpha = new Float_t[nentries];
4734 AliKink *kink = new AliKink();
4735 Int_t * usage = new Int_t[nentries];
4736 Float_t *zm = new Float_t[nentries];
4737 Float_t *z0 = new Float_t[nentries];
4738 Float_t *fim = new Float_t[nentries];
4739 Float_t *shared = new Float_t[nentries];
4740 Bool_t *circular = new Bool_t[nentries];
4741 Float_t *dca = new Float_t[nentries];
4742 //const AliESDVertex * primvertex = esd->GetVertex();
4744 // nentries = array->GetEntriesFast();
4749 for (Int_t i=0;i<nentries;i++){
4752 AliTPCseed* track = (AliTPCseed*)array->At(i);
4753 if (!track) continue;
4754 track->SetCircular(0);
4756 track->UpdatePoints();
4757 if (( track->GetPoints()[2]- track->GetPoints()[0])>5 && track->GetPoints()[3]>0.8){
4759 nclusters[i]=track->GetNumberOfClusters();
4760 alpha[i] = track->GetAlpha();
4761 new (&helixes[i]) AliHelix(*track);
4763 helixes[i].Evaluate(0,xyz);
4764 sign[i] = (track->GetC()>0) ? -1:1;
4767 if (track->GetProlongation(x,y,z)){
4769 fim[i] = alpha[i]+TMath::ATan2(y,x);
4772 zm[i] = track->GetZ();
4776 circular[i]= kFALSE;
4777 if (track->GetProlongation(0,y,z)) z0[i] = z;
4778 dca[i] = track->GetD(0,0);
4784 Int_t ncandidates =0;
4787 Double_t phase[2][2]={{0,0},{0,0}},radius[2]={0,0};
4790 // Find circling track
4792 for (Int_t i0=0;i0<nentries;i0++){
4793 AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
4794 if (!track0) continue;
4795 if (track0->GetNumberOfClusters()<40) continue;
4796 if (TMath::Abs(1./track0->GetC())>200) continue;
4797 for (Int_t i1=i0+1;i1<nentries;i1++){
4798 AliTPCseed * track1 = (AliTPCseed*)array->At(i1);
4799 if (!track1) continue;
4800 if (track1->GetNumberOfClusters()<40) continue;
4801 if ( TMath::Abs(track1->GetTgl()+track0->GetTgl())>0.1) continue;
4802 if (track0->GetBConstrain()&&track1->GetBConstrain()) continue;
4803 if (TMath::Abs(1./track1->GetC())>200) continue;
4804 if (track1->GetSigned1Pt()*track0->GetSigned1Pt()>0) continue;
4805 if (track1->GetTgl()*track0->GetTgl()>0) continue;
4806 if (TMath::Max(TMath::Abs(1./track0->GetC()),TMath::Abs(1./track1->GetC()))>190) continue;
4807 if (track0->GetBConstrain()&&track1->OneOverPt()<track0->OneOverPt()) continue; //returning - lower momenta
4808 if (track1->GetBConstrain()&&track0->OneOverPt()<track1->OneOverPt()) continue; //returning - lower momenta
4810 Float_t mindcar = TMath::Min(TMath::Abs(dca[i0]),TMath::Abs(dca[i1]));
4811 if (mindcar<5) continue;
4812 Float_t mindcaz = TMath::Min(TMath::Abs(z0[i0]-GetZ()),TMath::Abs(z0[i1]-GetZ()));
4813 if (mindcaz<5) continue;
4814 if (mindcar+mindcaz<20) continue;
4817 Float_t xc0 = helixes[i0].GetHelix(6);
4818 Float_t yc0 = helixes[i0].GetHelix(7);
4819 Float_t r0 = helixes[i0].GetHelix(8);
4820 Float_t xc1 = helixes[i1].GetHelix(6);
4821 Float_t yc1 = helixes[i1].GetHelix(7);
4822 Float_t r1 = helixes[i1].GetHelix(8);
4824 Float_t rmean = (r0+r1)*0.5;
4825 Float_t delta =TMath::Sqrt((xc1-xc0)*(xc1-xc0)+(yc1-yc0)*(yc1-yc0));
4826 //if (delta>30) continue;
4827 if (delta>rmean*0.25) continue;
4828 if (TMath::Abs(r0-r1)/rmean>0.3) continue;
4830 Int_t npoints = helixes[i0].GetRPHIintersections(helixes[i1], phase, radius,10);
4831 if (npoints==0) continue;
4832 helixes[i0].GetClosestPhases(helixes[i1], phase);
4836 Double_t hangles[3];
4837 helixes[i0].Evaluate(phase[0][0],xyz0);
4838 helixes[i1].Evaluate(phase[0][1],xyz1);
4840 helixes[i0].GetAngle(phase[0][0],helixes[i1],phase[0][1],hangles);
4841 Double_t deltah[2],deltabest;
4842 if (hangles[2]<2.8) continue;
4845 helixes[i0].ParabolicDCA(helixes[i1],phase[0][0],phase[0][1],radius[0],deltah[0],2);
4847 helixes[i0].ParabolicDCA(helixes[i1],phase[1][0],phase[1][1],radius[1],deltah[1],2);
4848 if (deltah[1]<deltah[0]) ibest=1;
4850 deltabest = TMath::Sqrt(deltah[ibest]);
4851 helixes[i0].Evaluate(phase[ibest][0],xyz0);
4852 helixes[i1].Evaluate(phase[ibest][1],xyz1);
4853 helixes[i0].GetAngle(phase[ibest][0],helixes[i1],phase[ibest][1],hangles);
4854 Double_t radiusbest = TMath::Sqrt(radius[ibest]);
4856 if (deltabest>6) continue;
4857 if (mindcar+mindcaz<40 && (hangles[2]<3.12||deltabest>3)) continue;
4858 Bool_t lsign =kFALSE;
4859 if (hangles[2]>3.06) lsign =kTRUE;
4862 circular[i0] = kTRUE;
4863 circular[i1] = kTRUE;
4864 if (track0->OneOverPt()<track1->OneOverPt()){
4865 track0->SetCircular(track0->GetCircular()+1);
4866 track1->SetCircular(track1->GetCircular()+2);
4869 track1->SetCircular(track1->GetCircular()+1);
4870 track0->SetCircular(track0->GetCircular()+2);
4873 if (lsign&&AliTPCReconstructor::StreamLevel()>1){
4875 Int_t lab0=track0->GetLabel();
4876 Int_t lab1=track1->GetLabel();
4877 TTreeSRedirector &cstream = *fDebugStreamer;
4878 cstream<<"Curling"<<
4885 "mindcar="<<mindcar<<
4886 "mindcaz="<<mindcaz<<
4889 "npoints="<<npoints<<
4890 "hangles0="<<hangles[0]<<
4891 "hangles2="<<hangles[2]<<
4896 "radius="<<radiusbest<<
4897 "deltabest="<<deltabest<<
4898 "phase0="<<phase[ibest][0]<<
4899 "phase1="<<phase[ibest][1]<<
4909 for (Int_t i =0;i<nentries;i++){
4910 if (sign[i]==0) continue;
4911 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
4918 Double_t cradius0 = 40*40;
4919 Double_t cradius1 = 270*270;
4922 Double_t cdist3=0.55;
4923 for (Int_t j =i+1;j<nentries;j++){
4925 if (sign[j]*sign[i]<1) continue;
4926 if ( (nclusters[i]+nclusters[j])>200) continue;
4927 if ( (nclusters[i]+nclusters[j])<80) continue;
4928 if ( TMath::Abs(zm[i]-zm[j])>60.) continue;
4929 if ( TMath::Abs(fim[i]-fim[j])>0.6 && TMath::Abs(fim[i]-fim[j])<5.7 ) continue;
4930 //AliTPCseed * track1 = (AliTPCseed*)array->At(j); Double_t phase[2][2],radius[2];
4931 Int_t npoints = helixes[i].GetRPHIintersections(helixes[j], phase, radius,20);
4932 if (npoints<1) continue;
4935 if (radius[0]<cradius0||radius[0]>cradius1) continue;
4938 if ( (radius[0]<cradius0||radius[0]>cradius1) && (radius[1]<cradius0||radius[1]>cradius1) ) continue;
4941 Double_t delta1=10000,delta2=10000;
4942 // cuts on the intersection radius
4943 helixes[i].LinearDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta1);
4944 if (radius[0]<20&&delta1<1) continue; //intersection at vertex
4945 if (radius[0]<10&&delta1<3) continue; //intersection at vertex
4947 helixes[i].LinearDCA(helixes[j],phase[1][0],phase[1][1],radius[1],delta2);
4948 if (radius[1]<20&&delta2<1) continue; //intersection at vertex
4949 if (radius[1]<10&&delta2<3) continue; //intersection at vertex
4952 Double_t distance1 = TMath::Min(delta1,delta2);
4953 if (distance1>cdist1) continue; // cut on DCA linear approximation
4955 npoints = helixes[i].GetRPHIintersections(helixes[j], phase, radius,20);
4956 helixes[i].ParabolicDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta1);
4957 if (radius[0]<20&&delta1<1) continue; //intersection at vertex
4958 if (radius[0]<10&&delta1<3) continue; //intersection at vertex
4961 helixes[i].ParabolicDCA(helixes[j],phase[1][0],phase[1][1],radius[1],delta2);
4962 if (radius[1]<20&&delta2<1) continue; //intersection at vertex
4963 if (radius[1]<10&&delta2<3) continue; //intersection at vertex
4965 distance1 = TMath::Min(delta1,delta2);
4968 rkink = TMath::Sqrt(radius[0]);
4971 rkink = TMath::Sqrt(radius[1]);
4973 if (distance1>cdist2) continue;
4976 AliTPCseed * track1 = (AliTPCseed*)array->At(j);
4979 Int_t row0 = GetRowNumber(rkink);
4980 if (row0<10) continue;
4981 if (row0>150) continue;
4984 Float_t dens00=-1,dens01=-1;
4985 Float_t dens10=-1,dens11=-1;
4987 Int_t found,foundable,ishared;
4988 track0->GetClusterStatistic(0,row0-5, found, foundable,ishared,kFALSE);
4989 if (foundable>5) dens00 = Float_t(found)/Float_t(foundable);
4990 track0->GetClusterStatistic(row0+5,155, found, foundable,ishared,kFALSE);
4991 if (foundable>5) dens01 = Float_t(found)/Float_t(foundable);
4993 track1->GetClusterStatistic(0,row0-5, found, foundable,ishared,kFALSE);
4994 if (foundable>10) dens10 = Float_t(found)/Float_t(foundable);
4995 track1->GetClusterStatistic(row0+5,155, found, foundable,ishared,kFALSE);
4996 if (foundable>10) dens11 = Float_t(found)/Float_t(foundable);
4998 if (dens00<dens10 && dens01<dens11) continue;
4999 if (dens00>dens10 && dens01>dens11) continue;
5000 if (TMath::Max(dens00,dens10)<0.1) continue;
5001 if (TMath::Max(dens01,dens11)<0.3) continue;
5003 if (TMath::Min(dens00,dens10)>0.6) continue;
5004 if (TMath::Min(dens01,dens11)>0.6) continue;
5007 AliTPCseed * ktrack0, *ktrack1;
5016 if (TMath::Abs(ktrack0->GetC())>5) continue; // cut on the curvature for mother particle
5017 AliExternalTrackParam paramm(*ktrack0);
5018 AliExternalTrackParam paramd(*ktrack1);
5019 if (row0>60&&ktrack1->GetReference().GetX()>90.)new (¶md) AliExternalTrackParam(ktrack1->GetReference());
5022 kink->SetMother(paramm);
5023 kink->SetDaughter(paramd);
5026 Float_t x[3] = { kink->GetPosition()[0],kink->GetPosition()[1],kink->GetPosition()[2]};
5028 fkParam->Transform0to1(x,index);
5029 fkParam->Transform1to2(x,index);
5030 row0 = GetRowNumber(x[0]);
5032 if (kink->GetR()<100) continue;
5033 if (kink->GetR()>240) continue;
5034 if (kink->GetPosition()[2]/kink->GetR()>AliTPCReconstructor::GetCtgRange()) continue; //out of fiducial volume
5035 if (kink->GetDistance()>cdist3) continue;
5036 Float_t dird = kink->GetDaughterP()[0]*kink->GetPosition()[0]+kink->GetDaughterP()[1]*kink->GetPosition()[1]; // rough direction estimate
5037 if (dird<0) continue;
5039 Float_t dirm = kink->GetMotherP()[0]*kink->GetPosition()[0]+kink->GetMotherP()[1]*kink->GetPosition()[1]; // rough direction estimate
5040 if (dirm<0) continue;
5041 Float_t mpt = TMath::Sqrt(kink->GetMotherP()[0]*kink->GetMotherP()[0]+kink->GetMotherP()[1]*kink->GetMotherP()[1]);
5042 if (mpt<0.2) continue;
5045 //for high momenta momentum not defined well in first iteration
5046 Double_t qt = TMath::Sin(kink->GetAngle(2))*ktrack1->GetP();
5047 if (qt>0.35) continue;
5050 kink->SetLabel(CookLabel(ktrack0,0.4,0,row0),0);
5051 kink->SetLabel(CookLabel(ktrack1,0.4,row0,160),1);
5053 kink->SetTPCDensity(dens00,0,0);
5054 kink->SetTPCDensity(dens01,0,1);
5055 kink->SetTPCDensity(dens10,1,0);
5056 kink->SetTPCDensity(dens11,1,1);
5057 kink->SetIndex(i,0);
5058 kink->SetIndex(j,1);
5061 kink->SetTPCDensity(dens10,0,0);
5062 kink->SetTPCDensity(dens11,0,1);
5063 kink->SetTPCDensity(dens00,1,0);
5064 kink->SetTPCDensity(dens01,1,1);
5065 kink->SetIndex(j,0);
5066 kink->SetIndex(i,1);
5069 if (mpt<1||kink->GetAngle(2)>0.1){
5070 // angle and densities not defined yet
5071 if (kink->GetTPCDensityFactor()<0.8) continue;
5072 if ((2-kink->GetTPCDensityFactor())*kink->GetDistance() >0.25) continue;
5073 if (kink->GetAngle(2)*ktrack0->GetP()<0.003) continue; //too small angle
5074 if (kink->GetAngle(2)>0.2&&kink->GetTPCDensityFactor()<1.15) continue;
5075 if (kink->GetAngle(2)>0.2&&kink->GetTPCDensity(0,1)>0.05) continue;
5077 Float_t criticalangle = track0->GetSigmaSnp2()+track0->GetSigmaTgl2();
5078 criticalangle+= track1->GetSigmaSnp2()+track1->GetSigmaTgl2();
5079 criticalangle= 3*TMath::Sqrt(criticalangle);
5080 if (criticalangle>0.02) criticalangle=0.02;
5081 if (kink->GetAngle(2)<criticalangle) continue;
5084 Int_t drow = Int_t(2.+0.5/(0.05+kink->GetAngle(2))); // overlap region defined
5085 Float_t shapesum =0;
5087 for ( Int_t row = row0-drow; row<row0+drow;row++){
5088 if (row<0) continue;
5089 if (row>155) continue;
5090 if (ktrack0->GetClusterPointer(row)){
5091 AliTPCTrackerPoint *point =ktrack0->GetTrackPoint(row);
5092 shapesum+=point->GetSigmaY()+point->GetSigmaZ();
5095 if (ktrack1->GetClusterPointer(row)){
5096 AliTPCTrackerPoint *point =ktrack1->GetTrackPoint(row);
5097 shapesum+=point->GetSigmaY()+point->GetSigmaZ();
5102 kink->SetShapeFactor(-1.);
5105 kink->SetShapeFactor(shapesum/sum);
5107 // esd->AddKink(kink);
5109 // kink->SetMother(paramm);
5110 //kink->SetDaughter(paramd);
5112 Double_t chi2P2 = paramm.GetParameter()[2]-paramd.GetParameter()[2];
5114 chi2P2/=paramm.GetCovariance()[5]+paramd.GetCovariance()[5];
5115 Double_t chi2P3 = paramm.GetParameter()[3]-paramd.GetParameter()[3];
5117 chi2P3/=paramm.GetCovariance()[9]+paramd.GetCovariance()[9];
5119 if (AliTPCReconstructor::StreamLevel()>1) {
5120 (*fDebugStreamer)<<"kinkLpt"<<
5128 if ( chi2P2+chi2P3<AliTPCReconstructor::GetRecoParam()->GetKinkAngleCutChi2(0)){
5132 kinks->AddLast(kink);
5138 // sort the kinks according quality - and refit them towards vertex
5140 Int_t nkinks = kinks->GetEntriesFast();
5141 Float_t *quality = new Float_t[nkinks];
5142 Int_t *indexes = new Int_t[nkinks];
5143 AliTPCseed *mothers = new AliTPCseed[nkinks];
5144 AliTPCseed *daughters = new AliTPCseed[nkinks];
5147 for (Int_t i=0;i<nkinks;i++){
5149 AliKink *kinkl = (AliKink*)kinks->At(i);
5151 // refit kinks towards vertex
5153 Int_t index0 = kinkl->GetIndex(0);
5154 Int_t index1 = kinkl->GetIndex(1);
5155 AliTPCseed * ktrack0 = (AliTPCseed*)array->At(index0);
5156 AliTPCseed * ktrack1 = (AliTPCseed*)array->At(index1);
5158 Int_t sumn=ktrack0->GetNumberOfClusters()+ktrack1->GetNumberOfClusters();
5160 // Refit Kink under if too small angle
5162 if (kinkl->GetAngle(2)<0.05){
5163 kinkl->SetTPCRow0(GetRowNumber(kinkl->GetR()));
5164 Int_t row0 = kinkl->GetTPCRow0();
5165 Int_t drow = Int_t(2.+0.5/(0.05+kinkl->GetAngle(2)));
5168 Int_t last = row0-drow;
5169 if (last<40) last=40;
5170 if (last<ktrack0->GetFirstPoint()+25) last = ktrack0->GetFirstPoint()+25;
5171 AliTPCseed* seed0 = ReSeed(ktrack0,last,kFALSE);
5174 Int_t first = row0+drow;
5175 if (first>130) first=130;
5176 if (first>ktrack1->GetLastPoint()-25) first = TMath::Max(ktrack1->GetLastPoint()-25,30);
5177 AliTPCseed* seed1 = ReSeed(ktrack1,first,kTRUE);
5179 if (seed0 && seed1){
5180 kinkl->SetStatus(1,8);
5181 if (RefitKink(*seed0,*seed1,*kinkl)) kinkl->SetStatus(1,9);
5182 row0 = GetRowNumber(kinkl->GetR());
5183 sumn = seed0->GetNumberOfClusters()+seed1->GetNumberOfClusters();
5184 mothers[i] = *seed0;
5185 daughters[i] = *seed1;
5188 delete kinks->RemoveAt(i);
5189 if (seed0) delete seed0;
5190 if (seed1) delete seed1;
5193 if (kinkl->GetDistance()>0.5 || kinkl->GetR()<110 || kinkl->GetR()>240) {
5194 delete kinks->RemoveAt(i);
5195 if (seed0) delete seed0;
5196 if (seed1) delete seed1;
5204 if (kinkl) quality[i] = 160*((0.1+kinkl->GetDistance())*(2.-kinkl->GetTPCDensityFactor()))/(sumn+40.); //the longest -clossest will win
5206 TMath::Sort(nkinks,quality,indexes,kFALSE);
5208 //remove double find kinks
5210 for (Int_t ikink0=1;ikink0<nkinks;ikink0++){
5211 AliKink * kink0 = (AliKink*) kinks->At(indexes[ikink0]);
5212 if (!kink0) continue;
5214 for (Int_t ikink1=0;ikink1<ikink0;ikink1++){
5215 kink0 = (AliKink*) kinks->At(indexes[ikink0]);
5216 if (!kink0) continue;
5217 AliKink * kink1 = (AliKink*) kinks->At(indexes[ikink1]);
5218 if (!kink1) continue;
5219 // if not close kink continue
5220 if (TMath::Abs(kink1->GetPosition()[2]-kink0->GetPosition()[2])>10) continue;
5221 if (TMath::Abs(kink1->GetPosition()[1]-kink0->GetPosition()[1])>10) continue;
5222 if (TMath::Abs(kink1->GetPosition()[0]-kink0->GetPosition()[0])>10) continue;
5224 AliTPCseed &mother0 = mothers[indexes[ikink0]];
5225 AliTPCseed &daughter0 = daughters[indexes[ikink0]];
5226 AliTPCseed &mother1 = mothers[indexes[ikink1]];
5227 AliTPCseed &daughter1 = daughters[indexes[ikink1]];
5228 Int_t row0 = (kink0->GetTPCRow0()+kink1->GetTPCRow0())/2;
5237 for (Int_t i=0;i<row0;i++){
5238 if (mother0.GetClusterIndex(i)>0 && mother1.GetClusterIndex(i)>0){
5241 if (mother0.GetClusterIndex(i)==mother1.GetClusterIndex(i)){
5248 for (Int_t i=row0;i<158;i++){
5249 if (daughter0.GetClusterIndex(i)>0 && daughter0.GetClusterIndex(i)>0){
5252 if (mother0.GetClusterIndex(i)==mother1.GetClusterIndex(i)){
5258 Float_t ratio = Float_t(same+1)/Float_t(both+1);
5259 Float_t ratiom = Float_t(samem+1)/Float_t(bothm+1);
5260 Float_t ratiod = Float_t(samed+1)/Float_t(bothd+1);
5261 if (ratio>0.3 && ratiom>0.5 &&ratiod>0.5) {
5262 Int_t sum0 = mother0.GetNumberOfClusters()+daughter0.GetNumberOfClusters();
5263 Int_t sum1 = mother1.GetNumberOfClusters()+daughter1.GetNumberOfClusters();
5265 shared[kink0->GetIndex(0)]= kTRUE;
5266 shared[kink0->GetIndex(1)]= kTRUE;
5267 delete kinks->RemoveAt(indexes[ikink0]);
5271 shared[kink1->GetIndex(0)]= kTRUE;
5272 shared[kink1->GetIndex(1)]= kTRUE;
5273 delete kinks->RemoveAt(indexes[ikink1]);
5280 for (Int_t i=0;i<nkinks;i++){
5281 AliKink * kinkl = (AliKink*) kinks->At(indexes[i]);
5282 if (!kinkl) continue;
5283 kinkl->SetTPCRow0(GetRowNumber(kinkl->GetR()));
5284 Int_t index0 = kinkl->GetIndex(0);
5285 Int_t index1 = kinkl->GetIndex(1);
5286 if (circular[index0]||(circular[index1]&&kinkl->GetDistance()>0.2)) continue;
5287 kinkl->SetMultiple(usage[index0],0);
5288 kinkl->SetMultiple(usage[index1],1);
5289 if (kinkl->GetMultiple()[0]+kinkl->GetMultiple()[1]>2) continue;
5290 if (kinkl->GetMultiple()[0]+kinkl->GetMultiple()[1]>0 && quality[indexes[i]]>0.2) continue;
5291 if (kinkl->GetMultiple()[0]+kinkl->GetMultiple()[1]>0 && kinkl->GetDistance()>0.2) continue;
5292 if (circular[index0]||(circular[index1]&&kinkl->GetDistance()>0.1)) continue;
5294 AliTPCseed * ktrack0 = (AliTPCseed*)array->At(index0);
5295 AliTPCseed * ktrack1 = (AliTPCseed*)array->At(index1);
5296 if (!ktrack0 || !ktrack1) continue;
5297 Int_t index = esd->AddKink(kinkl);
5300 if ( ktrack0->GetKinkIndex(0)==0 && ktrack1->GetKinkIndex(0)==0) { //best kink
5301 if (mothers[indexes[i]].GetNumberOfClusters()>20 && daughters[indexes[i]].GetNumberOfClusters()>20 && (mothers[indexes[i]].GetNumberOfClusters()+daughters[indexes[i]].GetNumberOfClusters())>100){
5302 *ktrack0 = mothers[indexes[i]];
5303 *ktrack1 = daughters[indexes[i]];
5307 ktrack0->SetKinkIndex(usage[index0],-(index+1));
5308 ktrack1->SetKinkIndex(usage[index1], (index+1));
5313 // Remove tracks corresponding to shared kink's
5315 for (Int_t i=0;i<nentries;i++){
5316 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5317 if (!track0) continue;
5318 if (track0->GetKinkIndex(0)!=0) continue;
5319 if (shared[i]) delete array->RemoveAt(i);
5324 RemoveUsed2(array,0.5,0.4,30);
5326 for (Int_t i=0;i<nentries;i++){
5327 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5328 if (!track0) continue;
5329 track0->CookdEdx(0.02,0.6);
5333 for (Int_t i=0;i<nentries;i++){
5334 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5335 if (!track0) continue;
5336 if (track0->Pt()<1.4) continue;
5337 //remove double high momenta tracks - overlapped with kink candidates
5340 for (Int_t icl=track0->GetFirstPoint();icl<track0->GetLastPoint(); icl++){
5341 if (track0->GetClusterPointer(icl)!=0){
5343 if (track0->GetClusterPointer(icl)->IsUsed(10)) ishared++;
5346 if (Float_t(ishared+1)/Float_t(all+1)>0.5) {
5347 delete array->RemoveAt(i);
5351 if (track0->GetKinkIndex(0)!=0) continue;
5352 if (track0->GetNumberOfClusters()<80) continue;
5354 AliTPCseed *pmother = new AliTPCseed();
5355 AliTPCseed *pdaughter = new AliTPCseed();
5356 AliKink *pkink = new AliKink;
5358 AliTPCseed & mother = *pmother;
5359 AliTPCseed & daughter = *pdaughter;
5360 AliKink & kinkl = *pkink;
5361 if (CheckKinkPoint(track0,mother,daughter, kinkl)){
5362 if (mother.GetNumberOfClusters()<30||daughter.GetNumberOfClusters()<20) {
5366 continue; //too short tracks
5368 if (mother.Pt()<1.4) {
5374 Int_t row0= kinkl.GetTPCRow0();
5375 if (kinkl.GetDistance()>0.5 || kinkl.GetR()<110. || kinkl.GetR()>240.) {
5382 Int_t index = esd->AddKink(&kinkl);
5383 mother.SetKinkIndex(0,-(index+1));
5384 daughter.SetKinkIndex(0,index+1);
5385 if (mother.GetNumberOfClusters()>50) {
5386 delete array->RemoveAt(i);
5387 array->AddAt(new AliTPCseed(mother),i);
5390 array->AddLast(new AliTPCseed(mother));
5392 array->AddLast(new AliTPCseed(daughter));
5393 for (Int_t icl=0;icl<row0;icl++) {
5394 if (mother.GetClusterPointer(icl)) mother.GetClusterPointer(icl)->Use(20);
5397 for (Int_t icl=row0;icl<158;icl++) {
5398 if (daughter.GetClusterPointer(icl)) daughter.GetClusterPointer(icl)->Use(20);
5407 delete [] daughters;
5429 AliInfo(Form("Ncandidates=\t%d\t%d\t%d\t%d\n",esd->GetNumberOfKinks(),ncandidates,ntracks,nall));
5434 Int_t AliTPCtrackerMI::RefitKink(AliTPCseed &mother, AliTPCseed &daughter, const AliESDkink &knk)
5437 // refit kink towards to the vertex
5440 AliKink &kink=(AliKink &)knk;
5442 Int_t row0 = GetRowNumber(kink.GetR());
5443 FollowProlongation(mother,0);
5444 mother.Reset(kFALSE);
5446 FollowProlongation(daughter,row0);
5447 daughter.Reset(kFALSE);
5448 FollowBackProlongation(daughter,158);
5449 daughter.Reset(kFALSE);
5450 Int_t first = TMath::Max(row0-20,30);
5451 Int_t last = TMath::Min(row0+20,140);
5453 const Int_t kNdiv =5;
5454 AliTPCseed param0[kNdiv]; // parameters along the track
5455 AliTPCseed param1[kNdiv]; // parameters along the track
5456 AliKink kinks[kNdiv]; // corresponding kink parameters
5459 for (Int_t irow=0; irow<kNdiv;irow++){
5460 rows[irow] = first +((last-first)*irow)/(kNdiv-1);
5462 // store parameters along the track
5464 for (Int_t irow=0;irow<kNdiv;irow++){
5465 FollowBackProlongation(mother, rows[irow]);
5466 FollowProlongation(daughter,rows[kNdiv-1-irow]);
5467 param0[irow] = mother;
5468 param1[kNdiv-1-irow] = daughter;
5472 for (Int_t irow=0; irow<kNdiv-1;irow++){
5473 if (param0[irow].GetNumberOfClusters()<kNdiv||param1[irow].GetNumberOfClusters()<kNdiv) continue;
5474 kinks[irow].SetMother(param0[irow]);
5475 kinks[irow].SetDaughter(param1[irow]);
5476 kinks[irow].Update();
5479 // choose kink with best "quality"
5481 Double_t mindist = 10000;
5482 for (Int_t irow=0;irow<kNdiv;irow++){
5483 if (param0[irow].GetNumberOfClusters()<20||param1[irow].GetNumberOfClusters()<20) continue;
5484 if (TMath::Abs(kinks[irow].GetR())>240.) continue;
5485 if (TMath::Abs(kinks[irow].GetR())<100.) continue;
5487 Float_t normdist = TMath::Abs(param0[irow].GetX()-kinks[irow].GetR())*(0.1+kink.GetDistance());
5488 normdist/= (param0[irow].GetNumberOfClusters()+param1[irow].GetNumberOfClusters()+40.);
5489 if (normdist < mindist){
5495 if (index==-1) return 0;
5498 param0[index].Reset(kTRUE);
5499 FollowProlongation(param0[index],0);
5501 mother = param0[index];
5502 daughter = param1[index]; // daughter in vertex
5504 kink.SetMother(mother);
5505 kink.SetDaughter(daughter);
5507 kink.SetTPCRow0(GetRowNumber(kink.GetR()));
5508 kink.SetTPCncls(param0[index].GetNumberOfClusters(),0);
5509 kink.SetTPCncls(param1[index].GetNumberOfClusters(),1);
5510 kink.SetLabel(CookLabel(&mother,0.4, 0,kink.GetTPCRow0()),0);
5511 kink.SetLabel(CookLabel(&daughter,0.4, kink.GetTPCRow0(),160),1);
5512 mother.SetLabel(kink.GetLabel(0));
5513 daughter.SetLabel(kink.GetLabel(1));
5519 void AliTPCtrackerMI::UpdateKinkQualityM(AliTPCseed * seed){
5521 // update Kink quality information for mother after back propagation
5523 if (seed->GetKinkIndex(0)>=0) return;
5524 for (Int_t ikink=0;ikink<3;ikink++){
5525 Int_t index = seed->GetKinkIndex(ikink);
5526 if (index>=0) break;
5527 index = TMath::Abs(index)-1;
5528 AliESDkink * kink = fEvent->GetKink(index);
5529 kink->SetTPCDensity(-1,0,0);
5530 kink->SetTPCDensity(1,0,1);
5532 Int_t row0 = kink->GetTPCRow0() - 2 - Int_t( 0.5/ (0.05+kink->GetAngle(2)));
5533 if (row0<15) row0=15;
5535 Int_t row1 = kink->GetTPCRow0() + 2 + Int_t( 0.5/ (0.05+kink->GetAngle(2)));
5536 if (row1>145) row1=145;
5538 Int_t found,foundable,shared;
5539 seed->GetClusterStatistic(0,row0, found, foundable,shared,kFALSE);
5540 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),0,0);
5541 seed->GetClusterStatistic(row1,155, found, foundable,shared,kFALSE);
5542 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),0,1);
5547 void AliTPCtrackerMI::UpdateKinkQualityD(AliTPCseed * seed){
5549 // update Kink quality information for daughter after refit
5551 if (seed->GetKinkIndex(0)<=0) return;
5552 for (Int_t ikink=0;ikink<3;ikink++){
5553 Int_t index = seed->GetKinkIndex(ikink);
5554 if (index<=0) break;
5555 index = TMath::Abs(index)-1;
5556 AliESDkink * kink = fEvent->GetKink(index);
5557 kink->SetTPCDensity(-1,1,0);
5558 kink->SetTPCDensity(-1,1,1);
5560 Int_t row0 = kink->GetTPCRow0() -2 - Int_t( 0.5/ (0.05+kink->GetAngle(2)));
5561 if (row0<15) row0=15;
5563 Int_t row1 = kink->GetTPCRow0() +2 + Int_t( 0.5/ (0.05+kink->GetAngle(2)));
5564 if (row1>145) row1=145;
5566 Int_t found,foundable,shared;
5567 seed->GetClusterStatistic(0,row0, found, foundable,shared,kFALSE);
5568 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),1,0);
5569 seed->GetClusterStatistic(row1,155, found, foundable,shared,kFALSE);
5570 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),1,1);
5576 Int_t AliTPCtrackerMI::CheckKinkPoint(AliTPCseed*seed,AliTPCseed &mother, AliTPCseed &daughter, const AliESDkink &knk)
5579 // check kink point for given track
5580 // if return value=0 kink point not found
5581 // otherwise seed0 correspond to mother particle
5582 // seed1 correspond to daughter particle
5583 // kink parameter of kink point
5584 AliKink &kink=(AliKink &)knk;
5586 Int_t middlerow = (seed->GetFirstPoint()+seed->GetLastPoint())/2;
5587 Int_t first = seed->GetFirstPoint();
5588 Int_t last = seed->GetLastPoint();
5589 if (last-first<20) return 0; // shortest length - 2*30 = 60 pad-rows
5592 AliTPCseed *seed1 = ReSeed(seed,middlerow+20, kTRUE); //middle of chamber
5593 if (!seed1) return 0;
5594 FollowProlongation(*seed1,seed->GetLastPoint()-20);
5595 seed1->Reset(kTRUE);
5596 FollowProlongation(*seed1,158);
5597 seed1->Reset(kTRUE);
5598 last = seed1->GetLastPoint();
5600 AliTPCseed *seed0 = new AliTPCseed(*seed);
5601 seed0->Reset(kFALSE);
5604 AliTPCseed param0[20]; // parameters along the track
5605 AliTPCseed param1[20]; // parameters along the track
5606 AliKink kinks[20]; // corresponding kink parameters
5608 for (Int_t irow=0; irow<20;irow++){
5609 rows[irow] = first +((last-first)*irow)/19;
5611 // store parameters along the track
5613 for (Int_t irow=0;irow<20;irow++){
5614 FollowBackProlongation(*seed0, rows[irow]);
5615 FollowProlongation(*seed1,rows[19-irow]);
5616 param0[irow] = *seed0;
5617 param1[19-irow] = *seed1;
5621 for (Int_t irow=0; irow<19;irow++){
5622 kinks[irow].SetMother(param0[irow]);
5623 kinks[irow].SetDaughter(param1[irow]);
5624 kinks[irow].Update();
5627 // choose kink with biggest change of angle
5629 Double_t maxchange= 0;
5630 for (Int_t irow=1;irow<19;irow++){
5631 if (TMath::Abs(kinks[irow].GetR())>240.) continue;
5632 if (TMath::Abs(kinks[irow].GetR())<110.) continue;
5633 Float_t quality = TMath::Abs(kinks[irow].GetAngle(2))/(3.+TMath::Abs(kinks[irow].GetR()-param0[irow].GetX()));
5634 if ( quality > maxchange){
5635 maxchange = quality;
5642 if (index<0) return 0;
5644 Int_t row0 = GetRowNumber(kinks[index].GetR()); //row 0 estimate
5645 seed0 = new AliTPCseed(param0[index]);
5646 seed1 = new AliTPCseed(param1[index]);
5647 seed0->Reset(kFALSE);
5648 seed1->Reset(kFALSE);
5649 seed0->ResetCovariance(10.);
5650 seed1->ResetCovariance(10.);
5651 FollowProlongation(*seed0,0);
5652 FollowBackProlongation(*seed1,158);
5653 mother = *seed0; // backup mother at position 0
5654 seed0->Reset(kFALSE);
5655 seed1->Reset(kFALSE);
5656 seed0->ResetCovariance(10.);
5657 seed1->ResetCovariance(10.);
5659 first = TMath::Max(row0-20,0);
5660 last = TMath::Min(row0+20,158);
5662 for (Int_t irow=0; irow<20;irow++){
5663 rows[irow] = first +((last-first)*irow)/19;
5665 // store parameters along the track
5667 for (Int_t irow=0;irow<20;irow++){
5668 FollowBackProlongation(*seed0, rows[irow]);
5669 FollowProlongation(*seed1,rows[19-irow]);
5670 param0[irow] = *seed0;
5671 param1[19-irow] = *seed1;
5675 for (Int_t irow=0; irow<19;irow++){
5676 kinks[irow].SetMother(param0[irow]);
5677 kinks[irow].SetDaughter(param1[irow]);
5678 // param0[irow].Dump();
5679 //param1[irow].Dump();
5680 kinks[irow].Update();
5683 // choose kink with biggest change of angle
5686 for (Int_t irow=0;irow<20;irow++){
5687 if (TMath::Abs(kinks[irow].GetR())>250.) continue;
5688 if (TMath::Abs(kinks[irow].GetR())<90.) continue;
5689 Float_t quality = TMath::Abs(kinks[irow].GetAngle(2))/(3.+TMath::Abs(kinks[irow].GetR()-param0[irow].GetX()));
5690 if ( quality > maxchange){
5691 maxchange = quality;
5698 if (index==-1 || param0[index].GetNumberOfClusters()+param1[index].GetNumberOfClusters()<100){
5704 // Float_t anglesigma = TMath::Sqrt(param0[index].fC22+param0[index].fC33+param1[index].fC22+param1[index].fC33);
5706 kink.SetMother(param0[index]);
5707 kink.SetDaughter(param1[index]);
5710 Double_t chi2P2 = param0[index].GetParameter()[2]-param1[index].GetParameter()[2];
5712 chi2P2/=param0[index].GetCovariance()[5]+param1[index].GetCovariance()[5];
5713 Double_t chi2P3 = param0[index].GetParameter()[3]-param1[index].GetParameter()[3];
5715 chi2P3/=param0[index].GetCovariance()[9]+param1[index].GetCovariance()[9];
5717 if (AliTPCReconstructor::StreamLevel()>1) {
5718 (*fDebugStreamer)<<"kinkHpt"<<
5721 "p0.="<<¶m0[index]<<
5722 "p1.="<<¶m1[index]<<
5726 if ( chi2P2+chi2P3<AliTPCReconstructor::GetRecoParam()->GetKinkAngleCutChi2(0)){
5733 row0 = GetRowNumber(kink.GetR());
5734 kink.SetTPCRow0(row0);
5735 kink.SetLabel(CookLabel(seed0,0.5,0,row0),0);
5736 kink.SetLabel(CookLabel(seed1,0.5,row0,158),1);
5737 kink.SetIndex(-10,0);
5738 kink.SetIndex(int(param0[index].GetNumberOfClusters()+param1[index].GetNumberOfClusters()),1);
5739 kink.SetTPCncls(param0[index].GetNumberOfClusters(),0);
5740 kink.SetTPCncls(param1[index].GetNumberOfClusters(),1);
5743 // new (&mother) AliTPCseed(param0[index]);
5744 daughter = param1[index];
5745 daughter.SetLabel(kink.GetLabel(1));
5746 param0[index].Reset(kTRUE);
5747 FollowProlongation(param0[index],0);
5748 mother = param0[index];
5749 mother.SetLabel(kink.GetLabel(0));
5750 if ( chi2P2+chi2P3<AliTPCReconstructor::GetRecoParam()->GetKinkAngleCutChi2(1)){
5762 AliTPCseed* AliTPCtrackerMI::ReSeed(AliTPCseed *t)
5765 // reseed - refit - track
5768 // Int_t last = fSectors->GetNRows()-1;
5770 if (fSectors == fOuterSec){
5771 first = TMath::Max(first, t->GetFirstPoint()-fInnerSec->GetNRows());
5775 first = t->GetFirstPoint();
5777 AliTPCseed * seed = MakeSeed(t,0.1,0.5,0.9);
5778 FollowBackProlongation(*t,fSectors->GetNRows()-1);
5780 FollowProlongation(*t,first);
5790 //_____________________________________________________________________________
5791 Int_t AliTPCtrackerMI::ReadSeeds(const TFile *inp) {
5792 //-----------------------------------------------------------------
5793 // This function reades track seeds.
5794 //-----------------------------------------------------------------
5795 TDirectory *savedir=gDirectory;
5797 TFile *in=(TFile*)inp;
5798 if (!in->IsOpen()) {
5799 cerr<<"AliTPCtrackerMI::ReadSeeds(): input file is not open !\n";
5804 TTree *seedTree=(TTree*)in->Get("Seeds");
5806 cerr<<"AliTPCtrackerMI::ReadSeeds(): ";
5807 cerr<<"can't get a tree with track seeds !\n";
5810 AliTPCtrack *seed=new AliTPCtrack;
5811 seedTree->SetBranchAddress("tracks",&seed);
5813 if (fSeeds==0) fSeeds=new TObjArray(15000);
5815 Int_t n=(Int_t)seedTree->GetEntries();
5816 for (Int_t i=0; i<n; i++) {
5817 seedTree->GetEvent(i);
5818 fSeeds->AddLast(new AliTPCseed(*seed/*,seed->GetAlpha()*/));
5827 Int_t AliTPCtrackerMI::Clusters2Tracks (AliESDEvent *const esd)
5830 if (fSeeds) DeleteSeeds();
5833 if (!fSeeds) return 1;
5835 if (AliTPCReconstructor::StreamLevel()>3) DumpClusters(0,fSeeds);
5841 //_____________________________________________________________________________
5842 Int_t AliTPCtrackerMI::Clusters2Tracks() {
5843 //-----------------------------------------------------------------
5844 // This is a track finder.
5845 //-----------------------------------------------------------------
5846 TDirectory *savedir=gDirectory;
5850 fSeeds = Tracking();
5853 Info("Clusters2Tracks","Time for tracking: \t");timer.Print();timer.Start();
5855 //activate again some tracks
5856 for (Int_t i=0; i<fSeeds->GetEntriesFast(); i++) {
5857 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
5859 Int_t nc=t.GetNumberOfClusters();
5861 delete fSeeds->RemoveAt(i);
5865 if (pt->GetRemoval()==10) {
5866 if (pt->GetDensityFirst(20)>0.8 || pt->GetDensityFirst(30)>0.8 || pt->GetDensityFirst(40)>0.7)
5867 pt->Desactivate(10); // make track again active
5869 pt->Desactivate(20);
5870 delete fSeeds->RemoveAt(i);
5875 RemoveUsed2(fSeeds,0.85,0.85,0);
5876 if (AliTPCReconstructor::GetRecoParam()->GetDoKinks()) FindKinks(fSeeds,fEvent);
5877 //FindCurling(fSeeds, fEvent,0);
5878 if (AliTPCReconstructor::StreamLevel()>5) FindMultiMC(fSeeds, fEvent,-1); // find multi found tracks
5879 RemoveUsed2(fSeeds,0.5,0.4,20);
5880 FindSplitted(fSeeds, fEvent,0); // find multi found tracks
5881 if (AliTPCReconstructor::StreamLevel()>5) FindMultiMC(fSeeds, fEvent,0); // find multi found tracks
5884 // // refit short tracks
5886 Int_t nseed=fSeeds->GetEntriesFast();
5889 for (Int_t i=0; i<nseed; i++) {
5890 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
5892 Int_t nc=t.GetNumberOfClusters();
5894 delete fSeeds->RemoveAt(i);
5897 CookLabel(pt,0.1); //For comparison only
5898 //if ((pt->IsActive() || (pt->fRemoval==10) )&& nc>50 &&pt->GetNumberOfClusters()>0.4*pt->fNFoundable){
5899 if ((pt->IsActive() || (pt->GetRemoval()==10) )){
5901 if (fDebug>0) cerr<<found<<'\r';
5905 delete fSeeds->RemoveAt(i);
5909 //RemoveOverlap(fSeeds,0.99,7,kTRUE);
5911 //RemoveUsed(fSeeds,0.9,0.9,6);
5913 nseed=fSeeds->GetEntriesFast();
5915 for (Int_t i=0; i<nseed; i++) {
5916 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
5918 Int_t nc=t.GetNumberOfClusters();
5920 delete fSeeds->RemoveAt(i);
5924 t.CookdEdx(0.02,0.6);
5925 // CheckKinkPoint(&t,0.05);
5926 //if ((pt->IsActive() || (pt->fRemoval==10) )&& nc>50 &&pt->GetNumberOfClusters()>0.4*pt->fNFoundable){
5927 if ((pt->IsActive() || (pt->GetRemoval()==10) )){
5935 delete fSeeds->RemoveAt(i);
5936 //AliTPCseed * seed1 = ReSeed(pt,0.05,0.5,1);
5938 // FollowProlongation(*seed1,0);
5939 // Int_t n = seed1->GetNumberOfClusters();
5940 // printf("fP4\t%f\t%f\n",seed1->GetC(),pt->GetC());
5941 // printf("fN\t%d\t%d\n", seed1->GetNumberOfClusters(),pt->GetNumberOfClusters());
5944 //AliTPCseed * seed2 = ReSeed(pt,0.95,0.5,0.05);
5948 SortTracks(fSeeds, 1);
5952 PrepareForBackProlongation(fSeeds,5.);
5953 PropagateBack(fSeeds);
5954 printf("Time for back propagation: \t");timer.Print();timer.Start();
5958 PrepareForProlongation(fSeeds,5.);
5959 PropagateForward2(fSeeds);
5961 printf("Time for FORWARD propagation: \t");timer.Print();timer.Start();
5962 // RemoveUsed(fSeeds,0.7,0.7,6);
5963 //RemoveOverlap(fSeeds,0.9,7,kTRUE);
5965 nseed=fSeeds->GetEntriesFast();
5967 for (Int_t i=0; i<nseed; i++) {
5968 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
5970 Int_t nc=t.GetNumberOfClusters();
5972 delete fSeeds->RemoveAt(i);
5975 t.CookdEdx(0.02,0.6);
5976 // CookLabel(pt,0.1); //For comparison only
5977 //if ((pt->IsActive() || (pt->fRemoval==10) )&& nc>50 &&pt->GetNumberOfClusters()>0.4*pt->fNFoundable){
5978 if ((pt->IsActive() || (pt->fRemoval==10) )){
5979 cerr<<found++<<'\r';
5982 delete fSeeds->RemoveAt(i);
5987 // fNTracks = found;
5989 Info("Clusters2Tracks","Time for overlap removal, track writing and dedx cooking: \t"); timer.Print();timer.Start();
5992 // cerr<<"Number of found tracks : "<<"\t"<<found<<endl;
5993 Info("Clusters2Tracks","Number of found tracks %d",found);
5995 // UnloadClusters();
6000 void AliTPCtrackerMI::Tracking(TObjArray * arr)
6003 // tracking of the seeds
6006 fSectors = fOuterSec;
6007 ParallelTracking(arr,150,63);
6008 fSectors = fOuterSec;
6009 ParallelTracking(arr,63,0);
6012 TObjArray * AliTPCtrackerMI::Tracking(Int_t seedtype, Int_t i1, Int_t i2, Float_t cuts[4], Float_t dy, Int_t dsec)
6017 TObjArray * arr = new TObjArray;
6019 fSectors = fOuterSec;
6022 for (Int_t sec=0;sec<fkNOS;sec++){
6023 if (seedtype==3) MakeSeeds3(arr,sec,i1,i2,cuts,dy, dsec);
6024 if (seedtype==4) MakeSeeds5(arr,sec,i1,i2,cuts,dy);
6025 if (seedtype==2) MakeSeeds2(arr,sec,i1,i2,cuts,dy);
6028 Info("Tracking","\nSeeding - %d\t%d\t%d\t%d\n",seedtype,i1,i2,arr->GetEntriesFast());
6040 TObjArray * AliTPCtrackerMI::Tracking()
6044 if (AliTPCReconstructor::GetRecoParam()->GetSpecialSeeding()) return TrackingSpecial();
6047 Int_t nup=fOuterSec->GetNRows()+fInnerSec->GetNRows();
6049 TObjArray * seeds = new TObjArray;
6051 Int_t fLastSeedRowSec=AliTPCReconstructor::GetRecoParam()->GetLastSeedRowSec();
6052 Int_t gapPrim = AliTPCReconstructor::GetRecoParam()->GetSeedGapPrim();
6053 Int_t gapSec = AliTPCReconstructor::GetRecoParam()->GetSeedGapSec();
6061 Float_t fnumber = 3.0;
6062 Float_t fdensity = 3.0;
6067 for (Int_t delta = 0; delta<18; delta+=gapPrim){
6071 arr = Tracking(3,nup-1-delta,nup-1-delta-gap,cuts,-1,1);
6072 SumTracks(seeds,arr);
6073 SignClusters(seeds,fnumber,fdensity);
6075 for (Int_t i=2;i<6;i+=2){
6076 // seed high pt tracks
6079 arr = Tracking(3,nup-i-delta,nup-i-delta-gap,cuts,-1,0);
6080 SumTracks(seeds,arr);
6081 SignClusters(seeds,fnumber,fdensity);
6086 // RemoveUsed(seeds,0.9,0.9,1);
6087 // UnsignClusters();
6088 // SignClusters(seeds,fnumber,fdensity);
6092 for (Int_t delta = 20; delta<120; delta+=gapPrim){
6094 // seed high pt tracks
6098 arr = Tracking(3,nup-delta,nup-delta-gap,cuts,-1);
6099 SumTracks(seeds,arr);
6100 SignClusters(seeds,fnumber,fdensity);
6105 arr = Tracking(3,nup-delta-5,nup-delta-5-gap,cuts,-1);
6106 SumTracks(seeds,arr);
6107 SignClusters(seeds,fnumber,fdensity);
6118 Info("Tracking()","\n\nPrimary seeding\t%d\n\n",seeds->GetEntriesFast());
6122 // RemoveUsed(seeds,0.75,0.75,1);
6124 //SignClusters(seeds,fnumber,fdensity);
6133 arr = Tracking(4,nup-1,nup-1-gap,cuts,-1);
6134 SumTracks(seeds,arr);
6135 SignClusters(seeds,fnumber,fdensity);
6137 arr = Tracking(4,nup-2,nup-2-gap,cuts,-1);
6138 SumTracks(seeds,arr);
6139 SignClusters(seeds,fnumber,fdensity);
6141 arr = Tracking(4,nup-3,nup-3-gap,cuts,-1);
6142 SumTracks(seeds,arr);
6143 SignClusters(seeds,fnumber,fdensity);
6145 arr = Tracking(4,nup-5,nup-5-gap,cuts,-1);
6146 SumTracks(seeds,arr);
6147 SignClusters(seeds,fnumber,fdensity);
6149 arr = Tracking(4,nup-7,nup-7-gap,cuts,-1);
6150 SumTracks(seeds,arr);
6151 SignClusters(seeds,fnumber,fdensity);
6154 arr = Tracking(4,nup-9,nup-9-gap,cuts,-1);
6155 SumTracks(seeds,arr);
6156 SignClusters(seeds,fnumber,fdensity);
6160 for (Int_t delta = 9; delta<30; delta+=gapSec){
6166 arr = Tracking(4,nup-1-delta,nup-1-delta-gap,cuts,-1);
6167 SumTracks(seeds,arr);
6168 SignClusters(seeds,fnumber,fdensity);
6170 arr = Tracking(4,nup-3-delta,nup-5-delta-gap,cuts,4);
6171 SumTracks(seeds,arr);
6172 SignClusters(seeds,fnumber,fdensity);
6185 for (Int_t delta = 30; delta<fLastSeedRowSec; delta+=gapSec){
6191 arr = Tracking(4,nup-1-delta,nup-1-delta-gap,cuts,-1);
6192 SumTracks(seeds,arr);
6193 SignClusters(seeds,fnumber,fdensity);
6195 arr = Tracking(4,nup-5-delta,nup-5-delta-gap,cuts,5 );
6196 SumTracks(seeds,arr);
6197 SignClusters(seeds,fnumber,fdensity);
6201 Info("Tracking()","\n\nSecondary seeding\t%d\n\n",seeds->GetEntriesFast());
6212 TObjArray * AliTPCtrackerMI::TrackingSpecial()
6215 // seeding adjusted for laser and cosmic tests - short tracks with big inclination angle
6216 // no primary vertex seeding tried
6220 Int_t nup=fOuterSec->GetNRows()+fInnerSec->GetNRows();
6222 TObjArray * seeds = new TObjArray;
6227 Float_t fnumber = 3.0;
6228 Float_t fdensity = 3.0;
6231 cuts[0] = AliTPCReconstructor::GetRecoParam()->GetMaxC(); // max curvature
6232 cuts[1] = 3.5; // max tan(phi) angle for seeding
6233 cuts[2] = 3.; // not used (cut on z primary vertex)
6234 cuts[3] = 3.5; // max tan(theta) angle for seeding
6236 for (Int_t delta = 0; nup-delta-gap-1>0; delta+=3){
6238 arr = Tracking(4,nup-1-delta,nup-1-delta-gap,cuts,-1);
6239 SumTracks(seeds,arr);
6240 SignClusters(seeds,fnumber,fdensity);
6244 Info("Tracking()","\n\nSecondary seeding\t%d\n\n",seeds->GetEntriesFast());
6255 void AliTPCtrackerMI::SumTracks(TObjArray *arr1,TObjArray *&arr2) const
6258 //sum tracks to common container
6259 //remove suspicious tracks
6260 Int_t nseed = arr2->GetEntriesFast();
6261 for (Int_t i=0;i<nseed;i++){
6262 AliTPCseed *pt=(AliTPCseed*)arr2->UncheckedAt(i);
6265 // remove tracks with too big curvature
6267 if (TMath::Abs(pt->GetC())>AliTPCReconstructor::GetRecoParam()->GetMaxC()){
6268 delete arr2->RemoveAt(i);
6271 // REMOVE VERY SHORT TRACKS
6272 if (pt->GetNumberOfClusters()<20){
6273 delete arr2->RemoveAt(i);
6276 // NORMAL ACTIVE TRACK
6277 if (pt->IsActive()){
6278 arr1->AddLast(arr2->RemoveAt(i));
6281 //remove not usable tracks
6282 if (pt->GetRemoval()!=10){
6283 delete arr2->RemoveAt(i);
6287 // ENABLE ONLY ENOUGH GOOD STOPPED TRACKS
6288 if (pt->GetDensityFirst(20)>0.8 || pt->GetDensityFirst(30)>0.8 || pt->GetDensityFirst(40)>0.7)
6289 arr1->AddLast(arr2->RemoveAt(i));
6291 delete arr2->RemoveAt(i);
6295 delete arr2; arr2 = 0;
6300 void AliTPCtrackerMI::ParallelTracking(TObjArray *const arr, Int_t rfirst, Int_t rlast)
6303 // try to track in parralel
6305 Int_t nseed=arr->GetEntriesFast();
6306 //prepare seeds for tracking
6307 for (Int_t i=0; i<nseed; i++) {
6308 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i), &t=*pt;
6310 if (!t.IsActive()) continue;
6311 // follow prolongation to the first layer
6312 if ( (fSectors ==fInnerSec) || (t.GetFirstPoint()-fkParam->GetNRowLow()>rfirst+1) )
6313 FollowProlongation(t, rfirst+1);
6318 for (Int_t nr=rfirst; nr>=rlast; nr--){
6319 if (nr<fInnerSec->GetNRows())
6320 fSectors = fInnerSec;
6322 fSectors = fOuterSec;
6323 // make indexes with the cluster tracks for given
6325 // find nearest cluster
6326 for (Int_t i=0; i<nseed; i++) {
6327 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i), &t=*pt;
6329 if (nr==80) pt->UpdateReference();
6330 if (!pt->IsActive()) continue;
6331 // if ( (fSectors ==fOuterSec) && (pt->fFirstPoint-fkParam->GetNRowLow())<nr) continue;
6332 if (pt->GetRelativeSector()>17) {
6335 UpdateClusters(t,nr);
6337 // prolonagate to the nearest cluster - if founded
6338 for (Int_t i=0; i<nseed; i++) {
6339 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
6341 if (!pt->IsActive()) continue;
6342 // if ((fSectors ==fOuterSec) && (pt->fFirstPoint-fkParam->GetNRowLow())<nr) continue;
6343 if (pt->GetRelativeSector()>17) {
6346 FollowToNextCluster(*pt,nr);
6351 void AliTPCtrackerMI::PrepareForBackProlongation(TObjArray *const arr,Float_t fac) const
6355 // if we use TPC track itself we have to "update" covariance
6357 Int_t nseed= arr->GetEntriesFast();
6358 for (Int_t i=0;i<nseed;i++){
6359 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
6363 //rotate to current local system at first accepted point
6364 Int_t index = pt->GetClusterIndex2(pt->GetFirstPoint());
6365 Int_t sec = (index&0xff000000)>>24;
6367 Float_t angle1 = fInnerSec->GetAlpha()*sec+fInnerSec->GetAlphaShift();
6368 if (angle1>TMath::Pi())
6369 angle1-=2.*TMath::Pi();
6370 Float_t angle2 = pt->GetAlpha();
6372 if (TMath::Abs(angle1-angle2)>0.001){
6373 pt->Rotate(angle1-angle2);
6374 //angle2 = pt->GetAlpha();
6375 //pt->fRelativeSector = pt->GetAlpha()/fInnerSec->GetAlpha();
6376 //if (pt->GetAlpha()<0)
6377 // pt->fRelativeSector+=18;
6378 //sec = pt->fRelativeSector;
6387 void AliTPCtrackerMI::PrepareForProlongation(TObjArray *const arr, Float_t fac) const
6391 // if we use TPC track itself we have to "update" covariance
6393 Int_t nseed= arr->GetEntriesFast();
6394 for (Int_t i=0;i<nseed;i++){
6395 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
6398 pt->SetFirstPoint(pt->GetLastPoint());
6406 Int_t AliTPCtrackerMI::PropagateBack(TObjArray *const arr)
6409 // make back propagation
6411 Int_t nseed= arr->GetEntriesFast();
6412 for (Int_t i=0;i<nseed;i++){
6413 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
6414 if (pt&& pt->GetKinkIndex(0)<=0) {
6415 //AliTPCseed *pt2 = new AliTPCseed(*pt);
6416 fSectors = fInnerSec;
6417 //FollowBackProlongation(*pt,fInnerSec->GetNRows()-1);
6418 //fSectors = fOuterSec;
6419 FollowBackProlongation(*pt,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1);
6420 //if (pt->GetNumberOfClusters()<(pt->fEsd->GetTPCclusters(0)) ){
6421 // Error("PropagateBack","Not prolonged track %d",pt->GetLabel());
6422 // FollowBackProlongation(*pt2,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1);
6425 if (pt&& pt->GetKinkIndex(0)>0) {
6426 AliESDkink * kink = fEvent->GetKink(pt->GetKinkIndex(0)-1);
6427 pt->SetFirstPoint(kink->GetTPCRow0());
6428 fSectors = fInnerSec;
6429 FollowBackProlongation(*pt,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1);
6437 Int_t AliTPCtrackerMI::PropagateForward2(TObjArray *const arr)
6440 // make forward propagation
6442 Int_t nseed= arr->GetEntriesFast();
6444 for (Int_t i=0;i<nseed;i++){
6445 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
6447 FollowProlongation(*pt,0);
6456 Int_t AliTPCtrackerMI::PropagateForward()
6459 // propagate track forward
6461 Int_t nseed = fSeeds->GetEntriesFast();
6462 for (Int_t i=0;i<nseed;i++){
6463 AliTPCseed *pt = (AliTPCseed*)fSeeds->UncheckedAt(i);
6465 AliTPCseed &t = *pt;
6466 Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
6467 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
6468 if (alpha < 0. ) alpha += 2.*TMath::Pi();
6469 t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN);
6473 fSectors = fOuterSec;
6474 ParallelTracking(fSeeds,fOuterSec->GetNRows()+fInnerSec->GetNRows()-1,fInnerSec->GetNRows());
6475 fSectors = fInnerSec;
6476 ParallelTracking(fSeeds,fInnerSec->GetNRows()-1,0);
6485 Int_t AliTPCtrackerMI::PropagateBack(AliTPCseed *const pt, Int_t row0, Int_t row1)
6488 // make back propagation, in between row0 and row1
6492 fSectors = fInnerSec;
6495 if (row1<fSectors->GetNRows())
6498 r1 = fSectors->GetNRows()-1;
6500 if (row0<fSectors->GetNRows()&& r1>0 )
6501 FollowBackProlongation(*pt,r1);
6502 if (row1<=fSectors->GetNRows())
6505 r1 = row1 - fSectors->GetNRows();
6506 if (r1<=0) return 0;
6507 if (r1>=fOuterSec->GetNRows()) return 0;
6508 fSectors = fOuterSec;
6509 return FollowBackProlongation(*pt,r1);
6517 void AliTPCtrackerMI::GetShape(AliTPCseed * seed, Int_t row)
6521 AliTPCClusterParam * clparam = AliTPCcalibDB::Instance()->GetClusterParam();
6522 Float_t zdrift = TMath::Abs((fkParam->GetZLength(0)-TMath::Abs(seed->GetZ())));
6523 Int_t type = (seed->GetSector() < fkParam->GetNSector()/2) ? 0: (row>126) ? 1:2;
6524 Double_t angulary = seed->GetSnp();
6526 if (TMath::Abs(angulary)>AliTPCReconstructor::GetMaxSnpTracker()) {
6527 angulary = TMath::Sign(AliTPCReconstructor::GetMaxSnpTracker(),angulary);
6530 angulary = angulary*angulary/((1.-angulary)*(1.+angulary));
6531 Double_t angularz = seed->GetTgl()*seed->GetTgl()*(1.+angulary);
6533 Double_t sigmay = clparam->GetRMS0(0,type,zdrift,TMath::Sqrt(TMath::Abs(angulary)));
6534 Double_t sigmaz = clparam->GetRMS0(1,type,zdrift,TMath::Sqrt(TMath::Abs(angularz)));
6535 seed->SetCurrentSigmaY2(sigmay*sigmay);
6536 seed->SetCurrentSigmaZ2(sigmaz*sigmaz);
6537 // Float_t sd2 = TMath::Abs((fkParam->GetZLength(0)-TMath::Abs(seed->GetZ())))*fkParam->GetDiffL()*fkParam->GetDiffL();
6538 // // Float_t padlength = fkParam->GetPadPitchLength(seed->fSector);
6539 // Float_t padlength = GetPadPitchLength(row);
6541 // Float_t sresy = (seed->GetSector() < fkParam->GetNSector()/2) ? 0.2 :0.3;
6542 // seed->SetCurrentSigmaY2(sd2+padlength*padlength*angulary/12.+sresy*sresy);
6544 // Float_t sresz = fkParam->GetZSigma();
6545 // seed->SetCurrentSigmaZ2(sd2+padlength*padlength*angularz*angularz*(1+angulary)/12.+sresz*sresz);
6547 Float_t wy = GetSigmaY(seed);
6548 Float_t wz = GetSigmaZ(seed);
6551 if (TMath::Abs(wy/seed->fCurrentSigmaY2-1)>0.0001 || TMath::Abs(wz/seed->fCurrentSigmaZ2-1)>0.0001 ){
6552 printf("problem\n");
6559 //__________________________________________________________________________
6560 void AliTPCtrackerMI::CookLabel(AliKalmanTrack *tk, Float_t wrong) const {
6561 //--------------------------------------------------------------------
6562 //This function "cooks" a track label. If label<0, this track is fake.
6563 //--------------------------------------------------------------------
6564 AliTPCseed * t = dynamic_cast<AliTPCseed*>(tk);
6566 printf("%s:%d wrong type \n",(char*)__FILE__,__LINE__);
6570 Int_t noc=t->GetNumberOfClusters();
6572 //printf("\nnot founded prolongation\n\n\n");
6578 AliTPCclusterMI *clusters[160];
6580 for (Int_t i=0;i<160;i++) {
6587 for (i=0; i<160 && current<noc; i++) {
6589 Int_t index=t->GetClusterIndex2(i);
6590 if (index<=0) continue;
6591 if (index&0x8000) continue;
6593 //clusters[current]=GetClusterMI(index);
6594 if (t->GetClusterPointer(i)){
6595 clusters[current]=t->GetClusterPointer(i);
6601 Int_t lab=123456789;
6602 for (i=0; i<noc; i++) {
6603 AliTPCclusterMI *c=clusters[i];
6605 lab=TMath::Abs(c->GetLabel(0));
6607 for (j=0; j<noc; j++) if (lb[j]==lab || mx[j]==0) break;
6613 for (i=0; i<noc; i++) if (mx[i]>max) {max=mx[i]; lab=lb[i];}
6615 for (i=0; i<noc; i++) {
6616 AliTPCclusterMI *c=clusters[i];
6618 if (TMath::Abs(c->GetLabel(1)) == lab ||
6619 TMath::Abs(c->GetLabel(2)) == lab ) max++;
6621 if (noc<=0) { lab=-1; return;}
6622 if ((1.- Float_t(max)/(noc)) > wrong) lab=-lab;
6625 Int_t tail=Int_t(0.10*noc);
6628 for (i=1; i<160&&ind<tail; i++) {
6629 // AliTPCclusterMI *c=clusters[noc-i];
6630 AliTPCclusterMI *c=clusters[i];
6632 if (lab == TMath::Abs(c->GetLabel(0)) ||
6633 lab == TMath::Abs(c->GetLabel(1)) ||
6634 lab == TMath::Abs(c->GetLabel(2))) max++;
6637 if (max < Int_t(0.5*tail)) lab=-lab;
6644 //delete[] clusters;
6648 //__________________________________________________________________________
6649 Int_t AliTPCtrackerMI::CookLabel(AliTPCseed *const t, Float_t wrong,Int_t first, Int_t last) const {
6650 //--------------------------------------------------------------------
6651 //This function "cooks" a track label. If label<0, this track is fake.
6652 //--------------------------------------------------------------------
6653 Int_t noc=t->GetNumberOfClusters();
6655 //printf("\nnot founded prolongation\n\n\n");
6661 AliTPCclusterMI *clusters[160];
6663 for (Int_t i=0;i<160;i++) {
6670 for (i=0; i<160 && current<noc; i++) {
6671 if (i<first) continue;
6672 if (i>last) continue;
6673 Int_t index=t->GetClusterIndex2(i);
6674 if (index<=0) continue;
6675 if (index&0x8000) continue;
6677 //clusters[current]=GetClusterMI(index);
6678 if (t->GetClusterPointer(i)){
6679 clusters[current]=t->GetClusterPointer(i);
6684 if (noc<5) return -1;
6685 Int_t lab=123456789;
6686 for (i=0; i<noc; i++) {
6687 AliTPCclusterMI *c=clusters[i];
6689 lab=TMath::Abs(c->GetLabel(0));
6691 for (j=0; j<noc; j++) if (lb[j]==lab || mx[j]==0) break;
6697 for (i=0; i<noc; i++) if (mx[i]>max) {max=mx[i]; lab=lb[i];}
6699 for (i=0; i<noc; i++) {
6700 AliTPCclusterMI *c=clusters[i];
6702 if (TMath::Abs(c->GetLabel(1)) == lab ||
6703 TMath::Abs(c->GetLabel(2)) == lab ) max++;
6705 if (noc<=0) { lab=-1; return -1;}
6706 if ((1.- Float_t(max)/(noc)) > wrong) lab=-lab;
6709 Int_t tail=Int_t(0.10*noc);
6712 for (i=1; i<160&&ind<tail; i++) {
6713 // AliTPCclusterMI *c=clusters[noc-i];
6714 AliTPCclusterMI *c=clusters[i];
6716 if (lab == TMath::Abs(c->GetLabel(0)) ||
6717 lab == TMath::Abs(c->GetLabel(1)) ||
6718 lab == TMath::Abs(c->GetLabel(2))) max++;
6721 if (max < Int_t(0.5*tail)) lab=-lab;
6724 // t->SetLabel(lab);
6728 //delete[] clusters;
6732 Int_t AliTPCtrackerMI::GetRowNumber(Double_t x[3]) const
6734 //return pad row number for given x vector
6735 Float_t phi = TMath::ATan2(x[1],x[0]);
6736 if(phi<0) phi=2.*TMath::Pi()+phi;
6737 // Get the local angle in the sector philoc
6738 const Float_t kRaddeg = 180/3.14159265358979312;
6739 Float_t phiangle = (Int_t (phi*kRaddeg/20.) + 0.5)*20./kRaddeg;
6740 Double_t localx = x[0]*TMath::Cos(phiangle)-x[1]*TMath::Sin(phiangle);
6741 return GetRowNumber(localx);
6746 void AliTPCtrackerMI::MakeBitmaps(AliTPCseed *t)
6748 //-----------------------------------------------------------------------
6749 // Fill the cluster and sharing bitmaps of the track
6750 //-----------------------------------------------------------------------
6752 Int_t firstpoint = 0;
6753 Int_t lastpoint = 159;
6754 AliTPCTrackerPoint *point;
6755 AliTPCclusterMI *cluster;
6757 for (int iter=firstpoint; iter<lastpoint; iter++) {
6758 // Change to cluster pointers to see if we have a cluster at given padrow
6759 cluster = t->GetClusterPointer(iter);
6761 t->SetClusterMapBit(iter, kTRUE);
6762 point = t->GetTrackPoint(iter);
6763 if (point->IsShared())
6764 t->SetSharedMapBit(iter,kTRUE);
6766 t->SetSharedMapBit(iter, kFALSE);
6769 t->SetClusterMapBit(iter, kFALSE);
6770 t->SetSharedMapBit(iter, kFALSE);
6775 Bool_t AliTPCtrackerMI::IsFindable(AliTPCseed & track){
6777 // return flag if there is findable cluster at given position
6780 Float_t z = track.GetZ();
6782 if (TMath::Abs(z)<(AliTPCReconstructor::GetCtgRange()*track.GetX()+kDeltaZ) &&
6783 TMath::Abs(z)<fkParam->GetZLength(0) &&
6784 (TMath::Abs(track.GetSnp())<AliTPCReconstructor::GetMaxSnpTracker()))
6790 void AliTPCtrackerMI::AddCovariance(AliTPCseed * seed){
6792 // Adding systematic error
6793 // !!!! the systematic error for element 4 is in 1/cm not in pt
6795 const Double_t *param = AliTPCReconstructor::GetRecoParam()->GetSystematicError();
6796 Double_t *covarIn= (Double_t*)seed->GetCovariance();
6798 for (Int_t i=0;i<15;i++) covar[i]=0;
6804 covar[0] = param[0]*param[0];
6805 covar[2] = param[1]*param[1];
6806 covar[5] = param[2]*param[2];
6807 covar[9] = param[3]*param[3];
6808 Double_t facC = AliTracker::GetBz()*kB2C;
6809 covar[14]= param[4]*param[4]*facC*facC;
6811 covar[1]=TMath::Sqrt((covar[0]*covar[2]))*covarIn[1]/TMath::Sqrt((covarIn[0]*covarIn[2]));
6813 covar[3]=TMath::Sqrt((covar[0]*covar[5]))*covarIn[3]/TMath::Sqrt((covarIn[0]*covarIn[5]));
6814 covar[4]=TMath::Sqrt((covar[2]*covar[5]))*covarIn[4]/TMath::Sqrt((covarIn[2]*covarIn[5]));
6816 covar[6]=TMath::Sqrt((covar[0]*covar[9]))*covarIn[6]/TMath::Sqrt((covarIn[0]*covarIn[9]));
6817 covar[7]=TMath::Sqrt((covar[2]*covar[9]))*covarIn[7]/TMath::Sqrt((covarIn[2]*covarIn[9]));
6818 covar[8]=TMath::Sqrt((covar[5]*covar[9]))*covarIn[8]/TMath::Sqrt((covarIn[5]*covarIn[9]));
6820 covar[10]=TMath::Sqrt((covar[0]*covar[14]))*covarIn[10]/TMath::Sqrt((covarIn[0]*covarIn[14]));
6821 covar[11]=TMath::Sqrt((covar[2]*covar[14]))*covarIn[11]/TMath::Sqrt((covarIn[2]*covarIn[14]));
6822 covar[12]=TMath::Sqrt((covar[5]*covar[14]))*covarIn[12]/TMath::Sqrt((covarIn[5]*covarIn[14]));
6823 covar[13]=TMath::Sqrt((covar[9]*covar[14]))*covarIn[13]/TMath::Sqrt((covarIn[9]*covarIn[14]));
6825 seed->AddCovariance(covar);