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>
116 #include <TGraphErrors.h>
118 #include "AliComplexCluster.h"
119 #include "AliESDEvent.h"
120 #include "AliESDtrack.h"
121 #include "AliESDVertex.h"
124 #include "AliHelix.h"
125 #include "AliRunLoader.h"
126 #include "AliTPCClustersRow.h"
127 #include "AliTPCParam.h"
128 #include "AliTPCReconstructor.h"
129 #include "AliTPCpolyTrack.h"
130 #include "AliTPCreco.h"
131 #include "AliTPCseed.h"
133 #include "AliTPCtrackerSector.h"
134 #include "AliTPCtrackerMI.h"
135 #include "TStopwatch.h"
136 #include "AliTPCReconstructor.h"
137 #include "AliAlignObj.h"
138 #include "AliTrackPointArray.h"
140 #include "AliTPCcalibDB.h"
141 #include "AliTPCcalibDButil.h"
142 #include "AliTPCTransform.h"
143 #include "AliTPCClusterParam.h"
147 ClassImp(AliTPCtrackerMI)
151 class AliTPCFastMath {
154 static Double_t FastAsin(Double_t x);
156 static Double_t fgFastAsin[20000]; //lookup table for fast asin computation
159 Double_t AliTPCFastMath::fgFastAsin[20000];
160 AliTPCFastMath gAliTPCFastMath; // needed to fill the LUT
162 AliTPCFastMath::AliTPCFastMath(){
164 // initialized lookup table;
165 for (Int_t i=0;i<10000;i++){
166 fgFastAsin[2*i] = TMath::ASin(i/10000.);
167 fgFastAsin[2*i+1] = (TMath::ASin((i+1)/10000.)-fgFastAsin[2*i]);
171 Double_t AliTPCFastMath::FastAsin(Double_t x){
173 // return asin using lookup table
175 Int_t index = int(x*10000);
176 return fgFastAsin[2*index]+(x*10000.-index)*fgFastAsin[2*index+1];
179 Int_t index = int(x*10000);
180 return -(fgFastAsin[2*index]+(x*10000.-index)*fgFastAsin[2*index+1]);
182 //__________________________________________________________________
183 AliTPCtrackerMI::AliTPCtrackerMI()
205 // default constructor
207 for (Int_t irow=0; irow<200; irow++){
214 //_____________________________________________________________________
218 Int_t AliTPCtrackerMI::UpdateTrack(AliTPCseed * track, Int_t accept){
220 //update track information using current cluster - track->fCurrentCluster
223 AliTPCclusterMI* c =track->GetCurrentCluster();
224 if (accept>0) track->SetCurrentClusterIndex1(track->GetCurrentClusterIndex1() | 0x8000); //sign not accepted clusters
226 UInt_t i = track->GetCurrentClusterIndex1();
228 Int_t sec=(i&0xff000000)>>24;
229 //Int_t row = (i&0x00ff0000)>>16;
230 track->SetRow((i&0x00ff0000)>>16);
231 track->SetSector(sec);
232 // Int_t index = i&0xFFFF;
233 if (sec>=fkParam->GetNInnerSector()) track->SetRow(track->GetRow()+fkParam->GetNRowLow());
234 track->SetClusterIndex2(track->GetRow(), i);
235 //track->fFirstPoint = row;
236 //if ( track->fLastPoint<row) track->fLastPoint =row;
237 // if (track->fRow<0 || track->fRow>160) {
238 // printf("problem\n");
240 if (track->GetFirstPoint()>track->GetRow())
241 track->SetFirstPoint(track->GetRow());
242 if (track->GetLastPoint()<track->GetRow())
243 track->SetLastPoint(track->GetRow());
246 track->SetClusterPointer(track->GetRow(),c);
249 Double_t angle2 = track->GetSnp()*track->GetSnp();
251 //SET NEW Track Point
253 if (angle2<1) //PH sometimes angle2 is very big. To be investigated...
255 angle2 = TMath::Sqrt(angle2/(1-angle2));
256 AliTPCTrackerPoint &point =*(track->GetTrackPoint(track->GetRow()));
258 point.SetSigmaY(c->GetSigmaY2()/track->GetCurrentSigmaY2());
259 point.SetSigmaZ(c->GetSigmaZ2()/track->GetCurrentSigmaZ2());
260 point.SetErrY(sqrt(track->GetErrorY2()));
261 point.SetErrZ(sqrt(track->GetErrorZ2()));
263 point.SetX(track->GetX());
264 point.SetY(track->GetY());
265 point.SetZ(track->GetZ());
266 point.SetAngleY(angle2);
267 point.SetAngleZ(track->GetTgl());
268 if (point.IsShared()){
269 track->SetErrorY2(track->GetErrorY2()*4);
270 track->SetErrorZ2(track->GetErrorZ2()*4);
274 Double_t chi2 = track->GetPredictedChi2(track->GetCurrentCluster());
276 // track->SetErrorY2(track->GetErrorY2()*1.3);
277 // track->SetErrorY2(track->GetErrorY2()+0.01);
278 // track->SetErrorZ2(track->GetErrorZ2()*1.3);
279 // track->SetErrorZ2(track->GetErrorZ2()+0.005);
281 if (accept>0) return 0;
282 if (track->GetNumberOfClusters()%20==0){
283 // if (track->fHelixIn){
284 // TClonesArray & larr = *(track->fHelixIn);
285 // Int_t ihelix = larr.GetEntriesFast();
286 // new(larr[ihelix]) AliHelix(*track) ;
289 track->SetNoCluster(0);
290 return track->Update(c,chi2,i);
295 Int_t AliTPCtrackerMI::AcceptCluster(AliTPCseed * seed, AliTPCclusterMI * cluster)
298 // decide according desired precision to accept given
299 // cluster for tracking
301 seed->GetProlongation(cluster->GetX(),yt,zt);
302 Double_t sy2=ErrY2(seed,cluster);
303 Double_t sz2=ErrZ2(seed,cluster);
305 Double_t sdistancey2 = sy2+seed->GetSigmaY2();
306 Double_t sdistancez2 = sz2+seed->GetSigmaZ2();
307 Double_t dy=seed->GetCurrentCluster()->GetY()-yt;
308 Double_t dz=seed->GetCurrentCluster()->GetZ()-zt;
309 Double_t rdistancey2 = (seed->GetCurrentCluster()->GetY()-yt)*
310 (seed->GetCurrentCluster()->GetY()-yt)/sdistancey2;
311 Double_t rdistancez2 = (seed->GetCurrentCluster()->GetZ()-zt)*
312 (seed->GetCurrentCluster()->GetZ()-zt)/sdistancez2;
314 Double_t rdistance2 = rdistancey2+rdistancez2;
317 if (AliTPCReconstructor::StreamLevel()>2 && seed->GetNumberOfClusters()>20) {
318 Float_t rmsy2 = seed->GetCurrentSigmaY2();
319 Float_t rmsz2 = seed->GetCurrentSigmaZ2();
320 Float_t rmsy2p30 = seed->GetCMeanSigmaY2p30();
321 Float_t rmsz2p30 = seed->GetCMeanSigmaZ2p30();
322 Float_t rmsy2p30R = seed->GetCMeanSigmaY2p30R();
323 Float_t rmsz2p30R = seed->GetCMeanSigmaZ2p30R();
324 AliExternalTrackParam param(*seed);
325 static TVectorD gcl(3),gtr(3);
327 param.GetXYZ(gcl.GetMatrixArray());
328 cluster->GetGlobalXYZ(gclf);
329 gcl[0]=gclf[0]; gcl[1]=gclf[1]; gcl[2]=gclf[2];
332 if (AliTPCReconstructor::StreamLevel()>2) {
333 (*fDebugStreamer)<<"ErrParam"<<
346 "rmsy2p30="<<rmsy2p30<<
347 "rmsz2p30="<<rmsz2p30<<
348 "rmsy2p30R="<<rmsy2p30R<<
349 "rmsz2p30R="<<rmsz2p30R<<
350 // normalize distance -
351 "rdisty="<<rdistancey2<<
352 "rdistz="<<rdistancez2<<
353 "rdist="<<rdistance2<< //
357 //return 0; // temporary
358 if (rdistance2>32) return 3;
361 if ((rdistancey2>9. || rdistancez2>9.) && cluster->GetType()==0)
362 return 2; //suspisiouce - will be changed
364 if ((rdistancey2>6.25 || rdistancez2>6.25) && cluster->GetType()>0)
365 // strict cut on overlaped cluster
366 return 2; //suspisiouce - will be changed
368 if ( (rdistancey2>1. || rdistancez2>6.25 )
369 && cluster->GetType()<0){
370 seed->SetNFoundable(seed->GetNFoundable()-1);
380 //_____________________________________________________________________________
381 AliTPCtrackerMI::AliTPCtrackerMI(const AliTPCParam *par):
383 fkNIS(par->GetNInnerSector()/2),
385 fkNOS(par->GetNOuterSector()/2),
402 //---------------------------------------------------------------------
403 // The main TPC tracker constructor
404 //---------------------------------------------------------------------
405 fInnerSec=new AliTPCtrackerSector[fkNIS];
406 fOuterSec=new AliTPCtrackerSector[fkNOS];
409 for (i=0; i<fkNIS; i++) fInnerSec[i].Setup(par,0);
410 for (i=0; i<fkNOS; i++) fOuterSec[i].Setup(par,1);
413 Int_t nrowlow = par->GetNRowLow();
414 Int_t nrowup = par->GetNRowUp();
417 for (i=0;i<nrowlow;i++){
418 fXRow[i] = par->GetPadRowRadiiLow(i);
419 fPadLength[i]= par->GetPadPitchLength(0,i);
420 fYMax[i] = fXRow[i]*TMath::Tan(0.5*par->GetInnerAngle());
424 for (i=0;i<nrowup;i++){
425 fXRow[i+nrowlow] = par->GetPadRowRadiiUp(i);
426 fPadLength[i+nrowlow] = par->GetPadPitchLength(60,i);
427 fYMax[i+nrowlow] = fXRow[i+nrowlow]*TMath::Tan(0.5*par->GetOuterAngle());
430 if (AliTPCReconstructor::StreamLevel()>0) {
431 fDebugStreamer = new TTreeSRedirector("TPCdebug.root");
434 //________________________________________________________________________
435 AliTPCtrackerMI::AliTPCtrackerMI(const AliTPCtrackerMI &t):
456 //------------------------------------
457 // dummy copy constructor
458 //------------------------------------------------------------------
460 for (Int_t irow=0; irow<200; irow++){
467 AliTPCtrackerMI & AliTPCtrackerMI::operator=(const AliTPCtrackerMI& /*r*/)
469 //------------------------------
471 //--------------------------------------------------------------
474 //_____________________________________________________________________________
475 AliTPCtrackerMI::~AliTPCtrackerMI() {
476 //------------------------------------------------------------------
477 // TPC tracker destructor
478 //------------------------------------------------------------------
485 if (fDebugStreamer) delete fDebugStreamer;
489 void AliTPCtrackerMI::FillESD(const TObjArray* arr)
493 //fill esds using updated tracks
496 // write tracks to the event
497 // store index of the track
498 Int_t nseed=arr->GetEntriesFast();
499 //FindKinks(arr,fEvent);
500 for (Int_t i=0; i<nseed; i++) {
501 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
505 if (AliTPCReconstructor::StreamLevel()>1) {
506 (*fDebugStreamer)<<"Track0"<<
510 // pt->PropagateTo(fkParam->GetInnerRadiusLow());
511 if (pt->GetKinkIndex(0)<=0){ //don't propagate daughter tracks
512 pt->PropagateTo(fkParam->GetInnerRadiusLow());
515 if (( pt->GetPoints()[2]- pt->GetPoints()[0])>5 && pt->GetPoints()[3]>0.8){
517 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
518 iotrack.SetTPCPoints(pt->GetPoints());
519 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
520 iotrack.SetV0Indexes(pt->GetV0Indexes());
521 // iotrack.SetTPCpid(pt->fTPCr);
522 //iotrack.SetTPCindex(i);
523 fEvent->AddTrack(&iotrack);
527 if ( (pt->GetNumberOfClusters()>70)&& (Float_t(pt->GetNumberOfClusters())/Float_t(pt->GetNFoundable()))>0.55) {
529 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
530 iotrack.SetTPCPoints(pt->GetPoints());
531 //iotrack.SetTPCindex(i);
532 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
533 iotrack.SetV0Indexes(pt->GetV0Indexes());
534 // iotrack.SetTPCpid(pt->fTPCr);
535 fEvent->AddTrack(&iotrack);
539 // short tracks - maybe decays
541 if ( (pt->GetNumberOfClusters()>30) && (Float_t(pt->GetNumberOfClusters())/Float_t(pt->GetNFoundable()))>0.70) {
542 Int_t found,foundable,shared;
543 pt->GetClusterStatistic(0,60,found, foundable,shared,kFALSE);
544 if ( (found>20) && (pt->GetNShared()/float(pt->GetNumberOfClusters())<0.2)){
546 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
547 //iotrack.SetTPCindex(i);
548 iotrack.SetTPCPoints(pt->GetPoints());
549 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
550 iotrack.SetV0Indexes(pt->GetV0Indexes());
551 //iotrack.SetTPCpid(pt->fTPCr);
552 fEvent->AddTrack(&iotrack);
557 if ( (pt->GetNumberOfClusters()>20) && (Float_t(pt->GetNumberOfClusters())/Float_t(pt->GetNFoundable()))>0.8) {
558 Int_t found,foundable,shared;
559 pt->GetClusterStatistic(0,60,found, foundable,shared,kFALSE);
560 if (found<20) continue;
561 if (pt->GetNShared()/float(pt->GetNumberOfClusters())>0.2) continue;
564 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
565 iotrack.SetTPCPoints(pt->GetPoints());
566 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
567 iotrack.SetV0Indexes(pt->GetV0Indexes());
568 //iotrack.SetTPCpid(pt->fTPCr);
569 //iotrack.SetTPCindex(i);
570 fEvent->AddTrack(&iotrack);
573 // short tracks - secondaties
575 if ( (pt->GetNumberOfClusters()>30) ) {
576 Int_t found,foundable,shared;
577 pt->GetClusterStatistic(128,158,found, foundable,shared,kFALSE);
578 if ( (found>20) && (pt->GetNShared()/float(pt->GetNumberOfClusters())<0.2) &&float(found)/float(foundable)>0.8){
580 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
581 iotrack.SetTPCPoints(pt->GetPoints());
582 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
583 iotrack.SetV0Indexes(pt->GetV0Indexes());
584 //iotrack.SetTPCpid(pt->fTPCr);
585 //iotrack.SetTPCindex(i);
586 fEvent->AddTrack(&iotrack);
591 if ( (pt->GetNumberOfClusters()>15)) {
592 Int_t found,foundable,shared;
593 pt->GetClusterStatistic(138,158,found, foundable,shared,kFALSE);
594 if (found<15) continue;
595 if (foundable<=0) continue;
596 if (pt->GetNShared()/float(pt->GetNumberOfClusters())>0.2) continue;
597 if (float(found)/float(foundable)<0.8) continue;
600 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
601 iotrack.SetTPCPoints(pt->GetPoints());
602 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
603 iotrack.SetV0Indexes(pt->GetV0Indexes());
604 // iotrack.SetTPCpid(pt->fTPCr);
605 //iotrack.SetTPCindex(i);
606 fEvent->AddTrack(&iotrack);
610 // >> account for suppressed tracks in the kink indices (RS)
611 int nESDtracks = fEvent->GetNumberOfTracks();
612 for (int it=nESDtracks;it--;) {
613 AliESDtrack* esdTr = fEvent->GetTrack(it);
614 if (!esdTr || !esdTr->GetKinkIndex(0)) continue;
615 for (int ik=0;ik<3;ik++) {
617 if (!(knkId=esdTr->GetKinkIndex(ik))) break; // no more kinks for this track
618 AliESDkink* kink = fEvent->GetKink(TMath::Abs(knkId)-1);
620 AliError(Form("ESDTrack%d refers to non-existing kink %d",it,TMath::Abs(knkId)-1));
623 kink->SetIndex(it, knkId<0 ? 0:1); // update track index of the kink: mother at 0, daughter at 1
626 // << account for suppressed tracks in the kink indices (RS)
627 AliInfo(Form("Number of filled ESDs-\t%d\n",fEvent->GetNumberOfTracks()));
635 Double_t AliTPCtrackerMI::ErrY2(AliTPCseed* seed, const AliTPCclusterMI * cl){
638 // Use calibrated cluster error from OCDB
640 AliTPCClusterParam * clparam = AliTPCcalibDB::Instance()->GetClusterParam();
642 Float_t z = TMath::Abs(fkParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
643 Int_t ctype = cl->GetType();
644 Int_t type = (cl->GetRow()<63) ? 0: (cl->GetRow()>126) ? 1:2;
645 Double_t angle = seed->GetSnp()*seed->GetSnp();
646 angle = TMath::Sqrt(TMath::Abs(angle/(1.-angle)));
647 Double_t erry2 = clparam->GetError0Par(0,type, z,angle);
649 erry2+=0.5; // edge cluster
652 seed->SetErrorY2(erry2);
656 //calculate look-up table at the beginning
657 // static Bool_t ginit = kFALSE;
658 // static Float_t gnoise1,gnoise2,gnoise3;
659 // static Float_t ggg1[10000];
660 // static Float_t ggg2[10000];
661 // static Float_t ggg3[10000];
662 // static Float_t glandau1[10000];
663 // static Float_t glandau2[10000];
664 // static Float_t glandau3[10000];
666 // static Float_t gcor01[500];
667 // static Float_t gcor02[500];
668 // static Float_t gcorp[500];
672 // if (ginit==kFALSE){
673 // for (Int_t i=1;i<500;i++){
674 // Float_t rsigma = float(i)/100.;
675 // gcor02[i] = TMath::Max(0.78 +TMath::Exp(7.4*(rsigma-1.2)),0.6);
676 // gcor01[i] = TMath::Max(0.72 +TMath::Exp(3.36*(rsigma-1.2)),0.6);
677 // gcorp[i] = TMath::Max(TMath::Power((rsigma+0.5),1.5),1.2);
681 // for (Int_t i=3;i<10000;i++){
685 // Float_t amp = float(i);
686 // Float_t padlength =0.75;
687 // gnoise1 = 0.0004/padlength;
688 // Float_t nel = 0.268*amp;
689 // Float_t nprim = 0.155*amp;
690 // ggg1[i] = fkParam->GetDiffT()*fkParam->GetDiffT()*(2+0.001*nel/(padlength*padlength))/nel;
691 // glandau1[i] = (2.+0.12*nprim)*0.5* (2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
692 // if (glandau1[i]>1) glandau1[i]=1;
693 // glandau1[i]*=padlength*padlength/12.;
697 // gnoise2 = 0.0004/padlength;
699 // nprim = 0.133*amp;
700 // ggg2[i] = fkParam->GetDiffT()*fkParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
701 // glandau2[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
702 // if (glandau2[i]>1) glandau2[i]=1;
703 // glandau2[i]*=padlength*padlength/12.;
708 // gnoise3 = 0.0004/padlength;
710 // nprim = 0.133*amp;
711 // ggg3[i] = fkParam->GetDiffT()*fkParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
712 // glandau3[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
713 // if (glandau3[i]>1) glandau3[i]=1;
714 // glandau3[i]*=padlength*padlength/12.;
722 // Int_t amp = int(TMath::Abs(cl->GetQ()));
724 // seed->SetErrorY2(1.);
728 // Float_t z = TMath::Abs(fkParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
729 // Int_t ctype = cl->GetType();
730 // Float_t padlength= GetPadPitchLength(seed->GetRow());
731 // Double_t angle2 = seed->GetSnp()*seed->GetSnp();
732 // angle2 = angle2/(1-angle2);
734 // //cluster "quality"
735 // Int_t rsigmay = int(100.*cl->GetSigmaY2()/(seed->GetCurrentSigmaY2()));
738 // if (fSectors==fInnerSec){
739 // snoise2 = gnoise1;
740 // res = ggg1[amp]*z+glandau1[amp]*angle2;
741 // if (ctype==0) res *= gcor01[rsigmay];
744 // res*= gcorp[rsigmay];
748 // if (padlength<1.1){
749 // snoise2 = gnoise2;
750 // res = ggg2[amp]*z+glandau2[amp]*angle2;
751 // if (ctype==0) res *= gcor02[rsigmay];
754 // res*= gcorp[rsigmay];
758 // snoise2 = gnoise3;
759 // res = ggg3[amp]*z+glandau3[amp]*angle2;
760 // if (ctype==0) res *= gcor02[rsigmay];
763 // res*= gcorp[rsigmay];
770 // res*=2.4; // overestimate error 2 times
774 // if (res<2*snoise2)
777 // seed->SetErrorY2(res);
785 Double_t AliTPCtrackerMI::ErrZ2(AliTPCseed* seed, const AliTPCclusterMI * cl){
788 // Use calibrated cluster error from OCDB
790 AliTPCClusterParam * clparam = AliTPCcalibDB::Instance()->GetClusterParam();
792 Float_t z = TMath::Abs(fkParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
793 Int_t ctype = cl->GetType();
794 Int_t type = (cl->GetRow()<63) ? 0: (cl->GetRow()>126) ? 1:2;
796 Double_t angle2 = seed->GetSnp()*seed->GetSnp();
797 angle2 = seed->GetTgl()*seed->GetTgl()*(1+angle2/(1-angle2));
798 Double_t angle = TMath::Sqrt(TMath::Abs(angle2));
799 Double_t errz2 = clparam->GetError0Par(1,type, z,angle);
801 errz2+=0.5; // edge cluster
804 seed->SetErrorZ2(errz2);
810 // //seed->SetErrorY2(0.1);
812 // //calculate look-up table at the beginning
813 // static Bool_t ginit = kFALSE;
814 // static Float_t gnoise1,gnoise2,gnoise3;
815 // static Float_t ggg1[10000];
816 // static Float_t ggg2[10000];
817 // static Float_t ggg3[10000];
818 // static Float_t glandau1[10000];
819 // static Float_t glandau2[10000];
820 // static Float_t glandau3[10000];
822 // static Float_t gcor01[1000];
823 // static Float_t gcor02[1000];
824 // static Float_t gcorp[1000];
828 // if (ginit==kFALSE){
829 // for (Int_t i=1;i<1000;i++){
830 // Float_t rsigma = float(i)/100.;
831 // gcor02[i] = TMath::Max(0.81 +TMath::Exp(6.8*(rsigma-1.2)),0.6);
832 // gcor01[i] = TMath::Max(0.72 +TMath::Exp(2.04*(rsigma-1.2)),0.6);
833 // gcorp[i] = TMath::Max(TMath::Power((rsigma+0.5),1.5),1.2);
837 // for (Int_t i=3;i<10000;i++){
841 // Float_t amp = float(i);
842 // Float_t padlength =0.75;
843 // gnoise1 = 0.0004/padlength;
844 // Float_t nel = 0.268*amp;
845 // Float_t nprim = 0.155*amp;
846 // ggg1[i] = fkParam->GetDiffT()*fkParam->GetDiffT()*(2+0.001*nel/(padlength*padlength))/nel;
847 // glandau1[i] = (2.+0.12*nprim)*0.5* (2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
848 // if (glandau1[i]>1) glandau1[i]=1;
849 // glandau1[i]*=padlength*padlength/12.;
853 // gnoise2 = 0.0004/padlength;
855 // nprim = 0.133*amp;
856 // ggg2[i] = fkParam->GetDiffT()*fkParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
857 // glandau2[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
858 // if (glandau2[i]>1) glandau2[i]=1;
859 // glandau2[i]*=padlength*padlength/12.;
864 // gnoise3 = 0.0004/padlength;
866 // nprim = 0.133*amp;
867 // ggg3[i] = fkParam->GetDiffT()*fkParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
868 // glandau3[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
869 // if (glandau3[i]>1) glandau3[i]=1;
870 // glandau3[i]*=padlength*padlength/12.;
878 // Int_t amp = int(TMath::Abs(cl->GetQ()));
880 // seed->SetErrorY2(1.);
884 // Float_t z = TMath::Abs(fkParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
885 // Int_t ctype = cl->GetType();
886 // Float_t padlength= GetPadPitchLength(seed->GetRow());
888 // Double_t angle2 = seed->GetSnp()*seed->GetSnp();
889 // // if (angle2<0.6) angle2 = 0.6;
890 // angle2 = seed->GetTgl()*seed->GetTgl()*(1+angle2/(1-angle2));
892 // //cluster "quality"
893 // Int_t rsigmaz = int(100.*cl->GetSigmaZ2()/(seed->GetCurrentSigmaZ2()));
896 // if (fSectors==fInnerSec){
897 // snoise2 = gnoise1;
898 // res = ggg1[amp]*z+glandau1[amp]*angle2;
899 // if (ctype==0) res *= gcor01[rsigmaz];
902 // res*= gcorp[rsigmaz];
906 // if (padlength<1.1){
907 // snoise2 = gnoise2;
908 // res = ggg2[amp]*z+glandau2[amp]*angle2;
909 // if (ctype==0) res *= gcor02[rsigmaz];
912 // res*= gcorp[rsigmaz];
916 // snoise2 = gnoise3;
917 // res = ggg3[amp]*z+glandau3[amp]*angle2;
918 // if (ctype==0) res *= gcor02[rsigmaz];
921 // res*= gcorp[rsigmaz];
930 // if ((ctype<0) &&<70){
935 // if (res<2*snoise2)
937 // if (res>3) res =3;
938 // seed->SetErrorZ2(res);
946 void AliTPCtrackerMI::RotateToLocal(AliTPCseed *seed)
948 //rotate to track "local coordinata
949 Float_t x = seed->GetX();
950 Float_t y = seed->GetY();
951 Float_t ymax = x*TMath::Tan(0.5*fSectors->GetAlpha());
954 seed->SetRelativeSector((seed->GetRelativeSector()+1) % fN);
955 if (!seed->Rotate(fSectors->GetAlpha()))
957 } else if (y <-ymax) {
958 seed->SetRelativeSector((seed->GetRelativeSector()-1+fN) % fN);
959 if (!seed->Rotate(-fSectors->GetAlpha()))
967 //_____________________________________________________________________________
968 Double_t AliTPCtrackerMI::F1old(Double_t x1,Double_t y1,
969 Double_t x2,Double_t y2,
970 Double_t x3,Double_t y3) const
972 //-----------------------------------------------------------------
973 // Initial approximation of the track curvature
974 //-----------------------------------------------------------------
975 Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
976 Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
977 (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
978 Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
979 (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
981 Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
982 if ( xr*xr+yr*yr<=0.00000000000001) return 100;
983 return -xr*yr/sqrt(xr*xr+yr*yr);
988 //_____________________________________________________________________________
989 Double_t AliTPCtrackerMI::F1(Double_t x1,Double_t y1,
990 Double_t x2,Double_t y2,
991 Double_t x3,Double_t y3) const
993 //-----------------------------------------------------------------
994 // Initial approximation of the track curvature
995 //-----------------------------------------------------------------
1001 Double_t det = x3*y2-x2*y3;
1002 if (TMath::Abs(det)<1e-10){
1006 Double_t u = 0.5* (x2*(x2-x3)+y2*(y2-y3))/det;
1007 Double_t x0 = x3*0.5-y3*u;
1008 Double_t y0 = y3*0.5+x3*u;
1009 Double_t c2 = 1/TMath::Sqrt(x0*x0+y0*y0);
1015 Double_t AliTPCtrackerMI::F2(Double_t x1,Double_t y1,
1016 Double_t x2,Double_t y2,
1017 Double_t x3,Double_t y3) const
1019 //-----------------------------------------------------------------
1020 // Initial approximation of the track curvature
1021 //-----------------------------------------------------------------
1027 Double_t det = x3*y2-x2*y3;
1028 if (TMath::Abs(det)<1e-10) {
1032 Double_t u = 0.5* (x2*(x2-x3)+y2*(y2-y3))/det;
1033 Double_t x0 = x3*0.5-y3*u;
1034 Double_t y0 = y3*0.5+x3*u;
1035 Double_t c2 = 1/TMath::Sqrt(x0*x0+y0*y0);
1044 //_____________________________________________________________________________
1045 Double_t AliTPCtrackerMI::F2old(Double_t x1,Double_t y1,
1046 Double_t x2,Double_t y2,
1047 Double_t x3,Double_t y3) const
1049 //-----------------------------------------------------------------
1050 // Initial approximation of the track curvature times center of curvature
1051 //-----------------------------------------------------------------
1052 Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
1053 Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
1054 (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
1055 Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
1056 (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
1058 Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
1060 return -a/(d*y1-b)*xr/sqrt(xr*xr+yr*yr);
1063 //_____________________________________________________________________________
1064 Double_t AliTPCtrackerMI::F3(Double_t x1,Double_t y1,
1065 Double_t x2,Double_t y2,
1066 Double_t z1,Double_t z2) const
1068 //-----------------------------------------------------------------
1069 // Initial approximation of the tangent of the track dip angle
1070 //-----------------------------------------------------------------
1071 return (z1 - z2)/sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
1075 Double_t AliTPCtrackerMI::F3n(Double_t x1,Double_t y1,
1076 Double_t x2,Double_t y2,
1077 Double_t z1,Double_t z2, Double_t c) const
1079 //-----------------------------------------------------------------
1080 // Initial approximation of the tangent of the track dip angle
1081 //-----------------------------------------------------------------
1085 //angle1 = (z1-z2)*c/(TMath::ASin(c*x1-ni)-TMath::ASin(c*x2-ni));
1087 Double_t d = TMath::Sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
1088 if (TMath::Abs(d*c*0.5)>1) return 0;
1089 // Double_t angle2 = TMath::ASin(d*c*0.5);
1090 // Double_t angle2 = AliTPCFastMath::FastAsin(d*c*0.5);
1091 Double_t angle2 = (d*c*0.5>0.1)? TMath::ASin(d*c*0.5): AliTPCFastMath::FastAsin(d*c*0.5);
1093 angle2 = (z1-z2)*c/(angle2*2.);
1097 Bool_t AliTPCtrackerMI::GetProlongation(Double_t x1, Double_t x2, Double_t x[5], Double_t &y, Double_t &z) const
1098 {//-----------------------------------------------------------------
1099 // This function find proloncation of a track to a reference plane x=x2.
1100 //-----------------------------------------------------------------
1104 if (TMath::Abs(x[4]*x1 - x[2]) >= 0.999) {
1108 Double_t c1=x[4]*x1 - x[2], r1=TMath::Sqrt((1.-c1)*(1.+c1));
1109 Double_t c2=x[4]*x2 - x[2], r2=TMath::Sqrt((1.-c2)*(1.+c2));
1113 Double_t dy = dx*(c1+c2)/(r1+r2);
1116 Double_t delta = x[4]*dx*(c1+c2)/(c1*r2 + c2*r1);
1118 if (TMath::Abs(delta)>0.01){
1119 dz = x[3]*TMath::ASin(delta)/x[4];
1121 dz = x[3]*AliTPCFastMath::FastAsin(delta)/x[4];
1124 //dz = x[3]*AliTPCFastMath::FastAsin(delta)/x[4];
1132 Int_t AliTPCtrackerMI::LoadClusters (TTree *const tree)
1137 return LoadClusters();
1141 Int_t AliTPCtrackerMI::LoadClusters(const TObjArray *arr)
1144 // load clusters to the memory
1145 AliTPCClustersRow *clrow = new AliTPCClustersRow("AliTPCclusterMI");
1146 Int_t lower = arr->LowerBound();
1147 Int_t entries = arr->GetEntriesFast();
1149 for (Int_t i=lower; i<entries; i++) {
1150 clrow = (AliTPCClustersRow*) arr->At(i);
1151 if(!clrow) continue;
1152 if(!clrow->GetArray()) continue;
1156 fkParam->AdjustSectorRow(clrow->GetID(),sec,row);
1158 for (Int_t icl=0; icl<clrow->GetArray()->GetEntriesFast(); icl++){
1159 Transform((AliTPCclusterMI*)(clrow->GetArray()->At(icl)));
1162 if (clrow->GetArray()->GetEntriesFast()<=0) continue;
1163 AliTPCtrackerRow * tpcrow=0;
1166 tpcrow = &(fInnerSec[sec%fkNIS][row]);
1170 tpcrow = &(fOuterSec[(sec-fkNIS*2)%fkNOS][row]);
1171 left = (sec-fkNIS*2)/fkNOS;
1174 tpcrow->SetN1(clrow->GetArray()->GetEntriesFast());
1175 tpcrow->SetClusters1(new AliTPCclusterMI[tpcrow->GetN1()]);
1176 for (Int_t j=0;j<tpcrow->GetN1();++j)
1177 tpcrow->SetCluster1(j, *(AliTPCclusterMI*)(clrow->GetArray()->At(j)));
1180 tpcrow->SetN2(clrow->GetArray()->GetEntriesFast());
1181 tpcrow->SetClusters2(new AliTPCclusterMI[tpcrow->GetN2()]);
1182 for (Int_t j=0;j<tpcrow->GetN2();++j)
1183 tpcrow->SetCluster2(j,*(AliTPCclusterMI*)(clrow->GetArray()->At(j)));
1185 clrow->GetArray()->Clear("C");
1194 Int_t AliTPCtrackerMI::LoadClusters(const TClonesArray *arr)
1197 // load clusters to the memory from one
1200 AliTPCclusterMI *clust=0;
1201 Int_t count[72][96] = { {0} , {0} };
1203 // loop over clusters
1204 for (Int_t icl=0; icl<arr->GetEntriesFast(); icl++) {
1205 clust = (AliTPCclusterMI*)arr->At(icl);
1206 if(!clust) continue;
1207 //printf("cluster: det %d, row %d \n", clust->GetDetector(),clust->GetRow());
1209 // transform clusters
1212 // count clusters per pad row
1213 count[clust->GetDetector()][clust->GetRow()]++;
1216 // insert clusters to sectors
1217 for (Int_t icl=0; icl<arr->GetEntriesFast(); icl++) {
1218 clust = (AliTPCclusterMI*)arr->At(icl);
1219 if(!clust) continue;
1221 Int_t sec = clust->GetDetector();
1222 Int_t row = clust->GetRow();
1224 // filter overlapping pad rows needed by HLT
1225 if(sec<fkNIS*2) { //IROCs
1226 if(row == 30) continue;
1229 if(row == 27 || row == 76) continue;
1235 fInnerSec[sec%fkNIS].InsertCluster(clust, count[sec][row], fkParam);
1238 left = (sec-fkNIS*2)/fkNOS;
1239 fOuterSec[(sec-fkNIS*2)%fkNOS].InsertCluster(clust, count[sec][row], fkParam);
1243 // Load functions must be called behind LoadCluster(TClonesArray*)
1245 //LoadOuterSectors();
1246 //LoadInnerSectors();
1252 Int_t AliTPCtrackerMI::LoadClusters()
1255 // load clusters to the memory
1256 AliTPCClustersRow *clrow= new AliTPCClustersRow("AliTPCclusterMI");
1258 // TTree * tree = fClustersArray.GetTree();
1260 TTree * tree = fInput;
1261 TBranch * br = tree->GetBranch("Segment");
1262 br->SetAddress(&clrow);
1264 Int_t j=Int_t(tree->GetEntries());
1265 for (Int_t i=0; i<j; i++) {
1269 fkParam->AdjustSectorRow(clrow->GetID(),sec,row);
1270 for (Int_t icl=0; icl<clrow->GetArray()->GetEntriesFast(); icl++){
1271 Transform((AliTPCclusterMI*)(clrow->GetArray()->At(icl)));
1274 AliTPCtrackerRow * tpcrow=0;
1277 tpcrow = &(fInnerSec[sec%fkNIS][row]);
1281 tpcrow = &(fOuterSec[(sec-fkNIS*2)%fkNOS][row]);
1282 left = (sec-fkNIS*2)/fkNOS;
1285 tpcrow->SetN1(clrow->GetArray()->GetEntriesFast());
1286 tpcrow->SetClusters1(new AliTPCclusterMI[tpcrow->GetN1()]);
1287 for (Int_t k=0;k<tpcrow->GetN1();++k)
1288 tpcrow->SetCluster1(k, *(AliTPCclusterMI*)(clrow->GetArray()->At(k)));
1291 tpcrow->SetN2(clrow->GetArray()->GetEntriesFast());
1292 tpcrow->SetClusters2(new AliTPCclusterMI[tpcrow->GetN2()]);
1293 for (Int_t k=0;k<tpcrow->GetN2();++k)
1294 tpcrow->SetCluster2(k,*(AliTPCclusterMI*)(clrow->GetArray()->At(k)));
1305 void AliTPCtrackerMI::UnloadClusters()
1308 // unload clusters from the memory
1310 Int_t nrows = fOuterSec->GetNRows();
1311 for (Int_t sec = 0;sec<fkNOS;sec++)
1312 for (Int_t row = 0;row<nrows;row++){
1313 AliTPCtrackerRow* tpcrow = &(fOuterSec[sec%fkNOS][row]);
1315 // if (tpcrow->fClusters1) delete []tpcrow->fClusters1;
1316 // if (tpcrow->fClusters2) delete []tpcrow->fClusters2;
1318 tpcrow->ResetClusters();
1321 nrows = fInnerSec->GetNRows();
1322 for (Int_t sec = 0;sec<fkNIS;sec++)
1323 for (Int_t row = 0;row<nrows;row++){
1324 AliTPCtrackerRow* tpcrow = &(fInnerSec[sec%fkNIS][row]);
1326 // if (tpcrow->fClusters1) delete []tpcrow->fClusters1;
1327 //if (tpcrow->fClusters2) delete []tpcrow->fClusters2;
1329 tpcrow->ResetClusters();
1335 void AliTPCtrackerMI::FillClusterArray(TObjArray* array) const{
1337 // Filling cluster to the array - For visualization purposes
1340 nrows = fOuterSec->GetNRows();
1341 for (Int_t sec = 0;sec<fkNOS;sec++)
1342 for (Int_t row = 0;row<nrows;row++){
1343 AliTPCtrackerRow* tpcrow = &(fOuterSec[sec%fkNOS][row]);
1344 if (!tpcrow) continue;
1345 for (Int_t icl = 0;icl<tpcrow->GetN();icl++){
1346 array->AddLast((TObject*)((*tpcrow)[icl]));
1349 nrows = fInnerSec->GetNRows();
1350 for (Int_t sec = 0;sec<fkNIS;sec++)
1351 for (Int_t row = 0;row<nrows;row++){
1352 AliTPCtrackerRow* tpcrow = &(fInnerSec[sec%fkNIS][row]);
1353 if (!tpcrow) continue;
1354 for (Int_t icl = 0;icl<tpcrow->GetN();icl++){
1355 array->AddLast((TObject*)(*tpcrow)[icl]);
1361 void AliTPCtrackerMI::Transform(AliTPCclusterMI * cluster){
1365 AliTPCcalibDB * calibDB = AliTPCcalibDB::Instance();
1366 AliTPCTransform *transform = calibDB->GetTransform() ;
1368 AliFatal("Tranformations not in calibDB");
1371 transform->SetCurrentRecoParam((AliTPCRecoParam*)AliTPCReconstructor::GetRecoParam());
1372 Double_t x[3]={cluster->GetRow(),cluster->GetPad(),cluster->GetTimeBin()};
1373 Int_t i[1]={cluster->GetDetector()};
1374 transform->Transform(x,i,0,1);
1375 // if (cluster->GetDetector()%36>17){
1380 // in debug mode check the transformation
1382 if (AliTPCReconstructor::StreamLevel()>2) {
1384 cluster->GetGlobalXYZ(gx);
1385 Int_t event = (fEvent==NULL)? 0: fEvent->GetEventNumberInFile();
1386 TTreeSRedirector &cstream = *fDebugStreamer;
1387 cstream<<"Transform"<<
1398 cluster->SetX(x[0]);
1399 cluster->SetY(x[1]);
1400 cluster->SetZ(x[2]);
1405 //if (!fkParam->IsGeoRead()) fkParam->ReadGeoMatrices();
1406 if (AliTPCReconstructor::GetRecoParam()->GetUseSectorAlignment() && (!calibDB->HasAlignmentOCDB())){
1407 TGeoHMatrix *mat = fkParam->GetClusterMatrix(cluster->GetDetector());
1409 Double_t pos[3]= {cluster->GetX(),cluster->GetY(),cluster->GetZ()};
1410 Double_t posC[3]={cluster->GetX(),cluster->GetY(),cluster->GetZ()};
1411 if (mat) mat->LocalToMaster(pos,posC);
1413 // chack Loading of Geo matrices from GeoManager - TEMPORARY FIX
1415 cluster->SetX(posC[0]);
1416 cluster->SetY(posC[1]);
1417 cluster->SetZ(posC[2]);
1421 //_____________________________________________________________________________
1422 Int_t AliTPCtrackerMI::LoadOuterSectors() {
1423 //-----------------------------------------------------------------
1424 // This function fills outer TPC sectors with clusters.
1425 //-----------------------------------------------------------------
1426 Int_t nrows = fOuterSec->GetNRows();
1428 for (Int_t sec = 0;sec<fkNOS;sec++)
1429 for (Int_t row = 0;row<nrows;row++){
1430 AliTPCtrackerRow* tpcrow = &(fOuterSec[sec%fkNOS][row]);
1431 Int_t sec2 = sec+2*fkNIS;
1433 Int_t ncl = tpcrow->GetN1();
1435 AliTPCclusterMI *c= (tpcrow->GetCluster1(ncl));
1436 index=(((sec2<<8)+row)<<16)+ncl;
1437 tpcrow->InsertCluster(c,index);
1440 ncl = tpcrow->GetN2();
1442 AliTPCclusterMI *c= (tpcrow->GetCluster2(ncl));
1443 index=((((sec2+fkNOS)<<8)+row)<<16)+ncl;
1444 tpcrow->InsertCluster(c,index);
1447 // write indexes for fast acces
1449 for (Int_t i=0;i<510;i++)
1450 tpcrow->SetFastCluster(i,-1);
1451 for (Int_t i=0;i<tpcrow->GetN();i++){
1452 Int_t zi = Int_t((*tpcrow)[i]->GetZ()+255.);
1453 tpcrow->SetFastCluster(zi,i); // write index
1456 for (Int_t i=0;i<510;i++){
1457 if (tpcrow->GetFastCluster(i)<0)
1458 tpcrow->SetFastCluster(i,last);
1460 last = tpcrow->GetFastCluster(i);
1469 //_____________________________________________________________________________
1470 Int_t AliTPCtrackerMI::LoadInnerSectors() {
1471 //-----------------------------------------------------------------
1472 // This function fills inner TPC sectors with clusters.
1473 //-----------------------------------------------------------------
1474 Int_t nrows = fInnerSec->GetNRows();
1476 for (Int_t sec = 0;sec<fkNIS;sec++)
1477 for (Int_t row = 0;row<nrows;row++){
1478 AliTPCtrackerRow* tpcrow = &(fInnerSec[sec%fkNIS][row]);
1481 Int_t ncl = tpcrow->GetN1();
1483 AliTPCclusterMI *c= (tpcrow->GetCluster1(ncl));
1484 index=(((sec<<8)+row)<<16)+ncl;
1485 tpcrow->InsertCluster(c,index);
1488 ncl = tpcrow->GetN2();
1490 AliTPCclusterMI *c= (tpcrow->GetCluster2(ncl));
1491 index=((((sec+fkNIS)<<8)+row)<<16)+ncl;
1492 tpcrow->InsertCluster(c,index);
1495 // write indexes for fast acces
1497 for (Int_t i=0;i<510;i++)
1498 tpcrow->SetFastCluster(i,-1);
1499 for (Int_t i=0;i<tpcrow->GetN();i++){
1500 Int_t zi = Int_t((*tpcrow)[i]->GetZ()+255.);
1501 tpcrow->SetFastCluster(zi,i); // write index
1504 for (Int_t i=0;i<510;i++){
1505 if (tpcrow->GetFastCluster(i)<0)
1506 tpcrow->SetFastCluster(i,last);
1508 last = tpcrow->GetFastCluster(i);
1520 //_________________________________________________________________________
1521 AliTPCclusterMI *AliTPCtrackerMI::GetClusterMI(Int_t index) const {
1522 //--------------------------------------------------------------------
1523 // Return pointer to a given cluster
1524 //--------------------------------------------------------------------
1525 if (index<0) return 0; // no cluster
1526 Int_t sec=(index&0xff000000)>>24;
1527 Int_t row=(index&0x00ff0000)>>16;
1528 Int_t ncl=(index&0x00007fff)>>00;
1530 const AliTPCtrackerRow * tpcrow=0;
1531 AliTPCclusterMI * clrow =0;
1533 if (sec<0 || sec>=fkNIS*4) {
1534 AliWarning(Form("Wrong sector %d",sec));
1539 AliTPCtrackerSector& tracksec = fInnerSec[sec%fkNIS];
1540 if (tracksec.GetNRows()<=row) return 0;
1541 tpcrow = &(tracksec[row]);
1542 if (tpcrow==0) return 0;
1545 if (tpcrow->GetN1()<=ncl) return 0;
1546 clrow = tpcrow->GetClusters1();
1549 if (tpcrow->GetN2()<=ncl) return 0;
1550 clrow = tpcrow->GetClusters2();
1554 AliTPCtrackerSector& tracksec = fOuterSec[(sec-fkNIS*2)%fkNOS];
1555 if (tracksec.GetNRows()<=row) return 0;
1556 tpcrow = &(tracksec[row]);
1557 if (tpcrow==0) return 0;
1559 if (sec-2*fkNIS<fkNOS) {
1560 if (tpcrow->GetN1()<=ncl) return 0;
1561 clrow = tpcrow->GetClusters1();
1564 if (tpcrow->GetN2()<=ncl) return 0;
1565 clrow = tpcrow->GetClusters2();
1569 return &(clrow[ncl]);
1575 Int_t AliTPCtrackerMI::FollowToNext(AliTPCseed& t, Int_t nr) {
1576 //-----------------------------------------------------------------
1577 // This function tries to find a track prolongation to next pad row
1578 //-----------------------------------------------------------------
1580 Double_t x= GetXrow(nr), ymax=GetMaxY(nr);
1583 AliTPCclusterMI *cl=0;
1584 Int_t tpcindex= t.GetClusterIndex2(nr);
1586 // update current shape info every 5 pad-row
1587 // if ( (nr%5==0) || t.GetNumberOfClusters()<2 || (t.fCurrentSigmaY2<0.0001) ){
1591 if (fIteration>0 && tpcindex>=-1){ //if we have already clusters
1593 if (tpcindex==-1) return 0; //track in dead zone
1595 cl = t.GetClusterPointer(nr);
1596 if ( (cl==0) ) cl = GetClusterMI(tpcindex);
1597 t.SetCurrentClusterIndex1(tpcindex);
1600 Int_t relativesector = ((tpcindex&0xff000000)>>24)%18; // if previously accepted cluster in different sector
1601 Float_t angle = relativesector*fSectors->GetAlpha()+fSectors->GetAlphaShift();
1603 if (angle<-TMath::Pi()) angle += 2*TMath::Pi();
1604 if (angle>=TMath::Pi()) angle -= 2*TMath::Pi();
1606 if (TMath::Abs(angle-t.GetAlpha())>0.001){
1607 Double_t rotation = angle-t.GetAlpha();
1608 t.SetRelativeSector(relativesector);
1609 if (!t.Rotate(rotation)) return 0;
1611 if (!t.PropagateTo(x)) return 0;
1613 t.SetCurrentCluster(cl);
1615 Int_t accept = AcceptCluster(&t,t.GetCurrentCluster());
1616 if ((tpcindex&0x8000)==0) accept =0;
1618 //if founded cluster is acceptible
1619 if (cl->IsUsed(11)) { // id cluster is shared inrease uncertainty
1620 t.SetErrorY2(t.GetErrorY2()+0.03);
1621 t.SetErrorZ2(t.GetErrorZ2()+0.03);
1622 t.SetErrorY2(t.GetErrorY2()*3);
1623 t.SetErrorZ2(t.GetErrorZ2()*3);
1625 t.SetNFoundable(t.GetNFoundable()+1);
1626 UpdateTrack(&t,accept);
1631 if (TMath::Abs(t.GetSnp())>AliTPCReconstructor::GetMaxSnpTracker()) return 0; // cut on angle
1632 if (fIteration>1 && IsFindable(t)){
1633 // not look for new cluster during refitting
1634 t.SetNFoundable(t.GetNFoundable()+1);
1639 // if (TMath::Abs(t.GetSnp())>0.95 || TMath::Abs(x*t.GetC()-t.GetEta())>0.95) return 0;// patch 28 fev 06
1640 Double_t y=t.GetYat(x);
1641 if (TMath::Abs(y)>ymax){
1643 t.SetRelativeSector((t.GetRelativeSector()+1) % fN);
1644 if (!t.Rotate(fSectors->GetAlpha()))
1646 } else if (y <-ymax) {
1647 t.SetRelativeSector((t.GetRelativeSector()-1+fN) % fN);
1648 if (!t.Rotate(-fSectors->GetAlpha()))
1654 if (!t.PropagateTo(x)) {
1655 if (fIteration==0) t.SetRemoval(10);
1659 Double_t z=t.GetZ();
1662 if (!IsActive(t.GetRelativeSector(),nr)) {
1664 t.SetClusterIndex2(nr,-1);
1667 //AliInfo(Form("A - Sector%d phi %f - alpha %f", t.fRelativeSector,y/x, t.GetAlpha()));
1668 Bool_t isActive = IsActive(t.GetRelativeSector(),nr);
1669 Bool_t isActive2 = (nr>=fInnerSec->GetNRows()) ? fOuterSec[t.GetRelativeSector()][nr-fInnerSec->GetNRows()].GetN()>0:fInnerSec[t.GetRelativeSector()][nr].GetN()>0;
1671 if (!isActive || !isActive2) return 0;
1673 const AliTPCtrackerRow &krow=GetRow(t.GetRelativeSector(),nr);
1674 if ( (t.GetSigmaY2()<0) || t.GetSigmaZ2()<0) return 0;
1676 Double_t roadz = 1.;
1678 if (TMath::Abs(TMath::Abs(y)-ymax)<krow.GetDeadZone()){
1680 t.SetClusterIndex2(nr,-1);
1686 // if (TMath::Abs(z)<(AliTPCReconstructor::GetCtgRange()*x+10) && TMath::Abs(z)<fkParam->GetZLength(0) && (TMath::Abs(t.GetSnp())<AliTPCReconstructor::GetMaxSnpTracker()))
1687 t.SetNFoundable(t.GetNFoundable()+1);
1693 // cl = krow.FindNearest2(y+10.,z,roady,roadz,index);
1694 cl = krow.FindNearest2(y,z,roady,roadz,index);
1695 if (cl) t.SetCurrentClusterIndex1(krow.GetIndex(index));
1698 t.SetCurrentCluster(cl);
1700 if (fIteration==2&&cl->IsUsed(10)) return 0;
1701 Int_t accept = AcceptCluster(&t,t.GetCurrentCluster());
1702 if (fIteration==2&&cl->IsUsed(11)) {
1703 t.SetErrorY2(t.GetErrorY2()+0.03);
1704 t.SetErrorZ2(t.GetErrorZ2()+0.03);
1705 t.SetErrorY2(t.GetErrorY2()*3);
1706 t.SetErrorZ2(t.GetErrorZ2()*3);
1709 if (t.fCurrentCluster->IsUsed(10)){
1714 if (t.fNShared>0.7*t.GetNumberOfClusters()) {
1720 if (accept<3) UpdateTrack(&t,accept);
1723 if ( fIteration==0 && t.GetNFoundable()*0.5 > t.GetNumberOfClusters()) t.SetRemoval(10);
1731 //_________________________________________________________________________
1732 Bool_t AliTPCtrackerMI::GetTrackPoint(Int_t index, AliTrackPoint &p ) const
1734 // Get track space point by index
1735 // return false in case the cluster doesn't exist
1736 AliTPCclusterMI *cl = GetClusterMI(index);
1737 if (!cl) return kFALSE;
1738 Int_t sector = (index&0xff000000)>>24;
1739 // Int_t row = (index&0x00ff0000)>>16;
1741 // xyz[0] = fkParam->GetPadRowRadii(sector,row);
1742 xyz[0] = cl->GetX();
1743 xyz[1] = cl->GetY();
1744 xyz[2] = cl->GetZ();
1746 fkParam->AdjustCosSin(sector,cos,sin);
1747 Float_t x = cos*xyz[0]-sin*xyz[1];
1748 Float_t y = cos*xyz[1]+sin*xyz[0];
1750 Float_t sigmaY2 = 0.027*cl->GetSigmaY2();
1751 if (sector < fkParam->GetNInnerSector()) sigmaY2 *= 2.07;
1752 Float_t sigmaZ2 = 0.066*cl->GetSigmaZ2();
1753 if (sector < fkParam->GetNInnerSector()) sigmaZ2 *= 1.77;
1754 cov[0] = sin*sin*sigmaY2;
1755 cov[1] = -sin*cos*sigmaY2;
1757 cov[3] = cos*cos*sigmaY2;
1760 p.SetXYZ(x,y,xyz[2],cov);
1761 AliGeomManager::ELayerID iLayer;
1763 if (sector < fkParam->GetNInnerSector()) {
1764 iLayer = AliGeomManager::kTPC1;
1768 iLayer = AliGeomManager::kTPC2;
1769 idet = sector - fkParam->GetNInnerSector();
1771 UShort_t volid = AliGeomManager::LayerToVolUID(iLayer,idet);
1772 p.SetVolumeID(volid);
1778 Int_t AliTPCtrackerMI::UpdateClusters(AliTPCseed& t, Int_t nr) {
1779 //-----------------------------------------------------------------
1780 // This function tries to find a track prolongation to next pad row
1781 //-----------------------------------------------------------------
1782 t.SetCurrentCluster(0);
1783 t.SetCurrentClusterIndex1(0);
1785 Double_t xt=t.GetX();
1786 Int_t row = GetRowNumber(xt)-1;
1787 Double_t ymax= GetMaxY(nr);
1789 if (row < nr) return 1; // don't prolongate if not information until now -
1790 // if (TMath::Abs(t.GetSnp())>0.9 && t.GetNumberOfClusters()>40. && fIteration!=2) {
1792 // return 0; // not prolongate strongly inclined tracks
1794 // if (TMath::Abs(t.GetSnp())>0.95) {
1796 // return 0; // not prolongate strongly inclined tracks
1797 // }// patch 28 fev 06
1799 Double_t x= GetXrow(nr);
1801 //t.PropagateTo(x+0.02);
1802 //t.PropagateTo(x+0.01);
1803 if (!t.PropagateTo(x)){
1810 if (TMath::Abs(y)>ymax){
1812 t.SetRelativeSector((t.GetRelativeSector()+1) % fN);
1813 if (!t.Rotate(fSectors->GetAlpha()))
1815 } else if (y <-ymax) {
1816 t.SetRelativeSector((t.GetRelativeSector()-1+fN) % fN);
1817 if (!t.Rotate(-fSectors->GetAlpha()))
1820 // if (!t.PropagateTo(x)){
1827 if (TMath::Abs(t.GetSnp())>AliTPCReconstructor::GetMaxSnpTracker()) return 0;
1829 if (!IsActive(t.GetRelativeSector(),nr)) {
1831 t.SetClusterIndex2(nr,-1);
1834 //AliInfo(Form("A - Sector%d phi %f - alpha %f", t.fRelativeSector,y/x, t.GetAlpha()));
1836 AliTPCtrackerRow &krow=GetRow(t.GetRelativeSector(),nr);
1838 if (TMath::Abs(TMath::Abs(y)-ymax)<krow.GetDeadZone()){
1840 t.SetClusterIndex2(nr,-1);
1846 // if (TMath::Abs(t.GetZ())<(AliTPCReconstructor::GetCtgRange()*t.GetX()+10) && (TMath::Abs(t.GetSnp())<AliTPCReconstructor::GetMaxSnpTracker()))
1847 if (IsFindable(t)) t.SetNFoundable(t.GetNFoundable()+1);
1853 if ( (nr%2==0) || t.GetNumberOfClusters()<2){
1854 // t.fCurrentSigmaY = GetSigmaY(&t);
1855 //t.fCurrentSigmaZ = GetSigmaZ(&t);
1859 AliTPCclusterMI *cl=0;
1862 Double_t roady = 1.;
1863 Double_t roadz = 1.;
1867 index = t.GetClusterIndex2(nr);
1868 if ( (index>0) && (index&0x8000)==0){
1869 cl = t.GetClusterPointer(nr);
1870 if ( (cl==0) && (index>0)) cl = GetClusterMI(index);
1871 t.SetCurrentClusterIndex1(index);
1873 t.SetCurrentCluster(cl);
1879 // if (index<0) return 0;
1880 UInt_t uindex = TMath::Abs(index);
1883 //cl = krow.FindNearest2(y+10,z,roady,roadz,uindex);
1884 cl = krow.FindNearest2(y,z,roady,roadz,uindex);
1887 if (cl) t.SetCurrentClusterIndex1(krow.GetIndex(uindex));
1888 t.SetCurrentCluster(cl);
1894 Int_t AliTPCtrackerMI::FollowToNextCluster(AliTPCseed & t, Int_t nr) {
1895 //-----------------------------------------------------------------
1896 // This function tries to find a track prolongation to next pad row
1897 //-----------------------------------------------------------------
1899 //update error according neighborhoud
1901 if (t.GetCurrentCluster()) {
1903 Int_t accept = AcceptCluster(&t,t.GetCurrentCluster());
1905 if (t.GetCurrentCluster()->IsUsed(10)){
1910 t.SetNShared(t.GetNShared()+1);
1911 if (t.GetNShared()>0.7*t.GetNumberOfClusters()) {
1916 if (fIteration>0) accept = 0;
1917 if (accept<3) UpdateTrack(&t,accept);
1921 if ( ( (t.GetSigmaY2()+t.GetSigmaZ2())>0.16)&& t.GetNumberOfClusters()>18) t.SetRemoval(10);
1922 if ( t.GetChi2()/t.GetNumberOfClusters()>6 &&t.GetNumberOfClusters()>18) t.SetRemoval(10);
1924 if (( (t.GetNFoundable()*0.5 > t.GetNumberOfClusters()) || t.GetNoCluster()>15)) t.SetRemoval(10);
1932 //_____________________________________________________________________________
1933 Int_t AliTPCtrackerMI::FollowProlongation(AliTPCseed& t, Int_t rf, Int_t step) {
1934 //-----------------------------------------------------------------
1935 // This function tries to find a track prolongation.
1936 //-----------------------------------------------------------------
1937 Double_t xt=t.GetX();
1939 Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
1940 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
1941 if (alpha < 0. ) alpha += 2.*TMath::Pi();
1943 t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN);
1945 Int_t first = GetRowNumber(xt)-1;
1946 for (Int_t nr= first; nr>=rf; nr-=step) {
1948 if (t.GetKinkIndexes()[0]>0){
1949 for (Int_t i=0;i<3;i++){
1950 Int_t index = t.GetKinkIndexes()[i];
1951 if (index==0) break;
1952 if (index<0) continue;
1954 AliKink * kink = (AliKink*)fEvent->GetKink(index-1);
1956 printf("PROBLEM\n");
1959 Int_t kinkrow = kink->GetTPCRow0()+2+Int_t(0.5/(0.05+kink->GetAngle(2)));
1961 AliExternalTrackParam paramd(t);
1962 kink->SetDaughter(paramd);
1963 kink->SetStatus(2,5);
1970 if (nr==80) t.UpdateReference();
1971 if (nr<fInnerSec->GetNRows())
1972 fSectors = fInnerSec;
1974 fSectors = fOuterSec;
1975 if (FollowToNext(t,nr)==0)
1988 Int_t AliTPCtrackerMI::FollowBackProlongation(AliTPCseed& t, Int_t rf) {
1989 //-----------------------------------------------------------------
1990 // This function tries to find a track prolongation.
1991 //-----------------------------------------------------------------
1993 Double_t xt=t.GetX();
1994 Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
1995 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
1996 if (alpha < 0. ) alpha += 2.*TMath::Pi();
1997 t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN);
1999 Int_t first = t.GetFirstPoint();
2000 if (first<GetRowNumber(xt)+1) first = GetRowNumber(xt)+1;
2002 if (first<0) first=0;
2003 for (Int_t nr=first; nr<=rf; nr++) {
2004 // if ( (TMath::Abs(t.GetSnp())>0.95)) break;//patch 28 fev 06
2005 if (t.GetKinkIndexes()[0]<0){
2006 for (Int_t i=0;i<3;i++){
2007 Int_t index = t.GetKinkIndexes()[i];
2008 if (index==0) break;
2009 if (index>0) continue;
2010 index = TMath::Abs(index);
2011 AliKink * kink = (AliKink*)fEvent->GetKink(index-1);
2013 printf("PROBLEM\n");
2016 Int_t kinkrow = kink->GetTPCRow0()-2-Int_t(0.5/(0.05+kink->GetAngle(2)));
2018 AliExternalTrackParam paramm(t);
2019 kink->SetMother(paramm);
2020 kink->SetStatus(2,1);
2027 if (nr<fInnerSec->GetNRows())
2028 fSectors = fInnerSec;
2030 fSectors = fOuterSec;
2041 Float_t AliTPCtrackerMI::OverlapFactor(AliTPCseed * s1, AliTPCseed * s2, Int_t &sum1, Int_t & sum2)
2049 Float_t dz2 =(s1->GetZ() - s2->GetZ());
2052 Float_t dy2 =TMath::Abs((s1->GetY() - s2->GetY()));
2054 Float_t distance = TMath::Sqrt(dz2+dy2);
2055 if (distance>4.) return 0; // if there are far away - not overlap - to reduce combinatorics
2058 Int_t firstpoint = TMath::Min(s1->GetFirstPoint(),s2->GetFirstPoint());
2059 Int_t lastpoint = TMath::Max(s1->GetLastPoint(),s2->GetLastPoint());
2064 if (firstpoint>lastpoint) {
2065 firstpoint =lastpoint;
2070 for (Int_t i=firstpoint-1;i<lastpoint+1;i++){
2071 if (s1->GetClusterIndex2(i)>0) sum1++;
2072 if (s2->GetClusterIndex2(i)>0) sum2++;
2073 if (s1->GetClusterIndex2(i)==s2->GetClusterIndex2(i) && s1->GetClusterIndex2(i)>0) {
2077 if (sum<5) return 0;
2079 Float_t summin = TMath::Min(sum1+1,sum2+1);
2080 Float_t ratio = (sum+1)/Float_t(summin);
2084 void AliTPCtrackerMI::SignShared(AliTPCseed * s1, AliTPCseed * s2)
2088 Float_t thetaCut = 0.2;//+10.*TMath::Sqrt(s1->GetSigmaTglZ()+ s2->GetSigmaTglZ());
2089 if (TMath::Abs(s1->GetTgl()-s2->GetTgl())>thetaCut) return;
2090 Float_t minCl = TMath::Min(s1->GetNumberOfClusters(),s2->GetNumberOfClusters());
2091 Int_t cutN0 = TMath::Max(5,TMath::Nint(0.1*minCl));
2096 //Int_t firstpoint = TMath::Max(s1->GetFirstPoint(),s2->GetFirstPoint());
2097 //Int_t lastpoint = TMath::Min(s1->GetLastPoint(),s2->GetLastPoint());
2098 Int_t firstpoint = 0;
2099 Int_t lastpoint = 160;
2101 // if (firstpoint>=lastpoint-5) return;;
2103 for (Int_t i=firstpoint;i<lastpoint;i++){
2104 // if ( (s1->GetClusterIndex2(i)&0xFFFF8FFF)==(s2->GetClusterIndex2(i)&0xFFFF8FFF) && s1->GetClusterIndex2(i)>0) {
2105 if ( (s1->GetClusterIndex2(i))==(s2->GetClusterIndex2(i)) && s1->GetClusterIndex2(i)>0) {
2107 s1->SetSharedMapBit(i, kTRUE);
2108 s2->SetSharedMapBit(i, kTRUE);
2110 if (s1->GetClusterIndex2(i)>0)
2111 s1->SetClusterMapBit(i, kTRUE);
2113 if (sumshared>cutN0){
2116 for (Int_t i=firstpoint;i<lastpoint;i++){
2117 // if ( (s1->GetClusterIndex2(i)&0xFFFF8FFF)==(s2->GetClusterIndex2(i)&0xFFFF8FFF) && s1->GetClusterIndex2(i)>0) {
2118 if ( (s1->GetClusterIndex2(i))==(s2->GetClusterIndex2(i)) && s1->GetClusterIndex2(i)>0) {
2119 AliTPCTrackerPoint *p1 = s1->GetTrackPoint(i);
2120 AliTPCTrackerPoint *p2 = s2->GetTrackPoint(i);;
2121 if (s1->IsActive()&&s2->IsActive()){
2122 p1->SetShared(kTRUE);
2123 p2->SetShared(kTRUE);
2129 if (sumshared>cutN0){
2130 for (Int_t i=0;i<4;i++){
2131 if (s1->GetOverlapLabel(3*i)==0){
2132 s1->SetOverlapLabel(3*i, s2->GetLabel());
2133 s1->SetOverlapLabel(3*i+1,sumshared);
2134 s1->SetOverlapLabel(3*i+2,s2->GetUniqueID());
2138 for (Int_t i=0;i<4;i++){
2139 if (s2->GetOverlapLabel(3*i)==0){
2140 s2->SetOverlapLabel(3*i, s1->GetLabel());
2141 s2->SetOverlapLabel(3*i+1,sumshared);
2142 s2->SetOverlapLabel(3*i+2,s1->GetUniqueID());
2149 void AliTPCtrackerMI::SignShared(TObjArray * arr)
2152 //sort trackss according sectors
2154 for (Int_t i=0; i<arr->GetEntriesFast(); i++) {
2155 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2157 //if (pt) RotateToLocal(pt);
2161 arr->Sort(); // sorting according relative sectors
2162 arr->Expand(arr->GetEntries());
2165 Int_t nseed=arr->GetEntriesFast();
2166 for (Int_t i=0; i<nseed; i++) {
2167 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2169 for (Int_t j=0;j<12;j++){
2170 pt->SetOverlapLabel(j,0);
2173 for (Int_t i=0; i<nseed; i++) {
2174 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2176 if (pt->GetRemoval()>10) continue;
2177 for (Int_t j=i+1; j<nseed; j++){
2178 AliTPCseed *pt2=(AliTPCseed*)arr->UncheckedAt(j);
2179 if (TMath::Abs(pt->GetRelativeSector()-pt2->GetRelativeSector())>1) continue;
2181 if (pt2->GetRemoval()<=10) {
2182 //if ( TMath::Abs(pt->GetRelativeSector()-pt2->GetRelativeSector())>0) break;
2190 void AliTPCtrackerMI::SortTracks(TObjArray * arr, Int_t mode) const
2193 //sort tracks in array according mode criteria
2194 Int_t nseed = arr->GetEntriesFast();
2195 for (Int_t i=0; i<nseed; i++) {
2196 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2207 void AliTPCtrackerMI::RemoveUsed2(TObjArray * arr, Float_t factor1, Float_t factor2, Int_t minimal)
2210 // Loop over all tracks and remove overlaped tracks (with lower quality)
2212 // 1. Unsign clusters
2213 // 2. Sort tracks according quality
2214 // Quality is defined by the number of cluster between first and last points
2216 // 3. Loop over tracks - decreasing quality order
2217 // a.) remove - If the fraction of shared cluster less than factor (1- n or 2)
2218 // b.) remove - If the minimal number of clusters less than minimal and not ITS
2219 // c.) if track accepted - sign clusters
2221 //Called in - AliTPCtrackerMI::Clusters2Tracks()
2222 // - AliTPCtrackerMI::PropagateBack()
2223 // - AliTPCtrackerMI::RefitInward()
2226 // factor1 - factor for constrained
2227 // factor2 - for non constrained tracks
2228 // if (Float_t(shared+1)/Float_t(found+1)>factor) - DELETE
2232 Int_t nseed = arr->GetEntriesFast();
2233 Float_t * quality = new Float_t[nseed];
2234 Int_t * indexes = new Int_t[nseed];
2238 for (Int_t i=0; i<nseed; i++) {
2239 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2244 pt->UpdatePoints(); //select first last max dens points
2245 Float_t * points = pt->GetPoints();
2246 if (points[3]<0.8) quality[i] =-1;
2247 quality[i] = (points[2]-points[0])+pt->GetNumberOfClusters();
2248 //prefer high momenta tracks if overlaps
2249 quality[i] *= TMath::Sqrt(TMath::Abs(pt->Pt())+0.5);
2251 TMath::Sort(nseed,quality,indexes);
2254 for (Int_t itrack=0; itrack<nseed; itrack++) {
2255 Int_t trackindex = indexes[itrack];
2256 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(trackindex);
2259 if (quality[trackindex]<0){
2260 delete arr->RemoveAt(trackindex);
2265 Int_t first = Int_t(pt->GetPoints()[0]);
2266 Int_t last = Int_t(pt->GetPoints()[2]);
2267 Double_t factor = (pt->GetBConstrain()) ? factor1: factor2;
2269 Int_t found,foundable,shared;
2270 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
2271 // pt->GetClusterStatistic(0,160, found, foundable,shared,kFALSE);
2272 Bool_t itsgold =kFALSE;
2275 if (pt->GetESD()->GetITSclusters(dummy)>4) itsgold= kTRUE;
2279 if (Float_t(shared+1)/Float_t(found+1)>factor){
2280 if (pt->GetKinkIndexes()[0]!=0) continue; //don't remove tracks - part of the kinks
2281 if( AliTPCReconstructor::StreamLevel()>3){
2282 TTreeSRedirector &cstream = *fDebugStreamer;
2283 cstream<<"RemoveUsed"<<
2284 "iter="<<fIteration<<
2288 delete arr->RemoveAt(trackindex);
2291 if (pt->GetNumberOfClusters()<50&&(found-0.5*shared)<minimal){ //remove short tracks
2292 if (pt->GetKinkIndexes()[0]!=0) continue; //don't remove tracks - part of the kinks
2293 if( AliTPCReconstructor::StreamLevel()>3){
2294 TTreeSRedirector &cstream = *fDebugStreamer;
2295 cstream<<"RemoveShort"<<
2296 "iter="<<fIteration<<
2300 delete arr->RemoveAt(trackindex);
2306 //if (sharedfactor>0.4) continue;
2307 if (pt->GetKinkIndexes()[0]>0) continue;
2308 //Remove tracks with undefined properties - seems
2309 if (pt->GetSigmaY2()<kAlmost0) continue; // ? what is the origin ?
2311 for (Int_t i=first; i<last; i++) {
2312 Int_t index=pt->GetClusterIndex2(i);
2313 // if (index<0 || index&0x8000 ) continue;
2314 if (index<0 || index&0x8000 ) continue;
2315 AliTPCclusterMI *c= pt->GetClusterPointer(i);
2322 Info("RemoveUsed","\n*****\nNumber of good tracks after shared removal\t%d\n",fNtracks);
2328 void AliTPCtrackerMI::DumpClusters(Int_t iter, TObjArray *trackArray)
2331 // Dump clusters after reco
2332 // signed and unsigned cluster can be visualized
2333 // 1. Unsign all cluster
2334 // 2. Sign all used clusters
2337 Int_t nseed = trackArray->GetEntries();
2338 for (Int_t i=0; i<nseed; i++){
2339 AliTPCseed *pt=(AliTPCseed*)trackArray->UncheckedAt(i);
2343 Bool_t isKink=pt->GetKinkIndex(0)!=0;
2344 for (Int_t j=0; j<160; ++j) {
2345 Int_t index=pt->GetClusterIndex2(j);
2346 if (index<0) continue;
2347 AliTPCclusterMI *c= pt->GetClusterPointer(j);
2349 if (isKink) c->Use(100); // kink
2350 c->Use(10); // by default usage 10
2355 for (Int_t sec=0;sec<fkNIS;sec++){
2356 for (Int_t row=0;row<fInnerSec->GetNRows();row++){
2357 AliTPCclusterMI *cl = fInnerSec[sec][row].GetClusters1();
2358 for (Int_t icl =0;icl< fInnerSec[sec][row].GetN1();icl++){
2359 Float_t gx[3]; cl[icl].GetGlobalXYZ(gx);
2360 (*fDebugStreamer)<<"clDump"<<
2368 cl = fInnerSec[sec][row].GetClusters2();
2369 for (Int_t icl =0;icl< fInnerSec[sec][row].GetN2();icl++){
2370 Float_t gx[3]; cl[icl].GetGlobalXYZ(gx);
2371 (*fDebugStreamer)<<"clDump"<<
2382 for (Int_t sec=0;sec<fkNOS;sec++){
2383 for (Int_t row=0;row<fOuterSec->GetNRows();row++){
2384 AliTPCclusterMI *cl = fOuterSec[sec][row].GetClusters1();
2385 for (Int_t icl =0;icl< fOuterSec[sec][row].GetN1();icl++){
2386 Float_t gx[3]; cl[icl].GetGlobalXYZ(gx);
2387 (*fDebugStreamer)<<"clDump"<<
2395 cl = fOuterSec[sec][row].GetClusters2();
2396 for (Int_t icl =0;icl< fOuterSec[sec][row].GetN2();icl++){
2397 Float_t gx[3]; cl[icl].GetGlobalXYZ(gx);
2398 (*fDebugStreamer)<<"clDump"<<
2410 void AliTPCtrackerMI::UnsignClusters()
2413 // loop over all clusters and unsign them
2416 for (Int_t sec=0;sec<fkNIS;sec++){
2417 for (Int_t row=0;row<fInnerSec->GetNRows();row++){
2418 AliTPCclusterMI *cl = fInnerSec[sec][row].GetClusters1();
2419 for (Int_t icl =0;icl< fInnerSec[sec][row].GetN1();icl++)
2420 // if (cl[icl].IsUsed(10))
2422 cl = fInnerSec[sec][row].GetClusters2();
2423 for (Int_t icl =0;icl< fInnerSec[sec][row].GetN2();icl++)
2424 //if (cl[icl].IsUsed(10))
2429 for (Int_t sec=0;sec<fkNOS;sec++){
2430 for (Int_t row=0;row<fOuterSec->GetNRows();row++){
2431 AliTPCclusterMI *cl = fOuterSec[sec][row].GetClusters1();
2432 for (Int_t icl =0;icl< fOuterSec[sec][row].GetN1();icl++)
2433 //if (cl[icl].IsUsed(10))
2435 cl = fOuterSec[sec][row].GetClusters2();
2436 for (Int_t icl =0;icl< fOuterSec[sec][row].GetN2();icl++)
2437 //if (cl[icl].IsUsed(10))
2446 void AliTPCtrackerMI::SignClusters(const TObjArray * arr, Float_t fnumber, Float_t fdensity)
2449 //sign clusters to be "used"
2451 // snumber and sdensity sign number of sigmas - bellow mean value to be accepted
2452 // loop over "primaries"
2466 Int_t nseed = arr->GetEntriesFast();
2467 for (Int_t i=0; i<nseed; i++) {
2468 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2472 if (!(pt->IsActive())) continue;
2473 Float_t dens = pt->GetNumberOfClusters()/Float_t(pt->GetNFoundable());
2474 if ( (dens>0.7) && (pt->GetNumberOfClusters()>70)){
2476 sumdens2+= dens*dens;
2477 sumn += pt->GetNumberOfClusters();
2478 sumn2 += pt->GetNumberOfClusters()*pt->GetNumberOfClusters();
2479 Float_t chi2 = pt->GetChi2()/pt->GetNumberOfClusters();
2482 sumchi2 +=chi2*chi2;
2487 Float_t mdensity = 0.9;
2488 Float_t meann = 130;
2489 Float_t meanchi = 1;
2490 Float_t sdensity = 0.1;
2491 Float_t smeann = 10;
2492 Float_t smeanchi =0.4;
2496 mdensity = sumdens/sum;
2498 meanchi = sumchi/sum;
2500 sdensity = sumdens2/sum-mdensity*mdensity;
2502 sdensity = TMath::Sqrt(sdensity);
2506 smeann = sumn2/sum-meann*meann;
2508 smeann = TMath::Sqrt(smeann);
2512 smeanchi = sumchi2/sum - meanchi*meanchi;
2514 smeanchi = TMath::Sqrt(smeanchi);
2520 //REMOVE SHORT DELTAS or tracks going out of sensitive volume of TPC
2522 for (Int_t i=0; i<nseed; i++) {
2523 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2527 if (pt->GetBSigned()) continue;
2528 if (pt->GetBConstrain()) continue;
2529 //if (!(pt->IsActive())) continue;
2531 Int_t found,foundable,shared;
2532 pt->GetClusterStatistic(0,160,found, foundable,shared);
2533 if (shared/float(found)>0.3) {
2534 if (shared/float(found)>0.9 ){
2535 //delete arr->RemoveAt(i);
2540 Bool_t isok =kFALSE;
2541 if ( (pt->GetNShared()/pt->GetNumberOfClusters()<0.5) &&pt->GetNumberOfClusters()>60)
2543 if ((TMath::Abs(1/pt->GetC())<100.) && (pt->GetNShared()/pt->GetNumberOfClusters()<0.7))
2545 if (TMath::Abs(pt->GetZ()/pt->GetX())>1.1)
2547 if ( (TMath::Abs(pt->GetSnp()>0.7) && pt->GetD(0,0)>60.))
2551 for (Int_t j=0; j<160; ++j) {
2552 Int_t index=pt->GetClusterIndex2(j);
2553 if (index<0) continue;
2554 AliTPCclusterMI *c= pt->GetClusterPointer(j);
2556 //if (!(c->IsUsed(10))) c->Use();
2563 Double_t maxchi = meanchi+2.*smeanchi;
2565 for (Int_t i=0; i<nseed; i++) {
2566 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2570 //if (!(pt->IsActive())) continue;
2571 if (pt->GetBSigned()) continue;
2572 Double_t chi = pt->GetChi2()/pt->GetNumberOfClusters();
2573 if (chi>maxchi) continue;
2576 Float_t dens = pt->GetNumberOfClusters()/Float_t(pt->GetNFoundable());
2578 //sign only tracks with enoug big density at the beginning
2580 if ((pt->GetDensityFirst(40)<0.75) && pt->GetNumberOfClusters()<meann) continue;
2583 Double_t mindens = TMath::Max(double(mdensity-sdensity*fdensity*bfactor),0.65);
2584 Double_t minn = TMath::Max(Int_t(meann-fnumber*smeann*bfactor),50);
2586 // if (pt->fBConstrain) mindens = TMath::Max(mdensity-sdensity*fdensity*bfactor,0.65);
2587 if ( (pt->GetRemoval()==10) && (pt->GetSnp()>0.8)&&(dens>mindens))
2590 if ((dens>mindens && pt->GetNumberOfClusters()>minn) && chi<maxchi ){
2591 //Int_t noc=pt->GetNumberOfClusters();
2592 pt->SetBSigned(kTRUE);
2593 for (Int_t j=0; j<160; ++j) {
2595 Int_t index=pt->GetClusterIndex2(j);
2596 if (index<0) continue;
2597 AliTPCclusterMI *c= pt->GetClusterPointer(j);
2599 // if (!(c->IsUsed(10))) c->Use();
2604 // gLastCheck = nseed;
2612 void AliTPCtrackerMI::StopNotActive(const TObjArray * arr, Int_t row0, Float_t th0, Float_t th1, Float_t th2) const
2614 // stop not active tracks
2615 // take th1 as threshold for number of founded to number of foundable on last 10 active rows
2616 // take th2 as threshold for number of founded to number of foundable on last 20 active rows
2617 Int_t nseed = arr->GetEntriesFast();
2619 for (Int_t i=0; i<nseed; i++) {
2620 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2624 if (!(pt->IsActive())) continue;
2625 StopNotActive(pt,row0,th0, th1,th2);
2631 void AliTPCtrackerMI::StopNotActive(AliTPCseed * seed, Int_t row0, Float_t th0, Float_t th1,
2634 // stop not active tracks
2635 // take th1 as threshold for number of founded to number of foundable on last 10 active rows
2636 // take th2 as threshold for number of founded to number of foundable on last 20 active rows
2639 Int_t foundable = 0;
2640 Int_t maxindex = seed->GetLastPoint(); //last foundable row
2641 if (seed->GetNFoundable()*th0 > seed->GetNumberOfClusters()) {
2642 seed->Desactivate(10) ;
2646 for (Int_t i=row0; i<maxindex; i++){
2647 Int_t index = seed->GetClusterIndex2(i);
2648 if (index!=-1) foundable++;
2650 if (foundable<=30) sumgood1++;
2651 if (foundable<=50) {
2658 if (foundable>=30.){
2659 if (sumgood1<(th1*30.)) seed->Desactivate(10);
2662 if (sumgood2<(th2*50.)) seed->Desactivate(10);
2666 Int_t AliTPCtrackerMI::RefitInward(AliESDEvent *event)
2669 // back propagation of ESD tracks
2672 if (!event) return 0;
2673 const Int_t kMaxFriendTracks=2000;
2675 // extract correction object for multiplicity dependence of dEdx
2676 TObjArray * gainCalibArray = AliTPCcalibDB::Instance()->GetTimeGainSplinesRun(event->GetRunNumber());
2678 AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ;
2680 AliFatal("Tranformations not in RefitInward");
2683 transform->SetCurrentRecoParam((AliTPCRecoParam*)AliTPCReconstructor::GetRecoParam());
2684 const AliTPCRecoParam * recoParam = AliTPCcalibDB::Instance()->GetTransform()->GetCurrentRecoParam();
2685 Int_t nContribut = event->GetNumberOfTracks();
2686 TGraphErrors * graphMultDependenceDeDx = 0x0;
2687 if (recoParam && recoParam->GetUseMultiplicityCorrectionDedx() && gainCalibArray) {
2688 if (recoParam->GetUseTotCharge()) {
2689 graphMultDependenceDeDx = (TGraphErrors *) gainCalibArray->FindObject("TGRAPHERRORS_MEANQTOT_MULTIPLICITYDEPENDENCE_BEAM_ALL");
2691 graphMultDependenceDeDx = (TGraphErrors *) gainCalibArray->FindObject("TGRAPHERRORS_MEANQMAX_MULTIPLICITYDEPENDENCE_BEAM_ALL");
2697 //PrepareForProlongation(fSeeds,1);
2698 PropagateForward2(fSeeds);
2699 RemoveUsed2(fSeeds,0.4,0.4,20);
2701 TObjArray arraySeed(fSeeds->GetEntries());
2702 for (Int_t i=0;i<fSeeds->GetEntries();i++) {
2703 arraySeed.AddAt(fSeeds->At(i),i);
2705 SignShared(&arraySeed);
2706 // FindCurling(fSeeds, event,2); // find multi found tracks
2707 FindSplitted(fSeeds, event,2); // find multi found tracks
2708 if (AliTPCReconstructor::StreamLevel()>5) FindMultiMC(fSeeds, fEvent,2); // find multi found tracks
2711 Int_t nseed = fSeeds->GetEntriesFast();
2712 for (Int_t i=0;i<nseed;i++){
2713 AliTPCseed * seed = (AliTPCseed*) fSeeds->UncheckedAt(i);
2714 if (!seed) continue;
2715 if (seed->GetKinkIndex(0)>0) UpdateKinkQualityD(seed); // update quality informations for kinks
2716 AliESDtrack *esd=event->GetTrack(i);
2718 if (seed->GetNumberOfClusters()<60 && seed->GetNumberOfClusters()<(esd->GetTPCclusters(0) -5)*0.8){
2719 AliExternalTrackParam paramIn;
2720 AliExternalTrackParam paramOut;
2721 Int_t ncl = seed->RefitTrack(seed,¶mIn,¶mOut);
2722 if (AliTPCReconstructor::StreamLevel()>2) {
2723 (*fDebugStreamer)<<"RecoverIn"<<
2727 "pout.="<<¶mOut<<
2732 seed->Set(paramIn.GetX(),paramIn.GetAlpha(),paramIn.GetParameter(),paramIn.GetCovariance());
2733 seed->SetNumberOfClusters(ncl);
2737 seed->PropagateTo(fkParam->GetInnerRadiusLow());
2738 seed->UpdatePoints();
2739 AddCovariance(seed);
2741 seed->CookdEdx(0.02,0.6);
2742 CookLabel(seed,0.1); //For comparison only
2743 esd->SetTPCClusterMap(seed->GetClusterMap());
2744 esd->SetTPCSharedMap(seed->GetSharedMap());
2746 if (AliTPCReconstructor::StreamLevel()>1 && seed!=0) {
2747 TTreeSRedirector &cstream = *fDebugStreamer;
2754 if (seed->GetNumberOfClusters()>15){
2755 esd->UpdateTrackParams(seed,AliESDtrack::kTPCrefit);
2756 esd->SetTPCPoints(seed->GetPoints());
2757 esd->SetTPCPointsF(seed->GetNFoundable());
2758 Int_t ndedx = seed->GetNCDEDX(0);
2759 Float_t sdedx = seed->GetSDEDX(0);
2760 Float_t dedx = seed->GetdEdx();
2761 // apply mutliplicity dependent dEdx correction if available
2762 if (graphMultDependenceDeDx) {
2763 Double_t corrGain = AliTPCcalibDButil::EvalGraphConst(graphMultDependenceDeDx, nContribut);
2764 dedx += (1 - corrGain)*50.; // MIP is normalized to 50
2766 esd->SetTPCsignal(dedx, sdedx, ndedx);
2768 // add seed to the esd track in Calib level
2770 Bool_t storeFriend = gRandom->Rndm()<(kMaxFriendTracks)/Float_t(nseed);
2771 if (AliTPCReconstructor::StreamLevel()>0 &&storeFriend){
2772 AliTPCseed * seedCopy = new AliTPCseed(*seed, kTRUE);
2773 esd->AddCalibObject(seedCopy);
2778 //printf("problem\n");
2781 //FindKinks(fSeeds,event);
2782 if (AliTPCReconstructor::StreamLevel()>3) DumpClusters(2,fSeeds);
2783 Info("RefitInward","Number of refitted tracks %d",ntracks);
2788 Int_t AliTPCtrackerMI::PropagateBack(AliESDEvent *event)
2791 // back propagation of ESD tracks
2793 if (!event) return 0;
2797 PropagateBack(fSeeds);
2798 RemoveUsed2(fSeeds,0.4,0.4,20);
2799 //FindCurling(fSeeds, fEvent,1);
2800 FindSplitted(fSeeds, event,1); // find multi found tracks
2801 if (AliTPCReconstructor::StreamLevel()>5) FindMultiMC(fSeeds, fEvent,1); // find multi found tracks
2804 Int_t nseed = fSeeds->GetEntriesFast();
2806 for (Int_t i=0;i<nseed;i++){
2807 AliTPCseed * seed = (AliTPCseed*) fSeeds->UncheckedAt(i);
2808 if (!seed) continue;
2809 if (seed->GetKinkIndex(0)<0) UpdateKinkQualityM(seed); // update quality informations for kinks
2810 seed->UpdatePoints();
2811 AddCovariance(seed);
2812 AliESDtrack *esd=event->GetTrack(i);
2813 if (seed->GetNumberOfClusters()<60 && seed->GetNumberOfClusters()<(esd->GetTPCclusters(0) -5)*0.8){
2814 AliExternalTrackParam paramIn;
2815 AliExternalTrackParam paramOut;
2816 Int_t ncl = seed->RefitTrack(seed,¶mIn,¶mOut);
2817 if (AliTPCReconstructor::StreamLevel()>2) {
2818 (*fDebugStreamer)<<"RecoverBack"<<
2822 "pout.="<<¶mOut<<
2827 seed->Set(paramOut.GetX(),paramOut.GetAlpha(),paramOut.GetParameter(),paramOut.GetCovariance());
2828 seed->SetNumberOfClusters(ncl);
2831 seed->CookdEdx(0.02,0.6);
2832 CookLabel(seed,0.1); //For comparison only
2833 if (seed->GetNumberOfClusters()>15){
2834 esd->UpdateTrackParams(seed,AliESDtrack::kTPCout);
2835 esd->SetTPCPoints(seed->GetPoints());
2836 esd->SetTPCPointsF(seed->GetNFoundable());
2837 Int_t ndedx = seed->GetNCDEDX(0);
2838 Float_t sdedx = seed->GetSDEDX(0);
2839 Float_t dedx = seed->GetdEdx();
2840 esd->SetTPCsignal(dedx, sdedx, ndedx);
2842 Int_t eventnumber = event->GetEventNumberInFile();// patch 28 fev 06
2843 // This is most likely NOT the event number you'd like to use. It has nothing to do with the 'real' event number
2844 if (AliTPCReconstructor::StreamLevel()>1 && esd) {
2845 (*fDebugStreamer)<<"Cback"<<
2848 "EventNrInFile="<<eventnumber<<
2853 if (AliTPCReconstructor::StreamLevel()>3) DumpClusters(1,fSeeds);
2854 //FindKinks(fSeeds,event);
2855 Info("PropagateBack","Number of back propagated tracks %d",ntracks);
2862 void AliTPCtrackerMI::DeleteSeeds()
2867 Int_t nseed = fSeeds->GetEntriesFast();
2868 for (Int_t i=0;i<nseed;i++){
2869 AliTPCseed * seed = (AliTPCseed*)fSeeds->At(i);
2870 if (seed) delete fSeeds->RemoveAt(i);
2877 void AliTPCtrackerMI::ReadSeeds(AliESDEvent *const event, Int_t direction)
2880 //read seeds from the event
2882 Int_t nentr=event->GetNumberOfTracks();
2884 Info("ReadSeeds", "Number of ESD tracks: %d\n", nentr);
2889 fSeeds = new TObjArray(nentr);
2893 for (Int_t i=0; i<nentr; i++) {
2894 AliESDtrack *esd=event->GetTrack(i);
2895 ULong_t status=esd->GetStatus();
2896 if (!(status&AliESDtrack::kTPCin)) continue;
2897 AliTPCtrack t(*esd);
2898 t.SetNumberOfClusters(0);
2899 // AliTPCseed *seed = new AliTPCseed(t,t.GetAlpha());
2900 AliTPCseed *seed = new AliTPCseed(t/*,t.GetAlpha()*/);
2901 seed->SetUniqueID(esd->GetID());
2902 AddCovariance(seed); //add systematic ucertainty
2903 for (Int_t ikink=0;ikink<3;ikink++) {
2904 Int_t index = esd->GetKinkIndex(ikink);
2905 seed->GetKinkIndexes()[ikink] = index;
2906 if (index==0) continue;
2907 index = TMath::Abs(index);
2908 AliESDkink * kink = fEvent->GetKink(index-1);
2909 if (kink&&esd->GetKinkIndex(ikink)<0){
2910 if ((status & AliESDtrack::kTRDrefit) != 0) kink->SetStatus(1,2);
2911 if ((status & AliESDtrack::kITSout) != 0) kink->SetStatus(1,0);
2913 if (kink&&esd->GetKinkIndex(ikink)>0){
2914 if ((status & AliESDtrack::kTRDrefit) != 0) kink->SetStatus(1,6);
2915 if ((status & AliESDtrack::kITSout) != 0) kink->SetStatus(1,4);
2919 if (((status&AliESDtrack::kITSout)==0)&&(direction==1)) seed->ResetCovariance(10.);
2920 if ( direction ==2 &&(status & AliESDtrack::kTRDrefit) == 0 ) seed->ResetCovariance(10.);
2921 //if ( direction ==2 && ((status & AliESDtrack::kTPCout) == 0) ) {
2922 // fSeeds->AddAt(0,i);
2926 if ( direction ==2 &&(status & AliESDtrack::kTRDrefit) > 0 ) {
2927 Double_t par0[5],par1[5],alpha,x;
2928 esd->GetInnerExternalParameters(alpha,x,par0);
2929 esd->GetExternalParameters(x,par1);
2930 Double_t delta1 = TMath::Abs(par0[4]-par1[4])/(0.000000001+TMath::Abs(par0[4]+par1[4]));
2931 Double_t delta2 = TMath::Abs(par0[3]-par1[3]);
2933 if (esd->GetTRDncls()>0) trdchi2 = esd->GetTRDchi2()/esd->GetTRDncls();
2934 //reset covariance if suspicious
2935 if ( (delta1>0.1) || (delta2>0.006) ||trdchi2>7.)
2936 seed->ResetCovariance(10.);
2941 // rotate to the local coordinate system
2943 fSectors=fInnerSec; fN=fkNIS;
2944 Double_t alpha=seed->GetAlpha() - fSectors->GetAlphaShift();
2945 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
2946 if (alpha < 0. ) alpha += 2.*TMath::Pi();
2947 Int_t ns=Int_t(alpha/fSectors->GetAlpha())%fN;
2948 alpha =ns*fSectors->GetAlpha() + fSectors->GetAlphaShift();
2949 if (alpha<-TMath::Pi()) alpha += 2*TMath::Pi();
2950 if (alpha>=TMath::Pi()) alpha -= 2*TMath::Pi();
2951 alpha-=seed->GetAlpha();
2952 if (!seed->Rotate(alpha)) {
2958 if (esd->GetKinkIndex(0)<=0){
2959 for (Int_t irow=0;irow<160;irow++){
2960 Int_t index = seed->GetClusterIndex2(irow);
2963 AliTPCclusterMI * cl = GetClusterMI(index);
2964 seed->SetClusterPointer(irow,cl);
2966 if ((index & 0x8000)==0){
2967 cl->Use(10); // accepted cluster
2969 cl->Use(6); // close cluster not accepted
2972 Info("ReadSeeds","Not found cluster");
2977 fSeeds->AddAt(seed,i);
2983 //_____________________________________________________________________________
2984 void AliTPCtrackerMI::MakeSeeds3(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, Float_t cuts[4],
2985 Float_t deltay, Int_t ddsec) {
2986 //-----------------------------------------------------------------
2987 // This function creates track seeds.
2988 // SEEDING WITH VERTEX CONSTRAIN
2989 //-----------------------------------------------------------------
2990 // cuts[0] - fP4 cut
2991 // cuts[1] - tan(phi) cut
2992 // cuts[2] - zvertex cut
2993 // cuts[3] - fP3 cut
3001 Double_t x[5], c[15];
3002 // Int_t di = i1-i2;
3004 AliTPCseed * seed = new AliTPCseed();
3005 Double_t alpha=fSectors->GetAlpha(), shift=fSectors->GetAlphaShift();
3006 Double_t cs=cos(alpha), sn=sin(alpha);
3008 // Double_t x1 =fOuterSec->GetX(i1);
3009 //Double_t xx2=fOuterSec->GetX(i2);
3011 Double_t x1 =GetXrow(i1);
3012 Double_t xx2=GetXrow(i2);
3014 Double_t x3=GetX(), y3=GetY(), z3=GetZ();
3016 Int_t imiddle = (i2+i1)/2; //middle pad row index
3017 Double_t xm = GetXrow(imiddle); // radius of middle pad-row
3018 const AliTPCtrackerRow& krm=GetRow(sec,imiddle); //middle pad -row
3022 const AliTPCtrackerRow& kr1=GetRow(ns,i1);
3023 Double_t ymax = GetMaxY(i1)-kr1.GetDeadZone()-1.5;
3024 Double_t ymaxm = GetMaxY(imiddle)-kr1.GetDeadZone()-1.5;
3027 // change cut on curvature if it can't reach this layer
3028 // maximal curvature set to reach it
3029 Double_t dvertexmax = TMath::Sqrt((x1-x3)*(x1-x3)+(ymax+5-y3)*(ymax+5-y3));
3030 if (dvertexmax*0.5*cuts[0]>0.85){
3031 cuts[0] = 0.85/(dvertexmax*0.5+1.);
3033 Double_t r2min = 1/(cuts[0]*cuts[0]); //minimal square of radius given by cut
3036 if (deltay>0) ddsec = 0;
3037 // loop over clusters
3038 for (Int_t is=0; is < kr1; is++) {
3040 if (kr1[is]->IsUsed(10)) continue;
3041 Double_t y1=kr1[is]->GetY(), z1=kr1[is]->GetZ();
3042 //if (TMath::Abs(y1)>ymax) continue;
3044 if (deltay>0 && TMath::Abs(ymax-TMath::Abs(y1))> deltay ) continue; // seed only at the edge
3046 // find possible directions
3047 Float_t anglez = (z1-z3)/(x1-x3);
3048 Float_t extraz = z1 - anglez*(x1-xx2); // extrapolated z
3051 //find rotation angles relative to line given by vertex and point 1
3052 Double_t dvertex2 = (x1-x3)*(x1-x3)+(y1-y3)*(y1-y3);
3053 Double_t dvertex = TMath::Sqrt(dvertex2);
3054 Double_t angle13 = TMath::ATan((y1-y3)/(x1-x3));
3055 Double_t cs13 = cos(-angle13), sn13 = sin(-angle13);
3058 // loop over 2 sectors
3064 Double_t dddz1=0; // direction of delta inclination in z axis
3071 for (Int_t dsec = dsec1; dsec<=dsec2;dsec++){
3072 Int_t sec2 = sec + dsec;
3074 // AliTPCtrackerRow& kr2 = fOuterSec[(sec2+fkNOS)%fkNOS][i2];
3075 //AliTPCtrackerRow& kr2m = fOuterSec[(sec2+fkNOS)%fkNOS][imiddle];
3076 AliTPCtrackerRow& kr2 = GetRow((sec2+fkNOS)%fkNOS,i2);
3077 AliTPCtrackerRow& kr2m = GetRow((sec2+fkNOS)%fkNOS,imiddle);
3078 Int_t index1 = TMath::Max(kr2.Find(extraz-0.6-dddz1*TMath::Abs(z1)*0.05)-1,0);
3079 Int_t index2 = TMath::Min(kr2.Find(extraz+0.6+dddz2*TMath::Abs(z1)*0.05)+1,kr2);
3081 // rotation angles to p1-p3
3082 Double_t cs13r = cos(-angle13+dsec*alpha)/dvertex, sn13r = sin(-angle13+dsec*alpha)/dvertex;
3083 Double_t x2, y2, z2;
3085 // Double_t dymax = maxangle*TMath::Abs(x1-xx2);
3088 Double_t dxx0 = (xx2-x3)*cs13r;
3089 Double_t dyy0 = (xx2-x3)*sn13r;
3090 for (Int_t js=index1; js < index2; js++) {
3091 const AliTPCclusterMI *kcl = kr2[js];
3092 if (kcl->IsUsed(10)) continue;
3094 //calcutate parameters
3096 Double_t yy0 = dyy0 +(kcl->GetY()-y3)*cs13r;
3098 if (TMath::Abs(yy0)<0.000001) continue;
3099 Double_t xx0 = dxx0 -(kcl->GetY()-y3)*sn13r;
3100 Double_t y0 = 0.5*(xx0*xx0+yy0*yy0-xx0)/yy0;
3101 Double_t r02 = (0.25+y0*y0)*dvertex2;
3102 //curvature (radius) cut
3103 if (r02<r2min) continue;
3107 Double_t c0 = 1/TMath::Sqrt(r02);
3111 //Double_t dfi0 = 2.*TMath::ASin(dvertex*c0*0.5);
3112 //Double_t dfi1 = 2.*TMath::ASin(TMath::Sqrt(yy0*yy0+(1-xx0)*(1-xx0))*dvertex*c0*0.5);
3113 Double_t dfi0 = 2.*AliTPCFastMath::FastAsin(dvertex*c0*0.5);
3114 Double_t dfi1 = 2.*AliTPCFastMath::FastAsin(TMath::Sqrt(yy0*yy0+(1-xx0)*(1-xx0))*dvertex*c0*0.5);
3117 Double_t z0 = kcl->GetZ();
3118 Double_t zzzz2 = z1-(z1-z3)*dfi1/dfi0;
3119 if (TMath::Abs(zzzz2-z0)>0.5) continue;
3122 Double_t dip = (z1-z0)*c0/dfi1;
3123 Double_t x0 = (0.5*cs13+y0*sn13)*dvertex*c0;
3134 x2= xx2*cs-y2*sn*dsec;
3135 y2=+xx2*sn*dsec+y2*cs;
3145 // do we have cluster at the middle ?
3147 GetProlongation(x1,xm,x,ym,zm);
3149 AliTPCclusterMI * cm=0;
3150 if (TMath::Abs(ym)-ymaxm<0){
3151 cm = krm.FindNearest2(ym,zm,1.0,0.6,dummy);
3152 if ((!cm) || (cm->IsUsed(10))) {
3157 // rotate y1 to system 0
3158 // get state vector in rotated system
3159 Double_t yr1 = (-0.5*sn13+y0*cs13)*dvertex*c0;
3160 Double_t xr2 = x0*cs+yr1*sn*dsec;
3161 Double_t xr[5]={kcl->GetY(),kcl->GetZ(), xr2, dip, c0};
3163 GetProlongation(xx2,xm,xr,ym,zm);
3164 if (TMath::Abs(ym)-ymaxm<0){
3165 cm = kr2m.FindNearest2(ym,zm,1.0,0.6,dummy);
3166 if ((!cm) || (cm->IsUsed(10))) {
3176 dym = ym - cm->GetY();
3177 dzm = zm - cm->GetZ();
3184 Double_t sy1=kr1[is]->GetSigmaY2()*2., sz1=kr1[is]->GetSigmaZ2()*2.;
3185 Double_t sy2=kcl->GetSigmaY2()*2., sz2=kcl->GetSigmaZ2()*2.;
3186 //Double_t sy3=400*3./12., sy=0.1, sz=0.1;
3187 Double_t sy3=25000*x[4]*x[4]+0.1, sy=0.1, sz=0.1;
3188 //Double_t sy3=25000*x[4]*x[4]*60+0.5, sy=0.1, sz=0.1;
3190 Double_t f40=(F1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
3191 Double_t f42=(F1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
3192 Double_t f43=(F1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
3193 Double_t f20=(F2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
3194 Double_t f22=(F2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
3195 Double_t f23=(F2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
3197 Double_t f30=(F3(x1,y1+sy,x2,y2,z1,z2)-x[3])/sy;
3198 Double_t f31=(F3(x1,y1,x2,y2,z1+sz,z2)-x[3])/sz;
3199 Double_t f32=(F3(x1,y1,x2,y2+sy,z1,z2)-x[3])/sy;
3200 Double_t f34=(F3(x1,y1,x2,y2,z1,z2+sz)-x[3])/sz;
3204 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3205 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3206 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3207 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3208 c[13]=f30*sy1*f40+f32*sy2*f42;
3209 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3211 // if (!BuildSeed(kr1[is],kcl,0,x1,x2,x3,x,c)) continue;
3213 UInt_t index=kr1.GetIndex(is);
3214 seed->~AliTPCseed(); // this does not set the pointer to 0...
3215 AliTPCseed *track=new(seed) AliTPCseed(x1, ns*alpha+shift, x, c, index);
3217 track->SetIsSeeding(kTRUE);
3218 track->SetSeed1(i1);
3219 track->SetSeed2(i2);
3220 track->SetSeedType(3);
3224 FollowProlongation(*track, (i1+i2)/2,1);
3225 Int_t foundable,found,shared;
3226 track->GetClusterStatistic((i1+i2)/2,i1, found, foundable, shared, kTRUE);
3227 if ((found<0.55*foundable) || shared>0.5*found || (track->GetSigmaY2()+track->GetSigmaZ2())>0.5){
3229 seed->~AliTPCseed();
3235 FollowProlongation(*track, i2,1);
3239 track->SetBConstrain(1);
3240 // track->fLastPoint = i1+fInnerSec->GetNRows(); // first cluster in track position
3241 track->SetLastPoint(i1); // first cluster in track position
3242 track->SetFirstPoint(track->GetLastPoint());
3244 if (track->GetNumberOfClusters()<(i1-i2)*0.5 ||
3245 track->GetNumberOfClusters() < track->GetNFoundable()*0.6 ||
3246 track->GetNShared()>0.4*track->GetNumberOfClusters() ) {
3248 seed->~AliTPCseed();
3252 // Z VERTEX CONDITION
3253 Double_t zv, bz=GetBz();
3254 if ( !track->GetZAt(0.,bz,zv) ) continue;
3255 if (TMath::Abs(zv-z3)>cuts[2]) {
3256 FollowProlongation(*track, TMath::Max(i2-20,0));
3257 if ( !track->GetZAt(0.,bz,zv) ) continue;
3258 if (TMath::Abs(zv-z3)>cuts[2]){
3259 FollowProlongation(*track, TMath::Max(i2-40,0));
3260 if ( !track->GetZAt(0.,bz,zv) ) continue;
3261 if (TMath::Abs(zv-z3)>cuts[2] &&(track->GetNumberOfClusters() > track->GetNFoundable()*0.7)){
3262 // make seed without constrain
3263 AliTPCseed * track2 = MakeSeed(track,0.2,0.5,1.);
3264 FollowProlongation(*track2, i2,1);
3265 track2->SetBConstrain(kFALSE);
3266 track2->SetSeedType(1);
3267 arr->AddLast(track2);
3269 seed->~AliTPCseed();
3274 seed->~AliTPCseed();
3281 track->SetSeedType(0);
3282 arr->AddLast(track);
3283 seed = new AliTPCseed;
3285 // don't consider other combinations
3286 if (track->GetNumberOfClusters() > track->GetNFoundable()*0.8)
3292 Info("MakeSeeds3","\nSeeding statistic:\t%d\t%d\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin,nout1,nout2);
3298 void AliTPCtrackerMI::MakeSeeds5(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, Float_t cuts[4],
3303 //-----------------------------------------------------------------
3304 // This function creates track seeds.
3305 //-----------------------------------------------------------------
3306 // cuts[0] - fP4 cut
3307 // cuts[1] - tan(phi) cut
3308 // cuts[2] - zvertex cut
3309 // cuts[3] - fP3 cut
3319 Double_t x[5], c[15];
3321 // make temporary seed
3322 AliTPCseed * seed = new AliTPCseed;
3323 Double_t alpha=fOuterSec->GetAlpha(), shift=fOuterSec->GetAlphaShift();
3324 // Double_t cs=cos(alpha), sn=sin(alpha);
3329 Double_t x1 = GetXrow(i1-1);
3330 const AliTPCtrackerRow& kr1=GetRow(sec,i1-1);
3331 Double_t y1max = GetMaxY(i1-1)-kr1.GetDeadZone()-1.5;
3333 Double_t x1p = GetXrow(i1);
3334 const AliTPCtrackerRow& kr1p=GetRow(sec,i1);
3336 Double_t x1m = GetXrow(i1-2);
3337 const AliTPCtrackerRow& kr1m=GetRow(sec,i1-2);
3340 //last 3 padrow for seeding
3341 AliTPCtrackerRow& kr3 = GetRow((sec+fkNOS)%fkNOS,i1-7);
3342 Double_t x3 = GetXrow(i1-7);
3343 // Double_t y3max= GetMaxY(i1-7)-kr3.fDeadZone-1.5;
3345 AliTPCtrackerRow& kr3p = GetRow((sec+fkNOS)%fkNOS,i1-6);
3346 Double_t x3p = GetXrow(i1-6);
3348 AliTPCtrackerRow& kr3m = GetRow((sec+fkNOS)%fkNOS,i1-8);
3349 Double_t x3m = GetXrow(i1-8);
3354 Int_t im = i1-4; //middle pad row index
3355 Double_t xm = GetXrow(im); // radius of middle pad-row
3356 const AliTPCtrackerRow& krm=GetRow(sec,im); //middle pad -row
3357 // Double_t ymmax = GetMaxY(im)-kr1.fDeadZone-1.5;
3360 Double_t deltax = x1-x3;
3361 Double_t dymax = deltax*cuts[1];
3362 Double_t dzmax = deltax*cuts[3];
3364 // loop over clusters
3365 for (Int_t is=0; is < kr1; is++) {
3367 if (kr1[is]->IsUsed(10)) continue;
3368 Double_t y1=kr1[is]->GetY(), z1=kr1[is]->GetZ();
3370 if (deltay>0 && TMath::Abs(y1max-TMath::Abs(y1))> deltay ) continue; // seed only at the edge
3372 Int_t index1 = TMath::Max(kr3.Find(z1-dzmax)-1,0);
3373 Int_t index2 = TMath::Min(kr3.Find(z1+dzmax)+1,kr3);
3379 for (Int_t js=index1; js < index2; js++) {
3380 const AliTPCclusterMI *kcl = kr3[js];
3381 if (kcl->IsUsed(10)) continue;
3383 // apply angular cuts
3384 if (TMath::Abs(y1-y3)>dymax) continue;
3387 if (TMath::Abs(z1-z3)>dzmax) continue;
3389 Double_t angley = (y1-y3)/(x1-x3);
3390 Double_t anglez = (z1-z3)/(x1-x3);
3392 Double_t erry = TMath::Abs(angley)*(x1-x1m)*0.5+0.5;
3393 Double_t errz = TMath::Abs(anglez)*(x1-x1m)*0.5+0.5;
3395 Double_t yyym = angley*(xm-x1)+y1;
3396 Double_t zzzm = anglez*(xm-x1)+z1;
3398 const AliTPCclusterMI *kcm = krm.FindNearest2(yyym,zzzm,erry,errz,index);
3400 if (kcm->IsUsed(10)) continue;
3402 erry = TMath::Abs(angley)*(x1-x1m)*0.4+0.5;
3403 errz = TMath::Abs(anglez)*(x1-x1m)*0.4+0.5;
3410 // look around first
3411 const AliTPCclusterMI *kc1m = kr1m.FindNearest2(angley*(x1m-x1)+y1,
3417 if (kc1m->IsUsed(10)) used++;
3419 const AliTPCclusterMI *kc1p = kr1p.FindNearest2(angley*(x1p-x1)+y1,
3425 if (kc1p->IsUsed(10)) used++;
3427 if (used>1) continue;
3428 if (found<1) continue;
3432 const AliTPCclusterMI *kc3m = kr3m.FindNearest2(angley*(x3m-x3)+y3,
3438 if (kc3m->IsUsed(10)) used++;
3442 const AliTPCclusterMI *kc3p = kr3p.FindNearest2(angley*(x3p-x3)+y3,
3448 if (kc3p->IsUsed(10)) used++;
3452 if (used>1) continue;
3453 if (found<3) continue;
3463 x[4]=F1(x1,y1,x2,y2,x3,y3);
3464 //if (TMath::Abs(x[4]) >= cuts[0]) continue;
3467 x[2]=F2(x1,y1,x2,y2,x3,y3);
3470 x[3]=F3n(x1,y1,x2,y2,z1,z2,x[4]);
3471 //if (TMath::Abs(x[3]) > cuts[3]) continue;
3475 Double_t sy1=0.1, sz1=0.1;
3476 Double_t sy2=0.1, sz2=0.1;
3477 Double_t sy3=0.1, sy=0.1, sz=0.1;
3479 Double_t f40=(F1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
3480 Double_t f42=(F1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
3481 Double_t f43=(F1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
3482 Double_t f20=(F2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
3483 Double_t f22=(F2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
3484 Double_t f23=(F2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
3486 Double_t f30=(F3(x1,y1+sy,x2,y2,z1,z2)-x[3])/sy;
3487 Double_t f31=(F3(x1,y1,x2,y2,z1+sz,z2)-x[3])/sz;
3488 Double_t f32=(F3(x1,y1,x2,y2+sy,z1,z2)-x[3])/sy;
3489 Double_t f34=(F3(x1,y1,x2,y2,z1,z2+sz)-x[3])/sz;
3493 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3494 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3495 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3496 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3497 c[13]=f30*sy1*f40+f32*sy2*f42;
3498 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3500 // if (!BuildSeed(kr1[is],kcl,0,x1,x2,x3,x,c)) continue;
3502 index=kr1.GetIndex(is);
3503 seed->~AliTPCseed();
3504 AliTPCseed *track=new(seed) AliTPCseed(x1, sec*alpha+shift, x, c, index);
3506 track->SetIsSeeding(kTRUE);
3509 FollowProlongation(*track, i1-7,1);
3510 if (track->GetNumberOfClusters() < track->GetNFoundable()*0.75 ||
3511 track->GetNShared()>0.6*track->GetNumberOfClusters() || ( track->GetSigmaY2()+ track->GetSigmaZ2())>0.6){
3513 seed->~AliTPCseed();
3519 FollowProlongation(*track, i2,1);
3520 track->SetBConstrain(0);
3521 track->SetLastPoint(i1+fInnerSec->GetNRows()); // first cluster in track position
3522 track->SetFirstPoint(track->GetLastPoint());
3524 if (track->GetNumberOfClusters()<(i1-i2)*0.5 ||
3525 track->GetNumberOfClusters()<track->GetNFoundable()*0.7 ||
3526 track->GetNShared()>2. || track->GetChi2()/track->GetNumberOfClusters()>6 || ( track->GetSigmaY2()+ track->GetSigmaZ2())>0.5 ) {
3528 seed->~AliTPCseed();
3533 FollowProlongation(*track, TMath::Max(i2-10,0),1);
3534 AliTPCseed * track2 = MakeSeed(track,0.2,0.5,0.9);
3535 FollowProlongation(*track2, i2,1);
3536 track2->SetBConstrain(kFALSE);
3537 track2->SetSeedType(4);
3538 arr->AddLast(track2);
3540 seed->~AliTPCseed();
3544 //arr->AddLast(track);
3545 //seed = new AliTPCseed;
3551 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);
3557 //_____________________________________________________________________________
3558 void AliTPCtrackerMI::MakeSeeds2(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, Float_t */*cuts[4]*/,
3559 Float_t deltay, Bool_t /*bconstrain*/) {
3560 //-----------------------------------------------------------------
3561 // This function creates track seeds - without vertex constraint
3562 //-----------------------------------------------------------------
3563 // cuts[0] - fP4 cut - not applied
3564 // cuts[1] - tan(phi) cut
3565 // cuts[2] - zvertex cut - not applied
3566 // cuts[3] - fP3 cut
3576 Double_t alpha=fOuterSec->GetAlpha(), shift=fOuterSec->GetAlphaShift();
3577 // Double_t cs=cos(alpha), sn=sin(alpha);
3578 Int_t row0 = (i1+i2)/2;
3579 Int_t drow = (i1-i2)/2;
3580 const AliTPCtrackerRow& kr0=fSectors[sec][row0];
3581 AliTPCtrackerRow * kr=0;
3583 AliTPCpolyTrack polytrack;
3584 Int_t nclusters=fSectors[sec][row0];
3585 AliTPCseed * seed = new AliTPCseed;
3590 for (Int_t is=0; is < nclusters; is++) { //LOOP over clusters
3592 Int_t nfoundable =0;
3593 for (Int_t iter =1; iter<2; iter++){ //iterations
3594 const AliTPCtrackerRow& krm=fSectors[sec][row0-iter];
3595 const AliTPCtrackerRow& krp=fSectors[sec][row0+iter];
3596 const AliTPCclusterMI * cl= kr0[is];
3598 if (cl->IsUsed(10)) {
3604 Double_t x = kr0.GetX();
3605 // Initialization of the polytrack
3610 Double_t y0= cl->GetY();
3611 Double_t z0= cl->GetZ();
3615 Double_t ymax = fSectors->GetMaxY(row0)-kr0.GetDeadZone()-1.5;
3616 if (deltay>0 && TMath::Abs(ymax-TMath::Abs(y0))> deltay ) continue; // seed only at the edge
3618 erry = (0.5)*cl->GetSigmaY2()/TMath::Sqrt(cl->GetQ())*6;
3619 errz = (0.5)*cl->GetSigmaZ2()/TMath::Sqrt(cl->GetQ())*6;
3620 polytrack.AddPoint(x,y0,z0,erry, errz);
3623 if (cl->IsUsed(10)) sumused++;
3626 Float_t roady = (5*TMath::Sqrt(cl->GetSigmaY2()+0.2)+1.)*iter;
3627 Float_t roadz = (5*TMath::Sqrt(cl->GetSigmaZ2()+0.2)+1.)*iter;
3630 AliTPCclusterMI * cl1 = krm.FindNearest(y0,z0,roady,roadz);
3631 if (cl1 && TMath::Abs(ymax-TMath::Abs(y0))) {
3632 erry = (0.5)*cl1->GetSigmaY2()/TMath::Sqrt(cl1->GetQ())*3;
3633 errz = (0.5)*cl1->GetSigmaZ2()/TMath::Sqrt(cl1->GetQ())*3;
3634 if (cl1->IsUsed(10)) sumused++;
3635 polytrack.AddPoint(x,cl1->GetY(),cl1->GetZ(),erry,errz);
3639 AliTPCclusterMI * cl2 = krp.FindNearest(y0,z0,roady,roadz);
3641 erry = (0.5)*cl2->GetSigmaY2()/TMath::Sqrt(cl2->GetQ())*3;
3642 errz = (0.5)*cl2->GetSigmaZ2()/TMath::Sqrt(cl2->GetQ())*3;
3643 if (cl2->IsUsed(10)) sumused++;
3644 polytrack.AddPoint(x,cl2->GetY(),cl2->GetZ(),erry,errz);
3647 if (sumused>0) continue;
3649 polytrack.UpdateParameters();
3655 nfoundable = polytrack.GetN();
3656 nfound = nfoundable;
3658 for (Int_t ddrow = iter+1; ddrow<drow;ddrow++){
3659 Float_t maxdist = 0.8*(1.+3./(ddrow));
3660 for (Int_t delta = -1;delta<=1;delta+=2){
3661 Int_t row = row0+ddrow*delta;
3662 kr = &(fSectors[sec][row]);
3663 Double_t xn = kr->GetX();
3664 Double_t ymax1 = fSectors->GetMaxY(row)-kr->GetDeadZone()-1.5;
3665 polytrack.GetFitPoint(xn,yn,zn);
3666 if (TMath::Abs(yn)>ymax1) continue;
3668 AliTPCclusterMI * cln = kr->FindNearest(yn,zn,roady,roadz);
3670 Float_t dist = TMath::Sqrt( (yn-cln->GetY())*(yn-cln->GetY())+(zn-cln->GetZ())*(zn-cln->GetZ()));
3673 erry = (dist+0.3)*cln->GetSigmaY2()/TMath::Sqrt(cln->GetQ())*(1.+1./(ddrow));
3674 errz = (dist+0.3)*cln->GetSigmaZ2()/TMath::Sqrt(cln->GetQ())*(1.+1./(ddrow));
3675 if (cln->IsUsed(10)) {
3676 // printf("used\n");
3684 polytrack.AddPoint(xn,cln->GetY(),cln->GetZ(),erry, errz);
3689 if ( (sumused>3) || (sumused>0.5*nfound) || (nfound<0.6*nfoundable)) break;
3690 polytrack.UpdateParameters();
3693 if ( (sumused>3) || (sumused>0.5*nfound)) {
3694 //printf("sumused %d\n",sumused);
3699 polytrack.GetFitDerivation(kr0.GetX(),dy,dz);
3700 AliTPCpolyTrack track2;
3702 polytrack.Refit(track2,0.5+TMath::Abs(dy)*0.3,0.4+TMath::Abs(dz)*0.3);
3703 if (track2.GetN()<0.5*nfoundable) continue;
3706 if ((nfound>0.6*nfoundable) &&( nfoundable>0.4*(i1-i2))) {
3708 // test seed with and without constrain
3709 for (Int_t constrain=0; constrain<=0;constrain++){
3710 // add polytrack candidate
3712 Double_t x[5], c[15];
3713 Double_t x1,x2,x3,y1,y2,y3,z1,z2,z3;
3714 track2.GetBoundaries(x3,x1);
3716 track2.GetFitPoint(x1,y1,z1);
3717 track2.GetFitPoint(x2,y2,z2);
3718 track2.GetFitPoint(x3,y3,z3);
3720 //is track pointing to the vertex ?
3723 polytrack.GetFitPoint(x0,y0,z0);
3736 x[4]=F1(x1,y1,x2,y2,x3,y3);
3738 // if (TMath::Abs(x[4]) >= cuts[0]) continue; //
3739 x[2]=F2(x1,y1,x2,y2,x3,y3);
3741 //if (TMath::Abs(x[4]*x1-x[2]) >= cuts[1]) continue;
3742 //x[3]=F3(x1,y1,x2,y2,z1,z2);
3743 x[3]=F3n(x1,y1,x3,y3,z1,z3,x[4]);
3744 //if (TMath::Abs(x[3]) > cuts[3]) continue;
3747 Double_t sy =0.1, sz =0.1;
3748 Double_t sy1=0.02, sz1=0.02;
3749 Double_t sy2=0.02, sz2=0.02;
3753 sy3=25000*x[4]*x[4]+0.1, sy=0.1, sz=0.1;
3756 Double_t f40=(F1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
3757 Double_t f42=(F1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
3758 Double_t f43=(F1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
3759 Double_t f20=(F2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
3760 Double_t f22=(F2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
3761 Double_t f23=(F2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
3763 Double_t f30=(F3(x1,y1+sy,x3,y3,z1,z3)-x[3])/sy;
3764 Double_t f31=(F3(x1,y1,x3,y3,z1+sz,z3)-x[3])/sz;
3765 Double_t f32=(F3(x1,y1,x3,y3+sy,z1,z3)-x[3])/sy;
3766 Double_t f34=(F3(x1,y1,x3,y3,z1,z3+sz)-x[3])/sz;
3771 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3772 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3773 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3774 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3775 c[13]=f30*sy1*f40+f32*sy2*f42;
3776 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3778 //Int_t row1 = fSectors->GetRowNumber(x1);
3779 Int_t row1 = GetRowNumber(x1);
3783 seed->~AliTPCseed();
3784 AliTPCseed *track=new(seed) AliTPCseed(x1,sec*alpha+shift,x,c,index);
3785 track->SetIsSeeding(kTRUE);
3786 Int_t rc=FollowProlongation(*track, i2);
3787 if (constrain) track->SetBConstrain(1);
3789 track->SetBConstrain(0);
3790 track->SetLastPoint(row1+fInnerSec->GetNRows()); // first cluster in track position
3791 track->SetFirstPoint(track->GetLastPoint());
3793 if (rc==0 || track->GetNumberOfClusters()<(i1-i2)*0.5 ||
3794 track->GetNumberOfClusters() < track->GetNFoundable()*0.6 ||
3795 track->GetNShared()>0.4*track->GetNumberOfClusters()) {
3798 seed->~AliTPCseed();
3801 arr->AddLast(track);
3802 seed = new AliTPCseed;
3806 } // if accepted seed
3809 Info("MakeSeeds2","\nSeeding statiistic:\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin3);
3815 AliTPCseed *AliTPCtrackerMI::MakeSeed(AliTPCseed *const track, Float_t r0, Float_t r1, Float_t r2)
3819 //reseed using track points
3820 Int_t p0 = int(r0*track->GetNumberOfClusters()); // point 0
3821 Int_t p1 = int(r1*track->GetNumberOfClusters());
3822 Int_t p2 = int(r2*track->GetNumberOfClusters()); // last point
3824 Double_t x0[3],x1[3],x2[3];
3825 for (Int_t i=0;i<3;i++){
3831 // find track position at given ratio of the length
3832 Int_t sec0=0, sec1=0, sec2=0;
3835 for (Int_t i=0;i<160;i++){
3836 if (track->GetClusterPointer(i)){
3838 AliTPCTrackerPoint *trpoint =track->GetTrackPoint(i);
3839 if ( (index<p0) || x0[0]<0 ){
3840 if (trpoint->GetX()>1){
3841 clindex = track->GetClusterIndex2(i);
3843 x0[0] = trpoint->GetX();
3844 x0[1] = trpoint->GetY();
3845 x0[2] = trpoint->GetZ();
3846 sec0 = ((clindex&0xff000000)>>24)%18;
3851 if ( (index<p1) &&(trpoint->GetX()>1)){
3852 clindex = track->GetClusterIndex2(i);
3854 x1[0] = trpoint->GetX();
3855 x1[1] = trpoint->GetY();
3856 x1[2] = trpoint->GetZ();
3857 sec1 = ((clindex&0xff000000)>>24)%18;
3860 if ( (index<p2) &&(trpoint->GetX()>1)){
3861 clindex = track->GetClusterIndex2(i);
3863 x2[0] = trpoint->GetX();
3864 x2[1] = trpoint->GetY();
3865 x2[2] = trpoint->GetZ();
3866 sec2 = ((clindex&0xff000000)>>24)%18;
3873 Double_t alpha, cs,sn, xx2,yy2;
3875 alpha = (sec1-sec2)*fSectors->GetAlpha();
3876 cs = TMath::Cos(alpha);
3877 sn = TMath::Sin(alpha);
3878 xx2= x1[0]*cs-x1[1]*sn;
3879 yy2= x1[0]*sn+x1[1]*cs;
3883 alpha = (sec0-sec2)*fSectors->GetAlpha();
3884 cs = TMath::Cos(alpha);
3885 sn = TMath::Sin(alpha);
3886 xx2= x0[0]*cs-x0[1]*sn;
3887 yy2= x0[0]*sn+x0[1]*cs;
3893 Double_t x[5],c[15];
3897 x[4]=F1(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]);
3898 // if (x[4]>1) return 0;
3899 x[2]=F2(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]);
3900 x[3]=F3n(x2[0],x2[1],x0[0],x0[1],x2[2],x0[2],x[4]);
3901 //if (TMath::Abs(x[3]) > 2.2) return 0;
3902 //if (TMath::Abs(x[2]) > 1.99) return 0;
3904 Double_t sy =0.1, sz =0.1;
3906 Double_t sy1=0.02+track->GetSigmaY2(), sz1=0.02+track->GetSigmaZ2();
3907 Double_t sy2=0.01+track->GetSigmaY2(), sz2=0.01+track->GetSigmaZ2();
3908 Double_t sy3=0.01+track->GetSigmaY2();
3910 Double_t f40=(F1(x2[0],x2[1]+sy,x1[0],x1[1],x0[0],x0[1])-x[4])/sy;
3911 Double_t f42=(F1(x2[0],x2[1],x1[0],x1[1]+sy,x0[0],x0[1])-x[4])/sy;
3912 Double_t f43=(F1(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]+sy)-x[4])/sy;
3913 Double_t f20=(F2(x2[0],x2[1]+sy,x1[0],x1[1],x0[0],x0[1])-x[2])/sy;
3914 Double_t f22=(F2(x2[0],x2[1],x1[0],x1[1]+sy,x0[0],x0[1])-x[2])/sy;
3915 Double_t f23=(F2(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]+sy)-x[2])/sy;
3917 Double_t f30=(F3(x2[0],x2[1]+sy,x0[0],x0[1],x2[2],x0[2])-x[3])/sy;
3918 Double_t f31=(F3(x2[0],x2[1],x0[0],x0[1],x2[2]+sz,x0[2])-x[3])/sz;
3919 Double_t f32=(F3(x2[0],x2[1],x0[0],x0[1]+sy,x2[2],x0[2])-x[3])/sy;
3920 Double_t f34=(F3(x2[0],x2[1],x0[0],x0[1],x2[2],x0[2]+sz)-x[3])/sz;
3925 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3926 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3927 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3928 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3929 c[13]=f30*sy1*f40+f32*sy2*f42;
3930 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3932 // Int_t row1 = fSectors->GetRowNumber(x2[0]);
3933 AliTPCseed *seed=new AliTPCseed(x2[0], sec2*fSectors->GetAlpha()+fSectors->GetAlphaShift(), x, c, 0);
3934 // Double_t y0,z0,y1,z1, y2,z2;
3935 //seed->GetProlongation(x0[0],y0,z0);
3936 // seed->GetProlongation(x1[0],y1,z1);
3937 //seed->GetProlongation(x2[0],y2,z2);
3939 seed->SetLastPoint(pp2);
3940 seed->SetFirstPoint(pp2);
3947 AliTPCseed *AliTPCtrackerMI::ReSeed(const AliTPCseed *track, Float_t r0, Float_t r1, Float_t r2)
3951 //reseed using founded clusters
3953 // Find the number of clusters
3954 Int_t nclusters = 0;
3955 for (Int_t irow=0;irow<160;irow++){
3956 if (track->GetClusterIndex(irow)>0) nclusters++;
3960 ipos[0] = TMath::Max(int(r0*nclusters),0); // point 0 cluster
3961 ipos[1] = TMath::Min(int(r1*nclusters),nclusters-1); //
3962 ipos[2] = TMath::Min(int(r2*nclusters),nclusters-1); // last point
3965 Double_t xyz[3][3]={{0}};
3966 Int_t row[3]={0},sec[3]={0,0,0};
3968 // find track row position at given ratio of the length
3970 for (Int_t irow=0;irow<160;irow++){
3971 if (track->GetClusterIndex2(irow)<0) continue;
3973 for (Int_t ipoint=0;ipoint<3;ipoint++){
3974 if (index<=ipos[ipoint]) row[ipoint] = irow;
3978 //Get cluster and sector position
3979 for (Int_t ipoint=0;ipoint<3;ipoint++){
3980 Int_t clindex = track->GetClusterIndex2(row[ipoint]);
3981 AliTPCclusterMI * cl = GetClusterMI(clindex);
3984 // AliTPCclusterMI * cl = GetClusterMI(clindex);
3987 sec[ipoint] = ((clindex&0xff000000)>>24)%18;
3988 xyz[ipoint][0] = GetXrow(row[ipoint]);
3989 xyz[ipoint][1] = cl->GetY();
3990 xyz[ipoint][2] = cl->GetZ();
3994 // Calculate seed state vector and covariance matrix
3996 Double_t alpha, cs,sn, xx2,yy2;
3998 alpha = (sec[1]-sec[2])*fSectors->GetAlpha();
3999 cs = TMath::Cos(alpha);
4000 sn = TMath::Sin(alpha);
4001 xx2= xyz[1][0]*cs-xyz[1][1]*sn;
4002 yy2= xyz[1][0]*sn+xyz[1][1]*cs;
4006 alpha = (sec[0]-sec[2])*fSectors->GetAlpha();
4007 cs = TMath::Cos(alpha);
4008 sn = TMath::Sin(alpha);
4009 xx2= xyz[0][0]*cs-xyz[0][1]*sn;
4010 yy2= xyz[0][0]*sn+xyz[0][1]*cs;
4016 Double_t x[5],c[15];
4020 x[4]=F1(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
4021 x[2]=F2(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
4022 x[3]=F3n(xyz[2][0],xyz[2][1],xyz[0][0],xyz[0][1],xyz[2][2],xyz[0][2],x[4]);
4024 Double_t sy =0.1, sz =0.1;
4026 Double_t sy1=0.2, sz1=0.2;
4027 Double_t sy2=0.2, sz2=0.2;
4030 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;
4031 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;
4032 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;
4033 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;
4034 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;
4035 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;
4037 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;
4038 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;
4039 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;
4040 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;
4045 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
4046 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
4047 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
4048 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
4049 c[13]=f30*sy1*f40+f32*sy2*f42;
4050 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
4052 // Int_t row1 = fSectors->GetRowNumber(xyz[2][0]);
4053 AliTPCseed *seed=new AliTPCseed(xyz[2][0], sec[2]*fSectors->GetAlpha()+fSectors->GetAlphaShift(), x, c, 0);
4054 seed->SetLastPoint(row[2]);
4055 seed->SetFirstPoint(row[2]);
4060 AliTPCseed *AliTPCtrackerMI::ReSeed(AliTPCseed *track,Int_t r0, Bool_t forward)
4064 //reseed using founded clusters
4067 Int_t row[3]={0,0,0};
4068 Int_t sec[3]={0,0,0};
4070 // forward direction
4072 for (Int_t irow=r0;irow<160;irow++){
4073 if (track->GetClusterIndex(irow)>0){
4078 for (Int_t irow=160;irow>r0;irow--){
4079 if (track->GetClusterIndex(irow)>0){
4084 for (Int_t irow=row[2]-15;irow>row[0];irow--){
4085 if (track->GetClusterIndex(irow)>0){
4093 for (Int_t irow=0;irow<r0;irow++){
4094 if (track->GetClusterIndex(irow)>0){
4099 for (Int_t irow=r0;irow>0;irow--){
4100 if (track->GetClusterIndex(irow)>0){
4105 for (Int_t irow=row[2]-15;irow>row[0];irow--){
4106 if (track->GetClusterIndex(irow)>0){
4113 if ((row[2]-row[0])<20) return 0;
4114 if (row[1]==0) return 0;
4117 //Get cluster and sector position
4118 for (Int_t ipoint=0;ipoint<3;ipoint++){
4119 Int_t clindex = track->GetClusterIndex2(row[ipoint]);
4120 AliTPCclusterMI * cl = GetClusterMI(clindex);
4123 // AliTPCclusterMI * cl = GetClusterMI(clindex);
4126 sec[ipoint] = ((clindex&0xff000000)>>24)%18;
4127 xyz[ipoint][0] = GetXrow(row[ipoint]);
4128 AliTPCTrackerPoint * point = track->GetTrackPoint(row[ipoint]);
4129 if (point&&ipoint<2){
4131 xyz[ipoint][1] = point->GetY();
4132 xyz[ipoint][2] = point->GetZ();
4135 xyz[ipoint][1] = cl->GetY();
4136 xyz[ipoint][2] = cl->GetZ();
4143 // Calculate seed state vector and covariance matrix
4145 Double_t alpha, cs,sn, xx2,yy2;
4147 alpha = (sec[1]-sec[2])*fSectors->GetAlpha();
4148 cs = TMath::Cos(alpha);
4149 sn = TMath::Sin(alpha);
4150 xx2= xyz[1][0]*cs-xyz[1][1]*sn;
4151 yy2= xyz[1][0]*sn+xyz[1][1]*cs;
4155 alpha = (sec[0]-sec[2])*fSectors->GetAlpha();
4156 cs = TMath::Cos(alpha);
4157 sn = TMath::Sin(alpha);
4158 xx2= xyz[0][0]*cs-xyz[0][1]*sn;
4159 yy2= xyz[0][0]*sn+xyz[0][1]*cs;
4165 Double_t x[5],c[15];
4169 x[4]=F1(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
4170 x[2]=F2(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
4171 x[3]=F3n(xyz[2][0],xyz[2][1],xyz[0][0],xyz[0][1],xyz[2][2],xyz[0][2],x[4]);
4173 Double_t sy =0.1, sz =0.1;
4175 Double_t sy1=0.2, sz1=0.2;
4176 Double_t sy2=0.2, sz2=0.2;
4179 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;
4180 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;
4181 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;
4182 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;
4183 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;
4184 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;
4186 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;
4187 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;
4188 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;
4189 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;
4194 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
4195 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
4196 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
4197 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
4198 c[13]=f30*sy1*f40+f32*sy2*f42;
4199 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
4201 // Int_t row1 = fSectors->GetRowNumber(xyz[2][0]);
4202 AliTPCseed *seed=new AliTPCseed(xyz[2][0], sec[2]*fSectors->GetAlpha()+fSectors->GetAlphaShift(), x, c, 0);
4203 seed->SetLastPoint(row[2]);
4204 seed->SetFirstPoint(row[2]);
4205 for (Int_t i=row[0];i<row[2];i++){
4206 seed->SetClusterIndex(i, track->GetClusterIndex(i));
4214 void AliTPCtrackerMI::FindMultiMC(const TObjArray * array, AliESDEvent */*esd*/, Int_t iter)
4217 // find multi tracks - THIS FUNCTION IS ONLY FOR DEBUG PURPOSES
4219 // Use AliTPCReconstructor::StreamLevel()>2 if you want to tune parameters - cuts
4221 // Two reasons to have multiple find tracks
4222 // 1. Curling tracks can be find more than once
4223 // 2. Splitted tracks
4224 // a.) Multiple seeding to increase tracking efficiency - (~ 100% reached)
4225 // b.) Edge effect on the sector boundaries
4228 // Algorithm done in 2 phases - because of CPU consumption
4229 // it is n^2 algorithm - for lead-lead 20000x20000 combination are investigated
4231 // Algorihm for curling tracks sign:
4232 // 1 phase -makes a very rough fast cuts to minimize combinatorics
4233 // a.) opposite sign
4234 // b.) one of the tracks - not pointing to the primary vertex -
4235 // c.) delta tan(theta)
4237 // 2 phase - calculates DCA between tracks - time consument
4242 // General cuts - for splitted tracks and for curling tracks
4244 const Float_t kMaxdPhi = 0.2; // maximal distance in phi
4246 // Curling tracks cuts
4251 Int_t nentries = array->GetEntriesFast();
4252 AliHelix *helixes = new AliHelix[nentries];
4253 Float_t *xm = new Float_t[nentries];
4254 Float_t *dz0 = new Float_t[nentries];
4255 Float_t *dz1 = new Float_t[nentries];
4261 // Find track COG in x direction - point with best defined parameters
4263 for (Int_t i=0;i<nentries;i++){
4264 AliTPCseed* track = (AliTPCseed*)array->At(i);
4265 if (!track) continue;
4266 track->SetCircular(0);
4267 new (&helixes[i]) AliHelix(*track);
4271 track->GetDZ(GetX(),GetY(),GetZ(),GetBz(),dz);
4274 for (Int_t icl=0; icl<160; icl++){
4275 AliTPCclusterMI * cl = track->GetClusterPointer(icl);
4281 if (ncl>0) xm[i]/=Float_t(ncl);
4284 for (Int_t i0=0;i0<nentries;i0++){
4285 AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
4286 if (!track0) continue;
4287 Float_t xc0 = helixes[i0].GetHelix(6);
4288 Float_t yc0 = helixes[i0].GetHelix(7);
4289 Float_t r0 = helixes[i0].GetHelix(8);
4290 Float_t rc0 = TMath::Sqrt(xc0*xc0+yc0*yc0);
4291 Float_t fi0 = TMath::ATan2(yc0,xc0);
4293 for (Int_t i1=i0+1;i1<nentries;i1++){
4294 AliTPCseed * track1 = (AliTPCseed*)array->At(i1);
4295 if (!track1) continue;
4296 Int_t lab0=track0->GetLabel();
4297 Int_t lab1=track1->GetLabel();
4298 if (TMath::Abs(lab0)!=TMath::Abs(lab1)) continue;
4300 Float_t xc1 = helixes[i1].GetHelix(6);
4301 Float_t yc1 = helixes[i1].GetHelix(7);
4302 Float_t r1 = helixes[i1].GetHelix(8);
4303 Float_t rc1 = TMath::Sqrt(xc1*xc1+yc1*yc1);
4304 Float_t fi1 = TMath::ATan2(yc1,xc1);
4306 Float_t dfi = fi0-fi1;
4309 if (dfi>1.5*TMath::Pi()) dfi-=TMath::Pi(); // take care about edge effect
4310 if (dfi<-1.5*TMath::Pi()) dfi+=TMath::Pi(); //
4311 if (TMath::Abs(dfi)>kMaxdPhi&&helixes[i0].GetHelix(4)*helixes[i1].GetHelix(4)<0){
4313 // if short tracks with undefined sign
4314 fi1 = -TMath::ATan2(yc1,-xc1);
4317 Float_t dtheta = TMath::Abs(track0->GetTgl()-track1->GetTgl())<TMath::Abs(track0->GetTgl()+track1->GetTgl())? track0->GetTgl()-track1->GetTgl():track0->GetTgl()+track1->GetTgl();
4320 // debug stream to tune "fast cuts"
4322 Double_t dist[3]; // distance at X
4323 Double_t mdist[3]={0,0,0}; // mean distance X+-40cm
4324 track0->GetDistance(track1,0.5*(xm[i0]+xm[i1])-40.,dist,AliTracker::GetBz());
4325 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[i]);
4326 track0->GetDistance(track1,0.5*(xm[i0]+xm[i1])+40.,dist,AliTracker::GetBz());
4327 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[i]);
4328 track0->GetDistance(track1,0.5*(xm[i0]+xm[i1]),dist,AliTracker::GetBz());
4329 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[i]);
4330 for (Int_t i=0;i<3;i++) mdist[i]*=0.33333;
4334 for (Int_t icl=0; icl<160; icl++){
4335 AliTPCclusterMI * cl0 = track0->GetClusterPointer(icl);
4336 AliTPCclusterMI * cl1 = track1->GetClusterPointer(icl);
4339 if (cl0==cl1) sums++;
4343 if (AliTPCReconstructor::StreamLevel()>5) {
4344 TTreeSRedirector &cstream = *fDebugStreamer;
4349 "Tr0.="<<track0<< // seed0
4350 "Tr1.="<<track1<< // seed1
4351 "h0.="<<&helixes[i0]<<
4352 "h1.="<<&helixes[i1]<<
4354 "sum="<<sum<< //the sum of rows with cl in both
4355 "sums="<<sums<< //the sum of shared clusters
4356 "xm0="<<xm[i0]<< // the center of track
4357 "xm1="<<xm[i1]<< // the x center of track
4358 // General cut variables
4359 "dfi="<<dfi<< // distance in fi angle
4360 "dtheta="<<dtheta<< // distance int theta angle
4366 "dist0="<<dist[0]<< //distance x
4367 "dist1="<<dist[1]<< //distance y
4368 "dist2="<<dist[2]<< //distance z
4369 "mdist0="<<mdist[0]<< //distance x
4370 "mdist1="<<mdist[1]<< //distance y
4371 "mdist2="<<mdist[2]<< //distance z
4387 if (AliTPCReconstructor::StreamLevel()>1) {
4388 AliInfo("Time for curling tracks removal DEBUGGING MC");
4395 void AliTPCtrackerMI::FindSplitted(TObjArray * array, AliESDEvent */*esd*/, Int_t /*iter*/){
4397 // Find Splitted tracks and remove the one with worst quality
4398 // Corresponding debug streamer to tune selections - "Splitted2"
4400 // 0. Sort tracks according quility
4401 // 1. Propagate the tracks to the reference radius
4402 // 2. Double_t loop to select close tracks (only to speed up process)
4403 // 3. Calculate cluster overlap ratio - and remove the track if bigger than a threshold
4404 // 4. Delete temporary parameters
4406 const Double_t xref=GetXrow(63); // reference radius -IROC/OROC boundary
4408 const Double_t kCutP1=10; // delta Z cut 10 cm
4409 const Double_t kCutP2=0.15; // delta snp(fi) cut 0.15
4410 const Double_t kCutP3=0.15; // delta tgl(theta) cut 0.15
4411 const Double_t kCutAlpha=0.15; // delta alpha cut
4412 Int_t firstpoint = 0;
4413 Int_t lastpoint = 160;
4415 Int_t nentries = array->GetEntriesFast();
4416 AliExternalTrackParam *params = new AliExternalTrackParam[nentries];
4422 //0. Sort tracks according quality
4423 //1. Propagate the ext. param to reference radius
4424 Int_t nseed = array->GetEntriesFast();
4425 if (nseed<=0) return;
4426 Float_t * quality = new Float_t[nseed];
4427 Int_t * indexes = new Int_t[nseed];
4428 for (Int_t i=0; i<nseed; i++) {
4429 AliTPCseed *pt=(AliTPCseed*)array->UncheckedAt(i);
4434 pt->UpdatePoints(); //select first last max dens points
4435 Float_t * points = pt->GetPoints();
4436 if (points[3]<0.8) quality[i] =-1;
4437 quality[i] = (points[2]-points[0])+pt->GetNumberOfClusters();
4438 //prefer high momenta tracks if overlaps
4439 quality[i] *= TMath::Sqrt(TMath::Abs(pt->Pt())+0.5);
4441 AliTracker::PropagateTrackToBxByBz(&(params[i]),xref,pt->GetMass(),5.,kTRUE);
4442 AliTracker::PropagateTrackToBxByBz(&(params[i]),xref,pt->GetMass(),1.,kTRUE);
4444 TMath::Sort(nseed,quality,indexes);
4446 // 3. Loop over pair of tracks
4448 for (Int_t i0=0; i0<nseed; i0++) {
4449 Int_t index0=indexes[i0];
4450 if (!(array->UncheckedAt(index0))) continue;
4451 AliTPCseed *s1 = (AliTPCseed*)array->UncheckedAt(index0);
4452 if (!s1->IsActive()) continue;
4453 AliExternalTrackParam &par0=params[index0];
4454 for (Int_t i1=i0+1; i1<nseed; i1++) {
4455 Int_t index1=indexes[i1];
4456 if (!(array->UncheckedAt(index1))) continue;
4457 AliTPCseed *s2 = (AliTPCseed*)array->UncheckedAt(index1);
4458 if (!s2->IsActive()) continue;
4459 if (s2->GetKinkIndexes()[0]!=0)
4460 if (s2->GetKinkIndexes()[0] == -s1->GetKinkIndexes()[0]) continue;
4461 AliExternalTrackParam &par1=params[index1];
4462 if (TMath::Abs(par0.GetParameter()[3]-par1.GetParameter()[3])>kCutP3) continue;
4463 if (TMath::Abs(par0.GetParameter()[1]-par1.GetParameter()[1])>kCutP1) continue;
4464 if (TMath::Abs(par0.GetParameter()[2]-par1.GetParameter()[2])>kCutP2) continue;
4465 Double_t dAlpha= TMath::Abs(par0.GetAlpha()-par1.GetAlpha());
4466 if (dAlpha>TMath::Pi()) dAlpha-=TMath::Pi();
4467 if (TMath::Abs(dAlpha)>kCutAlpha) continue;
4472 Int_t firstShared=lastpoint, lastShared=firstpoint;
4473 Int_t firstRow=lastpoint, lastRow=firstpoint;
4475 for (Int_t i=firstpoint;i<lastpoint;i++){
4476 if (s1->GetClusterIndex2(i)>0) nall0++;
4477 if (s2->GetClusterIndex2(i)>0) nall1++;
4478 if (s1->GetClusterIndex2(i)>0 && s2->GetClusterIndex2(i)>0) {
4479 if (i<firstRow) firstRow=i;
4480 if (i>lastRow) lastRow=i;
4482 if ( (s1->GetClusterIndex2(i))==(s2->GetClusterIndex2(i)) && s1->GetClusterIndex2(i)>0) {
4483 if (i<firstShared) firstShared=i;
4484 if (i>lastShared) lastShared=i;
4488 Double_t ratio0 = Float_t(sumShared)/Float_t(TMath::Min(nall0+1,nall1+1));
4489 Double_t ratio1 = Float_t(sumShared)/Float_t(TMath::Max(nall0+1,nall1+1));
4491 if( AliTPCReconstructor::StreamLevel()>1){
4492 TTreeSRedirector &cstream = *fDebugStreamer;
4493 Int_t n0=s1->GetNumberOfClusters();
4494 Int_t n1=s2->GetNumberOfClusters();
4495 Int_t n0F=s1->GetNFoundable();
4496 Int_t n1F=s2->GetNFoundable();
4497 Int_t lab0=s1->GetLabel();
4498 Int_t lab1=s2->GetLabel();
4500 cstream<<"Splitted2"<<
4501 "iter="<<fIteration<<
4502 "lab0="<<lab0<< // MC label if exist
4503 "lab1="<<lab1<< // MC label if exist
4506 "ratio0="<<ratio0<< // shared ratio
4507 "ratio1="<<ratio1<< // shared ratio
4508 "p0.="<<&par0<< // track parameters
4510 "s0.="<<s1<< // full seed
4512 "n0="<<n0<< // number of clusters track 0
4513 "n1="<<n1<< // number of clusters track 1
4514 "nall0="<<nall0<< // number of clusters track 0
4515 "nall1="<<nall1<< // number of clusters track 1
4516 "n0F="<<n0F<< // number of findable
4517 "n1F="<<n1F<< // number of findable
4518 "shared="<<sumShared<< // number of shared clusters
4519 "firstS="<<firstShared<< // first and the last shared row
4520 "lastS="<<lastShared<<
4521 "firstRow="<<firstRow<< // first and the last row with cluster
4522 "lastRow="<<lastRow<< //
4526 // remove track with lower quality
4528 if (ratio0>AliTPCReconstructor::GetRecoParam()->GetCutSharedClusters(0) ||
4529 ratio1>AliTPCReconstructor::GetRecoParam()->GetCutSharedClusters(1)){
4533 delete array->RemoveAt(index1);
4538 // 4. Delete temporary array
4548 void AliTPCtrackerMI::FindCurling(const TObjArray * array, AliESDEvent */*esd*/, Int_t iter)
4551 // find Curling tracks
4552 // Use AliTPCReconstructor::StreamLevel()>1 if you want to tune parameters - cuts
4555 // Algorithm done in 2 phases - because of CPU consumption
4556 // it is n^2 algorithm - for lead-lead 20000x20000 combination are investigated
4557 // see detal in MC part what can be used to cut
4561 const Float_t kMaxC = 400; // maximal curvature to of the track
4562 const Float_t kMaxdTheta = 0.15; // maximal distance in theta
4563 const Float_t kMaxdPhi = 0.15; // maximal distance in phi
4564 const Float_t kPtRatio = 0.3; // ratio between pt
4565 const Float_t kMinDCAR = 2.; // distance to the primary vertex in r - see cpipe cut
4568 // Curling tracks cuts
4571 const Float_t kMaxDeltaRMax = 40; // distance in outer radius
4572 const Float_t kMaxDeltaRMin = 5.; // distance in lower radius - see cpipe cut
4573 const Float_t kMinAngle = 2.9; // angle between tracks
4574 const Float_t kMaxDist = 5; // biggest distance
4576 // The cuts can be tuned using the "MC information stored in Multi tree ==> see FindMultiMC
4579 TCut csign("csign","Tr0.fP[4]*Tr1.fP[4]<0"); //opposite sign
4580 TCut cmax("cmax","abs(Tr0.GetC())>1/400");
4581 TCut cda("cda","sqrt(dtheta^2+dfi^2)<0.15");
4582 TCut ccratio("ccratio","abs((Tr0.fP[4]+Tr1.fP[4])/(abs(Tr0.fP[4])+abs(Tr1.fP[4])))<0.3");
4583 TCut cpipe("cpipe", "min(abs(r0-rc0),abs(r1-rc1))>5");
4585 TCut cdrmax("cdrmax","abs(abs(rc0+r0)-abs(rc1+r1))<40")
4586 TCut cdrmin("cdrmin","abs(abs(rc0+r0)-abs(rc1+r1))<10")
4588 Multi->Draw("dfi","iter==0"+csign+cmax+cda+ccratio); ~94% of curling tracks fulfill
4589 Multi->Draw("min(abs(r0-rc0),abs(r1-rc1))","iter==0&&abs(lab1)==abs(lab0)"+csign+cmax+cda+ccratio+cpipe+cdrmin+cdrmax); //80%
4591 Curling2->Draw("dfi","iter==0&&abs(lab0)==abs(lab1)"+csign+cmax+cdtheta+cdfi+ccratio)
4597 Int_t nentries = array->GetEntriesFast();
4598 AliHelix *helixes = new AliHelix[nentries];
4599 for (Int_t i=0;i<nentries;i++){
4600 AliTPCseed* track = (AliTPCseed*)array->At(i);
4601 if (!track) continue;
4602 track->SetCircular(0);
4603 new (&helixes[i]) AliHelix(*track);
4609 Double_t phase[2][2]={{0,0},{0,0}},radius[2]={0,0};
4615 for (Int_t i0=0;i0<nentries;i0++){
4616 AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
4617 if (!track0) continue;
4618 if (TMath::Abs(track0->GetC())<1/kMaxC) continue;
4619 Float_t xc0 = helixes[i0].GetHelix(6);
4620 Float_t yc0 = helixes[i0].GetHelix(7);
4621 Float_t r0 = helixes[i0].GetHelix(8);
4622 Float_t rc0 = TMath::Sqrt(xc0*xc0+yc0*yc0);
4623 Float_t fi0 = TMath::ATan2(yc0,xc0);
4625 for (Int_t i1=i0+1;i1<nentries;i1++){
4626 AliTPCseed * track1 = (AliTPCseed*)array->At(i1);
4627 if (!track1) continue;
4628 if (TMath::Abs(track1->GetC())<1/kMaxC) continue;
4629 Float_t xc1 = helixes[i1].GetHelix(6);
4630 Float_t yc1 = helixes[i1].GetHelix(7);
4631 Float_t r1 = helixes[i1].GetHelix(8);
4632 Float_t rc1 = TMath::Sqrt(xc1*xc1+yc1*yc1);
4633 Float_t fi1 = TMath::ATan2(yc1,xc1);
4635 Float_t dfi = fi0-fi1;
4638 if (dfi>1.5*TMath::Pi()) dfi-=TMath::Pi(); // take care about edge effect
4639 if (dfi<-1.5*TMath::Pi()) dfi+=TMath::Pi(); //
4640 Float_t dtheta = TMath::Abs(track0->GetTgl()-track1->GetTgl())<TMath::Abs(track0->GetTgl()+track1->GetTgl())? track0->GetTgl()-track1->GetTgl():track0->GetTgl()+track1->GetTgl();
4644 if (track0->GetBConstrain()&&track1->GetBConstrain()) continue; // not constrained
4645 if (track1->GetSigned1Pt()*track0->GetSigned1Pt()>0) continue; // not the same sign
4646 if ( TMath::Abs(track1->GetTgl()+track0->GetTgl())>kMaxdTheta) continue; //distance in the Theta
4647 if ( TMath::Abs(dfi)>kMaxdPhi) continue; //distance in phi
4648 if ( TMath::Sqrt(dfi*dfi+dtheta*dtheta)>kMaxdPhi) continue; //common angular offset
4650 Float_t pt0 = track0->GetSignedPt();
4651 Float_t pt1 = track1->GetSignedPt();
4652 if ((TMath::Abs(pt0+pt1)/(TMath::Abs(pt0)+TMath::Abs(pt1)))>kPtRatio) continue;
4653 if ((iter==1) && TMath::Abs(TMath::Abs(rc0+r0)-TMath::Abs(rc1+r1))>kMaxDeltaRMax) continue;
4654 if ((iter!=1) &&TMath::Abs(TMath::Abs(rc0-r0)-TMath::Abs(rc1-r1))>kMaxDeltaRMin) continue;
4655 if (TMath::Min(TMath::Abs(rc0-r0),TMath::Abs(rc1-r1))<kMinDCAR) continue;
4658 // Now find closest approach
4662 Int_t npoints = helixes[i0].GetRPHIintersections(helixes[i1], phase, radius,10);
4663 if (npoints==0) continue;
4664 helixes[i0].GetClosestPhases(helixes[i1], phase);
4668 Double_t hangles[3];
4669 helixes[i0].Evaluate(phase[0][0],xyz0);
4670 helixes[i1].Evaluate(phase[0][1],xyz1);
4672 helixes[i0].GetAngle(phase[0][0],helixes[i1],phase[0][1],hangles);
4673 Double_t deltah[2],deltabest;
4674 if (TMath::Abs(hangles[2])<kMinAngle) continue;
4678 helixes[i0].ParabolicDCA(helixes[i1],phase[0][0],phase[0][1],radius[0],deltah[0],2);
4680 helixes[i0].ParabolicDCA(helixes[i1],phase[1][0],phase[1][1],radius[1],deltah[1],2);
4681 if (deltah[1]<deltah[0]) ibest=1;
4683 deltabest = TMath::Sqrt(deltah[ibest]);
4684 helixes[i0].Evaluate(phase[ibest][0],xyz0);
4685 helixes[i1].Evaluate(phase[ibest][1],xyz1);
4686 helixes[i0].GetAngle(phase[ibest][0],helixes[i1],phase[ibest][1],hangles);
4687 Double_t radiusbest = TMath::Sqrt(radius[ibest]);
4689 if (deltabest>kMaxDist) continue;
4690 // if (mindcar+mindcaz<40 && (TMath::Abs(hangles[2])<kMinAngle ||deltabest>3)) continue;
4691 Bool_t sign =kFALSE;
4692 if (hangles[2]>kMinAngle) sign =kTRUE;
4695 // circular[i0] = kTRUE;
4696 // circular[i1] = kTRUE;
4697 if (track0->OneOverPt()<track1->OneOverPt()){
4698 track0->SetCircular(track0->GetCircular()+1);
4699 track1->SetCircular(track1->GetCircular()+2);
4702 track1->SetCircular(track1->GetCircular()+1);
4703 track0->SetCircular(track0->GetCircular()+2);
4706 if (AliTPCReconstructor::StreamLevel()>2){
4708 //debug stream to tune "fine" cuts
4709 Int_t lab0=track0->GetLabel();
4710 Int_t lab1=track1->GetLabel();
4711 TTreeSRedirector &cstream = *fDebugStreamer;
4712 cstream<<"Curling2"<<
4728 "npoints="<<npoints<<
4729 "hangles0="<<hangles[0]<<
4730 "hangles1="<<hangles[1]<<
4731 "hangles2="<<hangles[2]<<
4734 "radius="<<radiusbest<<
4735 "deltabest="<<deltabest<<
4736 "phase0="<<phase[ibest][0]<<
4737 "phase1="<<phase[ibest][1]<<
4745 if (AliTPCReconstructor::StreamLevel()>1) {
4746 AliInfo("Time for curling tracks removal");
4755 void AliTPCtrackerMI::FindKinks(TObjArray * array, AliESDEvent *esd)
4762 TObjArray *kinks= new TObjArray(10000);
4763 // TObjArray *v0s= new TObjArray(10000);
4764 Int_t nentries = array->GetEntriesFast();
4765 AliHelix *helixes = new AliHelix[nentries];
4766 Int_t *sign = new Int_t[nentries];
4767 Int_t *nclusters = new Int_t[nentries];
4768 Float_t *alpha = new Float_t[nentries];
4769 AliKink *kink = new AliKink();
4770 Int_t * usage = new Int_t[nentries];
4771 Float_t *zm = new Float_t[nentries];
4772 Float_t *z0 = new Float_t[nentries];
4773 Float_t *fim = new Float_t[nentries];
4774 Float_t *shared = new Float_t[nentries];
4775 Bool_t *circular = new Bool_t[nentries];
4776 Float_t *dca = new Float_t[nentries];
4777 //const AliESDVertex * primvertex = esd->GetVertex();
4779 // nentries = array->GetEntriesFast();
4784 for (Int_t i=0;i<nentries;i++){
4787 AliTPCseed* track = (AliTPCseed*)array->At(i);
4788 if (!track) continue;
4789 track->SetCircular(0);
4791 track->UpdatePoints();
4792 if (( track->GetPoints()[2]- track->GetPoints()[0])>5 && track->GetPoints()[3]>0.8){
4794 nclusters[i]=track->GetNumberOfClusters();
4795 alpha[i] = track->GetAlpha();
4796 new (&helixes[i]) AliHelix(*track);
4798 helixes[i].Evaluate(0,xyz);
4799 sign[i] = (track->GetC()>0) ? -1:1;
4802 if (track->GetProlongation(x,y,z)){
4804 fim[i] = alpha[i]+TMath::ATan2(y,x);
4807 zm[i] = track->GetZ();
4811 circular[i]= kFALSE;
4812 if (track->GetProlongation(0,y,z)) z0[i] = z;
4813 dca[i] = track->GetD(0,0);
4819 Int_t ncandidates =0;
4822 Double_t phase[2][2]={{0,0},{0,0}},radius[2]={0,0};
4825 // Find circling track
4827 for (Int_t i0=0;i0<nentries;i0++){
4828 AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
4829 if (!track0) continue;
4830 if (track0->GetNumberOfClusters()<40) continue;
4831 if (TMath::Abs(1./track0->GetC())>200) continue;
4832 for (Int_t i1=i0+1;i1<nentries;i1++){
4833 AliTPCseed * track1 = (AliTPCseed*)array->At(i1);
4834 if (!track1) continue;
4835 if (track1->GetNumberOfClusters()<40) continue;
4836 if ( TMath::Abs(track1->GetTgl()+track0->GetTgl())>0.1) continue;
4837 if (track0->GetBConstrain()&&track1->GetBConstrain()) continue;
4838 if (TMath::Abs(1./track1->GetC())>200) continue;
4839 if (track1->GetSigned1Pt()*track0->GetSigned1Pt()>0) continue;
4840 if (track1->GetTgl()*track0->GetTgl()>0) continue;
4841 if (TMath::Max(TMath::Abs(1./track0->GetC()),TMath::Abs(1./track1->GetC()))>190) continue;
4842 if (track0->GetBConstrain()&&track1->OneOverPt()<track0->OneOverPt()) continue; //returning - lower momenta
4843 if (track1->GetBConstrain()&&track0->OneOverPt()<track1->OneOverPt()) continue; //returning - lower momenta
4845 Float_t mindcar = TMath::Min(TMath::Abs(dca[i0]),TMath::Abs(dca[i1]));
4846 if (mindcar<5) continue;
4847 Float_t mindcaz = TMath::Min(TMath::Abs(z0[i0]-GetZ()),TMath::Abs(z0[i1]-GetZ()));
4848 if (mindcaz<5) continue;
4849 if (mindcar+mindcaz<20) continue;
4852 Float_t xc0 = helixes[i0].GetHelix(6);
4853 Float_t yc0 = helixes[i0].GetHelix(7);
4854 Float_t r0 = helixes[i0].GetHelix(8);
4855 Float_t xc1 = helixes[i1].GetHelix(6);
4856 Float_t yc1 = helixes[i1].GetHelix(7);
4857 Float_t r1 = helixes[i1].GetHelix(8);
4859 Float_t rmean = (r0+r1)*0.5;
4860 Float_t delta =TMath::Sqrt((xc1-xc0)*(xc1-xc0)+(yc1-yc0)*(yc1-yc0));
4861 //if (delta>30) continue;
4862 if (delta>rmean*0.25) continue;
4863 if (TMath::Abs(r0-r1)/rmean>0.3) continue;
4865 Int_t npoints = helixes[i0].GetRPHIintersections(helixes[i1], phase, radius,10);
4866 if (npoints==0) continue;
4867 helixes[i0].GetClosestPhases(helixes[i1], phase);
4871 Double_t hangles[3];
4872 helixes[i0].Evaluate(phase[0][0],xyz0);
4873 helixes[i1].Evaluate(phase[0][1],xyz1);
4875 helixes[i0].GetAngle(phase[0][0],helixes[i1],phase[0][1],hangles);
4876 Double_t deltah[2],deltabest;
4877 if (hangles[2]<2.8) continue;
4880 helixes[i0].ParabolicDCA(helixes[i1],phase[0][0],phase[0][1],radius[0],deltah[0],2);
4882 helixes[i0].ParabolicDCA(helixes[i1],phase[1][0],phase[1][1],radius[1],deltah[1],2);
4883 if (deltah[1]<deltah[0]) ibest=1;
4885 deltabest = TMath::Sqrt(deltah[ibest]);
4886 helixes[i0].Evaluate(phase[ibest][0],xyz0);
4887 helixes[i1].Evaluate(phase[ibest][1],xyz1);
4888 helixes[i0].GetAngle(phase[ibest][0],helixes[i1],phase[ibest][1],hangles);
4889 Double_t radiusbest = TMath::Sqrt(radius[ibest]);
4891 if (deltabest>6) continue;
4892 if (mindcar+mindcaz<40 && (hangles[2]<3.12||deltabest>3)) continue;
4893 Bool_t lsign =kFALSE;
4894 if (hangles[2]>3.06) lsign =kTRUE;
4897 circular[i0] = kTRUE;
4898 circular[i1] = kTRUE;
4899 if (track0->OneOverPt()<track1->OneOverPt()){
4900 track0->SetCircular(track0->GetCircular()+1);
4901 track1->SetCircular(track1->GetCircular()+2);
4904 track1->SetCircular(track1->GetCircular()+1);
4905 track0->SetCircular(track0->GetCircular()+2);
4908 if (lsign&&AliTPCReconstructor::StreamLevel()>1){
4910 Int_t lab0=track0->GetLabel();
4911 Int_t lab1=track1->GetLabel();
4912 TTreeSRedirector &cstream = *fDebugStreamer;
4913 cstream<<"Curling"<<
4920 "mindcar="<<mindcar<<
4921 "mindcaz="<<mindcaz<<
4924 "npoints="<<npoints<<
4925 "hangles0="<<hangles[0]<<
4926 "hangles2="<<hangles[2]<<
4931 "radius="<<radiusbest<<
4932 "deltabest="<<deltabest<<
4933 "phase0="<<phase[ibest][0]<<
4934 "phase1="<<phase[ibest][1]<<
4944 for (Int_t i =0;i<nentries;i++){
4945 if (sign[i]==0) continue;
4946 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
4953 Double_t cradius0 = 40*40;
4954 Double_t cradius1 = 270*270;
4957 Double_t cdist3=0.55;
4958 for (Int_t j =i+1;j<nentries;j++){
4960 if (sign[j]*sign[i]<1) continue;
4961 if ( (nclusters[i]+nclusters[j])>200) continue;
4962 if ( (nclusters[i]+nclusters[j])<80) continue;
4963 if ( TMath::Abs(zm[i]-zm[j])>60.) continue;
4964 if ( TMath::Abs(fim[i]-fim[j])>0.6 && TMath::Abs(fim[i]-fim[j])<5.7 ) continue;
4965 //AliTPCseed * track1 = (AliTPCseed*)array->At(j); Double_t phase[2][2],radius[2];
4966 Int_t npoints = helixes[i].GetRPHIintersections(helixes[j], phase, radius,20);
4967 if (npoints<1) continue;
4970 if (radius[0]<cradius0||radius[0]>cradius1) continue;
4973 if ( (radius[0]<cradius0||radius[0]>cradius1) && (radius[1]<cradius0||radius[1]>cradius1) ) continue;
4976 Double_t delta1=10000,delta2=10000;
4977 // cuts on the intersection radius
4978 helixes[i].LinearDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta1);
4979 if (radius[0]<20&&delta1<1) continue; //intersection at vertex
4980 if (radius[0]<10&&delta1<3) continue; //intersection at vertex
4982 helixes[i].LinearDCA(helixes[j],phase[1][0],phase[1][1],radius[1],delta2);
4983 if (radius[1]<20&&delta2<1) continue; //intersection at vertex
4984 if (radius[1]<10&&delta2<3) continue; //intersection at vertex
4987 Double_t distance1 = TMath::Min(delta1,delta2);
4988 if (distance1>cdist1) continue; // cut on DCA linear approximation
4990 npoints = helixes[i].GetRPHIintersections(helixes[j], phase, radius,20);
4991 helixes[i].ParabolicDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta1);
4992 if (radius[0]<20&&delta1<1) continue; //intersection at vertex
4993 if (radius[0]<10&&delta1<3) continue; //intersection at vertex
4996 helixes[i].ParabolicDCA(helixes[j],phase[1][0],phase[1][1],radius[1],delta2);
4997 if (radius[1]<20&&delta2<1) continue; //intersection at vertex
4998 if (radius[1]<10&&delta2<3) continue; //intersection at vertex
5000 distance1 = TMath::Min(delta1,delta2);
5003 rkink = TMath::Sqrt(radius[0]);
5006 rkink = TMath::Sqrt(radius[1]);
5008 if (distance1>cdist2) continue;
5011 AliTPCseed * track1 = (AliTPCseed*)array->At(j);
5014 Int_t row0 = GetRowNumber(rkink);
5015 if (row0<10) continue;
5016 if (row0>150) continue;
5019 Float_t dens00=-1,dens01=-1;
5020 Float_t dens10=-1,dens11=-1;
5022 Int_t found,foundable,ishared;
5023 track0->GetClusterStatistic(0,row0-5, found, foundable,ishared,kFALSE);
5024 if (foundable>5) dens00 = Float_t(found)/Float_t(foundable);
5025 track0->GetClusterStatistic(row0+5,155, found, foundable,ishared,kFALSE);
5026 if (foundable>5) dens01 = Float_t(found)/Float_t(foundable);
5028 track1->GetClusterStatistic(0,row0-5, found, foundable,ishared,kFALSE);
5029 if (foundable>10) dens10 = Float_t(found)/Float_t(foundable);
5030 track1->GetClusterStatistic(row0+5,155, found, foundable,ishared,kFALSE);
5031 if (foundable>10) dens11 = Float_t(found)/Float_t(foundable);
5033 if (dens00<dens10 && dens01<dens11) continue;
5034 if (dens00>dens10 && dens01>dens11) continue;
5035 if (TMath::Max(dens00,dens10)<0.1) continue;
5036 if (TMath::Max(dens01,dens11)<0.3) continue;
5038 if (TMath::Min(dens00,dens10)>0.6) continue;
5039 if (TMath::Min(dens01,dens11)>0.6) continue;
5042 AliTPCseed * ktrack0, *ktrack1;
5051 if (TMath::Abs(ktrack0->GetC())>5) continue; // cut on the curvature for mother particle
5052 AliExternalTrackParam paramm(*ktrack0);
5053 AliExternalTrackParam paramd(*ktrack1);
5054 if (row0>60&&ktrack1->GetReference().GetX()>90.)new (¶md) AliExternalTrackParam(ktrack1->GetReference());
5057 kink->SetMother(paramm);
5058 kink->SetDaughter(paramd);
5061 Float_t x[3] = { kink->GetPosition()[0],kink->GetPosition()[1],kink->GetPosition()[2]};
5063 fkParam->Transform0to1(x,index);
5064 fkParam->Transform1to2(x,index);
5065 row0 = GetRowNumber(x[0]);
5067 if (kink->GetR()<100) continue;
5068 if (kink->GetR()>240) continue;
5069 if (kink->GetPosition()[2]/kink->GetR()>AliTPCReconstructor::GetCtgRange()) continue; //out of fiducial volume
5070 if (kink->GetDistance()>cdist3) continue;
5071 Float_t dird = kink->GetDaughterP()[0]*kink->GetPosition()[0]+kink->GetDaughterP()[1]*kink->GetPosition()[1]; // rough direction estimate
5072 if (dird<0) continue;
5074 Float_t dirm = kink->GetMotherP()[0]*kink->GetPosition()[0]+kink->GetMotherP()[1]*kink->GetPosition()[1]; // rough direction estimate
5075 if (dirm<0) continue;
5076 Float_t mpt = TMath::Sqrt(kink->GetMotherP()[0]*kink->GetMotherP()[0]+kink->GetMotherP()[1]*kink->GetMotherP()[1]);
5077 if (mpt<0.2) continue;
5080 //for high momenta momentum not defined well in first iteration
5081 Double_t qt = TMath::Sin(kink->GetAngle(2))*ktrack1->GetP();
5082 if (qt>0.35) continue;
5085 kink->SetLabel(CookLabel(ktrack0,0.4,0,row0),0);
5086 kink->SetLabel(CookLabel(ktrack1,0.4,row0,160),1);
5088 kink->SetTPCDensity(dens00,0,0);
5089 kink->SetTPCDensity(dens01,0,1);
5090 kink->SetTPCDensity(dens10,1,0);
5091 kink->SetTPCDensity(dens11,1,1);
5092 kink->SetIndex(i,0);
5093 kink->SetIndex(j,1);
5096 kink->SetTPCDensity(dens10,0,0);
5097 kink->SetTPCDensity(dens11,0,1);
5098 kink->SetTPCDensity(dens00,1,0);
5099 kink->SetTPCDensity(dens01,1,1);
5100 kink->SetIndex(j,0);
5101 kink->SetIndex(i,1);
5104 if (mpt<1||kink->GetAngle(2)>0.1){
5105 // angle and densities not defined yet
5106 if (kink->GetTPCDensityFactor()<0.8) continue;
5107 if ((2-kink->GetTPCDensityFactor())*kink->GetDistance() >0.25) continue;
5108 if (kink->GetAngle(2)*ktrack0->GetP()<0.003) continue; //too small angle
5109 if (kink->GetAngle(2)>0.2&&kink->GetTPCDensityFactor()<1.15) continue;
5110 if (kink->GetAngle(2)>0.2&&kink->GetTPCDensity(0,1)>0.05) continue;
5112 Float_t criticalangle = track0->GetSigmaSnp2()+track0->GetSigmaTgl2();
5113 criticalangle+= track1->GetSigmaSnp2()+track1->GetSigmaTgl2();
5114 criticalangle= 3*TMath::Sqrt(criticalangle);
5115 if (criticalangle>0.02) criticalangle=0.02;
5116 if (kink->GetAngle(2)<criticalangle) continue;
5119 Int_t drow = Int_t(2.+0.5/(0.05+kink->GetAngle(2))); // overlap region defined
5120 Float_t shapesum =0;
5122 for ( Int_t row = row0-drow; row<row0+drow;row++){
5123 if (row<0) continue;
5124 if (row>155) continue;
5125 if (ktrack0->GetClusterPointer(row)){
5126 AliTPCTrackerPoint *point =ktrack0->GetTrackPoint(row);
5127 shapesum+=point->GetSigmaY()+point->GetSigmaZ();
5130 if (ktrack1->GetClusterPointer(row)){
5131 AliTPCTrackerPoint *point =ktrack1->GetTrackPoint(row);
5132 shapesum+=point->GetSigmaY()+point->GetSigmaZ();
5137 kink->SetShapeFactor(-1.);
5140 kink->SetShapeFactor(shapesum/sum);
5142 // esd->AddKink(kink);
5144 // kink->SetMother(paramm);
5145 //kink->SetDaughter(paramd);
5147 Double_t chi2P2 = paramm.GetParameter()[2]-paramd.GetParameter()[2];
5149 chi2P2/=paramm.GetCovariance()[5]+paramd.GetCovariance()[5];
5150 Double_t chi2P3 = paramm.GetParameter()[3]-paramd.GetParameter()[3];
5152 chi2P3/=paramm.GetCovariance()[9]+paramd.GetCovariance()[9];
5154 if (AliTPCReconstructor::StreamLevel()>1) {
5155 (*fDebugStreamer)<<"kinkLpt"<<
5163 if ( chi2P2+chi2P3<AliTPCReconstructor::GetRecoParam()->GetKinkAngleCutChi2(0)){
5167 kinks->AddLast(kink);
5173 // sort the kinks according quality - and refit them towards vertex
5175 Int_t nkinks = kinks->GetEntriesFast();
5176 Float_t *quality = new Float_t[nkinks];
5177 Int_t *indexes = new Int_t[nkinks];
5178 AliTPCseed *mothers = new AliTPCseed[nkinks];
5179 AliTPCseed *daughters = new AliTPCseed[nkinks];
5182 for (Int_t i=0;i<nkinks;i++){
5184 AliKink *kinkl = (AliKink*)kinks->At(i);
5186 // refit kinks towards vertex
5188 Int_t index0 = kinkl->GetIndex(0);
5189 Int_t index1 = kinkl->GetIndex(1);
5190 AliTPCseed * ktrack0 = (AliTPCseed*)array->At(index0);
5191 AliTPCseed * ktrack1 = (AliTPCseed*)array->At(index1);
5193 Int_t sumn=ktrack0->GetNumberOfClusters()+ktrack1->GetNumberOfClusters();
5195 // Refit Kink under if too small angle
5197 if (kinkl->GetAngle(2)<0.05){
5198 kinkl->SetTPCRow0(GetRowNumber(kinkl->GetR()));
5199 Int_t row0 = kinkl->GetTPCRow0();
5200 Int_t drow = Int_t(2.+0.5/(0.05+kinkl->GetAngle(2)));
5203 Int_t last = row0-drow;
5204 if (last<40) last=40;
5205 if (last<ktrack0->GetFirstPoint()+25) last = ktrack0->GetFirstPoint()+25;
5206 AliTPCseed* seed0 = ReSeed(ktrack0,last,kFALSE);
5209 Int_t first = row0+drow;
5210 if (first>130) first=130;
5211 if (first>ktrack1->GetLastPoint()-25) first = TMath::Max(ktrack1->GetLastPoint()-25,30);
5212 AliTPCseed* seed1 = ReSeed(ktrack1,first,kTRUE);
5214 if (seed0 && seed1){
5215 kinkl->SetStatus(1,8);
5216 if (RefitKink(*seed0,*seed1,*kinkl)) kinkl->SetStatus(1,9);
5217 row0 = GetRowNumber(kinkl->GetR());
5218 sumn = seed0->GetNumberOfClusters()+seed1->GetNumberOfClusters();
5219 mothers[i] = *seed0;
5220 daughters[i] = *seed1;
5223 delete kinks->RemoveAt(i);
5224 if (seed0) delete seed0;
5225 if (seed1) delete seed1;
5228 if (kinkl->GetDistance()>0.5 || kinkl->GetR()<110 || kinkl->GetR()>240) {
5229 delete kinks->RemoveAt(i);
5230 if (seed0) delete seed0;
5231 if (seed1) delete seed1;
5239 if (kinkl) quality[i] = 160*((0.1+kinkl->GetDistance())*(2.-kinkl->GetTPCDensityFactor()))/(sumn+40.); //the longest -clossest will win
5241 TMath::Sort(nkinks,quality,indexes,kFALSE);
5243 //remove double find kinks
5245 for (Int_t ikink0=1;ikink0<nkinks;ikink0++){
5246 AliKink * kink0 = (AliKink*) kinks->At(indexes[ikink0]);
5247 if (!kink0) continue;
5249 for (Int_t ikink1=0;ikink1<ikink0;ikink1++){
5250 kink0 = (AliKink*) kinks->At(indexes[ikink0]);
5251 if (!kink0) continue;
5252 AliKink * kink1 = (AliKink*) kinks->At(indexes[ikink1]);
5253 if (!kink1) continue;
5254 // if not close kink continue
5255 if (TMath::Abs(kink1->GetPosition()[2]-kink0->GetPosition()[2])>10) continue;
5256 if (TMath::Abs(kink1->GetPosition()[1]-kink0->GetPosition()[1])>10) continue;
5257 if (TMath::Abs(kink1->GetPosition()[0]-kink0->GetPosition()[0])>10) continue;
5259 AliTPCseed &mother0 = mothers[indexes[ikink0]];
5260 AliTPCseed &daughter0 = daughters[indexes[ikink0]];
5261 AliTPCseed &mother1 = mothers[indexes[ikink1]];
5262 AliTPCseed &daughter1 = daughters[indexes[ikink1]];
5263 Int_t row0 = (kink0->GetTPCRow0()+kink1->GetTPCRow0())/2;
5272 for (Int_t i=0;i<row0;i++){
5273 if (mother0.GetClusterIndex(i)>0 && mother1.GetClusterIndex(i)>0){
5276 if (mother0.GetClusterIndex(i)==mother1.GetClusterIndex(i)){
5283 for (Int_t i=row0;i<158;i++){
5284 if (daughter0.GetClusterIndex(i)>0 && daughter0.GetClusterIndex(i)>0){
5287 if (mother0.GetClusterIndex(i)==mother1.GetClusterIndex(i)){
5293 Float_t ratio = Float_t(same+1)/Float_t(both+1);
5294 Float_t ratiom = Float_t(samem+1)/Float_t(bothm+1);
5295 Float_t ratiod = Float_t(samed+1)/Float_t(bothd+1);
5296 if (ratio>0.3 && ratiom>0.5 &&ratiod>0.5) {
5297 Int_t sum0 = mother0.GetNumberOfClusters()+daughter0.GetNumberOfClusters();
5298 Int_t sum1 = mother1.GetNumberOfClusters()+daughter1.GetNumberOfClusters();
5300 shared[kink0->GetIndex(0)]= kTRUE;
5301 shared[kink0->GetIndex(1)]= kTRUE;
5302 delete kinks->RemoveAt(indexes[ikink0]);
5306 shared[kink1->GetIndex(0)]= kTRUE;
5307 shared[kink1->GetIndex(1)]= kTRUE;
5308 delete kinks->RemoveAt(indexes[ikink1]);
5315 for (Int_t i=0;i<nkinks;i++){
5316 AliKink * kinkl = (AliKink*) kinks->At(indexes[i]);
5317 if (!kinkl) continue;
5318 kinkl->SetTPCRow0(GetRowNumber(kinkl->GetR()));
5319 Int_t index0 = kinkl->GetIndex(0);
5320 Int_t index1 = kinkl->GetIndex(1);
5321 if (circular[index0]||(circular[index1]&&kinkl->GetDistance()>0.2)) continue;
5322 kinkl->SetMultiple(usage[index0],0);
5323 kinkl->SetMultiple(usage[index1],1);
5324 if (kinkl->GetMultiple()[0]+kinkl->GetMultiple()[1]>2) continue;
5325 if (kinkl->GetMultiple()[0]+kinkl->GetMultiple()[1]>0 && quality[indexes[i]]>0.2) continue;
5326 if (kinkl->GetMultiple()[0]+kinkl->GetMultiple()[1]>0 && kinkl->GetDistance()>0.2) continue;
5327 if (circular[index0]||(circular[index1]&&kinkl->GetDistance()>0.1)) continue;
5329 AliTPCseed * ktrack0 = (AliTPCseed*)array->At(index0);
5330 AliTPCseed * ktrack1 = (AliTPCseed*)array->At(index1);
5331 if (!ktrack0 || !ktrack1) continue;
5332 Int_t index = esd->AddKink(kinkl);
5335 if ( ktrack0->GetKinkIndex(0)==0 && ktrack1->GetKinkIndex(0)==0) { //best kink
5336 if (mothers[indexes[i]].GetNumberOfClusters()>20 && daughters[indexes[i]].GetNumberOfClusters()>20 && (mothers[indexes[i]].GetNumberOfClusters()+daughters[indexes[i]].GetNumberOfClusters())>100){
5337 *ktrack0 = mothers[indexes[i]];
5338 *ktrack1 = daughters[indexes[i]];
5342 ktrack0->SetKinkIndex(usage[index0],-(index+1));
5343 ktrack1->SetKinkIndex(usage[index1], (index+1));
5348 // Remove tracks corresponding to shared kink's
5350 for (Int_t i=0;i<nentries;i++){
5351 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5352 if (!track0) continue;
5353 if (track0->GetKinkIndex(0)!=0) continue;
5354 if (shared[i]) delete array->RemoveAt(i);
5359 RemoveUsed2(array,0.5,0.4,30);
5361 for (Int_t i=0;i<nentries;i++){
5362 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5363 if (!track0) continue;
5364 track0->CookdEdx(0.02,0.6);
5368 for (Int_t i=0;i<nentries;i++){
5369 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5370 if (!track0) continue;
5371 if (track0->Pt()<1.4) continue;
5372 //remove double high momenta tracks - overlapped with kink candidates
5375 for (Int_t icl=track0->GetFirstPoint();icl<track0->GetLastPoint(); icl++){
5376 if (track0->GetClusterPointer(icl)!=0){
5378 if (track0->GetClusterPointer(icl)->IsUsed(10)) ishared++;
5381 if (Float_t(ishared+1)/Float_t(all+1)>0.5) {
5382 delete array->RemoveAt(i);
5386 if (track0->GetKinkIndex(0)!=0) continue;
5387 if (track0->GetNumberOfClusters()<80) continue;
5389 AliTPCseed *pmother = new AliTPCseed();
5390 AliTPCseed *pdaughter = new AliTPCseed();
5391 AliKink *pkink = new AliKink;
5393 AliTPCseed & mother = *pmother;
5394 AliTPCseed & daughter = *pdaughter;
5395 AliKink & kinkl = *pkink;
5396 if (CheckKinkPoint(track0,mother,daughter, kinkl)){
5397 if (mother.GetNumberOfClusters()<30||daughter.GetNumberOfClusters()<20) {
5401 continue; //too short tracks
5403 if (mother.Pt()<1.4) {
5409 Int_t row0= kinkl.GetTPCRow0();
5410 if (kinkl.GetDistance()>0.5 || kinkl.GetR()<110. || kinkl.GetR()>240.) {
5417 Int_t index = esd->AddKink(&kinkl);
5418 mother.SetKinkIndex(0,-(index+1));
5419 daughter.SetKinkIndex(0,index+1);
5420 if (mother.GetNumberOfClusters()>50) {
5421 delete array->RemoveAt(i);
5422 array->AddAt(new AliTPCseed(mother),i);
5425 array->AddLast(new AliTPCseed(mother));
5427 array->AddLast(new AliTPCseed(daughter));
5428 for (Int_t icl=0;icl<row0;icl++) {
5429 if (mother.GetClusterPointer(icl)) mother.GetClusterPointer(icl)->Use(20);
5432 for (Int_t icl=row0;icl<158;icl++) {
5433 if (daughter.GetClusterPointer(icl)) daughter.GetClusterPointer(icl)->Use(20);
5442 delete [] daughters;
5464 AliInfo(Form("Ncandidates=\t%d\t%d\t%d\t%d\n",esd->GetNumberOfKinks(),ncandidates,ntracks,nall));
5469 Int_t AliTPCtrackerMI::RefitKink(AliTPCseed &mother, AliTPCseed &daughter, const AliESDkink &knk)
5472 // refit kink towards to the vertex
5475 AliKink &kink=(AliKink &)knk;
5477 Int_t row0 = GetRowNumber(kink.GetR());
5478 FollowProlongation(mother,0);
5479 mother.Reset(kFALSE);
5481 FollowProlongation(daughter,row0);
5482 daughter.Reset(kFALSE);
5483 FollowBackProlongation(daughter,158);
5484 daughter.Reset(kFALSE);
5485 Int_t first = TMath::Max(row0-20,30);
5486 Int_t last = TMath::Min(row0+20,140);
5488 const Int_t kNdiv =5;
5489 AliTPCseed param0[kNdiv]; // parameters along the track
5490 AliTPCseed param1[kNdiv]; // parameters along the track
5491 AliKink kinks[kNdiv]; // corresponding kink parameters
5494 for (Int_t irow=0; irow<kNdiv;irow++){
5495 rows[irow] = first +((last-first)*irow)/(kNdiv-1);
5497 // store parameters along the track
5499 for (Int_t irow=0;irow<kNdiv;irow++){
5500 FollowBackProlongation(mother, rows[irow]);
5501 FollowProlongation(daughter,rows[kNdiv-1-irow]);
5502 param0[irow] = mother;
5503 param1[kNdiv-1-irow] = daughter;
5507 for (Int_t irow=0; irow<kNdiv-1;irow++){
5508 if (param0[irow].GetNumberOfClusters()<kNdiv||param1[irow].GetNumberOfClusters()<kNdiv) continue;
5509 kinks[irow].SetMother(param0[irow]);
5510 kinks[irow].SetDaughter(param1[irow]);
5511 kinks[irow].Update();
5514 // choose kink with best "quality"
5516 Double_t mindist = 10000;
5517 for (Int_t irow=0;irow<kNdiv;irow++){
5518 if (param0[irow].GetNumberOfClusters()<20||param1[irow].GetNumberOfClusters()<20) continue;
5519 if (TMath::Abs(kinks[irow].GetR())>240.) continue;
5520 if (TMath::Abs(kinks[irow].GetR())<100.) continue;
5522 Float_t normdist = TMath::Abs(param0[irow].GetX()-kinks[irow].GetR())*(0.1+kink.GetDistance());
5523 normdist/= (param0[irow].GetNumberOfClusters()+param1[irow].GetNumberOfClusters()+40.);
5524 if (normdist < mindist){
5530 if (index==-1) return 0;
5533 param0[index].Reset(kTRUE);
5534 FollowProlongation(param0[index],0);
5536 mother = param0[index];
5537 daughter = param1[index]; // daughter in vertex
5539 kink.SetMother(mother);
5540 kink.SetDaughter(daughter);
5542 kink.SetTPCRow0(GetRowNumber(kink.GetR()));
5543 kink.SetTPCncls(param0[index].GetNumberOfClusters(),0);
5544 kink.SetTPCncls(param1[index].GetNumberOfClusters(),1);
5545 kink.SetLabel(CookLabel(&mother,0.4, 0,kink.GetTPCRow0()),0);
5546 kink.SetLabel(CookLabel(&daughter,0.4, kink.GetTPCRow0(),160),1);
5547 mother.SetLabel(kink.GetLabel(0));
5548 daughter.SetLabel(kink.GetLabel(1));
5554 void AliTPCtrackerMI::UpdateKinkQualityM(AliTPCseed * seed){
5556 // update Kink quality information for mother after back propagation
5558 if (seed->GetKinkIndex(0)>=0) return;
5559 for (Int_t ikink=0;ikink<3;ikink++){
5560 Int_t index = seed->GetKinkIndex(ikink);
5561 if (index>=0) break;
5562 index = TMath::Abs(index)-1;
5563 AliESDkink * kink = fEvent->GetKink(index);
5564 kink->SetTPCDensity(-1,0,0);
5565 kink->SetTPCDensity(1,0,1);
5567 Int_t row0 = kink->GetTPCRow0() - 2 - Int_t( 0.5/ (0.05+kink->GetAngle(2)));
5568 if (row0<15) row0=15;
5570 Int_t row1 = kink->GetTPCRow0() + 2 + Int_t( 0.5/ (0.05+kink->GetAngle(2)));
5571 if (row1>145) row1=145;
5573 Int_t found,foundable,shared;
5574 seed->GetClusterStatistic(0,row0, found, foundable,shared,kFALSE);
5575 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),0,0);
5576 seed->GetClusterStatistic(row1,155, found, foundable,shared,kFALSE);
5577 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),0,1);
5582 void AliTPCtrackerMI::UpdateKinkQualityD(AliTPCseed * seed){
5584 // update Kink quality information for daughter after refit
5586 if (seed->GetKinkIndex(0)<=0) return;
5587 for (Int_t ikink=0;ikink<3;ikink++){
5588 Int_t index = seed->GetKinkIndex(ikink);
5589 if (index<=0) break;
5590 index = TMath::Abs(index)-1;
5591 AliESDkink * kink = fEvent->GetKink(index);
5592 kink->SetTPCDensity(-1,1,0);
5593 kink->SetTPCDensity(-1,1,1);
5595 Int_t row0 = kink->GetTPCRow0() -2 - Int_t( 0.5/ (0.05+kink->GetAngle(2)));
5596 if (row0<15) row0=15;
5598 Int_t row1 = kink->GetTPCRow0() +2 + Int_t( 0.5/ (0.05+kink->GetAngle(2)));
5599 if (row1>145) row1=145;
5601 Int_t found,foundable,shared;
5602 seed->GetClusterStatistic(0,row0, found, foundable,shared,kFALSE);
5603 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),1,0);
5604 seed->GetClusterStatistic(row1,155, found, foundable,shared,kFALSE);
5605 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),1,1);
5611 Int_t AliTPCtrackerMI::CheckKinkPoint(AliTPCseed*seed,AliTPCseed &mother, AliTPCseed &daughter, const AliESDkink &knk)
5614 // check kink point for given track
5615 // if return value=0 kink point not found
5616 // otherwise seed0 correspond to mother particle
5617 // seed1 correspond to daughter particle
5618 // kink parameter of kink point
5619 AliKink &kink=(AliKink &)knk;
5621 Int_t middlerow = (seed->GetFirstPoint()+seed->GetLastPoint())/2;
5622 Int_t first = seed->GetFirstPoint();
5623 Int_t last = seed->GetLastPoint();
5624 if (last-first<20) return 0; // shortest length - 2*30 = 60 pad-rows
5627 AliTPCseed *seed1 = ReSeed(seed,middlerow+20, kTRUE); //middle of chamber
5628 if (!seed1) return 0;
5629 FollowProlongation(*seed1,seed->GetLastPoint()-20);
5630 seed1->Reset(kTRUE);
5631 FollowProlongation(*seed1,158);
5632 seed1->Reset(kTRUE);
5633 last = seed1->GetLastPoint();
5635 AliTPCseed *seed0 = new AliTPCseed(*seed);
5636 seed0->Reset(kFALSE);
5639 AliTPCseed param0[20]; // parameters along the track
5640 AliTPCseed param1[20]; // parameters along the track
5641 AliKink kinks[20]; // corresponding kink parameters
5643 for (Int_t irow=0; irow<20;irow++){
5644 rows[irow] = first +((last-first)*irow)/19;
5646 // store parameters along the track
5648 for (Int_t irow=0;irow<20;irow++){
5649 FollowBackProlongation(*seed0, rows[irow]);
5650 FollowProlongation(*seed1,rows[19-irow]);
5651 param0[irow] = *seed0;
5652 param1[19-irow] = *seed1;
5656 for (Int_t irow=0; irow<19;irow++){
5657 kinks[irow].SetMother(param0[irow]);
5658 kinks[irow].SetDaughter(param1[irow]);
5659 kinks[irow].Update();
5662 // choose kink with biggest change of angle
5664 Double_t maxchange= 0;
5665 for (Int_t irow=1;irow<19;irow++){
5666 if (TMath::Abs(kinks[irow].GetR())>240.) continue;
5667 if (TMath::Abs(kinks[irow].GetR())<110.) continue;
5668 Float_t quality = TMath::Abs(kinks[irow].GetAngle(2))/(3.+TMath::Abs(kinks[irow].GetR()-param0[irow].GetX()));
5669 if ( quality > maxchange){
5670 maxchange = quality;
5677 if (index<0) return 0;
5679 Int_t row0 = GetRowNumber(kinks[index].GetR()); //row 0 estimate
5680 seed0 = new AliTPCseed(param0[index]);
5681 seed1 = new AliTPCseed(param1[index]);
5682 seed0->Reset(kFALSE);
5683 seed1->Reset(kFALSE);
5684 seed0->ResetCovariance(10.);
5685 seed1->ResetCovariance(10.);
5686 FollowProlongation(*seed0,0);
5687 FollowBackProlongation(*seed1,158);
5688 mother = *seed0; // backup mother at position 0
5689 seed0->Reset(kFALSE);
5690 seed1->Reset(kFALSE);
5691 seed0->ResetCovariance(10.);
5692 seed1->ResetCovariance(10.);
5694 first = TMath::Max(row0-20,0);
5695 last = TMath::Min(row0+20,158);
5697 for (Int_t irow=0; irow<20;irow++){
5698 rows[irow] = first +((last-first)*irow)/19;
5700 // store parameters along the track
5702 for (Int_t irow=0;irow<20;irow++){
5703 FollowBackProlongation(*seed0, rows[irow]);
5704 FollowProlongation(*seed1,rows[19-irow]);
5705 param0[irow] = *seed0;
5706 param1[19-irow] = *seed1;
5710 for (Int_t irow=0; irow<19;irow++){
5711 kinks[irow].SetMother(param0[irow]);
5712 kinks[irow].SetDaughter(param1[irow]);
5713 // param0[irow].Dump();
5714 //param1[irow].Dump();
5715 kinks[irow].Update();
5718 // choose kink with biggest change of angle
5721 for (Int_t irow=0;irow<20;irow++){
5722 if (TMath::Abs(kinks[irow].GetR())>250.) continue;
5723 if (TMath::Abs(kinks[irow].GetR())<90.) continue;
5724 Float_t quality = TMath::Abs(kinks[irow].GetAngle(2))/(3.+TMath::Abs(kinks[irow].GetR()-param0[irow].GetX()));
5725 if ( quality > maxchange){
5726 maxchange = quality;
5733 if (index==-1 || param0[index].GetNumberOfClusters()+param1[index].GetNumberOfClusters()<100){
5739 // Float_t anglesigma = TMath::Sqrt(param0[index].fC22+param0[index].fC33+param1[index].fC22+param1[index].fC33);
5741 kink.SetMother(param0[index]);
5742 kink.SetDaughter(param1[index]);
5745 Double_t chi2P2 = param0[index].GetParameter()[2]-param1[index].GetParameter()[2];
5747 chi2P2/=param0[index].GetCovariance()[5]+param1[index].GetCovariance()[5];
5748 Double_t chi2P3 = param0[index].GetParameter()[3]-param1[index].GetParameter()[3];
5750 chi2P3/=param0[index].GetCovariance()[9]+param1[index].GetCovariance()[9];
5752 if (AliTPCReconstructor::StreamLevel()>1) {
5753 (*fDebugStreamer)<<"kinkHpt"<<
5756 "p0.="<<¶m0[index]<<
5757 "p1.="<<¶m1[index]<<
5761 if ( chi2P2+chi2P3<AliTPCReconstructor::GetRecoParam()->GetKinkAngleCutChi2(0)){
5768 row0 = GetRowNumber(kink.GetR());
5769 kink.SetTPCRow0(row0);
5770 kink.SetLabel(CookLabel(seed0,0.5,0,row0),0);
5771 kink.SetLabel(CookLabel(seed1,0.5,row0,158),1);
5772 kink.SetIndex(-10,0);
5773 kink.SetIndex(int(param0[index].GetNumberOfClusters()+param1[index].GetNumberOfClusters()),1);
5774 kink.SetTPCncls(param0[index].GetNumberOfClusters(),0);
5775 kink.SetTPCncls(param1[index].GetNumberOfClusters(),1);
5778 // new (&mother) AliTPCseed(param0[index]);
5779 daughter = param1[index];
5780 daughter.SetLabel(kink.GetLabel(1));
5781 param0[index].Reset(kTRUE);
5782 FollowProlongation(param0[index],0);
5783 mother = param0[index];
5784 mother.SetLabel(kink.GetLabel(0));
5785 if ( chi2P2+chi2P3<AliTPCReconstructor::GetRecoParam()->GetKinkAngleCutChi2(1)){
5797 AliTPCseed* AliTPCtrackerMI::ReSeed(AliTPCseed *t)
5800 // reseed - refit - track
5803 // Int_t last = fSectors->GetNRows()-1;
5805 if (fSectors == fOuterSec){
5806 first = TMath::Max(first, t->GetFirstPoint()-fInnerSec->GetNRows());
5810 first = t->GetFirstPoint();
5812 AliTPCseed * seed = MakeSeed(t,0.1,0.5,0.9);
5813 FollowBackProlongation(*t,fSectors->GetNRows()-1);
5815 FollowProlongation(*t,first);
5825 //_____________________________________________________________________________
5826 Int_t AliTPCtrackerMI::ReadSeeds(const TFile *inp) {
5827 //-----------------------------------------------------------------
5828 // This function reades track seeds.
5829 //-----------------------------------------------------------------
5830 TDirectory *savedir=gDirectory;
5832 TFile *in=(TFile*)inp;
5833 if (!in->IsOpen()) {
5834 cerr<<"AliTPCtrackerMI::ReadSeeds(): input file is not open !\n";
5839 TTree *seedTree=(TTree*)in->Get("Seeds");
5841 cerr<<"AliTPCtrackerMI::ReadSeeds(): ";
5842 cerr<<"can't get a tree with track seeds !\n";
5845 AliTPCtrack *seed=new AliTPCtrack;
5846 seedTree->SetBranchAddress("tracks",&seed);
5848 if (fSeeds==0) fSeeds=new TObjArray(15000);
5850 Int_t n=(Int_t)seedTree->GetEntries();
5851 for (Int_t i=0; i<n; i++) {
5852 seedTree->GetEvent(i);
5853 fSeeds->AddLast(new AliTPCseed(*seed/*,seed->GetAlpha()*/));
5862 Int_t AliTPCtrackerMI::Clusters2Tracks (AliESDEvent *const esd)
5866 if (fSeeds) DeleteSeeds();
5869 if (!fSeeds) return 1;
5871 if (AliTPCReconstructor::StreamLevel()>3) DumpClusters(0,fSeeds);
5877 //_____________________________________________________________________________
5878 Int_t AliTPCtrackerMI::Clusters2Tracks() {
5879 //-----------------------------------------------------------------
5880 // This is a track finder.
5881 //-----------------------------------------------------------------
5882 TDirectory *savedir=gDirectory;
5886 fSeeds = Tracking();
5889 Info("Clusters2Tracks","Time for tracking: \t");timer.Print();timer.Start();
5891 //activate again some tracks
5892 for (Int_t i=0; i<fSeeds->GetEntriesFast(); i++) {
5893 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
5895 Int_t nc=t.GetNumberOfClusters();
5897 delete fSeeds->RemoveAt(i);
5901 if (pt->GetRemoval()==10) {
5902 if (pt->GetDensityFirst(20)>0.8 || pt->GetDensityFirst(30)>0.8 || pt->GetDensityFirst(40)>0.7)
5903 pt->Desactivate(10); // make track again active
5905 pt->Desactivate(20);
5906 delete fSeeds->RemoveAt(i);
5911 RemoveUsed2(fSeeds,0.85,0.85,0);
5912 if (AliTPCReconstructor::GetRecoParam()->GetDoKinks()) FindKinks(fSeeds,fEvent);
5913 //FindCurling(fSeeds, fEvent,0);
5914 if (AliTPCReconstructor::StreamLevel()>5) FindMultiMC(fSeeds, fEvent,-1); // find multi found tracks
5915 RemoveUsed2(fSeeds,0.5,0.4,20);
5916 FindSplitted(fSeeds, fEvent,0); // find multi found tracks
5917 if (AliTPCReconstructor::StreamLevel()>5) FindMultiMC(fSeeds, fEvent,0); // find multi found tracks
5920 // // refit short tracks
5922 Int_t nseed=fSeeds->GetEntriesFast();
5925 for (Int_t i=0; i<nseed; i++) {
5926 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
5928 Int_t nc=t.GetNumberOfClusters();
5930 delete fSeeds->RemoveAt(i);
5933 CookLabel(pt,0.1); //For comparison only
5934 //if ((pt->IsActive() || (pt->fRemoval==10) )&& nc>50 &&pt->GetNumberOfClusters()>0.4*pt->fNFoundable){
5935 if ((pt->IsActive() || (pt->GetRemoval()==10) )){
5937 if (fDebug>0) cerr<<found<<'\r';
5941 delete fSeeds->RemoveAt(i);
5945 //RemoveOverlap(fSeeds,0.99,7,kTRUE);
5947 //RemoveUsed(fSeeds,0.9,0.9,6);
5949 nseed=fSeeds->GetEntriesFast();
5951 for (Int_t i=0; i<nseed; i++) {
5952 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
5954 Int_t nc=t.GetNumberOfClusters();
5956 delete fSeeds->RemoveAt(i);
5960 t.CookdEdx(0.02,0.6);
5961 // CheckKinkPoint(&t,0.05);
5962 //if ((pt->IsActive() || (pt->fRemoval==10) )&& nc>50 &&pt->GetNumberOfClusters()>0.4*pt->fNFoundable){
5963 if ((pt->IsActive() || (pt->GetRemoval()==10) )){
5971 delete fSeeds->RemoveAt(i);
5972 //AliTPCseed * seed1 = ReSeed(pt,0.05,0.5,1);
5974 // FollowProlongation(*seed1,0);
5975 // Int_t n = seed1->GetNumberOfClusters();
5976 // printf("fP4\t%f\t%f\n",seed1->GetC(),pt->GetC());
5977 // printf("fN\t%d\t%d\n", seed1->GetNumberOfClusters(),pt->GetNumberOfClusters());
5980 //AliTPCseed * seed2 = ReSeed(pt,0.95,0.5,0.05);
5984 SortTracks(fSeeds, 1);
5988 PrepareForBackProlongation(fSeeds,5.);
5989 PropagateBack(fSeeds);
5990 printf("Time for back propagation: \t");timer.Print();timer.Start();
5994 PrepareForProlongation(fSeeds,5.);
5995 PropagateForward2(fSeeds);
5997 printf("Time for FORWARD propagation: \t");timer.Print();timer.Start();
5998 // RemoveUsed(fSeeds,0.7,0.7,6);
5999 //RemoveOverlap(fSeeds,0.9,7,kTRUE);
6001 nseed=fSeeds->GetEntriesFast();
6003 for (Int_t i=0; i<nseed; i++) {
6004 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
6006 Int_t nc=t.GetNumberOfClusters();
6008 delete fSeeds->RemoveAt(i);
6011 t.CookdEdx(0.02,0.6);
6012 // CookLabel(pt,0.1); //For comparison only
6013 //if ((pt->IsActive() || (pt->fRemoval==10) )&& nc>50 &&pt->GetNumberOfClusters()>0.4*pt->fNFoundable){
6014 if ((pt->IsActive() || (pt->fRemoval==10) )){
6015 cerr<<found++<<'\r';
6018 delete fSeeds->RemoveAt(i);
6023 // fNTracks = found;
6025 Info("Clusters2Tracks","Time for overlap removal, track writing and dedx cooking: \t"); timer.Print();timer.Start();
6028 // cerr<<"Number of found tracks : "<<"\t"<<found<<endl;
6029 Info("Clusters2Tracks","Number of found tracks %d",found);
6031 // UnloadClusters();
6036 void AliTPCtrackerMI::Tracking(TObjArray * arr)
6039 // tracking of the seeds
6042 fSectors = fOuterSec;
6043 ParallelTracking(arr,150,63);
6044 fSectors = fOuterSec;
6045 ParallelTracking(arr,63,0);
6048 TObjArray * AliTPCtrackerMI::Tracking(Int_t seedtype, Int_t i1, Int_t i2, Float_t cuts[4], Float_t dy, Int_t dsec)
6053 TObjArray * arr = new TObjArray;
6055 fSectors = fOuterSec;
6058 for (Int_t sec=0;sec<fkNOS;sec++){
6059 if (seedtype==3) MakeSeeds3(arr,sec,i1,i2,cuts,dy, dsec);
6060 if (seedtype==4) MakeSeeds5(arr,sec,i1,i2,cuts,dy);
6061 if (seedtype==2) MakeSeeds2(arr,sec,i1,i2,cuts,dy);
6064 Info("Tracking","\nSeeding - %d\t%d\t%d\t%d\n",seedtype,i1,i2,arr->GetEntriesFast());
6076 TObjArray * AliTPCtrackerMI::Tracking()
6080 if (AliTPCReconstructor::GetRecoParam()->GetSpecialSeeding()) return TrackingSpecial();
6083 Int_t nup=fOuterSec->GetNRows()+fInnerSec->GetNRows();
6085 TObjArray * seeds = new TObjArray;
6087 Int_t fLastSeedRowSec=AliTPCReconstructor::GetRecoParam()->GetLastSeedRowSec();
6088 Int_t gapPrim = AliTPCReconstructor::GetRecoParam()->GetSeedGapPrim();
6089 Int_t gapSec = AliTPCReconstructor::GetRecoParam()->GetSeedGapSec();
6097 Float_t fnumber = 3.0;
6098 Float_t fdensity = 3.0;
6103 for (Int_t delta = 0; delta<18; delta+=gapPrim){
6107 arr = Tracking(3,nup-1-delta,nup-1-delta-gap,cuts,-1,1);
6108 SumTracks(seeds,arr);
6109 SignClusters(seeds,fnumber,fdensity);
6111 for (Int_t i=2;i<6;i+=2){
6112 // seed high pt tracks
6115 arr = Tracking(3,nup-i-delta,nup-i-delta-gap,cuts,-1,0);
6116 SumTracks(seeds,arr);
6117 SignClusters(seeds,fnumber,fdensity);
6122 // RemoveUsed(seeds,0.9,0.9,1);
6123 // UnsignClusters();
6124 // SignClusters(seeds,fnumber,fdensity);
6128 for (Int_t delta = 20; delta<120; delta+=gapPrim){
6130 // seed high pt tracks
6134 arr = Tracking(3,nup-delta,nup-delta-gap,cuts,-1);
6135 SumTracks(seeds,arr);
6136 SignClusters(seeds,fnumber,fdensity);
6141 arr = Tracking(3,nup-delta-5,nup-delta-5-gap,cuts,-1);
6142 SumTracks(seeds,arr);
6143 SignClusters(seeds,fnumber,fdensity);
6154 Info("Tracking()","\n\nPrimary seeding\t%d\n\n",seeds->GetEntriesFast());
6158 // RemoveUsed(seeds,0.75,0.75,1);
6160 //SignClusters(seeds,fnumber,fdensity);
6169 arr = Tracking(4,nup-1,nup-1-gap,cuts,-1);
6170 SumTracks(seeds,arr);
6171 SignClusters(seeds,fnumber,fdensity);
6173 arr = Tracking(4,nup-2,nup-2-gap,cuts,-1);
6174 SumTracks(seeds,arr);
6175 SignClusters(seeds,fnumber,fdensity);
6177 arr = Tracking(4,nup-3,nup-3-gap,cuts,-1);
6178 SumTracks(seeds,arr);
6179 SignClusters(seeds,fnumber,fdensity);
6181 arr = Tracking(4,nup-5,nup-5-gap,cuts,-1);
6182 SumTracks(seeds,arr);
6183 SignClusters(seeds,fnumber,fdensity);
6185 arr = Tracking(4,nup-7,nup-7-gap,cuts,-1);
6186 SumTracks(seeds,arr);
6187 SignClusters(seeds,fnumber,fdensity);
6190 arr = Tracking(4,nup-9,nup-9-gap,cuts,-1);
6191 SumTracks(seeds,arr);
6192 SignClusters(seeds,fnumber,fdensity);
6196 for (Int_t delta = 9; delta<30; delta+=gapSec){
6202 arr = Tracking(4,nup-1-delta,nup-1-delta-gap,cuts,-1);
6203 SumTracks(seeds,arr);
6204 SignClusters(seeds,fnumber,fdensity);
6206 arr = Tracking(4,nup-3-delta,nup-5-delta-gap,cuts,4);
6207 SumTracks(seeds,arr);
6208 SignClusters(seeds,fnumber,fdensity);
6221 for (Int_t delta = 30; delta<fLastSeedRowSec; delta+=gapSec){
6227 arr = Tracking(4,nup-1-delta,nup-1-delta-gap,cuts,-1);
6228 SumTracks(seeds,arr);
6229 SignClusters(seeds,fnumber,fdensity);
6231 arr = Tracking(4,nup-5-delta,nup-5-delta-gap,cuts,5 );
6232 SumTracks(seeds,arr);
6233 SignClusters(seeds,fnumber,fdensity);
6237 Info("Tracking()","\n\nSecondary seeding\t%d\n\n",seeds->GetEntriesFast());
6248 TObjArray * AliTPCtrackerMI::TrackingSpecial()
6251 // seeding adjusted for laser and cosmic tests - short tracks with big inclination angle
6252 // no primary vertex seeding tried
6256 Int_t nup=fOuterSec->GetNRows()+fInnerSec->GetNRows();
6258 TObjArray * seeds = new TObjArray;
6263 Float_t fnumber = 3.0;
6264 Float_t fdensity = 3.0;
6267 cuts[0] = AliTPCReconstructor::GetRecoParam()->GetMaxC(); // max curvature
6268 cuts[1] = 3.5; // max tan(phi) angle for seeding
6269 cuts[2] = 3.; // not used (cut on z primary vertex)
6270 cuts[3] = 3.5; // max tan(theta) angle for seeding
6272 for (Int_t delta = 0; nup-delta-gap-1>0; delta+=3){
6274 arr = Tracking(4,nup-1-delta,nup-1-delta-gap,cuts,-1);
6275 SumTracks(seeds,arr);
6276 SignClusters(seeds,fnumber,fdensity);
6280 Info("Tracking()","\n\nSecondary seeding\t%d\n\n",seeds->GetEntriesFast());
6291 void AliTPCtrackerMI::SumTracks(TObjArray *arr1,TObjArray *&arr2) const
6294 //sum tracks to common container
6295 //remove suspicious tracks
6296 Int_t nseed = arr2->GetEntriesFast();
6297 for (Int_t i=0;i<nseed;i++){
6298 AliTPCseed *pt=(AliTPCseed*)arr2->UncheckedAt(i);
6301 // remove tracks with too big curvature
6303 if (TMath::Abs(pt->GetC())>AliTPCReconstructor::GetRecoParam()->GetMaxC()){
6304 delete arr2->RemoveAt(i);
6307 // REMOVE VERY SHORT TRACKS
6308 if (pt->GetNumberOfClusters()<20){
6309 delete arr2->RemoveAt(i);
6312 // NORMAL ACTIVE TRACK
6313 if (pt->IsActive()){
6314 arr1->AddLast(arr2->RemoveAt(i));
6317 //remove not usable tracks
6318 if (pt->GetRemoval()!=10){
6319 delete arr2->RemoveAt(i);
6323 // ENABLE ONLY ENOUGH GOOD STOPPED TRACKS
6324 if (pt->GetDensityFirst(20)>0.8 || pt->GetDensityFirst(30)>0.8 || pt->GetDensityFirst(40)>0.7)
6325 arr1->AddLast(arr2->RemoveAt(i));
6327 delete arr2->RemoveAt(i);
6331 delete arr2; arr2 = 0;
6336 void AliTPCtrackerMI::ParallelTracking(TObjArray *const arr, Int_t rfirst, Int_t rlast)
6339 // try to track in parralel
6341 Int_t nseed=arr->GetEntriesFast();
6342 //prepare seeds for tracking
6343 for (Int_t i=0; i<nseed; i++) {
6344 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i), &t=*pt;
6346 if (!t.IsActive()) continue;
6347 // follow prolongation to the first layer
6348 if ( (fSectors ==fInnerSec) || (t.GetFirstPoint()-fkParam->GetNRowLow()>rfirst+1) )
6349 FollowProlongation(t, rfirst+1);
6354 for (Int_t nr=rfirst; nr>=rlast; nr--){
6355 if (nr<fInnerSec->GetNRows())
6356 fSectors = fInnerSec;
6358 fSectors = fOuterSec;
6359 // make indexes with the cluster tracks for given
6361 // find nearest cluster
6362 for (Int_t i=0; i<nseed; i++) {
6363 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i), &t=*pt;
6365 if (nr==80) pt->UpdateReference();
6366 if (!pt->IsActive()) continue;
6367 // if ( (fSectors ==fOuterSec) && (pt->fFirstPoint-fkParam->GetNRowLow())<nr) continue;
6368 if (pt->GetRelativeSector()>17) {
6371 UpdateClusters(t,nr);
6373 // prolonagate to the nearest cluster - if founded
6374 for (Int_t i=0; i<nseed; i++) {
6375 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
6377 if (!pt->IsActive()) continue;
6378 // if ((fSectors ==fOuterSec) && (pt->fFirstPoint-fkParam->GetNRowLow())<nr) continue;
6379 if (pt->GetRelativeSector()>17) {
6382 FollowToNextCluster(*pt,nr);
6387 void AliTPCtrackerMI::PrepareForBackProlongation(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);
6399 //rotate to current local system at first accepted point
6400 Int_t index = pt->GetClusterIndex2(pt->GetFirstPoint());
6401 Int_t sec = (index&0xff000000)>>24;
6403 Float_t angle1 = fInnerSec->GetAlpha()*sec+fInnerSec->GetAlphaShift();
6404 if (angle1>TMath::Pi())
6405 angle1-=2.*TMath::Pi();
6406 Float_t angle2 = pt->GetAlpha();
6408 if (TMath::Abs(angle1-angle2)>0.001){
6409 if (!pt->Rotate(angle1-angle2)) return;
6410 //angle2 = pt->GetAlpha();
6411 //pt->fRelativeSector = pt->GetAlpha()/fInnerSec->GetAlpha();
6412 //if (pt->GetAlpha()<0)
6413 // pt->fRelativeSector+=18;
6414 //sec = pt->fRelativeSector;
6423 void AliTPCtrackerMI::PrepareForProlongation(TObjArray *const arr, Float_t fac) const
6427 // if we use TPC track itself we have to "update" covariance
6429 Int_t nseed= arr->GetEntriesFast();
6430 for (Int_t i=0;i<nseed;i++){
6431 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
6434 pt->SetFirstPoint(pt->GetLastPoint());
6442 Int_t AliTPCtrackerMI::PropagateBack(TObjArray *const arr)
6445 // make back propagation
6447 Int_t nseed= arr->GetEntriesFast();
6448 for (Int_t i=0;i<nseed;i++){
6449 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
6450 if (pt&& pt->GetKinkIndex(0)<=0) {
6451 //AliTPCseed *pt2 = new AliTPCseed(*pt);
6452 fSectors = fInnerSec;
6453 //FollowBackProlongation(*pt,fInnerSec->GetNRows()-1);
6454 //fSectors = fOuterSec;
6455 FollowBackProlongation(*pt,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1);
6456 //if (pt->GetNumberOfClusters()<(pt->fEsd->GetTPCclusters(0)) ){
6457 // Error("PropagateBack","Not prolonged track %d",pt->GetLabel());
6458 // FollowBackProlongation(*pt2,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1);
6461 if (pt&& pt->GetKinkIndex(0)>0) {
6462 AliESDkink * kink = fEvent->GetKink(pt->GetKinkIndex(0)-1);
6463 pt->SetFirstPoint(kink->GetTPCRow0());
6464 fSectors = fInnerSec;
6465 FollowBackProlongation(*pt,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1);
6473 Int_t AliTPCtrackerMI::PropagateForward2(TObjArray *const arr)
6476 // make forward propagation
6478 Int_t nseed= arr->GetEntriesFast();
6480 for (Int_t i=0;i<nseed;i++){
6481 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
6483 FollowProlongation(*pt,0);
6492 Int_t AliTPCtrackerMI::PropagateForward()
6495 // propagate track forward
6497 Int_t nseed = fSeeds->GetEntriesFast();
6498 for (Int_t i=0;i<nseed;i++){
6499 AliTPCseed *pt = (AliTPCseed*)fSeeds->UncheckedAt(i);
6501 AliTPCseed &t = *pt;
6502 Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
6503 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
6504 if (alpha < 0. ) alpha += 2.*TMath::Pi();
6505 t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN);
6509 fSectors = fOuterSec;
6510 ParallelTracking(fSeeds,fOuterSec->GetNRows()+fInnerSec->GetNRows()-1,fInnerSec->GetNRows());
6511 fSectors = fInnerSec;
6512 ParallelTracking(fSeeds,fInnerSec->GetNRows()-1,0);
6521 Int_t AliTPCtrackerMI::PropagateBack(AliTPCseed *const pt, Int_t row0, Int_t row1)
6524 // make back propagation, in between row0 and row1
6528 fSectors = fInnerSec;
6531 if (row1<fSectors->GetNRows())
6534 r1 = fSectors->GetNRows()-1;
6536 if (row0<fSectors->GetNRows()&& r1>0 )
6537 FollowBackProlongation(*pt,r1);
6538 if (row1<=fSectors->GetNRows())
6541 r1 = row1 - fSectors->GetNRows();
6542 if (r1<=0) return 0;
6543 if (r1>=fOuterSec->GetNRows()) return 0;
6544 fSectors = fOuterSec;
6545 return FollowBackProlongation(*pt,r1);
6553 void AliTPCtrackerMI::GetShape(AliTPCseed * seed, Int_t row)
6557 AliTPCClusterParam * clparam = AliTPCcalibDB::Instance()->GetClusterParam();
6558 Float_t zdrift = TMath::Abs((fkParam->GetZLength(0)-TMath::Abs(seed->GetZ())));
6559 Int_t type = (seed->GetSector() < fkParam->GetNSector()/2) ? 0: (row>126) ? 1:2;
6560 Double_t angulary = seed->GetSnp();
6562 if (TMath::Abs(angulary)>AliTPCReconstructor::GetMaxSnpTracker()) {
6563 angulary = TMath::Sign(AliTPCReconstructor::GetMaxSnpTracker(),angulary);
6566 angulary = angulary*angulary/((1.-angulary)*(1.+angulary));
6567 Double_t angularz = seed->GetTgl()*seed->GetTgl()*(1.+angulary);
6569 Double_t sigmay = clparam->GetRMS0(0,type,zdrift,TMath::Sqrt(TMath::Abs(angulary)));
6570 Double_t sigmaz = clparam->GetRMS0(1,type,zdrift,TMath::Sqrt(TMath::Abs(angularz)));
6571 seed->SetCurrentSigmaY2(sigmay*sigmay);
6572 seed->SetCurrentSigmaZ2(sigmaz*sigmaz);
6573 // Float_t sd2 = TMath::Abs((fkParam->GetZLength(0)-TMath::Abs(seed->GetZ())))*fkParam->GetDiffL()*fkParam->GetDiffL();
6574 // // Float_t padlength = fkParam->GetPadPitchLength(seed->fSector);
6575 // Float_t padlength = GetPadPitchLength(row);
6577 // Float_t sresy = (seed->GetSector() < fkParam->GetNSector()/2) ? 0.2 :0.3;
6578 // seed->SetCurrentSigmaY2(sd2+padlength*padlength*angulary/12.+sresy*sresy);
6580 // Float_t sresz = fkParam->GetZSigma();
6581 // seed->SetCurrentSigmaZ2(sd2+padlength*padlength*angularz*angularz*(1+angulary)/12.+sresz*sresz);
6583 Float_t wy = GetSigmaY(seed);
6584 Float_t wz = GetSigmaZ(seed);
6587 if (TMath::Abs(wy/seed->fCurrentSigmaY2-1)>0.0001 || TMath::Abs(wz/seed->fCurrentSigmaZ2-1)>0.0001 ){
6588 printf("problem\n");
6595 //__________________________________________________________________________
6596 void AliTPCtrackerMI::CookLabel(AliKalmanTrack *tk, Float_t wrong) const {
6597 //--------------------------------------------------------------------
6598 //This function "cooks" a track label. If label<0, this track is fake.
6599 //--------------------------------------------------------------------
6600 AliTPCseed * t = dynamic_cast<AliTPCseed*>(tk);
6602 printf("%s:%d wrong type \n",(char*)__FILE__,__LINE__);
6606 Int_t noc=t->GetNumberOfClusters();
6608 //printf("\nnot founded prolongation\n\n\n");
6614 AliTPCclusterMI *clusters[160];
6616 for (Int_t i=0;i<160;i++) {
6623 for (i=0; i<160 && current<noc; i++) {
6625 Int_t index=t->GetClusterIndex2(i);
6626 if (index<=0) continue;
6627 if (index&0x8000) continue;
6629 //clusters[current]=GetClusterMI(index);
6630 if (t->GetClusterPointer(i)){
6631 clusters[current]=t->GetClusterPointer(i);
6637 Int_t lab=123456789;
6638 for (i=0; i<noc; i++) {
6639 AliTPCclusterMI *c=clusters[i];
6641 lab=TMath::Abs(c->GetLabel(0));
6643 for (j=0; j<noc; j++) if (lb[j]==lab || mx[j]==0) break;
6649 for (i=0; i<noc; i++) if (mx[i]>max) {max=mx[i]; lab=lb[i];}
6651 for (i=0; i<noc; i++) {
6652 AliTPCclusterMI *c=clusters[i];
6654 if (TMath::Abs(c->GetLabel(1)) == lab ||
6655 TMath::Abs(c->GetLabel(2)) == lab ) max++;
6657 if (noc<=0) { lab=-1; return;}
6658 if ((1.- Float_t(max)/(noc)) > wrong) lab=-lab;
6661 Int_t tail=Int_t(0.10*noc);
6664 for (i=1; i<160&&ind<tail; i++) {
6665 // AliTPCclusterMI *c=clusters[noc-i];
6666 AliTPCclusterMI *c=clusters[i];
6668 if (lab == TMath::Abs(c->GetLabel(0)) ||
6669 lab == TMath::Abs(c->GetLabel(1)) ||
6670 lab == TMath::Abs(c->GetLabel(2))) max++;
6673 if (max < Int_t(0.5*tail)) lab=-lab;
6680 //delete[] clusters;
6684 //__________________________________________________________________________
6685 Int_t AliTPCtrackerMI::CookLabel(AliTPCseed *const t, Float_t wrong,Int_t first, Int_t last) const {
6686 //--------------------------------------------------------------------
6687 //This function "cooks" a track label. If label<0, this track is fake.
6688 //--------------------------------------------------------------------
6689 Int_t noc=t->GetNumberOfClusters();
6691 //printf("\nnot founded prolongation\n\n\n");
6697 AliTPCclusterMI *clusters[160];
6699 for (Int_t i=0;i<160;i++) {
6706 for (i=0; i<160 && current<noc; i++) {
6707 if (i<first) continue;
6708 if (i>last) continue;
6709 Int_t index=t->GetClusterIndex2(i);
6710 if (index<=0) continue;
6711 if (index&0x8000) continue;
6713 //clusters[current]=GetClusterMI(index);
6714 if (t->GetClusterPointer(i)){
6715 clusters[current]=t->GetClusterPointer(i);
6720 //if (noc<5) return -1;
6721 Int_t lab=123456789;
6722 for (i=0; i<noc; i++) {
6723 AliTPCclusterMI *c=clusters[i];
6725 lab=TMath::Abs(c->GetLabel(0));
6727 for (j=0; j<noc; j++) if (lb[j]==lab || mx[j]==0) break;
6733 for (i=0; i<noc; i++) if (mx[i]>max) {max=mx[i]; lab=lb[i];}
6735 for (i=0; i<noc; i++) {
6736 AliTPCclusterMI *c=clusters[i];
6738 if (TMath::Abs(c->GetLabel(1)) == lab ||
6739 TMath::Abs(c->GetLabel(2)) == lab ) max++;
6741 if (noc<=0) { lab=-1; return -1;}
6742 if ((1.- Float_t(max)/(noc)) > wrong) lab=-lab;
6745 Int_t tail=Int_t(0.10*noc);
6748 for (i=1; i<160&&ind<tail; i++) {
6749 // AliTPCclusterMI *c=clusters[noc-i];
6750 AliTPCclusterMI *c=clusters[i];
6752 if (lab == TMath::Abs(c->GetLabel(0)) ||
6753 lab == TMath::Abs(c->GetLabel(1)) ||
6754 lab == TMath::Abs(c->GetLabel(2))) max++;
6757 if (max < Int_t(0.5*tail)) lab=-lab;
6760 // t->SetLabel(lab);
6764 //delete[] clusters;
6768 Int_t AliTPCtrackerMI::GetRowNumber(Double_t x[3]) const
6770 //return pad row number for given x vector
6771 Float_t phi = TMath::ATan2(x[1],x[0]);
6772 if(phi<0) phi=2.*TMath::Pi()+phi;
6773 // Get the local angle in the sector philoc
6774 const Float_t kRaddeg = 180/3.14159265358979312;
6775 Float_t phiangle = (Int_t (phi*kRaddeg/20.) + 0.5)*20./kRaddeg;
6776 Double_t localx = x[0]*TMath::Cos(phiangle)-x[1]*TMath::Sin(phiangle);
6777 return GetRowNumber(localx);
6782 void AliTPCtrackerMI::MakeBitmaps(AliTPCseed *t)
6784 //-----------------------------------------------------------------------
6785 // Fill the cluster and sharing bitmaps of the track
6786 //-----------------------------------------------------------------------
6788 Int_t firstpoint = 0;
6789 Int_t lastpoint = 159;
6790 AliTPCTrackerPoint *point;
6791 AliTPCclusterMI *cluster;
6793 for (int iter=firstpoint; iter<lastpoint; iter++) {
6794 // Change to cluster pointers to see if we have a cluster at given padrow
6795 cluster = t->GetClusterPointer(iter);
6797 t->SetClusterMapBit(iter, kTRUE);
6798 point = t->GetTrackPoint(iter);
6799 if (point->IsShared())
6800 t->SetSharedMapBit(iter,kTRUE);
6802 t->SetSharedMapBit(iter, kFALSE);
6805 t->SetClusterMapBit(iter, kFALSE);
6806 t->SetSharedMapBit(iter, kFALSE);
6811 Bool_t AliTPCtrackerMI::IsFindable(AliTPCseed & track){
6813 // return flag if there is findable cluster at given position
6816 Float_t z = track.GetZ();
6818 if (TMath::Abs(z)<(AliTPCReconstructor::GetCtgRange()*track.GetX()+kDeltaZ) &&
6819 TMath::Abs(z)<fkParam->GetZLength(0) &&
6820 (TMath::Abs(track.GetSnp())<AliTPCReconstructor::GetMaxSnpTracker()))
6826 void AliTPCtrackerMI::AddCovariance(AliTPCseed * seed){
6828 // Adding systematic error
6829 // !!!! the systematic error for element 4 is in 1/cm not in pt
6831 const Double_t *param = AliTPCReconstructor::GetRecoParam()->GetSystematicError();
6832 Double_t *covarIn= (Double_t*)seed->GetCovariance();
6834 for (Int_t i=0;i<15;i++) covar[i]=0;
6840 covar[0] = param[0]*param[0];
6841 covar[2] = param[1]*param[1];
6842 covar[5] = param[2]*param[2];
6843 covar[9] = param[3]*param[3];
6844 Double_t facC = AliTracker::GetBz()*kB2C;
6845 covar[14]= param[4]*param[4]*facC*facC;
6847 covar[1]=TMath::Sqrt((covar[0]*covar[2]))*covarIn[1]/TMath::Sqrt((covarIn[0]*covarIn[2]));
6849 covar[3]=TMath::Sqrt((covar[0]*covar[5]))*covarIn[3]/TMath::Sqrt((covarIn[0]*covarIn[5]));
6850 covar[4]=TMath::Sqrt((covar[2]*covar[5]))*covarIn[4]/TMath::Sqrt((covarIn[2]*covarIn[5]));
6852 covar[6]=TMath::Sqrt((covar[0]*covar[9]))*covarIn[6]/TMath::Sqrt((covarIn[0]*covarIn[9]));
6853 covar[7]=TMath::Sqrt((covar[2]*covar[9]))*covarIn[7]/TMath::Sqrt((covarIn[2]*covarIn[9]));
6854 covar[8]=TMath::Sqrt((covar[5]*covar[9]))*covarIn[8]/TMath::Sqrt((covarIn[5]*covarIn[9]));
6856 covar[10]=TMath::Sqrt((covar[0]*covar[14]))*covarIn[10]/TMath::Sqrt((covarIn[0]*covarIn[14]));
6857 covar[11]=TMath::Sqrt((covar[2]*covar[14]))*covarIn[11]/TMath::Sqrt((covarIn[2]*covarIn[14]));
6858 covar[12]=TMath::Sqrt((covar[5]*covar[14]))*covarIn[12]/TMath::Sqrt((covarIn[5]*covarIn[14]));
6859 covar[13]=TMath::Sqrt((covar[9]*covar[14]))*covarIn[13]/TMath::Sqrt((covarIn[9]*covarIn[14]));
6861 seed->AddCovariance(covar);