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()
206 // default constructor
208 for (Int_t irow=0; irow<200; irow++){
215 //_____________________________________________________________________
219 Int_t AliTPCtrackerMI::UpdateTrack(AliTPCseed * track, Int_t accept){
221 //update track information using current cluster - track->fCurrentCluster
224 AliTPCclusterMI* c =track->GetCurrentCluster();
225 if (accept > 0) //sign not accepted clusters
226 track->SetCurrentClusterIndex1(track->GetCurrentClusterIndex1() | 0x8000);
227 else // unsign accpeted clusters
228 track->SetCurrentClusterIndex1(track->GetCurrentClusterIndex1() & 0xffff7fff);
229 UInt_t i = track->GetCurrentClusterIndex1();
231 Int_t sec=(i&0xff000000)>>24;
232 //Int_t row = (i&0x00ff0000)>>16;
233 track->SetRow((i&0x00ff0000)>>16);
234 track->SetSector(sec);
235 // Int_t index = i&0xFFFF;
236 if (sec>=fkParam->GetNInnerSector()) track->SetRow(track->GetRow()+fkParam->GetNRowLow());
237 track->SetClusterIndex2(track->GetRow(), i);
238 //track->fFirstPoint = row;
239 //if ( track->fLastPoint<row) track->fLastPoint =row;
240 // if (track->fRow<0 || track->fRow>160) {
241 // printf("problem\n");
243 if (track->GetFirstPoint()>track->GetRow())
244 track->SetFirstPoint(track->GetRow());
245 if (track->GetLastPoint()<track->GetRow())
246 track->SetLastPoint(track->GetRow());
249 track->SetClusterPointer(track->GetRow(),c);
252 Double_t angle2 = track->GetSnp()*track->GetSnp();
254 //SET NEW Track Point
256 if (angle2<1) //PH sometimes angle2 is very big. To be investigated...
258 angle2 = TMath::Sqrt(angle2/(1-angle2));
259 AliTPCTrackerPoint &point =*(track->GetTrackPoint(track->GetRow()));
261 point.SetSigmaY(c->GetSigmaY2()/track->GetCurrentSigmaY2());
262 point.SetSigmaZ(c->GetSigmaZ2()/track->GetCurrentSigmaZ2());
263 point.SetErrY(sqrt(track->GetErrorY2()));
264 point.SetErrZ(sqrt(track->GetErrorZ2()));
266 point.SetX(track->GetX());
267 point.SetY(track->GetY());
268 point.SetZ(track->GetZ());
269 point.SetAngleY(angle2);
270 point.SetAngleZ(track->GetTgl());
271 if (point.IsShared()){
272 track->SetErrorY2(track->GetErrorY2()*4);
273 track->SetErrorZ2(track->GetErrorZ2()*4);
277 Double_t chi2 = track->GetPredictedChi2(track->GetCurrentCluster());
279 // track->SetErrorY2(track->GetErrorY2()*1.3);
280 // track->SetErrorY2(track->GetErrorY2()+0.01);
281 // track->SetErrorZ2(track->GetErrorZ2()*1.3);
282 // track->SetErrorZ2(track->GetErrorZ2()+0.005);
284 if (accept>0) return 0;
285 if (track->GetNumberOfClusters()%20==0){
286 // if (track->fHelixIn){
287 // TClonesArray & larr = *(track->fHelixIn);
288 // Int_t ihelix = larr.GetEntriesFast();
289 // new(larr[ihelix]) AliHelix(*track) ;
292 track->SetNoCluster(0);
293 return track->Update(c,chi2,i);
298 Int_t AliTPCtrackerMI::AcceptCluster(AliTPCseed * seed, AliTPCclusterMI * cluster)
301 // decide according desired precision to accept given
302 // cluster for tracking
304 seed->GetProlongation(cluster->GetX(),yt,zt);
305 Double_t sy2=ErrY2(seed,cluster);
306 Double_t sz2=ErrZ2(seed,cluster);
308 Double_t sdistancey2 = sy2+seed->GetSigmaY2();
309 Double_t sdistancez2 = sz2+seed->GetSigmaZ2();
310 Double_t dy=seed->GetCurrentCluster()->GetY()-yt;
311 Double_t dz=seed->GetCurrentCluster()->GetZ()-zt;
312 Double_t rdistancey2 = (seed->GetCurrentCluster()->GetY()-yt)*
313 (seed->GetCurrentCluster()->GetY()-yt)/sdistancey2;
314 Double_t rdistancez2 = (seed->GetCurrentCluster()->GetZ()-zt)*
315 (seed->GetCurrentCluster()->GetZ()-zt)/sdistancez2;
317 Double_t rdistance2 = rdistancey2+rdistancez2;
320 if (AliTPCReconstructor::StreamLevel()>2 && seed->GetNumberOfClusters()>20) {
321 Float_t rmsy2 = seed->GetCurrentSigmaY2();
322 Float_t rmsz2 = seed->GetCurrentSigmaZ2();
323 Float_t rmsy2p30 = seed->GetCMeanSigmaY2p30();
324 Float_t rmsz2p30 = seed->GetCMeanSigmaZ2p30();
325 Float_t rmsy2p30R = seed->GetCMeanSigmaY2p30R();
326 Float_t rmsz2p30R = seed->GetCMeanSigmaZ2p30R();
327 AliExternalTrackParam param(*seed);
328 static TVectorD gcl(3),gtr(3);
330 param.GetXYZ(gcl.GetMatrixArray());
331 cluster->GetGlobalXYZ(gclf);
332 gcl[0]=gclf[0]; gcl[1]=gclf[1]; gcl[2]=gclf[2];
335 if (AliTPCReconstructor::StreamLevel()>2) {
336 (*fDebugStreamer)<<"ErrParam"<<
349 "rmsy2p30="<<rmsy2p30<<
350 "rmsz2p30="<<rmsz2p30<<
351 "rmsy2p30R="<<rmsy2p30R<<
352 "rmsz2p30R="<<rmsz2p30R<<
353 // normalize distance -
354 "rdisty="<<rdistancey2<<
355 "rdistz="<<rdistancez2<<
356 "rdist="<<rdistance2<< //
360 //return 0; // temporary
361 if (rdistance2>32) return 3;
364 if ((rdistancey2>9. || rdistancez2>9.) && cluster->GetType()==0)
365 return 2; //suspisiouce - will be changed
367 if ((rdistancey2>6.25 || rdistancez2>6.25) && cluster->GetType()>0)
368 // strict cut on overlaped cluster
369 return 2; //suspisiouce - will be changed
371 if ( (rdistancey2>1. || rdistancez2>6.25 )
372 && cluster->GetType()<0){
373 seed->SetNFoundable(seed->GetNFoundable()-1);
377 if (AliTPCReconstructor::GetRecoParam()->GetUseHLTClusters() == 3 ||
378 AliTPCReconstructor::GetRecoParam()->GetUseHLTClusters() == 4 ) {
379 if(!AliTPCReconstructor::GetRecoParam()->GetUseOnePadCluster())
380 if (TMath::Abs(cluster->GetSigmaY2()) < kAlmost0)
390 //_____________________________________________________________________________
391 AliTPCtrackerMI::AliTPCtrackerMI(const AliTPCParam *par):
393 fkNIS(par->GetNInnerSector()/2),
395 fkNOS(par->GetNOuterSector()/2),
412 //---------------------------------------------------------------------
413 // The main TPC tracker constructor
414 //---------------------------------------------------------------------
415 fInnerSec=new AliTPCtrackerSector[fkNIS];
416 fOuterSec=new AliTPCtrackerSector[fkNOS];
419 for (i=0; i<fkNIS; i++) fInnerSec[i].Setup(par,0);
420 for (i=0; i<fkNOS; i++) fOuterSec[i].Setup(par,1);
423 Int_t nrowlow = par->GetNRowLow();
424 Int_t nrowup = par->GetNRowUp();
427 for (i=0;i<nrowlow;i++){
428 fXRow[i] = par->GetPadRowRadiiLow(i);
429 fPadLength[i]= par->GetPadPitchLength(0,i);
430 fYMax[i] = fXRow[i]*TMath::Tan(0.5*par->GetInnerAngle());
434 for (i=0;i<nrowup;i++){
435 fXRow[i+nrowlow] = par->GetPadRowRadiiUp(i);
436 fPadLength[i+nrowlow] = par->GetPadPitchLength(60,i);
437 fYMax[i+nrowlow] = fXRow[i+nrowlow]*TMath::Tan(0.5*par->GetOuterAngle());
440 if (AliTPCReconstructor::StreamLevel()>0) {
441 fDebugStreamer = new TTreeSRedirector("TPCdebug.root");
444 //________________________________________________________________________
445 AliTPCtrackerMI::AliTPCtrackerMI(const AliTPCtrackerMI &t):
466 //------------------------------------
467 // dummy copy constructor
468 //------------------------------------------------------------------
470 for (Int_t irow=0; irow<200; irow++){
477 AliTPCtrackerMI & AliTPCtrackerMI::operator=(const AliTPCtrackerMI& /*r*/)
479 //------------------------------
481 //--------------------------------------------------------------
484 //_____________________________________________________________________________
485 AliTPCtrackerMI::~AliTPCtrackerMI() {
486 //------------------------------------------------------------------
487 // TPC tracker destructor
488 //------------------------------------------------------------------
495 if (fDebugStreamer) delete fDebugStreamer;
499 void AliTPCtrackerMI::FillESD(const TObjArray* arr)
503 //fill esds using updated tracks
506 // write tracks to the event
507 // store index of the track
508 Int_t nseed=arr->GetEntriesFast();
509 //FindKinks(arr,fEvent);
510 for (Int_t i=0; i<nseed; i++) {
511 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
515 if (AliTPCReconstructor::StreamLevel()>1) {
516 (*fDebugStreamer)<<"Track0"<<
520 // pt->PropagateTo(fkParam->GetInnerRadiusLow());
521 if (pt->GetKinkIndex(0)<=0){ //don't propagate daughter tracks
522 pt->PropagateTo(fkParam->GetInnerRadiusLow());
525 if (( pt->GetPoints()[2]- pt->GetPoints()[0])>5 && pt->GetPoints()[3]>0.8){
527 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
528 iotrack.SetTPCPoints(pt->GetPoints());
529 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
530 iotrack.SetV0Indexes(pt->GetV0Indexes());
531 // iotrack.SetTPCpid(pt->fTPCr);
532 //iotrack.SetTPCindex(i);
533 fEvent->AddTrack(&iotrack);
537 if ( (pt->GetNumberOfClusters()>70)&& (Float_t(pt->GetNumberOfClusters())/Float_t(pt->GetNFoundable()))>0.55) {
539 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
540 iotrack.SetTPCPoints(pt->GetPoints());
541 //iotrack.SetTPCindex(i);
542 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
543 iotrack.SetV0Indexes(pt->GetV0Indexes());
544 // iotrack.SetTPCpid(pt->fTPCr);
545 fEvent->AddTrack(&iotrack);
549 // short tracks - maybe decays
551 if ( (pt->GetNumberOfClusters()>30) && (Float_t(pt->GetNumberOfClusters())/Float_t(pt->GetNFoundable()))>0.70) {
552 Int_t found,foundable,shared;
553 pt->GetClusterStatistic(0,60,found, foundable,shared,kFALSE);
554 if ( (found>20) && (pt->GetNShared()/float(pt->GetNumberOfClusters())<0.2)){
556 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
557 //iotrack.SetTPCindex(i);
558 iotrack.SetTPCPoints(pt->GetPoints());
559 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
560 iotrack.SetV0Indexes(pt->GetV0Indexes());
561 //iotrack.SetTPCpid(pt->fTPCr);
562 fEvent->AddTrack(&iotrack);
567 if ( (pt->GetNumberOfClusters()>20) && (Float_t(pt->GetNumberOfClusters())/Float_t(pt->GetNFoundable()))>0.8) {
568 Int_t found,foundable,shared;
569 pt->GetClusterStatistic(0,60,found, foundable,shared,kFALSE);
570 if (found<20) continue;
571 if (pt->GetNShared()/float(pt->GetNumberOfClusters())>0.2) continue;
574 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
575 iotrack.SetTPCPoints(pt->GetPoints());
576 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
577 iotrack.SetV0Indexes(pt->GetV0Indexes());
578 //iotrack.SetTPCpid(pt->fTPCr);
579 //iotrack.SetTPCindex(i);
580 fEvent->AddTrack(&iotrack);
583 // short tracks - secondaties
585 if ( (pt->GetNumberOfClusters()>30) ) {
586 Int_t found,foundable,shared;
587 pt->GetClusterStatistic(128,158,found, foundable,shared,kFALSE);
588 if ( (found>20) && (pt->GetNShared()/float(pt->GetNumberOfClusters())<0.2) &&float(found)/float(foundable)>0.8){
590 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
591 iotrack.SetTPCPoints(pt->GetPoints());
592 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
593 iotrack.SetV0Indexes(pt->GetV0Indexes());
594 //iotrack.SetTPCpid(pt->fTPCr);
595 //iotrack.SetTPCindex(i);
596 fEvent->AddTrack(&iotrack);
601 if ( (pt->GetNumberOfClusters()>15)) {
602 Int_t found,foundable,shared;
603 pt->GetClusterStatistic(138,158,found, foundable,shared,kFALSE);
604 if (found<15) continue;
605 if (foundable<=0) continue;
606 if (pt->GetNShared()/float(pt->GetNumberOfClusters())>0.2) continue;
607 if (float(found)/float(foundable)<0.8) continue;
610 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
611 iotrack.SetTPCPoints(pt->GetPoints());
612 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
613 iotrack.SetV0Indexes(pt->GetV0Indexes());
614 // iotrack.SetTPCpid(pt->fTPCr);
615 //iotrack.SetTPCindex(i);
616 fEvent->AddTrack(&iotrack);
620 // >> account for suppressed tracks in the kink indices (RS)
621 int nESDtracks = fEvent->GetNumberOfTracks();
622 for (int it=nESDtracks;it--;) {
623 AliESDtrack* esdTr = fEvent->GetTrack(it);
624 if (!esdTr || !esdTr->GetKinkIndex(0)) continue;
625 for (int ik=0;ik<3;ik++) {
627 if (!(knkId=esdTr->GetKinkIndex(ik))) break; // no more kinks for this track
628 AliESDkink* kink = fEvent->GetKink(TMath::Abs(knkId)-1);
630 AliError(Form("ESDTrack%d refers to non-existing kink %d",it,TMath::Abs(knkId)-1));
633 kink->SetIndex(it, knkId<0 ? 0:1); // update track index of the kink: mother at 0, daughter at 1
636 // << account for suppressed tracks in the kink indices (RS)
637 AliInfo(Form("Number of filled ESDs-\t%d\n",fEvent->GetNumberOfTracks()));
645 Double_t AliTPCtrackerMI::ErrY2(AliTPCseed* seed, const AliTPCclusterMI * cl){
648 // Use calibrated cluster error from OCDB
650 AliTPCClusterParam * clparam = AliTPCcalibDB::Instance()->GetClusterParam();
652 Float_t z = TMath::Abs(fkParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
653 Int_t ctype = cl->GetType();
654 Int_t type = (cl->GetRow()<63) ? 0: (cl->GetRow()>126) ? 1:2;
655 Double_t angle = seed->GetSnp()*seed->GetSnp();
656 angle = TMath::Sqrt(TMath::Abs(angle/(1.-angle)));
657 Double_t erry2 = clparam->GetError0Par(0,type, z,angle);
659 erry2+=0.5; // edge cluster
662 seed->SetErrorY2(erry2);
666 //calculate look-up table at the beginning
667 // static Bool_t ginit = kFALSE;
668 // static Float_t gnoise1,gnoise2,gnoise3;
669 // static Float_t ggg1[10000];
670 // static Float_t ggg2[10000];
671 // static Float_t ggg3[10000];
672 // static Float_t glandau1[10000];
673 // static Float_t glandau2[10000];
674 // static Float_t glandau3[10000];
676 // static Float_t gcor01[500];
677 // static Float_t gcor02[500];
678 // static Float_t gcorp[500];
682 // if (ginit==kFALSE){
683 // for (Int_t i=1;i<500;i++){
684 // Float_t rsigma = float(i)/100.;
685 // gcor02[i] = TMath::Max(0.78 +TMath::Exp(7.4*(rsigma-1.2)),0.6);
686 // gcor01[i] = TMath::Max(0.72 +TMath::Exp(3.36*(rsigma-1.2)),0.6);
687 // gcorp[i] = TMath::Max(TMath::Power((rsigma+0.5),1.5),1.2);
691 // for (Int_t i=3;i<10000;i++){
695 // Float_t amp = float(i);
696 // Float_t padlength =0.75;
697 // gnoise1 = 0.0004/padlength;
698 // Float_t nel = 0.268*amp;
699 // Float_t nprim = 0.155*amp;
700 // ggg1[i] = fkParam->GetDiffT()*fkParam->GetDiffT()*(2+0.001*nel/(padlength*padlength))/nel;
701 // glandau1[i] = (2.+0.12*nprim)*0.5* (2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
702 // if (glandau1[i]>1) glandau1[i]=1;
703 // glandau1[i]*=padlength*padlength/12.;
707 // gnoise2 = 0.0004/padlength;
709 // nprim = 0.133*amp;
710 // ggg2[i] = fkParam->GetDiffT()*fkParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
711 // glandau2[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
712 // if (glandau2[i]>1) glandau2[i]=1;
713 // glandau2[i]*=padlength*padlength/12.;
718 // gnoise3 = 0.0004/padlength;
720 // nprim = 0.133*amp;
721 // ggg3[i] = fkParam->GetDiffT()*fkParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
722 // glandau3[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
723 // if (glandau3[i]>1) glandau3[i]=1;
724 // glandau3[i]*=padlength*padlength/12.;
732 // Int_t amp = int(TMath::Abs(cl->GetQ()));
734 // seed->SetErrorY2(1.);
738 // Float_t z = TMath::Abs(fkParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
739 // Int_t ctype = cl->GetType();
740 // Float_t padlength= GetPadPitchLength(seed->GetRow());
741 // Double_t angle2 = seed->GetSnp()*seed->GetSnp();
742 // angle2 = angle2/(1-angle2);
744 // //cluster "quality"
745 // Int_t rsigmay = int(100.*cl->GetSigmaY2()/(seed->GetCurrentSigmaY2()));
748 // if (fSectors==fInnerSec){
749 // snoise2 = gnoise1;
750 // res = ggg1[amp]*z+glandau1[amp]*angle2;
751 // if (ctype==0) res *= gcor01[rsigmay];
754 // res*= gcorp[rsigmay];
758 // if (padlength<1.1){
759 // snoise2 = gnoise2;
760 // res = ggg2[amp]*z+glandau2[amp]*angle2;
761 // if (ctype==0) res *= gcor02[rsigmay];
764 // res*= gcorp[rsigmay];
768 // snoise2 = gnoise3;
769 // res = ggg3[amp]*z+glandau3[amp]*angle2;
770 // if (ctype==0) res *= gcor02[rsigmay];
773 // res*= gcorp[rsigmay];
780 // res*=2.4; // overestimate error 2 times
784 // if (res<2*snoise2)
787 // seed->SetErrorY2(res);
795 Double_t AliTPCtrackerMI::ErrZ2(AliTPCseed* seed, const AliTPCclusterMI * cl){
798 // Use calibrated cluster error from OCDB
800 AliTPCClusterParam * clparam = AliTPCcalibDB::Instance()->GetClusterParam();
802 Float_t z = TMath::Abs(fkParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
803 Int_t ctype = cl->GetType();
804 Int_t type = (cl->GetRow()<63) ? 0: (cl->GetRow()>126) ? 1:2;
806 Double_t angle2 = seed->GetSnp()*seed->GetSnp();
807 angle2 = seed->GetTgl()*seed->GetTgl()*(1+angle2/(1-angle2));
808 Double_t angle = TMath::Sqrt(TMath::Abs(angle2));
809 Double_t errz2 = clparam->GetError0Par(1,type, z,angle);
811 errz2+=0.5; // edge cluster
814 seed->SetErrorZ2(errz2);
820 // //seed->SetErrorY2(0.1);
822 // //calculate look-up table at the beginning
823 // static Bool_t ginit = kFALSE;
824 // static Float_t gnoise1,gnoise2,gnoise3;
825 // static Float_t ggg1[10000];
826 // static Float_t ggg2[10000];
827 // static Float_t ggg3[10000];
828 // static Float_t glandau1[10000];
829 // static Float_t glandau2[10000];
830 // static Float_t glandau3[10000];
832 // static Float_t gcor01[1000];
833 // static Float_t gcor02[1000];
834 // static Float_t gcorp[1000];
838 // if (ginit==kFALSE){
839 // for (Int_t i=1;i<1000;i++){
840 // Float_t rsigma = float(i)/100.;
841 // gcor02[i] = TMath::Max(0.81 +TMath::Exp(6.8*(rsigma-1.2)),0.6);
842 // gcor01[i] = TMath::Max(0.72 +TMath::Exp(2.04*(rsigma-1.2)),0.6);
843 // gcorp[i] = TMath::Max(TMath::Power((rsigma+0.5),1.5),1.2);
847 // for (Int_t i=3;i<10000;i++){
851 // Float_t amp = float(i);
852 // Float_t padlength =0.75;
853 // gnoise1 = 0.0004/padlength;
854 // Float_t nel = 0.268*amp;
855 // Float_t nprim = 0.155*amp;
856 // ggg1[i] = fkParam->GetDiffT()*fkParam->GetDiffT()*(2+0.001*nel/(padlength*padlength))/nel;
857 // glandau1[i] = (2.+0.12*nprim)*0.5* (2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
858 // if (glandau1[i]>1) glandau1[i]=1;
859 // glandau1[i]*=padlength*padlength/12.;
863 // gnoise2 = 0.0004/padlength;
865 // nprim = 0.133*amp;
866 // ggg2[i] = fkParam->GetDiffT()*fkParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
867 // glandau2[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
868 // if (glandau2[i]>1) glandau2[i]=1;
869 // glandau2[i]*=padlength*padlength/12.;
874 // gnoise3 = 0.0004/padlength;
876 // nprim = 0.133*amp;
877 // ggg3[i] = fkParam->GetDiffT()*fkParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
878 // glandau3[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
879 // if (glandau3[i]>1) glandau3[i]=1;
880 // glandau3[i]*=padlength*padlength/12.;
888 // Int_t amp = int(TMath::Abs(cl->GetQ()));
890 // seed->SetErrorY2(1.);
894 // Float_t z = TMath::Abs(fkParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
895 // Int_t ctype = cl->GetType();
896 // Float_t padlength= GetPadPitchLength(seed->GetRow());
898 // Double_t angle2 = seed->GetSnp()*seed->GetSnp();
899 // // if (angle2<0.6) angle2 = 0.6;
900 // angle2 = seed->GetTgl()*seed->GetTgl()*(1+angle2/(1-angle2));
902 // //cluster "quality"
903 // Int_t rsigmaz = int(100.*cl->GetSigmaZ2()/(seed->GetCurrentSigmaZ2()));
906 // if (fSectors==fInnerSec){
907 // snoise2 = gnoise1;
908 // res = ggg1[amp]*z+glandau1[amp]*angle2;
909 // if (ctype==0) res *= gcor01[rsigmaz];
912 // res*= gcorp[rsigmaz];
916 // if (padlength<1.1){
917 // snoise2 = gnoise2;
918 // res = ggg2[amp]*z+glandau2[amp]*angle2;
919 // if (ctype==0) res *= gcor02[rsigmaz];
922 // res*= gcorp[rsigmaz];
926 // snoise2 = gnoise3;
927 // res = ggg3[amp]*z+glandau3[amp]*angle2;
928 // if (ctype==0) res *= gcor02[rsigmaz];
931 // res*= gcorp[rsigmaz];
940 // if ((ctype<0) &&<70){
945 // if (res<2*snoise2)
947 // if (res>3) res =3;
948 // seed->SetErrorZ2(res);
956 void AliTPCtrackerMI::RotateToLocal(AliTPCseed *seed)
958 //rotate to track "local coordinata
959 Float_t x = seed->GetX();
960 Float_t y = seed->GetY();
961 Float_t ymax = x*TMath::Tan(0.5*fSectors->GetAlpha());
964 seed->SetRelativeSector((seed->GetRelativeSector()+1) % fN);
965 if (!seed->Rotate(fSectors->GetAlpha()))
967 } else if (y <-ymax) {
968 seed->SetRelativeSector((seed->GetRelativeSector()-1+fN) % fN);
969 if (!seed->Rotate(-fSectors->GetAlpha()))
977 //_____________________________________________________________________________
978 Double_t AliTPCtrackerMI::F1old(Double_t x1,Double_t y1,
979 Double_t x2,Double_t y2,
980 Double_t x3,Double_t y3) const
982 //-----------------------------------------------------------------
983 // Initial approximation of the track curvature
984 //-----------------------------------------------------------------
985 Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
986 Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
987 (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
988 Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
989 (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
991 Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
992 if ( xr*xr+yr*yr<=0.00000000000001) return 100;
993 return -xr*yr/sqrt(xr*xr+yr*yr);
998 //_____________________________________________________________________________
999 Double_t AliTPCtrackerMI::F1(Double_t x1,Double_t y1,
1000 Double_t x2,Double_t y2,
1001 Double_t x3,Double_t y3) const
1003 //-----------------------------------------------------------------
1004 // Initial approximation of the track curvature
1005 //-----------------------------------------------------------------
1011 Double_t det = x3*y2-x2*y3;
1012 if (TMath::Abs(det)<1e-10){
1016 Double_t u = 0.5* (x2*(x2-x3)+y2*(y2-y3))/det;
1017 Double_t x0 = x3*0.5-y3*u;
1018 Double_t y0 = y3*0.5+x3*u;
1019 Double_t c2 = 1/TMath::Sqrt(x0*x0+y0*y0);
1025 Double_t AliTPCtrackerMI::F2(Double_t x1,Double_t y1,
1026 Double_t x2,Double_t y2,
1027 Double_t x3,Double_t y3) const
1029 //-----------------------------------------------------------------
1030 // Initial approximation of the track curvature
1031 //-----------------------------------------------------------------
1037 Double_t det = x3*y2-x2*y3;
1038 if (TMath::Abs(det)<1e-10) {
1042 Double_t u = 0.5* (x2*(x2-x3)+y2*(y2-y3))/det;
1043 Double_t x0 = x3*0.5-y3*u;
1044 Double_t y0 = y3*0.5+x3*u;
1045 Double_t c2 = 1/TMath::Sqrt(x0*x0+y0*y0);
1054 //_____________________________________________________________________________
1055 Double_t AliTPCtrackerMI::F2old(Double_t x1,Double_t y1,
1056 Double_t x2,Double_t y2,
1057 Double_t x3,Double_t y3) const
1059 //-----------------------------------------------------------------
1060 // Initial approximation of the track curvature times center of curvature
1061 //-----------------------------------------------------------------
1062 Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
1063 Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
1064 (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
1065 Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
1066 (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
1068 Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
1070 return -a/(d*y1-b)*xr/sqrt(xr*xr+yr*yr);
1073 //_____________________________________________________________________________
1074 Double_t AliTPCtrackerMI::F3(Double_t x1,Double_t y1,
1075 Double_t x2,Double_t y2,
1076 Double_t z1,Double_t z2) const
1078 //-----------------------------------------------------------------
1079 // Initial approximation of the tangent of the track dip angle
1080 //-----------------------------------------------------------------
1081 return (z1 - z2)/sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
1085 Double_t AliTPCtrackerMI::F3n(Double_t x1,Double_t y1,
1086 Double_t x2,Double_t y2,
1087 Double_t z1,Double_t z2, Double_t c) const
1089 //-----------------------------------------------------------------
1090 // Initial approximation of the tangent of the track dip angle
1091 //-----------------------------------------------------------------
1095 //angle1 = (z1-z2)*c/(TMath::ASin(c*x1-ni)-TMath::ASin(c*x2-ni));
1097 Double_t d = TMath::Sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
1098 if (TMath::Abs(d*c*0.5)>1) return 0;
1099 // Double_t angle2 = TMath::ASin(d*c*0.5);
1100 // Double_t angle2 = AliTPCFastMath::FastAsin(d*c*0.5);
1101 Double_t angle2 = (d*c*0.5>0.1)? TMath::ASin(d*c*0.5): AliTPCFastMath::FastAsin(d*c*0.5);
1103 angle2 = (z1-z2)*c/(angle2*2.);
1107 Bool_t AliTPCtrackerMI::GetProlongation(Double_t x1, Double_t x2, Double_t x[5], Double_t &y, Double_t &z) const
1108 {//-----------------------------------------------------------------
1109 // This function find proloncation of a track to a reference plane x=x2.
1110 //-----------------------------------------------------------------
1114 if (TMath::Abs(x[4]*x1 - x[2]) >= 0.999) {
1118 Double_t c1=x[4]*x1 - x[2], r1=TMath::Sqrt((1.-c1)*(1.+c1));
1119 Double_t c2=x[4]*x2 - x[2], r2=TMath::Sqrt((1.-c2)*(1.+c2));
1123 Double_t dy = dx*(c1+c2)/(r1+r2);
1126 Double_t delta = x[4]*dx*(c1+c2)/(c1*r2 + c2*r1);
1128 if (TMath::Abs(delta)>0.01){
1129 dz = x[3]*TMath::ASin(delta)/x[4];
1131 dz = x[3]*AliTPCFastMath::FastAsin(delta)/x[4];
1134 //dz = x[3]*AliTPCFastMath::FastAsin(delta)/x[4];
1142 Int_t AliTPCtrackerMI::LoadClusters (TTree *const tree)
1147 return LoadClusters();
1151 Int_t AliTPCtrackerMI::LoadClusters(const TObjArray *arr)
1154 // load clusters to the memory
1155 AliTPCClustersRow *clrow = new AliTPCClustersRow("AliTPCclusterMI");
1156 Int_t lower = arr->LowerBound();
1157 Int_t entries = arr->GetEntriesFast();
1159 for (Int_t i=lower; i<entries; i++) {
1160 clrow = (AliTPCClustersRow*) arr->At(i);
1161 if(!clrow) continue;
1162 if(!clrow->GetArray()) continue;
1166 fkParam->AdjustSectorRow(clrow->GetID(),sec,row);
1168 for (Int_t icl=0; icl<clrow->GetArray()->GetEntriesFast(); icl++){
1169 Transform((AliTPCclusterMI*)(clrow->GetArray()->At(icl)));
1172 if (clrow->GetArray()->GetEntriesFast()<=0) continue;
1173 AliTPCtrackerRow * tpcrow=0;
1176 tpcrow = &(fInnerSec[sec%fkNIS][row]);
1180 tpcrow = &(fOuterSec[(sec-fkNIS*2)%fkNOS][row]);
1181 left = (sec-fkNIS*2)/fkNOS;
1184 tpcrow->SetN1(clrow->GetArray()->GetEntriesFast());
1185 tpcrow->SetClusters1(new AliTPCclusterMI[tpcrow->GetN1()]);
1186 for (Int_t j=0;j<tpcrow->GetN1();++j)
1187 tpcrow->SetCluster1(j, *(AliTPCclusterMI*)(clrow->GetArray()->At(j)));
1190 tpcrow->SetN2(clrow->GetArray()->GetEntriesFast());
1191 tpcrow->SetClusters2(new AliTPCclusterMI[tpcrow->GetN2()]);
1192 for (Int_t j=0;j<tpcrow->GetN2();++j)
1193 tpcrow->SetCluster2(j,*(AliTPCclusterMI*)(clrow->GetArray()->At(j)));
1195 clrow->GetArray()->Clear("C");
1204 Int_t AliTPCtrackerMI::LoadClusters(const TClonesArray *arr)
1207 // load clusters to the memory from one
1210 AliTPCclusterMI *clust=0;
1211 Int_t count[72][96] = { {0} , {0} };
1213 // loop over clusters
1214 for (Int_t icl=0; icl<arr->GetEntriesFast(); icl++) {
1215 clust = (AliTPCclusterMI*)arr->At(icl);
1216 if(!clust) continue;
1217 //printf("cluster: det %d, row %d \n", clust->GetDetector(),clust->GetRow());
1219 // transform clusters
1222 // count clusters per pad row
1223 count[clust->GetDetector()][clust->GetRow()]++;
1226 // insert clusters to sectors
1227 for (Int_t icl=0; icl<arr->GetEntriesFast(); icl++) {
1228 clust = (AliTPCclusterMI*)arr->At(icl);
1229 if(!clust) continue;
1231 Int_t sec = clust->GetDetector();
1232 Int_t row = clust->GetRow();
1234 // filter overlapping pad rows needed by HLT
1235 if(sec<fkNIS*2) { //IROCs
1236 if(row == 30) continue;
1239 if(row == 27 || row == 76) continue;
1245 fInnerSec[sec%fkNIS].InsertCluster(clust, count[sec][row], fkParam);
1248 left = (sec-fkNIS*2)/fkNOS;
1249 fOuterSec[(sec-fkNIS*2)%fkNOS].InsertCluster(clust, count[sec][row], fkParam);
1253 // Load functions must be called behind LoadCluster(TClonesArray*)
1255 //LoadOuterSectors();
1256 //LoadInnerSectors();
1262 Int_t AliTPCtrackerMI::LoadClusters()
1265 // load clusters to the memory
1266 AliTPCClustersRow *clrow= new AliTPCClustersRow("AliTPCclusterMI");
1268 // TTree * tree = fClustersArray.GetTree();
1270 TTree * tree = fInput;
1271 TBranch * br = tree->GetBranch("Segment");
1272 br->SetAddress(&clrow);
1274 // Conversion of pad, row coordinates in local tracking coords.
1275 // Could be skipped here; is already done in clusterfinder
1277 Int_t j=Int_t(tree->GetEntries());
1278 for (Int_t i=0; i<j; i++) {
1282 fkParam->AdjustSectorRow(clrow->GetID(),sec,row);
1283 for (Int_t icl=0; icl<clrow->GetArray()->GetEntriesFast(); icl++){
1284 Transform((AliTPCclusterMI*)(clrow->GetArray()->At(icl)));
1287 AliTPCtrackerRow * tpcrow=0;
1290 tpcrow = &(fInnerSec[sec%fkNIS][row]);
1294 tpcrow = &(fOuterSec[(sec-fkNIS*2)%fkNOS][row]);
1295 left = (sec-fkNIS*2)/fkNOS;
1298 tpcrow->SetN1(clrow->GetArray()->GetEntriesFast());
1299 tpcrow->SetClusters1(new AliTPCclusterMI[tpcrow->GetN1()]);
1300 for (Int_t k=0;k<tpcrow->GetN1();++k)
1301 tpcrow->SetCluster1(k, *(AliTPCclusterMI*)(clrow->GetArray()->At(k)));
1304 tpcrow->SetN2(clrow->GetArray()->GetEntriesFast());
1305 tpcrow->SetClusters2(new AliTPCclusterMI[tpcrow->GetN2()]);
1306 for (Int_t k=0;k<tpcrow->GetN2();++k)
1307 tpcrow->SetCluster2(k,*(AliTPCclusterMI*)(clrow->GetArray()->At(k)));
1318 void AliTPCtrackerMI::UnloadClusters()
1321 // unload clusters from the memory
1323 Int_t nrows = fOuterSec->GetNRows();
1324 for (Int_t sec = 0;sec<fkNOS;sec++)
1325 for (Int_t row = 0;row<nrows;row++){
1326 AliTPCtrackerRow* tpcrow = &(fOuterSec[sec%fkNOS][row]);
1328 // if (tpcrow->fClusters1) delete []tpcrow->fClusters1;
1329 // if (tpcrow->fClusters2) delete []tpcrow->fClusters2;
1331 tpcrow->ResetClusters();
1334 nrows = fInnerSec->GetNRows();
1335 for (Int_t sec = 0;sec<fkNIS;sec++)
1336 for (Int_t row = 0;row<nrows;row++){
1337 AliTPCtrackerRow* tpcrow = &(fInnerSec[sec%fkNIS][row]);
1339 // if (tpcrow->fClusters1) delete []tpcrow->fClusters1;
1340 //if (tpcrow->fClusters2) delete []tpcrow->fClusters2;
1342 tpcrow->ResetClusters();
1348 void AliTPCtrackerMI::FillClusterArray(TObjArray* array) const{
1350 // Filling cluster to the array - For visualization purposes
1353 nrows = fOuterSec->GetNRows();
1354 for (Int_t sec = 0;sec<fkNOS;sec++)
1355 for (Int_t row = 0;row<nrows;row++){
1356 AliTPCtrackerRow* tpcrow = &(fOuterSec[sec%fkNOS][row]);
1357 if (!tpcrow) continue;
1358 for (Int_t icl = 0;icl<tpcrow->GetN();icl++){
1359 array->AddLast((TObject*)((*tpcrow)[icl]));
1362 nrows = fInnerSec->GetNRows();
1363 for (Int_t sec = 0;sec<fkNIS;sec++)
1364 for (Int_t row = 0;row<nrows;row++){
1365 AliTPCtrackerRow* tpcrow = &(fInnerSec[sec%fkNIS][row]);
1366 if (!tpcrow) continue;
1367 for (Int_t icl = 0;icl<tpcrow->GetN();icl++){
1368 array->AddLast((TObject*)(*tpcrow)[icl]);
1374 void AliTPCtrackerMI::Transform(AliTPCclusterMI * cluster){
1378 AliTPCcalibDB * calibDB = AliTPCcalibDB::Instance();
1379 AliTPCTransform *transform = calibDB->GetTransform() ;
1381 AliFatal("Tranformations not in calibDB");
1384 transform->SetCurrentRecoParam((AliTPCRecoParam*)AliTPCReconstructor::GetRecoParam());
1385 Double_t x[3]={cluster->GetRow(),cluster->GetPad(),cluster->GetTimeBin()};
1386 Int_t i[1]={cluster->GetDetector()};
1387 transform->Transform(x,i,0,1);
1388 // if (cluster->GetDetector()%36>17){
1393 // in debug mode check the transformation
1395 if (AliTPCReconstructor::StreamLevel()>2) {
1397 cluster->GetGlobalXYZ(gx);
1398 Int_t event = (fEvent==NULL)? 0: fEvent->GetEventNumberInFile();
1399 TTreeSRedirector &cstream = *fDebugStreamer;
1400 cstream<<"Transform"<<
1411 cluster->SetX(x[0]);
1412 cluster->SetY(x[1]);
1413 cluster->SetZ(x[2]);
1418 //if (!fkParam->IsGeoRead()) fkParam->ReadGeoMatrices();
1419 if (AliTPCReconstructor::GetRecoParam()->GetUseSectorAlignment() && (!calibDB->HasAlignmentOCDB())){
1420 TGeoHMatrix *mat = fkParam->GetClusterMatrix(cluster->GetDetector());
1422 Double_t pos[3]= {cluster->GetX(),cluster->GetY(),cluster->GetZ()};
1423 Double_t posC[3]={cluster->GetX(),cluster->GetY(),cluster->GetZ()};
1424 if (mat) mat->LocalToMaster(pos,posC);
1426 // chack Loading of Geo matrices from GeoManager - TEMPORARY FIX
1428 cluster->SetX(posC[0]);
1429 cluster->SetY(posC[1]);
1430 cluster->SetZ(posC[2]);
1434 //_____________________________________________________________________________
1435 Int_t AliTPCtrackerMI::LoadOuterSectors() {
1436 //-----------------------------------------------------------------
1437 // This function fills outer TPC sectors with clusters.
1438 //-----------------------------------------------------------------
1439 Int_t nrows = fOuterSec->GetNRows();
1441 for (Int_t sec = 0;sec<fkNOS;sec++)
1442 for (Int_t row = 0;row<nrows;row++){
1443 AliTPCtrackerRow* tpcrow = &(fOuterSec[sec%fkNOS][row]);
1444 Int_t sec2 = sec+2*fkNIS;
1446 Int_t ncl = tpcrow->GetN1();
1448 AliTPCclusterMI *c= (tpcrow->GetCluster1(ncl));
1449 index=(((sec2<<8)+row)<<16)+ncl;
1450 tpcrow->InsertCluster(c,index);
1453 ncl = tpcrow->GetN2();
1455 AliTPCclusterMI *c= (tpcrow->GetCluster2(ncl));
1456 index=((((sec2+fkNOS)<<8)+row)<<16)+ncl;
1457 tpcrow->InsertCluster(c,index);
1460 // write indexes for fast acces
1462 for (Int_t i=0;i<510;i++)
1463 tpcrow->SetFastCluster(i,-1);
1464 for (Int_t i=0;i<tpcrow->GetN();i++){
1465 Int_t zi = Int_t((*tpcrow)[i]->GetZ()+255.);
1466 tpcrow->SetFastCluster(zi,i); // write index
1469 for (Int_t i=0;i<510;i++){
1470 if (tpcrow->GetFastCluster(i)<0)
1471 tpcrow->SetFastCluster(i,last);
1473 last = tpcrow->GetFastCluster(i);
1482 //_____________________________________________________________________________
1483 Int_t AliTPCtrackerMI::LoadInnerSectors() {
1484 //-----------------------------------------------------------------
1485 // This function fills inner TPC sectors with clusters.
1486 //-----------------------------------------------------------------
1487 Int_t nrows = fInnerSec->GetNRows();
1489 for (Int_t sec = 0;sec<fkNIS;sec++)
1490 for (Int_t row = 0;row<nrows;row++){
1491 AliTPCtrackerRow* tpcrow = &(fInnerSec[sec%fkNIS][row]);
1494 Int_t ncl = tpcrow->GetN1();
1496 AliTPCclusterMI *c= (tpcrow->GetCluster1(ncl));
1497 index=(((sec<<8)+row)<<16)+ncl;
1498 tpcrow->InsertCluster(c,index);
1501 ncl = tpcrow->GetN2();
1503 AliTPCclusterMI *c= (tpcrow->GetCluster2(ncl));
1504 index=((((sec+fkNIS)<<8)+row)<<16)+ncl;
1505 tpcrow->InsertCluster(c,index);
1508 // write indexes for fast acces
1510 for (Int_t i=0;i<510;i++)
1511 tpcrow->SetFastCluster(i,-1);
1512 for (Int_t i=0;i<tpcrow->GetN();i++){
1513 Int_t zi = Int_t((*tpcrow)[i]->GetZ()+255.);
1514 tpcrow->SetFastCluster(zi,i); // write index
1517 for (Int_t i=0;i<510;i++){
1518 if (tpcrow->GetFastCluster(i)<0)
1519 tpcrow->SetFastCluster(i,last);
1521 last = tpcrow->GetFastCluster(i);
1533 //_________________________________________________________________________
1534 AliTPCclusterMI *AliTPCtrackerMI::GetClusterMI(Int_t index) const {
1535 //--------------------------------------------------------------------
1536 // Return pointer to a given cluster
1537 //--------------------------------------------------------------------
1538 if (index<0) return 0; // no cluster
1539 Int_t sec=(index&0xff000000)>>24;
1540 Int_t row=(index&0x00ff0000)>>16;
1541 Int_t ncl=(index&0x00007fff)>>00;
1543 const AliTPCtrackerRow * tpcrow=0;
1544 AliTPCclusterMI * clrow =0;
1546 if (sec<0 || sec>=fkNIS*4) {
1547 AliWarning(Form("Wrong sector %d",sec));
1552 AliTPCtrackerSector& tracksec = fInnerSec[sec%fkNIS];
1553 if (tracksec.GetNRows()<=row) return 0;
1554 tpcrow = &(tracksec[row]);
1555 if (tpcrow==0) return 0;
1558 if (tpcrow->GetN1()<=ncl) return 0;
1559 clrow = tpcrow->GetClusters1();
1562 if (tpcrow->GetN2()<=ncl) return 0;
1563 clrow = tpcrow->GetClusters2();
1567 AliTPCtrackerSector& tracksec = fOuterSec[(sec-fkNIS*2)%fkNOS];
1568 if (tracksec.GetNRows()<=row) return 0;
1569 tpcrow = &(tracksec[row]);
1570 if (tpcrow==0) return 0;
1572 if (sec-2*fkNIS<fkNOS) {
1573 if (tpcrow->GetN1()<=ncl) return 0;
1574 clrow = tpcrow->GetClusters1();
1577 if (tpcrow->GetN2()<=ncl) return 0;
1578 clrow = tpcrow->GetClusters2();
1582 return &(clrow[ncl]);
1588 Int_t AliTPCtrackerMI::FollowToNext(AliTPCseed& t, Int_t nr) {
1589 //-----------------------------------------------------------------
1590 // This function tries to find a track prolongation to next pad row
1591 //-----------------------------------------------------------------
1593 Double_t x= GetXrow(nr), ymax=GetMaxY(nr);
1596 AliTPCclusterMI *cl=0;
1597 Int_t tpcindex= t.GetClusterIndex2(nr);
1599 // update current shape info every 5 pad-row
1600 // if ( (nr%5==0) || t.GetNumberOfClusters()<2 || (t.fCurrentSigmaY2<0.0001) ){
1604 if (fIteration>0 && tpcindex>=-1){ //if we have already clusters
1606 if (tpcindex==-1) return 0; //track in dead zone
1607 if (tpcindex >= 0){ //
1608 cl = t.GetClusterPointer(nr);
1609 if ( (cl==0) ) cl = GetClusterMI(tpcindex);
1610 t.SetCurrentClusterIndex1(tpcindex);
1613 Int_t relativesector = ((tpcindex&0xff000000)>>24)%18; // if previously accepted cluster in different sector
1614 Float_t angle = relativesector*fSectors->GetAlpha()+fSectors->GetAlphaShift();
1616 if (angle<-TMath::Pi()) angle += 2*TMath::Pi();
1617 if (angle>=TMath::Pi()) angle -= 2*TMath::Pi();
1619 if (TMath::Abs(angle-t.GetAlpha())>0.001){
1620 Double_t rotation = angle-t.GetAlpha();
1621 t.SetRelativeSector(relativesector);
1622 if (!t.Rotate(rotation)) return 0;
1624 if (!t.PropagateTo(x)) return 0;
1626 t.SetCurrentCluster(cl);
1628 Int_t accept = AcceptCluster(&t,t.GetCurrentCluster());
1629 if ((tpcindex&0x8000)==0) accept =0;
1631 //if founded cluster is acceptible
1632 if (cl->IsUsed(11)) { // id cluster is shared inrease uncertainty
1633 t.SetErrorY2(t.GetErrorY2()+0.03);
1634 t.SetErrorZ2(t.GetErrorZ2()+0.03);
1635 t.SetErrorY2(t.GetErrorY2()*3);
1636 t.SetErrorZ2(t.GetErrorZ2()*3);
1638 t.SetNFoundable(t.GetNFoundable()+1);
1639 UpdateTrack(&t,accept);
1642 else { // Remove old cluster from track
1643 t.SetClusterIndex(nr, -3);
1644 t.SetClusterPointer(nr, 0);
1648 if (TMath::Abs(t.GetSnp())>AliTPCReconstructor::GetMaxSnpTracker()) return 0; // cut on angle
1649 if (fIteration>1 && IsFindable(t)){
1650 // not look for new cluster during refitting
1651 t.SetNFoundable(t.GetNFoundable()+1);
1656 // if (TMath::Abs(t.GetSnp())>0.95 || TMath::Abs(x*t.GetC()-t.GetEta())>0.95) return 0;// patch 28 fev 06
1657 if (!t.PropagateTo(x)) {
1658 if (fIteration==0) t.SetRemoval(10);
1661 Double_t y = t.GetY();
1662 if (TMath::Abs(y)>ymax){
1664 t.SetRelativeSector((t.GetRelativeSector()+1) % fN);
1665 if (!t.Rotate(fSectors->GetAlpha()))
1667 } else if (y <-ymax) {
1668 t.SetRelativeSector((t.GetRelativeSector()-1+fN) % fN);
1669 if (!t.Rotate(-fSectors->GetAlpha()))
1672 if (!t.PropagateTo(x)) {
1673 if (fIteration==0) t.SetRemoval(10);
1679 Double_t z=t.GetZ();
1682 if (!IsActive(t.GetRelativeSector(),nr)) {
1684 t.SetClusterIndex2(nr,-1);
1687 //AliInfo(Form("A - Sector%d phi %f - alpha %f", t.fRelativeSector,y/x, t.GetAlpha()));
1688 Bool_t isActive = IsActive(t.GetRelativeSector(),nr);
1689 Bool_t isActive2 = (nr>=fInnerSec->GetNRows()) ? fOuterSec[t.GetRelativeSector()][nr-fInnerSec->GetNRows()].GetN()>0:fInnerSec[t.GetRelativeSector()][nr].GetN()>0;
1691 if (!isActive || !isActive2) return 0;
1693 const AliTPCtrackerRow &krow=GetRow(t.GetRelativeSector(),nr);
1694 if ( (t.GetSigmaY2()<0) || t.GetSigmaZ2()<0) return 0;
1696 Double_t roadz = 1.;
1698 if (TMath::Abs(TMath::Abs(y)-ymax)<krow.GetDeadZone()){
1700 t.SetClusterIndex2(nr,-1);
1706 // if (TMath::Abs(z)<(AliTPCReconstructor::GetCtgRange()*x+10) && TMath::Abs(z)<fkParam->GetZLength(0) && (TMath::Abs(t.GetSnp())<AliTPCReconstructor::GetMaxSnpTracker()))
1707 t.SetNFoundable(t.GetNFoundable()+1);
1713 // cl = krow.FindNearest2(y+10.,z,roady,roadz,index);
1714 cl = krow.FindNearest2(y,z,roady,roadz,index);
1715 if (cl) t.SetCurrentClusterIndex1(krow.GetIndex(index));
1718 t.SetCurrentCluster(cl);
1720 if (fIteration==2&&cl->IsUsed(10)) return 0;
1721 Int_t accept = AcceptCluster(&t,t.GetCurrentCluster());
1722 if (fIteration==2&&cl->IsUsed(11)) {
1723 t.SetErrorY2(t.GetErrorY2()+0.03);
1724 t.SetErrorZ2(t.GetErrorZ2()+0.03);
1725 t.SetErrorY2(t.GetErrorY2()*3);
1726 t.SetErrorZ2(t.GetErrorZ2()*3);
1729 if (t.fCurrentCluster->IsUsed(10)){
1734 if (t.fNShared>0.7*t.GetNumberOfClusters()) {
1740 if (accept<3) UpdateTrack(&t,accept);
1743 if ( fIteration==0 && t.GetNFoundable()*0.5 > t.GetNumberOfClusters()) t.SetRemoval(10);
1751 //_________________________________________________________________________
1752 Bool_t AliTPCtrackerMI::GetTrackPoint(Int_t index, AliTrackPoint &p ) const
1754 // Get track space point by index
1755 // return false in case the cluster doesn't exist
1756 AliTPCclusterMI *cl = GetClusterMI(index);
1757 if (!cl) return kFALSE;
1758 Int_t sector = (index&0xff000000)>>24;
1759 // Int_t row = (index&0x00ff0000)>>16;
1761 // xyz[0] = fkParam->GetPadRowRadii(sector,row);
1762 xyz[0] = cl->GetX();
1763 xyz[1] = cl->GetY();
1764 xyz[2] = cl->GetZ();
1766 fkParam->AdjustCosSin(sector,cos,sin);
1767 Float_t x = cos*xyz[0]-sin*xyz[1];
1768 Float_t y = cos*xyz[1]+sin*xyz[0];
1770 Float_t sigmaY2 = 0.027*cl->GetSigmaY2();
1771 if (sector < fkParam->GetNInnerSector()) sigmaY2 *= 2.07;
1772 Float_t sigmaZ2 = 0.066*cl->GetSigmaZ2();
1773 if (sector < fkParam->GetNInnerSector()) sigmaZ2 *= 1.77;
1774 cov[0] = sin*sin*sigmaY2;
1775 cov[1] = -sin*cos*sigmaY2;
1777 cov[3] = cos*cos*sigmaY2;
1780 p.SetXYZ(x,y,xyz[2],cov);
1781 AliGeomManager::ELayerID iLayer;
1783 if (sector < fkParam->GetNInnerSector()) {
1784 iLayer = AliGeomManager::kTPC1;
1788 iLayer = AliGeomManager::kTPC2;
1789 idet = sector - fkParam->GetNInnerSector();
1791 UShort_t volid = AliGeomManager::LayerToVolUID(iLayer,idet);
1792 p.SetVolumeID(volid);
1798 Int_t AliTPCtrackerMI::UpdateClusters(AliTPCseed& t, Int_t nr) {
1799 //-----------------------------------------------------------------
1800 // This function tries to find a track prolongation to next pad row
1801 //-----------------------------------------------------------------
1802 t.SetCurrentCluster(0);
1803 t.SetCurrentClusterIndex1(-3);
1805 Double_t xt=t.GetX();
1806 Int_t row = GetRowNumber(xt)-1;
1807 Double_t ymax= GetMaxY(nr);
1809 if (row < nr) return 1; // don't prolongate if not information until now -
1810 // if (TMath::Abs(t.GetSnp())>0.9 && t.GetNumberOfClusters()>40. && fIteration!=2) {
1812 // return 0; // not prolongate strongly inclined tracks
1814 // if (TMath::Abs(t.GetSnp())>0.95) {
1816 // return 0; // not prolongate strongly inclined tracks
1817 // }// patch 28 fev 06
1819 Double_t x= GetXrow(nr);
1821 //t.PropagateTo(x+0.02);
1822 //t.PropagateTo(x+0.01);
1823 if (!t.PropagateTo(x)){
1830 if (TMath::Abs(y)>ymax){
1832 t.SetRelativeSector((t.GetRelativeSector()+1) % fN);
1833 if (!t.Rotate(fSectors->GetAlpha()))
1835 } else if (y <-ymax) {
1836 t.SetRelativeSector((t.GetRelativeSector()-1+fN) % fN);
1837 if (!t.Rotate(-fSectors->GetAlpha()))
1840 // if (!t.PropagateTo(x)){
1847 if (TMath::Abs(t.GetSnp())>AliTPCReconstructor::GetMaxSnpTracker()) return 0;
1849 if (!IsActive(t.GetRelativeSector(),nr)) {
1851 t.SetClusterIndex2(nr,-1);
1854 //AliInfo(Form("A - Sector%d phi %f - alpha %f", t.fRelativeSector,y/x, t.GetAlpha()));
1856 AliTPCtrackerRow &krow=GetRow(t.GetRelativeSector(),nr);
1858 if (TMath::Abs(TMath::Abs(y)-ymax)<krow.GetDeadZone()){
1860 t.SetClusterIndex2(nr,-1);
1866 // if (TMath::Abs(t.GetZ())<(AliTPCReconstructor::GetCtgRange()*t.GetX()+10) && (TMath::Abs(t.GetSnp())<AliTPCReconstructor::GetMaxSnpTracker()))
1867 if (IsFindable(t)) t.SetNFoundable(t.GetNFoundable()+1);
1873 if ( (nr%2==0) || t.GetNumberOfClusters()<2){
1874 // t.fCurrentSigmaY = GetSigmaY(&t);
1875 //t.fCurrentSigmaZ = GetSigmaZ(&t);
1879 AliTPCclusterMI *cl=0;
1882 Double_t roady = 1.;
1883 Double_t roadz = 1.;
1887 index = t.GetClusterIndex2(nr);
1888 if ( (index >= 0) && (index&0x8000)==0){
1889 cl = t.GetClusterPointer(nr);
1890 if ( (cl==0) && (index >= 0)) cl = GetClusterMI(index);
1891 t.SetCurrentClusterIndex1(index);
1893 t.SetCurrentCluster(cl);
1899 // if (index<0) return 0;
1900 UInt_t uindex = TMath::Abs(index);
1903 //cl = krow.FindNearest2(y+10,z,roady,roadz,uindex);
1904 cl = krow.FindNearest2(y,z,roady,roadz,uindex);
1907 if (cl) t.SetCurrentClusterIndex1(krow.GetIndex(uindex));
1908 t.SetCurrentCluster(cl);
1914 Int_t AliTPCtrackerMI::FollowToNextCluster(AliTPCseed & t, Int_t nr) {
1915 //-----------------------------------------------------------------
1916 // This function tries to find a track prolongation to next pad row
1917 //-----------------------------------------------------------------
1919 //update error according neighborhoud
1921 if (t.GetCurrentCluster()) {
1923 Int_t accept = AcceptCluster(&t,t.GetCurrentCluster());
1925 if (t.GetCurrentCluster()->IsUsed(10)){
1930 t.SetNShared(t.GetNShared()+1);
1931 if (t.GetNShared()>0.7*t.GetNumberOfClusters()) {
1936 if (fIteration>0) accept = 0;
1937 if (accept<3) UpdateTrack(&t,accept);
1941 if ( t.GetNumberOfClusters()>18 && ( (t.GetSigmaY2()+t.GetSigmaZ2())>0.16)) t.SetRemoval(10);
1942 if ( t.GetNumberOfClusters()>18 && t.GetChi2()/t.GetNumberOfClusters()>6 ) t.SetRemoval(10);
1944 if (( (t.GetNFoundable()*0.5 > t.GetNumberOfClusters()) || t.GetNoCluster()>15)) t.SetRemoval(10);
1952 //_____________________________________________________________________________
1953 Int_t AliTPCtrackerMI::FollowProlongation(AliTPCseed& t, Int_t rf, Int_t step, Bool_t fromSeeds) {
1954 //-----------------------------------------------------------------
1955 // This function tries to find a track prolongation.
1956 //-----------------------------------------------------------------
1957 Double_t xt=t.GetX();
1959 Double_t alpha=t.GetAlpha();
1960 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
1961 if (alpha < 0. ) alpha += 2.*TMath::Pi();
1963 t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN);
1965 Int_t first = GetRowNumber(xt);
1970 for (Int_t nr= first; nr>=rf; nr-=step) {
1972 if (t.GetKinkIndexes()[0]>0){
1973 for (Int_t i=0;i<3;i++){
1974 Int_t index = t.GetKinkIndexes()[i];
1975 if (index==0) break;
1976 if (index<0) continue;
1978 AliKink * kink = (AliKink*)fEvent->GetKink(index-1);
1980 printf("PROBLEM\n");
1983 Int_t kinkrow = kink->GetTPCRow0()+2+Int_t(0.5/(0.05+kink->GetAngle(2)));
1985 AliExternalTrackParam paramd(t);
1986 kink->SetDaughter(paramd);
1987 kink->SetStatus(2,5);
1994 if (nr==80) t.UpdateReference();
1995 if (nr<fInnerSec->GetNRows())
1996 fSectors = fInnerSec;
1998 fSectors = fOuterSec;
1999 if (FollowToNext(t,nr)==0)
2012 Int_t AliTPCtrackerMI::FollowBackProlongation(AliTPCseed& t, Int_t rf, Bool_t fromSeeds) {
2013 //-----------------------------------------------------------------
2014 // This function tries to find a track prolongation.
2015 //-----------------------------------------------------------------
2017 Double_t xt=t.GetX();
2018 Double_t alpha=t.GetAlpha();
2019 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
2020 if (alpha < 0. ) alpha += 2.*TMath::Pi();
2021 t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN);
2023 Int_t first = t.GetFirstPoint();
2024 Int_t ri = GetRowNumber(xt);
2028 if (first<ri) first = ri;
2030 if (first<0) first=0;
2031 for (Int_t nr=first; nr<=rf; nr++) {
2032 // if ( (TMath::Abs(t.GetSnp())>0.95)) break;//patch 28 fev 06
2033 if (t.GetKinkIndexes()[0]<0){
2034 for (Int_t i=0;i<3;i++){
2035 Int_t index = t.GetKinkIndexes()[i];
2036 if (index==0) break;
2037 if (index>0) continue;
2038 index = TMath::Abs(index);
2039 AliKink * kink = (AliKink*)fEvent->GetKink(index-1);
2041 printf("PROBLEM\n");
2044 Int_t kinkrow = kink->GetTPCRow0()-2-Int_t(0.5/(0.05+kink->GetAngle(2)));
2046 AliExternalTrackParam paramm(t);
2047 kink->SetMother(paramm);
2048 kink->SetStatus(2,1);
2055 if (nr<fInnerSec->GetNRows())
2056 fSectors = fInnerSec;
2058 fSectors = fOuterSec;
2069 Float_t AliTPCtrackerMI::OverlapFactor(AliTPCseed * s1, AliTPCseed * s2, Int_t &sum1, Int_t & sum2)
2071 // overlapping factor
2077 Float_t dz2 =(s1->GetZ() - s2->GetZ());
2080 Float_t dy2 =TMath::Abs((s1->GetY() - s2->GetY()));
2082 Float_t distance = TMath::Sqrt(dz2+dy2);
2083 if (distance>4.) return 0; // if there are far away - not overlap - to reduce combinatorics
2086 Int_t firstpoint = TMath::Min(s1->GetFirstPoint(),s2->GetFirstPoint());
2087 Int_t lastpoint = TMath::Max(s1->GetLastPoint(),s2->GetLastPoint());
2092 if (firstpoint>lastpoint) {
2093 firstpoint =lastpoint;
2098 for (Int_t i=firstpoint-1;i<lastpoint+1;i++){
2099 if (s1->GetClusterIndex2(i)>0) sum1++;
2100 if (s2->GetClusterIndex2(i)>0) sum2++;
2101 if (s1->GetClusterIndex2(i)==s2->GetClusterIndex2(i) && s1->GetClusterIndex2(i)>0) {
2105 if (sum<5) return 0;
2107 Float_t summin = TMath::Min(sum1+1,sum2+1);
2108 Float_t ratio = (sum+1)/Float_t(summin);
2112 void AliTPCtrackerMI::SignShared(AliTPCseed * s1, AliTPCseed * s2)
2116 Float_t thetaCut = 0.2;//+10.*TMath::Sqrt(s1->GetSigmaTglZ()+ s2->GetSigmaTglZ());
2117 if (TMath::Abs(s1->GetTgl()-s2->GetTgl())>thetaCut) return;
2118 Float_t minCl = TMath::Min(s1->GetNumberOfClusters(),s2->GetNumberOfClusters());
2119 Int_t cutN0 = TMath::Max(5,TMath::Nint(0.1*minCl));
2124 //Int_t firstpoint = TMath::Max(s1->GetFirstPoint(),s2->GetFirstPoint());
2125 //Int_t lastpoint = TMath::Min(s1->GetLastPoint(),s2->GetLastPoint());
2126 Int_t firstpoint = 0;
2127 Int_t lastpoint = 160;
2129 // if (firstpoint>=lastpoint-5) return;;
2131 for (Int_t i=firstpoint;i<lastpoint;i++){
2132 // if ( (s1->GetClusterIndex2(i)&0xFFFF8FFF)==(s2->GetClusterIndex2(i)&0xFFFF8FFF) && s1->GetClusterIndex2(i)>0) {
2133 if ( (s1->GetClusterIndex2(i))==(s2->GetClusterIndex2(i)) && s1->GetClusterIndex2(i)>0) {
2137 if (sumshared>cutN0){
2140 for (Int_t i=firstpoint;i<lastpoint;i++){
2141 // if ( (s1->GetClusterIndex2(i)&0xFFFF8FFF)==(s2->GetClusterIndex2(i)&0xFFFF8FFF) && s1->GetClusterIndex2(i)>0) {
2142 if ( (s1->GetClusterIndex2(i))==(s2->GetClusterIndex2(i)) && s1->GetClusterIndex2(i)>0) {
2143 AliTPCTrackerPoint *p1 = s1->GetTrackPoint(i);
2144 AliTPCTrackerPoint *p2 = s2->GetTrackPoint(i);;
2145 if (s1->IsActive()&&s2->IsActive()){
2146 p1->SetShared(kTRUE);
2147 p2->SetShared(kTRUE);
2153 if (sumshared>cutN0){
2154 for (Int_t i=0;i<4;i++){
2155 if (s1->GetOverlapLabel(3*i)==0){
2156 s1->SetOverlapLabel(3*i, s2->GetLabel());
2157 s1->SetOverlapLabel(3*i+1,sumshared);
2158 s1->SetOverlapLabel(3*i+2,s2->GetUniqueID());
2162 for (Int_t i=0;i<4;i++){
2163 if (s2->GetOverlapLabel(3*i)==0){
2164 s2->SetOverlapLabel(3*i, s1->GetLabel());
2165 s2->SetOverlapLabel(3*i+1,sumshared);
2166 s2->SetOverlapLabel(3*i+2,s1->GetUniqueID());
2173 void AliTPCtrackerMI::SignShared(TObjArray * arr)
2176 //sort trackss according sectors
2178 for (Int_t i=0; i<arr->GetEntriesFast(); i++) {
2179 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2181 //if (pt) RotateToLocal(pt);
2185 arr->Sort(); // sorting according relative sectors
2186 arr->Expand(arr->GetEntries());
2189 Int_t nseed=arr->GetEntriesFast();
2190 for (Int_t i=0; i<nseed; i++) {
2191 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2193 for (Int_t j=0;j<12;j++){
2194 pt->SetOverlapLabel(j,0);
2197 for (Int_t i=0; i<nseed; i++) {
2198 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2200 if (pt->GetRemoval()>10) continue;
2201 for (Int_t j=i+1; j<nseed; j++){
2202 AliTPCseed *pt2=(AliTPCseed*)arr->UncheckedAt(j);
2203 if (TMath::Abs(pt->GetRelativeSector()-pt2->GetRelativeSector())>1) continue;
2205 if (pt2->GetRemoval()<=10) {
2206 //if ( TMath::Abs(pt->GetRelativeSector()-pt2->GetRelativeSector())>0) break;
2214 void AliTPCtrackerMI::SortTracks(TObjArray * arr, Int_t mode) const
2217 //sort tracks in array according mode criteria
2218 Int_t nseed = arr->GetEntriesFast();
2219 for (Int_t i=0; i<nseed; i++) {
2220 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2231 void AliTPCtrackerMI::RemoveUsed2(TObjArray * arr, Float_t factor1, Float_t factor2, Int_t minimal)
2234 // Loop over all tracks and remove overlaped tracks (with lower quality)
2236 // 1. Unsign clusters
2237 // 2. Sort tracks according quality
2238 // Quality is defined by the number of cluster between first and last points
2240 // 3. Loop over tracks - decreasing quality order
2241 // a.) remove - If the fraction of shared cluster less than factor (1- n or 2)
2242 // b.) remove - If the minimal number of clusters less than minimal and not ITS
2243 // c.) if track accepted - sign clusters
2245 //Called in - AliTPCtrackerMI::Clusters2Tracks()
2246 // - AliTPCtrackerMI::PropagateBack()
2247 // - AliTPCtrackerMI::RefitInward()
2250 // factor1 - factor for constrained
2251 // factor2 - for non constrained tracks
2252 // if (Float_t(shared+1)/Float_t(found+1)>factor) - DELETE
2256 Int_t nseed = arr->GetEntriesFast();
2257 Float_t * quality = new Float_t[nseed];
2258 Int_t * indexes = new Int_t[nseed];
2262 for (Int_t i=0; i<nseed; i++) {
2263 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2268 pt->UpdatePoints(); //select first last max dens points
2269 Float_t * points = pt->GetPoints();
2270 if (points[3]<0.8) quality[i] =-1;
2271 quality[i] = (points[2]-points[0])+pt->GetNumberOfClusters();
2272 //prefer high momenta tracks if overlaps
2273 quality[i] *= TMath::Sqrt(TMath::Abs(pt->Pt())+0.5);
2275 TMath::Sort(nseed,quality,indexes);
2278 for (Int_t itrack=0; itrack<nseed; itrack++) {
2279 Int_t trackindex = indexes[itrack];
2280 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(trackindex);
2283 if (quality[trackindex]<0){
2284 delete arr->RemoveAt(trackindex);
2289 Int_t first = Int_t(pt->GetPoints()[0]);
2290 Int_t last = Int_t(pt->GetPoints()[2]);
2291 Double_t factor = (pt->GetBConstrain()) ? factor1: factor2;
2293 Int_t found,foundable,shared;
2294 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
2295 // pt->GetClusterStatistic(0,160, found, foundable,shared,kFALSE);
2296 Bool_t itsgold =kFALSE;
2299 if (pt->GetESD()->GetITSclusters(dummy)>4) itsgold= kTRUE;
2303 if (Float_t(shared+1)/Float_t(found+1)>factor){
2304 if (pt->GetKinkIndexes()[0]!=0) continue; //don't remove tracks - part of the kinks
2305 if( AliTPCReconstructor::StreamLevel()>3){
2306 TTreeSRedirector &cstream = *fDebugStreamer;
2307 cstream<<"RemoveUsed"<<
2308 "iter="<<fIteration<<
2312 delete arr->RemoveAt(trackindex);
2315 if (pt->GetNumberOfClusters()<50&&(found-0.5*shared)<minimal){ //remove short tracks
2316 if (pt->GetKinkIndexes()[0]!=0) continue; //don't remove tracks - part of the kinks
2317 if( AliTPCReconstructor::StreamLevel()>3){
2318 TTreeSRedirector &cstream = *fDebugStreamer;
2319 cstream<<"RemoveShort"<<
2320 "iter="<<fIteration<<
2324 delete arr->RemoveAt(trackindex);
2330 //if (sharedfactor>0.4) continue;
2331 if (pt->GetKinkIndexes()[0]>0) continue;
2332 //Remove tracks with undefined properties - seems
2333 if (pt->GetSigmaY2()<kAlmost0) continue; // ? what is the origin ?
2335 for (Int_t i=first; i<last; i++) {
2336 Int_t index=pt->GetClusterIndex2(i);
2337 // if (index<0 || index&0x8000 ) continue;
2338 if (index<0 || index&0x8000 ) continue;
2339 AliTPCclusterMI *c= pt->GetClusterPointer(i);
2346 Info("RemoveUsed","\n*****\nNumber of good tracks after shared removal\t%d\n",fNtracks);
2352 void AliTPCtrackerMI::DumpClusters(Int_t iter, TObjArray *trackArray)
2355 // Dump clusters after reco
2356 // signed and unsigned cluster can be visualized
2357 // 1. Unsign all cluster
2358 // 2. Sign all used clusters
2361 Int_t nseed = trackArray->GetEntries();
2362 for (Int_t i=0; i<nseed; i++){
2363 AliTPCseed *pt=(AliTPCseed*)trackArray->UncheckedAt(i);
2367 Bool_t isKink=pt->GetKinkIndex(0)!=0;
2368 for (Int_t j=0; j<160; ++j) {
2369 Int_t index=pt->GetClusterIndex2(j);
2370 if (index<0) continue;
2371 AliTPCclusterMI *c= pt->GetClusterPointer(j);
2373 if (isKink) c->Use(100); // kink
2374 c->Use(10); // by default usage 10
2379 for (Int_t sec=0;sec<fkNIS;sec++){
2380 for (Int_t row=0;row<fInnerSec->GetNRows();row++){
2381 AliTPCclusterMI *cl = fInnerSec[sec][row].GetClusters1();
2382 for (Int_t icl =0;icl< fInnerSec[sec][row].GetN1();icl++){
2383 Float_t gx[3]; cl[icl].GetGlobalXYZ(gx);
2384 (*fDebugStreamer)<<"clDump"<<
2392 cl = fInnerSec[sec][row].GetClusters2();
2393 for (Int_t icl =0;icl< fInnerSec[sec][row].GetN2();icl++){
2394 Float_t gx[3]; cl[icl].GetGlobalXYZ(gx);
2395 (*fDebugStreamer)<<"clDump"<<
2406 for (Int_t sec=0;sec<fkNOS;sec++){
2407 for (Int_t row=0;row<fOuterSec->GetNRows();row++){
2408 AliTPCclusterMI *cl = fOuterSec[sec][row].GetClusters1();
2409 for (Int_t icl =0;icl< fOuterSec[sec][row].GetN1();icl++){
2410 Float_t gx[3]; cl[icl].GetGlobalXYZ(gx);
2411 (*fDebugStreamer)<<"clDump"<<
2419 cl = fOuterSec[sec][row].GetClusters2();
2420 for (Int_t icl =0;icl< fOuterSec[sec][row].GetN2();icl++){
2421 Float_t gx[3]; cl[icl].GetGlobalXYZ(gx);
2422 (*fDebugStreamer)<<"clDump"<<
2434 void AliTPCtrackerMI::UnsignClusters()
2437 // loop over all clusters and unsign them
2440 for (Int_t sec=0;sec<fkNIS;sec++){
2441 for (Int_t row=0;row<fInnerSec->GetNRows();row++){
2442 AliTPCclusterMI *cl = fInnerSec[sec][row].GetClusters1();
2443 for (Int_t icl =0;icl< fInnerSec[sec][row].GetN1();icl++)
2444 // if (cl[icl].IsUsed(10))
2446 cl = fInnerSec[sec][row].GetClusters2();
2447 for (Int_t icl =0;icl< fInnerSec[sec][row].GetN2();icl++)
2448 //if (cl[icl].IsUsed(10))
2453 for (Int_t sec=0;sec<fkNOS;sec++){
2454 for (Int_t row=0;row<fOuterSec->GetNRows();row++){
2455 AliTPCclusterMI *cl = fOuterSec[sec][row].GetClusters1();
2456 for (Int_t icl =0;icl< fOuterSec[sec][row].GetN1();icl++)
2457 //if (cl[icl].IsUsed(10))
2459 cl = fOuterSec[sec][row].GetClusters2();
2460 for (Int_t icl =0;icl< fOuterSec[sec][row].GetN2();icl++)
2461 //if (cl[icl].IsUsed(10))
2470 void AliTPCtrackerMI::SignClusters(const TObjArray * arr, Float_t fnumber, Float_t fdensity)
2473 //sign clusters to be "used"
2475 // snumber and sdensity sign number of sigmas - bellow mean value to be accepted
2476 // loop over "primaries"
2490 Int_t nseed = arr->GetEntriesFast();
2491 for (Int_t i=0; i<nseed; i++) {
2492 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2496 if (!(pt->IsActive())) continue;
2497 Float_t dens = pt->GetNumberOfClusters()/Float_t(pt->GetNFoundable());
2498 if ( (dens>0.7) && (pt->GetNumberOfClusters()>70)){
2500 sumdens2+= dens*dens;
2501 sumn += pt->GetNumberOfClusters();
2502 sumn2 += pt->GetNumberOfClusters()*pt->GetNumberOfClusters();
2503 Float_t chi2 = pt->GetChi2()/pt->GetNumberOfClusters();
2506 sumchi2 +=chi2*chi2;
2511 Float_t mdensity = 0.9;
2512 Float_t meann = 130;
2513 Float_t meanchi = 1;
2514 Float_t sdensity = 0.1;
2515 Float_t smeann = 10;
2516 Float_t smeanchi =0.4;
2520 mdensity = sumdens/sum;
2522 meanchi = sumchi/sum;
2524 sdensity = sumdens2/sum-mdensity*mdensity;
2526 sdensity = TMath::Sqrt(sdensity);
2530 smeann = sumn2/sum-meann*meann;
2532 smeann = TMath::Sqrt(smeann);
2536 smeanchi = sumchi2/sum - meanchi*meanchi;
2538 smeanchi = TMath::Sqrt(smeanchi);
2544 //REMOVE SHORT DELTAS or tracks going out of sensitive volume of TPC
2546 for (Int_t i=0; i<nseed; i++) {
2547 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2551 if (pt->GetBSigned()) continue;
2552 if (pt->GetBConstrain()) continue;
2553 //if (!(pt->IsActive())) continue;
2555 Int_t found,foundable,shared;
2556 pt->GetClusterStatistic(0,160,found, foundable,shared);
2557 if (shared/float(found)>0.3) {
2558 if (shared/float(found)>0.9 ){
2559 //delete arr->RemoveAt(i);
2564 Bool_t isok =kFALSE;
2565 if ( (pt->GetNShared()/pt->GetNumberOfClusters()<0.5) &&pt->GetNumberOfClusters()>60)
2567 if ((TMath::Abs(1/pt->GetC())<100.) && (pt->GetNShared()/pt->GetNumberOfClusters()<0.7))
2569 if (TMath::Abs(pt->GetZ()/pt->GetX())>1.1)
2571 if ( (TMath::Abs(pt->GetSnp()>0.7) && pt->GetD(0,0)>60.))
2575 for (Int_t j=0; j<160; ++j) {
2576 Int_t index=pt->GetClusterIndex2(j);
2577 if (index<0) continue;
2578 AliTPCclusterMI *c= pt->GetClusterPointer(j);
2580 //if (!(c->IsUsed(10))) c->Use();
2587 Double_t maxchi = meanchi+2.*smeanchi;
2589 for (Int_t i=0; i<nseed; i++) {
2590 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2594 //if (!(pt->IsActive())) continue;
2595 if (pt->GetBSigned()) continue;
2596 Double_t chi = pt->GetChi2()/pt->GetNumberOfClusters();
2597 if (chi>maxchi) continue;
2600 Float_t dens = pt->GetNumberOfClusters()/Float_t(pt->GetNFoundable());
2602 //sign only tracks with enoug big density at the beginning
2604 if ((pt->GetDensityFirst(40)<0.75) && pt->GetNumberOfClusters()<meann) continue;
2607 Double_t mindens = TMath::Max(double(mdensity-sdensity*fdensity*bfactor),0.65);
2608 Double_t minn = TMath::Max(Int_t(meann-fnumber*smeann*bfactor),50);
2610 // if (pt->fBConstrain) mindens = TMath::Max(mdensity-sdensity*fdensity*bfactor,0.65);
2611 if ( (pt->GetRemoval()==10) && (pt->GetSnp()>0.8)&&(dens>mindens))
2614 if ((dens>mindens && pt->GetNumberOfClusters()>minn) && chi<maxchi ){
2615 //Int_t noc=pt->GetNumberOfClusters();
2616 pt->SetBSigned(kTRUE);
2617 for (Int_t j=0; j<160; ++j) {
2619 Int_t index=pt->GetClusterIndex2(j);
2620 if (index<0) continue;
2621 AliTPCclusterMI *c= pt->GetClusterPointer(j);
2623 // if (!(c->IsUsed(10))) c->Use();
2628 // gLastCheck = nseed;
2636 void AliTPCtrackerMI::StopNotActive(const TObjArray * arr, Int_t row0, Float_t th0, Float_t th1, Float_t th2) const
2638 // stop not active tracks
2639 // take th1 as threshold for number of founded to number of foundable on last 10 active rows
2640 // take th2 as threshold for number of founded to number of foundable on last 20 active rows
2641 Int_t nseed = arr->GetEntriesFast();
2643 for (Int_t i=0; i<nseed; i++) {
2644 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2648 if (!(pt->IsActive())) continue;
2649 StopNotActive(pt,row0,th0, th1,th2);
2655 void AliTPCtrackerMI::StopNotActive(AliTPCseed * seed, Int_t row0, Float_t th0, Float_t th1,
2658 // stop not active tracks
2659 // take th1 as threshold for number of founded to number of foundable on last 10 active rows
2660 // take th2 as threshold for number of founded to number of foundable on last 20 active rows
2663 Int_t foundable = 0;
2664 Int_t maxindex = seed->GetLastPoint(); //last foundable row
2665 if (seed->GetNFoundable()*th0 > seed->GetNumberOfClusters()) {
2666 seed->Desactivate(10) ;
2670 for (Int_t i=row0; i<maxindex; i++){
2671 Int_t index = seed->GetClusterIndex2(i);
2672 if (index!=-1) foundable++;
2674 if (foundable<=30) sumgood1++;
2675 if (foundable<=50) {
2682 if (foundable>=30.){
2683 if (sumgood1<(th1*30.)) seed->Desactivate(10);
2686 if (sumgood2<(th2*50.)) seed->Desactivate(10);
2690 Int_t AliTPCtrackerMI::RefitInward(AliESDEvent *event)
2693 // back propagation of ESD tracks
2696 if (!event) return 0;
2697 const Int_t kMaxFriendTracks=2000;
2699 // extract correction object for multiplicity dependence of dEdx
2700 TObjArray * gainCalibArray = AliTPCcalibDB::Instance()->GetTimeGainSplinesRun(event->GetRunNumber());
2702 AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ;
2704 AliFatal("Tranformations not in RefitInward");
2707 transform->SetCurrentRecoParam((AliTPCRecoParam*)AliTPCReconstructor::GetRecoParam());
2708 const AliTPCRecoParam * recoParam = AliTPCcalibDB::Instance()->GetTransform()->GetCurrentRecoParam();
2709 Int_t nContribut = event->GetNumberOfTracks();
2710 TGraphErrors * graphMultDependenceDeDx = 0x0;
2711 if (recoParam && recoParam->GetUseMultiplicityCorrectionDedx() && gainCalibArray) {
2712 if (recoParam->GetUseTotCharge()) {
2713 graphMultDependenceDeDx = (TGraphErrors *) gainCalibArray->FindObject("TGRAPHERRORS_MEANQTOT_MULTIPLICITYDEPENDENCE_BEAM_ALL");
2715 graphMultDependenceDeDx = (TGraphErrors *) gainCalibArray->FindObject("TGRAPHERRORS_MEANQMAX_MULTIPLICITYDEPENDENCE_BEAM_ALL");
2721 //PrepareForProlongation(fSeeds,1);
2722 PropagateForward2(fSeeds);
2723 RemoveUsed2(fSeeds,0.4,0.4,20);
2725 TObjArray arraySeed(fSeeds->GetEntries());
2726 for (Int_t i=0;i<fSeeds->GetEntries();i++) {
2727 arraySeed.AddAt(fSeeds->At(i),i);
2729 SignShared(&arraySeed);
2730 // FindCurling(fSeeds, event,2); // find multi found tracks
2731 FindSplitted(fSeeds, event,2); // find multi found tracks
2732 if (AliTPCReconstructor::StreamLevel()>5) FindMultiMC(fSeeds, fEvent,2); // find multi found tracks
2735 Int_t nseed = fSeeds->GetEntriesFast();
2736 for (Int_t i=0;i<nseed;i++){
2737 AliTPCseed * seed = (AliTPCseed*) fSeeds->UncheckedAt(i);
2738 if (!seed) continue;
2739 if (seed->GetKinkIndex(0)>0) UpdateKinkQualityD(seed); // update quality informations for kinks
2740 AliESDtrack *esd=event->GetTrack(i);
2742 if (seed->GetNumberOfClusters()<60 && seed->GetNumberOfClusters()<(esd->GetTPCclusters(0) -5)*0.8){
2743 AliExternalTrackParam paramIn;
2744 AliExternalTrackParam paramOut;
2745 Int_t ncl = seed->RefitTrack(seed,¶mIn,¶mOut);
2746 if (AliTPCReconstructor::StreamLevel()>2) {
2747 (*fDebugStreamer)<<"RecoverIn"<<
2751 "pout.="<<¶mOut<<
2756 seed->Set(paramIn.GetX(),paramIn.GetAlpha(),paramIn.GetParameter(),paramIn.GetCovariance());
2757 seed->SetNumberOfClusters(ncl);
2761 seed->PropagateTo(fkParam->GetInnerRadiusLow());
2762 seed->UpdatePoints();
2763 AddCovariance(seed);
2764 MakeESDBitmaps(seed, esd);
2765 seed->CookdEdx(0.02,0.6);
2766 CookLabel(seed,0.1); //For comparison only
2768 if (AliTPCReconstructor::StreamLevel()>1 && seed!=0) {
2769 TTreeSRedirector &cstream = *fDebugStreamer;
2776 if (seed->GetNumberOfClusters()>15){
2777 esd->UpdateTrackParams(seed,AliESDtrack::kTPCrefit);
2778 esd->SetTPCPoints(seed->GetPoints());
2779 esd->SetTPCPointsF(seed->GetNFoundable());
2780 Int_t ndedx = seed->GetNCDEDX(0);
2781 Float_t sdedx = seed->GetSDEDX(0);
2782 Float_t dedx = seed->GetdEdx();
2783 // apply mutliplicity dependent dEdx correction if available
2784 if (graphMultDependenceDeDx) {
2785 Double_t corrGain = AliTPCcalibDButil::EvalGraphConst(graphMultDependenceDeDx, nContribut);
2786 dedx += (1 - corrGain)*50.; // MIP is normalized to 50
2788 esd->SetTPCsignal(dedx, sdedx, ndedx);
2790 // fill new dEdx information
2792 Double32_t signal[4];
2796 for(Int_t iarr=0;iarr<3;iarr++) {
2797 signal[iarr] = seed->GetDEDXregion(iarr+1);
2798 ncl[iarr] = seed->GetNCDEDX(iarr+1);
2799 nrows[iarr] = seed->GetNCDEDXInclThres(iarr+1);
2801 signal[3] = seed->GetDEDXregion(4);
2803 AliTPCdEdxInfo * infoTpcPid = new AliTPCdEdxInfo();
2804 infoTpcPid->SetTPCSignalRegionInfo(signal, ncl, nrows);
2805 esd->SetTPCdEdxInfo(infoTpcPid);
2807 // add seed to the esd track in Calib level
2809 Bool_t storeFriend = gRandom->Rndm()<(kMaxFriendTracks)/Float_t(nseed);
2810 if (AliTPCReconstructor::StreamLevel()>0 &&storeFriend){
2811 AliTPCseed * seedCopy = new AliTPCseed(*seed, kTRUE);
2812 esd->AddCalibObject(seedCopy);
2817 //printf("problem\n");
2820 //FindKinks(fSeeds,event);
2821 if (AliTPCReconstructor::StreamLevel()>3) DumpClusters(2,fSeeds);
2822 Info("RefitInward","Number of refitted tracks %d",ntracks);
2827 Int_t AliTPCtrackerMI::PropagateBack(AliESDEvent *event)
2830 // back propagation of ESD tracks
2832 if (!event) return 0;
2836 PropagateBack(fSeeds);
2837 RemoveUsed2(fSeeds,0.4,0.4,20);
2838 //FindCurling(fSeeds, fEvent,1);
2839 FindSplitted(fSeeds, event,1); // find multi found tracks
2840 if (AliTPCReconstructor::StreamLevel()>5) FindMultiMC(fSeeds, fEvent,1); // find multi found tracks
2843 Int_t nseed = fSeeds->GetEntriesFast();
2845 for (Int_t i=0;i<nseed;i++){
2846 AliTPCseed * seed = (AliTPCseed*) fSeeds->UncheckedAt(i);
2847 if (!seed) continue;
2848 if (seed->GetKinkIndex(0)<0) UpdateKinkQualityM(seed); // update quality informations for kinks
2849 seed->UpdatePoints();
2850 AddCovariance(seed);
2851 AliESDtrack *esd=event->GetTrack(i);
2852 if (!esd) continue; //never happen
2853 if (seed->GetNumberOfClusters()<60 && seed->GetNumberOfClusters()<(esd->GetTPCclusters(0) -5)*0.8){
2854 AliExternalTrackParam paramIn;
2855 AliExternalTrackParam paramOut;
2856 Int_t ncl = seed->RefitTrack(seed,¶mIn,¶mOut);
2857 if (AliTPCReconstructor::StreamLevel()>2) {
2858 (*fDebugStreamer)<<"RecoverBack"<<
2862 "pout.="<<¶mOut<<
2867 seed->Set(paramOut.GetX(),paramOut.GetAlpha(),paramOut.GetParameter(),paramOut.GetCovariance());
2868 seed->SetNumberOfClusters(ncl);
2871 seed->CookdEdx(0.02,0.6);
2872 CookLabel(seed,0.1); //For comparison only
2873 if (seed->GetNumberOfClusters()>15){
2874 esd->UpdateTrackParams(seed,AliESDtrack::kTPCout);
2875 esd->SetTPCPoints(seed->GetPoints());
2876 esd->SetTPCPointsF(seed->GetNFoundable());
2877 Int_t ndedx = seed->GetNCDEDX(0);
2878 Float_t sdedx = seed->GetSDEDX(0);
2879 Float_t dedx = seed->GetdEdx();
2880 esd->SetTPCsignal(dedx, sdedx, ndedx);
2882 Int_t eventnumber = event->GetEventNumberInFile();// patch 28 fev 06
2883 // This is most likely NOT the event number you'd like to use. It has nothing to do with the 'real' event number
2884 if (AliTPCReconstructor::StreamLevel()>1 && esd) {
2885 (*fDebugStreamer)<<"Cback"<<
2888 "EventNrInFile="<<eventnumber<<
2893 if (AliTPCReconstructor::StreamLevel()>3) DumpClusters(1,fSeeds);
2894 //FindKinks(fSeeds,event);
2895 Info("PropagateBack","Number of back propagated tracks %d",ntracks);
2902 void AliTPCtrackerMI::DeleteSeeds()
2907 Int_t nseed = fSeeds->GetEntriesFast();
2908 for (Int_t i=0;i<nseed;i++){
2909 AliTPCseed * seed = (AliTPCseed*)fSeeds->At(i);
2910 if (seed) delete fSeeds->RemoveAt(i);
2917 void AliTPCtrackerMI::ReadSeeds(const AliESDEvent *const event, Int_t direction)
2920 //read seeds from the event
2922 Int_t nentr=event->GetNumberOfTracks();
2924 Info("ReadSeeds", "Number of ESD tracks: %d\n", nentr);
2929 fSeeds = new TObjArray(nentr);
2933 for (Int_t i=0; i<nentr; i++) {
2934 AliESDtrack *esd=event->GetTrack(i);
2935 ULong_t status=esd->GetStatus();
2936 if (!(status&AliESDtrack::kTPCin)) continue;
2937 AliTPCtrack t(*esd);
2938 t.SetNumberOfClusters(0);
2939 // AliTPCseed *seed = new AliTPCseed(t,t.GetAlpha());
2940 AliTPCseed *seed = new AliTPCseed(t/*,t.GetAlpha()*/);
2941 seed->SetUniqueID(esd->GetID());
2942 AddCovariance(seed); //add systematic ucertainty
2943 for (Int_t ikink=0;ikink<3;ikink++) {
2944 Int_t index = esd->GetKinkIndex(ikink);
2945 seed->GetKinkIndexes()[ikink] = index;
2946 if (index==0) continue;
2947 index = TMath::Abs(index);
2948 AliESDkink * kink = fEvent->GetKink(index-1);
2949 if (kink&&esd->GetKinkIndex(ikink)<0){
2950 if ((status & AliESDtrack::kTRDrefit) != 0) kink->SetStatus(1,2);
2951 if ((status & AliESDtrack::kITSout) != 0) kink->SetStatus(1,0);
2953 if (kink&&esd->GetKinkIndex(ikink)>0){
2954 if ((status & AliESDtrack::kTRDrefit) != 0) kink->SetStatus(1,6);
2955 if ((status & AliESDtrack::kITSout) != 0) kink->SetStatus(1,4);
2959 if (((status&AliESDtrack::kITSout)==0)&&(direction==1)) seed->ResetCovariance(10.);
2960 if ( direction ==2 &&(status & AliESDtrack::kTRDrefit) == 0 ) seed->ResetCovariance(10.);
2961 //if ( direction ==2 && ((status & AliESDtrack::kTPCout) == 0) ) {
2962 // fSeeds->AddAt(0,i);
2966 if ( direction ==2 &&(status & AliESDtrack::kTRDrefit) > 0 ) {
2967 Double_t par0[5],par1[5],alpha,x;
2968 esd->GetInnerExternalParameters(alpha,x,par0);
2969 esd->GetExternalParameters(x,par1);
2970 Double_t delta1 = TMath::Abs(par0[4]-par1[4])/(0.000000001+TMath::Abs(par0[4]+par1[4]));
2971 Double_t delta2 = TMath::Abs(par0[3]-par1[3]);
2973 if (esd->GetTRDncls()>0) trdchi2 = esd->GetTRDchi2()/esd->GetTRDncls();
2974 //reset covariance if suspicious
2975 if ( (delta1>0.1) || (delta2>0.006) ||trdchi2>7.)
2976 seed->ResetCovariance(10.);
2981 // rotate to the local coordinate system
2983 fSectors=fInnerSec; fN=fkNIS;
2984 Double_t alpha=seed->GetAlpha();
2985 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
2986 if (alpha < 0. ) alpha += 2.*TMath::Pi();
2987 Int_t ns=Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN;
2988 alpha =ns*fSectors->GetAlpha() + fSectors->GetAlphaShift();
2989 alpha-=seed->GetAlpha();
2990 if (alpha<-TMath::Pi()) alpha += 2*TMath::Pi();
2991 if (alpha>=TMath::Pi()) alpha -= 2*TMath::Pi();
2992 if (TMath::Abs(alpha) > 0.001) { // This should not happen normally
2993 AliWarning(Form("Rotating track over %f",alpha));
2994 if (!seed->Rotate(alpha)) {
3001 if (esd->GetKinkIndex(0)<=0){
3002 for (Int_t irow=0;irow<160;irow++){
3003 Int_t index = seed->GetClusterIndex2(irow);
3006 AliTPCclusterMI * cl = GetClusterMI(index);
3007 seed->SetClusterPointer(irow,cl);
3009 if ((index & 0x8000)==0){
3010 cl->Use(10); // accepted cluster
3012 cl->Use(6); // close cluster not accepted
3015 Info("ReadSeeds","Not found cluster");
3020 fSeeds->AddAt(seed,i);
3026 //_____________________________________________________________________________
3027 void AliTPCtrackerMI::MakeSeeds3(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, Float_t cuts[4],
3028 Float_t deltay, Int_t ddsec) {
3029 //-----------------------------------------------------------------
3030 // This function creates track seeds.
3031 // SEEDING WITH VERTEX CONSTRAIN
3032 //-----------------------------------------------------------------
3033 // cuts[0] - fP4 cut
3034 // cuts[1] - tan(phi) cut
3035 // cuts[2] - zvertex cut
3036 // cuts[3] - fP3 cut
3044 Double_t x[5], c[15];
3045 // Int_t di = i1-i2;
3047 AliTPCseed * seed = new AliTPCseed();
3048 Double_t alpha=fSectors->GetAlpha(), shift=fSectors->GetAlphaShift();
3049 Double_t cs=cos(alpha), sn=sin(alpha);
3051 // Double_t x1 =fOuterSec->GetX(i1);
3052 //Double_t xx2=fOuterSec->GetX(i2);
3054 Double_t x1 =GetXrow(i1);
3055 Double_t xx2=GetXrow(i2);
3057 Double_t x3=GetX(), y3=GetY(), z3=GetZ();
3059 Int_t imiddle = (i2+i1)/2; //middle pad row index
3060 Double_t xm = GetXrow(imiddle); // radius of middle pad-row
3061 const AliTPCtrackerRow& krm=GetRow(sec,imiddle); //middle pad -row
3065 const AliTPCtrackerRow& kr1=GetRow(ns,i1);
3066 Double_t ymax = GetMaxY(i1)-kr1.GetDeadZone()-1.5;
3067 Double_t ymaxm = GetMaxY(imiddle)-kr1.GetDeadZone()-1.5;
3070 // change cut on curvature if it can't reach this layer
3071 // maximal curvature set to reach it
3072 Double_t dvertexmax = TMath::Sqrt((x1-x3)*(x1-x3)+(ymax+5-y3)*(ymax+5-y3));
3073 if (dvertexmax*0.5*cuts[0]>0.85){
3074 cuts[0] = 0.85/(dvertexmax*0.5+1.);
3076 Double_t r2min = 1/(cuts[0]*cuts[0]); //minimal square of radius given by cut
3079 if (deltay>0) ddsec = 0;
3080 // loop over clusters
3081 for (Int_t is=0; is < kr1; is++) {
3083 if (kr1[is]->IsUsed(10)) continue;
3084 Double_t y1=kr1[is]->GetY(), z1=kr1[is]->GetZ();
3085 //if (TMath::Abs(y1)>ymax) continue;
3087 if (deltay>0 && TMath::Abs(ymax-TMath::Abs(y1))> deltay ) continue; // seed only at the edge
3089 // find possible directions
3090 Float_t anglez = (z1-z3)/(x1-x3);
3091 Float_t extraz = z1 - anglez*(x1-xx2); // extrapolated z
3094 //find rotation angles relative to line given by vertex and point 1
3095 Double_t dvertex2 = (x1-x3)*(x1-x3)+(y1-y3)*(y1-y3);
3096 Double_t dvertex = TMath::Sqrt(dvertex2);
3097 Double_t angle13 = TMath::ATan((y1-y3)/(x1-x3));
3098 Double_t cs13 = cos(-angle13), sn13 = sin(-angle13);
3101 // loop over 2 sectors
3107 Double_t dddz1=0; // direction of delta inclination in z axis
3114 for (Int_t dsec = dsec1; dsec<=dsec2;dsec++){
3115 Int_t sec2 = sec + dsec;
3117 // AliTPCtrackerRow& kr2 = fOuterSec[(sec2+fkNOS)%fkNOS][i2];
3118 //AliTPCtrackerRow& kr2m = fOuterSec[(sec2+fkNOS)%fkNOS][imiddle];
3119 AliTPCtrackerRow& kr2 = GetRow((sec2+fkNOS)%fkNOS,i2);
3120 AliTPCtrackerRow& kr2m = GetRow((sec2+fkNOS)%fkNOS,imiddle);
3121 Int_t index1 = TMath::Max(kr2.Find(extraz-0.6-dddz1*TMath::Abs(z1)*0.05)-1,0);
3122 Int_t index2 = TMath::Min(kr2.Find(extraz+0.6+dddz2*TMath::Abs(z1)*0.05)+1,kr2);
3124 // rotation angles to p1-p3
3125 Double_t cs13r = cos(-angle13+dsec*alpha)/dvertex, sn13r = sin(-angle13+dsec*alpha)/dvertex;
3126 Double_t x2, y2, z2;
3128 // Double_t dymax = maxangle*TMath::Abs(x1-xx2);
3131 Double_t dxx0 = (xx2-x3)*cs13r;
3132 Double_t dyy0 = (xx2-x3)*sn13r;
3133 for (Int_t js=index1; js < index2; js++) {
3134 const AliTPCclusterMI *kcl = kr2[js];
3135 if (kcl->IsUsed(10)) continue;
3137 //calcutate parameters
3139 Double_t yy0 = dyy0 +(kcl->GetY()-y3)*cs13r;
3141 if (TMath::Abs(yy0)<0.000001) continue;
3142 Double_t xx0 = dxx0 -(kcl->GetY()-y3)*sn13r;
3143 Double_t y0 = 0.5*(xx0*xx0+yy0*yy0-xx0)/yy0;
3144 Double_t r02 = (0.25+y0*y0)*dvertex2;
3145 //curvature (radius) cut
3146 if (r02<r2min) continue;
3150 Double_t c0 = 1/TMath::Sqrt(r02);
3154 //Double_t dfi0 = 2.*TMath::ASin(dvertex*c0*0.5);
3155 //Double_t dfi1 = 2.*TMath::ASin(TMath::Sqrt(yy0*yy0+(1-xx0)*(1-xx0))*dvertex*c0*0.5);
3156 Double_t dfi0 = 2.*AliTPCFastMath::FastAsin(dvertex*c0*0.5);
3157 Double_t dfi1 = 2.*AliTPCFastMath::FastAsin(TMath::Sqrt(yy0*yy0+(1-xx0)*(1-xx0))*dvertex*c0*0.5);
3160 Double_t z0 = kcl->GetZ();
3161 Double_t zzzz2 = z1-(z1-z3)*dfi1/dfi0;
3162 if (TMath::Abs(zzzz2-z0)>0.5) continue;
3165 Double_t dip = (z1-z0)*c0/dfi1;
3166 Double_t x0 = (0.5*cs13+y0*sn13)*dvertex*c0;
3177 x2= xx2*cs-y2*sn*dsec;
3178 y2=+xx2*sn*dsec+y2*cs;
3188 // do we have cluster at the middle ?
3190 GetProlongation(x1,xm,x,ym,zm);
3192 AliTPCclusterMI * cm=0;
3193 if (TMath::Abs(ym)-ymaxm<0){
3194 cm = krm.FindNearest2(ym,zm,1.0,0.6,dummy);
3195 if ((!cm) || (cm->IsUsed(10))) {
3200 // rotate y1 to system 0
3201 // get state vector in rotated system
3202 Double_t yr1 = (-0.5*sn13+y0*cs13)*dvertex*c0;
3203 Double_t xr2 = x0*cs+yr1*sn*dsec;
3204 Double_t xr[5]={kcl->GetY(),kcl->GetZ(), xr2, dip, c0};
3206 GetProlongation(xx2,xm,xr,ym,zm);
3207 if (TMath::Abs(ym)-ymaxm<0){
3208 cm = kr2m.FindNearest2(ym,zm,1.0,0.6,dummy);
3209 if ((!cm) || (cm->IsUsed(10))) {
3219 dym = ym - cm->GetY();
3220 dzm = zm - cm->GetZ();
3227 Double_t sy1=kr1[is]->GetSigmaY2()*2., sz1=kr1[is]->GetSigmaZ2()*2.;
3228 Double_t sy2=kcl->GetSigmaY2()*2., sz2=kcl->GetSigmaZ2()*2.;
3229 //Double_t sy3=400*3./12., sy=0.1, sz=0.1;
3230 Double_t sy3=25000*x[4]*x[4]+0.1, sy=0.1, sz=0.1;
3231 //Double_t sy3=25000*x[4]*x[4]*60+0.5, sy=0.1, sz=0.1;
3233 Double_t f40=(F1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
3234 Double_t f42=(F1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
3235 Double_t f43=(F1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
3236 Double_t f20=(F2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
3237 Double_t f22=(F2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
3238 Double_t f23=(F2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
3240 Double_t f30=(F3(x1,y1+sy,x2,y2,z1,z2)-x[3])/sy;
3241 Double_t f31=(F3(x1,y1,x2,y2,z1+sz,z2)-x[3])/sz;
3242 Double_t f32=(F3(x1,y1,x2,y2+sy,z1,z2)-x[3])/sy;
3243 Double_t f34=(F3(x1,y1,x2,y2,z1,z2+sz)-x[3])/sz;
3247 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3248 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3249 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3250 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3251 c[13]=f30*sy1*f40+f32*sy2*f42;
3252 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3254 // if (!BuildSeed(kr1[is],kcl,0,x1,x2,x3,x,c)) continue;
3256 UInt_t index=kr1.GetIndex(is);
3257 seed->~AliTPCseed(); // this does not set the pointer to 0...
3258 AliTPCseed *track=new(seed) AliTPCseed(x1, ns*alpha+shift, x, c, index);
3260 track->SetIsSeeding(kTRUE);
3261 track->SetSeed1(i1);
3262 track->SetSeed2(i2);
3263 track->SetSeedType(3);
3267 FollowProlongation(*track, (i1+i2)/2,1);
3268 Int_t foundable,found,shared;
3269 track->GetClusterStatistic((i1+i2)/2,i1, found, foundable, shared, kTRUE);
3270 if ((found<0.55*foundable) || shared>0.5*found || (track->GetSigmaY2()+track->GetSigmaZ2())>0.5){
3272 seed->~AliTPCseed();
3278 FollowProlongation(*track, i2,1);
3282 track->SetBConstrain(1);
3283 // track->fLastPoint = i1+fInnerSec->GetNRows(); // first cluster in track position
3284 track->SetLastPoint(i1); // first cluster in track position
3285 track->SetFirstPoint(track->GetLastPoint());
3287 if (track->GetNumberOfClusters()<(i1-i2)*0.5 ||
3288 track->GetNumberOfClusters() < track->GetNFoundable()*0.6 ||
3289 track->GetNShared()>0.4*track->GetNumberOfClusters() ) {
3291 seed->~AliTPCseed();
3295 // Z VERTEX CONDITION
3296 Double_t zv, bz=GetBz();
3297 if ( !track->GetZAt(0.,bz,zv) ) continue;
3298 if (TMath::Abs(zv-z3)>cuts[2]) {
3299 FollowProlongation(*track, TMath::Max(i2-20,0));
3300 if ( !track->GetZAt(0.,bz,zv) ) continue;
3301 if (TMath::Abs(zv-z3)>cuts[2]){
3302 FollowProlongation(*track, TMath::Max(i2-40,0));
3303 if ( !track->GetZAt(0.,bz,zv) ) continue;
3304 if (TMath::Abs(zv-z3)>cuts[2] &&(track->GetNumberOfClusters() > track->GetNFoundable()*0.7)){
3305 // make seed without constrain
3306 AliTPCseed * track2 = MakeSeed(track,0.2,0.5,1.);
3307 FollowProlongation(*track2, i2,1);
3308 track2->SetBConstrain(kFALSE);
3309 track2->SetSeedType(1);
3310 arr->AddLast(track2);
3312 seed->~AliTPCseed();
3317 seed->~AliTPCseed();
3324 track->SetSeedType(0);
3325 arr->AddLast(track);
3326 seed = new AliTPCseed;
3328 // don't consider other combinations
3329 if (track->GetNumberOfClusters() > track->GetNFoundable()*0.8)
3335 Info("MakeSeeds3","\nSeeding statistic:\t%d\t%d\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin,nout1,nout2);
3341 void AliTPCtrackerMI::MakeSeeds5(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, Float_t cuts[4],
3346 //-----------------------------------------------------------------
3347 // This function creates track seeds.
3348 //-----------------------------------------------------------------
3349 // cuts[0] - fP4 cut
3350 // cuts[1] - tan(phi) cut
3351 // cuts[2] - zvertex cut
3352 // cuts[3] - fP3 cut
3362 Double_t x[5], c[15];
3364 // make temporary seed
3365 AliTPCseed * seed = new AliTPCseed;
3366 Double_t alpha=fOuterSec->GetAlpha(), shift=fOuterSec->GetAlphaShift();
3367 // Double_t cs=cos(alpha), sn=sin(alpha);
3372 Double_t x1 = GetXrow(i1-1);
3373 const AliTPCtrackerRow& kr1=GetRow(sec,i1-1);
3374 Double_t y1max = GetMaxY(i1-1)-kr1.GetDeadZone()-1.5;
3376 Double_t x1p = GetXrow(i1);
3377 const AliTPCtrackerRow& kr1p=GetRow(sec,i1);
3379 Double_t x1m = GetXrow(i1-2);
3380 const AliTPCtrackerRow& kr1m=GetRow(sec,i1-2);
3383 //last 3 padrow for seeding
3384 AliTPCtrackerRow& kr3 = GetRow((sec+fkNOS)%fkNOS,i1-7);
3385 Double_t x3 = GetXrow(i1-7);
3386 // Double_t y3max= GetMaxY(i1-7)-kr3.fDeadZone-1.5;
3388 AliTPCtrackerRow& kr3p = GetRow((sec+fkNOS)%fkNOS,i1-6);
3389 Double_t x3p = GetXrow(i1-6);
3391 AliTPCtrackerRow& kr3m = GetRow((sec+fkNOS)%fkNOS,i1-8);
3392 Double_t x3m = GetXrow(i1-8);
3397 Int_t im = i1-4; //middle pad row index
3398 Double_t xm = GetXrow(im); // radius of middle pad-row
3399 const AliTPCtrackerRow& krm=GetRow(sec,im); //middle pad -row
3400 // Double_t ymmax = GetMaxY(im)-kr1.fDeadZone-1.5;
3403 Double_t deltax = x1-x3;
3404 Double_t dymax = deltax*cuts[1];
3405 Double_t dzmax = deltax*cuts[3];
3407 // loop over clusters
3408 for (Int_t is=0; is < kr1; is++) {
3410 if (kr1[is]->IsUsed(10)) continue;
3411 Double_t y1=kr1[is]->GetY(), z1=kr1[is]->GetZ();
3413 if (deltay>0 && TMath::Abs(y1max-TMath::Abs(y1))> deltay ) continue; // seed only at the edge
3415 Int_t index1 = TMath::Max(kr3.Find(z1-dzmax)-1,0);
3416 Int_t index2 = TMath::Min(kr3.Find(z1+dzmax)+1,kr3);
3422 for (Int_t js=index1; js < index2; js++) {
3423 const AliTPCclusterMI *kcl = kr3[js];
3424 if (kcl->IsUsed(10)) continue;
3426 // apply angular cuts
3427 if (TMath::Abs(y1-y3)>dymax) continue;
3430 if (TMath::Abs(z1-z3)>dzmax) continue;
3432 Double_t angley = (y1-y3)/(x1-x3);
3433 Double_t anglez = (z1-z3)/(x1-x3);
3435 Double_t erry = TMath::Abs(angley)*(x1-x1m)*0.5+0.5;
3436 Double_t errz = TMath::Abs(anglez)*(x1-x1m)*0.5+0.5;
3438 Double_t yyym = angley*(xm-x1)+y1;
3439 Double_t zzzm = anglez*(xm-x1)+z1;
3441 const AliTPCclusterMI *kcm = krm.FindNearest2(yyym,zzzm,erry,errz,index);
3443 if (kcm->IsUsed(10)) continue;
3445 erry = TMath::Abs(angley)*(x1-x1m)*0.4+0.5;
3446 errz = TMath::Abs(anglez)*(x1-x1m)*0.4+0.5;
3453 // look around first
3454 const AliTPCclusterMI *kc1m = kr1m.FindNearest2(angley*(x1m-x1)+y1,
3460 if (kc1m->IsUsed(10)) used++;
3462 const AliTPCclusterMI *kc1p = kr1p.FindNearest2(angley*(x1p-x1)+y1,
3468 if (kc1p->IsUsed(10)) used++;
3470 if (used>1) continue;
3471 if (found<1) continue;
3475 const AliTPCclusterMI *kc3m = kr3m.FindNearest2(angley*(x3m-x3)+y3,
3481 if (kc3m->IsUsed(10)) used++;
3485 const AliTPCclusterMI *kc3p = kr3p.FindNearest2(angley*(x3p-x3)+y3,
3491 if (kc3p->IsUsed(10)) used++;
3495 if (used>1) continue;
3496 if (found<3) continue;
3506 x[4]=F1(x1,y1,x2,y2,x3,y3);
3507 //if (TMath::Abs(x[4]) >= cuts[0]) continue;
3510 x[2]=F2(x1,y1,x2,y2,x3,y3);
3513 x[3]=F3n(x1,y1,x2,y2,z1,z2,x[4]);
3514 //if (TMath::Abs(x[3]) > cuts[3]) continue;
3518 Double_t sy1=0.1, sz1=0.1;
3519 Double_t sy2=0.1, sz2=0.1;
3520 Double_t sy3=0.1, sy=0.1, sz=0.1;
3522 Double_t f40=(F1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
3523 Double_t f42=(F1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
3524 Double_t f43=(F1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
3525 Double_t f20=(F2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
3526 Double_t f22=(F2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
3527 Double_t f23=(F2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
3529 Double_t f30=(F3(x1,y1+sy,x2,y2,z1,z2)-x[3])/sy;
3530 Double_t f31=(F3(x1,y1,x2,y2,z1+sz,z2)-x[3])/sz;
3531 Double_t f32=(F3(x1,y1,x2,y2+sy,z1,z2)-x[3])/sy;
3532 Double_t f34=(F3(x1,y1,x2,y2,z1,z2+sz)-x[3])/sz;
3536 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3537 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3538 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3539 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3540 c[13]=f30*sy1*f40+f32*sy2*f42;
3541 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3543 // if (!BuildSeed(kr1[is],kcl,0,x1,x2,x3,x,c)) continue;
3545 index=kr1.GetIndex(is);
3546 seed->~AliTPCseed();
3547 AliTPCseed *track=new(seed) AliTPCseed(x1, sec*alpha+shift, x, c, index);
3549 track->SetIsSeeding(kTRUE);
3552 FollowProlongation(*track, i1-7,1);
3553 if (track->GetNumberOfClusters() < track->GetNFoundable()*0.75 ||
3554 track->GetNShared()>0.6*track->GetNumberOfClusters() || ( track->GetSigmaY2()+ track->GetSigmaZ2())>0.6){
3556 seed->~AliTPCseed();
3562 FollowProlongation(*track, i2,1);
3563 track->SetBConstrain(0);
3564 track->SetLastPoint(i1+fInnerSec->GetNRows()); // first cluster in track position
3565 track->SetFirstPoint(track->GetLastPoint());
3567 if (track->GetNumberOfClusters()<(i1-i2)*0.5 ||
3568 track->GetNumberOfClusters()<track->GetNFoundable()*0.7 ||
3569 track->GetNShared()>2. || track->GetChi2()/track->GetNumberOfClusters()>6 || ( track->GetSigmaY2()+ track->GetSigmaZ2())>0.5 ) {
3571 seed->~AliTPCseed();
3576 FollowProlongation(*track, TMath::Max(i2-10,0),1);
3577 AliTPCseed * track2 = MakeSeed(track,0.2,0.5,0.9);
3578 FollowProlongation(*track2, i2,1);
3579 track2->SetBConstrain(kFALSE);
3580 track2->SetSeedType(4);
3581 arr->AddLast(track2);
3583 seed->~AliTPCseed();
3587 //arr->AddLast(track);
3588 //seed = new AliTPCseed;
3594 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);
3600 //_____________________________________________________________________________
3601 void AliTPCtrackerMI::MakeSeeds2(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, Float_t */*cuts[4]*/,
3602 Float_t deltay, Bool_t /*bconstrain*/) {
3603 //-----------------------------------------------------------------
3604 // This function creates track seeds - without vertex constraint
3605 //-----------------------------------------------------------------
3606 // cuts[0] - fP4 cut - not applied
3607 // cuts[1] - tan(phi) cut
3608 // cuts[2] - zvertex cut - not applied
3609 // cuts[3] - fP3 cut
3619 Double_t alpha=fOuterSec->GetAlpha(), shift=fOuterSec->GetAlphaShift();
3620 // Double_t cs=cos(alpha), sn=sin(alpha);
3621 Int_t row0 = (i1+i2)/2;
3622 Int_t drow = (i1-i2)/2;
3623 const AliTPCtrackerRow& kr0=fSectors[sec][row0];
3624 AliTPCtrackerRow * kr=0;
3626 AliTPCpolyTrack polytrack;
3627 Int_t nclusters=fSectors[sec][row0];
3628 AliTPCseed * seed = new AliTPCseed;
3633 for (Int_t is=0; is < nclusters; is++) { //LOOP over clusters
3635 Int_t nfoundable =0;
3636 for (Int_t iter =1; iter<2; iter++){ //iterations
3637 const AliTPCtrackerRow& krm=fSectors[sec][row0-iter];
3638 const AliTPCtrackerRow& krp=fSectors[sec][row0+iter];
3639 const AliTPCclusterMI * cl= kr0[is];
3641 if (cl->IsUsed(10)) {
3647 Double_t x = kr0.GetX();
3648 // Initialization of the polytrack
3653 Double_t y0= cl->GetY();
3654 Double_t z0= cl->GetZ();
3658 Double_t ymax = fSectors->GetMaxY(row0)-kr0.GetDeadZone()-1.5;
3659 if (deltay>0 && TMath::Abs(ymax-TMath::Abs(y0))> deltay ) continue; // seed only at the edge
3661 erry = (0.5)*cl->GetSigmaY2()/TMath::Sqrt(cl->GetQ())*6;
3662 errz = (0.5)*cl->GetSigmaZ2()/TMath::Sqrt(cl->GetQ())*6;
3663 polytrack.AddPoint(x,y0,z0,erry, errz);
3666 if (cl->IsUsed(10)) sumused++;
3669 Float_t roady = (5*TMath::Sqrt(cl->GetSigmaY2()+0.2)+1.)*iter;
3670 Float_t roadz = (5*TMath::Sqrt(cl->GetSigmaZ2()+0.2)+1.)*iter;
3673 AliTPCclusterMI * cl1 = krm.FindNearest(y0,z0,roady,roadz);
3674 if (cl1 && TMath::Abs(ymax-TMath::Abs(y0))) {
3675 erry = (0.5)*cl1->GetSigmaY2()/TMath::Sqrt(cl1->GetQ())*3;
3676 errz = (0.5)*cl1->GetSigmaZ2()/TMath::Sqrt(cl1->GetQ())*3;
3677 if (cl1->IsUsed(10)) sumused++;
3678 polytrack.AddPoint(x,cl1->GetY(),cl1->GetZ(),erry,errz);
3682 AliTPCclusterMI * cl2 = krp.FindNearest(y0,z0,roady,roadz);
3684 erry = (0.5)*cl2->GetSigmaY2()/TMath::Sqrt(cl2->GetQ())*3;
3685 errz = (0.5)*cl2->GetSigmaZ2()/TMath::Sqrt(cl2->GetQ())*3;
3686 if (cl2->IsUsed(10)) sumused++;
3687 polytrack.AddPoint(x,cl2->GetY(),cl2->GetZ(),erry,errz);
3690 if (sumused>0) continue;
3692 polytrack.UpdateParameters();
3698 nfoundable = polytrack.GetN();
3699 nfound = nfoundable;
3701 for (Int_t ddrow = iter+1; ddrow<drow;ddrow++){
3702 Float_t maxdist = 0.8*(1.+3./(ddrow));
3703 for (Int_t delta = -1;delta<=1;delta+=2){
3704 Int_t row = row0+ddrow*delta;
3705 kr = &(fSectors[sec][row]);
3706 Double_t xn = kr->GetX();
3707 Double_t ymax1 = fSectors->GetMaxY(row)-kr->GetDeadZone()-1.5;
3708 polytrack.GetFitPoint(xn,yn,zn);
3709 if (TMath::Abs(yn)>ymax1) continue;
3711 AliTPCclusterMI * cln = kr->FindNearest(yn,zn,roady,roadz);
3713 Float_t dist = TMath::Sqrt( (yn-cln->GetY())*(yn-cln->GetY())+(zn-cln->GetZ())*(zn-cln->GetZ()));
3716 erry = (dist+0.3)*cln->GetSigmaY2()/TMath::Sqrt(cln->GetQ())*(1.+1./(ddrow));
3717 errz = (dist+0.3)*cln->GetSigmaZ2()/TMath::Sqrt(cln->GetQ())*(1.+1./(ddrow));
3718 if (cln->IsUsed(10)) {
3719 // printf("used\n");
3727 polytrack.AddPoint(xn,cln->GetY(),cln->GetZ(),erry, errz);
3732 if ( (sumused>3) || (sumused>0.5*nfound) || (nfound<0.6*nfoundable)) break;
3733 polytrack.UpdateParameters();
3736 if ( (sumused>3) || (sumused>0.5*nfound)) {
3737 //printf("sumused %d\n",sumused);
3742 polytrack.GetFitDerivation(kr0.GetX(),dy,dz);
3743 AliTPCpolyTrack track2;
3745 polytrack.Refit(track2,0.5+TMath::Abs(dy)*0.3,0.4+TMath::Abs(dz)*0.3);
3746 if (track2.GetN()<0.5*nfoundable) continue;
3749 if ((nfound>0.6*nfoundable) &&( nfoundable>0.4*(i1-i2))) {
3751 // test seed with and without constrain
3752 for (Int_t constrain=0; constrain<=0;constrain++){
3753 // add polytrack candidate
3755 Double_t x[5], c[15];
3756 Double_t x1,x2,x3,y1,y2,y3,z1,z2,z3;
3757 track2.GetBoundaries(x3,x1);
3759 track2.GetFitPoint(x1,y1,z1);
3760 track2.GetFitPoint(x2,y2,z2);
3761 track2.GetFitPoint(x3,y3,z3);
3763 //is track pointing to the vertex ?
3766 polytrack.GetFitPoint(x0,y0,z0);
3779 x[4]=F1(x1,y1,x2,y2,x3,y3);
3781 // if (TMath::Abs(x[4]) >= cuts[0]) continue; //
3782 x[2]=F2(x1,y1,x2,y2,x3,y3);
3784 //if (TMath::Abs(x[4]*x1-x[2]) >= cuts[1]) continue;
3785 //x[3]=F3(x1,y1,x2,y2,z1,z2);
3786 x[3]=F3n(x1,y1,x3,y3,z1,z3,x[4]);
3787 //if (TMath::Abs(x[3]) > cuts[3]) continue;
3790 Double_t sy =0.1, sz =0.1;
3791 Double_t sy1=0.02, sz1=0.02;
3792 Double_t sy2=0.02, sz2=0.02;
3796 sy3=25000*x[4]*x[4]+0.1, sy=0.1, sz=0.1;
3799 Double_t f40=(F1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
3800 Double_t f42=(F1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
3801 Double_t f43=(F1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
3802 Double_t f20=(F2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
3803 Double_t f22=(F2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
3804 Double_t f23=(F2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
3806 Double_t f30=(F3(x1,y1+sy,x3,y3,z1,z3)-x[3])/sy;
3807 Double_t f31=(F3(x1,y1,x3,y3,z1+sz,z3)-x[3])/sz;
3808 Double_t f32=(F3(x1,y1,x3,y3+sy,z1,z3)-x[3])/sy;
3809 Double_t f34=(F3(x1,y1,x3,y3,z1,z3+sz)-x[3])/sz;
3814 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3815 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3816 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3817 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3818 c[13]=f30*sy1*f40+f32*sy2*f42;
3819 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3821 //Int_t row1 = fSectors->GetRowNumber(x1);
3822 Int_t row1 = GetRowNumber(x1);
3826 seed->~AliTPCseed();
3827 AliTPCseed *track=new(seed) AliTPCseed(x1,sec*alpha+shift,x,c,index);
3828 track->SetIsSeeding(kTRUE);
3829 Int_t rc=FollowProlongation(*track, i2);
3830 if (constrain) track->SetBConstrain(1);
3832 track->SetBConstrain(0);
3833 track->SetLastPoint(row1+fInnerSec->GetNRows()); // first cluster in track position
3834 track->SetFirstPoint(track->GetLastPoint());
3836 if (rc==0 || track->GetNumberOfClusters()<(i1-i2)*0.5 ||
3837 track->GetNumberOfClusters() < track->GetNFoundable()*0.6 ||
3838 track->GetNShared()>0.4*track->GetNumberOfClusters()) {
3841 seed->~AliTPCseed();
3844 arr->AddLast(track);
3845 seed = new AliTPCseed;
3849 } // if accepted seed
3852 Info("MakeSeeds2","\nSeeding statiistic:\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin3);
3858 AliTPCseed *AliTPCtrackerMI::MakeSeed(AliTPCseed *const track, Float_t r0, Float_t r1, Float_t r2)
3862 //reseed using track points
3863 Int_t p0 = int(r0*track->GetNumberOfClusters()); // point 0
3864 Int_t p1 = int(r1*track->GetNumberOfClusters());
3865 Int_t p2 = int(r2*track->GetNumberOfClusters()); // last point
3867 Double_t x0[3],x1[3],x2[3];
3868 for (Int_t i=0;i<3;i++){
3874 // find track position at given ratio of the length
3875 Int_t sec0=0, sec1=0, sec2=0;
3878 for (Int_t i=0;i<160;i++){
3879 if (track->GetClusterPointer(i)){
3881 AliTPCTrackerPoint *trpoint =track->GetTrackPoint(i);
3882 if ( (index<p0) || x0[0]<0 ){
3883 if (trpoint->GetX()>1){
3884 clindex = track->GetClusterIndex2(i);
3886 x0[0] = trpoint->GetX();
3887 x0[1] = trpoint->GetY();
3888 x0[2] = trpoint->GetZ();
3889 sec0 = ((clindex&0xff000000)>>24)%18;
3894 if ( (index<p1) &&(trpoint->GetX()>1)){
3895 clindex = track->GetClusterIndex2(i);
3897 x1[0] = trpoint->GetX();
3898 x1[1] = trpoint->GetY();
3899 x1[2] = trpoint->GetZ();
3900 sec1 = ((clindex&0xff000000)>>24)%18;
3903 if ( (index<p2) &&(trpoint->GetX()>1)){
3904 clindex = track->GetClusterIndex2(i);
3906 x2[0] = trpoint->GetX();
3907 x2[1] = trpoint->GetY();
3908 x2[2] = trpoint->GetZ();
3909 sec2 = ((clindex&0xff000000)>>24)%18;
3916 Double_t alpha, cs,sn, xx2,yy2;
3918 alpha = (sec1-sec2)*fSectors->GetAlpha();
3919 cs = TMath::Cos(alpha);
3920 sn = TMath::Sin(alpha);
3921 xx2= x1[0]*cs-x1[1]*sn;
3922 yy2= x1[0]*sn+x1[1]*cs;
3926 alpha = (sec0-sec2)*fSectors->GetAlpha();
3927 cs = TMath::Cos(alpha);
3928 sn = TMath::Sin(alpha);
3929 xx2= x0[0]*cs-x0[1]*sn;
3930 yy2= x0[0]*sn+x0[1]*cs;
3936 Double_t x[5],c[15];
3940 x[4]=F1(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]);
3941 // if (x[4]>1) return 0;
3942 x[2]=F2(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]);
3943 x[3]=F3n(x2[0],x2[1],x0[0],x0[1],x2[2],x0[2],x[4]);
3944 //if (TMath::Abs(x[3]) > 2.2) return 0;
3945 //if (TMath::Abs(x[2]) > 1.99) return 0;
3947 Double_t sy =0.1, sz =0.1;
3949 Double_t sy1=0.02+track->GetSigmaY2(), sz1=0.02+track->GetSigmaZ2();
3950 Double_t sy2=0.01+track->GetSigmaY2(), sz2=0.01+track->GetSigmaZ2();
3951 Double_t sy3=0.01+track->GetSigmaY2();
3953 Double_t f40=(F1(x2[0],x2[1]+sy,x1[0],x1[1],x0[0],x0[1])-x[4])/sy;
3954 Double_t f42=(F1(x2[0],x2[1],x1[0],x1[1]+sy,x0[0],x0[1])-x[4])/sy;
3955 Double_t f43=(F1(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]+sy)-x[4])/sy;
3956 Double_t f20=(F2(x2[0],x2[1]+sy,x1[0],x1[1],x0[0],x0[1])-x[2])/sy;
3957 Double_t f22=(F2(x2[0],x2[1],x1[0],x1[1]+sy,x0[0],x0[1])-x[2])/sy;
3958 Double_t f23=(F2(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]+sy)-x[2])/sy;
3960 Double_t f30=(F3(x2[0],x2[1]+sy,x0[0],x0[1],x2[2],x0[2])-x[3])/sy;
3961 Double_t f31=(F3(x2[0],x2[1],x0[0],x0[1],x2[2]+sz,x0[2])-x[3])/sz;
3962 Double_t f32=(F3(x2[0],x2[1],x0[0],x0[1]+sy,x2[2],x0[2])-x[3])/sy;
3963 Double_t f34=(F3(x2[0],x2[1],x0[0],x0[1],x2[2],x0[2]+sz)-x[3])/sz;
3968 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3969 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3970 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3971 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3972 c[13]=f30*sy1*f40+f32*sy2*f42;
3973 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3975 // Int_t row1 = fSectors->GetRowNumber(x2[0]);
3976 AliTPCseed *seed=new AliTPCseed(x2[0], sec2*fSectors->GetAlpha()+fSectors->GetAlphaShift(), x, c, 0);
3977 // Double_t y0,z0,y1,z1, y2,z2;
3978 //seed->GetProlongation(x0[0],y0,z0);
3979 // seed->GetProlongation(x1[0],y1,z1);
3980 //seed->GetProlongation(x2[0],y2,z2);
3982 seed->SetLastPoint(pp2);
3983 seed->SetFirstPoint(pp2);
3990 AliTPCseed *AliTPCtrackerMI::ReSeed(const AliTPCseed *track, Float_t r0, Float_t r1, Float_t r2)
3994 //reseed using founded clusters
3996 // Find the number of clusters
3997 Int_t nclusters = 0;
3998 for (Int_t irow=0;irow<160;irow++){
3999 if (track->GetClusterIndex(irow)>0) nclusters++;
4003 ipos[0] = TMath::Max(int(r0*nclusters),0); // point 0 cluster
4004 ipos[1] = TMath::Min(int(r1*nclusters),nclusters-1); //
4005 ipos[2] = TMath::Min(int(r2*nclusters),nclusters-1); // last point
4008 Double_t xyz[3][3]={{0}};
4009 Int_t row[3]={0},sec[3]={0,0,0};
4011 // find track row position at given ratio of the length
4013 for (Int_t irow=0;irow<160;irow++){
4014 if (track->GetClusterIndex2(irow)<0) continue;
4016 for (Int_t ipoint=0;ipoint<3;ipoint++){
4017 if (index<=ipos[ipoint]) row[ipoint] = irow;
4021 //Get cluster and sector position
4022 for (Int_t ipoint=0;ipoint<3;ipoint++){
4023 Int_t clindex = track->GetClusterIndex2(row[ipoint]);
4024 AliTPCclusterMI * cl = GetClusterMI(clindex);
4027 // AliTPCclusterMI * cl = GetClusterMI(clindex);
4030 sec[ipoint] = ((clindex&0xff000000)>>24)%18;
4031 xyz[ipoint][0] = GetXrow(row[ipoint]);
4032 xyz[ipoint][1] = cl->GetY();
4033 xyz[ipoint][2] = cl->GetZ();
4037 // Calculate seed state vector and covariance matrix
4039 Double_t alpha, cs,sn, xx2,yy2;
4041 alpha = (sec[1]-sec[2])*fSectors->GetAlpha();
4042 cs = TMath::Cos(alpha);
4043 sn = TMath::Sin(alpha);
4044 xx2= xyz[1][0]*cs-xyz[1][1]*sn;
4045 yy2= xyz[1][0]*sn+xyz[1][1]*cs;
4049 alpha = (sec[0]-sec[2])*fSectors->GetAlpha();
4050 cs = TMath::Cos(alpha);
4051 sn = TMath::Sin(alpha);
4052 xx2= xyz[0][0]*cs-xyz[0][1]*sn;
4053 yy2= xyz[0][0]*sn+xyz[0][1]*cs;
4059 Double_t x[5],c[15];
4063 x[4]=F1(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
4064 x[2]=F2(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
4065 x[3]=F3n(xyz[2][0],xyz[2][1],xyz[0][0],xyz[0][1],xyz[2][2],xyz[0][2],x[4]);
4067 Double_t sy =0.1, sz =0.1;
4069 Double_t sy1=0.2, sz1=0.2;
4070 Double_t sy2=0.2, sz2=0.2;
4073 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;
4074 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;
4075 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;
4076 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;
4077 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;
4078 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;
4080 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;
4081 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;
4082 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;
4083 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;
4088 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
4089 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
4090 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
4091 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
4092 c[13]=f30*sy1*f40+f32*sy2*f42;
4093 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
4095 // Int_t row1 = fSectors->GetRowNumber(xyz[2][0]);
4096 AliTPCseed *seed=new AliTPCseed(xyz[2][0], sec[2]*fSectors->GetAlpha()+fSectors->GetAlphaShift(), x, c, 0);
4097 seed->SetLastPoint(row[2]);
4098 seed->SetFirstPoint(row[2]);
4103 AliTPCseed *AliTPCtrackerMI::ReSeed(AliTPCseed *track,Int_t r0, Bool_t forward)
4107 //reseed using founded clusters
4110 Int_t row[3]={0,0,0};
4111 Int_t sec[3]={0,0,0};
4113 // forward direction
4115 for (Int_t irow=r0;irow<160;irow++){
4116 if (track->GetClusterIndex(irow)>0){
4121 for (Int_t irow=160;irow>r0;irow--){
4122 if (track->GetClusterIndex(irow)>0){
4127 for (Int_t irow=row[2]-15;irow>row[0];irow--){
4128 if (track->GetClusterIndex(irow)>0){
4136 for (Int_t irow=0;irow<r0;irow++){
4137 if (track->GetClusterIndex(irow)>0){
4142 for (Int_t irow=r0;irow>0;irow--){
4143 if (track->GetClusterIndex(irow)>0){
4148 for (Int_t irow=row[2]-15;irow>row[0];irow--){
4149 if (track->GetClusterIndex(irow)>0){
4156 if ((row[2]-row[0])<20) return 0;
4157 if (row[1]==0) return 0;
4160 //Get cluster and sector position
4161 for (Int_t ipoint=0;ipoint<3;ipoint++){
4162 Int_t clindex = track->GetClusterIndex2(row[ipoint]);
4163 AliTPCclusterMI * cl = GetClusterMI(clindex);
4166 // AliTPCclusterMI * cl = GetClusterMI(clindex);
4169 sec[ipoint] = ((clindex&0xff000000)>>24)%18;
4170 xyz[ipoint][0] = GetXrow(row[ipoint]);
4171 AliTPCTrackerPoint * point = track->GetTrackPoint(row[ipoint]);
4172 if (point&&ipoint<2){
4174 xyz[ipoint][1] = point->GetY();
4175 xyz[ipoint][2] = point->GetZ();
4178 xyz[ipoint][1] = cl->GetY();
4179 xyz[ipoint][2] = cl->GetZ();
4186 // Calculate seed state vector and covariance matrix
4188 Double_t alpha, cs,sn, xx2,yy2;
4190 alpha = (sec[1]-sec[2])*fSectors->GetAlpha();
4191 cs = TMath::Cos(alpha);
4192 sn = TMath::Sin(alpha);
4193 xx2= xyz[1][0]*cs-xyz[1][1]*sn;
4194 yy2= xyz[1][0]*sn+xyz[1][1]*cs;
4198 alpha = (sec[0]-sec[2])*fSectors->GetAlpha();
4199 cs = TMath::Cos(alpha);
4200 sn = TMath::Sin(alpha);
4201 xx2= xyz[0][0]*cs-xyz[0][1]*sn;
4202 yy2= xyz[0][0]*sn+xyz[0][1]*cs;
4208 Double_t x[5],c[15];
4212 x[4]=F1(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
4213 x[2]=F2(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
4214 x[3]=F3n(xyz[2][0],xyz[2][1],xyz[0][0],xyz[0][1],xyz[2][2],xyz[0][2],x[4]);
4216 Double_t sy =0.1, sz =0.1;
4218 Double_t sy1=0.2, sz1=0.2;
4219 Double_t sy2=0.2, sz2=0.2;
4222 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;
4223 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;
4224 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;
4225 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;
4226 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;
4227 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;
4229 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;
4230 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;
4231 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;
4232 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;
4237 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
4238 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
4239 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
4240 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
4241 c[13]=f30*sy1*f40+f32*sy2*f42;
4242 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
4244 // Int_t row1 = fSectors->GetRowNumber(xyz[2][0]);
4245 AliTPCseed *seed=new AliTPCseed(xyz[2][0], sec[2]*fSectors->GetAlpha()+fSectors->GetAlphaShift(), x, c, 0);
4246 seed->SetLastPoint(row[2]);
4247 seed->SetFirstPoint(row[2]);
4248 for (Int_t i=row[0];i<row[2];i++){
4249 seed->SetClusterIndex(i, track->GetClusterIndex(i));
4257 void AliTPCtrackerMI::FindMultiMC(const TObjArray * array, AliESDEvent */*esd*/, Int_t iter)
4260 // find multi tracks - THIS FUNCTION IS ONLY FOR DEBUG PURPOSES
4262 // Use AliTPCReconstructor::StreamLevel()>2 if you want to tune parameters - cuts
4264 // Two reasons to have multiple find tracks
4265 // 1. Curling tracks can be find more than once
4266 // 2. Splitted tracks
4267 // a.) Multiple seeding to increase tracking efficiency - (~ 100% reached)
4268 // b.) Edge effect on the sector boundaries
4271 // Algorithm done in 2 phases - because of CPU consumption
4272 // it is n^2 algorithm - for lead-lead 20000x20000 combination are investigated
4274 // Algorihm for curling tracks sign:
4275 // 1 phase -makes a very rough fast cuts to minimize combinatorics
4276 // a.) opposite sign
4277 // b.) one of the tracks - not pointing to the primary vertex -
4278 // c.) delta tan(theta)
4280 // 2 phase - calculates DCA between tracks - time consument
4285 // General cuts - for splitted tracks and for curling tracks
4287 const Float_t kMaxdPhi = 0.2; // maximal distance in phi
4289 // Curling tracks cuts
4294 Int_t nentries = array->GetEntriesFast();
4295 AliHelix *helixes = new AliHelix[nentries];
4296 Float_t *xm = new Float_t[nentries];
4297 Float_t *dz0 = new Float_t[nentries];
4298 Float_t *dz1 = new Float_t[nentries];
4304 // Find track COG in x direction - point with best defined parameters
4306 for (Int_t i=0;i<nentries;i++){
4307 AliTPCseed* track = (AliTPCseed*)array->At(i);
4308 if (!track) continue;
4309 track->SetCircular(0);
4310 new (&helixes[i]) AliHelix(*track);
4314 track->GetDZ(GetX(),GetY(),GetZ(),GetBz(),dz);
4317 for (Int_t icl=0; icl<160; icl++){
4318 AliTPCclusterMI * cl = track->GetClusterPointer(icl);
4324 if (ncl>0) xm[i]/=Float_t(ncl);
4327 for (Int_t i0=0;i0<nentries;i0++){
4328 AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
4329 if (!track0) continue;
4330 Float_t xc0 = helixes[i0].GetHelix(6);
4331 Float_t yc0 = helixes[i0].GetHelix(7);
4332 Float_t r0 = helixes[i0].GetHelix(8);
4333 Float_t rc0 = TMath::Sqrt(xc0*xc0+yc0*yc0);
4334 Float_t fi0 = TMath::ATan2(yc0,xc0);
4336 for (Int_t i1=i0+1;i1<nentries;i1++){
4337 AliTPCseed * track1 = (AliTPCseed*)array->At(i1);
4338 if (!track1) continue;
4339 Int_t lab0=track0->GetLabel();
4340 Int_t lab1=track1->GetLabel();
4341 if (TMath::Abs(lab0)!=TMath::Abs(lab1)) continue;
4343 Float_t xc1 = helixes[i1].GetHelix(6);
4344 Float_t yc1 = helixes[i1].GetHelix(7);
4345 Float_t r1 = helixes[i1].GetHelix(8);
4346 Float_t rc1 = TMath::Sqrt(xc1*xc1+yc1*yc1);
4347 Float_t fi1 = TMath::ATan2(yc1,xc1);
4349 Float_t dfi = fi0-fi1;
4352 if (dfi>1.5*TMath::Pi()) dfi-=TMath::Pi(); // take care about edge effect
4353 if (dfi<-1.5*TMath::Pi()) dfi+=TMath::Pi(); //
4354 if (TMath::Abs(dfi)>kMaxdPhi&&helixes[i0].GetHelix(4)*helixes[i1].GetHelix(4)<0){
4356 // if short tracks with undefined sign
4357 fi1 = -TMath::ATan2(yc1,-xc1);
4360 Float_t dtheta = TMath::Abs(track0->GetTgl()-track1->GetTgl())<TMath::Abs(track0->GetTgl()+track1->GetTgl())? track0->GetTgl()-track1->GetTgl():track0->GetTgl()+track1->GetTgl();
4363 // debug stream to tune "fast cuts"
4365 Double_t dist[3]; // distance at X
4366 Double_t mdist[3]={0,0,0}; // mean distance X+-40cm
4367 track0->GetDistance(track1,0.5*(xm[i0]+xm[i1])-40.,dist,AliTracker::GetBz());
4368 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[i]);
4369 track0->GetDistance(track1,0.5*(xm[i0]+xm[i1])+40.,dist,AliTracker::GetBz());
4370 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[i]);
4371 track0->GetDistance(track1,0.5*(xm[i0]+xm[i1]),dist,AliTracker::GetBz());
4372 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[i]);
4373 for (Int_t i=0;i<3;i++) mdist[i]*=0.33333;
4377 for (Int_t icl=0; icl<160; icl++){
4378 AliTPCclusterMI * cl0 = track0->GetClusterPointer(icl);
4379 AliTPCclusterMI * cl1 = track1->GetClusterPointer(icl);
4382 if (cl0==cl1) sums++;
4386 if (AliTPCReconstructor::StreamLevel()>5) {
4387 TTreeSRedirector &cstream = *fDebugStreamer;
4392 "Tr0.="<<track0<< // seed0
4393 "Tr1.="<<track1<< // seed1
4394 "h0.="<<&helixes[i0]<<
4395 "h1.="<<&helixes[i1]<<
4397 "sum="<<sum<< //the sum of rows with cl in both
4398 "sums="<<sums<< //the sum of shared clusters
4399 "xm0="<<xm[i0]<< // the center of track
4400 "xm1="<<xm[i1]<< // the x center of track
4401 // General cut variables
4402 "dfi="<<dfi<< // distance in fi angle
4403 "dtheta="<<dtheta<< // distance int theta angle
4409 "dist0="<<dist[0]<< //distance x
4410 "dist1="<<dist[1]<< //distance y
4411 "dist2="<<dist[2]<< //distance z
4412 "mdist0="<<mdist[0]<< //distance x
4413 "mdist1="<<mdist[1]<< //distance y
4414 "mdist2="<<mdist[2]<< //distance z
4430 if (AliTPCReconstructor::StreamLevel()>1) {
4431 AliInfo("Time for curling tracks removal DEBUGGING MC");
4438 void AliTPCtrackerMI::FindSplitted(TObjArray * array, AliESDEvent */*esd*/, Int_t /*iter*/){
4440 // Find Splitted tracks and remove the one with worst quality
4441 // Corresponding debug streamer to tune selections - "Splitted2"
4443 // 0. Sort tracks according quility
4444 // 1. Propagate the tracks to the reference radius
4445 // 2. Double_t loop to select close tracks (only to speed up process)
4446 // 3. Calculate cluster overlap ratio - and remove the track if bigger than a threshold
4447 // 4. Delete temporary parameters
4449 const Double_t xref=GetXrow(63); // reference radius -IROC/OROC boundary
4451 const Double_t kCutP1=10; // delta Z cut 10 cm
4452 const Double_t kCutP2=0.15; // delta snp(fi) cut 0.15
4453 const Double_t kCutP3=0.15; // delta tgl(theta) cut 0.15
4454 const Double_t kCutAlpha=0.15; // delta alpha cut
4455 Int_t firstpoint = 0;
4456 Int_t lastpoint = 160;
4458 Int_t nentries = array->GetEntriesFast();
4459 AliExternalTrackParam *params = new AliExternalTrackParam[nentries];
4465 //0. Sort tracks according quality
4466 //1. Propagate the ext. param to reference radius
4467 Int_t nseed = array->GetEntriesFast();
4468 if (nseed<=0) return;
4469 Float_t * quality = new Float_t[nseed];
4470 Int_t * indexes = new Int_t[nseed];
4471 for (Int_t i=0; i<nseed; i++) {
4472 AliTPCseed *pt=(AliTPCseed*)array->UncheckedAt(i);
4477 pt->UpdatePoints(); //select first last max dens points
4478 Float_t * points = pt->GetPoints();
4479 if (points[3]<0.8) quality[i] =-1;
4480 quality[i] = (points[2]-points[0])+pt->GetNumberOfClusters();
4481 //prefer high momenta tracks if overlaps
4482 quality[i] *= TMath::Sqrt(TMath::Abs(pt->Pt())+0.5);
4484 AliTracker::PropagateTrackToBxByBz(&(params[i]),xref,pt->GetMass(),5.,kTRUE);
4485 AliTracker::PropagateTrackToBxByBz(&(params[i]),xref,pt->GetMass(),1.,kTRUE);
4487 TMath::Sort(nseed,quality,indexes);
4489 // 3. Loop over pair of tracks
4491 for (Int_t i0=0; i0<nseed; i0++) {
4492 Int_t index0=indexes[i0];
4493 if (!(array->UncheckedAt(index0))) continue;
4494 AliTPCseed *s1 = (AliTPCseed*)array->UncheckedAt(index0);
4495 if (!s1->IsActive()) continue;
4496 AliExternalTrackParam &par0=params[index0];
4497 for (Int_t i1=i0+1; i1<nseed; i1++) {
4498 Int_t index1=indexes[i1];
4499 if (!(array->UncheckedAt(index1))) continue;
4500 AliTPCseed *s2 = (AliTPCseed*)array->UncheckedAt(index1);
4501 if (!s2->IsActive()) continue;
4502 if (s2->GetKinkIndexes()[0]!=0)
4503 if (s2->GetKinkIndexes()[0] == -s1->GetKinkIndexes()[0]) continue;
4504 AliExternalTrackParam &par1=params[index1];
4505 if (TMath::Abs(par0.GetParameter()[3]-par1.GetParameter()[3])>kCutP3) continue;
4506 if (TMath::Abs(par0.GetParameter()[1]-par1.GetParameter()[1])>kCutP1) continue;
4507 if (TMath::Abs(par0.GetParameter()[2]-par1.GetParameter()[2])>kCutP2) continue;
4508 Double_t dAlpha= TMath::Abs(par0.GetAlpha()-par1.GetAlpha());
4509 if (dAlpha>TMath::Pi()) dAlpha-=TMath::Pi();
4510 if (TMath::Abs(dAlpha)>kCutAlpha) continue;
4515 Int_t firstShared=lastpoint, lastShared=firstpoint;
4516 Int_t firstRow=lastpoint, lastRow=firstpoint;
4518 for (Int_t i=firstpoint;i<lastpoint;i++){
4519 if (s1->GetClusterIndex2(i)>0) nall0++;
4520 if (s2->GetClusterIndex2(i)>0) nall1++;
4521 if (s1->GetClusterIndex2(i)>0 && s2->GetClusterIndex2(i)>0) {
4522 if (i<firstRow) firstRow=i;
4523 if (i>lastRow) lastRow=i;
4525 if ( (s1->GetClusterIndex2(i))==(s2->GetClusterIndex2(i)) && s1->GetClusterIndex2(i)>0) {
4526 if (i<firstShared) firstShared=i;
4527 if (i>lastShared) lastShared=i;
4531 Double_t ratio0 = Float_t(sumShared)/Float_t(TMath::Min(nall0+1,nall1+1));
4532 Double_t ratio1 = Float_t(sumShared)/Float_t(TMath::Max(nall0+1,nall1+1));
4534 if( AliTPCReconstructor::StreamLevel()>1){
4535 TTreeSRedirector &cstream = *fDebugStreamer;
4536 Int_t n0=s1->GetNumberOfClusters();
4537 Int_t n1=s2->GetNumberOfClusters();
4538 Int_t n0F=s1->GetNFoundable();
4539 Int_t n1F=s2->GetNFoundable();
4540 Int_t lab0=s1->GetLabel();
4541 Int_t lab1=s2->GetLabel();
4543 cstream<<"Splitted2"<<
4544 "iter="<<fIteration<<
4545 "lab0="<<lab0<< // MC label if exist
4546 "lab1="<<lab1<< // MC label if exist
4549 "ratio0="<<ratio0<< // shared ratio
4550 "ratio1="<<ratio1<< // shared ratio
4551 "p0.="<<&par0<< // track parameters
4553 "s0.="<<s1<< // full seed
4555 "n0="<<n0<< // number of clusters track 0
4556 "n1="<<n1<< // number of clusters track 1
4557 "nall0="<<nall0<< // number of clusters track 0
4558 "nall1="<<nall1<< // number of clusters track 1
4559 "n0F="<<n0F<< // number of findable
4560 "n1F="<<n1F<< // number of findable
4561 "shared="<<sumShared<< // number of shared clusters
4562 "firstS="<<firstShared<< // first and the last shared row
4563 "lastS="<<lastShared<<
4564 "firstRow="<<firstRow<< // first and the last row with cluster
4565 "lastRow="<<lastRow<< //
4569 // remove track with lower quality
4571 if (ratio0>AliTPCReconstructor::GetRecoParam()->GetCutSharedClusters(0) ||
4572 ratio1>AliTPCReconstructor::GetRecoParam()->GetCutSharedClusters(1)){
4576 delete array->RemoveAt(index1);
4581 // 4. Delete temporary array
4591 void AliTPCtrackerMI::FindCurling(const TObjArray * array, AliESDEvent */*esd*/, Int_t iter)
4594 // find Curling tracks
4595 // Use AliTPCReconstructor::StreamLevel()>1 if you want to tune parameters - cuts
4598 // Algorithm done in 2 phases - because of CPU consumption
4599 // it is n^2 algorithm - for lead-lead 20000x20000 combination are investigated
4600 // see detal in MC part what can be used to cut
4604 const Float_t kMaxC = 400; // maximal curvature to of the track
4605 const Float_t kMaxdTheta = 0.15; // maximal distance in theta
4606 const Float_t kMaxdPhi = 0.15; // maximal distance in phi
4607 const Float_t kPtRatio = 0.3; // ratio between pt
4608 const Float_t kMinDCAR = 2.; // distance to the primary vertex in r - see cpipe cut
4611 // Curling tracks cuts
4614 const Float_t kMaxDeltaRMax = 40; // distance in outer radius
4615 const Float_t kMaxDeltaRMin = 5.; // distance in lower radius - see cpipe cut
4616 const Float_t kMinAngle = 2.9; // angle between tracks
4617 const Float_t kMaxDist = 5; // biggest distance
4619 // The cuts can be tuned using the "MC information stored in Multi tree ==> see FindMultiMC
4622 TCut csign("csign","Tr0.fP[4]*Tr1.fP[4]<0"); //opposite sign
4623 TCut cmax("cmax","abs(Tr0.GetC())>1/400");
4624 TCut cda("cda","sqrt(dtheta^2+dfi^2)<0.15");
4625 TCut ccratio("ccratio","abs((Tr0.fP[4]+Tr1.fP[4])/(abs(Tr0.fP[4])+abs(Tr1.fP[4])))<0.3");
4626 TCut cpipe("cpipe", "min(abs(r0-rc0),abs(r1-rc1))>5");
4628 TCut cdrmax("cdrmax","abs(abs(rc0+r0)-abs(rc1+r1))<40")
4629 TCut cdrmin("cdrmin","abs(abs(rc0+r0)-abs(rc1+r1))<10")
4631 Multi->Draw("dfi","iter==0"+csign+cmax+cda+ccratio); ~94% of curling tracks fulfill
4632 Multi->Draw("min(abs(r0-rc0),abs(r1-rc1))","iter==0&&abs(lab1)==abs(lab0)"+csign+cmax+cda+ccratio+cpipe+cdrmin+cdrmax); //80%
4634 Curling2->Draw("dfi","iter==0&&abs(lab0)==abs(lab1)"+csign+cmax+cdtheta+cdfi+ccratio)
4640 Int_t nentries = array->GetEntriesFast();
4641 AliHelix *helixes = new AliHelix[nentries];
4642 for (Int_t i=0;i<nentries;i++){
4643 AliTPCseed* track = (AliTPCseed*)array->At(i);
4644 if (!track) continue;
4645 track->SetCircular(0);
4646 new (&helixes[i]) AliHelix(*track);
4652 Double_t phase[2][2]={{0,0},{0,0}},radius[2]={0,0};
4658 for (Int_t i0=0;i0<nentries;i0++){
4659 AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
4660 if (!track0) continue;
4661 if (TMath::Abs(track0->GetC())<1/kMaxC) continue;
4662 Float_t xc0 = helixes[i0].GetHelix(6);
4663 Float_t yc0 = helixes[i0].GetHelix(7);
4664 Float_t r0 = helixes[i0].GetHelix(8);
4665 Float_t rc0 = TMath::Sqrt(xc0*xc0+yc0*yc0);
4666 Float_t fi0 = TMath::ATan2(yc0,xc0);
4668 for (Int_t i1=i0+1;i1<nentries;i1++){
4669 AliTPCseed * track1 = (AliTPCseed*)array->At(i1);
4670 if (!track1) continue;
4671 if (TMath::Abs(track1->GetC())<1/kMaxC) continue;
4672 Float_t xc1 = helixes[i1].GetHelix(6);
4673 Float_t yc1 = helixes[i1].GetHelix(7);
4674 Float_t r1 = helixes[i1].GetHelix(8);
4675 Float_t rc1 = TMath::Sqrt(xc1*xc1+yc1*yc1);
4676 Float_t fi1 = TMath::ATan2(yc1,xc1);
4678 Float_t dfi = fi0-fi1;
4681 if (dfi>1.5*TMath::Pi()) dfi-=TMath::Pi(); // take care about edge effect
4682 if (dfi<-1.5*TMath::Pi()) dfi+=TMath::Pi(); //
4683 Float_t dtheta = TMath::Abs(track0->GetTgl()-track1->GetTgl())<TMath::Abs(track0->GetTgl()+track1->GetTgl())? track0->GetTgl()-track1->GetTgl():track0->GetTgl()+track1->GetTgl();
4687 if (track0->GetBConstrain()&&track1->GetBConstrain()) continue; // not constrained
4688 if (track1->GetSigned1Pt()*track0->GetSigned1Pt()>0) continue; // not the same sign
4689 if ( TMath::Abs(track1->GetTgl()+track0->GetTgl())>kMaxdTheta) continue; //distance in the Theta
4690 if ( TMath::Abs(dfi)>kMaxdPhi) continue; //distance in phi
4691 if ( TMath::Sqrt(dfi*dfi+dtheta*dtheta)>kMaxdPhi) continue; //common angular offset
4693 Float_t pt0 = track0->GetSignedPt();
4694 Float_t pt1 = track1->GetSignedPt();
4695 if ((TMath::Abs(pt0+pt1)/(TMath::Abs(pt0)+TMath::Abs(pt1)))>kPtRatio) continue;
4696 if ((iter==1) && TMath::Abs(TMath::Abs(rc0+r0)-TMath::Abs(rc1+r1))>kMaxDeltaRMax) continue;
4697 if ((iter!=1) &&TMath::Abs(TMath::Abs(rc0-r0)-TMath::Abs(rc1-r1))>kMaxDeltaRMin) continue;
4698 if (TMath::Min(TMath::Abs(rc0-r0),TMath::Abs(rc1-r1))<kMinDCAR) continue;
4701 // Now find closest approach
4705 Int_t npoints = helixes[i0].GetRPHIintersections(helixes[i1], phase, radius,10);
4706 if (npoints==0) continue;
4707 helixes[i0].GetClosestPhases(helixes[i1], phase);
4711 Double_t hangles[3];
4712 helixes[i0].Evaluate(phase[0][0],xyz0);
4713 helixes[i1].Evaluate(phase[0][1],xyz1);
4715 helixes[i0].GetAngle(phase[0][0],helixes[i1],phase[0][1],hangles);
4716 Double_t deltah[2],deltabest;
4717 if (TMath::Abs(hangles[2])<kMinAngle) continue;
4721 helixes[i0].ParabolicDCA(helixes[i1],phase[0][0],phase[0][1],radius[0],deltah[0],2);
4723 helixes[i0].ParabolicDCA(helixes[i1],phase[1][0],phase[1][1],radius[1],deltah[1],2);
4724 if (deltah[1]<deltah[0]) ibest=1;
4726 deltabest = TMath::Sqrt(deltah[ibest]);
4727 helixes[i0].Evaluate(phase[ibest][0],xyz0);
4728 helixes[i1].Evaluate(phase[ibest][1],xyz1);
4729 helixes[i0].GetAngle(phase[ibest][0],helixes[i1],phase[ibest][1],hangles);
4730 Double_t radiusbest = TMath::Sqrt(radius[ibest]);
4732 if (deltabest>kMaxDist) continue;
4733 // if (mindcar+mindcaz<40 && (TMath::Abs(hangles[2])<kMinAngle ||deltabest>3)) continue;
4734 Bool_t sign =kFALSE;
4735 if (hangles[2]>kMinAngle) sign =kTRUE;
4738 // circular[i0] = kTRUE;
4739 // circular[i1] = kTRUE;
4740 if (track0->OneOverPt()<track1->OneOverPt()){
4741 track0->SetCircular(track0->GetCircular()+1);
4742 track1->SetCircular(track1->GetCircular()+2);
4745 track1->SetCircular(track1->GetCircular()+1);
4746 track0->SetCircular(track0->GetCircular()+2);
4749 if (AliTPCReconstructor::StreamLevel()>2){
4751 //debug stream to tune "fine" cuts
4752 Int_t lab0=track0->GetLabel();
4753 Int_t lab1=track1->GetLabel();
4754 TTreeSRedirector &cstream = *fDebugStreamer;
4755 cstream<<"Curling2"<<
4771 "npoints="<<npoints<<
4772 "hangles0="<<hangles[0]<<
4773 "hangles1="<<hangles[1]<<
4774 "hangles2="<<hangles[2]<<
4777 "radius="<<radiusbest<<
4778 "deltabest="<<deltabest<<
4779 "phase0="<<phase[ibest][0]<<
4780 "phase1="<<phase[ibest][1]<<
4788 if (AliTPCReconstructor::StreamLevel()>1) {
4789 AliInfo("Time for curling tracks removal");
4798 void AliTPCtrackerMI::FindKinks(TObjArray * array, AliESDEvent *esd)
4805 TObjArray *kinks= new TObjArray(10000);
4806 // TObjArray *v0s= new TObjArray(10000);
4807 Int_t nentries = array->GetEntriesFast();
4808 AliHelix *helixes = new AliHelix[nentries];
4809 Int_t *sign = new Int_t[nentries];
4810 Int_t *nclusters = new Int_t[nentries];
4811 Float_t *alpha = new Float_t[nentries];
4812 AliKink *kink = new AliKink();
4813 Int_t * usage = new Int_t[nentries];
4814 Float_t *zm = new Float_t[nentries];
4815 Float_t *z0 = new Float_t[nentries];
4816 Float_t *fim = new Float_t[nentries];
4817 Float_t *shared = new Float_t[nentries];
4818 Bool_t *circular = new Bool_t[nentries];
4819 Float_t *dca = new Float_t[nentries];
4820 //const AliESDVertex * primvertex = esd->GetVertex();
4822 // nentries = array->GetEntriesFast();
4827 for (Int_t i=0;i<nentries;i++){
4830 AliTPCseed* track = (AliTPCseed*)array->At(i);
4831 if (!track) continue;
4832 track->SetCircular(0);
4834 track->UpdatePoints();
4835 if (( track->GetPoints()[2]- track->GetPoints()[0])>5 && track->GetPoints()[3]>0.8){
4837 nclusters[i]=track->GetNumberOfClusters();
4838 alpha[i] = track->GetAlpha();
4839 new (&helixes[i]) AliHelix(*track);
4841 helixes[i].Evaluate(0,xyz);
4842 sign[i] = (track->GetC()>0) ? -1:1;
4845 if (track->GetProlongation(x,y,z)){
4847 fim[i] = alpha[i]+TMath::ATan2(y,x);
4850 zm[i] = track->GetZ();
4854 circular[i]= kFALSE;
4855 if (track->GetProlongation(0,y,z)) z0[i] = z;
4856 dca[i] = track->GetD(0,0);
4862 Int_t ncandidates =0;
4865 Double_t phase[2][2]={{0,0},{0,0}},radius[2]={0,0};
4868 // Find circling track
4870 for (Int_t i0=0;i0<nentries;i0++){
4871 AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
4872 if (!track0) continue;
4873 if (track0->GetNumberOfClusters()<40) continue;
4874 if (TMath::Abs(1./track0->GetC())>200) continue;
4875 for (Int_t i1=i0+1;i1<nentries;i1++){
4876 AliTPCseed * track1 = (AliTPCseed*)array->At(i1);
4877 if (!track1) continue;
4878 if (track1->GetNumberOfClusters()<40) continue;
4879 if ( TMath::Abs(track1->GetTgl()+track0->GetTgl())>0.1) continue;
4880 if (track0->GetBConstrain()&&track1->GetBConstrain()) continue;
4881 if (TMath::Abs(1./track1->GetC())>200) continue;
4882 if (track1->GetSigned1Pt()*track0->GetSigned1Pt()>0) continue;
4883 if (track1->GetTgl()*track0->GetTgl()>0) continue;
4884 if (TMath::Max(TMath::Abs(1./track0->GetC()),TMath::Abs(1./track1->GetC()))>190) continue;
4885 if (track0->GetBConstrain()&&track1->OneOverPt()<track0->OneOverPt()) continue; //returning - lower momenta
4886 if (track1->GetBConstrain()&&track0->OneOverPt()<track1->OneOverPt()) continue; //returning - lower momenta
4888 Float_t mindcar = TMath::Min(TMath::Abs(dca[i0]),TMath::Abs(dca[i1]));
4889 if (mindcar<5) continue;
4890 Float_t mindcaz = TMath::Min(TMath::Abs(z0[i0]-GetZ()),TMath::Abs(z0[i1]-GetZ()));
4891 if (mindcaz<5) continue;
4892 if (mindcar+mindcaz<20) continue;
4895 Float_t xc0 = helixes[i0].GetHelix(6);
4896 Float_t yc0 = helixes[i0].GetHelix(7);
4897 Float_t r0 = helixes[i0].GetHelix(8);
4898 Float_t xc1 = helixes[i1].GetHelix(6);
4899 Float_t yc1 = helixes[i1].GetHelix(7);
4900 Float_t r1 = helixes[i1].GetHelix(8);
4902 Float_t rmean = (r0+r1)*0.5;
4903 Float_t delta =TMath::Sqrt((xc1-xc0)*(xc1-xc0)+(yc1-yc0)*(yc1-yc0));
4904 //if (delta>30) continue;
4905 if (delta>rmean*0.25) continue;
4906 if (TMath::Abs(r0-r1)/rmean>0.3) continue;
4908 Int_t npoints = helixes[i0].GetRPHIintersections(helixes[i1], phase, radius,10);
4909 if (npoints==0) continue;
4910 helixes[i0].GetClosestPhases(helixes[i1], phase);
4914 Double_t hangles[3];
4915 helixes[i0].Evaluate(phase[0][0],xyz0);
4916 helixes[i1].Evaluate(phase[0][1],xyz1);
4918 helixes[i0].GetAngle(phase[0][0],helixes[i1],phase[0][1],hangles);
4919 Double_t deltah[2],deltabest;
4920 if (hangles[2]<2.8) continue;
4923 helixes[i0].ParabolicDCA(helixes[i1],phase[0][0],phase[0][1],radius[0],deltah[0],2);
4925 helixes[i0].ParabolicDCA(helixes[i1],phase[1][0],phase[1][1],radius[1],deltah[1],2);
4926 if (deltah[1]<deltah[0]) ibest=1;
4928 deltabest = TMath::Sqrt(deltah[ibest]);
4929 helixes[i0].Evaluate(phase[ibest][0],xyz0);
4930 helixes[i1].Evaluate(phase[ibest][1],xyz1);
4931 helixes[i0].GetAngle(phase[ibest][0],helixes[i1],phase[ibest][1],hangles);
4932 Double_t radiusbest = TMath::Sqrt(radius[ibest]);
4934 if (deltabest>6) continue;
4935 if (mindcar+mindcaz<40 && (hangles[2]<3.12||deltabest>3)) continue;
4936 Bool_t lsign =kFALSE;
4937 if (hangles[2]>3.06) lsign =kTRUE;
4940 circular[i0] = kTRUE;
4941 circular[i1] = kTRUE;
4942 if (track0->OneOverPt()<track1->OneOverPt()){
4943 track0->SetCircular(track0->GetCircular()+1);
4944 track1->SetCircular(track1->GetCircular()+2);
4947 track1->SetCircular(track1->GetCircular()+1);
4948 track0->SetCircular(track0->GetCircular()+2);
4951 if (lsign&&AliTPCReconstructor::StreamLevel()>1){
4953 Int_t lab0=track0->GetLabel();
4954 Int_t lab1=track1->GetLabel();
4955 TTreeSRedirector &cstream = *fDebugStreamer;
4956 cstream<<"Curling"<<
4963 "mindcar="<<mindcar<<
4964 "mindcaz="<<mindcaz<<
4967 "npoints="<<npoints<<
4968 "hangles0="<<hangles[0]<<
4969 "hangles2="<<hangles[2]<<
4974 "radius="<<radiusbest<<
4975 "deltabest="<<deltabest<<
4976 "phase0="<<phase[ibest][0]<<
4977 "phase1="<<phase[ibest][1]<<
4987 for (Int_t i =0;i<nentries;i++){
4988 if (sign[i]==0) continue;
4989 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
4996 Double_t cradius0 = 40*40;
4997 Double_t cradius1 = 270*270;
5000 Double_t cdist3=0.55;
5001 for (Int_t j =i+1;j<nentries;j++){
5003 if (sign[j]*sign[i]<1) continue;
5004 if ( (nclusters[i]+nclusters[j])>200) continue;
5005 if ( (nclusters[i]+nclusters[j])<80) continue;
5006 if ( TMath::Abs(zm[i]-zm[j])>60.) continue;
5007 if ( TMath::Abs(fim[i]-fim[j])>0.6 && TMath::Abs(fim[i]-fim[j])<5.7 ) continue;
5008 //AliTPCseed * track1 = (AliTPCseed*)array->At(j); Double_t phase[2][2],radius[2];
5009 Int_t npoints = helixes[i].GetRPHIintersections(helixes[j], phase, radius,20);
5010 if (npoints<1) continue;
5013 if (radius[0]<cradius0||radius[0]>cradius1) continue;
5016 if ( (radius[0]<cradius0||radius[0]>cradius1) && (radius[1]<cradius0||radius[1]>cradius1) ) continue;
5019 Double_t delta1=10000,delta2=10000;
5020 // cuts on the intersection radius
5021 helixes[i].LinearDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta1);
5022 if (radius[0]<20&&delta1<1) continue; //intersection at vertex
5023 if (radius[0]<10&&delta1<3) continue; //intersection at vertex
5025 helixes[i].LinearDCA(helixes[j],phase[1][0],phase[1][1],radius[1],delta2);
5026 if (radius[1]<20&&delta2<1) continue; //intersection at vertex
5027 if (radius[1]<10&&delta2<3) continue; //intersection at vertex
5030 Double_t distance1 = TMath::Min(delta1,delta2);
5031 if (distance1>cdist1) continue; // cut on DCA linear approximation
5033 npoints = helixes[i].GetRPHIintersections(helixes[j], phase, radius,20);
5034 helixes[i].ParabolicDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta1);
5035 if (radius[0]<20&&delta1<1) continue; //intersection at vertex
5036 if (radius[0]<10&&delta1<3) continue; //intersection at vertex
5039 helixes[i].ParabolicDCA(helixes[j],phase[1][0],phase[1][1],radius[1],delta2);
5040 if (radius[1]<20&&delta2<1) continue; //intersection at vertex
5041 if (radius[1]<10&&delta2<3) continue; //intersection at vertex
5043 distance1 = TMath::Min(delta1,delta2);
5046 rkink = TMath::Sqrt(radius[0]);
5049 rkink = TMath::Sqrt(radius[1]);
5051 if (distance1>cdist2) continue;
5054 AliTPCseed * track1 = (AliTPCseed*)array->At(j);
5057 Int_t row0 = GetRowNumber(rkink);
5058 if (row0<10) continue;
5059 if (row0>150) continue;
5062 Float_t dens00=-1,dens01=-1;
5063 Float_t dens10=-1,dens11=-1;
5065 Int_t found,foundable,ishared;
5066 track0->GetClusterStatistic(0,row0-5, found, foundable,ishared,kFALSE);
5067 if (foundable>5) dens00 = Float_t(found)/Float_t(foundable);
5068 track0->GetClusterStatistic(row0+5,155, found, foundable,ishared,kFALSE);
5069 if (foundable>5) dens01 = Float_t(found)/Float_t(foundable);
5071 track1->GetClusterStatistic(0,row0-5, found, foundable,ishared,kFALSE);
5072 if (foundable>10) dens10 = Float_t(found)/Float_t(foundable);
5073 track1->GetClusterStatistic(row0+5,155, found, foundable,ishared,kFALSE);
5074 if (foundable>10) dens11 = Float_t(found)/Float_t(foundable);
5076 if (dens00<dens10 && dens01<dens11) continue;
5077 if (dens00>dens10 && dens01>dens11) continue;
5078 if (TMath::Max(dens00,dens10)<0.1) continue;
5079 if (TMath::Max(dens01,dens11)<0.3) continue;
5081 if (TMath::Min(dens00,dens10)>0.6) continue;
5082 if (TMath::Min(dens01,dens11)>0.6) continue;
5085 AliTPCseed * ktrack0, *ktrack1;
5094 if (TMath::Abs(ktrack0->GetC())>5) continue; // cut on the curvature for mother particle
5095 AliExternalTrackParam paramm(*ktrack0);
5096 AliExternalTrackParam paramd(*ktrack1);
5097 if (row0>60&&ktrack1->GetReference().GetX()>90.)new (¶md) AliExternalTrackParam(ktrack1->GetReference());
5100 kink->SetMother(paramm);
5101 kink->SetDaughter(paramd);
5104 Float_t x[3] = { kink->GetPosition()[0],kink->GetPosition()[1],kink->GetPosition()[2]};
5106 fkParam->Transform0to1(x,index);
5107 fkParam->Transform1to2(x,index);
5108 row0 = GetRowNumber(x[0]);
5110 if (kink->GetR()<100) continue;
5111 if (kink->GetR()>240) continue;
5112 if (kink->GetPosition()[2]/kink->GetR()>AliTPCReconstructor::GetCtgRange()) continue; //out of fiducial volume
5113 if (kink->GetDistance()>cdist3) continue;
5114 Float_t dird = kink->GetDaughterP()[0]*kink->GetPosition()[0]+kink->GetDaughterP()[1]*kink->GetPosition()[1]; // rough direction estimate
5115 if (dird<0) continue;
5117 Float_t dirm = kink->GetMotherP()[0]*kink->GetPosition()[0]+kink->GetMotherP()[1]*kink->GetPosition()[1]; // rough direction estimate
5118 if (dirm<0) continue;
5119 Float_t mpt = TMath::Sqrt(kink->GetMotherP()[0]*kink->GetMotherP()[0]+kink->GetMotherP()[1]*kink->GetMotherP()[1]);
5120 if (mpt<0.2) continue;
5123 //for high momenta momentum not defined well in first iteration
5124 Double_t qt = TMath::Sin(kink->GetAngle(2))*ktrack1->GetP();
5125 if (qt>0.35) continue;
5128 kink->SetLabel(CookLabel(ktrack0,0.4,0,row0),0);
5129 kink->SetLabel(CookLabel(ktrack1,0.4,row0,160),1);
5131 kink->SetTPCDensity(dens00,0,0);
5132 kink->SetTPCDensity(dens01,0,1);
5133 kink->SetTPCDensity(dens10,1,0);
5134 kink->SetTPCDensity(dens11,1,1);
5135 kink->SetIndex(i,0);
5136 kink->SetIndex(j,1);
5139 kink->SetTPCDensity(dens10,0,0);
5140 kink->SetTPCDensity(dens11,0,1);
5141 kink->SetTPCDensity(dens00,1,0);
5142 kink->SetTPCDensity(dens01,1,1);
5143 kink->SetIndex(j,0);
5144 kink->SetIndex(i,1);
5147 if (mpt<1||kink->GetAngle(2)>0.1){
5148 // angle and densities not defined yet
5149 if (kink->GetTPCDensityFactor()<0.8) continue;
5150 if ((2-kink->GetTPCDensityFactor())*kink->GetDistance() >0.25) continue;
5151 if (kink->GetAngle(2)*ktrack0->GetP()<0.003) continue; //too small angle
5152 if (kink->GetAngle(2)>0.2&&kink->GetTPCDensityFactor()<1.15) continue;
5153 if (kink->GetAngle(2)>0.2&&kink->GetTPCDensity(0,1)>0.05) continue;
5155 Float_t criticalangle = track0->GetSigmaSnp2()+track0->GetSigmaTgl2();
5156 criticalangle+= track1->GetSigmaSnp2()+track1->GetSigmaTgl2();
5157 criticalangle= 3*TMath::Sqrt(criticalangle);
5158 if (criticalangle>0.02) criticalangle=0.02;
5159 if (kink->GetAngle(2)<criticalangle) continue;
5162 Int_t drow = Int_t(2.+0.5/(0.05+kink->GetAngle(2))); // overlap region defined
5163 Float_t shapesum =0;
5165 for ( Int_t row = row0-drow; row<row0+drow;row++){
5166 if (row<0) continue;
5167 if (row>155) continue;
5168 if (ktrack0->GetClusterPointer(row)){
5169 AliTPCTrackerPoint *point =ktrack0->GetTrackPoint(row);
5170 shapesum+=point->GetSigmaY()+point->GetSigmaZ();
5173 if (ktrack1->GetClusterPointer(row)){
5174 AliTPCTrackerPoint *point =ktrack1->GetTrackPoint(row);
5175 shapesum+=point->GetSigmaY()+point->GetSigmaZ();
5180 kink->SetShapeFactor(-1.);
5183 kink->SetShapeFactor(shapesum/sum);
5185 // esd->AddKink(kink);
5187 // kink->SetMother(paramm);
5188 //kink->SetDaughter(paramd);
5190 Double_t chi2P2 = paramm.GetParameter()[2]-paramd.GetParameter()[2];
5192 chi2P2/=paramm.GetCovariance()[5]+paramd.GetCovariance()[5];
5193 Double_t chi2P3 = paramm.GetParameter()[3]-paramd.GetParameter()[3];
5195 chi2P3/=paramm.GetCovariance()[9]+paramd.GetCovariance()[9];
5197 if (AliTPCReconstructor::StreamLevel()>1) {
5198 (*fDebugStreamer)<<"kinkLpt"<<
5206 if ( chi2P2+chi2P3<AliTPCReconstructor::GetRecoParam()->GetKinkAngleCutChi2(0)){
5210 kinks->AddLast(kink);
5216 // sort the kinks according quality - and refit them towards vertex
5218 Int_t nkinks = kinks->GetEntriesFast();
5219 Float_t *quality = new Float_t[nkinks];
5220 Int_t *indexes = new Int_t[nkinks];
5221 AliTPCseed *mothers = new AliTPCseed[nkinks];
5222 AliTPCseed *daughters = new AliTPCseed[nkinks];
5225 for (Int_t i=0;i<nkinks;i++){
5227 AliKink *kinkl = (AliKink*)kinks->At(i);
5229 // refit kinks towards vertex
5231 Int_t index0 = kinkl->GetIndex(0);
5232 Int_t index1 = kinkl->GetIndex(1);
5233 AliTPCseed * ktrack0 = (AliTPCseed*)array->At(index0);
5234 AliTPCseed * ktrack1 = (AliTPCseed*)array->At(index1);
5236 Int_t sumn=ktrack0->GetNumberOfClusters()+ktrack1->GetNumberOfClusters();
5238 // Refit Kink under if too small angle
5240 if (kinkl->GetAngle(2)<0.05){
5241 kinkl->SetTPCRow0(GetRowNumber(kinkl->GetR()));
5242 Int_t row0 = kinkl->GetTPCRow0();
5243 Int_t drow = Int_t(2.+0.5/(0.05+kinkl->GetAngle(2)));
5246 Int_t last = row0-drow;
5247 if (last<40) last=40;
5248 if (last<ktrack0->GetFirstPoint()+25) last = ktrack0->GetFirstPoint()+25;
5249 AliTPCseed* seed0 = ReSeed(ktrack0,last,kFALSE);
5252 Int_t first = row0+drow;
5253 if (first>130) first=130;
5254 if (first>ktrack1->GetLastPoint()-25) first = TMath::Max(ktrack1->GetLastPoint()-25,30);
5255 AliTPCseed* seed1 = ReSeed(ktrack1,first,kTRUE);
5257 if (seed0 && seed1){
5258 kinkl->SetStatus(1,8);
5259 if (RefitKink(*seed0,*seed1,*kinkl)) kinkl->SetStatus(1,9);
5260 row0 = GetRowNumber(kinkl->GetR());
5261 sumn = seed0->GetNumberOfClusters()+seed1->GetNumberOfClusters();
5262 mothers[i] = *seed0;
5263 daughters[i] = *seed1;
5266 delete kinks->RemoveAt(i);
5267 if (seed0) delete seed0;
5268 if (seed1) delete seed1;
5271 if (kinkl->GetDistance()>0.5 || kinkl->GetR()<110 || kinkl->GetR()>240) {
5272 delete kinks->RemoveAt(i);
5273 if (seed0) delete seed0;
5274 if (seed1) delete seed1;
5282 if (kinkl) quality[i] = 160*((0.1+kinkl->GetDistance())*(2.-kinkl->GetTPCDensityFactor()))/(sumn+40.); //the longest -clossest will win
5284 TMath::Sort(nkinks,quality,indexes,kFALSE);
5286 //remove double find kinks
5288 for (Int_t ikink0=1;ikink0<nkinks;ikink0++){
5289 AliKink * kink0 = (AliKink*) kinks->At(indexes[ikink0]);
5290 if (!kink0) continue;
5292 for (Int_t ikink1=0;ikink1<ikink0;ikink1++){
5293 kink0 = (AliKink*) kinks->At(indexes[ikink0]);
5294 if (!kink0) continue;
5295 AliKink * kink1 = (AliKink*) kinks->At(indexes[ikink1]);
5296 if (!kink1) continue;
5297 // if not close kink continue
5298 if (TMath::Abs(kink1->GetPosition()[2]-kink0->GetPosition()[2])>10) continue;
5299 if (TMath::Abs(kink1->GetPosition()[1]-kink0->GetPosition()[1])>10) continue;
5300 if (TMath::Abs(kink1->GetPosition()[0]-kink0->GetPosition()[0])>10) continue;
5302 AliTPCseed &mother0 = mothers[indexes[ikink0]];
5303 AliTPCseed &daughter0 = daughters[indexes[ikink0]];
5304 AliTPCseed &mother1 = mothers[indexes[ikink1]];
5305 AliTPCseed &daughter1 = daughters[indexes[ikink1]];
5306 Int_t row0 = (kink0->GetTPCRow0()+kink1->GetTPCRow0())/2;
5315 for (Int_t i=0;i<row0;i++){
5316 if (mother0.GetClusterIndex(i)>0 && mother1.GetClusterIndex(i)>0){
5319 if (mother0.GetClusterIndex(i)==mother1.GetClusterIndex(i)){
5326 for (Int_t i=row0;i<158;i++){
5327 if (daughter0.GetClusterIndex(i)>0 && daughter0.GetClusterIndex(i)>0){
5330 if (mother0.GetClusterIndex(i)==mother1.GetClusterIndex(i)){
5336 Float_t ratio = Float_t(same+1)/Float_t(both+1);
5337 Float_t ratiom = Float_t(samem+1)/Float_t(bothm+1);
5338 Float_t ratiod = Float_t(samed+1)/Float_t(bothd+1);
5339 if (ratio>0.3 && ratiom>0.5 &&ratiod>0.5) {
5340 Int_t sum0 = mother0.GetNumberOfClusters()+daughter0.GetNumberOfClusters();
5341 Int_t sum1 = mother1.GetNumberOfClusters()+daughter1.GetNumberOfClusters();
5343 shared[kink0->GetIndex(0)]= kTRUE;
5344 shared[kink0->GetIndex(1)]= kTRUE;
5345 delete kinks->RemoveAt(indexes[ikink0]);
5349 shared[kink1->GetIndex(0)]= kTRUE;
5350 shared[kink1->GetIndex(1)]= kTRUE;
5351 delete kinks->RemoveAt(indexes[ikink1]);
5358 for (Int_t i=0;i<nkinks;i++){
5359 AliKink * kinkl = (AliKink*) kinks->At(indexes[i]);
5360 if (!kinkl) continue;
5361 kinkl->SetTPCRow0(GetRowNumber(kinkl->GetR()));
5362 Int_t index0 = kinkl->GetIndex(0);
5363 Int_t index1 = kinkl->GetIndex(1);
5364 if (circular[index0]||(circular[index1]&&kinkl->GetDistance()>0.2)) continue;
5365 kinkl->SetMultiple(usage[index0],0);
5366 kinkl->SetMultiple(usage[index1],1);
5367 if (kinkl->GetMultiple()[0]+kinkl->GetMultiple()[1]>2) continue;
5368 if (kinkl->GetMultiple()[0]+kinkl->GetMultiple()[1]>0 && quality[indexes[i]]>0.2) continue;
5369 if (kinkl->GetMultiple()[0]+kinkl->GetMultiple()[1]>0 && kinkl->GetDistance()>0.2) continue;
5370 if (circular[index0]||(circular[index1]&&kinkl->GetDistance()>0.1)) continue;
5372 AliTPCseed * ktrack0 = (AliTPCseed*)array->At(index0);
5373 AliTPCseed * ktrack1 = (AliTPCseed*)array->At(index1);
5374 if (!ktrack0 || !ktrack1) continue;
5375 Int_t index = esd->AddKink(kinkl);
5378 if ( ktrack0->GetKinkIndex(0)==0 && ktrack1->GetKinkIndex(0)==0) { //best kink
5379 if (mothers[indexes[i]].GetNumberOfClusters()>20 && daughters[indexes[i]].GetNumberOfClusters()>20 && (mothers[indexes[i]].GetNumberOfClusters()+daughters[indexes[i]].GetNumberOfClusters())>100){
5380 *ktrack0 = mothers[indexes[i]];
5381 *ktrack1 = daughters[indexes[i]];
5385 ktrack0->SetKinkIndex(usage[index0],-(index+1));
5386 ktrack1->SetKinkIndex(usage[index1], (index+1));
5391 // Remove tracks corresponding to shared kink's
5393 for (Int_t i=0;i<nentries;i++){
5394 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5395 if (!track0) continue;
5396 if (track0->GetKinkIndex(0)!=0) continue;
5397 if (shared[i]) delete array->RemoveAt(i);
5402 RemoveUsed2(array,0.5,0.4,30);
5404 for (Int_t i=0;i<nentries;i++){
5405 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5406 if (!track0) continue;
5407 track0->CookdEdx(0.02,0.6);
5411 for (Int_t i=0;i<nentries;i++){
5412 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5413 if (!track0) continue;
5414 if (track0->Pt()<1.4) continue;
5415 //remove double high momenta tracks - overlapped with kink candidates
5418 for (Int_t icl=track0->GetFirstPoint();icl<track0->GetLastPoint(); icl++){
5419 if (track0->GetClusterPointer(icl)!=0){
5421 if (track0->GetClusterPointer(icl)->IsUsed(10)) ishared++;
5424 if (Float_t(ishared+1)/Float_t(all+1)>0.5) {
5425 delete array->RemoveAt(i);
5429 if (track0->GetKinkIndex(0)!=0) continue;
5430 if (track0->GetNumberOfClusters()<80) continue;
5432 AliTPCseed *pmother = new AliTPCseed();
5433 AliTPCseed *pdaughter = new AliTPCseed();
5434 AliKink *pkink = new AliKink;
5436 AliTPCseed & mother = *pmother;
5437 AliTPCseed & daughter = *pdaughter;
5438 AliKink & kinkl = *pkink;
5439 if (CheckKinkPoint(track0,mother,daughter, kinkl)){
5440 if (mother.GetNumberOfClusters()<30||daughter.GetNumberOfClusters()<20) {
5444 continue; //too short tracks
5446 if (mother.Pt()<1.4) {
5452 Int_t row0= kinkl.GetTPCRow0();
5453 if (kinkl.GetDistance()>0.5 || kinkl.GetR()<110. || kinkl.GetR()>240.) {
5460 Int_t index = esd->AddKink(&kinkl);
5461 mother.SetKinkIndex(0,-(index+1));
5462 daughter.SetKinkIndex(0,index+1);
5463 if (mother.GetNumberOfClusters()>50) {
5464 delete array->RemoveAt(i);
5465 array->AddAt(new AliTPCseed(mother),i);
5468 array->AddLast(new AliTPCseed(mother));
5470 array->AddLast(new AliTPCseed(daughter));
5471 for (Int_t icl=0;icl<row0;icl++) {
5472 if (mother.GetClusterPointer(icl)) mother.GetClusterPointer(icl)->Use(20);
5475 for (Int_t icl=row0;icl<158;icl++) {
5476 if (daughter.GetClusterPointer(icl)) daughter.GetClusterPointer(icl)->Use(20);
5485 delete [] daughters;
5507 AliInfo(Form("Ncandidates=\t%d\t%d\t%d\t%d\n",esd->GetNumberOfKinks(),ncandidates,ntracks,nall));
5512 Int_t AliTPCtrackerMI::RefitKink(AliTPCseed &mother, AliTPCseed &daughter, const AliESDkink &knk)
5515 // refit kink towards to the vertex
5518 AliKink &kink=(AliKink &)knk;
5520 Int_t row0 = GetRowNumber(kink.GetR());
5521 FollowProlongation(mother,0);
5522 mother.Reset(kFALSE);
5524 FollowProlongation(daughter,row0);
5525 daughter.Reset(kFALSE);
5526 FollowBackProlongation(daughter,158);
5527 daughter.Reset(kFALSE);
5528 Int_t first = TMath::Max(row0-20,30);
5529 Int_t last = TMath::Min(row0+20,140);
5531 const Int_t kNdiv =5;
5532 AliTPCseed param0[kNdiv]; // parameters along the track
5533 AliTPCseed param1[kNdiv]; // parameters along the track
5534 AliKink kinks[kNdiv]; // corresponding kink parameters
5537 for (Int_t irow=0; irow<kNdiv;irow++){
5538 rows[irow] = first +((last-first)*irow)/(kNdiv-1);
5540 // store parameters along the track
5542 for (Int_t irow=0;irow<kNdiv;irow++){
5543 FollowBackProlongation(mother, rows[irow]);
5544 FollowProlongation(daughter,rows[kNdiv-1-irow]);
5545 param0[irow] = mother;
5546 param1[kNdiv-1-irow] = daughter;
5550 for (Int_t irow=0; irow<kNdiv-1;irow++){
5551 if (param0[irow].GetNumberOfClusters()<kNdiv||param1[irow].GetNumberOfClusters()<kNdiv) continue;
5552 kinks[irow].SetMother(param0[irow]);
5553 kinks[irow].SetDaughter(param1[irow]);
5554 kinks[irow].Update();
5557 // choose kink with best "quality"
5559 Double_t mindist = 10000;
5560 for (Int_t irow=0;irow<kNdiv;irow++){
5561 if (param0[irow].GetNumberOfClusters()<20||param1[irow].GetNumberOfClusters()<20) continue;
5562 if (TMath::Abs(kinks[irow].GetR())>240.) continue;
5563 if (TMath::Abs(kinks[irow].GetR())<100.) continue;
5565 Float_t normdist = TMath::Abs(param0[irow].GetX()-kinks[irow].GetR())*(0.1+kink.GetDistance());
5566 normdist/= (param0[irow].GetNumberOfClusters()+param1[irow].GetNumberOfClusters()+40.);
5567 if (normdist < mindist){
5573 if (index==-1) return 0;
5576 param0[index].Reset(kTRUE);
5577 FollowProlongation(param0[index],0);
5579 mother = param0[index];
5580 daughter = param1[index]; // daughter in vertex
5582 kink.SetMother(mother);
5583 kink.SetDaughter(daughter);
5585 kink.SetTPCRow0(GetRowNumber(kink.GetR()));
5586 kink.SetTPCncls(param0[index].GetNumberOfClusters(),0);
5587 kink.SetTPCncls(param1[index].GetNumberOfClusters(),1);
5588 kink.SetLabel(CookLabel(&mother,0.4, 0,kink.GetTPCRow0()),0);
5589 kink.SetLabel(CookLabel(&daughter,0.4, kink.GetTPCRow0(),160),1);
5590 mother.SetLabel(kink.GetLabel(0));
5591 daughter.SetLabel(kink.GetLabel(1));
5597 void AliTPCtrackerMI::UpdateKinkQualityM(AliTPCseed * seed){
5599 // update Kink quality information for mother after back propagation
5601 if (seed->GetKinkIndex(0)>=0) return;
5602 for (Int_t ikink=0;ikink<3;ikink++){
5603 Int_t index = seed->GetKinkIndex(ikink);
5604 if (index>=0) break;
5605 index = TMath::Abs(index)-1;
5606 AliESDkink * kink = fEvent->GetKink(index);
5607 kink->SetTPCDensity(-1,0,0);
5608 kink->SetTPCDensity(1,0,1);
5610 Int_t row0 = kink->GetTPCRow0() - 2 - Int_t( 0.5/ (0.05+kink->GetAngle(2)));
5611 if (row0<15) row0=15;
5613 Int_t row1 = kink->GetTPCRow0() + 2 + Int_t( 0.5/ (0.05+kink->GetAngle(2)));
5614 if (row1>145) row1=145;
5616 Int_t found,foundable,shared;
5617 seed->GetClusterStatistic(0,row0, found, foundable,shared,kFALSE);
5618 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),0,0);
5619 seed->GetClusterStatistic(row1,155, found, foundable,shared,kFALSE);
5620 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),0,1);
5625 void AliTPCtrackerMI::UpdateKinkQualityD(AliTPCseed * seed){
5627 // update Kink quality information for daughter after refit
5629 if (seed->GetKinkIndex(0)<=0) return;
5630 for (Int_t ikink=0;ikink<3;ikink++){
5631 Int_t index = seed->GetKinkIndex(ikink);
5632 if (index<=0) break;
5633 index = TMath::Abs(index)-1;
5634 AliESDkink * kink = fEvent->GetKink(index);
5635 kink->SetTPCDensity(-1,1,0);
5636 kink->SetTPCDensity(-1,1,1);
5638 Int_t row0 = kink->GetTPCRow0() -2 - Int_t( 0.5/ (0.05+kink->GetAngle(2)));
5639 if (row0<15) row0=15;
5641 Int_t row1 = kink->GetTPCRow0() +2 + Int_t( 0.5/ (0.05+kink->GetAngle(2)));
5642 if (row1>145) row1=145;
5644 Int_t found,foundable,shared;
5645 seed->GetClusterStatistic(0,row0, found, foundable,shared,kFALSE);
5646 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),1,0);
5647 seed->GetClusterStatistic(row1,155, found, foundable,shared,kFALSE);
5648 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),1,1);
5654 Int_t AliTPCtrackerMI::CheckKinkPoint(AliTPCseed*seed,AliTPCseed &mother, AliTPCseed &daughter, const AliESDkink &knk)
5657 // check kink point for given track
5658 // if return value=0 kink point not found
5659 // otherwise seed0 correspond to mother particle
5660 // seed1 correspond to daughter particle
5661 // kink parameter of kink point
5662 AliKink &kink=(AliKink &)knk;
5664 Int_t middlerow = (seed->GetFirstPoint()+seed->GetLastPoint())/2;
5665 Int_t first = seed->GetFirstPoint();
5666 Int_t last = seed->GetLastPoint();
5667 if (last-first<20) return 0; // shortest length - 2*30 = 60 pad-rows
5670 AliTPCseed *seed1 = ReSeed(seed,middlerow+20, kTRUE); //middle of chamber
5671 if (!seed1) return 0;
5672 FollowProlongation(*seed1,seed->GetLastPoint()-20);
5673 seed1->Reset(kTRUE);
5674 FollowProlongation(*seed1,158);
5675 seed1->Reset(kTRUE);
5676 last = seed1->GetLastPoint();
5678 AliTPCseed *seed0 = new AliTPCseed(*seed);
5679 seed0->Reset(kFALSE);
5682 AliTPCseed param0[20]; // parameters along the track
5683 AliTPCseed param1[20]; // parameters along the track
5684 AliKink kinks[20]; // corresponding kink parameters
5686 for (Int_t irow=0; irow<20;irow++){
5687 rows[irow] = first +((last-first)*irow)/19;
5689 // store parameters along the track
5691 for (Int_t irow=0;irow<20;irow++){
5692 FollowBackProlongation(*seed0, rows[irow]);
5693 FollowProlongation(*seed1,rows[19-irow]);
5694 param0[irow] = *seed0;
5695 param1[19-irow] = *seed1;
5699 for (Int_t irow=0; irow<19;irow++){
5700 kinks[irow].SetMother(param0[irow]);
5701 kinks[irow].SetDaughter(param1[irow]);
5702 kinks[irow].Update();
5705 // choose kink with biggest change of angle
5707 Double_t maxchange= 0;
5708 for (Int_t irow=1;irow<19;irow++){
5709 if (TMath::Abs(kinks[irow].GetR())>240.) continue;
5710 if (TMath::Abs(kinks[irow].GetR())<110.) continue;
5711 Float_t quality = TMath::Abs(kinks[irow].GetAngle(2))/(3.+TMath::Abs(kinks[irow].GetR()-param0[irow].GetX()));
5712 if ( quality > maxchange){
5713 maxchange = quality;
5720 if (index<0) return 0;
5722 Int_t row0 = GetRowNumber(kinks[index].GetR()); //row 0 estimate
5723 seed0 = new AliTPCseed(param0[index]);
5724 seed1 = new AliTPCseed(param1[index]);
5725 seed0->Reset(kFALSE);
5726 seed1->Reset(kFALSE);
5727 seed0->ResetCovariance(10.);
5728 seed1->ResetCovariance(10.);
5729 FollowProlongation(*seed0,0);
5730 FollowBackProlongation(*seed1,158);
5731 mother = *seed0; // backup mother at position 0
5732 seed0->Reset(kFALSE);
5733 seed1->Reset(kFALSE);
5734 seed0->ResetCovariance(10.);
5735 seed1->ResetCovariance(10.);
5737 first = TMath::Max(row0-20,0);
5738 last = TMath::Min(row0+20,158);
5740 for (Int_t irow=0; irow<20;irow++){
5741 rows[irow] = first +((last-first)*irow)/19;
5743 // store parameters along the track
5745 for (Int_t irow=0;irow<20;irow++){
5746 FollowBackProlongation(*seed0, rows[irow]);
5747 FollowProlongation(*seed1,rows[19-irow]);
5748 param0[irow] = *seed0;
5749 param1[19-irow] = *seed1;
5753 for (Int_t irow=0; irow<19;irow++){
5754 kinks[irow].SetMother(param0[irow]);
5755 kinks[irow].SetDaughter(param1[irow]);
5756 // param0[irow].Dump();
5757 //param1[irow].Dump();
5758 kinks[irow].Update();
5761 // choose kink with biggest change of angle
5764 for (Int_t irow=0;irow<20;irow++){
5765 if (TMath::Abs(kinks[irow].GetR())>250.) continue;
5766 if (TMath::Abs(kinks[irow].GetR())<90.) continue;
5767 Float_t quality = TMath::Abs(kinks[irow].GetAngle(2))/(3.+TMath::Abs(kinks[irow].GetR()-param0[irow].GetX()));
5768 if ( quality > maxchange){
5769 maxchange = quality;
5776 if (index==-1 || param0[index].GetNumberOfClusters()+param1[index].GetNumberOfClusters()<100){
5782 // Float_t anglesigma = TMath::Sqrt(param0[index].fC22+param0[index].fC33+param1[index].fC22+param1[index].fC33);
5784 kink.SetMother(param0[index]);
5785 kink.SetDaughter(param1[index]);
5788 Double_t chi2P2 = param0[index].GetParameter()[2]-param1[index].GetParameter()[2];
5790 chi2P2/=param0[index].GetCovariance()[5]+param1[index].GetCovariance()[5];
5791 Double_t chi2P3 = param0[index].GetParameter()[3]-param1[index].GetParameter()[3];
5793 chi2P3/=param0[index].GetCovariance()[9]+param1[index].GetCovariance()[9];
5795 if (AliTPCReconstructor::StreamLevel()>1) {
5796 (*fDebugStreamer)<<"kinkHpt"<<
5799 "p0.="<<¶m0[index]<<
5800 "p1.="<<¶m1[index]<<
5804 if ( chi2P2+chi2P3<AliTPCReconstructor::GetRecoParam()->GetKinkAngleCutChi2(0)){
5811 row0 = GetRowNumber(kink.GetR());
5812 kink.SetTPCRow0(row0);
5813 kink.SetLabel(CookLabel(seed0,0.5,0,row0),0);
5814 kink.SetLabel(CookLabel(seed1,0.5,row0,158),1);
5815 kink.SetIndex(-10,0);
5816 kink.SetIndex(int(param0[index].GetNumberOfClusters()+param1[index].GetNumberOfClusters()),1);
5817 kink.SetTPCncls(param0[index].GetNumberOfClusters(),0);
5818 kink.SetTPCncls(param1[index].GetNumberOfClusters(),1);
5821 // new (&mother) AliTPCseed(param0[index]);
5822 daughter = param1[index];
5823 daughter.SetLabel(kink.GetLabel(1));
5824 param0[index].Reset(kTRUE);
5825 FollowProlongation(param0[index],0);
5826 mother = param0[index];
5827 mother.SetLabel(kink.GetLabel(0));
5828 if ( chi2P2+chi2P3<AliTPCReconstructor::GetRecoParam()->GetKinkAngleCutChi2(1)){
5840 AliTPCseed* AliTPCtrackerMI::ReSeed(AliTPCseed *t)
5843 // reseed - refit - track
5846 // Int_t last = fSectors->GetNRows()-1;
5848 if (fSectors == fOuterSec){
5849 first = TMath::Max(first, t->GetFirstPoint()-fInnerSec->GetNRows());
5853 first = t->GetFirstPoint();
5855 AliTPCseed * seed = MakeSeed(t,0.1,0.5,0.9);
5856 FollowBackProlongation(*t,fSectors->GetNRows()-1);
5858 FollowProlongation(*t,first);
5868 //_____________________________________________________________________________
5869 Int_t AliTPCtrackerMI::ReadSeeds(const TFile *inp) {
5870 //-----------------------------------------------------------------
5871 // This function reades track seeds.
5872 //-----------------------------------------------------------------
5873 TDirectory *savedir=gDirectory;
5875 TFile *in=(TFile*)inp;
5876 if (!in->IsOpen()) {
5877 cerr<<"AliTPCtrackerMI::ReadSeeds(): input file is not open !\n";
5882 TTree *seedTree=(TTree*)in->Get("Seeds");
5884 cerr<<"AliTPCtrackerMI::ReadSeeds(): ";
5885 cerr<<"can't get a tree with track seeds !\n";
5888 AliTPCtrack *seed=new AliTPCtrack;
5889 seedTree->SetBranchAddress("tracks",&seed);
5891 if (fSeeds==0) fSeeds=new TObjArray(15000);
5893 Int_t n=(Int_t)seedTree->GetEntries();
5894 for (Int_t i=0; i<n; i++) {
5895 seedTree->GetEvent(i);
5896 fSeeds->AddLast(new AliTPCseed(*seed/*,seed->GetAlpha()*/));
5905 Int_t AliTPCtrackerMI::Clusters2Tracks (AliESDEvent *const esd)
5908 // clusters to tracks
5910 if (fSeeds) DeleteSeeds();
5913 AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ;
5914 transform->SetCurrentTimeStamp( esd->GetTimeStamp());
5915 transform->SetCurrentRun(esd->GetRunNumber());
5918 if (!fSeeds) return 1;
5920 if (AliTPCReconstructor::StreamLevel()>3) DumpClusters(0,fSeeds);
5926 //_____________________________________________________________________________
5927 Int_t AliTPCtrackerMI::Clusters2Tracks() {
5928 //-----------------------------------------------------------------
5929 // This is a track finder.
5930 //-----------------------------------------------------------------
5931 TDirectory *savedir=gDirectory;
5935 fSeeds = Tracking();
5938 Info("Clusters2Tracks","Time for tracking: \t");timer.Print();timer.Start();
5940 //activate again some tracks
5941 for (Int_t i=0; i<fSeeds->GetEntriesFast(); i++) {
5942 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
5944 Int_t nc=t.GetNumberOfClusters();
5946 delete fSeeds->RemoveAt(i);
5950 if (pt->GetRemoval()==10) {
5951 if (pt->GetDensityFirst(20)>0.8 || pt->GetDensityFirst(30)>0.8 || pt->GetDensityFirst(40)>0.7)
5952 pt->Desactivate(10); // make track again active // MvL: should be 0 ?
5954 pt->Desactivate(20);
5955 delete fSeeds->RemoveAt(i);
5960 RemoveUsed2(fSeeds,0.85,0.85,0);
5961 if (AliTPCReconstructor::GetRecoParam()->GetDoKinks()) FindKinks(fSeeds,fEvent);
5962 //FindCurling(fSeeds, fEvent,0);
5963 if (AliTPCReconstructor::StreamLevel()>5) FindMultiMC(fSeeds, fEvent,-1); // find multi found tracks
5964 RemoveUsed2(fSeeds,0.5,0.4,20);
5965 FindSplitted(fSeeds, fEvent,0); // find multi found tracks
5966 if (AliTPCReconstructor::StreamLevel()>5) FindMultiMC(fSeeds, fEvent,0); // find multi found tracks
5969 // // refit short tracks
5971 Int_t nseed=fSeeds->GetEntriesFast();
5974 for (Int_t i=0; i<nseed; i++) {
5975 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
5977 Int_t nc=t.GetNumberOfClusters();
5979 delete fSeeds->RemoveAt(i);
5982 CookLabel(pt,0.1); //For comparison only
5983 //if ((pt->IsActive() || (pt->fRemoval==10) )&& nc>50 &&pt->GetNumberOfClusters()>0.4*pt->fNFoundable){
5984 if ((pt->IsActive() || (pt->GetRemoval()==10) )){
5986 if (fDebug>0) cerr<<found<<'\r';
5990 delete fSeeds->RemoveAt(i);
5994 //RemoveOverlap(fSeeds,0.99,7,kTRUE);
5996 //RemoveUsed(fSeeds,0.9,0.9,6);
5998 nseed=fSeeds->GetEntriesFast();
6000 for (Int_t i=0; i<nseed; i++) {
6001 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
6003 Int_t nc=t.GetNumberOfClusters();
6005 delete fSeeds->RemoveAt(i);
6009 t.CookdEdx(0.02,0.6);
6010 // CheckKinkPoint(&t,0.05);
6011 //if ((pt->IsActive() || (pt->fRemoval==10) )&& nc>50 &&pt->GetNumberOfClusters()>0.4*pt->fNFoundable){
6012 if ((pt->IsActive() || (pt->GetRemoval()==10) )){
6020 delete fSeeds->RemoveAt(i);
6021 //AliTPCseed * seed1 = ReSeed(pt,0.05,0.5,1);
6023 // FollowProlongation(*seed1,0);
6024 // Int_t n = seed1->GetNumberOfClusters();
6025 // printf("fP4\t%f\t%f\n",seed1->GetC(),pt->GetC());
6026 // printf("fN\t%d\t%d\n", seed1->GetNumberOfClusters(),pt->GetNumberOfClusters());
6029 //AliTPCseed * seed2 = ReSeed(pt,0.95,0.5,0.05);
6033 SortTracks(fSeeds, 1);
6037 PrepareForBackProlongation(fSeeds,5.);
6038 PropagateBack(fSeeds);
6039 printf("Time for back propagation: \t");timer.Print();timer.Start();
6043 PrepareForProlongation(fSeeds,5.);
6044 PropagateForard2(fSeeds);
6046 printf("Time for FORWARD propagation: \t");timer.Print();timer.Start();
6047 // RemoveUsed(fSeeds,0.7,0.7,6);
6048 //RemoveOverlap(fSeeds,0.9,7,kTRUE);
6050 nseed=fSeeds->GetEntriesFast();
6052 for (Int_t i=0; i<nseed; i++) {
6053 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
6055 Int_t nc=t.GetNumberOfClusters();
6057 delete fSeeds->RemoveAt(i);
6060 t.CookdEdx(0.02,0.6);
6061 // CookLabel(pt,0.1); //For comparison only
6062 //if ((pt->IsActive() || (pt->fRemoval==10) )&& nc>50 &&pt->GetNumberOfClusters()>0.4*pt->fNFoundable){
6063 if ((pt->IsActive() || (pt->fRemoval==10) )){
6064 cerr<<found++<<'\r';
6067 delete fSeeds->RemoveAt(i);
6072 // fNTracks = found;
6074 Info("Clusters2Tracks","Time for overlap removal, track writing and dedx cooking: \t"); timer.Print();timer.Start();
6077 // cerr<<"Number of found tracks : "<<"\t"<<found<<endl;
6078 Info("Clusters2Tracks","Number of found tracks %d",found);
6080 // UnloadClusters();
6085 void AliTPCtrackerMI::Tracking(TObjArray * arr)
6088 // tracking of the seeds
6091 fSectors = fOuterSec;
6092 ParallelTracking(arr,150,63);
6093 fSectors = fOuterSec;
6094 ParallelTracking(arr,63,0);
6097 TObjArray * AliTPCtrackerMI::Tracking(Int_t seedtype, Int_t i1, Int_t i2, Float_t cuts[4], Float_t dy, Int_t dsec)
6102 TObjArray * arr = new TObjArray;
6104 fSectors = fOuterSec;
6107 for (Int_t sec=0;sec<fkNOS;sec++){
6108 if (seedtype==3) MakeSeeds3(arr,sec,i1,i2,cuts,dy, dsec);
6109 if (seedtype==4) MakeSeeds5(arr,sec,i1,i2,cuts,dy);
6110 if (seedtype==2) MakeSeeds2(arr,sec,i1,i2,cuts,dy);
6113 Info("Tracking","\nSeeding - %d\t%d\t%d\t%d\n",seedtype,i1,i2,arr->GetEntriesFast());
6125 TObjArray * AliTPCtrackerMI::Tracking()
6129 if (AliTPCReconstructor::GetRecoParam()->GetSpecialSeeding()) return TrackingSpecial();
6132 Int_t nup=fOuterSec->GetNRows()+fInnerSec->GetNRows();
6134 TObjArray * seeds = new TObjArray;
6136 Int_t fLastSeedRowSec=AliTPCReconstructor::GetRecoParam()->GetLastSeedRowSec();
6137 Int_t gapPrim = AliTPCReconstructor::GetRecoParam()->GetSeedGapPrim();
6138 Int_t gapSec = AliTPCReconstructor::GetRecoParam()->GetSeedGapSec();
6146 Float_t fnumber = 3.0;
6147 Float_t fdensity = 3.0;
6152 for (Int_t delta = 0; delta<18; delta+=gapPrim){
6156 arr = Tracking(3,nup-1-delta,nup-1-delta-gap,cuts,-1,1);
6157 SumTracks(seeds,arr);
6158 SignClusters(seeds,fnumber,fdensity);
6160 for (Int_t i=2;i<6;i+=2){
6161 // seed high pt tracks
6164 arr = Tracking(3,nup-i-delta,nup-i-delta-gap,cuts,-1,0);
6165 SumTracks(seeds,arr);
6166 SignClusters(seeds,fnumber,fdensity);
6171 // RemoveUsed(seeds,0.9,0.9,1);
6172 // UnsignClusters();
6173 // SignClusters(seeds,fnumber,fdensity);
6177 for (Int_t delta = 20; delta<120; delta+=gapPrim){
6179 // seed high pt tracks
6183 arr = Tracking(3,nup-delta,nup-delta-gap,cuts,-1);
6184 SumTracks(seeds,arr);
6185 SignClusters(seeds,fnumber,fdensity);
6190 arr = Tracking(3,nup-delta-5,nup-delta-5-gap,cuts,-1);
6191 SumTracks(seeds,arr);
6192 SignClusters(seeds,fnumber,fdensity);
6203 Info("Tracking()","\n\nPrimary seeding\t%d\n\n",seeds->GetEntriesFast());
6207 // RemoveUsed(seeds,0.75,0.75,1);
6209 //SignClusters(seeds,fnumber,fdensity);
6218 arr = Tracking(4,nup-1,nup-1-gap,cuts,-1);
6219 SumTracks(seeds,arr);
6220 SignClusters(seeds,fnumber,fdensity);
6222 arr = Tracking(4,nup-2,nup-2-gap,cuts,-1);
6223 SumTracks(seeds,arr);
6224 SignClusters(seeds,fnumber,fdensity);
6226 arr = Tracking(4,nup-3,nup-3-gap,cuts,-1);
6227 SumTracks(seeds,arr);
6228 SignClusters(seeds,fnumber,fdensity);
6230 arr = Tracking(4,nup-5,nup-5-gap,cuts,-1);
6231 SumTracks(seeds,arr);
6232 SignClusters(seeds,fnumber,fdensity);
6234 arr = Tracking(4,nup-7,nup-7-gap,cuts,-1);
6235 SumTracks(seeds,arr);
6236 SignClusters(seeds,fnumber,fdensity);
6239 arr = Tracking(4,nup-9,nup-9-gap,cuts,-1);
6240 SumTracks(seeds,arr);
6241 SignClusters(seeds,fnumber,fdensity);
6245 for (Int_t delta = 9; delta<30; delta+=gapSec){
6251 arr = Tracking(4,nup-1-delta,nup-1-delta-gap,cuts,-1);
6252 SumTracks(seeds,arr);
6253 SignClusters(seeds,fnumber,fdensity);
6255 arr = Tracking(4,nup-3-delta,nup-5-delta-gap,cuts,4);
6256 SumTracks(seeds,arr);
6257 SignClusters(seeds,fnumber,fdensity);
6270 for (Int_t delta = 30; delta<fLastSeedRowSec; delta+=gapSec){
6276 arr = Tracking(4,nup-1-delta,nup-1-delta-gap,cuts,-1);
6277 SumTracks(seeds,arr);
6278 SignClusters(seeds,fnumber,fdensity);
6280 arr = Tracking(4,nup-5-delta,nup-5-delta-gap,cuts,5 );
6281 SumTracks(seeds,arr);
6282 SignClusters(seeds,fnumber,fdensity);
6286 Info("Tracking()","\n\nSecondary seeding\t%d\n\n",seeds->GetEntriesFast());
6297 TObjArray * AliTPCtrackerMI::TrackingSpecial()
6300 // seeding adjusted for laser and cosmic tests - short tracks with big inclination angle
6301 // no primary vertex seeding tried
6305 Int_t nup=fOuterSec->GetNRows()+fInnerSec->GetNRows();
6307 TObjArray * seeds = new TObjArray;
6312 Float_t fnumber = 3.0;
6313 Float_t fdensity = 3.0;
6316 cuts[0] = AliTPCReconstructor::GetRecoParam()->GetMaxC(); // max curvature
6317 cuts[1] = 3.5; // max tan(phi) angle for seeding
6318 cuts[2] = 3.; // not used (cut on z primary vertex)
6319 cuts[3] = 3.5; // max tan(theta) angle for seeding
6321 for (Int_t delta = 0; nup-delta-gap-1>0; delta+=3){
6323 arr = Tracking(4,nup-1-delta,nup-1-delta-gap,cuts,-1);
6324 SumTracks(seeds,arr);
6325 SignClusters(seeds,fnumber,fdensity);
6329 Info("Tracking()","\n\nSecondary seeding\t%d\n\n",seeds->GetEntriesFast());
6340 void AliTPCtrackerMI::SumTracks(TObjArray *arr1,TObjArray *&arr2) const
6343 //sum tracks to common container
6344 //remove suspicious tracks
6345 Int_t nseed = arr2->GetEntriesFast();
6346 for (Int_t i=0;i<nseed;i++){
6347 AliTPCseed *pt=(AliTPCseed*)arr2->UncheckedAt(i);
6350 // remove tracks with too big curvature
6352 if (TMath::Abs(pt->GetC())>AliTPCReconstructor::GetRecoParam()->GetMaxC()){
6353 delete arr2->RemoveAt(i);
6356 // REMOVE VERY SHORT TRACKS
6357 if (pt->GetNumberOfClusters()<20){
6358 delete arr2->RemoveAt(i);
6361 // NORMAL ACTIVE TRACK
6362 if (pt->IsActive()){
6363 arr1->AddLast(arr2->RemoveAt(i));
6366 //remove not usable tracks
6367 if (pt->GetRemoval()!=10){
6368 delete arr2->RemoveAt(i);
6372 // ENABLE ONLY ENOUGH GOOD STOPPED TRACKS
6373 if (pt->GetDensityFirst(20)>0.8 || pt->GetDensityFirst(30)>0.8 || pt->GetDensityFirst(40)>0.7)
6374 arr1->AddLast(arr2->RemoveAt(i));
6376 delete arr2->RemoveAt(i);
6380 delete arr2; arr2 = 0;
6385 void AliTPCtrackerMI::ParallelTracking(TObjArray *const arr, Int_t rfirst, Int_t rlast)
6388 // try to track in parralel
6390 Int_t nseed=arr->GetEntriesFast();
6391 //prepare seeds for tracking
6392 for (Int_t i=0; i<nseed; i++) {
6393 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i), &t=*pt;
6395 if (!t.IsActive()) continue;
6396 // follow prolongation to the first layer
6397 if ( (fSectors ==fInnerSec) || (t.GetFirstPoint()-fkParam->GetNRowLow()>rfirst+1) )
6398 FollowProlongation(t, rfirst+1);
6403 for (Int_t nr=rfirst; nr>=rlast; nr--){
6404 if (nr<fInnerSec->GetNRows())
6405 fSectors = fInnerSec;
6407 fSectors = fOuterSec;
6408 // make indexes with the cluster tracks for given
6410 // find nearest cluster
6411 for (Int_t i=0; i<nseed; i++) {
6412 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i), &t=*pt;
6414 if (nr==80) pt->UpdateReference();
6415 if (!pt->IsActive()) continue;
6416 // if ( (fSectors ==fOuterSec) && (pt->fFirstPoint-fkParam->GetNRowLow())<nr) continue;
6417 if (pt->GetRelativeSector()>17) {
6420 UpdateClusters(t,nr);
6422 // prolonagate to the nearest cluster - if founded
6423 for (Int_t i=0; i<nseed; i++) {
6424 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
6426 if (!pt->IsActive()) continue;
6427 // if ((fSectors ==fOuterSec) && (pt->fFirstPoint-fkParam->GetNRowLow())<nr) continue;
6428 if (pt->GetRelativeSector()>17) {
6431 FollowToNextCluster(*pt,nr);
6436 void AliTPCtrackerMI::PrepareForBackProlongation(const TObjArray *const arr,Float_t fac) const
6440 // if we use TPC track itself we have to "update" covariance
6442 Int_t nseed= arr->GetEntriesFast();
6443 for (Int_t i=0;i<nseed;i++){
6444 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
6448 //rotate to current local system at first accepted point
6449 Int_t index = pt->GetClusterIndex2(pt->GetFirstPoint());
6450 Int_t sec = (index&0xff000000)>>24;
6452 Float_t angle1 = fInnerSec->GetAlpha()*sec+fInnerSec->GetAlphaShift();
6453 if (angle1>TMath::Pi())
6454 angle1-=2.*TMath::Pi();
6455 Float_t angle2 = pt->GetAlpha();
6457 if (TMath::Abs(angle1-angle2)>0.001){
6458 if (!pt->Rotate(angle1-angle2)) return;
6459 //angle2 = pt->GetAlpha();
6460 //pt->fRelativeSector = pt->GetAlpha()/fInnerSec->GetAlpha();
6461 //if (pt->GetAlpha()<0)
6462 // pt->fRelativeSector+=18;
6463 //sec = pt->fRelativeSector;
6472 void AliTPCtrackerMI::PrepareForProlongation(TObjArray *const arr, Float_t fac) const
6476 // if we use TPC track itself we have to "update" covariance
6478 Int_t nseed= arr->GetEntriesFast();
6479 for (Int_t i=0;i<nseed;i++){
6480 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
6483 pt->SetFirstPoint(pt->GetLastPoint());
6491 Int_t AliTPCtrackerMI::PropagateBack(const TObjArray *const arr)
6494 // make back propagation
6496 Int_t nseed= arr->GetEntriesFast();
6497 for (Int_t i=0;i<nseed;i++){
6498 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
6499 if (pt&& pt->GetKinkIndex(0)<=0) {
6500 //AliTPCseed *pt2 = new AliTPCseed(*pt);
6501 fSectors = fInnerSec;
6502 //FollowBackProlongation(*pt,fInnerSec->GetNRows()-1);
6503 //fSectors = fOuterSec;
6504 FollowBackProlongation(*pt,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1,1);
6505 //if (pt->GetNumberOfClusters()<(pt->fEsd->GetTPCclusters(0)) ){
6506 // Error("PropagateBack","Not prolonged track %d",pt->GetLabel());
6507 // FollowBackProlongation(*pt2,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1);
6510 if (pt&& pt->GetKinkIndex(0)>0) {
6511 AliESDkink * kink = fEvent->GetKink(pt->GetKinkIndex(0)-1);
6512 pt->SetFirstPoint(kink->GetTPCRow0());
6513 fSectors = fInnerSec;
6514 FollowBackProlongation(*pt,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1,1);
6522 Int_t AliTPCtrackerMI::PropagateForward2(const TObjArray *const arr)
6525 // make forward propagation
6527 Int_t nseed= arr->GetEntriesFast();
6529 for (Int_t i=0;i<nseed;i++){
6530 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
6532 FollowProlongation(*pt,0,1,1);
6541 Int_t AliTPCtrackerMI::PropagateForward()
6544 // propagate track forward
6546 Int_t nseed = fSeeds->GetEntriesFast();
6547 for (Int_t i=0;i<nseed;i++){
6548 AliTPCseed *pt = (AliTPCseed*)fSeeds->UncheckedAt(i);
6550 AliTPCseed &t = *pt;
6551 Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
6552 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
6553 if (alpha < 0. ) alpha += 2.*TMath::Pi();
6554 t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN);
6558 fSectors = fOuterSec;
6559 ParallelTracking(fSeeds,fOuterSec->GetNRows()+fInnerSec->GetNRows()-1,fInnerSec->GetNRows());
6560 fSectors = fInnerSec;
6561 ParallelTracking(fSeeds,fInnerSec->GetNRows()-1,0);
6570 Int_t AliTPCtrackerMI::PropagateBack(AliTPCseed *const pt, Int_t row0, Int_t row1)
6573 // make back propagation, in between row0 and row1
6577 fSectors = fInnerSec;
6580 if (row1<fSectors->GetNRows())
6583 r1 = fSectors->GetNRows()-1;
6585 if (row0<fSectors->GetNRows()&& r1>0 )
6586 FollowBackProlongation(*pt,r1);
6587 if (row1<=fSectors->GetNRows())
6590 r1 = row1 - fSectors->GetNRows();
6591 if (r1<=0) return 0;
6592 if (r1>=fOuterSec->GetNRows()) return 0;
6593 fSectors = fOuterSec;
6594 return FollowBackProlongation(*pt,r1);
6602 void AliTPCtrackerMI::GetShape(AliTPCseed * seed, Int_t row)
6604 // gets cluster shape
6606 AliTPCClusterParam * clparam = AliTPCcalibDB::Instance()->GetClusterParam();
6607 Float_t zdrift = TMath::Abs((fkParam->GetZLength(0)-TMath::Abs(seed->GetZ())));
6608 Int_t type = (seed->GetSector() < fkParam->GetNSector()/2) ? 0: (row>126) ? 1:2;
6609 Double_t angulary = seed->GetSnp();
6611 if (TMath::Abs(angulary)>AliTPCReconstructor::GetMaxSnpTracker()) {
6612 angulary = TMath::Sign(AliTPCReconstructor::GetMaxSnpTracker(),angulary);
6615 angulary = angulary*angulary/((1.-angulary)*(1.+angulary));
6616 Double_t angularz = seed->GetTgl()*seed->GetTgl()*(1.+angulary);
6618 Double_t sigmay = clparam->GetRMS0(0,type,zdrift,TMath::Sqrt(TMath::Abs(angulary)));
6619 Double_t sigmaz = clparam->GetRMS0(1,type,zdrift,TMath::Sqrt(TMath::Abs(angularz)));
6620 seed->SetCurrentSigmaY2(sigmay*sigmay);
6621 seed->SetCurrentSigmaZ2(sigmaz*sigmaz);
6622 // Float_t sd2 = TMath::Abs((fkParam->GetZLength(0)-TMath::Abs(seed->GetZ())))*fkParam->GetDiffL()*fkParam->GetDiffL();
6623 // // Float_t padlength = fkParam->GetPadPitchLength(seed->fSector);
6624 // Float_t padlength = GetPadPitchLength(row);
6626 // Float_t sresy = (seed->GetSector() < fkParam->GetNSector()/2) ? 0.2 :0.3;
6627 // seed->SetCurrentSigmaY2(sd2+padlength*padlength*angulary/12.+sresy*sresy);
6629 // Float_t sresz = fkParam->GetZSigma();
6630 // seed->SetCurrentSigmaZ2(sd2+padlength*padlength*angularz*angularz*(1+angulary)/12.+sresz*sresz);
6632 Float_t wy = GetSigmaY(seed);
6633 Float_t wz = GetSigmaZ(seed);
6636 if (TMath::Abs(wy/seed->fCurrentSigmaY2-1)>0.0001 || TMath::Abs(wz/seed->fCurrentSigmaZ2-1)>0.0001 ){
6637 printf("problem\n");
6644 //__________________________________________________________________________
6645 void AliTPCtrackerMI::CookLabel(AliKalmanTrack *tk, Float_t wrong) const {
6646 //--------------------------------------------------------------------
6647 //This function "cooks" a track label. If label<0, this track is fake.
6648 //--------------------------------------------------------------------
6649 AliTPCseed * t = dynamic_cast<AliTPCseed*>(tk);
6651 printf("%s:%d wrong type \n",(char*)__FILE__,__LINE__);
6655 Int_t noc=t->GetNumberOfClusters();
6657 //printf("\nnot founded prolongation\n\n\n");
6663 AliTPCclusterMI *clusters[160];
6665 for (Int_t i=0;i<160;i++) {
6672 for (i=0; i<160 && current<noc; i++) {
6674 Int_t index=t->GetClusterIndex2(i);
6675 if (index<=0) continue;
6676 if (index&0x8000) continue;
6678 //clusters[current]=GetClusterMI(index);
6679 if (t->GetClusterPointer(i)){
6680 clusters[current]=t->GetClusterPointer(i);
6686 Int_t lab=123456789;
6687 for (i=0; i<noc; i++) {
6688 AliTPCclusterMI *c=clusters[i];
6690 lab=TMath::Abs(c->GetLabel(0));
6692 for (j=0; j<noc; j++) if (lb[j]==lab || mx[j]==0) break;
6698 for (i=0; i<noc; i++) if (mx[i]>max) {max=mx[i]; lab=lb[i];}
6700 for (i=0; i<noc; i++) {
6701 AliTPCclusterMI *c=clusters[i];
6703 if (TMath::Abs(c->GetLabel(1)) == lab ||
6704 TMath::Abs(c->GetLabel(2)) == lab ) max++;
6706 if (noc<=0) { lab=-1; return;}
6707 if ((1.- Float_t(max)/(noc)) > wrong) lab=-lab;
6710 Int_t tail=Int_t(0.10*noc);
6713 for (i=1; i<160&&ind<tail; i++) {
6714 // AliTPCclusterMI *c=clusters[noc-i];
6715 AliTPCclusterMI *c=clusters[i];
6717 if (lab == TMath::Abs(c->GetLabel(0)) ||
6718 lab == TMath::Abs(c->GetLabel(1)) ||
6719 lab == TMath::Abs(c->GetLabel(2))) max++;
6722 if (max < Int_t(0.5*tail)) lab=-lab;
6729 //delete[] clusters;
6733 //__________________________________________________________________________
6734 Int_t AliTPCtrackerMI::CookLabel(AliTPCseed *const t, Float_t wrong,Int_t first, Int_t last) const {
6735 //--------------------------------------------------------------------
6736 //This function "cooks" a track label. If label<0, this track is fake.
6737 //--------------------------------------------------------------------
6738 Int_t noc=t->GetNumberOfClusters();
6740 //printf("\nnot founded prolongation\n\n\n");
6746 AliTPCclusterMI *clusters[160];
6748 for (Int_t i=0;i<160;i++) {
6755 for (i=0; i<160 && current<noc; i++) {
6756 if (i<first) continue;
6757 if (i>last) continue;
6758 Int_t index=t->GetClusterIndex2(i);
6759 if (index<=0) continue;
6760 if (index&0x8000) continue;
6762 //clusters[current]=GetClusterMI(index);
6763 if (t->GetClusterPointer(i)){
6764 clusters[current]=t->GetClusterPointer(i);
6769 //if (noc<5) return -1;
6770 Int_t lab=123456789;
6771 for (i=0; i<noc; i++) {
6772 AliTPCclusterMI *c=clusters[i];
6774 lab=TMath::Abs(c->GetLabel(0));
6776 for (j=0; j<noc; j++) if (lb[j]==lab || mx[j]==0) break;
6782 for (i=0; i<noc; i++) if (mx[i]>max) {max=mx[i]; lab=lb[i];}
6784 for (i=0; i<noc; i++) {
6785 AliTPCclusterMI *c=clusters[i];
6787 if (TMath::Abs(c->GetLabel(1)) == lab ||
6788 TMath::Abs(c->GetLabel(2)) == lab ) max++;
6790 if (noc<=0) { lab=-1; return -1;}
6791 if ((1.- Float_t(max)/(noc)) > wrong) lab=-lab;
6794 Int_t tail=Int_t(0.10*noc);
6797 for (i=1; i<160&&ind<tail; i++) {
6798 // AliTPCclusterMI *c=clusters[noc-i];
6799 AliTPCclusterMI *c=clusters[i];
6801 if (lab == TMath::Abs(c->GetLabel(0)) ||
6802 lab == TMath::Abs(c->GetLabel(1)) ||
6803 lab == TMath::Abs(c->GetLabel(2))) max++;
6806 if (max < Int_t(0.5*tail)) lab=-lab;
6809 // t->SetLabel(lab);
6813 //delete[] clusters;
6817 Int_t AliTPCtrackerMI::GetRowNumber(Double_t x[3]) const
6819 //return pad row number for given x vector
6820 Float_t phi = TMath::ATan2(x[1],x[0]);
6821 if(phi<0) phi=2.*TMath::Pi()+phi;
6822 // Get the local angle in the sector philoc
6823 const Float_t kRaddeg = 180/3.14159265358979312;
6824 Float_t phiangle = (Int_t (phi*kRaddeg/20.) + 0.5)*20./kRaddeg;
6825 Double_t localx = x[0]*TMath::Cos(phiangle)-x[1]*TMath::Sin(phiangle);
6826 return GetRowNumber(localx);
6831 void AliTPCtrackerMI::MakeESDBitmaps(AliTPCseed *t, AliESDtrack *esd)
6833 //-----------------------------------------------------------------------
6834 // Fill the cluster and sharing bitmaps of the track
6835 //-----------------------------------------------------------------------
6837 Int_t firstpoint = 0;
6838 Int_t lastpoint = 159;
6839 AliTPCTrackerPoint *point;
6840 AliTPCclusterMI *cluster;
6843 TBits clusterMap(159);
6844 TBits sharedMap(159);
6846 for (int iter=firstpoint; iter<lastpoint; iter++) {
6847 // Change to cluster pointers to see if we have a cluster at given padrow
6849 cluster = t->GetClusterPointer(iter);
6851 clusterMap.SetBitNumber(iter, kTRUE);
6852 point = t->GetTrackPoint(iter);
6853 if (point->IsShared())
6854 sharedMap.SetBitNumber(iter,kTRUE);
6856 if (t->GetClusterIndex(iter) >= 0 && (t->GetClusterIndex(iter) & 0x8000) == 0) {
6857 fitMap.SetBitNumber(iter, kTRUE);
6861 esd->SetTPCClusterMap(clusterMap);
6862 esd->SetTPCSharedMap(sharedMap);
6863 esd->SetTPCFitMap(fitMap);
6864 if (nclsf != t->GetNumberOfClusters())
6865 AliWarning(Form("Inconsistency between ncls %d and indices %d (found %d)",t->GetNumberOfClusters(),nclsf,esd->GetTPCClusterMap().CountBits()));
6868 Bool_t AliTPCtrackerMI::IsFindable(AliTPCseed & track){
6870 // return flag if there is findable cluster at given position
6873 Float_t z = track.GetZ();
6875 if (TMath::Abs(z)<(AliTPCReconstructor::GetCtgRange()*track.GetX()+kDeltaZ) &&
6876 TMath::Abs(z)<fkParam->GetZLength(0) &&
6877 (TMath::Abs(track.GetSnp())<AliTPCReconstructor::GetMaxSnpTracker()))
6883 void AliTPCtrackerMI::AddCovariance(AliTPCseed * seed){
6885 // Adding systematic error estimate to the covariance matrix
6886 // !!!! the systematic error for element 4 is in 1/cm not in pt
6888 const Double_t *param = AliTPCReconstructor::GetRecoParam()->GetSystematicError();
6890 // use only the diagonal part if not specified otherwise
6891 if (!AliTPCReconstructor::GetRecoParam()->GetUseSystematicCorrelation()) return AddCovarianceAdd(seed);
6893 Double_t *covarS= (Double_t*)seed->GetCovariance();
6894 Double_t factor[5]={1,1,1,1,1};
6895 Double_t facC = AliTracker::GetBz()*kB2C;
6896 factor[0]= TMath::Sqrt(TMath::Abs((covarS[0] + param[0]*param[0])/covarS[0]));
6897 factor[1]= TMath::Sqrt(TMath::Abs((covarS[2] + param[1]*param[1])/covarS[2]));
6898 factor[2]= TMath::Sqrt(TMath::Abs((covarS[5] + param[2]*param[2])/covarS[5]));
6899 factor[3]= TMath::Sqrt(TMath::Abs((covarS[9] + param[3]*param[3])/covarS[9]));
6900 factor[4]= TMath::Sqrt(TMath::Abs((covarS[14] + facC*facC*param[4]*param[4])/covarS[14]));
6906 for (Int_t i=0; i<5; i++){
6907 for (Int_t j=i; j<5; j++){
6908 Int_t index=seed->GetIndex(i,j);
6909 covarS[index]*=factor[i]*factor[j];
6915 void AliTPCtrackerMI::AddCovarianceAdd(AliTPCseed * seed){
6917 // Adding systematic error - as additive factor without correlation
6919 const Double_t *param = AliTPCReconstructor::GetRecoParam()->GetSystematicError();
6920 Double_t *covarIn= (Double_t*)seed->GetCovariance();
6922 for (Int_t i=0;i<15;i++) covar[i]=0;
6928 covar[0] = param[0]*param[0];
6929 covar[2] = param[1]*param[1];
6930 covar[5] = param[2]*param[2];
6931 covar[9] = param[3]*param[3];
6932 Double_t facC = AliTracker::GetBz()*kB2C;
6933 covar[14]= param[4]*param[4]*facC*facC;
6935 covar[1]=TMath::Sqrt((covar[0]*covar[2]))*covarIn[1]/TMath::Sqrt((covarIn[0]*covarIn[2]));
6937 covar[3]=TMath::Sqrt((covar[0]*covar[5]))*covarIn[3]/TMath::Sqrt((covarIn[0]*covarIn[5]));
6938 covar[4]=TMath::Sqrt((covar[2]*covar[5]))*covarIn[4]/TMath::Sqrt((covarIn[2]*covarIn[5]));
6940 covar[6]=TMath::Sqrt((covar[0]*covar[9]))*covarIn[6]/TMath::Sqrt((covarIn[0]*covarIn[9]));
6941 covar[7]=TMath::Sqrt((covar[2]*covar[9]))*covarIn[7]/TMath::Sqrt((covarIn[2]*covarIn[9]));
6942 covar[8]=TMath::Sqrt((covar[5]*covar[9]))*covarIn[8]/TMath::Sqrt((covarIn[5]*covarIn[9]));
6944 covar[10]=TMath::Sqrt((covar[0]*covar[14]))*covarIn[10]/TMath::Sqrt((covarIn[0]*covarIn[14]));
6945 covar[11]=TMath::Sqrt((covar[2]*covar[14]))*covarIn[11]/TMath::Sqrt((covarIn[2]*covarIn[14]));
6946 covar[12]=TMath::Sqrt((covar[5]*covar[14]))*covarIn[12]/TMath::Sqrt((covarIn[5]*covarIn[14]));
6947 covar[13]=TMath::Sqrt((covar[9]*covar[14]))*covarIn[13]/TMath::Sqrt((covarIn[9]*covarIn[14]));
6949 seed->AddCovariance(covar);