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"
144 #include "AliTPCdEdxInfo.h"
148 ClassImp(AliTPCtrackerMI)
152 class AliTPCFastMath {
155 static Double_t FastAsin(Double_t x);
157 static Double_t fgFastAsin[20000]; //lookup table for fast asin computation
160 Double_t AliTPCFastMath::fgFastAsin[20000];
161 AliTPCFastMath gAliTPCFastMath; // needed to fill the LUT
163 AliTPCFastMath::AliTPCFastMath(){
165 // initialized lookup table;
166 for (Int_t i=0;i<10000;i++){
167 fgFastAsin[2*i] = TMath::ASin(i/10000.);
168 fgFastAsin[2*i+1] = (TMath::ASin((i+1)/10000.)-fgFastAsin[2*i]);
172 Double_t AliTPCFastMath::FastAsin(Double_t x){
174 // return asin using lookup table
176 Int_t index = int(x*10000);
177 return fgFastAsin[2*index]+(x*10000.-index)*fgFastAsin[2*index+1];
180 Int_t index = int(x*10000);
181 return -(fgFastAsin[2*index]+(x*10000.-index)*fgFastAsin[2*index+1]);
183 //__________________________________________________________________
184 AliTPCtrackerMI::AliTPCtrackerMI()
207 // default constructor
209 for (Int_t irow=0; irow<200; irow++){
216 //_____________________________________________________________________
220 Int_t AliTPCtrackerMI::UpdateTrack(AliTPCseed * track, Int_t accept){
222 //update track information using current cluster - track->fCurrentCluster
225 AliTPCclusterMI* c =track->GetCurrentCluster();
226 if (accept > 0) //sign not accepted clusters
227 track->SetCurrentClusterIndex1(track->GetCurrentClusterIndex1() | 0x8000);
228 else // unsign accpeted clusters
229 track->SetCurrentClusterIndex1(track->GetCurrentClusterIndex1() & 0xffff7fff);
230 UInt_t i = track->GetCurrentClusterIndex1();
232 Int_t sec=(i&0xff000000)>>24;
233 //Int_t row = (i&0x00ff0000)>>16;
234 track->SetRow((i&0x00ff0000)>>16);
235 track->SetSector(sec);
236 // Int_t index = i&0xFFFF;
237 if (sec>=fkParam->GetNInnerSector()) track->SetRow(track->GetRow()+fkParam->GetNRowLow());
238 track->SetClusterIndex2(track->GetRow(), i);
239 //track->fFirstPoint = row;
240 //if ( track->fLastPoint<row) track->fLastPoint =row;
241 // if (track->fRow<0 || track->fRow>160) {
242 // printf("problem\n");
244 if (track->GetFirstPoint()>track->GetRow())
245 track->SetFirstPoint(track->GetRow());
246 if (track->GetLastPoint()<track->GetRow())
247 track->SetLastPoint(track->GetRow());
250 track->SetClusterPointer(track->GetRow(),c);
253 Double_t angle2 = track->GetSnp()*track->GetSnp();
255 //SET NEW Track Point
257 if (angle2<1) //PH sometimes angle2 is very big. To be investigated...
259 angle2 = TMath::Sqrt(angle2/(1-angle2));
260 AliTPCTrackerPoint &point =*(track->GetTrackPoint(track->GetRow()));
262 point.SetSigmaY(c->GetSigmaY2()/track->GetCurrentSigmaY2());
263 point.SetSigmaZ(c->GetSigmaZ2()/track->GetCurrentSigmaZ2());
264 point.SetErrY(sqrt(track->GetErrorY2()));
265 point.SetErrZ(sqrt(track->GetErrorZ2()));
267 point.SetX(track->GetX());
268 point.SetY(track->GetY());
269 point.SetZ(track->GetZ());
270 point.SetAngleY(angle2);
271 point.SetAngleZ(track->GetTgl());
272 if (point.IsShared()){
273 track->SetErrorY2(track->GetErrorY2()*4);
274 track->SetErrorZ2(track->GetErrorZ2()*4);
278 Double_t chi2 = track->GetPredictedChi2(track->GetCurrentCluster());
280 // track->SetErrorY2(track->GetErrorY2()*1.3);
281 // track->SetErrorY2(track->GetErrorY2()+0.01);
282 // track->SetErrorZ2(track->GetErrorZ2()*1.3);
283 // track->SetErrorZ2(track->GetErrorZ2()+0.005);
285 if (accept>0) return 0;
286 if (track->GetNumberOfClusters()%20==0){
287 // if (track->fHelixIn){
288 // TClonesArray & larr = *(track->fHelixIn);
289 // Int_t ihelix = larr.GetEntriesFast();
290 // new(larr[ihelix]) AliHelix(*track) ;
293 track->SetNoCluster(0);
294 return track->Update(c,chi2,i);
299 Int_t AliTPCtrackerMI::AcceptCluster(AliTPCseed * seed, AliTPCclusterMI * cluster)
302 // decide according desired precision to accept given
303 // cluster for tracking
305 seed->GetProlongation(cluster->GetX(),yt,zt);
306 Double_t sy2=ErrY2(seed,cluster);
307 Double_t sz2=ErrZ2(seed,cluster);
309 Double_t sdistancey2 = sy2+seed->GetSigmaY2();
310 Double_t sdistancez2 = sz2+seed->GetSigmaZ2();
311 Double_t dy=seed->GetCurrentCluster()->GetY()-yt;
312 Double_t dz=seed->GetCurrentCluster()->GetZ()-zt;
313 Double_t rdistancey2 = (seed->GetCurrentCluster()->GetY()-yt)*
314 (seed->GetCurrentCluster()->GetY()-yt)/sdistancey2;
315 Double_t rdistancez2 = (seed->GetCurrentCluster()->GetZ()-zt)*
316 (seed->GetCurrentCluster()->GetZ()-zt)/sdistancez2;
318 Double_t rdistance2 = rdistancey2+rdistancez2;
321 if (AliTPCReconstructor::StreamLevel()>2 && seed->GetNumberOfClusters()>20) {
322 Float_t rmsy2 = seed->GetCurrentSigmaY2();
323 Float_t rmsz2 = seed->GetCurrentSigmaZ2();
324 Float_t rmsy2p30 = seed->GetCMeanSigmaY2p30();
325 Float_t rmsz2p30 = seed->GetCMeanSigmaZ2p30();
326 Float_t rmsy2p30R = seed->GetCMeanSigmaY2p30R();
327 Float_t rmsz2p30R = seed->GetCMeanSigmaZ2p30R();
328 AliExternalTrackParam param(*seed);
329 static TVectorD gcl(3),gtr(3);
331 param.GetXYZ(gcl.GetMatrixArray());
332 cluster->GetGlobalXYZ(gclf);
333 gcl[0]=gclf[0]; gcl[1]=gclf[1]; gcl[2]=gclf[2];
336 if (AliTPCReconstructor::StreamLevel()>2) {
337 (*fDebugStreamer)<<"ErrParam"<<
350 "rmsy2p30="<<rmsy2p30<<
351 "rmsz2p30="<<rmsz2p30<<
352 "rmsy2p30R="<<rmsy2p30R<<
353 "rmsz2p30R="<<rmsz2p30R<<
354 // normalize distance -
355 "rdisty="<<rdistancey2<<
356 "rdistz="<<rdistancez2<<
357 "rdist="<<rdistance2<< //
361 //return 0; // temporary
362 if (rdistance2>32) return 3;
365 if ((rdistancey2>9. || rdistancez2>9.) && cluster->GetType()==0)
366 return 2; //suspisiouce - will be changed
368 if ((rdistancey2>6.25 || rdistancez2>6.25) && cluster->GetType()>0)
369 // strict cut on overlaped cluster
370 return 2; //suspisiouce - will be changed
372 if ( (rdistancey2>1. || rdistancez2>6.25 )
373 && cluster->GetType()<0){
374 seed->SetNFoundable(seed->GetNFoundable()-1);
378 if (fUseHLTClusters == 3 || fUseHLTClusters == 4) {
379 if(!AliTPCReconstructor::GetRecoParam()->GetUseHLTOnePadCluster()) {
380 if (TMath::Abs(cluster->GetSigmaY2()) < kAlmost0)
392 //_____________________________________________________________________________
393 AliTPCtrackerMI::AliTPCtrackerMI(const AliTPCParam *par):
395 fkNIS(par->GetNInnerSector()/2),
397 fkNOS(par->GetNOuterSector()/2),
416 //---------------------------------------------------------------------
417 // The main TPC tracker constructor
418 //---------------------------------------------------------------------
419 fInnerSec=new AliTPCtrackerSector[fkNIS];
420 fOuterSec=new AliTPCtrackerSector[fkNOS];
423 for (i=0; i<fkNIS; i++) fInnerSec[i].Setup(par,0);
424 for (i=0; i<fkNOS; i++) fOuterSec[i].Setup(par,1);
427 Int_t nrowlow = par->GetNRowLow();
428 Int_t nrowup = par->GetNRowUp();
431 for (i=0;i<nrowlow;i++){
432 fXRow[i] = par->GetPadRowRadiiLow(i);
433 fPadLength[i]= par->GetPadPitchLength(0,i);
434 fYMax[i] = fXRow[i]*TMath::Tan(0.5*par->GetInnerAngle());
438 for (i=0;i<nrowup;i++){
439 fXRow[i+nrowlow] = par->GetPadRowRadiiUp(i);
440 fPadLength[i+nrowlow] = par->GetPadPitchLength(60,i);
441 fYMax[i+nrowlow] = fXRow[i+nrowlow]*TMath::Tan(0.5*par->GetOuterAngle());
444 if (AliTPCReconstructor::StreamLevel()>0) {
445 fDebugStreamer = new TTreeSRedirector("TPCdebug.root");
448 //________________________________________________________________________
449 AliTPCtrackerMI::AliTPCtrackerMI(const AliTPCtrackerMI &t):
471 //------------------------------------
472 // dummy copy constructor
473 //------------------------------------------------------------------
475 for (Int_t irow=0; irow<200; irow++){
482 AliTPCtrackerMI & AliTPCtrackerMI::operator=(const AliTPCtrackerMI& /*r*/)
484 //------------------------------
486 //--------------------------------------------------------------
489 //_____________________________________________________________________________
490 AliTPCtrackerMI::~AliTPCtrackerMI() {
491 //------------------------------------------------------------------
492 // TPC tracker destructor
493 //------------------------------------------------------------------
500 if (fDebugStreamer) delete fDebugStreamer;
504 void AliTPCtrackerMI::FillESD(const TObjArray* arr)
508 //fill esds using updated tracks
511 // write tracks to the event
512 // store index of the track
513 Int_t nseed=arr->GetEntriesFast();
514 //FindKinks(arr,fEvent);
515 for (Int_t i=0; i<nseed; i++) {
516 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
520 if (AliTPCReconstructor::StreamLevel()>1) {
521 (*fDebugStreamer)<<"Track0"<<
525 // pt->PropagateTo(fkParam->GetInnerRadiusLow());
526 if (pt->GetKinkIndex(0)<=0){ //don't propagate daughter tracks
527 pt->PropagateTo(fkParam->GetInnerRadiusLow());
530 if (( pt->GetPoints()[2]- pt->GetPoints()[0])>5 && pt->GetPoints()[3]>0.8){
532 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
533 iotrack.SetTPCPoints(pt->GetPoints());
534 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
535 iotrack.SetV0Indexes(pt->GetV0Indexes());
536 // iotrack.SetTPCpid(pt->fTPCr);
537 //iotrack.SetTPCindex(i);
538 fEvent->AddTrack(&iotrack);
542 if ( (pt->GetNumberOfClusters()>70)&& (Float_t(pt->GetNumberOfClusters())/Float_t(pt->GetNFoundable()))>0.55) {
544 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
545 iotrack.SetTPCPoints(pt->GetPoints());
546 //iotrack.SetTPCindex(i);
547 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
548 iotrack.SetV0Indexes(pt->GetV0Indexes());
549 // iotrack.SetTPCpid(pt->fTPCr);
550 fEvent->AddTrack(&iotrack);
554 // short tracks - maybe decays
556 if ( (pt->GetNumberOfClusters()>30) && (Float_t(pt->GetNumberOfClusters())/Float_t(pt->GetNFoundable()))>0.70) {
557 Int_t found,foundable,shared;
558 pt->GetClusterStatistic(0,60,found, foundable,shared,kFALSE);
559 if ( (found>20) && (pt->GetNShared()/float(pt->GetNumberOfClusters())<0.2)){
561 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
562 //iotrack.SetTPCindex(i);
563 iotrack.SetTPCPoints(pt->GetPoints());
564 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
565 iotrack.SetV0Indexes(pt->GetV0Indexes());
566 //iotrack.SetTPCpid(pt->fTPCr);
567 fEvent->AddTrack(&iotrack);
572 if ( (pt->GetNumberOfClusters()>20) && (Float_t(pt->GetNumberOfClusters())/Float_t(pt->GetNFoundable()))>0.8) {
573 Int_t found,foundable,shared;
574 pt->GetClusterStatistic(0,60,found, foundable,shared,kFALSE);
575 if (found<20) continue;
576 if (pt->GetNShared()/float(pt->GetNumberOfClusters())>0.2) continue;
579 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
580 iotrack.SetTPCPoints(pt->GetPoints());
581 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
582 iotrack.SetV0Indexes(pt->GetV0Indexes());
583 //iotrack.SetTPCpid(pt->fTPCr);
584 //iotrack.SetTPCindex(i);
585 fEvent->AddTrack(&iotrack);
588 // short tracks - secondaties
590 if ( (pt->GetNumberOfClusters()>30) ) {
591 Int_t found,foundable,shared;
592 pt->GetClusterStatistic(128,158,found, foundable,shared,kFALSE);
593 if ( (found>20) && (pt->GetNShared()/float(pt->GetNumberOfClusters())<0.2) &&float(found)/float(foundable)>0.8){
595 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
596 iotrack.SetTPCPoints(pt->GetPoints());
597 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
598 iotrack.SetV0Indexes(pt->GetV0Indexes());
599 //iotrack.SetTPCpid(pt->fTPCr);
600 //iotrack.SetTPCindex(i);
601 fEvent->AddTrack(&iotrack);
606 if ( (pt->GetNumberOfClusters()>15)) {
607 Int_t found,foundable,shared;
608 pt->GetClusterStatistic(138,158,found, foundable,shared,kFALSE);
609 if (found<15) continue;
610 if (foundable<=0) continue;
611 if (pt->GetNShared()/float(pt->GetNumberOfClusters())>0.2) continue;
612 if (float(found)/float(foundable)<0.8) continue;
615 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
616 iotrack.SetTPCPoints(pt->GetPoints());
617 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
618 iotrack.SetV0Indexes(pt->GetV0Indexes());
619 // iotrack.SetTPCpid(pt->fTPCr);
620 //iotrack.SetTPCindex(i);
621 fEvent->AddTrack(&iotrack);
625 // >> account for suppressed tracks in the kink indices (RS)
626 int nESDtracks = fEvent->GetNumberOfTracks();
627 for (int it=nESDtracks;it--;) {
628 AliESDtrack* esdTr = fEvent->GetTrack(it);
629 if (!esdTr || !esdTr->GetKinkIndex(0)) continue;
630 for (int ik=0;ik<3;ik++) {
632 if (!(knkId=esdTr->GetKinkIndex(ik))) break; // no more kinks for this track
633 AliESDkink* kink = fEvent->GetKink(TMath::Abs(knkId)-1);
635 AliError(Form("ESDTrack%d refers to non-existing kink %d",it,TMath::Abs(knkId)-1));
638 kink->SetIndex(it, knkId<0 ? 0:1); // update track index of the kink: mother at 0, daughter at 1
641 // << account for suppressed tracks in the kink indices (RS)
642 AliInfo(Form("Number of filled ESDs-\t%d\n",fEvent->GetNumberOfTracks()));
650 Double_t AliTPCtrackerMI::ErrY2(AliTPCseed* seed, const AliTPCclusterMI * cl){
653 // Use calibrated cluster error from OCDB
655 AliTPCClusterParam * clparam = AliTPCcalibDB::Instance()->GetClusterParam();
657 Float_t z = TMath::Abs(fkParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
658 Int_t ctype = cl->GetType();
659 Int_t type = (cl->GetRow()<63) ? 0: (cl->GetRow()>126) ? 1:2;
660 Double_t angle = seed->GetSnp()*seed->GetSnp();
661 angle = TMath::Sqrt(TMath::Abs(angle/(1.-angle)));
662 Double_t erry2 = clparam->GetError0Par(0,type, z,angle);
664 erry2+=0.5; // edge cluster
667 seed->SetErrorY2(erry2);
671 //calculate look-up table at the beginning
672 // static Bool_t ginit = kFALSE;
673 // static Float_t gnoise1,gnoise2,gnoise3;
674 // static Float_t ggg1[10000];
675 // static Float_t ggg2[10000];
676 // static Float_t ggg3[10000];
677 // static Float_t glandau1[10000];
678 // static Float_t glandau2[10000];
679 // static Float_t glandau3[10000];
681 // static Float_t gcor01[500];
682 // static Float_t gcor02[500];
683 // static Float_t gcorp[500];
687 // if (ginit==kFALSE){
688 // for (Int_t i=1;i<500;i++){
689 // Float_t rsigma = float(i)/100.;
690 // gcor02[i] = TMath::Max(0.78 +TMath::Exp(7.4*(rsigma-1.2)),0.6);
691 // gcor01[i] = TMath::Max(0.72 +TMath::Exp(3.36*(rsigma-1.2)),0.6);
692 // gcorp[i] = TMath::Max(TMath::Power((rsigma+0.5),1.5),1.2);
696 // for (Int_t i=3;i<10000;i++){
700 // Float_t amp = float(i);
701 // Float_t padlength =0.75;
702 // gnoise1 = 0.0004/padlength;
703 // Float_t nel = 0.268*amp;
704 // Float_t nprim = 0.155*amp;
705 // ggg1[i] = fkParam->GetDiffT()*fkParam->GetDiffT()*(2+0.001*nel/(padlength*padlength))/nel;
706 // glandau1[i] = (2.+0.12*nprim)*0.5* (2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
707 // if (glandau1[i]>1) glandau1[i]=1;
708 // glandau1[i]*=padlength*padlength/12.;
712 // gnoise2 = 0.0004/padlength;
714 // nprim = 0.133*amp;
715 // ggg2[i] = fkParam->GetDiffT()*fkParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
716 // glandau2[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
717 // if (glandau2[i]>1) glandau2[i]=1;
718 // glandau2[i]*=padlength*padlength/12.;
723 // gnoise3 = 0.0004/padlength;
725 // nprim = 0.133*amp;
726 // ggg3[i] = fkParam->GetDiffT()*fkParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
727 // glandau3[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
728 // if (glandau3[i]>1) glandau3[i]=1;
729 // glandau3[i]*=padlength*padlength/12.;
737 // Int_t amp = int(TMath::Abs(cl->GetQ()));
739 // seed->SetErrorY2(1.);
743 // Float_t z = TMath::Abs(fkParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
744 // Int_t ctype = cl->GetType();
745 // Float_t padlength= GetPadPitchLength(seed->GetRow());
746 // Double_t angle2 = seed->GetSnp()*seed->GetSnp();
747 // angle2 = angle2/(1-angle2);
749 // //cluster "quality"
750 // Int_t rsigmay = int(100.*cl->GetSigmaY2()/(seed->GetCurrentSigmaY2()));
753 // if (fSectors==fInnerSec){
754 // snoise2 = gnoise1;
755 // res = ggg1[amp]*z+glandau1[amp]*angle2;
756 // if (ctype==0) res *= gcor01[rsigmay];
759 // res*= gcorp[rsigmay];
763 // if (padlength<1.1){
764 // snoise2 = gnoise2;
765 // res = ggg2[amp]*z+glandau2[amp]*angle2;
766 // if (ctype==0) res *= gcor02[rsigmay];
769 // res*= gcorp[rsigmay];
773 // snoise2 = gnoise3;
774 // res = ggg3[amp]*z+glandau3[amp]*angle2;
775 // if (ctype==0) res *= gcor02[rsigmay];
778 // res*= gcorp[rsigmay];
785 // res*=2.4; // overestimate error 2 times
789 // if (res<2*snoise2)
792 // seed->SetErrorY2(res);
800 Double_t AliTPCtrackerMI::ErrZ2(AliTPCseed* seed, const AliTPCclusterMI * cl){
803 // Use calibrated cluster error from OCDB
805 AliTPCClusterParam * clparam = AliTPCcalibDB::Instance()->GetClusterParam();
807 Float_t z = TMath::Abs(fkParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
808 Int_t ctype = cl->GetType();
809 Int_t type = (cl->GetRow()<63) ? 0: (cl->GetRow()>126) ? 1:2;
811 Double_t angle2 = seed->GetSnp()*seed->GetSnp();
812 angle2 = seed->GetTgl()*seed->GetTgl()*(1+angle2/(1-angle2));
813 Double_t angle = TMath::Sqrt(TMath::Abs(angle2));
814 Double_t errz2 = clparam->GetError0Par(1,type, z,angle);
816 errz2+=0.5; // edge cluster
819 seed->SetErrorZ2(errz2);
825 // //seed->SetErrorY2(0.1);
827 // //calculate look-up table at the beginning
828 // static Bool_t ginit = kFALSE;
829 // static Float_t gnoise1,gnoise2,gnoise3;
830 // static Float_t ggg1[10000];
831 // static Float_t ggg2[10000];
832 // static Float_t ggg3[10000];
833 // static Float_t glandau1[10000];
834 // static Float_t glandau2[10000];
835 // static Float_t glandau3[10000];
837 // static Float_t gcor01[1000];
838 // static Float_t gcor02[1000];
839 // static Float_t gcorp[1000];
843 // if (ginit==kFALSE){
844 // for (Int_t i=1;i<1000;i++){
845 // Float_t rsigma = float(i)/100.;
846 // gcor02[i] = TMath::Max(0.81 +TMath::Exp(6.8*(rsigma-1.2)),0.6);
847 // gcor01[i] = TMath::Max(0.72 +TMath::Exp(2.04*(rsigma-1.2)),0.6);
848 // gcorp[i] = TMath::Max(TMath::Power((rsigma+0.5),1.5),1.2);
852 // for (Int_t i=3;i<10000;i++){
856 // Float_t amp = float(i);
857 // Float_t padlength =0.75;
858 // gnoise1 = 0.0004/padlength;
859 // Float_t nel = 0.268*amp;
860 // Float_t nprim = 0.155*amp;
861 // ggg1[i] = fkParam->GetDiffT()*fkParam->GetDiffT()*(2+0.001*nel/(padlength*padlength))/nel;
862 // glandau1[i] = (2.+0.12*nprim)*0.5* (2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
863 // if (glandau1[i]>1) glandau1[i]=1;
864 // glandau1[i]*=padlength*padlength/12.;
868 // gnoise2 = 0.0004/padlength;
870 // nprim = 0.133*amp;
871 // ggg2[i] = fkParam->GetDiffT()*fkParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
872 // glandau2[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
873 // if (glandau2[i]>1) glandau2[i]=1;
874 // glandau2[i]*=padlength*padlength/12.;
879 // gnoise3 = 0.0004/padlength;
881 // nprim = 0.133*amp;
882 // ggg3[i] = fkParam->GetDiffT()*fkParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
883 // glandau3[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
884 // if (glandau3[i]>1) glandau3[i]=1;
885 // glandau3[i]*=padlength*padlength/12.;
893 // Int_t amp = int(TMath::Abs(cl->GetQ()));
895 // seed->SetErrorY2(1.);
899 // Float_t z = TMath::Abs(fkParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
900 // Int_t ctype = cl->GetType();
901 // Float_t padlength= GetPadPitchLength(seed->GetRow());
903 // Double_t angle2 = seed->GetSnp()*seed->GetSnp();
904 // // if (angle2<0.6) angle2 = 0.6;
905 // angle2 = seed->GetTgl()*seed->GetTgl()*(1+angle2/(1-angle2));
907 // //cluster "quality"
908 // Int_t rsigmaz = int(100.*cl->GetSigmaZ2()/(seed->GetCurrentSigmaZ2()));
911 // if (fSectors==fInnerSec){
912 // snoise2 = gnoise1;
913 // res = ggg1[amp]*z+glandau1[amp]*angle2;
914 // if (ctype==0) res *= gcor01[rsigmaz];
917 // res*= gcorp[rsigmaz];
921 // if (padlength<1.1){
922 // snoise2 = gnoise2;
923 // res = ggg2[amp]*z+glandau2[amp]*angle2;
924 // if (ctype==0) res *= gcor02[rsigmaz];
927 // res*= gcorp[rsigmaz];
931 // snoise2 = gnoise3;
932 // res = ggg3[amp]*z+glandau3[amp]*angle2;
933 // if (ctype==0) res *= gcor02[rsigmaz];
936 // res*= gcorp[rsigmaz];
945 // if ((ctype<0) &&<70){
950 // if (res<2*snoise2)
952 // if (res>3) res =3;
953 // seed->SetErrorZ2(res);
961 void AliTPCtrackerMI::RotateToLocal(AliTPCseed *seed)
963 //rotate to track "local coordinata
964 Float_t x = seed->GetX();
965 Float_t y = seed->GetY();
966 Float_t ymax = x*TMath::Tan(0.5*fSectors->GetAlpha());
969 seed->SetRelativeSector((seed->GetRelativeSector()+1) % fN);
970 if (!seed->Rotate(fSectors->GetAlpha()))
972 } else if (y <-ymax) {
973 seed->SetRelativeSector((seed->GetRelativeSector()-1+fN) % fN);
974 if (!seed->Rotate(-fSectors->GetAlpha()))
982 //_____________________________________________________________________________
983 Double_t AliTPCtrackerMI::F1old(Double_t x1,Double_t y1,
984 Double_t x2,Double_t y2,
985 Double_t x3,Double_t y3) const
987 //-----------------------------------------------------------------
988 // Initial approximation of the track curvature
989 //-----------------------------------------------------------------
990 Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
991 Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
992 (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
993 Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
994 (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
996 Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
997 if ( xr*xr+yr*yr<=0.00000000000001) return 100;
998 return -xr*yr/sqrt(xr*xr+yr*yr);
1003 //_____________________________________________________________________________
1004 Double_t AliTPCtrackerMI::F1(Double_t x1,Double_t y1,
1005 Double_t x2,Double_t y2,
1006 Double_t x3,Double_t y3) const
1008 //-----------------------------------------------------------------
1009 // Initial approximation of the track curvature
1010 //-----------------------------------------------------------------
1016 Double_t det = x3*y2-x2*y3;
1017 if (TMath::Abs(det)<1e-10){
1021 Double_t u = 0.5* (x2*(x2-x3)+y2*(y2-y3))/det;
1022 Double_t x0 = x3*0.5-y3*u;
1023 Double_t y0 = y3*0.5+x3*u;
1024 Double_t c2 = 1/TMath::Sqrt(x0*x0+y0*y0);
1030 Double_t AliTPCtrackerMI::F2(Double_t x1,Double_t y1,
1031 Double_t x2,Double_t y2,
1032 Double_t x3,Double_t y3) const
1034 //-----------------------------------------------------------------
1035 // Initial approximation of the track curvature
1036 //-----------------------------------------------------------------
1042 Double_t det = x3*y2-x2*y3;
1043 if (TMath::Abs(det)<1e-10) {
1047 Double_t u = 0.5* (x2*(x2-x3)+y2*(y2-y3))/det;
1048 Double_t x0 = x3*0.5-y3*u;
1049 Double_t y0 = y3*0.5+x3*u;
1050 Double_t c2 = 1/TMath::Sqrt(x0*x0+y0*y0);
1059 //_____________________________________________________________________________
1060 Double_t AliTPCtrackerMI::F2old(Double_t x1,Double_t y1,
1061 Double_t x2,Double_t y2,
1062 Double_t x3,Double_t y3) const
1064 //-----------------------------------------------------------------
1065 // Initial approximation of the track curvature times center of curvature
1066 //-----------------------------------------------------------------
1067 Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
1068 Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
1069 (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
1070 Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
1071 (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
1073 Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
1075 return -a/(d*y1-b)*xr/sqrt(xr*xr+yr*yr);
1078 //_____________________________________________________________________________
1079 Double_t AliTPCtrackerMI::F3(Double_t x1,Double_t y1,
1080 Double_t x2,Double_t y2,
1081 Double_t z1,Double_t z2) const
1083 //-----------------------------------------------------------------
1084 // Initial approximation of the tangent of the track dip angle
1085 //-----------------------------------------------------------------
1086 return (z1 - z2)/sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
1090 Double_t AliTPCtrackerMI::F3n(Double_t x1,Double_t y1,
1091 Double_t x2,Double_t y2,
1092 Double_t z1,Double_t z2, Double_t c) const
1094 //-----------------------------------------------------------------
1095 // Initial approximation of the tangent of the track dip angle
1096 //-----------------------------------------------------------------
1100 //angle1 = (z1-z2)*c/(TMath::ASin(c*x1-ni)-TMath::ASin(c*x2-ni));
1102 Double_t d = TMath::Sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
1103 if (TMath::Abs(d*c*0.5)>1) return 0;
1104 // Double_t angle2 = TMath::ASin(d*c*0.5);
1105 // Double_t angle2 = AliTPCFastMath::FastAsin(d*c*0.5);
1106 Double_t angle2 = (d*c*0.5>0.1)? TMath::ASin(d*c*0.5): AliTPCFastMath::FastAsin(d*c*0.5);
1108 angle2 = (z1-z2)*c/(angle2*2.);
1112 Bool_t AliTPCtrackerMI::GetProlongation(Double_t x1, Double_t x2, Double_t x[5], Double_t &y, Double_t &z) const
1113 {//-----------------------------------------------------------------
1114 // This function find proloncation of a track to a reference plane x=x2.
1115 //-----------------------------------------------------------------
1119 if (TMath::Abs(x[4]*x1 - x[2]) >= 0.999) {
1123 Double_t c1=x[4]*x1 - x[2], r1=TMath::Sqrt((1.-c1)*(1.+c1));
1124 Double_t c2=x[4]*x2 - x[2], r2=TMath::Sqrt((1.-c2)*(1.+c2));
1128 Double_t dy = dx*(c1+c2)/(r1+r2);
1131 Double_t delta = x[4]*dx*(c1+c2)/(c1*r2 + c2*r1);
1133 if (TMath::Abs(delta)>0.01){
1134 dz = x[3]*TMath::ASin(delta)/x[4];
1136 dz = x[3]*AliTPCFastMath::FastAsin(delta)/x[4];
1139 //dz = x[3]*AliTPCFastMath::FastAsin(delta)/x[4];
1147 Int_t AliTPCtrackerMI::LoadClusters (TTree *const tree)
1152 return LoadClusters();
1156 Int_t AliTPCtrackerMI::LoadClusters(const TObjArray *arr)
1159 // load clusters to the memory
1160 AliTPCClustersRow *clrow = new AliTPCClustersRow("AliTPCclusterMI");
1161 Int_t lower = arr->LowerBound();
1162 Int_t entries = arr->GetEntriesFast();
1164 for (Int_t i=lower; i<entries; i++) {
1165 clrow = (AliTPCClustersRow*) arr->At(i);
1166 if(!clrow) continue;
1167 if(!clrow->GetArray()) continue;
1171 fkParam->AdjustSectorRow(clrow->GetID(),sec,row);
1173 for (Int_t icl=0; icl<clrow->GetArray()->GetEntriesFast(); icl++){
1174 Transform((AliTPCclusterMI*)(clrow->GetArray()->At(icl)));
1177 if (clrow->GetArray()->GetEntriesFast()<=0) continue;
1178 AliTPCtrackerRow * tpcrow=0;
1181 tpcrow = &(fInnerSec[sec%fkNIS][row]);
1185 tpcrow = &(fOuterSec[(sec-fkNIS*2)%fkNOS][row]);
1186 left = (sec-fkNIS*2)/fkNOS;
1189 tpcrow->SetN1(clrow->GetArray()->GetEntriesFast());
1190 tpcrow->SetClusters1(new AliTPCclusterMI[tpcrow->GetN1()]);
1191 for (Int_t j=0;j<tpcrow->GetN1();++j)
1192 tpcrow->SetCluster1(j, *(AliTPCclusterMI*)(clrow->GetArray()->At(j)));
1195 tpcrow->SetN2(clrow->GetArray()->GetEntriesFast());
1196 tpcrow->SetClusters2(new AliTPCclusterMI[tpcrow->GetN2()]);
1197 for (Int_t j=0;j<tpcrow->GetN2();++j)
1198 tpcrow->SetCluster2(j,*(AliTPCclusterMI*)(clrow->GetArray()->At(j)));
1200 clrow->GetArray()->Clear("C");
1209 Int_t AliTPCtrackerMI::LoadClusters(const TClonesArray *arr)
1212 // load clusters to the memory from one
1215 AliTPCclusterMI *clust=0;
1216 Int_t count[72][96] = { {0} , {0} };
1218 // loop over clusters
1219 for (Int_t icl=0; icl<arr->GetEntriesFast(); icl++) {
1220 clust = (AliTPCclusterMI*)arr->At(icl);
1221 if(!clust) continue;
1222 //printf("cluster: det %d, row %d \n", clust->GetDetector(),clust->GetRow());
1224 // transform clusters
1227 // count clusters per pad row
1228 count[clust->GetDetector()][clust->GetRow()]++;
1231 // insert clusters to sectors
1232 for (Int_t icl=0; icl<arr->GetEntriesFast(); icl++) {
1233 clust = (AliTPCclusterMI*)arr->At(icl);
1234 if(!clust) continue;
1236 Int_t sec = clust->GetDetector();
1237 Int_t row = clust->GetRow();
1239 // filter overlapping pad rows needed by HLT
1240 if(sec<fkNIS*2) { //IROCs
1241 if(row == 30) continue;
1244 if(row == 27 || row == 76) continue;
1250 fInnerSec[sec%fkNIS].InsertCluster(clust, count[sec][row], fkParam);
1253 left = (sec-fkNIS*2)/fkNOS;
1254 fOuterSec[(sec-fkNIS*2)%fkNOS].InsertCluster(clust, count[sec][row], fkParam);
1258 // Load functions must be called behind LoadCluster(TClonesArray*)
1260 //LoadOuterSectors();
1261 //LoadInnerSectors();
1267 Int_t AliTPCtrackerMI::LoadClusters()
1270 // load clusters to the memory
1271 AliTPCClustersRow *clrow= new AliTPCClustersRow("AliTPCclusterMI");
1273 // TTree * tree = fClustersArray.GetTree();
1275 TTree * tree = fInput;
1276 TBranch * br = tree->GetBranch("Segment");
1277 br->SetAddress(&clrow);
1279 // Conversion of pad, row coordinates in local tracking coords.
1280 // Could be skipped here; is already done in clusterfinder
1282 Int_t j=Int_t(tree->GetEntries());
1283 for (Int_t i=0; i<j; i++) {
1287 fkParam->AdjustSectorRow(clrow->GetID(),sec,row);
1288 for (Int_t icl=0; icl<clrow->GetArray()->GetEntriesFast(); icl++){
1289 Transform((AliTPCclusterMI*)(clrow->GetArray()->At(icl)));
1292 AliTPCtrackerRow * tpcrow=0;
1295 tpcrow = &(fInnerSec[sec%fkNIS][row]);
1299 tpcrow = &(fOuterSec[(sec-fkNIS*2)%fkNOS][row]);
1300 left = (sec-fkNIS*2)/fkNOS;
1303 tpcrow->SetN1(clrow->GetArray()->GetEntriesFast());
1304 tpcrow->SetClusters1(new AliTPCclusterMI[tpcrow->GetN1()]);
1305 for (Int_t k=0;k<tpcrow->GetN1();++k)
1306 tpcrow->SetCluster1(k, *(AliTPCclusterMI*)(clrow->GetArray()->At(k)));
1309 tpcrow->SetN2(clrow->GetArray()->GetEntriesFast());
1310 tpcrow->SetClusters2(new AliTPCclusterMI[tpcrow->GetN2()]);
1311 for (Int_t k=0;k<tpcrow->GetN2();++k)
1312 tpcrow->SetCluster2(k,*(AliTPCclusterMI*)(clrow->GetArray()->At(k)));
1323 void AliTPCtrackerMI::UnloadClusters()
1326 // unload clusters from the memory
1328 Int_t nrows = fOuterSec->GetNRows();
1329 for (Int_t sec = 0;sec<fkNOS;sec++)
1330 for (Int_t row = 0;row<nrows;row++){
1331 AliTPCtrackerRow* tpcrow = &(fOuterSec[sec%fkNOS][row]);
1333 // if (tpcrow->fClusters1) delete []tpcrow->fClusters1;
1334 // if (tpcrow->fClusters2) delete []tpcrow->fClusters2;
1336 tpcrow->ResetClusters();
1339 nrows = fInnerSec->GetNRows();
1340 for (Int_t sec = 0;sec<fkNIS;sec++)
1341 for (Int_t row = 0;row<nrows;row++){
1342 AliTPCtrackerRow* tpcrow = &(fInnerSec[sec%fkNIS][row]);
1344 // if (tpcrow->fClusters1) delete []tpcrow->fClusters1;
1345 //if (tpcrow->fClusters2) delete []tpcrow->fClusters2;
1347 tpcrow->ResetClusters();
1353 void AliTPCtrackerMI::FillClusterArray(TObjArray* array) const{
1355 // Filling cluster to the array - For visualization purposes
1358 nrows = fOuterSec->GetNRows();
1359 for (Int_t sec = 0;sec<fkNOS;sec++)
1360 for (Int_t row = 0;row<nrows;row++){
1361 AliTPCtrackerRow* tpcrow = &(fOuterSec[sec%fkNOS][row]);
1362 if (!tpcrow) continue;
1363 for (Int_t icl = 0;icl<tpcrow->GetN();icl++){
1364 array->AddLast((TObject*)((*tpcrow)[icl]));
1367 nrows = fInnerSec->GetNRows();
1368 for (Int_t sec = 0;sec<fkNIS;sec++)
1369 for (Int_t row = 0;row<nrows;row++){
1370 AliTPCtrackerRow* tpcrow = &(fInnerSec[sec%fkNIS][row]);
1371 if (!tpcrow) continue;
1372 for (Int_t icl = 0;icl<tpcrow->GetN();icl++){
1373 array->AddLast((TObject*)(*tpcrow)[icl]);
1379 void AliTPCtrackerMI::Transform(AliTPCclusterMI * cluster){
1383 AliTPCcalibDB * calibDB = AliTPCcalibDB::Instance();
1384 AliTPCTransform *transform = calibDB->GetTransform() ;
1386 AliFatal("Tranformations not in calibDB");
1389 transform->SetCurrentRecoParam((AliTPCRecoParam*)AliTPCReconstructor::GetRecoParam());
1390 Double_t x[3]={cluster->GetRow(),cluster->GetPad(),cluster->GetTimeBin()};
1391 Int_t i[1]={cluster->GetDetector()};
1392 transform->Transform(x,i,0,1);
1393 // if (cluster->GetDetector()%36>17){
1398 // in debug mode check the transformation
1400 if (AliTPCReconstructor::StreamLevel()>2) {
1402 cluster->GetGlobalXYZ(gx);
1403 Int_t event = (fEvent==NULL)? 0: fEvent->GetEventNumberInFile();
1404 TTreeSRedirector &cstream = *fDebugStreamer;
1405 cstream<<"Transform"<<
1416 cluster->SetX(x[0]);
1417 cluster->SetY(x[1]);
1418 cluster->SetZ(x[2]);
1423 //if (!fkParam->IsGeoRead()) fkParam->ReadGeoMatrices();
1424 if (AliTPCReconstructor::GetRecoParam()->GetUseSectorAlignment() && (!calibDB->HasAlignmentOCDB())){
1425 TGeoHMatrix *mat = fkParam->GetClusterMatrix(cluster->GetDetector());
1427 Double_t pos[3]= {cluster->GetX(),cluster->GetY(),cluster->GetZ()};
1428 Double_t posC[3]={cluster->GetX(),cluster->GetY(),cluster->GetZ()};
1429 if (mat) mat->LocalToMaster(pos,posC);
1431 // chack Loading of Geo matrices from GeoManager - TEMPORARY FIX
1433 cluster->SetX(posC[0]);
1434 cluster->SetY(posC[1]);
1435 cluster->SetZ(posC[2]);
1439 //_____________________________________________________________________________
1440 Int_t AliTPCtrackerMI::LoadOuterSectors() {
1441 //-----------------------------------------------------------------
1442 // This function fills outer TPC sectors with clusters.
1443 //-----------------------------------------------------------------
1444 Int_t nrows = fOuterSec->GetNRows();
1446 for (Int_t sec = 0;sec<fkNOS;sec++)
1447 for (Int_t row = 0;row<nrows;row++){
1448 AliTPCtrackerRow* tpcrow = &(fOuterSec[sec%fkNOS][row]);
1449 Int_t sec2 = sec+2*fkNIS;
1451 Int_t ncl = tpcrow->GetN1();
1453 AliTPCclusterMI *c= (tpcrow->GetCluster1(ncl));
1454 index=(((sec2<<8)+row)<<16)+ncl;
1455 tpcrow->InsertCluster(c,index);
1458 ncl = tpcrow->GetN2();
1460 AliTPCclusterMI *c= (tpcrow->GetCluster2(ncl));
1461 index=((((sec2+fkNOS)<<8)+row)<<16)+ncl;
1462 tpcrow->InsertCluster(c,index);
1465 // write indexes for fast acces
1467 for (Int_t i=0;i<510;i++)
1468 tpcrow->SetFastCluster(i,-1);
1469 for (Int_t i=0;i<tpcrow->GetN();i++){
1470 Int_t zi = Int_t((*tpcrow)[i]->GetZ()+255.);
1471 tpcrow->SetFastCluster(zi,i); // write index
1474 for (Int_t i=0;i<510;i++){
1475 if (tpcrow->GetFastCluster(i)<0)
1476 tpcrow->SetFastCluster(i,last);
1478 last = tpcrow->GetFastCluster(i);
1487 //_____________________________________________________________________________
1488 Int_t AliTPCtrackerMI::LoadInnerSectors() {
1489 //-----------------------------------------------------------------
1490 // This function fills inner TPC sectors with clusters.
1491 //-----------------------------------------------------------------
1492 Int_t nrows = fInnerSec->GetNRows();
1494 for (Int_t sec = 0;sec<fkNIS;sec++)
1495 for (Int_t row = 0;row<nrows;row++){
1496 AliTPCtrackerRow* tpcrow = &(fInnerSec[sec%fkNIS][row]);
1499 Int_t ncl = tpcrow->GetN1();
1501 AliTPCclusterMI *c= (tpcrow->GetCluster1(ncl));
1502 index=(((sec<<8)+row)<<16)+ncl;
1503 tpcrow->InsertCluster(c,index);
1506 ncl = tpcrow->GetN2();
1508 AliTPCclusterMI *c= (tpcrow->GetCluster2(ncl));
1509 index=((((sec+fkNIS)<<8)+row)<<16)+ncl;
1510 tpcrow->InsertCluster(c,index);
1513 // write indexes for fast acces
1515 for (Int_t i=0;i<510;i++)
1516 tpcrow->SetFastCluster(i,-1);
1517 for (Int_t i=0;i<tpcrow->GetN();i++){
1518 Int_t zi = Int_t((*tpcrow)[i]->GetZ()+255.);
1519 tpcrow->SetFastCluster(zi,i); // write index
1522 for (Int_t i=0;i<510;i++){
1523 if (tpcrow->GetFastCluster(i)<0)
1524 tpcrow->SetFastCluster(i,last);
1526 last = tpcrow->GetFastCluster(i);
1538 //_________________________________________________________________________
1539 AliTPCclusterMI *AliTPCtrackerMI::GetClusterMI(Int_t index) const {
1540 //--------------------------------------------------------------------
1541 // Return pointer to a given cluster
1542 //--------------------------------------------------------------------
1543 if (index<0) return 0; // no cluster
1544 Int_t sec=(index&0xff000000)>>24;
1545 Int_t row=(index&0x00ff0000)>>16;
1546 Int_t ncl=(index&0x00007fff)>>00;
1548 const AliTPCtrackerRow * tpcrow=0;
1549 AliTPCclusterMI * clrow =0;
1551 if (sec<0 || sec>=fkNIS*4) {
1552 AliWarning(Form("Wrong sector %d",sec));
1557 AliTPCtrackerSector& tracksec = fInnerSec[sec%fkNIS];
1558 if (tracksec.GetNRows()<=row) return 0;
1559 tpcrow = &(tracksec[row]);
1560 if (tpcrow==0) return 0;
1563 if (tpcrow->GetN1()<=ncl) return 0;
1564 clrow = tpcrow->GetClusters1();
1567 if (tpcrow->GetN2()<=ncl) return 0;
1568 clrow = tpcrow->GetClusters2();
1572 AliTPCtrackerSector& tracksec = fOuterSec[(sec-fkNIS*2)%fkNOS];
1573 if (tracksec.GetNRows()<=row) return 0;
1574 tpcrow = &(tracksec[row]);
1575 if (tpcrow==0) return 0;
1577 if (sec-2*fkNIS<fkNOS) {
1578 if (tpcrow->GetN1()<=ncl) return 0;
1579 clrow = tpcrow->GetClusters1();
1582 if (tpcrow->GetN2()<=ncl) return 0;
1583 clrow = tpcrow->GetClusters2();
1587 return &(clrow[ncl]);
1593 Int_t AliTPCtrackerMI::FollowToNext(AliTPCseed& t, Int_t nr) {
1594 //-----------------------------------------------------------------
1595 // This function tries to find a track prolongation to next pad row
1596 //-----------------------------------------------------------------
1598 Double_t x= GetXrow(nr), ymax=GetMaxY(nr);
1601 AliTPCclusterMI *cl=0;
1602 Int_t tpcindex= t.GetClusterIndex2(nr);
1604 // update current shape info every 5 pad-row
1605 // if ( (nr%5==0) || t.GetNumberOfClusters()<2 || (t.fCurrentSigmaY2<0.0001) ){
1609 if (fIteration>0 && tpcindex>=-1){ //if we have already clusters
1611 if (tpcindex==-1) return 0; //track in dead zone
1612 if (tpcindex >= 0){ //
1613 cl = t.GetClusterPointer(nr);
1614 //if (cl==0) cl = GetClusterMI(tpcindex);
1615 if (!cl) cl = GetClusterMI(tpcindex);
1616 t.SetCurrentClusterIndex1(tpcindex);
1619 Int_t relativesector = ((tpcindex&0xff000000)>>24)%18; // if previously accepted cluster in different sector
1620 Float_t angle = relativesector*fSectors->GetAlpha()+fSectors->GetAlphaShift();
1622 if (angle<-TMath::Pi()) angle += 2*TMath::Pi();
1623 if (angle>=TMath::Pi()) angle -= 2*TMath::Pi();
1625 if (TMath::Abs(angle-t.GetAlpha())>0.001){
1626 Double_t rotation = angle-t.GetAlpha();
1627 t.SetRelativeSector(relativesector);
1628 if (!t.Rotate(rotation)) return 0;
1630 if (!t.PropagateTo(x)) return 0;
1632 t.SetCurrentCluster(cl);
1634 Int_t accept = AcceptCluster(&t,t.GetCurrentCluster());
1635 if ((tpcindex&0x8000)==0) accept =0;
1637 //if founded cluster is acceptible
1638 if (cl->IsUsed(11)) { // id cluster is shared inrease uncertainty
1639 t.SetErrorY2(t.GetErrorY2()+0.03);
1640 t.SetErrorZ2(t.GetErrorZ2()+0.03);
1641 t.SetErrorY2(t.GetErrorY2()*3);
1642 t.SetErrorZ2(t.GetErrorZ2()*3);
1644 t.SetNFoundable(t.GetNFoundable()+1);
1645 UpdateTrack(&t,accept);
1648 else { // Remove old cluster from track
1649 t.SetClusterIndex(nr, -3);
1650 t.SetClusterPointer(nr, 0);
1654 if (TMath::Abs(t.GetSnp())>AliTPCReconstructor::GetMaxSnpTracker()) return 0; // cut on angle
1655 if (fIteration>1 && IsFindable(t)){
1656 // not look for new cluster during refitting
1657 t.SetNFoundable(t.GetNFoundable()+1);
1662 // if (TMath::Abs(t.GetSnp())>0.95 || TMath::Abs(x*t.GetC()-t.GetEta())>0.95) return 0;// patch 28 fev 06
1663 if (!t.PropagateTo(x)) {
1664 if (fIteration==0) t.SetRemoval(10);
1667 Double_t y = t.GetY();
1668 if (TMath::Abs(y)>ymax){
1670 t.SetRelativeSector((t.GetRelativeSector()+1) % fN);
1671 if (!t.Rotate(fSectors->GetAlpha()))
1673 } else if (y <-ymax) {
1674 t.SetRelativeSector((t.GetRelativeSector()-1+fN) % fN);
1675 if (!t.Rotate(-fSectors->GetAlpha()))
1678 if (!t.PropagateTo(x)) {
1679 if (fIteration==0) t.SetRemoval(10);
1685 Double_t z=t.GetZ();
1688 if (!IsActive(t.GetRelativeSector(),nr)) {
1690 t.SetClusterIndex2(nr,-1);
1693 //AliInfo(Form("A - Sector%d phi %f - alpha %f", t.fRelativeSector,y/x, t.GetAlpha()));
1694 Bool_t isActive = IsActive(t.GetRelativeSector(),nr);
1695 Bool_t isActive2 = (nr>=fInnerSec->GetNRows()) ? fOuterSec[t.GetRelativeSector()][nr-fInnerSec->GetNRows()].GetN()>0:fInnerSec[t.GetRelativeSector()][nr].GetN()>0;
1697 if (!isActive || !isActive2) return 0;
1699 const AliTPCtrackerRow &krow=GetRow(t.GetRelativeSector(),nr);
1700 if ( (t.GetSigmaY2()<0) || t.GetSigmaZ2()<0) return 0;
1702 Double_t roadz = 1.;
1704 if (TMath::Abs(TMath::Abs(y)-ymax)<krow.GetDeadZone()){
1706 t.SetClusterIndex2(nr,-1);
1712 // if (TMath::Abs(z)<(AliTPCReconstructor::GetCtgRange()*x+10) && TMath::Abs(z)<fkParam->GetZLength(0) && (TMath::Abs(t.GetSnp())<AliTPCReconstructor::GetMaxSnpTracker()))
1713 t.SetNFoundable(t.GetNFoundable()+1);
1719 // cl = krow.FindNearest2(y+10.,z,roady,roadz,index);
1720 cl = krow.FindNearest2(y,z,roady,roadz,index);
1721 if (cl) t.SetCurrentClusterIndex1(krow.GetIndex(index));
1724 t.SetCurrentCluster(cl);
1726 if (fIteration==2&&cl->IsUsed(10)) return 0;
1727 Int_t accept = AcceptCluster(&t,t.GetCurrentCluster());
1728 if (fIteration==2&&cl->IsUsed(11)) {
1729 t.SetErrorY2(t.GetErrorY2()+0.03);
1730 t.SetErrorZ2(t.GetErrorZ2()+0.03);
1731 t.SetErrorY2(t.GetErrorY2()*3);
1732 t.SetErrorZ2(t.GetErrorZ2()*3);
1735 if (t.fCurrentCluster->IsUsed(10)){
1740 if (t.fNShared>0.7*t.GetNumberOfClusters()) {
1746 if (accept<3) UpdateTrack(&t,accept);
1749 if ( fIteration==0 && t.GetNFoundable()*0.5 > t.GetNumberOfClusters()) t.SetRemoval(10);
1757 //_________________________________________________________________________
1758 Bool_t AliTPCtrackerMI::GetTrackPoint(Int_t index, AliTrackPoint &p ) const
1760 // Get track space point by index
1761 // return false in case the cluster doesn't exist
1762 AliTPCclusterMI *cl = GetClusterMI(index);
1763 if (!cl) return kFALSE;
1764 Int_t sector = (index&0xff000000)>>24;
1765 // Int_t row = (index&0x00ff0000)>>16;
1767 // xyz[0] = fkParam->GetPadRowRadii(sector,row);
1768 xyz[0] = cl->GetX();
1769 xyz[1] = cl->GetY();
1770 xyz[2] = cl->GetZ();
1772 fkParam->AdjustCosSin(sector,cos,sin);
1773 Float_t x = cos*xyz[0]-sin*xyz[1];
1774 Float_t y = cos*xyz[1]+sin*xyz[0];
1776 Float_t sigmaY2 = 0.027*cl->GetSigmaY2();
1777 if (sector < fkParam->GetNInnerSector()) sigmaY2 *= 2.07;
1778 Float_t sigmaZ2 = 0.066*cl->GetSigmaZ2();
1779 if (sector < fkParam->GetNInnerSector()) sigmaZ2 *= 1.77;
1780 cov[0] = sin*sin*sigmaY2;
1781 cov[1] = -sin*cos*sigmaY2;
1783 cov[3] = cos*cos*sigmaY2;
1786 p.SetXYZ(x,y,xyz[2],cov);
1787 AliGeomManager::ELayerID iLayer;
1789 if (sector < fkParam->GetNInnerSector()) {
1790 iLayer = AliGeomManager::kTPC1;
1794 iLayer = AliGeomManager::kTPC2;
1795 idet = sector - fkParam->GetNInnerSector();
1797 UShort_t volid = AliGeomManager::LayerToVolUID(iLayer,idet);
1798 p.SetVolumeID(volid);
1804 Int_t AliTPCtrackerMI::UpdateClusters(AliTPCseed& t, Int_t nr) {
1805 //-----------------------------------------------------------------
1806 // This function tries to find a track prolongation to next pad row
1807 //-----------------------------------------------------------------
1808 t.SetCurrentCluster(0);
1809 t.SetCurrentClusterIndex1(-3);
1811 Double_t xt=t.GetX();
1812 Int_t row = GetRowNumber(xt)-1;
1813 Double_t ymax= GetMaxY(nr);
1815 if (row < nr) return 1; // don't prolongate if not information until now -
1816 // if (TMath::Abs(t.GetSnp())>0.9 && t.GetNumberOfClusters()>40. && fIteration!=2) {
1818 // return 0; // not prolongate strongly inclined tracks
1820 // if (TMath::Abs(t.GetSnp())>0.95) {
1822 // return 0; // not prolongate strongly inclined tracks
1823 // }// patch 28 fev 06
1825 Double_t x= GetXrow(nr);
1827 //t.PropagateTo(x+0.02);
1828 //t.PropagateTo(x+0.01);
1829 if (!t.PropagateTo(x)){
1836 if (TMath::Abs(y)>ymax){
1838 t.SetRelativeSector((t.GetRelativeSector()+1) % fN);
1839 if (!t.Rotate(fSectors->GetAlpha()))
1841 } else if (y <-ymax) {
1842 t.SetRelativeSector((t.GetRelativeSector()-1+fN) % fN);
1843 if (!t.Rotate(-fSectors->GetAlpha()))
1846 // if (!t.PropagateTo(x)){
1853 if (TMath::Abs(t.GetSnp())>AliTPCReconstructor::GetMaxSnpTracker()) return 0;
1855 if (!IsActive(t.GetRelativeSector(),nr)) {
1857 t.SetClusterIndex2(nr,-1);
1860 //AliInfo(Form("A - Sector%d phi %f - alpha %f", t.fRelativeSector,y/x, t.GetAlpha()));
1862 AliTPCtrackerRow &krow=GetRow(t.GetRelativeSector(),nr);
1864 if (TMath::Abs(TMath::Abs(y)-ymax)<krow.GetDeadZone()){
1866 t.SetClusterIndex2(nr,-1);
1872 // if (TMath::Abs(t.GetZ())<(AliTPCReconstructor::GetCtgRange()*t.GetX()+10) && (TMath::Abs(t.GetSnp())<AliTPCReconstructor::GetMaxSnpTracker()))
1873 if (IsFindable(t)) t.SetNFoundable(t.GetNFoundable()+1);
1879 if ( (nr%2==0) || t.GetNumberOfClusters()<2){
1880 // t.fCurrentSigmaY = GetSigmaY(&t);
1881 //t.fCurrentSigmaZ = GetSigmaZ(&t);
1885 AliTPCclusterMI *cl=0;
1888 Double_t roady = 1.;
1889 Double_t roadz = 1.;
1893 index = t.GetClusterIndex2(nr);
1894 if ( (index >= 0) && (index&0x8000)==0){
1895 cl = t.GetClusterPointer(nr);
1896 if ( (cl==0) && (index >= 0)) cl = GetClusterMI(index);
1897 t.SetCurrentClusterIndex1(index);
1899 t.SetCurrentCluster(cl);
1905 // if (index<0) return 0;
1906 UInt_t uindex = TMath::Abs(index);
1909 //cl = krow.FindNearest2(y+10,z,roady,roadz,uindex);
1910 cl = krow.FindNearest2(y,z,roady,roadz,uindex);
1913 if (cl) t.SetCurrentClusterIndex1(krow.GetIndex(uindex));
1914 t.SetCurrentCluster(cl);
1920 Int_t AliTPCtrackerMI::FollowToNextCluster(AliTPCseed & t, Int_t nr) {
1921 //-----------------------------------------------------------------
1922 // This function tries to find a track prolongation to next pad row
1923 //-----------------------------------------------------------------
1925 //update error according neighborhoud
1927 if (t.GetCurrentCluster()) {
1929 Int_t accept = AcceptCluster(&t,t.GetCurrentCluster());
1931 if (t.GetCurrentCluster()->IsUsed(10)){
1936 t.SetNShared(t.GetNShared()+1);
1937 if (t.GetNShared()>0.7*t.GetNumberOfClusters()) {
1942 if (fIteration>0) accept = 0;
1943 if (accept<3) UpdateTrack(&t,accept);
1947 if ( t.GetNumberOfClusters()>18 && ( (t.GetSigmaY2()+t.GetSigmaZ2())>0.16)) t.SetRemoval(10);
1948 if ( t.GetNumberOfClusters()>18 && t.GetChi2()/t.GetNumberOfClusters()>6 ) t.SetRemoval(10);
1950 if (( (t.GetNFoundable()*0.5 > t.GetNumberOfClusters()) || t.GetNoCluster()>15)) t.SetRemoval(10);
1958 //_____________________________________________________________________________
1959 Int_t AliTPCtrackerMI::FollowProlongation(AliTPCseed& t, Int_t rf, Int_t step, Bool_t fromSeeds) {
1960 //-----------------------------------------------------------------
1961 // This function tries to find a track prolongation.
1962 //-----------------------------------------------------------------
1963 Double_t xt=t.GetX();
1965 Double_t alpha=t.GetAlpha();
1966 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
1967 if (alpha < 0. ) alpha += 2.*TMath::Pi();
1969 t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN);
1971 Int_t first = GetRowNumber(xt);
1976 for (Int_t nr= first; nr>=rf; nr-=step) {
1978 if (t.GetKinkIndexes()[0]>0){
1979 for (Int_t i=0;i<3;i++){
1980 Int_t index = t.GetKinkIndexes()[i];
1981 if (index==0) break;
1982 if (index<0) continue;
1984 AliKink * kink = (AliKink*)fEvent->GetKink(index-1);
1986 printf("PROBLEM\n");
1989 Int_t kinkrow = kink->GetTPCRow0()+2+Int_t(0.5/(0.05+kink->GetAngle(2)));
1991 AliExternalTrackParam paramd(t);
1992 kink->SetDaughter(paramd);
1993 kink->SetStatus(2,5);
2000 if (nr==80) t.UpdateReference();
2001 if (nr<fInnerSec->GetNRows())
2002 fSectors = fInnerSec;
2004 fSectors = fOuterSec;
2005 if (FollowToNext(t,nr)==0)
2018 Int_t AliTPCtrackerMI::FollowBackProlongation(AliTPCseed& t, Int_t rf, Bool_t fromSeeds) {
2019 //-----------------------------------------------------------------
2020 // This function tries to find a track prolongation.
2021 //-----------------------------------------------------------------
2023 Double_t xt=t.GetX();
2024 Double_t alpha=t.GetAlpha();
2025 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
2026 if (alpha < 0. ) alpha += 2.*TMath::Pi();
2027 t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN);
2029 Int_t first = t.GetFirstPoint();
2030 Int_t ri = GetRowNumber(xt);
2034 if (first<ri) first = ri;
2036 if (first<0) first=0;
2037 for (Int_t nr=first; nr<=rf; nr++) {
2038 // if ( (TMath::Abs(t.GetSnp())>0.95)) break;//patch 28 fev 06
2039 if (t.GetKinkIndexes()[0]<0){
2040 for (Int_t i=0;i<3;i++){
2041 Int_t index = t.GetKinkIndexes()[i];
2042 if (index==0) break;
2043 if (index>0) continue;
2044 index = TMath::Abs(index);
2045 AliKink * kink = (AliKink*)fEvent->GetKink(index-1);
2047 printf("PROBLEM\n");
2050 Int_t kinkrow = kink->GetTPCRow0()-2-Int_t(0.5/(0.05+kink->GetAngle(2)));
2052 AliExternalTrackParam paramm(t);
2053 kink->SetMother(paramm);
2054 kink->SetStatus(2,1);
2061 if (nr<fInnerSec->GetNRows())
2062 fSectors = fInnerSec;
2064 fSectors = fOuterSec;
2075 Float_t AliTPCtrackerMI::OverlapFactor(AliTPCseed * s1, AliTPCseed * s2, Int_t &sum1, Int_t & sum2)
2077 // overlapping factor
2083 Float_t dz2 =(s1->GetZ() - s2->GetZ());
2086 Float_t dy2 =TMath::Abs((s1->GetY() - s2->GetY()));
2088 Float_t distance = TMath::Sqrt(dz2+dy2);
2089 if (distance>4.) return 0; // if there are far away - not overlap - to reduce combinatorics
2092 Int_t firstpoint = TMath::Min(s1->GetFirstPoint(),s2->GetFirstPoint());
2093 Int_t lastpoint = TMath::Max(s1->GetLastPoint(),s2->GetLastPoint());
2098 if (firstpoint>lastpoint) {
2099 firstpoint =lastpoint;
2104 for (Int_t i=firstpoint-1;i<lastpoint+1;i++){
2105 if (s1->GetClusterIndex2(i)>0) sum1++;
2106 if (s2->GetClusterIndex2(i)>0) sum2++;
2107 if (s1->GetClusterIndex2(i)==s2->GetClusterIndex2(i) && s1->GetClusterIndex2(i)>0) {
2111 if (sum<5) return 0;
2113 Float_t summin = TMath::Min(sum1+1,sum2+1);
2114 Float_t ratio = (sum+1)/Float_t(summin);
2118 void AliTPCtrackerMI::SignShared(AliTPCseed * s1, AliTPCseed * s2)
2122 Float_t thetaCut = 0.2;//+10.*TMath::Sqrt(s1->GetSigmaTglZ()+ s2->GetSigmaTglZ());
2123 if (TMath::Abs(s1->GetTgl()-s2->GetTgl())>thetaCut) return;
2124 Float_t minCl = TMath::Min(s1->GetNumberOfClusters(),s2->GetNumberOfClusters());
2125 Int_t cutN0 = TMath::Max(5,TMath::Nint(0.1*minCl));
2130 //Int_t firstpoint = TMath::Max(s1->GetFirstPoint(),s2->GetFirstPoint());
2131 //Int_t lastpoint = TMath::Min(s1->GetLastPoint(),s2->GetLastPoint());
2132 Int_t firstpoint = 0;
2133 Int_t lastpoint = 160;
2135 // if (firstpoint>=lastpoint-5) return;;
2137 for (Int_t i=firstpoint;i<lastpoint;i++){
2138 // if ( (s1->GetClusterIndex2(i)&0xFFFF8FFF)==(s2->GetClusterIndex2(i)&0xFFFF8FFF) && s1->GetClusterIndex2(i)>0) {
2139 if ( (s1->GetClusterIndex2(i))==(s2->GetClusterIndex2(i)) && s1->GetClusterIndex2(i)>0) {
2143 if (sumshared>cutN0){
2146 for (Int_t i=firstpoint;i<lastpoint;i++){
2147 // if ( (s1->GetClusterIndex2(i)&0xFFFF8FFF)==(s2->GetClusterIndex2(i)&0xFFFF8FFF) && s1->GetClusterIndex2(i)>0) {
2148 if ( (s1->GetClusterIndex2(i))==(s2->GetClusterIndex2(i)) && s1->GetClusterIndex2(i)>0) {
2149 AliTPCTrackerPoint *p1 = s1->GetTrackPoint(i);
2150 AliTPCTrackerPoint *p2 = s2->GetTrackPoint(i);;
2151 if (s1->IsActive()&&s2->IsActive()){
2152 p1->SetShared(kTRUE);
2153 p2->SetShared(kTRUE);
2159 if (sumshared>cutN0){
2160 for (Int_t i=0;i<4;i++){
2161 if (s1->GetOverlapLabel(3*i)==0){
2162 s1->SetOverlapLabel(3*i, s2->GetLabel());
2163 s1->SetOverlapLabel(3*i+1,sumshared);
2164 s1->SetOverlapLabel(3*i+2,s2->GetUniqueID());
2168 for (Int_t i=0;i<4;i++){
2169 if (s2->GetOverlapLabel(3*i)==0){
2170 s2->SetOverlapLabel(3*i, s1->GetLabel());
2171 s2->SetOverlapLabel(3*i+1,sumshared);
2172 s2->SetOverlapLabel(3*i+2,s1->GetUniqueID());
2179 void AliTPCtrackerMI::SignShared(TObjArray * arr)
2182 //sort trackss according sectors
2184 for (Int_t i=0; i<arr->GetEntriesFast(); i++) {
2185 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2187 //if (pt) RotateToLocal(pt);
2191 arr->Sort(); // sorting according relative sectors
2192 arr->Expand(arr->GetEntries());
2195 Int_t nseed=arr->GetEntriesFast();
2196 for (Int_t i=0; i<nseed; i++) {
2197 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2199 for (Int_t j=0;j<12;j++){
2200 pt->SetOverlapLabel(j,0);
2203 for (Int_t i=0; i<nseed; i++) {
2204 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2206 if (pt->GetRemoval()>10) continue;
2207 for (Int_t j=i+1; j<nseed; j++){
2208 AliTPCseed *pt2=(AliTPCseed*)arr->UncheckedAt(j);
2209 if (TMath::Abs(pt->GetRelativeSector()-pt2->GetRelativeSector())>1) continue;
2211 if (pt2->GetRemoval()<=10) {
2212 //if ( TMath::Abs(pt->GetRelativeSector()-pt2->GetRelativeSector())>0) break;
2220 void AliTPCtrackerMI::SortTracks(TObjArray * arr, Int_t mode) const
2223 //sort tracks in array according mode criteria
2224 Int_t nseed = arr->GetEntriesFast();
2225 for (Int_t i=0; i<nseed; i++) {
2226 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2237 void AliTPCtrackerMI::RemoveUsed2(TObjArray * arr, Float_t factor1, Float_t factor2, Int_t minimal)
2240 // Loop over all tracks and remove overlaped tracks (with lower quality)
2242 // 1. Unsign clusters
2243 // 2. Sort tracks according quality
2244 // Quality is defined by the number of cluster between first and last points
2246 // 3. Loop over tracks - decreasing quality order
2247 // a.) remove - If the fraction of shared cluster less than factor (1- n or 2)
2248 // b.) remove - If the minimal number of clusters less than minimal and not ITS
2249 // c.) if track accepted - sign clusters
2251 //Called in - AliTPCtrackerMI::Clusters2Tracks()
2252 // - AliTPCtrackerMI::PropagateBack()
2253 // - AliTPCtrackerMI::RefitInward()
2256 // factor1 - factor for constrained
2257 // factor2 - for non constrained tracks
2258 // if (Float_t(shared+1)/Float_t(found+1)>factor) - DELETE
2262 Int_t nseed = arr->GetEntriesFast();
2263 Float_t * quality = new Float_t[nseed];
2264 Int_t * indexes = new Int_t[nseed];
2268 for (Int_t i=0; i<nseed; i++) {
2269 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2274 pt->UpdatePoints(); //select first last max dens points
2275 Float_t * points = pt->GetPoints();
2276 if (points[3]<0.8) quality[i] =-1;
2277 quality[i] = (points[2]-points[0])+pt->GetNumberOfClusters();
2278 //prefer high momenta tracks if overlaps
2279 quality[i] *= TMath::Sqrt(TMath::Abs(pt->Pt())+0.5);
2281 TMath::Sort(nseed,quality,indexes);
2284 for (Int_t itrack=0; itrack<nseed; itrack++) {
2285 Int_t trackindex = indexes[itrack];
2286 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(trackindex);
2289 if (quality[trackindex]<0){
2290 delete arr->RemoveAt(trackindex);
2295 Int_t first = Int_t(pt->GetPoints()[0]);
2296 Int_t last = Int_t(pt->GetPoints()[2]);
2297 Double_t factor = (pt->GetBConstrain()) ? factor1: factor2;
2299 Int_t found,foundable,shared;
2300 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
2301 // pt->GetClusterStatistic(0,160, found, foundable,shared,kFALSE);
2302 Bool_t itsgold =kFALSE;
2305 if (pt->GetESD()->GetITSclusters(dummy)>4) itsgold= kTRUE;
2309 if (Float_t(shared+1)/Float_t(found+1)>factor){
2310 if (pt->GetKinkIndexes()[0]!=0) continue; //don't remove tracks - part of the kinks
2311 if( AliTPCReconstructor::StreamLevel()>3){
2312 TTreeSRedirector &cstream = *fDebugStreamer;
2313 cstream<<"RemoveUsed"<<
2314 "iter="<<fIteration<<
2318 delete arr->RemoveAt(trackindex);
2321 if (pt->GetNumberOfClusters()<50&&(found-0.5*shared)<minimal){ //remove short tracks
2322 if (pt->GetKinkIndexes()[0]!=0) continue; //don't remove tracks - part of the kinks
2323 if( AliTPCReconstructor::StreamLevel()>3){
2324 TTreeSRedirector &cstream = *fDebugStreamer;
2325 cstream<<"RemoveShort"<<
2326 "iter="<<fIteration<<
2330 delete arr->RemoveAt(trackindex);
2336 //if (sharedfactor>0.4) continue;
2337 if (pt->GetKinkIndexes()[0]>0) continue;
2338 //Remove tracks with undefined properties - seems
2339 if (pt->GetSigmaY2()<kAlmost0) continue; // ? what is the origin ?
2341 for (Int_t i=first; i<last; i++) {
2342 Int_t index=pt->GetClusterIndex2(i);
2343 // if (index<0 || index&0x8000 ) continue;
2344 if (index<0 || index&0x8000 ) continue;
2345 AliTPCclusterMI *c= pt->GetClusterPointer(i);
2352 Info("RemoveUsed","\n*****\nNumber of good tracks after shared removal\t%d\n",fNtracks);
2358 void AliTPCtrackerMI::DumpClusters(Int_t iter, TObjArray *trackArray)
2361 // Dump clusters after reco
2362 // signed and unsigned cluster can be visualized
2363 // 1. Unsign all cluster
2364 // 2. Sign all used clusters
2367 Int_t nseed = trackArray->GetEntries();
2368 for (Int_t i=0; i<nseed; i++){
2369 AliTPCseed *pt=(AliTPCseed*)trackArray->UncheckedAt(i);
2373 Bool_t isKink=pt->GetKinkIndex(0)!=0;
2374 for (Int_t j=0; j<160; ++j) {
2375 Int_t index=pt->GetClusterIndex2(j);
2376 if (index<0) continue;
2377 AliTPCclusterMI *c= pt->GetClusterPointer(j);
2379 if (isKink) c->Use(100); // kink
2380 c->Use(10); // by default usage 10
2385 for (Int_t sec=0;sec<fkNIS;sec++){
2386 for (Int_t row=0;row<fInnerSec->GetNRows();row++){
2387 AliTPCclusterMI *cl = fInnerSec[sec][row].GetClusters1();
2388 for (Int_t icl =0;icl< fInnerSec[sec][row].GetN1();icl++){
2389 Float_t gx[3]; cl[icl].GetGlobalXYZ(gx);
2390 (*fDebugStreamer)<<"clDump"<<
2398 cl = fInnerSec[sec][row].GetClusters2();
2399 for (Int_t icl =0;icl< fInnerSec[sec][row].GetN2();icl++){
2400 Float_t gx[3]; cl[icl].GetGlobalXYZ(gx);
2401 (*fDebugStreamer)<<"clDump"<<
2412 for (Int_t sec=0;sec<fkNOS;sec++){
2413 for (Int_t row=0;row<fOuterSec->GetNRows();row++){
2414 AliTPCclusterMI *cl = fOuterSec[sec][row].GetClusters1();
2415 for (Int_t icl =0;icl< fOuterSec[sec][row].GetN1();icl++){
2416 Float_t gx[3]; cl[icl].GetGlobalXYZ(gx);
2417 (*fDebugStreamer)<<"clDump"<<
2425 cl = fOuterSec[sec][row].GetClusters2();
2426 for (Int_t icl =0;icl< fOuterSec[sec][row].GetN2();icl++){
2427 Float_t gx[3]; cl[icl].GetGlobalXYZ(gx);
2428 (*fDebugStreamer)<<"clDump"<<
2440 void AliTPCtrackerMI::UnsignClusters()
2443 // loop over all clusters and unsign them
2446 for (Int_t sec=0;sec<fkNIS;sec++){
2447 for (Int_t row=0;row<fInnerSec->GetNRows();row++){
2448 AliTPCclusterMI *cl = fInnerSec[sec][row].GetClusters1();
2449 for (Int_t icl =0;icl< fInnerSec[sec][row].GetN1();icl++)
2450 // if (cl[icl].IsUsed(10))
2452 cl = fInnerSec[sec][row].GetClusters2();
2453 for (Int_t icl =0;icl< fInnerSec[sec][row].GetN2();icl++)
2454 //if (cl[icl].IsUsed(10))
2459 for (Int_t sec=0;sec<fkNOS;sec++){
2460 for (Int_t row=0;row<fOuterSec->GetNRows();row++){
2461 AliTPCclusterMI *cl = fOuterSec[sec][row].GetClusters1();
2462 for (Int_t icl =0;icl< fOuterSec[sec][row].GetN1();icl++)
2463 //if (cl[icl].IsUsed(10))
2465 cl = fOuterSec[sec][row].GetClusters2();
2466 for (Int_t icl =0;icl< fOuterSec[sec][row].GetN2();icl++)
2467 //if (cl[icl].IsUsed(10))
2476 void AliTPCtrackerMI::SignClusters(const TObjArray * arr, Float_t fnumber, Float_t fdensity)
2479 //sign clusters to be "used"
2481 // snumber and sdensity sign number of sigmas - bellow mean value to be accepted
2482 // loop over "primaries"
2496 Int_t nseed = arr->GetEntriesFast();
2497 for (Int_t i=0; i<nseed; i++) {
2498 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2502 if (!(pt->IsActive())) continue;
2503 Float_t dens = pt->GetNumberOfClusters()/Float_t(pt->GetNFoundable());
2504 if ( (dens>0.7) && (pt->GetNumberOfClusters()>70)){
2506 sumdens2+= dens*dens;
2507 sumn += pt->GetNumberOfClusters();
2508 sumn2 += pt->GetNumberOfClusters()*pt->GetNumberOfClusters();
2509 Float_t chi2 = pt->GetChi2()/pt->GetNumberOfClusters();
2512 sumchi2 +=chi2*chi2;
2517 Float_t mdensity = 0.9;
2518 Float_t meann = 130;
2519 Float_t meanchi = 1;
2520 Float_t sdensity = 0.1;
2521 Float_t smeann = 10;
2522 Float_t smeanchi =0.4;
2526 mdensity = sumdens/sum;
2528 meanchi = sumchi/sum;
2530 sdensity = sumdens2/sum-mdensity*mdensity;
2532 sdensity = TMath::Sqrt(sdensity);
2536 smeann = sumn2/sum-meann*meann;
2538 smeann = TMath::Sqrt(smeann);
2542 smeanchi = sumchi2/sum - meanchi*meanchi;
2544 smeanchi = TMath::Sqrt(smeanchi);
2550 //REMOVE SHORT DELTAS or tracks going out of sensitive volume of TPC
2552 for (Int_t i=0; i<nseed; i++) {
2553 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2557 if (pt->GetBSigned()) continue;
2558 if (pt->GetBConstrain()) continue;
2559 //if (!(pt->IsActive())) continue;
2561 Int_t found,foundable,shared;
2562 pt->GetClusterStatistic(0,160,found, foundable,shared);
2563 if (shared/float(found)>0.3) {
2564 if (shared/float(found)>0.9 ){
2565 //delete arr->RemoveAt(i);
2570 Bool_t isok =kFALSE;
2571 if ( (pt->GetNShared()/pt->GetNumberOfClusters()<0.5) &&pt->GetNumberOfClusters()>60)
2573 if ((TMath::Abs(1/pt->GetC())<100.) && (pt->GetNShared()/pt->GetNumberOfClusters()<0.7))
2575 if (TMath::Abs(pt->GetZ()/pt->GetX())>1.1)
2577 if ( (TMath::Abs(pt->GetSnp()>0.7) && pt->GetD(0,0)>60.))
2581 for (Int_t j=0; j<160; ++j) {
2582 Int_t index=pt->GetClusterIndex2(j);
2583 if (index<0) continue;
2584 AliTPCclusterMI *c= pt->GetClusterPointer(j);
2586 //if (!(c->IsUsed(10))) c->Use();
2593 Double_t maxchi = meanchi+2.*smeanchi;
2595 for (Int_t i=0; i<nseed; i++) {
2596 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2600 //if (!(pt->IsActive())) continue;
2601 if (pt->GetBSigned()) continue;
2602 Double_t chi = pt->GetChi2()/pt->GetNumberOfClusters();
2603 if (chi>maxchi) continue;
2606 Float_t dens = pt->GetNumberOfClusters()/Float_t(pt->GetNFoundable());
2608 //sign only tracks with enoug big density at the beginning
2610 if ((pt->GetDensityFirst(40)<0.75) && pt->GetNumberOfClusters()<meann) continue;
2613 Double_t mindens = TMath::Max(double(mdensity-sdensity*fdensity*bfactor),0.65);
2614 Double_t minn = TMath::Max(Int_t(meann-fnumber*smeann*bfactor),50);
2616 // if (pt->fBConstrain) mindens = TMath::Max(mdensity-sdensity*fdensity*bfactor,0.65);
2617 if ( (pt->GetRemoval()==10) && (pt->GetSnp()>0.8)&&(dens>mindens))
2620 if ((dens>mindens && pt->GetNumberOfClusters()>minn) && chi<maxchi ){
2621 //Int_t noc=pt->GetNumberOfClusters();
2622 pt->SetBSigned(kTRUE);
2623 for (Int_t j=0; j<160; ++j) {
2625 Int_t index=pt->GetClusterIndex2(j);
2626 if (index<0) continue;
2627 AliTPCclusterMI *c= pt->GetClusterPointer(j);
2629 // if (!(c->IsUsed(10))) c->Use();
2634 // gLastCheck = nseed;
2642 void AliTPCtrackerMI::StopNotActive(const TObjArray * arr, Int_t row0, Float_t th0, Float_t th1, Float_t th2) const
2644 // stop not active tracks
2645 // take th1 as threshold for number of founded to number of foundable on last 10 active rows
2646 // take th2 as threshold for number of founded to number of foundable on last 20 active rows
2647 Int_t nseed = arr->GetEntriesFast();
2649 for (Int_t i=0; i<nseed; i++) {
2650 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2654 if (!(pt->IsActive())) continue;
2655 StopNotActive(pt,row0,th0, th1,th2);
2661 void AliTPCtrackerMI::StopNotActive(AliTPCseed * seed, Int_t row0, Float_t th0, Float_t th1,
2664 // stop not active tracks
2665 // take th1 as threshold for number of founded to number of foundable on last 10 active rows
2666 // take th2 as threshold for number of founded to number of foundable on last 20 active rows
2669 Int_t foundable = 0;
2670 Int_t maxindex = seed->GetLastPoint(); //last foundable row
2671 if (seed->GetNFoundable()*th0 > seed->GetNumberOfClusters()) {
2672 seed->Desactivate(10) ;
2676 for (Int_t i=row0; i<maxindex; i++){
2677 Int_t index = seed->GetClusterIndex2(i);
2678 if (index!=-1) foundable++;
2680 if (foundable<=30) sumgood1++;
2681 if (foundable<=50) {
2688 if (foundable>=30.){
2689 if (sumgood1<(th1*30.)) seed->Desactivate(10);
2692 if (sumgood2<(th2*50.)) seed->Desactivate(10);
2696 Int_t AliTPCtrackerMI::RefitInward(AliESDEvent *event)
2699 // back propagation of ESD tracks
2702 if (!event) return 0;
2703 const Int_t kMaxFriendTracks=2000;
2705 // extract correction object for multiplicity dependence of dEdx
2706 TObjArray * gainCalibArray = AliTPCcalibDB::Instance()->GetTimeGainSplinesRun(event->GetRunNumber());
2708 AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ;
2710 AliFatal("Tranformations not in RefitInward");
2713 transform->SetCurrentRecoParam((AliTPCRecoParam*)AliTPCReconstructor::GetRecoParam());
2714 const AliTPCRecoParam * recoParam = AliTPCcalibDB::Instance()->GetTransform()->GetCurrentRecoParam();
2715 Int_t nContribut = event->GetNumberOfTracks();
2716 TGraphErrors * graphMultDependenceDeDx = 0x0;
2717 if (recoParam && recoParam->GetUseMultiplicityCorrectionDedx() && gainCalibArray) {
2718 if (recoParam->GetUseTotCharge()) {
2719 graphMultDependenceDeDx = (TGraphErrors *) gainCalibArray->FindObject("TGRAPHERRORS_MEANQTOT_MULTIPLICITYDEPENDENCE_BEAM_ALL");
2721 graphMultDependenceDeDx = (TGraphErrors *) gainCalibArray->FindObject("TGRAPHERRORS_MEANQMAX_MULTIPLICITYDEPENDENCE_BEAM_ALL");
2727 //PrepareForProlongation(fSeeds,1);
2728 PropagateForward2(fSeeds);
2729 RemoveUsed2(fSeeds,0.4,0.4,20);
2731 TObjArray arraySeed(fSeeds->GetEntries());
2732 for (Int_t i=0;i<fSeeds->GetEntries();i++) {
2733 arraySeed.AddAt(fSeeds->At(i),i);
2735 SignShared(&arraySeed);
2736 // FindCurling(fSeeds, event,2); // find multi found tracks
2737 FindSplitted(fSeeds, event,2); // find multi found tracks
2738 if (AliTPCReconstructor::StreamLevel()>5) FindMultiMC(fSeeds, fEvent,2); // find multi found tracks
2741 Int_t nseed = fSeeds->GetEntriesFast();
2742 for (Int_t i=0;i<nseed;i++){
2743 AliTPCseed * seed = (AliTPCseed*) fSeeds->UncheckedAt(i);
2744 if (!seed) continue;
2745 if (seed->GetKinkIndex(0)>0) UpdateKinkQualityD(seed); // update quality informations for kinks
2746 AliESDtrack *esd=event->GetTrack(i);
2748 if (seed->GetNumberOfClusters()<60 && seed->GetNumberOfClusters()<(esd->GetTPCclusters(0) -5)*0.8){
2749 AliExternalTrackParam paramIn;
2750 AliExternalTrackParam paramOut;
2751 Int_t ncl = seed->RefitTrack(seed,¶mIn,¶mOut);
2752 if (AliTPCReconstructor::StreamLevel()>2) {
2753 (*fDebugStreamer)<<"RecoverIn"<<
2757 "pout.="<<¶mOut<<
2762 seed->Set(paramIn.GetX(),paramIn.GetAlpha(),paramIn.GetParameter(),paramIn.GetCovariance());
2763 seed->SetNumberOfClusters(ncl);
2767 seed->PropagateTo(fkParam->GetInnerRadiusLow());
2768 seed->UpdatePoints();
2769 AddCovariance(seed);
2770 MakeESDBitmaps(seed, esd);
2771 seed->CookdEdx(0.02,0.6);
2772 CookLabel(seed,0.1); //For comparison only
2774 if (AliTPCReconstructor::StreamLevel()>1 && seed!=0) {
2775 TTreeSRedirector &cstream = *fDebugStreamer;
2782 if (seed->GetNumberOfClusters()>15){
2783 esd->UpdateTrackParams(seed,AliESDtrack::kTPCrefit);
2784 esd->SetTPCPoints(seed->GetPoints());
2785 esd->SetTPCPointsF(seed->GetNFoundable());
2786 Int_t ndedx = seed->GetNCDEDX(0);
2787 Float_t sdedx = seed->GetSDEDX(0);
2788 Float_t dedx = seed->GetdEdx();
2789 // apply mutliplicity dependent dEdx correction if available
2790 if (graphMultDependenceDeDx) {
2791 Double_t corrGain = AliTPCcalibDButil::EvalGraphConst(graphMultDependenceDeDx, nContribut);
2792 dedx += (1 - corrGain)*50.; // MIP is normalized to 50
2794 esd->SetTPCsignal(dedx, sdedx, ndedx);
2796 // fill new dEdx information
2798 Double32_t signal[4];
2802 for(Int_t iarr=0;iarr<3;iarr++) {
2803 signal[iarr] = seed->GetDEDXregion(iarr+1);
2804 ncl[iarr] = seed->GetNCDEDX(iarr+1);
2805 nrows[iarr] = seed->GetNCDEDXInclThres(iarr+1);
2807 signal[3] = seed->GetDEDXregion(4);
2809 AliTPCdEdxInfo * infoTpcPid = new AliTPCdEdxInfo();
2810 infoTpcPid->SetTPCSignalRegionInfo(signal, ncl, nrows);
2811 esd->SetTPCdEdxInfo(infoTpcPid);
2813 // add seed to the esd track in Calib level
2815 Bool_t storeFriend = gRandom->Rndm()<(kMaxFriendTracks)/Float_t(nseed);
2816 if (AliTPCReconstructor::StreamLevel()>0 &&storeFriend){
2817 AliTPCseed * seedCopy = new AliTPCseed(*seed, kTRUE);
2818 esd->AddCalibObject(seedCopy);
2823 //printf("problem\n");
2826 //FindKinks(fSeeds,event);
2827 if (AliTPCReconstructor::StreamLevel()>3) DumpClusters(2,fSeeds);
2828 Info("RefitInward","Number of refitted tracks %d",ntracks);
2833 Int_t AliTPCtrackerMI::PropagateBack(AliESDEvent *event)
2836 // back propagation of ESD tracks
2838 if (!event) return 0;
2842 PropagateBack(fSeeds);
2843 RemoveUsed2(fSeeds,0.4,0.4,20);
2844 //FindCurling(fSeeds, fEvent,1);
2845 FindSplitted(fSeeds, event,1); // find multi found tracks
2846 if (AliTPCReconstructor::StreamLevel()>5) FindMultiMC(fSeeds, fEvent,1); // find multi found tracks
2849 Int_t nseed = fSeeds->GetEntriesFast();
2851 for (Int_t i=0;i<nseed;i++){
2852 AliTPCseed * seed = (AliTPCseed*) fSeeds->UncheckedAt(i);
2853 if (!seed) continue;
2854 if (seed->GetKinkIndex(0)<0) UpdateKinkQualityM(seed); // update quality informations for kinks
2855 seed->UpdatePoints();
2856 AddCovariance(seed);
2857 AliESDtrack *esd=event->GetTrack(i);
2858 if (!esd) continue; //never happen
2859 if (seed->GetNumberOfClusters()<60 && seed->GetNumberOfClusters()<(esd->GetTPCclusters(0) -5)*0.8){
2860 AliExternalTrackParam paramIn;
2861 AliExternalTrackParam paramOut;
2862 Int_t ncl = seed->RefitTrack(seed,¶mIn,¶mOut);
2863 if (AliTPCReconstructor::StreamLevel()>2) {
2864 (*fDebugStreamer)<<"RecoverBack"<<
2868 "pout.="<<¶mOut<<
2873 seed->Set(paramOut.GetX(),paramOut.GetAlpha(),paramOut.GetParameter(),paramOut.GetCovariance());
2874 seed->SetNumberOfClusters(ncl);
2877 seed->CookdEdx(0.02,0.6);
2878 CookLabel(seed,0.1); //For comparison only
2879 if (seed->GetNumberOfClusters()>15){
2880 esd->UpdateTrackParams(seed,AliESDtrack::kTPCout);
2881 esd->SetTPCPoints(seed->GetPoints());
2882 esd->SetTPCPointsF(seed->GetNFoundable());
2883 Int_t ndedx = seed->GetNCDEDX(0);
2884 Float_t sdedx = seed->GetSDEDX(0);
2885 Float_t dedx = seed->GetdEdx();
2886 esd->SetTPCsignal(dedx, sdedx, ndedx);
2888 Int_t eventnumber = event->GetEventNumberInFile();// patch 28 fev 06
2889 // This is most likely NOT the event number you'd like to use. It has nothing to do with the 'real' event number
2890 if (AliTPCReconstructor::StreamLevel()>1 && esd) {
2891 (*fDebugStreamer)<<"Cback"<<
2894 "EventNrInFile="<<eventnumber<<
2899 if (AliTPCReconstructor::StreamLevel()>3) DumpClusters(1,fSeeds);
2900 //FindKinks(fSeeds,event);
2901 Info("PropagateBack","Number of back propagated tracks %d",ntracks);
2908 void AliTPCtrackerMI::DeleteSeeds()
2913 Int_t nseed = fSeeds->GetEntriesFast();
2914 for (Int_t i=0;i<nseed;i++){
2915 AliTPCseed * seed = (AliTPCseed*)fSeeds->At(i);
2916 if (seed) delete fSeeds->RemoveAt(i);
2923 void AliTPCtrackerMI::ReadSeeds(const AliESDEvent *const event, Int_t direction)
2926 //read seeds from the event
2928 Int_t nentr=event->GetNumberOfTracks();
2930 Info("ReadSeeds", "Number of ESD tracks: %d\n", nentr);
2935 fSeeds = new TObjArray(nentr);
2939 for (Int_t i=0; i<nentr; i++) {
2940 AliESDtrack *esd=event->GetTrack(i);
2941 ULong_t status=esd->GetStatus();
2942 if (!(status&AliESDtrack::kTPCin)) continue;
2943 AliTPCtrack t(*esd);
2944 t.SetNumberOfClusters(0);
2945 // AliTPCseed *seed = new AliTPCseed(t,t.GetAlpha());
2946 AliTPCseed *seed = new AliTPCseed(t/*,t.GetAlpha()*/);
2947 seed->SetUniqueID(esd->GetID());
2948 AddCovariance(seed); //add systematic ucertainty
2949 for (Int_t ikink=0;ikink<3;ikink++) {
2950 Int_t index = esd->GetKinkIndex(ikink);
2951 seed->GetKinkIndexes()[ikink] = index;
2952 if (index==0) continue;
2953 index = TMath::Abs(index);
2954 AliESDkink * kink = fEvent->GetKink(index-1);
2955 if (kink&&esd->GetKinkIndex(ikink)<0){
2956 if ((status & AliESDtrack::kTRDrefit) != 0) kink->SetStatus(1,2);
2957 if ((status & AliESDtrack::kITSout) != 0) kink->SetStatus(1,0);
2959 if (kink&&esd->GetKinkIndex(ikink)>0){
2960 if ((status & AliESDtrack::kTRDrefit) != 0) kink->SetStatus(1,6);
2961 if ((status & AliESDtrack::kITSout) != 0) kink->SetStatus(1,4);
2965 if (((status&AliESDtrack::kITSout)==0)&&(direction==1)) seed->ResetCovariance(10.);
2966 if ( direction ==2 &&(status & AliESDtrack::kTRDrefit) == 0 ) seed->ResetCovariance(10.);
2967 //if ( direction ==2 && ((status & AliESDtrack::kTPCout) == 0) ) {
2968 // fSeeds->AddAt(0,i);
2972 if ( direction ==2 &&(status & AliESDtrack::kTRDrefit) > 0 ) {
2973 Double_t par0[5],par1[5],alpha,x;
2974 esd->GetInnerExternalParameters(alpha,x,par0);
2975 esd->GetExternalParameters(x,par1);
2976 Double_t delta1 = TMath::Abs(par0[4]-par1[4])/(0.000000001+TMath::Abs(par0[4]+par1[4]));
2977 Double_t delta2 = TMath::Abs(par0[3]-par1[3]);
2979 if (esd->GetTRDncls()>0) trdchi2 = esd->GetTRDchi2()/esd->GetTRDncls();
2980 //reset covariance if suspicious
2981 if ( (delta1>0.1) || (delta2>0.006) ||trdchi2>7.)
2982 seed->ResetCovariance(10.);
2987 // rotate to the local coordinate system
2989 fSectors=fInnerSec; fN=fkNIS;
2990 Double_t alpha=seed->GetAlpha();
2991 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
2992 if (alpha < 0. ) alpha += 2.*TMath::Pi();
2993 Int_t ns=Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN;
2994 alpha =ns*fSectors->GetAlpha() + fSectors->GetAlphaShift();
2995 alpha-=seed->GetAlpha();
2996 if (alpha<-TMath::Pi()) alpha += 2*TMath::Pi();
2997 if (alpha>=TMath::Pi()) alpha -= 2*TMath::Pi();
2998 if (TMath::Abs(alpha) > 0.001) { // This should not happen normally
2999 AliWarning(Form("Rotating track over %f",alpha));
3000 if (!seed->Rotate(alpha)) {
3007 if (esd->GetKinkIndex(0)<=0){
3008 for (Int_t irow=0;irow<160;irow++){
3009 Int_t index = seed->GetClusterIndex2(irow);
3012 AliTPCclusterMI * cl = GetClusterMI(index);
3013 seed->SetClusterPointer(irow,cl);
3015 if ((index & 0x8000)==0){
3016 cl->Use(10); // accepted cluster
3018 cl->Use(6); // close cluster not accepted
3021 Info("ReadSeeds","Not found cluster");
3026 fSeeds->AddAt(seed,i);
3032 //_____________________________________________________________________________
3033 void AliTPCtrackerMI::MakeSeeds3(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, Float_t cuts[4],
3034 Float_t deltay, Int_t ddsec) {
3035 //-----------------------------------------------------------------
3036 // This function creates track seeds.
3037 // SEEDING WITH VERTEX CONSTRAIN
3038 //-----------------------------------------------------------------
3039 // cuts[0] - fP4 cut
3040 // cuts[1] - tan(phi) cut
3041 // cuts[2] - zvertex cut
3042 // cuts[3] - fP3 cut
3050 Double_t x[5], c[15];
3051 // Int_t di = i1-i2;
3053 AliTPCseed * seed = new AliTPCseed();
3054 Double_t alpha=fSectors->GetAlpha(), shift=fSectors->GetAlphaShift();
3055 Double_t cs=cos(alpha), sn=sin(alpha);
3057 // Double_t x1 =fOuterSec->GetX(i1);
3058 //Double_t xx2=fOuterSec->GetX(i2);
3060 Double_t x1 =GetXrow(i1);
3061 Double_t xx2=GetXrow(i2);
3063 Double_t x3=GetX(), y3=GetY(), z3=GetZ();
3065 Int_t imiddle = (i2+i1)/2; //middle pad row index
3066 Double_t xm = GetXrow(imiddle); // radius of middle pad-row
3067 const AliTPCtrackerRow& krm=GetRow(sec,imiddle); //middle pad -row
3071 const AliTPCtrackerRow& kr1=GetRow(ns,i1);
3072 Double_t ymax = GetMaxY(i1)-kr1.GetDeadZone()-1.5;
3073 Double_t ymaxm = GetMaxY(imiddle)-kr1.GetDeadZone()-1.5;
3076 // change cut on curvature if it can't reach this layer
3077 // maximal curvature set to reach it
3078 Double_t dvertexmax = TMath::Sqrt((x1-x3)*(x1-x3)+(ymax+5-y3)*(ymax+5-y3));
3079 if (dvertexmax*0.5*cuts[0]>0.85){
3080 cuts[0] = 0.85/(dvertexmax*0.5+1.);
3082 Double_t r2min = 1/(cuts[0]*cuts[0]); //minimal square of radius given by cut
3085 if (deltay>0) ddsec = 0;
3086 // loop over clusters
3087 for (Int_t is=0; is < kr1; is++) {
3089 if (kr1[is]->IsUsed(10)) continue;
3090 Double_t y1=kr1[is]->GetY(), z1=kr1[is]->GetZ();
3091 //if (TMath::Abs(y1)>ymax) continue;
3093 if (deltay>0 && TMath::Abs(ymax-TMath::Abs(y1))> deltay ) continue; // seed only at the edge
3095 // find possible directions
3096 Float_t anglez = (z1-z3)/(x1-x3);
3097 Float_t extraz = z1 - anglez*(x1-xx2); // extrapolated z
3100 //find rotation angles relative to line given by vertex and point 1
3101 Double_t dvertex2 = (x1-x3)*(x1-x3)+(y1-y3)*(y1-y3);
3102 Double_t dvertex = TMath::Sqrt(dvertex2);
3103 Double_t angle13 = TMath::ATan((y1-y3)/(x1-x3));
3104 Double_t cs13 = cos(-angle13), sn13 = sin(-angle13);
3107 // loop over 2 sectors
3113 Double_t dddz1=0; // direction of delta inclination in z axis
3120 for (Int_t dsec = dsec1; dsec<=dsec2;dsec++){
3121 Int_t sec2 = sec + dsec;
3123 // AliTPCtrackerRow& kr2 = fOuterSec[(sec2+fkNOS)%fkNOS][i2];
3124 //AliTPCtrackerRow& kr2m = fOuterSec[(sec2+fkNOS)%fkNOS][imiddle];
3125 AliTPCtrackerRow& kr2 = GetRow((sec2+fkNOS)%fkNOS,i2);
3126 AliTPCtrackerRow& kr2m = GetRow((sec2+fkNOS)%fkNOS,imiddle);
3127 Int_t index1 = TMath::Max(kr2.Find(extraz-0.6-dddz1*TMath::Abs(z1)*0.05)-1,0);
3128 Int_t index2 = TMath::Min(kr2.Find(extraz+0.6+dddz2*TMath::Abs(z1)*0.05)+1,kr2);
3130 // rotation angles to p1-p3
3131 Double_t cs13r = cos(-angle13+dsec*alpha)/dvertex, sn13r = sin(-angle13+dsec*alpha)/dvertex;
3132 Double_t x2, y2, z2;
3134 // Double_t dymax = maxangle*TMath::Abs(x1-xx2);
3137 Double_t dxx0 = (xx2-x3)*cs13r;
3138 Double_t dyy0 = (xx2-x3)*sn13r;
3139 for (Int_t js=index1; js < index2; js++) {
3140 const AliTPCclusterMI *kcl = kr2[js];
3141 if (kcl->IsUsed(10)) continue;
3143 //calcutate parameters
3145 Double_t yy0 = dyy0 +(kcl->GetY()-y3)*cs13r;
3147 if (TMath::Abs(yy0)<0.000001) continue;
3148 Double_t xx0 = dxx0 -(kcl->GetY()-y3)*sn13r;
3149 Double_t y0 = 0.5*(xx0*xx0+yy0*yy0-xx0)/yy0;
3150 Double_t r02 = (0.25+y0*y0)*dvertex2;
3151 //curvature (radius) cut
3152 if (r02<r2min) continue;
3156 Double_t c0 = 1/TMath::Sqrt(r02);
3160 //Double_t dfi0 = 2.*TMath::ASin(dvertex*c0*0.5);
3161 //Double_t dfi1 = 2.*TMath::ASin(TMath::Sqrt(yy0*yy0+(1-xx0)*(1-xx0))*dvertex*c0*0.5);
3162 Double_t dfi0 = 2.*AliTPCFastMath::FastAsin(dvertex*c0*0.5);
3163 Double_t dfi1 = 2.*AliTPCFastMath::FastAsin(TMath::Sqrt(yy0*yy0+(1-xx0)*(1-xx0))*dvertex*c0*0.5);
3166 Double_t z0 = kcl->GetZ();
3167 Double_t zzzz2 = z1-(z1-z3)*dfi1/dfi0;
3168 if (TMath::Abs(zzzz2-z0)>0.5) continue;
3171 Double_t dip = (z1-z0)*c0/dfi1;
3172 Double_t x0 = (0.5*cs13+y0*sn13)*dvertex*c0;
3183 x2= xx2*cs-y2*sn*dsec;
3184 y2=+xx2*sn*dsec+y2*cs;
3194 // do we have cluster at the middle ?
3196 GetProlongation(x1,xm,x,ym,zm);
3198 AliTPCclusterMI * cm=0;
3199 if (TMath::Abs(ym)-ymaxm<0){
3200 cm = krm.FindNearest2(ym,zm,1.0,0.6,dummy);
3201 if ((!cm) || (cm->IsUsed(10))) {
3206 // rotate y1 to system 0
3207 // get state vector in rotated system
3208 Double_t yr1 = (-0.5*sn13+y0*cs13)*dvertex*c0;
3209 Double_t xr2 = x0*cs+yr1*sn*dsec;
3210 Double_t xr[5]={kcl->GetY(),kcl->GetZ(), xr2, dip, c0};
3212 GetProlongation(xx2,xm,xr,ym,zm);
3213 if (TMath::Abs(ym)-ymaxm<0){
3214 cm = kr2m.FindNearest2(ym,zm,1.0,0.6,dummy);
3215 if ((!cm) || (cm->IsUsed(10))) {
3225 dym = ym - cm->GetY();
3226 dzm = zm - cm->GetZ();
3233 Double_t sy1=kr1[is]->GetSigmaY2()*2., sz1=kr1[is]->GetSigmaZ2()*2.;
3234 Double_t sy2=kcl->GetSigmaY2()*2., sz2=kcl->GetSigmaZ2()*2.;
3235 //Double_t sy3=400*3./12., sy=0.1, sz=0.1;
3236 Double_t sy3=25000*x[4]*x[4]+0.1, sy=0.1, sz=0.1;
3237 //Double_t sy3=25000*x[4]*x[4]*60+0.5, sy=0.1, sz=0.1;
3239 Double_t f40=(F1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
3240 Double_t f42=(F1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
3241 Double_t f43=(F1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
3242 Double_t f20=(F2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
3243 Double_t f22=(F2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
3244 Double_t f23=(F2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
3246 Double_t f30=(F3(x1,y1+sy,x2,y2,z1,z2)-x[3])/sy;
3247 Double_t f31=(F3(x1,y1,x2,y2,z1+sz,z2)-x[3])/sz;
3248 Double_t f32=(F3(x1,y1,x2,y2+sy,z1,z2)-x[3])/sy;
3249 Double_t f34=(F3(x1,y1,x2,y2,z1,z2+sz)-x[3])/sz;
3253 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3254 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3255 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3256 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3257 c[13]=f30*sy1*f40+f32*sy2*f42;
3258 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3260 // if (!BuildSeed(kr1[is],kcl,0,x1,x2,x3,x,c)) continue;
3262 UInt_t index=kr1.GetIndex(is);
3263 seed->~AliTPCseed(); // this does not set the pointer to 0...
3264 AliTPCseed *track=new(seed) AliTPCseed(x1, ns*alpha+shift, x, c, index);
3266 track->SetIsSeeding(kTRUE);
3267 track->SetSeed1(i1);
3268 track->SetSeed2(i2);
3269 track->SetSeedType(3);
3273 FollowProlongation(*track, (i1+i2)/2,1);
3274 Int_t foundable,found,shared;
3275 track->GetClusterStatistic((i1+i2)/2,i1, found, foundable, shared, kTRUE);
3276 if ((found<0.55*foundable) || shared>0.5*found || (track->GetSigmaY2()+track->GetSigmaZ2())>0.5){
3278 seed->~AliTPCseed();
3284 FollowProlongation(*track, i2,1);
3288 track->SetBConstrain(1);
3289 // track->fLastPoint = i1+fInnerSec->GetNRows(); // first cluster in track position
3290 track->SetLastPoint(i1); // first cluster in track position
3291 track->SetFirstPoint(track->GetLastPoint());
3293 if (track->GetNumberOfClusters()<(i1-i2)*0.5 ||
3294 track->GetNumberOfClusters() < track->GetNFoundable()*0.6 ||
3295 track->GetNShared()>0.4*track->GetNumberOfClusters() ) {
3297 seed->~AliTPCseed();
3301 // Z VERTEX CONDITION
3302 Double_t zv, bz=GetBz();
3303 if ( !track->GetZAt(0.,bz,zv) ) continue;
3304 if (TMath::Abs(zv-z3)>cuts[2]) {
3305 FollowProlongation(*track, TMath::Max(i2-20,0));
3306 if ( !track->GetZAt(0.,bz,zv) ) continue;
3307 if (TMath::Abs(zv-z3)>cuts[2]){
3308 FollowProlongation(*track, TMath::Max(i2-40,0));
3309 if ( !track->GetZAt(0.,bz,zv) ) continue;
3310 if (TMath::Abs(zv-z3)>cuts[2] &&(track->GetNumberOfClusters() > track->GetNFoundable()*0.7)){
3311 // make seed without constrain
3312 AliTPCseed * track2 = MakeSeed(track,0.2,0.5,1.);
3313 FollowProlongation(*track2, i2,1);
3314 track2->SetBConstrain(kFALSE);
3315 track2->SetSeedType(1);
3316 arr->AddLast(track2);
3318 seed->~AliTPCseed();
3323 seed->~AliTPCseed();
3330 track->SetSeedType(0);
3331 arr->AddLast(track);
3332 seed = new AliTPCseed;
3334 // don't consider other combinations
3335 if (track->GetNumberOfClusters() > track->GetNFoundable()*0.8)
3341 Info("MakeSeeds3","\nSeeding statistic:\t%d\t%d\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin,nout1,nout2);
3347 void AliTPCtrackerMI::MakeSeeds5(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, Float_t cuts[4],
3352 //-----------------------------------------------------------------
3353 // This function creates track seeds.
3354 //-----------------------------------------------------------------
3355 // cuts[0] - fP4 cut
3356 // cuts[1] - tan(phi) cut
3357 // cuts[2] - zvertex cut
3358 // cuts[3] - fP3 cut
3368 Double_t x[5], c[15];
3370 // make temporary seed
3371 AliTPCseed * seed = new AliTPCseed;
3372 Double_t alpha=fOuterSec->GetAlpha(), shift=fOuterSec->GetAlphaShift();
3373 // Double_t cs=cos(alpha), sn=sin(alpha);
3378 Double_t x1 = GetXrow(i1-1);
3379 const AliTPCtrackerRow& kr1=GetRow(sec,i1-1);
3380 Double_t y1max = GetMaxY(i1-1)-kr1.GetDeadZone()-1.5;
3382 Double_t x1p = GetXrow(i1);
3383 const AliTPCtrackerRow& kr1p=GetRow(sec,i1);
3385 Double_t x1m = GetXrow(i1-2);
3386 const AliTPCtrackerRow& kr1m=GetRow(sec,i1-2);
3389 //last 3 padrow for seeding
3390 AliTPCtrackerRow& kr3 = GetRow((sec+fkNOS)%fkNOS,i1-7);
3391 Double_t x3 = GetXrow(i1-7);
3392 // Double_t y3max= GetMaxY(i1-7)-kr3.fDeadZone-1.5;
3394 AliTPCtrackerRow& kr3p = GetRow((sec+fkNOS)%fkNOS,i1-6);
3395 Double_t x3p = GetXrow(i1-6);
3397 AliTPCtrackerRow& kr3m = GetRow((sec+fkNOS)%fkNOS,i1-8);
3398 Double_t x3m = GetXrow(i1-8);
3403 Int_t im = i1-4; //middle pad row index
3404 Double_t xm = GetXrow(im); // radius of middle pad-row
3405 const AliTPCtrackerRow& krm=GetRow(sec,im); //middle pad -row
3406 // Double_t ymmax = GetMaxY(im)-kr1.fDeadZone-1.5;
3409 Double_t deltax = x1-x3;
3410 Double_t dymax = deltax*cuts[1];
3411 Double_t dzmax = deltax*cuts[3];
3413 // loop over clusters
3414 for (Int_t is=0; is < kr1; is++) {
3416 if (kr1[is]->IsUsed(10)) continue;
3417 Double_t y1=kr1[is]->GetY(), z1=kr1[is]->GetZ();
3419 if (deltay>0 && TMath::Abs(y1max-TMath::Abs(y1))> deltay ) continue; // seed only at the edge
3421 Int_t index1 = TMath::Max(kr3.Find(z1-dzmax)-1,0);
3422 Int_t index2 = TMath::Min(kr3.Find(z1+dzmax)+1,kr3);
3428 for (Int_t js=index1; js < index2; js++) {
3429 const AliTPCclusterMI *kcl = kr3[js];
3430 if (kcl->IsUsed(10)) continue;
3432 // apply angular cuts
3433 if (TMath::Abs(y1-y3)>dymax) continue;
3436 if (TMath::Abs(z1-z3)>dzmax) continue;
3438 Double_t angley = (y1-y3)/(x1-x3);
3439 Double_t anglez = (z1-z3)/(x1-x3);
3441 Double_t erry = TMath::Abs(angley)*(x1-x1m)*0.5+0.5;
3442 Double_t errz = TMath::Abs(anglez)*(x1-x1m)*0.5+0.5;
3444 Double_t yyym = angley*(xm-x1)+y1;
3445 Double_t zzzm = anglez*(xm-x1)+z1;
3447 const AliTPCclusterMI *kcm = krm.FindNearest2(yyym,zzzm,erry,errz,index);
3449 if (kcm->IsUsed(10)) continue;
3451 erry = TMath::Abs(angley)*(x1-x1m)*0.4+0.5;
3452 errz = TMath::Abs(anglez)*(x1-x1m)*0.4+0.5;
3459 // look around first
3460 const AliTPCclusterMI *kc1m = kr1m.FindNearest2(angley*(x1m-x1)+y1,
3466 if (kc1m->IsUsed(10)) used++;
3468 const AliTPCclusterMI *kc1p = kr1p.FindNearest2(angley*(x1p-x1)+y1,
3474 if (kc1p->IsUsed(10)) used++;
3476 if (used>1) continue;
3477 if (found<1) continue;
3481 const AliTPCclusterMI *kc3m = kr3m.FindNearest2(angley*(x3m-x3)+y3,
3487 if (kc3m->IsUsed(10)) used++;
3491 const AliTPCclusterMI *kc3p = kr3p.FindNearest2(angley*(x3p-x3)+y3,
3497 if (kc3p->IsUsed(10)) used++;
3501 if (used>1) continue;
3502 if (found<3) continue;
3512 x[4]=F1(x1,y1,x2,y2,x3,y3);
3513 //if (TMath::Abs(x[4]) >= cuts[0]) continue;
3516 x[2]=F2(x1,y1,x2,y2,x3,y3);
3519 x[3]=F3n(x1,y1,x2,y2,z1,z2,x[4]);
3520 //if (TMath::Abs(x[3]) > cuts[3]) continue;
3524 Double_t sy1=0.1, sz1=0.1;
3525 Double_t sy2=0.1, sz2=0.1;
3526 Double_t sy3=0.1, sy=0.1, sz=0.1;
3528 Double_t f40=(F1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
3529 Double_t f42=(F1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
3530 Double_t f43=(F1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
3531 Double_t f20=(F2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
3532 Double_t f22=(F2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
3533 Double_t f23=(F2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
3535 Double_t f30=(F3(x1,y1+sy,x2,y2,z1,z2)-x[3])/sy;
3536 Double_t f31=(F3(x1,y1,x2,y2,z1+sz,z2)-x[3])/sz;
3537 Double_t f32=(F3(x1,y1,x2,y2+sy,z1,z2)-x[3])/sy;
3538 Double_t f34=(F3(x1,y1,x2,y2,z1,z2+sz)-x[3])/sz;
3542 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3543 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3544 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3545 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3546 c[13]=f30*sy1*f40+f32*sy2*f42;
3547 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3549 // if (!BuildSeed(kr1[is],kcl,0,x1,x2,x3,x,c)) continue;
3551 index=kr1.GetIndex(is);
3552 seed->~AliTPCseed();
3553 AliTPCseed *track=new(seed) AliTPCseed(x1, sec*alpha+shift, x, c, index);
3555 track->SetIsSeeding(kTRUE);
3558 FollowProlongation(*track, i1-7,1);
3559 if (track->GetNumberOfClusters() < track->GetNFoundable()*0.75 ||
3560 track->GetNShared()>0.6*track->GetNumberOfClusters() || ( track->GetSigmaY2()+ track->GetSigmaZ2())>0.6){
3562 seed->~AliTPCseed();
3568 FollowProlongation(*track, i2,1);
3569 track->SetBConstrain(0);
3570 track->SetLastPoint(i1+fInnerSec->GetNRows()); // first cluster in track position
3571 track->SetFirstPoint(track->GetLastPoint());
3573 if (track->GetNumberOfClusters()<(i1-i2)*0.5 ||
3574 track->GetNumberOfClusters()<track->GetNFoundable()*0.7 ||
3575 track->GetNShared()>2. || track->GetChi2()/track->GetNumberOfClusters()>6 || ( track->GetSigmaY2()+ track->GetSigmaZ2())>0.5 ) {
3577 seed->~AliTPCseed();
3582 FollowProlongation(*track, TMath::Max(i2-10,0),1);
3583 AliTPCseed * track2 = MakeSeed(track,0.2,0.5,0.9);
3584 FollowProlongation(*track2, i2,1);
3585 track2->SetBConstrain(kFALSE);
3586 track2->SetSeedType(4);
3587 arr->AddLast(track2);
3589 seed->~AliTPCseed();
3593 //arr->AddLast(track);
3594 //seed = new AliTPCseed;
3600 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);
3606 //_____________________________________________________________________________
3607 void AliTPCtrackerMI::MakeSeeds2(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, Float_t */*cuts[4]*/,
3608 Float_t deltay, Bool_t /*bconstrain*/) {
3609 //-----------------------------------------------------------------
3610 // This function creates track seeds - without vertex constraint
3611 //-----------------------------------------------------------------
3612 // cuts[0] - fP4 cut - not applied
3613 // cuts[1] - tan(phi) cut
3614 // cuts[2] - zvertex cut - not applied
3615 // cuts[3] - fP3 cut
3625 Double_t alpha=fOuterSec->GetAlpha(), shift=fOuterSec->GetAlphaShift();
3626 // Double_t cs=cos(alpha), sn=sin(alpha);
3627 Int_t row0 = (i1+i2)/2;
3628 Int_t drow = (i1-i2)/2;
3629 const AliTPCtrackerRow& kr0=fSectors[sec][row0];
3630 AliTPCtrackerRow * kr=0;
3632 AliTPCpolyTrack polytrack;
3633 Int_t nclusters=fSectors[sec][row0];
3634 AliTPCseed * seed = new AliTPCseed;
3639 for (Int_t is=0; is < nclusters; is++) { //LOOP over clusters
3641 Int_t nfoundable =0;
3642 for (Int_t iter =1; iter<2; iter++){ //iterations
3643 const AliTPCtrackerRow& krm=fSectors[sec][row0-iter];
3644 const AliTPCtrackerRow& krp=fSectors[sec][row0+iter];
3645 const AliTPCclusterMI * cl= kr0[is];
3647 if (cl->IsUsed(10)) {
3653 Double_t x = kr0.GetX();
3654 // Initialization of the polytrack
3659 Double_t y0= cl->GetY();
3660 Double_t z0= cl->GetZ();
3664 Double_t ymax = fSectors->GetMaxY(row0)-kr0.GetDeadZone()-1.5;
3665 if (deltay>0 && TMath::Abs(ymax-TMath::Abs(y0))> deltay ) continue; // seed only at the edge
3667 erry = (0.5)*cl->GetSigmaY2()/TMath::Sqrt(cl->GetQ())*6;
3668 errz = (0.5)*cl->GetSigmaZ2()/TMath::Sqrt(cl->GetQ())*6;
3669 polytrack.AddPoint(x,y0,z0,erry, errz);
3672 if (cl->IsUsed(10)) sumused++;
3675 Float_t roady = (5*TMath::Sqrt(cl->GetSigmaY2()+0.2)+1.)*iter;
3676 Float_t roadz = (5*TMath::Sqrt(cl->GetSigmaZ2()+0.2)+1.)*iter;
3679 AliTPCclusterMI * cl1 = krm.FindNearest(y0,z0,roady,roadz);
3680 if (cl1 && TMath::Abs(ymax-TMath::Abs(y0))) {
3681 erry = (0.5)*cl1->GetSigmaY2()/TMath::Sqrt(cl1->GetQ())*3;
3682 errz = (0.5)*cl1->GetSigmaZ2()/TMath::Sqrt(cl1->GetQ())*3;
3683 if (cl1->IsUsed(10)) sumused++;
3684 polytrack.AddPoint(x,cl1->GetY(),cl1->GetZ(),erry,errz);
3688 AliTPCclusterMI * cl2 = krp.FindNearest(y0,z0,roady,roadz);
3690 erry = (0.5)*cl2->GetSigmaY2()/TMath::Sqrt(cl2->GetQ())*3;
3691 errz = (0.5)*cl2->GetSigmaZ2()/TMath::Sqrt(cl2->GetQ())*3;
3692 if (cl2->IsUsed(10)) sumused++;
3693 polytrack.AddPoint(x,cl2->GetY(),cl2->GetZ(),erry,errz);
3696 if (sumused>0) continue;
3698 polytrack.UpdateParameters();
3704 nfoundable = polytrack.GetN();
3705 nfound = nfoundable;
3707 for (Int_t ddrow = iter+1; ddrow<drow;ddrow++){
3708 Float_t maxdist = 0.8*(1.+3./(ddrow));
3709 for (Int_t delta = -1;delta<=1;delta+=2){
3710 Int_t row = row0+ddrow*delta;
3711 kr = &(fSectors[sec][row]);
3712 Double_t xn = kr->GetX();
3713 Double_t ymax1 = fSectors->GetMaxY(row)-kr->GetDeadZone()-1.5;
3714 polytrack.GetFitPoint(xn,yn,zn);
3715 if (TMath::Abs(yn)>ymax1) continue;
3717 AliTPCclusterMI * cln = kr->FindNearest(yn,zn,roady,roadz);
3719 Float_t dist = TMath::Sqrt( (yn-cln->GetY())*(yn-cln->GetY())+(zn-cln->GetZ())*(zn-cln->GetZ()));
3722 erry = (dist+0.3)*cln->GetSigmaY2()/TMath::Sqrt(cln->GetQ())*(1.+1./(ddrow));
3723 errz = (dist+0.3)*cln->GetSigmaZ2()/TMath::Sqrt(cln->GetQ())*(1.+1./(ddrow));
3724 if (cln->IsUsed(10)) {
3725 // printf("used\n");
3733 polytrack.AddPoint(xn,cln->GetY(),cln->GetZ(),erry, errz);
3738 if ( (sumused>3) || (sumused>0.5*nfound) || (nfound<0.6*nfoundable)) break;
3739 polytrack.UpdateParameters();
3742 if ( (sumused>3) || (sumused>0.5*nfound)) {
3743 //printf("sumused %d\n",sumused);
3748 polytrack.GetFitDerivation(kr0.GetX(),dy,dz);
3749 AliTPCpolyTrack track2;
3751 polytrack.Refit(track2,0.5+TMath::Abs(dy)*0.3,0.4+TMath::Abs(dz)*0.3);
3752 if (track2.GetN()<0.5*nfoundable) continue;
3755 if ((nfound>0.6*nfoundable) &&( nfoundable>0.4*(i1-i2))) {
3757 // test seed with and without constrain
3758 for (Int_t constrain=0; constrain<=0;constrain++){
3759 // add polytrack candidate
3761 Double_t x[5], c[15];
3762 Double_t x1,x2,x3,y1,y2,y3,z1,z2,z3;
3763 track2.GetBoundaries(x3,x1);
3765 track2.GetFitPoint(x1,y1,z1);
3766 track2.GetFitPoint(x2,y2,z2);
3767 track2.GetFitPoint(x3,y3,z3);
3769 //is track pointing to the vertex ?
3772 polytrack.GetFitPoint(x0,y0,z0);
3785 x[4]=F1(x1,y1,x2,y2,x3,y3);
3787 // if (TMath::Abs(x[4]) >= cuts[0]) continue; //
3788 x[2]=F2(x1,y1,x2,y2,x3,y3);
3790 //if (TMath::Abs(x[4]*x1-x[2]) >= cuts[1]) continue;
3791 //x[3]=F3(x1,y1,x2,y2,z1,z2);
3792 x[3]=F3n(x1,y1,x3,y3,z1,z3,x[4]);
3793 //if (TMath::Abs(x[3]) > cuts[3]) continue;
3796 Double_t sy =0.1, sz =0.1;
3797 Double_t sy1=0.02, sz1=0.02;
3798 Double_t sy2=0.02, sz2=0.02;
3802 sy3=25000*x[4]*x[4]+0.1, sy=0.1, sz=0.1;
3805 Double_t f40=(F1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
3806 Double_t f42=(F1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
3807 Double_t f43=(F1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
3808 Double_t f20=(F2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
3809 Double_t f22=(F2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
3810 Double_t f23=(F2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
3812 Double_t f30=(F3(x1,y1+sy,x3,y3,z1,z3)-x[3])/sy;
3813 Double_t f31=(F3(x1,y1,x3,y3,z1+sz,z3)-x[3])/sz;
3814 Double_t f32=(F3(x1,y1,x3,y3+sy,z1,z3)-x[3])/sy;
3815 Double_t f34=(F3(x1,y1,x3,y3,z1,z3+sz)-x[3])/sz;
3820 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3821 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3822 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3823 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3824 c[13]=f30*sy1*f40+f32*sy2*f42;
3825 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3827 //Int_t row1 = fSectors->GetRowNumber(x1);
3828 Int_t row1 = GetRowNumber(x1);
3832 seed->~AliTPCseed();
3833 AliTPCseed *track=new(seed) AliTPCseed(x1,sec*alpha+shift,x,c,index);
3834 track->SetIsSeeding(kTRUE);
3835 Int_t rc=FollowProlongation(*track, i2);
3836 if (constrain) track->SetBConstrain(1);
3838 track->SetBConstrain(0);
3839 track->SetLastPoint(row1+fInnerSec->GetNRows()); // first cluster in track position
3840 track->SetFirstPoint(track->GetLastPoint());
3842 if (rc==0 || track->GetNumberOfClusters()<(i1-i2)*0.5 ||
3843 track->GetNumberOfClusters() < track->GetNFoundable()*0.6 ||
3844 track->GetNShared()>0.4*track->GetNumberOfClusters()) {
3847 seed->~AliTPCseed();
3850 arr->AddLast(track);
3851 seed = new AliTPCseed;
3855 } // if accepted seed
3858 Info("MakeSeeds2","\nSeeding statiistic:\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin3);
3864 AliTPCseed *AliTPCtrackerMI::MakeSeed(AliTPCseed *const track, Float_t r0, Float_t r1, Float_t r2)
3868 //reseed using track points
3869 Int_t p0 = int(r0*track->GetNumberOfClusters()); // point 0
3870 Int_t p1 = int(r1*track->GetNumberOfClusters());
3871 Int_t p2 = int(r2*track->GetNumberOfClusters()); // last point
3873 Double_t x0[3],x1[3],x2[3];
3874 for (Int_t i=0;i<3;i++){
3880 // find track position at given ratio of the length
3881 Int_t sec0=0, sec1=0, sec2=0;
3884 for (Int_t i=0;i<160;i++){
3885 if (track->GetClusterPointer(i)){
3887 AliTPCTrackerPoint *trpoint =track->GetTrackPoint(i);
3888 if ( (index<p0) || x0[0]<0 ){
3889 if (trpoint->GetX()>1){
3890 clindex = track->GetClusterIndex2(i);
3892 x0[0] = trpoint->GetX();
3893 x0[1] = trpoint->GetY();
3894 x0[2] = trpoint->GetZ();
3895 sec0 = ((clindex&0xff000000)>>24)%18;
3900 if ( (index<p1) &&(trpoint->GetX()>1)){
3901 clindex = track->GetClusterIndex2(i);
3903 x1[0] = trpoint->GetX();
3904 x1[1] = trpoint->GetY();
3905 x1[2] = trpoint->GetZ();
3906 sec1 = ((clindex&0xff000000)>>24)%18;
3909 if ( (index<p2) &&(trpoint->GetX()>1)){
3910 clindex = track->GetClusterIndex2(i);
3912 x2[0] = trpoint->GetX();
3913 x2[1] = trpoint->GetY();
3914 x2[2] = trpoint->GetZ();
3915 sec2 = ((clindex&0xff000000)>>24)%18;
3922 Double_t alpha, cs,sn, xx2,yy2;
3924 alpha = (sec1-sec2)*fSectors->GetAlpha();
3925 cs = TMath::Cos(alpha);
3926 sn = TMath::Sin(alpha);
3927 xx2= x1[0]*cs-x1[1]*sn;
3928 yy2= x1[0]*sn+x1[1]*cs;
3932 alpha = (sec0-sec2)*fSectors->GetAlpha();
3933 cs = TMath::Cos(alpha);
3934 sn = TMath::Sin(alpha);
3935 xx2= x0[0]*cs-x0[1]*sn;
3936 yy2= x0[0]*sn+x0[1]*cs;
3942 Double_t x[5],c[15];
3946 x[4]=F1(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]);
3947 // if (x[4]>1) return 0;
3948 x[2]=F2(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]);
3949 x[3]=F3n(x2[0],x2[1],x0[0],x0[1],x2[2],x0[2],x[4]);
3950 //if (TMath::Abs(x[3]) > 2.2) return 0;
3951 //if (TMath::Abs(x[2]) > 1.99) return 0;
3953 Double_t sy =0.1, sz =0.1;
3955 Double_t sy1=0.02+track->GetSigmaY2(), sz1=0.02+track->GetSigmaZ2();
3956 Double_t sy2=0.01+track->GetSigmaY2(), sz2=0.01+track->GetSigmaZ2();
3957 Double_t sy3=0.01+track->GetSigmaY2();
3959 Double_t f40=(F1(x2[0],x2[1]+sy,x1[0],x1[1],x0[0],x0[1])-x[4])/sy;
3960 Double_t f42=(F1(x2[0],x2[1],x1[0],x1[1]+sy,x0[0],x0[1])-x[4])/sy;
3961 Double_t f43=(F1(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]+sy)-x[4])/sy;
3962 Double_t f20=(F2(x2[0],x2[1]+sy,x1[0],x1[1],x0[0],x0[1])-x[2])/sy;
3963 Double_t f22=(F2(x2[0],x2[1],x1[0],x1[1]+sy,x0[0],x0[1])-x[2])/sy;
3964 Double_t f23=(F2(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]+sy)-x[2])/sy;
3966 Double_t f30=(F3(x2[0],x2[1]+sy,x0[0],x0[1],x2[2],x0[2])-x[3])/sy;
3967 Double_t f31=(F3(x2[0],x2[1],x0[0],x0[1],x2[2]+sz,x0[2])-x[3])/sz;
3968 Double_t f32=(F3(x2[0],x2[1],x0[0],x0[1]+sy,x2[2],x0[2])-x[3])/sy;
3969 Double_t f34=(F3(x2[0],x2[1],x0[0],x0[1],x2[2],x0[2]+sz)-x[3])/sz;
3974 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3975 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3976 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3977 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3978 c[13]=f30*sy1*f40+f32*sy2*f42;
3979 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3981 // Int_t row1 = fSectors->GetRowNumber(x2[0]);
3982 AliTPCseed *seed=new AliTPCseed(x2[0], sec2*fSectors->GetAlpha()+fSectors->GetAlphaShift(), x, c, 0);
3983 // Double_t y0,z0,y1,z1, y2,z2;
3984 //seed->GetProlongation(x0[0],y0,z0);
3985 // seed->GetProlongation(x1[0],y1,z1);
3986 //seed->GetProlongation(x2[0],y2,z2);
3988 seed->SetLastPoint(pp2);
3989 seed->SetFirstPoint(pp2);
3996 AliTPCseed *AliTPCtrackerMI::ReSeed(const AliTPCseed *track, Float_t r0, Float_t r1, Float_t r2)
4000 //reseed using founded clusters
4002 // Find the number of clusters
4003 Int_t nclusters = 0;
4004 for (Int_t irow=0;irow<160;irow++){
4005 if (track->GetClusterIndex(irow)>0) nclusters++;
4009 ipos[0] = TMath::Max(int(r0*nclusters),0); // point 0 cluster
4010 ipos[1] = TMath::Min(int(r1*nclusters),nclusters-1); //
4011 ipos[2] = TMath::Min(int(r2*nclusters),nclusters-1); // last point
4014 Double_t xyz[3][3]={{0}};
4015 Int_t row[3]={0},sec[3]={0,0,0};
4017 // find track row position at given ratio of the length
4019 for (Int_t irow=0;irow<160;irow++){
4020 if (track->GetClusterIndex2(irow)<0) continue;
4022 for (Int_t ipoint=0;ipoint<3;ipoint++){
4023 if (index<=ipos[ipoint]) row[ipoint] = irow;
4027 //Get cluster and sector position
4028 for (Int_t ipoint=0;ipoint<3;ipoint++){
4029 Int_t clindex = track->GetClusterIndex2(row[ipoint]);
4030 AliTPCclusterMI * cl = GetClusterMI(clindex);
4033 // AliTPCclusterMI * cl = GetClusterMI(clindex);
4036 sec[ipoint] = ((clindex&0xff000000)>>24)%18;
4037 xyz[ipoint][0] = GetXrow(row[ipoint]);
4038 xyz[ipoint][1] = cl->GetY();
4039 xyz[ipoint][2] = cl->GetZ();
4043 // Calculate seed state vector and covariance matrix
4045 Double_t alpha, cs,sn, xx2,yy2;
4047 alpha = (sec[1]-sec[2])*fSectors->GetAlpha();
4048 cs = TMath::Cos(alpha);
4049 sn = TMath::Sin(alpha);
4050 xx2= xyz[1][0]*cs-xyz[1][1]*sn;
4051 yy2= xyz[1][0]*sn+xyz[1][1]*cs;
4055 alpha = (sec[0]-sec[2])*fSectors->GetAlpha();
4056 cs = TMath::Cos(alpha);
4057 sn = TMath::Sin(alpha);
4058 xx2= xyz[0][0]*cs-xyz[0][1]*sn;
4059 yy2= xyz[0][0]*sn+xyz[0][1]*cs;
4065 Double_t x[5],c[15];
4069 x[4]=F1(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
4070 x[2]=F2(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
4071 x[3]=F3n(xyz[2][0],xyz[2][1],xyz[0][0],xyz[0][1],xyz[2][2],xyz[0][2],x[4]);
4073 Double_t sy =0.1, sz =0.1;
4075 Double_t sy1=0.2, sz1=0.2;
4076 Double_t sy2=0.2, sz2=0.2;
4079 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;
4080 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;
4081 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;
4082 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;
4083 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;
4084 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;
4086 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;
4087 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;
4088 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;
4089 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;
4094 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
4095 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
4096 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
4097 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
4098 c[13]=f30*sy1*f40+f32*sy2*f42;
4099 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
4101 // Int_t row1 = fSectors->GetRowNumber(xyz[2][0]);
4102 AliTPCseed *seed=new AliTPCseed(xyz[2][0], sec[2]*fSectors->GetAlpha()+fSectors->GetAlphaShift(), x, c, 0);
4103 seed->SetLastPoint(row[2]);
4104 seed->SetFirstPoint(row[2]);
4109 AliTPCseed *AliTPCtrackerMI::ReSeed(AliTPCseed *track,Int_t r0, Bool_t forward)
4113 //reseed using founded clusters
4116 Int_t row[3]={0,0,0};
4117 Int_t sec[3]={0,0,0};
4119 // forward direction
4121 for (Int_t irow=r0;irow<160;irow++){
4122 if (track->GetClusterIndex(irow)>0){
4127 for (Int_t irow=160;irow>r0;irow--){
4128 if (track->GetClusterIndex(irow)>0){
4133 for (Int_t irow=row[2]-15;irow>row[0];irow--){
4134 if (track->GetClusterIndex(irow)>0){
4142 for (Int_t irow=0;irow<r0;irow++){
4143 if (track->GetClusterIndex(irow)>0){
4148 for (Int_t irow=r0;irow>0;irow--){
4149 if (track->GetClusterIndex(irow)>0){
4154 for (Int_t irow=row[2]-15;irow>row[0];irow--){
4155 if (track->GetClusterIndex(irow)>0){
4162 if ((row[2]-row[0])<20) return 0;
4163 if (row[1]==0) return 0;
4166 //Get cluster and sector position
4167 for (Int_t ipoint=0;ipoint<3;ipoint++){
4168 Int_t clindex = track->GetClusterIndex2(row[ipoint]);
4169 AliTPCclusterMI * cl = GetClusterMI(clindex);
4172 // AliTPCclusterMI * cl = GetClusterMI(clindex);
4175 sec[ipoint] = ((clindex&0xff000000)>>24)%18;
4176 xyz[ipoint][0] = GetXrow(row[ipoint]);
4177 AliTPCTrackerPoint * point = track->GetTrackPoint(row[ipoint]);
4178 if (point&&ipoint<2){
4180 xyz[ipoint][1] = point->GetY();
4181 xyz[ipoint][2] = point->GetZ();
4184 xyz[ipoint][1] = cl->GetY();
4185 xyz[ipoint][2] = cl->GetZ();
4192 // Calculate seed state vector and covariance matrix
4194 Double_t alpha, cs,sn, xx2,yy2;
4196 alpha = (sec[1]-sec[2])*fSectors->GetAlpha();
4197 cs = TMath::Cos(alpha);
4198 sn = TMath::Sin(alpha);
4199 xx2= xyz[1][0]*cs-xyz[1][1]*sn;
4200 yy2= xyz[1][0]*sn+xyz[1][1]*cs;
4204 alpha = (sec[0]-sec[2])*fSectors->GetAlpha();
4205 cs = TMath::Cos(alpha);
4206 sn = TMath::Sin(alpha);
4207 xx2= xyz[0][0]*cs-xyz[0][1]*sn;
4208 yy2= xyz[0][0]*sn+xyz[0][1]*cs;
4214 Double_t x[5],c[15];
4218 x[4]=F1(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
4219 x[2]=F2(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
4220 x[3]=F3n(xyz[2][0],xyz[2][1],xyz[0][0],xyz[0][1],xyz[2][2],xyz[0][2],x[4]);
4222 Double_t sy =0.1, sz =0.1;
4224 Double_t sy1=0.2, sz1=0.2;
4225 Double_t sy2=0.2, sz2=0.2;
4228 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;
4229 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;
4230 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;
4231 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;
4232 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;
4233 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;
4235 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;
4236 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;
4237 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;
4238 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;
4243 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
4244 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
4245 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
4246 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
4247 c[13]=f30*sy1*f40+f32*sy2*f42;
4248 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
4250 // Int_t row1 = fSectors->GetRowNumber(xyz[2][0]);
4251 AliTPCseed *seed=new AliTPCseed(xyz[2][0], sec[2]*fSectors->GetAlpha()+fSectors->GetAlphaShift(), x, c, 0);
4252 seed->SetLastPoint(row[2]);
4253 seed->SetFirstPoint(row[2]);
4254 for (Int_t i=row[0];i<row[2];i++){
4255 seed->SetClusterIndex(i, track->GetClusterIndex(i));
4263 void AliTPCtrackerMI::FindMultiMC(const TObjArray * array, AliESDEvent */*esd*/, Int_t iter)
4266 // find multi tracks - THIS FUNCTION IS ONLY FOR DEBUG PURPOSES
4268 // Use AliTPCReconstructor::StreamLevel()>2 if you want to tune parameters - cuts
4270 // Two reasons to have multiple find tracks
4271 // 1. Curling tracks can be find more than once
4272 // 2. Splitted tracks
4273 // a.) Multiple seeding to increase tracking efficiency - (~ 100% reached)
4274 // b.) Edge effect on the sector boundaries
4277 // Algorithm done in 2 phases - because of CPU consumption
4278 // it is n^2 algorithm - for lead-lead 20000x20000 combination are investigated
4280 // Algorihm for curling tracks sign:
4281 // 1 phase -makes a very rough fast cuts to minimize combinatorics
4282 // a.) opposite sign
4283 // b.) one of the tracks - not pointing to the primary vertex -
4284 // c.) delta tan(theta)
4286 // 2 phase - calculates DCA between tracks - time consument
4291 // General cuts - for splitted tracks and for curling tracks
4293 const Float_t kMaxdPhi = 0.2; // maximal distance in phi
4295 // Curling tracks cuts
4300 Int_t nentries = array->GetEntriesFast();
4301 AliHelix *helixes = new AliHelix[nentries];
4302 Float_t *xm = new Float_t[nentries];
4303 Float_t *dz0 = new Float_t[nentries];
4304 Float_t *dz1 = new Float_t[nentries];
4310 // Find track COG in x direction - point with best defined parameters
4312 for (Int_t i=0;i<nentries;i++){
4313 AliTPCseed* track = (AliTPCseed*)array->At(i);
4314 if (!track) continue;
4315 track->SetCircular(0);
4316 new (&helixes[i]) AliHelix(*track);
4320 track->GetDZ(GetX(),GetY(),GetZ(),GetBz(),dz);
4323 for (Int_t icl=0; icl<160; icl++){
4324 AliTPCclusterMI * cl = track->GetClusterPointer(icl);
4330 if (ncl>0) xm[i]/=Float_t(ncl);
4333 for (Int_t i0=0;i0<nentries;i0++){
4334 AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
4335 if (!track0) continue;
4336 Float_t xc0 = helixes[i0].GetHelix(6);
4337 Float_t yc0 = helixes[i0].GetHelix(7);
4338 Float_t r0 = helixes[i0].GetHelix(8);
4339 Float_t rc0 = TMath::Sqrt(xc0*xc0+yc0*yc0);
4340 Float_t fi0 = TMath::ATan2(yc0,xc0);
4342 for (Int_t i1=i0+1;i1<nentries;i1++){
4343 AliTPCseed * track1 = (AliTPCseed*)array->At(i1);
4344 if (!track1) continue;
4345 Int_t lab0=track0->GetLabel();
4346 Int_t lab1=track1->GetLabel();
4347 if (TMath::Abs(lab0)!=TMath::Abs(lab1)) continue;
4349 Float_t xc1 = helixes[i1].GetHelix(6);
4350 Float_t yc1 = helixes[i1].GetHelix(7);
4351 Float_t r1 = helixes[i1].GetHelix(8);
4352 Float_t rc1 = TMath::Sqrt(xc1*xc1+yc1*yc1);
4353 Float_t fi1 = TMath::ATan2(yc1,xc1);
4355 Float_t dfi = fi0-fi1;
4358 if (dfi>1.5*TMath::Pi()) dfi-=TMath::Pi(); // take care about edge effect
4359 if (dfi<-1.5*TMath::Pi()) dfi+=TMath::Pi(); //
4360 if (TMath::Abs(dfi)>kMaxdPhi&&helixes[i0].GetHelix(4)*helixes[i1].GetHelix(4)<0){
4362 // if short tracks with undefined sign
4363 fi1 = -TMath::ATan2(yc1,-xc1);
4366 Float_t dtheta = TMath::Abs(track0->GetTgl()-track1->GetTgl())<TMath::Abs(track0->GetTgl()+track1->GetTgl())? track0->GetTgl()-track1->GetTgl():track0->GetTgl()+track1->GetTgl();
4369 // debug stream to tune "fast cuts"
4371 Double_t dist[3]; // distance at X
4372 Double_t mdist[3]={0,0,0}; // mean distance X+-40cm
4373 track0->GetDistance(track1,0.5*(xm[i0]+xm[i1])-40.,dist,AliTracker::GetBz());
4374 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[i]);
4375 track0->GetDistance(track1,0.5*(xm[i0]+xm[i1])+40.,dist,AliTracker::GetBz());
4376 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[i]);
4377 track0->GetDistance(track1,0.5*(xm[i0]+xm[i1]),dist,AliTracker::GetBz());
4378 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[i]);
4379 for (Int_t i=0;i<3;i++) mdist[i]*=0.33333;
4383 for (Int_t icl=0; icl<160; icl++){
4384 AliTPCclusterMI * cl0 = track0->GetClusterPointer(icl);
4385 AliTPCclusterMI * cl1 = track1->GetClusterPointer(icl);
4388 if (cl0==cl1) sums++;
4392 if (AliTPCReconstructor::StreamLevel()>5) {
4393 TTreeSRedirector &cstream = *fDebugStreamer;
4398 "Tr0.="<<track0<< // seed0
4399 "Tr1.="<<track1<< // seed1
4400 "h0.="<<&helixes[i0]<<
4401 "h1.="<<&helixes[i1]<<
4403 "sum="<<sum<< //the sum of rows with cl in both
4404 "sums="<<sums<< //the sum of shared clusters
4405 "xm0="<<xm[i0]<< // the center of track
4406 "xm1="<<xm[i1]<< // the x center of track
4407 // General cut variables
4408 "dfi="<<dfi<< // distance in fi angle
4409 "dtheta="<<dtheta<< // distance int theta angle
4415 "dist0="<<dist[0]<< //distance x
4416 "dist1="<<dist[1]<< //distance y
4417 "dist2="<<dist[2]<< //distance z
4418 "mdist0="<<mdist[0]<< //distance x
4419 "mdist1="<<mdist[1]<< //distance y
4420 "mdist2="<<mdist[2]<< //distance z
4436 if (AliTPCReconstructor::StreamLevel()>1) {
4437 AliInfo("Time for curling tracks removal DEBUGGING MC");
4444 void AliTPCtrackerMI::FindSplitted(TObjArray * array, AliESDEvent */*esd*/, Int_t /*iter*/){
4446 // Find Splitted tracks and remove the one with worst quality
4447 // Corresponding debug streamer to tune selections - "Splitted2"
4449 // 0. Sort tracks according quility
4450 // 1. Propagate the tracks to the reference radius
4451 // 2. Double_t loop to select close tracks (only to speed up process)
4452 // 3. Calculate cluster overlap ratio - and remove the track if bigger than a threshold
4453 // 4. Delete temporary parameters
4455 const Double_t xref=GetXrow(63); // reference radius -IROC/OROC boundary
4457 const Double_t kCutP1=10; // delta Z cut 10 cm
4458 const Double_t kCutP2=0.15; // delta snp(fi) cut 0.15
4459 const Double_t kCutP3=0.15; // delta tgl(theta) cut 0.15
4460 const Double_t kCutAlpha=0.15; // delta alpha cut
4461 Int_t firstpoint = 0;
4462 Int_t lastpoint = 160;
4464 Int_t nentries = array->GetEntriesFast();
4465 AliExternalTrackParam *params = new AliExternalTrackParam[nentries];
4471 //0. Sort tracks according quality
4472 //1. Propagate the ext. param to reference radius
4473 Int_t nseed = array->GetEntriesFast();
4474 if (nseed<=0) return;
4475 Float_t * quality = new Float_t[nseed];
4476 Int_t * indexes = new Int_t[nseed];
4477 for (Int_t i=0; i<nseed; i++) {
4478 AliTPCseed *pt=(AliTPCseed*)array->UncheckedAt(i);
4483 pt->UpdatePoints(); //select first last max dens points
4484 Float_t * points = pt->GetPoints();
4485 if (points[3]<0.8) quality[i] =-1;
4486 quality[i] = (points[2]-points[0])+pt->GetNumberOfClusters();
4487 //prefer high momenta tracks if overlaps
4488 quality[i] *= TMath::Sqrt(TMath::Abs(pt->Pt())+0.5);
4490 AliTracker::PropagateTrackToBxByBz(&(params[i]),xref,pt->GetMass(),5.,kTRUE);
4491 AliTracker::PropagateTrackToBxByBz(&(params[i]),xref,pt->GetMass(),1.,kTRUE);
4493 TMath::Sort(nseed,quality,indexes);
4495 // 3. Loop over pair of tracks
4497 for (Int_t i0=0; i0<nseed; i0++) {
4498 Int_t index0=indexes[i0];
4499 if (!(array->UncheckedAt(index0))) continue;
4500 AliTPCseed *s1 = (AliTPCseed*)array->UncheckedAt(index0);
4501 if (!s1->IsActive()) continue;
4502 AliExternalTrackParam &par0=params[index0];
4503 for (Int_t i1=i0+1; i1<nseed; i1++) {
4504 Int_t index1=indexes[i1];
4505 if (!(array->UncheckedAt(index1))) continue;
4506 AliTPCseed *s2 = (AliTPCseed*)array->UncheckedAt(index1);
4507 if (!s2->IsActive()) continue;
4508 if (s2->GetKinkIndexes()[0]!=0)
4509 if (s2->GetKinkIndexes()[0] == -s1->GetKinkIndexes()[0]) continue;
4510 AliExternalTrackParam &par1=params[index1];
4511 if (TMath::Abs(par0.GetParameter()[3]-par1.GetParameter()[3])>kCutP3) continue;
4512 if (TMath::Abs(par0.GetParameter()[1]-par1.GetParameter()[1])>kCutP1) continue;
4513 if (TMath::Abs(par0.GetParameter()[2]-par1.GetParameter()[2])>kCutP2) continue;
4514 Double_t dAlpha= TMath::Abs(par0.GetAlpha()-par1.GetAlpha());
4515 if (dAlpha>TMath::Pi()) dAlpha-=TMath::Pi();
4516 if (TMath::Abs(dAlpha)>kCutAlpha) continue;
4521 Int_t firstShared=lastpoint, lastShared=firstpoint;
4522 Int_t firstRow=lastpoint, lastRow=firstpoint;
4524 for (Int_t i=firstpoint;i<lastpoint;i++){
4525 if (s1->GetClusterIndex2(i)>0) nall0++;
4526 if (s2->GetClusterIndex2(i)>0) nall1++;
4527 if (s1->GetClusterIndex2(i)>0 && s2->GetClusterIndex2(i)>0) {
4528 if (i<firstRow) firstRow=i;
4529 if (i>lastRow) lastRow=i;
4531 if ( (s1->GetClusterIndex2(i))==(s2->GetClusterIndex2(i)) && s1->GetClusterIndex2(i)>0) {
4532 if (i<firstShared) firstShared=i;
4533 if (i>lastShared) lastShared=i;
4537 Double_t ratio0 = Float_t(sumShared)/Float_t(TMath::Min(nall0+1,nall1+1));
4538 Double_t ratio1 = Float_t(sumShared)/Float_t(TMath::Max(nall0+1,nall1+1));
4540 if( AliTPCReconstructor::StreamLevel()>1){
4541 TTreeSRedirector &cstream = *fDebugStreamer;
4542 Int_t n0=s1->GetNumberOfClusters();
4543 Int_t n1=s2->GetNumberOfClusters();
4544 Int_t n0F=s1->GetNFoundable();
4545 Int_t n1F=s2->GetNFoundable();
4546 Int_t lab0=s1->GetLabel();
4547 Int_t lab1=s2->GetLabel();
4549 cstream<<"Splitted2"<<
4550 "iter="<<fIteration<<
4551 "lab0="<<lab0<< // MC label if exist
4552 "lab1="<<lab1<< // MC label if exist
4555 "ratio0="<<ratio0<< // shared ratio
4556 "ratio1="<<ratio1<< // shared ratio
4557 "p0.="<<&par0<< // track parameters
4559 "s0.="<<s1<< // full seed
4561 "n0="<<n0<< // number of clusters track 0
4562 "n1="<<n1<< // number of clusters track 1
4563 "nall0="<<nall0<< // number of clusters track 0
4564 "nall1="<<nall1<< // number of clusters track 1
4565 "n0F="<<n0F<< // number of findable
4566 "n1F="<<n1F<< // number of findable
4567 "shared="<<sumShared<< // number of shared clusters
4568 "firstS="<<firstShared<< // first and the last shared row
4569 "lastS="<<lastShared<<
4570 "firstRow="<<firstRow<< // first and the last row with cluster
4571 "lastRow="<<lastRow<< //
4575 // remove track with lower quality
4577 if (ratio0>AliTPCReconstructor::GetRecoParam()->GetCutSharedClusters(0) ||
4578 ratio1>AliTPCReconstructor::GetRecoParam()->GetCutSharedClusters(1)){
4582 delete array->RemoveAt(index1);
4587 // 4. Delete temporary array
4597 void AliTPCtrackerMI::FindCurling(const TObjArray * array, AliESDEvent */*esd*/, Int_t iter)
4600 // find Curling tracks
4601 // Use AliTPCReconstructor::StreamLevel()>1 if you want to tune parameters - cuts
4604 // Algorithm done in 2 phases - because of CPU consumption
4605 // it is n^2 algorithm - for lead-lead 20000x20000 combination are investigated
4606 // see detal in MC part what can be used to cut
4610 const Float_t kMaxC = 400; // maximal curvature to of the track
4611 const Float_t kMaxdTheta = 0.15; // maximal distance in theta
4612 const Float_t kMaxdPhi = 0.15; // maximal distance in phi
4613 const Float_t kPtRatio = 0.3; // ratio between pt
4614 const Float_t kMinDCAR = 2.; // distance to the primary vertex in r - see cpipe cut
4617 // Curling tracks cuts
4620 const Float_t kMaxDeltaRMax = 40; // distance in outer radius
4621 const Float_t kMaxDeltaRMin = 5.; // distance in lower radius - see cpipe cut
4622 const Float_t kMinAngle = 2.9; // angle between tracks
4623 const Float_t kMaxDist = 5; // biggest distance
4625 // The cuts can be tuned using the "MC information stored in Multi tree ==> see FindMultiMC
4628 TCut csign("csign","Tr0.fP[4]*Tr1.fP[4]<0"); //opposite sign
4629 TCut cmax("cmax","abs(Tr0.GetC())>1/400");
4630 TCut cda("cda","sqrt(dtheta^2+dfi^2)<0.15");
4631 TCut ccratio("ccratio","abs((Tr0.fP[4]+Tr1.fP[4])/(abs(Tr0.fP[4])+abs(Tr1.fP[4])))<0.3");
4632 TCut cpipe("cpipe", "min(abs(r0-rc0),abs(r1-rc1))>5");
4634 TCut cdrmax("cdrmax","abs(abs(rc0+r0)-abs(rc1+r1))<40")
4635 TCut cdrmin("cdrmin","abs(abs(rc0+r0)-abs(rc1+r1))<10")
4637 Multi->Draw("dfi","iter==0"+csign+cmax+cda+ccratio); ~94% of curling tracks fulfill
4638 Multi->Draw("min(abs(r0-rc0),abs(r1-rc1))","iter==0&&abs(lab1)==abs(lab0)"+csign+cmax+cda+ccratio+cpipe+cdrmin+cdrmax); //80%
4640 Curling2->Draw("dfi","iter==0&&abs(lab0)==abs(lab1)"+csign+cmax+cdtheta+cdfi+ccratio)
4646 Int_t nentries = array->GetEntriesFast();
4647 AliHelix *helixes = new AliHelix[nentries];
4648 for (Int_t i=0;i<nentries;i++){
4649 AliTPCseed* track = (AliTPCseed*)array->At(i);
4650 if (!track) continue;
4651 track->SetCircular(0);
4652 new (&helixes[i]) AliHelix(*track);
4658 Double_t phase[2][2]={{0,0},{0,0}},radius[2]={0,0};
4664 for (Int_t i0=0;i0<nentries;i0++){
4665 AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
4666 if (!track0) continue;
4667 if (TMath::Abs(track0->GetC())<1/kMaxC) continue;
4668 Float_t xc0 = helixes[i0].GetHelix(6);
4669 Float_t yc0 = helixes[i0].GetHelix(7);
4670 Float_t r0 = helixes[i0].GetHelix(8);
4671 Float_t rc0 = TMath::Sqrt(xc0*xc0+yc0*yc0);
4672 Float_t fi0 = TMath::ATan2(yc0,xc0);
4674 for (Int_t i1=i0+1;i1<nentries;i1++){
4675 AliTPCseed * track1 = (AliTPCseed*)array->At(i1);
4676 if (!track1) continue;
4677 if (TMath::Abs(track1->GetC())<1/kMaxC) continue;
4678 Float_t xc1 = helixes[i1].GetHelix(6);
4679 Float_t yc1 = helixes[i1].GetHelix(7);
4680 Float_t r1 = helixes[i1].GetHelix(8);
4681 Float_t rc1 = TMath::Sqrt(xc1*xc1+yc1*yc1);
4682 Float_t fi1 = TMath::ATan2(yc1,xc1);
4684 Float_t dfi = fi0-fi1;
4687 if (dfi>1.5*TMath::Pi()) dfi-=TMath::Pi(); // take care about edge effect
4688 if (dfi<-1.5*TMath::Pi()) dfi+=TMath::Pi(); //
4689 Float_t dtheta = TMath::Abs(track0->GetTgl()-track1->GetTgl())<TMath::Abs(track0->GetTgl()+track1->GetTgl())? track0->GetTgl()-track1->GetTgl():track0->GetTgl()+track1->GetTgl();
4693 if (track0->GetBConstrain()&&track1->GetBConstrain()) continue; // not constrained
4694 if (track1->GetSigned1Pt()*track0->GetSigned1Pt()>0) continue; // not the same sign
4695 if ( TMath::Abs(track1->GetTgl()+track0->GetTgl())>kMaxdTheta) continue; //distance in the Theta
4696 if ( TMath::Abs(dfi)>kMaxdPhi) continue; //distance in phi
4697 if ( TMath::Sqrt(dfi*dfi+dtheta*dtheta)>kMaxdPhi) continue; //common angular offset
4699 Float_t pt0 = track0->GetSignedPt();
4700 Float_t pt1 = track1->GetSignedPt();
4701 if ((TMath::Abs(pt0+pt1)/(TMath::Abs(pt0)+TMath::Abs(pt1)))>kPtRatio) continue;
4702 if ((iter==1) && TMath::Abs(TMath::Abs(rc0+r0)-TMath::Abs(rc1+r1))>kMaxDeltaRMax) continue;
4703 if ((iter!=1) &&TMath::Abs(TMath::Abs(rc0-r0)-TMath::Abs(rc1-r1))>kMaxDeltaRMin) continue;
4704 if (TMath::Min(TMath::Abs(rc0-r0),TMath::Abs(rc1-r1))<kMinDCAR) continue;
4707 // Now find closest approach
4711 Int_t npoints = helixes[i0].GetRPHIintersections(helixes[i1], phase, radius,10);
4712 if (npoints==0) continue;
4713 helixes[i0].GetClosestPhases(helixes[i1], phase);
4717 Double_t hangles[3];
4718 helixes[i0].Evaluate(phase[0][0],xyz0);
4719 helixes[i1].Evaluate(phase[0][1],xyz1);
4721 helixes[i0].GetAngle(phase[0][0],helixes[i1],phase[0][1],hangles);
4722 Double_t deltah[2],deltabest;
4723 if (TMath::Abs(hangles[2])<kMinAngle) continue;
4727 helixes[i0].ParabolicDCA(helixes[i1],phase[0][0],phase[0][1],radius[0],deltah[0],2);
4729 helixes[i0].ParabolicDCA(helixes[i1],phase[1][0],phase[1][1],radius[1],deltah[1],2);
4730 if (deltah[1]<deltah[0]) ibest=1;
4732 deltabest = TMath::Sqrt(deltah[ibest]);
4733 helixes[i0].Evaluate(phase[ibest][0],xyz0);
4734 helixes[i1].Evaluate(phase[ibest][1],xyz1);
4735 helixes[i0].GetAngle(phase[ibest][0],helixes[i1],phase[ibest][1],hangles);
4736 Double_t radiusbest = TMath::Sqrt(radius[ibest]);
4738 if (deltabest>kMaxDist) continue;
4739 // if (mindcar+mindcaz<40 && (TMath::Abs(hangles[2])<kMinAngle ||deltabest>3)) continue;
4740 Bool_t sign =kFALSE;
4741 if (hangles[2]>kMinAngle) sign =kTRUE;
4744 // circular[i0] = kTRUE;
4745 // circular[i1] = kTRUE;
4746 if (track0->OneOverPt()<track1->OneOverPt()){
4747 track0->SetCircular(track0->GetCircular()+1);
4748 track1->SetCircular(track1->GetCircular()+2);
4751 track1->SetCircular(track1->GetCircular()+1);
4752 track0->SetCircular(track0->GetCircular()+2);
4755 if (AliTPCReconstructor::StreamLevel()>2){
4757 //debug stream to tune "fine" cuts
4758 Int_t lab0=track0->GetLabel();
4759 Int_t lab1=track1->GetLabel();
4760 TTreeSRedirector &cstream = *fDebugStreamer;
4761 cstream<<"Curling2"<<
4777 "npoints="<<npoints<<
4778 "hangles0="<<hangles[0]<<
4779 "hangles1="<<hangles[1]<<
4780 "hangles2="<<hangles[2]<<
4783 "radius="<<radiusbest<<
4784 "deltabest="<<deltabest<<
4785 "phase0="<<phase[ibest][0]<<
4786 "phase1="<<phase[ibest][1]<<
4794 if (AliTPCReconstructor::StreamLevel()>1) {
4795 AliInfo("Time for curling tracks removal");
4804 void AliTPCtrackerMI::FindKinks(TObjArray * array, AliESDEvent *esd)
4811 TObjArray *kinks= new TObjArray(10000);
4812 // TObjArray *v0s= new TObjArray(10000);
4813 Int_t nentries = array->GetEntriesFast();
4814 AliHelix *helixes = new AliHelix[nentries];
4815 Int_t *sign = new Int_t[nentries];
4816 Int_t *nclusters = new Int_t[nentries];
4817 Float_t *alpha = new Float_t[nentries];
4818 AliKink *kink = new AliKink();
4819 Int_t * usage = new Int_t[nentries];
4820 Float_t *zm = new Float_t[nentries];
4821 Float_t *z0 = new Float_t[nentries];
4822 Float_t *fim = new Float_t[nentries];
4823 Float_t *shared = new Float_t[nentries];
4824 Bool_t *circular = new Bool_t[nentries];
4825 Float_t *dca = new Float_t[nentries];
4826 //const AliESDVertex * primvertex = esd->GetVertex();
4828 // nentries = array->GetEntriesFast();
4833 for (Int_t i=0;i<nentries;i++){
4836 AliTPCseed* track = (AliTPCseed*)array->At(i);
4837 if (!track) continue;
4838 track->SetCircular(0);
4840 track->UpdatePoints();
4841 if (( track->GetPoints()[2]- track->GetPoints()[0])>5 && track->GetPoints()[3]>0.8){
4843 nclusters[i]=track->GetNumberOfClusters();
4844 alpha[i] = track->GetAlpha();
4845 new (&helixes[i]) AliHelix(*track);
4847 helixes[i].Evaluate(0,xyz);
4848 sign[i] = (track->GetC()>0) ? -1:1;
4851 if (track->GetProlongation(x,y,z)){
4853 fim[i] = alpha[i]+TMath::ATan2(y,x);
4856 zm[i] = track->GetZ();
4860 circular[i]= kFALSE;
4861 if (track->GetProlongation(0,y,z)) z0[i] = z;
4862 dca[i] = track->GetD(0,0);
4868 Int_t ncandidates =0;
4871 Double_t phase[2][2]={{0,0},{0,0}},radius[2]={0,0};
4874 // Find circling track
4876 for (Int_t i0=0;i0<nentries;i0++){
4877 AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
4878 if (!track0) continue;
4879 if (track0->GetNumberOfClusters()<40) continue;
4880 if (TMath::Abs(1./track0->GetC())>200) continue;
4881 for (Int_t i1=i0+1;i1<nentries;i1++){
4882 AliTPCseed * track1 = (AliTPCseed*)array->At(i1);
4883 if (!track1) continue;
4884 if (track1->GetNumberOfClusters()<40) continue;
4885 if ( TMath::Abs(track1->GetTgl()+track0->GetTgl())>0.1) continue;
4886 if (track0->GetBConstrain()&&track1->GetBConstrain()) continue;
4887 if (TMath::Abs(1./track1->GetC())>200) continue;
4888 if (track1->GetSigned1Pt()*track0->GetSigned1Pt()>0) continue;
4889 if (track1->GetTgl()*track0->GetTgl()>0) continue;
4890 if (TMath::Max(TMath::Abs(1./track0->GetC()),TMath::Abs(1./track1->GetC()))>190) continue;
4891 if (track0->GetBConstrain()&&track1->OneOverPt()<track0->OneOverPt()) continue; //returning - lower momenta
4892 if (track1->GetBConstrain()&&track0->OneOverPt()<track1->OneOverPt()) continue; //returning - lower momenta
4894 Float_t mindcar = TMath::Min(TMath::Abs(dca[i0]),TMath::Abs(dca[i1]));
4895 if (mindcar<5) continue;
4896 Float_t mindcaz = TMath::Min(TMath::Abs(z0[i0]-GetZ()),TMath::Abs(z0[i1]-GetZ()));
4897 if (mindcaz<5) continue;
4898 if (mindcar+mindcaz<20) continue;
4901 Float_t xc0 = helixes[i0].GetHelix(6);
4902 Float_t yc0 = helixes[i0].GetHelix(7);
4903 Float_t r0 = helixes[i0].GetHelix(8);
4904 Float_t xc1 = helixes[i1].GetHelix(6);
4905 Float_t yc1 = helixes[i1].GetHelix(7);
4906 Float_t r1 = helixes[i1].GetHelix(8);
4908 Float_t rmean = (r0+r1)*0.5;
4909 Float_t delta =TMath::Sqrt((xc1-xc0)*(xc1-xc0)+(yc1-yc0)*(yc1-yc0));
4910 //if (delta>30) continue;
4911 if (delta>rmean*0.25) continue;
4912 if (TMath::Abs(r0-r1)/rmean>0.3) continue;
4914 Int_t npoints = helixes[i0].GetRPHIintersections(helixes[i1], phase, radius,10);
4915 if (npoints==0) continue;
4916 helixes[i0].GetClosestPhases(helixes[i1], phase);
4920 Double_t hangles[3];
4921 helixes[i0].Evaluate(phase[0][0],xyz0);
4922 helixes[i1].Evaluate(phase[0][1],xyz1);
4924 helixes[i0].GetAngle(phase[0][0],helixes[i1],phase[0][1],hangles);
4925 Double_t deltah[2],deltabest;
4926 if (hangles[2]<2.8) continue;
4929 helixes[i0].ParabolicDCA(helixes[i1],phase[0][0],phase[0][1],radius[0],deltah[0],2);
4931 helixes[i0].ParabolicDCA(helixes[i1],phase[1][0],phase[1][1],radius[1],deltah[1],2);
4932 if (deltah[1]<deltah[0]) ibest=1;
4934 deltabest = TMath::Sqrt(deltah[ibest]);
4935 helixes[i0].Evaluate(phase[ibest][0],xyz0);
4936 helixes[i1].Evaluate(phase[ibest][1],xyz1);
4937 helixes[i0].GetAngle(phase[ibest][0],helixes[i1],phase[ibest][1],hangles);
4938 Double_t radiusbest = TMath::Sqrt(radius[ibest]);
4940 if (deltabest>6) continue;
4941 if (mindcar+mindcaz<40 && (hangles[2]<3.12||deltabest>3)) continue;
4942 Bool_t lsign =kFALSE;
4943 if (hangles[2]>3.06) lsign =kTRUE;
4946 circular[i0] = kTRUE;
4947 circular[i1] = kTRUE;
4948 if (track0->OneOverPt()<track1->OneOverPt()){
4949 track0->SetCircular(track0->GetCircular()+1);
4950 track1->SetCircular(track1->GetCircular()+2);
4953 track1->SetCircular(track1->GetCircular()+1);
4954 track0->SetCircular(track0->GetCircular()+2);
4957 if (lsign&&AliTPCReconstructor::StreamLevel()>1){
4959 Int_t lab0=track0->GetLabel();
4960 Int_t lab1=track1->GetLabel();
4961 TTreeSRedirector &cstream = *fDebugStreamer;
4962 cstream<<"Curling"<<
4969 "mindcar="<<mindcar<<
4970 "mindcaz="<<mindcaz<<
4973 "npoints="<<npoints<<
4974 "hangles0="<<hangles[0]<<
4975 "hangles2="<<hangles[2]<<
4980 "radius="<<radiusbest<<
4981 "deltabest="<<deltabest<<
4982 "phase0="<<phase[ibest][0]<<
4983 "phase1="<<phase[ibest][1]<<
4993 for (Int_t i =0;i<nentries;i++){
4994 if (sign[i]==0) continue;
4995 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5002 Double_t cradius0 = 40*40;
5003 Double_t cradius1 = 270*270;
5006 Double_t cdist3=0.55;
5007 for (Int_t j =i+1;j<nentries;j++){
5009 if (sign[j]*sign[i]<1) continue;
5010 if ( (nclusters[i]+nclusters[j])>200) continue;
5011 if ( (nclusters[i]+nclusters[j])<80) continue;
5012 if ( TMath::Abs(zm[i]-zm[j])>60.) continue;
5013 if ( TMath::Abs(fim[i]-fim[j])>0.6 && TMath::Abs(fim[i]-fim[j])<5.7 ) continue;
5014 //AliTPCseed * track1 = (AliTPCseed*)array->At(j); Double_t phase[2][2],radius[2];
5015 Int_t npoints = helixes[i].GetRPHIintersections(helixes[j], phase, radius,20);
5016 if (npoints<1) continue;
5019 if (radius[0]<cradius0||radius[0]>cradius1) continue;
5022 if ( (radius[0]<cradius0||radius[0]>cradius1) && (radius[1]<cradius0||radius[1]>cradius1) ) continue;
5025 Double_t delta1=10000,delta2=10000;
5026 // cuts on the intersection radius
5027 helixes[i].LinearDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta1);
5028 if (radius[0]<20&&delta1<1) continue; //intersection at vertex
5029 if (radius[0]<10&&delta1<3) continue; //intersection at vertex
5031 helixes[i].LinearDCA(helixes[j],phase[1][0],phase[1][1],radius[1],delta2);
5032 if (radius[1]<20&&delta2<1) continue; //intersection at vertex
5033 if (radius[1]<10&&delta2<3) continue; //intersection at vertex
5036 Double_t distance1 = TMath::Min(delta1,delta2);
5037 if (distance1>cdist1) continue; // cut on DCA linear approximation
5039 npoints = helixes[i].GetRPHIintersections(helixes[j], phase, radius,20);
5040 helixes[i].ParabolicDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta1);
5041 if (radius[0]<20&&delta1<1) continue; //intersection at vertex
5042 if (radius[0]<10&&delta1<3) continue; //intersection at vertex
5045 helixes[i].ParabolicDCA(helixes[j],phase[1][0],phase[1][1],radius[1],delta2);
5046 if (radius[1]<20&&delta2<1) continue; //intersection at vertex
5047 if (radius[1]<10&&delta2<3) continue; //intersection at vertex
5049 distance1 = TMath::Min(delta1,delta2);
5052 rkink = TMath::Sqrt(radius[0]);
5055 rkink = TMath::Sqrt(radius[1]);
5057 if (distance1>cdist2) continue;
5060 AliTPCseed * track1 = (AliTPCseed*)array->At(j);
5063 Int_t row0 = GetRowNumber(rkink);
5064 if (row0<10) continue;
5065 if (row0>150) continue;
5068 Float_t dens00=-1,dens01=-1;
5069 Float_t dens10=-1,dens11=-1;
5071 Int_t found,foundable,ishared;
5072 track0->GetClusterStatistic(0,row0-5, found, foundable,ishared,kFALSE);
5073 if (foundable>5) dens00 = Float_t(found)/Float_t(foundable);
5074 track0->GetClusterStatistic(row0+5,155, found, foundable,ishared,kFALSE);
5075 if (foundable>5) dens01 = Float_t(found)/Float_t(foundable);
5077 track1->GetClusterStatistic(0,row0-5, found, foundable,ishared,kFALSE);
5078 if (foundable>10) dens10 = Float_t(found)/Float_t(foundable);
5079 track1->GetClusterStatistic(row0+5,155, found, foundable,ishared,kFALSE);
5080 if (foundable>10) dens11 = Float_t(found)/Float_t(foundable);
5082 if (dens00<dens10 && dens01<dens11) continue;
5083 if (dens00>dens10 && dens01>dens11) continue;
5084 if (TMath::Max(dens00,dens10)<0.1) continue;
5085 if (TMath::Max(dens01,dens11)<0.3) continue;
5087 if (TMath::Min(dens00,dens10)>0.6) continue;
5088 if (TMath::Min(dens01,dens11)>0.6) continue;
5091 AliTPCseed * ktrack0, *ktrack1;
5100 if (TMath::Abs(ktrack0->GetC())>5) continue; // cut on the curvature for mother particle
5101 AliExternalTrackParam paramm(*ktrack0);
5102 AliExternalTrackParam paramd(*ktrack1);
5103 if (row0>60&&ktrack1->GetReference().GetX()>90.)new (¶md) AliExternalTrackParam(ktrack1->GetReference());
5106 kink->SetMother(paramm);
5107 kink->SetDaughter(paramd);
5110 Float_t x[3] = { kink->GetPosition()[0],kink->GetPosition()[1],kink->GetPosition()[2]};
5112 fkParam->Transform0to1(x,index);
5113 fkParam->Transform1to2(x,index);
5114 row0 = GetRowNumber(x[0]);
5116 if (kink->GetR()<100) continue;
5117 if (kink->GetR()>240) continue;
5118 if (kink->GetPosition()[2]/kink->GetR()>AliTPCReconstructor::GetCtgRange()) continue; //out of fiducial volume
5119 if (kink->GetDistance()>cdist3) continue;
5120 Float_t dird = kink->GetDaughterP()[0]*kink->GetPosition()[0]+kink->GetDaughterP()[1]*kink->GetPosition()[1]; // rough direction estimate
5121 if (dird<0) continue;
5123 Float_t dirm = kink->GetMotherP()[0]*kink->GetPosition()[0]+kink->GetMotherP()[1]*kink->GetPosition()[1]; // rough direction estimate
5124 if (dirm<0) continue;
5125 Float_t mpt = TMath::Sqrt(kink->GetMotherP()[0]*kink->GetMotherP()[0]+kink->GetMotherP()[1]*kink->GetMotherP()[1]);
5126 if (mpt<0.2) continue;
5129 //for high momenta momentum not defined well in first iteration
5130 Double_t qt = TMath::Sin(kink->GetAngle(2))*ktrack1->GetP();
5131 if (qt>0.35) continue;
5134 kink->SetLabel(CookLabel(ktrack0,0.4,0,row0),0);
5135 kink->SetLabel(CookLabel(ktrack1,0.4,row0,160),1);
5137 kink->SetTPCDensity(dens00,0,0);
5138 kink->SetTPCDensity(dens01,0,1);
5139 kink->SetTPCDensity(dens10,1,0);
5140 kink->SetTPCDensity(dens11,1,1);
5141 kink->SetIndex(i,0);
5142 kink->SetIndex(j,1);
5145 kink->SetTPCDensity(dens10,0,0);
5146 kink->SetTPCDensity(dens11,0,1);
5147 kink->SetTPCDensity(dens00,1,0);
5148 kink->SetTPCDensity(dens01,1,1);
5149 kink->SetIndex(j,0);
5150 kink->SetIndex(i,1);
5153 if (mpt<1||kink->GetAngle(2)>0.1){
5154 // angle and densities not defined yet
5155 if (kink->GetTPCDensityFactor()<0.8) continue;
5156 if ((2-kink->GetTPCDensityFactor())*kink->GetDistance() >0.25) continue;
5157 if (kink->GetAngle(2)*ktrack0->GetP()<0.003) continue; //too small angle
5158 if (kink->GetAngle(2)>0.2&&kink->GetTPCDensityFactor()<1.15) continue;
5159 if (kink->GetAngle(2)>0.2&&kink->GetTPCDensity(0,1)>0.05) continue;
5161 Float_t criticalangle = track0->GetSigmaSnp2()+track0->GetSigmaTgl2();
5162 criticalangle+= track1->GetSigmaSnp2()+track1->GetSigmaTgl2();
5163 criticalangle= 3*TMath::Sqrt(criticalangle);
5164 if (criticalangle>0.02) criticalangle=0.02;
5165 if (kink->GetAngle(2)<criticalangle) continue;
5168 Int_t drow = Int_t(2.+0.5/(0.05+kink->GetAngle(2))); // overlap region defined
5169 Float_t shapesum =0;
5171 for ( Int_t row = row0-drow; row<row0+drow;row++){
5172 if (row<0) continue;
5173 if (row>155) continue;
5174 if (ktrack0->GetClusterPointer(row)){
5175 AliTPCTrackerPoint *point =ktrack0->GetTrackPoint(row);
5176 shapesum+=point->GetSigmaY()+point->GetSigmaZ();
5179 if (ktrack1->GetClusterPointer(row)){
5180 AliTPCTrackerPoint *point =ktrack1->GetTrackPoint(row);
5181 shapesum+=point->GetSigmaY()+point->GetSigmaZ();
5186 kink->SetShapeFactor(-1.);
5189 kink->SetShapeFactor(shapesum/sum);
5191 // esd->AddKink(kink);
5193 // kink->SetMother(paramm);
5194 //kink->SetDaughter(paramd);
5196 Double_t chi2P2 = paramm.GetParameter()[2]-paramd.GetParameter()[2];
5198 chi2P2/=paramm.GetCovariance()[5]+paramd.GetCovariance()[5];
5199 Double_t chi2P3 = paramm.GetParameter()[3]-paramd.GetParameter()[3];
5201 chi2P3/=paramm.GetCovariance()[9]+paramd.GetCovariance()[9];
5203 if (AliTPCReconstructor::StreamLevel()>1) {
5204 (*fDebugStreamer)<<"kinkLpt"<<
5212 if ( chi2P2+chi2P3<AliTPCReconstructor::GetRecoParam()->GetKinkAngleCutChi2(0)){
5216 kinks->AddLast(kink);
5222 // sort the kinks according quality - and refit them towards vertex
5224 Int_t nkinks = kinks->GetEntriesFast();
5225 Float_t *quality = new Float_t[nkinks];
5226 Int_t *indexes = new Int_t[nkinks];
5227 AliTPCseed *mothers = new AliTPCseed[nkinks];
5228 AliTPCseed *daughters = new AliTPCseed[nkinks];
5231 for (Int_t i=0;i<nkinks;i++){
5233 AliKink *kinkl = (AliKink*)kinks->At(i);
5235 // refit kinks towards vertex
5237 Int_t index0 = kinkl->GetIndex(0);
5238 Int_t index1 = kinkl->GetIndex(1);
5239 AliTPCseed * ktrack0 = (AliTPCseed*)array->At(index0);
5240 AliTPCseed * ktrack1 = (AliTPCseed*)array->At(index1);
5242 Int_t sumn=ktrack0->GetNumberOfClusters()+ktrack1->GetNumberOfClusters();
5244 // Refit Kink under if too small angle
5246 if (kinkl->GetAngle(2)<0.05){
5247 kinkl->SetTPCRow0(GetRowNumber(kinkl->GetR()));
5248 Int_t row0 = kinkl->GetTPCRow0();
5249 Int_t drow = Int_t(2.+0.5/(0.05+kinkl->GetAngle(2)));
5252 Int_t last = row0-drow;
5253 if (last<40) last=40;
5254 if (last<ktrack0->GetFirstPoint()+25) last = ktrack0->GetFirstPoint()+25;
5255 AliTPCseed* seed0 = ReSeed(ktrack0,last,kFALSE);
5258 Int_t first = row0+drow;
5259 if (first>130) first=130;
5260 if (first>ktrack1->GetLastPoint()-25) first = TMath::Max(ktrack1->GetLastPoint()-25,30);
5261 AliTPCseed* seed1 = ReSeed(ktrack1,first,kTRUE);
5263 if (seed0 && seed1){
5264 kinkl->SetStatus(1,8);
5265 if (RefitKink(*seed0,*seed1,*kinkl)) kinkl->SetStatus(1,9);
5266 row0 = GetRowNumber(kinkl->GetR());
5267 sumn = seed0->GetNumberOfClusters()+seed1->GetNumberOfClusters();
5268 mothers[i] = *seed0;
5269 daughters[i] = *seed1;
5272 delete kinks->RemoveAt(i);
5273 if (seed0) delete seed0;
5274 if (seed1) delete seed1;
5277 if (kinkl->GetDistance()>0.5 || kinkl->GetR()<110 || kinkl->GetR()>240) {
5278 delete kinks->RemoveAt(i);
5279 if (seed0) delete seed0;
5280 if (seed1) delete seed1;
5288 if (kinkl) quality[i] = 160*((0.1+kinkl->GetDistance())*(2.-kinkl->GetTPCDensityFactor()))/(sumn+40.); //the longest -clossest will win
5290 TMath::Sort(nkinks,quality,indexes,kFALSE);
5292 //remove double find kinks
5294 for (Int_t ikink0=1;ikink0<nkinks;ikink0++){
5295 AliKink * kink0 = (AliKink*) kinks->At(indexes[ikink0]);
5296 if (!kink0) continue;
5298 for (Int_t ikink1=0;ikink1<ikink0;ikink1++){
5299 kink0 = (AliKink*) kinks->At(indexes[ikink0]);
5300 if (!kink0) continue;
5301 AliKink * kink1 = (AliKink*) kinks->At(indexes[ikink1]);
5302 if (!kink1) continue;
5303 // if not close kink continue
5304 if (TMath::Abs(kink1->GetPosition()[2]-kink0->GetPosition()[2])>10) continue;
5305 if (TMath::Abs(kink1->GetPosition()[1]-kink0->GetPosition()[1])>10) continue;
5306 if (TMath::Abs(kink1->GetPosition()[0]-kink0->GetPosition()[0])>10) continue;
5308 AliTPCseed &mother0 = mothers[indexes[ikink0]];
5309 AliTPCseed &daughter0 = daughters[indexes[ikink0]];
5310 AliTPCseed &mother1 = mothers[indexes[ikink1]];
5311 AliTPCseed &daughter1 = daughters[indexes[ikink1]];
5312 Int_t row0 = (kink0->GetTPCRow0()+kink1->GetTPCRow0())/2;
5321 for (Int_t i=0;i<row0;i++){
5322 if (mother0.GetClusterIndex(i)>0 && mother1.GetClusterIndex(i)>0){
5325 if (mother0.GetClusterIndex(i)==mother1.GetClusterIndex(i)){
5332 for (Int_t i=row0;i<158;i++){
5333 if (daughter0.GetClusterIndex(i)>0 && daughter0.GetClusterIndex(i)>0){
5336 if (mother0.GetClusterIndex(i)==mother1.GetClusterIndex(i)){
5342 Float_t ratio = Float_t(same+1)/Float_t(both+1);
5343 Float_t ratiom = Float_t(samem+1)/Float_t(bothm+1);
5344 Float_t ratiod = Float_t(samed+1)/Float_t(bothd+1);
5345 if (ratio>0.3 && ratiom>0.5 &&ratiod>0.5) {
5346 Int_t sum0 = mother0.GetNumberOfClusters()+daughter0.GetNumberOfClusters();
5347 Int_t sum1 = mother1.GetNumberOfClusters()+daughter1.GetNumberOfClusters();
5349 shared[kink0->GetIndex(0)]= kTRUE;
5350 shared[kink0->GetIndex(1)]= kTRUE;
5351 delete kinks->RemoveAt(indexes[ikink0]);
5355 shared[kink1->GetIndex(0)]= kTRUE;
5356 shared[kink1->GetIndex(1)]= kTRUE;
5357 delete kinks->RemoveAt(indexes[ikink1]);
5364 for (Int_t i=0;i<nkinks;i++){
5365 AliKink * kinkl = (AliKink*) kinks->At(indexes[i]);
5366 if (!kinkl) continue;
5367 kinkl->SetTPCRow0(GetRowNumber(kinkl->GetR()));
5368 Int_t index0 = kinkl->GetIndex(0);
5369 Int_t index1 = kinkl->GetIndex(1);
5370 if (circular[index0]||(circular[index1]&&kinkl->GetDistance()>0.2)) continue;
5371 kinkl->SetMultiple(usage[index0],0);
5372 kinkl->SetMultiple(usage[index1],1);
5373 if (kinkl->GetMultiple()[0]+kinkl->GetMultiple()[1]>2) continue;
5374 if (kinkl->GetMultiple()[0]+kinkl->GetMultiple()[1]>0 && quality[indexes[i]]>0.2) continue;
5375 if (kinkl->GetMultiple()[0]+kinkl->GetMultiple()[1]>0 && kinkl->GetDistance()>0.2) continue;
5376 if (circular[index0]||(circular[index1]&&kinkl->GetDistance()>0.1)) continue;
5378 AliTPCseed * ktrack0 = (AliTPCseed*)array->At(index0);
5379 AliTPCseed * ktrack1 = (AliTPCseed*)array->At(index1);
5380 if (!ktrack0 || !ktrack1) continue;
5381 Int_t index = esd->AddKink(kinkl);
5384 if ( ktrack0->GetKinkIndex(0)==0 && ktrack1->GetKinkIndex(0)==0) { //best kink
5385 if (mothers[indexes[i]].GetNumberOfClusters()>20 && daughters[indexes[i]].GetNumberOfClusters()>20 && (mothers[indexes[i]].GetNumberOfClusters()+daughters[indexes[i]].GetNumberOfClusters())>100){
5386 *ktrack0 = mothers[indexes[i]];
5387 *ktrack1 = daughters[indexes[i]];
5391 ktrack0->SetKinkIndex(usage[index0],-(index+1));
5392 ktrack1->SetKinkIndex(usage[index1], (index+1));
5397 // Remove tracks corresponding to shared kink's
5399 for (Int_t i=0;i<nentries;i++){
5400 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5401 if (!track0) continue;
5402 if (track0->GetKinkIndex(0)!=0) continue;
5403 if (shared[i]) delete array->RemoveAt(i);
5408 RemoveUsed2(array,0.5,0.4,30);
5410 for (Int_t i=0;i<nentries;i++){
5411 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5412 if (!track0) continue;
5413 track0->CookdEdx(0.02,0.6);
5417 for (Int_t i=0;i<nentries;i++){
5418 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5419 if (!track0) continue;
5420 if (track0->Pt()<1.4) continue;
5421 //remove double high momenta tracks - overlapped with kink candidates
5424 for (Int_t icl=track0->GetFirstPoint();icl<track0->GetLastPoint(); icl++){
5425 if (track0->GetClusterPointer(icl)!=0){
5427 if (track0->GetClusterPointer(icl)->IsUsed(10)) ishared++;
5430 if (Float_t(ishared+1)/Float_t(all+1)>0.5) {
5431 delete array->RemoveAt(i);
5435 if (track0->GetKinkIndex(0)!=0) continue;
5436 if (track0->GetNumberOfClusters()<80) continue;
5438 AliTPCseed *pmother = new AliTPCseed();
5439 AliTPCseed *pdaughter = new AliTPCseed();
5440 AliKink *pkink = new AliKink;
5442 AliTPCseed & mother = *pmother;
5443 AliTPCseed & daughter = *pdaughter;
5444 AliKink & kinkl = *pkink;
5445 if (CheckKinkPoint(track0,mother,daughter, kinkl)){
5446 if (mother.GetNumberOfClusters()<30||daughter.GetNumberOfClusters()<20) {
5450 continue; //too short tracks
5452 if (mother.Pt()<1.4) {
5458 Int_t row0= kinkl.GetTPCRow0();
5459 if (kinkl.GetDistance()>0.5 || kinkl.GetR()<110. || kinkl.GetR()>240.) {
5466 Int_t index = esd->AddKink(&kinkl);
5467 mother.SetKinkIndex(0,-(index+1));
5468 daughter.SetKinkIndex(0,index+1);
5469 if (mother.GetNumberOfClusters()>50) {
5470 delete array->RemoveAt(i);
5471 array->AddAt(new AliTPCseed(mother),i);
5474 array->AddLast(new AliTPCseed(mother));
5476 array->AddLast(new AliTPCseed(daughter));
5477 for (Int_t icl=0;icl<row0;icl++) {
5478 if (mother.GetClusterPointer(icl)) mother.GetClusterPointer(icl)->Use(20);
5481 for (Int_t icl=row0;icl<158;icl++) {
5482 if (daughter.GetClusterPointer(icl)) daughter.GetClusterPointer(icl)->Use(20);
5491 delete [] daughters;
5513 AliInfo(Form("Ncandidates=\t%d\t%d\t%d\t%d\n",esd->GetNumberOfKinks(),ncandidates,ntracks,nall));
5518 Int_t AliTPCtrackerMI::RefitKink(AliTPCseed &mother, AliTPCseed &daughter, const AliESDkink &knk)
5521 // refit kink towards to the vertex
5524 AliKink &kink=(AliKink &)knk;
5526 Int_t row0 = GetRowNumber(kink.GetR());
5527 FollowProlongation(mother,0);
5528 mother.Reset(kFALSE);
5530 FollowProlongation(daughter,row0);
5531 daughter.Reset(kFALSE);
5532 FollowBackProlongation(daughter,158);
5533 daughter.Reset(kFALSE);
5534 Int_t first = TMath::Max(row0-20,30);
5535 Int_t last = TMath::Min(row0+20,140);
5537 const Int_t kNdiv =5;
5538 AliTPCseed param0[kNdiv]; // parameters along the track
5539 AliTPCseed param1[kNdiv]; // parameters along the track
5540 AliKink kinks[kNdiv]; // corresponding kink parameters
5543 for (Int_t irow=0; irow<kNdiv;irow++){
5544 rows[irow] = first +((last-first)*irow)/(kNdiv-1);
5546 // store parameters along the track
5548 for (Int_t irow=0;irow<kNdiv;irow++){
5549 FollowBackProlongation(mother, rows[irow]);
5550 FollowProlongation(daughter,rows[kNdiv-1-irow]);
5551 param0[irow] = mother;
5552 param1[kNdiv-1-irow] = daughter;
5556 for (Int_t irow=0; irow<kNdiv-1;irow++){
5557 if (param0[irow].GetNumberOfClusters()<kNdiv||param1[irow].GetNumberOfClusters()<kNdiv) continue;
5558 kinks[irow].SetMother(param0[irow]);
5559 kinks[irow].SetDaughter(param1[irow]);
5560 kinks[irow].Update();
5563 // choose kink with best "quality"
5565 Double_t mindist = 10000;
5566 for (Int_t irow=0;irow<kNdiv;irow++){
5567 if (param0[irow].GetNumberOfClusters()<20||param1[irow].GetNumberOfClusters()<20) continue;
5568 if (TMath::Abs(kinks[irow].GetR())>240.) continue;
5569 if (TMath::Abs(kinks[irow].GetR())<100.) continue;
5571 Float_t normdist = TMath::Abs(param0[irow].GetX()-kinks[irow].GetR())*(0.1+kink.GetDistance());
5572 normdist/= (param0[irow].GetNumberOfClusters()+param1[irow].GetNumberOfClusters()+40.);
5573 if (normdist < mindist){
5579 if (index==-1) return 0;
5582 param0[index].Reset(kTRUE);
5583 FollowProlongation(param0[index],0);
5585 mother = param0[index];
5586 daughter = param1[index]; // daughter in vertex
5588 kink.SetMother(mother);
5589 kink.SetDaughter(daughter);
5591 kink.SetTPCRow0(GetRowNumber(kink.GetR()));
5592 kink.SetTPCncls(param0[index].GetNumberOfClusters(),0);
5593 kink.SetTPCncls(param1[index].GetNumberOfClusters(),1);
5594 kink.SetLabel(CookLabel(&mother,0.4, 0,kink.GetTPCRow0()),0);
5595 kink.SetLabel(CookLabel(&daughter,0.4, kink.GetTPCRow0(),160),1);
5596 mother.SetLabel(kink.GetLabel(0));
5597 daughter.SetLabel(kink.GetLabel(1));
5603 void AliTPCtrackerMI::UpdateKinkQualityM(AliTPCseed * seed){
5605 // update Kink quality information for mother after back propagation
5607 if (seed->GetKinkIndex(0)>=0) return;
5608 for (Int_t ikink=0;ikink<3;ikink++){
5609 Int_t index = seed->GetKinkIndex(ikink);
5610 if (index>=0) break;
5611 index = TMath::Abs(index)-1;
5612 AliESDkink * kink = fEvent->GetKink(index);
5613 kink->SetTPCDensity(-1,0,0);
5614 kink->SetTPCDensity(1,0,1);
5616 Int_t row0 = kink->GetTPCRow0() - 2 - Int_t( 0.5/ (0.05+kink->GetAngle(2)));
5617 if (row0<15) row0=15;
5619 Int_t row1 = kink->GetTPCRow0() + 2 + Int_t( 0.5/ (0.05+kink->GetAngle(2)));
5620 if (row1>145) row1=145;
5622 Int_t found,foundable,shared;
5623 seed->GetClusterStatistic(0,row0, found, foundable,shared,kFALSE);
5624 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),0,0);
5625 seed->GetClusterStatistic(row1,155, found, foundable,shared,kFALSE);
5626 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),0,1);
5631 void AliTPCtrackerMI::UpdateKinkQualityD(AliTPCseed * seed){
5633 // update Kink quality information for daughter after refit
5635 if (seed->GetKinkIndex(0)<=0) return;
5636 for (Int_t ikink=0;ikink<3;ikink++){
5637 Int_t index = seed->GetKinkIndex(ikink);
5638 if (index<=0) break;
5639 index = TMath::Abs(index)-1;
5640 AliESDkink * kink = fEvent->GetKink(index);
5641 kink->SetTPCDensity(-1,1,0);
5642 kink->SetTPCDensity(-1,1,1);
5644 Int_t row0 = kink->GetTPCRow0() -2 - Int_t( 0.5/ (0.05+kink->GetAngle(2)));
5645 if (row0<15) row0=15;
5647 Int_t row1 = kink->GetTPCRow0() +2 + Int_t( 0.5/ (0.05+kink->GetAngle(2)));
5648 if (row1>145) row1=145;
5650 Int_t found,foundable,shared;
5651 seed->GetClusterStatistic(0,row0, found, foundable,shared,kFALSE);
5652 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),1,0);
5653 seed->GetClusterStatistic(row1,155, found, foundable,shared,kFALSE);
5654 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),1,1);
5660 Int_t AliTPCtrackerMI::CheckKinkPoint(AliTPCseed*seed,AliTPCseed &mother, AliTPCseed &daughter, const AliESDkink &knk)
5663 // check kink point for given track
5664 // if return value=0 kink point not found
5665 // otherwise seed0 correspond to mother particle
5666 // seed1 correspond to daughter particle
5667 // kink parameter of kink point
5668 AliKink &kink=(AliKink &)knk;
5670 Int_t middlerow = (seed->GetFirstPoint()+seed->GetLastPoint())/2;
5671 Int_t first = seed->GetFirstPoint();
5672 Int_t last = seed->GetLastPoint();
5673 if (last-first<20) return 0; // shortest length - 2*30 = 60 pad-rows
5676 AliTPCseed *seed1 = ReSeed(seed,middlerow+20, kTRUE); //middle of chamber
5677 if (!seed1) return 0;
5678 FollowProlongation(*seed1,seed->GetLastPoint()-20);
5679 seed1->Reset(kTRUE);
5680 FollowProlongation(*seed1,158);
5681 seed1->Reset(kTRUE);
5682 last = seed1->GetLastPoint();
5684 AliTPCseed *seed0 = new AliTPCseed(*seed);
5685 seed0->Reset(kFALSE);
5688 AliTPCseed param0[20]; // parameters along the track
5689 AliTPCseed param1[20]; // parameters along the track
5690 AliKink kinks[20]; // corresponding kink parameters
5692 for (Int_t irow=0; irow<20;irow++){
5693 rows[irow] = first +((last-first)*irow)/19;
5695 // store parameters along the track
5697 for (Int_t irow=0;irow<20;irow++){
5698 FollowBackProlongation(*seed0, rows[irow]);
5699 FollowProlongation(*seed1,rows[19-irow]);
5700 param0[irow] = *seed0;
5701 param1[19-irow] = *seed1;
5705 for (Int_t irow=0; irow<19;irow++){
5706 kinks[irow].SetMother(param0[irow]);
5707 kinks[irow].SetDaughter(param1[irow]);
5708 kinks[irow].Update();
5711 // choose kink with biggest change of angle
5713 Double_t maxchange= 0;
5714 for (Int_t irow=1;irow<19;irow++){
5715 if (TMath::Abs(kinks[irow].GetR())>240.) continue;
5716 if (TMath::Abs(kinks[irow].GetR())<110.) continue;
5717 Float_t quality = TMath::Abs(kinks[irow].GetAngle(2))/(3.+TMath::Abs(kinks[irow].GetR()-param0[irow].GetX()));
5718 if ( quality > maxchange){
5719 maxchange = quality;
5726 if (index<0) return 0;
5728 Int_t row0 = GetRowNumber(kinks[index].GetR()); //row 0 estimate
5729 seed0 = new AliTPCseed(param0[index]);
5730 seed1 = new AliTPCseed(param1[index]);
5731 seed0->Reset(kFALSE);
5732 seed1->Reset(kFALSE);
5733 seed0->ResetCovariance(10.);
5734 seed1->ResetCovariance(10.);
5735 FollowProlongation(*seed0,0);
5736 FollowBackProlongation(*seed1,158);
5737 mother = *seed0; // backup mother at position 0
5738 seed0->Reset(kFALSE);
5739 seed1->Reset(kFALSE);
5740 seed0->ResetCovariance(10.);
5741 seed1->ResetCovariance(10.);
5743 first = TMath::Max(row0-20,0);
5744 last = TMath::Min(row0+20,158);
5746 for (Int_t irow=0; irow<20;irow++){
5747 rows[irow] = first +((last-first)*irow)/19;
5749 // store parameters along the track
5751 for (Int_t irow=0;irow<20;irow++){
5752 FollowBackProlongation(*seed0, rows[irow]);
5753 FollowProlongation(*seed1,rows[19-irow]);
5754 param0[irow] = *seed0;
5755 param1[19-irow] = *seed1;
5759 for (Int_t irow=0; irow<19;irow++){
5760 kinks[irow].SetMother(param0[irow]);
5761 kinks[irow].SetDaughter(param1[irow]);
5762 // param0[irow].Dump();
5763 //param1[irow].Dump();
5764 kinks[irow].Update();
5767 // choose kink with biggest change of angle
5770 for (Int_t irow=0;irow<20;irow++){
5771 if (TMath::Abs(kinks[irow].GetR())>250.) continue;
5772 if (TMath::Abs(kinks[irow].GetR())<90.) continue;
5773 Float_t quality = TMath::Abs(kinks[irow].GetAngle(2))/(3.+TMath::Abs(kinks[irow].GetR()-param0[irow].GetX()));
5774 if ( quality > maxchange){
5775 maxchange = quality;
5782 if (index==-1 || param0[index].GetNumberOfClusters()+param1[index].GetNumberOfClusters()<100){
5788 // Float_t anglesigma = TMath::Sqrt(param0[index].fC22+param0[index].fC33+param1[index].fC22+param1[index].fC33);
5790 kink.SetMother(param0[index]);
5791 kink.SetDaughter(param1[index]);
5794 Double_t chi2P2 = param0[index].GetParameter()[2]-param1[index].GetParameter()[2];
5796 chi2P2/=param0[index].GetCovariance()[5]+param1[index].GetCovariance()[5];
5797 Double_t chi2P3 = param0[index].GetParameter()[3]-param1[index].GetParameter()[3];
5799 chi2P3/=param0[index].GetCovariance()[9]+param1[index].GetCovariance()[9];
5801 if (AliTPCReconstructor::StreamLevel()>1) {
5802 (*fDebugStreamer)<<"kinkHpt"<<
5805 "p0.="<<¶m0[index]<<
5806 "p1.="<<¶m1[index]<<
5810 if ( chi2P2+chi2P3<AliTPCReconstructor::GetRecoParam()->GetKinkAngleCutChi2(0)){
5817 row0 = GetRowNumber(kink.GetR());
5818 kink.SetTPCRow0(row0);
5819 kink.SetLabel(CookLabel(seed0,0.5,0,row0),0);
5820 kink.SetLabel(CookLabel(seed1,0.5,row0,158),1);
5821 kink.SetIndex(-10,0);
5822 kink.SetIndex(int(param0[index].GetNumberOfClusters()+param1[index].GetNumberOfClusters()),1);
5823 kink.SetTPCncls(param0[index].GetNumberOfClusters(),0);
5824 kink.SetTPCncls(param1[index].GetNumberOfClusters(),1);
5827 // new (&mother) AliTPCseed(param0[index]);
5828 daughter = param1[index];
5829 daughter.SetLabel(kink.GetLabel(1));
5830 param0[index].Reset(kTRUE);
5831 FollowProlongation(param0[index],0);
5832 mother = param0[index];
5833 mother.SetLabel(kink.GetLabel(0));
5834 if ( chi2P2+chi2P3<AliTPCReconstructor::GetRecoParam()->GetKinkAngleCutChi2(1)){
5846 AliTPCseed* AliTPCtrackerMI::ReSeed(AliTPCseed *t)
5849 // reseed - refit - track
5852 // Int_t last = fSectors->GetNRows()-1;
5854 if (fSectors == fOuterSec){
5855 first = TMath::Max(first, t->GetFirstPoint()-fInnerSec->GetNRows());
5859 first = t->GetFirstPoint();
5861 AliTPCseed * seed = MakeSeed(t,0.1,0.5,0.9);
5862 FollowBackProlongation(*t,fSectors->GetNRows()-1);
5864 FollowProlongation(*t,first);
5874 //_____________________________________________________________________________
5875 Int_t AliTPCtrackerMI::ReadSeeds(const TFile *inp) {
5876 //-----------------------------------------------------------------
5877 // This function reades track seeds.
5878 //-----------------------------------------------------------------
5879 TDirectory *savedir=gDirectory;
5881 TFile *in=(TFile*)inp;
5882 if (!in->IsOpen()) {
5883 cerr<<"AliTPCtrackerMI::ReadSeeds(): input file is not open !\n";
5888 TTree *seedTree=(TTree*)in->Get("Seeds");
5890 cerr<<"AliTPCtrackerMI::ReadSeeds(): ";
5891 cerr<<"can't get a tree with track seeds !\n";
5894 AliTPCtrack *seed=new AliTPCtrack;
5895 seedTree->SetBranchAddress("tracks",&seed);
5897 if (fSeeds==0) fSeeds=new TObjArray(15000);
5899 Int_t n=(Int_t)seedTree->GetEntries();
5900 for (Int_t i=0; i<n; i++) {
5901 seedTree->GetEvent(i);
5902 fSeeds->AddLast(new AliTPCseed(*seed/*,seed->GetAlpha()*/));
5911 Int_t AliTPCtrackerMI::Clusters2Tracks (AliESDEvent *const esd)
5914 // clusters to tracks
5916 if (fSeeds) DeleteSeeds();
5919 AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ;
5920 transform->SetCurrentTimeStamp( esd->GetTimeStamp());
5921 transform->SetCurrentRun(esd->GetRunNumber());
5924 if (!fSeeds) return 1;
5926 if (AliTPCReconstructor::StreamLevel()>3) DumpClusters(0,fSeeds);
5932 //_____________________________________________________________________________
5933 Int_t AliTPCtrackerMI::Clusters2Tracks() {
5934 //-----------------------------------------------------------------
5935 // This is a track finder.
5936 //-----------------------------------------------------------------
5937 TDirectory *savedir=gDirectory;
5941 fSeeds = Tracking();
5944 Info("Clusters2Tracks","Time for tracking: \t");timer.Print();timer.Start();
5946 //activate again some tracks
5947 for (Int_t i=0; i<fSeeds->GetEntriesFast(); i++) {
5948 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
5950 Int_t nc=t.GetNumberOfClusters();
5952 delete fSeeds->RemoveAt(i);
5956 if (pt->GetRemoval()==10) {
5957 if (pt->GetDensityFirst(20)>0.8 || pt->GetDensityFirst(30)>0.8 || pt->GetDensityFirst(40)>0.7)
5958 pt->Desactivate(10); // make track again active // MvL: should be 0 ?
5960 pt->Desactivate(20);
5961 delete fSeeds->RemoveAt(i);
5966 RemoveUsed2(fSeeds,0.85,0.85,0);
5967 if (AliTPCReconstructor::GetRecoParam()->GetDoKinks()) FindKinks(fSeeds,fEvent);
5968 //FindCurling(fSeeds, fEvent,0);
5969 if (AliTPCReconstructor::StreamLevel()>5) FindMultiMC(fSeeds, fEvent,-1); // find multi found tracks
5970 RemoveUsed2(fSeeds,0.5,0.4,20);
5971 FindSplitted(fSeeds, fEvent,0); // find multi found tracks
5972 if (AliTPCReconstructor::StreamLevel()>5) FindMultiMC(fSeeds, fEvent,0); // find multi found tracks
5975 // // refit short tracks
5977 Int_t nseed=fSeeds->GetEntriesFast();
5980 for (Int_t i=0; i<nseed; i++) {
5981 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
5983 Int_t nc=t.GetNumberOfClusters();
5985 delete fSeeds->RemoveAt(i);
5988 CookLabel(pt,0.1); //For comparison only
5989 //if ((pt->IsActive() || (pt->fRemoval==10) )&& nc>50 &&pt->GetNumberOfClusters()>0.4*pt->fNFoundable){
5990 if ((pt->IsActive() || (pt->GetRemoval()==10) )){
5992 if (fDebug>0) cerr<<found<<'\r';
5996 delete fSeeds->RemoveAt(i);
6000 //RemoveOverlap(fSeeds,0.99,7,kTRUE);
6002 //RemoveUsed(fSeeds,0.9,0.9,6);
6004 nseed=fSeeds->GetEntriesFast();
6006 for (Int_t i=0; i<nseed; i++) {
6007 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
6009 Int_t nc=t.GetNumberOfClusters();
6011 delete fSeeds->RemoveAt(i);
6015 t.CookdEdx(0.02,0.6);
6016 // CheckKinkPoint(&t,0.05);
6017 //if ((pt->IsActive() || (pt->fRemoval==10) )&& nc>50 &&pt->GetNumberOfClusters()>0.4*pt->fNFoundable){
6018 if ((pt->IsActive() || (pt->GetRemoval()==10) )){
6026 delete fSeeds->RemoveAt(i);
6027 //AliTPCseed * seed1 = ReSeed(pt,0.05,0.5,1);
6029 // FollowProlongation(*seed1,0);
6030 // Int_t n = seed1->GetNumberOfClusters();
6031 // printf("fP4\t%f\t%f\n",seed1->GetC(),pt->GetC());
6032 // printf("fN\t%d\t%d\n", seed1->GetNumberOfClusters(),pt->GetNumberOfClusters());
6035 //AliTPCseed * seed2 = ReSeed(pt,0.95,0.5,0.05);
6039 SortTracks(fSeeds, 1);
6043 PrepareForBackProlongation(fSeeds,5.);
6044 PropagateBack(fSeeds);
6045 printf("Time for back propagation: \t");timer.Print();timer.Start();
6049 PrepareForProlongation(fSeeds,5.);
6050 PropagateForard2(fSeeds);
6052 printf("Time for FORWARD propagation: \t");timer.Print();timer.Start();
6053 // RemoveUsed(fSeeds,0.7,0.7,6);
6054 //RemoveOverlap(fSeeds,0.9,7,kTRUE);
6056 nseed=fSeeds->GetEntriesFast();
6058 for (Int_t i=0; i<nseed; i++) {
6059 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
6061 Int_t nc=t.GetNumberOfClusters();
6063 delete fSeeds->RemoveAt(i);
6066 t.CookdEdx(0.02,0.6);
6067 // CookLabel(pt,0.1); //For comparison only
6068 //if ((pt->IsActive() || (pt->fRemoval==10) )&& nc>50 &&pt->GetNumberOfClusters()>0.4*pt->fNFoundable){
6069 if ((pt->IsActive() || (pt->fRemoval==10) )){
6070 cerr<<found++<<'\r';
6073 delete fSeeds->RemoveAt(i);
6078 // fNTracks = found;
6080 Info("Clusters2Tracks","Time for overlap removal, track writing and dedx cooking: \t"); timer.Print();timer.Start();
6083 // cerr<<"Number of found tracks : "<<"\t"<<found<<endl;
6084 Info("Clusters2Tracks","Number of found tracks %d",found);
6086 // UnloadClusters();
6091 void AliTPCtrackerMI::Tracking(TObjArray * arr)
6094 // tracking of the seeds
6097 fSectors = fOuterSec;
6098 ParallelTracking(arr,150,63);
6099 fSectors = fOuterSec;
6100 ParallelTracking(arr,63,0);
6103 TObjArray * AliTPCtrackerMI::Tracking(Int_t seedtype, Int_t i1, Int_t i2, Float_t cuts[4], Float_t dy, Int_t dsec)
6108 TObjArray * arr = new TObjArray;
6110 fSectors = fOuterSec;
6113 for (Int_t sec=0;sec<fkNOS;sec++){
6114 if (seedtype==3) MakeSeeds3(arr,sec,i1,i2,cuts,dy, dsec);
6115 if (seedtype==4) MakeSeeds5(arr,sec,i1,i2,cuts,dy);
6116 if (seedtype==2) MakeSeeds2(arr,sec,i1,i2,cuts,dy);
6119 Info("Tracking","\nSeeding - %d\t%d\t%d\t%d\n",seedtype,i1,i2,arr->GetEntriesFast());
6131 TObjArray * AliTPCtrackerMI::Tracking()
6135 if (AliTPCReconstructor::GetRecoParam()->GetSpecialSeeding()) return TrackingSpecial();
6138 Int_t nup=fOuterSec->GetNRows()+fInnerSec->GetNRows();
6140 TObjArray * seeds = new TObjArray;
6142 Int_t fLastSeedRowSec=AliTPCReconstructor::GetRecoParam()->GetLastSeedRowSec();
6143 Int_t gapPrim = AliTPCReconstructor::GetRecoParam()->GetSeedGapPrim();
6144 Int_t gapSec = AliTPCReconstructor::GetRecoParam()->GetSeedGapSec();
6152 Float_t fnumber = 3.0;
6153 Float_t fdensity = 3.0;
6158 for (Int_t delta = 0; delta<18; delta+=gapPrim){
6162 arr = Tracking(3,nup-1-delta,nup-1-delta-gap,cuts,-1,1);
6163 SumTracks(seeds,arr);
6164 SignClusters(seeds,fnumber,fdensity);
6166 for (Int_t i=2;i<6;i+=2){
6167 // seed high pt tracks
6170 arr = Tracking(3,nup-i-delta,nup-i-delta-gap,cuts,-1,0);
6171 SumTracks(seeds,arr);
6172 SignClusters(seeds,fnumber,fdensity);
6177 // RemoveUsed(seeds,0.9,0.9,1);
6178 // UnsignClusters();
6179 // SignClusters(seeds,fnumber,fdensity);
6183 for (Int_t delta = 20; delta<120; delta+=gapPrim){
6185 // seed high pt tracks
6189 arr = Tracking(3,nup-delta,nup-delta-gap,cuts,-1);
6190 SumTracks(seeds,arr);
6191 SignClusters(seeds,fnumber,fdensity);
6196 arr = Tracking(3,nup-delta-5,nup-delta-5-gap,cuts,-1);
6197 SumTracks(seeds,arr);
6198 SignClusters(seeds,fnumber,fdensity);
6209 Info("Tracking()","\n\nPrimary seeding\t%d\n\n",seeds->GetEntriesFast());
6213 // RemoveUsed(seeds,0.75,0.75,1);
6215 //SignClusters(seeds,fnumber,fdensity);
6224 arr = Tracking(4,nup-1,nup-1-gap,cuts,-1);
6225 SumTracks(seeds,arr);
6226 SignClusters(seeds,fnumber,fdensity);
6228 arr = Tracking(4,nup-2,nup-2-gap,cuts,-1);
6229 SumTracks(seeds,arr);
6230 SignClusters(seeds,fnumber,fdensity);
6232 arr = Tracking(4,nup-3,nup-3-gap,cuts,-1);
6233 SumTracks(seeds,arr);
6234 SignClusters(seeds,fnumber,fdensity);
6236 arr = Tracking(4,nup-5,nup-5-gap,cuts,-1);
6237 SumTracks(seeds,arr);
6238 SignClusters(seeds,fnumber,fdensity);
6240 arr = Tracking(4,nup-7,nup-7-gap,cuts,-1);
6241 SumTracks(seeds,arr);
6242 SignClusters(seeds,fnumber,fdensity);
6245 arr = Tracking(4,nup-9,nup-9-gap,cuts,-1);
6246 SumTracks(seeds,arr);
6247 SignClusters(seeds,fnumber,fdensity);
6251 for (Int_t delta = 9; delta<30; delta+=gapSec){
6257 arr = Tracking(4,nup-1-delta,nup-1-delta-gap,cuts,-1);
6258 SumTracks(seeds,arr);
6259 SignClusters(seeds,fnumber,fdensity);
6261 arr = Tracking(4,nup-3-delta,nup-5-delta-gap,cuts,4);
6262 SumTracks(seeds,arr);
6263 SignClusters(seeds,fnumber,fdensity);
6276 for (Int_t delta = 30; delta<fLastSeedRowSec; delta+=gapSec){
6282 arr = Tracking(4,nup-1-delta,nup-1-delta-gap,cuts,-1);
6283 SumTracks(seeds,arr);
6284 SignClusters(seeds,fnumber,fdensity);
6286 arr = Tracking(4,nup-5-delta,nup-5-delta-gap,cuts,5 );
6287 SumTracks(seeds,arr);
6288 SignClusters(seeds,fnumber,fdensity);
6292 Info("Tracking()","\n\nSecondary seeding\t%d\n\n",seeds->GetEntriesFast());
6303 TObjArray * AliTPCtrackerMI::TrackingSpecial()
6306 // seeding adjusted for laser and cosmic tests - short tracks with big inclination angle
6307 // no primary vertex seeding tried
6311 Int_t nup=fOuterSec->GetNRows()+fInnerSec->GetNRows();
6313 TObjArray * seeds = new TObjArray;
6318 Float_t fnumber = 3.0;
6319 Float_t fdensity = 3.0;
6322 cuts[0] = AliTPCReconstructor::GetRecoParam()->GetMaxC(); // max curvature
6323 cuts[1] = 3.5; // max tan(phi) angle for seeding
6324 cuts[2] = 3.; // not used (cut on z primary vertex)
6325 cuts[3] = 3.5; // max tan(theta) angle for seeding
6327 for (Int_t delta = 0; nup-delta-gap-1>0; delta+=3){
6329 arr = Tracking(4,nup-1-delta,nup-1-delta-gap,cuts,-1);
6330 SumTracks(seeds,arr);
6331 SignClusters(seeds,fnumber,fdensity);
6335 Info("Tracking()","\n\nSecondary seeding\t%d\n\n",seeds->GetEntriesFast());
6346 void AliTPCtrackerMI::SumTracks(TObjArray *arr1,TObjArray *&arr2) const
6349 //sum tracks to common container
6350 //remove suspicious tracks
6351 Int_t nseed = arr2->GetEntriesFast();
6352 for (Int_t i=0;i<nseed;i++){
6353 AliTPCseed *pt=(AliTPCseed*)arr2->UncheckedAt(i);
6356 // remove tracks with too big curvature
6358 if (TMath::Abs(pt->GetC())>AliTPCReconstructor::GetRecoParam()->GetMaxC()){
6359 delete arr2->RemoveAt(i);
6362 // REMOVE VERY SHORT TRACKS
6363 if (pt->GetNumberOfClusters()<20){
6364 delete arr2->RemoveAt(i);
6367 // NORMAL ACTIVE TRACK
6368 if (pt->IsActive()){
6369 arr1->AddLast(arr2->RemoveAt(i));
6372 //remove not usable tracks
6373 if (pt->GetRemoval()!=10){
6374 delete arr2->RemoveAt(i);
6378 // ENABLE ONLY ENOUGH GOOD STOPPED TRACKS
6379 if (pt->GetDensityFirst(20)>0.8 || pt->GetDensityFirst(30)>0.8 || pt->GetDensityFirst(40)>0.7)
6380 arr1->AddLast(arr2->RemoveAt(i));
6382 delete arr2->RemoveAt(i);
6386 delete arr2; arr2 = 0;
6391 void AliTPCtrackerMI::ParallelTracking(TObjArray *const arr, Int_t rfirst, Int_t rlast)
6394 // try to track in parralel
6396 Int_t nseed=arr->GetEntriesFast();
6397 //prepare seeds for tracking
6398 for (Int_t i=0; i<nseed; i++) {
6399 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i), &t=*pt;
6401 if (!t.IsActive()) continue;
6402 // follow prolongation to the first layer
6403 if ( (fSectors ==fInnerSec) || (t.GetFirstPoint()-fkParam->GetNRowLow()>rfirst+1) )
6404 FollowProlongation(t, rfirst+1);
6409 for (Int_t nr=rfirst; nr>=rlast; nr--){
6410 if (nr<fInnerSec->GetNRows())
6411 fSectors = fInnerSec;
6413 fSectors = fOuterSec;
6414 // make indexes with the cluster tracks for given
6416 // find nearest cluster
6417 for (Int_t i=0; i<nseed; i++) {
6418 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i), &t=*pt;
6420 if (nr==80) pt->UpdateReference();
6421 if (!pt->IsActive()) continue;
6422 // if ( (fSectors ==fOuterSec) && (pt->fFirstPoint-fkParam->GetNRowLow())<nr) continue;
6423 if (pt->GetRelativeSector()>17) {
6426 UpdateClusters(t,nr);
6428 // prolonagate to the nearest cluster - if founded
6429 for (Int_t i=0; i<nseed; i++) {
6430 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
6432 if (!pt->IsActive()) continue;
6433 // if ((fSectors ==fOuterSec) && (pt->fFirstPoint-fkParam->GetNRowLow())<nr) continue;
6434 if (pt->GetRelativeSector()>17) {
6437 FollowToNextCluster(*pt,nr);
6442 void AliTPCtrackerMI::PrepareForBackProlongation(const TObjArray *const arr,Float_t fac) const
6446 // if we use TPC track itself we have to "update" covariance
6448 Int_t nseed= arr->GetEntriesFast();
6449 for (Int_t i=0;i<nseed;i++){
6450 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
6454 //rotate to current local system at first accepted point
6455 Int_t index = pt->GetClusterIndex2(pt->GetFirstPoint());
6456 Int_t sec = (index&0xff000000)>>24;
6458 Float_t angle1 = fInnerSec->GetAlpha()*sec+fInnerSec->GetAlphaShift();
6459 if (angle1>TMath::Pi())
6460 angle1-=2.*TMath::Pi();
6461 Float_t angle2 = pt->GetAlpha();
6463 if (TMath::Abs(angle1-angle2)>0.001){
6464 if (!pt->Rotate(angle1-angle2)) return;
6465 //angle2 = pt->GetAlpha();
6466 //pt->fRelativeSector = pt->GetAlpha()/fInnerSec->GetAlpha();
6467 //if (pt->GetAlpha()<0)
6468 // pt->fRelativeSector+=18;
6469 //sec = pt->fRelativeSector;
6478 void AliTPCtrackerMI::PrepareForProlongation(TObjArray *const arr, Float_t fac) const
6482 // if we use TPC track itself we have to "update" covariance
6484 Int_t nseed= arr->GetEntriesFast();
6485 for (Int_t i=0;i<nseed;i++){
6486 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
6489 pt->SetFirstPoint(pt->GetLastPoint());
6497 Int_t AliTPCtrackerMI::PropagateBack(const TObjArray *const arr)
6500 // make back propagation
6502 Int_t nseed= arr->GetEntriesFast();
6503 for (Int_t i=0;i<nseed;i++){
6504 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
6505 if (pt&& pt->GetKinkIndex(0)<=0) {
6506 //AliTPCseed *pt2 = new AliTPCseed(*pt);
6507 fSectors = fInnerSec;
6508 //FollowBackProlongation(*pt,fInnerSec->GetNRows()-1);
6509 //fSectors = fOuterSec;
6510 FollowBackProlongation(*pt,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1,1);
6511 //if (pt->GetNumberOfClusters()<(pt->fEsd->GetTPCclusters(0)) ){
6512 // Error("PropagateBack","Not prolonged track %d",pt->GetLabel());
6513 // FollowBackProlongation(*pt2,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1);
6516 if (pt&& pt->GetKinkIndex(0)>0) {
6517 AliESDkink * kink = fEvent->GetKink(pt->GetKinkIndex(0)-1);
6518 pt->SetFirstPoint(kink->GetTPCRow0());
6519 fSectors = fInnerSec;
6520 FollowBackProlongation(*pt,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1,1);
6528 Int_t AliTPCtrackerMI::PropagateForward2(const TObjArray *const arr)
6531 // make forward propagation
6533 Int_t nseed= arr->GetEntriesFast();
6535 for (Int_t i=0;i<nseed;i++){
6536 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
6538 FollowProlongation(*pt,0,1,1);
6547 Int_t AliTPCtrackerMI::PropagateForward()
6550 // propagate track forward
6552 Int_t nseed = fSeeds->GetEntriesFast();
6553 for (Int_t i=0;i<nseed;i++){
6554 AliTPCseed *pt = (AliTPCseed*)fSeeds->UncheckedAt(i);
6556 AliTPCseed &t = *pt;
6557 Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
6558 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
6559 if (alpha < 0. ) alpha += 2.*TMath::Pi();
6560 t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN);
6564 fSectors = fOuterSec;
6565 ParallelTracking(fSeeds,fOuterSec->GetNRows()+fInnerSec->GetNRows()-1,fInnerSec->GetNRows());
6566 fSectors = fInnerSec;
6567 ParallelTracking(fSeeds,fInnerSec->GetNRows()-1,0);
6576 Int_t AliTPCtrackerMI::PropagateBack(AliTPCseed *const pt, Int_t row0, Int_t row1)
6579 // make back propagation, in between row0 and row1
6583 fSectors = fInnerSec;
6586 if (row1<fSectors->GetNRows())
6589 r1 = fSectors->GetNRows()-1;
6591 if (row0<fSectors->GetNRows()&& r1>0 )
6592 FollowBackProlongation(*pt,r1);
6593 if (row1<=fSectors->GetNRows())
6596 r1 = row1 - fSectors->GetNRows();
6597 if (r1<=0) return 0;
6598 if (r1>=fOuterSec->GetNRows()) return 0;
6599 fSectors = fOuterSec;
6600 return FollowBackProlongation(*pt,r1);
6608 void AliTPCtrackerMI::GetShape(AliTPCseed * seed, Int_t row)
6610 // gets cluster shape
6612 AliTPCClusterParam * clparam = AliTPCcalibDB::Instance()->GetClusterParam();
6613 Float_t zdrift = TMath::Abs((fkParam->GetZLength(0)-TMath::Abs(seed->GetZ())));
6614 Int_t type = (seed->GetSector() < fkParam->GetNSector()/2) ? 0: (row>126) ? 1:2;
6615 Double_t angulary = seed->GetSnp();
6617 if (TMath::Abs(angulary)>AliTPCReconstructor::GetMaxSnpTracker()) {
6618 angulary = TMath::Sign(AliTPCReconstructor::GetMaxSnpTracker(),angulary);
6621 angulary = angulary*angulary/((1.-angulary)*(1.+angulary));
6622 Double_t angularz = seed->GetTgl()*seed->GetTgl()*(1.+angulary);
6624 Double_t sigmay = clparam->GetRMS0(0,type,zdrift,TMath::Sqrt(TMath::Abs(angulary)));
6625 Double_t sigmaz = clparam->GetRMS0(1,type,zdrift,TMath::Sqrt(TMath::Abs(angularz)));
6626 seed->SetCurrentSigmaY2(sigmay*sigmay);
6627 seed->SetCurrentSigmaZ2(sigmaz*sigmaz);
6628 // Float_t sd2 = TMath::Abs((fkParam->GetZLength(0)-TMath::Abs(seed->GetZ())))*fkParam->GetDiffL()*fkParam->GetDiffL();
6629 // // Float_t padlength = fkParam->GetPadPitchLength(seed->fSector);
6630 // Float_t padlength = GetPadPitchLength(row);
6632 // Float_t sresy = (seed->GetSector() < fkParam->GetNSector()/2) ? 0.2 :0.3;
6633 // seed->SetCurrentSigmaY2(sd2+padlength*padlength*angulary/12.+sresy*sresy);
6635 // Float_t sresz = fkParam->GetZSigma();
6636 // seed->SetCurrentSigmaZ2(sd2+padlength*padlength*angularz*angularz*(1+angulary)/12.+sresz*sresz);
6638 Float_t wy = GetSigmaY(seed);
6639 Float_t wz = GetSigmaZ(seed);
6642 if (TMath::Abs(wy/seed->fCurrentSigmaY2-1)>0.0001 || TMath::Abs(wz/seed->fCurrentSigmaZ2-1)>0.0001 ){
6643 printf("problem\n");
6650 //__________________________________________________________________________
6651 void AliTPCtrackerMI::CookLabel(AliKalmanTrack *tk, Float_t wrong) const {
6652 //--------------------------------------------------------------------
6653 //This function "cooks" a track label. If label<0, this track is fake.
6654 //--------------------------------------------------------------------
6655 AliTPCseed * t = dynamic_cast<AliTPCseed*>(tk);
6657 printf("%s:%d wrong type \n",(char*)__FILE__,__LINE__);
6661 Int_t noc=t->GetNumberOfClusters();
6663 //printf("\nnot founded prolongation\n\n\n");
6669 AliTPCclusterMI *clusters[160];
6671 for (Int_t i=0;i<160;i++) {
6678 for (i=0; i<160 && current<noc; i++) {
6680 Int_t index=t->GetClusterIndex2(i);
6681 if (index<=0) continue;
6682 if (index&0x8000) continue;
6684 //clusters[current]=GetClusterMI(index);
6685 if (t->GetClusterPointer(i)){
6686 clusters[current]=t->GetClusterPointer(i);
6692 Int_t lab=123456789;
6693 for (i=0; i<noc; i++) {
6694 AliTPCclusterMI *c=clusters[i];
6696 lab=TMath::Abs(c->GetLabel(0));
6698 for (j=0; j<noc; j++) if (lb[j]==lab || mx[j]==0) break;
6704 for (i=0; i<noc; i++) if (mx[i]>max) {max=mx[i]; lab=lb[i];}
6706 for (i=0; i<noc; i++) {
6707 AliTPCclusterMI *c=clusters[i];
6709 if (TMath::Abs(c->GetLabel(1)) == lab ||
6710 TMath::Abs(c->GetLabel(2)) == lab ) max++;
6712 if (noc<=0) { lab=-1; return;}
6713 if ((1.- Float_t(max)/(noc)) > wrong) lab=-lab;
6716 Int_t tail=Int_t(0.10*noc);
6719 for (i=1; i<160&&ind<tail; i++) {
6720 // AliTPCclusterMI *c=clusters[noc-i];
6721 AliTPCclusterMI *c=clusters[i];
6723 if (lab == TMath::Abs(c->GetLabel(0)) ||
6724 lab == TMath::Abs(c->GetLabel(1)) ||
6725 lab == TMath::Abs(c->GetLabel(2))) max++;
6728 if (max < Int_t(0.5*tail)) lab=-lab;
6735 //delete[] clusters;
6739 //__________________________________________________________________________
6740 Int_t AliTPCtrackerMI::CookLabel(AliTPCseed *const t, Float_t wrong,Int_t first, Int_t last) const {
6741 //--------------------------------------------------------------------
6742 //This function "cooks" a track label. If label<0, this track is fake.
6743 //--------------------------------------------------------------------
6744 Int_t noc=t->GetNumberOfClusters();
6746 //printf("\nnot founded prolongation\n\n\n");
6752 AliTPCclusterMI *clusters[160];
6754 for (Int_t i=0;i<160;i++) {
6761 for (i=0; i<160 && current<noc; i++) {
6762 if (i<first) continue;
6763 if (i>last) continue;
6764 Int_t index=t->GetClusterIndex2(i);
6765 if (index<=0) continue;
6766 if (index&0x8000) continue;
6768 //clusters[current]=GetClusterMI(index);
6769 if (t->GetClusterPointer(i)){
6770 clusters[current]=t->GetClusterPointer(i);
6775 //if (noc<5) return -1;
6776 Int_t lab=123456789;
6777 for (i=0; i<noc; i++) {
6778 AliTPCclusterMI *c=clusters[i];
6780 lab=TMath::Abs(c->GetLabel(0));
6782 for (j=0; j<noc; j++) if (lb[j]==lab || mx[j]==0) break;
6788 for (i=0; i<noc; i++) if (mx[i]>max) {max=mx[i]; lab=lb[i];}
6790 for (i=0; i<noc; i++) {
6791 AliTPCclusterMI *c=clusters[i];
6793 if (TMath::Abs(c->GetLabel(1)) == lab ||
6794 TMath::Abs(c->GetLabel(2)) == lab ) max++;
6796 if (noc<=0) { lab=-1; return -1;}
6797 if ((1.- Float_t(max)/(noc)) > wrong) lab=-lab;
6800 Int_t tail=Int_t(0.10*noc);
6803 for (i=1; i<160&&ind<tail; i++) {
6804 // AliTPCclusterMI *c=clusters[noc-i];
6805 AliTPCclusterMI *c=clusters[i];
6807 if (lab == TMath::Abs(c->GetLabel(0)) ||
6808 lab == TMath::Abs(c->GetLabel(1)) ||
6809 lab == TMath::Abs(c->GetLabel(2))) max++;
6812 if (max < Int_t(0.5*tail)) lab=-lab;
6815 // t->SetLabel(lab);
6819 //delete[] clusters;
6823 Int_t AliTPCtrackerMI::GetRowNumber(Double_t x[3]) const
6825 //return pad row number for given x vector
6826 Float_t phi = TMath::ATan2(x[1],x[0]);
6827 if(phi<0) phi=2.*TMath::Pi()+phi;
6828 // Get the local angle in the sector philoc
6829 const Float_t kRaddeg = 180/3.14159265358979312;
6830 Float_t phiangle = (Int_t (phi*kRaddeg/20.) + 0.5)*20./kRaddeg;
6831 Double_t localx = x[0]*TMath::Cos(phiangle)-x[1]*TMath::Sin(phiangle);
6832 return GetRowNumber(localx);
6837 void AliTPCtrackerMI::MakeESDBitmaps(AliTPCseed *t, AliESDtrack *esd)
6839 //-----------------------------------------------------------------------
6840 // Fill the cluster and sharing bitmaps of the track
6841 //-----------------------------------------------------------------------
6843 Int_t firstpoint = 0;
6844 Int_t lastpoint = 159;
6845 AliTPCTrackerPoint *point;
6846 AliTPCclusterMI *cluster;
6849 TBits clusterMap(159);
6850 TBits sharedMap(159);
6852 for (int iter=firstpoint; iter<lastpoint; iter++) {
6853 // Change to cluster pointers to see if we have a cluster at given padrow
6855 cluster = t->GetClusterPointer(iter);
6857 clusterMap.SetBitNumber(iter, kTRUE);
6858 point = t->GetTrackPoint(iter);
6859 if (point->IsShared())
6860 sharedMap.SetBitNumber(iter,kTRUE);
6862 if (t->GetClusterIndex(iter) >= 0 && (t->GetClusterIndex(iter) & 0x8000) == 0) {
6863 fitMap.SetBitNumber(iter, kTRUE);
6867 esd->SetTPCClusterMap(clusterMap);
6868 esd->SetTPCSharedMap(sharedMap);
6869 esd->SetTPCFitMap(fitMap);
6870 if (nclsf != t->GetNumberOfClusters())
6871 AliWarning(Form("Inconsistency between ncls %d and indices %d (found %d)",t->GetNumberOfClusters(),nclsf,esd->GetTPCClusterMap().CountBits()));
6874 Bool_t AliTPCtrackerMI::IsFindable(AliTPCseed & track){
6876 // return flag if there is findable cluster at given position
6879 Float_t z = track.GetZ();
6881 if (TMath::Abs(z)<(AliTPCReconstructor::GetCtgRange()*track.GetX()+kDeltaZ) &&
6882 TMath::Abs(z)<fkParam->GetZLength(0) &&
6883 (TMath::Abs(track.GetSnp())<AliTPCReconstructor::GetMaxSnpTracker()))
6889 void AliTPCtrackerMI::AddCovariance(AliTPCseed * seed){
6891 // Adding systematic error estimate to the covariance matrix
6892 // !!!! the systematic error for element 4 is in 1/cm not in pt
6894 const Double_t *param = AliTPCReconstructor::GetRecoParam()->GetSystematicError();
6896 // use only the diagonal part if not specified otherwise
6897 if (!AliTPCReconstructor::GetRecoParam()->GetUseSystematicCorrelation()) return AddCovarianceAdd(seed);
6899 Double_t *covarS= (Double_t*)seed->GetCovariance();
6900 Double_t factor[5]={1,1,1,1,1};
6901 Double_t facC = AliTracker::GetBz()*kB2C;
6902 factor[0]= TMath::Sqrt(TMath::Abs((covarS[0] + param[0]*param[0])/covarS[0]));
6903 factor[1]= TMath::Sqrt(TMath::Abs((covarS[2] + param[1]*param[1])/covarS[2]));
6904 factor[2]= TMath::Sqrt(TMath::Abs((covarS[5] + param[2]*param[2])/covarS[5]));
6905 factor[3]= TMath::Sqrt(TMath::Abs((covarS[9] + param[3]*param[3])/covarS[9]));
6906 factor[4]= TMath::Sqrt(TMath::Abs((covarS[14] + facC*facC*param[4]*param[4])/covarS[14]));
6912 for (Int_t i=0; i<5; i++){
6913 for (Int_t j=i; j<5; j++){
6914 Int_t index=seed->GetIndex(i,j);
6915 covarS[index]*=factor[i]*factor[j];
6921 void AliTPCtrackerMI::AddCovarianceAdd(AliTPCseed * seed){
6923 // Adding systematic error - as additive factor without correlation
6925 const Double_t *param = AliTPCReconstructor::GetRecoParam()->GetSystematicError();
6926 Double_t *covarIn= (Double_t*)seed->GetCovariance();
6928 for (Int_t i=0;i<15;i++) covar[i]=0;
6934 covar[0] = param[0]*param[0];
6935 covar[2] = param[1]*param[1];
6936 covar[5] = param[2]*param[2];
6937 covar[9] = param[3]*param[3];
6938 Double_t facC = AliTracker::GetBz()*kB2C;
6939 covar[14]= param[4]*param[4]*facC*facC;
6941 covar[1]=TMath::Sqrt((covar[0]*covar[2]))*covarIn[1]/TMath::Sqrt((covarIn[0]*covarIn[2]));
6943 covar[3]=TMath::Sqrt((covar[0]*covar[5]))*covarIn[3]/TMath::Sqrt((covarIn[0]*covarIn[5]));
6944 covar[4]=TMath::Sqrt((covar[2]*covar[5]))*covarIn[4]/TMath::Sqrt((covarIn[2]*covarIn[5]));
6946 covar[6]=TMath::Sqrt((covar[0]*covar[9]))*covarIn[6]/TMath::Sqrt((covarIn[0]*covarIn[9]));
6947 covar[7]=TMath::Sqrt((covar[2]*covar[9]))*covarIn[7]/TMath::Sqrt((covarIn[2]*covarIn[9]));
6948 covar[8]=TMath::Sqrt((covar[5]*covar[9]))*covarIn[8]/TMath::Sqrt((covarIn[5]*covarIn[9]));
6950 covar[10]=TMath::Sqrt((covar[0]*covar[14]))*covarIn[10]/TMath::Sqrt((covarIn[0]*covarIn[14]));
6951 covar[11]=TMath::Sqrt((covar[2]*covar[14]))*covarIn[11]/TMath::Sqrt((covarIn[2]*covarIn[14]));
6952 covar[12]=TMath::Sqrt((covar[5]*covar[14]))*covarIn[12]/TMath::Sqrt((covarIn[5]*covarIn[14]));
6953 covar[13]=TMath::Sqrt((covar[9]*covar[14]))*covarIn[13]/TMath::Sqrt((covarIn[9]*covarIn[14]));
6955 seed->AddCovariance(covar);