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)
124 ClassImp(AliTPCtrackerRow)
127 class AliTPCFastMath {
130 static Double_t FastAsin(Double_t x);
132 static Double_t fgFastAsin[20000]; //lookup table for fast asin computation
135 Double_t AliTPCFastMath::fgFastAsin[20000];
136 AliTPCFastMath gAliTPCFastMath; // needed to fill the LUT
138 AliTPCFastMath::AliTPCFastMath(){
140 // initialized lookup table;
141 for (Int_t i=0;i<10000;i++){
142 fgFastAsin[2*i] = TMath::ASin(i/10000.);
143 fgFastAsin[2*i+1] = (TMath::ASin((i+1)/10000.)-fgFastAsin[2*i]);
147 Double_t AliTPCFastMath::FastAsin(Double_t x){
149 // return asin using lookup table
151 Int_t index = int(x*10000);
152 return fgFastAsin[2*index]+(x*10000.-index)*fgFastAsin[2*index+1];
155 Int_t index = int(x*10000);
156 return -(fgFastAsin[2*index]+(x*10000.-index)*fgFastAsin[2*index+1]);
158 //__________________________________________________________________
159 AliTPCtrackerMI::AliTPCtrackerMI()
181 // default constructor
184 //_____________________________________________________________________
188 Int_t AliTPCtrackerMI::UpdateTrack(AliTPCseed * track, Int_t accept){
190 //update track information using current cluster - track->fCurrentCluster
193 AliTPCclusterMI* c =track->GetCurrentCluster();
194 if (accept>0) track->SetCurrentClusterIndex1(track->GetCurrentClusterIndex1() | 0x8000); //sign not accepted clusters
196 UInt_t i = track->GetCurrentClusterIndex1();
198 Int_t sec=(i&0xff000000)>>24;
199 //Int_t row = (i&0x00ff0000)>>16;
200 track->SetRow((i&0x00ff0000)>>16);
201 track->SetSector(sec);
202 // Int_t index = i&0xFFFF;
203 if (sec>=fParam->GetNInnerSector()) track->SetRow(track->GetRow()+fParam->GetNRowLow());
204 track->SetClusterIndex2(track->GetRow(), i);
205 //track->fFirstPoint = row;
206 //if ( track->fLastPoint<row) track->fLastPoint =row;
207 // if (track->fRow<0 || track->fRow>160) {
208 // printf("problem\n");
210 if (track->GetFirstPoint()>track->GetRow())
211 track->SetFirstPoint(track->GetRow());
212 if (track->GetLastPoint()<track->GetRow())
213 track->SetLastPoint(track->GetRow());
216 track->SetClusterPointer(track->GetRow(),c);
219 Double_t angle2 = track->GetSnp()*track->GetSnp();
221 //SET NEW Track Point
223 if (angle2<1) //PH sometimes angle2 is very big. To be investigated...
225 angle2 = TMath::Sqrt(angle2/(1-angle2));
226 AliTPCTrackerPoint &point =*(track->GetTrackPoint(track->GetRow()));
228 point.SetSigmaY(c->GetSigmaY2()/track->GetCurrentSigmaY2());
229 point.SetSigmaZ(c->GetSigmaZ2()/track->GetCurrentSigmaZ2());
230 point.SetErrY(sqrt(track->GetErrorY2()));
231 point.SetErrZ(sqrt(track->GetErrorZ2()));
233 point.SetX(track->GetX());
234 point.SetY(track->GetY());
235 point.SetZ(track->GetZ());
236 point.SetAngleY(angle2);
237 point.SetAngleZ(track->GetTgl());
238 if (point.IsShared()){
239 track->SetErrorY2(track->GetErrorY2()*4);
240 track->SetErrorZ2(track->GetErrorZ2()*4);
244 Double_t chi2 = track->GetPredictedChi2(track->GetCurrentCluster());
246 // track->SetErrorY2(track->GetErrorY2()*1.3);
247 // track->SetErrorY2(track->GetErrorY2()+0.01);
248 // track->SetErrorZ2(track->GetErrorZ2()*1.3);
249 // track->SetErrorZ2(track->GetErrorZ2()+0.005);
251 if (accept>0) return 0;
252 if (track->GetNumberOfClusters()%20==0){
253 // if (track->fHelixIn){
254 // TClonesArray & larr = *(track->fHelixIn);
255 // Int_t ihelix = larr.GetEntriesFast();
256 // new(larr[ihelix]) AliHelix(*track) ;
259 track->SetNoCluster(0);
260 return track->Update(c,chi2,i);
265 Int_t AliTPCtrackerMI::AcceptCluster(AliTPCseed * seed, AliTPCclusterMI * cluster)
268 // decide according desired precision to accept given
269 // cluster for tracking
270 Double_t sy2=ErrY2(seed,cluster);
271 Double_t sz2=ErrZ2(seed,cluster);
273 Double_t sdistancey2 = sy2+seed->GetSigmaY2();
274 Double_t sdistancez2 = sz2+seed->GetSigmaZ2();
276 Double_t rdistancey2 = (seed->GetCurrentCluster()->GetY()-seed->GetY())*
277 (seed->GetCurrentCluster()->GetY()-seed->GetY())/sdistancey2;
278 Double_t rdistancez2 = (seed->GetCurrentCluster()->GetZ()-seed->GetZ())*
279 (seed->GetCurrentCluster()->GetZ()-seed->GetZ())/sdistancez2;
281 Double_t rdistance2 = rdistancey2+rdistancez2;
284 if (AliTPCReconstructor::StreamLevel()>5) {
285 Float_t rmsy2 = seed->GetCurrentSigmaY2();
286 Float_t rmsz2 = seed->GetCurrentSigmaZ2();
287 Float_t rmsy2p30 = seed->GetCMeanSigmaY2p30();
288 Float_t rmsz2p30 = seed->GetCMeanSigmaZ2p30();
289 Float_t rmsy2p30R = seed->GetCMeanSigmaY2p30R();
290 Float_t rmsz2p30R = seed->GetCMeanSigmaZ2p30R();
292 AliExternalTrackParam param(*seed);
293 (*fDebugStreamer)<<"ErrParam"<<
300 "rmsy2p30="<<rmsy2p30<<
301 "rmsz2p30="<<rmsz2p30<<
302 "rmsy2p30R="<<rmsy2p30R<<
303 "rmsz2p30R="<<rmsz2p30R<<
307 if (rdistance2>16) return 3;
310 if ((rdistancey2>9. || rdistancez2>9.) && cluster->GetType()==0)
311 return 2; //suspisiouce - will be changed
313 if ((rdistancey2>6.25 || rdistancez2>6.25) && cluster->GetType()>0)
314 // strict cut on overlaped cluster
315 return 2; //suspisiouce - will be changed
317 if ( (rdistancey2>1. || rdistancez2>6.25 )
318 && cluster->GetType()<0){
319 seed->SetNFoundable(seed->GetNFoundable()-1);
328 //_____________________________________________________________________________
329 AliTPCtrackerMI::AliTPCtrackerMI(const AliTPCParam *par):
331 fkNIS(par->GetNInnerSector()/2),
333 fkNOS(par->GetNOuterSector()/2),
350 //---------------------------------------------------------------------
351 // The main TPC tracker constructor
352 //---------------------------------------------------------------------
353 fInnerSec=new AliTPCSector[fkNIS];
354 fOuterSec=new AliTPCSector[fkNOS];
357 for (i=0; i<fkNIS; i++) fInnerSec[i].Setup(par,0);
358 for (i=0; i<fkNOS; i++) fOuterSec[i].Setup(par,1);
361 Int_t nrowlow = par->GetNRowLow();
362 Int_t nrowup = par->GetNRowUp();
365 for (Int_t i=0;i<nrowlow;i++){
366 fXRow[i] = par->GetPadRowRadiiLow(i);
367 fPadLength[i]= par->GetPadPitchLength(0,i);
368 fYMax[i] = fXRow[i]*TMath::Tan(0.5*par->GetInnerAngle());
372 for (Int_t i=0;i<nrowup;i++){
373 fXRow[i+nrowlow] = par->GetPadRowRadiiUp(i);
374 fPadLength[i+nrowlow] = par->GetPadPitchLength(60,i);
375 fYMax[i+nrowlow] = fXRow[i+nrowlow]*TMath::Tan(0.5*par->GetOuterAngle());
378 fDebugStreamer = new TTreeSRedirector("TPCdebug.root");
380 //________________________________________________________________________
381 AliTPCtrackerMI::AliTPCtrackerMI(const AliTPCtrackerMI &t):
402 //------------------------------------
403 // dummy copy constructor
404 //------------------------------------------------------------------
407 AliTPCtrackerMI & AliTPCtrackerMI::operator=(const AliTPCtrackerMI& /*r*/){
408 //------------------------------
410 //--------------------------------------------------------------
413 //_____________________________________________________________________________
414 AliTPCtrackerMI::~AliTPCtrackerMI() {
415 //------------------------------------------------------------------
416 // TPC tracker destructor
417 //------------------------------------------------------------------
424 if (fDebugStreamer) delete fDebugStreamer;
428 void AliTPCtrackerMI::FillESD(TObjArray* arr)
432 //fill esds using updated tracks
434 // write tracks to the event
435 // store index of the track
436 Int_t nseed=arr->GetEntriesFast();
437 //FindKinks(arr,fEvent);
438 for (Int_t i=0; i<nseed; i++) {
439 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
442 // pt->PropagateTo(fParam->GetInnerRadiusLow());
443 if (pt->GetKinkIndex(0)<=0){ //don't propagate daughter tracks
444 pt->PropagateTo(fParam->GetInnerRadiusLow());
447 if (( pt->GetPoints()[2]- pt->GetPoints()[0])>5 && pt->GetPoints()[3]>0.8){
449 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
450 iotrack.SetTPCPoints(pt->GetPoints());
451 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
452 iotrack.SetV0Indexes(pt->GetV0Indexes());
453 // iotrack.SetTPCpid(pt->fTPCr);
454 //iotrack.SetTPCindex(i);
455 fEvent->AddTrack(&iotrack);
459 if ( (pt->GetNumberOfClusters()>70)&& (Float_t(pt->GetNumberOfClusters())/Float_t(pt->GetNFoundable()))>0.55) {
461 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
462 iotrack.SetTPCPoints(pt->GetPoints());
463 //iotrack.SetTPCindex(i);
464 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
465 iotrack.SetV0Indexes(pt->GetV0Indexes());
466 // iotrack.SetTPCpid(pt->fTPCr);
467 fEvent->AddTrack(&iotrack);
471 // short tracks - maybe decays
473 if ( (pt->GetNumberOfClusters()>30) && (Float_t(pt->GetNumberOfClusters())/Float_t(pt->GetNFoundable()))>0.70) {
474 Int_t found,foundable,shared;
475 pt->GetClusterStatistic(0,60,found, foundable,shared,kFALSE);
476 if ( (found>20) && (pt->GetNShared()/float(pt->GetNumberOfClusters())<0.2)){
478 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
479 //iotrack.SetTPCindex(i);
480 iotrack.SetTPCPoints(pt->GetPoints());
481 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
482 iotrack.SetV0Indexes(pt->GetV0Indexes());
483 //iotrack.SetTPCpid(pt->fTPCr);
484 fEvent->AddTrack(&iotrack);
489 if ( (pt->GetNumberOfClusters()>20) && (Float_t(pt->GetNumberOfClusters())/Float_t(pt->GetNFoundable()))>0.8) {
490 Int_t found,foundable,shared;
491 pt->GetClusterStatistic(0,60,found, foundable,shared,kFALSE);
492 if (found<20) continue;
493 if (pt->GetNShared()/float(pt->GetNumberOfClusters())>0.2) continue;
496 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
497 iotrack.SetTPCPoints(pt->GetPoints());
498 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
499 iotrack.SetV0Indexes(pt->GetV0Indexes());
500 //iotrack.SetTPCpid(pt->fTPCr);
501 //iotrack.SetTPCindex(i);
502 fEvent->AddTrack(&iotrack);
505 // short tracks - secondaties
507 if ( (pt->GetNumberOfClusters()>30) ) {
508 Int_t found,foundable,shared;
509 pt->GetClusterStatistic(128,158,found, foundable,shared,kFALSE);
510 if ( (found>20) && (pt->GetNShared()/float(pt->GetNumberOfClusters())<0.2) &&float(found)/float(foundable)>0.8){
512 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
513 iotrack.SetTPCPoints(pt->GetPoints());
514 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
515 iotrack.SetV0Indexes(pt->GetV0Indexes());
516 //iotrack.SetTPCpid(pt->fTPCr);
517 //iotrack.SetTPCindex(i);
518 fEvent->AddTrack(&iotrack);
523 if ( (pt->GetNumberOfClusters()>15)) {
524 Int_t found,foundable,shared;
525 pt->GetClusterStatistic(138,158,found, foundable,shared,kFALSE);
526 if (found<15) continue;
527 if (foundable<=0) continue;
528 if (pt->GetNShared()/float(pt->GetNumberOfClusters())>0.2) continue;
529 if (float(found)/float(foundable)<0.8) continue;
532 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
533 iotrack.SetTPCPoints(pt->GetPoints());
534 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
535 iotrack.SetV0Indexes(pt->GetV0Indexes());
536 // iotrack.SetTPCpid(pt->fTPCr);
537 //iotrack.SetTPCindex(i);
538 fEvent->AddTrack(&iotrack);
543 printf("Number of filled ESDs-\t%d\n",fEvent->GetNumberOfTracks());
550 Double_t AliTPCtrackerMI::ErrY2(AliTPCseed* seed, AliTPCclusterMI * cl){
553 // Use calibrated cluster error from OCDB
555 AliTPCClusterParam * clparam = AliTPCcalibDB::Instance()->GetClusterParam();
557 Float_t z = TMath::Abs(fParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
558 Int_t ctype = cl->GetType();
559 Int_t type = (cl->GetRow()<63) ? 0: (cl->GetRow()>126) ? 1:2;
560 Double_t angle = seed->GetSnp()*seed->GetSnp();
561 angle = TMath::Sqrt(TMath::Abs(angle/(1.-angle)));
562 Double_t erry2 = clparam->GetError0Par(0,type, z,angle);
564 erry2+=0.5; // edge cluster
567 seed->SetErrorY2(erry2);
571 //calculate look-up table at the beginning
572 // static Bool_t ginit = kFALSE;
573 // static Float_t gnoise1,gnoise2,gnoise3;
574 // static Float_t ggg1[10000];
575 // static Float_t ggg2[10000];
576 // static Float_t ggg3[10000];
577 // static Float_t glandau1[10000];
578 // static Float_t glandau2[10000];
579 // static Float_t glandau3[10000];
581 // static Float_t gcor01[500];
582 // static Float_t gcor02[500];
583 // static Float_t gcorp[500];
587 // if (ginit==kFALSE){
588 // for (Int_t i=1;i<500;i++){
589 // Float_t rsigma = float(i)/100.;
590 // gcor02[i] = TMath::Max(0.78 +TMath::Exp(7.4*(rsigma-1.2)),0.6);
591 // gcor01[i] = TMath::Max(0.72 +TMath::Exp(3.36*(rsigma-1.2)),0.6);
592 // gcorp[i] = TMath::Max(TMath::Power((rsigma+0.5),1.5),1.2);
596 // for (Int_t i=3;i<10000;i++){
600 // Float_t amp = float(i);
601 // Float_t padlength =0.75;
602 // gnoise1 = 0.0004/padlength;
603 // Float_t nel = 0.268*amp;
604 // Float_t nprim = 0.155*amp;
605 // ggg1[i] = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.001*nel/(padlength*padlength))/nel;
606 // glandau1[i] = (2.+0.12*nprim)*0.5* (2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
607 // if (glandau1[i]>1) glandau1[i]=1;
608 // glandau1[i]*=padlength*padlength/12.;
612 // gnoise2 = 0.0004/padlength;
614 // nprim = 0.133*amp;
615 // ggg2[i] = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
616 // glandau2[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
617 // if (glandau2[i]>1) glandau2[i]=1;
618 // glandau2[i]*=padlength*padlength/12.;
623 // gnoise3 = 0.0004/padlength;
625 // nprim = 0.133*amp;
626 // ggg3[i] = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
627 // glandau3[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
628 // if (glandau3[i]>1) glandau3[i]=1;
629 // glandau3[i]*=padlength*padlength/12.;
637 // Int_t amp = int(TMath::Abs(cl->GetQ()));
639 // seed->SetErrorY2(1.);
643 // Float_t z = TMath::Abs(fParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
644 // Int_t ctype = cl->GetType();
645 // Float_t padlength= GetPadPitchLength(seed->GetRow());
646 // Double_t angle2 = seed->GetSnp()*seed->GetSnp();
647 // angle2 = angle2/(1-angle2);
649 // //cluster "quality"
650 // Int_t rsigmay = int(100.*cl->GetSigmaY2()/(seed->GetCurrentSigmaY2()));
653 // if (fSectors==fInnerSec){
654 // snoise2 = gnoise1;
655 // res = ggg1[amp]*z+glandau1[amp]*angle2;
656 // if (ctype==0) res *= gcor01[rsigmay];
659 // res*= gcorp[rsigmay];
663 // if (padlength<1.1){
664 // snoise2 = gnoise2;
665 // res = ggg2[amp]*z+glandau2[amp]*angle2;
666 // if (ctype==0) res *= gcor02[rsigmay];
669 // res*= gcorp[rsigmay];
673 // snoise2 = gnoise3;
674 // res = ggg3[amp]*z+glandau3[amp]*angle2;
675 // if (ctype==0) res *= gcor02[rsigmay];
678 // res*= gcorp[rsigmay];
685 // res*=2.4; // overestimate error 2 times
689 // if (res<2*snoise2)
692 // seed->SetErrorY2(res);
700 Double_t AliTPCtrackerMI::ErrZ2(AliTPCseed* seed, AliTPCclusterMI * cl){
703 // Use calibrated cluster error from OCDB
705 AliTPCClusterParam * clparam = AliTPCcalibDB::Instance()->GetClusterParam();
707 Float_t z = TMath::Abs(fParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
708 Int_t ctype = cl->GetType();
709 Int_t type = (cl->GetRow()<63) ? 0: (cl->GetRow()>126) ? 1:2;
711 Double_t angle2 = seed->GetSnp()*seed->GetSnp();
712 angle2 = seed->GetTgl()*seed->GetTgl()*(1+angle2/(1-angle2));
713 Double_t angle = TMath::Sqrt(TMath::Abs(angle2));
714 Double_t errz2 = clparam->GetError0Par(1,type, z,angle);
716 errz2+=0.5; // edge cluster
719 seed->SetErrorZ2(errz2);
725 // //seed->SetErrorY2(0.1);
727 // //calculate look-up table at the beginning
728 // static Bool_t ginit = kFALSE;
729 // static Float_t gnoise1,gnoise2,gnoise3;
730 // static Float_t ggg1[10000];
731 // static Float_t ggg2[10000];
732 // static Float_t ggg3[10000];
733 // static Float_t glandau1[10000];
734 // static Float_t glandau2[10000];
735 // static Float_t glandau3[10000];
737 // static Float_t gcor01[1000];
738 // static Float_t gcor02[1000];
739 // static Float_t gcorp[1000];
743 // if (ginit==kFALSE){
744 // for (Int_t i=1;i<1000;i++){
745 // Float_t rsigma = float(i)/100.;
746 // gcor02[i] = TMath::Max(0.81 +TMath::Exp(6.8*(rsigma-1.2)),0.6);
747 // gcor01[i] = TMath::Max(0.72 +TMath::Exp(2.04*(rsigma-1.2)),0.6);
748 // gcorp[i] = TMath::Max(TMath::Power((rsigma+0.5),1.5),1.2);
752 // for (Int_t i=3;i<10000;i++){
756 // Float_t amp = float(i);
757 // Float_t padlength =0.75;
758 // gnoise1 = 0.0004/padlength;
759 // Float_t nel = 0.268*amp;
760 // Float_t nprim = 0.155*amp;
761 // ggg1[i] = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.001*nel/(padlength*padlength))/nel;
762 // glandau1[i] = (2.+0.12*nprim)*0.5* (2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
763 // if (glandau1[i]>1) glandau1[i]=1;
764 // glandau1[i]*=padlength*padlength/12.;
768 // gnoise2 = 0.0004/padlength;
770 // nprim = 0.133*amp;
771 // ggg2[i] = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
772 // glandau2[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
773 // if (glandau2[i]>1) glandau2[i]=1;
774 // glandau2[i]*=padlength*padlength/12.;
779 // gnoise3 = 0.0004/padlength;
781 // nprim = 0.133*amp;
782 // ggg3[i] = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
783 // glandau3[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
784 // if (glandau3[i]>1) glandau3[i]=1;
785 // glandau3[i]*=padlength*padlength/12.;
793 // Int_t amp = int(TMath::Abs(cl->GetQ()));
795 // seed->SetErrorY2(1.);
799 // Float_t z = TMath::Abs(fParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
800 // Int_t ctype = cl->GetType();
801 // Float_t padlength= GetPadPitchLength(seed->GetRow());
803 // Double_t angle2 = seed->GetSnp()*seed->GetSnp();
804 // // if (angle2<0.6) angle2 = 0.6;
805 // angle2 = seed->GetTgl()*seed->GetTgl()*(1+angle2/(1-angle2));
807 // //cluster "quality"
808 // Int_t rsigmaz = int(100.*cl->GetSigmaZ2()/(seed->GetCurrentSigmaZ2()));
811 // if (fSectors==fInnerSec){
812 // snoise2 = gnoise1;
813 // res = ggg1[amp]*z+glandau1[amp]*angle2;
814 // if (ctype==0) res *= gcor01[rsigmaz];
817 // res*= gcorp[rsigmaz];
821 // if (padlength<1.1){
822 // snoise2 = gnoise2;
823 // res = ggg2[amp]*z+glandau2[amp]*angle2;
824 // if (ctype==0) res *= gcor02[rsigmaz];
827 // res*= gcorp[rsigmaz];
831 // snoise2 = gnoise3;
832 // res = ggg3[amp]*z+glandau3[amp]*angle2;
833 // if (ctype==0) res *= gcor02[rsigmaz];
836 // res*= gcorp[rsigmaz];
845 // if ((ctype<0) &&<70){
850 // if (res<2*snoise2)
852 // if (res>3) res =3;
853 // seed->SetErrorZ2(res);
861 void AliTPCtrackerMI::RotateToLocal(AliTPCseed *seed)
863 //rotate to track "local coordinata
864 Float_t x = seed->GetX();
865 Float_t y = seed->GetY();
866 Float_t ymax = x*TMath::Tan(0.5*fSectors->GetAlpha());
869 seed->SetRelativeSector((seed->GetRelativeSector()+1) % fN);
870 if (!seed->Rotate(fSectors->GetAlpha()))
872 } else if (y <-ymax) {
873 seed->SetRelativeSector((seed->GetRelativeSector()-1+fN) % fN);
874 if (!seed->Rotate(-fSectors->GetAlpha()))
882 //_____________________________________________________________________________
883 Double_t AliTPCtrackerMI::F1old(Double_t x1,Double_t y1,
884 Double_t x2,Double_t y2,
885 Double_t x3,Double_t y3)
887 //-----------------------------------------------------------------
888 // Initial approximation of the track curvature
889 //-----------------------------------------------------------------
890 Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
891 Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
892 (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
893 Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
894 (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
896 Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
897 if ( xr*xr+yr*yr<=0.00000000000001) return 100;
898 return -xr*yr/sqrt(xr*xr+yr*yr);
903 //_____________________________________________________________________________
904 Double_t AliTPCtrackerMI::F1(Double_t x1,Double_t y1,
905 Double_t x2,Double_t y2,
906 Double_t x3,Double_t y3)
908 //-----------------------------------------------------------------
909 // Initial approximation of the track curvature
910 //-----------------------------------------------------------------
916 Double_t det = x3*y2-x2*y3;
921 Double_t u = 0.5* (x2*(x2-x3)+y2*(y2-y3))/det;
922 Double_t x0 = x3*0.5-y3*u;
923 Double_t y0 = y3*0.5+x3*u;
924 Double_t c2 = 1/TMath::Sqrt(x0*x0+y0*y0);
930 Double_t AliTPCtrackerMI::F2(Double_t x1,Double_t y1,
931 Double_t x2,Double_t y2,
932 Double_t x3,Double_t y3)
934 //-----------------------------------------------------------------
935 // Initial approximation of the track curvature
936 //-----------------------------------------------------------------
942 Double_t det = x3*y2-x2*y3;
947 Double_t u = 0.5* (x2*(x2-x3)+y2*(y2-y3))/det;
948 Double_t x0 = x3*0.5-y3*u;
949 Double_t y0 = y3*0.5+x3*u;
950 Double_t c2 = 1/TMath::Sqrt(x0*x0+y0*y0);
959 //_____________________________________________________________________________
960 Double_t AliTPCtrackerMI::F2old(Double_t x1,Double_t y1,
961 Double_t x2,Double_t y2,
962 Double_t x3,Double_t y3)
964 //-----------------------------------------------------------------
965 // Initial approximation of the track curvature times center of curvature
966 //-----------------------------------------------------------------
967 Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
968 Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
969 (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
970 Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
971 (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
973 Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
975 return -a/(d*y1-b)*xr/sqrt(xr*xr+yr*yr);
978 //_____________________________________________________________________________
979 Double_t AliTPCtrackerMI::F3(Double_t x1,Double_t y1,
980 Double_t x2,Double_t y2,
981 Double_t z1,Double_t z2)
983 //-----------------------------------------------------------------
984 // Initial approximation of the tangent of the track dip angle
985 //-----------------------------------------------------------------
986 return (z1 - z2)/sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
990 Double_t AliTPCtrackerMI::F3n(Double_t x1,Double_t y1,
991 Double_t x2,Double_t y2,
992 Double_t z1,Double_t z2, Double_t c)
994 //-----------------------------------------------------------------
995 // Initial approximation of the tangent of the track dip angle
996 //-----------------------------------------------------------------
1000 //angle1 = (z1-z2)*c/(TMath::ASin(c*x1-ni)-TMath::ASin(c*x2-ni));
1002 Double_t d = TMath::Sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
1003 if (TMath::Abs(d*c*0.5)>1) return 0;
1004 // Double_t angle2 = TMath::ASin(d*c*0.5);
1005 // Double_t angle2 = AliTPCFastMath::FastAsin(d*c*0.5);
1006 Double_t angle2 = (d*c*0.5>0.1)? TMath::ASin(d*c*0.5): AliTPCFastMath::FastAsin(d*c*0.5);
1008 angle2 = (z1-z2)*c/(angle2*2.);
1012 Bool_t AliTPCtrackerMI::GetProlongation(Double_t x1, Double_t x2, Double_t x[5], Double_t &y, Double_t &z)
1013 {//-----------------------------------------------------------------
1014 // This function find proloncation of a track to a reference plane x=x2.
1015 //-----------------------------------------------------------------
1019 if (TMath::Abs(x[4]*x1 - x[2]) >= 0.999) {
1023 Double_t c1=x[4]*x1 - x[2], r1=sqrt(1.- c1*c1);
1024 Double_t c2=x[4]*x2 - x[2], r2=sqrt(1.- c2*c2);
1028 Double_t dy = dx*(c1+c2)/(r1+r2);
1031 Double_t delta = x[4]*dx*(c1+c2)/(c1*r2 + c2*r1);
1033 if (TMath::Abs(delta)>0.01){
1034 dz = x[3]*TMath::ASin(delta)/x[4];
1036 dz = x[3]*AliTPCFastMath::FastAsin(delta)/x[4];
1039 //dz = x[3]*AliTPCFastMath::FastAsin(delta)/x[4];
1047 Int_t AliTPCtrackerMI::LoadClusters (TTree *tree)
1052 return LoadClusters();
1056 Int_t AliTPCtrackerMI::LoadClusters(TObjArray *arr)
1059 // load clusters to the memory
1060 AliTPCClustersRow *clrow = 0x0;
1061 Int_t lower = arr->LowerBound();
1062 Int_t entries = arr->GetEntriesFast();
1064 for (Int_t i=lower; i<entries; i++) {
1065 clrow = (AliTPCClustersRow*) arr->At(i);
1066 if(!clrow) continue;
1067 if(!clrow->GetArray()) continue;
1071 fParam->AdjustSectorRow(clrow->GetID(),sec,row);
1073 for (Int_t icl=0; icl<clrow->GetArray()->GetEntriesFast(); icl++){
1074 Transform((AliTPCclusterMI*)(clrow->GetArray()->At(icl)));
1077 if (clrow->GetArray()->GetEntriesFast()<=0) continue;
1078 AliTPCtrackerRow * tpcrow=0;
1081 tpcrow = &(fInnerSec[sec%fkNIS][row]);
1085 tpcrow = &(fOuterSec[(sec-fkNIS*2)%fkNOS][row]);
1086 left = (sec-fkNIS*2)/fkNOS;
1089 tpcrow->SetN1(clrow->GetArray()->GetEntriesFast());
1090 tpcrow->SetClusters1(new AliTPCclusterMI[tpcrow->GetN1()]);
1091 for (Int_t i=0;i<tpcrow->GetN1();i++)
1092 tpcrow->SetCluster1(i, *(AliTPCclusterMI*)(clrow->GetArray()->At(i)));
1095 tpcrow->SetN2(clrow->GetArray()->GetEntriesFast());
1096 tpcrow->SetClusters2(new AliTPCclusterMI[tpcrow->GetN2()]);
1097 for (Int_t i=0;i<tpcrow->GetN2();i++)
1098 tpcrow->SetCluster2(i,*(AliTPCclusterMI*)(clrow->GetArray()->At(i)));
1109 Int_t AliTPCtrackerMI::LoadClusters()
1112 // load clusters to the memory
1113 AliTPCClustersRow *clrow= new AliTPCClustersRow;
1114 clrow->SetClass("AliTPCclusterMI");
1116 clrow->GetArray()->ExpandCreateFast(10000);
1118 // TTree * tree = fClustersArray.GetTree();
1120 TTree * tree = fInput;
1121 TBranch * br = tree->GetBranch("Segment");
1122 br->SetAddress(&clrow);
1124 Int_t j=Int_t(tree->GetEntries());
1125 for (Int_t i=0; i<j; i++) {
1129 fParam->AdjustSectorRow(clrow->GetID(),sec,row);
1130 for (Int_t icl=0; icl<clrow->GetArray()->GetEntriesFast(); icl++){
1131 Transform((AliTPCclusterMI*)(clrow->GetArray()->At(icl)));
1134 AliTPCtrackerRow * tpcrow=0;
1137 tpcrow = &(fInnerSec[sec%fkNIS][row]);
1141 tpcrow = &(fOuterSec[(sec-fkNIS*2)%fkNOS][row]);
1142 left = (sec-fkNIS*2)/fkNOS;
1145 tpcrow->SetN1(clrow->GetArray()->GetEntriesFast());
1146 tpcrow->SetClusters1(new AliTPCclusterMI[tpcrow->GetN1()]);
1147 for (Int_t i=0;i<tpcrow->GetN1();i++)
1148 tpcrow->SetCluster1(i, *(AliTPCclusterMI*)(clrow->GetArray()->At(i)));
1151 tpcrow->SetN2(clrow->GetArray()->GetEntriesFast());
1152 tpcrow->SetClusters2(new AliTPCclusterMI[tpcrow->GetN2()]);
1153 for (Int_t i=0;i<tpcrow->GetN2();i++)
1154 tpcrow->SetCluster2(i,*(AliTPCclusterMI*)(clrow->GetArray()->At(i)));
1165 void AliTPCtrackerMI::UnloadClusters()
1168 // unload clusters from the memory
1170 Int_t nrows = fOuterSec->GetNRows();
1171 for (Int_t sec = 0;sec<fkNOS;sec++)
1172 for (Int_t row = 0;row<nrows;row++){
1173 AliTPCtrackerRow* tpcrow = &(fOuterSec[sec%fkNOS][row]);
1175 // if (tpcrow->fClusters1) delete []tpcrow->fClusters1;
1176 // if (tpcrow->fClusters2) delete []tpcrow->fClusters2;
1178 tpcrow->ResetClusters();
1181 nrows = fInnerSec->GetNRows();
1182 for (Int_t sec = 0;sec<fkNIS;sec++)
1183 for (Int_t row = 0;row<nrows;row++){
1184 AliTPCtrackerRow* tpcrow = &(fInnerSec[sec%fkNIS][row]);
1186 // if (tpcrow->fClusters1) delete []tpcrow->fClusters1;
1187 //if (tpcrow->fClusters2) delete []tpcrow->fClusters2;
1189 tpcrow->ResetClusters();
1195 void AliTPCtrackerMI::Transform(AliTPCclusterMI * cluster){
1199 AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ;
1201 AliFatal("Tranformations not in calibDB");
1203 Double_t x[3]={cluster->GetRow(),cluster->GetPad(),cluster->GetTimeBin()};
1204 Int_t i[1]={cluster->GetDetector()};
1205 transform->Transform(x,i,0,1);
1206 if (!AliTPCReconstructor::GetRecoParam()->GetBYMirror()){
1207 if (cluster->GetDetector()%36>17){
1213 // in debug mode check the transformation
1215 if (AliTPCReconstructor::StreamLevel()>1) {
1217 cluster->GetGlobalXYZ(gx);
1218 TTreeSRedirector &cstream = *fDebugStreamer;
1219 cstream<<"Transform"<<
1229 cluster->SetX(x[0]);
1230 cluster->SetY(x[1]);
1231 cluster->SetZ(x[2]);
1237 //if (!fParam->IsGeoRead()) fParam->ReadGeoMatrices();
1238 TGeoHMatrix *mat = fParam->GetClusterMatrix(cluster->GetDetector());
1240 Double_t pos[3]= {cluster->GetX(),cluster->GetY(),cluster->GetZ()};
1241 Double_t posC[3]={cluster->GetX(),cluster->GetY(),cluster->GetZ()};
1242 if (mat) mat->LocalToMaster(pos,posC);
1244 // chack Loading of Geo matrices from GeoManager - TEMPORARY FIX
1246 cluster->SetX(posC[0]);
1247 cluster->SetY(posC[1]);
1248 cluster->SetZ(posC[2]);
1251 //_____________________________________________________________________________
1252 Int_t AliTPCtrackerMI::LoadOuterSectors() {
1253 //-----------------------------------------------------------------
1254 // This function fills outer TPC sectors with clusters.
1255 //-----------------------------------------------------------------
1256 Int_t nrows = fOuterSec->GetNRows();
1258 for (Int_t sec = 0;sec<fkNOS;sec++)
1259 for (Int_t row = 0;row<nrows;row++){
1260 AliTPCtrackerRow* tpcrow = &(fOuterSec[sec%fkNOS][row]);
1261 Int_t sec2 = sec+2*fkNIS;
1263 Int_t ncl = tpcrow->GetN1();
1265 AliTPCclusterMI *c= (tpcrow->GetCluster1(ncl));
1266 index=(((sec2<<8)+row)<<16)+ncl;
1267 tpcrow->InsertCluster(c,index);
1270 ncl = tpcrow->GetN2();
1272 AliTPCclusterMI *c= (tpcrow->GetCluster2(ncl));
1273 index=((((sec2+fkNOS)<<8)+row)<<16)+ncl;
1274 tpcrow->InsertCluster(c,index);
1277 // write indexes for fast acces
1279 for (Int_t i=0;i<510;i++)
1280 tpcrow->SetFastCluster(i,-1);
1281 for (Int_t i=0;i<tpcrow->GetN();i++){
1282 Int_t zi = Int_t((*tpcrow)[i]->GetZ()+255.);
1283 tpcrow->SetFastCluster(zi,i); // write index
1286 for (Int_t i=0;i<510;i++){
1287 if (tpcrow->GetFastCluster(i)<0)
1288 tpcrow->SetFastCluster(i,last);
1290 last = tpcrow->GetFastCluster(i);
1299 //_____________________________________________________________________________
1300 Int_t AliTPCtrackerMI::LoadInnerSectors() {
1301 //-----------------------------------------------------------------
1302 // This function fills inner TPC sectors with clusters.
1303 //-----------------------------------------------------------------
1304 Int_t nrows = fInnerSec->GetNRows();
1306 for (Int_t sec = 0;sec<fkNIS;sec++)
1307 for (Int_t row = 0;row<nrows;row++){
1308 AliTPCtrackerRow* tpcrow = &(fInnerSec[sec%fkNIS][row]);
1311 Int_t ncl = tpcrow->GetN1();
1313 AliTPCclusterMI *c= (tpcrow->GetCluster1(ncl));
1314 index=(((sec<<8)+row)<<16)+ncl;
1315 tpcrow->InsertCluster(c,index);
1318 ncl = tpcrow->GetN2();
1320 AliTPCclusterMI *c= (tpcrow->GetCluster2(ncl));
1321 index=((((sec+fkNIS)<<8)+row)<<16)+ncl;
1322 tpcrow->InsertCluster(c,index);
1325 // write indexes for fast acces
1327 for (Int_t i=0;i<510;i++)
1328 tpcrow->SetFastCluster(i,-1);
1329 for (Int_t i=0;i<tpcrow->GetN();i++){
1330 Int_t zi = Int_t((*tpcrow)[i]->GetZ()+255.);
1331 tpcrow->SetFastCluster(zi,i); // write index
1334 for (Int_t i=0;i<510;i++){
1335 if (tpcrow->GetFastCluster(i)<0)
1336 tpcrow->SetFastCluster(i,last);
1338 last = tpcrow->GetFastCluster(i);
1350 //_________________________________________________________________________
1351 AliTPCclusterMI *AliTPCtrackerMI::GetClusterMI(Int_t index) const {
1352 //--------------------------------------------------------------------
1353 // Return pointer to a given cluster
1354 //--------------------------------------------------------------------
1355 if (index<0) return 0; // no cluster
1356 Int_t sec=(index&0xff000000)>>24;
1357 Int_t row=(index&0x00ff0000)>>16;
1358 Int_t ncl=(index&0x00007fff)>>00;
1360 const AliTPCtrackerRow * tpcrow=0;
1361 AliTPCclusterMI * clrow =0;
1363 if (sec<0 || sec>=fkNIS*4) {
1364 AliWarning(Form("Wrong sector %d",sec));
1369 tpcrow = &(fInnerSec[sec%fkNIS][row]);
1370 if (tpcrow==0) return 0;
1373 if (tpcrow->GetN1()<=ncl) return 0;
1374 clrow = tpcrow->GetClusters1();
1377 if (tpcrow->GetN2()<=ncl) return 0;
1378 clrow = tpcrow->GetClusters2();
1382 tpcrow = &(fOuterSec[(sec-fkNIS*2)%fkNOS][row]);
1383 if (tpcrow==0) return 0;
1385 if (sec-2*fkNIS<fkNOS) {
1386 if (tpcrow->GetN1()<=ncl) return 0;
1387 clrow = tpcrow->GetClusters1();
1390 if (tpcrow->GetN2()<=ncl) return 0;
1391 clrow = tpcrow->GetClusters2();
1395 return &(clrow[ncl]);
1401 Int_t AliTPCtrackerMI::FollowToNext(AliTPCseed& t, Int_t nr) {
1402 //-----------------------------------------------------------------
1403 // This function tries to find a track prolongation to next pad row
1404 //-----------------------------------------------------------------
1406 Double_t x= GetXrow(nr), ymax=GetMaxY(nr);
1407 AliTPCclusterMI *cl=0;
1408 Int_t tpcindex= t.GetClusterIndex2(nr);
1410 // update current shape info every 5 pad-row
1411 // if ( (nr%5==0) || t.GetNumberOfClusters()<2 || (t.fCurrentSigmaY2<0.0001) ){
1415 if (fIteration>0 && tpcindex>=-1){ //if we have already clusters
1417 if (tpcindex==-1) return 0; //track in dead zone
1419 cl = t.GetClusterPointer(nr);
1420 if ( (cl==0) ) cl = GetClusterMI(tpcindex);
1421 t.SetCurrentClusterIndex1(tpcindex);
1424 Int_t relativesector = ((tpcindex&0xff000000)>>24)%18; // if previously accepted cluster in different sector
1425 Float_t angle = relativesector*fSectors->GetAlpha()+fSectors->GetAlphaShift();
1427 if (angle<-TMath::Pi()) angle += 2*TMath::Pi();
1428 if (angle>=TMath::Pi()) angle -= 2*TMath::Pi();
1430 if (TMath::Abs(angle-t.GetAlpha())>0.001){
1431 Double_t rotation = angle-t.GetAlpha();
1432 t.SetRelativeSector(relativesector);
1433 if (!t.Rotate(rotation)) return 0;
1435 if (!t.PropagateTo(x)) return 0;
1437 t.SetCurrentCluster(cl);
1439 Int_t accept = AcceptCluster(&t,t.GetCurrentCluster());
1440 if ((tpcindex&0x8000)==0) accept =0;
1442 //if founded cluster is acceptible
1443 if (cl->IsUsed(11)) { // id cluster is shared inrease uncertainty
1444 t.SetErrorY2(t.GetErrorY2()+0.03);
1445 t.SetErrorZ2(t.GetErrorZ2()+0.03);
1446 t.SetErrorY2(t.GetErrorY2()*3);
1447 t.SetErrorZ2(t.GetErrorZ2()*3);
1449 t.SetNFoundable(t.GetNFoundable()+1);
1450 UpdateTrack(&t,accept);
1455 if (TMath::Abs(t.GetSnp())>AliTPCReconstructor::GetMaxSnpTracker()) return 0; // cut on angle
1457 // not look for new cluster during refitting
1458 t.SetNFoundable(t.GetNFoundable()+1);
1463 // if (TMath::Abs(t.GetSnp())>0.95 || TMath::Abs(x*t.GetC()-t.GetEta())>0.95) return 0;// patch 28 fev 06
1464 Double_t y=t.GetYat(x);
1465 if (TMath::Abs(y)>ymax){
1467 t.SetRelativeSector((t.GetRelativeSector()+1) % fN);
1468 if (!t.Rotate(fSectors->GetAlpha()))
1470 } else if (y <-ymax) {
1471 t.SetRelativeSector((t.GetRelativeSector()-1+fN) % fN);
1472 if (!t.Rotate(-fSectors->GetAlpha()))
1478 if (!t.PropagateTo(x)) {
1479 if (fIteration==0) t.SetRemoval(10);
1483 Double_t z=t.GetZ();
1486 if (!IsActive(t.GetRelativeSector(),nr)) {
1488 t.SetClusterIndex2(nr,-1);
1491 //AliInfo(Form("A - Sector%d phi %f - alpha %f", t.fRelativeSector,y/x, t.GetAlpha()));
1492 Bool_t isActive = IsActive(t.GetRelativeSector(),nr);
1493 Bool_t isActive2 = (nr>=fInnerSec->GetNRows()) ? fOuterSec[t.GetRelativeSector()][nr-fInnerSec->GetNRows()].GetN()>0:fInnerSec[t.GetRelativeSector()][nr].GetN()>0;
1495 if (!isActive || !isActive2) return 0;
1497 const AliTPCtrackerRow &krow=GetRow(t.GetRelativeSector(),nr);
1498 if ( (t.GetSigmaY2()<0) || t.GetSigmaZ2()<0) return 0;
1500 Double_t roadz = 1.;
1502 if (TMath::Abs(TMath::Abs(y)-ymax)<krow.GetDeadZone()){
1504 t.SetClusterIndex2(nr,-1);
1509 if (TMath::Abs(z)<(AliTPCReconstructor::GetCtgRange()*x+10) && TMath::Abs(z)<fParam->GetZLength(0) && (TMath::Abs(t.GetSnp())<AliTPCReconstructor::GetMaxSnpTracker()))
1510 t.SetNFoundable(t.GetNFoundable()+1);
1516 // cl = krow.FindNearest2(y+10.,z,roady,roadz,index);
1517 cl = krow.FindNearest2(y,z,roady,roadz,index);
1518 if (cl) t.SetCurrentClusterIndex1(krow.GetIndex(index));
1521 t.SetCurrentCluster(cl);
1523 if (fIteration==2&&cl->IsUsed(10)) return 0;
1524 Int_t accept = AcceptCluster(&t,t.GetCurrentCluster());
1525 if (fIteration==2&&cl->IsUsed(11)) {
1526 t.SetErrorY2(t.GetErrorY2()+0.03);
1527 t.SetErrorZ2(t.GetErrorZ2()+0.03);
1528 t.SetErrorY2(t.GetErrorY2()*3);
1529 t.SetErrorZ2(t.GetErrorZ2()*3);
1532 if (t.fCurrentCluster->IsUsed(10)){
1537 if (t.fNShared>0.7*t.GetNumberOfClusters()) {
1543 if (accept<3) UpdateTrack(&t,accept);
1546 if ( fIteration==0 && t.GetNFoundable()*0.5 > t.GetNumberOfClusters()) t.SetRemoval(10);
1552 Int_t AliTPCtrackerMI::FollowToNextFast(AliTPCseed& t, Int_t nr) {
1553 //-----------------------------------------------------------------
1554 // This function tries to find a track prolongation to next pad row
1555 //-----------------------------------------------------------------
1557 Double_t x= GetXrow(nr), ymax=GetMaxY(nr);
1559 if (!t.GetProlongation(x,y,z)) {
1565 if (TMath::Abs(y)>ymax){
1568 t.SetRelativeSector((t.GetRelativeSector()+1) % fN);
1569 if (!t.Rotate(fSectors->GetAlpha()))
1571 } else if (y <-ymax) {
1572 t.SetRelativeSector((t.GetRelativeSector()-1+fN) % fN);
1573 if (!t.Rotate(-fSectors->GetAlpha()))
1576 if (!t.PropagateTo(x)) {
1579 t.GetProlongation(x,y,z);
1582 // update current shape info every 2 pad-row
1583 if ( (nr%2==0) || t.GetNumberOfClusters()<2 || (t.GetCurrentSigmaY2()<0.0001) ){
1584 // t.fCurrentSigmaY = GetSigmaY(&t);
1585 //t.fCurrentSigmaZ = GetSigmaZ(&t);
1589 AliTPCclusterMI *cl=0;
1594 const AliTPCtrackerRow &krow=GetRow(t.GetRelativeSector(),nr);
1595 if ( (t.GetSigmaY2()<0) || t.GetSigmaZ2()<0) return 0;
1597 Double_t roadz = 1.;
1600 if (TMath::Abs(TMath::Abs(y)-ymax)<krow.GetDeadZone()){
1602 t.SetClusterIndex2(row,-1);
1607 if (TMath::Abs(z)>(AliTPCReconstructor::GetCtgRange()*x+10)) t.SetClusterIndex2(row,-1);
1611 if ((cl==0)&&(krow)) {
1612 // cl = krow.FindNearest2(y+10,z,roady,roadz,index);
1613 cl = krow.FindNearest2(y,z,roady,roadz,index);
1615 if (cl) t.SetCurrentClusterIndex1(krow.GetIndex(index));
1619 t.SetCurrentCluster(cl);
1620 // Int_t accept = AcceptCluster(&t,t.fCurrentCluster);
1622 t.SetClusterIndex2(row,index);
1623 t.SetClusterPointer(row, cl);
1631 //_________________________________________________________________________
1632 Bool_t AliTPCtrackerMI::GetTrackPoint(Int_t index, AliTrackPoint &p ) const
1634 // Get track space point by index
1635 // return false in case the cluster doesn't exist
1636 AliTPCclusterMI *cl = GetClusterMI(index);
1637 if (!cl) return kFALSE;
1638 Int_t sector = (index&0xff000000)>>24;
1639 // Int_t row = (index&0x00ff0000)>>16;
1641 // xyz[0] = fParam->GetPadRowRadii(sector,row);
1642 xyz[0] = cl->GetX();
1643 xyz[1] = cl->GetY();
1644 xyz[2] = cl->GetZ();
1646 fParam->AdjustCosSin(sector,cos,sin);
1647 Float_t x = cos*xyz[0]-sin*xyz[1];
1648 Float_t y = cos*xyz[1]+sin*xyz[0];
1650 Float_t sigmaY2 = 0.027*cl->GetSigmaY2();
1651 if (sector < fParam->GetNInnerSector()) sigmaY2 *= 2.07;
1652 Float_t sigmaZ2 = 0.066*cl->GetSigmaZ2();
1653 if (sector < fParam->GetNInnerSector()) sigmaZ2 *= 1.77;
1654 cov[0] = sin*sin*sigmaY2;
1655 cov[1] = -sin*cos*sigmaY2;
1657 cov[3] = cos*cos*sigmaY2;
1660 p.SetXYZ(x,y,xyz[2],cov);
1661 AliGeomManager::ELayerID iLayer;
1663 if (sector < fParam->GetNInnerSector()) {
1664 iLayer = AliGeomManager::kTPC1;
1668 iLayer = AliGeomManager::kTPC2;
1669 idet = sector - fParam->GetNInnerSector();
1671 UShort_t volid = AliGeomManager::LayerToVolUID(iLayer,idet);
1672 p.SetVolumeID(volid);
1678 Int_t AliTPCtrackerMI::UpdateClusters(AliTPCseed& t, Int_t nr) {
1679 //-----------------------------------------------------------------
1680 // This function tries to find a track prolongation to next pad row
1681 //-----------------------------------------------------------------
1682 t.SetCurrentCluster(0);
1683 t.SetCurrentClusterIndex1(0);
1685 Double_t xt=t.GetX();
1686 Int_t row = GetRowNumber(xt)-1;
1687 Double_t ymax= GetMaxY(nr);
1689 if (row < nr) return 1; // don't prolongate if not information until now -
1690 // if (TMath::Abs(t.GetSnp())>0.9 && t.GetNumberOfClusters()>40. && fIteration!=2) {
1692 // return 0; // not prolongate strongly inclined tracks
1694 // if (TMath::Abs(t.GetSnp())>0.95) {
1696 // return 0; // not prolongate strongly inclined tracks
1697 // }// patch 28 fev 06
1699 Double_t x= GetXrow(nr);
1701 //t.PropagateTo(x+0.02);
1702 //t.PropagateTo(x+0.01);
1703 if (!t.PropagateTo(x)){
1710 if (TMath::Abs(y)>ymax){
1712 t.SetRelativeSector((t.GetRelativeSector()+1) % fN);
1713 if (!t.Rotate(fSectors->GetAlpha()))
1715 } else if (y <-ymax) {
1716 t.SetRelativeSector((t.GetRelativeSector()-1+fN) % fN);
1717 if (!t.Rotate(-fSectors->GetAlpha()))
1720 // if (!t.PropagateTo(x)){
1727 if (TMath::Abs(t.GetSnp())>AliTPCReconstructor::GetMaxSnpTracker()) return 0;
1729 if (!IsActive(t.GetRelativeSector(),nr)) {
1731 t.SetClusterIndex2(nr,-1);
1734 //AliInfo(Form("A - Sector%d phi %f - alpha %f", t.fRelativeSector,y/x, t.GetAlpha()));
1736 AliTPCtrackerRow &krow=GetRow(t.GetRelativeSector(),nr);
1738 if (TMath::Abs(TMath::Abs(y)-ymax)<krow.GetDeadZone()){
1740 t.SetClusterIndex2(nr,-1);
1745 if (TMath::Abs(t.GetZ())<(AliTPCReconstructor::GetCtgRange()*t.GetX()+10) && (TMath::Abs(t.GetSnp())<AliTPCReconstructor::GetMaxSnpTracker())) t.SetNFoundable(t.GetNFoundable()+1);
1751 if ( (nr%2==0) || t.GetNumberOfClusters()<2){
1752 // t.fCurrentSigmaY = GetSigmaY(&t);
1753 //t.fCurrentSigmaZ = GetSigmaZ(&t);
1757 AliTPCclusterMI *cl=0;
1760 Double_t roady = 1.;
1761 Double_t roadz = 1.;
1765 index = t.GetClusterIndex2(nr);
1766 if ( (index>0) && (index&0x8000)==0){
1767 cl = t.GetClusterPointer(nr);
1768 if ( (cl==0) && (index>0)) cl = GetClusterMI(index);
1769 t.SetCurrentClusterIndex1(index);
1771 t.SetCurrentCluster(cl);
1777 // if (index<0) return 0;
1778 UInt_t uindex = TMath::Abs(index);
1781 //cl = krow.FindNearest2(y+10,z,roady,roadz,uindex);
1782 cl = krow.FindNearest2(y,z,roady,roadz,uindex);
1785 if (cl) t.SetCurrentClusterIndex1(krow.GetIndex(uindex));
1786 t.SetCurrentCluster(cl);
1792 Int_t AliTPCtrackerMI::FollowToNextCluster(AliTPCseed & t, Int_t nr) {
1793 //-----------------------------------------------------------------
1794 // This function tries to find a track prolongation to next pad row
1795 //-----------------------------------------------------------------
1797 //update error according neighborhoud
1799 if (t.GetCurrentCluster()) {
1801 Int_t accept = AcceptCluster(&t,t.GetCurrentCluster());
1803 if (t.GetCurrentCluster()->IsUsed(10)){
1808 t.SetNShared(t.GetNShared()+1);
1809 if (t.GetNShared()>0.7*t.GetNumberOfClusters()) {
1814 if (fIteration>0) accept = 0;
1815 if (accept<3) UpdateTrack(&t,accept);
1819 if ( ( (t.GetSigmaY2()+t.GetSigmaZ2())>0.16)&& t.GetNumberOfClusters()>18) t.SetRemoval(10);
1820 if ( t.GetChi2()/t.GetNumberOfClusters()>6 &&t.GetNumberOfClusters()>18) t.SetRemoval(10);
1822 if (( (t.GetNFoundable()*0.5 > t.GetNumberOfClusters()) || t.GetNoCluster()>15)) t.SetRemoval(10);
1830 //_____________________________________________________________________________
1831 Int_t AliTPCtrackerMI::FollowProlongation(AliTPCseed& t, Int_t rf, Int_t step) {
1832 //-----------------------------------------------------------------
1833 // This function tries to find a track prolongation.
1834 //-----------------------------------------------------------------
1835 Double_t xt=t.GetX();
1837 Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
1838 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
1839 if (alpha < 0. ) alpha += 2.*TMath::Pi();
1841 t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN);
1843 Int_t first = GetRowNumber(xt)-1;
1844 for (Int_t nr= first; nr>=rf; nr-=step) {
1846 if (t.GetKinkIndexes()[0]>0){
1847 for (Int_t i=0;i<3;i++){
1848 Int_t index = t.GetKinkIndexes()[i];
1849 if (index==0) break;
1850 if (index<0) continue;
1852 AliKink * kink = (AliKink*)fEvent->GetKink(index-1);
1854 printf("PROBLEM\n");
1857 Int_t kinkrow = kink->GetTPCRow0()+2+Int_t(0.5/(0.05+kink->GetAngle(2)));
1859 AliExternalTrackParam paramd(t);
1860 kink->SetDaughter(paramd);
1861 kink->SetStatus(2,5);
1868 if (nr==80) t.UpdateReference();
1869 if (nr<fInnerSec->GetNRows())
1870 fSectors = fInnerSec;
1872 fSectors = fOuterSec;
1873 if (FollowToNext(t,nr)==0)
1882 //_____________________________________________________________________________
1883 Int_t AliTPCtrackerMI::FollowProlongationFast(AliTPCseed& t, Int_t rf, Int_t step) {
1884 //-----------------------------------------------------------------
1885 // This function tries to find a track prolongation.
1886 //-----------------------------------------------------------------
1887 Double_t xt=t.GetX();
1889 Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
1890 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
1891 if (alpha < 0. ) alpha += 2.*TMath::Pi();
1892 t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN);
1894 for (Int_t nr=GetRowNumber(xt)-1; nr>=rf; nr-=step) {
1896 if (FollowToNextFast(t,nr)==0)
1897 if (!t.IsActive()) return 0;
1907 Int_t AliTPCtrackerMI::FollowBackProlongation(AliTPCseed& t, Int_t rf) {
1908 //-----------------------------------------------------------------
1909 // This function tries to find a track prolongation.
1910 //-----------------------------------------------------------------
1912 Double_t xt=t.GetX();
1913 Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
1914 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
1915 if (alpha < 0. ) alpha += 2.*TMath::Pi();
1916 t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN);
1918 Int_t first = t.GetFirstPoint();
1919 if (first<GetRowNumber(xt)+1) first = GetRowNumber(xt)+1;
1921 if (first<0) first=0;
1922 for (Int_t nr=first; nr<=rf; nr++) {
1923 // if ( (TMath::Abs(t.GetSnp())>0.95)) break;//patch 28 fev 06
1924 if (t.GetKinkIndexes()[0]<0){
1925 for (Int_t i=0;i<3;i++){
1926 Int_t index = t.GetKinkIndexes()[i];
1927 if (index==0) break;
1928 if (index>0) continue;
1929 index = TMath::Abs(index);
1930 AliKink * kink = (AliKink*)fEvent->GetKink(index-1);
1932 printf("PROBLEM\n");
1935 Int_t kinkrow = kink->GetTPCRow0()-2-Int_t(0.5/(0.05+kink->GetAngle(2)));
1937 AliExternalTrackParam paramm(t);
1938 kink->SetMother(paramm);
1939 kink->SetStatus(2,1);
1946 if (nr<fInnerSec->GetNRows())
1947 fSectors = fInnerSec;
1949 fSectors = fOuterSec;
1960 Float_t AliTPCtrackerMI::OverlapFactor(AliTPCseed * s1, AliTPCseed * s2, Int_t &sum1, Int_t & sum2)
1968 Float_t dz2 =(s1->GetZ() - s2->GetZ());
1971 Float_t dy2 =TMath::Abs((s1->GetY() - s2->GetY()));
1973 Float_t distance = TMath::Sqrt(dz2+dy2);
1974 if (distance>4.) return 0; // if there are far away - not overlap - to reduce combinatorics
1977 Int_t firstpoint = TMath::Min(s1->GetFirstPoint(),s2->GetFirstPoint());
1978 Int_t lastpoint = TMath::Max(s1->GetLastPoint(),s2->GetLastPoint());
1983 if (firstpoint>lastpoint) {
1984 firstpoint =lastpoint;
1989 for (Int_t i=firstpoint-1;i<lastpoint+1;i++){
1990 if (s1->GetClusterIndex2(i)>0) sum1++;
1991 if (s2->GetClusterIndex2(i)>0) sum2++;
1992 if (s1->GetClusterIndex2(i)==s2->GetClusterIndex2(i) && s1->GetClusterIndex2(i)>0) {
1996 if (sum<5) return 0;
1998 Float_t summin = TMath::Min(sum1+1,sum2+1);
1999 Float_t ratio = (sum+1)/Float_t(summin);
2003 void AliTPCtrackerMI::SignShared(AliTPCseed * s1, AliTPCseed * s2)
2007 Float_t thetaCut = 0.2;//+10.*TMath::Sqrt(s1->GetSigmaTglZ()+ s2->GetSigmaTglZ());
2008 if (TMath::Abs(s1->GetTgl()-s2->GetTgl())>thetaCut) return;
2009 Float_t minCl = TMath::Min(s1->GetNumberOfClusters(),s2->GetNumberOfClusters());
2010 Int_t cutN0 = TMath::Max(5,TMath::Nint(0.1*minCl));
2015 //Int_t firstpoint = TMath::Max(s1->GetFirstPoint(),s2->GetFirstPoint());
2016 //Int_t lastpoint = TMath::Min(s1->GetLastPoint(),s2->GetLastPoint());
2017 Int_t firstpoint = 0;
2018 Int_t lastpoint = 160;
2020 // if (firstpoint>=lastpoint-5) return;;
2022 for (Int_t i=firstpoint;i<lastpoint;i++){
2023 // if ( (s1->GetClusterIndex2(i)&0xFFFF8FFF)==(s2->GetClusterIndex2(i)&0xFFFF8FFF) && s1->GetClusterIndex2(i)>0) {
2024 if ( (s1->GetClusterIndex2(i))==(s2->GetClusterIndex2(i)) && s1->GetClusterIndex2(i)>0) {
2026 s1->SetSharedMapBit(i, kTRUE);
2027 s2->SetSharedMapBit(i, kTRUE);
2029 if (s1->GetClusterIndex2(i)>0)
2030 s1->SetClusterMapBit(i, kTRUE);
2032 if (sumshared>cutN0){
2035 for (Int_t i=firstpoint;i<lastpoint;i++){
2036 // if ( (s1->GetClusterIndex2(i)&0xFFFF8FFF)==(s2->GetClusterIndex2(i)&0xFFFF8FFF) && s1->GetClusterIndex2(i)>0) {
2037 if ( (s1->GetClusterIndex2(i))==(s2->GetClusterIndex2(i)) && s1->GetClusterIndex2(i)>0) {
2038 AliTPCTrackerPoint *p1 = s1->GetTrackPoint(i);
2039 AliTPCTrackerPoint *p2 = s2->GetTrackPoint(i);;
2040 if (s1->IsActive()&&s2->IsActive()){
2041 p1->SetShared(kTRUE);
2042 p2->SetShared(kTRUE);
2048 if (sumshared>cutN0){
2049 for (Int_t i=0;i<4;i++){
2050 if (s1->GetOverlapLabel(3*i)==0){
2051 s1->SetOverlapLabel(3*i, s2->GetLabel());
2052 s1->SetOverlapLabel(3*i+1,sumshared);
2053 s1->SetOverlapLabel(3*i+2,s2->GetUniqueID());
2057 for (Int_t i=0;i<4;i++){
2058 if (s2->GetOverlapLabel(3*i)==0){
2059 s2->SetOverlapLabel(3*i, s1->GetLabel());
2060 s2->SetOverlapLabel(3*i+1,sumshared);
2061 s2->SetOverlapLabel(3*i+2,s1->GetUniqueID());
2068 void AliTPCtrackerMI::SignShared(TObjArray * arr)
2071 //sort trackss according sectors
2073 for (Int_t i=0; i<arr->GetEntriesFast(); i++) {
2074 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2076 //if (pt) RotateToLocal(pt);
2080 arr->Sort(); // sorting according relative sectors
2081 arr->Expand(arr->GetEntries());
2084 Int_t nseed=arr->GetEntriesFast();
2085 for (Int_t i=0; i<nseed; i++) {
2086 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2088 for (Int_t j=0;j<=12;j++){
2089 pt->SetOverlapLabel(j,0);
2092 for (Int_t i=0; i<nseed; i++) {
2093 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2095 if (pt->GetRemoval()>10) continue;
2096 for (Int_t j=i+1; j<nseed; j++){
2097 AliTPCseed *pt2=(AliTPCseed*)arr->UncheckedAt(j);
2098 if (TMath::Abs(pt->GetRelativeSector()-pt2->GetRelativeSector())>1) continue;
2100 if (pt2->GetRemoval()<=10) {
2101 //if ( TMath::Abs(pt->GetRelativeSector()-pt2->GetRelativeSector())>0) break;
2109 void AliTPCtrackerMI::SortTracks(TObjArray * arr, Int_t mode) const
2112 //sort tracks in array according mode criteria
2113 Int_t nseed = arr->GetEntriesFast();
2114 for (Int_t i=0; i<nseed; i++) {
2115 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2126 void AliTPCtrackerMI::RemoveUsed2(TObjArray * arr, Float_t factor1, Float_t factor2, Int_t minimal)
2129 // Loop over all tracks and remove overlaped tracks (with lower quality)
2131 // 1. Unsign clusters
2132 // 2. Sort tracks according quality
2133 // Quality is defined by the number of cluster between first and last points
2135 // 3. Loop over tracks - decreasing quality order
2136 // a.) remove - If the fraction of shared cluster less than factor (1- n or 2)
2137 // b.) remove - If the minimal number of clusters less than minimal and not ITS
2138 // c.) if track accepted - sign clusters
2140 //Called in - AliTPCtrackerMI::Clusters2Tracks()
2141 // - AliTPCtrackerMI::PropagateBack()
2142 // - AliTPCtrackerMI::RefitInward()
2148 Int_t nseed = arr->GetEntriesFast();
2149 Float_t * quality = new Float_t[nseed];
2150 Int_t * indexes = new Int_t[nseed];
2154 for (Int_t i=0; i<nseed; i++) {
2155 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2160 pt->UpdatePoints(); //select first last max dens points
2161 Float_t * points = pt->GetPoints();
2162 if (points[3]<0.8) quality[i] =-1;
2163 quality[i] = (points[2]-points[0])+pt->GetNumberOfClusters();
2164 //prefer high momenta tracks if overlaps
2165 quality[i] *= TMath::Sqrt(TMath::Abs(pt->Pt())+0.5);
2167 TMath::Sort(nseed,quality,indexes);
2170 for (Int_t itrack=0; itrack<nseed; itrack++) {
2171 Int_t trackindex = indexes[itrack];
2172 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(trackindex);
2175 if (quality[trackindex]<0){
2177 delete arr->RemoveAt(trackindex);
2180 arr->RemoveAt(trackindex);
2186 Int_t first = Int_t(pt->GetPoints()[0]);
2187 Int_t last = Int_t(pt->GetPoints()[2]);
2188 Double_t factor = (pt->GetBConstrain()) ? factor1: factor2;
2190 Int_t found,foundable,shared;
2191 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
2192 // pt->GetClusterStatistic(0,160, found, foundable,shared,kFALSE);
2193 Bool_t itsgold =kFALSE;
2196 if (pt->GetESD()->GetITSclusters(dummy)>4) itsgold= kTRUE;
2200 if (Float_t(shared+1)/Float_t(found+1)>factor){
2201 if (pt->GetKinkIndexes()[0]!=0) continue; //don't remove tracks - part of the kinks
2202 delete arr->RemoveAt(trackindex);
2205 if (pt->GetNumberOfClusters()<50&&(found-0.5*shared)<minimal){ //remove short tracks
2206 if (pt->GetKinkIndexes()[0]!=0) continue; //don't remove tracks - part of the kinks
2207 delete arr->RemoveAt(trackindex);
2213 //if (sharedfactor>0.4) continue;
2214 if (pt->GetKinkIndexes()[0]>0) continue;
2215 //Remove tracks with undefined properties - seems
2216 if (pt->GetSigmaY2()<kAlmost0) continue; // ? what is the origin ?
2218 for (Int_t i=first; i<last; i++) {
2219 Int_t index=pt->GetClusterIndex2(i);
2220 // if (index<0 || index&0x8000 ) continue;
2221 if (index<0 || index&0x8000 ) continue;
2222 AliTPCclusterMI *c= pt->GetClusterPointer(i);
2229 Info("RemoveUsed","\n*****\nNumber of good tracks after shared removal\t%d\n",fNtracks);
2235 void AliTPCtrackerMI::UnsignClusters()
2238 // loop over all clusters and unsign them
2241 for (Int_t sec=0;sec<fkNIS;sec++){
2242 for (Int_t row=0;row<fInnerSec->GetNRows();row++){
2243 AliTPCclusterMI *cl = fInnerSec[sec][row].GetClusters1();
2244 for (Int_t icl =0;icl< fInnerSec[sec][row].GetN1();icl++)
2245 // if (cl[icl].IsUsed(10))
2247 cl = fInnerSec[sec][row].GetClusters2();
2248 for (Int_t icl =0;icl< fInnerSec[sec][row].GetN2();icl++)
2249 //if (cl[icl].IsUsed(10))
2254 for (Int_t sec=0;sec<fkNOS;sec++){
2255 for (Int_t row=0;row<fOuterSec->GetNRows();row++){
2256 AliTPCclusterMI *cl = fOuterSec[sec][row].GetClusters1();
2257 for (Int_t icl =0;icl< fOuterSec[sec][row].GetN1();icl++)
2258 //if (cl[icl].IsUsed(10))
2260 cl = fOuterSec[sec][row].GetClusters2();
2261 for (Int_t icl =0;icl< fOuterSec[sec][row].GetN2();icl++)
2262 //if (cl[icl].IsUsed(10))
2271 void AliTPCtrackerMI::SignClusters(TObjArray * arr, Float_t fnumber, Float_t fdensity)
2274 //sign clusters to be "used"
2276 // snumber and sdensity sign number of sigmas - bellow mean value to be accepted
2277 // loop over "primaries"
2291 Int_t nseed = arr->GetEntriesFast();
2292 for (Int_t i=0; i<nseed; i++) {
2293 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2297 if (!(pt->IsActive())) continue;
2298 Float_t dens = pt->GetNumberOfClusters()/Float_t(pt->GetNFoundable());
2299 if ( (dens>0.7) && (pt->GetNumberOfClusters()>70)){
2301 sumdens2+= dens*dens;
2302 sumn += pt->GetNumberOfClusters();
2303 sumn2 += pt->GetNumberOfClusters()*pt->GetNumberOfClusters();
2304 Float_t chi2 = pt->GetChi2()/pt->GetNumberOfClusters();
2307 sumchi2 +=chi2*chi2;
2312 Float_t mdensity = 0.9;
2313 Float_t meann = 130;
2314 Float_t meanchi = 1;
2315 Float_t sdensity = 0.1;
2316 Float_t smeann = 10;
2317 Float_t smeanchi =0.4;
2321 mdensity = sumdens/sum;
2323 meanchi = sumchi/sum;
2325 sdensity = sumdens2/sum-mdensity*mdensity;
2327 sdensity = TMath::Sqrt(sdensity);
2331 smeann = sumn2/sum-meann*meann;
2333 smeann = TMath::Sqrt(smeann);
2337 smeanchi = sumchi2/sum - meanchi*meanchi;
2339 smeanchi = TMath::Sqrt(smeanchi);
2345 //REMOVE SHORT DELTAS or tracks going out of sensitive volume of TPC
2347 for (Int_t i=0; i<nseed; i++) {
2348 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2352 if (pt->GetBSigned()) continue;
2353 if (pt->GetBConstrain()) continue;
2354 //if (!(pt->IsActive())) continue;
2356 Int_t found,foundable,shared;
2357 pt->GetClusterStatistic(0,160,found, foundable,shared);
2358 if (shared/float(found)>0.3) {
2359 if (shared/float(found)>0.9 ){
2360 //delete arr->RemoveAt(i);
2365 Bool_t isok =kFALSE;
2366 if ( (pt->GetNShared()/pt->GetNumberOfClusters()<0.5) &&pt->GetNumberOfClusters()>60)
2368 if ((TMath::Abs(1/pt->GetC())<100.) && (pt->GetNShared()/pt->GetNumberOfClusters()<0.7))
2370 if (TMath::Abs(pt->GetZ()/pt->GetX())>1.1)
2372 if ( (TMath::Abs(pt->GetSnp()>0.7) && pt->GetD(0,0)>60.))
2376 for (Int_t i=0; i<160; i++) {
2377 Int_t index=pt->GetClusterIndex2(i);
2378 if (index<0) continue;
2379 AliTPCclusterMI *c= pt->GetClusterPointer(i);
2381 //if (!(c->IsUsed(10))) c->Use();
2388 Double_t maxchi = meanchi+2.*smeanchi;
2390 for (Int_t i=0; i<nseed; i++) {
2391 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2395 //if (!(pt->IsActive())) continue;
2396 if (pt->GetBSigned()) continue;
2397 Double_t chi = pt->GetChi2()/pt->GetNumberOfClusters();
2398 if (chi>maxchi) continue;
2401 Float_t dens = pt->GetNumberOfClusters()/Float_t(pt->GetNFoundable());
2403 //sign only tracks with enoug big density at the beginning
2405 if ((pt->GetDensityFirst(40)<0.75) && pt->GetNumberOfClusters()<meann) continue;
2408 Double_t mindens = TMath::Max(double(mdensity-sdensity*fdensity*bfactor),0.65);
2409 Double_t minn = TMath::Max(Int_t(meann-fnumber*smeann*bfactor),50);
2411 // if (pt->fBConstrain) mindens = TMath::Max(mdensity-sdensity*fdensity*bfactor,0.65);
2412 if ( (pt->GetRemoval()==10) && (pt->GetSnp()>0.8)&&(dens>mindens))
2415 if ((dens>mindens && pt->GetNumberOfClusters()>minn) && chi<maxchi ){
2416 //Int_t noc=pt->GetNumberOfClusters();
2417 pt->SetBSigned(kTRUE);
2418 for (Int_t i=0; i<160; i++) {
2420 Int_t index=pt->GetClusterIndex2(i);
2421 if (index<0) continue;
2422 AliTPCclusterMI *c= pt->GetClusterPointer(i);
2424 // if (!(c->IsUsed(10))) c->Use();
2429 // gLastCheck = nseed;
2437 void AliTPCtrackerMI::StopNotActive(TObjArray * arr, Int_t row0, Float_t th0, Float_t th1, Float_t th2) const
2439 // stop not active tracks
2440 // take th1 as threshold for number of founded to number of foundable on last 10 active rows
2441 // take th2 as threshold for number of founded to number of foundable on last 20 active rows
2442 Int_t nseed = arr->GetEntriesFast();
2444 for (Int_t i=0; i<nseed; i++) {
2445 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2449 if (!(pt->IsActive())) continue;
2450 StopNotActive(pt,row0,th0, th1,th2);
2456 void AliTPCtrackerMI::StopNotActive(AliTPCseed * seed, Int_t row0, Float_t th0, Float_t th1,
2459 // stop not active tracks
2460 // take th1 as threshold for number of founded to number of foundable on last 10 active rows
2461 // take th2 as threshold for number of founded to number of foundable on last 20 active rows
2464 Int_t foundable = 0;
2465 Int_t maxindex = seed->GetLastPoint(); //last foundable row
2466 if (seed->GetNFoundable()*th0 > seed->GetNumberOfClusters()) {
2467 seed->Desactivate(10) ;
2471 for (Int_t i=row0; i<maxindex; i++){
2472 Int_t index = seed->GetClusterIndex2(i);
2473 if (index!=-1) foundable++;
2475 if (foundable<=30) sumgood1++;
2476 if (foundable<=50) {
2483 if (foundable>=30.){
2484 if (sumgood1<(th1*30.)) seed->Desactivate(10);
2487 if (sumgood2<(th2*50.)) seed->Desactivate(10);
2491 Int_t AliTPCtrackerMI::RefitInward(AliESDEvent *event)
2494 // back propagation of ESD tracks
2497 const Int_t kMaxFriendTracks=2000;
2501 //PrepareForProlongation(fSeeds,1);
2502 PropagateForward2(fSeeds);
2503 RemoveUsed2(fSeeds,0.4,0.4,20);
2505 TObjArray arraySeed(fSeeds->GetEntries());
2506 for (Int_t i=0;i<fSeeds->GetEntries();i++) {
2507 arraySeed.AddAt(fSeeds->At(i),i);
2509 SignShared(&arraySeed);
2510 // FindCurling(fSeeds, event,2); // find multi found tracks
2511 FindSplitted(fSeeds, event,2); // find multi found tracks
2514 Int_t nseed = fSeeds->GetEntriesFast();
2515 for (Int_t i=0;i<nseed;i++){
2516 AliTPCseed * seed = (AliTPCseed*) fSeeds->UncheckedAt(i);
2517 if (!seed) continue;
2518 if (seed->GetKinkIndex(0)>0) UpdateKinkQualityD(seed); // update quality informations for kinks
2520 seed->PropagateTo(fParam->GetInnerRadiusLow());
2521 seed->UpdatePoints();
2523 AliESDtrack *esd=event->GetTrack(i);
2524 seed->CookdEdx(0.02,0.6);
2525 CookLabel(seed,0.1); //For comparison only
2526 esd->SetTPCClusterMap(seed->GetClusterMap());
2527 esd->SetTPCSharedMap(seed->GetSharedMap());
2529 if (AliTPCReconstructor::StreamLevel()>1 && seed!=0&&esd!=0) {
2530 TTreeSRedirector &cstream = *fDebugStreamer;
2537 if (seed->GetNumberOfClusters()>15){
2538 esd->UpdateTrackParams(seed,AliESDtrack::kTPCrefit);
2539 esd->SetTPCPoints(seed->GetPoints());
2540 esd->SetTPCPointsF(seed->GetNFoundable());
2541 Int_t ndedx = seed->GetNCDEDX(0)+seed->GetNCDEDX(1)+seed->GetNCDEDX(2)+seed->GetNCDEDX(3);
2542 Float_t sdedx = (seed->GetSDEDX(0)+seed->GetSDEDX(1)+seed->GetSDEDX(2)+seed->GetSDEDX(3))*0.25;
2543 Float_t dedx = seed->GetdEdx();
2544 esd->SetTPCsignal(dedx, sdedx, ndedx);
2546 // add seed to the esd track in Calib level
2548 Bool_t storeFriend = gRandom->Rndm()<(kMaxFriendTracks)/Float_t(nseed);
2549 if (AliTPCReconstructor::StreamLevel()>0 &&storeFriend){
2550 AliTPCseed * seedCopy = new AliTPCseed(*seed, kTRUE);
2551 esd->AddCalibObject(seedCopy);
2556 //printf("problem\n");
2559 //FindKinks(fSeeds,event);
2560 Info("RefitInward","Number of refitted tracks %d",ntracks);
2566 Int_t AliTPCtrackerMI::PropagateBack(AliESDEvent *event)
2569 // back propagation of ESD tracks
2575 PropagateBack(fSeeds);
2576 RemoveUsed2(fSeeds,0.4,0.4,20);
2577 //FindCurling(fSeeds, fEvent,1);
2578 FindSplitted(fSeeds, event,1); // find multi found tracks
2581 Int_t nseed = fSeeds->GetEntriesFast();
2583 for (Int_t i=0;i<nseed;i++){
2584 AliTPCseed * seed = (AliTPCseed*) fSeeds->UncheckedAt(i);
2585 if (!seed) continue;
2586 if (seed->GetKinkIndex(0)<0) UpdateKinkQualityM(seed); // update quality informations for kinks
2587 seed->UpdatePoints();
2588 AliESDtrack *esd=event->GetTrack(i);
2589 seed->CookdEdx(0.02,0.6);
2590 CookLabel(seed,0.1); //For comparison only
2591 if (seed->GetNumberOfClusters()>15){
2592 esd->UpdateTrackParams(seed,AliESDtrack::kTPCout);
2593 esd->SetTPCPoints(seed->GetPoints());
2594 esd->SetTPCPointsF(seed->GetNFoundable());
2595 Int_t ndedx = seed->GetNCDEDX(0)+seed->GetNCDEDX(1)+seed->GetNCDEDX(2)+seed->GetNCDEDX(3);
2596 Float_t sdedx = (seed->GetSDEDX(0)+seed->GetSDEDX(1)+seed->GetSDEDX(2)+seed->GetSDEDX(3))*0.25;
2597 Float_t dedx = seed->GetdEdx();
2598 esd->SetTPCsignal(dedx, sdedx, ndedx);
2600 Int_t eventnumber = event->GetEventNumberInFile();// patch 28 fev 06
2601 // This is most likely NOT the event number you'd like to use. It has nothing to do with the 'real' event number
2602 if (AliTPCReconstructor::StreamLevel()>1) {
2603 (*fDebugStreamer)<<"Cback"<<
2605 "EventNrInFile="<<eventnumber<<
2606 "\n"; // patch 28 fev 06
2610 //FindKinks(fSeeds,event);
2611 Info("PropagateBack","Number of back propagated tracks %d",ntracks);
2618 void AliTPCtrackerMI::DeleteSeeds()
2623 Int_t nseed = fSeeds->GetEntriesFast();
2624 for (Int_t i=0;i<nseed;i++){
2625 AliTPCseed * seed = (AliTPCseed*)fSeeds->At(i);
2626 if (seed) delete fSeeds->RemoveAt(i);
2633 void AliTPCtrackerMI::ReadSeeds(AliESDEvent *event, Int_t direction)
2636 //read seeds from the event
2638 Int_t nentr=event->GetNumberOfTracks();
2640 Info("ReadSeeds", "Number of ESD tracks: %d\n", nentr);
2645 fSeeds = new TObjArray(nentr);
2649 for (Int_t i=0; i<nentr; i++) {
2650 AliESDtrack *esd=event->GetTrack(i);
2651 ULong_t status=esd->GetStatus();
2652 if (!(status&AliESDtrack::kTPCin)) continue;
2653 AliTPCtrack t(*esd);
2654 t.SetNumberOfClusters(0);
2655 // AliTPCseed *seed = new AliTPCseed(t,t.GetAlpha());
2656 AliTPCseed *seed = new AliTPCseed(t/*,t.GetAlpha()*/);
2657 seed->SetUniqueID(esd->GetID());
2658 for (Int_t ikink=0;ikink<3;ikink++) {
2659 Int_t index = esd->GetKinkIndex(ikink);
2660 seed->GetKinkIndexes()[ikink] = index;
2661 if (index==0) continue;
2662 index = TMath::Abs(index);
2663 AliESDkink * kink = fEvent->GetKink(index-1);
2664 if (kink&&esd->GetKinkIndex(ikink)<0){
2665 if ((status & AliESDtrack::kTRDrefit) != 0) kink->SetStatus(1,2);
2666 if ((status & AliESDtrack::kITSout) != 0) kink->SetStatus(1,0);
2668 if (kink&&esd->GetKinkIndex(ikink)>0){
2669 if ((status & AliESDtrack::kTRDrefit) != 0) kink->SetStatus(1,6);
2670 if ((status & AliESDtrack::kITSout) != 0) kink->SetStatus(1,4);
2674 if (((status&AliESDtrack::kITSout)==0)&&(direction==1)) seed->ResetCovariance(10.);
2675 if ( direction ==2 &&(status & AliESDtrack::kTRDrefit) == 0 ) seed->ResetCovariance(10.);
2676 if ( direction ==2 && ((status & AliESDtrack::kTPCout) == 0) ) {
2681 if ( direction ==2 &&(status & AliESDtrack::kTRDrefit) > 0 ) {
2682 Double_t par0[5],par1[5],alpha,x;
2683 esd->GetInnerExternalParameters(alpha,x,par0);
2684 esd->GetExternalParameters(x,par1);
2685 Double_t delta1 = TMath::Abs(par0[4]-par1[4])/(0.000000001+TMath::Abs(par0[4]+par1[4]));
2686 Double_t delta2 = TMath::Abs(par0[3]-par1[3]);
2688 if (esd->GetTRDncls()>0) trdchi2 = esd->GetTRDchi2()/esd->GetTRDncls();
2689 //reset covariance if suspicious
2690 if ( (delta1>0.1) || (delta2>0.006) ||trdchi2>7.)
2691 seed->ResetCovariance(10.);
2696 // rotate to the local coordinate system
2698 fSectors=fInnerSec; fN=fkNIS;
2699 Double_t alpha=seed->GetAlpha() - fSectors->GetAlphaShift();
2700 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
2701 if (alpha < 0. ) alpha += 2.*TMath::Pi();
2702 Int_t ns=Int_t(alpha/fSectors->GetAlpha())%fN;
2703 alpha =ns*fSectors->GetAlpha() + fSectors->GetAlphaShift();
2704 if (alpha<-TMath::Pi()) alpha += 2*TMath::Pi();
2705 if (alpha>=TMath::Pi()) alpha -= 2*TMath::Pi();
2706 alpha-=seed->GetAlpha();
2707 if (!seed->Rotate(alpha)) {
2713 if (esd->GetKinkIndex(0)<=0){
2714 for (Int_t irow=0;irow<160;irow++){
2715 Int_t index = seed->GetClusterIndex2(irow);
2718 AliTPCclusterMI * cl = GetClusterMI(index);
2719 seed->SetClusterPointer(irow,cl);
2721 if ((index & 0x8000)==0){
2722 cl->Use(10); // accepted cluster
2724 cl->Use(6); // close cluster not accepted
2727 Info("ReadSeeds","Not found cluster");
2732 fSeeds->AddAt(seed,i);
2738 //_____________________________________________________________________________
2739 void AliTPCtrackerMI::MakeSeeds3(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, Float_t cuts[4],
2740 Float_t deltay, Int_t ddsec) {
2741 //-----------------------------------------------------------------
2742 // This function creates track seeds.
2743 // SEEDING WITH VERTEX CONSTRAIN
2744 //-----------------------------------------------------------------
2745 // cuts[0] - fP4 cut
2746 // cuts[1] - tan(phi) cut
2747 // cuts[2] - zvertex cut
2748 // cuts[3] - fP3 cut
2756 Double_t x[5], c[15];
2757 // Int_t di = i1-i2;
2759 AliTPCseed * seed = new AliTPCseed();
2760 Double_t alpha=fSectors->GetAlpha(), shift=fSectors->GetAlphaShift();
2761 Double_t cs=cos(alpha), sn=sin(alpha);
2763 // Double_t x1 =fOuterSec->GetX(i1);
2764 //Double_t xx2=fOuterSec->GetX(i2);
2766 Double_t x1 =GetXrow(i1);
2767 Double_t xx2=GetXrow(i2);
2769 Double_t x3=GetX(), y3=GetY(), z3=GetZ();
2771 Int_t imiddle = (i2+i1)/2; //middle pad row index
2772 Double_t xm = GetXrow(imiddle); // radius of middle pad-row
2773 const AliTPCtrackerRow& krm=GetRow(sec,imiddle); //middle pad -row
2777 const AliTPCtrackerRow& kr1=GetRow(ns,i1);
2778 Double_t ymax = GetMaxY(i1)-kr1.GetDeadZone()-1.5;
2779 Double_t ymaxm = GetMaxY(imiddle)-kr1.GetDeadZone()-1.5;
2782 // change cut on curvature if it can't reach this layer
2783 // maximal curvature set to reach it
2784 Double_t dvertexmax = TMath::Sqrt((x1-x3)*(x1-x3)+(ymax+5-y3)*(ymax+5-y3));
2785 if (dvertexmax*0.5*cuts[0]>0.85){
2786 cuts[0] = 0.85/(dvertexmax*0.5+1.);
2788 Double_t r2min = 1/(cuts[0]*cuts[0]); //minimal square of radius given by cut
2791 if (deltay>0) ddsec = 0;
2792 // loop over clusters
2793 for (Int_t is=0; is < kr1; is++) {
2795 if (kr1[is]->IsUsed(10)) continue;
2796 Double_t y1=kr1[is]->GetY(), z1=kr1[is]->GetZ();
2797 //if (TMath::Abs(y1)>ymax) continue;
2799 if (deltay>0 && TMath::Abs(ymax-TMath::Abs(y1))> deltay ) continue; // seed only at the edge
2801 // find possible directions
2802 Float_t anglez = (z1-z3)/(x1-x3);
2803 Float_t extraz = z1 - anglez*(x1-xx2); // extrapolated z
2806 //find rotation angles relative to line given by vertex and point 1
2807 Double_t dvertex2 = (x1-x3)*(x1-x3)+(y1-y3)*(y1-y3);
2808 Double_t dvertex = TMath::Sqrt(dvertex2);
2809 Double_t angle13 = TMath::ATan((y1-y3)/(x1-x3));
2810 Double_t cs13 = cos(-angle13), sn13 = sin(-angle13);
2813 // loop over 2 sectors
2819 Double_t dddz1=0; // direction of delta inclination in z axis
2826 for (Int_t dsec = dsec1; dsec<=dsec2;dsec++){
2827 Int_t sec2 = sec + dsec;
2829 // AliTPCtrackerRow& kr2 = fOuterSec[(sec2+fkNOS)%fkNOS][i2];
2830 //AliTPCtrackerRow& kr2m = fOuterSec[(sec2+fkNOS)%fkNOS][imiddle];
2831 AliTPCtrackerRow& kr2 = GetRow((sec2+fkNOS)%fkNOS,i2);
2832 AliTPCtrackerRow& kr2m = GetRow((sec2+fkNOS)%fkNOS,imiddle);
2833 Int_t index1 = TMath::Max(kr2.Find(extraz-0.6-dddz1*TMath::Abs(z1)*0.05)-1,0);
2834 Int_t index2 = TMath::Min(kr2.Find(extraz+0.6+dddz2*TMath::Abs(z1)*0.05)+1,kr2);
2836 // rotation angles to p1-p3
2837 Double_t cs13r = cos(-angle13+dsec*alpha)/dvertex, sn13r = sin(-angle13+dsec*alpha)/dvertex;
2838 Double_t x2, y2, z2;
2840 // Double_t dymax = maxangle*TMath::Abs(x1-xx2);
2843 Double_t dxx0 = (xx2-x3)*cs13r;
2844 Double_t dyy0 = (xx2-x3)*sn13r;
2845 for (Int_t js=index1; js < index2; js++) {
2846 const AliTPCclusterMI *kcl = kr2[js];
2847 if (kcl->IsUsed(10)) continue;
2849 //calcutate parameters
2851 Double_t yy0 = dyy0 +(kcl->GetY()-y3)*cs13r;
2853 if (TMath::Abs(yy0)<0.000001) continue;
2854 Double_t xx0 = dxx0 -(kcl->GetY()-y3)*sn13r;
2855 Double_t y0 = 0.5*(xx0*xx0+yy0*yy0-xx0)/yy0;
2856 Double_t r02 = (0.25+y0*y0)*dvertex2;
2857 //curvature (radius) cut
2858 if (r02<r2min) continue;
2862 Double_t c0 = 1/TMath::Sqrt(r02);
2866 //Double_t dfi0 = 2.*TMath::ASin(dvertex*c0*0.5);
2867 //Double_t dfi1 = 2.*TMath::ASin(TMath::Sqrt(yy0*yy0+(1-xx0)*(1-xx0))*dvertex*c0*0.5);
2868 Double_t dfi0 = 2.*AliTPCFastMath::FastAsin(dvertex*c0*0.5);
2869 Double_t dfi1 = 2.*AliTPCFastMath::FastAsin(TMath::Sqrt(yy0*yy0+(1-xx0)*(1-xx0))*dvertex*c0*0.5);
2872 Double_t z0 = kcl->GetZ();
2873 Double_t zzzz2 = z1-(z1-z3)*dfi1/dfi0;
2874 if (TMath::Abs(zzzz2-z0)>0.5) continue;
2877 Double_t dip = (z1-z0)*c0/dfi1;
2878 Double_t x0 = (0.5*cs13+y0*sn13)*dvertex*c0;
2889 x2= xx2*cs-y2*sn*dsec;
2890 y2=+xx2*sn*dsec+y2*cs;
2900 // do we have cluster at the middle ?
2902 GetProlongation(x1,xm,x,ym,zm);
2904 AliTPCclusterMI * cm=0;
2905 if (TMath::Abs(ym)-ymaxm<0){
2906 cm = krm.FindNearest2(ym,zm,1.0,0.6,dummy);
2907 if ((!cm) || (cm->IsUsed(10))) {
2912 // rotate y1 to system 0
2913 // get state vector in rotated system
2914 Double_t yr1 = (-0.5*sn13+y0*cs13)*dvertex*c0;
2915 Double_t xr2 = x0*cs+yr1*sn*dsec;
2916 Double_t xr[5]={kcl->GetY(),kcl->GetZ(), xr2, dip, c0};
2918 GetProlongation(xx2,xm,xr,ym,zm);
2919 if (TMath::Abs(ym)-ymaxm<0){
2920 cm = kr2m.FindNearest2(ym,zm,1.0,0.6,dummy);
2921 if ((!cm) || (cm->IsUsed(10))) {
2931 dym = ym - cm->GetY();
2932 dzm = zm - cm->GetZ();
2939 Double_t sy1=kr1[is]->GetSigmaY2()*2., sz1=kr1[is]->GetSigmaZ2()*2.;
2940 Double_t sy2=kcl->GetSigmaY2()*2., sz2=kcl->GetSigmaZ2()*2.;
2941 //Double_t sy3=400*3./12., sy=0.1, sz=0.1;
2942 Double_t sy3=25000*x[4]*x[4]+0.1, sy=0.1, sz=0.1;
2943 //Double_t sy3=25000*x[4]*x[4]*60+0.5, sy=0.1, sz=0.1;
2945 Double_t f40=(F1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
2946 Double_t f42=(F1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
2947 Double_t f43=(F1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
2948 Double_t f20=(F2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
2949 Double_t f22=(F2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
2950 Double_t f23=(F2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
2952 Double_t f30=(F3(x1,y1+sy,x2,y2,z1,z2)-x[3])/sy;
2953 Double_t f31=(F3(x1,y1,x2,y2,z1+sz,z2)-x[3])/sz;
2954 Double_t f32=(F3(x1,y1,x2,y2+sy,z1,z2)-x[3])/sy;
2955 Double_t f34=(F3(x1,y1,x2,y2,z1,z2+sz)-x[3])/sz;
2959 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
2960 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
2961 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
2962 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
2963 c[13]=f30*sy1*f40+f32*sy2*f42;
2964 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
2966 // if (!BuildSeed(kr1[is],kcl,0,x1,x2,x3,x,c)) continue;
2968 UInt_t index=kr1.GetIndex(is);
2969 seed->~AliTPCseed(); // this does not set the pointer to 0...
2970 AliTPCseed *track=new(seed) AliTPCseed(x1, ns*alpha+shift, x, c, index);
2972 track->SetIsSeeding(kTRUE);
2973 track->SetSeed1(i1);
2974 track->SetSeed2(i2);
2975 track->SetSeedType(3);
2979 FollowProlongation(*track, (i1+i2)/2,1);
2980 Int_t foundable,found,shared;
2981 track->GetClusterStatistic((i1+i2)/2,i1, found, foundable, shared, kTRUE);
2982 if ((found<0.55*foundable) || shared>0.5*found || (track->GetSigmaY2()+track->GetSigmaZ2())>0.5){
2984 seed->~AliTPCseed();
2990 FollowProlongation(*track, i2,1);
2994 track->SetBConstrain(1);
2995 // track->fLastPoint = i1+fInnerSec->GetNRows(); // first cluster in track position
2996 track->SetLastPoint(i1); // first cluster in track position
2997 track->SetFirstPoint(track->GetLastPoint());
2999 if (track->GetNumberOfClusters()<(i1-i2)*0.5 ||
3000 track->GetNumberOfClusters() < track->GetNFoundable()*0.6 ||
3001 track->GetNShared()>0.4*track->GetNumberOfClusters() ) {
3003 seed->~AliTPCseed();
3007 // Z VERTEX CONDITION
3008 Double_t zv, bz=GetBz();
3009 if ( !track->GetZAt(0.,bz,zv) ) continue;
3010 if (TMath::Abs(zv-z3)>cuts[2]) {
3011 FollowProlongation(*track, TMath::Max(i2-20,0));
3012 if ( !track->GetZAt(0.,bz,zv) ) continue;
3013 if (TMath::Abs(zv-z3)>cuts[2]){
3014 FollowProlongation(*track, TMath::Max(i2-40,0));
3015 if ( !track->GetZAt(0.,bz,zv) ) continue;
3016 if (TMath::Abs(zv-z3)>cuts[2] &&(track->GetNumberOfClusters() > track->GetNFoundable()*0.7)){
3017 // make seed without constrain
3018 AliTPCseed * track2 = MakeSeed(track,0.2,0.5,1.);
3019 FollowProlongation(*track2, i2,1);
3020 track2->SetBConstrain(kFALSE);
3021 track2->SetSeedType(1);
3022 arr->AddLast(track2);
3024 seed->~AliTPCseed();
3029 seed->~AliTPCseed();
3036 track->SetSeedType(0);
3037 arr->AddLast(track);
3038 seed = new AliTPCseed;
3040 // don't consider other combinations
3041 if (track->GetNumberOfClusters() > track->GetNFoundable()*0.8)
3047 Info("MakeSeeds3","\nSeeding statistic:\t%d\t%d\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin,nout1,nout2);
3053 void AliTPCtrackerMI::MakeSeeds5(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, Float_t cuts[4],
3058 //-----------------------------------------------------------------
3059 // This function creates track seeds.
3060 //-----------------------------------------------------------------
3061 // cuts[0] - fP4 cut
3062 // cuts[1] - tan(phi) cut
3063 // cuts[2] - zvertex cut
3064 // cuts[3] - fP3 cut
3074 Double_t x[5], c[15];
3076 // make temporary seed
3077 AliTPCseed * seed = new AliTPCseed;
3078 Double_t alpha=fOuterSec->GetAlpha(), shift=fOuterSec->GetAlphaShift();
3079 // Double_t cs=cos(alpha), sn=sin(alpha);
3084 Double_t x1 = GetXrow(i1-1);
3085 const AliTPCtrackerRow& kr1=GetRow(sec,i1-1);
3086 Double_t y1max = GetMaxY(i1-1)-kr1.GetDeadZone()-1.5;
3088 Double_t x1p = GetXrow(i1);
3089 const AliTPCtrackerRow& kr1p=GetRow(sec,i1);
3091 Double_t x1m = GetXrow(i1-2);
3092 const AliTPCtrackerRow& kr1m=GetRow(sec,i1-2);
3095 //last 3 padrow for seeding
3096 AliTPCtrackerRow& kr3 = GetRow((sec+fkNOS)%fkNOS,i1-7);
3097 Double_t x3 = GetXrow(i1-7);
3098 // Double_t y3max= GetMaxY(i1-7)-kr3.fDeadZone-1.5;
3100 AliTPCtrackerRow& kr3p = GetRow((sec+fkNOS)%fkNOS,i1-6);
3101 Double_t x3p = GetXrow(i1-6);
3103 AliTPCtrackerRow& kr3m = GetRow((sec+fkNOS)%fkNOS,i1-8);
3104 Double_t x3m = GetXrow(i1-8);
3109 Int_t im = i1-4; //middle pad row index
3110 Double_t xm = GetXrow(im); // radius of middle pad-row
3111 const AliTPCtrackerRow& krm=GetRow(sec,im); //middle pad -row
3112 // Double_t ymmax = GetMaxY(im)-kr1.fDeadZone-1.5;
3115 Double_t deltax = x1-x3;
3116 Double_t dymax = deltax*cuts[1];
3117 Double_t dzmax = deltax*cuts[3];
3119 // loop over clusters
3120 for (Int_t is=0; is < kr1; is++) {
3122 if (kr1[is]->IsUsed(10)) continue;
3123 Double_t y1=kr1[is]->GetY(), z1=kr1[is]->GetZ();
3125 if (deltay>0 && TMath::Abs(y1max-TMath::Abs(y1))> deltay ) continue; // seed only at the edge
3127 Int_t index1 = TMath::Max(kr3.Find(z1-dzmax)-1,0);
3128 Int_t index2 = TMath::Min(kr3.Find(z1+dzmax)+1,kr3);
3134 for (Int_t js=index1; js < index2; js++) {
3135 const AliTPCclusterMI *kcl = kr3[js];
3136 if (kcl->IsUsed(10)) continue;
3138 // apply angular cuts
3139 if (TMath::Abs(y1-y3)>dymax) continue;
3142 if (TMath::Abs(z1-z3)>dzmax) continue;
3144 Double_t angley = (y1-y3)/(x1-x3);
3145 Double_t anglez = (z1-z3)/(x1-x3);
3147 Double_t erry = TMath::Abs(angley)*(x1-x1m)*0.5+0.5;
3148 Double_t errz = TMath::Abs(anglez)*(x1-x1m)*0.5+0.5;
3150 Double_t yyym = angley*(xm-x1)+y1;
3151 Double_t zzzm = anglez*(xm-x1)+z1;
3153 const AliTPCclusterMI *kcm = krm.FindNearest2(yyym,zzzm,erry,errz,index);
3155 if (kcm->IsUsed(10)) continue;
3157 erry = TMath::Abs(angley)*(x1-x1m)*0.4+0.5;
3158 errz = TMath::Abs(anglez)*(x1-x1m)*0.4+0.5;
3165 // look around first
3166 const AliTPCclusterMI *kc1m = kr1m.FindNearest2(angley*(x1m-x1)+y1,
3172 if (kc1m->IsUsed(10)) used++;
3174 const AliTPCclusterMI *kc1p = kr1p.FindNearest2(angley*(x1p-x1)+y1,
3180 if (kc1p->IsUsed(10)) used++;
3182 if (used>1) continue;
3183 if (found<1) continue;
3187 const AliTPCclusterMI *kc3m = kr3m.FindNearest2(angley*(x3m-x3)+y3,
3193 if (kc3m->IsUsed(10)) used++;
3197 const AliTPCclusterMI *kc3p = kr3p.FindNearest2(angley*(x3p-x3)+y3,
3203 if (kc3p->IsUsed(10)) used++;
3207 if (used>1) continue;
3208 if (found<3) continue;
3218 x[4]=F1(x1,y1,x2,y2,x3,y3);
3219 //if (TMath::Abs(x[4]) >= cuts[0]) continue;
3222 x[2]=F2(x1,y1,x2,y2,x3,y3);
3225 x[3]=F3n(x1,y1,x2,y2,z1,z2,x[4]);
3226 //if (TMath::Abs(x[3]) > cuts[3]) continue;
3230 Double_t sy1=0.1, sz1=0.1;
3231 Double_t sy2=0.1, sz2=0.1;
3232 Double_t sy3=0.1, sy=0.1, sz=0.1;
3234 Double_t f40=(F1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
3235 Double_t f42=(F1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
3236 Double_t f43=(F1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
3237 Double_t f20=(F2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
3238 Double_t f22=(F2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
3239 Double_t f23=(F2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
3241 Double_t f30=(F3(x1,y1+sy,x2,y2,z1,z2)-x[3])/sy;
3242 Double_t f31=(F3(x1,y1,x2,y2,z1+sz,z2)-x[3])/sz;
3243 Double_t f32=(F3(x1,y1,x2,y2+sy,z1,z2)-x[3])/sy;
3244 Double_t f34=(F3(x1,y1,x2,y2,z1,z2+sz)-x[3])/sz;
3248 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3249 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3250 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3251 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3252 c[13]=f30*sy1*f40+f32*sy2*f42;
3253 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3255 // if (!BuildSeed(kr1[is],kcl,0,x1,x2,x3,x,c)) continue;
3257 UInt_t index=kr1.GetIndex(is);
3258 seed->~AliTPCseed();
3259 AliTPCseed *track=new(seed) AliTPCseed(x1, sec*alpha+shift, x, c, index);
3261 track->SetIsSeeding(kTRUE);
3264 FollowProlongation(*track, i1-7,1);
3265 if (track->GetNumberOfClusters() < track->GetNFoundable()*0.75 ||
3266 track->GetNShared()>0.6*track->GetNumberOfClusters() || ( track->GetSigmaY2()+ track->GetSigmaZ2())>0.6){
3268 seed->~AliTPCseed();
3274 FollowProlongation(*track, i2,1);
3275 track->SetBConstrain(0);
3276 track->SetLastPoint(i1+fInnerSec->GetNRows()); // first cluster in track position
3277 track->SetFirstPoint(track->GetLastPoint());
3279 if (track->GetNumberOfClusters()<(i1-i2)*0.5 ||
3280 track->GetNumberOfClusters()<track->GetNFoundable()*0.7 ||
3281 track->GetNShared()>2. || track->GetChi2()/track->GetNumberOfClusters()>6 || ( track->GetSigmaY2()+ track->GetSigmaZ2())>0.5 ) {
3283 seed->~AliTPCseed();
3288 FollowProlongation(*track, TMath::Max(i2-10,0),1);
3289 AliTPCseed * track2 = MakeSeed(track,0.2,0.5,0.9);
3290 FollowProlongation(*track2, i2,1);
3291 track2->SetBConstrain(kFALSE);
3292 track2->SetSeedType(4);
3293 arr->AddLast(track2);
3295 seed->~AliTPCseed();
3299 //arr->AddLast(track);
3300 //seed = new AliTPCseed;
3306 Info("MakeSeeds5","\nSeeding statiistic:\t%d\t%d\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin,nout1,nout2,nout3);
3312 //_____________________________________________________________________________
3313 void AliTPCtrackerMI::MakeSeeds2(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, Float_t */*cuts[4]*/,
3314 Float_t deltay, Bool_t /*bconstrain*/) {
3315 //-----------------------------------------------------------------
3316 // This function creates track seeds - without vertex constraint
3317 //-----------------------------------------------------------------
3318 // cuts[0] - fP4 cut - not applied
3319 // cuts[1] - tan(phi) cut
3320 // cuts[2] - zvertex cut - not applied
3321 // cuts[3] - fP3 cut
3331 Double_t alpha=fOuterSec->GetAlpha(), shift=fOuterSec->GetAlphaShift();
3332 // Double_t cs=cos(alpha), sn=sin(alpha);
3333 Int_t row0 = (i1+i2)/2;
3334 Int_t drow = (i1-i2)/2;
3335 const AliTPCtrackerRow& kr0=fSectors[sec][row0];
3336 AliTPCtrackerRow * kr=0;
3338 AliTPCpolyTrack polytrack;
3339 Int_t nclusters=fSectors[sec][row0];
3340 AliTPCseed * seed = new AliTPCseed;
3345 for (Int_t is=0; is < nclusters; is++) { //LOOP over clusters
3347 Int_t nfoundable =0;
3348 for (Int_t iter =1; iter<2; iter++){ //iterations
3349 const AliTPCtrackerRow& krm=fSectors[sec][row0-iter];
3350 const AliTPCtrackerRow& krp=fSectors[sec][row0+iter];
3351 const AliTPCclusterMI * cl= kr0[is];
3353 if (cl->IsUsed(10)) {
3359 Double_t x = kr0.GetX();
3360 // Initialization of the polytrack
3365 Double_t y0= cl->GetY();
3366 Double_t z0= cl->GetZ();
3370 Double_t ymax = fSectors->GetMaxY(row0)-kr0.GetDeadZone()-1.5;
3371 if (deltay>0 && TMath::Abs(ymax-TMath::Abs(y0))> deltay ) continue; // seed only at the edge
3373 erry = (0.5)*cl->GetSigmaY2()/TMath::Sqrt(cl->GetQ())*6;
3374 errz = (0.5)*cl->GetSigmaZ2()/TMath::Sqrt(cl->GetQ())*6;
3375 polytrack.AddPoint(x,y0,z0,erry, errz);
3378 if (cl->IsUsed(10)) sumused++;
3381 Float_t roady = (5*TMath::Sqrt(cl->GetSigmaY2()+0.2)+1.)*iter;
3382 Float_t roadz = (5*TMath::Sqrt(cl->GetSigmaZ2()+0.2)+1.)*iter;
3385 AliTPCclusterMI * cl1 = krm.FindNearest(y0,z0,roady,roadz);
3386 if (cl1 && TMath::Abs(ymax-TMath::Abs(y0))) {
3387 erry = (0.5)*cl1->GetSigmaY2()/TMath::Sqrt(cl1->GetQ())*3;
3388 errz = (0.5)*cl1->GetSigmaZ2()/TMath::Sqrt(cl1->GetQ())*3;
3389 if (cl1->IsUsed(10)) sumused++;
3390 polytrack.AddPoint(x,cl1->GetY(),cl1->GetZ(),erry,errz);
3394 AliTPCclusterMI * cl2 = krp.FindNearest(y0,z0,roady,roadz);
3396 erry = (0.5)*cl2->GetSigmaY2()/TMath::Sqrt(cl2->GetQ())*3;
3397 errz = (0.5)*cl2->GetSigmaZ2()/TMath::Sqrt(cl2->GetQ())*3;
3398 if (cl2->IsUsed(10)) sumused++;
3399 polytrack.AddPoint(x,cl2->GetY(),cl2->GetZ(),erry,errz);
3402 if (sumused>0) continue;
3404 polytrack.UpdateParameters();
3410 nfoundable = polytrack.GetN();
3411 nfound = nfoundable;
3413 for (Int_t ddrow = iter+1; ddrow<drow;ddrow++){
3414 Float_t maxdist = 0.8*(1.+3./(ddrow));
3415 for (Int_t delta = -1;delta<=1;delta+=2){
3416 Int_t row = row0+ddrow*delta;
3417 kr = &(fSectors[sec][row]);
3418 Double_t xn = kr->GetX();
3419 Double_t ymax = fSectors->GetMaxY(row)-kr->GetDeadZone()-1.5;
3420 polytrack.GetFitPoint(xn,yn,zn);
3421 if (TMath::Abs(yn)>ymax) continue;
3423 AliTPCclusterMI * cln = kr->FindNearest(yn,zn,roady,roadz);
3425 Float_t dist = TMath::Sqrt( (yn-cln->GetY())*(yn-cln->GetY())+(zn-cln->GetZ())*(zn-cln->GetZ()));
3428 erry = (dist+0.3)*cln->GetSigmaY2()/TMath::Sqrt(cln->GetQ())*(1.+1./(ddrow));
3429 errz = (dist+0.3)*cln->GetSigmaZ2()/TMath::Sqrt(cln->GetQ())*(1.+1./(ddrow));
3430 if (cln->IsUsed(10)) {
3431 // printf("used\n");
3439 polytrack.AddPoint(xn,cln->GetY(),cln->GetZ(),erry, errz);
3444 if ( (sumused>3) || (sumused>0.5*nfound) || (nfound<0.6*nfoundable)) break;
3445 polytrack.UpdateParameters();
3448 if ( (sumused>3) || (sumused>0.5*nfound)) {
3449 //printf("sumused %d\n",sumused);
3454 polytrack.GetFitDerivation(kr0.GetX(),dy,dz);
3455 AliTPCpolyTrack track2;
3457 polytrack.Refit(track2,0.5+TMath::Abs(dy)*0.3,0.4+TMath::Abs(dz)*0.3);
3458 if (track2.GetN()<0.5*nfoundable) continue;
3461 if ((nfound>0.6*nfoundable) &&( nfoundable>0.4*(i1-i2))) {
3463 // test seed with and without constrain
3464 for (Int_t constrain=0; constrain<=0;constrain++){
3465 // add polytrack candidate
3467 Double_t x[5], c[15];
3468 Double_t x1,x2,x3,y1,y2,y3,z1,z2,z3;
3469 track2.GetBoundaries(x3,x1);
3471 track2.GetFitPoint(x1,y1,z1);
3472 track2.GetFitPoint(x2,y2,z2);
3473 track2.GetFitPoint(x3,y3,z3);
3475 //is track pointing to the vertex ?
3478 polytrack.GetFitPoint(x0,y0,z0);
3491 x[4]=F1(x1,y1,x2,y2,x3,y3);
3493 // if (TMath::Abs(x[4]) >= cuts[0]) continue; //
3494 x[2]=F2(x1,y1,x2,y2,x3,y3);
3496 //if (TMath::Abs(x[4]*x1-x[2]) >= cuts[1]) continue;
3497 //x[3]=F3(x1,y1,x2,y2,z1,z2);
3498 x[3]=F3n(x1,y1,x3,y3,z1,z3,x[4]);
3499 //if (TMath::Abs(x[3]) > cuts[3]) continue;
3502 Double_t sy =0.1, sz =0.1;
3503 Double_t sy1=0.02, sz1=0.02;
3504 Double_t sy2=0.02, sz2=0.02;
3508 sy3=25000*x[4]*x[4]+0.1, sy=0.1, sz=0.1;
3511 Double_t f40=(F1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
3512 Double_t f42=(F1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
3513 Double_t f43=(F1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
3514 Double_t f20=(F2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
3515 Double_t f22=(F2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
3516 Double_t f23=(F2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
3518 Double_t f30=(F3(x1,y1+sy,x3,y3,z1,z3)-x[3])/sy;
3519 Double_t f31=(F3(x1,y1,x3,y3,z1+sz,z3)-x[3])/sz;
3520 Double_t f32=(F3(x1,y1,x3,y3+sy,z1,z3)-x[3])/sy;
3521 Double_t f34=(F3(x1,y1,x3,y3,z1,z3+sz)-x[3])/sz;
3526 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3527 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3528 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3529 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3530 c[13]=f30*sy1*f40+f32*sy2*f42;
3531 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3533 //Int_t row1 = fSectors->GetRowNumber(x1);
3534 Int_t row1 = GetRowNumber(x1);
3538 seed->~AliTPCseed();
3539 AliTPCseed *track=new(seed) AliTPCseed(x1,sec*alpha+shift,x,c,index);
3540 track->SetIsSeeding(kTRUE);
3541 Int_t rc=FollowProlongation(*track, i2);
3542 if (constrain) track->SetBConstrain(1);
3544 track->SetBConstrain(0);
3545 track->SetLastPoint(row1+fInnerSec->GetNRows()); // first cluster in track position
3546 track->SetFirstPoint(track->GetLastPoint());
3548 if (rc==0 || track->GetNumberOfClusters()<(i1-i2)*0.5 ||
3549 track->GetNumberOfClusters() < track->GetNFoundable()*0.6 ||
3550 track->GetNShared()>0.4*track->GetNumberOfClusters()) {
3553 seed->~AliTPCseed();
3556 arr->AddLast(track);
3557 seed = new AliTPCseed;
3561 } // if accepted seed
3564 Info("MakeSeeds2","\nSeeding statiistic:\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin3);
3570 AliTPCseed *AliTPCtrackerMI::MakeSeed(AliTPCseed *track, Float_t r0, Float_t r1, Float_t r2)
3574 //reseed using track points
3575 Int_t p0 = int(r0*track->GetNumberOfClusters()); // point 0
3576 Int_t p1 = int(r1*track->GetNumberOfClusters());
3577 Int_t p2 = int(r2*track->GetNumberOfClusters()); // last point
3579 Double_t x0[3],x1[3],x2[3];
3580 for (Int_t i=0;i<3;i++){
3586 // find track position at given ratio of the length
3587 Int_t sec0=0, sec1=0, sec2=0;
3590 for (Int_t i=0;i<160;i++){
3591 if (track->GetClusterPointer(i)){
3593 AliTPCTrackerPoint *trpoint =track->GetTrackPoint(i);
3594 if ( (index<p0) || x0[0]<0 ){
3595 if (trpoint->GetX()>1){
3596 clindex = track->GetClusterIndex2(i);
3598 x0[0] = trpoint->GetX();
3599 x0[1] = trpoint->GetY();
3600 x0[2] = trpoint->GetZ();
3601 sec0 = ((clindex&0xff000000)>>24)%18;
3606 if ( (index<p1) &&(trpoint->GetX()>1)){
3607 clindex = track->GetClusterIndex2(i);
3609 x1[0] = trpoint->GetX();
3610 x1[1] = trpoint->GetY();
3611 x1[2] = trpoint->GetZ();
3612 sec1 = ((clindex&0xff000000)>>24)%18;
3615 if ( (index<p2) &&(trpoint->GetX()>1)){
3616 clindex = track->GetClusterIndex2(i);
3618 x2[0] = trpoint->GetX();
3619 x2[1] = trpoint->GetY();
3620 x2[2] = trpoint->GetZ();
3621 sec2 = ((clindex&0xff000000)>>24)%18;
3628 Double_t alpha, cs,sn, xx2,yy2;
3630 alpha = (sec1-sec2)*fSectors->GetAlpha();
3631 cs = TMath::Cos(alpha);
3632 sn = TMath::Sin(alpha);
3633 xx2= x1[0]*cs-x1[1]*sn;
3634 yy2= x1[0]*sn+x1[1]*cs;
3638 alpha = (sec0-sec2)*fSectors->GetAlpha();
3639 cs = TMath::Cos(alpha);
3640 sn = TMath::Sin(alpha);
3641 xx2= x0[0]*cs-x0[1]*sn;
3642 yy2= x0[0]*sn+x0[1]*cs;
3648 Double_t x[5],c[15];
3652 x[4]=F1(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]);
3653 // if (x[4]>1) return 0;
3654 x[2]=F2(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]);
3655 x[3]=F3n(x2[0],x2[1],x0[0],x0[1],x2[2],x0[2],x[4]);
3656 //if (TMath::Abs(x[3]) > 2.2) return 0;
3657 //if (TMath::Abs(x[2]) > 1.99) return 0;
3659 Double_t sy =0.1, sz =0.1;
3661 Double_t sy1=0.02+track->GetSigmaY2(), sz1=0.02+track->GetSigmaZ2();
3662 Double_t sy2=0.01+track->GetSigmaY2(), sz2=0.01+track->GetSigmaZ2();
3663 Double_t sy3=0.01+track->GetSigmaY2();
3665 Double_t f40=(F1(x2[0],x2[1]+sy,x1[0],x1[1],x0[0],x0[1])-x[4])/sy;
3666 Double_t f42=(F1(x2[0],x2[1],x1[0],x1[1]+sy,x0[0],x0[1])-x[4])/sy;
3667 Double_t f43=(F1(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]+sy)-x[4])/sy;
3668 Double_t f20=(F2(x2[0],x2[1]+sy,x1[0],x1[1],x0[0],x0[1])-x[2])/sy;
3669 Double_t f22=(F2(x2[0],x2[1],x1[0],x1[1]+sy,x0[0],x0[1])-x[2])/sy;
3670 Double_t f23=(F2(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]+sy)-x[2])/sy;
3672 Double_t f30=(F3(x2[0],x2[1]+sy,x0[0],x0[1],x2[2],x0[2])-x[3])/sy;
3673 Double_t f31=(F3(x2[0],x2[1],x0[0],x0[1],x2[2]+sz,x0[2])-x[3])/sz;
3674 Double_t f32=(F3(x2[0],x2[1],x0[0],x0[1]+sy,x2[2],x0[2])-x[3])/sy;
3675 Double_t f34=(F3(x2[0],x2[1],x0[0],x0[1],x2[2],x0[2]+sz)-x[3])/sz;
3680 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3681 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3682 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3683 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3684 c[13]=f30*sy1*f40+f32*sy2*f42;
3685 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3687 // Int_t row1 = fSectors->GetRowNumber(x2[0]);
3688 AliTPCseed *seed=new AliTPCseed(x2[0], sec2*fSectors->GetAlpha()+fSectors->GetAlphaShift(), x, c, 0);
3689 // Double_t y0,z0,y1,z1, y2,z2;
3690 //seed->GetProlongation(x0[0],y0,z0);
3691 // seed->GetProlongation(x1[0],y1,z1);
3692 //seed->GetProlongation(x2[0],y2,z2);
3694 seed->SetLastPoint(pp2);
3695 seed->SetFirstPoint(pp2);
3702 AliTPCseed *AliTPCtrackerMI::ReSeed(AliTPCseed *track, Float_t r0, Float_t r1, Float_t r2)
3706 //reseed using founded clusters
3708 // Find the number of clusters
3709 Int_t nclusters = 0;
3710 for (Int_t irow=0;irow<160;irow++){
3711 if (track->GetClusterIndex(irow)>0) nclusters++;
3715 ipos[0] = TMath::Max(int(r0*nclusters),0); // point 0 cluster
3716 ipos[1] = TMath::Min(int(r1*nclusters),nclusters-1); //
3717 ipos[2] = TMath::Min(int(r2*nclusters),nclusters-1); // last point
3721 Int_t row[3],sec[3]={0,0,0};
3723 // find track row position at given ratio of the length
3725 for (Int_t irow=0;irow<160;irow++){
3726 if (track->GetClusterIndex2(irow)<0) continue;
3728 for (Int_t ipoint=0;ipoint<3;ipoint++){
3729 if (index<=ipos[ipoint]) row[ipoint] = irow;
3733 //Get cluster and sector position
3734 for (Int_t ipoint=0;ipoint<3;ipoint++){
3735 Int_t clindex = track->GetClusterIndex2(row[ipoint]);
3736 AliTPCclusterMI * cl = GetClusterMI(clindex);
3739 // AliTPCclusterMI * cl = GetClusterMI(clindex);
3742 sec[ipoint] = ((clindex&0xff000000)>>24)%18;
3743 xyz[ipoint][0] = GetXrow(row[ipoint]);
3744 xyz[ipoint][1] = cl->GetY();
3745 xyz[ipoint][2] = cl->GetZ();
3749 // Calculate seed state vector and covariance matrix
3751 Double_t alpha, cs,sn, xx2,yy2;
3753 alpha = (sec[1]-sec[2])*fSectors->GetAlpha();
3754 cs = TMath::Cos(alpha);
3755 sn = TMath::Sin(alpha);
3756 xx2= xyz[1][0]*cs-xyz[1][1]*sn;
3757 yy2= xyz[1][0]*sn+xyz[1][1]*cs;
3761 alpha = (sec[0]-sec[2])*fSectors->GetAlpha();
3762 cs = TMath::Cos(alpha);
3763 sn = TMath::Sin(alpha);
3764 xx2= xyz[0][0]*cs-xyz[0][1]*sn;
3765 yy2= xyz[0][0]*sn+xyz[0][1]*cs;
3771 Double_t x[5],c[15];
3775 x[4]=F1(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
3776 x[2]=F2(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
3777 x[3]=F3n(xyz[2][0],xyz[2][1],xyz[0][0],xyz[0][1],xyz[2][2],xyz[0][2],x[4]);
3779 Double_t sy =0.1, sz =0.1;
3781 Double_t sy1=0.2, sz1=0.2;
3782 Double_t sy2=0.2, sz2=0.2;
3785 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;
3786 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;
3787 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;
3788 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;
3789 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;
3790 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;
3792 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;
3793 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;
3794 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;
3795 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;
3800 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3801 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3802 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3803 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3804 c[13]=f30*sy1*f40+f32*sy2*f42;
3805 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3807 // Int_t row1 = fSectors->GetRowNumber(xyz[2][0]);
3808 AliTPCseed *seed=new AliTPCseed(xyz[2][0], sec[2]*fSectors->GetAlpha()+fSectors->GetAlphaShift(), x, c, 0);
3809 seed->SetLastPoint(row[2]);
3810 seed->SetFirstPoint(row[2]);
3815 AliTPCseed *AliTPCtrackerMI::ReSeed(AliTPCseed *track,Int_t r0, Bool_t forward)
3819 //reseed using founded clusters
3822 Int_t row[3]={0,0,0};
3823 Int_t sec[3]={0,0,0};
3825 // forward direction
3827 for (Int_t irow=r0;irow<160;irow++){
3828 if (track->GetClusterIndex(irow)>0){
3833 for (Int_t irow=160;irow>r0;irow--){
3834 if (track->GetClusterIndex(irow)>0){
3839 for (Int_t irow=row[2]-15;irow>row[0];irow--){
3840 if (track->GetClusterIndex(irow)>0){
3848 for (Int_t irow=0;irow<r0;irow++){
3849 if (track->GetClusterIndex(irow)>0){
3854 for (Int_t irow=r0;irow>0;irow--){
3855 if (track->GetClusterIndex(irow)>0){
3860 for (Int_t irow=row[2]-15;irow>row[0];irow--){
3861 if (track->GetClusterIndex(irow)>0){
3868 if ((row[2]-row[0])<20) return 0;
3869 if (row[1]==0) return 0;
3872 //Get cluster and sector position
3873 for (Int_t ipoint=0;ipoint<3;ipoint++){
3874 Int_t clindex = track->GetClusterIndex2(row[ipoint]);
3875 AliTPCclusterMI * cl = GetClusterMI(clindex);
3878 // AliTPCclusterMI * cl = GetClusterMI(clindex);
3881 sec[ipoint] = ((clindex&0xff000000)>>24)%18;
3882 xyz[ipoint][0] = GetXrow(row[ipoint]);
3883 AliTPCTrackerPoint * point = track->GetTrackPoint(row[ipoint]);
3884 if (point&&ipoint<2){
3886 xyz[ipoint][1] = point->GetY();
3887 xyz[ipoint][2] = point->GetZ();
3890 xyz[ipoint][1] = cl->GetY();
3891 xyz[ipoint][2] = cl->GetZ();
3898 // Calculate seed state vector and covariance matrix
3900 Double_t alpha, cs,sn, xx2,yy2;
3902 alpha = (sec[1]-sec[2])*fSectors->GetAlpha();
3903 cs = TMath::Cos(alpha);
3904 sn = TMath::Sin(alpha);
3905 xx2= xyz[1][0]*cs-xyz[1][1]*sn;
3906 yy2= xyz[1][0]*sn+xyz[1][1]*cs;
3910 alpha = (sec[0]-sec[2])*fSectors->GetAlpha();
3911 cs = TMath::Cos(alpha);
3912 sn = TMath::Sin(alpha);
3913 xx2= xyz[0][0]*cs-xyz[0][1]*sn;
3914 yy2= xyz[0][0]*sn+xyz[0][1]*cs;
3920 Double_t x[5],c[15];
3924 x[4]=F1(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
3925 x[2]=F2(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
3926 x[3]=F3n(xyz[2][0],xyz[2][1],xyz[0][0],xyz[0][1],xyz[2][2],xyz[0][2],x[4]);
3928 Double_t sy =0.1, sz =0.1;
3930 Double_t sy1=0.2, sz1=0.2;
3931 Double_t sy2=0.2, sz2=0.2;
3934 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;
3935 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;
3936 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;
3937 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;
3938 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;
3939 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;
3941 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;
3942 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;
3943 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;
3944 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;
3949 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3950 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3951 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3952 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3953 c[13]=f30*sy1*f40+f32*sy2*f42;
3954 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3956 // Int_t row1 = fSectors->GetRowNumber(xyz[2][0]);
3957 AliTPCseed *seed=new AliTPCseed(xyz[2][0], sec[2]*fSectors->GetAlpha()+fSectors->GetAlphaShift(), x, c, 0);
3958 seed->SetLastPoint(row[2]);
3959 seed->SetFirstPoint(row[2]);
3960 for (Int_t i=row[0];i<row[2];i++){
3961 seed->SetClusterIndex(i, track->GetClusterIndex(i));
3969 void AliTPCtrackerMI::FindMultiMC(TObjArray * array, AliESDEvent */*esd*/, Int_t iter)
3972 // find multi tracks - THIS FUNCTION IS ONLY FOR DEBUG PURPOSES
3974 // Use AliTPCReconstructor::StreamLevel()>2 if you want to tune parameters - cuts
3976 // Two reasons to have multiple find tracks
3977 // 1. Curling tracks can be find more than once
3978 // 2. Splitted tracks
3979 // a.) Multiple seeding to increase tracking efficiency - (~ 100% reached)
3980 // b.) Edge effect on the sector boundaries
3983 // Algorithm done in 2 phases - because of CPU consumption
3984 // it is n^2 algorithm - for lead-lead 20000x20000 combination are investigated
3986 // Algorihm for curling tracks sign:
3987 // 1 phase -makes a very rough fast cuts to minimize combinatorics
3988 // a.) opposite sign
3989 // b.) one of the tracks - not pointing to the primary vertex -
3990 // c.) delta tan(theta)
3992 // 2 phase - calculates DCA between tracks - time consument
3997 // General cuts - for splitted tracks and for curling tracks
3999 const Float_t kMaxdPhi = 0.2; // maximal distance in phi
4001 // Curling tracks cuts
4006 Int_t nentries = array->GetEntriesFast();
4007 AliHelix *helixes = new AliHelix[nentries];
4008 Float_t *xm = new Float_t[nentries];
4009 Float_t *dz0 = new Float_t[nentries];
4010 Float_t *dz1 = new Float_t[nentries];
4016 // Find track COG in x direction - point with best defined parameters
4018 for (Int_t i=0;i<nentries;i++){
4019 AliTPCseed* track = (AliTPCseed*)array->At(i);
4020 if (!track) continue;
4021 track->SetCircular(0);
4022 new (&helixes[i]) AliHelix(*track);
4026 track->GetDZ(GetX(),GetY(),GetZ(),GetBz(),dz);
4029 for (Int_t icl=0; icl<160; icl++){
4030 AliTPCclusterMI * cl = track->GetClusterPointer(icl);
4036 if (ncl>0) xm[i]/=Float_t(ncl);
4038 TTreeSRedirector &cstream = *fDebugStreamer;
4040 for (Int_t i0=0;i0<nentries;i0++){
4041 AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
4042 if (!track0) continue;
4043 Float_t xc0 = helixes[i0].GetHelix(6);
4044 Float_t yc0 = helixes[i0].GetHelix(7);
4045 Float_t r0 = helixes[i0].GetHelix(8);
4046 Float_t rc0 = TMath::Sqrt(xc0*xc0+yc0*yc0);
4047 Float_t fi0 = TMath::ATan2(yc0,xc0);
4049 for (Int_t i1=i0+1;i1<nentries;i1++){
4050 AliTPCseed * track1 = (AliTPCseed*)array->At(i1);
4051 if (!track1) continue;
4052 Int_t lab0=track0->GetLabel();
4053 Int_t lab1=track1->GetLabel();
4054 if (TMath::Abs(lab0)!=TMath::Abs(lab1)) continue;
4056 Float_t xc1 = helixes[i1].GetHelix(6);
4057 Float_t yc1 = helixes[i1].GetHelix(7);
4058 Float_t r1 = helixes[i1].GetHelix(8);
4059 Float_t rc1 = TMath::Sqrt(xc1*xc1+yc1*yc1);
4060 Float_t fi1 = TMath::ATan2(yc1,xc1);
4062 Float_t dfi = fi0-fi1;
4065 if (dfi>1.5*TMath::Pi()) dfi-=TMath::Pi(); // take care about edge effect
4066 if (dfi<-1.5*TMath::Pi()) dfi+=TMath::Pi(); //
4067 if (TMath::Abs(dfi)>kMaxdPhi&&helixes[i0].GetHelix(4)*helixes[i1].GetHelix(4)<0){
4069 // if short tracks with undefined sign
4070 fi1 = -TMath::ATan2(yc1,-xc1);
4073 Float_t dtheta = TMath::Abs(track0->GetTgl()-track1->GetTgl())<TMath::Abs(track0->GetTgl()+track1->GetTgl())? track0->GetTgl()-track1->GetTgl():track0->GetTgl()+track1->GetTgl();
4076 // debug stream to tune "fast cuts"
4078 Double_t dist[3]; // distance at X
4079 Double_t mdist[3]={0,0,0}; // mean distance X+-40cm
4080 track0->GetDistance(track1,0.5*(xm[i0]+xm[i1])-40.,dist,AliTracker::GetBz());
4081 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[i]);
4082 track0->GetDistance(track1,0.5*(xm[i0]+xm[i1])+40.,dist,AliTracker::GetBz());
4083 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[i]);
4084 track0->GetDistance(track1,0.5*(xm[i0]+xm[i1]),dist,AliTracker::GetBz());
4085 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[i]);
4086 for (Int_t i=0;i<3;i++) mdist[i]*=0.33333;
4090 for (Int_t icl=0; icl<160; icl++){
4091 AliTPCclusterMI * cl0 = track0->GetClusterPointer(icl);
4092 AliTPCclusterMI * cl1 = track1->GetClusterPointer(icl);
4095 if (cl0==cl1) sums++;
4103 "Tr0.="<<track0<< // seed0
4104 "Tr1.="<<track1<< // seed1
4105 "h0.="<<&helixes[i0]<<
4106 "h1.="<<&helixes[i1]<<
4108 "sum="<<sum<< //the sum of rows with cl in both
4109 "sums="<<sums<< //the sum of shared clusters
4110 "xm0="<<xm[i0]<< // the center of track
4111 "xm1="<<xm[i1]<< // the x center of track
4112 // General cut variables
4113 "dfi="<<dfi<< // distance in fi angle
4114 "dtheta="<<dtheta<< // distance int theta angle
4120 "dist0="<<dist[0]<< //distance x
4121 "dist1="<<dist[1]<< //distance y
4122 "dist2="<<dist[2]<< //distance z
4123 "mdist0="<<mdist[0]<< //distance x
4124 "mdist1="<<mdist[1]<< //distance y
4125 "mdist2="<<mdist[2]<< //distance z
4138 if (AliTPCReconstructor::StreamLevel()>1) {
4139 AliInfo("Time for curling tracks removal DEBUGGING MC");
4145 void AliTPCtrackerMI::FindSplitted(TObjArray * array, AliESDEvent */*esd*/, Int_t iter)
4149 // Two reasons to have multiple find tracks
4150 // 1. Curling tracks can be find more than once
4151 // 2. Splitted tracks
4152 // a.) Multiple seeding to increase tracking efficiency - (~ 100% reached)
4153 // b.) Edge effect on the sector boundaries
4155 // This function tries to find tracks closed in the parametric space
4157 // cut logic if distance is bigger than cut continue - Do Nothing
4158 const Float_t kMaxdTheta = 0.05; // maximal distance in theta
4159 const Float_t kMaxdPhi = 0.05; // maximal deistance in phi
4160 const Float_t kdelta = 40.; //delta r to calculate track distance
4162 // const Float_t kMaxDist0 = 1.; // maximal distance 0
4163 //const Float_t kMaxDist1 = 0.3; // maximal distance 1 - cut if track in separate rows
4166 TCut csec("csec","abs(Tr0.fRelativeSector-Tr1.fRelativeSector)<2");
4167 TCut cdtheta("cdtheta","abs(dtheta)<0.05");
4172 Int_t nentries = array->GetEntriesFast();
4173 AliHelix *helixes = new AliHelix[nentries];
4174 Float_t *xm = new Float_t[nentries];
4180 //Sort tracks according quality
4182 Int_t nseed = array->GetEntriesFast();
4183 Float_t * quality = new Float_t[nseed];
4184 Int_t * indexes = new Int_t[nseed];
4185 for (Int_t i=0; i<nseed; i++) {
4186 AliTPCseed *pt=(AliTPCseed*)array->UncheckedAt(i);
4191 pt->UpdatePoints(); //select first last max dens points
4192 Float_t * points = pt->GetPoints();
4193 if (points[3]<0.8) quality[i] =-1;
4194 quality[i] = (points[2]-points[0])+pt->GetNumberOfClusters();
4195 //prefer high momenta tracks if overlaps
4196 quality[i] *= TMath::Sqrt(TMath::Abs(pt->Pt())+0.5);
4198 TMath::Sort(nseed,quality,indexes);
4202 // Find track COG in x direction - point with best defined parameters
4204 for (Int_t i=0;i<nentries;i++){
4205 AliTPCseed* track = (AliTPCseed*)array->At(i);
4206 if (!track) continue;
4207 track->SetCircular(0);
4208 new (&helixes[i]) AliHelix(*track);
4211 for (Int_t icl=0; icl<160; icl++){
4212 AliTPCclusterMI * cl = track->GetClusterPointer(icl);
4218 if (ncl>0) xm[i]/=Float_t(ncl);
4220 TTreeSRedirector &cstream = *fDebugStreamer;
4222 for (Int_t is0=0;is0<nentries;is0++){
4223 Int_t i0 = indexes[is0];
4224 AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
4225 if (!track0) continue;
4226 if (track0->GetKinkIndexes()[0]!=0) continue;
4227 Float_t xc0 = helixes[i0].GetHelix(6);
4228 Float_t yc0 = helixes[i0].GetHelix(7);
4229 Float_t fi0 = TMath::ATan2(yc0,xc0);
4231 for (Int_t is1=is0+1;is1<nentries;is1++){
4232 Int_t i1 = indexes[is1];
4233 AliTPCseed * track1 = (AliTPCseed*)array->At(i1);
4234 if (!track1) continue;
4236 if (TMath::Abs(track0->GetRelativeSector()-track1->GetRelativeSector())>1) continue;
4237 if (track1->GetKinkIndexes()[0]>0 &&track0->GetKinkIndexes()[0]<0) continue;
4238 if (track1->GetKinkIndexes()[0]!=0) continue;
4240 Float_t dtheta = TMath::Abs(track0->GetTgl()-track1->GetTgl())<TMath::Abs(track0->GetTgl()+track1->GetTgl())? track0->GetTgl()-track1->GetTgl():track0->GetTgl()+track1->GetTgl();
4241 if (TMath::Abs(dtheta)>kMaxdTheta) continue;
4243 Float_t xc1 = helixes[i1].GetHelix(6);
4244 Float_t yc1 = helixes[i1].GetHelix(7);
4245 Float_t fi1 = TMath::ATan2(yc1,xc1);
4247 Float_t dfi = fi0-fi1;
4248 if (dfi>1.5*TMath::Pi()) dfi-=TMath::Pi(); // take care about edge effect
4249 if (dfi<-1.5*TMath::Pi()) dfi+=TMath::Pi(); //
4250 if (TMath::Abs(dfi)>kMaxdPhi&&helixes[i0].GetHelix(4)*helixes[i1].GetHelix(4)<0){
4252 // if short tracks with undefined sign
4253 fi1 = -TMath::ATan2(yc1,-xc1);
4256 if (TMath::Abs(dfi)>kMaxdPhi) continue;
4263 for (Int_t icl=0; icl<160; icl++){
4264 Int_t index0=track0->GetClusterIndex2(icl);
4265 Int_t index1=track1->GetClusterIndex2(icl);
4266 Bool_t used0 = (index0>0 && !(index0&0x8000));
4267 Bool_t used1 = (index1>0 && !(index1&0x8000));
4269 if (used0) sum0++; // used cluster0
4270 if (used1) sum1++; // used clusters1
4271 if (used0&&used1) sum++;
4272 if (index0==index1 && used0 && used1) sums++;
4276 if (sums<10) continue;
4277 if (sum<40) continue;
4278 if (sums/Float_t(TMath::Min(sum0,sum1))<0.5) continue;
4280 Double_t dist[5][4]; // distance at X
4281 Double_t mdist[4]={0,0,0,0}; // mean distance on range +-delta
4285 track0->GetDistance(track1,xm[i0],dist[0],AliTracker::GetBz());
4286 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[0][i]);
4287 track0->GetDistance(track1,xm[i1],dist[1],AliTracker::GetBz());
4288 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[1][i]);
4290 track0->GetDistance(track1,TMath::Min(xm[i1],xm[i0])-kdelta,dist[2],AliTracker::GetBz());
4291 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[2][i]);
4292 track0->GetDistance(track1,TMath::Max(xm[i1],xm[i0])+kdelta,dist[3],AliTracker::GetBz());
4293 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[3][i]);
4295 track0->GetDistance(track1,(xm[i1]+xm[i0])*0.5,dist[4],AliTracker::GetBz());
4296 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[4][i]);
4297 for (Int_t i=0;i<3;i++) mdist[i]*=0.2;
4300 Int_t lab0=track0->GetLabel();
4301 Int_t lab1=track1->GetLabel();
4302 cstream<<"Splitted"<<
4306 "Tr0.="<<track0<< // seed0
4307 "Tr1.="<<track1<< // seed1
4308 "h0.="<<&helixes[i0]<<
4309 "h1.="<<&helixes[i1]<<
4311 "sum="<<sum<< //the sum of rows with cl in both
4312 "sum0="<<sum0<< //the sum of rows with cl in 0 track
4313 "sum1="<<sum1<< //the sum of rows with cl in 1 track
4314 "sums="<<sums<< //the sum of shared clusters
4315 "xm0="<<xm[i0]<< // the center of track
4316 "xm1="<<xm[i1]<< // the x center of track
4317 // General cut variables
4318 "dfi="<<dfi<< // distance in fi angle
4319 "dtheta="<<dtheta<< // distance int theta angle
4322 "dist0="<<dist[4][0]<< //distance x
4323 "dist1="<<dist[4][1]<< //distance y
4324 "dist2="<<dist[4][2]<< //distance z
4325 "mdist0="<<mdist[0]<< //distance x
4326 "mdist1="<<mdist[1]<< //distance y
4327 "mdist2="<<mdist[2]<< //distance z
4330 delete array->RemoveAt(i1);
4335 AliInfo("Time for splitted tracks removal");
4341 void AliTPCtrackerMI::FindCurling(TObjArray * array, AliESDEvent */*esd*/, Int_t iter)
4344 // find Curling tracks
4345 // Use AliTPCReconstructor::StreamLevel()>1 if you want to tune parameters - cuts
4348 // Algorithm done in 2 phases - because of CPU consumption
4349 // it is n^2 algorithm - for lead-lead 20000x20000 combination are investigated
4350 // see detal in MC part what can be used to cut
4354 const Float_t kMaxC = 400; // maximal curvature to of the track
4355 const Float_t kMaxdTheta = 0.15; // maximal distance in theta
4356 const Float_t kMaxdPhi = 0.15; // maximal distance in phi
4357 const Float_t kPtRatio = 0.3; // ratio between pt
4358 const Float_t kMinDCAR = 2.; // distance to the primary vertex in r - see cpipe cut
4361 // Curling tracks cuts
4364 const Float_t kMaxDeltaRMax = 40; // distance in outer radius
4365 const Float_t kMaxDeltaRMin = 5.; // distance in lower radius - see cpipe cut
4366 const Float_t kMinAngle = 2.9; // angle between tracks
4367 const Float_t kMaxDist = 5; // biggest distance
4369 // The cuts can be tuned using the "MC information stored in Multi tree ==> see FindMultiMC
4372 TCut csign("csign","Tr0.fP[4]*Tr1.fP[4]<0"); //opposite sign
4373 TCut cmax("cmax","abs(Tr0.GetC())>1/400");
4374 TCut cda("cda","sqrt(dtheta^2+dfi^2)<0.15");
4375 TCut ccratio("ccratio","abs((Tr0.fP[4]+Tr1.fP[4])/(abs(Tr0.fP[4])+abs(Tr1.fP[4])))<0.3");
4376 TCut cpipe("cpipe", "min(abs(r0-rc0),abs(r1-rc1))>5");
4378 TCut cdrmax("cdrmax","abs(abs(rc0+r0)-abs(rc1+r1))<40")
4379 TCut cdrmin("cdrmin","abs(abs(rc0+r0)-abs(rc1+r1))<10")
4381 Multi->Draw("dfi","iter==0"+csign+cmax+cda+ccratio); ~94% of curling tracks fulfill
4382 Multi->Draw("min(abs(r0-rc0),abs(r1-rc1))","iter==0&&abs(lab1)==abs(lab0)"+csign+cmax+cda+ccratio+cpipe+cdrmin+cdrmax); //80%
4384 Curling2->Draw("dfi","iter==0&&abs(lab0)==abs(lab1)"+csign+cmax+cdtheta+cdfi+ccratio)
4390 Int_t nentries = array->GetEntriesFast();
4391 AliHelix *helixes = new AliHelix[nentries];
4392 for (Int_t i=0;i<nentries;i++){
4393 AliTPCseed* track = (AliTPCseed*)array->At(i);
4394 if (!track) continue;
4395 track->SetCircular(0);
4396 new (&helixes[i]) AliHelix(*track);
4402 Double_t phase[2][2],radius[2];
4406 TTreeSRedirector &cstream = *fDebugStreamer;
4408 for (Int_t i0=0;i0<nentries;i0++){
4409 AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
4410 if (!track0) continue;
4411 if (TMath::Abs(track0->GetC())<1/kMaxC) continue;
4412 Float_t xc0 = helixes[i0].GetHelix(6);
4413 Float_t yc0 = helixes[i0].GetHelix(7);
4414 Float_t r0 = helixes[i0].GetHelix(8);
4415 Float_t rc0 = TMath::Sqrt(xc0*xc0+yc0*yc0);
4416 Float_t fi0 = TMath::ATan2(yc0,xc0);
4418 for (Int_t i1=i0+1;i1<nentries;i1++){
4419 AliTPCseed * track1 = (AliTPCseed*)array->At(i1);
4420 if (!track1) continue;
4421 if (TMath::Abs(track1->GetC())<1/kMaxC) continue;
4422 Float_t xc1 = helixes[i1].GetHelix(6);
4423 Float_t yc1 = helixes[i1].GetHelix(7);
4424 Float_t r1 = helixes[i1].GetHelix(8);
4425 Float_t rc1 = TMath::Sqrt(xc1*xc1+yc1*yc1);
4426 Float_t fi1 = TMath::ATan2(yc1,xc1);
4428 Float_t dfi = fi0-fi1;
4431 if (dfi>1.5*TMath::Pi()) dfi-=TMath::Pi(); // take care about edge effect
4432 if (dfi<-1.5*TMath::Pi()) dfi+=TMath::Pi(); //
4433 Float_t dtheta = TMath::Abs(track0->GetTgl()-track1->GetTgl())<TMath::Abs(track0->GetTgl()+track1->GetTgl())? track0->GetTgl()-track1->GetTgl():track0->GetTgl()+track1->GetTgl();
4437 if (track0->GetBConstrain()&&track1->GetBConstrain()) continue; // not constrained
4438 if (track1->GetSigned1Pt()*track0->GetSigned1Pt()>0) continue; // not the same sign
4439 if ( TMath::Abs(track1->GetTgl()+track0->GetTgl())>kMaxdTheta) continue; //distance in the Theta
4440 if ( TMath::Abs(dfi)>kMaxdPhi) continue; //distance in phi
4441 if ( TMath::Sqrt(dfi*dfi+dtheta*dtheta)>kMaxdPhi) continue; //common angular offset
4443 Float_t pt0 = track0->GetSignedPt();
4444 Float_t pt1 = track1->GetSignedPt();
4445 if ((TMath::Abs(pt0+pt1)/(TMath::Abs(pt0)+TMath::Abs(pt1)))>kPtRatio) continue;
4446 if ((iter==1) && TMath::Abs(TMath::Abs(rc0+r0)-TMath::Abs(rc1+r1))>kMaxDeltaRMax) continue;
4447 if ((iter!=1) &&TMath::Abs(TMath::Abs(rc0-r0)-TMath::Abs(rc1-r1))>kMaxDeltaRMin) continue;
4448 if (TMath::Min(TMath::Abs(rc0-r0),TMath::Abs(rc1-r1))<kMinDCAR) continue;
4451 // Now find closest approach
4455 Int_t npoints = helixes[i0].GetRPHIintersections(helixes[i1], phase, radius,10);
4456 if (npoints==0) continue;
4457 helixes[i0].GetClosestPhases(helixes[i1], phase);
4461 Double_t hangles[3];
4462 helixes[i0].Evaluate(phase[0][0],xyz0);
4463 helixes[i1].Evaluate(phase[0][1],xyz1);
4465 helixes[i0].GetAngle(phase[0][0],helixes[i1],phase[0][1],hangles);
4466 Double_t deltah[2],deltabest;
4467 if (TMath::Abs(hangles[2])<kMinAngle) continue;
4471 helixes[i0].ParabolicDCA(helixes[i1],phase[0][0],phase[0][1],radius[0],deltah[0],2);
4473 helixes[i0].ParabolicDCA(helixes[i1],phase[1][0],phase[1][1],radius[1],deltah[1],2);
4474 if (deltah[1]<deltah[0]) ibest=1;
4476 deltabest = TMath::Sqrt(deltah[ibest]);
4477 helixes[i0].Evaluate(phase[ibest][0],xyz0);
4478 helixes[i1].Evaluate(phase[ibest][1],xyz1);
4479 helixes[i0].GetAngle(phase[ibest][0],helixes[i1],phase[ibest][1],hangles);
4480 Double_t radiusbest = TMath::Sqrt(radius[ibest]);
4482 if (deltabest>kMaxDist) continue;
4483 // if (mindcar+mindcaz<40 && (TMath::Abs(hangles[2])<kMinAngle ||deltabest>3)) continue;
4484 Bool_t sign =kFALSE;
4485 if (hangles[2]>kMinAngle) sign =kTRUE;
4488 // circular[i0] = kTRUE;
4489 // circular[i1] = kTRUE;
4490 if (track0->OneOverPt()<track1->OneOverPt()){
4491 track0->SetCircular(track0->GetCircular()+1);
4492 track1->SetCircular(track1->GetCircular()+2);
4495 track1->SetCircular(track1->GetCircular()+1);
4496 track0->SetCircular(track0->GetCircular()+2);
4499 if (AliTPCReconstructor::StreamLevel()>1){
4501 //debug stream to tune "fine" cuts
4502 Int_t lab0=track0->GetLabel();
4503 Int_t lab1=track1->GetLabel();
4504 cstream<<"Curling2"<<
4520 "npoints="<<npoints<<
4521 "hangles0="<<hangles[0]<<
4522 "hangles1="<<hangles[1]<<
4523 "hangles2="<<hangles[2]<<
4526 "radius="<<radiusbest<<
4527 "deltabest="<<deltabest<<
4528 "phase0="<<phase[ibest][0]<<
4529 "phase1="<<phase[ibest][1]<<
4537 if (AliTPCReconstructor::StreamLevel()>1) {
4538 AliInfo("Time for curling tracks removal");
4547 void AliTPCtrackerMI::FindKinks(TObjArray * array, AliESDEvent *esd)
4554 TObjArray *kinks= new TObjArray(10000);
4555 // TObjArray *v0s= new TObjArray(10000);
4556 Int_t nentries = array->GetEntriesFast();
4557 AliHelix *helixes = new AliHelix[nentries];
4558 Int_t *sign = new Int_t[nentries];
4559 Int_t *nclusters = new Int_t[nentries];
4560 Float_t *alpha = new Float_t[nentries];
4561 AliKink *kink = new AliKink();
4562 Int_t * usage = new Int_t[nentries];
4563 Float_t *zm = new Float_t[nentries];
4564 Float_t *z0 = new Float_t[nentries];
4565 Float_t *fim = new Float_t[nentries];
4566 Float_t *shared = new Float_t[nentries];
4567 Bool_t *circular = new Bool_t[nentries];
4568 Float_t *dca = new Float_t[nentries];
4569 //const AliESDVertex * primvertex = esd->GetVertex();
4571 // nentries = array->GetEntriesFast();
4576 for (Int_t i=0;i<nentries;i++){
4579 AliTPCseed* track = (AliTPCseed*)array->At(i);
4580 if (!track) continue;
4581 track->SetCircular(0);
4583 track->UpdatePoints();
4584 if (( track->GetPoints()[2]- track->GetPoints()[0])>5 && track->GetPoints()[3]>0.8){
4586 nclusters[i]=track->GetNumberOfClusters();
4587 alpha[i] = track->GetAlpha();
4588 new (&helixes[i]) AliHelix(*track);
4590 helixes[i].Evaluate(0,xyz);
4591 sign[i] = (track->GetC()>0) ? -1:1;
4594 if (track->GetProlongation(x,y,z)){
4596 fim[i] = alpha[i]+TMath::ATan2(y,x);
4599 zm[i] = track->GetZ();
4603 circular[i]= kFALSE;
4604 if (track->GetProlongation(0,y,z)) z0[i] = z;
4605 dca[i] = track->GetD(0,0);
4611 Int_t ncandidates =0;
4614 Double_t phase[2][2],radius[2];
4617 // Find circling track
4618 TTreeSRedirector &cstream = *fDebugStreamer;
4620 for (Int_t i0=0;i0<nentries;i0++){
4621 AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
4622 if (!track0) continue;
4623 if (track0->GetNumberOfClusters()<40) continue;
4624 if (TMath::Abs(1./track0->GetC())>200) continue;
4625 for (Int_t i1=i0+1;i1<nentries;i1++){
4626 AliTPCseed * track1 = (AliTPCseed*)array->At(i1);
4627 if (!track1) continue;
4628 if (track1->GetNumberOfClusters()<40) continue;
4629 if ( TMath::Abs(track1->GetTgl()+track0->GetTgl())>0.1) continue;
4630 if (track0->GetBConstrain()&&track1->GetBConstrain()) continue;
4631 if (TMath::Abs(1./track1->GetC())>200) continue;
4632 if (track1->GetSigned1Pt()*track0->GetSigned1Pt()>0) continue;
4633 if (track1->GetTgl()*track0->GetTgl()>0) continue;
4634 if (TMath::Max(TMath::Abs(1./track0->GetC()),TMath::Abs(1./track1->GetC()))>190) continue;
4635 if (track0->GetBConstrain()&&track1->OneOverPt()<track0->OneOverPt()) continue; //returning - lower momenta
4636 if (track1->GetBConstrain()&&track0->OneOverPt()<track1->OneOverPt()) continue; //returning - lower momenta
4638 Float_t mindcar = TMath::Min(TMath::Abs(dca[i0]),TMath::Abs(dca[i1]));
4639 if (mindcar<5) continue;
4640 Float_t mindcaz = TMath::Min(TMath::Abs(z0[i0]-GetZ()),TMath::Abs(z0[i1]-GetZ()));
4641 if (mindcaz<5) continue;
4642 if (mindcar+mindcaz<20) continue;
4645 Float_t xc0 = helixes[i0].GetHelix(6);
4646 Float_t yc0 = helixes[i0].GetHelix(7);
4647 Float_t r0 = helixes[i0].GetHelix(8);
4648 Float_t xc1 = helixes[i1].GetHelix(6);
4649 Float_t yc1 = helixes[i1].GetHelix(7);
4650 Float_t r1 = helixes[i1].GetHelix(8);
4652 Float_t rmean = (r0+r1)*0.5;
4653 Float_t delta =TMath::Sqrt((xc1-xc0)*(xc1-xc0)+(yc1-yc0)*(yc1-yc0));
4654 //if (delta>30) continue;
4655 if (delta>rmean*0.25) continue;
4656 if (TMath::Abs(r0-r1)/rmean>0.3) continue;
4658 Int_t npoints = helixes[i0].GetRPHIintersections(helixes[i1], phase, radius,10);
4659 if (npoints==0) continue;
4660 helixes[i0].GetClosestPhases(helixes[i1], phase);
4664 Double_t hangles[3];
4665 helixes[i0].Evaluate(phase[0][0],xyz0);
4666 helixes[i1].Evaluate(phase[0][1],xyz1);
4668 helixes[i0].GetAngle(phase[0][0],helixes[i1],phase[0][1],hangles);
4669 Double_t deltah[2],deltabest;
4670 if (hangles[2]<2.8) continue;
4673 helixes[i0].ParabolicDCA(helixes[i1],phase[0][0],phase[0][1],radius[0],deltah[0],2);
4675 helixes[i0].ParabolicDCA(helixes[i1],phase[1][0],phase[1][1],radius[1],deltah[1],2);
4676 if (deltah[1]<deltah[0]) ibest=1;
4678 deltabest = TMath::Sqrt(deltah[ibest]);
4679 helixes[i0].Evaluate(phase[ibest][0],xyz0);
4680 helixes[i1].Evaluate(phase[ibest][1],xyz1);
4681 helixes[i0].GetAngle(phase[ibest][0],helixes[i1],phase[ibest][1],hangles);
4682 Double_t radiusbest = TMath::Sqrt(radius[ibest]);
4684 if (deltabest>6) continue;
4685 if (mindcar+mindcaz<40 && (hangles[2]<3.12||deltabest>3)) continue;
4686 Bool_t sign =kFALSE;
4687 if (hangles[2]>3.06) sign =kTRUE;
4690 circular[i0] = kTRUE;
4691 circular[i1] = kTRUE;
4692 if (track0->OneOverPt()<track1->OneOverPt()){
4693 track0->SetCircular(track0->GetCircular()+1);
4694 track1->SetCircular(track1->GetCircular()+2);
4697 track1->SetCircular(track1->GetCircular()+1);
4698 track0->SetCircular(track0->GetCircular()+2);
4701 if (sign&&AliTPCReconstructor::StreamLevel()>1){
4703 Int_t lab0=track0->GetLabel();
4704 Int_t lab1=track1->GetLabel();
4705 cstream<<"Curling"<<
4712 "mindcar="<<mindcar<<
4713 "mindcaz="<<mindcaz<<
4716 "npoints="<<npoints<<
4717 "hangles0="<<hangles[0]<<
4718 "hangles2="<<hangles[2]<<
4723 "radius="<<radiusbest<<
4724 "deltabest="<<deltabest<<
4725 "phase0="<<phase[ibest][0]<<
4726 "phase1="<<phase[ibest][1]<<
4736 for (Int_t i =0;i<nentries;i++){
4737 if (sign[i]==0) continue;
4738 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
4745 Double_t cradius0 = 40*40;
4746 Double_t cradius1 = 270*270;
4749 Double_t cdist3=0.55;
4750 for (Int_t j =i+1;j<nentries;j++){
4752 if (sign[j]*sign[i]<1) continue;
4753 if ( (nclusters[i]+nclusters[j])>200) continue;
4754 if ( (nclusters[i]+nclusters[j])<80) continue;
4755 if ( TMath::Abs(zm[i]-zm[j])>60.) continue;
4756 if ( TMath::Abs(fim[i]-fim[j])>0.6 && TMath::Abs(fim[i]-fim[j])<5.7 ) continue;
4757 //AliTPCseed * track1 = (AliTPCseed*)array->At(j); Double_t phase[2][2],radius[2];
4758 Int_t npoints = helixes[i].GetRPHIintersections(helixes[j], phase, radius,20);
4759 if (npoints<1) continue;
4762 if (radius[0]<cradius0||radius[0]>cradius1) continue;
4765 if ( (radius[0]<cradius0||radius[0]>cradius1) && (radius[1]<cradius0||radius[1]>cradius1) ) continue;
4768 Double_t delta1=10000,delta2=10000;
4769 // cuts on the intersection radius
4770 helixes[i].LinearDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta1);
4771 if (radius[0]<20&&delta1<1) continue; //intersection at vertex
4772 if (radius[0]<10&&delta1<3) continue; //intersection at vertex
4774 helixes[i].LinearDCA(helixes[j],phase[1][0],phase[1][1],radius[1],delta2);
4775 if (radius[1]<20&&delta2<1) continue; //intersection at vertex
4776 if (radius[1]<10&&delta2<3) continue; //intersection at vertex
4779 Double_t distance1 = TMath::Min(delta1,delta2);
4780 if (distance1>cdist1) continue; // cut on DCA linear approximation
4782 npoints = helixes[i].GetRPHIintersections(helixes[j], phase, radius,20);
4783 helixes[i].ParabolicDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta1);
4784 if (radius[0]<20&&delta1<1) continue; //intersection at vertex
4785 if (radius[0]<10&&delta1<3) continue; //intersection at vertex
4788 helixes[i].ParabolicDCA(helixes[j],phase[1][0],phase[1][1],radius[1],delta2);
4789 if (radius[1]<20&&delta2<1) continue; //intersection at vertex
4790 if (radius[1]<10&&delta2<3) continue; //intersection at vertex
4792 distance1 = TMath::Min(delta1,delta2);
4795 rkink = TMath::Sqrt(radius[0]);
4798 rkink = TMath::Sqrt(radius[1]);
4800 if (distance1>cdist2) continue;
4803 AliTPCseed * track1 = (AliTPCseed*)array->At(j);
4806 Int_t row0 = GetRowNumber(rkink);
4807 if (row0<10) continue;
4808 if (row0>150) continue;
4811 Float_t dens00=-1,dens01=-1;
4812 Float_t dens10=-1,dens11=-1;
4814 Int_t found,foundable,shared;
4815 track0->GetClusterStatistic(0,row0-5, found, foundable,shared,kFALSE);
4816 if (foundable>5) dens00 = Float_t(found)/Float_t(foundable);
4817 track0->GetClusterStatistic(row0+5,155, found, foundable,shared,kFALSE);
4818 if (foundable>5) dens01 = Float_t(found)/Float_t(foundable);
4820 track1->GetClusterStatistic(0,row0-5, found, foundable,shared,kFALSE);
4821 if (foundable>10) dens10 = Float_t(found)/Float_t(foundable);
4822 track1->GetClusterStatistic(row0+5,155, found, foundable,shared,kFALSE);
4823 if (foundable>10) dens11 = Float_t(found)/Float_t(foundable);
4825 if (dens00<dens10 && dens01<dens11) continue;
4826 if (dens00>dens10 && dens01>dens11) continue;
4827 if (TMath::Max(dens00,dens10)<0.1) continue;
4828 if (TMath::Max(dens01,dens11)<0.3) continue;
4830 if (TMath::Min(dens00,dens10)>0.6) continue;
4831 if (TMath::Min(dens01,dens11)>0.6) continue;
4834 AliTPCseed * ktrack0, *ktrack1;
4843 if (TMath::Abs(ktrack0->GetC())>5) continue; // cut on the curvature for mother particle
4844 AliExternalTrackParam paramm(*ktrack0);
4845 AliExternalTrackParam paramd(*ktrack1);
4846 if (row0>60&&ktrack1->GetReference().GetX()>90.)new (¶md) AliExternalTrackParam(ktrack1->GetReference());
4849 kink->SetMother(paramm);
4850 kink->SetDaughter(paramd);
4853 Float_t x[3] = { kink->GetPosition()[0],kink->GetPosition()[1],kink->GetPosition()[2]};
4855 fParam->Transform0to1(x,index);
4856 fParam->Transform1to2(x,index);
4857 row0 = GetRowNumber(x[0]);
4859 if (kink->GetR()<100) continue;
4860 if (kink->GetR()>240) continue;
4861 if (kink->GetPosition()[2]/kink->GetR()>AliTPCReconstructor::GetCtgRange()) continue; //out of fiducial volume
4862 if (kink->GetDistance()>cdist3) continue;
4863 Float_t dird = kink->GetDaughterP()[0]*kink->GetPosition()[0]+kink->GetDaughterP()[1]*kink->GetPosition()[1]; // rough direction estimate
4864 if (dird<0) continue;
4866 Float_t dirm = kink->GetMotherP()[0]*kink->GetPosition()[0]+kink->GetMotherP()[1]*kink->GetPosition()[1]; // rough direction estimate
4867 if (dirm<0) continue;
4868 Float_t mpt = TMath::Sqrt(kink->GetMotherP()[0]*kink->GetMotherP()[0]+kink->GetMotherP()[1]*kink->GetMotherP()[1]);
4869 if (mpt<0.2) continue;
4872 //for high momenta momentum not defined well in first iteration
4873 Double_t qt = TMath::Sin(kink->GetAngle(2))*ktrack1->GetP();
4874 if (qt>0.35) continue;
4877 kink->SetLabel(CookLabel(ktrack0,0.4,0,row0),0);
4878 kink->SetLabel(CookLabel(ktrack1,0.4,row0,160),1);
4880 kink->SetTPCDensity(dens00,0,0);
4881 kink->SetTPCDensity(dens01,0,1);
4882 kink->SetTPCDensity(dens10,1,0);
4883 kink->SetTPCDensity(dens11,1,1);
4884 kink->SetIndex(i,0);
4885 kink->SetIndex(j,1);
4888 kink->SetTPCDensity(dens10,0,0);
4889 kink->SetTPCDensity(dens11,0,1);
4890 kink->SetTPCDensity(dens00,1,0);
4891 kink->SetTPCDensity(dens01,1,1);
4892 kink->SetIndex(j,0);
4893 kink->SetIndex(i,1);
4896 if (mpt<1||kink->GetAngle(2)>0.1){
4897 // angle and densities not defined yet
4898 if (kink->GetTPCDensityFactor()<0.8) continue;
4899 if ((2-kink->GetTPCDensityFactor())*kink->GetDistance() >0.25) continue;
4900 if (kink->GetAngle(2)*ktrack0->GetP()<0.003) continue; //too small angle
4901 if (kink->GetAngle(2)>0.2&&kink->GetTPCDensityFactor()<1.15) continue;
4902 if (kink->GetAngle(2)>0.2&&kink->GetTPCDensity(0,1)>0.05) continue;
4904 Float_t criticalangle = track0->GetSigmaSnp2()+track0->GetSigmaTgl2();
4905 criticalangle+= track1->GetSigmaSnp2()+track1->GetSigmaTgl2();
4906 criticalangle= 3*TMath::Sqrt(criticalangle);
4907 if (criticalangle>0.02) criticalangle=0.02;
4908 if (kink->GetAngle(2)<criticalangle) continue;
4911 Int_t drow = Int_t(2.+0.5/(0.05+kink->GetAngle(2))); // overlap region defined
4912 Float_t shapesum =0;
4914 for ( Int_t row = row0-drow; row<row0+drow;row++){
4915 if (row<0) continue;
4916 if (row>155) continue;
4917 if (ktrack0->GetClusterPointer(row)){
4918 AliTPCTrackerPoint *point =ktrack0->GetTrackPoint(row);
4919 shapesum+=point->GetSigmaY()+point->GetSigmaZ();
4922 if (ktrack1->GetClusterPointer(row)){
4923 AliTPCTrackerPoint *point =ktrack1->GetTrackPoint(row);
4924 shapesum+=point->GetSigmaY()+point->GetSigmaZ();
4929 kink->SetShapeFactor(-1.);
4932 kink->SetShapeFactor(shapesum/sum);
4934 // esd->AddKink(kink);
4935 kinks->AddLast(kink);
4941 // sort the kinks according quality - and refit them towards vertex
4943 Int_t nkinks = kinks->GetEntriesFast();
4944 Float_t *quality = new Float_t[nkinks];
4945 Int_t *indexes = new Int_t[nkinks];
4946 AliTPCseed *mothers = new AliTPCseed[nkinks];
4947 AliTPCseed *daughters = new AliTPCseed[nkinks];
4950 for (Int_t i=0;i<nkinks;i++){
4952 AliKink *kink = (AliKink*)kinks->At(i);
4954 // refit kinks towards vertex
4956 Int_t index0 = kink->GetIndex(0);
4957 Int_t index1 = kink->GetIndex(1);
4958 AliTPCseed * ktrack0 = (AliTPCseed*)array->At(index0);
4959 AliTPCseed * ktrack1 = (AliTPCseed*)array->At(index1);
4961 Int_t sumn=ktrack0->GetNumberOfClusters()+ktrack1->GetNumberOfClusters();
4963 // Refit Kink under if too small angle
4965 if (kink->GetAngle(2)<0.05){
4966 kink->SetTPCRow0(GetRowNumber(kink->GetR()));
4967 Int_t row0 = kink->GetTPCRow0();
4968 Int_t drow = Int_t(2.+0.5/(0.05+kink->GetAngle(2)));
4971 Int_t last = row0-drow;
4972 if (last<40) last=40;
4973 if (last<ktrack0->GetFirstPoint()+25) last = ktrack0->GetFirstPoint()+25;
4974 AliTPCseed* seed0 = ReSeed(ktrack0,last,kFALSE);
4977 Int_t first = row0+drow;
4978 if (first>130) first=130;
4979 if (first>ktrack1->GetLastPoint()-25) first = TMath::Max(ktrack1->GetLastPoint()-25,30);
4980 AliTPCseed* seed1 = ReSeed(ktrack1,first,kTRUE);
4982 if (seed0 && seed1){
4983 kink->SetStatus(1,8);
4984 if (RefitKink(*seed0,*seed1,*kink)) kink->SetStatus(1,9);
4985 row0 = GetRowNumber(kink->GetR());
4986 sumn = seed0->GetNumberOfClusters()+seed1->GetNumberOfClusters();
4987 mothers[i] = *seed0;
4988 daughters[i] = *seed1;
4991 delete kinks->RemoveAt(i);
4992 if (seed0) delete seed0;
4993 if (seed1) delete seed1;
4996 if (kink->GetDistance()>0.5 || kink->GetR()<110 || kink->GetR()>240) {
4997 delete kinks->RemoveAt(i);
4998 if (seed0) delete seed0;
4999 if (seed1) delete seed1;
5007 if (kink) quality[i] = 160*((0.1+kink->GetDistance())*(2.-kink->GetTPCDensityFactor()))/(sumn+40.); //the longest -clossest will win
5009 TMath::Sort(nkinks,quality,indexes,kFALSE);
5011 //remove double find kinks
5013 for (Int_t ikink0=1;ikink0<nkinks;ikink0++){
5014 AliKink * kink0 = (AliKink*) kinks->At(indexes[ikink0]);
5015 if (!kink0) continue;
5017 for (Int_t ikink1=0;ikink1<ikink0;ikink1++){
5018 if (!kink0) continue;
5019 AliKink * kink1 = (AliKink*) kinks->At(indexes[ikink1]);
5020 if (!kink1) continue;
5021 // if not close kink continue
5022 if (TMath::Abs(kink1->GetPosition()[2]-kink0->GetPosition()[2])>10) continue;
5023 if (TMath::Abs(kink1->GetPosition()[1]-kink0->GetPosition()[1])>10) continue;
5024 if (TMath::Abs(kink1->GetPosition()[0]-kink0->GetPosition()[0])>10) continue;
5026 AliTPCseed &mother0 = mothers[indexes[ikink0]];
5027 AliTPCseed &daughter0 = daughters[indexes[ikink0]];
5028 AliTPCseed &mother1 = mothers[indexes[ikink1]];
5029 AliTPCseed &daughter1 = daughters[indexes[ikink1]];
5030 Int_t row0 = (kink0->GetTPCRow0()+kink1->GetTPCRow0())/2;
5039 for (Int_t i=0;i<row0;i++){
5040 if (mother0.GetClusterIndex(i)>0 && mother1.GetClusterIndex(i)>0){
5043 if (mother0.GetClusterIndex(i)==mother1.GetClusterIndex(i)){
5050 for (Int_t i=row0;i<158;i++){
5051 if (daughter0.GetClusterIndex(i)>0 && daughter0.GetClusterIndex(i)>0){
5054 if (mother0.GetClusterIndex(i)==mother1.GetClusterIndex(i)){
5060 Float_t ratio = Float_t(same+1)/Float_t(both+1);
5061 Float_t ratiom = Float_t(samem+1)/Float_t(bothm+1);
5062 Float_t ratiod = Float_t(samed+1)/Float_t(bothd+1);
5063 if (ratio>0.3 && ratiom>0.5 &&ratiod>0.5) {
5064 Int_t sum0 = mother0.GetNumberOfClusters()+daughter0.GetNumberOfClusters();
5065 Int_t sum1 = mother1.GetNumberOfClusters()+daughter1.GetNumberOfClusters();
5067 shared[kink0->GetIndex(0)]= kTRUE;
5068 shared[kink0->GetIndex(1)]= kTRUE;
5069 delete kinks->RemoveAt(indexes[ikink0]);
5072 shared[kink1->GetIndex(0)]= kTRUE;
5073 shared[kink1->GetIndex(1)]= kTRUE;
5074 delete kinks->RemoveAt(indexes[ikink1]);
5081 for (Int_t i=0;i<nkinks;i++){
5082 AliKink * kink = (AliKink*) kinks->At(indexes[i]);
5083 if (!kink) continue;
5084 kink->SetTPCRow0(GetRowNumber(kink->GetR()));
5085 Int_t index0 = kink->GetIndex(0);
5086 Int_t index1 = kink->GetIndex(1);
5087 if (circular[index0]||circular[index1]&&kink->GetDistance()>0.2) continue;
5088 kink->SetMultiple(usage[index0],0);
5089 kink->SetMultiple(usage[index1],1);
5090 if (kink->GetMultiple()[0]+kink->GetMultiple()[1]>2) continue;
5091 if (kink->GetMultiple()[0]+kink->GetMultiple()[1]>0 && quality[indexes[i]]>0.2) continue;
5092 if (kink->GetMultiple()[0]+kink->GetMultiple()[1]>0 && kink->GetDistance()>0.2) continue;
5093 if (circular[index0]||circular[index1]&&kink->GetDistance()>0.1) continue;
5095 AliTPCseed * ktrack0 = (AliTPCseed*)array->At(index0);
5096 AliTPCseed * ktrack1 = (AliTPCseed*)array->At(index1);
5097 if (!ktrack0 || !ktrack1) continue;
5098 Int_t index = esd->AddKink(kink);
5101 if ( ktrack0->GetKinkIndex(0)==0 && ktrack1->GetKinkIndex(0)==0) { //best kink
5102 if (mothers[indexes[i]].GetNumberOfClusters()>20 && daughters[indexes[i]].GetNumberOfClusters()>20 && (mothers[indexes[i]].GetNumberOfClusters()+daughters[indexes[i]].GetNumberOfClusters())>100){
5103 *ktrack0 = mothers[indexes[i]];
5104 *ktrack1 = daughters[indexes[i]];
5108 ktrack0->SetKinkIndex(usage[index0],-(index+1));
5109 ktrack1->SetKinkIndex(usage[index1], (index+1));
5114 // Remove tracks corresponding to shared kink's
5116 for (Int_t i=0;i<nentries;i++){
5117 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5118 if (!track0) continue;
5119 if (track0->GetKinkIndex(0)!=0) continue;
5120 if (shared[i]) delete array->RemoveAt(i);
5125 RemoveUsed2(array,0.5,0.4,30);
5127 for (Int_t i=0;i<nentries;i++){
5128 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5129 if (!track0) continue;
5130 track0->CookdEdx(0.02,0.6);
5134 for (Int_t i=0;i<nentries;i++){
5135 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5136 if (!track0) continue;
5137 if (track0->Pt()<1.4) continue;
5138 //remove double high momenta tracks - overlapped with kink candidates
5141 for (Int_t icl=track0->GetFirstPoint();icl<track0->GetLastPoint(); icl++){
5142 if (track0->GetClusterPointer(icl)!=0){
5144 if (track0->GetClusterPointer(icl)->IsUsed(10)) shared++;
5147 if (Float_t(shared+1)/Float_t(all+1)>0.5) {
5148 delete array->RemoveAt(i);
5152 if (track0->GetKinkIndex(0)!=0) continue;
5153 if (track0->GetNumberOfClusters()<80) continue;
5155 AliTPCseed *pmother = new AliTPCseed();
5156 AliTPCseed *pdaughter = new AliTPCseed();
5157 AliKink *pkink = new AliKink;
5159 AliTPCseed & mother = *pmother;
5160 AliTPCseed & daughter = *pdaughter;
5161 AliKink & kink = *pkink;
5162 if (CheckKinkPoint(track0,mother,daughter, kink)){
5163 if (mother.GetNumberOfClusters()<30||daughter.GetNumberOfClusters()<20) {
5167 continue; //too short tracks
5169 if (mother.Pt()<1.4) {
5175 Int_t row0= kink.GetTPCRow0();
5176 if (kink.GetDistance()>0.5 || kink.GetR()<110. || kink.GetR()>240.) {
5183 Int_t index = esd->AddKink(&kink);
5184 mother.SetKinkIndex(0,-(index+1));
5185 daughter.SetKinkIndex(0,index+1);
5186 if (mother.GetNumberOfClusters()>50) {
5187 delete array->RemoveAt(i);
5188 array->AddAt(new AliTPCseed(mother),i);
5191 array->AddLast(new AliTPCseed(mother));
5193 array->AddLast(new AliTPCseed(daughter));
5194 for (Int_t icl=0;icl<row0;icl++) {
5195 if (mother.GetClusterPointer(icl)) mother.GetClusterPointer(icl)->Use(20);
5198 for (Int_t icl=row0;icl<158;icl++) {
5199 if (daughter.GetClusterPointer(icl)) daughter.GetClusterPointer(icl)->Use(20);
5208 delete [] daughters;
5230 printf("Ncandidates=\t%d\t%d\t%d\t%d\n",esd->GetNumberOfKinks(),ncandidates,ntracks,nall);
5234 void AliTPCtrackerMI::FindV0s(TObjArray * array, AliESDEvent *esd)
5240 TObjArray *tpcv0s = new TObjArray(100000);
5241 Int_t nentries = array->GetEntriesFast();
5242 AliHelix *helixes = new AliHelix[nentries];
5243 Int_t *sign = new Int_t[nentries];
5244 Float_t *alpha = new Float_t[nentries];
5245 Float_t *z0 = new Float_t[nentries];
5246 Float_t *dca = new Float_t[nentries];
5247 Float_t *sdcar = new Float_t[nentries];
5248 Float_t *cdcar = new Float_t[nentries];
5249 Float_t *pulldcar = new Float_t[nentries];
5250 Float_t *pulldcaz = new Float_t[nentries];
5251 Float_t *pulldca = new Float_t[nentries];
5252 Bool_t *isPrim = new Bool_t[nentries];
5253 const AliESDVertex * primvertex = esd->GetVertex();
5254 Double_t zvertex = primvertex->GetZv();
5256 // nentries = array->GetEntriesFast();
5258 for (Int_t i=0;i<nentries;i++){
5261 AliTPCseed* track = (AliTPCseed*)array->At(i);
5262 if (!track) continue;
5263 track->GetV0Indexes()[0] = 0; //rest v0 indexes
5264 track->GetV0Indexes()[1] = 0; //rest v0 indexes
5265 track->GetV0Indexes()[2] = 0; //rest v0 indexes
5267 alpha[i] = track->GetAlpha();
5268 new (&helixes[i]) AliHelix(*track);
5270 helixes[i].Evaluate(0,xyz);
5271 sign[i] = (track->GetC()>0) ? -1:1;
5275 if (track->GetProlongation(0,y,z)) z0[i] = z;
5276 dca[i] = track->GetD(0,0);
5278 // dca error parrameterezation + pulls
5280 sdcar[i] = TMath::Sqrt(0.150*0.150+(100*track->GetC())*(100*track->GetC()));
5281 if (TMath::Abs(track->GetTgl())>1) sdcar[i]*=2.5;
5282 cdcar[i] = TMath::Exp((TMath::Abs(track->GetC())-0.0106)*525.3);
5283 pulldcar[i] = (dca[i]-cdcar[i])/sdcar[i];
5284 pulldcaz[i] = (z0[i]-zvertex)/sdcar[i];
5285 pulldca[i] = TMath::Sqrt(pulldcar[i]*pulldcar[i]+pulldcaz[i]*pulldcaz[i]);
5286 if (track->TPCrPID(1)+track->TPCrPID(2)+track->TPCrPID(3)>0.5) {
5287 if (pulldca[i]<3.) isPrim[i]=kTRUE; //pion, muon and Kaon 3 sigma cut
5289 if (track->TPCrPID(4)>0.5) {
5290 if (pulldca[i]<0.5) isPrim[i]=kTRUE; //proton 0.5 sigma cut
5292 if (track->TPCrPID(0)>0.4) {
5293 isPrim[i]=kFALSE; //electron no sigma cut
5300 Int_t ncandidates =0;
5303 Double_t phase[2][2],radius[2];
5309 TTreeSRedirector &cstream = *fDebugStreamer;
5310 Float_t fprimvertex[3]={GetX(),GetY(),GetZ()};
5312 Double_t cradius0 = 10*10;
5313 Double_t cradius1 = 200*200;
5316 Double_t cpointAngle = 0.95;
5318 Double_t delta[2]={10000,10000};
5319 for (Int_t i =0;i<nentries;i++){
5320 if (sign[i]==0) continue;
5321 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5322 if (!track0) continue;
5323 if (AliTPCReconstructor::StreamLevel()>1){
5328 "zvertex="<<zvertex<<
5329 "sdcar0="<<sdcar[i]<<
5330 "cdcar0="<<cdcar[i]<<
5331 "pulldcar0="<<pulldcar[i]<<
5332 "pulldcaz0="<<pulldcaz[i]<<
5333 "pulldca0="<<pulldca[i]<<
5334 "isPrim="<<isPrim[i]<<
5338 if (track0->GetSigned1Pt()<0) continue;
5339 if (track0->GetKinkIndex(0)>0||isPrim[i]) continue; //daughter kink
5341 if (TMath::Abs(helixes[i].GetHelix(4))<0.000000001) continue;
5346 for (Int_t j =0;j<nentries;j++){
5347 AliTPCseed * track1 = (AliTPCseed*)array->At(j);
5348 if (!track1) continue;
5349 if (track1->GetKinkIndex(0)>0 || isPrim[j]) continue; //daughter kink
5350 if (sign[j]*sign[i]>0) continue;
5351 if (TMath::Abs(helixes[j].GetHelix(4))<0.000001) continue;
5352 if (track0->GetCircular()+track1->GetCircular()>1) continue; //circling -returning track
5355 // DCA to prim vertex cut
5361 Int_t npoints = helixes[i].GetRPHIintersections(helixes[j], phase, radius,cdist2);
5362 if (npoints<1) continue;
5366 if (radius[0]<cradius0||radius[0]>cradius1) continue;
5367 helixes[i].LinearDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta[0]);
5368 if (delta[0]>cdist1) continue;
5371 if (TMath::Max(radius[0],radius[1])<cradius0|| TMath::Min(radius[0],radius[1])>cradius1) continue;
5372 helixes[i].LinearDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta[0]);
5373 helixes[i].LinearDCA(helixes[j],phase[1][0],phase[1][1],radius[1],delta[1]);
5374 if (delta[1]<delta[0]) iclosest=1;
5375 if (delta[iclosest]>cdist1) continue;
5377 helixes[i].ParabolicDCA(helixes[j],phase[iclosest][0],phase[iclosest][1],radius[iclosest],delta[iclosest]);
5378 if (radius[iclosest]<cradius0 || radius[iclosest]>cradius1 || delta[iclosest]>cdist1) continue;
5380 Double_t pointAngle = helixes[i].GetPointAngle(helixes[j],phase[iclosest],fprimvertex);
5381 if (pointAngle<cpointAngle) continue;
5383 Bool_t isGamma = kFALSE;
5384 vertex.SetParamP(*track0); //track0 - plus
5385 vertex.SetParamN(*track1); //track1 - minus
5386 vertex.Update(fprimvertex);
5387 if (track0->TPCrPID(0)>0.3&&track1->TPCrPID(0)>0.3&&vertex.GetAnglep()[2]<0.15) isGamma=kTRUE; // gamma conversion candidate
5388 Double_t pointAngle2 = vertex.GetV0CosineOfPointingAngle();
5390 if (vertex.GetV0CosineOfPointingAngle()<cpointAngle && (!isGamma)) continue;// point angle cut
5391 //Bo: if (vertex.GetDist2()>2&&(!isGamma)) continue; // point angle cut
5392 if (vertex.GetDcaV0Daughters()>2&&(!isGamma)) continue;//Bo: // point angle cut
5393 Float_t sigmae = 0.15*0.15;
5394 if (vertex.GetRr()<80)
5395 sigmae += (sdcar[i]*sdcar[i]+sdcar[j]*sdcar[j])*(1.-vertex.GetRr()/80.)*(1.-vertex.GetRr()/80.);
5396 sigmae+= TMath::Sqrt(sigmae);
5397 //Bo: if (vertex.GetDist2()/sigmae>3.&&(!isGamma)) continue;
5398 if (vertex.GetDcaV0Daughters()/sigmae>3.&&(!isGamma)) continue;
5399 Float_t densb0=0,densb1=0,densa0=0,densa1=0;
5400 Int_t row0 = GetRowNumber(vertex.GetRr());
5402 //Bo: if (vertex.GetDist2()>0.2) continue;
5403 if (vertex.GetDcaV0Daughters()>0.2) continue;
5404 densb0 = track0->Density2(0,row0-5);
5405 densb1 = track1->Density2(0,row0-5);
5406 if (densb0>0.3|| densb1>0.3) continue; //clusters before vertex
5407 densa0 = track0->Density2(row0+5,row0+40);
5408 densa1 = track1->Density2(row0+5,row0+40);
5409 if ((densa0<0.4|| densa1<0.4)&&(!isGamma)) continue; //missing clusters after vertex
5412 densa0 = track0->Density2(0,40); //cluster density
5413 densa1 = track1->Density2(0,40); //cluster density
5414 if ((vertex.GetRr()<80&&densa0+densa1<1.)&&(!isGamma)) continue;
5416 //Bo: vertex.SetLab(0,track0->GetLabel());
5417 //Bo: vertex.SetLab(1,track1->GetLabel());
5418 vertex.SetChi2After((densa0+densa1)*0.5);
5419 vertex.SetChi2Before((densb0+densb1)*0.5);
5420 vertex.SetIndex(0,i);
5421 vertex.SetIndex(1,j);
5422 //Bo: vertex.SetStatus(1); // TPC v0 candidate
5423 vertex.SetOnFlyStatus(2);//Bo: // TPC v0 candidate
5424 //Bo: vertex.SetRp(track0->TPCrPIDs());
5425 //Bo: vertex.SetRm(track1->TPCrPIDs());
5426 tpcv0s->AddLast(new AliESDv0(vertex));
5429 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
5430 Double_t radiusm= (delta[0]<delta[1])? TMath::Sqrt(radius[0]):TMath::Sqrt(radius[1]);
5431 Double_t deltam= (delta[0]<delta[1])? TMath::Sqrt(delta[0]):TMath::Sqrt(delta[1]);
5432 if (AliTPCReconstructor::StreamLevel()>1) {
5433 Int_t lab0=track0->GetLabel();
5434 Int_t lab1=track1->GetLabel();
5435 Char_t c0=track0->GetCircular();
5436 Char_t c1=track1->GetCircular();
5439 "vertex.="<<&vertex<<
5442 "Helix0.="<<&helixes[i]<<
5445 "Helix1.="<<&helixes[j]<<
5446 "pointAngle="<<pointAngle<<
5447 "pointAngle2="<<pointAngle2<<
5452 "zvertex="<<zvertex<<
5455 "npoints="<<npoints<<
5456 "radius0="<<radius[0]<<
5457 "delta0="<<delta[0]<<
5458 "radius1="<<radius[1]<<
5459 "delta1="<<delta[1]<<
5460 "radiusm="<<radiusm<<
5462 "sdcar0="<<sdcar[i]<<
5463 "sdcar1="<<sdcar[j]<<
5464 "cdcar0="<<cdcar[i]<<
5465 "cdcar1="<<cdcar[j]<<
5466 "pulldcar0="<<pulldcar[i]<<
5467 "pulldcar1="<<pulldcar[j]<<
5468 "pulldcaz0="<<pulldcaz[i]<<
5469 "pulldcaz1="<<pulldcaz[j]<<
5470 "pulldca0="<<pulldca[i]<<
5471 "pulldca1="<<pulldca[j]<<
5481 Float_t *quality = new Float_t[ncandidates];
5482 Int_t *indexes = new Int_t[ncandidates];
5484 for (Int_t i=0;i<ncandidates;i++){
5486 AliESDv0 *v0 = (AliESDv0*)tpcv0s->At(i);
5487 quality[i] = 1./(1.00001-v0->GetV0CosineOfPointingAngle()); //base point angle
5488 // quality[i] /= (0.5+v0->GetDist2());
5489 // quality[i] *= v0->GetChi2After(); //density factor
5491 Int_t index0 = v0->GetIndex(0);
5492 Int_t index1 = v0->GetIndex(1);
5493 //Bo: Double_t minpulldca = TMath::Min(2.+pulldca[v0->GetIndex(0)],(2.+pulldca[v0->GetIndex(1)]) ); //pull
5494 Double_t minpulldca = TMath::Min(2.+pulldca[index0],(2.+pulldca[index1]) );//Bo:
5498 AliTPCseed * track0 = (AliTPCseed*)array->At(index0);
5499 AliTPCseed * track1 = (AliTPCseed*)array->At(index1);
5500 if (track0->TPCrPID(0)>0.3&&track1->TPCrPID(0)>0.3&&v0->GetAnglep()[2]<0.15) quality[i]+=1000000; // gamma conversion candidate
5501 if (track0->TPCrPID(4)>0.9||track1->TPCrPID(4)>0.9&&minpulldca>4) quality[i]*=10; // lambda candidate candidate
5504 TMath::Sort(ncandidates,quality,indexes,kTRUE);
5507 for (Int_t i=0;i<ncandidates;i++){
5508 AliESDv0 * v0 = (AliESDv0*)tpcv0s->At(indexes[i]);
5510 Int_t index0 = v0->GetIndex(0);
5511 Int_t index1 = v0->GetIndex(1);
5512 AliTPCseed * track0 = (AliTPCseed*)array->At(index0);
5513 AliTPCseed * track1 = (AliTPCseed*)array->At(index1);
5514 if (!track0||!track1) {
5518 Bool_t accept =kTRUE; //default accept
5519 Int_t *v0indexes0 = track0->GetV0Indexes();
5520 Int_t *v0indexes1 = track1->GetV0Indexes();
5522 Int_t order0 = (v0indexes0[0]!=0) ? 1:0;
5523 Int_t order1 = (v0indexes1[0]!=0) ? 1:0;
5524 if (v0indexes0[1]!=0) order0 =2;
5525 if (v0indexes1[1]!=0) order1 =2;
5527 if (v0indexes0[2]!=0) {order0=3; accept=kFALSE;}
5528 if (v0indexes0[2]!=0) {order1=3; accept=kFALSE;}
5530 AliESDv0 * v02 = v0;
5532 //Bo: v0->SetOrder(0,order0);
5533 //Bo: v0->SetOrder(1,order1);
5534 //Bo: v0->SetOrder(1,order0+order1);
5535 v0->SetOnFlyStatus(kTRUE);
5536 Int_t index = esd->AddV0(v0);
5537 v02 = esd->GetV0(index);
5538 v0indexes0[order0]=index;
5539 v0indexes1[order1]=index;
5543 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
5544 if (AliTPCReconstructor::StreamLevel()>1) {
5545 Int_t lab0=track0->GetLabel();
5546 Int_t lab1=track1->GetLabel();
5555 "dca0="<<dca[index0]<<
5556 "dca1="<<dca[index1]<<
5560 "quality="<<quality[i]<<
5561 "pulldca0="<<pulldca[index0]<<
5562 "pulldca1="<<pulldca[index1]<<
5586 printf("TPC V0 finder : naccepted\t%d\tncandidates\t%d\tntracks\t%d\tnall\t%d\n",naccepted,ncandidates,ntracks,nall);
5590 Int_t AliTPCtrackerMI::RefitKink(AliTPCseed &mother, AliTPCseed &daughter, AliESDkink &knk)
5593 // refit kink towards to the vertex
5596 AliKink &kink=(AliKink &)knk;
5598 Int_t row0 = GetRowNumber(kink.GetR());
5599 FollowProlongation(mother,0);
5600 mother.Reset(kFALSE);
5602 FollowProlongation(daughter,row0);
5603 daughter.Reset(kFALSE);
5604 FollowBackProlongation(daughter,158);
5605 daughter.Reset(kFALSE);
5606 Int_t first = TMath::Max(row0-20,30);
5607 Int_t last = TMath::Min(row0+20,140);
5609 const Int_t kNdiv =5;
5610 AliTPCseed param0[kNdiv]; // parameters along the track
5611 AliTPCseed param1[kNdiv]; // parameters along the track
5612 AliKink kinks[kNdiv]; // corresponding kink parameters
5615 for (Int_t irow=0; irow<kNdiv;irow++){
5616 rows[irow] = first +((last-first)*irow)/(kNdiv-1);
5618 // store parameters along the track
5620 for (Int_t irow=0;irow<kNdiv;irow++){
5621 FollowBackProlongation(mother, rows[irow]);
5622 FollowProlongation(daughter,rows[kNdiv-1-irow]);
5623 param0[irow] = mother;
5624 param1[kNdiv-1-irow] = daughter;
5628 for (Int_t irow=0; irow<kNdiv-1;irow++){
5629 if (param0[irow].GetNumberOfClusters()<kNdiv||param1[irow].GetNumberOfClusters()<kNdiv) continue;
5630 kinks[irow].SetMother(param0[irow]);
5631 kinks[irow].SetDaughter(param1[irow]);
5632 kinks[irow].Update();
5635 // choose kink with best "quality"
5637 Double_t mindist = 10000;
5638 for (Int_t irow=0;irow<kNdiv;irow++){
5639 if (param0[irow].GetNumberOfClusters()<20||param1[irow].GetNumberOfClusters()<20) continue;
5640 if (TMath::Abs(kinks[irow].GetR())>240.) continue;
5641 if (TMath::Abs(kinks[irow].GetR())<100.) continue;
5643 Float_t normdist = TMath::Abs(param0[irow].GetX()-kinks[irow].GetR())*(0.1+kink.GetDistance());
5644 normdist/= (param0[irow].GetNumberOfClusters()+param1[irow].GetNumberOfClusters()+40.);
5645 if (normdist < mindist){
5651 if (index==-1) return 0;
5654 param0[index].Reset(kTRUE);
5655 FollowProlongation(param0[index],0);
5657 mother = param0[index];
5658 daughter = param1[index]; // daughter in vertex
5660 kink.SetMother(mother);
5661 kink.SetDaughter(daughter);
5663 kink.SetTPCRow0(GetRowNumber(kink.GetR()));
5664 kink.SetTPCncls(param0[index].GetNumberOfClusters(),0);
5665 kink.SetTPCncls(param1[index].GetNumberOfClusters(),1);
5666 kink.SetLabel(CookLabel(&mother,0.4, 0,kink.GetTPCRow0()),0);
5667 kink.SetLabel(CookLabel(&daughter,0.4, kink.GetTPCRow0(),160),1);
5668 mother.SetLabel(kink.GetLabel(0));
5669 daughter.SetLabel(kink.GetLabel(1));
5675 void AliTPCtrackerMI::UpdateKinkQualityM(AliTPCseed * seed){
5677 // update Kink quality information for mother after back propagation
5679 if (seed->GetKinkIndex(0)>=0) return;
5680 for (Int_t ikink=0;ikink<3;ikink++){
5681 Int_t index = seed->GetKinkIndex(ikink);
5682 if (index>=0) break;
5683 index = TMath::Abs(index)-1;
5684 AliESDkink * kink = fEvent->GetKink(index);
5685 kink->SetTPCDensity(-1,0,0);
5686 kink->SetTPCDensity(1,0,1);
5688 Int_t row0 = kink->GetTPCRow0() - 2 - Int_t( 0.5/ (0.05+kink->GetAngle(2)));
5689 if (row0<15) row0=15;
5691 Int_t row1 = kink->GetTPCRow0() + 2 + Int_t( 0.5/ (0.05+kink->GetAngle(2)));
5692 if (row1>145) row1=145;
5694 Int_t found,foundable,shared;
5695 seed->GetClusterStatistic(0,row0, found, foundable,shared,kFALSE);
5696 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),0,0);
5697 seed->GetClusterStatistic(row1,155, found, foundable,shared,kFALSE);
5698 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),0,1);
5703 void AliTPCtrackerMI::UpdateKinkQualityD(AliTPCseed * seed){
5705 // update Kink quality information for daughter after refit
5707 if (seed->GetKinkIndex(0)<=0) return;
5708 for (Int_t ikink=0;ikink<3;ikink++){
5709 Int_t index = seed->GetKinkIndex(ikink);
5710 if (index<=0) break;
5711 index = TMath::Abs(index)-1;
5712 AliESDkink * kink = fEvent->GetKink(index);
5713 kink->SetTPCDensity(-1,1,0);
5714 kink->SetTPCDensity(-1,1,1);
5716 Int_t row0 = kink->GetTPCRow0() -2 - Int_t( 0.5/ (0.05+kink->GetAngle(2)));
5717 if (row0<15) row0=15;
5719 Int_t row1 = kink->GetTPCRow0() +2 + Int_t( 0.5/ (0.05+kink->GetAngle(2)));
5720 if (row1>145) row1=145;
5722 Int_t found,foundable,shared;
5723 seed->GetClusterStatistic(0,row0, found, foundable,shared,kFALSE);
5724 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),1,0);
5725 seed->GetClusterStatistic(row1,155, found, foundable,shared,kFALSE);
5726 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),1,1);
5732 Int_t AliTPCtrackerMI::CheckKinkPoint(AliTPCseed*seed,AliTPCseed &mother, AliTPCseed &daughter, AliESDkink &knk)
5735 // check kink point for given track
5736 // if return value=0 kink point not found
5737 // otherwise seed0 correspond to mother particle
5738 // seed1 correspond to daughter particle
5739 // kink parameter of kink point
5740 AliKink &kink=(AliKink &)knk;
5742 Int_t middlerow = (seed->GetFirstPoint()+seed->GetLastPoint())/2;
5743 Int_t first = seed->GetFirstPoint();
5744 Int_t last = seed->GetLastPoint();
5745 if (last-first<20) return 0; // shortest length - 2*30 = 60 pad-rows
5748 AliTPCseed *seed1 = ReSeed(seed,middlerow+20, kTRUE); //middle of chamber
5749 if (!seed1) return 0;
5750 FollowProlongation(*seed1,seed->GetLastPoint()-20);
5751 seed1->Reset(kTRUE);
5752 FollowProlongation(*seed1,158);
5753 seed1->Reset(kTRUE);
5754 last = seed1->GetLastPoint();
5756 AliTPCseed *seed0 = new AliTPCseed(*seed);
5757 seed0->Reset(kFALSE);
5760 AliTPCseed param0[20]; // parameters along the track
5761 AliTPCseed param1[20]; // parameters along the track
5762 AliKink kinks[20]; // corresponding kink parameters
5764 for (Int_t irow=0; irow<20;irow++){
5765 rows[irow] = first +((last-first)*irow)/19;
5767 // store parameters along the track
5769 for (Int_t irow=0;irow<20;irow++){
5770 FollowBackProlongation(*seed0, rows[irow]);
5771 FollowProlongation(*seed1,rows[19-irow]);
5772 param0[irow] = *seed0;
5773 param1[19-irow] = *seed1;
5777 for (Int_t irow=0; irow<19;irow++){
5778 kinks[irow].SetMother(param0[irow]);
5779 kinks[irow].SetDaughter(param1[irow]);
5780 kinks[irow].Update();
5783 // choose kink with biggest change of angle
5785 Double_t maxchange= 0;
5786 for (Int_t irow=1;irow<19;irow++){
5787 if (TMath::Abs(kinks[irow].GetR())>240.) continue;
5788 if (TMath::Abs(kinks[irow].GetR())<110.) continue;
5789 Float_t quality = TMath::Abs(kinks[irow].GetAngle(2))/(3.+TMath::Abs(kinks[irow].GetR()-param0[irow].GetX()));
5790 if ( quality > maxchange){
5791 maxchange = quality;
5798 if (index<0) return 0;
5800 Int_t row0 = GetRowNumber(kinks[index].GetR()); //row 0 estimate
5801 seed0 = new AliTPCseed(param0[index]);
5802 seed1 = new AliTPCseed(param1[index]);
5803 seed0->Reset(kFALSE);
5804 seed1->Reset(kFALSE);
5805 seed0->ResetCovariance(10.);
5806 seed1->ResetCovariance(10.);
5807 FollowProlongation(*seed0,0);
5808 FollowBackProlongation(*seed1,158);
5809 mother = *seed0; // backup mother at position 0
5810 seed0->Reset(kFALSE);
5811 seed1->Reset(kFALSE);
5812 seed0->ResetCovariance(10.);
5813 seed1->ResetCovariance(10.);
5815 first = TMath::Max(row0-20,0);
5816 last = TMath::Min(row0+20,158);
5818 for (Int_t irow=0; irow<20;irow++){
5819 rows[irow] = first +((last-first)*irow)/19;
5821 // store parameters along the track
5823 for (Int_t irow=0;irow<20;irow++){
5824 FollowBackProlongation(*seed0, rows[irow]);
5825 FollowProlongation(*seed1,rows[19-irow]);
5826 param0[irow] = *seed0;
5827 param1[19-irow] = *seed1;
5831 for (Int_t irow=0; irow<19;irow++){
5832 kinks[irow].SetMother(param0[irow]);
5833 kinks[irow].SetDaughter(param1[irow]);
5834 // param0[irow].Dump();
5835 //param1[irow].Dump();
5836 kinks[irow].Update();
5839 // choose kink with biggest change of angle
5842 for (Int_t irow=0;irow<20;irow++){
5843 if (TMath::Abs(kinks[irow].GetR())>250.) continue;
5844 if (TMath::Abs(kinks[irow].GetR())<90.) continue;
5845 Float_t quality = TMath::Abs(kinks[irow].GetAngle(2))/(3.+TMath::Abs(kinks[irow].GetR()-param0[irow].GetX()));
5846 if ( quality > maxchange){
5847 maxchange = quality;
5854 if (index==-1 || param0[index].GetNumberOfClusters()+param1[index].GetNumberOfClusters()<100){
5859 // Float_t anglesigma = TMath::Sqrt(param0[index].fC22+param0[index].fC33+param1[index].fC22+param1[index].fC33);
5861 kink.SetMother(param0[index]);
5862 kink.SetDaughter(param1[index]);
5864 row0 = GetRowNumber(kink.GetR());
5865 kink.SetTPCRow0(row0);
5866 kink.SetLabel(CookLabel(seed0,0.5,0,row0),0);
5867 kink.SetLabel(CookLabel(seed1,0.5,row0,158),1);
5868 kink.SetIndex(-10,0);
5869 kink.SetIndex(int(param0[index].GetNumberOfClusters()+param1[index].GetNumberOfClusters()),1);
5870 kink.SetTPCncls(param0[index].GetNumberOfClusters(),0);
5871 kink.SetTPCncls(param1[index].GetNumberOfClusters(),1);
5874 // new (&mother) AliTPCseed(param0[index]);
5875 daughter = param1[index];
5876 daughter.SetLabel(kink.GetLabel(1));
5877 param0[index].Reset(kTRUE);
5878 FollowProlongation(param0[index],0);
5879 mother = param0[index];
5880 mother.SetLabel(kink.GetLabel(0));
5890 AliTPCseed* AliTPCtrackerMI::ReSeed(AliTPCseed *t)
5893 // reseed - refit - track
5896 // Int_t last = fSectors->GetNRows()-1;
5898 if (fSectors == fOuterSec){
5899 first = TMath::Max(first, t->GetFirstPoint()-fInnerSec->GetNRows());
5903 first = t->GetFirstPoint();
5905 AliTPCseed * seed = MakeSeed(t,0.1,0.5,0.9);
5906 FollowBackProlongation(*t,fSectors->GetNRows()-1);
5908 FollowProlongation(*t,first);
5918 //_____________________________________________________________________________
5919 Int_t AliTPCtrackerMI::ReadSeeds(const TFile *inp) {
5920 //-----------------------------------------------------------------
5921 // This function reades track seeds.
5922 //-----------------------------------------------------------------
5923 TDirectory *savedir=gDirectory;
5925 TFile *in=(TFile*)inp;
5926 if (!in->IsOpen()) {
5927 cerr<<"AliTPCtrackerMI::ReadSeeds(): input file is not open !\n";
5932 TTree *seedTree=(TTree*)in->Get("Seeds");
5934 cerr<<"AliTPCtrackerMI::ReadSeeds(): ";
5935 cerr<<"can't get a tree with track seeds !\n";
5938 AliTPCtrack *seed=new AliTPCtrack;
5939 seedTree->SetBranchAddress("tracks",&seed);
5941 if (fSeeds==0) fSeeds=new TObjArray(15000);
5943 Int_t n=(Int_t)seedTree->GetEntries();
5944 for (Int_t i=0; i<n; i++) {
5945 seedTree->GetEvent(i);
5946 fSeeds->AddLast(new AliTPCseed(*seed/*,seed->GetAlpha()*/));
5955 Int_t AliTPCtrackerMI::Clusters2Tracks (AliESDEvent *esd)
5958 if (fSeeds) DeleteSeeds();
5961 if (!fSeeds) return 1;
5968 //_____________________________________________________________________________
5969 Int_t AliTPCtrackerMI::Clusters2Tracks() {
5970 //-----------------------------------------------------------------
5971 // This is a track finder.
5972 //-----------------------------------------------------------------
5973 TDirectory *savedir=gDirectory;
5977 fSeeds = Tracking();
5980 Info("Clusters2Tracks","Time for tracking: \t");timer.Print();timer.Start();
5982 //activate again some tracks
5983 for (Int_t i=0; i<fSeeds->GetEntriesFast(); i++) {
5984 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
5986 Int_t nc=t.GetNumberOfClusters();
5988 delete fSeeds->RemoveAt(i);
5992 if (pt->GetRemoval()==10) {
5993 if (pt->GetDensityFirst(20)>0.8 || pt->GetDensityFirst(30)>0.8 || pt->GetDensityFirst(40)>0.7)
5994 pt->Desactivate(10); // make track again active
5996 pt->Desactivate(20);
5997 delete fSeeds->RemoveAt(i);
6002 RemoveUsed2(fSeeds,0.85,0.85,0);
6003 if (AliTPCReconstructor::GetRecoParam()->GetDoKinks()) FindKinks(fSeeds,fEvent);
6004 //FindCurling(fSeeds, fEvent,0);
6005 if (AliTPCReconstructor::StreamLevel()>2) FindMultiMC(fSeeds, fEvent,0); // find multi found tracks
6006 RemoveUsed2(fSeeds,0.5,0.4,20);
6007 FindSplitted(fSeeds, fEvent,0); // find multi found tracks
6010 // // refit short tracks
6012 Int_t nseed=fSeeds->GetEntriesFast();
6015 for (Int_t i=0; i<nseed; i++) {
6016 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
6018 Int_t nc=t.GetNumberOfClusters();
6020 delete fSeeds->RemoveAt(i);
6023 CookLabel(pt,0.1); //For comparison only
6024 //if ((pt->IsActive() || (pt->fRemoval==10) )&& nc>50 &&pt->GetNumberOfClusters()>0.4*pt->fNFoundable){
6025 if ((pt->IsActive() || (pt->GetRemoval()==10) )){
6027 if (fDebug>0) cerr<<found<<'\r';
6031 delete fSeeds->RemoveAt(i);
6035 //RemoveOverlap(fSeeds,0.99,7,kTRUE);
6037 //RemoveUsed(fSeeds,0.9,0.9,6);
6039 nseed=fSeeds->GetEntriesFast();
6041 for (Int_t i=0; i<nseed; i++) {
6042 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
6044 Int_t nc=t.GetNumberOfClusters();
6046 delete fSeeds->RemoveAt(i);
6050 t.CookdEdx(0.02,0.6);
6051 // CheckKinkPoint(&t,0.05);
6052 //if ((pt->IsActive() || (pt->fRemoval==10) )&& nc>50 &&pt->GetNumberOfClusters()>0.4*pt->fNFoundable){
6053 if ((pt->IsActive() || (pt->GetRemoval()==10) )){
6061 delete fSeeds->RemoveAt(i);
6062 //AliTPCseed * seed1 = ReSeed(pt,0.05,0.5,1);
6064 // FollowProlongation(*seed1,0);
6065 // Int_t n = seed1->GetNumberOfClusters();
6066 // printf("fP4\t%f\t%f\n",seed1->GetC(),pt->GetC());
6067 // printf("fN\t%d\t%d\n", seed1->GetNumberOfClusters(),pt->GetNumberOfClusters());
6070 //AliTPCseed * seed2 = ReSeed(pt,0.95,0.5,0.05);
6074 SortTracks(fSeeds, 1);
6078 PrepareForBackProlongation(fSeeds,5.);
6079 PropagateBack(fSeeds);
6080 printf("Time for back propagation: \t");timer.Print();timer.Start();
6084 PrepareForProlongation(fSeeds,5.);
6085 PropagateForward2(fSeeds);
6087 printf("Time for FORWARD propagation: \t");timer.Print();timer.Start();
6088 // RemoveUsed(fSeeds,0.7,0.7,6);
6089 //RemoveOverlap(fSeeds,0.9,7,kTRUE);
6091 nseed=fSeeds->GetEntriesFast();
6093 for (Int_t i=0; i<nseed; i++) {
6094 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
6096 Int_t nc=t.GetNumberOfClusters();
6098 delete fSeeds->RemoveAt(i);
6101 t.CookdEdx(0.02,0.6);
6102 // CookLabel(pt,0.1); //For comparison only
6103 //if ((pt->IsActive() || (pt->fRemoval==10) )&& nc>50 &&pt->GetNumberOfClusters()>0.4*pt->fNFoundable){
6104 if ((pt->IsActive() || (pt->fRemoval==10) )){
6105 cerr<<found++<<'\r';
6108 delete fSeeds->RemoveAt(i);
6113 // fNTracks = found;
6115 Info("Clusters2Tracks","Time for overlap removal, track writing and dedx cooking: \t"); timer.Print();timer.Start();
6118 // cerr<<"Number of found tracks : "<<"\t"<<found<<endl;
6119 Info("Clusters2Tracks","Number of found tracks %d",found);
6121 // UnloadClusters();
6126 void AliTPCtrackerMI::Tracking(TObjArray * arr)
6129 // tracking of the seeds
6132 fSectors = fOuterSec;
6133 ParallelTracking(arr,150,63);
6134 fSectors = fOuterSec;
6135 ParallelTracking(arr,63,0);
6138 TObjArray * AliTPCtrackerMI::Tracking(Int_t seedtype, Int_t i1, Int_t i2, Float_t cuts[4], Float_t dy, Int_t dsec)
6143 TObjArray * arr = new TObjArray;
6145 fSectors = fOuterSec;
6148 for (Int_t sec=0;sec<fkNOS;sec++){
6149 if (seedtype==3) MakeSeeds3(arr,sec,i1,i2,cuts,dy, dsec);
6150 if (seedtype==4) MakeSeeds5(arr,sec,i1,i2,cuts,dy);
6151 if (seedtype==2) MakeSeeds2(arr,sec,i1,i2,cuts,dy);
6154 Info("Tracking","\nSeeding - %d\t%d\t%d\t%d\n",seedtype,i1,i2,arr->GetEntriesFast());
6166 TObjArray * AliTPCtrackerMI::Tracking()
6170 if (AliTPCReconstructor::GetRecoParam()->GetSpecialSeeding()) return TrackingSpecial();
6173 Int_t nup=fOuterSec->GetNRows()+fInnerSec->GetNRows();
6175 TObjArray * seeds = new TObjArray;
6184 Float_t fnumber = 3.0;
6185 Float_t fdensity = 3.0;
6190 for (Int_t delta = 0; delta<18; delta+=6){
6194 arr = Tracking(3,nup-1-delta,nup-1-delta-gap,cuts,-1,1);
6195 SumTracks(seeds,arr);
6196 SignClusters(seeds,fnumber,fdensity);
6198 for (Int_t i=2;i<6;i+=2){
6199 // seed high pt tracks
6202 arr = Tracking(3,nup-i-delta,nup-i-delta-gap,cuts,-1,0);
6203 SumTracks(seeds,arr);
6204 SignClusters(seeds,fnumber,fdensity);
6209 // RemoveUsed(seeds,0.9,0.9,1);
6210 // UnsignClusters();
6211 // SignClusters(seeds,fnumber,fdensity);
6215 for (Int_t delta = 20; delta<120; delta+=10){
6217 // seed high pt tracks
6221 arr = Tracking(3,nup-delta,nup-delta-gap,cuts,-1);
6222 SumTracks(seeds,arr);
6223 SignClusters(seeds,fnumber,fdensity);
6228 arr = Tracking(3,nup-delta-5,nup-delta-5-gap,cuts,-1);
6229 SumTracks(seeds,arr);
6230 SignClusters(seeds,fnumber,fdensity);
6241 Info("Tracking()","\n\nPrimary seeding\t%d\n\n",seeds->GetEntriesFast());
6245 // RemoveUsed(seeds,0.75,0.75,1);
6247 //SignClusters(seeds,fnumber,fdensity);
6256 arr = Tracking(4,nup-1,nup-1-gap,cuts,-1);
6257 SumTracks(seeds,arr);
6258 SignClusters(seeds,fnumber,fdensity);
6260 arr = Tracking(4,nup-2,nup-2-gap,cuts,-1);
6261 SumTracks(seeds,arr);
6262 SignClusters(seeds,fnumber,fdensity);
6264 arr = Tracking(4,nup-3,nup-3-gap,cuts,-1);
6265 SumTracks(seeds,arr);
6266 SignClusters(seeds,fnumber,fdensity);
6270 for (Int_t delta = 3; delta<30; delta+=5){
6276 arr = Tracking(4,nup-1-delta,nup-1-delta-gap,cuts,-1);
6277 SumTracks(seeds,arr);
6278 SignClusters(seeds,fnumber,fdensity);
6280 arr = Tracking(4,nup-3-delta,nup-5-delta-gap,cuts,4);
6281 SumTracks(seeds,arr);
6282 SignClusters(seeds,fnumber,fdensity);
6293 Int_t fLastSeedRowSec=AliTPCReconstructor::GetRecoParam()->GetLastSeedRowSec();
6296 for (Int_t delta = 30; delta<fLastSeedRowSec; delta+=10){
6302 arr = Tracking(4,nup-1-delta,nup-1-delta-gap,cuts,-1);
6303 SumTracks(seeds,arr);
6304 SignClusters(seeds,fnumber,fdensity);
6306 arr = Tracking(4,nup-5-delta,nup-5-delta-gap,cuts,5 );
6307 SumTracks(seeds,arr);
6308 SignClusters(seeds,fnumber,fdensity);
6312 Info("Tracking()","\n\nSecondary seeding\t%d\n\n",seeds->GetEntriesFast());
6323 TObjArray * AliTPCtrackerMI::TrackingSpecial()
6326 // seeding adjusted for laser and cosmic tests - short tracks with big inclination angle
6327 // no primary vertex seeding tried
6331 Int_t nup=fOuterSec->GetNRows()+fInnerSec->GetNRows();
6333 TObjArray * seeds = new TObjArray;
6338 Float_t fnumber = 3.0;
6339 Float_t fdensity = 3.0;
6342 cuts[0] = AliTPCReconstructor::GetRecoParam()->GetMaxC(); // max curvature
6343 cuts[1] = 3.5; // max tan(phi) angle for seeding
6344 cuts[2] = 3.; // not used (cut on z primary vertex)
6345 cuts[3] = 3.5; // max tan(theta) angle for seeding
6347 for (Int_t delta = 0; nup-delta-gap-1>0; delta+=3){
6349 arr = Tracking(4,nup-1-delta,nup-1-delta-gap,cuts,-1);
6350 SumTracks(seeds,arr);
6351 SignClusters(seeds,fnumber,fdensity);
6355 Info("Tracking()","\n\nSecondary seeding\t%d\n\n",seeds->GetEntriesFast());
6366 void AliTPCtrackerMI::SumTracks(TObjArray *arr1,TObjArray *&arr2) const
6369 //sum tracks to common container
6370 //remove suspicious tracks
6371 Int_t nseed = arr2->GetEntriesFast();
6372 for (Int_t i=0;i<nseed;i++){
6373 AliTPCseed *pt=(AliTPCseed*)arr2->UncheckedAt(i);
6376 // remove tracks with too big curvature
6378 if (TMath::Abs(pt->GetC())>AliTPCReconstructor::GetRecoParam()->GetMaxC()){
6379 delete arr2->RemoveAt(i);
6382 // REMOVE VERY SHORT TRACKS
6383 if (pt->GetNumberOfClusters()<20){
6384 delete arr2->RemoveAt(i);
6387 // NORMAL ACTIVE TRACK
6388 if (pt->IsActive()){
6389 arr1->AddLast(arr2->RemoveAt(i));
6392 //remove not usable tracks
6393 if (pt->GetRemoval()!=10){
6394 delete arr2->RemoveAt(i);
6398 // ENABLE ONLY ENOUGH GOOD STOPPED TRACKS
6399 if (pt->GetDensityFirst(20)>0.8 || pt->GetDensityFirst(30)>0.8 || pt->GetDensityFirst(40)>0.7)
6400 arr1->AddLast(arr2->RemoveAt(i));
6402 delete arr2->RemoveAt(i);
6406 delete arr2; arr2 = 0;
6411 void AliTPCtrackerMI::ParallelTracking(TObjArray * arr, Int_t rfirst, Int_t rlast)
6414 // try to track in parralel
6416 Int_t nseed=arr->GetEntriesFast();
6417 //prepare seeds for tracking
6418 for (Int_t i=0; i<nseed; i++) {
6419 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i), &t=*pt;
6421 if (!t.IsActive()) continue;
6422 // follow prolongation to the first layer
6423 if ( (fSectors ==fInnerSec) || (t.GetFirstPoint()-fParam->GetNRowLow()>rfirst+1) )
6424 FollowProlongation(t, rfirst+1);
6429 for (Int_t nr=rfirst; nr>=rlast; nr--){
6430 if (nr<fInnerSec->GetNRows())
6431 fSectors = fInnerSec;
6433 fSectors = fOuterSec;
6434 // make indexes with the cluster tracks for given
6436 // find nearest cluster
6437 for (Int_t i=0; i<nseed; i++) {
6438 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i), &t=*pt;
6440 if (nr==80) pt->UpdateReference();
6441 if (!pt->IsActive()) continue;
6442 // if ( (fSectors ==fOuterSec) && (pt->fFirstPoint-fParam->GetNRowLow())<nr) continue;
6443 if (pt->GetRelativeSector()>17) {
6446 UpdateClusters(t,nr);
6448 // prolonagate to the nearest cluster - if founded
6449 for (Int_t i=0; i<nseed; i++) {
6450 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
6452 if (!pt->IsActive()) continue;
6453 // if ((fSectors ==fOuterSec) && (pt->fFirstPoint-fParam->GetNRowLow())<nr) continue;
6454 if (pt->GetRelativeSector()>17) {
6457 FollowToNextCluster(*pt,nr);
6462 void AliTPCtrackerMI::PrepareForBackProlongation(TObjArray * arr,Float_t fac) const
6466 // if we use TPC track itself we have to "update" covariance
6468 Int_t nseed= arr->GetEntriesFast();
6469 for (Int_t i=0;i<nseed;i++){
6470 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
6474 //rotate to current local system at first accepted point
6475 Int_t index = pt->GetClusterIndex2(pt->GetFirstPoint());
6476 Int_t sec = (index&0xff000000)>>24;
6478 Float_t angle1 = fInnerSec->GetAlpha()*sec+fInnerSec->GetAlphaShift();
6479 if (angle1>TMath::Pi())
6480 angle1-=2.*TMath::Pi();
6481 Float_t angle2 = pt->GetAlpha();
6483 if (TMath::Abs(angle1-angle2)>0.001){
6484 pt->Rotate(angle1-angle2);
6485 //angle2 = pt->GetAlpha();
6486 //pt->fRelativeSector = pt->GetAlpha()/fInnerSec->GetAlpha();
6487 //if (pt->GetAlpha()<0)
6488 // pt->fRelativeSector+=18;
6489 //sec = pt->fRelativeSector;
6498 void AliTPCtrackerMI::PrepareForProlongation(TObjArray * arr, Float_t fac) const
6502 // if we use TPC track itself we have to "update" covariance
6504 Int_t nseed= arr->GetEntriesFast();
6505 for (Int_t i=0;i<nseed;i++){
6506 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
6509 pt->SetFirstPoint(pt->GetLastPoint());
6517 Int_t AliTPCtrackerMI::PropagateBack(TObjArray * arr)
6520 // make back propagation
6522 Int_t nseed= arr->GetEntriesFast();
6523 for (Int_t i=0;i<nseed;i++){
6524 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
6525 if (pt&& pt->GetKinkIndex(0)<=0) {
6526 //AliTPCseed *pt2 = new AliTPCseed(*pt);
6527 fSectors = fInnerSec;
6528 //FollowBackProlongation(*pt,fInnerSec->GetNRows()-1);
6529 //fSectors = fOuterSec;
6530 FollowBackProlongation(*pt,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1);
6531 //if (pt->GetNumberOfClusters()<(pt->fEsd->GetTPCclusters(0)) ){
6532 // Error("PropagateBack","Not prolonged track %d",pt->GetLabel());
6533 // FollowBackProlongation(*pt2,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1);
6536 if (pt&& pt->GetKinkIndex(0)>0) {
6537 AliESDkink * kink = fEvent->GetKink(pt->GetKinkIndex(0)-1);
6538 pt->SetFirstPoint(kink->GetTPCRow0());
6539 fSectors = fInnerSec;
6540 FollowBackProlongation(*pt,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1);
6548 Int_t AliTPCtrackerMI::PropagateForward2(TObjArray * arr)
6551 // make forward propagation
6553 Int_t nseed= arr->GetEntriesFast();
6555 for (Int_t i=0;i<nseed;i++){
6556 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
6558 FollowProlongation(*pt,0);
6567 Int_t AliTPCtrackerMI::PropagateForward()
6570 // propagate track forward
6572 Int_t nseed = fSeeds->GetEntriesFast();
6573 for (Int_t i=0;i<nseed;i++){
6574 AliTPCseed *pt = (AliTPCseed*)fSeeds->UncheckedAt(i);
6576 AliTPCseed &t = *pt;
6577 Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
6578 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
6579 if (alpha < 0. ) alpha += 2.*TMath::Pi();
6580 t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN);
6584 fSectors = fOuterSec;
6585 ParallelTracking(fSeeds,fOuterSec->GetNRows()+fInnerSec->GetNRows()-1,fInnerSec->GetNRows());
6586 fSectors = fInnerSec;
6587 ParallelTracking(fSeeds,fInnerSec->GetNRows()-1,0);
6596 Int_t AliTPCtrackerMI::PropagateBack(AliTPCseed * pt, Int_t row0, Int_t row1)
6599 // make back propagation, in between row0 and row1
6603 fSectors = fInnerSec;
6606 if (row1<fSectors->GetNRows())
6609 r1 = fSectors->GetNRows()-1;
6611 if (row0<fSectors->GetNRows()&& r1>0 )
6612 FollowBackProlongation(*pt,r1);
6613 if (row1<=fSectors->GetNRows())
6616 r1 = row1 - fSectors->GetNRows();
6617 if (r1<=0) return 0;
6618 if (r1>=fOuterSec->GetNRows()) return 0;
6619 fSectors = fOuterSec;
6620 return FollowBackProlongation(*pt,r1);
6628 void AliTPCtrackerMI::GetShape(AliTPCseed * seed, Int_t row)
6632 AliTPCClusterParam * clparam = AliTPCcalibDB::Instance()->GetClusterParam();
6633 Float_t zdrift = TMath::Abs((fParam->GetZLength(0)-TMath::Abs(seed->GetZ())));
6634 Int_t type = (seed->GetSector() < fParam->GetNSector()/2) ? 0: (row>126) ? 1:2;
6635 Double_t angulary = seed->GetSnp();
6636 angulary = angulary*angulary/(1.-angulary*angulary);
6637 Double_t angularz = seed->GetTgl()*seed->GetTgl()*(1.+angulary);
6639 Double_t sigmay = clparam->GetRMS0(0,type,zdrift,TMath::Sqrt(TMath::Abs(angulary)));
6640 Double_t sigmaz = clparam->GetRMS0(1,type,zdrift,TMath::Sqrt(TMath::Abs(angularz)));
6641 seed->SetCurrentSigmaY2(sigmay*sigmay);
6642 seed->SetCurrentSigmaZ2(sigmaz*sigmaz);
6643 // Float_t sd2 = TMath::Abs((fParam->GetZLength(0)-TMath::Abs(seed->GetZ())))*fParam->GetDiffL()*fParam->GetDiffL();
6644 // // Float_t padlength = fParam->GetPadPitchLength(seed->fSector);
6645 // Float_t padlength = GetPadPitchLength(row);
6647 // Float_t sresy = (seed->GetSector() < fParam->GetNSector()/2) ? 0.2 :0.3;
6648 // seed->SetCurrentSigmaY2(sd2+padlength*padlength*angulary/12.+sresy*sresy);
6650 // Float_t sresz = fParam->GetZSigma();
6651 // seed->SetCurrentSigmaZ2(sd2+padlength*padlength*angularz*angularz*(1+angulary)/12.+sresz*sresz);
6653 Float_t wy = GetSigmaY(seed);
6654 Float_t wz = GetSigmaZ(seed);
6657 if (TMath::Abs(wy/seed->fCurrentSigmaY2-1)>0.0001 || TMath::Abs(wz/seed->fCurrentSigmaZ2-1)>0.0001 ){
6658 printf("problem\n");
6665 //__________________________________________________________________________
6666 void AliTPCtrackerMI::CookLabel(AliKalmanTrack *tk, Float_t wrong) const {
6667 //--------------------------------------------------------------------
6668 //This function "cooks" a track label. If label<0, this track is fake.
6669 //--------------------------------------------------------------------
6670 AliTPCseed * t = dynamic_cast<AliTPCseed*>(tk);
6672 printf("%s:%d wrong type \n",(char*)__FILE__,__LINE__);
6676 Int_t noc=t->GetNumberOfClusters();
6678 //printf("\nnot founded prolongation\n\n\n");
6684 AliTPCclusterMI *clusters[160];
6686 for (Int_t i=0;i<160;i++) {
6693 for (i=0; i<160 && current<noc; i++) {
6695 Int_t index=t->GetClusterIndex2(i);
6696 if (index<=0) continue;
6697 if (index&0x8000) continue;
6699 //clusters[current]=GetClusterMI(index);
6700 if (t->GetClusterPointer(i)){
6701 clusters[current]=t->GetClusterPointer(i);
6707 Int_t lab=123456789;
6708 for (i=0; i<noc; i++) {
6709 AliTPCclusterMI *c=clusters[i];
6711 lab=TMath::Abs(c->GetLabel(0));
6713 for (j=0; j<noc; j++) if (lb[j]==lab || mx[j]==0) break;
6719 for (i=0; i<noc; i++) if (mx[i]>max) {max=mx[i]; lab=lb[i];}
6721 for (i=0; i<noc; i++) {
6722 AliTPCclusterMI *c=clusters[i];
6724 if (TMath::Abs(c->GetLabel(1)) == lab ||
6725 TMath::Abs(c->GetLabel(2)) == lab ) max++;
6728 if ((1.- Float_t(max)/noc) > wrong) lab=-lab;
6731 Int_t tail=Int_t(0.10*noc);
6734 for (i=1; i<=160&&ind<tail; i++) {
6735 // AliTPCclusterMI *c=clusters[noc-i];
6736 AliTPCclusterMI *c=clusters[i];
6738 if (lab == TMath::Abs(c->GetLabel(0)) ||
6739 lab == TMath::Abs(c->GetLabel(1)) ||
6740 lab == TMath::Abs(c->GetLabel(2))) max++;
6743 if (max < Int_t(0.5*tail)) lab=-lab;
6750 //delete[] clusters;
6754 //__________________________________________________________________________
6755 Int_t AliTPCtrackerMI::CookLabel(AliTPCseed *t, Float_t wrong,Int_t first, Int_t last) const {
6756 //--------------------------------------------------------------------
6757 //This function "cooks" a track label. If label<0, this track is fake.
6758 //--------------------------------------------------------------------
6759 Int_t noc=t->GetNumberOfClusters();
6761 //printf("\nnot founded prolongation\n\n\n");
6767 AliTPCclusterMI *clusters[160];
6769 for (Int_t i=0;i<160;i++) {
6776 for (i=0; i<160 && current<noc; i++) {
6777 if (i<first) continue;
6778 if (i>last) continue;
6779 Int_t index=t->GetClusterIndex2(i);
6780 if (index<=0) continue;
6781 if (index&0x8000) continue;
6783 //clusters[current]=GetClusterMI(index);
6784 if (t->GetClusterPointer(i)){
6785 clusters[current]=t->GetClusterPointer(i);
6790 if (noc<5) return -1;
6791 Int_t lab=123456789;
6792 for (i=0; i<noc; i++) {
6793 AliTPCclusterMI *c=clusters[i];
6795 lab=TMath::Abs(c->GetLabel(0));
6797 for (j=0; j<noc; j++) if (lb[j]==lab || mx[j]==0) break;
6803 for (i=0; i<noc; i++) if (mx[i]>max) {max=mx[i]; lab=lb[i];}
6805 for (i=0; i<noc; i++) {
6806 AliTPCclusterMI *c=clusters[i];
6808 if (TMath::Abs(c->GetLabel(1)) == lab ||
6809 TMath::Abs(c->GetLabel(2)) == lab ) max++;
6812 if ((1.- Float_t(max)/noc) > wrong) lab=-lab;
6815 Int_t tail=Int_t(0.10*noc);
6818 for (i=1; i<=160&&ind<tail; i++) {
6819 // AliTPCclusterMI *c=clusters[noc-i];
6820 AliTPCclusterMI *c=clusters[i];
6822 if (lab == TMath::Abs(c->GetLabel(0)) ||
6823 lab == TMath::Abs(c->GetLabel(1)) ||
6824 lab == TMath::Abs(c->GetLabel(2))) max++;
6827 if (max < Int_t(0.5*tail)) lab=-lab;
6830 // t->SetLabel(lab);
6834 //delete[] clusters;
6838 Int_t AliTPCtrackerMI::AliTPCSector::GetRowNumber(Double_t x) const
6840 //return pad row number for this x
6843 r=fRow[fN-1].GetX();
6844 if (x > r) return fN;
6846 if (x < r) return -1;
6847 return Int_t((x-r)/fPadPitchLength + 0.5);}
6849 r=fRow[fN-1].GetX();
6850 if (x > r) return fN;
6852 if (x < r) return -1;
6853 Double_t r1=fRow[64].GetX();
6855 return Int_t((x-r)/f1PadPitchLength + 0.5);}
6857 return (Int_t((x-r1)/f2PadPitchLength + 0.5)+64);}
6861 Int_t AliTPCtrackerMI::GetRowNumber(Double_t x[3]) const
6863 //return pad row number for given x vector
6864 Float_t phi = TMath::ATan2(x[1],x[0]);
6865 if(phi<0) phi=2.*TMath::Pi()+phi;
6866 // Get the local angle in the sector philoc
6867 const Float_t kRaddeg = 180/3.14159265358979312;
6868 Float_t phiangle = (Int_t (phi*kRaddeg/20.) + 0.5)*20./kRaddeg;
6869 Double_t localx = x[0]*TMath::Cos(phiangle)-x[1]*TMath::Sin(phiangle);
6870 return GetRowNumber(localx);
6873 //_________________________________________________________________________
6874 void AliTPCtrackerMI::AliTPCSector::Setup(const AliTPCParam *par, Int_t f) {
6875 //-----------------------------------------------------------------------
6876 // Setup inner sector
6877 //-----------------------------------------------------------------------
6879 fAlpha=par->GetInnerAngle();
6880 fAlphaShift=par->GetInnerAngleShift();
6881 fPadPitchWidth=par->GetInnerPadPitchWidth();
6882 fPadPitchLength=par->GetInnerPadPitchLength();
6883 fN=par->GetNRowLow();
6884 if(fRow)delete [] fRow;fRow = 0;
6885 fRow=new AliTPCtrackerRow[fN];
6886 for (Int_t i=0; i<fN; i++) {
6887 fRow[i].SetX(par->GetPadRowRadiiLow(i));
6888 fRow[i].SetDeadZone(1.5); //1.5 cm of dead zone
6891 fAlpha=par->GetOuterAngle();
6892 fAlphaShift=par->GetOuterAngleShift();
6893 fPadPitchWidth = par->GetOuterPadPitchWidth();
6894 fPadPitchLength = par->GetOuter1PadPitchLength();
6895 f1PadPitchLength = par->GetOuter1PadPitchLength();
6896 f2PadPitchLength = par->GetOuter2PadPitchLength();
6897 fN=par->GetNRowUp();
6898 if(fRow)delete [] fRow;fRow = 0;
6899 fRow=new AliTPCtrackerRow[fN];
6900 for (Int_t i=0; i<fN; i++) {
6901 fRow[i].SetX(par->GetPadRowRadiiUp(i));
6902 fRow[i].SetDeadZone(1.5); // 1.5 cm of dead zone
6907 AliTPCtrackerRow::AliTPCtrackerRow():
6917 // default constructor
6921 AliTPCtrackerRow::~AliTPCtrackerRow(){
6923 for (Int_t i = 0; i < fN1; i++)
6924 fClusters1[i].~AliTPCclusterMI();
6925 delete [] fClusters1;
6926 for (Int_t i = 0; i < fN2; i++)
6927 fClusters2[i].~AliTPCclusterMI();
6928 delete [] fClusters2;
6933 //_________________________________________________________________________
6935 AliTPCtrackerRow::InsertCluster(const AliTPCclusterMI* c, UInt_t index) {
6936 //-----------------------------------------------------------------------
6937 // Insert a cluster into this pad row in accordence with its y-coordinate
6938 //-----------------------------------------------------------------------
6939 if (fN==kMaxClusterPerRow) {
6940 //AliInfo("AliTPCtrackerRow::InsertCluster(): Too many clusters");
6944 //AliInfo("AliTPCtrackerRow::InsertCluster(): Too many clusters !");
6947 if (fN==0) {fIndex[0]=index; fClusters[fN++]=c; return;}
6948 Int_t i=Find(c->GetZ());
6949 memmove(fClusters+i+1 ,fClusters+i,(fN-i)*sizeof(AliTPCclusterMI*));
6950 memmove(fIndex +i+1 ,fIndex +i,(fN-i)*sizeof(UInt_t));
6951 fIndex[i]=index; fClusters[i]=c; fN++;
6954 void AliTPCtrackerRow::ResetClusters() {
6957 // MvL: Need to call destructors for AliTPCclusterMI, to delete fInfo
6958 for (Int_t i = 0; i < fN1; i++)
6959 fClusters1[i].~AliTPCclusterMI();
6960 delete [] fClusters1; fClusters1=0;
6961 for (Int_t i = 0; i < fN2; i++)
6962 fClusters2[i].~AliTPCclusterMI();
6963 delete [] fClusters2; fClusters2=0;
6968 //delete[] fClusterArray;
6974 //___________________________________________________________________
6975 Int_t AliTPCtrackerRow::Find(Double_t z) const {
6976 //-----------------------------------------------------------------------
6977 // Return the index of the nearest cluster
6978 //-----------------------------------------------------------------------
6979 if (fN==0) return 0;
6980 if (z <= fClusters[0]->GetZ()) return 0;
6981 if (z > fClusters[fN-1]->GetZ()) return fN;
6982 Int_t b=0, e=fN-1, m=(b+e)/2;
6983 for (; b<e; m=(b+e)/2) {
6984 if (z > fClusters[m]->GetZ()) b=m+1;
6992 //___________________________________________________________________
6993 AliTPCclusterMI * AliTPCtrackerRow::FindNearest(Double_t y, Double_t z, Double_t roady, Double_t roadz) const {
6994 //-----------------------------------------------------------------------
6995 // Return the index of the nearest cluster in z y
6996 //-----------------------------------------------------------------------
6997 Float_t maxdistance = roady*roady + roadz*roadz;
6999 AliTPCclusterMI *cl =0;
7000 for (Int_t i=Find(z-roadz); i<fN; i++) {
7001 AliTPCclusterMI *c=(AliTPCclusterMI*)(fClusters[i]);
7002 if (c->GetZ() > z+roadz) break;
7003 if ( (c->GetY()-y) > roady ) continue;
7004 Float_t distance = (c->GetZ()-z)*(c->GetZ()-z)+(c->GetY()-y)*(c->GetY()-y);
7005 if (maxdistance>distance) {
7006 maxdistance = distance;
7013 AliTPCclusterMI * AliTPCtrackerRow::FindNearest2(Double_t y, Double_t z, Double_t roady, Double_t roadz,UInt_t & index) const
7015 //-----------------------------------------------------------------------
7016 // Return the index of the nearest cluster in z y
7017 //-----------------------------------------------------------------------
7018 Float_t maxdistance = roady*roady + roadz*roadz;
7019 AliTPCclusterMI *cl =0;
7021 //PH Check boundaries. 510 is the size of fFastCluster
7022 Int_t iz1 = Int_t(z-roadz+254.5);
7023 if (iz1<0 || iz1>=510) return cl;
7024 iz1 = TMath::Max(GetFastCluster(iz1)-1,0);
7025 Int_t iz2 = Int_t(z+roadz+255.5);
7026 if (iz2<0 || iz2>=510) return cl;
7027 iz2 = TMath::Min(GetFastCluster(iz2)+1,fN);
7028 Bool_t skipUsed = !(AliTPCReconstructor::GetRecoParam()->GetClusterSharing());
7029 //FindNearest3(y,z,roady,roadz,index);
7030 // for (Int_t i=Find(z-roadz); i<fN; i++) {
7031 for (Int_t i=iz1; i<iz2; i++) {
7032 AliTPCclusterMI *c=(AliTPCclusterMI*)(fClusters[i]);
7033 if (c->GetZ() > z+roadz) break;
7034 if ( c->GetY()-y > roady ) continue;
7035 if ( y-c->GetY() > roady ) continue;
7036 if (skipUsed && c->IsUsed(11)) continue;
7037 Float_t distance = (c->GetZ()-z)*(c->GetZ()-z)+(c->GetY()-y)*(c->GetY()-y);
7038 if (maxdistance>distance) {
7039 maxdistance = distance;
7042 //roady = TMath::Sqrt(maxdistance);
7049 void AliTPCtrackerRow::SetFastCluster(Int_t i, Short_t cl){
7051 // Set cluster info for fast navigation
7062 void AliTPCtrackerMI::MakeBitmaps(AliTPCseed *t)
7064 //-----------------------------------------------------------------------
7065 // Fill the cluster and sharing bitmaps of the track
7066 //-----------------------------------------------------------------------
7068 Int_t firstpoint = 0;
7069 Int_t lastpoint = 159;
7070 AliTPCTrackerPoint *point;
7072 for (int iter=firstpoint; iter<lastpoint; iter++) {
7073 point = t->GetTrackPoint(iter);
7075 t->SetClusterMapBit(iter, kTRUE);
7076 if (point->IsShared())
7077 t->SetSharedMapBit(iter,kTRUE);
7079 t->SetSharedMapBit(iter, kFALSE);
7082 t->SetClusterMapBit(iter, kFALSE);
7083 t->SetSharedMapBit(iter, kFALSE);