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);
72 // The debug level - different procedure produce tree for numerical debugging
73 // To enable them set AliTPCReconstructor::SetStreamLevel(n); where nis bigger 1
77 // Adding systematic errors to the covariance:
79 // The systematic errors due to the misalignment and miscalibration are added to the covariance matrix
80 // of the tracks (not to the clusters as they are dependent):
81 // The parameters form AliTPCRecoParam are used AliTPCRecoParam::GetSystematicError
82 // The systematic errors are expressed there in RMS - position (cm), angle (rad), curvature (1/cm)
83 // The default values are 0.
85 // The sytematic errors are added to the covariance matrix in following places:
87 // 1. During fisrt itteration - AliTPCtrackerMI::FillESD
88 // 2. Second iteration -
89 // 2.a ITS->TPC - AliTPCtrackerMI::ReadSeeds
90 // 2.b TPC->TRD - AliTPCtrackerMI::PropagateBack
91 // 3. Third iteration -
92 // 3.a TRD->TPC - AliTPCtrackerMI::ReadSeeds
93 // 3.b TPC->ITS - AliTPCtrackerMI::RefitInward
95 // There are several places in the code which can be numerically debuged
96 // This code is keeped in order to enable code development and to check the calibration implementtion
98 // 1. ErrParam stream - dump information about
100 // 2.a) cluster error estimate
101 // 3.a) cluster shape estimate
104 // Debug streamer levels:
106 //-------------------------------------------------------
111 #include "Riostream.h"
112 #include <TClonesArray.h>
114 #include <TObjArray.h>
117 #include "AliComplexCluster.h"
118 #include "AliESDEvent.h"
119 #include "AliESDtrack.h"
120 #include "AliESDVertex.h"
123 #include "AliHelix.h"
124 #include "AliRunLoader.h"
125 #include "AliTPCClustersRow.h"
126 #include "AliTPCParam.h"
127 #include "AliTPCReconstructor.h"
128 #include "AliTPCpolyTrack.h"
129 #include "AliTPCreco.h"
130 #include "AliTPCseed.h"
132 #include "AliTPCtrackerSector.h"
133 #include "AliTPCtrackerMI.h"
134 #include "TStopwatch.h"
135 #include "AliTPCReconstructor.h"
136 #include "AliAlignObj.h"
137 #include "AliTrackPointArray.h"
139 #include "AliTPCcalibDB.h"
140 #include "AliTPCTransform.h"
141 #include "AliTPCClusterParam.h"
145 ClassImp(AliTPCtrackerMI)
149 class AliTPCFastMath {
152 static Double_t FastAsin(Double_t x);
154 static Double_t fgFastAsin[20000]; //lookup table for fast asin computation
157 Double_t AliTPCFastMath::fgFastAsin[20000];
158 AliTPCFastMath gAliTPCFastMath; // needed to fill the LUT
160 AliTPCFastMath::AliTPCFastMath(){
162 // initialized lookup table;
163 for (Int_t i=0;i<10000;i++){
164 fgFastAsin[2*i] = TMath::ASin(i/10000.);
165 fgFastAsin[2*i+1] = (TMath::ASin((i+1)/10000.)-fgFastAsin[2*i]);
169 Double_t AliTPCFastMath::FastAsin(Double_t x){
171 // return asin using lookup table
173 Int_t index = int(x*10000);
174 return fgFastAsin[2*index]+(x*10000.-index)*fgFastAsin[2*index+1];
177 Int_t index = int(x*10000);
178 return -(fgFastAsin[2*index]+(x*10000.-index)*fgFastAsin[2*index+1]);
180 //__________________________________________________________________
181 AliTPCtrackerMI::AliTPCtrackerMI()
203 // default constructor
206 //_____________________________________________________________________
210 Int_t AliTPCtrackerMI::UpdateTrack(AliTPCseed * track, Int_t accept){
212 //update track information using current cluster - track->fCurrentCluster
215 AliTPCclusterMI* c =track->GetCurrentCluster();
216 if (accept>0) track->SetCurrentClusterIndex1(track->GetCurrentClusterIndex1() | 0x8000); //sign not accepted clusters
218 UInt_t i = track->GetCurrentClusterIndex1();
220 Int_t sec=(i&0xff000000)>>24;
221 //Int_t row = (i&0x00ff0000)>>16;
222 track->SetRow((i&0x00ff0000)>>16);
223 track->SetSector(sec);
224 // Int_t index = i&0xFFFF;
225 if (sec>=fkParam->GetNInnerSector()) track->SetRow(track->GetRow()+fkParam->GetNRowLow());
226 track->SetClusterIndex2(track->GetRow(), i);
227 //track->fFirstPoint = row;
228 //if ( track->fLastPoint<row) track->fLastPoint =row;
229 // if (track->fRow<0 || track->fRow>160) {
230 // printf("problem\n");
232 if (track->GetFirstPoint()>track->GetRow())
233 track->SetFirstPoint(track->GetRow());
234 if (track->GetLastPoint()<track->GetRow())
235 track->SetLastPoint(track->GetRow());
238 track->SetClusterPointer(track->GetRow(),c);
241 Double_t angle2 = track->GetSnp()*track->GetSnp();
243 //SET NEW Track Point
245 if (angle2<1) //PH sometimes angle2 is very big. To be investigated...
247 angle2 = TMath::Sqrt(angle2/(1-angle2));
248 AliTPCTrackerPoint &point =*(track->GetTrackPoint(track->GetRow()));
250 point.SetSigmaY(c->GetSigmaY2()/track->GetCurrentSigmaY2());
251 point.SetSigmaZ(c->GetSigmaZ2()/track->GetCurrentSigmaZ2());
252 point.SetErrY(sqrt(track->GetErrorY2()));
253 point.SetErrZ(sqrt(track->GetErrorZ2()));
255 point.SetX(track->GetX());
256 point.SetY(track->GetY());
257 point.SetZ(track->GetZ());
258 point.SetAngleY(angle2);
259 point.SetAngleZ(track->GetTgl());
260 if (point.IsShared()){
261 track->SetErrorY2(track->GetErrorY2()*4);
262 track->SetErrorZ2(track->GetErrorZ2()*4);
266 Double_t chi2 = track->GetPredictedChi2(track->GetCurrentCluster());
268 // track->SetErrorY2(track->GetErrorY2()*1.3);
269 // track->SetErrorY2(track->GetErrorY2()+0.01);
270 // track->SetErrorZ2(track->GetErrorZ2()*1.3);
271 // track->SetErrorZ2(track->GetErrorZ2()+0.005);
273 if (accept>0) return 0;
274 if (track->GetNumberOfClusters()%20==0){
275 // if (track->fHelixIn){
276 // TClonesArray & larr = *(track->fHelixIn);
277 // Int_t ihelix = larr.GetEntriesFast();
278 // new(larr[ihelix]) AliHelix(*track) ;
281 track->SetNoCluster(0);
282 return track->Update(c,chi2,i);
287 Int_t AliTPCtrackerMI::AcceptCluster(AliTPCseed * seed, AliTPCclusterMI * cluster)
290 // decide according desired precision to accept given
291 // cluster for tracking
293 seed->GetProlongation(cluster->GetX(),yt,zt);
294 Double_t sy2=ErrY2(seed,cluster);
295 Double_t sz2=ErrZ2(seed,cluster);
297 Double_t sdistancey2 = sy2+seed->GetSigmaY2();
298 Double_t sdistancez2 = sz2+seed->GetSigmaZ2();
299 Double_t dy=seed->GetCurrentCluster()->GetY()-yt;
300 Double_t dz=seed->GetCurrentCluster()->GetY()-zt;
301 Double_t rdistancey2 = (seed->GetCurrentCluster()->GetY()-yt)*
302 (seed->GetCurrentCluster()->GetY()-yt)/sdistancey2;
303 Double_t rdistancez2 = (seed->GetCurrentCluster()->GetZ()-zt)*
304 (seed->GetCurrentCluster()->GetZ()-zt)/sdistancez2;
306 Double_t rdistance2 = rdistancey2+rdistancez2;
309 if (AliTPCReconstructor::StreamLevel()>2 && seed->GetNumberOfClusters()>20) {
310 Float_t rmsy2 = seed->GetCurrentSigmaY2();
311 Float_t rmsz2 = seed->GetCurrentSigmaZ2();
312 Float_t rmsy2p30 = seed->GetCMeanSigmaY2p30();
313 Float_t rmsz2p30 = seed->GetCMeanSigmaZ2p30();
314 Float_t rmsy2p30R = seed->GetCMeanSigmaY2p30R();
315 Float_t rmsz2p30R = seed->GetCMeanSigmaZ2p30R();
316 AliExternalTrackParam param(*seed);
317 static TVectorD gcl(3),gtr(3);
319 param.GetXYZ(gcl.GetMatrixArray());
320 cluster->GetGlobalXYZ(gclf);
321 gcl[0]=gclf[0]; gcl[1]=gclf[1]; gcl[2]=gclf[2];
324 if (AliTPCReconstructor::StreamLevel()>2) {
325 (*fDebugStreamer)<<"ErrParam"<<
338 "rmsy2p30="<<rmsy2p30<<
339 "rmsz2p30="<<rmsz2p30<<
340 "rmsy2p30R="<<rmsy2p30R<<
341 "rmsz2p30R="<<rmsz2p30R<<
345 //return 0; // temporary
346 if (rdistance2>32) return 3;
349 if ((rdistancey2>9. || rdistancez2>9.) && cluster->GetType()==0)
350 return 2; //suspisiouce - will be changed
352 if ((rdistancey2>6.25 || rdistancez2>6.25) && cluster->GetType()>0)
353 // strict cut on overlaped cluster
354 return 2; //suspisiouce - will be changed
356 if ( (rdistancey2>1. || rdistancez2>6.25 )
357 && cluster->GetType()<0){
358 seed->SetNFoundable(seed->GetNFoundable()-1);
368 //_____________________________________________________________________________
369 AliTPCtrackerMI::AliTPCtrackerMI(const AliTPCParam *par):
371 fkNIS(par->GetNInnerSector()/2),
373 fkNOS(par->GetNOuterSector()/2),
390 //---------------------------------------------------------------------
391 // The main TPC tracker constructor
392 //---------------------------------------------------------------------
393 fInnerSec=new AliTPCtrackerSector[fkNIS];
394 fOuterSec=new AliTPCtrackerSector[fkNOS];
397 for (i=0; i<fkNIS; i++) fInnerSec[i].Setup(par,0);
398 for (i=0; i<fkNOS; i++) fOuterSec[i].Setup(par,1);
401 Int_t nrowlow = par->GetNRowLow();
402 Int_t nrowup = par->GetNRowUp();
405 for (i=0;i<nrowlow;i++){
406 fXRow[i] = par->GetPadRowRadiiLow(i);
407 fPadLength[i]= par->GetPadPitchLength(0,i);
408 fYMax[i] = fXRow[i]*TMath::Tan(0.5*par->GetInnerAngle());
412 for (i=0;i<nrowup;i++){
413 fXRow[i+nrowlow] = par->GetPadRowRadiiUp(i);
414 fPadLength[i+nrowlow] = par->GetPadPitchLength(60,i);
415 fYMax[i+nrowlow] = fXRow[i+nrowlow]*TMath::Tan(0.5*par->GetOuterAngle());
418 if (AliTPCReconstructor::StreamLevel()>0) {
419 fDebugStreamer = new TTreeSRedirector("TPCdebug.root");
422 //________________________________________________________________________
423 AliTPCtrackerMI::AliTPCtrackerMI(const AliTPCtrackerMI &t):
444 //------------------------------------
445 // dummy copy constructor
446 //------------------------------------------------------------------
449 AliTPCtrackerMI & AliTPCtrackerMI::operator=(const AliTPCtrackerMI& /*r*/)
451 //------------------------------
453 //--------------------------------------------------------------
456 //_____________________________________________________________________________
457 AliTPCtrackerMI::~AliTPCtrackerMI() {
458 //------------------------------------------------------------------
459 // TPC tracker destructor
460 //------------------------------------------------------------------
467 if (fDebugStreamer) delete fDebugStreamer;
471 void AliTPCtrackerMI::FillESD(const TObjArray* arr)
475 //fill esds using updated tracks
477 // write tracks to the event
478 // store index of the track
479 Int_t nseed=arr->GetEntriesFast();
480 //FindKinks(arr,fEvent);
481 for (Int_t i=0; i<nseed; i++) {
482 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
486 if (AliTPCReconstructor::StreamLevel()>1) {
487 (*fDebugStreamer)<<"Track0"<<
491 // pt->PropagateTo(fkParam->GetInnerRadiusLow());
492 if (pt->GetKinkIndex(0)<=0){ //don't propagate daughter tracks
493 pt->PropagateTo(fkParam->GetInnerRadiusLow());
496 if (( pt->GetPoints()[2]- pt->GetPoints()[0])>5 && pt->GetPoints()[3]>0.8){
498 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
499 iotrack.SetTPCPoints(pt->GetPoints());
500 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
501 iotrack.SetV0Indexes(pt->GetV0Indexes());
502 // iotrack.SetTPCpid(pt->fTPCr);
503 //iotrack.SetTPCindex(i);
504 fEvent->AddTrack(&iotrack);
508 if ( (pt->GetNumberOfClusters()>70)&& (Float_t(pt->GetNumberOfClusters())/Float_t(pt->GetNFoundable()))>0.55) {
510 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
511 iotrack.SetTPCPoints(pt->GetPoints());
512 //iotrack.SetTPCindex(i);
513 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
514 iotrack.SetV0Indexes(pt->GetV0Indexes());
515 // iotrack.SetTPCpid(pt->fTPCr);
516 fEvent->AddTrack(&iotrack);
520 // short tracks - maybe decays
522 if ( (pt->GetNumberOfClusters()>30) && (Float_t(pt->GetNumberOfClusters())/Float_t(pt->GetNFoundable()))>0.70) {
523 Int_t found,foundable,shared;
524 pt->GetClusterStatistic(0,60,found, foundable,shared,kFALSE);
525 if ( (found>20) && (pt->GetNShared()/float(pt->GetNumberOfClusters())<0.2)){
527 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
528 //iotrack.SetTPCindex(i);
529 iotrack.SetTPCPoints(pt->GetPoints());
530 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
531 iotrack.SetV0Indexes(pt->GetV0Indexes());
532 //iotrack.SetTPCpid(pt->fTPCr);
533 fEvent->AddTrack(&iotrack);
538 if ( (pt->GetNumberOfClusters()>20) && (Float_t(pt->GetNumberOfClusters())/Float_t(pt->GetNFoundable()))>0.8) {
539 Int_t found,foundable,shared;
540 pt->GetClusterStatistic(0,60,found, foundable,shared,kFALSE);
541 if (found<20) continue;
542 if (pt->GetNShared()/float(pt->GetNumberOfClusters())>0.2) continue;
545 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
546 iotrack.SetTPCPoints(pt->GetPoints());
547 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
548 iotrack.SetV0Indexes(pt->GetV0Indexes());
549 //iotrack.SetTPCpid(pt->fTPCr);
550 //iotrack.SetTPCindex(i);
551 fEvent->AddTrack(&iotrack);
554 // short tracks - secondaties
556 if ( (pt->GetNumberOfClusters()>30) ) {
557 Int_t found,foundable,shared;
558 pt->GetClusterStatistic(128,158,found, foundable,shared,kFALSE);
559 if ( (found>20) && (pt->GetNShared()/float(pt->GetNumberOfClusters())<0.2) &&float(found)/float(foundable)>0.8){
561 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
562 iotrack.SetTPCPoints(pt->GetPoints());
563 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
564 iotrack.SetV0Indexes(pt->GetV0Indexes());
565 //iotrack.SetTPCpid(pt->fTPCr);
566 //iotrack.SetTPCindex(i);
567 fEvent->AddTrack(&iotrack);
572 if ( (pt->GetNumberOfClusters()>15)) {
573 Int_t found,foundable,shared;
574 pt->GetClusterStatistic(138,158,found, foundable,shared,kFALSE);
575 if (found<15) continue;
576 if (foundable<=0) continue;
577 if (pt->GetNShared()/float(pt->GetNumberOfClusters())>0.2) continue;
578 if (float(found)/float(foundable)<0.8) continue;
581 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
582 iotrack.SetTPCPoints(pt->GetPoints());
583 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
584 iotrack.SetV0Indexes(pt->GetV0Indexes());
585 // iotrack.SetTPCpid(pt->fTPCr);
586 //iotrack.SetTPCindex(i);
587 fEvent->AddTrack(&iotrack);
591 // >> account for suppressed tracks in the kink indices (RS)
592 int nESDtracks = fEvent->GetNumberOfTracks();
593 for (int it=nESDtracks;it--;) {
594 AliESDtrack* esdTr = fEvent->GetTrack(it);
595 if (!esdTr || !esdTr->GetKinkIndex(0)) continue;
596 for (int ik=0;ik<3;ik++) {
598 if (!(knkId=esdTr->GetKinkIndex(ik))) break; // no more kinks for this track
599 AliESDkink* kink = fEvent->GetKink(TMath::Abs(knkId)-1);
601 AliError(Form("ESDTrack%d refers to non-existing kink %d",it,TMath::Abs(knkId)-1));
604 kink->SetIndex(it, knkId<0 ? 0:1); // update track index of the kink: mother at 0, daughter at 1
607 // << account for suppressed tracks in the kink indices (RS)
609 printf("Number of filled ESDs-\t%d\n",fEvent->GetNumberOfTracks());
616 Double_t AliTPCtrackerMI::ErrY2(AliTPCseed* seed, const AliTPCclusterMI * cl){
619 // Use calibrated cluster error from OCDB
621 AliTPCClusterParam * clparam = AliTPCcalibDB::Instance()->GetClusterParam();
623 Float_t z = TMath::Abs(fkParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
624 Int_t ctype = cl->GetType();
625 Int_t type = (cl->GetRow()<63) ? 0: (cl->GetRow()>126) ? 1:2;
626 Double_t angle = seed->GetSnp()*seed->GetSnp();
627 angle = TMath::Sqrt(TMath::Abs(angle/(1.-angle)));
628 Double_t erry2 = clparam->GetError0Par(0,type, z,angle);
630 erry2+=0.5; // edge cluster
633 seed->SetErrorY2(erry2);
637 //calculate look-up table at the beginning
638 // static Bool_t ginit = kFALSE;
639 // static Float_t gnoise1,gnoise2,gnoise3;
640 // static Float_t ggg1[10000];
641 // static Float_t ggg2[10000];
642 // static Float_t ggg3[10000];
643 // static Float_t glandau1[10000];
644 // static Float_t glandau2[10000];
645 // static Float_t glandau3[10000];
647 // static Float_t gcor01[500];
648 // static Float_t gcor02[500];
649 // static Float_t gcorp[500];
653 // if (ginit==kFALSE){
654 // for (Int_t i=1;i<500;i++){
655 // Float_t rsigma = float(i)/100.;
656 // gcor02[i] = TMath::Max(0.78 +TMath::Exp(7.4*(rsigma-1.2)),0.6);
657 // gcor01[i] = TMath::Max(0.72 +TMath::Exp(3.36*(rsigma-1.2)),0.6);
658 // gcorp[i] = TMath::Max(TMath::Power((rsigma+0.5),1.5),1.2);
662 // for (Int_t i=3;i<10000;i++){
666 // Float_t amp = float(i);
667 // Float_t padlength =0.75;
668 // gnoise1 = 0.0004/padlength;
669 // Float_t nel = 0.268*amp;
670 // Float_t nprim = 0.155*amp;
671 // ggg1[i] = fkParam->GetDiffT()*fkParam->GetDiffT()*(2+0.001*nel/(padlength*padlength))/nel;
672 // glandau1[i] = (2.+0.12*nprim)*0.5* (2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
673 // if (glandau1[i]>1) glandau1[i]=1;
674 // glandau1[i]*=padlength*padlength/12.;
678 // gnoise2 = 0.0004/padlength;
680 // nprim = 0.133*amp;
681 // ggg2[i] = fkParam->GetDiffT()*fkParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
682 // glandau2[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
683 // if (glandau2[i]>1) glandau2[i]=1;
684 // glandau2[i]*=padlength*padlength/12.;
689 // gnoise3 = 0.0004/padlength;
691 // nprim = 0.133*amp;
692 // ggg3[i] = fkParam->GetDiffT()*fkParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
693 // glandau3[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
694 // if (glandau3[i]>1) glandau3[i]=1;
695 // glandau3[i]*=padlength*padlength/12.;
703 // Int_t amp = int(TMath::Abs(cl->GetQ()));
705 // seed->SetErrorY2(1.);
709 // Float_t z = TMath::Abs(fkParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
710 // Int_t ctype = cl->GetType();
711 // Float_t padlength= GetPadPitchLength(seed->GetRow());
712 // Double_t angle2 = seed->GetSnp()*seed->GetSnp();
713 // angle2 = angle2/(1-angle2);
715 // //cluster "quality"
716 // Int_t rsigmay = int(100.*cl->GetSigmaY2()/(seed->GetCurrentSigmaY2()));
719 // if (fSectors==fInnerSec){
720 // snoise2 = gnoise1;
721 // res = ggg1[amp]*z+glandau1[amp]*angle2;
722 // if (ctype==0) res *= gcor01[rsigmay];
725 // res*= gcorp[rsigmay];
729 // if (padlength<1.1){
730 // snoise2 = gnoise2;
731 // res = ggg2[amp]*z+glandau2[amp]*angle2;
732 // if (ctype==0) res *= gcor02[rsigmay];
735 // res*= gcorp[rsigmay];
739 // snoise2 = gnoise3;
740 // res = ggg3[amp]*z+glandau3[amp]*angle2;
741 // if (ctype==0) res *= gcor02[rsigmay];
744 // res*= gcorp[rsigmay];
751 // res*=2.4; // overestimate error 2 times
755 // if (res<2*snoise2)
758 // seed->SetErrorY2(res);
766 Double_t AliTPCtrackerMI::ErrZ2(AliTPCseed* seed, const AliTPCclusterMI * cl){
769 // Use calibrated cluster error from OCDB
771 AliTPCClusterParam * clparam = AliTPCcalibDB::Instance()->GetClusterParam();
773 Float_t z = TMath::Abs(fkParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
774 Int_t ctype = cl->GetType();
775 Int_t type = (cl->GetRow()<63) ? 0: (cl->GetRow()>126) ? 1:2;
777 Double_t angle2 = seed->GetSnp()*seed->GetSnp();
778 angle2 = seed->GetTgl()*seed->GetTgl()*(1+angle2/(1-angle2));
779 Double_t angle = TMath::Sqrt(TMath::Abs(angle2));
780 Double_t errz2 = clparam->GetError0Par(1,type, z,angle);
782 errz2+=0.5; // edge cluster
785 seed->SetErrorZ2(errz2);
791 // //seed->SetErrorY2(0.1);
793 // //calculate look-up table at the beginning
794 // static Bool_t ginit = kFALSE;
795 // static Float_t gnoise1,gnoise2,gnoise3;
796 // static Float_t ggg1[10000];
797 // static Float_t ggg2[10000];
798 // static Float_t ggg3[10000];
799 // static Float_t glandau1[10000];
800 // static Float_t glandau2[10000];
801 // static Float_t glandau3[10000];
803 // static Float_t gcor01[1000];
804 // static Float_t gcor02[1000];
805 // static Float_t gcorp[1000];
809 // if (ginit==kFALSE){
810 // for (Int_t i=1;i<1000;i++){
811 // Float_t rsigma = float(i)/100.;
812 // gcor02[i] = TMath::Max(0.81 +TMath::Exp(6.8*(rsigma-1.2)),0.6);
813 // gcor01[i] = TMath::Max(0.72 +TMath::Exp(2.04*(rsigma-1.2)),0.6);
814 // gcorp[i] = TMath::Max(TMath::Power((rsigma+0.5),1.5),1.2);
818 // for (Int_t i=3;i<10000;i++){
822 // Float_t amp = float(i);
823 // Float_t padlength =0.75;
824 // gnoise1 = 0.0004/padlength;
825 // Float_t nel = 0.268*amp;
826 // Float_t nprim = 0.155*amp;
827 // ggg1[i] = fkParam->GetDiffT()*fkParam->GetDiffT()*(2+0.001*nel/(padlength*padlength))/nel;
828 // glandau1[i] = (2.+0.12*nprim)*0.5* (2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
829 // if (glandau1[i]>1) glandau1[i]=1;
830 // glandau1[i]*=padlength*padlength/12.;
834 // gnoise2 = 0.0004/padlength;
836 // nprim = 0.133*amp;
837 // ggg2[i] = fkParam->GetDiffT()*fkParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
838 // glandau2[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
839 // if (glandau2[i]>1) glandau2[i]=1;
840 // glandau2[i]*=padlength*padlength/12.;
845 // gnoise3 = 0.0004/padlength;
847 // nprim = 0.133*amp;
848 // ggg3[i] = fkParam->GetDiffT()*fkParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
849 // glandau3[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
850 // if (glandau3[i]>1) glandau3[i]=1;
851 // glandau3[i]*=padlength*padlength/12.;
859 // Int_t amp = int(TMath::Abs(cl->GetQ()));
861 // seed->SetErrorY2(1.);
865 // Float_t z = TMath::Abs(fkParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
866 // Int_t ctype = cl->GetType();
867 // Float_t padlength= GetPadPitchLength(seed->GetRow());
869 // Double_t angle2 = seed->GetSnp()*seed->GetSnp();
870 // // if (angle2<0.6) angle2 = 0.6;
871 // angle2 = seed->GetTgl()*seed->GetTgl()*(1+angle2/(1-angle2));
873 // //cluster "quality"
874 // Int_t rsigmaz = int(100.*cl->GetSigmaZ2()/(seed->GetCurrentSigmaZ2()));
877 // if (fSectors==fInnerSec){
878 // snoise2 = gnoise1;
879 // res = ggg1[amp]*z+glandau1[amp]*angle2;
880 // if (ctype==0) res *= gcor01[rsigmaz];
883 // res*= gcorp[rsigmaz];
887 // if (padlength<1.1){
888 // snoise2 = gnoise2;
889 // res = ggg2[amp]*z+glandau2[amp]*angle2;
890 // if (ctype==0) res *= gcor02[rsigmaz];
893 // res*= gcorp[rsigmaz];
897 // snoise2 = gnoise3;
898 // res = ggg3[amp]*z+glandau3[amp]*angle2;
899 // if (ctype==0) res *= gcor02[rsigmaz];
902 // res*= gcorp[rsigmaz];
911 // if ((ctype<0) &&<70){
916 // if (res<2*snoise2)
918 // if (res>3) res =3;
919 // seed->SetErrorZ2(res);
927 void AliTPCtrackerMI::RotateToLocal(AliTPCseed *seed)
929 //rotate to track "local coordinata
930 Float_t x = seed->GetX();
931 Float_t y = seed->GetY();
932 Float_t ymax = x*TMath::Tan(0.5*fSectors->GetAlpha());
935 seed->SetRelativeSector((seed->GetRelativeSector()+1) % fN);
936 if (!seed->Rotate(fSectors->GetAlpha()))
938 } else if (y <-ymax) {
939 seed->SetRelativeSector((seed->GetRelativeSector()-1+fN) % fN);
940 if (!seed->Rotate(-fSectors->GetAlpha()))
948 //_____________________________________________________________________________
949 Double_t AliTPCtrackerMI::F1old(Double_t x1,Double_t y1,
950 Double_t x2,Double_t y2,
951 Double_t x3,Double_t y3) const
953 //-----------------------------------------------------------------
954 // Initial approximation of the track curvature
955 //-----------------------------------------------------------------
956 Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
957 Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
958 (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
959 Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
960 (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
962 Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
963 if ( xr*xr+yr*yr<=0.00000000000001) return 100;
964 return -xr*yr/sqrt(xr*xr+yr*yr);
969 //_____________________________________________________________________________
970 Double_t AliTPCtrackerMI::F1(Double_t x1,Double_t y1,
971 Double_t x2,Double_t y2,
972 Double_t x3,Double_t y3) const
974 //-----------------------------------------------------------------
975 // Initial approximation of the track curvature
976 //-----------------------------------------------------------------
982 Double_t det = x3*y2-x2*y3;
983 if (TMath::Abs(det)<1e-16){
987 Double_t u = 0.5* (x2*(x2-x3)+y2*(y2-y3))/det;
988 Double_t x0 = x3*0.5-y3*u;
989 Double_t y0 = y3*0.5+x3*u;
990 Double_t c2 = 1/TMath::Sqrt(x0*x0+y0*y0);
996 Double_t AliTPCtrackerMI::F2(Double_t x1,Double_t y1,
997 Double_t x2,Double_t y2,
998 Double_t x3,Double_t y3) const
1000 //-----------------------------------------------------------------
1001 // Initial approximation of the track curvature
1002 //-----------------------------------------------------------------
1008 Double_t det = x3*y2-x2*y3;
1009 if (TMath::Abs(det)<1e-16)) {
1013 Double_t u = 0.5* (x2*(x2-x3)+y2*(y2-y3))/det;
1014 Double_t x0 = x3*0.5-y3*u;
1015 Double_t y0 = y3*0.5+x3*u;
1016 Double_t c2 = 1/TMath::Sqrt(x0*x0+y0*y0);
1025 //_____________________________________________________________________________
1026 Double_t AliTPCtrackerMI::F2old(Double_t x1,Double_t y1,
1027 Double_t x2,Double_t y2,
1028 Double_t x3,Double_t y3) const
1030 //-----------------------------------------------------------------
1031 // Initial approximation of the track curvature times center of curvature
1032 //-----------------------------------------------------------------
1033 Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
1034 Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
1035 (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
1036 Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
1037 (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
1039 Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
1041 return -a/(d*y1-b)*xr/sqrt(xr*xr+yr*yr);
1044 //_____________________________________________________________________________
1045 Double_t AliTPCtrackerMI::F3(Double_t x1,Double_t y1,
1046 Double_t x2,Double_t y2,
1047 Double_t z1,Double_t z2) const
1049 //-----------------------------------------------------------------
1050 // Initial approximation of the tangent of the track dip angle
1051 //-----------------------------------------------------------------
1052 return (z1 - z2)/sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
1056 Double_t AliTPCtrackerMI::F3n(Double_t x1,Double_t y1,
1057 Double_t x2,Double_t y2,
1058 Double_t z1,Double_t z2, Double_t c) const
1060 //-----------------------------------------------------------------
1061 // Initial approximation of the tangent of the track dip angle
1062 //-----------------------------------------------------------------
1066 //angle1 = (z1-z2)*c/(TMath::ASin(c*x1-ni)-TMath::ASin(c*x2-ni));
1068 Double_t d = TMath::Sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
1069 if (TMath::Abs(d*c*0.5)>1) return 0;
1070 // Double_t angle2 = TMath::ASin(d*c*0.5);
1071 // Double_t angle2 = AliTPCFastMath::FastAsin(d*c*0.5);
1072 Double_t angle2 = (d*c*0.5>0.1)? TMath::ASin(d*c*0.5): AliTPCFastMath::FastAsin(d*c*0.5);
1074 angle2 = (z1-z2)*c/(angle2*2.);
1078 Bool_t AliTPCtrackerMI::GetProlongation(Double_t x1, Double_t x2, Double_t x[5], Double_t &y, Double_t &z) const
1079 {//-----------------------------------------------------------------
1080 // This function find proloncation of a track to a reference plane x=x2.
1081 //-----------------------------------------------------------------
1085 if (TMath::Abs(x[4]*x1 - x[2]) >= 0.999) {
1089 Double_t c1=x[4]*x1 - x[2], r1=TMath::Sqrt((1.-c1)*(1.+c1));
1090 Double_t c2=x[4]*x2 - x[2], r2=TMath::Sqrt((1.-c2)*(1.+c2));
1094 Double_t dy = dx*(c1+c2)/(r1+r2);
1097 Double_t delta = x[4]*dx*(c1+c2)/(c1*r2 + c2*r1);
1099 if (TMath::Abs(delta)>0.01){
1100 dz = x[3]*TMath::ASin(delta)/x[4];
1102 dz = x[3]*AliTPCFastMath::FastAsin(delta)/x[4];
1105 //dz = x[3]*AliTPCFastMath::FastAsin(delta)/x[4];
1113 Int_t AliTPCtrackerMI::LoadClusters (TTree *const tree)
1118 return LoadClusters();
1122 Int_t AliTPCtrackerMI::LoadClusters(const TObjArray *arr)
1125 // load clusters to the memory
1126 AliTPCClustersRow *clrow = 0x0;
1127 Int_t lower = arr->LowerBound();
1128 Int_t entries = arr->GetEntriesFast();
1130 for (Int_t i=lower; i<entries; i++) {
1131 clrow = (AliTPCClustersRow*) arr->At(i);
1132 if(!clrow) continue;
1133 if(!clrow->GetArray()) continue;
1137 fkParam->AdjustSectorRow(clrow->GetID(),sec,row);
1139 for (Int_t icl=0; icl<clrow->GetArray()->GetEntriesFast(); icl++){
1140 Transform((AliTPCclusterMI*)(clrow->GetArray()->At(icl)));
1143 if (clrow->GetArray()->GetEntriesFast()<=0) continue;
1144 AliTPCtrackerRow * tpcrow=0;
1147 tpcrow = &(fInnerSec[sec%fkNIS][row]);
1151 tpcrow = &(fOuterSec[(sec-fkNIS*2)%fkNOS][row]);
1152 left = (sec-fkNIS*2)/fkNOS;
1155 tpcrow->SetN1(clrow->GetArray()->GetEntriesFast());
1156 tpcrow->SetClusters1(new AliTPCclusterMI[tpcrow->GetN1()]);
1157 for (Int_t j=0;j<tpcrow->GetN1();++j)
1158 tpcrow->SetCluster1(j, *(AliTPCclusterMI*)(clrow->GetArray()->At(j)));
1161 tpcrow->SetN2(clrow->GetArray()->GetEntriesFast());
1162 tpcrow->SetClusters2(new AliTPCclusterMI[tpcrow->GetN2()]);
1163 for (Int_t j=0;j<tpcrow->GetN2();++j)
1164 tpcrow->SetCluster2(j,*(AliTPCclusterMI*)(clrow->GetArray()->At(j)));
1174 Int_t AliTPCtrackerMI::LoadClusters(const TClonesArray *arr)
1177 // load clusters to the memory from one
1180 AliTPCclusterMI *clust=0;
1181 Int_t count[72][96] = { {0} , {0} };
1183 // loop over clusters
1184 for (Int_t icl=0; icl<arr->GetEntriesFast(); icl++) {
1185 clust = (AliTPCclusterMI*)arr->At(icl);
1186 if(!clust) continue;
1187 //printf("cluster: det %d, row %d \n", clust->GetDetector(),clust->GetRow());
1189 // transform clusters
1192 // count clusters per pad row
1193 count[clust->GetDetector()][clust->GetRow()]++;
1196 // insert clusters to sectors
1197 for (Int_t icl=0; icl<arr->GetEntriesFast(); icl++) {
1198 clust = (AliTPCclusterMI*)arr->At(icl);
1199 if(!clust) continue;
1201 Int_t sec = clust->GetDetector();
1202 Int_t row = clust->GetRow();
1204 // filter overlapping pad rows needed by HLT
1205 if(sec<fkNIS*2) { //IROCs
1206 if(row == 30) continue;
1209 if(row == 27 || row == 76) continue;
1215 fInnerSec[sec%fkNIS].InsertCluster(clust, count[sec][row], fkParam);
1218 left = (sec-fkNIS*2)/fkNOS;
1219 fOuterSec[(sec-fkNIS*2)%fkNOS].InsertCluster(clust, count[sec][row], fkParam);
1223 // Load functions must be called behind LoadCluster(TClonesArray*)
1225 //LoadOuterSectors();
1226 //LoadInnerSectors();
1232 Int_t AliTPCtrackerMI::LoadClusters()
1235 // load clusters to the memory
1236 AliTPCClustersRow *clrow= new AliTPCClustersRow;
1237 clrow->SetClass("AliTPCclusterMI");
1239 clrow->GetArray()->ExpandCreateFast(10000);
1241 // TTree * tree = fClustersArray.GetTree();
1243 TTree * tree = fInput;
1244 TBranch * br = tree->GetBranch("Segment");
1245 br->SetAddress(&clrow);
1247 Int_t j=Int_t(tree->GetEntries());
1248 for (Int_t i=0; i<j; i++) {
1252 fkParam->AdjustSectorRow(clrow->GetID(),sec,row);
1253 for (Int_t icl=0; icl<clrow->GetArray()->GetEntriesFast(); icl++){
1254 Transform((AliTPCclusterMI*)(clrow->GetArray()->At(icl)));
1257 AliTPCtrackerRow * tpcrow=0;
1260 tpcrow = &(fInnerSec[sec%fkNIS][row]);
1264 tpcrow = &(fOuterSec[(sec-fkNIS*2)%fkNOS][row]);
1265 left = (sec-fkNIS*2)/fkNOS;
1268 tpcrow->SetN1(clrow->GetArray()->GetEntriesFast());
1269 tpcrow->SetClusters1(new AliTPCclusterMI[tpcrow->GetN1()]);
1270 for (Int_t k=0;k<tpcrow->GetN1();++k)
1271 tpcrow->SetCluster1(k, *(AliTPCclusterMI*)(clrow->GetArray()->At(k)));
1274 tpcrow->SetN2(clrow->GetArray()->GetEntriesFast());
1275 tpcrow->SetClusters2(new AliTPCclusterMI[tpcrow->GetN2()]);
1276 for (Int_t k=0;k<tpcrow->GetN2();++k)
1277 tpcrow->SetCluster2(k,*(AliTPCclusterMI*)(clrow->GetArray()->At(k)));
1288 void AliTPCtrackerMI::UnloadClusters()
1291 // unload clusters from the memory
1293 Int_t nrows = fOuterSec->GetNRows();
1294 for (Int_t sec = 0;sec<fkNOS;sec++)
1295 for (Int_t row = 0;row<nrows;row++){
1296 AliTPCtrackerRow* tpcrow = &(fOuterSec[sec%fkNOS][row]);
1298 // if (tpcrow->fClusters1) delete []tpcrow->fClusters1;
1299 // if (tpcrow->fClusters2) delete []tpcrow->fClusters2;
1301 tpcrow->ResetClusters();
1304 nrows = fInnerSec->GetNRows();
1305 for (Int_t sec = 0;sec<fkNIS;sec++)
1306 for (Int_t row = 0;row<nrows;row++){
1307 AliTPCtrackerRow* tpcrow = &(fInnerSec[sec%fkNIS][row]);
1309 // if (tpcrow->fClusters1) delete []tpcrow->fClusters1;
1310 //if (tpcrow->fClusters2) delete []tpcrow->fClusters2;
1312 tpcrow->ResetClusters();
1318 void AliTPCtrackerMI::FillClusterArray(TObjArray* array) const{
1320 // Filling cluster to the array - For visualization purposes
1323 nrows = fOuterSec->GetNRows();
1324 for (Int_t sec = 0;sec<fkNOS;sec++)
1325 for (Int_t row = 0;row<nrows;row++){
1326 AliTPCtrackerRow* tpcrow = &(fOuterSec[sec%fkNOS][row]);
1327 if (!tpcrow) continue;
1328 for (Int_t icl = 0;icl<tpcrow->GetN();icl++){
1329 array->AddLast((TObject*)((*tpcrow)[icl]));
1332 nrows = fInnerSec->GetNRows();
1333 for (Int_t sec = 0;sec<fkNIS;sec++)
1334 for (Int_t row = 0;row<nrows;row++){
1335 AliTPCtrackerRow* tpcrow = &(fInnerSec[sec%fkNIS][row]);
1336 if (!tpcrow) continue;
1337 for (Int_t icl = 0;icl<tpcrow->GetN();icl++){
1338 array->AddLast((TObject*)(*tpcrow)[icl]);
1344 void AliTPCtrackerMI::Transform(AliTPCclusterMI * cluster){
1348 AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ;
1350 AliFatal("Tranformations not in calibDB");
1352 transform->SetCurrentRecoParam((AliTPCRecoParam*)AliTPCReconstructor::GetRecoParam());
1353 Double_t x[3]={cluster->GetRow(),cluster->GetPad(),cluster->GetTimeBin()};
1354 Int_t i[1]={cluster->GetDetector()};
1355 transform->Transform(x,i,0,1);
1356 // if (cluster->GetDetector()%36>17){
1361 // in debug mode check the transformation
1363 if (AliTPCReconstructor::StreamLevel()>2) {
1365 cluster->GetGlobalXYZ(gx);
1366 Int_t event = (fEvent==NULL)? 0: fEvent->GetEventNumberInFile();
1367 TTreeSRedirector &cstream = *fDebugStreamer;
1368 cstream<<"Transform"<<
1379 cluster->SetX(x[0]);
1380 cluster->SetY(x[1]);
1381 cluster->SetZ(x[2]);
1387 //if (!fkParam->IsGeoRead()) fkParam->ReadGeoMatrices();
1388 TGeoHMatrix *mat = fkParam->GetClusterMatrix(cluster->GetDetector());
1390 Double_t pos[3]= {cluster->GetX(),cluster->GetY(),cluster->GetZ()};
1391 Double_t posC[3]={cluster->GetX(),cluster->GetY(),cluster->GetZ()};
1392 if (mat) mat->LocalToMaster(pos,posC);
1394 // chack Loading of Geo matrices from GeoManager - TEMPORARY FIX
1396 cluster->SetX(posC[0]);
1397 cluster->SetY(posC[1]);
1398 cluster->SetZ(posC[2]);
1401 //_____________________________________________________________________________
1402 Int_t AliTPCtrackerMI::LoadOuterSectors() {
1403 //-----------------------------------------------------------------
1404 // This function fills outer TPC sectors with clusters.
1405 //-----------------------------------------------------------------
1406 Int_t nrows = fOuterSec->GetNRows();
1408 for (Int_t sec = 0;sec<fkNOS;sec++)
1409 for (Int_t row = 0;row<nrows;row++){
1410 AliTPCtrackerRow* tpcrow = &(fOuterSec[sec%fkNOS][row]);
1411 Int_t sec2 = sec+2*fkNIS;
1413 Int_t ncl = tpcrow->GetN1();
1415 AliTPCclusterMI *c= (tpcrow->GetCluster1(ncl));
1416 index=(((sec2<<8)+row)<<16)+ncl;
1417 tpcrow->InsertCluster(c,index);
1420 ncl = tpcrow->GetN2();
1422 AliTPCclusterMI *c= (tpcrow->GetCluster2(ncl));
1423 index=((((sec2+fkNOS)<<8)+row)<<16)+ncl;
1424 tpcrow->InsertCluster(c,index);
1427 // write indexes for fast acces
1429 for (Int_t i=0;i<510;i++)
1430 tpcrow->SetFastCluster(i,-1);
1431 for (Int_t i=0;i<tpcrow->GetN();i++){
1432 Int_t zi = Int_t((*tpcrow)[i]->GetZ()+255.);
1433 tpcrow->SetFastCluster(zi,i); // write index
1436 for (Int_t i=0;i<510;i++){
1437 if (tpcrow->GetFastCluster(i)<0)
1438 tpcrow->SetFastCluster(i,last);
1440 last = tpcrow->GetFastCluster(i);
1449 //_____________________________________________________________________________
1450 Int_t AliTPCtrackerMI::LoadInnerSectors() {
1451 //-----------------------------------------------------------------
1452 // This function fills inner TPC sectors with clusters.
1453 //-----------------------------------------------------------------
1454 Int_t nrows = fInnerSec->GetNRows();
1456 for (Int_t sec = 0;sec<fkNIS;sec++)
1457 for (Int_t row = 0;row<nrows;row++){
1458 AliTPCtrackerRow* tpcrow = &(fInnerSec[sec%fkNIS][row]);
1461 Int_t ncl = tpcrow->GetN1();
1463 AliTPCclusterMI *c= (tpcrow->GetCluster1(ncl));
1464 index=(((sec<<8)+row)<<16)+ncl;
1465 tpcrow->InsertCluster(c,index);
1468 ncl = tpcrow->GetN2();
1470 AliTPCclusterMI *c= (tpcrow->GetCluster2(ncl));
1471 index=((((sec+fkNIS)<<8)+row)<<16)+ncl;
1472 tpcrow->InsertCluster(c,index);
1475 // write indexes for fast acces
1477 for (Int_t i=0;i<510;i++)
1478 tpcrow->SetFastCluster(i,-1);
1479 for (Int_t i=0;i<tpcrow->GetN();i++){
1480 Int_t zi = Int_t((*tpcrow)[i]->GetZ()+255.);
1481 tpcrow->SetFastCluster(zi,i); // write index
1484 for (Int_t i=0;i<510;i++){
1485 if (tpcrow->GetFastCluster(i)<0)
1486 tpcrow->SetFastCluster(i,last);
1488 last = tpcrow->GetFastCluster(i);
1500 //_________________________________________________________________________
1501 AliTPCclusterMI *AliTPCtrackerMI::GetClusterMI(Int_t index) const {
1502 //--------------------------------------------------------------------
1503 // Return pointer to a given cluster
1504 //--------------------------------------------------------------------
1505 if (index<0) return 0; // no cluster
1506 Int_t sec=(index&0xff000000)>>24;
1507 Int_t row=(index&0x00ff0000)>>16;
1508 Int_t ncl=(index&0x00007fff)>>00;
1510 const AliTPCtrackerRow * tpcrow=0;
1511 AliTPCclusterMI * clrow =0;
1513 if (sec<0 || sec>=fkNIS*4) {
1514 AliWarning(Form("Wrong sector %d",sec));
1519 AliTPCtrackerSector& tracksec = fInnerSec[sec%fkNIS];
1520 if (tracksec.GetNRows()<=row) return 0;
1521 tpcrow = &(tracksec[row]);
1522 if (tpcrow==0) return 0;
1525 if (tpcrow->GetN1()<=ncl) return 0;
1526 clrow = tpcrow->GetClusters1();
1529 if (tpcrow->GetN2()<=ncl) return 0;
1530 clrow = tpcrow->GetClusters2();
1534 AliTPCtrackerSector& tracksec = fOuterSec[(sec-fkNIS*2)%fkNOS];
1535 if (tracksec.GetNRows()<=row) return 0;
1536 tpcrow = &(tracksec[row]);
1537 if (tpcrow==0) return 0;
1539 if (sec-2*fkNIS<fkNOS) {
1540 if (tpcrow->GetN1()<=ncl) return 0;
1541 clrow = tpcrow->GetClusters1();
1544 if (tpcrow->GetN2()<=ncl) return 0;
1545 clrow = tpcrow->GetClusters2();
1549 return &(clrow[ncl]);
1555 Int_t AliTPCtrackerMI::FollowToNext(AliTPCseed& t, Int_t nr) {
1556 //-----------------------------------------------------------------
1557 // This function tries to find a track prolongation to next pad row
1558 //-----------------------------------------------------------------
1560 Double_t x= GetXrow(nr), ymax=GetMaxY(nr);
1563 AliTPCclusterMI *cl=0;
1564 Int_t tpcindex= t.GetClusterIndex2(nr);
1566 // update current shape info every 5 pad-row
1567 // if ( (nr%5==0) || t.GetNumberOfClusters()<2 || (t.fCurrentSigmaY2<0.0001) ){
1571 if (fIteration>0 && tpcindex>=-1){ //if we have already clusters
1573 if (tpcindex==-1) return 0; //track in dead zone
1575 cl = t.GetClusterPointer(nr);
1576 if ( (cl==0) ) cl = GetClusterMI(tpcindex);
1577 t.SetCurrentClusterIndex1(tpcindex);
1580 Int_t relativesector = ((tpcindex&0xff000000)>>24)%18; // if previously accepted cluster in different sector
1581 Float_t angle = relativesector*fSectors->GetAlpha()+fSectors->GetAlphaShift();
1583 if (angle<-TMath::Pi()) angle += 2*TMath::Pi();
1584 if (angle>=TMath::Pi()) angle -= 2*TMath::Pi();
1586 if (TMath::Abs(angle-t.GetAlpha())>0.001){
1587 Double_t rotation = angle-t.GetAlpha();
1588 t.SetRelativeSector(relativesector);
1589 if (!t.Rotate(rotation)) return 0;
1591 if (!t.PropagateTo(x)) return 0;
1593 t.SetCurrentCluster(cl);
1595 Int_t accept = AcceptCluster(&t,t.GetCurrentCluster());
1596 if ((tpcindex&0x8000)==0) accept =0;
1598 //if founded cluster is acceptible
1599 if (cl->IsUsed(11)) { // id cluster is shared inrease uncertainty
1600 t.SetErrorY2(t.GetErrorY2()+0.03);
1601 t.SetErrorZ2(t.GetErrorZ2()+0.03);
1602 t.SetErrorY2(t.GetErrorY2()*3);
1603 t.SetErrorZ2(t.GetErrorZ2()*3);
1605 t.SetNFoundable(t.GetNFoundable()+1);
1606 UpdateTrack(&t,accept);
1611 if (TMath::Abs(t.GetSnp())>AliTPCReconstructor::GetMaxSnpTracker()) return 0; // cut on angle
1612 if (fIteration>1 && IsFindable(t)){
1613 // not look for new cluster during refitting
1614 t.SetNFoundable(t.GetNFoundable()+1);
1619 // if (TMath::Abs(t.GetSnp())>0.95 || TMath::Abs(x*t.GetC()-t.GetEta())>0.95) return 0;// patch 28 fev 06
1620 Double_t y=t.GetYat(x);
1621 if (TMath::Abs(y)>ymax){
1623 t.SetRelativeSector((t.GetRelativeSector()+1) % fN);
1624 if (!t.Rotate(fSectors->GetAlpha()))
1626 } else if (y <-ymax) {
1627 t.SetRelativeSector((t.GetRelativeSector()-1+fN) % fN);
1628 if (!t.Rotate(-fSectors->GetAlpha()))
1634 if (!t.PropagateTo(x)) {
1635 if (fIteration==0) t.SetRemoval(10);
1639 Double_t z=t.GetZ();
1642 if (!IsActive(t.GetRelativeSector(),nr)) {
1644 t.SetClusterIndex2(nr,-1);
1647 //AliInfo(Form("A - Sector%d phi %f - alpha %f", t.fRelativeSector,y/x, t.GetAlpha()));
1648 Bool_t isActive = IsActive(t.GetRelativeSector(),nr);
1649 Bool_t isActive2 = (nr>=fInnerSec->GetNRows()) ? fOuterSec[t.GetRelativeSector()][nr-fInnerSec->GetNRows()].GetN()>0:fInnerSec[t.GetRelativeSector()][nr].GetN()>0;
1651 if (!isActive || !isActive2) return 0;
1653 const AliTPCtrackerRow &krow=GetRow(t.GetRelativeSector(),nr);
1654 if ( (t.GetSigmaY2()<0) || t.GetSigmaZ2()<0) return 0;
1656 Double_t roadz = 1.;
1658 if (TMath::Abs(TMath::Abs(y)-ymax)<krow.GetDeadZone()){
1660 t.SetClusterIndex2(nr,-1);
1666 // if (TMath::Abs(z)<(AliTPCReconstructor::GetCtgRange()*x+10) && TMath::Abs(z)<fkParam->GetZLength(0) && (TMath::Abs(t.GetSnp())<AliTPCReconstructor::GetMaxSnpTracker()))
1667 t.SetNFoundable(t.GetNFoundable()+1);
1673 // cl = krow.FindNearest2(y+10.,z,roady,roadz,index);
1674 cl = krow.FindNearest2(y,z,roady,roadz,index);
1675 if (cl) t.SetCurrentClusterIndex1(krow.GetIndex(index));
1678 t.SetCurrentCluster(cl);
1680 if (fIteration==2&&cl->IsUsed(10)) return 0;
1681 Int_t accept = AcceptCluster(&t,t.GetCurrentCluster());
1682 if (fIteration==2&&cl->IsUsed(11)) {
1683 t.SetErrorY2(t.GetErrorY2()+0.03);
1684 t.SetErrorZ2(t.GetErrorZ2()+0.03);
1685 t.SetErrorY2(t.GetErrorY2()*3);
1686 t.SetErrorZ2(t.GetErrorZ2()*3);
1689 if (t.fCurrentCluster->IsUsed(10)){
1694 if (t.fNShared>0.7*t.GetNumberOfClusters()) {
1700 if (accept<3) UpdateTrack(&t,accept);
1703 if ( fIteration==0 && t.GetNFoundable()*0.5 > t.GetNumberOfClusters()) t.SetRemoval(10);
1711 //_________________________________________________________________________
1712 Bool_t AliTPCtrackerMI::GetTrackPoint(Int_t index, AliTrackPoint &p ) const
1714 // Get track space point by index
1715 // return false in case the cluster doesn't exist
1716 AliTPCclusterMI *cl = GetClusterMI(index);
1717 if (!cl) return kFALSE;
1718 Int_t sector = (index&0xff000000)>>24;
1719 // Int_t row = (index&0x00ff0000)>>16;
1721 // xyz[0] = fkParam->GetPadRowRadii(sector,row);
1722 xyz[0] = cl->GetX();
1723 xyz[1] = cl->GetY();
1724 xyz[2] = cl->GetZ();
1726 fkParam->AdjustCosSin(sector,cos,sin);
1727 Float_t x = cos*xyz[0]-sin*xyz[1];
1728 Float_t y = cos*xyz[1]+sin*xyz[0];
1730 Float_t sigmaY2 = 0.027*cl->GetSigmaY2();
1731 if (sector < fkParam->GetNInnerSector()) sigmaY2 *= 2.07;
1732 Float_t sigmaZ2 = 0.066*cl->GetSigmaZ2();
1733 if (sector < fkParam->GetNInnerSector()) sigmaZ2 *= 1.77;
1734 cov[0] = sin*sin*sigmaY2;
1735 cov[1] = -sin*cos*sigmaY2;
1737 cov[3] = cos*cos*sigmaY2;
1740 p.SetXYZ(x,y,xyz[2],cov);
1741 AliGeomManager::ELayerID iLayer;
1743 if (sector < fkParam->GetNInnerSector()) {
1744 iLayer = AliGeomManager::kTPC1;
1748 iLayer = AliGeomManager::kTPC2;
1749 idet = sector - fkParam->GetNInnerSector();
1751 UShort_t volid = AliGeomManager::LayerToVolUID(iLayer,idet);
1752 p.SetVolumeID(volid);
1758 Int_t AliTPCtrackerMI::UpdateClusters(AliTPCseed& t, Int_t nr) {
1759 //-----------------------------------------------------------------
1760 // This function tries to find a track prolongation to next pad row
1761 //-----------------------------------------------------------------
1762 t.SetCurrentCluster(0);
1763 t.SetCurrentClusterIndex1(0);
1765 Double_t xt=t.GetX();
1766 Int_t row = GetRowNumber(xt)-1;
1767 Double_t ymax= GetMaxY(nr);
1769 if (row < nr) return 1; // don't prolongate if not information until now -
1770 // if (TMath::Abs(t.GetSnp())>0.9 && t.GetNumberOfClusters()>40. && fIteration!=2) {
1772 // return 0; // not prolongate strongly inclined tracks
1774 // if (TMath::Abs(t.GetSnp())>0.95) {
1776 // return 0; // not prolongate strongly inclined tracks
1777 // }// patch 28 fev 06
1779 Double_t x= GetXrow(nr);
1781 //t.PropagateTo(x+0.02);
1782 //t.PropagateTo(x+0.01);
1783 if (!t.PropagateTo(x)){
1790 if (TMath::Abs(y)>ymax){
1792 t.SetRelativeSector((t.GetRelativeSector()+1) % fN);
1793 if (!t.Rotate(fSectors->GetAlpha()))
1795 } else if (y <-ymax) {
1796 t.SetRelativeSector((t.GetRelativeSector()-1+fN) % fN);
1797 if (!t.Rotate(-fSectors->GetAlpha()))
1800 // if (!t.PropagateTo(x)){
1807 if (TMath::Abs(t.GetSnp())>AliTPCReconstructor::GetMaxSnpTracker()) return 0;
1809 if (!IsActive(t.GetRelativeSector(),nr)) {
1811 t.SetClusterIndex2(nr,-1);
1814 //AliInfo(Form("A - Sector%d phi %f - alpha %f", t.fRelativeSector,y/x, t.GetAlpha()));
1816 AliTPCtrackerRow &krow=GetRow(t.GetRelativeSector(),nr);
1818 if (TMath::Abs(TMath::Abs(y)-ymax)<krow.GetDeadZone()){
1820 t.SetClusterIndex2(nr,-1);
1826 // if (TMath::Abs(t.GetZ())<(AliTPCReconstructor::GetCtgRange()*t.GetX()+10) && (TMath::Abs(t.GetSnp())<AliTPCReconstructor::GetMaxSnpTracker()))
1827 if (IsFindable(t)) t.SetNFoundable(t.GetNFoundable()+1);
1833 if ( (nr%2==0) || t.GetNumberOfClusters()<2){
1834 // t.fCurrentSigmaY = GetSigmaY(&t);
1835 //t.fCurrentSigmaZ = GetSigmaZ(&t);
1839 AliTPCclusterMI *cl=0;
1842 Double_t roady = 1.;
1843 Double_t roadz = 1.;
1847 index = t.GetClusterIndex2(nr);
1848 if ( (index>0) && (index&0x8000)==0){
1849 cl = t.GetClusterPointer(nr);
1850 if ( (cl==0) && (index>0)) cl = GetClusterMI(index);
1851 t.SetCurrentClusterIndex1(index);
1853 t.SetCurrentCluster(cl);
1859 // if (index<0) return 0;
1860 UInt_t uindex = TMath::Abs(index);
1863 //cl = krow.FindNearest2(y+10,z,roady,roadz,uindex);
1864 cl = krow.FindNearest2(y,z,roady,roadz,uindex);
1867 if (cl) t.SetCurrentClusterIndex1(krow.GetIndex(uindex));
1868 t.SetCurrentCluster(cl);
1874 Int_t AliTPCtrackerMI::FollowToNextCluster(AliTPCseed & t, Int_t nr) {
1875 //-----------------------------------------------------------------
1876 // This function tries to find a track prolongation to next pad row
1877 //-----------------------------------------------------------------
1879 //update error according neighborhoud
1881 if (t.GetCurrentCluster()) {
1883 Int_t accept = AcceptCluster(&t,t.GetCurrentCluster());
1885 if (t.GetCurrentCluster()->IsUsed(10)){
1890 t.SetNShared(t.GetNShared()+1);
1891 if (t.GetNShared()>0.7*t.GetNumberOfClusters()) {
1896 if (fIteration>0) accept = 0;
1897 if (accept<3) UpdateTrack(&t,accept);
1901 if ( ( (t.GetSigmaY2()+t.GetSigmaZ2())>0.16)&& t.GetNumberOfClusters()>18) t.SetRemoval(10);
1902 if ( t.GetChi2()/t.GetNumberOfClusters()>6 &&t.GetNumberOfClusters()>18) t.SetRemoval(10);
1904 if (( (t.GetNFoundable()*0.5 > t.GetNumberOfClusters()) || t.GetNoCluster()>15)) t.SetRemoval(10);
1912 //_____________________________________________________________________________
1913 Int_t AliTPCtrackerMI::FollowProlongation(AliTPCseed& t, Int_t rf, Int_t step) {
1914 //-----------------------------------------------------------------
1915 // This function tries to find a track prolongation.
1916 //-----------------------------------------------------------------
1917 Double_t xt=t.GetX();
1919 Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
1920 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
1921 if (alpha < 0. ) alpha += 2.*TMath::Pi();
1923 t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN);
1925 Int_t first = GetRowNumber(xt)-1;
1926 for (Int_t nr= first; nr>=rf; nr-=step) {
1928 if (t.GetKinkIndexes()[0]>0){
1929 for (Int_t i=0;i<3;i++){
1930 Int_t index = t.GetKinkIndexes()[i];
1931 if (index==0) break;
1932 if (index<0) continue;
1934 AliKink * kink = (AliKink*)fEvent->GetKink(index-1);
1936 printf("PROBLEM\n");
1939 Int_t kinkrow = kink->GetTPCRow0()+2+Int_t(0.5/(0.05+kink->GetAngle(2)));
1941 AliExternalTrackParam paramd(t);
1942 kink->SetDaughter(paramd);
1943 kink->SetStatus(2,5);
1950 if (nr==80) t.UpdateReference();
1951 if (nr<fInnerSec->GetNRows())
1952 fSectors = fInnerSec;
1954 fSectors = fOuterSec;
1955 if (FollowToNext(t,nr)==0)
1968 Int_t AliTPCtrackerMI::FollowBackProlongation(AliTPCseed& t, Int_t rf) {
1969 //-----------------------------------------------------------------
1970 // This function tries to find a track prolongation.
1971 //-----------------------------------------------------------------
1973 Double_t xt=t.GetX();
1974 Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
1975 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
1976 if (alpha < 0. ) alpha += 2.*TMath::Pi();
1977 t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN);
1979 Int_t first = t.GetFirstPoint();
1980 if (first<GetRowNumber(xt)+1) first = GetRowNumber(xt)+1;
1982 if (first<0) first=0;
1983 for (Int_t nr=first; nr<=rf; nr++) {
1984 // if ( (TMath::Abs(t.GetSnp())>0.95)) break;//patch 28 fev 06
1985 if (t.GetKinkIndexes()[0]<0){
1986 for (Int_t i=0;i<3;i++){
1987 Int_t index = t.GetKinkIndexes()[i];
1988 if (index==0) break;
1989 if (index>0) continue;
1990 index = TMath::Abs(index);
1991 AliKink * kink = (AliKink*)fEvent->GetKink(index-1);
1993 printf("PROBLEM\n");
1996 Int_t kinkrow = kink->GetTPCRow0()-2-Int_t(0.5/(0.05+kink->GetAngle(2)));
1998 AliExternalTrackParam paramm(t);
1999 kink->SetMother(paramm);
2000 kink->SetStatus(2,1);
2007 if (nr<fInnerSec->GetNRows())
2008 fSectors = fInnerSec;
2010 fSectors = fOuterSec;
2021 Float_t AliTPCtrackerMI::OverlapFactor(AliTPCseed * s1, AliTPCseed * s2, Int_t &sum1, Int_t & sum2)
2029 Float_t dz2 =(s1->GetZ() - s2->GetZ());
2032 Float_t dy2 =TMath::Abs((s1->GetY() - s2->GetY()));
2034 Float_t distance = TMath::Sqrt(dz2+dy2);
2035 if (distance>4.) return 0; // if there are far away - not overlap - to reduce combinatorics
2038 Int_t firstpoint = TMath::Min(s1->GetFirstPoint(),s2->GetFirstPoint());
2039 Int_t lastpoint = TMath::Max(s1->GetLastPoint(),s2->GetLastPoint());
2044 if (firstpoint>lastpoint) {
2045 firstpoint =lastpoint;
2050 for (Int_t i=firstpoint-1;i<lastpoint+1;i++){
2051 if (s1->GetClusterIndex2(i)>0) sum1++;
2052 if (s2->GetClusterIndex2(i)>0) sum2++;
2053 if (s1->GetClusterIndex2(i)==s2->GetClusterIndex2(i) && s1->GetClusterIndex2(i)>0) {
2057 if (sum<5) return 0;
2059 Float_t summin = TMath::Min(sum1+1,sum2+1);
2060 Float_t ratio = (sum+1)/Float_t(summin);
2064 void AliTPCtrackerMI::SignShared(AliTPCseed * s1, AliTPCseed * s2)
2068 Float_t thetaCut = 0.2;//+10.*TMath::Sqrt(s1->GetSigmaTglZ()+ s2->GetSigmaTglZ());
2069 if (TMath::Abs(s1->GetTgl()-s2->GetTgl())>thetaCut) return;
2070 Float_t minCl = TMath::Min(s1->GetNumberOfClusters(),s2->GetNumberOfClusters());
2071 Int_t cutN0 = TMath::Max(5,TMath::Nint(0.1*minCl));
2076 //Int_t firstpoint = TMath::Max(s1->GetFirstPoint(),s2->GetFirstPoint());
2077 //Int_t lastpoint = TMath::Min(s1->GetLastPoint(),s2->GetLastPoint());
2078 Int_t firstpoint = 0;
2079 Int_t lastpoint = 160;
2081 // if (firstpoint>=lastpoint-5) return;;
2083 for (Int_t i=firstpoint;i<lastpoint;i++){
2084 // if ( (s1->GetClusterIndex2(i)&0xFFFF8FFF)==(s2->GetClusterIndex2(i)&0xFFFF8FFF) && s1->GetClusterIndex2(i)>0) {
2085 if ( (s1->GetClusterIndex2(i))==(s2->GetClusterIndex2(i)) && s1->GetClusterIndex2(i)>0) {
2087 s1->SetSharedMapBit(i, kTRUE);
2088 s2->SetSharedMapBit(i, kTRUE);
2090 if (s1->GetClusterIndex2(i)>0)
2091 s1->SetClusterMapBit(i, kTRUE);
2093 if (sumshared>cutN0){
2096 for (Int_t i=firstpoint;i<lastpoint;i++){
2097 // if ( (s1->GetClusterIndex2(i)&0xFFFF8FFF)==(s2->GetClusterIndex2(i)&0xFFFF8FFF) && s1->GetClusterIndex2(i)>0) {
2098 if ( (s1->GetClusterIndex2(i))==(s2->GetClusterIndex2(i)) && s1->GetClusterIndex2(i)>0) {
2099 AliTPCTrackerPoint *p1 = s1->GetTrackPoint(i);
2100 AliTPCTrackerPoint *p2 = s2->GetTrackPoint(i);;
2101 if (s1->IsActive()&&s2->IsActive()){
2102 p1->SetShared(kTRUE);
2103 p2->SetShared(kTRUE);
2109 if (sumshared>cutN0){
2110 for (Int_t i=0;i<4;i++){
2111 if (s1->GetOverlapLabel(3*i)==0){
2112 s1->SetOverlapLabel(3*i, s2->GetLabel());
2113 s1->SetOverlapLabel(3*i+1,sumshared);
2114 s1->SetOverlapLabel(3*i+2,s2->GetUniqueID());
2118 for (Int_t i=0;i<4;i++){
2119 if (s2->GetOverlapLabel(3*i)==0){
2120 s2->SetOverlapLabel(3*i, s1->GetLabel());
2121 s2->SetOverlapLabel(3*i+1,sumshared);
2122 s2->SetOverlapLabel(3*i+2,s1->GetUniqueID());
2129 void AliTPCtrackerMI::SignShared(TObjArray * arr)
2132 //sort trackss according sectors
2134 for (Int_t i=0; i<arr->GetEntriesFast(); i++) {
2135 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2137 //if (pt) RotateToLocal(pt);
2141 arr->Sort(); // sorting according relative sectors
2142 arr->Expand(arr->GetEntries());
2145 Int_t nseed=arr->GetEntriesFast();
2146 for (Int_t i=0; i<nseed; i++) {
2147 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2149 for (Int_t j=0;j<=12;j++){
2150 pt->SetOverlapLabel(j,0);
2153 for (Int_t i=0; i<nseed; i++) {
2154 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2156 if (pt->GetRemoval()>10) continue;
2157 for (Int_t j=i+1; j<nseed; j++){
2158 AliTPCseed *pt2=(AliTPCseed*)arr->UncheckedAt(j);
2159 if (TMath::Abs(pt->GetRelativeSector()-pt2->GetRelativeSector())>1) continue;
2161 if (pt2->GetRemoval()<=10) {
2162 //if ( TMath::Abs(pt->GetRelativeSector()-pt2->GetRelativeSector())>0) break;
2170 void AliTPCtrackerMI::SortTracks(TObjArray * arr, Int_t mode) const
2173 //sort tracks in array according mode criteria
2174 Int_t nseed = arr->GetEntriesFast();
2175 for (Int_t i=0; i<nseed; i++) {
2176 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2187 void AliTPCtrackerMI::RemoveUsed2(TObjArray * arr, Float_t factor1, Float_t factor2, Int_t minimal)
2190 // Loop over all tracks and remove overlaped tracks (with lower quality)
2192 // 1. Unsign clusters
2193 // 2. Sort tracks according quality
2194 // Quality is defined by the number of cluster between first and last points
2196 // 3. Loop over tracks - decreasing quality order
2197 // a.) remove - If the fraction of shared cluster less than factor (1- n or 2)
2198 // b.) remove - If the minimal number of clusters less than minimal and not ITS
2199 // c.) if track accepted - sign clusters
2201 //Called in - AliTPCtrackerMI::Clusters2Tracks()
2202 // - AliTPCtrackerMI::PropagateBack()
2203 // - AliTPCtrackerMI::RefitInward()
2206 // factor1 - factor for constrained
2207 // factor2 - for non constrained tracks
2208 // if (Float_t(shared+1)/Float_t(found+1)>factor) - DELETE
2212 Int_t nseed = arr->GetEntriesFast();
2213 Float_t * quality = new Float_t[nseed];
2214 Int_t * indexes = new Int_t[nseed];
2218 for (Int_t i=0; i<nseed; i++) {
2219 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2224 pt->UpdatePoints(); //select first last max dens points
2225 Float_t * points = pt->GetPoints();
2226 if (points[3]<0.8) quality[i] =-1;
2227 quality[i] = (points[2]-points[0])+pt->GetNumberOfClusters();
2228 //prefer high momenta tracks if overlaps
2229 quality[i] *= TMath::Sqrt(TMath::Abs(pt->Pt())+0.5);
2231 TMath::Sort(nseed,quality,indexes);
2234 for (Int_t itrack=0; itrack<nseed; itrack++) {
2235 Int_t trackindex = indexes[itrack];
2236 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(trackindex);
2239 if (quality[trackindex]<0){
2241 delete arr->RemoveAt(trackindex);
2244 arr->RemoveAt(trackindex);
2250 Int_t first = Int_t(pt->GetPoints()[0]);
2251 Int_t last = Int_t(pt->GetPoints()[2]);
2252 Double_t factor = (pt->GetBConstrain()) ? factor1: factor2;
2254 Int_t found,foundable,shared;
2255 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
2256 // pt->GetClusterStatistic(0,160, found, foundable,shared,kFALSE);
2257 Bool_t itsgold =kFALSE;
2260 if (pt->GetESD()->GetITSclusters(dummy)>4) itsgold= kTRUE;
2264 if (Float_t(shared+1)/Float_t(found+1)>factor){
2265 if (pt->GetKinkIndexes()[0]!=0) continue; //don't remove tracks - part of the kinks
2266 if( AliTPCReconstructor::StreamLevel()>3){
2267 TTreeSRedirector &cstream = *fDebugStreamer;
2268 cstream<<"RemoveUsed"<<
2269 "iter="<<fIteration<<
2273 delete arr->RemoveAt(trackindex);
2276 if (pt->GetNumberOfClusters()<50&&(found-0.5*shared)<minimal){ //remove short tracks
2277 if (pt->GetKinkIndexes()[0]!=0) continue; //don't remove tracks - part of the kinks
2278 if( AliTPCReconstructor::StreamLevel()>3){
2279 TTreeSRedirector &cstream = *fDebugStreamer;
2280 cstream<<"RemoveShort"<<
2281 "iter="<<fIteration<<
2285 delete arr->RemoveAt(trackindex);
2291 //if (sharedfactor>0.4) continue;
2292 if (pt->GetKinkIndexes()[0]>0) continue;
2293 //Remove tracks with undefined properties - seems
2294 if (pt->GetSigmaY2()<kAlmost0) continue; // ? what is the origin ?
2296 for (Int_t i=first; i<last; i++) {
2297 Int_t index=pt->GetClusterIndex2(i);
2298 // if (index<0 || index&0x8000 ) continue;
2299 if (index<0 || index&0x8000 ) continue;
2300 AliTPCclusterMI *c= pt->GetClusterPointer(i);
2307 Info("RemoveUsed","\n*****\nNumber of good tracks after shared removal\t%d\n",fNtracks);
2313 void AliTPCtrackerMI::UnsignClusters()
2316 // loop over all clusters and unsign them
2319 for (Int_t sec=0;sec<fkNIS;sec++){
2320 for (Int_t row=0;row<fInnerSec->GetNRows();row++){
2321 AliTPCclusterMI *cl = fInnerSec[sec][row].GetClusters1();
2322 for (Int_t icl =0;icl< fInnerSec[sec][row].GetN1();icl++)
2323 // if (cl[icl].IsUsed(10))
2325 cl = fInnerSec[sec][row].GetClusters2();
2326 for (Int_t icl =0;icl< fInnerSec[sec][row].GetN2();icl++)
2327 //if (cl[icl].IsUsed(10))
2332 for (Int_t sec=0;sec<fkNOS;sec++){
2333 for (Int_t row=0;row<fOuterSec->GetNRows();row++){
2334 AliTPCclusterMI *cl = fOuterSec[sec][row].GetClusters1();
2335 for (Int_t icl =0;icl< fOuterSec[sec][row].GetN1();icl++)
2336 //if (cl[icl].IsUsed(10))
2338 cl = fOuterSec[sec][row].GetClusters2();
2339 for (Int_t icl =0;icl< fOuterSec[sec][row].GetN2();icl++)
2340 //if (cl[icl].IsUsed(10))
2349 void AliTPCtrackerMI::SignClusters(const TObjArray * arr, Float_t fnumber, Float_t fdensity)
2352 //sign clusters to be "used"
2354 // snumber and sdensity sign number of sigmas - bellow mean value to be accepted
2355 // loop over "primaries"
2369 Int_t nseed = arr->GetEntriesFast();
2370 for (Int_t i=0; i<nseed; i++) {
2371 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2375 if (!(pt->IsActive())) continue;
2376 Float_t dens = pt->GetNumberOfClusters()/Float_t(pt->GetNFoundable());
2377 if ( (dens>0.7) && (pt->GetNumberOfClusters()>70)){
2379 sumdens2+= dens*dens;
2380 sumn += pt->GetNumberOfClusters();
2381 sumn2 += pt->GetNumberOfClusters()*pt->GetNumberOfClusters();
2382 Float_t chi2 = pt->GetChi2()/pt->GetNumberOfClusters();
2385 sumchi2 +=chi2*chi2;
2390 Float_t mdensity = 0.9;
2391 Float_t meann = 130;
2392 Float_t meanchi = 1;
2393 Float_t sdensity = 0.1;
2394 Float_t smeann = 10;
2395 Float_t smeanchi =0.4;
2399 mdensity = sumdens/sum;
2401 meanchi = sumchi/sum;
2403 sdensity = sumdens2/sum-mdensity*mdensity;
2405 sdensity = TMath::Sqrt(sdensity);
2409 smeann = sumn2/sum-meann*meann;
2411 smeann = TMath::Sqrt(smeann);
2415 smeanchi = sumchi2/sum - meanchi*meanchi;
2417 smeanchi = TMath::Sqrt(smeanchi);
2423 //REMOVE SHORT DELTAS or tracks going out of sensitive volume of TPC
2425 for (Int_t i=0; i<nseed; i++) {
2426 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2430 if (pt->GetBSigned()) continue;
2431 if (pt->GetBConstrain()) continue;
2432 //if (!(pt->IsActive())) continue;
2434 Int_t found,foundable,shared;
2435 pt->GetClusterStatistic(0,160,found, foundable,shared);
2436 if (shared/float(found)>0.3) {
2437 if (shared/float(found)>0.9 ){
2438 //delete arr->RemoveAt(i);
2443 Bool_t isok =kFALSE;
2444 if ( (pt->GetNShared()/pt->GetNumberOfClusters()<0.5) &&pt->GetNumberOfClusters()>60)
2446 if ((TMath::Abs(1/pt->GetC())<100.) && (pt->GetNShared()/pt->GetNumberOfClusters()<0.7))
2448 if (TMath::Abs(pt->GetZ()/pt->GetX())>1.1)
2450 if ( (TMath::Abs(pt->GetSnp()>0.7) && pt->GetD(0,0)>60.))
2454 for (Int_t j=0; j<160; ++j) {
2455 Int_t index=pt->GetClusterIndex2(j);
2456 if (index<0) continue;
2457 AliTPCclusterMI *c= pt->GetClusterPointer(j);
2459 //if (!(c->IsUsed(10))) c->Use();
2466 Double_t maxchi = meanchi+2.*smeanchi;
2468 for (Int_t i=0; i<nseed; i++) {
2469 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2473 //if (!(pt->IsActive())) continue;
2474 if (pt->GetBSigned()) continue;
2475 Double_t chi = pt->GetChi2()/pt->GetNumberOfClusters();
2476 if (chi>maxchi) continue;
2479 Float_t dens = pt->GetNumberOfClusters()/Float_t(pt->GetNFoundable());
2481 //sign only tracks with enoug big density at the beginning
2483 if ((pt->GetDensityFirst(40)<0.75) && pt->GetNumberOfClusters()<meann) continue;
2486 Double_t mindens = TMath::Max(double(mdensity-sdensity*fdensity*bfactor),0.65);
2487 Double_t minn = TMath::Max(Int_t(meann-fnumber*smeann*bfactor),50);
2489 // if (pt->fBConstrain) mindens = TMath::Max(mdensity-sdensity*fdensity*bfactor,0.65);
2490 if ( (pt->GetRemoval()==10) && (pt->GetSnp()>0.8)&&(dens>mindens))
2493 if ((dens>mindens && pt->GetNumberOfClusters()>minn) && chi<maxchi ){
2494 //Int_t noc=pt->GetNumberOfClusters();
2495 pt->SetBSigned(kTRUE);
2496 for (Int_t j=0; j<160; ++j) {
2498 Int_t index=pt->GetClusterIndex2(j);
2499 if (index<0) continue;
2500 AliTPCclusterMI *c= pt->GetClusterPointer(j);
2502 // if (!(c->IsUsed(10))) c->Use();
2507 // gLastCheck = nseed;
2515 void AliTPCtrackerMI::StopNotActive(const TObjArray * arr, Int_t row0, Float_t th0, Float_t th1, Float_t th2) const
2517 // stop not active tracks
2518 // take th1 as threshold for number of founded to number of foundable on last 10 active rows
2519 // take th2 as threshold for number of founded to number of foundable on last 20 active rows
2520 Int_t nseed = arr->GetEntriesFast();
2522 for (Int_t i=0; i<nseed; i++) {
2523 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2527 if (!(pt->IsActive())) continue;
2528 StopNotActive(pt,row0,th0, th1,th2);
2534 void AliTPCtrackerMI::StopNotActive(AliTPCseed * seed, Int_t row0, Float_t th0, Float_t th1,
2537 // stop not active tracks
2538 // take th1 as threshold for number of founded to number of foundable on last 10 active rows
2539 // take th2 as threshold for number of founded to number of foundable on last 20 active rows
2542 Int_t foundable = 0;
2543 Int_t maxindex = seed->GetLastPoint(); //last foundable row
2544 if (seed->GetNFoundable()*th0 > seed->GetNumberOfClusters()) {
2545 seed->Desactivate(10) ;
2549 for (Int_t i=row0; i<maxindex; i++){
2550 Int_t index = seed->GetClusterIndex2(i);
2551 if (index!=-1) foundable++;
2553 if (foundable<=30) sumgood1++;
2554 if (foundable<=50) {
2561 if (foundable>=30.){
2562 if (sumgood1<(th1*30.)) seed->Desactivate(10);
2565 if (sumgood2<(th2*50.)) seed->Desactivate(10);
2569 Int_t AliTPCtrackerMI::RefitInward(AliESDEvent *event)
2572 // back propagation of ESD tracks
2575 const Int_t kMaxFriendTracks=2000;
2579 //PrepareForProlongation(fSeeds,1);
2580 PropagateForward2(fSeeds);
2581 RemoveUsed2(fSeeds,0.4,0.4,20);
2583 TObjArray arraySeed(fSeeds->GetEntries());
2584 for (Int_t i=0;i<fSeeds->GetEntries();i++) {
2585 arraySeed.AddAt(fSeeds->At(i),i);
2587 SignShared(&arraySeed);
2588 // FindCurling(fSeeds, event,2); // find multi found tracks
2589 FindSplitted(fSeeds, event,2); // find multi found tracks
2590 if (AliTPCReconstructor::StreamLevel()>5) FindMultiMC(fSeeds, fEvent,2); // find multi found tracks
2593 Int_t nseed = fSeeds->GetEntriesFast();
2594 for (Int_t i=0;i<nseed;i++){
2595 AliTPCseed * seed = (AliTPCseed*) fSeeds->UncheckedAt(i);
2596 if (!seed) continue;
2597 if (seed->GetKinkIndex(0)>0) UpdateKinkQualityD(seed); // update quality informations for kinks
2598 AliESDtrack *esd=event->GetTrack(i);
2600 if (seed->GetNumberOfClusters()<60 && seed->GetNumberOfClusters()<(esd->GetTPCclusters(0) -5)*0.8){
2601 AliExternalTrackParam paramIn;
2602 AliExternalTrackParam paramOut;
2603 Int_t ncl = seed->RefitTrack(seed,¶mIn,¶mOut);
2604 if (AliTPCReconstructor::StreamLevel()>2) {
2605 (*fDebugStreamer)<<"RecoverIn"<<
2609 "pout.="<<¶mOut<<
2614 seed->Set(paramIn.GetX(),paramIn.GetAlpha(),paramIn.GetParameter(),paramIn.GetCovariance());
2615 seed->SetNumberOfClusters(ncl);
2619 seed->PropagateTo(fkParam->GetInnerRadiusLow());
2620 seed->UpdatePoints();
2621 AddCovariance(seed);
2623 seed->CookdEdx(0.02,0.6);
2624 CookLabel(seed,0.1); //For comparison only
2625 esd->SetTPCClusterMap(seed->GetClusterMap());
2626 esd->SetTPCSharedMap(seed->GetSharedMap());
2628 if (AliTPCReconstructor::StreamLevel()>1 && seed!=0&&esd!=0) {
2629 TTreeSRedirector &cstream = *fDebugStreamer;
2636 if (seed->GetNumberOfClusters()>15){
2637 esd->UpdateTrackParams(seed,AliESDtrack::kTPCrefit);
2638 esd->SetTPCPoints(seed->GetPoints());
2639 esd->SetTPCPointsF(seed->GetNFoundable());
2640 Int_t ndedx = seed->GetNCDEDX(0);
2641 Float_t sdedx = seed->GetSDEDX(0);
2642 Float_t dedx = seed->GetdEdx();
2643 esd->SetTPCsignal(dedx, sdedx, ndedx);
2645 // add seed to the esd track in Calib level
2647 Bool_t storeFriend = gRandom->Rndm()<(kMaxFriendTracks)/Float_t(nseed);
2648 if (AliTPCReconstructor::StreamLevel()>0 &&storeFriend){
2649 AliTPCseed * seedCopy = new AliTPCseed(*seed, kTRUE);
2650 esd->AddCalibObject(seedCopy);
2655 //printf("problem\n");
2658 //FindKinks(fSeeds,event);
2659 Info("RefitInward","Number of refitted tracks %d",ntracks);
2664 Int_t AliTPCtrackerMI::PropagateBack(AliESDEvent *event)
2667 // back propagation of ESD tracks
2673 PropagateBack(fSeeds);
2674 RemoveUsed2(fSeeds,0.4,0.4,20);
2675 //FindCurling(fSeeds, fEvent,1);
2676 FindSplitted(fSeeds, event,1); // find multi found tracks
2677 if (AliTPCReconstructor::StreamLevel()>5) FindMultiMC(fSeeds, fEvent,1); // find multi found tracks
2680 Int_t nseed = fSeeds->GetEntriesFast();
2682 for (Int_t i=0;i<nseed;i++){
2683 AliTPCseed * seed = (AliTPCseed*) fSeeds->UncheckedAt(i);
2684 if (!seed) continue;
2685 if (seed->GetKinkIndex(0)<0) UpdateKinkQualityM(seed); // update quality informations for kinks
2686 seed->UpdatePoints();
2687 AddCovariance(seed);
2688 AliESDtrack *esd=event->GetTrack(i);
2689 if (seed->GetNumberOfClusters()<60 && seed->GetNumberOfClusters()<(esd->GetTPCclusters(0) -5)*0.8){
2690 AliExternalTrackParam paramIn;
2691 AliExternalTrackParam paramOut;
2692 Int_t ncl = seed->RefitTrack(seed,¶mIn,¶mOut);
2693 if (AliTPCReconstructor::StreamLevel()>2) {
2694 (*fDebugStreamer)<<"RecoverBack"<<
2698 "pout.="<<¶mOut<<
2703 seed->Set(paramOut.GetX(),paramOut.GetAlpha(),paramOut.GetParameter(),paramOut.GetCovariance());
2704 seed->SetNumberOfClusters(ncl);
2707 seed->CookdEdx(0.02,0.6);
2708 CookLabel(seed,0.1); //For comparison only
2709 if (seed->GetNumberOfClusters()>15){
2710 esd->UpdateTrackParams(seed,AliESDtrack::kTPCout);
2711 esd->SetTPCPoints(seed->GetPoints());
2712 esd->SetTPCPointsF(seed->GetNFoundable());
2713 Int_t ndedx = seed->GetNCDEDX(0);
2714 Float_t sdedx = seed->GetSDEDX(0);
2715 Float_t dedx = seed->GetdEdx();
2716 esd->SetTPCsignal(dedx, sdedx, ndedx);
2718 Int_t eventnumber = event->GetEventNumberInFile();// patch 28 fev 06
2719 // This is most likely NOT the event number you'd like to use. It has nothing to do with the 'real' event number
2720 if (AliTPCReconstructor::StreamLevel()>1 && esd) {
2721 (*fDebugStreamer)<<"Cback"<<
2724 "EventNrInFile="<<eventnumber<<
2729 //FindKinks(fSeeds,event);
2730 Info("PropagateBack","Number of back propagated tracks %d",ntracks);
2737 void AliTPCtrackerMI::DeleteSeeds()
2742 Int_t nseed = fSeeds->GetEntriesFast();
2743 for (Int_t i=0;i<nseed;i++){
2744 AliTPCseed * seed = (AliTPCseed*)fSeeds->At(i);
2745 if (seed) delete fSeeds->RemoveAt(i);
2752 void AliTPCtrackerMI::ReadSeeds(AliESDEvent *const event, Int_t direction)
2755 //read seeds from the event
2757 Int_t nentr=event->GetNumberOfTracks();
2759 Info("ReadSeeds", "Number of ESD tracks: %d\n", nentr);
2764 fSeeds = new TObjArray(nentr);
2768 for (Int_t i=0; i<nentr; i++) {
2769 AliESDtrack *esd=event->GetTrack(i);
2770 ULong_t status=esd->GetStatus();
2771 if (!(status&AliESDtrack::kTPCin)) continue;
2772 AliTPCtrack t(*esd);
2773 t.SetNumberOfClusters(0);
2774 // AliTPCseed *seed = new AliTPCseed(t,t.GetAlpha());
2775 AliTPCseed *seed = new AliTPCseed(t/*,t.GetAlpha()*/);
2776 seed->SetUniqueID(esd->GetID());
2777 AddCovariance(seed); //add systematic ucertainty
2778 for (Int_t ikink=0;ikink<3;ikink++) {
2779 Int_t index = esd->GetKinkIndex(ikink);
2780 seed->GetKinkIndexes()[ikink] = index;
2781 if (index==0) continue;
2782 index = TMath::Abs(index);
2783 AliESDkink * kink = fEvent->GetKink(index-1);
2784 if (kink&&esd->GetKinkIndex(ikink)<0){
2785 if ((status & AliESDtrack::kTRDrefit) != 0) kink->SetStatus(1,2);
2786 if ((status & AliESDtrack::kITSout) != 0) kink->SetStatus(1,0);
2788 if (kink&&esd->GetKinkIndex(ikink)>0){
2789 if ((status & AliESDtrack::kTRDrefit) != 0) kink->SetStatus(1,6);
2790 if ((status & AliESDtrack::kITSout) != 0) kink->SetStatus(1,4);
2794 if (((status&AliESDtrack::kITSout)==0)&&(direction==1)) seed->ResetCovariance(10.);
2795 if ( direction ==2 &&(status & AliESDtrack::kTRDrefit) == 0 ) seed->ResetCovariance(10.);
2796 //if ( direction ==2 && ((status & AliESDtrack::kTPCout) == 0) ) {
2797 // fSeeds->AddAt(0,i);
2801 if ( direction ==2 &&(status & AliESDtrack::kTRDrefit) > 0 ) {
2802 Double_t par0[5],par1[5],alpha,x;
2803 esd->GetInnerExternalParameters(alpha,x,par0);
2804 esd->GetExternalParameters(x,par1);
2805 Double_t delta1 = TMath::Abs(par0[4]-par1[4])/(0.000000001+TMath::Abs(par0[4]+par1[4]));
2806 Double_t delta2 = TMath::Abs(par0[3]-par1[3]);
2808 if (esd->GetTRDncls()>0) trdchi2 = esd->GetTRDchi2()/esd->GetTRDncls();
2809 //reset covariance if suspicious
2810 if ( (delta1>0.1) || (delta2>0.006) ||trdchi2>7.)
2811 seed->ResetCovariance(10.);
2816 // rotate to the local coordinate system
2818 fSectors=fInnerSec; fN=fkNIS;
2819 Double_t alpha=seed->GetAlpha() - fSectors->GetAlphaShift();
2820 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
2821 if (alpha < 0. ) alpha += 2.*TMath::Pi();
2822 Int_t ns=Int_t(alpha/fSectors->GetAlpha())%fN;
2823 alpha =ns*fSectors->GetAlpha() + fSectors->GetAlphaShift();
2824 if (alpha<-TMath::Pi()) alpha += 2*TMath::Pi();
2825 if (alpha>=TMath::Pi()) alpha -= 2*TMath::Pi();
2826 alpha-=seed->GetAlpha();
2827 if (!seed->Rotate(alpha)) {
2833 if (esd->GetKinkIndex(0)<=0){
2834 for (Int_t irow=0;irow<160;irow++){
2835 Int_t index = seed->GetClusterIndex2(irow);
2838 AliTPCclusterMI * cl = GetClusterMI(index);
2839 seed->SetClusterPointer(irow,cl);
2841 if ((index & 0x8000)==0){
2842 cl->Use(10); // accepted cluster
2844 cl->Use(6); // close cluster not accepted
2847 Info("ReadSeeds","Not found cluster");
2852 fSeeds->AddAt(seed,i);
2858 //_____________________________________________________________________________
2859 void AliTPCtrackerMI::MakeSeeds3(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, Float_t cuts[4],
2860 Float_t deltay, Int_t ddsec) {
2861 //-----------------------------------------------------------------
2862 // This function creates track seeds.
2863 // SEEDING WITH VERTEX CONSTRAIN
2864 //-----------------------------------------------------------------
2865 // cuts[0] - fP4 cut
2866 // cuts[1] - tan(phi) cut
2867 // cuts[2] - zvertex cut
2868 // cuts[3] - fP3 cut
2876 Double_t x[5], c[15];
2877 // Int_t di = i1-i2;
2879 AliTPCseed * seed = new AliTPCseed();
2880 Double_t alpha=fSectors->GetAlpha(), shift=fSectors->GetAlphaShift();
2881 Double_t cs=cos(alpha), sn=sin(alpha);
2883 // Double_t x1 =fOuterSec->GetX(i1);
2884 //Double_t xx2=fOuterSec->GetX(i2);
2886 Double_t x1 =GetXrow(i1);
2887 Double_t xx2=GetXrow(i2);
2889 Double_t x3=GetX(), y3=GetY(), z3=GetZ();
2891 Int_t imiddle = (i2+i1)/2; //middle pad row index
2892 Double_t xm = GetXrow(imiddle); // radius of middle pad-row
2893 const AliTPCtrackerRow& krm=GetRow(sec,imiddle); //middle pad -row
2897 const AliTPCtrackerRow& kr1=GetRow(ns,i1);
2898 Double_t ymax = GetMaxY(i1)-kr1.GetDeadZone()-1.5;
2899 Double_t ymaxm = GetMaxY(imiddle)-kr1.GetDeadZone()-1.5;
2902 // change cut on curvature if it can't reach this layer
2903 // maximal curvature set to reach it
2904 Double_t dvertexmax = TMath::Sqrt((x1-x3)*(x1-x3)+(ymax+5-y3)*(ymax+5-y3));
2905 if (dvertexmax*0.5*cuts[0]>0.85){
2906 cuts[0] = 0.85/(dvertexmax*0.5+1.);
2908 Double_t r2min = 1/(cuts[0]*cuts[0]); //minimal square of radius given by cut
2911 if (deltay>0) ddsec = 0;
2912 // loop over clusters
2913 for (Int_t is=0; is < kr1; is++) {
2915 if (kr1[is]->IsUsed(10)) continue;
2916 Double_t y1=kr1[is]->GetY(), z1=kr1[is]->GetZ();
2917 //if (TMath::Abs(y1)>ymax) continue;
2919 if (deltay>0 && TMath::Abs(ymax-TMath::Abs(y1))> deltay ) continue; // seed only at the edge
2921 // find possible directions
2922 Float_t anglez = (z1-z3)/(x1-x3);
2923 Float_t extraz = z1 - anglez*(x1-xx2); // extrapolated z
2926 //find rotation angles relative to line given by vertex and point 1
2927 Double_t dvertex2 = (x1-x3)*(x1-x3)+(y1-y3)*(y1-y3);
2928 Double_t dvertex = TMath::Sqrt(dvertex2);
2929 Double_t angle13 = TMath::ATan((y1-y3)/(x1-x3));
2930 Double_t cs13 = cos(-angle13), sn13 = sin(-angle13);
2933 // loop over 2 sectors
2939 Double_t dddz1=0; // direction of delta inclination in z axis
2946 for (Int_t dsec = dsec1; dsec<=dsec2;dsec++){
2947 Int_t sec2 = sec + dsec;
2949 // AliTPCtrackerRow& kr2 = fOuterSec[(sec2+fkNOS)%fkNOS][i2];
2950 //AliTPCtrackerRow& kr2m = fOuterSec[(sec2+fkNOS)%fkNOS][imiddle];
2951 AliTPCtrackerRow& kr2 = GetRow((sec2+fkNOS)%fkNOS,i2);
2952 AliTPCtrackerRow& kr2m = GetRow((sec2+fkNOS)%fkNOS,imiddle);
2953 Int_t index1 = TMath::Max(kr2.Find(extraz-0.6-dddz1*TMath::Abs(z1)*0.05)-1,0);
2954 Int_t index2 = TMath::Min(kr2.Find(extraz+0.6+dddz2*TMath::Abs(z1)*0.05)+1,kr2);
2956 // rotation angles to p1-p3
2957 Double_t cs13r = cos(-angle13+dsec*alpha)/dvertex, sn13r = sin(-angle13+dsec*alpha)/dvertex;
2958 Double_t x2, y2, z2;
2960 // Double_t dymax = maxangle*TMath::Abs(x1-xx2);
2963 Double_t dxx0 = (xx2-x3)*cs13r;
2964 Double_t dyy0 = (xx2-x3)*sn13r;
2965 for (Int_t js=index1; js < index2; js++) {
2966 const AliTPCclusterMI *kcl = kr2[js];
2967 if (kcl->IsUsed(10)) continue;
2969 //calcutate parameters
2971 Double_t yy0 = dyy0 +(kcl->GetY()-y3)*cs13r;
2973 if (TMath::Abs(yy0)<0.000001) continue;
2974 Double_t xx0 = dxx0 -(kcl->GetY()-y3)*sn13r;
2975 Double_t y0 = 0.5*(xx0*xx0+yy0*yy0-xx0)/yy0;
2976 Double_t r02 = (0.25+y0*y0)*dvertex2;
2977 //curvature (radius) cut
2978 if (r02<r2min) continue;
2982 Double_t c0 = 1/TMath::Sqrt(r02);
2986 //Double_t dfi0 = 2.*TMath::ASin(dvertex*c0*0.5);
2987 //Double_t dfi1 = 2.*TMath::ASin(TMath::Sqrt(yy0*yy0+(1-xx0)*(1-xx0))*dvertex*c0*0.5);
2988 Double_t dfi0 = 2.*AliTPCFastMath::FastAsin(dvertex*c0*0.5);
2989 Double_t dfi1 = 2.*AliTPCFastMath::FastAsin(TMath::Sqrt(yy0*yy0+(1-xx0)*(1-xx0))*dvertex*c0*0.5);
2992 Double_t z0 = kcl->GetZ();
2993 Double_t zzzz2 = z1-(z1-z3)*dfi1/dfi0;
2994 if (TMath::Abs(zzzz2-z0)>0.5) continue;
2997 Double_t dip = (z1-z0)*c0/dfi1;
2998 Double_t x0 = (0.5*cs13+y0*sn13)*dvertex*c0;
3009 x2= xx2*cs-y2*sn*dsec;
3010 y2=+xx2*sn*dsec+y2*cs;
3020 // do we have cluster at the middle ?
3022 GetProlongation(x1,xm,x,ym,zm);
3024 AliTPCclusterMI * cm=0;
3025 if (TMath::Abs(ym)-ymaxm<0){
3026 cm = krm.FindNearest2(ym,zm,1.0,0.6,dummy);
3027 if ((!cm) || (cm->IsUsed(10))) {
3032 // rotate y1 to system 0
3033 // get state vector in rotated system
3034 Double_t yr1 = (-0.5*sn13+y0*cs13)*dvertex*c0;
3035 Double_t xr2 = x0*cs+yr1*sn*dsec;
3036 Double_t xr[5]={kcl->GetY(),kcl->GetZ(), xr2, dip, c0};
3038 GetProlongation(xx2,xm,xr,ym,zm);
3039 if (TMath::Abs(ym)-ymaxm<0){
3040 cm = kr2m.FindNearest2(ym,zm,1.0,0.6,dummy);
3041 if ((!cm) || (cm->IsUsed(10))) {
3051 dym = ym - cm->GetY();
3052 dzm = zm - cm->GetZ();
3059 Double_t sy1=kr1[is]->GetSigmaY2()*2., sz1=kr1[is]->GetSigmaZ2()*2.;
3060 Double_t sy2=kcl->GetSigmaY2()*2., sz2=kcl->GetSigmaZ2()*2.;
3061 //Double_t sy3=400*3./12., sy=0.1, sz=0.1;
3062 Double_t sy3=25000*x[4]*x[4]+0.1, sy=0.1, sz=0.1;
3063 //Double_t sy3=25000*x[4]*x[4]*60+0.5, sy=0.1, sz=0.1;
3065 Double_t f40=(F1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
3066 Double_t f42=(F1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
3067 Double_t f43=(F1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
3068 Double_t f20=(F2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
3069 Double_t f22=(F2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
3070 Double_t f23=(F2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
3072 Double_t f30=(F3(x1,y1+sy,x2,y2,z1,z2)-x[3])/sy;
3073 Double_t f31=(F3(x1,y1,x2,y2,z1+sz,z2)-x[3])/sz;
3074 Double_t f32=(F3(x1,y1,x2,y2+sy,z1,z2)-x[3])/sy;
3075 Double_t f34=(F3(x1,y1,x2,y2,z1,z2+sz)-x[3])/sz;
3079 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3080 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3081 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3082 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3083 c[13]=f30*sy1*f40+f32*sy2*f42;
3084 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3086 // if (!BuildSeed(kr1[is],kcl,0,x1,x2,x3,x,c)) continue;
3088 UInt_t index=kr1.GetIndex(is);
3089 seed->~AliTPCseed(); // this does not set the pointer to 0...
3090 AliTPCseed *track=new(seed) AliTPCseed(x1, ns*alpha+shift, x, c, index);
3092 track->SetIsSeeding(kTRUE);
3093 track->SetSeed1(i1);
3094 track->SetSeed2(i2);
3095 track->SetSeedType(3);
3099 FollowProlongation(*track, (i1+i2)/2,1);
3100 Int_t foundable,found,shared;
3101 track->GetClusterStatistic((i1+i2)/2,i1, found, foundable, shared, kTRUE);
3102 if ((found<0.55*foundable) || shared>0.5*found || (track->GetSigmaY2()+track->GetSigmaZ2())>0.5){
3104 seed->~AliTPCseed();
3110 FollowProlongation(*track, i2,1);
3114 track->SetBConstrain(1);
3115 // track->fLastPoint = i1+fInnerSec->GetNRows(); // first cluster in track position
3116 track->SetLastPoint(i1); // first cluster in track position
3117 track->SetFirstPoint(track->GetLastPoint());
3119 if (track->GetNumberOfClusters()<(i1-i2)*0.5 ||
3120 track->GetNumberOfClusters() < track->GetNFoundable()*0.6 ||
3121 track->GetNShared()>0.4*track->GetNumberOfClusters() ) {
3123 seed->~AliTPCseed();
3127 // Z VERTEX CONDITION
3128 Double_t zv, bz=GetBz();
3129 if ( !track->GetZAt(0.,bz,zv) ) continue;
3130 if (TMath::Abs(zv-z3)>cuts[2]) {
3131 FollowProlongation(*track, TMath::Max(i2-20,0));
3132 if ( !track->GetZAt(0.,bz,zv) ) continue;
3133 if (TMath::Abs(zv-z3)>cuts[2]){
3134 FollowProlongation(*track, TMath::Max(i2-40,0));
3135 if ( !track->GetZAt(0.,bz,zv) ) continue;
3136 if (TMath::Abs(zv-z3)>cuts[2] &&(track->GetNumberOfClusters() > track->GetNFoundable()*0.7)){
3137 // make seed without constrain
3138 AliTPCseed * track2 = MakeSeed(track,0.2,0.5,1.);
3139 FollowProlongation(*track2, i2,1);
3140 track2->SetBConstrain(kFALSE);
3141 track2->SetSeedType(1);
3142 arr->AddLast(track2);
3144 seed->~AliTPCseed();
3149 seed->~AliTPCseed();
3156 track->SetSeedType(0);
3157 arr->AddLast(track);
3158 seed = new AliTPCseed;
3160 // don't consider other combinations
3161 if (track->GetNumberOfClusters() > track->GetNFoundable()*0.8)
3167 Info("MakeSeeds3","\nSeeding statistic:\t%d\t%d\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin,nout1,nout2);
3173 void AliTPCtrackerMI::MakeSeeds5(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, Float_t cuts[4],
3178 //-----------------------------------------------------------------
3179 // This function creates track seeds.
3180 //-----------------------------------------------------------------
3181 // cuts[0] - fP4 cut
3182 // cuts[1] - tan(phi) cut
3183 // cuts[2] - zvertex cut
3184 // cuts[3] - fP3 cut
3194 Double_t x[5], c[15];
3196 // make temporary seed
3197 AliTPCseed * seed = new AliTPCseed;
3198 Double_t alpha=fOuterSec->GetAlpha(), shift=fOuterSec->GetAlphaShift();
3199 // Double_t cs=cos(alpha), sn=sin(alpha);
3204 Double_t x1 = GetXrow(i1-1);
3205 const AliTPCtrackerRow& kr1=GetRow(sec,i1-1);
3206 Double_t y1max = GetMaxY(i1-1)-kr1.GetDeadZone()-1.5;
3208 Double_t x1p = GetXrow(i1);
3209 const AliTPCtrackerRow& kr1p=GetRow(sec,i1);
3211 Double_t x1m = GetXrow(i1-2);
3212 const AliTPCtrackerRow& kr1m=GetRow(sec,i1-2);
3215 //last 3 padrow for seeding
3216 AliTPCtrackerRow& kr3 = GetRow((sec+fkNOS)%fkNOS,i1-7);
3217 Double_t x3 = GetXrow(i1-7);
3218 // Double_t y3max= GetMaxY(i1-7)-kr3.fDeadZone-1.5;
3220 AliTPCtrackerRow& kr3p = GetRow((sec+fkNOS)%fkNOS,i1-6);
3221 Double_t x3p = GetXrow(i1-6);
3223 AliTPCtrackerRow& kr3m = GetRow((sec+fkNOS)%fkNOS,i1-8);
3224 Double_t x3m = GetXrow(i1-8);
3229 Int_t im = i1-4; //middle pad row index
3230 Double_t xm = GetXrow(im); // radius of middle pad-row
3231 const AliTPCtrackerRow& krm=GetRow(sec,im); //middle pad -row
3232 // Double_t ymmax = GetMaxY(im)-kr1.fDeadZone-1.5;
3235 Double_t deltax = x1-x3;
3236 Double_t dymax = deltax*cuts[1];
3237 Double_t dzmax = deltax*cuts[3];
3239 // loop over clusters
3240 for (Int_t is=0; is < kr1; is++) {
3242 if (kr1[is]->IsUsed(10)) continue;
3243 Double_t y1=kr1[is]->GetY(), z1=kr1[is]->GetZ();
3245 if (deltay>0 && TMath::Abs(y1max-TMath::Abs(y1))> deltay ) continue; // seed only at the edge
3247 Int_t index1 = TMath::Max(kr3.Find(z1-dzmax)-1,0);
3248 Int_t index2 = TMath::Min(kr3.Find(z1+dzmax)+1,kr3);
3254 for (Int_t js=index1; js < index2; js++) {
3255 const AliTPCclusterMI *kcl = kr3[js];
3256 if (kcl->IsUsed(10)) continue;
3258 // apply angular cuts
3259 if (TMath::Abs(y1-y3)>dymax) continue;
3262 if (TMath::Abs(z1-z3)>dzmax) continue;
3264 Double_t angley = (y1-y3)/(x1-x3);
3265 Double_t anglez = (z1-z3)/(x1-x3);
3267 Double_t erry = TMath::Abs(angley)*(x1-x1m)*0.5+0.5;
3268 Double_t errz = TMath::Abs(anglez)*(x1-x1m)*0.5+0.5;
3270 Double_t yyym = angley*(xm-x1)+y1;
3271 Double_t zzzm = anglez*(xm-x1)+z1;
3273 const AliTPCclusterMI *kcm = krm.FindNearest2(yyym,zzzm,erry,errz,index);
3275 if (kcm->IsUsed(10)) continue;
3277 erry = TMath::Abs(angley)*(x1-x1m)*0.4+0.5;
3278 errz = TMath::Abs(anglez)*(x1-x1m)*0.4+0.5;
3285 // look around first
3286 const AliTPCclusterMI *kc1m = kr1m.FindNearest2(angley*(x1m-x1)+y1,
3292 if (kc1m->IsUsed(10)) used++;
3294 const AliTPCclusterMI *kc1p = kr1p.FindNearest2(angley*(x1p-x1)+y1,
3300 if (kc1p->IsUsed(10)) used++;
3302 if (used>1) continue;
3303 if (found<1) continue;
3307 const AliTPCclusterMI *kc3m = kr3m.FindNearest2(angley*(x3m-x3)+y3,
3313 if (kc3m->IsUsed(10)) used++;
3317 const AliTPCclusterMI *kc3p = kr3p.FindNearest2(angley*(x3p-x3)+y3,
3323 if (kc3p->IsUsed(10)) used++;
3327 if (used>1) continue;
3328 if (found<3) continue;
3338 x[4]=F1(x1,y1,x2,y2,x3,y3);
3339 //if (TMath::Abs(x[4]) >= cuts[0]) continue;
3342 x[2]=F2(x1,y1,x2,y2,x3,y3);
3345 x[3]=F3n(x1,y1,x2,y2,z1,z2,x[4]);
3346 //if (TMath::Abs(x[3]) > cuts[3]) continue;
3350 Double_t sy1=0.1, sz1=0.1;
3351 Double_t sy2=0.1, sz2=0.1;
3352 Double_t sy3=0.1, sy=0.1, sz=0.1;
3354 Double_t f40=(F1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
3355 Double_t f42=(F1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
3356 Double_t f43=(F1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
3357 Double_t f20=(F2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
3358 Double_t f22=(F2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
3359 Double_t f23=(F2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
3361 Double_t f30=(F3(x1,y1+sy,x2,y2,z1,z2)-x[3])/sy;
3362 Double_t f31=(F3(x1,y1,x2,y2,z1+sz,z2)-x[3])/sz;
3363 Double_t f32=(F3(x1,y1,x2,y2+sy,z1,z2)-x[3])/sy;
3364 Double_t f34=(F3(x1,y1,x2,y2,z1,z2+sz)-x[3])/sz;
3368 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3369 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3370 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3371 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3372 c[13]=f30*sy1*f40+f32*sy2*f42;
3373 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3375 // if (!BuildSeed(kr1[is],kcl,0,x1,x2,x3,x,c)) continue;
3377 index=kr1.GetIndex(is);
3378 seed->~AliTPCseed();
3379 AliTPCseed *track=new(seed) AliTPCseed(x1, sec*alpha+shift, x, c, index);
3381 track->SetIsSeeding(kTRUE);
3384 FollowProlongation(*track, i1-7,1);
3385 if (track->GetNumberOfClusters() < track->GetNFoundable()*0.75 ||
3386 track->GetNShared()>0.6*track->GetNumberOfClusters() || ( track->GetSigmaY2()+ track->GetSigmaZ2())>0.6){
3388 seed->~AliTPCseed();
3394 FollowProlongation(*track, i2,1);
3395 track->SetBConstrain(0);
3396 track->SetLastPoint(i1+fInnerSec->GetNRows()); // first cluster in track position
3397 track->SetFirstPoint(track->GetLastPoint());
3399 if (track->GetNumberOfClusters()<(i1-i2)*0.5 ||
3400 track->GetNumberOfClusters()<track->GetNFoundable()*0.7 ||
3401 track->GetNShared()>2. || track->GetChi2()/track->GetNumberOfClusters()>6 || ( track->GetSigmaY2()+ track->GetSigmaZ2())>0.5 ) {
3403 seed->~AliTPCseed();
3408 FollowProlongation(*track, TMath::Max(i2-10,0),1);
3409 AliTPCseed * track2 = MakeSeed(track,0.2,0.5,0.9);
3410 FollowProlongation(*track2, i2,1);
3411 track2->SetBConstrain(kFALSE);
3412 track2->SetSeedType(4);
3413 arr->AddLast(track2);
3415 seed->~AliTPCseed();
3419 //arr->AddLast(track);
3420 //seed = new AliTPCseed;
3426 Info("MakeSeeds5","\nSeeding statiistic:\t%d\t%d\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin,nout1,nout2,nout3);
3432 //_____________________________________________________________________________
3433 void AliTPCtrackerMI::MakeSeeds2(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, Float_t */*cuts[4]*/,
3434 Float_t deltay, Bool_t /*bconstrain*/) {
3435 //-----------------------------------------------------------------
3436 // This function creates track seeds - without vertex constraint
3437 //-----------------------------------------------------------------
3438 // cuts[0] - fP4 cut - not applied
3439 // cuts[1] - tan(phi) cut
3440 // cuts[2] - zvertex cut - not applied
3441 // cuts[3] - fP3 cut
3451 Double_t alpha=fOuterSec->GetAlpha(), shift=fOuterSec->GetAlphaShift();
3452 // Double_t cs=cos(alpha), sn=sin(alpha);
3453 Int_t row0 = (i1+i2)/2;
3454 Int_t drow = (i1-i2)/2;
3455 const AliTPCtrackerRow& kr0=fSectors[sec][row0];
3456 AliTPCtrackerRow * kr=0;
3458 AliTPCpolyTrack polytrack;
3459 Int_t nclusters=fSectors[sec][row0];
3460 AliTPCseed * seed = new AliTPCseed;
3465 for (Int_t is=0; is < nclusters; is++) { //LOOP over clusters
3467 Int_t nfoundable =0;
3468 for (Int_t iter =1; iter<2; iter++){ //iterations
3469 const AliTPCtrackerRow& krm=fSectors[sec][row0-iter];
3470 const AliTPCtrackerRow& krp=fSectors[sec][row0+iter];
3471 const AliTPCclusterMI * cl= kr0[is];
3473 if (cl->IsUsed(10)) {
3479 Double_t x = kr0.GetX();
3480 // Initialization of the polytrack
3485 Double_t y0= cl->GetY();
3486 Double_t z0= cl->GetZ();
3490 Double_t ymax = fSectors->GetMaxY(row0)-kr0.GetDeadZone()-1.5;
3491 if (deltay>0 && TMath::Abs(ymax-TMath::Abs(y0))> deltay ) continue; // seed only at the edge
3493 erry = (0.5)*cl->GetSigmaY2()/TMath::Sqrt(cl->GetQ())*6;
3494 errz = (0.5)*cl->GetSigmaZ2()/TMath::Sqrt(cl->GetQ())*6;
3495 polytrack.AddPoint(x,y0,z0,erry, errz);
3498 if (cl->IsUsed(10)) sumused++;
3501 Float_t roady = (5*TMath::Sqrt(cl->GetSigmaY2()+0.2)+1.)*iter;
3502 Float_t roadz = (5*TMath::Sqrt(cl->GetSigmaZ2()+0.2)+1.)*iter;
3505 AliTPCclusterMI * cl1 = krm.FindNearest(y0,z0,roady,roadz);
3506 if (cl1 && TMath::Abs(ymax-TMath::Abs(y0))) {
3507 erry = (0.5)*cl1->GetSigmaY2()/TMath::Sqrt(cl1->GetQ())*3;
3508 errz = (0.5)*cl1->GetSigmaZ2()/TMath::Sqrt(cl1->GetQ())*3;
3509 if (cl1->IsUsed(10)) sumused++;
3510 polytrack.AddPoint(x,cl1->GetY(),cl1->GetZ(),erry,errz);
3514 AliTPCclusterMI * cl2 = krp.FindNearest(y0,z0,roady,roadz);
3516 erry = (0.5)*cl2->GetSigmaY2()/TMath::Sqrt(cl2->GetQ())*3;
3517 errz = (0.5)*cl2->GetSigmaZ2()/TMath::Sqrt(cl2->GetQ())*3;
3518 if (cl2->IsUsed(10)) sumused++;
3519 polytrack.AddPoint(x,cl2->GetY(),cl2->GetZ(),erry,errz);
3522 if (sumused>0) continue;
3524 polytrack.UpdateParameters();
3530 nfoundable = polytrack.GetN();
3531 nfound = nfoundable;
3533 for (Int_t ddrow = iter+1; ddrow<drow;ddrow++){
3534 Float_t maxdist = 0.8*(1.+3./(ddrow));
3535 for (Int_t delta = -1;delta<=1;delta+=2){
3536 Int_t row = row0+ddrow*delta;
3537 kr = &(fSectors[sec][row]);
3538 Double_t xn = kr->GetX();
3539 Double_t ymax1 = fSectors->GetMaxY(row)-kr->GetDeadZone()-1.5;
3540 polytrack.GetFitPoint(xn,yn,zn);
3541 if (TMath::Abs(yn)>ymax1) continue;
3543 AliTPCclusterMI * cln = kr->FindNearest(yn,zn,roady,roadz);
3545 Float_t dist = TMath::Sqrt( (yn-cln->GetY())*(yn-cln->GetY())+(zn-cln->GetZ())*(zn-cln->GetZ()));
3548 erry = (dist+0.3)*cln->GetSigmaY2()/TMath::Sqrt(cln->GetQ())*(1.+1./(ddrow));
3549 errz = (dist+0.3)*cln->GetSigmaZ2()/TMath::Sqrt(cln->GetQ())*(1.+1./(ddrow));
3550 if (cln->IsUsed(10)) {
3551 // printf("used\n");
3559 polytrack.AddPoint(xn,cln->GetY(),cln->GetZ(),erry, errz);
3564 if ( (sumused>3) || (sumused>0.5*nfound) || (nfound<0.6*nfoundable)) break;
3565 polytrack.UpdateParameters();
3568 if ( (sumused>3) || (sumused>0.5*nfound)) {
3569 //printf("sumused %d\n",sumused);
3574 polytrack.GetFitDerivation(kr0.GetX(),dy,dz);
3575 AliTPCpolyTrack track2;
3577 polytrack.Refit(track2,0.5+TMath::Abs(dy)*0.3,0.4+TMath::Abs(dz)*0.3);
3578 if (track2.GetN()<0.5*nfoundable) continue;
3581 if ((nfound>0.6*nfoundable) &&( nfoundable>0.4*(i1-i2))) {
3583 // test seed with and without constrain
3584 for (Int_t constrain=0; constrain<=0;constrain++){
3585 // add polytrack candidate
3587 Double_t x[5], c[15];
3588 Double_t x1,x2,x3,y1,y2,y3,z1,z2,z3;
3589 track2.GetBoundaries(x3,x1);
3591 track2.GetFitPoint(x1,y1,z1);
3592 track2.GetFitPoint(x2,y2,z2);
3593 track2.GetFitPoint(x3,y3,z3);
3595 //is track pointing to the vertex ?
3598 polytrack.GetFitPoint(x0,y0,z0);
3611 x[4]=F1(x1,y1,x2,y2,x3,y3);
3613 // if (TMath::Abs(x[4]) >= cuts[0]) continue; //
3614 x[2]=F2(x1,y1,x2,y2,x3,y3);
3616 //if (TMath::Abs(x[4]*x1-x[2]) >= cuts[1]) continue;
3617 //x[3]=F3(x1,y1,x2,y2,z1,z2);
3618 x[3]=F3n(x1,y1,x3,y3,z1,z3,x[4]);
3619 //if (TMath::Abs(x[3]) > cuts[3]) continue;
3622 Double_t sy =0.1, sz =0.1;
3623 Double_t sy1=0.02, sz1=0.02;
3624 Double_t sy2=0.02, sz2=0.02;
3628 sy3=25000*x[4]*x[4]+0.1, sy=0.1, sz=0.1;
3631 Double_t f40=(F1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
3632 Double_t f42=(F1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
3633 Double_t f43=(F1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
3634 Double_t f20=(F2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
3635 Double_t f22=(F2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
3636 Double_t f23=(F2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
3638 Double_t f30=(F3(x1,y1+sy,x3,y3,z1,z3)-x[3])/sy;
3639 Double_t f31=(F3(x1,y1,x3,y3,z1+sz,z3)-x[3])/sz;
3640 Double_t f32=(F3(x1,y1,x3,y3+sy,z1,z3)-x[3])/sy;
3641 Double_t f34=(F3(x1,y1,x3,y3,z1,z3+sz)-x[3])/sz;
3646 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3647 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3648 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3649 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3650 c[13]=f30*sy1*f40+f32*sy2*f42;
3651 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3653 //Int_t row1 = fSectors->GetRowNumber(x1);
3654 Int_t row1 = GetRowNumber(x1);
3658 seed->~AliTPCseed();
3659 AliTPCseed *track=new(seed) AliTPCseed(x1,sec*alpha+shift,x,c,index);
3660 track->SetIsSeeding(kTRUE);
3661 Int_t rc=FollowProlongation(*track, i2);
3662 if (constrain) track->SetBConstrain(1);
3664 track->SetBConstrain(0);
3665 track->SetLastPoint(row1+fInnerSec->GetNRows()); // first cluster in track position
3666 track->SetFirstPoint(track->GetLastPoint());
3668 if (rc==0 || track->GetNumberOfClusters()<(i1-i2)*0.5 ||
3669 track->GetNumberOfClusters() < track->GetNFoundable()*0.6 ||
3670 track->GetNShared()>0.4*track->GetNumberOfClusters()) {
3673 seed->~AliTPCseed();
3676 arr->AddLast(track);
3677 seed = new AliTPCseed;
3681 } // if accepted seed
3684 Info("MakeSeeds2","\nSeeding statiistic:\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin3);
3690 AliTPCseed *AliTPCtrackerMI::MakeSeed(AliTPCseed *const track, Float_t r0, Float_t r1, Float_t r2)
3694 //reseed using track points
3695 Int_t p0 = int(r0*track->GetNumberOfClusters()); // point 0
3696 Int_t p1 = int(r1*track->GetNumberOfClusters());
3697 Int_t p2 = int(r2*track->GetNumberOfClusters()); // last point
3699 Double_t x0[3],x1[3],x2[3];
3700 for (Int_t i=0;i<3;i++){
3706 // find track position at given ratio of the length
3707 Int_t sec0=0, sec1=0, sec2=0;
3710 for (Int_t i=0;i<160;i++){
3711 if (track->GetClusterPointer(i)){
3713 AliTPCTrackerPoint *trpoint =track->GetTrackPoint(i);
3714 if ( (index<p0) || x0[0]<0 ){
3715 if (trpoint->GetX()>1){
3716 clindex = track->GetClusterIndex2(i);
3718 x0[0] = trpoint->GetX();
3719 x0[1] = trpoint->GetY();
3720 x0[2] = trpoint->GetZ();
3721 sec0 = ((clindex&0xff000000)>>24)%18;
3726 if ( (index<p1) &&(trpoint->GetX()>1)){
3727 clindex = track->GetClusterIndex2(i);
3729 x1[0] = trpoint->GetX();
3730 x1[1] = trpoint->GetY();
3731 x1[2] = trpoint->GetZ();
3732 sec1 = ((clindex&0xff000000)>>24)%18;
3735 if ( (index<p2) &&(trpoint->GetX()>1)){
3736 clindex = track->GetClusterIndex2(i);
3738 x2[0] = trpoint->GetX();
3739 x2[1] = trpoint->GetY();
3740 x2[2] = trpoint->GetZ();
3741 sec2 = ((clindex&0xff000000)>>24)%18;
3748 Double_t alpha, cs,sn, xx2,yy2;
3750 alpha = (sec1-sec2)*fSectors->GetAlpha();
3751 cs = TMath::Cos(alpha);
3752 sn = TMath::Sin(alpha);
3753 xx2= x1[0]*cs-x1[1]*sn;
3754 yy2= x1[0]*sn+x1[1]*cs;
3758 alpha = (sec0-sec2)*fSectors->GetAlpha();
3759 cs = TMath::Cos(alpha);
3760 sn = TMath::Sin(alpha);
3761 xx2= x0[0]*cs-x0[1]*sn;
3762 yy2= x0[0]*sn+x0[1]*cs;
3768 Double_t x[5],c[15];
3772 x[4]=F1(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]);
3773 // if (x[4]>1) return 0;
3774 x[2]=F2(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]);
3775 x[3]=F3n(x2[0],x2[1],x0[0],x0[1],x2[2],x0[2],x[4]);
3776 //if (TMath::Abs(x[3]) > 2.2) return 0;
3777 //if (TMath::Abs(x[2]) > 1.99) return 0;
3779 Double_t sy =0.1, sz =0.1;
3781 Double_t sy1=0.02+track->GetSigmaY2(), sz1=0.02+track->GetSigmaZ2();
3782 Double_t sy2=0.01+track->GetSigmaY2(), sz2=0.01+track->GetSigmaZ2();
3783 Double_t sy3=0.01+track->GetSigmaY2();
3785 Double_t f40=(F1(x2[0],x2[1]+sy,x1[0],x1[1],x0[0],x0[1])-x[4])/sy;
3786 Double_t f42=(F1(x2[0],x2[1],x1[0],x1[1]+sy,x0[0],x0[1])-x[4])/sy;
3787 Double_t f43=(F1(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]+sy)-x[4])/sy;
3788 Double_t f20=(F2(x2[0],x2[1]+sy,x1[0],x1[1],x0[0],x0[1])-x[2])/sy;
3789 Double_t f22=(F2(x2[0],x2[1],x1[0],x1[1]+sy,x0[0],x0[1])-x[2])/sy;
3790 Double_t f23=(F2(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]+sy)-x[2])/sy;
3792 Double_t f30=(F3(x2[0],x2[1]+sy,x0[0],x0[1],x2[2],x0[2])-x[3])/sy;
3793 Double_t f31=(F3(x2[0],x2[1],x0[0],x0[1],x2[2]+sz,x0[2])-x[3])/sz;
3794 Double_t f32=(F3(x2[0],x2[1],x0[0],x0[1]+sy,x2[2],x0[2])-x[3])/sy;
3795 Double_t f34=(F3(x2[0],x2[1],x0[0],x0[1],x2[2],x0[2]+sz)-x[3])/sz;
3800 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3801 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3802 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3803 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3804 c[13]=f30*sy1*f40+f32*sy2*f42;
3805 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3807 // Int_t row1 = fSectors->GetRowNumber(x2[0]);
3808 AliTPCseed *seed=new AliTPCseed(x2[0], sec2*fSectors->GetAlpha()+fSectors->GetAlphaShift(), x, c, 0);
3809 // Double_t y0,z0,y1,z1, y2,z2;
3810 //seed->GetProlongation(x0[0],y0,z0);
3811 // seed->GetProlongation(x1[0],y1,z1);
3812 //seed->GetProlongation(x2[0],y2,z2);
3814 seed->SetLastPoint(pp2);
3815 seed->SetFirstPoint(pp2);
3822 AliTPCseed *AliTPCtrackerMI::ReSeed(const AliTPCseed *track, Float_t r0, Float_t r1, Float_t r2)
3826 //reseed using founded clusters
3828 // Find the number of clusters
3829 Int_t nclusters = 0;
3830 for (Int_t irow=0;irow<160;irow++){
3831 if (track->GetClusterIndex(irow)>0) nclusters++;
3835 ipos[0] = TMath::Max(int(r0*nclusters),0); // point 0 cluster
3836 ipos[1] = TMath::Min(int(r1*nclusters),nclusters-1); //
3837 ipos[2] = TMath::Min(int(r2*nclusters),nclusters-1); // last point
3841 Int_t row[3],sec[3]={0,0,0};
3843 // find track row position at given ratio of the length
3845 for (Int_t irow=0;irow<160;irow++){
3846 if (track->GetClusterIndex2(irow)<0) continue;
3848 for (Int_t ipoint=0;ipoint<3;ipoint++){
3849 if (index<=ipos[ipoint]) row[ipoint] = irow;
3853 //Get cluster and sector position
3854 for (Int_t ipoint=0;ipoint<3;ipoint++){
3855 Int_t clindex = track->GetClusterIndex2(row[ipoint]);
3856 AliTPCclusterMI * cl = GetClusterMI(clindex);
3859 // AliTPCclusterMI * cl = GetClusterMI(clindex);
3862 sec[ipoint] = ((clindex&0xff000000)>>24)%18;
3863 xyz[ipoint][0] = GetXrow(row[ipoint]);
3864 xyz[ipoint][1] = cl->GetY();
3865 xyz[ipoint][2] = cl->GetZ();
3869 // Calculate seed state vector and covariance matrix
3871 Double_t alpha, cs,sn, xx2,yy2;
3873 alpha = (sec[1]-sec[2])*fSectors->GetAlpha();
3874 cs = TMath::Cos(alpha);
3875 sn = TMath::Sin(alpha);
3876 xx2= xyz[1][0]*cs-xyz[1][1]*sn;
3877 yy2= xyz[1][0]*sn+xyz[1][1]*cs;
3881 alpha = (sec[0]-sec[2])*fSectors->GetAlpha();
3882 cs = TMath::Cos(alpha);
3883 sn = TMath::Sin(alpha);
3884 xx2= xyz[0][0]*cs-xyz[0][1]*sn;
3885 yy2= xyz[0][0]*sn+xyz[0][1]*cs;
3891 Double_t x[5],c[15];
3895 x[4]=F1(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
3896 x[2]=F2(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
3897 x[3]=F3n(xyz[2][0],xyz[2][1],xyz[0][0],xyz[0][1],xyz[2][2],xyz[0][2],x[4]);
3899 Double_t sy =0.1, sz =0.1;
3901 Double_t sy1=0.2, sz1=0.2;
3902 Double_t sy2=0.2, sz2=0.2;
3905 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;
3906 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;
3907 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;
3908 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;
3909 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;
3910 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;
3912 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;
3913 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;
3914 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;
3915 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;
3920 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3921 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3922 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3923 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3924 c[13]=f30*sy1*f40+f32*sy2*f42;
3925 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3927 // Int_t row1 = fSectors->GetRowNumber(xyz[2][0]);
3928 AliTPCseed *seed=new AliTPCseed(xyz[2][0], sec[2]*fSectors->GetAlpha()+fSectors->GetAlphaShift(), x, c, 0);
3929 seed->SetLastPoint(row[2]);
3930 seed->SetFirstPoint(row[2]);
3935 AliTPCseed *AliTPCtrackerMI::ReSeed(AliTPCseed *track,Int_t r0, Bool_t forward)
3939 //reseed using founded clusters
3942 Int_t row[3]={0,0,0};
3943 Int_t sec[3]={0,0,0};
3945 // forward direction
3947 for (Int_t irow=r0;irow<160;irow++){
3948 if (track->GetClusterIndex(irow)>0){
3953 for (Int_t irow=160;irow>r0;irow--){
3954 if (track->GetClusterIndex(irow)>0){
3959 for (Int_t irow=row[2]-15;irow>row[0];irow--){
3960 if (track->GetClusterIndex(irow)>0){
3968 for (Int_t irow=0;irow<r0;irow++){
3969 if (track->GetClusterIndex(irow)>0){
3974 for (Int_t irow=r0;irow>0;irow--){
3975 if (track->GetClusterIndex(irow)>0){
3980 for (Int_t irow=row[2]-15;irow>row[0];irow--){
3981 if (track->GetClusterIndex(irow)>0){
3988 if ((row[2]-row[0])<20) return 0;
3989 if (row[1]==0) return 0;
3992 //Get cluster and sector position
3993 for (Int_t ipoint=0;ipoint<3;ipoint++){
3994 Int_t clindex = track->GetClusterIndex2(row[ipoint]);
3995 AliTPCclusterMI * cl = GetClusterMI(clindex);
3998 // AliTPCclusterMI * cl = GetClusterMI(clindex);
4001 sec[ipoint] = ((clindex&0xff000000)>>24)%18;
4002 xyz[ipoint][0] = GetXrow(row[ipoint]);
4003 AliTPCTrackerPoint * point = track->GetTrackPoint(row[ipoint]);
4004 if (point&&ipoint<2){
4006 xyz[ipoint][1] = point->GetY();
4007 xyz[ipoint][2] = point->GetZ();
4010 xyz[ipoint][1] = cl->GetY();
4011 xyz[ipoint][2] = cl->GetZ();
4018 // Calculate seed state vector and covariance matrix
4020 Double_t alpha, cs,sn, xx2,yy2;
4022 alpha = (sec[1]-sec[2])*fSectors->GetAlpha();
4023 cs = TMath::Cos(alpha);
4024 sn = TMath::Sin(alpha);
4025 xx2= xyz[1][0]*cs-xyz[1][1]*sn;
4026 yy2= xyz[1][0]*sn+xyz[1][1]*cs;
4030 alpha = (sec[0]-sec[2])*fSectors->GetAlpha();
4031 cs = TMath::Cos(alpha);
4032 sn = TMath::Sin(alpha);
4033 xx2= xyz[0][0]*cs-xyz[0][1]*sn;
4034 yy2= xyz[0][0]*sn+xyz[0][1]*cs;
4040 Double_t x[5],c[15];
4044 x[4]=F1(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
4045 x[2]=F2(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
4046 x[3]=F3n(xyz[2][0],xyz[2][1],xyz[0][0],xyz[0][1],xyz[2][2],xyz[0][2],x[4]);
4048 Double_t sy =0.1, sz =0.1;
4050 Double_t sy1=0.2, sz1=0.2;
4051 Double_t sy2=0.2, sz2=0.2;
4054 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;
4055 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;
4056 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;
4057 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;
4058 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;
4059 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;
4061 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;
4062 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;
4063 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;
4064 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;
4069 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
4070 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
4071 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
4072 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
4073 c[13]=f30*sy1*f40+f32*sy2*f42;
4074 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
4076 // Int_t row1 = fSectors->GetRowNumber(xyz[2][0]);
4077 AliTPCseed *seed=new AliTPCseed(xyz[2][0], sec[2]*fSectors->GetAlpha()+fSectors->GetAlphaShift(), x, c, 0);
4078 seed->SetLastPoint(row[2]);
4079 seed->SetFirstPoint(row[2]);
4080 for (Int_t i=row[0];i<row[2];i++){
4081 seed->SetClusterIndex(i, track->GetClusterIndex(i));
4089 void AliTPCtrackerMI::FindMultiMC(const TObjArray * array, AliESDEvent */*esd*/, Int_t iter)
4092 // find multi tracks - THIS FUNCTION IS ONLY FOR DEBUG PURPOSES
4094 // Use AliTPCReconstructor::StreamLevel()>2 if you want to tune parameters - cuts
4096 // Two reasons to have multiple find tracks
4097 // 1. Curling tracks can be find more than once
4098 // 2. Splitted tracks
4099 // a.) Multiple seeding to increase tracking efficiency - (~ 100% reached)
4100 // b.) Edge effect on the sector boundaries
4103 // Algorithm done in 2 phases - because of CPU consumption
4104 // it is n^2 algorithm - for lead-lead 20000x20000 combination are investigated
4106 // Algorihm for curling tracks sign:
4107 // 1 phase -makes a very rough fast cuts to minimize combinatorics
4108 // a.) opposite sign
4109 // b.) one of the tracks - not pointing to the primary vertex -
4110 // c.) delta tan(theta)
4112 // 2 phase - calculates DCA between tracks - time consument
4117 // General cuts - for splitted tracks and for curling tracks
4119 const Float_t kMaxdPhi = 0.2; // maximal distance in phi
4121 // Curling tracks cuts
4126 Int_t nentries = array->GetEntriesFast();
4127 AliHelix *helixes = new AliHelix[nentries];
4128 Float_t *xm = new Float_t[nentries];
4129 Float_t *dz0 = new Float_t[nentries];
4130 Float_t *dz1 = new Float_t[nentries];
4136 // Find track COG in x direction - point with best defined parameters
4138 for (Int_t i=0;i<nentries;i++){
4139 AliTPCseed* track = (AliTPCseed*)array->At(i);
4140 if (!track) continue;
4141 track->SetCircular(0);
4142 new (&helixes[i]) AliHelix(*track);
4146 track->GetDZ(GetX(),GetY(),GetZ(),GetBz(),dz);
4149 for (Int_t icl=0; icl<160; icl++){
4150 AliTPCclusterMI * cl = track->GetClusterPointer(icl);
4156 if (ncl>0) xm[i]/=Float_t(ncl);
4159 for (Int_t i0=0;i0<nentries;i0++){
4160 AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
4161 if (!track0) continue;
4162 Float_t xc0 = helixes[i0].GetHelix(6);
4163 Float_t yc0 = helixes[i0].GetHelix(7);
4164 Float_t r0 = helixes[i0].GetHelix(8);
4165 Float_t rc0 = TMath::Sqrt(xc0*xc0+yc0*yc0);
4166 Float_t fi0 = TMath::ATan2(yc0,xc0);
4168 for (Int_t i1=i0+1;i1<nentries;i1++){
4169 AliTPCseed * track1 = (AliTPCseed*)array->At(i1);
4170 if (!track1) continue;
4171 Int_t lab0=track0->GetLabel();
4172 Int_t lab1=track1->GetLabel();
4173 if (TMath::Abs(lab0)!=TMath::Abs(lab1)) continue;
4175 Float_t xc1 = helixes[i1].GetHelix(6);
4176 Float_t yc1 = helixes[i1].GetHelix(7);
4177 Float_t r1 = helixes[i1].GetHelix(8);
4178 Float_t rc1 = TMath::Sqrt(xc1*xc1+yc1*yc1);
4179 Float_t fi1 = TMath::ATan2(yc1,xc1);
4181 Float_t dfi = fi0-fi1;
4184 if (dfi>1.5*TMath::Pi()) dfi-=TMath::Pi(); // take care about edge effect
4185 if (dfi<-1.5*TMath::Pi()) dfi+=TMath::Pi(); //
4186 if (TMath::Abs(dfi)>kMaxdPhi&&helixes[i0].GetHelix(4)*helixes[i1].GetHelix(4)<0){
4188 // if short tracks with undefined sign
4189 fi1 = -TMath::ATan2(yc1,-xc1);
4192 Float_t dtheta = TMath::Abs(track0->GetTgl()-track1->GetTgl())<TMath::Abs(track0->GetTgl()+track1->GetTgl())? track0->GetTgl()-track1->GetTgl():track0->GetTgl()+track1->GetTgl();
4195 // debug stream to tune "fast cuts"
4197 Double_t dist[3]; // distance at X
4198 Double_t mdist[3]={0,0,0}; // mean distance X+-40cm
4199 track0->GetDistance(track1,0.5*(xm[i0]+xm[i1])-40.,dist,AliTracker::GetBz());
4200 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[i]);
4201 track0->GetDistance(track1,0.5*(xm[i0]+xm[i1])+40.,dist,AliTracker::GetBz());
4202 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[i]);
4203 track0->GetDistance(track1,0.5*(xm[i0]+xm[i1]),dist,AliTracker::GetBz());
4204 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[i]);
4205 for (Int_t i=0;i<3;i++) mdist[i]*=0.33333;
4209 for (Int_t icl=0; icl<160; icl++){
4210 AliTPCclusterMI * cl0 = track0->GetClusterPointer(icl);
4211 AliTPCclusterMI * cl1 = track1->GetClusterPointer(icl);
4214 if (cl0==cl1) sums++;
4218 if (AliTPCReconstructor::StreamLevel()>5) {
4219 TTreeSRedirector &cstream = *fDebugStreamer;
4224 "Tr0.="<<track0<< // seed0
4225 "Tr1.="<<track1<< // seed1
4226 "h0.="<<&helixes[i0]<<
4227 "h1.="<<&helixes[i1]<<
4229 "sum="<<sum<< //the sum of rows with cl in both
4230 "sums="<<sums<< //the sum of shared clusters
4231 "xm0="<<xm[i0]<< // the center of track
4232 "xm1="<<xm[i1]<< // the x center of track
4233 // General cut variables
4234 "dfi="<<dfi<< // distance in fi angle
4235 "dtheta="<<dtheta<< // distance int theta angle
4241 "dist0="<<dist[0]<< //distance x
4242 "dist1="<<dist[1]<< //distance y
4243 "dist2="<<dist[2]<< //distance z
4244 "mdist0="<<mdist[0]<< //distance x
4245 "mdist1="<<mdist[1]<< //distance y
4246 "mdist2="<<mdist[2]<< //distance z
4260 if (AliTPCReconstructor::StreamLevel()>1) {
4261 AliInfo("Time for curling tracks removal DEBUGGING MC");
4268 void AliTPCtrackerMI::FindSplitted(TObjArray * array, AliESDEvent */*esd*/, Int_t /*iter*/){
4270 // Find Splitted tracks and remove the one with worst quality
4271 // Corresponding debug streamer to tune selections - "Splitted2"
4273 // 0. Sort tracks according quility
4274 // 1. Propagate the tracks to the reference radius
4275 // 2. Double_t loop to select close tracks (only to speed up process)
4276 // 3. Calculate cluster overlap ratio - and remove the track if bigger than a threshold
4277 // 4. Delete temporary parameters
4279 const Double_t xref=GetXrow(63); // reference radius -IROC/OROC boundary
4281 const Double_t kCutP1=10; // delta Z cut 10 cm
4282 const Double_t kCutP2=0.15; // delta snp(fi) cut 0.15
4283 const Double_t kCutP3=0.15; // delta tgl(theta) cut 0.15
4284 const Double_t kCutAlpha=0.15; // delta alpha cut
4285 Int_t firstpoint = 0;
4286 Int_t lastpoint = 160;
4288 Int_t nentries = array->GetEntriesFast();
4289 AliExternalTrackParam *params = new AliExternalTrackParam[nentries];
4295 //0. Sort tracks according quality
4296 //1. Propagate the ext. param to reference radius
4297 Int_t nseed = array->GetEntriesFast();
4298 Float_t * quality = new Float_t[nseed];
4299 Int_t * indexes = new Int_t[nseed];
4300 for (Int_t i=0; i<nseed; i++) {
4301 AliTPCseed *pt=(AliTPCseed*)array->UncheckedAt(i);
4306 pt->UpdatePoints(); //select first last max dens points
4307 Float_t * points = pt->GetPoints();
4308 if (points[3]<0.8) quality[i] =-1;
4309 quality[i] = (points[2]-points[0])+pt->GetNumberOfClusters();
4310 //prefer high momenta tracks if overlaps
4311 quality[i] *= TMath::Sqrt(TMath::Abs(pt->Pt())+0.5);
4313 AliTracker::PropagateTrackToBxByBz(&(params[i]),xref,pt->GetMass(),5.,kTRUE);
4314 AliTracker::PropagateTrackToBxByBz(&(params[i]),xref,pt->GetMass(),1.,kTRUE);
4316 TMath::Sort(nseed,quality,indexes);
4318 // 3. Loop over pair of tracks
4320 for (Int_t i0=0; i0<nseed; i0++) {
4321 Int_t index0=indexes[i0];
4322 if (!(array->UncheckedAt(index0))) continue;
4323 AliTPCseed *s1 = (AliTPCseed*)array->UncheckedAt(index0);
4324 if (!s1->IsActive()) continue;
4325 AliExternalTrackParam &par0=params[index0];
4326 for (Int_t i1=i0+1; i1<nseed; i1++) {
4327 Int_t index1=indexes[i1];
4328 if (!(array->UncheckedAt(index1))) continue;
4329 AliTPCseed *s2 = (AliTPCseed*)array->UncheckedAt(index1);
4330 if (!s2->IsActive()) continue;
4331 if (s2->GetKinkIndexes()[0]!=0)
4332 if (s2->GetKinkIndexes()[0] == -s1->GetKinkIndexes()[0]) continue;
4333 AliExternalTrackParam &par1=params[index1];
4334 if (TMath::Abs(par0.GetParameter()[3]-par1.GetParameter()[3])>kCutP3) continue;
4335 if (TMath::Abs(par0.GetParameter()[1]-par1.GetParameter()[1])>kCutP1) continue;
4336 if (TMath::Abs(par0.GetParameter()[2]-par1.GetParameter()[2])>kCutP2) continue;
4337 Double_t dAlpha= TMath::Abs(par0.GetAlpha()-par1.GetAlpha());
4338 if (dAlpha>TMath::Pi()) dAlpha-=TMath::Pi();
4339 if (TMath::Abs(dAlpha)>kCutAlpha) continue;
4344 Int_t firstShared=lastpoint, lastShared=firstpoint;
4345 Int_t firstRow=lastpoint, lastRow=firstpoint;
4347 for (Int_t i=firstpoint;i<lastpoint;i++){
4348 if (s1->GetClusterIndex2(i)>0) nall0++;
4349 if (s2->GetClusterIndex2(i)>0) nall1++;
4350 if (s1->GetClusterIndex2(i)>0 && s2->GetClusterIndex2(i)>0) {
4351 if (i<firstRow) firstRow=i;
4352 if (i>lastRow) lastRow=i;
4354 if ( (s1->GetClusterIndex2(i))==(s2->GetClusterIndex2(i)) && s1->GetClusterIndex2(i)>0) {
4355 if (i<firstShared) firstShared=i;
4356 if (i>lastShared) lastShared=i;
4360 Double_t ratio0 = Float_t(sumShared)/Float_t(TMath::Min(nall0+1,nall1+1));
4361 Double_t ratio1 = Float_t(sumShared)/Float_t(TMath::Max(nall0+1,nall1+1));
4363 if( AliTPCReconstructor::StreamLevel()>1){
4364 TTreeSRedirector &cstream = *fDebugStreamer;
4365 Int_t n0=s1->GetNumberOfClusters();
4366 Int_t n1=s2->GetNumberOfClusters();
4367 Int_t n0F=s1->GetNFoundable();
4368 Int_t n1F=s2->GetNFoundable();
4369 Int_t lab0=s1->GetLabel();
4370 Int_t lab1=s2->GetLabel();
4372 cstream<<"Splitted2"<<
4373 "iter="<<fIteration<<
4374 "lab0="<<lab0<< // MC label if exist
4375 "lab1="<<lab1<< // MC label if exist
4378 "ratio0="<<ratio0<< // shared ratio
4379 "ratio1="<<ratio1<< // shared ratio
4380 "p0.="<<&par0<< // track parameters
4382 "s0.="<<s1<< // full seed
4384 "n0="<<n0<< // number of clusters track 0
4385 "n1="<<n1<< // number of clusters track 1
4386 "nall0="<<nall0<< // number of clusters track 0
4387 "nall1="<<nall1<< // number of clusters track 1
4388 "n0F="<<n0F<< // number of findable
4389 "n1F="<<n1F<< // number of findable
4390 "shared="<<sumShared<< // number of shared clusters
4391 "firstS="<<firstShared<< // first and the last shared row
4392 "lastS="<<lastShared<<
4393 "firstRow="<<firstRow<< // first and the last row with cluster
4394 "lastRow="<<lastRow<< //
4398 // remove track with lower quality
4400 if (ratio0>AliTPCReconstructor::GetRecoParam()->GetCutSharedClusters(0) ||
4401 ratio1>AliTPCReconstructor::GetRecoParam()->GetCutSharedClusters(1)){
4405 delete array->RemoveAt(index1);
4410 // 4. Delete temporary array
4417 void AliTPCtrackerMI::FindCurling(const TObjArray * array, AliESDEvent */*esd*/, Int_t iter)
4420 // find Curling tracks
4421 // Use AliTPCReconstructor::StreamLevel()>1 if you want to tune parameters - cuts
4424 // Algorithm done in 2 phases - because of CPU consumption
4425 // it is n^2 algorithm - for lead-lead 20000x20000 combination are investigated
4426 // see detal in MC part what can be used to cut
4430 const Float_t kMaxC = 400; // maximal curvature to of the track
4431 const Float_t kMaxdTheta = 0.15; // maximal distance in theta
4432 const Float_t kMaxdPhi = 0.15; // maximal distance in phi
4433 const Float_t kPtRatio = 0.3; // ratio between pt
4434 const Float_t kMinDCAR = 2.; // distance to the primary vertex in r - see cpipe cut
4437 // Curling tracks cuts
4440 const Float_t kMaxDeltaRMax = 40; // distance in outer radius
4441 const Float_t kMaxDeltaRMin = 5.; // distance in lower radius - see cpipe cut
4442 const Float_t kMinAngle = 2.9; // angle between tracks
4443 const Float_t kMaxDist = 5; // biggest distance
4445 // The cuts can be tuned using the "MC information stored in Multi tree ==> see FindMultiMC
4448 TCut csign("csign","Tr0.fP[4]*Tr1.fP[4]<0"); //opposite sign
4449 TCut cmax("cmax","abs(Tr0.GetC())>1/400");
4450 TCut cda("cda","sqrt(dtheta^2+dfi^2)<0.15");
4451 TCut ccratio("ccratio","abs((Tr0.fP[4]+Tr1.fP[4])/(abs(Tr0.fP[4])+abs(Tr1.fP[4])))<0.3");
4452 TCut cpipe("cpipe", "min(abs(r0-rc0),abs(r1-rc1))>5");
4454 TCut cdrmax("cdrmax","abs(abs(rc0+r0)-abs(rc1+r1))<40")
4455 TCut cdrmin("cdrmin","abs(abs(rc0+r0)-abs(rc1+r1))<10")
4457 Multi->Draw("dfi","iter==0"+csign+cmax+cda+ccratio); ~94% of curling tracks fulfill
4458 Multi->Draw("min(abs(r0-rc0),abs(r1-rc1))","iter==0&&abs(lab1)==abs(lab0)"+csign+cmax+cda+ccratio+cpipe+cdrmin+cdrmax); //80%
4460 Curling2->Draw("dfi","iter==0&&abs(lab0)==abs(lab1)"+csign+cmax+cdtheta+cdfi+ccratio)
4466 Int_t nentries = array->GetEntriesFast();
4467 AliHelix *helixes = new AliHelix[nentries];
4468 for (Int_t i=0;i<nentries;i++){
4469 AliTPCseed* track = (AliTPCseed*)array->At(i);
4470 if (!track) continue;
4471 track->SetCircular(0);
4472 new (&helixes[i]) AliHelix(*track);
4478 Double_t phase[2][2],radius[2];
4483 for (Int_t i0=0;i0<nentries;i0++){
4484 AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
4485 if (!track0) continue;
4486 if (TMath::Abs(track0->GetC())<1/kMaxC) continue;
4487 Float_t xc0 = helixes[i0].GetHelix(6);
4488 Float_t yc0 = helixes[i0].GetHelix(7);
4489 Float_t r0 = helixes[i0].GetHelix(8);
4490 Float_t rc0 = TMath::Sqrt(xc0*xc0+yc0*yc0);
4491 Float_t fi0 = TMath::ATan2(yc0,xc0);
4493 for (Int_t i1=i0+1;i1<nentries;i1++){
4494 AliTPCseed * track1 = (AliTPCseed*)array->At(i1);
4495 if (!track1) continue;
4496 if (TMath::Abs(track1->GetC())<1/kMaxC) continue;
4497 Float_t xc1 = helixes[i1].GetHelix(6);
4498 Float_t yc1 = helixes[i1].GetHelix(7);
4499 Float_t r1 = helixes[i1].GetHelix(8);
4500 Float_t rc1 = TMath::Sqrt(xc1*xc1+yc1*yc1);
4501 Float_t fi1 = TMath::ATan2(yc1,xc1);
4503 Float_t dfi = fi0-fi1;
4506 if (dfi>1.5*TMath::Pi()) dfi-=TMath::Pi(); // take care about edge effect
4507 if (dfi<-1.5*TMath::Pi()) dfi+=TMath::Pi(); //
4508 Float_t dtheta = TMath::Abs(track0->GetTgl()-track1->GetTgl())<TMath::Abs(track0->GetTgl()+track1->GetTgl())? track0->GetTgl()-track1->GetTgl():track0->GetTgl()+track1->GetTgl();
4512 if (track0->GetBConstrain()&&track1->GetBConstrain()) continue; // not constrained
4513 if (track1->GetSigned1Pt()*track0->GetSigned1Pt()>0) continue; // not the same sign
4514 if ( TMath::Abs(track1->GetTgl()+track0->GetTgl())>kMaxdTheta) continue; //distance in the Theta
4515 if ( TMath::Abs(dfi)>kMaxdPhi) continue; //distance in phi
4516 if ( TMath::Sqrt(dfi*dfi+dtheta*dtheta)>kMaxdPhi) continue; //common angular offset
4518 Float_t pt0 = track0->GetSignedPt();
4519 Float_t pt1 = track1->GetSignedPt();
4520 if ((TMath::Abs(pt0+pt1)/(TMath::Abs(pt0)+TMath::Abs(pt1)))>kPtRatio) continue;
4521 if ((iter==1) && TMath::Abs(TMath::Abs(rc0+r0)-TMath::Abs(rc1+r1))>kMaxDeltaRMax) continue;
4522 if ((iter!=1) &&TMath::Abs(TMath::Abs(rc0-r0)-TMath::Abs(rc1-r1))>kMaxDeltaRMin) continue;
4523 if (TMath::Min(TMath::Abs(rc0-r0),TMath::Abs(rc1-r1))<kMinDCAR) continue;
4526 // Now find closest approach
4530 Int_t npoints = helixes[i0].GetRPHIintersections(helixes[i1], phase, radius,10);
4531 if (npoints==0) continue;
4532 helixes[i0].GetClosestPhases(helixes[i1], phase);
4536 Double_t hangles[3];
4537 helixes[i0].Evaluate(phase[0][0],xyz0);
4538 helixes[i1].Evaluate(phase[0][1],xyz1);
4540 helixes[i0].GetAngle(phase[0][0],helixes[i1],phase[0][1],hangles);
4541 Double_t deltah[2],deltabest;
4542 if (TMath::Abs(hangles[2])<kMinAngle) continue;
4546 helixes[i0].ParabolicDCA(helixes[i1],phase[0][0],phase[0][1],radius[0],deltah[0],2);
4548 helixes[i0].ParabolicDCA(helixes[i1],phase[1][0],phase[1][1],radius[1],deltah[1],2);
4549 if (deltah[1]<deltah[0]) ibest=1;
4551 deltabest = TMath::Sqrt(deltah[ibest]);
4552 helixes[i0].Evaluate(phase[ibest][0],xyz0);
4553 helixes[i1].Evaluate(phase[ibest][1],xyz1);
4554 helixes[i0].GetAngle(phase[ibest][0],helixes[i1],phase[ibest][1],hangles);
4555 Double_t radiusbest = TMath::Sqrt(radius[ibest]);
4557 if (deltabest>kMaxDist) continue;
4558 // if (mindcar+mindcaz<40 && (TMath::Abs(hangles[2])<kMinAngle ||deltabest>3)) continue;
4559 Bool_t sign =kFALSE;
4560 if (hangles[2]>kMinAngle) sign =kTRUE;
4563 // circular[i0] = kTRUE;
4564 // circular[i1] = kTRUE;
4565 if (track0->OneOverPt()<track1->OneOverPt()){
4566 track0->SetCircular(track0->GetCircular()+1);
4567 track1->SetCircular(track1->GetCircular()+2);
4570 track1->SetCircular(track1->GetCircular()+1);
4571 track0->SetCircular(track0->GetCircular()+2);
4574 if (AliTPCReconstructor::StreamLevel()>2){
4576 //debug stream to tune "fine" cuts
4577 Int_t lab0=track0->GetLabel();
4578 Int_t lab1=track1->GetLabel();
4579 TTreeSRedirector &cstream = *fDebugStreamer;
4580 cstream<<"Curling2"<<
4596 "npoints="<<npoints<<
4597 "hangles0="<<hangles[0]<<
4598 "hangles1="<<hangles[1]<<
4599 "hangles2="<<hangles[2]<<
4602 "radius="<<radiusbest<<
4603 "deltabest="<<deltabest<<
4604 "phase0="<<phase[ibest][0]<<
4605 "phase1="<<phase[ibest][1]<<
4613 if (AliTPCReconstructor::StreamLevel()>1) {
4614 AliInfo("Time for curling tracks removal");
4623 void AliTPCtrackerMI::FindKinks(TObjArray * array, AliESDEvent *esd)
4630 TObjArray *kinks= new TObjArray(10000);
4631 // TObjArray *v0s= new TObjArray(10000);
4632 Int_t nentries = array->GetEntriesFast();
4633 AliHelix *helixes = new AliHelix[nentries];
4634 Int_t *sign = new Int_t[nentries];
4635 Int_t *nclusters = new Int_t[nentries];
4636 Float_t *alpha = new Float_t[nentries];
4637 AliKink *kink = new AliKink();
4638 Int_t * usage = new Int_t[nentries];
4639 Float_t *zm = new Float_t[nentries];
4640 Float_t *z0 = new Float_t[nentries];
4641 Float_t *fim = new Float_t[nentries];
4642 Float_t *shared = new Float_t[nentries];
4643 Bool_t *circular = new Bool_t[nentries];
4644 Float_t *dca = new Float_t[nentries];
4645 //const AliESDVertex * primvertex = esd->GetVertex();
4647 // nentries = array->GetEntriesFast();
4652 for (Int_t i=0;i<nentries;i++){
4655 AliTPCseed* track = (AliTPCseed*)array->At(i);
4656 if (!track) continue;
4657 track->SetCircular(0);
4659 track->UpdatePoints();
4660 if (( track->GetPoints()[2]- track->GetPoints()[0])>5 && track->GetPoints()[3]>0.8){
4662 nclusters[i]=track->GetNumberOfClusters();
4663 alpha[i] = track->GetAlpha();
4664 new (&helixes[i]) AliHelix(*track);
4666 helixes[i].Evaluate(0,xyz);
4667 sign[i] = (track->GetC()>0) ? -1:1;
4670 if (track->GetProlongation(x,y,z)){
4672 fim[i] = alpha[i]+TMath::ATan2(y,x);
4675 zm[i] = track->GetZ();
4679 circular[i]= kFALSE;
4680 if (track->GetProlongation(0,y,z)) z0[i] = z;
4681 dca[i] = track->GetD(0,0);
4687 Int_t ncandidates =0;
4690 Double_t phase[2][2],radius[2];
4693 // Find circling track
4695 for (Int_t i0=0;i0<nentries;i0++){
4696 AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
4697 if (!track0) continue;
4698 if (track0->GetNumberOfClusters()<40) continue;
4699 if (TMath::Abs(1./track0->GetC())>200) continue;
4700 for (Int_t i1=i0+1;i1<nentries;i1++){
4701 AliTPCseed * track1 = (AliTPCseed*)array->At(i1);
4702 if (!track1) continue;
4703 if (track1->GetNumberOfClusters()<40) continue;
4704 if ( TMath::Abs(track1->GetTgl()+track0->GetTgl())>0.1) continue;
4705 if (track0->GetBConstrain()&&track1->GetBConstrain()) continue;
4706 if (TMath::Abs(1./track1->GetC())>200) continue;
4707 if (track1->GetSigned1Pt()*track0->GetSigned1Pt()>0) continue;
4708 if (track1->GetTgl()*track0->GetTgl()>0) continue;
4709 if (TMath::Max(TMath::Abs(1./track0->GetC()),TMath::Abs(1./track1->GetC()))>190) continue;
4710 if (track0->GetBConstrain()&&track1->OneOverPt()<track0->OneOverPt()) continue; //returning - lower momenta
4711 if (track1->GetBConstrain()&&track0->OneOverPt()<track1->OneOverPt()) continue; //returning - lower momenta
4713 Float_t mindcar = TMath::Min(TMath::Abs(dca[i0]),TMath::Abs(dca[i1]));
4714 if (mindcar<5) continue;
4715 Float_t mindcaz = TMath::Min(TMath::Abs(z0[i0]-GetZ()),TMath::Abs(z0[i1]-GetZ()));
4716 if (mindcaz<5) continue;
4717 if (mindcar+mindcaz<20) continue;
4720 Float_t xc0 = helixes[i0].GetHelix(6);
4721 Float_t yc0 = helixes[i0].GetHelix(7);
4722 Float_t r0 = helixes[i0].GetHelix(8);
4723 Float_t xc1 = helixes[i1].GetHelix(6);
4724 Float_t yc1 = helixes[i1].GetHelix(7);
4725 Float_t r1 = helixes[i1].GetHelix(8);
4727 Float_t rmean = (r0+r1)*0.5;
4728 Float_t delta =TMath::Sqrt((xc1-xc0)*(xc1-xc0)+(yc1-yc0)*(yc1-yc0));
4729 //if (delta>30) continue;
4730 if (delta>rmean*0.25) continue;
4731 if (TMath::Abs(r0-r1)/rmean>0.3) continue;
4733 Int_t npoints = helixes[i0].GetRPHIintersections(helixes[i1], phase, radius,10);
4734 if (npoints==0) continue;
4735 helixes[i0].GetClosestPhases(helixes[i1], phase);
4739 Double_t hangles[3];
4740 helixes[i0].Evaluate(phase[0][0],xyz0);
4741 helixes[i1].Evaluate(phase[0][1],xyz1);
4743 helixes[i0].GetAngle(phase[0][0],helixes[i1],phase[0][1],hangles);
4744 Double_t deltah[2],deltabest;
4745 if (hangles[2]<2.8) continue;
4748 helixes[i0].ParabolicDCA(helixes[i1],phase[0][0],phase[0][1],radius[0],deltah[0],2);
4750 helixes[i0].ParabolicDCA(helixes[i1],phase[1][0],phase[1][1],radius[1],deltah[1],2);
4751 if (deltah[1]<deltah[0]) ibest=1;
4753 deltabest = TMath::Sqrt(deltah[ibest]);
4754 helixes[i0].Evaluate(phase[ibest][0],xyz0);
4755 helixes[i1].Evaluate(phase[ibest][1],xyz1);
4756 helixes[i0].GetAngle(phase[ibest][0],helixes[i1],phase[ibest][1],hangles);
4757 Double_t radiusbest = TMath::Sqrt(radius[ibest]);
4759 if (deltabest>6) continue;
4760 if (mindcar+mindcaz<40 && (hangles[2]<3.12||deltabest>3)) continue;
4761 Bool_t lsign =kFALSE;
4762 if (hangles[2]>3.06) lsign =kTRUE;
4765 circular[i0] = kTRUE;
4766 circular[i1] = kTRUE;
4767 if (track0->OneOverPt()<track1->OneOverPt()){
4768 track0->SetCircular(track0->GetCircular()+1);
4769 track1->SetCircular(track1->GetCircular()+2);
4772 track1->SetCircular(track1->GetCircular()+1);
4773 track0->SetCircular(track0->GetCircular()+2);
4776 if (lsign&&AliTPCReconstructor::StreamLevel()>1){
4778 Int_t lab0=track0->GetLabel();
4779 Int_t lab1=track1->GetLabel();
4780 TTreeSRedirector &cstream = *fDebugStreamer;
4781 cstream<<"Curling"<<
4788 "mindcar="<<mindcar<<
4789 "mindcaz="<<mindcaz<<
4792 "npoints="<<npoints<<
4793 "hangles0="<<hangles[0]<<
4794 "hangles2="<<hangles[2]<<
4799 "radius="<<radiusbest<<
4800 "deltabest="<<deltabest<<
4801 "phase0="<<phase[ibest][0]<<
4802 "phase1="<<phase[ibest][1]<<
4812 for (Int_t i =0;i<nentries;i++){
4813 if (sign[i]==0) continue;
4814 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
4821 Double_t cradius0 = 40*40;
4822 Double_t cradius1 = 270*270;
4825 Double_t cdist3=0.55;
4826 for (Int_t j =i+1;j<nentries;j++){
4828 if (sign[j]*sign[i]<1) continue;
4829 if ( (nclusters[i]+nclusters[j])>200) continue;
4830 if ( (nclusters[i]+nclusters[j])<80) continue;
4831 if ( TMath::Abs(zm[i]-zm[j])>60.) continue;
4832 if ( TMath::Abs(fim[i]-fim[j])>0.6 && TMath::Abs(fim[i]-fim[j])<5.7 ) continue;
4833 //AliTPCseed * track1 = (AliTPCseed*)array->At(j); Double_t phase[2][2],radius[2];
4834 Int_t npoints = helixes[i].GetRPHIintersections(helixes[j], phase, radius,20);
4835 if (npoints<1) continue;
4838 if (radius[0]<cradius0||radius[0]>cradius1) continue;
4841 if ( (radius[0]<cradius0||radius[0]>cradius1) && (radius[1]<cradius0||radius[1]>cradius1) ) continue;
4844 Double_t delta1=10000,delta2=10000;
4845 // cuts on the intersection radius
4846 helixes[i].LinearDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta1);
4847 if (radius[0]<20&&delta1<1) continue; //intersection at vertex
4848 if (radius[0]<10&&delta1<3) continue; //intersection at vertex
4850 helixes[i].LinearDCA(helixes[j],phase[1][0],phase[1][1],radius[1],delta2);
4851 if (radius[1]<20&&delta2<1) continue; //intersection at vertex
4852 if (radius[1]<10&&delta2<3) continue; //intersection at vertex
4855 Double_t distance1 = TMath::Min(delta1,delta2);
4856 if (distance1>cdist1) continue; // cut on DCA linear approximation
4858 npoints = helixes[i].GetRPHIintersections(helixes[j], phase, radius,20);
4859 helixes[i].ParabolicDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta1);
4860 if (radius[0]<20&&delta1<1) continue; //intersection at vertex
4861 if (radius[0]<10&&delta1<3) continue; //intersection at vertex
4864 helixes[i].ParabolicDCA(helixes[j],phase[1][0],phase[1][1],radius[1],delta2);
4865 if (radius[1]<20&&delta2<1) continue; //intersection at vertex
4866 if (radius[1]<10&&delta2<3) continue; //intersection at vertex
4868 distance1 = TMath::Min(delta1,delta2);
4871 rkink = TMath::Sqrt(radius[0]);
4874 rkink = TMath::Sqrt(radius[1]);
4876 if (distance1>cdist2) continue;
4879 AliTPCseed * track1 = (AliTPCseed*)array->At(j);
4882 Int_t row0 = GetRowNumber(rkink);
4883 if (row0<10) continue;
4884 if (row0>150) continue;
4887 Float_t dens00=-1,dens01=-1;
4888 Float_t dens10=-1,dens11=-1;
4890 Int_t found,foundable,ishared;
4891 track0->GetClusterStatistic(0,row0-5, found, foundable,ishared,kFALSE);
4892 if (foundable>5) dens00 = Float_t(found)/Float_t(foundable);
4893 track0->GetClusterStatistic(row0+5,155, found, foundable,ishared,kFALSE);
4894 if (foundable>5) dens01 = Float_t(found)/Float_t(foundable);
4896 track1->GetClusterStatistic(0,row0-5, found, foundable,ishared,kFALSE);
4897 if (foundable>10) dens10 = Float_t(found)/Float_t(foundable);
4898 track1->GetClusterStatistic(row0+5,155, found, foundable,ishared,kFALSE);
4899 if (foundable>10) dens11 = Float_t(found)/Float_t(foundable);
4901 if (dens00<dens10 && dens01<dens11) continue;
4902 if (dens00>dens10 && dens01>dens11) continue;
4903 if (TMath::Max(dens00,dens10)<0.1) continue;
4904 if (TMath::Max(dens01,dens11)<0.3) continue;
4906 if (TMath::Min(dens00,dens10)>0.6) continue;
4907 if (TMath::Min(dens01,dens11)>0.6) continue;
4910 AliTPCseed * ktrack0, *ktrack1;
4919 if (TMath::Abs(ktrack0->GetC())>5) continue; // cut on the curvature for mother particle
4920 AliExternalTrackParam paramm(*ktrack0);
4921 AliExternalTrackParam paramd(*ktrack1);
4922 if (row0>60&&ktrack1->GetReference().GetX()>90.)new (¶md) AliExternalTrackParam(ktrack1->GetReference());
4925 kink->SetMother(paramm);
4926 kink->SetDaughter(paramd);
4929 Float_t x[3] = { kink->GetPosition()[0],kink->GetPosition()[1],kink->GetPosition()[2]};
4931 fkParam->Transform0to1(x,index);
4932 fkParam->Transform1to2(x,index);
4933 row0 = GetRowNumber(x[0]);
4935 if (kink->GetR()<100) continue;
4936 if (kink->GetR()>240) continue;
4937 if (kink->GetPosition()[2]/kink->GetR()>AliTPCReconstructor::GetCtgRange()) continue; //out of fiducial volume
4938 if (kink->GetDistance()>cdist3) continue;
4939 Float_t dird = kink->GetDaughterP()[0]*kink->GetPosition()[0]+kink->GetDaughterP()[1]*kink->GetPosition()[1]; // rough direction estimate
4940 if (dird<0) continue;
4942 Float_t dirm = kink->GetMotherP()[0]*kink->GetPosition()[0]+kink->GetMotherP()[1]*kink->GetPosition()[1]; // rough direction estimate
4943 if (dirm<0) continue;
4944 Float_t mpt = TMath::Sqrt(kink->GetMotherP()[0]*kink->GetMotherP()[0]+kink->GetMotherP()[1]*kink->GetMotherP()[1]);
4945 if (mpt<0.2) continue;
4948 //for high momenta momentum not defined well in first iteration
4949 Double_t qt = TMath::Sin(kink->GetAngle(2))*ktrack1->GetP();
4950 if (qt>0.35) continue;
4953 kink->SetLabel(CookLabel(ktrack0,0.4,0,row0),0);
4954 kink->SetLabel(CookLabel(ktrack1,0.4,row0,160),1);
4956 kink->SetTPCDensity(dens00,0,0);
4957 kink->SetTPCDensity(dens01,0,1);
4958 kink->SetTPCDensity(dens10,1,0);
4959 kink->SetTPCDensity(dens11,1,1);
4960 kink->SetIndex(i,0);
4961 kink->SetIndex(j,1);
4964 kink->SetTPCDensity(dens10,0,0);
4965 kink->SetTPCDensity(dens11,0,1);
4966 kink->SetTPCDensity(dens00,1,0);
4967 kink->SetTPCDensity(dens01,1,1);
4968 kink->SetIndex(j,0);
4969 kink->SetIndex(i,1);
4972 if (mpt<1||kink->GetAngle(2)>0.1){
4973 // angle and densities not defined yet
4974 if (kink->GetTPCDensityFactor()<0.8) continue;
4975 if ((2-kink->GetTPCDensityFactor())*kink->GetDistance() >0.25) continue;
4976 if (kink->GetAngle(2)*ktrack0->GetP()<0.003) continue; //too small angle
4977 if (kink->GetAngle(2)>0.2&&kink->GetTPCDensityFactor()<1.15) continue;
4978 if (kink->GetAngle(2)>0.2&&kink->GetTPCDensity(0,1)>0.05) continue;
4980 Float_t criticalangle = track0->GetSigmaSnp2()+track0->GetSigmaTgl2();
4981 criticalangle+= track1->GetSigmaSnp2()+track1->GetSigmaTgl2();
4982 criticalangle= 3*TMath::Sqrt(criticalangle);
4983 if (criticalangle>0.02) criticalangle=0.02;
4984 if (kink->GetAngle(2)<criticalangle) continue;
4987 Int_t drow = Int_t(2.+0.5/(0.05+kink->GetAngle(2))); // overlap region defined
4988 Float_t shapesum =0;
4990 for ( Int_t row = row0-drow; row<row0+drow;row++){
4991 if (row<0) continue;
4992 if (row>155) continue;
4993 if (ktrack0->GetClusterPointer(row)){
4994 AliTPCTrackerPoint *point =ktrack0->GetTrackPoint(row);
4995 shapesum+=point->GetSigmaY()+point->GetSigmaZ();
4998 if (ktrack1->GetClusterPointer(row)){
4999 AliTPCTrackerPoint *point =ktrack1->GetTrackPoint(row);
5000 shapesum+=point->GetSigmaY()+point->GetSigmaZ();
5005 kink->SetShapeFactor(-1.);
5008 kink->SetShapeFactor(shapesum/sum);
5010 // esd->AddKink(kink);
5012 // kink->SetMother(paramm);
5013 //kink->SetDaughter(paramd);
5015 Double_t chi2P2 = paramm.GetParameter()[2]-paramd.GetParameter()[2];
5017 chi2P2/=paramm.GetCovariance()[5]+paramd.GetCovariance()[5];
5018 Double_t chi2P3 = paramm.GetParameter()[3]-paramd.GetParameter()[3];
5020 chi2P3/=paramm.GetCovariance()[9]+paramd.GetCovariance()[9];
5022 if (AliTPCReconstructor::StreamLevel()>1) {
5023 (*fDebugStreamer)<<"kinkLpt"<<
5031 if ( chi2P2+chi2P3<AliTPCReconstructor::GetRecoParam()->GetKinkAngleCutChi2(0)){
5035 kinks->AddLast(kink);
5041 // sort the kinks according quality - and refit them towards vertex
5043 Int_t nkinks = kinks->GetEntriesFast();
5044 Float_t *quality = new Float_t[nkinks];
5045 Int_t *indexes = new Int_t[nkinks];
5046 AliTPCseed *mothers = new AliTPCseed[nkinks];
5047 AliTPCseed *daughters = new AliTPCseed[nkinks];
5050 for (Int_t i=0;i<nkinks;i++){
5052 AliKink *kinkl = (AliKink*)kinks->At(i);
5054 // refit kinks towards vertex
5056 Int_t index0 = kinkl->GetIndex(0);
5057 Int_t index1 = kinkl->GetIndex(1);
5058 AliTPCseed * ktrack0 = (AliTPCseed*)array->At(index0);
5059 AliTPCseed * ktrack1 = (AliTPCseed*)array->At(index1);
5061 Int_t sumn=ktrack0->GetNumberOfClusters()+ktrack1->GetNumberOfClusters();
5063 // Refit Kink under if too small angle
5065 if (kinkl->GetAngle(2)<0.05){
5066 kinkl->SetTPCRow0(GetRowNumber(kinkl->GetR()));
5067 Int_t row0 = kinkl->GetTPCRow0();
5068 Int_t drow = Int_t(2.+0.5/(0.05+kinkl->GetAngle(2)));
5071 Int_t last = row0-drow;
5072 if (last<40) last=40;
5073 if (last<ktrack0->GetFirstPoint()+25) last = ktrack0->GetFirstPoint()+25;
5074 AliTPCseed* seed0 = ReSeed(ktrack0,last,kFALSE);
5077 Int_t first = row0+drow;
5078 if (first>130) first=130;
5079 if (first>ktrack1->GetLastPoint()-25) first = TMath::Max(ktrack1->GetLastPoint()-25,30);
5080 AliTPCseed* seed1 = ReSeed(ktrack1,first,kTRUE);
5082 if (seed0 && seed1){
5083 kinkl->SetStatus(1,8);
5084 if (RefitKink(*seed0,*seed1,*kinkl)) kinkl->SetStatus(1,9);
5085 row0 = GetRowNumber(kinkl->GetR());
5086 sumn = seed0->GetNumberOfClusters()+seed1->GetNumberOfClusters();
5087 mothers[i] = *seed0;
5088 daughters[i] = *seed1;
5091 delete kinks->RemoveAt(i);
5092 if (seed0) delete seed0;
5093 if (seed1) delete seed1;
5096 if (kinkl->GetDistance()>0.5 || kinkl->GetR()<110 || kinkl->GetR()>240) {
5097 delete kinks->RemoveAt(i);
5098 if (seed0) delete seed0;
5099 if (seed1) delete seed1;
5107 if (kinkl) quality[i] = 160*((0.1+kinkl->GetDistance())*(2.-kinkl->GetTPCDensityFactor()))/(sumn+40.); //the longest -clossest will win
5109 TMath::Sort(nkinks,quality,indexes,kFALSE);
5111 //remove double find kinks
5113 for (Int_t ikink0=1;ikink0<nkinks;ikink0++){
5114 AliKink * kink0 = (AliKink*) kinks->At(indexes[ikink0]);
5115 if (!kink0) continue;
5117 for (Int_t ikink1=0;ikink1<ikink0;ikink1++){
5118 if (!kink0) continue;
5119 AliKink * kink1 = (AliKink*) kinks->At(indexes[ikink1]);
5120 if (!kink1) continue;
5121 // if not close kink continue
5122 if (TMath::Abs(kink1->GetPosition()[2]-kink0->GetPosition()[2])>10) continue;
5123 if (TMath::Abs(kink1->GetPosition()[1]-kink0->GetPosition()[1])>10) continue;
5124 if (TMath::Abs(kink1->GetPosition()[0]-kink0->GetPosition()[0])>10) continue;
5126 AliTPCseed &mother0 = mothers[indexes[ikink0]];
5127 AliTPCseed &daughter0 = daughters[indexes[ikink0]];
5128 AliTPCseed &mother1 = mothers[indexes[ikink1]];
5129 AliTPCseed &daughter1 = daughters[indexes[ikink1]];
5130 Int_t row0 = (kink0->GetTPCRow0()+kink1->GetTPCRow0())/2;
5139 for (Int_t i=0;i<row0;i++){
5140 if (mother0.GetClusterIndex(i)>0 && mother1.GetClusterIndex(i)>0){
5143 if (mother0.GetClusterIndex(i)==mother1.GetClusterIndex(i)){
5150 for (Int_t i=row0;i<158;i++){
5151 if (daughter0.GetClusterIndex(i)>0 && daughter0.GetClusterIndex(i)>0){
5154 if (mother0.GetClusterIndex(i)==mother1.GetClusterIndex(i)){
5160 Float_t ratio = Float_t(same+1)/Float_t(both+1);
5161 Float_t ratiom = Float_t(samem+1)/Float_t(bothm+1);
5162 Float_t ratiod = Float_t(samed+1)/Float_t(bothd+1);
5163 if (ratio>0.3 && ratiom>0.5 &&ratiod>0.5) {
5164 Int_t sum0 = mother0.GetNumberOfClusters()+daughter0.GetNumberOfClusters();
5165 Int_t sum1 = mother1.GetNumberOfClusters()+daughter1.GetNumberOfClusters();
5167 shared[kink0->GetIndex(0)]= kTRUE;
5168 shared[kink0->GetIndex(1)]= kTRUE;
5169 delete kinks->RemoveAt(indexes[ikink0]);
5172 shared[kink1->GetIndex(0)]= kTRUE;
5173 shared[kink1->GetIndex(1)]= kTRUE;
5174 delete kinks->RemoveAt(indexes[ikink1]);
5181 for (Int_t i=0;i<nkinks;i++){
5182 AliKink * kinkl = (AliKink*) kinks->At(indexes[i]);
5183 if (!kinkl) continue;
5184 kinkl->SetTPCRow0(GetRowNumber(kinkl->GetR()));
5185 Int_t index0 = kinkl->GetIndex(0);
5186 Int_t index1 = kinkl->GetIndex(1);
5187 if (circular[index0]||(circular[index1]&&kinkl->GetDistance()>0.2)) continue;
5188 kinkl->SetMultiple(usage[index0],0);
5189 kinkl->SetMultiple(usage[index1],1);
5190 if (kinkl->GetMultiple()[0]+kinkl->GetMultiple()[1]>2) continue;
5191 if (kinkl->GetMultiple()[0]+kinkl->GetMultiple()[1]>0 && quality[indexes[i]]>0.2) continue;
5192 if (kinkl->GetMultiple()[0]+kinkl->GetMultiple()[1]>0 && kinkl->GetDistance()>0.2) continue;
5193 if (circular[index0]||(circular[index1]&&kinkl->GetDistance()>0.1)) continue;
5195 AliTPCseed * ktrack0 = (AliTPCseed*)array->At(index0);
5196 AliTPCseed * ktrack1 = (AliTPCseed*)array->At(index1);
5197 if (!ktrack0 || !ktrack1) continue;
5198 Int_t index = esd->AddKink(kinkl);
5201 if ( ktrack0->GetKinkIndex(0)==0 && ktrack1->GetKinkIndex(0)==0) { //best kink
5202 if (mothers[indexes[i]].GetNumberOfClusters()>20 && daughters[indexes[i]].GetNumberOfClusters()>20 && (mothers[indexes[i]].GetNumberOfClusters()+daughters[indexes[i]].GetNumberOfClusters())>100){
5203 *ktrack0 = mothers[indexes[i]];
5204 *ktrack1 = daughters[indexes[i]];
5208 ktrack0->SetKinkIndex(usage[index0],-(index+1));
5209 ktrack1->SetKinkIndex(usage[index1], (index+1));
5214 // Remove tracks corresponding to shared kink's
5216 for (Int_t i=0;i<nentries;i++){
5217 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5218 if (!track0) continue;
5219 if (track0->GetKinkIndex(0)!=0) continue;
5220 if (shared[i]) delete array->RemoveAt(i);
5225 RemoveUsed2(array,0.5,0.4,30);
5227 for (Int_t i=0;i<nentries;i++){
5228 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5229 if (!track0) continue;
5230 track0->CookdEdx(0.02,0.6);
5234 for (Int_t i=0;i<nentries;i++){
5235 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5236 if (!track0) continue;
5237 if (track0->Pt()<1.4) continue;
5238 //remove double high momenta tracks - overlapped with kink candidates
5241 for (Int_t icl=track0->GetFirstPoint();icl<track0->GetLastPoint(); icl++){
5242 if (track0->GetClusterPointer(icl)!=0){
5244 if (track0->GetClusterPointer(icl)->IsUsed(10)) ishared++;
5247 if (Float_t(ishared+1)/Float_t(all+1)>0.5) {
5248 delete array->RemoveAt(i);
5252 if (track0->GetKinkIndex(0)!=0) continue;
5253 if (track0->GetNumberOfClusters()<80) continue;
5255 AliTPCseed *pmother = new AliTPCseed();
5256 AliTPCseed *pdaughter = new AliTPCseed();
5257 AliKink *pkink = new AliKink;
5259 AliTPCseed & mother = *pmother;
5260 AliTPCseed & daughter = *pdaughter;
5261 AliKink & kinkl = *pkink;
5262 if (CheckKinkPoint(track0,mother,daughter, kinkl)){
5263 if (mother.GetNumberOfClusters()<30||daughter.GetNumberOfClusters()<20) {
5267 continue; //too short tracks
5269 if (mother.Pt()<1.4) {
5275 Int_t row0= kinkl.GetTPCRow0();
5276 if (kinkl.GetDistance()>0.5 || kinkl.GetR()<110. || kinkl.GetR()>240.) {
5283 Int_t index = esd->AddKink(&kinkl);
5284 mother.SetKinkIndex(0,-(index+1));
5285 daughter.SetKinkIndex(0,index+1);
5286 if (mother.GetNumberOfClusters()>50) {
5287 delete array->RemoveAt(i);
5288 array->AddAt(new AliTPCseed(mother),i);
5291 array->AddLast(new AliTPCseed(mother));
5293 array->AddLast(new AliTPCseed(daughter));
5294 for (Int_t icl=0;icl<row0;icl++) {
5295 if (mother.GetClusterPointer(icl)) mother.GetClusterPointer(icl)->Use(20);
5298 for (Int_t icl=row0;icl<158;icl++) {
5299 if (daughter.GetClusterPointer(icl)) daughter.GetClusterPointer(icl)->Use(20);
5308 delete [] daughters;
5330 printf("Ncandidates=\t%d\t%d\t%d\t%d\n",esd->GetNumberOfKinks(),ncandidates,ntracks,nall);
5334 void AliTPCtrackerMI::FindV0s(const TObjArray * array, AliESDEvent *const esd)
5340 TObjArray *tpcv0s = new TObjArray(100000);
5341 Int_t nentries = array->GetEntriesFast();
5342 AliHelix *helixes = new AliHelix[nentries];
5343 Int_t *sign = new Int_t[nentries];
5344 Float_t *alpha = new Float_t[nentries];
5345 Float_t *z0 = new Float_t[nentries];
5346 Float_t *dca = new Float_t[nentries];
5347 Float_t *sdcar = new Float_t[nentries];
5348 Float_t *cdcar = new Float_t[nentries];
5349 Float_t *pulldcar = new Float_t[nentries];
5350 Float_t *pulldcaz = new Float_t[nentries];
5351 Float_t *pulldca = new Float_t[nentries];
5352 Bool_t *isPrim = new Bool_t[nentries];
5353 const AliESDVertex * primvertex = esd->GetVertex();
5354 Double_t zvertex = primvertex->GetZv();
5356 // nentries = array->GetEntriesFast();
5358 for (Int_t i=0;i<nentries;i++){
5361 AliTPCseed* track = (AliTPCseed*)array->At(i);
5362 if (!track) continue;
5363 track->GetV0Indexes()[0] = 0; //rest v0 indexes
5364 track->GetV0Indexes()[1] = 0; //rest v0 indexes
5365 track->GetV0Indexes()[2] = 0; //rest v0 indexes
5367 alpha[i] = track->GetAlpha();
5368 new (&helixes[i]) AliHelix(*track);
5370 helixes[i].Evaluate(0,xyz);
5371 sign[i] = (track->GetC()>0) ? -1:1;
5375 if (track->GetProlongation(0,y,z)) z0[i] = z;
5376 dca[i] = track->GetD(0,0);
5378 // dca error parrameterezation + pulls
5380 sdcar[i] = TMath::Sqrt(0.150*0.150+(100*track->GetC())*(100*track->GetC()));
5381 if (TMath::Abs(track->GetTgl())>1) sdcar[i]*=2.5;
5382 cdcar[i] = TMath::Exp((TMath::Abs(track->GetC())-0.0106)*525.3);
5383 pulldcar[i] = (dca[i]-cdcar[i])/sdcar[i];
5384 pulldcaz[i] = (z0[i]-zvertex)/sdcar[i];
5385 pulldca[i] = TMath::Sqrt(pulldcar[i]*pulldcar[i]+pulldcaz[i]*pulldcaz[i]);
5386 if (track->TPCrPID(1)+track->TPCrPID(2)+track->TPCrPID(3)>0.5) {
5387 if (pulldca[i]<3.) isPrim[i]=kTRUE; //pion, muon and Kaon 3 sigma cut
5389 if (track->TPCrPID(4)>0.5) {
5390 if (pulldca[i]<0.5) isPrim[i]=kTRUE; //proton 0.5 sigma cut
5392 if (track->TPCrPID(0)>0.4) {
5393 isPrim[i]=kFALSE; //electron no sigma cut
5400 Int_t ncandidates =0;
5403 Double_t phase[2][2],radius[2];
5409 Float_t fprimvertex[3]={GetX(),GetY(),GetZ()};
5411 Double_t cradius0 = 10*10;
5412 Double_t cradius1 = 200*200;
5415 Double_t cpointAngle = 0.95;
5417 Double_t delta[2]={10000,10000};
5418 for (Int_t i =0;i<nentries;i++){
5419 if (sign[i]==0) continue;
5420 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5421 if (!track0) continue;
5422 if (AliTPCReconstructor::StreamLevel()>1){
5423 TTreeSRedirector &cstream = *fDebugStreamer;
5428 "zvertex="<<zvertex<<
5429 "sdcar0="<<sdcar[i]<<
5430 "cdcar0="<<cdcar[i]<<
5431 "pulldcar0="<<pulldcar[i]<<
5432 "pulldcaz0="<<pulldcaz[i]<<
5433 "pulldca0="<<pulldca[i]<<
5434 "isPrim="<<isPrim[i]<<
5438 if (track0->GetSigned1Pt()<0) continue;
5439 if (track0->GetKinkIndex(0)>0||isPrim[i]) continue; //daughter kink
5441 if (TMath::Abs(helixes[i].GetHelix(4))<0.000000001) continue;
5446 for (Int_t j =0;j<nentries;j++){
5447 AliTPCseed * track1 = (AliTPCseed*)array->At(j);
5448 if (!track1) continue;
5449 if (track1->GetKinkIndex(0)>0 || isPrim[j]) continue; //daughter kink
5450 if (sign[j]*sign[i]>0) continue;
5451 if (TMath::Abs(helixes[j].GetHelix(4))<0.000001) continue;
5452 if (track0->GetCircular()+track1->GetCircular()>1) continue; //circling -returning track
5455 // DCA to prim vertex cut
5461 Int_t npoints = helixes[i].GetRPHIintersections(helixes[j], phase, radius,cdist2);
5462 if (npoints<1) continue;
5466 if (radius[0]<cradius0||radius[0]>cradius1) continue;
5467 helixes[i].LinearDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta[0]);
5468 if (delta[0]>cdist1) continue;
5471 if (TMath::Max(radius[0],radius[1])<cradius0|| TMath::Min(radius[0],radius[1])>cradius1) continue;
5472 helixes[i].LinearDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta[0]);
5473 helixes[i].LinearDCA(helixes[j],phase[1][0],phase[1][1],radius[1],delta[1]);
5474 if (delta[1]<delta[0]) iclosest=1;
5475 if (delta[iclosest]>cdist1) continue;
5477 helixes[i].ParabolicDCA(helixes[j],phase[iclosest][0],phase[iclosest][1],radius[iclosest],delta[iclosest]);
5478 if (radius[iclosest]<cradius0 || radius[iclosest]>cradius1 || delta[iclosest]>cdist1) continue;
5480 Double_t pointAngle = helixes[i].GetPointAngle(helixes[j],phase[iclosest],fprimvertex);
5481 if (pointAngle<cpointAngle) continue;
5483 Bool_t isGamma = kFALSE;
5484 vertex.SetParamP(*track0); //track0 - plus
5485 vertex.SetParamN(*track1); //track1 - minus
5486 vertex.Update(fprimvertex);
5487 if (track0->TPCrPID(0)>0.3&&track1->TPCrPID(0)>0.3&&vertex.GetAnglep()[2]<0.15) isGamma=kTRUE; // gamma conversion candidate
5488 Double_t pointAngle2 = vertex.GetV0CosineOfPointingAngle();
5490 if (vertex.GetV0CosineOfPointingAngle()<cpointAngle && (!isGamma)) continue;// point angle cut
5491 //Bo: if (vertex.GetDist2()>2&&(!isGamma)) continue; // point angle cut
5492 if (vertex.GetDcaV0Daughters()>2&&(!isGamma)) continue;//Bo: // point angle cut
5493 Float_t sigmae = 0.15*0.15;
5494 if (vertex.GetRr()<80)
5495 sigmae += (sdcar[i]*sdcar[i]+sdcar[j]*sdcar[j])*(1.-vertex.GetRr()/80.)*(1.-vertex.GetRr()/80.);
5496 sigmae+= TMath::Sqrt(sigmae);
5497 //Bo: if (vertex.GetDist2()/sigmae>3.&&(!isGamma)) continue;
5498 if (vertex.GetDcaV0Daughters()/sigmae>3.&&(!isGamma)) continue;
5499 Float_t densb0=0,densb1=0,densa0=0,densa1=0;
5500 Int_t row0 = GetRowNumber(vertex.GetRr());
5502 //Bo: if (vertex.GetDist2()>0.2) continue;
5503 if (vertex.GetDcaV0Daughters()>0.2) continue;
5504 densb0 = track0->Density2(0,row0-5);
5505 densb1 = track1->Density2(0,row0-5);
5506 if (densb0>0.3|| densb1>0.3) continue; //clusters before vertex
5507 densa0 = track0->Density2(row0+5,row0+40);
5508 densa1 = track1->Density2(row0+5,row0+40);
5509 if ((densa0<0.4|| densa1<0.4)&&(!isGamma)) continue; //missing clusters after vertex
5512 densa0 = track0->Density2(0,40); //cluster density
5513 densa1 = track1->Density2(0,40); //cluster density
5514 if ((vertex.GetRr()<80&&densa0+densa1<1.)&&(!isGamma)) continue;
5516 //Bo: vertex.SetLab(0,track0->GetLabel());
5517 //Bo: vertex.SetLab(1,track1->GetLabel());
5518 vertex.SetChi2After((densa0+densa1)*0.5);
5519 vertex.SetChi2Before((densb0+densb1)*0.5);
5520 vertex.SetIndex(0,i);
5521 vertex.SetIndex(1,j);
5522 //Bo: vertex.SetStatus(1); // TPC v0 candidate
5523 vertex.SetOnFlyStatus(2);//Bo: // TPC v0 candidate
5524 //Bo: vertex.SetRp(track0->TPCrPIDs());
5525 //Bo: vertex.SetRm(track1->TPCrPIDs());
5526 tpcv0s->AddLast(new AliESDv0(vertex));
5529 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
5530 Double_t radiusm= (delta[0]<delta[1])? TMath::Sqrt(radius[0]):TMath::Sqrt(radius[1]);
5531 Double_t deltam= (delta[0]<delta[1])? TMath::Sqrt(delta[0]):TMath::Sqrt(delta[1]);
5532 if (AliTPCReconstructor::StreamLevel()>1) {
5533 Int_t lab0=track0->GetLabel();
5534 Int_t lab1=track1->GetLabel();
5535 Char_t c0=track0->GetCircular();
5536 Char_t c1=track1->GetCircular();
5537 TTreeSRedirector &cstream = *fDebugStreamer;
5540 "vertex.="<<&vertex<<
5543 "Helix0.="<<&helixes[i]<<
5546 "Helix1.="<<&helixes[j]<<
5547 "pointAngle="<<pointAngle<<
5548 "pointAngle2="<<pointAngle2<<
5553 "zvertex="<<zvertex<<
5556 "npoints="<<npoints<<
5557 "radius0="<<radius[0]<<
5558 "delta0="<<delta[0]<<
5559 "radius1="<<radius[1]<<
5560 "delta1="<<delta[1]<<
5561 "radiusm="<<radiusm<<
5563 "sdcar0="<<sdcar[i]<<
5564 "sdcar1="<<sdcar[j]<<
5565 "cdcar0="<<cdcar[i]<<
5566 "cdcar1="<<cdcar[j]<<
5567 "pulldcar0="<<pulldcar[i]<<
5568 "pulldcar1="<<pulldcar[j]<<
5569 "pulldcaz0="<<pulldcaz[i]<<
5570 "pulldcaz1="<<pulldcaz[j]<<
5571 "pulldca0="<<pulldca[i]<<
5572 "pulldca1="<<pulldca[j]<<
5582 Float_t *quality = new Float_t[ncandidates];
5583 Int_t *indexes = new Int_t[ncandidates];
5585 for (Int_t i=0;i<ncandidates;i++){
5587 AliESDv0 *v0 = (AliESDv0*)tpcv0s->At(i);
5588 quality[i] = 1./(1.00001-v0->GetV0CosineOfPointingAngle()); //base point angle
5589 // quality[i] /= (0.5+v0->GetDist2());
5590 // quality[i] *= v0->GetChi2After(); //density factor
5592 Int_t index0 = v0->GetIndex(0);
5593 Int_t index1 = v0->GetIndex(1);
5594 //Bo: Double_t minpulldca = TMath::Min(2.+pulldca[v0->GetIndex(0)],(2.+pulldca[v0->GetIndex(1)]) ); //pull
5595 Double_t minpulldca = TMath::Min(2.+pulldca[index0],(2.+pulldca[index1]) );//Bo:
5599 AliTPCseed * track0 = (AliTPCseed*)array->At(index0);
5600 AliTPCseed * track1 = (AliTPCseed*)array->At(index1);
5601 if (track0->TPCrPID(0)>0.3&&track1->TPCrPID(0)>0.3&&v0->GetAnglep()[2]<0.15) quality[i]+=1000000; // gamma conversion candidate
5602 if (track0->TPCrPID(4)>0.9||(track1->TPCrPID(4)>0.9&&minpulldca>4)) quality[i]*=10; // lambda candidate candidate
5605 TMath::Sort(ncandidates,quality,indexes,kTRUE);
5608 for (Int_t i=0;i<ncandidates;i++){
5609 AliESDv0 * v0 = (AliESDv0*)tpcv0s->At(indexes[i]);
5611 Int_t index0 = v0->GetIndex(0);
5612 Int_t index1 = v0->GetIndex(1);
5613 AliTPCseed * track0 = (AliTPCseed*)array->At(index0);
5614 AliTPCseed * track1 = (AliTPCseed*)array->At(index1);
5615 if (!track0||!track1) {
5619 Bool_t accept =kTRUE; //default accept
5620 Int_t *v0indexes0 = track0->GetV0Indexes();
5621 Int_t *v0indexes1 = track1->GetV0Indexes();
5623 Int_t order0 = (v0indexes0[0]!=0) ? 1:0;
5624 Int_t order1 = (v0indexes1[0]!=0) ? 1:0;
5625 if (v0indexes0[1]!=0) order0 =2;
5626 if (v0indexes1[1]!=0) order1 =2;
5628 if (v0indexes0[2]!=0) {order0=3; accept=kFALSE;}
5629 if (v0indexes0[2]!=0) {order1=3; accept=kFALSE;}
5631 AliESDv0 * v02 = v0;
5633 //Bo: v0->SetOrder(0,order0);
5634 //Bo: v0->SetOrder(1,order1);
5635 //Bo: v0->SetOrder(1,order0+order1);
5636 v0->SetOnFlyStatus(kTRUE);
5637 Int_t index = esd->AddV0(v0);
5638 v02 = esd->GetV0(index);
5639 v0indexes0[order0]=index;
5640 v0indexes1[order1]=index;
5644 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
5645 if (AliTPCReconstructor::StreamLevel()>1) {
5646 Int_t lab0=track0->GetLabel();
5647 Int_t lab1=track1->GetLabel();
5648 TTreeSRedirector &cstream = *fDebugStreamer;
5657 "dca0="<<dca[index0]<<
5658 "dca1="<<dca[index1]<<
5662 "quality="<<quality[i]<<
5663 "pulldca0="<<pulldca[index0]<<
5664 "pulldca1="<<pulldca[index1]<<
5688 printf("TPC V0 finder : naccepted\t%d\tncandidates\t%d\tntracks\t%d\tnall\t%d\n",naccepted,ncandidates,ntracks,nall);
5692 Int_t AliTPCtrackerMI::RefitKink(AliTPCseed &mother, AliTPCseed &daughter, const AliESDkink &knk)
5695 // refit kink towards to the vertex
5698 AliKink &kink=(AliKink &)knk;
5700 Int_t row0 = GetRowNumber(kink.GetR());
5701 FollowProlongation(mother,0);
5702 mother.Reset(kFALSE);
5704 FollowProlongation(daughter,row0);
5705 daughter.Reset(kFALSE);
5706 FollowBackProlongation(daughter,158);
5707 daughter.Reset(kFALSE);
5708 Int_t first = TMath::Max(row0-20,30);
5709 Int_t last = TMath::Min(row0+20,140);
5711 const Int_t kNdiv =5;
5712 AliTPCseed param0[kNdiv]; // parameters along the track
5713 AliTPCseed param1[kNdiv]; // parameters along the track
5714 AliKink kinks[kNdiv]; // corresponding kink parameters
5717 for (Int_t irow=0; irow<kNdiv;irow++){
5718 rows[irow] = first +((last-first)*irow)/(kNdiv-1);
5720 // store parameters along the track
5722 for (Int_t irow=0;irow<kNdiv;irow++){
5723 FollowBackProlongation(mother, rows[irow]);
5724 FollowProlongation(daughter,rows[kNdiv-1-irow]);
5725 param0[irow] = mother;
5726 param1[kNdiv-1-irow] = daughter;
5730 for (Int_t irow=0; irow<kNdiv-1;irow++){
5731 if (param0[irow].GetNumberOfClusters()<kNdiv||param1[irow].GetNumberOfClusters()<kNdiv) continue;
5732 kinks[irow].SetMother(param0[irow]);
5733 kinks[irow].SetDaughter(param1[irow]);
5734 kinks[irow].Update();
5737 // choose kink with best "quality"
5739 Double_t mindist = 10000;
5740 for (Int_t irow=0;irow<kNdiv;irow++){
5741 if (param0[irow].GetNumberOfClusters()<20||param1[irow].GetNumberOfClusters()<20) continue;
5742 if (TMath::Abs(kinks[irow].GetR())>240.) continue;
5743 if (TMath::Abs(kinks[irow].GetR())<100.) continue;
5745 Float_t normdist = TMath::Abs(param0[irow].GetX()-kinks[irow].GetR())*(0.1+kink.GetDistance());
5746 normdist/= (param0[irow].GetNumberOfClusters()+param1[irow].GetNumberOfClusters()+40.);
5747 if (normdist < mindist){
5753 if (index==-1) return 0;
5756 param0[index].Reset(kTRUE);
5757 FollowProlongation(param0[index],0);
5759 mother = param0[index];
5760 daughter = param1[index]; // daughter in vertex
5762 kink.SetMother(mother);
5763 kink.SetDaughter(daughter);
5765 kink.SetTPCRow0(GetRowNumber(kink.GetR()));
5766 kink.SetTPCncls(param0[index].GetNumberOfClusters(),0);
5767 kink.SetTPCncls(param1[index].GetNumberOfClusters(),1);
5768 kink.SetLabel(CookLabel(&mother,0.4, 0,kink.GetTPCRow0()),0);
5769 kink.SetLabel(CookLabel(&daughter,0.4, kink.GetTPCRow0(),160),1);
5770 mother.SetLabel(kink.GetLabel(0));
5771 daughter.SetLabel(kink.GetLabel(1));
5777 void AliTPCtrackerMI::UpdateKinkQualityM(AliTPCseed * seed){
5779 // update Kink quality information for mother after back propagation
5781 if (seed->GetKinkIndex(0)>=0) return;
5782 for (Int_t ikink=0;ikink<3;ikink++){
5783 Int_t index = seed->GetKinkIndex(ikink);
5784 if (index>=0) break;
5785 index = TMath::Abs(index)-1;
5786 AliESDkink * kink = fEvent->GetKink(index);
5787 kink->SetTPCDensity(-1,0,0);
5788 kink->SetTPCDensity(1,0,1);
5790 Int_t row0 = kink->GetTPCRow0() - 2 - Int_t( 0.5/ (0.05+kink->GetAngle(2)));
5791 if (row0<15) row0=15;
5793 Int_t row1 = kink->GetTPCRow0() + 2 + Int_t( 0.5/ (0.05+kink->GetAngle(2)));
5794 if (row1>145) row1=145;
5796 Int_t found,foundable,shared;
5797 seed->GetClusterStatistic(0,row0, found, foundable,shared,kFALSE);
5798 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),0,0);
5799 seed->GetClusterStatistic(row1,155, found, foundable,shared,kFALSE);
5800 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),0,1);
5805 void AliTPCtrackerMI::UpdateKinkQualityD(AliTPCseed * seed){
5807 // update Kink quality information for daughter after refit
5809 if (seed->GetKinkIndex(0)<=0) return;
5810 for (Int_t ikink=0;ikink<3;ikink++){
5811 Int_t index = seed->GetKinkIndex(ikink);
5812 if (index<=0) break;
5813 index = TMath::Abs(index)-1;
5814 AliESDkink * kink = fEvent->GetKink(index);
5815 kink->SetTPCDensity(-1,1,0);
5816 kink->SetTPCDensity(-1,1,1);
5818 Int_t row0 = kink->GetTPCRow0() -2 - Int_t( 0.5/ (0.05+kink->GetAngle(2)));
5819 if (row0<15) row0=15;
5821 Int_t row1 = kink->GetTPCRow0() +2 + Int_t( 0.5/ (0.05+kink->GetAngle(2)));
5822 if (row1>145) row1=145;
5824 Int_t found,foundable,shared;
5825 seed->GetClusterStatistic(0,row0, found, foundable,shared,kFALSE);
5826 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),1,0);
5827 seed->GetClusterStatistic(row1,155, found, foundable,shared,kFALSE);
5828 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),1,1);
5834 Int_t AliTPCtrackerMI::CheckKinkPoint(AliTPCseed*seed,AliTPCseed &mother, AliTPCseed &daughter, const AliESDkink &knk)
5837 // check kink point for given track
5838 // if return value=0 kink point not found
5839 // otherwise seed0 correspond to mother particle
5840 // seed1 correspond to daughter particle
5841 // kink parameter of kink point
5842 AliKink &kink=(AliKink &)knk;
5844 Int_t middlerow = (seed->GetFirstPoint()+seed->GetLastPoint())/2;
5845 Int_t first = seed->GetFirstPoint();
5846 Int_t last = seed->GetLastPoint();
5847 if (last-first<20) return 0; // shortest length - 2*30 = 60 pad-rows
5850 AliTPCseed *seed1 = ReSeed(seed,middlerow+20, kTRUE); //middle of chamber
5851 if (!seed1) return 0;
5852 FollowProlongation(*seed1,seed->GetLastPoint()-20);
5853 seed1->Reset(kTRUE);
5854 FollowProlongation(*seed1,158);
5855 seed1->Reset(kTRUE);
5856 last = seed1->GetLastPoint();
5858 AliTPCseed *seed0 = new AliTPCseed(*seed);
5859 seed0->Reset(kFALSE);
5862 AliTPCseed param0[20]; // parameters along the track
5863 AliTPCseed param1[20]; // parameters along the track
5864 AliKink kinks[20]; // corresponding kink parameters
5866 for (Int_t irow=0; irow<20;irow++){
5867 rows[irow] = first +((last-first)*irow)/19;
5869 // store parameters along the track
5871 for (Int_t irow=0;irow<20;irow++){
5872 FollowBackProlongation(*seed0, rows[irow]);
5873 FollowProlongation(*seed1,rows[19-irow]);
5874 param0[irow] = *seed0;
5875 param1[19-irow] = *seed1;
5879 for (Int_t irow=0; irow<19;irow++){
5880 kinks[irow].SetMother(param0[irow]);
5881 kinks[irow].SetDaughter(param1[irow]);
5882 kinks[irow].Update();
5885 // choose kink with biggest change of angle
5887 Double_t maxchange= 0;
5888 for (Int_t irow=1;irow<19;irow++){
5889 if (TMath::Abs(kinks[irow].GetR())>240.) continue;
5890 if (TMath::Abs(kinks[irow].GetR())<110.) continue;
5891 Float_t quality = TMath::Abs(kinks[irow].GetAngle(2))/(3.+TMath::Abs(kinks[irow].GetR()-param0[irow].GetX()));
5892 if ( quality > maxchange){
5893 maxchange = quality;
5900 if (index<0) return 0;
5902 Int_t row0 = GetRowNumber(kinks[index].GetR()); //row 0 estimate
5903 seed0 = new AliTPCseed(param0[index]);
5904 seed1 = new AliTPCseed(param1[index]);
5905 seed0->Reset(kFALSE);
5906 seed1->Reset(kFALSE);
5907 seed0->ResetCovariance(10.);
5908 seed1->ResetCovariance(10.);
5909 FollowProlongation(*seed0,0);
5910 FollowBackProlongation(*seed1,158);
5911 mother = *seed0; // backup mother at position 0
5912 seed0->Reset(kFALSE);
5913 seed1->Reset(kFALSE);
5914 seed0->ResetCovariance(10.);
5915 seed1->ResetCovariance(10.);
5917 first = TMath::Max(row0-20,0);
5918 last = TMath::Min(row0+20,158);
5920 for (Int_t irow=0; irow<20;irow++){
5921 rows[irow] = first +((last-first)*irow)/19;
5923 // store parameters along the track
5925 for (Int_t irow=0;irow<20;irow++){
5926 FollowBackProlongation(*seed0, rows[irow]);
5927 FollowProlongation(*seed1,rows[19-irow]);
5928 param0[irow] = *seed0;
5929 param1[19-irow] = *seed1;
5933 for (Int_t irow=0; irow<19;irow++){
5934 kinks[irow].SetMother(param0[irow]);
5935 kinks[irow].SetDaughter(param1[irow]);
5936 // param0[irow].Dump();
5937 //param1[irow].Dump();
5938 kinks[irow].Update();
5941 // choose kink with biggest change of angle
5944 for (Int_t irow=0;irow<20;irow++){
5945 if (TMath::Abs(kinks[irow].GetR())>250.) continue;
5946 if (TMath::Abs(kinks[irow].GetR())<90.) continue;
5947 Float_t quality = TMath::Abs(kinks[irow].GetAngle(2))/(3.+TMath::Abs(kinks[irow].GetR()-param0[irow].GetX()));
5948 if ( quality > maxchange){
5949 maxchange = quality;
5956 if (index==-1 || param0[index].GetNumberOfClusters()+param1[index].GetNumberOfClusters()<100){
5962 // Float_t anglesigma = TMath::Sqrt(param0[index].fC22+param0[index].fC33+param1[index].fC22+param1[index].fC33);
5964 kink.SetMother(param0[index]);
5965 kink.SetDaughter(param1[index]);
5968 Double_t chi2P2 = param0[index].GetParameter()[2]-param1[index].GetParameter()[2];
5970 chi2P2/=param0[index].GetCovariance()[5]+param1[index].GetCovariance()[5];
5971 Double_t chi2P3 = param0[index].GetParameter()[3]-param1[index].GetParameter()[3];
5973 chi2P3/=param0[index].GetCovariance()[9]+param1[index].GetCovariance()[9];
5975 if (AliTPCReconstructor::StreamLevel()>1) {
5976 (*fDebugStreamer)<<"kinkHpt"<<
5979 "p0.="<<¶m0[index]<<
5980 "p1.="<<¶m1[index]<<
5984 if ( chi2P2+chi2P3<AliTPCReconstructor::GetRecoParam()->GetKinkAngleCutChi2(0)){
5991 row0 = GetRowNumber(kink.GetR());
5992 kink.SetTPCRow0(row0);
5993 kink.SetLabel(CookLabel(seed0,0.5,0,row0),0);
5994 kink.SetLabel(CookLabel(seed1,0.5,row0,158),1);
5995 kink.SetIndex(-10,0);
5996 kink.SetIndex(int(param0[index].GetNumberOfClusters()+param1[index].GetNumberOfClusters()),1);
5997 kink.SetTPCncls(param0[index].GetNumberOfClusters(),0);
5998 kink.SetTPCncls(param1[index].GetNumberOfClusters(),1);
6001 // new (&mother) AliTPCseed(param0[index]);
6002 daughter = param1[index];
6003 daughter.SetLabel(kink.GetLabel(1));
6004 param0[index].Reset(kTRUE);
6005 FollowProlongation(param0[index],0);
6006 mother = param0[index];
6007 mother.SetLabel(kink.GetLabel(0));
6008 if ( chi2P2+chi2P3<AliTPCReconstructor::GetRecoParam()->GetKinkAngleCutChi2(1)){
6020 AliTPCseed* AliTPCtrackerMI::ReSeed(AliTPCseed *t)
6023 // reseed - refit - track
6026 // Int_t last = fSectors->GetNRows()-1;
6028 if (fSectors == fOuterSec){
6029 first = TMath::Max(first, t->GetFirstPoint()-fInnerSec->GetNRows());
6033 first = t->GetFirstPoint();
6035 AliTPCseed * seed = MakeSeed(t,0.1,0.5,0.9);
6036 FollowBackProlongation(*t,fSectors->GetNRows()-1);
6038 FollowProlongation(*t,first);
6048 //_____________________________________________________________________________
6049 Int_t AliTPCtrackerMI::ReadSeeds(const TFile *inp) {
6050 //-----------------------------------------------------------------
6051 // This function reades track seeds.
6052 //-----------------------------------------------------------------
6053 TDirectory *savedir=gDirectory;
6055 TFile *in=(TFile*)inp;
6056 if (!in->IsOpen()) {
6057 cerr<<"AliTPCtrackerMI::ReadSeeds(): input file is not open !\n";
6062 TTree *seedTree=(TTree*)in->Get("Seeds");
6064 cerr<<"AliTPCtrackerMI::ReadSeeds(): ";
6065 cerr<<"can't get a tree with track seeds !\n";
6068 AliTPCtrack *seed=new AliTPCtrack;
6069 seedTree->SetBranchAddress("tracks",&seed);
6071 if (fSeeds==0) fSeeds=new TObjArray(15000);
6073 Int_t n=(Int_t)seedTree->GetEntries();
6074 for (Int_t i=0; i<n; i++) {
6075 seedTree->GetEvent(i);
6076 fSeeds->AddLast(new AliTPCseed(*seed/*,seed->GetAlpha()*/));
6085 Int_t AliTPCtrackerMI::Clusters2Tracks (AliESDEvent *const esd)
6088 if (fSeeds) DeleteSeeds();
6091 if (!fSeeds) return 1;
6098 //_____________________________________________________________________________
6099 Int_t AliTPCtrackerMI::Clusters2Tracks() {
6100 //-----------------------------------------------------------------
6101 // This is a track finder.
6102 //-----------------------------------------------------------------
6103 TDirectory *savedir=gDirectory;
6107 fSeeds = Tracking();
6110 Info("Clusters2Tracks","Time for tracking: \t");timer.Print();timer.Start();
6112 //activate again some tracks
6113 for (Int_t i=0; i<fSeeds->GetEntriesFast(); i++) {
6114 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
6116 Int_t nc=t.GetNumberOfClusters();
6118 delete fSeeds->RemoveAt(i);
6122 if (pt->GetRemoval()==10) {
6123 if (pt->GetDensityFirst(20)>0.8 || pt->GetDensityFirst(30)>0.8 || pt->GetDensityFirst(40)>0.7)
6124 pt->Desactivate(10); // make track again active
6126 pt->Desactivate(20);
6127 delete fSeeds->RemoveAt(i);
6132 RemoveUsed2(fSeeds,0.85,0.85,0);
6133 if (AliTPCReconstructor::GetRecoParam()->GetDoKinks()) FindKinks(fSeeds,fEvent);
6134 //FindCurling(fSeeds, fEvent,0);
6135 if (AliTPCReconstructor::StreamLevel()>5) FindMultiMC(fSeeds, fEvent,-1); // find multi found tracks
6136 RemoveUsed2(fSeeds,0.5,0.4,20);
6137 FindSplitted(fSeeds, fEvent,0); // find multi found tracks
6138 if (AliTPCReconstructor::StreamLevel()>5) FindMultiMC(fSeeds, fEvent,0); // find multi found tracks
6141 // // refit short tracks
6143 Int_t nseed=fSeeds->GetEntriesFast();
6146 for (Int_t i=0; i<nseed; i++) {
6147 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
6149 Int_t nc=t.GetNumberOfClusters();
6151 delete fSeeds->RemoveAt(i);
6154 CookLabel(pt,0.1); //For comparison only
6155 //if ((pt->IsActive() || (pt->fRemoval==10) )&& nc>50 &&pt->GetNumberOfClusters()>0.4*pt->fNFoundable){
6156 if ((pt->IsActive() || (pt->GetRemoval()==10) )){
6158 if (fDebug>0) cerr<<found<<'\r';
6162 delete fSeeds->RemoveAt(i);
6166 //RemoveOverlap(fSeeds,0.99,7,kTRUE);
6168 //RemoveUsed(fSeeds,0.9,0.9,6);
6170 nseed=fSeeds->GetEntriesFast();
6172 for (Int_t i=0; i<nseed; i++) {
6173 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
6175 Int_t nc=t.GetNumberOfClusters();
6177 delete fSeeds->RemoveAt(i);
6181 t.CookdEdx(0.02,0.6);
6182 // CheckKinkPoint(&t,0.05);
6183 //if ((pt->IsActive() || (pt->fRemoval==10) )&& nc>50 &&pt->GetNumberOfClusters()>0.4*pt->fNFoundable){
6184 if ((pt->IsActive() || (pt->GetRemoval()==10) )){
6192 delete fSeeds->RemoveAt(i);
6193 //AliTPCseed * seed1 = ReSeed(pt,0.05,0.5,1);
6195 // FollowProlongation(*seed1,0);
6196 // Int_t n = seed1->GetNumberOfClusters();
6197 // printf("fP4\t%f\t%f\n",seed1->GetC(),pt->GetC());
6198 // printf("fN\t%d\t%d\n", seed1->GetNumberOfClusters(),pt->GetNumberOfClusters());
6201 //AliTPCseed * seed2 = ReSeed(pt,0.95,0.5,0.05);
6205 SortTracks(fSeeds, 1);
6209 PrepareForBackProlongation(fSeeds,5.);
6210 PropagateBack(fSeeds);
6211 printf("Time for back propagation: \t");timer.Print();timer.Start();
6215 PrepareForProlongation(fSeeds,5.);
6216 PropagateForward2(fSeeds);
6218 printf("Time for FORWARD propagation: \t");timer.Print();timer.Start();
6219 // RemoveUsed(fSeeds,0.7,0.7,6);
6220 //RemoveOverlap(fSeeds,0.9,7,kTRUE);
6222 nseed=fSeeds->GetEntriesFast();
6224 for (Int_t i=0; i<nseed; i++) {
6225 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
6227 Int_t nc=t.GetNumberOfClusters();
6229 delete fSeeds->RemoveAt(i);
6232 t.CookdEdx(0.02,0.6);
6233 // CookLabel(pt,0.1); //For comparison only
6234 //if ((pt->IsActive() || (pt->fRemoval==10) )&& nc>50 &&pt->GetNumberOfClusters()>0.4*pt->fNFoundable){
6235 if ((pt->IsActive() || (pt->fRemoval==10) )){
6236 cerr<<found++<<'\r';
6239 delete fSeeds->RemoveAt(i);
6244 // fNTracks = found;
6246 Info("Clusters2Tracks","Time for overlap removal, track writing and dedx cooking: \t"); timer.Print();timer.Start();
6249 // cerr<<"Number of found tracks : "<<"\t"<<found<<endl;
6250 Info("Clusters2Tracks","Number of found tracks %d",found);
6252 // UnloadClusters();
6257 void AliTPCtrackerMI::Tracking(TObjArray * arr)
6260 // tracking of the seeds
6263 fSectors = fOuterSec;
6264 ParallelTracking(arr,150,63);
6265 fSectors = fOuterSec;
6266 ParallelTracking(arr,63,0);
6269 TObjArray * AliTPCtrackerMI::Tracking(Int_t seedtype, Int_t i1, Int_t i2, Float_t cuts[4], Float_t dy, Int_t dsec)
6274 TObjArray * arr = new TObjArray;
6276 fSectors = fOuterSec;
6279 for (Int_t sec=0;sec<fkNOS;sec++){
6280 if (seedtype==3) MakeSeeds3(arr,sec,i1,i2,cuts,dy, dsec);
6281 if (seedtype==4) MakeSeeds5(arr,sec,i1,i2,cuts,dy);
6282 if (seedtype==2) MakeSeeds2(arr,sec,i1,i2,cuts,dy);
6285 Info("Tracking","\nSeeding - %d\t%d\t%d\t%d\n",seedtype,i1,i2,arr->GetEntriesFast());
6297 TObjArray * AliTPCtrackerMI::Tracking()
6301 if (AliTPCReconstructor::GetRecoParam()->GetSpecialSeeding()) return TrackingSpecial();
6304 Int_t nup=fOuterSec->GetNRows()+fInnerSec->GetNRows();
6306 TObjArray * seeds = new TObjArray;
6315 Float_t fnumber = 3.0;
6316 Float_t fdensity = 3.0;
6321 for (Int_t delta = 0; delta<18; delta+=6){
6325 arr = Tracking(3,nup-1-delta,nup-1-delta-gap,cuts,-1,1);
6326 SumTracks(seeds,arr);
6327 SignClusters(seeds,fnumber,fdensity);
6329 for (Int_t i=2;i<6;i+=2){
6330 // seed high pt tracks
6333 arr = Tracking(3,nup-i-delta,nup-i-delta-gap,cuts,-1,0);
6334 SumTracks(seeds,arr);
6335 SignClusters(seeds,fnumber,fdensity);
6340 // RemoveUsed(seeds,0.9,0.9,1);
6341 // UnsignClusters();
6342 // SignClusters(seeds,fnumber,fdensity);
6346 for (Int_t delta = 20; delta<120; delta+=10){
6348 // seed high pt tracks
6352 arr = Tracking(3,nup-delta,nup-delta-gap,cuts,-1);
6353 SumTracks(seeds,arr);
6354 SignClusters(seeds,fnumber,fdensity);
6359 arr = Tracking(3,nup-delta-5,nup-delta-5-gap,cuts,-1);
6360 SumTracks(seeds,arr);
6361 SignClusters(seeds,fnumber,fdensity);
6372 Info("Tracking()","\n\nPrimary seeding\t%d\n\n",seeds->GetEntriesFast());
6376 // RemoveUsed(seeds,0.75,0.75,1);
6378 //SignClusters(seeds,fnumber,fdensity);
6387 arr = Tracking(4,nup-1,nup-1-gap,cuts,-1);
6388 SumTracks(seeds,arr);
6389 SignClusters(seeds,fnumber,fdensity);
6391 arr = Tracking(4,nup-2,nup-2-gap,cuts,-1);
6392 SumTracks(seeds,arr);
6393 SignClusters(seeds,fnumber,fdensity);
6395 arr = Tracking(4,nup-3,nup-3-gap,cuts,-1);
6396 SumTracks(seeds,arr);
6397 SignClusters(seeds,fnumber,fdensity);
6401 for (Int_t delta = 3; delta<30; delta+=5){
6407 arr = Tracking(4,nup-1-delta,nup-1-delta-gap,cuts,-1);
6408 SumTracks(seeds,arr);
6409 SignClusters(seeds,fnumber,fdensity);
6411 arr = Tracking(4,nup-3-delta,nup-5-delta-gap,cuts,4);
6412 SumTracks(seeds,arr);
6413 SignClusters(seeds,fnumber,fdensity);
6424 Int_t fLastSeedRowSec=AliTPCReconstructor::GetRecoParam()->GetLastSeedRowSec();
6427 for (Int_t delta = 30; delta<fLastSeedRowSec; delta+=10){
6433 arr = Tracking(4,nup-1-delta,nup-1-delta-gap,cuts,-1);
6434 SumTracks(seeds,arr);
6435 SignClusters(seeds,fnumber,fdensity);
6437 arr = Tracking(4,nup-5-delta,nup-5-delta-gap,cuts,5 );
6438 SumTracks(seeds,arr);
6439 SignClusters(seeds,fnumber,fdensity);
6443 Info("Tracking()","\n\nSecondary seeding\t%d\n\n",seeds->GetEntriesFast());
6454 TObjArray * AliTPCtrackerMI::TrackingSpecial()
6457 // seeding adjusted for laser and cosmic tests - short tracks with big inclination angle
6458 // no primary vertex seeding tried
6462 Int_t nup=fOuterSec->GetNRows()+fInnerSec->GetNRows();
6464 TObjArray * seeds = new TObjArray;
6469 Float_t fnumber = 3.0;
6470 Float_t fdensity = 3.0;
6473 cuts[0] = AliTPCReconstructor::GetRecoParam()->GetMaxC(); // max curvature
6474 cuts[1] = 3.5; // max tan(phi) angle for seeding
6475 cuts[2] = 3.; // not used (cut on z primary vertex)
6476 cuts[3] = 3.5; // max tan(theta) angle for seeding
6478 for (Int_t delta = 0; nup-delta-gap-1>0; delta+=3){
6480 arr = Tracking(4,nup-1-delta,nup-1-delta-gap,cuts,-1);
6481 SumTracks(seeds,arr);
6482 SignClusters(seeds,fnumber,fdensity);
6486 Info("Tracking()","\n\nSecondary seeding\t%d\n\n",seeds->GetEntriesFast());
6497 void AliTPCtrackerMI::SumTracks(TObjArray *arr1,TObjArray *&arr2) const
6500 //sum tracks to common container
6501 //remove suspicious tracks
6502 Int_t nseed = arr2->GetEntriesFast();
6503 for (Int_t i=0;i<nseed;i++){
6504 AliTPCseed *pt=(AliTPCseed*)arr2->UncheckedAt(i);
6507 // remove tracks with too big curvature
6509 if (TMath::Abs(pt->GetC())>AliTPCReconstructor::GetRecoParam()->GetMaxC()){
6510 delete arr2->RemoveAt(i);
6513 // REMOVE VERY SHORT TRACKS
6514 if (pt->GetNumberOfClusters()<20){
6515 delete arr2->RemoveAt(i);
6518 // NORMAL ACTIVE TRACK
6519 if (pt->IsActive()){
6520 arr1->AddLast(arr2->RemoveAt(i));
6523 //remove not usable tracks
6524 if (pt->GetRemoval()!=10){
6525 delete arr2->RemoveAt(i);
6529 // ENABLE ONLY ENOUGH GOOD STOPPED TRACKS
6530 if (pt->GetDensityFirst(20)>0.8 || pt->GetDensityFirst(30)>0.8 || pt->GetDensityFirst(40)>0.7)
6531 arr1->AddLast(arr2->RemoveAt(i));
6533 delete arr2->RemoveAt(i);
6537 delete arr2; arr2 = 0;
6542 void AliTPCtrackerMI::ParallelTracking(TObjArray *const arr, Int_t rfirst, Int_t rlast)
6545 // try to track in parralel
6547 Int_t nseed=arr->GetEntriesFast();
6548 //prepare seeds for tracking
6549 for (Int_t i=0; i<nseed; i++) {
6550 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i), &t=*pt;
6552 if (!t.IsActive()) continue;
6553 // follow prolongation to the first layer
6554 if ( (fSectors ==fInnerSec) || (t.GetFirstPoint()-fkParam->GetNRowLow()>rfirst+1) )
6555 FollowProlongation(t, rfirst+1);
6560 for (Int_t nr=rfirst; nr>=rlast; nr--){
6561 if (nr<fInnerSec->GetNRows())
6562 fSectors = fInnerSec;
6564 fSectors = fOuterSec;
6565 // make indexes with the cluster tracks for given
6567 // find nearest cluster
6568 for (Int_t i=0; i<nseed; i++) {
6569 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i), &t=*pt;
6571 if (nr==80) pt->UpdateReference();
6572 if (!pt->IsActive()) continue;
6573 // if ( (fSectors ==fOuterSec) && (pt->fFirstPoint-fkParam->GetNRowLow())<nr) continue;
6574 if (pt->GetRelativeSector()>17) {
6577 UpdateClusters(t,nr);
6579 // prolonagate to the nearest cluster - if founded
6580 for (Int_t i=0; i<nseed; i++) {
6581 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
6583 if (!pt->IsActive()) continue;
6584 // if ((fSectors ==fOuterSec) && (pt->fFirstPoint-fkParam->GetNRowLow())<nr) continue;
6585 if (pt->GetRelativeSector()>17) {
6588 FollowToNextCluster(*pt,nr);
6593 void AliTPCtrackerMI::PrepareForBackProlongation(TObjArray *const arr,Float_t fac) const
6597 // if we use TPC track itself we have to "update" covariance
6599 Int_t nseed= arr->GetEntriesFast();
6600 for (Int_t i=0;i<nseed;i++){
6601 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
6605 //rotate to current local system at first accepted point
6606 Int_t index = pt->GetClusterIndex2(pt->GetFirstPoint());
6607 Int_t sec = (index&0xff000000)>>24;
6609 Float_t angle1 = fInnerSec->GetAlpha()*sec+fInnerSec->GetAlphaShift();
6610 if (angle1>TMath::Pi())
6611 angle1-=2.*TMath::Pi();
6612 Float_t angle2 = pt->GetAlpha();
6614 if (TMath::Abs(angle1-angle2)>0.001){
6615 pt->Rotate(angle1-angle2);
6616 //angle2 = pt->GetAlpha();
6617 //pt->fRelativeSector = pt->GetAlpha()/fInnerSec->GetAlpha();
6618 //if (pt->GetAlpha()<0)
6619 // pt->fRelativeSector+=18;
6620 //sec = pt->fRelativeSector;
6629 void AliTPCtrackerMI::PrepareForProlongation(TObjArray *const arr, Float_t fac) const
6633 // if we use TPC track itself we have to "update" covariance
6635 Int_t nseed= arr->GetEntriesFast();
6636 for (Int_t i=0;i<nseed;i++){
6637 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
6640 pt->SetFirstPoint(pt->GetLastPoint());
6648 Int_t AliTPCtrackerMI::PropagateBack(TObjArray *const arr)
6651 // make back propagation
6653 Int_t nseed= arr->GetEntriesFast();
6654 for (Int_t i=0;i<nseed;i++){
6655 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
6656 if (pt&& pt->GetKinkIndex(0)<=0) {
6657 //AliTPCseed *pt2 = new AliTPCseed(*pt);
6658 fSectors = fInnerSec;
6659 //FollowBackProlongation(*pt,fInnerSec->GetNRows()-1);
6660 //fSectors = fOuterSec;
6661 FollowBackProlongation(*pt,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1);
6662 //if (pt->GetNumberOfClusters()<(pt->fEsd->GetTPCclusters(0)) ){
6663 // Error("PropagateBack","Not prolonged track %d",pt->GetLabel());
6664 // FollowBackProlongation(*pt2,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1);
6667 if (pt&& pt->GetKinkIndex(0)>0) {
6668 AliESDkink * kink = fEvent->GetKink(pt->GetKinkIndex(0)-1);
6669 pt->SetFirstPoint(kink->GetTPCRow0());
6670 fSectors = fInnerSec;
6671 FollowBackProlongation(*pt,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1);
6679 Int_t AliTPCtrackerMI::PropagateForward2(TObjArray *const arr)
6682 // make forward propagation
6684 Int_t nseed= arr->GetEntriesFast();
6686 for (Int_t i=0;i<nseed;i++){
6687 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
6689 FollowProlongation(*pt,0);
6698 Int_t AliTPCtrackerMI::PropagateForward()
6701 // propagate track forward
6703 Int_t nseed = fSeeds->GetEntriesFast();
6704 for (Int_t i=0;i<nseed;i++){
6705 AliTPCseed *pt = (AliTPCseed*)fSeeds->UncheckedAt(i);
6707 AliTPCseed &t = *pt;
6708 Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
6709 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
6710 if (alpha < 0. ) alpha += 2.*TMath::Pi();
6711 t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN);
6715 fSectors = fOuterSec;
6716 ParallelTracking(fSeeds,fOuterSec->GetNRows()+fInnerSec->GetNRows()-1,fInnerSec->GetNRows());
6717 fSectors = fInnerSec;
6718 ParallelTracking(fSeeds,fInnerSec->GetNRows()-1,0);
6727 Int_t AliTPCtrackerMI::PropagateBack(AliTPCseed *const pt, Int_t row0, Int_t row1)
6730 // make back propagation, in between row0 and row1
6734 fSectors = fInnerSec;
6737 if (row1<fSectors->GetNRows())
6740 r1 = fSectors->GetNRows()-1;
6742 if (row0<fSectors->GetNRows()&& r1>0 )
6743 FollowBackProlongation(*pt,r1);
6744 if (row1<=fSectors->GetNRows())
6747 r1 = row1 - fSectors->GetNRows();
6748 if (r1<=0) return 0;
6749 if (r1>=fOuterSec->GetNRows()) return 0;
6750 fSectors = fOuterSec;
6751 return FollowBackProlongation(*pt,r1);
6759 void AliTPCtrackerMI::GetShape(AliTPCseed * seed, Int_t row)
6763 AliTPCClusterParam * clparam = AliTPCcalibDB::Instance()->GetClusterParam();
6764 Float_t zdrift = TMath::Abs((fkParam->GetZLength(0)-TMath::Abs(seed->GetZ())));
6765 Int_t type = (seed->GetSector() < fkParam->GetNSector()/2) ? 0: (row>126) ? 1:2;
6766 Double_t angulary = seed->GetSnp();
6767 angulary = angulary*angulary/((1.-angulary)*(1.+angulary));
6768 Double_t angularz = seed->GetTgl()*seed->GetTgl()*(1.+angulary);
6770 Double_t sigmay = clparam->GetRMS0(0,type,zdrift,TMath::Sqrt(TMath::Abs(angulary)));
6771 Double_t sigmaz = clparam->GetRMS0(1,type,zdrift,TMath::Sqrt(TMath::Abs(angularz)));
6772 seed->SetCurrentSigmaY2(sigmay*sigmay);
6773 seed->SetCurrentSigmaZ2(sigmaz*sigmaz);
6774 // Float_t sd2 = TMath::Abs((fkParam->GetZLength(0)-TMath::Abs(seed->GetZ())))*fkParam->GetDiffL()*fkParam->GetDiffL();
6775 // // Float_t padlength = fkParam->GetPadPitchLength(seed->fSector);
6776 // Float_t padlength = GetPadPitchLength(row);
6778 // Float_t sresy = (seed->GetSector() < fkParam->GetNSector()/2) ? 0.2 :0.3;
6779 // seed->SetCurrentSigmaY2(sd2+padlength*padlength*angulary/12.+sresy*sresy);
6781 // Float_t sresz = fkParam->GetZSigma();
6782 // seed->SetCurrentSigmaZ2(sd2+padlength*padlength*angularz*angularz*(1+angulary)/12.+sresz*sresz);
6784 Float_t wy = GetSigmaY(seed);
6785 Float_t wz = GetSigmaZ(seed);
6788 if (TMath::Abs(wy/seed->fCurrentSigmaY2-1)>0.0001 || TMath::Abs(wz/seed->fCurrentSigmaZ2-1)>0.0001 ){
6789 printf("problem\n");
6796 //__________________________________________________________________________
6797 void AliTPCtrackerMI::CookLabel(AliKalmanTrack *tk, Float_t wrong) const {
6798 //--------------------------------------------------------------------
6799 //This function "cooks" a track label. If label<0, this track is fake.
6800 //--------------------------------------------------------------------
6801 AliTPCseed * t = dynamic_cast<AliTPCseed*>(tk);
6803 printf("%s:%d wrong type \n",(char*)__FILE__,__LINE__);
6807 Int_t noc=t->GetNumberOfClusters();
6809 //printf("\nnot founded prolongation\n\n\n");
6815 AliTPCclusterMI *clusters[160];
6817 for (Int_t i=0;i<160;i++) {
6824 for (i=0; i<160 && current<noc; i++) {
6826 Int_t index=t->GetClusterIndex2(i);
6827 if (index<=0) continue;
6828 if (index&0x8000) continue;
6830 //clusters[current]=GetClusterMI(index);
6831 if (t->GetClusterPointer(i)){
6832 clusters[current]=t->GetClusterPointer(i);
6838 Int_t lab=123456789;
6839 for (i=0; i<noc; i++) {
6840 AliTPCclusterMI *c=clusters[i];
6842 lab=TMath::Abs(c->GetLabel(0));
6844 for (j=0; j<noc; j++) if (lb[j]==lab || mx[j]==0) break;
6850 for (i=0; i<noc; i++) if (mx[i]>max) {max=mx[i]; lab=lb[i];}
6852 for (i=0; i<noc; i++) {
6853 AliTPCclusterMI *c=clusters[i];
6855 if (TMath::Abs(c->GetLabel(1)) == lab ||
6856 TMath::Abs(c->GetLabel(2)) == lab ) max++;
6859 if ((1.- Float_t(max)/noc) > wrong) lab=-lab;
6862 Int_t tail=Int_t(0.10*noc);
6865 for (i=1; i<=160&&ind<tail; i++) {
6866 // AliTPCclusterMI *c=clusters[noc-i];
6867 AliTPCclusterMI *c=clusters[i];
6869 if (lab == TMath::Abs(c->GetLabel(0)) ||
6870 lab == TMath::Abs(c->GetLabel(1)) ||
6871 lab == TMath::Abs(c->GetLabel(2))) max++;
6874 if (max < Int_t(0.5*tail)) lab=-lab;
6881 //delete[] clusters;
6885 //__________________________________________________________________________
6886 Int_t AliTPCtrackerMI::CookLabel(AliTPCseed *const t, Float_t wrong,Int_t first, Int_t last) const {
6887 //--------------------------------------------------------------------
6888 //This function "cooks" a track label. If label<0, this track is fake.
6889 //--------------------------------------------------------------------
6890 Int_t noc=t->GetNumberOfClusters();
6892 //printf("\nnot founded prolongation\n\n\n");
6898 AliTPCclusterMI *clusters[160];
6900 for (Int_t i=0;i<160;i++) {
6907 for (i=0; i<160 && current<noc; i++) {
6908 if (i<first) continue;
6909 if (i>last) continue;
6910 Int_t index=t->GetClusterIndex2(i);
6911 if (index<=0) continue;
6912 if (index&0x8000) continue;
6914 //clusters[current]=GetClusterMI(index);
6915 if (t->GetClusterPointer(i)){
6916 clusters[current]=t->GetClusterPointer(i);
6921 if (noc<5) return -1;
6922 Int_t lab=123456789;
6923 for (i=0; i<noc; i++) {
6924 AliTPCclusterMI *c=clusters[i];
6926 lab=TMath::Abs(c->GetLabel(0));
6928 for (j=0; j<noc; j++) if (lb[j]==lab || mx[j]==0) break;
6934 for (i=0; i<noc; i++) if (mx[i]>max) {max=mx[i]; lab=lb[i];}
6936 for (i=0; i<noc; i++) {
6937 AliTPCclusterMI *c=clusters[i];
6939 if (TMath::Abs(c->GetLabel(1)) == lab ||
6940 TMath::Abs(c->GetLabel(2)) == lab ) max++;
6943 if ((1.- Float_t(max)/noc) > wrong) lab=-lab;
6946 Int_t tail=Int_t(0.10*noc);
6949 for (i=1; i<=160&&ind<tail; i++) {
6950 // AliTPCclusterMI *c=clusters[noc-i];
6951 AliTPCclusterMI *c=clusters[i];
6953 if (lab == TMath::Abs(c->GetLabel(0)) ||
6954 lab == TMath::Abs(c->GetLabel(1)) ||
6955 lab == TMath::Abs(c->GetLabel(2))) max++;
6958 if (max < Int_t(0.5*tail)) lab=-lab;
6961 // t->SetLabel(lab);
6965 //delete[] clusters;
6969 Int_t AliTPCtrackerMI::GetRowNumber(Double_t x[3]) const
6971 //return pad row number for given x vector
6972 Float_t phi = TMath::ATan2(x[1],x[0]);
6973 if(phi<0) phi=2.*TMath::Pi()+phi;
6974 // Get the local angle in the sector philoc
6975 const Float_t kRaddeg = 180/3.14159265358979312;
6976 Float_t phiangle = (Int_t (phi*kRaddeg/20.) + 0.5)*20./kRaddeg;
6977 Double_t localx = x[0]*TMath::Cos(phiangle)-x[1]*TMath::Sin(phiangle);
6978 return GetRowNumber(localx);
6983 void AliTPCtrackerMI::MakeBitmaps(AliTPCseed *t)
6985 //-----------------------------------------------------------------------
6986 // Fill the cluster and sharing bitmaps of the track
6987 //-----------------------------------------------------------------------
6989 Int_t firstpoint = 0;
6990 Int_t lastpoint = 159;
6991 AliTPCTrackerPoint *point;
6992 AliTPCclusterMI *cluster;
6994 for (int iter=firstpoint; iter<lastpoint; iter++) {
6995 // Change to cluster pointers to see if we have a cluster at given padrow
6996 cluster = t->GetClusterPointer(iter);
6998 t->SetClusterMapBit(iter, kTRUE);
6999 point = t->GetTrackPoint(iter);
7000 if (point->IsShared())
7001 t->SetSharedMapBit(iter,kTRUE);
7003 t->SetSharedMapBit(iter, kFALSE);
7006 t->SetClusterMapBit(iter, kFALSE);
7007 t->SetSharedMapBit(iter, kFALSE);
7012 Bool_t AliTPCtrackerMI::IsFindable(AliTPCseed & track){
7014 // return flag if there is findable cluster at given position
7017 Float_t z = track.GetZ();
7019 if (TMath::Abs(z)<(AliTPCReconstructor::GetCtgRange()*track.GetX()+kDeltaZ) &&
7020 TMath::Abs(z)<fkParam->GetZLength(0) &&
7021 (TMath::Abs(track.GetSnp())<AliTPCReconstructor::GetMaxSnpTracker()))
7027 void AliTPCtrackerMI::AddCovariance(AliTPCseed * seed){
7029 // Adding systematic error
7030 // !!!! the systematic error for element 4 is in 1/cm not in pt
7032 const Double_t *param = AliTPCReconstructor::GetRecoParam()->GetSystematicError();
7034 for (Int_t i=0;i<15;i++) covar[i]=0;
7040 covar[0] = param[0]*param[0];
7041 covar[2] = param[1]*param[1];
7042 covar[5] = param[2]*param[2];
7043 covar[9] = param[3]*param[3];
7044 Double_t facC = AliTracker::GetBz()*kB2C;
7045 covar[14]= param[4]*param[4]*facC*facC;
7046 seed->AddCovariance(covar);