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); where nis bigger 0
72 // The debug level - different procedure produce tree for numerical debugging
73 // To enable them set AliTPCReconstructor::SetStreamLevel(n); where nis bigger 1
75 // There are several places in the code which can be numerically debuged
76 // This code is keeped in order to enable code development and to check the calibration implementtion
78 // 1. ErrParam stream (Log level 9) - dump information about
80 // 2.a) cluster error estimate
81 // 3.a) cluster shape estimate
84 //-------------------------------------------------------
89 #include "Riostream.h"
90 #include <TClonesArray.h>
92 #include <TObjArray.h>
95 #include "AliComplexCluster.h"
96 #include "AliESDEvent.h"
97 #include "AliESDtrack.h"
98 #include "AliESDVertex.h"
101 #include "AliHelix.h"
102 #include "AliRunLoader.h"
103 #include "AliTPCClustersRow.h"
104 #include "AliTPCParam.h"
105 #include "AliTPCReconstructor.h"
106 #include "AliTPCpolyTrack.h"
107 #include "AliTPCreco.h"
108 #include "AliTPCseed.h"
110 #include "AliTPCtrackerSector.h"
111 #include "AliTPCtrackerMI.h"
112 #include "TStopwatch.h"
113 #include "AliTPCReconstructor.h"
115 #include "TTreeStream.h"
116 #include "AliAlignObj.h"
117 #include "AliTrackPointArray.h"
119 #include "AliTPCcalibDB.h"
120 #include "AliTPCTransform.h"
121 #include "AliTPCClusterParam.h"
125 ClassImp(AliTPCtrackerMI)
129 class AliTPCFastMath {
132 static Double_t FastAsin(Double_t x);
134 static Double_t fgFastAsin[20000]; //lookup table for fast asin computation
137 Double_t AliTPCFastMath::fgFastAsin[20000];
138 AliTPCFastMath gAliTPCFastMath; // needed to fill the LUT
140 AliTPCFastMath::AliTPCFastMath(){
142 // initialized lookup table;
143 for (Int_t i=0;i<10000;i++){
144 fgFastAsin[2*i] = TMath::ASin(i/10000.);
145 fgFastAsin[2*i+1] = (TMath::ASin((i+1)/10000.)-fgFastAsin[2*i]);
149 Double_t AliTPCFastMath::FastAsin(Double_t x){
151 // return asin using lookup table
153 Int_t index = int(x*10000);
154 return fgFastAsin[2*index]+(x*10000.-index)*fgFastAsin[2*index+1];
157 Int_t index = int(x*10000);
158 return -(fgFastAsin[2*index]+(x*10000.-index)*fgFastAsin[2*index+1]);
160 //__________________________________________________________________
161 AliTPCtrackerMI::AliTPCtrackerMI()
183 // default constructor
186 //_____________________________________________________________________
190 Int_t AliTPCtrackerMI::UpdateTrack(AliTPCseed * track, Int_t accept){
192 //update track information using current cluster - track->fCurrentCluster
195 AliTPCclusterMI* c =track->GetCurrentCluster();
196 if (accept>0) track->SetCurrentClusterIndex1(track->GetCurrentClusterIndex1() | 0x8000); //sign not accepted clusters
198 UInt_t i = track->GetCurrentClusterIndex1();
200 Int_t sec=(i&0xff000000)>>24;
201 //Int_t row = (i&0x00ff0000)>>16;
202 track->SetRow((i&0x00ff0000)>>16);
203 track->SetSector(sec);
204 // Int_t index = i&0xFFFF;
205 if (sec>=fParam->GetNInnerSector()) track->SetRow(track->GetRow()+fParam->GetNRowLow());
206 track->SetClusterIndex2(track->GetRow(), i);
207 //track->fFirstPoint = row;
208 //if ( track->fLastPoint<row) track->fLastPoint =row;
209 // if (track->fRow<0 || track->fRow>160) {
210 // printf("problem\n");
212 if (track->GetFirstPoint()>track->GetRow())
213 track->SetFirstPoint(track->GetRow());
214 if (track->GetLastPoint()<track->GetRow())
215 track->SetLastPoint(track->GetRow());
218 track->SetClusterPointer(track->GetRow(),c);
221 Double_t angle2 = track->GetSnp()*track->GetSnp();
223 //SET NEW Track Point
225 if (angle2<1) //PH sometimes angle2 is very big. To be investigated...
227 angle2 = TMath::Sqrt(angle2/(1-angle2));
228 AliTPCTrackerPoint &point =*(track->GetTrackPoint(track->GetRow()));
230 point.SetSigmaY(c->GetSigmaY2()/track->GetCurrentSigmaY2());
231 point.SetSigmaZ(c->GetSigmaZ2()/track->GetCurrentSigmaZ2());
232 point.SetErrY(sqrt(track->GetErrorY2()));
233 point.SetErrZ(sqrt(track->GetErrorZ2()));
235 point.SetX(track->GetX());
236 point.SetY(track->GetY());
237 point.SetZ(track->GetZ());
238 point.SetAngleY(angle2);
239 point.SetAngleZ(track->GetTgl());
240 if (point.IsShared()){
241 track->SetErrorY2(track->GetErrorY2()*4);
242 track->SetErrorZ2(track->GetErrorZ2()*4);
246 Double_t chi2 = track->GetPredictedChi2(track->GetCurrentCluster());
248 // track->SetErrorY2(track->GetErrorY2()*1.3);
249 // track->SetErrorY2(track->GetErrorY2()+0.01);
250 // track->SetErrorZ2(track->GetErrorZ2()*1.3);
251 // track->SetErrorZ2(track->GetErrorZ2()+0.005);
253 if (accept>0) return 0;
254 if (track->GetNumberOfClusters()%20==0){
255 // if (track->fHelixIn){
256 // TClonesArray & larr = *(track->fHelixIn);
257 // Int_t ihelix = larr.GetEntriesFast();
258 // new(larr[ihelix]) AliHelix(*track) ;
261 track->SetNoCluster(0);
262 return track->Update(c,chi2,i);
267 Int_t AliTPCtrackerMI::AcceptCluster(AliTPCseed * seed, AliTPCclusterMI * cluster)
270 // decide according desired precision to accept given
271 // cluster for tracking
272 Double_t sy2=ErrY2(seed,cluster);
273 Double_t sz2=ErrZ2(seed,cluster);
275 Double_t sdistancey2 = sy2+seed->GetSigmaY2();
276 Double_t sdistancez2 = sz2+seed->GetSigmaZ2();
278 Double_t rdistancey2 = (seed->GetCurrentCluster()->GetY()-seed->GetY())*
279 (seed->GetCurrentCluster()->GetY()-seed->GetY())/sdistancey2;
280 Double_t rdistancez2 = (seed->GetCurrentCluster()->GetZ()-seed->GetZ())*
281 (seed->GetCurrentCluster()->GetZ()-seed->GetZ())/sdistancez2;
283 Double_t rdistance2 = rdistancey2+rdistancez2;
286 if (AliTPCReconstructor::StreamLevel()>5) {
287 Float_t rmsy2 = seed->GetCurrentSigmaY2();
288 Float_t rmsz2 = seed->GetCurrentSigmaZ2();
289 Float_t rmsy2p30 = seed->GetCMeanSigmaY2p30();
290 Float_t rmsz2p30 = seed->GetCMeanSigmaZ2p30();
291 Float_t rmsy2p30R = seed->GetCMeanSigmaY2p30R();
292 Float_t rmsz2p30R = seed->GetCMeanSigmaZ2p30R();
294 AliExternalTrackParam param(*seed);
295 (*fDebugStreamer)<<"ErrParam"<<
302 "rmsy2p30="<<rmsy2p30<<
303 "rmsz2p30="<<rmsz2p30<<
304 "rmsy2p30R="<<rmsy2p30R<<
305 "rmsz2p30R="<<rmsz2p30R<<
309 if (rdistance2>16) return 3;
312 if ((rdistancey2>9. || rdistancez2>9.) && cluster->GetType()==0)
313 return 2; //suspisiouce - will be changed
315 if ((rdistancey2>6.25 || rdistancez2>6.25) && cluster->GetType()>0)
316 // strict cut on overlaped cluster
317 return 2; //suspisiouce - will be changed
319 if ( (rdistancey2>1. || rdistancez2>6.25 )
320 && cluster->GetType()<0){
321 seed->SetNFoundable(seed->GetNFoundable()-1);
330 //_____________________________________________________________________________
331 AliTPCtrackerMI::AliTPCtrackerMI(const AliTPCParam *par):
333 fkNIS(par->GetNInnerSector()/2),
335 fkNOS(par->GetNOuterSector()/2),
352 //---------------------------------------------------------------------
353 // The main TPC tracker constructor
354 //---------------------------------------------------------------------
355 fInnerSec=new AliTPCtrackerSector[fkNIS];
356 fOuterSec=new AliTPCtrackerSector[fkNOS];
359 for (i=0; i<fkNIS; i++) fInnerSec[i].Setup(par,0);
360 for (i=0; i<fkNOS; i++) fOuterSec[i].Setup(par,1);
363 Int_t nrowlow = par->GetNRowLow();
364 Int_t nrowup = par->GetNRowUp();
367 for (Int_t i=0;i<nrowlow;i++){
368 fXRow[i] = par->GetPadRowRadiiLow(i);
369 fPadLength[i]= par->GetPadPitchLength(0,i);
370 fYMax[i] = fXRow[i]*TMath::Tan(0.5*par->GetInnerAngle());
374 for (Int_t i=0;i<nrowup;i++){
375 fXRow[i+nrowlow] = par->GetPadRowRadiiUp(i);
376 fPadLength[i+nrowlow] = par->GetPadPitchLength(60,i);
377 fYMax[i+nrowlow] = fXRow[i+nrowlow]*TMath::Tan(0.5*par->GetOuterAngle());
380 fDebugStreamer = new TTreeSRedirector("TPCdebug.root");
382 //________________________________________________________________________
383 AliTPCtrackerMI::AliTPCtrackerMI(const AliTPCtrackerMI &t):
404 //------------------------------------
405 // dummy copy constructor
406 //------------------------------------------------------------------
409 AliTPCtrackerMI & AliTPCtrackerMI::operator=(const AliTPCtrackerMI& /*r*/){
410 //------------------------------
412 //--------------------------------------------------------------
415 //_____________________________________________________________________________
416 AliTPCtrackerMI::~AliTPCtrackerMI() {
417 //------------------------------------------------------------------
418 // TPC tracker destructor
419 //------------------------------------------------------------------
426 if (fDebugStreamer) delete fDebugStreamer;
430 void AliTPCtrackerMI::FillESD(TObjArray* arr)
434 //fill esds using updated tracks
436 // write tracks to the event
437 // store index of the track
438 Int_t nseed=arr->GetEntriesFast();
439 //FindKinks(arr,fEvent);
440 for (Int_t i=0; i<nseed; i++) {
441 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
444 // pt->PropagateTo(fParam->GetInnerRadiusLow());
445 if (pt->GetKinkIndex(0)<=0){ //don't propagate daughter tracks
446 pt->PropagateTo(fParam->GetInnerRadiusLow());
449 if (( pt->GetPoints()[2]- pt->GetPoints()[0])>5 && pt->GetPoints()[3]>0.8){
451 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
452 iotrack.SetTPCPoints(pt->GetPoints());
453 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
454 iotrack.SetV0Indexes(pt->GetV0Indexes());
455 // iotrack.SetTPCpid(pt->fTPCr);
456 //iotrack.SetTPCindex(i);
457 fEvent->AddTrack(&iotrack);
461 if ( (pt->GetNumberOfClusters()>70)&& (Float_t(pt->GetNumberOfClusters())/Float_t(pt->GetNFoundable()))>0.55) {
463 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
464 iotrack.SetTPCPoints(pt->GetPoints());
465 //iotrack.SetTPCindex(i);
466 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
467 iotrack.SetV0Indexes(pt->GetV0Indexes());
468 // iotrack.SetTPCpid(pt->fTPCr);
469 fEvent->AddTrack(&iotrack);
473 // short tracks - maybe decays
475 if ( (pt->GetNumberOfClusters()>30) && (Float_t(pt->GetNumberOfClusters())/Float_t(pt->GetNFoundable()))>0.70) {
476 Int_t found,foundable,shared;
477 pt->GetClusterStatistic(0,60,found, foundable,shared,kFALSE);
478 if ( (found>20) && (pt->GetNShared()/float(pt->GetNumberOfClusters())<0.2)){
480 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
481 //iotrack.SetTPCindex(i);
482 iotrack.SetTPCPoints(pt->GetPoints());
483 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
484 iotrack.SetV0Indexes(pt->GetV0Indexes());
485 //iotrack.SetTPCpid(pt->fTPCr);
486 fEvent->AddTrack(&iotrack);
491 if ( (pt->GetNumberOfClusters()>20) && (Float_t(pt->GetNumberOfClusters())/Float_t(pt->GetNFoundable()))>0.8) {
492 Int_t found,foundable,shared;
493 pt->GetClusterStatistic(0,60,found, foundable,shared,kFALSE);
494 if (found<20) continue;
495 if (pt->GetNShared()/float(pt->GetNumberOfClusters())>0.2) continue;
498 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
499 iotrack.SetTPCPoints(pt->GetPoints());
500 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
501 iotrack.SetV0Indexes(pt->GetV0Indexes());
502 //iotrack.SetTPCpid(pt->fTPCr);
503 //iotrack.SetTPCindex(i);
504 fEvent->AddTrack(&iotrack);
507 // short tracks - secondaties
509 if ( (pt->GetNumberOfClusters()>30) ) {
510 Int_t found,foundable,shared;
511 pt->GetClusterStatistic(128,158,found, foundable,shared,kFALSE);
512 if ( (found>20) && (pt->GetNShared()/float(pt->GetNumberOfClusters())<0.2) &&float(found)/float(foundable)>0.8){
514 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
515 iotrack.SetTPCPoints(pt->GetPoints());
516 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
517 iotrack.SetV0Indexes(pt->GetV0Indexes());
518 //iotrack.SetTPCpid(pt->fTPCr);
519 //iotrack.SetTPCindex(i);
520 fEvent->AddTrack(&iotrack);
525 if ( (pt->GetNumberOfClusters()>15)) {
526 Int_t found,foundable,shared;
527 pt->GetClusterStatistic(138,158,found, foundable,shared,kFALSE);
528 if (found<15) continue;
529 if (foundable<=0) continue;
530 if (pt->GetNShared()/float(pt->GetNumberOfClusters())>0.2) continue;
531 if (float(found)/float(foundable)<0.8) continue;
534 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
535 iotrack.SetTPCPoints(pt->GetPoints());
536 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
537 iotrack.SetV0Indexes(pt->GetV0Indexes());
538 // iotrack.SetTPCpid(pt->fTPCr);
539 //iotrack.SetTPCindex(i);
540 fEvent->AddTrack(&iotrack);
545 printf("Number of filled ESDs-\t%d\n",fEvent->GetNumberOfTracks());
552 Double_t AliTPCtrackerMI::ErrY2(AliTPCseed* seed, AliTPCclusterMI * cl){
555 // Use calibrated cluster error from OCDB
557 AliTPCClusterParam * clparam = AliTPCcalibDB::Instance()->GetClusterParam();
559 Float_t z = TMath::Abs(fParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
560 Int_t ctype = cl->GetType();
561 Int_t type = (cl->GetRow()<63) ? 0: (cl->GetRow()>126) ? 1:2;
562 Double_t angle = seed->GetSnp()*seed->GetSnp();
563 angle = TMath::Sqrt(TMath::Abs(angle/(1.-angle)));
564 Double_t erry2 = clparam->GetError0Par(0,type, z,angle);
566 erry2+=0.5; // edge cluster
569 seed->SetErrorY2(erry2);
573 //calculate look-up table at the beginning
574 // static Bool_t ginit = kFALSE;
575 // static Float_t gnoise1,gnoise2,gnoise3;
576 // static Float_t ggg1[10000];
577 // static Float_t ggg2[10000];
578 // static Float_t ggg3[10000];
579 // static Float_t glandau1[10000];
580 // static Float_t glandau2[10000];
581 // static Float_t glandau3[10000];
583 // static Float_t gcor01[500];
584 // static Float_t gcor02[500];
585 // static Float_t gcorp[500];
589 // if (ginit==kFALSE){
590 // for (Int_t i=1;i<500;i++){
591 // Float_t rsigma = float(i)/100.;
592 // gcor02[i] = TMath::Max(0.78 +TMath::Exp(7.4*(rsigma-1.2)),0.6);
593 // gcor01[i] = TMath::Max(0.72 +TMath::Exp(3.36*(rsigma-1.2)),0.6);
594 // gcorp[i] = TMath::Max(TMath::Power((rsigma+0.5),1.5),1.2);
598 // for (Int_t i=3;i<10000;i++){
602 // Float_t amp = float(i);
603 // Float_t padlength =0.75;
604 // gnoise1 = 0.0004/padlength;
605 // Float_t nel = 0.268*amp;
606 // Float_t nprim = 0.155*amp;
607 // ggg1[i] = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.001*nel/(padlength*padlength))/nel;
608 // glandau1[i] = (2.+0.12*nprim)*0.5* (2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
609 // if (glandau1[i]>1) glandau1[i]=1;
610 // glandau1[i]*=padlength*padlength/12.;
614 // gnoise2 = 0.0004/padlength;
616 // nprim = 0.133*amp;
617 // ggg2[i] = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
618 // glandau2[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
619 // if (glandau2[i]>1) glandau2[i]=1;
620 // glandau2[i]*=padlength*padlength/12.;
625 // gnoise3 = 0.0004/padlength;
627 // nprim = 0.133*amp;
628 // ggg3[i] = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
629 // glandau3[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
630 // if (glandau3[i]>1) glandau3[i]=1;
631 // glandau3[i]*=padlength*padlength/12.;
639 // Int_t amp = int(TMath::Abs(cl->GetQ()));
641 // seed->SetErrorY2(1.);
645 // Float_t z = TMath::Abs(fParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
646 // Int_t ctype = cl->GetType();
647 // Float_t padlength= GetPadPitchLength(seed->GetRow());
648 // Double_t angle2 = seed->GetSnp()*seed->GetSnp();
649 // angle2 = angle2/(1-angle2);
651 // //cluster "quality"
652 // Int_t rsigmay = int(100.*cl->GetSigmaY2()/(seed->GetCurrentSigmaY2()));
655 // if (fSectors==fInnerSec){
656 // snoise2 = gnoise1;
657 // res = ggg1[amp]*z+glandau1[amp]*angle2;
658 // if (ctype==0) res *= gcor01[rsigmay];
661 // res*= gcorp[rsigmay];
665 // if (padlength<1.1){
666 // snoise2 = gnoise2;
667 // res = ggg2[amp]*z+glandau2[amp]*angle2;
668 // if (ctype==0) res *= gcor02[rsigmay];
671 // res*= gcorp[rsigmay];
675 // snoise2 = gnoise3;
676 // res = ggg3[amp]*z+glandau3[amp]*angle2;
677 // if (ctype==0) res *= gcor02[rsigmay];
680 // res*= gcorp[rsigmay];
687 // res*=2.4; // overestimate error 2 times
691 // if (res<2*snoise2)
694 // seed->SetErrorY2(res);
702 Double_t AliTPCtrackerMI::ErrZ2(AliTPCseed* seed, AliTPCclusterMI * cl){
705 // Use calibrated cluster error from OCDB
707 AliTPCClusterParam * clparam = AliTPCcalibDB::Instance()->GetClusterParam();
709 Float_t z = TMath::Abs(fParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
710 Int_t ctype = cl->GetType();
711 Int_t type = (cl->GetRow()<63) ? 0: (cl->GetRow()>126) ? 1:2;
713 Double_t angle2 = seed->GetSnp()*seed->GetSnp();
714 angle2 = seed->GetTgl()*seed->GetTgl()*(1+angle2/(1-angle2));
715 Double_t angle = TMath::Sqrt(TMath::Abs(angle2));
716 Double_t errz2 = clparam->GetError0Par(1,type, z,angle);
718 errz2+=0.5; // edge cluster
721 seed->SetErrorZ2(errz2);
727 // //seed->SetErrorY2(0.1);
729 // //calculate look-up table at the beginning
730 // static Bool_t ginit = kFALSE;
731 // static Float_t gnoise1,gnoise2,gnoise3;
732 // static Float_t ggg1[10000];
733 // static Float_t ggg2[10000];
734 // static Float_t ggg3[10000];
735 // static Float_t glandau1[10000];
736 // static Float_t glandau2[10000];
737 // static Float_t glandau3[10000];
739 // static Float_t gcor01[1000];
740 // static Float_t gcor02[1000];
741 // static Float_t gcorp[1000];
745 // if (ginit==kFALSE){
746 // for (Int_t i=1;i<1000;i++){
747 // Float_t rsigma = float(i)/100.;
748 // gcor02[i] = TMath::Max(0.81 +TMath::Exp(6.8*(rsigma-1.2)),0.6);
749 // gcor01[i] = TMath::Max(0.72 +TMath::Exp(2.04*(rsigma-1.2)),0.6);
750 // gcorp[i] = TMath::Max(TMath::Power((rsigma+0.5),1.5),1.2);
754 // for (Int_t i=3;i<10000;i++){
758 // Float_t amp = float(i);
759 // Float_t padlength =0.75;
760 // gnoise1 = 0.0004/padlength;
761 // Float_t nel = 0.268*amp;
762 // Float_t nprim = 0.155*amp;
763 // ggg1[i] = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.001*nel/(padlength*padlength))/nel;
764 // glandau1[i] = (2.+0.12*nprim)*0.5* (2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
765 // if (glandau1[i]>1) glandau1[i]=1;
766 // glandau1[i]*=padlength*padlength/12.;
770 // gnoise2 = 0.0004/padlength;
772 // nprim = 0.133*amp;
773 // ggg2[i] = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
774 // glandau2[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
775 // if (glandau2[i]>1) glandau2[i]=1;
776 // glandau2[i]*=padlength*padlength/12.;
781 // gnoise3 = 0.0004/padlength;
783 // nprim = 0.133*amp;
784 // ggg3[i] = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
785 // glandau3[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
786 // if (glandau3[i]>1) glandau3[i]=1;
787 // glandau3[i]*=padlength*padlength/12.;
795 // Int_t amp = int(TMath::Abs(cl->GetQ()));
797 // seed->SetErrorY2(1.);
801 // Float_t z = TMath::Abs(fParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
802 // Int_t ctype = cl->GetType();
803 // Float_t padlength= GetPadPitchLength(seed->GetRow());
805 // Double_t angle2 = seed->GetSnp()*seed->GetSnp();
806 // // if (angle2<0.6) angle2 = 0.6;
807 // angle2 = seed->GetTgl()*seed->GetTgl()*(1+angle2/(1-angle2));
809 // //cluster "quality"
810 // Int_t rsigmaz = int(100.*cl->GetSigmaZ2()/(seed->GetCurrentSigmaZ2()));
813 // if (fSectors==fInnerSec){
814 // snoise2 = gnoise1;
815 // res = ggg1[amp]*z+glandau1[amp]*angle2;
816 // if (ctype==0) res *= gcor01[rsigmaz];
819 // res*= gcorp[rsigmaz];
823 // if (padlength<1.1){
824 // snoise2 = gnoise2;
825 // res = ggg2[amp]*z+glandau2[amp]*angle2;
826 // if (ctype==0) res *= gcor02[rsigmaz];
829 // res*= gcorp[rsigmaz];
833 // snoise2 = gnoise3;
834 // res = ggg3[amp]*z+glandau3[amp]*angle2;
835 // if (ctype==0) res *= gcor02[rsigmaz];
838 // res*= gcorp[rsigmaz];
847 // if ((ctype<0) &&<70){
852 // if (res<2*snoise2)
854 // if (res>3) res =3;
855 // seed->SetErrorZ2(res);
863 void AliTPCtrackerMI::RotateToLocal(AliTPCseed *seed)
865 //rotate to track "local coordinata
866 Float_t x = seed->GetX();
867 Float_t y = seed->GetY();
868 Float_t ymax = x*TMath::Tan(0.5*fSectors->GetAlpha());
871 seed->SetRelativeSector((seed->GetRelativeSector()+1) % fN);
872 if (!seed->Rotate(fSectors->GetAlpha()))
874 } else if (y <-ymax) {
875 seed->SetRelativeSector((seed->GetRelativeSector()-1+fN) % fN);
876 if (!seed->Rotate(-fSectors->GetAlpha()))
884 //_____________________________________________________________________________
885 Double_t AliTPCtrackerMI::F1old(Double_t x1,Double_t y1,
886 Double_t x2,Double_t y2,
887 Double_t x3,Double_t y3)
889 //-----------------------------------------------------------------
890 // Initial approximation of the track curvature
891 //-----------------------------------------------------------------
892 Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
893 Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
894 (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
895 Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
896 (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
898 Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
899 if ( xr*xr+yr*yr<=0.00000000000001) return 100;
900 return -xr*yr/sqrt(xr*xr+yr*yr);
905 //_____________________________________________________________________________
906 Double_t AliTPCtrackerMI::F1(Double_t x1,Double_t y1,
907 Double_t x2,Double_t y2,
908 Double_t x3,Double_t y3)
910 //-----------------------------------------------------------------
911 // Initial approximation of the track curvature
912 //-----------------------------------------------------------------
918 Double_t det = x3*y2-x2*y3;
923 Double_t u = 0.5* (x2*(x2-x3)+y2*(y2-y3))/det;
924 Double_t x0 = x3*0.5-y3*u;
925 Double_t y0 = y3*0.5+x3*u;
926 Double_t c2 = 1/TMath::Sqrt(x0*x0+y0*y0);
932 Double_t AliTPCtrackerMI::F2(Double_t x1,Double_t y1,
933 Double_t x2,Double_t y2,
934 Double_t x3,Double_t y3)
936 //-----------------------------------------------------------------
937 // Initial approximation of the track curvature
938 //-----------------------------------------------------------------
944 Double_t det = x3*y2-x2*y3;
949 Double_t u = 0.5* (x2*(x2-x3)+y2*(y2-y3))/det;
950 Double_t x0 = x3*0.5-y3*u;
951 Double_t y0 = y3*0.5+x3*u;
952 Double_t c2 = 1/TMath::Sqrt(x0*x0+y0*y0);
961 //_____________________________________________________________________________
962 Double_t AliTPCtrackerMI::F2old(Double_t x1,Double_t y1,
963 Double_t x2,Double_t y2,
964 Double_t x3,Double_t y3)
966 //-----------------------------------------------------------------
967 // Initial approximation of the track curvature times center of curvature
968 //-----------------------------------------------------------------
969 Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
970 Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
971 (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
972 Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
973 (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
975 Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
977 return -a/(d*y1-b)*xr/sqrt(xr*xr+yr*yr);
980 //_____________________________________________________________________________
981 Double_t AliTPCtrackerMI::F3(Double_t x1,Double_t y1,
982 Double_t x2,Double_t y2,
983 Double_t z1,Double_t z2)
985 //-----------------------------------------------------------------
986 // Initial approximation of the tangent of the track dip angle
987 //-----------------------------------------------------------------
988 return (z1 - z2)/sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
992 Double_t AliTPCtrackerMI::F3n(Double_t x1,Double_t y1,
993 Double_t x2,Double_t y2,
994 Double_t z1,Double_t z2, Double_t c)
996 //-----------------------------------------------------------------
997 // Initial approximation of the tangent of the track dip angle
998 //-----------------------------------------------------------------
1002 //angle1 = (z1-z2)*c/(TMath::ASin(c*x1-ni)-TMath::ASin(c*x2-ni));
1004 Double_t d = TMath::Sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
1005 if (TMath::Abs(d*c*0.5)>1) return 0;
1006 // Double_t angle2 = TMath::ASin(d*c*0.5);
1007 // Double_t angle2 = AliTPCFastMath::FastAsin(d*c*0.5);
1008 Double_t angle2 = (d*c*0.5>0.1)? TMath::ASin(d*c*0.5): AliTPCFastMath::FastAsin(d*c*0.5);
1010 angle2 = (z1-z2)*c/(angle2*2.);
1014 Bool_t AliTPCtrackerMI::GetProlongation(Double_t x1, Double_t x2, Double_t x[5], Double_t &y, Double_t &z)
1015 {//-----------------------------------------------------------------
1016 // This function find proloncation of a track to a reference plane x=x2.
1017 //-----------------------------------------------------------------
1021 if (TMath::Abs(x[4]*x1 - x[2]) >= 0.999) {
1025 Double_t c1=x[4]*x1 - x[2], r1=sqrt(1.- c1*c1);
1026 Double_t c2=x[4]*x2 - x[2], r2=sqrt(1.- c2*c2);
1030 Double_t dy = dx*(c1+c2)/(r1+r2);
1033 Double_t delta = x[4]*dx*(c1+c2)/(c1*r2 + c2*r1);
1035 if (TMath::Abs(delta)>0.01){
1036 dz = x[3]*TMath::ASin(delta)/x[4];
1038 dz = x[3]*AliTPCFastMath::FastAsin(delta)/x[4];
1041 //dz = x[3]*AliTPCFastMath::FastAsin(delta)/x[4];
1049 Int_t AliTPCtrackerMI::LoadClusters (TTree *tree)
1054 return LoadClusters();
1058 Int_t AliTPCtrackerMI::LoadClusters(TObjArray *arr)
1061 // load clusters to the memory
1062 AliTPCClustersRow *clrow = 0x0;
1063 Int_t lower = arr->LowerBound();
1064 Int_t entries = arr->GetEntriesFast();
1066 for (Int_t i=lower; i<entries; i++) {
1067 clrow = (AliTPCClustersRow*) arr->At(i);
1068 if(!clrow) continue;
1069 if(!clrow->GetArray()) continue;
1073 fParam->AdjustSectorRow(clrow->GetID(),sec,row);
1075 for (Int_t icl=0; icl<clrow->GetArray()->GetEntriesFast(); icl++){
1076 Transform((AliTPCclusterMI*)(clrow->GetArray()->At(icl)));
1079 if (clrow->GetArray()->GetEntriesFast()<=0) continue;
1080 AliTPCtrackerRow * tpcrow=0;
1083 tpcrow = &(fInnerSec[sec%fkNIS][row]);
1087 tpcrow = &(fOuterSec[(sec-fkNIS*2)%fkNOS][row]);
1088 left = (sec-fkNIS*2)/fkNOS;
1091 tpcrow->SetN1(clrow->GetArray()->GetEntriesFast());
1092 tpcrow->SetClusters1(new AliTPCclusterMI[tpcrow->GetN1()]);
1093 for (Int_t i=0;i<tpcrow->GetN1();i++)
1094 tpcrow->SetCluster1(i, *(AliTPCclusterMI*)(clrow->GetArray()->At(i)));
1097 tpcrow->SetN2(clrow->GetArray()->GetEntriesFast());
1098 tpcrow->SetClusters2(new AliTPCclusterMI[tpcrow->GetN2()]);
1099 for (Int_t i=0;i<tpcrow->GetN2();i++)
1100 tpcrow->SetCluster2(i,*(AliTPCclusterMI*)(clrow->GetArray()->At(i)));
1111 Int_t AliTPCtrackerMI::LoadClusters()
1114 // load clusters to the memory
1115 AliTPCClustersRow *clrow= new AliTPCClustersRow;
1116 clrow->SetClass("AliTPCclusterMI");
1118 clrow->GetArray()->ExpandCreateFast(10000);
1120 // TTree * tree = fClustersArray.GetTree();
1122 TTree * tree = fInput;
1123 TBranch * br = tree->GetBranch("Segment");
1124 br->SetAddress(&clrow);
1126 Int_t j=Int_t(tree->GetEntries());
1127 for (Int_t i=0; i<j; i++) {
1131 fParam->AdjustSectorRow(clrow->GetID(),sec,row);
1132 for (Int_t icl=0; icl<clrow->GetArray()->GetEntriesFast(); icl++){
1133 Transform((AliTPCclusterMI*)(clrow->GetArray()->At(icl)));
1136 AliTPCtrackerRow * tpcrow=0;
1139 tpcrow = &(fInnerSec[sec%fkNIS][row]);
1143 tpcrow = &(fOuterSec[(sec-fkNIS*2)%fkNOS][row]);
1144 left = (sec-fkNIS*2)/fkNOS;
1147 tpcrow->SetN1(clrow->GetArray()->GetEntriesFast());
1148 tpcrow->SetClusters1(new AliTPCclusterMI[tpcrow->GetN1()]);
1149 for (Int_t i=0;i<tpcrow->GetN1();i++)
1150 tpcrow->SetCluster1(i, *(AliTPCclusterMI*)(clrow->GetArray()->At(i)));
1153 tpcrow->SetN2(clrow->GetArray()->GetEntriesFast());
1154 tpcrow->SetClusters2(new AliTPCclusterMI[tpcrow->GetN2()]);
1155 for (Int_t i=0;i<tpcrow->GetN2();i++)
1156 tpcrow->SetCluster2(i,*(AliTPCclusterMI*)(clrow->GetArray()->At(i)));
1167 void AliTPCtrackerMI::UnloadClusters()
1170 // unload clusters from the memory
1172 Int_t nrows = fOuterSec->GetNRows();
1173 for (Int_t sec = 0;sec<fkNOS;sec++)
1174 for (Int_t row = 0;row<nrows;row++){
1175 AliTPCtrackerRow* tpcrow = &(fOuterSec[sec%fkNOS][row]);
1177 // if (tpcrow->fClusters1) delete []tpcrow->fClusters1;
1178 // if (tpcrow->fClusters2) delete []tpcrow->fClusters2;
1180 tpcrow->ResetClusters();
1183 nrows = fInnerSec->GetNRows();
1184 for (Int_t sec = 0;sec<fkNIS;sec++)
1185 for (Int_t row = 0;row<nrows;row++){
1186 AliTPCtrackerRow* tpcrow = &(fInnerSec[sec%fkNIS][row]);
1188 // if (tpcrow->fClusters1) delete []tpcrow->fClusters1;
1189 //if (tpcrow->fClusters2) delete []tpcrow->fClusters2;
1191 tpcrow->ResetClusters();
1197 void AliTPCtrackerMI::FillClusterArray(TObjArray* array) const{
1199 // Filling cluster to the array - For visualization purposes
1202 nrows = fOuterSec->GetNRows();
1203 for (Int_t sec = 0;sec<fkNOS;sec++)
1204 for (Int_t row = 0;row<nrows;row++){
1205 AliTPCtrackerRow* tpcrow = &(fOuterSec[sec%fkNOS][row]);
1206 if (!tpcrow) continue;
1207 for (Int_t icl = 0;icl<tpcrow->GetN();icl++){
1208 array->AddLast((TObject*)((*tpcrow)[icl]));
1211 nrows = fInnerSec->GetNRows();
1212 for (Int_t sec = 0;sec<fkNIS;sec++)
1213 for (Int_t row = 0;row<nrows;row++){
1214 AliTPCtrackerRow* tpcrow = &(fInnerSec[sec%fkNIS][row]);
1215 if (!tpcrow) continue;
1216 for (Int_t icl = 0;icl<tpcrow->GetN();icl++){
1217 array->AddLast((TObject*)(*tpcrow)[icl]);
1223 void AliTPCtrackerMI::Transform(AliTPCclusterMI * cluster){
1227 AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ;
1229 AliFatal("Tranformations not in calibDB");
1231 Double_t x[3]={cluster->GetRow(),cluster->GetPad(),cluster->GetTimeBin()};
1232 Int_t i[1]={cluster->GetDetector()};
1233 transform->Transform(x,i,0,1);
1234 if (!AliTPCReconstructor::GetRecoParam()->GetBYMirror()){
1235 if (cluster->GetDetector()%36>17){
1241 // in debug mode check the transformation
1243 if (AliTPCReconstructor::StreamLevel()>1) {
1245 cluster->GetGlobalXYZ(gx);
1246 Int_t event = (fEvent==NULL)? 0: fEvent->GetEventNumberInFile();
1247 TTreeSRedirector &cstream = *fDebugStreamer;
1248 cstream<<"Transform"<<
1259 cluster->SetX(x[0]);
1260 cluster->SetY(x[1]);
1261 cluster->SetZ(x[2]);
1267 //if (!fParam->IsGeoRead()) fParam->ReadGeoMatrices();
1268 TGeoHMatrix *mat = fParam->GetClusterMatrix(cluster->GetDetector());
1270 Double_t pos[3]= {cluster->GetX(),cluster->GetY(),cluster->GetZ()};
1271 Double_t posC[3]={cluster->GetX(),cluster->GetY(),cluster->GetZ()};
1272 if (mat) mat->LocalToMaster(pos,posC);
1274 // chack Loading of Geo matrices from GeoManager - TEMPORARY FIX
1276 cluster->SetX(posC[0]);
1277 cluster->SetY(posC[1]);
1278 cluster->SetZ(posC[2]);
1281 //_____________________________________________________________________________
1282 Int_t AliTPCtrackerMI::LoadOuterSectors() {
1283 //-----------------------------------------------------------------
1284 // This function fills outer TPC sectors with clusters.
1285 //-----------------------------------------------------------------
1286 Int_t nrows = fOuterSec->GetNRows();
1288 for (Int_t sec = 0;sec<fkNOS;sec++)
1289 for (Int_t row = 0;row<nrows;row++){
1290 AliTPCtrackerRow* tpcrow = &(fOuterSec[sec%fkNOS][row]);
1291 Int_t sec2 = sec+2*fkNIS;
1293 Int_t ncl = tpcrow->GetN1();
1295 AliTPCclusterMI *c= (tpcrow->GetCluster1(ncl));
1296 index=(((sec2<<8)+row)<<16)+ncl;
1297 tpcrow->InsertCluster(c,index);
1300 ncl = tpcrow->GetN2();
1302 AliTPCclusterMI *c= (tpcrow->GetCluster2(ncl));
1303 index=((((sec2+fkNOS)<<8)+row)<<16)+ncl;
1304 tpcrow->InsertCluster(c,index);
1307 // write indexes for fast acces
1309 for (Int_t i=0;i<510;i++)
1310 tpcrow->SetFastCluster(i,-1);
1311 for (Int_t i=0;i<tpcrow->GetN();i++){
1312 Int_t zi = Int_t((*tpcrow)[i]->GetZ()+255.);
1313 tpcrow->SetFastCluster(zi,i); // write index
1316 for (Int_t i=0;i<510;i++){
1317 if (tpcrow->GetFastCluster(i)<0)
1318 tpcrow->SetFastCluster(i,last);
1320 last = tpcrow->GetFastCluster(i);
1329 //_____________________________________________________________________________
1330 Int_t AliTPCtrackerMI::LoadInnerSectors() {
1331 //-----------------------------------------------------------------
1332 // This function fills inner TPC sectors with clusters.
1333 //-----------------------------------------------------------------
1334 Int_t nrows = fInnerSec->GetNRows();
1336 for (Int_t sec = 0;sec<fkNIS;sec++)
1337 for (Int_t row = 0;row<nrows;row++){
1338 AliTPCtrackerRow* tpcrow = &(fInnerSec[sec%fkNIS][row]);
1341 Int_t ncl = tpcrow->GetN1();
1343 AliTPCclusterMI *c= (tpcrow->GetCluster1(ncl));
1344 index=(((sec<<8)+row)<<16)+ncl;
1345 tpcrow->InsertCluster(c,index);
1348 ncl = tpcrow->GetN2();
1350 AliTPCclusterMI *c= (tpcrow->GetCluster2(ncl));
1351 index=((((sec+fkNIS)<<8)+row)<<16)+ncl;
1352 tpcrow->InsertCluster(c,index);
1355 // write indexes for fast acces
1357 for (Int_t i=0;i<510;i++)
1358 tpcrow->SetFastCluster(i,-1);
1359 for (Int_t i=0;i<tpcrow->GetN();i++){
1360 Int_t zi = Int_t((*tpcrow)[i]->GetZ()+255.);
1361 tpcrow->SetFastCluster(zi,i); // write index
1364 for (Int_t i=0;i<510;i++){
1365 if (tpcrow->GetFastCluster(i)<0)
1366 tpcrow->SetFastCluster(i,last);
1368 last = tpcrow->GetFastCluster(i);
1380 //_________________________________________________________________________
1381 AliTPCclusterMI *AliTPCtrackerMI::GetClusterMI(Int_t index) const {
1382 //--------------------------------------------------------------------
1383 // Return pointer to a given cluster
1384 //--------------------------------------------------------------------
1385 if (index<0) return 0; // no cluster
1386 Int_t sec=(index&0xff000000)>>24;
1387 Int_t row=(index&0x00ff0000)>>16;
1388 Int_t ncl=(index&0x00007fff)>>00;
1390 const AliTPCtrackerRow * tpcrow=0;
1391 AliTPCclusterMI * clrow =0;
1393 if (sec<0 || sec>=fkNIS*4) {
1394 AliWarning(Form("Wrong sector %d",sec));
1399 tpcrow = &(fInnerSec[sec%fkNIS][row]);
1400 if (tpcrow==0) return 0;
1403 if (tpcrow->GetN1()<=ncl) return 0;
1404 clrow = tpcrow->GetClusters1();
1407 if (tpcrow->GetN2()<=ncl) return 0;
1408 clrow = tpcrow->GetClusters2();
1412 tpcrow = &(fOuterSec[(sec-fkNIS*2)%fkNOS][row]);
1413 if (tpcrow==0) return 0;
1415 if (sec-2*fkNIS<fkNOS) {
1416 if (tpcrow->GetN1()<=ncl) return 0;
1417 clrow = tpcrow->GetClusters1();
1420 if (tpcrow->GetN2()<=ncl) return 0;
1421 clrow = tpcrow->GetClusters2();
1425 return &(clrow[ncl]);
1431 Int_t AliTPCtrackerMI::FollowToNext(AliTPCseed& t, Int_t nr) {
1432 //-----------------------------------------------------------------
1433 // This function tries to find a track prolongation to next pad row
1434 //-----------------------------------------------------------------
1436 Double_t x= GetXrow(nr), ymax=GetMaxY(nr);
1437 AliTPCclusterMI *cl=0;
1438 Int_t tpcindex= t.GetClusterIndex2(nr);
1440 // update current shape info every 5 pad-row
1441 // if ( (nr%5==0) || t.GetNumberOfClusters()<2 || (t.fCurrentSigmaY2<0.0001) ){
1445 if (fIteration>0 && tpcindex>=-1){ //if we have already clusters
1447 if (tpcindex==-1) return 0; //track in dead zone
1449 cl = t.GetClusterPointer(nr);
1450 if ( (cl==0) ) cl = GetClusterMI(tpcindex);
1451 t.SetCurrentClusterIndex1(tpcindex);
1454 Int_t relativesector = ((tpcindex&0xff000000)>>24)%18; // if previously accepted cluster in different sector
1455 Float_t angle = relativesector*fSectors->GetAlpha()+fSectors->GetAlphaShift();
1457 if (angle<-TMath::Pi()) angle += 2*TMath::Pi();
1458 if (angle>=TMath::Pi()) angle -= 2*TMath::Pi();
1460 if (TMath::Abs(angle-t.GetAlpha())>0.001){
1461 Double_t rotation = angle-t.GetAlpha();
1462 t.SetRelativeSector(relativesector);
1463 if (!t.Rotate(rotation)) return 0;
1465 if (!t.PropagateTo(x)) return 0;
1467 t.SetCurrentCluster(cl);
1469 Int_t accept = AcceptCluster(&t,t.GetCurrentCluster());
1470 if ((tpcindex&0x8000)==0) accept =0;
1472 //if founded cluster is acceptible
1473 if (cl->IsUsed(11)) { // id cluster is shared inrease uncertainty
1474 t.SetErrorY2(t.GetErrorY2()+0.03);
1475 t.SetErrorZ2(t.GetErrorZ2()+0.03);
1476 t.SetErrorY2(t.GetErrorY2()*3);
1477 t.SetErrorZ2(t.GetErrorZ2()*3);
1479 t.SetNFoundable(t.GetNFoundable()+1);
1480 UpdateTrack(&t,accept);
1485 if (TMath::Abs(t.GetSnp())>AliTPCReconstructor::GetMaxSnpTracker()) return 0; // cut on angle
1487 // not look for new cluster during refitting
1488 t.SetNFoundable(t.GetNFoundable()+1);
1493 // if (TMath::Abs(t.GetSnp())>0.95 || TMath::Abs(x*t.GetC()-t.GetEta())>0.95) return 0;// patch 28 fev 06
1494 Double_t y=t.GetYat(x);
1495 if (TMath::Abs(y)>ymax){
1497 t.SetRelativeSector((t.GetRelativeSector()+1) % fN);
1498 if (!t.Rotate(fSectors->GetAlpha()))
1500 } else if (y <-ymax) {
1501 t.SetRelativeSector((t.GetRelativeSector()-1+fN) % fN);
1502 if (!t.Rotate(-fSectors->GetAlpha()))
1508 if (!t.PropagateTo(x)) {
1509 if (fIteration==0) t.SetRemoval(10);
1513 Double_t z=t.GetZ();
1516 if (!IsActive(t.GetRelativeSector(),nr)) {
1518 t.SetClusterIndex2(nr,-1);
1521 //AliInfo(Form("A - Sector%d phi %f - alpha %f", t.fRelativeSector,y/x, t.GetAlpha()));
1522 Bool_t isActive = IsActive(t.GetRelativeSector(),nr);
1523 Bool_t isActive2 = (nr>=fInnerSec->GetNRows()) ? fOuterSec[t.GetRelativeSector()][nr-fInnerSec->GetNRows()].GetN()>0:fInnerSec[t.GetRelativeSector()][nr].GetN()>0;
1525 if (!isActive || !isActive2) return 0;
1527 const AliTPCtrackerRow &krow=GetRow(t.GetRelativeSector(),nr);
1528 if ( (t.GetSigmaY2()<0) || t.GetSigmaZ2()<0) return 0;
1530 Double_t roadz = 1.;
1532 if (TMath::Abs(TMath::Abs(y)-ymax)<krow.GetDeadZone()){
1534 t.SetClusterIndex2(nr,-1);
1539 if (TMath::Abs(z)<(AliTPCReconstructor::GetCtgRange()*x+10) && TMath::Abs(z)<fParam->GetZLength(0) && (TMath::Abs(t.GetSnp())<AliTPCReconstructor::GetMaxSnpTracker()))
1540 t.SetNFoundable(t.GetNFoundable()+1);
1546 // cl = krow.FindNearest2(y+10.,z,roady,roadz,index);
1547 cl = krow.FindNearest2(y,z,roady,roadz,index);
1548 if (cl) t.SetCurrentClusterIndex1(krow.GetIndex(index));
1551 t.SetCurrentCluster(cl);
1553 if (fIteration==2&&cl->IsUsed(10)) return 0;
1554 Int_t accept = AcceptCluster(&t,t.GetCurrentCluster());
1555 if (fIteration==2&&cl->IsUsed(11)) {
1556 t.SetErrorY2(t.GetErrorY2()+0.03);
1557 t.SetErrorZ2(t.GetErrorZ2()+0.03);
1558 t.SetErrorY2(t.GetErrorY2()*3);
1559 t.SetErrorZ2(t.GetErrorZ2()*3);
1562 if (t.fCurrentCluster->IsUsed(10)){
1567 if (t.fNShared>0.7*t.GetNumberOfClusters()) {
1573 if (accept<3) UpdateTrack(&t,accept);
1576 if ( fIteration==0 && t.GetNFoundable()*0.5 > t.GetNumberOfClusters()) t.SetRemoval(10);
1582 Int_t AliTPCtrackerMI::FollowToNextFast(AliTPCseed& t, Int_t nr) {
1583 //-----------------------------------------------------------------
1584 // This function tries to find a track prolongation to next pad row
1585 //-----------------------------------------------------------------
1587 Double_t x= GetXrow(nr), ymax=GetMaxY(nr);
1589 if (!t.GetProlongation(x,y,z)) {
1595 if (TMath::Abs(y)>ymax){
1598 t.SetRelativeSector((t.GetRelativeSector()+1) % fN);
1599 if (!t.Rotate(fSectors->GetAlpha()))
1601 } else if (y <-ymax) {
1602 t.SetRelativeSector((t.GetRelativeSector()-1+fN) % fN);
1603 if (!t.Rotate(-fSectors->GetAlpha()))
1606 if (!t.PropagateTo(x)) {
1609 t.GetProlongation(x,y,z);
1612 // update current shape info every 2 pad-row
1613 if ( (nr%2==0) || t.GetNumberOfClusters()<2 || (t.GetCurrentSigmaY2()<0.0001) ){
1614 // t.fCurrentSigmaY = GetSigmaY(&t);
1615 //t.fCurrentSigmaZ = GetSigmaZ(&t);
1619 AliTPCclusterMI *cl=0;
1624 const AliTPCtrackerRow &krow=GetRow(t.GetRelativeSector(),nr);
1625 if ( (t.GetSigmaY2()<0) || t.GetSigmaZ2()<0) return 0;
1627 Double_t roadz = 1.;
1630 if (TMath::Abs(TMath::Abs(y)-ymax)<krow.GetDeadZone()){
1632 t.SetClusterIndex2(row,-1);
1637 if (TMath::Abs(z)>(AliTPCReconstructor::GetCtgRange()*x+10)) t.SetClusterIndex2(row,-1);
1641 if ((cl==0)&&(krow)) {
1642 // cl = krow.FindNearest2(y+10,z,roady,roadz,index);
1643 cl = krow.FindNearest2(y,z,roady,roadz,index);
1645 if (cl) t.SetCurrentClusterIndex1(krow.GetIndex(index));
1649 t.SetCurrentCluster(cl);
1650 // Int_t accept = AcceptCluster(&t,t.fCurrentCluster);
1652 t.SetClusterIndex2(row,index);
1653 t.SetClusterPointer(row, cl);
1661 //_________________________________________________________________________
1662 Bool_t AliTPCtrackerMI::GetTrackPoint(Int_t index, AliTrackPoint &p ) const
1664 // Get track space point by index
1665 // return false in case the cluster doesn't exist
1666 AliTPCclusterMI *cl = GetClusterMI(index);
1667 if (!cl) return kFALSE;
1668 Int_t sector = (index&0xff000000)>>24;
1669 // Int_t row = (index&0x00ff0000)>>16;
1671 // xyz[0] = fParam->GetPadRowRadii(sector,row);
1672 xyz[0] = cl->GetX();
1673 xyz[1] = cl->GetY();
1674 xyz[2] = cl->GetZ();
1676 fParam->AdjustCosSin(sector,cos,sin);
1677 Float_t x = cos*xyz[0]-sin*xyz[1];
1678 Float_t y = cos*xyz[1]+sin*xyz[0];
1680 Float_t sigmaY2 = 0.027*cl->GetSigmaY2();
1681 if (sector < fParam->GetNInnerSector()) sigmaY2 *= 2.07;
1682 Float_t sigmaZ2 = 0.066*cl->GetSigmaZ2();
1683 if (sector < fParam->GetNInnerSector()) sigmaZ2 *= 1.77;
1684 cov[0] = sin*sin*sigmaY2;
1685 cov[1] = -sin*cos*sigmaY2;
1687 cov[3] = cos*cos*sigmaY2;
1690 p.SetXYZ(x,y,xyz[2],cov);
1691 AliGeomManager::ELayerID iLayer;
1693 if (sector < fParam->GetNInnerSector()) {
1694 iLayer = AliGeomManager::kTPC1;
1698 iLayer = AliGeomManager::kTPC2;
1699 idet = sector - fParam->GetNInnerSector();
1701 UShort_t volid = AliGeomManager::LayerToVolUID(iLayer,idet);
1702 p.SetVolumeID(volid);
1708 Int_t AliTPCtrackerMI::UpdateClusters(AliTPCseed& t, Int_t nr) {
1709 //-----------------------------------------------------------------
1710 // This function tries to find a track prolongation to next pad row
1711 //-----------------------------------------------------------------
1712 t.SetCurrentCluster(0);
1713 t.SetCurrentClusterIndex1(0);
1715 Double_t xt=t.GetX();
1716 Int_t row = GetRowNumber(xt)-1;
1717 Double_t ymax= GetMaxY(nr);
1719 if (row < nr) return 1; // don't prolongate if not information until now -
1720 // if (TMath::Abs(t.GetSnp())>0.9 && t.GetNumberOfClusters()>40. && fIteration!=2) {
1722 // return 0; // not prolongate strongly inclined tracks
1724 // if (TMath::Abs(t.GetSnp())>0.95) {
1726 // return 0; // not prolongate strongly inclined tracks
1727 // }// patch 28 fev 06
1729 Double_t x= GetXrow(nr);
1731 //t.PropagateTo(x+0.02);
1732 //t.PropagateTo(x+0.01);
1733 if (!t.PropagateTo(x)){
1740 if (TMath::Abs(y)>ymax){
1742 t.SetRelativeSector((t.GetRelativeSector()+1) % fN);
1743 if (!t.Rotate(fSectors->GetAlpha()))
1745 } else if (y <-ymax) {
1746 t.SetRelativeSector((t.GetRelativeSector()-1+fN) % fN);
1747 if (!t.Rotate(-fSectors->GetAlpha()))
1750 // if (!t.PropagateTo(x)){
1757 if (TMath::Abs(t.GetSnp())>AliTPCReconstructor::GetMaxSnpTracker()) return 0;
1759 if (!IsActive(t.GetRelativeSector(),nr)) {
1761 t.SetClusterIndex2(nr,-1);
1764 //AliInfo(Form("A - Sector%d phi %f - alpha %f", t.fRelativeSector,y/x, t.GetAlpha()));
1766 AliTPCtrackerRow &krow=GetRow(t.GetRelativeSector(),nr);
1768 if (TMath::Abs(TMath::Abs(y)-ymax)<krow.GetDeadZone()){
1770 t.SetClusterIndex2(nr,-1);
1775 if (TMath::Abs(t.GetZ())<(AliTPCReconstructor::GetCtgRange()*t.GetX()+10) && (TMath::Abs(t.GetSnp())<AliTPCReconstructor::GetMaxSnpTracker())) t.SetNFoundable(t.GetNFoundable()+1);
1781 if ( (nr%2==0) || t.GetNumberOfClusters()<2){
1782 // t.fCurrentSigmaY = GetSigmaY(&t);
1783 //t.fCurrentSigmaZ = GetSigmaZ(&t);
1787 AliTPCclusterMI *cl=0;
1790 Double_t roady = 1.;
1791 Double_t roadz = 1.;
1795 index = t.GetClusterIndex2(nr);
1796 if ( (index>0) && (index&0x8000)==0){
1797 cl = t.GetClusterPointer(nr);
1798 if ( (cl==0) && (index>0)) cl = GetClusterMI(index);
1799 t.SetCurrentClusterIndex1(index);
1801 t.SetCurrentCluster(cl);
1807 // if (index<0) return 0;
1808 UInt_t uindex = TMath::Abs(index);
1811 //cl = krow.FindNearest2(y+10,z,roady,roadz,uindex);
1812 cl = krow.FindNearest2(y,z,roady,roadz,uindex);
1815 if (cl) t.SetCurrentClusterIndex1(krow.GetIndex(uindex));
1816 t.SetCurrentCluster(cl);
1822 Int_t AliTPCtrackerMI::FollowToNextCluster(AliTPCseed & t, Int_t nr) {
1823 //-----------------------------------------------------------------
1824 // This function tries to find a track prolongation to next pad row
1825 //-----------------------------------------------------------------
1827 //update error according neighborhoud
1829 if (t.GetCurrentCluster()) {
1831 Int_t accept = AcceptCluster(&t,t.GetCurrentCluster());
1833 if (t.GetCurrentCluster()->IsUsed(10)){
1838 t.SetNShared(t.GetNShared()+1);
1839 if (t.GetNShared()>0.7*t.GetNumberOfClusters()) {
1844 if (fIteration>0) accept = 0;
1845 if (accept<3) UpdateTrack(&t,accept);
1849 if ( ( (t.GetSigmaY2()+t.GetSigmaZ2())>0.16)&& t.GetNumberOfClusters()>18) t.SetRemoval(10);
1850 if ( t.GetChi2()/t.GetNumberOfClusters()>6 &&t.GetNumberOfClusters()>18) t.SetRemoval(10);
1852 if (( (t.GetNFoundable()*0.5 > t.GetNumberOfClusters()) || t.GetNoCluster()>15)) t.SetRemoval(10);
1860 //_____________________________________________________________________________
1861 Int_t AliTPCtrackerMI::FollowProlongation(AliTPCseed& t, Int_t rf, Int_t step) {
1862 //-----------------------------------------------------------------
1863 // This function tries to find a track prolongation.
1864 //-----------------------------------------------------------------
1865 Double_t xt=t.GetX();
1867 Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
1868 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
1869 if (alpha < 0. ) alpha += 2.*TMath::Pi();
1871 t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN);
1873 Int_t first = GetRowNumber(xt)-1;
1874 for (Int_t nr= first; nr>=rf; nr-=step) {
1876 if (t.GetKinkIndexes()[0]>0){
1877 for (Int_t i=0;i<3;i++){
1878 Int_t index = t.GetKinkIndexes()[i];
1879 if (index==0) break;
1880 if (index<0) continue;
1882 AliKink * kink = (AliKink*)fEvent->GetKink(index-1);
1884 printf("PROBLEM\n");
1887 Int_t kinkrow = kink->GetTPCRow0()+2+Int_t(0.5/(0.05+kink->GetAngle(2)));
1889 AliExternalTrackParam paramd(t);
1890 kink->SetDaughter(paramd);
1891 kink->SetStatus(2,5);
1898 if (nr==80) t.UpdateReference();
1899 if (nr<fInnerSec->GetNRows())
1900 fSectors = fInnerSec;
1902 fSectors = fOuterSec;
1903 if (FollowToNext(t,nr)==0)
1912 //_____________________________________________________________________________
1913 Int_t AliTPCtrackerMI::FollowProlongationFast(AliTPCseed& t, Int_t rf, Int_t step) {
1914 //-----------------------------------------------------------------
1915 // This function tries to find a track prolongation.
1916 //-----------------------------------------------------------------
1917 Double_t xt=t.GetX();
1919 Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
1920 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
1921 if (alpha < 0. ) alpha += 2.*TMath::Pi();
1922 t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN);
1924 for (Int_t nr=GetRowNumber(xt)-1; nr>=rf; nr-=step) {
1926 if (FollowToNextFast(t,nr)==0)
1927 if (!t.IsActive()) return 0;
1937 Int_t AliTPCtrackerMI::FollowBackProlongation(AliTPCseed& t, Int_t rf) {
1938 //-----------------------------------------------------------------
1939 // This function tries to find a track prolongation.
1940 //-----------------------------------------------------------------
1942 Double_t xt=t.GetX();
1943 Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
1944 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
1945 if (alpha < 0. ) alpha += 2.*TMath::Pi();
1946 t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN);
1948 Int_t first = t.GetFirstPoint();
1949 if (first<GetRowNumber(xt)+1) first = GetRowNumber(xt)+1;
1951 if (first<0) first=0;
1952 for (Int_t nr=first; nr<=rf; nr++) {
1953 // if ( (TMath::Abs(t.GetSnp())>0.95)) break;//patch 28 fev 06
1954 if (t.GetKinkIndexes()[0]<0){
1955 for (Int_t i=0;i<3;i++){
1956 Int_t index = t.GetKinkIndexes()[i];
1957 if (index==0) break;
1958 if (index>0) continue;
1959 index = TMath::Abs(index);
1960 AliKink * kink = (AliKink*)fEvent->GetKink(index-1);
1962 printf("PROBLEM\n");
1965 Int_t kinkrow = kink->GetTPCRow0()-2-Int_t(0.5/(0.05+kink->GetAngle(2)));
1967 AliExternalTrackParam paramm(t);
1968 kink->SetMother(paramm);
1969 kink->SetStatus(2,1);
1976 if (nr<fInnerSec->GetNRows())
1977 fSectors = fInnerSec;
1979 fSectors = fOuterSec;
1990 Float_t AliTPCtrackerMI::OverlapFactor(AliTPCseed * s1, AliTPCseed * s2, Int_t &sum1, Int_t & sum2)
1998 Float_t dz2 =(s1->GetZ() - s2->GetZ());
2001 Float_t dy2 =TMath::Abs((s1->GetY() - s2->GetY()));
2003 Float_t distance = TMath::Sqrt(dz2+dy2);
2004 if (distance>4.) return 0; // if there are far away - not overlap - to reduce combinatorics
2007 Int_t firstpoint = TMath::Min(s1->GetFirstPoint(),s2->GetFirstPoint());
2008 Int_t lastpoint = TMath::Max(s1->GetLastPoint(),s2->GetLastPoint());
2013 if (firstpoint>lastpoint) {
2014 firstpoint =lastpoint;
2019 for (Int_t i=firstpoint-1;i<lastpoint+1;i++){
2020 if (s1->GetClusterIndex2(i)>0) sum1++;
2021 if (s2->GetClusterIndex2(i)>0) sum2++;
2022 if (s1->GetClusterIndex2(i)==s2->GetClusterIndex2(i) && s1->GetClusterIndex2(i)>0) {
2026 if (sum<5) return 0;
2028 Float_t summin = TMath::Min(sum1+1,sum2+1);
2029 Float_t ratio = (sum+1)/Float_t(summin);
2033 void AliTPCtrackerMI::SignShared(AliTPCseed * s1, AliTPCseed * s2)
2037 Float_t thetaCut = 0.2;//+10.*TMath::Sqrt(s1->GetSigmaTglZ()+ s2->GetSigmaTglZ());
2038 if (TMath::Abs(s1->GetTgl()-s2->GetTgl())>thetaCut) return;
2039 Float_t minCl = TMath::Min(s1->GetNumberOfClusters(),s2->GetNumberOfClusters());
2040 Int_t cutN0 = TMath::Max(5,TMath::Nint(0.1*minCl));
2045 //Int_t firstpoint = TMath::Max(s1->GetFirstPoint(),s2->GetFirstPoint());
2046 //Int_t lastpoint = TMath::Min(s1->GetLastPoint(),s2->GetLastPoint());
2047 Int_t firstpoint = 0;
2048 Int_t lastpoint = 160;
2050 // if (firstpoint>=lastpoint-5) return;;
2052 for (Int_t i=firstpoint;i<lastpoint;i++){
2053 // if ( (s1->GetClusterIndex2(i)&0xFFFF8FFF)==(s2->GetClusterIndex2(i)&0xFFFF8FFF) && s1->GetClusterIndex2(i)>0) {
2054 if ( (s1->GetClusterIndex2(i))==(s2->GetClusterIndex2(i)) && s1->GetClusterIndex2(i)>0) {
2056 s1->SetSharedMapBit(i, kTRUE);
2057 s2->SetSharedMapBit(i, kTRUE);
2059 if (s1->GetClusterIndex2(i)>0)
2060 s1->SetClusterMapBit(i, kTRUE);
2062 if (sumshared>cutN0){
2065 for (Int_t i=firstpoint;i<lastpoint;i++){
2066 // if ( (s1->GetClusterIndex2(i)&0xFFFF8FFF)==(s2->GetClusterIndex2(i)&0xFFFF8FFF) && s1->GetClusterIndex2(i)>0) {
2067 if ( (s1->GetClusterIndex2(i))==(s2->GetClusterIndex2(i)) && s1->GetClusterIndex2(i)>0) {
2068 AliTPCTrackerPoint *p1 = s1->GetTrackPoint(i);
2069 AliTPCTrackerPoint *p2 = s2->GetTrackPoint(i);;
2070 if (s1->IsActive()&&s2->IsActive()){
2071 p1->SetShared(kTRUE);
2072 p2->SetShared(kTRUE);
2078 if (sumshared>cutN0){
2079 for (Int_t i=0;i<4;i++){
2080 if (s1->GetOverlapLabel(3*i)==0){
2081 s1->SetOverlapLabel(3*i, s2->GetLabel());
2082 s1->SetOverlapLabel(3*i+1,sumshared);
2083 s1->SetOverlapLabel(3*i+2,s2->GetUniqueID());
2087 for (Int_t i=0;i<4;i++){
2088 if (s2->GetOverlapLabel(3*i)==0){
2089 s2->SetOverlapLabel(3*i, s1->GetLabel());
2090 s2->SetOverlapLabel(3*i+1,sumshared);
2091 s2->SetOverlapLabel(3*i+2,s1->GetUniqueID());
2098 void AliTPCtrackerMI::SignShared(TObjArray * arr)
2101 //sort trackss according sectors
2103 for (Int_t i=0; i<arr->GetEntriesFast(); i++) {
2104 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2106 //if (pt) RotateToLocal(pt);
2110 arr->Sort(); // sorting according relative sectors
2111 arr->Expand(arr->GetEntries());
2114 Int_t nseed=arr->GetEntriesFast();
2115 for (Int_t i=0; i<nseed; i++) {
2116 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2118 for (Int_t j=0;j<=12;j++){
2119 pt->SetOverlapLabel(j,0);
2122 for (Int_t i=0; i<nseed; i++) {
2123 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2125 if (pt->GetRemoval()>10) continue;
2126 for (Int_t j=i+1; j<nseed; j++){
2127 AliTPCseed *pt2=(AliTPCseed*)arr->UncheckedAt(j);
2128 if (TMath::Abs(pt->GetRelativeSector()-pt2->GetRelativeSector())>1) continue;
2130 if (pt2->GetRemoval()<=10) {
2131 //if ( TMath::Abs(pt->GetRelativeSector()-pt2->GetRelativeSector())>0) break;
2139 void AliTPCtrackerMI::SortTracks(TObjArray * arr, Int_t mode) const
2142 //sort tracks in array according mode criteria
2143 Int_t nseed = arr->GetEntriesFast();
2144 for (Int_t i=0; i<nseed; i++) {
2145 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2156 void AliTPCtrackerMI::RemoveUsed2(TObjArray * arr, Float_t factor1, Float_t factor2, Int_t minimal)
2159 // Loop over all tracks and remove overlaped tracks (with lower quality)
2161 // 1. Unsign clusters
2162 // 2. Sort tracks according quality
2163 // Quality is defined by the number of cluster between first and last points
2165 // 3. Loop over tracks - decreasing quality order
2166 // a.) remove - If the fraction of shared cluster less than factor (1- n or 2)
2167 // b.) remove - If the minimal number of clusters less than minimal and not ITS
2168 // c.) if track accepted - sign clusters
2170 //Called in - AliTPCtrackerMI::Clusters2Tracks()
2171 // - AliTPCtrackerMI::PropagateBack()
2172 // - AliTPCtrackerMI::RefitInward()
2178 Int_t nseed = arr->GetEntriesFast();
2179 Float_t * quality = new Float_t[nseed];
2180 Int_t * indexes = new Int_t[nseed];
2184 for (Int_t i=0; i<nseed; i++) {
2185 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2190 pt->UpdatePoints(); //select first last max dens points
2191 Float_t * points = pt->GetPoints();
2192 if (points[3]<0.8) quality[i] =-1;
2193 quality[i] = (points[2]-points[0])+pt->GetNumberOfClusters();
2194 //prefer high momenta tracks if overlaps
2195 quality[i] *= TMath::Sqrt(TMath::Abs(pt->Pt())+0.5);
2197 TMath::Sort(nseed,quality,indexes);
2200 for (Int_t itrack=0; itrack<nseed; itrack++) {
2201 Int_t trackindex = indexes[itrack];
2202 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(trackindex);
2205 if (quality[trackindex]<0){
2207 delete arr->RemoveAt(trackindex);
2210 arr->RemoveAt(trackindex);
2216 Int_t first = Int_t(pt->GetPoints()[0]);
2217 Int_t last = Int_t(pt->GetPoints()[2]);
2218 Double_t factor = (pt->GetBConstrain()) ? factor1: factor2;
2220 Int_t found,foundable,shared;
2221 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
2222 // pt->GetClusterStatistic(0,160, found, foundable,shared,kFALSE);
2223 Bool_t itsgold =kFALSE;
2226 if (pt->GetESD()->GetITSclusters(dummy)>4) itsgold= kTRUE;
2230 if (Float_t(shared+1)/Float_t(found+1)>factor){
2231 if (pt->GetKinkIndexes()[0]!=0) continue; //don't remove tracks - part of the kinks
2232 delete arr->RemoveAt(trackindex);
2235 if (pt->GetNumberOfClusters()<50&&(found-0.5*shared)<minimal){ //remove short tracks
2236 if (pt->GetKinkIndexes()[0]!=0) continue; //don't remove tracks - part of the kinks
2237 delete arr->RemoveAt(trackindex);
2243 //if (sharedfactor>0.4) continue;
2244 if (pt->GetKinkIndexes()[0]>0) continue;
2245 //Remove tracks with undefined properties - seems
2246 if (pt->GetSigmaY2()<kAlmost0) continue; // ? what is the origin ?
2248 for (Int_t i=first; i<last; i++) {
2249 Int_t index=pt->GetClusterIndex2(i);
2250 // if (index<0 || index&0x8000 ) continue;
2251 if (index<0 || index&0x8000 ) continue;
2252 AliTPCclusterMI *c= pt->GetClusterPointer(i);
2259 Info("RemoveUsed","\n*****\nNumber of good tracks after shared removal\t%d\n",fNtracks);
2265 void AliTPCtrackerMI::UnsignClusters()
2268 // loop over all clusters and unsign them
2271 for (Int_t sec=0;sec<fkNIS;sec++){
2272 for (Int_t row=0;row<fInnerSec->GetNRows();row++){
2273 AliTPCclusterMI *cl = fInnerSec[sec][row].GetClusters1();
2274 for (Int_t icl =0;icl< fInnerSec[sec][row].GetN1();icl++)
2275 // if (cl[icl].IsUsed(10))
2277 cl = fInnerSec[sec][row].GetClusters2();
2278 for (Int_t icl =0;icl< fInnerSec[sec][row].GetN2();icl++)
2279 //if (cl[icl].IsUsed(10))
2284 for (Int_t sec=0;sec<fkNOS;sec++){
2285 for (Int_t row=0;row<fOuterSec->GetNRows();row++){
2286 AliTPCclusterMI *cl = fOuterSec[sec][row].GetClusters1();
2287 for (Int_t icl =0;icl< fOuterSec[sec][row].GetN1();icl++)
2288 //if (cl[icl].IsUsed(10))
2290 cl = fOuterSec[sec][row].GetClusters2();
2291 for (Int_t icl =0;icl< fOuterSec[sec][row].GetN2();icl++)
2292 //if (cl[icl].IsUsed(10))
2301 void AliTPCtrackerMI::SignClusters(TObjArray * arr, Float_t fnumber, Float_t fdensity)
2304 //sign clusters to be "used"
2306 // snumber and sdensity sign number of sigmas - bellow mean value to be accepted
2307 // loop over "primaries"
2321 Int_t nseed = arr->GetEntriesFast();
2322 for (Int_t i=0; i<nseed; i++) {
2323 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2327 if (!(pt->IsActive())) continue;
2328 Float_t dens = pt->GetNumberOfClusters()/Float_t(pt->GetNFoundable());
2329 if ( (dens>0.7) && (pt->GetNumberOfClusters()>70)){
2331 sumdens2+= dens*dens;
2332 sumn += pt->GetNumberOfClusters();
2333 sumn2 += pt->GetNumberOfClusters()*pt->GetNumberOfClusters();
2334 Float_t chi2 = pt->GetChi2()/pt->GetNumberOfClusters();
2337 sumchi2 +=chi2*chi2;
2342 Float_t mdensity = 0.9;
2343 Float_t meann = 130;
2344 Float_t meanchi = 1;
2345 Float_t sdensity = 0.1;
2346 Float_t smeann = 10;
2347 Float_t smeanchi =0.4;
2351 mdensity = sumdens/sum;
2353 meanchi = sumchi/sum;
2355 sdensity = sumdens2/sum-mdensity*mdensity;
2357 sdensity = TMath::Sqrt(sdensity);
2361 smeann = sumn2/sum-meann*meann;
2363 smeann = TMath::Sqrt(smeann);
2367 smeanchi = sumchi2/sum - meanchi*meanchi;
2369 smeanchi = TMath::Sqrt(smeanchi);
2375 //REMOVE SHORT DELTAS or tracks going out of sensitive volume of TPC
2377 for (Int_t i=0; i<nseed; i++) {
2378 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2382 if (pt->GetBSigned()) continue;
2383 if (pt->GetBConstrain()) continue;
2384 //if (!(pt->IsActive())) continue;
2386 Int_t found,foundable,shared;
2387 pt->GetClusterStatistic(0,160,found, foundable,shared);
2388 if (shared/float(found)>0.3) {
2389 if (shared/float(found)>0.9 ){
2390 //delete arr->RemoveAt(i);
2395 Bool_t isok =kFALSE;
2396 if ( (pt->GetNShared()/pt->GetNumberOfClusters()<0.5) &&pt->GetNumberOfClusters()>60)
2398 if ((TMath::Abs(1/pt->GetC())<100.) && (pt->GetNShared()/pt->GetNumberOfClusters()<0.7))
2400 if (TMath::Abs(pt->GetZ()/pt->GetX())>1.1)
2402 if ( (TMath::Abs(pt->GetSnp()>0.7) && pt->GetD(0,0)>60.))
2406 for (Int_t i=0; i<160; i++) {
2407 Int_t index=pt->GetClusterIndex2(i);
2408 if (index<0) continue;
2409 AliTPCclusterMI *c= pt->GetClusterPointer(i);
2411 //if (!(c->IsUsed(10))) c->Use();
2418 Double_t maxchi = meanchi+2.*smeanchi;
2420 for (Int_t i=0; i<nseed; i++) {
2421 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2425 //if (!(pt->IsActive())) continue;
2426 if (pt->GetBSigned()) continue;
2427 Double_t chi = pt->GetChi2()/pt->GetNumberOfClusters();
2428 if (chi>maxchi) continue;
2431 Float_t dens = pt->GetNumberOfClusters()/Float_t(pt->GetNFoundable());
2433 //sign only tracks with enoug big density at the beginning
2435 if ((pt->GetDensityFirst(40)<0.75) && pt->GetNumberOfClusters()<meann) continue;
2438 Double_t mindens = TMath::Max(double(mdensity-sdensity*fdensity*bfactor),0.65);
2439 Double_t minn = TMath::Max(Int_t(meann-fnumber*smeann*bfactor),50);
2441 // if (pt->fBConstrain) mindens = TMath::Max(mdensity-sdensity*fdensity*bfactor,0.65);
2442 if ( (pt->GetRemoval()==10) && (pt->GetSnp()>0.8)&&(dens>mindens))
2445 if ((dens>mindens && pt->GetNumberOfClusters()>minn) && chi<maxchi ){
2446 //Int_t noc=pt->GetNumberOfClusters();
2447 pt->SetBSigned(kTRUE);
2448 for (Int_t i=0; i<160; i++) {
2450 Int_t index=pt->GetClusterIndex2(i);
2451 if (index<0) continue;
2452 AliTPCclusterMI *c= pt->GetClusterPointer(i);
2454 // if (!(c->IsUsed(10))) c->Use();
2459 // gLastCheck = nseed;
2467 void AliTPCtrackerMI::StopNotActive(TObjArray * arr, Int_t row0, Float_t th0, Float_t th1, Float_t th2) const
2469 // stop not active tracks
2470 // take th1 as threshold for number of founded to number of foundable on last 10 active rows
2471 // take th2 as threshold for number of founded to number of foundable on last 20 active rows
2472 Int_t nseed = arr->GetEntriesFast();
2474 for (Int_t i=0; i<nseed; i++) {
2475 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2479 if (!(pt->IsActive())) continue;
2480 StopNotActive(pt,row0,th0, th1,th2);
2486 void AliTPCtrackerMI::StopNotActive(AliTPCseed * seed, Int_t row0, Float_t th0, Float_t th1,
2489 // stop not active tracks
2490 // take th1 as threshold for number of founded to number of foundable on last 10 active rows
2491 // take th2 as threshold for number of founded to number of foundable on last 20 active rows
2494 Int_t foundable = 0;
2495 Int_t maxindex = seed->GetLastPoint(); //last foundable row
2496 if (seed->GetNFoundable()*th0 > seed->GetNumberOfClusters()) {
2497 seed->Desactivate(10) ;
2501 for (Int_t i=row0; i<maxindex; i++){
2502 Int_t index = seed->GetClusterIndex2(i);
2503 if (index!=-1) foundable++;
2505 if (foundable<=30) sumgood1++;
2506 if (foundable<=50) {
2513 if (foundable>=30.){
2514 if (sumgood1<(th1*30.)) seed->Desactivate(10);
2517 if (sumgood2<(th2*50.)) seed->Desactivate(10);
2521 Int_t AliTPCtrackerMI::RefitInward(AliESDEvent *event)
2524 // back propagation of ESD tracks
2527 const Int_t kMaxFriendTracks=2000;
2531 //PrepareForProlongation(fSeeds,1);
2532 PropagateForward2(fSeeds);
2533 RemoveUsed2(fSeeds,0.4,0.4,20);
2535 TObjArray arraySeed(fSeeds->GetEntries());
2536 for (Int_t i=0;i<fSeeds->GetEntries();i++) {
2537 arraySeed.AddAt(fSeeds->At(i),i);
2539 SignShared(&arraySeed);
2540 // FindCurling(fSeeds, event,2); // find multi found tracks
2541 FindSplitted(fSeeds, event,2); // find multi found tracks
2544 Int_t nseed = fSeeds->GetEntriesFast();
2545 for (Int_t i=0;i<nseed;i++){
2546 AliTPCseed * seed = (AliTPCseed*) fSeeds->UncheckedAt(i);
2547 if (!seed) continue;
2548 if (seed->GetKinkIndex(0)>0) UpdateKinkQualityD(seed); // update quality informations for kinks
2550 seed->PropagateTo(fParam->GetInnerRadiusLow());
2551 seed->UpdatePoints();
2553 AliESDtrack *esd=event->GetTrack(i);
2554 seed->CookdEdx(0.02,0.6);
2555 CookLabel(seed,0.1); //For comparison only
2556 esd->SetTPCClusterMap(seed->GetClusterMap());
2557 esd->SetTPCSharedMap(seed->GetSharedMap());
2559 if (AliTPCReconstructor::StreamLevel()>1 && seed!=0&&esd!=0) {
2560 TTreeSRedirector &cstream = *fDebugStreamer;
2567 if (seed->GetNumberOfClusters()>15){
2568 esd->UpdateTrackParams(seed,AliESDtrack::kTPCrefit);
2569 esd->SetTPCPoints(seed->GetPoints());
2570 esd->SetTPCPointsF(seed->GetNFoundable());
2571 Int_t ndedx = seed->GetNCDEDX(0)+seed->GetNCDEDX(1)+seed->GetNCDEDX(2)+seed->GetNCDEDX(3);
2572 Float_t sdedx = (seed->GetSDEDX(0)+seed->GetSDEDX(1)+seed->GetSDEDX(2)+seed->GetSDEDX(3))*0.25;
2573 Float_t dedx = seed->GetdEdx();
2574 esd->SetTPCsignal(dedx, sdedx, ndedx);
2576 // add seed to the esd track in Calib level
2578 Bool_t storeFriend = gRandom->Rndm()<(kMaxFriendTracks)/Float_t(nseed);
2579 if (AliTPCReconstructor::StreamLevel()>0 &&storeFriend){
2580 AliTPCseed * seedCopy = new AliTPCseed(*seed, kTRUE);
2581 esd->AddCalibObject(seedCopy);
2586 //printf("problem\n");
2589 //FindKinks(fSeeds,event);
2590 Info("RefitInward","Number of refitted tracks %d",ntracks);
2595 Int_t AliTPCtrackerMI::PropagateBack(AliESDEvent *event)
2598 // back propagation of ESD tracks
2604 PropagateBack(fSeeds);
2605 RemoveUsed2(fSeeds,0.4,0.4,20);
2606 //FindCurling(fSeeds, fEvent,1);
2607 FindSplitted(fSeeds, event,1); // find multi found tracks
2610 Int_t nseed = fSeeds->GetEntriesFast();
2612 for (Int_t i=0;i<nseed;i++){
2613 AliTPCseed * seed = (AliTPCseed*) fSeeds->UncheckedAt(i);
2614 if (!seed) continue;
2615 if (seed->GetKinkIndex(0)<0) UpdateKinkQualityM(seed); // update quality informations for kinks
2616 seed->UpdatePoints();
2617 AliESDtrack *esd=event->GetTrack(i);
2618 seed->CookdEdx(0.02,0.6);
2619 CookLabel(seed,0.1); //For comparison only
2620 if (seed->GetNumberOfClusters()>15){
2621 esd->UpdateTrackParams(seed,AliESDtrack::kTPCout);
2622 esd->SetTPCPoints(seed->GetPoints());
2623 esd->SetTPCPointsF(seed->GetNFoundable());
2624 Int_t ndedx = seed->GetNCDEDX(0)+seed->GetNCDEDX(1)+seed->GetNCDEDX(2)+seed->GetNCDEDX(3);
2625 Float_t sdedx = (seed->GetSDEDX(0)+seed->GetSDEDX(1)+seed->GetSDEDX(2)+seed->GetSDEDX(3))*0.25;
2626 Float_t dedx = seed->GetdEdx();
2627 esd->SetTPCsignal(dedx, sdedx, ndedx);
2629 Int_t eventnumber = event->GetEventNumberInFile();// patch 28 fev 06
2630 // This is most likely NOT the event number you'd like to use. It has nothing to do with the 'real' event number
2631 if (AliTPCReconstructor::StreamLevel()>1) {
2632 (*fDebugStreamer)<<"Cback"<<
2634 "EventNrInFile="<<eventnumber<<
2635 "\n"; // patch 28 fev 06
2639 //FindKinks(fSeeds,event);
2640 Info("PropagateBack","Number of back propagated tracks %d",ntracks);
2647 void AliTPCtrackerMI::DeleteSeeds()
2652 Int_t nseed = fSeeds->GetEntriesFast();
2653 for (Int_t i=0;i<nseed;i++){
2654 AliTPCseed * seed = (AliTPCseed*)fSeeds->At(i);
2655 if (seed) delete fSeeds->RemoveAt(i);
2662 void AliTPCtrackerMI::ReadSeeds(AliESDEvent *event, Int_t direction)
2665 //read seeds from the event
2667 Int_t nentr=event->GetNumberOfTracks();
2669 Info("ReadSeeds", "Number of ESD tracks: %d\n", nentr);
2674 fSeeds = new TObjArray(nentr);
2678 for (Int_t i=0; i<nentr; i++) {
2679 AliESDtrack *esd=event->GetTrack(i);
2680 ULong_t status=esd->GetStatus();
2681 if (!(status&AliESDtrack::kTPCin)) continue;
2682 AliTPCtrack t(*esd);
2683 t.SetNumberOfClusters(0);
2684 // AliTPCseed *seed = new AliTPCseed(t,t.GetAlpha());
2685 AliTPCseed *seed = new AliTPCseed(t/*,t.GetAlpha()*/);
2686 seed->SetUniqueID(esd->GetID());
2687 for (Int_t ikink=0;ikink<3;ikink++) {
2688 Int_t index = esd->GetKinkIndex(ikink);
2689 seed->GetKinkIndexes()[ikink] = index;
2690 if (index==0) continue;
2691 index = TMath::Abs(index);
2692 AliESDkink * kink = fEvent->GetKink(index-1);
2693 if (kink&&esd->GetKinkIndex(ikink)<0){
2694 if ((status & AliESDtrack::kTRDrefit) != 0) kink->SetStatus(1,2);
2695 if ((status & AliESDtrack::kITSout) != 0) kink->SetStatus(1,0);
2697 if (kink&&esd->GetKinkIndex(ikink)>0){
2698 if ((status & AliESDtrack::kTRDrefit) != 0) kink->SetStatus(1,6);
2699 if ((status & AliESDtrack::kITSout) != 0) kink->SetStatus(1,4);
2703 if (((status&AliESDtrack::kITSout)==0)&&(direction==1)) seed->ResetCovariance(10.);
2704 if ( direction ==2 &&(status & AliESDtrack::kTRDrefit) == 0 ) seed->ResetCovariance(10.);
2705 if ( direction ==2 && ((status & AliESDtrack::kTPCout) == 0) ) {
2710 if ( direction ==2 &&(status & AliESDtrack::kTRDrefit) > 0 ) {
2711 Double_t par0[5],par1[5],alpha,x;
2712 esd->GetInnerExternalParameters(alpha,x,par0);
2713 esd->GetExternalParameters(x,par1);
2714 Double_t delta1 = TMath::Abs(par0[4]-par1[4])/(0.000000001+TMath::Abs(par0[4]+par1[4]));
2715 Double_t delta2 = TMath::Abs(par0[3]-par1[3]);
2717 if (esd->GetTRDncls()>0) trdchi2 = esd->GetTRDchi2()/esd->GetTRDncls();
2718 //reset covariance if suspicious
2719 if ( (delta1>0.1) || (delta2>0.006) ||trdchi2>7.)
2720 seed->ResetCovariance(10.);
2725 // rotate to the local coordinate system
2727 fSectors=fInnerSec; fN=fkNIS;
2728 Double_t alpha=seed->GetAlpha() - fSectors->GetAlphaShift();
2729 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
2730 if (alpha < 0. ) alpha += 2.*TMath::Pi();
2731 Int_t ns=Int_t(alpha/fSectors->GetAlpha())%fN;
2732 alpha =ns*fSectors->GetAlpha() + fSectors->GetAlphaShift();
2733 if (alpha<-TMath::Pi()) alpha += 2*TMath::Pi();
2734 if (alpha>=TMath::Pi()) alpha -= 2*TMath::Pi();
2735 alpha-=seed->GetAlpha();
2736 if (!seed->Rotate(alpha)) {
2742 if (esd->GetKinkIndex(0)<=0){
2743 for (Int_t irow=0;irow<160;irow++){
2744 Int_t index = seed->GetClusterIndex2(irow);
2747 AliTPCclusterMI * cl = GetClusterMI(index);
2748 seed->SetClusterPointer(irow,cl);
2750 if ((index & 0x8000)==0){
2751 cl->Use(10); // accepted cluster
2753 cl->Use(6); // close cluster not accepted
2756 Info("ReadSeeds","Not found cluster");
2761 fSeeds->AddAt(seed,i);
2767 //_____________________________________________________________________________
2768 void AliTPCtrackerMI::MakeSeeds3(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, Float_t cuts[4],
2769 Float_t deltay, Int_t ddsec) {
2770 //-----------------------------------------------------------------
2771 // This function creates track seeds.
2772 // SEEDING WITH VERTEX CONSTRAIN
2773 //-----------------------------------------------------------------
2774 // cuts[0] - fP4 cut
2775 // cuts[1] - tan(phi) cut
2776 // cuts[2] - zvertex cut
2777 // cuts[3] - fP3 cut
2785 Double_t x[5], c[15];
2786 // Int_t di = i1-i2;
2788 AliTPCseed * seed = new AliTPCseed();
2789 Double_t alpha=fSectors->GetAlpha(), shift=fSectors->GetAlphaShift();
2790 Double_t cs=cos(alpha), sn=sin(alpha);
2792 // Double_t x1 =fOuterSec->GetX(i1);
2793 //Double_t xx2=fOuterSec->GetX(i2);
2795 Double_t x1 =GetXrow(i1);
2796 Double_t xx2=GetXrow(i2);
2798 Double_t x3=GetX(), y3=GetY(), z3=GetZ();
2800 Int_t imiddle = (i2+i1)/2; //middle pad row index
2801 Double_t xm = GetXrow(imiddle); // radius of middle pad-row
2802 const AliTPCtrackerRow& krm=GetRow(sec,imiddle); //middle pad -row
2806 const AliTPCtrackerRow& kr1=GetRow(ns,i1);
2807 Double_t ymax = GetMaxY(i1)-kr1.GetDeadZone()-1.5;
2808 Double_t ymaxm = GetMaxY(imiddle)-kr1.GetDeadZone()-1.5;
2811 // change cut on curvature if it can't reach this layer
2812 // maximal curvature set to reach it
2813 Double_t dvertexmax = TMath::Sqrt((x1-x3)*(x1-x3)+(ymax+5-y3)*(ymax+5-y3));
2814 if (dvertexmax*0.5*cuts[0]>0.85){
2815 cuts[0] = 0.85/(dvertexmax*0.5+1.);
2817 Double_t r2min = 1/(cuts[0]*cuts[0]); //minimal square of radius given by cut
2820 if (deltay>0) ddsec = 0;
2821 // loop over clusters
2822 for (Int_t is=0; is < kr1; is++) {
2824 if (kr1[is]->IsUsed(10)) continue;
2825 Double_t y1=kr1[is]->GetY(), z1=kr1[is]->GetZ();
2826 //if (TMath::Abs(y1)>ymax) continue;
2828 if (deltay>0 && TMath::Abs(ymax-TMath::Abs(y1))> deltay ) continue; // seed only at the edge
2830 // find possible directions
2831 Float_t anglez = (z1-z3)/(x1-x3);
2832 Float_t extraz = z1 - anglez*(x1-xx2); // extrapolated z
2835 //find rotation angles relative to line given by vertex and point 1
2836 Double_t dvertex2 = (x1-x3)*(x1-x3)+(y1-y3)*(y1-y3);
2837 Double_t dvertex = TMath::Sqrt(dvertex2);
2838 Double_t angle13 = TMath::ATan((y1-y3)/(x1-x3));
2839 Double_t cs13 = cos(-angle13), sn13 = sin(-angle13);
2842 // loop over 2 sectors
2848 Double_t dddz1=0; // direction of delta inclination in z axis
2855 for (Int_t dsec = dsec1; dsec<=dsec2;dsec++){
2856 Int_t sec2 = sec + dsec;
2858 // AliTPCtrackerRow& kr2 = fOuterSec[(sec2+fkNOS)%fkNOS][i2];
2859 //AliTPCtrackerRow& kr2m = fOuterSec[(sec2+fkNOS)%fkNOS][imiddle];
2860 AliTPCtrackerRow& kr2 = GetRow((sec2+fkNOS)%fkNOS,i2);
2861 AliTPCtrackerRow& kr2m = GetRow((sec2+fkNOS)%fkNOS,imiddle);
2862 Int_t index1 = TMath::Max(kr2.Find(extraz-0.6-dddz1*TMath::Abs(z1)*0.05)-1,0);
2863 Int_t index2 = TMath::Min(kr2.Find(extraz+0.6+dddz2*TMath::Abs(z1)*0.05)+1,kr2);
2865 // rotation angles to p1-p3
2866 Double_t cs13r = cos(-angle13+dsec*alpha)/dvertex, sn13r = sin(-angle13+dsec*alpha)/dvertex;
2867 Double_t x2, y2, z2;
2869 // Double_t dymax = maxangle*TMath::Abs(x1-xx2);
2872 Double_t dxx0 = (xx2-x3)*cs13r;
2873 Double_t dyy0 = (xx2-x3)*sn13r;
2874 for (Int_t js=index1; js < index2; js++) {
2875 const AliTPCclusterMI *kcl = kr2[js];
2876 if (kcl->IsUsed(10)) continue;
2878 //calcutate parameters
2880 Double_t yy0 = dyy0 +(kcl->GetY()-y3)*cs13r;
2882 if (TMath::Abs(yy0)<0.000001) continue;
2883 Double_t xx0 = dxx0 -(kcl->GetY()-y3)*sn13r;
2884 Double_t y0 = 0.5*(xx0*xx0+yy0*yy0-xx0)/yy0;
2885 Double_t r02 = (0.25+y0*y0)*dvertex2;
2886 //curvature (radius) cut
2887 if (r02<r2min) continue;
2891 Double_t c0 = 1/TMath::Sqrt(r02);
2895 //Double_t dfi0 = 2.*TMath::ASin(dvertex*c0*0.5);
2896 //Double_t dfi1 = 2.*TMath::ASin(TMath::Sqrt(yy0*yy0+(1-xx0)*(1-xx0))*dvertex*c0*0.5);
2897 Double_t dfi0 = 2.*AliTPCFastMath::FastAsin(dvertex*c0*0.5);
2898 Double_t dfi1 = 2.*AliTPCFastMath::FastAsin(TMath::Sqrt(yy0*yy0+(1-xx0)*(1-xx0))*dvertex*c0*0.5);
2901 Double_t z0 = kcl->GetZ();
2902 Double_t zzzz2 = z1-(z1-z3)*dfi1/dfi0;
2903 if (TMath::Abs(zzzz2-z0)>0.5) continue;
2906 Double_t dip = (z1-z0)*c0/dfi1;
2907 Double_t x0 = (0.5*cs13+y0*sn13)*dvertex*c0;
2918 x2= xx2*cs-y2*sn*dsec;
2919 y2=+xx2*sn*dsec+y2*cs;
2929 // do we have cluster at the middle ?
2931 GetProlongation(x1,xm,x,ym,zm);
2933 AliTPCclusterMI * cm=0;
2934 if (TMath::Abs(ym)-ymaxm<0){
2935 cm = krm.FindNearest2(ym,zm,1.0,0.6,dummy);
2936 if ((!cm) || (cm->IsUsed(10))) {
2941 // rotate y1 to system 0
2942 // get state vector in rotated system
2943 Double_t yr1 = (-0.5*sn13+y0*cs13)*dvertex*c0;
2944 Double_t xr2 = x0*cs+yr1*sn*dsec;
2945 Double_t xr[5]={kcl->GetY(),kcl->GetZ(), xr2, dip, c0};
2947 GetProlongation(xx2,xm,xr,ym,zm);
2948 if (TMath::Abs(ym)-ymaxm<0){
2949 cm = kr2m.FindNearest2(ym,zm,1.0,0.6,dummy);
2950 if ((!cm) || (cm->IsUsed(10))) {
2960 dym = ym - cm->GetY();
2961 dzm = zm - cm->GetZ();
2968 Double_t sy1=kr1[is]->GetSigmaY2()*2., sz1=kr1[is]->GetSigmaZ2()*2.;
2969 Double_t sy2=kcl->GetSigmaY2()*2., sz2=kcl->GetSigmaZ2()*2.;
2970 //Double_t sy3=400*3./12., sy=0.1, sz=0.1;
2971 Double_t sy3=25000*x[4]*x[4]+0.1, sy=0.1, sz=0.1;
2972 //Double_t sy3=25000*x[4]*x[4]*60+0.5, sy=0.1, sz=0.1;
2974 Double_t f40=(F1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
2975 Double_t f42=(F1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
2976 Double_t f43=(F1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
2977 Double_t f20=(F2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
2978 Double_t f22=(F2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
2979 Double_t f23=(F2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
2981 Double_t f30=(F3(x1,y1+sy,x2,y2,z1,z2)-x[3])/sy;
2982 Double_t f31=(F3(x1,y1,x2,y2,z1+sz,z2)-x[3])/sz;
2983 Double_t f32=(F3(x1,y1,x2,y2+sy,z1,z2)-x[3])/sy;
2984 Double_t f34=(F3(x1,y1,x2,y2,z1,z2+sz)-x[3])/sz;
2988 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
2989 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
2990 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
2991 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
2992 c[13]=f30*sy1*f40+f32*sy2*f42;
2993 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
2995 // if (!BuildSeed(kr1[is],kcl,0,x1,x2,x3,x,c)) continue;
2997 UInt_t index=kr1.GetIndex(is);
2998 seed->~AliTPCseed(); // this does not set the pointer to 0...
2999 AliTPCseed *track=new(seed) AliTPCseed(x1, ns*alpha+shift, x, c, index);
3001 track->SetIsSeeding(kTRUE);
3002 track->SetSeed1(i1);
3003 track->SetSeed2(i2);
3004 track->SetSeedType(3);
3008 FollowProlongation(*track, (i1+i2)/2,1);
3009 Int_t foundable,found,shared;
3010 track->GetClusterStatistic((i1+i2)/2,i1, found, foundable, shared, kTRUE);
3011 if ((found<0.55*foundable) || shared>0.5*found || (track->GetSigmaY2()+track->GetSigmaZ2())>0.5){
3013 seed->~AliTPCseed();
3019 FollowProlongation(*track, i2,1);
3023 track->SetBConstrain(1);
3024 // track->fLastPoint = i1+fInnerSec->GetNRows(); // first cluster in track position
3025 track->SetLastPoint(i1); // first cluster in track position
3026 track->SetFirstPoint(track->GetLastPoint());
3028 if (track->GetNumberOfClusters()<(i1-i2)*0.5 ||
3029 track->GetNumberOfClusters() < track->GetNFoundable()*0.6 ||
3030 track->GetNShared()>0.4*track->GetNumberOfClusters() ) {
3032 seed->~AliTPCseed();
3036 // Z VERTEX CONDITION
3037 Double_t zv, bz=GetBz();
3038 if ( !track->GetZAt(0.,bz,zv) ) continue;
3039 if (TMath::Abs(zv-z3)>cuts[2]) {
3040 FollowProlongation(*track, TMath::Max(i2-20,0));
3041 if ( !track->GetZAt(0.,bz,zv) ) continue;
3042 if (TMath::Abs(zv-z3)>cuts[2]){
3043 FollowProlongation(*track, TMath::Max(i2-40,0));
3044 if ( !track->GetZAt(0.,bz,zv) ) continue;
3045 if (TMath::Abs(zv-z3)>cuts[2] &&(track->GetNumberOfClusters() > track->GetNFoundable()*0.7)){
3046 // make seed without constrain
3047 AliTPCseed * track2 = MakeSeed(track,0.2,0.5,1.);
3048 FollowProlongation(*track2, i2,1);
3049 track2->SetBConstrain(kFALSE);
3050 track2->SetSeedType(1);
3051 arr->AddLast(track2);
3053 seed->~AliTPCseed();
3058 seed->~AliTPCseed();
3065 track->SetSeedType(0);
3066 arr->AddLast(track);
3067 seed = new AliTPCseed;
3069 // don't consider other combinations
3070 if (track->GetNumberOfClusters() > track->GetNFoundable()*0.8)
3076 Info("MakeSeeds3","\nSeeding statistic:\t%d\t%d\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin,nout1,nout2);
3082 void AliTPCtrackerMI::MakeSeeds5(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, Float_t cuts[4],
3087 //-----------------------------------------------------------------
3088 // This function creates track seeds.
3089 //-----------------------------------------------------------------
3090 // cuts[0] - fP4 cut
3091 // cuts[1] - tan(phi) cut
3092 // cuts[2] - zvertex cut
3093 // cuts[3] - fP3 cut
3103 Double_t x[5], c[15];
3105 // make temporary seed
3106 AliTPCseed * seed = new AliTPCseed;
3107 Double_t alpha=fOuterSec->GetAlpha(), shift=fOuterSec->GetAlphaShift();
3108 // Double_t cs=cos(alpha), sn=sin(alpha);
3113 Double_t x1 = GetXrow(i1-1);
3114 const AliTPCtrackerRow& kr1=GetRow(sec,i1-1);
3115 Double_t y1max = GetMaxY(i1-1)-kr1.GetDeadZone()-1.5;
3117 Double_t x1p = GetXrow(i1);
3118 const AliTPCtrackerRow& kr1p=GetRow(sec,i1);
3120 Double_t x1m = GetXrow(i1-2);
3121 const AliTPCtrackerRow& kr1m=GetRow(sec,i1-2);
3124 //last 3 padrow for seeding
3125 AliTPCtrackerRow& kr3 = GetRow((sec+fkNOS)%fkNOS,i1-7);
3126 Double_t x3 = GetXrow(i1-7);
3127 // Double_t y3max= GetMaxY(i1-7)-kr3.fDeadZone-1.5;
3129 AliTPCtrackerRow& kr3p = GetRow((sec+fkNOS)%fkNOS,i1-6);
3130 Double_t x3p = GetXrow(i1-6);
3132 AliTPCtrackerRow& kr3m = GetRow((sec+fkNOS)%fkNOS,i1-8);
3133 Double_t x3m = GetXrow(i1-8);
3138 Int_t im = i1-4; //middle pad row index
3139 Double_t xm = GetXrow(im); // radius of middle pad-row
3140 const AliTPCtrackerRow& krm=GetRow(sec,im); //middle pad -row
3141 // Double_t ymmax = GetMaxY(im)-kr1.fDeadZone-1.5;
3144 Double_t deltax = x1-x3;
3145 Double_t dymax = deltax*cuts[1];
3146 Double_t dzmax = deltax*cuts[3];
3148 // loop over clusters
3149 for (Int_t is=0; is < kr1; is++) {
3151 if (kr1[is]->IsUsed(10)) continue;
3152 Double_t y1=kr1[is]->GetY(), z1=kr1[is]->GetZ();
3154 if (deltay>0 && TMath::Abs(y1max-TMath::Abs(y1))> deltay ) continue; // seed only at the edge
3156 Int_t index1 = TMath::Max(kr3.Find(z1-dzmax)-1,0);
3157 Int_t index2 = TMath::Min(kr3.Find(z1+dzmax)+1,kr3);
3163 for (Int_t js=index1; js < index2; js++) {
3164 const AliTPCclusterMI *kcl = kr3[js];
3165 if (kcl->IsUsed(10)) continue;
3167 // apply angular cuts
3168 if (TMath::Abs(y1-y3)>dymax) continue;
3171 if (TMath::Abs(z1-z3)>dzmax) continue;
3173 Double_t angley = (y1-y3)/(x1-x3);
3174 Double_t anglez = (z1-z3)/(x1-x3);
3176 Double_t erry = TMath::Abs(angley)*(x1-x1m)*0.5+0.5;
3177 Double_t errz = TMath::Abs(anglez)*(x1-x1m)*0.5+0.5;
3179 Double_t yyym = angley*(xm-x1)+y1;
3180 Double_t zzzm = anglez*(xm-x1)+z1;
3182 const AliTPCclusterMI *kcm = krm.FindNearest2(yyym,zzzm,erry,errz,index);
3184 if (kcm->IsUsed(10)) continue;
3186 erry = TMath::Abs(angley)*(x1-x1m)*0.4+0.5;
3187 errz = TMath::Abs(anglez)*(x1-x1m)*0.4+0.5;
3194 // look around first
3195 const AliTPCclusterMI *kc1m = kr1m.FindNearest2(angley*(x1m-x1)+y1,
3201 if (kc1m->IsUsed(10)) used++;
3203 const AliTPCclusterMI *kc1p = kr1p.FindNearest2(angley*(x1p-x1)+y1,
3209 if (kc1p->IsUsed(10)) used++;
3211 if (used>1) continue;
3212 if (found<1) continue;
3216 const AliTPCclusterMI *kc3m = kr3m.FindNearest2(angley*(x3m-x3)+y3,
3222 if (kc3m->IsUsed(10)) used++;
3226 const AliTPCclusterMI *kc3p = kr3p.FindNearest2(angley*(x3p-x3)+y3,
3232 if (kc3p->IsUsed(10)) used++;
3236 if (used>1) continue;
3237 if (found<3) continue;
3247 x[4]=F1(x1,y1,x2,y2,x3,y3);
3248 //if (TMath::Abs(x[4]) >= cuts[0]) continue;
3251 x[2]=F2(x1,y1,x2,y2,x3,y3);
3254 x[3]=F3n(x1,y1,x2,y2,z1,z2,x[4]);
3255 //if (TMath::Abs(x[3]) > cuts[3]) continue;
3259 Double_t sy1=0.1, sz1=0.1;
3260 Double_t sy2=0.1, sz2=0.1;
3261 Double_t sy3=0.1, sy=0.1, sz=0.1;
3263 Double_t f40=(F1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
3264 Double_t f42=(F1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
3265 Double_t f43=(F1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
3266 Double_t f20=(F2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
3267 Double_t f22=(F2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
3268 Double_t f23=(F2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
3270 Double_t f30=(F3(x1,y1+sy,x2,y2,z1,z2)-x[3])/sy;
3271 Double_t f31=(F3(x1,y1,x2,y2,z1+sz,z2)-x[3])/sz;
3272 Double_t f32=(F3(x1,y1,x2,y2+sy,z1,z2)-x[3])/sy;
3273 Double_t f34=(F3(x1,y1,x2,y2,z1,z2+sz)-x[3])/sz;
3277 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3278 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3279 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3280 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3281 c[13]=f30*sy1*f40+f32*sy2*f42;
3282 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3284 // if (!BuildSeed(kr1[is],kcl,0,x1,x2,x3,x,c)) continue;
3286 UInt_t index=kr1.GetIndex(is);
3287 seed->~AliTPCseed();
3288 AliTPCseed *track=new(seed) AliTPCseed(x1, sec*alpha+shift, x, c, index);
3290 track->SetIsSeeding(kTRUE);
3293 FollowProlongation(*track, i1-7,1);
3294 if (track->GetNumberOfClusters() < track->GetNFoundable()*0.75 ||
3295 track->GetNShared()>0.6*track->GetNumberOfClusters() || ( track->GetSigmaY2()+ track->GetSigmaZ2())>0.6){
3297 seed->~AliTPCseed();
3303 FollowProlongation(*track, i2,1);
3304 track->SetBConstrain(0);
3305 track->SetLastPoint(i1+fInnerSec->GetNRows()); // first cluster in track position
3306 track->SetFirstPoint(track->GetLastPoint());
3308 if (track->GetNumberOfClusters()<(i1-i2)*0.5 ||
3309 track->GetNumberOfClusters()<track->GetNFoundable()*0.7 ||
3310 track->GetNShared()>2. || track->GetChi2()/track->GetNumberOfClusters()>6 || ( track->GetSigmaY2()+ track->GetSigmaZ2())>0.5 ) {
3312 seed->~AliTPCseed();
3317 FollowProlongation(*track, TMath::Max(i2-10,0),1);
3318 AliTPCseed * track2 = MakeSeed(track,0.2,0.5,0.9);
3319 FollowProlongation(*track2, i2,1);
3320 track2->SetBConstrain(kFALSE);
3321 track2->SetSeedType(4);
3322 arr->AddLast(track2);
3324 seed->~AliTPCseed();
3328 //arr->AddLast(track);
3329 //seed = new AliTPCseed;
3335 Info("MakeSeeds5","\nSeeding statiistic:\t%d\t%d\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin,nout1,nout2,nout3);
3341 //_____________________________________________________________________________
3342 void AliTPCtrackerMI::MakeSeeds2(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, Float_t */*cuts[4]*/,
3343 Float_t deltay, Bool_t /*bconstrain*/) {
3344 //-----------------------------------------------------------------
3345 // This function creates track seeds - without vertex constraint
3346 //-----------------------------------------------------------------
3347 // cuts[0] - fP4 cut - not applied
3348 // cuts[1] - tan(phi) cut
3349 // cuts[2] - zvertex cut - not applied
3350 // cuts[3] - fP3 cut
3360 Double_t alpha=fOuterSec->GetAlpha(), shift=fOuterSec->GetAlphaShift();
3361 // Double_t cs=cos(alpha), sn=sin(alpha);
3362 Int_t row0 = (i1+i2)/2;
3363 Int_t drow = (i1-i2)/2;
3364 const AliTPCtrackerRow& kr0=fSectors[sec][row0];
3365 AliTPCtrackerRow * kr=0;
3367 AliTPCpolyTrack polytrack;
3368 Int_t nclusters=fSectors[sec][row0];
3369 AliTPCseed * seed = new AliTPCseed;
3374 for (Int_t is=0; is < nclusters; is++) { //LOOP over clusters
3376 Int_t nfoundable =0;
3377 for (Int_t iter =1; iter<2; iter++){ //iterations
3378 const AliTPCtrackerRow& krm=fSectors[sec][row0-iter];
3379 const AliTPCtrackerRow& krp=fSectors[sec][row0+iter];
3380 const AliTPCclusterMI * cl= kr0[is];
3382 if (cl->IsUsed(10)) {
3388 Double_t x = kr0.GetX();
3389 // Initialization of the polytrack
3394 Double_t y0= cl->GetY();
3395 Double_t z0= cl->GetZ();
3399 Double_t ymax = fSectors->GetMaxY(row0)-kr0.GetDeadZone()-1.5;
3400 if (deltay>0 && TMath::Abs(ymax-TMath::Abs(y0))> deltay ) continue; // seed only at the edge
3402 erry = (0.5)*cl->GetSigmaY2()/TMath::Sqrt(cl->GetQ())*6;
3403 errz = (0.5)*cl->GetSigmaZ2()/TMath::Sqrt(cl->GetQ())*6;
3404 polytrack.AddPoint(x,y0,z0,erry, errz);
3407 if (cl->IsUsed(10)) sumused++;
3410 Float_t roady = (5*TMath::Sqrt(cl->GetSigmaY2()+0.2)+1.)*iter;
3411 Float_t roadz = (5*TMath::Sqrt(cl->GetSigmaZ2()+0.2)+1.)*iter;
3414 AliTPCclusterMI * cl1 = krm.FindNearest(y0,z0,roady,roadz);
3415 if (cl1 && TMath::Abs(ymax-TMath::Abs(y0))) {
3416 erry = (0.5)*cl1->GetSigmaY2()/TMath::Sqrt(cl1->GetQ())*3;
3417 errz = (0.5)*cl1->GetSigmaZ2()/TMath::Sqrt(cl1->GetQ())*3;
3418 if (cl1->IsUsed(10)) sumused++;
3419 polytrack.AddPoint(x,cl1->GetY(),cl1->GetZ(),erry,errz);
3423 AliTPCclusterMI * cl2 = krp.FindNearest(y0,z0,roady,roadz);
3425 erry = (0.5)*cl2->GetSigmaY2()/TMath::Sqrt(cl2->GetQ())*3;
3426 errz = (0.5)*cl2->GetSigmaZ2()/TMath::Sqrt(cl2->GetQ())*3;
3427 if (cl2->IsUsed(10)) sumused++;
3428 polytrack.AddPoint(x,cl2->GetY(),cl2->GetZ(),erry,errz);
3431 if (sumused>0) continue;
3433 polytrack.UpdateParameters();
3439 nfoundable = polytrack.GetN();
3440 nfound = nfoundable;
3442 for (Int_t ddrow = iter+1; ddrow<drow;ddrow++){
3443 Float_t maxdist = 0.8*(1.+3./(ddrow));
3444 for (Int_t delta = -1;delta<=1;delta+=2){
3445 Int_t row = row0+ddrow*delta;
3446 kr = &(fSectors[sec][row]);
3447 Double_t xn = kr->GetX();
3448 Double_t ymax = fSectors->GetMaxY(row)-kr->GetDeadZone()-1.5;
3449 polytrack.GetFitPoint(xn,yn,zn);
3450 if (TMath::Abs(yn)>ymax) continue;
3452 AliTPCclusterMI * cln = kr->FindNearest(yn,zn,roady,roadz);
3454 Float_t dist = TMath::Sqrt( (yn-cln->GetY())*(yn-cln->GetY())+(zn-cln->GetZ())*(zn-cln->GetZ()));
3457 erry = (dist+0.3)*cln->GetSigmaY2()/TMath::Sqrt(cln->GetQ())*(1.+1./(ddrow));
3458 errz = (dist+0.3)*cln->GetSigmaZ2()/TMath::Sqrt(cln->GetQ())*(1.+1./(ddrow));
3459 if (cln->IsUsed(10)) {
3460 // printf("used\n");
3468 polytrack.AddPoint(xn,cln->GetY(),cln->GetZ(),erry, errz);
3473 if ( (sumused>3) || (sumused>0.5*nfound) || (nfound<0.6*nfoundable)) break;
3474 polytrack.UpdateParameters();
3477 if ( (sumused>3) || (sumused>0.5*nfound)) {
3478 //printf("sumused %d\n",sumused);
3483 polytrack.GetFitDerivation(kr0.GetX(),dy,dz);
3484 AliTPCpolyTrack track2;
3486 polytrack.Refit(track2,0.5+TMath::Abs(dy)*0.3,0.4+TMath::Abs(dz)*0.3);
3487 if (track2.GetN()<0.5*nfoundable) continue;
3490 if ((nfound>0.6*nfoundable) &&( nfoundable>0.4*(i1-i2))) {
3492 // test seed with and without constrain
3493 for (Int_t constrain=0; constrain<=0;constrain++){
3494 // add polytrack candidate
3496 Double_t x[5], c[15];
3497 Double_t x1,x2,x3,y1,y2,y3,z1,z2,z3;
3498 track2.GetBoundaries(x3,x1);
3500 track2.GetFitPoint(x1,y1,z1);
3501 track2.GetFitPoint(x2,y2,z2);
3502 track2.GetFitPoint(x3,y3,z3);
3504 //is track pointing to the vertex ?
3507 polytrack.GetFitPoint(x0,y0,z0);
3520 x[4]=F1(x1,y1,x2,y2,x3,y3);
3522 // if (TMath::Abs(x[4]) >= cuts[0]) continue; //
3523 x[2]=F2(x1,y1,x2,y2,x3,y3);
3525 //if (TMath::Abs(x[4]*x1-x[2]) >= cuts[1]) continue;
3526 //x[3]=F3(x1,y1,x2,y2,z1,z2);
3527 x[3]=F3n(x1,y1,x3,y3,z1,z3,x[4]);
3528 //if (TMath::Abs(x[3]) > cuts[3]) continue;
3531 Double_t sy =0.1, sz =0.1;
3532 Double_t sy1=0.02, sz1=0.02;
3533 Double_t sy2=0.02, sz2=0.02;
3537 sy3=25000*x[4]*x[4]+0.1, sy=0.1, sz=0.1;
3540 Double_t f40=(F1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
3541 Double_t f42=(F1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
3542 Double_t f43=(F1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
3543 Double_t f20=(F2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
3544 Double_t f22=(F2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
3545 Double_t f23=(F2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
3547 Double_t f30=(F3(x1,y1+sy,x3,y3,z1,z3)-x[3])/sy;
3548 Double_t f31=(F3(x1,y1,x3,y3,z1+sz,z3)-x[3])/sz;
3549 Double_t f32=(F3(x1,y1,x3,y3+sy,z1,z3)-x[3])/sy;
3550 Double_t f34=(F3(x1,y1,x3,y3,z1,z3+sz)-x[3])/sz;
3555 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3556 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3557 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3558 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3559 c[13]=f30*sy1*f40+f32*sy2*f42;
3560 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3562 //Int_t row1 = fSectors->GetRowNumber(x1);
3563 Int_t row1 = GetRowNumber(x1);
3567 seed->~AliTPCseed();
3568 AliTPCseed *track=new(seed) AliTPCseed(x1,sec*alpha+shift,x,c,index);
3569 track->SetIsSeeding(kTRUE);
3570 Int_t rc=FollowProlongation(*track, i2);
3571 if (constrain) track->SetBConstrain(1);
3573 track->SetBConstrain(0);
3574 track->SetLastPoint(row1+fInnerSec->GetNRows()); // first cluster in track position
3575 track->SetFirstPoint(track->GetLastPoint());
3577 if (rc==0 || track->GetNumberOfClusters()<(i1-i2)*0.5 ||
3578 track->GetNumberOfClusters() < track->GetNFoundable()*0.6 ||
3579 track->GetNShared()>0.4*track->GetNumberOfClusters()) {
3582 seed->~AliTPCseed();
3585 arr->AddLast(track);
3586 seed = new AliTPCseed;
3590 } // if accepted seed
3593 Info("MakeSeeds2","\nSeeding statiistic:\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin3);
3599 AliTPCseed *AliTPCtrackerMI::MakeSeed(AliTPCseed *track, Float_t r0, Float_t r1, Float_t r2)
3603 //reseed using track points
3604 Int_t p0 = int(r0*track->GetNumberOfClusters()); // point 0
3605 Int_t p1 = int(r1*track->GetNumberOfClusters());
3606 Int_t p2 = int(r2*track->GetNumberOfClusters()); // last point
3608 Double_t x0[3],x1[3],x2[3];
3609 for (Int_t i=0;i<3;i++){
3615 // find track position at given ratio of the length
3616 Int_t sec0=0, sec1=0, sec2=0;
3619 for (Int_t i=0;i<160;i++){
3620 if (track->GetClusterPointer(i)){
3622 AliTPCTrackerPoint *trpoint =track->GetTrackPoint(i);
3623 if ( (index<p0) || x0[0]<0 ){
3624 if (trpoint->GetX()>1){
3625 clindex = track->GetClusterIndex2(i);
3627 x0[0] = trpoint->GetX();
3628 x0[1] = trpoint->GetY();
3629 x0[2] = trpoint->GetZ();
3630 sec0 = ((clindex&0xff000000)>>24)%18;
3635 if ( (index<p1) &&(trpoint->GetX()>1)){
3636 clindex = track->GetClusterIndex2(i);
3638 x1[0] = trpoint->GetX();
3639 x1[1] = trpoint->GetY();
3640 x1[2] = trpoint->GetZ();
3641 sec1 = ((clindex&0xff000000)>>24)%18;
3644 if ( (index<p2) &&(trpoint->GetX()>1)){
3645 clindex = track->GetClusterIndex2(i);
3647 x2[0] = trpoint->GetX();
3648 x2[1] = trpoint->GetY();
3649 x2[2] = trpoint->GetZ();
3650 sec2 = ((clindex&0xff000000)>>24)%18;
3657 Double_t alpha, cs,sn, xx2,yy2;
3659 alpha = (sec1-sec2)*fSectors->GetAlpha();
3660 cs = TMath::Cos(alpha);
3661 sn = TMath::Sin(alpha);
3662 xx2= x1[0]*cs-x1[1]*sn;
3663 yy2= x1[0]*sn+x1[1]*cs;
3667 alpha = (sec0-sec2)*fSectors->GetAlpha();
3668 cs = TMath::Cos(alpha);
3669 sn = TMath::Sin(alpha);
3670 xx2= x0[0]*cs-x0[1]*sn;
3671 yy2= x0[0]*sn+x0[1]*cs;
3677 Double_t x[5],c[15];
3681 x[4]=F1(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]);
3682 // if (x[4]>1) return 0;
3683 x[2]=F2(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]);
3684 x[3]=F3n(x2[0],x2[1],x0[0],x0[1],x2[2],x0[2],x[4]);
3685 //if (TMath::Abs(x[3]) > 2.2) return 0;
3686 //if (TMath::Abs(x[2]) > 1.99) return 0;
3688 Double_t sy =0.1, sz =0.1;
3690 Double_t sy1=0.02+track->GetSigmaY2(), sz1=0.02+track->GetSigmaZ2();
3691 Double_t sy2=0.01+track->GetSigmaY2(), sz2=0.01+track->GetSigmaZ2();
3692 Double_t sy3=0.01+track->GetSigmaY2();
3694 Double_t f40=(F1(x2[0],x2[1]+sy,x1[0],x1[1],x0[0],x0[1])-x[4])/sy;
3695 Double_t f42=(F1(x2[0],x2[1],x1[0],x1[1]+sy,x0[0],x0[1])-x[4])/sy;
3696 Double_t f43=(F1(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]+sy)-x[4])/sy;
3697 Double_t f20=(F2(x2[0],x2[1]+sy,x1[0],x1[1],x0[0],x0[1])-x[2])/sy;
3698 Double_t f22=(F2(x2[0],x2[1],x1[0],x1[1]+sy,x0[0],x0[1])-x[2])/sy;
3699 Double_t f23=(F2(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]+sy)-x[2])/sy;
3701 Double_t f30=(F3(x2[0],x2[1]+sy,x0[0],x0[1],x2[2],x0[2])-x[3])/sy;
3702 Double_t f31=(F3(x2[0],x2[1],x0[0],x0[1],x2[2]+sz,x0[2])-x[3])/sz;
3703 Double_t f32=(F3(x2[0],x2[1],x0[0],x0[1]+sy,x2[2],x0[2])-x[3])/sy;
3704 Double_t f34=(F3(x2[0],x2[1],x0[0],x0[1],x2[2],x0[2]+sz)-x[3])/sz;
3709 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3710 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3711 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3712 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3713 c[13]=f30*sy1*f40+f32*sy2*f42;
3714 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3716 // Int_t row1 = fSectors->GetRowNumber(x2[0]);
3717 AliTPCseed *seed=new AliTPCseed(x2[0], sec2*fSectors->GetAlpha()+fSectors->GetAlphaShift(), x, c, 0);
3718 // Double_t y0,z0,y1,z1, y2,z2;
3719 //seed->GetProlongation(x0[0],y0,z0);
3720 // seed->GetProlongation(x1[0],y1,z1);
3721 //seed->GetProlongation(x2[0],y2,z2);
3723 seed->SetLastPoint(pp2);
3724 seed->SetFirstPoint(pp2);
3731 AliTPCseed *AliTPCtrackerMI::ReSeed(AliTPCseed *track, Float_t r0, Float_t r1, Float_t r2)
3735 //reseed using founded clusters
3737 // Find the number of clusters
3738 Int_t nclusters = 0;
3739 for (Int_t irow=0;irow<160;irow++){
3740 if (track->GetClusterIndex(irow)>0) nclusters++;
3744 ipos[0] = TMath::Max(int(r0*nclusters),0); // point 0 cluster
3745 ipos[1] = TMath::Min(int(r1*nclusters),nclusters-1); //
3746 ipos[2] = TMath::Min(int(r2*nclusters),nclusters-1); // last point
3750 Int_t row[3],sec[3]={0,0,0};
3752 // find track row position at given ratio of the length
3754 for (Int_t irow=0;irow<160;irow++){
3755 if (track->GetClusterIndex2(irow)<0) continue;
3757 for (Int_t ipoint=0;ipoint<3;ipoint++){
3758 if (index<=ipos[ipoint]) row[ipoint] = irow;
3762 //Get cluster and sector position
3763 for (Int_t ipoint=0;ipoint<3;ipoint++){
3764 Int_t clindex = track->GetClusterIndex2(row[ipoint]);
3765 AliTPCclusterMI * cl = GetClusterMI(clindex);
3768 // AliTPCclusterMI * cl = GetClusterMI(clindex);
3771 sec[ipoint] = ((clindex&0xff000000)>>24)%18;
3772 xyz[ipoint][0] = GetXrow(row[ipoint]);
3773 xyz[ipoint][1] = cl->GetY();
3774 xyz[ipoint][2] = cl->GetZ();
3778 // Calculate seed state vector and covariance matrix
3780 Double_t alpha, cs,sn, xx2,yy2;
3782 alpha = (sec[1]-sec[2])*fSectors->GetAlpha();
3783 cs = TMath::Cos(alpha);
3784 sn = TMath::Sin(alpha);
3785 xx2= xyz[1][0]*cs-xyz[1][1]*sn;
3786 yy2= xyz[1][0]*sn+xyz[1][1]*cs;
3790 alpha = (sec[0]-sec[2])*fSectors->GetAlpha();
3791 cs = TMath::Cos(alpha);
3792 sn = TMath::Sin(alpha);
3793 xx2= xyz[0][0]*cs-xyz[0][1]*sn;
3794 yy2= xyz[0][0]*sn+xyz[0][1]*cs;
3800 Double_t x[5],c[15];
3804 x[4]=F1(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
3805 x[2]=F2(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
3806 x[3]=F3n(xyz[2][0],xyz[2][1],xyz[0][0],xyz[0][1],xyz[2][2],xyz[0][2],x[4]);
3808 Double_t sy =0.1, sz =0.1;
3810 Double_t sy1=0.2, sz1=0.2;
3811 Double_t sy2=0.2, sz2=0.2;
3814 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;
3815 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;
3816 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;
3817 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;
3818 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;
3819 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;
3821 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;
3822 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;
3823 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;
3824 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;
3829 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3830 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3831 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3832 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3833 c[13]=f30*sy1*f40+f32*sy2*f42;
3834 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3836 // Int_t row1 = fSectors->GetRowNumber(xyz[2][0]);
3837 AliTPCseed *seed=new AliTPCseed(xyz[2][0], sec[2]*fSectors->GetAlpha()+fSectors->GetAlphaShift(), x, c, 0);
3838 seed->SetLastPoint(row[2]);
3839 seed->SetFirstPoint(row[2]);
3844 AliTPCseed *AliTPCtrackerMI::ReSeed(AliTPCseed *track,Int_t r0, Bool_t forward)
3848 //reseed using founded clusters
3851 Int_t row[3]={0,0,0};
3852 Int_t sec[3]={0,0,0};
3854 // forward direction
3856 for (Int_t irow=r0;irow<160;irow++){
3857 if (track->GetClusterIndex(irow)>0){
3862 for (Int_t irow=160;irow>r0;irow--){
3863 if (track->GetClusterIndex(irow)>0){
3868 for (Int_t irow=row[2]-15;irow>row[0];irow--){
3869 if (track->GetClusterIndex(irow)>0){
3877 for (Int_t irow=0;irow<r0;irow++){
3878 if (track->GetClusterIndex(irow)>0){
3883 for (Int_t irow=r0;irow>0;irow--){
3884 if (track->GetClusterIndex(irow)>0){
3889 for (Int_t irow=row[2]-15;irow>row[0];irow--){
3890 if (track->GetClusterIndex(irow)>0){
3897 if ((row[2]-row[0])<20) return 0;
3898 if (row[1]==0) return 0;
3901 //Get cluster and sector position
3902 for (Int_t ipoint=0;ipoint<3;ipoint++){
3903 Int_t clindex = track->GetClusterIndex2(row[ipoint]);
3904 AliTPCclusterMI * cl = GetClusterMI(clindex);
3907 // AliTPCclusterMI * cl = GetClusterMI(clindex);
3910 sec[ipoint] = ((clindex&0xff000000)>>24)%18;
3911 xyz[ipoint][0] = GetXrow(row[ipoint]);
3912 AliTPCTrackerPoint * point = track->GetTrackPoint(row[ipoint]);
3913 if (point&&ipoint<2){
3915 xyz[ipoint][1] = point->GetY();
3916 xyz[ipoint][2] = point->GetZ();
3919 xyz[ipoint][1] = cl->GetY();
3920 xyz[ipoint][2] = cl->GetZ();
3927 // Calculate seed state vector and covariance matrix
3929 Double_t alpha, cs,sn, xx2,yy2;
3931 alpha = (sec[1]-sec[2])*fSectors->GetAlpha();
3932 cs = TMath::Cos(alpha);
3933 sn = TMath::Sin(alpha);
3934 xx2= xyz[1][0]*cs-xyz[1][1]*sn;
3935 yy2= xyz[1][0]*sn+xyz[1][1]*cs;
3939 alpha = (sec[0]-sec[2])*fSectors->GetAlpha();
3940 cs = TMath::Cos(alpha);
3941 sn = TMath::Sin(alpha);
3942 xx2= xyz[0][0]*cs-xyz[0][1]*sn;
3943 yy2= xyz[0][0]*sn+xyz[0][1]*cs;
3949 Double_t x[5],c[15];
3953 x[4]=F1(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
3954 x[2]=F2(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
3955 x[3]=F3n(xyz[2][0],xyz[2][1],xyz[0][0],xyz[0][1],xyz[2][2],xyz[0][2],x[4]);
3957 Double_t sy =0.1, sz =0.1;
3959 Double_t sy1=0.2, sz1=0.2;
3960 Double_t sy2=0.2, sz2=0.2;
3963 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;
3964 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;
3965 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;
3966 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;
3967 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;
3968 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;
3970 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;
3971 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;
3972 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;
3973 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;
3978 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3979 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3980 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3981 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3982 c[13]=f30*sy1*f40+f32*sy2*f42;
3983 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3985 // Int_t row1 = fSectors->GetRowNumber(xyz[2][0]);
3986 AliTPCseed *seed=new AliTPCseed(xyz[2][0], sec[2]*fSectors->GetAlpha()+fSectors->GetAlphaShift(), x, c, 0);
3987 seed->SetLastPoint(row[2]);
3988 seed->SetFirstPoint(row[2]);
3989 for (Int_t i=row[0];i<row[2];i++){
3990 seed->SetClusterIndex(i, track->GetClusterIndex(i));
3998 void AliTPCtrackerMI::FindMultiMC(TObjArray * array, AliESDEvent */*esd*/, Int_t iter)
4001 // find multi tracks - THIS FUNCTION IS ONLY FOR DEBUG PURPOSES
4003 // Use AliTPCReconstructor::StreamLevel()>2 if you want to tune parameters - cuts
4005 // Two reasons to have multiple find tracks
4006 // 1. Curling tracks can be find more than once
4007 // 2. Splitted tracks
4008 // a.) Multiple seeding to increase tracking efficiency - (~ 100% reached)
4009 // b.) Edge effect on the sector boundaries
4012 // Algorithm done in 2 phases - because of CPU consumption
4013 // it is n^2 algorithm - for lead-lead 20000x20000 combination are investigated
4015 // Algorihm for curling tracks sign:
4016 // 1 phase -makes a very rough fast cuts to minimize combinatorics
4017 // a.) opposite sign
4018 // b.) one of the tracks - not pointing to the primary vertex -
4019 // c.) delta tan(theta)
4021 // 2 phase - calculates DCA between tracks - time consument
4026 // General cuts - for splitted tracks and for curling tracks
4028 const Float_t kMaxdPhi = 0.2; // maximal distance in phi
4030 // Curling tracks cuts
4035 Int_t nentries = array->GetEntriesFast();
4036 AliHelix *helixes = new AliHelix[nentries];
4037 Float_t *xm = new Float_t[nentries];
4038 Float_t *dz0 = new Float_t[nentries];
4039 Float_t *dz1 = new Float_t[nentries];
4045 // Find track COG in x direction - point with best defined parameters
4047 for (Int_t i=0;i<nentries;i++){
4048 AliTPCseed* track = (AliTPCseed*)array->At(i);
4049 if (!track) continue;
4050 track->SetCircular(0);
4051 new (&helixes[i]) AliHelix(*track);
4055 track->GetDZ(GetX(),GetY(),GetZ(),GetBz(),dz);
4058 for (Int_t icl=0; icl<160; icl++){
4059 AliTPCclusterMI * cl = track->GetClusterPointer(icl);
4065 if (ncl>0) xm[i]/=Float_t(ncl);
4067 TTreeSRedirector &cstream = *fDebugStreamer;
4069 for (Int_t i0=0;i0<nentries;i0++){
4070 AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
4071 if (!track0) continue;
4072 Float_t xc0 = helixes[i0].GetHelix(6);
4073 Float_t yc0 = helixes[i0].GetHelix(7);
4074 Float_t r0 = helixes[i0].GetHelix(8);
4075 Float_t rc0 = TMath::Sqrt(xc0*xc0+yc0*yc0);
4076 Float_t fi0 = TMath::ATan2(yc0,xc0);
4078 for (Int_t i1=i0+1;i1<nentries;i1++){
4079 AliTPCseed * track1 = (AliTPCseed*)array->At(i1);
4080 if (!track1) continue;
4081 Int_t lab0=track0->GetLabel();
4082 Int_t lab1=track1->GetLabel();
4083 if (TMath::Abs(lab0)!=TMath::Abs(lab1)) continue;
4085 Float_t xc1 = helixes[i1].GetHelix(6);
4086 Float_t yc1 = helixes[i1].GetHelix(7);
4087 Float_t r1 = helixes[i1].GetHelix(8);
4088 Float_t rc1 = TMath::Sqrt(xc1*xc1+yc1*yc1);
4089 Float_t fi1 = TMath::ATan2(yc1,xc1);
4091 Float_t dfi = fi0-fi1;
4094 if (dfi>1.5*TMath::Pi()) dfi-=TMath::Pi(); // take care about edge effect
4095 if (dfi<-1.5*TMath::Pi()) dfi+=TMath::Pi(); //
4096 if (TMath::Abs(dfi)>kMaxdPhi&&helixes[i0].GetHelix(4)*helixes[i1].GetHelix(4)<0){
4098 // if short tracks with undefined sign
4099 fi1 = -TMath::ATan2(yc1,-xc1);
4102 Float_t dtheta = TMath::Abs(track0->GetTgl()-track1->GetTgl())<TMath::Abs(track0->GetTgl()+track1->GetTgl())? track0->GetTgl()-track1->GetTgl():track0->GetTgl()+track1->GetTgl();
4105 // debug stream to tune "fast cuts"
4107 Double_t dist[3]; // distance at X
4108 Double_t mdist[3]={0,0,0}; // mean distance X+-40cm
4109 track0->GetDistance(track1,0.5*(xm[i0]+xm[i1])-40.,dist,AliTracker::GetBz());
4110 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[i]);
4111 track0->GetDistance(track1,0.5*(xm[i0]+xm[i1])+40.,dist,AliTracker::GetBz());
4112 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[i]);
4113 track0->GetDistance(track1,0.5*(xm[i0]+xm[i1]),dist,AliTracker::GetBz());
4114 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[i]);
4115 for (Int_t i=0;i<3;i++) mdist[i]*=0.33333;
4119 for (Int_t icl=0; icl<160; icl++){
4120 AliTPCclusterMI * cl0 = track0->GetClusterPointer(icl);
4121 AliTPCclusterMI * cl1 = track1->GetClusterPointer(icl);
4124 if (cl0==cl1) sums++;
4132 "Tr0.="<<track0<< // seed0
4133 "Tr1.="<<track1<< // seed1
4134 "h0.="<<&helixes[i0]<<
4135 "h1.="<<&helixes[i1]<<
4137 "sum="<<sum<< //the sum of rows with cl in both
4138 "sums="<<sums<< //the sum of shared clusters
4139 "xm0="<<xm[i0]<< // the center of track
4140 "xm1="<<xm[i1]<< // the x center of track
4141 // General cut variables
4142 "dfi="<<dfi<< // distance in fi angle
4143 "dtheta="<<dtheta<< // distance int theta angle
4149 "dist0="<<dist[0]<< //distance x
4150 "dist1="<<dist[1]<< //distance y
4151 "dist2="<<dist[2]<< //distance z
4152 "mdist0="<<mdist[0]<< //distance x
4153 "mdist1="<<mdist[1]<< //distance y
4154 "mdist2="<<mdist[2]<< //distance z
4167 if (AliTPCReconstructor::StreamLevel()>1) {
4168 AliInfo("Time for curling tracks removal DEBUGGING MC");
4174 void AliTPCtrackerMI::FindSplitted(TObjArray * array, AliESDEvent */*esd*/, Int_t iter)
4178 // Two reasons to have multiple find tracks
4179 // 1. Curling tracks can be find more than once
4180 // 2. Splitted tracks
4181 // a.) Multiple seeding to increase tracking efficiency - (~ 100% reached)
4182 // b.) Edge effect on the sector boundaries
4184 // This function tries to find tracks closed in the parametric space
4186 // cut logic if distance is bigger than cut continue - Do Nothing
4187 const Float_t kMaxdTheta = 0.05; // maximal distance in theta
4188 const Float_t kMaxdPhi = 0.05; // maximal deistance in phi
4189 const Float_t kdelta = 40.; //delta r to calculate track distance
4191 // const Float_t kMaxDist0 = 1.; // maximal distance 0
4192 //const Float_t kMaxDist1 = 0.3; // maximal distance 1 - cut if track in separate rows
4195 TCut csec("csec","abs(Tr0.fRelativeSector-Tr1.fRelativeSector)<2");
4196 TCut cdtheta("cdtheta","abs(dtheta)<0.05");
4201 Int_t nentries = array->GetEntriesFast();
4202 AliHelix *helixes = new AliHelix[nentries];
4203 Float_t *xm = new Float_t[nentries];
4209 //Sort tracks according quality
4211 Int_t nseed = array->GetEntriesFast();
4212 Float_t * quality = new Float_t[nseed];
4213 Int_t * indexes = new Int_t[nseed];
4214 for (Int_t i=0; i<nseed; i++) {
4215 AliTPCseed *pt=(AliTPCseed*)array->UncheckedAt(i);
4220 pt->UpdatePoints(); //select first last max dens points
4221 Float_t * points = pt->GetPoints();
4222 if (points[3]<0.8) quality[i] =-1;
4223 quality[i] = (points[2]-points[0])+pt->GetNumberOfClusters();
4224 //prefer high momenta tracks if overlaps
4225 quality[i] *= TMath::Sqrt(TMath::Abs(pt->Pt())+0.5);
4227 TMath::Sort(nseed,quality,indexes);
4231 // Find track COG in x direction - point with best defined parameters
4233 for (Int_t i=0;i<nentries;i++){
4234 AliTPCseed* track = (AliTPCseed*)array->At(i);
4235 if (!track) continue;
4236 track->SetCircular(0);
4237 new (&helixes[i]) AliHelix(*track);
4240 for (Int_t icl=0; icl<160; icl++){
4241 AliTPCclusterMI * cl = track->GetClusterPointer(icl);
4247 if (ncl>0) xm[i]/=Float_t(ncl);
4249 TTreeSRedirector &cstream = *fDebugStreamer;
4251 for (Int_t is0=0;is0<nentries;is0++){
4252 Int_t i0 = indexes[is0];
4253 AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
4254 if (!track0) continue;
4255 if (track0->GetKinkIndexes()[0]!=0) continue;
4256 Float_t xc0 = helixes[i0].GetHelix(6);
4257 Float_t yc0 = helixes[i0].GetHelix(7);
4258 Float_t fi0 = TMath::ATan2(yc0,xc0);
4260 for (Int_t is1=is0+1;is1<nentries;is1++){
4261 Int_t i1 = indexes[is1];
4262 AliTPCseed * track1 = (AliTPCseed*)array->At(i1);
4263 if (!track1) continue;
4265 if (TMath::Abs(track0->GetRelativeSector()-track1->GetRelativeSector())>1) continue;
4266 if (track1->GetKinkIndexes()[0]>0 &&track0->GetKinkIndexes()[0]<0) continue;
4267 if (track1->GetKinkIndexes()[0]!=0) continue;
4269 Float_t dtheta = TMath::Abs(track0->GetTgl()-track1->GetTgl())<TMath::Abs(track0->GetTgl()+track1->GetTgl())? track0->GetTgl()-track1->GetTgl():track0->GetTgl()+track1->GetTgl();
4270 if (TMath::Abs(dtheta)>kMaxdTheta) continue;
4272 Float_t xc1 = helixes[i1].GetHelix(6);
4273 Float_t yc1 = helixes[i1].GetHelix(7);
4274 Float_t fi1 = TMath::ATan2(yc1,xc1);
4276 Float_t dfi = fi0-fi1;
4277 if (dfi>1.5*TMath::Pi()) dfi-=TMath::Pi(); // take care about edge effect
4278 if (dfi<-1.5*TMath::Pi()) dfi+=TMath::Pi(); //
4279 if (TMath::Abs(dfi)>kMaxdPhi&&helixes[i0].GetHelix(4)*helixes[i1].GetHelix(4)<0){
4281 // if short tracks with undefined sign
4282 fi1 = -TMath::ATan2(yc1,-xc1);
4285 if (TMath::Abs(dfi)>kMaxdPhi) continue;
4292 for (Int_t icl=0; icl<160; icl++){
4293 Int_t index0=track0->GetClusterIndex2(icl);
4294 Int_t index1=track1->GetClusterIndex2(icl);
4295 Bool_t used0 = (index0>0 && !(index0&0x8000));
4296 Bool_t used1 = (index1>0 && !(index1&0x8000));
4298 if (used0) sum0++; // used cluster0
4299 if (used1) sum1++; // used clusters1
4300 if (used0&&used1) sum++;
4301 if (index0==index1 && used0 && used1) sums++;
4305 if (sums<10) continue;
4306 if (sum<40) continue;
4307 if (sums/Float_t(TMath::Min(sum0,sum1))<0.5) continue;
4309 Double_t dist[5][4]; // distance at X
4310 Double_t mdist[4]={0,0,0,0}; // mean distance on range +-delta
4314 track0->GetDistance(track1,xm[i0],dist[0],AliTracker::GetBz());
4315 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[0][i]);
4316 track0->GetDistance(track1,xm[i1],dist[1],AliTracker::GetBz());
4317 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[1][i]);
4319 track0->GetDistance(track1,TMath::Min(xm[i1],xm[i0])-kdelta,dist[2],AliTracker::GetBz());
4320 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[2][i]);
4321 track0->GetDistance(track1,TMath::Max(xm[i1],xm[i0])+kdelta,dist[3],AliTracker::GetBz());
4322 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[3][i]);
4324 track0->GetDistance(track1,(xm[i1]+xm[i0])*0.5,dist[4],AliTracker::GetBz());
4325 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[4][i]);
4326 for (Int_t i=0;i<3;i++) mdist[i]*=0.2;
4329 Int_t lab0=track0->GetLabel();
4330 Int_t lab1=track1->GetLabel();
4331 cstream<<"Splitted"<<
4335 "Tr0.="<<track0<< // seed0
4336 "Tr1.="<<track1<< // seed1
4337 "h0.="<<&helixes[i0]<<
4338 "h1.="<<&helixes[i1]<<
4340 "sum="<<sum<< //the sum of rows with cl in both
4341 "sum0="<<sum0<< //the sum of rows with cl in 0 track
4342 "sum1="<<sum1<< //the sum of rows with cl in 1 track
4343 "sums="<<sums<< //the sum of shared clusters
4344 "xm0="<<xm[i0]<< // the center of track
4345 "xm1="<<xm[i1]<< // the x center of track
4346 // General cut variables
4347 "dfi="<<dfi<< // distance in fi angle
4348 "dtheta="<<dtheta<< // distance int theta angle
4351 "dist0="<<dist[4][0]<< //distance x
4352 "dist1="<<dist[4][1]<< //distance y
4353 "dist2="<<dist[4][2]<< //distance z
4354 "mdist0="<<mdist[0]<< //distance x
4355 "mdist1="<<mdist[1]<< //distance y
4356 "mdist2="<<mdist[2]<< //distance z
4359 delete array->RemoveAt(i1);
4364 AliInfo("Time for splitted tracks removal");
4370 void AliTPCtrackerMI::FindCurling(TObjArray * array, AliESDEvent */*esd*/, Int_t iter)
4373 // find Curling tracks
4374 // Use AliTPCReconstructor::StreamLevel()>1 if you want to tune parameters - cuts
4377 // Algorithm done in 2 phases - because of CPU consumption
4378 // it is n^2 algorithm - for lead-lead 20000x20000 combination are investigated
4379 // see detal in MC part what can be used to cut
4383 const Float_t kMaxC = 400; // maximal curvature to of the track
4384 const Float_t kMaxdTheta = 0.15; // maximal distance in theta
4385 const Float_t kMaxdPhi = 0.15; // maximal distance in phi
4386 const Float_t kPtRatio = 0.3; // ratio between pt
4387 const Float_t kMinDCAR = 2.; // distance to the primary vertex in r - see cpipe cut
4390 // Curling tracks cuts
4393 const Float_t kMaxDeltaRMax = 40; // distance in outer radius
4394 const Float_t kMaxDeltaRMin = 5.; // distance in lower radius - see cpipe cut
4395 const Float_t kMinAngle = 2.9; // angle between tracks
4396 const Float_t kMaxDist = 5; // biggest distance
4398 // The cuts can be tuned using the "MC information stored in Multi tree ==> see FindMultiMC
4401 TCut csign("csign","Tr0.fP[4]*Tr1.fP[4]<0"); //opposite sign
4402 TCut cmax("cmax","abs(Tr0.GetC())>1/400");
4403 TCut cda("cda","sqrt(dtheta^2+dfi^2)<0.15");
4404 TCut ccratio("ccratio","abs((Tr0.fP[4]+Tr1.fP[4])/(abs(Tr0.fP[4])+abs(Tr1.fP[4])))<0.3");
4405 TCut cpipe("cpipe", "min(abs(r0-rc0),abs(r1-rc1))>5");
4407 TCut cdrmax("cdrmax","abs(abs(rc0+r0)-abs(rc1+r1))<40")
4408 TCut cdrmin("cdrmin","abs(abs(rc0+r0)-abs(rc1+r1))<10")
4410 Multi->Draw("dfi","iter==0"+csign+cmax+cda+ccratio); ~94% of curling tracks fulfill
4411 Multi->Draw("min(abs(r0-rc0),abs(r1-rc1))","iter==0&&abs(lab1)==abs(lab0)"+csign+cmax+cda+ccratio+cpipe+cdrmin+cdrmax); //80%
4413 Curling2->Draw("dfi","iter==0&&abs(lab0)==abs(lab1)"+csign+cmax+cdtheta+cdfi+ccratio)
4419 Int_t nentries = array->GetEntriesFast();
4420 AliHelix *helixes = new AliHelix[nentries];
4421 for (Int_t i=0;i<nentries;i++){
4422 AliTPCseed* track = (AliTPCseed*)array->At(i);
4423 if (!track) continue;
4424 track->SetCircular(0);
4425 new (&helixes[i]) AliHelix(*track);
4431 Double_t phase[2][2],radius[2];
4435 TTreeSRedirector &cstream = *fDebugStreamer;
4437 for (Int_t i0=0;i0<nentries;i0++){
4438 AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
4439 if (!track0) continue;
4440 if (TMath::Abs(track0->GetC())<1/kMaxC) continue;
4441 Float_t xc0 = helixes[i0].GetHelix(6);
4442 Float_t yc0 = helixes[i0].GetHelix(7);
4443 Float_t r0 = helixes[i0].GetHelix(8);
4444 Float_t rc0 = TMath::Sqrt(xc0*xc0+yc0*yc0);
4445 Float_t fi0 = TMath::ATan2(yc0,xc0);
4447 for (Int_t i1=i0+1;i1<nentries;i1++){
4448 AliTPCseed * track1 = (AliTPCseed*)array->At(i1);
4449 if (!track1) continue;
4450 if (TMath::Abs(track1->GetC())<1/kMaxC) continue;
4451 Float_t xc1 = helixes[i1].GetHelix(6);
4452 Float_t yc1 = helixes[i1].GetHelix(7);
4453 Float_t r1 = helixes[i1].GetHelix(8);
4454 Float_t rc1 = TMath::Sqrt(xc1*xc1+yc1*yc1);
4455 Float_t fi1 = TMath::ATan2(yc1,xc1);
4457 Float_t dfi = fi0-fi1;
4460 if (dfi>1.5*TMath::Pi()) dfi-=TMath::Pi(); // take care about edge effect
4461 if (dfi<-1.5*TMath::Pi()) dfi+=TMath::Pi(); //
4462 Float_t dtheta = TMath::Abs(track0->GetTgl()-track1->GetTgl())<TMath::Abs(track0->GetTgl()+track1->GetTgl())? track0->GetTgl()-track1->GetTgl():track0->GetTgl()+track1->GetTgl();
4466 if (track0->GetBConstrain()&&track1->GetBConstrain()) continue; // not constrained
4467 if (track1->GetSigned1Pt()*track0->GetSigned1Pt()>0) continue; // not the same sign
4468 if ( TMath::Abs(track1->GetTgl()+track0->GetTgl())>kMaxdTheta) continue; //distance in the Theta
4469 if ( TMath::Abs(dfi)>kMaxdPhi) continue; //distance in phi
4470 if ( TMath::Sqrt(dfi*dfi+dtheta*dtheta)>kMaxdPhi) continue; //common angular offset
4472 Float_t pt0 = track0->GetSignedPt();
4473 Float_t pt1 = track1->GetSignedPt();
4474 if ((TMath::Abs(pt0+pt1)/(TMath::Abs(pt0)+TMath::Abs(pt1)))>kPtRatio) continue;
4475 if ((iter==1) && TMath::Abs(TMath::Abs(rc0+r0)-TMath::Abs(rc1+r1))>kMaxDeltaRMax) continue;
4476 if ((iter!=1) &&TMath::Abs(TMath::Abs(rc0-r0)-TMath::Abs(rc1-r1))>kMaxDeltaRMin) continue;
4477 if (TMath::Min(TMath::Abs(rc0-r0),TMath::Abs(rc1-r1))<kMinDCAR) continue;
4480 // Now find closest approach
4484 Int_t npoints = helixes[i0].GetRPHIintersections(helixes[i1], phase, radius,10);
4485 if (npoints==0) continue;
4486 helixes[i0].GetClosestPhases(helixes[i1], phase);
4490 Double_t hangles[3];
4491 helixes[i0].Evaluate(phase[0][0],xyz0);
4492 helixes[i1].Evaluate(phase[0][1],xyz1);
4494 helixes[i0].GetAngle(phase[0][0],helixes[i1],phase[0][1],hangles);
4495 Double_t deltah[2],deltabest;
4496 if (TMath::Abs(hangles[2])<kMinAngle) continue;
4500 helixes[i0].ParabolicDCA(helixes[i1],phase[0][0],phase[0][1],radius[0],deltah[0],2);
4502 helixes[i0].ParabolicDCA(helixes[i1],phase[1][0],phase[1][1],radius[1],deltah[1],2);
4503 if (deltah[1]<deltah[0]) ibest=1;
4505 deltabest = TMath::Sqrt(deltah[ibest]);
4506 helixes[i0].Evaluate(phase[ibest][0],xyz0);
4507 helixes[i1].Evaluate(phase[ibest][1],xyz1);
4508 helixes[i0].GetAngle(phase[ibest][0],helixes[i1],phase[ibest][1],hangles);
4509 Double_t radiusbest = TMath::Sqrt(radius[ibest]);
4511 if (deltabest>kMaxDist) continue;
4512 // if (mindcar+mindcaz<40 && (TMath::Abs(hangles[2])<kMinAngle ||deltabest>3)) continue;
4513 Bool_t sign =kFALSE;
4514 if (hangles[2]>kMinAngle) sign =kTRUE;
4517 // circular[i0] = kTRUE;
4518 // circular[i1] = kTRUE;
4519 if (track0->OneOverPt()<track1->OneOverPt()){
4520 track0->SetCircular(track0->GetCircular()+1);
4521 track1->SetCircular(track1->GetCircular()+2);
4524 track1->SetCircular(track1->GetCircular()+1);
4525 track0->SetCircular(track0->GetCircular()+2);
4528 if (AliTPCReconstructor::StreamLevel()>1){
4530 //debug stream to tune "fine" cuts
4531 Int_t lab0=track0->GetLabel();
4532 Int_t lab1=track1->GetLabel();
4533 cstream<<"Curling2"<<
4549 "npoints="<<npoints<<
4550 "hangles0="<<hangles[0]<<
4551 "hangles1="<<hangles[1]<<
4552 "hangles2="<<hangles[2]<<
4555 "radius="<<radiusbest<<
4556 "deltabest="<<deltabest<<
4557 "phase0="<<phase[ibest][0]<<
4558 "phase1="<<phase[ibest][1]<<
4566 if (AliTPCReconstructor::StreamLevel()>1) {
4567 AliInfo("Time for curling tracks removal");
4576 void AliTPCtrackerMI::FindKinks(TObjArray * array, AliESDEvent *esd)
4583 TObjArray *kinks= new TObjArray(10000);
4584 // TObjArray *v0s= new TObjArray(10000);
4585 Int_t nentries = array->GetEntriesFast();
4586 AliHelix *helixes = new AliHelix[nentries];
4587 Int_t *sign = new Int_t[nentries];
4588 Int_t *nclusters = new Int_t[nentries];
4589 Float_t *alpha = new Float_t[nentries];
4590 AliKink *kink = new AliKink();
4591 Int_t * usage = new Int_t[nentries];
4592 Float_t *zm = new Float_t[nentries];
4593 Float_t *z0 = new Float_t[nentries];
4594 Float_t *fim = new Float_t[nentries];
4595 Float_t *shared = new Float_t[nentries];
4596 Bool_t *circular = new Bool_t[nentries];
4597 Float_t *dca = new Float_t[nentries];
4598 //const AliESDVertex * primvertex = esd->GetVertex();
4600 // nentries = array->GetEntriesFast();
4605 for (Int_t i=0;i<nentries;i++){
4608 AliTPCseed* track = (AliTPCseed*)array->At(i);
4609 if (!track) continue;
4610 track->SetCircular(0);
4612 track->UpdatePoints();
4613 if (( track->GetPoints()[2]- track->GetPoints()[0])>5 && track->GetPoints()[3]>0.8){
4615 nclusters[i]=track->GetNumberOfClusters();
4616 alpha[i] = track->GetAlpha();
4617 new (&helixes[i]) AliHelix(*track);
4619 helixes[i].Evaluate(0,xyz);
4620 sign[i] = (track->GetC()>0) ? -1:1;
4623 if (track->GetProlongation(x,y,z)){
4625 fim[i] = alpha[i]+TMath::ATan2(y,x);
4628 zm[i] = track->GetZ();
4632 circular[i]= kFALSE;
4633 if (track->GetProlongation(0,y,z)) z0[i] = z;
4634 dca[i] = track->GetD(0,0);
4640 Int_t ncandidates =0;
4643 Double_t phase[2][2],radius[2];
4646 // Find circling track
4647 TTreeSRedirector &cstream = *fDebugStreamer;
4649 for (Int_t i0=0;i0<nentries;i0++){
4650 AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
4651 if (!track0) continue;
4652 if (track0->GetNumberOfClusters()<40) continue;
4653 if (TMath::Abs(1./track0->GetC())>200) continue;
4654 for (Int_t i1=i0+1;i1<nentries;i1++){
4655 AliTPCseed * track1 = (AliTPCseed*)array->At(i1);
4656 if (!track1) continue;
4657 if (track1->GetNumberOfClusters()<40) continue;
4658 if ( TMath::Abs(track1->GetTgl()+track0->GetTgl())>0.1) continue;
4659 if (track0->GetBConstrain()&&track1->GetBConstrain()) continue;
4660 if (TMath::Abs(1./track1->GetC())>200) continue;
4661 if (track1->GetSigned1Pt()*track0->GetSigned1Pt()>0) continue;
4662 if (track1->GetTgl()*track0->GetTgl()>0) continue;
4663 if (TMath::Max(TMath::Abs(1./track0->GetC()),TMath::Abs(1./track1->GetC()))>190) continue;
4664 if (track0->GetBConstrain()&&track1->OneOverPt()<track0->OneOverPt()) continue; //returning - lower momenta
4665 if (track1->GetBConstrain()&&track0->OneOverPt()<track1->OneOverPt()) continue; //returning - lower momenta
4667 Float_t mindcar = TMath::Min(TMath::Abs(dca[i0]),TMath::Abs(dca[i1]));
4668 if (mindcar<5) continue;
4669 Float_t mindcaz = TMath::Min(TMath::Abs(z0[i0]-GetZ()),TMath::Abs(z0[i1]-GetZ()));
4670 if (mindcaz<5) continue;
4671 if (mindcar+mindcaz<20) continue;
4674 Float_t xc0 = helixes[i0].GetHelix(6);
4675 Float_t yc0 = helixes[i0].GetHelix(7);
4676 Float_t r0 = helixes[i0].GetHelix(8);
4677 Float_t xc1 = helixes[i1].GetHelix(6);
4678 Float_t yc1 = helixes[i1].GetHelix(7);
4679 Float_t r1 = helixes[i1].GetHelix(8);
4681 Float_t rmean = (r0+r1)*0.5;
4682 Float_t delta =TMath::Sqrt((xc1-xc0)*(xc1-xc0)+(yc1-yc0)*(yc1-yc0));
4683 //if (delta>30) continue;
4684 if (delta>rmean*0.25) continue;
4685 if (TMath::Abs(r0-r1)/rmean>0.3) continue;
4687 Int_t npoints = helixes[i0].GetRPHIintersections(helixes[i1], phase, radius,10);
4688 if (npoints==0) continue;
4689 helixes[i0].GetClosestPhases(helixes[i1], phase);
4693 Double_t hangles[3];
4694 helixes[i0].Evaluate(phase[0][0],xyz0);
4695 helixes[i1].Evaluate(phase[0][1],xyz1);
4697 helixes[i0].GetAngle(phase[0][0],helixes[i1],phase[0][1],hangles);
4698 Double_t deltah[2],deltabest;
4699 if (hangles[2]<2.8) continue;
4702 helixes[i0].ParabolicDCA(helixes[i1],phase[0][0],phase[0][1],radius[0],deltah[0],2);
4704 helixes[i0].ParabolicDCA(helixes[i1],phase[1][0],phase[1][1],radius[1],deltah[1],2);
4705 if (deltah[1]<deltah[0]) ibest=1;
4707 deltabest = TMath::Sqrt(deltah[ibest]);
4708 helixes[i0].Evaluate(phase[ibest][0],xyz0);
4709 helixes[i1].Evaluate(phase[ibest][1],xyz1);
4710 helixes[i0].GetAngle(phase[ibest][0],helixes[i1],phase[ibest][1],hangles);
4711 Double_t radiusbest = TMath::Sqrt(radius[ibest]);
4713 if (deltabest>6) continue;
4714 if (mindcar+mindcaz<40 && (hangles[2]<3.12||deltabest>3)) continue;
4715 Bool_t sign =kFALSE;
4716 if (hangles[2]>3.06) sign =kTRUE;
4719 circular[i0] = kTRUE;
4720 circular[i1] = kTRUE;
4721 if (track0->OneOverPt()<track1->OneOverPt()){
4722 track0->SetCircular(track0->GetCircular()+1);
4723 track1->SetCircular(track1->GetCircular()+2);
4726 track1->SetCircular(track1->GetCircular()+1);
4727 track0->SetCircular(track0->GetCircular()+2);
4730 if (sign&&AliTPCReconstructor::StreamLevel()>1){
4732 Int_t lab0=track0->GetLabel();
4733 Int_t lab1=track1->GetLabel();
4734 cstream<<"Curling"<<
4741 "mindcar="<<mindcar<<
4742 "mindcaz="<<mindcaz<<
4745 "npoints="<<npoints<<
4746 "hangles0="<<hangles[0]<<
4747 "hangles2="<<hangles[2]<<
4752 "radius="<<radiusbest<<
4753 "deltabest="<<deltabest<<
4754 "phase0="<<phase[ibest][0]<<
4755 "phase1="<<phase[ibest][1]<<
4765 for (Int_t i =0;i<nentries;i++){
4766 if (sign[i]==0) continue;
4767 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
4774 Double_t cradius0 = 40*40;
4775 Double_t cradius1 = 270*270;
4778 Double_t cdist3=0.55;
4779 for (Int_t j =i+1;j<nentries;j++){
4781 if (sign[j]*sign[i]<1) continue;
4782 if ( (nclusters[i]+nclusters[j])>200) continue;
4783 if ( (nclusters[i]+nclusters[j])<80) continue;
4784 if ( TMath::Abs(zm[i]-zm[j])>60.) continue;
4785 if ( TMath::Abs(fim[i]-fim[j])>0.6 && TMath::Abs(fim[i]-fim[j])<5.7 ) continue;
4786 //AliTPCseed * track1 = (AliTPCseed*)array->At(j); Double_t phase[2][2],radius[2];
4787 Int_t npoints = helixes[i].GetRPHIintersections(helixes[j], phase, radius,20);
4788 if (npoints<1) continue;
4791 if (radius[0]<cradius0||radius[0]>cradius1) continue;
4794 if ( (radius[0]<cradius0||radius[0]>cradius1) && (radius[1]<cradius0||radius[1]>cradius1) ) continue;
4797 Double_t delta1=10000,delta2=10000;
4798 // cuts on the intersection radius
4799 helixes[i].LinearDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta1);
4800 if (radius[0]<20&&delta1<1) continue; //intersection at vertex
4801 if (radius[0]<10&&delta1<3) continue; //intersection at vertex
4803 helixes[i].LinearDCA(helixes[j],phase[1][0],phase[1][1],radius[1],delta2);
4804 if (radius[1]<20&&delta2<1) continue; //intersection at vertex
4805 if (radius[1]<10&&delta2<3) continue; //intersection at vertex
4808 Double_t distance1 = TMath::Min(delta1,delta2);
4809 if (distance1>cdist1) continue; // cut on DCA linear approximation
4811 npoints = helixes[i].GetRPHIintersections(helixes[j], phase, radius,20);
4812 helixes[i].ParabolicDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta1);
4813 if (radius[0]<20&&delta1<1) continue; //intersection at vertex
4814 if (radius[0]<10&&delta1<3) continue; //intersection at vertex
4817 helixes[i].ParabolicDCA(helixes[j],phase[1][0],phase[1][1],radius[1],delta2);
4818 if (radius[1]<20&&delta2<1) continue; //intersection at vertex
4819 if (radius[1]<10&&delta2<3) continue; //intersection at vertex
4821 distance1 = TMath::Min(delta1,delta2);
4824 rkink = TMath::Sqrt(radius[0]);
4827 rkink = TMath::Sqrt(radius[1]);
4829 if (distance1>cdist2) continue;
4832 AliTPCseed * track1 = (AliTPCseed*)array->At(j);
4835 Int_t row0 = GetRowNumber(rkink);
4836 if (row0<10) continue;
4837 if (row0>150) continue;
4840 Float_t dens00=-1,dens01=-1;
4841 Float_t dens10=-1,dens11=-1;
4843 Int_t found,foundable,shared;
4844 track0->GetClusterStatistic(0,row0-5, found, foundable,shared,kFALSE);
4845 if (foundable>5) dens00 = Float_t(found)/Float_t(foundable);
4846 track0->GetClusterStatistic(row0+5,155, found, foundable,shared,kFALSE);
4847 if (foundable>5) dens01 = Float_t(found)/Float_t(foundable);
4849 track1->GetClusterStatistic(0,row0-5, found, foundable,shared,kFALSE);
4850 if (foundable>10) dens10 = Float_t(found)/Float_t(foundable);
4851 track1->GetClusterStatistic(row0+5,155, found, foundable,shared,kFALSE);
4852 if (foundable>10) dens11 = Float_t(found)/Float_t(foundable);
4854 if (dens00<dens10 && dens01<dens11) continue;
4855 if (dens00>dens10 && dens01>dens11) continue;
4856 if (TMath::Max(dens00,dens10)<0.1) continue;
4857 if (TMath::Max(dens01,dens11)<0.3) continue;
4859 if (TMath::Min(dens00,dens10)>0.6) continue;
4860 if (TMath::Min(dens01,dens11)>0.6) continue;
4863 AliTPCseed * ktrack0, *ktrack1;
4872 if (TMath::Abs(ktrack0->GetC())>5) continue; // cut on the curvature for mother particle
4873 AliExternalTrackParam paramm(*ktrack0);
4874 AliExternalTrackParam paramd(*ktrack1);
4875 if (row0>60&&ktrack1->GetReference().GetX()>90.)new (¶md) AliExternalTrackParam(ktrack1->GetReference());
4878 kink->SetMother(paramm);
4879 kink->SetDaughter(paramd);
4882 Float_t x[3] = { kink->GetPosition()[0],kink->GetPosition()[1],kink->GetPosition()[2]};
4884 fParam->Transform0to1(x,index);
4885 fParam->Transform1to2(x,index);
4886 row0 = GetRowNumber(x[0]);
4888 if (kink->GetR()<100) continue;
4889 if (kink->GetR()>240) continue;
4890 if (kink->GetPosition()[2]/kink->GetR()>AliTPCReconstructor::GetCtgRange()) continue; //out of fiducial volume
4891 if (kink->GetDistance()>cdist3) continue;
4892 Float_t dird = kink->GetDaughterP()[0]*kink->GetPosition()[0]+kink->GetDaughterP()[1]*kink->GetPosition()[1]; // rough direction estimate
4893 if (dird<0) continue;
4895 Float_t dirm = kink->GetMotherP()[0]*kink->GetPosition()[0]+kink->GetMotherP()[1]*kink->GetPosition()[1]; // rough direction estimate
4896 if (dirm<0) continue;
4897 Float_t mpt = TMath::Sqrt(kink->GetMotherP()[0]*kink->GetMotherP()[0]+kink->GetMotherP()[1]*kink->GetMotherP()[1]);
4898 if (mpt<0.2) continue;
4901 //for high momenta momentum not defined well in first iteration
4902 Double_t qt = TMath::Sin(kink->GetAngle(2))*ktrack1->GetP();
4903 if (qt>0.35) continue;
4906 kink->SetLabel(CookLabel(ktrack0,0.4,0,row0),0);
4907 kink->SetLabel(CookLabel(ktrack1,0.4,row0,160),1);
4909 kink->SetTPCDensity(dens00,0,0);
4910 kink->SetTPCDensity(dens01,0,1);
4911 kink->SetTPCDensity(dens10,1,0);
4912 kink->SetTPCDensity(dens11,1,1);
4913 kink->SetIndex(i,0);
4914 kink->SetIndex(j,1);
4917 kink->SetTPCDensity(dens10,0,0);
4918 kink->SetTPCDensity(dens11,0,1);
4919 kink->SetTPCDensity(dens00,1,0);
4920 kink->SetTPCDensity(dens01,1,1);
4921 kink->SetIndex(j,0);
4922 kink->SetIndex(i,1);
4925 if (mpt<1||kink->GetAngle(2)>0.1){
4926 // angle and densities not defined yet
4927 if (kink->GetTPCDensityFactor()<0.8) continue;
4928 if ((2-kink->GetTPCDensityFactor())*kink->GetDistance() >0.25) continue;
4929 if (kink->GetAngle(2)*ktrack0->GetP()<0.003) continue; //too small angle
4930 if (kink->GetAngle(2)>0.2&&kink->GetTPCDensityFactor()<1.15) continue;
4931 if (kink->GetAngle(2)>0.2&&kink->GetTPCDensity(0,1)>0.05) continue;
4933 Float_t criticalangle = track0->GetSigmaSnp2()+track0->GetSigmaTgl2();
4934 criticalangle+= track1->GetSigmaSnp2()+track1->GetSigmaTgl2();
4935 criticalangle= 3*TMath::Sqrt(criticalangle);
4936 if (criticalangle>0.02) criticalangle=0.02;
4937 if (kink->GetAngle(2)<criticalangle) continue;
4940 Int_t drow = Int_t(2.+0.5/(0.05+kink->GetAngle(2))); // overlap region defined
4941 Float_t shapesum =0;
4943 for ( Int_t row = row0-drow; row<row0+drow;row++){
4944 if (row<0) continue;
4945 if (row>155) continue;
4946 if (ktrack0->GetClusterPointer(row)){
4947 AliTPCTrackerPoint *point =ktrack0->GetTrackPoint(row);
4948 shapesum+=point->GetSigmaY()+point->GetSigmaZ();
4951 if (ktrack1->GetClusterPointer(row)){
4952 AliTPCTrackerPoint *point =ktrack1->GetTrackPoint(row);
4953 shapesum+=point->GetSigmaY()+point->GetSigmaZ();
4958 kink->SetShapeFactor(-1.);
4961 kink->SetShapeFactor(shapesum/sum);
4963 // esd->AddKink(kink);
4964 kinks->AddLast(kink);
4970 // sort the kinks according quality - and refit them towards vertex
4972 Int_t nkinks = kinks->GetEntriesFast();
4973 Float_t *quality = new Float_t[nkinks];
4974 Int_t *indexes = new Int_t[nkinks];
4975 AliTPCseed *mothers = new AliTPCseed[nkinks];
4976 AliTPCseed *daughters = new AliTPCseed[nkinks];
4979 for (Int_t i=0;i<nkinks;i++){
4981 AliKink *kink = (AliKink*)kinks->At(i);
4983 // refit kinks towards vertex
4985 Int_t index0 = kink->GetIndex(0);
4986 Int_t index1 = kink->GetIndex(1);
4987 AliTPCseed * ktrack0 = (AliTPCseed*)array->At(index0);
4988 AliTPCseed * ktrack1 = (AliTPCseed*)array->At(index1);
4990 Int_t sumn=ktrack0->GetNumberOfClusters()+ktrack1->GetNumberOfClusters();
4992 // Refit Kink under if too small angle
4994 if (kink->GetAngle(2)<0.05){
4995 kink->SetTPCRow0(GetRowNumber(kink->GetR()));
4996 Int_t row0 = kink->GetTPCRow0();
4997 Int_t drow = Int_t(2.+0.5/(0.05+kink->GetAngle(2)));
5000 Int_t last = row0-drow;
5001 if (last<40) last=40;
5002 if (last<ktrack0->GetFirstPoint()+25) last = ktrack0->GetFirstPoint()+25;
5003 AliTPCseed* seed0 = ReSeed(ktrack0,last,kFALSE);
5006 Int_t first = row0+drow;
5007 if (first>130) first=130;
5008 if (first>ktrack1->GetLastPoint()-25) first = TMath::Max(ktrack1->GetLastPoint()-25,30);
5009 AliTPCseed* seed1 = ReSeed(ktrack1,first,kTRUE);
5011 if (seed0 && seed1){
5012 kink->SetStatus(1,8);
5013 if (RefitKink(*seed0,*seed1,*kink)) kink->SetStatus(1,9);
5014 row0 = GetRowNumber(kink->GetR());
5015 sumn = seed0->GetNumberOfClusters()+seed1->GetNumberOfClusters();
5016 mothers[i] = *seed0;
5017 daughters[i] = *seed1;
5020 delete kinks->RemoveAt(i);
5021 if (seed0) delete seed0;
5022 if (seed1) delete seed1;
5025 if (kink->GetDistance()>0.5 || kink->GetR()<110 || kink->GetR()>240) {
5026 delete kinks->RemoveAt(i);
5027 if (seed0) delete seed0;
5028 if (seed1) delete seed1;
5036 if (kink) quality[i] = 160*((0.1+kink->GetDistance())*(2.-kink->GetTPCDensityFactor()))/(sumn+40.); //the longest -clossest will win
5038 TMath::Sort(nkinks,quality,indexes,kFALSE);
5040 //remove double find kinks
5042 for (Int_t ikink0=1;ikink0<nkinks;ikink0++){
5043 AliKink * kink0 = (AliKink*) kinks->At(indexes[ikink0]);
5044 if (!kink0) continue;
5046 for (Int_t ikink1=0;ikink1<ikink0;ikink1++){
5047 if (!kink0) continue;
5048 AliKink * kink1 = (AliKink*) kinks->At(indexes[ikink1]);
5049 if (!kink1) continue;
5050 // if not close kink continue
5051 if (TMath::Abs(kink1->GetPosition()[2]-kink0->GetPosition()[2])>10) continue;
5052 if (TMath::Abs(kink1->GetPosition()[1]-kink0->GetPosition()[1])>10) continue;
5053 if (TMath::Abs(kink1->GetPosition()[0]-kink0->GetPosition()[0])>10) continue;
5055 AliTPCseed &mother0 = mothers[indexes[ikink0]];
5056 AliTPCseed &daughter0 = daughters[indexes[ikink0]];
5057 AliTPCseed &mother1 = mothers[indexes[ikink1]];
5058 AliTPCseed &daughter1 = daughters[indexes[ikink1]];
5059 Int_t row0 = (kink0->GetTPCRow0()+kink1->GetTPCRow0())/2;
5068 for (Int_t i=0;i<row0;i++){
5069 if (mother0.GetClusterIndex(i)>0 && mother1.GetClusterIndex(i)>0){
5072 if (mother0.GetClusterIndex(i)==mother1.GetClusterIndex(i)){
5079 for (Int_t i=row0;i<158;i++){
5080 if (daughter0.GetClusterIndex(i)>0 && daughter0.GetClusterIndex(i)>0){
5083 if (mother0.GetClusterIndex(i)==mother1.GetClusterIndex(i)){
5089 Float_t ratio = Float_t(same+1)/Float_t(both+1);
5090 Float_t ratiom = Float_t(samem+1)/Float_t(bothm+1);
5091 Float_t ratiod = Float_t(samed+1)/Float_t(bothd+1);
5092 if (ratio>0.3 && ratiom>0.5 &&ratiod>0.5) {
5093 Int_t sum0 = mother0.GetNumberOfClusters()+daughter0.GetNumberOfClusters();
5094 Int_t sum1 = mother1.GetNumberOfClusters()+daughter1.GetNumberOfClusters();
5096 shared[kink0->GetIndex(0)]= kTRUE;
5097 shared[kink0->GetIndex(1)]= kTRUE;
5098 delete kinks->RemoveAt(indexes[ikink0]);
5101 shared[kink1->GetIndex(0)]= kTRUE;
5102 shared[kink1->GetIndex(1)]= kTRUE;
5103 delete kinks->RemoveAt(indexes[ikink1]);
5110 for (Int_t i=0;i<nkinks;i++){
5111 AliKink * kink = (AliKink*) kinks->At(indexes[i]);
5112 if (!kink) continue;
5113 kink->SetTPCRow0(GetRowNumber(kink->GetR()));
5114 Int_t index0 = kink->GetIndex(0);
5115 Int_t index1 = kink->GetIndex(1);
5116 if (circular[index0]||circular[index1]&&kink->GetDistance()>0.2) continue;
5117 kink->SetMultiple(usage[index0],0);
5118 kink->SetMultiple(usage[index1],1);
5119 if (kink->GetMultiple()[0]+kink->GetMultiple()[1]>2) continue;
5120 if (kink->GetMultiple()[0]+kink->GetMultiple()[1]>0 && quality[indexes[i]]>0.2) continue;
5121 if (kink->GetMultiple()[0]+kink->GetMultiple()[1]>0 && kink->GetDistance()>0.2) continue;
5122 if (circular[index0]||circular[index1]&&kink->GetDistance()>0.1) continue;
5124 AliTPCseed * ktrack0 = (AliTPCseed*)array->At(index0);
5125 AliTPCseed * ktrack1 = (AliTPCseed*)array->At(index1);
5126 if (!ktrack0 || !ktrack1) continue;
5127 Int_t index = esd->AddKink(kink);
5130 if ( ktrack0->GetKinkIndex(0)==0 && ktrack1->GetKinkIndex(0)==0) { //best kink
5131 if (mothers[indexes[i]].GetNumberOfClusters()>20 && daughters[indexes[i]].GetNumberOfClusters()>20 && (mothers[indexes[i]].GetNumberOfClusters()+daughters[indexes[i]].GetNumberOfClusters())>100){
5132 *ktrack0 = mothers[indexes[i]];
5133 *ktrack1 = daughters[indexes[i]];
5137 ktrack0->SetKinkIndex(usage[index0],-(index+1));
5138 ktrack1->SetKinkIndex(usage[index1], (index+1));
5143 // Remove tracks corresponding to shared kink's
5145 for (Int_t i=0;i<nentries;i++){
5146 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5147 if (!track0) continue;
5148 if (track0->GetKinkIndex(0)!=0) continue;
5149 if (shared[i]) delete array->RemoveAt(i);
5154 RemoveUsed2(array,0.5,0.4,30);
5156 for (Int_t i=0;i<nentries;i++){
5157 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5158 if (!track0) continue;
5159 track0->CookdEdx(0.02,0.6);
5163 for (Int_t i=0;i<nentries;i++){
5164 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5165 if (!track0) continue;
5166 if (track0->Pt()<1.4) continue;
5167 //remove double high momenta tracks - overlapped with kink candidates
5170 for (Int_t icl=track0->GetFirstPoint();icl<track0->GetLastPoint(); icl++){
5171 if (track0->GetClusterPointer(icl)!=0){
5173 if (track0->GetClusterPointer(icl)->IsUsed(10)) shared++;
5176 if (Float_t(shared+1)/Float_t(all+1)>0.5) {
5177 delete array->RemoveAt(i);
5181 if (track0->GetKinkIndex(0)!=0) continue;
5182 if (track0->GetNumberOfClusters()<80) continue;
5184 AliTPCseed *pmother = new AliTPCseed();
5185 AliTPCseed *pdaughter = new AliTPCseed();
5186 AliKink *pkink = new AliKink;
5188 AliTPCseed & mother = *pmother;
5189 AliTPCseed & daughter = *pdaughter;
5190 AliKink & kink = *pkink;
5191 if (CheckKinkPoint(track0,mother,daughter, kink)){
5192 if (mother.GetNumberOfClusters()<30||daughter.GetNumberOfClusters()<20) {
5196 continue; //too short tracks
5198 if (mother.Pt()<1.4) {
5204 Int_t row0= kink.GetTPCRow0();
5205 if (kink.GetDistance()>0.5 || kink.GetR()<110. || kink.GetR()>240.) {
5212 Int_t index = esd->AddKink(&kink);
5213 mother.SetKinkIndex(0,-(index+1));
5214 daughter.SetKinkIndex(0,index+1);
5215 if (mother.GetNumberOfClusters()>50) {
5216 delete array->RemoveAt(i);
5217 array->AddAt(new AliTPCseed(mother),i);
5220 array->AddLast(new AliTPCseed(mother));
5222 array->AddLast(new AliTPCseed(daughter));
5223 for (Int_t icl=0;icl<row0;icl++) {
5224 if (mother.GetClusterPointer(icl)) mother.GetClusterPointer(icl)->Use(20);
5227 for (Int_t icl=row0;icl<158;icl++) {
5228 if (daughter.GetClusterPointer(icl)) daughter.GetClusterPointer(icl)->Use(20);
5237 delete [] daughters;
5259 printf("Ncandidates=\t%d\t%d\t%d\t%d\n",esd->GetNumberOfKinks(),ncandidates,ntracks,nall);
5263 void AliTPCtrackerMI::FindV0s(TObjArray * array, AliESDEvent *esd)
5269 TObjArray *tpcv0s = new TObjArray(100000);
5270 Int_t nentries = array->GetEntriesFast();
5271 AliHelix *helixes = new AliHelix[nentries];
5272 Int_t *sign = new Int_t[nentries];
5273 Float_t *alpha = new Float_t[nentries];
5274 Float_t *z0 = new Float_t[nentries];
5275 Float_t *dca = new Float_t[nentries];
5276 Float_t *sdcar = new Float_t[nentries];
5277 Float_t *cdcar = new Float_t[nentries];
5278 Float_t *pulldcar = new Float_t[nentries];
5279 Float_t *pulldcaz = new Float_t[nentries];
5280 Float_t *pulldca = new Float_t[nentries];
5281 Bool_t *isPrim = new Bool_t[nentries];
5282 const AliESDVertex * primvertex = esd->GetVertex();
5283 Double_t zvertex = primvertex->GetZv();
5285 // nentries = array->GetEntriesFast();
5287 for (Int_t i=0;i<nentries;i++){
5290 AliTPCseed* track = (AliTPCseed*)array->At(i);
5291 if (!track) continue;
5292 track->GetV0Indexes()[0] = 0; //rest v0 indexes
5293 track->GetV0Indexes()[1] = 0; //rest v0 indexes
5294 track->GetV0Indexes()[2] = 0; //rest v0 indexes
5296 alpha[i] = track->GetAlpha();
5297 new (&helixes[i]) AliHelix(*track);
5299 helixes[i].Evaluate(0,xyz);
5300 sign[i] = (track->GetC()>0) ? -1:1;
5304 if (track->GetProlongation(0,y,z)) z0[i] = z;
5305 dca[i] = track->GetD(0,0);
5307 // dca error parrameterezation + pulls
5309 sdcar[i] = TMath::Sqrt(0.150*0.150+(100*track->GetC())*(100*track->GetC()));
5310 if (TMath::Abs(track->GetTgl())>1) sdcar[i]*=2.5;
5311 cdcar[i] = TMath::Exp((TMath::Abs(track->GetC())-0.0106)*525.3);
5312 pulldcar[i] = (dca[i]-cdcar[i])/sdcar[i];
5313 pulldcaz[i] = (z0[i]-zvertex)/sdcar[i];
5314 pulldca[i] = TMath::Sqrt(pulldcar[i]*pulldcar[i]+pulldcaz[i]*pulldcaz[i]);
5315 if (track->TPCrPID(1)+track->TPCrPID(2)+track->TPCrPID(3)>0.5) {
5316 if (pulldca[i]<3.) isPrim[i]=kTRUE; //pion, muon and Kaon 3 sigma cut
5318 if (track->TPCrPID(4)>0.5) {
5319 if (pulldca[i]<0.5) isPrim[i]=kTRUE; //proton 0.5 sigma cut
5321 if (track->TPCrPID(0)>0.4) {
5322 isPrim[i]=kFALSE; //electron no sigma cut
5329 Int_t ncandidates =0;
5332 Double_t phase[2][2],radius[2];
5338 TTreeSRedirector &cstream = *fDebugStreamer;
5339 Float_t fprimvertex[3]={GetX(),GetY(),GetZ()};
5341 Double_t cradius0 = 10*10;
5342 Double_t cradius1 = 200*200;
5345 Double_t cpointAngle = 0.95;
5347 Double_t delta[2]={10000,10000};
5348 for (Int_t i =0;i<nentries;i++){
5349 if (sign[i]==0) continue;
5350 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5351 if (!track0) continue;
5352 if (AliTPCReconstructor::StreamLevel()>1){
5357 "zvertex="<<zvertex<<
5358 "sdcar0="<<sdcar[i]<<
5359 "cdcar0="<<cdcar[i]<<
5360 "pulldcar0="<<pulldcar[i]<<
5361 "pulldcaz0="<<pulldcaz[i]<<
5362 "pulldca0="<<pulldca[i]<<
5363 "isPrim="<<isPrim[i]<<
5367 if (track0->GetSigned1Pt()<0) continue;
5368 if (track0->GetKinkIndex(0)>0||isPrim[i]) continue; //daughter kink
5370 if (TMath::Abs(helixes[i].GetHelix(4))<0.000000001) continue;
5375 for (Int_t j =0;j<nentries;j++){
5376 AliTPCseed * track1 = (AliTPCseed*)array->At(j);
5377 if (!track1) continue;
5378 if (track1->GetKinkIndex(0)>0 || isPrim[j]) continue; //daughter kink
5379 if (sign[j]*sign[i]>0) continue;
5380 if (TMath::Abs(helixes[j].GetHelix(4))<0.000001) continue;
5381 if (track0->GetCircular()+track1->GetCircular()>1) continue; //circling -returning track
5384 // DCA to prim vertex cut
5390 Int_t npoints = helixes[i].GetRPHIintersections(helixes[j], phase, radius,cdist2);
5391 if (npoints<1) continue;
5395 if (radius[0]<cradius0||radius[0]>cradius1) continue;
5396 helixes[i].LinearDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta[0]);
5397 if (delta[0]>cdist1) continue;
5400 if (TMath::Max(radius[0],radius[1])<cradius0|| TMath::Min(radius[0],radius[1])>cradius1) continue;
5401 helixes[i].LinearDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta[0]);
5402 helixes[i].LinearDCA(helixes[j],phase[1][0],phase[1][1],radius[1],delta[1]);
5403 if (delta[1]<delta[0]) iclosest=1;
5404 if (delta[iclosest]>cdist1) continue;
5406 helixes[i].ParabolicDCA(helixes[j],phase[iclosest][0],phase[iclosest][1],radius[iclosest],delta[iclosest]);
5407 if (radius[iclosest]<cradius0 || radius[iclosest]>cradius1 || delta[iclosest]>cdist1) continue;
5409 Double_t pointAngle = helixes[i].GetPointAngle(helixes[j],phase[iclosest],fprimvertex);
5410 if (pointAngle<cpointAngle) continue;
5412 Bool_t isGamma = kFALSE;
5413 vertex.SetParamP(*track0); //track0 - plus
5414 vertex.SetParamN(*track1); //track1 - minus
5415 vertex.Update(fprimvertex);
5416 if (track0->TPCrPID(0)>0.3&&track1->TPCrPID(0)>0.3&&vertex.GetAnglep()[2]<0.15) isGamma=kTRUE; // gamma conversion candidate
5417 Double_t pointAngle2 = vertex.GetV0CosineOfPointingAngle();
5419 if (vertex.GetV0CosineOfPointingAngle()<cpointAngle && (!isGamma)) continue;// point angle cut
5420 //Bo: if (vertex.GetDist2()>2&&(!isGamma)) continue; // point angle cut
5421 if (vertex.GetDcaV0Daughters()>2&&(!isGamma)) continue;//Bo: // point angle cut
5422 Float_t sigmae = 0.15*0.15;
5423 if (vertex.GetRr()<80)
5424 sigmae += (sdcar[i]*sdcar[i]+sdcar[j]*sdcar[j])*(1.-vertex.GetRr()/80.)*(1.-vertex.GetRr()/80.);
5425 sigmae+= TMath::Sqrt(sigmae);
5426 //Bo: if (vertex.GetDist2()/sigmae>3.&&(!isGamma)) continue;
5427 if (vertex.GetDcaV0Daughters()/sigmae>3.&&(!isGamma)) continue;
5428 Float_t densb0=0,densb1=0,densa0=0,densa1=0;
5429 Int_t row0 = GetRowNumber(vertex.GetRr());
5431 //Bo: if (vertex.GetDist2()>0.2) continue;
5432 if (vertex.GetDcaV0Daughters()>0.2) continue;
5433 densb0 = track0->Density2(0,row0-5);
5434 densb1 = track1->Density2(0,row0-5);
5435 if (densb0>0.3|| densb1>0.3) continue; //clusters before vertex
5436 densa0 = track0->Density2(row0+5,row0+40);
5437 densa1 = track1->Density2(row0+5,row0+40);
5438 if ((densa0<0.4|| densa1<0.4)&&(!isGamma)) continue; //missing clusters after vertex
5441 densa0 = track0->Density2(0,40); //cluster density
5442 densa1 = track1->Density2(0,40); //cluster density
5443 if ((vertex.GetRr()<80&&densa0+densa1<1.)&&(!isGamma)) continue;
5445 //Bo: vertex.SetLab(0,track0->GetLabel());
5446 //Bo: vertex.SetLab(1,track1->GetLabel());
5447 vertex.SetChi2After((densa0+densa1)*0.5);
5448 vertex.SetChi2Before((densb0+densb1)*0.5);
5449 vertex.SetIndex(0,i);
5450 vertex.SetIndex(1,j);
5451 //Bo: vertex.SetStatus(1); // TPC v0 candidate
5452 vertex.SetOnFlyStatus(2);//Bo: // TPC v0 candidate
5453 //Bo: vertex.SetRp(track0->TPCrPIDs());
5454 //Bo: vertex.SetRm(track1->TPCrPIDs());
5455 tpcv0s->AddLast(new AliESDv0(vertex));
5458 Int_t eventNr = esd->GetEventNumberInFile(); // This is most likely NOT the event number you'd like to use. It has nothing to do with the 'real' event number
5459 Double_t radiusm= (delta[0]<delta[1])? TMath::Sqrt(radius[0]):TMath::Sqrt(radius[1]);
5460 Double_t deltam= (delta[0]<delta[1])? TMath::Sqrt(delta[0]):TMath::Sqrt(delta[1]);
5461 if (AliTPCReconstructor::StreamLevel()>1) {
5462 Int_t lab0=track0->GetLabel();
5463 Int_t lab1=track1->GetLabel();
5464 Char_t c0=track0->GetCircular();
5465 Char_t c1=track1->GetCircular();
5468 "vertex.="<<&vertex<<
5471 "Helix0.="<<&helixes[i]<<
5474 "Helix1.="<<&helixes[j]<<
5475 "pointAngle="<<pointAngle<<
5476 "pointAngle2="<<pointAngle2<<
5481 "zvertex="<<zvertex<<
5484 "npoints="<<npoints<<
5485 "radius0="<<radius[0]<<
5486 "delta0="<<delta[0]<<
5487 "radius1="<<radius[1]<<
5488 "delta1="<<delta[1]<<
5489 "radiusm="<<radiusm<<
5491 "sdcar0="<<sdcar[i]<<
5492 "sdcar1="<<sdcar[j]<<
5493 "cdcar0="<<cdcar[i]<<
5494 "cdcar1="<<cdcar[j]<<
5495 "pulldcar0="<<pulldcar[i]<<
5496 "pulldcar1="<<pulldcar[j]<<
5497 "pulldcaz0="<<pulldcaz[i]<<
5498 "pulldcaz1="<<pulldcaz[j]<<
5499 "pulldca0="<<pulldca[i]<<
5500 "pulldca1="<<pulldca[j]<<
5510 Float_t *quality = new Float_t[ncandidates];
5511 Int_t *indexes = new Int_t[ncandidates];
5513 for (Int_t i=0;i<ncandidates;i++){
5515 AliESDv0 *v0 = (AliESDv0*)tpcv0s->At(i);
5516 quality[i] = 1./(1.00001-v0->GetV0CosineOfPointingAngle()); //base point angle
5517 // quality[i] /= (0.5+v0->GetDist2());
5518 // quality[i] *= v0->GetChi2After(); //density factor
5520 Int_t index0 = v0->GetIndex(0);
5521 Int_t index1 = v0->GetIndex(1);
5522 //Bo: Double_t minpulldca = TMath::Min(2.+pulldca[v0->GetIndex(0)],(2.+pulldca[v0->GetIndex(1)]) ); //pull
5523 Double_t minpulldca = TMath::Min(2.+pulldca[index0],(2.+pulldca[index1]) );//Bo:
5527 AliTPCseed * track0 = (AliTPCseed*)array->At(index0);
5528 AliTPCseed * track1 = (AliTPCseed*)array->At(index1);
5529 if (track0->TPCrPID(0)>0.3&&track1->TPCrPID(0)>0.3&&v0->GetAnglep()[2]<0.15) quality[i]+=1000000; // gamma conversion candidate
5530 if (track0->TPCrPID(4)>0.9||track1->TPCrPID(4)>0.9&&minpulldca>4) quality[i]*=10; // lambda candidate candidate
5533 TMath::Sort(ncandidates,quality,indexes,kTRUE);
5536 for (Int_t i=0;i<ncandidates;i++){
5537 AliESDv0 * v0 = (AliESDv0*)tpcv0s->At(indexes[i]);
5539 Int_t index0 = v0->GetIndex(0);
5540 Int_t index1 = v0->GetIndex(1);
5541 AliTPCseed * track0 = (AliTPCseed*)array->At(index0);
5542 AliTPCseed * track1 = (AliTPCseed*)array->At(index1);
5543 if (!track0||!track1) {
5547 Bool_t accept =kTRUE; //default accept
5548 Int_t *v0indexes0 = track0->GetV0Indexes();
5549 Int_t *v0indexes1 = track1->GetV0Indexes();
5551 Int_t order0 = (v0indexes0[0]!=0) ? 1:0;
5552 Int_t order1 = (v0indexes1[0]!=0) ? 1:0;
5553 if (v0indexes0[1]!=0) order0 =2;
5554 if (v0indexes1[1]!=0) order1 =2;
5556 if (v0indexes0[2]!=0) {order0=3; accept=kFALSE;}
5557 if (v0indexes0[2]!=0) {order1=3; accept=kFALSE;}
5559 AliESDv0 * v02 = v0;
5561 //Bo: v0->SetOrder(0,order0);
5562 //Bo: v0->SetOrder(1,order1);
5563 //Bo: v0->SetOrder(1,order0+order1);
5564 v0->SetOnFlyStatus(kTRUE);
5565 Int_t index = esd->AddV0(v0);
5566 v02 = esd->GetV0(index);
5567 v0indexes0[order0]=index;
5568 v0indexes1[order1]=index;
5572 Int_t eventNr = esd->GetEventNumberInFile(); // This is most likely NOT the event number you'd like to use. It has nothing to do with the 'real' event number
5573 if (AliTPCReconstructor::StreamLevel()>1) {
5574 Int_t lab0=track0->GetLabel();
5575 Int_t lab1=track1->GetLabel();
5584 "dca0="<<dca[index0]<<
5585 "dca1="<<dca[index1]<<
5589 "quality="<<quality[i]<<
5590 "pulldca0="<<pulldca[index0]<<
5591 "pulldca1="<<pulldca[index1]<<
5615 printf("TPC V0 finder : naccepted\t%d\tncandidates\t%d\tntracks\t%d\tnall\t%d\n",naccepted,ncandidates,ntracks,nall);
5619 Int_t AliTPCtrackerMI::RefitKink(AliTPCseed &mother, AliTPCseed &daughter, AliESDkink &knk)
5622 // refit kink towards to the vertex
5625 AliKink &kink=(AliKink &)knk;
5627 Int_t row0 = GetRowNumber(kink.GetR());
5628 FollowProlongation(mother,0);
5629 mother.Reset(kFALSE);
5631 FollowProlongation(daughter,row0);
5632 daughter.Reset(kFALSE);
5633 FollowBackProlongation(daughter,158);
5634 daughter.Reset(kFALSE);
5635 Int_t first = TMath::Max(row0-20,30);
5636 Int_t last = TMath::Min(row0+20,140);
5638 const Int_t kNdiv =5;
5639 AliTPCseed param0[kNdiv]; // parameters along the track
5640 AliTPCseed param1[kNdiv]; // parameters along the track
5641 AliKink kinks[kNdiv]; // corresponding kink parameters
5644 for (Int_t irow=0; irow<kNdiv;irow++){
5645 rows[irow] = first +((last-first)*irow)/(kNdiv-1);
5647 // store parameters along the track
5649 for (Int_t irow=0;irow<kNdiv;irow++){
5650 FollowBackProlongation(mother, rows[irow]);
5651 FollowProlongation(daughter,rows[kNdiv-1-irow]);
5652 param0[irow] = mother;
5653 param1[kNdiv-1-irow] = daughter;
5657 for (Int_t irow=0; irow<kNdiv-1;irow++){
5658 if (param0[irow].GetNumberOfClusters()<kNdiv||param1[irow].GetNumberOfClusters()<kNdiv) continue;
5659 kinks[irow].SetMother(param0[irow]);
5660 kinks[irow].SetDaughter(param1[irow]);
5661 kinks[irow].Update();
5664 // choose kink with best "quality"
5666 Double_t mindist = 10000;
5667 for (Int_t irow=0;irow<kNdiv;irow++){
5668 if (param0[irow].GetNumberOfClusters()<20||param1[irow].GetNumberOfClusters()<20) continue;
5669 if (TMath::Abs(kinks[irow].GetR())>240.) continue;
5670 if (TMath::Abs(kinks[irow].GetR())<100.) continue;
5672 Float_t normdist = TMath::Abs(param0[irow].GetX()-kinks[irow].GetR())*(0.1+kink.GetDistance());
5673 normdist/= (param0[irow].GetNumberOfClusters()+param1[irow].GetNumberOfClusters()+40.);
5674 if (normdist < mindist){
5680 if (index==-1) return 0;
5683 param0[index].Reset(kTRUE);
5684 FollowProlongation(param0[index],0);
5686 mother = param0[index];
5687 daughter = param1[index]; // daughter in vertex
5689 kink.SetMother(mother);
5690 kink.SetDaughter(daughter);
5692 kink.SetTPCRow0(GetRowNumber(kink.GetR()));
5693 kink.SetTPCncls(param0[index].GetNumberOfClusters(),0);
5694 kink.SetTPCncls(param1[index].GetNumberOfClusters(),1);
5695 kink.SetLabel(CookLabel(&mother,0.4, 0,kink.GetTPCRow0()),0);
5696 kink.SetLabel(CookLabel(&daughter,0.4, kink.GetTPCRow0(),160),1);
5697 mother.SetLabel(kink.GetLabel(0));
5698 daughter.SetLabel(kink.GetLabel(1));
5704 void AliTPCtrackerMI::UpdateKinkQualityM(AliTPCseed * seed){
5706 // update Kink quality information for mother after back propagation
5708 if (seed->GetKinkIndex(0)>=0) return;
5709 for (Int_t ikink=0;ikink<3;ikink++){
5710 Int_t index = seed->GetKinkIndex(ikink);
5711 if (index>=0) break;
5712 index = TMath::Abs(index)-1;
5713 AliESDkink * kink = fEvent->GetKink(index);
5714 kink->SetTPCDensity(-1,0,0);
5715 kink->SetTPCDensity(1,0,1);
5717 Int_t row0 = kink->GetTPCRow0() - 2 - Int_t( 0.5/ (0.05+kink->GetAngle(2)));
5718 if (row0<15) row0=15;
5720 Int_t row1 = kink->GetTPCRow0() + 2 + Int_t( 0.5/ (0.05+kink->GetAngle(2)));
5721 if (row1>145) row1=145;
5723 Int_t found,foundable,shared;
5724 seed->GetClusterStatistic(0,row0, found, foundable,shared,kFALSE);
5725 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),0,0);
5726 seed->GetClusterStatistic(row1,155, found, foundable,shared,kFALSE);
5727 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),0,1);
5732 void AliTPCtrackerMI::UpdateKinkQualityD(AliTPCseed * seed){
5734 // update Kink quality information for daughter after refit
5736 if (seed->GetKinkIndex(0)<=0) return;
5737 for (Int_t ikink=0;ikink<3;ikink++){
5738 Int_t index = seed->GetKinkIndex(ikink);
5739 if (index<=0) break;
5740 index = TMath::Abs(index)-1;
5741 AliESDkink * kink = fEvent->GetKink(index);
5742 kink->SetTPCDensity(-1,1,0);
5743 kink->SetTPCDensity(-1,1,1);
5745 Int_t row0 = kink->GetTPCRow0() -2 - Int_t( 0.5/ (0.05+kink->GetAngle(2)));
5746 if (row0<15) row0=15;
5748 Int_t row1 = kink->GetTPCRow0() +2 + Int_t( 0.5/ (0.05+kink->GetAngle(2)));
5749 if (row1>145) row1=145;
5751 Int_t found,foundable,shared;
5752 seed->GetClusterStatistic(0,row0, found, foundable,shared,kFALSE);
5753 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),1,0);
5754 seed->GetClusterStatistic(row1,155, found, foundable,shared,kFALSE);
5755 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),1,1);
5761 Int_t AliTPCtrackerMI::CheckKinkPoint(AliTPCseed*seed,AliTPCseed &mother, AliTPCseed &daughter, AliESDkink &knk)
5764 // check kink point for given track
5765 // if return value=0 kink point not found
5766 // otherwise seed0 correspond to mother particle
5767 // seed1 correspond to daughter particle
5768 // kink parameter of kink point
5769 AliKink &kink=(AliKink &)knk;
5771 Int_t middlerow = (seed->GetFirstPoint()+seed->GetLastPoint())/2;
5772 Int_t first = seed->GetFirstPoint();
5773 Int_t last = seed->GetLastPoint();
5774 if (last-first<20) return 0; // shortest length - 2*30 = 60 pad-rows
5777 AliTPCseed *seed1 = ReSeed(seed,middlerow+20, kTRUE); //middle of chamber
5778 if (!seed1) return 0;
5779 FollowProlongation(*seed1,seed->GetLastPoint()-20);
5780 seed1->Reset(kTRUE);
5781 FollowProlongation(*seed1,158);
5782 seed1->Reset(kTRUE);
5783 last = seed1->GetLastPoint();
5785 AliTPCseed *seed0 = new AliTPCseed(*seed);
5786 seed0->Reset(kFALSE);
5789 AliTPCseed param0[20]; // parameters along the track
5790 AliTPCseed param1[20]; // parameters along the track
5791 AliKink kinks[20]; // corresponding kink parameters
5793 for (Int_t irow=0; irow<20;irow++){
5794 rows[irow] = first +((last-first)*irow)/19;
5796 // store parameters along the track
5798 for (Int_t irow=0;irow<20;irow++){
5799 FollowBackProlongation(*seed0, rows[irow]);
5800 FollowProlongation(*seed1,rows[19-irow]);
5801 param0[irow] = *seed0;
5802 param1[19-irow] = *seed1;
5806 for (Int_t irow=0; irow<19;irow++){
5807 kinks[irow].SetMother(param0[irow]);
5808 kinks[irow].SetDaughter(param1[irow]);
5809 kinks[irow].Update();
5812 // choose kink with biggest change of angle
5814 Double_t maxchange= 0;
5815 for (Int_t irow=1;irow<19;irow++){
5816 if (TMath::Abs(kinks[irow].GetR())>240.) continue;
5817 if (TMath::Abs(kinks[irow].GetR())<110.) continue;
5818 Float_t quality = TMath::Abs(kinks[irow].GetAngle(2))/(3.+TMath::Abs(kinks[irow].GetR()-param0[irow].GetX()));
5819 if ( quality > maxchange){
5820 maxchange = quality;
5827 if (index<0) return 0;
5829 Int_t row0 = GetRowNumber(kinks[index].GetR()); //row 0 estimate
5830 seed0 = new AliTPCseed(param0[index]);
5831 seed1 = new AliTPCseed(param1[index]);
5832 seed0->Reset(kFALSE);
5833 seed1->Reset(kFALSE);
5834 seed0->ResetCovariance(10.);
5835 seed1->ResetCovariance(10.);
5836 FollowProlongation(*seed0,0);
5837 FollowBackProlongation(*seed1,158);
5838 mother = *seed0; // backup mother at position 0
5839 seed0->Reset(kFALSE);
5840 seed1->Reset(kFALSE);
5841 seed0->ResetCovariance(10.);
5842 seed1->ResetCovariance(10.);
5844 first = TMath::Max(row0-20,0);
5845 last = TMath::Min(row0+20,158);
5847 for (Int_t irow=0; irow<20;irow++){
5848 rows[irow] = first +((last-first)*irow)/19;
5850 // store parameters along the track
5852 for (Int_t irow=0;irow<20;irow++){
5853 FollowBackProlongation(*seed0, rows[irow]);
5854 FollowProlongation(*seed1,rows[19-irow]);
5855 param0[irow] = *seed0;
5856 param1[19-irow] = *seed1;
5860 for (Int_t irow=0; irow<19;irow++){
5861 kinks[irow].SetMother(param0[irow]);
5862 kinks[irow].SetDaughter(param1[irow]);
5863 // param0[irow].Dump();
5864 //param1[irow].Dump();
5865 kinks[irow].Update();
5868 // choose kink with biggest change of angle
5871 for (Int_t irow=0;irow<20;irow++){
5872 if (TMath::Abs(kinks[irow].GetR())>250.) continue;
5873 if (TMath::Abs(kinks[irow].GetR())<90.) continue;
5874 Float_t quality = TMath::Abs(kinks[irow].GetAngle(2))/(3.+TMath::Abs(kinks[irow].GetR()-param0[irow].GetX()));
5875 if ( quality > maxchange){
5876 maxchange = quality;
5883 if (index==-1 || param0[index].GetNumberOfClusters()+param1[index].GetNumberOfClusters()<100){
5888 // Float_t anglesigma = TMath::Sqrt(param0[index].fC22+param0[index].fC33+param1[index].fC22+param1[index].fC33);
5890 kink.SetMother(param0[index]);
5891 kink.SetDaughter(param1[index]);
5893 row0 = GetRowNumber(kink.GetR());
5894 kink.SetTPCRow0(row0);
5895 kink.SetLabel(CookLabel(seed0,0.5,0,row0),0);
5896 kink.SetLabel(CookLabel(seed1,0.5,row0,158),1);
5897 kink.SetIndex(-10,0);
5898 kink.SetIndex(int(param0[index].GetNumberOfClusters()+param1[index].GetNumberOfClusters()),1);
5899 kink.SetTPCncls(param0[index].GetNumberOfClusters(),0);
5900 kink.SetTPCncls(param1[index].GetNumberOfClusters(),1);
5903 // new (&mother) AliTPCseed(param0[index]);
5904 daughter = param1[index];
5905 daughter.SetLabel(kink.GetLabel(1));
5906 param0[index].Reset(kTRUE);
5907 FollowProlongation(param0[index],0);
5908 mother = param0[index];
5909 mother.SetLabel(kink.GetLabel(0));
5919 AliTPCseed* AliTPCtrackerMI::ReSeed(AliTPCseed *t)
5922 // reseed - refit - track
5925 // Int_t last = fSectors->GetNRows()-1;
5927 if (fSectors == fOuterSec){
5928 first = TMath::Max(first, t->GetFirstPoint()-fInnerSec->GetNRows());
5932 first = t->GetFirstPoint();
5934 AliTPCseed * seed = MakeSeed(t,0.1,0.5,0.9);
5935 FollowBackProlongation(*t,fSectors->GetNRows()-1);
5937 FollowProlongation(*t,first);
5947 //_____________________________________________________________________________
5948 Int_t AliTPCtrackerMI::ReadSeeds(const TFile *inp) {
5949 //-----------------------------------------------------------------
5950 // This function reades track seeds.
5951 //-----------------------------------------------------------------
5952 TDirectory *savedir=gDirectory;
5954 TFile *in=(TFile*)inp;
5955 if (!in->IsOpen()) {
5956 cerr<<"AliTPCtrackerMI::ReadSeeds(): input file is not open !\n";
5961 TTree *seedTree=(TTree*)in->Get("Seeds");
5963 cerr<<"AliTPCtrackerMI::ReadSeeds(): ";
5964 cerr<<"can't get a tree with track seeds !\n";
5967 AliTPCtrack *seed=new AliTPCtrack;
5968 seedTree->SetBranchAddress("tracks",&seed);
5970 if (fSeeds==0) fSeeds=new TObjArray(15000);
5972 Int_t n=(Int_t)seedTree->GetEntries();
5973 for (Int_t i=0; i<n; i++) {
5974 seedTree->GetEvent(i);
5975 fSeeds->AddLast(new AliTPCseed(*seed/*,seed->GetAlpha()*/));
5984 Int_t AliTPCtrackerMI::Clusters2Tracks (AliESDEvent *esd)
5987 if (fSeeds) DeleteSeeds();
5990 if (!fSeeds) return 1;
5997 //_____________________________________________________________________________
5998 Int_t AliTPCtrackerMI::Clusters2Tracks() {
5999 //-----------------------------------------------------------------
6000 // This is a track finder.
6001 //-----------------------------------------------------------------
6002 TDirectory *savedir=gDirectory;
6006 fSeeds = Tracking();
6009 Info("Clusters2Tracks","Time for tracking: \t");timer.Print();timer.Start();
6011 //activate again some tracks
6012 for (Int_t i=0; i<fSeeds->GetEntriesFast(); i++) {
6013 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
6015 Int_t nc=t.GetNumberOfClusters();
6017 delete fSeeds->RemoveAt(i);
6021 if (pt->GetRemoval()==10) {
6022 if (pt->GetDensityFirst(20)>0.8 || pt->GetDensityFirst(30)>0.8 || pt->GetDensityFirst(40)>0.7)
6023 pt->Desactivate(10); // make track again active
6025 pt->Desactivate(20);
6026 delete fSeeds->RemoveAt(i);
6031 RemoveUsed2(fSeeds,0.85,0.85,0);
6032 if (AliTPCReconstructor::GetRecoParam()->GetDoKinks()) FindKinks(fSeeds,fEvent);
6033 //FindCurling(fSeeds, fEvent,0);
6034 if (AliTPCReconstructor::StreamLevel()>2) FindMultiMC(fSeeds, fEvent,0); // find multi found tracks
6035 RemoveUsed2(fSeeds,0.5,0.4,20);
6036 FindSplitted(fSeeds, fEvent,0); // find multi found tracks
6039 // // refit short tracks
6041 Int_t nseed=fSeeds->GetEntriesFast();
6044 for (Int_t i=0; i<nseed; i++) {
6045 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
6047 Int_t nc=t.GetNumberOfClusters();
6049 delete fSeeds->RemoveAt(i);
6052 CookLabel(pt,0.1); //For comparison only
6053 //if ((pt->IsActive() || (pt->fRemoval==10) )&& nc>50 &&pt->GetNumberOfClusters()>0.4*pt->fNFoundable){
6054 if ((pt->IsActive() || (pt->GetRemoval()==10) )){
6056 if (fDebug>0) cerr<<found<<'\r';
6060 delete fSeeds->RemoveAt(i);
6064 //RemoveOverlap(fSeeds,0.99,7,kTRUE);
6066 //RemoveUsed(fSeeds,0.9,0.9,6);
6068 nseed=fSeeds->GetEntriesFast();
6070 for (Int_t i=0; i<nseed; i++) {
6071 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
6073 Int_t nc=t.GetNumberOfClusters();
6075 delete fSeeds->RemoveAt(i);
6079 t.CookdEdx(0.02,0.6);
6080 // CheckKinkPoint(&t,0.05);
6081 //if ((pt->IsActive() || (pt->fRemoval==10) )&& nc>50 &&pt->GetNumberOfClusters()>0.4*pt->fNFoundable){
6082 if ((pt->IsActive() || (pt->GetRemoval()==10) )){
6090 delete fSeeds->RemoveAt(i);
6091 //AliTPCseed * seed1 = ReSeed(pt,0.05,0.5,1);
6093 // FollowProlongation(*seed1,0);
6094 // Int_t n = seed1->GetNumberOfClusters();
6095 // printf("fP4\t%f\t%f\n",seed1->GetC(),pt->GetC());
6096 // printf("fN\t%d\t%d\n", seed1->GetNumberOfClusters(),pt->GetNumberOfClusters());
6099 //AliTPCseed * seed2 = ReSeed(pt,0.95,0.5,0.05);
6103 SortTracks(fSeeds, 1);
6107 PrepareForBackProlongation(fSeeds,5.);
6108 PropagateBack(fSeeds);
6109 printf("Time for back propagation: \t");timer.Print();timer.Start();
6113 PrepareForProlongation(fSeeds,5.);
6114 PropagateForward2(fSeeds);
6116 printf("Time for FORWARD propagation: \t");timer.Print();timer.Start();
6117 // RemoveUsed(fSeeds,0.7,0.7,6);
6118 //RemoveOverlap(fSeeds,0.9,7,kTRUE);
6120 nseed=fSeeds->GetEntriesFast();
6122 for (Int_t i=0; i<nseed; i++) {
6123 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
6125 Int_t nc=t.GetNumberOfClusters();
6127 delete fSeeds->RemoveAt(i);
6130 t.CookdEdx(0.02,0.6);
6131 // CookLabel(pt,0.1); //For comparison only
6132 //if ((pt->IsActive() || (pt->fRemoval==10) )&& nc>50 &&pt->GetNumberOfClusters()>0.4*pt->fNFoundable){
6133 if ((pt->IsActive() || (pt->fRemoval==10) )){
6134 cerr<<found++<<'\r';
6137 delete fSeeds->RemoveAt(i);
6142 // fNTracks = found;
6144 Info("Clusters2Tracks","Time for overlap removal, track writing and dedx cooking: \t"); timer.Print();timer.Start();
6147 // cerr<<"Number of found tracks : "<<"\t"<<found<<endl;
6148 Info("Clusters2Tracks","Number of found tracks %d",found);
6150 // UnloadClusters();
6155 void AliTPCtrackerMI::Tracking(TObjArray * arr)
6158 // tracking of the seeds
6161 fSectors = fOuterSec;
6162 ParallelTracking(arr,150,63);
6163 fSectors = fOuterSec;
6164 ParallelTracking(arr,63,0);
6167 TObjArray * AliTPCtrackerMI::Tracking(Int_t seedtype, Int_t i1, Int_t i2, Float_t cuts[4], Float_t dy, Int_t dsec)
6172 TObjArray * arr = new TObjArray;
6174 fSectors = fOuterSec;
6177 for (Int_t sec=0;sec<fkNOS;sec++){
6178 if (seedtype==3) MakeSeeds3(arr,sec,i1,i2,cuts,dy, dsec);
6179 if (seedtype==4) MakeSeeds5(arr,sec,i1,i2,cuts,dy);
6180 if (seedtype==2) MakeSeeds2(arr,sec,i1,i2,cuts,dy);
6183 Info("Tracking","\nSeeding - %d\t%d\t%d\t%d\n",seedtype,i1,i2,arr->GetEntriesFast());
6195 TObjArray * AliTPCtrackerMI::Tracking()
6199 if (AliTPCReconstructor::GetRecoParam()->GetSpecialSeeding()) return TrackingSpecial();
6202 Int_t nup=fOuterSec->GetNRows()+fInnerSec->GetNRows();
6204 TObjArray * seeds = new TObjArray;
6213 Float_t fnumber = 3.0;
6214 Float_t fdensity = 3.0;
6219 for (Int_t delta = 0; delta<18; delta+=6){
6223 arr = Tracking(3,nup-1-delta,nup-1-delta-gap,cuts,-1,1);
6224 SumTracks(seeds,arr);
6225 SignClusters(seeds,fnumber,fdensity);
6227 for (Int_t i=2;i<6;i+=2){
6228 // seed high pt tracks
6231 arr = Tracking(3,nup-i-delta,nup-i-delta-gap,cuts,-1,0);
6232 SumTracks(seeds,arr);
6233 SignClusters(seeds,fnumber,fdensity);
6238 // RemoveUsed(seeds,0.9,0.9,1);
6239 // UnsignClusters();
6240 // SignClusters(seeds,fnumber,fdensity);
6244 for (Int_t delta = 20; delta<120; delta+=10){
6246 // seed high pt tracks
6250 arr = Tracking(3,nup-delta,nup-delta-gap,cuts,-1);
6251 SumTracks(seeds,arr);
6252 SignClusters(seeds,fnumber,fdensity);
6257 arr = Tracking(3,nup-delta-5,nup-delta-5-gap,cuts,-1);
6258 SumTracks(seeds,arr);
6259 SignClusters(seeds,fnumber,fdensity);
6270 Info("Tracking()","\n\nPrimary seeding\t%d\n\n",seeds->GetEntriesFast());
6274 // RemoveUsed(seeds,0.75,0.75,1);
6276 //SignClusters(seeds,fnumber,fdensity);
6285 arr = Tracking(4,nup-1,nup-1-gap,cuts,-1);
6286 SumTracks(seeds,arr);
6287 SignClusters(seeds,fnumber,fdensity);
6289 arr = Tracking(4,nup-2,nup-2-gap,cuts,-1);
6290 SumTracks(seeds,arr);
6291 SignClusters(seeds,fnumber,fdensity);
6293 arr = Tracking(4,nup-3,nup-3-gap,cuts,-1);
6294 SumTracks(seeds,arr);
6295 SignClusters(seeds,fnumber,fdensity);
6299 for (Int_t delta = 3; delta<30; delta+=5){
6305 arr = Tracking(4,nup-1-delta,nup-1-delta-gap,cuts,-1);
6306 SumTracks(seeds,arr);
6307 SignClusters(seeds,fnumber,fdensity);
6309 arr = Tracking(4,nup-3-delta,nup-5-delta-gap,cuts,4);
6310 SumTracks(seeds,arr);
6311 SignClusters(seeds,fnumber,fdensity);
6322 Int_t fLastSeedRowSec=AliTPCReconstructor::GetRecoParam()->GetLastSeedRowSec();
6325 for (Int_t delta = 30; delta<fLastSeedRowSec; delta+=10){
6331 arr = Tracking(4,nup-1-delta,nup-1-delta-gap,cuts,-1);
6332 SumTracks(seeds,arr);
6333 SignClusters(seeds,fnumber,fdensity);
6335 arr = Tracking(4,nup-5-delta,nup-5-delta-gap,cuts,5 );
6336 SumTracks(seeds,arr);
6337 SignClusters(seeds,fnumber,fdensity);
6341 Info("Tracking()","\n\nSecondary seeding\t%d\n\n",seeds->GetEntriesFast());
6352 TObjArray * AliTPCtrackerMI::TrackingSpecial()
6355 // seeding adjusted for laser and cosmic tests - short tracks with big inclination angle
6356 // no primary vertex seeding tried
6360 Int_t nup=fOuterSec->GetNRows()+fInnerSec->GetNRows();
6362 TObjArray * seeds = new TObjArray;
6367 Float_t fnumber = 3.0;
6368 Float_t fdensity = 3.0;
6371 cuts[0] = AliTPCReconstructor::GetRecoParam()->GetMaxC(); // max curvature
6372 cuts[1] = 3.5; // max tan(phi) angle for seeding
6373 cuts[2] = 3.; // not used (cut on z primary vertex)
6374 cuts[3] = 3.5; // max tan(theta) angle for seeding
6376 for (Int_t delta = 0; nup-delta-gap-1>0; delta+=3){
6378 arr = Tracking(4,nup-1-delta,nup-1-delta-gap,cuts,-1);
6379 SumTracks(seeds,arr);
6380 SignClusters(seeds,fnumber,fdensity);
6384 Info("Tracking()","\n\nSecondary seeding\t%d\n\n",seeds->GetEntriesFast());
6395 void AliTPCtrackerMI::SumTracks(TObjArray *arr1,TObjArray *&arr2) const
6398 //sum tracks to common container
6399 //remove suspicious tracks
6400 Int_t nseed = arr2->GetEntriesFast();
6401 for (Int_t i=0;i<nseed;i++){
6402 AliTPCseed *pt=(AliTPCseed*)arr2->UncheckedAt(i);
6405 // remove tracks with too big curvature
6407 if (TMath::Abs(pt->GetC())>AliTPCReconstructor::GetRecoParam()->GetMaxC()){
6408 delete arr2->RemoveAt(i);
6411 // REMOVE VERY SHORT TRACKS
6412 if (pt->GetNumberOfClusters()<20){
6413 delete arr2->RemoveAt(i);
6416 // NORMAL ACTIVE TRACK
6417 if (pt->IsActive()){
6418 arr1->AddLast(arr2->RemoveAt(i));
6421 //remove not usable tracks
6422 if (pt->GetRemoval()!=10){
6423 delete arr2->RemoveAt(i);
6427 // ENABLE ONLY ENOUGH GOOD STOPPED TRACKS
6428 if (pt->GetDensityFirst(20)>0.8 || pt->GetDensityFirst(30)>0.8 || pt->GetDensityFirst(40)>0.7)
6429 arr1->AddLast(arr2->RemoveAt(i));
6431 delete arr2->RemoveAt(i);
6435 delete arr2; arr2 = 0;
6440 void AliTPCtrackerMI::ParallelTracking(TObjArray * arr, Int_t rfirst, Int_t rlast)
6443 // try to track in parralel
6445 Int_t nseed=arr->GetEntriesFast();
6446 //prepare seeds for tracking
6447 for (Int_t i=0; i<nseed; i++) {
6448 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i), &t=*pt;
6450 if (!t.IsActive()) continue;
6451 // follow prolongation to the first layer
6452 if ( (fSectors ==fInnerSec) || (t.GetFirstPoint()-fParam->GetNRowLow()>rfirst+1) )
6453 FollowProlongation(t, rfirst+1);
6458 for (Int_t nr=rfirst; nr>=rlast; nr--){
6459 if (nr<fInnerSec->GetNRows())
6460 fSectors = fInnerSec;
6462 fSectors = fOuterSec;
6463 // make indexes with the cluster tracks for given
6465 // find nearest cluster
6466 for (Int_t i=0; i<nseed; i++) {
6467 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i), &t=*pt;
6469 if (nr==80) pt->UpdateReference();
6470 if (!pt->IsActive()) continue;
6471 // if ( (fSectors ==fOuterSec) && (pt->fFirstPoint-fParam->GetNRowLow())<nr) continue;
6472 if (pt->GetRelativeSector()>17) {
6475 UpdateClusters(t,nr);
6477 // prolonagate to the nearest cluster - if founded
6478 for (Int_t i=0; i<nseed; i++) {
6479 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
6481 if (!pt->IsActive()) continue;
6482 // if ((fSectors ==fOuterSec) && (pt->fFirstPoint-fParam->GetNRowLow())<nr) continue;
6483 if (pt->GetRelativeSector()>17) {
6486 FollowToNextCluster(*pt,nr);
6491 void AliTPCtrackerMI::PrepareForBackProlongation(TObjArray * arr,Float_t fac) const
6495 // if we use TPC track itself we have to "update" covariance
6497 Int_t nseed= arr->GetEntriesFast();
6498 for (Int_t i=0;i<nseed;i++){
6499 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
6503 //rotate to current local system at first accepted point
6504 Int_t index = pt->GetClusterIndex2(pt->GetFirstPoint());
6505 Int_t sec = (index&0xff000000)>>24;
6507 Float_t angle1 = fInnerSec->GetAlpha()*sec+fInnerSec->GetAlphaShift();
6508 if (angle1>TMath::Pi())
6509 angle1-=2.*TMath::Pi();
6510 Float_t angle2 = pt->GetAlpha();
6512 if (TMath::Abs(angle1-angle2)>0.001){
6513 pt->Rotate(angle1-angle2);
6514 //angle2 = pt->GetAlpha();
6515 //pt->fRelativeSector = pt->GetAlpha()/fInnerSec->GetAlpha();
6516 //if (pt->GetAlpha()<0)
6517 // pt->fRelativeSector+=18;
6518 //sec = pt->fRelativeSector;
6527 void AliTPCtrackerMI::PrepareForProlongation(TObjArray * arr, Float_t fac) const
6531 // if we use TPC track itself we have to "update" covariance
6533 Int_t nseed= arr->GetEntriesFast();
6534 for (Int_t i=0;i<nseed;i++){
6535 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
6538 pt->SetFirstPoint(pt->GetLastPoint());
6546 Int_t AliTPCtrackerMI::PropagateBack(TObjArray * arr)
6549 // make back propagation
6551 Int_t nseed= arr->GetEntriesFast();
6552 for (Int_t i=0;i<nseed;i++){
6553 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
6554 if (pt&& pt->GetKinkIndex(0)<=0) {
6555 //AliTPCseed *pt2 = new AliTPCseed(*pt);
6556 fSectors = fInnerSec;
6557 //FollowBackProlongation(*pt,fInnerSec->GetNRows()-1);
6558 //fSectors = fOuterSec;
6559 FollowBackProlongation(*pt,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1);
6560 //if (pt->GetNumberOfClusters()<(pt->fEsd->GetTPCclusters(0)) ){
6561 // Error("PropagateBack","Not prolonged track %d",pt->GetLabel());
6562 // FollowBackProlongation(*pt2,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1);
6565 if (pt&& pt->GetKinkIndex(0)>0) {
6566 AliESDkink * kink = fEvent->GetKink(pt->GetKinkIndex(0)-1);
6567 pt->SetFirstPoint(kink->GetTPCRow0());
6568 fSectors = fInnerSec;
6569 FollowBackProlongation(*pt,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1);
6577 Int_t AliTPCtrackerMI::PropagateForward2(TObjArray * arr)
6580 // make forward propagation
6582 Int_t nseed= arr->GetEntriesFast();
6584 for (Int_t i=0;i<nseed;i++){
6585 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
6587 FollowProlongation(*pt,0);
6596 Int_t AliTPCtrackerMI::PropagateForward()
6599 // propagate track forward
6601 Int_t nseed = fSeeds->GetEntriesFast();
6602 for (Int_t i=0;i<nseed;i++){
6603 AliTPCseed *pt = (AliTPCseed*)fSeeds->UncheckedAt(i);
6605 AliTPCseed &t = *pt;
6606 Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
6607 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
6608 if (alpha < 0. ) alpha += 2.*TMath::Pi();
6609 t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN);
6613 fSectors = fOuterSec;
6614 ParallelTracking(fSeeds,fOuterSec->GetNRows()+fInnerSec->GetNRows()-1,fInnerSec->GetNRows());
6615 fSectors = fInnerSec;
6616 ParallelTracking(fSeeds,fInnerSec->GetNRows()-1,0);
6625 Int_t AliTPCtrackerMI::PropagateBack(AliTPCseed * pt, Int_t row0, Int_t row1)
6628 // make back propagation, in between row0 and row1
6632 fSectors = fInnerSec;
6635 if (row1<fSectors->GetNRows())
6638 r1 = fSectors->GetNRows()-1;
6640 if (row0<fSectors->GetNRows()&& r1>0 )
6641 FollowBackProlongation(*pt,r1);
6642 if (row1<=fSectors->GetNRows())
6645 r1 = row1 - fSectors->GetNRows();
6646 if (r1<=0) return 0;
6647 if (r1>=fOuterSec->GetNRows()) return 0;
6648 fSectors = fOuterSec;
6649 return FollowBackProlongation(*pt,r1);
6657 void AliTPCtrackerMI::GetShape(AliTPCseed * seed, Int_t row)
6661 AliTPCClusterParam * clparam = AliTPCcalibDB::Instance()->GetClusterParam();
6662 Float_t zdrift = TMath::Abs((fParam->GetZLength(0)-TMath::Abs(seed->GetZ())));
6663 Int_t type = (seed->GetSector() < fParam->GetNSector()/2) ? 0: (row>126) ? 1:2;
6664 Double_t angulary = seed->GetSnp();
6665 angulary = angulary*angulary/(1.-angulary*angulary);
6666 Double_t angularz = seed->GetTgl()*seed->GetTgl()*(1.+angulary);
6668 Double_t sigmay = clparam->GetRMS0(0,type,zdrift,TMath::Sqrt(TMath::Abs(angulary)));
6669 Double_t sigmaz = clparam->GetRMS0(1,type,zdrift,TMath::Sqrt(TMath::Abs(angularz)));
6670 seed->SetCurrentSigmaY2(sigmay*sigmay);
6671 seed->SetCurrentSigmaZ2(sigmaz*sigmaz);
6672 // Float_t sd2 = TMath::Abs((fParam->GetZLength(0)-TMath::Abs(seed->GetZ())))*fParam->GetDiffL()*fParam->GetDiffL();
6673 // // Float_t padlength = fParam->GetPadPitchLength(seed->fSector);
6674 // Float_t padlength = GetPadPitchLength(row);
6676 // Float_t sresy = (seed->GetSector() < fParam->GetNSector()/2) ? 0.2 :0.3;
6677 // seed->SetCurrentSigmaY2(sd2+padlength*padlength*angulary/12.+sresy*sresy);
6679 // Float_t sresz = fParam->GetZSigma();
6680 // seed->SetCurrentSigmaZ2(sd2+padlength*padlength*angularz*angularz*(1+angulary)/12.+sresz*sresz);
6682 Float_t wy = GetSigmaY(seed);
6683 Float_t wz = GetSigmaZ(seed);
6686 if (TMath::Abs(wy/seed->fCurrentSigmaY2-1)>0.0001 || TMath::Abs(wz/seed->fCurrentSigmaZ2-1)>0.0001 ){
6687 printf("problem\n");
6694 //__________________________________________________________________________
6695 void AliTPCtrackerMI::CookLabel(AliKalmanTrack *tk, Float_t wrong) const {
6696 //--------------------------------------------------------------------
6697 //This function "cooks" a track label. If label<0, this track is fake.
6698 //--------------------------------------------------------------------
6699 AliTPCseed * t = dynamic_cast<AliTPCseed*>(tk);
6701 printf("%s:%d wrong type \n",(char*)__FILE__,__LINE__);
6705 Int_t noc=t->GetNumberOfClusters();
6707 //printf("\nnot founded prolongation\n\n\n");
6713 AliTPCclusterMI *clusters[160];
6715 for (Int_t i=0;i<160;i++) {
6722 for (i=0; i<160 && current<noc; i++) {
6724 Int_t index=t->GetClusterIndex2(i);
6725 if (index<=0) continue;
6726 if (index&0x8000) continue;
6728 //clusters[current]=GetClusterMI(index);
6729 if (t->GetClusterPointer(i)){
6730 clusters[current]=t->GetClusterPointer(i);
6736 Int_t lab=123456789;
6737 for (i=0; i<noc; i++) {
6738 AliTPCclusterMI *c=clusters[i];
6740 lab=TMath::Abs(c->GetLabel(0));
6742 for (j=0; j<noc; j++) if (lb[j]==lab || mx[j]==0) break;
6748 for (i=0; i<noc; i++) if (mx[i]>max) {max=mx[i]; lab=lb[i];}
6750 for (i=0; i<noc; i++) {
6751 AliTPCclusterMI *c=clusters[i];
6753 if (TMath::Abs(c->GetLabel(1)) == lab ||
6754 TMath::Abs(c->GetLabel(2)) == lab ) max++;
6757 if ((1.- Float_t(max)/noc) > wrong) lab=-lab;
6760 Int_t tail=Int_t(0.10*noc);
6763 for (i=1; i<=160&&ind<tail; i++) {
6764 // AliTPCclusterMI *c=clusters[noc-i];
6765 AliTPCclusterMI *c=clusters[i];
6767 if (lab == TMath::Abs(c->GetLabel(0)) ||
6768 lab == TMath::Abs(c->GetLabel(1)) ||
6769 lab == TMath::Abs(c->GetLabel(2))) max++;
6772 if (max < Int_t(0.5*tail)) lab=-lab;
6779 //delete[] clusters;
6783 //__________________________________________________________________________
6784 Int_t AliTPCtrackerMI::CookLabel(AliTPCseed *t, Float_t wrong,Int_t first, Int_t last) const {
6785 //--------------------------------------------------------------------
6786 //This function "cooks" a track label. If label<0, this track is fake.
6787 //--------------------------------------------------------------------
6788 Int_t noc=t->GetNumberOfClusters();
6790 //printf("\nnot founded prolongation\n\n\n");
6796 AliTPCclusterMI *clusters[160];
6798 for (Int_t i=0;i<160;i++) {
6805 for (i=0; i<160 && current<noc; i++) {
6806 if (i<first) continue;
6807 if (i>last) continue;
6808 Int_t index=t->GetClusterIndex2(i);
6809 if (index<=0) continue;
6810 if (index&0x8000) continue;
6812 //clusters[current]=GetClusterMI(index);
6813 if (t->GetClusterPointer(i)){
6814 clusters[current]=t->GetClusterPointer(i);
6819 if (noc<5) return -1;
6820 Int_t lab=123456789;
6821 for (i=0; i<noc; i++) {
6822 AliTPCclusterMI *c=clusters[i];
6824 lab=TMath::Abs(c->GetLabel(0));
6826 for (j=0; j<noc; j++) if (lb[j]==lab || mx[j]==0) break;
6832 for (i=0; i<noc; i++) if (mx[i]>max) {max=mx[i]; lab=lb[i];}
6834 for (i=0; i<noc; i++) {
6835 AliTPCclusterMI *c=clusters[i];
6837 if (TMath::Abs(c->GetLabel(1)) == lab ||
6838 TMath::Abs(c->GetLabel(2)) == lab ) max++;
6841 if ((1.- Float_t(max)/noc) > wrong) lab=-lab;
6844 Int_t tail=Int_t(0.10*noc);
6847 for (i=1; i<=160&&ind<tail; i++) {
6848 // AliTPCclusterMI *c=clusters[noc-i];
6849 AliTPCclusterMI *c=clusters[i];
6851 if (lab == TMath::Abs(c->GetLabel(0)) ||
6852 lab == TMath::Abs(c->GetLabel(1)) ||
6853 lab == TMath::Abs(c->GetLabel(2))) max++;
6856 if (max < Int_t(0.5*tail)) lab=-lab;
6859 // t->SetLabel(lab);
6863 //delete[] clusters;
6867 Int_t AliTPCtrackerMI::GetRowNumber(Double_t x[3]) const
6869 //return pad row number for given x vector
6870 Float_t phi = TMath::ATan2(x[1],x[0]);
6871 if(phi<0) phi=2.*TMath::Pi()+phi;
6872 // Get the local angle in the sector philoc
6873 const Float_t kRaddeg = 180/3.14159265358979312;
6874 Float_t phiangle = (Int_t (phi*kRaddeg/20.) + 0.5)*20./kRaddeg;
6875 Double_t localx = x[0]*TMath::Cos(phiangle)-x[1]*TMath::Sin(phiangle);
6876 return GetRowNumber(localx);
6881 void AliTPCtrackerMI::MakeBitmaps(AliTPCseed *t)
6883 //-----------------------------------------------------------------------
6884 // Fill the cluster and sharing bitmaps of the track
6885 //-----------------------------------------------------------------------
6887 Int_t firstpoint = 0;
6888 Int_t lastpoint = 159;
6889 AliTPCTrackerPoint *point;
6891 for (int iter=firstpoint; iter<lastpoint; iter++) {
6892 point = t->GetTrackPoint(iter);
6894 t->SetClusterMapBit(iter, kTRUE);
6895 if (point->IsShared())
6896 t->SetSharedMapBit(iter,kTRUE);
6898 t->SetSharedMapBit(iter, kFALSE);
6901 t->SetClusterMapBit(iter, kFALSE);
6902 t->SetSharedMapBit(iter, kFALSE);