1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
17 //-------------------------------------------------------
18 // Implementation of the TPC tracker
20 // Origin: Marian Ivanov Marian.Ivanov@cern.ch
22 // AliTPC parallel tracker
24 // The track fitting is based on Kalaman filtering approach
26 // The track finding steps:
27 // 1. Seeding - with and without vertex constraint
28 // - seeding with vertex constain done at first n^2 proble
29 // - seeding without vertex constraint n^3 problem
30 // 2. Tracking - follow prolongation road - find cluster - update kalman track
32 // The seeding and tracking is repeated several times, in different seeding region.
33 // This approach enables to find the track which cannot be seeded in some region of TPC
34 // This can happen because of low momenta (track do not reach outer radius), or track is currently in the ded region between sectors, or the track is for the moment overlapped with other track (seed quality is poor) ...
36 // With this approach we reach almost 100 % efficiency also for high occupancy events.
37 // (If the seeding efficiency in a region is about 90 % than with logical or of several
38 // regions we will reach 100% (in theory - supposing independence)
40 // Repeating several seeding - tracking procedures some of the tracks can be find
43 // The procedures to remove multi find tacks are impremented:
44 // RemoveUsed2 - fast procedure n problem -
45 // Algorithm - Sorting tracks according quality
46 // remove tracks with some shared fraction
47 // Sharing in respect to all tacks
48 // Signing clusters in gold region
49 // FindSplitted - slower algorithm n^2
50 // Sort the tracks according quality
51 // Loop over pair of tracks
52 // If overlap with other track bigger than threshold - remove track
54 // FindCurling - Finds the pair of tracks which are curling
55 // - About 10% of tracks can be find with this procedure
56 // The combinatorial background is too big to be used in High
57 // multiplicity environment
58 // - n^2 problem - Slow procedure - currently it is disabled because of
61 // The number of splitted tracks can be reduced disabling the sharing of the cluster.
62 // tpcRecoParam-> SetClusterSharing(kFALSE);
63 // IT IS HIGHLY non recomended to use it in high flux enviroonment
64 // Even using this switch some tracks can be found more than once
65 // (because of multiple seeding and low quality tracks which will not cross full chamber)
68 // The tracker itself can be debugged - the information about tracks can be stored in several // phases of the reconstruction
69 // To enable storage of the TPC tracks in the ESD friend track
70 // use AliTPCReconstructor::SetStreamLevel(n); where nis bigger 0
72 // The debug level - different procedure produce tree for numerical debugging
73 // To enable them set AliTPCReconstructor::SetStreamLevel(n); where nis bigger 1
75 // There are several places in the code which can be numerically debuged
76 // This code is keeped in order to enable code development and to check the calibration implementtion
78 // 1. ErrParam stream (Log level 9) - dump information about
80 // 2.a) cluster error estimate
81 // 3.a) cluster shape estimate
84 //-------------------------------------------------------
89 #include "Riostream.h"
90 #include <TClonesArray.h>
92 #include <TObjArray.h>
95 #include "AliComplexCluster.h"
96 #include "AliESDEvent.h"
97 #include "AliESDtrack.h"
98 #include "AliESDVertex.h"
101 #include "AliHelix.h"
102 #include "AliRunLoader.h"
103 #include "AliTPCClustersRow.h"
104 #include "AliTPCParam.h"
105 #include "AliTPCReconstructor.h"
106 #include "AliTPCpolyTrack.h"
107 #include "AliTPCreco.h"
108 #include "AliTPCseed.h"
110 #include "AliTPCtrackerSector.h"
111 #include "AliTPCtrackerMI.h"
112 #include "TStopwatch.h"
113 #include "AliTPCReconstructor.h"
115 #include "TTreeStream.h"
116 #include "AliAlignObj.h"
117 #include "AliTrackPointArray.h"
119 #include "AliTPCcalibDB.h"
120 #include "AliTPCTransform.h"
121 #include "AliTPCClusterParam.h"
125 ClassImp(AliTPCtrackerMI)
129 class AliTPCFastMath {
132 static Double_t FastAsin(Double_t x);
134 static Double_t fgFastAsin[20000]; //lookup table for fast asin computation
137 Double_t AliTPCFastMath::fgFastAsin[20000];
138 AliTPCFastMath gAliTPCFastMath; // needed to fill the LUT
140 AliTPCFastMath::AliTPCFastMath(){
142 // initialized lookup table;
143 for (Int_t i=0;i<10000;i++){
144 fgFastAsin[2*i] = TMath::ASin(i/10000.);
145 fgFastAsin[2*i+1] = (TMath::ASin((i+1)/10000.)-fgFastAsin[2*i]);
149 Double_t AliTPCFastMath::FastAsin(Double_t x){
151 // return asin using lookup table
153 Int_t index = int(x*10000);
154 return fgFastAsin[2*index]+(x*10000.-index)*fgFastAsin[2*index+1];
157 Int_t index = int(x*10000);
158 return -(fgFastAsin[2*index]+(x*10000.-index)*fgFastAsin[2*index+1]);
160 //__________________________________________________________________
161 AliTPCtrackerMI::AliTPCtrackerMI()
183 // default constructor
186 //_____________________________________________________________________
190 Int_t AliTPCtrackerMI::UpdateTrack(AliTPCseed * track, Int_t accept){
192 //update track information using current cluster - track->fCurrentCluster
195 AliTPCclusterMI* c =track->GetCurrentCluster();
196 if (accept>0) track->SetCurrentClusterIndex1(track->GetCurrentClusterIndex1() | 0x8000); //sign not accepted clusters
198 UInt_t i = track->GetCurrentClusterIndex1();
200 Int_t sec=(i&0xff000000)>>24;
201 //Int_t row = (i&0x00ff0000)>>16;
202 track->SetRow((i&0x00ff0000)>>16);
203 track->SetSector(sec);
204 // Int_t index = i&0xFFFF;
205 if (sec>=fParam->GetNInnerSector()) track->SetRow(track->GetRow()+fParam->GetNRowLow());
206 track->SetClusterIndex2(track->GetRow(), i);
207 //track->fFirstPoint = row;
208 //if ( track->fLastPoint<row) track->fLastPoint =row;
209 // if (track->fRow<0 || track->fRow>160) {
210 // printf("problem\n");
212 if (track->GetFirstPoint()>track->GetRow())
213 track->SetFirstPoint(track->GetRow());
214 if (track->GetLastPoint()<track->GetRow())
215 track->SetLastPoint(track->GetRow());
218 track->SetClusterPointer(track->GetRow(),c);
221 Double_t angle2 = track->GetSnp()*track->GetSnp();
223 //SET NEW Track Point
225 if (angle2<1) //PH sometimes angle2 is very big. To be investigated...
227 angle2 = TMath::Sqrt(angle2/(1-angle2));
228 AliTPCTrackerPoint &point =*(track->GetTrackPoint(track->GetRow()));
230 point.SetSigmaY(c->GetSigmaY2()/track->GetCurrentSigmaY2());
231 point.SetSigmaZ(c->GetSigmaZ2()/track->GetCurrentSigmaZ2());
232 point.SetErrY(sqrt(track->GetErrorY2()));
233 point.SetErrZ(sqrt(track->GetErrorZ2()));
235 point.SetX(track->GetX());
236 point.SetY(track->GetY());
237 point.SetZ(track->GetZ());
238 point.SetAngleY(angle2);
239 point.SetAngleZ(track->GetTgl());
240 if (point.IsShared()){
241 track->SetErrorY2(track->GetErrorY2()*4);
242 track->SetErrorZ2(track->GetErrorZ2()*4);
246 Double_t chi2 = track->GetPredictedChi2(track->GetCurrentCluster());
248 // track->SetErrorY2(track->GetErrorY2()*1.3);
249 // track->SetErrorY2(track->GetErrorY2()+0.01);
250 // track->SetErrorZ2(track->GetErrorZ2()*1.3);
251 // track->SetErrorZ2(track->GetErrorZ2()+0.005);
253 if (accept>0) return 0;
254 if (track->GetNumberOfClusters()%20==0){
255 // if (track->fHelixIn){
256 // TClonesArray & larr = *(track->fHelixIn);
257 // Int_t ihelix = larr.GetEntriesFast();
258 // new(larr[ihelix]) AliHelix(*track) ;
261 track->SetNoCluster(0);
262 return track->Update(c,chi2,i);
267 Int_t AliTPCtrackerMI::AcceptCluster(AliTPCseed * seed, AliTPCclusterMI * cluster)
270 // decide according desired precision to accept given
271 // cluster for tracking
272 Double_t sy2=ErrY2(seed,cluster);
273 Double_t sz2=ErrZ2(seed,cluster);
275 Double_t sdistancey2 = sy2+seed->GetSigmaY2();
276 Double_t sdistancez2 = sz2+seed->GetSigmaZ2();
278 Double_t rdistancey2 = (seed->GetCurrentCluster()->GetY()-seed->GetY())*
279 (seed->GetCurrentCluster()->GetY()-seed->GetY())/sdistancey2;
280 Double_t rdistancez2 = (seed->GetCurrentCluster()->GetZ()-seed->GetZ())*
281 (seed->GetCurrentCluster()->GetZ()-seed->GetZ())/sdistancez2;
283 Double_t rdistance2 = rdistancey2+rdistancez2;
286 if (AliTPCReconstructor::StreamLevel()>5) {
287 Float_t rmsy2 = seed->GetCurrentSigmaY2();
288 Float_t rmsz2 = seed->GetCurrentSigmaZ2();
289 Float_t rmsy2p30 = seed->GetCMeanSigmaY2p30();
290 Float_t rmsz2p30 = seed->GetCMeanSigmaZ2p30();
291 Float_t rmsy2p30R = seed->GetCMeanSigmaY2p30R();
292 Float_t rmsz2p30R = seed->GetCMeanSigmaZ2p30R();
294 AliExternalTrackParam param(*seed);
295 static TVectorD gcl(3),gtr(3);
297 param.GetXYZ(gcl.GetMatrixArray());
298 cluster->GetGlobalXYZ(gclf);
299 gcl[0]=gclf[0]; gcl[1]=gclf[1]; gcl[2]=gclf[2];
300 (*fDebugStreamer)<<"ErrParam"<<
309 "rmsy2p30="<<rmsy2p30<<
310 "rmsz2p30="<<rmsz2p30<<
311 "rmsy2p30R="<<rmsy2p30R<<
312 "rmsz2p30R="<<rmsz2p30R<<
316 if (rdistance2>16) return 3;
319 if ((rdistancey2>9. || rdistancez2>9.) && cluster->GetType()==0)
320 return 2; //suspisiouce - will be changed
322 if ((rdistancey2>6.25 || rdistancez2>6.25) && cluster->GetType()>0)
323 // strict cut on overlaped cluster
324 return 2; //suspisiouce - will be changed
326 if ( (rdistancey2>1. || rdistancez2>6.25 )
327 && cluster->GetType()<0){
328 seed->SetNFoundable(seed->GetNFoundable()-1);
337 //_____________________________________________________________________________
338 AliTPCtrackerMI::AliTPCtrackerMI(const AliTPCParam *par):
340 fkNIS(par->GetNInnerSector()/2),
342 fkNOS(par->GetNOuterSector()/2),
359 //---------------------------------------------------------------------
360 // The main TPC tracker constructor
361 //---------------------------------------------------------------------
362 fInnerSec=new AliTPCtrackerSector[fkNIS];
363 fOuterSec=new AliTPCtrackerSector[fkNOS];
366 for (i=0; i<fkNIS; i++) fInnerSec[i].Setup(par,0);
367 for (i=0; i<fkNOS; i++) fOuterSec[i].Setup(par,1);
370 Int_t nrowlow = par->GetNRowLow();
371 Int_t nrowup = par->GetNRowUp();
374 for (Int_t i=0;i<nrowlow;i++){
375 fXRow[i] = par->GetPadRowRadiiLow(i);
376 fPadLength[i]= par->GetPadPitchLength(0,i);
377 fYMax[i] = fXRow[i]*TMath::Tan(0.5*par->GetInnerAngle());
381 for (Int_t i=0;i<nrowup;i++){
382 fXRow[i+nrowlow] = par->GetPadRowRadiiUp(i);
383 fPadLength[i+nrowlow] = par->GetPadPitchLength(60,i);
384 fYMax[i+nrowlow] = fXRow[i+nrowlow]*TMath::Tan(0.5*par->GetOuterAngle());
387 fDebugStreamer = new TTreeSRedirector("TPCdebug.root");
389 //________________________________________________________________________
390 AliTPCtrackerMI::AliTPCtrackerMI(const AliTPCtrackerMI &t):
411 //------------------------------------
412 // dummy copy constructor
413 //------------------------------------------------------------------
416 AliTPCtrackerMI & AliTPCtrackerMI::operator=(const AliTPCtrackerMI& /*r*/){
417 //------------------------------
419 //--------------------------------------------------------------
422 //_____________________________________________________________________________
423 AliTPCtrackerMI::~AliTPCtrackerMI() {
424 //------------------------------------------------------------------
425 // TPC tracker destructor
426 //------------------------------------------------------------------
433 if (fDebugStreamer) delete fDebugStreamer;
437 void AliTPCtrackerMI::FillESD(TObjArray* arr)
441 //fill esds using updated tracks
443 // write tracks to the event
444 // store index of the track
445 Int_t nseed=arr->GetEntriesFast();
446 //FindKinks(arr,fEvent);
447 for (Int_t i=0; i<nseed; i++) {
448 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
451 // pt->PropagateTo(fParam->GetInnerRadiusLow());
452 if (pt->GetKinkIndex(0)<=0){ //don't propagate daughter tracks
453 pt->PropagateTo(fParam->GetInnerRadiusLow());
456 if (( pt->GetPoints()[2]- pt->GetPoints()[0])>5 && pt->GetPoints()[3]>0.8){
458 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
459 iotrack.SetTPCPoints(pt->GetPoints());
460 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
461 iotrack.SetV0Indexes(pt->GetV0Indexes());
462 // iotrack.SetTPCpid(pt->fTPCr);
463 //iotrack.SetTPCindex(i);
464 fEvent->AddTrack(&iotrack);
468 if ( (pt->GetNumberOfClusters()>70)&& (Float_t(pt->GetNumberOfClusters())/Float_t(pt->GetNFoundable()))>0.55) {
470 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
471 iotrack.SetTPCPoints(pt->GetPoints());
472 //iotrack.SetTPCindex(i);
473 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
474 iotrack.SetV0Indexes(pt->GetV0Indexes());
475 // iotrack.SetTPCpid(pt->fTPCr);
476 fEvent->AddTrack(&iotrack);
480 // short tracks - maybe decays
482 if ( (pt->GetNumberOfClusters()>30) && (Float_t(pt->GetNumberOfClusters())/Float_t(pt->GetNFoundable()))>0.70) {
483 Int_t found,foundable,shared;
484 pt->GetClusterStatistic(0,60,found, foundable,shared,kFALSE);
485 if ( (found>20) && (pt->GetNShared()/float(pt->GetNumberOfClusters())<0.2)){
487 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
488 //iotrack.SetTPCindex(i);
489 iotrack.SetTPCPoints(pt->GetPoints());
490 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
491 iotrack.SetV0Indexes(pt->GetV0Indexes());
492 //iotrack.SetTPCpid(pt->fTPCr);
493 fEvent->AddTrack(&iotrack);
498 if ( (pt->GetNumberOfClusters()>20) && (Float_t(pt->GetNumberOfClusters())/Float_t(pt->GetNFoundable()))>0.8) {
499 Int_t found,foundable,shared;
500 pt->GetClusterStatistic(0,60,found, foundable,shared,kFALSE);
501 if (found<20) continue;
502 if (pt->GetNShared()/float(pt->GetNumberOfClusters())>0.2) continue;
505 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
506 iotrack.SetTPCPoints(pt->GetPoints());
507 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
508 iotrack.SetV0Indexes(pt->GetV0Indexes());
509 //iotrack.SetTPCpid(pt->fTPCr);
510 //iotrack.SetTPCindex(i);
511 fEvent->AddTrack(&iotrack);
514 // short tracks - secondaties
516 if ( (pt->GetNumberOfClusters()>30) ) {
517 Int_t found,foundable,shared;
518 pt->GetClusterStatistic(128,158,found, foundable,shared,kFALSE);
519 if ( (found>20) && (pt->GetNShared()/float(pt->GetNumberOfClusters())<0.2) &&float(found)/float(foundable)>0.8){
521 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
522 iotrack.SetTPCPoints(pt->GetPoints());
523 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
524 iotrack.SetV0Indexes(pt->GetV0Indexes());
525 //iotrack.SetTPCpid(pt->fTPCr);
526 //iotrack.SetTPCindex(i);
527 fEvent->AddTrack(&iotrack);
532 if ( (pt->GetNumberOfClusters()>15)) {
533 Int_t found,foundable,shared;
534 pt->GetClusterStatistic(138,158,found, foundable,shared,kFALSE);
535 if (found<15) continue;
536 if (foundable<=0) continue;
537 if (pt->GetNShared()/float(pt->GetNumberOfClusters())>0.2) continue;
538 if (float(found)/float(foundable)<0.8) continue;
541 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
542 iotrack.SetTPCPoints(pt->GetPoints());
543 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
544 iotrack.SetV0Indexes(pt->GetV0Indexes());
545 // iotrack.SetTPCpid(pt->fTPCr);
546 //iotrack.SetTPCindex(i);
547 fEvent->AddTrack(&iotrack);
552 printf("Number of filled ESDs-\t%d\n",fEvent->GetNumberOfTracks());
559 Double_t AliTPCtrackerMI::ErrY2(AliTPCseed* seed, AliTPCclusterMI * cl){
562 // Use calibrated cluster error from OCDB
564 AliTPCClusterParam * clparam = AliTPCcalibDB::Instance()->GetClusterParam();
566 Float_t z = TMath::Abs(fParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
567 Int_t ctype = cl->GetType();
568 Int_t type = (cl->GetRow()<63) ? 0: (cl->GetRow()>126) ? 1:2;
569 Double_t angle = seed->GetSnp()*seed->GetSnp();
570 angle = TMath::Sqrt(TMath::Abs(angle/(1.-angle)));
571 Double_t erry2 = clparam->GetError0Par(0,type, z,angle);
573 erry2+=0.5; // edge cluster
576 seed->SetErrorY2(erry2);
580 //calculate look-up table at the beginning
581 // static Bool_t ginit = kFALSE;
582 // static Float_t gnoise1,gnoise2,gnoise3;
583 // static Float_t ggg1[10000];
584 // static Float_t ggg2[10000];
585 // static Float_t ggg3[10000];
586 // static Float_t glandau1[10000];
587 // static Float_t glandau2[10000];
588 // static Float_t glandau3[10000];
590 // static Float_t gcor01[500];
591 // static Float_t gcor02[500];
592 // static Float_t gcorp[500];
596 // if (ginit==kFALSE){
597 // for (Int_t i=1;i<500;i++){
598 // Float_t rsigma = float(i)/100.;
599 // gcor02[i] = TMath::Max(0.78 +TMath::Exp(7.4*(rsigma-1.2)),0.6);
600 // gcor01[i] = TMath::Max(0.72 +TMath::Exp(3.36*(rsigma-1.2)),0.6);
601 // gcorp[i] = TMath::Max(TMath::Power((rsigma+0.5),1.5),1.2);
605 // for (Int_t i=3;i<10000;i++){
609 // Float_t amp = float(i);
610 // Float_t padlength =0.75;
611 // gnoise1 = 0.0004/padlength;
612 // Float_t nel = 0.268*amp;
613 // Float_t nprim = 0.155*amp;
614 // ggg1[i] = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.001*nel/(padlength*padlength))/nel;
615 // glandau1[i] = (2.+0.12*nprim)*0.5* (2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
616 // if (glandau1[i]>1) glandau1[i]=1;
617 // glandau1[i]*=padlength*padlength/12.;
621 // gnoise2 = 0.0004/padlength;
623 // nprim = 0.133*amp;
624 // ggg2[i] = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
625 // glandau2[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
626 // if (glandau2[i]>1) glandau2[i]=1;
627 // glandau2[i]*=padlength*padlength/12.;
632 // gnoise3 = 0.0004/padlength;
634 // nprim = 0.133*amp;
635 // ggg3[i] = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
636 // glandau3[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
637 // if (glandau3[i]>1) glandau3[i]=1;
638 // glandau3[i]*=padlength*padlength/12.;
646 // Int_t amp = int(TMath::Abs(cl->GetQ()));
648 // seed->SetErrorY2(1.);
652 // Float_t z = TMath::Abs(fParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
653 // Int_t ctype = cl->GetType();
654 // Float_t padlength= GetPadPitchLength(seed->GetRow());
655 // Double_t angle2 = seed->GetSnp()*seed->GetSnp();
656 // angle2 = angle2/(1-angle2);
658 // //cluster "quality"
659 // Int_t rsigmay = int(100.*cl->GetSigmaY2()/(seed->GetCurrentSigmaY2()));
662 // if (fSectors==fInnerSec){
663 // snoise2 = gnoise1;
664 // res = ggg1[amp]*z+glandau1[amp]*angle2;
665 // if (ctype==0) res *= gcor01[rsigmay];
668 // res*= gcorp[rsigmay];
672 // if (padlength<1.1){
673 // snoise2 = gnoise2;
674 // res = ggg2[amp]*z+glandau2[amp]*angle2;
675 // if (ctype==0) res *= gcor02[rsigmay];
678 // res*= gcorp[rsigmay];
682 // snoise2 = gnoise3;
683 // res = ggg3[amp]*z+glandau3[amp]*angle2;
684 // if (ctype==0) res *= gcor02[rsigmay];
687 // res*= gcorp[rsigmay];
694 // res*=2.4; // overestimate error 2 times
698 // if (res<2*snoise2)
701 // seed->SetErrorY2(res);
709 Double_t AliTPCtrackerMI::ErrZ2(AliTPCseed* seed, AliTPCclusterMI * cl){
712 // Use calibrated cluster error from OCDB
714 AliTPCClusterParam * clparam = AliTPCcalibDB::Instance()->GetClusterParam();
716 Float_t z = TMath::Abs(fParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
717 Int_t ctype = cl->GetType();
718 Int_t type = (cl->GetRow()<63) ? 0: (cl->GetRow()>126) ? 1:2;
720 Double_t angle2 = seed->GetSnp()*seed->GetSnp();
721 angle2 = seed->GetTgl()*seed->GetTgl()*(1+angle2/(1-angle2));
722 Double_t angle = TMath::Sqrt(TMath::Abs(angle2));
723 Double_t errz2 = clparam->GetError0Par(1,type, z,angle);
725 errz2+=0.5; // edge cluster
728 seed->SetErrorZ2(errz2);
734 // //seed->SetErrorY2(0.1);
736 // //calculate look-up table at the beginning
737 // static Bool_t ginit = kFALSE;
738 // static Float_t gnoise1,gnoise2,gnoise3;
739 // static Float_t ggg1[10000];
740 // static Float_t ggg2[10000];
741 // static Float_t ggg3[10000];
742 // static Float_t glandau1[10000];
743 // static Float_t glandau2[10000];
744 // static Float_t glandau3[10000];
746 // static Float_t gcor01[1000];
747 // static Float_t gcor02[1000];
748 // static Float_t gcorp[1000];
752 // if (ginit==kFALSE){
753 // for (Int_t i=1;i<1000;i++){
754 // Float_t rsigma = float(i)/100.;
755 // gcor02[i] = TMath::Max(0.81 +TMath::Exp(6.8*(rsigma-1.2)),0.6);
756 // gcor01[i] = TMath::Max(0.72 +TMath::Exp(2.04*(rsigma-1.2)),0.6);
757 // gcorp[i] = TMath::Max(TMath::Power((rsigma+0.5),1.5),1.2);
761 // for (Int_t i=3;i<10000;i++){
765 // Float_t amp = float(i);
766 // Float_t padlength =0.75;
767 // gnoise1 = 0.0004/padlength;
768 // Float_t nel = 0.268*amp;
769 // Float_t nprim = 0.155*amp;
770 // ggg1[i] = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.001*nel/(padlength*padlength))/nel;
771 // glandau1[i] = (2.+0.12*nprim)*0.5* (2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
772 // if (glandau1[i]>1) glandau1[i]=1;
773 // glandau1[i]*=padlength*padlength/12.;
777 // gnoise2 = 0.0004/padlength;
779 // nprim = 0.133*amp;
780 // ggg2[i] = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
781 // glandau2[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
782 // if (glandau2[i]>1) glandau2[i]=1;
783 // glandau2[i]*=padlength*padlength/12.;
788 // gnoise3 = 0.0004/padlength;
790 // nprim = 0.133*amp;
791 // ggg3[i] = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
792 // glandau3[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
793 // if (glandau3[i]>1) glandau3[i]=1;
794 // glandau3[i]*=padlength*padlength/12.;
802 // Int_t amp = int(TMath::Abs(cl->GetQ()));
804 // seed->SetErrorY2(1.);
808 // Float_t z = TMath::Abs(fParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
809 // Int_t ctype = cl->GetType();
810 // Float_t padlength= GetPadPitchLength(seed->GetRow());
812 // Double_t angle2 = seed->GetSnp()*seed->GetSnp();
813 // // if (angle2<0.6) angle2 = 0.6;
814 // angle2 = seed->GetTgl()*seed->GetTgl()*(1+angle2/(1-angle2));
816 // //cluster "quality"
817 // Int_t rsigmaz = int(100.*cl->GetSigmaZ2()/(seed->GetCurrentSigmaZ2()));
820 // if (fSectors==fInnerSec){
821 // snoise2 = gnoise1;
822 // res = ggg1[amp]*z+glandau1[amp]*angle2;
823 // if (ctype==0) res *= gcor01[rsigmaz];
826 // res*= gcorp[rsigmaz];
830 // if (padlength<1.1){
831 // snoise2 = gnoise2;
832 // res = ggg2[amp]*z+glandau2[amp]*angle2;
833 // if (ctype==0) res *= gcor02[rsigmaz];
836 // res*= gcorp[rsigmaz];
840 // snoise2 = gnoise3;
841 // res = ggg3[amp]*z+glandau3[amp]*angle2;
842 // if (ctype==0) res *= gcor02[rsigmaz];
845 // res*= gcorp[rsigmaz];
854 // if ((ctype<0) &&<70){
859 // if (res<2*snoise2)
861 // if (res>3) res =3;
862 // seed->SetErrorZ2(res);
870 void AliTPCtrackerMI::RotateToLocal(AliTPCseed *seed)
872 //rotate to track "local coordinata
873 Float_t x = seed->GetX();
874 Float_t y = seed->GetY();
875 Float_t ymax = x*TMath::Tan(0.5*fSectors->GetAlpha());
878 seed->SetRelativeSector((seed->GetRelativeSector()+1) % fN);
879 if (!seed->Rotate(fSectors->GetAlpha()))
881 } else if (y <-ymax) {
882 seed->SetRelativeSector((seed->GetRelativeSector()-1+fN) % fN);
883 if (!seed->Rotate(-fSectors->GetAlpha()))
891 //_____________________________________________________________________________
892 Double_t AliTPCtrackerMI::F1old(Double_t x1,Double_t y1,
893 Double_t x2,Double_t y2,
894 Double_t x3,Double_t y3)
896 //-----------------------------------------------------------------
897 // Initial approximation of the track curvature
898 //-----------------------------------------------------------------
899 Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
900 Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
901 (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
902 Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
903 (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
905 Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
906 if ( xr*xr+yr*yr<=0.00000000000001) return 100;
907 return -xr*yr/sqrt(xr*xr+yr*yr);
912 //_____________________________________________________________________________
913 Double_t AliTPCtrackerMI::F1(Double_t x1,Double_t y1,
914 Double_t x2,Double_t y2,
915 Double_t x3,Double_t y3)
917 //-----------------------------------------------------------------
918 // Initial approximation of the track curvature
919 //-----------------------------------------------------------------
925 Double_t det = x3*y2-x2*y3;
930 Double_t u = 0.5* (x2*(x2-x3)+y2*(y2-y3))/det;
931 Double_t x0 = x3*0.5-y3*u;
932 Double_t y0 = y3*0.5+x3*u;
933 Double_t c2 = 1/TMath::Sqrt(x0*x0+y0*y0);
939 Double_t AliTPCtrackerMI::F2(Double_t x1,Double_t y1,
940 Double_t x2,Double_t y2,
941 Double_t x3,Double_t y3)
943 //-----------------------------------------------------------------
944 // Initial approximation of the track curvature
945 //-----------------------------------------------------------------
951 Double_t det = x3*y2-x2*y3;
956 Double_t u = 0.5* (x2*(x2-x3)+y2*(y2-y3))/det;
957 Double_t x0 = x3*0.5-y3*u;
958 Double_t y0 = y3*0.5+x3*u;
959 Double_t c2 = 1/TMath::Sqrt(x0*x0+y0*y0);
968 //_____________________________________________________________________________
969 Double_t AliTPCtrackerMI::F2old(Double_t x1,Double_t y1,
970 Double_t x2,Double_t y2,
971 Double_t x3,Double_t y3)
973 //-----------------------------------------------------------------
974 // Initial approximation of the track curvature times center of curvature
975 //-----------------------------------------------------------------
976 Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
977 Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
978 (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
979 Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
980 (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
982 Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
984 return -a/(d*y1-b)*xr/sqrt(xr*xr+yr*yr);
987 //_____________________________________________________________________________
988 Double_t AliTPCtrackerMI::F3(Double_t x1,Double_t y1,
989 Double_t x2,Double_t y2,
990 Double_t z1,Double_t z2)
992 //-----------------------------------------------------------------
993 // Initial approximation of the tangent of the track dip angle
994 //-----------------------------------------------------------------
995 return (z1 - z2)/sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
999 Double_t AliTPCtrackerMI::F3n(Double_t x1,Double_t y1,
1000 Double_t x2,Double_t y2,
1001 Double_t z1,Double_t z2, Double_t c)
1003 //-----------------------------------------------------------------
1004 // Initial approximation of the tangent of the track dip angle
1005 //-----------------------------------------------------------------
1009 //angle1 = (z1-z2)*c/(TMath::ASin(c*x1-ni)-TMath::ASin(c*x2-ni));
1011 Double_t d = TMath::Sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
1012 if (TMath::Abs(d*c*0.5)>1) return 0;
1013 // Double_t angle2 = TMath::ASin(d*c*0.5);
1014 // Double_t angle2 = AliTPCFastMath::FastAsin(d*c*0.5);
1015 Double_t angle2 = (d*c*0.5>0.1)? TMath::ASin(d*c*0.5): AliTPCFastMath::FastAsin(d*c*0.5);
1017 angle2 = (z1-z2)*c/(angle2*2.);
1021 Bool_t AliTPCtrackerMI::GetProlongation(Double_t x1, Double_t x2, Double_t x[5], Double_t &y, Double_t &z)
1022 {//-----------------------------------------------------------------
1023 // This function find proloncation of a track to a reference plane x=x2.
1024 //-----------------------------------------------------------------
1028 if (TMath::Abs(x[4]*x1 - x[2]) >= 0.999) {
1032 Double_t c1=x[4]*x1 - x[2], r1=sqrt(1.- c1*c1);
1033 Double_t c2=x[4]*x2 - x[2], r2=sqrt(1.- c2*c2);
1037 Double_t dy = dx*(c1+c2)/(r1+r2);
1040 Double_t delta = x[4]*dx*(c1+c2)/(c1*r2 + c2*r1);
1042 if (TMath::Abs(delta)>0.01){
1043 dz = x[3]*TMath::ASin(delta)/x[4];
1045 dz = x[3]*AliTPCFastMath::FastAsin(delta)/x[4];
1048 //dz = x[3]*AliTPCFastMath::FastAsin(delta)/x[4];
1056 Int_t AliTPCtrackerMI::LoadClusters (TTree *tree)
1061 return LoadClusters();
1065 Int_t AliTPCtrackerMI::LoadClusters(TObjArray *arr)
1068 // load clusters to the memory
1069 AliTPCClustersRow *clrow = 0x0;
1070 Int_t lower = arr->LowerBound();
1071 Int_t entries = arr->GetEntriesFast();
1073 for (Int_t i=lower; i<entries; i++) {
1074 clrow = (AliTPCClustersRow*) arr->At(i);
1075 if(!clrow) continue;
1076 if(!clrow->GetArray()) continue;
1080 fParam->AdjustSectorRow(clrow->GetID(),sec,row);
1082 for (Int_t icl=0; icl<clrow->GetArray()->GetEntriesFast(); icl++){
1083 Transform((AliTPCclusterMI*)(clrow->GetArray()->At(icl)));
1086 if (clrow->GetArray()->GetEntriesFast()<=0) continue;
1087 AliTPCtrackerRow * tpcrow=0;
1090 tpcrow = &(fInnerSec[sec%fkNIS][row]);
1094 tpcrow = &(fOuterSec[(sec-fkNIS*2)%fkNOS][row]);
1095 left = (sec-fkNIS*2)/fkNOS;
1098 tpcrow->SetN1(clrow->GetArray()->GetEntriesFast());
1099 tpcrow->SetClusters1(new AliTPCclusterMI[tpcrow->GetN1()]);
1100 for (Int_t i=0;i<tpcrow->GetN1();i++)
1101 tpcrow->SetCluster1(i, *(AliTPCclusterMI*)(clrow->GetArray()->At(i)));
1104 tpcrow->SetN2(clrow->GetArray()->GetEntriesFast());
1105 tpcrow->SetClusters2(new AliTPCclusterMI[tpcrow->GetN2()]);
1106 for (Int_t i=0;i<tpcrow->GetN2();i++)
1107 tpcrow->SetCluster2(i,*(AliTPCclusterMI*)(clrow->GetArray()->At(i)));
1118 Int_t AliTPCtrackerMI::LoadClusters()
1121 // load clusters to the memory
1122 AliTPCClustersRow *clrow= new AliTPCClustersRow;
1123 clrow->SetClass("AliTPCclusterMI");
1125 clrow->GetArray()->ExpandCreateFast(10000);
1127 // TTree * tree = fClustersArray.GetTree();
1129 TTree * tree = fInput;
1130 TBranch * br = tree->GetBranch("Segment");
1131 br->SetAddress(&clrow);
1133 Int_t j=Int_t(tree->GetEntries());
1134 for (Int_t i=0; i<j; i++) {
1138 fParam->AdjustSectorRow(clrow->GetID(),sec,row);
1139 for (Int_t icl=0; icl<clrow->GetArray()->GetEntriesFast(); icl++){
1140 Transform((AliTPCclusterMI*)(clrow->GetArray()->At(icl)));
1143 AliTPCtrackerRow * tpcrow=0;
1146 tpcrow = &(fInnerSec[sec%fkNIS][row]);
1150 tpcrow = &(fOuterSec[(sec-fkNIS*2)%fkNOS][row]);
1151 left = (sec-fkNIS*2)/fkNOS;
1154 tpcrow->SetN1(clrow->GetArray()->GetEntriesFast());
1155 tpcrow->SetClusters1(new AliTPCclusterMI[tpcrow->GetN1()]);
1156 for (Int_t i=0;i<tpcrow->GetN1();i++)
1157 tpcrow->SetCluster1(i, *(AliTPCclusterMI*)(clrow->GetArray()->At(i)));
1160 tpcrow->SetN2(clrow->GetArray()->GetEntriesFast());
1161 tpcrow->SetClusters2(new AliTPCclusterMI[tpcrow->GetN2()]);
1162 for (Int_t i=0;i<tpcrow->GetN2();i++)
1163 tpcrow->SetCluster2(i,*(AliTPCclusterMI*)(clrow->GetArray()->At(i)));
1174 void AliTPCtrackerMI::UnloadClusters()
1177 // unload clusters from the memory
1179 Int_t nrows = fOuterSec->GetNRows();
1180 for (Int_t sec = 0;sec<fkNOS;sec++)
1181 for (Int_t row = 0;row<nrows;row++){
1182 AliTPCtrackerRow* tpcrow = &(fOuterSec[sec%fkNOS][row]);
1184 // if (tpcrow->fClusters1) delete []tpcrow->fClusters1;
1185 // if (tpcrow->fClusters2) delete []tpcrow->fClusters2;
1187 tpcrow->ResetClusters();
1190 nrows = fInnerSec->GetNRows();
1191 for (Int_t sec = 0;sec<fkNIS;sec++)
1192 for (Int_t row = 0;row<nrows;row++){
1193 AliTPCtrackerRow* tpcrow = &(fInnerSec[sec%fkNIS][row]);
1195 // if (tpcrow->fClusters1) delete []tpcrow->fClusters1;
1196 //if (tpcrow->fClusters2) delete []tpcrow->fClusters2;
1198 tpcrow->ResetClusters();
1204 void AliTPCtrackerMI::FillClusterArray(TObjArray* array) const{
1206 // Filling cluster to the array - For visualization purposes
1209 nrows = fOuterSec->GetNRows();
1210 for (Int_t sec = 0;sec<fkNOS;sec++)
1211 for (Int_t row = 0;row<nrows;row++){
1212 AliTPCtrackerRow* tpcrow = &(fOuterSec[sec%fkNOS][row]);
1213 if (!tpcrow) continue;
1214 for (Int_t icl = 0;icl<tpcrow->GetN();icl++){
1215 array->AddLast((TObject*)((*tpcrow)[icl]));
1218 nrows = fInnerSec->GetNRows();
1219 for (Int_t sec = 0;sec<fkNIS;sec++)
1220 for (Int_t row = 0;row<nrows;row++){
1221 AliTPCtrackerRow* tpcrow = &(fInnerSec[sec%fkNIS][row]);
1222 if (!tpcrow) continue;
1223 for (Int_t icl = 0;icl<tpcrow->GetN();icl++){
1224 array->AddLast((TObject*)(*tpcrow)[icl]);
1230 void AliTPCtrackerMI::Transform(AliTPCclusterMI * cluster){
1234 AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ;
1236 AliFatal("Tranformations not in calibDB");
1238 Double_t x[3]={cluster->GetRow(),cluster->GetPad(),cluster->GetTimeBin()};
1239 Int_t i[1]={cluster->GetDetector()};
1240 transform->Transform(x,i,0,1);
1241 if (!AliTPCReconstructor::GetRecoParam()->GetBYMirror()){
1242 if (cluster->GetDetector()%36>17){
1248 // in debug mode check the transformation
1250 if (AliTPCReconstructor::StreamLevel()>1) {
1252 cluster->GetGlobalXYZ(gx);
1253 Int_t event = (fEvent==NULL)? 0: fEvent->GetEventNumberInFile();
1254 TTreeSRedirector &cstream = *fDebugStreamer;
1255 cstream<<"Transform"<<
1266 cluster->SetX(x[0]);
1267 cluster->SetY(x[1]);
1268 cluster->SetZ(x[2]);
1274 //if (!fParam->IsGeoRead()) fParam->ReadGeoMatrices();
1275 TGeoHMatrix *mat = fParam->GetClusterMatrix(cluster->GetDetector());
1277 Double_t pos[3]= {cluster->GetX(),cluster->GetY(),cluster->GetZ()};
1278 Double_t posC[3]={cluster->GetX(),cluster->GetY(),cluster->GetZ()};
1279 if (mat) mat->LocalToMaster(pos,posC);
1281 // chack Loading of Geo matrices from GeoManager - TEMPORARY FIX
1283 cluster->SetX(posC[0]);
1284 cluster->SetY(posC[1]);
1285 cluster->SetZ(posC[2]);
1288 //_____________________________________________________________________________
1289 Int_t AliTPCtrackerMI::LoadOuterSectors() {
1290 //-----------------------------------------------------------------
1291 // This function fills outer TPC sectors with clusters.
1292 //-----------------------------------------------------------------
1293 Int_t nrows = fOuterSec->GetNRows();
1295 for (Int_t sec = 0;sec<fkNOS;sec++)
1296 for (Int_t row = 0;row<nrows;row++){
1297 AliTPCtrackerRow* tpcrow = &(fOuterSec[sec%fkNOS][row]);
1298 Int_t sec2 = sec+2*fkNIS;
1300 Int_t ncl = tpcrow->GetN1();
1302 AliTPCclusterMI *c= (tpcrow->GetCluster1(ncl));
1303 index=(((sec2<<8)+row)<<16)+ncl;
1304 tpcrow->InsertCluster(c,index);
1307 ncl = tpcrow->GetN2();
1309 AliTPCclusterMI *c= (tpcrow->GetCluster2(ncl));
1310 index=((((sec2+fkNOS)<<8)+row)<<16)+ncl;
1311 tpcrow->InsertCluster(c,index);
1314 // write indexes for fast acces
1316 for (Int_t i=0;i<510;i++)
1317 tpcrow->SetFastCluster(i,-1);
1318 for (Int_t i=0;i<tpcrow->GetN();i++){
1319 Int_t zi = Int_t((*tpcrow)[i]->GetZ()+255.);
1320 tpcrow->SetFastCluster(zi,i); // write index
1323 for (Int_t i=0;i<510;i++){
1324 if (tpcrow->GetFastCluster(i)<0)
1325 tpcrow->SetFastCluster(i,last);
1327 last = tpcrow->GetFastCluster(i);
1336 //_____________________________________________________________________________
1337 Int_t AliTPCtrackerMI::LoadInnerSectors() {
1338 //-----------------------------------------------------------------
1339 // This function fills inner TPC sectors with clusters.
1340 //-----------------------------------------------------------------
1341 Int_t nrows = fInnerSec->GetNRows();
1343 for (Int_t sec = 0;sec<fkNIS;sec++)
1344 for (Int_t row = 0;row<nrows;row++){
1345 AliTPCtrackerRow* tpcrow = &(fInnerSec[sec%fkNIS][row]);
1348 Int_t ncl = tpcrow->GetN1();
1350 AliTPCclusterMI *c= (tpcrow->GetCluster1(ncl));
1351 index=(((sec<<8)+row)<<16)+ncl;
1352 tpcrow->InsertCluster(c,index);
1355 ncl = tpcrow->GetN2();
1357 AliTPCclusterMI *c= (tpcrow->GetCluster2(ncl));
1358 index=((((sec+fkNIS)<<8)+row)<<16)+ncl;
1359 tpcrow->InsertCluster(c,index);
1362 // write indexes for fast acces
1364 for (Int_t i=0;i<510;i++)
1365 tpcrow->SetFastCluster(i,-1);
1366 for (Int_t i=0;i<tpcrow->GetN();i++){
1367 Int_t zi = Int_t((*tpcrow)[i]->GetZ()+255.);
1368 tpcrow->SetFastCluster(zi,i); // write index
1371 for (Int_t i=0;i<510;i++){
1372 if (tpcrow->GetFastCluster(i)<0)
1373 tpcrow->SetFastCluster(i,last);
1375 last = tpcrow->GetFastCluster(i);
1387 //_________________________________________________________________________
1388 AliTPCclusterMI *AliTPCtrackerMI::GetClusterMI(Int_t index) const {
1389 //--------------------------------------------------------------------
1390 // Return pointer to a given cluster
1391 //--------------------------------------------------------------------
1392 if (index<0) return 0; // no cluster
1393 Int_t sec=(index&0xff000000)>>24;
1394 Int_t row=(index&0x00ff0000)>>16;
1395 Int_t ncl=(index&0x00007fff)>>00;
1397 const AliTPCtrackerRow * tpcrow=0;
1398 AliTPCclusterMI * clrow =0;
1400 if (sec<0 || sec>=fkNIS*4) {
1401 AliWarning(Form("Wrong sector %d",sec));
1406 tpcrow = &(fInnerSec[sec%fkNIS][row]);
1407 if (tpcrow==0) return 0;
1410 if (tpcrow->GetN1()<=ncl) return 0;
1411 clrow = tpcrow->GetClusters1();
1414 if (tpcrow->GetN2()<=ncl) return 0;
1415 clrow = tpcrow->GetClusters2();
1419 tpcrow = &(fOuterSec[(sec-fkNIS*2)%fkNOS][row]);
1420 if (tpcrow==0) return 0;
1422 if (sec-2*fkNIS<fkNOS) {
1423 if (tpcrow->GetN1()<=ncl) return 0;
1424 clrow = tpcrow->GetClusters1();
1427 if (tpcrow->GetN2()<=ncl) return 0;
1428 clrow = tpcrow->GetClusters2();
1432 return &(clrow[ncl]);
1438 Int_t AliTPCtrackerMI::FollowToNext(AliTPCseed& t, Int_t nr) {
1439 //-----------------------------------------------------------------
1440 // This function tries to find a track prolongation to next pad row
1441 //-----------------------------------------------------------------
1443 Double_t x= GetXrow(nr), ymax=GetMaxY(nr);
1444 AliTPCclusterMI *cl=0;
1445 Int_t tpcindex= t.GetClusterIndex2(nr);
1447 // update current shape info every 5 pad-row
1448 // if ( (nr%5==0) || t.GetNumberOfClusters()<2 || (t.fCurrentSigmaY2<0.0001) ){
1452 if (fIteration>0 && tpcindex>=-1){ //if we have already clusters
1454 if (tpcindex==-1) return 0; //track in dead zone
1456 cl = t.GetClusterPointer(nr);
1457 if ( (cl==0) ) cl = GetClusterMI(tpcindex);
1458 t.SetCurrentClusterIndex1(tpcindex);
1461 Int_t relativesector = ((tpcindex&0xff000000)>>24)%18; // if previously accepted cluster in different sector
1462 Float_t angle = relativesector*fSectors->GetAlpha()+fSectors->GetAlphaShift();
1464 if (angle<-TMath::Pi()) angle += 2*TMath::Pi();
1465 if (angle>=TMath::Pi()) angle -= 2*TMath::Pi();
1467 if (TMath::Abs(angle-t.GetAlpha())>0.001){
1468 Double_t rotation = angle-t.GetAlpha();
1469 t.SetRelativeSector(relativesector);
1470 if (!t.Rotate(rotation)) return 0;
1472 if (!t.PropagateTo(x)) return 0;
1474 t.SetCurrentCluster(cl);
1476 Int_t accept = AcceptCluster(&t,t.GetCurrentCluster());
1477 if ((tpcindex&0x8000)==0) accept =0;
1479 //if founded cluster is acceptible
1480 if (cl->IsUsed(11)) { // id cluster is shared inrease uncertainty
1481 t.SetErrorY2(t.GetErrorY2()+0.03);
1482 t.SetErrorZ2(t.GetErrorZ2()+0.03);
1483 t.SetErrorY2(t.GetErrorY2()*3);
1484 t.SetErrorZ2(t.GetErrorZ2()*3);
1486 t.SetNFoundable(t.GetNFoundable()+1);
1487 UpdateTrack(&t,accept);
1492 if (TMath::Abs(t.GetSnp())>AliTPCReconstructor::GetMaxSnpTracker()) return 0; // cut on angle
1494 // not look for new cluster during refitting
1495 t.SetNFoundable(t.GetNFoundable()+1);
1500 // if (TMath::Abs(t.GetSnp())>0.95 || TMath::Abs(x*t.GetC()-t.GetEta())>0.95) return 0;// patch 28 fev 06
1501 Double_t y=t.GetYat(x);
1502 if (TMath::Abs(y)>ymax){
1504 t.SetRelativeSector((t.GetRelativeSector()+1) % fN);
1505 if (!t.Rotate(fSectors->GetAlpha()))
1507 } else if (y <-ymax) {
1508 t.SetRelativeSector((t.GetRelativeSector()-1+fN) % fN);
1509 if (!t.Rotate(-fSectors->GetAlpha()))
1515 if (!t.PropagateTo(x)) {
1516 if (fIteration==0) t.SetRemoval(10);
1520 Double_t z=t.GetZ();
1523 if (!IsActive(t.GetRelativeSector(),nr)) {
1525 t.SetClusterIndex2(nr,-1);
1528 //AliInfo(Form("A - Sector%d phi %f - alpha %f", t.fRelativeSector,y/x, t.GetAlpha()));
1529 Bool_t isActive = IsActive(t.GetRelativeSector(),nr);
1530 Bool_t isActive2 = (nr>=fInnerSec->GetNRows()) ? fOuterSec[t.GetRelativeSector()][nr-fInnerSec->GetNRows()].GetN()>0:fInnerSec[t.GetRelativeSector()][nr].GetN()>0;
1532 if (!isActive || !isActive2) return 0;
1534 const AliTPCtrackerRow &krow=GetRow(t.GetRelativeSector(),nr);
1535 if ( (t.GetSigmaY2()<0) || t.GetSigmaZ2()<0) return 0;
1537 Double_t roadz = 1.;
1539 if (TMath::Abs(TMath::Abs(y)-ymax)<krow.GetDeadZone()){
1541 t.SetClusterIndex2(nr,-1);
1546 if (TMath::Abs(z)<(AliTPCReconstructor::GetCtgRange()*x+10) && TMath::Abs(z)<fParam->GetZLength(0) && (TMath::Abs(t.GetSnp())<AliTPCReconstructor::GetMaxSnpTracker()))
1547 t.SetNFoundable(t.GetNFoundable()+1);
1553 // cl = krow.FindNearest2(y+10.,z,roady,roadz,index);
1554 cl = krow.FindNearest2(y,z,roady,roadz,index);
1555 if (cl) t.SetCurrentClusterIndex1(krow.GetIndex(index));
1558 t.SetCurrentCluster(cl);
1560 if (fIteration==2&&cl->IsUsed(10)) return 0;
1561 Int_t accept = AcceptCluster(&t,t.GetCurrentCluster());
1562 if (fIteration==2&&cl->IsUsed(11)) {
1563 t.SetErrorY2(t.GetErrorY2()+0.03);
1564 t.SetErrorZ2(t.GetErrorZ2()+0.03);
1565 t.SetErrorY2(t.GetErrorY2()*3);
1566 t.SetErrorZ2(t.GetErrorZ2()*3);
1569 if (t.fCurrentCluster->IsUsed(10)){
1574 if (t.fNShared>0.7*t.GetNumberOfClusters()) {
1580 if (accept<3) UpdateTrack(&t,accept);
1583 if ( fIteration==0 && t.GetNFoundable()*0.5 > t.GetNumberOfClusters()) t.SetRemoval(10);
1589 Int_t AliTPCtrackerMI::FollowToNextFast(AliTPCseed& t, Int_t nr) {
1590 //-----------------------------------------------------------------
1591 // This function tries to find a track prolongation to next pad row
1592 //-----------------------------------------------------------------
1594 Double_t x= GetXrow(nr), ymax=GetMaxY(nr);
1596 if (!t.GetProlongation(x,y,z)) {
1602 if (TMath::Abs(y)>ymax){
1605 t.SetRelativeSector((t.GetRelativeSector()+1) % fN);
1606 if (!t.Rotate(fSectors->GetAlpha()))
1608 } else if (y <-ymax) {
1609 t.SetRelativeSector((t.GetRelativeSector()-1+fN) % fN);
1610 if (!t.Rotate(-fSectors->GetAlpha()))
1613 if (!t.PropagateTo(x)) {
1616 t.GetProlongation(x,y,z);
1619 // update current shape info every 2 pad-row
1620 if ( (nr%2==0) || t.GetNumberOfClusters()<2 || (t.GetCurrentSigmaY2()<0.0001) ){
1621 // t.fCurrentSigmaY = GetSigmaY(&t);
1622 //t.fCurrentSigmaZ = GetSigmaZ(&t);
1626 AliTPCclusterMI *cl=0;
1631 const AliTPCtrackerRow &krow=GetRow(t.GetRelativeSector(),nr);
1632 if ( (t.GetSigmaY2()<0) || t.GetSigmaZ2()<0) return 0;
1634 Double_t roadz = 1.;
1637 if (TMath::Abs(TMath::Abs(y)-ymax)<krow.GetDeadZone()){
1639 t.SetClusterIndex2(row,-1);
1644 if (TMath::Abs(z)>(AliTPCReconstructor::GetCtgRange()*x+10)) t.SetClusterIndex2(row,-1);
1648 if ((cl==0)&&(krow)) {
1649 // cl = krow.FindNearest2(y+10,z,roady,roadz,index);
1650 cl = krow.FindNearest2(y,z,roady,roadz,index);
1652 if (cl) t.SetCurrentClusterIndex1(krow.GetIndex(index));
1656 t.SetCurrentCluster(cl);
1657 // Int_t accept = AcceptCluster(&t,t.fCurrentCluster);
1659 t.SetClusterIndex2(row,index);
1660 t.SetClusterPointer(row, cl);
1668 //_________________________________________________________________________
1669 Bool_t AliTPCtrackerMI::GetTrackPoint(Int_t index, AliTrackPoint &p ) const
1671 // Get track space point by index
1672 // return false in case the cluster doesn't exist
1673 AliTPCclusterMI *cl = GetClusterMI(index);
1674 if (!cl) return kFALSE;
1675 Int_t sector = (index&0xff000000)>>24;
1676 // Int_t row = (index&0x00ff0000)>>16;
1678 // xyz[0] = fParam->GetPadRowRadii(sector,row);
1679 xyz[0] = cl->GetX();
1680 xyz[1] = cl->GetY();
1681 xyz[2] = cl->GetZ();
1683 fParam->AdjustCosSin(sector,cos,sin);
1684 Float_t x = cos*xyz[0]-sin*xyz[1];
1685 Float_t y = cos*xyz[1]+sin*xyz[0];
1687 Float_t sigmaY2 = 0.027*cl->GetSigmaY2();
1688 if (sector < fParam->GetNInnerSector()) sigmaY2 *= 2.07;
1689 Float_t sigmaZ2 = 0.066*cl->GetSigmaZ2();
1690 if (sector < fParam->GetNInnerSector()) sigmaZ2 *= 1.77;
1691 cov[0] = sin*sin*sigmaY2;
1692 cov[1] = -sin*cos*sigmaY2;
1694 cov[3] = cos*cos*sigmaY2;
1697 p.SetXYZ(x,y,xyz[2],cov);
1698 AliGeomManager::ELayerID iLayer;
1700 if (sector < fParam->GetNInnerSector()) {
1701 iLayer = AliGeomManager::kTPC1;
1705 iLayer = AliGeomManager::kTPC2;
1706 idet = sector - fParam->GetNInnerSector();
1708 UShort_t volid = AliGeomManager::LayerToVolUID(iLayer,idet);
1709 p.SetVolumeID(volid);
1715 Int_t AliTPCtrackerMI::UpdateClusters(AliTPCseed& t, Int_t nr) {
1716 //-----------------------------------------------------------------
1717 // This function tries to find a track prolongation to next pad row
1718 //-----------------------------------------------------------------
1719 t.SetCurrentCluster(0);
1720 t.SetCurrentClusterIndex1(0);
1722 Double_t xt=t.GetX();
1723 Int_t row = GetRowNumber(xt)-1;
1724 Double_t ymax= GetMaxY(nr);
1726 if (row < nr) return 1; // don't prolongate if not information until now -
1727 // if (TMath::Abs(t.GetSnp())>0.9 && t.GetNumberOfClusters()>40. && fIteration!=2) {
1729 // return 0; // not prolongate strongly inclined tracks
1731 // if (TMath::Abs(t.GetSnp())>0.95) {
1733 // return 0; // not prolongate strongly inclined tracks
1734 // }// patch 28 fev 06
1736 Double_t x= GetXrow(nr);
1738 //t.PropagateTo(x+0.02);
1739 //t.PropagateTo(x+0.01);
1740 if (!t.PropagateTo(x)){
1747 if (TMath::Abs(y)>ymax){
1749 t.SetRelativeSector((t.GetRelativeSector()+1) % fN);
1750 if (!t.Rotate(fSectors->GetAlpha()))
1752 } else if (y <-ymax) {
1753 t.SetRelativeSector((t.GetRelativeSector()-1+fN) % fN);
1754 if (!t.Rotate(-fSectors->GetAlpha()))
1757 // if (!t.PropagateTo(x)){
1764 if (TMath::Abs(t.GetSnp())>AliTPCReconstructor::GetMaxSnpTracker()) return 0;
1766 if (!IsActive(t.GetRelativeSector(),nr)) {
1768 t.SetClusterIndex2(nr,-1);
1771 //AliInfo(Form("A - Sector%d phi %f - alpha %f", t.fRelativeSector,y/x, t.GetAlpha()));
1773 AliTPCtrackerRow &krow=GetRow(t.GetRelativeSector(),nr);
1775 if (TMath::Abs(TMath::Abs(y)-ymax)<krow.GetDeadZone()){
1777 t.SetClusterIndex2(nr,-1);
1782 if (TMath::Abs(t.GetZ())<(AliTPCReconstructor::GetCtgRange()*t.GetX()+10) && (TMath::Abs(t.GetSnp())<AliTPCReconstructor::GetMaxSnpTracker())) t.SetNFoundable(t.GetNFoundable()+1);
1788 if ( (nr%2==0) || t.GetNumberOfClusters()<2){
1789 // t.fCurrentSigmaY = GetSigmaY(&t);
1790 //t.fCurrentSigmaZ = GetSigmaZ(&t);
1794 AliTPCclusterMI *cl=0;
1797 Double_t roady = 1.;
1798 Double_t roadz = 1.;
1802 index = t.GetClusterIndex2(nr);
1803 if ( (index>0) && (index&0x8000)==0){
1804 cl = t.GetClusterPointer(nr);
1805 if ( (cl==0) && (index>0)) cl = GetClusterMI(index);
1806 t.SetCurrentClusterIndex1(index);
1808 t.SetCurrentCluster(cl);
1814 // if (index<0) return 0;
1815 UInt_t uindex = TMath::Abs(index);
1818 //cl = krow.FindNearest2(y+10,z,roady,roadz,uindex);
1819 cl = krow.FindNearest2(y,z,roady,roadz,uindex);
1822 if (cl) t.SetCurrentClusterIndex1(krow.GetIndex(uindex));
1823 t.SetCurrentCluster(cl);
1829 Int_t AliTPCtrackerMI::FollowToNextCluster(AliTPCseed & t, Int_t nr) {
1830 //-----------------------------------------------------------------
1831 // This function tries to find a track prolongation to next pad row
1832 //-----------------------------------------------------------------
1834 //update error according neighborhoud
1836 if (t.GetCurrentCluster()) {
1838 Int_t accept = AcceptCluster(&t,t.GetCurrentCluster());
1840 if (t.GetCurrentCluster()->IsUsed(10)){
1845 t.SetNShared(t.GetNShared()+1);
1846 if (t.GetNShared()>0.7*t.GetNumberOfClusters()) {
1851 if (fIteration>0) accept = 0;
1852 if (accept<3) UpdateTrack(&t,accept);
1856 if ( ( (t.GetSigmaY2()+t.GetSigmaZ2())>0.16)&& t.GetNumberOfClusters()>18) t.SetRemoval(10);
1857 if ( t.GetChi2()/t.GetNumberOfClusters()>6 &&t.GetNumberOfClusters()>18) t.SetRemoval(10);
1859 if (( (t.GetNFoundable()*0.5 > t.GetNumberOfClusters()) || t.GetNoCluster()>15)) t.SetRemoval(10);
1867 //_____________________________________________________________________________
1868 Int_t AliTPCtrackerMI::FollowProlongation(AliTPCseed& t, Int_t rf, Int_t step) {
1869 //-----------------------------------------------------------------
1870 // This function tries to find a track prolongation.
1871 //-----------------------------------------------------------------
1872 Double_t xt=t.GetX();
1874 Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
1875 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
1876 if (alpha < 0. ) alpha += 2.*TMath::Pi();
1878 t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN);
1880 Int_t first = GetRowNumber(xt)-1;
1881 for (Int_t nr= first; nr>=rf; nr-=step) {
1883 if (t.GetKinkIndexes()[0]>0){
1884 for (Int_t i=0;i<3;i++){
1885 Int_t index = t.GetKinkIndexes()[i];
1886 if (index==0) break;
1887 if (index<0) continue;
1889 AliKink * kink = (AliKink*)fEvent->GetKink(index-1);
1891 printf("PROBLEM\n");
1894 Int_t kinkrow = kink->GetTPCRow0()+2+Int_t(0.5/(0.05+kink->GetAngle(2)));
1896 AliExternalTrackParam paramd(t);
1897 kink->SetDaughter(paramd);
1898 kink->SetStatus(2,5);
1905 if (nr==80) t.UpdateReference();
1906 if (nr<fInnerSec->GetNRows())
1907 fSectors = fInnerSec;
1909 fSectors = fOuterSec;
1910 if (FollowToNext(t,nr)==0)
1919 //_____________________________________________________________________________
1920 Int_t AliTPCtrackerMI::FollowProlongationFast(AliTPCseed& t, Int_t rf, Int_t step) {
1921 //-----------------------------------------------------------------
1922 // This function tries to find a track prolongation.
1923 //-----------------------------------------------------------------
1924 Double_t xt=t.GetX();
1926 Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
1927 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
1928 if (alpha < 0. ) alpha += 2.*TMath::Pi();
1929 t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN);
1931 for (Int_t nr=GetRowNumber(xt)-1; nr>=rf; nr-=step) {
1933 if (FollowToNextFast(t,nr)==0)
1934 if (!t.IsActive()) return 0;
1944 Int_t AliTPCtrackerMI::FollowBackProlongation(AliTPCseed& t, Int_t rf) {
1945 //-----------------------------------------------------------------
1946 // This function tries to find a track prolongation.
1947 //-----------------------------------------------------------------
1949 Double_t xt=t.GetX();
1950 Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
1951 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
1952 if (alpha < 0. ) alpha += 2.*TMath::Pi();
1953 t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN);
1955 Int_t first = t.GetFirstPoint();
1956 if (first<GetRowNumber(xt)+1) first = GetRowNumber(xt)+1;
1958 if (first<0) first=0;
1959 for (Int_t nr=first; nr<=rf; nr++) {
1960 // if ( (TMath::Abs(t.GetSnp())>0.95)) break;//patch 28 fev 06
1961 if (t.GetKinkIndexes()[0]<0){
1962 for (Int_t i=0;i<3;i++){
1963 Int_t index = t.GetKinkIndexes()[i];
1964 if (index==0) break;
1965 if (index>0) continue;
1966 index = TMath::Abs(index);
1967 AliKink * kink = (AliKink*)fEvent->GetKink(index-1);
1969 printf("PROBLEM\n");
1972 Int_t kinkrow = kink->GetTPCRow0()-2-Int_t(0.5/(0.05+kink->GetAngle(2)));
1974 AliExternalTrackParam paramm(t);
1975 kink->SetMother(paramm);
1976 kink->SetStatus(2,1);
1983 if (nr<fInnerSec->GetNRows())
1984 fSectors = fInnerSec;
1986 fSectors = fOuterSec;
1997 Float_t AliTPCtrackerMI::OverlapFactor(AliTPCseed * s1, AliTPCseed * s2, Int_t &sum1, Int_t & sum2)
2005 Float_t dz2 =(s1->GetZ() - s2->GetZ());
2008 Float_t dy2 =TMath::Abs((s1->GetY() - s2->GetY()));
2010 Float_t distance = TMath::Sqrt(dz2+dy2);
2011 if (distance>4.) return 0; // if there are far away - not overlap - to reduce combinatorics
2014 Int_t firstpoint = TMath::Min(s1->GetFirstPoint(),s2->GetFirstPoint());
2015 Int_t lastpoint = TMath::Max(s1->GetLastPoint(),s2->GetLastPoint());
2020 if (firstpoint>lastpoint) {
2021 firstpoint =lastpoint;
2026 for (Int_t i=firstpoint-1;i<lastpoint+1;i++){
2027 if (s1->GetClusterIndex2(i)>0) sum1++;
2028 if (s2->GetClusterIndex2(i)>0) sum2++;
2029 if (s1->GetClusterIndex2(i)==s2->GetClusterIndex2(i) && s1->GetClusterIndex2(i)>0) {
2033 if (sum<5) return 0;
2035 Float_t summin = TMath::Min(sum1+1,sum2+1);
2036 Float_t ratio = (sum+1)/Float_t(summin);
2040 void AliTPCtrackerMI::SignShared(AliTPCseed * s1, AliTPCseed * s2)
2044 Float_t thetaCut = 0.2;//+10.*TMath::Sqrt(s1->GetSigmaTglZ()+ s2->GetSigmaTglZ());
2045 if (TMath::Abs(s1->GetTgl()-s2->GetTgl())>thetaCut) return;
2046 Float_t minCl = TMath::Min(s1->GetNumberOfClusters(),s2->GetNumberOfClusters());
2047 Int_t cutN0 = TMath::Max(5,TMath::Nint(0.1*minCl));
2052 //Int_t firstpoint = TMath::Max(s1->GetFirstPoint(),s2->GetFirstPoint());
2053 //Int_t lastpoint = TMath::Min(s1->GetLastPoint(),s2->GetLastPoint());
2054 Int_t firstpoint = 0;
2055 Int_t lastpoint = 160;
2057 // if (firstpoint>=lastpoint-5) return;;
2059 for (Int_t i=firstpoint;i<lastpoint;i++){
2060 // if ( (s1->GetClusterIndex2(i)&0xFFFF8FFF)==(s2->GetClusterIndex2(i)&0xFFFF8FFF) && s1->GetClusterIndex2(i)>0) {
2061 if ( (s1->GetClusterIndex2(i))==(s2->GetClusterIndex2(i)) && s1->GetClusterIndex2(i)>0) {
2063 s1->SetSharedMapBit(i, kTRUE);
2064 s2->SetSharedMapBit(i, kTRUE);
2066 if (s1->GetClusterIndex2(i)>0)
2067 s1->SetClusterMapBit(i, kTRUE);
2069 if (sumshared>cutN0){
2072 for (Int_t i=firstpoint;i<lastpoint;i++){
2073 // if ( (s1->GetClusterIndex2(i)&0xFFFF8FFF)==(s2->GetClusterIndex2(i)&0xFFFF8FFF) && s1->GetClusterIndex2(i)>0) {
2074 if ( (s1->GetClusterIndex2(i))==(s2->GetClusterIndex2(i)) && s1->GetClusterIndex2(i)>0) {
2075 AliTPCTrackerPoint *p1 = s1->GetTrackPoint(i);
2076 AliTPCTrackerPoint *p2 = s2->GetTrackPoint(i);;
2077 if (s1->IsActive()&&s2->IsActive()){
2078 p1->SetShared(kTRUE);
2079 p2->SetShared(kTRUE);
2085 if (sumshared>cutN0){
2086 for (Int_t i=0;i<4;i++){
2087 if (s1->GetOverlapLabel(3*i)==0){
2088 s1->SetOverlapLabel(3*i, s2->GetLabel());
2089 s1->SetOverlapLabel(3*i+1,sumshared);
2090 s1->SetOverlapLabel(3*i+2,s2->GetUniqueID());
2094 for (Int_t i=0;i<4;i++){
2095 if (s2->GetOverlapLabel(3*i)==0){
2096 s2->SetOverlapLabel(3*i, s1->GetLabel());
2097 s2->SetOverlapLabel(3*i+1,sumshared);
2098 s2->SetOverlapLabel(3*i+2,s1->GetUniqueID());
2105 void AliTPCtrackerMI::SignShared(TObjArray * arr)
2108 //sort trackss according sectors
2110 for (Int_t i=0; i<arr->GetEntriesFast(); i++) {
2111 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2113 //if (pt) RotateToLocal(pt);
2117 arr->Sort(); // sorting according relative sectors
2118 arr->Expand(arr->GetEntries());
2121 Int_t nseed=arr->GetEntriesFast();
2122 for (Int_t i=0; i<nseed; i++) {
2123 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2125 for (Int_t j=0;j<=12;j++){
2126 pt->SetOverlapLabel(j,0);
2129 for (Int_t i=0; i<nseed; i++) {
2130 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2132 if (pt->GetRemoval()>10) continue;
2133 for (Int_t j=i+1; j<nseed; j++){
2134 AliTPCseed *pt2=(AliTPCseed*)arr->UncheckedAt(j);
2135 if (TMath::Abs(pt->GetRelativeSector()-pt2->GetRelativeSector())>1) continue;
2137 if (pt2->GetRemoval()<=10) {
2138 //if ( TMath::Abs(pt->GetRelativeSector()-pt2->GetRelativeSector())>0) break;
2146 void AliTPCtrackerMI::SortTracks(TObjArray * arr, Int_t mode) const
2149 //sort tracks in array according mode criteria
2150 Int_t nseed = arr->GetEntriesFast();
2151 for (Int_t i=0; i<nseed; i++) {
2152 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2163 void AliTPCtrackerMI::RemoveUsed2(TObjArray * arr, Float_t factor1, Float_t factor2, Int_t minimal)
2166 // Loop over all tracks and remove overlaped tracks (with lower quality)
2168 // 1. Unsign clusters
2169 // 2. Sort tracks according quality
2170 // Quality is defined by the number of cluster between first and last points
2172 // 3. Loop over tracks - decreasing quality order
2173 // a.) remove - If the fraction of shared cluster less than factor (1- n or 2)
2174 // b.) remove - If the minimal number of clusters less than minimal and not ITS
2175 // c.) if track accepted - sign clusters
2177 //Called in - AliTPCtrackerMI::Clusters2Tracks()
2178 // - AliTPCtrackerMI::PropagateBack()
2179 // - AliTPCtrackerMI::RefitInward()
2185 Int_t nseed = arr->GetEntriesFast();
2186 Float_t * quality = new Float_t[nseed];
2187 Int_t * indexes = new Int_t[nseed];
2191 for (Int_t i=0; i<nseed; i++) {
2192 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2197 pt->UpdatePoints(); //select first last max dens points
2198 Float_t * points = pt->GetPoints();
2199 if (points[3]<0.8) quality[i] =-1;
2200 quality[i] = (points[2]-points[0])+pt->GetNumberOfClusters();
2201 //prefer high momenta tracks if overlaps
2202 quality[i] *= TMath::Sqrt(TMath::Abs(pt->Pt())+0.5);
2204 TMath::Sort(nseed,quality,indexes);
2207 for (Int_t itrack=0; itrack<nseed; itrack++) {
2208 Int_t trackindex = indexes[itrack];
2209 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(trackindex);
2212 if (quality[trackindex]<0){
2214 delete arr->RemoveAt(trackindex);
2217 arr->RemoveAt(trackindex);
2223 Int_t first = Int_t(pt->GetPoints()[0]);
2224 Int_t last = Int_t(pt->GetPoints()[2]);
2225 Double_t factor = (pt->GetBConstrain()) ? factor1: factor2;
2227 Int_t found,foundable,shared;
2228 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
2229 // pt->GetClusterStatistic(0,160, found, foundable,shared,kFALSE);
2230 Bool_t itsgold =kFALSE;
2233 if (pt->GetESD()->GetITSclusters(dummy)>4) itsgold= kTRUE;
2237 if (Float_t(shared+1)/Float_t(found+1)>factor){
2238 if (pt->GetKinkIndexes()[0]!=0) continue; //don't remove tracks - part of the kinks
2239 delete arr->RemoveAt(trackindex);
2242 if (pt->GetNumberOfClusters()<50&&(found-0.5*shared)<minimal){ //remove short tracks
2243 if (pt->GetKinkIndexes()[0]!=0) continue; //don't remove tracks - part of the kinks
2244 delete arr->RemoveAt(trackindex);
2250 //if (sharedfactor>0.4) continue;
2251 if (pt->GetKinkIndexes()[0]>0) continue;
2252 //Remove tracks with undefined properties - seems
2253 if (pt->GetSigmaY2()<kAlmost0) continue; // ? what is the origin ?
2255 for (Int_t i=first; i<last; i++) {
2256 Int_t index=pt->GetClusterIndex2(i);
2257 // if (index<0 || index&0x8000 ) continue;
2258 if (index<0 || index&0x8000 ) continue;
2259 AliTPCclusterMI *c= pt->GetClusterPointer(i);
2266 Info("RemoveUsed","\n*****\nNumber of good tracks after shared removal\t%d\n",fNtracks);
2272 void AliTPCtrackerMI::UnsignClusters()
2275 // loop over all clusters and unsign them
2278 for (Int_t sec=0;sec<fkNIS;sec++){
2279 for (Int_t row=0;row<fInnerSec->GetNRows();row++){
2280 AliTPCclusterMI *cl = fInnerSec[sec][row].GetClusters1();
2281 for (Int_t icl =0;icl< fInnerSec[sec][row].GetN1();icl++)
2282 // if (cl[icl].IsUsed(10))
2284 cl = fInnerSec[sec][row].GetClusters2();
2285 for (Int_t icl =0;icl< fInnerSec[sec][row].GetN2();icl++)
2286 //if (cl[icl].IsUsed(10))
2291 for (Int_t sec=0;sec<fkNOS;sec++){
2292 for (Int_t row=0;row<fOuterSec->GetNRows();row++){
2293 AliTPCclusterMI *cl = fOuterSec[sec][row].GetClusters1();
2294 for (Int_t icl =0;icl< fOuterSec[sec][row].GetN1();icl++)
2295 //if (cl[icl].IsUsed(10))
2297 cl = fOuterSec[sec][row].GetClusters2();
2298 for (Int_t icl =0;icl< fOuterSec[sec][row].GetN2();icl++)
2299 //if (cl[icl].IsUsed(10))
2308 void AliTPCtrackerMI::SignClusters(TObjArray * arr, Float_t fnumber, Float_t fdensity)
2311 //sign clusters to be "used"
2313 // snumber and sdensity sign number of sigmas - bellow mean value to be accepted
2314 // loop over "primaries"
2328 Int_t nseed = arr->GetEntriesFast();
2329 for (Int_t i=0; i<nseed; i++) {
2330 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2334 if (!(pt->IsActive())) continue;
2335 Float_t dens = pt->GetNumberOfClusters()/Float_t(pt->GetNFoundable());
2336 if ( (dens>0.7) && (pt->GetNumberOfClusters()>70)){
2338 sumdens2+= dens*dens;
2339 sumn += pt->GetNumberOfClusters();
2340 sumn2 += pt->GetNumberOfClusters()*pt->GetNumberOfClusters();
2341 Float_t chi2 = pt->GetChi2()/pt->GetNumberOfClusters();
2344 sumchi2 +=chi2*chi2;
2349 Float_t mdensity = 0.9;
2350 Float_t meann = 130;
2351 Float_t meanchi = 1;
2352 Float_t sdensity = 0.1;
2353 Float_t smeann = 10;
2354 Float_t smeanchi =0.4;
2358 mdensity = sumdens/sum;
2360 meanchi = sumchi/sum;
2362 sdensity = sumdens2/sum-mdensity*mdensity;
2364 sdensity = TMath::Sqrt(sdensity);
2368 smeann = sumn2/sum-meann*meann;
2370 smeann = TMath::Sqrt(smeann);
2374 smeanchi = sumchi2/sum - meanchi*meanchi;
2376 smeanchi = TMath::Sqrt(smeanchi);
2382 //REMOVE SHORT DELTAS or tracks going out of sensitive volume of TPC
2384 for (Int_t i=0; i<nseed; i++) {
2385 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2389 if (pt->GetBSigned()) continue;
2390 if (pt->GetBConstrain()) continue;
2391 //if (!(pt->IsActive())) continue;
2393 Int_t found,foundable,shared;
2394 pt->GetClusterStatistic(0,160,found, foundable,shared);
2395 if (shared/float(found)>0.3) {
2396 if (shared/float(found)>0.9 ){
2397 //delete arr->RemoveAt(i);
2402 Bool_t isok =kFALSE;
2403 if ( (pt->GetNShared()/pt->GetNumberOfClusters()<0.5) &&pt->GetNumberOfClusters()>60)
2405 if ((TMath::Abs(1/pt->GetC())<100.) && (pt->GetNShared()/pt->GetNumberOfClusters()<0.7))
2407 if (TMath::Abs(pt->GetZ()/pt->GetX())>1.1)
2409 if ( (TMath::Abs(pt->GetSnp()>0.7) && pt->GetD(0,0)>60.))
2413 for (Int_t i=0; i<160; i++) {
2414 Int_t index=pt->GetClusterIndex2(i);
2415 if (index<0) continue;
2416 AliTPCclusterMI *c= pt->GetClusterPointer(i);
2418 //if (!(c->IsUsed(10))) c->Use();
2425 Double_t maxchi = meanchi+2.*smeanchi;
2427 for (Int_t i=0; i<nseed; i++) {
2428 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2432 //if (!(pt->IsActive())) continue;
2433 if (pt->GetBSigned()) continue;
2434 Double_t chi = pt->GetChi2()/pt->GetNumberOfClusters();
2435 if (chi>maxchi) continue;
2438 Float_t dens = pt->GetNumberOfClusters()/Float_t(pt->GetNFoundable());
2440 //sign only tracks with enoug big density at the beginning
2442 if ((pt->GetDensityFirst(40)<0.75) && pt->GetNumberOfClusters()<meann) continue;
2445 Double_t mindens = TMath::Max(double(mdensity-sdensity*fdensity*bfactor),0.65);
2446 Double_t minn = TMath::Max(Int_t(meann-fnumber*smeann*bfactor),50);
2448 // if (pt->fBConstrain) mindens = TMath::Max(mdensity-sdensity*fdensity*bfactor,0.65);
2449 if ( (pt->GetRemoval()==10) && (pt->GetSnp()>0.8)&&(dens>mindens))
2452 if ((dens>mindens && pt->GetNumberOfClusters()>minn) && chi<maxchi ){
2453 //Int_t noc=pt->GetNumberOfClusters();
2454 pt->SetBSigned(kTRUE);
2455 for (Int_t i=0; i<160; i++) {
2457 Int_t index=pt->GetClusterIndex2(i);
2458 if (index<0) continue;
2459 AliTPCclusterMI *c= pt->GetClusterPointer(i);
2461 // if (!(c->IsUsed(10))) c->Use();
2466 // gLastCheck = nseed;
2474 void AliTPCtrackerMI::StopNotActive(TObjArray * arr, Int_t row0, Float_t th0, Float_t th1, Float_t th2) const
2476 // stop not active tracks
2477 // take th1 as threshold for number of founded to number of foundable on last 10 active rows
2478 // take th2 as threshold for number of founded to number of foundable on last 20 active rows
2479 Int_t nseed = arr->GetEntriesFast();
2481 for (Int_t i=0; i<nseed; i++) {
2482 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2486 if (!(pt->IsActive())) continue;
2487 StopNotActive(pt,row0,th0, th1,th2);
2493 void AliTPCtrackerMI::StopNotActive(AliTPCseed * seed, Int_t row0, Float_t th0, Float_t th1,
2496 // stop not active tracks
2497 // take th1 as threshold for number of founded to number of foundable on last 10 active rows
2498 // take th2 as threshold for number of founded to number of foundable on last 20 active rows
2501 Int_t foundable = 0;
2502 Int_t maxindex = seed->GetLastPoint(); //last foundable row
2503 if (seed->GetNFoundable()*th0 > seed->GetNumberOfClusters()) {
2504 seed->Desactivate(10) ;
2508 for (Int_t i=row0; i<maxindex; i++){
2509 Int_t index = seed->GetClusterIndex2(i);
2510 if (index!=-1) foundable++;
2512 if (foundable<=30) sumgood1++;
2513 if (foundable<=50) {
2520 if (foundable>=30.){
2521 if (sumgood1<(th1*30.)) seed->Desactivate(10);
2524 if (sumgood2<(th2*50.)) seed->Desactivate(10);
2528 Int_t AliTPCtrackerMI::RefitInward(AliESDEvent *event)
2531 // back propagation of ESD tracks
2534 const Int_t kMaxFriendTracks=2000;
2538 //PrepareForProlongation(fSeeds,1);
2539 PropagateForward2(fSeeds);
2540 RemoveUsed2(fSeeds,0.4,0.4,20);
2542 TObjArray arraySeed(fSeeds->GetEntries());
2543 for (Int_t i=0;i<fSeeds->GetEntries();i++) {
2544 arraySeed.AddAt(fSeeds->At(i),i);
2546 SignShared(&arraySeed);
2547 // FindCurling(fSeeds, event,2); // find multi found tracks
2548 FindSplitted(fSeeds, event,2); // find multi found tracks
2551 Int_t nseed = fSeeds->GetEntriesFast();
2552 for (Int_t i=0;i<nseed;i++){
2553 AliTPCseed * seed = (AliTPCseed*) fSeeds->UncheckedAt(i);
2554 if (!seed) continue;
2555 if (seed->GetKinkIndex(0)>0) UpdateKinkQualityD(seed); // update quality informations for kinks
2557 seed->PropagateTo(fParam->GetInnerRadiusLow());
2558 seed->UpdatePoints();
2560 AliESDtrack *esd=event->GetTrack(i);
2561 seed->CookdEdx(0.02,0.6);
2562 CookLabel(seed,0.1); //For comparison only
2563 esd->SetTPCClusterMap(seed->GetClusterMap());
2564 esd->SetTPCSharedMap(seed->GetSharedMap());
2566 if (AliTPCReconstructor::StreamLevel()>1 && seed!=0&&esd!=0) {
2567 TTreeSRedirector &cstream = *fDebugStreamer;
2574 if (seed->GetNumberOfClusters()>15){
2575 esd->UpdateTrackParams(seed,AliESDtrack::kTPCrefit);
2576 esd->SetTPCPoints(seed->GetPoints());
2577 esd->SetTPCPointsF(seed->GetNFoundable());
2578 Int_t ndedx = seed->GetNCDEDX(0)+seed->GetNCDEDX(1)+seed->GetNCDEDX(2)+seed->GetNCDEDX(3);
2579 Float_t sdedx = (seed->GetSDEDX(0)+seed->GetSDEDX(1)+seed->GetSDEDX(2)+seed->GetSDEDX(3))*0.25;
2580 Float_t dedx = seed->GetdEdx();
2581 esd->SetTPCsignal(dedx, sdedx, ndedx);
2583 // add seed to the esd track in Calib level
2585 Bool_t storeFriend = gRandom->Rndm()<(kMaxFriendTracks)/Float_t(nseed);
2586 if (AliTPCReconstructor::StreamLevel()>0 &&storeFriend){
2587 AliTPCseed * seedCopy = new AliTPCseed(*seed, kTRUE);
2588 esd->AddCalibObject(seedCopy);
2593 //printf("problem\n");
2596 //FindKinks(fSeeds,event);
2597 Info("RefitInward","Number of refitted tracks %d",ntracks);
2602 Int_t AliTPCtrackerMI::PropagateBack(AliESDEvent *event)
2605 // back propagation of ESD tracks
2611 PropagateBack(fSeeds);
2612 RemoveUsed2(fSeeds,0.4,0.4,20);
2613 //FindCurling(fSeeds, fEvent,1);
2614 FindSplitted(fSeeds, event,1); // find multi found tracks
2617 Int_t nseed = fSeeds->GetEntriesFast();
2619 for (Int_t i=0;i<nseed;i++){
2620 AliTPCseed * seed = (AliTPCseed*) fSeeds->UncheckedAt(i);
2621 if (!seed) continue;
2622 if (seed->GetKinkIndex(0)<0) UpdateKinkQualityM(seed); // update quality informations for kinks
2623 seed->UpdatePoints();
2624 AliESDtrack *esd=event->GetTrack(i);
2625 seed->CookdEdx(0.02,0.6);
2626 CookLabel(seed,0.1); //For comparison only
2627 if (seed->GetNumberOfClusters()>15){
2628 esd->UpdateTrackParams(seed,AliESDtrack::kTPCout);
2629 esd->SetTPCPoints(seed->GetPoints());
2630 esd->SetTPCPointsF(seed->GetNFoundable());
2631 Int_t ndedx = seed->GetNCDEDX(0)+seed->GetNCDEDX(1)+seed->GetNCDEDX(2)+seed->GetNCDEDX(3);
2632 Float_t sdedx = (seed->GetSDEDX(0)+seed->GetSDEDX(1)+seed->GetSDEDX(2)+seed->GetSDEDX(3))*0.25;
2633 Float_t dedx = seed->GetdEdx();
2634 esd->SetTPCsignal(dedx, sdedx, ndedx);
2636 Int_t eventnumber = event->GetEventNumberInFile();// patch 28 fev 06
2637 // This is most likely NOT the event number you'd like to use. It has nothing to do with the 'real' event number
2638 if (AliTPCReconstructor::StreamLevel()>1) {
2639 (*fDebugStreamer)<<"Cback"<<
2641 "EventNrInFile="<<eventnumber<<
2642 "\n"; // patch 28 fev 06
2646 //FindKinks(fSeeds,event);
2647 Info("PropagateBack","Number of back propagated tracks %d",ntracks);
2654 void AliTPCtrackerMI::DeleteSeeds()
2659 Int_t nseed = fSeeds->GetEntriesFast();
2660 for (Int_t i=0;i<nseed;i++){
2661 AliTPCseed * seed = (AliTPCseed*)fSeeds->At(i);
2662 if (seed) delete fSeeds->RemoveAt(i);
2669 void AliTPCtrackerMI::ReadSeeds(AliESDEvent *event, Int_t direction)
2672 //read seeds from the event
2674 Int_t nentr=event->GetNumberOfTracks();
2676 Info("ReadSeeds", "Number of ESD tracks: %d\n", nentr);
2681 fSeeds = new TObjArray(nentr);
2685 for (Int_t i=0; i<nentr; i++) {
2686 AliESDtrack *esd=event->GetTrack(i);
2687 ULong_t status=esd->GetStatus();
2688 if (!(status&AliESDtrack::kTPCin)) continue;
2689 AliTPCtrack t(*esd);
2690 t.SetNumberOfClusters(0);
2691 // AliTPCseed *seed = new AliTPCseed(t,t.GetAlpha());
2692 AliTPCseed *seed = new AliTPCseed(t/*,t.GetAlpha()*/);
2693 seed->SetUniqueID(esd->GetID());
2694 for (Int_t ikink=0;ikink<3;ikink++) {
2695 Int_t index = esd->GetKinkIndex(ikink);
2696 seed->GetKinkIndexes()[ikink] = index;
2697 if (index==0) continue;
2698 index = TMath::Abs(index);
2699 AliESDkink * kink = fEvent->GetKink(index-1);
2700 if (kink&&esd->GetKinkIndex(ikink)<0){
2701 if ((status & AliESDtrack::kTRDrefit) != 0) kink->SetStatus(1,2);
2702 if ((status & AliESDtrack::kITSout) != 0) kink->SetStatus(1,0);
2704 if (kink&&esd->GetKinkIndex(ikink)>0){
2705 if ((status & AliESDtrack::kTRDrefit) != 0) kink->SetStatus(1,6);
2706 if ((status & AliESDtrack::kITSout) != 0) kink->SetStatus(1,4);
2710 if (((status&AliESDtrack::kITSout)==0)&&(direction==1)) seed->ResetCovariance(10.);
2711 if ( direction ==2 &&(status & AliESDtrack::kTRDrefit) == 0 ) seed->ResetCovariance(10.);
2712 if ( direction ==2 && ((status & AliESDtrack::kTPCout) == 0) ) {
2717 if ( direction ==2 &&(status & AliESDtrack::kTRDrefit) > 0 ) {
2718 Double_t par0[5],par1[5],alpha,x;
2719 esd->GetInnerExternalParameters(alpha,x,par0);
2720 esd->GetExternalParameters(x,par1);
2721 Double_t delta1 = TMath::Abs(par0[4]-par1[4])/(0.000000001+TMath::Abs(par0[4]+par1[4]));
2722 Double_t delta2 = TMath::Abs(par0[3]-par1[3]);
2724 if (esd->GetTRDncls()>0) trdchi2 = esd->GetTRDchi2()/esd->GetTRDncls();
2725 //reset covariance if suspicious
2726 if ( (delta1>0.1) || (delta2>0.006) ||trdchi2>7.)
2727 seed->ResetCovariance(10.);
2732 // rotate to the local coordinate system
2734 fSectors=fInnerSec; fN=fkNIS;
2735 Double_t alpha=seed->GetAlpha() - fSectors->GetAlphaShift();
2736 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
2737 if (alpha < 0. ) alpha += 2.*TMath::Pi();
2738 Int_t ns=Int_t(alpha/fSectors->GetAlpha())%fN;
2739 alpha =ns*fSectors->GetAlpha() + fSectors->GetAlphaShift();
2740 if (alpha<-TMath::Pi()) alpha += 2*TMath::Pi();
2741 if (alpha>=TMath::Pi()) alpha -= 2*TMath::Pi();
2742 alpha-=seed->GetAlpha();
2743 if (!seed->Rotate(alpha)) {
2749 if (esd->GetKinkIndex(0)<=0){
2750 for (Int_t irow=0;irow<160;irow++){
2751 Int_t index = seed->GetClusterIndex2(irow);
2754 AliTPCclusterMI * cl = GetClusterMI(index);
2755 seed->SetClusterPointer(irow,cl);
2757 if ((index & 0x8000)==0){
2758 cl->Use(10); // accepted cluster
2760 cl->Use(6); // close cluster not accepted
2763 Info("ReadSeeds","Not found cluster");
2768 fSeeds->AddAt(seed,i);
2774 //_____________________________________________________________________________
2775 void AliTPCtrackerMI::MakeSeeds3(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, Float_t cuts[4],
2776 Float_t deltay, Int_t ddsec) {
2777 //-----------------------------------------------------------------
2778 // This function creates track seeds.
2779 // SEEDING WITH VERTEX CONSTRAIN
2780 //-----------------------------------------------------------------
2781 // cuts[0] - fP4 cut
2782 // cuts[1] - tan(phi) cut
2783 // cuts[2] - zvertex cut
2784 // cuts[3] - fP3 cut
2792 Double_t x[5], c[15];
2793 // Int_t di = i1-i2;
2795 AliTPCseed * seed = new AliTPCseed();
2796 Double_t alpha=fSectors->GetAlpha(), shift=fSectors->GetAlphaShift();
2797 Double_t cs=cos(alpha), sn=sin(alpha);
2799 // Double_t x1 =fOuterSec->GetX(i1);
2800 //Double_t xx2=fOuterSec->GetX(i2);
2802 Double_t x1 =GetXrow(i1);
2803 Double_t xx2=GetXrow(i2);
2805 Double_t x3=GetX(), y3=GetY(), z3=GetZ();
2807 Int_t imiddle = (i2+i1)/2; //middle pad row index
2808 Double_t xm = GetXrow(imiddle); // radius of middle pad-row
2809 const AliTPCtrackerRow& krm=GetRow(sec,imiddle); //middle pad -row
2813 const AliTPCtrackerRow& kr1=GetRow(ns,i1);
2814 Double_t ymax = GetMaxY(i1)-kr1.GetDeadZone()-1.5;
2815 Double_t ymaxm = GetMaxY(imiddle)-kr1.GetDeadZone()-1.5;
2818 // change cut on curvature if it can't reach this layer
2819 // maximal curvature set to reach it
2820 Double_t dvertexmax = TMath::Sqrt((x1-x3)*(x1-x3)+(ymax+5-y3)*(ymax+5-y3));
2821 if (dvertexmax*0.5*cuts[0]>0.85){
2822 cuts[0] = 0.85/(dvertexmax*0.5+1.);
2824 Double_t r2min = 1/(cuts[0]*cuts[0]); //minimal square of radius given by cut
2827 if (deltay>0) ddsec = 0;
2828 // loop over clusters
2829 for (Int_t is=0; is < kr1; is++) {
2831 if (kr1[is]->IsUsed(10)) continue;
2832 Double_t y1=kr1[is]->GetY(), z1=kr1[is]->GetZ();
2833 //if (TMath::Abs(y1)>ymax) continue;
2835 if (deltay>0 && TMath::Abs(ymax-TMath::Abs(y1))> deltay ) continue; // seed only at the edge
2837 // find possible directions
2838 Float_t anglez = (z1-z3)/(x1-x3);
2839 Float_t extraz = z1 - anglez*(x1-xx2); // extrapolated z
2842 //find rotation angles relative to line given by vertex and point 1
2843 Double_t dvertex2 = (x1-x3)*(x1-x3)+(y1-y3)*(y1-y3);
2844 Double_t dvertex = TMath::Sqrt(dvertex2);
2845 Double_t angle13 = TMath::ATan((y1-y3)/(x1-x3));
2846 Double_t cs13 = cos(-angle13), sn13 = sin(-angle13);
2849 // loop over 2 sectors
2855 Double_t dddz1=0; // direction of delta inclination in z axis
2862 for (Int_t dsec = dsec1; dsec<=dsec2;dsec++){
2863 Int_t sec2 = sec + dsec;
2865 // AliTPCtrackerRow& kr2 = fOuterSec[(sec2+fkNOS)%fkNOS][i2];
2866 //AliTPCtrackerRow& kr2m = fOuterSec[(sec2+fkNOS)%fkNOS][imiddle];
2867 AliTPCtrackerRow& kr2 = GetRow((sec2+fkNOS)%fkNOS,i2);
2868 AliTPCtrackerRow& kr2m = GetRow((sec2+fkNOS)%fkNOS,imiddle);
2869 Int_t index1 = TMath::Max(kr2.Find(extraz-0.6-dddz1*TMath::Abs(z1)*0.05)-1,0);
2870 Int_t index2 = TMath::Min(kr2.Find(extraz+0.6+dddz2*TMath::Abs(z1)*0.05)+1,kr2);
2872 // rotation angles to p1-p3
2873 Double_t cs13r = cos(-angle13+dsec*alpha)/dvertex, sn13r = sin(-angle13+dsec*alpha)/dvertex;
2874 Double_t x2, y2, z2;
2876 // Double_t dymax = maxangle*TMath::Abs(x1-xx2);
2879 Double_t dxx0 = (xx2-x3)*cs13r;
2880 Double_t dyy0 = (xx2-x3)*sn13r;
2881 for (Int_t js=index1; js < index2; js++) {
2882 const AliTPCclusterMI *kcl = kr2[js];
2883 if (kcl->IsUsed(10)) continue;
2885 //calcutate parameters
2887 Double_t yy0 = dyy0 +(kcl->GetY()-y3)*cs13r;
2889 if (TMath::Abs(yy0)<0.000001) continue;
2890 Double_t xx0 = dxx0 -(kcl->GetY()-y3)*sn13r;
2891 Double_t y0 = 0.5*(xx0*xx0+yy0*yy0-xx0)/yy0;
2892 Double_t r02 = (0.25+y0*y0)*dvertex2;
2893 //curvature (radius) cut
2894 if (r02<r2min) continue;
2898 Double_t c0 = 1/TMath::Sqrt(r02);
2902 //Double_t dfi0 = 2.*TMath::ASin(dvertex*c0*0.5);
2903 //Double_t dfi1 = 2.*TMath::ASin(TMath::Sqrt(yy0*yy0+(1-xx0)*(1-xx0))*dvertex*c0*0.5);
2904 Double_t dfi0 = 2.*AliTPCFastMath::FastAsin(dvertex*c0*0.5);
2905 Double_t dfi1 = 2.*AliTPCFastMath::FastAsin(TMath::Sqrt(yy0*yy0+(1-xx0)*(1-xx0))*dvertex*c0*0.5);
2908 Double_t z0 = kcl->GetZ();
2909 Double_t zzzz2 = z1-(z1-z3)*dfi1/dfi0;
2910 if (TMath::Abs(zzzz2-z0)>0.5) continue;
2913 Double_t dip = (z1-z0)*c0/dfi1;
2914 Double_t x0 = (0.5*cs13+y0*sn13)*dvertex*c0;
2925 x2= xx2*cs-y2*sn*dsec;
2926 y2=+xx2*sn*dsec+y2*cs;
2936 // do we have cluster at the middle ?
2938 GetProlongation(x1,xm,x,ym,zm);
2940 AliTPCclusterMI * cm=0;
2941 if (TMath::Abs(ym)-ymaxm<0){
2942 cm = krm.FindNearest2(ym,zm,1.0,0.6,dummy);
2943 if ((!cm) || (cm->IsUsed(10))) {
2948 // rotate y1 to system 0
2949 // get state vector in rotated system
2950 Double_t yr1 = (-0.5*sn13+y0*cs13)*dvertex*c0;
2951 Double_t xr2 = x0*cs+yr1*sn*dsec;
2952 Double_t xr[5]={kcl->GetY(),kcl->GetZ(), xr2, dip, c0};
2954 GetProlongation(xx2,xm,xr,ym,zm);
2955 if (TMath::Abs(ym)-ymaxm<0){
2956 cm = kr2m.FindNearest2(ym,zm,1.0,0.6,dummy);
2957 if ((!cm) || (cm->IsUsed(10))) {
2967 dym = ym - cm->GetY();
2968 dzm = zm - cm->GetZ();
2975 Double_t sy1=kr1[is]->GetSigmaY2()*2., sz1=kr1[is]->GetSigmaZ2()*2.;
2976 Double_t sy2=kcl->GetSigmaY2()*2., sz2=kcl->GetSigmaZ2()*2.;
2977 //Double_t sy3=400*3./12., sy=0.1, sz=0.1;
2978 Double_t sy3=25000*x[4]*x[4]+0.1, sy=0.1, sz=0.1;
2979 //Double_t sy3=25000*x[4]*x[4]*60+0.5, sy=0.1, sz=0.1;
2981 Double_t f40=(F1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
2982 Double_t f42=(F1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
2983 Double_t f43=(F1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
2984 Double_t f20=(F2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
2985 Double_t f22=(F2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
2986 Double_t f23=(F2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
2988 Double_t f30=(F3(x1,y1+sy,x2,y2,z1,z2)-x[3])/sy;
2989 Double_t f31=(F3(x1,y1,x2,y2,z1+sz,z2)-x[3])/sz;
2990 Double_t f32=(F3(x1,y1,x2,y2+sy,z1,z2)-x[3])/sy;
2991 Double_t f34=(F3(x1,y1,x2,y2,z1,z2+sz)-x[3])/sz;
2995 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
2996 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
2997 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
2998 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
2999 c[13]=f30*sy1*f40+f32*sy2*f42;
3000 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3002 // if (!BuildSeed(kr1[is],kcl,0,x1,x2,x3,x,c)) continue;
3004 UInt_t index=kr1.GetIndex(is);
3005 seed->~AliTPCseed(); // this does not set the pointer to 0...
3006 AliTPCseed *track=new(seed) AliTPCseed(x1, ns*alpha+shift, x, c, index);
3008 track->SetIsSeeding(kTRUE);
3009 track->SetSeed1(i1);
3010 track->SetSeed2(i2);
3011 track->SetSeedType(3);
3015 FollowProlongation(*track, (i1+i2)/2,1);
3016 Int_t foundable,found,shared;
3017 track->GetClusterStatistic((i1+i2)/2,i1, found, foundable, shared, kTRUE);
3018 if ((found<0.55*foundable) || shared>0.5*found || (track->GetSigmaY2()+track->GetSigmaZ2())>0.5){
3020 seed->~AliTPCseed();
3026 FollowProlongation(*track, i2,1);
3030 track->SetBConstrain(1);
3031 // track->fLastPoint = i1+fInnerSec->GetNRows(); // first cluster in track position
3032 track->SetLastPoint(i1); // first cluster in track position
3033 track->SetFirstPoint(track->GetLastPoint());
3035 if (track->GetNumberOfClusters()<(i1-i2)*0.5 ||
3036 track->GetNumberOfClusters() < track->GetNFoundable()*0.6 ||
3037 track->GetNShared()>0.4*track->GetNumberOfClusters() ) {
3039 seed->~AliTPCseed();
3043 // Z VERTEX CONDITION
3044 Double_t zv, bz=GetBz();
3045 if ( !track->GetZAt(0.,bz,zv) ) continue;
3046 if (TMath::Abs(zv-z3)>cuts[2]) {
3047 FollowProlongation(*track, TMath::Max(i2-20,0));
3048 if ( !track->GetZAt(0.,bz,zv) ) continue;
3049 if (TMath::Abs(zv-z3)>cuts[2]){
3050 FollowProlongation(*track, TMath::Max(i2-40,0));
3051 if ( !track->GetZAt(0.,bz,zv) ) continue;
3052 if (TMath::Abs(zv-z3)>cuts[2] &&(track->GetNumberOfClusters() > track->GetNFoundable()*0.7)){
3053 // make seed without constrain
3054 AliTPCseed * track2 = MakeSeed(track,0.2,0.5,1.);
3055 FollowProlongation(*track2, i2,1);
3056 track2->SetBConstrain(kFALSE);
3057 track2->SetSeedType(1);
3058 arr->AddLast(track2);
3060 seed->~AliTPCseed();
3065 seed->~AliTPCseed();
3072 track->SetSeedType(0);
3073 arr->AddLast(track);
3074 seed = new AliTPCseed;
3076 // don't consider other combinations
3077 if (track->GetNumberOfClusters() > track->GetNFoundable()*0.8)
3083 Info("MakeSeeds3","\nSeeding statistic:\t%d\t%d\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin,nout1,nout2);
3089 void AliTPCtrackerMI::MakeSeeds5(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, Float_t cuts[4],
3094 //-----------------------------------------------------------------
3095 // This function creates track seeds.
3096 //-----------------------------------------------------------------
3097 // cuts[0] - fP4 cut
3098 // cuts[1] - tan(phi) cut
3099 // cuts[2] - zvertex cut
3100 // cuts[3] - fP3 cut
3110 Double_t x[5], c[15];
3112 // make temporary seed
3113 AliTPCseed * seed = new AliTPCseed;
3114 Double_t alpha=fOuterSec->GetAlpha(), shift=fOuterSec->GetAlphaShift();
3115 // Double_t cs=cos(alpha), sn=sin(alpha);
3120 Double_t x1 = GetXrow(i1-1);
3121 const AliTPCtrackerRow& kr1=GetRow(sec,i1-1);
3122 Double_t y1max = GetMaxY(i1-1)-kr1.GetDeadZone()-1.5;
3124 Double_t x1p = GetXrow(i1);
3125 const AliTPCtrackerRow& kr1p=GetRow(sec,i1);
3127 Double_t x1m = GetXrow(i1-2);
3128 const AliTPCtrackerRow& kr1m=GetRow(sec,i1-2);
3131 //last 3 padrow for seeding
3132 AliTPCtrackerRow& kr3 = GetRow((sec+fkNOS)%fkNOS,i1-7);
3133 Double_t x3 = GetXrow(i1-7);
3134 // Double_t y3max= GetMaxY(i1-7)-kr3.fDeadZone-1.5;
3136 AliTPCtrackerRow& kr3p = GetRow((sec+fkNOS)%fkNOS,i1-6);
3137 Double_t x3p = GetXrow(i1-6);
3139 AliTPCtrackerRow& kr3m = GetRow((sec+fkNOS)%fkNOS,i1-8);
3140 Double_t x3m = GetXrow(i1-8);
3145 Int_t im = i1-4; //middle pad row index
3146 Double_t xm = GetXrow(im); // radius of middle pad-row
3147 const AliTPCtrackerRow& krm=GetRow(sec,im); //middle pad -row
3148 // Double_t ymmax = GetMaxY(im)-kr1.fDeadZone-1.5;
3151 Double_t deltax = x1-x3;
3152 Double_t dymax = deltax*cuts[1];
3153 Double_t dzmax = deltax*cuts[3];
3155 // loop over clusters
3156 for (Int_t is=0; is < kr1; is++) {
3158 if (kr1[is]->IsUsed(10)) continue;
3159 Double_t y1=kr1[is]->GetY(), z1=kr1[is]->GetZ();
3161 if (deltay>0 && TMath::Abs(y1max-TMath::Abs(y1))> deltay ) continue; // seed only at the edge
3163 Int_t index1 = TMath::Max(kr3.Find(z1-dzmax)-1,0);
3164 Int_t index2 = TMath::Min(kr3.Find(z1+dzmax)+1,kr3);
3170 for (Int_t js=index1; js < index2; js++) {
3171 const AliTPCclusterMI *kcl = kr3[js];
3172 if (kcl->IsUsed(10)) continue;
3174 // apply angular cuts
3175 if (TMath::Abs(y1-y3)>dymax) continue;
3178 if (TMath::Abs(z1-z3)>dzmax) continue;
3180 Double_t angley = (y1-y3)/(x1-x3);
3181 Double_t anglez = (z1-z3)/(x1-x3);
3183 Double_t erry = TMath::Abs(angley)*(x1-x1m)*0.5+0.5;
3184 Double_t errz = TMath::Abs(anglez)*(x1-x1m)*0.5+0.5;
3186 Double_t yyym = angley*(xm-x1)+y1;
3187 Double_t zzzm = anglez*(xm-x1)+z1;
3189 const AliTPCclusterMI *kcm = krm.FindNearest2(yyym,zzzm,erry,errz,index);
3191 if (kcm->IsUsed(10)) continue;
3193 erry = TMath::Abs(angley)*(x1-x1m)*0.4+0.5;
3194 errz = TMath::Abs(anglez)*(x1-x1m)*0.4+0.5;
3201 // look around first
3202 const AliTPCclusterMI *kc1m = kr1m.FindNearest2(angley*(x1m-x1)+y1,
3208 if (kc1m->IsUsed(10)) used++;
3210 const AliTPCclusterMI *kc1p = kr1p.FindNearest2(angley*(x1p-x1)+y1,
3216 if (kc1p->IsUsed(10)) used++;
3218 if (used>1) continue;
3219 if (found<1) continue;
3223 const AliTPCclusterMI *kc3m = kr3m.FindNearest2(angley*(x3m-x3)+y3,
3229 if (kc3m->IsUsed(10)) used++;
3233 const AliTPCclusterMI *kc3p = kr3p.FindNearest2(angley*(x3p-x3)+y3,
3239 if (kc3p->IsUsed(10)) used++;
3243 if (used>1) continue;
3244 if (found<3) continue;
3254 x[4]=F1(x1,y1,x2,y2,x3,y3);
3255 //if (TMath::Abs(x[4]) >= cuts[0]) continue;
3258 x[2]=F2(x1,y1,x2,y2,x3,y3);
3261 x[3]=F3n(x1,y1,x2,y2,z1,z2,x[4]);
3262 //if (TMath::Abs(x[3]) > cuts[3]) continue;
3266 Double_t sy1=0.1, sz1=0.1;
3267 Double_t sy2=0.1, sz2=0.1;
3268 Double_t sy3=0.1, sy=0.1, sz=0.1;
3270 Double_t f40=(F1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
3271 Double_t f42=(F1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
3272 Double_t f43=(F1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
3273 Double_t f20=(F2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
3274 Double_t f22=(F2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
3275 Double_t f23=(F2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
3277 Double_t f30=(F3(x1,y1+sy,x2,y2,z1,z2)-x[3])/sy;
3278 Double_t f31=(F3(x1,y1,x2,y2,z1+sz,z2)-x[3])/sz;
3279 Double_t f32=(F3(x1,y1,x2,y2+sy,z1,z2)-x[3])/sy;
3280 Double_t f34=(F3(x1,y1,x2,y2,z1,z2+sz)-x[3])/sz;
3284 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3285 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3286 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3287 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3288 c[13]=f30*sy1*f40+f32*sy2*f42;
3289 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3291 // if (!BuildSeed(kr1[is],kcl,0,x1,x2,x3,x,c)) continue;
3293 UInt_t index=kr1.GetIndex(is);
3294 seed->~AliTPCseed();
3295 AliTPCseed *track=new(seed) AliTPCseed(x1, sec*alpha+shift, x, c, index);
3297 track->SetIsSeeding(kTRUE);
3300 FollowProlongation(*track, i1-7,1);
3301 if (track->GetNumberOfClusters() < track->GetNFoundable()*0.75 ||
3302 track->GetNShared()>0.6*track->GetNumberOfClusters() || ( track->GetSigmaY2()+ track->GetSigmaZ2())>0.6){
3304 seed->~AliTPCseed();
3310 FollowProlongation(*track, i2,1);
3311 track->SetBConstrain(0);
3312 track->SetLastPoint(i1+fInnerSec->GetNRows()); // first cluster in track position
3313 track->SetFirstPoint(track->GetLastPoint());
3315 if (track->GetNumberOfClusters()<(i1-i2)*0.5 ||
3316 track->GetNumberOfClusters()<track->GetNFoundable()*0.7 ||
3317 track->GetNShared()>2. || track->GetChi2()/track->GetNumberOfClusters()>6 || ( track->GetSigmaY2()+ track->GetSigmaZ2())>0.5 ) {
3319 seed->~AliTPCseed();
3324 FollowProlongation(*track, TMath::Max(i2-10,0),1);
3325 AliTPCseed * track2 = MakeSeed(track,0.2,0.5,0.9);
3326 FollowProlongation(*track2, i2,1);
3327 track2->SetBConstrain(kFALSE);
3328 track2->SetSeedType(4);
3329 arr->AddLast(track2);
3331 seed->~AliTPCseed();
3335 //arr->AddLast(track);
3336 //seed = new AliTPCseed;
3342 Info("MakeSeeds5","\nSeeding statiistic:\t%d\t%d\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin,nout1,nout2,nout3);
3348 //_____________________________________________________________________________
3349 void AliTPCtrackerMI::MakeSeeds2(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, Float_t */*cuts[4]*/,
3350 Float_t deltay, Bool_t /*bconstrain*/) {
3351 //-----------------------------------------------------------------
3352 // This function creates track seeds - without vertex constraint
3353 //-----------------------------------------------------------------
3354 // cuts[0] - fP4 cut - not applied
3355 // cuts[1] - tan(phi) cut
3356 // cuts[2] - zvertex cut - not applied
3357 // cuts[3] - fP3 cut
3367 Double_t alpha=fOuterSec->GetAlpha(), shift=fOuterSec->GetAlphaShift();
3368 // Double_t cs=cos(alpha), sn=sin(alpha);
3369 Int_t row0 = (i1+i2)/2;
3370 Int_t drow = (i1-i2)/2;
3371 const AliTPCtrackerRow& kr0=fSectors[sec][row0];
3372 AliTPCtrackerRow * kr=0;
3374 AliTPCpolyTrack polytrack;
3375 Int_t nclusters=fSectors[sec][row0];
3376 AliTPCseed * seed = new AliTPCseed;
3381 for (Int_t is=0; is < nclusters; is++) { //LOOP over clusters
3383 Int_t nfoundable =0;
3384 for (Int_t iter =1; iter<2; iter++){ //iterations
3385 const AliTPCtrackerRow& krm=fSectors[sec][row0-iter];
3386 const AliTPCtrackerRow& krp=fSectors[sec][row0+iter];
3387 const AliTPCclusterMI * cl= kr0[is];
3389 if (cl->IsUsed(10)) {
3395 Double_t x = kr0.GetX();
3396 // Initialization of the polytrack
3401 Double_t y0= cl->GetY();
3402 Double_t z0= cl->GetZ();
3406 Double_t ymax = fSectors->GetMaxY(row0)-kr0.GetDeadZone()-1.5;
3407 if (deltay>0 && TMath::Abs(ymax-TMath::Abs(y0))> deltay ) continue; // seed only at the edge
3409 erry = (0.5)*cl->GetSigmaY2()/TMath::Sqrt(cl->GetQ())*6;
3410 errz = (0.5)*cl->GetSigmaZ2()/TMath::Sqrt(cl->GetQ())*6;
3411 polytrack.AddPoint(x,y0,z0,erry, errz);
3414 if (cl->IsUsed(10)) sumused++;
3417 Float_t roady = (5*TMath::Sqrt(cl->GetSigmaY2()+0.2)+1.)*iter;
3418 Float_t roadz = (5*TMath::Sqrt(cl->GetSigmaZ2()+0.2)+1.)*iter;
3421 AliTPCclusterMI * cl1 = krm.FindNearest(y0,z0,roady,roadz);
3422 if (cl1 && TMath::Abs(ymax-TMath::Abs(y0))) {
3423 erry = (0.5)*cl1->GetSigmaY2()/TMath::Sqrt(cl1->GetQ())*3;
3424 errz = (0.5)*cl1->GetSigmaZ2()/TMath::Sqrt(cl1->GetQ())*3;
3425 if (cl1->IsUsed(10)) sumused++;
3426 polytrack.AddPoint(x,cl1->GetY(),cl1->GetZ(),erry,errz);
3430 AliTPCclusterMI * cl2 = krp.FindNearest(y0,z0,roady,roadz);
3432 erry = (0.5)*cl2->GetSigmaY2()/TMath::Sqrt(cl2->GetQ())*3;
3433 errz = (0.5)*cl2->GetSigmaZ2()/TMath::Sqrt(cl2->GetQ())*3;
3434 if (cl2->IsUsed(10)) sumused++;
3435 polytrack.AddPoint(x,cl2->GetY(),cl2->GetZ(),erry,errz);
3438 if (sumused>0) continue;
3440 polytrack.UpdateParameters();
3446 nfoundable = polytrack.GetN();
3447 nfound = nfoundable;
3449 for (Int_t ddrow = iter+1; ddrow<drow;ddrow++){
3450 Float_t maxdist = 0.8*(1.+3./(ddrow));
3451 for (Int_t delta = -1;delta<=1;delta+=2){
3452 Int_t row = row0+ddrow*delta;
3453 kr = &(fSectors[sec][row]);
3454 Double_t xn = kr->GetX();
3455 Double_t ymax = fSectors->GetMaxY(row)-kr->GetDeadZone()-1.5;
3456 polytrack.GetFitPoint(xn,yn,zn);
3457 if (TMath::Abs(yn)>ymax) continue;
3459 AliTPCclusterMI * cln = kr->FindNearest(yn,zn,roady,roadz);
3461 Float_t dist = TMath::Sqrt( (yn-cln->GetY())*(yn-cln->GetY())+(zn-cln->GetZ())*(zn-cln->GetZ()));
3464 erry = (dist+0.3)*cln->GetSigmaY2()/TMath::Sqrt(cln->GetQ())*(1.+1./(ddrow));
3465 errz = (dist+0.3)*cln->GetSigmaZ2()/TMath::Sqrt(cln->GetQ())*(1.+1./(ddrow));
3466 if (cln->IsUsed(10)) {
3467 // printf("used\n");
3475 polytrack.AddPoint(xn,cln->GetY(),cln->GetZ(),erry, errz);
3480 if ( (sumused>3) || (sumused>0.5*nfound) || (nfound<0.6*nfoundable)) break;
3481 polytrack.UpdateParameters();
3484 if ( (sumused>3) || (sumused>0.5*nfound)) {
3485 //printf("sumused %d\n",sumused);
3490 polytrack.GetFitDerivation(kr0.GetX(),dy,dz);
3491 AliTPCpolyTrack track2;
3493 polytrack.Refit(track2,0.5+TMath::Abs(dy)*0.3,0.4+TMath::Abs(dz)*0.3);
3494 if (track2.GetN()<0.5*nfoundable) continue;
3497 if ((nfound>0.6*nfoundable) &&( nfoundable>0.4*(i1-i2))) {
3499 // test seed with and without constrain
3500 for (Int_t constrain=0; constrain<=0;constrain++){
3501 // add polytrack candidate
3503 Double_t x[5], c[15];
3504 Double_t x1,x2,x3,y1,y2,y3,z1,z2,z3;
3505 track2.GetBoundaries(x3,x1);
3507 track2.GetFitPoint(x1,y1,z1);
3508 track2.GetFitPoint(x2,y2,z2);
3509 track2.GetFitPoint(x3,y3,z3);
3511 //is track pointing to the vertex ?
3514 polytrack.GetFitPoint(x0,y0,z0);
3527 x[4]=F1(x1,y1,x2,y2,x3,y3);
3529 // if (TMath::Abs(x[4]) >= cuts[0]) continue; //
3530 x[2]=F2(x1,y1,x2,y2,x3,y3);
3532 //if (TMath::Abs(x[4]*x1-x[2]) >= cuts[1]) continue;
3533 //x[3]=F3(x1,y1,x2,y2,z1,z2);
3534 x[3]=F3n(x1,y1,x3,y3,z1,z3,x[4]);
3535 //if (TMath::Abs(x[3]) > cuts[3]) continue;
3538 Double_t sy =0.1, sz =0.1;
3539 Double_t sy1=0.02, sz1=0.02;
3540 Double_t sy2=0.02, sz2=0.02;
3544 sy3=25000*x[4]*x[4]+0.1, sy=0.1, sz=0.1;
3547 Double_t f40=(F1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
3548 Double_t f42=(F1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
3549 Double_t f43=(F1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
3550 Double_t f20=(F2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
3551 Double_t f22=(F2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
3552 Double_t f23=(F2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
3554 Double_t f30=(F3(x1,y1+sy,x3,y3,z1,z3)-x[3])/sy;
3555 Double_t f31=(F3(x1,y1,x3,y3,z1+sz,z3)-x[3])/sz;
3556 Double_t f32=(F3(x1,y1,x3,y3+sy,z1,z3)-x[3])/sy;
3557 Double_t f34=(F3(x1,y1,x3,y3,z1,z3+sz)-x[3])/sz;
3562 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3563 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3564 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3565 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3566 c[13]=f30*sy1*f40+f32*sy2*f42;
3567 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3569 //Int_t row1 = fSectors->GetRowNumber(x1);
3570 Int_t row1 = GetRowNumber(x1);
3574 seed->~AliTPCseed();
3575 AliTPCseed *track=new(seed) AliTPCseed(x1,sec*alpha+shift,x,c,index);
3576 track->SetIsSeeding(kTRUE);
3577 Int_t rc=FollowProlongation(*track, i2);
3578 if (constrain) track->SetBConstrain(1);
3580 track->SetBConstrain(0);
3581 track->SetLastPoint(row1+fInnerSec->GetNRows()); // first cluster in track position
3582 track->SetFirstPoint(track->GetLastPoint());
3584 if (rc==0 || track->GetNumberOfClusters()<(i1-i2)*0.5 ||
3585 track->GetNumberOfClusters() < track->GetNFoundable()*0.6 ||
3586 track->GetNShared()>0.4*track->GetNumberOfClusters()) {
3589 seed->~AliTPCseed();
3592 arr->AddLast(track);
3593 seed = new AliTPCseed;
3597 } // if accepted seed
3600 Info("MakeSeeds2","\nSeeding statiistic:\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin3);
3606 AliTPCseed *AliTPCtrackerMI::MakeSeed(AliTPCseed *track, Float_t r0, Float_t r1, Float_t r2)
3610 //reseed using track points
3611 Int_t p0 = int(r0*track->GetNumberOfClusters()); // point 0
3612 Int_t p1 = int(r1*track->GetNumberOfClusters());
3613 Int_t p2 = int(r2*track->GetNumberOfClusters()); // last point
3615 Double_t x0[3],x1[3],x2[3];
3616 for (Int_t i=0;i<3;i++){
3622 // find track position at given ratio of the length
3623 Int_t sec0=0, sec1=0, sec2=0;
3626 for (Int_t i=0;i<160;i++){
3627 if (track->GetClusterPointer(i)){
3629 AliTPCTrackerPoint *trpoint =track->GetTrackPoint(i);
3630 if ( (index<p0) || x0[0]<0 ){
3631 if (trpoint->GetX()>1){
3632 clindex = track->GetClusterIndex2(i);
3634 x0[0] = trpoint->GetX();
3635 x0[1] = trpoint->GetY();
3636 x0[2] = trpoint->GetZ();
3637 sec0 = ((clindex&0xff000000)>>24)%18;
3642 if ( (index<p1) &&(trpoint->GetX()>1)){
3643 clindex = track->GetClusterIndex2(i);
3645 x1[0] = trpoint->GetX();
3646 x1[1] = trpoint->GetY();
3647 x1[2] = trpoint->GetZ();
3648 sec1 = ((clindex&0xff000000)>>24)%18;
3651 if ( (index<p2) &&(trpoint->GetX()>1)){
3652 clindex = track->GetClusterIndex2(i);
3654 x2[0] = trpoint->GetX();
3655 x2[1] = trpoint->GetY();
3656 x2[2] = trpoint->GetZ();
3657 sec2 = ((clindex&0xff000000)>>24)%18;
3664 Double_t alpha, cs,sn, xx2,yy2;
3666 alpha = (sec1-sec2)*fSectors->GetAlpha();
3667 cs = TMath::Cos(alpha);
3668 sn = TMath::Sin(alpha);
3669 xx2= x1[0]*cs-x1[1]*sn;
3670 yy2= x1[0]*sn+x1[1]*cs;
3674 alpha = (sec0-sec2)*fSectors->GetAlpha();
3675 cs = TMath::Cos(alpha);
3676 sn = TMath::Sin(alpha);
3677 xx2= x0[0]*cs-x0[1]*sn;
3678 yy2= x0[0]*sn+x0[1]*cs;
3684 Double_t x[5],c[15];
3688 x[4]=F1(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]);
3689 // if (x[4]>1) return 0;
3690 x[2]=F2(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]);
3691 x[3]=F3n(x2[0],x2[1],x0[0],x0[1],x2[2],x0[2],x[4]);
3692 //if (TMath::Abs(x[3]) > 2.2) return 0;
3693 //if (TMath::Abs(x[2]) > 1.99) return 0;
3695 Double_t sy =0.1, sz =0.1;
3697 Double_t sy1=0.02+track->GetSigmaY2(), sz1=0.02+track->GetSigmaZ2();
3698 Double_t sy2=0.01+track->GetSigmaY2(), sz2=0.01+track->GetSigmaZ2();
3699 Double_t sy3=0.01+track->GetSigmaY2();
3701 Double_t f40=(F1(x2[0],x2[1]+sy,x1[0],x1[1],x0[0],x0[1])-x[4])/sy;
3702 Double_t f42=(F1(x2[0],x2[1],x1[0],x1[1]+sy,x0[0],x0[1])-x[4])/sy;
3703 Double_t f43=(F1(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]+sy)-x[4])/sy;
3704 Double_t f20=(F2(x2[0],x2[1]+sy,x1[0],x1[1],x0[0],x0[1])-x[2])/sy;
3705 Double_t f22=(F2(x2[0],x2[1],x1[0],x1[1]+sy,x0[0],x0[1])-x[2])/sy;
3706 Double_t f23=(F2(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]+sy)-x[2])/sy;
3708 Double_t f30=(F3(x2[0],x2[1]+sy,x0[0],x0[1],x2[2],x0[2])-x[3])/sy;
3709 Double_t f31=(F3(x2[0],x2[1],x0[0],x0[1],x2[2]+sz,x0[2])-x[3])/sz;
3710 Double_t f32=(F3(x2[0],x2[1],x0[0],x0[1]+sy,x2[2],x0[2])-x[3])/sy;
3711 Double_t f34=(F3(x2[0],x2[1],x0[0],x0[1],x2[2],x0[2]+sz)-x[3])/sz;
3716 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3717 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3718 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3719 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3720 c[13]=f30*sy1*f40+f32*sy2*f42;
3721 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3723 // Int_t row1 = fSectors->GetRowNumber(x2[0]);
3724 AliTPCseed *seed=new AliTPCseed(x2[0], sec2*fSectors->GetAlpha()+fSectors->GetAlphaShift(), x, c, 0);
3725 // Double_t y0,z0,y1,z1, y2,z2;
3726 //seed->GetProlongation(x0[0],y0,z0);
3727 // seed->GetProlongation(x1[0],y1,z1);
3728 //seed->GetProlongation(x2[0],y2,z2);
3730 seed->SetLastPoint(pp2);
3731 seed->SetFirstPoint(pp2);
3738 AliTPCseed *AliTPCtrackerMI::ReSeed(AliTPCseed *track, Float_t r0, Float_t r1, Float_t r2)
3742 //reseed using founded clusters
3744 // Find the number of clusters
3745 Int_t nclusters = 0;
3746 for (Int_t irow=0;irow<160;irow++){
3747 if (track->GetClusterIndex(irow)>0) nclusters++;
3751 ipos[0] = TMath::Max(int(r0*nclusters),0); // point 0 cluster
3752 ipos[1] = TMath::Min(int(r1*nclusters),nclusters-1); //
3753 ipos[2] = TMath::Min(int(r2*nclusters),nclusters-1); // last point
3757 Int_t row[3],sec[3]={0,0,0};
3759 // find track row position at given ratio of the length
3761 for (Int_t irow=0;irow<160;irow++){
3762 if (track->GetClusterIndex2(irow)<0) continue;
3764 for (Int_t ipoint=0;ipoint<3;ipoint++){
3765 if (index<=ipos[ipoint]) row[ipoint] = irow;
3769 //Get cluster and sector position
3770 for (Int_t ipoint=0;ipoint<3;ipoint++){
3771 Int_t clindex = track->GetClusterIndex2(row[ipoint]);
3772 AliTPCclusterMI * cl = GetClusterMI(clindex);
3775 // AliTPCclusterMI * cl = GetClusterMI(clindex);
3778 sec[ipoint] = ((clindex&0xff000000)>>24)%18;
3779 xyz[ipoint][0] = GetXrow(row[ipoint]);
3780 xyz[ipoint][1] = cl->GetY();
3781 xyz[ipoint][2] = cl->GetZ();
3785 // Calculate seed state vector and covariance matrix
3787 Double_t alpha, cs,sn, xx2,yy2;
3789 alpha = (sec[1]-sec[2])*fSectors->GetAlpha();
3790 cs = TMath::Cos(alpha);
3791 sn = TMath::Sin(alpha);
3792 xx2= xyz[1][0]*cs-xyz[1][1]*sn;
3793 yy2= xyz[1][0]*sn+xyz[1][1]*cs;
3797 alpha = (sec[0]-sec[2])*fSectors->GetAlpha();
3798 cs = TMath::Cos(alpha);
3799 sn = TMath::Sin(alpha);
3800 xx2= xyz[0][0]*cs-xyz[0][1]*sn;
3801 yy2= xyz[0][0]*sn+xyz[0][1]*cs;
3807 Double_t x[5],c[15];
3811 x[4]=F1(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
3812 x[2]=F2(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
3813 x[3]=F3n(xyz[2][0],xyz[2][1],xyz[0][0],xyz[0][1],xyz[2][2],xyz[0][2],x[4]);
3815 Double_t sy =0.1, sz =0.1;
3817 Double_t sy1=0.2, sz1=0.2;
3818 Double_t sy2=0.2, sz2=0.2;
3821 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;
3822 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;
3823 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;
3824 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;
3825 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;
3826 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;
3828 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;
3829 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;
3830 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;
3831 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;
3836 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3837 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3838 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3839 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3840 c[13]=f30*sy1*f40+f32*sy2*f42;
3841 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3843 // Int_t row1 = fSectors->GetRowNumber(xyz[2][0]);
3844 AliTPCseed *seed=new AliTPCseed(xyz[2][0], sec[2]*fSectors->GetAlpha()+fSectors->GetAlphaShift(), x, c, 0);
3845 seed->SetLastPoint(row[2]);
3846 seed->SetFirstPoint(row[2]);
3851 AliTPCseed *AliTPCtrackerMI::ReSeed(AliTPCseed *track,Int_t r0, Bool_t forward)
3855 //reseed using founded clusters
3858 Int_t row[3]={0,0,0};
3859 Int_t sec[3]={0,0,0};
3861 // forward direction
3863 for (Int_t irow=r0;irow<160;irow++){
3864 if (track->GetClusterIndex(irow)>0){
3869 for (Int_t irow=160;irow>r0;irow--){
3870 if (track->GetClusterIndex(irow)>0){
3875 for (Int_t irow=row[2]-15;irow>row[0];irow--){
3876 if (track->GetClusterIndex(irow)>0){
3884 for (Int_t irow=0;irow<r0;irow++){
3885 if (track->GetClusterIndex(irow)>0){
3890 for (Int_t irow=r0;irow>0;irow--){
3891 if (track->GetClusterIndex(irow)>0){
3896 for (Int_t irow=row[2]-15;irow>row[0];irow--){
3897 if (track->GetClusterIndex(irow)>0){
3904 if ((row[2]-row[0])<20) return 0;
3905 if (row[1]==0) return 0;
3908 //Get cluster and sector position
3909 for (Int_t ipoint=0;ipoint<3;ipoint++){
3910 Int_t clindex = track->GetClusterIndex2(row[ipoint]);
3911 AliTPCclusterMI * cl = GetClusterMI(clindex);
3914 // AliTPCclusterMI * cl = GetClusterMI(clindex);
3917 sec[ipoint] = ((clindex&0xff000000)>>24)%18;
3918 xyz[ipoint][0] = GetXrow(row[ipoint]);
3919 AliTPCTrackerPoint * point = track->GetTrackPoint(row[ipoint]);
3920 if (point&&ipoint<2){
3922 xyz[ipoint][1] = point->GetY();
3923 xyz[ipoint][2] = point->GetZ();
3926 xyz[ipoint][1] = cl->GetY();
3927 xyz[ipoint][2] = cl->GetZ();
3934 // Calculate seed state vector and covariance matrix
3936 Double_t alpha, cs,sn, xx2,yy2;
3938 alpha = (sec[1]-sec[2])*fSectors->GetAlpha();
3939 cs = TMath::Cos(alpha);
3940 sn = TMath::Sin(alpha);
3941 xx2= xyz[1][0]*cs-xyz[1][1]*sn;
3942 yy2= xyz[1][0]*sn+xyz[1][1]*cs;
3946 alpha = (sec[0]-sec[2])*fSectors->GetAlpha();
3947 cs = TMath::Cos(alpha);
3948 sn = TMath::Sin(alpha);
3949 xx2= xyz[0][0]*cs-xyz[0][1]*sn;
3950 yy2= xyz[0][0]*sn+xyz[0][1]*cs;
3956 Double_t x[5],c[15];
3960 x[4]=F1(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
3961 x[2]=F2(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
3962 x[3]=F3n(xyz[2][0],xyz[2][1],xyz[0][0],xyz[0][1],xyz[2][2],xyz[0][2],x[4]);
3964 Double_t sy =0.1, sz =0.1;
3966 Double_t sy1=0.2, sz1=0.2;
3967 Double_t sy2=0.2, sz2=0.2;
3970 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;
3971 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;
3972 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;
3973 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;
3974 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;
3975 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;
3977 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;
3978 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;
3979 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;
3980 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;
3985 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3986 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3987 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3988 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3989 c[13]=f30*sy1*f40+f32*sy2*f42;
3990 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3992 // Int_t row1 = fSectors->GetRowNumber(xyz[2][0]);
3993 AliTPCseed *seed=new AliTPCseed(xyz[2][0], sec[2]*fSectors->GetAlpha()+fSectors->GetAlphaShift(), x, c, 0);
3994 seed->SetLastPoint(row[2]);
3995 seed->SetFirstPoint(row[2]);
3996 for (Int_t i=row[0];i<row[2];i++){
3997 seed->SetClusterIndex(i, track->GetClusterIndex(i));
4005 void AliTPCtrackerMI::FindMultiMC(TObjArray * array, AliESDEvent */*esd*/, Int_t iter)
4008 // find multi tracks - THIS FUNCTION IS ONLY FOR DEBUG PURPOSES
4010 // Use AliTPCReconstructor::StreamLevel()>2 if you want to tune parameters - cuts
4012 // Two reasons to have multiple find tracks
4013 // 1. Curling tracks can be find more than once
4014 // 2. Splitted tracks
4015 // a.) Multiple seeding to increase tracking efficiency - (~ 100% reached)
4016 // b.) Edge effect on the sector boundaries
4019 // Algorithm done in 2 phases - because of CPU consumption
4020 // it is n^2 algorithm - for lead-lead 20000x20000 combination are investigated
4022 // Algorihm for curling tracks sign:
4023 // 1 phase -makes a very rough fast cuts to minimize combinatorics
4024 // a.) opposite sign
4025 // b.) one of the tracks - not pointing to the primary vertex -
4026 // c.) delta tan(theta)
4028 // 2 phase - calculates DCA between tracks - time consument
4033 // General cuts - for splitted tracks and for curling tracks
4035 const Float_t kMaxdPhi = 0.2; // maximal distance in phi
4037 // Curling tracks cuts
4042 Int_t nentries = array->GetEntriesFast();
4043 AliHelix *helixes = new AliHelix[nentries];
4044 Float_t *xm = new Float_t[nentries];
4045 Float_t *dz0 = new Float_t[nentries];
4046 Float_t *dz1 = new Float_t[nentries];
4052 // Find track COG in x direction - point with best defined parameters
4054 for (Int_t i=0;i<nentries;i++){
4055 AliTPCseed* track = (AliTPCseed*)array->At(i);
4056 if (!track) continue;
4057 track->SetCircular(0);
4058 new (&helixes[i]) AliHelix(*track);
4062 track->GetDZ(GetX(),GetY(),GetZ(),GetBz(),dz);
4065 for (Int_t icl=0; icl<160; icl++){
4066 AliTPCclusterMI * cl = track->GetClusterPointer(icl);
4072 if (ncl>0) xm[i]/=Float_t(ncl);
4074 TTreeSRedirector &cstream = *fDebugStreamer;
4076 for (Int_t i0=0;i0<nentries;i0++){
4077 AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
4078 if (!track0) continue;
4079 Float_t xc0 = helixes[i0].GetHelix(6);
4080 Float_t yc0 = helixes[i0].GetHelix(7);
4081 Float_t r0 = helixes[i0].GetHelix(8);
4082 Float_t rc0 = TMath::Sqrt(xc0*xc0+yc0*yc0);
4083 Float_t fi0 = TMath::ATan2(yc0,xc0);
4085 for (Int_t i1=i0+1;i1<nentries;i1++){
4086 AliTPCseed * track1 = (AliTPCseed*)array->At(i1);
4087 if (!track1) continue;
4088 Int_t lab0=track0->GetLabel();
4089 Int_t lab1=track1->GetLabel();
4090 if (TMath::Abs(lab0)!=TMath::Abs(lab1)) continue;
4092 Float_t xc1 = helixes[i1].GetHelix(6);
4093 Float_t yc1 = helixes[i1].GetHelix(7);
4094 Float_t r1 = helixes[i1].GetHelix(8);
4095 Float_t rc1 = TMath::Sqrt(xc1*xc1+yc1*yc1);
4096 Float_t fi1 = TMath::ATan2(yc1,xc1);
4098 Float_t dfi = fi0-fi1;
4101 if (dfi>1.5*TMath::Pi()) dfi-=TMath::Pi(); // take care about edge effect
4102 if (dfi<-1.5*TMath::Pi()) dfi+=TMath::Pi(); //
4103 if (TMath::Abs(dfi)>kMaxdPhi&&helixes[i0].GetHelix(4)*helixes[i1].GetHelix(4)<0){
4105 // if short tracks with undefined sign
4106 fi1 = -TMath::ATan2(yc1,-xc1);
4109 Float_t dtheta = TMath::Abs(track0->GetTgl()-track1->GetTgl())<TMath::Abs(track0->GetTgl()+track1->GetTgl())? track0->GetTgl()-track1->GetTgl():track0->GetTgl()+track1->GetTgl();
4112 // debug stream to tune "fast cuts"
4114 Double_t dist[3]; // distance at X
4115 Double_t mdist[3]={0,0,0}; // mean distance X+-40cm
4116 track0->GetDistance(track1,0.5*(xm[i0]+xm[i1])-40.,dist,AliTracker::GetBz());
4117 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[i]);
4118 track0->GetDistance(track1,0.5*(xm[i0]+xm[i1])+40.,dist,AliTracker::GetBz());
4119 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[i]);
4120 track0->GetDistance(track1,0.5*(xm[i0]+xm[i1]),dist,AliTracker::GetBz());
4121 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[i]);
4122 for (Int_t i=0;i<3;i++) mdist[i]*=0.33333;
4126 for (Int_t icl=0; icl<160; icl++){
4127 AliTPCclusterMI * cl0 = track0->GetClusterPointer(icl);
4128 AliTPCclusterMI * cl1 = track1->GetClusterPointer(icl);
4131 if (cl0==cl1) sums++;
4139 "Tr0.="<<track0<< // seed0
4140 "Tr1.="<<track1<< // seed1
4141 "h0.="<<&helixes[i0]<<
4142 "h1.="<<&helixes[i1]<<
4144 "sum="<<sum<< //the sum of rows with cl in both
4145 "sums="<<sums<< //the sum of shared clusters
4146 "xm0="<<xm[i0]<< // the center of track
4147 "xm1="<<xm[i1]<< // the x center of track
4148 // General cut variables
4149 "dfi="<<dfi<< // distance in fi angle
4150 "dtheta="<<dtheta<< // distance int theta angle
4156 "dist0="<<dist[0]<< //distance x
4157 "dist1="<<dist[1]<< //distance y
4158 "dist2="<<dist[2]<< //distance z
4159 "mdist0="<<mdist[0]<< //distance x
4160 "mdist1="<<mdist[1]<< //distance y
4161 "mdist2="<<mdist[2]<< //distance z
4174 if (AliTPCReconstructor::StreamLevel()>1) {
4175 AliInfo("Time for curling tracks removal DEBUGGING MC");
4181 void AliTPCtrackerMI::FindSplitted(TObjArray * array, AliESDEvent */*esd*/, Int_t iter)
4185 // Two reasons to have multiple find tracks
4186 // 1. Curling tracks can be find more than once
4187 // 2. Splitted tracks
4188 // a.) Multiple seeding to increase tracking efficiency - (~ 100% reached)
4189 // b.) Edge effect on the sector boundaries
4191 // This function tries to find tracks closed in the parametric space
4193 // cut logic if distance is bigger than cut continue - Do Nothing
4194 const Float_t kMaxdTheta = 0.05; // maximal distance in theta
4195 const Float_t kMaxdPhi = 0.05; // maximal deistance in phi
4196 const Float_t kdelta = 40.; //delta r to calculate track distance
4198 // const Float_t kMaxDist0 = 1.; // maximal distance 0
4199 //const Float_t kMaxDist1 = 0.3; // maximal distance 1 - cut if track in separate rows
4202 TCut csec("csec","abs(Tr0.fRelativeSector-Tr1.fRelativeSector)<2");
4203 TCut cdtheta("cdtheta","abs(dtheta)<0.05");
4208 Int_t nentries = array->GetEntriesFast();
4209 AliHelix *helixes = new AliHelix[nentries];
4210 Float_t *xm = new Float_t[nentries];
4216 //Sort tracks according quality
4218 Int_t nseed = array->GetEntriesFast();
4219 Float_t * quality = new Float_t[nseed];
4220 Int_t * indexes = new Int_t[nseed];
4221 for (Int_t i=0; i<nseed; i++) {
4222 AliTPCseed *pt=(AliTPCseed*)array->UncheckedAt(i);
4227 pt->UpdatePoints(); //select first last max dens points
4228 Float_t * points = pt->GetPoints();
4229 if (points[3]<0.8) quality[i] =-1;
4230 quality[i] = (points[2]-points[0])+pt->GetNumberOfClusters();
4231 //prefer high momenta tracks if overlaps
4232 quality[i] *= TMath::Sqrt(TMath::Abs(pt->Pt())+0.5);
4234 TMath::Sort(nseed,quality,indexes);
4238 // Find track COG in x direction - point with best defined parameters
4240 for (Int_t i=0;i<nentries;i++){
4241 AliTPCseed* track = (AliTPCseed*)array->At(i);
4242 if (!track) continue;
4243 track->SetCircular(0);
4244 new (&helixes[i]) AliHelix(*track);
4247 for (Int_t icl=0; icl<160; icl++){
4248 AliTPCclusterMI * cl = track->GetClusterPointer(icl);
4254 if (ncl>0) xm[i]/=Float_t(ncl);
4256 TTreeSRedirector &cstream = *fDebugStreamer;
4258 for (Int_t is0=0;is0<nentries;is0++){
4259 Int_t i0 = indexes[is0];
4260 AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
4261 if (!track0) continue;
4262 if (track0->GetKinkIndexes()[0]!=0) continue;
4263 Float_t xc0 = helixes[i0].GetHelix(6);
4264 Float_t yc0 = helixes[i0].GetHelix(7);
4265 Float_t fi0 = TMath::ATan2(yc0,xc0);
4267 for (Int_t is1=is0+1;is1<nentries;is1++){
4268 Int_t i1 = indexes[is1];
4269 AliTPCseed * track1 = (AliTPCseed*)array->At(i1);
4270 if (!track1) continue;
4272 if (TMath::Abs(track0->GetRelativeSector()-track1->GetRelativeSector())>1) continue;
4273 if (track1->GetKinkIndexes()[0]>0 &&track0->GetKinkIndexes()[0]<0) continue;
4274 if (track1->GetKinkIndexes()[0]!=0) continue;
4276 Float_t dtheta = TMath::Abs(track0->GetTgl()-track1->GetTgl())<TMath::Abs(track0->GetTgl()+track1->GetTgl())? track0->GetTgl()-track1->GetTgl():track0->GetTgl()+track1->GetTgl();
4277 if (TMath::Abs(dtheta)>kMaxdTheta) continue;
4279 Float_t xc1 = helixes[i1].GetHelix(6);
4280 Float_t yc1 = helixes[i1].GetHelix(7);
4281 Float_t fi1 = TMath::ATan2(yc1,xc1);
4283 Float_t dfi = fi0-fi1;
4284 if (dfi>1.5*TMath::Pi()) dfi-=TMath::Pi(); // take care about edge effect
4285 if (dfi<-1.5*TMath::Pi()) dfi+=TMath::Pi(); //
4286 if (TMath::Abs(dfi)>kMaxdPhi&&helixes[i0].GetHelix(4)*helixes[i1].GetHelix(4)<0){
4288 // if short tracks with undefined sign
4289 fi1 = -TMath::ATan2(yc1,-xc1);
4292 if (TMath::Abs(dfi)>kMaxdPhi) continue;
4299 for (Int_t icl=0; icl<160; icl++){
4300 Int_t index0=track0->GetClusterIndex2(icl);
4301 Int_t index1=track1->GetClusterIndex2(icl);
4302 Bool_t used0 = (index0>0 && !(index0&0x8000));
4303 Bool_t used1 = (index1>0 && !(index1&0x8000));
4305 if (used0) sum0++; // used cluster0
4306 if (used1) sum1++; // used clusters1
4307 if (used0&&used1) sum++;
4308 if (index0==index1 && used0 && used1) sums++;
4312 if (sums<10) continue;
4313 if (sum<40) continue;
4314 if (sums/Float_t(TMath::Min(sum0,sum1))<0.5) continue;
4316 Double_t dist[5][4]; // distance at X
4317 Double_t mdist[4]={0,0,0,0}; // mean distance on range +-delta
4321 track0->GetDistance(track1,xm[i0],dist[0],AliTracker::GetBz());
4322 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[0][i]);
4323 track0->GetDistance(track1,xm[i1],dist[1],AliTracker::GetBz());
4324 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[1][i]);
4326 track0->GetDistance(track1,TMath::Min(xm[i1],xm[i0])-kdelta,dist[2],AliTracker::GetBz());
4327 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[2][i]);
4328 track0->GetDistance(track1,TMath::Max(xm[i1],xm[i0])+kdelta,dist[3],AliTracker::GetBz());
4329 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[3][i]);
4331 track0->GetDistance(track1,(xm[i1]+xm[i0])*0.5,dist[4],AliTracker::GetBz());
4332 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[4][i]);
4333 for (Int_t i=0;i<3;i++) mdist[i]*=0.2;
4336 Int_t lab0=track0->GetLabel();
4337 Int_t lab1=track1->GetLabel();
4338 cstream<<"Splitted"<<
4342 "Tr0.="<<track0<< // seed0
4343 "Tr1.="<<track1<< // seed1
4344 "h0.="<<&helixes[i0]<<
4345 "h1.="<<&helixes[i1]<<
4347 "sum="<<sum<< //the sum of rows with cl in both
4348 "sum0="<<sum0<< //the sum of rows with cl in 0 track
4349 "sum1="<<sum1<< //the sum of rows with cl in 1 track
4350 "sums="<<sums<< //the sum of shared clusters
4351 "xm0="<<xm[i0]<< // the center of track
4352 "xm1="<<xm[i1]<< // the x center of track
4353 // General cut variables
4354 "dfi="<<dfi<< // distance in fi angle
4355 "dtheta="<<dtheta<< // distance int theta angle
4358 "dist0="<<dist[4][0]<< //distance x
4359 "dist1="<<dist[4][1]<< //distance y
4360 "dist2="<<dist[4][2]<< //distance z
4361 "mdist0="<<mdist[0]<< //distance x
4362 "mdist1="<<mdist[1]<< //distance y
4363 "mdist2="<<mdist[2]<< //distance z
4366 delete array->RemoveAt(i1);
4371 AliInfo("Time for splitted tracks removal");
4377 void AliTPCtrackerMI::FindCurling(TObjArray * array, AliESDEvent */*esd*/, Int_t iter)
4380 // find Curling tracks
4381 // Use AliTPCReconstructor::StreamLevel()>1 if you want to tune parameters - cuts
4384 // Algorithm done in 2 phases - because of CPU consumption
4385 // it is n^2 algorithm - for lead-lead 20000x20000 combination are investigated
4386 // see detal in MC part what can be used to cut
4390 const Float_t kMaxC = 400; // maximal curvature to of the track
4391 const Float_t kMaxdTheta = 0.15; // maximal distance in theta
4392 const Float_t kMaxdPhi = 0.15; // maximal distance in phi
4393 const Float_t kPtRatio = 0.3; // ratio between pt
4394 const Float_t kMinDCAR = 2.; // distance to the primary vertex in r - see cpipe cut
4397 // Curling tracks cuts
4400 const Float_t kMaxDeltaRMax = 40; // distance in outer radius
4401 const Float_t kMaxDeltaRMin = 5.; // distance in lower radius - see cpipe cut
4402 const Float_t kMinAngle = 2.9; // angle between tracks
4403 const Float_t kMaxDist = 5; // biggest distance
4405 // The cuts can be tuned using the "MC information stored in Multi tree ==> see FindMultiMC
4408 TCut csign("csign","Tr0.fP[4]*Tr1.fP[4]<0"); //opposite sign
4409 TCut cmax("cmax","abs(Tr0.GetC())>1/400");
4410 TCut cda("cda","sqrt(dtheta^2+dfi^2)<0.15");
4411 TCut ccratio("ccratio","abs((Tr0.fP[4]+Tr1.fP[4])/(abs(Tr0.fP[4])+abs(Tr1.fP[4])))<0.3");
4412 TCut cpipe("cpipe", "min(abs(r0-rc0),abs(r1-rc1))>5");
4414 TCut cdrmax("cdrmax","abs(abs(rc0+r0)-abs(rc1+r1))<40")
4415 TCut cdrmin("cdrmin","abs(abs(rc0+r0)-abs(rc1+r1))<10")
4417 Multi->Draw("dfi","iter==0"+csign+cmax+cda+ccratio); ~94% of curling tracks fulfill
4418 Multi->Draw("min(abs(r0-rc0),abs(r1-rc1))","iter==0&&abs(lab1)==abs(lab0)"+csign+cmax+cda+ccratio+cpipe+cdrmin+cdrmax); //80%
4420 Curling2->Draw("dfi","iter==0&&abs(lab0)==abs(lab1)"+csign+cmax+cdtheta+cdfi+ccratio)
4426 Int_t nentries = array->GetEntriesFast();
4427 AliHelix *helixes = new AliHelix[nentries];
4428 for (Int_t i=0;i<nentries;i++){
4429 AliTPCseed* track = (AliTPCseed*)array->At(i);
4430 if (!track) continue;
4431 track->SetCircular(0);
4432 new (&helixes[i]) AliHelix(*track);
4438 Double_t phase[2][2],radius[2];
4442 TTreeSRedirector &cstream = *fDebugStreamer;
4444 for (Int_t i0=0;i0<nentries;i0++){
4445 AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
4446 if (!track0) continue;
4447 if (TMath::Abs(track0->GetC())<1/kMaxC) continue;
4448 Float_t xc0 = helixes[i0].GetHelix(6);
4449 Float_t yc0 = helixes[i0].GetHelix(7);
4450 Float_t r0 = helixes[i0].GetHelix(8);
4451 Float_t rc0 = TMath::Sqrt(xc0*xc0+yc0*yc0);
4452 Float_t fi0 = TMath::ATan2(yc0,xc0);
4454 for (Int_t i1=i0+1;i1<nentries;i1++){
4455 AliTPCseed * track1 = (AliTPCseed*)array->At(i1);
4456 if (!track1) continue;
4457 if (TMath::Abs(track1->GetC())<1/kMaxC) continue;
4458 Float_t xc1 = helixes[i1].GetHelix(6);
4459 Float_t yc1 = helixes[i1].GetHelix(7);
4460 Float_t r1 = helixes[i1].GetHelix(8);
4461 Float_t rc1 = TMath::Sqrt(xc1*xc1+yc1*yc1);
4462 Float_t fi1 = TMath::ATan2(yc1,xc1);
4464 Float_t dfi = fi0-fi1;
4467 if (dfi>1.5*TMath::Pi()) dfi-=TMath::Pi(); // take care about edge effect
4468 if (dfi<-1.5*TMath::Pi()) dfi+=TMath::Pi(); //
4469 Float_t dtheta = TMath::Abs(track0->GetTgl()-track1->GetTgl())<TMath::Abs(track0->GetTgl()+track1->GetTgl())? track0->GetTgl()-track1->GetTgl():track0->GetTgl()+track1->GetTgl();
4473 if (track0->GetBConstrain()&&track1->GetBConstrain()) continue; // not constrained
4474 if (track1->GetSigned1Pt()*track0->GetSigned1Pt()>0) continue; // not the same sign
4475 if ( TMath::Abs(track1->GetTgl()+track0->GetTgl())>kMaxdTheta) continue; //distance in the Theta
4476 if ( TMath::Abs(dfi)>kMaxdPhi) continue; //distance in phi
4477 if ( TMath::Sqrt(dfi*dfi+dtheta*dtheta)>kMaxdPhi) continue; //common angular offset
4479 Float_t pt0 = track0->GetSignedPt();
4480 Float_t pt1 = track1->GetSignedPt();
4481 if ((TMath::Abs(pt0+pt1)/(TMath::Abs(pt0)+TMath::Abs(pt1)))>kPtRatio) continue;
4482 if ((iter==1) && TMath::Abs(TMath::Abs(rc0+r0)-TMath::Abs(rc1+r1))>kMaxDeltaRMax) continue;
4483 if ((iter!=1) &&TMath::Abs(TMath::Abs(rc0-r0)-TMath::Abs(rc1-r1))>kMaxDeltaRMin) continue;
4484 if (TMath::Min(TMath::Abs(rc0-r0),TMath::Abs(rc1-r1))<kMinDCAR) continue;
4487 // Now find closest approach
4491 Int_t npoints = helixes[i0].GetRPHIintersections(helixes[i1], phase, radius,10);
4492 if (npoints==0) continue;
4493 helixes[i0].GetClosestPhases(helixes[i1], phase);
4497 Double_t hangles[3];
4498 helixes[i0].Evaluate(phase[0][0],xyz0);
4499 helixes[i1].Evaluate(phase[0][1],xyz1);
4501 helixes[i0].GetAngle(phase[0][0],helixes[i1],phase[0][1],hangles);
4502 Double_t deltah[2],deltabest;
4503 if (TMath::Abs(hangles[2])<kMinAngle) continue;
4507 helixes[i0].ParabolicDCA(helixes[i1],phase[0][0],phase[0][1],radius[0],deltah[0],2);
4509 helixes[i0].ParabolicDCA(helixes[i1],phase[1][0],phase[1][1],radius[1],deltah[1],2);
4510 if (deltah[1]<deltah[0]) ibest=1;
4512 deltabest = TMath::Sqrt(deltah[ibest]);
4513 helixes[i0].Evaluate(phase[ibest][0],xyz0);
4514 helixes[i1].Evaluate(phase[ibest][1],xyz1);
4515 helixes[i0].GetAngle(phase[ibest][0],helixes[i1],phase[ibest][1],hangles);
4516 Double_t radiusbest = TMath::Sqrt(radius[ibest]);
4518 if (deltabest>kMaxDist) continue;
4519 // if (mindcar+mindcaz<40 && (TMath::Abs(hangles[2])<kMinAngle ||deltabest>3)) continue;
4520 Bool_t sign =kFALSE;
4521 if (hangles[2]>kMinAngle) sign =kTRUE;
4524 // circular[i0] = kTRUE;
4525 // circular[i1] = kTRUE;
4526 if (track0->OneOverPt()<track1->OneOverPt()){
4527 track0->SetCircular(track0->GetCircular()+1);
4528 track1->SetCircular(track1->GetCircular()+2);
4531 track1->SetCircular(track1->GetCircular()+1);
4532 track0->SetCircular(track0->GetCircular()+2);
4535 if (AliTPCReconstructor::StreamLevel()>1){
4537 //debug stream to tune "fine" cuts
4538 Int_t lab0=track0->GetLabel();
4539 Int_t lab1=track1->GetLabel();
4540 cstream<<"Curling2"<<
4556 "npoints="<<npoints<<
4557 "hangles0="<<hangles[0]<<
4558 "hangles1="<<hangles[1]<<
4559 "hangles2="<<hangles[2]<<
4562 "radius="<<radiusbest<<
4563 "deltabest="<<deltabest<<
4564 "phase0="<<phase[ibest][0]<<
4565 "phase1="<<phase[ibest][1]<<
4573 if (AliTPCReconstructor::StreamLevel()>1) {
4574 AliInfo("Time for curling tracks removal");
4583 void AliTPCtrackerMI::FindKinks(TObjArray * array, AliESDEvent *esd)
4590 TObjArray *kinks= new TObjArray(10000);
4591 // TObjArray *v0s= new TObjArray(10000);
4592 Int_t nentries = array->GetEntriesFast();
4593 AliHelix *helixes = new AliHelix[nentries];
4594 Int_t *sign = new Int_t[nentries];
4595 Int_t *nclusters = new Int_t[nentries];
4596 Float_t *alpha = new Float_t[nentries];
4597 AliKink *kink = new AliKink();
4598 Int_t * usage = new Int_t[nentries];
4599 Float_t *zm = new Float_t[nentries];
4600 Float_t *z0 = new Float_t[nentries];
4601 Float_t *fim = new Float_t[nentries];
4602 Float_t *shared = new Float_t[nentries];
4603 Bool_t *circular = new Bool_t[nentries];
4604 Float_t *dca = new Float_t[nentries];
4605 //const AliESDVertex * primvertex = esd->GetVertex();
4607 // nentries = array->GetEntriesFast();
4612 for (Int_t i=0;i<nentries;i++){
4615 AliTPCseed* track = (AliTPCseed*)array->At(i);
4616 if (!track) continue;
4617 track->SetCircular(0);
4619 track->UpdatePoints();
4620 if (( track->GetPoints()[2]- track->GetPoints()[0])>5 && track->GetPoints()[3]>0.8){
4622 nclusters[i]=track->GetNumberOfClusters();
4623 alpha[i] = track->GetAlpha();
4624 new (&helixes[i]) AliHelix(*track);
4626 helixes[i].Evaluate(0,xyz);
4627 sign[i] = (track->GetC()>0) ? -1:1;
4630 if (track->GetProlongation(x,y,z)){
4632 fim[i] = alpha[i]+TMath::ATan2(y,x);
4635 zm[i] = track->GetZ();
4639 circular[i]= kFALSE;
4640 if (track->GetProlongation(0,y,z)) z0[i] = z;
4641 dca[i] = track->GetD(0,0);
4647 Int_t ncandidates =0;
4650 Double_t phase[2][2],radius[2];
4653 // Find circling track
4654 TTreeSRedirector &cstream = *fDebugStreamer;
4656 for (Int_t i0=0;i0<nentries;i0++){
4657 AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
4658 if (!track0) continue;
4659 if (track0->GetNumberOfClusters()<40) continue;
4660 if (TMath::Abs(1./track0->GetC())>200) continue;
4661 for (Int_t i1=i0+1;i1<nentries;i1++){
4662 AliTPCseed * track1 = (AliTPCseed*)array->At(i1);
4663 if (!track1) continue;
4664 if (track1->GetNumberOfClusters()<40) continue;
4665 if ( TMath::Abs(track1->GetTgl()+track0->GetTgl())>0.1) continue;
4666 if (track0->GetBConstrain()&&track1->GetBConstrain()) continue;
4667 if (TMath::Abs(1./track1->GetC())>200) continue;
4668 if (track1->GetSigned1Pt()*track0->GetSigned1Pt()>0) continue;
4669 if (track1->GetTgl()*track0->GetTgl()>0) continue;
4670 if (TMath::Max(TMath::Abs(1./track0->GetC()),TMath::Abs(1./track1->GetC()))>190) continue;
4671 if (track0->GetBConstrain()&&track1->OneOverPt()<track0->OneOverPt()) continue; //returning - lower momenta
4672 if (track1->GetBConstrain()&&track0->OneOverPt()<track1->OneOverPt()) continue; //returning - lower momenta
4674 Float_t mindcar = TMath::Min(TMath::Abs(dca[i0]),TMath::Abs(dca[i1]));
4675 if (mindcar<5) continue;
4676 Float_t mindcaz = TMath::Min(TMath::Abs(z0[i0]-GetZ()),TMath::Abs(z0[i1]-GetZ()));
4677 if (mindcaz<5) continue;
4678 if (mindcar+mindcaz<20) continue;
4681 Float_t xc0 = helixes[i0].GetHelix(6);
4682 Float_t yc0 = helixes[i0].GetHelix(7);
4683 Float_t r0 = helixes[i0].GetHelix(8);
4684 Float_t xc1 = helixes[i1].GetHelix(6);
4685 Float_t yc1 = helixes[i1].GetHelix(7);
4686 Float_t r1 = helixes[i1].GetHelix(8);
4688 Float_t rmean = (r0+r1)*0.5;
4689 Float_t delta =TMath::Sqrt((xc1-xc0)*(xc1-xc0)+(yc1-yc0)*(yc1-yc0));
4690 //if (delta>30) continue;
4691 if (delta>rmean*0.25) continue;
4692 if (TMath::Abs(r0-r1)/rmean>0.3) continue;
4694 Int_t npoints = helixes[i0].GetRPHIintersections(helixes[i1], phase, radius,10);
4695 if (npoints==0) continue;
4696 helixes[i0].GetClosestPhases(helixes[i1], phase);
4700 Double_t hangles[3];
4701 helixes[i0].Evaluate(phase[0][0],xyz0);
4702 helixes[i1].Evaluate(phase[0][1],xyz1);
4704 helixes[i0].GetAngle(phase[0][0],helixes[i1],phase[0][1],hangles);
4705 Double_t deltah[2],deltabest;
4706 if (hangles[2]<2.8) continue;
4709 helixes[i0].ParabolicDCA(helixes[i1],phase[0][0],phase[0][1],radius[0],deltah[0],2);
4711 helixes[i0].ParabolicDCA(helixes[i1],phase[1][0],phase[1][1],radius[1],deltah[1],2);
4712 if (deltah[1]<deltah[0]) ibest=1;
4714 deltabest = TMath::Sqrt(deltah[ibest]);
4715 helixes[i0].Evaluate(phase[ibest][0],xyz0);
4716 helixes[i1].Evaluate(phase[ibest][1],xyz1);
4717 helixes[i0].GetAngle(phase[ibest][0],helixes[i1],phase[ibest][1],hangles);
4718 Double_t radiusbest = TMath::Sqrt(radius[ibest]);
4720 if (deltabest>6) continue;
4721 if (mindcar+mindcaz<40 && (hangles[2]<3.12||deltabest>3)) continue;
4722 Bool_t sign =kFALSE;
4723 if (hangles[2]>3.06) sign =kTRUE;
4726 circular[i0] = kTRUE;
4727 circular[i1] = kTRUE;
4728 if (track0->OneOverPt()<track1->OneOverPt()){
4729 track0->SetCircular(track0->GetCircular()+1);
4730 track1->SetCircular(track1->GetCircular()+2);
4733 track1->SetCircular(track1->GetCircular()+1);
4734 track0->SetCircular(track0->GetCircular()+2);
4737 if (sign&&AliTPCReconstructor::StreamLevel()>1){
4739 Int_t lab0=track0->GetLabel();
4740 Int_t lab1=track1->GetLabel();
4741 cstream<<"Curling"<<
4748 "mindcar="<<mindcar<<
4749 "mindcaz="<<mindcaz<<
4752 "npoints="<<npoints<<
4753 "hangles0="<<hangles[0]<<
4754 "hangles2="<<hangles[2]<<
4759 "radius="<<radiusbest<<
4760 "deltabest="<<deltabest<<
4761 "phase0="<<phase[ibest][0]<<
4762 "phase1="<<phase[ibest][1]<<
4772 for (Int_t i =0;i<nentries;i++){
4773 if (sign[i]==0) continue;
4774 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
4781 Double_t cradius0 = 40*40;
4782 Double_t cradius1 = 270*270;
4785 Double_t cdist3=0.55;
4786 for (Int_t j =i+1;j<nentries;j++){
4788 if (sign[j]*sign[i]<1) continue;
4789 if ( (nclusters[i]+nclusters[j])>200) continue;
4790 if ( (nclusters[i]+nclusters[j])<80) continue;
4791 if ( TMath::Abs(zm[i]-zm[j])>60.) continue;
4792 if ( TMath::Abs(fim[i]-fim[j])>0.6 && TMath::Abs(fim[i]-fim[j])<5.7 ) continue;
4793 //AliTPCseed * track1 = (AliTPCseed*)array->At(j); Double_t phase[2][2],radius[2];
4794 Int_t npoints = helixes[i].GetRPHIintersections(helixes[j], phase, radius,20);
4795 if (npoints<1) continue;
4798 if (radius[0]<cradius0||radius[0]>cradius1) continue;
4801 if ( (radius[0]<cradius0||radius[0]>cradius1) && (radius[1]<cradius0||radius[1]>cradius1) ) continue;
4804 Double_t delta1=10000,delta2=10000;
4805 // cuts on the intersection radius
4806 helixes[i].LinearDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta1);
4807 if (radius[0]<20&&delta1<1) continue; //intersection at vertex
4808 if (radius[0]<10&&delta1<3) continue; //intersection at vertex
4810 helixes[i].LinearDCA(helixes[j],phase[1][0],phase[1][1],radius[1],delta2);
4811 if (radius[1]<20&&delta2<1) continue; //intersection at vertex
4812 if (radius[1]<10&&delta2<3) continue; //intersection at vertex
4815 Double_t distance1 = TMath::Min(delta1,delta2);
4816 if (distance1>cdist1) continue; // cut on DCA linear approximation
4818 npoints = helixes[i].GetRPHIintersections(helixes[j], phase, radius,20);
4819 helixes[i].ParabolicDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta1);
4820 if (radius[0]<20&&delta1<1) continue; //intersection at vertex
4821 if (radius[0]<10&&delta1<3) continue; //intersection at vertex
4824 helixes[i].ParabolicDCA(helixes[j],phase[1][0],phase[1][1],radius[1],delta2);
4825 if (radius[1]<20&&delta2<1) continue; //intersection at vertex
4826 if (radius[1]<10&&delta2<3) continue; //intersection at vertex
4828 distance1 = TMath::Min(delta1,delta2);
4831 rkink = TMath::Sqrt(radius[0]);
4834 rkink = TMath::Sqrt(radius[1]);
4836 if (distance1>cdist2) continue;
4839 AliTPCseed * track1 = (AliTPCseed*)array->At(j);
4842 Int_t row0 = GetRowNumber(rkink);
4843 if (row0<10) continue;
4844 if (row0>150) continue;
4847 Float_t dens00=-1,dens01=-1;
4848 Float_t dens10=-1,dens11=-1;
4850 Int_t found,foundable,shared;
4851 track0->GetClusterStatistic(0,row0-5, found, foundable,shared,kFALSE);
4852 if (foundable>5) dens00 = Float_t(found)/Float_t(foundable);
4853 track0->GetClusterStatistic(row0+5,155, found, foundable,shared,kFALSE);
4854 if (foundable>5) dens01 = Float_t(found)/Float_t(foundable);
4856 track1->GetClusterStatistic(0,row0-5, found, foundable,shared,kFALSE);
4857 if (foundable>10) dens10 = Float_t(found)/Float_t(foundable);
4858 track1->GetClusterStatistic(row0+5,155, found, foundable,shared,kFALSE);
4859 if (foundable>10) dens11 = Float_t(found)/Float_t(foundable);
4861 if (dens00<dens10 && dens01<dens11) continue;
4862 if (dens00>dens10 && dens01>dens11) continue;
4863 if (TMath::Max(dens00,dens10)<0.1) continue;
4864 if (TMath::Max(dens01,dens11)<0.3) continue;
4866 if (TMath::Min(dens00,dens10)>0.6) continue;
4867 if (TMath::Min(dens01,dens11)>0.6) continue;
4870 AliTPCseed * ktrack0, *ktrack1;
4879 if (TMath::Abs(ktrack0->GetC())>5) continue; // cut on the curvature for mother particle
4880 AliExternalTrackParam paramm(*ktrack0);
4881 AliExternalTrackParam paramd(*ktrack1);
4882 if (row0>60&&ktrack1->GetReference().GetX()>90.)new (¶md) AliExternalTrackParam(ktrack1->GetReference());
4885 kink->SetMother(paramm);
4886 kink->SetDaughter(paramd);
4889 Float_t x[3] = { kink->GetPosition()[0],kink->GetPosition()[1],kink->GetPosition()[2]};
4891 fParam->Transform0to1(x,index);
4892 fParam->Transform1to2(x,index);
4893 row0 = GetRowNumber(x[0]);
4895 if (kink->GetR()<100) continue;
4896 if (kink->GetR()>240) continue;
4897 if (kink->GetPosition()[2]/kink->GetR()>AliTPCReconstructor::GetCtgRange()) continue; //out of fiducial volume
4898 if (kink->GetDistance()>cdist3) continue;
4899 Float_t dird = kink->GetDaughterP()[0]*kink->GetPosition()[0]+kink->GetDaughterP()[1]*kink->GetPosition()[1]; // rough direction estimate
4900 if (dird<0) continue;
4902 Float_t dirm = kink->GetMotherP()[0]*kink->GetPosition()[0]+kink->GetMotherP()[1]*kink->GetPosition()[1]; // rough direction estimate
4903 if (dirm<0) continue;
4904 Float_t mpt = TMath::Sqrt(kink->GetMotherP()[0]*kink->GetMotherP()[0]+kink->GetMotherP()[1]*kink->GetMotherP()[1]);
4905 if (mpt<0.2) continue;
4908 //for high momenta momentum not defined well in first iteration
4909 Double_t qt = TMath::Sin(kink->GetAngle(2))*ktrack1->GetP();
4910 if (qt>0.35) continue;
4913 kink->SetLabel(CookLabel(ktrack0,0.4,0,row0),0);
4914 kink->SetLabel(CookLabel(ktrack1,0.4,row0,160),1);
4916 kink->SetTPCDensity(dens00,0,0);
4917 kink->SetTPCDensity(dens01,0,1);
4918 kink->SetTPCDensity(dens10,1,0);
4919 kink->SetTPCDensity(dens11,1,1);
4920 kink->SetIndex(i,0);
4921 kink->SetIndex(j,1);
4924 kink->SetTPCDensity(dens10,0,0);
4925 kink->SetTPCDensity(dens11,0,1);
4926 kink->SetTPCDensity(dens00,1,0);
4927 kink->SetTPCDensity(dens01,1,1);
4928 kink->SetIndex(j,0);
4929 kink->SetIndex(i,1);
4932 if (mpt<1||kink->GetAngle(2)>0.1){
4933 // angle and densities not defined yet
4934 if (kink->GetTPCDensityFactor()<0.8) continue;
4935 if ((2-kink->GetTPCDensityFactor())*kink->GetDistance() >0.25) continue;
4936 if (kink->GetAngle(2)*ktrack0->GetP()<0.003) continue; //too small angle
4937 if (kink->GetAngle(2)>0.2&&kink->GetTPCDensityFactor()<1.15) continue;
4938 if (kink->GetAngle(2)>0.2&&kink->GetTPCDensity(0,1)>0.05) continue;
4940 Float_t criticalangle = track0->GetSigmaSnp2()+track0->GetSigmaTgl2();
4941 criticalangle+= track1->GetSigmaSnp2()+track1->GetSigmaTgl2();
4942 criticalangle= 3*TMath::Sqrt(criticalangle);
4943 if (criticalangle>0.02) criticalangle=0.02;
4944 if (kink->GetAngle(2)<criticalangle) continue;
4947 Int_t drow = Int_t(2.+0.5/(0.05+kink->GetAngle(2))); // overlap region defined
4948 Float_t shapesum =0;
4950 for ( Int_t row = row0-drow; row<row0+drow;row++){
4951 if (row<0) continue;
4952 if (row>155) continue;
4953 if (ktrack0->GetClusterPointer(row)){
4954 AliTPCTrackerPoint *point =ktrack0->GetTrackPoint(row);
4955 shapesum+=point->GetSigmaY()+point->GetSigmaZ();
4958 if (ktrack1->GetClusterPointer(row)){
4959 AliTPCTrackerPoint *point =ktrack1->GetTrackPoint(row);
4960 shapesum+=point->GetSigmaY()+point->GetSigmaZ();
4965 kink->SetShapeFactor(-1.);
4968 kink->SetShapeFactor(shapesum/sum);
4970 // esd->AddKink(kink);
4971 kinks->AddLast(kink);
4977 // sort the kinks according quality - and refit them towards vertex
4979 Int_t nkinks = kinks->GetEntriesFast();
4980 Float_t *quality = new Float_t[nkinks];
4981 Int_t *indexes = new Int_t[nkinks];
4982 AliTPCseed *mothers = new AliTPCseed[nkinks];
4983 AliTPCseed *daughters = new AliTPCseed[nkinks];
4986 for (Int_t i=0;i<nkinks;i++){
4988 AliKink *kink = (AliKink*)kinks->At(i);
4990 // refit kinks towards vertex
4992 Int_t index0 = kink->GetIndex(0);
4993 Int_t index1 = kink->GetIndex(1);
4994 AliTPCseed * ktrack0 = (AliTPCseed*)array->At(index0);
4995 AliTPCseed * ktrack1 = (AliTPCseed*)array->At(index1);
4997 Int_t sumn=ktrack0->GetNumberOfClusters()+ktrack1->GetNumberOfClusters();
4999 // Refit Kink under if too small angle
5001 if (kink->GetAngle(2)<0.05){
5002 kink->SetTPCRow0(GetRowNumber(kink->GetR()));
5003 Int_t row0 = kink->GetTPCRow0();
5004 Int_t drow = Int_t(2.+0.5/(0.05+kink->GetAngle(2)));
5007 Int_t last = row0-drow;
5008 if (last<40) last=40;
5009 if (last<ktrack0->GetFirstPoint()+25) last = ktrack0->GetFirstPoint()+25;
5010 AliTPCseed* seed0 = ReSeed(ktrack0,last,kFALSE);
5013 Int_t first = row0+drow;
5014 if (first>130) first=130;
5015 if (first>ktrack1->GetLastPoint()-25) first = TMath::Max(ktrack1->GetLastPoint()-25,30);
5016 AliTPCseed* seed1 = ReSeed(ktrack1,first,kTRUE);
5018 if (seed0 && seed1){
5019 kink->SetStatus(1,8);
5020 if (RefitKink(*seed0,*seed1,*kink)) kink->SetStatus(1,9);
5021 row0 = GetRowNumber(kink->GetR());
5022 sumn = seed0->GetNumberOfClusters()+seed1->GetNumberOfClusters();
5023 mothers[i] = *seed0;
5024 daughters[i] = *seed1;
5027 delete kinks->RemoveAt(i);
5028 if (seed0) delete seed0;
5029 if (seed1) delete seed1;
5032 if (kink->GetDistance()>0.5 || kink->GetR()<110 || kink->GetR()>240) {
5033 delete kinks->RemoveAt(i);
5034 if (seed0) delete seed0;
5035 if (seed1) delete seed1;
5043 if (kink) quality[i] = 160*((0.1+kink->GetDistance())*(2.-kink->GetTPCDensityFactor()))/(sumn+40.); //the longest -clossest will win
5045 TMath::Sort(nkinks,quality,indexes,kFALSE);
5047 //remove double find kinks
5049 for (Int_t ikink0=1;ikink0<nkinks;ikink0++){
5050 AliKink * kink0 = (AliKink*) kinks->At(indexes[ikink0]);
5051 if (!kink0) continue;
5053 for (Int_t ikink1=0;ikink1<ikink0;ikink1++){
5054 if (!kink0) continue;
5055 AliKink * kink1 = (AliKink*) kinks->At(indexes[ikink1]);
5056 if (!kink1) continue;
5057 // if not close kink continue
5058 if (TMath::Abs(kink1->GetPosition()[2]-kink0->GetPosition()[2])>10) continue;
5059 if (TMath::Abs(kink1->GetPosition()[1]-kink0->GetPosition()[1])>10) continue;
5060 if (TMath::Abs(kink1->GetPosition()[0]-kink0->GetPosition()[0])>10) continue;
5062 AliTPCseed &mother0 = mothers[indexes[ikink0]];
5063 AliTPCseed &daughter0 = daughters[indexes[ikink0]];
5064 AliTPCseed &mother1 = mothers[indexes[ikink1]];
5065 AliTPCseed &daughter1 = daughters[indexes[ikink1]];
5066 Int_t row0 = (kink0->GetTPCRow0()+kink1->GetTPCRow0())/2;
5075 for (Int_t i=0;i<row0;i++){
5076 if (mother0.GetClusterIndex(i)>0 && mother1.GetClusterIndex(i)>0){
5079 if (mother0.GetClusterIndex(i)==mother1.GetClusterIndex(i)){
5086 for (Int_t i=row0;i<158;i++){
5087 if (daughter0.GetClusterIndex(i)>0 && daughter0.GetClusterIndex(i)>0){
5090 if (mother0.GetClusterIndex(i)==mother1.GetClusterIndex(i)){
5096 Float_t ratio = Float_t(same+1)/Float_t(both+1);
5097 Float_t ratiom = Float_t(samem+1)/Float_t(bothm+1);
5098 Float_t ratiod = Float_t(samed+1)/Float_t(bothd+1);
5099 if (ratio>0.3 && ratiom>0.5 &&ratiod>0.5) {
5100 Int_t sum0 = mother0.GetNumberOfClusters()+daughter0.GetNumberOfClusters();
5101 Int_t sum1 = mother1.GetNumberOfClusters()+daughter1.GetNumberOfClusters();
5103 shared[kink0->GetIndex(0)]= kTRUE;
5104 shared[kink0->GetIndex(1)]= kTRUE;
5105 delete kinks->RemoveAt(indexes[ikink0]);
5108 shared[kink1->GetIndex(0)]= kTRUE;
5109 shared[kink1->GetIndex(1)]= kTRUE;
5110 delete kinks->RemoveAt(indexes[ikink1]);
5117 for (Int_t i=0;i<nkinks;i++){
5118 AliKink * kink = (AliKink*) kinks->At(indexes[i]);
5119 if (!kink) continue;
5120 kink->SetTPCRow0(GetRowNumber(kink->GetR()));
5121 Int_t index0 = kink->GetIndex(0);
5122 Int_t index1 = kink->GetIndex(1);
5123 if (circular[index0]||circular[index1]&&kink->GetDistance()>0.2) continue;
5124 kink->SetMultiple(usage[index0],0);
5125 kink->SetMultiple(usage[index1],1);
5126 if (kink->GetMultiple()[0]+kink->GetMultiple()[1]>2) continue;
5127 if (kink->GetMultiple()[0]+kink->GetMultiple()[1]>0 && quality[indexes[i]]>0.2) continue;
5128 if (kink->GetMultiple()[0]+kink->GetMultiple()[1]>0 && kink->GetDistance()>0.2) continue;
5129 if (circular[index0]||circular[index1]&&kink->GetDistance()>0.1) continue;
5131 AliTPCseed * ktrack0 = (AliTPCseed*)array->At(index0);
5132 AliTPCseed * ktrack1 = (AliTPCseed*)array->At(index1);
5133 if (!ktrack0 || !ktrack1) continue;
5134 Int_t index = esd->AddKink(kink);
5137 if ( ktrack0->GetKinkIndex(0)==0 && ktrack1->GetKinkIndex(0)==0) { //best kink
5138 if (mothers[indexes[i]].GetNumberOfClusters()>20 && daughters[indexes[i]].GetNumberOfClusters()>20 && (mothers[indexes[i]].GetNumberOfClusters()+daughters[indexes[i]].GetNumberOfClusters())>100){
5139 *ktrack0 = mothers[indexes[i]];
5140 *ktrack1 = daughters[indexes[i]];
5144 ktrack0->SetKinkIndex(usage[index0],-(index+1));
5145 ktrack1->SetKinkIndex(usage[index1], (index+1));
5150 // Remove tracks corresponding to shared kink's
5152 for (Int_t i=0;i<nentries;i++){
5153 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5154 if (!track0) continue;
5155 if (track0->GetKinkIndex(0)!=0) continue;
5156 if (shared[i]) delete array->RemoveAt(i);
5161 RemoveUsed2(array,0.5,0.4,30);
5163 for (Int_t i=0;i<nentries;i++){
5164 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5165 if (!track0) continue;
5166 track0->CookdEdx(0.02,0.6);
5170 for (Int_t i=0;i<nentries;i++){
5171 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5172 if (!track0) continue;
5173 if (track0->Pt()<1.4) continue;
5174 //remove double high momenta tracks - overlapped with kink candidates
5177 for (Int_t icl=track0->GetFirstPoint();icl<track0->GetLastPoint(); icl++){
5178 if (track0->GetClusterPointer(icl)!=0){
5180 if (track0->GetClusterPointer(icl)->IsUsed(10)) shared++;
5183 if (Float_t(shared+1)/Float_t(all+1)>0.5) {
5184 delete array->RemoveAt(i);
5188 if (track0->GetKinkIndex(0)!=0) continue;
5189 if (track0->GetNumberOfClusters()<80) continue;
5191 AliTPCseed *pmother = new AliTPCseed();
5192 AliTPCseed *pdaughter = new AliTPCseed();
5193 AliKink *pkink = new AliKink;
5195 AliTPCseed & mother = *pmother;
5196 AliTPCseed & daughter = *pdaughter;
5197 AliKink & kink = *pkink;
5198 if (CheckKinkPoint(track0,mother,daughter, kink)){
5199 if (mother.GetNumberOfClusters()<30||daughter.GetNumberOfClusters()<20) {
5203 continue; //too short tracks
5205 if (mother.Pt()<1.4) {
5211 Int_t row0= kink.GetTPCRow0();
5212 if (kink.GetDistance()>0.5 || kink.GetR()<110. || kink.GetR()>240.) {
5219 Int_t index = esd->AddKink(&kink);
5220 mother.SetKinkIndex(0,-(index+1));
5221 daughter.SetKinkIndex(0,index+1);
5222 if (mother.GetNumberOfClusters()>50) {
5223 delete array->RemoveAt(i);
5224 array->AddAt(new AliTPCseed(mother),i);
5227 array->AddLast(new AliTPCseed(mother));
5229 array->AddLast(new AliTPCseed(daughter));
5230 for (Int_t icl=0;icl<row0;icl++) {
5231 if (mother.GetClusterPointer(icl)) mother.GetClusterPointer(icl)->Use(20);
5234 for (Int_t icl=row0;icl<158;icl++) {
5235 if (daughter.GetClusterPointer(icl)) daughter.GetClusterPointer(icl)->Use(20);
5244 delete [] daughters;
5266 printf("Ncandidates=\t%d\t%d\t%d\t%d\n",esd->GetNumberOfKinks(),ncandidates,ntracks,nall);
5270 void AliTPCtrackerMI::FindV0s(TObjArray * array, AliESDEvent *esd)
5276 TObjArray *tpcv0s = new TObjArray(100000);
5277 Int_t nentries = array->GetEntriesFast();
5278 AliHelix *helixes = new AliHelix[nentries];
5279 Int_t *sign = new Int_t[nentries];
5280 Float_t *alpha = new Float_t[nentries];
5281 Float_t *z0 = new Float_t[nentries];
5282 Float_t *dca = new Float_t[nentries];
5283 Float_t *sdcar = new Float_t[nentries];
5284 Float_t *cdcar = new Float_t[nentries];
5285 Float_t *pulldcar = new Float_t[nentries];
5286 Float_t *pulldcaz = new Float_t[nentries];
5287 Float_t *pulldca = new Float_t[nentries];
5288 Bool_t *isPrim = new Bool_t[nentries];
5289 const AliESDVertex * primvertex = esd->GetVertex();
5290 Double_t zvertex = primvertex->GetZv();
5292 // nentries = array->GetEntriesFast();
5294 for (Int_t i=0;i<nentries;i++){
5297 AliTPCseed* track = (AliTPCseed*)array->At(i);
5298 if (!track) continue;
5299 track->GetV0Indexes()[0] = 0; //rest v0 indexes
5300 track->GetV0Indexes()[1] = 0; //rest v0 indexes
5301 track->GetV0Indexes()[2] = 0; //rest v0 indexes
5303 alpha[i] = track->GetAlpha();
5304 new (&helixes[i]) AliHelix(*track);
5306 helixes[i].Evaluate(0,xyz);
5307 sign[i] = (track->GetC()>0) ? -1:1;
5311 if (track->GetProlongation(0,y,z)) z0[i] = z;
5312 dca[i] = track->GetD(0,0);
5314 // dca error parrameterezation + pulls
5316 sdcar[i] = TMath::Sqrt(0.150*0.150+(100*track->GetC())*(100*track->GetC()));
5317 if (TMath::Abs(track->GetTgl())>1) sdcar[i]*=2.5;
5318 cdcar[i] = TMath::Exp((TMath::Abs(track->GetC())-0.0106)*525.3);
5319 pulldcar[i] = (dca[i]-cdcar[i])/sdcar[i];
5320 pulldcaz[i] = (z0[i]-zvertex)/sdcar[i];
5321 pulldca[i] = TMath::Sqrt(pulldcar[i]*pulldcar[i]+pulldcaz[i]*pulldcaz[i]);
5322 if (track->TPCrPID(1)+track->TPCrPID(2)+track->TPCrPID(3)>0.5) {
5323 if (pulldca[i]<3.) isPrim[i]=kTRUE; //pion, muon and Kaon 3 sigma cut
5325 if (track->TPCrPID(4)>0.5) {
5326 if (pulldca[i]<0.5) isPrim[i]=kTRUE; //proton 0.5 sigma cut
5328 if (track->TPCrPID(0)>0.4) {
5329 isPrim[i]=kFALSE; //electron no sigma cut
5336 Int_t ncandidates =0;
5339 Double_t phase[2][2],radius[2];
5345 TTreeSRedirector &cstream = *fDebugStreamer;
5346 Float_t fprimvertex[3]={GetX(),GetY(),GetZ()};
5348 Double_t cradius0 = 10*10;
5349 Double_t cradius1 = 200*200;
5352 Double_t cpointAngle = 0.95;
5354 Double_t delta[2]={10000,10000};
5355 for (Int_t i =0;i<nentries;i++){
5356 if (sign[i]==0) continue;
5357 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5358 if (!track0) continue;
5359 if (AliTPCReconstructor::StreamLevel()>1){
5364 "zvertex="<<zvertex<<
5365 "sdcar0="<<sdcar[i]<<
5366 "cdcar0="<<cdcar[i]<<
5367 "pulldcar0="<<pulldcar[i]<<
5368 "pulldcaz0="<<pulldcaz[i]<<
5369 "pulldca0="<<pulldca[i]<<
5370 "isPrim="<<isPrim[i]<<
5374 if (track0->GetSigned1Pt()<0) continue;
5375 if (track0->GetKinkIndex(0)>0||isPrim[i]) continue; //daughter kink
5377 if (TMath::Abs(helixes[i].GetHelix(4))<0.000000001) continue;
5382 for (Int_t j =0;j<nentries;j++){
5383 AliTPCseed * track1 = (AliTPCseed*)array->At(j);
5384 if (!track1) continue;
5385 if (track1->GetKinkIndex(0)>0 || isPrim[j]) continue; //daughter kink
5386 if (sign[j]*sign[i]>0) continue;
5387 if (TMath::Abs(helixes[j].GetHelix(4))<0.000001) continue;
5388 if (track0->GetCircular()+track1->GetCircular()>1) continue; //circling -returning track
5391 // DCA to prim vertex cut
5397 Int_t npoints = helixes[i].GetRPHIintersections(helixes[j], phase, radius,cdist2);
5398 if (npoints<1) continue;
5402 if (radius[0]<cradius0||radius[0]>cradius1) continue;
5403 helixes[i].LinearDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta[0]);
5404 if (delta[0]>cdist1) continue;
5407 if (TMath::Max(radius[0],radius[1])<cradius0|| TMath::Min(radius[0],radius[1])>cradius1) continue;
5408 helixes[i].LinearDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta[0]);
5409 helixes[i].LinearDCA(helixes[j],phase[1][0],phase[1][1],radius[1],delta[1]);
5410 if (delta[1]<delta[0]) iclosest=1;
5411 if (delta[iclosest]>cdist1) continue;
5413 helixes[i].ParabolicDCA(helixes[j],phase[iclosest][0],phase[iclosest][1],radius[iclosest],delta[iclosest]);
5414 if (radius[iclosest]<cradius0 || radius[iclosest]>cradius1 || delta[iclosest]>cdist1) continue;
5416 Double_t pointAngle = helixes[i].GetPointAngle(helixes[j],phase[iclosest],fprimvertex);
5417 if (pointAngle<cpointAngle) continue;
5419 Bool_t isGamma = kFALSE;
5420 vertex.SetParamP(*track0); //track0 - plus
5421 vertex.SetParamN(*track1); //track1 - minus
5422 vertex.Update(fprimvertex);
5423 if (track0->TPCrPID(0)>0.3&&track1->TPCrPID(0)>0.3&&vertex.GetAnglep()[2]<0.15) isGamma=kTRUE; // gamma conversion candidate
5424 Double_t pointAngle2 = vertex.GetV0CosineOfPointingAngle();
5426 if (vertex.GetV0CosineOfPointingAngle()<cpointAngle && (!isGamma)) continue;// point angle cut
5427 //Bo: if (vertex.GetDist2()>2&&(!isGamma)) continue; // point angle cut
5428 if (vertex.GetDcaV0Daughters()>2&&(!isGamma)) continue;//Bo: // point angle cut
5429 Float_t sigmae = 0.15*0.15;
5430 if (vertex.GetRr()<80)
5431 sigmae += (sdcar[i]*sdcar[i]+sdcar[j]*sdcar[j])*(1.-vertex.GetRr()/80.)*(1.-vertex.GetRr()/80.);
5432 sigmae+= TMath::Sqrt(sigmae);
5433 //Bo: if (vertex.GetDist2()/sigmae>3.&&(!isGamma)) continue;
5434 if (vertex.GetDcaV0Daughters()/sigmae>3.&&(!isGamma)) continue;
5435 Float_t densb0=0,densb1=0,densa0=0,densa1=0;
5436 Int_t row0 = GetRowNumber(vertex.GetRr());
5438 //Bo: if (vertex.GetDist2()>0.2) continue;
5439 if (vertex.GetDcaV0Daughters()>0.2) continue;
5440 densb0 = track0->Density2(0,row0-5);
5441 densb1 = track1->Density2(0,row0-5);
5442 if (densb0>0.3|| densb1>0.3) continue; //clusters before vertex
5443 densa0 = track0->Density2(row0+5,row0+40);
5444 densa1 = track1->Density2(row0+5,row0+40);
5445 if ((densa0<0.4|| densa1<0.4)&&(!isGamma)) continue; //missing clusters after vertex
5448 densa0 = track0->Density2(0,40); //cluster density
5449 densa1 = track1->Density2(0,40); //cluster density
5450 if ((vertex.GetRr()<80&&densa0+densa1<1.)&&(!isGamma)) continue;
5452 //Bo: vertex.SetLab(0,track0->GetLabel());
5453 //Bo: vertex.SetLab(1,track1->GetLabel());
5454 vertex.SetChi2After((densa0+densa1)*0.5);
5455 vertex.SetChi2Before((densb0+densb1)*0.5);
5456 vertex.SetIndex(0,i);
5457 vertex.SetIndex(1,j);
5458 //Bo: vertex.SetStatus(1); // TPC v0 candidate
5459 vertex.SetOnFlyStatus(2);//Bo: // TPC v0 candidate
5460 //Bo: vertex.SetRp(track0->TPCrPIDs());
5461 //Bo: vertex.SetRm(track1->TPCrPIDs());
5462 tpcv0s->AddLast(new AliESDv0(vertex));
5465 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
5466 Double_t radiusm= (delta[0]<delta[1])? TMath::Sqrt(radius[0]):TMath::Sqrt(radius[1]);
5467 Double_t deltam= (delta[0]<delta[1])? TMath::Sqrt(delta[0]):TMath::Sqrt(delta[1]);
5468 if (AliTPCReconstructor::StreamLevel()>1) {
5469 Int_t lab0=track0->GetLabel();
5470 Int_t lab1=track1->GetLabel();
5471 Char_t c0=track0->GetCircular();
5472 Char_t c1=track1->GetCircular();
5475 "vertex.="<<&vertex<<
5478 "Helix0.="<<&helixes[i]<<
5481 "Helix1.="<<&helixes[j]<<
5482 "pointAngle="<<pointAngle<<
5483 "pointAngle2="<<pointAngle2<<
5488 "zvertex="<<zvertex<<
5491 "npoints="<<npoints<<
5492 "radius0="<<radius[0]<<
5493 "delta0="<<delta[0]<<
5494 "radius1="<<radius[1]<<
5495 "delta1="<<delta[1]<<
5496 "radiusm="<<radiusm<<
5498 "sdcar0="<<sdcar[i]<<
5499 "sdcar1="<<sdcar[j]<<
5500 "cdcar0="<<cdcar[i]<<
5501 "cdcar1="<<cdcar[j]<<
5502 "pulldcar0="<<pulldcar[i]<<
5503 "pulldcar1="<<pulldcar[j]<<
5504 "pulldcaz0="<<pulldcaz[i]<<
5505 "pulldcaz1="<<pulldcaz[j]<<
5506 "pulldca0="<<pulldca[i]<<
5507 "pulldca1="<<pulldca[j]<<
5517 Float_t *quality = new Float_t[ncandidates];
5518 Int_t *indexes = new Int_t[ncandidates];
5520 for (Int_t i=0;i<ncandidates;i++){
5522 AliESDv0 *v0 = (AliESDv0*)tpcv0s->At(i);
5523 quality[i] = 1./(1.00001-v0->GetV0CosineOfPointingAngle()); //base point angle
5524 // quality[i] /= (0.5+v0->GetDist2());
5525 // quality[i] *= v0->GetChi2After(); //density factor
5527 Int_t index0 = v0->GetIndex(0);
5528 Int_t index1 = v0->GetIndex(1);
5529 //Bo: Double_t minpulldca = TMath::Min(2.+pulldca[v0->GetIndex(0)],(2.+pulldca[v0->GetIndex(1)]) ); //pull
5530 Double_t minpulldca = TMath::Min(2.+pulldca[index0],(2.+pulldca[index1]) );//Bo:
5534 AliTPCseed * track0 = (AliTPCseed*)array->At(index0);
5535 AliTPCseed * track1 = (AliTPCseed*)array->At(index1);
5536 if (track0->TPCrPID(0)>0.3&&track1->TPCrPID(0)>0.3&&v0->GetAnglep()[2]<0.15) quality[i]+=1000000; // gamma conversion candidate
5537 if (track0->TPCrPID(4)>0.9||track1->TPCrPID(4)>0.9&&minpulldca>4) quality[i]*=10; // lambda candidate candidate
5540 TMath::Sort(ncandidates,quality,indexes,kTRUE);
5543 for (Int_t i=0;i<ncandidates;i++){
5544 AliESDv0 * v0 = (AliESDv0*)tpcv0s->At(indexes[i]);
5546 Int_t index0 = v0->GetIndex(0);
5547 Int_t index1 = v0->GetIndex(1);
5548 AliTPCseed * track0 = (AliTPCseed*)array->At(index0);
5549 AliTPCseed * track1 = (AliTPCseed*)array->At(index1);
5550 if (!track0||!track1) {
5554 Bool_t accept =kTRUE; //default accept
5555 Int_t *v0indexes0 = track0->GetV0Indexes();
5556 Int_t *v0indexes1 = track1->GetV0Indexes();
5558 Int_t order0 = (v0indexes0[0]!=0) ? 1:0;
5559 Int_t order1 = (v0indexes1[0]!=0) ? 1:0;
5560 if (v0indexes0[1]!=0) order0 =2;
5561 if (v0indexes1[1]!=0) order1 =2;
5563 if (v0indexes0[2]!=0) {order0=3; accept=kFALSE;}
5564 if (v0indexes0[2]!=0) {order1=3; accept=kFALSE;}
5566 AliESDv0 * v02 = v0;
5568 //Bo: v0->SetOrder(0,order0);
5569 //Bo: v0->SetOrder(1,order1);
5570 //Bo: v0->SetOrder(1,order0+order1);
5571 v0->SetOnFlyStatus(kTRUE);
5572 Int_t index = esd->AddV0(v0);
5573 v02 = esd->GetV0(index);
5574 v0indexes0[order0]=index;
5575 v0indexes1[order1]=index;
5579 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
5580 if (AliTPCReconstructor::StreamLevel()>1) {
5581 Int_t lab0=track0->GetLabel();
5582 Int_t lab1=track1->GetLabel();
5591 "dca0="<<dca[index0]<<
5592 "dca1="<<dca[index1]<<
5596 "quality="<<quality[i]<<
5597 "pulldca0="<<pulldca[index0]<<
5598 "pulldca1="<<pulldca[index1]<<
5622 printf("TPC V0 finder : naccepted\t%d\tncandidates\t%d\tntracks\t%d\tnall\t%d\n",naccepted,ncandidates,ntracks,nall);
5626 Int_t AliTPCtrackerMI::RefitKink(AliTPCseed &mother, AliTPCseed &daughter, AliESDkink &knk)
5629 // refit kink towards to the vertex
5632 AliKink &kink=(AliKink &)knk;
5634 Int_t row0 = GetRowNumber(kink.GetR());
5635 FollowProlongation(mother,0);
5636 mother.Reset(kFALSE);
5638 FollowProlongation(daughter,row0);
5639 daughter.Reset(kFALSE);
5640 FollowBackProlongation(daughter,158);
5641 daughter.Reset(kFALSE);
5642 Int_t first = TMath::Max(row0-20,30);
5643 Int_t last = TMath::Min(row0+20,140);
5645 const Int_t kNdiv =5;
5646 AliTPCseed param0[kNdiv]; // parameters along the track
5647 AliTPCseed param1[kNdiv]; // parameters along the track
5648 AliKink kinks[kNdiv]; // corresponding kink parameters
5651 for (Int_t irow=0; irow<kNdiv;irow++){
5652 rows[irow] = first +((last-first)*irow)/(kNdiv-1);
5654 // store parameters along the track
5656 for (Int_t irow=0;irow<kNdiv;irow++){
5657 FollowBackProlongation(mother, rows[irow]);
5658 FollowProlongation(daughter,rows[kNdiv-1-irow]);
5659 param0[irow] = mother;
5660 param1[kNdiv-1-irow] = daughter;
5664 for (Int_t irow=0; irow<kNdiv-1;irow++){
5665 if (param0[irow].GetNumberOfClusters()<kNdiv||param1[irow].GetNumberOfClusters()<kNdiv) continue;
5666 kinks[irow].SetMother(param0[irow]);
5667 kinks[irow].SetDaughter(param1[irow]);
5668 kinks[irow].Update();
5671 // choose kink with best "quality"
5673 Double_t mindist = 10000;
5674 for (Int_t irow=0;irow<kNdiv;irow++){
5675 if (param0[irow].GetNumberOfClusters()<20||param1[irow].GetNumberOfClusters()<20) continue;
5676 if (TMath::Abs(kinks[irow].GetR())>240.) continue;
5677 if (TMath::Abs(kinks[irow].GetR())<100.) continue;
5679 Float_t normdist = TMath::Abs(param0[irow].GetX()-kinks[irow].GetR())*(0.1+kink.GetDistance());
5680 normdist/= (param0[irow].GetNumberOfClusters()+param1[irow].GetNumberOfClusters()+40.);
5681 if (normdist < mindist){
5687 if (index==-1) return 0;
5690 param0[index].Reset(kTRUE);
5691 FollowProlongation(param0[index],0);
5693 mother = param0[index];
5694 daughter = param1[index]; // daughter in vertex
5696 kink.SetMother(mother);
5697 kink.SetDaughter(daughter);
5699 kink.SetTPCRow0(GetRowNumber(kink.GetR()));
5700 kink.SetTPCncls(param0[index].GetNumberOfClusters(),0);
5701 kink.SetTPCncls(param1[index].GetNumberOfClusters(),1);
5702 kink.SetLabel(CookLabel(&mother,0.4, 0,kink.GetTPCRow0()),0);
5703 kink.SetLabel(CookLabel(&daughter,0.4, kink.GetTPCRow0(),160),1);
5704 mother.SetLabel(kink.GetLabel(0));
5705 daughter.SetLabel(kink.GetLabel(1));
5711 void AliTPCtrackerMI::UpdateKinkQualityM(AliTPCseed * seed){
5713 // update Kink quality information for mother after back propagation
5715 if (seed->GetKinkIndex(0)>=0) return;
5716 for (Int_t ikink=0;ikink<3;ikink++){
5717 Int_t index = seed->GetKinkIndex(ikink);
5718 if (index>=0) break;
5719 index = TMath::Abs(index)-1;
5720 AliESDkink * kink = fEvent->GetKink(index);
5721 kink->SetTPCDensity(-1,0,0);
5722 kink->SetTPCDensity(1,0,1);
5724 Int_t row0 = kink->GetTPCRow0() - 2 - Int_t( 0.5/ (0.05+kink->GetAngle(2)));
5725 if (row0<15) row0=15;
5727 Int_t row1 = kink->GetTPCRow0() + 2 + Int_t( 0.5/ (0.05+kink->GetAngle(2)));
5728 if (row1>145) row1=145;
5730 Int_t found,foundable,shared;
5731 seed->GetClusterStatistic(0,row0, found, foundable,shared,kFALSE);
5732 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),0,0);
5733 seed->GetClusterStatistic(row1,155, found, foundable,shared,kFALSE);
5734 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),0,1);
5739 void AliTPCtrackerMI::UpdateKinkQualityD(AliTPCseed * seed){
5741 // update Kink quality information for daughter after refit
5743 if (seed->GetKinkIndex(0)<=0) return;
5744 for (Int_t ikink=0;ikink<3;ikink++){
5745 Int_t index = seed->GetKinkIndex(ikink);
5746 if (index<=0) break;
5747 index = TMath::Abs(index)-1;
5748 AliESDkink * kink = fEvent->GetKink(index);
5749 kink->SetTPCDensity(-1,1,0);
5750 kink->SetTPCDensity(-1,1,1);
5752 Int_t row0 = kink->GetTPCRow0() -2 - Int_t( 0.5/ (0.05+kink->GetAngle(2)));
5753 if (row0<15) row0=15;
5755 Int_t row1 = kink->GetTPCRow0() +2 + Int_t( 0.5/ (0.05+kink->GetAngle(2)));
5756 if (row1>145) row1=145;
5758 Int_t found,foundable,shared;
5759 seed->GetClusterStatistic(0,row0, found, foundable,shared,kFALSE);
5760 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),1,0);
5761 seed->GetClusterStatistic(row1,155, found, foundable,shared,kFALSE);
5762 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),1,1);
5768 Int_t AliTPCtrackerMI::CheckKinkPoint(AliTPCseed*seed,AliTPCseed &mother, AliTPCseed &daughter, AliESDkink &knk)
5771 // check kink point for given track
5772 // if return value=0 kink point not found
5773 // otherwise seed0 correspond to mother particle
5774 // seed1 correspond to daughter particle
5775 // kink parameter of kink point
5776 AliKink &kink=(AliKink &)knk;
5778 Int_t middlerow = (seed->GetFirstPoint()+seed->GetLastPoint())/2;
5779 Int_t first = seed->GetFirstPoint();
5780 Int_t last = seed->GetLastPoint();
5781 if (last-first<20) return 0; // shortest length - 2*30 = 60 pad-rows
5784 AliTPCseed *seed1 = ReSeed(seed,middlerow+20, kTRUE); //middle of chamber
5785 if (!seed1) return 0;
5786 FollowProlongation(*seed1,seed->GetLastPoint()-20);
5787 seed1->Reset(kTRUE);
5788 FollowProlongation(*seed1,158);
5789 seed1->Reset(kTRUE);
5790 last = seed1->GetLastPoint();
5792 AliTPCseed *seed0 = new AliTPCseed(*seed);
5793 seed0->Reset(kFALSE);
5796 AliTPCseed param0[20]; // parameters along the track
5797 AliTPCseed param1[20]; // parameters along the track
5798 AliKink kinks[20]; // corresponding kink parameters
5800 for (Int_t irow=0; irow<20;irow++){
5801 rows[irow] = first +((last-first)*irow)/19;
5803 // store parameters along the track
5805 for (Int_t irow=0;irow<20;irow++){
5806 FollowBackProlongation(*seed0, rows[irow]);
5807 FollowProlongation(*seed1,rows[19-irow]);
5808 param0[irow] = *seed0;
5809 param1[19-irow] = *seed1;
5813 for (Int_t irow=0; irow<19;irow++){
5814 kinks[irow].SetMother(param0[irow]);
5815 kinks[irow].SetDaughter(param1[irow]);
5816 kinks[irow].Update();
5819 // choose kink with biggest change of angle
5821 Double_t maxchange= 0;
5822 for (Int_t irow=1;irow<19;irow++){
5823 if (TMath::Abs(kinks[irow].GetR())>240.) continue;
5824 if (TMath::Abs(kinks[irow].GetR())<110.) continue;
5825 Float_t quality = TMath::Abs(kinks[irow].GetAngle(2))/(3.+TMath::Abs(kinks[irow].GetR()-param0[irow].GetX()));
5826 if ( quality > maxchange){
5827 maxchange = quality;
5834 if (index<0) return 0;
5836 Int_t row0 = GetRowNumber(kinks[index].GetR()); //row 0 estimate
5837 seed0 = new AliTPCseed(param0[index]);
5838 seed1 = new AliTPCseed(param1[index]);
5839 seed0->Reset(kFALSE);
5840 seed1->Reset(kFALSE);
5841 seed0->ResetCovariance(10.);
5842 seed1->ResetCovariance(10.);
5843 FollowProlongation(*seed0,0);
5844 FollowBackProlongation(*seed1,158);
5845 mother = *seed0; // backup mother at position 0
5846 seed0->Reset(kFALSE);
5847 seed1->Reset(kFALSE);
5848 seed0->ResetCovariance(10.);
5849 seed1->ResetCovariance(10.);
5851 first = TMath::Max(row0-20,0);
5852 last = TMath::Min(row0+20,158);
5854 for (Int_t irow=0; irow<20;irow++){
5855 rows[irow] = first +((last-first)*irow)/19;
5857 // store parameters along the track
5859 for (Int_t irow=0;irow<20;irow++){
5860 FollowBackProlongation(*seed0, rows[irow]);
5861 FollowProlongation(*seed1,rows[19-irow]);
5862 param0[irow] = *seed0;
5863 param1[19-irow] = *seed1;
5867 for (Int_t irow=0; irow<19;irow++){
5868 kinks[irow].SetMother(param0[irow]);
5869 kinks[irow].SetDaughter(param1[irow]);
5870 // param0[irow].Dump();
5871 //param1[irow].Dump();
5872 kinks[irow].Update();
5875 // choose kink with biggest change of angle
5878 for (Int_t irow=0;irow<20;irow++){
5879 if (TMath::Abs(kinks[irow].GetR())>250.) continue;
5880 if (TMath::Abs(kinks[irow].GetR())<90.) continue;
5881 Float_t quality = TMath::Abs(kinks[irow].GetAngle(2))/(3.+TMath::Abs(kinks[irow].GetR()-param0[irow].GetX()));
5882 if ( quality > maxchange){
5883 maxchange = quality;
5890 if (index==-1 || param0[index].GetNumberOfClusters()+param1[index].GetNumberOfClusters()<100){
5895 // Float_t anglesigma = TMath::Sqrt(param0[index].fC22+param0[index].fC33+param1[index].fC22+param1[index].fC33);
5897 kink.SetMother(param0[index]);
5898 kink.SetDaughter(param1[index]);
5900 row0 = GetRowNumber(kink.GetR());
5901 kink.SetTPCRow0(row0);
5902 kink.SetLabel(CookLabel(seed0,0.5,0,row0),0);
5903 kink.SetLabel(CookLabel(seed1,0.5,row0,158),1);
5904 kink.SetIndex(-10,0);
5905 kink.SetIndex(int(param0[index].GetNumberOfClusters()+param1[index].GetNumberOfClusters()),1);
5906 kink.SetTPCncls(param0[index].GetNumberOfClusters(),0);
5907 kink.SetTPCncls(param1[index].GetNumberOfClusters(),1);
5910 // new (&mother) AliTPCseed(param0[index]);
5911 daughter = param1[index];
5912 daughter.SetLabel(kink.GetLabel(1));
5913 param0[index].Reset(kTRUE);
5914 FollowProlongation(param0[index],0);
5915 mother = param0[index];
5916 mother.SetLabel(kink.GetLabel(0));
5926 AliTPCseed* AliTPCtrackerMI::ReSeed(AliTPCseed *t)
5929 // reseed - refit - track
5932 // Int_t last = fSectors->GetNRows()-1;
5934 if (fSectors == fOuterSec){
5935 first = TMath::Max(first, t->GetFirstPoint()-fInnerSec->GetNRows());
5939 first = t->GetFirstPoint();
5941 AliTPCseed * seed = MakeSeed(t,0.1,0.5,0.9);
5942 FollowBackProlongation(*t,fSectors->GetNRows()-1);
5944 FollowProlongation(*t,first);
5954 //_____________________________________________________________________________
5955 Int_t AliTPCtrackerMI::ReadSeeds(const TFile *inp) {
5956 //-----------------------------------------------------------------
5957 // This function reades track seeds.
5958 //-----------------------------------------------------------------
5959 TDirectory *savedir=gDirectory;
5961 TFile *in=(TFile*)inp;
5962 if (!in->IsOpen()) {
5963 cerr<<"AliTPCtrackerMI::ReadSeeds(): input file is not open !\n";
5968 TTree *seedTree=(TTree*)in->Get("Seeds");
5970 cerr<<"AliTPCtrackerMI::ReadSeeds(): ";
5971 cerr<<"can't get a tree with track seeds !\n";
5974 AliTPCtrack *seed=new AliTPCtrack;
5975 seedTree->SetBranchAddress("tracks",&seed);
5977 if (fSeeds==0) fSeeds=new TObjArray(15000);
5979 Int_t n=(Int_t)seedTree->GetEntries();
5980 for (Int_t i=0; i<n; i++) {
5981 seedTree->GetEvent(i);
5982 fSeeds->AddLast(new AliTPCseed(*seed/*,seed->GetAlpha()*/));
5991 Int_t AliTPCtrackerMI::Clusters2Tracks (AliESDEvent *esd)
5994 if (fSeeds) DeleteSeeds();
5997 if (!fSeeds) return 1;
6004 //_____________________________________________________________________________
6005 Int_t AliTPCtrackerMI::Clusters2Tracks() {
6006 //-----------------------------------------------------------------
6007 // This is a track finder.
6008 //-----------------------------------------------------------------
6009 TDirectory *savedir=gDirectory;
6013 fSeeds = Tracking();
6016 Info("Clusters2Tracks","Time for tracking: \t");timer.Print();timer.Start();
6018 //activate again some tracks
6019 for (Int_t i=0; i<fSeeds->GetEntriesFast(); i++) {
6020 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
6022 Int_t nc=t.GetNumberOfClusters();
6024 delete fSeeds->RemoveAt(i);
6028 if (pt->GetRemoval()==10) {
6029 if (pt->GetDensityFirst(20)>0.8 || pt->GetDensityFirst(30)>0.8 || pt->GetDensityFirst(40)>0.7)
6030 pt->Desactivate(10); // make track again active
6032 pt->Desactivate(20);
6033 delete fSeeds->RemoveAt(i);
6038 RemoveUsed2(fSeeds,0.85,0.85,0);
6039 if (AliTPCReconstructor::GetRecoParam()->GetDoKinks()) FindKinks(fSeeds,fEvent);
6040 //FindCurling(fSeeds, fEvent,0);
6041 if (AliTPCReconstructor::StreamLevel()>2) FindMultiMC(fSeeds, fEvent,0); // find multi found tracks
6042 RemoveUsed2(fSeeds,0.5,0.4,20);
6043 FindSplitted(fSeeds, fEvent,0); // find multi found tracks
6046 // // refit short tracks
6048 Int_t nseed=fSeeds->GetEntriesFast();
6051 for (Int_t i=0; i<nseed; i++) {
6052 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
6054 Int_t nc=t.GetNumberOfClusters();
6056 delete fSeeds->RemoveAt(i);
6059 CookLabel(pt,0.1); //For comparison only
6060 //if ((pt->IsActive() || (pt->fRemoval==10) )&& nc>50 &&pt->GetNumberOfClusters()>0.4*pt->fNFoundable){
6061 if ((pt->IsActive() || (pt->GetRemoval()==10) )){
6063 if (fDebug>0) cerr<<found<<'\r';
6067 delete fSeeds->RemoveAt(i);
6071 //RemoveOverlap(fSeeds,0.99,7,kTRUE);
6073 //RemoveUsed(fSeeds,0.9,0.9,6);
6075 nseed=fSeeds->GetEntriesFast();
6077 for (Int_t i=0; i<nseed; i++) {
6078 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
6080 Int_t nc=t.GetNumberOfClusters();
6082 delete fSeeds->RemoveAt(i);
6086 t.CookdEdx(0.02,0.6);
6087 // CheckKinkPoint(&t,0.05);
6088 //if ((pt->IsActive() || (pt->fRemoval==10) )&& nc>50 &&pt->GetNumberOfClusters()>0.4*pt->fNFoundable){
6089 if ((pt->IsActive() || (pt->GetRemoval()==10) )){
6097 delete fSeeds->RemoveAt(i);
6098 //AliTPCseed * seed1 = ReSeed(pt,0.05,0.5,1);
6100 // FollowProlongation(*seed1,0);
6101 // Int_t n = seed1->GetNumberOfClusters();
6102 // printf("fP4\t%f\t%f\n",seed1->GetC(),pt->GetC());
6103 // printf("fN\t%d\t%d\n", seed1->GetNumberOfClusters(),pt->GetNumberOfClusters());
6106 //AliTPCseed * seed2 = ReSeed(pt,0.95,0.5,0.05);
6110 SortTracks(fSeeds, 1);
6114 PrepareForBackProlongation(fSeeds,5.);
6115 PropagateBack(fSeeds);
6116 printf("Time for back propagation: \t");timer.Print();timer.Start();
6120 PrepareForProlongation(fSeeds,5.);
6121 PropagateForward2(fSeeds);
6123 printf("Time for FORWARD propagation: \t");timer.Print();timer.Start();
6124 // RemoveUsed(fSeeds,0.7,0.7,6);
6125 //RemoveOverlap(fSeeds,0.9,7,kTRUE);
6127 nseed=fSeeds->GetEntriesFast();
6129 for (Int_t i=0; i<nseed; i++) {
6130 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
6132 Int_t nc=t.GetNumberOfClusters();
6134 delete fSeeds->RemoveAt(i);
6137 t.CookdEdx(0.02,0.6);
6138 // CookLabel(pt,0.1); //For comparison only
6139 //if ((pt->IsActive() || (pt->fRemoval==10) )&& nc>50 &&pt->GetNumberOfClusters()>0.4*pt->fNFoundable){
6140 if ((pt->IsActive() || (pt->fRemoval==10) )){
6141 cerr<<found++<<'\r';
6144 delete fSeeds->RemoveAt(i);
6149 // fNTracks = found;
6151 Info("Clusters2Tracks","Time for overlap removal, track writing and dedx cooking: \t"); timer.Print();timer.Start();
6154 // cerr<<"Number of found tracks : "<<"\t"<<found<<endl;
6155 Info("Clusters2Tracks","Number of found tracks %d",found);
6157 // UnloadClusters();
6162 void AliTPCtrackerMI::Tracking(TObjArray * arr)
6165 // tracking of the seeds
6168 fSectors = fOuterSec;
6169 ParallelTracking(arr,150,63);
6170 fSectors = fOuterSec;
6171 ParallelTracking(arr,63,0);
6174 TObjArray * AliTPCtrackerMI::Tracking(Int_t seedtype, Int_t i1, Int_t i2, Float_t cuts[4], Float_t dy, Int_t dsec)
6179 TObjArray * arr = new TObjArray;
6181 fSectors = fOuterSec;
6184 for (Int_t sec=0;sec<fkNOS;sec++){
6185 if (seedtype==3) MakeSeeds3(arr,sec,i1,i2,cuts,dy, dsec);
6186 if (seedtype==4) MakeSeeds5(arr,sec,i1,i2,cuts,dy);
6187 if (seedtype==2) MakeSeeds2(arr,sec,i1,i2,cuts,dy);
6190 Info("Tracking","\nSeeding - %d\t%d\t%d\t%d\n",seedtype,i1,i2,arr->GetEntriesFast());
6202 TObjArray * AliTPCtrackerMI::Tracking()
6206 if (AliTPCReconstructor::GetRecoParam()->GetSpecialSeeding()) return TrackingSpecial();
6209 Int_t nup=fOuterSec->GetNRows()+fInnerSec->GetNRows();
6211 TObjArray * seeds = new TObjArray;
6220 Float_t fnumber = 3.0;
6221 Float_t fdensity = 3.0;
6226 for (Int_t delta = 0; delta<18; delta+=6){
6230 arr = Tracking(3,nup-1-delta,nup-1-delta-gap,cuts,-1,1);
6231 SumTracks(seeds,arr);
6232 SignClusters(seeds,fnumber,fdensity);
6234 for (Int_t i=2;i<6;i+=2){
6235 // seed high pt tracks
6238 arr = Tracking(3,nup-i-delta,nup-i-delta-gap,cuts,-1,0);
6239 SumTracks(seeds,arr);
6240 SignClusters(seeds,fnumber,fdensity);
6245 // RemoveUsed(seeds,0.9,0.9,1);
6246 // UnsignClusters();
6247 // SignClusters(seeds,fnumber,fdensity);
6251 for (Int_t delta = 20; delta<120; delta+=10){
6253 // seed high pt tracks
6257 arr = Tracking(3,nup-delta,nup-delta-gap,cuts,-1);
6258 SumTracks(seeds,arr);
6259 SignClusters(seeds,fnumber,fdensity);
6264 arr = Tracking(3,nup-delta-5,nup-delta-5-gap,cuts,-1);
6265 SumTracks(seeds,arr);
6266 SignClusters(seeds,fnumber,fdensity);
6277 Info("Tracking()","\n\nPrimary seeding\t%d\n\n",seeds->GetEntriesFast());
6281 // RemoveUsed(seeds,0.75,0.75,1);
6283 //SignClusters(seeds,fnumber,fdensity);
6292 arr = Tracking(4,nup-1,nup-1-gap,cuts,-1);
6293 SumTracks(seeds,arr);
6294 SignClusters(seeds,fnumber,fdensity);
6296 arr = Tracking(4,nup-2,nup-2-gap,cuts,-1);
6297 SumTracks(seeds,arr);
6298 SignClusters(seeds,fnumber,fdensity);
6300 arr = Tracking(4,nup-3,nup-3-gap,cuts,-1);
6301 SumTracks(seeds,arr);
6302 SignClusters(seeds,fnumber,fdensity);
6306 for (Int_t delta = 3; delta<30; delta+=5){
6312 arr = Tracking(4,nup-1-delta,nup-1-delta-gap,cuts,-1);
6313 SumTracks(seeds,arr);
6314 SignClusters(seeds,fnumber,fdensity);
6316 arr = Tracking(4,nup-3-delta,nup-5-delta-gap,cuts,4);
6317 SumTracks(seeds,arr);
6318 SignClusters(seeds,fnumber,fdensity);
6329 Int_t fLastSeedRowSec=AliTPCReconstructor::GetRecoParam()->GetLastSeedRowSec();
6332 for (Int_t delta = 30; delta<fLastSeedRowSec; delta+=10){
6338 arr = Tracking(4,nup-1-delta,nup-1-delta-gap,cuts,-1);
6339 SumTracks(seeds,arr);
6340 SignClusters(seeds,fnumber,fdensity);
6342 arr = Tracking(4,nup-5-delta,nup-5-delta-gap,cuts,5 );
6343 SumTracks(seeds,arr);
6344 SignClusters(seeds,fnumber,fdensity);
6348 Info("Tracking()","\n\nSecondary seeding\t%d\n\n",seeds->GetEntriesFast());
6359 TObjArray * AliTPCtrackerMI::TrackingSpecial()
6362 // seeding adjusted for laser and cosmic tests - short tracks with big inclination angle
6363 // no primary vertex seeding tried
6367 Int_t nup=fOuterSec->GetNRows()+fInnerSec->GetNRows();
6369 TObjArray * seeds = new TObjArray;
6374 Float_t fnumber = 3.0;
6375 Float_t fdensity = 3.0;
6378 cuts[0] = AliTPCReconstructor::GetRecoParam()->GetMaxC(); // max curvature
6379 cuts[1] = 3.5; // max tan(phi) angle for seeding
6380 cuts[2] = 3.; // not used (cut on z primary vertex)
6381 cuts[3] = 3.5; // max tan(theta) angle for seeding
6383 for (Int_t delta = 0; nup-delta-gap-1>0; delta+=3){
6385 arr = Tracking(4,nup-1-delta,nup-1-delta-gap,cuts,-1);
6386 SumTracks(seeds,arr);
6387 SignClusters(seeds,fnumber,fdensity);
6391 Info("Tracking()","\n\nSecondary seeding\t%d\n\n",seeds->GetEntriesFast());
6402 void AliTPCtrackerMI::SumTracks(TObjArray *arr1,TObjArray *&arr2) const
6405 //sum tracks to common container
6406 //remove suspicious tracks
6407 Int_t nseed = arr2->GetEntriesFast();
6408 for (Int_t i=0;i<nseed;i++){
6409 AliTPCseed *pt=(AliTPCseed*)arr2->UncheckedAt(i);
6412 // remove tracks with too big curvature
6414 if (TMath::Abs(pt->GetC())>AliTPCReconstructor::GetRecoParam()->GetMaxC()){
6415 delete arr2->RemoveAt(i);
6418 // REMOVE VERY SHORT TRACKS
6419 if (pt->GetNumberOfClusters()<20){
6420 delete arr2->RemoveAt(i);
6423 // NORMAL ACTIVE TRACK
6424 if (pt->IsActive()){
6425 arr1->AddLast(arr2->RemoveAt(i));
6428 //remove not usable tracks
6429 if (pt->GetRemoval()!=10){
6430 delete arr2->RemoveAt(i);
6434 // ENABLE ONLY ENOUGH GOOD STOPPED TRACKS
6435 if (pt->GetDensityFirst(20)>0.8 || pt->GetDensityFirst(30)>0.8 || pt->GetDensityFirst(40)>0.7)
6436 arr1->AddLast(arr2->RemoveAt(i));
6438 delete arr2->RemoveAt(i);
6442 delete arr2; arr2 = 0;
6447 void AliTPCtrackerMI::ParallelTracking(TObjArray * arr, Int_t rfirst, Int_t rlast)
6450 // try to track in parralel
6452 Int_t nseed=arr->GetEntriesFast();
6453 //prepare seeds for tracking
6454 for (Int_t i=0; i<nseed; i++) {
6455 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i), &t=*pt;
6457 if (!t.IsActive()) continue;
6458 // follow prolongation to the first layer
6459 if ( (fSectors ==fInnerSec) || (t.GetFirstPoint()-fParam->GetNRowLow()>rfirst+1) )
6460 FollowProlongation(t, rfirst+1);
6465 for (Int_t nr=rfirst; nr>=rlast; nr--){
6466 if (nr<fInnerSec->GetNRows())
6467 fSectors = fInnerSec;
6469 fSectors = fOuterSec;
6470 // make indexes with the cluster tracks for given
6472 // find nearest cluster
6473 for (Int_t i=0; i<nseed; i++) {
6474 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i), &t=*pt;
6476 if (nr==80) pt->UpdateReference();
6477 if (!pt->IsActive()) continue;
6478 // if ( (fSectors ==fOuterSec) && (pt->fFirstPoint-fParam->GetNRowLow())<nr) continue;
6479 if (pt->GetRelativeSector()>17) {
6482 UpdateClusters(t,nr);
6484 // prolonagate to the nearest cluster - if founded
6485 for (Int_t i=0; i<nseed; i++) {
6486 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
6488 if (!pt->IsActive()) continue;
6489 // if ((fSectors ==fOuterSec) && (pt->fFirstPoint-fParam->GetNRowLow())<nr) continue;
6490 if (pt->GetRelativeSector()>17) {
6493 FollowToNextCluster(*pt,nr);
6498 void AliTPCtrackerMI::PrepareForBackProlongation(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);
6510 //rotate to current local system at first accepted point
6511 Int_t index = pt->GetClusterIndex2(pt->GetFirstPoint());
6512 Int_t sec = (index&0xff000000)>>24;
6514 Float_t angle1 = fInnerSec->GetAlpha()*sec+fInnerSec->GetAlphaShift();
6515 if (angle1>TMath::Pi())
6516 angle1-=2.*TMath::Pi();
6517 Float_t angle2 = pt->GetAlpha();
6519 if (TMath::Abs(angle1-angle2)>0.001){
6520 pt->Rotate(angle1-angle2);
6521 //angle2 = pt->GetAlpha();
6522 //pt->fRelativeSector = pt->GetAlpha()/fInnerSec->GetAlpha();
6523 //if (pt->GetAlpha()<0)
6524 // pt->fRelativeSector+=18;
6525 //sec = pt->fRelativeSector;
6534 void AliTPCtrackerMI::PrepareForProlongation(TObjArray * arr, Float_t fac) const
6538 // if we use TPC track itself we have to "update" covariance
6540 Int_t nseed= arr->GetEntriesFast();
6541 for (Int_t i=0;i<nseed;i++){
6542 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
6545 pt->SetFirstPoint(pt->GetLastPoint());
6553 Int_t AliTPCtrackerMI::PropagateBack(TObjArray * arr)
6556 // make back propagation
6558 Int_t nseed= arr->GetEntriesFast();
6559 for (Int_t i=0;i<nseed;i++){
6560 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
6561 if (pt&& pt->GetKinkIndex(0)<=0) {
6562 //AliTPCseed *pt2 = new AliTPCseed(*pt);
6563 fSectors = fInnerSec;
6564 //FollowBackProlongation(*pt,fInnerSec->GetNRows()-1);
6565 //fSectors = fOuterSec;
6566 FollowBackProlongation(*pt,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1);
6567 //if (pt->GetNumberOfClusters()<(pt->fEsd->GetTPCclusters(0)) ){
6568 // Error("PropagateBack","Not prolonged track %d",pt->GetLabel());
6569 // FollowBackProlongation(*pt2,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1);
6572 if (pt&& pt->GetKinkIndex(0)>0) {
6573 AliESDkink * kink = fEvent->GetKink(pt->GetKinkIndex(0)-1);
6574 pt->SetFirstPoint(kink->GetTPCRow0());
6575 fSectors = fInnerSec;
6576 FollowBackProlongation(*pt,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1);
6584 Int_t AliTPCtrackerMI::PropagateForward2(TObjArray * arr)
6587 // make forward propagation
6589 Int_t nseed= arr->GetEntriesFast();
6591 for (Int_t i=0;i<nseed;i++){
6592 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
6594 FollowProlongation(*pt,0);
6603 Int_t AliTPCtrackerMI::PropagateForward()
6606 // propagate track forward
6608 Int_t nseed = fSeeds->GetEntriesFast();
6609 for (Int_t i=0;i<nseed;i++){
6610 AliTPCseed *pt = (AliTPCseed*)fSeeds->UncheckedAt(i);
6612 AliTPCseed &t = *pt;
6613 Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
6614 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
6615 if (alpha < 0. ) alpha += 2.*TMath::Pi();
6616 t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN);
6620 fSectors = fOuterSec;
6621 ParallelTracking(fSeeds,fOuterSec->GetNRows()+fInnerSec->GetNRows()-1,fInnerSec->GetNRows());
6622 fSectors = fInnerSec;
6623 ParallelTracking(fSeeds,fInnerSec->GetNRows()-1,0);
6632 Int_t AliTPCtrackerMI::PropagateBack(AliTPCseed * pt, Int_t row0, Int_t row1)
6635 // make back propagation, in between row0 and row1
6639 fSectors = fInnerSec;
6642 if (row1<fSectors->GetNRows())
6645 r1 = fSectors->GetNRows()-1;
6647 if (row0<fSectors->GetNRows()&& r1>0 )
6648 FollowBackProlongation(*pt,r1);
6649 if (row1<=fSectors->GetNRows())
6652 r1 = row1 - fSectors->GetNRows();
6653 if (r1<=0) return 0;
6654 if (r1>=fOuterSec->GetNRows()) return 0;
6655 fSectors = fOuterSec;
6656 return FollowBackProlongation(*pt,r1);
6664 void AliTPCtrackerMI::GetShape(AliTPCseed * seed, Int_t row)
6668 AliTPCClusterParam * clparam = AliTPCcalibDB::Instance()->GetClusterParam();
6669 Float_t zdrift = TMath::Abs((fParam->GetZLength(0)-TMath::Abs(seed->GetZ())));
6670 Int_t type = (seed->GetSector() < fParam->GetNSector()/2) ? 0: (row>126) ? 1:2;
6671 Double_t angulary = seed->GetSnp();
6672 angulary = angulary*angulary/(1.-angulary*angulary);
6673 Double_t angularz = seed->GetTgl()*seed->GetTgl()*(1.+angulary);
6675 Double_t sigmay = clparam->GetRMS0(0,type,zdrift,TMath::Sqrt(TMath::Abs(angulary)));
6676 Double_t sigmaz = clparam->GetRMS0(1,type,zdrift,TMath::Sqrt(TMath::Abs(angularz)));
6677 seed->SetCurrentSigmaY2(sigmay*sigmay);
6678 seed->SetCurrentSigmaZ2(sigmaz*sigmaz);
6679 // Float_t sd2 = TMath::Abs((fParam->GetZLength(0)-TMath::Abs(seed->GetZ())))*fParam->GetDiffL()*fParam->GetDiffL();
6680 // // Float_t padlength = fParam->GetPadPitchLength(seed->fSector);
6681 // Float_t padlength = GetPadPitchLength(row);
6683 // Float_t sresy = (seed->GetSector() < fParam->GetNSector()/2) ? 0.2 :0.3;
6684 // seed->SetCurrentSigmaY2(sd2+padlength*padlength*angulary/12.+sresy*sresy);
6686 // Float_t sresz = fParam->GetZSigma();
6687 // seed->SetCurrentSigmaZ2(sd2+padlength*padlength*angularz*angularz*(1+angulary)/12.+sresz*sresz);
6689 Float_t wy = GetSigmaY(seed);
6690 Float_t wz = GetSigmaZ(seed);
6693 if (TMath::Abs(wy/seed->fCurrentSigmaY2-1)>0.0001 || TMath::Abs(wz/seed->fCurrentSigmaZ2-1)>0.0001 ){
6694 printf("problem\n");
6701 //__________________________________________________________________________
6702 void AliTPCtrackerMI::CookLabel(AliKalmanTrack *tk, Float_t wrong) const {
6703 //--------------------------------------------------------------------
6704 //This function "cooks" a track label. If label<0, this track is fake.
6705 //--------------------------------------------------------------------
6706 AliTPCseed * t = dynamic_cast<AliTPCseed*>(tk);
6708 printf("%s:%d wrong type \n",(char*)__FILE__,__LINE__);
6712 Int_t noc=t->GetNumberOfClusters();
6714 //printf("\nnot founded prolongation\n\n\n");
6720 AliTPCclusterMI *clusters[160];
6722 for (Int_t i=0;i<160;i++) {
6729 for (i=0; i<160 && current<noc; i++) {
6731 Int_t index=t->GetClusterIndex2(i);
6732 if (index<=0) continue;
6733 if (index&0x8000) continue;
6735 //clusters[current]=GetClusterMI(index);
6736 if (t->GetClusterPointer(i)){
6737 clusters[current]=t->GetClusterPointer(i);
6743 Int_t lab=123456789;
6744 for (i=0; i<noc; i++) {
6745 AliTPCclusterMI *c=clusters[i];
6747 lab=TMath::Abs(c->GetLabel(0));
6749 for (j=0; j<noc; j++) if (lb[j]==lab || mx[j]==0) break;
6755 for (i=0; i<noc; i++) if (mx[i]>max) {max=mx[i]; lab=lb[i];}
6757 for (i=0; i<noc; i++) {
6758 AliTPCclusterMI *c=clusters[i];
6760 if (TMath::Abs(c->GetLabel(1)) == lab ||
6761 TMath::Abs(c->GetLabel(2)) == lab ) max++;
6764 if ((1.- Float_t(max)/noc) > wrong) lab=-lab;
6767 Int_t tail=Int_t(0.10*noc);
6770 for (i=1; i<=160&&ind<tail; i++) {
6771 // AliTPCclusterMI *c=clusters[noc-i];
6772 AliTPCclusterMI *c=clusters[i];
6774 if (lab == TMath::Abs(c->GetLabel(0)) ||
6775 lab == TMath::Abs(c->GetLabel(1)) ||
6776 lab == TMath::Abs(c->GetLabel(2))) max++;
6779 if (max < Int_t(0.5*tail)) lab=-lab;
6786 //delete[] clusters;
6790 //__________________________________________________________________________
6791 Int_t AliTPCtrackerMI::CookLabel(AliTPCseed *t, Float_t wrong,Int_t first, Int_t last) const {
6792 //--------------------------------------------------------------------
6793 //This function "cooks" a track label. If label<0, this track is fake.
6794 //--------------------------------------------------------------------
6795 Int_t noc=t->GetNumberOfClusters();
6797 //printf("\nnot founded prolongation\n\n\n");
6803 AliTPCclusterMI *clusters[160];
6805 for (Int_t i=0;i<160;i++) {
6812 for (i=0; i<160 && current<noc; i++) {
6813 if (i<first) continue;
6814 if (i>last) continue;
6815 Int_t index=t->GetClusterIndex2(i);
6816 if (index<=0) continue;
6817 if (index&0x8000) continue;
6819 //clusters[current]=GetClusterMI(index);
6820 if (t->GetClusterPointer(i)){
6821 clusters[current]=t->GetClusterPointer(i);
6826 if (noc<5) return -1;
6827 Int_t lab=123456789;
6828 for (i=0; i<noc; i++) {
6829 AliTPCclusterMI *c=clusters[i];
6831 lab=TMath::Abs(c->GetLabel(0));
6833 for (j=0; j<noc; j++) if (lb[j]==lab || mx[j]==0) break;
6839 for (i=0; i<noc; i++) if (mx[i]>max) {max=mx[i]; lab=lb[i];}
6841 for (i=0; i<noc; i++) {
6842 AliTPCclusterMI *c=clusters[i];
6844 if (TMath::Abs(c->GetLabel(1)) == lab ||
6845 TMath::Abs(c->GetLabel(2)) == lab ) max++;
6848 if ((1.- Float_t(max)/noc) > wrong) lab=-lab;
6851 Int_t tail=Int_t(0.10*noc);
6854 for (i=1; i<=160&&ind<tail; i++) {
6855 // AliTPCclusterMI *c=clusters[noc-i];
6856 AliTPCclusterMI *c=clusters[i];
6858 if (lab == TMath::Abs(c->GetLabel(0)) ||
6859 lab == TMath::Abs(c->GetLabel(1)) ||
6860 lab == TMath::Abs(c->GetLabel(2))) max++;
6863 if (max < Int_t(0.5*tail)) lab=-lab;
6866 // t->SetLabel(lab);
6870 //delete[] clusters;
6874 Int_t AliTPCtrackerMI::GetRowNumber(Double_t x[3]) const
6876 //return pad row number for given x vector
6877 Float_t phi = TMath::ATan2(x[1],x[0]);
6878 if(phi<0) phi=2.*TMath::Pi()+phi;
6879 // Get the local angle in the sector philoc
6880 const Float_t kRaddeg = 180/3.14159265358979312;
6881 Float_t phiangle = (Int_t (phi*kRaddeg/20.) + 0.5)*20./kRaddeg;
6882 Double_t localx = x[0]*TMath::Cos(phiangle)-x[1]*TMath::Sin(phiangle);
6883 return GetRowNumber(localx);
6888 void AliTPCtrackerMI::MakeBitmaps(AliTPCseed *t)
6890 //-----------------------------------------------------------------------
6891 // Fill the cluster and sharing bitmaps of the track
6892 //-----------------------------------------------------------------------
6894 Int_t firstpoint = 0;
6895 Int_t lastpoint = 159;
6896 AliTPCTrackerPoint *point;
6898 for (int iter=firstpoint; iter<lastpoint; iter++) {
6899 point = t->GetTrackPoint(iter);
6901 t->SetClusterMapBit(iter, kTRUE);
6902 if (point->IsShared())
6903 t->SetSharedMapBit(iter,kTRUE);
6905 t->SetSharedMapBit(iter, kFALSE);
6908 t->SetClusterMapBit(iter, kFALSE);
6909 t->SetSharedMapBit(iter, kFALSE);