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)
128 class AliTPCFastMath {
131 static Double_t FastAsin(Double_t x);
133 static Double_t fgFastAsin[20000]; //lookup table for fast asin computation
136 Double_t AliTPCFastMath::fgFastAsin[20000];
137 AliTPCFastMath gAliTPCFastMath; // needed to fill the LUT
139 AliTPCFastMath::AliTPCFastMath(){
141 // initialized lookup table;
142 for (Int_t i=0;i<10000;i++){
143 fgFastAsin[2*i] = TMath::ASin(i/10000.);
144 fgFastAsin[2*i+1] = (TMath::ASin((i+1)/10000.)-fgFastAsin[2*i]);
148 Double_t AliTPCFastMath::FastAsin(Double_t x){
150 // return asin using lookup table
152 Int_t index = int(x*10000);
153 return fgFastAsin[2*index]+(x*10000.-index)*fgFastAsin[2*index+1];
156 Int_t index = int(x*10000);
157 return -(fgFastAsin[2*index]+(x*10000.-index)*fgFastAsin[2*index+1]);
159 //__________________________________________________________________
160 AliTPCtrackerMI::AliTPCtrackerMI()
182 // default constructor
185 //_____________________________________________________________________
189 Int_t AliTPCtrackerMI::UpdateTrack(AliTPCseed * track, Int_t accept){
191 //update track information using current cluster - track->fCurrentCluster
194 AliTPCclusterMI* c =track->GetCurrentCluster();
195 if (accept>0) track->SetCurrentClusterIndex1(track->GetCurrentClusterIndex1() | 0x8000); //sign not accepted clusters
197 UInt_t i = track->GetCurrentClusterIndex1();
199 Int_t sec=(i&0xff000000)>>24;
200 //Int_t row = (i&0x00ff0000)>>16;
201 track->SetRow((i&0x00ff0000)>>16);
202 track->SetSector(sec);
203 // Int_t index = i&0xFFFF;
204 if (sec>=fParam->GetNInnerSector()) track->SetRow(track->GetRow()+fParam->GetNRowLow());
205 track->SetClusterIndex2(track->GetRow(), i);
206 //track->fFirstPoint = row;
207 //if ( track->fLastPoint<row) track->fLastPoint =row;
208 // if (track->fRow<0 || track->fRow>160) {
209 // printf("problem\n");
211 if (track->GetFirstPoint()>track->GetRow())
212 track->SetFirstPoint(track->GetRow());
213 if (track->GetLastPoint()<track->GetRow())
214 track->SetLastPoint(track->GetRow());
217 track->SetClusterPointer(track->GetRow(),c);
220 Double_t angle2 = track->GetSnp()*track->GetSnp();
222 //SET NEW Track Point
224 if (angle2<1) //PH sometimes angle2 is very big. To be investigated...
226 angle2 = TMath::Sqrt(angle2/(1-angle2));
227 AliTPCTrackerPoint &point =*(track->GetTrackPoint(track->GetRow()));
229 point.SetSigmaY(c->GetSigmaY2()/track->GetCurrentSigmaY2());
230 point.SetSigmaZ(c->GetSigmaZ2()/track->GetCurrentSigmaZ2());
231 point.SetErrY(sqrt(track->GetErrorY2()));
232 point.SetErrZ(sqrt(track->GetErrorZ2()));
234 point.SetX(track->GetX());
235 point.SetY(track->GetY());
236 point.SetZ(track->GetZ());
237 point.SetAngleY(angle2);
238 point.SetAngleZ(track->GetTgl());
239 if (point.IsShared()){
240 track->SetErrorY2(track->GetErrorY2()*4);
241 track->SetErrorZ2(track->GetErrorZ2()*4);
245 Double_t chi2 = track->GetPredictedChi2(track->GetCurrentCluster());
247 // track->SetErrorY2(track->GetErrorY2()*1.3);
248 // track->SetErrorY2(track->GetErrorY2()+0.01);
249 // track->SetErrorZ2(track->GetErrorZ2()*1.3);
250 // track->SetErrorZ2(track->GetErrorZ2()+0.005);
252 if (accept>0) return 0;
253 if (track->GetNumberOfClusters()%20==0){
254 // if (track->fHelixIn){
255 // TClonesArray & larr = *(track->fHelixIn);
256 // Int_t ihelix = larr.GetEntriesFast();
257 // new(larr[ihelix]) AliHelix(*track) ;
260 track->SetNoCluster(0);
261 return track->Update(c,chi2,i);
266 Int_t AliTPCtrackerMI::AcceptCluster(AliTPCseed * seed, AliTPCclusterMI * cluster)
269 // decide according desired precision to accept given
270 // cluster for tracking
271 Double_t sy2=ErrY2(seed,cluster);
272 Double_t sz2=ErrZ2(seed,cluster);
274 Double_t sdistancey2 = sy2+seed->GetSigmaY2();
275 Double_t sdistancez2 = sz2+seed->GetSigmaZ2();
277 Double_t rdistancey2 = (seed->GetCurrentCluster()->GetY()-seed->GetY())*
278 (seed->GetCurrentCluster()->GetY()-seed->GetY())/sdistancey2;
279 Double_t rdistancez2 = (seed->GetCurrentCluster()->GetZ()-seed->GetZ())*
280 (seed->GetCurrentCluster()->GetZ()-seed->GetZ())/sdistancez2;
282 Double_t rdistance2 = rdistancey2+rdistancez2;
285 if (AliTPCReconstructor::StreamLevel()>5) {
286 Float_t rmsy2 = seed->GetCurrentSigmaY2();
287 Float_t rmsz2 = seed->GetCurrentSigmaZ2();
288 Float_t rmsy2p30 = seed->GetCMeanSigmaY2p30();
289 Float_t rmsz2p30 = seed->GetCMeanSigmaZ2p30();
290 Float_t rmsy2p30R = seed->GetCMeanSigmaY2p30R();
291 Float_t rmsz2p30R = seed->GetCMeanSigmaZ2p30R();
293 AliExternalTrackParam param(*seed);
294 (*fDebugStreamer)<<"ErrParam"<<
301 "rmsy2p30="<<rmsy2p30<<
302 "rmsz2p30="<<rmsz2p30<<
303 "rmsy2p30R="<<rmsy2p30R<<
304 "rmsz2p30R="<<rmsz2p30R<<
308 if (rdistance2>16) return 3;
311 if ((rdistancey2>9. || rdistancez2>9.) && cluster->GetType()==0)
312 return 2; //suspisiouce - will be changed
314 if ((rdistancey2>6.25 || rdistancez2>6.25) && cluster->GetType()>0)
315 // strict cut on overlaped cluster
316 return 2; //suspisiouce - will be changed
318 if ( (rdistancey2>1. || rdistancez2>6.25 )
319 && cluster->GetType()<0){
320 seed->SetNFoundable(seed->GetNFoundable()-1);
329 //_____________________________________________________________________________
330 AliTPCtrackerMI::AliTPCtrackerMI(const AliTPCParam *par):
332 fkNIS(par->GetNInnerSector()/2),
334 fkNOS(par->GetNOuterSector()/2),
351 //---------------------------------------------------------------------
352 // The main TPC tracker constructor
353 //---------------------------------------------------------------------
354 fInnerSec=new AliTPCtrackerSector[fkNIS];
355 fOuterSec=new AliTPCtrackerSector[fkNOS];
358 for (i=0; i<fkNIS; i++) fInnerSec[i].Setup(par,0);
359 for (i=0; i<fkNOS; i++) fOuterSec[i].Setup(par,1);
362 Int_t nrowlow = par->GetNRowLow();
363 Int_t nrowup = par->GetNRowUp();
366 for (Int_t i=0;i<nrowlow;i++){
367 fXRow[i] = par->GetPadRowRadiiLow(i);
368 fPadLength[i]= par->GetPadPitchLength(0,i);
369 fYMax[i] = fXRow[i]*TMath::Tan(0.5*par->GetInnerAngle());
373 for (Int_t i=0;i<nrowup;i++){
374 fXRow[i+nrowlow] = par->GetPadRowRadiiUp(i);
375 fPadLength[i+nrowlow] = par->GetPadPitchLength(60,i);
376 fYMax[i+nrowlow] = fXRow[i+nrowlow]*TMath::Tan(0.5*par->GetOuterAngle());
379 fDebugStreamer = new TTreeSRedirector("TPCdebug.root");
381 //________________________________________________________________________
382 AliTPCtrackerMI::AliTPCtrackerMI(const AliTPCtrackerMI &t):
403 //------------------------------------
404 // dummy copy constructor
405 //------------------------------------------------------------------
408 AliTPCtrackerMI & AliTPCtrackerMI::operator=(const AliTPCtrackerMI& /*r*/){
409 //------------------------------
411 //--------------------------------------------------------------
414 //_____________________________________________________________________________
415 AliTPCtrackerMI::~AliTPCtrackerMI() {
416 //------------------------------------------------------------------
417 // TPC tracker destructor
418 //------------------------------------------------------------------
425 if (fDebugStreamer) delete fDebugStreamer;
429 void AliTPCtrackerMI::FillESD(TObjArray* arr)
433 //fill esds using updated tracks
435 // write tracks to the event
436 // store index of the track
437 Int_t nseed=arr->GetEntriesFast();
438 //FindKinks(arr,fEvent);
439 for (Int_t i=0; i<nseed; i++) {
440 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
443 // pt->PropagateTo(fParam->GetInnerRadiusLow());
444 if (pt->GetKinkIndex(0)<=0){ //don't propagate daughter tracks
445 pt->PropagateTo(fParam->GetInnerRadiusLow());
448 if (( pt->GetPoints()[2]- pt->GetPoints()[0])>5 && pt->GetPoints()[3]>0.8){
450 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
451 iotrack.SetTPCPoints(pt->GetPoints());
452 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
453 iotrack.SetV0Indexes(pt->GetV0Indexes());
454 // iotrack.SetTPCpid(pt->fTPCr);
455 //iotrack.SetTPCindex(i);
456 fEvent->AddTrack(&iotrack);
460 if ( (pt->GetNumberOfClusters()>70)&& (Float_t(pt->GetNumberOfClusters())/Float_t(pt->GetNFoundable()))>0.55) {
462 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
463 iotrack.SetTPCPoints(pt->GetPoints());
464 //iotrack.SetTPCindex(i);
465 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
466 iotrack.SetV0Indexes(pt->GetV0Indexes());
467 // iotrack.SetTPCpid(pt->fTPCr);
468 fEvent->AddTrack(&iotrack);
472 // short tracks - maybe decays
474 if ( (pt->GetNumberOfClusters()>30) && (Float_t(pt->GetNumberOfClusters())/Float_t(pt->GetNFoundable()))>0.70) {
475 Int_t found,foundable,shared;
476 pt->GetClusterStatistic(0,60,found, foundable,shared,kFALSE);
477 if ( (found>20) && (pt->GetNShared()/float(pt->GetNumberOfClusters())<0.2)){
479 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
480 //iotrack.SetTPCindex(i);
481 iotrack.SetTPCPoints(pt->GetPoints());
482 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
483 iotrack.SetV0Indexes(pt->GetV0Indexes());
484 //iotrack.SetTPCpid(pt->fTPCr);
485 fEvent->AddTrack(&iotrack);
490 if ( (pt->GetNumberOfClusters()>20) && (Float_t(pt->GetNumberOfClusters())/Float_t(pt->GetNFoundable()))>0.8) {
491 Int_t found,foundable,shared;
492 pt->GetClusterStatistic(0,60,found, foundable,shared,kFALSE);
493 if (found<20) continue;
494 if (pt->GetNShared()/float(pt->GetNumberOfClusters())>0.2) continue;
497 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
498 iotrack.SetTPCPoints(pt->GetPoints());
499 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
500 iotrack.SetV0Indexes(pt->GetV0Indexes());
501 //iotrack.SetTPCpid(pt->fTPCr);
502 //iotrack.SetTPCindex(i);
503 fEvent->AddTrack(&iotrack);
506 // short tracks - secondaties
508 if ( (pt->GetNumberOfClusters()>30) ) {
509 Int_t found,foundable,shared;
510 pt->GetClusterStatistic(128,158,found, foundable,shared,kFALSE);
511 if ( (found>20) && (pt->GetNShared()/float(pt->GetNumberOfClusters())<0.2) &&float(found)/float(foundable)>0.8){
513 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
514 iotrack.SetTPCPoints(pt->GetPoints());
515 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
516 iotrack.SetV0Indexes(pt->GetV0Indexes());
517 //iotrack.SetTPCpid(pt->fTPCr);
518 //iotrack.SetTPCindex(i);
519 fEvent->AddTrack(&iotrack);
524 if ( (pt->GetNumberOfClusters()>15)) {
525 Int_t found,foundable,shared;
526 pt->GetClusterStatistic(138,158,found, foundable,shared,kFALSE);
527 if (found<15) continue;
528 if (foundable<=0) continue;
529 if (pt->GetNShared()/float(pt->GetNumberOfClusters())>0.2) continue;
530 if (float(found)/float(foundable)<0.8) continue;
533 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
534 iotrack.SetTPCPoints(pt->GetPoints());
535 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
536 iotrack.SetV0Indexes(pt->GetV0Indexes());
537 // iotrack.SetTPCpid(pt->fTPCr);
538 //iotrack.SetTPCindex(i);
539 fEvent->AddTrack(&iotrack);
544 printf("Number of filled ESDs-\t%d\n",fEvent->GetNumberOfTracks());
551 Double_t AliTPCtrackerMI::ErrY2(AliTPCseed* seed, AliTPCclusterMI * cl){
554 // Use calibrated cluster error from OCDB
556 AliTPCClusterParam * clparam = AliTPCcalibDB::Instance()->GetClusterParam();
558 Float_t z = TMath::Abs(fParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
559 Int_t ctype = cl->GetType();
560 Int_t type = (cl->GetRow()<63) ? 0: (cl->GetRow()>126) ? 1:2;
561 Double_t angle = seed->GetSnp()*seed->GetSnp();
562 angle = TMath::Sqrt(TMath::Abs(angle/(1.-angle)));
563 Double_t erry2 = clparam->GetError0Par(0,type, z,angle);
565 erry2+=0.5; // edge cluster
568 seed->SetErrorY2(erry2);
572 //calculate look-up table at the beginning
573 // static Bool_t ginit = kFALSE;
574 // static Float_t gnoise1,gnoise2,gnoise3;
575 // static Float_t ggg1[10000];
576 // static Float_t ggg2[10000];
577 // static Float_t ggg3[10000];
578 // static Float_t glandau1[10000];
579 // static Float_t glandau2[10000];
580 // static Float_t glandau3[10000];
582 // static Float_t gcor01[500];
583 // static Float_t gcor02[500];
584 // static Float_t gcorp[500];
588 // if (ginit==kFALSE){
589 // for (Int_t i=1;i<500;i++){
590 // Float_t rsigma = float(i)/100.;
591 // gcor02[i] = TMath::Max(0.78 +TMath::Exp(7.4*(rsigma-1.2)),0.6);
592 // gcor01[i] = TMath::Max(0.72 +TMath::Exp(3.36*(rsigma-1.2)),0.6);
593 // gcorp[i] = TMath::Max(TMath::Power((rsigma+0.5),1.5),1.2);
597 // for (Int_t i=3;i<10000;i++){
601 // Float_t amp = float(i);
602 // Float_t padlength =0.75;
603 // gnoise1 = 0.0004/padlength;
604 // Float_t nel = 0.268*amp;
605 // Float_t nprim = 0.155*amp;
606 // ggg1[i] = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.001*nel/(padlength*padlength))/nel;
607 // glandau1[i] = (2.+0.12*nprim)*0.5* (2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
608 // if (glandau1[i]>1) glandau1[i]=1;
609 // glandau1[i]*=padlength*padlength/12.;
613 // gnoise2 = 0.0004/padlength;
615 // nprim = 0.133*amp;
616 // ggg2[i] = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
617 // glandau2[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
618 // if (glandau2[i]>1) glandau2[i]=1;
619 // glandau2[i]*=padlength*padlength/12.;
624 // gnoise3 = 0.0004/padlength;
626 // nprim = 0.133*amp;
627 // ggg3[i] = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
628 // glandau3[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
629 // if (glandau3[i]>1) glandau3[i]=1;
630 // glandau3[i]*=padlength*padlength/12.;
638 // Int_t amp = int(TMath::Abs(cl->GetQ()));
640 // seed->SetErrorY2(1.);
644 // Float_t z = TMath::Abs(fParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
645 // Int_t ctype = cl->GetType();
646 // Float_t padlength= GetPadPitchLength(seed->GetRow());
647 // Double_t angle2 = seed->GetSnp()*seed->GetSnp();
648 // angle2 = angle2/(1-angle2);
650 // //cluster "quality"
651 // Int_t rsigmay = int(100.*cl->GetSigmaY2()/(seed->GetCurrentSigmaY2()));
654 // if (fSectors==fInnerSec){
655 // snoise2 = gnoise1;
656 // res = ggg1[amp]*z+glandau1[amp]*angle2;
657 // if (ctype==0) res *= gcor01[rsigmay];
660 // res*= gcorp[rsigmay];
664 // if (padlength<1.1){
665 // snoise2 = gnoise2;
666 // res = ggg2[amp]*z+glandau2[amp]*angle2;
667 // if (ctype==0) res *= gcor02[rsigmay];
670 // res*= gcorp[rsigmay];
674 // snoise2 = gnoise3;
675 // res = ggg3[amp]*z+glandau3[amp]*angle2;
676 // if (ctype==0) res *= gcor02[rsigmay];
679 // res*= gcorp[rsigmay];
686 // res*=2.4; // overestimate error 2 times
690 // if (res<2*snoise2)
693 // seed->SetErrorY2(res);
701 Double_t AliTPCtrackerMI::ErrZ2(AliTPCseed* seed, AliTPCclusterMI * cl){
704 // Use calibrated cluster error from OCDB
706 AliTPCClusterParam * clparam = AliTPCcalibDB::Instance()->GetClusterParam();
708 Float_t z = TMath::Abs(fParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
709 Int_t ctype = cl->GetType();
710 Int_t type = (cl->GetRow()<63) ? 0: (cl->GetRow()>126) ? 1:2;
712 Double_t angle2 = seed->GetSnp()*seed->GetSnp();
713 angle2 = seed->GetTgl()*seed->GetTgl()*(1+angle2/(1-angle2));
714 Double_t angle = TMath::Sqrt(TMath::Abs(angle2));
715 Double_t errz2 = clparam->GetError0Par(1,type, z,angle);
717 errz2+=0.5; // edge cluster
720 seed->SetErrorZ2(errz2);
726 // //seed->SetErrorY2(0.1);
728 // //calculate look-up table at the beginning
729 // static Bool_t ginit = kFALSE;
730 // static Float_t gnoise1,gnoise2,gnoise3;
731 // static Float_t ggg1[10000];
732 // static Float_t ggg2[10000];
733 // static Float_t ggg3[10000];
734 // static Float_t glandau1[10000];
735 // static Float_t glandau2[10000];
736 // static Float_t glandau3[10000];
738 // static Float_t gcor01[1000];
739 // static Float_t gcor02[1000];
740 // static Float_t gcorp[1000];
744 // if (ginit==kFALSE){
745 // for (Int_t i=1;i<1000;i++){
746 // Float_t rsigma = float(i)/100.;
747 // gcor02[i] = TMath::Max(0.81 +TMath::Exp(6.8*(rsigma-1.2)),0.6);
748 // gcor01[i] = TMath::Max(0.72 +TMath::Exp(2.04*(rsigma-1.2)),0.6);
749 // gcorp[i] = TMath::Max(TMath::Power((rsigma+0.5),1.5),1.2);
753 // for (Int_t i=3;i<10000;i++){
757 // Float_t amp = float(i);
758 // Float_t padlength =0.75;
759 // gnoise1 = 0.0004/padlength;
760 // Float_t nel = 0.268*amp;
761 // Float_t nprim = 0.155*amp;
762 // ggg1[i] = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.001*nel/(padlength*padlength))/nel;
763 // glandau1[i] = (2.+0.12*nprim)*0.5* (2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
764 // if (glandau1[i]>1) glandau1[i]=1;
765 // glandau1[i]*=padlength*padlength/12.;
769 // gnoise2 = 0.0004/padlength;
771 // nprim = 0.133*amp;
772 // ggg2[i] = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
773 // glandau2[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
774 // if (glandau2[i]>1) glandau2[i]=1;
775 // glandau2[i]*=padlength*padlength/12.;
780 // gnoise3 = 0.0004/padlength;
782 // nprim = 0.133*amp;
783 // ggg3[i] = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
784 // glandau3[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
785 // if (glandau3[i]>1) glandau3[i]=1;
786 // glandau3[i]*=padlength*padlength/12.;
794 // Int_t amp = int(TMath::Abs(cl->GetQ()));
796 // seed->SetErrorY2(1.);
800 // Float_t z = TMath::Abs(fParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
801 // Int_t ctype = cl->GetType();
802 // Float_t padlength= GetPadPitchLength(seed->GetRow());
804 // Double_t angle2 = seed->GetSnp()*seed->GetSnp();
805 // // if (angle2<0.6) angle2 = 0.6;
806 // angle2 = seed->GetTgl()*seed->GetTgl()*(1+angle2/(1-angle2));
808 // //cluster "quality"
809 // Int_t rsigmaz = int(100.*cl->GetSigmaZ2()/(seed->GetCurrentSigmaZ2()));
812 // if (fSectors==fInnerSec){
813 // snoise2 = gnoise1;
814 // res = ggg1[amp]*z+glandau1[amp]*angle2;
815 // if (ctype==0) res *= gcor01[rsigmaz];
818 // res*= gcorp[rsigmaz];
822 // if (padlength<1.1){
823 // snoise2 = gnoise2;
824 // res = ggg2[amp]*z+glandau2[amp]*angle2;
825 // if (ctype==0) res *= gcor02[rsigmaz];
828 // res*= gcorp[rsigmaz];
832 // snoise2 = gnoise3;
833 // res = ggg3[amp]*z+glandau3[amp]*angle2;
834 // if (ctype==0) res *= gcor02[rsigmaz];
837 // res*= gcorp[rsigmaz];
846 // if ((ctype<0) &&<70){
851 // if (res<2*snoise2)
853 // if (res>3) res =3;
854 // seed->SetErrorZ2(res);
862 void AliTPCtrackerMI::RotateToLocal(AliTPCseed *seed)
864 //rotate to track "local coordinata
865 Float_t x = seed->GetX();
866 Float_t y = seed->GetY();
867 Float_t ymax = x*TMath::Tan(0.5*fSectors->GetAlpha());
870 seed->SetRelativeSector((seed->GetRelativeSector()+1) % fN);
871 if (!seed->Rotate(fSectors->GetAlpha()))
873 } else if (y <-ymax) {
874 seed->SetRelativeSector((seed->GetRelativeSector()-1+fN) % fN);
875 if (!seed->Rotate(-fSectors->GetAlpha()))
883 //_____________________________________________________________________________
884 Double_t AliTPCtrackerMI::F1old(Double_t x1,Double_t y1,
885 Double_t x2,Double_t y2,
886 Double_t x3,Double_t y3)
888 //-----------------------------------------------------------------
889 // Initial approximation of the track curvature
890 //-----------------------------------------------------------------
891 Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
892 Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
893 (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
894 Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
895 (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
897 Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
898 if ( xr*xr+yr*yr<=0.00000000000001) return 100;
899 return -xr*yr/sqrt(xr*xr+yr*yr);
904 //_____________________________________________________________________________
905 Double_t AliTPCtrackerMI::F1(Double_t x1,Double_t y1,
906 Double_t x2,Double_t y2,
907 Double_t x3,Double_t y3)
909 //-----------------------------------------------------------------
910 // Initial approximation of the track curvature
911 //-----------------------------------------------------------------
917 Double_t det = x3*y2-x2*y3;
922 Double_t u = 0.5* (x2*(x2-x3)+y2*(y2-y3))/det;
923 Double_t x0 = x3*0.5-y3*u;
924 Double_t y0 = y3*0.5+x3*u;
925 Double_t c2 = 1/TMath::Sqrt(x0*x0+y0*y0);
931 Double_t AliTPCtrackerMI::F2(Double_t x1,Double_t y1,
932 Double_t x2,Double_t y2,
933 Double_t x3,Double_t y3)
935 //-----------------------------------------------------------------
936 // Initial approximation of the track curvature
937 //-----------------------------------------------------------------
943 Double_t det = x3*y2-x2*y3;
948 Double_t u = 0.5* (x2*(x2-x3)+y2*(y2-y3))/det;
949 Double_t x0 = x3*0.5-y3*u;
950 Double_t y0 = y3*0.5+x3*u;
951 Double_t c2 = 1/TMath::Sqrt(x0*x0+y0*y0);
960 //_____________________________________________________________________________
961 Double_t AliTPCtrackerMI::F2old(Double_t x1,Double_t y1,
962 Double_t x2,Double_t y2,
963 Double_t x3,Double_t y3)
965 //-----------------------------------------------------------------
966 // Initial approximation of the track curvature times center of curvature
967 //-----------------------------------------------------------------
968 Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
969 Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
970 (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
971 Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
972 (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
974 Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
976 return -a/(d*y1-b)*xr/sqrt(xr*xr+yr*yr);
979 //_____________________________________________________________________________
980 Double_t AliTPCtrackerMI::F3(Double_t x1,Double_t y1,
981 Double_t x2,Double_t y2,
982 Double_t z1,Double_t z2)
984 //-----------------------------------------------------------------
985 // Initial approximation of the tangent of the track dip angle
986 //-----------------------------------------------------------------
987 return (z1 - z2)/sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
991 Double_t AliTPCtrackerMI::F3n(Double_t x1,Double_t y1,
992 Double_t x2,Double_t y2,
993 Double_t z1,Double_t z2, Double_t c)
995 //-----------------------------------------------------------------
996 // Initial approximation of the tangent of the track dip angle
997 //-----------------------------------------------------------------
1001 //angle1 = (z1-z2)*c/(TMath::ASin(c*x1-ni)-TMath::ASin(c*x2-ni));
1003 Double_t d = TMath::Sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
1004 if (TMath::Abs(d*c*0.5)>1) return 0;
1005 // Double_t angle2 = TMath::ASin(d*c*0.5);
1006 // Double_t angle2 = AliTPCFastMath::FastAsin(d*c*0.5);
1007 Double_t angle2 = (d*c*0.5>0.1)? TMath::ASin(d*c*0.5): AliTPCFastMath::FastAsin(d*c*0.5);
1009 angle2 = (z1-z2)*c/(angle2*2.);
1013 Bool_t AliTPCtrackerMI::GetProlongation(Double_t x1, Double_t x2, Double_t x[5], Double_t &y, Double_t &z)
1014 {//-----------------------------------------------------------------
1015 // This function find proloncation of a track to a reference plane x=x2.
1016 //-----------------------------------------------------------------
1020 if (TMath::Abs(x[4]*x1 - x[2]) >= 0.999) {
1024 Double_t c1=x[4]*x1 - x[2], r1=sqrt(1.- c1*c1);
1025 Double_t c2=x[4]*x2 - x[2], r2=sqrt(1.- c2*c2);
1029 Double_t dy = dx*(c1+c2)/(r1+r2);
1032 Double_t delta = x[4]*dx*(c1+c2)/(c1*r2 + c2*r1);
1034 if (TMath::Abs(delta)>0.01){
1035 dz = x[3]*TMath::ASin(delta)/x[4];
1037 dz = x[3]*AliTPCFastMath::FastAsin(delta)/x[4];
1040 //dz = x[3]*AliTPCFastMath::FastAsin(delta)/x[4];
1048 Int_t AliTPCtrackerMI::LoadClusters (TTree *tree)
1053 return LoadClusters();
1057 Int_t AliTPCtrackerMI::LoadClusters(TObjArray *arr)
1060 // load clusters to the memory
1061 AliTPCClustersRow *clrow = 0x0;
1062 Int_t lower = arr->LowerBound();
1063 Int_t entries = arr->GetEntriesFast();
1065 for (Int_t i=lower; i<entries; i++) {
1066 clrow = (AliTPCClustersRow*) arr->At(i);
1067 if(!clrow) continue;
1068 if(!clrow->GetArray()) continue;
1072 fParam->AdjustSectorRow(clrow->GetID(),sec,row);
1074 for (Int_t icl=0; icl<clrow->GetArray()->GetEntriesFast(); icl++){
1075 Transform((AliTPCclusterMI*)(clrow->GetArray()->At(icl)));
1078 if (clrow->GetArray()->GetEntriesFast()<=0) continue;
1079 AliTPCtrackerRow * tpcrow=0;
1082 tpcrow = &(fInnerSec[sec%fkNIS][row]);
1086 tpcrow = &(fOuterSec[(sec-fkNIS*2)%fkNOS][row]);
1087 left = (sec-fkNIS*2)/fkNOS;
1090 tpcrow->SetN1(clrow->GetArray()->GetEntriesFast());
1091 tpcrow->SetClusters1(new AliTPCclusterMI[tpcrow->GetN1()]);
1092 for (Int_t i=0;i<tpcrow->GetN1();i++)
1093 tpcrow->SetCluster1(i, *(AliTPCclusterMI*)(clrow->GetArray()->At(i)));
1096 tpcrow->SetN2(clrow->GetArray()->GetEntriesFast());
1097 tpcrow->SetClusters2(new AliTPCclusterMI[tpcrow->GetN2()]);
1098 for (Int_t i=0;i<tpcrow->GetN2();i++)
1099 tpcrow->SetCluster2(i,*(AliTPCclusterMI*)(clrow->GetArray()->At(i)));
1110 Int_t AliTPCtrackerMI::LoadClusters()
1113 // load clusters to the memory
1114 AliTPCClustersRow *clrow= new AliTPCClustersRow;
1115 clrow->SetClass("AliTPCclusterMI");
1117 clrow->GetArray()->ExpandCreateFast(10000);
1119 // TTree * tree = fClustersArray.GetTree();
1121 TTree * tree = fInput;
1122 TBranch * br = tree->GetBranch("Segment");
1123 br->SetAddress(&clrow);
1125 Int_t j=Int_t(tree->GetEntries());
1126 for (Int_t i=0; i<j; i++) {
1130 fParam->AdjustSectorRow(clrow->GetID(),sec,row);
1131 for (Int_t icl=0; icl<clrow->GetArray()->GetEntriesFast(); icl++){
1132 Transform((AliTPCclusterMI*)(clrow->GetArray()->At(icl)));
1135 AliTPCtrackerRow * tpcrow=0;
1138 tpcrow = &(fInnerSec[sec%fkNIS][row]);
1142 tpcrow = &(fOuterSec[(sec-fkNIS*2)%fkNOS][row]);
1143 left = (sec-fkNIS*2)/fkNOS;
1146 tpcrow->SetN1(clrow->GetArray()->GetEntriesFast());
1147 tpcrow->SetClusters1(new AliTPCclusterMI[tpcrow->GetN1()]);
1148 for (Int_t i=0;i<tpcrow->GetN1();i++)
1149 tpcrow->SetCluster1(i, *(AliTPCclusterMI*)(clrow->GetArray()->At(i)));
1152 tpcrow->SetN2(clrow->GetArray()->GetEntriesFast());
1153 tpcrow->SetClusters2(new AliTPCclusterMI[tpcrow->GetN2()]);
1154 for (Int_t i=0;i<tpcrow->GetN2();i++)
1155 tpcrow->SetCluster2(i,*(AliTPCclusterMI*)(clrow->GetArray()->At(i)));
1166 void AliTPCtrackerMI::UnloadClusters()
1169 // unload clusters from the memory
1171 Int_t nrows = fOuterSec->GetNRows();
1172 for (Int_t sec = 0;sec<fkNOS;sec++)
1173 for (Int_t row = 0;row<nrows;row++){
1174 AliTPCtrackerRow* tpcrow = &(fOuterSec[sec%fkNOS][row]);
1176 // if (tpcrow->fClusters1) delete []tpcrow->fClusters1;
1177 // if (tpcrow->fClusters2) delete []tpcrow->fClusters2;
1179 tpcrow->ResetClusters();
1182 nrows = fInnerSec->GetNRows();
1183 for (Int_t sec = 0;sec<fkNIS;sec++)
1184 for (Int_t row = 0;row<nrows;row++){
1185 AliTPCtrackerRow* tpcrow = &(fInnerSec[sec%fkNIS][row]);
1187 // if (tpcrow->fClusters1) delete []tpcrow->fClusters1;
1188 //if (tpcrow->fClusters2) delete []tpcrow->fClusters2;
1190 tpcrow->ResetClusters();
1196 void AliTPCtrackerMI::Transform(AliTPCclusterMI * cluster){
1200 AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ;
1202 AliFatal("Tranformations not in calibDB");
1204 Double_t x[3]={cluster->GetRow(),cluster->GetPad(),cluster->GetTimeBin()};
1205 Int_t i[1]={cluster->GetDetector()};
1206 transform->Transform(x,i,0,1);
1207 if (!AliTPCReconstructor::GetRecoParam()->GetBYMirror()){
1208 if (cluster->GetDetector()%36>17){
1214 // in debug mode check the transformation
1216 if (AliTPCReconstructor::StreamLevel()>1) {
1218 cluster->GetGlobalXYZ(gx);
1219 TTreeSRedirector &cstream = *fDebugStreamer;
1220 cstream<<"Transform"<<
1230 cluster->SetX(x[0]);
1231 cluster->SetY(x[1]);
1232 cluster->SetZ(x[2]);
1238 //if (!fParam->IsGeoRead()) fParam->ReadGeoMatrices();
1239 TGeoHMatrix *mat = fParam->GetClusterMatrix(cluster->GetDetector());
1241 Double_t pos[3]= {cluster->GetX(),cluster->GetY(),cluster->GetZ()};
1242 Double_t posC[3]={cluster->GetX(),cluster->GetY(),cluster->GetZ()};
1243 if (mat) mat->LocalToMaster(pos,posC);
1245 // chack Loading of Geo matrices from GeoManager - TEMPORARY FIX
1247 cluster->SetX(posC[0]);
1248 cluster->SetY(posC[1]);
1249 cluster->SetZ(posC[2]);
1252 //_____________________________________________________________________________
1253 Int_t AliTPCtrackerMI::LoadOuterSectors() {
1254 //-----------------------------------------------------------------
1255 // This function fills outer TPC sectors with clusters.
1256 //-----------------------------------------------------------------
1257 Int_t nrows = fOuterSec->GetNRows();
1259 for (Int_t sec = 0;sec<fkNOS;sec++)
1260 for (Int_t row = 0;row<nrows;row++){
1261 AliTPCtrackerRow* tpcrow = &(fOuterSec[sec%fkNOS][row]);
1262 Int_t sec2 = sec+2*fkNIS;
1264 Int_t ncl = tpcrow->GetN1();
1266 AliTPCclusterMI *c= (tpcrow->GetCluster1(ncl));
1267 index=(((sec2<<8)+row)<<16)+ncl;
1268 tpcrow->InsertCluster(c,index);
1271 ncl = tpcrow->GetN2();
1273 AliTPCclusterMI *c= (tpcrow->GetCluster2(ncl));
1274 index=((((sec2+fkNOS)<<8)+row)<<16)+ncl;
1275 tpcrow->InsertCluster(c,index);
1278 // write indexes for fast acces
1280 for (Int_t i=0;i<510;i++)
1281 tpcrow->SetFastCluster(i,-1);
1282 for (Int_t i=0;i<tpcrow->GetN();i++){
1283 Int_t zi = Int_t((*tpcrow)[i]->GetZ()+255.);
1284 tpcrow->SetFastCluster(zi,i); // write index
1287 for (Int_t i=0;i<510;i++){
1288 if (tpcrow->GetFastCluster(i)<0)
1289 tpcrow->SetFastCluster(i,last);
1291 last = tpcrow->GetFastCluster(i);
1300 //_____________________________________________________________________________
1301 Int_t AliTPCtrackerMI::LoadInnerSectors() {
1302 //-----------------------------------------------------------------
1303 // This function fills inner TPC sectors with clusters.
1304 //-----------------------------------------------------------------
1305 Int_t nrows = fInnerSec->GetNRows();
1307 for (Int_t sec = 0;sec<fkNIS;sec++)
1308 for (Int_t row = 0;row<nrows;row++){
1309 AliTPCtrackerRow* tpcrow = &(fInnerSec[sec%fkNIS][row]);
1312 Int_t ncl = tpcrow->GetN1();
1314 AliTPCclusterMI *c= (tpcrow->GetCluster1(ncl));
1315 index=(((sec<<8)+row)<<16)+ncl;
1316 tpcrow->InsertCluster(c,index);
1319 ncl = tpcrow->GetN2();
1321 AliTPCclusterMI *c= (tpcrow->GetCluster2(ncl));
1322 index=((((sec+fkNIS)<<8)+row)<<16)+ncl;
1323 tpcrow->InsertCluster(c,index);
1326 // write indexes for fast acces
1328 for (Int_t i=0;i<510;i++)
1329 tpcrow->SetFastCluster(i,-1);
1330 for (Int_t i=0;i<tpcrow->GetN();i++){
1331 Int_t zi = Int_t((*tpcrow)[i]->GetZ()+255.);
1332 tpcrow->SetFastCluster(zi,i); // write index
1335 for (Int_t i=0;i<510;i++){
1336 if (tpcrow->GetFastCluster(i)<0)
1337 tpcrow->SetFastCluster(i,last);
1339 last = tpcrow->GetFastCluster(i);
1351 //_________________________________________________________________________
1352 AliTPCclusterMI *AliTPCtrackerMI::GetClusterMI(Int_t index) const {
1353 //--------------------------------------------------------------------
1354 // Return pointer to a given cluster
1355 //--------------------------------------------------------------------
1356 if (index<0) return 0; // no cluster
1357 Int_t sec=(index&0xff000000)>>24;
1358 Int_t row=(index&0x00ff0000)>>16;
1359 Int_t ncl=(index&0x00007fff)>>00;
1361 const AliTPCtrackerRow * tpcrow=0;
1362 AliTPCclusterMI * clrow =0;
1364 if (sec<0 || sec>=fkNIS*4) {
1365 AliWarning(Form("Wrong sector %d",sec));
1370 tpcrow = &(fInnerSec[sec%fkNIS][row]);
1371 if (tpcrow==0) return 0;
1374 if (tpcrow->GetN1()<=ncl) return 0;
1375 clrow = tpcrow->GetClusters1();
1378 if (tpcrow->GetN2()<=ncl) return 0;
1379 clrow = tpcrow->GetClusters2();
1383 tpcrow = &(fOuterSec[(sec-fkNIS*2)%fkNOS][row]);
1384 if (tpcrow==0) return 0;
1386 if (sec-2*fkNIS<fkNOS) {
1387 if (tpcrow->GetN1()<=ncl) return 0;
1388 clrow = tpcrow->GetClusters1();
1391 if (tpcrow->GetN2()<=ncl) return 0;
1392 clrow = tpcrow->GetClusters2();
1396 return &(clrow[ncl]);
1402 Int_t AliTPCtrackerMI::FollowToNext(AliTPCseed& t, Int_t nr) {
1403 //-----------------------------------------------------------------
1404 // This function tries to find a track prolongation to next pad row
1405 //-----------------------------------------------------------------
1407 Double_t x= GetXrow(nr), ymax=GetMaxY(nr);
1408 AliTPCclusterMI *cl=0;
1409 Int_t tpcindex= t.GetClusterIndex2(nr);
1411 // update current shape info every 5 pad-row
1412 // if ( (nr%5==0) || t.GetNumberOfClusters()<2 || (t.fCurrentSigmaY2<0.0001) ){
1416 if (fIteration>0 && tpcindex>=-1){ //if we have already clusters
1418 if (tpcindex==-1) return 0; //track in dead zone
1420 cl = t.GetClusterPointer(nr);
1421 if ( (cl==0) ) cl = GetClusterMI(tpcindex);
1422 t.SetCurrentClusterIndex1(tpcindex);
1425 Int_t relativesector = ((tpcindex&0xff000000)>>24)%18; // if previously accepted cluster in different sector
1426 Float_t angle = relativesector*fSectors->GetAlpha()+fSectors->GetAlphaShift();
1428 if (angle<-TMath::Pi()) angle += 2*TMath::Pi();
1429 if (angle>=TMath::Pi()) angle -= 2*TMath::Pi();
1431 if (TMath::Abs(angle-t.GetAlpha())>0.001){
1432 Double_t rotation = angle-t.GetAlpha();
1433 t.SetRelativeSector(relativesector);
1434 if (!t.Rotate(rotation)) return 0;
1436 if (!t.PropagateTo(x)) return 0;
1438 t.SetCurrentCluster(cl);
1440 Int_t accept = AcceptCluster(&t,t.GetCurrentCluster());
1441 if ((tpcindex&0x8000)==0) accept =0;
1443 //if founded cluster is acceptible
1444 if (cl->IsUsed(11)) { // id cluster is shared inrease uncertainty
1445 t.SetErrorY2(t.GetErrorY2()+0.03);
1446 t.SetErrorZ2(t.GetErrorZ2()+0.03);
1447 t.SetErrorY2(t.GetErrorY2()*3);
1448 t.SetErrorZ2(t.GetErrorZ2()*3);
1450 t.SetNFoundable(t.GetNFoundable()+1);
1451 UpdateTrack(&t,accept);
1456 if (TMath::Abs(t.GetSnp())>AliTPCReconstructor::GetMaxSnpTracker()) return 0; // cut on angle
1458 // not look for new cluster during refitting
1459 t.SetNFoundable(t.GetNFoundable()+1);
1464 // if (TMath::Abs(t.GetSnp())>0.95 || TMath::Abs(x*t.GetC()-t.GetEta())>0.95) return 0;// patch 28 fev 06
1465 Double_t y=t.GetYat(x);
1466 if (TMath::Abs(y)>ymax){
1468 t.SetRelativeSector((t.GetRelativeSector()+1) % fN);
1469 if (!t.Rotate(fSectors->GetAlpha()))
1471 } else if (y <-ymax) {
1472 t.SetRelativeSector((t.GetRelativeSector()-1+fN) % fN);
1473 if (!t.Rotate(-fSectors->GetAlpha()))
1479 if (!t.PropagateTo(x)) {
1480 if (fIteration==0) t.SetRemoval(10);
1484 Double_t z=t.GetZ();
1487 if (!IsActive(t.GetRelativeSector(),nr)) {
1489 t.SetClusterIndex2(nr,-1);
1492 //AliInfo(Form("A - Sector%d phi %f - alpha %f", t.fRelativeSector,y/x, t.GetAlpha()));
1493 Bool_t isActive = IsActive(t.GetRelativeSector(),nr);
1494 Bool_t isActive2 = (nr>=fInnerSec->GetNRows()) ? fOuterSec[t.GetRelativeSector()][nr-fInnerSec->GetNRows()].GetN()>0:fInnerSec[t.GetRelativeSector()][nr].GetN()>0;
1496 if (!isActive || !isActive2) return 0;
1498 const AliTPCtrackerRow &krow=GetRow(t.GetRelativeSector(),nr);
1499 if ( (t.GetSigmaY2()<0) || t.GetSigmaZ2()<0) return 0;
1501 Double_t roadz = 1.;
1503 if (TMath::Abs(TMath::Abs(y)-ymax)<krow.GetDeadZone()){
1505 t.SetClusterIndex2(nr,-1);
1510 if (TMath::Abs(z)<(AliTPCReconstructor::GetCtgRange()*x+10) && TMath::Abs(z)<fParam->GetZLength(0) && (TMath::Abs(t.GetSnp())<AliTPCReconstructor::GetMaxSnpTracker()))
1511 t.SetNFoundable(t.GetNFoundable()+1);
1517 // cl = krow.FindNearest2(y+10.,z,roady,roadz,index);
1518 cl = krow.FindNearest2(y,z,roady,roadz,index);
1519 if (cl) t.SetCurrentClusterIndex1(krow.GetIndex(index));
1522 t.SetCurrentCluster(cl);
1524 if (fIteration==2&&cl->IsUsed(10)) return 0;
1525 Int_t accept = AcceptCluster(&t,t.GetCurrentCluster());
1526 if (fIteration==2&&cl->IsUsed(11)) {
1527 t.SetErrorY2(t.GetErrorY2()+0.03);
1528 t.SetErrorZ2(t.GetErrorZ2()+0.03);
1529 t.SetErrorY2(t.GetErrorY2()*3);
1530 t.SetErrorZ2(t.GetErrorZ2()*3);
1533 if (t.fCurrentCluster->IsUsed(10)){
1538 if (t.fNShared>0.7*t.GetNumberOfClusters()) {
1544 if (accept<3) UpdateTrack(&t,accept);
1547 if ( fIteration==0 && t.GetNFoundable()*0.5 > t.GetNumberOfClusters()) t.SetRemoval(10);
1553 Int_t AliTPCtrackerMI::FollowToNextFast(AliTPCseed& t, Int_t nr) {
1554 //-----------------------------------------------------------------
1555 // This function tries to find a track prolongation to next pad row
1556 //-----------------------------------------------------------------
1558 Double_t x= GetXrow(nr), ymax=GetMaxY(nr);
1560 if (!t.GetProlongation(x,y,z)) {
1566 if (TMath::Abs(y)>ymax){
1569 t.SetRelativeSector((t.GetRelativeSector()+1) % fN);
1570 if (!t.Rotate(fSectors->GetAlpha()))
1572 } else if (y <-ymax) {
1573 t.SetRelativeSector((t.GetRelativeSector()-1+fN) % fN);
1574 if (!t.Rotate(-fSectors->GetAlpha()))
1577 if (!t.PropagateTo(x)) {
1580 t.GetProlongation(x,y,z);
1583 // update current shape info every 2 pad-row
1584 if ( (nr%2==0) || t.GetNumberOfClusters()<2 || (t.GetCurrentSigmaY2()<0.0001) ){
1585 // t.fCurrentSigmaY = GetSigmaY(&t);
1586 //t.fCurrentSigmaZ = GetSigmaZ(&t);
1590 AliTPCclusterMI *cl=0;
1595 const AliTPCtrackerRow &krow=GetRow(t.GetRelativeSector(),nr);
1596 if ( (t.GetSigmaY2()<0) || t.GetSigmaZ2()<0) return 0;
1598 Double_t roadz = 1.;
1601 if (TMath::Abs(TMath::Abs(y)-ymax)<krow.GetDeadZone()){
1603 t.SetClusterIndex2(row,-1);
1608 if (TMath::Abs(z)>(AliTPCReconstructor::GetCtgRange()*x+10)) t.SetClusterIndex2(row,-1);
1612 if ((cl==0)&&(krow)) {
1613 // cl = krow.FindNearest2(y+10,z,roady,roadz,index);
1614 cl = krow.FindNearest2(y,z,roady,roadz,index);
1616 if (cl) t.SetCurrentClusterIndex1(krow.GetIndex(index));
1620 t.SetCurrentCluster(cl);
1621 // Int_t accept = AcceptCluster(&t,t.fCurrentCluster);
1623 t.SetClusterIndex2(row,index);
1624 t.SetClusterPointer(row, cl);
1632 //_________________________________________________________________________
1633 Bool_t AliTPCtrackerMI::GetTrackPoint(Int_t index, AliTrackPoint &p ) const
1635 // Get track space point by index
1636 // return false in case the cluster doesn't exist
1637 AliTPCclusterMI *cl = GetClusterMI(index);
1638 if (!cl) return kFALSE;
1639 Int_t sector = (index&0xff000000)>>24;
1640 // Int_t row = (index&0x00ff0000)>>16;
1642 // xyz[0] = fParam->GetPadRowRadii(sector,row);
1643 xyz[0] = cl->GetX();
1644 xyz[1] = cl->GetY();
1645 xyz[2] = cl->GetZ();
1647 fParam->AdjustCosSin(sector,cos,sin);
1648 Float_t x = cos*xyz[0]-sin*xyz[1];
1649 Float_t y = cos*xyz[1]+sin*xyz[0];
1651 Float_t sigmaY2 = 0.027*cl->GetSigmaY2();
1652 if (sector < fParam->GetNInnerSector()) sigmaY2 *= 2.07;
1653 Float_t sigmaZ2 = 0.066*cl->GetSigmaZ2();
1654 if (sector < fParam->GetNInnerSector()) sigmaZ2 *= 1.77;
1655 cov[0] = sin*sin*sigmaY2;
1656 cov[1] = -sin*cos*sigmaY2;
1658 cov[3] = cos*cos*sigmaY2;
1661 p.SetXYZ(x,y,xyz[2],cov);
1662 AliGeomManager::ELayerID iLayer;
1664 if (sector < fParam->GetNInnerSector()) {
1665 iLayer = AliGeomManager::kTPC1;
1669 iLayer = AliGeomManager::kTPC2;
1670 idet = sector - fParam->GetNInnerSector();
1672 UShort_t volid = AliGeomManager::LayerToVolUID(iLayer,idet);
1673 p.SetVolumeID(volid);
1679 Int_t AliTPCtrackerMI::UpdateClusters(AliTPCseed& t, Int_t nr) {
1680 //-----------------------------------------------------------------
1681 // This function tries to find a track prolongation to next pad row
1682 //-----------------------------------------------------------------
1683 t.SetCurrentCluster(0);
1684 t.SetCurrentClusterIndex1(0);
1686 Double_t xt=t.GetX();
1687 Int_t row = GetRowNumber(xt)-1;
1688 Double_t ymax= GetMaxY(nr);
1690 if (row < nr) return 1; // don't prolongate if not information until now -
1691 // if (TMath::Abs(t.GetSnp())>0.9 && t.GetNumberOfClusters()>40. && fIteration!=2) {
1693 // return 0; // not prolongate strongly inclined tracks
1695 // if (TMath::Abs(t.GetSnp())>0.95) {
1697 // return 0; // not prolongate strongly inclined tracks
1698 // }// patch 28 fev 06
1700 Double_t x= GetXrow(nr);
1702 //t.PropagateTo(x+0.02);
1703 //t.PropagateTo(x+0.01);
1704 if (!t.PropagateTo(x)){
1711 if (TMath::Abs(y)>ymax){
1713 t.SetRelativeSector((t.GetRelativeSector()+1) % fN);
1714 if (!t.Rotate(fSectors->GetAlpha()))
1716 } else if (y <-ymax) {
1717 t.SetRelativeSector((t.GetRelativeSector()-1+fN) % fN);
1718 if (!t.Rotate(-fSectors->GetAlpha()))
1721 // if (!t.PropagateTo(x)){
1728 if (TMath::Abs(t.GetSnp())>AliTPCReconstructor::GetMaxSnpTracker()) return 0;
1730 if (!IsActive(t.GetRelativeSector(),nr)) {
1732 t.SetClusterIndex2(nr,-1);
1735 //AliInfo(Form("A - Sector%d phi %f - alpha %f", t.fRelativeSector,y/x, t.GetAlpha()));
1737 AliTPCtrackerRow &krow=GetRow(t.GetRelativeSector(),nr);
1739 if (TMath::Abs(TMath::Abs(y)-ymax)<krow.GetDeadZone()){
1741 t.SetClusterIndex2(nr,-1);
1746 if (TMath::Abs(t.GetZ())<(AliTPCReconstructor::GetCtgRange()*t.GetX()+10) && (TMath::Abs(t.GetSnp())<AliTPCReconstructor::GetMaxSnpTracker())) t.SetNFoundable(t.GetNFoundable()+1);
1752 if ( (nr%2==0) || t.GetNumberOfClusters()<2){
1753 // t.fCurrentSigmaY = GetSigmaY(&t);
1754 //t.fCurrentSigmaZ = GetSigmaZ(&t);
1758 AliTPCclusterMI *cl=0;
1761 Double_t roady = 1.;
1762 Double_t roadz = 1.;
1766 index = t.GetClusterIndex2(nr);
1767 if ( (index>0) && (index&0x8000)==0){
1768 cl = t.GetClusterPointer(nr);
1769 if ( (cl==0) && (index>0)) cl = GetClusterMI(index);
1770 t.SetCurrentClusterIndex1(index);
1772 t.SetCurrentCluster(cl);
1778 // if (index<0) return 0;
1779 UInt_t uindex = TMath::Abs(index);
1782 //cl = krow.FindNearest2(y+10,z,roady,roadz,uindex);
1783 cl = krow.FindNearest2(y,z,roady,roadz,uindex);
1786 if (cl) t.SetCurrentClusterIndex1(krow.GetIndex(uindex));
1787 t.SetCurrentCluster(cl);
1793 Int_t AliTPCtrackerMI::FollowToNextCluster(AliTPCseed & t, Int_t nr) {
1794 //-----------------------------------------------------------------
1795 // This function tries to find a track prolongation to next pad row
1796 //-----------------------------------------------------------------
1798 //update error according neighborhoud
1800 if (t.GetCurrentCluster()) {
1802 Int_t accept = AcceptCluster(&t,t.GetCurrentCluster());
1804 if (t.GetCurrentCluster()->IsUsed(10)){
1809 t.SetNShared(t.GetNShared()+1);
1810 if (t.GetNShared()>0.7*t.GetNumberOfClusters()) {
1815 if (fIteration>0) accept = 0;
1816 if (accept<3) UpdateTrack(&t,accept);
1820 if ( ( (t.GetSigmaY2()+t.GetSigmaZ2())>0.16)&& t.GetNumberOfClusters()>18) t.SetRemoval(10);
1821 if ( t.GetChi2()/t.GetNumberOfClusters()>6 &&t.GetNumberOfClusters()>18) t.SetRemoval(10);
1823 if (( (t.GetNFoundable()*0.5 > t.GetNumberOfClusters()) || t.GetNoCluster()>15)) t.SetRemoval(10);
1831 //_____________________________________________________________________________
1832 Int_t AliTPCtrackerMI::FollowProlongation(AliTPCseed& t, Int_t rf, Int_t step) {
1833 //-----------------------------------------------------------------
1834 // This function tries to find a track prolongation.
1835 //-----------------------------------------------------------------
1836 Double_t xt=t.GetX();
1838 Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
1839 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
1840 if (alpha < 0. ) alpha += 2.*TMath::Pi();
1842 t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN);
1844 Int_t first = GetRowNumber(xt)-1;
1845 for (Int_t nr= first; nr>=rf; nr-=step) {
1847 if (t.GetKinkIndexes()[0]>0){
1848 for (Int_t i=0;i<3;i++){
1849 Int_t index = t.GetKinkIndexes()[i];
1850 if (index==0) break;
1851 if (index<0) continue;
1853 AliKink * kink = (AliKink*)fEvent->GetKink(index-1);
1855 printf("PROBLEM\n");
1858 Int_t kinkrow = kink->GetTPCRow0()+2+Int_t(0.5/(0.05+kink->GetAngle(2)));
1860 AliExternalTrackParam paramd(t);
1861 kink->SetDaughter(paramd);
1862 kink->SetStatus(2,5);
1869 if (nr==80) t.UpdateReference();
1870 if (nr<fInnerSec->GetNRows())
1871 fSectors = fInnerSec;
1873 fSectors = fOuterSec;
1874 if (FollowToNext(t,nr)==0)
1883 //_____________________________________________________________________________
1884 Int_t AliTPCtrackerMI::FollowProlongationFast(AliTPCseed& t, Int_t rf, Int_t step) {
1885 //-----------------------------------------------------------------
1886 // This function tries to find a track prolongation.
1887 //-----------------------------------------------------------------
1888 Double_t xt=t.GetX();
1890 Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
1891 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
1892 if (alpha < 0. ) alpha += 2.*TMath::Pi();
1893 t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN);
1895 for (Int_t nr=GetRowNumber(xt)-1; nr>=rf; nr-=step) {
1897 if (FollowToNextFast(t,nr)==0)
1898 if (!t.IsActive()) return 0;
1908 Int_t AliTPCtrackerMI::FollowBackProlongation(AliTPCseed& t, Int_t rf) {
1909 //-----------------------------------------------------------------
1910 // This function tries to find a track prolongation.
1911 //-----------------------------------------------------------------
1913 Double_t xt=t.GetX();
1914 Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
1915 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
1916 if (alpha < 0. ) alpha += 2.*TMath::Pi();
1917 t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN);
1919 Int_t first = t.GetFirstPoint();
1920 if (first<GetRowNumber(xt)+1) first = GetRowNumber(xt)+1;
1922 if (first<0) first=0;
1923 for (Int_t nr=first; nr<=rf; nr++) {
1924 // if ( (TMath::Abs(t.GetSnp())>0.95)) break;//patch 28 fev 06
1925 if (t.GetKinkIndexes()[0]<0){
1926 for (Int_t i=0;i<3;i++){
1927 Int_t index = t.GetKinkIndexes()[i];
1928 if (index==0) break;
1929 if (index>0) continue;
1930 index = TMath::Abs(index);
1931 AliKink * kink = (AliKink*)fEvent->GetKink(index-1);
1933 printf("PROBLEM\n");
1936 Int_t kinkrow = kink->GetTPCRow0()-2-Int_t(0.5/(0.05+kink->GetAngle(2)));
1938 AliExternalTrackParam paramm(t);
1939 kink->SetMother(paramm);
1940 kink->SetStatus(2,1);
1947 if (nr<fInnerSec->GetNRows())
1948 fSectors = fInnerSec;
1950 fSectors = fOuterSec;
1961 Float_t AliTPCtrackerMI::OverlapFactor(AliTPCseed * s1, AliTPCseed * s2, Int_t &sum1, Int_t & sum2)
1969 Float_t dz2 =(s1->GetZ() - s2->GetZ());
1972 Float_t dy2 =TMath::Abs((s1->GetY() - s2->GetY()));
1974 Float_t distance = TMath::Sqrt(dz2+dy2);
1975 if (distance>4.) return 0; // if there are far away - not overlap - to reduce combinatorics
1978 Int_t firstpoint = TMath::Min(s1->GetFirstPoint(),s2->GetFirstPoint());
1979 Int_t lastpoint = TMath::Max(s1->GetLastPoint(),s2->GetLastPoint());
1984 if (firstpoint>lastpoint) {
1985 firstpoint =lastpoint;
1990 for (Int_t i=firstpoint-1;i<lastpoint+1;i++){
1991 if (s1->GetClusterIndex2(i)>0) sum1++;
1992 if (s2->GetClusterIndex2(i)>0) sum2++;
1993 if (s1->GetClusterIndex2(i)==s2->GetClusterIndex2(i) && s1->GetClusterIndex2(i)>0) {
1997 if (sum<5) return 0;
1999 Float_t summin = TMath::Min(sum1+1,sum2+1);
2000 Float_t ratio = (sum+1)/Float_t(summin);
2004 void AliTPCtrackerMI::SignShared(AliTPCseed * s1, AliTPCseed * s2)
2008 Float_t thetaCut = 0.2;//+10.*TMath::Sqrt(s1->GetSigmaTglZ()+ s2->GetSigmaTglZ());
2009 if (TMath::Abs(s1->GetTgl()-s2->GetTgl())>thetaCut) return;
2010 Float_t minCl = TMath::Min(s1->GetNumberOfClusters(),s2->GetNumberOfClusters());
2011 Int_t cutN0 = TMath::Max(5,TMath::Nint(0.1*minCl));
2016 //Int_t firstpoint = TMath::Max(s1->GetFirstPoint(),s2->GetFirstPoint());
2017 //Int_t lastpoint = TMath::Min(s1->GetLastPoint(),s2->GetLastPoint());
2018 Int_t firstpoint = 0;
2019 Int_t lastpoint = 160;
2021 // if (firstpoint>=lastpoint-5) return;;
2023 for (Int_t i=firstpoint;i<lastpoint;i++){
2024 // if ( (s1->GetClusterIndex2(i)&0xFFFF8FFF)==(s2->GetClusterIndex2(i)&0xFFFF8FFF) && s1->GetClusterIndex2(i)>0) {
2025 if ( (s1->GetClusterIndex2(i))==(s2->GetClusterIndex2(i)) && s1->GetClusterIndex2(i)>0) {
2027 s1->SetSharedMapBit(i, kTRUE);
2028 s2->SetSharedMapBit(i, kTRUE);
2030 if (s1->GetClusterIndex2(i)>0)
2031 s1->SetClusterMapBit(i, kTRUE);
2033 if (sumshared>cutN0){
2036 for (Int_t i=firstpoint;i<lastpoint;i++){
2037 // if ( (s1->GetClusterIndex2(i)&0xFFFF8FFF)==(s2->GetClusterIndex2(i)&0xFFFF8FFF) && s1->GetClusterIndex2(i)>0) {
2038 if ( (s1->GetClusterIndex2(i))==(s2->GetClusterIndex2(i)) && s1->GetClusterIndex2(i)>0) {
2039 AliTPCTrackerPoint *p1 = s1->GetTrackPoint(i);
2040 AliTPCTrackerPoint *p2 = s2->GetTrackPoint(i);;
2041 if (s1->IsActive()&&s2->IsActive()){
2042 p1->SetShared(kTRUE);
2043 p2->SetShared(kTRUE);
2049 if (sumshared>cutN0){
2050 for (Int_t i=0;i<4;i++){
2051 if (s1->GetOverlapLabel(3*i)==0){
2052 s1->SetOverlapLabel(3*i, s2->GetLabel());
2053 s1->SetOverlapLabel(3*i+1,sumshared);
2054 s1->SetOverlapLabel(3*i+2,s2->GetUniqueID());
2058 for (Int_t i=0;i<4;i++){
2059 if (s2->GetOverlapLabel(3*i)==0){
2060 s2->SetOverlapLabel(3*i, s1->GetLabel());
2061 s2->SetOverlapLabel(3*i+1,sumshared);
2062 s2->SetOverlapLabel(3*i+2,s1->GetUniqueID());
2069 void AliTPCtrackerMI::SignShared(TObjArray * arr)
2072 //sort trackss according sectors
2074 for (Int_t i=0; i<arr->GetEntriesFast(); i++) {
2075 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2077 //if (pt) RotateToLocal(pt);
2081 arr->Sort(); // sorting according relative sectors
2082 arr->Expand(arr->GetEntries());
2085 Int_t nseed=arr->GetEntriesFast();
2086 for (Int_t i=0; i<nseed; i++) {
2087 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2089 for (Int_t j=0;j<=12;j++){
2090 pt->SetOverlapLabel(j,0);
2093 for (Int_t i=0; i<nseed; i++) {
2094 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2096 if (pt->GetRemoval()>10) continue;
2097 for (Int_t j=i+1; j<nseed; j++){
2098 AliTPCseed *pt2=(AliTPCseed*)arr->UncheckedAt(j);
2099 if (TMath::Abs(pt->GetRelativeSector()-pt2->GetRelativeSector())>1) continue;
2101 if (pt2->GetRemoval()<=10) {
2102 //if ( TMath::Abs(pt->GetRelativeSector()-pt2->GetRelativeSector())>0) break;
2110 void AliTPCtrackerMI::SortTracks(TObjArray * arr, Int_t mode) const
2113 //sort tracks in array according mode criteria
2114 Int_t nseed = arr->GetEntriesFast();
2115 for (Int_t i=0; i<nseed; i++) {
2116 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2127 void AliTPCtrackerMI::RemoveUsed2(TObjArray * arr, Float_t factor1, Float_t factor2, Int_t minimal)
2130 // Loop over all tracks and remove overlaped tracks (with lower quality)
2132 // 1. Unsign clusters
2133 // 2. Sort tracks according quality
2134 // Quality is defined by the number of cluster between first and last points
2136 // 3. Loop over tracks - decreasing quality order
2137 // a.) remove - If the fraction of shared cluster less than factor (1- n or 2)
2138 // b.) remove - If the minimal number of clusters less than minimal and not ITS
2139 // c.) if track accepted - sign clusters
2141 //Called in - AliTPCtrackerMI::Clusters2Tracks()
2142 // - AliTPCtrackerMI::PropagateBack()
2143 // - AliTPCtrackerMI::RefitInward()
2149 Int_t nseed = arr->GetEntriesFast();
2150 Float_t * quality = new Float_t[nseed];
2151 Int_t * indexes = new Int_t[nseed];
2155 for (Int_t i=0; i<nseed; i++) {
2156 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2161 pt->UpdatePoints(); //select first last max dens points
2162 Float_t * points = pt->GetPoints();
2163 if (points[3]<0.8) quality[i] =-1;
2164 quality[i] = (points[2]-points[0])+pt->GetNumberOfClusters();
2165 //prefer high momenta tracks if overlaps
2166 quality[i] *= TMath::Sqrt(TMath::Abs(pt->Pt())+0.5);
2168 TMath::Sort(nseed,quality,indexes);
2171 for (Int_t itrack=0; itrack<nseed; itrack++) {
2172 Int_t trackindex = indexes[itrack];
2173 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(trackindex);
2176 if (quality[trackindex]<0){
2178 delete arr->RemoveAt(trackindex);
2181 arr->RemoveAt(trackindex);
2187 Int_t first = Int_t(pt->GetPoints()[0]);
2188 Int_t last = Int_t(pt->GetPoints()[2]);
2189 Double_t factor = (pt->GetBConstrain()) ? factor1: factor2;
2191 Int_t found,foundable,shared;
2192 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
2193 // pt->GetClusterStatistic(0,160, found, foundable,shared,kFALSE);
2194 Bool_t itsgold =kFALSE;
2197 if (pt->GetESD()->GetITSclusters(dummy)>4) itsgold= kTRUE;
2201 if (Float_t(shared+1)/Float_t(found+1)>factor){
2202 if (pt->GetKinkIndexes()[0]!=0) continue; //don't remove tracks - part of the kinks
2203 delete arr->RemoveAt(trackindex);
2206 if (pt->GetNumberOfClusters()<50&&(found-0.5*shared)<minimal){ //remove short tracks
2207 if (pt->GetKinkIndexes()[0]!=0) continue; //don't remove tracks - part of the kinks
2208 delete arr->RemoveAt(trackindex);
2214 //if (sharedfactor>0.4) continue;
2215 if (pt->GetKinkIndexes()[0]>0) continue;
2216 //Remove tracks with undefined properties - seems
2217 if (pt->GetSigmaY2()<kAlmost0) continue; // ? what is the origin ?
2219 for (Int_t i=first; i<last; i++) {
2220 Int_t index=pt->GetClusterIndex2(i);
2221 // if (index<0 || index&0x8000 ) continue;
2222 if (index<0 || index&0x8000 ) continue;
2223 AliTPCclusterMI *c= pt->GetClusterPointer(i);
2230 Info("RemoveUsed","\n*****\nNumber of good tracks after shared removal\t%d\n",fNtracks);
2236 void AliTPCtrackerMI::UnsignClusters()
2239 // loop over all clusters and unsign them
2242 for (Int_t sec=0;sec<fkNIS;sec++){
2243 for (Int_t row=0;row<fInnerSec->GetNRows();row++){
2244 AliTPCclusterMI *cl = fInnerSec[sec][row].GetClusters1();
2245 for (Int_t icl =0;icl< fInnerSec[sec][row].GetN1();icl++)
2246 // if (cl[icl].IsUsed(10))
2248 cl = fInnerSec[sec][row].GetClusters2();
2249 for (Int_t icl =0;icl< fInnerSec[sec][row].GetN2();icl++)
2250 //if (cl[icl].IsUsed(10))
2255 for (Int_t sec=0;sec<fkNOS;sec++){
2256 for (Int_t row=0;row<fOuterSec->GetNRows();row++){
2257 AliTPCclusterMI *cl = fOuterSec[sec][row].GetClusters1();
2258 for (Int_t icl =0;icl< fOuterSec[sec][row].GetN1();icl++)
2259 //if (cl[icl].IsUsed(10))
2261 cl = fOuterSec[sec][row].GetClusters2();
2262 for (Int_t icl =0;icl< fOuterSec[sec][row].GetN2();icl++)
2263 //if (cl[icl].IsUsed(10))
2272 void AliTPCtrackerMI::SignClusters(TObjArray * arr, Float_t fnumber, Float_t fdensity)
2275 //sign clusters to be "used"
2277 // snumber and sdensity sign number of sigmas - bellow mean value to be accepted
2278 // loop over "primaries"
2292 Int_t nseed = arr->GetEntriesFast();
2293 for (Int_t i=0; i<nseed; i++) {
2294 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2298 if (!(pt->IsActive())) continue;
2299 Float_t dens = pt->GetNumberOfClusters()/Float_t(pt->GetNFoundable());
2300 if ( (dens>0.7) && (pt->GetNumberOfClusters()>70)){
2302 sumdens2+= dens*dens;
2303 sumn += pt->GetNumberOfClusters();
2304 sumn2 += pt->GetNumberOfClusters()*pt->GetNumberOfClusters();
2305 Float_t chi2 = pt->GetChi2()/pt->GetNumberOfClusters();
2308 sumchi2 +=chi2*chi2;
2313 Float_t mdensity = 0.9;
2314 Float_t meann = 130;
2315 Float_t meanchi = 1;
2316 Float_t sdensity = 0.1;
2317 Float_t smeann = 10;
2318 Float_t smeanchi =0.4;
2322 mdensity = sumdens/sum;
2324 meanchi = sumchi/sum;
2326 sdensity = sumdens2/sum-mdensity*mdensity;
2328 sdensity = TMath::Sqrt(sdensity);
2332 smeann = sumn2/sum-meann*meann;
2334 smeann = TMath::Sqrt(smeann);
2338 smeanchi = sumchi2/sum - meanchi*meanchi;
2340 smeanchi = TMath::Sqrt(smeanchi);
2346 //REMOVE SHORT DELTAS or tracks going out of sensitive volume of TPC
2348 for (Int_t i=0; i<nseed; i++) {
2349 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2353 if (pt->GetBSigned()) continue;
2354 if (pt->GetBConstrain()) continue;
2355 //if (!(pt->IsActive())) continue;
2357 Int_t found,foundable,shared;
2358 pt->GetClusterStatistic(0,160,found, foundable,shared);
2359 if (shared/float(found)>0.3) {
2360 if (shared/float(found)>0.9 ){
2361 //delete arr->RemoveAt(i);
2366 Bool_t isok =kFALSE;
2367 if ( (pt->GetNShared()/pt->GetNumberOfClusters()<0.5) &&pt->GetNumberOfClusters()>60)
2369 if ((TMath::Abs(1/pt->GetC())<100.) && (pt->GetNShared()/pt->GetNumberOfClusters()<0.7))
2371 if (TMath::Abs(pt->GetZ()/pt->GetX())>1.1)
2373 if ( (TMath::Abs(pt->GetSnp()>0.7) && pt->GetD(0,0)>60.))
2377 for (Int_t i=0; i<160; i++) {
2378 Int_t index=pt->GetClusterIndex2(i);
2379 if (index<0) continue;
2380 AliTPCclusterMI *c= pt->GetClusterPointer(i);
2382 //if (!(c->IsUsed(10))) c->Use();
2389 Double_t maxchi = meanchi+2.*smeanchi;
2391 for (Int_t i=0; i<nseed; i++) {
2392 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2396 //if (!(pt->IsActive())) continue;
2397 if (pt->GetBSigned()) continue;
2398 Double_t chi = pt->GetChi2()/pt->GetNumberOfClusters();
2399 if (chi>maxchi) continue;
2402 Float_t dens = pt->GetNumberOfClusters()/Float_t(pt->GetNFoundable());
2404 //sign only tracks with enoug big density at the beginning
2406 if ((pt->GetDensityFirst(40)<0.75) && pt->GetNumberOfClusters()<meann) continue;
2409 Double_t mindens = TMath::Max(double(mdensity-sdensity*fdensity*bfactor),0.65);
2410 Double_t minn = TMath::Max(Int_t(meann-fnumber*smeann*bfactor),50);
2412 // if (pt->fBConstrain) mindens = TMath::Max(mdensity-sdensity*fdensity*bfactor,0.65);
2413 if ( (pt->GetRemoval()==10) && (pt->GetSnp()>0.8)&&(dens>mindens))
2416 if ((dens>mindens && pt->GetNumberOfClusters()>minn) && chi<maxchi ){
2417 //Int_t noc=pt->GetNumberOfClusters();
2418 pt->SetBSigned(kTRUE);
2419 for (Int_t i=0; i<160; i++) {
2421 Int_t index=pt->GetClusterIndex2(i);
2422 if (index<0) continue;
2423 AliTPCclusterMI *c= pt->GetClusterPointer(i);
2425 // if (!(c->IsUsed(10))) c->Use();
2430 // gLastCheck = nseed;
2438 void AliTPCtrackerMI::StopNotActive(TObjArray * arr, Int_t row0, Float_t th0, Float_t th1, Float_t th2) const
2440 // stop not active tracks
2441 // take th1 as threshold for number of founded to number of foundable on last 10 active rows
2442 // take th2 as threshold for number of founded to number of foundable on last 20 active rows
2443 Int_t nseed = arr->GetEntriesFast();
2445 for (Int_t i=0; i<nseed; i++) {
2446 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2450 if (!(pt->IsActive())) continue;
2451 StopNotActive(pt,row0,th0, th1,th2);
2457 void AliTPCtrackerMI::StopNotActive(AliTPCseed * seed, Int_t row0, Float_t th0, Float_t th1,
2460 // stop not active tracks
2461 // take th1 as threshold for number of founded to number of foundable on last 10 active rows
2462 // take th2 as threshold for number of founded to number of foundable on last 20 active rows
2465 Int_t foundable = 0;
2466 Int_t maxindex = seed->GetLastPoint(); //last foundable row
2467 if (seed->GetNFoundable()*th0 > seed->GetNumberOfClusters()) {
2468 seed->Desactivate(10) ;
2472 for (Int_t i=row0; i<maxindex; i++){
2473 Int_t index = seed->GetClusterIndex2(i);
2474 if (index!=-1) foundable++;
2476 if (foundable<=30) sumgood1++;
2477 if (foundable<=50) {
2484 if (foundable>=30.){
2485 if (sumgood1<(th1*30.)) seed->Desactivate(10);
2488 if (sumgood2<(th2*50.)) seed->Desactivate(10);
2492 Int_t AliTPCtrackerMI::RefitInward(AliESDEvent *event)
2495 // back propagation of ESD tracks
2498 const Int_t kMaxFriendTracks=2000;
2502 //PrepareForProlongation(fSeeds,1);
2503 PropagateForward2(fSeeds);
2504 RemoveUsed2(fSeeds,0.4,0.4,20);
2506 TObjArray arraySeed(fSeeds->GetEntries());
2507 for (Int_t i=0;i<fSeeds->GetEntries();i++) {
2508 arraySeed.AddAt(fSeeds->At(i),i);
2510 SignShared(&arraySeed);
2511 // FindCurling(fSeeds, event,2); // find multi found tracks
2512 FindSplitted(fSeeds, event,2); // find multi found tracks
2515 Int_t nseed = fSeeds->GetEntriesFast();
2516 for (Int_t i=0;i<nseed;i++){
2517 AliTPCseed * seed = (AliTPCseed*) fSeeds->UncheckedAt(i);
2518 if (!seed) continue;
2519 if (seed->GetKinkIndex(0)>0) UpdateKinkQualityD(seed); // update quality informations for kinks
2521 seed->PropagateTo(fParam->GetInnerRadiusLow());
2522 seed->UpdatePoints();
2524 AliESDtrack *esd=event->GetTrack(i);
2525 seed->CookdEdx(0.02,0.6);
2526 CookLabel(seed,0.1); //For comparison only
2527 esd->SetTPCClusterMap(seed->GetClusterMap());
2528 esd->SetTPCSharedMap(seed->GetSharedMap());
2530 if (AliTPCReconstructor::StreamLevel()>1 && seed!=0&&esd!=0) {
2531 TTreeSRedirector &cstream = *fDebugStreamer;
2538 if (seed->GetNumberOfClusters()>15){
2539 esd->UpdateTrackParams(seed,AliESDtrack::kTPCrefit);
2540 esd->SetTPCPoints(seed->GetPoints());
2541 esd->SetTPCPointsF(seed->GetNFoundable());
2542 Int_t ndedx = seed->GetNCDEDX(0)+seed->GetNCDEDX(1)+seed->GetNCDEDX(2)+seed->GetNCDEDX(3);
2543 Float_t sdedx = (seed->GetSDEDX(0)+seed->GetSDEDX(1)+seed->GetSDEDX(2)+seed->GetSDEDX(3))*0.25;
2544 Float_t dedx = seed->GetdEdx();
2545 esd->SetTPCsignal(dedx, sdedx, ndedx);
2547 // add seed to the esd track in Calib level
2549 Bool_t storeFriend = gRandom->Rndm()<(kMaxFriendTracks)/Float_t(nseed);
2550 if (AliTPCReconstructor::StreamLevel()>0 &&storeFriend){
2551 AliTPCseed * seedCopy = new AliTPCseed(*seed, kTRUE);
2552 esd->AddCalibObject(seedCopy);
2557 //printf("problem\n");
2560 //FindKinks(fSeeds,event);
2561 Info("RefitInward","Number of refitted tracks %d",ntracks);
2567 Int_t AliTPCtrackerMI::PropagateBack(AliESDEvent *event)
2570 // back propagation of ESD tracks
2576 PropagateBack(fSeeds);
2577 RemoveUsed2(fSeeds,0.4,0.4,20);
2578 //FindCurling(fSeeds, fEvent,1);
2579 FindSplitted(fSeeds, event,1); // find multi found tracks
2582 Int_t nseed = fSeeds->GetEntriesFast();
2584 for (Int_t i=0;i<nseed;i++){
2585 AliTPCseed * seed = (AliTPCseed*) fSeeds->UncheckedAt(i);
2586 if (!seed) continue;
2587 if (seed->GetKinkIndex(0)<0) UpdateKinkQualityM(seed); // update quality informations for kinks
2588 seed->UpdatePoints();
2589 AliESDtrack *esd=event->GetTrack(i);
2590 seed->CookdEdx(0.02,0.6);
2591 CookLabel(seed,0.1); //For comparison only
2592 if (seed->GetNumberOfClusters()>15){
2593 esd->UpdateTrackParams(seed,AliESDtrack::kTPCout);
2594 esd->SetTPCPoints(seed->GetPoints());
2595 esd->SetTPCPointsF(seed->GetNFoundable());
2596 Int_t ndedx = seed->GetNCDEDX(0)+seed->GetNCDEDX(1)+seed->GetNCDEDX(2)+seed->GetNCDEDX(3);
2597 Float_t sdedx = (seed->GetSDEDX(0)+seed->GetSDEDX(1)+seed->GetSDEDX(2)+seed->GetSDEDX(3))*0.25;
2598 Float_t dedx = seed->GetdEdx();
2599 esd->SetTPCsignal(dedx, sdedx, ndedx);
2601 Int_t eventnumber = event->GetEventNumberInFile();// patch 28 fev 06
2602 // This is most likely NOT the event number you'd like to use. It has nothing to do with the 'real' event number
2603 if (AliTPCReconstructor::StreamLevel()>1) {
2604 (*fDebugStreamer)<<"Cback"<<
2606 "EventNrInFile="<<eventnumber<<
2607 "\n"; // patch 28 fev 06
2611 //FindKinks(fSeeds,event);
2612 Info("PropagateBack","Number of back propagated tracks %d",ntracks);
2619 void AliTPCtrackerMI::DeleteSeeds()
2624 Int_t nseed = fSeeds->GetEntriesFast();
2625 for (Int_t i=0;i<nseed;i++){
2626 AliTPCseed * seed = (AliTPCseed*)fSeeds->At(i);
2627 if (seed) delete fSeeds->RemoveAt(i);
2634 void AliTPCtrackerMI::ReadSeeds(AliESDEvent *event, Int_t direction)
2637 //read seeds from the event
2639 Int_t nentr=event->GetNumberOfTracks();
2641 Info("ReadSeeds", "Number of ESD tracks: %d\n", nentr);
2646 fSeeds = new TObjArray(nentr);
2650 for (Int_t i=0; i<nentr; i++) {
2651 AliESDtrack *esd=event->GetTrack(i);
2652 ULong_t status=esd->GetStatus();
2653 if (!(status&AliESDtrack::kTPCin)) continue;
2654 AliTPCtrack t(*esd);
2655 t.SetNumberOfClusters(0);
2656 // AliTPCseed *seed = new AliTPCseed(t,t.GetAlpha());
2657 AliTPCseed *seed = new AliTPCseed(t/*,t.GetAlpha()*/);
2658 seed->SetUniqueID(esd->GetID());
2659 for (Int_t ikink=0;ikink<3;ikink++) {
2660 Int_t index = esd->GetKinkIndex(ikink);
2661 seed->GetKinkIndexes()[ikink] = index;
2662 if (index==0) continue;
2663 index = TMath::Abs(index);
2664 AliESDkink * kink = fEvent->GetKink(index-1);
2665 if (kink&&esd->GetKinkIndex(ikink)<0){
2666 if ((status & AliESDtrack::kTRDrefit) != 0) kink->SetStatus(1,2);
2667 if ((status & AliESDtrack::kITSout) != 0) kink->SetStatus(1,0);
2669 if (kink&&esd->GetKinkIndex(ikink)>0){
2670 if ((status & AliESDtrack::kTRDrefit) != 0) kink->SetStatus(1,6);
2671 if ((status & AliESDtrack::kITSout) != 0) kink->SetStatus(1,4);
2675 if (((status&AliESDtrack::kITSout)==0)&&(direction==1)) seed->ResetCovariance(10.);
2676 if ( direction ==2 &&(status & AliESDtrack::kTRDrefit) == 0 ) seed->ResetCovariance(10.);
2677 if ( direction ==2 && ((status & AliESDtrack::kTPCout) == 0) ) {
2682 if ( direction ==2 &&(status & AliESDtrack::kTRDrefit) > 0 ) {
2683 Double_t par0[5],par1[5],alpha,x;
2684 esd->GetInnerExternalParameters(alpha,x,par0);
2685 esd->GetExternalParameters(x,par1);
2686 Double_t delta1 = TMath::Abs(par0[4]-par1[4])/(0.000000001+TMath::Abs(par0[4]+par1[4]));
2687 Double_t delta2 = TMath::Abs(par0[3]-par1[3]);
2689 if (esd->GetTRDncls()>0) trdchi2 = esd->GetTRDchi2()/esd->GetTRDncls();
2690 //reset covariance if suspicious
2691 if ( (delta1>0.1) || (delta2>0.006) ||trdchi2>7.)
2692 seed->ResetCovariance(10.);
2697 // rotate to the local coordinate system
2699 fSectors=fInnerSec; fN=fkNIS;
2700 Double_t alpha=seed->GetAlpha() - fSectors->GetAlphaShift();
2701 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
2702 if (alpha < 0. ) alpha += 2.*TMath::Pi();
2703 Int_t ns=Int_t(alpha/fSectors->GetAlpha())%fN;
2704 alpha =ns*fSectors->GetAlpha() + fSectors->GetAlphaShift();
2705 if (alpha<-TMath::Pi()) alpha += 2*TMath::Pi();
2706 if (alpha>=TMath::Pi()) alpha -= 2*TMath::Pi();
2707 alpha-=seed->GetAlpha();
2708 if (!seed->Rotate(alpha)) {
2714 if (esd->GetKinkIndex(0)<=0){
2715 for (Int_t irow=0;irow<160;irow++){
2716 Int_t index = seed->GetClusterIndex2(irow);
2719 AliTPCclusterMI * cl = GetClusterMI(index);
2720 seed->SetClusterPointer(irow,cl);
2722 if ((index & 0x8000)==0){
2723 cl->Use(10); // accepted cluster
2725 cl->Use(6); // close cluster not accepted
2728 Info("ReadSeeds","Not found cluster");
2733 fSeeds->AddAt(seed,i);
2739 //_____________________________________________________________________________
2740 void AliTPCtrackerMI::MakeSeeds3(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, Float_t cuts[4],
2741 Float_t deltay, Int_t ddsec) {
2742 //-----------------------------------------------------------------
2743 // This function creates track seeds.
2744 // SEEDING WITH VERTEX CONSTRAIN
2745 //-----------------------------------------------------------------
2746 // cuts[0] - fP4 cut
2747 // cuts[1] - tan(phi) cut
2748 // cuts[2] - zvertex cut
2749 // cuts[3] - fP3 cut
2757 Double_t x[5], c[15];
2758 // Int_t di = i1-i2;
2760 AliTPCseed * seed = new AliTPCseed();
2761 Double_t alpha=fSectors->GetAlpha(), shift=fSectors->GetAlphaShift();
2762 Double_t cs=cos(alpha), sn=sin(alpha);
2764 // Double_t x1 =fOuterSec->GetX(i1);
2765 //Double_t xx2=fOuterSec->GetX(i2);
2767 Double_t x1 =GetXrow(i1);
2768 Double_t xx2=GetXrow(i2);
2770 Double_t x3=GetX(), y3=GetY(), z3=GetZ();
2772 Int_t imiddle = (i2+i1)/2; //middle pad row index
2773 Double_t xm = GetXrow(imiddle); // radius of middle pad-row
2774 const AliTPCtrackerRow& krm=GetRow(sec,imiddle); //middle pad -row
2778 const AliTPCtrackerRow& kr1=GetRow(ns,i1);
2779 Double_t ymax = GetMaxY(i1)-kr1.GetDeadZone()-1.5;
2780 Double_t ymaxm = GetMaxY(imiddle)-kr1.GetDeadZone()-1.5;
2783 // change cut on curvature if it can't reach this layer
2784 // maximal curvature set to reach it
2785 Double_t dvertexmax = TMath::Sqrt((x1-x3)*(x1-x3)+(ymax+5-y3)*(ymax+5-y3));
2786 if (dvertexmax*0.5*cuts[0]>0.85){
2787 cuts[0] = 0.85/(dvertexmax*0.5+1.);
2789 Double_t r2min = 1/(cuts[0]*cuts[0]); //minimal square of radius given by cut
2792 if (deltay>0) ddsec = 0;
2793 // loop over clusters
2794 for (Int_t is=0; is < kr1; is++) {
2796 if (kr1[is]->IsUsed(10)) continue;
2797 Double_t y1=kr1[is]->GetY(), z1=kr1[is]->GetZ();
2798 //if (TMath::Abs(y1)>ymax) continue;
2800 if (deltay>0 && TMath::Abs(ymax-TMath::Abs(y1))> deltay ) continue; // seed only at the edge
2802 // find possible directions
2803 Float_t anglez = (z1-z3)/(x1-x3);
2804 Float_t extraz = z1 - anglez*(x1-xx2); // extrapolated z
2807 //find rotation angles relative to line given by vertex and point 1
2808 Double_t dvertex2 = (x1-x3)*(x1-x3)+(y1-y3)*(y1-y3);
2809 Double_t dvertex = TMath::Sqrt(dvertex2);
2810 Double_t angle13 = TMath::ATan((y1-y3)/(x1-x3));
2811 Double_t cs13 = cos(-angle13), sn13 = sin(-angle13);
2814 // loop over 2 sectors
2820 Double_t dddz1=0; // direction of delta inclination in z axis
2827 for (Int_t dsec = dsec1; dsec<=dsec2;dsec++){
2828 Int_t sec2 = sec + dsec;
2830 // AliTPCtrackerRow& kr2 = fOuterSec[(sec2+fkNOS)%fkNOS][i2];
2831 //AliTPCtrackerRow& kr2m = fOuterSec[(sec2+fkNOS)%fkNOS][imiddle];
2832 AliTPCtrackerRow& kr2 = GetRow((sec2+fkNOS)%fkNOS,i2);
2833 AliTPCtrackerRow& kr2m = GetRow((sec2+fkNOS)%fkNOS,imiddle);
2834 Int_t index1 = TMath::Max(kr2.Find(extraz-0.6-dddz1*TMath::Abs(z1)*0.05)-1,0);
2835 Int_t index2 = TMath::Min(kr2.Find(extraz+0.6+dddz2*TMath::Abs(z1)*0.05)+1,kr2);
2837 // rotation angles to p1-p3
2838 Double_t cs13r = cos(-angle13+dsec*alpha)/dvertex, sn13r = sin(-angle13+dsec*alpha)/dvertex;
2839 Double_t x2, y2, z2;
2841 // Double_t dymax = maxangle*TMath::Abs(x1-xx2);
2844 Double_t dxx0 = (xx2-x3)*cs13r;
2845 Double_t dyy0 = (xx2-x3)*sn13r;
2846 for (Int_t js=index1; js < index2; js++) {
2847 const AliTPCclusterMI *kcl = kr2[js];
2848 if (kcl->IsUsed(10)) continue;
2850 //calcutate parameters
2852 Double_t yy0 = dyy0 +(kcl->GetY()-y3)*cs13r;
2854 if (TMath::Abs(yy0)<0.000001) continue;
2855 Double_t xx0 = dxx0 -(kcl->GetY()-y3)*sn13r;
2856 Double_t y0 = 0.5*(xx0*xx0+yy0*yy0-xx0)/yy0;
2857 Double_t r02 = (0.25+y0*y0)*dvertex2;
2858 //curvature (radius) cut
2859 if (r02<r2min) continue;
2863 Double_t c0 = 1/TMath::Sqrt(r02);
2867 //Double_t dfi0 = 2.*TMath::ASin(dvertex*c0*0.5);
2868 //Double_t dfi1 = 2.*TMath::ASin(TMath::Sqrt(yy0*yy0+(1-xx0)*(1-xx0))*dvertex*c0*0.5);
2869 Double_t dfi0 = 2.*AliTPCFastMath::FastAsin(dvertex*c0*0.5);
2870 Double_t dfi1 = 2.*AliTPCFastMath::FastAsin(TMath::Sqrt(yy0*yy0+(1-xx0)*(1-xx0))*dvertex*c0*0.5);
2873 Double_t z0 = kcl->GetZ();
2874 Double_t zzzz2 = z1-(z1-z3)*dfi1/dfi0;
2875 if (TMath::Abs(zzzz2-z0)>0.5) continue;
2878 Double_t dip = (z1-z0)*c0/dfi1;
2879 Double_t x0 = (0.5*cs13+y0*sn13)*dvertex*c0;
2890 x2= xx2*cs-y2*sn*dsec;
2891 y2=+xx2*sn*dsec+y2*cs;
2901 // do we have cluster at the middle ?
2903 GetProlongation(x1,xm,x,ym,zm);
2905 AliTPCclusterMI * cm=0;
2906 if (TMath::Abs(ym)-ymaxm<0){
2907 cm = krm.FindNearest2(ym,zm,1.0,0.6,dummy);
2908 if ((!cm) || (cm->IsUsed(10))) {
2913 // rotate y1 to system 0
2914 // get state vector in rotated system
2915 Double_t yr1 = (-0.5*sn13+y0*cs13)*dvertex*c0;
2916 Double_t xr2 = x0*cs+yr1*sn*dsec;
2917 Double_t xr[5]={kcl->GetY(),kcl->GetZ(), xr2, dip, c0};
2919 GetProlongation(xx2,xm,xr,ym,zm);
2920 if (TMath::Abs(ym)-ymaxm<0){
2921 cm = kr2m.FindNearest2(ym,zm,1.0,0.6,dummy);
2922 if ((!cm) || (cm->IsUsed(10))) {
2932 dym = ym - cm->GetY();
2933 dzm = zm - cm->GetZ();
2940 Double_t sy1=kr1[is]->GetSigmaY2()*2., sz1=kr1[is]->GetSigmaZ2()*2.;
2941 Double_t sy2=kcl->GetSigmaY2()*2., sz2=kcl->GetSigmaZ2()*2.;
2942 //Double_t sy3=400*3./12., sy=0.1, sz=0.1;
2943 Double_t sy3=25000*x[4]*x[4]+0.1, sy=0.1, sz=0.1;
2944 //Double_t sy3=25000*x[4]*x[4]*60+0.5, sy=0.1, sz=0.1;
2946 Double_t f40=(F1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
2947 Double_t f42=(F1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
2948 Double_t f43=(F1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
2949 Double_t f20=(F2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
2950 Double_t f22=(F2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
2951 Double_t f23=(F2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
2953 Double_t f30=(F3(x1,y1+sy,x2,y2,z1,z2)-x[3])/sy;
2954 Double_t f31=(F3(x1,y1,x2,y2,z1+sz,z2)-x[3])/sz;
2955 Double_t f32=(F3(x1,y1,x2,y2+sy,z1,z2)-x[3])/sy;
2956 Double_t f34=(F3(x1,y1,x2,y2,z1,z2+sz)-x[3])/sz;
2960 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
2961 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
2962 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
2963 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
2964 c[13]=f30*sy1*f40+f32*sy2*f42;
2965 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
2967 // if (!BuildSeed(kr1[is],kcl,0,x1,x2,x3,x,c)) continue;
2969 UInt_t index=kr1.GetIndex(is);
2970 seed->~AliTPCseed(); // this does not set the pointer to 0...
2971 AliTPCseed *track=new(seed) AliTPCseed(x1, ns*alpha+shift, x, c, index);
2973 track->SetIsSeeding(kTRUE);
2974 track->SetSeed1(i1);
2975 track->SetSeed2(i2);
2976 track->SetSeedType(3);
2980 FollowProlongation(*track, (i1+i2)/2,1);
2981 Int_t foundable,found,shared;
2982 track->GetClusterStatistic((i1+i2)/2,i1, found, foundable, shared, kTRUE);
2983 if ((found<0.55*foundable) || shared>0.5*found || (track->GetSigmaY2()+track->GetSigmaZ2())>0.5){
2985 seed->~AliTPCseed();
2991 FollowProlongation(*track, i2,1);
2995 track->SetBConstrain(1);
2996 // track->fLastPoint = i1+fInnerSec->GetNRows(); // first cluster in track position
2997 track->SetLastPoint(i1); // first cluster in track position
2998 track->SetFirstPoint(track->GetLastPoint());
3000 if (track->GetNumberOfClusters()<(i1-i2)*0.5 ||
3001 track->GetNumberOfClusters() < track->GetNFoundable()*0.6 ||
3002 track->GetNShared()>0.4*track->GetNumberOfClusters() ) {
3004 seed->~AliTPCseed();
3008 // Z VERTEX CONDITION
3009 Double_t zv, bz=GetBz();
3010 if ( !track->GetZAt(0.,bz,zv) ) continue;
3011 if (TMath::Abs(zv-z3)>cuts[2]) {
3012 FollowProlongation(*track, TMath::Max(i2-20,0));
3013 if ( !track->GetZAt(0.,bz,zv) ) continue;
3014 if (TMath::Abs(zv-z3)>cuts[2]){
3015 FollowProlongation(*track, TMath::Max(i2-40,0));
3016 if ( !track->GetZAt(0.,bz,zv) ) continue;
3017 if (TMath::Abs(zv-z3)>cuts[2] &&(track->GetNumberOfClusters() > track->GetNFoundable()*0.7)){
3018 // make seed without constrain
3019 AliTPCseed * track2 = MakeSeed(track,0.2,0.5,1.);
3020 FollowProlongation(*track2, i2,1);
3021 track2->SetBConstrain(kFALSE);
3022 track2->SetSeedType(1);
3023 arr->AddLast(track2);
3025 seed->~AliTPCseed();
3030 seed->~AliTPCseed();
3037 track->SetSeedType(0);
3038 arr->AddLast(track);
3039 seed = new AliTPCseed;
3041 // don't consider other combinations
3042 if (track->GetNumberOfClusters() > track->GetNFoundable()*0.8)
3048 Info("MakeSeeds3","\nSeeding statistic:\t%d\t%d\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin,nout1,nout2);
3054 void AliTPCtrackerMI::MakeSeeds5(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, Float_t cuts[4],
3059 //-----------------------------------------------------------------
3060 // This function creates track seeds.
3061 //-----------------------------------------------------------------
3062 // cuts[0] - fP4 cut
3063 // cuts[1] - tan(phi) cut
3064 // cuts[2] - zvertex cut
3065 // cuts[3] - fP3 cut
3075 Double_t x[5], c[15];
3077 // make temporary seed
3078 AliTPCseed * seed = new AliTPCseed;
3079 Double_t alpha=fOuterSec->GetAlpha(), shift=fOuterSec->GetAlphaShift();
3080 // Double_t cs=cos(alpha), sn=sin(alpha);
3085 Double_t x1 = GetXrow(i1-1);
3086 const AliTPCtrackerRow& kr1=GetRow(sec,i1-1);
3087 Double_t y1max = GetMaxY(i1-1)-kr1.GetDeadZone()-1.5;
3089 Double_t x1p = GetXrow(i1);
3090 const AliTPCtrackerRow& kr1p=GetRow(sec,i1);
3092 Double_t x1m = GetXrow(i1-2);
3093 const AliTPCtrackerRow& kr1m=GetRow(sec,i1-2);
3096 //last 3 padrow for seeding
3097 AliTPCtrackerRow& kr3 = GetRow((sec+fkNOS)%fkNOS,i1-7);
3098 Double_t x3 = GetXrow(i1-7);
3099 // Double_t y3max= GetMaxY(i1-7)-kr3.fDeadZone-1.5;
3101 AliTPCtrackerRow& kr3p = GetRow((sec+fkNOS)%fkNOS,i1-6);
3102 Double_t x3p = GetXrow(i1-6);
3104 AliTPCtrackerRow& kr3m = GetRow((sec+fkNOS)%fkNOS,i1-8);
3105 Double_t x3m = GetXrow(i1-8);
3110 Int_t im = i1-4; //middle pad row index
3111 Double_t xm = GetXrow(im); // radius of middle pad-row
3112 const AliTPCtrackerRow& krm=GetRow(sec,im); //middle pad -row
3113 // Double_t ymmax = GetMaxY(im)-kr1.fDeadZone-1.5;
3116 Double_t deltax = x1-x3;
3117 Double_t dymax = deltax*cuts[1];
3118 Double_t dzmax = deltax*cuts[3];
3120 // loop over clusters
3121 for (Int_t is=0; is < kr1; is++) {
3123 if (kr1[is]->IsUsed(10)) continue;
3124 Double_t y1=kr1[is]->GetY(), z1=kr1[is]->GetZ();
3126 if (deltay>0 && TMath::Abs(y1max-TMath::Abs(y1))> deltay ) continue; // seed only at the edge
3128 Int_t index1 = TMath::Max(kr3.Find(z1-dzmax)-1,0);
3129 Int_t index2 = TMath::Min(kr3.Find(z1+dzmax)+1,kr3);
3135 for (Int_t js=index1; js < index2; js++) {
3136 const AliTPCclusterMI *kcl = kr3[js];
3137 if (kcl->IsUsed(10)) continue;
3139 // apply angular cuts
3140 if (TMath::Abs(y1-y3)>dymax) continue;
3143 if (TMath::Abs(z1-z3)>dzmax) continue;
3145 Double_t angley = (y1-y3)/(x1-x3);
3146 Double_t anglez = (z1-z3)/(x1-x3);
3148 Double_t erry = TMath::Abs(angley)*(x1-x1m)*0.5+0.5;
3149 Double_t errz = TMath::Abs(anglez)*(x1-x1m)*0.5+0.5;
3151 Double_t yyym = angley*(xm-x1)+y1;
3152 Double_t zzzm = anglez*(xm-x1)+z1;
3154 const AliTPCclusterMI *kcm = krm.FindNearest2(yyym,zzzm,erry,errz,index);
3156 if (kcm->IsUsed(10)) continue;
3158 erry = TMath::Abs(angley)*(x1-x1m)*0.4+0.5;
3159 errz = TMath::Abs(anglez)*(x1-x1m)*0.4+0.5;
3166 // look around first
3167 const AliTPCclusterMI *kc1m = kr1m.FindNearest2(angley*(x1m-x1)+y1,
3173 if (kc1m->IsUsed(10)) used++;
3175 const AliTPCclusterMI *kc1p = kr1p.FindNearest2(angley*(x1p-x1)+y1,
3181 if (kc1p->IsUsed(10)) used++;
3183 if (used>1) continue;
3184 if (found<1) continue;
3188 const AliTPCclusterMI *kc3m = kr3m.FindNearest2(angley*(x3m-x3)+y3,
3194 if (kc3m->IsUsed(10)) used++;
3198 const AliTPCclusterMI *kc3p = kr3p.FindNearest2(angley*(x3p-x3)+y3,
3204 if (kc3p->IsUsed(10)) used++;
3208 if (used>1) continue;
3209 if (found<3) continue;
3219 x[4]=F1(x1,y1,x2,y2,x3,y3);
3220 //if (TMath::Abs(x[4]) >= cuts[0]) continue;
3223 x[2]=F2(x1,y1,x2,y2,x3,y3);
3226 x[3]=F3n(x1,y1,x2,y2,z1,z2,x[4]);
3227 //if (TMath::Abs(x[3]) > cuts[3]) continue;
3231 Double_t sy1=0.1, sz1=0.1;
3232 Double_t sy2=0.1, sz2=0.1;
3233 Double_t sy3=0.1, sy=0.1, sz=0.1;
3235 Double_t f40=(F1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
3236 Double_t f42=(F1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
3237 Double_t f43=(F1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
3238 Double_t f20=(F2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
3239 Double_t f22=(F2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
3240 Double_t f23=(F2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
3242 Double_t f30=(F3(x1,y1+sy,x2,y2,z1,z2)-x[3])/sy;
3243 Double_t f31=(F3(x1,y1,x2,y2,z1+sz,z2)-x[3])/sz;
3244 Double_t f32=(F3(x1,y1,x2,y2+sy,z1,z2)-x[3])/sy;
3245 Double_t f34=(F3(x1,y1,x2,y2,z1,z2+sz)-x[3])/sz;
3249 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3250 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3251 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3252 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3253 c[13]=f30*sy1*f40+f32*sy2*f42;
3254 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3256 // if (!BuildSeed(kr1[is],kcl,0,x1,x2,x3,x,c)) continue;
3258 UInt_t index=kr1.GetIndex(is);
3259 seed->~AliTPCseed();
3260 AliTPCseed *track=new(seed) AliTPCseed(x1, sec*alpha+shift, x, c, index);
3262 track->SetIsSeeding(kTRUE);
3265 FollowProlongation(*track, i1-7,1);
3266 if (track->GetNumberOfClusters() < track->GetNFoundable()*0.75 ||
3267 track->GetNShared()>0.6*track->GetNumberOfClusters() || ( track->GetSigmaY2()+ track->GetSigmaZ2())>0.6){
3269 seed->~AliTPCseed();
3275 FollowProlongation(*track, i2,1);
3276 track->SetBConstrain(0);
3277 track->SetLastPoint(i1+fInnerSec->GetNRows()); // first cluster in track position
3278 track->SetFirstPoint(track->GetLastPoint());
3280 if (track->GetNumberOfClusters()<(i1-i2)*0.5 ||
3281 track->GetNumberOfClusters()<track->GetNFoundable()*0.7 ||
3282 track->GetNShared()>2. || track->GetChi2()/track->GetNumberOfClusters()>6 || ( track->GetSigmaY2()+ track->GetSigmaZ2())>0.5 ) {
3284 seed->~AliTPCseed();
3289 FollowProlongation(*track, TMath::Max(i2-10,0),1);
3290 AliTPCseed * track2 = MakeSeed(track,0.2,0.5,0.9);
3291 FollowProlongation(*track2, i2,1);
3292 track2->SetBConstrain(kFALSE);
3293 track2->SetSeedType(4);
3294 arr->AddLast(track2);
3296 seed->~AliTPCseed();
3300 //arr->AddLast(track);
3301 //seed = new AliTPCseed;
3307 Info("MakeSeeds5","\nSeeding statiistic:\t%d\t%d\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin,nout1,nout2,nout3);
3313 //_____________________________________________________________________________
3314 void AliTPCtrackerMI::MakeSeeds2(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, Float_t */*cuts[4]*/,
3315 Float_t deltay, Bool_t /*bconstrain*/) {
3316 //-----------------------------------------------------------------
3317 // This function creates track seeds - without vertex constraint
3318 //-----------------------------------------------------------------
3319 // cuts[0] - fP4 cut - not applied
3320 // cuts[1] - tan(phi) cut
3321 // cuts[2] - zvertex cut - not applied
3322 // cuts[3] - fP3 cut
3332 Double_t alpha=fOuterSec->GetAlpha(), shift=fOuterSec->GetAlphaShift();
3333 // Double_t cs=cos(alpha), sn=sin(alpha);
3334 Int_t row0 = (i1+i2)/2;
3335 Int_t drow = (i1-i2)/2;
3336 const AliTPCtrackerRow& kr0=fSectors[sec][row0];
3337 AliTPCtrackerRow * kr=0;
3339 AliTPCpolyTrack polytrack;
3340 Int_t nclusters=fSectors[sec][row0];
3341 AliTPCseed * seed = new AliTPCseed;
3346 for (Int_t is=0; is < nclusters; is++) { //LOOP over clusters
3348 Int_t nfoundable =0;
3349 for (Int_t iter =1; iter<2; iter++){ //iterations
3350 const AliTPCtrackerRow& krm=fSectors[sec][row0-iter];
3351 const AliTPCtrackerRow& krp=fSectors[sec][row0+iter];
3352 const AliTPCclusterMI * cl= kr0[is];
3354 if (cl->IsUsed(10)) {
3360 Double_t x = kr0.GetX();
3361 // Initialization of the polytrack
3366 Double_t y0= cl->GetY();
3367 Double_t z0= cl->GetZ();
3371 Double_t ymax = fSectors->GetMaxY(row0)-kr0.GetDeadZone()-1.5;
3372 if (deltay>0 && TMath::Abs(ymax-TMath::Abs(y0))> deltay ) continue; // seed only at the edge
3374 erry = (0.5)*cl->GetSigmaY2()/TMath::Sqrt(cl->GetQ())*6;
3375 errz = (0.5)*cl->GetSigmaZ2()/TMath::Sqrt(cl->GetQ())*6;
3376 polytrack.AddPoint(x,y0,z0,erry, errz);
3379 if (cl->IsUsed(10)) sumused++;
3382 Float_t roady = (5*TMath::Sqrt(cl->GetSigmaY2()+0.2)+1.)*iter;
3383 Float_t roadz = (5*TMath::Sqrt(cl->GetSigmaZ2()+0.2)+1.)*iter;
3386 AliTPCclusterMI * cl1 = krm.FindNearest(y0,z0,roady,roadz);
3387 if (cl1 && TMath::Abs(ymax-TMath::Abs(y0))) {
3388 erry = (0.5)*cl1->GetSigmaY2()/TMath::Sqrt(cl1->GetQ())*3;
3389 errz = (0.5)*cl1->GetSigmaZ2()/TMath::Sqrt(cl1->GetQ())*3;
3390 if (cl1->IsUsed(10)) sumused++;
3391 polytrack.AddPoint(x,cl1->GetY(),cl1->GetZ(),erry,errz);
3395 AliTPCclusterMI * cl2 = krp.FindNearest(y0,z0,roady,roadz);
3397 erry = (0.5)*cl2->GetSigmaY2()/TMath::Sqrt(cl2->GetQ())*3;
3398 errz = (0.5)*cl2->GetSigmaZ2()/TMath::Sqrt(cl2->GetQ())*3;
3399 if (cl2->IsUsed(10)) sumused++;
3400 polytrack.AddPoint(x,cl2->GetY(),cl2->GetZ(),erry,errz);
3403 if (sumused>0) continue;
3405 polytrack.UpdateParameters();
3411 nfoundable = polytrack.GetN();
3412 nfound = nfoundable;
3414 for (Int_t ddrow = iter+1; ddrow<drow;ddrow++){
3415 Float_t maxdist = 0.8*(1.+3./(ddrow));
3416 for (Int_t delta = -1;delta<=1;delta+=2){
3417 Int_t row = row0+ddrow*delta;
3418 kr = &(fSectors[sec][row]);
3419 Double_t xn = kr->GetX();
3420 Double_t ymax = fSectors->GetMaxY(row)-kr->GetDeadZone()-1.5;
3421 polytrack.GetFitPoint(xn,yn,zn);
3422 if (TMath::Abs(yn)>ymax) continue;
3424 AliTPCclusterMI * cln = kr->FindNearest(yn,zn,roady,roadz);
3426 Float_t dist = TMath::Sqrt( (yn-cln->GetY())*(yn-cln->GetY())+(zn-cln->GetZ())*(zn-cln->GetZ()));
3429 erry = (dist+0.3)*cln->GetSigmaY2()/TMath::Sqrt(cln->GetQ())*(1.+1./(ddrow));
3430 errz = (dist+0.3)*cln->GetSigmaZ2()/TMath::Sqrt(cln->GetQ())*(1.+1./(ddrow));
3431 if (cln->IsUsed(10)) {
3432 // printf("used\n");
3440 polytrack.AddPoint(xn,cln->GetY(),cln->GetZ(),erry, errz);
3445 if ( (sumused>3) || (sumused>0.5*nfound) || (nfound<0.6*nfoundable)) break;
3446 polytrack.UpdateParameters();
3449 if ( (sumused>3) || (sumused>0.5*nfound)) {
3450 //printf("sumused %d\n",sumused);
3455 polytrack.GetFitDerivation(kr0.GetX(),dy,dz);
3456 AliTPCpolyTrack track2;
3458 polytrack.Refit(track2,0.5+TMath::Abs(dy)*0.3,0.4+TMath::Abs(dz)*0.3);
3459 if (track2.GetN()<0.5*nfoundable) continue;
3462 if ((nfound>0.6*nfoundable) &&( nfoundable>0.4*(i1-i2))) {
3464 // test seed with and without constrain
3465 for (Int_t constrain=0; constrain<=0;constrain++){
3466 // add polytrack candidate
3468 Double_t x[5], c[15];
3469 Double_t x1,x2,x3,y1,y2,y3,z1,z2,z3;
3470 track2.GetBoundaries(x3,x1);
3472 track2.GetFitPoint(x1,y1,z1);
3473 track2.GetFitPoint(x2,y2,z2);
3474 track2.GetFitPoint(x3,y3,z3);
3476 //is track pointing to the vertex ?
3479 polytrack.GetFitPoint(x0,y0,z0);
3492 x[4]=F1(x1,y1,x2,y2,x3,y3);
3494 // if (TMath::Abs(x[4]) >= cuts[0]) continue; //
3495 x[2]=F2(x1,y1,x2,y2,x3,y3);
3497 //if (TMath::Abs(x[4]*x1-x[2]) >= cuts[1]) continue;
3498 //x[3]=F3(x1,y1,x2,y2,z1,z2);
3499 x[3]=F3n(x1,y1,x3,y3,z1,z3,x[4]);
3500 //if (TMath::Abs(x[3]) > cuts[3]) continue;
3503 Double_t sy =0.1, sz =0.1;
3504 Double_t sy1=0.02, sz1=0.02;
3505 Double_t sy2=0.02, sz2=0.02;
3509 sy3=25000*x[4]*x[4]+0.1, sy=0.1, sz=0.1;
3512 Double_t f40=(F1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
3513 Double_t f42=(F1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
3514 Double_t f43=(F1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
3515 Double_t f20=(F2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
3516 Double_t f22=(F2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
3517 Double_t f23=(F2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
3519 Double_t f30=(F3(x1,y1+sy,x3,y3,z1,z3)-x[3])/sy;
3520 Double_t f31=(F3(x1,y1,x3,y3,z1+sz,z3)-x[3])/sz;
3521 Double_t f32=(F3(x1,y1,x3,y3+sy,z1,z3)-x[3])/sy;
3522 Double_t f34=(F3(x1,y1,x3,y3,z1,z3+sz)-x[3])/sz;
3527 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3528 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3529 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3530 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3531 c[13]=f30*sy1*f40+f32*sy2*f42;
3532 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3534 //Int_t row1 = fSectors->GetRowNumber(x1);
3535 Int_t row1 = GetRowNumber(x1);
3539 seed->~AliTPCseed();
3540 AliTPCseed *track=new(seed) AliTPCseed(x1,sec*alpha+shift,x,c,index);
3541 track->SetIsSeeding(kTRUE);
3542 Int_t rc=FollowProlongation(*track, i2);
3543 if (constrain) track->SetBConstrain(1);
3545 track->SetBConstrain(0);
3546 track->SetLastPoint(row1+fInnerSec->GetNRows()); // first cluster in track position
3547 track->SetFirstPoint(track->GetLastPoint());
3549 if (rc==0 || track->GetNumberOfClusters()<(i1-i2)*0.5 ||
3550 track->GetNumberOfClusters() < track->GetNFoundable()*0.6 ||
3551 track->GetNShared()>0.4*track->GetNumberOfClusters()) {
3554 seed->~AliTPCseed();
3557 arr->AddLast(track);
3558 seed = new AliTPCseed;
3562 } // if accepted seed
3565 Info("MakeSeeds2","\nSeeding statiistic:\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin3);
3571 AliTPCseed *AliTPCtrackerMI::MakeSeed(AliTPCseed *track, Float_t r0, Float_t r1, Float_t r2)
3575 //reseed using track points
3576 Int_t p0 = int(r0*track->GetNumberOfClusters()); // point 0
3577 Int_t p1 = int(r1*track->GetNumberOfClusters());
3578 Int_t p2 = int(r2*track->GetNumberOfClusters()); // last point
3580 Double_t x0[3],x1[3],x2[3];
3581 for (Int_t i=0;i<3;i++){
3587 // find track position at given ratio of the length
3588 Int_t sec0=0, sec1=0, sec2=0;
3591 for (Int_t i=0;i<160;i++){
3592 if (track->GetClusterPointer(i)){
3594 AliTPCTrackerPoint *trpoint =track->GetTrackPoint(i);
3595 if ( (index<p0) || x0[0]<0 ){
3596 if (trpoint->GetX()>1){
3597 clindex = track->GetClusterIndex2(i);
3599 x0[0] = trpoint->GetX();
3600 x0[1] = trpoint->GetY();
3601 x0[2] = trpoint->GetZ();
3602 sec0 = ((clindex&0xff000000)>>24)%18;
3607 if ( (index<p1) &&(trpoint->GetX()>1)){
3608 clindex = track->GetClusterIndex2(i);
3610 x1[0] = trpoint->GetX();
3611 x1[1] = trpoint->GetY();
3612 x1[2] = trpoint->GetZ();
3613 sec1 = ((clindex&0xff000000)>>24)%18;
3616 if ( (index<p2) &&(trpoint->GetX()>1)){
3617 clindex = track->GetClusterIndex2(i);
3619 x2[0] = trpoint->GetX();
3620 x2[1] = trpoint->GetY();
3621 x2[2] = trpoint->GetZ();
3622 sec2 = ((clindex&0xff000000)>>24)%18;
3629 Double_t alpha, cs,sn, xx2,yy2;
3631 alpha = (sec1-sec2)*fSectors->GetAlpha();
3632 cs = TMath::Cos(alpha);
3633 sn = TMath::Sin(alpha);
3634 xx2= x1[0]*cs-x1[1]*sn;
3635 yy2= x1[0]*sn+x1[1]*cs;
3639 alpha = (sec0-sec2)*fSectors->GetAlpha();
3640 cs = TMath::Cos(alpha);
3641 sn = TMath::Sin(alpha);
3642 xx2= x0[0]*cs-x0[1]*sn;
3643 yy2= x0[0]*sn+x0[1]*cs;
3649 Double_t x[5],c[15];
3653 x[4]=F1(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]);
3654 // if (x[4]>1) return 0;
3655 x[2]=F2(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]);
3656 x[3]=F3n(x2[0],x2[1],x0[0],x0[1],x2[2],x0[2],x[4]);
3657 //if (TMath::Abs(x[3]) > 2.2) return 0;
3658 //if (TMath::Abs(x[2]) > 1.99) return 0;
3660 Double_t sy =0.1, sz =0.1;
3662 Double_t sy1=0.02+track->GetSigmaY2(), sz1=0.02+track->GetSigmaZ2();
3663 Double_t sy2=0.01+track->GetSigmaY2(), sz2=0.01+track->GetSigmaZ2();
3664 Double_t sy3=0.01+track->GetSigmaY2();
3666 Double_t f40=(F1(x2[0],x2[1]+sy,x1[0],x1[1],x0[0],x0[1])-x[4])/sy;
3667 Double_t f42=(F1(x2[0],x2[1],x1[0],x1[1]+sy,x0[0],x0[1])-x[4])/sy;
3668 Double_t f43=(F1(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]+sy)-x[4])/sy;
3669 Double_t f20=(F2(x2[0],x2[1]+sy,x1[0],x1[1],x0[0],x0[1])-x[2])/sy;
3670 Double_t f22=(F2(x2[0],x2[1],x1[0],x1[1]+sy,x0[0],x0[1])-x[2])/sy;
3671 Double_t f23=(F2(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]+sy)-x[2])/sy;
3673 Double_t f30=(F3(x2[0],x2[1]+sy,x0[0],x0[1],x2[2],x0[2])-x[3])/sy;
3674 Double_t f31=(F3(x2[0],x2[1],x0[0],x0[1],x2[2]+sz,x0[2])-x[3])/sz;
3675 Double_t f32=(F3(x2[0],x2[1],x0[0],x0[1]+sy,x2[2],x0[2])-x[3])/sy;
3676 Double_t f34=(F3(x2[0],x2[1],x0[0],x0[1],x2[2],x0[2]+sz)-x[3])/sz;
3681 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3682 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3683 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3684 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3685 c[13]=f30*sy1*f40+f32*sy2*f42;
3686 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3688 // Int_t row1 = fSectors->GetRowNumber(x2[0]);
3689 AliTPCseed *seed=new AliTPCseed(x2[0], sec2*fSectors->GetAlpha()+fSectors->GetAlphaShift(), x, c, 0);
3690 // Double_t y0,z0,y1,z1, y2,z2;
3691 //seed->GetProlongation(x0[0],y0,z0);
3692 // seed->GetProlongation(x1[0],y1,z1);
3693 //seed->GetProlongation(x2[0],y2,z2);
3695 seed->SetLastPoint(pp2);
3696 seed->SetFirstPoint(pp2);
3703 AliTPCseed *AliTPCtrackerMI::ReSeed(AliTPCseed *track, Float_t r0, Float_t r1, Float_t r2)
3707 //reseed using founded clusters
3709 // Find the number of clusters
3710 Int_t nclusters = 0;
3711 for (Int_t irow=0;irow<160;irow++){
3712 if (track->GetClusterIndex(irow)>0) nclusters++;
3716 ipos[0] = TMath::Max(int(r0*nclusters),0); // point 0 cluster
3717 ipos[1] = TMath::Min(int(r1*nclusters),nclusters-1); //
3718 ipos[2] = TMath::Min(int(r2*nclusters),nclusters-1); // last point
3722 Int_t row[3],sec[3]={0,0,0};
3724 // find track row position at given ratio of the length
3726 for (Int_t irow=0;irow<160;irow++){
3727 if (track->GetClusterIndex2(irow)<0) continue;
3729 for (Int_t ipoint=0;ipoint<3;ipoint++){
3730 if (index<=ipos[ipoint]) row[ipoint] = irow;
3734 //Get cluster and sector position
3735 for (Int_t ipoint=0;ipoint<3;ipoint++){
3736 Int_t clindex = track->GetClusterIndex2(row[ipoint]);
3737 AliTPCclusterMI * cl = GetClusterMI(clindex);
3740 // AliTPCclusterMI * cl = GetClusterMI(clindex);
3743 sec[ipoint] = ((clindex&0xff000000)>>24)%18;
3744 xyz[ipoint][0] = GetXrow(row[ipoint]);
3745 xyz[ipoint][1] = cl->GetY();
3746 xyz[ipoint][2] = cl->GetZ();
3750 // Calculate seed state vector and covariance matrix
3752 Double_t alpha, cs,sn, xx2,yy2;
3754 alpha = (sec[1]-sec[2])*fSectors->GetAlpha();
3755 cs = TMath::Cos(alpha);
3756 sn = TMath::Sin(alpha);
3757 xx2= xyz[1][0]*cs-xyz[1][1]*sn;
3758 yy2= xyz[1][0]*sn+xyz[1][1]*cs;
3762 alpha = (sec[0]-sec[2])*fSectors->GetAlpha();
3763 cs = TMath::Cos(alpha);
3764 sn = TMath::Sin(alpha);
3765 xx2= xyz[0][0]*cs-xyz[0][1]*sn;
3766 yy2= xyz[0][0]*sn+xyz[0][1]*cs;
3772 Double_t x[5],c[15];
3776 x[4]=F1(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
3777 x[2]=F2(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
3778 x[3]=F3n(xyz[2][0],xyz[2][1],xyz[0][0],xyz[0][1],xyz[2][2],xyz[0][2],x[4]);
3780 Double_t sy =0.1, sz =0.1;
3782 Double_t sy1=0.2, sz1=0.2;
3783 Double_t sy2=0.2, sz2=0.2;
3786 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;
3787 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;
3788 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;
3789 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;
3790 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;
3791 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;
3793 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;
3794 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;
3795 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;
3796 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;
3801 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3802 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3803 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3804 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3805 c[13]=f30*sy1*f40+f32*sy2*f42;
3806 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3808 // Int_t row1 = fSectors->GetRowNumber(xyz[2][0]);
3809 AliTPCseed *seed=new AliTPCseed(xyz[2][0], sec[2]*fSectors->GetAlpha()+fSectors->GetAlphaShift(), x, c, 0);
3810 seed->SetLastPoint(row[2]);
3811 seed->SetFirstPoint(row[2]);
3816 AliTPCseed *AliTPCtrackerMI::ReSeed(AliTPCseed *track,Int_t r0, Bool_t forward)
3820 //reseed using founded clusters
3823 Int_t row[3]={0,0,0};
3824 Int_t sec[3]={0,0,0};
3826 // forward direction
3828 for (Int_t irow=r0;irow<160;irow++){
3829 if (track->GetClusterIndex(irow)>0){
3834 for (Int_t irow=160;irow>r0;irow--){
3835 if (track->GetClusterIndex(irow)>0){
3840 for (Int_t irow=row[2]-15;irow>row[0];irow--){
3841 if (track->GetClusterIndex(irow)>0){
3849 for (Int_t irow=0;irow<r0;irow++){
3850 if (track->GetClusterIndex(irow)>0){
3855 for (Int_t irow=r0;irow>0;irow--){
3856 if (track->GetClusterIndex(irow)>0){
3861 for (Int_t irow=row[2]-15;irow>row[0];irow--){
3862 if (track->GetClusterIndex(irow)>0){
3869 if ((row[2]-row[0])<20) return 0;
3870 if (row[1]==0) return 0;
3873 //Get cluster and sector position
3874 for (Int_t ipoint=0;ipoint<3;ipoint++){
3875 Int_t clindex = track->GetClusterIndex2(row[ipoint]);
3876 AliTPCclusterMI * cl = GetClusterMI(clindex);
3879 // AliTPCclusterMI * cl = GetClusterMI(clindex);
3882 sec[ipoint] = ((clindex&0xff000000)>>24)%18;
3883 xyz[ipoint][0] = GetXrow(row[ipoint]);
3884 AliTPCTrackerPoint * point = track->GetTrackPoint(row[ipoint]);
3885 if (point&&ipoint<2){
3887 xyz[ipoint][1] = point->GetY();
3888 xyz[ipoint][2] = point->GetZ();
3891 xyz[ipoint][1] = cl->GetY();
3892 xyz[ipoint][2] = cl->GetZ();
3899 // Calculate seed state vector and covariance matrix
3901 Double_t alpha, cs,sn, xx2,yy2;
3903 alpha = (sec[1]-sec[2])*fSectors->GetAlpha();
3904 cs = TMath::Cos(alpha);
3905 sn = TMath::Sin(alpha);
3906 xx2= xyz[1][0]*cs-xyz[1][1]*sn;
3907 yy2= xyz[1][0]*sn+xyz[1][1]*cs;
3911 alpha = (sec[0]-sec[2])*fSectors->GetAlpha();
3912 cs = TMath::Cos(alpha);
3913 sn = TMath::Sin(alpha);
3914 xx2= xyz[0][0]*cs-xyz[0][1]*sn;
3915 yy2= xyz[0][0]*sn+xyz[0][1]*cs;
3921 Double_t x[5],c[15];
3925 x[4]=F1(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
3926 x[2]=F2(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
3927 x[3]=F3n(xyz[2][0],xyz[2][1],xyz[0][0],xyz[0][1],xyz[2][2],xyz[0][2],x[4]);
3929 Double_t sy =0.1, sz =0.1;
3931 Double_t sy1=0.2, sz1=0.2;
3932 Double_t sy2=0.2, sz2=0.2;
3935 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;
3936 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;
3937 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;
3938 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;
3939 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;
3940 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;
3942 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;
3943 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;
3944 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;
3945 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;
3950 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3951 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3952 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3953 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3954 c[13]=f30*sy1*f40+f32*sy2*f42;
3955 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3957 // Int_t row1 = fSectors->GetRowNumber(xyz[2][0]);
3958 AliTPCseed *seed=new AliTPCseed(xyz[2][0], sec[2]*fSectors->GetAlpha()+fSectors->GetAlphaShift(), x, c, 0);
3959 seed->SetLastPoint(row[2]);
3960 seed->SetFirstPoint(row[2]);
3961 for (Int_t i=row[0];i<row[2];i++){
3962 seed->SetClusterIndex(i, track->GetClusterIndex(i));
3970 void AliTPCtrackerMI::FindMultiMC(TObjArray * array, AliESDEvent */*esd*/, Int_t iter)
3973 // find multi tracks - THIS FUNCTION IS ONLY FOR DEBUG PURPOSES
3975 // Use AliTPCReconstructor::StreamLevel()>2 if you want to tune parameters - cuts
3977 // Two reasons to have multiple find tracks
3978 // 1. Curling tracks can be find more than once
3979 // 2. Splitted tracks
3980 // a.) Multiple seeding to increase tracking efficiency - (~ 100% reached)
3981 // b.) Edge effect on the sector boundaries
3984 // Algorithm done in 2 phases - because of CPU consumption
3985 // it is n^2 algorithm - for lead-lead 20000x20000 combination are investigated
3987 // Algorihm for curling tracks sign:
3988 // 1 phase -makes a very rough fast cuts to minimize combinatorics
3989 // a.) opposite sign
3990 // b.) one of the tracks - not pointing to the primary vertex -
3991 // c.) delta tan(theta)
3993 // 2 phase - calculates DCA between tracks - time consument
3998 // General cuts - for splitted tracks and for curling tracks
4000 const Float_t kMaxdPhi = 0.2; // maximal distance in phi
4002 // Curling tracks cuts
4007 Int_t nentries = array->GetEntriesFast();
4008 AliHelix *helixes = new AliHelix[nentries];
4009 Float_t *xm = new Float_t[nentries];
4010 Float_t *dz0 = new Float_t[nentries];
4011 Float_t *dz1 = new Float_t[nentries];
4017 // Find track COG in x direction - point with best defined parameters
4019 for (Int_t i=0;i<nentries;i++){
4020 AliTPCseed* track = (AliTPCseed*)array->At(i);
4021 if (!track) continue;
4022 track->SetCircular(0);
4023 new (&helixes[i]) AliHelix(*track);
4027 track->GetDZ(GetX(),GetY(),GetZ(),GetBz(),dz);
4030 for (Int_t icl=0; icl<160; icl++){
4031 AliTPCclusterMI * cl = track->GetClusterPointer(icl);
4037 if (ncl>0) xm[i]/=Float_t(ncl);
4039 TTreeSRedirector &cstream = *fDebugStreamer;
4041 for (Int_t i0=0;i0<nentries;i0++){
4042 AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
4043 if (!track0) continue;
4044 Float_t xc0 = helixes[i0].GetHelix(6);
4045 Float_t yc0 = helixes[i0].GetHelix(7);
4046 Float_t r0 = helixes[i0].GetHelix(8);
4047 Float_t rc0 = TMath::Sqrt(xc0*xc0+yc0*yc0);
4048 Float_t fi0 = TMath::ATan2(yc0,xc0);
4050 for (Int_t i1=i0+1;i1<nentries;i1++){
4051 AliTPCseed * track1 = (AliTPCseed*)array->At(i1);
4052 if (!track1) continue;
4053 Int_t lab0=track0->GetLabel();
4054 Int_t lab1=track1->GetLabel();
4055 if (TMath::Abs(lab0)!=TMath::Abs(lab1)) continue;
4057 Float_t xc1 = helixes[i1].GetHelix(6);
4058 Float_t yc1 = helixes[i1].GetHelix(7);
4059 Float_t r1 = helixes[i1].GetHelix(8);
4060 Float_t rc1 = TMath::Sqrt(xc1*xc1+yc1*yc1);
4061 Float_t fi1 = TMath::ATan2(yc1,xc1);
4063 Float_t dfi = fi0-fi1;
4066 if (dfi>1.5*TMath::Pi()) dfi-=TMath::Pi(); // take care about edge effect
4067 if (dfi<-1.5*TMath::Pi()) dfi+=TMath::Pi(); //
4068 if (TMath::Abs(dfi)>kMaxdPhi&&helixes[i0].GetHelix(4)*helixes[i1].GetHelix(4)<0){
4070 // if short tracks with undefined sign
4071 fi1 = -TMath::ATan2(yc1,-xc1);
4074 Float_t dtheta = TMath::Abs(track0->GetTgl()-track1->GetTgl())<TMath::Abs(track0->GetTgl()+track1->GetTgl())? track0->GetTgl()-track1->GetTgl():track0->GetTgl()+track1->GetTgl();
4077 // debug stream to tune "fast cuts"
4079 Double_t dist[3]; // distance at X
4080 Double_t mdist[3]={0,0,0}; // mean distance X+-40cm
4081 track0->GetDistance(track1,0.5*(xm[i0]+xm[i1])-40.,dist,AliTracker::GetBz());
4082 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[i]);
4083 track0->GetDistance(track1,0.5*(xm[i0]+xm[i1])+40.,dist,AliTracker::GetBz());
4084 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[i]);
4085 track0->GetDistance(track1,0.5*(xm[i0]+xm[i1]),dist,AliTracker::GetBz());
4086 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[i]);
4087 for (Int_t i=0;i<3;i++) mdist[i]*=0.33333;
4091 for (Int_t icl=0; icl<160; icl++){
4092 AliTPCclusterMI * cl0 = track0->GetClusterPointer(icl);
4093 AliTPCclusterMI * cl1 = track1->GetClusterPointer(icl);
4096 if (cl0==cl1) sums++;
4104 "Tr0.="<<track0<< // seed0
4105 "Tr1.="<<track1<< // seed1
4106 "h0.="<<&helixes[i0]<<
4107 "h1.="<<&helixes[i1]<<
4109 "sum="<<sum<< //the sum of rows with cl in both
4110 "sums="<<sums<< //the sum of shared clusters
4111 "xm0="<<xm[i0]<< // the center of track
4112 "xm1="<<xm[i1]<< // the x center of track
4113 // General cut variables
4114 "dfi="<<dfi<< // distance in fi angle
4115 "dtheta="<<dtheta<< // distance int theta angle
4121 "dist0="<<dist[0]<< //distance x
4122 "dist1="<<dist[1]<< //distance y
4123 "dist2="<<dist[2]<< //distance z
4124 "mdist0="<<mdist[0]<< //distance x
4125 "mdist1="<<mdist[1]<< //distance y
4126 "mdist2="<<mdist[2]<< //distance z
4139 if (AliTPCReconstructor::StreamLevel()>1) {
4140 AliInfo("Time for curling tracks removal DEBUGGING MC");
4146 void AliTPCtrackerMI::FindSplitted(TObjArray * array, AliESDEvent */*esd*/, Int_t iter)
4150 // Two reasons to have multiple find tracks
4151 // 1. Curling tracks can be find more than once
4152 // 2. Splitted tracks
4153 // a.) Multiple seeding to increase tracking efficiency - (~ 100% reached)
4154 // b.) Edge effect on the sector boundaries
4156 // This function tries to find tracks closed in the parametric space
4158 // cut logic if distance is bigger than cut continue - Do Nothing
4159 const Float_t kMaxdTheta = 0.05; // maximal distance in theta
4160 const Float_t kMaxdPhi = 0.05; // maximal deistance in phi
4161 const Float_t kdelta = 40.; //delta r to calculate track distance
4163 // const Float_t kMaxDist0 = 1.; // maximal distance 0
4164 //const Float_t kMaxDist1 = 0.3; // maximal distance 1 - cut if track in separate rows
4167 TCut csec("csec","abs(Tr0.fRelativeSector-Tr1.fRelativeSector)<2");
4168 TCut cdtheta("cdtheta","abs(dtheta)<0.05");
4173 Int_t nentries = array->GetEntriesFast();
4174 AliHelix *helixes = new AliHelix[nentries];
4175 Float_t *xm = new Float_t[nentries];
4181 //Sort tracks according quality
4183 Int_t nseed = array->GetEntriesFast();
4184 Float_t * quality = new Float_t[nseed];
4185 Int_t * indexes = new Int_t[nseed];
4186 for (Int_t i=0; i<nseed; i++) {
4187 AliTPCseed *pt=(AliTPCseed*)array->UncheckedAt(i);
4192 pt->UpdatePoints(); //select first last max dens points
4193 Float_t * points = pt->GetPoints();
4194 if (points[3]<0.8) quality[i] =-1;
4195 quality[i] = (points[2]-points[0])+pt->GetNumberOfClusters();
4196 //prefer high momenta tracks if overlaps
4197 quality[i] *= TMath::Sqrt(TMath::Abs(pt->Pt())+0.5);
4199 TMath::Sort(nseed,quality,indexes);
4203 // Find track COG in x direction - point with best defined parameters
4205 for (Int_t i=0;i<nentries;i++){
4206 AliTPCseed* track = (AliTPCseed*)array->At(i);
4207 if (!track) continue;
4208 track->SetCircular(0);
4209 new (&helixes[i]) AliHelix(*track);
4212 for (Int_t icl=0; icl<160; icl++){
4213 AliTPCclusterMI * cl = track->GetClusterPointer(icl);
4219 if (ncl>0) xm[i]/=Float_t(ncl);
4221 TTreeSRedirector &cstream = *fDebugStreamer;
4223 for (Int_t is0=0;is0<nentries;is0++){
4224 Int_t i0 = indexes[is0];
4225 AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
4226 if (!track0) continue;
4227 if (track0->GetKinkIndexes()[0]!=0) continue;
4228 Float_t xc0 = helixes[i0].GetHelix(6);
4229 Float_t yc0 = helixes[i0].GetHelix(7);
4230 Float_t fi0 = TMath::ATan2(yc0,xc0);
4232 for (Int_t is1=is0+1;is1<nentries;is1++){
4233 Int_t i1 = indexes[is1];
4234 AliTPCseed * track1 = (AliTPCseed*)array->At(i1);
4235 if (!track1) continue;
4237 if (TMath::Abs(track0->GetRelativeSector()-track1->GetRelativeSector())>1) continue;
4238 if (track1->GetKinkIndexes()[0]>0 &&track0->GetKinkIndexes()[0]<0) continue;
4239 if (track1->GetKinkIndexes()[0]!=0) continue;
4241 Float_t dtheta = TMath::Abs(track0->GetTgl()-track1->GetTgl())<TMath::Abs(track0->GetTgl()+track1->GetTgl())? track0->GetTgl()-track1->GetTgl():track0->GetTgl()+track1->GetTgl();
4242 if (TMath::Abs(dtheta)>kMaxdTheta) continue;
4244 Float_t xc1 = helixes[i1].GetHelix(6);
4245 Float_t yc1 = helixes[i1].GetHelix(7);
4246 Float_t fi1 = TMath::ATan2(yc1,xc1);
4248 Float_t dfi = fi0-fi1;
4249 if (dfi>1.5*TMath::Pi()) dfi-=TMath::Pi(); // take care about edge effect
4250 if (dfi<-1.5*TMath::Pi()) dfi+=TMath::Pi(); //
4251 if (TMath::Abs(dfi)>kMaxdPhi&&helixes[i0].GetHelix(4)*helixes[i1].GetHelix(4)<0){
4253 // if short tracks with undefined sign
4254 fi1 = -TMath::ATan2(yc1,-xc1);
4257 if (TMath::Abs(dfi)>kMaxdPhi) continue;
4264 for (Int_t icl=0; icl<160; icl++){
4265 Int_t index0=track0->GetClusterIndex2(icl);
4266 Int_t index1=track1->GetClusterIndex2(icl);
4267 Bool_t used0 = (index0>0 && !(index0&0x8000));
4268 Bool_t used1 = (index1>0 && !(index1&0x8000));
4270 if (used0) sum0++; // used cluster0
4271 if (used1) sum1++; // used clusters1
4272 if (used0&&used1) sum++;
4273 if (index0==index1 && used0 && used1) sums++;
4277 if (sums<10) continue;
4278 if (sum<40) continue;
4279 if (sums/Float_t(TMath::Min(sum0,sum1))<0.5) continue;
4281 Double_t dist[5][4]; // distance at X
4282 Double_t mdist[4]={0,0,0,0}; // mean distance on range +-delta
4286 track0->GetDistance(track1,xm[i0],dist[0],AliTracker::GetBz());
4287 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[0][i]);
4288 track0->GetDistance(track1,xm[i1],dist[1],AliTracker::GetBz());
4289 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[1][i]);
4291 track0->GetDistance(track1,TMath::Min(xm[i1],xm[i0])-kdelta,dist[2],AliTracker::GetBz());
4292 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[2][i]);
4293 track0->GetDistance(track1,TMath::Max(xm[i1],xm[i0])+kdelta,dist[3],AliTracker::GetBz());
4294 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[3][i]);
4296 track0->GetDistance(track1,(xm[i1]+xm[i0])*0.5,dist[4],AliTracker::GetBz());
4297 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[4][i]);
4298 for (Int_t i=0;i<3;i++) mdist[i]*=0.2;
4301 Int_t lab0=track0->GetLabel();
4302 Int_t lab1=track1->GetLabel();
4303 cstream<<"Splitted"<<
4307 "Tr0.="<<track0<< // seed0
4308 "Tr1.="<<track1<< // seed1
4309 "h0.="<<&helixes[i0]<<
4310 "h1.="<<&helixes[i1]<<
4312 "sum="<<sum<< //the sum of rows with cl in both
4313 "sum0="<<sum0<< //the sum of rows with cl in 0 track
4314 "sum1="<<sum1<< //the sum of rows with cl in 1 track
4315 "sums="<<sums<< //the sum of shared clusters
4316 "xm0="<<xm[i0]<< // the center of track
4317 "xm1="<<xm[i1]<< // the x center of track
4318 // General cut variables
4319 "dfi="<<dfi<< // distance in fi angle
4320 "dtheta="<<dtheta<< // distance int theta angle
4323 "dist0="<<dist[4][0]<< //distance x
4324 "dist1="<<dist[4][1]<< //distance y
4325 "dist2="<<dist[4][2]<< //distance z
4326 "mdist0="<<mdist[0]<< //distance x
4327 "mdist1="<<mdist[1]<< //distance y
4328 "mdist2="<<mdist[2]<< //distance z
4331 delete array->RemoveAt(i1);
4336 AliInfo("Time for splitted tracks removal");
4342 void AliTPCtrackerMI::FindCurling(TObjArray * array, AliESDEvent */*esd*/, Int_t iter)
4345 // find Curling tracks
4346 // Use AliTPCReconstructor::StreamLevel()>1 if you want to tune parameters - cuts
4349 // Algorithm done in 2 phases - because of CPU consumption
4350 // it is n^2 algorithm - for lead-lead 20000x20000 combination are investigated
4351 // see detal in MC part what can be used to cut
4355 const Float_t kMaxC = 400; // maximal curvature to of the track
4356 const Float_t kMaxdTheta = 0.15; // maximal distance in theta
4357 const Float_t kMaxdPhi = 0.15; // maximal distance in phi
4358 const Float_t kPtRatio = 0.3; // ratio between pt
4359 const Float_t kMinDCAR = 2.; // distance to the primary vertex in r - see cpipe cut
4362 // Curling tracks cuts
4365 const Float_t kMaxDeltaRMax = 40; // distance in outer radius
4366 const Float_t kMaxDeltaRMin = 5.; // distance in lower radius - see cpipe cut
4367 const Float_t kMinAngle = 2.9; // angle between tracks
4368 const Float_t kMaxDist = 5; // biggest distance
4370 // The cuts can be tuned using the "MC information stored in Multi tree ==> see FindMultiMC
4373 TCut csign("csign","Tr0.fP[4]*Tr1.fP[4]<0"); //opposite sign
4374 TCut cmax("cmax","abs(Tr0.GetC())>1/400");
4375 TCut cda("cda","sqrt(dtheta^2+dfi^2)<0.15");
4376 TCut ccratio("ccratio","abs((Tr0.fP[4]+Tr1.fP[4])/(abs(Tr0.fP[4])+abs(Tr1.fP[4])))<0.3");
4377 TCut cpipe("cpipe", "min(abs(r0-rc0),abs(r1-rc1))>5");
4379 TCut cdrmax("cdrmax","abs(abs(rc0+r0)-abs(rc1+r1))<40")
4380 TCut cdrmin("cdrmin","abs(abs(rc0+r0)-abs(rc1+r1))<10")
4382 Multi->Draw("dfi","iter==0"+csign+cmax+cda+ccratio); ~94% of curling tracks fulfill
4383 Multi->Draw("min(abs(r0-rc0),abs(r1-rc1))","iter==0&&abs(lab1)==abs(lab0)"+csign+cmax+cda+ccratio+cpipe+cdrmin+cdrmax); //80%
4385 Curling2->Draw("dfi","iter==0&&abs(lab0)==abs(lab1)"+csign+cmax+cdtheta+cdfi+ccratio)
4391 Int_t nentries = array->GetEntriesFast();
4392 AliHelix *helixes = new AliHelix[nentries];
4393 for (Int_t i=0;i<nentries;i++){
4394 AliTPCseed* track = (AliTPCseed*)array->At(i);
4395 if (!track) continue;
4396 track->SetCircular(0);
4397 new (&helixes[i]) AliHelix(*track);
4403 Double_t phase[2][2],radius[2];
4407 TTreeSRedirector &cstream = *fDebugStreamer;
4409 for (Int_t i0=0;i0<nentries;i0++){
4410 AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
4411 if (!track0) continue;
4412 if (TMath::Abs(track0->GetC())<1/kMaxC) continue;
4413 Float_t xc0 = helixes[i0].GetHelix(6);
4414 Float_t yc0 = helixes[i0].GetHelix(7);
4415 Float_t r0 = helixes[i0].GetHelix(8);
4416 Float_t rc0 = TMath::Sqrt(xc0*xc0+yc0*yc0);
4417 Float_t fi0 = TMath::ATan2(yc0,xc0);
4419 for (Int_t i1=i0+1;i1<nentries;i1++){
4420 AliTPCseed * track1 = (AliTPCseed*)array->At(i1);
4421 if (!track1) continue;
4422 if (TMath::Abs(track1->GetC())<1/kMaxC) continue;
4423 Float_t xc1 = helixes[i1].GetHelix(6);
4424 Float_t yc1 = helixes[i1].GetHelix(7);
4425 Float_t r1 = helixes[i1].GetHelix(8);
4426 Float_t rc1 = TMath::Sqrt(xc1*xc1+yc1*yc1);
4427 Float_t fi1 = TMath::ATan2(yc1,xc1);
4429 Float_t dfi = fi0-fi1;
4432 if (dfi>1.5*TMath::Pi()) dfi-=TMath::Pi(); // take care about edge effect
4433 if (dfi<-1.5*TMath::Pi()) dfi+=TMath::Pi(); //
4434 Float_t dtheta = TMath::Abs(track0->GetTgl()-track1->GetTgl())<TMath::Abs(track0->GetTgl()+track1->GetTgl())? track0->GetTgl()-track1->GetTgl():track0->GetTgl()+track1->GetTgl();
4438 if (track0->GetBConstrain()&&track1->GetBConstrain()) continue; // not constrained
4439 if (track1->GetSigned1Pt()*track0->GetSigned1Pt()>0) continue; // not the same sign
4440 if ( TMath::Abs(track1->GetTgl()+track0->GetTgl())>kMaxdTheta) continue; //distance in the Theta
4441 if ( TMath::Abs(dfi)>kMaxdPhi) continue; //distance in phi
4442 if ( TMath::Sqrt(dfi*dfi+dtheta*dtheta)>kMaxdPhi) continue; //common angular offset
4444 Float_t pt0 = track0->GetSignedPt();
4445 Float_t pt1 = track1->GetSignedPt();
4446 if ((TMath::Abs(pt0+pt1)/(TMath::Abs(pt0)+TMath::Abs(pt1)))>kPtRatio) continue;
4447 if ((iter==1) && TMath::Abs(TMath::Abs(rc0+r0)-TMath::Abs(rc1+r1))>kMaxDeltaRMax) continue;
4448 if ((iter!=1) &&TMath::Abs(TMath::Abs(rc0-r0)-TMath::Abs(rc1-r1))>kMaxDeltaRMin) continue;
4449 if (TMath::Min(TMath::Abs(rc0-r0),TMath::Abs(rc1-r1))<kMinDCAR) continue;
4452 // Now find closest approach
4456 Int_t npoints = helixes[i0].GetRPHIintersections(helixes[i1], phase, radius,10);
4457 if (npoints==0) continue;
4458 helixes[i0].GetClosestPhases(helixes[i1], phase);
4462 Double_t hangles[3];
4463 helixes[i0].Evaluate(phase[0][0],xyz0);
4464 helixes[i1].Evaluate(phase[0][1],xyz1);
4466 helixes[i0].GetAngle(phase[0][0],helixes[i1],phase[0][1],hangles);
4467 Double_t deltah[2],deltabest;
4468 if (TMath::Abs(hangles[2])<kMinAngle) continue;
4472 helixes[i0].ParabolicDCA(helixes[i1],phase[0][0],phase[0][1],radius[0],deltah[0],2);
4474 helixes[i0].ParabolicDCA(helixes[i1],phase[1][0],phase[1][1],radius[1],deltah[1],2);
4475 if (deltah[1]<deltah[0]) ibest=1;
4477 deltabest = TMath::Sqrt(deltah[ibest]);
4478 helixes[i0].Evaluate(phase[ibest][0],xyz0);
4479 helixes[i1].Evaluate(phase[ibest][1],xyz1);
4480 helixes[i0].GetAngle(phase[ibest][0],helixes[i1],phase[ibest][1],hangles);
4481 Double_t radiusbest = TMath::Sqrt(radius[ibest]);
4483 if (deltabest>kMaxDist) continue;
4484 // if (mindcar+mindcaz<40 && (TMath::Abs(hangles[2])<kMinAngle ||deltabest>3)) continue;
4485 Bool_t sign =kFALSE;
4486 if (hangles[2]>kMinAngle) sign =kTRUE;
4489 // circular[i0] = kTRUE;
4490 // circular[i1] = kTRUE;
4491 if (track0->OneOverPt()<track1->OneOverPt()){
4492 track0->SetCircular(track0->GetCircular()+1);
4493 track1->SetCircular(track1->GetCircular()+2);
4496 track1->SetCircular(track1->GetCircular()+1);
4497 track0->SetCircular(track0->GetCircular()+2);
4500 if (AliTPCReconstructor::StreamLevel()>1){
4502 //debug stream to tune "fine" cuts
4503 Int_t lab0=track0->GetLabel();
4504 Int_t lab1=track1->GetLabel();
4505 cstream<<"Curling2"<<
4521 "npoints="<<npoints<<
4522 "hangles0="<<hangles[0]<<
4523 "hangles1="<<hangles[1]<<
4524 "hangles2="<<hangles[2]<<
4527 "radius="<<radiusbest<<
4528 "deltabest="<<deltabest<<
4529 "phase0="<<phase[ibest][0]<<
4530 "phase1="<<phase[ibest][1]<<
4538 if (AliTPCReconstructor::StreamLevel()>1) {
4539 AliInfo("Time for curling tracks removal");
4548 void AliTPCtrackerMI::FindKinks(TObjArray * array, AliESDEvent *esd)
4555 TObjArray *kinks= new TObjArray(10000);
4556 // TObjArray *v0s= new TObjArray(10000);
4557 Int_t nentries = array->GetEntriesFast();
4558 AliHelix *helixes = new AliHelix[nentries];
4559 Int_t *sign = new Int_t[nentries];
4560 Int_t *nclusters = new Int_t[nentries];
4561 Float_t *alpha = new Float_t[nentries];
4562 AliKink *kink = new AliKink();
4563 Int_t * usage = new Int_t[nentries];
4564 Float_t *zm = new Float_t[nentries];
4565 Float_t *z0 = new Float_t[nentries];
4566 Float_t *fim = new Float_t[nentries];
4567 Float_t *shared = new Float_t[nentries];
4568 Bool_t *circular = new Bool_t[nentries];
4569 Float_t *dca = new Float_t[nentries];
4570 //const AliESDVertex * primvertex = esd->GetVertex();
4572 // nentries = array->GetEntriesFast();
4577 for (Int_t i=0;i<nentries;i++){
4580 AliTPCseed* track = (AliTPCseed*)array->At(i);
4581 if (!track) continue;
4582 track->SetCircular(0);
4584 track->UpdatePoints();
4585 if (( track->GetPoints()[2]- track->GetPoints()[0])>5 && track->GetPoints()[3]>0.8){
4587 nclusters[i]=track->GetNumberOfClusters();
4588 alpha[i] = track->GetAlpha();
4589 new (&helixes[i]) AliHelix(*track);
4591 helixes[i].Evaluate(0,xyz);
4592 sign[i] = (track->GetC()>0) ? -1:1;
4595 if (track->GetProlongation(x,y,z)){
4597 fim[i] = alpha[i]+TMath::ATan2(y,x);
4600 zm[i] = track->GetZ();
4604 circular[i]= kFALSE;
4605 if (track->GetProlongation(0,y,z)) z0[i] = z;
4606 dca[i] = track->GetD(0,0);
4612 Int_t ncandidates =0;
4615 Double_t phase[2][2],radius[2];
4618 // Find circling track
4619 TTreeSRedirector &cstream = *fDebugStreamer;
4621 for (Int_t i0=0;i0<nentries;i0++){
4622 AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
4623 if (!track0) continue;
4624 if (track0->GetNumberOfClusters()<40) continue;
4625 if (TMath::Abs(1./track0->GetC())>200) continue;
4626 for (Int_t i1=i0+1;i1<nentries;i1++){
4627 AliTPCseed * track1 = (AliTPCseed*)array->At(i1);
4628 if (!track1) continue;
4629 if (track1->GetNumberOfClusters()<40) continue;
4630 if ( TMath::Abs(track1->GetTgl()+track0->GetTgl())>0.1) continue;
4631 if (track0->GetBConstrain()&&track1->GetBConstrain()) continue;
4632 if (TMath::Abs(1./track1->GetC())>200) continue;
4633 if (track1->GetSigned1Pt()*track0->GetSigned1Pt()>0) continue;
4634 if (track1->GetTgl()*track0->GetTgl()>0) continue;
4635 if (TMath::Max(TMath::Abs(1./track0->GetC()),TMath::Abs(1./track1->GetC()))>190) continue;
4636 if (track0->GetBConstrain()&&track1->OneOverPt()<track0->OneOverPt()) continue; //returning - lower momenta
4637 if (track1->GetBConstrain()&&track0->OneOverPt()<track1->OneOverPt()) continue; //returning - lower momenta
4639 Float_t mindcar = TMath::Min(TMath::Abs(dca[i0]),TMath::Abs(dca[i1]));
4640 if (mindcar<5) continue;
4641 Float_t mindcaz = TMath::Min(TMath::Abs(z0[i0]-GetZ()),TMath::Abs(z0[i1]-GetZ()));
4642 if (mindcaz<5) continue;
4643 if (mindcar+mindcaz<20) continue;
4646 Float_t xc0 = helixes[i0].GetHelix(6);
4647 Float_t yc0 = helixes[i0].GetHelix(7);
4648 Float_t r0 = helixes[i0].GetHelix(8);
4649 Float_t xc1 = helixes[i1].GetHelix(6);
4650 Float_t yc1 = helixes[i1].GetHelix(7);
4651 Float_t r1 = helixes[i1].GetHelix(8);
4653 Float_t rmean = (r0+r1)*0.5;
4654 Float_t delta =TMath::Sqrt((xc1-xc0)*(xc1-xc0)+(yc1-yc0)*(yc1-yc0));
4655 //if (delta>30) continue;
4656 if (delta>rmean*0.25) continue;
4657 if (TMath::Abs(r0-r1)/rmean>0.3) continue;
4659 Int_t npoints = helixes[i0].GetRPHIintersections(helixes[i1], phase, radius,10);
4660 if (npoints==0) continue;
4661 helixes[i0].GetClosestPhases(helixes[i1], phase);
4665 Double_t hangles[3];
4666 helixes[i0].Evaluate(phase[0][0],xyz0);
4667 helixes[i1].Evaluate(phase[0][1],xyz1);
4669 helixes[i0].GetAngle(phase[0][0],helixes[i1],phase[0][1],hangles);
4670 Double_t deltah[2],deltabest;
4671 if (hangles[2]<2.8) continue;
4674 helixes[i0].ParabolicDCA(helixes[i1],phase[0][0],phase[0][1],radius[0],deltah[0],2);
4676 helixes[i0].ParabolicDCA(helixes[i1],phase[1][0],phase[1][1],radius[1],deltah[1],2);
4677 if (deltah[1]<deltah[0]) ibest=1;
4679 deltabest = TMath::Sqrt(deltah[ibest]);
4680 helixes[i0].Evaluate(phase[ibest][0],xyz0);
4681 helixes[i1].Evaluate(phase[ibest][1],xyz1);
4682 helixes[i0].GetAngle(phase[ibest][0],helixes[i1],phase[ibest][1],hangles);
4683 Double_t radiusbest = TMath::Sqrt(radius[ibest]);
4685 if (deltabest>6) continue;
4686 if (mindcar+mindcaz<40 && (hangles[2]<3.12||deltabest>3)) continue;
4687 Bool_t sign =kFALSE;
4688 if (hangles[2]>3.06) sign =kTRUE;
4691 circular[i0] = kTRUE;
4692 circular[i1] = kTRUE;
4693 if (track0->OneOverPt()<track1->OneOverPt()){
4694 track0->SetCircular(track0->GetCircular()+1);
4695 track1->SetCircular(track1->GetCircular()+2);
4698 track1->SetCircular(track1->GetCircular()+1);
4699 track0->SetCircular(track0->GetCircular()+2);
4702 if (sign&&AliTPCReconstructor::StreamLevel()>1){
4704 Int_t lab0=track0->GetLabel();
4705 Int_t lab1=track1->GetLabel();
4706 cstream<<"Curling"<<
4713 "mindcar="<<mindcar<<
4714 "mindcaz="<<mindcaz<<
4717 "npoints="<<npoints<<
4718 "hangles0="<<hangles[0]<<
4719 "hangles2="<<hangles[2]<<
4724 "radius="<<radiusbest<<
4725 "deltabest="<<deltabest<<
4726 "phase0="<<phase[ibest][0]<<
4727 "phase1="<<phase[ibest][1]<<
4737 for (Int_t i =0;i<nentries;i++){
4738 if (sign[i]==0) continue;
4739 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
4746 Double_t cradius0 = 40*40;
4747 Double_t cradius1 = 270*270;
4750 Double_t cdist3=0.55;
4751 for (Int_t j =i+1;j<nentries;j++){
4753 if (sign[j]*sign[i]<1) continue;
4754 if ( (nclusters[i]+nclusters[j])>200) continue;
4755 if ( (nclusters[i]+nclusters[j])<80) continue;
4756 if ( TMath::Abs(zm[i]-zm[j])>60.) continue;
4757 if ( TMath::Abs(fim[i]-fim[j])>0.6 && TMath::Abs(fim[i]-fim[j])<5.7 ) continue;
4758 //AliTPCseed * track1 = (AliTPCseed*)array->At(j); Double_t phase[2][2],radius[2];
4759 Int_t npoints = helixes[i].GetRPHIintersections(helixes[j], phase, radius,20);
4760 if (npoints<1) continue;
4763 if (radius[0]<cradius0||radius[0]>cradius1) continue;
4766 if ( (radius[0]<cradius0||radius[0]>cradius1) && (radius[1]<cradius0||radius[1]>cradius1) ) continue;
4769 Double_t delta1=10000,delta2=10000;
4770 // cuts on the intersection radius
4771 helixes[i].LinearDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta1);
4772 if (radius[0]<20&&delta1<1) continue; //intersection at vertex
4773 if (radius[0]<10&&delta1<3) continue; //intersection at vertex
4775 helixes[i].LinearDCA(helixes[j],phase[1][0],phase[1][1],radius[1],delta2);
4776 if (radius[1]<20&&delta2<1) continue; //intersection at vertex
4777 if (radius[1]<10&&delta2<3) continue; //intersection at vertex
4780 Double_t distance1 = TMath::Min(delta1,delta2);
4781 if (distance1>cdist1) continue; // cut on DCA linear approximation
4783 npoints = helixes[i].GetRPHIintersections(helixes[j], phase, radius,20);
4784 helixes[i].ParabolicDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta1);
4785 if (radius[0]<20&&delta1<1) continue; //intersection at vertex
4786 if (radius[0]<10&&delta1<3) continue; //intersection at vertex
4789 helixes[i].ParabolicDCA(helixes[j],phase[1][0],phase[1][1],radius[1],delta2);
4790 if (radius[1]<20&&delta2<1) continue; //intersection at vertex
4791 if (radius[1]<10&&delta2<3) continue; //intersection at vertex
4793 distance1 = TMath::Min(delta1,delta2);
4796 rkink = TMath::Sqrt(radius[0]);
4799 rkink = TMath::Sqrt(radius[1]);
4801 if (distance1>cdist2) continue;
4804 AliTPCseed * track1 = (AliTPCseed*)array->At(j);
4807 Int_t row0 = GetRowNumber(rkink);
4808 if (row0<10) continue;
4809 if (row0>150) continue;
4812 Float_t dens00=-1,dens01=-1;
4813 Float_t dens10=-1,dens11=-1;
4815 Int_t found,foundable,shared;
4816 track0->GetClusterStatistic(0,row0-5, found, foundable,shared,kFALSE);
4817 if (foundable>5) dens00 = Float_t(found)/Float_t(foundable);
4818 track0->GetClusterStatistic(row0+5,155, found, foundable,shared,kFALSE);
4819 if (foundable>5) dens01 = Float_t(found)/Float_t(foundable);
4821 track1->GetClusterStatistic(0,row0-5, found, foundable,shared,kFALSE);
4822 if (foundable>10) dens10 = Float_t(found)/Float_t(foundable);
4823 track1->GetClusterStatistic(row0+5,155, found, foundable,shared,kFALSE);
4824 if (foundable>10) dens11 = Float_t(found)/Float_t(foundable);
4826 if (dens00<dens10 && dens01<dens11) continue;
4827 if (dens00>dens10 && dens01>dens11) continue;
4828 if (TMath::Max(dens00,dens10)<0.1) continue;
4829 if (TMath::Max(dens01,dens11)<0.3) continue;
4831 if (TMath::Min(dens00,dens10)>0.6) continue;
4832 if (TMath::Min(dens01,dens11)>0.6) continue;
4835 AliTPCseed * ktrack0, *ktrack1;
4844 if (TMath::Abs(ktrack0->GetC())>5) continue; // cut on the curvature for mother particle
4845 AliExternalTrackParam paramm(*ktrack0);
4846 AliExternalTrackParam paramd(*ktrack1);
4847 if (row0>60&&ktrack1->GetReference().GetX()>90.)new (¶md) AliExternalTrackParam(ktrack1->GetReference());
4850 kink->SetMother(paramm);
4851 kink->SetDaughter(paramd);
4854 Float_t x[3] = { kink->GetPosition()[0],kink->GetPosition()[1],kink->GetPosition()[2]};
4856 fParam->Transform0to1(x,index);
4857 fParam->Transform1to2(x,index);
4858 row0 = GetRowNumber(x[0]);
4860 if (kink->GetR()<100) continue;
4861 if (kink->GetR()>240) continue;
4862 if (kink->GetPosition()[2]/kink->GetR()>AliTPCReconstructor::GetCtgRange()) continue; //out of fiducial volume
4863 if (kink->GetDistance()>cdist3) continue;
4864 Float_t dird = kink->GetDaughterP()[0]*kink->GetPosition()[0]+kink->GetDaughterP()[1]*kink->GetPosition()[1]; // rough direction estimate
4865 if (dird<0) continue;
4867 Float_t dirm = kink->GetMotherP()[0]*kink->GetPosition()[0]+kink->GetMotherP()[1]*kink->GetPosition()[1]; // rough direction estimate
4868 if (dirm<0) continue;
4869 Float_t mpt = TMath::Sqrt(kink->GetMotherP()[0]*kink->GetMotherP()[0]+kink->GetMotherP()[1]*kink->GetMotherP()[1]);
4870 if (mpt<0.2) continue;
4873 //for high momenta momentum not defined well in first iteration
4874 Double_t qt = TMath::Sin(kink->GetAngle(2))*ktrack1->GetP();
4875 if (qt>0.35) continue;
4878 kink->SetLabel(CookLabel(ktrack0,0.4,0,row0),0);
4879 kink->SetLabel(CookLabel(ktrack1,0.4,row0,160),1);
4881 kink->SetTPCDensity(dens00,0,0);
4882 kink->SetTPCDensity(dens01,0,1);
4883 kink->SetTPCDensity(dens10,1,0);
4884 kink->SetTPCDensity(dens11,1,1);
4885 kink->SetIndex(i,0);
4886 kink->SetIndex(j,1);
4889 kink->SetTPCDensity(dens10,0,0);
4890 kink->SetTPCDensity(dens11,0,1);
4891 kink->SetTPCDensity(dens00,1,0);
4892 kink->SetTPCDensity(dens01,1,1);
4893 kink->SetIndex(j,0);
4894 kink->SetIndex(i,1);
4897 if (mpt<1||kink->GetAngle(2)>0.1){
4898 // angle and densities not defined yet
4899 if (kink->GetTPCDensityFactor()<0.8) continue;
4900 if ((2-kink->GetTPCDensityFactor())*kink->GetDistance() >0.25) continue;
4901 if (kink->GetAngle(2)*ktrack0->GetP()<0.003) continue; //too small angle
4902 if (kink->GetAngle(2)>0.2&&kink->GetTPCDensityFactor()<1.15) continue;
4903 if (kink->GetAngle(2)>0.2&&kink->GetTPCDensity(0,1)>0.05) continue;
4905 Float_t criticalangle = track0->GetSigmaSnp2()+track0->GetSigmaTgl2();
4906 criticalangle+= track1->GetSigmaSnp2()+track1->GetSigmaTgl2();
4907 criticalangle= 3*TMath::Sqrt(criticalangle);
4908 if (criticalangle>0.02) criticalangle=0.02;
4909 if (kink->GetAngle(2)<criticalangle) continue;
4912 Int_t drow = Int_t(2.+0.5/(0.05+kink->GetAngle(2))); // overlap region defined
4913 Float_t shapesum =0;
4915 for ( Int_t row = row0-drow; row<row0+drow;row++){
4916 if (row<0) continue;
4917 if (row>155) continue;
4918 if (ktrack0->GetClusterPointer(row)){
4919 AliTPCTrackerPoint *point =ktrack0->GetTrackPoint(row);
4920 shapesum+=point->GetSigmaY()+point->GetSigmaZ();
4923 if (ktrack1->GetClusterPointer(row)){
4924 AliTPCTrackerPoint *point =ktrack1->GetTrackPoint(row);
4925 shapesum+=point->GetSigmaY()+point->GetSigmaZ();
4930 kink->SetShapeFactor(-1.);
4933 kink->SetShapeFactor(shapesum/sum);
4935 // esd->AddKink(kink);
4936 kinks->AddLast(kink);
4942 // sort the kinks according quality - and refit them towards vertex
4944 Int_t nkinks = kinks->GetEntriesFast();
4945 Float_t *quality = new Float_t[nkinks];
4946 Int_t *indexes = new Int_t[nkinks];
4947 AliTPCseed *mothers = new AliTPCseed[nkinks];
4948 AliTPCseed *daughters = new AliTPCseed[nkinks];
4951 for (Int_t i=0;i<nkinks;i++){
4953 AliKink *kink = (AliKink*)kinks->At(i);
4955 // refit kinks towards vertex
4957 Int_t index0 = kink->GetIndex(0);
4958 Int_t index1 = kink->GetIndex(1);
4959 AliTPCseed * ktrack0 = (AliTPCseed*)array->At(index0);
4960 AliTPCseed * ktrack1 = (AliTPCseed*)array->At(index1);
4962 Int_t sumn=ktrack0->GetNumberOfClusters()+ktrack1->GetNumberOfClusters();
4964 // Refit Kink under if too small angle
4966 if (kink->GetAngle(2)<0.05){
4967 kink->SetTPCRow0(GetRowNumber(kink->GetR()));
4968 Int_t row0 = kink->GetTPCRow0();
4969 Int_t drow = Int_t(2.+0.5/(0.05+kink->GetAngle(2)));
4972 Int_t last = row0-drow;
4973 if (last<40) last=40;
4974 if (last<ktrack0->GetFirstPoint()+25) last = ktrack0->GetFirstPoint()+25;
4975 AliTPCseed* seed0 = ReSeed(ktrack0,last,kFALSE);
4978 Int_t first = row0+drow;
4979 if (first>130) first=130;
4980 if (first>ktrack1->GetLastPoint()-25) first = TMath::Max(ktrack1->GetLastPoint()-25,30);
4981 AliTPCseed* seed1 = ReSeed(ktrack1,first,kTRUE);
4983 if (seed0 && seed1){
4984 kink->SetStatus(1,8);
4985 if (RefitKink(*seed0,*seed1,*kink)) kink->SetStatus(1,9);
4986 row0 = GetRowNumber(kink->GetR());
4987 sumn = seed0->GetNumberOfClusters()+seed1->GetNumberOfClusters();
4988 mothers[i] = *seed0;
4989 daughters[i] = *seed1;
4992 delete kinks->RemoveAt(i);
4993 if (seed0) delete seed0;
4994 if (seed1) delete seed1;
4997 if (kink->GetDistance()>0.5 || kink->GetR()<110 || kink->GetR()>240) {
4998 delete kinks->RemoveAt(i);
4999 if (seed0) delete seed0;
5000 if (seed1) delete seed1;
5008 if (kink) quality[i] = 160*((0.1+kink->GetDistance())*(2.-kink->GetTPCDensityFactor()))/(sumn+40.); //the longest -clossest will win
5010 TMath::Sort(nkinks,quality,indexes,kFALSE);
5012 //remove double find kinks
5014 for (Int_t ikink0=1;ikink0<nkinks;ikink0++){
5015 AliKink * kink0 = (AliKink*) kinks->At(indexes[ikink0]);
5016 if (!kink0) continue;
5018 for (Int_t ikink1=0;ikink1<ikink0;ikink1++){
5019 if (!kink0) continue;
5020 AliKink * kink1 = (AliKink*) kinks->At(indexes[ikink1]);
5021 if (!kink1) continue;
5022 // if not close kink continue
5023 if (TMath::Abs(kink1->GetPosition()[2]-kink0->GetPosition()[2])>10) continue;
5024 if (TMath::Abs(kink1->GetPosition()[1]-kink0->GetPosition()[1])>10) continue;
5025 if (TMath::Abs(kink1->GetPosition()[0]-kink0->GetPosition()[0])>10) continue;
5027 AliTPCseed &mother0 = mothers[indexes[ikink0]];
5028 AliTPCseed &daughter0 = daughters[indexes[ikink0]];
5029 AliTPCseed &mother1 = mothers[indexes[ikink1]];
5030 AliTPCseed &daughter1 = daughters[indexes[ikink1]];
5031 Int_t row0 = (kink0->GetTPCRow0()+kink1->GetTPCRow0())/2;
5040 for (Int_t i=0;i<row0;i++){
5041 if (mother0.GetClusterIndex(i)>0 && mother1.GetClusterIndex(i)>0){
5044 if (mother0.GetClusterIndex(i)==mother1.GetClusterIndex(i)){
5051 for (Int_t i=row0;i<158;i++){
5052 if (daughter0.GetClusterIndex(i)>0 && daughter0.GetClusterIndex(i)>0){
5055 if (mother0.GetClusterIndex(i)==mother1.GetClusterIndex(i)){
5061 Float_t ratio = Float_t(same+1)/Float_t(both+1);
5062 Float_t ratiom = Float_t(samem+1)/Float_t(bothm+1);
5063 Float_t ratiod = Float_t(samed+1)/Float_t(bothd+1);
5064 if (ratio>0.3 && ratiom>0.5 &&ratiod>0.5) {
5065 Int_t sum0 = mother0.GetNumberOfClusters()+daughter0.GetNumberOfClusters();
5066 Int_t sum1 = mother1.GetNumberOfClusters()+daughter1.GetNumberOfClusters();
5068 shared[kink0->GetIndex(0)]= kTRUE;
5069 shared[kink0->GetIndex(1)]= kTRUE;
5070 delete kinks->RemoveAt(indexes[ikink0]);
5073 shared[kink1->GetIndex(0)]= kTRUE;
5074 shared[kink1->GetIndex(1)]= kTRUE;
5075 delete kinks->RemoveAt(indexes[ikink1]);
5082 for (Int_t i=0;i<nkinks;i++){
5083 AliKink * kink = (AliKink*) kinks->At(indexes[i]);
5084 if (!kink) continue;
5085 kink->SetTPCRow0(GetRowNumber(kink->GetR()));
5086 Int_t index0 = kink->GetIndex(0);
5087 Int_t index1 = kink->GetIndex(1);
5088 if (circular[index0]||circular[index1]&&kink->GetDistance()>0.2) continue;
5089 kink->SetMultiple(usage[index0],0);
5090 kink->SetMultiple(usage[index1],1);
5091 if (kink->GetMultiple()[0]+kink->GetMultiple()[1]>2) continue;
5092 if (kink->GetMultiple()[0]+kink->GetMultiple()[1]>0 && quality[indexes[i]]>0.2) continue;
5093 if (kink->GetMultiple()[0]+kink->GetMultiple()[1]>0 && kink->GetDistance()>0.2) continue;
5094 if (circular[index0]||circular[index1]&&kink->GetDistance()>0.1) continue;
5096 AliTPCseed * ktrack0 = (AliTPCseed*)array->At(index0);
5097 AliTPCseed * ktrack1 = (AliTPCseed*)array->At(index1);
5098 if (!ktrack0 || !ktrack1) continue;
5099 Int_t index = esd->AddKink(kink);
5102 if ( ktrack0->GetKinkIndex(0)==0 && ktrack1->GetKinkIndex(0)==0) { //best kink
5103 if (mothers[indexes[i]].GetNumberOfClusters()>20 && daughters[indexes[i]].GetNumberOfClusters()>20 && (mothers[indexes[i]].GetNumberOfClusters()+daughters[indexes[i]].GetNumberOfClusters())>100){
5104 *ktrack0 = mothers[indexes[i]];
5105 *ktrack1 = daughters[indexes[i]];
5109 ktrack0->SetKinkIndex(usage[index0],-(index+1));
5110 ktrack1->SetKinkIndex(usage[index1], (index+1));
5115 // Remove tracks corresponding to shared kink's
5117 for (Int_t i=0;i<nentries;i++){
5118 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5119 if (!track0) continue;
5120 if (track0->GetKinkIndex(0)!=0) continue;
5121 if (shared[i]) delete array->RemoveAt(i);
5126 RemoveUsed2(array,0.5,0.4,30);
5128 for (Int_t i=0;i<nentries;i++){
5129 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5130 if (!track0) continue;
5131 track0->CookdEdx(0.02,0.6);
5135 for (Int_t i=0;i<nentries;i++){
5136 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5137 if (!track0) continue;
5138 if (track0->Pt()<1.4) continue;
5139 //remove double high momenta tracks - overlapped with kink candidates
5142 for (Int_t icl=track0->GetFirstPoint();icl<track0->GetLastPoint(); icl++){
5143 if (track0->GetClusterPointer(icl)!=0){
5145 if (track0->GetClusterPointer(icl)->IsUsed(10)) shared++;
5148 if (Float_t(shared+1)/Float_t(all+1)>0.5) {
5149 delete array->RemoveAt(i);
5153 if (track0->GetKinkIndex(0)!=0) continue;
5154 if (track0->GetNumberOfClusters()<80) continue;
5156 AliTPCseed *pmother = new AliTPCseed();
5157 AliTPCseed *pdaughter = new AliTPCseed();
5158 AliKink *pkink = new AliKink;
5160 AliTPCseed & mother = *pmother;
5161 AliTPCseed & daughter = *pdaughter;
5162 AliKink & kink = *pkink;
5163 if (CheckKinkPoint(track0,mother,daughter, kink)){
5164 if (mother.GetNumberOfClusters()<30||daughter.GetNumberOfClusters()<20) {
5168 continue; //too short tracks
5170 if (mother.Pt()<1.4) {
5176 Int_t row0= kink.GetTPCRow0();
5177 if (kink.GetDistance()>0.5 || kink.GetR()<110. || kink.GetR()>240.) {
5184 Int_t index = esd->AddKink(&kink);
5185 mother.SetKinkIndex(0,-(index+1));
5186 daughter.SetKinkIndex(0,index+1);
5187 if (mother.GetNumberOfClusters()>50) {
5188 delete array->RemoveAt(i);
5189 array->AddAt(new AliTPCseed(mother),i);
5192 array->AddLast(new AliTPCseed(mother));
5194 array->AddLast(new AliTPCseed(daughter));
5195 for (Int_t icl=0;icl<row0;icl++) {
5196 if (mother.GetClusterPointer(icl)) mother.GetClusterPointer(icl)->Use(20);
5199 for (Int_t icl=row0;icl<158;icl++) {
5200 if (daughter.GetClusterPointer(icl)) daughter.GetClusterPointer(icl)->Use(20);
5209 delete [] daughters;
5231 printf("Ncandidates=\t%d\t%d\t%d\t%d\n",esd->GetNumberOfKinks(),ncandidates,ntracks,nall);
5235 void AliTPCtrackerMI::FindV0s(TObjArray * array, AliESDEvent *esd)
5241 TObjArray *tpcv0s = new TObjArray(100000);
5242 Int_t nentries = array->GetEntriesFast();
5243 AliHelix *helixes = new AliHelix[nentries];
5244 Int_t *sign = new Int_t[nentries];
5245 Float_t *alpha = new Float_t[nentries];
5246 Float_t *z0 = new Float_t[nentries];
5247 Float_t *dca = new Float_t[nentries];
5248 Float_t *sdcar = new Float_t[nentries];
5249 Float_t *cdcar = new Float_t[nentries];
5250 Float_t *pulldcar = new Float_t[nentries];
5251 Float_t *pulldcaz = new Float_t[nentries];
5252 Float_t *pulldca = new Float_t[nentries];
5253 Bool_t *isPrim = new Bool_t[nentries];
5254 const AliESDVertex * primvertex = esd->GetVertex();
5255 Double_t zvertex = primvertex->GetZv();
5257 // nentries = array->GetEntriesFast();
5259 for (Int_t i=0;i<nentries;i++){
5262 AliTPCseed* track = (AliTPCseed*)array->At(i);
5263 if (!track) continue;
5264 track->GetV0Indexes()[0] = 0; //rest v0 indexes
5265 track->GetV0Indexes()[1] = 0; //rest v0 indexes
5266 track->GetV0Indexes()[2] = 0; //rest v0 indexes
5268 alpha[i] = track->GetAlpha();
5269 new (&helixes[i]) AliHelix(*track);
5271 helixes[i].Evaluate(0,xyz);
5272 sign[i] = (track->GetC()>0) ? -1:1;
5276 if (track->GetProlongation(0,y,z)) z0[i] = z;
5277 dca[i] = track->GetD(0,0);
5279 // dca error parrameterezation + pulls
5281 sdcar[i] = TMath::Sqrt(0.150*0.150+(100*track->GetC())*(100*track->GetC()));
5282 if (TMath::Abs(track->GetTgl())>1) sdcar[i]*=2.5;
5283 cdcar[i] = TMath::Exp((TMath::Abs(track->GetC())-0.0106)*525.3);
5284 pulldcar[i] = (dca[i]-cdcar[i])/sdcar[i];
5285 pulldcaz[i] = (z0[i]-zvertex)/sdcar[i];
5286 pulldca[i] = TMath::Sqrt(pulldcar[i]*pulldcar[i]+pulldcaz[i]*pulldcaz[i]);
5287 if (track->TPCrPID(1)+track->TPCrPID(2)+track->TPCrPID(3)>0.5) {
5288 if (pulldca[i]<3.) isPrim[i]=kTRUE; //pion, muon and Kaon 3 sigma cut
5290 if (track->TPCrPID(4)>0.5) {
5291 if (pulldca[i]<0.5) isPrim[i]=kTRUE; //proton 0.5 sigma cut
5293 if (track->TPCrPID(0)>0.4) {
5294 isPrim[i]=kFALSE; //electron no sigma cut
5301 Int_t ncandidates =0;
5304 Double_t phase[2][2],radius[2];
5310 TTreeSRedirector &cstream = *fDebugStreamer;
5311 Float_t fprimvertex[3]={GetX(),GetY(),GetZ()};
5313 Double_t cradius0 = 10*10;
5314 Double_t cradius1 = 200*200;
5317 Double_t cpointAngle = 0.95;
5319 Double_t delta[2]={10000,10000};
5320 for (Int_t i =0;i<nentries;i++){
5321 if (sign[i]==0) continue;
5322 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5323 if (!track0) continue;
5324 if (AliTPCReconstructor::StreamLevel()>1){
5329 "zvertex="<<zvertex<<
5330 "sdcar0="<<sdcar[i]<<
5331 "cdcar0="<<cdcar[i]<<
5332 "pulldcar0="<<pulldcar[i]<<
5333 "pulldcaz0="<<pulldcaz[i]<<
5334 "pulldca0="<<pulldca[i]<<
5335 "isPrim="<<isPrim[i]<<
5339 if (track0->GetSigned1Pt()<0) continue;
5340 if (track0->GetKinkIndex(0)>0||isPrim[i]) continue; //daughter kink
5342 if (TMath::Abs(helixes[i].GetHelix(4))<0.000000001) continue;
5347 for (Int_t j =0;j<nentries;j++){
5348 AliTPCseed * track1 = (AliTPCseed*)array->At(j);
5349 if (!track1) continue;
5350 if (track1->GetKinkIndex(0)>0 || isPrim[j]) continue; //daughter kink
5351 if (sign[j]*sign[i]>0) continue;
5352 if (TMath::Abs(helixes[j].GetHelix(4))<0.000001) continue;
5353 if (track0->GetCircular()+track1->GetCircular()>1) continue; //circling -returning track
5356 // DCA to prim vertex cut
5362 Int_t npoints = helixes[i].GetRPHIintersections(helixes[j], phase, radius,cdist2);
5363 if (npoints<1) continue;
5367 if (radius[0]<cradius0||radius[0]>cradius1) continue;
5368 helixes[i].LinearDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta[0]);
5369 if (delta[0]>cdist1) continue;
5372 if (TMath::Max(radius[0],radius[1])<cradius0|| TMath::Min(radius[0],radius[1])>cradius1) continue;
5373 helixes[i].LinearDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta[0]);
5374 helixes[i].LinearDCA(helixes[j],phase[1][0],phase[1][1],radius[1],delta[1]);
5375 if (delta[1]<delta[0]) iclosest=1;
5376 if (delta[iclosest]>cdist1) continue;
5378 helixes[i].ParabolicDCA(helixes[j],phase[iclosest][0],phase[iclosest][1],radius[iclosest],delta[iclosest]);
5379 if (radius[iclosest]<cradius0 || radius[iclosest]>cradius1 || delta[iclosest]>cdist1) continue;
5381 Double_t pointAngle = helixes[i].GetPointAngle(helixes[j],phase[iclosest],fprimvertex);
5382 if (pointAngle<cpointAngle) continue;
5384 Bool_t isGamma = kFALSE;
5385 vertex.SetParamP(*track0); //track0 - plus
5386 vertex.SetParamN(*track1); //track1 - minus
5387 vertex.Update(fprimvertex);
5388 if (track0->TPCrPID(0)>0.3&&track1->TPCrPID(0)>0.3&&vertex.GetAnglep()[2]<0.15) isGamma=kTRUE; // gamma conversion candidate
5389 Double_t pointAngle2 = vertex.GetV0CosineOfPointingAngle();
5391 if (vertex.GetV0CosineOfPointingAngle()<cpointAngle && (!isGamma)) continue;// point angle cut
5392 //Bo: if (vertex.GetDist2()>2&&(!isGamma)) continue; // point angle cut
5393 if (vertex.GetDcaV0Daughters()>2&&(!isGamma)) continue;//Bo: // point angle cut
5394 Float_t sigmae = 0.15*0.15;
5395 if (vertex.GetRr()<80)
5396 sigmae += (sdcar[i]*sdcar[i]+sdcar[j]*sdcar[j])*(1.-vertex.GetRr()/80.)*(1.-vertex.GetRr()/80.);
5397 sigmae+= TMath::Sqrt(sigmae);
5398 //Bo: if (vertex.GetDist2()/sigmae>3.&&(!isGamma)) continue;
5399 if (vertex.GetDcaV0Daughters()/sigmae>3.&&(!isGamma)) continue;
5400 Float_t densb0=0,densb1=0,densa0=0,densa1=0;
5401 Int_t row0 = GetRowNumber(vertex.GetRr());
5403 //Bo: if (vertex.GetDist2()>0.2) continue;
5404 if (vertex.GetDcaV0Daughters()>0.2) continue;
5405 densb0 = track0->Density2(0,row0-5);
5406 densb1 = track1->Density2(0,row0-5);
5407 if (densb0>0.3|| densb1>0.3) continue; //clusters before vertex
5408 densa0 = track0->Density2(row0+5,row0+40);
5409 densa1 = track1->Density2(row0+5,row0+40);
5410 if ((densa0<0.4|| densa1<0.4)&&(!isGamma)) continue; //missing clusters after vertex
5413 densa0 = track0->Density2(0,40); //cluster density
5414 densa1 = track1->Density2(0,40); //cluster density
5415 if ((vertex.GetRr()<80&&densa0+densa1<1.)&&(!isGamma)) continue;
5417 //Bo: vertex.SetLab(0,track0->GetLabel());
5418 //Bo: vertex.SetLab(1,track1->GetLabel());
5419 vertex.SetChi2After((densa0+densa1)*0.5);
5420 vertex.SetChi2Before((densb0+densb1)*0.5);
5421 vertex.SetIndex(0,i);
5422 vertex.SetIndex(1,j);
5423 //Bo: vertex.SetStatus(1); // TPC v0 candidate
5424 vertex.SetOnFlyStatus(2);//Bo: // TPC v0 candidate
5425 //Bo: vertex.SetRp(track0->TPCrPIDs());
5426 //Bo: vertex.SetRm(track1->TPCrPIDs());
5427 tpcv0s->AddLast(new AliESDv0(vertex));
5430 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
5431 Double_t radiusm= (delta[0]<delta[1])? TMath::Sqrt(radius[0]):TMath::Sqrt(radius[1]);
5432 Double_t deltam= (delta[0]<delta[1])? TMath::Sqrt(delta[0]):TMath::Sqrt(delta[1]);
5433 if (AliTPCReconstructor::StreamLevel()>1) {
5434 Int_t lab0=track0->GetLabel();
5435 Int_t lab1=track1->GetLabel();
5436 Char_t c0=track0->GetCircular();
5437 Char_t c1=track1->GetCircular();
5440 "vertex.="<<&vertex<<
5443 "Helix0.="<<&helixes[i]<<
5446 "Helix1.="<<&helixes[j]<<
5447 "pointAngle="<<pointAngle<<
5448 "pointAngle2="<<pointAngle2<<
5453 "zvertex="<<zvertex<<
5456 "npoints="<<npoints<<
5457 "radius0="<<radius[0]<<
5458 "delta0="<<delta[0]<<
5459 "radius1="<<radius[1]<<
5460 "delta1="<<delta[1]<<
5461 "radiusm="<<radiusm<<
5463 "sdcar0="<<sdcar[i]<<
5464 "sdcar1="<<sdcar[j]<<
5465 "cdcar0="<<cdcar[i]<<
5466 "cdcar1="<<cdcar[j]<<
5467 "pulldcar0="<<pulldcar[i]<<
5468 "pulldcar1="<<pulldcar[j]<<
5469 "pulldcaz0="<<pulldcaz[i]<<
5470 "pulldcaz1="<<pulldcaz[j]<<
5471 "pulldca0="<<pulldca[i]<<
5472 "pulldca1="<<pulldca[j]<<
5482 Float_t *quality = new Float_t[ncandidates];
5483 Int_t *indexes = new Int_t[ncandidates];
5485 for (Int_t i=0;i<ncandidates;i++){
5487 AliESDv0 *v0 = (AliESDv0*)tpcv0s->At(i);
5488 quality[i] = 1./(1.00001-v0->GetV0CosineOfPointingAngle()); //base point angle
5489 // quality[i] /= (0.5+v0->GetDist2());
5490 // quality[i] *= v0->GetChi2After(); //density factor
5492 Int_t index0 = v0->GetIndex(0);
5493 Int_t index1 = v0->GetIndex(1);
5494 //Bo: Double_t minpulldca = TMath::Min(2.+pulldca[v0->GetIndex(0)],(2.+pulldca[v0->GetIndex(1)]) ); //pull
5495 Double_t minpulldca = TMath::Min(2.+pulldca[index0],(2.+pulldca[index1]) );//Bo:
5499 AliTPCseed * track0 = (AliTPCseed*)array->At(index0);
5500 AliTPCseed * track1 = (AliTPCseed*)array->At(index1);
5501 if (track0->TPCrPID(0)>0.3&&track1->TPCrPID(0)>0.3&&v0->GetAnglep()[2]<0.15) quality[i]+=1000000; // gamma conversion candidate
5502 if (track0->TPCrPID(4)>0.9||track1->TPCrPID(4)>0.9&&minpulldca>4) quality[i]*=10; // lambda candidate candidate
5505 TMath::Sort(ncandidates,quality,indexes,kTRUE);
5508 for (Int_t i=0;i<ncandidates;i++){
5509 AliESDv0 * v0 = (AliESDv0*)tpcv0s->At(indexes[i]);
5511 Int_t index0 = v0->GetIndex(0);
5512 Int_t index1 = v0->GetIndex(1);
5513 AliTPCseed * track0 = (AliTPCseed*)array->At(index0);
5514 AliTPCseed * track1 = (AliTPCseed*)array->At(index1);
5515 if (!track0||!track1) {
5519 Bool_t accept =kTRUE; //default accept
5520 Int_t *v0indexes0 = track0->GetV0Indexes();
5521 Int_t *v0indexes1 = track1->GetV0Indexes();
5523 Int_t order0 = (v0indexes0[0]!=0) ? 1:0;
5524 Int_t order1 = (v0indexes1[0]!=0) ? 1:0;
5525 if (v0indexes0[1]!=0) order0 =2;
5526 if (v0indexes1[1]!=0) order1 =2;
5528 if (v0indexes0[2]!=0) {order0=3; accept=kFALSE;}
5529 if (v0indexes0[2]!=0) {order1=3; accept=kFALSE;}
5531 AliESDv0 * v02 = v0;
5533 //Bo: v0->SetOrder(0,order0);
5534 //Bo: v0->SetOrder(1,order1);
5535 //Bo: v0->SetOrder(1,order0+order1);
5536 v0->SetOnFlyStatus(kTRUE);
5537 Int_t index = esd->AddV0(v0);
5538 v02 = esd->GetV0(index);
5539 v0indexes0[order0]=index;
5540 v0indexes1[order1]=index;
5544 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
5545 if (AliTPCReconstructor::StreamLevel()>1) {
5546 Int_t lab0=track0->GetLabel();
5547 Int_t lab1=track1->GetLabel();
5556 "dca0="<<dca[index0]<<
5557 "dca1="<<dca[index1]<<
5561 "quality="<<quality[i]<<
5562 "pulldca0="<<pulldca[index0]<<
5563 "pulldca1="<<pulldca[index1]<<
5587 printf("TPC V0 finder : naccepted\t%d\tncandidates\t%d\tntracks\t%d\tnall\t%d\n",naccepted,ncandidates,ntracks,nall);
5591 Int_t AliTPCtrackerMI::RefitKink(AliTPCseed &mother, AliTPCseed &daughter, AliESDkink &knk)
5594 // refit kink towards to the vertex
5597 AliKink &kink=(AliKink &)knk;
5599 Int_t row0 = GetRowNumber(kink.GetR());
5600 FollowProlongation(mother,0);
5601 mother.Reset(kFALSE);
5603 FollowProlongation(daughter,row0);
5604 daughter.Reset(kFALSE);
5605 FollowBackProlongation(daughter,158);
5606 daughter.Reset(kFALSE);
5607 Int_t first = TMath::Max(row0-20,30);
5608 Int_t last = TMath::Min(row0+20,140);
5610 const Int_t kNdiv =5;
5611 AliTPCseed param0[kNdiv]; // parameters along the track
5612 AliTPCseed param1[kNdiv]; // parameters along the track
5613 AliKink kinks[kNdiv]; // corresponding kink parameters
5616 for (Int_t irow=0; irow<kNdiv;irow++){
5617 rows[irow] = first +((last-first)*irow)/(kNdiv-1);
5619 // store parameters along the track
5621 for (Int_t irow=0;irow<kNdiv;irow++){
5622 FollowBackProlongation(mother, rows[irow]);
5623 FollowProlongation(daughter,rows[kNdiv-1-irow]);
5624 param0[irow] = mother;
5625 param1[kNdiv-1-irow] = daughter;
5629 for (Int_t irow=0; irow<kNdiv-1;irow++){
5630 if (param0[irow].GetNumberOfClusters()<kNdiv||param1[irow].GetNumberOfClusters()<kNdiv) continue;
5631 kinks[irow].SetMother(param0[irow]);
5632 kinks[irow].SetDaughter(param1[irow]);
5633 kinks[irow].Update();
5636 // choose kink with best "quality"
5638 Double_t mindist = 10000;
5639 for (Int_t irow=0;irow<kNdiv;irow++){
5640 if (param0[irow].GetNumberOfClusters()<20||param1[irow].GetNumberOfClusters()<20) continue;
5641 if (TMath::Abs(kinks[irow].GetR())>240.) continue;
5642 if (TMath::Abs(kinks[irow].GetR())<100.) continue;
5644 Float_t normdist = TMath::Abs(param0[irow].GetX()-kinks[irow].GetR())*(0.1+kink.GetDistance());
5645 normdist/= (param0[irow].GetNumberOfClusters()+param1[irow].GetNumberOfClusters()+40.);
5646 if (normdist < mindist){
5652 if (index==-1) return 0;
5655 param0[index].Reset(kTRUE);
5656 FollowProlongation(param0[index],0);
5658 mother = param0[index];
5659 daughter = param1[index]; // daughter in vertex
5661 kink.SetMother(mother);
5662 kink.SetDaughter(daughter);
5664 kink.SetTPCRow0(GetRowNumber(kink.GetR()));
5665 kink.SetTPCncls(param0[index].GetNumberOfClusters(),0);
5666 kink.SetTPCncls(param1[index].GetNumberOfClusters(),1);
5667 kink.SetLabel(CookLabel(&mother,0.4, 0,kink.GetTPCRow0()),0);
5668 kink.SetLabel(CookLabel(&daughter,0.4, kink.GetTPCRow0(),160),1);
5669 mother.SetLabel(kink.GetLabel(0));
5670 daughter.SetLabel(kink.GetLabel(1));
5676 void AliTPCtrackerMI::UpdateKinkQualityM(AliTPCseed * seed){
5678 // update Kink quality information for mother after back propagation
5680 if (seed->GetKinkIndex(0)>=0) return;
5681 for (Int_t ikink=0;ikink<3;ikink++){
5682 Int_t index = seed->GetKinkIndex(ikink);
5683 if (index>=0) break;
5684 index = TMath::Abs(index)-1;
5685 AliESDkink * kink = fEvent->GetKink(index);
5686 kink->SetTPCDensity(-1,0,0);
5687 kink->SetTPCDensity(1,0,1);
5689 Int_t row0 = kink->GetTPCRow0() - 2 - Int_t( 0.5/ (0.05+kink->GetAngle(2)));
5690 if (row0<15) row0=15;
5692 Int_t row1 = kink->GetTPCRow0() + 2 + Int_t( 0.5/ (0.05+kink->GetAngle(2)));
5693 if (row1>145) row1=145;
5695 Int_t found,foundable,shared;
5696 seed->GetClusterStatistic(0,row0, found, foundable,shared,kFALSE);
5697 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),0,0);
5698 seed->GetClusterStatistic(row1,155, found, foundable,shared,kFALSE);
5699 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),0,1);
5704 void AliTPCtrackerMI::UpdateKinkQualityD(AliTPCseed * seed){
5706 // update Kink quality information for daughter after refit
5708 if (seed->GetKinkIndex(0)<=0) return;
5709 for (Int_t ikink=0;ikink<3;ikink++){
5710 Int_t index = seed->GetKinkIndex(ikink);
5711 if (index<=0) break;
5712 index = TMath::Abs(index)-1;
5713 AliESDkink * kink = fEvent->GetKink(index);
5714 kink->SetTPCDensity(-1,1,0);
5715 kink->SetTPCDensity(-1,1,1);
5717 Int_t row0 = kink->GetTPCRow0() -2 - Int_t( 0.5/ (0.05+kink->GetAngle(2)));
5718 if (row0<15) row0=15;
5720 Int_t row1 = kink->GetTPCRow0() +2 + Int_t( 0.5/ (0.05+kink->GetAngle(2)));
5721 if (row1>145) row1=145;
5723 Int_t found,foundable,shared;
5724 seed->GetClusterStatistic(0,row0, found, foundable,shared,kFALSE);
5725 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),1,0);
5726 seed->GetClusterStatistic(row1,155, found, foundable,shared,kFALSE);
5727 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),1,1);
5733 Int_t AliTPCtrackerMI::CheckKinkPoint(AliTPCseed*seed,AliTPCseed &mother, AliTPCseed &daughter, AliESDkink &knk)
5736 // check kink point for given track
5737 // if return value=0 kink point not found
5738 // otherwise seed0 correspond to mother particle
5739 // seed1 correspond to daughter particle
5740 // kink parameter of kink point
5741 AliKink &kink=(AliKink &)knk;
5743 Int_t middlerow = (seed->GetFirstPoint()+seed->GetLastPoint())/2;
5744 Int_t first = seed->GetFirstPoint();
5745 Int_t last = seed->GetLastPoint();
5746 if (last-first<20) return 0; // shortest length - 2*30 = 60 pad-rows
5749 AliTPCseed *seed1 = ReSeed(seed,middlerow+20, kTRUE); //middle of chamber
5750 if (!seed1) return 0;
5751 FollowProlongation(*seed1,seed->GetLastPoint()-20);
5752 seed1->Reset(kTRUE);
5753 FollowProlongation(*seed1,158);
5754 seed1->Reset(kTRUE);
5755 last = seed1->GetLastPoint();
5757 AliTPCseed *seed0 = new AliTPCseed(*seed);
5758 seed0->Reset(kFALSE);
5761 AliTPCseed param0[20]; // parameters along the track
5762 AliTPCseed param1[20]; // parameters along the track
5763 AliKink kinks[20]; // corresponding kink parameters
5765 for (Int_t irow=0; irow<20;irow++){
5766 rows[irow] = first +((last-first)*irow)/19;
5768 // store parameters along the track
5770 for (Int_t irow=0;irow<20;irow++){
5771 FollowBackProlongation(*seed0, rows[irow]);
5772 FollowProlongation(*seed1,rows[19-irow]);
5773 param0[irow] = *seed0;
5774 param1[19-irow] = *seed1;
5778 for (Int_t irow=0; irow<19;irow++){
5779 kinks[irow].SetMother(param0[irow]);
5780 kinks[irow].SetDaughter(param1[irow]);
5781 kinks[irow].Update();
5784 // choose kink with biggest change of angle
5786 Double_t maxchange= 0;
5787 for (Int_t irow=1;irow<19;irow++){
5788 if (TMath::Abs(kinks[irow].GetR())>240.) continue;
5789 if (TMath::Abs(kinks[irow].GetR())<110.) continue;
5790 Float_t quality = TMath::Abs(kinks[irow].GetAngle(2))/(3.+TMath::Abs(kinks[irow].GetR()-param0[irow].GetX()));
5791 if ( quality > maxchange){
5792 maxchange = quality;
5799 if (index<0) return 0;
5801 Int_t row0 = GetRowNumber(kinks[index].GetR()); //row 0 estimate
5802 seed0 = new AliTPCseed(param0[index]);
5803 seed1 = new AliTPCseed(param1[index]);
5804 seed0->Reset(kFALSE);
5805 seed1->Reset(kFALSE);
5806 seed0->ResetCovariance(10.);
5807 seed1->ResetCovariance(10.);
5808 FollowProlongation(*seed0,0);
5809 FollowBackProlongation(*seed1,158);
5810 mother = *seed0; // backup mother at position 0
5811 seed0->Reset(kFALSE);
5812 seed1->Reset(kFALSE);
5813 seed0->ResetCovariance(10.);
5814 seed1->ResetCovariance(10.);
5816 first = TMath::Max(row0-20,0);
5817 last = TMath::Min(row0+20,158);
5819 for (Int_t irow=0; irow<20;irow++){
5820 rows[irow] = first +((last-first)*irow)/19;
5822 // store parameters along the track
5824 for (Int_t irow=0;irow<20;irow++){
5825 FollowBackProlongation(*seed0, rows[irow]);
5826 FollowProlongation(*seed1,rows[19-irow]);
5827 param0[irow] = *seed0;
5828 param1[19-irow] = *seed1;
5832 for (Int_t irow=0; irow<19;irow++){
5833 kinks[irow].SetMother(param0[irow]);
5834 kinks[irow].SetDaughter(param1[irow]);
5835 // param0[irow].Dump();
5836 //param1[irow].Dump();
5837 kinks[irow].Update();
5840 // choose kink with biggest change of angle
5843 for (Int_t irow=0;irow<20;irow++){
5844 if (TMath::Abs(kinks[irow].GetR())>250.) continue;
5845 if (TMath::Abs(kinks[irow].GetR())<90.) continue;
5846 Float_t quality = TMath::Abs(kinks[irow].GetAngle(2))/(3.+TMath::Abs(kinks[irow].GetR()-param0[irow].GetX()));
5847 if ( quality > maxchange){
5848 maxchange = quality;
5855 if (index==-1 || param0[index].GetNumberOfClusters()+param1[index].GetNumberOfClusters()<100){
5860 // Float_t anglesigma = TMath::Sqrt(param0[index].fC22+param0[index].fC33+param1[index].fC22+param1[index].fC33);
5862 kink.SetMother(param0[index]);
5863 kink.SetDaughter(param1[index]);
5865 row0 = GetRowNumber(kink.GetR());
5866 kink.SetTPCRow0(row0);
5867 kink.SetLabel(CookLabel(seed0,0.5,0,row0),0);
5868 kink.SetLabel(CookLabel(seed1,0.5,row0,158),1);
5869 kink.SetIndex(-10,0);
5870 kink.SetIndex(int(param0[index].GetNumberOfClusters()+param1[index].GetNumberOfClusters()),1);
5871 kink.SetTPCncls(param0[index].GetNumberOfClusters(),0);
5872 kink.SetTPCncls(param1[index].GetNumberOfClusters(),1);
5875 // new (&mother) AliTPCseed(param0[index]);
5876 daughter = param1[index];
5877 daughter.SetLabel(kink.GetLabel(1));
5878 param0[index].Reset(kTRUE);
5879 FollowProlongation(param0[index],0);
5880 mother = param0[index];
5881 mother.SetLabel(kink.GetLabel(0));
5891 AliTPCseed* AliTPCtrackerMI::ReSeed(AliTPCseed *t)
5894 // reseed - refit - track
5897 // Int_t last = fSectors->GetNRows()-1;
5899 if (fSectors == fOuterSec){
5900 first = TMath::Max(first, t->GetFirstPoint()-fInnerSec->GetNRows());
5904 first = t->GetFirstPoint();
5906 AliTPCseed * seed = MakeSeed(t,0.1,0.5,0.9);
5907 FollowBackProlongation(*t,fSectors->GetNRows()-1);
5909 FollowProlongation(*t,first);
5919 //_____________________________________________________________________________
5920 Int_t AliTPCtrackerMI::ReadSeeds(const TFile *inp) {
5921 //-----------------------------------------------------------------
5922 // This function reades track seeds.
5923 //-----------------------------------------------------------------
5924 TDirectory *savedir=gDirectory;
5926 TFile *in=(TFile*)inp;
5927 if (!in->IsOpen()) {
5928 cerr<<"AliTPCtrackerMI::ReadSeeds(): input file is not open !\n";
5933 TTree *seedTree=(TTree*)in->Get("Seeds");
5935 cerr<<"AliTPCtrackerMI::ReadSeeds(): ";
5936 cerr<<"can't get a tree with track seeds !\n";
5939 AliTPCtrack *seed=new AliTPCtrack;
5940 seedTree->SetBranchAddress("tracks",&seed);
5942 if (fSeeds==0) fSeeds=new TObjArray(15000);
5944 Int_t n=(Int_t)seedTree->GetEntries();
5945 for (Int_t i=0; i<n; i++) {
5946 seedTree->GetEvent(i);
5947 fSeeds->AddLast(new AliTPCseed(*seed/*,seed->GetAlpha()*/));
5956 Int_t AliTPCtrackerMI::Clusters2Tracks (AliESDEvent *esd)
5959 if (fSeeds) DeleteSeeds();
5962 if (!fSeeds) return 1;
5969 //_____________________________________________________________________________
5970 Int_t AliTPCtrackerMI::Clusters2Tracks() {
5971 //-----------------------------------------------------------------
5972 // This is a track finder.
5973 //-----------------------------------------------------------------
5974 TDirectory *savedir=gDirectory;
5978 fSeeds = Tracking();
5981 Info("Clusters2Tracks","Time for tracking: \t");timer.Print();timer.Start();
5983 //activate again some tracks
5984 for (Int_t i=0; i<fSeeds->GetEntriesFast(); i++) {
5985 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
5987 Int_t nc=t.GetNumberOfClusters();
5989 delete fSeeds->RemoveAt(i);
5993 if (pt->GetRemoval()==10) {
5994 if (pt->GetDensityFirst(20)>0.8 || pt->GetDensityFirst(30)>0.8 || pt->GetDensityFirst(40)>0.7)
5995 pt->Desactivate(10); // make track again active
5997 pt->Desactivate(20);
5998 delete fSeeds->RemoveAt(i);
6003 RemoveUsed2(fSeeds,0.85,0.85,0);
6004 if (AliTPCReconstructor::GetRecoParam()->GetDoKinks()) FindKinks(fSeeds,fEvent);
6005 //FindCurling(fSeeds, fEvent,0);
6006 if (AliTPCReconstructor::StreamLevel()>2) FindMultiMC(fSeeds, fEvent,0); // find multi found tracks
6007 RemoveUsed2(fSeeds,0.5,0.4,20);
6008 FindSplitted(fSeeds, fEvent,0); // find multi found tracks
6011 // // refit short tracks
6013 Int_t nseed=fSeeds->GetEntriesFast();
6016 for (Int_t i=0; i<nseed; i++) {
6017 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
6019 Int_t nc=t.GetNumberOfClusters();
6021 delete fSeeds->RemoveAt(i);
6024 CookLabel(pt,0.1); //For comparison only
6025 //if ((pt->IsActive() || (pt->fRemoval==10) )&& nc>50 &&pt->GetNumberOfClusters()>0.4*pt->fNFoundable){
6026 if ((pt->IsActive() || (pt->GetRemoval()==10) )){
6028 if (fDebug>0) cerr<<found<<'\r';
6032 delete fSeeds->RemoveAt(i);
6036 //RemoveOverlap(fSeeds,0.99,7,kTRUE);
6038 //RemoveUsed(fSeeds,0.9,0.9,6);
6040 nseed=fSeeds->GetEntriesFast();
6042 for (Int_t i=0; i<nseed; i++) {
6043 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
6045 Int_t nc=t.GetNumberOfClusters();
6047 delete fSeeds->RemoveAt(i);
6051 t.CookdEdx(0.02,0.6);
6052 // CheckKinkPoint(&t,0.05);
6053 //if ((pt->IsActive() || (pt->fRemoval==10) )&& nc>50 &&pt->GetNumberOfClusters()>0.4*pt->fNFoundable){
6054 if ((pt->IsActive() || (pt->GetRemoval()==10) )){
6062 delete fSeeds->RemoveAt(i);
6063 //AliTPCseed * seed1 = ReSeed(pt,0.05,0.5,1);
6065 // FollowProlongation(*seed1,0);
6066 // Int_t n = seed1->GetNumberOfClusters();
6067 // printf("fP4\t%f\t%f\n",seed1->GetC(),pt->GetC());
6068 // printf("fN\t%d\t%d\n", seed1->GetNumberOfClusters(),pt->GetNumberOfClusters());
6071 //AliTPCseed * seed2 = ReSeed(pt,0.95,0.5,0.05);
6075 SortTracks(fSeeds, 1);
6079 PrepareForBackProlongation(fSeeds,5.);
6080 PropagateBack(fSeeds);
6081 printf("Time for back propagation: \t");timer.Print();timer.Start();
6085 PrepareForProlongation(fSeeds,5.);
6086 PropagateForward2(fSeeds);
6088 printf("Time for FORWARD propagation: \t");timer.Print();timer.Start();
6089 // RemoveUsed(fSeeds,0.7,0.7,6);
6090 //RemoveOverlap(fSeeds,0.9,7,kTRUE);
6092 nseed=fSeeds->GetEntriesFast();
6094 for (Int_t i=0; i<nseed; i++) {
6095 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
6097 Int_t nc=t.GetNumberOfClusters();
6099 delete fSeeds->RemoveAt(i);
6102 t.CookdEdx(0.02,0.6);
6103 // CookLabel(pt,0.1); //For comparison only
6104 //if ((pt->IsActive() || (pt->fRemoval==10) )&& nc>50 &&pt->GetNumberOfClusters()>0.4*pt->fNFoundable){
6105 if ((pt->IsActive() || (pt->fRemoval==10) )){
6106 cerr<<found++<<'\r';
6109 delete fSeeds->RemoveAt(i);
6114 // fNTracks = found;
6116 Info("Clusters2Tracks","Time for overlap removal, track writing and dedx cooking: \t"); timer.Print();timer.Start();
6119 // cerr<<"Number of found tracks : "<<"\t"<<found<<endl;
6120 Info("Clusters2Tracks","Number of found tracks %d",found);
6122 // UnloadClusters();
6127 void AliTPCtrackerMI::Tracking(TObjArray * arr)
6130 // tracking of the seeds
6133 fSectors = fOuterSec;
6134 ParallelTracking(arr,150,63);
6135 fSectors = fOuterSec;
6136 ParallelTracking(arr,63,0);
6139 TObjArray * AliTPCtrackerMI::Tracking(Int_t seedtype, Int_t i1, Int_t i2, Float_t cuts[4], Float_t dy, Int_t dsec)
6144 TObjArray * arr = new TObjArray;
6146 fSectors = fOuterSec;
6149 for (Int_t sec=0;sec<fkNOS;sec++){
6150 if (seedtype==3) MakeSeeds3(arr,sec,i1,i2,cuts,dy, dsec);
6151 if (seedtype==4) MakeSeeds5(arr,sec,i1,i2,cuts,dy);
6152 if (seedtype==2) MakeSeeds2(arr,sec,i1,i2,cuts,dy);
6155 Info("Tracking","\nSeeding - %d\t%d\t%d\t%d\n",seedtype,i1,i2,arr->GetEntriesFast());
6167 TObjArray * AliTPCtrackerMI::Tracking()
6171 if (AliTPCReconstructor::GetRecoParam()->GetSpecialSeeding()) return TrackingSpecial();
6174 Int_t nup=fOuterSec->GetNRows()+fInnerSec->GetNRows();
6176 TObjArray * seeds = new TObjArray;
6185 Float_t fnumber = 3.0;
6186 Float_t fdensity = 3.0;
6191 for (Int_t delta = 0; delta<18; delta+=6){
6195 arr = Tracking(3,nup-1-delta,nup-1-delta-gap,cuts,-1,1);
6196 SumTracks(seeds,arr);
6197 SignClusters(seeds,fnumber,fdensity);
6199 for (Int_t i=2;i<6;i+=2){
6200 // seed high pt tracks
6203 arr = Tracking(3,nup-i-delta,nup-i-delta-gap,cuts,-1,0);
6204 SumTracks(seeds,arr);
6205 SignClusters(seeds,fnumber,fdensity);
6210 // RemoveUsed(seeds,0.9,0.9,1);
6211 // UnsignClusters();
6212 // SignClusters(seeds,fnumber,fdensity);
6216 for (Int_t delta = 20; delta<120; delta+=10){
6218 // seed high pt tracks
6222 arr = Tracking(3,nup-delta,nup-delta-gap,cuts,-1);
6223 SumTracks(seeds,arr);
6224 SignClusters(seeds,fnumber,fdensity);
6229 arr = Tracking(3,nup-delta-5,nup-delta-5-gap,cuts,-1);
6230 SumTracks(seeds,arr);
6231 SignClusters(seeds,fnumber,fdensity);
6242 Info("Tracking()","\n\nPrimary seeding\t%d\n\n",seeds->GetEntriesFast());
6246 // RemoveUsed(seeds,0.75,0.75,1);
6248 //SignClusters(seeds,fnumber,fdensity);
6257 arr = Tracking(4,nup-1,nup-1-gap,cuts,-1);
6258 SumTracks(seeds,arr);
6259 SignClusters(seeds,fnumber,fdensity);
6261 arr = Tracking(4,nup-2,nup-2-gap,cuts,-1);
6262 SumTracks(seeds,arr);
6263 SignClusters(seeds,fnumber,fdensity);
6265 arr = Tracking(4,nup-3,nup-3-gap,cuts,-1);
6266 SumTracks(seeds,arr);
6267 SignClusters(seeds,fnumber,fdensity);
6271 for (Int_t delta = 3; delta<30; delta+=5){
6277 arr = Tracking(4,nup-1-delta,nup-1-delta-gap,cuts,-1);
6278 SumTracks(seeds,arr);
6279 SignClusters(seeds,fnumber,fdensity);
6281 arr = Tracking(4,nup-3-delta,nup-5-delta-gap,cuts,4);
6282 SumTracks(seeds,arr);
6283 SignClusters(seeds,fnumber,fdensity);
6294 Int_t fLastSeedRowSec=AliTPCReconstructor::GetRecoParam()->GetLastSeedRowSec();
6297 for (Int_t delta = 30; delta<fLastSeedRowSec; delta+=10){
6303 arr = Tracking(4,nup-1-delta,nup-1-delta-gap,cuts,-1);
6304 SumTracks(seeds,arr);
6305 SignClusters(seeds,fnumber,fdensity);
6307 arr = Tracking(4,nup-5-delta,nup-5-delta-gap,cuts,5 );
6308 SumTracks(seeds,arr);
6309 SignClusters(seeds,fnumber,fdensity);
6313 Info("Tracking()","\n\nSecondary seeding\t%d\n\n",seeds->GetEntriesFast());
6324 TObjArray * AliTPCtrackerMI::TrackingSpecial()
6327 // seeding adjusted for laser and cosmic tests - short tracks with big inclination angle
6328 // no primary vertex seeding tried
6332 Int_t nup=fOuterSec->GetNRows()+fInnerSec->GetNRows();
6334 TObjArray * seeds = new TObjArray;
6339 Float_t fnumber = 3.0;
6340 Float_t fdensity = 3.0;
6343 cuts[0] = AliTPCReconstructor::GetRecoParam()->GetMaxC(); // max curvature
6344 cuts[1] = 3.5; // max tan(phi) angle for seeding
6345 cuts[2] = 3.; // not used (cut on z primary vertex)
6346 cuts[3] = 3.5; // max tan(theta) angle for seeding
6348 for (Int_t delta = 0; nup-delta-gap-1>0; delta+=3){
6350 arr = Tracking(4,nup-1-delta,nup-1-delta-gap,cuts,-1);
6351 SumTracks(seeds,arr);
6352 SignClusters(seeds,fnumber,fdensity);
6356 Info("Tracking()","\n\nSecondary seeding\t%d\n\n",seeds->GetEntriesFast());
6367 void AliTPCtrackerMI::SumTracks(TObjArray *arr1,TObjArray *&arr2) const
6370 //sum tracks to common container
6371 //remove suspicious tracks
6372 Int_t nseed = arr2->GetEntriesFast();
6373 for (Int_t i=0;i<nseed;i++){
6374 AliTPCseed *pt=(AliTPCseed*)arr2->UncheckedAt(i);
6377 // remove tracks with too big curvature
6379 if (TMath::Abs(pt->GetC())>AliTPCReconstructor::GetRecoParam()->GetMaxC()){
6380 delete arr2->RemoveAt(i);
6383 // REMOVE VERY SHORT TRACKS
6384 if (pt->GetNumberOfClusters()<20){
6385 delete arr2->RemoveAt(i);
6388 // NORMAL ACTIVE TRACK
6389 if (pt->IsActive()){
6390 arr1->AddLast(arr2->RemoveAt(i));
6393 //remove not usable tracks
6394 if (pt->GetRemoval()!=10){
6395 delete arr2->RemoveAt(i);
6399 // ENABLE ONLY ENOUGH GOOD STOPPED TRACKS
6400 if (pt->GetDensityFirst(20)>0.8 || pt->GetDensityFirst(30)>0.8 || pt->GetDensityFirst(40)>0.7)
6401 arr1->AddLast(arr2->RemoveAt(i));
6403 delete arr2->RemoveAt(i);
6407 delete arr2; arr2 = 0;
6412 void AliTPCtrackerMI::ParallelTracking(TObjArray * arr, Int_t rfirst, Int_t rlast)
6415 // try to track in parralel
6417 Int_t nseed=arr->GetEntriesFast();
6418 //prepare seeds for tracking
6419 for (Int_t i=0; i<nseed; i++) {
6420 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i), &t=*pt;
6422 if (!t.IsActive()) continue;
6423 // follow prolongation to the first layer
6424 if ( (fSectors ==fInnerSec) || (t.GetFirstPoint()-fParam->GetNRowLow()>rfirst+1) )
6425 FollowProlongation(t, rfirst+1);
6430 for (Int_t nr=rfirst; nr>=rlast; nr--){
6431 if (nr<fInnerSec->GetNRows())
6432 fSectors = fInnerSec;
6434 fSectors = fOuterSec;
6435 // make indexes with the cluster tracks for given
6437 // find nearest cluster
6438 for (Int_t i=0; i<nseed; i++) {
6439 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i), &t=*pt;
6441 if (nr==80) pt->UpdateReference();
6442 if (!pt->IsActive()) continue;
6443 // if ( (fSectors ==fOuterSec) && (pt->fFirstPoint-fParam->GetNRowLow())<nr) continue;
6444 if (pt->GetRelativeSector()>17) {
6447 UpdateClusters(t,nr);
6449 // prolonagate to the nearest cluster - if founded
6450 for (Int_t i=0; i<nseed; i++) {
6451 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
6453 if (!pt->IsActive()) continue;
6454 // if ((fSectors ==fOuterSec) && (pt->fFirstPoint-fParam->GetNRowLow())<nr) continue;
6455 if (pt->GetRelativeSector()>17) {
6458 FollowToNextCluster(*pt,nr);
6463 void AliTPCtrackerMI::PrepareForBackProlongation(TObjArray * arr,Float_t fac) const
6467 // if we use TPC track itself we have to "update" covariance
6469 Int_t nseed= arr->GetEntriesFast();
6470 for (Int_t i=0;i<nseed;i++){
6471 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
6475 //rotate to current local system at first accepted point
6476 Int_t index = pt->GetClusterIndex2(pt->GetFirstPoint());
6477 Int_t sec = (index&0xff000000)>>24;
6479 Float_t angle1 = fInnerSec->GetAlpha()*sec+fInnerSec->GetAlphaShift();
6480 if (angle1>TMath::Pi())
6481 angle1-=2.*TMath::Pi();
6482 Float_t angle2 = pt->GetAlpha();
6484 if (TMath::Abs(angle1-angle2)>0.001){
6485 pt->Rotate(angle1-angle2);
6486 //angle2 = pt->GetAlpha();
6487 //pt->fRelativeSector = pt->GetAlpha()/fInnerSec->GetAlpha();
6488 //if (pt->GetAlpha()<0)
6489 // pt->fRelativeSector+=18;
6490 //sec = pt->fRelativeSector;
6499 void AliTPCtrackerMI::PrepareForProlongation(TObjArray * arr, Float_t fac) const
6503 // if we use TPC track itself we have to "update" covariance
6505 Int_t nseed= arr->GetEntriesFast();
6506 for (Int_t i=0;i<nseed;i++){
6507 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
6510 pt->SetFirstPoint(pt->GetLastPoint());
6518 Int_t AliTPCtrackerMI::PropagateBack(TObjArray * arr)
6521 // make back propagation
6523 Int_t nseed= arr->GetEntriesFast();
6524 for (Int_t i=0;i<nseed;i++){
6525 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
6526 if (pt&& pt->GetKinkIndex(0)<=0) {
6527 //AliTPCseed *pt2 = new AliTPCseed(*pt);
6528 fSectors = fInnerSec;
6529 //FollowBackProlongation(*pt,fInnerSec->GetNRows()-1);
6530 //fSectors = fOuterSec;
6531 FollowBackProlongation(*pt,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1);
6532 //if (pt->GetNumberOfClusters()<(pt->fEsd->GetTPCclusters(0)) ){
6533 // Error("PropagateBack","Not prolonged track %d",pt->GetLabel());
6534 // FollowBackProlongation(*pt2,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1);
6537 if (pt&& pt->GetKinkIndex(0)>0) {
6538 AliESDkink * kink = fEvent->GetKink(pt->GetKinkIndex(0)-1);
6539 pt->SetFirstPoint(kink->GetTPCRow0());
6540 fSectors = fInnerSec;
6541 FollowBackProlongation(*pt,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1);
6549 Int_t AliTPCtrackerMI::PropagateForward2(TObjArray * arr)
6552 // make forward propagation
6554 Int_t nseed= arr->GetEntriesFast();
6556 for (Int_t i=0;i<nseed;i++){
6557 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
6559 FollowProlongation(*pt,0);
6568 Int_t AliTPCtrackerMI::PropagateForward()
6571 // propagate track forward
6573 Int_t nseed = fSeeds->GetEntriesFast();
6574 for (Int_t i=0;i<nseed;i++){
6575 AliTPCseed *pt = (AliTPCseed*)fSeeds->UncheckedAt(i);
6577 AliTPCseed &t = *pt;
6578 Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
6579 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
6580 if (alpha < 0. ) alpha += 2.*TMath::Pi();
6581 t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN);
6585 fSectors = fOuterSec;
6586 ParallelTracking(fSeeds,fOuterSec->GetNRows()+fInnerSec->GetNRows()-1,fInnerSec->GetNRows());
6587 fSectors = fInnerSec;
6588 ParallelTracking(fSeeds,fInnerSec->GetNRows()-1,0);
6597 Int_t AliTPCtrackerMI::PropagateBack(AliTPCseed * pt, Int_t row0, Int_t row1)
6600 // make back propagation, in between row0 and row1
6604 fSectors = fInnerSec;
6607 if (row1<fSectors->GetNRows())
6610 r1 = fSectors->GetNRows()-1;
6612 if (row0<fSectors->GetNRows()&& r1>0 )
6613 FollowBackProlongation(*pt,r1);
6614 if (row1<=fSectors->GetNRows())
6617 r1 = row1 - fSectors->GetNRows();
6618 if (r1<=0) return 0;
6619 if (r1>=fOuterSec->GetNRows()) return 0;
6620 fSectors = fOuterSec;
6621 return FollowBackProlongation(*pt,r1);
6629 void AliTPCtrackerMI::GetShape(AliTPCseed * seed, Int_t row)
6633 AliTPCClusterParam * clparam = AliTPCcalibDB::Instance()->GetClusterParam();
6634 Float_t zdrift = TMath::Abs((fParam->GetZLength(0)-TMath::Abs(seed->GetZ())));
6635 Int_t type = (seed->GetSector() < fParam->GetNSector()/2) ? 0: (row>126) ? 1:2;
6636 Double_t angulary = seed->GetSnp();
6637 angulary = angulary*angulary/(1.-angulary*angulary);
6638 Double_t angularz = seed->GetTgl()*seed->GetTgl()*(1.+angulary);
6640 Double_t sigmay = clparam->GetRMS0(0,type,zdrift,TMath::Sqrt(TMath::Abs(angulary)));
6641 Double_t sigmaz = clparam->GetRMS0(1,type,zdrift,TMath::Sqrt(TMath::Abs(angularz)));
6642 seed->SetCurrentSigmaY2(sigmay*sigmay);
6643 seed->SetCurrentSigmaZ2(sigmaz*sigmaz);
6644 // Float_t sd2 = TMath::Abs((fParam->GetZLength(0)-TMath::Abs(seed->GetZ())))*fParam->GetDiffL()*fParam->GetDiffL();
6645 // // Float_t padlength = fParam->GetPadPitchLength(seed->fSector);
6646 // Float_t padlength = GetPadPitchLength(row);
6648 // Float_t sresy = (seed->GetSector() < fParam->GetNSector()/2) ? 0.2 :0.3;
6649 // seed->SetCurrentSigmaY2(sd2+padlength*padlength*angulary/12.+sresy*sresy);
6651 // Float_t sresz = fParam->GetZSigma();
6652 // seed->SetCurrentSigmaZ2(sd2+padlength*padlength*angularz*angularz*(1+angulary)/12.+sresz*sresz);
6654 Float_t wy = GetSigmaY(seed);
6655 Float_t wz = GetSigmaZ(seed);
6658 if (TMath::Abs(wy/seed->fCurrentSigmaY2-1)>0.0001 || TMath::Abs(wz/seed->fCurrentSigmaZ2-1)>0.0001 ){
6659 printf("problem\n");
6666 //__________________________________________________________________________
6667 void AliTPCtrackerMI::CookLabel(AliKalmanTrack *tk, Float_t wrong) const {
6668 //--------------------------------------------------------------------
6669 //This function "cooks" a track label. If label<0, this track is fake.
6670 //--------------------------------------------------------------------
6671 AliTPCseed * t = dynamic_cast<AliTPCseed*>(tk);
6673 printf("%s:%d wrong type \n",(char*)__FILE__,__LINE__);
6677 Int_t noc=t->GetNumberOfClusters();
6679 //printf("\nnot founded prolongation\n\n\n");
6685 AliTPCclusterMI *clusters[160];
6687 for (Int_t i=0;i<160;i++) {
6694 for (i=0; i<160 && current<noc; i++) {
6696 Int_t index=t->GetClusterIndex2(i);
6697 if (index<=0) continue;
6698 if (index&0x8000) continue;
6700 //clusters[current]=GetClusterMI(index);
6701 if (t->GetClusterPointer(i)){
6702 clusters[current]=t->GetClusterPointer(i);
6708 Int_t lab=123456789;
6709 for (i=0; i<noc; i++) {
6710 AliTPCclusterMI *c=clusters[i];
6712 lab=TMath::Abs(c->GetLabel(0));
6714 for (j=0; j<noc; j++) if (lb[j]==lab || mx[j]==0) break;
6720 for (i=0; i<noc; i++) if (mx[i]>max) {max=mx[i]; lab=lb[i];}
6722 for (i=0; i<noc; i++) {
6723 AliTPCclusterMI *c=clusters[i];
6725 if (TMath::Abs(c->GetLabel(1)) == lab ||
6726 TMath::Abs(c->GetLabel(2)) == lab ) max++;
6729 if ((1.- Float_t(max)/noc) > wrong) lab=-lab;
6732 Int_t tail=Int_t(0.10*noc);
6735 for (i=1; i<=160&&ind<tail; i++) {
6736 // AliTPCclusterMI *c=clusters[noc-i];
6737 AliTPCclusterMI *c=clusters[i];
6739 if (lab == TMath::Abs(c->GetLabel(0)) ||
6740 lab == TMath::Abs(c->GetLabel(1)) ||
6741 lab == TMath::Abs(c->GetLabel(2))) max++;
6744 if (max < Int_t(0.5*tail)) lab=-lab;
6751 //delete[] clusters;
6755 //__________________________________________________________________________
6756 Int_t AliTPCtrackerMI::CookLabel(AliTPCseed *t, Float_t wrong,Int_t first, Int_t last) const {
6757 //--------------------------------------------------------------------
6758 //This function "cooks" a track label. If label<0, this track is fake.
6759 //--------------------------------------------------------------------
6760 Int_t noc=t->GetNumberOfClusters();
6762 //printf("\nnot founded prolongation\n\n\n");
6768 AliTPCclusterMI *clusters[160];
6770 for (Int_t i=0;i<160;i++) {
6777 for (i=0; i<160 && current<noc; i++) {
6778 if (i<first) continue;
6779 if (i>last) continue;
6780 Int_t index=t->GetClusterIndex2(i);
6781 if (index<=0) continue;
6782 if (index&0x8000) continue;
6784 //clusters[current]=GetClusterMI(index);
6785 if (t->GetClusterPointer(i)){
6786 clusters[current]=t->GetClusterPointer(i);
6791 if (noc<5) return -1;
6792 Int_t lab=123456789;
6793 for (i=0; i<noc; i++) {
6794 AliTPCclusterMI *c=clusters[i];
6796 lab=TMath::Abs(c->GetLabel(0));
6798 for (j=0; j<noc; j++) if (lb[j]==lab || mx[j]==0) break;
6804 for (i=0; i<noc; i++) if (mx[i]>max) {max=mx[i]; lab=lb[i];}
6806 for (i=0; i<noc; i++) {
6807 AliTPCclusterMI *c=clusters[i];
6809 if (TMath::Abs(c->GetLabel(1)) == lab ||
6810 TMath::Abs(c->GetLabel(2)) == lab ) max++;
6813 if ((1.- Float_t(max)/noc) > wrong) lab=-lab;
6816 Int_t tail=Int_t(0.10*noc);
6819 for (i=1; i<=160&&ind<tail; i++) {
6820 // AliTPCclusterMI *c=clusters[noc-i];
6821 AliTPCclusterMI *c=clusters[i];
6823 if (lab == TMath::Abs(c->GetLabel(0)) ||
6824 lab == TMath::Abs(c->GetLabel(1)) ||
6825 lab == TMath::Abs(c->GetLabel(2))) max++;
6828 if (max < Int_t(0.5*tail)) lab=-lab;
6831 // t->SetLabel(lab);
6835 //delete[] clusters;
6839 Int_t AliTPCtrackerMI::GetRowNumber(Double_t x[3]) const
6841 //return pad row number for given x vector
6842 Float_t phi = TMath::ATan2(x[1],x[0]);
6843 if(phi<0) phi=2.*TMath::Pi()+phi;
6844 // Get the local angle in the sector philoc
6845 const Float_t kRaddeg = 180/3.14159265358979312;
6846 Float_t phiangle = (Int_t (phi*kRaddeg/20.) + 0.5)*20./kRaddeg;
6847 Double_t localx = x[0]*TMath::Cos(phiangle)-x[1]*TMath::Sin(phiangle);
6848 return GetRowNumber(localx);
6853 void AliTPCtrackerMI::MakeBitmaps(AliTPCseed *t)
6855 //-----------------------------------------------------------------------
6856 // Fill the cluster and sharing bitmaps of the track
6857 //-----------------------------------------------------------------------
6859 Int_t firstpoint = 0;
6860 Int_t lastpoint = 159;
6861 AliTPCTrackerPoint *point;
6863 for (int iter=firstpoint; iter<lastpoint; iter++) {
6864 point = t->GetTrackPoint(iter);
6866 t->SetClusterMapBit(iter, kTRUE);
6867 if (point->IsShared())
6868 t->SetSharedMapBit(iter,kTRUE);
6870 t->SetSharedMapBit(iter, kFALSE);
6873 t->SetClusterMapBit(iter, kFALSE);
6874 t->SetSharedMapBit(iter, kFALSE);