1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
17 //-------------------------------------------------------
18 // Implementation of the TPC tracker
20 // Origin: Marian Ivanov Marian.Ivanov@cern.ch
22 // AliTPC parallel tracker
24 // The track fitting is based on Kalaman filtering approach
26 // The track finding steps:
27 // 1. Seeding - with and without vertex constraint
28 // - seeding with vertex constain done at first n^2 proble
29 // - seeding without vertex constraint n^3 problem
30 // 2. Tracking - follow prolongation road - find cluster - update kalman track
32 // The seeding and tracking is repeated several times, in different seeding region.
33 // This approach enables to find the track which cannot be seeded in some region of TPC
34 // This can happen because of low momenta (track do not reach outer radius), or track is currently in the ded region between sectors, or the track is for the moment overlapped with other track (seed quality is poor) ...
36 // With this approach we reach almost 100 % efficiency also for high occupancy events.
37 // (If the seeding efficiency in a region is about 90 % than with logical or of several
38 // regions we will reach 100% (in theory - supposing independence)
40 // Repeating several seeding - tracking procedures some of the tracks can be find
43 // The procedures to remove multi find tacks are impremented:
44 // RemoveUsed2 - fast procedure n problem -
45 // Algorithm - Sorting tracks according quality
46 // remove tracks with some shared fraction
47 // Sharing in respect to all tacks
48 // Signing clusters in gold region
49 // FindSplitted - slower algorithm n^2
50 // Sort the tracks according quality
51 // Loop over pair of tracks
52 // If overlap with other track bigger than threshold - remove track
54 // FindCurling - Finds the pair of tracks which are curling
55 // - About 10% of tracks can be find with this procedure
56 // The combinatorial background is too big to be used in High
57 // multiplicity environment
58 // - n^2 problem - Slow procedure - currently it is disabled because of
61 // The number of splitted tracks can be reduced disabling the sharing of the cluster.
62 // tpcRecoParam-> SetClusterSharing(kFALSE);
63 // IT IS HIGHLY non recomended to use it in high flux enviroonment
64 // Even using this switch some tracks can be found more than once
65 // (because of multiple seeding and low quality tracks which will not cross full chamber)
68 // The tracker itself can be debugged - the information about tracks can be stored in several // phases of the reconstruction
69 // To enable storage of the TPC tracks in the ESD friend track
70 // use AliTPCReconstructor::SetStreamLevel(n);
72 // The debug level - different procedure produce tree for numerical debugging
73 // To enable them set AliTPCReconstructor::SetStreamLevel(n); where nis bigger 1
77 // Adding systematic errors to the covariance:
79 // The systematic errors due to the misalignment and miscalibration are added to the covariance matrix
80 // of the tracks (not to the clusters as they are dependent):
81 // The parameters form AliTPCRecoParam are used AliTPCRecoParam::GetSystematicError
82 // The systematic errors are expressed there in RMS - position (cm), angle (rad), curvature (1/cm)
83 // The default values are 0.
85 // The sytematic errors are added to the covariance matrix in following places:
87 // 1. During fisrt itteration - AliTPCtrackerMI::FillESD
88 // 2. Second iteration -
89 // 2.a ITS->TPC - AliTPCtrackerMI::ReadSeeds
90 // 2.b TPC->TRD - AliTPCtrackerMI::PropagateBack
91 // 3. Third iteration -
92 // 3.a TRD->TPC - AliTPCtrackerMI::ReadSeeds
93 // 3.b TPC->ITS - AliTPCtrackerMI::RefitInward
95 // There are several places in the code which can be numerically debuged
96 // This code is keeped in order to enable code development and to check the calibration implementtion
98 // 1. ErrParam stream - dump information about
100 // 2.a) cluster error estimate
101 // 3.a) cluster shape estimate
104 // Debug streamer levels:
106 //-------------------------------------------------------
111 #include "Riostream.h"
112 #include <TClonesArray.h>
114 #include <TObjArray.h>
116 #include <TGraphErrors.h>
118 #include "AliComplexCluster.h"
119 #include "AliESDEvent.h"
120 #include "AliESDtrack.h"
121 #include "AliESDVertex.h"
124 #include "AliHelix.h"
125 #include "AliRunLoader.h"
126 #include "AliTPCClustersRow.h"
127 #include "AliTPCParam.h"
128 #include "AliTPCReconstructor.h"
129 #include "AliTPCpolyTrack.h"
130 #include "AliTPCreco.h"
131 #include "AliTPCseed.h"
133 #include "AliTPCtrackerSector.h"
134 #include "AliTPCtrackerMI.h"
135 #include "TStopwatch.h"
136 #include "AliTPCReconstructor.h"
137 #include "AliAlignObj.h"
138 #include "AliTrackPointArray.h"
140 #include "AliTPCcalibDB.h"
141 #include "AliTPCcalibDButil.h"
142 #include "AliTPCTransform.h"
143 #include "AliTPCClusterParam.h"
147 ClassImp(AliTPCtrackerMI)
151 class AliTPCFastMath {
154 static Double_t FastAsin(Double_t x);
156 static Double_t fgFastAsin[20000]; //lookup table for fast asin computation
159 Double_t AliTPCFastMath::fgFastAsin[20000];
160 AliTPCFastMath gAliTPCFastMath; // needed to fill the LUT
162 AliTPCFastMath::AliTPCFastMath(){
164 // initialized lookup table;
165 for (Int_t i=0;i<10000;i++){
166 fgFastAsin[2*i] = TMath::ASin(i/10000.);
167 fgFastAsin[2*i+1] = (TMath::ASin((i+1)/10000.)-fgFastAsin[2*i]);
171 Double_t AliTPCFastMath::FastAsin(Double_t x){
173 // return asin using lookup table
175 Int_t index = int(x*10000);
176 return fgFastAsin[2*index]+(x*10000.-index)*fgFastAsin[2*index+1];
179 Int_t index = int(x*10000);
180 return -(fgFastAsin[2*index]+(x*10000.-index)*fgFastAsin[2*index+1]);
182 //__________________________________________________________________
183 AliTPCtrackerMI::AliTPCtrackerMI()
205 // default constructor
207 for (Int_t irow=0; irow<200; irow++){
214 //_____________________________________________________________________
218 Int_t AliTPCtrackerMI::UpdateTrack(AliTPCseed * track, Int_t accept){
220 //update track information using current cluster - track->fCurrentCluster
223 AliTPCclusterMI* c =track->GetCurrentCluster();
224 if (accept > 0) //sign not accepted clusters
225 track->SetCurrentClusterIndex1(track->GetCurrentClusterIndex1() | 0x8000);
226 else // unsign accpeted clusters
227 track->SetCurrentClusterIndex1(track->GetCurrentClusterIndex1() & 0xffff7fff);
228 UInt_t i = track->GetCurrentClusterIndex1();
230 Int_t sec=(i&0xff000000)>>24;
231 //Int_t row = (i&0x00ff0000)>>16;
232 track->SetRow((i&0x00ff0000)>>16);
233 track->SetSector(sec);
234 // Int_t index = i&0xFFFF;
235 if (sec>=fkParam->GetNInnerSector()) track->SetRow(track->GetRow()+fkParam->GetNRowLow());
236 track->SetClusterIndex2(track->GetRow(), i);
237 //track->fFirstPoint = row;
238 //if ( track->fLastPoint<row) track->fLastPoint =row;
239 // if (track->fRow<0 || track->fRow>160) {
240 // printf("problem\n");
242 if (track->GetFirstPoint()>track->GetRow())
243 track->SetFirstPoint(track->GetRow());
244 if (track->GetLastPoint()<track->GetRow())
245 track->SetLastPoint(track->GetRow());
248 track->SetClusterPointer(track->GetRow(),c);
251 Double_t angle2 = track->GetSnp()*track->GetSnp();
253 //SET NEW Track Point
255 if (angle2<1) //PH sometimes angle2 is very big. To be investigated...
257 angle2 = TMath::Sqrt(angle2/(1-angle2));
258 AliTPCTrackerPoint &point =*(track->GetTrackPoint(track->GetRow()));
260 point.SetSigmaY(c->GetSigmaY2()/track->GetCurrentSigmaY2());
261 point.SetSigmaZ(c->GetSigmaZ2()/track->GetCurrentSigmaZ2());
262 point.SetErrY(sqrt(track->GetErrorY2()));
263 point.SetErrZ(sqrt(track->GetErrorZ2()));
265 point.SetX(track->GetX());
266 point.SetY(track->GetY());
267 point.SetZ(track->GetZ());
268 point.SetAngleY(angle2);
269 point.SetAngleZ(track->GetTgl());
270 if (point.IsShared()){
271 track->SetErrorY2(track->GetErrorY2()*4);
272 track->SetErrorZ2(track->GetErrorZ2()*4);
276 Double_t chi2 = track->GetPredictedChi2(track->GetCurrentCluster());
278 // track->SetErrorY2(track->GetErrorY2()*1.3);
279 // track->SetErrorY2(track->GetErrorY2()+0.01);
280 // track->SetErrorZ2(track->GetErrorZ2()*1.3);
281 // track->SetErrorZ2(track->GetErrorZ2()+0.005);
283 if (accept>0) return 0;
284 if (track->GetNumberOfClusters()%20==0){
285 // if (track->fHelixIn){
286 // TClonesArray & larr = *(track->fHelixIn);
287 // Int_t ihelix = larr.GetEntriesFast();
288 // new(larr[ihelix]) AliHelix(*track) ;
291 track->SetNoCluster(0);
292 return track->Update(c,chi2,i);
297 Int_t AliTPCtrackerMI::AcceptCluster(AliTPCseed * seed, AliTPCclusterMI * cluster)
300 // decide according desired precision to accept given
301 // cluster for tracking
303 seed->GetProlongation(cluster->GetX(),yt,zt);
304 Double_t sy2=ErrY2(seed,cluster);
305 Double_t sz2=ErrZ2(seed,cluster);
307 Double_t sdistancey2 = sy2+seed->GetSigmaY2();
308 Double_t sdistancez2 = sz2+seed->GetSigmaZ2();
309 Double_t dy=seed->GetCurrentCluster()->GetY()-yt;
310 Double_t dz=seed->GetCurrentCluster()->GetZ()-zt;
311 Double_t rdistancey2 = (seed->GetCurrentCluster()->GetY()-yt)*
312 (seed->GetCurrentCluster()->GetY()-yt)/sdistancey2;
313 Double_t rdistancez2 = (seed->GetCurrentCluster()->GetZ()-zt)*
314 (seed->GetCurrentCluster()->GetZ()-zt)/sdistancez2;
316 Double_t rdistance2 = rdistancey2+rdistancez2;
319 if (AliTPCReconstructor::StreamLevel()>2 && seed->GetNumberOfClusters()>20) {
320 Float_t rmsy2 = seed->GetCurrentSigmaY2();
321 Float_t rmsz2 = seed->GetCurrentSigmaZ2();
322 Float_t rmsy2p30 = seed->GetCMeanSigmaY2p30();
323 Float_t rmsz2p30 = seed->GetCMeanSigmaZ2p30();
324 Float_t rmsy2p30R = seed->GetCMeanSigmaY2p30R();
325 Float_t rmsz2p30R = seed->GetCMeanSigmaZ2p30R();
326 AliExternalTrackParam param(*seed);
327 static TVectorD gcl(3),gtr(3);
329 param.GetXYZ(gcl.GetMatrixArray());
330 cluster->GetGlobalXYZ(gclf);
331 gcl[0]=gclf[0]; gcl[1]=gclf[1]; gcl[2]=gclf[2];
334 if (AliTPCReconstructor::StreamLevel()>2) {
335 (*fDebugStreamer)<<"ErrParam"<<
348 "rmsy2p30="<<rmsy2p30<<
349 "rmsz2p30="<<rmsz2p30<<
350 "rmsy2p30R="<<rmsy2p30R<<
351 "rmsz2p30R="<<rmsz2p30R<<
352 // normalize distance -
353 "rdisty="<<rdistancey2<<
354 "rdistz="<<rdistancez2<<
355 "rdist="<<rdistance2<< //
359 //return 0; // temporary
360 if (rdistance2>32) return 3;
363 if ((rdistancey2>9. || rdistancez2>9.) && cluster->GetType()==0)
364 return 2; //suspisiouce - will be changed
366 if ((rdistancey2>6.25 || rdistancez2>6.25) && cluster->GetType()>0)
367 // strict cut on overlaped cluster
368 return 2; //suspisiouce - will be changed
370 if ( (rdistancey2>1. || rdistancez2>6.25 )
371 && cluster->GetType()<0){
372 seed->SetNFoundable(seed->GetNFoundable()-1);
376 if (AliTPCReconstructor::GetRecoParam()->GetUseHLTClusters() == 3 ||
377 AliTPCReconstructor::GetRecoParam()->GetUseHLTClusters() == 4 ) {
378 if(!AliTPCReconstructor::GetRecoParam()->GetUseOnePadCluster())
379 if (TMath::Abs(cluster->GetSigmaY2()) < kAlmost0)
389 //_____________________________________________________________________________
390 AliTPCtrackerMI::AliTPCtrackerMI(const AliTPCParam *par):
392 fkNIS(par->GetNInnerSector()/2),
394 fkNOS(par->GetNOuterSector()/2),
411 //---------------------------------------------------------------------
412 // The main TPC tracker constructor
413 //---------------------------------------------------------------------
414 fInnerSec=new AliTPCtrackerSector[fkNIS];
415 fOuterSec=new AliTPCtrackerSector[fkNOS];
418 for (i=0; i<fkNIS; i++) fInnerSec[i].Setup(par,0);
419 for (i=0; i<fkNOS; i++) fOuterSec[i].Setup(par,1);
422 Int_t nrowlow = par->GetNRowLow();
423 Int_t nrowup = par->GetNRowUp();
426 for (i=0;i<nrowlow;i++){
427 fXRow[i] = par->GetPadRowRadiiLow(i);
428 fPadLength[i]= par->GetPadPitchLength(0,i);
429 fYMax[i] = fXRow[i]*TMath::Tan(0.5*par->GetInnerAngle());
433 for (i=0;i<nrowup;i++){
434 fXRow[i+nrowlow] = par->GetPadRowRadiiUp(i);
435 fPadLength[i+nrowlow] = par->GetPadPitchLength(60,i);
436 fYMax[i+nrowlow] = fXRow[i+nrowlow]*TMath::Tan(0.5*par->GetOuterAngle());
439 if (AliTPCReconstructor::StreamLevel()>0) {
440 fDebugStreamer = new TTreeSRedirector("TPCdebug.root");
443 //________________________________________________________________________
444 AliTPCtrackerMI::AliTPCtrackerMI(const AliTPCtrackerMI &t):
465 //------------------------------------
466 // dummy copy constructor
467 //------------------------------------------------------------------
469 for (Int_t irow=0; irow<200; irow++){
476 AliTPCtrackerMI & AliTPCtrackerMI::operator=(const AliTPCtrackerMI& /*r*/)
478 //------------------------------
480 //--------------------------------------------------------------
483 //_____________________________________________________________________________
484 AliTPCtrackerMI::~AliTPCtrackerMI() {
485 //------------------------------------------------------------------
486 // TPC tracker destructor
487 //------------------------------------------------------------------
494 if (fDebugStreamer) delete fDebugStreamer;
498 void AliTPCtrackerMI::FillESD(const TObjArray* arr)
502 //fill esds using updated tracks
505 // write tracks to the event
506 // store index of the track
507 Int_t nseed=arr->GetEntriesFast();
508 //FindKinks(arr,fEvent);
509 for (Int_t i=0; i<nseed; i++) {
510 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
514 if (AliTPCReconstructor::StreamLevel()>1) {
515 (*fDebugStreamer)<<"Track0"<<
519 // pt->PropagateTo(fkParam->GetInnerRadiusLow());
520 if (pt->GetKinkIndex(0)<=0){ //don't propagate daughter tracks
521 pt->PropagateTo(fkParam->GetInnerRadiusLow());
524 if (( pt->GetPoints()[2]- pt->GetPoints()[0])>5 && pt->GetPoints()[3]>0.8){
526 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
527 iotrack.SetTPCPoints(pt->GetPoints());
528 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
529 iotrack.SetV0Indexes(pt->GetV0Indexes());
530 // iotrack.SetTPCpid(pt->fTPCr);
531 //iotrack.SetTPCindex(i);
532 fEvent->AddTrack(&iotrack);
536 if ( (pt->GetNumberOfClusters()>70)&& (Float_t(pt->GetNumberOfClusters())/Float_t(pt->GetNFoundable()))>0.55) {
538 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
539 iotrack.SetTPCPoints(pt->GetPoints());
540 //iotrack.SetTPCindex(i);
541 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
542 iotrack.SetV0Indexes(pt->GetV0Indexes());
543 // iotrack.SetTPCpid(pt->fTPCr);
544 fEvent->AddTrack(&iotrack);
548 // short tracks - maybe decays
550 if ( (pt->GetNumberOfClusters()>30) && (Float_t(pt->GetNumberOfClusters())/Float_t(pt->GetNFoundable()))>0.70) {
551 Int_t found,foundable,shared;
552 pt->GetClusterStatistic(0,60,found, foundable,shared,kFALSE);
553 if ( (found>20) && (pt->GetNShared()/float(pt->GetNumberOfClusters())<0.2)){
555 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
556 //iotrack.SetTPCindex(i);
557 iotrack.SetTPCPoints(pt->GetPoints());
558 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
559 iotrack.SetV0Indexes(pt->GetV0Indexes());
560 //iotrack.SetTPCpid(pt->fTPCr);
561 fEvent->AddTrack(&iotrack);
566 if ( (pt->GetNumberOfClusters()>20) && (Float_t(pt->GetNumberOfClusters())/Float_t(pt->GetNFoundable()))>0.8) {
567 Int_t found,foundable,shared;
568 pt->GetClusterStatistic(0,60,found, foundable,shared,kFALSE);
569 if (found<20) continue;
570 if (pt->GetNShared()/float(pt->GetNumberOfClusters())>0.2) continue;
573 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
574 iotrack.SetTPCPoints(pt->GetPoints());
575 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
576 iotrack.SetV0Indexes(pt->GetV0Indexes());
577 //iotrack.SetTPCpid(pt->fTPCr);
578 //iotrack.SetTPCindex(i);
579 fEvent->AddTrack(&iotrack);
582 // short tracks - secondaties
584 if ( (pt->GetNumberOfClusters()>30) ) {
585 Int_t found,foundable,shared;
586 pt->GetClusterStatistic(128,158,found, foundable,shared,kFALSE);
587 if ( (found>20) && (pt->GetNShared()/float(pt->GetNumberOfClusters())<0.2) &&float(found)/float(foundable)>0.8){
589 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
590 iotrack.SetTPCPoints(pt->GetPoints());
591 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
592 iotrack.SetV0Indexes(pt->GetV0Indexes());
593 //iotrack.SetTPCpid(pt->fTPCr);
594 //iotrack.SetTPCindex(i);
595 fEvent->AddTrack(&iotrack);
600 if ( (pt->GetNumberOfClusters()>15)) {
601 Int_t found,foundable,shared;
602 pt->GetClusterStatistic(138,158,found, foundable,shared,kFALSE);
603 if (found<15) continue;
604 if (foundable<=0) continue;
605 if (pt->GetNShared()/float(pt->GetNumberOfClusters())>0.2) continue;
606 if (float(found)/float(foundable)<0.8) continue;
609 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
610 iotrack.SetTPCPoints(pt->GetPoints());
611 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
612 iotrack.SetV0Indexes(pt->GetV0Indexes());
613 // iotrack.SetTPCpid(pt->fTPCr);
614 //iotrack.SetTPCindex(i);
615 fEvent->AddTrack(&iotrack);
619 // >> account for suppressed tracks in the kink indices (RS)
620 int nESDtracks = fEvent->GetNumberOfTracks();
621 for (int it=nESDtracks;it--;) {
622 AliESDtrack* esdTr = fEvent->GetTrack(it);
623 if (!esdTr || !esdTr->GetKinkIndex(0)) continue;
624 for (int ik=0;ik<3;ik++) {
626 if (!(knkId=esdTr->GetKinkIndex(ik))) break; // no more kinks for this track
627 AliESDkink* kink = fEvent->GetKink(TMath::Abs(knkId)-1);
629 AliError(Form("ESDTrack%d refers to non-existing kink %d",it,TMath::Abs(knkId)-1));
632 kink->SetIndex(it, knkId<0 ? 0:1); // update track index of the kink: mother at 0, daughter at 1
635 // << account for suppressed tracks in the kink indices (RS)
636 AliInfo(Form("Number of filled ESDs-\t%d\n",fEvent->GetNumberOfTracks()));
644 Double_t AliTPCtrackerMI::ErrY2(AliTPCseed* seed, const AliTPCclusterMI * cl){
647 // Use calibrated cluster error from OCDB
649 AliTPCClusterParam * clparam = AliTPCcalibDB::Instance()->GetClusterParam();
651 Float_t z = TMath::Abs(fkParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
652 Int_t ctype = cl->GetType();
653 Int_t type = (cl->GetRow()<63) ? 0: (cl->GetRow()>126) ? 1:2;
654 Double_t angle = seed->GetSnp()*seed->GetSnp();
655 angle = TMath::Sqrt(TMath::Abs(angle/(1.-angle)));
656 Double_t erry2 = clparam->GetError0Par(0,type, z,angle);
658 erry2+=0.5; // edge cluster
661 seed->SetErrorY2(erry2);
665 //calculate look-up table at the beginning
666 // static Bool_t ginit = kFALSE;
667 // static Float_t gnoise1,gnoise2,gnoise3;
668 // static Float_t ggg1[10000];
669 // static Float_t ggg2[10000];
670 // static Float_t ggg3[10000];
671 // static Float_t glandau1[10000];
672 // static Float_t glandau2[10000];
673 // static Float_t glandau3[10000];
675 // static Float_t gcor01[500];
676 // static Float_t gcor02[500];
677 // static Float_t gcorp[500];
681 // if (ginit==kFALSE){
682 // for (Int_t i=1;i<500;i++){
683 // Float_t rsigma = float(i)/100.;
684 // gcor02[i] = TMath::Max(0.78 +TMath::Exp(7.4*(rsigma-1.2)),0.6);
685 // gcor01[i] = TMath::Max(0.72 +TMath::Exp(3.36*(rsigma-1.2)),0.6);
686 // gcorp[i] = TMath::Max(TMath::Power((rsigma+0.5),1.5),1.2);
690 // for (Int_t i=3;i<10000;i++){
694 // Float_t amp = float(i);
695 // Float_t padlength =0.75;
696 // gnoise1 = 0.0004/padlength;
697 // Float_t nel = 0.268*amp;
698 // Float_t nprim = 0.155*amp;
699 // ggg1[i] = fkParam->GetDiffT()*fkParam->GetDiffT()*(2+0.001*nel/(padlength*padlength))/nel;
700 // glandau1[i] = (2.+0.12*nprim)*0.5* (2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
701 // if (glandau1[i]>1) glandau1[i]=1;
702 // glandau1[i]*=padlength*padlength/12.;
706 // gnoise2 = 0.0004/padlength;
708 // nprim = 0.133*amp;
709 // ggg2[i] = fkParam->GetDiffT()*fkParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
710 // glandau2[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
711 // if (glandau2[i]>1) glandau2[i]=1;
712 // glandau2[i]*=padlength*padlength/12.;
717 // gnoise3 = 0.0004/padlength;
719 // nprim = 0.133*amp;
720 // ggg3[i] = fkParam->GetDiffT()*fkParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
721 // glandau3[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
722 // if (glandau3[i]>1) glandau3[i]=1;
723 // glandau3[i]*=padlength*padlength/12.;
731 // Int_t amp = int(TMath::Abs(cl->GetQ()));
733 // seed->SetErrorY2(1.);
737 // Float_t z = TMath::Abs(fkParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
738 // Int_t ctype = cl->GetType();
739 // Float_t padlength= GetPadPitchLength(seed->GetRow());
740 // Double_t angle2 = seed->GetSnp()*seed->GetSnp();
741 // angle2 = angle2/(1-angle2);
743 // //cluster "quality"
744 // Int_t rsigmay = int(100.*cl->GetSigmaY2()/(seed->GetCurrentSigmaY2()));
747 // if (fSectors==fInnerSec){
748 // snoise2 = gnoise1;
749 // res = ggg1[amp]*z+glandau1[amp]*angle2;
750 // if (ctype==0) res *= gcor01[rsigmay];
753 // res*= gcorp[rsigmay];
757 // if (padlength<1.1){
758 // snoise2 = gnoise2;
759 // res = ggg2[amp]*z+glandau2[amp]*angle2;
760 // if (ctype==0) res *= gcor02[rsigmay];
763 // res*= gcorp[rsigmay];
767 // snoise2 = gnoise3;
768 // res = ggg3[amp]*z+glandau3[amp]*angle2;
769 // if (ctype==0) res *= gcor02[rsigmay];
772 // res*= gcorp[rsigmay];
779 // res*=2.4; // overestimate error 2 times
783 // if (res<2*snoise2)
786 // seed->SetErrorY2(res);
794 Double_t AliTPCtrackerMI::ErrZ2(AliTPCseed* seed, const AliTPCclusterMI * cl){
797 // Use calibrated cluster error from OCDB
799 AliTPCClusterParam * clparam = AliTPCcalibDB::Instance()->GetClusterParam();
801 Float_t z = TMath::Abs(fkParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
802 Int_t ctype = cl->GetType();
803 Int_t type = (cl->GetRow()<63) ? 0: (cl->GetRow()>126) ? 1:2;
805 Double_t angle2 = seed->GetSnp()*seed->GetSnp();
806 angle2 = seed->GetTgl()*seed->GetTgl()*(1+angle2/(1-angle2));
807 Double_t angle = TMath::Sqrt(TMath::Abs(angle2));
808 Double_t errz2 = clparam->GetError0Par(1,type, z,angle);
810 errz2+=0.5; // edge cluster
813 seed->SetErrorZ2(errz2);
819 // //seed->SetErrorY2(0.1);
821 // //calculate look-up table at the beginning
822 // static Bool_t ginit = kFALSE;
823 // static Float_t gnoise1,gnoise2,gnoise3;
824 // static Float_t ggg1[10000];
825 // static Float_t ggg2[10000];
826 // static Float_t ggg3[10000];
827 // static Float_t glandau1[10000];
828 // static Float_t glandau2[10000];
829 // static Float_t glandau3[10000];
831 // static Float_t gcor01[1000];
832 // static Float_t gcor02[1000];
833 // static Float_t gcorp[1000];
837 // if (ginit==kFALSE){
838 // for (Int_t i=1;i<1000;i++){
839 // Float_t rsigma = float(i)/100.;
840 // gcor02[i] = TMath::Max(0.81 +TMath::Exp(6.8*(rsigma-1.2)),0.6);
841 // gcor01[i] = TMath::Max(0.72 +TMath::Exp(2.04*(rsigma-1.2)),0.6);
842 // gcorp[i] = TMath::Max(TMath::Power((rsigma+0.5),1.5),1.2);
846 // for (Int_t i=3;i<10000;i++){
850 // Float_t amp = float(i);
851 // Float_t padlength =0.75;
852 // gnoise1 = 0.0004/padlength;
853 // Float_t nel = 0.268*amp;
854 // Float_t nprim = 0.155*amp;
855 // ggg1[i] = fkParam->GetDiffT()*fkParam->GetDiffT()*(2+0.001*nel/(padlength*padlength))/nel;
856 // glandau1[i] = (2.+0.12*nprim)*0.5* (2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
857 // if (glandau1[i]>1) glandau1[i]=1;
858 // glandau1[i]*=padlength*padlength/12.;
862 // gnoise2 = 0.0004/padlength;
864 // nprim = 0.133*amp;
865 // ggg2[i] = fkParam->GetDiffT()*fkParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
866 // glandau2[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
867 // if (glandau2[i]>1) glandau2[i]=1;
868 // glandau2[i]*=padlength*padlength/12.;
873 // gnoise3 = 0.0004/padlength;
875 // nprim = 0.133*amp;
876 // ggg3[i] = fkParam->GetDiffT()*fkParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
877 // glandau3[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
878 // if (glandau3[i]>1) glandau3[i]=1;
879 // glandau3[i]*=padlength*padlength/12.;
887 // Int_t amp = int(TMath::Abs(cl->GetQ()));
889 // seed->SetErrorY2(1.);
893 // Float_t z = TMath::Abs(fkParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
894 // Int_t ctype = cl->GetType();
895 // Float_t padlength= GetPadPitchLength(seed->GetRow());
897 // Double_t angle2 = seed->GetSnp()*seed->GetSnp();
898 // // if (angle2<0.6) angle2 = 0.6;
899 // angle2 = seed->GetTgl()*seed->GetTgl()*(1+angle2/(1-angle2));
901 // //cluster "quality"
902 // Int_t rsigmaz = int(100.*cl->GetSigmaZ2()/(seed->GetCurrentSigmaZ2()));
905 // if (fSectors==fInnerSec){
906 // snoise2 = gnoise1;
907 // res = ggg1[amp]*z+glandau1[amp]*angle2;
908 // if (ctype==0) res *= gcor01[rsigmaz];
911 // res*= gcorp[rsigmaz];
915 // if (padlength<1.1){
916 // snoise2 = gnoise2;
917 // res = ggg2[amp]*z+glandau2[amp]*angle2;
918 // if (ctype==0) res *= gcor02[rsigmaz];
921 // res*= gcorp[rsigmaz];
925 // snoise2 = gnoise3;
926 // res = ggg3[amp]*z+glandau3[amp]*angle2;
927 // if (ctype==0) res *= gcor02[rsigmaz];
930 // res*= gcorp[rsigmaz];
939 // if ((ctype<0) &&<70){
944 // if (res<2*snoise2)
946 // if (res>3) res =3;
947 // seed->SetErrorZ2(res);
955 void AliTPCtrackerMI::RotateToLocal(AliTPCseed *seed)
957 //rotate to track "local coordinata
958 Float_t x = seed->GetX();
959 Float_t y = seed->GetY();
960 Float_t ymax = x*TMath::Tan(0.5*fSectors->GetAlpha());
963 seed->SetRelativeSector((seed->GetRelativeSector()+1) % fN);
964 if (!seed->Rotate(fSectors->GetAlpha()))
966 } else if (y <-ymax) {
967 seed->SetRelativeSector((seed->GetRelativeSector()-1+fN) % fN);
968 if (!seed->Rotate(-fSectors->GetAlpha()))
976 //_____________________________________________________________________________
977 Double_t AliTPCtrackerMI::F1old(Double_t x1,Double_t y1,
978 Double_t x2,Double_t y2,
979 Double_t x3,Double_t y3) const
981 //-----------------------------------------------------------------
982 // Initial approximation of the track curvature
983 //-----------------------------------------------------------------
984 Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
985 Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
986 (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
987 Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
988 (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
990 Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
991 if ( xr*xr+yr*yr<=0.00000000000001) return 100;
992 return -xr*yr/sqrt(xr*xr+yr*yr);
997 //_____________________________________________________________________________
998 Double_t AliTPCtrackerMI::F1(Double_t x1,Double_t y1,
999 Double_t x2,Double_t y2,
1000 Double_t x3,Double_t y3) const
1002 //-----------------------------------------------------------------
1003 // Initial approximation of the track curvature
1004 //-----------------------------------------------------------------
1010 Double_t det = x3*y2-x2*y3;
1011 if (TMath::Abs(det)<1e-10){
1015 Double_t u = 0.5* (x2*(x2-x3)+y2*(y2-y3))/det;
1016 Double_t x0 = x3*0.5-y3*u;
1017 Double_t y0 = y3*0.5+x3*u;
1018 Double_t c2 = 1/TMath::Sqrt(x0*x0+y0*y0);
1024 Double_t AliTPCtrackerMI::F2(Double_t x1,Double_t y1,
1025 Double_t x2,Double_t y2,
1026 Double_t x3,Double_t y3) const
1028 //-----------------------------------------------------------------
1029 // Initial approximation of the track curvature
1030 //-----------------------------------------------------------------
1036 Double_t det = x3*y2-x2*y3;
1037 if (TMath::Abs(det)<1e-10) {
1041 Double_t u = 0.5* (x2*(x2-x3)+y2*(y2-y3))/det;
1042 Double_t x0 = x3*0.5-y3*u;
1043 Double_t y0 = y3*0.5+x3*u;
1044 Double_t c2 = 1/TMath::Sqrt(x0*x0+y0*y0);
1053 //_____________________________________________________________________________
1054 Double_t AliTPCtrackerMI::F2old(Double_t x1,Double_t y1,
1055 Double_t x2,Double_t y2,
1056 Double_t x3,Double_t y3) const
1058 //-----------------------------------------------------------------
1059 // Initial approximation of the track curvature times center of curvature
1060 //-----------------------------------------------------------------
1061 Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
1062 Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
1063 (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
1064 Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
1065 (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
1067 Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
1069 return -a/(d*y1-b)*xr/sqrt(xr*xr+yr*yr);
1072 //_____________________________________________________________________________
1073 Double_t AliTPCtrackerMI::F3(Double_t x1,Double_t y1,
1074 Double_t x2,Double_t y2,
1075 Double_t z1,Double_t z2) const
1077 //-----------------------------------------------------------------
1078 // Initial approximation of the tangent of the track dip angle
1079 //-----------------------------------------------------------------
1080 return (z1 - z2)/sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
1084 Double_t AliTPCtrackerMI::F3n(Double_t x1,Double_t y1,
1085 Double_t x2,Double_t y2,
1086 Double_t z1,Double_t z2, Double_t c) const
1088 //-----------------------------------------------------------------
1089 // Initial approximation of the tangent of the track dip angle
1090 //-----------------------------------------------------------------
1094 //angle1 = (z1-z2)*c/(TMath::ASin(c*x1-ni)-TMath::ASin(c*x2-ni));
1096 Double_t d = TMath::Sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
1097 if (TMath::Abs(d*c*0.5)>1) return 0;
1098 // Double_t angle2 = TMath::ASin(d*c*0.5);
1099 // Double_t angle2 = AliTPCFastMath::FastAsin(d*c*0.5);
1100 Double_t angle2 = (d*c*0.5>0.1)? TMath::ASin(d*c*0.5): AliTPCFastMath::FastAsin(d*c*0.5);
1102 angle2 = (z1-z2)*c/(angle2*2.);
1106 Bool_t AliTPCtrackerMI::GetProlongation(Double_t x1, Double_t x2, Double_t x[5], Double_t &y, Double_t &z) const
1107 {//-----------------------------------------------------------------
1108 // This function find proloncation of a track to a reference plane x=x2.
1109 //-----------------------------------------------------------------
1113 if (TMath::Abs(x[4]*x1 - x[2]) >= 0.999) {
1117 Double_t c1=x[4]*x1 - x[2], r1=TMath::Sqrt((1.-c1)*(1.+c1));
1118 Double_t c2=x[4]*x2 - x[2], r2=TMath::Sqrt((1.-c2)*(1.+c2));
1122 Double_t dy = dx*(c1+c2)/(r1+r2);
1125 Double_t delta = x[4]*dx*(c1+c2)/(c1*r2 + c2*r1);
1127 if (TMath::Abs(delta)>0.01){
1128 dz = x[3]*TMath::ASin(delta)/x[4];
1130 dz = x[3]*AliTPCFastMath::FastAsin(delta)/x[4];
1133 //dz = x[3]*AliTPCFastMath::FastAsin(delta)/x[4];
1141 Int_t AliTPCtrackerMI::LoadClusters (TTree *const tree)
1146 return LoadClusters();
1150 Int_t AliTPCtrackerMI::LoadClusters(const TObjArray *arr)
1153 // load clusters to the memory
1154 AliTPCClustersRow *clrow = new AliTPCClustersRow("AliTPCclusterMI");
1155 Int_t lower = arr->LowerBound();
1156 Int_t entries = arr->GetEntriesFast();
1158 for (Int_t i=lower; i<entries; i++) {
1159 clrow = (AliTPCClustersRow*) arr->At(i);
1160 if(!clrow) continue;
1161 if(!clrow->GetArray()) continue;
1165 fkParam->AdjustSectorRow(clrow->GetID(),sec,row);
1167 for (Int_t icl=0; icl<clrow->GetArray()->GetEntriesFast(); icl++){
1168 Transform((AliTPCclusterMI*)(clrow->GetArray()->At(icl)));
1171 if (clrow->GetArray()->GetEntriesFast()<=0) continue;
1172 AliTPCtrackerRow * tpcrow=0;
1175 tpcrow = &(fInnerSec[sec%fkNIS][row]);
1179 tpcrow = &(fOuterSec[(sec-fkNIS*2)%fkNOS][row]);
1180 left = (sec-fkNIS*2)/fkNOS;
1183 tpcrow->SetN1(clrow->GetArray()->GetEntriesFast());
1184 tpcrow->SetClusters1(new AliTPCclusterMI[tpcrow->GetN1()]);
1185 for (Int_t j=0;j<tpcrow->GetN1();++j)
1186 tpcrow->SetCluster1(j, *(AliTPCclusterMI*)(clrow->GetArray()->At(j)));
1189 tpcrow->SetN2(clrow->GetArray()->GetEntriesFast());
1190 tpcrow->SetClusters2(new AliTPCclusterMI[tpcrow->GetN2()]);
1191 for (Int_t j=0;j<tpcrow->GetN2();++j)
1192 tpcrow->SetCluster2(j,*(AliTPCclusterMI*)(clrow->GetArray()->At(j)));
1194 clrow->GetArray()->Clear("C");
1203 Int_t AliTPCtrackerMI::LoadClusters(const TClonesArray *arr)
1206 // load clusters to the memory from one
1209 AliTPCclusterMI *clust=0;
1210 Int_t count[72][96] = { {0} , {0} };
1212 // loop over clusters
1213 for (Int_t icl=0; icl<arr->GetEntriesFast(); icl++) {
1214 clust = (AliTPCclusterMI*)arr->At(icl);
1215 if(!clust) continue;
1216 //printf("cluster: det %d, row %d \n", clust->GetDetector(),clust->GetRow());
1218 // transform clusters
1221 // count clusters per pad row
1222 count[clust->GetDetector()][clust->GetRow()]++;
1225 // insert clusters to sectors
1226 for (Int_t icl=0; icl<arr->GetEntriesFast(); icl++) {
1227 clust = (AliTPCclusterMI*)arr->At(icl);
1228 if(!clust) continue;
1230 Int_t sec = clust->GetDetector();
1231 Int_t row = clust->GetRow();
1233 // filter overlapping pad rows needed by HLT
1234 if(sec<fkNIS*2) { //IROCs
1235 if(row == 30) continue;
1238 if(row == 27 || row == 76) continue;
1244 fInnerSec[sec%fkNIS].InsertCluster(clust, count[sec][row], fkParam);
1247 left = (sec-fkNIS*2)/fkNOS;
1248 fOuterSec[(sec-fkNIS*2)%fkNOS].InsertCluster(clust, count[sec][row], fkParam);
1252 // Load functions must be called behind LoadCluster(TClonesArray*)
1254 //LoadOuterSectors();
1255 //LoadInnerSectors();
1261 Int_t AliTPCtrackerMI::LoadClusters()
1264 // load clusters to the memory
1265 AliTPCClustersRow *clrow= new AliTPCClustersRow("AliTPCclusterMI");
1267 // TTree * tree = fClustersArray.GetTree();
1269 TTree * tree = fInput;
1270 TBranch * br = tree->GetBranch("Segment");
1271 br->SetAddress(&clrow);
1273 Int_t j=Int_t(tree->GetEntries());
1274 for (Int_t i=0; i<j; i++) {
1278 fkParam->AdjustSectorRow(clrow->GetID(),sec,row);
1279 for (Int_t icl=0; icl<clrow->GetArray()->GetEntriesFast(); icl++){
1280 Transform((AliTPCclusterMI*)(clrow->GetArray()->At(icl)));
1283 AliTPCtrackerRow * tpcrow=0;
1286 tpcrow = &(fInnerSec[sec%fkNIS][row]);
1290 tpcrow = &(fOuterSec[(sec-fkNIS*2)%fkNOS][row]);
1291 left = (sec-fkNIS*2)/fkNOS;
1294 tpcrow->SetN1(clrow->GetArray()->GetEntriesFast());
1295 tpcrow->SetClusters1(new AliTPCclusterMI[tpcrow->GetN1()]);
1296 for (Int_t k=0;k<tpcrow->GetN1();++k)
1297 tpcrow->SetCluster1(k, *(AliTPCclusterMI*)(clrow->GetArray()->At(k)));
1300 tpcrow->SetN2(clrow->GetArray()->GetEntriesFast());
1301 tpcrow->SetClusters2(new AliTPCclusterMI[tpcrow->GetN2()]);
1302 for (Int_t k=0;k<tpcrow->GetN2();++k)
1303 tpcrow->SetCluster2(k,*(AliTPCclusterMI*)(clrow->GetArray()->At(k)));
1314 void AliTPCtrackerMI::UnloadClusters()
1317 // unload clusters from the memory
1319 Int_t nrows = fOuterSec->GetNRows();
1320 for (Int_t sec = 0;sec<fkNOS;sec++)
1321 for (Int_t row = 0;row<nrows;row++){
1322 AliTPCtrackerRow* tpcrow = &(fOuterSec[sec%fkNOS][row]);
1324 // if (tpcrow->fClusters1) delete []tpcrow->fClusters1;
1325 // if (tpcrow->fClusters2) delete []tpcrow->fClusters2;
1327 tpcrow->ResetClusters();
1330 nrows = fInnerSec->GetNRows();
1331 for (Int_t sec = 0;sec<fkNIS;sec++)
1332 for (Int_t row = 0;row<nrows;row++){
1333 AliTPCtrackerRow* tpcrow = &(fInnerSec[sec%fkNIS][row]);
1335 // if (tpcrow->fClusters1) delete []tpcrow->fClusters1;
1336 //if (tpcrow->fClusters2) delete []tpcrow->fClusters2;
1338 tpcrow->ResetClusters();
1344 void AliTPCtrackerMI::FillClusterArray(TObjArray* array) const{
1346 // Filling cluster to the array - For visualization purposes
1349 nrows = fOuterSec->GetNRows();
1350 for (Int_t sec = 0;sec<fkNOS;sec++)
1351 for (Int_t row = 0;row<nrows;row++){
1352 AliTPCtrackerRow* tpcrow = &(fOuterSec[sec%fkNOS][row]);
1353 if (!tpcrow) continue;
1354 for (Int_t icl = 0;icl<tpcrow->GetN();icl++){
1355 array->AddLast((TObject*)((*tpcrow)[icl]));
1358 nrows = fInnerSec->GetNRows();
1359 for (Int_t sec = 0;sec<fkNIS;sec++)
1360 for (Int_t row = 0;row<nrows;row++){
1361 AliTPCtrackerRow* tpcrow = &(fInnerSec[sec%fkNIS][row]);
1362 if (!tpcrow) continue;
1363 for (Int_t icl = 0;icl<tpcrow->GetN();icl++){
1364 array->AddLast((TObject*)(*tpcrow)[icl]);
1370 void AliTPCtrackerMI::Transform(AliTPCclusterMI * cluster){
1374 AliTPCcalibDB * calibDB = AliTPCcalibDB::Instance();
1375 AliTPCTransform *transform = calibDB->GetTransform() ;
1377 AliFatal("Tranformations not in calibDB");
1380 transform->SetCurrentRecoParam((AliTPCRecoParam*)AliTPCReconstructor::GetRecoParam());
1381 Double_t x[3]={cluster->GetRow(),cluster->GetPad(),cluster->GetTimeBin()};
1382 Int_t i[1]={cluster->GetDetector()};
1383 transform->Transform(x,i,0,1);
1384 // if (cluster->GetDetector()%36>17){
1389 // in debug mode check the transformation
1391 if (AliTPCReconstructor::StreamLevel()>2) {
1393 cluster->GetGlobalXYZ(gx);
1394 Int_t event = (fEvent==NULL)? 0: fEvent->GetEventNumberInFile();
1395 TTreeSRedirector &cstream = *fDebugStreamer;
1396 cstream<<"Transform"<<
1407 cluster->SetX(x[0]);
1408 cluster->SetY(x[1]);
1409 cluster->SetZ(x[2]);
1414 //if (!fkParam->IsGeoRead()) fkParam->ReadGeoMatrices();
1415 if (AliTPCReconstructor::GetRecoParam()->GetUseSectorAlignment() && (!calibDB->HasAlignmentOCDB())){
1416 TGeoHMatrix *mat = fkParam->GetClusterMatrix(cluster->GetDetector());
1418 Double_t pos[3]= {cluster->GetX(),cluster->GetY(),cluster->GetZ()};
1419 Double_t posC[3]={cluster->GetX(),cluster->GetY(),cluster->GetZ()};
1420 if (mat) mat->LocalToMaster(pos,posC);
1422 // chack Loading of Geo matrices from GeoManager - TEMPORARY FIX
1424 cluster->SetX(posC[0]);
1425 cluster->SetY(posC[1]);
1426 cluster->SetZ(posC[2]);
1430 //_____________________________________________________________________________
1431 Int_t AliTPCtrackerMI::LoadOuterSectors() {
1432 //-----------------------------------------------------------------
1433 // This function fills outer TPC sectors with clusters.
1434 //-----------------------------------------------------------------
1435 Int_t nrows = fOuterSec->GetNRows();
1437 for (Int_t sec = 0;sec<fkNOS;sec++)
1438 for (Int_t row = 0;row<nrows;row++){
1439 AliTPCtrackerRow* tpcrow = &(fOuterSec[sec%fkNOS][row]);
1440 Int_t sec2 = sec+2*fkNIS;
1442 Int_t ncl = tpcrow->GetN1();
1444 AliTPCclusterMI *c= (tpcrow->GetCluster1(ncl));
1445 index=(((sec2<<8)+row)<<16)+ncl;
1446 tpcrow->InsertCluster(c,index);
1449 ncl = tpcrow->GetN2();
1451 AliTPCclusterMI *c= (tpcrow->GetCluster2(ncl));
1452 index=((((sec2+fkNOS)<<8)+row)<<16)+ncl;
1453 tpcrow->InsertCluster(c,index);
1456 // write indexes for fast acces
1458 for (Int_t i=0;i<510;i++)
1459 tpcrow->SetFastCluster(i,-1);
1460 for (Int_t i=0;i<tpcrow->GetN();i++){
1461 Int_t zi = Int_t((*tpcrow)[i]->GetZ()+255.);
1462 tpcrow->SetFastCluster(zi,i); // write index
1465 for (Int_t i=0;i<510;i++){
1466 if (tpcrow->GetFastCluster(i)<0)
1467 tpcrow->SetFastCluster(i,last);
1469 last = tpcrow->GetFastCluster(i);
1478 //_____________________________________________________________________________
1479 Int_t AliTPCtrackerMI::LoadInnerSectors() {
1480 //-----------------------------------------------------------------
1481 // This function fills inner TPC sectors with clusters.
1482 //-----------------------------------------------------------------
1483 Int_t nrows = fInnerSec->GetNRows();
1485 for (Int_t sec = 0;sec<fkNIS;sec++)
1486 for (Int_t row = 0;row<nrows;row++){
1487 AliTPCtrackerRow* tpcrow = &(fInnerSec[sec%fkNIS][row]);
1490 Int_t ncl = tpcrow->GetN1();
1492 AliTPCclusterMI *c= (tpcrow->GetCluster1(ncl));
1493 index=(((sec<<8)+row)<<16)+ncl;
1494 tpcrow->InsertCluster(c,index);
1497 ncl = tpcrow->GetN2();
1499 AliTPCclusterMI *c= (tpcrow->GetCluster2(ncl));
1500 index=((((sec+fkNIS)<<8)+row)<<16)+ncl;
1501 tpcrow->InsertCluster(c,index);
1504 // write indexes for fast acces
1506 for (Int_t i=0;i<510;i++)
1507 tpcrow->SetFastCluster(i,-1);
1508 for (Int_t i=0;i<tpcrow->GetN();i++){
1509 Int_t zi = Int_t((*tpcrow)[i]->GetZ()+255.);
1510 tpcrow->SetFastCluster(zi,i); // write index
1513 for (Int_t i=0;i<510;i++){
1514 if (tpcrow->GetFastCluster(i)<0)
1515 tpcrow->SetFastCluster(i,last);
1517 last = tpcrow->GetFastCluster(i);
1529 //_________________________________________________________________________
1530 AliTPCclusterMI *AliTPCtrackerMI::GetClusterMI(Int_t index) const {
1531 //--------------------------------------------------------------------
1532 // Return pointer to a given cluster
1533 //--------------------------------------------------------------------
1534 if (index<0) return 0; // no cluster
1535 Int_t sec=(index&0xff000000)>>24;
1536 Int_t row=(index&0x00ff0000)>>16;
1537 Int_t ncl=(index&0x00007fff)>>00;
1539 const AliTPCtrackerRow * tpcrow=0;
1540 AliTPCclusterMI * clrow =0;
1542 if (sec<0 || sec>=fkNIS*4) {
1543 AliWarning(Form("Wrong sector %d",sec));
1548 AliTPCtrackerSector& tracksec = fInnerSec[sec%fkNIS];
1549 if (tracksec.GetNRows()<=row) return 0;
1550 tpcrow = &(tracksec[row]);
1551 if (tpcrow==0) return 0;
1554 if (tpcrow->GetN1()<=ncl) return 0;
1555 clrow = tpcrow->GetClusters1();
1558 if (tpcrow->GetN2()<=ncl) return 0;
1559 clrow = tpcrow->GetClusters2();
1563 AliTPCtrackerSector& tracksec = fOuterSec[(sec-fkNIS*2)%fkNOS];
1564 if (tracksec.GetNRows()<=row) return 0;
1565 tpcrow = &(tracksec[row]);
1566 if (tpcrow==0) return 0;
1568 if (sec-2*fkNIS<fkNOS) {
1569 if (tpcrow->GetN1()<=ncl) return 0;
1570 clrow = tpcrow->GetClusters1();
1573 if (tpcrow->GetN2()<=ncl) return 0;
1574 clrow = tpcrow->GetClusters2();
1578 return &(clrow[ncl]);
1584 Int_t AliTPCtrackerMI::FollowToNext(AliTPCseed& t, Int_t nr) {
1585 //-----------------------------------------------------------------
1586 // This function tries to find a track prolongation to next pad row
1587 //-----------------------------------------------------------------
1589 Double_t x= GetXrow(nr), ymax=GetMaxY(nr);
1592 AliTPCclusterMI *cl=0;
1593 Int_t tpcindex= t.GetClusterIndex2(nr);
1595 // update current shape info every 5 pad-row
1596 // if ( (nr%5==0) || t.GetNumberOfClusters()<2 || (t.fCurrentSigmaY2<0.0001) ){
1600 if (fIteration>0 && tpcindex>=-1){ //if we have already clusters
1602 if (tpcindex==-1) return 0; //track in dead zone
1603 if (tpcindex >= 0){ //
1604 cl = t.GetClusterPointer(nr);
1605 if ( (cl==0) ) cl = GetClusterMI(tpcindex);
1606 t.SetCurrentClusterIndex1(tpcindex);
1609 Int_t relativesector = ((tpcindex&0xff000000)>>24)%18; // if previously accepted cluster in different sector
1610 Float_t angle = relativesector*fSectors->GetAlpha()+fSectors->GetAlphaShift();
1612 if (angle<-TMath::Pi()) angle += 2*TMath::Pi();
1613 if (angle>=TMath::Pi()) angle -= 2*TMath::Pi();
1615 if (TMath::Abs(angle-t.GetAlpha())>0.001){
1616 Double_t rotation = angle-t.GetAlpha();
1617 t.SetRelativeSector(relativesector);
1618 if (!t.Rotate(rotation)) return 0;
1620 if (!t.PropagateTo(x)) return 0;
1622 t.SetCurrentCluster(cl);
1624 Int_t accept = AcceptCluster(&t,t.GetCurrentCluster());
1625 if ((tpcindex&0x8000)==0) accept =0;
1627 //if founded cluster is acceptible
1628 if (cl->IsUsed(11)) { // id cluster is shared inrease uncertainty
1629 t.SetErrorY2(t.GetErrorY2()+0.03);
1630 t.SetErrorZ2(t.GetErrorZ2()+0.03);
1631 t.SetErrorY2(t.GetErrorY2()*3);
1632 t.SetErrorZ2(t.GetErrorZ2()*3);
1634 t.SetNFoundable(t.GetNFoundable()+1);
1635 UpdateTrack(&t,accept);
1638 else { // Remove old cluster from track
1639 t.SetClusterIndex(nr, -3);
1640 t.SetClusterPointer(nr, 0);
1644 if (TMath::Abs(t.GetSnp())>AliTPCReconstructor::GetMaxSnpTracker()) return 0; // cut on angle
1645 if (fIteration>1 && IsFindable(t)){
1646 // not look for new cluster during refitting
1647 t.SetNFoundable(t.GetNFoundable()+1);
1652 // if (TMath::Abs(t.GetSnp())>0.95 || TMath::Abs(x*t.GetC()-t.GetEta())>0.95) return 0;// patch 28 fev 06
1653 if (!t.PropagateTo(x)) {
1654 if (fIteration==0) t.SetRemoval(10);
1657 Double_t y = t.GetY();
1658 if (TMath::Abs(y)>ymax){
1660 t.SetRelativeSector((t.GetRelativeSector()+1) % fN);
1661 if (!t.Rotate(fSectors->GetAlpha()))
1663 } else if (y <-ymax) {
1664 t.SetRelativeSector((t.GetRelativeSector()-1+fN) % fN);
1665 if (!t.Rotate(-fSectors->GetAlpha()))
1668 if (!t.PropagateTo(x)) {
1669 if (fIteration==0) t.SetRemoval(10);
1675 Double_t z=t.GetZ();
1678 if (!IsActive(t.GetRelativeSector(),nr)) {
1680 t.SetClusterIndex2(nr,-1);
1683 //AliInfo(Form("A - Sector%d phi %f - alpha %f", t.fRelativeSector,y/x, t.GetAlpha()));
1684 Bool_t isActive = IsActive(t.GetRelativeSector(),nr);
1685 Bool_t isActive2 = (nr>=fInnerSec->GetNRows()) ? fOuterSec[t.GetRelativeSector()][nr-fInnerSec->GetNRows()].GetN()>0:fInnerSec[t.GetRelativeSector()][nr].GetN()>0;
1687 if (!isActive || !isActive2) return 0;
1689 const AliTPCtrackerRow &krow=GetRow(t.GetRelativeSector(),nr);
1690 if ( (t.GetSigmaY2()<0) || t.GetSigmaZ2()<0) return 0;
1692 Double_t roadz = 1.;
1694 if (TMath::Abs(TMath::Abs(y)-ymax)<krow.GetDeadZone()){
1696 t.SetClusterIndex2(nr,-1);
1702 // if (TMath::Abs(z)<(AliTPCReconstructor::GetCtgRange()*x+10) && TMath::Abs(z)<fkParam->GetZLength(0) && (TMath::Abs(t.GetSnp())<AliTPCReconstructor::GetMaxSnpTracker()))
1703 t.SetNFoundable(t.GetNFoundable()+1);
1709 // cl = krow.FindNearest2(y+10.,z,roady,roadz,index);
1710 cl = krow.FindNearest2(y,z,roady,roadz,index);
1711 if (cl) t.SetCurrentClusterIndex1(krow.GetIndex(index));
1714 t.SetCurrentCluster(cl);
1716 if (fIteration==2&&cl->IsUsed(10)) return 0;
1717 Int_t accept = AcceptCluster(&t,t.GetCurrentCluster());
1718 if (fIteration==2&&cl->IsUsed(11)) {
1719 t.SetErrorY2(t.GetErrorY2()+0.03);
1720 t.SetErrorZ2(t.GetErrorZ2()+0.03);
1721 t.SetErrorY2(t.GetErrorY2()*3);
1722 t.SetErrorZ2(t.GetErrorZ2()*3);
1725 if (t.fCurrentCluster->IsUsed(10)){
1730 if (t.fNShared>0.7*t.GetNumberOfClusters()) {
1736 if (accept<3) UpdateTrack(&t,accept);
1739 if ( fIteration==0 && t.GetNFoundable()*0.5 > t.GetNumberOfClusters()) t.SetRemoval(10);
1747 //_________________________________________________________________________
1748 Bool_t AliTPCtrackerMI::GetTrackPoint(Int_t index, AliTrackPoint &p ) const
1750 // Get track space point by index
1751 // return false in case the cluster doesn't exist
1752 AliTPCclusterMI *cl = GetClusterMI(index);
1753 if (!cl) return kFALSE;
1754 Int_t sector = (index&0xff000000)>>24;
1755 // Int_t row = (index&0x00ff0000)>>16;
1757 // xyz[0] = fkParam->GetPadRowRadii(sector,row);
1758 xyz[0] = cl->GetX();
1759 xyz[1] = cl->GetY();
1760 xyz[2] = cl->GetZ();
1762 fkParam->AdjustCosSin(sector,cos,sin);
1763 Float_t x = cos*xyz[0]-sin*xyz[1];
1764 Float_t y = cos*xyz[1]+sin*xyz[0];
1766 Float_t sigmaY2 = 0.027*cl->GetSigmaY2();
1767 if (sector < fkParam->GetNInnerSector()) sigmaY2 *= 2.07;
1768 Float_t sigmaZ2 = 0.066*cl->GetSigmaZ2();
1769 if (sector < fkParam->GetNInnerSector()) sigmaZ2 *= 1.77;
1770 cov[0] = sin*sin*sigmaY2;
1771 cov[1] = -sin*cos*sigmaY2;
1773 cov[3] = cos*cos*sigmaY2;
1776 p.SetXYZ(x,y,xyz[2],cov);
1777 AliGeomManager::ELayerID iLayer;
1779 if (sector < fkParam->GetNInnerSector()) {
1780 iLayer = AliGeomManager::kTPC1;
1784 iLayer = AliGeomManager::kTPC2;
1785 idet = sector - fkParam->GetNInnerSector();
1787 UShort_t volid = AliGeomManager::LayerToVolUID(iLayer,idet);
1788 p.SetVolumeID(volid);
1794 Int_t AliTPCtrackerMI::UpdateClusters(AliTPCseed& t, Int_t nr) {
1795 //-----------------------------------------------------------------
1796 // This function tries to find a track prolongation to next pad row
1797 //-----------------------------------------------------------------
1798 t.SetCurrentCluster(0);
1799 t.SetCurrentClusterIndex1(-3);
1801 Double_t xt=t.GetX();
1802 Int_t row = GetRowNumber(xt)-1;
1803 Double_t ymax= GetMaxY(nr);
1805 if (row < nr) return 1; // don't prolongate if not information until now -
1806 // if (TMath::Abs(t.GetSnp())>0.9 && t.GetNumberOfClusters()>40. && fIteration!=2) {
1808 // return 0; // not prolongate strongly inclined tracks
1810 // if (TMath::Abs(t.GetSnp())>0.95) {
1812 // return 0; // not prolongate strongly inclined tracks
1813 // }// patch 28 fev 06
1815 Double_t x= GetXrow(nr);
1817 //t.PropagateTo(x+0.02);
1818 //t.PropagateTo(x+0.01);
1819 if (!t.PropagateTo(x)){
1826 if (TMath::Abs(y)>ymax){
1828 t.SetRelativeSector((t.GetRelativeSector()+1) % fN);
1829 if (!t.Rotate(fSectors->GetAlpha()))
1831 } else if (y <-ymax) {
1832 t.SetRelativeSector((t.GetRelativeSector()-1+fN) % fN);
1833 if (!t.Rotate(-fSectors->GetAlpha()))
1836 // if (!t.PropagateTo(x)){
1843 if (TMath::Abs(t.GetSnp())>AliTPCReconstructor::GetMaxSnpTracker()) return 0;
1845 if (!IsActive(t.GetRelativeSector(),nr)) {
1847 t.SetClusterIndex2(nr,-1);
1850 //AliInfo(Form("A - Sector%d phi %f - alpha %f", t.fRelativeSector,y/x, t.GetAlpha()));
1852 AliTPCtrackerRow &krow=GetRow(t.GetRelativeSector(),nr);
1854 if (TMath::Abs(TMath::Abs(y)-ymax)<krow.GetDeadZone()){
1856 t.SetClusterIndex2(nr,-1);
1862 // if (TMath::Abs(t.GetZ())<(AliTPCReconstructor::GetCtgRange()*t.GetX()+10) && (TMath::Abs(t.GetSnp())<AliTPCReconstructor::GetMaxSnpTracker()))
1863 if (IsFindable(t)) t.SetNFoundable(t.GetNFoundable()+1);
1869 if ( (nr%2==0) || t.GetNumberOfClusters()<2){
1870 // t.fCurrentSigmaY = GetSigmaY(&t);
1871 //t.fCurrentSigmaZ = GetSigmaZ(&t);
1875 AliTPCclusterMI *cl=0;
1878 Double_t roady = 1.;
1879 Double_t roadz = 1.;
1883 index = t.GetClusterIndex2(nr);
1884 if ( (index>0) && (index&0x8000)==0){
1885 cl = t.GetClusterPointer(nr);
1886 if ( (cl==0) && (index>0)) cl = GetClusterMI(index);
1887 t.SetCurrentClusterIndex1(index);
1889 t.SetCurrentCluster(cl);
1895 // if (index<0) return 0;
1896 UInt_t uindex = TMath::Abs(index);
1899 //cl = krow.FindNearest2(y+10,z,roady,roadz,uindex);
1900 cl = krow.FindNearest2(y,z,roady,roadz,uindex);
1903 if (cl) t.SetCurrentClusterIndex1(krow.GetIndex(uindex));
1904 t.SetCurrentCluster(cl);
1910 Int_t AliTPCtrackerMI::FollowToNextCluster(AliTPCseed & t, Int_t nr) {
1911 //-----------------------------------------------------------------
1912 // This function tries to find a track prolongation to next pad row
1913 //-----------------------------------------------------------------
1915 //update error according neighborhoud
1917 if (t.GetCurrentCluster()) {
1919 Int_t accept = AcceptCluster(&t,t.GetCurrentCluster());
1921 if (t.GetCurrentCluster()->IsUsed(10)){
1926 t.SetNShared(t.GetNShared()+1);
1927 if (t.GetNShared()>0.7*t.GetNumberOfClusters()) {
1932 if (fIteration>0) accept = 0;
1933 if (accept<3) UpdateTrack(&t,accept);
1937 if ( t.GetNumberOfClusters()>18 && ( (t.GetSigmaY2()+t.GetSigmaZ2())>0.16)) t.SetRemoval(10);
1938 if ( t.GetNumberOfClusters()>18 && t.GetChi2()/t.GetNumberOfClusters()>6 ) t.SetRemoval(10);
1940 if (( (t.GetNFoundable()*0.5 > t.GetNumberOfClusters()) || t.GetNoCluster()>15)) t.SetRemoval(10);
1948 //_____________________________________________________________________________
1949 Int_t AliTPCtrackerMI::FollowProlongation(AliTPCseed& t, Int_t rf, Int_t step, Bool_t fromSeeds) {
1950 //-----------------------------------------------------------------
1951 // This function tries to find a track prolongation.
1952 //-----------------------------------------------------------------
1953 Double_t xt=t.GetX();
1955 Double_t alpha=t.GetAlpha();
1956 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
1957 if (alpha < 0. ) alpha += 2.*TMath::Pi();
1959 t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha())%fN);
1961 Int_t first = GetRowNumber(xt);
1964 for (Int_t nr= first; nr>=rf; nr-=step) {
1966 if (t.GetKinkIndexes()[0]>0){
1967 for (Int_t i=0;i<3;i++){
1968 Int_t index = t.GetKinkIndexes()[i];
1969 if (index==0) break;
1970 if (index<0) continue;
1972 AliKink * kink = (AliKink*)fEvent->GetKink(index-1);
1974 printf("PROBLEM\n");
1977 Int_t kinkrow = kink->GetTPCRow0()+2+Int_t(0.5/(0.05+kink->GetAngle(2)));
1979 AliExternalTrackParam paramd(t);
1980 kink->SetDaughter(paramd);
1981 kink->SetStatus(2,5);
1988 if (nr==80) t.UpdateReference();
1989 if (nr<fInnerSec->GetNRows())
1990 fSectors = fInnerSec;
1992 fSectors = fOuterSec;
1993 if (FollowToNext(t,nr)==0)
2006 Int_t AliTPCtrackerMI::FollowBackProlongation(AliTPCseed& t, Int_t rf, Bool_t fromSeeds) {
2007 //-----------------------------------------------------------------
2008 // This function tries to find a track prolongation.
2009 //-----------------------------------------------------------------
2011 Double_t xt=t.GetX();
2012 Double_t alpha=t.GetAlpha();
2013 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
2014 if (alpha < 0. ) alpha += 2.*TMath::Pi();
2015 t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha())%fN);
2017 Int_t first = t.GetFirstPoint();
2018 Int_t ri = GetRowNumber(xt);
2022 if (first<ri) first = ri;
2024 if (first<0) first=0;
2025 for (Int_t nr=first; nr<=rf; nr++) {
2026 // if ( (TMath::Abs(t.GetSnp())>0.95)) break;//patch 28 fev 06
2027 if (t.GetKinkIndexes()[0]<0){
2028 for (Int_t i=0;i<3;i++){
2029 Int_t index = t.GetKinkIndexes()[i];
2030 if (index==0) break;
2031 if (index>0) continue;
2032 index = TMath::Abs(index);
2033 AliKink * kink = (AliKink*)fEvent->GetKink(index-1);
2035 printf("PROBLEM\n");
2038 Int_t kinkrow = kink->GetTPCRow0()-2-Int_t(0.5/(0.05+kink->GetAngle(2)));
2040 AliExternalTrackParam paramm(t);
2041 kink->SetMother(paramm);
2042 kink->SetStatus(2,1);
2049 if (nr<fInnerSec->GetNRows())
2050 fSectors = fInnerSec;
2052 fSectors = fOuterSec;
2063 Float_t AliTPCtrackerMI::OverlapFactor(AliTPCseed * s1, AliTPCseed * s2, Int_t &sum1, Int_t & sum2)
2065 // overlapping factor
2071 Float_t dz2 =(s1->GetZ() - s2->GetZ());
2074 Float_t dy2 =TMath::Abs((s1->GetY() - s2->GetY()));
2076 Float_t distance = TMath::Sqrt(dz2+dy2);
2077 if (distance>4.) return 0; // if there are far away - not overlap - to reduce combinatorics
2080 Int_t firstpoint = TMath::Min(s1->GetFirstPoint(),s2->GetFirstPoint());
2081 Int_t lastpoint = TMath::Max(s1->GetLastPoint(),s2->GetLastPoint());
2086 if (firstpoint>lastpoint) {
2087 firstpoint =lastpoint;
2092 for (Int_t i=firstpoint-1;i<lastpoint+1;i++){
2093 if (s1->GetClusterIndex2(i)>0) sum1++;
2094 if (s2->GetClusterIndex2(i)>0) sum2++;
2095 if (s1->GetClusterIndex2(i)==s2->GetClusterIndex2(i) && s1->GetClusterIndex2(i)>0) {
2099 if (sum<5) return 0;
2101 Float_t summin = TMath::Min(sum1+1,sum2+1);
2102 Float_t ratio = (sum+1)/Float_t(summin);
2106 void AliTPCtrackerMI::SignShared(AliTPCseed * s1, AliTPCseed * s2)
2110 Float_t thetaCut = 0.2;//+10.*TMath::Sqrt(s1->GetSigmaTglZ()+ s2->GetSigmaTglZ());
2111 if (TMath::Abs(s1->GetTgl()-s2->GetTgl())>thetaCut) return;
2112 Float_t minCl = TMath::Min(s1->GetNumberOfClusters(),s2->GetNumberOfClusters());
2113 Int_t cutN0 = TMath::Max(5,TMath::Nint(0.1*minCl));
2118 //Int_t firstpoint = TMath::Max(s1->GetFirstPoint(),s2->GetFirstPoint());
2119 //Int_t lastpoint = TMath::Min(s1->GetLastPoint(),s2->GetLastPoint());
2120 Int_t firstpoint = 0;
2121 Int_t lastpoint = 160;
2123 // if (firstpoint>=lastpoint-5) return;;
2125 for (Int_t i=firstpoint;i<lastpoint;i++){
2126 // if ( (s1->GetClusterIndex2(i)&0xFFFF8FFF)==(s2->GetClusterIndex2(i)&0xFFFF8FFF) && s1->GetClusterIndex2(i)>0) {
2127 if ( (s1->GetClusterIndex2(i))==(s2->GetClusterIndex2(i)) && s1->GetClusterIndex2(i)>0) {
2129 s1->SetSharedMapBit(i, kTRUE);
2130 s2->SetSharedMapBit(i, kTRUE);
2132 if (s1->GetClusterIndex2(i)>0)
2133 s1->SetClusterMapBit(i, kTRUE);
2135 if (sumshared>cutN0){
2138 for (Int_t i=firstpoint;i<lastpoint;i++){
2139 // if ( (s1->GetClusterIndex2(i)&0xFFFF8FFF)==(s2->GetClusterIndex2(i)&0xFFFF8FFF) && s1->GetClusterIndex2(i)>0) {
2140 if ( (s1->GetClusterIndex2(i))==(s2->GetClusterIndex2(i)) && s1->GetClusterIndex2(i)>0) {
2141 AliTPCTrackerPoint *p1 = s1->GetTrackPoint(i);
2142 AliTPCTrackerPoint *p2 = s2->GetTrackPoint(i);;
2143 if (s1->IsActive()&&s2->IsActive()){
2144 p1->SetShared(kTRUE);
2145 p2->SetShared(kTRUE);
2151 if (sumshared>cutN0){
2152 for (Int_t i=0;i<4;i++){
2153 if (s1->GetOverlapLabel(3*i)==0){
2154 s1->SetOverlapLabel(3*i, s2->GetLabel());
2155 s1->SetOverlapLabel(3*i+1,sumshared);
2156 s1->SetOverlapLabel(3*i+2,s2->GetUniqueID());
2160 for (Int_t i=0;i<4;i++){
2161 if (s2->GetOverlapLabel(3*i)==0){
2162 s2->SetOverlapLabel(3*i, s1->GetLabel());
2163 s2->SetOverlapLabel(3*i+1,sumshared);
2164 s2->SetOverlapLabel(3*i+2,s1->GetUniqueID());
2171 void AliTPCtrackerMI::SignShared(TObjArray * arr)
2174 //sort trackss according sectors
2176 for (Int_t i=0; i<arr->GetEntriesFast(); i++) {
2177 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2179 //if (pt) RotateToLocal(pt);
2183 arr->Sort(); // sorting according relative sectors
2184 arr->Expand(arr->GetEntries());
2187 Int_t nseed=arr->GetEntriesFast();
2188 for (Int_t i=0; i<nseed; i++) {
2189 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2191 for (Int_t j=0;j<12;j++){
2192 pt->SetOverlapLabel(j,0);
2195 for (Int_t i=0; i<nseed; i++) {
2196 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2198 if (pt->GetRemoval()>10) continue;
2199 for (Int_t j=i+1; j<nseed; j++){
2200 AliTPCseed *pt2=(AliTPCseed*)arr->UncheckedAt(j);
2201 if (TMath::Abs(pt->GetRelativeSector()-pt2->GetRelativeSector())>1) continue;
2203 if (pt2->GetRemoval()<=10) {
2204 //if ( TMath::Abs(pt->GetRelativeSector()-pt2->GetRelativeSector())>0) break;
2212 void AliTPCtrackerMI::SortTracks(TObjArray * arr, Int_t mode) const
2215 //sort tracks in array according mode criteria
2216 Int_t nseed = arr->GetEntriesFast();
2217 for (Int_t i=0; i<nseed; i++) {
2218 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2229 void AliTPCtrackerMI::RemoveUsed2(TObjArray * arr, Float_t factor1, Float_t factor2, Int_t minimal)
2232 // Loop over all tracks and remove overlaped tracks (with lower quality)
2234 // 1. Unsign clusters
2235 // 2. Sort tracks according quality
2236 // Quality is defined by the number of cluster between first and last points
2238 // 3. Loop over tracks - decreasing quality order
2239 // a.) remove - If the fraction of shared cluster less than factor (1- n or 2)
2240 // b.) remove - If the minimal number of clusters less than minimal and not ITS
2241 // c.) if track accepted - sign clusters
2243 //Called in - AliTPCtrackerMI::Clusters2Tracks()
2244 // - AliTPCtrackerMI::PropagateBack()
2245 // - AliTPCtrackerMI::RefitInward()
2248 // factor1 - factor for constrained
2249 // factor2 - for non constrained tracks
2250 // if (Float_t(shared+1)/Float_t(found+1)>factor) - DELETE
2254 Int_t nseed = arr->GetEntriesFast();
2255 Float_t * quality = new Float_t[nseed];
2256 Int_t * indexes = new Int_t[nseed];
2260 for (Int_t i=0; i<nseed; i++) {
2261 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2266 pt->UpdatePoints(); //select first last max dens points
2267 Float_t * points = pt->GetPoints();
2268 if (points[3]<0.8) quality[i] =-1;
2269 quality[i] = (points[2]-points[0])+pt->GetNumberOfClusters();
2270 //prefer high momenta tracks if overlaps
2271 quality[i] *= TMath::Sqrt(TMath::Abs(pt->Pt())+0.5);
2273 TMath::Sort(nseed,quality,indexes);
2276 for (Int_t itrack=0; itrack<nseed; itrack++) {
2277 Int_t trackindex = indexes[itrack];
2278 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(trackindex);
2281 if (quality[trackindex]<0){
2282 delete arr->RemoveAt(trackindex);
2287 Int_t first = Int_t(pt->GetPoints()[0]);
2288 Int_t last = Int_t(pt->GetPoints()[2]);
2289 Double_t factor = (pt->GetBConstrain()) ? factor1: factor2;
2291 Int_t found,foundable,shared;
2292 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
2293 // pt->GetClusterStatistic(0,160, found, foundable,shared,kFALSE);
2294 Bool_t itsgold =kFALSE;
2297 if (pt->GetESD()->GetITSclusters(dummy)>4) itsgold= kTRUE;
2301 if (Float_t(shared+1)/Float_t(found+1)>factor){
2302 if (pt->GetKinkIndexes()[0]!=0) continue; //don't remove tracks - part of the kinks
2303 if( AliTPCReconstructor::StreamLevel()>3){
2304 TTreeSRedirector &cstream = *fDebugStreamer;
2305 cstream<<"RemoveUsed"<<
2306 "iter="<<fIteration<<
2310 delete arr->RemoveAt(trackindex);
2313 if (pt->GetNumberOfClusters()<50&&(found-0.5*shared)<minimal){ //remove short tracks
2314 if (pt->GetKinkIndexes()[0]!=0) continue; //don't remove tracks - part of the kinks
2315 if( AliTPCReconstructor::StreamLevel()>3){
2316 TTreeSRedirector &cstream = *fDebugStreamer;
2317 cstream<<"RemoveShort"<<
2318 "iter="<<fIteration<<
2322 delete arr->RemoveAt(trackindex);
2328 //if (sharedfactor>0.4) continue;
2329 if (pt->GetKinkIndexes()[0]>0) continue;
2330 //Remove tracks with undefined properties - seems
2331 if (pt->GetSigmaY2()<kAlmost0) continue; // ? what is the origin ?
2333 for (Int_t i=first; i<last; i++) {
2334 Int_t index=pt->GetClusterIndex2(i);
2335 // if (index<0 || index&0x8000 ) continue;
2336 if (index<0 || index&0x8000 ) continue;
2337 AliTPCclusterMI *c= pt->GetClusterPointer(i);
2344 Info("RemoveUsed","\n*****\nNumber of good tracks after shared removal\t%d\n",fNtracks);
2350 void AliTPCtrackerMI::DumpClusters(Int_t iter, TObjArray *trackArray)
2353 // Dump clusters after reco
2354 // signed and unsigned cluster can be visualized
2355 // 1. Unsign all cluster
2356 // 2. Sign all used clusters
2359 Int_t nseed = trackArray->GetEntries();
2360 for (Int_t i=0; i<nseed; i++){
2361 AliTPCseed *pt=(AliTPCseed*)trackArray->UncheckedAt(i);
2365 Bool_t isKink=pt->GetKinkIndex(0)!=0;
2366 for (Int_t j=0; j<160; ++j) {
2367 Int_t index=pt->GetClusterIndex2(j);
2368 if (index<0) continue;
2369 AliTPCclusterMI *c= pt->GetClusterPointer(j);
2371 if (isKink) c->Use(100); // kink
2372 c->Use(10); // by default usage 10
2377 for (Int_t sec=0;sec<fkNIS;sec++){
2378 for (Int_t row=0;row<fInnerSec->GetNRows();row++){
2379 AliTPCclusterMI *cl = fInnerSec[sec][row].GetClusters1();
2380 for (Int_t icl =0;icl< fInnerSec[sec][row].GetN1();icl++){
2381 Float_t gx[3]; cl[icl].GetGlobalXYZ(gx);
2382 (*fDebugStreamer)<<"clDump"<<
2390 cl = fInnerSec[sec][row].GetClusters2();
2391 for (Int_t icl =0;icl< fInnerSec[sec][row].GetN2();icl++){
2392 Float_t gx[3]; cl[icl].GetGlobalXYZ(gx);
2393 (*fDebugStreamer)<<"clDump"<<
2404 for (Int_t sec=0;sec<fkNOS;sec++){
2405 for (Int_t row=0;row<fOuterSec->GetNRows();row++){
2406 AliTPCclusterMI *cl = fOuterSec[sec][row].GetClusters1();
2407 for (Int_t icl =0;icl< fOuterSec[sec][row].GetN1();icl++){
2408 Float_t gx[3]; cl[icl].GetGlobalXYZ(gx);
2409 (*fDebugStreamer)<<"clDump"<<
2417 cl = fOuterSec[sec][row].GetClusters2();
2418 for (Int_t icl =0;icl< fOuterSec[sec][row].GetN2();icl++){
2419 Float_t gx[3]; cl[icl].GetGlobalXYZ(gx);
2420 (*fDebugStreamer)<<"clDump"<<
2432 void AliTPCtrackerMI::UnsignClusters()
2435 // loop over all clusters and unsign them
2438 for (Int_t sec=0;sec<fkNIS;sec++){
2439 for (Int_t row=0;row<fInnerSec->GetNRows();row++){
2440 AliTPCclusterMI *cl = fInnerSec[sec][row].GetClusters1();
2441 for (Int_t icl =0;icl< fInnerSec[sec][row].GetN1();icl++)
2442 // if (cl[icl].IsUsed(10))
2444 cl = fInnerSec[sec][row].GetClusters2();
2445 for (Int_t icl =0;icl< fInnerSec[sec][row].GetN2();icl++)
2446 //if (cl[icl].IsUsed(10))
2451 for (Int_t sec=0;sec<fkNOS;sec++){
2452 for (Int_t row=0;row<fOuterSec->GetNRows();row++){
2453 AliTPCclusterMI *cl = fOuterSec[sec][row].GetClusters1();
2454 for (Int_t icl =0;icl< fOuterSec[sec][row].GetN1();icl++)
2455 //if (cl[icl].IsUsed(10))
2457 cl = fOuterSec[sec][row].GetClusters2();
2458 for (Int_t icl =0;icl< fOuterSec[sec][row].GetN2();icl++)
2459 //if (cl[icl].IsUsed(10))
2468 void AliTPCtrackerMI::SignClusters(const TObjArray * arr, Float_t fnumber, Float_t fdensity)
2471 //sign clusters to be "used"
2473 // snumber and sdensity sign number of sigmas - bellow mean value to be accepted
2474 // loop over "primaries"
2488 Int_t nseed = arr->GetEntriesFast();
2489 for (Int_t i=0; i<nseed; i++) {
2490 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2494 if (!(pt->IsActive())) continue;
2495 Float_t dens = pt->GetNumberOfClusters()/Float_t(pt->GetNFoundable());
2496 if ( (dens>0.7) && (pt->GetNumberOfClusters()>70)){
2498 sumdens2+= dens*dens;
2499 sumn += pt->GetNumberOfClusters();
2500 sumn2 += pt->GetNumberOfClusters()*pt->GetNumberOfClusters();
2501 Float_t chi2 = pt->GetChi2()/pt->GetNumberOfClusters();
2504 sumchi2 +=chi2*chi2;
2509 Float_t mdensity = 0.9;
2510 Float_t meann = 130;
2511 Float_t meanchi = 1;
2512 Float_t sdensity = 0.1;
2513 Float_t smeann = 10;
2514 Float_t smeanchi =0.4;
2518 mdensity = sumdens/sum;
2520 meanchi = sumchi/sum;
2522 sdensity = sumdens2/sum-mdensity*mdensity;
2524 sdensity = TMath::Sqrt(sdensity);
2528 smeann = sumn2/sum-meann*meann;
2530 smeann = TMath::Sqrt(smeann);
2534 smeanchi = sumchi2/sum - meanchi*meanchi;
2536 smeanchi = TMath::Sqrt(smeanchi);
2542 //REMOVE SHORT DELTAS or tracks going out of sensitive volume of TPC
2544 for (Int_t i=0; i<nseed; i++) {
2545 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2549 if (pt->GetBSigned()) continue;
2550 if (pt->GetBConstrain()) continue;
2551 //if (!(pt->IsActive())) continue;
2553 Int_t found,foundable,shared;
2554 pt->GetClusterStatistic(0,160,found, foundable,shared);
2555 if (shared/float(found)>0.3) {
2556 if (shared/float(found)>0.9 ){
2557 //delete arr->RemoveAt(i);
2562 Bool_t isok =kFALSE;
2563 if ( (pt->GetNShared()/pt->GetNumberOfClusters()<0.5) &&pt->GetNumberOfClusters()>60)
2565 if ((TMath::Abs(1/pt->GetC())<100.) && (pt->GetNShared()/pt->GetNumberOfClusters()<0.7))
2567 if (TMath::Abs(pt->GetZ()/pt->GetX())>1.1)
2569 if ( (TMath::Abs(pt->GetSnp()>0.7) && pt->GetD(0,0)>60.))
2573 for (Int_t j=0; j<160; ++j) {
2574 Int_t index=pt->GetClusterIndex2(j);
2575 if (index<0) continue;
2576 AliTPCclusterMI *c= pt->GetClusterPointer(j);
2578 //if (!(c->IsUsed(10))) c->Use();
2585 Double_t maxchi = meanchi+2.*smeanchi;
2587 for (Int_t i=0; i<nseed; i++) {
2588 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2592 //if (!(pt->IsActive())) continue;
2593 if (pt->GetBSigned()) continue;
2594 Double_t chi = pt->GetChi2()/pt->GetNumberOfClusters();
2595 if (chi>maxchi) continue;
2598 Float_t dens = pt->GetNumberOfClusters()/Float_t(pt->GetNFoundable());
2600 //sign only tracks with enoug big density at the beginning
2602 if ((pt->GetDensityFirst(40)<0.75) && pt->GetNumberOfClusters()<meann) continue;
2605 Double_t mindens = TMath::Max(double(mdensity-sdensity*fdensity*bfactor),0.65);
2606 Double_t minn = TMath::Max(Int_t(meann-fnumber*smeann*bfactor),50);
2608 // if (pt->fBConstrain) mindens = TMath::Max(mdensity-sdensity*fdensity*bfactor,0.65);
2609 if ( (pt->GetRemoval()==10) && (pt->GetSnp()>0.8)&&(dens>mindens))
2612 if ((dens>mindens && pt->GetNumberOfClusters()>minn) && chi<maxchi ){
2613 //Int_t noc=pt->GetNumberOfClusters();
2614 pt->SetBSigned(kTRUE);
2615 for (Int_t j=0; j<160; ++j) {
2617 Int_t index=pt->GetClusterIndex2(j);
2618 if (index<0) continue;
2619 AliTPCclusterMI *c= pt->GetClusterPointer(j);
2621 // if (!(c->IsUsed(10))) c->Use();
2626 // gLastCheck = nseed;
2634 void AliTPCtrackerMI::StopNotActive(const TObjArray * arr, Int_t row0, Float_t th0, Float_t th1, Float_t th2) const
2636 // stop not active tracks
2637 // take th1 as threshold for number of founded to number of foundable on last 10 active rows
2638 // take th2 as threshold for number of founded to number of foundable on last 20 active rows
2639 Int_t nseed = arr->GetEntriesFast();
2641 for (Int_t i=0; i<nseed; i++) {
2642 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2646 if (!(pt->IsActive())) continue;
2647 StopNotActive(pt,row0,th0, th1,th2);
2653 void AliTPCtrackerMI::StopNotActive(AliTPCseed * seed, Int_t row0, Float_t th0, Float_t th1,
2656 // stop not active tracks
2657 // take th1 as threshold for number of founded to number of foundable on last 10 active rows
2658 // take th2 as threshold for number of founded to number of foundable on last 20 active rows
2661 Int_t foundable = 0;
2662 Int_t maxindex = seed->GetLastPoint(); //last foundable row
2663 if (seed->GetNFoundable()*th0 > seed->GetNumberOfClusters()) {
2664 seed->Desactivate(10) ;
2668 for (Int_t i=row0; i<maxindex; i++){
2669 Int_t index = seed->GetClusterIndex2(i);
2670 if (index!=-1) foundable++;
2672 if (foundable<=30) sumgood1++;
2673 if (foundable<=50) {
2680 if (foundable>=30.){
2681 if (sumgood1<(th1*30.)) seed->Desactivate(10);
2684 if (sumgood2<(th2*50.)) seed->Desactivate(10);
2688 Int_t AliTPCtrackerMI::RefitInward(AliESDEvent *event)
2691 // back propagation of ESD tracks
2694 if (!event) return 0;
2695 const Int_t kMaxFriendTracks=2000;
2697 // extract correction object for multiplicity dependence of dEdx
2698 TObjArray * gainCalibArray = AliTPCcalibDB::Instance()->GetTimeGainSplinesRun(event->GetRunNumber());
2700 AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ;
2702 AliFatal("Tranformations not in RefitInward");
2705 transform->SetCurrentRecoParam((AliTPCRecoParam*)AliTPCReconstructor::GetRecoParam());
2706 const AliTPCRecoParam * recoParam = AliTPCcalibDB::Instance()->GetTransform()->GetCurrentRecoParam();
2707 Int_t nContribut = event->GetNumberOfTracks();
2708 TGraphErrors * graphMultDependenceDeDx = 0x0;
2709 if (recoParam && recoParam->GetUseMultiplicityCorrectionDedx() && gainCalibArray) {
2710 if (recoParam->GetUseTotCharge()) {
2711 graphMultDependenceDeDx = (TGraphErrors *) gainCalibArray->FindObject("TGRAPHERRORS_MEANQTOT_MULTIPLICITYDEPENDENCE_BEAM_ALL");
2713 graphMultDependenceDeDx = (TGraphErrors *) gainCalibArray->FindObject("TGRAPHERRORS_MEANQMAX_MULTIPLICITYDEPENDENCE_BEAM_ALL");
2719 //PrepareForProlongation(fSeeds,1);
2720 PropagateForward2(fSeeds);
2721 RemoveUsed2(fSeeds,0.4,0.4,20);
2723 TObjArray arraySeed(fSeeds->GetEntries());
2724 for (Int_t i=0;i<fSeeds->GetEntries();i++) {
2725 arraySeed.AddAt(fSeeds->At(i),i);
2727 SignShared(&arraySeed);
2728 // FindCurling(fSeeds, event,2); // find multi found tracks
2729 FindSplitted(fSeeds, event,2); // find multi found tracks
2730 if (AliTPCReconstructor::StreamLevel()>5) FindMultiMC(fSeeds, fEvent,2); // find multi found tracks
2733 Int_t nseed = fSeeds->GetEntriesFast();
2734 for (Int_t i=0;i<nseed;i++){
2735 AliTPCseed * seed = (AliTPCseed*) fSeeds->UncheckedAt(i);
2736 if (!seed) continue;
2737 if (seed->GetKinkIndex(0)>0) UpdateKinkQualityD(seed); // update quality informations for kinks
2738 AliESDtrack *esd=event->GetTrack(i);
2740 if (seed->GetNumberOfClusters()<60 && seed->GetNumberOfClusters()<(esd->GetTPCclusters(0) -5)*0.8){
2741 AliInfo("Refitting track");
2742 AliExternalTrackParam paramIn;
2743 AliExternalTrackParam paramOut;
2744 Int_t ncl = seed->RefitTrack(seed,¶mIn,¶mOut);
2745 if (AliTPCReconstructor::StreamLevel()>2) {
2746 (*fDebugStreamer)<<"RecoverIn"<<
2750 "pout.="<<¶mOut<<
2755 seed->Set(paramIn.GetX(),paramIn.GetAlpha(),paramIn.GetParameter(),paramIn.GetCovariance());
2756 seed->SetNumberOfClusters(ncl);
2760 seed->PropagateTo(fkParam->GetInnerRadiusLow());
2761 seed->UpdatePoints();
2762 AddCovariance(seed);
2764 seed->CookdEdx(0.02,0.6);
2765 CookLabel(seed,0.1); //For comparison only
2766 esd->SetTPCClusterMap(seed->GetClusterMap());
2767 esd->SetTPCSharedMap(seed->GetSharedMap());
2769 if (AliTPCReconstructor::StreamLevel()>1 && seed!=0) {
2770 TTreeSRedirector &cstream = *fDebugStreamer;
2777 if (seed->GetNumberOfClusters()>15){
2778 esd->UpdateTrackParams(seed,AliESDtrack::kTPCrefit);
2779 esd->SetTPCPoints(seed->GetPoints());
2780 esd->SetTPCPointsF(seed->GetNFoundable());
2781 Int_t ndedx = seed->GetNCDEDX(0);
2782 Float_t sdedx = seed->GetSDEDX(0);
2783 Float_t dedx = seed->GetdEdx();
2784 // apply mutliplicity dependent dEdx correction if available
2785 if (graphMultDependenceDeDx) {
2786 Double_t corrGain = AliTPCcalibDButil::EvalGraphConst(graphMultDependenceDeDx, nContribut);
2787 dedx += (1 - corrGain)*50.; // MIP is normalized to 50
2789 esd->SetTPCsignal(dedx, sdedx, ndedx);
2791 // add seed to the esd track in Calib level
2793 Bool_t storeFriend = gRandom->Rndm()<(kMaxFriendTracks)/Float_t(nseed);
2794 if (AliTPCReconstructor::StreamLevel()>0 &&storeFriend){
2795 AliTPCseed * seedCopy = new AliTPCseed(*seed, kTRUE);
2796 esd->AddCalibObject(seedCopy);
2801 //printf("problem\n");
2804 //FindKinks(fSeeds,event);
2805 if (AliTPCReconstructor::StreamLevel()>3) DumpClusters(2,fSeeds);
2806 Info("RefitInward","Number of refitted tracks %d",ntracks);
2811 Int_t AliTPCtrackerMI::PropagateBack(AliESDEvent *event)
2814 // back propagation of ESD tracks
2816 if (!event) return 0;
2820 PropagateBack(fSeeds);
2821 RemoveUsed2(fSeeds,0.4,0.4,20);
2822 //FindCurling(fSeeds, fEvent,1);
2823 FindSplitted(fSeeds, event,1); // find multi found tracks
2824 if (AliTPCReconstructor::StreamLevel()>5) FindMultiMC(fSeeds, fEvent,1); // find multi found tracks
2827 Int_t nseed = fSeeds->GetEntriesFast();
2829 for (Int_t i=0;i<nseed;i++){
2830 AliTPCseed * seed = (AliTPCseed*) fSeeds->UncheckedAt(i);
2831 if (!seed) continue;
2832 if (seed->GetKinkIndex(0)<0) UpdateKinkQualityM(seed); // update quality informations for kinks
2833 seed->UpdatePoints();
2834 AddCovariance(seed);
2835 AliESDtrack *esd=event->GetTrack(i);
2836 if (!esd) continue; //never happen
2837 if (seed->GetNumberOfClusters()<60 && seed->GetNumberOfClusters()<(esd->GetTPCclusters(0) -5)*0.8){
2838 AliExternalTrackParam paramIn;
2839 AliExternalTrackParam paramOut;
2840 Int_t ncl = seed->RefitTrack(seed,¶mIn,¶mOut);
2841 if (AliTPCReconstructor::StreamLevel()>2) {
2842 (*fDebugStreamer)<<"RecoverBack"<<
2846 "pout.="<<¶mOut<<
2851 seed->Set(paramOut.GetX(),paramOut.GetAlpha(),paramOut.GetParameter(),paramOut.GetCovariance());
2852 seed->SetNumberOfClusters(ncl);
2855 seed->CookdEdx(0.02,0.6);
2856 CookLabel(seed,0.1); //For comparison only
2857 if (seed->GetNumberOfClusters()>15){
2858 esd->UpdateTrackParams(seed,AliESDtrack::kTPCout);
2859 esd->SetTPCPoints(seed->GetPoints());
2860 esd->SetTPCPointsF(seed->GetNFoundable());
2861 Int_t ndedx = seed->GetNCDEDX(0);
2862 Float_t sdedx = seed->GetSDEDX(0);
2863 Float_t dedx = seed->GetdEdx();
2864 esd->SetTPCsignal(dedx, sdedx, ndedx);
2866 Int_t eventnumber = event->GetEventNumberInFile();// patch 28 fev 06
2867 // This is most likely NOT the event number you'd like to use. It has nothing to do with the 'real' event number
2868 if (AliTPCReconstructor::StreamLevel()>1 && esd) {
2869 (*fDebugStreamer)<<"Cback"<<
2872 "EventNrInFile="<<eventnumber<<
2877 if (AliTPCReconstructor::StreamLevel()>3) DumpClusters(1,fSeeds);
2878 //FindKinks(fSeeds,event);
2879 Info("PropagateBack","Number of back propagated tracks %d",ntracks);
2886 void AliTPCtrackerMI::DeleteSeeds()
2891 Int_t nseed = fSeeds->GetEntriesFast();
2892 for (Int_t i=0;i<nseed;i++){
2893 AliTPCseed * seed = (AliTPCseed*)fSeeds->At(i);
2894 if (seed) delete fSeeds->RemoveAt(i);
2901 void AliTPCtrackerMI::ReadSeeds(const AliESDEvent *const event, Int_t direction)
2904 //read seeds from the event
2906 Int_t nentr=event->GetNumberOfTracks();
2908 Info("ReadSeeds", "Number of ESD tracks: %d\n", nentr);
2913 fSeeds = new TObjArray(nentr);
2917 for (Int_t i=0; i<nentr; i++) {
2918 AliESDtrack *esd=event->GetTrack(i);
2919 ULong_t status=esd->GetStatus();
2920 if (!(status&AliESDtrack::kTPCin)) continue;
2921 AliTPCtrack t(*esd);
2922 t.SetNumberOfClusters(0);
2923 // AliTPCseed *seed = new AliTPCseed(t,t.GetAlpha());
2924 AliTPCseed *seed = new AliTPCseed(t/*,t.GetAlpha()*/);
2925 seed->SetUniqueID(esd->GetID());
2926 AddCovariance(seed); //add systematic ucertainty
2927 for (Int_t ikink=0;ikink<3;ikink++) {
2928 Int_t index = esd->GetKinkIndex(ikink);
2929 seed->GetKinkIndexes()[ikink] = index;
2930 if (index==0) continue;
2931 index = TMath::Abs(index);
2932 AliESDkink * kink = fEvent->GetKink(index-1);
2933 if (kink&&esd->GetKinkIndex(ikink)<0){
2934 if ((status & AliESDtrack::kTRDrefit) != 0) kink->SetStatus(1,2);
2935 if ((status & AliESDtrack::kITSout) != 0) kink->SetStatus(1,0);
2937 if (kink&&esd->GetKinkIndex(ikink)>0){
2938 if ((status & AliESDtrack::kTRDrefit) != 0) kink->SetStatus(1,6);
2939 if ((status & AliESDtrack::kITSout) != 0) kink->SetStatus(1,4);
2943 if (((status&AliESDtrack::kITSout)==0)&&(direction==1)) seed->ResetCovariance(10.);
2944 if ( direction ==2 &&(status & AliESDtrack::kTRDrefit) == 0 ) seed->ResetCovariance(10.);
2945 //if ( direction ==2 && ((status & AliESDtrack::kTPCout) == 0) ) {
2946 // fSeeds->AddAt(0,i);
2950 if ( direction ==2 &&(status & AliESDtrack::kTRDrefit) > 0 ) {
2951 Double_t par0[5],par1[5],alpha,x;
2952 esd->GetInnerExternalParameters(alpha,x,par0);
2953 esd->GetExternalParameters(x,par1);
2954 Double_t delta1 = TMath::Abs(par0[4]-par1[4])/(0.000000001+TMath::Abs(par0[4]+par1[4]));
2955 Double_t delta2 = TMath::Abs(par0[3]-par1[3]);
2957 if (esd->GetTRDncls()>0) trdchi2 = esd->GetTRDchi2()/esd->GetTRDncls();
2958 //reset covariance if suspicious
2959 if ( (delta1>0.1) || (delta2>0.006) ||trdchi2>7.)
2960 seed->ResetCovariance(10.);
2965 // rotate to the local coordinate system
2967 fSectors=fInnerSec; fN=fkNIS;
2968 Double_t alpha=seed->GetAlpha();
2969 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
2970 if (alpha < 0. ) alpha += 2.*TMath::Pi();
2971 Int_t ns=Int_t(alpha/fSectors->GetAlpha())%fN;
2972 alpha =ns*fSectors->GetAlpha() + fSectors->GetAlphaShift();
2973 alpha-=seed->GetAlpha();
2974 if (alpha<-TMath::Pi()) alpha += 2*TMath::Pi();
2975 if (alpha>=TMath::Pi()) alpha -= 2*TMath::Pi();
2976 if (TMath::Abs(alpha) > 0.001) { // This should not happen normally
2977 AliWarning(Form("Rotating track over %f",alpha));
2978 if (!seed->Rotate(alpha)) {
2985 if (esd->GetKinkIndex(0)<=0){
2986 for (Int_t irow=0;irow<160;irow++){
2987 Int_t index = seed->GetClusterIndex2(irow);
2990 AliTPCclusterMI * cl = GetClusterMI(index);
2991 seed->SetClusterPointer(irow,cl);
2993 if ((index & 0x8000)==0){
2994 cl->Use(10); // accepted cluster
2996 cl->Use(6); // close cluster not accepted
2999 Info("ReadSeeds","Not found cluster");
3004 fSeeds->AddAt(seed,i);
3010 //_____________________________________________________________________________
3011 void AliTPCtrackerMI::MakeSeeds3(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, Float_t cuts[4],
3012 Float_t deltay, Int_t ddsec) {
3013 //-----------------------------------------------------------------
3014 // This function creates track seeds.
3015 // SEEDING WITH VERTEX CONSTRAIN
3016 //-----------------------------------------------------------------
3017 // cuts[0] - fP4 cut
3018 // cuts[1] - tan(phi) cut
3019 // cuts[2] - zvertex cut
3020 // cuts[3] - fP3 cut
3028 Double_t x[5], c[15];
3029 // Int_t di = i1-i2;
3031 AliTPCseed * seed = new AliTPCseed();
3032 Double_t alpha=fSectors->GetAlpha(), shift=fSectors->GetAlphaShift();
3033 Double_t cs=cos(alpha), sn=sin(alpha);
3035 // Double_t x1 =fOuterSec->GetX(i1);
3036 //Double_t xx2=fOuterSec->GetX(i2);
3038 Double_t x1 =GetXrow(i1);
3039 Double_t xx2=GetXrow(i2);
3041 Double_t x3=GetX(), y3=GetY(), z3=GetZ();
3043 Int_t imiddle = (i2+i1)/2; //middle pad row index
3044 Double_t xm = GetXrow(imiddle); // radius of middle pad-row
3045 const AliTPCtrackerRow& krm=GetRow(sec,imiddle); //middle pad -row
3049 const AliTPCtrackerRow& kr1=GetRow(ns,i1);
3050 Double_t ymax = GetMaxY(i1)-kr1.GetDeadZone()-1.5;
3051 Double_t ymaxm = GetMaxY(imiddle)-kr1.GetDeadZone()-1.5;
3054 // change cut on curvature if it can't reach this layer
3055 // maximal curvature set to reach it
3056 Double_t dvertexmax = TMath::Sqrt((x1-x3)*(x1-x3)+(ymax+5-y3)*(ymax+5-y3));
3057 if (dvertexmax*0.5*cuts[0]>0.85){
3058 cuts[0] = 0.85/(dvertexmax*0.5+1.);
3060 Double_t r2min = 1/(cuts[0]*cuts[0]); //minimal square of radius given by cut
3063 if (deltay>0) ddsec = 0;
3064 // loop over clusters
3065 for (Int_t is=0; is < kr1; is++) {
3067 if (kr1[is]->IsUsed(10)) continue;
3068 Double_t y1=kr1[is]->GetY(), z1=kr1[is]->GetZ();
3069 //if (TMath::Abs(y1)>ymax) continue;
3071 if (deltay>0 && TMath::Abs(ymax-TMath::Abs(y1))> deltay ) continue; // seed only at the edge
3073 // find possible directions
3074 Float_t anglez = (z1-z3)/(x1-x3);
3075 Float_t extraz = z1 - anglez*(x1-xx2); // extrapolated z
3078 //find rotation angles relative to line given by vertex and point 1
3079 Double_t dvertex2 = (x1-x3)*(x1-x3)+(y1-y3)*(y1-y3);
3080 Double_t dvertex = TMath::Sqrt(dvertex2);
3081 Double_t angle13 = TMath::ATan((y1-y3)/(x1-x3));
3082 Double_t cs13 = cos(-angle13), sn13 = sin(-angle13);
3085 // loop over 2 sectors
3091 Double_t dddz1=0; // direction of delta inclination in z axis
3098 for (Int_t dsec = dsec1; dsec<=dsec2;dsec++){
3099 Int_t sec2 = sec + dsec;
3101 // AliTPCtrackerRow& kr2 = fOuterSec[(sec2+fkNOS)%fkNOS][i2];
3102 //AliTPCtrackerRow& kr2m = fOuterSec[(sec2+fkNOS)%fkNOS][imiddle];
3103 AliTPCtrackerRow& kr2 = GetRow((sec2+fkNOS)%fkNOS,i2);
3104 AliTPCtrackerRow& kr2m = GetRow((sec2+fkNOS)%fkNOS,imiddle);
3105 Int_t index1 = TMath::Max(kr2.Find(extraz-0.6-dddz1*TMath::Abs(z1)*0.05)-1,0);
3106 Int_t index2 = TMath::Min(kr2.Find(extraz+0.6+dddz2*TMath::Abs(z1)*0.05)+1,kr2);
3108 // rotation angles to p1-p3
3109 Double_t cs13r = cos(-angle13+dsec*alpha)/dvertex, sn13r = sin(-angle13+dsec*alpha)/dvertex;
3110 Double_t x2, y2, z2;
3112 // Double_t dymax = maxangle*TMath::Abs(x1-xx2);
3115 Double_t dxx0 = (xx2-x3)*cs13r;
3116 Double_t dyy0 = (xx2-x3)*sn13r;
3117 for (Int_t js=index1; js < index2; js++) {
3118 const AliTPCclusterMI *kcl = kr2[js];
3119 if (kcl->IsUsed(10)) continue;
3121 //calcutate parameters
3123 Double_t yy0 = dyy0 +(kcl->GetY()-y3)*cs13r;
3125 if (TMath::Abs(yy0)<0.000001) continue;
3126 Double_t xx0 = dxx0 -(kcl->GetY()-y3)*sn13r;
3127 Double_t y0 = 0.5*(xx0*xx0+yy0*yy0-xx0)/yy0;
3128 Double_t r02 = (0.25+y0*y0)*dvertex2;
3129 //curvature (radius) cut
3130 if (r02<r2min) continue;
3134 Double_t c0 = 1/TMath::Sqrt(r02);
3138 //Double_t dfi0 = 2.*TMath::ASin(dvertex*c0*0.5);
3139 //Double_t dfi1 = 2.*TMath::ASin(TMath::Sqrt(yy0*yy0+(1-xx0)*(1-xx0))*dvertex*c0*0.5);
3140 Double_t dfi0 = 2.*AliTPCFastMath::FastAsin(dvertex*c0*0.5);
3141 Double_t dfi1 = 2.*AliTPCFastMath::FastAsin(TMath::Sqrt(yy0*yy0+(1-xx0)*(1-xx0))*dvertex*c0*0.5);
3144 Double_t z0 = kcl->GetZ();
3145 Double_t zzzz2 = z1-(z1-z3)*dfi1/dfi0;
3146 if (TMath::Abs(zzzz2-z0)>0.5) continue;
3149 Double_t dip = (z1-z0)*c0/dfi1;
3150 Double_t x0 = (0.5*cs13+y0*sn13)*dvertex*c0;
3161 x2= xx2*cs-y2*sn*dsec;
3162 y2=+xx2*sn*dsec+y2*cs;
3172 // do we have cluster at the middle ?
3174 GetProlongation(x1,xm,x,ym,zm);
3176 AliTPCclusterMI * cm=0;
3177 if (TMath::Abs(ym)-ymaxm<0){
3178 cm = krm.FindNearest2(ym,zm,1.0,0.6,dummy);
3179 if ((!cm) || (cm->IsUsed(10))) {
3184 // rotate y1 to system 0
3185 // get state vector in rotated system
3186 Double_t yr1 = (-0.5*sn13+y0*cs13)*dvertex*c0;
3187 Double_t xr2 = x0*cs+yr1*sn*dsec;
3188 Double_t xr[5]={kcl->GetY(),kcl->GetZ(), xr2, dip, c0};
3190 GetProlongation(xx2,xm,xr,ym,zm);
3191 if (TMath::Abs(ym)-ymaxm<0){
3192 cm = kr2m.FindNearest2(ym,zm,1.0,0.6,dummy);
3193 if ((!cm) || (cm->IsUsed(10))) {
3203 dym = ym - cm->GetY();
3204 dzm = zm - cm->GetZ();
3211 Double_t sy1=kr1[is]->GetSigmaY2()*2., sz1=kr1[is]->GetSigmaZ2()*2.;
3212 Double_t sy2=kcl->GetSigmaY2()*2., sz2=kcl->GetSigmaZ2()*2.;
3213 //Double_t sy3=400*3./12., sy=0.1, sz=0.1;
3214 Double_t sy3=25000*x[4]*x[4]+0.1, sy=0.1, sz=0.1;
3215 //Double_t sy3=25000*x[4]*x[4]*60+0.5, sy=0.1, sz=0.1;
3217 Double_t f40=(F1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
3218 Double_t f42=(F1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
3219 Double_t f43=(F1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
3220 Double_t f20=(F2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
3221 Double_t f22=(F2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
3222 Double_t f23=(F2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
3224 Double_t f30=(F3(x1,y1+sy,x2,y2,z1,z2)-x[3])/sy;
3225 Double_t f31=(F3(x1,y1,x2,y2,z1+sz,z2)-x[3])/sz;
3226 Double_t f32=(F3(x1,y1,x2,y2+sy,z1,z2)-x[3])/sy;
3227 Double_t f34=(F3(x1,y1,x2,y2,z1,z2+sz)-x[3])/sz;
3231 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3232 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3233 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3234 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3235 c[13]=f30*sy1*f40+f32*sy2*f42;
3236 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3238 // if (!BuildSeed(kr1[is],kcl,0,x1,x2,x3,x,c)) continue;
3240 UInt_t index=kr1.GetIndex(is);
3241 seed->~AliTPCseed(); // this does not set the pointer to 0...
3242 AliTPCseed *track=new(seed) AliTPCseed(x1, ns*alpha+shift, x, c, index);
3244 track->SetIsSeeding(kTRUE);
3245 track->SetSeed1(i1);
3246 track->SetSeed2(i2);
3247 track->SetSeedType(3);
3251 FollowProlongation(*track, (i1+i2)/2,1);
3252 Int_t foundable,found,shared;
3253 track->GetClusterStatistic((i1+i2)/2,i1, found, foundable, shared, kTRUE);
3254 if ((found<0.55*foundable) || shared>0.5*found || (track->GetSigmaY2()+track->GetSigmaZ2())>0.5){
3256 seed->~AliTPCseed();
3262 FollowProlongation(*track, i2,1);
3266 track->SetBConstrain(1);
3267 // track->fLastPoint = i1+fInnerSec->GetNRows(); // first cluster in track position
3268 track->SetLastPoint(i1); // first cluster in track position
3269 track->SetFirstPoint(track->GetLastPoint());
3271 if (track->GetNumberOfClusters()<(i1-i2)*0.5 ||
3272 track->GetNumberOfClusters() < track->GetNFoundable()*0.6 ||
3273 track->GetNShared()>0.4*track->GetNumberOfClusters() ) {
3275 seed->~AliTPCseed();
3279 // Z VERTEX CONDITION
3280 Double_t zv, bz=GetBz();
3281 if ( !track->GetZAt(0.,bz,zv) ) continue;
3282 if (TMath::Abs(zv-z3)>cuts[2]) {
3283 FollowProlongation(*track, TMath::Max(i2-20,0));
3284 if ( !track->GetZAt(0.,bz,zv) ) continue;
3285 if (TMath::Abs(zv-z3)>cuts[2]){
3286 FollowProlongation(*track, TMath::Max(i2-40,0));
3287 if ( !track->GetZAt(0.,bz,zv) ) continue;
3288 if (TMath::Abs(zv-z3)>cuts[2] &&(track->GetNumberOfClusters() > track->GetNFoundable()*0.7)){
3289 // make seed without constrain
3290 AliTPCseed * track2 = MakeSeed(track,0.2,0.5,1.);
3291 FollowProlongation(*track2, i2,1);
3292 track2->SetBConstrain(kFALSE);
3293 track2->SetSeedType(1);
3294 arr->AddLast(track2);
3296 seed->~AliTPCseed();
3301 seed->~AliTPCseed();
3308 track->SetSeedType(0);
3309 arr->AddLast(track);
3310 seed = new AliTPCseed;
3312 // don't consider other combinations
3313 if (track->GetNumberOfClusters() > track->GetNFoundable()*0.8)
3319 Info("MakeSeeds3","\nSeeding statistic:\t%d\t%d\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin,nout1,nout2);
3325 void AliTPCtrackerMI::MakeSeeds5(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, Float_t cuts[4],
3330 //-----------------------------------------------------------------
3331 // This function creates track seeds.
3332 //-----------------------------------------------------------------
3333 // cuts[0] - fP4 cut
3334 // cuts[1] - tan(phi) cut
3335 // cuts[2] - zvertex cut
3336 // cuts[3] - fP3 cut
3346 Double_t x[5], c[15];
3348 // make temporary seed
3349 AliTPCseed * seed = new AliTPCseed;
3350 Double_t alpha=fOuterSec->GetAlpha(), shift=fOuterSec->GetAlphaShift();
3351 // Double_t cs=cos(alpha), sn=sin(alpha);
3356 Double_t x1 = GetXrow(i1-1);
3357 const AliTPCtrackerRow& kr1=GetRow(sec,i1-1);
3358 Double_t y1max = GetMaxY(i1-1)-kr1.GetDeadZone()-1.5;
3360 Double_t x1p = GetXrow(i1);
3361 const AliTPCtrackerRow& kr1p=GetRow(sec,i1);
3363 Double_t x1m = GetXrow(i1-2);
3364 const AliTPCtrackerRow& kr1m=GetRow(sec,i1-2);
3367 //last 3 padrow for seeding
3368 AliTPCtrackerRow& kr3 = GetRow((sec+fkNOS)%fkNOS,i1-7);
3369 Double_t x3 = GetXrow(i1-7);
3370 // Double_t y3max= GetMaxY(i1-7)-kr3.fDeadZone-1.5;
3372 AliTPCtrackerRow& kr3p = GetRow((sec+fkNOS)%fkNOS,i1-6);
3373 Double_t x3p = GetXrow(i1-6);
3375 AliTPCtrackerRow& kr3m = GetRow((sec+fkNOS)%fkNOS,i1-8);
3376 Double_t x3m = GetXrow(i1-8);
3381 Int_t im = i1-4; //middle pad row index
3382 Double_t xm = GetXrow(im); // radius of middle pad-row
3383 const AliTPCtrackerRow& krm=GetRow(sec,im); //middle pad -row
3384 // Double_t ymmax = GetMaxY(im)-kr1.fDeadZone-1.5;
3387 Double_t deltax = x1-x3;
3388 Double_t dymax = deltax*cuts[1];
3389 Double_t dzmax = deltax*cuts[3];
3391 // loop over clusters
3392 for (Int_t is=0; is < kr1; is++) {
3394 if (kr1[is]->IsUsed(10)) continue;
3395 Double_t y1=kr1[is]->GetY(), z1=kr1[is]->GetZ();
3397 if (deltay>0 && TMath::Abs(y1max-TMath::Abs(y1))> deltay ) continue; // seed only at the edge
3399 Int_t index1 = TMath::Max(kr3.Find(z1-dzmax)-1,0);
3400 Int_t index2 = TMath::Min(kr3.Find(z1+dzmax)+1,kr3);
3406 for (Int_t js=index1; js < index2; js++) {
3407 const AliTPCclusterMI *kcl = kr3[js];
3408 if (kcl->IsUsed(10)) continue;
3410 // apply angular cuts
3411 if (TMath::Abs(y1-y3)>dymax) continue;
3414 if (TMath::Abs(z1-z3)>dzmax) continue;
3416 Double_t angley = (y1-y3)/(x1-x3);
3417 Double_t anglez = (z1-z3)/(x1-x3);
3419 Double_t erry = TMath::Abs(angley)*(x1-x1m)*0.5+0.5;
3420 Double_t errz = TMath::Abs(anglez)*(x1-x1m)*0.5+0.5;
3422 Double_t yyym = angley*(xm-x1)+y1;
3423 Double_t zzzm = anglez*(xm-x1)+z1;
3425 const AliTPCclusterMI *kcm = krm.FindNearest2(yyym,zzzm,erry,errz,index);
3427 if (kcm->IsUsed(10)) continue;
3429 erry = TMath::Abs(angley)*(x1-x1m)*0.4+0.5;
3430 errz = TMath::Abs(anglez)*(x1-x1m)*0.4+0.5;
3437 // look around first
3438 const AliTPCclusterMI *kc1m = kr1m.FindNearest2(angley*(x1m-x1)+y1,
3444 if (kc1m->IsUsed(10)) used++;
3446 const AliTPCclusterMI *kc1p = kr1p.FindNearest2(angley*(x1p-x1)+y1,
3452 if (kc1p->IsUsed(10)) used++;
3454 if (used>1) continue;
3455 if (found<1) continue;
3459 const AliTPCclusterMI *kc3m = kr3m.FindNearest2(angley*(x3m-x3)+y3,
3465 if (kc3m->IsUsed(10)) used++;
3469 const AliTPCclusterMI *kc3p = kr3p.FindNearest2(angley*(x3p-x3)+y3,
3475 if (kc3p->IsUsed(10)) used++;
3479 if (used>1) continue;
3480 if (found<3) continue;
3490 x[4]=F1(x1,y1,x2,y2,x3,y3);
3491 //if (TMath::Abs(x[4]) >= cuts[0]) continue;
3494 x[2]=F2(x1,y1,x2,y2,x3,y3);
3497 x[3]=F3n(x1,y1,x2,y2,z1,z2,x[4]);
3498 //if (TMath::Abs(x[3]) > cuts[3]) continue;
3502 Double_t sy1=0.1, sz1=0.1;
3503 Double_t sy2=0.1, sz2=0.1;
3504 Double_t sy3=0.1, sy=0.1, sz=0.1;
3506 Double_t f40=(F1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
3507 Double_t f42=(F1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
3508 Double_t f43=(F1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
3509 Double_t f20=(F2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
3510 Double_t f22=(F2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
3511 Double_t f23=(F2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
3513 Double_t f30=(F3(x1,y1+sy,x2,y2,z1,z2)-x[3])/sy;
3514 Double_t f31=(F3(x1,y1,x2,y2,z1+sz,z2)-x[3])/sz;
3515 Double_t f32=(F3(x1,y1,x2,y2+sy,z1,z2)-x[3])/sy;
3516 Double_t f34=(F3(x1,y1,x2,y2,z1,z2+sz)-x[3])/sz;
3520 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3521 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3522 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3523 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3524 c[13]=f30*sy1*f40+f32*sy2*f42;
3525 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3527 // if (!BuildSeed(kr1[is],kcl,0,x1,x2,x3,x,c)) continue;
3529 index=kr1.GetIndex(is);
3530 seed->~AliTPCseed();
3531 AliTPCseed *track=new(seed) AliTPCseed(x1, sec*alpha+shift, x, c, index);
3533 track->SetIsSeeding(kTRUE);
3536 FollowProlongation(*track, i1-7,1);
3537 if (track->GetNumberOfClusters() < track->GetNFoundable()*0.75 ||
3538 track->GetNShared()>0.6*track->GetNumberOfClusters() || ( track->GetSigmaY2()+ track->GetSigmaZ2())>0.6){
3540 seed->~AliTPCseed();
3546 FollowProlongation(*track, i2,1);
3547 track->SetBConstrain(0);
3548 track->SetLastPoint(i1+fInnerSec->GetNRows()); // first cluster in track position
3549 track->SetFirstPoint(track->GetLastPoint());
3551 if (track->GetNumberOfClusters()<(i1-i2)*0.5 ||
3552 track->GetNumberOfClusters()<track->GetNFoundable()*0.7 ||
3553 track->GetNShared()>2. || track->GetChi2()/track->GetNumberOfClusters()>6 || ( track->GetSigmaY2()+ track->GetSigmaZ2())>0.5 ) {
3555 seed->~AliTPCseed();
3560 FollowProlongation(*track, TMath::Max(i2-10,0),1);
3561 AliTPCseed * track2 = MakeSeed(track,0.2,0.5,0.9);
3562 FollowProlongation(*track2, i2,1);
3563 track2->SetBConstrain(kFALSE);
3564 track2->SetSeedType(4);
3565 arr->AddLast(track2);
3567 seed->~AliTPCseed();
3571 //arr->AddLast(track);
3572 //seed = new AliTPCseed;
3578 Info("MakeSeeds5","\nSeeding statiistic:\t%d\t%d\t%d\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin,nout1,nout2,nout3);
3584 //_____________________________________________________________________________
3585 void AliTPCtrackerMI::MakeSeeds2(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, Float_t */*cuts[4]*/,
3586 Float_t deltay, Bool_t /*bconstrain*/) {
3587 //-----------------------------------------------------------------
3588 // This function creates track seeds - without vertex constraint
3589 //-----------------------------------------------------------------
3590 // cuts[0] - fP4 cut - not applied
3591 // cuts[1] - tan(phi) cut
3592 // cuts[2] - zvertex cut - not applied
3593 // cuts[3] - fP3 cut
3603 Double_t alpha=fOuterSec->GetAlpha(), shift=fOuterSec->GetAlphaShift();
3604 // Double_t cs=cos(alpha), sn=sin(alpha);
3605 Int_t row0 = (i1+i2)/2;
3606 Int_t drow = (i1-i2)/2;
3607 const AliTPCtrackerRow& kr0=fSectors[sec][row0];
3608 AliTPCtrackerRow * kr=0;
3610 AliTPCpolyTrack polytrack;
3611 Int_t nclusters=fSectors[sec][row0];
3612 AliTPCseed * seed = new AliTPCseed;
3617 for (Int_t is=0; is < nclusters; is++) { //LOOP over clusters
3619 Int_t nfoundable =0;
3620 for (Int_t iter =1; iter<2; iter++){ //iterations
3621 const AliTPCtrackerRow& krm=fSectors[sec][row0-iter];
3622 const AliTPCtrackerRow& krp=fSectors[sec][row0+iter];
3623 const AliTPCclusterMI * cl= kr0[is];
3625 if (cl->IsUsed(10)) {
3631 Double_t x = kr0.GetX();
3632 // Initialization of the polytrack
3637 Double_t y0= cl->GetY();
3638 Double_t z0= cl->GetZ();
3642 Double_t ymax = fSectors->GetMaxY(row0)-kr0.GetDeadZone()-1.5;
3643 if (deltay>0 && TMath::Abs(ymax-TMath::Abs(y0))> deltay ) continue; // seed only at the edge
3645 erry = (0.5)*cl->GetSigmaY2()/TMath::Sqrt(cl->GetQ())*6;
3646 errz = (0.5)*cl->GetSigmaZ2()/TMath::Sqrt(cl->GetQ())*6;
3647 polytrack.AddPoint(x,y0,z0,erry, errz);
3650 if (cl->IsUsed(10)) sumused++;
3653 Float_t roady = (5*TMath::Sqrt(cl->GetSigmaY2()+0.2)+1.)*iter;
3654 Float_t roadz = (5*TMath::Sqrt(cl->GetSigmaZ2()+0.2)+1.)*iter;
3657 AliTPCclusterMI * cl1 = krm.FindNearest(y0,z0,roady,roadz);
3658 if (cl1 && TMath::Abs(ymax-TMath::Abs(y0))) {
3659 erry = (0.5)*cl1->GetSigmaY2()/TMath::Sqrt(cl1->GetQ())*3;
3660 errz = (0.5)*cl1->GetSigmaZ2()/TMath::Sqrt(cl1->GetQ())*3;
3661 if (cl1->IsUsed(10)) sumused++;
3662 polytrack.AddPoint(x,cl1->GetY(),cl1->GetZ(),erry,errz);
3666 AliTPCclusterMI * cl2 = krp.FindNearest(y0,z0,roady,roadz);
3668 erry = (0.5)*cl2->GetSigmaY2()/TMath::Sqrt(cl2->GetQ())*3;
3669 errz = (0.5)*cl2->GetSigmaZ2()/TMath::Sqrt(cl2->GetQ())*3;
3670 if (cl2->IsUsed(10)) sumused++;
3671 polytrack.AddPoint(x,cl2->GetY(),cl2->GetZ(),erry,errz);
3674 if (sumused>0) continue;
3676 polytrack.UpdateParameters();
3682 nfoundable = polytrack.GetN();
3683 nfound = nfoundable;
3685 for (Int_t ddrow = iter+1; ddrow<drow;ddrow++){
3686 Float_t maxdist = 0.8*(1.+3./(ddrow));
3687 for (Int_t delta = -1;delta<=1;delta+=2){
3688 Int_t row = row0+ddrow*delta;
3689 kr = &(fSectors[sec][row]);
3690 Double_t xn = kr->GetX();
3691 Double_t ymax1 = fSectors->GetMaxY(row)-kr->GetDeadZone()-1.5;
3692 polytrack.GetFitPoint(xn,yn,zn);
3693 if (TMath::Abs(yn)>ymax1) continue;
3695 AliTPCclusterMI * cln = kr->FindNearest(yn,zn,roady,roadz);
3697 Float_t dist = TMath::Sqrt( (yn-cln->GetY())*(yn-cln->GetY())+(zn-cln->GetZ())*(zn-cln->GetZ()));
3700 erry = (dist+0.3)*cln->GetSigmaY2()/TMath::Sqrt(cln->GetQ())*(1.+1./(ddrow));
3701 errz = (dist+0.3)*cln->GetSigmaZ2()/TMath::Sqrt(cln->GetQ())*(1.+1./(ddrow));
3702 if (cln->IsUsed(10)) {
3703 // printf("used\n");
3711 polytrack.AddPoint(xn,cln->GetY(),cln->GetZ(),erry, errz);
3716 if ( (sumused>3) || (sumused>0.5*nfound) || (nfound<0.6*nfoundable)) break;
3717 polytrack.UpdateParameters();
3720 if ( (sumused>3) || (sumused>0.5*nfound)) {
3721 //printf("sumused %d\n",sumused);
3726 polytrack.GetFitDerivation(kr0.GetX(),dy,dz);
3727 AliTPCpolyTrack track2;
3729 polytrack.Refit(track2,0.5+TMath::Abs(dy)*0.3,0.4+TMath::Abs(dz)*0.3);
3730 if (track2.GetN()<0.5*nfoundable) continue;
3733 if ((nfound>0.6*nfoundable) &&( nfoundable>0.4*(i1-i2))) {
3735 // test seed with and without constrain
3736 for (Int_t constrain=0; constrain<=0;constrain++){
3737 // add polytrack candidate
3739 Double_t x[5], c[15];
3740 Double_t x1,x2,x3,y1,y2,y3,z1,z2,z3;
3741 track2.GetBoundaries(x3,x1);
3743 track2.GetFitPoint(x1,y1,z1);
3744 track2.GetFitPoint(x2,y2,z2);
3745 track2.GetFitPoint(x3,y3,z3);
3747 //is track pointing to the vertex ?
3750 polytrack.GetFitPoint(x0,y0,z0);
3763 x[4]=F1(x1,y1,x2,y2,x3,y3);
3765 // if (TMath::Abs(x[4]) >= cuts[0]) continue; //
3766 x[2]=F2(x1,y1,x2,y2,x3,y3);
3768 //if (TMath::Abs(x[4]*x1-x[2]) >= cuts[1]) continue;
3769 //x[3]=F3(x1,y1,x2,y2,z1,z2);
3770 x[3]=F3n(x1,y1,x3,y3,z1,z3,x[4]);
3771 //if (TMath::Abs(x[3]) > cuts[3]) continue;
3774 Double_t sy =0.1, sz =0.1;
3775 Double_t sy1=0.02, sz1=0.02;
3776 Double_t sy2=0.02, sz2=0.02;
3780 sy3=25000*x[4]*x[4]+0.1, sy=0.1, sz=0.1;
3783 Double_t f40=(F1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
3784 Double_t f42=(F1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
3785 Double_t f43=(F1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
3786 Double_t f20=(F2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
3787 Double_t f22=(F2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
3788 Double_t f23=(F2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
3790 Double_t f30=(F3(x1,y1+sy,x3,y3,z1,z3)-x[3])/sy;
3791 Double_t f31=(F3(x1,y1,x3,y3,z1+sz,z3)-x[3])/sz;
3792 Double_t f32=(F3(x1,y1,x3,y3+sy,z1,z3)-x[3])/sy;
3793 Double_t f34=(F3(x1,y1,x3,y3,z1,z3+sz)-x[3])/sz;
3798 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3799 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3800 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3801 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3802 c[13]=f30*sy1*f40+f32*sy2*f42;
3803 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3805 //Int_t row1 = fSectors->GetRowNumber(x1);
3806 Int_t row1 = GetRowNumber(x1);
3810 seed->~AliTPCseed();
3811 AliTPCseed *track=new(seed) AliTPCseed(x1,sec*alpha+shift,x,c,index);
3812 track->SetIsSeeding(kTRUE);
3813 Int_t rc=FollowProlongation(*track, i2);
3814 if (constrain) track->SetBConstrain(1);
3816 track->SetBConstrain(0);
3817 track->SetLastPoint(row1+fInnerSec->GetNRows()); // first cluster in track position
3818 track->SetFirstPoint(track->GetLastPoint());
3820 if (rc==0 || track->GetNumberOfClusters()<(i1-i2)*0.5 ||
3821 track->GetNumberOfClusters() < track->GetNFoundable()*0.6 ||
3822 track->GetNShared()>0.4*track->GetNumberOfClusters()) {
3825 seed->~AliTPCseed();
3828 arr->AddLast(track);
3829 seed = new AliTPCseed;
3833 } // if accepted seed
3836 Info("MakeSeeds2","\nSeeding statiistic:\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin3);
3842 AliTPCseed *AliTPCtrackerMI::MakeSeed(AliTPCseed *const track, Float_t r0, Float_t r1, Float_t r2)
3846 //reseed using track points
3847 Int_t p0 = int(r0*track->GetNumberOfClusters()); // point 0
3848 Int_t p1 = int(r1*track->GetNumberOfClusters());
3849 Int_t p2 = int(r2*track->GetNumberOfClusters()); // last point
3851 Double_t x0[3],x1[3],x2[3];
3852 for (Int_t i=0;i<3;i++){
3858 // find track position at given ratio of the length
3859 Int_t sec0=0, sec1=0, sec2=0;
3862 for (Int_t i=0;i<160;i++){
3863 if (track->GetClusterPointer(i)){
3865 AliTPCTrackerPoint *trpoint =track->GetTrackPoint(i);
3866 if ( (index<p0) || x0[0]<0 ){
3867 if (trpoint->GetX()>1){
3868 clindex = track->GetClusterIndex2(i);
3870 x0[0] = trpoint->GetX();
3871 x0[1] = trpoint->GetY();
3872 x0[2] = trpoint->GetZ();
3873 sec0 = ((clindex&0xff000000)>>24)%18;
3878 if ( (index<p1) &&(trpoint->GetX()>1)){
3879 clindex = track->GetClusterIndex2(i);
3881 x1[0] = trpoint->GetX();
3882 x1[1] = trpoint->GetY();
3883 x1[2] = trpoint->GetZ();
3884 sec1 = ((clindex&0xff000000)>>24)%18;
3887 if ( (index<p2) &&(trpoint->GetX()>1)){
3888 clindex = track->GetClusterIndex2(i);
3890 x2[0] = trpoint->GetX();
3891 x2[1] = trpoint->GetY();
3892 x2[2] = trpoint->GetZ();
3893 sec2 = ((clindex&0xff000000)>>24)%18;
3900 Double_t alpha, cs,sn, xx2,yy2;
3902 alpha = (sec1-sec2)*fSectors->GetAlpha();
3903 cs = TMath::Cos(alpha);
3904 sn = TMath::Sin(alpha);
3905 xx2= x1[0]*cs-x1[1]*sn;
3906 yy2= x1[0]*sn+x1[1]*cs;
3910 alpha = (sec0-sec2)*fSectors->GetAlpha();
3911 cs = TMath::Cos(alpha);
3912 sn = TMath::Sin(alpha);
3913 xx2= x0[0]*cs-x0[1]*sn;
3914 yy2= x0[0]*sn+x0[1]*cs;
3920 Double_t x[5],c[15];
3924 x[4]=F1(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]);
3925 // if (x[4]>1) return 0;
3926 x[2]=F2(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]);
3927 x[3]=F3n(x2[0],x2[1],x0[0],x0[1],x2[2],x0[2],x[4]);
3928 //if (TMath::Abs(x[3]) > 2.2) return 0;
3929 //if (TMath::Abs(x[2]) > 1.99) return 0;
3931 Double_t sy =0.1, sz =0.1;
3933 Double_t sy1=0.02+track->GetSigmaY2(), sz1=0.02+track->GetSigmaZ2();
3934 Double_t sy2=0.01+track->GetSigmaY2(), sz2=0.01+track->GetSigmaZ2();
3935 Double_t sy3=0.01+track->GetSigmaY2();
3937 Double_t f40=(F1(x2[0],x2[1]+sy,x1[0],x1[1],x0[0],x0[1])-x[4])/sy;
3938 Double_t f42=(F1(x2[0],x2[1],x1[0],x1[1]+sy,x0[0],x0[1])-x[4])/sy;
3939 Double_t f43=(F1(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]+sy)-x[4])/sy;
3940 Double_t f20=(F2(x2[0],x2[1]+sy,x1[0],x1[1],x0[0],x0[1])-x[2])/sy;
3941 Double_t f22=(F2(x2[0],x2[1],x1[0],x1[1]+sy,x0[0],x0[1])-x[2])/sy;
3942 Double_t f23=(F2(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]+sy)-x[2])/sy;
3944 Double_t f30=(F3(x2[0],x2[1]+sy,x0[0],x0[1],x2[2],x0[2])-x[3])/sy;
3945 Double_t f31=(F3(x2[0],x2[1],x0[0],x0[1],x2[2]+sz,x0[2])-x[3])/sz;
3946 Double_t f32=(F3(x2[0],x2[1],x0[0],x0[1]+sy,x2[2],x0[2])-x[3])/sy;
3947 Double_t f34=(F3(x2[0],x2[1],x0[0],x0[1],x2[2],x0[2]+sz)-x[3])/sz;
3952 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3953 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3954 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3955 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3956 c[13]=f30*sy1*f40+f32*sy2*f42;
3957 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3959 // Int_t row1 = fSectors->GetRowNumber(x2[0]);
3960 AliTPCseed *seed=new AliTPCseed(x2[0], sec2*fSectors->GetAlpha()+fSectors->GetAlphaShift(), x, c, 0);
3961 // Double_t y0,z0,y1,z1, y2,z2;
3962 //seed->GetProlongation(x0[0],y0,z0);
3963 // seed->GetProlongation(x1[0],y1,z1);
3964 //seed->GetProlongation(x2[0],y2,z2);
3966 seed->SetLastPoint(pp2);
3967 seed->SetFirstPoint(pp2);
3974 AliTPCseed *AliTPCtrackerMI::ReSeed(const AliTPCseed *track, Float_t r0, Float_t r1, Float_t r2)
3978 //reseed using founded clusters
3980 // Find the number of clusters
3981 Int_t nclusters = 0;
3982 for (Int_t irow=0;irow<160;irow++){
3983 if (track->GetClusterIndex(irow)>0) nclusters++;
3987 ipos[0] = TMath::Max(int(r0*nclusters),0); // point 0 cluster
3988 ipos[1] = TMath::Min(int(r1*nclusters),nclusters-1); //
3989 ipos[2] = TMath::Min(int(r2*nclusters),nclusters-1); // last point
3992 Double_t xyz[3][3]={{0}};
3993 Int_t row[3]={0},sec[3]={0,0,0};
3995 // find track row position at given ratio of the length
3997 for (Int_t irow=0;irow<160;irow++){
3998 if (track->GetClusterIndex2(irow)<0) continue;
4000 for (Int_t ipoint=0;ipoint<3;ipoint++){
4001 if (index<=ipos[ipoint]) row[ipoint] = irow;
4005 //Get cluster and sector position
4006 for (Int_t ipoint=0;ipoint<3;ipoint++){
4007 Int_t clindex = track->GetClusterIndex2(row[ipoint]);
4008 AliTPCclusterMI * cl = GetClusterMI(clindex);
4011 // AliTPCclusterMI * cl = GetClusterMI(clindex);
4014 sec[ipoint] = ((clindex&0xff000000)>>24)%18;
4015 xyz[ipoint][0] = GetXrow(row[ipoint]);
4016 xyz[ipoint][1] = cl->GetY();
4017 xyz[ipoint][2] = cl->GetZ();
4021 // Calculate seed state vector and covariance matrix
4023 Double_t alpha, cs,sn, xx2,yy2;
4025 alpha = (sec[1]-sec[2])*fSectors->GetAlpha();
4026 cs = TMath::Cos(alpha);
4027 sn = TMath::Sin(alpha);
4028 xx2= xyz[1][0]*cs-xyz[1][1]*sn;
4029 yy2= xyz[1][0]*sn+xyz[1][1]*cs;
4033 alpha = (sec[0]-sec[2])*fSectors->GetAlpha();
4034 cs = TMath::Cos(alpha);
4035 sn = TMath::Sin(alpha);
4036 xx2= xyz[0][0]*cs-xyz[0][1]*sn;
4037 yy2= xyz[0][0]*sn+xyz[0][1]*cs;
4043 Double_t x[5],c[15];
4047 x[4]=F1(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
4048 x[2]=F2(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
4049 x[3]=F3n(xyz[2][0],xyz[2][1],xyz[0][0],xyz[0][1],xyz[2][2],xyz[0][2],x[4]);
4051 Double_t sy =0.1, sz =0.1;
4053 Double_t sy1=0.2, sz1=0.2;
4054 Double_t sy2=0.2, sz2=0.2;
4057 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;
4058 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;
4059 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;
4060 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;
4061 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;
4062 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;
4064 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;
4065 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;
4066 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;
4067 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;
4072 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
4073 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
4074 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
4075 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
4076 c[13]=f30*sy1*f40+f32*sy2*f42;
4077 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
4079 // Int_t row1 = fSectors->GetRowNumber(xyz[2][0]);
4080 AliTPCseed *seed=new AliTPCseed(xyz[2][0], sec[2]*fSectors->GetAlpha()+fSectors->GetAlphaShift(), x, c, 0);
4081 seed->SetLastPoint(row[2]);
4082 seed->SetFirstPoint(row[2]);
4087 AliTPCseed *AliTPCtrackerMI::ReSeed(AliTPCseed *track,Int_t r0, Bool_t forward)
4091 //reseed using founded clusters
4094 Int_t row[3]={0,0,0};
4095 Int_t sec[3]={0,0,0};
4097 // forward direction
4099 for (Int_t irow=r0;irow<160;irow++){
4100 if (track->GetClusterIndex(irow)>0){
4105 for (Int_t irow=160;irow>r0;irow--){
4106 if (track->GetClusterIndex(irow)>0){
4111 for (Int_t irow=row[2]-15;irow>row[0];irow--){
4112 if (track->GetClusterIndex(irow)>0){
4120 for (Int_t irow=0;irow<r0;irow++){
4121 if (track->GetClusterIndex(irow)>0){
4126 for (Int_t irow=r0;irow>0;irow--){
4127 if (track->GetClusterIndex(irow)>0){
4132 for (Int_t irow=row[2]-15;irow>row[0];irow--){
4133 if (track->GetClusterIndex(irow)>0){
4140 if ((row[2]-row[0])<20) return 0;
4141 if (row[1]==0) return 0;
4144 //Get cluster and sector position
4145 for (Int_t ipoint=0;ipoint<3;ipoint++){
4146 Int_t clindex = track->GetClusterIndex2(row[ipoint]);
4147 AliTPCclusterMI * cl = GetClusterMI(clindex);
4150 // AliTPCclusterMI * cl = GetClusterMI(clindex);
4153 sec[ipoint] = ((clindex&0xff000000)>>24)%18;
4154 xyz[ipoint][0] = GetXrow(row[ipoint]);
4155 AliTPCTrackerPoint * point = track->GetTrackPoint(row[ipoint]);
4156 if (point&&ipoint<2){
4158 xyz[ipoint][1] = point->GetY();
4159 xyz[ipoint][2] = point->GetZ();
4162 xyz[ipoint][1] = cl->GetY();
4163 xyz[ipoint][2] = cl->GetZ();
4170 // Calculate seed state vector and covariance matrix
4172 Double_t alpha, cs,sn, xx2,yy2;
4174 alpha = (sec[1]-sec[2])*fSectors->GetAlpha();
4175 cs = TMath::Cos(alpha);
4176 sn = TMath::Sin(alpha);
4177 xx2= xyz[1][0]*cs-xyz[1][1]*sn;
4178 yy2= xyz[1][0]*sn+xyz[1][1]*cs;
4182 alpha = (sec[0]-sec[2])*fSectors->GetAlpha();
4183 cs = TMath::Cos(alpha);
4184 sn = TMath::Sin(alpha);
4185 xx2= xyz[0][0]*cs-xyz[0][1]*sn;
4186 yy2= xyz[0][0]*sn+xyz[0][1]*cs;
4192 Double_t x[5],c[15];
4196 x[4]=F1(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
4197 x[2]=F2(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
4198 x[3]=F3n(xyz[2][0],xyz[2][1],xyz[0][0],xyz[0][1],xyz[2][2],xyz[0][2],x[4]);
4200 Double_t sy =0.1, sz =0.1;
4202 Double_t sy1=0.2, sz1=0.2;
4203 Double_t sy2=0.2, sz2=0.2;
4206 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;
4207 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;
4208 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;
4209 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;
4210 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;
4211 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;
4213 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;
4214 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;
4215 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;
4216 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;
4221 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
4222 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
4223 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
4224 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
4225 c[13]=f30*sy1*f40+f32*sy2*f42;
4226 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
4228 // Int_t row1 = fSectors->GetRowNumber(xyz[2][0]);
4229 AliTPCseed *seed=new AliTPCseed(xyz[2][0], sec[2]*fSectors->GetAlpha()+fSectors->GetAlphaShift(), x, c, 0);
4230 seed->SetLastPoint(row[2]);
4231 seed->SetFirstPoint(row[2]);
4232 for (Int_t i=row[0];i<row[2];i++){
4233 seed->SetClusterIndex(i, track->GetClusterIndex(i));
4241 void AliTPCtrackerMI::FindMultiMC(const TObjArray * array, AliESDEvent */*esd*/, Int_t iter)
4244 // find multi tracks - THIS FUNCTION IS ONLY FOR DEBUG PURPOSES
4246 // Use AliTPCReconstructor::StreamLevel()>2 if you want to tune parameters - cuts
4248 // Two reasons to have multiple find tracks
4249 // 1. Curling tracks can be find more than once
4250 // 2. Splitted tracks
4251 // a.) Multiple seeding to increase tracking efficiency - (~ 100% reached)
4252 // b.) Edge effect on the sector boundaries
4255 // Algorithm done in 2 phases - because of CPU consumption
4256 // it is n^2 algorithm - for lead-lead 20000x20000 combination are investigated
4258 // Algorihm for curling tracks sign:
4259 // 1 phase -makes a very rough fast cuts to minimize combinatorics
4260 // a.) opposite sign
4261 // b.) one of the tracks - not pointing to the primary vertex -
4262 // c.) delta tan(theta)
4264 // 2 phase - calculates DCA between tracks - time consument
4269 // General cuts - for splitted tracks and for curling tracks
4271 const Float_t kMaxdPhi = 0.2; // maximal distance in phi
4273 // Curling tracks cuts
4278 Int_t nentries = array->GetEntriesFast();
4279 AliHelix *helixes = new AliHelix[nentries];
4280 Float_t *xm = new Float_t[nentries];
4281 Float_t *dz0 = new Float_t[nentries];
4282 Float_t *dz1 = new Float_t[nentries];
4288 // Find track COG in x direction - point with best defined parameters
4290 for (Int_t i=0;i<nentries;i++){
4291 AliTPCseed* track = (AliTPCseed*)array->At(i);
4292 if (!track) continue;
4293 track->SetCircular(0);
4294 new (&helixes[i]) AliHelix(*track);
4298 track->GetDZ(GetX(),GetY(),GetZ(),GetBz(),dz);
4301 for (Int_t icl=0; icl<160; icl++){
4302 AliTPCclusterMI * cl = track->GetClusterPointer(icl);
4308 if (ncl>0) xm[i]/=Float_t(ncl);
4311 for (Int_t i0=0;i0<nentries;i0++){
4312 AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
4313 if (!track0) continue;
4314 Float_t xc0 = helixes[i0].GetHelix(6);
4315 Float_t yc0 = helixes[i0].GetHelix(7);
4316 Float_t r0 = helixes[i0].GetHelix(8);
4317 Float_t rc0 = TMath::Sqrt(xc0*xc0+yc0*yc0);
4318 Float_t fi0 = TMath::ATan2(yc0,xc0);
4320 for (Int_t i1=i0+1;i1<nentries;i1++){
4321 AliTPCseed * track1 = (AliTPCseed*)array->At(i1);
4322 if (!track1) continue;
4323 Int_t lab0=track0->GetLabel();
4324 Int_t lab1=track1->GetLabel();
4325 if (TMath::Abs(lab0)!=TMath::Abs(lab1)) continue;
4327 Float_t xc1 = helixes[i1].GetHelix(6);
4328 Float_t yc1 = helixes[i1].GetHelix(7);
4329 Float_t r1 = helixes[i1].GetHelix(8);
4330 Float_t rc1 = TMath::Sqrt(xc1*xc1+yc1*yc1);
4331 Float_t fi1 = TMath::ATan2(yc1,xc1);
4333 Float_t dfi = fi0-fi1;
4336 if (dfi>1.5*TMath::Pi()) dfi-=TMath::Pi(); // take care about edge effect
4337 if (dfi<-1.5*TMath::Pi()) dfi+=TMath::Pi(); //
4338 if (TMath::Abs(dfi)>kMaxdPhi&&helixes[i0].GetHelix(4)*helixes[i1].GetHelix(4)<0){
4340 // if short tracks with undefined sign
4341 fi1 = -TMath::ATan2(yc1,-xc1);
4344 Float_t dtheta = TMath::Abs(track0->GetTgl()-track1->GetTgl())<TMath::Abs(track0->GetTgl()+track1->GetTgl())? track0->GetTgl()-track1->GetTgl():track0->GetTgl()+track1->GetTgl();
4347 // debug stream to tune "fast cuts"
4349 Double_t dist[3]; // distance at X
4350 Double_t mdist[3]={0,0,0}; // mean distance X+-40cm
4351 track0->GetDistance(track1,0.5*(xm[i0]+xm[i1])-40.,dist,AliTracker::GetBz());
4352 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[i]);
4353 track0->GetDistance(track1,0.5*(xm[i0]+xm[i1])+40.,dist,AliTracker::GetBz());
4354 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[i]);
4355 track0->GetDistance(track1,0.5*(xm[i0]+xm[i1]),dist,AliTracker::GetBz());
4356 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[i]);
4357 for (Int_t i=0;i<3;i++) mdist[i]*=0.33333;
4361 for (Int_t icl=0; icl<160; icl++){
4362 AliTPCclusterMI * cl0 = track0->GetClusterPointer(icl);
4363 AliTPCclusterMI * cl1 = track1->GetClusterPointer(icl);
4366 if (cl0==cl1) sums++;
4370 if (AliTPCReconstructor::StreamLevel()>5) {
4371 TTreeSRedirector &cstream = *fDebugStreamer;
4376 "Tr0.="<<track0<< // seed0
4377 "Tr1.="<<track1<< // seed1
4378 "h0.="<<&helixes[i0]<<
4379 "h1.="<<&helixes[i1]<<
4381 "sum="<<sum<< //the sum of rows with cl in both
4382 "sums="<<sums<< //the sum of shared clusters
4383 "xm0="<<xm[i0]<< // the center of track
4384 "xm1="<<xm[i1]<< // the x center of track
4385 // General cut variables
4386 "dfi="<<dfi<< // distance in fi angle
4387 "dtheta="<<dtheta<< // distance int theta angle
4393 "dist0="<<dist[0]<< //distance x
4394 "dist1="<<dist[1]<< //distance y
4395 "dist2="<<dist[2]<< //distance z
4396 "mdist0="<<mdist[0]<< //distance x
4397 "mdist1="<<mdist[1]<< //distance y
4398 "mdist2="<<mdist[2]<< //distance z
4414 if (AliTPCReconstructor::StreamLevel()>1) {
4415 AliInfo("Time for curling tracks removal DEBUGGING MC");
4422 void AliTPCtrackerMI::FindSplitted(TObjArray * array, AliESDEvent */*esd*/, Int_t /*iter*/){
4424 // Find Splitted tracks and remove the one with worst quality
4425 // Corresponding debug streamer to tune selections - "Splitted2"
4427 // 0. Sort tracks according quility
4428 // 1. Propagate the tracks to the reference radius
4429 // 2. Double_t loop to select close tracks (only to speed up process)
4430 // 3. Calculate cluster overlap ratio - and remove the track if bigger than a threshold
4431 // 4. Delete temporary parameters
4433 const Double_t xref=GetXrow(63); // reference radius -IROC/OROC boundary
4435 const Double_t kCutP1=10; // delta Z cut 10 cm
4436 const Double_t kCutP2=0.15; // delta snp(fi) cut 0.15
4437 const Double_t kCutP3=0.15; // delta tgl(theta) cut 0.15
4438 const Double_t kCutAlpha=0.15; // delta alpha cut
4439 Int_t firstpoint = 0;
4440 Int_t lastpoint = 160;
4442 Int_t nentries = array->GetEntriesFast();
4443 AliExternalTrackParam *params = new AliExternalTrackParam[nentries];
4449 //0. Sort tracks according quality
4450 //1. Propagate the ext. param to reference radius
4451 Int_t nseed = array->GetEntriesFast();
4452 if (nseed<=0) return;
4453 Float_t * quality = new Float_t[nseed];
4454 Int_t * indexes = new Int_t[nseed];
4455 for (Int_t i=0; i<nseed; i++) {
4456 AliTPCseed *pt=(AliTPCseed*)array->UncheckedAt(i);
4461 pt->UpdatePoints(); //select first last max dens points
4462 Float_t * points = pt->GetPoints();
4463 if (points[3]<0.8) quality[i] =-1;
4464 quality[i] = (points[2]-points[0])+pt->GetNumberOfClusters();
4465 //prefer high momenta tracks if overlaps
4466 quality[i] *= TMath::Sqrt(TMath::Abs(pt->Pt())+0.5);
4468 AliTracker::PropagateTrackToBxByBz(&(params[i]),xref,pt->GetMass(),5.,kTRUE);
4469 AliTracker::PropagateTrackToBxByBz(&(params[i]),xref,pt->GetMass(),1.,kTRUE);
4471 TMath::Sort(nseed,quality,indexes);
4473 // 3. Loop over pair of tracks
4475 for (Int_t i0=0; i0<nseed; i0++) {
4476 Int_t index0=indexes[i0];
4477 if (!(array->UncheckedAt(index0))) continue;
4478 AliTPCseed *s1 = (AliTPCseed*)array->UncheckedAt(index0);
4479 if (!s1->IsActive()) continue;
4480 AliExternalTrackParam &par0=params[index0];
4481 for (Int_t i1=i0+1; i1<nseed; i1++) {
4482 Int_t index1=indexes[i1];
4483 if (!(array->UncheckedAt(index1))) continue;
4484 AliTPCseed *s2 = (AliTPCseed*)array->UncheckedAt(index1);
4485 if (!s2->IsActive()) continue;
4486 if (s2->GetKinkIndexes()[0]!=0)
4487 if (s2->GetKinkIndexes()[0] == -s1->GetKinkIndexes()[0]) continue;
4488 AliExternalTrackParam &par1=params[index1];
4489 if (TMath::Abs(par0.GetParameter()[3]-par1.GetParameter()[3])>kCutP3) continue;
4490 if (TMath::Abs(par0.GetParameter()[1]-par1.GetParameter()[1])>kCutP1) continue;
4491 if (TMath::Abs(par0.GetParameter()[2]-par1.GetParameter()[2])>kCutP2) continue;
4492 Double_t dAlpha= TMath::Abs(par0.GetAlpha()-par1.GetAlpha());
4493 if (dAlpha>TMath::Pi()) dAlpha-=TMath::Pi();
4494 if (TMath::Abs(dAlpha)>kCutAlpha) continue;
4499 Int_t firstShared=lastpoint, lastShared=firstpoint;
4500 Int_t firstRow=lastpoint, lastRow=firstpoint;
4502 for (Int_t i=firstpoint;i<lastpoint;i++){
4503 if (s1->GetClusterIndex2(i)>0) nall0++;
4504 if (s2->GetClusterIndex2(i)>0) nall1++;
4505 if (s1->GetClusterIndex2(i)>0 && s2->GetClusterIndex2(i)>0) {
4506 if (i<firstRow) firstRow=i;
4507 if (i>lastRow) lastRow=i;
4509 if ( (s1->GetClusterIndex2(i))==(s2->GetClusterIndex2(i)) && s1->GetClusterIndex2(i)>0) {
4510 if (i<firstShared) firstShared=i;
4511 if (i>lastShared) lastShared=i;
4515 Double_t ratio0 = Float_t(sumShared)/Float_t(TMath::Min(nall0+1,nall1+1));
4516 Double_t ratio1 = Float_t(sumShared)/Float_t(TMath::Max(nall0+1,nall1+1));
4518 if( AliTPCReconstructor::StreamLevel()>1){
4519 TTreeSRedirector &cstream = *fDebugStreamer;
4520 Int_t n0=s1->GetNumberOfClusters();
4521 Int_t n1=s2->GetNumberOfClusters();
4522 Int_t n0F=s1->GetNFoundable();
4523 Int_t n1F=s2->GetNFoundable();
4524 Int_t lab0=s1->GetLabel();
4525 Int_t lab1=s2->GetLabel();
4527 cstream<<"Splitted2"<<
4528 "iter="<<fIteration<<
4529 "lab0="<<lab0<< // MC label if exist
4530 "lab1="<<lab1<< // MC label if exist
4533 "ratio0="<<ratio0<< // shared ratio
4534 "ratio1="<<ratio1<< // shared ratio
4535 "p0.="<<&par0<< // track parameters
4537 "s0.="<<s1<< // full seed
4539 "n0="<<n0<< // number of clusters track 0
4540 "n1="<<n1<< // number of clusters track 1
4541 "nall0="<<nall0<< // number of clusters track 0
4542 "nall1="<<nall1<< // number of clusters track 1
4543 "n0F="<<n0F<< // number of findable
4544 "n1F="<<n1F<< // number of findable
4545 "shared="<<sumShared<< // number of shared clusters
4546 "firstS="<<firstShared<< // first and the last shared row
4547 "lastS="<<lastShared<<
4548 "firstRow="<<firstRow<< // first and the last row with cluster
4549 "lastRow="<<lastRow<< //
4553 // remove track with lower quality
4555 if (ratio0>AliTPCReconstructor::GetRecoParam()->GetCutSharedClusters(0) ||
4556 ratio1>AliTPCReconstructor::GetRecoParam()->GetCutSharedClusters(1)){
4560 delete array->RemoveAt(index1);
4565 // 4. Delete temporary array
4575 void AliTPCtrackerMI::FindCurling(const TObjArray * array, AliESDEvent */*esd*/, Int_t iter)
4578 // find Curling tracks
4579 // Use AliTPCReconstructor::StreamLevel()>1 if you want to tune parameters - cuts
4582 // Algorithm done in 2 phases - because of CPU consumption
4583 // it is n^2 algorithm - for lead-lead 20000x20000 combination are investigated
4584 // see detal in MC part what can be used to cut
4588 const Float_t kMaxC = 400; // maximal curvature to of the track
4589 const Float_t kMaxdTheta = 0.15; // maximal distance in theta
4590 const Float_t kMaxdPhi = 0.15; // maximal distance in phi
4591 const Float_t kPtRatio = 0.3; // ratio between pt
4592 const Float_t kMinDCAR = 2.; // distance to the primary vertex in r - see cpipe cut
4595 // Curling tracks cuts
4598 const Float_t kMaxDeltaRMax = 40; // distance in outer radius
4599 const Float_t kMaxDeltaRMin = 5.; // distance in lower radius - see cpipe cut
4600 const Float_t kMinAngle = 2.9; // angle between tracks
4601 const Float_t kMaxDist = 5; // biggest distance
4603 // The cuts can be tuned using the "MC information stored in Multi tree ==> see FindMultiMC
4606 TCut csign("csign","Tr0.fP[4]*Tr1.fP[4]<0"); //opposite sign
4607 TCut cmax("cmax","abs(Tr0.GetC())>1/400");
4608 TCut cda("cda","sqrt(dtheta^2+dfi^2)<0.15");
4609 TCut ccratio("ccratio","abs((Tr0.fP[4]+Tr1.fP[4])/(abs(Tr0.fP[4])+abs(Tr1.fP[4])))<0.3");
4610 TCut cpipe("cpipe", "min(abs(r0-rc0),abs(r1-rc1))>5");
4612 TCut cdrmax("cdrmax","abs(abs(rc0+r0)-abs(rc1+r1))<40")
4613 TCut cdrmin("cdrmin","abs(abs(rc0+r0)-abs(rc1+r1))<10")
4615 Multi->Draw("dfi","iter==0"+csign+cmax+cda+ccratio); ~94% of curling tracks fulfill
4616 Multi->Draw("min(abs(r0-rc0),abs(r1-rc1))","iter==0&&abs(lab1)==abs(lab0)"+csign+cmax+cda+ccratio+cpipe+cdrmin+cdrmax); //80%
4618 Curling2->Draw("dfi","iter==0&&abs(lab0)==abs(lab1)"+csign+cmax+cdtheta+cdfi+ccratio)
4624 Int_t nentries = array->GetEntriesFast();
4625 AliHelix *helixes = new AliHelix[nentries];
4626 for (Int_t i=0;i<nentries;i++){
4627 AliTPCseed* track = (AliTPCseed*)array->At(i);
4628 if (!track) continue;
4629 track->SetCircular(0);
4630 new (&helixes[i]) AliHelix(*track);
4636 Double_t phase[2][2]={{0,0},{0,0}},radius[2]={0,0};
4642 for (Int_t i0=0;i0<nentries;i0++){
4643 AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
4644 if (!track0) continue;
4645 if (TMath::Abs(track0->GetC())<1/kMaxC) continue;
4646 Float_t xc0 = helixes[i0].GetHelix(6);
4647 Float_t yc0 = helixes[i0].GetHelix(7);
4648 Float_t r0 = helixes[i0].GetHelix(8);
4649 Float_t rc0 = TMath::Sqrt(xc0*xc0+yc0*yc0);
4650 Float_t fi0 = TMath::ATan2(yc0,xc0);
4652 for (Int_t i1=i0+1;i1<nentries;i1++){
4653 AliTPCseed * track1 = (AliTPCseed*)array->At(i1);
4654 if (!track1) continue;
4655 if (TMath::Abs(track1->GetC())<1/kMaxC) continue;
4656 Float_t xc1 = helixes[i1].GetHelix(6);
4657 Float_t yc1 = helixes[i1].GetHelix(7);
4658 Float_t r1 = helixes[i1].GetHelix(8);
4659 Float_t rc1 = TMath::Sqrt(xc1*xc1+yc1*yc1);
4660 Float_t fi1 = TMath::ATan2(yc1,xc1);
4662 Float_t dfi = fi0-fi1;
4665 if (dfi>1.5*TMath::Pi()) dfi-=TMath::Pi(); // take care about edge effect
4666 if (dfi<-1.5*TMath::Pi()) dfi+=TMath::Pi(); //
4667 Float_t dtheta = TMath::Abs(track0->GetTgl()-track1->GetTgl())<TMath::Abs(track0->GetTgl()+track1->GetTgl())? track0->GetTgl()-track1->GetTgl():track0->GetTgl()+track1->GetTgl();
4671 if (track0->GetBConstrain()&&track1->GetBConstrain()) continue; // not constrained
4672 if (track1->GetSigned1Pt()*track0->GetSigned1Pt()>0) continue; // not the same sign
4673 if ( TMath::Abs(track1->GetTgl()+track0->GetTgl())>kMaxdTheta) continue; //distance in the Theta
4674 if ( TMath::Abs(dfi)>kMaxdPhi) continue; //distance in phi
4675 if ( TMath::Sqrt(dfi*dfi+dtheta*dtheta)>kMaxdPhi) continue; //common angular offset
4677 Float_t pt0 = track0->GetSignedPt();
4678 Float_t pt1 = track1->GetSignedPt();
4679 if ((TMath::Abs(pt0+pt1)/(TMath::Abs(pt0)+TMath::Abs(pt1)))>kPtRatio) continue;
4680 if ((iter==1) && TMath::Abs(TMath::Abs(rc0+r0)-TMath::Abs(rc1+r1))>kMaxDeltaRMax) continue;
4681 if ((iter!=1) &&TMath::Abs(TMath::Abs(rc0-r0)-TMath::Abs(rc1-r1))>kMaxDeltaRMin) continue;
4682 if (TMath::Min(TMath::Abs(rc0-r0),TMath::Abs(rc1-r1))<kMinDCAR) continue;
4685 // Now find closest approach
4689 Int_t npoints = helixes[i0].GetRPHIintersections(helixes[i1], phase, radius,10);
4690 if (npoints==0) continue;
4691 helixes[i0].GetClosestPhases(helixes[i1], phase);
4695 Double_t hangles[3];
4696 helixes[i0].Evaluate(phase[0][0],xyz0);
4697 helixes[i1].Evaluate(phase[0][1],xyz1);
4699 helixes[i0].GetAngle(phase[0][0],helixes[i1],phase[0][1],hangles);
4700 Double_t deltah[2],deltabest;
4701 if (TMath::Abs(hangles[2])<kMinAngle) continue;
4705 helixes[i0].ParabolicDCA(helixes[i1],phase[0][0],phase[0][1],radius[0],deltah[0],2);
4707 helixes[i0].ParabolicDCA(helixes[i1],phase[1][0],phase[1][1],radius[1],deltah[1],2);
4708 if (deltah[1]<deltah[0]) ibest=1;
4710 deltabest = TMath::Sqrt(deltah[ibest]);
4711 helixes[i0].Evaluate(phase[ibest][0],xyz0);
4712 helixes[i1].Evaluate(phase[ibest][1],xyz1);
4713 helixes[i0].GetAngle(phase[ibest][0],helixes[i1],phase[ibest][1],hangles);
4714 Double_t radiusbest = TMath::Sqrt(radius[ibest]);
4716 if (deltabest>kMaxDist) continue;
4717 // if (mindcar+mindcaz<40 && (TMath::Abs(hangles[2])<kMinAngle ||deltabest>3)) continue;
4718 Bool_t sign =kFALSE;
4719 if (hangles[2]>kMinAngle) sign =kTRUE;
4722 // circular[i0] = kTRUE;
4723 // circular[i1] = kTRUE;
4724 if (track0->OneOverPt()<track1->OneOverPt()){
4725 track0->SetCircular(track0->GetCircular()+1);
4726 track1->SetCircular(track1->GetCircular()+2);
4729 track1->SetCircular(track1->GetCircular()+1);
4730 track0->SetCircular(track0->GetCircular()+2);
4733 if (AliTPCReconstructor::StreamLevel()>2){
4735 //debug stream to tune "fine" cuts
4736 Int_t lab0=track0->GetLabel();
4737 Int_t lab1=track1->GetLabel();
4738 TTreeSRedirector &cstream = *fDebugStreamer;
4739 cstream<<"Curling2"<<
4755 "npoints="<<npoints<<
4756 "hangles0="<<hangles[0]<<
4757 "hangles1="<<hangles[1]<<
4758 "hangles2="<<hangles[2]<<
4761 "radius="<<radiusbest<<
4762 "deltabest="<<deltabest<<
4763 "phase0="<<phase[ibest][0]<<
4764 "phase1="<<phase[ibest][1]<<
4772 if (AliTPCReconstructor::StreamLevel()>1) {
4773 AliInfo("Time for curling tracks removal");
4782 void AliTPCtrackerMI::FindKinks(TObjArray * array, AliESDEvent *esd)
4789 TObjArray *kinks= new TObjArray(10000);
4790 // TObjArray *v0s= new TObjArray(10000);
4791 Int_t nentries = array->GetEntriesFast();
4792 AliHelix *helixes = new AliHelix[nentries];
4793 Int_t *sign = new Int_t[nentries];
4794 Int_t *nclusters = new Int_t[nentries];
4795 Float_t *alpha = new Float_t[nentries];
4796 AliKink *kink = new AliKink();
4797 Int_t * usage = new Int_t[nentries];
4798 Float_t *zm = new Float_t[nentries];
4799 Float_t *z0 = new Float_t[nentries];
4800 Float_t *fim = new Float_t[nentries];
4801 Float_t *shared = new Float_t[nentries];
4802 Bool_t *circular = new Bool_t[nentries];
4803 Float_t *dca = new Float_t[nentries];
4804 //const AliESDVertex * primvertex = esd->GetVertex();
4806 // nentries = array->GetEntriesFast();
4811 for (Int_t i=0;i<nentries;i++){
4814 AliTPCseed* track = (AliTPCseed*)array->At(i);
4815 if (!track) continue;
4816 track->SetCircular(0);
4818 track->UpdatePoints();
4819 if (( track->GetPoints()[2]- track->GetPoints()[0])>5 && track->GetPoints()[3]>0.8){
4821 nclusters[i]=track->GetNumberOfClusters();
4822 alpha[i] = track->GetAlpha();
4823 new (&helixes[i]) AliHelix(*track);
4825 helixes[i].Evaluate(0,xyz);
4826 sign[i] = (track->GetC()>0) ? -1:1;
4829 if (track->GetProlongation(x,y,z)){
4831 fim[i] = alpha[i]+TMath::ATan2(y,x);
4834 zm[i] = track->GetZ();
4838 circular[i]= kFALSE;
4839 if (track->GetProlongation(0,y,z)) z0[i] = z;
4840 dca[i] = track->GetD(0,0);
4846 Int_t ncandidates =0;
4849 Double_t phase[2][2]={{0,0},{0,0}},radius[2]={0,0};
4852 // Find circling track
4854 for (Int_t i0=0;i0<nentries;i0++){
4855 AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
4856 if (!track0) continue;
4857 if (track0->GetNumberOfClusters()<40) continue;
4858 if (TMath::Abs(1./track0->GetC())>200) continue;
4859 for (Int_t i1=i0+1;i1<nentries;i1++){
4860 AliTPCseed * track1 = (AliTPCseed*)array->At(i1);
4861 if (!track1) continue;
4862 if (track1->GetNumberOfClusters()<40) continue;
4863 if ( TMath::Abs(track1->GetTgl()+track0->GetTgl())>0.1) continue;
4864 if (track0->GetBConstrain()&&track1->GetBConstrain()) continue;
4865 if (TMath::Abs(1./track1->GetC())>200) continue;
4866 if (track1->GetSigned1Pt()*track0->GetSigned1Pt()>0) continue;
4867 if (track1->GetTgl()*track0->GetTgl()>0) continue;
4868 if (TMath::Max(TMath::Abs(1./track0->GetC()),TMath::Abs(1./track1->GetC()))>190) continue;
4869 if (track0->GetBConstrain()&&track1->OneOverPt()<track0->OneOverPt()) continue; //returning - lower momenta
4870 if (track1->GetBConstrain()&&track0->OneOverPt()<track1->OneOverPt()) continue; //returning - lower momenta
4872 Float_t mindcar = TMath::Min(TMath::Abs(dca[i0]),TMath::Abs(dca[i1]));
4873 if (mindcar<5) continue;
4874 Float_t mindcaz = TMath::Min(TMath::Abs(z0[i0]-GetZ()),TMath::Abs(z0[i1]-GetZ()));
4875 if (mindcaz<5) continue;
4876 if (mindcar+mindcaz<20) continue;
4879 Float_t xc0 = helixes[i0].GetHelix(6);
4880 Float_t yc0 = helixes[i0].GetHelix(7);
4881 Float_t r0 = helixes[i0].GetHelix(8);
4882 Float_t xc1 = helixes[i1].GetHelix(6);
4883 Float_t yc1 = helixes[i1].GetHelix(7);
4884 Float_t r1 = helixes[i1].GetHelix(8);
4886 Float_t rmean = (r0+r1)*0.5;
4887 Float_t delta =TMath::Sqrt((xc1-xc0)*(xc1-xc0)+(yc1-yc0)*(yc1-yc0));
4888 //if (delta>30) continue;
4889 if (delta>rmean*0.25) continue;
4890 if (TMath::Abs(r0-r1)/rmean>0.3) continue;
4892 Int_t npoints = helixes[i0].GetRPHIintersections(helixes[i1], phase, radius,10);
4893 if (npoints==0) continue;
4894 helixes[i0].GetClosestPhases(helixes[i1], phase);
4898 Double_t hangles[3];
4899 helixes[i0].Evaluate(phase[0][0],xyz0);
4900 helixes[i1].Evaluate(phase[0][1],xyz1);
4902 helixes[i0].GetAngle(phase[0][0],helixes[i1],phase[0][1],hangles);
4903 Double_t deltah[2],deltabest;
4904 if (hangles[2]<2.8) continue;
4907 helixes[i0].ParabolicDCA(helixes[i1],phase[0][0],phase[0][1],radius[0],deltah[0],2);
4909 helixes[i0].ParabolicDCA(helixes[i1],phase[1][0],phase[1][1],radius[1],deltah[1],2);
4910 if (deltah[1]<deltah[0]) ibest=1;
4912 deltabest = TMath::Sqrt(deltah[ibest]);
4913 helixes[i0].Evaluate(phase[ibest][0],xyz0);
4914 helixes[i1].Evaluate(phase[ibest][1],xyz1);
4915 helixes[i0].GetAngle(phase[ibest][0],helixes[i1],phase[ibest][1],hangles);
4916 Double_t radiusbest = TMath::Sqrt(radius[ibest]);
4918 if (deltabest>6) continue;
4919 if (mindcar+mindcaz<40 && (hangles[2]<3.12||deltabest>3)) continue;
4920 Bool_t lsign =kFALSE;
4921 if (hangles[2]>3.06) lsign =kTRUE;
4924 circular[i0] = kTRUE;
4925 circular[i1] = kTRUE;
4926 if (track0->OneOverPt()<track1->OneOverPt()){
4927 track0->SetCircular(track0->GetCircular()+1);
4928 track1->SetCircular(track1->GetCircular()+2);
4931 track1->SetCircular(track1->GetCircular()+1);
4932 track0->SetCircular(track0->GetCircular()+2);
4935 if (lsign&&AliTPCReconstructor::StreamLevel()>1){
4937 Int_t lab0=track0->GetLabel();
4938 Int_t lab1=track1->GetLabel();
4939 TTreeSRedirector &cstream = *fDebugStreamer;
4940 cstream<<"Curling"<<
4947 "mindcar="<<mindcar<<
4948 "mindcaz="<<mindcaz<<
4951 "npoints="<<npoints<<
4952 "hangles0="<<hangles[0]<<
4953 "hangles2="<<hangles[2]<<
4958 "radius="<<radiusbest<<
4959 "deltabest="<<deltabest<<
4960 "phase0="<<phase[ibest][0]<<
4961 "phase1="<<phase[ibest][1]<<
4971 for (Int_t i =0;i<nentries;i++){
4972 if (sign[i]==0) continue;
4973 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
4980 Double_t cradius0 = 40*40;
4981 Double_t cradius1 = 270*270;
4984 Double_t cdist3=0.55;
4985 for (Int_t j =i+1;j<nentries;j++){
4987 if (sign[j]*sign[i]<1) continue;
4988 if ( (nclusters[i]+nclusters[j])>200) continue;
4989 if ( (nclusters[i]+nclusters[j])<80) continue;
4990 if ( TMath::Abs(zm[i]-zm[j])>60.) continue;
4991 if ( TMath::Abs(fim[i]-fim[j])>0.6 && TMath::Abs(fim[i]-fim[j])<5.7 ) continue;
4992 //AliTPCseed * track1 = (AliTPCseed*)array->At(j); Double_t phase[2][2],radius[2];
4993 Int_t npoints = helixes[i].GetRPHIintersections(helixes[j], phase, radius,20);
4994 if (npoints<1) continue;
4997 if (radius[0]<cradius0||radius[0]>cradius1) continue;
5000 if ( (radius[0]<cradius0||radius[0]>cradius1) && (radius[1]<cradius0||radius[1]>cradius1) ) continue;
5003 Double_t delta1=10000,delta2=10000;
5004 // cuts on the intersection radius
5005 helixes[i].LinearDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta1);
5006 if (radius[0]<20&&delta1<1) continue; //intersection at vertex
5007 if (radius[0]<10&&delta1<3) continue; //intersection at vertex
5009 helixes[i].LinearDCA(helixes[j],phase[1][0],phase[1][1],radius[1],delta2);
5010 if (radius[1]<20&&delta2<1) continue; //intersection at vertex
5011 if (radius[1]<10&&delta2<3) continue; //intersection at vertex
5014 Double_t distance1 = TMath::Min(delta1,delta2);
5015 if (distance1>cdist1) continue; // cut on DCA linear approximation
5017 npoints = helixes[i].GetRPHIintersections(helixes[j], phase, radius,20);
5018 helixes[i].ParabolicDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta1);
5019 if (radius[0]<20&&delta1<1) continue; //intersection at vertex
5020 if (radius[0]<10&&delta1<3) continue; //intersection at vertex
5023 helixes[i].ParabolicDCA(helixes[j],phase[1][0],phase[1][1],radius[1],delta2);
5024 if (radius[1]<20&&delta2<1) continue; //intersection at vertex
5025 if (radius[1]<10&&delta2<3) continue; //intersection at vertex
5027 distance1 = TMath::Min(delta1,delta2);
5030 rkink = TMath::Sqrt(radius[0]);
5033 rkink = TMath::Sqrt(radius[1]);
5035 if (distance1>cdist2) continue;
5038 AliTPCseed * track1 = (AliTPCseed*)array->At(j);
5041 Int_t row0 = GetRowNumber(rkink);
5042 if (row0<10) continue;
5043 if (row0>150) continue;
5046 Float_t dens00=-1,dens01=-1;
5047 Float_t dens10=-1,dens11=-1;
5049 Int_t found,foundable,ishared;
5050 track0->GetClusterStatistic(0,row0-5, found, foundable,ishared,kFALSE);
5051 if (foundable>5) dens00 = Float_t(found)/Float_t(foundable);
5052 track0->GetClusterStatistic(row0+5,155, found, foundable,ishared,kFALSE);
5053 if (foundable>5) dens01 = Float_t(found)/Float_t(foundable);
5055 track1->GetClusterStatistic(0,row0-5, found, foundable,ishared,kFALSE);
5056 if (foundable>10) dens10 = Float_t(found)/Float_t(foundable);
5057 track1->GetClusterStatistic(row0+5,155, found, foundable,ishared,kFALSE);
5058 if (foundable>10) dens11 = Float_t(found)/Float_t(foundable);
5060 if (dens00<dens10 && dens01<dens11) continue;
5061 if (dens00>dens10 && dens01>dens11) continue;
5062 if (TMath::Max(dens00,dens10)<0.1) continue;
5063 if (TMath::Max(dens01,dens11)<0.3) continue;
5065 if (TMath::Min(dens00,dens10)>0.6) continue;
5066 if (TMath::Min(dens01,dens11)>0.6) continue;
5069 AliTPCseed * ktrack0, *ktrack1;
5078 if (TMath::Abs(ktrack0->GetC())>5) continue; // cut on the curvature for mother particle
5079 AliExternalTrackParam paramm(*ktrack0);
5080 AliExternalTrackParam paramd(*ktrack1);
5081 if (row0>60&&ktrack1->GetReference().GetX()>90.)new (¶md) AliExternalTrackParam(ktrack1->GetReference());
5084 kink->SetMother(paramm);
5085 kink->SetDaughter(paramd);
5088 Float_t x[3] = { kink->GetPosition()[0],kink->GetPosition()[1],kink->GetPosition()[2]};
5090 fkParam->Transform0to1(x,index);
5091 fkParam->Transform1to2(x,index);
5092 row0 = GetRowNumber(x[0]);
5094 if (kink->GetR()<100) continue;
5095 if (kink->GetR()>240) continue;
5096 if (kink->GetPosition()[2]/kink->GetR()>AliTPCReconstructor::GetCtgRange()) continue; //out of fiducial volume
5097 if (kink->GetDistance()>cdist3) continue;
5098 Float_t dird = kink->GetDaughterP()[0]*kink->GetPosition()[0]+kink->GetDaughterP()[1]*kink->GetPosition()[1]; // rough direction estimate
5099 if (dird<0) continue;
5101 Float_t dirm = kink->GetMotherP()[0]*kink->GetPosition()[0]+kink->GetMotherP()[1]*kink->GetPosition()[1]; // rough direction estimate
5102 if (dirm<0) continue;
5103 Float_t mpt = TMath::Sqrt(kink->GetMotherP()[0]*kink->GetMotherP()[0]+kink->GetMotherP()[1]*kink->GetMotherP()[1]);
5104 if (mpt<0.2) continue;
5107 //for high momenta momentum not defined well in first iteration
5108 Double_t qt = TMath::Sin(kink->GetAngle(2))*ktrack1->GetP();
5109 if (qt>0.35) continue;
5112 kink->SetLabel(CookLabel(ktrack0,0.4,0,row0),0);
5113 kink->SetLabel(CookLabel(ktrack1,0.4,row0,160),1);
5115 kink->SetTPCDensity(dens00,0,0);
5116 kink->SetTPCDensity(dens01,0,1);
5117 kink->SetTPCDensity(dens10,1,0);
5118 kink->SetTPCDensity(dens11,1,1);
5119 kink->SetIndex(i,0);
5120 kink->SetIndex(j,1);
5123 kink->SetTPCDensity(dens10,0,0);
5124 kink->SetTPCDensity(dens11,0,1);
5125 kink->SetTPCDensity(dens00,1,0);
5126 kink->SetTPCDensity(dens01,1,1);
5127 kink->SetIndex(j,0);
5128 kink->SetIndex(i,1);
5131 if (mpt<1||kink->GetAngle(2)>0.1){
5132 // angle and densities not defined yet
5133 if (kink->GetTPCDensityFactor()<0.8) continue;
5134 if ((2-kink->GetTPCDensityFactor())*kink->GetDistance() >0.25) continue;
5135 if (kink->GetAngle(2)*ktrack0->GetP()<0.003) continue; //too small angle
5136 if (kink->GetAngle(2)>0.2&&kink->GetTPCDensityFactor()<1.15) continue;
5137 if (kink->GetAngle(2)>0.2&&kink->GetTPCDensity(0,1)>0.05) continue;
5139 Float_t criticalangle = track0->GetSigmaSnp2()+track0->GetSigmaTgl2();
5140 criticalangle+= track1->GetSigmaSnp2()+track1->GetSigmaTgl2();
5141 criticalangle= 3*TMath::Sqrt(criticalangle);
5142 if (criticalangle>0.02) criticalangle=0.02;
5143 if (kink->GetAngle(2)<criticalangle) continue;
5146 Int_t drow = Int_t(2.+0.5/(0.05+kink->GetAngle(2))); // overlap region defined
5147 Float_t shapesum =0;
5149 for ( Int_t row = row0-drow; row<row0+drow;row++){
5150 if (row<0) continue;
5151 if (row>155) continue;
5152 if (ktrack0->GetClusterPointer(row)){
5153 AliTPCTrackerPoint *point =ktrack0->GetTrackPoint(row);
5154 shapesum+=point->GetSigmaY()+point->GetSigmaZ();
5157 if (ktrack1->GetClusterPointer(row)){
5158 AliTPCTrackerPoint *point =ktrack1->GetTrackPoint(row);
5159 shapesum+=point->GetSigmaY()+point->GetSigmaZ();
5164 kink->SetShapeFactor(-1.);
5167 kink->SetShapeFactor(shapesum/sum);
5169 // esd->AddKink(kink);
5171 // kink->SetMother(paramm);
5172 //kink->SetDaughter(paramd);
5174 Double_t chi2P2 = paramm.GetParameter()[2]-paramd.GetParameter()[2];
5176 chi2P2/=paramm.GetCovariance()[5]+paramd.GetCovariance()[5];
5177 Double_t chi2P3 = paramm.GetParameter()[3]-paramd.GetParameter()[3];
5179 chi2P3/=paramm.GetCovariance()[9]+paramd.GetCovariance()[9];
5181 if (AliTPCReconstructor::StreamLevel()>1) {
5182 (*fDebugStreamer)<<"kinkLpt"<<
5190 if ( chi2P2+chi2P3<AliTPCReconstructor::GetRecoParam()->GetKinkAngleCutChi2(0)){
5194 kinks->AddLast(kink);
5200 // sort the kinks according quality - and refit them towards vertex
5202 Int_t nkinks = kinks->GetEntriesFast();
5203 Float_t *quality = new Float_t[nkinks];
5204 Int_t *indexes = new Int_t[nkinks];
5205 AliTPCseed *mothers = new AliTPCseed[nkinks];
5206 AliTPCseed *daughters = new AliTPCseed[nkinks];
5209 for (Int_t i=0;i<nkinks;i++){
5211 AliKink *kinkl = (AliKink*)kinks->At(i);
5213 // refit kinks towards vertex
5215 Int_t index0 = kinkl->GetIndex(0);
5216 Int_t index1 = kinkl->GetIndex(1);
5217 AliTPCseed * ktrack0 = (AliTPCseed*)array->At(index0);
5218 AliTPCseed * ktrack1 = (AliTPCseed*)array->At(index1);
5220 Int_t sumn=ktrack0->GetNumberOfClusters()+ktrack1->GetNumberOfClusters();
5222 // Refit Kink under if too small angle
5224 if (kinkl->GetAngle(2)<0.05){
5225 kinkl->SetTPCRow0(GetRowNumber(kinkl->GetR()));
5226 Int_t row0 = kinkl->GetTPCRow0();
5227 Int_t drow = Int_t(2.+0.5/(0.05+kinkl->GetAngle(2)));
5230 Int_t last = row0-drow;
5231 if (last<40) last=40;
5232 if (last<ktrack0->GetFirstPoint()+25) last = ktrack0->GetFirstPoint()+25;
5233 AliTPCseed* seed0 = ReSeed(ktrack0,last,kFALSE);
5236 Int_t first = row0+drow;
5237 if (first>130) first=130;
5238 if (first>ktrack1->GetLastPoint()-25) first = TMath::Max(ktrack1->GetLastPoint()-25,30);
5239 AliTPCseed* seed1 = ReSeed(ktrack1,first,kTRUE);
5241 if (seed0 && seed1){
5242 kinkl->SetStatus(1,8);
5243 if (RefitKink(*seed0,*seed1,*kinkl)) kinkl->SetStatus(1,9);
5244 row0 = GetRowNumber(kinkl->GetR());
5245 sumn = seed0->GetNumberOfClusters()+seed1->GetNumberOfClusters();
5246 mothers[i] = *seed0;
5247 daughters[i] = *seed1;
5250 delete kinks->RemoveAt(i);
5251 if (seed0) delete seed0;
5252 if (seed1) delete seed1;
5255 if (kinkl->GetDistance()>0.5 || kinkl->GetR()<110 || kinkl->GetR()>240) {
5256 delete kinks->RemoveAt(i);
5257 if (seed0) delete seed0;
5258 if (seed1) delete seed1;
5266 if (kinkl) quality[i] = 160*((0.1+kinkl->GetDistance())*(2.-kinkl->GetTPCDensityFactor()))/(sumn+40.); //the longest -clossest will win
5268 TMath::Sort(nkinks,quality,indexes,kFALSE);
5270 //remove double find kinks
5272 for (Int_t ikink0=1;ikink0<nkinks;ikink0++){
5273 AliKink * kink0 = (AliKink*) kinks->At(indexes[ikink0]);
5274 if (!kink0) continue;
5276 for (Int_t ikink1=0;ikink1<ikink0;ikink1++){
5277 kink0 = (AliKink*) kinks->At(indexes[ikink0]);
5278 if (!kink0) continue;
5279 AliKink * kink1 = (AliKink*) kinks->At(indexes[ikink1]);
5280 if (!kink1) continue;
5281 // if not close kink continue
5282 if (TMath::Abs(kink1->GetPosition()[2]-kink0->GetPosition()[2])>10) continue;
5283 if (TMath::Abs(kink1->GetPosition()[1]-kink0->GetPosition()[1])>10) continue;
5284 if (TMath::Abs(kink1->GetPosition()[0]-kink0->GetPosition()[0])>10) continue;
5286 AliTPCseed &mother0 = mothers[indexes[ikink0]];
5287 AliTPCseed &daughter0 = daughters[indexes[ikink0]];
5288 AliTPCseed &mother1 = mothers[indexes[ikink1]];
5289 AliTPCseed &daughter1 = daughters[indexes[ikink1]];
5290 Int_t row0 = (kink0->GetTPCRow0()+kink1->GetTPCRow0())/2;
5299 for (Int_t i=0;i<row0;i++){
5300 if (mother0.GetClusterIndex(i)>0 && mother1.GetClusterIndex(i)>0){
5303 if (mother0.GetClusterIndex(i)==mother1.GetClusterIndex(i)){
5310 for (Int_t i=row0;i<158;i++){
5311 if (daughter0.GetClusterIndex(i)>0 && daughter0.GetClusterIndex(i)>0){
5314 if (mother0.GetClusterIndex(i)==mother1.GetClusterIndex(i)){
5320 Float_t ratio = Float_t(same+1)/Float_t(both+1);
5321 Float_t ratiom = Float_t(samem+1)/Float_t(bothm+1);
5322 Float_t ratiod = Float_t(samed+1)/Float_t(bothd+1);
5323 if (ratio>0.3 && ratiom>0.5 &&ratiod>0.5) {
5324 Int_t sum0 = mother0.GetNumberOfClusters()+daughter0.GetNumberOfClusters();
5325 Int_t sum1 = mother1.GetNumberOfClusters()+daughter1.GetNumberOfClusters();
5327 shared[kink0->GetIndex(0)]= kTRUE;
5328 shared[kink0->GetIndex(1)]= kTRUE;
5329 delete kinks->RemoveAt(indexes[ikink0]);
5333 shared[kink1->GetIndex(0)]= kTRUE;
5334 shared[kink1->GetIndex(1)]= kTRUE;
5335 delete kinks->RemoveAt(indexes[ikink1]);
5342 for (Int_t i=0;i<nkinks;i++){
5343 AliKink * kinkl = (AliKink*) kinks->At(indexes[i]);
5344 if (!kinkl) continue;
5345 kinkl->SetTPCRow0(GetRowNumber(kinkl->GetR()));
5346 Int_t index0 = kinkl->GetIndex(0);
5347 Int_t index1 = kinkl->GetIndex(1);
5348 if (circular[index0]||(circular[index1]&&kinkl->GetDistance()>0.2)) continue;
5349 kinkl->SetMultiple(usage[index0],0);
5350 kinkl->SetMultiple(usage[index1],1);
5351 if (kinkl->GetMultiple()[0]+kinkl->GetMultiple()[1]>2) continue;
5352 if (kinkl->GetMultiple()[0]+kinkl->GetMultiple()[1]>0 && quality[indexes[i]]>0.2) continue;
5353 if (kinkl->GetMultiple()[0]+kinkl->GetMultiple()[1]>0 && kinkl->GetDistance()>0.2) continue;
5354 if (circular[index0]||(circular[index1]&&kinkl->GetDistance()>0.1)) continue;
5356 AliTPCseed * ktrack0 = (AliTPCseed*)array->At(index0);
5357 AliTPCseed * ktrack1 = (AliTPCseed*)array->At(index1);
5358 if (!ktrack0 || !ktrack1) continue;
5359 Int_t index = esd->AddKink(kinkl);
5362 if ( ktrack0->GetKinkIndex(0)==0 && ktrack1->GetKinkIndex(0)==0) { //best kink
5363 if (mothers[indexes[i]].GetNumberOfClusters()>20 && daughters[indexes[i]].GetNumberOfClusters()>20 && (mothers[indexes[i]].GetNumberOfClusters()+daughters[indexes[i]].GetNumberOfClusters())>100){
5364 *ktrack0 = mothers[indexes[i]];
5365 *ktrack1 = daughters[indexes[i]];
5369 ktrack0->SetKinkIndex(usage[index0],-(index+1));
5370 ktrack1->SetKinkIndex(usage[index1], (index+1));
5375 // Remove tracks corresponding to shared kink's
5377 for (Int_t i=0;i<nentries;i++){
5378 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5379 if (!track0) continue;
5380 if (track0->GetKinkIndex(0)!=0) continue;
5381 if (shared[i]) delete array->RemoveAt(i);
5386 RemoveUsed2(array,0.5,0.4,30);
5388 for (Int_t i=0;i<nentries;i++){
5389 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5390 if (!track0) continue;
5391 track0->CookdEdx(0.02,0.6);
5395 for (Int_t i=0;i<nentries;i++){
5396 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5397 if (!track0) continue;
5398 if (track0->Pt()<1.4) continue;
5399 //remove double high momenta tracks - overlapped with kink candidates
5402 for (Int_t icl=track0->GetFirstPoint();icl<track0->GetLastPoint(); icl++){
5403 if (track0->GetClusterPointer(icl)!=0){
5405 if (track0->GetClusterPointer(icl)->IsUsed(10)) ishared++;
5408 if (Float_t(ishared+1)/Float_t(all+1)>0.5) {
5409 delete array->RemoveAt(i);
5413 if (track0->GetKinkIndex(0)!=0) continue;
5414 if (track0->GetNumberOfClusters()<80) continue;
5416 AliTPCseed *pmother = new AliTPCseed();
5417 AliTPCseed *pdaughter = new AliTPCseed();
5418 AliKink *pkink = new AliKink;
5420 AliTPCseed & mother = *pmother;
5421 AliTPCseed & daughter = *pdaughter;
5422 AliKink & kinkl = *pkink;
5423 if (CheckKinkPoint(track0,mother,daughter, kinkl)){
5424 if (mother.GetNumberOfClusters()<30||daughter.GetNumberOfClusters()<20) {
5428 continue; //too short tracks
5430 if (mother.Pt()<1.4) {
5436 Int_t row0= kinkl.GetTPCRow0();
5437 if (kinkl.GetDistance()>0.5 || kinkl.GetR()<110. || kinkl.GetR()>240.) {
5444 Int_t index = esd->AddKink(&kinkl);
5445 mother.SetKinkIndex(0,-(index+1));
5446 daughter.SetKinkIndex(0,index+1);
5447 if (mother.GetNumberOfClusters()>50) {
5448 delete array->RemoveAt(i);
5449 array->AddAt(new AliTPCseed(mother),i);
5452 array->AddLast(new AliTPCseed(mother));
5454 array->AddLast(new AliTPCseed(daughter));
5455 for (Int_t icl=0;icl<row0;icl++) {
5456 if (mother.GetClusterPointer(icl)) mother.GetClusterPointer(icl)->Use(20);
5459 for (Int_t icl=row0;icl<158;icl++) {
5460 if (daughter.GetClusterPointer(icl)) daughter.GetClusterPointer(icl)->Use(20);
5469 delete [] daughters;
5491 AliInfo(Form("Ncandidates=\t%d\t%d\t%d\t%d\n",esd->GetNumberOfKinks(),ncandidates,ntracks,nall));
5496 Int_t AliTPCtrackerMI::RefitKink(AliTPCseed &mother, AliTPCseed &daughter, const AliESDkink &knk)
5499 // refit kink towards to the vertex
5502 AliKink &kink=(AliKink &)knk;
5504 Int_t row0 = GetRowNumber(kink.GetR());
5505 FollowProlongation(mother,0);
5506 mother.Reset(kFALSE);
5508 FollowProlongation(daughter,row0);
5509 daughter.Reset(kFALSE);
5510 FollowBackProlongation(daughter,158);
5511 daughter.Reset(kFALSE);
5512 Int_t first = TMath::Max(row0-20,30);
5513 Int_t last = TMath::Min(row0+20,140);
5515 const Int_t kNdiv =5;
5516 AliTPCseed param0[kNdiv]; // parameters along the track
5517 AliTPCseed param1[kNdiv]; // parameters along the track
5518 AliKink kinks[kNdiv]; // corresponding kink parameters
5521 for (Int_t irow=0; irow<kNdiv;irow++){
5522 rows[irow] = first +((last-first)*irow)/(kNdiv-1);
5524 // store parameters along the track
5526 for (Int_t irow=0;irow<kNdiv;irow++){
5527 FollowBackProlongation(mother, rows[irow]);
5528 FollowProlongation(daughter,rows[kNdiv-1-irow]);
5529 param0[irow] = mother;
5530 param1[kNdiv-1-irow] = daughter;
5534 for (Int_t irow=0; irow<kNdiv-1;irow++){
5535 if (param0[irow].GetNumberOfClusters()<kNdiv||param1[irow].GetNumberOfClusters()<kNdiv) continue;
5536 kinks[irow].SetMother(param0[irow]);
5537 kinks[irow].SetDaughter(param1[irow]);
5538 kinks[irow].Update();
5541 // choose kink with best "quality"
5543 Double_t mindist = 10000;
5544 for (Int_t irow=0;irow<kNdiv;irow++){
5545 if (param0[irow].GetNumberOfClusters()<20||param1[irow].GetNumberOfClusters()<20) continue;
5546 if (TMath::Abs(kinks[irow].GetR())>240.) continue;
5547 if (TMath::Abs(kinks[irow].GetR())<100.) continue;
5549 Float_t normdist = TMath::Abs(param0[irow].GetX()-kinks[irow].GetR())*(0.1+kink.GetDistance());
5550 normdist/= (param0[irow].GetNumberOfClusters()+param1[irow].GetNumberOfClusters()+40.);
5551 if (normdist < mindist){
5557 if (index==-1) return 0;
5560 param0[index].Reset(kTRUE);
5561 FollowProlongation(param0[index],0);
5563 mother = param0[index];
5564 daughter = param1[index]; // daughter in vertex
5566 kink.SetMother(mother);
5567 kink.SetDaughter(daughter);
5569 kink.SetTPCRow0(GetRowNumber(kink.GetR()));
5570 kink.SetTPCncls(param0[index].GetNumberOfClusters(),0);
5571 kink.SetTPCncls(param1[index].GetNumberOfClusters(),1);
5572 kink.SetLabel(CookLabel(&mother,0.4, 0,kink.GetTPCRow0()),0);
5573 kink.SetLabel(CookLabel(&daughter,0.4, kink.GetTPCRow0(),160),1);
5574 mother.SetLabel(kink.GetLabel(0));
5575 daughter.SetLabel(kink.GetLabel(1));
5581 void AliTPCtrackerMI::UpdateKinkQualityM(AliTPCseed * seed){
5583 // update Kink quality information for mother after back propagation
5585 if (seed->GetKinkIndex(0)>=0) return;
5586 for (Int_t ikink=0;ikink<3;ikink++){
5587 Int_t index = seed->GetKinkIndex(ikink);
5588 if (index>=0) break;
5589 index = TMath::Abs(index)-1;
5590 AliESDkink * kink = fEvent->GetKink(index);
5591 kink->SetTPCDensity(-1,0,0);
5592 kink->SetTPCDensity(1,0,1);
5594 Int_t row0 = kink->GetTPCRow0() - 2 - Int_t( 0.5/ (0.05+kink->GetAngle(2)));
5595 if (row0<15) row0=15;
5597 Int_t row1 = kink->GetTPCRow0() + 2 + Int_t( 0.5/ (0.05+kink->GetAngle(2)));
5598 if (row1>145) row1=145;
5600 Int_t found,foundable,shared;
5601 seed->GetClusterStatistic(0,row0, found, foundable,shared,kFALSE);
5602 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),0,0);
5603 seed->GetClusterStatistic(row1,155, found, foundable,shared,kFALSE);
5604 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),0,1);
5609 void AliTPCtrackerMI::UpdateKinkQualityD(AliTPCseed * seed){
5611 // update Kink quality information for daughter after refit
5613 if (seed->GetKinkIndex(0)<=0) return;
5614 for (Int_t ikink=0;ikink<3;ikink++){
5615 Int_t index = seed->GetKinkIndex(ikink);
5616 if (index<=0) break;
5617 index = TMath::Abs(index)-1;
5618 AliESDkink * kink = fEvent->GetKink(index);
5619 kink->SetTPCDensity(-1,1,0);
5620 kink->SetTPCDensity(-1,1,1);
5622 Int_t row0 = kink->GetTPCRow0() -2 - Int_t( 0.5/ (0.05+kink->GetAngle(2)));
5623 if (row0<15) row0=15;
5625 Int_t row1 = kink->GetTPCRow0() +2 + Int_t( 0.5/ (0.05+kink->GetAngle(2)));
5626 if (row1>145) row1=145;
5628 Int_t found,foundable,shared;
5629 seed->GetClusterStatistic(0,row0, found, foundable,shared,kFALSE);
5630 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),1,0);
5631 seed->GetClusterStatistic(row1,155, found, foundable,shared,kFALSE);
5632 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),1,1);
5638 Int_t AliTPCtrackerMI::CheckKinkPoint(AliTPCseed*seed,AliTPCseed &mother, AliTPCseed &daughter, const AliESDkink &knk)
5641 // check kink point for given track
5642 // if return value=0 kink point not found
5643 // otherwise seed0 correspond to mother particle
5644 // seed1 correspond to daughter particle
5645 // kink parameter of kink point
5646 AliKink &kink=(AliKink &)knk;
5648 Int_t middlerow = (seed->GetFirstPoint()+seed->GetLastPoint())/2;
5649 Int_t first = seed->GetFirstPoint();
5650 Int_t last = seed->GetLastPoint();
5651 if (last-first<20) return 0; // shortest length - 2*30 = 60 pad-rows
5654 AliTPCseed *seed1 = ReSeed(seed,middlerow+20, kTRUE); //middle of chamber
5655 if (!seed1) return 0;
5656 FollowProlongation(*seed1,seed->GetLastPoint()-20);
5657 seed1->Reset(kTRUE);
5658 FollowProlongation(*seed1,158);
5659 seed1->Reset(kTRUE);
5660 last = seed1->GetLastPoint();
5662 AliTPCseed *seed0 = new AliTPCseed(*seed);
5663 seed0->Reset(kFALSE);
5666 AliTPCseed param0[20]; // parameters along the track
5667 AliTPCseed param1[20]; // parameters along the track
5668 AliKink kinks[20]; // corresponding kink parameters
5670 for (Int_t irow=0; irow<20;irow++){
5671 rows[irow] = first +((last-first)*irow)/19;
5673 // store parameters along the track
5675 for (Int_t irow=0;irow<20;irow++){
5676 FollowBackProlongation(*seed0, rows[irow]);
5677 FollowProlongation(*seed1,rows[19-irow]);
5678 param0[irow] = *seed0;
5679 param1[19-irow] = *seed1;
5683 for (Int_t irow=0; irow<19;irow++){
5684 kinks[irow].SetMother(param0[irow]);
5685 kinks[irow].SetDaughter(param1[irow]);
5686 kinks[irow].Update();
5689 // choose kink with biggest change of angle
5691 Double_t maxchange= 0;
5692 for (Int_t irow=1;irow<19;irow++){
5693 if (TMath::Abs(kinks[irow].GetR())>240.) continue;
5694 if (TMath::Abs(kinks[irow].GetR())<110.) continue;
5695 Float_t quality = TMath::Abs(kinks[irow].GetAngle(2))/(3.+TMath::Abs(kinks[irow].GetR()-param0[irow].GetX()));
5696 if ( quality > maxchange){
5697 maxchange = quality;
5704 if (index<0) return 0;
5706 Int_t row0 = GetRowNumber(kinks[index].GetR()); //row 0 estimate
5707 seed0 = new AliTPCseed(param0[index]);
5708 seed1 = new AliTPCseed(param1[index]);
5709 seed0->Reset(kFALSE);
5710 seed1->Reset(kFALSE);
5711 seed0->ResetCovariance(10.);
5712 seed1->ResetCovariance(10.);
5713 FollowProlongation(*seed0,0);
5714 FollowBackProlongation(*seed1,158);
5715 mother = *seed0; // backup mother at position 0
5716 seed0->Reset(kFALSE);
5717 seed1->Reset(kFALSE);
5718 seed0->ResetCovariance(10.);
5719 seed1->ResetCovariance(10.);
5721 first = TMath::Max(row0-20,0);
5722 last = TMath::Min(row0+20,158);
5724 for (Int_t irow=0; irow<20;irow++){
5725 rows[irow] = first +((last-first)*irow)/19;
5727 // store parameters along the track
5729 for (Int_t irow=0;irow<20;irow++){
5730 FollowBackProlongation(*seed0, rows[irow]);
5731 FollowProlongation(*seed1,rows[19-irow]);
5732 param0[irow] = *seed0;
5733 param1[19-irow] = *seed1;
5737 for (Int_t irow=0; irow<19;irow++){
5738 kinks[irow].SetMother(param0[irow]);
5739 kinks[irow].SetDaughter(param1[irow]);
5740 // param0[irow].Dump();
5741 //param1[irow].Dump();
5742 kinks[irow].Update();
5745 // choose kink with biggest change of angle
5748 for (Int_t irow=0;irow<20;irow++){
5749 if (TMath::Abs(kinks[irow].GetR())>250.) continue;
5750 if (TMath::Abs(kinks[irow].GetR())<90.) continue;
5751 Float_t quality = TMath::Abs(kinks[irow].GetAngle(2))/(3.+TMath::Abs(kinks[irow].GetR()-param0[irow].GetX()));
5752 if ( quality > maxchange){
5753 maxchange = quality;
5760 if (index==-1 || param0[index].GetNumberOfClusters()+param1[index].GetNumberOfClusters()<100){
5766 // Float_t anglesigma = TMath::Sqrt(param0[index].fC22+param0[index].fC33+param1[index].fC22+param1[index].fC33);
5768 kink.SetMother(param0[index]);
5769 kink.SetDaughter(param1[index]);
5772 Double_t chi2P2 = param0[index].GetParameter()[2]-param1[index].GetParameter()[2];
5774 chi2P2/=param0[index].GetCovariance()[5]+param1[index].GetCovariance()[5];
5775 Double_t chi2P3 = param0[index].GetParameter()[3]-param1[index].GetParameter()[3];
5777 chi2P3/=param0[index].GetCovariance()[9]+param1[index].GetCovariance()[9];
5779 if (AliTPCReconstructor::StreamLevel()>1) {
5780 (*fDebugStreamer)<<"kinkHpt"<<
5783 "p0.="<<¶m0[index]<<
5784 "p1.="<<¶m1[index]<<
5788 if ( chi2P2+chi2P3<AliTPCReconstructor::GetRecoParam()->GetKinkAngleCutChi2(0)){
5795 row0 = GetRowNumber(kink.GetR());
5796 kink.SetTPCRow0(row0);
5797 kink.SetLabel(CookLabel(seed0,0.5,0,row0),0);
5798 kink.SetLabel(CookLabel(seed1,0.5,row0,158),1);
5799 kink.SetIndex(-10,0);
5800 kink.SetIndex(int(param0[index].GetNumberOfClusters()+param1[index].GetNumberOfClusters()),1);
5801 kink.SetTPCncls(param0[index].GetNumberOfClusters(),0);
5802 kink.SetTPCncls(param1[index].GetNumberOfClusters(),1);
5805 // new (&mother) AliTPCseed(param0[index]);
5806 daughter = param1[index];
5807 daughter.SetLabel(kink.GetLabel(1));
5808 param0[index].Reset(kTRUE);
5809 FollowProlongation(param0[index],0);
5810 mother = param0[index];
5811 mother.SetLabel(kink.GetLabel(0));
5812 if ( chi2P2+chi2P3<AliTPCReconstructor::GetRecoParam()->GetKinkAngleCutChi2(1)){
5824 AliTPCseed* AliTPCtrackerMI::ReSeed(AliTPCseed *t)
5827 // reseed - refit - track
5830 // Int_t last = fSectors->GetNRows()-1;
5832 if (fSectors == fOuterSec){
5833 first = TMath::Max(first, t->GetFirstPoint()-fInnerSec->GetNRows());
5837 first = t->GetFirstPoint();
5839 AliTPCseed * seed = MakeSeed(t,0.1,0.5,0.9);
5840 FollowBackProlongation(*t,fSectors->GetNRows()-1);
5842 FollowProlongation(*t,first);
5852 //_____________________________________________________________________________
5853 Int_t AliTPCtrackerMI::ReadSeeds(const TFile *inp) {
5854 //-----------------------------------------------------------------
5855 // This function reades track seeds.
5856 //-----------------------------------------------------------------
5857 TDirectory *savedir=gDirectory;
5859 TFile *in=(TFile*)inp;
5860 if (!in->IsOpen()) {
5861 cerr<<"AliTPCtrackerMI::ReadSeeds(): input file is not open !\n";
5866 TTree *seedTree=(TTree*)in->Get("Seeds");
5868 cerr<<"AliTPCtrackerMI::ReadSeeds(): ";
5869 cerr<<"can't get a tree with track seeds !\n";
5872 AliTPCtrack *seed=new AliTPCtrack;
5873 seedTree->SetBranchAddress("tracks",&seed);
5875 if (fSeeds==0) fSeeds=new TObjArray(15000);
5877 Int_t n=(Int_t)seedTree->GetEntries();
5878 for (Int_t i=0; i<n; i++) {
5879 seedTree->GetEvent(i);
5880 fSeeds->AddLast(new AliTPCseed(*seed/*,seed->GetAlpha()*/));
5889 Int_t AliTPCtrackerMI::Clusters2Tracks (AliESDEvent *const esd)
5892 // clusters to tracks
5894 if (fSeeds) DeleteSeeds();
5897 AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ;
5898 transform->SetCurrentTimeStamp( esd->GetTimeStamp());
5899 transform->SetCurrentRun(esd->GetRunNumber());
5902 if (!fSeeds) return 1;
5904 if (AliTPCReconstructor::StreamLevel()>3) DumpClusters(0,fSeeds);
5910 //_____________________________________________________________________________
5911 Int_t AliTPCtrackerMI::Clusters2Tracks() {
5912 //-----------------------------------------------------------------
5913 // This is a track finder.
5914 //-----------------------------------------------------------------
5915 TDirectory *savedir=gDirectory;
5919 fSeeds = Tracking();
5922 Info("Clusters2Tracks","Time for tracking: \t");timer.Print();timer.Start();
5924 //activate again some tracks
5925 for (Int_t i=0; i<fSeeds->GetEntriesFast(); i++) {
5926 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
5928 Int_t nc=t.GetNumberOfClusters();
5930 delete fSeeds->RemoveAt(i);
5934 if (pt->GetRemoval()==10) {
5935 if (pt->GetDensityFirst(20)>0.8 || pt->GetDensityFirst(30)>0.8 || pt->GetDensityFirst(40)>0.7)
5936 pt->Desactivate(10); // make track again active
5938 pt->Desactivate(20);
5939 delete fSeeds->RemoveAt(i);
5944 RemoveUsed2(fSeeds,0.85,0.85,0);
5945 if (AliTPCReconstructor::GetRecoParam()->GetDoKinks()) FindKinks(fSeeds,fEvent);
5946 //FindCurling(fSeeds, fEvent,0);
5947 if (AliTPCReconstructor::StreamLevel()>5) FindMultiMC(fSeeds, fEvent,-1); // find multi found tracks
5948 RemoveUsed2(fSeeds,0.5,0.4,20);
5949 FindSplitted(fSeeds, fEvent,0); // find multi found tracks
5950 if (AliTPCReconstructor::StreamLevel()>5) FindMultiMC(fSeeds, fEvent,0); // find multi found tracks
5953 // // refit short tracks
5955 Int_t nseed=fSeeds->GetEntriesFast();
5958 for (Int_t i=0; i<nseed; i++) {
5959 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
5961 Int_t nc=t.GetNumberOfClusters();
5963 delete fSeeds->RemoveAt(i);
5966 CookLabel(pt,0.1); //For comparison only
5967 //if ((pt->IsActive() || (pt->fRemoval==10) )&& nc>50 &&pt->GetNumberOfClusters()>0.4*pt->fNFoundable){
5968 if ((pt->IsActive() || (pt->GetRemoval()==10) )){
5970 if (fDebug>0) cerr<<found<<'\r';
5974 delete fSeeds->RemoveAt(i);
5978 //RemoveOverlap(fSeeds,0.99,7,kTRUE);
5980 //RemoveUsed(fSeeds,0.9,0.9,6);
5982 nseed=fSeeds->GetEntriesFast();
5984 for (Int_t i=0; i<nseed; i++) {
5985 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
5987 Int_t nc=t.GetNumberOfClusters();
5989 delete fSeeds->RemoveAt(i);
5993 t.CookdEdx(0.02,0.6);
5994 // CheckKinkPoint(&t,0.05);
5995 //if ((pt->IsActive() || (pt->fRemoval==10) )&& nc>50 &&pt->GetNumberOfClusters()>0.4*pt->fNFoundable){
5996 if ((pt->IsActive() || (pt->GetRemoval()==10) )){
6004 delete fSeeds->RemoveAt(i);
6005 //AliTPCseed * seed1 = ReSeed(pt,0.05,0.5,1);
6007 // FollowProlongation(*seed1,0);
6008 // Int_t n = seed1->GetNumberOfClusters();
6009 // printf("fP4\t%f\t%f\n",seed1->GetC(),pt->GetC());
6010 // printf("fN\t%d\t%d\n", seed1->GetNumberOfClusters(),pt->GetNumberOfClusters());
6013 //AliTPCseed * seed2 = ReSeed(pt,0.95,0.5,0.05);
6017 SortTracks(fSeeds, 1);
6021 PrepareForBackProlongation(fSeeds,5.);
6022 PropagateBack(fSeeds);
6023 printf("Time for back propagation: \t");timer.Print();timer.Start();
6027 PrepareForProlongation(fSeeds,5.);
6028 PropagateForard2(fSeeds);
6030 printf("Time for FORWARD propagation: \t");timer.Print();timer.Start();
6031 // RemoveUsed(fSeeds,0.7,0.7,6);
6032 //RemoveOverlap(fSeeds,0.9,7,kTRUE);
6034 nseed=fSeeds->GetEntriesFast();
6036 for (Int_t i=0; i<nseed; i++) {
6037 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
6039 Int_t nc=t.GetNumberOfClusters();
6041 delete fSeeds->RemoveAt(i);
6044 t.CookdEdx(0.02,0.6);
6045 // CookLabel(pt,0.1); //For comparison only
6046 //if ((pt->IsActive() || (pt->fRemoval==10) )&& nc>50 &&pt->GetNumberOfClusters()>0.4*pt->fNFoundable){
6047 if ((pt->IsActive() || (pt->fRemoval==10) )){
6048 cerr<<found++<<'\r';
6051 delete fSeeds->RemoveAt(i);
6056 // fNTracks = found;
6058 Info("Clusters2Tracks","Time for overlap removal, track writing and dedx cooking: \t"); timer.Print();timer.Start();
6061 // cerr<<"Number of found tracks : "<<"\t"<<found<<endl;
6062 Info("Clusters2Tracks","Number of found tracks %d",found);
6064 // UnloadClusters();
6069 void AliTPCtrackerMI::Tracking(TObjArray * arr)
6072 // tracking of the seeds
6075 fSectors = fOuterSec;
6076 ParallelTracking(arr,150,63);
6077 fSectors = fOuterSec;
6078 ParallelTracking(arr,63,0);
6081 TObjArray * AliTPCtrackerMI::Tracking(Int_t seedtype, Int_t i1, Int_t i2, Float_t cuts[4], Float_t dy, Int_t dsec)
6086 TObjArray * arr = new TObjArray;
6088 fSectors = fOuterSec;
6091 for (Int_t sec=0;sec<fkNOS;sec++){
6092 if (seedtype==3) MakeSeeds3(arr,sec,i1,i2,cuts,dy, dsec);
6093 if (seedtype==4) MakeSeeds5(arr,sec,i1,i2,cuts,dy);
6094 if (seedtype==2) MakeSeeds2(arr,sec,i1,i2,cuts,dy);
6097 Info("Tracking","\nSeeding - %d\t%d\t%d\t%d\n",seedtype,i1,i2,arr->GetEntriesFast());
6109 TObjArray * AliTPCtrackerMI::Tracking()
6113 if (AliTPCReconstructor::GetRecoParam()->GetSpecialSeeding()) return TrackingSpecial();
6116 Int_t nup=fOuterSec->GetNRows()+fInnerSec->GetNRows();
6118 TObjArray * seeds = new TObjArray;
6120 Int_t fLastSeedRowSec=AliTPCReconstructor::GetRecoParam()->GetLastSeedRowSec();
6121 Int_t gapPrim = AliTPCReconstructor::GetRecoParam()->GetSeedGapPrim();
6122 Int_t gapSec = AliTPCReconstructor::GetRecoParam()->GetSeedGapSec();
6130 Float_t fnumber = 3.0;
6131 Float_t fdensity = 3.0;
6136 for (Int_t delta = 0; delta<18; delta+=gapPrim){
6140 arr = Tracking(3,nup-1-delta,nup-1-delta-gap,cuts,-1,1);
6141 SumTracks(seeds,arr);
6142 SignClusters(seeds,fnumber,fdensity);
6144 for (Int_t i=2;i<6;i+=2){
6145 // seed high pt tracks
6148 arr = Tracking(3,nup-i-delta,nup-i-delta-gap,cuts,-1,0);
6149 SumTracks(seeds,arr);
6150 SignClusters(seeds,fnumber,fdensity);
6155 // RemoveUsed(seeds,0.9,0.9,1);
6156 // UnsignClusters();
6157 // SignClusters(seeds,fnumber,fdensity);
6161 for (Int_t delta = 20; delta<120; delta+=gapPrim){
6163 // seed high pt tracks
6167 arr = Tracking(3,nup-delta,nup-delta-gap,cuts,-1);
6168 SumTracks(seeds,arr);
6169 SignClusters(seeds,fnumber,fdensity);
6174 arr = Tracking(3,nup-delta-5,nup-delta-5-gap,cuts,-1);
6175 SumTracks(seeds,arr);
6176 SignClusters(seeds,fnumber,fdensity);
6187 Info("Tracking()","\n\nPrimary seeding\t%d\n\n",seeds->GetEntriesFast());
6191 // RemoveUsed(seeds,0.75,0.75,1);
6193 //SignClusters(seeds,fnumber,fdensity);
6202 arr = Tracking(4,nup-1,nup-1-gap,cuts,-1);
6203 SumTracks(seeds,arr);
6204 SignClusters(seeds,fnumber,fdensity);
6206 arr = Tracking(4,nup-2,nup-2-gap,cuts,-1);
6207 SumTracks(seeds,arr);
6208 SignClusters(seeds,fnumber,fdensity);
6210 arr = Tracking(4,nup-3,nup-3-gap,cuts,-1);
6211 SumTracks(seeds,arr);
6212 SignClusters(seeds,fnumber,fdensity);
6214 arr = Tracking(4,nup-5,nup-5-gap,cuts,-1);
6215 SumTracks(seeds,arr);
6216 SignClusters(seeds,fnumber,fdensity);
6218 arr = Tracking(4,nup-7,nup-7-gap,cuts,-1);
6219 SumTracks(seeds,arr);
6220 SignClusters(seeds,fnumber,fdensity);
6223 arr = Tracking(4,nup-9,nup-9-gap,cuts,-1);
6224 SumTracks(seeds,arr);
6225 SignClusters(seeds,fnumber,fdensity);
6229 for (Int_t delta = 9; delta<30; delta+=gapSec){
6235 arr = Tracking(4,nup-1-delta,nup-1-delta-gap,cuts,-1);
6236 SumTracks(seeds,arr);
6237 SignClusters(seeds,fnumber,fdensity);
6239 arr = Tracking(4,nup-3-delta,nup-5-delta-gap,cuts,4);
6240 SumTracks(seeds,arr);
6241 SignClusters(seeds,fnumber,fdensity);
6254 for (Int_t delta = 30; delta<fLastSeedRowSec; delta+=gapSec){
6260 arr = Tracking(4,nup-1-delta,nup-1-delta-gap,cuts,-1);
6261 SumTracks(seeds,arr);
6262 SignClusters(seeds,fnumber,fdensity);
6264 arr = Tracking(4,nup-5-delta,nup-5-delta-gap,cuts,5 );
6265 SumTracks(seeds,arr);
6266 SignClusters(seeds,fnumber,fdensity);
6270 Info("Tracking()","\n\nSecondary seeding\t%d\n\n",seeds->GetEntriesFast());
6281 TObjArray * AliTPCtrackerMI::TrackingSpecial()
6284 // seeding adjusted for laser and cosmic tests - short tracks with big inclination angle
6285 // no primary vertex seeding tried
6289 Int_t nup=fOuterSec->GetNRows()+fInnerSec->GetNRows();
6291 TObjArray * seeds = new TObjArray;
6296 Float_t fnumber = 3.0;
6297 Float_t fdensity = 3.0;
6300 cuts[0] = AliTPCReconstructor::GetRecoParam()->GetMaxC(); // max curvature
6301 cuts[1] = 3.5; // max tan(phi) angle for seeding
6302 cuts[2] = 3.; // not used (cut on z primary vertex)
6303 cuts[3] = 3.5; // max tan(theta) angle for seeding
6305 for (Int_t delta = 0; nup-delta-gap-1>0; delta+=3){
6307 arr = Tracking(4,nup-1-delta,nup-1-delta-gap,cuts,-1);
6308 SumTracks(seeds,arr);
6309 SignClusters(seeds,fnumber,fdensity);
6313 Info("Tracking()","\n\nSecondary seeding\t%d\n\n",seeds->GetEntriesFast());
6324 void AliTPCtrackerMI::SumTracks(TObjArray *arr1,TObjArray *&arr2) const
6327 //sum tracks to common container
6328 //remove suspicious tracks
6329 Int_t nseed = arr2->GetEntriesFast();
6330 for (Int_t i=0;i<nseed;i++){
6331 AliTPCseed *pt=(AliTPCseed*)arr2->UncheckedAt(i);
6334 // remove tracks with too big curvature
6336 if (TMath::Abs(pt->GetC())>AliTPCReconstructor::GetRecoParam()->GetMaxC()){
6337 delete arr2->RemoveAt(i);
6340 // REMOVE VERY SHORT TRACKS
6341 if (pt->GetNumberOfClusters()<20){
6342 delete arr2->RemoveAt(i);
6345 // NORMAL ACTIVE TRACK
6346 if (pt->IsActive()){
6347 arr1->AddLast(arr2->RemoveAt(i));
6350 //remove not usable tracks
6351 if (pt->GetRemoval()!=10){
6352 delete arr2->RemoveAt(i);
6356 // ENABLE ONLY ENOUGH GOOD STOPPED TRACKS
6357 if (pt->GetDensityFirst(20)>0.8 || pt->GetDensityFirst(30)>0.8 || pt->GetDensityFirst(40)>0.7)
6358 arr1->AddLast(arr2->RemoveAt(i));
6360 delete arr2->RemoveAt(i);
6364 delete arr2; arr2 = 0;
6369 void AliTPCtrackerMI::ParallelTracking(TObjArray *const arr, Int_t rfirst, Int_t rlast)
6372 // try to track in parralel
6374 Int_t nseed=arr->GetEntriesFast();
6375 //prepare seeds for tracking
6376 for (Int_t i=0; i<nseed; i++) {
6377 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i), &t=*pt;
6379 if (!t.IsActive()) continue;
6380 // follow prolongation to the first layer
6381 if ( (fSectors ==fInnerSec) || (t.GetFirstPoint()-fkParam->GetNRowLow()>rfirst+1) )
6382 FollowProlongation(t, rfirst+1);
6387 for (Int_t nr=rfirst; nr>=rlast; nr--){
6388 if (nr<fInnerSec->GetNRows())
6389 fSectors = fInnerSec;
6391 fSectors = fOuterSec;
6392 // make indexes with the cluster tracks for given
6394 // find nearest cluster
6395 for (Int_t i=0; i<nseed; i++) {
6396 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i), &t=*pt;
6398 if (nr==80) pt->UpdateReference();
6399 if (!pt->IsActive()) continue;
6400 // if ( (fSectors ==fOuterSec) && (pt->fFirstPoint-fkParam->GetNRowLow())<nr) continue;
6401 if (pt->GetRelativeSector()>17) {
6404 UpdateClusters(t,nr);
6406 // prolonagate to the nearest cluster - if founded
6407 for (Int_t i=0; i<nseed; i++) {
6408 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
6410 if (!pt->IsActive()) continue;
6411 // if ((fSectors ==fOuterSec) && (pt->fFirstPoint-fkParam->GetNRowLow())<nr) continue;
6412 if (pt->GetRelativeSector()>17) {
6415 FollowToNextCluster(*pt,nr);
6420 void AliTPCtrackerMI::PrepareForBackProlongation(const TObjArray *const arr,Float_t fac) const
6424 // if we use TPC track itself we have to "update" covariance
6426 Int_t nseed= arr->GetEntriesFast();
6427 for (Int_t i=0;i<nseed;i++){
6428 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
6432 //rotate to current local system at first accepted point
6433 Int_t index = pt->GetClusterIndex2(pt->GetFirstPoint());
6434 Int_t sec = (index&0xff000000)>>24;
6436 Float_t angle1 = fInnerSec->GetAlpha()*sec+fInnerSec->GetAlphaShift();
6437 if (angle1>TMath::Pi())
6438 angle1-=2.*TMath::Pi();
6439 Float_t angle2 = pt->GetAlpha();
6441 if (TMath::Abs(angle1-angle2)>0.001){
6442 if (!pt->Rotate(angle1-angle2)) return;
6443 //angle2 = pt->GetAlpha();
6444 //pt->fRelativeSector = pt->GetAlpha()/fInnerSec->GetAlpha();
6445 //if (pt->GetAlpha()<0)
6446 // pt->fRelativeSector+=18;
6447 //sec = pt->fRelativeSector;
6456 void AliTPCtrackerMI::PrepareForProlongation(TObjArray *const arr, Float_t fac) const
6460 // if we use TPC track itself we have to "update" covariance
6462 Int_t nseed= arr->GetEntriesFast();
6463 for (Int_t i=0;i<nseed;i++){
6464 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
6467 pt->SetFirstPoint(pt->GetLastPoint());
6475 Int_t AliTPCtrackerMI::PropagateBack(const TObjArray *const arr)
6478 // make back propagation
6480 Int_t nseed= arr->GetEntriesFast();
6481 for (Int_t i=0;i<nseed;i++){
6482 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
6483 if (pt&& pt->GetKinkIndex(0)<=0) {
6484 //AliTPCseed *pt2 = new AliTPCseed(*pt);
6485 fSectors = fInnerSec;
6486 //FollowBackProlongation(*pt,fInnerSec->GetNRows()-1);
6487 //fSectors = fOuterSec;
6488 FollowBackProlongation(*pt,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1,1);
6489 //if (pt->GetNumberOfClusters()<(pt->fEsd->GetTPCclusters(0)) ){
6490 // Error("PropagateBack","Not prolonged track %d",pt->GetLabel());
6491 // FollowBackProlongation(*pt2,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1);
6494 if (pt&& pt->GetKinkIndex(0)>0) {
6495 AliESDkink * kink = fEvent->GetKink(pt->GetKinkIndex(0)-1);
6496 pt->SetFirstPoint(kink->GetTPCRow0());
6497 fSectors = fInnerSec;
6498 FollowBackProlongation(*pt,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1,1);
6506 Int_t AliTPCtrackerMI::PropagateForward2(const TObjArray *const arr)
6509 // make forward propagation
6511 Int_t nseed= arr->GetEntriesFast();
6513 for (Int_t i=0;i<nseed;i++){
6514 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
6516 FollowProlongation(*pt,0,1,1);
6525 Int_t AliTPCtrackerMI::PropagateForward()
6528 // propagate track forward
6530 Int_t nseed = fSeeds->GetEntriesFast();
6531 for (Int_t i=0;i<nseed;i++){
6532 AliTPCseed *pt = (AliTPCseed*)fSeeds->UncheckedAt(i);
6534 AliTPCseed &t = *pt;
6535 Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
6536 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
6537 if (alpha < 0. ) alpha += 2.*TMath::Pi();
6538 t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN);
6542 fSectors = fOuterSec;
6543 ParallelTracking(fSeeds,fOuterSec->GetNRows()+fInnerSec->GetNRows()-1,fInnerSec->GetNRows());
6544 fSectors = fInnerSec;
6545 ParallelTracking(fSeeds,fInnerSec->GetNRows()-1,0);
6554 Int_t AliTPCtrackerMI::PropagateBack(AliTPCseed *const pt, Int_t row0, Int_t row1)
6557 // make back propagation, in between row0 and row1
6561 fSectors = fInnerSec;
6564 if (row1<fSectors->GetNRows())
6567 r1 = fSectors->GetNRows()-1;
6569 if (row0<fSectors->GetNRows()&& r1>0 )
6570 FollowBackProlongation(*pt,r1);
6571 if (row1<=fSectors->GetNRows())
6574 r1 = row1 - fSectors->GetNRows();
6575 if (r1<=0) return 0;
6576 if (r1>=fOuterSec->GetNRows()) return 0;
6577 fSectors = fOuterSec;
6578 return FollowBackProlongation(*pt,r1);
6586 void AliTPCtrackerMI::GetShape(AliTPCseed * seed, Int_t row)
6588 // gets cluster shape
6590 AliTPCClusterParam * clparam = AliTPCcalibDB::Instance()->GetClusterParam();
6591 Float_t zdrift = TMath::Abs((fkParam->GetZLength(0)-TMath::Abs(seed->GetZ())));
6592 Int_t type = (seed->GetSector() < fkParam->GetNSector()/2) ? 0: (row>126) ? 1:2;
6593 Double_t angulary = seed->GetSnp();
6595 if (TMath::Abs(angulary)>AliTPCReconstructor::GetMaxSnpTracker()) {
6596 angulary = TMath::Sign(AliTPCReconstructor::GetMaxSnpTracker(),angulary);
6599 angulary = angulary*angulary/((1.-angulary)*(1.+angulary));
6600 Double_t angularz = seed->GetTgl()*seed->GetTgl()*(1.+angulary);
6602 Double_t sigmay = clparam->GetRMS0(0,type,zdrift,TMath::Sqrt(TMath::Abs(angulary)));
6603 Double_t sigmaz = clparam->GetRMS0(1,type,zdrift,TMath::Sqrt(TMath::Abs(angularz)));
6604 seed->SetCurrentSigmaY2(sigmay*sigmay);
6605 seed->SetCurrentSigmaZ2(sigmaz*sigmaz);
6606 // Float_t sd2 = TMath::Abs((fkParam->GetZLength(0)-TMath::Abs(seed->GetZ())))*fkParam->GetDiffL()*fkParam->GetDiffL();
6607 // // Float_t padlength = fkParam->GetPadPitchLength(seed->fSector);
6608 // Float_t padlength = GetPadPitchLength(row);
6610 // Float_t sresy = (seed->GetSector() < fkParam->GetNSector()/2) ? 0.2 :0.3;
6611 // seed->SetCurrentSigmaY2(sd2+padlength*padlength*angulary/12.+sresy*sresy);
6613 // Float_t sresz = fkParam->GetZSigma();
6614 // seed->SetCurrentSigmaZ2(sd2+padlength*padlength*angularz*angularz*(1+angulary)/12.+sresz*sresz);
6616 Float_t wy = GetSigmaY(seed);
6617 Float_t wz = GetSigmaZ(seed);
6620 if (TMath::Abs(wy/seed->fCurrentSigmaY2-1)>0.0001 || TMath::Abs(wz/seed->fCurrentSigmaZ2-1)>0.0001 ){
6621 printf("problem\n");
6628 //__________________________________________________________________________
6629 void AliTPCtrackerMI::CookLabel(AliKalmanTrack *tk, Float_t wrong) const {
6630 //--------------------------------------------------------------------
6631 //This function "cooks" a track label. If label<0, this track is fake.
6632 //--------------------------------------------------------------------
6633 AliTPCseed * t = dynamic_cast<AliTPCseed*>(tk);
6635 printf("%s:%d wrong type \n",(char*)__FILE__,__LINE__);
6639 Int_t noc=t->GetNumberOfClusters();
6641 //printf("\nnot founded prolongation\n\n\n");
6647 AliTPCclusterMI *clusters[160];
6649 for (Int_t i=0;i<160;i++) {
6656 for (i=0; i<160 && current<noc; i++) {
6658 Int_t index=t->GetClusterIndex2(i);
6659 if (index<=0) continue;
6660 if (index&0x8000) continue;
6662 //clusters[current]=GetClusterMI(index);
6663 if (t->GetClusterPointer(i)){
6664 clusters[current]=t->GetClusterPointer(i);
6670 Int_t lab=123456789;
6671 for (i=0; i<noc; i++) {
6672 AliTPCclusterMI *c=clusters[i];
6674 lab=TMath::Abs(c->GetLabel(0));
6676 for (j=0; j<noc; j++) if (lb[j]==lab || mx[j]==0) break;
6682 for (i=0; i<noc; i++) if (mx[i]>max) {max=mx[i]; lab=lb[i];}
6684 for (i=0; i<noc; i++) {
6685 AliTPCclusterMI *c=clusters[i];
6687 if (TMath::Abs(c->GetLabel(1)) == lab ||
6688 TMath::Abs(c->GetLabel(2)) == lab ) max++;
6690 if (noc<=0) { lab=-1; return;}
6691 if ((1.- Float_t(max)/(noc)) > wrong) lab=-lab;
6694 Int_t tail=Int_t(0.10*noc);
6697 for (i=1; i<160&&ind<tail; i++) {
6698 // AliTPCclusterMI *c=clusters[noc-i];
6699 AliTPCclusterMI *c=clusters[i];
6701 if (lab == TMath::Abs(c->GetLabel(0)) ||
6702 lab == TMath::Abs(c->GetLabel(1)) ||
6703 lab == TMath::Abs(c->GetLabel(2))) max++;
6706 if (max < Int_t(0.5*tail)) lab=-lab;
6713 //delete[] clusters;
6717 //__________________________________________________________________________
6718 Int_t AliTPCtrackerMI::CookLabel(AliTPCseed *const t, Float_t wrong,Int_t first, Int_t last) const {
6719 //--------------------------------------------------------------------
6720 //This function "cooks" a track label. If label<0, this track is fake.
6721 //--------------------------------------------------------------------
6722 Int_t noc=t->GetNumberOfClusters();
6724 //printf("\nnot founded prolongation\n\n\n");
6730 AliTPCclusterMI *clusters[160];
6732 for (Int_t i=0;i<160;i++) {
6739 for (i=0; i<160 && current<noc; i++) {
6740 if (i<first) continue;
6741 if (i>last) continue;
6742 Int_t index=t->GetClusterIndex2(i);
6743 if (index<=0) continue;
6744 if (index&0x8000) continue;
6746 //clusters[current]=GetClusterMI(index);
6747 if (t->GetClusterPointer(i)){
6748 clusters[current]=t->GetClusterPointer(i);
6753 //if (noc<5) return -1;
6754 Int_t lab=123456789;
6755 for (i=0; i<noc; i++) {
6756 AliTPCclusterMI *c=clusters[i];
6758 lab=TMath::Abs(c->GetLabel(0));
6760 for (j=0; j<noc; j++) if (lb[j]==lab || mx[j]==0) break;
6766 for (i=0; i<noc; i++) if (mx[i]>max) {max=mx[i]; lab=lb[i];}
6768 for (i=0; i<noc; i++) {
6769 AliTPCclusterMI *c=clusters[i];
6771 if (TMath::Abs(c->GetLabel(1)) == lab ||
6772 TMath::Abs(c->GetLabel(2)) == lab ) max++;
6774 if (noc<=0) { lab=-1; return -1;}
6775 if ((1.- Float_t(max)/(noc)) > wrong) lab=-lab;
6778 Int_t tail=Int_t(0.10*noc);
6781 for (i=1; i<160&&ind<tail; i++) {
6782 // AliTPCclusterMI *c=clusters[noc-i];
6783 AliTPCclusterMI *c=clusters[i];
6785 if (lab == TMath::Abs(c->GetLabel(0)) ||
6786 lab == TMath::Abs(c->GetLabel(1)) ||
6787 lab == TMath::Abs(c->GetLabel(2))) max++;
6790 if (max < Int_t(0.5*tail)) lab=-lab;
6793 // t->SetLabel(lab);
6797 //delete[] clusters;
6801 Int_t AliTPCtrackerMI::GetRowNumber(Double_t x[3]) const
6803 //return pad row number for given x vector
6804 Float_t phi = TMath::ATan2(x[1],x[0]);
6805 if(phi<0) phi=2.*TMath::Pi()+phi;
6806 // Get the local angle in the sector philoc
6807 const Float_t kRaddeg = 180/3.14159265358979312;
6808 Float_t phiangle = (Int_t (phi*kRaddeg/20.) + 0.5)*20./kRaddeg;
6809 Double_t localx = x[0]*TMath::Cos(phiangle)-x[1]*TMath::Sin(phiangle);
6810 return GetRowNumber(localx);
6815 void AliTPCtrackerMI::MakeBitmaps(AliTPCseed *t)
6817 //-----------------------------------------------------------------------
6818 // Fill the cluster and sharing bitmaps of the track
6819 //-----------------------------------------------------------------------
6821 Int_t firstpoint = 0;
6822 Int_t lastpoint = 159;
6823 AliTPCTrackerPoint *point;
6824 AliTPCclusterMI *cluster;
6826 for (int iter=firstpoint; iter<lastpoint; iter++) {
6827 // Change to cluster pointers to see if we have a cluster at given padrow
6828 cluster = t->GetClusterPointer(iter);
6830 t->SetClusterMapBit(iter, kTRUE);
6831 point = t->GetTrackPoint(iter);
6832 if (point->IsShared())
6833 t->SetSharedMapBit(iter,kTRUE);
6835 t->SetSharedMapBit(iter, kFALSE);
6838 t->SetClusterMapBit(iter, kFALSE);
6839 t->SetSharedMapBit(iter, kFALSE);
6844 Bool_t AliTPCtrackerMI::IsFindable(AliTPCseed & track){
6846 // return flag if there is findable cluster at given position
6849 Float_t z = track.GetZ();
6851 if (TMath::Abs(z)<(AliTPCReconstructor::GetCtgRange()*track.GetX()+kDeltaZ) &&
6852 TMath::Abs(z)<fkParam->GetZLength(0) &&
6853 (TMath::Abs(track.GetSnp())<AliTPCReconstructor::GetMaxSnpTracker()))
6859 void AliTPCtrackerMI::AddCovariance(AliTPCseed * seed){
6861 // Adding systematic error estimate to the covariance matrix
6862 // !!!! the systematic error for element 4 is in 1/cm not in pt
6864 const Double_t *param = AliTPCReconstructor::GetRecoParam()->GetSystematicError();
6866 // use only the diagonal part if not specified otherwise
6867 if (!AliTPCReconstructor::GetRecoParam()->GetUseSystematicCorrelation()) return AddCovarianceAdd(seed);
6869 Double_t *covarS= (Double_t*)seed->GetCovariance();
6870 Double_t factor[5]={1,1,1,1,1};
6871 Double_t facC = AliTracker::GetBz()*kB2C;
6872 factor[0]= TMath::Sqrt(TMath::Abs((covarS[0] + param[0]*param[0])/covarS[0]));
6873 factor[1]= TMath::Sqrt(TMath::Abs((covarS[2] + param[1]*param[1])/covarS[2]));
6874 factor[2]= TMath::Sqrt(TMath::Abs((covarS[5] + param[2]*param[2])/covarS[5]));
6875 factor[3]= TMath::Sqrt(TMath::Abs((covarS[9] + param[3]*param[3])/covarS[9]));
6876 factor[4]= TMath::Sqrt(TMath::Abs((covarS[14] + facC*facC*param[4]*param[4])/covarS[14]));
6882 for (Int_t i=0; i<5; i++){
6883 for (Int_t j=i; j<5; j++){
6884 Int_t index=seed->GetIndex(i,j);
6885 covarS[index]*=factor[i]*factor[j];
6891 void AliTPCtrackerMI::AddCovarianceAdd(AliTPCseed * seed){
6893 // Adding systematic error - as additive factor without correlation
6895 const Double_t *param = AliTPCReconstructor::GetRecoParam()->GetSystematicError();
6896 Double_t *covarIn= (Double_t*)seed->GetCovariance();
6898 for (Int_t i=0;i<15;i++) covar[i]=0;
6904 covar[0] = param[0]*param[0];
6905 covar[2] = param[1]*param[1];
6906 covar[5] = param[2]*param[2];
6907 covar[9] = param[3]*param[3];
6908 Double_t facC = AliTracker::GetBz()*kB2C;
6909 covar[14]= param[4]*param[4]*facC*facC;
6911 covar[1]=TMath::Sqrt((covar[0]*covar[2]))*covarIn[1]/TMath::Sqrt((covarIn[0]*covarIn[2]));
6913 covar[3]=TMath::Sqrt((covar[0]*covar[5]))*covarIn[3]/TMath::Sqrt((covarIn[0]*covarIn[5]));
6914 covar[4]=TMath::Sqrt((covar[2]*covar[5]))*covarIn[4]/TMath::Sqrt((covarIn[2]*covarIn[5]));
6916 covar[6]=TMath::Sqrt((covar[0]*covar[9]))*covarIn[6]/TMath::Sqrt((covarIn[0]*covarIn[9]));
6917 covar[7]=TMath::Sqrt((covar[2]*covar[9]))*covarIn[7]/TMath::Sqrt((covarIn[2]*covarIn[9]));
6918 covar[8]=TMath::Sqrt((covar[5]*covar[9]))*covarIn[8]/TMath::Sqrt((covarIn[5]*covarIn[9]));
6920 covar[10]=TMath::Sqrt((covar[0]*covar[14]))*covarIn[10]/TMath::Sqrt((covarIn[0]*covarIn[14]));
6921 covar[11]=TMath::Sqrt((covar[2]*covar[14]))*covarIn[11]/TMath::Sqrt((covarIn[2]*covarIn[14]));
6922 covar[12]=TMath::Sqrt((covar[5]*covar[14]))*covarIn[12]/TMath::Sqrt((covarIn[5]*covarIn[14]));
6923 covar[13]=TMath::Sqrt((covar[9]*covar[14]))*covarIn[13]/TMath::Sqrt((covarIn[9]*covarIn[14]));
6925 seed->AddCovariance(covar);