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"
109 #include "AliTPCtrackerMI.h"
110 #include "TStopwatch.h"
111 #include "AliTPCReconstructor.h"
113 #include "TTreeStream.h"
114 #include "AliAlignObj.h"
115 #include "AliTrackPointArray.h"
117 #include "AliTPCcalibDB.h"
118 #include "AliTPCTransform.h"
119 #include "AliTPCClusterParam.h"
123 ClassImp(AliTPCtrackerMI)
126 class AliTPCFastMath {
129 static Double_t FastAsin(Double_t x);
131 static Double_t fgFastAsin[20000]; //lookup table for fast asin computation
134 Double_t AliTPCFastMath::fgFastAsin[20000];
135 AliTPCFastMath gAliTPCFastMath; // needed to fill the LUT
137 AliTPCFastMath::AliTPCFastMath(){
139 // initialized lookup table;
140 for (Int_t i=0;i<10000;i++){
141 fgFastAsin[2*i] = TMath::ASin(i/10000.);
142 fgFastAsin[2*i+1] = (TMath::ASin((i+1)/10000.)-fgFastAsin[2*i]);
146 Double_t AliTPCFastMath::FastAsin(Double_t x){
148 // return asin using lookup table
150 Int_t index = int(x*10000);
151 return fgFastAsin[2*index]+(x*10000.-index)*fgFastAsin[2*index+1];
154 Int_t index = int(x*10000);
155 return -(fgFastAsin[2*index]+(x*10000.-index)*fgFastAsin[2*index+1]);
157 //__________________________________________________________________
158 AliTPCtrackerMI::AliTPCtrackerMI()
180 // default constructor
183 //_____________________________________________________________________
187 Int_t AliTPCtrackerMI::UpdateTrack(AliTPCseed * track, Int_t accept){
189 //update track information using current cluster - track->fCurrentCluster
192 AliTPCclusterMI* c =track->GetCurrentCluster();
193 if (accept>0) track->SetCurrentClusterIndex1(track->GetCurrentClusterIndex1() | 0x8000); //sign not accepted clusters
195 UInt_t i = track->GetCurrentClusterIndex1();
197 Int_t sec=(i&0xff000000)>>24;
198 //Int_t row = (i&0x00ff0000)>>16;
199 track->SetRow((i&0x00ff0000)>>16);
200 track->SetSector(sec);
201 // Int_t index = i&0xFFFF;
202 if (sec>=fParam->GetNInnerSector()) track->SetRow(track->GetRow()+fParam->GetNRowLow());
203 track->SetClusterIndex2(track->GetRow(), i);
204 //track->fFirstPoint = row;
205 //if ( track->fLastPoint<row) track->fLastPoint =row;
206 // if (track->fRow<0 || track->fRow>160) {
207 // printf("problem\n");
209 if (track->GetFirstPoint()>track->GetRow())
210 track->SetFirstPoint(track->GetRow());
211 if (track->GetLastPoint()<track->GetRow())
212 track->SetLastPoint(track->GetRow());
215 track->SetClusterPointer(track->GetRow(),c);
218 Double_t angle2 = track->GetSnp()*track->GetSnp();
220 //SET NEW Track Point
222 if (angle2<1) //PH sometimes angle2 is very big. To be investigated...
224 angle2 = TMath::Sqrt(angle2/(1-angle2));
225 AliTPCTrackerPoint &point =*(track->GetTrackPoint(track->GetRow()));
227 point.SetSigmaY(c->GetSigmaY2()/track->GetCurrentSigmaY2());
228 point.SetSigmaZ(c->GetSigmaZ2()/track->GetCurrentSigmaZ2());
229 point.SetErrY(sqrt(track->GetErrorY2()));
230 point.SetErrZ(sqrt(track->GetErrorZ2()));
232 point.SetX(track->GetX());
233 point.SetY(track->GetY());
234 point.SetZ(track->GetZ());
235 point.SetAngleY(angle2);
236 point.SetAngleZ(track->GetTgl());
237 if (point.IsShared()){
238 track->SetErrorY2(track->GetErrorY2()*4);
239 track->SetErrorZ2(track->GetErrorZ2()*4);
243 Double_t chi2 = track->GetPredictedChi2(track->GetCurrentCluster());
245 // track->SetErrorY2(track->GetErrorY2()*1.3);
246 // track->SetErrorY2(track->GetErrorY2()+0.01);
247 // track->SetErrorZ2(track->GetErrorZ2()*1.3);
248 // track->SetErrorZ2(track->GetErrorZ2()+0.005);
250 if (accept>0) return 0;
251 if (track->GetNumberOfClusters()%20==0){
252 // if (track->fHelixIn){
253 // TClonesArray & larr = *(track->fHelixIn);
254 // Int_t ihelix = larr.GetEntriesFast();
255 // new(larr[ihelix]) AliHelix(*track) ;
258 track->SetNoCluster(0);
259 return track->Update(c,chi2,i);
264 Int_t AliTPCtrackerMI::AcceptCluster(AliTPCseed * seed, AliTPCclusterMI * cluster)
267 // decide according desired precision to accept given
268 // cluster for tracking
269 Double_t sy2=ErrY2(seed,cluster);
270 Double_t sz2=ErrZ2(seed,cluster);
272 Double_t sdistancey2 = sy2+seed->GetSigmaY2();
273 Double_t sdistancez2 = sz2+seed->GetSigmaZ2();
275 Double_t rdistancey2 = (seed->GetCurrentCluster()->GetY()-seed->GetY())*
276 (seed->GetCurrentCluster()->GetY()-seed->GetY())/sdistancey2;
277 Double_t rdistancez2 = (seed->GetCurrentCluster()->GetZ()-seed->GetZ())*
278 (seed->GetCurrentCluster()->GetZ()-seed->GetZ())/sdistancez2;
280 Double_t rdistance2 = rdistancey2+rdistancez2;
283 if (AliTPCReconstructor::StreamLevel()>5) {
284 Float_t rmsy2 = seed->GetCurrentSigmaY2();
285 Float_t rmsz2 = seed->GetCurrentSigmaZ2();
286 Float_t rmsy2p30 = seed->GetCMeanSigmaY2p30();
287 Float_t rmsz2p30 = seed->GetCMeanSigmaZ2p30();
288 Float_t rmsy2p2 = seed->GetCMeanSigmaY2p2();
289 Float_t rmsz2p2 = seed->GetCMeanSigmaZ2p2();
291 AliExternalTrackParam param(*seed);
292 (*fDebugStreamer)<<"ErrParam"<<
299 "rmsy2p30="<<rmsy2p30<<
300 "rmsz2p30="<<rmsz2p30<<
301 "rmsy2p2="<<rmsy2p2<<
302 "rmsz2p2="<<rmsz2p2<<
306 if (rdistance2>16) return 3;
309 if ((rdistancey2>9. || rdistancez2>9.) && cluster->GetType()==0)
310 return 2; //suspisiouce - will be changed
312 if ((rdistancey2>6.25 || rdistancez2>6.25) && cluster->GetType()>0)
313 // strict cut on overlaped cluster
314 return 2; //suspisiouce - will be changed
316 if ( (rdistancey2>1. || rdistancez2>6.25 )
317 && cluster->GetType()<0){
318 seed->SetNFoundable(seed->GetNFoundable()-1);
327 //_____________________________________________________________________________
328 AliTPCtrackerMI::AliTPCtrackerMI(const AliTPCParam *par):
330 fkNIS(par->GetNInnerSector()/2),
332 fkNOS(par->GetNOuterSector()/2),
349 //---------------------------------------------------------------------
350 // The main TPC tracker constructor
351 //---------------------------------------------------------------------
352 fInnerSec=new AliTPCSector[fkNIS];
353 fOuterSec=new AliTPCSector[fkNOS];
356 for (i=0; i<fkNIS; i++) fInnerSec[i].Setup(par,0);
357 for (i=0; i<fkNOS; i++) fOuterSec[i].Setup(par,1);
360 Int_t nrowlow = par->GetNRowLow();
361 Int_t nrowup = par->GetNRowUp();
364 for (Int_t i=0;i<nrowlow;i++){
365 fXRow[i] = par->GetPadRowRadiiLow(i);
366 fPadLength[i]= par->GetPadPitchLength(0,i);
367 fYMax[i] = fXRow[i]*TMath::Tan(0.5*par->GetInnerAngle());
371 for (Int_t i=0;i<nrowup;i++){
372 fXRow[i+nrowlow] = par->GetPadRowRadiiUp(i);
373 fPadLength[i+nrowlow] = par->GetPadPitchLength(60,i);
374 fYMax[i+nrowlow] = fXRow[i+nrowlow]*TMath::Tan(0.5*par->GetOuterAngle());
377 fDebugStreamer = new TTreeSRedirector("TPCdebug.root");
379 //________________________________________________________________________
380 AliTPCtrackerMI::AliTPCtrackerMI(const AliTPCtrackerMI &t):
401 //------------------------------------
402 // dummy copy constructor
403 //------------------------------------------------------------------
406 AliTPCtrackerMI & AliTPCtrackerMI::operator=(const AliTPCtrackerMI& /*r*/){
407 //------------------------------
409 //--------------------------------------------------------------
412 //_____________________________________________________________________________
413 AliTPCtrackerMI::~AliTPCtrackerMI() {
414 //------------------------------------------------------------------
415 // TPC tracker destructor
416 //------------------------------------------------------------------
423 if (fDebugStreamer) delete fDebugStreamer;
427 void AliTPCtrackerMI::FillESD(TObjArray* arr)
431 //fill esds using updated tracks
433 // write tracks to the event
434 // store index of the track
435 Int_t nseed=arr->GetEntriesFast();
436 //FindKinks(arr,fEvent);
437 for (Int_t i=0; i<nseed; i++) {
438 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
441 // pt->PropagateTo(fParam->GetInnerRadiusLow());
442 if (pt->GetKinkIndex(0)<=0){ //don't propagate daughter tracks
443 pt->PropagateTo(fParam->GetInnerRadiusLow());
446 if (( pt->GetPoints()[2]- pt->GetPoints()[0])>5 && pt->GetPoints()[3]>0.8){
448 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
449 iotrack.SetTPCPoints(pt->GetPoints());
450 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
451 iotrack.SetV0Indexes(pt->GetV0Indexes());
452 // iotrack.SetTPCpid(pt->fTPCr);
453 //iotrack.SetTPCindex(i);
454 fEvent->AddTrack(&iotrack);
458 if ( (pt->GetNumberOfClusters()>70)&& (Float_t(pt->GetNumberOfClusters())/Float_t(pt->GetNFoundable()))>0.55) {
460 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
461 iotrack.SetTPCPoints(pt->GetPoints());
462 //iotrack.SetTPCindex(i);
463 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
464 iotrack.SetV0Indexes(pt->GetV0Indexes());
465 // iotrack.SetTPCpid(pt->fTPCr);
466 fEvent->AddTrack(&iotrack);
470 // short tracks - maybe decays
472 if ( (pt->GetNumberOfClusters()>30) && (Float_t(pt->GetNumberOfClusters())/Float_t(pt->GetNFoundable()))>0.70) {
473 Int_t found,foundable,shared;
474 pt->GetClusterStatistic(0,60,found, foundable,shared,kFALSE);
475 if ( (found>20) && (pt->GetNShared()/float(pt->GetNumberOfClusters())<0.2)){
477 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
478 //iotrack.SetTPCindex(i);
479 iotrack.SetTPCPoints(pt->GetPoints());
480 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
481 iotrack.SetV0Indexes(pt->GetV0Indexes());
482 //iotrack.SetTPCpid(pt->fTPCr);
483 fEvent->AddTrack(&iotrack);
488 if ( (pt->GetNumberOfClusters()>20) && (Float_t(pt->GetNumberOfClusters())/Float_t(pt->GetNFoundable()))>0.8) {
489 Int_t found,foundable,shared;
490 pt->GetClusterStatistic(0,60,found, foundable,shared,kFALSE);
491 if (found<20) continue;
492 if (pt->GetNShared()/float(pt->GetNumberOfClusters())>0.2) continue;
495 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
496 iotrack.SetTPCPoints(pt->GetPoints());
497 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
498 iotrack.SetV0Indexes(pt->GetV0Indexes());
499 //iotrack.SetTPCpid(pt->fTPCr);
500 //iotrack.SetTPCindex(i);
501 fEvent->AddTrack(&iotrack);
504 // short tracks - secondaties
506 if ( (pt->GetNumberOfClusters()>30) ) {
507 Int_t found,foundable,shared;
508 pt->GetClusterStatistic(128,158,found, foundable,shared,kFALSE);
509 if ( (found>20) && (pt->GetNShared()/float(pt->GetNumberOfClusters())<0.2) &&float(found)/float(foundable)>0.8){
511 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
512 iotrack.SetTPCPoints(pt->GetPoints());
513 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
514 iotrack.SetV0Indexes(pt->GetV0Indexes());
515 //iotrack.SetTPCpid(pt->fTPCr);
516 //iotrack.SetTPCindex(i);
517 fEvent->AddTrack(&iotrack);
522 if ( (pt->GetNumberOfClusters()>15)) {
523 Int_t found,foundable,shared;
524 pt->GetClusterStatistic(138,158,found, foundable,shared,kFALSE);
525 if (found<15) continue;
526 if (foundable<=0) continue;
527 if (pt->GetNShared()/float(pt->GetNumberOfClusters())>0.2) continue;
528 if (float(found)/float(foundable)<0.8) continue;
531 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
532 iotrack.SetTPCPoints(pt->GetPoints());
533 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
534 iotrack.SetV0Indexes(pt->GetV0Indexes());
535 // iotrack.SetTPCpid(pt->fTPCr);
536 //iotrack.SetTPCindex(i);
537 fEvent->AddTrack(&iotrack);
542 printf("Number of filled ESDs-\t%d\n",fEvent->GetNumberOfTracks());
549 Double_t AliTPCtrackerMI::ErrY2(AliTPCseed* seed, AliTPCclusterMI * cl){
552 // Use calibrated cluster error from OCDB
554 AliTPCClusterParam * clparam = AliTPCcalibDB::Instance()->GetClusterParam();
556 Float_t z = TMath::Abs(fParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
557 Int_t ctype = cl->GetType();
558 Int_t type = (cl->GetRow()<63) ? 0: (cl->GetRow()>126) ? 1:2;
559 Double_t angle = seed->GetSnp()*seed->GetSnp();
560 angle = TMath::Sqrt(TMath::Abs(angle/(1.-angle)));
561 Double_t erry2 = clparam->GetError0Par(0,type, z,angle);
563 erry2+=0.5; // edge cluster
566 seed->SetErrorY2(erry2);
570 //calculate look-up table at the beginning
571 // static Bool_t ginit = kFALSE;
572 // static Float_t gnoise1,gnoise2,gnoise3;
573 // static Float_t ggg1[10000];
574 // static Float_t ggg2[10000];
575 // static Float_t ggg3[10000];
576 // static Float_t glandau1[10000];
577 // static Float_t glandau2[10000];
578 // static Float_t glandau3[10000];
580 // static Float_t gcor01[500];
581 // static Float_t gcor02[500];
582 // static Float_t gcorp[500];
586 // if (ginit==kFALSE){
587 // for (Int_t i=1;i<500;i++){
588 // Float_t rsigma = float(i)/100.;
589 // gcor02[i] = TMath::Max(0.78 +TMath::Exp(7.4*(rsigma-1.2)),0.6);
590 // gcor01[i] = TMath::Max(0.72 +TMath::Exp(3.36*(rsigma-1.2)),0.6);
591 // gcorp[i] = TMath::Max(TMath::Power((rsigma+0.5),1.5),1.2);
595 // for (Int_t i=3;i<10000;i++){
599 // Float_t amp = float(i);
600 // Float_t padlength =0.75;
601 // gnoise1 = 0.0004/padlength;
602 // Float_t nel = 0.268*amp;
603 // Float_t nprim = 0.155*amp;
604 // ggg1[i] = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.001*nel/(padlength*padlength))/nel;
605 // glandau1[i] = (2.+0.12*nprim)*0.5* (2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
606 // if (glandau1[i]>1) glandau1[i]=1;
607 // glandau1[i]*=padlength*padlength/12.;
611 // gnoise2 = 0.0004/padlength;
613 // nprim = 0.133*amp;
614 // ggg2[i] = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
615 // glandau2[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
616 // if (glandau2[i]>1) glandau2[i]=1;
617 // glandau2[i]*=padlength*padlength/12.;
622 // gnoise3 = 0.0004/padlength;
624 // nprim = 0.133*amp;
625 // ggg3[i] = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
626 // glandau3[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
627 // if (glandau3[i]>1) glandau3[i]=1;
628 // glandau3[i]*=padlength*padlength/12.;
636 // Int_t amp = int(TMath::Abs(cl->GetQ()));
638 // seed->SetErrorY2(1.);
642 // Float_t z = TMath::Abs(fParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
643 // Int_t ctype = cl->GetType();
644 // Float_t padlength= GetPadPitchLength(seed->GetRow());
645 // Double_t angle2 = seed->GetSnp()*seed->GetSnp();
646 // angle2 = angle2/(1-angle2);
648 // //cluster "quality"
649 // Int_t rsigmay = int(100.*cl->GetSigmaY2()/(seed->GetCurrentSigmaY2()));
652 // if (fSectors==fInnerSec){
653 // snoise2 = gnoise1;
654 // res = ggg1[amp]*z+glandau1[amp]*angle2;
655 // if (ctype==0) res *= gcor01[rsigmay];
658 // res*= gcorp[rsigmay];
662 // if (padlength<1.1){
663 // snoise2 = gnoise2;
664 // res = ggg2[amp]*z+glandau2[amp]*angle2;
665 // if (ctype==0) res *= gcor02[rsigmay];
668 // res*= gcorp[rsigmay];
672 // snoise2 = gnoise3;
673 // res = ggg3[amp]*z+glandau3[amp]*angle2;
674 // if (ctype==0) res *= gcor02[rsigmay];
677 // res*= gcorp[rsigmay];
684 // res*=2.4; // overestimate error 2 times
688 // if (res<2*snoise2)
691 // seed->SetErrorY2(res);
699 Double_t AliTPCtrackerMI::ErrZ2(AliTPCseed* seed, AliTPCclusterMI * cl){
702 // Use calibrated cluster error from OCDB
704 AliTPCClusterParam * clparam = AliTPCcalibDB::Instance()->GetClusterParam();
706 Float_t z = TMath::Abs(fParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
707 Int_t ctype = cl->GetType();
708 Int_t type = (cl->GetRow()<63) ? 0: (cl->GetRow()>126) ? 1:2;
710 Double_t angle2 = seed->GetSnp()*seed->GetSnp();
711 angle2 = seed->GetTgl()*seed->GetTgl()*(1+angle2/(1-angle2));
712 Double_t angle = TMath::Sqrt(TMath::Abs(angle2));
713 Double_t errz2 = clparam->GetError0Par(1,type, z,angle);
715 errz2+=0.5; // edge cluster
718 seed->SetErrorZ2(errz2);
724 // //seed->SetErrorY2(0.1);
726 // //calculate look-up table at the beginning
727 // static Bool_t ginit = kFALSE;
728 // static Float_t gnoise1,gnoise2,gnoise3;
729 // static Float_t ggg1[10000];
730 // static Float_t ggg2[10000];
731 // static Float_t ggg3[10000];
732 // static Float_t glandau1[10000];
733 // static Float_t glandau2[10000];
734 // static Float_t glandau3[10000];
736 // static Float_t gcor01[1000];
737 // static Float_t gcor02[1000];
738 // static Float_t gcorp[1000];
742 // if (ginit==kFALSE){
743 // for (Int_t i=1;i<1000;i++){
744 // Float_t rsigma = float(i)/100.;
745 // gcor02[i] = TMath::Max(0.81 +TMath::Exp(6.8*(rsigma-1.2)),0.6);
746 // gcor01[i] = TMath::Max(0.72 +TMath::Exp(2.04*(rsigma-1.2)),0.6);
747 // gcorp[i] = TMath::Max(TMath::Power((rsigma+0.5),1.5),1.2);
751 // for (Int_t i=3;i<10000;i++){
755 // Float_t amp = float(i);
756 // Float_t padlength =0.75;
757 // gnoise1 = 0.0004/padlength;
758 // Float_t nel = 0.268*amp;
759 // Float_t nprim = 0.155*amp;
760 // ggg1[i] = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.001*nel/(padlength*padlength))/nel;
761 // glandau1[i] = (2.+0.12*nprim)*0.5* (2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
762 // if (glandau1[i]>1) glandau1[i]=1;
763 // glandau1[i]*=padlength*padlength/12.;
767 // gnoise2 = 0.0004/padlength;
769 // nprim = 0.133*amp;
770 // ggg2[i] = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
771 // glandau2[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
772 // if (glandau2[i]>1) glandau2[i]=1;
773 // glandau2[i]*=padlength*padlength/12.;
778 // gnoise3 = 0.0004/padlength;
780 // nprim = 0.133*amp;
781 // ggg3[i] = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
782 // glandau3[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
783 // if (glandau3[i]>1) glandau3[i]=1;
784 // glandau3[i]*=padlength*padlength/12.;
792 // Int_t amp = int(TMath::Abs(cl->GetQ()));
794 // seed->SetErrorY2(1.);
798 // Float_t z = TMath::Abs(fParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
799 // Int_t ctype = cl->GetType();
800 // Float_t padlength= GetPadPitchLength(seed->GetRow());
802 // Double_t angle2 = seed->GetSnp()*seed->GetSnp();
803 // // if (angle2<0.6) angle2 = 0.6;
804 // angle2 = seed->GetTgl()*seed->GetTgl()*(1+angle2/(1-angle2));
806 // //cluster "quality"
807 // Int_t rsigmaz = int(100.*cl->GetSigmaZ2()/(seed->GetCurrentSigmaZ2()));
810 // if (fSectors==fInnerSec){
811 // snoise2 = gnoise1;
812 // res = ggg1[amp]*z+glandau1[amp]*angle2;
813 // if (ctype==0) res *= gcor01[rsigmaz];
816 // res*= gcorp[rsigmaz];
820 // if (padlength<1.1){
821 // snoise2 = gnoise2;
822 // res = ggg2[amp]*z+glandau2[amp]*angle2;
823 // if (ctype==0) res *= gcor02[rsigmaz];
826 // res*= gcorp[rsigmaz];
830 // snoise2 = gnoise3;
831 // res = ggg3[amp]*z+glandau3[amp]*angle2;
832 // if (ctype==0) res *= gcor02[rsigmaz];
835 // res*= gcorp[rsigmaz];
844 // if ((ctype<0) &&<70){
849 // if (res<2*snoise2)
851 // if (res>3) res =3;
852 // seed->SetErrorZ2(res);
860 void AliTPCtrackerMI::RotateToLocal(AliTPCseed *seed)
862 //rotate to track "local coordinata
863 Float_t x = seed->GetX();
864 Float_t y = seed->GetY();
865 Float_t ymax = x*TMath::Tan(0.5*fSectors->GetAlpha());
868 seed->SetRelativeSector((seed->GetRelativeSector()+1) % fN);
869 if (!seed->Rotate(fSectors->GetAlpha()))
871 } else if (y <-ymax) {
872 seed->SetRelativeSector((seed->GetRelativeSector()-1+fN) % fN);
873 if (!seed->Rotate(-fSectors->GetAlpha()))
881 //_____________________________________________________________________________
882 Double_t AliTPCtrackerMI::F1old(Double_t x1,Double_t y1,
883 Double_t x2,Double_t y2,
884 Double_t x3,Double_t y3)
886 //-----------------------------------------------------------------
887 // Initial approximation of the track curvature
888 //-----------------------------------------------------------------
889 Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
890 Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
891 (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
892 Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
893 (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
895 Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
896 if ( xr*xr+yr*yr<=0.00000000000001) return 100;
897 return -xr*yr/sqrt(xr*xr+yr*yr);
902 //_____________________________________________________________________________
903 Double_t AliTPCtrackerMI::F1(Double_t x1,Double_t y1,
904 Double_t x2,Double_t y2,
905 Double_t x3,Double_t y3)
907 //-----------------------------------------------------------------
908 // Initial approximation of the track curvature
909 //-----------------------------------------------------------------
915 Double_t det = x3*y2-x2*y3;
920 Double_t u = 0.5* (x2*(x2-x3)+y2*(y2-y3))/det;
921 Double_t x0 = x3*0.5-y3*u;
922 Double_t y0 = y3*0.5+x3*u;
923 Double_t c2 = 1/TMath::Sqrt(x0*x0+y0*y0);
929 Double_t AliTPCtrackerMI::F2(Double_t x1,Double_t y1,
930 Double_t x2,Double_t y2,
931 Double_t x3,Double_t y3)
933 //-----------------------------------------------------------------
934 // Initial approximation of the track curvature
935 //-----------------------------------------------------------------
941 Double_t det = x3*y2-x2*y3;
946 Double_t u = 0.5* (x2*(x2-x3)+y2*(y2-y3))/det;
947 Double_t x0 = x3*0.5-y3*u;
948 Double_t y0 = y3*0.5+x3*u;
949 Double_t c2 = 1/TMath::Sqrt(x0*x0+y0*y0);
958 //_____________________________________________________________________________
959 Double_t AliTPCtrackerMI::F2old(Double_t x1,Double_t y1,
960 Double_t x2,Double_t y2,
961 Double_t x3,Double_t y3)
963 //-----------------------------------------------------------------
964 // Initial approximation of the track curvature times center of curvature
965 //-----------------------------------------------------------------
966 Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
967 Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
968 (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
969 Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
970 (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
972 Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
974 return -a/(d*y1-b)*xr/sqrt(xr*xr+yr*yr);
977 //_____________________________________________________________________________
978 Double_t AliTPCtrackerMI::F3(Double_t x1,Double_t y1,
979 Double_t x2,Double_t y2,
980 Double_t z1,Double_t z2)
982 //-----------------------------------------------------------------
983 // Initial approximation of the tangent of the track dip angle
984 //-----------------------------------------------------------------
985 return (z1 - z2)/sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
989 Double_t AliTPCtrackerMI::F3n(Double_t x1,Double_t y1,
990 Double_t x2,Double_t y2,
991 Double_t z1,Double_t z2, Double_t c)
993 //-----------------------------------------------------------------
994 // Initial approximation of the tangent of the track dip angle
995 //-----------------------------------------------------------------
999 //angle1 = (z1-z2)*c/(TMath::ASin(c*x1-ni)-TMath::ASin(c*x2-ni));
1001 Double_t d = TMath::Sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
1002 if (TMath::Abs(d*c*0.5)>1) return 0;
1003 // Double_t angle2 = TMath::ASin(d*c*0.5);
1004 // Double_t angle2 = AliTPCFastMath::FastAsin(d*c*0.5);
1005 Double_t angle2 = (d*c*0.5>0.1)? TMath::ASin(d*c*0.5): AliTPCFastMath::FastAsin(d*c*0.5);
1007 angle2 = (z1-z2)*c/(angle2*2.);
1011 Bool_t AliTPCtrackerMI::GetProlongation(Double_t x1, Double_t x2, Double_t x[5], Double_t &y, Double_t &z)
1012 {//-----------------------------------------------------------------
1013 // This function find proloncation of a track to a reference plane x=x2.
1014 //-----------------------------------------------------------------
1018 if (TMath::Abs(x[4]*x1 - x[2]) >= 0.999) {
1022 Double_t c1=x[4]*x1 - x[2], r1=sqrt(1.- c1*c1);
1023 Double_t c2=x[4]*x2 - x[2], r2=sqrt(1.- c2*c2);
1027 Double_t dy = dx*(c1+c2)/(r1+r2);
1030 Double_t delta = x[4]*dx*(c1+c2)/(c1*r2 + c2*r1);
1032 if (TMath::Abs(delta)>0.01){
1033 dz = x[3]*TMath::ASin(delta)/x[4];
1035 dz = x[3]*AliTPCFastMath::FastAsin(delta)/x[4];
1038 //dz = x[3]*AliTPCFastMath::FastAsin(delta)/x[4];
1046 Int_t AliTPCtrackerMI::LoadClusters (TTree *tree)
1051 return LoadClusters();
1054 Int_t AliTPCtrackerMI::LoadClusters()
1057 // load clusters to the memory
1058 AliTPCClustersRow *clrow= new AliTPCClustersRow;
1059 clrow->SetClass("AliTPCclusterMI");
1061 clrow->GetArray()->ExpandCreateFast(10000);
1063 // TTree * tree = fClustersArray.GetTree();
1065 TTree * tree = fInput;
1066 TBranch * br = tree->GetBranch("Segment");
1067 br->SetAddress(&clrow);
1069 Int_t j=Int_t(tree->GetEntries());
1070 for (Int_t i=0; i<j; i++) {
1074 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 AliTPCRow * tpcrow=0;
1082 tpcrow = &(fInnerSec[sec%fkNIS][row]);
1086 tpcrow = &(fOuterSec[(sec-fkNIS*2)%fkNOS][row]);
1087 left = (sec-fkNIS*2)/fkNOS;
1090 tpcrow->SetN1(clrow->GetArray()->GetEntriesFast());
1091 tpcrow->SetClusters1(new AliTPCclusterMI[tpcrow->GetN1()]);
1092 for (Int_t i=0;i<tpcrow->GetN1();i++)
1093 tpcrow->SetCluster1(i, *(AliTPCclusterMI*)(clrow->GetArray()->At(i)));
1096 tpcrow->SetN2(clrow->GetArray()->GetEntriesFast());
1097 tpcrow->SetClusters2(new AliTPCclusterMI[tpcrow->GetN2()]);
1098 for (Int_t i=0;i<tpcrow->GetN2();i++)
1099 tpcrow->SetCluster2(i,*(AliTPCclusterMI*)(clrow->GetArray()->At(i)));
1110 void AliTPCtrackerMI::UnloadClusters()
1113 // unload clusters from the memory
1115 Int_t nrows = fOuterSec->GetNRows();
1116 for (Int_t sec = 0;sec<fkNOS;sec++)
1117 for (Int_t row = 0;row<nrows;row++){
1118 AliTPCRow* tpcrow = &(fOuterSec[sec%fkNOS][row]);
1120 // if (tpcrow->fClusters1) delete []tpcrow->fClusters1;
1121 // if (tpcrow->fClusters2) delete []tpcrow->fClusters2;
1123 tpcrow->ResetClusters();
1126 nrows = fInnerSec->GetNRows();
1127 for (Int_t sec = 0;sec<fkNIS;sec++)
1128 for (Int_t row = 0;row<nrows;row++){
1129 AliTPCRow* tpcrow = &(fInnerSec[sec%fkNIS][row]);
1131 // if (tpcrow->fClusters1) delete []tpcrow->fClusters1;
1132 //if (tpcrow->fClusters2) delete []tpcrow->fClusters2;
1134 tpcrow->ResetClusters();
1140 void AliTPCtrackerMI::Transform(AliTPCclusterMI * cluster){
1144 AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ;
1146 AliFatal("Tranformations not in calibDB");
1148 Double_t x[3]={cluster->GetRow(),cluster->GetPad(),cluster->GetTimeBin()};
1149 Int_t i[1]={cluster->GetDetector()};
1150 transform->Transform(x,i,0,1);
1151 if (!AliTPCReconstructor::GetRecoParam()->GetBYMirror()){
1152 if (cluster->GetDetector()%36>17){
1158 // in debug mode check the transformation
1160 if (AliTPCReconstructor::StreamLevel()>1) {
1161 TTreeSRedirector &cstream = *fDebugStreamer;
1162 cstream<<"Transform"<<
1169 cluster->SetX(x[0]);
1170 cluster->SetY(x[1]);
1171 cluster->SetZ(x[2]);
1177 //if (!fParam->IsGeoRead()) fParam->ReadGeoMatrices();
1178 TGeoHMatrix *mat = fParam->GetClusterMatrix(cluster->GetDetector());
1180 Double_t pos[3]= {cluster->GetX(),cluster->GetY(),cluster->GetZ()};
1181 Double_t posC[3]={cluster->GetX(),cluster->GetY(),cluster->GetZ()};
1182 if (mat) mat->LocalToMaster(pos,posC);
1184 // chack Loading of Geo matrices from GeoManager - TEMPORARY FIX
1186 cluster->SetX(posC[0]);
1187 cluster->SetY(posC[1]);
1188 cluster->SetZ(posC[2]);
1191 //_____________________________________________________________________________
1192 Int_t AliTPCtrackerMI::LoadOuterSectors() {
1193 //-----------------------------------------------------------------
1194 // This function fills outer TPC sectors with clusters.
1195 //-----------------------------------------------------------------
1196 Int_t nrows = fOuterSec->GetNRows();
1198 for (Int_t sec = 0;sec<fkNOS;sec++)
1199 for (Int_t row = 0;row<nrows;row++){
1200 AliTPCRow* tpcrow = &(fOuterSec[sec%fkNOS][row]);
1201 Int_t sec2 = sec+2*fkNIS;
1203 Int_t ncl = tpcrow->GetN1();
1205 AliTPCclusterMI *c= (tpcrow->GetCluster1(ncl));
1206 index=(((sec2<<8)+row)<<16)+ncl;
1207 tpcrow->InsertCluster(c,index);
1210 ncl = tpcrow->GetN2();
1212 AliTPCclusterMI *c= (tpcrow->GetCluster2(ncl));
1213 index=((((sec2+fkNOS)<<8)+row)<<16)+ncl;
1214 tpcrow->InsertCluster(c,index);
1217 // write indexes for fast acces
1219 for (Int_t i=0;i<510;i++)
1220 tpcrow->SetFastCluster(i,-1);
1221 for (Int_t i=0;i<tpcrow->GetN();i++){
1222 Int_t zi = Int_t((*tpcrow)[i]->GetZ()+255.);
1223 tpcrow->SetFastCluster(zi,i); // write index
1226 for (Int_t i=0;i<510;i++){
1227 if (tpcrow->GetFastCluster(i)<0)
1228 tpcrow->SetFastCluster(i,last);
1230 last = tpcrow->GetFastCluster(i);
1239 //_____________________________________________________________________________
1240 Int_t AliTPCtrackerMI::LoadInnerSectors() {
1241 //-----------------------------------------------------------------
1242 // This function fills inner TPC sectors with clusters.
1243 //-----------------------------------------------------------------
1244 Int_t nrows = fInnerSec->GetNRows();
1246 for (Int_t sec = 0;sec<fkNIS;sec++)
1247 for (Int_t row = 0;row<nrows;row++){
1248 AliTPCRow* tpcrow = &(fInnerSec[sec%fkNIS][row]);
1251 Int_t ncl = tpcrow->GetN1();
1253 AliTPCclusterMI *c= (tpcrow->GetCluster1(ncl));
1254 index=(((sec<<8)+row)<<16)+ncl;
1255 tpcrow->InsertCluster(c,index);
1258 ncl = tpcrow->GetN2();
1260 AliTPCclusterMI *c= (tpcrow->GetCluster2(ncl));
1261 index=((((sec+fkNIS)<<8)+row)<<16)+ncl;
1262 tpcrow->InsertCluster(c,index);
1265 // write indexes for fast acces
1267 for (Int_t i=0;i<510;i++)
1268 tpcrow->SetFastCluster(i,-1);
1269 for (Int_t i=0;i<tpcrow->GetN();i++){
1270 Int_t zi = Int_t((*tpcrow)[i]->GetZ()+255.);
1271 tpcrow->SetFastCluster(zi,i); // write index
1274 for (Int_t i=0;i<510;i++){
1275 if (tpcrow->GetFastCluster(i)<0)
1276 tpcrow->SetFastCluster(i,last);
1278 last = tpcrow->GetFastCluster(i);
1290 //_________________________________________________________________________
1291 AliTPCclusterMI *AliTPCtrackerMI::GetClusterMI(Int_t index) const {
1292 //--------------------------------------------------------------------
1293 // Return pointer to a given cluster
1294 //--------------------------------------------------------------------
1295 if (index<0) return 0; // no cluster
1296 Int_t sec=(index&0xff000000)>>24;
1297 Int_t row=(index&0x00ff0000)>>16;
1298 Int_t ncl=(index&0x00007fff)>>00;
1300 const AliTPCRow * tpcrow=0;
1301 AliTPCclusterMI * clrow =0;
1303 if (sec<0 || sec>=fkNIS*4) {
1304 AliWarning(Form("Wrong sector %d",sec));
1309 tpcrow = &(fInnerSec[sec%fkNIS][row]);
1310 if (tpcrow==0) return 0;
1313 if (tpcrow->GetN1()<=ncl) return 0;
1314 clrow = tpcrow->GetClusters1();
1317 if (tpcrow->GetN2()<=ncl) return 0;
1318 clrow = tpcrow->GetClusters2();
1322 tpcrow = &(fOuterSec[(sec-fkNIS*2)%fkNOS][row]);
1323 if (tpcrow==0) return 0;
1325 if (sec-2*fkNIS<fkNOS) {
1326 if (tpcrow->GetN1()<=ncl) return 0;
1327 clrow = tpcrow->GetClusters1();
1330 if (tpcrow->GetN2()<=ncl) return 0;
1331 clrow = tpcrow->GetClusters2();
1335 return &(clrow[ncl]);
1341 Int_t AliTPCtrackerMI::FollowToNext(AliTPCseed& t, Int_t nr) {
1342 //-----------------------------------------------------------------
1343 // This function tries to find a track prolongation to next pad row
1344 //-----------------------------------------------------------------
1346 Double_t x= GetXrow(nr), ymax=GetMaxY(nr);
1347 AliTPCclusterMI *cl=0;
1348 Int_t tpcindex= t.GetClusterIndex2(nr);
1350 // update current shape info every 5 pad-row
1351 // if ( (nr%5==0) || t.GetNumberOfClusters()<2 || (t.fCurrentSigmaY2<0.0001) ){
1355 if (fIteration>0 && tpcindex>=-1){ //if we have already clusters
1357 if (tpcindex==-1) return 0; //track in dead zone
1359 cl = t.GetClusterPointer(nr);
1360 if ( (cl==0) ) cl = GetClusterMI(tpcindex);
1361 t.SetCurrentClusterIndex1(tpcindex);
1364 Int_t relativesector = ((tpcindex&0xff000000)>>24)%18; // if previously accepted cluster in different sector
1365 Float_t angle = relativesector*fSectors->GetAlpha()+fSectors->GetAlphaShift();
1367 if (angle<-TMath::Pi()) angle += 2*TMath::Pi();
1368 if (angle>=TMath::Pi()) angle -= 2*TMath::Pi();
1370 if (TMath::Abs(angle-t.GetAlpha())>0.001){
1371 Double_t rotation = angle-t.GetAlpha();
1372 t.SetRelativeSector(relativesector);
1373 if (!t.Rotate(rotation)) return 0;
1375 if (!t.PropagateTo(x)) return 0;
1377 t.SetCurrentCluster(cl);
1379 Int_t accept = AcceptCluster(&t,t.GetCurrentCluster());
1380 if ((tpcindex&0x8000)==0) accept =0;
1382 //if founded cluster is acceptible
1383 if (cl->IsUsed(11)) { // id cluster is shared inrease uncertainty
1384 t.SetErrorY2(t.GetErrorY2()+0.03);
1385 t.SetErrorZ2(t.GetErrorZ2()+0.03);
1386 t.SetErrorY2(t.GetErrorY2()*3);
1387 t.SetErrorZ2(t.GetErrorZ2()*3);
1389 t.SetNFoundable(t.GetNFoundable()+1);
1390 UpdateTrack(&t,accept);
1395 if (TMath::Abs(t.GetSnp())>AliTPCReconstructor::GetMaxSnpTracker()) return 0; // cut on angle
1397 // not look for new cluster during refitting
1398 t.SetNFoundable(t.GetNFoundable()+1);
1403 // if (TMath::Abs(t.GetSnp())>0.95 || TMath::Abs(x*t.GetC()-t.GetEta())>0.95) return 0;// patch 28 fev 06
1404 Double_t y=t.GetYat(x);
1405 if (TMath::Abs(y)>ymax){
1407 t.SetRelativeSector((t.GetRelativeSector()+1) % fN);
1408 if (!t.Rotate(fSectors->GetAlpha()))
1410 } else if (y <-ymax) {
1411 t.SetRelativeSector((t.GetRelativeSector()-1+fN) % fN);
1412 if (!t.Rotate(-fSectors->GetAlpha()))
1418 if (!t.PropagateTo(x)) {
1419 if (fIteration==0) t.SetRemoval(10);
1423 Double_t z=t.GetZ();
1426 if (!IsActive(t.GetRelativeSector(),nr)) {
1428 t.SetClusterIndex2(nr,-1);
1431 //AliInfo(Form("A - Sector%d phi %f - alpha %f", t.fRelativeSector,y/x, t.GetAlpha()));
1432 Bool_t isActive = IsActive(t.GetRelativeSector(),nr);
1433 Bool_t isActive2 = (nr>=fInnerSec->GetNRows()) ? fOuterSec[t.GetRelativeSector()][nr-fInnerSec->GetNRows()].GetN()>0:fInnerSec[t.GetRelativeSector()][nr].GetN()>0;
1435 if (!isActive || !isActive2) return 0;
1437 const AliTPCRow &krow=GetRow(t.GetRelativeSector(),nr);
1438 if ( (t.GetSigmaY2()<0) || t.GetSigmaZ2()<0) return 0;
1440 Double_t roadz = 1.;
1442 if (TMath::Abs(TMath::Abs(y)-ymax)<krow.GetDeadZone()){
1444 t.SetClusterIndex2(nr,-1);
1449 if (TMath::Abs(z)<(AliTPCReconstructor::GetCtgRange()*x+10) && TMath::Abs(z)<fParam->GetZLength(0) && (TMath::Abs(t.GetSnp())<AliTPCReconstructor::GetMaxSnpTracker()))
1450 t.SetNFoundable(t.GetNFoundable()+1);
1456 // cl = krow.FindNearest2(y+10.,z,roady,roadz,index);
1457 cl = krow.FindNearest2(y,z,roady,roadz,index);
1458 if (cl) t.SetCurrentClusterIndex1(krow.GetIndex(index));
1461 t.SetCurrentCluster(cl);
1463 if (fIteration==2&&cl->IsUsed(10)) return 0;
1464 Int_t accept = AcceptCluster(&t,t.GetCurrentCluster());
1465 if (fIteration==2&&cl->IsUsed(11)) {
1466 t.SetErrorY2(t.GetErrorY2()+0.03);
1467 t.SetErrorZ2(t.GetErrorZ2()+0.03);
1468 t.SetErrorY2(t.GetErrorY2()*3);
1469 t.SetErrorZ2(t.GetErrorZ2()*3);
1472 if (t.fCurrentCluster->IsUsed(10)){
1477 if (t.fNShared>0.7*t.GetNumberOfClusters()) {
1483 if (accept<3) UpdateTrack(&t,accept);
1486 if ( fIteration==0 && t.GetNFoundable()*0.5 > t.GetNumberOfClusters()) t.SetRemoval(10);
1492 Int_t AliTPCtrackerMI::FollowToNextFast(AliTPCseed& t, Int_t nr) {
1493 //-----------------------------------------------------------------
1494 // This function tries to find a track prolongation to next pad row
1495 //-----------------------------------------------------------------
1497 Double_t x= GetXrow(nr), ymax=GetMaxY(nr);
1499 if (!t.GetProlongation(x,y,z)) {
1505 if (TMath::Abs(y)>ymax){
1508 t.SetRelativeSector((t.GetRelativeSector()+1) % fN);
1509 if (!t.Rotate(fSectors->GetAlpha()))
1511 } else if (y <-ymax) {
1512 t.SetRelativeSector((t.GetRelativeSector()-1+fN) % fN);
1513 if (!t.Rotate(-fSectors->GetAlpha()))
1516 if (!t.PropagateTo(x)) {
1519 t.GetProlongation(x,y,z);
1522 // update current shape info every 3 pad-row
1523 if ( (nr%6==0) || t.GetNumberOfClusters()<2 || (t.GetCurrentSigmaY2()<0.0001) ){
1524 // t.fCurrentSigmaY = GetSigmaY(&t);
1525 //t.fCurrentSigmaZ = GetSigmaZ(&t);
1529 AliTPCclusterMI *cl=0;
1534 const AliTPCRow &krow=GetRow(t.GetRelativeSector(),nr);
1535 if ( (t.GetSigmaY2()<0) || t.GetSigmaZ2()<0) return 0;
1537 Double_t roadz = 1.;
1540 if (TMath::Abs(TMath::Abs(y)-ymax)<krow.GetDeadZone()){
1542 t.SetClusterIndex2(row,-1);
1547 if (TMath::Abs(z)>(AliTPCReconstructor::GetCtgRange()*x+10)) t.SetClusterIndex2(row,-1);
1551 if ((cl==0)&&(krow)) {
1552 // cl = krow.FindNearest2(y+10,z,roady,roadz,index);
1553 cl = krow.FindNearest2(y,z,roady,roadz,index);
1555 if (cl) t.SetCurrentClusterIndex1(krow.GetIndex(index));
1559 t.SetCurrentCluster(cl);
1560 // Int_t accept = AcceptCluster(&t,t.fCurrentCluster);
1562 t.SetClusterIndex2(row,index);
1563 t.SetClusterPointer(row, cl);
1571 //_________________________________________________________________________
1572 Bool_t AliTPCtrackerMI::GetTrackPoint(Int_t index, AliTrackPoint &p ) const
1574 // Get track space point by index
1575 // return false in case the cluster doesn't exist
1576 AliTPCclusterMI *cl = GetClusterMI(index);
1577 if (!cl) return kFALSE;
1578 Int_t sector = (index&0xff000000)>>24;
1579 // Int_t row = (index&0x00ff0000)>>16;
1581 // xyz[0] = fParam->GetPadRowRadii(sector,row);
1582 xyz[0] = cl->GetX();
1583 xyz[1] = cl->GetY();
1584 xyz[2] = cl->GetZ();
1586 fParam->AdjustCosSin(sector,cos,sin);
1587 Float_t x = cos*xyz[0]-sin*xyz[1];
1588 Float_t y = cos*xyz[1]+sin*xyz[0];
1590 Float_t sigmaY2 = 0.027*cl->GetSigmaY2();
1591 if (sector < fParam->GetNInnerSector()) sigmaY2 *= 2.07;
1592 Float_t sigmaZ2 = 0.066*cl->GetSigmaZ2();
1593 if (sector < fParam->GetNInnerSector()) sigmaZ2 *= 1.77;
1594 cov[0] = sin*sin*sigmaY2;
1595 cov[1] = -sin*cos*sigmaY2;
1597 cov[3] = cos*cos*sigmaY2;
1600 p.SetXYZ(x,y,xyz[2],cov);
1601 AliGeomManager::ELayerID iLayer;
1603 if (sector < fParam->GetNInnerSector()) {
1604 iLayer = AliGeomManager::kTPC1;
1608 iLayer = AliGeomManager::kTPC2;
1609 idet = sector - fParam->GetNInnerSector();
1611 UShort_t volid = AliGeomManager::LayerToVolUID(iLayer,idet);
1612 p.SetVolumeID(volid);
1618 Int_t AliTPCtrackerMI::UpdateClusters(AliTPCseed& t, Int_t nr) {
1619 //-----------------------------------------------------------------
1620 // This function tries to find a track prolongation to next pad row
1621 //-----------------------------------------------------------------
1622 t.SetCurrentCluster(0);
1623 t.SetCurrentClusterIndex1(0);
1625 Double_t xt=t.GetX();
1626 Int_t row = GetRowNumber(xt)-1;
1627 Double_t ymax= GetMaxY(nr);
1629 if (row < nr) return 1; // don't prolongate if not information until now -
1630 // if (TMath::Abs(t.GetSnp())>0.9 && t.GetNumberOfClusters()>40. && fIteration!=2) {
1632 // return 0; // not prolongate strongly inclined tracks
1634 // if (TMath::Abs(t.GetSnp())>0.95) {
1636 // return 0; // not prolongate strongly inclined tracks
1637 // }// patch 28 fev 06
1639 Double_t x= GetXrow(nr);
1641 //t.PropagateTo(x+0.02);
1642 //t.PropagateTo(x+0.01);
1643 if (!t.PropagateTo(x)){
1650 if (TMath::Abs(y)>ymax){
1652 t.SetRelativeSector((t.GetRelativeSector()+1) % fN);
1653 if (!t.Rotate(fSectors->GetAlpha()))
1655 } else if (y <-ymax) {
1656 t.SetRelativeSector((t.GetRelativeSector()-1+fN) % fN);
1657 if (!t.Rotate(-fSectors->GetAlpha()))
1660 // if (!t.PropagateTo(x)){
1667 if (TMath::Abs(t.GetSnp())>AliTPCReconstructor::GetMaxSnpTracker()) return 0;
1669 if (!IsActive(t.GetRelativeSector(),nr)) {
1671 t.SetClusterIndex2(nr,-1);
1674 //AliInfo(Form("A - Sector%d phi %f - alpha %f", t.fRelativeSector,y/x, t.GetAlpha()));
1676 AliTPCRow &krow=GetRow(t.GetRelativeSector(),nr);
1678 if (TMath::Abs(TMath::Abs(y)-ymax)<krow.GetDeadZone()){
1680 t.SetClusterIndex2(nr,-1);
1685 if (TMath::Abs(t.GetZ())<(AliTPCReconstructor::GetCtgRange()*t.GetX()+10) && (TMath::Abs(t.GetSnp())<AliTPCReconstructor::GetMaxSnpTracker())) t.SetNFoundable(t.GetNFoundable()+1);
1691 if ( (nr%6==0) || t.GetNumberOfClusters()<2){
1692 // t.fCurrentSigmaY = GetSigmaY(&t);
1693 //t.fCurrentSigmaZ = GetSigmaZ(&t);
1697 AliTPCclusterMI *cl=0;
1700 Double_t roady = 1.;
1701 Double_t roadz = 1.;
1705 index = t.GetClusterIndex2(nr);
1706 if ( (index>0) && (index&0x8000)==0){
1707 cl = t.GetClusterPointer(nr);
1708 if ( (cl==0) && (index>0)) cl = GetClusterMI(index);
1709 t.SetCurrentClusterIndex1(index);
1711 t.SetCurrentCluster(cl);
1717 // if (index<0) return 0;
1718 UInt_t uindex = TMath::Abs(index);
1721 //cl = krow.FindNearest2(y+10,z,roady,roadz,uindex);
1722 cl = krow.FindNearest2(y,z,roady,roadz,uindex);
1725 if (cl) t.SetCurrentClusterIndex1(krow.GetIndex(uindex));
1726 t.SetCurrentCluster(cl);
1732 Int_t AliTPCtrackerMI::FollowToNextCluster(AliTPCseed & t, Int_t nr) {
1733 //-----------------------------------------------------------------
1734 // This function tries to find a track prolongation to next pad row
1735 //-----------------------------------------------------------------
1737 //update error according neighborhoud
1739 if (t.GetCurrentCluster()) {
1741 Int_t accept = AcceptCluster(&t,t.GetCurrentCluster());
1743 if (t.GetCurrentCluster()->IsUsed(10)){
1748 t.SetNShared(t.GetNShared()+1);
1749 if (t.GetNShared()>0.7*t.GetNumberOfClusters()) {
1754 if (fIteration>0) accept = 0;
1755 if (accept<3) UpdateTrack(&t,accept);
1759 if ( ( (t.GetSigmaY2()+t.GetSigmaZ2())>0.16)&& t.GetNumberOfClusters()>18) t.SetRemoval(10);
1760 if ( t.GetChi2()/t.GetNumberOfClusters()>6 &&t.GetNumberOfClusters()>18) t.SetRemoval(10);
1762 if (( (t.GetNFoundable()*0.5 > t.GetNumberOfClusters()) || t.GetNoCluster()>15)) t.SetRemoval(10);
1770 //_____________________________________________________________________________
1771 Int_t AliTPCtrackerMI::FollowProlongation(AliTPCseed& t, Int_t rf, Int_t step) {
1772 //-----------------------------------------------------------------
1773 // This function tries to find a track prolongation.
1774 //-----------------------------------------------------------------
1775 Double_t xt=t.GetX();
1777 Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
1778 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
1779 if (alpha < 0. ) alpha += 2.*TMath::Pi();
1781 t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN);
1783 Int_t first = GetRowNumber(xt)-1;
1784 for (Int_t nr= first; nr>=rf; nr-=step) {
1786 if (t.GetKinkIndexes()[0]>0){
1787 for (Int_t i=0;i<3;i++){
1788 Int_t index = t.GetKinkIndexes()[i];
1789 if (index==0) break;
1790 if (index<0) continue;
1792 AliKink * kink = (AliKink*)fEvent->GetKink(index-1);
1794 printf("PROBLEM\n");
1797 Int_t kinkrow = kink->GetTPCRow0()+2+Int_t(0.5/(0.05+kink->GetAngle(2)));
1799 AliExternalTrackParam paramd(t);
1800 kink->SetDaughter(paramd);
1801 kink->SetStatus(2,5);
1808 if (nr==80) t.UpdateReference();
1809 if (nr<fInnerSec->GetNRows())
1810 fSectors = fInnerSec;
1812 fSectors = fOuterSec;
1813 if (FollowToNext(t,nr)==0)
1822 //_____________________________________________________________________________
1823 Int_t AliTPCtrackerMI::FollowProlongationFast(AliTPCseed& t, Int_t rf, Int_t step) {
1824 //-----------------------------------------------------------------
1825 // This function tries to find a track prolongation.
1826 //-----------------------------------------------------------------
1827 Double_t xt=t.GetX();
1829 Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
1830 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
1831 if (alpha < 0. ) alpha += 2.*TMath::Pi();
1832 t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN);
1834 for (Int_t nr=GetRowNumber(xt)-1; nr>=rf; nr-=step) {
1836 if (FollowToNextFast(t,nr)==0)
1837 if (!t.IsActive()) return 0;
1847 Int_t AliTPCtrackerMI::FollowBackProlongation(AliTPCseed& t, Int_t rf) {
1848 //-----------------------------------------------------------------
1849 // This function tries to find a track prolongation.
1850 //-----------------------------------------------------------------
1852 Double_t xt=t.GetX();
1853 Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
1854 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
1855 if (alpha < 0. ) alpha += 2.*TMath::Pi();
1856 t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN);
1858 Int_t first = t.GetFirstPoint();
1859 if (first<GetRowNumber(xt)+1) first = GetRowNumber(xt)+1;
1861 if (first<0) first=0;
1862 for (Int_t nr=first; nr<=rf; nr++) {
1863 // if ( (TMath::Abs(t.GetSnp())>0.95)) break;//patch 28 fev 06
1864 if (t.GetKinkIndexes()[0]<0){
1865 for (Int_t i=0;i<3;i++){
1866 Int_t index = t.GetKinkIndexes()[i];
1867 if (index==0) break;
1868 if (index>0) continue;
1869 index = TMath::Abs(index);
1870 AliKink * kink = (AliKink*)fEvent->GetKink(index-1);
1872 printf("PROBLEM\n");
1875 Int_t kinkrow = kink->GetTPCRow0()-2-Int_t(0.5/(0.05+kink->GetAngle(2)));
1877 AliExternalTrackParam paramm(t);
1878 kink->SetMother(paramm);
1879 kink->SetStatus(2,1);
1886 if (nr<fInnerSec->GetNRows())
1887 fSectors = fInnerSec;
1889 fSectors = fOuterSec;
1900 Float_t AliTPCtrackerMI::OverlapFactor(AliTPCseed * s1, AliTPCseed * s2, Int_t &sum1, Int_t & sum2)
1908 Float_t dz2 =(s1->GetZ() - s2->GetZ());
1911 Float_t dy2 =TMath::Abs((s1->GetY() - s2->GetY()));
1913 Float_t distance = TMath::Sqrt(dz2+dy2);
1914 if (distance>4.) return 0; // if there are far away - not overlap - to reduce combinatorics
1917 Int_t firstpoint = TMath::Min(s1->GetFirstPoint(),s2->GetFirstPoint());
1918 Int_t lastpoint = TMath::Max(s1->GetLastPoint(),s2->GetLastPoint());
1923 if (firstpoint>lastpoint) {
1924 firstpoint =lastpoint;
1929 for (Int_t i=firstpoint-1;i<lastpoint+1;i++){
1930 if (s1->GetClusterIndex2(i)>0) sum1++;
1931 if (s2->GetClusterIndex2(i)>0) sum2++;
1932 if (s1->GetClusterIndex2(i)==s2->GetClusterIndex2(i) && s1->GetClusterIndex2(i)>0) {
1936 if (sum<5) return 0;
1938 Float_t summin = TMath::Min(sum1+1,sum2+1);
1939 Float_t ratio = (sum+1)/Float_t(summin);
1943 void AliTPCtrackerMI::SignShared(AliTPCseed * s1, AliTPCseed * s2)
1947 Float_t thetaCut = 0.2;//+10.*TMath::Sqrt(s1->GetSigmaTglZ()+ s2->GetSigmaTglZ());
1948 if (TMath::Abs(s1->GetTgl()-s2->GetTgl())>thetaCut) return;
1949 Float_t minCl = TMath::Min(s1->GetNumberOfClusters(),s2->GetNumberOfClusters());
1950 Int_t cutN0 = TMath::Max(5,TMath::Nint(0.1*minCl));
1955 //Int_t firstpoint = TMath::Max(s1->GetFirstPoint(),s2->GetFirstPoint());
1956 //Int_t lastpoint = TMath::Min(s1->GetLastPoint(),s2->GetLastPoint());
1957 Int_t firstpoint = 0;
1958 Int_t lastpoint = 160;
1960 // if (firstpoint>=lastpoint-5) return;;
1962 for (Int_t i=firstpoint;i<lastpoint;i++){
1963 // if ( (s1->GetClusterIndex2(i)&0xFFFF8FFF)==(s2->GetClusterIndex2(i)&0xFFFF8FFF) && s1->GetClusterIndex2(i)>0) {
1964 if ( (s1->GetClusterIndex2(i))==(s2->GetClusterIndex2(i)) && s1->GetClusterIndex2(i)>0) {
1966 s1->SetSharedMapBit(i, kTRUE);
1967 s2->SetSharedMapBit(i, kTRUE);
1969 if (s1->GetClusterIndex2(i)>0)
1970 s1->SetClusterMapBit(i, kTRUE);
1972 if (sumshared>cutN0){
1975 for (Int_t i=firstpoint;i<lastpoint;i++){
1976 // if ( (s1->GetClusterIndex2(i)&0xFFFF8FFF)==(s2->GetClusterIndex2(i)&0xFFFF8FFF) && s1->GetClusterIndex2(i)>0) {
1977 if ( (s1->GetClusterIndex2(i))==(s2->GetClusterIndex2(i)) && s1->GetClusterIndex2(i)>0) {
1978 AliTPCTrackerPoint *p1 = s1->GetTrackPoint(i);
1979 AliTPCTrackerPoint *p2 = s2->GetTrackPoint(i);;
1980 if (s1->IsActive()&&s2->IsActive()){
1981 p1->SetShared(kTRUE);
1982 p2->SetShared(kTRUE);
1988 if (sumshared>cutN0){
1989 for (Int_t i=0;i<4;i++){
1990 if (s1->GetOverlapLabel(3*i)==0){
1991 s1->SetOverlapLabel(3*i, s2->GetLabel());
1992 s1->SetOverlapLabel(3*i+1,sumshared);
1993 s1->SetOverlapLabel(3*i+2,s2->GetUniqueID());
1997 for (Int_t i=0;i<4;i++){
1998 if (s2->GetOverlapLabel(3*i)==0){
1999 s2->SetOverlapLabel(3*i, s1->GetLabel());
2000 s2->SetOverlapLabel(3*i+1,sumshared);
2001 s2->SetOverlapLabel(3*i+2,s1->GetUniqueID());
2008 void AliTPCtrackerMI::SignShared(TObjArray * arr)
2011 //sort trackss according sectors
2013 for (Int_t i=0; i<arr->GetEntriesFast(); i++) {
2014 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2016 //if (pt) RotateToLocal(pt);
2020 arr->Sort(); // sorting according relative sectors
2021 arr->Expand(arr->GetEntries());
2024 Int_t nseed=arr->GetEntriesFast();
2025 for (Int_t i=0; i<nseed; i++) {
2026 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2028 for (Int_t j=0;j<=12;j++){
2029 pt->SetOverlapLabel(j,0);
2032 for (Int_t i=0; i<nseed; i++) {
2033 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2035 if (pt->GetRemoval()>10) continue;
2036 for (Int_t j=i+1; j<nseed; j++){
2037 AliTPCseed *pt2=(AliTPCseed*)arr->UncheckedAt(j);
2038 if (TMath::Abs(pt->GetRelativeSector()-pt2->GetRelativeSector())>1) continue;
2040 if (pt2->GetRemoval()<=10) {
2041 //if ( TMath::Abs(pt->GetRelativeSector()-pt2->GetRelativeSector())>0) break;
2049 void AliTPCtrackerMI::SortTracks(TObjArray * arr, Int_t mode) const
2052 //sort tracks in array according mode criteria
2053 Int_t nseed = arr->GetEntriesFast();
2054 for (Int_t i=0; i<nseed; i++) {
2055 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2066 void AliTPCtrackerMI::RemoveUsed2(TObjArray * arr, Float_t factor1, Float_t factor2, Int_t minimal)
2069 // Loop over all tracks and remove overlaped tracks (with lower quality)
2071 // 1. Unsign clusters
2072 // 2. Sort tracks according quality
2073 // Quality is defined by the number of cluster between first and last points
2075 // 3. Loop over tracks - decreasing quality order
2076 // a.) remove - If the fraction of shared cluster less than factor (1- n or 2)
2077 // b.) remove - If the minimal number of clusters less than minimal and not ITS
2078 // c.) if track accepted - sign clusters
2080 //Called in - AliTPCtrackerMI::Clusters2Tracks()
2081 // - AliTPCtrackerMI::PropagateBack()
2082 // - AliTPCtrackerMI::RefitInward()
2088 Int_t nseed = arr->GetEntriesFast();
2089 Float_t * quality = new Float_t[nseed];
2090 Int_t * indexes = new Int_t[nseed];
2094 for (Int_t i=0; i<nseed; i++) {
2095 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2100 pt->UpdatePoints(); //select first last max dens points
2101 Float_t * points = pt->GetPoints();
2102 if (points[3]<0.8) quality[i] =-1;
2103 quality[i] = (points[2]-points[0])+pt->GetNumberOfClusters();
2104 //prefer high momenta tracks if overlaps
2105 quality[i] *= TMath::Sqrt(TMath::Abs(pt->Pt())+0.5);
2107 TMath::Sort(nseed,quality,indexes);
2110 for (Int_t itrack=0; itrack<nseed; itrack++) {
2111 Int_t trackindex = indexes[itrack];
2112 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(trackindex);
2115 if (quality[trackindex]<0){
2117 delete arr->RemoveAt(trackindex);
2120 arr->RemoveAt(trackindex);
2126 Int_t first = Int_t(pt->GetPoints()[0]);
2127 Int_t last = Int_t(pt->GetPoints()[2]);
2128 Double_t factor = (pt->GetBConstrain()) ? factor1: factor2;
2130 Int_t found,foundable,shared;
2131 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
2132 // pt->GetClusterStatistic(0,160, found, foundable,shared,kFALSE);
2133 Bool_t itsgold =kFALSE;
2136 if (pt->GetESD()->GetITSclusters(dummy)>4) itsgold= kTRUE;
2140 if (Float_t(shared+1)/Float_t(found+1)>factor){
2141 if (pt->GetKinkIndexes()[0]!=0) continue; //don't remove tracks - part of the kinks
2142 delete arr->RemoveAt(trackindex);
2145 if (pt->GetNumberOfClusters()<50&&(found-0.5*shared)<minimal){ //remove short tracks
2146 if (pt->GetKinkIndexes()[0]!=0) continue; //don't remove tracks - part of the kinks
2147 delete arr->RemoveAt(trackindex);
2153 //if (sharedfactor>0.4) continue;
2154 if (pt->GetKinkIndexes()[0]>0) continue;
2155 //Remove tracks with undefined properties - seems
2156 if (pt->GetSigmaY2()<kAlmost0) continue; // ? what is the origin ?
2158 for (Int_t i=first; i<last; i++) {
2159 Int_t index=pt->GetClusterIndex2(i);
2160 // if (index<0 || index&0x8000 ) continue;
2161 if (index<0 || index&0x8000 ) continue;
2162 AliTPCclusterMI *c= pt->GetClusterPointer(i);
2169 Info("RemoveUsed","\n*****\nNumber of good tracks after shared removal\t%d\n",fNtracks);
2175 void AliTPCtrackerMI::UnsignClusters()
2178 // loop over all clusters and unsign them
2181 for (Int_t sec=0;sec<fkNIS;sec++){
2182 for (Int_t row=0;row<fInnerSec->GetNRows();row++){
2183 AliTPCclusterMI *cl = fInnerSec[sec][row].GetClusters1();
2184 for (Int_t icl =0;icl< fInnerSec[sec][row].GetN1();icl++)
2185 // if (cl[icl].IsUsed(10))
2187 cl = fInnerSec[sec][row].GetClusters2();
2188 for (Int_t icl =0;icl< fInnerSec[sec][row].GetN2();icl++)
2189 //if (cl[icl].IsUsed(10))
2194 for (Int_t sec=0;sec<fkNOS;sec++){
2195 for (Int_t row=0;row<fOuterSec->GetNRows();row++){
2196 AliTPCclusterMI *cl = fOuterSec[sec][row].GetClusters1();
2197 for (Int_t icl =0;icl< fOuterSec[sec][row].GetN1();icl++)
2198 //if (cl[icl].IsUsed(10))
2200 cl = fOuterSec[sec][row].GetClusters2();
2201 for (Int_t icl =0;icl< fOuterSec[sec][row].GetN2();icl++)
2202 //if (cl[icl].IsUsed(10))
2211 void AliTPCtrackerMI::SignClusters(TObjArray * arr, Float_t fnumber, Float_t fdensity)
2214 //sign clusters to be "used"
2216 // snumber and sdensity sign number of sigmas - bellow mean value to be accepted
2217 // loop over "primaries"
2231 Int_t nseed = arr->GetEntriesFast();
2232 for (Int_t i=0; i<nseed; i++) {
2233 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2237 if (!(pt->IsActive())) continue;
2238 Float_t dens = pt->GetNumberOfClusters()/Float_t(pt->GetNFoundable());
2239 if ( (dens>0.7) && (pt->GetNumberOfClusters()>70)){
2241 sumdens2+= dens*dens;
2242 sumn += pt->GetNumberOfClusters();
2243 sumn2 += pt->GetNumberOfClusters()*pt->GetNumberOfClusters();
2244 Float_t chi2 = pt->GetChi2()/pt->GetNumberOfClusters();
2247 sumchi2 +=chi2*chi2;
2252 Float_t mdensity = 0.9;
2253 Float_t meann = 130;
2254 Float_t meanchi = 1;
2255 Float_t sdensity = 0.1;
2256 Float_t smeann = 10;
2257 Float_t smeanchi =0.4;
2261 mdensity = sumdens/sum;
2263 meanchi = sumchi/sum;
2265 sdensity = sumdens2/sum-mdensity*mdensity;
2267 sdensity = TMath::Sqrt(sdensity);
2271 smeann = sumn2/sum-meann*meann;
2273 smeann = TMath::Sqrt(smeann);
2277 smeanchi = sumchi2/sum - meanchi*meanchi;
2279 smeanchi = TMath::Sqrt(smeanchi);
2285 //REMOVE SHORT DELTAS or tracks going out of sensitive volume of TPC
2287 for (Int_t i=0; i<nseed; i++) {
2288 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2292 if (pt->GetBSigned()) continue;
2293 if (pt->GetBConstrain()) continue;
2294 //if (!(pt->IsActive())) continue;
2296 Int_t found,foundable,shared;
2297 pt->GetClusterStatistic(0,160,found, foundable,shared);
2298 if (shared/float(found)>0.3) {
2299 if (shared/float(found)>0.9 ){
2300 //delete arr->RemoveAt(i);
2305 Bool_t isok =kFALSE;
2306 if ( (pt->GetNShared()/pt->GetNumberOfClusters()<0.5) &&pt->GetNumberOfClusters()>60)
2308 if ((TMath::Abs(1/pt->GetC())<100.) && (pt->GetNShared()/pt->GetNumberOfClusters()<0.7))
2310 if (TMath::Abs(pt->GetZ()/pt->GetX())>1.1)
2312 if ( (TMath::Abs(pt->GetSnp()>0.7) && pt->GetD(0,0)>60.))
2316 for (Int_t i=0; i<160; i++) {
2317 Int_t index=pt->GetClusterIndex2(i);
2318 if (index<0) continue;
2319 AliTPCclusterMI *c= pt->GetClusterPointer(i);
2321 //if (!(c->IsUsed(10))) c->Use();
2328 Double_t maxchi = meanchi+2.*smeanchi;
2330 for (Int_t i=0; i<nseed; i++) {
2331 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2335 //if (!(pt->IsActive())) continue;
2336 if (pt->GetBSigned()) continue;
2337 Double_t chi = pt->GetChi2()/pt->GetNumberOfClusters();
2338 if (chi>maxchi) continue;
2341 Float_t dens = pt->GetNumberOfClusters()/Float_t(pt->GetNFoundable());
2343 //sign only tracks with enoug big density at the beginning
2345 if ((pt->GetDensityFirst(40)<0.75) && pt->GetNumberOfClusters()<meann) continue;
2348 Double_t mindens = TMath::Max(double(mdensity-sdensity*fdensity*bfactor),0.65);
2349 Double_t minn = TMath::Max(Int_t(meann-fnumber*smeann*bfactor),50);
2351 // if (pt->fBConstrain) mindens = TMath::Max(mdensity-sdensity*fdensity*bfactor,0.65);
2352 if ( (pt->GetRemoval()==10) && (pt->GetSnp()>0.8)&&(dens>mindens))
2355 if ((dens>mindens && pt->GetNumberOfClusters()>minn) && chi<maxchi ){
2356 //Int_t noc=pt->GetNumberOfClusters();
2357 pt->SetBSigned(kTRUE);
2358 for (Int_t i=0; i<160; i++) {
2360 Int_t index=pt->GetClusterIndex2(i);
2361 if (index<0) continue;
2362 AliTPCclusterMI *c= pt->GetClusterPointer(i);
2364 // if (!(c->IsUsed(10))) c->Use();
2369 // gLastCheck = nseed;
2377 void AliTPCtrackerMI::StopNotActive(TObjArray * arr, Int_t row0, Float_t th0, Float_t th1, Float_t th2) const
2379 // stop not active tracks
2380 // take th1 as threshold for number of founded to number of foundable on last 10 active rows
2381 // take th2 as threshold for number of founded to number of foundable on last 20 active rows
2382 Int_t nseed = arr->GetEntriesFast();
2384 for (Int_t i=0; i<nseed; i++) {
2385 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2389 if (!(pt->IsActive())) continue;
2390 StopNotActive(pt,row0,th0, th1,th2);
2396 void AliTPCtrackerMI::StopNotActive(AliTPCseed * seed, Int_t row0, Float_t th0, Float_t th1,
2399 // stop not active tracks
2400 // take th1 as threshold for number of founded to number of foundable on last 10 active rows
2401 // take th2 as threshold for number of founded to number of foundable on last 20 active rows
2404 Int_t foundable = 0;
2405 Int_t maxindex = seed->GetLastPoint(); //last foundable row
2406 if (seed->GetNFoundable()*th0 > seed->GetNumberOfClusters()) {
2407 seed->Desactivate(10) ;
2411 for (Int_t i=row0; i<maxindex; i++){
2412 Int_t index = seed->GetClusterIndex2(i);
2413 if (index!=-1) foundable++;
2415 if (foundable<=30) sumgood1++;
2416 if (foundable<=50) {
2423 if (foundable>=30.){
2424 if (sumgood1<(th1*30.)) seed->Desactivate(10);
2427 if (sumgood2<(th2*50.)) seed->Desactivate(10);
2431 Int_t AliTPCtrackerMI::RefitInward(AliESDEvent *event)
2434 // back propagation of ESD tracks
2437 const Int_t kMaxFriendTracks=2000;
2441 //PrepareForProlongation(fSeeds,1);
2442 PropagateForward2(fSeeds);
2443 RemoveUsed2(fSeeds,0.4,0.4,20);
2445 TObjArray arraySeed(fSeeds->GetEntries());
2446 for (Int_t i=0;i<fSeeds->GetEntries();i++) {
2447 arraySeed.AddAt(fSeeds->At(i),i);
2449 SignShared(&arraySeed);
2450 // FindCurling(fSeeds, event,2); // find multi found tracks
2451 FindSplitted(fSeeds, event,2); // find multi found tracks
2454 Int_t nseed = fSeeds->GetEntriesFast();
2455 for (Int_t i=0;i<nseed;i++){
2456 AliTPCseed * seed = (AliTPCseed*) fSeeds->UncheckedAt(i);
2457 if (!seed) continue;
2458 if (seed->GetKinkIndex(0)>0) UpdateKinkQualityD(seed); // update quality informations for kinks
2460 seed->PropagateTo(fParam->GetInnerRadiusLow());
2461 seed->UpdatePoints();
2463 AliESDtrack *esd=event->GetTrack(i);
2464 seed->CookdEdx(0.02,0.6);
2465 CookLabel(seed,0.1); //For comparison only
2466 esd->SetTPCClusterMap(seed->GetClusterMap());
2467 esd->SetTPCSharedMap(seed->GetSharedMap());
2469 if (AliTPCReconstructor::StreamLevel()>1 && seed!=0&&esd!=0) {
2470 TTreeSRedirector &cstream = *fDebugStreamer;
2477 if (seed->GetNumberOfClusters()>15){
2478 esd->UpdateTrackParams(seed,AliESDtrack::kTPCrefit);
2479 esd->SetTPCPoints(seed->GetPoints());
2480 esd->SetTPCPointsF(seed->GetNFoundable());
2481 Int_t ndedx = seed->GetNCDEDX(0)+seed->GetNCDEDX(1)+seed->GetNCDEDX(2)+seed->GetNCDEDX(3);
2482 Float_t sdedx = (seed->GetSDEDX(0)+seed->GetSDEDX(1)+seed->GetSDEDX(2)+seed->GetSDEDX(3))*0.25;
2483 Float_t dedx = seed->GetdEdx();
2484 esd->SetTPCsignal(dedx, sdedx, ndedx);
2486 // add seed to the esd track in Calib level
2488 Bool_t storeFriend = gRandom->Rndm()<(kMaxFriendTracks)/Float_t(nseed);
2489 if (AliTPCReconstructor::StreamLevel()>0 &&storeFriend){
2490 AliTPCseed * seedCopy = new AliTPCseed(*seed, kTRUE);
2491 esd->AddCalibObject(seedCopy);
2496 //printf("problem\n");
2499 //FindKinks(fSeeds,event);
2500 Info("RefitInward","Number of refitted tracks %d",ntracks);
2506 Int_t AliTPCtrackerMI::PropagateBack(AliESDEvent *event)
2509 // back propagation of ESD tracks
2515 PropagateBack(fSeeds);
2516 RemoveUsed2(fSeeds,0.4,0.4,20);
2517 //FindCurling(fSeeds, fEvent,1);
2518 FindSplitted(fSeeds, event,1); // find multi found tracks
2521 Int_t nseed = fSeeds->GetEntriesFast();
2523 for (Int_t i=0;i<nseed;i++){
2524 AliTPCseed * seed = (AliTPCseed*) fSeeds->UncheckedAt(i);
2525 if (!seed) continue;
2526 if (seed->GetKinkIndex(0)<0) UpdateKinkQualityM(seed); // update quality informations for kinks
2527 seed->UpdatePoints();
2528 AliESDtrack *esd=event->GetTrack(i);
2529 seed->CookdEdx(0.02,0.6);
2530 CookLabel(seed,0.1); //For comparison only
2531 if (seed->GetNumberOfClusters()>15){
2532 esd->UpdateTrackParams(seed,AliESDtrack::kTPCout);
2533 esd->SetTPCPoints(seed->GetPoints());
2534 esd->SetTPCPointsF(seed->GetNFoundable());
2535 Int_t ndedx = seed->GetNCDEDX(0)+seed->GetNCDEDX(1)+seed->GetNCDEDX(2)+seed->GetNCDEDX(3);
2536 Float_t sdedx = (seed->GetSDEDX(0)+seed->GetSDEDX(1)+seed->GetSDEDX(2)+seed->GetSDEDX(3))*0.25;
2537 Float_t dedx = seed->GetdEdx();
2538 esd->SetTPCsignal(dedx, sdedx, ndedx);
2540 Int_t eventnumber = event->GetEventNumberInFile();// patch 28 fev 06
2541 // This is most likely NOT the event number you'd like to use. It has nothing to do with the 'real' event number
2542 if (AliTPCReconstructor::StreamLevel()>1) {
2543 (*fDebugStreamer)<<"Cback"<<
2545 "EventNrInFile="<<eventnumber<<
2546 "\n"; // patch 28 fev 06
2550 //FindKinks(fSeeds,event);
2551 Info("PropagateBack","Number of back propagated tracks %d",ntracks);
2558 void AliTPCtrackerMI::DeleteSeeds()
2563 Int_t nseed = fSeeds->GetEntriesFast();
2564 for (Int_t i=0;i<nseed;i++){
2565 AliTPCseed * seed = (AliTPCseed*)fSeeds->At(i);
2566 if (seed) delete fSeeds->RemoveAt(i);
2573 void AliTPCtrackerMI::ReadSeeds(AliESDEvent *event, Int_t direction)
2576 //read seeds from the event
2578 Int_t nentr=event->GetNumberOfTracks();
2580 Info("ReadSeeds", "Number of ESD tracks: %d\n", nentr);
2585 fSeeds = new TObjArray(nentr);
2589 for (Int_t i=0; i<nentr; i++) {
2590 AliESDtrack *esd=event->GetTrack(i);
2591 ULong_t status=esd->GetStatus();
2592 if (!(status&AliESDtrack::kTPCin)) continue;
2593 AliTPCtrack t(*esd);
2594 t.SetNumberOfClusters(0);
2595 // AliTPCseed *seed = new AliTPCseed(t,t.GetAlpha());
2596 AliTPCseed *seed = new AliTPCseed(t/*,t.GetAlpha()*/);
2597 seed->SetUniqueID(esd->GetID());
2598 for (Int_t ikink=0;ikink<3;ikink++) {
2599 Int_t index = esd->GetKinkIndex(ikink);
2600 seed->GetKinkIndexes()[ikink] = index;
2601 if (index==0) continue;
2602 index = TMath::Abs(index);
2603 AliESDkink * kink = fEvent->GetKink(index-1);
2604 if (kink&&esd->GetKinkIndex(ikink)<0){
2605 if ((status & AliESDtrack::kTRDrefit) != 0) kink->SetStatus(1,2);
2606 if ((status & AliESDtrack::kITSout) != 0) kink->SetStatus(1,0);
2608 if (kink&&esd->GetKinkIndex(ikink)>0){
2609 if ((status & AliESDtrack::kTRDrefit) != 0) kink->SetStatus(1,6);
2610 if ((status & AliESDtrack::kITSout) != 0) kink->SetStatus(1,4);
2614 if (((status&AliESDtrack::kITSout)==0)&&(direction==1)) seed->ResetCovariance(10.);
2615 if ( direction ==2 &&(status & AliESDtrack::kTRDrefit) == 0 ) seed->ResetCovariance(10.);
2616 if ( direction ==2 && ((status & AliESDtrack::kTPCout) == 0) ) {
2621 if ( direction ==2 &&(status & AliESDtrack::kTRDrefit) > 0 ) {
2622 Double_t par0[5],par1[5],alpha,x;
2623 esd->GetInnerExternalParameters(alpha,x,par0);
2624 esd->GetExternalParameters(x,par1);
2625 Double_t delta1 = TMath::Abs(par0[4]-par1[4])/(0.000000001+TMath::Abs(par0[4]+par1[4]));
2626 Double_t delta2 = TMath::Abs(par0[3]-par1[3]);
2628 if (esd->GetTRDncls()>0) trdchi2 = esd->GetTRDchi2()/esd->GetTRDncls();
2629 //reset covariance if suspicious
2630 if ( (delta1>0.1) || (delta2>0.006) ||trdchi2>7.)
2631 seed->ResetCovariance(10.);
2636 // rotate to the local coordinate system
2638 fSectors=fInnerSec; fN=fkNIS;
2639 Double_t alpha=seed->GetAlpha() - fSectors->GetAlphaShift();
2640 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
2641 if (alpha < 0. ) alpha += 2.*TMath::Pi();
2642 Int_t ns=Int_t(alpha/fSectors->GetAlpha())%fN;
2643 alpha =ns*fSectors->GetAlpha() + fSectors->GetAlphaShift();
2644 if (alpha<-TMath::Pi()) alpha += 2*TMath::Pi();
2645 if (alpha>=TMath::Pi()) alpha -= 2*TMath::Pi();
2646 alpha-=seed->GetAlpha();
2647 if (!seed->Rotate(alpha)) {
2653 if (esd->GetKinkIndex(0)<=0){
2654 for (Int_t irow=0;irow<160;irow++){
2655 Int_t index = seed->GetClusterIndex2(irow);
2658 AliTPCclusterMI * cl = GetClusterMI(index);
2659 seed->SetClusterPointer(irow,cl);
2661 if ((index & 0x8000)==0){
2662 cl->Use(10); // accepted cluster
2664 cl->Use(6); // close cluster not accepted
2667 Info("ReadSeeds","Not found cluster");
2672 fSeeds->AddAt(seed,i);
2678 //_____________________________________________________________________________
2679 void AliTPCtrackerMI::MakeSeeds3(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, Float_t cuts[4],
2680 Float_t deltay, Int_t ddsec) {
2681 //-----------------------------------------------------------------
2682 // This function creates track seeds.
2683 // SEEDING WITH VERTEX CONSTRAIN
2684 //-----------------------------------------------------------------
2685 // cuts[0] - fP4 cut
2686 // cuts[1] - tan(phi) cut
2687 // cuts[2] - zvertex cut
2688 // cuts[3] - fP3 cut
2696 Double_t x[5], c[15];
2697 // Int_t di = i1-i2;
2699 AliTPCseed * seed = new AliTPCseed();
2700 Double_t alpha=fSectors->GetAlpha(), shift=fSectors->GetAlphaShift();
2701 Double_t cs=cos(alpha), sn=sin(alpha);
2703 // Double_t x1 =fOuterSec->GetX(i1);
2704 //Double_t xx2=fOuterSec->GetX(i2);
2706 Double_t x1 =GetXrow(i1);
2707 Double_t xx2=GetXrow(i2);
2709 Double_t x3=GetX(), y3=GetY(), z3=GetZ();
2711 Int_t imiddle = (i2+i1)/2; //middle pad row index
2712 Double_t xm = GetXrow(imiddle); // radius of middle pad-row
2713 const AliTPCRow& krm=GetRow(sec,imiddle); //middle pad -row
2717 const AliTPCRow& kr1=GetRow(ns,i1);
2718 Double_t ymax = GetMaxY(i1)-kr1.GetDeadZone()-1.5;
2719 Double_t ymaxm = GetMaxY(imiddle)-kr1.GetDeadZone()-1.5;
2722 // change cut on curvature if it can't reach this layer
2723 // maximal curvature set to reach it
2724 Double_t dvertexmax = TMath::Sqrt((x1-x3)*(x1-x3)+(ymax+5-y3)*(ymax+5-y3));
2725 if (dvertexmax*0.5*cuts[0]>0.85){
2726 cuts[0] = 0.85/(dvertexmax*0.5+1.);
2728 Double_t r2min = 1/(cuts[0]*cuts[0]); //minimal square of radius given by cut
2731 if (deltay>0) ddsec = 0;
2732 // loop over clusters
2733 for (Int_t is=0; is < kr1; is++) {
2735 if (kr1[is]->IsUsed(10)) continue;
2736 Double_t y1=kr1[is]->GetY(), z1=kr1[is]->GetZ();
2737 //if (TMath::Abs(y1)>ymax) continue;
2739 if (deltay>0 && TMath::Abs(ymax-TMath::Abs(y1))> deltay ) continue; // seed only at the edge
2741 // find possible directions
2742 Float_t anglez = (z1-z3)/(x1-x3);
2743 Float_t extraz = z1 - anglez*(x1-xx2); // extrapolated z
2746 //find rotation angles relative to line given by vertex and point 1
2747 Double_t dvertex2 = (x1-x3)*(x1-x3)+(y1-y3)*(y1-y3);
2748 Double_t dvertex = TMath::Sqrt(dvertex2);
2749 Double_t angle13 = TMath::ATan((y1-y3)/(x1-x3));
2750 Double_t cs13 = cos(-angle13), sn13 = sin(-angle13);
2753 // loop over 2 sectors
2759 Double_t dddz1=0; // direction of delta inclination in z axis
2766 for (Int_t dsec = dsec1; dsec<=dsec2;dsec++){
2767 Int_t sec2 = sec + dsec;
2769 // AliTPCRow& kr2 = fOuterSec[(sec2+fkNOS)%fkNOS][i2];
2770 //AliTPCRow& kr2m = fOuterSec[(sec2+fkNOS)%fkNOS][imiddle];
2771 AliTPCRow& kr2 = GetRow((sec2+fkNOS)%fkNOS,i2);
2772 AliTPCRow& kr2m = GetRow((sec2+fkNOS)%fkNOS,imiddle);
2773 Int_t index1 = TMath::Max(kr2.Find(extraz-0.6-dddz1*TMath::Abs(z1)*0.05)-1,0);
2774 Int_t index2 = TMath::Min(kr2.Find(extraz+0.6+dddz2*TMath::Abs(z1)*0.05)+1,kr2);
2776 // rotation angles to p1-p3
2777 Double_t cs13r = cos(-angle13+dsec*alpha)/dvertex, sn13r = sin(-angle13+dsec*alpha)/dvertex;
2778 Double_t x2, y2, z2;
2780 // Double_t dymax = maxangle*TMath::Abs(x1-xx2);
2783 Double_t dxx0 = (xx2-x3)*cs13r;
2784 Double_t dyy0 = (xx2-x3)*sn13r;
2785 for (Int_t js=index1; js < index2; js++) {
2786 const AliTPCclusterMI *kcl = kr2[js];
2787 if (kcl->IsUsed(10)) continue;
2789 //calcutate parameters
2791 Double_t yy0 = dyy0 +(kcl->GetY()-y3)*cs13r;
2793 if (TMath::Abs(yy0)<0.000001) continue;
2794 Double_t xx0 = dxx0 -(kcl->GetY()-y3)*sn13r;
2795 Double_t y0 = 0.5*(xx0*xx0+yy0*yy0-xx0)/yy0;
2796 Double_t r02 = (0.25+y0*y0)*dvertex2;
2797 //curvature (radius) cut
2798 if (r02<r2min) continue;
2802 Double_t c0 = 1/TMath::Sqrt(r02);
2806 //Double_t dfi0 = 2.*TMath::ASin(dvertex*c0*0.5);
2807 //Double_t dfi1 = 2.*TMath::ASin(TMath::Sqrt(yy0*yy0+(1-xx0)*(1-xx0))*dvertex*c0*0.5);
2808 Double_t dfi0 = 2.*AliTPCFastMath::FastAsin(dvertex*c0*0.5);
2809 Double_t dfi1 = 2.*AliTPCFastMath::FastAsin(TMath::Sqrt(yy0*yy0+(1-xx0)*(1-xx0))*dvertex*c0*0.5);
2812 Double_t z0 = kcl->GetZ();
2813 Double_t zzzz2 = z1-(z1-z3)*dfi1/dfi0;
2814 if (TMath::Abs(zzzz2-z0)>0.5) continue;
2817 Double_t dip = (z1-z0)*c0/dfi1;
2818 Double_t x0 = (0.5*cs13+y0*sn13)*dvertex*c0;
2829 x2= xx2*cs-y2*sn*dsec;
2830 y2=+xx2*sn*dsec+y2*cs;
2840 // do we have cluster at the middle ?
2842 GetProlongation(x1,xm,x,ym,zm);
2844 AliTPCclusterMI * cm=0;
2845 if (TMath::Abs(ym)-ymaxm<0){
2846 cm = krm.FindNearest2(ym,zm,1.0,0.6,dummy);
2847 if ((!cm) || (cm->IsUsed(10))) {
2852 // rotate y1 to system 0
2853 // get state vector in rotated system
2854 Double_t yr1 = (-0.5*sn13+y0*cs13)*dvertex*c0;
2855 Double_t xr2 = x0*cs+yr1*sn*dsec;
2856 Double_t xr[5]={kcl->GetY(),kcl->GetZ(), xr2, dip, c0};
2858 GetProlongation(xx2,xm,xr,ym,zm);
2859 if (TMath::Abs(ym)-ymaxm<0){
2860 cm = kr2m.FindNearest2(ym,zm,1.0,0.6,dummy);
2861 if ((!cm) || (cm->IsUsed(10))) {
2871 dym = ym - cm->GetY();
2872 dzm = zm - cm->GetZ();
2879 Double_t sy1=kr1[is]->GetSigmaY2()*2., sz1=kr1[is]->GetSigmaZ2()*2.;
2880 Double_t sy2=kcl->GetSigmaY2()*2., sz2=kcl->GetSigmaZ2()*2.;
2881 //Double_t sy3=400*3./12., sy=0.1, sz=0.1;
2882 Double_t sy3=25000*x[4]*x[4]+0.1, sy=0.1, sz=0.1;
2883 //Double_t sy3=25000*x[4]*x[4]*60+0.5, sy=0.1, sz=0.1;
2885 Double_t f40=(F1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
2886 Double_t f42=(F1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
2887 Double_t f43=(F1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
2888 Double_t f20=(F2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
2889 Double_t f22=(F2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
2890 Double_t f23=(F2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
2892 Double_t f30=(F3(x1,y1+sy,x2,y2,z1,z2)-x[3])/sy;
2893 Double_t f31=(F3(x1,y1,x2,y2,z1+sz,z2)-x[3])/sz;
2894 Double_t f32=(F3(x1,y1,x2,y2+sy,z1,z2)-x[3])/sy;
2895 Double_t f34=(F3(x1,y1,x2,y2,z1,z2+sz)-x[3])/sz;
2899 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
2900 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
2901 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
2902 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
2903 c[13]=f30*sy1*f40+f32*sy2*f42;
2904 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
2906 // if (!BuildSeed(kr1[is],kcl,0,x1,x2,x3,x,c)) continue;
2908 UInt_t index=kr1.GetIndex(is);
2909 seed->~AliTPCseed(); // this does not set the pointer to 0...
2910 AliTPCseed *track=new(seed) AliTPCseed(x1, ns*alpha+shift, x, c, index);
2912 track->SetIsSeeding(kTRUE);
2913 track->SetSeed1(i1);
2914 track->SetSeed2(i2);
2915 track->SetSeedType(3);
2919 FollowProlongation(*track, (i1+i2)/2,1);
2920 Int_t foundable,found,shared;
2921 track->GetClusterStatistic((i1+i2)/2,i1, found, foundable, shared, kTRUE);
2922 if ((found<0.55*foundable) || shared>0.5*found || (track->GetSigmaY2()+track->GetSigmaZ2())>0.5){
2924 seed->~AliTPCseed();
2930 FollowProlongation(*track, i2,1);
2934 track->SetBConstrain(1);
2935 // track->fLastPoint = i1+fInnerSec->GetNRows(); // first cluster in track position
2936 track->SetLastPoint(i1); // first cluster in track position
2937 track->SetFirstPoint(track->GetLastPoint());
2939 if (track->GetNumberOfClusters()<(i1-i2)*0.5 ||
2940 track->GetNumberOfClusters() < track->GetNFoundable()*0.6 ||
2941 track->GetNShared()>0.4*track->GetNumberOfClusters() ) {
2943 seed->~AliTPCseed();
2947 // Z VERTEX CONDITION
2948 Double_t zv, bz=GetBz();
2949 if ( !track->GetZAt(0.,bz,zv) ) continue;
2950 if (TMath::Abs(zv-z3)>cuts[2]) {
2951 FollowProlongation(*track, TMath::Max(i2-20,0));
2952 if ( !track->GetZAt(0.,bz,zv) ) continue;
2953 if (TMath::Abs(zv-z3)>cuts[2]){
2954 FollowProlongation(*track, TMath::Max(i2-40,0));
2955 if ( !track->GetZAt(0.,bz,zv) ) continue;
2956 if (TMath::Abs(zv-z3)>cuts[2] &&(track->GetNumberOfClusters() > track->GetNFoundable()*0.7)){
2957 // make seed without constrain
2958 AliTPCseed * track2 = MakeSeed(track,0.2,0.5,1.);
2959 FollowProlongation(*track2, i2,1);
2960 track2->SetBConstrain(kFALSE);
2961 track2->SetSeedType(1);
2962 arr->AddLast(track2);
2964 seed->~AliTPCseed();
2969 seed->~AliTPCseed();
2976 track->SetSeedType(0);
2977 arr->AddLast(track);
2978 seed = new AliTPCseed;
2980 // don't consider other combinations
2981 if (track->GetNumberOfClusters() > track->GetNFoundable()*0.8)
2987 Info("MakeSeeds3","\nSeeding statistic:\t%d\t%d\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin,nout1,nout2);
2993 void AliTPCtrackerMI::MakeSeeds5(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, Float_t cuts[4],
2998 //-----------------------------------------------------------------
2999 // This function creates track seeds.
3000 //-----------------------------------------------------------------
3001 // cuts[0] - fP4 cut
3002 // cuts[1] - tan(phi) cut
3003 // cuts[2] - zvertex cut
3004 // cuts[3] - fP3 cut
3014 Double_t x[5], c[15];
3016 // make temporary seed
3017 AliTPCseed * seed = new AliTPCseed;
3018 Double_t alpha=fOuterSec->GetAlpha(), shift=fOuterSec->GetAlphaShift();
3019 // Double_t cs=cos(alpha), sn=sin(alpha);
3024 Double_t x1 = GetXrow(i1-1);
3025 const AliTPCRow& kr1=GetRow(sec,i1-1);
3026 Double_t y1max = GetMaxY(i1-1)-kr1.GetDeadZone()-1.5;
3028 Double_t x1p = GetXrow(i1);
3029 const AliTPCRow& kr1p=GetRow(sec,i1);
3031 Double_t x1m = GetXrow(i1-2);
3032 const AliTPCRow& kr1m=GetRow(sec,i1-2);
3035 //last 3 padrow for seeding
3036 AliTPCRow& kr3 = GetRow((sec+fkNOS)%fkNOS,i1-7);
3037 Double_t x3 = GetXrow(i1-7);
3038 // Double_t y3max= GetMaxY(i1-7)-kr3.fDeadZone-1.5;
3040 AliTPCRow& kr3p = GetRow((sec+fkNOS)%fkNOS,i1-6);
3041 Double_t x3p = GetXrow(i1-6);
3043 AliTPCRow& kr3m = GetRow((sec+fkNOS)%fkNOS,i1-8);
3044 Double_t x3m = GetXrow(i1-8);
3049 Int_t im = i1-4; //middle pad row index
3050 Double_t xm = GetXrow(im); // radius of middle pad-row
3051 const AliTPCRow& krm=GetRow(sec,im); //middle pad -row
3052 // Double_t ymmax = GetMaxY(im)-kr1.fDeadZone-1.5;
3055 Double_t deltax = x1-x3;
3056 Double_t dymax = deltax*cuts[1];
3057 Double_t dzmax = deltax*cuts[3];
3059 // loop over clusters
3060 for (Int_t is=0; is < kr1; is++) {
3062 if (kr1[is]->IsUsed(10)) continue;
3063 Double_t y1=kr1[is]->GetY(), z1=kr1[is]->GetZ();
3065 if (deltay>0 && TMath::Abs(y1max-TMath::Abs(y1))> deltay ) continue; // seed only at the edge
3067 Int_t index1 = TMath::Max(kr3.Find(z1-dzmax)-1,0);
3068 Int_t index2 = TMath::Min(kr3.Find(z1+dzmax)+1,kr3);
3074 for (Int_t js=index1; js < index2; js++) {
3075 const AliTPCclusterMI *kcl = kr3[js];
3076 if (kcl->IsUsed(10)) continue;
3078 // apply angular cuts
3079 if (TMath::Abs(y1-y3)>dymax) continue;
3082 if (TMath::Abs(z1-z3)>dzmax) continue;
3084 Double_t angley = (y1-y3)/(x1-x3);
3085 Double_t anglez = (z1-z3)/(x1-x3);
3087 Double_t erry = TMath::Abs(angley)*(x1-x1m)*0.5+0.5;
3088 Double_t errz = TMath::Abs(anglez)*(x1-x1m)*0.5+0.5;
3090 Double_t yyym = angley*(xm-x1)+y1;
3091 Double_t zzzm = anglez*(xm-x1)+z1;
3093 const AliTPCclusterMI *kcm = krm.FindNearest2(yyym,zzzm,erry,errz,index);
3095 if (kcm->IsUsed(10)) continue;
3097 erry = TMath::Abs(angley)*(x1-x1m)*0.4+0.5;
3098 errz = TMath::Abs(anglez)*(x1-x1m)*0.4+0.5;
3105 // look around first
3106 const AliTPCclusterMI *kc1m = kr1m.FindNearest2(angley*(x1m-x1)+y1,
3112 if (kc1m->IsUsed(10)) used++;
3114 const AliTPCclusterMI *kc1p = kr1p.FindNearest2(angley*(x1p-x1)+y1,
3120 if (kc1p->IsUsed(10)) used++;
3122 if (used>1) continue;
3123 if (found<1) continue;
3127 const AliTPCclusterMI *kc3m = kr3m.FindNearest2(angley*(x3m-x3)+y3,
3133 if (kc3m->IsUsed(10)) used++;
3137 const AliTPCclusterMI *kc3p = kr3p.FindNearest2(angley*(x3p-x3)+y3,
3143 if (kc3p->IsUsed(10)) used++;
3147 if (used>1) continue;
3148 if (found<3) continue;
3158 x[4]=F1(x1,y1,x2,y2,x3,y3);
3159 //if (TMath::Abs(x[4]) >= cuts[0]) continue;
3162 x[2]=F2(x1,y1,x2,y2,x3,y3);
3165 x[3]=F3n(x1,y1,x2,y2,z1,z2,x[4]);
3166 //if (TMath::Abs(x[3]) > cuts[3]) continue;
3170 Double_t sy1=0.1, sz1=0.1;
3171 Double_t sy2=0.1, sz2=0.1;
3172 Double_t sy3=0.1, sy=0.1, sz=0.1;
3174 Double_t f40=(F1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
3175 Double_t f42=(F1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
3176 Double_t f43=(F1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
3177 Double_t f20=(F2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
3178 Double_t f22=(F2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
3179 Double_t f23=(F2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
3181 Double_t f30=(F3(x1,y1+sy,x2,y2,z1,z2)-x[3])/sy;
3182 Double_t f31=(F3(x1,y1,x2,y2,z1+sz,z2)-x[3])/sz;
3183 Double_t f32=(F3(x1,y1,x2,y2+sy,z1,z2)-x[3])/sy;
3184 Double_t f34=(F3(x1,y1,x2,y2,z1,z2+sz)-x[3])/sz;
3188 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3189 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3190 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3191 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3192 c[13]=f30*sy1*f40+f32*sy2*f42;
3193 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3195 // if (!BuildSeed(kr1[is],kcl,0,x1,x2,x3,x,c)) continue;
3197 UInt_t index=kr1.GetIndex(is);
3198 seed->~AliTPCseed();
3199 AliTPCseed *track=new(seed) AliTPCseed(x1, sec*alpha+shift, x, c, index);
3201 track->SetIsSeeding(kTRUE);
3204 FollowProlongation(*track, i1-7,1);
3205 if (track->GetNumberOfClusters() < track->GetNFoundable()*0.75 ||
3206 track->GetNShared()>0.6*track->GetNumberOfClusters() || ( track->GetSigmaY2()+ track->GetSigmaZ2())>0.6){
3208 seed->~AliTPCseed();
3214 FollowProlongation(*track, i2,1);
3215 track->SetBConstrain(0);
3216 track->SetLastPoint(i1+fInnerSec->GetNRows()); // first cluster in track position
3217 track->SetFirstPoint(track->GetLastPoint());
3219 if (track->GetNumberOfClusters()<(i1-i2)*0.5 ||
3220 track->GetNumberOfClusters()<track->GetNFoundable()*0.7 ||
3221 track->GetNShared()>2. || track->GetChi2()/track->GetNumberOfClusters()>6 || ( track->GetSigmaY2()+ track->GetSigmaZ2())>0.5 ) {
3223 seed->~AliTPCseed();
3228 FollowProlongation(*track, TMath::Max(i2-10,0),1);
3229 AliTPCseed * track2 = MakeSeed(track,0.2,0.5,0.9);
3230 FollowProlongation(*track2, i2,1);
3231 track2->SetBConstrain(kFALSE);
3232 track2->SetSeedType(4);
3233 arr->AddLast(track2);
3235 seed->~AliTPCseed();
3239 //arr->AddLast(track);
3240 //seed = new AliTPCseed;
3246 Info("MakeSeeds5","\nSeeding statiistic:\t%d\t%d\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin,nout1,nout2,nout3);
3252 //_____________________________________________________________________________
3253 void AliTPCtrackerMI::MakeSeeds2(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, Float_t */*cuts[4]*/,
3254 Float_t deltay, Bool_t /*bconstrain*/) {
3255 //-----------------------------------------------------------------
3256 // This function creates track seeds - without vertex constraint
3257 //-----------------------------------------------------------------
3258 // cuts[0] - fP4 cut - not applied
3259 // cuts[1] - tan(phi) cut
3260 // cuts[2] - zvertex cut - not applied
3261 // cuts[3] - fP3 cut
3271 Double_t alpha=fOuterSec->GetAlpha(), shift=fOuterSec->GetAlphaShift();
3272 // Double_t cs=cos(alpha), sn=sin(alpha);
3273 Int_t row0 = (i1+i2)/2;
3274 Int_t drow = (i1-i2)/2;
3275 const AliTPCRow& kr0=fSectors[sec][row0];
3278 AliTPCpolyTrack polytrack;
3279 Int_t nclusters=fSectors[sec][row0];
3280 AliTPCseed * seed = new AliTPCseed;
3285 for (Int_t is=0; is < nclusters; is++) { //LOOP over clusters
3287 Int_t nfoundable =0;
3288 for (Int_t iter =1; iter<2; iter++){ //iterations
3289 const AliTPCRow& krm=fSectors[sec][row0-iter];
3290 const AliTPCRow& krp=fSectors[sec][row0+iter];
3291 const AliTPCclusterMI * cl= kr0[is];
3293 if (cl->IsUsed(10)) {
3299 Double_t x = kr0.GetX();
3300 // Initialization of the polytrack
3305 Double_t y0= cl->GetY();
3306 Double_t z0= cl->GetZ();
3310 Double_t ymax = fSectors->GetMaxY(row0)-kr0.GetDeadZone()-1.5;
3311 if (deltay>0 && TMath::Abs(ymax-TMath::Abs(y0))> deltay ) continue; // seed only at the edge
3313 erry = (0.5)*cl->GetSigmaY2()/TMath::Sqrt(cl->GetQ())*6;
3314 errz = (0.5)*cl->GetSigmaZ2()/TMath::Sqrt(cl->GetQ())*6;
3315 polytrack.AddPoint(x,y0,z0,erry, errz);
3318 if (cl->IsUsed(10)) sumused++;
3321 Float_t roady = (5*TMath::Sqrt(cl->GetSigmaY2()+0.2)+1.)*iter;
3322 Float_t roadz = (5*TMath::Sqrt(cl->GetSigmaZ2()+0.2)+1.)*iter;
3325 AliTPCclusterMI * cl1 = krm.FindNearest(y0,z0,roady,roadz);
3326 if (cl1 && TMath::Abs(ymax-TMath::Abs(y0))) {
3327 erry = (0.5)*cl1->GetSigmaY2()/TMath::Sqrt(cl1->GetQ())*3;
3328 errz = (0.5)*cl1->GetSigmaZ2()/TMath::Sqrt(cl1->GetQ())*3;
3329 if (cl1->IsUsed(10)) sumused++;
3330 polytrack.AddPoint(x,cl1->GetY(),cl1->GetZ(),erry,errz);
3334 AliTPCclusterMI * cl2 = krp.FindNearest(y0,z0,roady,roadz);
3336 erry = (0.5)*cl2->GetSigmaY2()/TMath::Sqrt(cl2->GetQ())*3;
3337 errz = (0.5)*cl2->GetSigmaZ2()/TMath::Sqrt(cl2->GetQ())*3;
3338 if (cl2->IsUsed(10)) sumused++;
3339 polytrack.AddPoint(x,cl2->GetY(),cl2->GetZ(),erry,errz);
3342 if (sumused>0) continue;
3344 polytrack.UpdateParameters();
3350 nfoundable = polytrack.GetN();
3351 nfound = nfoundable;
3353 for (Int_t ddrow = iter+1; ddrow<drow;ddrow++){
3354 Float_t maxdist = 0.8*(1.+3./(ddrow));
3355 for (Int_t delta = -1;delta<=1;delta+=2){
3356 Int_t row = row0+ddrow*delta;
3357 kr = &(fSectors[sec][row]);
3358 Double_t xn = kr->GetX();
3359 Double_t ymax = fSectors->GetMaxY(row)-kr->GetDeadZone()-1.5;
3360 polytrack.GetFitPoint(xn,yn,zn);
3361 if (TMath::Abs(yn)>ymax) continue;
3363 AliTPCclusterMI * cln = kr->FindNearest(yn,zn,roady,roadz);
3365 Float_t dist = TMath::Sqrt( (yn-cln->GetY())*(yn-cln->GetY())+(zn-cln->GetZ())*(zn-cln->GetZ()));
3368 erry = (dist+0.3)*cln->GetSigmaY2()/TMath::Sqrt(cln->GetQ())*(1.+1./(ddrow));
3369 errz = (dist+0.3)*cln->GetSigmaZ2()/TMath::Sqrt(cln->GetQ())*(1.+1./(ddrow));
3370 if (cln->IsUsed(10)) {
3371 // printf("used\n");
3379 polytrack.AddPoint(xn,cln->GetY(),cln->GetZ(),erry, errz);
3384 if ( (sumused>3) || (sumused>0.5*nfound) || (nfound<0.6*nfoundable)) break;
3385 polytrack.UpdateParameters();
3388 if ( (sumused>3) || (sumused>0.5*nfound)) {
3389 //printf("sumused %d\n",sumused);
3394 polytrack.GetFitDerivation(kr0.GetX(),dy,dz);
3395 AliTPCpolyTrack track2;
3397 polytrack.Refit(track2,0.5+TMath::Abs(dy)*0.3,0.4+TMath::Abs(dz)*0.3);
3398 if (track2.GetN()<0.5*nfoundable) continue;
3401 if ((nfound>0.6*nfoundable) &&( nfoundable>0.4*(i1-i2))) {
3403 // test seed with and without constrain
3404 for (Int_t constrain=0; constrain<=0;constrain++){
3405 // add polytrack candidate
3407 Double_t x[5], c[15];
3408 Double_t x1,x2,x3,y1,y2,y3,z1,z2,z3;
3409 track2.GetBoundaries(x3,x1);
3411 track2.GetFitPoint(x1,y1,z1);
3412 track2.GetFitPoint(x2,y2,z2);
3413 track2.GetFitPoint(x3,y3,z3);
3415 //is track pointing to the vertex ?
3418 polytrack.GetFitPoint(x0,y0,z0);
3431 x[4]=F1(x1,y1,x2,y2,x3,y3);
3433 // if (TMath::Abs(x[4]) >= cuts[0]) continue; //
3434 x[2]=F2(x1,y1,x2,y2,x3,y3);
3436 //if (TMath::Abs(x[4]*x1-x[2]) >= cuts[1]) continue;
3437 //x[3]=F3(x1,y1,x2,y2,z1,z2);
3438 x[3]=F3n(x1,y1,x3,y3,z1,z3,x[4]);
3439 //if (TMath::Abs(x[3]) > cuts[3]) continue;
3442 Double_t sy =0.1, sz =0.1;
3443 Double_t sy1=0.02, sz1=0.02;
3444 Double_t sy2=0.02, sz2=0.02;
3448 sy3=25000*x[4]*x[4]+0.1, sy=0.1, sz=0.1;
3451 Double_t f40=(F1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
3452 Double_t f42=(F1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
3453 Double_t f43=(F1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
3454 Double_t f20=(F2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
3455 Double_t f22=(F2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
3456 Double_t f23=(F2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
3458 Double_t f30=(F3(x1,y1+sy,x3,y3,z1,z3)-x[3])/sy;
3459 Double_t f31=(F3(x1,y1,x3,y3,z1+sz,z3)-x[3])/sz;
3460 Double_t f32=(F3(x1,y1,x3,y3+sy,z1,z3)-x[3])/sy;
3461 Double_t f34=(F3(x1,y1,x3,y3,z1,z3+sz)-x[3])/sz;
3466 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3467 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3468 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3469 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3470 c[13]=f30*sy1*f40+f32*sy2*f42;
3471 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3473 //Int_t row1 = fSectors->GetRowNumber(x1);
3474 Int_t row1 = GetRowNumber(x1);
3478 seed->~AliTPCseed();
3479 AliTPCseed *track=new(seed) AliTPCseed(x1,sec*alpha+shift,x,c,index);
3480 track->SetIsSeeding(kTRUE);
3481 Int_t rc=FollowProlongation(*track, i2);
3482 if (constrain) track->SetBConstrain(1);
3484 track->SetBConstrain(0);
3485 track->SetLastPoint(row1+fInnerSec->GetNRows()); // first cluster in track position
3486 track->SetFirstPoint(track->GetLastPoint());
3488 if (rc==0 || track->GetNumberOfClusters()<(i1-i2)*0.5 ||
3489 track->GetNumberOfClusters() < track->GetNFoundable()*0.6 ||
3490 track->GetNShared()>0.4*track->GetNumberOfClusters()) {
3493 seed->~AliTPCseed();
3496 arr->AddLast(track);
3497 seed = new AliTPCseed;
3501 } // if accepted seed
3504 Info("MakeSeeds2","\nSeeding statiistic:\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin3);
3510 AliTPCseed *AliTPCtrackerMI::MakeSeed(AliTPCseed *track, Float_t r0, Float_t r1, Float_t r2)
3514 //reseed using track points
3515 Int_t p0 = int(r0*track->GetNumberOfClusters()); // point 0
3516 Int_t p1 = int(r1*track->GetNumberOfClusters());
3517 Int_t p2 = int(r2*track->GetNumberOfClusters()); // last point
3519 Double_t x0[3],x1[3],x2[3];
3520 for (Int_t i=0;i<3;i++){
3526 // find track position at given ratio of the length
3527 Int_t sec0=0, sec1=0, sec2=0;
3530 for (Int_t i=0;i<160;i++){
3531 if (track->GetClusterPointer(i)){
3533 AliTPCTrackerPoint *trpoint =track->GetTrackPoint(i);
3534 if ( (index<p0) || x0[0]<0 ){
3535 if (trpoint->GetX()>1){
3536 clindex = track->GetClusterIndex2(i);
3538 x0[0] = trpoint->GetX();
3539 x0[1] = trpoint->GetY();
3540 x0[2] = trpoint->GetZ();
3541 sec0 = ((clindex&0xff000000)>>24)%18;
3546 if ( (index<p1) &&(trpoint->GetX()>1)){
3547 clindex = track->GetClusterIndex2(i);
3549 x1[0] = trpoint->GetX();
3550 x1[1] = trpoint->GetY();
3551 x1[2] = trpoint->GetZ();
3552 sec1 = ((clindex&0xff000000)>>24)%18;
3555 if ( (index<p2) &&(trpoint->GetX()>1)){
3556 clindex = track->GetClusterIndex2(i);
3558 x2[0] = trpoint->GetX();
3559 x2[1] = trpoint->GetY();
3560 x2[2] = trpoint->GetZ();
3561 sec2 = ((clindex&0xff000000)>>24)%18;
3568 Double_t alpha, cs,sn, xx2,yy2;
3570 alpha = (sec1-sec2)*fSectors->GetAlpha();
3571 cs = TMath::Cos(alpha);
3572 sn = TMath::Sin(alpha);
3573 xx2= x1[0]*cs-x1[1]*sn;
3574 yy2= x1[0]*sn+x1[1]*cs;
3578 alpha = (sec0-sec2)*fSectors->GetAlpha();
3579 cs = TMath::Cos(alpha);
3580 sn = TMath::Sin(alpha);
3581 xx2= x0[0]*cs-x0[1]*sn;
3582 yy2= x0[0]*sn+x0[1]*cs;
3588 Double_t x[5],c[15];
3592 x[4]=F1(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]);
3593 // if (x[4]>1) return 0;
3594 x[2]=F2(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]);
3595 x[3]=F3n(x2[0],x2[1],x0[0],x0[1],x2[2],x0[2],x[4]);
3596 //if (TMath::Abs(x[3]) > 2.2) return 0;
3597 //if (TMath::Abs(x[2]) > 1.99) return 0;
3599 Double_t sy =0.1, sz =0.1;
3601 Double_t sy1=0.02+track->GetSigmaY2(), sz1=0.02+track->GetSigmaZ2();
3602 Double_t sy2=0.01+track->GetSigmaY2(), sz2=0.01+track->GetSigmaZ2();
3603 Double_t sy3=0.01+track->GetSigmaY2();
3605 Double_t f40=(F1(x2[0],x2[1]+sy,x1[0],x1[1],x0[0],x0[1])-x[4])/sy;
3606 Double_t f42=(F1(x2[0],x2[1],x1[0],x1[1]+sy,x0[0],x0[1])-x[4])/sy;
3607 Double_t f43=(F1(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]+sy)-x[4])/sy;
3608 Double_t f20=(F2(x2[0],x2[1]+sy,x1[0],x1[1],x0[0],x0[1])-x[2])/sy;
3609 Double_t f22=(F2(x2[0],x2[1],x1[0],x1[1]+sy,x0[0],x0[1])-x[2])/sy;
3610 Double_t f23=(F2(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]+sy)-x[2])/sy;
3612 Double_t f30=(F3(x2[0],x2[1]+sy,x0[0],x0[1],x2[2],x0[2])-x[3])/sy;
3613 Double_t f31=(F3(x2[0],x2[1],x0[0],x0[1],x2[2]+sz,x0[2])-x[3])/sz;
3614 Double_t f32=(F3(x2[0],x2[1],x0[0],x0[1]+sy,x2[2],x0[2])-x[3])/sy;
3615 Double_t f34=(F3(x2[0],x2[1],x0[0],x0[1],x2[2],x0[2]+sz)-x[3])/sz;
3620 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3621 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3622 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3623 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3624 c[13]=f30*sy1*f40+f32*sy2*f42;
3625 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3627 // Int_t row1 = fSectors->GetRowNumber(x2[0]);
3628 AliTPCseed *seed=new AliTPCseed(x2[0], sec2*fSectors->GetAlpha()+fSectors->GetAlphaShift(), x, c, 0);
3629 // Double_t y0,z0,y1,z1, y2,z2;
3630 //seed->GetProlongation(x0[0],y0,z0);
3631 // seed->GetProlongation(x1[0],y1,z1);
3632 //seed->GetProlongation(x2[0],y2,z2);
3634 seed->SetLastPoint(pp2);
3635 seed->SetFirstPoint(pp2);
3642 AliTPCseed *AliTPCtrackerMI::ReSeed(AliTPCseed *track, Float_t r0, Float_t r1, Float_t r2)
3646 //reseed using founded clusters
3648 // Find the number of clusters
3649 Int_t nclusters = 0;
3650 for (Int_t irow=0;irow<160;irow++){
3651 if (track->GetClusterIndex(irow)>0) nclusters++;
3655 ipos[0] = TMath::Max(int(r0*nclusters),0); // point 0 cluster
3656 ipos[1] = TMath::Min(int(r1*nclusters),nclusters-1); //
3657 ipos[2] = TMath::Min(int(r2*nclusters),nclusters-1); // last point
3661 Int_t row[3],sec[3]={0,0,0};
3663 // find track row position at given ratio of the length
3665 for (Int_t irow=0;irow<160;irow++){
3666 if (track->GetClusterIndex2(irow)<0) continue;
3668 for (Int_t ipoint=0;ipoint<3;ipoint++){
3669 if (index<=ipos[ipoint]) row[ipoint] = irow;
3673 //Get cluster and sector position
3674 for (Int_t ipoint=0;ipoint<3;ipoint++){
3675 Int_t clindex = track->GetClusterIndex2(row[ipoint]);
3676 AliTPCclusterMI * cl = GetClusterMI(clindex);
3679 // AliTPCclusterMI * cl = GetClusterMI(clindex);
3682 sec[ipoint] = ((clindex&0xff000000)>>24)%18;
3683 xyz[ipoint][0] = GetXrow(row[ipoint]);
3684 xyz[ipoint][1] = cl->GetY();
3685 xyz[ipoint][2] = cl->GetZ();
3689 // Calculate seed state vector and covariance matrix
3691 Double_t alpha, cs,sn, xx2,yy2;
3693 alpha = (sec[1]-sec[2])*fSectors->GetAlpha();
3694 cs = TMath::Cos(alpha);
3695 sn = TMath::Sin(alpha);
3696 xx2= xyz[1][0]*cs-xyz[1][1]*sn;
3697 yy2= xyz[1][0]*sn+xyz[1][1]*cs;
3701 alpha = (sec[0]-sec[2])*fSectors->GetAlpha();
3702 cs = TMath::Cos(alpha);
3703 sn = TMath::Sin(alpha);
3704 xx2= xyz[0][0]*cs-xyz[0][1]*sn;
3705 yy2= xyz[0][0]*sn+xyz[0][1]*cs;
3711 Double_t x[5],c[15];
3715 x[4]=F1(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
3716 x[2]=F2(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
3717 x[3]=F3n(xyz[2][0],xyz[2][1],xyz[0][0],xyz[0][1],xyz[2][2],xyz[0][2],x[4]);
3719 Double_t sy =0.1, sz =0.1;
3721 Double_t sy1=0.2, sz1=0.2;
3722 Double_t sy2=0.2, sz2=0.2;
3725 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;
3726 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;
3727 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;
3728 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;
3729 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;
3730 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;
3732 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;
3733 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;
3734 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;
3735 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;
3740 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3741 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3742 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3743 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3744 c[13]=f30*sy1*f40+f32*sy2*f42;
3745 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3747 // Int_t row1 = fSectors->GetRowNumber(xyz[2][0]);
3748 AliTPCseed *seed=new AliTPCseed(xyz[2][0], sec[2]*fSectors->GetAlpha()+fSectors->GetAlphaShift(), x, c, 0);
3749 seed->SetLastPoint(row[2]);
3750 seed->SetFirstPoint(row[2]);
3755 AliTPCseed *AliTPCtrackerMI::ReSeed(AliTPCseed *track,Int_t r0, Bool_t forward)
3759 //reseed using founded clusters
3762 Int_t row[3]={0,0,0};
3763 Int_t sec[3]={0,0,0};
3765 // forward direction
3767 for (Int_t irow=r0;irow<160;irow++){
3768 if (track->GetClusterIndex(irow)>0){
3773 for (Int_t irow=160;irow>r0;irow--){
3774 if (track->GetClusterIndex(irow)>0){
3779 for (Int_t irow=row[2]-15;irow>row[0];irow--){
3780 if (track->GetClusterIndex(irow)>0){
3788 for (Int_t irow=0;irow<r0;irow++){
3789 if (track->GetClusterIndex(irow)>0){
3794 for (Int_t irow=r0;irow>0;irow--){
3795 if (track->GetClusterIndex(irow)>0){
3800 for (Int_t irow=row[2]-15;irow>row[0];irow--){
3801 if (track->GetClusterIndex(irow)>0){
3808 if ((row[2]-row[0])<20) return 0;
3809 if (row[1]==0) return 0;
3812 //Get cluster and sector position
3813 for (Int_t ipoint=0;ipoint<3;ipoint++){
3814 Int_t clindex = track->GetClusterIndex2(row[ipoint]);
3815 AliTPCclusterMI * cl = GetClusterMI(clindex);
3818 // AliTPCclusterMI * cl = GetClusterMI(clindex);
3821 sec[ipoint] = ((clindex&0xff000000)>>24)%18;
3822 xyz[ipoint][0] = GetXrow(row[ipoint]);
3823 AliTPCTrackerPoint * point = track->GetTrackPoint(row[ipoint]);
3824 if (point&&ipoint<2){
3826 xyz[ipoint][1] = point->GetY();
3827 xyz[ipoint][2] = point->GetZ();
3830 xyz[ipoint][1] = cl->GetY();
3831 xyz[ipoint][2] = cl->GetZ();
3838 // Calculate seed state vector and covariance matrix
3840 Double_t alpha, cs,sn, xx2,yy2;
3842 alpha = (sec[1]-sec[2])*fSectors->GetAlpha();
3843 cs = TMath::Cos(alpha);
3844 sn = TMath::Sin(alpha);
3845 xx2= xyz[1][0]*cs-xyz[1][1]*sn;
3846 yy2= xyz[1][0]*sn+xyz[1][1]*cs;
3850 alpha = (sec[0]-sec[2])*fSectors->GetAlpha();
3851 cs = TMath::Cos(alpha);
3852 sn = TMath::Sin(alpha);
3853 xx2= xyz[0][0]*cs-xyz[0][1]*sn;
3854 yy2= xyz[0][0]*sn+xyz[0][1]*cs;
3860 Double_t x[5],c[15];
3864 x[4]=F1(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
3865 x[2]=F2(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
3866 x[3]=F3n(xyz[2][0],xyz[2][1],xyz[0][0],xyz[0][1],xyz[2][2],xyz[0][2],x[4]);
3868 Double_t sy =0.1, sz =0.1;
3870 Double_t sy1=0.2, sz1=0.2;
3871 Double_t sy2=0.2, sz2=0.2;
3874 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;
3875 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;
3876 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;
3877 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;
3878 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;
3879 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;
3881 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;
3882 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;
3883 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;
3884 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;
3889 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3890 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3891 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3892 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3893 c[13]=f30*sy1*f40+f32*sy2*f42;
3894 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3896 // Int_t row1 = fSectors->GetRowNumber(xyz[2][0]);
3897 AliTPCseed *seed=new AliTPCseed(xyz[2][0], sec[2]*fSectors->GetAlpha()+fSectors->GetAlphaShift(), x, c, 0);
3898 seed->SetLastPoint(row[2]);
3899 seed->SetFirstPoint(row[2]);
3900 for (Int_t i=row[0];i<row[2];i++){
3901 seed->SetClusterIndex(i, track->GetClusterIndex(i));
3909 void AliTPCtrackerMI::FindMultiMC(TObjArray * array, AliESDEvent */*esd*/, Int_t iter)
3912 // find multi tracks - THIS FUNCTION IS ONLY FOR DEBUG PURPOSES
3914 // Use AliTPCReconstructor::StreamLevel()>2 if you want to tune parameters - cuts
3916 // Two reasons to have multiple find tracks
3917 // 1. Curling tracks can be find more than once
3918 // 2. Splitted tracks
3919 // a.) Multiple seeding to increase tracking efficiency - (~ 100% reached)
3920 // b.) Edge effect on the sector boundaries
3923 // Algorithm done in 2 phases - because of CPU consumption
3924 // it is n^2 algorithm - for lead-lead 20000x20000 combination are investigated
3926 // Algorihm for curling tracks sign:
3927 // 1 phase -makes a very rough fast cuts to minimize combinatorics
3928 // a.) opposite sign
3929 // b.) one of the tracks - not pointing to the primary vertex -
3930 // c.) delta tan(theta)
3932 // 2 phase - calculates DCA between tracks - time consument
3937 // General cuts - for splitted tracks and for curling tracks
3939 const Float_t kMaxdPhi = 0.2; // maximal distance in phi
3941 // Curling tracks cuts
3946 Int_t nentries = array->GetEntriesFast();
3947 AliHelix *helixes = new AliHelix[nentries];
3948 Float_t *xm = new Float_t[nentries];
3949 Float_t *dz0 = new Float_t[nentries];
3950 Float_t *dz1 = new Float_t[nentries];
3956 // Find track COG in x direction - point with best defined parameters
3958 for (Int_t i=0;i<nentries;i++){
3959 AliTPCseed* track = (AliTPCseed*)array->At(i);
3960 if (!track) continue;
3961 track->SetCircular(0);
3962 new (&helixes[i]) AliHelix(*track);
3966 track->GetDZ(GetX(),GetY(),GetZ(),GetBz(),dz);
3969 for (Int_t icl=0; icl<160; icl++){
3970 AliTPCclusterMI * cl = track->GetClusterPointer(icl);
3976 if (ncl>0) xm[i]/=Float_t(ncl);
3978 TTreeSRedirector &cstream = *fDebugStreamer;
3980 for (Int_t i0=0;i0<nentries;i0++){
3981 AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
3982 if (!track0) continue;
3983 Float_t xc0 = helixes[i0].GetHelix(6);
3984 Float_t yc0 = helixes[i0].GetHelix(7);
3985 Float_t r0 = helixes[i0].GetHelix(8);
3986 Float_t rc0 = TMath::Sqrt(xc0*xc0+yc0*yc0);
3987 Float_t fi0 = TMath::ATan2(yc0,xc0);
3989 for (Int_t i1=i0+1;i1<nentries;i1++){
3990 AliTPCseed * track1 = (AliTPCseed*)array->At(i1);
3991 if (!track1) continue;
3992 Int_t lab0=track0->GetLabel();
3993 Int_t lab1=track1->GetLabel();
3994 if (TMath::Abs(lab0)!=TMath::Abs(lab1)) continue;
3996 Float_t xc1 = helixes[i1].GetHelix(6);
3997 Float_t yc1 = helixes[i1].GetHelix(7);
3998 Float_t r1 = helixes[i1].GetHelix(8);
3999 Float_t rc1 = TMath::Sqrt(xc1*xc1+yc1*yc1);
4000 Float_t fi1 = TMath::ATan2(yc1,xc1);
4002 Float_t dfi = fi0-fi1;
4005 if (dfi>1.5*TMath::Pi()) dfi-=TMath::Pi(); // take care about edge effect
4006 if (dfi<-1.5*TMath::Pi()) dfi+=TMath::Pi(); //
4007 if (TMath::Abs(dfi)>kMaxdPhi&&helixes[i0].GetHelix(4)*helixes[i1].GetHelix(4)<0){
4009 // if short tracks with undefined sign
4010 fi1 = -TMath::ATan2(yc1,-xc1);
4013 Float_t dtheta = TMath::Abs(track0->GetTgl()-track1->GetTgl())<TMath::Abs(track0->GetTgl()+track1->GetTgl())? track0->GetTgl()-track1->GetTgl():track0->GetTgl()+track1->GetTgl();
4016 // debug stream to tune "fast cuts"
4018 Double_t dist[3]; // distance at X
4019 Double_t mdist[3]={0,0,0}; // mean distance X+-40cm
4020 track0->GetDistance(track1,0.5*(xm[i0]+xm[i1])-40.,dist,AliTracker::GetBz());
4021 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[i]);
4022 track0->GetDistance(track1,0.5*(xm[i0]+xm[i1])+40.,dist,AliTracker::GetBz());
4023 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[i]);
4024 track0->GetDistance(track1,0.5*(xm[i0]+xm[i1]),dist,AliTracker::GetBz());
4025 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[i]);
4026 for (Int_t i=0;i<3;i++) mdist[i]*=0.33333;
4030 for (Int_t icl=0; icl<160; icl++){
4031 AliTPCclusterMI * cl0 = track0->GetClusterPointer(icl);
4032 AliTPCclusterMI * cl1 = track1->GetClusterPointer(icl);
4035 if (cl0==cl1) sums++;
4043 "Tr0.="<<track0<< // seed0
4044 "Tr1.="<<track1<< // seed1
4045 "h0.="<<&helixes[i0]<<
4046 "h1.="<<&helixes[i1]<<
4048 "sum="<<sum<< //the sum of rows with cl in both
4049 "sums="<<sums<< //the sum of shared clusters
4050 "xm0="<<xm[i0]<< // the center of track
4051 "xm1="<<xm[i1]<< // the x center of track
4052 // General cut variables
4053 "dfi="<<dfi<< // distance in fi angle
4054 "dtheta="<<dtheta<< // distance int theta angle
4060 "dist0="<<dist[0]<< //distance x
4061 "dist1="<<dist[1]<< //distance y
4062 "dist2="<<dist[2]<< //distance z
4063 "mdist0="<<mdist[0]<< //distance x
4064 "mdist1="<<mdist[1]<< //distance y
4065 "mdist2="<<mdist[2]<< //distance z
4078 if (AliTPCReconstructor::StreamLevel()>1) {
4079 AliInfo("Time for curling tracks removal DEBUGGING MC");
4085 void AliTPCtrackerMI::FindSplitted(TObjArray * array, AliESDEvent */*esd*/, Int_t iter)
4089 // Two reasons to have multiple find tracks
4090 // 1. Curling tracks can be find more than once
4091 // 2. Splitted tracks
4092 // a.) Multiple seeding to increase tracking efficiency - (~ 100% reached)
4093 // b.) Edge effect on the sector boundaries
4095 // This function tries to find tracks closed in the parametric space
4097 // cut logic if distance is bigger than cut continue - Do Nothing
4098 const Float_t kMaxdTheta = 0.05; // maximal distance in theta
4099 const Float_t kMaxdPhi = 0.05; // maximal deistance in phi
4100 const Float_t kdelta = 40.; //delta r to calculate track distance
4102 // const Float_t kMaxDist0 = 1.; // maximal distance 0
4103 //const Float_t kMaxDist1 = 0.3; // maximal distance 1 - cut if track in separate rows
4106 TCut csec("csec","abs(Tr0.fRelativeSector-Tr1.fRelativeSector)<2");
4107 TCut cdtheta("cdtheta","abs(dtheta)<0.05");
4112 Int_t nentries = array->GetEntriesFast();
4113 AliHelix *helixes = new AliHelix[nentries];
4114 Float_t *xm = new Float_t[nentries];
4120 //Sort tracks according quality
4122 Int_t nseed = array->GetEntriesFast();
4123 Float_t * quality = new Float_t[nseed];
4124 Int_t * indexes = new Int_t[nseed];
4125 for (Int_t i=0; i<nseed; i++) {
4126 AliTPCseed *pt=(AliTPCseed*)array->UncheckedAt(i);
4131 pt->UpdatePoints(); //select first last max dens points
4132 Float_t * points = pt->GetPoints();
4133 if (points[3]<0.8) quality[i] =-1;
4134 quality[i] = (points[2]-points[0])+pt->GetNumberOfClusters();
4135 //prefer high momenta tracks if overlaps
4136 quality[i] *= TMath::Sqrt(TMath::Abs(pt->Pt())+0.5);
4138 TMath::Sort(nseed,quality,indexes);
4142 // Find track COG in x direction - point with best defined parameters
4144 for (Int_t i=0;i<nentries;i++){
4145 AliTPCseed* track = (AliTPCseed*)array->At(i);
4146 if (!track) continue;
4147 track->SetCircular(0);
4148 new (&helixes[i]) AliHelix(*track);
4151 for (Int_t icl=0; icl<160; icl++){
4152 AliTPCclusterMI * cl = track->GetClusterPointer(icl);
4158 if (ncl>0) xm[i]/=Float_t(ncl);
4160 TTreeSRedirector &cstream = *fDebugStreamer;
4162 for (Int_t is0=0;is0<nentries;is0++){
4163 Int_t i0 = indexes[is0];
4164 AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
4165 if (!track0) continue;
4166 if (track0->GetKinkIndexes()[0]!=0) continue;
4167 Float_t xc0 = helixes[i0].GetHelix(6);
4168 Float_t yc0 = helixes[i0].GetHelix(7);
4169 Float_t fi0 = TMath::ATan2(yc0,xc0);
4171 for (Int_t is1=is0+1;is1<nentries;is1++){
4172 Int_t i1 = indexes[is1];
4173 AliTPCseed * track1 = (AliTPCseed*)array->At(i1);
4174 if (!track1) continue;
4176 if (TMath::Abs(track0->GetRelativeSector()-track1->GetRelativeSector())>1) continue;
4177 if (track1->GetKinkIndexes()[0]>0 &&track0->GetKinkIndexes()[0]<0) continue;
4178 if (track1->GetKinkIndexes()[0]!=0) continue;
4180 Float_t dtheta = TMath::Abs(track0->GetTgl()-track1->GetTgl())<TMath::Abs(track0->GetTgl()+track1->GetTgl())? track0->GetTgl()-track1->GetTgl():track0->GetTgl()+track1->GetTgl();
4181 if (TMath::Abs(dtheta)>kMaxdTheta) continue;
4183 Float_t xc1 = helixes[i1].GetHelix(6);
4184 Float_t yc1 = helixes[i1].GetHelix(7);
4185 Float_t fi1 = TMath::ATan2(yc1,xc1);
4187 Float_t dfi = fi0-fi1;
4188 if (dfi>1.5*TMath::Pi()) dfi-=TMath::Pi(); // take care about edge effect
4189 if (dfi<-1.5*TMath::Pi()) dfi+=TMath::Pi(); //
4190 if (TMath::Abs(dfi)>kMaxdPhi&&helixes[i0].GetHelix(4)*helixes[i1].GetHelix(4)<0){
4192 // if short tracks with undefined sign
4193 fi1 = -TMath::ATan2(yc1,-xc1);
4196 if (TMath::Abs(dfi)>kMaxdPhi) continue;
4203 for (Int_t icl=0; icl<160; icl++){
4204 Int_t index0=track0->GetClusterIndex2(icl);
4205 Int_t index1=track1->GetClusterIndex2(icl);
4206 Bool_t used0 = (index0>0 && !(index0&0x8000));
4207 Bool_t used1 = (index1>0 && !(index1&0x8000));
4209 if (used0) sum0++; // used cluster0
4210 if (used1) sum1++; // used clusters1
4211 if (used0&&used1) sum++;
4212 if (index0==index1 && used0 && used1) sums++;
4216 if (sums<10) continue;
4217 if (sum<40) continue;
4218 if (sums/Float_t(TMath::Min(sum0,sum1))<0.5) continue;
4220 Double_t dist[5][4]; // distance at X
4221 Double_t mdist[4]={0,0,0,0}; // mean distance on range +-delta
4225 track0->GetDistance(track1,xm[i0],dist[0],AliTracker::GetBz());
4226 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[0][i]);
4227 track0->GetDistance(track1,xm[i1],dist[1],AliTracker::GetBz());
4228 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[1][i]);
4230 track0->GetDistance(track1,TMath::Min(xm[i1],xm[i0])-kdelta,dist[2],AliTracker::GetBz());
4231 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[2][i]);
4232 track0->GetDistance(track1,TMath::Max(xm[i1],xm[i0])+kdelta,dist[3],AliTracker::GetBz());
4233 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[3][i]);
4235 track0->GetDistance(track1,(xm[i1]+xm[i0])*0.5,dist[4],AliTracker::GetBz());
4236 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[4][i]);
4237 for (Int_t i=0;i<3;i++) mdist[i]*=0.2;
4240 Int_t lab0=track0->GetLabel();
4241 Int_t lab1=track1->GetLabel();
4242 cstream<<"Splitted"<<
4246 "Tr0.="<<track0<< // seed0
4247 "Tr1.="<<track1<< // seed1
4248 "h0.="<<&helixes[i0]<<
4249 "h1.="<<&helixes[i1]<<
4251 "sum="<<sum<< //the sum of rows with cl in both
4252 "sum0="<<sum0<< //the sum of rows with cl in 0 track
4253 "sum1="<<sum1<< //the sum of rows with cl in 1 track
4254 "sums="<<sums<< //the sum of shared clusters
4255 "xm0="<<xm[i0]<< // the center of track
4256 "xm1="<<xm[i1]<< // the x center of track
4257 // General cut variables
4258 "dfi="<<dfi<< // distance in fi angle
4259 "dtheta="<<dtheta<< // distance int theta angle
4262 "dist0="<<dist[4][0]<< //distance x
4263 "dist1="<<dist[4][1]<< //distance y
4264 "dist2="<<dist[4][2]<< //distance z
4265 "mdist0="<<mdist[0]<< //distance x
4266 "mdist1="<<mdist[1]<< //distance y
4267 "mdist2="<<mdist[2]<< //distance z
4270 delete array->RemoveAt(i1);
4275 AliInfo("Time for splitted tracks removal");
4281 void AliTPCtrackerMI::FindCurling(TObjArray * array, AliESDEvent */*esd*/, Int_t iter)
4284 // find Curling tracks
4285 // Use AliTPCReconstructor::StreamLevel()>1 if you want to tune parameters - cuts
4288 // Algorithm done in 2 phases - because of CPU consumption
4289 // it is n^2 algorithm - for lead-lead 20000x20000 combination are investigated
4290 // see detal in MC part what can be used to cut
4294 const Float_t kMaxC = 400; // maximal curvature to of the track
4295 const Float_t kMaxdTheta = 0.15; // maximal distance in theta
4296 const Float_t kMaxdPhi = 0.15; // maximal distance in phi
4297 const Float_t kPtRatio = 0.3; // ratio between pt
4298 const Float_t kMinDCAR = 2.; // distance to the primary vertex in r - see cpipe cut
4301 // Curling tracks cuts
4304 const Float_t kMaxDeltaRMax = 40; // distance in outer radius
4305 const Float_t kMaxDeltaRMin = 5.; // distance in lower radius - see cpipe cut
4306 const Float_t kMinAngle = 2.9; // angle between tracks
4307 const Float_t kMaxDist = 5; // biggest distance
4309 // The cuts can be tuned using the "MC information stored in Multi tree ==> see FindMultiMC
4312 TCut csign("csign","Tr0.fP[4]*Tr1.fP[4]<0"); //opposite sign
4313 TCut cmax("cmax","abs(Tr0.GetC())>1/400");
4314 TCut cda("cda","sqrt(dtheta^2+dfi^2)<0.15");
4315 TCut ccratio("ccratio","abs((Tr0.fP[4]+Tr1.fP[4])/(abs(Tr0.fP[4])+abs(Tr1.fP[4])))<0.3");
4316 TCut cpipe("cpipe", "min(abs(r0-rc0),abs(r1-rc1))>5");
4318 TCut cdrmax("cdrmax","abs(abs(rc0+r0)-abs(rc1+r1))<40")
4319 TCut cdrmin("cdrmin","abs(abs(rc0+r0)-abs(rc1+r1))<10")
4321 Multi->Draw("dfi","iter==0"+csign+cmax+cda+ccratio); ~94% of curling tracks fulfill
4322 Multi->Draw("min(abs(r0-rc0),abs(r1-rc1))","iter==0&&abs(lab1)==abs(lab0)"+csign+cmax+cda+ccratio+cpipe+cdrmin+cdrmax); //80%
4324 Curling2->Draw("dfi","iter==0&&abs(lab0)==abs(lab1)"+csign+cmax+cdtheta+cdfi+ccratio)
4330 Int_t nentries = array->GetEntriesFast();
4331 AliHelix *helixes = new AliHelix[nentries];
4332 for (Int_t i=0;i<nentries;i++){
4333 AliTPCseed* track = (AliTPCseed*)array->At(i);
4334 if (!track) continue;
4335 track->SetCircular(0);
4336 new (&helixes[i]) AliHelix(*track);
4342 Double_t phase[2][2],radius[2];
4346 TTreeSRedirector &cstream = *fDebugStreamer;
4348 for (Int_t i0=0;i0<nentries;i0++){
4349 AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
4350 if (!track0) continue;
4351 if (TMath::Abs(track0->GetC())<1/kMaxC) continue;
4352 Float_t xc0 = helixes[i0].GetHelix(6);
4353 Float_t yc0 = helixes[i0].GetHelix(7);
4354 Float_t r0 = helixes[i0].GetHelix(8);
4355 Float_t rc0 = TMath::Sqrt(xc0*xc0+yc0*yc0);
4356 Float_t fi0 = TMath::ATan2(yc0,xc0);
4358 for (Int_t i1=i0+1;i1<nentries;i1++){
4359 AliTPCseed * track1 = (AliTPCseed*)array->At(i1);
4360 if (!track1) continue;
4361 if (TMath::Abs(track1->GetC())<1/kMaxC) continue;
4362 Float_t xc1 = helixes[i1].GetHelix(6);
4363 Float_t yc1 = helixes[i1].GetHelix(7);
4364 Float_t r1 = helixes[i1].GetHelix(8);
4365 Float_t rc1 = TMath::Sqrt(xc1*xc1+yc1*yc1);
4366 Float_t fi1 = TMath::ATan2(yc1,xc1);
4368 Float_t dfi = fi0-fi1;
4371 if (dfi>1.5*TMath::Pi()) dfi-=TMath::Pi(); // take care about edge effect
4372 if (dfi<-1.5*TMath::Pi()) dfi+=TMath::Pi(); //
4373 Float_t dtheta = TMath::Abs(track0->GetTgl()-track1->GetTgl())<TMath::Abs(track0->GetTgl()+track1->GetTgl())? track0->GetTgl()-track1->GetTgl():track0->GetTgl()+track1->GetTgl();
4377 if (track0->GetBConstrain()&&track1->GetBConstrain()) continue; // not constrained
4378 if (track1->GetSigned1Pt()*track0->GetSigned1Pt()>0) continue; // not the same sign
4379 if ( TMath::Abs(track1->GetTgl()+track0->GetTgl())>kMaxdTheta) continue; //distance in the Theta
4380 if ( TMath::Abs(dfi)>kMaxdPhi) continue; //distance in phi
4381 if ( TMath::Sqrt(dfi*dfi+dtheta*dtheta)>kMaxdPhi) continue; //common angular offset
4383 Float_t pt0 = track0->GetSignedPt();
4384 Float_t pt1 = track1->GetSignedPt();
4385 if ((TMath::Abs(pt0+pt1)/(TMath::Abs(pt0)+TMath::Abs(pt1)))>kPtRatio) continue;
4386 if ((iter==1) && TMath::Abs(TMath::Abs(rc0+r0)-TMath::Abs(rc1+r1))>kMaxDeltaRMax) continue;
4387 if ((iter!=1) &&TMath::Abs(TMath::Abs(rc0-r0)-TMath::Abs(rc1-r1))>kMaxDeltaRMin) continue;
4388 if (TMath::Min(TMath::Abs(rc0-r0),TMath::Abs(rc1-r1))<kMinDCAR) continue;
4391 // Now find closest approach
4395 Int_t npoints = helixes[i0].GetRPHIintersections(helixes[i1], phase, radius,10);
4396 if (npoints==0) continue;
4397 helixes[i0].GetClosestPhases(helixes[i1], phase);
4401 Double_t hangles[3];
4402 helixes[i0].Evaluate(phase[0][0],xyz0);
4403 helixes[i1].Evaluate(phase[0][1],xyz1);
4405 helixes[i0].GetAngle(phase[0][0],helixes[i1],phase[0][1],hangles);
4406 Double_t deltah[2],deltabest;
4407 if (TMath::Abs(hangles[2])<kMinAngle) continue;
4411 helixes[i0].ParabolicDCA(helixes[i1],phase[0][0],phase[0][1],radius[0],deltah[0],2);
4413 helixes[i0].ParabolicDCA(helixes[i1],phase[1][0],phase[1][1],radius[1],deltah[1],2);
4414 if (deltah[1]<deltah[0]) ibest=1;
4416 deltabest = TMath::Sqrt(deltah[ibest]);
4417 helixes[i0].Evaluate(phase[ibest][0],xyz0);
4418 helixes[i1].Evaluate(phase[ibest][1],xyz1);
4419 helixes[i0].GetAngle(phase[ibest][0],helixes[i1],phase[ibest][1],hangles);
4420 Double_t radiusbest = TMath::Sqrt(radius[ibest]);
4422 if (deltabest>kMaxDist) continue;
4423 // if (mindcar+mindcaz<40 && (TMath::Abs(hangles[2])<kMinAngle ||deltabest>3)) continue;
4424 Bool_t sign =kFALSE;
4425 if (hangles[2]>kMinAngle) sign =kTRUE;
4428 // circular[i0] = kTRUE;
4429 // circular[i1] = kTRUE;
4430 if (track0->OneOverPt()<track1->OneOverPt()){
4431 track0->SetCircular(track0->GetCircular()+1);
4432 track1->SetCircular(track1->GetCircular()+2);
4435 track1->SetCircular(track1->GetCircular()+1);
4436 track0->SetCircular(track0->GetCircular()+2);
4439 if (AliTPCReconstructor::StreamLevel()>1){
4441 //debug stream to tune "fine" cuts
4442 Int_t lab0=track0->GetLabel();
4443 Int_t lab1=track1->GetLabel();
4444 cstream<<"Curling2"<<
4460 "npoints="<<npoints<<
4461 "hangles0="<<hangles[0]<<
4462 "hangles1="<<hangles[1]<<
4463 "hangles2="<<hangles[2]<<
4466 "radius="<<radiusbest<<
4467 "deltabest="<<deltabest<<
4468 "phase0="<<phase[ibest][0]<<
4469 "phase1="<<phase[ibest][1]<<
4477 if (AliTPCReconstructor::StreamLevel()>1) {
4478 AliInfo("Time for curling tracks removal");
4487 void AliTPCtrackerMI::FindKinks(TObjArray * array, AliESDEvent *esd)
4494 TObjArray *kinks= new TObjArray(10000);
4495 // TObjArray *v0s= new TObjArray(10000);
4496 Int_t nentries = array->GetEntriesFast();
4497 AliHelix *helixes = new AliHelix[nentries];
4498 Int_t *sign = new Int_t[nentries];
4499 Int_t *nclusters = new Int_t[nentries];
4500 Float_t *alpha = new Float_t[nentries];
4501 AliKink *kink = new AliKink();
4502 Int_t * usage = new Int_t[nentries];
4503 Float_t *zm = new Float_t[nentries];
4504 Float_t *z0 = new Float_t[nentries];
4505 Float_t *fim = new Float_t[nentries];
4506 Float_t *shared = new Float_t[nentries];
4507 Bool_t *circular = new Bool_t[nentries];
4508 Float_t *dca = new Float_t[nentries];
4509 //const AliESDVertex * primvertex = esd->GetVertex();
4511 // nentries = array->GetEntriesFast();
4516 for (Int_t i=0;i<nentries;i++){
4519 AliTPCseed* track = (AliTPCseed*)array->At(i);
4520 if (!track) continue;
4521 track->SetCircular(0);
4523 track->UpdatePoints();
4524 if (( track->GetPoints()[2]- track->GetPoints()[0])>5 && track->GetPoints()[3]>0.8){
4526 nclusters[i]=track->GetNumberOfClusters();
4527 alpha[i] = track->GetAlpha();
4528 new (&helixes[i]) AliHelix(*track);
4530 helixes[i].Evaluate(0,xyz);
4531 sign[i] = (track->GetC()>0) ? -1:1;
4534 if (track->GetProlongation(x,y,z)){
4536 fim[i] = alpha[i]+TMath::ATan2(y,x);
4539 zm[i] = track->GetZ();
4543 circular[i]= kFALSE;
4544 if (track->GetProlongation(0,y,z)) z0[i] = z;
4545 dca[i] = track->GetD(0,0);
4551 Int_t ncandidates =0;
4554 Double_t phase[2][2],radius[2];
4557 // Find circling track
4558 TTreeSRedirector &cstream = *fDebugStreamer;
4560 for (Int_t i0=0;i0<nentries;i0++){
4561 AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
4562 if (!track0) continue;
4563 if (track0->GetNumberOfClusters()<40) continue;
4564 if (TMath::Abs(1./track0->GetC())>200) continue;
4565 for (Int_t i1=i0+1;i1<nentries;i1++){
4566 AliTPCseed * track1 = (AliTPCseed*)array->At(i1);
4567 if (!track1) continue;
4568 if (track1->GetNumberOfClusters()<40) continue;
4569 if ( TMath::Abs(track1->GetTgl()+track0->GetTgl())>0.1) continue;
4570 if (track0->GetBConstrain()&&track1->GetBConstrain()) continue;
4571 if (TMath::Abs(1./track1->GetC())>200) continue;
4572 if (track1->GetSigned1Pt()*track0->GetSigned1Pt()>0) continue;
4573 if (track1->GetTgl()*track0->GetTgl()>0) continue;
4574 if (TMath::Max(TMath::Abs(1./track0->GetC()),TMath::Abs(1./track1->GetC()))>190) continue;
4575 if (track0->GetBConstrain()&&track1->OneOverPt()<track0->OneOverPt()) continue; //returning - lower momenta
4576 if (track1->GetBConstrain()&&track0->OneOverPt()<track1->OneOverPt()) continue; //returning - lower momenta
4578 Float_t mindcar = TMath::Min(TMath::Abs(dca[i0]),TMath::Abs(dca[i1]));
4579 if (mindcar<5) continue;
4580 Float_t mindcaz = TMath::Min(TMath::Abs(z0[i0]-GetZ()),TMath::Abs(z0[i1]-GetZ()));
4581 if (mindcaz<5) continue;
4582 if (mindcar+mindcaz<20) continue;
4585 Float_t xc0 = helixes[i0].GetHelix(6);
4586 Float_t yc0 = helixes[i0].GetHelix(7);
4587 Float_t r0 = helixes[i0].GetHelix(8);
4588 Float_t xc1 = helixes[i1].GetHelix(6);
4589 Float_t yc1 = helixes[i1].GetHelix(7);
4590 Float_t r1 = helixes[i1].GetHelix(8);
4592 Float_t rmean = (r0+r1)*0.5;
4593 Float_t delta =TMath::Sqrt((xc1-xc0)*(xc1-xc0)+(yc1-yc0)*(yc1-yc0));
4594 //if (delta>30) continue;
4595 if (delta>rmean*0.25) continue;
4596 if (TMath::Abs(r0-r1)/rmean>0.3) continue;
4598 Int_t npoints = helixes[i0].GetRPHIintersections(helixes[i1], phase, radius,10);
4599 if (npoints==0) continue;
4600 helixes[i0].GetClosestPhases(helixes[i1], phase);
4604 Double_t hangles[3];
4605 helixes[i0].Evaluate(phase[0][0],xyz0);
4606 helixes[i1].Evaluate(phase[0][1],xyz1);
4608 helixes[i0].GetAngle(phase[0][0],helixes[i1],phase[0][1],hangles);
4609 Double_t deltah[2],deltabest;
4610 if (hangles[2]<2.8) continue;
4613 helixes[i0].ParabolicDCA(helixes[i1],phase[0][0],phase[0][1],radius[0],deltah[0],2);
4615 helixes[i0].ParabolicDCA(helixes[i1],phase[1][0],phase[1][1],radius[1],deltah[1],2);
4616 if (deltah[1]<deltah[0]) ibest=1;
4618 deltabest = TMath::Sqrt(deltah[ibest]);
4619 helixes[i0].Evaluate(phase[ibest][0],xyz0);
4620 helixes[i1].Evaluate(phase[ibest][1],xyz1);
4621 helixes[i0].GetAngle(phase[ibest][0],helixes[i1],phase[ibest][1],hangles);
4622 Double_t radiusbest = TMath::Sqrt(radius[ibest]);
4624 if (deltabest>6) continue;
4625 if (mindcar+mindcaz<40 && (hangles[2]<3.12||deltabest>3)) continue;
4626 Bool_t sign =kFALSE;
4627 if (hangles[2]>3.06) sign =kTRUE;
4630 circular[i0] = kTRUE;
4631 circular[i1] = kTRUE;
4632 if (track0->OneOverPt()<track1->OneOverPt()){
4633 track0->SetCircular(track0->GetCircular()+1);
4634 track1->SetCircular(track1->GetCircular()+2);
4637 track1->SetCircular(track1->GetCircular()+1);
4638 track0->SetCircular(track0->GetCircular()+2);
4641 if (sign&&AliTPCReconstructor::StreamLevel()>1){
4643 Int_t lab0=track0->GetLabel();
4644 Int_t lab1=track1->GetLabel();
4645 cstream<<"Curling"<<
4652 "mindcar="<<mindcar<<
4653 "mindcaz="<<mindcaz<<
4656 "npoints="<<npoints<<
4657 "hangles0="<<hangles[0]<<
4658 "hangles2="<<hangles[2]<<
4663 "radius="<<radiusbest<<
4664 "deltabest="<<deltabest<<
4665 "phase0="<<phase[ibest][0]<<
4666 "phase1="<<phase[ibest][1]<<
4676 for (Int_t i =0;i<nentries;i++){
4677 if (sign[i]==0) continue;
4678 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
4685 Double_t cradius0 = 40*40;
4686 Double_t cradius1 = 270*270;
4689 Double_t cdist3=0.55;
4690 for (Int_t j =i+1;j<nentries;j++){
4692 if (sign[j]*sign[i]<1) continue;
4693 if ( (nclusters[i]+nclusters[j])>200) continue;
4694 if ( (nclusters[i]+nclusters[j])<80) continue;
4695 if ( TMath::Abs(zm[i]-zm[j])>60.) continue;
4696 if ( TMath::Abs(fim[i]-fim[j])>0.6 && TMath::Abs(fim[i]-fim[j])<5.7 ) continue;
4697 //AliTPCseed * track1 = (AliTPCseed*)array->At(j); Double_t phase[2][2],radius[2];
4698 Int_t npoints = helixes[i].GetRPHIintersections(helixes[j], phase, radius,20);
4699 if (npoints<1) continue;
4702 if (radius[0]<cradius0||radius[0]>cradius1) continue;
4705 if ( (radius[0]<cradius0||radius[0]>cradius1) && (radius[1]<cradius0||radius[1]>cradius1) ) continue;
4708 Double_t delta1=10000,delta2=10000;
4709 // cuts on the intersection radius
4710 helixes[i].LinearDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta1);
4711 if (radius[0]<20&&delta1<1) continue; //intersection at vertex
4712 if (radius[0]<10&&delta1<3) continue; //intersection at vertex
4714 helixes[i].LinearDCA(helixes[j],phase[1][0],phase[1][1],radius[1],delta2);
4715 if (radius[1]<20&&delta2<1) continue; //intersection at vertex
4716 if (radius[1]<10&&delta2<3) continue; //intersection at vertex
4719 Double_t distance1 = TMath::Min(delta1,delta2);
4720 if (distance1>cdist1) continue; // cut on DCA linear approximation
4722 npoints = helixes[i].GetRPHIintersections(helixes[j], phase, radius,20);
4723 helixes[i].ParabolicDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta1);
4724 if (radius[0]<20&&delta1<1) continue; //intersection at vertex
4725 if (radius[0]<10&&delta1<3) continue; //intersection at vertex
4728 helixes[i].ParabolicDCA(helixes[j],phase[1][0],phase[1][1],radius[1],delta2);
4729 if (radius[1]<20&&delta2<1) continue; //intersection at vertex
4730 if (radius[1]<10&&delta2<3) continue; //intersection at vertex
4732 distance1 = TMath::Min(delta1,delta2);
4735 rkink = TMath::Sqrt(radius[0]);
4738 rkink = TMath::Sqrt(radius[1]);
4740 if (distance1>cdist2) continue;
4743 AliTPCseed * track1 = (AliTPCseed*)array->At(j);
4746 Int_t row0 = GetRowNumber(rkink);
4747 if (row0<10) continue;
4748 if (row0>150) continue;
4751 Float_t dens00=-1,dens01=-1;
4752 Float_t dens10=-1,dens11=-1;
4754 Int_t found,foundable,shared;
4755 track0->GetClusterStatistic(0,row0-5, found, foundable,shared,kFALSE);
4756 if (foundable>5) dens00 = Float_t(found)/Float_t(foundable);
4757 track0->GetClusterStatistic(row0+5,155, found, foundable,shared,kFALSE);
4758 if (foundable>5) dens01 = Float_t(found)/Float_t(foundable);
4760 track1->GetClusterStatistic(0,row0-5, found, foundable,shared,kFALSE);
4761 if (foundable>10) dens10 = Float_t(found)/Float_t(foundable);
4762 track1->GetClusterStatistic(row0+5,155, found, foundable,shared,kFALSE);
4763 if (foundable>10) dens11 = Float_t(found)/Float_t(foundable);
4765 if (dens00<dens10 && dens01<dens11) continue;
4766 if (dens00>dens10 && dens01>dens11) continue;
4767 if (TMath::Max(dens00,dens10)<0.1) continue;
4768 if (TMath::Max(dens01,dens11)<0.3) continue;
4770 if (TMath::Min(dens00,dens10)>0.6) continue;
4771 if (TMath::Min(dens01,dens11)>0.6) continue;
4774 AliTPCseed * ktrack0, *ktrack1;
4783 if (TMath::Abs(ktrack0->GetC())>5) continue; // cut on the curvature for mother particle
4784 AliExternalTrackParam paramm(*ktrack0);
4785 AliExternalTrackParam paramd(*ktrack1);
4786 if (row0>60&&ktrack1->GetReference().GetX()>90.)new (¶md) AliExternalTrackParam(ktrack1->GetReference());
4789 kink->SetMother(paramm);
4790 kink->SetDaughter(paramd);
4793 Float_t x[3] = { kink->GetPosition()[0],kink->GetPosition()[1],kink->GetPosition()[2]};
4795 fParam->Transform0to1(x,index);
4796 fParam->Transform1to2(x,index);
4797 row0 = GetRowNumber(x[0]);
4799 if (kink->GetR()<100) continue;
4800 if (kink->GetR()>240) continue;
4801 if (kink->GetPosition()[2]/kink->GetR()>AliTPCReconstructor::GetCtgRange()) continue; //out of fiducial volume
4802 if (kink->GetDistance()>cdist3) continue;
4803 Float_t dird = kink->GetDaughterP()[0]*kink->GetPosition()[0]+kink->GetDaughterP()[1]*kink->GetPosition()[1]; // rough direction estimate
4804 if (dird<0) continue;
4806 Float_t dirm = kink->GetMotherP()[0]*kink->GetPosition()[0]+kink->GetMotherP()[1]*kink->GetPosition()[1]; // rough direction estimate
4807 if (dirm<0) continue;
4808 Float_t mpt = TMath::Sqrt(kink->GetMotherP()[0]*kink->GetMotherP()[0]+kink->GetMotherP()[1]*kink->GetMotherP()[1]);
4809 if (mpt<0.2) continue;
4812 //for high momenta momentum not defined well in first iteration
4813 Double_t qt = TMath::Sin(kink->GetAngle(2))*ktrack1->GetP();
4814 if (qt>0.35) continue;
4817 kink->SetLabel(CookLabel(ktrack0,0.4,0,row0),0);
4818 kink->SetLabel(CookLabel(ktrack1,0.4,row0,160),1);
4820 kink->SetTPCDensity(dens00,0,0);
4821 kink->SetTPCDensity(dens01,0,1);
4822 kink->SetTPCDensity(dens10,1,0);
4823 kink->SetTPCDensity(dens11,1,1);
4824 kink->SetIndex(i,0);
4825 kink->SetIndex(j,1);
4828 kink->SetTPCDensity(dens10,0,0);
4829 kink->SetTPCDensity(dens11,0,1);
4830 kink->SetTPCDensity(dens00,1,0);
4831 kink->SetTPCDensity(dens01,1,1);
4832 kink->SetIndex(j,0);
4833 kink->SetIndex(i,1);
4836 if (mpt<1||kink->GetAngle(2)>0.1){
4837 // angle and densities not defined yet
4838 if (kink->GetTPCDensityFactor()<0.8) continue;
4839 if ((2-kink->GetTPCDensityFactor())*kink->GetDistance() >0.25) continue;
4840 if (kink->GetAngle(2)*ktrack0->GetP()<0.003) continue; //too small angle
4841 if (kink->GetAngle(2)>0.2&&kink->GetTPCDensityFactor()<1.15) continue;
4842 if (kink->GetAngle(2)>0.2&&kink->GetTPCDensity(0,1)>0.05) continue;
4844 Float_t criticalangle = track0->GetSigmaSnp2()+track0->GetSigmaTgl2();
4845 criticalangle+= track1->GetSigmaSnp2()+track1->GetSigmaTgl2();
4846 criticalangle= 3*TMath::Sqrt(criticalangle);
4847 if (criticalangle>0.02) criticalangle=0.02;
4848 if (kink->GetAngle(2)<criticalangle) continue;
4851 Int_t drow = Int_t(2.+0.5/(0.05+kink->GetAngle(2))); // overlap region defined
4852 Float_t shapesum =0;
4854 for ( Int_t row = row0-drow; row<row0+drow;row++){
4855 if (row<0) continue;
4856 if (row>155) continue;
4857 if (ktrack0->GetClusterPointer(row)){
4858 AliTPCTrackerPoint *point =ktrack0->GetTrackPoint(row);
4859 shapesum+=point->GetSigmaY()+point->GetSigmaZ();
4862 if (ktrack1->GetClusterPointer(row)){
4863 AliTPCTrackerPoint *point =ktrack1->GetTrackPoint(row);
4864 shapesum+=point->GetSigmaY()+point->GetSigmaZ();
4869 kink->SetShapeFactor(-1.);
4872 kink->SetShapeFactor(shapesum/sum);
4874 // esd->AddKink(kink);
4875 kinks->AddLast(kink);
4881 // sort the kinks according quality - and refit them towards vertex
4883 Int_t nkinks = kinks->GetEntriesFast();
4884 Float_t *quality = new Float_t[nkinks];
4885 Int_t *indexes = new Int_t[nkinks];
4886 AliTPCseed *mothers = new AliTPCseed[nkinks];
4887 AliTPCseed *daughters = new AliTPCseed[nkinks];
4890 for (Int_t i=0;i<nkinks;i++){
4892 AliKink *kink = (AliKink*)kinks->At(i);
4894 // refit kinks towards vertex
4896 Int_t index0 = kink->GetIndex(0);
4897 Int_t index1 = kink->GetIndex(1);
4898 AliTPCseed * ktrack0 = (AliTPCseed*)array->At(index0);
4899 AliTPCseed * ktrack1 = (AliTPCseed*)array->At(index1);
4901 Int_t sumn=ktrack0->GetNumberOfClusters()+ktrack1->GetNumberOfClusters();
4903 // Refit Kink under if too small angle
4905 if (kink->GetAngle(2)<0.05){
4906 kink->SetTPCRow0(GetRowNumber(kink->GetR()));
4907 Int_t row0 = kink->GetTPCRow0();
4908 Int_t drow = Int_t(2.+0.5/(0.05+kink->GetAngle(2)));
4911 Int_t last = row0-drow;
4912 if (last<40) last=40;
4913 if (last<ktrack0->GetFirstPoint()+25) last = ktrack0->GetFirstPoint()+25;
4914 AliTPCseed* seed0 = ReSeed(ktrack0,last,kFALSE);
4917 Int_t first = row0+drow;
4918 if (first>130) first=130;
4919 if (first>ktrack1->GetLastPoint()-25) first = TMath::Max(ktrack1->GetLastPoint()-25,30);
4920 AliTPCseed* seed1 = ReSeed(ktrack1,first,kTRUE);
4922 if (seed0 && seed1){
4923 kink->SetStatus(1,8);
4924 if (RefitKink(*seed0,*seed1,*kink)) kink->SetStatus(1,9);
4925 row0 = GetRowNumber(kink->GetR());
4926 sumn = seed0->GetNumberOfClusters()+seed1->GetNumberOfClusters();
4927 mothers[i] = *seed0;
4928 daughters[i] = *seed1;
4931 delete kinks->RemoveAt(i);
4932 if (seed0) delete seed0;
4933 if (seed1) delete seed1;
4936 if (kink->GetDistance()>0.5 || kink->GetR()<110 || kink->GetR()>240) {
4937 delete kinks->RemoveAt(i);
4938 if (seed0) delete seed0;
4939 if (seed1) delete seed1;
4947 if (kink) quality[i] = 160*((0.1+kink->GetDistance())*(2.-kink->GetTPCDensityFactor()))/(sumn+40.); //the longest -clossest will win
4949 TMath::Sort(nkinks,quality,indexes,kFALSE);
4951 //remove double find kinks
4953 for (Int_t ikink0=1;ikink0<nkinks;ikink0++){
4954 AliKink * kink0 = (AliKink*) kinks->At(indexes[ikink0]);
4955 if (!kink0) continue;
4957 for (Int_t ikink1=0;ikink1<ikink0;ikink1++){
4958 if (!kink0) continue;
4959 AliKink * kink1 = (AliKink*) kinks->At(indexes[ikink1]);
4960 if (!kink1) continue;
4961 // if not close kink continue
4962 if (TMath::Abs(kink1->GetPosition()[2]-kink0->GetPosition()[2])>10) continue;
4963 if (TMath::Abs(kink1->GetPosition()[1]-kink0->GetPosition()[1])>10) continue;
4964 if (TMath::Abs(kink1->GetPosition()[0]-kink0->GetPosition()[0])>10) continue;
4966 AliTPCseed &mother0 = mothers[indexes[ikink0]];
4967 AliTPCseed &daughter0 = daughters[indexes[ikink0]];
4968 AliTPCseed &mother1 = mothers[indexes[ikink1]];
4969 AliTPCseed &daughter1 = daughters[indexes[ikink1]];
4970 Int_t row0 = (kink0->GetTPCRow0()+kink1->GetTPCRow0())/2;
4979 for (Int_t i=0;i<row0;i++){
4980 if (mother0.GetClusterIndex(i)>0 && mother1.GetClusterIndex(i)>0){
4983 if (mother0.GetClusterIndex(i)==mother1.GetClusterIndex(i)){
4990 for (Int_t i=row0;i<158;i++){
4991 if (daughter0.GetClusterIndex(i)>0 && daughter0.GetClusterIndex(i)>0){
4994 if (mother0.GetClusterIndex(i)==mother1.GetClusterIndex(i)){
5000 Float_t ratio = Float_t(same+1)/Float_t(both+1);
5001 Float_t ratiom = Float_t(samem+1)/Float_t(bothm+1);
5002 Float_t ratiod = Float_t(samed+1)/Float_t(bothd+1);
5003 if (ratio>0.3 && ratiom>0.5 &&ratiod>0.5) {
5004 Int_t sum0 = mother0.GetNumberOfClusters()+daughter0.GetNumberOfClusters();
5005 Int_t sum1 = mother1.GetNumberOfClusters()+daughter1.GetNumberOfClusters();
5007 shared[kink0->GetIndex(0)]= kTRUE;
5008 shared[kink0->GetIndex(1)]= kTRUE;
5009 delete kinks->RemoveAt(indexes[ikink0]);
5012 shared[kink1->GetIndex(0)]= kTRUE;
5013 shared[kink1->GetIndex(1)]= kTRUE;
5014 delete kinks->RemoveAt(indexes[ikink1]);
5021 for (Int_t i=0;i<nkinks;i++){
5022 AliKink * kink = (AliKink*) kinks->At(indexes[i]);
5023 if (!kink) continue;
5024 kink->SetTPCRow0(GetRowNumber(kink->GetR()));
5025 Int_t index0 = kink->GetIndex(0);
5026 Int_t index1 = kink->GetIndex(1);
5027 if (circular[index0]||circular[index1]&&kink->GetDistance()>0.2) continue;
5028 kink->SetMultiple(usage[index0],0);
5029 kink->SetMultiple(usage[index1],1);
5030 if (kink->GetMultiple()[0]+kink->GetMultiple()[1]>2) continue;
5031 if (kink->GetMultiple()[0]+kink->GetMultiple()[1]>0 && quality[indexes[i]]>0.2) continue;
5032 if (kink->GetMultiple()[0]+kink->GetMultiple()[1]>0 && kink->GetDistance()>0.2) continue;
5033 if (circular[index0]||circular[index1]&&kink->GetDistance()>0.1) continue;
5035 AliTPCseed * ktrack0 = (AliTPCseed*)array->At(index0);
5036 AliTPCseed * ktrack1 = (AliTPCseed*)array->At(index1);
5037 if (!ktrack0 || !ktrack1) continue;
5038 Int_t index = esd->AddKink(kink);
5041 if ( ktrack0->GetKinkIndex(0)==0 && ktrack1->GetKinkIndex(0)==0) { //best kink
5042 if (mothers[indexes[i]].GetNumberOfClusters()>20 && daughters[indexes[i]].GetNumberOfClusters()>20 && (mothers[indexes[i]].GetNumberOfClusters()+daughters[indexes[i]].GetNumberOfClusters())>100){
5043 *ktrack0 = mothers[indexes[i]];
5044 *ktrack1 = daughters[indexes[i]];
5048 ktrack0->SetKinkIndex(usage[index0],-(index+1));
5049 ktrack1->SetKinkIndex(usage[index1], (index+1));
5054 // Remove tracks corresponding to shared kink's
5056 for (Int_t i=0;i<nentries;i++){
5057 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5058 if (!track0) continue;
5059 if (track0->GetKinkIndex(0)!=0) continue;
5060 if (shared[i]) delete array->RemoveAt(i);
5065 RemoveUsed2(array,0.5,0.4,30);
5067 for (Int_t i=0;i<nentries;i++){
5068 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5069 if (!track0) continue;
5070 track0->CookdEdx(0.02,0.6);
5074 for (Int_t i=0;i<nentries;i++){
5075 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5076 if (!track0) continue;
5077 if (track0->Pt()<1.4) continue;
5078 //remove double high momenta tracks - overlapped with kink candidates
5081 for (Int_t icl=track0->GetFirstPoint();icl<track0->GetLastPoint(); icl++){
5082 if (track0->GetClusterPointer(icl)!=0){
5084 if (track0->GetClusterPointer(icl)->IsUsed(10)) shared++;
5087 if (Float_t(shared+1)/Float_t(all+1)>0.5) {
5088 delete array->RemoveAt(i);
5092 if (track0->GetKinkIndex(0)!=0) continue;
5093 if (track0->GetNumberOfClusters()<80) continue;
5095 AliTPCseed *pmother = new AliTPCseed();
5096 AliTPCseed *pdaughter = new AliTPCseed();
5097 AliKink *pkink = new AliKink;
5099 AliTPCseed & mother = *pmother;
5100 AliTPCseed & daughter = *pdaughter;
5101 AliKink & kink = *pkink;
5102 if (CheckKinkPoint(track0,mother,daughter, kink)){
5103 if (mother.GetNumberOfClusters()<30||daughter.GetNumberOfClusters()<20) {
5107 continue; //too short tracks
5109 if (mother.Pt()<1.4) {
5115 Int_t row0= kink.GetTPCRow0();
5116 if (kink.GetDistance()>0.5 || kink.GetR()<110. || kink.GetR()>240.) {
5123 Int_t index = esd->AddKink(&kink);
5124 mother.SetKinkIndex(0,-(index+1));
5125 daughter.SetKinkIndex(0,index+1);
5126 if (mother.GetNumberOfClusters()>50) {
5127 delete array->RemoveAt(i);
5128 array->AddAt(new AliTPCseed(mother),i);
5131 array->AddLast(new AliTPCseed(mother));
5133 array->AddLast(new AliTPCseed(daughter));
5134 for (Int_t icl=0;icl<row0;icl++) {
5135 if (mother.GetClusterPointer(icl)) mother.GetClusterPointer(icl)->Use(20);
5138 for (Int_t icl=row0;icl<158;icl++) {
5139 if (daughter.GetClusterPointer(icl)) daughter.GetClusterPointer(icl)->Use(20);
5148 delete [] daughters;
5170 printf("Ncandidates=\t%d\t%d\t%d\t%d\n",esd->GetNumberOfKinks(),ncandidates,ntracks,nall);
5174 void AliTPCtrackerMI::FindV0s(TObjArray * array, AliESDEvent *esd)
5180 TObjArray *tpcv0s = new TObjArray(100000);
5181 Int_t nentries = array->GetEntriesFast();
5182 AliHelix *helixes = new AliHelix[nentries];
5183 Int_t *sign = new Int_t[nentries];
5184 Float_t *alpha = new Float_t[nentries];
5185 Float_t *z0 = new Float_t[nentries];
5186 Float_t *dca = new Float_t[nentries];
5187 Float_t *sdcar = new Float_t[nentries];
5188 Float_t *cdcar = new Float_t[nentries];
5189 Float_t *pulldcar = new Float_t[nentries];
5190 Float_t *pulldcaz = new Float_t[nentries];
5191 Float_t *pulldca = new Float_t[nentries];
5192 Bool_t *isPrim = new Bool_t[nentries];
5193 const AliESDVertex * primvertex = esd->GetVertex();
5194 Double_t zvertex = primvertex->GetZv();
5196 // nentries = array->GetEntriesFast();
5198 for (Int_t i=0;i<nentries;i++){
5201 AliTPCseed* track = (AliTPCseed*)array->At(i);
5202 if (!track) continue;
5203 track->GetV0Indexes()[0] = 0; //rest v0 indexes
5204 track->GetV0Indexes()[1] = 0; //rest v0 indexes
5205 track->GetV0Indexes()[2] = 0; //rest v0 indexes
5207 alpha[i] = track->GetAlpha();
5208 new (&helixes[i]) AliHelix(*track);
5210 helixes[i].Evaluate(0,xyz);
5211 sign[i] = (track->GetC()>0) ? -1:1;
5215 if (track->GetProlongation(0,y,z)) z0[i] = z;
5216 dca[i] = track->GetD(0,0);
5218 // dca error parrameterezation + pulls
5220 sdcar[i] = TMath::Sqrt(0.150*0.150+(100*track->GetC())*(100*track->GetC()));
5221 if (TMath::Abs(track->GetTgl())>1) sdcar[i]*=2.5;
5222 cdcar[i] = TMath::Exp((TMath::Abs(track->GetC())-0.0106)*525.3);
5223 pulldcar[i] = (dca[i]-cdcar[i])/sdcar[i];
5224 pulldcaz[i] = (z0[i]-zvertex)/sdcar[i];
5225 pulldca[i] = TMath::Sqrt(pulldcar[i]*pulldcar[i]+pulldcaz[i]*pulldcaz[i]);
5226 if (track->TPCrPID(1)+track->TPCrPID(2)+track->TPCrPID(3)>0.5) {
5227 if (pulldca[i]<3.) isPrim[i]=kTRUE; //pion, muon and Kaon 3 sigma cut
5229 if (track->TPCrPID(4)>0.5) {
5230 if (pulldca[i]<0.5) isPrim[i]=kTRUE; //proton 0.5 sigma cut
5232 if (track->TPCrPID(0)>0.4) {
5233 isPrim[i]=kFALSE; //electron no sigma cut
5240 Int_t ncandidates =0;
5243 Double_t phase[2][2],radius[2];
5249 TTreeSRedirector &cstream = *fDebugStreamer;
5250 Float_t fprimvertex[3]={GetX(),GetY(),GetZ()};
5252 Double_t cradius0 = 10*10;
5253 Double_t cradius1 = 200*200;
5256 Double_t cpointAngle = 0.95;
5258 Double_t delta[2]={10000,10000};
5259 for (Int_t i =0;i<nentries;i++){
5260 if (sign[i]==0) continue;
5261 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5262 if (!track0) continue;
5263 if (AliTPCReconstructor::StreamLevel()>1){
5268 "zvertex="<<zvertex<<
5269 "sdcar0="<<sdcar[i]<<
5270 "cdcar0="<<cdcar[i]<<
5271 "pulldcar0="<<pulldcar[i]<<
5272 "pulldcaz0="<<pulldcaz[i]<<
5273 "pulldca0="<<pulldca[i]<<
5274 "isPrim="<<isPrim[i]<<
5278 if (track0->GetSigned1Pt()<0) continue;
5279 if (track0->GetKinkIndex(0)>0||isPrim[i]) continue; //daughter kink
5281 if (TMath::Abs(helixes[i].GetHelix(4))<0.000000001) continue;
5286 for (Int_t j =0;j<nentries;j++){
5287 AliTPCseed * track1 = (AliTPCseed*)array->At(j);
5288 if (!track1) continue;
5289 if (track1->GetKinkIndex(0)>0 || isPrim[j]) continue; //daughter kink
5290 if (sign[j]*sign[i]>0) continue;
5291 if (TMath::Abs(helixes[j].GetHelix(4))<0.000001) continue;
5292 if (track0->GetCircular()+track1->GetCircular()>1) continue; //circling -returning track
5295 // DCA to prim vertex cut
5301 Int_t npoints = helixes[i].GetRPHIintersections(helixes[j], phase, radius,cdist2);
5302 if (npoints<1) continue;
5306 if (radius[0]<cradius0||radius[0]>cradius1) continue;
5307 helixes[i].LinearDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta[0]);
5308 if (delta[0]>cdist1) continue;
5311 if (TMath::Max(radius[0],radius[1])<cradius0|| TMath::Min(radius[0],radius[1])>cradius1) continue;
5312 helixes[i].LinearDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta[0]);
5313 helixes[i].LinearDCA(helixes[j],phase[1][0],phase[1][1],radius[1],delta[1]);
5314 if (delta[1]<delta[0]) iclosest=1;
5315 if (delta[iclosest]>cdist1) continue;
5317 helixes[i].ParabolicDCA(helixes[j],phase[iclosest][0],phase[iclosest][1],radius[iclosest],delta[iclosest]);
5318 if (radius[iclosest]<cradius0 || radius[iclosest]>cradius1 || delta[iclosest]>cdist1) continue;
5320 Double_t pointAngle = helixes[i].GetPointAngle(helixes[j],phase[iclosest],fprimvertex);
5321 if (pointAngle<cpointAngle) continue;
5323 Bool_t isGamma = kFALSE;
5324 vertex.SetParamP(*track0); //track0 - plus
5325 vertex.SetParamN(*track1); //track1 - minus
5326 vertex.Update(fprimvertex);
5327 if (track0->TPCrPID(0)>0.3&&track1->TPCrPID(0)>0.3&&vertex.GetAnglep()[2]<0.15) isGamma=kTRUE; // gamma conversion candidate
5328 Double_t pointAngle2 = vertex.GetV0CosineOfPointingAngle();
5330 if (vertex.GetV0CosineOfPointingAngle()<cpointAngle && (!isGamma)) continue;// point angle cut
5331 //Bo: if (vertex.GetDist2()>2&&(!isGamma)) continue; // point angle cut
5332 if (vertex.GetDcaV0Daughters()>2&&(!isGamma)) continue;//Bo: // point angle cut
5333 Float_t sigmae = 0.15*0.15;
5334 if (vertex.GetRr()<80)
5335 sigmae += (sdcar[i]*sdcar[i]+sdcar[j]*sdcar[j])*(1.-vertex.GetRr()/80.)*(1.-vertex.GetRr()/80.);
5336 sigmae+= TMath::Sqrt(sigmae);
5337 //Bo: if (vertex.GetDist2()/sigmae>3.&&(!isGamma)) continue;
5338 if (vertex.GetDcaV0Daughters()/sigmae>3.&&(!isGamma)) continue;
5339 Float_t densb0=0,densb1=0,densa0=0,densa1=0;
5340 Int_t row0 = GetRowNumber(vertex.GetRr());
5342 //Bo: if (vertex.GetDist2()>0.2) continue;
5343 if (vertex.GetDcaV0Daughters()>0.2) continue;
5344 densb0 = track0->Density2(0,row0-5);
5345 densb1 = track1->Density2(0,row0-5);
5346 if (densb0>0.3|| densb1>0.3) continue; //clusters before vertex
5347 densa0 = track0->Density2(row0+5,row0+40);
5348 densa1 = track1->Density2(row0+5,row0+40);
5349 if ((densa0<0.4|| densa1<0.4)&&(!isGamma)) continue; //missing clusters after vertex
5352 densa0 = track0->Density2(0,40); //cluster density
5353 densa1 = track1->Density2(0,40); //cluster density
5354 if ((vertex.GetRr()<80&&densa0+densa1<1.)&&(!isGamma)) continue;
5356 //Bo: vertex.SetLab(0,track0->GetLabel());
5357 //Bo: vertex.SetLab(1,track1->GetLabel());
5358 vertex.SetChi2After((densa0+densa1)*0.5);
5359 vertex.SetChi2Before((densb0+densb1)*0.5);
5360 vertex.SetIndex(0,i);
5361 vertex.SetIndex(1,j);
5362 //Bo: vertex.SetStatus(1); // TPC v0 candidate
5363 vertex.SetOnFlyStatus(2);//Bo: // TPC v0 candidate
5364 //Bo: vertex.SetRp(track0->TPCrPIDs());
5365 //Bo: vertex.SetRm(track1->TPCrPIDs());
5366 tpcv0s->AddLast(new AliESDv0(vertex));
5369 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
5370 Double_t radiusm= (delta[0]<delta[1])? TMath::Sqrt(radius[0]):TMath::Sqrt(radius[1]);
5371 Double_t deltam= (delta[0]<delta[1])? TMath::Sqrt(delta[0]):TMath::Sqrt(delta[1]);
5372 if (AliTPCReconstructor::StreamLevel()>1) {
5373 Int_t lab0=track0->GetLabel();
5374 Int_t lab1=track1->GetLabel();
5375 Char_t c0=track0->GetCircular();
5376 Char_t c1=track1->GetCircular();
5379 "vertex.="<<&vertex<<
5382 "Helix0.="<<&helixes[i]<<
5385 "Helix1.="<<&helixes[j]<<
5386 "pointAngle="<<pointAngle<<
5387 "pointAngle2="<<pointAngle2<<
5392 "zvertex="<<zvertex<<
5395 "npoints="<<npoints<<
5396 "radius0="<<radius[0]<<
5397 "delta0="<<delta[0]<<
5398 "radius1="<<radius[1]<<
5399 "delta1="<<delta[1]<<
5400 "radiusm="<<radiusm<<
5402 "sdcar0="<<sdcar[i]<<
5403 "sdcar1="<<sdcar[j]<<
5404 "cdcar0="<<cdcar[i]<<
5405 "cdcar1="<<cdcar[j]<<
5406 "pulldcar0="<<pulldcar[i]<<
5407 "pulldcar1="<<pulldcar[j]<<
5408 "pulldcaz0="<<pulldcaz[i]<<
5409 "pulldcaz1="<<pulldcaz[j]<<
5410 "pulldca0="<<pulldca[i]<<
5411 "pulldca1="<<pulldca[j]<<
5421 Float_t *quality = new Float_t[ncandidates];
5422 Int_t *indexes = new Int_t[ncandidates];
5424 for (Int_t i=0;i<ncandidates;i++){
5426 AliESDv0 *v0 = (AliESDv0*)tpcv0s->At(i);
5427 quality[i] = 1./(1.00001-v0->GetV0CosineOfPointingAngle()); //base point angle
5428 // quality[i] /= (0.5+v0->GetDist2());
5429 // quality[i] *= v0->GetChi2After(); //density factor
5431 Int_t index0 = v0->GetIndex(0);
5432 Int_t index1 = v0->GetIndex(1);
5433 //Bo: Double_t minpulldca = TMath::Min(2.+pulldca[v0->GetIndex(0)],(2.+pulldca[v0->GetIndex(1)]) ); //pull
5434 Double_t minpulldca = TMath::Min(2.+pulldca[index0],(2.+pulldca[index1]) );//Bo:
5438 AliTPCseed * track0 = (AliTPCseed*)array->At(index0);
5439 AliTPCseed * track1 = (AliTPCseed*)array->At(index1);
5440 if (track0->TPCrPID(0)>0.3&&track1->TPCrPID(0)>0.3&&v0->GetAnglep()[2]<0.15) quality[i]+=1000000; // gamma conversion candidate
5441 if (track0->TPCrPID(4)>0.9||track1->TPCrPID(4)>0.9&&minpulldca>4) quality[i]*=10; // lambda candidate candidate
5444 TMath::Sort(ncandidates,quality,indexes,kTRUE);
5447 for (Int_t i=0;i<ncandidates;i++){
5448 AliESDv0 * v0 = (AliESDv0*)tpcv0s->At(indexes[i]);
5450 Int_t index0 = v0->GetIndex(0);
5451 Int_t index1 = v0->GetIndex(1);
5452 AliTPCseed * track0 = (AliTPCseed*)array->At(index0);
5453 AliTPCseed * track1 = (AliTPCseed*)array->At(index1);
5454 if (!track0||!track1) {
5458 Bool_t accept =kTRUE; //default accept
5459 Int_t *v0indexes0 = track0->GetV0Indexes();
5460 Int_t *v0indexes1 = track1->GetV0Indexes();
5462 Int_t order0 = (v0indexes0[0]!=0) ? 1:0;
5463 Int_t order1 = (v0indexes1[0]!=0) ? 1:0;
5464 if (v0indexes0[1]!=0) order0 =2;
5465 if (v0indexes1[1]!=0) order1 =2;
5467 if (v0indexes0[2]!=0) {order0=3; accept=kFALSE;}
5468 if (v0indexes0[2]!=0) {order1=3; accept=kFALSE;}
5470 AliESDv0 * v02 = v0;
5472 //Bo: v0->SetOrder(0,order0);
5473 //Bo: v0->SetOrder(1,order1);
5474 //Bo: v0->SetOrder(1,order0+order1);
5475 v0->SetOnFlyStatus(kTRUE);
5476 Int_t index = esd->AddV0(v0);
5477 v02 = esd->GetV0(index);
5478 v0indexes0[order0]=index;
5479 v0indexes1[order1]=index;
5483 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
5484 if (AliTPCReconstructor::StreamLevel()>1) {
5485 Int_t lab0=track0->GetLabel();
5486 Int_t lab1=track1->GetLabel();
5495 "dca0="<<dca[index0]<<
5496 "dca1="<<dca[index1]<<
5500 "quality="<<quality[i]<<
5501 "pulldca0="<<pulldca[index0]<<
5502 "pulldca1="<<pulldca[index1]<<
5526 printf("TPC V0 finder : naccepted\t%d\tncandidates\t%d\tntracks\t%d\tnall\t%d\n",naccepted,ncandidates,ntracks,nall);
5530 Int_t AliTPCtrackerMI::RefitKink(AliTPCseed &mother, AliTPCseed &daughter, AliESDkink &knk)
5533 // refit kink towards to the vertex
5536 AliKink &kink=(AliKink &)knk;
5538 Int_t row0 = GetRowNumber(kink.GetR());
5539 FollowProlongation(mother,0);
5540 mother.Reset(kFALSE);
5542 FollowProlongation(daughter,row0);
5543 daughter.Reset(kFALSE);
5544 FollowBackProlongation(daughter,158);
5545 daughter.Reset(kFALSE);
5546 Int_t first = TMath::Max(row0-20,30);
5547 Int_t last = TMath::Min(row0+20,140);
5549 const Int_t kNdiv =5;
5550 AliTPCseed param0[kNdiv]; // parameters along the track
5551 AliTPCseed param1[kNdiv]; // parameters along the track
5552 AliKink kinks[kNdiv]; // corresponding kink parameters
5555 for (Int_t irow=0; irow<kNdiv;irow++){
5556 rows[irow] = first +((last-first)*irow)/(kNdiv-1);
5558 // store parameters along the track
5560 for (Int_t irow=0;irow<kNdiv;irow++){
5561 FollowBackProlongation(mother, rows[irow]);
5562 FollowProlongation(daughter,rows[kNdiv-1-irow]);
5563 param0[irow] = mother;
5564 param1[kNdiv-1-irow] = daughter;
5568 for (Int_t irow=0; irow<kNdiv-1;irow++){
5569 if (param0[irow].GetNumberOfClusters()<kNdiv||param1[irow].GetNumberOfClusters()<kNdiv) continue;
5570 kinks[irow].SetMother(param0[irow]);
5571 kinks[irow].SetDaughter(param1[irow]);
5572 kinks[irow].Update();
5575 // choose kink with best "quality"
5577 Double_t mindist = 10000;
5578 for (Int_t irow=0;irow<kNdiv;irow++){
5579 if (param0[irow].GetNumberOfClusters()<20||param1[irow].GetNumberOfClusters()<20) continue;
5580 if (TMath::Abs(kinks[irow].GetR())>240.) continue;
5581 if (TMath::Abs(kinks[irow].GetR())<100.) continue;
5583 Float_t normdist = TMath::Abs(param0[irow].GetX()-kinks[irow].GetR())*(0.1+kink.GetDistance());
5584 normdist/= (param0[irow].GetNumberOfClusters()+param1[irow].GetNumberOfClusters()+40.);
5585 if (normdist < mindist){
5591 if (index==-1) return 0;
5594 param0[index].Reset(kTRUE);
5595 FollowProlongation(param0[index],0);
5597 mother = param0[index];
5598 daughter = param1[index]; // daughter in vertex
5600 kink.SetMother(mother);
5601 kink.SetDaughter(daughter);
5603 kink.SetTPCRow0(GetRowNumber(kink.GetR()));
5604 kink.SetTPCncls(param0[index].GetNumberOfClusters(),0);
5605 kink.SetTPCncls(param1[index].GetNumberOfClusters(),1);
5606 kink.SetLabel(CookLabel(&mother,0.4, 0,kink.GetTPCRow0()),0);
5607 kink.SetLabel(CookLabel(&daughter,0.4, kink.GetTPCRow0(),160),1);
5608 mother.SetLabel(kink.GetLabel(0));
5609 daughter.SetLabel(kink.GetLabel(1));
5615 void AliTPCtrackerMI::UpdateKinkQualityM(AliTPCseed * seed){
5617 // update Kink quality information for mother after back propagation
5619 if (seed->GetKinkIndex(0)>=0) return;
5620 for (Int_t ikink=0;ikink<3;ikink++){
5621 Int_t index = seed->GetKinkIndex(ikink);
5622 if (index>=0) break;
5623 index = TMath::Abs(index)-1;
5624 AliESDkink * kink = fEvent->GetKink(index);
5625 kink->SetTPCDensity(-1,0,0);
5626 kink->SetTPCDensity(1,0,1);
5628 Int_t row0 = kink->GetTPCRow0() - 2 - Int_t( 0.5/ (0.05+kink->GetAngle(2)));
5629 if (row0<15) row0=15;
5631 Int_t row1 = kink->GetTPCRow0() + 2 + Int_t( 0.5/ (0.05+kink->GetAngle(2)));
5632 if (row1>145) row1=145;
5634 Int_t found,foundable,shared;
5635 seed->GetClusterStatistic(0,row0, found, foundable,shared,kFALSE);
5636 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),0,0);
5637 seed->GetClusterStatistic(row1,155, found, foundable,shared,kFALSE);
5638 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),0,1);
5643 void AliTPCtrackerMI::UpdateKinkQualityD(AliTPCseed * seed){
5645 // update Kink quality information for daughter after refit
5647 if (seed->GetKinkIndex(0)<=0) return;
5648 for (Int_t ikink=0;ikink<3;ikink++){
5649 Int_t index = seed->GetKinkIndex(ikink);
5650 if (index<=0) break;
5651 index = TMath::Abs(index)-1;
5652 AliESDkink * kink = fEvent->GetKink(index);
5653 kink->SetTPCDensity(-1,1,0);
5654 kink->SetTPCDensity(-1,1,1);
5656 Int_t row0 = kink->GetTPCRow0() -2 - Int_t( 0.5/ (0.05+kink->GetAngle(2)));
5657 if (row0<15) row0=15;
5659 Int_t row1 = kink->GetTPCRow0() +2 + Int_t( 0.5/ (0.05+kink->GetAngle(2)));
5660 if (row1>145) row1=145;
5662 Int_t found,foundable,shared;
5663 seed->GetClusterStatistic(0,row0, found, foundable,shared,kFALSE);
5664 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),1,0);
5665 seed->GetClusterStatistic(row1,155, found, foundable,shared,kFALSE);
5666 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),1,1);
5672 Int_t AliTPCtrackerMI::CheckKinkPoint(AliTPCseed*seed,AliTPCseed &mother, AliTPCseed &daughter, AliESDkink &knk)
5675 // check kink point for given track
5676 // if return value=0 kink point not found
5677 // otherwise seed0 correspond to mother particle
5678 // seed1 correspond to daughter particle
5679 // kink parameter of kink point
5680 AliKink &kink=(AliKink &)knk;
5682 Int_t middlerow = (seed->GetFirstPoint()+seed->GetLastPoint())/2;
5683 Int_t first = seed->GetFirstPoint();
5684 Int_t last = seed->GetLastPoint();
5685 if (last-first<20) return 0; // shortest length - 2*30 = 60 pad-rows
5688 AliTPCseed *seed1 = ReSeed(seed,middlerow+20, kTRUE); //middle of chamber
5689 if (!seed1) return 0;
5690 FollowProlongation(*seed1,seed->GetLastPoint()-20);
5691 seed1->Reset(kTRUE);
5692 FollowProlongation(*seed1,158);
5693 seed1->Reset(kTRUE);
5694 last = seed1->GetLastPoint();
5696 AliTPCseed *seed0 = new AliTPCseed(*seed);
5697 seed0->Reset(kFALSE);
5700 AliTPCseed param0[20]; // parameters along the track
5701 AliTPCseed param1[20]; // parameters along the track
5702 AliKink kinks[20]; // corresponding kink parameters
5704 for (Int_t irow=0; irow<20;irow++){
5705 rows[irow] = first +((last-first)*irow)/19;
5707 // store parameters along the track
5709 for (Int_t irow=0;irow<20;irow++){
5710 FollowBackProlongation(*seed0, rows[irow]);
5711 FollowProlongation(*seed1,rows[19-irow]);
5712 param0[irow] = *seed0;
5713 param1[19-irow] = *seed1;
5717 for (Int_t irow=0; irow<19;irow++){
5718 kinks[irow].SetMother(param0[irow]);
5719 kinks[irow].SetDaughter(param1[irow]);
5720 kinks[irow].Update();
5723 // choose kink with biggest change of angle
5725 Double_t maxchange= 0;
5726 for (Int_t irow=1;irow<19;irow++){
5727 if (TMath::Abs(kinks[irow].GetR())>240.) continue;
5728 if (TMath::Abs(kinks[irow].GetR())<110.) continue;
5729 Float_t quality = TMath::Abs(kinks[irow].GetAngle(2))/(3.+TMath::Abs(kinks[irow].GetR()-param0[irow].GetX()));
5730 if ( quality > maxchange){
5731 maxchange = quality;
5738 if (index<0) return 0;
5740 Int_t row0 = GetRowNumber(kinks[index].GetR()); //row 0 estimate
5741 seed0 = new AliTPCseed(param0[index]);
5742 seed1 = new AliTPCseed(param1[index]);
5743 seed0->Reset(kFALSE);
5744 seed1->Reset(kFALSE);
5745 seed0->ResetCovariance(10.);
5746 seed1->ResetCovariance(10.);
5747 FollowProlongation(*seed0,0);
5748 FollowBackProlongation(*seed1,158);
5749 mother = *seed0; // backup mother at position 0
5750 seed0->Reset(kFALSE);
5751 seed1->Reset(kFALSE);
5752 seed0->ResetCovariance(10.);
5753 seed1->ResetCovariance(10.);
5755 first = TMath::Max(row0-20,0);
5756 last = TMath::Min(row0+20,158);
5758 for (Int_t irow=0; irow<20;irow++){
5759 rows[irow] = first +((last-first)*irow)/19;
5761 // store parameters along the track
5763 for (Int_t irow=0;irow<20;irow++){
5764 FollowBackProlongation(*seed0, rows[irow]);
5765 FollowProlongation(*seed1,rows[19-irow]);
5766 param0[irow] = *seed0;
5767 param1[19-irow] = *seed1;
5771 for (Int_t irow=0; irow<19;irow++){
5772 kinks[irow].SetMother(param0[irow]);
5773 kinks[irow].SetDaughter(param1[irow]);
5774 // param0[irow].Dump();
5775 //param1[irow].Dump();
5776 kinks[irow].Update();
5779 // choose kink with biggest change of angle
5782 for (Int_t irow=0;irow<20;irow++){
5783 if (TMath::Abs(kinks[irow].GetR())>250.) continue;
5784 if (TMath::Abs(kinks[irow].GetR())<90.) continue;
5785 Float_t quality = TMath::Abs(kinks[irow].GetAngle(2))/(3.+TMath::Abs(kinks[irow].GetR()-param0[irow].GetX()));
5786 if ( quality > maxchange){
5787 maxchange = quality;
5794 if (index==-1 || param0[index].GetNumberOfClusters()+param1[index].GetNumberOfClusters()<100){
5799 // Float_t anglesigma = TMath::Sqrt(param0[index].fC22+param0[index].fC33+param1[index].fC22+param1[index].fC33);
5801 kink.SetMother(param0[index]);
5802 kink.SetDaughter(param1[index]);
5804 row0 = GetRowNumber(kink.GetR());
5805 kink.SetTPCRow0(row0);
5806 kink.SetLabel(CookLabel(seed0,0.5,0,row0),0);
5807 kink.SetLabel(CookLabel(seed1,0.5,row0,158),1);
5808 kink.SetIndex(-10,0);
5809 kink.SetIndex(int(param0[index].GetNumberOfClusters()+param1[index].GetNumberOfClusters()),1);
5810 kink.SetTPCncls(param0[index].GetNumberOfClusters(),0);
5811 kink.SetTPCncls(param1[index].GetNumberOfClusters(),1);
5814 // new (&mother) AliTPCseed(param0[index]);
5815 daughter = param1[index];
5816 daughter.SetLabel(kink.GetLabel(1));
5817 param0[index].Reset(kTRUE);
5818 FollowProlongation(param0[index],0);
5819 mother = param0[index];
5820 mother.SetLabel(kink.GetLabel(0));
5830 AliTPCseed* AliTPCtrackerMI::ReSeed(AliTPCseed *t)
5833 // reseed - refit - track
5836 // Int_t last = fSectors->GetNRows()-1;
5838 if (fSectors == fOuterSec){
5839 first = TMath::Max(first, t->GetFirstPoint()-fInnerSec->GetNRows());
5843 first = t->GetFirstPoint();
5845 AliTPCseed * seed = MakeSeed(t,0.1,0.5,0.9);
5846 FollowBackProlongation(*t,fSectors->GetNRows()-1);
5848 FollowProlongation(*t,first);
5858 //_____________________________________________________________________________
5859 Int_t AliTPCtrackerMI::ReadSeeds(const TFile *inp) {
5860 //-----------------------------------------------------------------
5861 // This function reades track seeds.
5862 //-----------------------------------------------------------------
5863 TDirectory *savedir=gDirectory;
5865 TFile *in=(TFile*)inp;
5866 if (!in->IsOpen()) {
5867 cerr<<"AliTPCtrackerMI::ReadSeeds(): input file is not open !\n";
5872 TTree *seedTree=(TTree*)in->Get("Seeds");
5874 cerr<<"AliTPCtrackerMI::ReadSeeds(): ";
5875 cerr<<"can't get a tree with track seeds !\n";
5878 AliTPCtrack *seed=new AliTPCtrack;
5879 seedTree->SetBranchAddress("tracks",&seed);
5881 if (fSeeds==0) fSeeds=new TObjArray(15000);
5883 Int_t n=(Int_t)seedTree->GetEntries();
5884 for (Int_t i=0; i<n; i++) {
5885 seedTree->GetEvent(i);
5886 fSeeds->AddLast(new AliTPCseed(*seed/*,seed->GetAlpha()*/));
5895 Int_t AliTPCtrackerMI::Clusters2Tracks (AliESDEvent *esd)
5898 if (fSeeds) DeleteSeeds();
5901 if (!fSeeds) return 1;
5908 //_____________________________________________________________________________
5909 Int_t AliTPCtrackerMI::Clusters2Tracks() {
5910 //-----------------------------------------------------------------
5911 // This is a track finder.
5912 //-----------------------------------------------------------------
5913 TDirectory *savedir=gDirectory;
5917 fSeeds = Tracking();
5920 Info("Clusters2Tracks","Time for tracking: \t");timer.Print();timer.Start();
5922 //activate again some tracks
5923 for (Int_t i=0; i<fSeeds->GetEntriesFast(); i++) {
5924 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
5926 Int_t nc=t.GetNumberOfClusters();
5928 delete fSeeds->RemoveAt(i);
5932 if (pt->GetRemoval()==10) {
5933 if (pt->GetDensityFirst(20)>0.8 || pt->GetDensityFirst(30)>0.8 || pt->GetDensityFirst(40)>0.7)
5934 pt->Desactivate(10); // make track again active
5936 pt->Desactivate(20);
5937 delete fSeeds->RemoveAt(i);
5942 RemoveUsed2(fSeeds,0.85,0.85,0);
5943 if (AliTPCReconstructor::GetRecoParam()->GetDoKinks()) FindKinks(fSeeds,fEvent);
5944 //FindCurling(fSeeds, fEvent,0);
5945 if (AliTPCReconstructor::StreamLevel()>2) FindMultiMC(fSeeds, fEvent,0); // find multi found tracks
5946 RemoveUsed2(fSeeds,0.5,0.4,20);
5947 FindSplitted(fSeeds, fEvent,0); // find multi found tracks
5950 // // refit short tracks
5952 Int_t nseed=fSeeds->GetEntriesFast();
5955 for (Int_t i=0; i<nseed; i++) {
5956 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
5958 Int_t nc=t.GetNumberOfClusters();
5960 delete fSeeds->RemoveAt(i);
5963 CookLabel(pt,0.1); //For comparison only
5964 //if ((pt->IsActive() || (pt->fRemoval==10) )&& nc>50 &&pt->GetNumberOfClusters()>0.4*pt->fNFoundable){
5965 if ((pt->IsActive() || (pt->GetRemoval()==10) )){
5967 if (fDebug>0) cerr<<found<<'\r';
5971 delete fSeeds->RemoveAt(i);
5975 //RemoveOverlap(fSeeds,0.99,7,kTRUE);
5977 //RemoveUsed(fSeeds,0.9,0.9,6);
5979 nseed=fSeeds->GetEntriesFast();
5981 for (Int_t i=0; i<nseed; i++) {
5982 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
5984 Int_t nc=t.GetNumberOfClusters();
5986 delete fSeeds->RemoveAt(i);
5990 t.CookdEdx(0.02,0.6);
5991 // CheckKinkPoint(&t,0.05);
5992 //if ((pt->IsActive() || (pt->fRemoval==10) )&& nc>50 &&pt->GetNumberOfClusters()>0.4*pt->fNFoundable){
5993 if ((pt->IsActive() || (pt->GetRemoval()==10) )){
6001 delete fSeeds->RemoveAt(i);
6002 //AliTPCseed * seed1 = ReSeed(pt,0.05,0.5,1);
6004 // FollowProlongation(*seed1,0);
6005 // Int_t n = seed1->GetNumberOfClusters();
6006 // printf("fP4\t%f\t%f\n",seed1->GetC(),pt->GetC());
6007 // printf("fN\t%d\t%d\n", seed1->GetNumberOfClusters(),pt->GetNumberOfClusters());
6010 //AliTPCseed * seed2 = ReSeed(pt,0.95,0.5,0.05);
6014 SortTracks(fSeeds, 1);
6018 PrepareForBackProlongation(fSeeds,5.);
6019 PropagateBack(fSeeds);
6020 printf("Time for back propagation: \t");timer.Print();timer.Start();
6024 PrepareForProlongation(fSeeds,5.);
6025 PropagateForward2(fSeeds);
6027 printf("Time for FORWARD propagation: \t");timer.Print();timer.Start();
6028 // RemoveUsed(fSeeds,0.7,0.7,6);
6029 //RemoveOverlap(fSeeds,0.9,7,kTRUE);
6031 nseed=fSeeds->GetEntriesFast();
6033 for (Int_t i=0; i<nseed; i++) {
6034 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
6036 Int_t nc=t.GetNumberOfClusters();
6038 delete fSeeds->RemoveAt(i);
6041 t.CookdEdx(0.02,0.6);
6042 // CookLabel(pt,0.1); //For comparison only
6043 //if ((pt->IsActive() || (pt->fRemoval==10) )&& nc>50 &&pt->GetNumberOfClusters()>0.4*pt->fNFoundable){
6044 if ((pt->IsActive() || (pt->fRemoval==10) )){
6045 cerr<<found++<<'\r';
6048 delete fSeeds->RemoveAt(i);
6053 // fNTracks = found;
6055 Info("Clusters2Tracks","Time for overlap removal, track writing and dedx cooking: \t"); timer.Print();timer.Start();
6058 // cerr<<"Number of found tracks : "<<"\t"<<found<<endl;
6059 Info("Clusters2Tracks","Number of found tracks %d",found);
6061 // UnloadClusters();
6066 void AliTPCtrackerMI::Tracking(TObjArray * arr)
6069 // tracking of the seeds
6072 fSectors = fOuterSec;
6073 ParallelTracking(arr,150,63);
6074 fSectors = fOuterSec;
6075 ParallelTracking(arr,63,0);
6078 TObjArray * AliTPCtrackerMI::Tracking(Int_t seedtype, Int_t i1, Int_t i2, Float_t cuts[4], Float_t dy, Int_t dsec)
6083 TObjArray * arr = new TObjArray;
6085 fSectors = fOuterSec;
6088 for (Int_t sec=0;sec<fkNOS;sec++){
6089 if (seedtype==3) MakeSeeds3(arr,sec,i1,i2,cuts,dy, dsec);
6090 if (seedtype==4) MakeSeeds5(arr,sec,i1,i2,cuts,dy);
6091 if (seedtype==2) MakeSeeds2(arr,sec,i1,i2,cuts,dy);
6094 Info("Tracking","\nSeeding - %d\t%d\t%d\t%d\n",seedtype,i1,i2,arr->GetEntriesFast());
6106 TObjArray * AliTPCtrackerMI::Tracking()
6110 if (AliTPCReconstructor::GetRecoParam()->GetSpecialSeeding()) return TrackingSpecial();
6113 Int_t nup=fOuterSec->GetNRows()+fInnerSec->GetNRows();
6115 TObjArray * seeds = new TObjArray;
6124 Float_t fnumber = 3.0;
6125 Float_t fdensity = 3.0;
6130 for (Int_t delta = 0; delta<18; delta+=6){
6134 arr = Tracking(3,nup-1-delta,nup-1-delta-gap,cuts,-1,1);
6135 SumTracks(seeds,arr);
6136 SignClusters(seeds,fnumber,fdensity);
6138 for (Int_t i=2;i<6;i+=2){
6139 // seed high pt tracks
6142 arr = Tracking(3,nup-i-delta,nup-i-delta-gap,cuts,-1,0);
6143 SumTracks(seeds,arr);
6144 SignClusters(seeds,fnumber,fdensity);
6149 // RemoveUsed(seeds,0.9,0.9,1);
6150 // UnsignClusters();
6151 // SignClusters(seeds,fnumber,fdensity);
6155 for (Int_t delta = 20; delta<120; delta+=10){
6157 // seed high pt tracks
6161 arr = Tracking(3,nup-delta,nup-delta-gap,cuts,-1);
6162 SumTracks(seeds,arr);
6163 SignClusters(seeds,fnumber,fdensity);
6168 arr = Tracking(3,nup-delta-5,nup-delta-5-gap,cuts,-1);
6169 SumTracks(seeds,arr);
6170 SignClusters(seeds,fnumber,fdensity);
6181 Info("Tracking()","\n\nPrimary seeding\t%d\n\n",seeds->GetEntriesFast());
6185 // RemoveUsed(seeds,0.75,0.75,1);
6187 //SignClusters(seeds,fnumber,fdensity);
6196 arr = Tracking(4,nup-1,nup-1-gap,cuts,-1);
6197 SumTracks(seeds,arr);
6198 SignClusters(seeds,fnumber,fdensity);
6200 arr = Tracking(4,nup-2,nup-2-gap,cuts,-1);
6201 SumTracks(seeds,arr);
6202 SignClusters(seeds,fnumber,fdensity);
6204 arr = Tracking(4,nup-3,nup-3-gap,cuts,-1);
6205 SumTracks(seeds,arr);
6206 SignClusters(seeds,fnumber,fdensity);
6210 for (Int_t delta = 3; delta<30; delta+=5){
6216 arr = Tracking(4,nup-1-delta,nup-1-delta-gap,cuts,-1);
6217 SumTracks(seeds,arr);
6218 SignClusters(seeds,fnumber,fdensity);
6220 arr = Tracking(4,nup-3-delta,nup-5-delta-gap,cuts,4);
6221 SumTracks(seeds,arr);
6222 SignClusters(seeds,fnumber,fdensity);
6233 Int_t fLastSeedRowSec=AliTPCReconstructor::GetRecoParam()->GetLastSeedRowSec();
6236 for (Int_t delta = 30; delta<fLastSeedRowSec; delta+=10){
6242 arr = Tracking(4,nup-1-delta,nup-1-delta-gap,cuts,-1);
6243 SumTracks(seeds,arr);
6244 SignClusters(seeds,fnumber,fdensity);
6246 arr = Tracking(4,nup-5-delta,nup-5-delta-gap,cuts,5 );
6247 SumTracks(seeds,arr);
6248 SignClusters(seeds,fnumber,fdensity);
6252 Info("Tracking()","\n\nSecondary seeding\t%d\n\n",seeds->GetEntriesFast());
6263 TObjArray * AliTPCtrackerMI::TrackingSpecial()
6266 // seeding adjusted for laser and cosmic tests - short tracks with big inclination angle
6267 // no primary vertex seeding tried
6271 Int_t nup=fOuterSec->GetNRows()+fInnerSec->GetNRows();
6273 TObjArray * seeds = new TObjArray;
6278 Float_t fnumber = 3.0;
6279 Float_t fdensity = 3.0;
6282 cuts[0] = AliTPCReconstructor::GetRecoParam()->GetMaxC(); // max curvature
6283 cuts[1] = 3.5; // max tan(phi) angle for seeding
6284 cuts[2] = 3.; // not used (cut on z primary vertex)
6285 cuts[3] = 3.5; // max tan(theta) angle for seeding
6287 for (Int_t delta = 0; nup-delta-gap-1>0; delta+=3){
6289 arr = Tracking(4,nup-1-delta,nup-1-delta-gap,cuts,-1);
6290 SumTracks(seeds,arr);
6291 SignClusters(seeds,fnumber,fdensity);
6295 Info("Tracking()","\n\nSecondary seeding\t%d\n\n",seeds->GetEntriesFast());
6306 void AliTPCtrackerMI::SumTracks(TObjArray *arr1,TObjArray *&arr2) const
6309 //sum tracks to common container
6310 //remove suspicious tracks
6311 Int_t nseed = arr2->GetEntriesFast();
6312 for (Int_t i=0;i<nseed;i++){
6313 AliTPCseed *pt=(AliTPCseed*)arr2->UncheckedAt(i);
6316 // remove tracks with too big curvature
6318 if (TMath::Abs(pt->GetC())>AliTPCReconstructor::GetRecoParam()->GetMaxC()){
6319 delete arr2->RemoveAt(i);
6322 // REMOVE VERY SHORT TRACKS
6323 if (pt->GetNumberOfClusters()<20){
6324 delete arr2->RemoveAt(i);
6327 // NORMAL ACTIVE TRACK
6328 if (pt->IsActive()){
6329 arr1->AddLast(arr2->RemoveAt(i));
6332 //remove not usable tracks
6333 if (pt->GetRemoval()!=10){
6334 delete arr2->RemoveAt(i);
6338 // ENABLE ONLY ENOUGH GOOD STOPPED TRACKS
6339 if (pt->GetDensityFirst(20)>0.8 || pt->GetDensityFirst(30)>0.8 || pt->GetDensityFirst(40)>0.7)
6340 arr1->AddLast(arr2->RemoveAt(i));
6342 delete arr2->RemoveAt(i);
6346 delete arr2; arr2 = 0;
6351 void AliTPCtrackerMI::ParallelTracking(TObjArray * arr, Int_t rfirst, Int_t rlast)
6354 // try to track in parralel
6356 Int_t nseed=arr->GetEntriesFast();
6357 //prepare seeds for tracking
6358 for (Int_t i=0; i<nseed; i++) {
6359 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i), &t=*pt;
6361 if (!t.IsActive()) continue;
6362 // follow prolongation to the first layer
6363 if ( (fSectors ==fInnerSec) || (t.GetFirstPoint()-fParam->GetNRowLow()>rfirst+1) )
6364 FollowProlongation(t, rfirst+1);
6369 for (Int_t nr=rfirst; nr>=rlast; nr--){
6370 if (nr<fInnerSec->GetNRows())
6371 fSectors = fInnerSec;
6373 fSectors = fOuterSec;
6374 // make indexes with the cluster tracks for given
6376 // find nearest cluster
6377 for (Int_t i=0; i<nseed; i++) {
6378 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i), &t=*pt;
6380 if (nr==80) pt->UpdateReference();
6381 if (!pt->IsActive()) continue;
6382 // if ( (fSectors ==fOuterSec) && (pt->fFirstPoint-fParam->GetNRowLow())<nr) continue;
6383 if (pt->GetRelativeSector()>17) {
6386 UpdateClusters(t,nr);
6388 // prolonagate to the nearest cluster - if founded
6389 for (Int_t i=0; i<nseed; i++) {
6390 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
6392 if (!pt->IsActive()) continue;
6393 // if ((fSectors ==fOuterSec) && (pt->fFirstPoint-fParam->GetNRowLow())<nr) continue;
6394 if (pt->GetRelativeSector()>17) {
6397 FollowToNextCluster(*pt,nr);
6402 void AliTPCtrackerMI::PrepareForBackProlongation(TObjArray * arr,Float_t fac) const
6406 // if we use TPC track itself we have to "update" covariance
6408 Int_t nseed= arr->GetEntriesFast();
6409 for (Int_t i=0;i<nseed;i++){
6410 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
6414 //rotate to current local system at first accepted point
6415 Int_t index = pt->GetClusterIndex2(pt->GetFirstPoint());
6416 Int_t sec = (index&0xff000000)>>24;
6418 Float_t angle1 = fInnerSec->GetAlpha()*sec+fInnerSec->GetAlphaShift();
6419 if (angle1>TMath::Pi())
6420 angle1-=2.*TMath::Pi();
6421 Float_t angle2 = pt->GetAlpha();
6423 if (TMath::Abs(angle1-angle2)>0.001){
6424 pt->Rotate(angle1-angle2);
6425 //angle2 = pt->GetAlpha();
6426 //pt->fRelativeSector = pt->GetAlpha()/fInnerSec->GetAlpha();
6427 //if (pt->GetAlpha()<0)
6428 // pt->fRelativeSector+=18;
6429 //sec = pt->fRelativeSector;
6438 void AliTPCtrackerMI::PrepareForProlongation(TObjArray * arr, Float_t fac) const
6442 // if we use TPC track itself we have to "update" covariance
6444 Int_t nseed= arr->GetEntriesFast();
6445 for (Int_t i=0;i<nseed;i++){
6446 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
6449 pt->SetFirstPoint(pt->GetLastPoint());
6457 Int_t AliTPCtrackerMI::PropagateBack(TObjArray * arr)
6460 // make back propagation
6462 Int_t nseed= arr->GetEntriesFast();
6463 for (Int_t i=0;i<nseed;i++){
6464 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
6465 if (pt&& pt->GetKinkIndex(0)<=0) {
6466 //AliTPCseed *pt2 = new AliTPCseed(*pt);
6467 fSectors = fInnerSec;
6468 //FollowBackProlongation(*pt,fInnerSec->GetNRows()-1);
6469 //fSectors = fOuterSec;
6470 FollowBackProlongation(*pt,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1);
6471 //if (pt->GetNumberOfClusters()<(pt->fEsd->GetTPCclusters(0)) ){
6472 // Error("PropagateBack","Not prolonged track %d",pt->GetLabel());
6473 // FollowBackProlongation(*pt2,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1);
6476 if (pt&& pt->GetKinkIndex(0)>0) {
6477 AliESDkink * kink = fEvent->GetKink(pt->GetKinkIndex(0)-1);
6478 pt->SetFirstPoint(kink->GetTPCRow0());
6479 fSectors = fInnerSec;
6480 FollowBackProlongation(*pt,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1);
6488 Int_t AliTPCtrackerMI::PropagateForward2(TObjArray * arr)
6491 // make forward propagation
6493 Int_t nseed= arr->GetEntriesFast();
6495 for (Int_t i=0;i<nseed;i++){
6496 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
6498 FollowProlongation(*pt,0);
6507 Int_t AliTPCtrackerMI::PropagateForward()
6510 // propagate track forward
6512 Int_t nseed = fSeeds->GetEntriesFast();
6513 for (Int_t i=0;i<nseed;i++){
6514 AliTPCseed *pt = (AliTPCseed*)fSeeds->UncheckedAt(i);
6516 AliTPCseed &t = *pt;
6517 Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
6518 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
6519 if (alpha < 0. ) alpha += 2.*TMath::Pi();
6520 t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN);
6524 fSectors = fOuterSec;
6525 ParallelTracking(fSeeds,fOuterSec->GetNRows()+fInnerSec->GetNRows()-1,fInnerSec->GetNRows());
6526 fSectors = fInnerSec;
6527 ParallelTracking(fSeeds,fInnerSec->GetNRows()-1,0);
6536 Int_t AliTPCtrackerMI::PropagateBack(AliTPCseed * pt, Int_t row0, Int_t row1)
6539 // make back propagation, in between row0 and row1
6543 fSectors = fInnerSec;
6546 if (row1<fSectors->GetNRows())
6549 r1 = fSectors->GetNRows()-1;
6551 if (row0<fSectors->GetNRows()&& r1>0 )
6552 FollowBackProlongation(*pt,r1);
6553 if (row1<=fSectors->GetNRows())
6556 r1 = row1 - fSectors->GetNRows();
6557 if (r1<=0) return 0;
6558 if (r1>=fOuterSec->GetNRows()) return 0;
6559 fSectors = fOuterSec;
6560 return FollowBackProlongation(*pt,r1);
6568 void AliTPCtrackerMI::GetShape(AliTPCseed * seed, Int_t row)
6572 AliTPCClusterParam * clparam = AliTPCcalibDB::Instance()->GetClusterParam();
6573 Float_t zdrift = TMath::Abs((fParam->GetZLength(0)-TMath::Abs(seed->GetZ())));
6574 Int_t type = (seed->GetSector() < fParam->GetNSector()/2) ? 0: (row>126) ? 1:2;
6575 Double_t angulary = seed->GetSnp();
6576 angulary = angulary*angulary/(1.-angulary*angulary);
6577 Double_t angularz = seed->GetTgl()*seed->GetTgl()*(1.+angulary);
6579 Double_t sigmay = clparam->GetRMS0(0,type,zdrift,TMath::Sqrt(TMath::Abs(angulary)));
6580 Double_t sigmaz = clparam->GetRMS0(1,type,zdrift,TMath::Sqrt(TMath::Abs(angularz)));
6581 seed->SetCurrentSigmaY2(sigmay*sigmay);
6582 seed->SetCurrentSigmaZ2(sigmaz*sigmaz);
6583 // Float_t sd2 = TMath::Abs((fParam->GetZLength(0)-TMath::Abs(seed->GetZ())))*fParam->GetDiffL()*fParam->GetDiffL();
6584 // // Float_t padlength = fParam->GetPadPitchLength(seed->fSector);
6585 // Float_t padlength = GetPadPitchLength(row);
6587 // Float_t sresy = (seed->GetSector() < fParam->GetNSector()/2) ? 0.2 :0.3;
6588 // seed->SetCurrentSigmaY2(sd2+padlength*padlength*angulary/12.+sresy*sresy);
6590 // Float_t sresz = fParam->GetZSigma();
6591 // seed->SetCurrentSigmaZ2(sd2+padlength*padlength*angularz*angularz*(1+angulary)/12.+sresz*sresz);
6593 Float_t wy = GetSigmaY(seed);
6594 Float_t wz = GetSigmaZ(seed);
6597 if (TMath::Abs(wy/seed->fCurrentSigmaY2-1)>0.0001 || TMath::Abs(wz/seed->fCurrentSigmaZ2-1)>0.0001 ){
6598 printf("problem\n");
6605 //__________________________________________________________________________
6606 void AliTPCtrackerMI::CookLabel(AliKalmanTrack *tk, Float_t wrong) const {
6607 //--------------------------------------------------------------------
6608 //This function "cooks" a track label. If label<0, this track is fake.
6609 //--------------------------------------------------------------------
6610 AliTPCseed * t = dynamic_cast<AliTPCseed*>(tk);
6612 printf("%s:%d wrong type \n",(char*)__FILE__,__LINE__);
6616 Int_t noc=t->GetNumberOfClusters();
6618 //printf("\nnot founded prolongation\n\n\n");
6624 AliTPCclusterMI *clusters[160];
6626 for (Int_t i=0;i<160;i++) {
6633 for (i=0; i<160 && current<noc; i++) {
6635 Int_t index=t->GetClusterIndex2(i);
6636 if (index<=0) continue;
6637 if (index&0x8000) continue;
6639 //clusters[current]=GetClusterMI(index);
6640 if (t->GetClusterPointer(i)){
6641 clusters[current]=t->GetClusterPointer(i);
6647 Int_t lab=123456789;
6648 for (i=0; i<noc; i++) {
6649 AliTPCclusterMI *c=clusters[i];
6651 lab=TMath::Abs(c->GetLabel(0));
6653 for (j=0; j<noc; j++) if (lb[j]==lab || mx[j]==0) break;
6659 for (i=0; i<noc; i++) if (mx[i]>max) {max=mx[i]; lab=lb[i];}
6661 for (i=0; i<noc; i++) {
6662 AliTPCclusterMI *c=clusters[i];
6664 if (TMath::Abs(c->GetLabel(1)) == lab ||
6665 TMath::Abs(c->GetLabel(2)) == lab ) max++;
6668 if ((1.- Float_t(max)/noc) > wrong) lab=-lab;
6671 Int_t tail=Int_t(0.10*noc);
6674 for (i=1; i<=160&&ind<tail; i++) {
6675 // AliTPCclusterMI *c=clusters[noc-i];
6676 AliTPCclusterMI *c=clusters[i];
6678 if (lab == TMath::Abs(c->GetLabel(0)) ||
6679 lab == TMath::Abs(c->GetLabel(1)) ||
6680 lab == TMath::Abs(c->GetLabel(2))) max++;
6683 if (max < Int_t(0.5*tail)) lab=-lab;
6690 //delete[] clusters;
6694 //__________________________________________________________________________
6695 Int_t AliTPCtrackerMI::CookLabel(AliTPCseed *t, Float_t wrong,Int_t first, Int_t last) const {
6696 //--------------------------------------------------------------------
6697 //This function "cooks" a track label. If label<0, this track is fake.
6698 //--------------------------------------------------------------------
6699 Int_t noc=t->GetNumberOfClusters();
6701 //printf("\nnot founded prolongation\n\n\n");
6707 AliTPCclusterMI *clusters[160];
6709 for (Int_t i=0;i<160;i++) {
6716 for (i=0; i<160 && current<noc; i++) {
6717 if (i<first) continue;
6718 if (i>last) continue;
6719 Int_t index=t->GetClusterIndex2(i);
6720 if (index<=0) continue;
6721 if (index&0x8000) continue;
6723 //clusters[current]=GetClusterMI(index);
6724 if (t->GetClusterPointer(i)){
6725 clusters[current]=t->GetClusterPointer(i);
6730 if (noc<5) return -1;
6731 Int_t lab=123456789;
6732 for (i=0; i<noc; i++) {
6733 AliTPCclusterMI *c=clusters[i];
6735 lab=TMath::Abs(c->GetLabel(0));
6737 for (j=0; j<noc; j++) if (lb[j]==lab || mx[j]==0) break;
6743 for (i=0; i<noc; i++) if (mx[i]>max) {max=mx[i]; lab=lb[i];}
6745 for (i=0; i<noc; i++) {
6746 AliTPCclusterMI *c=clusters[i];
6748 if (TMath::Abs(c->GetLabel(1)) == lab ||
6749 TMath::Abs(c->GetLabel(2)) == lab ) max++;
6752 if ((1.- Float_t(max)/noc) > wrong) lab=-lab;
6755 Int_t tail=Int_t(0.10*noc);
6758 for (i=1; i<=160&&ind<tail; i++) {
6759 // AliTPCclusterMI *c=clusters[noc-i];
6760 AliTPCclusterMI *c=clusters[i];
6762 if (lab == TMath::Abs(c->GetLabel(0)) ||
6763 lab == TMath::Abs(c->GetLabel(1)) ||
6764 lab == TMath::Abs(c->GetLabel(2))) max++;
6767 if (max < Int_t(0.5*tail)) lab=-lab;
6770 // t->SetLabel(lab);
6774 //delete[] clusters;
6778 Int_t AliTPCtrackerMI::AliTPCSector::GetRowNumber(Double_t x) const
6780 //return pad row number for this x
6783 r=fRow[fN-1].GetX();
6784 if (x > r) return fN;
6786 if (x < r) return -1;
6787 return Int_t((x-r)/fPadPitchLength + 0.5);}
6789 r=fRow[fN-1].GetX();
6790 if (x > r) return fN;
6792 if (x < r) return -1;
6793 Double_t r1=fRow[64].GetX();
6795 return Int_t((x-r)/f1PadPitchLength + 0.5);}
6797 return (Int_t((x-r1)/f2PadPitchLength + 0.5)+64);}
6801 Int_t AliTPCtrackerMI::GetRowNumber(Double_t x[3]) const
6803 //return pad row number for given x vector
6804 Float_t phi = TMath::ATan2(x[1],x[0]);
6805 if(phi<0) phi=2.*TMath::Pi()+phi;
6806 // Get the local angle in the sector philoc
6807 const Float_t kRaddeg = 180/3.14159265358979312;
6808 Float_t phiangle = (Int_t (phi*kRaddeg/20.) + 0.5)*20./kRaddeg;
6809 Double_t localx = x[0]*TMath::Cos(phiangle)-x[1]*TMath::Sin(phiangle);
6810 return GetRowNumber(localx);
6813 //_________________________________________________________________________
6814 void AliTPCtrackerMI::AliTPCSector::Setup(const AliTPCParam *par, Int_t f) {
6815 //-----------------------------------------------------------------------
6816 // Setup inner sector
6817 //-----------------------------------------------------------------------
6819 fAlpha=par->GetInnerAngle();
6820 fAlphaShift=par->GetInnerAngleShift();
6821 fPadPitchWidth=par->GetInnerPadPitchWidth();
6822 fPadPitchLength=par->GetInnerPadPitchLength();
6823 fN=par->GetNRowLow();
6824 if(fRow)delete [] fRow;fRow = 0;
6825 fRow=new AliTPCRow[fN];
6826 for (Int_t i=0; i<fN; i++) {
6827 fRow[i].SetX(par->GetPadRowRadiiLow(i));
6828 fRow[i].SetDeadZone(1.5); //1.5 cm of dead zone
6831 fAlpha=par->GetOuterAngle();
6832 fAlphaShift=par->GetOuterAngleShift();
6833 fPadPitchWidth = par->GetOuterPadPitchWidth();
6834 fPadPitchLength = par->GetOuter1PadPitchLength();
6835 f1PadPitchLength = par->GetOuter1PadPitchLength();
6836 f2PadPitchLength = par->GetOuter2PadPitchLength();
6837 fN=par->GetNRowUp();
6838 if(fRow)delete [] fRow;fRow = 0;
6839 fRow=new AliTPCRow[fN];
6840 for (Int_t i=0; i<fN; i++) {
6841 fRow[i].SetX(par->GetPadRowRadiiUp(i));
6842 fRow[i].SetDeadZone(1.5); // 1.5 cm of dead zone
6847 AliTPCtrackerMI::AliTPCRow::AliTPCRow():
6857 // default constructor
6861 AliTPCtrackerMI::AliTPCRow::~AliTPCRow(){
6863 for (Int_t i = 0; i < fN1; i++)
6864 fClusters1[i].~AliTPCclusterMI();
6865 delete [] fClusters1;
6866 for (Int_t i = 0; i < fN2; i++)
6867 fClusters2[i].~AliTPCclusterMI();
6868 delete [] fClusters2;
6873 //_________________________________________________________________________
6875 AliTPCtrackerMI::AliTPCRow::InsertCluster(const AliTPCclusterMI* c, UInt_t index) {
6876 //-----------------------------------------------------------------------
6877 // Insert a cluster into this pad row in accordence with its y-coordinate
6878 //-----------------------------------------------------------------------
6879 if (fN==kMaxClusterPerRow) {
6880 //AliInfo("AliTPCRow::InsertCluster(): Too many clusters");
6884 //AliInfo("AliTPCRow::InsertCluster(): Too many clusters !");
6887 if (fN==0) {fIndex[0]=index; fClusters[fN++]=c; return;}
6888 Int_t i=Find(c->GetZ());
6889 memmove(fClusters+i+1 ,fClusters+i,(fN-i)*sizeof(AliTPCclusterMI*));
6890 memmove(fIndex +i+1 ,fIndex +i,(fN-i)*sizeof(UInt_t));
6891 fIndex[i]=index; fClusters[i]=c; fN++;
6894 void AliTPCtrackerMI::AliTPCRow::ResetClusters() {
6897 // MvL: Need to call destructors for AliTPCclusterMI, to delete fInfo
6898 for (Int_t i = 0; i < fN1; i++)
6899 fClusters1[i].~AliTPCclusterMI();
6900 delete [] fClusters1; fClusters1=0;
6901 for (Int_t i = 0; i < fN2; i++)
6902 fClusters2[i].~AliTPCclusterMI();
6903 delete [] fClusters2; fClusters2=0;
6908 //delete[] fClusterArray;
6914 //___________________________________________________________________
6915 Int_t AliTPCtrackerMI::AliTPCRow::Find(Double_t z) const {
6916 //-----------------------------------------------------------------------
6917 // Return the index of the nearest cluster
6918 //-----------------------------------------------------------------------
6919 if (fN==0) return 0;
6920 if (z <= fClusters[0]->GetZ()) return 0;
6921 if (z > fClusters[fN-1]->GetZ()) return fN;
6922 Int_t b=0, e=fN-1, m=(b+e)/2;
6923 for (; b<e; m=(b+e)/2) {
6924 if (z > fClusters[m]->GetZ()) b=m+1;
6932 //___________________________________________________________________
6933 AliTPCclusterMI * AliTPCtrackerMI::AliTPCRow::FindNearest(Double_t y, Double_t z, Double_t roady, Double_t roadz) const {
6934 //-----------------------------------------------------------------------
6935 // Return the index of the nearest cluster in z y
6936 //-----------------------------------------------------------------------
6937 Float_t maxdistance = roady*roady + roadz*roadz;
6939 AliTPCclusterMI *cl =0;
6940 for (Int_t i=Find(z-roadz); i<fN; i++) {
6941 AliTPCclusterMI *c=(AliTPCclusterMI*)(fClusters[i]);
6942 if (c->GetZ() > z+roadz) break;
6943 if ( (c->GetY()-y) > roady ) continue;
6944 Float_t distance = (c->GetZ()-z)*(c->GetZ()-z)+(c->GetY()-y)*(c->GetY()-y);
6945 if (maxdistance>distance) {
6946 maxdistance = distance;
6953 AliTPCclusterMI * AliTPCtrackerMI::AliTPCRow::FindNearest2(Double_t y, Double_t z, Double_t roady, Double_t roadz,UInt_t & index) const
6955 //-----------------------------------------------------------------------
6956 // Return the index of the nearest cluster in z y
6957 //-----------------------------------------------------------------------
6958 Float_t maxdistance = roady*roady + roadz*roadz;
6959 AliTPCclusterMI *cl =0;
6961 //PH Check boundaries. 510 is the size of fFastCluster
6962 Int_t iz1 = Int_t(z-roadz+254.5);
6963 if (iz1<0 || iz1>=510) return cl;
6964 iz1 = TMath::Max(GetFastCluster(iz1)-1,0);
6965 Int_t iz2 = Int_t(z+roadz+255.5);
6966 if (iz2<0 || iz2>=510) return cl;
6967 iz2 = TMath::Min(GetFastCluster(iz2)+1,fN);
6968 Bool_t skipUsed = !(AliTPCReconstructor::GetRecoParam()->GetClusterSharing());
6969 //FindNearest3(y,z,roady,roadz,index);
6970 // for (Int_t i=Find(z-roadz); i<fN; i++) {
6971 for (Int_t i=iz1; i<iz2; i++) {
6972 AliTPCclusterMI *c=(AliTPCclusterMI*)(fClusters[i]);
6973 if (c->GetZ() > z+roadz) break;
6974 if ( c->GetY()-y > roady ) continue;
6975 if ( y-c->GetY() > roady ) continue;
6976 if (skipUsed && c->IsUsed(11)) continue;
6977 Float_t distance = (c->GetZ()-z)*(c->GetZ()-z)+(c->GetY()-y)*(c->GetY()-y);
6978 if (maxdistance>distance) {
6979 maxdistance = distance;
6982 //roady = TMath::Sqrt(maxdistance);
6989 void AliTPCtrackerMI::AliTPCRow::SetFastCluster(Int_t i, Short_t cl){
6991 // Set cluster info for fast navigation
7002 void AliTPCtrackerMI::MakeBitmaps(AliTPCseed *t)
7004 //-----------------------------------------------------------------------
7005 // Fill the cluster and sharing bitmaps of the track
7006 //-----------------------------------------------------------------------
7008 Int_t firstpoint = 0;
7009 Int_t lastpoint = 159;
7010 AliTPCTrackerPoint *point;
7012 for (int iter=firstpoint; iter<lastpoint; iter++) {
7013 point = t->GetTrackPoint(iter);
7015 t->SetClusterMapBit(iter, kTRUE);
7016 if (point->IsShared())
7017 t->SetSharedMapBit(iter,kTRUE);
7019 t->SetSharedMapBit(iter, kFALSE);
7022 t->SetClusterMapBit(iter, kFALSE);
7023 t->SetSharedMapBit(iter, kFALSE);