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 t.SetCurrentClusterIndex1(tpcindex);
1618 Int_t relativesector = ((tpcindex&0xff000000)>>24)%18; // if previously accepted cluster in different sector
1619 Float_t angle = relativesector*fSectors->GetAlpha()+fSectors->GetAlphaShift();
1621 if (angle<-TMath::Pi()) angle += 2*TMath::Pi();
1622 if (angle>=TMath::Pi()) angle -= 2*TMath::Pi();
1624 if (TMath::Abs(angle-t.GetAlpha())>0.001){
1625 Double_t rotation = angle-t.GetAlpha();
1626 t.SetRelativeSector(relativesector);
1627 if (!t.Rotate(rotation)) return 0;
1629 if (!t.PropagateTo(x)) return 0;
1631 t.SetCurrentCluster(cl);
1633 Int_t accept = AcceptCluster(&t,t.GetCurrentCluster());
1634 if ((tpcindex&0x8000)==0) accept =0;
1636 //if founded cluster is acceptible
1637 if (cl->IsUsed(11)) { // id cluster is shared inrease uncertainty
1638 t.SetErrorY2(t.GetErrorY2()+0.03);
1639 t.SetErrorZ2(t.GetErrorZ2()+0.03);
1640 t.SetErrorY2(t.GetErrorY2()*3);
1641 t.SetErrorZ2(t.GetErrorZ2()*3);
1643 t.SetNFoundable(t.GetNFoundable()+1);
1644 UpdateTrack(&t,accept);
1647 else { // Remove old cluster from track
1648 t.SetClusterIndex(nr, -3);
1649 t.SetClusterPointer(nr, 0);
1653 if (TMath::Abs(t.GetSnp())>AliTPCReconstructor::GetMaxSnpTracker()) return 0; // cut on angle
1654 if (fIteration>1 && IsFindable(t)){
1655 // not look for new cluster during refitting
1656 t.SetNFoundable(t.GetNFoundable()+1);
1661 // if (TMath::Abs(t.GetSnp())>0.95 || TMath::Abs(x*t.GetC()-t.GetEta())>0.95) return 0;// patch 28 fev 06
1662 if (!t.PropagateTo(x)) {
1663 if (fIteration==0) t.SetRemoval(10);
1666 Double_t y = t.GetY();
1667 if (TMath::Abs(y)>ymax){
1669 t.SetRelativeSector((t.GetRelativeSector()+1) % fN);
1670 if (!t.Rotate(fSectors->GetAlpha()))
1672 } else if (y <-ymax) {
1673 t.SetRelativeSector((t.GetRelativeSector()-1+fN) % fN);
1674 if (!t.Rotate(-fSectors->GetAlpha()))
1677 if (!t.PropagateTo(x)) {
1678 if (fIteration==0) t.SetRemoval(10);
1684 Double_t z=t.GetZ();
1687 if (!IsActive(t.GetRelativeSector(),nr)) {
1689 t.SetClusterIndex2(nr,-1);
1692 //AliInfo(Form("A - Sector%d phi %f - alpha %f", t.fRelativeSector,y/x, t.GetAlpha()));
1693 Bool_t isActive = IsActive(t.GetRelativeSector(),nr);
1694 Bool_t isActive2 = (nr>=fInnerSec->GetNRows()) ? fOuterSec[t.GetRelativeSector()][nr-fInnerSec->GetNRows()].GetN()>0:fInnerSec[t.GetRelativeSector()][nr].GetN()>0;
1696 if (!isActive || !isActive2) return 0;
1698 const AliTPCtrackerRow &krow=GetRow(t.GetRelativeSector(),nr);
1699 if ( (t.GetSigmaY2()<0) || t.GetSigmaZ2()<0) return 0;
1701 Double_t roadz = 1.;
1703 if (TMath::Abs(TMath::Abs(y)-ymax)<krow.GetDeadZone()){
1705 t.SetClusterIndex2(nr,-1);
1711 // if (TMath::Abs(z)<(AliTPCReconstructor::GetCtgRange()*x+10) && TMath::Abs(z)<fkParam->GetZLength(0) && (TMath::Abs(t.GetSnp())<AliTPCReconstructor::GetMaxSnpTracker()))
1712 t.SetNFoundable(t.GetNFoundable()+1);
1718 // cl = krow.FindNearest2(y+10.,z,roady,roadz,index);
1719 cl = krow.FindNearest2(y,z,roady,roadz,index);
1720 if (cl) t.SetCurrentClusterIndex1(krow.GetIndex(index));
1723 t.SetCurrentCluster(cl);
1725 if (fIteration==2&&cl->IsUsed(10)) return 0;
1726 Int_t accept = AcceptCluster(&t,t.GetCurrentCluster());
1727 if (fIteration==2&&cl->IsUsed(11)) {
1728 t.SetErrorY2(t.GetErrorY2()+0.03);
1729 t.SetErrorZ2(t.GetErrorZ2()+0.03);
1730 t.SetErrorY2(t.GetErrorY2()*3);
1731 t.SetErrorZ2(t.GetErrorZ2()*3);
1734 if (t.fCurrentCluster->IsUsed(10)){
1739 if (t.fNShared>0.7*t.GetNumberOfClusters()) {
1745 if (accept<3) UpdateTrack(&t,accept);
1748 if ( fIteration==0 && t.GetNFoundable()*0.5 > t.GetNumberOfClusters()) t.SetRemoval(10);
1756 //_________________________________________________________________________
1757 Bool_t AliTPCtrackerMI::GetTrackPoint(Int_t index, AliTrackPoint &p ) const
1759 // Get track space point by index
1760 // return false in case the cluster doesn't exist
1761 AliTPCclusterMI *cl = GetClusterMI(index);
1762 if (!cl) return kFALSE;
1763 Int_t sector = (index&0xff000000)>>24;
1764 // Int_t row = (index&0x00ff0000)>>16;
1766 // xyz[0] = fkParam->GetPadRowRadii(sector,row);
1767 xyz[0] = cl->GetX();
1768 xyz[1] = cl->GetY();
1769 xyz[2] = cl->GetZ();
1771 fkParam->AdjustCosSin(sector,cos,sin);
1772 Float_t x = cos*xyz[0]-sin*xyz[1];
1773 Float_t y = cos*xyz[1]+sin*xyz[0];
1775 Float_t sigmaY2 = 0.027*cl->GetSigmaY2();
1776 if (sector < fkParam->GetNInnerSector()) sigmaY2 *= 2.07;
1777 Float_t sigmaZ2 = 0.066*cl->GetSigmaZ2();
1778 if (sector < fkParam->GetNInnerSector()) sigmaZ2 *= 1.77;
1779 cov[0] = sin*sin*sigmaY2;
1780 cov[1] = -sin*cos*sigmaY2;
1782 cov[3] = cos*cos*sigmaY2;
1785 p.SetXYZ(x,y,xyz[2],cov);
1786 AliGeomManager::ELayerID iLayer;
1788 if (sector < fkParam->GetNInnerSector()) {
1789 iLayer = AliGeomManager::kTPC1;
1793 iLayer = AliGeomManager::kTPC2;
1794 idet = sector - fkParam->GetNInnerSector();
1796 UShort_t volid = AliGeomManager::LayerToVolUID(iLayer,idet);
1797 p.SetVolumeID(volid);
1803 Int_t AliTPCtrackerMI::UpdateClusters(AliTPCseed& t, Int_t nr) {
1804 //-----------------------------------------------------------------
1805 // This function tries to find a track prolongation to next pad row
1806 //-----------------------------------------------------------------
1807 t.SetCurrentCluster(0);
1808 t.SetCurrentClusterIndex1(-3);
1810 Double_t xt=t.GetX();
1811 Int_t row = GetRowNumber(xt)-1;
1812 Double_t ymax= GetMaxY(nr);
1814 if (row < nr) return 1; // don't prolongate if not information until now -
1815 // if (TMath::Abs(t.GetSnp())>0.9 && t.GetNumberOfClusters()>40. && fIteration!=2) {
1817 // return 0; // not prolongate strongly inclined tracks
1819 // if (TMath::Abs(t.GetSnp())>0.95) {
1821 // return 0; // not prolongate strongly inclined tracks
1822 // }// patch 28 fev 06
1824 Double_t x= GetXrow(nr);
1826 //t.PropagateTo(x+0.02);
1827 //t.PropagateTo(x+0.01);
1828 if (!t.PropagateTo(x)){
1835 if (TMath::Abs(y)>ymax){
1837 t.SetRelativeSector((t.GetRelativeSector()+1) % fN);
1838 if (!t.Rotate(fSectors->GetAlpha()))
1840 } else if (y <-ymax) {
1841 t.SetRelativeSector((t.GetRelativeSector()-1+fN) % fN);
1842 if (!t.Rotate(-fSectors->GetAlpha()))
1845 // if (!t.PropagateTo(x)){
1852 if (TMath::Abs(t.GetSnp())>AliTPCReconstructor::GetMaxSnpTracker()) return 0;
1854 if (!IsActive(t.GetRelativeSector(),nr)) {
1856 t.SetClusterIndex2(nr,-1);
1859 //AliInfo(Form("A - Sector%d phi %f - alpha %f", t.fRelativeSector,y/x, t.GetAlpha()));
1861 AliTPCtrackerRow &krow=GetRow(t.GetRelativeSector(),nr);
1863 if (TMath::Abs(TMath::Abs(y)-ymax)<krow.GetDeadZone()){
1865 t.SetClusterIndex2(nr,-1);
1871 // if (TMath::Abs(t.GetZ())<(AliTPCReconstructor::GetCtgRange()*t.GetX()+10) && (TMath::Abs(t.GetSnp())<AliTPCReconstructor::GetMaxSnpTracker()))
1872 if (IsFindable(t)) t.SetNFoundable(t.GetNFoundable()+1);
1878 if ( (nr%2==0) || t.GetNumberOfClusters()<2){
1879 // t.fCurrentSigmaY = GetSigmaY(&t);
1880 //t.fCurrentSigmaZ = GetSigmaZ(&t);
1884 AliTPCclusterMI *cl=0;
1887 Double_t roady = 1.;
1888 Double_t roadz = 1.;
1892 index = t.GetClusterIndex2(nr);
1893 if ( (index >= 0) && (index&0x8000)==0){
1894 cl = t.GetClusterPointer(nr);
1895 if ( (cl==0) && (index >= 0)) cl = GetClusterMI(index);
1896 t.SetCurrentClusterIndex1(index);
1898 t.SetCurrentCluster(cl);
1904 // if (index<0) return 0;
1905 UInt_t uindex = TMath::Abs(index);
1908 //cl = krow.FindNearest2(y+10,z,roady,roadz,uindex);
1909 cl = krow.FindNearest2(y,z,roady,roadz,uindex);
1912 if (cl) t.SetCurrentClusterIndex1(krow.GetIndex(uindex));
1913 t.SetCurrentCluster(cl);
1919 Int_t AliTPCtrackerMI::FollowToNextCluster(AliTPCseed & t, Int_t nr) {
1920 //-----------------------------------------------------------------
1921 // This function tries to find a track prolongation to next pad row
1922 //-----------------------------------------------------------------
1924 //update error according neighborhoud
1926 if (t.GetCurrentCluster()) {
1928 Int_t accept = AcceptCluster(&t,t.GetCurrentCluster());
1930 if (t.GetCurrentCluster()->IsUsed(10)){
1935 t.SetNShared(t.GetNShared()+1);
1936 if (t.GetNShared()>0.7*t.GetNumberOfClusters()) {
1941 if (fIteration>0) accept = 0;
1942 if (accept<3) UpdateTrack(&t,accept);
1946 if ( t.GetNumberOfClusters()>18 && ( (t.GetSigmaY2()+t.GetSigmaZ2())>0.16)) t.SetRemoval(10);
1947 if ( t.GetNumberOfClusters()>18 && t.GetChi2()/t.GetNumberOfClusters()>6 ) t.SetRemoval(10);
1949 if (( (t.GetNFoundable()*0.5 > t.GetNumberOfClusters()) || t.GetNoCluster()>15)) t.SetRemoval(10);
1957 //_____________________________________________________________________________
1958 Int_t AliTPCtrackerMI::FollowProlongation(AliTPCseed& t, Int_t rf, Int_t step, Bool_t fromSeeds) {
1959 //-----------------------------------------------------------------
1960 // This function tries to find a track prolongation.
1961 //-----------------------------------------------------------------
1962 Double_t xt=t.GetX();
1964 Double_t alpha=t.GetAlpha();
1965 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
1966 if (alpha < 0. ) alpha += 2.*TMath::Pi();
1968 t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN);
1970 Int_t first = GetRowNumber(xt);
1975 for (Int_t nr= first; nr>=rf; nr-=step) {
1977 if (t.GetKinkIndexes()[0]>0){
1978 for (Int_t i=0;i<3;i++){
1979 Int_t index = t.GetKinkIndexes()[i];
1980 if (index==0) break;
1981 if (index<0) continue;
1983 AliKink * kink = (AliKink*)fEvent->GetKink(index-1);
1985 printf("PROBLEM\n");
1988 Int_t kinkrow = kink->GetTPCRow0()+2+Int_t(0.5/(0.05+kink->GetAngle(2)));
1990 AliExternalTrackParam paramd(t);
1991 kink->SetDaughter(paramd);
1992 kink->SetStatus(2,5);
1999 if (nr==80) t.UpdateReference();
2000 if (nr<fInnerSec->GetNRows())
2001 fSectors = fInnerSec;
2003 fSectors = fOuterSec;
2004 if (FollowToNext(t,nr)==0)
2017 Int_t AliTPCtrackerMI::FollowBackProlongation(AliTPCseed& t, Int_t rf, Bool_t fromSeeds) {
2018 //-----------------------------------------------------------------
2019 // This function tries to find a track prolongation.
2020 //-----------------------------------------------------------------
2022 Double_t xt=t.GetX();
2023 Double_t alpha=t.GetAlpha();
2024 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
2025 if (alpha < 0. ) alpha += 2.*TMath::Pi();
2026 t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN);
2028 Int_t first = t.GetFirstPoint();
2029 Int_t ri = GetRowNumber(xt);
2033 if (first<ri) first = ri;
2035 if (first<0) first=0;
2036 for (Int_t nr=first; nr<=rf; nr++) {
2037 // if ( (TMath::Abs(t.GetSnp())>0.95)) break;//patch 28 fev 06
2038 if (t.GetKinkIndexes()[0]<0){
2039 for (Int_t i=0;i<3;i++){
2040 Int_t index = t.GetKinkIndexes()[i];
2041 if (index==0) break;
2042 if (index>0) continue;
2043 index = TMath::Abs(index);
2044 AliKink * kink = (AliKink*)fEvent->GetKink(index-1);
2046 printf("PROBLEM\n");
2049 Int_t kinkrow = kink->GetTPCRow0()-2-Int_t(0.5/(0.05+kink->GetAngle(2)));
2051 AliExternalTrackParam paramm(t);
2052 kink->SetMother(paramm);
2053 kink->SetStatus(2,1);
2060 if (nr<fInnerSec->GetNRows())
2061 fSectors = fInnerSec;
2063 fSectors = fOuterSec;
2074 Float_t AliTPCtrackerMI::OverlapFactor(AliTPCseed * s1, AliTPCseed * s2, Int_t &sum1, Int_t & sum2)
2076 // overlapping factor
2082 Float_t dz2 =(s1->GetZ() - s2->GetZ());
2085 Float_t dy2 =TMath::Abs((s1->GetY() - s2->GetY()));
2087 Float_t distance = TMath::Sqrt(dz2+dy2);
2088 if (distance>4.) return 0; // if there are far away - not overlap - to reduce combinatorics
2091 Int_t firstpoint = TMath::Min(s1->GetFirstPoint(),s2->GetFirstPoint());
2092 Int_t lastpoint = TMath::Max(s1->GetLastPoint(),s2->GetLastPoint());
2097 if (firstpoint>lastpoint) {
2098 firstpoint =lastpoint;
2103 for (Int_t i=firstpoint-1;i<lastpoint+1;i++){
2104 if (s1->GetClusterIndex2(i)>0) sum1++;
2105 if (s2->GetClusterIndex2(i)>0) sum2++;
2106 if (s1->GetClusterIndex2(i)==s2->GetClusterIndex2(i) && s1->GetClusterIndex2(i)>0) {
2110 if (sum<5) return 0;
2112 Float_t summin = TMath::Min(sum1+1,sum2+1);
2113 Float_t ratio = (sum+1)/Float_t(summin);
2117 void AliTPCtrackerMI::SignShared(AliTPCseed * s1, AliTPCseed * s2)
2121 Float_t thetaCut = 0.2;//+10.*TMath::Sqrt(s1->GetSigmaTglZ()+ s2->GetSigmaTglZ());
2122 if (TMath::Abs(s1->GetTgl()-s2->GetTgl())>thetaCut) return;
2123 Float_t minCl = TMath::Min(s1->GetNumberOfClusters(),s2->GetNumberOfClusters());
2124 Int_t cutN0 = TMath::Max(5,TMath::Nint(0.1*minCl));
2129 //Int_t firstpoint = TMath::Max(s1->GetFirstPoint(),s2->GetFirstPoint());
2130 //Int_t lastpoint = TMath::Min(s1->GetLastPoint(),s2->GetLastPoint());
2131 Int_t firstpoint = 0;
2132 Int_t lastpoint = 160;
2134 // if (firstpoint>=lastpoint-5) return;;
2136 for (Int_t i=firstpoint;i<lastpoint;i++){
2137 // if ( (s1->GetClusterIndex2(i)&0xFFFF8FFF)==(s2->GetClusterIndex2(i)&0xFFFF8FFF) && s1->GetClusterIndex2(i)>0) {
2138 if ( (s1->GetClusterIndex2(i))==(s2->GetClusterIndex2(i)) && s1->GetClusterIndex2(i)>0) {
2142 if (sumshared>cutN0){
2145 for (Int_t i=firstpoint;i<lastpoint;i++){
2146 // if ( (s1->GetClusterIndex2(i)&0xFFFF8FFF)==(s2->GetClusterIndex2(i)&0xFFFF8FFF) && s1->GetClusterIndex2(i)>0) {
2147 if ( (s1->GetClusterIndex2(i))==(s2->GetClusterIndex2(i)) && s1->GetClusterIndex2(i)>0) {
2148 AliTPCTrackerPoint *p1 = s1->GetTrackPoint(i);
2149 AliTPCTrackerPoint *p2 = s2->GetTrackPoint(i);;
2150 if (s1->IsActive()&&s2->IsActive()){
2151 p1->SetShared(kTRUE);
2152 p2->SetShared(kTRUE);
2158 if (sumshared>cutN0){
2159 for (Int_t i=0;i<4;i++){
2160 if (s1->GetOverlapLabel(3*i)==0){
2161 s1->SetOverlapLabel(3*i, s2->GetLabel());
2162 s1->SetOverlapLabel(3*i+1,sumshared);
2163 s1->SetOverlapLabel(3*i+2,s2->GetUniqueID());
2167 for (Int_t i=0;i<4;i++){
2168 if (s2->GetOverlapLabel(3*i)==0){
2169 s2->SetOverlapLabel(3*i, s1->GetLabel());
2170 s2->SetOverlapLabel(3*i+1,sumshared);
2171 s2->SetOverlapLabel(3*i+2,s1->GetUniqueID());
2178 void AliTPCtrackerMI::SignShared(TObjArray * arr)
2181 //sort trackss according sectors
2183 for (Int_t i=0; i<arr->GetEntriesFast(); i++) {
2184 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2186 //if (pt) RotateToLocal(pt);
2190 arr->Sort(); // sorting according relative sectors
2191 arr->Expand(arr->GetEntries());
2194 Int_t nseed=arr->GetEntriesFast();
2195 for (Int_t i=0; i<nseed; i++) {
2196 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2198 for (Int_t j=0;j<12;j++){
2199 pt->SetOverlapLabel(j,0);
2202 for (Int_t i=0; i<nseed; i++) {
2203 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2205 if (pt->GetRemoval()>10) continue;
2206 for (Int_t j=i+1; j<nseed; j++){
2207 AliTPCseed *pt2=(AliTPCseed*)arr->UncheckedAt(j);
2208 if (TMath::Abs(pt->GetRelativeSector()-pt2->GetRelativeSector())>1) continue;
2210 if (pt2->GetRemoval()<=10) {
2211 //if ( TMath::Abs(pt->GetRelativeSector()-pt2->GetRelativeSector())>0) break;
2219 void AliTPCtrackerMI::SortTracks(TObjArray * arr, Int_t mode) const
2222 //sort tracks in array according mode criteria
2223 Int_t nseed = arr->GetEntriesFast();
2224 for (Int_t i=0; i<nseed; i++) {
2225 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2236 void AliTPCtrackerMI::RemoveUsed2(TObjArray * arr, Float_t factor1, Float_t factor2, Int_t minimal)
2239 // Loop over all tracks and remove overlaped tracks (with lower quality)
2241 // 1. Unsign clusters
2242 // 2. Sort tracks according quality
2243 // Quality is defined by the number of cluster between first and last points
2245 // 3. Loop over tracks - decreasing quality order
2246 // a.) remove - If the fraction of shared cluster less than factor (1- n or 2)
2247 // b.) remove - If the minimal number of clusters less than minimal and not ITS
2248 // c.) if track accepted - sign clusters
2250 //Called in - AliTPCtrackerMI::Clusters2Tracks()
2251 // - AliTPCtrackerMI::PropagateBack()
2252 // - AliTPCtrackerMI::RefitInward()
2255 // factor1 - factor for constrained
2256 // factor2 - for non constrained tracks
2257 // if (Float_t(shared+1)/Float_t(found+1)>factor) - DELETE
2261 Int_t nseed = arr->GetEntriesFast();
2262 Float_t * quality = new Float_t[nseed];
2263 Int_t * indexes = new Int_t[nseed];
2267 for (Int_t i=0; i<nseed; i++) {
2268 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2273 pt->UpdatePoints(); //select first last max dens points
2274 Float_t * points = pt->GetPoints();
2275 if (points[3]<0.8) quality[i] =-1;
2276 quality[i] = (points[2]-points[0])+pt->GetNumberOfClusters();
2277 //prefer high momenta tracks if overlaps
2278 quality[i] *= TMath::Sqrt(TMath::Abs(pt->Pt())+0.5);
2280 TMath::Sort(nseed,quality,indexes);
2283 for (Int_t itrack=0; itrack<nseed; itrack++) {
2284 Int_t trackindex = indexes[itrack];
2285 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(trackindex);
2288 if (quality[trackindex]<0){
2289 delete arr->RemoveAt(trackindex);
2294 Int_t first = Int_t(pt->GetPoints()[0]);
2295 Int_t last = Int_t(pt->GetPoints()[2]);
2296 Double_t factor = (pt->GetBConstrain()) ? factor1: factor2;
2298 Int_t found,foundable,shared;
2299 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
2300 // pt->GetClusterStatistic(0,160, found, foundable,shared,kFALSE);
2301 Bool_t itsgold =kFALSE;
2304 if (pt->GetESD()->GetITSclusters(dummy)>4) itsgold= kTRUE;
2308 if (Float_t(shared+1)/Float_t(found+1)>factor){
2309 if (pt->GetKinkIndexes()[0]!=0) continue; //don't remove tracks - part of the kinks
2310 if( AliTPCReconstructor::StreamLevel()>3){
2311 TTreeSRedirector &cstream = *fDebugStreamer;
2312 cstream<<"RemoveUsed"<<
2313 "iter="<<fIteration<<
2317 delete arr->RemoveAt(trackindex);
2320 if (pt->GetNumberOfClusters()<50&&(found-0.5*shared)<minimal){ //remove short tracks
2321 if (pt->GetKinkIndexes()[0]!=0) continue; //don't remove tracks - part of the kinks
2322 if( AliTPCReconstructor::StreamLevel()>3){
2323 TTreeSRedirector &cstream = *fDebugStreamer;
2324 cstream<<"RemoveShort"<<
2325 "iter="<<fIteration<<
2329 delete arr->RemoveAt(trackindex);
2335 //if (sharedfactor>0.4) continue;
2336 if (pt->GetKinkIndexes()[0]>0) continue;
2337 //Remove tracks with undefined properties - seems
2338 if (pt->GetSigmaY2()<kAlmost0) continue; // ? what is the origin ?
2340 for (Int_t i=first; i<last; i++) {
2341 Int_t index=pt->GetClusterIndex2(i);
2342 // if (index<0 || index&0x8000 ) continue;
2343 if (index<0 || index&0x8000 ) continue;
2344 AliTPCclusterMI *c= pt->GetClusterPointer(i);
2351 Info("RemoveUsed","\n*****\nNumber of good tracks after shared removal\t%d\n",fNtracks);
2357 void AliTPCtrackerMI::DumpClusters(Int_t iter, TObjArray *trackArray)
2360 // Dump clusters after reco
2361 // signed and unsigned cluster can be visualized
2362 // 1. Unsign all cluster
2363 // 2. Sign all used clusters
2366 Int_t nseed = trackArray->GetEntries();
2367 for (Int_t i=0; i<nseed; i++){
2368 AliTPCseed *pt=(AliTPCseed*)trackArray->UncheckedAt(i);
2372 Bool_t isKink=pt->GetKinkIndex(0)!=0;
2373 for (Int_t j=0; j<160; ++j) {
2374 Int_t index=pt->GetClusterIndex2(j);
2375 if (index<0) continue;
2376 AliTPCclusterMI *c= pt->GetClusterPointer(j);
2378 if (isKink) c->Use(100); // kink
2379 c->Use(10); // by default usage 10
2384 for (Int_t sec=0;sec<fkNIS;sec++){
2385 for (Int_t row=0;row<fInnerSec->GetNRows();row++){
2386 AliTPCclusterMI *cl = fInnerSec[sec][row].GetClusters1();
2387 for (Int_t icl =0;icl< fInnerSec[sec][row].GetN1();icl++){
2388 Float_t gx[3]; cl[icl].GetGlobalXYZ(gx);
2389 (*fDebugStreamer)<<"clDump"<<
2397 cl = fInnerSec[sec][row].GetClusters2();
2398 for (Int_t icl =0;icl< fInnerSec[sec][row].GetN2();icl++){
2399 Float_t gx[3]; cl[icl].GetGlobalXYZ(gx);
2400 (*fDebugStreamer)<<"clDump"<<
2411 for (Int_t sec=0;sec<fkNOS;sec++){
2412 for (Int_t row=0;row<fOuterSec->GetNRows();row++){
2413 AliTPCclusterMI *cl = fOuterSec[sec][row].GetClusters1();
2414 for (Int_t icl =0;icl< fOuterSec[sec][row].GetN1();icl++){
2415 Float_t gx[3]; cl[icl].GetGlobalXYZ(gx);
2416 (*fDebugStreamer)<<"clDump"<<
2424 cl = fOuterSec[sec][row].GetClusters2();
2425 for (Int_t icl =0;icl< fOuterSec[sec][row].GetN2();icl++){
2426 Float_t gx[3]; cl[icl].GetGlobalXYZ(gx);
2427 (*fDebugStreamer)<<"clDump"<<
2439 void AliTPCtrackerMI::UnsignClusters()
2442 // loop over all clusters and unsign them
2445 for (Int_t sec=0;sec<fkNIS;sec++){
2446 for (Int_t row=0;row<fInnerSec->GetNRows();row++){
2447 AliTPCclusterMI *cl = fInnerSec[sec][row].GetClusters1();
2448 for (Int_t icl =0;icl< fInnerSec[sec][row].GetN1();icl++)
2449 // if (cl[icl].IsUsed(10))
2451 cl = fInnerSec[sec][row].GetClusters2();
2452 for (Int_t icl =0;icl< fInnerSec[sec][row].GetN2();icl++)
2453 //if (cl[icl].IsUsed(10))
2458 for (Int_t sec=0;sec<fkNOS;sec++){
2459 for (Int_t row=0;row<fOuterSec->GetNRows();row++){
2460 AliTPCclusterMI *cl = fOuterSec[sec][row].GetClusters1();
2461 for (Int_t icl =0;icl< fOuterSec[sec][row].GetN1();icl++)
2462 //if (cl[icl].IsUsed(10))
2464 cl = fOuterSec[sec][row].GetClusters2();
2465 for (Int_t icl =0;icl< fOuterSec[sec][row].GetN2();icl++)
2466 //if (cl[icl].IsUsed(10))
2475 void AliTPCtrackerMI::SignClusters(const TObjArray * arr, Float_t fnumber, Float_t fdensity)
2478 //sign clusters to be "used"
2480 // snumber and sdensity sign number of sigmas - bellow mean value to be accepted
2481 // loop over "primaries"
2495 Int_t nseed = arr->GetEntriesFast();
2496 for (Int_t i=0; i<nseed; i++) {
2497 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2501 if (!(pt->IsActive())) continue;
2502 Float_t dens = pt->GetNumberOfClusters()/Float_t(pt->GetNFoundable());
2503 if ( (dens>0.7) && (pt->GetNumberOfClusters()>70)){
2505 sumdens2+= dens*dens;
2506 sumn += pt->GetNumberOfClusters();
2507 sumn2 += pt->GetNumberOfClusters()*pt->GetNumberOfClusters();
2508 Float_t chi2 = pt->GetChi2()/pt->GetNumberOfClusters();
2511 sumchi2 +=chi2*chi2;
2516 Float_t mdensity = 0.9;
2517 Float_t meann = 130;
2518 Float_t meanchi = 1;
2519 Float_t sdensity = 0.1;
2520 Float_t smeann = 10;
2521 Float_t smeanchi =0.4;
2525 mdensity = sumdens/sum;
2527 meanchi = sumchi/sum;
2529 sdensity = sumdens2/sum-mdensity*mdensity;
2531 sdensity = TMath::Sqrt(sdensity);
2535 smeann = sumn2/sum-meann*meann;
2537 smeann = TMath::Sqrt(smeann);
2541 smeanchi = sumchi2/sum - meanchi*meanchi;
2543 smeanchi = TMath::Sqrt(smeanchi);
2549 //REMOVE SHORT DELTAS or tracks going out of sensitive volume of TPC
2551 for (Int_t i=0; i<nseed; i++) {
2552 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2556 if (pt->GetBSigned()) continue;
2557 if (pt->GetBConstrain()) continue;
2558 //if (!(pt->IsActive())) continue;
2560 Int_t found,foundable,shared;
2561 pt->GetClusterStatistic(0,160,found, foundable,shared);
2562 if (shared/float(found)>0.3) {
2563 if (shared/float(found)>0.9 ){
2564 //delete arr->RemoveAt(i);
2569 Bool_t isok =kFALSE;
2570 if ( (pt->GetNShared()/pt->GetNumberOfClusters()<0.5) &&pt->GetNumberOfClusters()>60)
2572 if ((TMath::Abs(1/pt->GetC())<100.) && (pt->GetNShared()/pt->GetNumberOfClusters()<0.7))
2574 if (TMath::Abs(pt->GetZ()/pt->GetX())>1.1)
2576 if ( (TMath::Abs(pt->GetSnp()>0.7) && pt->GetD(0,0)>60.))
2580 for (Int_t j=0; j<160; ++j) {
2581 Int_t index=pt->GetClusterIndex2(j);
2582 if (index<0) continue;
2583 AliTPCclusterMI *c= pt->GetClusterPointer(j);
2585 //if (!(c->IsUsed(10))) c->Use();
2592 Double_t maxchi = meanchi+2.*smeanchi;
2594 for (Int_t i=0; i<nseed; i++) {
2595 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2599 //if (!(pt->IsActive())) continue;
2600 if (pt->GetBSigned()) continue;
2601 Double_t chi = pt->GetChi2()/pt->GetNumberOfClusters();
2602 if (chi>maxchi) continue;
2605 Float_t dens = pt->GetNumberOfClusters()/Float_t(pt->GetNFoundable());
2607 //sign only tracks with enoug big density at the beginning
2609 if ((pt->GetDensityFirst(40)<0.75) && pt->GetNumberOfClusters()<meann) continue;
2612 Double_t mindens = TMath::Max(double(mdensity-sdensity*fdensity*bfactor),0.65);
2613 Double_t minn = TMath::Max(Int_t(meann-fnumber*smeann*bfactor),50);
2615 // if (pt->fBConstrain) mindens = TMath::Max(mdensity-sdensity*fdensity*bfactor,0.65);
2616 if ( (pt->GetRemoval()==10) && (pt->GetSnp()>0.8)&&(dens>mindens))
2619 if ((dens>mindens && pt->GetNumberOfClusters()>minn) && chi<maxchi ){
2620 //Int_t noc=pt->GetNumberOfClusters();
2621 pt->SetBSigned(kTRUE);
2622 for (Int_t j=0; j<160; ++j) {
2624 Int_t index=pt->GetClusterIndex2(j);
2625 if (index<0) continue;
2626 AliTPCclusterMI *c= pt->GetClusterPointer(j);
2628 // if (!(c->IsUsed(10))) c->Use();
2633 // gLastCheck = nseed;
2641 void AliTPCtrackerMI::StopNotActive(const TObjArray * arr, Int_t row0, Float_t th0, Float_t th1, Float_t th2) const
2643 // stop not active tracks
2644 // take th1 as threshold for number of founded to number of foundable on last 10 active rows
2645 // take th2 as threshold for number of founded to number of foundable on last 20 active rows
2646 Int_t nseed = arr->GetEntriesFast();
2648 for (Int_t i=0; i<nseed; i++) {
2649 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2653 if (!(pt->IsActive())) continue;
2654 StopNotActive(pt,row0,th0, th1,th2);
2660 void AliTPCtrackerMI::StopNotActive(AliTPCseed * seed, Int_t row0, Float_t th0, Float_t th1,
2663 // stop not active tracks
2664 // take th1 as threshold for number of founded to number of foundable on last 10 active rows
2665 // take th2 as threshold for number of founded to number of foundable on last 20 active rows
2668 Int_t foundable = 0;
2669 Int_t maxindex = seed->GetLastPoint(); //last foundable row
2670 if (seed->GetNFoundable()*th0 > seed->GetNumberOfClusters()) {
2671 seed->Desactivate(10) ;
2675 for (Int_t i=row0; i<maxindex; i++){
2676 Int_t index = seed->GetClusterIndex2(i);
2677 if (index!=-1) foundable++;
2679 if (foundable<=30) sumgood1++;
2680 if (foundable<=50) {
2687 if (foundable>=30.){
2688 if (sumgood1<(th1*30.)) seed->Desactivate(10);
2691 if (sumgood2<(th2*50.)) seed->Desactivate(10);
2695 Int_t AliTPCtrackerMI::RefitInward(AliESDEvent *event)
2698 // back propagation of ESD tracks
2701 if (!event) return 0;
2702 const Int_t kMaxFriendTracks=2000;
2704 // extract correction object for multiplicity dependence of dEdx
2705 TObjArray * gainCalibArray = AliTPCcalibDB::Instance()->GetTimeGainSplinesRun(event->GetRunNumber());
2707 AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ;
2709 AliFatal("Tranformations not in RefitInward");
2712 transform->SetCurrentRecoParam((AliTPCRecoParam*)AliTPCReconstructor::GetRecoParam());
2713 const AliTPCRecoParam * recoParam = AliTPCcalibDB::Instance()->GetTransform()->GetCurrentRecoParam();
2714 Int_t nContribut = event->GetNumberOfTracks();
2715 TGraphErrors * graphMultDependenceDeDx = 0x0;
2716 if (recoParam && recoParam->GetUseMultiplicityCorrectionDedx() && gainCalibArray) {
2717 if (recoParam->GetUseTotCharge()) {
2718 graphMultDependenceDeDx = (TGraphErrors *) gainCalibArray->FindObject("TGRAPHERRORS_MEANQTOT_MULTIPLICITYDEPENDENCE_BEAM_ALL");
2720 graphMultDependenceDeDx = (TGraphErrors *) gainCalibArray->FindObject("TGRAPHERRORS_MEANQMAX_MULTIPLICITYDEPENDENCE_BEAM_ALL");
2726 //PrepareForProlongation(fSeeds,1);
2727 PropagateForward2(fSeeds);
2728 RemoveUsed2(fSeeds,0.4,0.4,20);
2730 TObjArray arraySeed(fSeeds->GetEntries());
2731 for (Int_t i=0;i<fSeeds->GetEntries();i++) {
2732 arraySeed.AddAt(fSeeds->At(i),i);
2734 SignShared(&arraySeed);
2735 // FindCurling(fSeeds, event,2); // find multi found tracks
2736 FindSplitted(fSeeds, event,2); // find multi found tracks
2737 if (AliTPCReconstructor::StreamLevel()>5) FindMultiMC(fSeeds, fEvent,2); // find multi found tracks
2740 Int_t nseed = fSeeds->GetEntriesFast();
2741 for (Int_t i=0;i<nseed;i++){
2742 AliTPCseed * seed = (AliTPCseed*) fSeeds->UncheckedAt(i);
2743 if (!seed) continue;
2744 if (seed->GetKinkIndex(0)>0) UpdateKinkQualityD(seed); // update quality informations for kinks
2745 AliESDtrack *esd=event->GetTrack(i);
2747 if (seed->GetNumberOfClusters()<60 && seed->GetNumberOfClusters()<(esd->GetTPCclusters(0) -5)*0.8){
2748 AliExternalTrackParam paramIn;
2749 AliExternalTrackParam paramOut;
2750 Int_t ncl = seed->RefitTrack(seed,¶mIn,¶mOut);
2751 if (AliTPCReconstructor::StreamLevel()>2) {
2752 (*fDebugStreamer)<<"RecoverIn"<<
2756 "pout.="<<¶mOut<<
2761 seed->Set(paramIn.GetX(),paramIn.GetAlpha(),paramIn.GetParameter(),paramIn.GetCovariance());
2762 seed->SetNumberOfClusters(ncl);
2766 seed->PropagateTo(fkParam->GetInnerRadiusLow());
2767 seed->UpdatePoints();
2768 AddCovariance(seed);
2769 MakeESDBitmaps(seed, esd);
2770 seed->CookdEdx(0.02,0.6);
2771 CookLabel(seed,0.1); //For comparison only
2773 if (AliTPCReconstructor::StreamLevel()>1 && seed!=0) {
2774 TTreeSRedirector &cstream = *fDebugStreamer;
2781 if (seed->GetNumberOfClusters()>15){
2782 esd->UpdateTrackParams(seed,AliESDtrack::kTPCrefit);
2783 esd->SetTPCPoints(seed->GetPoints());
2784 esd->SetTPCPointsF(seed->GetNFoundable());
2785 Int_t ndedx = seed->GetNCDEDX(0);
2786 Float_t sdedx = seed->GetSDEDX(0);
2787 Float_t dedx = seed->GetdEdx();
2788 // apply mutliplicity dependent dEdx correction if available
2789 if (graphMultDependenceDeDx) {
2790 Double_t corrGain = AliTPCcalibDButil::EvalGraphConst(graphMultDependenceDeDx, nContribut);
2791 dedx += (1 - corrGain)*50.; // MIP is normalized to 50
2793 esd->SetTPCsignal(dedx, sdedx, ndedx);
2795 // fill new dEdx information
2797 Double32_t signal[4];
2801 for(Int_t iarr=0;iarr<3;iarr++) {
2802 signal[iarr] = seed->GetDEDXregion(iarr+1);
2803 ncl[iarr] = seed->GetNCDEDX(iarr+1);
2804 nrows[iarr] = seed->GetNCDEDXInclThres(iarr+1);
2806 signal[3] = seed->GetDEDXregion(4);
2808 AliTPCdEdxInfo * infoTpcPid = new AliTPCdEdxInfo();
2809 infoTpcPid->SetTPCSignalRegionInfo(signal, ncl, nrows);
2810 esd->SetTPCdEdxInfo(infoTpcPid);
2812 // add seed to the esd track in Calib level
2814 Bool_t storeFriend = gRandom->Rndm()<(kMaxFriendTracks)/Float_t(nseed);
2815 if (AliTPCReconstructor::StreamLevel()>0 &&storeFriend){
2816 AliTPCseed * seedCopy = new AliTPCseed(*seed, kTRUE);
2817 esd->AddCalibObject(seedCopy);
2822 //printf("problem\n");
2825 //FindKinks(fSeeds,event);
2826 if (AliTPCReconstructor::StreamLevel()>3) DumpClusters(2,fSeeds);
2827 Info("RefitInward","Number of refitted tracks %d",ntracks);
2832 Int_t AliTPCtrackerMI::PropagateBack(AliESDEvent *event)
2835 // back propagation of ESD tracks
2837 if (!event) return 0;
2841 PropagateBack(fSeeds);
2842 RemoveUsed2(fSeeds,0.4,0.4,20);
2843 //FindCurling(fSeeds, fEvent,1);
2844 FindSplitted(fSeeds, event,1); // find multi found tracks
2845 if (AliTPCReconstructor::StreamLevel()>5) FindMultiMC(fSeeds, fEvent,1); // find multi found tracks
2848 Int_t nseed = fSeeds->GetEntriesFast();
2850 for (Int_t i=0;i<nseed;i++){
2851 AliTPCseed * seed = (AliTPCseed*) fSeeds->UncheckedAt(i);
2852 if (!seed) continue;
2853 if (seed->GetKinkIndex(0)<0) UpdateKinkQualityM(seed); // update quality informations for kinks
2854 seed->UpdatePoints();
2855 AddCovariance(seed);
2856 AliESDtrack *esd=event->GetTrack(i);
2857 if (!esd) continue; //never happen
2858 if (seed->GetNumberOfClusters()<60 && seed->GetNumberOfClusters()<(esd->GetTPCclusters(0) -5)*0.8){
2859 AliExternalTrackParam paramIn;
2860 AliExternalTrackParam paramOut;
2861 Int_t ncl = seed->RefitTrack(seed,¶mIn,¶mOut);
2862 if (AliTPCReconstructor::StreamLevel()>2) {
2863 (*fDebugStreamer)<<"RecoverBack"<<
2867 "pout.="<<¶mOut<<
2872 seed->Set(paramOut.GetX(),paramOut.GetAlpha(),paramOut.GetParameter(),paramOut.GetCovariance());
2873 seed->SetNumberOfClusters(ncl);
2876 seed->CookdEdx(0.02,0.6);
2877 CookLabel(seed,0.1); //For comparison only
2878 if (seed->GetNumberOfClusters()>15){
2879 esd->UpdateTrackParams(seed,AliESDtrack::kTPCout);
2880 esd->SetTPCPoints(seed->GetPoints());
2881 esd->SetTPCPointsF(seed->GetNFoundable());
2882 Int_t ndedx = seed->GetNCDEDX(0);
2883 Float_t sdedx = seed->GetSDEDX(0);
2884 Float_t dedx = seed->GetdEdx();
2885 esd->SetTPCsignal(dedx, sdedx, ndedx);
2887 Int_t eventnumber = event->GetEventNumberInFile();// patch 28 fev 06
2888 // This is most likely NOT the event number you'd like to use. It has nothing to do with the 'real' event number
2889 if (AliTPCReconstructor::StreamLevel()>1 && esd) {
2890 (*fDebugStreamer)<<"Cback"<<
2893 "EventNrInFile="<<eventnumber<<
2898 if (AliTPCReconstructor::StreamLevel()>3) DumpClusters(1,fSeeds);
2899 //FindKinks(fSeeds,event);
2900 Info("PropagateBack","Number of back propagated tracks %d",ntracks);
2907 void AliTPCtrackerMI::DeleteSeeds()
2912 Int_t nseed = fSeeds->GetEntriesFast();
2913 for (Int_t i=0;i<nseed;i++){
2914 AliTPCseed * seed = (AliTPCseed*)fSeeds->At(i);
2915 if (seed) delete fSeeds->RemoveAt(i);
2922 void AliTPCtrackerMI::ReadSeeds(const AliESDEvent *const event, Int_t direction)
2925 //read seeds from the event
2927 Int_t nentr=event->GetNumberOfTracks();
2929 Info("ReadSeeds", "Number of ESD tracks: %d\n", nentr);
2934 fSeeds = new TObjArray(nentr);
2938 for (Int_t i=0; i<nentr; i++) {
2939 AliESDtrack *esd=event->GetTrack(i);
2940 ULong_t status=esd->GetStatus();
2941 if (!(status&AliESDtrack::kTPCin)) continue;
2942 AliTPCtrack t(*esd);
2943 t.SetNumberOfClusters(0);
2944 // AliTPCseed *seed = new AliTPCseed(t,t.GetAlpha());
2945 AliTPCseed *seed = new AliTPCseed(t/*,t.GetAlpha()*/);
2946 seed->SetUniqueID(esd->GetID());
2947 AddCovariance(seed); //add systematic ucertainty
2948 for (Int_t ikink=0;ikink<3;ikink++) {
2949 Int_t index = esd->GetKinkIndex(ikink);
2950 seed->GetKinkIndexes()[ikink] = index;
2951 if (index==0) continue;
2952 index = TMath::Abs(index);
2953 AliESDkink * kink = fEvent->GetKink(index-1);
2954 if (kink&&esd->GetKinkIndex(ikink)<0){
2955 if ((status & AliESDtrack::kTRDrefit) != 0) kink->SetStatus(1,2);
2956 if ((status & AliESDtrack::kITSout) != 0) kink->SetStatus(1,0);
2958 if (kink&&esd->GetKinkIndex(ikink)>0){
2959 if ((status & AliESDtrack::kTRDrefit) != 0) kink->SetStatus(1,6);
2960 if ((status & AliESDtrack::kITSout) != 0) kink->SetStatus(1,4);
2964 if (((status&AliESDtrack::kITSout)==0)&&(direction==1)) seed->ResetCovariance(10.);
2965 if ( direction ==2 &&(status & AliESDtrack::kTRDrefit) == 0 ) seed->ResetCovariance(10.);
2966 //if ( direction ==2 && ((status & AliESDtrack::kTPCout) == 0) ) {
2967 // fSeeds->AddAt(0,i);
2971 if ( direction ==2 &&(status & AliESDtrack::kTRDrefit) > 0 ) {
2972 Double_t par0[5],par1[5],alpha,x;
2973 esd->GetInnerExternalParameters(alpha,x,par0);
2974 esd->GetExternalParameters(x,par1);
2975 Double_t delta1 = TMath::Abs(par0[4]-par1[4])/(0.000000001+TMath::Abs(par0[4]+par1[4]));
2976 Double_t delta2 = TMath::Abs(par0[3]-par1[3]);
2978 if (esd->GetTRDncls()>0) trdchi2 = esd->GetTRDchi2()/esd->GetTRDncls();
2979 //reset covariance if suspicious
2980 if ( (delta1>0.1) || (delta2>0.006) ||trdchi2>7.)
2981 seed->ResetCovariance(10.);
2986 // rotate to the local coordinate system
2988 fSectors=fInnerSec; fN=fkNIS;
2989 Double_t alpha=seed->GetAlpha();
2990 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
2991 if (alpha < 0. ) alpha += 2.*TMath::Pi();
2992 Int_t ns=Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN;
2993 alpha =ns*fSectors->GetAlpha() + fSectors->GetAlphaShift();
2994 alpha-=seed->GetAlpha();
2995 if (alpha<-TMath::Pi()) alpha += 2*TMath::Pi();
2996 if (alpha>=TMath::Pi()) alpha -= 2*TMath::Pi();
2997 if (TMath::Abs(alpha) > 0.001) { // This should not happen normally
2998 AliWarning(Form("Rotating track over %f",alpha));
2999 if (!seed->Rotate(alpha)) {
3006 if (esd->GetKinkIndex(0)<=0){
3007 for (Int_t irow=0;irow<160;irow++){
3008 Int_t index = seed->GetClusterIndex2(irow);
3011 AliTPCclusterMI * cl = GetClusterMI(index);
3012 seed->SetClusterPointer(irow,cl);
3014 if ((index & 0x8000)==0){
3015 cl->Use(10); // accepted cluster
3017 cl->Use(6); // close cluster not accepted
3020 Info("ReadSeeds","Not found cluster");
3025 fSeeds->AddAt(seed,i);
3031 //_____________________________________________________________________________
3032 void AliTPCtrackerMI::MakeSeeds3(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, Float_t cuts[4],
3033 Float_t deltay, Int_t ddsec) {
3034 //-----------------------------------------------------------------
3035 // This function creates track seeds.
3036 // SEEDING WITH VERTEX CONSTRAIN
3037 //-----------------------------------------------------------------
3038 // cuts[0] - fP4 cut
3039 // cuts[1] - tan(phi) cut
3040 // cuts[2] - zvertex cut
3041 // cuts[3] - fP3 cut
3049 Double_t x[5], c[15];
3050 // Int_t di = i1-i2;
3052 AliTPCseed * seed = new AliTPCseed();
3053 Double_t alpha=fSectors->GetAlpha(), shift=fSectors->GetAlphaShift();
3054 Double_t cs=cos(alpha), sn=sin(alpha);
3056 // Double_t x1 =fOuterSec->GetX(i1);
3057 //Double_t xx2=fOuterSec->GetX(i2);
3059 Double_t x1 =GetXrow(i1);
3060 Double_t xx2=GetXrow(i2);
3062 Double_t x3=GetX(), y3=GetY(), z3=GetZ();
3064 Int_t imiddle = (i2+i1)/2; //middle pad row index
3065 Double_t xm = GetXrow(imiddle); // radius of middle pad-row
3066 const AliTPCtrackerRow& krm=GetRow(sec,imiddle); //middle pad -row
3070 const AliTPCtrackerRow& kr1=GetRow(ns,i1);
3071 Double_t ymax = GetMaxY(i1)-kr1.GetDeadZone()-1.5;
3072 Double_t ymaxm = GetMaxY(imiddle)-kr1.GetDeadZone()-1.5;
3075 // change cut on curvature if it can't reach this layer
3076 // maximal curvature set to reach it
3077 Double_t dvertexmax = TMath::Sqrt((x1-x3)*(x1-x3)+(ymax+5-y3)*(ymax+5-y3));
3078 if (dvertexmax*0.5*cuts[0]>0.85){
3079 cuts[0] = 0.85/(dvertexmax*0.5+1.);
3081 Double_t r2min = 1/(cuts[0]*cuts[0]); //minimal square of radius given by cut
3084 if (deltay>0) ddsec = 0;
3085 // loop over clusters
3086 for (Int_t is=0; is < kr1; is++) {
3088 if (kr1[is]->IsUsed(10)) continue;
3089 Double_t y1=kr1[is]->GetY(), z1=kr1[is]->GetZ();
3090 //if (TMath::Abs(y1)>ymax) continue;
3092 if (deltay>0 && TMath::Abs(ymax-TMath::Abs(y1))> deltay ) continue; // seed only at the edge
3094 // find possible directions
3095 Float_t anglez = (z1-z3)/(x1-x3);
3096 Float_t extraz = z1 - anglez*(x1-xx2); // extrapolated z
3099 //find rotation angles relative to line given by vertex and point 1
3100 Double_t dvertex2 = (x1-x3)*(x1-x3)+(y1-y3)*(y1-y3);
3101 Double_t dvertex = TMath::Sqrt(dvertex2);
3102 Double_t angle13 = TMath::ATan((y1-y3)/(x1-x3));
3103 Double_t cs13 = cos(-angle13), sn13 = sin(-angle13);
3106 // loop over 2 sectors
3112 Double_t dddz1=0; // direction of delta inclination in z axis
3119 for (Int_t dsec = dsec1; dsec<=dsec2;dsec++){
3120 Int_t sec2 = sec + dsec;
3122 // AliTPCtrackerRow& kr2 = fOuterSec[(sec2+fkNOS)%fkNOS][i2];
3123 //AliTPCtrackerRow& kr2m = fOuterSec[(sec2+fkNOS)%fkNOS][imiddle];
3124 AliTPCtrackerRow& kr2 = GetRow((sec2+fkNOS)%fkNOS,i2);
3125 AliTPCtrackerRow& kr2m = GetRow((sec2+fkNOS)%fkNOS,imiddle);
3126 Int_t index1 = TMath::Max(kr2.Find(extraz-0.6-dddz1*TMath::Abs(z1)*0.05)-1,0);
3127 Int_t index2 = TMath::Min(kr2.Find(extraz+0.6+dddz2*TMath::Abs(z1)*0.05)+1,kr2);
3129 // rotation angles to p1-p3
3130 Double_t cs13r = cos(-angle13+dsec*alpha)/dvertex, sn13r = sin(-angle13+dsec*alpha)/dvertex;
3131 Double_t x2, y2, z2;
3133 // Double_t dymax = maxangle*TMath::Abs(x1-xx2);
3136 Double_t dxx0 = (xx2-x3)*cs13r;
3137 Double_t dyy0 = (xx2-x3)*sn13r;
3138 for (Int_t js=index1; js < index2; js++) {
3139 const AliTPCclusterMI *kcl = kr2[js];
3140 if (kcl->IsUsed(10)) continue;
3142 //calcutate parameters
3144 Double_t yy0 = dyy0 +(kcl->GetY()-y3)*cs13r;
3146 if (TMath::Abs(yy0)<0.000001) continue;
3147 Double_t xx0 = dxx0 -(kcl->GetY()-y3)*sn13r;
3148 Double_t y0 = 0.5*(xx0*xx0+yy0*yy0-xx0)/yy0;
3149 Double_t r02 = (0.25+y0*y0)*dvertex2;
3150 //curvature (radius) cut
3151 if (r02<r2min) continue;
3155 Double_t c0 = 1/TMath::Sqrt(r02);
3159 //Double_t dfi0 = 2.*TMath::ASin(dvertex*c0*0.5);
3160 //Double_t dfi1 = 2.*TMath::ASin(TMath::Sqrt(yy0*yy0+(1-xx0)*(1-xx0))*dvertex*c0*0.5);
3161 Double_t dfi0 = 2.*AliTPCFastMath::FastAsin(dvertex*c0*0.5);
3162 Double_t dfi1 = 2.*AliTPCFastMath::FastAsin(TMath::Sqrt(yy0*yy0+(1-xx0)*(1-xx0))*dvertex*c0*0.5);
3165 Double_t z0 = kcl->GetZ();
3166 Double_t zzzz2 = z1-(z1-z3)*dfi1/dfi0;
3167 if (TMath::Abs(zzzz2-z0)>0.5) continue;
3170 Double_t dip = (z1-z0)*c0/dfi1;
3171 Double_t x0 = (0.5*cs13+y0*sn13)*dvertex*c0;
3182 x2= xx2*cs-y2*sn*dsec;
3183 y2=+xx2*sn*dsec+y2*cs;
3193 // do we have cluster at the middle ?
3195 GetProlongation(x1,xm,x,ym,zm);
3197 AliTPCclusterMI * cm=0;
3198 if (TMath::Abs(ym)-ymaxm<0){
3199 cm = krm.FindNearest2(ym,zm,1.0,0.6,dummy);
3200 if ((!cm) || (cm->IsUsed(10))) {
3205 // rotate y1 to system 0
3206 // get state vector in rotated system
3207 Double_t yr1 = (-0.5*sn13+y0*cs13)*dvertex*c0;
3208 Double_t xr2 = x0*cs+yr1*sn*dsec;
3209 Double_t xr[5]={kcl->GetY(),kcl->GetZ(), xr2, dip, c0};
3211 GetProlongation(xx2,xm,xr,ym,zm);
3212 if (TMath::Abs(ym)-ymaxm<0){
3213 cm = kr2m.FindNearest2(ym,zm,1.0,0.6,dummy);
3214 if ((!cm) || (cm->IsUsed(10))) {
3224 dym = ym - cm->GetY();
3225 dzm = zm - cm->GetZ();
3232 Double_t sy1=kr1[is]->GetSigmaY2()*2., sz1=kr1[is]->GetSigmaZ2()*2.;
3233 Double_t sy2=kcl->GetSigmaY2()*2., sz2=kcl->GetSigmaZ2()*2.;
3234 //Double_t sy3=400*3./12., sy=0.1, sz=0.1;
3235 Double_t sy3=25000*x[4]*x[4]+0.1, sy=0.1, sz=0.1;
3236 //Double_t sy3=25000*x[4]*x[4]*60+0.5, sy=0.1, sz=0.1;
3238 Double_t f40=(F1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
3239 Double_t f42=(F1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
3240 Double_t f43=(F1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
3241 Double_t f20=(F2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
3242 Double_t f22=(F2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
3243 Double_t f23=(F2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
3245 Double_t f30=(F3(x1,y1+sy,x2,y2,z1,z2)-x[3])/sy;
3246 Double_t f31=(F3(x1,y1,x2,y2,z1+sz,z2)-x[3])/sz;
3247 Double_t f32=(F3(x1,y1,x2,y2+sy,z1,z2)-x[3])/sy;
3248 Double_t f34=(F3(x1,y1,x2,y2,z1,z2+sz)-x[3])/sz;
3252 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3253 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3254 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3255 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3256 c[13]=f30*sy1*f40+f32*sy2*f42;
3257 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3259 // if (!BuildSeed(kr1[is],kcl,0,x1,x2,x3,x,c)) continue;
3261 UInt_t index=kr1.GetIndex(is);
3262 seed->~AliTPCseed(); // this does not set the pointer to 0...
3263 AliTPCseed *track=new(seed) AliTPCseed(x1, ns*alpha+shift, x, c, index);
3265 track->SetIsSeeding(kTRUE);
3266 track->SetSeed1(i1);
3267 track->SetSeed2(i2);
3268 track->SetSeedType(3);
3272 FollowProlongation(*track, (i1+i2)/2,1);
3273 Int_t foundable,found,shared;
3274 track->GetClusterStatistic((i1+i2)/2,i1, found, foundable, shared, kTRUE);
3275 if ((found<0.55*foundable) || shared>0.5*found || (track->GetSigmaY2()+track->GetSigmaZ2())>0.5){
3277 seed->~AliTPCseed();
3283 FollowProlongation(*track, i2,1);
3287 track->SetBConstrain(1);
3288 // track->fLastPoint = i1+fInnerSec->GetNRows(); // first cluster in track position
3289 track->SetLastPoint(i1); // first cluster in track position
3290 track->SetFirstPoint(track->GetLastPoint());
3292 if (track->GetNumberOfClusters()<(i1-i2)*0.5 ||
3293 track->GetNumberOfClusters() < track->GetNFoundable()*0.6 ||
3294 track->GetNShared()>0.4*track->GetNumberOfClusters() ) {
3296 seed->~AliTPCseed();
3300 // Z VERTEX CONDITION
3301 Double_t zv, bz=GetBz();
3302 if ( !track->GetZAt(0.,bz,zv) ) continue;
3303 if (TMath::Abs(zv-z3)>cuts[2]) {
3304 FollowProlongation(*track, TMath::Max(i2-20,0));
3305 if ( !track->GetZAt(0.,bz,zv) ) continue;
3306 if (TMath::Abs(zv-z3)>cuts[2]){
3307 FollowProlongation(*track, TMath::Max(i2-40,0));
3308 if ( !track->GetZAt(0.,bz,zv) ) continue;
3309 if (TMath::Abs(zv-z3)>cuts[2] &&(track->GetNumberOfClusters() > track->GetNFoundable()*0.7)){
3310 // make seed without constrain
3311 AliTPCseed * track2 = MakeSeed(track,0.2,0.5,1.);
3312 FollowProlongation(*track2, i2,1);
3313 track2->SetBConstrain(kFALSE);
3314 track2->SetSeedType(1);
3315 arr->AddLast(track2);
3317 seed->~AliTPCseed();
3322 seed->~AliTPCseed();
3329 track->SetSeedType(0);
3330 arr->AddLast(track);
3331 seed = new AliTPCseed;
3333 // don't consider other combinations
3334 if (track->GetNumberOfClusters() > track->GetNFoundable()*0.8)
3340 Info("MakeSeeds3","\nSeeding statistic:\t%d\t%d\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin,nout1,nout2);
3346 void AliTPCtrackerMI::MakeSeeds5(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, Float_t cuts[4],
3351 //-----------------------------------------------------------------
3352 // This function creates track seeds.
3353 //-----------------------------------------------------------------
3354 // cuts[0] - fP4 cut
3355 // cuts[1] - tan(phi) cut
3356 // cuts[2] - zvertex cut
3357 // cuts[3] - fP3 cut
3367 Double_t x[5], c[15];
3369 // make temporary seed
3370 AliTPCseed * seed = new AliTPCseed;
3371 Double_t alpha=fOuterSec->GetAlpha(), shift=fOuterSec->GetAlphaShift();
3372 // Double_t cs=cos(alpha), sn=sin(alpha);
3377 Double_t x1 = GetXrow(i1-1);
3378 const AliTPCtrackerRow& kr1=GetRow(sec,i1-1);
3379 Double_t y1max = GetMaxY(i1-1)-kr1.GetDeadZone()-1.5;
3381 Double_t x1p = GetXrow(i1);
3382 const AliTPCtrackerRow& kr1p=GetRow(sec,i1);
3384 Double_t x1m = GetXrow(i1-2);
3385 const AliTPCtrackerRow& kr1m=GetRow(sec,i1-2);
3388 //last 3 padrow for seeding
3389 AliTPCtrackerRow& kr3 = GetRow((sec+fkNOS)%fkNOS,i1-7);
3390 Double_t x3 = GetXrow(i1-7);
3391 // Double_t y3max= GetMaxY(i1-7)-kr3.fDeadZone-1.5;
3393 AliTPCtrackerRow& kr3p = GetRow((sec+fkNOS)%fkNOS,i1-6);
3394 Double_t x3p = GetXrow(i1-6);
3396 AliTPCtrackerRow& kr3m = GetRow((sec+fkNOS)%fkNOS,i1-8);
3397 Double_t x3m = GetXrow(i1-8);
3402 Int_t im = i1-4; //middle pad row index
3403 Double_t xm = GetXrow(im); // radius of middle pad-row
3404 const AliTPCtrackerRow& krm=GetRow(sec,im); //middle pad -row
3405 // Double_t ymmax = GetMaxY(im)-kr1.fDeadZone-1.5;
3408 Double_t deltax = x1-x3;
3409 Double_t dymax = deltax*cuts[1];
3410 Double_t dzmax = deltax*cuts[3];
3412 // loop over clusters
3413 for (Int_t is=0; is < kr1; is++) {
3415 if (kr1[is]->IsUsed(10)) continue;
3416 Double_t y1=kr1[is]->GetY(), z1=kr1[is]->GetZ();
3418 if (deltay>0 && TMath::Abs(y1max-TMath::Abs(y1))> deltay ) continue; // seed only at the edge
3420 Int_t index1 = TMath::Max(kr3.Find(z1-dzmax)-1,0);
3421 Int_t index2 = TMath::Min(kr3.Find(z1+dzmax)+1,kr3);
3427 for (Int_t js=index1; js < index2; js++) {
3428 const AliTPCclusterMI *kcl = kr3[js];
3429 if (kcl->IsUsed(10)) continue;
3431 // apply angular cuts
3432 if (TMath::Abs(y1-y3)>dymax) continue;
3435 if (TMath::Abs(z1-z3)>dzmax) continue;
3437 Double_t angley = (y1-y3)/(x1-x3);
3438 Double_t anglez = (z1-z3)/(x1-x3);
3440 Double_t erry = TMath::Abs(angley)*(x1-x1m)*0.5+0.5;
3441 Double_t errz = TMath::Abs(anglez)*(x1-x1m)*0.5+0.5;
3443 Double_t yyym = angley*(xm-x1)+y1;
3444 Double_t zzzm = anglez*(xm-x1)+z1;
3446 const AliTPCclusterMI *kcm = krm.FindNearest2(yyym,zzzm,erry,errz,index);
3448 if (kcm->IsUsed(10)) continue;
3450 erry = TMath::Abs(angley)*(x1-x1m)*0.4+0.5;
3451 errz = TMath::Abs(anglez)*(x1-x1m)*0.4+0.5;
3458 // look around first
3459 const AliTPCclusterMI *kc1m = kr1m.FindNearest2(angley*(x1m-x1)+y1,
3465 if (kc1m->IsUsed(10)) used++;
3467 const AliTPCclusterMI *kc1p = kr1p.FindNearest2(angley*(x1p-x1)+y1,
3473 if (kc1p->IsUsed(10)) used++;
3475 if (used>1) continue;
3476 if (found<1) continue;
3480 const AliTPCclusterMI *kc3m = kr3m.FindNearest2(angley*(x3m-x3)+y3,
3486 if (kc3m->IsUsed(10)) used++;
3490 const AliTPCclusterMI *kc3p = kr3p.FindNearest2(angley*(x3p-x3)+y3,
3496 if (kc3p->IsUsed(10)) used++;
3500 if (used>1) continue;
3501 if (found<3) continue;
3511 x[4]=F1(x1,y1,x2,y2,x3,y3);
3512 //if (TMath::Abs(x[4]) >= cuts[0]) continue;
3515 x[2]=F2(x1,y1,x2,y2,x3,y3);
3518 x[3]=F3n(x1,y1,x2,y2,z1,z2,x[4]);
3519 //if (TMath::Abs(x[3]) > cuts[3]) continue;
3523 Double_t sy1=0.1, sz1=0.1;
3524 Double_t sy2=0.1, sz2=0.1;
3525 Double_t sy3=0.1, sy=0.1, sz=0.1;
3527 Double_t f40=(F1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
3528 Double_t f42=(F1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
3529 Double_t f43=(F1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
3530 Double_t f20=(F2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
3531 Double_t f22=(F2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
3532 Double_t f23=(F2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
3534 Double_t f30=(F3(x1,y1+sy,x2,y2,z1,z2)-x[3])/sy;
3535 Double_t f31=(F3(x1,y1,x2,y2,z1+sz,z2)-x[3])/sz;
3536 Double_t f32=(F3(x1,y1,x2,y2+sy,z1,z2)-x[3])/sy;
3537 Double_t f34=(F3(x1,y1,x2,y2,z1,z2+sz)-x[3])/sz;
3541 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3542 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3543 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3544 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3545 c[13]=f30*sy1*f40+f32*sy2*f42;
3546 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3548 // if (!BuildSeed(kr1[is],kcl,0,x1,x2,x3,x,c)) continue;
3550 index=kr1.GetIndex(is);
3551 seed->~AliTPCseed();
3552 AliTPCseed *track=new(seed) AliTPCseed(x1, sec*alpha+shift, x, c, index);
3554 track->SetIsSeeding(kTRUE);
3557 FollowProlongation(*track, i1-7,1);
3558 if (track->GetNumberOfClusters() < track->GetNFoundable()*0.75 ||
3559 track->GetNShared()>0.6*track->GetNumberOfClusters() || ( track->GetSigmaY2()+ track->GetSigmaZ2())>0.6){
3561 seed->~AliTPCseed();
3567 FollowProlongation(*track, i2,1);
3568 track->SetBConstrain(0);
3569 track->SetLastPoint(i1+fInnerSec->GetNRows()); // first cluster in track position
3570 track->SetFirstPoint(track->GetLastPoint());
3572 if (track->GetNumberOfClusters()<(i1-i2)*0.5 ||
3573 track->GetNumberOfClusters()<track->GetNFoundable()*0.7 ||
3574 track->GetNShared()>2. || track->GetChi2()/track->GetNumberOfClusters()>6 || ( track->GetSigmaY2()+ track->GetSigmaZ2())>0.5 ) {
3576 seed->~AliTPCseed();
3581 FollowProlongation(*track, TMath::Max(i2-10,0),1);
3582 AliTPCseed * track2 = MakeSeed(track,0.2,0.5,0.9);
3583 FollowProlongation(*track2, i2,1);
3584 track2->SetBConstrain(kFALSE);
3585 track2->SetSeedType(4);
3586 arr->AddLast(track2);
3588 seed->~AliTPCseed();
3592 //arr->AddLast(track);
3593 //seed = new AliTPCseed;
3599 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);
3605 //_____________________________________________________________________________
3606 void AliTPCtrackerMI::MakeSeeds2(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, Float_t */*cuts[4]*/,
3607 Float_t deltay, Bool_t /*bconstrain*/) {
3608 //-----------------------------------------------------------------
3609 // This function creates track seeds - without vertex constraint
3610 //-----------------------------------------------------------------
3611 // cuts[0] - fP4 cut - not applied
3612 // cuts[1] - tan(phi) cut
3613 // cuts[2] - zvertex cut - not applied
3614 // cuts[3] - fP3 cut
3624 Double_t alpha=fOuterSec->GetAlpha(), shift=fOuterSec->GetAlphaShift();
3625 // Double_t cs=cos(alpha), sn=sin(alpha);
3626 Int_t row0 = (i1+i2)/2;
3627 Int_t drow = (i1-i2)/2;
3628 const AliTPCtrackerRow& kr0=fSectors[sec][row0];
3629 AliTPCtrackerRow * kr=0;
3631 AliTPCpolyTrack polytrack;
3632 Int_t nclusters=fSectors[sec][row0];
3633 AliTPCseed * seed = new AliTPCseed;
3638 for (Int_t is=0; is < nclusters; is++) { //LOOP over clusters
3640 Int_t nfoundable =0;
3641 for (Int_t iter =1; iter<2; iter++){ //iterations
3642 const AliTPCtrackerRow& krm=fSectors[sec][row0-iter];
3643 const AliTPCtrackerRow& krp=fSectors[sec][row0+iter];
3644 const AliTPCclusterMI * cl= kr0[is];
3646 if (cl->IsUsed(10)) {
3652 Double_t x = kr0.GetX();
3653 // Initialization of the polytrack
3658 Double_t y0= cl->GetY();
3659 Double_t z0= cl->GetZ();
3663 Double_t ymax = fSectors->GetMaxY(row0)-kr0.GetDeadZone()-1.5;
3664 if (deltay>0 && TMath::Abs(ymax-TMath::Abs(y0))> deltay ) continue; // seed only at the edge
3666 erry = (0.5)*cl->GetSigmaY2()/TMath::Sqrt(cl->GetQ())*6;
3667 errz = (0.5)*cl->GetSigmaZ2()/TMath::Sqrt(cl->GetQ())*6;
3668 polytrack.AddPoint(x,y0,z0,erry, errz);
3671 if (cl->IsUsed(10)) sumused++;
3674 Float_t roady = (5*TMath::Sqrt(cl->GetSigmaY2()+0.2)+1.)*iter;
3675 Float_t roadz = (5*TMath::Sqrt(cl->GetSigmaZ2()+0.2)+1.)*iter;
3678 AliTPCclusterMI * cl1 = krm.FindNearest(y0,z0,roady,roadz);
3679 if (cl1 && TMath::Abs(ymax-TMath::Abs(y0))) {
3680 erry = (0.5)*cl1->GetSigmaY2()/TMath::Sqrt(cl1->GetQ())*3;
3681 errz = (0.5)*cl1->GetSigmaZ2()/TMath::Sqrt(cl1->GetQ())*3;
3682 if (cl1->IsUsed(10)) sumused++;
3683 polytrack.AddPoint(x,cl1->GetY(),cl1->GetZ(),erry,errz);
3687 AliTPCclusterMI * cl2 = krp.FindNearest(y0,z0,roady,roadz);
3689 erry = (0.5)*cl2->GetSigmaY2()/TMath::Sqrt(cl2->GetQ())*3;
3690 errz = (0.5)*cl2->GetSigmaZ2()/TMath::Sqrt(cl2->GetQ())*3;
3691 if (cl2->IsUsed(10)) sumused++;
3692 polytrack.AddPoint(x,cl2->GetY(),cl2->GetZ(),erry,errz);
3695 if (sumused>0) continue;
3697 polytrack.UpdateParameters();
3703 nfoundable = polytrack.GetN();
3704 nfound = nfoundable;
3706 for (Int_t ddrow = iter+1; ddrow<drow;ddrow++){
3707 Float_t maxdist = 0.8*(1.+3./(ddrow));
3708 for (Int_t delta = -1;delta<=1;delta+=2){
3709 Int_t row = row0+ddrow*delta;
3710 kr = &(fSectors[sec][row]);
3711 Double_t xn = kr->GetX();
3712 Double_t ymax1 = fSectors->GetMaxY(row)-kr->GetDeadZone()-1.5;
3713 polytrack.GetFitPoint(xn,yn,zn);
3714 if (TMath::Abs(yn)>ymax1) continue;
3716 AliTPCclusterMI * cln = kr->FindNearest(yn,zn,roady,roadz);
3718 Float_t dist = TMath::Sqrt( (yn-cln->GetY())*(yn-cln->GetY())+(zn-cln->GetZ())*(zn-cln->GetZ()));
3721 erry = (dist+0.3)*cln->GetSigmaY2()/TMath::Sqrt(cln->GetQ())*(1.+1./(ddrow));
3722 errz = (dist+0.3)*cln->GetSigmaZ2()/TMath::Sqrt(cln->GetQ())*(1.+1./(ddrow));
3723 if (cln->IsUsed(10)) {
3724 // printf("used\n");
3732 polytrack.AddPoint(xn,cln->GetY(),cln->GetZ(),erry, errz);
3737 if ( (sumused>3) || (sumused>0.5*nfound) || (nfound<0.6*nfoundable)) break;
3738 polytrack.UpdateParameters();
3741 if ( (sumused>3) || (sumused>0.5*nfound)) {
3742 //printf("sumused %d\n",sumused);
3747 polytrack.GetFitDerivation(kr0.GetX(),dy,dz);
3748 AliTPCpolyTrack track2;
3750 polytrack.Refit(track2,0.5+TMath::Abs(dy)*0.3,0.4+TMath::Abs(dz)*0.3);
3751 if (track2.GetN()<0.5*nfoundable) continue;
3754 if ((nfound>0.6*nfoundable) &&( nfoundable>0.4*(i1-i2))) {
3756 // test seed with and without constrain
3757 for (Int_t constrain=0; constrain<=0;constrain++){
3758 // add polytrack candidate
3760 Double_t x[5], c[15];
3761 Double_t x1,x2,x3,y1,y2,y3,z1,z2,z3;
3762 track2.GetBoundaries(x3,x1);
3764 track2.GetFitPoint(x1,y1,z1);
3765 track2.GetFitPoint(x2,y2,z2);
3766 track2.GetFitPoint(x3,y3,z3);
3768 //is track pointing to the vertex ?
3771 polytrack.GetFitPoint(x0,y0,z0);
3784 x[4]=F1(x1,y1,x2,y2,x3,y3);
3786 // if (TMath::Abs(x[4]) >= cuts[0]) continue; //
3787 x[2]=F2(x1,y1,x2,y2,x3,y3);
3789 //if (TMath::Abs(x[4]*x1-x[2]) >= cuts[1]) continue;
3790 //x[3]=F3(x1,y1,x2,y2,z1,z2);
3791 x[3]=F3n(x1,y1,x3,y3,z1,z3,x[4]);
3792 //if (TMath::Abs(x[3]) > cuts[3]) continue;
3795 Double_t sy =0.1, sz =0.1;
3796 Double_t sy1=0.02, sz1=0.02;
3797 Double_t sy2=0.02, sz2=0.02;
3801 sy3=25000*x[4]*x[4]+0.1, sy=0.1, sz=0.1;
3804 Double_t f40=(F1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
3805 Double_t f42=(F1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
3806 Double_t f43=(F1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
3807 Double_t f20=(F2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
3808 Double_t f22=(F2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
3809 Double_t f23=(F2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
3811 Double_t f30=(F3(x1,y1+sy,x3,y3,z1,z3)-x[3])/sy;
3812 Double_t f31=(F3(x1,y1,x3,y3,z1+sz,z3)-x[3])/sz;
3813 Double_t f32=(F3(x1,y1,x3,y3+sy,z1,z3)-x[3])/sy;
3814 Double_t f34=(F3(x1,y1,x3,y3,z1,z3+sz)-x[3])/sz;
3819 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3820 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3821 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3822 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3823 c[13]=f30*sy1*f40+f32*sy2*f42;
3824 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3826 //Int_t row1 = fSectors->GetRowNumber(x1);
3827 Int_t row1 = GetRowNumber(x1);
3831 seed->~AliTPCseed();
3832 AliTPCseed *track=new(seed) AliTPCseed(x1,sec*alpha+shift,x,c,index);
3833 track->SetIsSeeding(kTRUE);
3834 Int_t rc=FollowProlongation(*track, i2);
3835 if (constrain) track->SetBConstrain(1);
3837 track->SetBConstrain(0);
3838 track->SetLastPoint(row1+fInnerSec->GetNRows()); // first cluster in track position
3839 track->SetFirstPoint(track->GetLastPoint());
3841 if (rc==0 || track->GetNumberOfClusters()<(i1-i2)*0.5 ||
3842 track->GetNumberOfClusters() < track->GetNFoundable()*0.6 ||
3843 track->GetNShared()>0.4*track->GetNumberOfClusters()) {
3846 seed->~AliTPCseed();
3849 arr->AddLast(track);
3850 seed = new AliTPCseed;
3854 } // if accepted seed
3857 Info("MakeSeeds2","\nSeeding statiistic:\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin3);
3863 AliTPCseed *AliTPCtrackerMI::MakeSeed(AliTPCseed *const track, Float_t r0, Float_t r1, Float_t r2)
3867 //reseed using track points
3868 Int_t p0 = int(r0*track->GetNumberOfClusters()); // point 0
3869 Int_t p1 = int(r1*track->GetNumberOfClusters());
3870 Int_t p2 = int(r2*track->GetNumberOfClusters()); // last point
3872 Double_t x0[3],x1[3],x2[3];
3873 for (Int_t i=0;i<3;i++){
3879 // find track position at given ratio of the length
3880 Int_t sec0=0, sec1=0, sec2=0;
3883 for (Int_t i=0;i<160;i++){
3884 if (track->GetClusterPointer(i)){
3886 AliTPCTrackerPoint *trpoint =track->GetTrackPoint(i);
3887 if ( (index<p0) || x0[0]<0 ){
3888 if (trpoint->GetX()>1){
3889 clindex = track->GetClusterIndex2(i);
3891 x0[0] = trpoint->GetX();
3892 x0[1] = trpoint->GetY();
3893 x0[2] = trpoint->GetZ();
3894 sec0 = ((clindex&0xff000000)>>24)%18;
3899 if ( (index<p1) &&(trpoint->GetX()>1)){
3900 clindex = track->GetClusterIndex2(i);
3902 x1[0] = trpoint->GetX();
3903 x1[1] = trpoint->GetY();
3904 x1[2] = trpoint->GetZ();
3905 sec1 = ((clindex&0xff000000)>>24)%18;
3908 if ( (index<p2) &&(trpoint->GetX()>1)){
3909 clindex = track->GetClusterIndex2(i);
3911 x2[0] = trpoint->GetX();
3912 x2[1] = trpoint->GetY();
3913 x2[2] = trpoint->GetZ();
3914 sec2 = ((clindex&0xff000000)>>24)%18;
3921 Double_t alpha, cs,sn, xx2,yy2;
3923 alpha = (sec1-sec2)*fSectors->GetAlpha();
3924 cs = TMath::Cos(alpha);
3925 sn = TMath::Sin(alpha);
3926 xx2= x1[0]*cs-x1[1]*sn;
3927 yy2= x1[0]*sn+x1[1]*cs;
3931 alpha = (sec0-sec2)*fSectors->GetAlpha();
3932 cs = TMath::Cos(alpha);
3933 sn = TMath::Sin(alpha);
3934 xx2= x0[0]*cs-x0[1]*sn;
3935 yy2= x0[0]*sn+x0[1]*cs;
3941 Double_t x[5],c[15];
3945 x[4]=F1(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]);
3946 // if (x[4]>1) return 0;
3947 x[2]=F2(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]);
3948 x[3]=F3n(x2[0],x2[1],x0[0],x0[1],x2[2],x0[2],x[4]);
3949 //if (TMath::Abs(x[3]) > 2.2) return 0;
3950 //if (TMath::Abs(x[2]) > 1.99) return 0;
3952 Double_t sy =0.1, sz =0.1;
3954 Double_t sy1=0.02+track->GetSigmaY2(), sz1=0.02+track->GetSigmaZ2();
3955 Double_t sy2=0.01+track->GetSigmaY2(), sz2=0.01+track->GetSigmaZ2();
3956 Double_t sy3=0.01+track->GetSigmaY2();
3958 Double_t f40=(F1(x2[0],x2[1]+sy,x1[0],x1[1],x0[0],x0[1])-x[4])/sy;
3959 Double_t f42=(F1(x2[0],x2[1],x1[0],x1[1]+sy,x0[0],x0[1])-x[4])/sy;
3960 Double_t f43=(F1(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]+sy)-x[4])/sy;
3961 Double_t f20=(F2(x2[0],x2[1]+sy,x1[0],x1[1],x0[0],x0[1])-x[2])/sy;
3962 Double_t f22=(F2(x2[0],x2[1],x1[0],x1[1]+sy,x0[0],x0[1])-x[2])/sy;
3963 Double_t f23=(F2(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]+sy)-x[2])/sy;
3965 Double_t f30=(F3(x2[0],x2[1]+sy,x0[0],x0[1],x2[2],x0[2])-x[3])/sy;
3966 Double_t f31=(F3(x2[0],x2[1],x0[0],x0[1],x2[2]+sz,x0[2])-x[3])/sz;
3967 Double_t f32=(F3(x2[0],x2[1],x0[0],x0[1]+sy,x2[2],x0[2])-x[3])/sy;
3968 Double_t f34=(F3(x2[0],x2[1],x0[0],x0[1],x2[2],x0[2]+sz)-x[3])/sz;
3973 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3974 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3975 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3976 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3977 c[13]=f30*sy1*f40+f32*sy2*f42;
3978 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3980 // Int_t row1 = fSectors->GetRowNumber(x2[0]);
3981 AliTPCseed *seed=new AliTPCseed(x2[0], sec2*fSectors->GetAlpha()+fSectors->GetAlphaShift(), x, c, 0);
3982 // Double_t y0,z0,y1,z1, y2,z2;
3983 //seed->GetProlongation(x0[0],y0,z0);
3984 // seed->GetProlongation(x1[0],y1,z1);
3985 //seed->GetProlongation(x2[0],y2,z2);
3987 seed->SetLastPoint(pp2);
3988 seed->SetFirstPoint(pp2);
3995 AliTPCseed *AliTPCtrackerMI::ReSeed(const AliTPCseed *track, Float_t r0, Float_t r1, Float_t r2)
3999 //reseed using founded clusters
4001 // Find the number of clusters
4002 Int_t nclusters = 0;
4003 for (Int_t irow=0;irow<160;irow++){
4004 if (track->GetClusterIndex(irow)>0) nclusters++;
4008 ipos[0] = TMath::Max(int(r0*nclusters),0); // point 0 cluster
4009 ipos[1] = TMath::Min(int(r1*nclusters),nclusters-1); //
4010 ipos[2] = TMath::Min(int(r2*nclusters),nclusters-1); // last point
4013 Double_t xyz[3][3]={{0}};
4014 Int_t row[3]={0},sec[3]={0,0,0};
4016 // find track row position at given ratio of the length
4018 for (Int_t irow=0;irow<160;irow++){
4019 if (track->GetClusterIndex2(irow)<0) continue;
4021 for (Int_t ipoint=0;ipoint<3;ipoint++){
4022 if (index<=ipos[ipoint]) row[ipoint] = irow;
4026 //Get cluster and sector position
4027 for (Int_t ipoint=0;ipoint<3;ipoint++){
4028 Int_t clindex = track->GetClusterIndex2(row[ipoint]);
4029 AliTPCclusterMI * cl = GetClusterMI(clindex);
4032 // AliTPCclusterMI * cl = GetClusterMI(clindex);
4035 sec[ipoint] = ((clindex&0xff000000)>>24)%18;
4036 xyz[ipoint][0] = GetXrow(row[ipoint]);
4037 xyz[ipoint][1] = cl->GetY();
4038 xyz[ipoint][2] = cl->GetZ();
4042 // Calculate seed state vector and covariance matrix
4044 Double_t alpha, cs,sn, xx2,yy2;
4046 alpha = (sec[1]-sec[2])*fSectors->GetAlpha();
4047 cs = TMath::Cos(alpha);
4048 sn = TMath::Sin(alpha);
4049 xx2= xyz[1][0]*cs-xyz[1][1]*sn;
4050 yy2= xyz[1][0]*sn+xyz[1][1]*cs;
4054 alpha = (sec[0]-sec[2])*fSectors->GetAlpha();
4055 cs = TMath::Cos(alpha);
4056 sn = TMath::Sin(alpha);
4057 xx2= xyz[0][0]*cs-xyz[0][1]*sn;
4058 yy2= xyz[0][0]*sn+xyz[0][1]*cs;
4064 Double_t x[5],c[15];
4068 x[4]=F1(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
4069 x[2]=F2(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
4070 x[3]=F3n(xyz[2][0],xyz[2][1],xyz[0][0],xyz[0][1],xyz[2][2],xyz[0][2],x[4]);
4072 Double_t sy =0.1, sz =0.1;
4074 Double_t sy1=0.2, sz1=0.2;
4075 Double_t sy2=0.2, sz2=0.2;
4078 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;
4079 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;
4080 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;
4081 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;
4082 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;
4083 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;
4085 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;
4086 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;
4087 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;
4088 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;
4093 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
4094 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
4095 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
4096 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
4097 c[13]=f30*sy1*f40+f32*sy2*f42;
4098 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
4100 // Int_t row1 = fSectors->GetRowNumber(xyz[2][0]);
4101 AliTPCseed *seed=new AliTPCseed(xyz[2][0], sec[2]*fSectors->GetAlpha()+fSectors->GetAlphaShift(), x, c, 0);
4102 seed->SetLastPoint(row[2]);
4103 seed->SetFirstPoint(row[2]);
4108 AliTPCseed *AliTPCtrackerMI::ReSeed(AliTPCseed *track,Int_t r0, Bool_t forward)
4112 //reseed using founded clusters
4115 Int_t row[3]={0,0,0};
4116 Int_t sec[3]={0,0,0};
4118 // forward direction
4120 for (Int_t irow=r0;irow<160;irow++){
4121 if (track->GetClusterIndex(irow)>0){
4126 for (Int_t irow=160;irow>r0;irow--){
4127 if (track->GetClusterIndex(irow)>0){
4132 for (Int_t irow=row[2]-15;irow>row[0];irow--){
4133 if (track->GetClusterIndex(irow)>0){
4141 for (Int_t irow=0;irow<r0;irow++){
4142 if (track->GetClusterIndex(irow)>0){
4147 for (Int_t irow=r0;irow>0;irow--){
4148 if (track->GetClusterIndex(irow)>0){
4153 for (Int_t irow=row[2]-15;irow>row[0];irow--){
4154 if (track->GetClusterIndex(irow)>0){
4161 if ((row[2]-row[0])<20) return 0;
4162 if (row[1]==0) return 0;
4165 //Get cluster and sector position
4166 for (Int_t ipoint=0;ipoint<3;ipoint++){
4167 Int_t clindex = track->GetClusterIndex2(row[ipoint]);
4168 AliTPCclusterMI * cl = GetClusterMI(clindex);
4171 // AliTPCclusterMI * cl = GetClusterMI(clindex);
4174 sec[ipoint] = ((clindex&0xff000000)>>24)%18;
4175 xyz[ipoint][0] = GetXrow(row[ipoint]);
4176 AliTPCTrackerPoint * point = track->GetTrackPoint(row[ipoint]);
4177 if (point&&ipoint<2){
4179 xyz[ipoint][1] = point->GetY();
4180 xyz[ipoint][2] = point->GetZ();
4183 xyz[ipoint][1] = cl->GetY();
4184 xyz[ipoint][2] = cl->GetZ();
4191 // Calculate seed state vector and covariance matrix
4193 Double_t alpha, cs,sn, xx2,yy2;
4195 alpha = (sec[1]-sec[2])*fSectors->GetAlpha();
4196 cs = TMath::Cos(alpha);
4197 sn = TMath::Sin(alpha);
4198 xx2= xyz[1][0]*cs-xyz[1][1]*sn;
4199 yy2= xyz[1][0]*sn+xyz[1][1]*cs;
4203 alpha = (sec[0]-sec[2])*fSectors->GetAlpha();
4204 cs = TMath::Cos(alpha);
4205 sn = TMath::Sin(alpha);
4206 xx2= xyz[0][0]*cs-xyz[0][1]*sn;
4207 yy2= xyz[0][0]*sn+xyz[0][1]*cs;
4213 Double_t x[5],c[15];
4217 x[4]=F1(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
4218 x[2]=F2(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
4219 x[3]=F3n(xyz[2][0],xyz[2][1],xyz[0][0],xyz[0][1],xyz[2][2],xyz[0][2],x[4]);
4221 Double_t sy =0.1, sz =0.1;
4223 Double_t sy1=0.2, sz1=0.2;
4224 Double_t sy2=0.2, sz2=0.2;
4227 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;
4228 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;
4229 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;
4230 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;
4231 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;
4232 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;
4234 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;
4235 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;
4236 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;
4237 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;
4242 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
4243 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
4244 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
4245 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
4246 c[13]=f30*sy1*f40+f32*sy2*f42;
4247 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
4249 // Int_t row1 = fSectors->GetRowNumber(xyz[2][0]);
4250 AliTPCseed *seed=new AliTPCseed(xyz[2][0], sec[2]*fSectors->GetAlpha()+fSectors->GetAlphaShift(), x, c, 0);
4251 seed->SetLastPoint(row[2]);
4252 seed->SetFirstPoint(row[2]);
4253 for (Int_t i=row[0];i<row[2];i++){
4254 seed->SetClusterIndex(i, track->GetClusterIndex(i));
4262 void AliTPCtrackerMI::FindMultiMC(const TObjArray * array, AliESDEvent */*esd*/, Int_t iter)
4265 // find multi tracks - THIS FUNCTION IS ONLY FOR DEBUG PURPOSES
4267 // Use AliTPCReconstructor::StreamLevel()>2 if you want to tune parameters - cuts
4269 // Two reasons to have multiple find tracks
4270 // 1. Curling tracks can be find more than once
4271 // 2. Splitted tracks
4272 // a.) Multiple seeding to increase tracking efficiency - (~ 100% reached)
4273 // b.) Edge effect on the sector boundaries
4276 // Algorithm done in 2 phases - because of CPU consumption
4277 // it is n^2 algorithm - for lead-lead 20000x20000 combination are investigated
4279 // Algorihm for curling tracks sign:
4280 // 1 phase -makes a very rough fast cuts to minimize combinatorics
4281 // a.) opposite sign
4282 // b.) one of the tracks - not pointing to the primary vertex -
4283 // c.) delta tan(theta)
4285 // 2 phase - calculates DCA between tracks - time consument
4290 // General cuts - for splitted tracks and for curling tracks
4292 const Float_t kMaxdPhi = 0.2; // maximal distance in phi
4294 // Curling tracks cuts
4299 Int_t nentries = array->GetEntriesFast();
4300 AliHelix *helixes = new AliHelix[nentries];
4301 Float_t *xm = new Float_t[nentries];
4302 Float_t *dz0 = new Float_t[nentries];
4303 Float_t *dz1 = new Float_t[nentries];
4309 // Find track COG in x direction - point with best defined parameters
4311 for (Int_t i=0;i<nentries;i++){
4312 AliTPCseed* track = (AliTPCseed*)array->At(i);
4313 if (!track) continue;
4314 track->SetCircular(0);
4315 new (&helixes[i]) AliHelix(*track);
4319 track->GetDZ(GetX(),GetY(),GetZ(),GetBz(),dz);
4322 for (Int_t icl=0; icl<160; icl++){
4323 AliTPCclusterMI * cl = track->GetClusterPointer(icl);
4329 if (ncl>0) xm[i]/=Float_t(ncl);
4332 for (Int_t i0=0;i0<nentries;i0++){
4333 AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
4334 if (!track0) continue;
4335 Float_t xc0 = helixes[i0].GetHelix(6);
4336 Float_t yc0 = helixes[i0].GetHelix(7);
4337 Float_t r0 = helixes[i0].GetHelix(8);
4338 Float_t rc0 = TMath::Sqrt(xc0*xc0+yc0*yc0);
4339 Float_t fi0 = TMath::ATan2(yc0,xc0);
4341 for (Int_t i1=i0+1;i1<nentries;i1++){
4342 AliTPCseed * track1 = (AliTPCseed*)array->At(i1);
4343 if (!track1) continue;
4344 Int_t lab0=track0->GetLabel();
4345 Int_t lab1=track1->GetLabel();
4346 if (TMath::Abs(lab0)!=TMath::Abs(lab1)) continue;
4348 Float_t xc1 = helixes[i1].GetHelix(6);
4349 Float_t yc1 = helixes[i1].GetHelix(7);
4350 Float_t r1 = helixes[i1].GetHelix(8);
4351 Float_t rc1 = TMath::Sqrt(xc1*xc1+yc1*yc1);
4352 Float_t fi1 = TMath::ATan2(yc1,xc1);
4354 Float_t dfi = fi0-fi1;
4357 if (dfi>1.5*TMath::Pi()) dfi-=TMath::Pi(); // take care about edge effect
4358 if (dfi<-1.5*TMath::Pi()) dfi+=TMath::Pi(); //
4359 if (TMath::Abs(dfi)>kMaxdPhi&&helixes[i0].GetHelix(4)*helixes[i1].GetHelix(4)<0){
4361 // if short tracks with undefined sign
4362 fi1 = -TMath::ATan2(yc1,-xc1);
4365 Float_t dtheta = TMath::Abs(track0->GetTgl()-track1->GetTgl())<TMath::Abs(track0->GetTgl()+track1->GetTgl())? track0->GetTgl()-track1->GetTgl():track0->GetTgl()+track1->GetTgl();
4368 // debug stream to tune "fast cuts"
4370 Double_t dist[3]; // distance at X
4371 Double_t mdist[3]={0,0,0}; // mean distance X+-40cm
4372 track0->GetDistance(track1,0.5*(xm[i0]+xm[i1])-40.,dist,AliTracker::GetBz());
4373 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[i]);
4374 track0->GetDistance(track1,0.5*(xm[i0]+xm[i1])+40.,dist,AliTracker::GetBz());
4375 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[i]);
4376 track0->GetDistance(track1,0.5*(xm[i0]+xm[i1]),dist,AliTracker::GetBz());
4377 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[i]);
4378 for (Int_t i=0;i<3;i++) mdist[i]*=0.33333;
4382 for (Int_t icl=0; icl<160; icl++){
4383 AliTPCclusterMI * cl0 = track0->GetClusterPointer(icl);
4384 AliTPCclusterMI * cl1 = track1->GetClusterPointer(icl);
4387 if (cl0==cl1) sums++;
4391 if (AliTPCReconstructor::StreamLevel()>5) {
4392 TTreeSRedirector &cstream = *fDebugStreamer;
4397 "Tr0.="<<track0<< // seed0
4398 "Tr1.="<<track1<< // seed1
4399 "h0.="<<&helixes[i0]<<
4400 "h1.="<<&helixes[i1]<<
4402 "sum="<<sum<< //the sum of rows with cl in both
4403 "sums="<<sums<< //the sum of shared clusters
4404 "xm0="<<xm[i0]<< // the center of track
4405 "xm1="<<xm[i1]<< // the x center of track
4406 // General cut variables
4407 "dfi="<<dfi<< // distance in fi angle
4408 "dtheta="<<dtheta<< // distance int theta angle
4414 "dist0="<<dist[0]<< //distance x
4415 "dist1="<<dist[1]<< //distance y
4416 "dist2="<<dist[2]<< //distance z
4417 "mdist0="<<mdist[0]<< //distance x
4418 "mdist1="<<mdist[1]<< //distance y
4419 "mdist2="<<mdist[2]<< //distance z
4435 if (AliTPCReconstructor::StreamLevel()>1) {
4436 AliInfo("Time for curling tracks removal DEBUGGING MC");
4443 void AliTPCtrackerMI::FindSplitted(TObjArray * array, AliESDEvent */*esd*/, Int_t /*iter*/){
4445 // Find Splitted tracks and remove the one with worst quality
4446 // Corresponding debug streamer to tune selections - "Splitted2"
4448 // 0. Sort tracks according quility
4449 // 1. Propagate the tracks to the reference radius
4450 // 2. Double_t loop to select close tracks (only to speed up process)
4451 // 3. Calculate cluster overlap ratio - and remove the track if bigger than a threshold
4452 // 4. Delete temporary parameters
4454 const Double_t xref=GetXrow(63); // reference radius -IROC/OROC boundary
4456 const Double_t kCutP1=10; // delta Z cut 10 cm
4457 const Double_t kCutP2=0.15; // delta snp(fi) cut 0.15
4458 const Double_t kCutP3=0.15; // delta tgl(theta) cut 0.15
4459 const Double_t kCutAlpha=0.15; // delta alpha cut
4460 Int_t firstpoint = 0;
4461 Int_t lastpoint = 160;
4463 Int_t nentries = array->GetEntriesFast();
4464 AliExternalTrackParam *params = new AliExternalTrackParam[nentries];
4470 //0. Sort tracks according quality
4471 //1. Propagate the ext. param to reference radius
4472 Int_t nseed = array->GetEntriesFast();
4473 if (nseed<=0) return;
4474 Float_t * quality = new Float_t[nseed];
4475 Int_t * indexes = new Int_t[nseed];
4476 for (Int_t i=0; i<nseed; i++) {
4477 AliTPCseed *pt=(AliTPCseed*)array->UncheckedAt(i);
4482 pt->UpdatePoints(); //select first last max dens points
4483 Float_t * points = pt->GetPoints();
4484 if (points[3]<0.8) quality[i] =-1;
4485 quality[i] = (points[2]-points[0])+pt->GetNumberOfClusters();
4486 //prefer high momenta tracks if overlaps
4487 quality[i] *= TMath::Sqrt(TMath::Abs(pt->Pt())+0.5);
4489 AliTracker::PropagateTrackToBxByBz(&(params[i]),xref,pt->GetMass(),5.,kTRUE);
4490 AliTracker::PropagateTrackToBxByBz(&(params[i]),xref,pt->GetMass(),1.,kTRUE);
4492 TMath::Sort(nseed,quality,indexes);
4494 // 3. Loop over pair of tracks
4496 for (Int_t i0=0; i0<nseed; i0++) {
4497 Int_t index0=indexes[i0];
4498 if (!(array->UncheckedAt(index0))) continue;
4499 AliTPCseed *s1 = (AliTPCseed*)array->UncheckedAt(index0);
4500 if (!s1->IsActive()) continue;
4501 AliExternalTrackParam &par0=params[index0];
4502 for (Int_t i1=i0+1; i1<nseed; i1++) {
4503 Int_t index1=indexes[i1];
4504 if (!(array->UncheckedAt(index1))) continue;
4505 AliTPCseed *s2 = (AliTPCseed*)array->UncheckedAt(index1);
4506 if (!s2->IsActive()) continue;
4507 if (s2->GetKinkIndexes()[0]!=0)
4508 if (s2->GetKinkIndexes()[0] == -s1->GetKinkIndexes()[0]) continue;
4509 AliExternalTrackParam &par1=params[index1];
4510 if (TMath::Abs(par0.GetParameter()[3]-par1.GetParameter()[3])>kCutP3) continue;
4511 if (TMath::Abs(par0.GetParameter()[1]-par1.GetParameter()[1])>kCutP1) continue;
4512 if (TMath::Abs(par0.GetParameter()[2]-par1.GetParameter()[2])>kCutP2) continue;
4513 Double_t dAlpha= TMath::Abs(par0.GetAlpha()-par1.GetAlpha());
4514 if (dAlpha>TMath::Pi()) dAlpha-=TMath::Pi();
4515 if (TMath::Abs(dAlpha)>kCutAlpha) continue;
4520 Int_t firstShared=lastpoint, lastShared=firstpoint;
4521 Int_t firstRow=lastpoint, lastRow=firstpoint;
4523 for (Int_t i=firstpoint;i<lastpoint;i++){
4524 if (s1->GetClusterIndex2(i)>0) nall0++;
4525 if (s2->GetClusterIndex2(i)>0) nall1++;
4526 if (s1->GetClusterIndex2(i)>0 && s2->GetClusterIndex2(i)>0) {
4527 if (i<firstRow) firstRow=i;
4528 if (i>lastRow) lastRow=i;
4530 if ( (s1->GetClusterIndex2(i))==(s2->GetClusterIndex2(i)) && s1->GetClusterIndex2(i)>0) {
4531 if (i<firstShared) firstShared=i;
4532 if (i>lastShared) lastShared=i;
4536 Double_t ratio0 = Float_t(sumShared)/Float_t(TMath::Min(nall0+1,nall1+1));
4537 Double_t ratio1 = Float_t(sumShared)/Float_t(TMath::Max(nall0+1,nall1+1));
4539 if( AliTPCReconstructor::StreamLevel()>1){
4540 TTreeSRedirector &cstream = *fDebugStreamer;
4541 Int_t n0=s1->GetNumberOfClusters();
4542 Int_t n1=s2->GetNumberOfClusters();
4543 Int_t n0F=s1->GetNFoundable();
4544 Int_t n1F=s2->GetNFoundable();
4545 Int_t lab0=s1->GetLabel();
4546 Int_t lab1=s2->GetLabel();
4548 cstream<<"Splitted2"<<
4549 "iter="<<fIteration<<
4550 "lab0="<<lab0<< // MC label if exist
4551 "lab1="<<lab1<< // MC label if exist
4554 "ratio0="<<ratio0<< // shared ratio
4555 "ratio1="<<ratio1<< // shared ratio
4556 "p0.="<<&par0<< // track parameters
4558 "s0.="<<s1<< // full seed
4560 "n0="<<n0<< // number of clusters track 0
4561 "n1="<<n1<< // number of clusters track 1
4562 "nall0="<<nall0<< // number of clusters track 0
4563 "nall1="<<nall1<< // number of clusters track 1
4564 "n0F="<<n0F<< // number of findable
4565 "n1F="<<n1F<< // number of findable
4566 "shared="<<sumShared<< // number of shared clusters
4567 "firstS="<<firstShared<< // first and the last shared row
4568 "lastS="<<lastShared<<
4569 "firstRow="<<firstRow<< // first and the last row with cluster
4570 "lastRow="<<lastRow<< //
4574 // remove track with lower quality
4576 if (ratio0>AliTPCReconstructor::GetRecoParam()->GetCutSharedClusters(0) ||
4577 ratio1>AliTPCReconstructor::GetRecoParam()->GetCutSharedClusters(1)){
4581 delete array->RemoveAt(index1);
4586 // 4. Delete temporary array
4596 void AliTPCtrackerMI::FindCurling(const TObjArray * array, AliESDEvent */*esd*/, Int_t iter)
4599 // find Curling tracks
4600 // Use AliTPCReconstructor::StreamLevel()>1 if you want to tune parameters - cuts
4603 // Algorithm done in 2 phases - because of CPU consumption
4604 // it is n^2 algorithm - for lead-lead 20000x20000 combination are investigated
4605 // see detal in MC part what can be used to cut
4609 const Float_t kMaxC = 400; // maximal curvature to of the track
4610 const Float_t kMaxdTheta = 0.15; // maximal distance in theta
4611 const Float_t kMaxdPhi = 0.15; // maximal distance in phi
4612 const Float_t kPtRatio = 0.3; // ratio between pt
4613 const Float_t kMinDCAR = 2.; // distance to the primary vertex in r - see cpipe cut
4616 // Curling tracks cuts
4619 const Float_t kMaxDeltaRMax = 40; // distance in outer radius
4620 const Float_t kMaxDeltaRMin = 5.; // distance in lower radius - see cpipe cut
4621 const Float_t kMinAngle = 2.9; // angle between tracks
4622 const Float_t kMaxDist = 5; // biggest distance
4624 // The cuts can be tuned using the "MC information stored in Multi tree ==> see FindMultiMC
4627 TCut csign("csign","Tr0.fP[4]*Tr1.fP[4]<0"); //opposite sign
4628 TCut cmax("cmax","abs(Tr0.GetC())>1/400");
4629 TCut cda("cda","sqrt(dtheta^2+dfi^2)<0.15");
4630 TCut ccratio("ccratio","abs((Tr0.fP[4]+Tr1.fP[4])/(abs(Tr0.fP[4])+abs(Tr1.fP[4])))<0.3");
4631 TCut cpipe("cpipe", "min(abs(r0-rc0),abs(r1-rc1))>5");
4633 TCut cdrmax("cdrmax","abs(abs(rc0+r0)-abs(rc1+r1))<40")
4634 TCut cdrmin("cdrmin","abs(abs(rc0+r0)-abs(rc1+r1))<10")
4636 Multi->Draw("dfi","iter==0"+csign+cmax+cda+ccratio); ~94% of curling tracks fulfill
4637 Multi->Draw("min(abs(r0-rc0),abs(r1-rc1))","iter==0&&abs(lab1)==abs(lab0)"+csign+cmax+cda+ccratio+cpipe+cdrmin+cdrmax); //80%
4639 Curling2->Draw("dfi","iter==0&&abs(lab0)==abs(lab1)"+csign+cmax+cdtheta+cdfi+ccratio)
4645 Int_t nentries = array->GetEntriesFast();
4646 AliHelix *helixes = new AliHelix[nentries];
4647 for (Int_t i=0;i<nentries;i++){
4648 AliTPCseed* track = (AliTPCseed*)array->At(i);
4649 if (!track) continue;
4650 track->SetCircular(0);
4651 new (&helixes[i]) AliHelix(*track);
4657 Double_t phase[2][2]={{0,0},{0,0}},radius[2]={0,0};
4663 for (Int_t i0=0;i0<nentries;i0++){
4664 AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
4665 if (!track0) continue;
4666 if (TMath::Abs(track0->GetC())<1/kMaxC) continue;
4667 Float_t xc0 = helixes[i0].GetHelix(6);
4668 Float_t yc0 = helixes[i0].GetHelix(7);
4669 Float_t r0 = helixes[i0].GetHelix(8);
4670 Float_t rc0 = TMath::Sqrt(xc0*xc0+yc0*yc0);
4671 Float_t fi0 = TMath::ATan2(yc0,xc0);
4673 for (Int_t i1=i0+1;i1<nentries;i1++){
4674 AliTPCseed * track1 = (AliTPCseed*)array->At(i1);
4675 if (!track1) continue;
4676 if (TMath::Abs(track1->GetC())<1/kMaxC) continue;
4677 Float_t xc1 = helixes[i1].GetHelix(6);
4678 Float_t yc1 = helixes[i1].GetHelix(7);
4679 Float_t r1 = helixes[i1].GetHelix(8);
4680 Float_t rc1 = TMath::Sqrt(xc1*xc1+yc1*yc1);
4681 Float_t fi1 = TMath::ATan2(yc1,xc1);
4683 Float_t dfi = fi0-fi1;
4686 if (dfi>1.5*TMath::Pi()) dfi-=TMath::Pi(); // take care about edge effect
4687 if (dfi<-1.5*TMath::Pi()) dfi+=TMath::Pi(); //
4688 Float_t dtheta = TMath::Abs(track0->GetTgl()-track1->GetTgl())<TMath::Abs(track0->GetTgl()+track1->GetTgl())? track0->GetTgl()-track1->GetTgl():track0->GetTgl()+track1->GetTgl();
4692 if (track0->GetBConstrain()&&track1->GetBConstrain()) continue; // not constrained
4693 if (track1->GetSigned1Pt()*track0->GetSigned1Pt()>0) continue; // not the same sign
4694 if ( TMath::Abs(track1->GetTgl()+track0->GetTgl())>kMaxdTheta) continue; //distance in the Theta
4695 if ( TMath::Abs(dfi)>kMaxdPhi) continue; //distance in phi
4696 if ( TMath::Sqrt(dfi*dfi+dtheta*dtheta)>kMaxdPhi) continue; //common angular offset
4698 Float_t pt0 = track0->GetSignedPt();
4699 Float_t pt1 = track1->GetSignedPt();
4700 if ((TMath::Abs(pt0+pt1)/(TMath::Abs(pt0)+TMath::Abs(pt1)))>kPtRatio) continue;
4701 if ((iter==1) && TMath::Abs(TMath::Abs(rc0+r0)-TMath::Abs(rc1+r1))>kMaxDeltaRMax) continue;
4702 if ((iter!=1) &&TMath::Abs(TMath::Abs(rc0-r0)-TMath::Abs(rc1-r1))>kMaxDeltaRMin) continue;
4703 if (TMath::Min(TMath::Abs(rc0-r0),TMath::Abs(rc1-r1))<kMinDCAR) continue;
4706 // Now find closest approach
4710 Int_t npoints = helixes[i0].GetRPHIintersections(helixes[i1], phase, radius,10);
4711 if (npoints==0) continue;
4712 helixes[i0].GetClosestPhases(helixes[i1], phase);
4716 Double_t hangles[3];
4717 helixes[i0].Evaluate(phase[0][0],xyz0);
4718 helixes[i1].Evaluate(phase[0][1],xyz1);
4720 helixes[i0].GetAngle(phase[0][0],helixes[i1],phase[0][1],hangles);
4721 Double_t deltah[2],deltabest;
4722 if (TMath::Abs(hangles[2])<kMinAngle) continue;
4726 helixes[i0].ParabolicDCA(helixes[i1],phase[0][0],phase[0][1],radius[0],deltah[0],2);
4728 helixes[i0].ParabolicDCA(helixes[i1],phase[1][0],phase[1][1],radius[1],deltah[1],2);
4729 if (deltah[1]<deltah[0]) ibest=1;
4731 deltabest = TMath::Sqrt(deltah[ibest]);
4732 helixes[i0].Evaluate(phase[ibest][0],xyz0);
4733 helixes[i1].Evaluate(phase[ibest][1],xyz1);
4734 helixes[i0].GetAngle(phase[ibest][0],helixes[i1],phase[ibest][1],hangles);
4735 Double_t radiusbest = TMath::Sqrt(radius[ibest]);
4737 if (deltabest>kMaxDist) continue;
4738 // if (mindcar+mindcaz<40 && (TMath::Abs(hangles[2])<kMinAngle ||deltabest>3)) continue;
4739 Bool_t sign =kFALSE;
4740 if (hangles[2]>kMinAngle) sign =kTRUE;
4743 // circular[i0] = kTRUE;
4744 // circular[i1] = kTRUE;
4745 if (track0->OneOverPt()<track1->OneOverPt()){
4746 track0->SetCircular(track0->GetCircular()+1);
4747 track1->SetCircular(track1->GetCircular()+2);
4750 track1->SetCircular(track1->GetCircular()+1);
4751 track0->SetCircular(track0->GetCircular()+2);
4754 if (AliTPCReconstructor::StreamLevel()>2){
4756 //debug stream to tune "fine" cuts
4757 Int_t lab0=track0->GetLabel();
4758 Int_t lab1=track1->GetLabel();
4759 TTreeSRedirector &cstream = *fDebugStreamer;
4760 cstream<<"Curling2"<<
4776 "npoints="<<npoints<<
4777 "hangles0="<<hangles[0]<<
4778 "hangles1="<<hangles[1]<<
4779 "hangles2="<<hangles[2]<<
4782 "radius="<<radiusbest<<
4783 "deltabest="<<deltabest<<
4784 "phase0="<<phase[ibest][0]<<
4785 "phase1="<<phase[ibest][1]<<
4793 if (AliTPCReconstructor::StreamLevel()>1) {
4794 AliInfo("Time for curling tracks removal");
4803 void AliTPCtrackerMI::FindKinks(TObjArray * array, AliESDEvent *esd)
4810 TObjArray *kinks= new TObjArray(10000);
4811 // TObjArray *v0s= new TObjArray(10000);
4812 Int_t nentries = array->GetEntriesFast();
4813 AliHelix *helixes = new AliHelix[nentries];
4814 Int_t *sign = new Int_t[nentries];
4815 Int_t *nclusters = new Int_t[nentries];
4816 Float_t *alpha = new Float_t[nentries];
4817 AliKink *kink = new AliKink();
4818 Int_t * usage = new Int_t[nentries];
4819 Float_t *zm = new Float_t[nentries];
4820 Float_t *z0 = new Float_t[nentries];
4821 Float_t *fim = new Float_t[nentries];
4822 Float_t *shared = new Float_t[nentries];
4823 Bool_t *circular = new Bool_t[nentries];
4824 Float_t *dca = new Float_t[nentries];
4825 //const AliESDVertex * primvertex = esd->GetVertex();
4827 // nentries = array->GetEntriesFast();
4832 for (Int_t i=0;i<nentries;i++){
4835 AliTPCseed* track = (AliTPCseed*)array->At(i);
4836 if (!track) continue;
4837 track->SetCircular(0);
4839 track->UpdatePoints();
4840 if (( track->GetPoints()[2]- track->GetPoints()[0])>5 && track->GetPoints()[3]>0.8){
4842 nclusters[i]=track->GetNumberOfClusters();
4843 alpha[i] = track->GetAlpha();
4844 new (&helixes[i]) AliHelix(*track);
4846 helixes[i].Evaluate(0,xyz);
4847 sign[i] = (track->GetC()>0) ? -1:1;
4850 if (track->GetProlongation(x,y,z)){
4852 fim[i] = alpha[i]+TMath::ATan2(y,x);
4855 zm[i] = track->GetZ();
4859 circular[i]= kFALSE;
4860 if (track->GetProlongation(0,y,z)) z0[i] = z;
4861 dca[i] = track->GetD(0,0);
4867 Int_t ncandidates =0;
4870 Double_t phase[2][2]={{0,0},{0,0}},radius[2]={0,0};
4873 // Find circling track
4875 for (Int_t i0=0;i0<nentries;i0++){
4876 AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
4877 if (!track0) continue;
4878 if (track0->GetNumberOfClusters()<40) continue;
4879 if (TMath::Abs(1./track0->GetC())>200) continue;
4880 for (Int_t i1=i0+1;i1<nentries;i1++){
4881 AliTPCseed * track1 = (AliTPCseed*)array->At(i1);
4882 if (!track1) continue;
4883 if (track1->GetNumberOfClusters()<40) continue;
4884 if ( TMath::Abs(track1->GetTgl()+track0->GetTgl())>0.1) continue;
4885 if (track0->GetBConstrain()&&track1->GetBConstrain()) continue;
4886 if (TMath::Abs(1./track1->GetC())>200) continue;
4887 if (track1->GetSigned1Pt()*track0->GetSigned1Pt()>0) continue;
4888 if (track1->GetTgl()*track0->GetTgl()>0) continue;
4889 if (TMath::Max(TMath::Abs(1./track0->GetC()),TMath::Abs(1./track1->GetC()))>190) continue;
4890 if (track0->GetBConstrain()&&track1->OneOverPt()<track0->OneOverPt()) continue; //returning - lower momenta
4891 if (track1->GetBConstrain()&&track0->OneOverPt()<track1->OneOverPt()) continue; //returning - lower momenta
4893 Float_t mindcar = TMath::Min(TMath::Abs(dca[i0]),TMath::Abs(dca[i1]));
4894 if (mindcar<5) continue;
4895 Float_t mindcaz = TMath::Min(TMath::Abs(z0[i0]-GetZ()),TMath::Abs(z0[i1]-GetZ()));
4896 if (mindcaz<5) continue;
4897 if (mindcar+mindcaz<20) continue;
4900 Float_t xc0 = helixes[i0].GetHelix(6);
4901 Float_t yc0 = helixes[i0].GetHelix(7);
4902 Float_t r0 = helixes[i0].GetHelix(8);
4903 Float_t xc1 = helixes[i1].GetHelix(6);
4904 Float_t yc1 = helixes[i1].GetHelix(7);
4905 Float_t r1 = helixes[i1].GetHelix(8);
4907 Float_t rmean = (r0+r1)*0.5;
4908 Float_t delta =TMath::Sqrt((xc1-xc0)*(xc1-xc0)+(yc1-yc0)*(yc1-yc0));
4909 //if (delta>30) continue;
4910 if (delta>rmean*0.25) continue;
4911 if (TMath::Abs(r0-r1)/rmean>0.3) continue;
4913 Int_t npoints = helixes[i0].GetRPHIintersections(helixes[i1], phase, radius,10);
4914 if (npoints==0) continue;
4915 helixes[i0].GetClosestPhases(helixes[i1], phase);
4919 Double_t hangles[3];
4920 helixes[i0].Evaluate(phase[0][0],xyz0);
4921 helixes[i1].Evaluate(phase[0][1],xyz1);
4923 helixes[i0].GetAngle(phase[0][0],helixes[i1],phase[0][1],hangles);
4924 Double_t deltah[2],deltabest;
4925 if (hangles[2]<2.8) continue;
4928 helixes[i0].ParabolicDCA(helixes[i1],phase[0][0],phase[0][1],radius[0],deltah[0],2);
4930 helixes[i0].ParabolicDCA(helixes[i1],phase[1][0],phase[1][1],radius[1],deltah[1],2);
4931 if (deltah[1]<deltah[0]) ibest=1;
4933 deltabest = TMath::Sqrt(deltah[ibest]);
4934 helixes[i0].Evaluate(phase[ibest][0],xyz0);
4935 helixes[i1].Evaluate(phase[ibest][1],xyz1);
4936 helixes[i0].GetAngle(phase[ibest][0],helixes[i1],phase[ibest][1],hangles);
4937 Double_t radiusbest = TMath::Sqrt(radius[ibest]);
4939 if (deltabest>6) continue;
4940 if (mindcar+mindcaz<40 && (hangles[2]<3.12||deltabest>3)) continue;
4941 Bool_t lsign =kFALSE;
4942 if (hangles[2]>3.06) lsign =kTRUE;
4945 circular[i0] = kTRUE;
4946 circular[i1] = kTRUE;
4947 if (track0->OneOverPt()<track1->OneOverPt()){
4948 track0->SetCircular(track0->GetCircular()+1);
4949 track1->SetCircular(track1->GetCircular()+2);
4952 track1->SetCircular(track1->GetCircular()+1);
4953 track0->SetCircular(track0->GetCircular()+2);
4956 if (lsign&&AliTPCReconstructor::StreamLevel()>1){
4958 Int_t lab0=track0->GetLabel();
4959 Int_t lab1=track1->GetLabel();
4960 TTreeSRedirector &cstream = *fDebugStreamer;
4961 cstream<<"Curling"<<
4968 "mindcar="<<mindcar<<
4969 "mindcaz="<<mindcaz<<
4972 "npoints="<<npoints<<
4973 "hangles0="<<hangles[0]<<
4974 "hangles2="<<hangles[2]<<
4979 "radius="<<radiusbest<<
4980 "deltabest="<<deltabest<<
4981 "phase0="<<phase[ibest][0]<<
4982 "phase1="<<phase[ibest][1]<<
4992 for (Int_t i =0;i<nentries;i++){
4993 if (sign[i]==0) continue;
4994 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5001 Double_t cradius0 = 40*40;
5002 Double_t cradius1 = 270*270;
5005 Double_t cdist3=0.55;
5006 for (Int_t j =i+1;j<nentries;j++){
5008 if (sign[j]*sign[i]<1) continue;
5009 if ( (nclusters[i]+nclusters[j])>200) continue;
5010 if ( (nclusters[i]+nclusters[j])<80) continue;
5011 if ( TMath::Abs(zm[i]-zm[j])>60.) continue;
5012 if ( TMath::Abs(fim[i]-fim[j])>0.6 && TMath::Abs(fim[i]-fim[j])<5.7 ) continue;
5013 //AliTPCseed * track1 = (AliTPCseed*)array->At(j); Double_t phase[2][2],radius[2];
5014 Int_t npoints = helixes[i].GetRPHIintersections(helixes[j], phase, radius,20);
5015 if (npoints<1) continue;
5018 if (radius[0]<cradius0||radius[0]>cradius1) continue;
5021 if ( (radius[0]<cradius0||radius[0]>cradius1) && (radius[1]<cradius0||radius[1]>cradius1) ) continue;
5024 Double_t delta1=10000,delta2=10000;
5025 // cuts on the intersection radius
5026 helixes[i].LinearDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta1);
5027 if (radius[0]<20&&delta1<1) continue; //intersection at vertex
5028 if (radius[0]<10&&delta1<3) continue; //intersection at vertex
5030 helixes[i].LinearDCA(helixes[j],phase[1][0],phase[1][1],radius[1],delta2);
5031 if (radius[1]<20&&delta2<1) continue; //intersection at vertex
5032 if (radius[1]<10&&delta2<3) continue; //intersection at vertex
5035 Double_t distance1 = TMath::Min(delta1,delta2);
5036 if (distance1>cdist1) continue; // cut on DCA linear approximation
5038 npoints = helixes[i].GetRPHIintersections(helixes[j], phase, radius,20);
5039 helixes[i].ParabolicDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta1);
5040 if (radius[0]<20&&delta1<1) continue; //intersection at vertex
5041 if (radius[0]<10&&delta1<3) continue; //intersection at vertex
5044 helixes[i].ParabolicDCA(helixes[j],phase[1][0],phase[1][1],radius[1],delta2);
5045 if (radius[1]<20&&delta2<1) continue; //intersection at vertex
5046 if (radius[1]<10&&delta2<3) continue; //intersection at vertex
5048 distance1 = TMath::Min(delta1,delta2);
5051 rkink = TMath::Sqrt(radius[0]);
5054 rkink = TMath::Sqrt(radius[1]);
5056 if (distance1>cdist2) continue;
5059 AliTPCseed * track1 = (AliTPCseed*)array->At(j);
5062 Int_t row0 = GetRowNumber(rkink);
5063 if (row0<10) continue;
5064 if (row0>150) continue;
5067 Float_t dens00=-1,dens01=-1;
5068 Float_t dens10=-1,dens11=-1;
5070 Int_t found,foundable,ishared;
5071 track0->GetClusterStatistic(0,row0-5, found, foundable,ishared,kFALSE);
5072 if (foundable>5) dens00 = Float_t(found)/Float_t(foundable);
5073 track0->GetClusterStatistic(row0+5,155, found, foundable,ishared,kFALSE);
5074 if (foundable>5) dens01 = Float_t(found)/Float_t(foundable);
5076 track1->GetClusterStatistic(0,row0-5, found, foundable,ishared,kFALSE);
5077 if (foundable>10) dens10 = Float_t(found)/Float_t(foundable);
5078 track1->GetClusterStatistic(row0+5,155, found, foundable,ishared,kFALSE);
5079 if (foundable>10) dens11 = Float_t(found)/Float_t(foundable);
5081 if (dens00<dens10 && dens01<dens11) continue;
5082 if (dens00>dens10 && dens01>dens11) continue;
5083 if (TMath::Max(dens00,dens10)<0.1) continue;
5084 if (TMath::Max(dens01,dens11)<0.3) continue;
5086 if (TMath::Min(dens00,dens10)>0.6) continue;
5087 if (TMath::Min(dens01,dens11)>0.6) continue;
5090 AliTPCseed * ktrack0, *ktrack1;
5099 if (TMath::Abs(ktrack0->GetC())>5) continue; // cut on the curvature for mother particle
5100 AliExternalTrackParam paramm(*ktrack0);
5101 AliExternalTrackParam paramd(*ktrack1);
5102 if (row0>60&&ktrack1->GetReference().GetX()>90.)new (¶md) AliExternalTrackParam(ktrack1->GetReference());
5105 kink->SetMother(paramm);
5106 kink->SetDaughter(paramd);
5109 Float_t x[3] = { kink->GetPosition()[0],kink->GetPosition()[1],kink->GetPosition()[2]};
5111 fkParam->Transform0to1(x,index);
5112 fkParam->Transform1to2(x,index);
5113 row0 = GetRowNumber(x[0]);
5115 if (kink->GetR()<100) continue;
5116 if (kink->GetR()>240) continue;
5117 if (kink->GetPosition()[2]/kink->GetR()>AliTPCReconstructor::GetCtgRange()) continue; //out of fiducial volume
5118 if (kink->GetDistance()>cdist3) continue;
5119 Float_t dird = kink->GetDaughterP()[0]*kink->GetPosition()[0]+kink->GetDaughterP()[1]*kink->GetPosition()[1]; // rough direction estimate
5120 if (dird<0) continue;
5122 Float_t dirm = kink->GetMotherP()[0]*kink->GetPosition()[0]+kink->GetMotherP()[1]*kink->GetPosition()[1]; // rough direction estimate
5123 if (dirm<0) continue;
5124 Float_t mpt = TMath::Sqrt(kink->GetMotherP()[0]*kink->GetMotherP()[0]+kink->GetMotherP()[1]*kink->GetMotherP()[1]);
5125 if (mpt<0.2) continue;
5128 //for high momenta momentum not defined well in first iteration
5129 Double_t qt = TMath::Sin(kink->GetAngle(2))*ktrack1->GetP();
5130 if (qt>0.35) continue;
5133 kink->SetLabel(CookLabel(ktrack0,0.4,0,row0),0);
5134 kink->SetLabel(CookLabel(ktrack1,0.4,row0,160),1);
5136 kink->SetTPCDensity(dens00,0,0);
5137 kink->SetTPCDensity(dens01,0,1);
5138 kink->SetTPCDensity(dens10,1,0);
5139 kink->SetTPCDensity(dens11,1,1);
5140 kink->SetIndex(i,0);
5141 kink->SetIndex(j,1);
5144 kink->SetTPCDensity(dens10,0,0);
5145 kink->SetTPCDensity(dens11,0,1);
5146 kink->SetTPCDensity(dens00,1,0);
5147 kink->SetTPCDensity(dens01,1,1);
5148 kink->SetIndex(j,0);
5149 kink->SetIndex(i,1);
5152 if (mpt<1||kink->GetAngle(2)>0.1){
5153 // angle and densities not defined yet
5154 if (kink->GetTPCDensityFactor()<0.8) continue;
5155 if ((2-kink->GetTPCDensityFactor())*kink->GetDistance() >0.25) continue;
5156 if (kink->GetAngle(2)*ktrack0->GetP()<0.003) continue; //too small angle
5157 if (kink->GetAngle(2)>0.2&&kink->GetTPCDensityFactor()<1.15) continue;
5158 if (kink->GetAngle(2)>0.2&&kink->GetTPCDensity(0,1)>0.05) continue;
5160 Float_t criticalangle = track0->GetSigmaSnp2()+track0->GetSigmaTgl2();
5161 criticalangle+= track1->GetSigmaSnp2()+track1->GetSigmaTgl2();
5162 criticalangle= 3*TMath::Sqrt(criticalangle);
5163 if (criticalangle>0.02) criticalangle=0.02;
5164 if (kink->GetAngle(2)<criticalangle) continue;
5167 Int_t drow = Int_t(2.+0.5/(0.05+kink->GetAngle(2))); // overlap region defined
5168 Float_t shapesum =0;
5170 for ( Int_t row = row0-drow; row<row0+drow;row++){
5171 if (row<0) continue;
5172 if (row>155) continue;
5173 if (ktrack0->GetClusterPointer(row)){
5174 AliTPCTrackerPoint *point =ktrack0->GetTrackPoint(row);
5175 shapesum+=point->GetSigmaY()+point->GetSigmaZ();
5178 if (ktrack1->GetClusterPointer(row)){
5179 AliTPCTrackerPoint *point =ktrack1->GetTrackPoint(row);
5180 shapesum+=point->GetSigmaY()+point->GetSigmaZ();
5185 kink->SetShapeFactor(-1.);
5188 kink->SetShapeFactor(shapesum/sum);
5190 // esd->AddKink(kink);
5192 // kink->SetMother(paramm);
5193 //kink->SetDaughter(paramd);
5195 Double_t chi2P2 = paramm.GetParameter()[2]-paramd.GetParameter()[2];
5197 chi2P2/=paramm.GetCovariance()[5]+paramd.GetCovariance()[5];
5198 Double_t chi2P3 = paramm.GetParameter()[3]-paramd.GetParameter()[3];
5200 chi2P3/=paramm.GetCovariance()[9]+paramd.GetCovariance()[9];
5202 if (AliTPCReconstructor::StreamLevel()>1) {
5203 (*fDebugStreamer)<<"kinkLpt"<<
5211 if ( chi2P2+chi2P3<AliTPCReconstructor::GetRecoParam()->GetKinkAngleCutChi2(0)){
5215 kinks->AddLast(kink);
5221 // sort the kinks according quality - and refit them towards vertex
5223 Int_t nkinks = kinks->GetEntriesFast();
5224 Float_t *quality = new Float_t[nkinks];
5225 Int_t *indexes = new Int_t[nkinks];
5226 AliTPCseed *mothers = new AliTPCseed[nkinks];
5227 AliTPCseed *daughters = new AliTPCseed[nkinks];
5230 for (Int_t i=0;i<nkinks;i++){
5232 AliKink *kinkl = (AliKink*)kinks->At(i);
5234 // refit kinks towards vertex
5236 Int_t index0 = kinkl->GetIndex(0);
5237 Int_t index1 = kinkl->GetIndex(1);
5238 AliTPCseed * ktrack0 = (AliTPCseed*)array->At(index0);
5239 AliTPCseed * ktrack1 = (AliTPCseed*)array->At(index1);
5241 Int_t sumn=ktrack0->GetNumberOfClusters()+ktrack1->GetNumberOfClusters();
5243 // Refit Kink under if too small angle
5245 if (kinkl->GetAngle(2)<0.05){
5246 kinkl->SetTPCRow0(GetRowNumber(kinkl->GetR()));
5247 Int_t row0 = kinkl->GetTPCRow0();
5248 Int_t drow = Int_t(2.+0.5/(0.05+kinkl->GetAngle(2)));
5251 Int_t last = row0-drow;
5252 if (last<40) last=40;
5253 if (last<ktrack0->GetFirstPoint()+25) last = ktrack0->GetFirstPoint()+25;
5254 AliTPCseed* seed0 = ReSeed(ktrack0,last,kFALSE);
5257 Int_t first = row0+drow;
5258 if (first>130) first=130;
5259 if (first>ktrack1->GetLastPoint()-25) first = TMath::Max(ktrack1->GetLastPoint()-25,30);
5260 AliTPCseed* seed1 = ReSeed(ktrack1,first,kTRUE);
5262 if (seed0 && seed1){
5263 kinkl->SetStatus(1,8);
5264 if (RefitKink(*seed0,*seed1,*kinkl)) kinkl->SetStatus(1,9);
5265 row0 = GetRowNumber(kinkl->GetR());
5266 sumn = seed0->GetNumberOfClusters()+seed1->GetNumberOfClusters();
5267 mothers[i] = *seed0;
5268 daughters[i] = *seed1;
5271 delete kinks->RemoveAt(i);
5272 if (seed0) delete seed0;
5273 if (seed1) delete seed1;
5276 if (kinkl->GetDistance()>0.5 || kinkl->GetR()<110 || kinkl->GetR()>240) {
5277 delete kinks->RemoveAt(i);
5278 if (seed0) delete seed0;
5279 if (seed1) delete seed1;
5287 if (kinkl) quality[i] = 160*((0.1+kinkl->GetDistance())*(2.-kinkl->GetTPCDensityFactor()))/(sumn+40.); //the longest -clossest will win
5289 TMath::Sort(nkinks,quality,indexes,kFALSE);
5291 //remove double find kinks
5293 for (Int_t ikink0=1;ikink0<nkinks;ikink0++){
5294 AliKink * kink0 = (AliKink*) kinks->At(indexes[ikink0]);
5295 if (!kink0) continue;
5297 for (Int_t ikink1=0;ikink1<ikink0;ikink1++){
5298 kink0 = (AliKink*) kinks->At(indexes[ikink0]);
5299 if (!kink0) continue;
5300 AliKink * kink1 = (AliKink*) kinks->At(indexes[ikink1]);
5301 if (!kink1) continue;
5302 // if not close kink continue
5303 if (TMath::Abs(kink1->GetPosition()[2]-kink0->GetPosition()[2])>10) continue;
5304 if (TMath::Abs(kink1->GetPosition()[1]-kink0->GetPosition()[1])>10) continue;
5305 if (TMath::Abs(kink1->GetPosition()[0]-kink0->GetPosition()[0])>10) continue;
5307 AliTPCseed &mother0 = mothers[indexes[ikink0]];
5308 AliTPCseed &daughter0 = daughters[indexes[ikink0]];
5309 AliTPCseed &mother1 = mothers[indexes[ikink1]];
5310 AliTPCseed &daughter1 = daughters[indexes[ikink1]];
5311 Int_t row0 = (kink0->GetTPCRow0()+kink1->GetTPCRow0())/2;
5320 for (Int_t i=0;i<row0;i++){
5321 if (mother0.GetClusterIndex(i)>0 && mother1.GetClusterIndex(i)>0){
5324 if (mother0.GetClusterIndex(i)==mother1.GetClusterIndex(i)){
5331 for (Int_t i=row0;i<158;i++){
5332 if (daughter0.GetClusterIndex(i)>0 && daughter0.GetClusterIndex(i)>0){
5335 if (mother0.GetClusterIndex(i)==mother1.GetClusterIndex(i)){
5341 Float_t ratio = Float_t(same+1)/Float_t(both+1);
5342 Float_t ratiom = Float_t(samem+1)/Float_t(bothm+1);
5343 Float_t ratiod = Float_t(samed+1)/Float_t(bothd+1);
5344 if (ratio>0.3 && ratiom>0.5 &&ratiod>0.5) {
5345 Int_t sum0 = mother0.GetNumberOfClusters()+daughter0.GetNumberOfClusters();
5346 Int_t sum1 = mother1.GetNumberOfClusters()+daughter1.GetNumberOfClusters();
5348 shared[kink0->GetIndex(0)]= kTRUE;
5349 shared[kink0->GetIndex(1)]= kTRUE;
5350 delete kinks->RemoveAt(indexes[ikink0]);
5354 shared[kink1->GetIndex(0)]= kTRUE;
5355 shared[kink1->GetIndex(1)]= kTRUE;
5356 delete kinks->RemoveAt(indexes[ikink1]);
5363 for (Int_t i=0;i<nkinks;i++){
5364 AliKink * kinkl = (AliKink*) kinks->At(indexes[i]);
5365 if (!kinkl) continue;
5366 kinkl->SetTPCRow0(GetRowNumber(kinkl->GetR()));
5367 Int_t index0 = kinkl->GetIndex(0);
5368 Int_t index1 = kinkl->GetIndex(1);
5369 if (circular[index0]||(circular[index1]&&kinkl->GetDistance()>0.2)) continue;
5370 kinkl->SetMultiple(usage[index0],0);
5371 kinkl->SetMultiple(usage[index1],1);
5372 if (kinkl->GetMultiple()[0]+kinkl->GetMultiple()[1]>2) continue;
5373 if (kinkl->GetMultiple()[0]+kinkl->GetMultiple()[1]>0 && quality[indexes[i]]>0.2) continue;
5374 if (kinkl->GetMultiple()[0]+kinkl->GetMultiple()[1]>0 && kinkl->GetDistance()>0.2) continue;
5375 if (circular[index0]||(circular[index1]&&kinkl->GetDistance()>0.1)) continue;
5377 AliTPCseed * ktrack0 = (AliTPCseed*)array->At(index0);
5378 AliTPCseed * ktrack1 = (AliTPCseed*)array->At(index1);
5379 if (!ktrack0 || !ktrack1) continue;
5380 Int_t index = esd->AddKink(kinkl);
5383 if ( ktrack0->GetKinkIndex(0)==0 && ktrack1->GetKinkIndex(0)==0) { //best kink
5384 if (mothers[indexes[i]].GetNumberOfClusters()>20 && daughters[indexes[i]].GetNumberOfClusters()>20 && (mothers[indexes[i]].GetNumberOfClusters()+daughters[indexes[i]].GetNumberOfClusters())>100){
5385 *ktrack0 = mothers[indexes[i]];
5386 *ktrack1 = daughters[indexes[i]];
5390 ktrack0->SetKinkIndex(usage[index0],-(index+1));
5391 ktrack1->SetKinkIndex(usage[index1], (index+1));
5396 // Remove tracks corresponding to shared kink's
5398 for (Int_t i=0;i<nentries;i++){
5399 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5400 if (!track0) continue;
5401 if (track0->GetKinkIndex(0)!=0) continue;
5402 if (shared[i]) delete array->RemoveAt(i);
5407 RemoveUsed2(array,0.5,0.4,30);
5409 for (Int_t i=0;i<nentries;i++){
5410 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5411 if (!track0) continue;
5412 track0->CookdEdx(0.02,0.6);
5416 for (Int_t i=0;i<nentries;i++){
5417 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5418 if (!track0) continue;
5419 if (track0->Pt()<1.4) continue;
5420 //remove double high momenta tracks - overlapped with kink candidates
5423 for (Int_t icl=track0->GetFirstPoint();icl<track0->GetLastPoint(); icl++){
5424 if (track0->GetClusterPointer(icl)!=0){
5426 if (track0->GetClusterPointer(icl)->IsUsed(10)) ishared++;
5429 if (Float_t(ishared+1)/Float_t(all+1)>0.5) {
5430 delete array->RemoveAt(i);
5434 if (track0->GetKinkIndex(0)!=0) continue;
5435 if (track0->GetNumberOfClusters()<80) continue;
5437 AliTPCseed *pmother = new AliTPCseed();
5438 AliTPCseed *pdaughter = new AliTPCseed();
5439 AliKink *pkink = new AliKink;
5441 AliTPCseed & mother = *pmother;
5442 AliTPCseed & daughter = *pdaughter;
5443 AliKink & kinkl = *pkink;
5444 if (CheckKinkPoint(track0,mother,daughter, kinkl)){
5445 if (mother.GetNumberOfClusters()<30||daughter.GetNumberOfClusters()<20) {
5449 continue; //too short tracks
5451 if (mother.Pt()<1.4) {
5457 Int_t row0= kinkl.GetTPCRow0();
5458 if (kinkl.GetDistance()>0.5 || kinkl.GetR()<110. || kinkl.GetR()>240.) {
5465 Int_t index = esd->AddKink(&kinkl);
5466 mother.SetKinkIndex(0,-(index+1));
5467 daughter.SetKinkIndex(0,index+1);
5468 if (mother.GetNumberOfClusters()>50) {
5469 delete array->RemoveAt(i);
5470 array->AddAt(new AliTPCseed(mother),i);
5473 array->AddLast(new AliTPCseed(mother));
5475 array->AddLast(new AliTPCseed(daughter));
5476 for (Int_t icl=0;icl<row0;icl++) {
5477 if (mother.GetClusterPointer(icl)) mother.GetClusterPointer(icl)->Use(20);
5480 for (Int_t icl=row0;icl<158;icl++) {
5481 if (daughter.GetClusterPointer(icl)) daughter.GetClusterPointer(icl)->Use(20);
5490 delete [] daughters;
5512 AliInfo(Form("Ncandidates=\t%d\t%d\t%d\t%d\n",esd->GetNumberOfKinks(),ncandidates,ntracks,nall));
5517 Int_t AliTPCtrackerMI::RefitKink(AliTPCseed &mother, AliTPCseed &daughter, const AliESDkink &knk)
5520 // refit kink towards to the vertex
5523 AliKink &kink=(AliKink &)knk;
5525 Int_t row0 = GetRowNumber(kink.GetR());
5526 FollowProlongation(mother,0);
5527 mother.Reset(kFALSE);
5529 FollowProlongation(daughter,row0);
5530 daughter.Reset(kFALSE);
5531 FollowBackProlongation(daughter,158);
5532 daughter.Reset(kFALSE);
5533 Int_t first = TMath::Max(row0-20,30);
5534 Int_t last = TMath::Min(row0+20,140);
5536 const Int_t kNdiv =5;
5537 AliTPCseed param0[kNdiv]; // parameters along the track
5538 AliTPCseed param1[kNdiv]; // parameters along the track
5539 AliKink kinks[kNdiv]; // corresponding kink parameters
5542 for (Int_t irow=0; irow<kNdiv;irow++){
5543 rows[irow] = first +((last-first)*irow)/(kNdiv-1);
5545 // store parameters along the track
5547 for (Int_t irow=0;irow<kNdiv;irow++){
5548 FollowBackProlongation(mother, rows[irow]);
5549 FollowProlongation(daughter,rows[kNdiv-1-irow]);
5550 param0[irow] = mother;
5551 param1[kNdiv-1-irow] = daughter;
5555 for (Int_t irow=0; irow<kNdiv-1;irow++){
5556 if (param0[irow].GetNumberOfClusters()<kNdiv||param1[irow].GetNumberOfClusters()<kNdiv) continue;
5557 kinks[irow].SetMother(param0[irow]);
5558 kinks[irow].SetDaughter(param1[irow]);
5559 kinks[irow].Update();
5562 // choose kink with best "quality"
5564 Double_t mindist = 10000;
5565 for (Int_t irow=0;irow<kNdiv;irow++){
5566 if (param0[irow].GetNumberOfClusters()<20||param1[irow].GetNumberOfClusters()<20) continue;
5567 if (TMath::Abs(kinks[irow].GetR())>240.) continue;
5568 if (TMath::Abs(kinks[irow].GetR())<100.) continue;
5570 Float_t normdist = TMath::Abs(param0[irow].GetX()-kinks[irow].GetR())*(0.1+kink.GetDistance());
5571 normdist/= (param0[irow].GetNumberOfClusters()+param1[irow].GetNumberOfClusters()+40.);
5572 if (normdist < mindist){
5578 if (index==-1) return 0;
5581 param0[index].Reset(kTRUE);
5582 FollowProlongation(param0[index],0);
5584 mother = param0[index];
5585 daughter = param1[index]; // daughter in vertex
5587 kink.SetMother(mother);
5588 kink.SetDaughter(daughter);
5590 kink.SetTPCRow0(GetRowNumber(kink.GetR()));
5591 kink.SetTPCncls(param0[index].GetNumberOfClusters(),0);
5592 kink.SetTPCncls(param1[index].GetNumberOfClusters(),1);
5593 kink.SetLabel(CookLabel(&mother,0.4, 0,kink.GetTPCRow0()),0);
5594 kink.SetLabel(CookLabel(&daughter,0.4, kink.GetTPCRow0(),160),1);
5595 mother.SetLabel(kink.GetLabel(0));
5596 daughter.SetLabel(kink.GetLabel(1));
5602 void AliTPCtrackerMI::UpdateKinkQualityM(AliTPCseed * seed){
5604 // update Kink quality information for mother after back propagation
5606 if (seed->GetKinkIndex(0)>=0) return;
5607 for (Int_t ikink=0;ikink<3;ikink++){
5608 Int_t index = seed->GetKinkIndex(ikink);
5609 if (index>=0) break;
5610 index = TMath::Abs(index)-1;
5611 AliESDkink * kink = fEvent->GetKink(index);
5612 kink->SetTPCDensity(-1,0,0);
5613 kink->SetTPCDensity(1,0,1);
5615 Int_t row0 = kink->GetTPCRow0() - 2 - Int_t( 0.5/ (0.05+kink->GetAngle(2)));
5616 if (row0<15) row0=15;
5618 Int_t row1 = kink->GetTPCRow0() + 2 + Int_t( 0.5/ (0.05+kink->GetAngle(2)));
5619 if (row1>145) row1=145;
5621 Int_t found,foundable,shared;
5622 seed->GetClusterStatistic(0,row0, found, foundable,shared,kFALSE);
5623 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),0,0);
5624 seed->GetClusterStatistic(row1,155, found, foundable,shared,kFALSE);
5625 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),0,1);
5630 void AliTPCtrackerMI::UpdateKinkQualityD(AliTPCseed * seed){
5632 // update Kink quality information for daughter after refit
5634 if (seed->GetKinkIndex(0)<=0) return;
5635 for (Int_t ikink=0;ikink<3;ikink++){
5636 Int_t index = seed->GetKinkIndex(ikink);
5637 if (index<=0) break;
5638 index = TMath::Abs(index)-1;
5639 AliESDkink * kink = fEvent->GetKink(index);
5640 kink->SetTPCDensity(-1,1,0);
5641 kink->SetTPCDensity(-1,1,1);
5643 Int_t row0 = kink->GetTPCRow0() -2 - Int_t( 0.5/ (0.05+kink->GetAngle(2)));
5644 if (row0<15) row0=15;
5646 Int_t row1 = kink->GetTPCRow0() +2 + Int_t( 0.5/ (0.05+kink->GetAngle(2)));
5647 if (row1>145) row1=145;
5649 Int_t found,foundable,shared;
5650 seed->GetClusterStatistic(0,row0, found, foundable,shared,kFALSE);
5651 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),1,0);
5652 seed->GetClusterStatistic(row1,155, found, foundable,shared,kFALSE);
5653 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),1,1);
5659 Int_t AliTPCtrackerMI::CheckKinkPoint(AliTPCseed*seed,AliTPCseed &mother, AliTPCseed &daughter, const AliESDkink &knk)
5662 // check kink point for given track
5663 // if return value=0 kink point not found
5664 // otherwise seed0 correspond to mother particle
5665 // seed1 correspond to daughter particle
5666 // kink parameter of kink point
5667 AliKink &kink=(AliKink &)knk;
5669 Int_t middlerow = (seed->GetFirstPoint()+seed->GetLastPoint())/2;
5670 Int_t first = seed->GetFirstPoint();
5671 Int_t last = seed->GetLastPoint();
5672 if (last-first<20) return 0; // shortest length - 2*30 = 60 pad-rows
5675 AliTPCseed *seed1 = ReSeed(seed,middlerow+20, kTRUE); //middle of chamber
5676 if (!seed1) return 0;
5677 FollowProlongation(*seed1,seed->GetLastPoint()-20);
5678 seed1->Reset(kTRUE);
5679 FollowProlongation(*seed1,158);
5680 seed1->Reset(kTRUE);
5681 last = seed1->GetLastPoint();
5683 AliTPCseed *seed0 = new AliTPCseed(*seed);
5684 seed0->Reset(kFALSE);
5687 AliTPCseed param0[20]; // parameters along the track
5688 AliTPCseed param1[20]; // parameters along the track
5689 AliKink kinks[20]; // corresponding kink parameters
5691 for (Int_t irow=0; irow<20;irow++){
5692 rows[irow] = first +((last-first)*irow)/19;
5694 // store parameters along the track
5696 for (Int_t irow=0;irow<20;irow++){
5697 FollowBackProlongation(*seed0, rows[irow]);
5698 FollowProlongation(*seed1,rows[19-irow]);
5699 param0[irow] = *seed0;
5700 param1[19-irow] = *seed1;
5704 for (Int_t irow=0; irow<19;irow++){
5705 kinks[irow].SetMother(param0[irow]);
5706 kinks[irow].SetDaughter(param1[irow]);
5707 kinks[irow].Update();
5710 // choose kink with biggest change of angle
5712 Double_t maxchange= 0;
5713 for (Int_t irow=1;irow<19;irow++){
5714 if (TMath::Abs(kinks[irow].GetR())>240.) continue;
5715 if (TMath::Abs(kinks[irow].GetR())<110.) continue;
5716 Float_t quality = TMath::Abs(kinks[irow].GetAngle(2))/(3.+TMath::Abs(kinks[irow].GetR()-param0[irow].GetX()));
5717 if ( quality > maxchange){
5718 maxchange = quality;
5725 if (index<0) return 0;
5727 Int_t row0 = GetRowNumber(kinks[index].GetR()); //row 0 estimate
5728 seed0 = new AliTPCseed(param0[index]);
5729 seed1 = new AliTPCseed(param1[index]);
5730 seed0->Reset(kFALSE);
5731 seed1->Reset(kFALSE);
5732 seed0->ResetCovariance(10.);
5733 seed1->ResetCovariance(10.);
5734 FollowProlongation(*seed0,0);
5735 FollowBackProlongation(*seed1,158);
5736 mother = *seed0; // backup mother at position 0
5737 seed0->Reset(kFALSE);
5738 seed1->Reset(kFALSE);
5739 seed0->ResetCovariance(10.);
5740 seed1->ResetCovariance(10.);
5742 first = TMath::Max(row0-20,0);
5743 last = TMath::Min(row0+20,158);
5745 for (Int_t irow=0; irow<20;irow++){
5746 rows[irow] = first +((last-first)*irow)/19;
5748 // store parameters along the track
5750 for (Int_t irow=0;irow<20;irow++){
5751 FollowBackProlongation(*seed0, rows[irow]);
5752 FollowProlongation(*seed1,rows[19-irow]);
5753 param0[irow] = *seed0;
5754 param1[19-irow] = *seed1;
5758 for (Int_t irow=0; irow<19;irow++){
5759 kinks[irow].SetMother(param0[irow]);
5760 kinks[irow].SetDaughter(param1[irow]);
5761 // param0[irow].Dump();
5762 //param1[irow].Dump();
5763 kinks[irow].Update();
5766 // choose kink with biggest change of angle
5769 for (Int_t irow=0;irow<20;irow++){
5770 if (TMath::Abs(kinks[irow].GetR())>250.) continue;
5771 if (TMath::Abs(kinks[irow].GetR())<90.) continue;
5772 Float_t quality = TMath::Abs(kinks[irow].GetAngle(2))/(3.+TMath::Abs(kinks[irow].GetR()-param0[irow].GetX()));
5773 if ( quality > maxchange){
5774 maxchange = quality;
5781 if (index==-1 || param0[index].GetNumberOfClusters()+param1[index].GetNumberOfClusters()<100){
5787 // Float_t anglesigma = TMath::Sqrt(param0[index].fC22+param0[index].fC33+param1[index].fC22+param1[index].fC33);
5789 kink.SetMother(param0[index]);
5790 kink.SetDaughter(param1[index]);
5793 Double_t chi2P2 = param0[index].GetParameter()[2]-param1[index].GetParameter()[2];
5795 chi2P2/=param0[index].GetCovariance()[5]+param1[index].GetCovariance()[5];
5796 Double_t chi2P3 = param0[index].GetParameter()[3]-param1[index].GetParameter()[3];
5798 chi2P3/=param0[index].GetCovariance()[9]+param1[index].GetCovariance()[9];
5800 if (AliTPCReconstructor::StreamLevel()>1) {
5801 (*fDebugStreamer)<<"kinkHpt"<<
5804 "p0.="<<¶m0[index]<<
5805 "p1.="<<¶m1[index]<<
5809 if ( chi2P2+chi2P3<AliTPCReconstructor::GetRecoParam()->GetKinkAngleCutChi2(0)){
5816 row0 = GetRowNumber(kink.GetR());
5817 kink.SetTPCRow0(row0);
5818 kink.SetLabel(CookLabel(seed0,0.5,0,row0),0);
5819 kink.SetLabel(CookLabel(seed1,0.5,row0,158),1);
5820 kink.SetIndex(-10,0);
5821 kink.SetIndex(int(param0[index].GetNumberOfClusters()+param1[index].GetNumberOfClusters()),1);
5822 kink.SetTPCncls(param0[index].GetNumberOfClusters(),0);
5823 kink.SetTPCncls(param1[index].GetNumberOfClusters(),1);
5826 // new (&mother) AliTPCseed(param0[index]);
5827 daughter = param1[index];
5828 daughter.SetLabel(kink.GetLabel(1));
5829 param0[index].Reset(kTRUE);
5830 FollowProlongation(param0[index],0);
5831 mother = param0[index];
5832 mother.SetLabel(kink.GetLabel(0));
5833 if ( chi2P2+chi2P3<AliTPCReconstructor::GetRecoParam()->GetKinkAngleCutChi2(1)){
5845 AliTPCseed* AliTPCtrackerMI::ReSeed(AliTPCseed *t)
5848 // reseed - refit - track
5851 // Int_t last = fSectors->GetNRows()-1;
5853 if (fSectors == fOuterSec){
5854 first = TMath::Max(first, t->GetFirstPoint()-fInnerSec->GetNRows());
5858 first = t->GetFirstPoint();
5860 AliTPCseed * seed = MakeSeed(t,0.1,0.5,0.9);
5861 FollowBackProlongation(*t,fSectors->GetNRows()-1);
5863 FollowProlongation(*t,first);
5873 //_____________________________________________________________________________
5874 Int_t AliTPCtrackerMI::ReadSeeds(const TFile *inp) {
5875 //-----------------------------------------------------------------
5876 // This function reades track seeds.
5877 //-----------------------------------------------------------------
5878 TDirectory *savedir=gDirectory;
5880 TFile *in=(TFile*)inp;
5881 if (!in->IsOpen()) {
5882 cerr<<"AliTPCtrackerMI::ReadSeeds(): input file is not open !\n";
5887 TTree *seedTree=(TTree*)in->Get("Seeds");
5889 cerr<<"AliTPCtrackerMI::ReadSeeds(): ";
5890 cerr<<"can't get a tree with track seeds !\n";
5893 AliTPCtrack *seed=new AliTPCtrack;
5894 seedTree->SetBranchAddress("tracks",&seed);
5896 if (fSeeds==0) fSeeds=new TObjArray(15000);
5898 Int_t n=(Int_t)seedTree->GetEntries();
5899 for (Int_t i=0; i<n; i++) {
5900 seedTree->GetEvent(i);
5901 fSeeds->AddLast(new AliTPCseed(*seed/*,seed->GetAlpha()*/));
5910 Int_t AliTPCtrackerMI::Clusters2Tracks (AliESDEvent *const esd)
5913 // clusters to tracks
5915 if (fSeeds) DeleteSeeds();
5918 AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ;
5919 transform->SetCurrentTimeStamp( esd->GetTimeStamp());
5920 transform->SetCurrentRun(esd->GetRunNumber());
5923 if (!fSeeds) return 1;
5925 if (AliTPCReconstructor::StreamLevel()>3) DumpClusters(0,fSeeds);
5931 //_____________________________________________________________________________
5932 Int_t AliTPCtrackerMI::Clusters2Tracks() {
5933 //-----------------------------------------------------------------
5934 // This is a track finder.
5935 //-----------------------------------------------------------------
5936 TDirectory *savedir=gDirectory;
5940 fSeeds = Tracking();
5943 Info("Clusters2Tracks","Time for tracking: \t");timer.Print();timer.Start();
5945 //activate again some tracks
5946 for (Int_t i=0; i<fSeeds->GetEntriesFast(); i++) {
5947 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
5949 Int_t nc=t.GetNumberOfClusters();
5951 delete fSeeds->RemoveAt(i);
5955 if (pt->GetRemoval()==10) {
5956 if (pt->GetDensityFirst(20)>0.8 || pt->GetDensityFirst(30)>0.8 || pt->GetDensityFirst(40)>0.7)
5957 pt->Desactivate(10); // make track again active // MvL: should be 0 ?
5959 pt->Desactivate(20);
5960 delete fSeeds->RemoveAt(i);
5965 RemoveUsed2(fSeeds,0.85,0.85,0);
5966 if (AliTPCReconstructor::GetRecoParam()->GetDoKinks()) FindKinks(fSeeds,fEvent);
5967 //FindCurling(fSeeds, fEvent,0);
5968 if (AliTPCReconstructor::StreamLevel()>5) FindMultiMC(fSeeds, fEvent,-1); // find multi found tracks
5969 RemoveUsed2(fSeeds,0.5,0.4,20);
5970 FindSplitted(fSeeds, fEvent,0); // find multi found tracks
5971 if (AliTPCReconstructor::StreamLevel()>5) FindMultiMC(fSeeds, fEvent,0); // find multi found tracks
5974 // // refit short tracks
5976 Int_t nseed=fSeeds->GetEntriesFast();
5979 for (Int_t i=0; i<nseed; i++) {
5980 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
5982 Int_t nc=t.GetNumberOfClusters();
5984 delete fSeeds->RemoveAt(i);
5987 CookLabel(pt,0.1); //For comparison only
5988 //if ((pt->IsActive() || (pt->fRemoval==10) )&& nc>50 &&pt->GetNumberOfClusters()>0.4*pt->fNFoundable){
5989 if ((pt->IsActive() || (pt->GetRemoval()==10) )){
5991 if (fDebug>0) cerr<<found<<'\r';
5995 delete fSeeds->RemoveAt(i);
5999 //RemoveOverlap(fSeeds,0.99,7,kTRUE);
6001 //RemoveUsed(fSeeds,0.9,0.9,6);
6003 nseed=fSeeds->GetEntriesFast();
6005 for (Int_t i=0; i<nseed; i++) {
6006 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
6008 Int_t nc=t.GetNumberOfClusters();
6010 delete fSeeds->RemoveAt(i);
6014 t.CookdEdx(0.02,0.6);
6015 // CheckKinkPoint(&t,0.05);
6016 //if ((pt->IsActive() || (pt->fRemoval==10) )&& nc>50 &&pt->GetNumberOfClusters()>0.4*pt->fNFoundable){
6017 if ((pt->IsActive() || (pt->GetRemoval()==10) )){
6025 delete fSeeds->RemoveAt(i);
6026 //AliTPCseed * seed1 = ReSeed(pt,0.05,0.5,1);
6028 // FollowProlongation(*seed1,0);
6029 // Int_t n = seed1->GetNumberOfClusters();
6030 // printf("fP4\t%f\t%f\n",seed1->GetC(),pt->GetC());
6031 // printf("fN\t%d\t%d\n", seed1->GetNumberOfClusters(),pt->GetNumberOfClusters());
6034 //AliTPCseed * seed2 = ReSeed(pt,0.95,0.5,0.05);
6038 SortTracks(fSeeds, 1);
6042 PrepareForBackProlongation(fSeeds,5.);
6043 PropagateBack(fSeeds);
6044 printf("Time for back propagation: \t");timer.Print();timer.Start();
6048 PrepareForProlongation(fSeeds,5.);
6049 PropagateForard2(fSeeds);
6051 printf("Time for FORWARD propagation: \t");timer.Print();timer.Start();
6052 // RemoveUsed(fSeeds,0.7,0.7,6);
6053 //RemoveOverlap(fSeeds,0.9,7,kTRUE);
6055 nseed=fSeeds->GetEntriesFast();
6057 for (Int_t i=0; i<nseed; i++) {
6058 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
6060 Int_t nc=t.GetNumberOfClusters();
6062 delete fSeeds->RemoveAt(i);
6065 t.CookdEdx(0.02,0.6);
6066 // CookLabel(pt,0.1); //For comparison only
6067 //if ((pt->IsActive() || (pt->fRemoval==10) )&& nc>50 &&pt->GetNumberOfClusters()>0.4*pt->fNFoundable){
6068 if ((pt->IsActive() || (pt->fRemoval==10) )){
6069 cerr<<found++<<'\r';
6072 delete fSeeds->RemoveAt(i);
6077 // fNTracks = found;
6079 Info("Clusters2Tracks","Time for overlap removal, track writing and dedx cooking: \t"); timer.Print();timer.Start();
6082 // cerr<<"Number of found tracks : "<<"\t"<<found<<endl;
6083 Info("Clusters2Tracks","Number of found tracks %d",found);
6085 // UnloadClusters();
6090 void AliTPCtrackerMI::Tracking(TObjArray * arr)
6093 // tracking of the seeds
6096 fSectors = fOuterSec;
6097 ParallelTracking(arr,150,63);
6098 fSectors = fOuterSec;
6099 ParallelTracking(arr,63,0);
6102 TObjArray * AliTPCtrackerMI::Tracking(Int_t seedtype, Int_t i1, Int_t i2, Float_t cuts[4], Float_t dy, Int_t dsec)
6107 TObjArray * arr = new TObjArray;
6109 fSectors = fOuterSec;
6112 for (Int_t sec=0;sec<fkNOS;sec++){
6113 if (seedtype==3) MakeSeeds3(arr,sec,i1,i2,cuts,dy, dsec);
6114 if (seedtype==4) MakeSeeds5(arr,sec,i1,i2,cuts,dy);
6115 if (seedtype==2) MakeSeeds2(arr,sec,i1,i2,cuts,dy);
6118 Info("Tracking","\nSeeding - %d\t%d\t%d\t%d\n",seedtype,i1,i2,arr->GetEntriesFast());
6130 TObjArray * AliTPCtrackerMI::Tracking()
6134 if (AliTPCReconstructor::GetRecoParam()->GetSpecialSeeding()) return TrackingSpecial();
6137 Int_t nup=fOuterSec->GetNRows()+fInnerSec->GetNRows();
6139 TObjArray * seeds = new TObjArray;
6141 Int_t fLastSeedRowSec=AliTPCReconstructor::GetRecoParam()->GetLastSeedRowSec();
6142 Int_t gapPrim = AliTPCReconstructor::GetRecoParam()->GetSeedGapPrim();
6143 Int_t gapSec = AliTPCReconstructor::GetRecoParam()->GetSeedGapSec();
6151 Float_t fnumber = 3.0;
6152 Float_t fdensity = 3.0;
6157 for (Int_t delta = 0; delta<18; delta+=gapPrim){
6161 arr = Tracking(3,nup-1-delta,nup-1-delta-gap,cuts,-1,1);
6162 SumTracks(seeds,arr);
6163 SignClusters(seeds,fnumber,fdensity);
6165 for (Int_t i=2;i<6;i+=2){
6166 // seed high pt tracks
6169 arr = Tracking(3,nup-i-delta,nup-i-delta-gap,cuts,-1,0);
6170 SumTracks(seeds,arr);
6171 SignClusters(seeds,fnumber,fdensity);
6176 // RemoveUsed(seeds,0.9,0.9,1);
6177 // UnsignClusters();
6178 // SignClusters(seeds,fnumber,fdensity);
6182 for (Int_t delta = 20; delta<120; delta+=gapPrim){
6184 // seed high pt tracks
6188 arr = Tracking(3,nup-delta,nup-delta-gap,cuts,-1);
6189 SumTracks(seeds,arr);
6190 SignClusters(seeds,fnumber,fdensity);
6195 arr = Tracking(3,nup-delta-5,nup-delta-5-gap,cuts,-1);
6196 SumTracks(seeds,arr);
6197 SignClusters(seeds,fnumber,fdensity);
6208 Info("Tracking()","\n\nPrimary seeding\t%d\n\n",seeds->GetEntriesFast());
6212 // RemoveUsed(seeds,0.75,0.75,1);
6214 //SignClusters(seeds,fnumber,fdensity);
6223 arr = Tracking(4,nup-1,nup-1-gap,cuts,-1);
6224 SumTracks(seeds,arr);
6225 SignClusters(seeds,fnumber,fdensity);
6227 arr = Tracking(4,nup-2,nup-2-gap,cuts,-1);
6228 SumTracks(seeds,arr);
6229 SignClusters(seeds,fnumber,fdensity);
6231 arr = Tracking(4,nup-3,nup-3-gap,cuts,-1);
6232 SumTracks(seeds,arr);
6233 SignClusters(seeds,fnumber,fdensity);
6235 arr = Tracking(4,nup-5,nup-5-gap,cuts,-1);
6236 SumTracks(seeds,arr);
6237 SignClusters(seeds,fnumber,fdensity);
6239 arr = Tracking(4,nup-7,nup-7-gap,cuts,-1);
6240 SumTracks(seeds,arr);
6241 SignClusters(seeds,fnumber,fdensity);
6244 arr = Tracking(4,nup-9,nup-9-gap,cuts,-1);
6245 SumTracks(seeds,arr);
6246 SignClusters(seeds,fnumber,fdensity);
6250 for (Int_t delta = 9; delta<30; delta+=gapSec){
6256 arr = Tracking(4,nup-1-delta,nup-1-delta-gap,cuts,-1);
6257 SumTracks(seeds,arr);
6258 SignClusters(seeds,fnumber,fdensity);
6260 arr = Tracking(4,nup-3-delta,nup-5-delta-gap,cuts,4);
6261 SumTracks(seeds,arr);
6262 SignClusters(seeds,fnumber,fdensity);
6275 for (Int_t delta = 30; delta<fLastSeedRowSec; delta+=gapSec){
6281 arr = Tracking(4,nup-1-delta,nup-1-delta-gap,cuts,-1);
6282 SumTracks(seeds,arr);
6283 SignClusters(seeds,fnumber,fdensity);
6285 arr = Tracking(4,nup-5-delta,nup-5-delta-gap,cuts,5 );
6286 SumTracks(seeds,arr);
6287 SignClusters(seeds,fnumber,fdensity);
6291 Info("Tracking()","\n\nSecondary seeding\t%d\n\n",seeds->GetEntriesFast());
6302 TObjArray * AliTPCtrackerMI::TrackingSpecial()
6305 // seeding adjusted for laser and cosmic tests - short tracks with big inclination angle
6306 // no primary vertex seeding tried
6310 Int_t nup=fOuterSec->GetNRows()+fInnerSec->GetNRows();
6312 TObjArray * seeds = new TObjArray;
6317 Float_t fnumber = 3.0;
6318 Float_t fdensity = 3.0;
6321 cuts[0] = AliTPCReconstructor::GetRecoParam()->GetMaxC(); // max curvature
6322 cuts[1] = 3.5; // max tan(phi) angle for seeding
6323 cuts[2] = 3.; // not used (cut on z primary vertex)
6324 cuts[3] = 3.5; // max tan(theta) angle for seeding
6326 for (Int_t delta = 0; nup-delta-gap-1>0; delta+=3){
6328 arr = Tracking(4,nup-1-delta,nup-1-delta-gap,cuts,-1);
6329 SumTracks(seeds,arr);
6330 SignClusters(seeds,fnumber,fdensity);
6334 Info("Tracking()","\n\nSecondary seeding\t%d\n\n",seeds->GetEntriesFast());
6345 void AliTPCtrackerMI::SumTracks(TObjArray *arr1,TObjArray *&arr2) const
6348 //sum tracks to common container
6349 //remove suspicious tracks
6350 Int_t nseed = arr2->GetEntriesFast();
6351 for (Int_t i=0;i<nseed;i++){
6352 AliTPCseed *pt=(AliTPCseed*)arr2->UncheckedAt(i);
6355 // remove tracks with too big curvature
6357 if (TMath::Abs(pt->GetC())>AliTPCReconstructor::GetRecoParam()->GetMaxC()){
6358 delete arr2->RemoveAt(i);
6361 // REMOVE VERY SHORT TRACKS
6362 if (pt->GetNumberOfClusters()<20){
6363 delete arr2->RemoveAt(i);
6366 // NORMAL ACTIVE TRACK
6367 if (pt->IsActive()){
6368 arr1->AddLast(arr2->RemoveAt(i));
6371 //remove not usable tracks
6372 if (pt->GetRemoval()!=10){
6373 delete arr2->RemoveAt(i);
6377 // ENABLE ONLY ENOUGH GOOD STOPPED TRACKS
6378 if (pt->GetDensityFirst(20)>0.8 || pt->GetDensityFirst(30)>0.8 || pt->GetDensityFirst(40)>0.7)
6379 arr1->AddLast(arr2->RemoveAt(i));
6381 delete arr2->RemoveAt(i);
6385 delete arr2; arr2 = 0;
6390 void AliTPCtrackerMI::ParallelTracking(TObjArray *const arr, Int_t rfirst, Int_t rlast)
6393 // try to track in parralel
6395 Int_t nseed=arr->GetEntriesFast();
6396 //prepare seeds for tracking
6397 for (Int_t i=0; i<nseed; i++) {
6398 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i), &t=*pt;
6400 if (!t.IsActive()) continue;
6401 // follow prolongation to the first layer
6402 if ( (fSectors ==fInnerSec) || (t.GetFirstPoint()-fkParam->GetNRowLow()>rfirst+1) )
6403 FollowProlongation(t, rfirst+1);
6408 for (Int_t nr=rfirst; nr>=rlast; nr--){
6409 if (nr<fInnerSec->GetNRows())
6410 fSectors = fInnerSec;
6412 fSectors = fOuterSec;
6413 // make indexes with the cluster tracks for given
6415 // find nearest cluster
6416 for (Int_t i=0; i<nseed; i++) {
6417 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i), &t=*pt;
6419 if (nr==80) pt->UpdateReference();
6420 if (!pt->IsActive()) continue;
6421 // if ( (fSectors ==fOuterSec) && (pt->fFirstPoint-fkParam->GetNRowLow())<nr) continue;
6422 if (pt->GetRelativeSector()>17) {
6425 UpdateClusters(t,nr);
6427 // prolonagate to the nearest cluster - if founded
6428 for (Int_t i=0; i<nseed; i++) {
6429 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
6431 if (!pt->IsActive()) continue;
6432 // if ((fSectors ==fOuterSec) && (pt->fFirstPoint-fkParam->GetNRowLow())<nr) continue;
6433 if (pt->GetRelativeSector()>17) {
6436 FollowToNextCluster(*pt,nr);
6441 void AliTPCtrackerMI::PrepareForBackProlongation(const TObjArray *const arr,Float_t fac) const
6445 // if we use TPC track itself we have to "update" covariance
6447 Int_t nseed= arr->GetEntriesFast();
6448 for (Int_t i=0;i<nseed;i++){
6449 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
6453 //rotate to current local system at first accepted point
6454 Int_t index = pt->GetClusterIndex2(pt->GetFirstPoint());
6455 Int_t sec = (index&0xff000000)>>24;
6457 Float_t angle1 = fInnerSec->GetAlpha()*sec+fInnerSec->GetAlphaShift();
6458 if (angle1>TMath::Pi())
6459 angle1-=2.*TMath::Pi();
6460 Float_t angle2 = pt->GetAlpha();
6462 if (TMath::Abs(angle1-angle2)>0.001){
6463 if (!pt->Rotate(angle1-angle2)) return;
6464 //angle2 = pt->GetAlpha();
6465 //pt->fRelativeSector = pt->GetAlpha()/fInnerSec->GetAlpha();
6466 //if (pt->GetAlpha()<0)
6467 // pt->fRelativeSector+=18;
6468 //sec = pt->fRelativeSector;
6477 void AliTPCtrackerMI::PrepareForProlongation(TObjArray *const arr, Float_t fac) const
6481 // if we use TPC track itself we have to "update" covariance
6483 Int_t nseed= arr->GetEntriesFast();
6484 for (Int_t i=0;i<nseed;i++){
6485 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
6488 pt->SetFirstPoint(pt->GetLastPoint());
6496 Int_t AliTPCtrackerMI::PropagateBack(const TObjArray *const arr)
6499 // make back propagation
6501 Int_t nseed= arr->GetEntriesFast();
6502 for (Int_t i=0;i<nseed;i++){
6503 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
6504 if (pt&& pt->GetKinkIndex(0)<=0) {
6505 //AliTPCseed *pt2 = new AliTPCseed(*pt);
6506 fSectors = fInnerSec;
6507 //FollowBackProlongation(*pt,fInnerSec->GetNRows()-1);
6508 //fSectors = fOuterSec;
6509 FollowBackProlongation(*pt,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1,1);
6510 //if (pt->GetNumberOfClusters()<(pt->fEsd->GetTPCclusters(0)) ){
6511 // Error("PropagateBack","Not prolonged track %d",pt->GetLabel());
6512 // FollowBackProlongation(*pt2,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1);
6515 if (pt&& pt->GetKinkIndex(0)>0) {
6516 AliESDkink * kink = fEvent->GetKink(pt->GetKinkIndex(0)-1);
6517 pt->SetFirstPoint(kink->GetTPCRow0());
6518 fSectors = fInnerSec;
6519 FollowBackProlongation(*pt,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1,1);
6527 Int_t AliTPCtrackerMI::PropagateForward2(const TObjArray *const arr)
6530 // make forward propagation
6532 Int_t nseed= arr->GetEntriesFast();
6534 for (Int_t i=0;i<nseed;i++){
6535 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
6537 FollowProlongation(*pt,0,1,1);
6546 Int_t AliTPCtrackerMI::PropagateForward()
6549 // propagate track forward
6551 Int_t nseed = fSeeds->GetEntriesFast();
6552 for (Int_t i=0;i<nseed;i++){
6553 AliTPCseed *pt = (AliTPCseed*)fSeeds->UncheckedAt(i);
6555 AliTPCseed &t = *pt;
6556 Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
6557 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
6558 if (alpha < 0. ) alpha += 2.*TMath::Pi();
6559 t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN);
6563 fSectors = fOuterSec;
6564 ParallelTracking(fSeeds,fOuterSec->GetNRows()+fInnerSec->GetNRows()-1,fInnerSec->GetNRows());
6565 fSectors = fInnerSec;
6566 ParallelTracking(fSeeds,fInnerSec->GetNRows()-1,0);
6575 Int_t AliTPCtrackerMI::PropagateBack(AliTPCseed *const pt, Int_t row0, Int_t row1)
6578 // make back propagation, in between row0 and row1
6582 fSectors = fInnerSec;
6585 if (row1<fSectors->GetNRows())
6588 r1 = fSectors->GetNRows()-1;
6590 if (row0<fSectors->GetNRows()&& r1>0 )
6591 FollowBackProlongation(*pt,r1);
6592 if (row1<=fSectors->GetNRows())
6595 r1 = row1 - fSectors->GetNRows();
6596 if (r1<=0) return 0;
6597 if (r1>=fOuterSec->GetNRows()) return 0;
6598 fSectors = fOuterSec;
6599 return FollowBackProlongation(*pt,r1);
6607 void AliTPCtrackerMI::GetShape(AliTPCseed * seed, Int_t row)
6609 // gets cluster shape
6611 AliTPCClusterParam * clparam = AliTPCcalibDB::Instance()->GetClusterParam();
6612 Float_t zdrift = TMath::Abs((fkParam->GetZLength(0)-TMath::Abs(seed->GetZ())));
6613 Int_t type = (seed->GetSector() < fkParam->GetNSector()/2) ? 0: (row>126) ? 1:2;
6614 Double_t angulary = seed->GetSnp();
6616 if (TMath::Abs(angulary)>AliTPCReconstructor::GetMaxSnpTracker()) {
6617 angulary = TMath::Sign(AliTPCReconstructor::GetMaxSnpTracker(),angulary);
6620 angulary = angulary*angulary/((1.-angulary)*(1.+angulary));
6621 Double_t angularz = seed->GetTgl()*seed->GetTgl()*(1.+angulary);
6623 Double_t sigmay = clparam->GetRMS0(0,type,zdrift,TMath::Sqrt(TMath::Abs(angulary)));
6624 Double_t sigmaz = clparam->GetRMS0(1,type,zdrift,TMath::Sqrt(TMath::Abs(angularz)));
6625 seed->SetCurrentSigmaY2(sigmay*sigmay);
6626 seed->SetCurrentSigmaZ2(sigmaz*sigmaz);
6627 // Float_t sd2 = TMath::Abs((fkParam->GetZLength(0)-TMath::Abs(seed->GetZ())))*fkParam->GetDiffL()*fkParam->GetDiffL();
6628 // // Float_t padlength = fkParam->GetPadPitchLength(seed->fSector);
6629 // Float_t padlength = GetPadPitchLength(row);
6631 // Float_t sresy = (seed->GetSector() < fkParam->GetNSector()/2) ? 0.2 :0.3;
6632 // seed->SetCurrentSigmaY2(sd2+padlength*padlength*angulary/12.+sresy*sresy);
6634 // Float_t sresz = fkParam->GetZSigma();
6635 // seed->SetCurrentSigmaZ2(sd2+padlength*padlength*angularz*angularz*(1+angulary)/12.+sresz*sresz);
6637 Float_t wy = GetSigmaY(seed);
6638 Float_t wz = GetSigmaZ(seed);
6641 if (TMath::Abs(wy/seed->fCurrentSigmaY2-1)>0.0001 || TMath::Abs(wz/seed->fCurrentSigmaZ2-1)>0.0001 ){
6642 printf("problem\n");
6649 //__________________________________________________________________________
6650 void AliTPCtrackerMI::CookLabel(AliKalmanTrack *tk, Float_t wrong) const {
6651 //--------------------------------------------------------------------
6652 //This function "cooks" a track label. If label<0, this track is fake.
6653 //--------------------------------------------------------------------
6654 AliTPCseed * t = dynamic_cast<AliTPCseed*>(tk);
6656 printf("%s:%d wrong type \n",(char*)__FILE__,__LINE__);
6660 Int_t noc=t->GetNumberOfClusters();
6662 //printf("\nnot founded prolongation\n\n\n");
6668 AliTPCclusterMI *clusters[160];
6670 for (Int_t i=0;i<160;i++) {
6677 for (i=0; i<160 && current<noc; i++) {
6679 Int_t index=t->GetClusterIndex2(i);
6680 if (index<=0) continue;
6681 if (index&0x8000) continue;
6683 //clusters[current]=GetClusterMI(index);
6684 if (t->GetClusterPointer(i)){
6685 clusters[current]=t->GetClusterPointer(i);
6691 Int_t lab=123456789;
6692 for (i=0; i<noc; i++) {
6693 AliTPCclusterMI *c=clusters[i];
6695 lab=TMath::Abs(c->GetLabel(0));
6697 for (j=0; j<noc; j++) if (lb[j]==lab || mx[j]==0) break;
6703 for (i=0; i<noc; i++) if (mx[i]>max) {max=mx[i]; lab=lb[i];}
6705 for (i=0; i<noc; i++) {
6706 AliTPCclusterMI *c=clusters[i];
6708 if (TMath::Abs(c->GetLabel(1)) == lab ||
6709 TMath::Abs(c->GetLabel(2)) == lab ) max++;
6711 if (noc<=0) { lab=-1; return;}
6712 if ((1.- Float_t(max)/(noc)) > wrong) lab=-lab;
6715 Int_t tail=Int_t(0.10*noc);
6718 for (i=1; i<160&&ind<tail; i++) {
6719 // AliTPCclusterMI *c=clusters[noc-i];
6720 AliTPCclusterMI *c=clusters[i];
6722 if (lab == TMath::Abs(c->GetLabel(0)) ||
6723 lab == TMath::Abs(c->GetLabel(1)) ||
6724 lab == TMath::Abs(c->GetLabel(2))) max++;
6727 if (max < Int_t(0.5*tail)) lab=-lab;
6734 //delete[] clusters;
6738 //__________________________________________________________________________
6739 Int_t AliTPCtrackerMI::CookLabel(AliTPCseed *const t, Float_t wrong,Int_t first, Int_t last) const {
6740 //--------------------------------------------------------------------
6741 //This function "cooks" a track label. If label<0, this track is fake.
6742 //--------------------------------------------------------------------
6743 Int_t noc=t->GetNumberOfClusters();
6745 //printf("\nnot founded prolongation\n\n\n");
6751 AliTPCclusterMI *clusters[160];
6753 for (Int_t i=0;i<160;i++) {
6760 for (i=0; i<160 && current<noc; i++) {
6761 if (i<first) continue;
6762 if (i>last) continue;
6763 Int_t index=t->GetClusterIndex2(i);
6764 if (index<=0) continue;
6765 if (index&0x8000) continue;
6767 //clusters[current]=GetClusterMI(index);
6768 if (t->GetClusterPointer(i)){
6769 clusters[current]=t->GetClusterPointer(i);
6774 //if (noc<5) return -1;
6775 Int_t lab=123456789;
6776 for (i=0; i<noc; i++) {
6777 AliTPCclusterMI *c=clusters[i];
6779 lab=TMath::Abs(c->GetLabel(0));
6781 for (j=0; j<noc; j++) if (lb[j]==lab || mx[j]==0) break;
6787 for (i=0; i<noc; i++) if (mx[i]>max) {max=mx[i]; lab=lb[i];}
6789 for (i=0; i<noc; i++) {
6790 AliTPCclusterMI *c=clusters[i];
6792 if (TMath::Abs(c->GetLabel(1)) == lab ||
6793 TMath::Abs(c->GetLabel(2)) == lab ) max++;
6795 if (noc<=0) { lab=-1; return -1;}
6796 if ((1.- Float_t(max)/(noc)) > wrong) lab=-lab;
6799 Int_t tail=Int_t(0.10*noc);
6802 for (i=1; i<160&&ind<tail; i++) {
6803 // AliTPCclusterMI *c=clusters[noc-i];
6804 AliTPCclusterMI *c=clusters[i];
6806 if (lab == TMath::Abs(c->GetLabel(0)) ||
6807 lab == TMath::Abs(c->GetLabel(1)) ||
6808 lab == TMath::Abs(c->GetLabel(2))) max++;
6811 if (max < Int_t(0.5*tail)) lab=-lab;
6814 // t->SetLabel(lab);
6818 //delete[] clusters;
6822 Int_t AliTPCtrackerMI::GetRowNumber(Double_t x[3]) const
6824 //return pad row number for given x vector
6825 Float_t phi = TMath::ATan2(x[1],x[0]);
6826 if(phi<0) phi=2.*TMath::Pi()+phi;
6827 // Get the local angle in the sector philoc
6828 const Float_t kRaddeg = 180/3.14159265358979312;
6829 Float_t phiangle = (Int_t (phi*kRaddeg/20.) + 0.5)*20./kRaddeg;
6830 Double_t localx = x[0]*TMath::Cos(phiangle)-x[1]*TMath::Sin(phiangle);
6831 return GetRowNumber(localx);
6836 void AliTPCtrackerMI::MakeESDBitmaps(AliTPCseed *t, AliESDtrack *esd)
6838 //-----------------------------------------------------------------------
6839 // Fill the cluster and sharing bitmaps of the track
6840 //-----------------------------------------------------------------------
6842 Int_t firstpoint = 0;
6843 Int_t lastpoint = 159;
6844 AliTPCTrackerPoint *point;
6845 AliTPCclusterMI *cluster;
6848 TBits clusterMap(159);
6849 TBits sharedMap(159);
6851 for (int iter=firstpoint; iter<lastpoint; iter++) {
6852 // Change to cluster pointers to see if we have a cluster at given padrow
6854 cluster = t->GetClusterPointer(iter);
6856 clusterMap.SetBitNumber(iter, kTRUE);
6857 point = t->GetTrackPoint(iter);
6858 if (point->IsShared())
6859 sharedMap.SetBitNumber(iter,kTRUE);
6861 if (t->GetClusterIndex(iter) >= 0 && (t->GetClusterIndex(iter) & 0x8000) == 0) {
6862 fitMap.SetBitNumber(iter, kTRUE);
6866 esd->SetTPCClusterMap(clusterMap);
6867 esd->SetTPCSharedMap(sharedMap);
6868 esd->SetTPCFitMap(fitMap);
6869 if (nclsf != t->GetNumberOfClusters())
6870 AliWarning(Form("Inconsistency between ncls %d and indices %d (found %d)",t->GetNumberOfClusters(),nclsf,esd->GetTPCClusterMap().CountBits()));
6873 Bool_t AliTPCtrackerMI::IsFindable(AliTPCseed & track){
6875 // return flag if there is findable cluster at given position
6878 Float_t z = track.GetZ();
6880 if (TMath::Abs(z)<(AliTPCReconstructor::GetCtgRange()*track.GetX()+kDeltaZ) &&
6881 TMath::Abs(z)<fkParam->GetZLength(0) &&
6882 (TMath::Abs(track.GetSnp())<AliTPCReconstructor::GetMaxSnpTracker()))
6888 void AliTPCtrackerMI::AddCovariance(AliTPCseed * seed){
6890 // Adding systematic error estimate to the covariance matrix
6891 // !!!! the systematic error for element 4 is in 1/cm not in pt
6893 const Double_t *param = AliTPCReconstructor::GetRecoParam()->GetSystematicError();
6895 // use only the diagonal part if not specified otherwise
6896 if (!AliTPCReconstructor::GetRecoParam()->GetUseSystematicCorrelation()) return AddCovarianceAdd(seed);
6898 Double_t *covarS= (Double_t*)seed->GetCovariance();
6899 Double_t factor[5]={1,1,1,1,1};
6900 Double_t facC = AliTracker::GetBz()*kB2C;
6901 factor[0]= TMath::Sqrt(TMath::Abs((covarS[0] + param[0]*param[0])/covarS[0]));
6902 factor[1]= TMath::Sqrt(TMath::Abs((covarS[2] + param[1]*param[1])/covarS[2]));
6903 factor[2]= TMath::Sqrt(TMath::Abs((covarS[5] + param[2]*param[2])/covarS[5]));
6904 factor[3]= TMath::Sqrt(TMath::Abs((covarS[9] + param[3]*param[3])/covarS[9]));
6905 factor[4]= TMath::Sqrt(TMath::Abs((covarS[14] + facC*facC*param[4]*param[4])/covarS[14]));
6911 for (Int_t i=0; i<5; i++){
6912 for (Int_t j=i; j<5; j++){
6913 Int_t index=seed->GetIndex(i,j);
6914 covarS[index]*=factor[i]*factor[j];
6920 void AliTPCtrackerMI::AddCovarianceAdd(AliTPCseed * seed){
6922 // Adding systematic error - as additive factor without correlation
6924 const Double_t *param = AliTPCReconstructor::GetRecoParam()->GetSystematicError();
6925 Double_t *covarIn= (Double_t*)seed->GetCovariance();
6927 for (Int_t i=0;i<15;i++) covar[i]=0;
6933 covar[0] = param[0]*param[0];
6934 covar[2] = param[1]*param[1];
6935 covar[5] = param[2]*param[2];
6936 covar[9] = param[3]*param[3];
6937 Double_t facC = AliTracker::GetBz()*kB2C;
6938 covar[14]= param[4]*param[4]*facC*facC;
6940 covar[1]=TMath::Sqrt((covar[0]*covar[2]))*covarIn[1]/TMath::Sqrt((covarIn[0]*covarIn[2]));
6942 covar[3]=TMath::Sqrt((covar[0]*covar[5]))*covarIn[3]/TMath::Sqrt((covarIn[0]*covarIn[5]));
6943 covar[4]=TMath::Sqrt((covar[2]*covar[5]))*covarIn[4]/TMath::Sqrt((covarIn[2]*covarIn[5]));
6945 covar[6]=TMath::Sqrt((covar[0]*covar[9]))*covarIn[6]/TMath::Sqrt((covarIn[0]*covarIn[9]));
6946 covar[7]=TMath::Sqrt((covar[2]*covar[9]))*covarIn[7]/TMath::Sqrt((covarIn[2]*covarIn[9]));
6947 covar[8]=TMath::Sqrt((covar[5]*covar[9]))*covarIn[8]/TMath::Sqrt((covarIn[5]*covarIn[9]));
6949 covar[10]=TMath::Sqrt((covar[0]*covar[14]))*covarIn[10]/TMath::Sqrt((covarIn[0]*covarIn[14]));
6950 covar[11]=TMath::Sqrt((covar[2]*covar[14]))*covarIn[11]/TMath::Sqrt((covarIn[2]*covarIn[14]));
6951 covar[12]=TMath::Sqrt((covar[5]*covar[14]))*covarIn[12]/TMath::Sqrt((covarIn[5]*covarIn[14]));
6952 covar[13]=TMath::Sqrt((covar[9]*covar[14]))*covarIn[13]/TMath::Sqrt((covarIn[9]*covarIn[14]));
6954 seed->AddCovariance(covar);