1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
17 //-------------------------------------------------------
18 // Implementation of the TPC tracker
20 // Origin: Marian Ivanov Marian.Ivanov@cern.ch
22 // AliTPC parallel tracker
24 // The track fitting is based on Kalaman filtering approach
26 // The track finding steps:
27 // 1. Seeding - with and without vertex constraint
28 // - seeding with vertex constain done at first n^2 proble
29 // - seeding without vertex constraint n^3 problem
30 // 2. Tracking - follow prolongation road - find cluster - update kalman track
32 // The seeding and tracking is repeated several times, in different seeding region.
33 // This approach enables to find the track which cannot be seeded in some region of TPC
34 // This can happen because of low momenta (track do not reach outer radius), or track is currently in the ded region between sectors, or the track is for the moment overlapped with other track (seed quality is poor) ...
36 // With this approach we reach almost 100 % efficiency also for high occupancy events.
37 // (If the seeding efficiency in a region is about 90 % than with logical or of several
38 // regions we will reach 100% (in theory - supposing independence)
40 // Repeating several seeding - tracking procedures some of the tracks can be find
43 // The procedures to remove multi find tacks are impremented:
44 // RemoveUsed2 - fast procedure n problem -
45 // Algorithm - Sorting tracks according quality
46 // remove tracks with some shared fraction
47 // Sharing in respect to all tacks
48 // Signing clusters in gold region
49 // FindSplitted - slower algorithm n^2
50 // Sort the tracks according quality
51 // Loop over pair of tracks
52 // If overlap with other track bigger than threshold - remove track
54 // FindCurling - Finds the pair of tracks which are curling
55 // - About 10% of tracks can be find with this procedure
56 // The combinatorial background is too big to be used in High
57 // multiplicity environment
58 // - n^2 problem - Slow procedure - currently it is disabled because of
61 // The number of splitted tracks can be reduced disabling the sharing of the cluster.
62 // tpcRecoParam-> SetClusterSharing(kFALSE);
63 // IT IS HIGHLY non recomended to use it in high flux enviroonment
64 // Even using this switch some tracks can be found more than once
65 // (because of multiple seeding and low quality tracks which will not cross full chamber)
68 // The tracker itself can be debugged - the information about tracks can be stored in several // phases of the reconstruction
69 // To enable storage of the TPC tracks in the ESD friend track
70 // use AliTPCReconstructor::SetStreamLevel(n); where nis bigger 0
72 // The debug level - different procedure produce tree for numerical debugging
73 // To enable them set AliTPCReconstructor::SetStreamLevel(n); where nis bigger 1
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 (Log level 9) - dump information about
100 // 2.a) cluster error estimate
101 // 3.a) cluster shape estimate
104 //-------------------------------------------------------
109 #include "Riostream.h"
110 #include <TClonesArray.h>
112 #include <TObjArray.h>
115 #include "AliComplexCluster.h"
116 #include "AliESDEvent.h"
117 #include "AliESDtrack.h"
118 #include "AliESDVertex.h"
121 #include "AliHelix.h"
122 #include "AliRunLoader.h"
123 #include "AliTPCClustersRow.h"
124 #include "AliTPCParam.h"
125 #include "AliTPCReconstructor.h"
126 #include "AliTPCpolyTrack.h"
127 #include "AliTPCreco.h"
128 #include "AliTPCseed.h"
130 #include "AliTPCtrackerSector.h"
131 #include "AliTPCtrackerMI.h"
132 #include "TStopwatch.h"
133 #include "AliTPCReconstructor.h"
135 #include "TTreeStream.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()>1 && 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()>0) {
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;
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;
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)
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 *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()>1) {
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()>15){
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()>15){
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()>2) 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()>0) {
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()>2) 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()>0) {
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 *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 *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()>0) {
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()>2){
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()>1){
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);
5011 kinks->AddLast(kink);
5017 // sort the kinks according quality - and refit them towards vertex
5019 Int_t nkinks = kinks->GetEntriesFast();
5020 Float_t *quality = new Float_t[nkinks];
5021 Int_t *indexes = new Int_t[nkinks];
5022 AliTPCseed *mothers = new AliTPCseed[nkinks];
5023 AliTPCseed *daughters = new AliTPCseed[nkinks];
5026 for (Int_t i=0;i<nkinks;i++){
5028 AliKink *kinkl = (AliKink*)kinks->At(i);
5030 // refit kinks towards vertex
5032 Int_t index0 = kinkl->GetIndex(0);
5033 Int_t index1 = kinkl->GetIndex(1);
5034 AliTPCseed * ktrack0 = (AliTPCseed*)array->At(index0);
5035 AliTPCseed * ktrack1 = (AliTPCseed*)array->At(index1);
5037 Int_t sumn=ktrack0->GetNumberOfClusters()+ktrack1->GetNumberOfClusters();
5039 // Refit Kink under if too small angle
5041 if (kinkl->GetAngle(2)<0.05){
5042 kinkl->SetTPCRow0(GetRowNumber(kinkl->GetR()));
5043 Int_t row0 = kinkl->GetTPCRow0();
5044 Int_t drow = Int_t(2.+0.5/(0.05+kinkl->GetAngle(2)));
5047 Int_t last = row0-drow;
5048 if (last<40) last=40;
5049 if (last<ktrack0->GetFirstPoint()+25) last = ktrack0->GetFirstPoint()+25;
5050 AliTPCseed* seed0 = ReSeed(ktrack0,last,kFALSE);
5053 Int_t first = row0+drow;
5054 if (first>130) first=130;
5055 if (first>ktrack1->GetLastPoint()-25) first = TMath::Max(ktrack1->GetLastPoint()-25,30);
5056 AliTPCseed* seed1 = ReSeed(ktrack1,first,kTRUE);
5058 if (seed0 && seed1){
5059 kinkl->SetStatus(1,8);
5060 if (RefitKink(*seed0,*seed1,*kinkl)) kinkl->SetStatus(1,9);
5061 row0 = GetRowNumber(kinkl->GetR());
5062 sumn = seed0->GetNumberOfClusters()+seed1->GetNumberOfClusters();
5063 mothers[i] = *seed0;
5064 daughters[i] = *seed1;
5067 delete kinks->RemoveAt(i);
5068 if (seed0) delete seed0;
5069 if (seed1) delete seed1;
5072 if (kinkl->GetDistance()>0.5 || kinkl->GetR()<110 || kinkl->GetR()>240) {
5073 delete kinks->RemoveAt(i);
5074 if (seed0) delete seed0;
5075 if (seed1) delete seed1;
5083 if (kinkl) quality[i] = 160*((0.1+kinkl->GetDistance())*(2.-kinkl->GetTPCDensityFactor()))/(sumn+40.); //the longest -clossest will win
5085 TMath::Sort(nkinks,quality,indexes,kFALSE);
5087 //remove double find kinks
5089 for (Int_t ikink0=1;ikink0<nkinks;ikink0++){
5090 AliKink * kink0 = (AliKink*) kinks->At(indexes[ikink0]);
5091 if (!kink0) continue;
5093 for (Int_t ikink1=0;ikink1<ikink0;ikink1++){
5094 if (!kink0) continue;
5095 AliKink * kink1 = (AliKink*) kinks->At(indexes[ikink1]);
5096 if (!kink1) continue;
5097 // if not close kink continue
5098 if (TMath::Abs(kink1->GetPosition()[2]-kink0->GetPosition()[2])>10) continue;
5099 if (TMath::Abs(kink1->GetPosition()[1]-kink0->GetPosition()[1])>10) continue;
5100 if (TMath::Abs(kink1->GetPosition()[0]-kink0->GetPosition()[0])>10) continue;
5102 AliTPCseed &mother0 = mothers[indexes[ikink0]];
5103 AliTPCseed &daughter0 = daughters[indexes[ikink0]];
5104 AliTPCseed &mother1 = mothers[indexes[ikink1]];
5105 AliTPCseed &daughter1 = daughters[indexes[ikink1]];
5106 Int_t row0 = (kink0->GetTPCRow0()+kink1->GetTPCRow0())/2;
5115 for (Int_t i=0;i<row0;i++){
5116 if (mother0.GetClusterIndex(i)>0 && mother1.GetClusterIndex(i)>0){
5119 if (mother0.GetClusterIndex(i)==mother1.GetClusterIndex(i)){
5126 for (Int_t i=row0;i<158;i++){
5127 if (daughter0.GetClusterIndex(i)>0 && daughter0.GetClusterIndex(i)>0){
5130 if (mother0.GetClusterIndex(i)==mother1.GetClusterIndex(i)){
5136 Float_t ratio = Float_t(same+1)/Float_t(both+1);
5137 Float_t ratiom = Float_t(samem+1)/Float_t(bothm+1);
5138 Float_t ratiod = Float_t(samed+1)/Float_t(bothd+1);
5139 if (ratio>0.3 && ratiom>0.5 &&ratiod>0.5) {
5140 Int_t sum0 = mother0.GetNumberOfClusters()+daughter0.GetNumberOfClusters();
5141 Int_t sum1 = mother1.GetNumberOfClusters()+daughter1.GetNumberOfClusters();
5143 shared[kink0->GetIndex(0)]= kTRUE;
5144 shared[kink0->GetIndex(1)]= kTRUE;
5145 delete kinks->RemoveAt(indexes[ikink0]);
5148 shared[kink1->GetIndex(0)]= kTRUE;
5149 shared[kink1->GetIndex(1)]= kTRUE;
5150 delete kinks->RemoveAt(indexes[ikink1]);
5157 for (Int_t i=0;i<nkinks;i++){
5158 AliKink * kinkl = (AliKink*) kinks->At(indexes[i]);
5159 if (!kinkl) continue;
5160 kinkl->SetTPCRow0(GetRowNumber(kinkl->GetR()));
5161 Int_t index0 = kinkl->GetIndex(0);
5162 Int_t index1 = kinkl->GetIndex(1);
5163 if (circular[index0]||(circular[index1]&&kinkl->GetDistance()>0.2)) continue;
5164 kinkl->SetMultiple(usage[index0],0);
5165 kinkl->SetMultiple(usage[index1],1);
5166 if (kinkl->GetMultiple()[0]+kinkl->GetMultiple()[1]>2) continue;
5167 if (kinkl->GetMultiple()[0]+kinkl->GetMultiple()[1]>0 && quality[indexes[i]]>0.2) continue;
5168 if (kinkl->GetMultiple()[0]+kinkl->GetMultiple()[1]>0 && kinkl->GetDistance()>0.2) continue;
5169 if (circular[index0]||(circular[index1]&&kinkl->GetDistance()>0.1)) continue;
5171 AliTPCseed * ktrack0 = (AliTPCseed*)array->At(index0);
5172 AliTPCseed * ktrack1 = (AliTPCseed*)array->At(index1);
5173 if (!ktrack0 || !ktrack1) continue;
5174 Int_t index = esd->AddKink(kinkl);
5177 if ( ktrack0->GetKinkIndex(0)==0 && ktrack1->GetKinkIndex(0)==0) { //best kink
5178 if (mothers[indexes[i]].GetNumberOfClusters()>20 && daughters[indexes[i]].GetNumberOfClusters()>20 && (mothers[indexes[i]].GetNumberOfClusters()+daughters[indexes[i]].GetNumberOfClusters())>100){
5179 *ktrack0 = mothers[indexes[i]];
5180 *ktrack1 = daughters[indexes[i]];
5184 ktrack0->SetKinkIndex(usage[index0],-(index+1));
5185 ktrack1->SetKinkIndex(usage[index1], (index+1));
5190 // Remove tracks corresponding to shared kink's
5192 for (Int_t i=0;i<nentries;i++){
5193 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5194 if (!track0) continue;
5195 if (track0->GetKinkIndex(0)!=0) continue;
5196 if (shared[i]) delete array->RemoveAt(i);
5201 RemoveUsed2(array,0.5,0.4,30);
5203 for (Int_t i=0;i<nentries;i++){
5204 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5205 if (!track0) continue;
5206 track0->CookdEdx(0.02,0.6);
5210 for (Int_t i=0;i<nentries;i++){
5211 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5212 if (!track0) continue;
5213 if (track0->Pt()<1.4) continue;
5214 //remove double high momenta tracks - overlapped with kink candidates
5217 for (Int_t icl=track0->GetFirstPoint();icl<track0->GetLastPoint(); icl++){
5218 if (track0->GetClusterPointer(icl)!=0){
5220 if (track0->GetClusterPointer(icl)->IsUsed(10)) ishared++;
5223 if (Float_t(ishared+1)/Float_t(all+1)>0.5) {
5224 delete array->RemoveAt(i);
5228 if (track0->GetKinkIndex(0)!=0) continue;
5229 if (track0->GetNumberOfClusters()<80) continue;
5231 AliTPCseed *pmother = new AliTPCseed();
5232 AliTPCseed *pdaughter = new AliTPCseed();
5233 AliKink *pkink = new AliKink;
5235 AliTPCseed & mother = *pmother;
5236 AliTPCseed & daughter = *pdaughter;
5237 AliKink & kinkl = *pkink;
5238 if (CheckKinkPoint(track0,mother,daughter, kinkl)){
5239 if (mother.GetNumberOfClusters()<30||daughter.GetNumberOfClusters()<20) {
5243 continue; //too short tracks
5245 if (mother.Pt()<1.4) {
5251 Int_t row0= kinkl.GetTPCRow0();
5252 if (kinkl.GetDistance()>0.5 || kinkl.GetR()<110. || kinkl.GetR()>240.) {
5259 Int_t index = esd->AddKink(&kinkl);
5260 mother.SetKinkIndex(0,-(index+1));
5261 daughter.SetKinkIndex(0,index+1);
5262 if (mother.GetNumberOfClusters()>50) {
5263 delete array->RemoveAt(i);
5264 array->AddAt(new AliTPCseed(mother),i);
5267 array->AddLast(new AliTPCseed(mother));
5269 array->AddLast(new AliTPCseed(daughter));
5270 for (Int_t icl=0;icl<row0;icl++) {
5271 if (mother.GetClusterPointer(icl)) mother.GetClusterPointer(icl)->Use(20);
5274 for (Int_t icl=row0;icl<158;icl++) {
5275 if (daughter.GetClusterPointer(icl)) daughter.GetClusterPointer(icl)->Use(20);
5284 delete [] daughters;
5306 printf("Ncandidates=\t%d\t%d\t%d\t%d\n",esd->GetNumberOfKinks(),ncandidates,ntracks,nall);
5310 void AliTPCtrackerMI::FindV0s(const TObjArray * array, AliESDEvent *esd)
5316 TObjArray *tpcv0s = new TObjArray(100000);
5317 Int_t nentries = array->GetEntriesFast();
5318 AliHelix *helixes = new AliHelix[nentries];
5319 Int_t *sign = new Int_t[nentries];
5320 Float_t *alpha = new Float_t[nentries];
5321 Float_t *z0 = new Float_t[nentries];
5322 Float_t *dca = new Float_t[nentries];
5323 Float_t *sdcar = new Float_t[nentries];
5324 Float_t *cdcar = new Float_t[nentries];
5325 Float_t *pulldcar = new Float_t[nentries];
5326 Float_t *pulldcaz = new Float_t[nentries];
5327 Float_t *pulldca = new Float_t[nentries];
5328 Bool_t *isPrim = new Bool_t[nentries];
5329 const AliESDVertex * primvertex = esd->GetVertex();
5330 Double_t zvertex = primvertex->GetZv();
5332 // nentries = array->GetEntriesFast();
5334 for (Int_t i=0;i<nentries;i++){
5337 AliTPCseed* track = (AliTPCseed*)array->At(i);
5338 if (!track) continue;
5339 track->GetV0Indexes()[0] = 0; //rest v0 indexes
5340 track->GetV0Indexes()[1] = 0; //rest v0 indexes
5341 track->GetV0Indexes()[2] = 0; //rest v0 indexes
5343 alpha[i] = track->GetAlpha();
5344 new (&helixes[i]) AliHelix(*track);
5346 helixes[i].Evaluate(0,xyz);
5347 sign[i] = (track->GetC()>0) ? -1:1;
5351 if (track->GetProlongation(0,y,z)) z0[i] = z;
5352 dca[i] = track->GetD(0,0);
5354 // dca error parrameterezation + pulls
5356 sdcar[i] = TMath::Sqrt(0.150*0.150+(100*track->GetC())*(100*track->GetC()));
5357 if (TMath::Abs(track->GetTgl())>1) sdcar[i]*=2.5;
5358 cdcar[i] = TMath::Exp((TMath::Abs(track->GetC())-0.0106)*525.3);
5359 pulldcar[i] = (dca[i]-cdcar[i])/sdcar[i];
5360 pulldcaz[i] = (z0[i]-zvertex)/sdcar[i];
5361 pulldca[i] = TMath::Sqrt(pulldcar[i]*pulldcar[i]+pulldcaz[i]*pulldcaz[i]);
5362 if (track->TPCrPID(1)+track->TPCrPID(2)+track->TPCrPID(3)>0.5) {
5363 if (pulldca[i]<3.) isPrim[i]=kTRUE; //pion, muon and Kaon 3 sigma cut
5365 if (track->TPCrPID(4)>0.5) {
5366 if (pulldca[i]<0.5) isPrim[i]=kTRUE; //proton 0.5 sigma cut
5368 if (track->TPCrPID(0)>0.4) {
5369 isPrim[i]=kFALSE; //electron no sigma cut
5376 Int_t ncandidates =0;
5379 Double_t phase[2][2],radius[2];
5385 Float_t fprimvertex[3]={GetX(),GetY(),GetZ()};
5387 Double_t cradius0 = 10*10;
5388 Double_t cradius1 = 200*200;
5391 Double_t cpointAngle = 0.95;
5393 Double_t delta[2]={10000,10000};
5394 for (Int_t i =0;i<nentries;i++){
5395 if (sign[i]==0) continue;
5396 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5397 if (!track0) continue;
5398 if (AliTPCReconstructor::StreamLevel()>1){
5399 TTreeSRedirector &cstream = *fDebugStreamer;
5404 "zvertex="<<zvertex<<
5405 "sdcar0="<<sdcar[i]<<
5406 "cdcar0="<<cdcar[i]<<
5407 "pulldcar0="<<pulldcar[i]<<
5408 "pulldcaz0="<<pulldcaz[i]<<
5409 "pulldca0="<<pulldca[i]<<
5410 "isPrim="<<isPrim[i]<<
5414 if (track0->GetSigned1Pt()<0) continue;
5415 if (track0->GetKinkIndex(0)>0||isPrim[i]) continue; //daughter kink
5417 if (TMath::Abs(helixes[i].GetHelix(4))<0.000000001) continue;
5422 for (Int_t j =0;j<nentries;j++){
5423 AliTPCseed * track1 = (AliTPCseed*)array->At(j);
5424 if (!track1) continue;
5425 if (track1->GetKinkIndex(0)>0 || isPrim[j]) continue; //daughter kink
5426 if (sign[j]*sign[i]>0) continue;
5427 if (TMath::Abs(helixes[j].GetHelix(4))<0.000001) continue;
5428 if (track0->GetCircular()+track1->GetCircular()>1) continue; //circling -returning track
5431 // DCA to prim vertex cut
5437 Int_t npoints = helixes[i].GetRPHIintersections(helixes[j], phase, radius,cdist2);
5438 if (npoints<1) continue;
5442 if (radius[0]<cradius0||radius[0]>cradius1) continue;
5443 helixes[i].LinearDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta[0]);
5444 if (delta[0]>cdist1) continue;
5447 if (TMath::Max(radius[0],radius[1])<cradius0|| TMath::Min(radius[0],radius[1])>cradius1) continue;
5448 helixes[i].LinearDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta[0]);
5449 helixes[i].LinearDCA(helixes[j],phase[1][0],phase[1][1],radius[1],delta[1]);
5450 if (delta[1]<delta[0]) iclosest=1;
5451 if (delta[iclosest]>cdist1) continue;
5453 helixes[i].ParabolicDCA(helixes[j],phase[iclosest][0],phase[iclosest][1],radius[iclosest],delta[iclosest]);
5454 if (radius[iclosest]<cradius0 || radius[iclosest]>cradius1 || delta[iclosest]>cdist1) continue;
5456 Double_t pointAngle = helixes[i].GetPointAngle(helixes[j],phase[iclosest],fprimvertex);
5457 if (pointAngle<cpointAngle) continue;
5459 Bool_t isGamma = kFALSE;
5460 vertex.SetParamP(*track0); //track0 - plus
5461 vertex.SetParamN(*track1); //track1 - minus
5462 vertex.Update(fprimvertex);
5463 if (track0->TPCrPID(0)>0.3&&track1->TPCrPID(0)>0.3&&vertex.GetAnglep()[2]<0.15) isGamma=kTRUE; // gamma conversion candidate
5464 Double_t pointAngle2 = vertex.GetV0CosineOfPointingAngle();
5466 if (vertex.GetV0CosineOfPointingAngle()<cpointAngle && (!isGamma)) continue;// point angle cut
5467 //Bo: if (vertex.GetDist2()>2&&(!isGamma)) continue; // point angle cut
5468 if (vertex.GetDcaV0Daughters()>2&&(!isGamma)) continue;//Bo: // point angle cut
5469 Float_t sigmae = 0.15*0.15;
5470 if (vertex.GetRr()<80)
5471 sigmae += (sdcar[i]*sdcar[i]+sdcar[j]*sdcar[j])*(1.-vertex.GetRr()/80.)*(1.-vertex.GetRr()/80.);
5472 sigmae+= TMath::Sqrt(sigmae);
5473 //Bo: if (vertex.GetDist2()/sigmae>3.&&(!isGamma)) continue;
5474 if (vertex.GetDcaV0Daughters()/sigmae>3.&&(!isGamma)) continue;
5475 Float_t densb0=0,densb1=0,densa0=0,densa1=0;
5476 Int_t row0 = GetRowNumber(vertex.GetRr());
5478 //Bo: if (vertex.GetDist2()>0.2) continue;
5479 if (vertex.GetDcaV0Daughters()>0.2) continue;
5480 densb0 = track0->Density2(0,row0-5);
5481 densb1 = track1->Density2(0,row0-5);
5482 if (densb0>0.3|| densb1>0.3) continue; //clusters before vertex
5483 densa0 = track0->Density2(row0+5,row0+40);
5484 densa1 = track1->Density2(row0+5,row0+40);
5485 if ((densa0<0.4|| densa1<0.4)&&(!isGamma)) continue; //missing clusters after vertex
5488 densa0 = track0->Density2(0,40); //cluster density
5489 densa1 = track1->Density2(0,40); //cluster density
5490 if ((vertex.GetRr()<80&&densa0+densa1<1.)&&(!isGamma)) continue;
5492 //Bo: vertex.SetLab(0,track0->GetLabel());
5493 //Bo: vertex.SetLab(1,track1->GetLabel());
5494 vertex.SetChi2After((densa0+densa1)*0.5);
5495 vertex.SetChi2Before((densb0+densb1)*0.5);
5496 vertex.SetIndex(0,i);
5497 vertex.SetIndex(1,j);
5498 //Bo: vertex.SetStatus(1); // TPC v0 candidate
5499 vertex.SetOnFlyStatus(2);//Bo: // TPC v0 candidate
5500 //Bo: vertex.SetRp(track0->TPCrPIDs());
5501 //Bo: vertex.SetRm(track1->TPCrPIDs());
5502 tpcv0s->AddLast(new AliESDv0(vertex));
5505 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
5506 Double_t radiusm= (delta[0]<delta[1])? TMath::Sqrt(radius[0]):TMath::Sqrt(radius[1]);
5507 Double_t deltam= (delta[0]<delta[1])? TMath::Sqrt(delta[0]):TMath::Sqrt(delta[1]);
5508 if (AliTPCReconstructor::StreamLevel()>1) {
5509 Int_t lab0=track0->GetLabel();
5510 Int_t lab1=track1->GetLabel();
5511 Char_t c0=track0->GetCircular();
5512 Char_t c1=track1->GetCircular();
5513 TTreeSRedirector &cstream = *fDebugStreamer;
5516 "vertex.="<<&vertex<<
5519 "Helix0.="<<&helixes[i]<<
5522 "Helix1.="<<&helixes[j]<<
5523 "pointAngle="<<pointAngle<<
5524 "pointAngle2="<<pointAngle2<<
5529 "zvertex="<<zvertex<<
5532 "npoints="<<npoints<<
5533 "radius0="<<radius[0]<<
5534 "delta0="<<delta[0]<<
5535 "radius1="<<radius[1]<<
5536 "delta1="<<delta[1]<<
5537 "radiusm="<<radiusm<<
5539 "sdcar0="<<sdcar[i]<<
5540 "sdcar1="<<sdcar[j]<<
5541 "cdcar0="<<cdcar[i]<<
5542 "cdcar1="<<cdcar[j]<<
5543 "pulldcar0="<<pulldcar[i]<<
5544 "pulldcar1="<<pulldcar[j]<<
5545 "pulldcaz0="<<pulldcaz[i]<<
5546 "pulldcaz1="<<pulldcaz[j]<<
5547 "pulldca0="<<pulldca[i]<<
5548 "pulldca1="<<pulldca[j]<<
5558 Float_t *quality = new Float_t[ncandidates];
5559 Int_t *indexes = new Int_t[ncandidates];
5561 for (Int_t i=0;i<ncandidates;i++){
5563 AliESDv0 *v0 = (AliESDv0*)tpcv0s->At(i);
5564 quality[i] = 1./(1.00001-v0->GetV0CosineOfPointingAngle()); //base point angle
5565 // quality[i] /= (0.5+v0->GetDist2());
5566 // quality[i] *= v0->GetChi2After(); //density factor
5568 Int_t index0 = v0->GetIndex(0);
5569 Int_t index1 = v0->GetIndex(1);
5570 //Bo: Double_t minpulldca = TMath::Min(2.+pulldca[v0->GetIndex(0)],(2.+pulldca[v0->GetIndex(1)]) ); //pull
5571 Double_t minpulldca = TMath::Min(2.+pulldca[index0],(2.+pulldca[index1]) );//Bo:
5575 AliTPCseed * track0 = (AliTPCseed*)array->At(index0);
5576 AliTPCseed * track1 = (AliTPCseed*)array->At(index1);
5577 if (track0->TPCrPID(0)>0.3&&track1->TPCrPID(0)>0.3&&v0->GetAnglep()[2]<0.15) quality[i]+=1000000; // gamma conversion candidate
5578 if (track0->TPCrPID(4)>0.9||(track1->TPCrPID(4)>0.9&&minpulldca>4)) quality[i]*=10; // lambda candidate candidate
5581 TMath::Sort(ncandidates,quality,indexes,kTRUE);
5584 for (Int_t i=0;i<ncandidates;i++){
5585 AliESDv0 * v0 = (AliESDv0*)tpcv0s->At(indexes[i]);
5587 Int_t index0 = v0->GetIndex(0);
5588 Int_t index1 = v0->GetIndex(1);
5589 AliTPCseed * track0 = (AliTPCseed*)array->At(index0);
5590 AliTPCseed * track1 = (AliTPCseed*)array->At(index1);
5591 if (!track0||!track1) {
5595 Bool_t accept =kTRUE; //default accept
5596 Int_t *v0indexes0 = track0->GetV0Indexes();
5597 Int_t *v0indexes1 = track1->GetV0Indexes();
5599 Int_t order0 = (v0indexes0[0]!=0) ? 1:0;
5600 Int_t order1 = (v0indexes1[0]!=0) ? 1:0;
5601 if (v0indexes0[1]!=0) order0 =2;
5602 if (v0indexes1[1]!=0) order1 =2;
5604 if (v0indexes0[2]!=0) {order0=3; accept=kFALSE;}
5605 if (v0indexes0[2]!=0) {order1=3; accept=kFALSE;}
5607 AliESDv0 * v02 = v0;
5609 //Bo: v0->SetOrder(0,order0);
5610 //Bo: v0->SetOrder(1,order1);
5611 //Bo: v0->SetOrder(1,order0+order1);
5612 v0->SetOnFlyStatus(kTRUE);
5613 Int_t index = esd->AddV0(v0);
5614 v02 = esd->GetV0(index);
5615 v0indexes0[order0]=index;
5616 v0indexes1[order1]=index;
5620 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
5621 if (AliTPCReconstructor::StreamLevel()>1) {
5622 Int_t lab0=track0->GetLabel();
5623 Int_t lab1=track1->GetLabel();
5624 TTreeSRedirector &cstream = *fDebugStreamer;
5633 "dca0="<<dca[index0]<<
5634 "dca1="<<dca[index1]<<
5638 "quality="<<quality[i]<<
5639 "pulldca0="<<pulldca[index0]<<
5640 "pulldca1="<<pulldca[index1]<<
5664 printf("TPC V0 finder : naccepted\t%d\tncandidates\t%d\tntracks\t%d\tnall\t%d\n",naccepted,ncandidates,ntracks,nall);
5668 Int_t AliTPCtrackerMI::RefitKink(AliTPCseed &mother, AliTPCseed &daughter, AliESDkink &knk)
5671 // refit kink towards to the vertex
5674 AliKink &kink=(AliKink &)knk;
5676 Int_t row0 = GetRowNumber(kink.GetR());
5677 FollowProlongation(mother,0);
5678 mother.Reset(kFALSE);
5680 FollowProlongation(daughter,row0);
5681 daughter.Reset(kFALSE);
5682 FollowBackProlongation(daughter,158);
5683 daughter.Reset(kFALSE);
5684 Int_t first = TMath::Max(row0-20,30);
5685 Int_t last = TMath::Min(row0+20,140);
5687 const Int_t kNdiv =5;
5688 AliTPCseed param0[kNdiv]; // parameters along the track
5689 AliTPCseed param1[kNdiv]; // parameters along the track
5690 AliKink kinks[kNdiv]; // corresponding kink parameters
5693 for (Int_t irow=0; irow<kNdiv;irow++){
5694 rows[irow] = first +((last-first)*irow)/(kNdiv-1);
5696 // store parameters along the track
5698 for (Int_t irow=0;irow<kNdiv;irow++){
5699 FollowBackProlongation(mother, rows[irow]);
5700 FollowProlongation(daughter,rows[kNdiv-1-irow]);
5701 param0[irow] = mother;
5702 param1[kNdiv-1-irow] = daughter;
5706 for (Int_t irow=0; irow<kNdiv-1;irow++){
5707 if (param0[irow].GetNumberOfClusters()<kNdiv||param1[irow].GetNumberOfClusters()<kNdiv) continue;
5708 kinks[irow].SetMother(param0[irow]);
5709 kinks[irow].SetDaughter(param1[irow]);
5710 kinks[irow].Update();
5713 // choose kink with best "quality"
5715 Double_t mindist = 10000;
5716 for (Int_t irow=0;irow<kNdiv;irow++){
5717 if (param0[irow].GetNumberOfClusters()<20||param1[irow].GetNumberOfClusters()<20) continue;
5718 if (TMath::Abs(kinks[irow].GetR())>240.) continue;
5719 if (TMath::Abs(kinks[irow].GetR())<100.) continue;
5721 Float_t normdist = TMath::Abs(param0[irow].GetX()-kinks[irow].GetR())*(0.1+kink.GetDistance());
5722 normdist/= (param0[irow].GetNumberOfClusters()+param1[irow].GetNumberOfClusters()+40.);
5723 if (normdist < mindist){
5729 if (index==-1) return 0;
5732 param0[index].Reset(kTRUE);
5733 FollowProlongation(param0[index],0);
5735 mother = param0[index];
5736 daughter = param1[index]; // daughter in vertex
5738 kink.SetMother(mother);
5739 kink.SetDaughter(daughter);
5741 kink.SetTPCRow0(GetRowNumber(kink.GetR()));
5742 kink.SetTPCncls(param0[index].GetNumberOfClusters(),0);
5743 kink.SetTPCncls(param1[index].GetNumberOfClusters(),1);
5744 kink.SetLabel(CookLabel(&mother,0.4, 0,kink.GetTPCRow0()),0);
5745 kink.SetLabel(CookLabel(&daughter,0.4, kink.GetTPCRow0(),160),1);
5746 mother.SetLabel(kink.GetLabel(0));
5747 daughter.SetLabel(kink.GetLabel(1));
5753 void AliTPCtrackerMI::UpdateKinkQualityM(AliTPCseed * seed){
5755 // update Kink quality information for mother after back propagation
5757 if (seed->GetKinkIndex(0)>=0) return;
5758 for (Int_t ikink=0;ikink<3;ikink++){
5759 Int_t index = seed->GetKinkIndex(ikink);
5760 if (index>=0) break;
5761 index = TMath::Abs(index)-1;
5762 AliESDkink * kink = fEvent->GetKink(index);
5763 kink->SetTPCDensity(-1,0,0);
5764 kink->SetTPCDensity(1,0,1);
5766 Int_t row0 = kink->GetTPCRow0() - 2 - Int_t( 0.5/ (0.05+kink->GetAngle(2)));
5767 if (row0<15) row0=15;
5769 Int_t row1 = kink->GetTPCRow0() + 2 + Int_t( 0.5/ (0.05+kink->GetAngle(2)));
5770 if (row1>145) row1=145;
5772 Int_t found,foundable,shared;
5773 seed->GetClusterStatistic(0,row0, found, foundable,shared,kFALSE);
5774 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),0,0);
5775 seed->GetClusterStatistic(row1,155, found, foundable,shared,kFALSE);
5776 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),0,1);
5781 void AliTPCtrackerMI::UpdateKinkQualityD(AliTPCseed * seed){
5783 // update Kink quality information for daughter after refit
5785 if (seed->GetKinkIndex(0)<=0) return;
5786 for (Int_t ikink=0;ikink<3;ikink++){
5787 Int_t index = seed->GetKinkIndex(ikink);
5788 if (index<=0) break;
5789 index = TMath::Abs(index)-1;
5790 AliESDkink * kink = fEvent->GetKink(index);
5791 kink->SetTPCDensity(-1,1,0);
5792 kink->SetTPCDensity(-1,1,1);
5794 Int_t row0 = kink->GetTPCRow0() -2 - Int_t( 0.5/ (0.05+kink->GetAngle(2)));
5795 if (row0<15) row0=15;
5797 Int_t row1 = kink->GetTPCRow0() +2 + Int_t( 0.5/ (0.05+kink->GetAngle(2)));
5798 if (row1>145) row1=145;
5800 Int_t found,foundable,shared;
5801 seed->GetClusterStatistic(0,row0, found, foundable,shared,kFALSE);
5802 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),1,0);
5803 seed->GetClusterStatistic(row1,155, found, foundable,shared,kFALSE);
5804 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),1,1);
5810 Int_t AliTPCtrackerMI::CheckKinkPoint(AliTPCseed*seed,AliTPCseed &mother, AliTPCseed &daughter, AliESDkink &knk)
5813 // check kink point for given track
5814 // if return value=0 kink point not found
5815 // otherwise seed0 correspond to mother particle
5816 // seed1 correspond to daughter particle
5817 // kink parameter of kink point
5818 AliKink &kink=(AliKink &)knk;
5820 Int_t middlerow = (seed->GetFirstPoint()+seed->GetLastPoint())/2;
5821 Int_t first = seed->GetFirstPoint();
5822 Int_t last = seed->GetLastPoint();
5823 if (last-first<20) return 0; // shortest length - 2*30 = 60 pad-rows
5826 AliTPCseed *seed1 = ReSeed(seed,middlerow+20, kTRUE); //middle of chamber
5827 if (!seed1) return 0;
5828 FollowProlongation(*seed1,seed->GetLastPoint()-20);
5829 seed1->Reset(kTRUE);
5830 FollowProlongation(*seed1,158);
5831 seed1->Reset(kTRUE);
5832 last = seed1->GetLastPoint();
5834 AliTPCseed *seed0 = new AliTPCseed(*seed);
5835 seed0->Reset(kFALSE);
5838 AliTPCseed param0[20]; // parameters along the track
5839 AliTPCseed param1[20]; // parameters along the track
5840 AliKink kinks[20]; // corresponding kink parameters
5842 for (Int_t irow=0; irow<20;irow++){
5843 rows[irow] = first +((last-first)*irow)/19;
5845 // store parameters along the track
5847 for (Int_t irow=0;irow<20;irow++){
5848 FollowBackProlongation(*seed0, rows[irow]);
5849 FollowProlongation(*seed1,rows[19-irow]);
5850 param0[irow] = *seed0;
5851 param1[19-irow] = *seed1;
5855 for (Int_t irow=0; irow<19;irow++){
5856 kinks[irow].SetMother(param0[irow]);
5857 kinks[irow].SetDaughter(param1[irow]);
5858 kinks[irow].Update();
5861 // choose kink with biggest change of angle
5863 Double_t maxchange= 0;
5864 for (Int_t irow=1;irow<19;irow++){
5865 if (TMath::Abs(kinks[irow].GetR())>240.) continue;
5866 if (TMath::Abs(kinks[irow].GetR())<110.) continue;
5867 Float_t quality = TMath::Abs(kinks[irow].GetAngle(2))/(3.+TMath::Abs(kinks[irow].GetR()-param0[irow].GetX()));
5868 if ( quality > maxchange){
5869 maxchange = quality;
5876 if (index<0) return 0;
5878 Int_t row0 = GetRowNumber(kinks[index].GetR()); //row 0 estimate
5879 seed0 = new AliTPCseed(param0[index]);
5880 seed1 = new AliTPCseed(param1[index]);
5881 seed0->Reset(kFALSE);
5882 seed1->Reset(kFALSE);
5883 seed0->ResetCovariance(10.);
5884 seed1->ResetCovariance(10.);
5885 FollowProlongation(*seed0,0);
5886 FollowBackProlongation(*seed1,158);
5887 mother = *seed0; // backup mother at position 0
5888 seed0->Reset(kFALSE);
5889 seed1->Reset(kFALSE);
5890 seed0->ResetCovariance(10.);
5891 seed1->ResetCovariance(10.);
5893 first = TMath::Max(row0-20,0);
5894 last = TMath::Min(row0+20,158);
5896 for (Int_t irow=0; irow<20;irow++){
5897 rows[irow] = first +((last-first)*irow)/19;
5899 // store parameters along the track
5901 for (Int_t irow=0;irow<20;irow++){
5902 FollowBackProlongation(*seed0, rows[irow]);
5903 FollowProlongation(*seed1,rows[19-irow]);
5904 param0[irow] = *seed0;
5905 param1[19-irow] = *seed1;
5909 for (Int_t irow=0; irow<19;irow++){
5910 kinks[irow].SetMother(param0[irow]);
5911 kinks[irow].SetDaughter(param1[irow]);
5912 // param0[irow].Dump();
5913 //param1[irow].Dump();
5914 kinks[irow].Update();
5917 // choose kink with biggest change of angle
5920 for (Int_t irow=0;irow<20;irow++){
5921 if (TMath::Abs(kinks[irow].GetR())>250.) continue;
5922 if (TMath::Abs(kinks[irow].GetR())<90.) continue;
5923 Float_t quality = TMath::Abs(kinks[irow].GetAngle(2))/(3.+TMath::Abs(kinks[irow].GetR()-param0[irow].GetX()));
5924 if ( quality > maxchange){
5925 maxchange = quality;
5932 if (index==-1 || param0[index].GetNumberOfClusters()+param1[index].GetNumberOfClusters()<100){
5937 // Float_t anglesigma = TMath::Sqrt(param0[index].fC22+param0[index].fC33+param1[index].fC22+param1[index].fC33);
5939 kink.SetMother(param0[index]);
5940 kink.SetDaughter(param1[index]);
5942 row0 = GetRowNumber(kink.GetR());
5943 kink.SetTPCRow0(row0);
5944 kink.SetLabel(CookLabel(seed0,0.5,0,row0),0);
5945 kink.SetLabel(CookLabel(seed1,0.5,row0,158),1);
5946 kink.SetIndex(-10,0);
5947 kink.SetIndex(int(param0[index].GetNumberOfClusters()+param1[index].GetNumberOfClusters()),1);
5948 kink.SetTPCncls(param0[index].GetNumberOfClusters(),0);
5949 kink.SetTPCncls(param1[index].GetNumberOfClusters(),1);
5952 // new (&mother) AliTPCseed(param0[index]);
5953 daughter = param1[index];
5954 daughter.SetLabel(kink.GetLabel(1));
5955 param0[index].Reset(kTRUE);
5956 FollowProlongation(param0[index],0);
5957 mother = param0[index];
5958 mother.SetLabel(kink.GetLabel(0));
5968 AliTPCseed* AliTPCtrackerMI::ReSeed(AliTPCseed *t)
5971 // reseed - refit - track
5974 // Int_t last = fSectors->GetNRows()-1;
5976 if (fSectors == fOuterSec){
5977 first = TMath::Max(first, t->GetFirstPoint()-fInnerSec->GetNRows());
5981 first = t->GetFirstPoint();
5983 AliTPCseed * seed = MakeSeed(t,0.1,0.5,0.9);
5984 FollowBackProlongation(*t,fSectors->GetNRows()-1);
5986 FollowProlongation(*t,first);
5996 //_____________________________________________________________________________
5997 Int_t AliTPCtrackerMI::ReadSeeds(const TFile *inp) {
5998 //-----------------------------------------------------------------
5999 // This function reades track seeds.
6000 //-----------------------------------------------------------------
6001 TDirectory *savedir=gDirectory;
6003 TFile *in=(TFile*)inp;
6004 if (!in->IsOpen()) {
6005 cerr<<"AliTPCtrackerMI::ReadSeeds(): input file is not open !\n";
6010 TTree *seedTree=(TTree*)in->Get("Seeds");
6012 cerr<<"AliTPCtrackerMI::ReadSeeds(): ";
6013 cerr<<"can't get a tree with track seeds !\n";
6016 AliTPCtrack *seed=new AliTPCtrack;
6017 seedTree->SetBranchAddress("tracks",&seed);
6019 if (fSeeds==0) fSeeds=new TObjArray(15000);
6021 Int_t n=(Int_t)seedTree->GetEntries();
6022 for (Int_t i=0; i<n; i++) {
6023 seedTree->GetEvent(i);
6024 fSeeds->AddLast(new AliTPCseed(*seed/*,seed->GetAlpha()*/));
6033 Int_t AliTPCtrackerMI::Clusters2Tracks (AliESDEvent *esd)
6036 if (fSeeds) DeleteSeeds();
6039 if (!fSeeds) return 1;
6046 //_____________________________________________________________________________
6047 Int_t AliTPCtrackerMI::Clusters2Tracks() {
6048 //-----------------------------------------------------------------
6049 // This is a track finder.
6050 //-----------------------------------------------------------------
6051 TDirectory *savedir=gDirectory;
6055 fSeeds = Tracking();
6058 Info("Clusters2Tracks","Time for tracking: \t");timer.Print();timer.Start();
6060 //activate again some tracks
6061 for (Int_t i=0; i<fSeeds->GetEntriesFast(); i++) {
6062 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
6064 Int_t nc=t.GetNumberOfClusters();
6066 delete fSeeds->RemoveAt(i);
6070 if (pt->GetRemoval()==10) {
6071 if (pt->GetDensityFirst(20)>0.8 || pt->GetDensityFirst(30)>0.8 || pt->GetDensityFirst(40)>0.7)
6072 pt->Desactivate(10); // make track again active
6074 pt->Desactivate(20);
6075 delete fSeeds->RemoveAt(i);
6080 RemoveUsed2(fSeeds,0.85,0.85,0);
6081 if (AliTPCReconstructor::GetRecoParam()->GetDoKinks()) FindKinks(fSeeds,fEvent);
6082 //FindCurling(fSeeds, fEvent,0);
6083 if (AliTPCReconstructor::StreamLevel()>2) FindMultiMC(fSeeds, fEvent,-1); // find multi found tracks
6084 RemoveUsed2(fSeeds,0.5,0.4,20);
6085 FindSplitted(fSeeds, fEvent,0); // find multi found tracks
6086 if (AliTPCReconstructor::StreamLevel()>2) FindMultiMC(fSeeds, fEvent,0); // find multi found tracks
6089 // // refit short tracks
6091 Int_t nseed=fSeeds->GetEntriesFast();
6094 for (Int_t i=0; i<nseed; i++) {
6095 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
6097 Int_t nc=t.GetNumberOfClusters();
6099 delete fSeeds->RemoveAt(i);
6102 CookLabel(pt,0.1); //For comparison only
6103 //if ((pt->IsActive() || (pt->fRemoval==10) )&& nc>50 &&pt->GetNumberOfClusters()>0.4*pt->fNFoundable){
6104 if ((pt->IsActive() || (pt->GetRemoval()==10) )){
6106 if (fDebug>0) cerr<<found<<'\r';
6110 delete fSeeds->RemoveAt(i);
6114 //RemoveOverlap(fSeeds,0.99,7,kTRUE);
6116 //RemoveUsed(fSeeds,0.9,0.9,6);
6118 nseed=fSeeds->GetEntriesFast();
6120 for (Int_t i=0; i<nseed; i++) {
6121 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
6123 Int_t nc=t.GetNumberOfClusters();
6125 delete fSeeds->RemoveAt(i);
6129 t.CookdEdx(0.02,0.6);
6130 // CheckKinkPoint(&t,0.05);
6131 //if ((pt->IsActive() || (pt->fRemoval==10) )&& nc>50 &&pt->GetNumberOfClusters()>0.4*pt->fNFoundable){
6132 if ((pt->IsActive() || (pt->GetRemoval()==10) )){
6140 delete fSeeds->RemoveAt(i);
6141 //AliTPCseed * seed1 = ReSeed(pt,0.05,0.5,1);
6143 // FollowProlongation(*seed1,0);
6144 // Int_t n = seed1->GetNumberOfClusters();
6145 // printf("fP4\t%f\t%f\n",seed1->GetC(),pt->GetC());
6146 // printf("fN\t%d\t%d\n", seed1->GetNumberOfClusters(),pt->GetNumberOfClusters());
6149 //AliTPCseed * seed2 = ReSeed(pt,0.95,0.5,0.05);
6153 SortTracks(fSeeds, 1);
6157 PrepareForBackProlongation(fSeeds,5.);
6158 PropagateBack(fSeeds);
6159 printf("Time for back propagation: \t");timer.Print();timer.Start();
6163 PrepareForProlongation(fSeeds,5.);
6164 PropagateForward2(fSeeds);
6166 printf("Time for FORWARD propagation: \t");timer.Print();timer.Start();
6167 // RemoveUsed(fSeeds,0.7,0.7,6);
6168 //RemoveOverlap(fSeeds,0.9,7,kTRUE);
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);
6180 t.CookdEdx(0.02,0.6);
6181 // CookLabel(pt,0.1); //For comparison only
6182 //if ((pt->IsActive() || (pt->fRemoval==10) )&& nc>50 &&pt->GetNumberOfClusters()>0.4*pt->fNFoundable){
6183 if ((pt->IsActive() || (pt->fRemoval==10) )){
6184 cerr<<found++<<'\r';
6187 delete fSeeds->RemoveAt(i);
6192 // fNTracks = found;
6194 Info("Clusters2Tracks","Time for overlap removal, track writing and dedx cooking: \t"); timer.Print();timer.Start();
6197 // cerr<<"Number of found tracks : "<<"\t"<<found<<endl;
6198 Info("Clusters2Tracks","Number of found tracks %d",found);
6200 // UnloadClusters();
6205 void AliTPCtrackerMI::Tracking(TObjArray * arr)
6208 // tracking of the seeds
6211 fSectors = fOuterSec;
6212 ParallelTracking(arr,150,63);
6213 fSectors = fOuterSec;
6214 ParallelTracking(arr,63,0);
6217 TObjArray * AliTPCtrackerMI::Tracking(Int_t seedtype, Int_t i1, Int_t i2, Float_t cuts[4], Float_t dy, Int_t dsec)
6222 TObjArray * arr = new TObjArray;
6224 fSectors = fOuterSec;
6227 for (Int_t sec=0;sec<fkNOS;sec++){
6228 if (seedtype==3) MakeSeeds3(arr,sec,i1,i2,cuts,dy, dsec);
6229 if (seedtype==4) MakeSeeds5(arr,sec,i1,i2,cuts,dy);
6230 if (seedtype==2) MakeSeeds2(arr,sec,i1,i2,cuts,dy);
6233 Info("Tracking","\nSeeding - %d\t%d\t%d\t%d\n",seedtype,i1,i2,arr->GetEntriesFast());
6245 TObjArray * AliTPCtrackerMI::Tracking()
6249 if (AliTPCReconstructor::GetRecoParam()->GetSpecialSeeding()) return TrackingSpecial();
6252 Int_t nup=fOuterSec->GetNRows()+fInnerSec->GetNRows();
6254 TObjArray * seeds = new TObjArray;
6263 Float_t fnumber = 3.0;
6264 Float_t fdensity = 3.0;
6269 for (Int_t delta = 0; delta<18; delta+=6){
6273 arr = Tracking(3,nup-1-delta,nup-1-delta-gap,cuts,-1,1);
6274 SumTracks(seeds,arr);
6275 SignClusters(seeds,fnumber,fdensity);
6277 for (Int_t i=2;i<6;i+=2){
6278 // seed high pt tracks
6281 arr = Tracking(3,nup-i-delta,nup-i-delta-gap,cuts,-1,0);
6282 SumTracks(seeds,arr);
6283 SignClusters(seeds,fnumber,fdensity);
6288 // RemoveUsed(seeds,0.9,0.9,1);
6289 // UnsignClusters();
6290 // SignClusters(seeds,fnumber,fdensity);
6294 for (Int_t delta = 20; delta<120; delta+=10){
6296 // seed high pt tracks
6300 arr = Tracking(3,nup-delta,nup-delta-gap,cuts,-1);
6301 SumTracks(seeds,arr);
6302 SignClusters(seeds,fnumber,fdensity);
6307 arr = Tracking(3,nup-delta-5,nup-delta-5-gap,cuts,-1);
6308 SumTracks(seeds,arr);
6309 SignClusters(seeds,fnumber,fdensity);
6320 Info("Tracking()","\n\nPrimary seeding\t%d\n\n",seeds->GetEntriesFast());
6324 // RemoveUsed(seeds,0.75,0.75,1);
6326 //SignClusters(seeds,fnumber,fdensity);
6335 arr = Tracking(4,nup-1,nup-1-gap,cuts,-1);
6336 SumTracks(seeds,arr);
6337 SignClusters(seeds,fnumber,fdensity);
6339 arr = Tracking(4,nup-2,nup-2-gap,cuts,-1);
6340 SumTracks(seeds,arr);
6341 SignClusters(seeds,fnumber,fdensity);
6343 arr = Tracking(4,nup-3,nup-3-gap,cuts,-1);
6344 SumTracks(seeds,arr);
6345 SignClusters(seeds,fnumber,fdensity);
6349 for (Int_t delta = 3; delta<30; delta+=5){
6355 arr = Tracking(4,nup-1-delta,nup-1-delta-gap,cuts,-1);
6356 SumTracks(seeds,arr);
6357 SignClusters(seeds,fnumber,fdensity);
6359 arr = Tracking(4,nup-3-delta,nup-5-delta-gap,cuts,4);
6360 SumTracks(seeds,arr);
6361 SignClusters(seeds,fnumber,fdensity);
6372 Int_t fLastSeedRowSec=AliTPCReconstructor::GetRecoParam()->GetLastSeedRowSec();
6375 for (Int_t delta = 30; delta<fLastSeedRowSec; delta+=10){
6381 arr = Tracking(4,nup-1-delta,nup-1-delta-gap,cuts,-1);
6382 SumTracks(seeds,arr);
6383 SignClusters(seeds,fnumber,fdensity);
6385 arr = Tracking(4,nup-5-delta,nup-5-delta-gap,cuts,5 );
6386 SumTracks(seeds,arr);
6387 SignClusters(seeds,fnumber,fdensity);
6391 Info("Tracking()","\n\nSecondary seeding\t%d\n\n",seeds->GetEntriesFast());
6402 TObjArray * AliTPCtrackerMI::TrackingSpecial()
6405 // seeding adjusted for laser and cosmic tests - short tracks with big inclination angle
6406 // no primary vertex seeding tried
6410 Int_t nup=fOuterSec->GetNRows()+fInnerSec->GetNRows();
6412 TObjArray * seeds = new TObjArray;
6417 Float_t fnumber = 3.0;
6418 Float_t fdensity = 3.0;
6421 cuts[0] = AliTPCReconstructor::GetRecoParam()->GetMaxC(); // max curvature
6422 cuts[1] = 3.5; // max tan(phi) angle for seeding
6423 cuts[2] = 3.; // not used (cut on z primary vertex)
6424 cuts[3] = 3.5; // max tan(theta) angle for seeding
6426 for (Int_t delta = 0; nup-delta-gap-1>0; delta+=3){
6428 arr = Tracking(4,nup-1-delta,nup-1-delta-gap,cuts,-1);
6429 SumTracks(seeds,arr);
6430 SignClusters(seeds,fnumber,fdensity);
6434 Info("Tracking()","\n\nSecondary seeding\t%d\n\n",seeds->GetEntriesFast());
6445 void AliTPCtrackerMI::SumTracks(TObjArray *arr1,TObjArray *&arr2) const
6448 //sum tracks to common container
6449 //remove suspicious tracks
6450 Int_t nseed = arr2->GetEntriesFast();
6451 for (Int_t i=0;i<nseed;i++){
6452 AliTPCseed *pt=(AliTPCseed*)arr2->UncheckedAt(i);
6455 // remove tracks with too big curvature
6457 if (TMath::Abs(pt->GetC())>AliTPCReconstructor::GetRecoParam()->GetMaxC()){
6458 delete arr2->RemoveAt(i);
6461 // REMOVE VERY SHORT TRACKS
6462 if (pt->GetNumberOfClusters()<20){
6463 delete arr2->RemoveAt(i);
6466 // NORMAL ACTIVE TRACK
6467 if (pt->IsActive()){
6468 arr1->AddLast(arr2->RemoveAt(i));
6471 //remove not usable tracks
6472 if (pt->GetRemoval()!=10){
6473 delete arr2->RemoveAt(i);
6477 // ENABLE ONLY ENOUGH GOOD STOPPED TRACKS
6478 if (pt->GetDensityFirst(20)>0.8 || pt->GetDensityFirst(30)>0.8 || pt->GetDensityFirst(40)>0.7)
6479 arr1->AddLast(arr2->RemoveAt(i));
6481 delete arr2->RemoveAt(i);
6485 delete arr2; arr2 = 0;
6490 void AliTPCtrackerMI::ParallelTracking(TObjArray * arr, Int_t rfirst, Int_t rlast)
6493 // try to track in parralel
6495 Int_t nseed=arr->GetEntriesFast();
6496 //prepare seeds for tracking
6497 for (Int_t i=0; i<nseed; i++) {
6498 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i), &t=*pt;
6500 if (!t.IsActive()) continue;
6501 // follow prolongation to the first layer
6502 if ( (fSectors ==fInnerSec) || (t.GetFirstPoint()-fkParam->GetNRowLow()>rfirst+1) )
6503 FollowProlongation(t, rfirst+1);
6508 for (Int_t nr=rfirst; nr>=rlast; nr--){
6509 if (nr<fInnerSec->GetNRows())
6510 fSectors = fInnerSec;
6512 fSectors = fOuterSec;
6513 // make indexes with the cluster tracks for given
6515 // find nearest cluster
6516 for (Int_t i=0; i<nseed; i++) {
6517 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i), &t=*pt;
6519 if (nr==80) pt->UpdateReference();
6520 if (!pt->IsActive()) continue;
6521 // if ( (fSectors ==fOuterSec) && (pt->fFirstPoint-fkParam->GetNRowLow())<nr) continue;
6522 if (pt->GetRelativeSector()>17) {
6525 UpdateClusters(t,nr);
6527 // prolonagate to the nearest cluster - if founded
6528 for (Int_t i=0; i<nseed; i++) {
6529 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
6531 if (!pt->IsActive()) continue;
6532 // if ((fSectors ==fOuterSec) && (pt->fFirstPoint-fkParam->GetNRowLow())<nr) continue;
6533 if (pt->GetRelativeSector()>17) {
6536 FollowToNextCluster(*pt,nr);
6541 void AliTPCtrackerMI::PrepareForBackProlongation(TObjArray * arr,Float_t fac) const
6545 // if we use TPC track itself we have to "update" covariance
6547 Int_t nseed= arr->GetEntriesFast();
6548 for (Int_t i=0;i<nseed;i++){
6549 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
6553 //rotate to current local system at first accepted point
6554 Int_t index = pt->GetClusterIndex2(pt->GetFirstPoint());
6555 Int_t sec = (index&0xff000000)>>24;
6557 Float_t angle1 = fInnerSec->GetAlpha()*sec+fInnerSec->GetAlphaShift();
6558 if (angle1>TMath::Pi())
6559 angle1-=2.*TMath::Pi();
6560 Float_t angle2 = pt->GetAlpha();
6562 if (TMath::Abs(angle1-angle2)>0.001){
6563 pt->Rotate(angle1-angle2);
6564 //angle2 = pt->GetAlpha();
6565 //pt->fRelativeSector = pt->GetAlpha()/fInnerSec->GetAlpha();
6566 //if (pt->GetAlpha()<0)
6567 // pt->fRelativeSector+=18;
6568 //sec = pt->fRelativeSector;
6577 void AliTPCtrackerMI::PrepareForProlongation(TObjArray * arr, Float_t fac) const
6581 // if we use TPC track itself we have to "update" covariance
6583 Int_t nseed= arr->GetEntriesFast();
6584 for (Int_t i=0;i<nseed;i++){
6585 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
6588 pt->SetFirstPoint(pt->GetLastPoint());
6596 Int_t AliTPCtrackerMI::PropagateBack(TObjArray * arr)
6599 // make back propagation
6601 Int_t nseed= arr->GetEntriesFast();
6602 for (Int_t i=0;i<nseed;i++){
6603 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
6604 if (pt&& pt->GetKinkIndex(0)<=0) {
6605 //AliTPCseed *pt2 = new AliTPCseed(*pt);
6606 fSectors = fInnerSec;
6607 //FollowBackProlongation(*pt,fInnerSec->GetNRows()-1);
6608 //fSectors = fOuterSec;
6609 FollowBackProlongation(*pt,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1);
6610 //if (pt->GetNumberOfClusters()<(pt->fEsd->GetTPCclusters(0)) ){
6611 // Error("PropagateBack","Not prolonged track %d",pt->GetLabel());
6612 // FollowBackProlongation(*pt2,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1);
6615 if (pt&& pt->GetKinkIndex(0)>0) {
6616 AliESDkink * kink = fEvent->GetKink(pt->GetKinkIndex(0)-1);
6617 pt->SetFirstPoint(kink->GetTPCRow0());
6618 fSectors = fInnerSec;
6619 FollowBackProlongation(*pt,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1);
6627 Int_t AliTPCtrackerMI::PropagateForward2(TObjArray * arr)
6630 // make forward propagation
6632 Int_t nseed= arr->GetEntriesFast();
6634 for (Int_t i=0;i<nseed;i++){
6635 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
6637 FollowProlongation(*pt,0);
6646 Int_t AliTPCtrackerMI::PropagateForward()
6649 // propagate track forward
6651 Int_t nseed = fSeeds->GetEntriesFast();
6652 for (Int_t i=0;i<nseed;i++){
6653 AliTPCseed *pt = (AliTPCseed*)fSeeds->UncheckedAt(i);
6655 AliTPCseed &t = *pt;
6656 Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
6657 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
6658 if (alpha < 0. ) alpha += 2.*TMath::Pi();
6659 t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN);
6663 fSectors = fOuterSec;
6664 ParallelTracking(fSeeds,fOuterSec->GetNRows()+fInnerSec->GetNRows()-1,fInnerSec->GetNRows());
6665 fSectors = fInnerSec;
6666 ParallelTracking(fSeeds,fInnerSec->GetNRows()-1,0);
6675 Int_t AliTPCtrackerMI::PropagateBack(AliTPCseed * pt, Int_t row0, Int_t row1)
6678 // make back propagation, in between row0 and row1
6682 fSectors = fInnerSec;
6685 if (row1<fSectors->GetNRows())
6688 r1 = fSectors->GetNRows()-1;
6690 if (row0<fSectors->GetNRows()&& r1>0 )
6691 FollowBackProlongation(*pt,r1);
6692 if (row1<=fSectors->GetNRows())
6695 r1 = row1 - fSectors->GetNRows();
6696 if (r1<=0) return 0;
6697 if (r1>=fOuterSec->GetNRows()) return 0;
6698 fSectors = fOuterSec;
6699 return FollowBackProlongation(*pt,r1);
6707 void AliTPCtrackerMI::GetShape(AliTPCseed * seed, Int_t row)
6711 AliTPCClusterParam * clparam = AliTPCcalibDB::Instance()->GetClusterParam();
6712 Float_t zdrift = TMath::Abs((fkParam->GetZLength(0)-TMath::Abs(seed->GetZ())));
6713 Int_t type = (seed->GetSector() < fkParam->GetNSector()/2) ? 0: (row>126) ? 1:2;
6714 Double_t angulary = seed->GetSnp();
6715 angulary = angulary*angulary/((1.-angulary)*(1.+angulary));
6716 Double_t angularz = seed->GetTgl()*seed->GetTgl()*(1.+angulary);
6718 Double_t sigmay = clparam->GetRMS0(0,type,zdrift,TMath::Sqrt(TMath::Abs(angulary)));
6719 Double_t sigmaz = clparam->GetRMS0(1,type,zdrift,TMath::Sqrt(TMath::Abs(angularz)));
6720 seed->SetCurrentSigmaY2(sigmay*sigmay);
6721 seed->SetCurrentSigmaZ2(sigmaz*sigmaz);
6722 // Float_t sd2 = TMath::Abs((fkParam->GetZLength(0)-TMath::Abs(seed->GetZ())))*fkParam->GetDiffL()*fkParam->GetDiffL();
6723 // // Float_t padlength = fkParam->GetPadPitchLength(seed->fSector);
6724 // Float_t padlength = GetPadPitchLength(row);
6726 // Float_t sresy = (seed->GetSector() < fkParam->GetNSector()/2) ? 0.2 :0.3;
6727 // seed->SetCurrentSigmaY2(sd2+padlength*padlength*angulary/12.+sresy*sresy);
6729 // Float_t sresz = fkParam->GetZSigma();
6730 // seed->SetCurrentSigmaZ2(sd2+padlength*padlength*angularz*angularz*(1+angulary)/12.+sresz*sresz);
6732 Float_t wy = GetSigmaY(seed);
6733 Float_t wz = GetSigmaZ(seed);
6736 if (TMath::Abs(wy/seed->fCurrentSigmaY2-1)>0.0001 || TMath::Abs(wz/seed->fCurrentSigmaZ2-1)>0.0001 ){
6737 printf("problem\n");
6744 //__________________________________________________________________________
6745 void AliTPCtrackerMI::CookLabel(AliKalmanTrack *tk, Float_t wrong) const {
6746 //--------------------------------------------------------------------
6747 //This function "cooks" a track label. If label<0, this track is fake.
6748 //--------------------------------------------------------------------
6749 AliTPCseed * t = dynamic_cast<AliTPCseed*>(tk);
6751 printf("%s:%d wrong type \n",(char*)__FILE__,__LINE__);
6755 Int_t noc=t->GetNumberOfClusters();
6757 //printf("\nnot founded prolongation\n\n\n");
6763 AliTPCclusterMI *clusters[160];
6765 for (Int_t i=0;i<160;i++) {
6772 for (i=0; i<160 && current<noc; i++) {
6774 Int_t index=t->GetClusterIndex2(i);
6775 if (index<=0) continue;
6776 if (index&0x8000) continue;
6778 //clusters[current]=GetClusterMI(index);
6779 if (t->GetClusterPointer(i)){
6780 clusters[current]=t->GetClusterPointer(i);
6786 Int_t lab=123456789;
6787 for (i=0; i<noc; i++) {
6788 AliTPCclusterMI *c=clusters[i];
6790 lab=TMath::Abs(c->GetLabel(0));
6792 for (j=0; j<noc; j++) if (lb[j]==lab || mx[j]==0) break;
6798 for (i=0; i<noc; i++) if (mx[i]>max) {max=mx[i]; lab=lb[i];}
6800 for (i=0; i<noc; i++) {
6801 AliTPCclusterMI *c=clusters[i];
6803 if (TMath::Abs(c->GetLabel(1)) == lab ||
6804 TMath::Abs(c->GetLabel(2)) == lab ) max++;
6807 if ((1.- Float_t(max)/noc) > wrong) lab=-lab;
6810 Int_t tail=Int_t(0.10*noc);
6813 for (i=1; i<=160&&ind<tail; i++) {
6814 // AliTPCclusterMI *c=clusters[noc-i];
6815 AliTPCclusterMI *c=clusters[i];
6817 if (lab == TMath::Abs(c->GetLabel(0)) ||
6818 lab == TMath::Abs(c->GetLabel(1)) ||
6819 lab == TMath::Abs(c->GetLabel(2))) max++;
6822 if (max < Int_t(0.5*tail)) lab=-lab;
6829 //delete[] clusters;
6833 //__________________________________________________________________________
6834 Int_t AliTPCtrackerMI::CookLabel(AliTPCseed *t, Float_t wrong,Int_t first, Int_t last) const {
6835 //--------------------------------------------------------------------
6836 //This function "cooks" a track label. If label<0, this track is fake.
6837 //--------------------------------------------------------------------
6838 Int_t noc=t->GetNumberOfClusters();
6840 //printf("\nnot founded prolongation\n\n\n");
6846 AliTPCclusterMI *clusters[160];
6848 for (Int_t i=0;i<160;i++) {
6855 for (i=0; i<160 && current<noc; i++) {
6856 if (i<first) continue;
6857 if (i>last) continue;
6858 Int_t index=t->GetClusterIndex2(i);
6859 if (index<=0) continue;
6860 if (index&0x8000) continue;
6862 //clusters[current]=GetClusterMI(index);
6863 if (t->GetClusterPointer(i)){
6864 clusters[current]=t->GetClusterPointer(i);
6869 if (noc<5) return -1;
6870 Int_t lab=123456789;
6871 for (i=0; i<noc; i++) {
6872 AliTPCclusterMI *c=clusters[i];
6874 lab=TMath::Abs(c->GetLabel(0));
6876 for (j=0; j<noc; j++) if (lb[j]==lab || mx[j]==0) break;
6882 for (i=0; i<noc; i++) if (mx[i]>max) {max=mx[i]; lab=lb[i];}
6884 for (i=0; i<noc; i++) {
6885 AliTPCclusterMI *c=clusters[i];
6887 if (TMath::Abs(c->GetLabel(1)) == lab ||
6888 TMath::Abs(c->GetLabel(2)) == lab ) max++;
6891 if ((1.- Float_t(max)/noc) > wrong) lab=-lab;
6894 Int_t tail=Int_t(0.10*noc);
6897 for (i=1; i<=160&&ind<tail; i++) {
6898 // AliTPCclusterMI *c=clusters[noc-i];
6899 AliTPCclusterMI *c=clusters[i];
6901 if (lab == TMath::Abs(c->GetLabel(0)) ||
6902 lab == TMath::Abs(c->GetLabel(1)) ||
6903 lab == TMath::Abs(c->GetLabel(2))) max++;
6906 if (max < Int_t(0.5*tail)) lab=-lab;
6909 // t->SetLabel(lab);
6913 //delete[] clusters;
6917 Int_t AliTPCtrackerMI::GetRowNumber(Double_t x[3]) const
6919 //return pad row number for given x vector
6920 Float_t phi = TMath::ATan2(x[1],x[0]);
6921 if(phi<0) phi=2.*TMath::Pi()+phi;
6922 // Get the local angle in the sector philoc
6923 const Float_t kRaddeg = 180/3.14159265358979312;
6924 Float_t phiangle = (Int_t (phi*kRaddeg/20.) + 0.5)*20./kRaddeg;
6925 Double_t localx = x[0]*TMath::Cos(phiangle)-x[1]*TMath::Sin(phiangle);
6926 return GetRowNumber(localx);
6931 void AliTPCtrackerMI::MakeBitmaps(AliTPCseed *t)
6933 //-----------------------------------------------------------------------
6934 // Fill the cluster and sharing bitmaps of the track
6935 //-----------------------------------------------------------------------
6937 Int_t firstpoint = 0;
6938 Int_t lastpoint = 159;
6939 AliTPCTrackerPoint *point;
6940 AliTPCclusterMI *cluster;
6942 for (int iter=firstpoint; iter<lastpoint; iter++) {
6943 // Change to cluster pointers to see if we have a cluster at given padrow
6944 cluster = t->GetClusterPointer(iter);
6946 t->SetClusterMapBit(iter, kTRUE);
6947 point = t->GetTrackPoint(iter);
6948 if (point->IsShared())
6949 t->SetSharedMapBit(iter,kTRUE);
6951 t->SetSharedMapBit(iter, kFALSE);
6954 t->SetClusterMapBit(iter, kFALSE);
6955 t->SetSharedMapBit(iter, kFALSE);
6960 Bool_t AliTPCtrackerMI::IsFindable(AliTPCseed & track){
6962 // return flag if there is findable cluster at given position
6965 Float_t z = track.GetZ();
6967 if (TMath::Abs(z)<(AliTPCReconstructor::GetCtgRange()*track.GetX()+kDeltaZ) &&
6968 TMath::Abs(z)<fkParam->GetZLength(0) &&
6969 (TMath::Abs(track.GetSnp())<AliTPCReconstructor::GetMaxSnpTracker()))
6975 void AliTPCtrackerMI::AddCovariance(AliTPCseed * seed){
6977 // Adding systematic error
6978 // !!!! the systematic error for element 4 is in 1/cm not in pt
6980 const Double_t *param = AliTPCReconstructor::GetRecoParam()->GetSystematicError();
6982 for (Int_t i=0;i<15;i++) covar[i]=0;
6988 covar[0] = param[0]*param[0];
6989 covar[2] = param[1]*param[1];
6990 covar[5] = param[2]*param[2];
6991 covar[9] = param[3]*param[3];
6992 Double_t facC = AliTracker::GetBz()*kB2C;
6993 covar[14]= param[4]*param[4]*facC*facC;
6994 seed->AddCovariance(covar);