1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
17 //-------------------------------------------------------
18 // Implementation of the TPC tracker
20 // Origin: Marian Ivanov Marian.Ivanov@cern.ch
22 // AliTPC parallel tracker
24 // The track fitting is based on Kalaman filtering approach
26 // The track finding steps:
27 // 1. Seeding - with and without vertex constraint
28 // - seeding with vertex constain done at first n^2 proble
29 // - seeding without vertex constraint n^3 problem
30 // 2. Tracking - follow prolongation road - find cluster - update kalman track
32 // The seeding and tracking is repeated several times, in different seeding region.
33 // This approach enables to find the track which cannot be seeded in some region of TPC
34 // This can happen because of low momenta (track do not reach outer radius), or track is currently in the ded region between sectors, or the track is for the moment overlapped with other track (seed quality is poor) ...
36 // With this approach we reach almost 100 % efficiency also for high occupancy events.
37 // (If the seeding efficiency in a region is about 90 % than with logical or of several
38 // regions we will reach 100% (in theory - supposing independence)
40 // Repeating several seeding - tracking procedures some of the tracks can be find
43 // The procedures to remove multi find tacks are impremented:
44 // RemoveUsed2 - fast procedure n problem -
45 // Algorithm - Sorting tracks according quality
46 // remove tracks with some shared fraction
47 // Sharing in respect to all tacks
48 // Signing clusters in gold region
49 // FindSplitted - slower algorithm n^2
50 // Sort the tracks according quality
51 // Loop over pair of tracks
52 // If overlap with other track bigger than threshold - remove track
54 // FindCurling - Finds the pair of tracks which are curling
55 // - About 10% of tracks can be find with this procedure
56 // The combinatorial background is too big to be used in High
57 // multiplicity environment
58 // - n^2 problem - Slow procedure - currently it is disabled because of
61 // The number of splitted tracks can be reduced disabling the sharing of the cluster.
62 // tpcRecoParam-> SetClusterSharing(kFALSE);
63 // IT IS HIGHLY non recomended to use it in high flux enviroonment
64 // Even using this switch some tracks can be found more than once
65 // (because of multiple seeding and low quality tracks which will not cross full chamber)
68 // The tracker itself can be debugged - the information about tracks can be stored in several // phases of the reconstruction
69 // To enable storage of the TPC tracks in the ESD friend track
70 // use AliTPCReconstructor::SetStreamLevel(n);
72 // The debug level - different procedure produce tree for numerical debugging
73 // To enable them set AliTPCReconstructor::SetStreamLevel(n); where nis bigger 1
77 // Adding systematic errors to the covariance:
79 // The systematic errors due to the misalignment and miscalibration are added to the covariance matrix
80 // of the tracks (not to the clusters as they are dependent):
81 // The parameters form AliTPCRecoParam are used AliTPCRecoParam::GetSystematicError
82 // The systematic errors are expressed there in RMS - position (cm), angle (rad), curvature (1/cm)
83 // The default values are 0.
85 // The sytematic errors are added to the covariance matrix in following places:
87 // 1. During fisrt itteration - AliTPCtrackerMI::FillESD
88 // 2. Second iteration -
89 // 2.a ITS->TPC - AliTPCtrackerMI::ReadSeeds
90 // 2.b TPC->TRD - AliTPCtrackerMI::PropagateBack
91 // 3. Third iteration -
92 // 3.a TRD->TPC - AliTPCtrackerMI::ReadSeeds
93 // 3.b TPC->ITS - AliTPCtrackerMI::RefitInward
95 // There are several places in the code which can be numerically debuged
96 // This code is keeped in order to enable code development and to check the calibration implementtion
98 // 1. ErrParam stream - dump information about
100 // 2.a) cluster error estimate
101 // 3.a) cluster shape estimate
104 // Debug streamer levels:
106 //-------------------------------------------------------
111 #include "Riostream.h"
112 #include <TClonesArray.h>
114 #include <TObjArray.h>
116 #include <TGraphErrors.h>
118 #include "AliComplexCluster.h"
119 #include "AliESDEvent.h"
120 #include "AliESDtrack.h"
121 #include "AliESDVertex.h"
124 #include "AliHelix.h"
125 #include "AliRunLoader.h"
126 #include "AliTPCClustersRow.h"
127 #include "AliTPCParam.h"
128 #include "AliTPCReconstructor.h"
129 #include "AliTPCpolyTrack.h"
130 #include "AliTPCreco.h"
131 #include "AliTPCseed.h"
133 #include "AliTPCtrackerSector.h"
134 #include "AliTPCtrackerMI.h"
135 #include "TStopwatch.h"
136 #include "AliTPCReconstructor.h"
137 #include "AliAlignObj.h"
138 #include "AliTrackPointArray.h"
140 #include "AliTPCcalibDB.h"
141 #include "AliTPCcalibDButil.h"
142 #include "AliTPCTransform.h"
143 #include "AliTPCClusterParam.h"
147 ClassImp(AliTPCtrackerMI)
151 class AliTPCFastMath {
154 static Double_t FastAsin(Double_t x);
156 static Double_t fgFastAsin[20000]; //lookup table for fast asin computation
159 Double_t AliTPCFastMath::fgFastAsin[20000];
160 AliTPCFastMath gAliTPCFastMath; // needed to fill the LUT
162 AliTPCFastMath::AliTPCFastMath(){
164 // initialized lookup table;
165 for (Int_t i=0;i<10000;i++){
166 fgFastAsin[2*i] = TMath::ASin(i/10000.);
167 fgFastAsin[2*i+1] = (TMath::ASin((i+1)/10000.)-fgFastAsin[2*i]);
171 Double_t AliTPCFastMath::FastAsin(Double_t x){
173 // return asin using lookup table
175 Int_t index = int(x*10000);
176 return fgFastAsin[2*index]+(x*10000.-index)*fgFastAsin[2*index+1];
179 Int_t index = int(x*10000);
180 return -(fgFastAsin[2*index]+(x*10000.-index)*fgFastAsin[2*index+1]);
182 //__________________________________________________________________
183 AliTPCtrackerMI::AliTPCtrackerMI()
205 // default constructor
207 for (Int_t irow=0; irow<200; irow++){
214 //_____________________________________________________________________
218 Int_t AliTPCtrackerMI::UpdateTrack(AliTPCseed * track, Int_t accept){
220 //update track information using current cluster - track->fCurrentCluster
223 AliTPCclusterMI* c =track->GetCurrentCluster();
224 if (accept>0) track->SetCurrentClusterIndex1(track->GetCurrentClusterIndex1() | 0x8000); //sign not accepted clusters
226 UInt_t i = track->GetCurrentClusterIndex1();
228 Int_t sec=(i&0xff000000)>>24;
229 //Int_t row = (i&0x00ff0000)>>16;
230 track->SetRow((i&0x00ff0000)>>16);
231 track->SetSector(sec);
232 // Int_t index = i&0xFFFF;
233 if (sec>=fkParam->GetNInnerSector()) track->SetRow(track->GetRow()+fkParam->GetNRowLow());
234 track->SetClusterIndex2(track->GetRow(), i);
235 //track->fFirstPoint = row;
236 //if ( track->fLastPoint<row) track->fLastPoint =row;
237 // if (track->fRow<0 || track->fRow>160) {
238 // printf("problem\n");
240 if (track->GetFirstPoint()>track->GetRow())
241 track->SetFirstPoint(track->GetRow());
242 if (track->GetLastPoint()<track->GetRow())
243 track->SetLastPoint(track->GetRow());
246 track->SetClusterPointer(track->GetRow(),c);
249 Double_t angle2 = track->GetSnp()*track->GetSnp();
251 //SET NEW Track Point
253 if (angle2<1) //PH sometimes angle2 is very big. To be investigated...
255 angle2 = TMath::Sqrt(angle2/(1-angle2));
256 AliTPCTrackerPoint &point =*(track->GetTrackPoint(track->GetRow()));
258 point.SetSigmaY(c->GetSigmaY2()/track->GetCurrentSigmaY2());
259 point.SetSigmaZ(c->GetSigmaZ2()/track->GetCurrentSigmaZ2());
260 point.SetErrY(sqrt(track->GetErrorY2()));
261 point.SetErrZ(sqrt(track->GetErrorZ2()));
263 point.SetX(track->GetX());
264 point.SetY(track->GetY());
265 point.SetZ(track->GetZ());
266 point.SetAngleY(angle2);
267 point.SetAngleZ(track->GetTgl());
268 if (point.IsShared()){
269 track->SetErrorY2(track->GetErrorY2()*4);
270 track->SetErrorZ2(track->GetErrorZ2()*4);
274 Double_t chi2 = track->GetPredictedChi2(track->GetCurrentCluster());
276 // track->SetErrorY2(track->GetErrorY2()*1.3);
277 // track->SetErrorY2(track->GetErrorY2()+0.01);
278 // track->SetErrorZ2(track->GetErrorZ2()*1.3);
279 // track->SetErrorZ2(track->GetErrorZ2()+0.005);
281 if (accept>0) return 0;
282 if (track->GetNumberOfClusters()%20==0){
283 // if (track->fHelixIn){
284 // TClonesArray & larr = *(track->fHelixIn);
285 // Int_t ihelix = larr.GetEntriesFast();
286 // new(larr[ihelix]) AliHelix(*track) ;
289 track->SetNoCluster(0);
290 return track->Update(c,chi2,i);
295 Int_t AliTPCtrackerMI::AcceptCluster(AliTPCseed * seed, AliTPCclusterMI * cluster)
298 // decide according desired precision to accept given
299 // cluster for tracking
301 seed->GetProlongation(cluster->GetX(),yt,zt);
302 Double_t sy2=ErrY2(seed,cluster);
303 Double_t sz2=ErrZ2(seed,cluster);
305 Double_t sdistancey2 = sy2+seed->GetSigmaY2();
306 Double_t sdistancez2 = sz2+seed->GetSigmaZ2();
307 Double_t dy=seed->GetCurrentCluster()->GetY()-yt;
308 Double_t dz=seed->GetCurrentCluster()->GetZ()-zt;
309 Double_t rdistancey2 = (seed->GetCurrentCluster()->GetY()-yt)*
310 (seed->GetCurrentCluster()->GetY()-yt)/sdistancey2;
311 Double_t rdistancez2 = (seed->GetCurrentCluster()->GetZ()-zt)*
312 (seed->GetCurrentCluster()->GetZ()-zt)/sdistancez2;
314 Double_t rdistance2 = rdistancey2+rdistancez2;
317 if (AliTPCReconstructor::StreamLevel()>2 && seed->GetNumberOfClusters()>20) {
318 Float_t rmsy2 = seed->GetCurrentSigmaY2();
319 Float_t rmsz2 = seed->GetCurrentSigmaZ2();
320 Float_t rmsy2p30 = seed->GetCMeanSigmaY2p30();
321 Float_t rmsz2p30 = seed->GetCMeanSigmaZ2p30();
322 Float_t rmsy2p30R = seed->GetCMeanSigmaY2p30R();
323 Float_t rmsz2p30R = seed->GetCMeanSigmaZ2p30R();
324 AliExternalTrackParam param(*seed);
325 static TVectorD gcl(3),gtr(3);
327 param.GetXYZ(gcl.GetMatrixArray());
328 cluster->GetGlobalXYZ(gclf);
329 gcl[0]=gclf[0]; gcl[1]=gclf[1]; gcl[2]=gclf[2];
332 if (AliTPCReconstructor::StreamLevel()>2) {
333 (*fDebugStreamer)<<"ErrParam"<<
346 "rmsy2p30="<<rmsy2p30<<
347 "rmsz2p30="<<rmsz2p30<<
348 "rmsy2p30R="<<rmsy2p30R<<
349 "rmsz2p30R="<<rmsz2p30R<<
350 // normalize distance -
351 "rdisty="<<rdistancey2<<
352 "rdistz="<<rdistancez2<<
353 "rdist="<<rdistance2<< //
357 //return 0; // temporary
358 if (rdistance2>32) return 3;
361 if ((rdistancey2>9. || rdistancez2>9.) && cluster->GetType()==0)
362 return 2; //suspisiouce - will be changed
364 if ((rdistancey2>6.25 || rdistancez2>6.25) && cluster->GetType()>0)
365 // strict cut on overlaped cluster
366 return 2; //suspisiouce - will be changed
368 if ( (rdistancey2>1. || rdistancez2>6.25 )
369 && cluster->GetType()<0){
370 seed->SetNFoundable(seed->GetNFoundable()-1);
380 //_____________________________________________________________________________
381 AliTPCtrackerMI::AliTPCtrackerMI(const AliTPCParam *par):
383 fkNIS(par->GetNInnerSector()/2),
385 fkNOS(par->GetNOuterSector()/2),
402 //---------------------------------------------------------------------
403 // The main TPC tracker constructor
404 //---------------------------------------------------------------------
405 fInnerSec=new AliTPCtrackerSector[fkNIS];
406 fOuterSec=new AliTPCtrackerSector[fkNOS];
409 for (i=0; i<fkNIS; i++) fInnerSec[i].Setup(par,0);
410 for (i=0; i<fkNOS; i++) fOuterSec[i].Setup(par,1);
413 Int_t nrowlow = par->GetNRowLow();
414 Int_t nrowup = par->GetNRowUp();
417 for (i=0;i<nrowlow;i++){
418 fXRow[i] = par->GetPadRowRadiiLow(i);
419 fPadLength[i]= par->GetPadPitchLength(0,i);
420 fYMax[i] = fXRow[i]*TMath::Tan(0.5*par->GetInnerAngle());
424 for (i=0;i<nrowup;i++){
425 fXRow[i+nrowlow] = par->GetPadRowRadiiUp(i);
426 fPadLength[i+nrowlow] = par->GetPadPitchLength(60,i);
427 fYMax[i+nrowlow] = fXRow[i+nrowlow]*TMath::Tan(0.5*par->GetOuterAngle());
430 if (AliTPCReconstructor::StreamLevel()>0) {
431 fDebugStreamer = new TTreeSRedirector("TPCdebug.root");
434 //________________________________________________________________________
435 AliTPCtrackerMI::AliTPCtrackerMI(const AliTPCtrackerMI &t):
456 //------------------------------------
457 // dummy copy constructor
458 //------------------------------------------------------------------
461 AliTPCtrackerMI & AliTPCtrackerMI::operator=(const AliTPCtrackerMI& /*r*/)
463 //------------------------------
465 //--------------------------------------------------------------
468 //_____________________________________________________________________________
469 AliTPCtrackerMI::~AliTPCtrackerMI() {
470 //------------------------------------------------------------------
471 // TPC tracker destructor
472 //------------------------------------------------------------------
479 if (fDebugStreamer) delete fDebugStreamer;
483 void AliTPCtrackerMI::FillESD(const TObjArray* arr)
487 //fill esds using updated tracks
490 // write tracks to the event
491 // store index of the track
492 Int_t nseed=arr->GetEntriesFast();
493 //FindKinks(arr,fEvent);
494 for (Int_t i=0; i<nseed; i++) {
495 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
499 if (AliTPCReconstructor::StreamLevel()>1) {
500 (*fDebugStreamer)<<"Track0"<<
504 // pt->PropagateTo(fkParam->GetInnerRadiusLow());
505 if (pt->GetKinkIndex(0)<=0){ //don't propagate daughter tracks
506 pt->PropagateTo(fkParam->GetInnerRadiusLow());
509 if (( pt->GetPoints()[2]- pt->GetPoints()[0])>5 && pt->GetPoints()[3]>0.8){
511 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
512 iotrack.SetTPCPoints(pt->GetPoints());
513 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
514 iotrack.SetV0Indexes(pt->GetV0Indexes());
515 // iotrack.SetTPCpid(pt->fTPCr);
516 //iotrack.SetTPCindex(i);
517 fEvent->AddTrack(&iotrack);
521 if ( (pt->GetNumberOfClusters()>70)&& (Float_t(pt->GetNumberOfClusters())/Float_t(pt->GetNFoundable()))>0.55) {
523 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
524 iotrack.SetTPCPoints(pt->GetPoints());
525 //iotrack.SetTPCindex(i);
526 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
527 iotrack.SetV0Indexes(pt->GetV0Indexes());
528 // iotrack.SetTPCpid(pt->fTPCr);
529 fEvent->AddTrack(&iotrack);
533 // short tracks - maybe decays
535 if ( (pt->GetNumberOfClusters()>30) && (Float_t(pt->GetNumberOfClusters())/Float_t(pt->GetNFoundable()))>0.70) {
536 Int_t found,foundable,shared;
537 pt->GetClusterStatistic(0,60,found, foundable,shared,kFALSE);
538 if ( (found>20) && (pt->GetNShared()/float(pt->GetNumberOfClusters())<0.2)){
540 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
541 //iotrack.SetTPCindex(i);
542 iotrack.SetTPCPoints(pt->GetPoints());
543 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
544 iotrack.SetV0Indexes(pt->GetV0Indexes());
545 //iotrack.SetTPCpid(pt->fTPCr);
546 fEvent->AddTrack(&iotrack);
551 if ( (pt->GetNumberOfClusters()>20) && (Float_t(pt->GetNumberOfClusters())/Float_t(pt->GetNFoundable()))>0.8) {
552 Int_t found,foundable,shared;
553 pt->GetClusterStatistic(0,60,found, foundable,shared,kFALSE);
554 if (found<20) continue;
555 if (pt->GetNShared()/float(pt->GetNumberOfClusters())>0.2) continue;
558 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
559 iotrack.SetTPCPoints(pt->GetPoints());
560 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
561 iotrack.SetV0Indexes(pt->GetV0Indexes());
562 //iotrack.SetTPCpid(pt->fTPCr);
563 //iotrack.SetTPCindex(i);
564 fEvent->AddTrack(&iotrack);
567 // short tracks - secondaties
569 if ( (pt->GetNumberOfClusters()>30) ) {
570 Int_t found,foundable,shared;
571 pt->GetClusterStatistic(128,158,found, foundable,shared,kFALSE);
572 if ( (found>20) && (pt->GetNShared()/float(pt->GetNumberOfClusters())<0.2) &&float(found)/float(foundable)>0.8){
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);
585 if ( (pt->GetNumberOfClusters()>15)) {
586 Int_t found,foundable,shared;
587 pt->GetClusterStatistic(138,158,found, foundable,shared,kFALSE);
588 if (found<15) continue;
589 if (foundable<=0) continue;
590 if (pt->GetNShared()/float(pt->GetNumberOfClusters())>0.2) continue;
591 if (float(found)/float(foundable)<0.8) continue;
594 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
595 iotrack.SetTPCPoints(pt->GetPoints());
596 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
597 iotrack.SetV0Indexes(pt->GetV0Indexes());
598 // iotrack.SetTPCpid(pt->fTPCr);
599 //iotrack.SetTPCindex(i);
600 fEvent->AddTrack(&iotrack);
604 // >> account for suppressed tracks in the kink indices (RS)
605 int nESDtracks = fEvent->GetNumberOfTracks();
606 for (int it=nESDtracks;it--;) {
607 AliESDtrack* esdTr = fEvent->GetTrack(it);
608 if (!esdTr || !esdTr->GetKinkIndex(0)) continue;
609 for (int ik=0;ik<3;ik++) {
611 if (!(knkId=esdTr->GetKinkIndex(ik))) break; // no more kinks for this track
612 AliESDkink* kink = fEvent->GetKink(TMath::Abs(knkId)-1);
614 AliError(Form("ESDTrack%d refers to non-existing kink %d",it,TMath::Abs(knkId)-1));
617 kink->SetIndex(it, knkId<0 ? 0:1); // update track index of the kink: mother at 0, daughter at 1
620 // << account for suppressed tracks in the kink indices (RS)
621 AliInfo(Form("Number of filled ESDs-\t%d\n",fEvent->GetNumberOfTracks()));
629 Double_t AliTPCtrackerMI::ErrY2(AliTPCseed* seed, const AliTPCclusterMI * cl){
632 // Use calibrated cluster error from OCDB
634 AliTPCClusterParam * clparam = AliTPCcalibDB::Instance()->GetClusterParam();
636 Float_t z = TMath::Abs(fkParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
637 Int_t ctype = cl->GetType();
638 Int_t type = (cl->GetRow()<63) ? 0: (cl->GetRow()>126) ? 1:2;
639 Double_t angle = seed->GetSnp()*seed->GetSnp();
640 angle = TMath::Sqrt(TMath::Abs(angle/(1.-angle)));
641 Double_t erry2 = clparam->GetError0Par(0,type, z,angle);
643 erry2+=0.5; // edge cluster
646 seed->SetErrorY2(erry2);
650 //calculate look-up table at the beginning
651 // static Bool_t ginit = kFALSE;
652 // static Float_t gnoise1,gnoise2,gnoise3;
653 // static Float_t ggg1[10000];
654 // static Float_t ggg2[10000];
655 // static Float_t ggg3[10000];
656 // static Float_t glandau1[10000];
657 // static Float_t glandau2[10000];
658 // static Float_t glandau3[10000];
660 // static Float_t gcor01[500];
661 // static Float_t gcor02[500];
662 // static Float_t gcorp[500];
666 // if (ginit==kFALSE){
667 // for (Int_t i=1;i<500;i++){
668 // Float_t rsigma = float(i)/100.;
669 // gcor02[i] = TMath::Max(0.78 +TMath::Exp(7.4*(rsigma-1.2)),0.6);
670 // gcor01[i] = TMath::Max(0.72 +TMath::Exp(3.36*(rsigma-1.2)),0.6);
671 // gcorp[i] = TMath::Max(TMath::Power((rsigma+0.5),1.5),1.2);
675 // for (Int_t i=3;i<10000;i++){
679 // Float_t amp = float(i);
680 // Float_t padlength =0.75;
681 // gnoise1 = 0.0004/padlength;
682 // Float_t nel = 0.268*amp;
683 // Float_t nprim = 0.155*amp;
684 // ggg1[i] = fkParam->GetDiffT()*fkParam->GetDiffT()*(2+0.001*nel/(padlength*padlength))/nel;
685 // glandau1[i] = (2.+0.12*nprim)*0.5* (2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
686 // if (glandau1[i]>1) glandau1[i]=1;
687 // glandau1[i]*=padlength*padlength/12.;
691 // gnoise2 = 0.0004/padlength;
693 // nprim = 0.133*amp;
694 // ggg2[i] = fkParam->GetDiffT()*fkParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
695 // glandau2[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
696 // if (glandau2[i]>1) glandau2[i]=1;
697 // glandau2[i]*=padlength*padlength/12.;
702 // gnoise3 = 0.0004/padlength;
704 // nprim = 0.133*amp;
705 // ggg3[i] = fkParam->GetDiffT()*fkParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
706 // glandau3[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
707 // if (glandau3[i]>1) glandau3[i]=1;
708 // glandau3[i]*=padlength*padlength/12.;
716 // Int_t amp = int(TMath::Abs(cl->GetQ()));
718 // seed->SetErrorY2(1.);
722 // Float_t z = TMath::Abs(fkParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
723 // Int_t ctype = cl->GetType();
724 // Float_t padlength= GetPadPitchLength(seed->GetRow());
725 // Double_t angle2 = seed->GetSnp()*seed->GetSnp();
726 // angle2 = angle2/(1-angle2);
728 // //cluster "quality"
729 // Int_t rsigmay = int(100.*cl->GetSigmaY2()/(seed->GetCurrentSigmaY2()));
732 // if (fSectors==fInnerSec){
733 // snoise2 = gnoise1;
734 // res = ggg1[amp]*z+glandau1[amp]*angle2;
735 // if (ctype==0) res *= gcor01[rsigmay];
738 // res*= gcorp[rsigmay];
742 // if (padlength<1.1){
743 // snoise2 = gnoise2;
744 // res = ggg2[amp]*z+glandau2[amp]*angle2;
745 // if (ctype==0) res *= gcor02[rsigmay];
748 // res*= gcorp[rsigmay];
752 // snoise2 = gnoise3;
753 // res = ggg3[amp]*z+glandau3[amp]*angle2;
754 // if (ctype==0) res *= gcor02[rsigmay];
757 // res*= gcorp[rsigmay];
764 // res*=2.4; // overestimate error 2 times
768 // if (res<2*snoise2)
771 // seed->SetErrorY2(res);
779 Double_t AliTPCtrackerMI::ErrZ2(AliTPCseed* seed, const AliTPCclusterMI * cl){
782 // Use calibrated cluster error from OCDB
784 AliTPCClusterParam * clparam = AliTPCcalibDB::Instance()->GetClusterParam();
786 Float_t z = TMath::Abs(fkParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
787 Int_t ctype = cl->GetType();
788 Int_t type = (cl->GetRow()<63) ? 0: (cl->GetRow()>126) ? 1:2;
790 Double_t angle2 = seed->GetSnp()*seed->GetSnp();
791 angle2 = seed->GetTgl()*seed->GetTgl()*(1+angle2/(1-angle2));
792 Double_t angle = TMath::Sqrt(TMath::Abs(angle2));
793 Double_t errz2 = clparam->GetError0Par(1,type, z,angle);
795 errz2+=0.5; // edge cluster
798 seed->SetErrorZ2(errz2);
804 // //seed->SetErrorY2(0.1);
806 // //calculate look-up table at the beginning
807 // static Bool_t ginit = kFALSE;
808 // static Float_t gnoise1,gnoise2,gnoise3;
809 // static Float_t ggg1[10000];
810 // static Float_t ggg2[10000];
811 // static Float_t ggg3[10000];
812 // static Float_t glandau1[10000];
813 // static Float_t glandau2[10000];
814 // static Float_t glandau3[10000];
816 // static Float_t gcor01[1000];
817 // static Float_t gcor02[1000];
818 // static Float_t gcorp[1000];
822 // if (ginit==kFALSE){
823 // for (Int_t i=1;i<1000;i++){
824 // Float_t rsigma = float(i)/100.;
825 // gcor02[i] = TMath::Max(0.81 +TMath::Exp(6.8*(rsigma-1.2)),0.6);
826 // gcor01[i] = TMath::Max(0.72 +TMath::Exp(2.04*(rsigma-1.2)),0.6);
827 // gcorp[i] = TMath::Max(TMath::Power((rsigma+0.5),1.5),1.2);
831 // for (Int_t i=3;i<10000;i++){
835 // Float_t amp = float(i);
836 // Float_t padlength =0.75;
837 // gnoise1 = 0.0004/padlength;
838 // Float_t nel = 0.268*amp;
839 // Float_t nprim = 0.155*amp;
840 // ggg1[i] = fkParam->GetDiffT()*fkParam->GetDiffT()*(2+0.001*nel/(padlength*padlength))/nel;
841 // glandau1[i] = (2.+0.12*nprim)*0.5* (2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
842 // if (glandau1[i]>1) glandau1[i]=1;
843 // glandau1[i]*=padlength*padlength/12.;
847 // gnoise2 = 0.0004/padlength;
849 // nprim = 0.133*amp;
850 // ggg2[i] = fkParam->GetDiffT()*fkParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
851 // glandau2[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
852 // if (glandau2[i]>1) glandau2[i]=1;
853 // glandau2[i]*=padlength*padlength/12.;
858 // gnoise3 = 0.0004/padlength;
860 // nprim = 0.133*amp;
861 // ggg3[i] = fkParam->GetDiffT()*fkParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
862 // glandau3[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
863 // if (glandau3[i]>1) glandau3[i]=1;
864 // glandau3[i]*=padlength*padlength/12.;
872 // Int_t amp = int(TMath::Abs(cl->GetQ()));
874 // seed->SetErrorY2(1.);
878 // Float_t z = TMath::Abs(fkParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
879 // Int_t ctype = cl->GetType();
880 // Float_t padlength= GetPadPitchLength(seed->GetRow());
882 // Double_t angle2 = seed->GetSnp()*seed->GetSnp();
883 // // if (angle2<0.6) angle2 = 0.6;
884 // angle2 = seed->GetTgl()*seed->GetTgl()*(1+angle2/(1-angle2));
886 // //cluster "quality"
887 // Int_t rsigmaz = int(100.*cl->GetSigmaZ2()/(seed->GetCurrentSigmaZ2()));
890 // if (fSectors==fInnerSec){
891 // snoise2 = gnoise1;
892 // res = ggg1[amp]*z+glandau1[amp]*angle2;
893 // if (ctype==0) res *= gcor01[rsigmaz];
896 // res*= gcorp[rsigmaz];
900 // if (padlength<1.1){
901 // snoise2 = gnoise2;
902 // res = ggg2[amp]*z+glandau2[amp]*angle2;
903 // if (ctype==0) res *= gcor02[rsigmaz];
906 // res*= gcorp[rsigmaz];
910 // snoise2 = gnoise3;
911 // res = ggg3[amp]*z+glandau3[amp]*angle2;
912 // if (ctype==0) res *= gcor02[rsigmaz];
915 // res*= gcorp[rsigmaz];
924 // if ((ctype<0) &&<70){
929 // if (res<2*snoise2)
931 // if (res>3) res =3;
932 // seed->SetErrorZ2(res);
940 void AliTPCtrackerMI::RotateToLocal(AliTPCseed *seed)
942 //rotate to track "local coordinata
943 Float_t x = seed->GetX();
944 Float_t y = seed->GetY();
945 Float_t ymax = x*TMath::Tan(0.5*fSectors->GetAlpha());
948 seed->SetRelativeSector((seed->GetRelativeSector()+1) % fN);
949 if (!seed->Rotate(fSectors->GetAlpha()))
951 } else if (y <-ymax) {
952 seed->SetRelativeSector((seed->GetRelativeSector()-1+fN) % fN);
953 if (!seed->Rotate(-fSectors->GetAlpha()))
961 //_____________________________________________________________________________
962 Double_t AliTPCtrackerMI::F1old(Double_t x1,Double_t y1,
963 Double_t x2,Double_t y2,
964 Double_t x3,Double_t y3) const
966 //-----------------------------------------------------------------
967 // Initial approximation of the track curvature
968 //-----------------------------------------------------------------
969 Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
970 Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
971 (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
972 Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
973 (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
975 Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
976 if ( xr*xr+yr*yr<=0.00000000000001) return 100;
977 return -xr*yr/sqrt(xr*xr+yr*yr);
982 //_____________________________________________________________________________
983 Double_t AliTPCtrackerMI::F1(Double_t x1,Double_t y1,
984 Double_t x2,Double_t y2,
985 Double_t x3,Double_t y3) const
987 //-----------------------------------------------------------------
988 // Initial approximation of the track curvature
989 //-----------------------------------------------------------------
995 Double_t det = x3*y2-x2*y3;
996 if (TMath::Abs(det)<1e-10){
1000 Double_t u = 0.5* (x2*(x2-x3)+y2*(y2-y3))/det;
1001 Double_t x0 = x3*0.5-y3*u;
1002 Double_t y0 = y3*0.5+x3*u;
1003 Double_t c2 = 1/TMath::Sqrt(x0*x0+y0*y0);
1009 Double_t AliTPCtrackerMI::F2(Double_t x1,Double_t y1,
1010 Double_t x2,Double_t y2,
1011 Double_t x3,Double_t y3) const
1013 //-----------------------------------------------------------------
1014 // Initial approximation of the track curvature
1015 //-----------------------------------------------------------------
1021 Double_t det = x3*y2-x2*y3;
1022 if (TMath::Abs(det)<1e-10) {
1026 Double_t u = 0.5* (x2*(x2-x3)+y2*(y2-y3))/det;
1027 Double_t x0 = x3*0.5-y3*u;
1028 Double_t y0 = y3*0.5+x3*u;
1029 Double_t c2 = 1/TMath::Sqrt(x0*x0+y0*y0);
1038 //_____________________________________________________________________________
1039 Double_t AliTPCtrackerMI::F2old(Double_t x1,Double_t y1,
1040 Double_t x2,Double_t y2,
1041 Double_t x3,Double_t y3) const
1043 //-----------------------------------------------------------------
1044 // Initial approximation of the track curvature times center of curvature
1045 //-----------------------------------------------------------------
1046 Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
1047 Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
1048 (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
1049 Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
1050 (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
1052 Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
1054 return -a/(d*y1-b)*xr/sqrt(xr*xr+yr*yr);
1057 //_____________________________________________________________________________
1058 Double_t AliTPCtrackerMI::F3(Double_t x1,Double_t y1,
1059 Double_t x2,Double_t y2,
1060 Double_t z1,Double_t z2) const
1062 //-----------------------------------------------------------------
1063 // Initial approximation of the tangent of the track dip angle
1064 //-----------------------------------------------------------------
1065 return (z1 - z2)/sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
1069 Double_t AliTPCtrackerMI::F3n(Double_t x1,Double_t y1,
1070 Double_t x2,Double_t y2,
1071 Double_t z1,Double_t z2, Double_t c) const
1073 //-----------------------------------------------------------------
1074 // Initial approximation of the tangent of the track dip angle
1075 //-----------------------------------------------------------------
1079 //angle1 = (z1-z2)*c/(TMath::ASin(c*x1-ni)-TMath::ASin(c*x2-ni));
1081 Double_t d = TMath::Sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
1082 if (TMath::Abs(d*c*0.5)>1) return 0;
1083 // Double_t angle2 = TMath::ASin(d*c*0.5);
1084 // Double_t angle2 = AliTPCFastMath::FastAsin(d*c*0.5);
1085 Double_t angle2 = (d*c*0.5>0.1)? TMath::ASin(d*c*0.5): AliTPCFastMath::FastAsin(d*c*0.5);
1087 angle2 = (z1-z2)*c/(angle2*2.);
1091 Bool_t AliTPCtrackerMI::GetProlongation(Double_t x1, Double_t x2, Double_t x[5], Double_t &y, Double_t &z) const
1092 {//-----------------------------------------------------------------
1093 // This function find proloncation of a track to a reference plane x=x2.
1094 //-----------------------------------------------------------------
1098 if (TMath::Abs(x[4]*x1 - x[2]) >= 0.999) {
1102 Double_t c1=x[4]*x1 - x[2], r1=TMath::Sqrt((1.-c1)*(1.+c1));
1103 Double_t c2=x[4]*x2 - x[2], r2=TMath::Sqrt((1.-c2)*(1.+c2));
1107 Double_t dy = dx*(c1+c2)/(r1+r2);
1110 Double_t delta = x[4]*dx*(c1+c2)/(c1*r2 + c2*r1);
1112 if (TMath::Abs(delta)>0.01){
1113 dz = x[3]*TMath::ASin(delta)/x[4];
1115 dz = x[3]*AliTPCFastMath::FastAsin(delta)/x[4];
1118 //dz = x[3]*AliTPCFastMath::FastAsin(delta)/x[4];
1126 Int_t AliTPCtrackerMI::LoadClusters (TTree *const tree)
1131 return LoadClusters();
1135 Int_t AliTPCtrackerMI::LoadClusters(const TObjArray *arr)
1138 // load clusters to the memory
1139 AliTPCClustersRow *clrow = new AliTPCClustersRow("AliTPCclusterMI");
1140 Int_t lower = arr->LowerBound();
1141 Int_t entries = arr->GetEntriesFast();
1143 for (Int_t i=lower; i<entries; i++) {
1144 clrow = (AliTPCClustersRow*) arr->At(i);
1145 if(!clrow) continue;
1146 if(!clrow->GetArray()) continue;
1150 fkParam->AdjustSectorRow(clrow->GetID(),sec,row);
1152 for (Int_t icl=0; icl<clrow->GetArray()->GetEntriesFast(); icl++){
1153 Transform((AliTPCclusterMI*)(clrow->GetArray()->At(icl)));
1156 if (clrow->GetArray()->GetEntriesFast()<=0) continue;
1157 AliTPCtrackerRow * tpcrow=0;
1160 tpcrow = &(fInnerSec[sec%fkNIS][row]);
1164 tpcrow = &(fOuterSec[(sec-fkNIS*2)%fkNOS][row]);
1165 left = (sec-fkNIS*2)/fkNOS;
1168 tpcrow->SetN1(clrow->GetArray()->GetEntriesFast());
1169 tpcrow->SetClusters1(new AliTPCclusterMI[tpcrow->GetN1()]);
1170 for (Int_t j=0;j<tpcrow->GetN1();++j)
1171 tpcrow->SetCluster1(j, *(AliTPCclusterMI*)(clrow->GetArray()->At(j)));
1174 tpcrow->SetN2(clrow->GetArray()->GetEntriesFast());
1175 tpcrow->SetClusters2(new AliTPCclusterMI[tpcrow->GetN2()]);
1176 for (Int_t j=0;j<tpcrow->GetN2();++j)
1177 tpcrow->SetCluster2(j,*(AliTPCclusterMI*)(clrow->GetArray()->At(j)));
1179 clrow->GetArray()->Clear("C");
1188 Int_t AliTPCtrackerMI::LoadClusters(const TClonesArray *arr)
1191 // load clusters to the memory from one
1194 AliTPCclusterMI *clust=0;
1195 Int_t count[72][96] = { {0} , {0} };
1197 // loop over clusters
1198 for (Int_t icl=0; icl<arr->GetEntriesFast(); icl++) {
1199 clust = (AliTPCclusterMI*)arr->At(icl);
1200 if(!clust) continue;
1201 //printf("cluster: det %d, row %d \n", clust->GetDetector(),clust->GetRow());
1203 // transform clusters
1206 // count clusters per pad row
1207 count[clust->GetDetector()][clust->GetRow()]++;
1210 // insert clusters to sectors
1211 for (Int_t icl=0; icl<arr->GetEntriesFast(); icl++) {
1212 clust = (AliTPCclusterMI*)arr->At(icl);
1213 if(!clust) continue;
1215 Int_t sec = clust->GetDetector();
1216 Int_t row = clust->GetRow();
1218 // filter overlapping pad rows needed by HLT
1219 if(sec<fkNIS*2) { //IROCs
1220 if(row == 30) continue;
1223 if(row == 27 || row == 76) continue;
1229 fInnerSec[sec%fkNIS].InsertCluster(clust, count[sec][row], fkParam);
1232 left = (sec-fkNIS*2)/fkNOS;
1233 fOuterSec[(sec-fkNIS*2)%fkNOS].InsertCluster(clust, count[sec][row], fkParam);
1237 // Load functions must be called behind LoadCluster(TClonesArray*)
1239 //LoadOuterSectors();
1240 //LoadInnerSectors();
1246 Int_t AliTPCtrackerMI::LoadClusters()
1249 // load clusters to the memory
1250 AliTPCClustersRow *clrow= new AliTPCClustersRow("AliTPCclusterMI");
1252 // TTree * tree = fClustersArray.GetTree();
1254 TTree * tree = fInput;
1255 TBranch * br = tree->GetBranch("Segment");
1256 br->SetAddress(&clrow);
1258 Int_t j=Int_t(tree->GetEntries());
1259 for (Int_t i=0; i<j; i++) {
1263 fkParam->AdjustSectorRow(clrow->GetID(),sec,row);
1264 for (Int_t icl=0; icl<clrow->GetArray()->GetEntriesFast(); icl++){
1265 Transform((AliTPCclusterMI*)(clrow->GetArray()->At(icl)));
1268 AliTPCtrackerRow * tpcrow=0;
1271 tpcrow = &(fInnerSec[sec%fkNIS][row]);
1275 tpcrow = &(fOuterSec[(sec-fkNIS*2)%fkNOS][row]);
1276 left = (sec-fkNIS*2)/fkNOS;
1279 tpcrow->SetN1(clrow->GetArray()->GetEntriesFast());
1280 tpcrow->SetClusters1(new AliTPCclusterMI[tpcrow->GetN1()]);
1281 for (Int_t k=0;k<tpcrow->GetN1();++k)
1282 tpcrow->SetCluster1(k, *(AliTPCclusterMI*)(clrow->GetArray()->At(k)));
1285 tpcrow->SetN2(clrow->GetArray()->GetEntriesFast());
1286 tpcrow->SetClusters2(new AliTPCclusterMI[tpcrow->GetN2()]);
1287 for (Int_t k=0;k<tpcrow->GetN2();++k)
1288 tpcrow->SetCluster2(k,*(AliTPCclusterMI*)(clrow->GetArray()->At(k)));
1299 void AliTPCtrackerMI::UnloadClusters()
1302 // unload clusters from the memory
1304 Int_t nrows = fOuterSec->GetNRows();
1305 for (Int_t sec = 0;sec<fkNOS;sec++)
1306 for (Int_t row = 0;row<nrows;row++){
1307 AliTPCtrackerRow* tpcrow = &(fOuterSec[sec%fkNOS][row]);
1309 // if (tpcrow->fClusters1) delete []tpcrow->fClusters1;
1310 // if (tpcrow->fClusters2) delete []tpcrow->fClusters2;
1312 tpcrow->ResetClusters();
1315 nrows = fInnerSec->GetNRows();
1316 for (Int_t sec = 0;sec<fkNIS;sec++)
1317 for (Int_t row = 0;row<nrows;row++){
1318 AliTPCtrackerRow* tpcrow = &(fInnerSec[sec%fkNIS][row]);
1320 // if (tpcrow->fClusters1) delete []tpcrow->fClusters1;
1321 //if (tpcrow->fClusters2) delete []tpcrow->fClusters2;
1323 tpcrow->ResetClusters();
1329 void AliTPCtrackerMI::FillClusterArray(TObjArray* array) const{
1331 // Filling cluster to the array - For visualization purposes
1334 nrows = fOuterSec->GetNRows();
1335 for (Int_t sec = 0;sec<fkNOS;sec++)
1336 for (Int_t row = 0;row<nrows;row++){
1337 AliTPCtrackerRow* tpcrow = &(fOuterSec[sec%fkNOS][row]);
1338 if (!tpcrow) continue;
1339 for (Int_t icl = 0;icl<tpcrow->GetN();icl++){
1340 array->AddLast((TObject*)((*tpcrow)[icl]));
1343 nrows = fInnerSec->GetNRows();
1344 for (Int_t sec = 0;sec<fkNIS;sec++)
1345 for (Int_t row = 0;row<nrows;row++){
1346 AliTPCtrackerRow* tpcrow = &(fInnerSec[sec%fkNIS][row]);
1347 if (!tpcrow) continue;
1348 for (Int_t icl = 0;icl<tpcrow->GetN();icl++){
1349 array->AddLast((TObject*)(*tpcrow)[icl]);
1355 void AliTPCtrackerMI::Transform(AliTPCclusterMI * cluster){
1359 AliTPCcalibDB * calibDB = AliTPCcalibDB::Instance();
1360 AliTPCTransform *transform = calibDB->GetTransform() ;
1362 AliFatal("Tranformations not in calibDB");
1365 transform->SetCurrentRecoParam((AliTPCRecoParam*)AliTPCReconstructor::GetRecoParam());
1366 Double_t x[3]={cluster->GetRow(),cluster->GetPad(),cluster->GetTimeBin()};
1367 Int_t i[1]={cluster->GetDetector()};
1368 transform->Transform(x,i,0,1);
1369 // if (cluster->GetDetector()%36>17){
1374 // in debug mode check the transformation
1376 if (AliTPCReconstructor::StreamLevel()>2) {
1378 cluster->GetGlobalXYZ(gx);
1379 Int_t event = (fEvent==NULL)? 0: fEvent->GetEventNumberInFile();
1380 TTreeSRedirector &cstream = *fDebugStreamer;
1381 cstream<<"Transform"<<
1392 cluster->SetX(x[0]);
1393 cluster->SetY(x[1]);
1394 cluster->SetZ(x[2]);
1399 //if (!fkParam->IsGeoRead()) fkParam->ReadGeoMatrices();
1400 if (AliTPCReconstructor::GetRecoParam()->GetUseSectorAlignment() && (!calibDB->HasAlignmentOCDB())){
1401 TGeoHMatrix *mat = fkParam->GetClusterMatrix(cluster->GetDetector());
1403 Double_t pos[3]= {cluster->GetX(),cluster->GetY(),cluster->GetZ()};
1404 Double_t posC[3]={cluster->GetX(),cluster->GetY(),cluster->GetZ()};
1405 if (mat) mat->LocalToMaster(pos,posC);
1407 // chack Loading of Geo matrices from GeoManager - TEMPORARY FIX
1409 cluster->SetX(posC[0]);
1410 cluster->SetY(posC[1]);
1411 cluster->SetZ(posC[2]);
1415 //_____________________________________________________________________________
1416 Int_t AliTPCtrackerMI::LoadOuterSectors() {
1417 //-----------------------------------------------------------------
1418 // This function fills outer TPC sectors with clusters.
1419 //-----------------------------------------------------------------
1420 Int_t nrows = fOuterSec->GetNRows();
1422 for (Int_t sec = 0;sec<fkNOS;sec++)
1423 for (Int_t row = 0;row<nrows;row++){
1424 AliTPCtrackerRow* tpcrow = &(fOuterSec[sec%fkNOS][row]);
1425 Int_t sec2 = sec+2*fkNIS;
1427 Int_t ncl = tpcrow->GetN1();
1429 AliTPCclusterMI *c= (tpcrow->GetCluster1(ncl));
1430 index=(((sec2<<8)+row)<<16)+ncl;
1431 tpcrow->InsertCluster(c,index);
1434 ncl = tpcrow->GetN2();
1436 AliTPCclusterMI *c= (tpcrow->GetCluster2(ncl));
1437 index=((((sec2+fkNOS)<<8)+row)<<16)+ncl;
1438 tpcrow->InsertCluster(c,index);
1441 // write indexes for fast acces
1443 for (Int_t i=0;i<510;i++)
1444 tpcrow->SetFastCluster(i,-1);
1445 for (Int_t i=0;i<tpcrow->GetN();i++){
1446 Int_t zi = Int_t((*tpcrow)[i]->GetZ()+255.);
1447 tpcrow->SetFastCluster(zi,i); // write index
1450 for (Int_t i=0;i<510;i++){
1451 if (tpcrow->GetFastCluster(i)<0)
1452 tpcrow->SetFastCluster(i,last);
1454 last = tpcrow->GetFastCluster(i);
1463 //_____________________________________________________________________________
1464 Int_t AliTPCtrackerMI::LoadInnerSectors() {
1465 //-----------------------------------------------------------------
1466 // This function fills inner TPC sectors with clusters.
1467 //-----------------------------------------------------------------
1468 Int_t nrows = fInnerSec->GetNRows();
1470 for (Int_t sec = 0;sec<fkNIS;sec++)
1471 for (Int_t row = 0;row<nrows;row++){
1472 AliTPCtrackerRow* tpcrow = &(fInnerSec[sec%fkNIS][row]);
1475 Int_t ncl = tpcrow->GetN1();
1477 AliTPCclusterMI *c= (tpcrow->GetCluster1(ncl));
1478 index=(((sec<<8)+row)<<16)+ncl;
1479 tpcrow->InsertCluster(c,index);
1482 ncl = tpcrow->GetN2();
1484 AliTPCclusterMI *c= (tpcrow->GetCluster2(ncl));
1485 index=((((sec+fkNIS)<<8)+row)<<16)+ncl;
1486 tpcrow->InsertCluster(c,index);
1489 // write indexes for fast acces
1491 for (Int_t i=0;i<510;i++)
1492 tpcrow->SetFastCluster(i,-1);
1493 for (Int_t i=0;i<tpcrow->GetN();i++){
1494 Int_t zi = Int_t((*tpcrow)[i]->GetZ()+255.);
1495 tpcrow->SetFastCluster(zi,i); // write index
1498 for (Int_t i=0;i<510;i++){
1499 if (tpcrow->GetFastCluster(i)<0)
1500 tpcrow->SetFastCluster(i,last);
1502 last = tpcrow->GetFastCluster(i);
1514 //_________________________________________________________________________
1515 AliTPCclusterMI *AliTPCtrackerMI::GetClusterMI(Int_t index) const {
1516 //--------------------------------------------------------------------
1517 // Return pointer to a given cluster
1518 //--------------------------------------------------------------------
1519 if (index<0) return 0; // no cluster
1520 Int_t sec=(index&0xff000000)>>24;
1521 Int_t row=(index&0x00ff0000)>>16;
1522 Int_t ncl=(index&0x00007fff)>>00;
1524 const AliTPCtrackerRow * tpcrow=0;
1525 AliTPCclusterMI * clrow =0;
1527 if (sec<0 || sec>=fkNIS*4) {
1528 AliWarning(Form("Wrong sector %d",sec));
1533 AliTPCtrackerSector& tracksec = fInnerSec[sec%fkNIS];
1534 if (tracksec.GetNRows()<=row) return 0;
1535 tpcrow = &(tracksec[row]);
1536 if (tpcrow==0) return 0;
1539 if (tpcrow->GetN1()<=ncl) return 0;
1540 clrow = tpcrow->GetClusters1();
1543 if (tpcrow->GetN2()<=ncl) return 0;
1544 clrow = tpcrow->GetClusters2();
1548 AliTPCtrackerSector& tracksec = fOuterSec[(sec-fkNIS*2)%fkNOS];
1549 if (tracksec.GetNRows()<=row) return 0;
1550 tpcrow = &(tracksec[row]);
1551 if (tpcrow==0) return 0;
1553 if (sec-2*fkNIS<fkNOS) {
1554 if (tpcrow->GetN1()<=ncl) return 0;
1555 clrow = tpcrow->GetClusters1();
1558 if (tpcrow->GetN2()<=ncl) return 0;
1559 clrow = tpcrow->GetClusters2();
1563 return &(clrow[ncl]);
1569 Int_t AliTPCtrackerMI::FollowToNext(AliTPCseed& t, Int_t nr) {
1570 //-----------------------------------------------------------------
1571 // This function tries to find a track prolongation to next pad row
1572 //-----------------------------------------------------------------
1574 Double_t x= GetXrow(nr), ymax=GetMaxY(nr);
1577 AliTPCclusterMI *cl=0;
1578 Int_t tpcindex= t.GetClusterIndex2(nr);
1580 // update current shape info every 5 pad-row
1581 // if ( (nr%5==0) || t.GetNumberOfClusters()<2 || (t.fCurrentSigmaY2<0.0001) ){
1585 if (fIteration>0 && tpcindex>=-1){ //if we have already clusters
1587 if (tpcindex==-1) return 0; //track in dead zone
1589 cl = t.GetClusterPointer(nr);
1590 if ( (cl==0) ) cl = GetClusterMI(tpcindex);
1591 t.SetCurrentClusterIndex1(tpcindex);
1594 Int_t relativesector = ((tpcindex&0xff000000)>>24)%18; // if previously accepted cluster in different sector
1595 Float_t angle = relativesector*fSectors->GetAlpha()+fSectors->GetAlphaShift();
1597 if (angle<-TMath::Pi()) angle += 2*TMath::Pi();
1598 if (angle>=TMath::Pi()) angle -= 2*TMath::Pi();
1600 if (TMath::Abs(angle-t.GetAlpha())>0.001){
1601 Double_t rotation = angle-t.GetAlpha();
1602 t.SetRelativeSector(relativesector);
1603 if (!t.Rotate(rotation)) return 0;
1605 if (!t.PropagateTo(x)) return 0;
1607 t.SetCurrentCluster(cl);
1609 Int_t accept = AcceptCluster(&t,t.GetCurrentCluster());
1610 if ((tpcindex&0x8000)==0) accept =0;
1612 //if founded cluster is acceptible
1613 if (cl->IsUsed(11)) { // id cluster is shared inrease uncertainty
1614 t.SetErrorY2(t.GetErrorY2()+0.03);
1615 t.SetErrorZ2(t.GetErrorZ2()+0.03);
1616 t.SetErrorY2(t.GetErrorY2()*3);
1617 t.SetErrorZ2(t.GetErrorZ2()*3);
1619 t.SetNFoundable(t.GetNFoundable()+1);
1620 UpdateTrack(&t,accept);
1625 if (TMath::Abs(t.GetSnp())>AliTPCReconstructor::GetMaxSnpTracker()) return 0; // cut on angle
1626 if (fIteration>1 && IsFindable(t)){
1627 // not look for new cluster during refitting
1628 t.SetNFoundable(t.GetNFoundable()+1);
1633 // if (TMath::Abs(t.GetSnp())>0.95 || TMath::Abs(x*t.GetC()-t.GetEta())>0.95) return 0;// patch 28 fev 06
1634 Double_t y=t.GetYat(x);
1635 if (TMath::Abs(y)>ymax){
1637 t.SetRelativeSector((t.GetRelativeSector()+1) % fN);
1638 if (!t.Rotate(fSectors->GetAlpha()))
1640 } else if (y <-ymax) {
1641 t.SetRelativeSector((t.GetRelativeSector()-1+fN) % fN);
1642 if (!t.Rotate(-fSectors->GetAlpha()))
1648 if (!t.PropagateTo(x)) {
1649 if (fIteration==0) t.SetRemoval(10);
1653 Double_t z=t.GetZ();
1656 if (!IsActive(t.GetRelativeSector(),nr)) {
1658 t.SetClusterIndex2(nr,-1);
1661 //AliInfo(Form("A - Sector%d phi %f - alpha %f", t.fRelativeSector,y/x, t.GetAlpha()));
1662 Bool_t isActive = IsActive(t.GetRelativeSector(),nr);
1663 Bool_t isActive2 = (nr>=fInnerSec->GetNRows()) ? fOuterSec[t.GetRelativeSector()][nr-fInnerSec->GetNRows()].GetN()>0:fInnerSec[t.GetRelativeSector()][nr].GetN()>0;
1665 if (!isActive || !isActive2) return 0;
1667 const AliTPCtrackerRow &krow=GetRow(t.GetRelativeSector(),nr);
1668 if ( (t.GetSigmaY2()<0) || t.GetSigmaZ2()<0) return 0;
1670 Double_t roadz = 1.;
1672 if (TMath::Abs(TMath::Abs(y)-ymax)<krow.GetDeadZone()){
1674 t.SetClusterIndex2(nr,-1);
1680 // if (TMath::Abs(z)<(AliTPCReconstructor::GetCtgRange()*x+10) && TMath::Abs(z)<fkParam->GetZLength(0) && (TMath::Abs(t.GetSnp())<AliTPCReconstructor::GetMaxSnpTracker()))
1681 t.SetNFoundable(t.GetNFoundable()+1);
1687 // cl = krow.FindNearest2(y+10.,z,roady,roadz,index);
1688 cl = krow.FindNearest2(y,z,roady,roadz,index);
1689 if (cl) t.SetCurrentClusterIndex1(krow.GetIndex(index));
1692 t.SetCurrentCluster(cl);
1694 if (fIteration==2&&cl->IsUsed(10)) return 0;
1695 Int_t accept = AcceptCluster(&t,t.GetCurrentCluster());
1696 if (fIteration==2&&cl->IsUsed(11)) {
1697 t.SetErrorY2(t.GetErrorY2()+0.03);
1698 t.SetErrorZ2(t.GetErrorZ2()+0.03);
1699 t.SetErrorY2(t.GetErrorY2()*3);
1700 t.SetErrorZ2(t.GetErrorZ2()*3);
1703 if (t.fCurrentCluster->IsUsed(10)){
1708 if (t.fNShared>0.7*t.GetNumberOfClusters()) {
1714 if (accept<3) UpdateTrack(&t,accept);
1717 if ( fIteration==0 && t.GetNFoundable()*0.5 > t.GetNumberOfClusters()) t.SetRemoval(10);
1725 //_________________________________________________________________________
1726 Bool_t AliTPCtrackerMI::GetTrackPoint(Int_t index, AliTrackPoint &p ) const
1728 // Get track space point by index
1729 // return false in case the cluster doesn't exist
1730 AliTPCclusterMI *cl = GetClusterMI(index);
1731 if (!cl) return kFALSE;
1732 Int_t sector = (index&0xff000000)>>24;
1733 // Int_t row = (index&0x00ff0000)>>16;
1735 // xyz[0] = fkParam->GetPadRowRadii(sector,row);
1736 xyz[0] = cl->GetX();
1737 xyz[1] = cl->GetY();
1738 xyz[2] = cl->GetZ();
1740 fkParam->AdjustCosSin(sector,cos,sin);
1741 Float_t x = cos*xyz[0]-sin*xyz[1];
1742 Float_t y = cos*xyz[1]+sin*xyz[0];
1744 Float_t sigmaY2 = 0.027*cl->GetSigmaY2();
1745 if (sector < fkParam->GetNInnerSector()) sigmaY2 *= 2.07;
1746 Float_t sigmaZ2 = 0.066*cl->GetSigmaZ2();
1747 if (sector < fkParam->GetNInnerSector()) sigmaZ2 *= 1.77;
1748 cov[0] = sin*sin*sigmaY2;
1749 cov[1] = -sin*cos*sigmaY2;
1751 cov[3] = cos*cos*sigmaY2;
1754 p.SetXYZ(x,y,xyz[2],cov);
1755 AliGeomManager::ELayerID iLayer;
1757 if (sector < fkParam->GetNInnerSector()) {
1758 iLayer = AliGeomManager::kTPC1;
1762 iLayer = AliGeomManager::kTPC2;
1763 idet = sector - fkParam->GetNInnerSector();
1765 UShort_t volid = AliGeomManager::LayerToVolUID(iLayer,idet);
1766 p.SetVolumeID(volid);
1772 Int_t AliTPCtrackerMI::UpdateClusters(AliTPCseed& t, Int_t nr) {
1773 //-----------------------------------------------------------------
1774 // This function tries to find a track prolongation to next pad row
1775 //-----------------------------------------------------------------
1776 t.SetCurrentCluster(0);
1777 t.SetCurrentClusterIndex1(0);
1779 Double_t xt=t.GetX();
1780 Int_t row = GetRowNumber(xt)-1;
1781 Double_t ymax= GetMaxY(nr);
1783 if (row < nr) return 1; // don't prolongate if not information until now -
1784 // if (TMath::Abs(t.GetSnp())>0.9 && t.GetNumberOfClusters()>40. && fIteration!=2) {
1786 // return 0; // not prolongate strongly inclined tracks
1788 // if (TMath::Abs(t.GetSnp())>0.95) {
1790 // return 0; // not prolongate strongly inclined tracks
1791 // }// patch 28 fev 06
1793 Double_t x= GetXrow(nr);
1795 //t.PropagateTo(x+0.02);
1796 //t.PropagateTo(x+0.01);
1797 if (!t.PropagateTo(x)){
1804 if (TMath::Abs(y)>ymax){
1806 t.SetRelativeSector((t.GetRelativeSector()+1) % fN);
1807 if (!t.Rotate(fSectors->GetAlpha()))
1809 } else if (y <-ymax) {
1810 t.SetRelativeSector((t.GetRelativeSector()-1+fN) % fN);
1811 if (!t.Rotate(-fSectors->GetAlpha()))
1814 // if (!t.PropagateTo(x)){
1821 if (TMath::Abs(t.GetSnp())>AliTPCReconstructor::GetMaxSnpTracker()) return 0;
1823 if (!IsActive(t.GetRelativeSector(),nr)) {
1825 t.SetClusterIndex2(nr,-1);
1828 //AliInfo(Form("A - Sector%d phi %f - alpha %f", t.fRelativeSector,y/x, t.GetAlpha()));
1830 AliTPCtrackerRow &krow=GetRow(t.GetRelativeSector(),nr);
1832 if (TMath::Abs(TMath::Abs(y)-ymax)<krow.GetDeadZone()){
1834 t.SetClusterIndex2(nr,-1);
1840 // if (TMath::Abs(t.GetZ())<(AliTPCReconstructor::GetCtgRange()*t.GetX()+10) && (TMath::Abs(t.GetSnp())<AliTPCReconstructor::GetMaxSnpTracker()))
1841 if (IsFindable(t)) t.SetNFoundable(t.GetNFoundable()+1);
1847 if ( (nr%2==0) || t.GetNumberOfClusters()<2){
1848 // t.fCurrentSigmaY = GetSigmaY(&t);
1849 //t.fCurrentSigmaZ = GetSigmaZ(&t);
1853 AliTPCclusterMI *cl=0;
1856 Double_t roady = 1.;
1857 Double_t roadz = 1.;
1861 index = t.GetClusterIndex2(nr);
1862 if ( (index>0) && (index&0x8000)==0){
1863 cl = t.GetClusterPointer(nr);
1864 if ( (cl==0) && (index>0)) cl = GetClusterMI(index);
1865 t.SetCurrentClusterIndex1(index);
1867 t.SetCurrentCluster(cl);
1873 // if (index<0) return 0;
1874 UInt_t uindex = TMath::Abs(index);
1877 //cl = krow.FindNearest2(y+10,z,roady,roadz,uindex);
1878 cl = krow.FindNearest2(y,z,roady,roadz,uindex);
1881 if (cl) t.SetCurrentClusterIndex1(krow.GetIndex(uindex));
1882 t.SetCurrentCluster(cl);
1888 Int_t AliTPCtrackerMI::FollowToNextCluster(AliTPCseed & t, Int_t nr) {
1889 //-----------------------------------------------------------------
1890 // This function tries to find a track prolongation to next pad row
1891 //-----------------------------------------------------------------
1893 //update error according neighborhoud
1895 if (t.GetCurrentCluster()) {
1897 Int_t accept = AcceptCluster(&t,t.GetCurrentCluster());
1899 if (t.GetCurrentCluster()->IsUsed(10)){
1904 t.SetNShared(t.GetNShared()+1);
1905 if (t.GetNShared()>0.7*t.GetNumberOfClusters()) {
1910 if (fIteration>0) accept = 0;
1911 if (accept<3) UpdateTrack(&t,accept);
1915 if ( ( (t.GetSigmaY2()+t.GetSigmaZ2())>0.16)&& t.GetNumberOfClusters()>18) t.SetRemoval(10);
1916 if ( t.GetChi2()/t.GetNumberOfClusters()>6 &&t.GetNumberOfClusters()>18) t.SetRemoval(10);
1918 if (( (t.GetNFoundable()*0.5 > t.GetNumberOfClusters()) || t.GetNoCluster()>15)) t.SetRemoval(10);
1926 //_____________________________________________________________________________
1927 Int_t AliTPCtrackerMI::FollowProlongation(AliTPCseed& t, Int_t rf, Int_t step) {
1928 //-----------------------------------------------------------------
1929 // This function tries to find a track prolongation.
1930 //-----------------------------------------------------------------
1931 Double_t xt=t.GetX();
1933 Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
1934 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
1935 if (alpha < 0. ) alpha += 2.*TMath::Pi();
1937 t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN);
1939 Int_t first = GetRowNumber(xt)-1;
1940 for (Int_t nr= first; nr>=rf; nr-=step) {
1942 if (t.GetKinkIndexes()[0]>0){
1943 for (Int_t i=0;i<3;i++){
1944 Int_t index = t.GetKinkIndexes()[i];
1945 if (index==0) break;
1946 if (index<0) continue;
1948 AliKink * kink = (AliKink*)fEvent->GetKink(index-1);
1950 printf("PROBLEM\n");
1953 Int_t kinkrow = kink->GetTPCRow0()+2+Int_t(0.5/(0.05+kink->GetAngle(2)));
1955 AliExternalTrackParam paramd(t);
1956 kink->SetDaughter(paramd);
1957 kink->SetStatus(2,5);
1964 if (nr==80) t.UpdateReference();
1965 if (nr<fInnerSec->GetNRows())
1966 fSectors = fInnerSec;
1968 fSectors = fOuterSec;
1969 if (FollowToNext(t,nr)==0)
1982 Int_t AliTPCtrackerMI::FollowBackProlongation(AliTPCseed& t, Int_t rf) {
1983 //-----------------------------------------------------------------
1984 // This function tries to find a track prolongation.
1985 //-----------------------------------------------------------------
1987 Double_t xt=t.GetX();
1988 Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
1989 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
1990 if (alpha < 0. ) alpha += 2.*TMath::Pi();
1991 t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN);
1993 Int_t first = t.GetFirstPoint();
1994 if (first<GetRowNumber(xt)+1) first = GetRowNumber(xt)+1;
1996 if (first<0) first=0;
1997 for (Int_t nr=first; nr<=rf; nr++) {
1998 // if ( (TMath::Abs(t.GetSnp())>0.95)) break;//patch 28 fev 06
1999 if (t.GetKinkIndexes()[0]<0){
2000 for (Int_t i=0;i<3;i++){
2001 Int_t index = t.GetKinkIndexes()[i];
2002 if (index==0) break;
2003 if (index>0) continue;
2004 index = TMath::Abs(index);
2005 AliKink * kink = (AliKink*)fEvent->GetKink(index-1);
2007 printf("PROBLEM\n");
2010 Int_t kinkrow = kink->GetTPCRow0()-2-Int_t(0.5/(0.05+kink->GetAngle(2)));
2012 AliExternalTrackParam paramm(t);
2013 kink->SetMother(paramm);
2014 kink->SetStatus(2,1);
2021 if (nr<fInnerSec->GetNRows())
2022 fSectors = fInnerSec;
2024 fSectors = fOuterSec;
2035 Float_t AliTPCtrackerMI::OverlapFactor(AliTPCseed * s1, AliTPCseed * s2, Int_t &sum1, Int_t & sum2)
2043 Float_t dz2 =(s1->GetZ() - s2->GetZ());
2046 Float_t dy2 =TMath::Abs((s1->GetY() - s2->GetY()));
2048 Float_t distance = TMath::Sqrt(dz2+dy2);
2049 if (distance>4.) return 0; // if there are far away - not overlap - to reduce combinatorics
2052 Int_t firstpoint = TMath::Min(s1->GetFirstPoint(),s2->GetFirstPoint());
2053 Int_t lastpoint = TMath::Max(s1->GetLastPoint(),s2->GetLastPoint());
2058 if (firstpoint>lastpoint) {
2059 firstpoint =lastpoint;
2064 for (Int_t i=firstpoint-1;i<lastpoint+1;i++){
2065 if (s1->GetClusterIndex2(i)>0) sum1++;
2066 if (s2->GetClusterIndex2(i)>0) sum2++;
2067 if (s1->GetClusterIndex2(i)==s2->GetClusterIndex2(i) && s1->GetClusterIndex2(i)>0) {
2071 if (sum<5) return 0;
2073 Float_t summin = TMath::Min(sum1+1,sum2+1);
2074 Float_t ratio = (sum+1)/Float_t(summin);
2078 void AliTPCtrackerMI::SignShared(AliTPCseed * s1, AliTPCseed * s2)
2082 Float_t thetaCut = 0.2;//+10.*TMath::Sqrt(s1->GetSigmaTglZ()+ s2->GetSigmaTglZ());
2083 if (TMath::Abs(s1->GetTgl()-s2->GetTgl())>thetaCut) return;
2084 Float_t minCl = TMath::Min(s1->GetNumberOfClusters(),s2->GetNumberOfClusters());
2085 Int_t cutN0 = TMath::Max(5,TMath::Nint(0.1*minCl));
2090 //Int_t firstpoint = TMath::Max(s1->GetFirstPoint(),s2->GetFirstPoint());
2091 //Int_t lastpoint = TMath::Min(s1->GetLastPoint(),s2->GetLastPoint());
2092 Int_t firstpoint = 0;
2093 Int_t lastpoint = 160;
2095 // if (firstpoint>=lastpoint-5) return;;
2097 for (Int_t i=firstpoint;i<lastpoint;i++){
2098 // if ( (s1->GetClusterIndex2(i)&0xFFFF8FFF)==(s2->GetClusterIndex2(i)&0xFFFF8FFF) && s1->GetClusterIndex2(i)>0) {
2099 if ( (s1->GetClusterIndex2(i))==(s2->GetClusterIndex2(i)) && s1->GetClusterIndex2(i)>0) {
2101 s1->SetSharedMapBit(i, kTRUE);
2102 s2->SetSharedMapBit(i, kTRUE);
2104 if (s1->GetClusterIndex2(i)>0)
2105 s1->SetClusterMapBit(i, kTRUE);
2107 if (sumshared>cutN0){
2110 for (Int_t i=firstpoint;i<lastpoint;i++){
2111 // if ( (s1->GetClusterIndex2(i)&0xFFFF8FFF)==(s2->GetClusterIndex2(i)&0xFFFF8FFF) && s1->GetClusterIndex2(i)>0) {
2112 if ( (s1->GetClusterIndex2(i))==(s2->GetClusterIndex2(i)) && s1->GetClusterIndex2(i)>0) {
2113 AliTPCTrackerPoint *p1 = s1->GetTrackPoint(i);
2114 AliTPCTrackerPoint *p2 = s2->GetTrackPoint(i);;
2115 if (s1->IsActive()&&s2->IsActive()){
2116 p1->SetShared(kTRUE);
2117 p2->SetShared(kTRUE);
2123 if (sumshared>cutN0){
2124 for (Int_t i=0;i<4;i++){
2125 if (s1->GetOverlapLabel(3*i)==0){
2126 s1->SetOverlapLabel(3*i, s2->GetLabel());
2127 s1->SetOverlapLabel(3*i+1,sumshared);
2128 s1->SetOverlapLabel(3*i+2,s2->GetUniqueID());
2132 for (Int_t i=0;i<4;i++){
2133 if (s2->GetOverlapLabel(3*i)==0){
2134 s2->SetOverlapLabel(3*i, s1->GetLabel());
2135 s2->SetOverlapLabel(3*i+1,sumshared);
2136 s2->SetOverlapLabel(3*i+2,s1->GetUniqueID());
2143 void AliTPCtrackerMI::SignShared(TObjArray * arr)
2146 //sort trackss according sectors
2148 for (Int_t i=0; i<arr->GetEntriesFast(); i++) {
2149 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2151 //if (pt) RotateToLocal(pt);
2155 arr->Sort(); // sorting according relative sectors
2156 arr->Expand(arr->GetEntries());
2159 Int_t nseed=arr->GetEntriesFast();
2160 for (Int_t i=0; i<nseed; i++) {
2161 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2163 for (Int_t j=0;j<12;j++){
2164 pt->SetOverlapLabel(j,0);
2167 for (Int_t i=0; i<nseed; i++) {
2168 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2170 if (pt->GetRemoval()>10) continue;
2171 for (Int_t j=i+1; j<nseed; j++){
2172 AliTPCseed *pt2=(AliTPCseed*)arr->UncheckedAt(j);
2173 if (TMath::Abs(pt->GetRelativeSector()-pt2->GetRelativeSector())>1) continue;
2175 if (pt2->GetRemoval()<=10) {
2176 //if ( TMath::Abs(pt->GetRelativeSector()-pt2->GetRelativeSector())>0) break;
2184 void AliTPCtrackerMI::SortTracks(TObjArray * arr, Int_t mode) const
2187 //sort tracks in array according mode criteria
2188 Int_t nseed = arr->GetEntriesFast();
2189 for (Int_t i=0; i<nseed; i++) {
2190 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2201 void AliTPCtrackerMI::RemoveUsed2(TObjArray * arr, Float_t factor1, Float_t factor2, Int_t minimal)
2204 // Loop over all tracks and remove overlaped tracks (with lower quality)
2206 // 1. Unsign clusters
2207 // 2. Sort tracks according quality
2208 // Quality is defined by the number of cluster between first and last points
2210 // 3. Loop over tracks - decreasing quality order
2211 // a.) remove - If the fraction of shared cluster less than factor (1- n or 2)
2212 // b.) remove - If the minimal number of clusters less than minimal and not ITS
2213 // c.) if track accepted - sign clusters
2215 //Called in - AliTPCtrackerMI::Clusters2Tracks()
2216 // - AliTPCtrackerMI::PropagateBack()
2217 // - AliTPCtrackerMI::RefitInward()
2220 // factor1 - factor for constrained
2221 // factor2 - for non constrained tracks
2222 // if (Float_t(shared+1)/Float_t(found+1)>factor) - DELETE
2226 Int_t nseed = arr->GetEntriesFast();
2227 Float_t * quality = new Float_t[nseed];
2228 Int_t * indexes = new Int_t[nseed];
2232 for (Int_t i=0; i<nseed; i++) {
2233 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2238 pt->UpdatePoints(); //select first last max dens points
2239 Float_t * points = pt->GetPoints();
2240 if (points[3]<0.8) quality[i] =-1;
2241 quality[i] = (points[2]-points[0])+pt->GetNumberOfClusters();
2242 //prefer high momenta tracks if overlaps
2243 quality[i] *= TMath::Sqrt(TMath::Abs(pt->Pt())+0.5);
2245 TMath::Sort(nseed,quality,indexes);
2248 for (Int_t itrack=0; itrack<nseed; itrack++) {
2249 Int_t trackindex = indexes[itrack];
2250 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(trackindex);
2253 if (quality[trackindex]<0){
2254 delete arr->RemoveAt(trackindex);
2259 Int_t first = Int_t(pt->GetPoints()[0]);
2260 Int_t last = Int_t(pt->GetPoints()[2]);
2261 Double_t factor = (pt->GetBConstrain()) ? factor1: factor2;
2263 Int_t found,foundable,shared;
2264 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
2265 // pt->GetClusterStatistic(0,160, found, foundable,shared,kFALSE);
2266 Bool_t itsgold =kFALSE;
2269 if (pt->GetESD()->GetITSclusters(dummy)>4) itsgold= kTRUE;
2273 if (Float_t(shared+1)/Float_t(found+1)>factor){
2274 if (pt->GetKinkIndexes()[0]!=0) continue; //don't remove tracks - part of the kinks
2275 if( AliTPCReconstructor::StreamLevel()>3){
2276 TTreeSRedirector &cstream = *fDebugStreamer;
2277 cstream<<"RemoveUsed"<<
2278 "iter="<<fIteration<<
2282 delete arr->RemoveAt(trackindex);
2285 if (pt->GetNumberOfClusters()<50&&(found-0.5*shared)<minimal){ //remove short tracks
2286 if (pt->GetKinkIndexes()[0]!=0) continue; //don't remove tracks - part of the kinks
2287 if( AliTPCReconstructor::StreamLevel()>3){
2288 TTreeSRedirector &cstream = *fDebugStreamer;
2289 cstream<<"RemoveShort"<<
2290 "iter="<<fIteration<<
2294 delete arr->RemoveAt(trackindex);
2300 //if (sharedfactor>0.4) continue;
2301 if (pt->GetKinkIndexes()[0]>0) continue;
2302 //Remove tracks with undefined properties - seems
2303 if (pt->GetSigmaY2()<kAlmost0) continue; // ? what is the origin ?
2305 for (Int_t i=first; i<last; i++) {
2306 Int_t index=pt->GetClusterIndex2(i);
2307 // if (index<0 || index&0x8000 ) continue;
2308 if (index<0 || index&0x8000 ) continue;
2309 AliTPCclusterMI *c= pt->GetClusterPointer(i);
2316 Info("RemoveUsed","\n*****\nNumber of good tracks after shared removal\t%d\n",fNtracks);
2322 void AliTPCtrackerMI::DumpClusters(Int_t iter, TObjArray *trackArray)
2325 // Dump clusters after reco
2326 // signed and unsigned cluster can be visualized
2327 // 1. Unsign all cluster
2328 // 2. Sign all used clusters
2331 Int_t nseed = trackArray->GetEntries();
2332 for (Int_t i=0; i<nseed; i++){
2333 AliTPCseed *pt=(AliTPCseed*)trackArray->UncheckedAt(i);
2337 Bool_t isKink=pt->GetKinkIndex(0)!=0;
2338 for (Int_t j=0; j<160; ++j) {
2339 Int_t index=pt->GetClusterIndex2(j);
2340 if (index<0) continue;
2341 AliTPCclusterMI *c= pt->GetClusterPointer(j);
2343 if (isKink) c->Use(100); // kink
2344 c->Use(10); // by default usage 10
2349 for (Int_t sec=0;sec<fkNIS;sec++){
2350 for (Int_t row=0;row<fInnerSec->GetNRows();row++){
2351 AliTPCclusterMI *cl = fInnerSec[sec][row].GetClusters1();
2352 for (Int_t icl =0;icl< fInnerSec[sec][row].GetN1();icl++){
2353 Float_t gx[3]; cl[icl].GetGlobalXYZ(gx);
2354 (*fDebugStreamer)<<"clDump"<<
2362 cl = fInnerSec[sec][row].GetClusters2();
2363 for (Int_t icl =0;icl< fInnerSec[sec][row].GetN2();icl++){
2364 Float_t gx[3]; cl[icl].GetGlobalXYZ(gx);
2365 (*fDebugStreamer)<<"clDump"<<
2376 for (Int_t sec=0;sec<fkNOS;sec++){
2377 for (Int_t row=0;row<fOuterSec->GetNRows();row++){
2378 AliTPCclusterMI *cl = fOuterSec[sec][row].GetClusters1();
2379 for (Int_t icl =0;icl< fOuterSec[sec][row].GetN1();icl++){
2380 Float_t gx[3]; cl[icl].GetGlobalXYZ(gx);
2381 (*fDebugStreamer)<<"clDump"<<
2389 cl = fOuterSec[sec][row].GetClusters2();
2390 for (Int_t icl =0;icl< fOuterSec[sec][row].GetN2();icl++){
2391 Float_t gx[3]; cl[icl].GetGlobalXYZ(gx);
2392 (*fDebugStreamer)<<"clDump"<<
2404 void AliTPCtrackerMI::UnsignClusters()
2407 // loop over all clusters and unsign them
2410 for (Int_t sec=0;sec<fkNIS;sec++){
2411 for (Int_t row=0;row<fInnerSec->GetNRows();row++){
2412 AliTPCclusterMI *cl = fInnerSec[sec][row].GetClusters1();
2413 for (Int_t icl =0;icl< fInnerSec[sec][row].GetN1();icl++)
2414 // if (cl[icl].IsUsed(10))
2416 cl = fInnerSec[sec][row].GetClusters2();
2417 for (Int_t icl =0;icl< fInnerSec[sec][row].GetN2();icl++)
2418 //if (cl[icl].IsUsed(10))
2423 for (Int_t sec=0;sec<fkNOS;sec++){
2424 for (Int_t row=0;row<fOuterSec->GetNRows();row++){
2425 AliTPCclusterMI *cl = fOuterSec[sec][row].GetClusters1();
2426 for (Int_t icl =0;icl< fOuterSec[sec][row].GetN1();icl++)
2427 //if (cl[icl].IsUsed(10))
2429 cl = fOuterSec[sec][row].GetClusters2();
2430 for (Int_t icl =0;icl< fOuterSec[sec][row].GetN2();icl++)
2431 //if (cl[icl].IsUsed(10))
2440 void AliTPCtrackerMI::SignClusters(const TObjArray * arr, Float_t fnumber, Float_t fdensity)
2443 //sign clusters to be "used"
2445 // snumber and sdensity sign number of sigmas - bellow mean value to be accepted
2446 // loop over "primaries"
2460 Int_t nseed = arr->GetEntriesFast();
2461 for (Int_t i=0; i<nseed; i++) {
2462 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2466 if (!(pt->IsActive())) continue;
2467 Float_t dens = pt->GetNumberOfClusters()/Float_t(pt->GetNFoundable());
2468 if ( (dens>0.7) && (pt->GetNumberOfClusters()>70)){
2470 sumdens2+= dens*dens;
2471 sumn += pt->GetNumberOfClusters();
2472 sumn2 += pt->GetNumberOfClusters()*pt->GetNumberOfClusters();
2473 Float_t chi2 = pt->GetChi2()/pt->GetNumberOfClusters();
2476 sumchi2 +=chi2*chi2;
2481 Float_t mdensity = 0.9;
2482 Float_t meann = 130;
2483 Float_t meanchi = 1;
2484 Float_t sdensity = 0.1;
2485 Float_t smeann = 10;
2486 Float_t smeanchi =0.4;
2490 mdensity = sumdens/sum;
2492 meanchi = sumchi/sum;
2494 sdensity = sumdens2/sum-mdensity*mdensity;
2496 sdensity = TMath::Sqrt(sdensity);
2500 smeann = sumn2/sum-meann*meann;
2502 smeann = TMath::Sqrt(smeann);
2506 smeanchi = sumchi2/sum - meanchi*meanchi;
2508 smeanchi = TMath::Sqrt(smeanchi);
2514 //REMOVE SHORT DELTAS or tracks going out of sensitive volume of TPC
2516 for (Int_t i=0; i<nseed; i++) {
2517 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2521 if (pt->GetBSigned()) continue;
2522 if (pt->GetBConstrain()) continue;
2523 //if (!(pt->IsActive())) continue;
2525 Int_t found,foundable,shared;
2526 pt->GetClusterStatistic(0,160,found, foundable,shared);
2527 if (shared/float(found)>0.3) {
2528 if (shared/float(found)>0.9 ){
2529 //delete arr->RemoveAt(i);
2534 Bool_t isok =kFALSE;
2535 if ( (pt->GetNShared()/pt->GetNumberOfClusters()<0.5) &&pt->GetNumberOfClusters()>60)
2537 if ((TMath::Abs(1/pt->GetC())<100.) && (pt->GetNShared()/pt->GetNumberOfClusters()<0.7))
2539 if (TMath::Abs(pt->GetZ()/pt->GetX())>1.1)
2541 if ( (TMath::Abs(pt->GetSnp()>0.7) && pt->GetD(0,0)>60.))
2545 for (Int_t j=0; j<160; ++j) {
2546 Int_t index=pt->GetClusterIndex2(j);
2547 if (index<0) continue;
2548 AliTPCclusterMI *c= pt->GetClusterPointer(j);
2550 //if (!(c->IsUsed(10))) c->Use();
2557 Double_t maxchi = meanchi+2.*smeanchi;
2559 for (Int_t i=0; i<nseed; i++) {
2560 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2564 //if (!(pt->IsActive())) continue;
2565 if (pt->GetBSigned()) continue;
2566 Double_t chi = pt->GetChi2()/pt->GetNumberOfClusters();
2567 if (chi>maxchi) continue;
2570 Float_t dens = pt->GetNumberOfClusters()/Float_t(pt->GetNFoundable());
2572 //sign only tracks with enoug big density at the beginning
2574 if ((pt->GetDensityFirst(40)<0.75) && pt->GetNumberOfClusters()<meann) continue;
2577 Double_t mindens = TMath::Max(double(mdensity-sdensity*fdensity*bfactor),0.65);
2578 Double_t minn = TMath::Max(Int_t(meann-fnumber*smeann*bfactor),50);
2580 // if (pt->fBConstrain) mindens = TMath::Max(mdensity-sdensity*fdensity*bfactor,0.65);
2581 if ( (pt->GetRemoval()==10) && (pt->GetSnp()>0.8)&&(dens>mindens))
2584 if ((dens>mindens && pt->GetNumberOfClusters()>minn) && chi<maxchi ){
2585 //Int_t noc=pt->GetNumberOfClusters();
2586 pt->SetBSigned(kTRUE);
2587 for (Int_t j=0; j<160; ++j) {
2589 Int_t index=pt->GetClusterIndex2(j);
2590 if (index<0) continue;
2591 AliTPCclusterMI *c= pt->GetClusterPointer(j);
2593 // if (!(c->IsUsed(10))) c->Use();
2598 // gLastCheck = nseed;
2606 void AliTPCtrackerMI::StopNotActive(const TObjArray * arr, Int_t row0, Float_t th0, Float_t th1, Float_t th2) const
2608 // stop not active tracks
2609 // take th1 as threshold for number of founded to number of foundable on last 10 active rows
2610 // take th2 as threshold for number of founded to number of foundable on last 20 active rows
2611 Int_t nseed = arr->GetEntriesFast();
2613 for (Int_t i=0; i<nseed; i++) {
2614 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2618 if (!(pt->IsActive())) continue;
2619 StopNotActive(pt,row0,th0, th1,th2);
2625 void AliTPCtrackerMI::StopNotActive(AliTPCseed * seed, Int_t row0, Float_t th0, Float_t th1,
2628 // stop not active tracks
2629 // take th1 as threshold for number of founded to number of foundable on last 10 active rows
2630 // take th2 as threshold for number of founded to number of foundable on last 20 active rows
2633 Int_t foundable = 0;
2634 Int_t maxindex = seed->GetLastPoint(); //last foundable row
2635 if (seed->GetNFoundable()*th0 > seed->GetNumberOfClusters()) {
2636 seed->Desactivate(10) ;
2640 for (Int_t i=row0; i<maxindex; i++){
2641 Int_t index = seed->GetClusterIndex2(i);
2642 if (index!=-1) foundable++;
2644 if (foundable<=30) sumgood1++;
2645 if (foundable<=50) {
2652 if (foundable>=30.){
2653 if (sumgood1<(th1*30.)) seed->Desactivate(10);
2656 if (sumgood2<(th2*50.)) seed->Desactivate(10);
2660 Int_t AliTPCtrackerMI::RefitInward(AliESDEvent *event)
2663 // back propagation of ESD tracks
2666 if (!event) return 0;
2667 const Int_t kMaxFriendTracks=2000;
2669 // extract correction object for multiplicity dependence of dEdx
2670 TObjArray * gainCalibArray = AliTPCcalibDB::Instance()->GetTimeGainSplinesRun(event->GetRunNumber());
2672 AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ;
2674 AliFatal("Tranformations not in RefitInward");
2677 transform->SetCurrentRecoParam((AliTPCRecoParam*)AliTPCReconstructor::GetRecoParam());
2678 const AliTPCRecoParam * recoParam = AliTPCcalibDB::Instance()->GetTransform()->GetCurrentRecoParam();
2680 TGraphErrors * graphMultDependenceDeDx = 0x0;
2681 if (recoParam && recoParam->GetUseMultiplicityCorrectionDedx() && gainCalibArray) {
2682 if (recoParam->GetUseTotCharge()) {
2683 graphMultDependenceDeDx = (TGraphErrors *) gainCalibArray->FindObject("TGRAPHERRORS_MEANQTOT_MULTIPLICITYDEPENDENCE_BEAM_ALL");
2685 graphMultDependenceDeDx = (TGraphErrors *) gainCalibArray->FindObject("TGRAPHERRORS_MEANQMAX_MULTIPLICITYDEPENDENCE_BEAM_ALL");
2691 //PrepareForProlongation(fSeeds,1);
2692 PropagateForward2(fSeeds);
2693 RemoveUsed2(fSeeds,0.4,0.4,20);
2695 TObjArray arraySeed(fSeeds->GetEntries());
2696 for (Int_t i=0;i<fSeeds->GetEntries();i++) {
2697 arraySeed.AddAt(fSeeds->At(i),i);
2699 SignShared(&arraySeed);
2700 // FindCurling(fSeeds, event,2); // find multi found tracks
2701 FindSplitted(fSeeds, event,2); // find multi found tracks
2702 if (AliTPCReconstructor::StreamLevel()>5) FindMultiMC(fSeeds, fEvent,2); // find multi found tracks
2705 Int_t nseed = fSeeds->GetEntriesFast();
2706 for (Int_t i=0;i<nseed;i++){
2707 AliTPCseed * seed = (AliTPCseed*) fSeeds->UncheckedAt(i);
2708 if (!seed) continue;
2709 if (seed->GetKinkIndex(0)>0) UpdateKinkQualityD(seed); // update quality informations for kinks
2710 AliESDtrack *esd=event->GetTrack(i);
2712 if (seed->GetNumberOfClusters()<60 && seed->GetNumberOfClusters()<(esd->GetTPCclusters(0) -5)*0.8){
2713 AliExternalTrackParam paramIn;
2714 AliExternalTrackParam paramOut;
2715 Int_t ncl = seed->RefitTrack(seed,¶mIn,¶mOut);
2716 if (AliTPCReconstructor::StreamLevel()>2) {
2717 (*fDebugStreamer)<<"RecoverIn"<<
2721 "pout.="<<¶mOut<<
2726 seed->Set(paramIn.GetX(),paramIn.GetAlpha(),paramIn.GetParameter(),paramIn.GetCovariance());
2727 seed->SetNumberOfClusters(ncl);
2731 seed->PropagateTo(fkParam->GetInnerRadiusLow());
2732 seed->UpdatePoints();
2733 AddCovariance(seed);
2735 seed->CookdEdx(0.02,0.6);
2736 CookLabel(seed,0.1); //For comparison only
2737 esd->SetTPCClusterMap(seed->GetClusterMap());
2738 esd->SetTPCSharedMap(seed->GetSharedMap());
2740 if (AliTPCReconstructor::StreamLevel()>1 && seed!=0) {
2741 TTreeSRedirector &cstream = *fDebugStreamer;
2748 if (seed->GetNumberOfClusters()>15){
2749 esd->UpdateTrackParams(seed,AliESDtrack::kTPCrefit);
2750 esd->SetTPCPoints(seed->GetPoints());
2751 esd->SetTPCPointsF(seed->GetNFoundable());
2752 Int_t ndedx = seed->GetNCDEDX(0);
2753 Float_t sdedx = seed->GetSDEDX(0);
2754 Float_t dedx = seed->GetdEdx();
2755 // apply mutliplicity dependent dEdx correction if available
2756 if (graphMultDependenceDeDx) {
2757 Int_t nContribut = event->GetPrimaryVertexTPC()->GetNContributors();
2758 Double_t corrGain = AliTPCcalibDButil::EvalGraphConst(graphMultDependenceDeDx, nContribut);
2761 esd->SetTPCsignal(dedx, sdedx, ndedx);
2763 // add seed to the esd track in Calib level
2765 Bool_t storeFriend = gRandom->Rndm()<(kMaxFriendTracks)/Float_t(nseed);
2766 if (AliTPCReconstructor::StreamLevel()>0 &&storeFriend){
2767 AliTPCseed * seedCopy = new AliTPCseed(*seed, kTRUE);
2768 esd->AddCalibObject(seedCopy);
2773 //printf("problem\n");
2776 //FindKinks(fSeeds,event);
2777 if (AliTPCReconstructor::StreamLevel()>3) DumpClusters(2,fSeeds);
2778 Info("RefitInward","Number of refitted tracks %d",ntracks);
2783 Int_t AliTPCtrackerMI::PropagateBack(AliESDEvent *event)
2786 // back propagation of ESD tracks
2788 if (!event) return 0;
2792 PropagateBack(fSeeds);
2793 RemoveUsed2(fSeeds,0.4,0.4,20);
2794 //FindCurling(fSeeds, fEvent,1);
2795 FindSplitted(fSeeds, event,1); // find multi found tracks
2796 if (AliTPCReconstructor::StreamLevel()>5) FindMultiMC(fSeeds, fEvent,1); // find multi found tracks
2799 Int_t nseed = fSeeds->GetEntriesFast();
2801 for (Int_t i=0;i<nseed;i++){
2802 AliTPCseed * seed = (AliTPCseed*) fSeeds->UncheckedAt(i);
2803 if (!seed) continue;
2804 if (seed->GetKinkIndex(0)<0) UpdateKinkQualityM(seed); // update quality informations for kinks
2805 seed->UpdatePoints();
2806 AddCovariance(seed);
2807 AliESDtrack *esd=event->GetTrack(i);
2808 if (seed->GetNumberOfClusters()<60 && seed->GetNumberOfClusters()<(esd->GetTPCclusters(0) -5)*0.8){
2809 AliExternalTrackParam paramIn;
2810 AliExternalTrackParam paramOut;
2811 Int_t ncl = seed->RefitTrack(seed,¶mIn,¶mOut);
2812 if (AliTPCReconstructor::StreamLevel()>2) {
2813 (*fDebugStreamer)<<"RecoverBack"<<
2817 "pout.="<<¶mOut<<
2822 seed->Set(paramOut.GetX(),paramOut.GetAlpha(),paramOut.GetParameter(),paramOut.GetCovariance());
2823 seed->SetNumberOfClusters(ncl);
2826 seed->CookdEdx(0.02,0.6);
2827 CookLabel(seed,0.1); //For comparison only
2828 if (seed->GetNumberOfClusters()>15){
2829 esd->UpdateTrackParams(seed,AliESDtrack::kTPCout);
2830 esd->SetTPCPoints(seed->GetPoints());
2831 esd->SetTPCPointsF(seed->GetNFoundable());
2832 Int_t ndedx = seed->GetNCDEDX(0);
2833 Float_t sdedx = seed->GetSDEDX(0);
2834 Float_t dedx = seed->GetdEdx();
2835 esd->SetTPCsignal(dedx, sdedx, ndedx);
2837 Int_t eventnumber = event->GetEventNumberInFile();// patch 28 fev 06
2838 // This is most likely NOT the event number you'd like to use. It has nothing to do with the 'real' event number
2839 if (AliTPCReconstructor::StreamLevel()>1 && esd) {
2840 (*fDebugStreamer)<<"Cback"<<
2843 "EventNrInFile="<<eventnumber<<
2848 if (AliTPCReconstructor::StreamLevel()>3) DumpClusters(1,fSeeds);
2849 //FindKinks(fSeeds,event);
2850 Info("PropagateBack","Number of back propagated tracks %d",ntracks);
2857 void AliTPCtrackerMI::DeleteSeeds()
2862 Int_t nseed = fSeeds->GetEntriesFast();
2863 for (Int_t i=0;i<nseed;i++){
2864 AliTPCseed * seed = (AliTPCseed*)fSeeds->At(i);
2865 if (seed) delete fSeeds->RemoveAt(i);
2872 void AliTPCtrackerMI::ReadSeeds(AliESDEvent *const event, Int_t direction)
2875 //read seeds from the event
2877 Int_t nentr=event->GetNumberOfTracks();
2879 Info("ReadSeeds", "Number of ESD tracks: %d\n", nentr);
2884 fSeeds = new TObjArray(nentr);
2888 for (Int_t i=0; i<nentr; i++) {
2889 AliESDtrack *esd=event->GetTrack(i);
2890 ULong_t status=esd->GetStatus();
2891 if (!(status&AliESDtrack::kTPCin)) continue;
2892 AliTPCtrack t(*esd);
2893 t.SetNumberOfClusters(0);
2894 // AliTPCseed *seed = new AliTPCseed(t,t.GetAlpha());
2895 AliTPCseed *seed = new AliTPCseed(t/*,t.GetAlpha()*/);
2896 seed->SetUniqueID(esd->GetID());
2897 AddCovariance(seed); //add systematic ucertainty
2898 for (Int_t ikink=0;ikink<3;ikink++) {
2899 Int_t index = esd->GetKinkIndex(ikink);
2900 seed->GetKinkIndexes()[ikink] = index;
2901 if (index==0) continue;
2902 index = TMath::Abs(index);
2903 AliESDkink * kink = fEvent->GetKink(index-1);
2904 if (kink&&esd->GetKinkIndex(ikink)<0){
2905 if ((status & AliESDtrack::kTRDrefit) != 0) kink->SetStatus(1,2);
2906 if ((status & AliESDtrack::kITSout) != 0) kink->SetStatus(1,0);
2908 if (kink&&esd->GetKinkIndex(ikink)>0){
2909 if ((status & AliESDtrack::kTRDrefit) != 0) kink->SetStatus(1,6);
2910 if ((status & AliESDtrack::kITSout) != 0) kink->SetStatus(1,4);
2914 if (((status&AliESDtrack::kITSout)==0)&&(direction==1)) seed->ResetCovariance(10.);
2915 if ( direction ==2 &&(status & AliESDtrack::kTRDrefit) == 0 ) seed->ResetCovariance(10.);
2916 //if ( direction ==2 && ((status & AliESDtrack::kTPCout) == 0) ) {
2917 // fSeeds->AddAt(0,i);
2921 if ( direction ==2 &&(status & AliESDtrack::kTRDrefit) > 0 ) {
2922 Double_t par0[5],par1[5],alpha,x;
2923 esd->GetInnerExternalParameters(alpha,x,par0);
2924 esd->GetExternalParameters(x,par1);
2925 Double_t delta1 = TMath::Abs(par0[4]-par1[4])/(0.000000001+TMath::Abs(par0[4]+par1[4]));
2926 Double_t delta2 = TMath::Abs(par0[3]-par1[3]);
2928 if (esd->GetTRDncls()>0) trdchi2 = esd->GetTRDchi2()/esd->GetTRDncls();
2929 //reset covariance if suspicious
2930 if ( (delta1>0.1) || (delta2>0.006) ||trdchi2>7.)
2931 seed->ResetCovariance(10.);
2936 // rotate to the local coordinate system
2938 fSectors=fInnerSec; fN=fkNIS;
2939 Double_t alpha=seed->GetAlpha() - fSectors->GetAlphaShift();
2940 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
2941 if (alpha < 0. ) alpha += 2.*TMath::Pi();
2942 Int_t ns=Int_t(alpha/fSectors->GetAlpha())%fN;
2943 alpha =ns*fSectors->GetAlpha() + fSectors->GetAlphaShift();
2944 if (alpha<-TMath::Pi()) alpha += 2*TMath::Pi();
2945 if (alpha>=TMath::Pi()) alpha -= 2*TMath::Pi();
2946 alpha-=seed->GetAlpha();
2947 if (!seed->Rotate(alpha)) {
2953 if (esd->GetKinkIndex(0)<=0){
2954 for (Int_t irow=0;irow<160;irow++){
2955 Int_t index = seed->GetClusterIndex2(irow);
2958 AliTPCclusterMI * cl = GetClusterMI(index);
2959 seed->SetClusterPointer(irow,cl);
2961 if ((index & 0x8000)==0){
2962 cl->Use(10); // accepted cluster
2964 cl->Use(6); // close cluster not accepted
2967 Info("ReadSeeds","Not found cluster");
2972 fSeeds->AddAt(seed,i);
2978 //_____________________________________________________________________________
2979 void AliTPCtrackerMI::MakeSeeds3(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, Float_t cuts[4],
2980 Float_t deltay, Int_t ddsec) {
2981 //-----------------------------------------------------------------
2982 // This function creates track seeds.
2983 // SEEDING WITH VERTEX CONSTRAIN
2984 //-----------------------------------------------------------------
2985 // cuts[0] - fP4 cut
2986 // cuts[1] - tan(phi) cut
2987 // cuts[2] - zvertex cut
2988 // cuts[3] - fP3 cut
2996 Double_t x[5], c[15];
2997 // Int_t di = i1-i2;
2999 AliTPCseed * seed = new AliTPCseed();
3000 Double_t alpha=fSectors->GetAlpha(), shift=fSectors->GetAlphaShift();
3001 Double_t cs=cos(alpha), sn=sin(alpha);
3003 // Double_t x1 =fOuterSec->GetX(i1);
3004 //Double_t xx2=fOuterSec->GetX(i2);
3006 Double_t x1 =GetXrow(i1);
3007 Double_t xx2=GetXrow(i2);
3009 Double_t x3=GetX(), y3=GetY(), z3=GetZ();
3011 Int_t imiddle = (i2+i1)/2; //middle pad row index
3012 Double_t xm = GetXrow(imiddle); // radius of middle pad-row
3013 const AliTPCtrackerRow& krm=GetRow(sec,imiddle); //middle pad -row
3017 const AliTPCtrackerRow& kr1=GetRow(ns,i1);
3018 Double_t ymax = GetMaxY(i1)-kr1.GetDeadZone()-1.5;
3019 Double_t ymaxm = GetMaxY(imiddle)-kr1.GetDeadZone()-1.5;
3022 // change cut on curvature if it can't reach this layer
3023 // maximal curvature set to reach it
3024 Double_t dvertexmax = TMath::Sqrt((x1-x3)*(x1-x3)+(ymax+5-y3)*(ymax+5-y3));
3025 if (dvertexmax*0.5*cuts[0]>0.85){
3026 cuts[0] = 0.85/(dvertexmax*0.5+1.);
3028 Double_t r2min = 1/(cuts[0]*cuts[0]); //minimal square of radius given by cut
3031 if (deltay>0) ddsec = 0;
3032 // loop over clusters
3033 for (Int_t is=0; is < kr1; is++) {
3035 if (kr1[is]->IsUsed(10)) continue;
3036 Double_t y1=kr1[is]->GetY(), z1=kr1[is]->GetZ();
3037 //if (TMath::Abs(y1)>ymax) continue;
3039 if (deltay>0 && TMath::Abs(ymax-TMath::Abs(y1))> deltay ) continue; // seed only at the edge
3041 // find possible directions
3042 Float_t anglez = (z1-z3)/(x1-x3);
3043 Float_t extraz = z1 - anglez*(x1-xx2); // extrapolated z
3046 //find rotation angles relative to line given by vertex and point 1
3047 Double_t dvertex2 = (x1-x3)*(x1-x3)+(y1-y3)*(y1-y3);
3048 Double_t dvertex = TMath::Sqrt(dvertex2);
3049 Double_t angle13 = TMath::ATan((y1-y3)/(x1-x3));
3050 Double_t cs13 = cos(-angle13), sn13 = sin(-angle13);
3053 // loop over 2 sectors
3059 Double_t dddz1=0; // direction of delta inclination in z axis
3066 for (Int_t dsec = dsec1; dsec<=dsec2;dsec++){
3067 Int_t sec2 = sec + dsec;
3069 // AliTPCtrackerRow& kr2 = fOuterSec[(sec2+fkNOS)%fkNOS][i2];
3070 //AliTPCtrackerRow& kr2m = fOuterSec[(sec2+fkNOS)%fkNOS][imiddle];
3071 AliTPCtrackerRow& kr2 = GetRow((sec2+fkNOS)%fkNOS,i2);
3072 AliTPCtrackerRow& kr2m = GetRow((sec2+fkNOS)%fkNOS,imiddle);
3073 Int_t index1 = TMath::Max(kr2.Find(extraz-0.6-dddz1*TMath::Abs(z1)*0.05)-1,0);
3074 Int_t index2 = TMath::Min(kr2.Find(extraz+0.6+dddz2*TMath::Abs(z1)*0.05)+1,kr2);
3076 // rotation angles to p1-p3
3077 Double_t cs13r = cos(-angle13+dsec*alpha)/dvertex, sn13r = sin(-angle13+dsec*alpha)/dvertex;
3078 Double_t x2, y2, z2;
3080 // Double_t dymax = maxangle*TMath::Abs(x1-xx2);
3083 Double_t dxx0 = (xx2-x3)*cs13r;
3084 Double_t dyy0 = (xx2-x3)*sn13r;
3085 for (Int_t js=index1; js < index2; js++) {
3086 const AliTPCclusterMI *kcl = kr2[js];
3087 if (kcl->IsUsed(10)) continue;
3089 //calcutate parameters
3091 Double_t yy0 = dyy0 +(kcl->GetY()-y3)*cs13r;
3093 if (TMath::Abs(yy0)<0.000001) continue;
3094 Double_t xx0 = dxx0 -(kcl->GetY()-y3)*sn13r;
3095 Double_t y0 = 0.5*(xx0*xx0+yy0*yy0-xx0)/yy0;
3096 Double_t r02 = (0.25+y0*y0)*dvertex2;
3097 //curvature (radius) cut
3098 if (r02<r2min) continue;
3102 Double_t c0 = 1/TMath::Sqrt(r02);
3106 //Double_t dfi0 = 2.*TMath::ASin(dvertex*c0*0.5);
3107 //Double_t dfi1 = 2.*TMath::ASin(TMath::Sqrt(yy0*yy0+(1-xx0)*(1-xx0))*dvertex*c0*0.5);
3108 Double_t dfi0 = 2.*AliTPCFastMath::FastAsin(dvertex*c0*0.5);
3109 Double_t dfi1 = 2.*AliTPCFastMath::FastAsin(TMath::Sqrt(yy0*yy0+(1-xx0)*(1-xx0))*dvertex*c0*0.5);
3112 Double_t z0 = kcl->GetZ();
3113 Double_t zzzz2 = z1-(z1-z3)*dfi1/dfi0;
3114 if (TMath::Abs(zzzz2-z0)>0.5) continue;
3117 Double_t dip = (z1-z0)*c0/dfi1;
3118 Double_t x0 = (0.5*cs13+y0*sn13)*dvertex*c0;
3129 x2= xx2*cs-y2*sn*dsec;
3130 y2=+xx2*sn*dsec+y2*cs;
3140 // do we have cluster at the middle ?
3142 GetProlongation(x1,xm,x,ym,zm);
3144 AliTPCclusterMI * cm=0;
3145 if (TMath::Abs(ym)-ymaxm<0){
3146 cm = krm.FindNearest2(ym,zm,1.0,0.6,dummy);
3147 if ((!cm) || (cm->IsUsed(10))) {
3152 // rotate y1 to system 0
3153 // get state vector in rotated system
3154 Double_t yr1 = (-0.5*sn13+y0*cs13)*dvertex*c0;
3155 Double_t xr2 = x0*cs+yr1*sn*dsec;
3156 Double_t xr[5]={kcl->GetY(),kcl->GetZ(), xr2, dip, c0};
3158 GetProlongation(xx2,xm,xr,ym,zm);
3159 if (TMath::Abs(ym)-ymaxm<0){
3160 cm = kr2m.FindNearest2(ym,zm,1.0,0.6,dummy);
3161 if ((!cm) || (cm->IsUsed(10))) {
3171 dym = ym - cm->GetY();
3172 dzm = zm - cm->GetZ();
3179 Double_t sy1=kr1[is]->GetSigmaY2()*2., sz1=kr1[is]->GetSigmaZ2()*2.;
3180 Double_t sy2=kcl->GetSigmaY2()*2., sz2=kcl->GetSigmaZ2()*2.;
3181 //Double_t sy3=400*3./12., sy=0.1, sz=0.1;
3182 Double_t sy3=25000*x[4]*x[4]+0.1, sy=0.1, sz=0.1;
3183 //Double_t sy3=25000*x[4]*x[4]*60+0.5, sy=0.1, sz=0.1;
3185 Double_t f40=(F1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
3186 Double_t f42=(F1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
3187 Double_t f43=(F1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
3188 Double_t f20=(F2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
3189 Double_t f22=(F2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
3190 Double_t f23=(F2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
3192 Double_t f30=(F3(x1,y1+sy,x2,y2,z1,z2)-x[3])/sy;
3193 Double_t f31=(F3(x1,y1,x2,y2,z1+sz,z2)-x[3])/sz;
3194 Double_t f32=(F3(x1,y1,x2,y2+sy,z1,z2)-x[3])/sy;
3195 Double_t f34=(F3(x1,y1,x2,y2,z1,z2+sz)-x[3])/sz;
3199 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3200 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3201 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3202 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3203 c[13]=f30*sy1*f40+f32*sy2*f42;
3204 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3206 // if (!BuildSeed(kr1[is],kcl,0,x1,x2,x3,x,c)) continue;
3208 UInt_t index=kr1.GetIndex(is);
3209 seed->~AliTPCseed(); // this does not set the pointer to 0...
3210 AliTPCseed *track=new(seed) AliTPCseed(x1, ns*alpha+shift, x, c, index);
3212 track->SetIsSeeding(kTRUE);
3213 track->SetSeed1(i1);
3214 track->SetSeed2(i2);
3215 track->SetSeedType(3);
3219 FollowProlongation(*track, (i1+i2)/2,1);
3220 Int_t foundable,found,shared;
3221 track->GetClusterStatistic((i1+i2)/2,i1, found, foundable, shared, kTRUE);
3222 if ((found<0.55*foundable) || shared>0.5*found || (track->GetSigmaY2()+track->GetSigmaZ2())>0.5){
3224 seed->~AliTPCseed();
3230 FollowProlongation(*track, i2,1);
3234 track->SetBConstrain(1);
3235 // track->fLastPoint = i1+fInnerSec->GetNRows(); // first cluster in track position
3236 track->SetLastPoint(i1); // first cluster in track position
3237 track->SetFirstPoint(track->GetLastPoint());
3239 if (track->GetNumberOfClusters()<(i1-i2)*0.5 ||
3240 track->GetNumberOfClusters() < track->GetNFoundable()*0.6 ||
3241 track->GetNShared()>0.4*track->GetNumberOfClusters() ) {
3243 seed->~AliTPCseed();
3247 // Z VERTEX CONDITION
3248 Double_t zv, bz=GetBz();
3249 if ( !track->GetZAt(0.,bz,zv) ) continue;
3250 if (TMath::Abs(zv-z3)>cuts[2]) {
3251 FollowProlongation(*track, TMath::Max(i2-20,0));
3252 if ( !track->GetZAt(0.,bz,zv) ) continue;
3253 if (TMath::Abs(zv-z3)>cuts[2]){
3254 FollowProlongation(*track, TMath::Max(i2-40,0));
3255 if ( !track->GetZAt(0.,bz,zv) ) continue;
3256 if (TMath::Abs(zv-z3)>cuts[2] &&(track->GetNumberOfClusters() > track->GetNFoundable()*0.7)){
3257 // make seed without constrain
3258 AliTPCseed * track2 = MakeSeed(track,0.2,0.5,1.);
3259 FollowProlongation(*track2, i2,1);
3260 track2->SetBConstrain(kFALSE);
3261 track2->SetSeedType(1);
3262 arr->AddLast(track2);
3264 seed->~AliTPCseed();
3269 seed->~AliTPCseed();
3276 track->SetSeedType(0);
3277 arr->AddLast(track);
3278 seed = new AliTPCseed;
3280 // don't consider other combinations
3281 if (track->GetNumberOfClusters() > track->GetNFoundable()*0.8)
3287 Info("MakeSeeds3","\nSeeding statistic:\t%d\t%d\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin,nout1,nout2);
3293 void AliTPCtrackerMI::MakeSeeds5(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, Float_t cuts[4],
3298 //-----------------------------------------------------------------
3299 // This function creates track seeds.
3300 //-----------------------------------------------------------------
3301 // cuts[0] - fP4 cut
3302 // cuts[1] - tan(phi) cut
3303 // cuts[2] - zvertex cut
3304 // cuts[3] - fP3 cut
3314 Double_t x[5], c[15];
3316 // make temporary seed
3317 AliTPCseed * seed = new AliTPCseed;
3318 Double_t alpha=fOuterSec->GetAlpha(), shift=fOuterSec->GetAlphaShift();
3319 // Double_t cs=cos(alpha), sn=sin(alpha);
3324 Double_t x1 = GetXrow(i1-1);
3325 const AliTPCtrackerRow& kr1=GetRow(sec,i1-1);
3326 Double_t y1max = GetMaxY(i1-1)-kr1.GetDeadZone()-1.5;
3328 Double_t x1p = GetXrow(i1);
3329 const AliTPCtrackerRow& kr1p=GetRow(sec,i1);
3331 Double_t x1m = GetXrow(i1-2);
3332 const AliTPCtrackerRow& kr1m=GetRow(sec,i1-2);
3335 //last 3 padrow for seeding
3336 AliTPCtrackerRow& kr3 = GetRow((sec+fkNOS)%fkNOS,i1-7);
3337 Double_t x3 = GetXrow(i1-7);
3338 // Double_t y3max= GetMaxY(i1-7)-kr3.fDeadZone-1.5;
3340 AliTPCtrackerRow& kr3p = GetRow((sec+fkNOS)%fkNOS,i1-6);
3341 Double_t x3p = GetXrow(i1-6);
3343 AliTPCtrackerRow& kr3m = GetRow((sec+fkNOS)%fkNOS,i1-8);
3344 Double_t x3m = GetXrow(i1-8);
3349 Int_t im = i1-4; //middle pad row index
3350 Double_t xm = GetXrow(im); // radius of middle pad-row
3351 const AliTPCtrackerRow& krm=GetRow(sec,im); //middle pad -row
3352 // Double_t ymmax = GetMaxY(im)-kr1.fDeadZone-1.5;
3355 Double_t deltax = x1-x3;
3356 Double_t dymax = deltax*cuts[1];
3357 Double_t dzmax = deltax*cuts[3];
3359 // loop over clusters
3360 for (Int_t is=0; is < kr1; is++) {
3362 if (kr1[is]->IsUsed(10)) continue;
3363 Double_t y1=kr1[is]->GetY(), z1=kr1[is]->GetZ();
3365 if (deltay>0 && TMath::Abs(y1max-TMath::Abs(y1))> deltay ) continue; // seed only at the edge
3367 Int_t index1 = TMath::Max(kr3.Find(z1-dzmax)-1,0);
3368 Int_t index2 = TMath::Min(kr3.Find(z1+dzmax)+1,kr3);
3374 for (Int_t js=index1; js < index2; js++) {
3375 const AliTPCclusterMI *kcl = kr3[js];
3376 if (kcl->IsUsed(10)) continue;
3378 // apply angular cuts
3379 if (TMath::Abs(y1-y3)>dymax) continue;
3382 if (TMath::Abs(z1-z3)>dzmax) continue;
3384 Double_t angley = (y1-y3)/(x1-x3);
3385 Double_t anglez = (z1-z3)/(x1-x3);
3387 Double_t erry = TMath::Abs(angley)*(x1-x1m)*0.5+0.5;
3388 Double_t errz = TMath::Abs(anglez)*(x1-x1m)*0.5+0.5;
3390 Double_t yyym = angley*(xm-x1)+y1;
3391 Double_t zzzm = anglez*(xm-x1)+z1;
3393 const AliTPCclusterMI *kcm = krm.FindNearest2(yyym,zzzm,erry,errz,index);
3395 if (kcm->IsUsed(10)) continue;
3397 erry = TMath::Abs(angley)*(x1-x1m)*0.4+0.5;
3398 errz = TMath::Abs(anglez)*(x1-x1m)*0.4+0.5;
3405 // look around first
3406 const AliTPCclusterMI *kc1m = kr1m.FindNearest2(angley*(x1m-x1)+y1,
3412 if (kc1m->IsUsed(10)) used++;
3414 const AliTPCclusterMI *kc1p = kr1p.FindNearest2(angley*(x1p-x1)+y1,
3420 if (kc1p->IsUsed(10)) used++;
3422 if (used>1) continue;
3423 if (found<1) continue;
3427 const AliTPCclusterMI *kc3m = kr3m.FindNearest2(angley*(x3m-x3)+y3,
3433 if (kc3m->IsUsed(10)) used++;
3437 const AliTPCclusterMI *kc3p = kr3p.FindNearest2(angley*(x3p-x3)+y3,
3443 if (kc3p->IsUsed(10)) used++;
3447 if (used>1) continue;
3448 if (found<3) continue;
3458 x[4]=F1(x1,y1,x2,y2,x3,y3);
3459 //if (TMath::Abs(x[4]) >= cuts[0]) continue;
3462 x[2]=F2(x1,y1,x2,y2,x3,y3);
3465 x[3]=F3n(x1,y1,x2,y2,z1,z2,x[4]);
3466 //if (TMath::Abs(x[3]) > cuts[3]) continue;
3470 Double_t sy1=0.1, sz1=0.1;
3471 Double_t sy2=0.1, sz2=0.1;
3472 Double_t sy3=0.1, sy=0.1, sz=0.1;
3474 Double_t f40=(F1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
3475 Double_t f42=(F1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
3476 Double_t f43=(F1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
3477 Double_t f20=(F2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
3478 Double_t f22=(F2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
3479 Double_t f23=(F2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
3481 Double_t f30=(F3(x1,y1+sy,x2,y2,z1,z2)-x[3])/sy;
3482 Double_t f31=(F3(x1,y1,x2,y2,z1+sz,z2)-x[3])/sz;
3483 Double_t f32=(F3(x1,y1,x2,y2+sy,z1,z2)-x[3])/sy;
3484 Double_t f34=(F3(x1,y1,x2,y2,z1,z2+sz)-x[3])/sz;
3488 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3489 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3490 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3491 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3492 c[13]=f30*sy1*f40+f32*sy2*f42;
3493 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3495 // if (!BuildSeed(kr1[is],kcl,0,x1,x2,x3,x,c)) continue;
3497 index=kr1.GetIndex(is);
3498 seed->~AliTPCseed();
3499 AliTPCseed *track=new(seed) AliTPCseed(x1, sec*alpha+shift, x, c, index);
3501 track->SetIsSeeding(kTRUE);
3504 FollowProlongation(*track, i1-7,1);
3505 if (track->GetNumberOfClusters() < track->GetNFoundable()*0.75 ||
3506 track->GetNShared()>0.6*track->GetNumberOfClusters() || ( track->GetSigmaY2()+ track->GetSigmaZ2())>0.6){
3508 seed->~AliTPCseed();
3514 FollowProlongation(*track, i2,1);
3515 track->SetBConstrain(0);
3516 track->SetLastPoint(i1+fInnerSec->GetNRows()); // first cluster in track position
3517 track->SetFirstPoint(track->GetLastPoint());
3519 if (track->GetNumberOfClusters()<(i1-i2)*0.5 ||
3520 track->GetNumberOfClusters()<track->GetNFoundable()*0.7 ||
3521 track->GetNShared()>2. || track->GetChi2()/track->GetNumberOfClusters()>6 || ( track->GetSigmaY2()+ track->GetSigmaZ2())>0.5 ) {
3523 seed->~AliTPCseed();
3528 FollowProlongation(*track, TMath::Max(i2-10,0),1);
3529 AliTPCseed * track2 = MakeSeed(track,0.2,0.5,0.9);
3530 FollowProlongation(*track2, i2,1);
3531 track2->SetBConstrain(kFALSE);
3532 track2->SetSeedType(4);
3533 arr->AddLast(track2);
3535 seed->~AliTPCseed();
3539 //arr->AddLast(track);
3540 //seed = new AliTPCseed;
3546 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);
3552 //_____________________________________________________________________________
3553 void AliTPCtrackerMI::MakeSeeds2(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, Float_t */*cuts[4]*/,
3554 Float_t deltay, Bool_t /*bconstrain*/) {
3555 //-----------------------------------------------------------------
3556 // This function creates track seeds - without vertex constraint
3557 //-----------------------------------------------------------------
3558 // cuts[0] - fP4 cut - not applied
3559 // cuts[1] - tan(phi) cut
3560 // cuts[2] - zvertex cut - not applied
3561 // cuts[3] - fP3 cut
3571 Double_t alpha=fOuterSec->GetAlpha(), shift=fOuterSec->GetAlphaShift();
3572 // Double_t cs=cos(alpha), sn=sin(alpha);
3573 Int_t row0 = (i1+i2)/2;
3574 Int_t drow = (i1-i2)/2;
3575 const AliTPCtrackerRow& kr0=fSectors[sec][row0];
3576 AliTPCtrackerRow * kr=0;
3578 AliTPCpolyTrack polytrack;
3579 Int_t nclusters=fSectors[sec][row0];
3580 AliTPCseed * seed = new AliTPCseed;
3585 for (Int_t is=0; is < nclusters; is++) { //LOOP over clusters
3587 Int_t nfoundable =0;
3588 for (Int_t iter =1; iter<2; iter++){ //iterations
3589 const AliTPCtrackerRow& krm=fSectors[sec][row0-iter];
3590 const AliTPCtrackerRow& krp=fSectors[sec][row0+iter];
3591 const AliTPCclusterMI * cl= kr0[is];
3593 if (cl->IsUsed(10)) {
3599 Double_t x = kr0.GetX();
3600 // Initialization of the polytrack
3605 Double_t y0= cl->GetY();
3606 Double_t z0= cl->GetZ();
3610 Double_t ymax = fSectors->GetMaxY(row0)-kr0.GetDeadZone()-1.5;
3611 if (deltay>0 && TMath::Abs(ymax-TMath::Abs(y0))> deltay ) continue; // seed only at the edge
3613 erry = (0.5)*cl->GetSigmaY2()/TMath::Sqrt(cl->GetQ())*6;
3614 errz = (0.5)*cl->GetSigmaZ2()/TMath::Sqrt(cl->GetQ())*6;
3615 polytrack.AddPoint(x,y0,z0,erry, errz);
3618 if (cl->IsUsed(10)) sumused++;
3621 Float_t roady = (5*TMath::Sqrt(cl->GetSigmaY2()+0.2)+1.)*iter;
3622 Float_t roadz = (5*TMath::Sqrt(cl->GetSigmaZ2()+0.2)+1.)*iter;
3625 AliTPCclusterMI * cl1 = krm.FindNearest(y0,z0,roady,roadz);
3626 if (cl1 && TMath::Abs(ymax-TMath::Abs(y0))) {
3627 erry = (0.5)*cl1->GetSigmaY2()/TMath::Sqrt(cl1->GetQ())*3;
3628 errz = (0.5)*cl1->GetSigmaZ2()/TMath::Sqrt(cl1->GetQ())*3;
3629 if (cl1->IsUsed(10)) sumused++;
3630 polytrack.AddPoint(x,cl1->GetY(),cl1->GetZ(),erry,errz);
3634 AliTPCclusterMI * cl2 = krp.FindNearest(y0,z0,roady,roadz);
3636 erry = (0.5)*cl2->GetSigmaY2()/TMath::Sqrt(cl2->GetQ())*3;
3637 errz = (0.5)*cl2->GetSigmaZ2()/TMath::Sqrt(cl2->GetQ())*3;
3638 if (cl2->IsUsed(10)) sumused++;
3639 polytrack.AddPoint(x,cl2->GetY(),cl2->GetZ(),erry,errz);
3642 if (sumused>0) continue;
3644 polytrack.UpdateParameters();
3650 nfoundable = polytrack.GetN();
3651 nfound = nfoundable;
3653 for (Int_t ddrow = iter+1; ddrow<drow;ddrow++){
3654 Float_t maxdist = 0.8*(1.+3./(ddrow));
3655 for (Int_t delta = -1;delta<=1;delta+=2){
3656 Int_t row = row0+ddrow*delta;
3657 kr = &(fSectors[sec][row]);
3658 Double_t xn = kr->GetX();
3659 Double_t ymax1 = fSectors->GetMaxY(row)-kr->GetDeadZone()-1.5;
3660 polytrack.GetFitPoint(xn,yn,zn);
3661 if (TMath::Abs(yn)>ymax1) continue;
3663 AliTPCclusterMI * cln = kr->FindNearest(yn,zn,roady,roadz);
3665 Float_t dist = TMath::Sqrt( (yn-cln->GetY())*(yn-cln->GetY())+(zn-cln->GetZ())*(zn-cln->GetZ()));
3668 erry = (dist+0.3)*cln->GetSigmaY2()/TMath::Sqrt(cln->GetQ())*(1.+1./(ddrow));
3669 errz = (dist+0.3)*cln->GetSigmaZ2()/TMath::Sqrt(cln->GetQ())*(1.+1./(ddrow));
3670 if (cln->IsUsed(10)) {
3671 // printf("used\n");
3679 polytrack.AddPoint(xn,cln->GetY(),cln->GetZ(),erry, errz);
3684 if ( (sumused>3) || (sumused>0.5*nfound) || (nfound<0.6*nfoundable)) break;
3685 polytrack.UpdateParameters();
3688 if ( (sumused>3) || (sumused>0.5*nfound)) {
3689 //printf("sumused %d\n",sumused);
3694 polytrack.GetFitDerivation(kr0.GetX(),dy,dz);
3695 AliTPCpolyTrack track2;
3697 polytrack.Refit(track2,0.5+TMath::Abs(dy)*0.3,0.4+TMath::Abs(dz)*0.3);
3698 if (track2.GetN()<0.5*nfoundable) continue;
3701 if ((nfound>0.6*nfoundable) &&( nfoundable>0.4*(i1-i2))) {
3703 // test seed with and without constrain
3704 for (Int_t constrain=0; constrain<=0;constrain++){
3705 // add polytrack candidate
3707 Double_t x[5], c[15];
3708 Double_t x1,x2,x3,y1,y2,y3,z1,z2,z3;
3709 track2.GetBoundaries(x3,x1);
3711 track2.GetFitPoint(x1,y1,z1);
3712 track2.GetFitPoint(x2,y2,z2);
3713 track2.GetFitPoint(x3,y3,z3);
3715 //is track pointing to the vertex ?
3718 polytrack.GetFitPoint(x0,y0,z0);
3731 x[4]=F1(x1,y1,x2,y2,x3,y3);
3733 // if (TMath::Abs(x[4]) >= cuts[0]) continue; //
3734 x[2]=F2(x1,y1,x2,y2,x3,y3);
3736 //if (TMath::Abs(x[4]*x1-x[2]) >= cuts[1]) continue;
3737 //x[3]=F3(x1,y1,x2,y2,z1,z2);
3738 x[3]=F3n(x1,y1,x3,y3,z1,z3,x[4]);
3739 //if (TMath::Abs(x[3]) > cuts[3]) continue;
3742 Double_t sy =0.1, sz =0.1;
3743 Double_t sy1=0.02, sz1=0.02;
3744 Double_t sy2=0.02, sz2=0.02;
3748 sy3=25000*x[4]*x[4]+0.1, sy=0.1, sz=0.1;
3751 Double_t f40=(F1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
3752 Double_t f42=(F1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
3753 Double_t f43=(F1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
3754 Double_t f20=(F2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
3755 Double_t f22=(F2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
3756 Double_t f23=(F2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
3758 Double_t f30=(F3(x1,y1+sy,x3,y3,z1,z3)-x[3])/sy;
3759 Double_t f31=(F3(x1,y1,x3,y3,z1+sz,z3)-x[3])/sz;
3760 Double_t f32=(F3(x1,y1,x3,y3+sy,z1,z3)-x[3])/sy;
3761 Double_t f34=(F3(x1,y1,x3,y3,z1,z3+sz)-x[3])/sz;
3766 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3767 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3768 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3769 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3770 c[13]=f30*sy1*f40+f32*sy2*f42;
3771 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3773 //Int_t row1 = fSectors->GetRowNumber(x1);
3774 Int_t row1 = GetRowNumber(x1);
3778 seed->~AliTPCseed();
3779 AliTPCseed *track=new(seed) AliTPCseed(x1,sec*alpha+shift,x,c,index);
3780 track->SetIsSeeding(kTRUE);
3781 Int_t rc=FollowProlongation(*track, i2);
3782 if (constrain) track->SetBConstrain(1);
3784 track->SetBConstrain(0);
3785 track->SetLastPoint(row1+fInnerSec->GetNRows()); // first cluster in track position
3786 track->SetFirstPoint(track->GetLastPoint());
3788 if (rc==0 || track->GetNumberOfClusters()<(i1-i2)*0.5 ||
3789 track->GetNumberOfClusters() < track->GetNFoundable()*0.6 ||
3790 track->GetNShared()>0.4*track->GetNumberOfClusters()) {
3793 seed->~AliTPCseed();
3796 arr->AddLast(track);
3797 seed = new AliTPCseed;
3801 } // if accepted seed
3804 Info("MakeSeeds2","\nSeeding statiistic:\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin3);
3810 AliTPCseed *AliTPCtrackerMI::MakeSeed(AliTPCseed *const track, Float_t r0, Float_t r1, Float_t r2)
3814 //reseed using track points
3815 Int_t p0 = int(r0*track->GetNumberOfClusters()); // point 0
3816 Int_t p1 = int(r1*track->GetNumberOfClusters());
3817 Int_t p2 = int(r2*track->GetNumberOfClusters()); // last point
3819 Double_t x0[3],x1[3],x2[3];
3820 for (Int_t i=0;i<3;i++){
3826 // find track position at given ratio of the length
3827 Int_t sec0=0, sec1=0, sec2=0;
3830 for (Int_t i=0;i<160;i++){
3831 if (track->GetClusterPointer(i)){
3833 AliTPCTrackerPoint *trpoint =track->GetTrackPoint(i);
3834 if ( (index<p0) || x0[0]<0 ){
3835 if (trpoint->GetX()>1){
3836 clindex = track->GetClusterIndex2(i);
3838 x0[0] = trpoint->GetX();
3839 x0[1] = trpoint->GetY();
3840 x0[2] = trpoint->GetZ();
3841 sec0 = ((clindex&0xff000000)>>24)%18;
3846 if ( (index<p1) &&(trpoint->GetX()>1)){
3847 clindex = track->GetClusterIndex2(i);
3849 x1[0] = trpoint->GetX();
3850 x1[1] = trpoint->GetY();
3851 x1[2] = trpoint->GetZ();
3852 sec1 = ((clindex&0xff000000)>>24)%18;
3855 if ( (index<p2) &&(trpoint->GetX()>1)){
3856 clindex = track->GetClusterIndex2(i);
3858 x2[0] = trpoint->GetX();
3859 x2[1] = trpoint->GetY();
3860 x2[2] = trpoint->GetZ();
3861 sec2 = ((clindex&0xff000000)>>24)%18;
3868 Double_t alpha, cs,sn, xx2,yy2;
3870 alpha = (sec1-sec2)*fSectors->GetAlpha();
3871 cs = TMath::Cos(alpha);
3872 sn = TMath::Sin(alpha);
3873 xx2= x1[0]*cs-x1[1]*sn;
3874 yy2= x1[0]*sn+x1[1]*cs;
3878 alpha = (sec0-sec2)*fSectors->GetAlpha();
3879 cs = TMath::Cos(alpha);
3880 sn = TMath::Sin(alpha);
3881 xx2= x0[0]*cs-x0[1]*sn;
3882 yy2= x0[0]*sn+x0[1]*cs;
3888 Double_t x[5],c[15];
3892 x[4]=F1(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]);
3893 // if (x[4]>1) return 0;
3894 x[2]=F2(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]);
3895 x[3]=F3n(x2[0],x2[1],x0[0],x0[1],x2[2],x0[2],x[4]);
3896 //if (TMath::Abs(x[3]) > 2.2) return 0;
3897 //if (TMath::Abs(x[2]) > 1.99) return 0;
3899 Double_t sy =0.1, sz =0.1;
3901 Double_t sy1=0.02+track->GetSigmaY2(), sz1=0.02+track->GetSigmaZ2();
3902 Double_t sy2=0.01+track->GetSigmaY2(), sz2=0.01+track->GetSigmaZ2();
3903 Double_t sy3=0.01+track->GetSigmaY2();
3905 Double_t f40=(F1(x2[0],x2[1]+sy,x1[0],x1[1],x0[0],x0[1])-x[4])/sy;
3906 Double_t f42=(F1(x2[0],x2[1],x1[0],x1[1]+sy,x0[0],x0[1])-x[4])/sy;
3907 Double_t f43=(F1(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]+sy)-x[4])/sy;
3908 Double_t f20=(F2(x2[0],x2[1]+sy,x1[0],x1[1],x0[0],x0[1])-x[2])/sy;
3909 Double_t f22=(F2(x2[0],x2[1],x1[0],x1[1]+sy,x0[0],x0[1])-x[2])/sy;
3910 Double_t f23=(F2(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]+sy)-x[2])/sy;
3912 Double_t f30=(F3(x2[0],x2[1]+sy,x0[0],x0[1],x2[2],x0[2])-x[3])/sy;
3913 Double_t f31=(F3(x2[0],x2[1],x0[0],x0[1],x2[2]+sz,x0[2])-x[3])/sz;
3914 Double_t f32=(F3(x2[0],x2[1],x0[0],x0[1]+sy,x2[2],x0[2])-x[3])/sy;
3915 Double_t f34=(F3(x2[0],x2[1],x0[0],x0[1],x2[2],x0[2]+sz)-x[3])/sz;
3920 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3921 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3922 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3923 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3924 c[13]=f30*sy1*f40+f32*sy2*f42;
3925 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3927 // Int_t row1 = fSectors->GetRowNumber(x2[0]);
3928 AliTPCseed *seed=new AliTPCseed(x2[0], sec2*fSectors->GetAlpha()+fSectors->GetAlphaShift(), x, c, 0);
3929 // Double_t y0,z0,y1,z1, y2,z2;
3930 //seed->GetProlongation(x0[0],y0,z0);
3931 // seed->GetProlongation(x1[0],y1,z1);
3932 //seed->GetProlongation(x2[0],y2,z2);
3934 seed->SetLastPoint(pp2);
3935 seed->SetFirstPoint(pp2);
3942 AliTPCseed *AliTPCtrackerMI::ReSeed(const AliTPCseed *track, Float_t r0, Float_t r1, Float_t r2)
3946 //reseed using founded clusters
3948 // Find the number of clusters
3949 Int_t nclusters = 0;
3950 for (Int_t irow=0;irow<160;irow++){
3951 if (track->GetClusterIndex(irow)>0) nclusters++;
3955 ipos[0] = TMath::Max(int(r0*nclusters),0); // point 0 cluster
3956 ipos[1] = TMath::Min(int(r1*nclusters),nclusters-1); //
3957 ipos[2] = TMath::Min(int(r2*nclusters),nclusters-1); // last point
3960 Double_t xyz[3][3]={{0}};
3961 Int_t row[3]={0},sec[3]={0,0,0};
3963 // find track row position at given ratio of the length
3965 for (Int_t irow=0;irow<160;irow++){
3966 if (track->GetClusterIndex2(irow)<0) continue;
3968 for (Int_t ipoint=0;ipoint<3;ipoint++){
3969 if (index<=ipos[ipoint]) row[ipoint] = irow;
3973 //Get cluster and sector position
3974 for (Int_t ipoint=0;ipoint<3;ipoint++){
3975 Int_t clindex = track->GetClusterIndex2(row[ipoint]);
3976 AliTPCclusterMI * cl = GetClusterMI(clindex);
3979 // AliTPCclusterMI * cl = GetClusterMI(clindex);
3982 sec[ipoint] = ((clindex&0xff000000)>>24)%18;
3983 xyz[ipoint][0] = GetXrow(row[ipoint]);
3984 xyz[ipoint][1] = cl->GetY();
3985 xyz[ipoint][2] = cl->GetZ();
3989 // Calculate seed state vector and covariance matrix
3991 Double_t alpha, cs,sn, xx2,yy2;
3993 alpha = (sec[1]-sec[2])*fSectors->GetAlpha();
3994 cs = TMath::Cos(alpha);
3995 sn = TMath::Sin(alpha);
3996 xx2= xyz[1][0]*cs-xyz[1][1]*sn;
3997 yy2= xyz[1][0]*sn+xyz[1][1]*cs;
4001 alpha = (sec[0]-sec[2])*fSectors->GetAlpha();
4002 cs = TMath::Cos(alpha);
4003 sn = TMath::Sin(alpha);
4004 xx2= xyz[0][0]*cs-xyz[0][1]*sn;
4005 yy2= xyz[0][0]*sn+xyz[0][1]*cs;
4011 Double_t x[5],c[15];
4015 x[4]=F1(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
4016 x[2]=F2(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
4017 x[3]=F3n(xyz[2][0],xyz[2][1],xyz[0][0],xyz[0][1],xyz[2][2],xyz[0][2],x[4]);
4019 Double_t sy =0.1, sz =0.1;
4021 Double_t sy1=0.2, sz1=0.2;
4022 Double_t sy2=0.2, sz2=0.2;
4025 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;
4026 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;
4027 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;
4028 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;
4029 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;
4030 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;
4032 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;
4033 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;
4034 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;
4035 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;
4040 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
4041 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
4042 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
4043 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
4044 c[13]=f30*sy1*f40+f32*sy2*f42;
4045 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
4047 // Int_t row1 = fSectors->GetRowNumber(xyz[2][0]);
4048 AliTPCseed *seed=new AliTPCseed(xyz[2][0], sec[2]*fSectors->GetAlpha()+fSectors->GetAlphaShift(), x, c, 0);
4049 seed->SetLastPoint(row[2]);
4050 seed->SetFirstPoint(row[2]);
4055 AliTPCseed *AliTPCtrackerMI::ReSeed(AliTPCseed *track,Int_t r0, Bool_t forward)
4059 //reseed using founded clusters
4062 Int_t row[3]={0,0,0};
4063 Int_t sec[3]={0,0,0};
4065 // forward direction
4067 for (Int_t irow=r0;irow<160;irow++){
4068 if (track->GetClusterIndex(irow)>0){
4073 for (Int_t irow=160;irow>r0;irow--){
4074 if (track->GetClusterIndex(irow)>0){
4079 for (Int_t irow=row[2]-15;irow>row[0];irow--){
4080 if (track->GetClusterIndex(irow)>0){
4088 for (Int_t irow=0;irow<r0;irow++){
4089 if (track->GetClusterIndex(irow)>0){
4094 for (Int_t irow=r0;irow>0;irow--){
4095 if (track->GetClusterIndex(irow)>0){
4100 for (Int_t irow=row[2]-15;irow>row[0];irow--){
4101 if (track->GetClusterIndex(irow)>0){
4108 if ((row[2]-row[0])<20) return 0;
4109 if (row[1]==0) return 0;
4112 //Get cluster and sector position
4113 for (Int_t ipoint=0;ipoint<3;ipoint++){
4114 Int_t clindex = track->GetClusterIndex2(row[ipoint]);
4115 AliTPCclusterMI * cl = GetClusterMI(clindex);
4118 // AliTPCclusterMI * cl = GetClusterMI(clindex);
4121 sec[ipoint] = ((clindex&0xff000000)>>24)%18;
4122 xyz[ipoint][0] = GetXrow(row[ipoint]);
4123 AliTPCTrackerPoint * point = track->GetTrackPoint(row[ipoint]);
4124 if (point&&ipoint<2){
4126 xyz[ipoint][1] = point->GetY();
4127 xyz[ipoint][2] = point->GetZ();
4130 xyz[ipoint][1] = cl->GetY();
4131 xyz[ipoint][2] = cl->GetZ();
4138 // Calculate seed state vector and covariance matrix
4140 Double_t alpha, cs,sn, xx2,yy2;
4142 alpha = (sec[1]-sec[2])*fSectors->GetAlpha();
4143 cs = TMath::Cos(alpha);
4144 sn = TMath::Sin(alpha);
4145 xx2= xyz[1][0]*cs-xyz[1][1]*sn;
4146 yy2= xyz[1][0]*sn+xyz[1][1]*cs;
4150 alpha = (sec[0]-sec[2])*fSectors->GetAlpha();
4151 cs = TMath::Cos(alpha);
4152 sn = TMath::Sin(alpha);
4153 xx2= xyz[0][0]*cs-xyz[0][1]*sn;
4154 yy2= xyz[0][0]*sn+xyz[0][1]*cs;
4160 Double_t x[5],c[15];
4164 x[4]=F1(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
4165 x[2]=F2(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
4166 x[3]=F3n(xyz[2][0],xyz[2][1],xyz[0][0],xyz[0][1],xyz[2][2],xyz[0][2],x[4]);
4168 Double_t sy =0.1, sz =0.1;
4170 Double_t sy1=0.2, sz1=0.2;
4171 Double_t sy2=0.2, sz2=0.2;
4174 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;
4175 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;
4176 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;
4177 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;
4178 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;
4179 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;
4181 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;
4182 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;
4183 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;
4184 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;
4189 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
4190 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
4191 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
4192 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
4193 c[13]=f30*sy1*f40+f32*sy2*f42;
4194 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
4196 // Int_t row1 = fSectors->GetRowNumber(xyz[2][0]);
4197 AliTPCseed *seed=new AliTPCseed(xyz[2][0], sec[2]*fSectors->GetAlpha()+fSectors->GetAlphaShift(), x, c, 0);
4198 seed->SetLastPoint(row[2]);
4199 seed->SetFirstPoint(row[2]);
4200 for (Int_t i=row[0];i<row[2];i++){
4201 seed->SetClusterIndex(i, track->GetClusterIndex(i));
4209 void AliTPCtrackerMI::FindMultiMC(const TObjArray * array, AliESDEvent */*esd*/, Int_t iter)
4212 // find multi tracks - THIS FUNCTION IS ONLY FOR DEBUG PURPOSES
4214 // Use AliTPCReconstructor::StreamLevel()>2 if you want to tune parameters - cuts
4216 // Two reasons to have multiple find tracks
4217 // 1. Curling tracks can be find more than once
4218 // 2. Splitted tracks
4219 // a.) Multiple seeding to increase tracking efficiency - (~ 100% reached)
4220 // b.) Edge effect on the sector boundaries
4223 // Algorithm done in 2 phases - because of CPU consumption
4224 // it is n^2 algorithm - for lead-lead 20000x20000 combination are investigated
4226 // Algorihm for curling tracks sign:
4227 // 1 phase -makes a very rough fast cuts to minimize combinatorics
4228 // a.) opposite sign
4229 // b.) one of the tracks - not pointing to the primary vertex -
4230 // c.) delta tan(theta)
4232 // 2 phase - calculates DCA between tracks - time consument
4237 // General cuts - for splitted tracks and for curling tracks
4239 const Float_t kMaxdPhi = 0.2; // maximal distance in phi
4241 // Curling tracks cuts
4246 Int_t nentries = array->GetEntriesFast();
4247 AliHelix *helixes = new AliHelix[nentries];
4248 Float_t *xm = new Float_t[nentries];
4249 Float_t *dz0 = new Float_t[nentries];
4250 Float_t *dz1 = new Float_t[nentries];
4256 // Find track COG in x direction - point with best defined parameters
4258 for (Int_t i=0;i<nentries;i++){
4259 AliTPCseed* track = (AliTPCseed*)array->At(i);
4260 if (!track) continue;
4261 track->SetCircular(0);
4262 new (&helixes[i]) AliHelix(*track);
4266 track->GetDZ(GetX(),GetY(),GetZ(),GetBz(),dz);
4269 for (Int_t icl=0; icl<160; icl++){
4270 AliTPCclusterMI * cl = track->GetClusterPointer(icl);
4276 if (ncl>0) xm[i]/=Float_t(ncl);
4279 for (Int_t i0=0;i0<nentries;i0++){
4280 AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
4281 if (!track0) continue;
4282 Float_t xc0 = helixes[i0].GetHelix(6);
4283 Float_t yc0 = helixes[i0].GetHelix(7);
4284 Float_t r0 = helixes[i0].GetHelix(8);
4285 Float_t rc0 = TMath::Sqrt(xc0*xc0+yc0*yc0);
4286 Float_t fi0 = TMath::ATan2(yc0,xc0);
4288 for (Int_t i1=i0+1;i1<nentries;i1++){
4289 AliTPCseed * track1 = (AliTPCseed*)array->At(i1);
4290 if (!track1) continue;
4291 Int_t lab0=track0->GetLabel();
4292 Int_t lab1=track1->GetLabel();
4293 if (TMath::Abs(lab0)!=TMath::Abs(lab1)) continue;
4295 Float_t xc1 = helixes[i1].GetHelix(6);
4296 Float_t yc1 = helixes[i1].GetHelix(7);
4297 Float_t r1 = helixes[i1].GetHelix(8);
4298 Float_t rc1 = TMath::Sqrt(xc1*xc1+yc1*yc1);
4299 Float_t fi1 = TMath::ATan2(yc1,xc1);
4301 Float_t dfi = fi0-fi1;
4304 if (dfi>1.5*TMath::Pi()) dfi-=TMath::Pi(); // take care about edge effect
4305 if (dfi<-1.5*TMath::Pi()) dfi+=TMath::Pi(); //
4306 if (TMath::Abs(dfi)>kMaxdPhi&&helixes[i0].GetHelix(4)*helixes[i1].GetHelix(4)<0){
4308 // if short tracks with undefined sign
4309 fi1 = -TMath::ATan2(yc1,-xc1);
4312 Float_t dtheta = TMath::Abs(track0->GetTgl()-track1->GetTgl())<TMath::Abs(track0->GetTgl()+track1->GetTgl())? track0->GetTgl()-track1->GetTgl():track0->GetTgl()+track1->GetTgl();
4315 // debug stream to tune "fast cuts"
4317 Double_t dist[3]; // distance at X
4318 Double_t mdist[3]={0,0,0}; // mean distance X+-40cm
4319 track0->GetDistance(track1,0.5*(xm[i0]+xm[i1])-40.,dist,AliTracker::GetBz());
4320 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[i]);
4321 track0->GetDistance(track1,0.5*(xm[i0]+xm[i1])+40.,dist,AliTracker::GetBz());
4322 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[i]);
4323 track0->GetDistance(track1,0.5*(xm[i0]+xm[i1]),dist,AliTracker::GetBz());
4324 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[i]);
4325 for (Int_t i=0;i<3;i++) mdist[i]*=0.33333;
4329 for (Int_t icl=0; icl<160; icl++){
4330 AliTPCclusterMI * cl0 = track0->GetClusterPointer(icl);
4331 AliTPCclusterMI * cl1 = track1->GetClusterPointer(icl);
4334 if (cl0==cl1) sums++;
4338 if (AliTPCReconstructor::StreamLevel()>5) {
4339 TTreeSRedirector &cstream = *fDebugStreamer;
4344 "Tr0.="<<track0<< // seed0
4345 "Tr1.="<<track1<< // seed1
4346 "h0.="<<&helixes[i0]<<
4347 "h1.="<<&helixes[i1]<<
4349 "sum="<<sum<< //the sum of rows with cl in both
4350 "sums="<<sums<< //the sum of shared clusters
4351 "xm0="<<xm[i0]<< // the center of track
4352 "xm1="<<xm[i1]<< // the x center of track
4353 // General cut variables
4354 "dfi="<<dfi<< // distance in fi angle
4355 "dtheta="<<dtheta<< // distance int theta angle
4361 "dist0="<<dist[0]<< //distance x
4362 "dist1="<<dist[1]<< //distance y
4363 "dist2="<<dist[2]<< //distance z
4364 "mdist0="<<mdist[0]<< //distance x
4365 "mdist1="<<mdist[1]<< //distance y
4366 "mdist2="<<mdist[2]<< //distance z
4382 if (AliTPCReconstructor::StreamLevel()>1) {
4383 AliInfo("Time for curling tracks removal DEBUGGING MC");
4390 void AliTPCtrackerMI::FindSplitted(TObjArray * array, AliESDEvent */*esd*/, Int_t /*iter*/){
4392 // Find Splitted tracks and remove the one with worst quality
4393 // Corresponding debug streamer to tune selections - "Splitted2"
4395 // 0. Sort tracks according quility
4396 // 1. Propagate the tracks to the reference radius
4397 // 2. Double_t loop to select close tracks (only to speed up process)
4398 // 3. Calculate cluster overlap ratio - and remove the track if bigger than a threshold
4399 // 4. Delete temporary parameters
4401 const Double_t xref=GetXrow(63); // reference radius -IROC/OROC boundary
4403 const Double_t kCutP1=10; // delta Z cut 10 cm
4404 const Double_t kCutP2=0.15; // delta snp(fi) cut 0.15
4405 const Double_t kCutP3=0.15; // delta tgl(theta) cut 0.15
4406 const Double_t kCutAlpha=0.15; // delta alpha cut
4407 Int_t firstpoint = 0;
4408 Int_t lastpoint = 160;
4410 Int_t nentries = array->GetEntriesFast();
4411 AliExternalTrackParam *params = new AliExternalTrackParam[nentries];
4417 //0. Sort tracks according quality
4418 //1. Propagate the ext. param to reference radius
4419 Int_t nseed = array->GetEntriesFast();
4420 if (nseed<=0) return;
4421 Float_t * quality = new Float_t[nseed];
4422 Int_t * indexes = new Int_t[nseed];
4423 for (Int_t i=0; i<nseed; i++) {
4424 AliTPCseed *pt=(AliTPCseed*)array->UncheckedAt(i);
4429 pt->UpdatePoints(); //select first last max dens points
4430 Float_t * points = pt->GetPoints();
4431 if (points[3]<0.8) quality[i] =-1;
4432 quality[i] = (points[2]-points[0])+pt->GetNumberOfClusters();
4433 //prefer high momenta tracks if overlaps
4434 quality[i] *= TMath::Sqrt(TMath::Abs(pt->Pt())+0.5);
4436 AliTracker::PropagateTrackToBxByBz(&(params[i]),xref,pt->GetMass(),5.,kTRUE);
4437 AliTracker::PropagateTrackToBxByBz(&(params[i]),xref,pt->GetMass(),1.,kTRUE);
4439 TMath::Sort(nseed,quality,indexes);
4441 // 3. Loop over pair of tracks
4443 for (Int_t i0=0; i0<nseed; i0++) {
4444 Int_t index0=indexes[i0];
4445 if (!(array->UncheckedAt(index0))) continue;
4446 AliTPCseed *s1 = (AliTPCseed*)array->UncheckedAt(index0);
4447 if (!s1->IsActive()) continue;
4448 AliExternalTrackParam &par0=params[index0];
4449 for (Int_t i1=i0+1; i1<nseed; i1++) {
4450 Int_t index1=indexes[i1];
4451 if (!(array->UncheckedAt(index1))) continue;
4452 AliTPCseed *s2 = (AliTPCseed*)array->UncheckedAt(index1);
4453 if (!s2->IsActive()) continue;
4454 if (s2->GetKinkIndexes()[0]!=0)
4455 if (s2->GetKinkIndexes()[0] == -s1->GetKinkIndexes()[0]) continue;
4456 AliExternalTrackParam &par1=params[index1];
4457 if (TMath::Abs(par0.GetParameter()[3]-par1.GetParameter()[3])>kCutP3) continue;
4458 if (TMath::Abs(par0.GetParameter()[1]-par1.GetParameter()[1])>kCutP1) continue;
4459 if (TMath::Abs(par0.GetParameter()[2]-par1.GetParameter()[2])>kCutP2) continue;
4460 Double_t dAlpha= TMath::Abs(par0.GetAlpha()-par1.GetAlpha());
4461 if (dAlpha>TMath::Pi()) dAlpha-=TMath::Pi();
4462 if (TMath::Abs(dAlpha)>kCutAlpha) continue;
4467 Int_t firstShared=lastpoint, lastShared=firstpoint;
4468 Int_t firstRow=lastpoint, lastRow=firstpoint;
4470 for (Int_t i=firstpoint;i<lastpoint;i++){
4471 if (s1->GetClusterIndex2(i)>0) nall0++;
4472 if (s2->GetClusterIndex2(i)>0) nall1++;
4473 if (s1->GetClusterIndex2(i)>0 && s2->GetClusterIndex2(i)>0) {
4474 if (i<firstRow) firstRow=i;
4475 if (i>lastRow) lastRow=i;
4477 if ( (s1->GetClusterIndex2(i))==(s2->GetClusterIndex2(i)) && s1->GetClusterIndex2(i)>0) {
4478 if (i<firstShared) firstShared=i;
4479 if (i>lastShared) lastShared=i;
4483 Double_t ratio0 = Float_t(sumShared)/Float_t(TMath::Min(nall0+1,nall1+1));
4484 Double_t ratio1 = Float_t(sumShared)/Float_t(TMath::Max(nall0+1,nall1+1));
4486 if( AliTPCReconstructor::StreamLevel()>1){
4487 TTreeSRedirector &cstream = *fDebugStreamer;
4488 Int_t n0=s1->GetNumberOfClusters();
4489 Int_t n1=s2->GetNumberOfClusters();
4490 Int_t n0F=s1->GetNFoundable();
4491 Int_t n1F=s2->GetNFoundable();
4492 Int_t lab0=s1->GetLabel();
4493 Int_t lab1=s2->GetLabel();
4495 cstream<<"Splitted2"<<
4496 "iter="<<fIteration<<
4497 "lab0="<<lab0<< // MC label if exist
4498 "lab1="<<lab1<< // MC label if exist
4501 "ratio0="<<ratio0<< // shared ratio
4502 "ratio1="<<ratio1<< // shared ratio
4503 "p0.="<<&par0<< // track parameters
4505 "s0.="<<s1<< // full seed
4507 "n0="<<n0<< // number of clusters track 0
4508 "n1="<<n1<< // number of clusters track 1
4509 "nall0="<<nall0<< // number of clusters track 0
4510 "nall1="<<nall1<< // number of clusters track 1
4511 "n0F="<<n0F<< // number of findable
4512 "n1F="<<n1F<< // number of findable
4513 "shared="<<sumShared<< // number of shared clusters
4514 "firstS="<<firstShared<< // first and the last shared row
4515 "lastS="<<lastShared<<
4516 "firstRow="<<firstRow<< // first and the last row with cluster
4517 "lastRow="<<lastRow<< //
4521 // remove track with lower quality
4523 if (ratio0>AliTPCReconstructor::GetRecoParam()->GetCutSharedClusters(0) ||
4524 ratio1>AliTPCReconstructor::GetRecoParam()->GetCutSharedClusters(1)){
4528 delete array->RemoveAt(index1);
4533 // 4. Delete temporary array
4543 void AliTPCtrackerMI::FindCurling(const TObjArray * array, AliESDEvent */*esd*/, Int_t iter)
4546 // find Curling tracks
4547 // Use AliTPCReconstructor::StreamLevel()>1 if you want to tune parameters - cuts
4550 // Algorithm done in 2 phases - because of CPU consumption
4551 // it is n^2 algorithm - for lead-lead 20000x20000 combination are investigated
4552 // see detal in MC part what can be used to cut
4556 const Float_t kMaxC = 400; // maximal curvature to of the track
4557 const Float_t kMaxdTheta = 0.15; // maximal distance in theta
4558 const Float_t kMaxdPhi = 0.15; // maximal distance in phi
4559 const Float_t kPtRatio = 0.3; // ratio between pt
4560 const Float_t kMinDCAR = 2.; // distance to the primary vertex in r - see cpipe cut
4563 // Curling tracks cuts
4566 const Float_t kMaxDeltaRMax = 40; // distance in outer radius
4567 const Float_t kMaxDeltaRMin = 5.; // distance in lower radius - see cpipe cut
4568 const Float_t kMinAngle = 2.9; // angle between tracks
4569 const Float_t kMaxDist = 5; // biggest distance
4571 // The cuts can be tuned using the "MC information stored in Multi tree ==> see FindMultiMC
4574 TCut csign("csign","Tr0.fP[4]*Tr1.fP[4]<0"); //opposite sign
4575 TCut cmax("cmax","abs(Tr0.GetC())>1/400");
4576 TCut cda("cda","sqrt(dtheta^2+dfi^2)<0.15");
4577 TCut ccratio("ccratio","abs((Tr0.fP[4]+Tr1.fP[4])/(abs(Tr0.fP[4])+abs(Tr1.fP[4])))<0.3");
4578 TCut cpipe("cpipe", "min(abs(r0-rc0),abs(r1-rc1))>5");
4580 TCut cdrmax("cdrmax","abs(abs(rc0+r0)-abs(rc1+r1))<40")
4581 TCut cdrmin("cdrmin","abs(abs(rc0+r0)-abs(rc1+r1))<10")
4583 Multi->Draw("dfi","iter==0"+csign+cmax+cda+ccratio); ~94% of curling tracks fulfill
4584 Multi->Draw("min(abs(r0-rc0),abs(r1-rc1))","iter==0&&abs(lab1)==abs(lab0)"+csign+cmax+cda+ccratio+cpipe+cdrmin+cdrmax); //80%
4586 Curling2->Draw("dfi","iter==0&&abs(lab0)==abs(lab1)"+csign+cmax+cdtheta+cdfi+ccratio)
4592 Int_t nentries = array->GetEntriesFast();
4593 AliHelix *helixes = new AliHelix[nentries];
4594 for (Int_t i=0;i<nentries;i++){
4595 AliTPCseed* track = (AliTPCseed*)array->At(i);
4596 if (!track) continue;
4597 track->SetCircular(0);
4598 new (&helixes[i]) AliHelix(*track);
4604 Double_t phase[2][2]={{0,0},{0,0}},radius[2]={0,0};
4610 for (Int_t i0=0;i0<nentries;i0++){
4611 AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
4612 if (!track0) continue;
4613 if (TMath::Abs(track0->GetC())<1/kMaxC) continue;
4614 Float_t xc0 = helixes[i0].GetHelix(6);
4615 Float_t yc0 = helixes[i0].GetHelix(7);
4616 Float_t r0 = helixes[i0].GetHelix(8);
4617 Float_t rc0 = TMath::Sqrt(xc0*xc0+yc0*yc0);
4618 Float_t fi0 = TMath::ATan2(yc0,xc0);
4620 for (Int_t i1=i0+1;i1<nentries;i1++){
4621 AliTPCseed * track1 = (AliTPCseed*)array->At(i1);
4622 if (!track1) continue;
4623 if (TMath::Abs(track1->GetC())<1/kMaxC) continue;
4624 Float_t xc1 = helixes[i1].GetHelix(6);
4625 Float_t yc1 = helixes[i1].GetHelix(7);
4626 Float_t r1 = helixes[i1].GetHelix(8);
4627 Float_t rc1 = TMath::Sqrt(xc1*xc1+yc1*yc1);
4628 Float_t fi1 = TMath::ATan2(yc1,xc1);
4630 Float_t dfi = fi0-fi1;
4633 if (dfi>1.5*TMath::Pi()) dfi-=TMath::Pi(); // take care about edge effect
4634 if (dfi<-1.5*TMath::Pi()) dfi+=TMath::Pi(); //
4635 Float_t dtheta = TMath::Abs(track0->GetTgl()-track1->GetTgl())<TMath::Abs(track0->GetTgl()+track1->GetTgl())? track0->GetTgl()-track1->GetTgl():track0->GetTgl()+track1->GetTgl();
4639 if (track0->GetBConstrain()&&track1->GetBConstrain()) continue; // not constrained
4640 if (track1->GetSigned1Pt()*track0->GetSigned1Pt()>0) continue; // not the same sign
4641 if ( TMath::Abs(track1->GetTgl()+track0->GetTgl())>kMaxdTheta) continue; //distance in the Theta
4642 if ( TMath::Abs(dfi)>kMaxdPhi) continue; //distance in phi
4643 if ( TMath::Sqrt(dfi*dfi+dtheta*dtheta)>kMaxdPhi) continue; //common angular offset
4645 Float_t pt0 = track0->GetSignedPt();
4646 Float_t pt1 = track1->GetSignedPt();
4647 if ((TMath::Abs(pt0+pt1)/(TMath::Abs(pt0)+TMath::Abs(pt1)))>kPtRatio) continue;
4648 if ((iter==1) && TMath::Abs(TMath::Abs(rc0+r0)-TMath::Abs(rc1+r1))>kMaxDeltaRMax) continue;
4649 if ((iter!=1) &&TMath::Abs(TMath::Abs(rc0-r0)-TMath::Abs(rc1-r1))>kMaxDeltaRMin) continue;
4650 if (TMath::Min(TMath::Abs(rc0-r0),TMath::Abs(rc1-r1))<kMinDCAR) continue;
4653 // Now find closest approach
4657 Int_t npoints = helixes[i0].GetRPHIintersections(helixes[i1], phase, radius,10);
4658 if (npoints==0) continue;
4659 helixes[i0].GetClosestPhases(helixes[i1], phase);
4663 Double_t hangles[3];
4664 helixes[i0].Evaluate(phase[0][0],xyz0);
4665 helixes[i1].Evaluate(phase[0][1],xyz1);
4667 helixes[i0].GetAngle(phase[0][0],helixes[i1],phase[0][1],hangles);
4668 Double_t deltah[2],deltabest;
4669 if (TMath::Abs(hangles[2])<kMinAngle) continue;
4673 helixes[i0].ParabolicDCA(helixes[i1],phase[0][0],phase[0][1],radius[0],deltah[0],2);
4675 helixes[i0].ParabolicDCA(helixes[i1],phase[1][0],phase[1][1],radius[1],deltah[1],2);
4676 if (deltah[1]<deltah[0]) ibest=1;
4678 deltabest = TMath::Sqrt(deltah[ibest]);
4679 helixes[i0].Evaluate(phase[ibest][0],xyz0);
4680 helixes[i1].Evaluate(phase[ibest][1],xyz1);
4681 helixes[i0].GetAngle(phase[ibest][0],helixes[i1],phase[ibest][1],hangles);
4682 Double_t radiusbest = TMath::Sqrt(radius[ibest]);
4684 if (deltabest>kMaxDist) continue;
4685 // if (mindcar+mindcaz<40 && (TMath::Abs(hangles[2])<kMinAngle ||deltabest>3)) continue;
4686 Bool_t sign =kFALSE;
4687 if (hangles[2]>kMinAngle) sign =kTRUE;
4690 // circular[i0] = kTRUE;
4691 // circular[i1] = kTRUE;
4692 if (track0->OneOverPt()<track1->OneOverPt()){
4693 track0->SetCircular(track0->GetCircular()+1);
4694 track1->SetCircular(track1->GetCircular()+2);
4697 track1->SetCircular(track1->GetCircular()+1);
4698 track0->SetCircular(track0->GetCircular()+2);
4701 if (AliTPCReconstructor::StreamLevel()>2){
4703 //debug stream to tune "fine" cuts
4704 Int_t lab0=track0->GetLabel();
4705 Int_t lab1=track1->GetLabel();
4706 TTreeSRedirector &cstream = *fDebugStreamer;
4707 cstream<<"Curling2"<<
4723 "npoints="<<npoints<<
4724 "hangles0="<<hangles[0]<<
4725 "hangles1="<<hangles[1]<<
4726 "hangles2="<<hangles[2]<<
4729 "radius="<<radiusbest<<
4730 "deltabest="<<deltabest<<
4731 "phase0="<<phase[ibest][0]<<
4732 "phase1="<<phase[ibest][1]<<
4740 if (AliTPCReconstructor::StreamLevel()>1) {
4741 AliInfo("Time for curling tracks removal");
4750 void AliTPCtrackerMI::FindKinks(TObjArray * array, AliESDEvent *esd)
4757 TObjArray *kinks= new TObjArray(10000);
4758 // TObjArray *v0s= new TObjArray(10000);
4759 Int_t nentries = array->GetEntriesFast();
4760 AliHelix *helixes = new AliHelix[nentries];
4761 Int_t *sign = new Int_t[nentries];
4762 Int_t *nclusters = new Int_t[nentries];
4763 Float_t *alpha = new Float_t[nentries];
4764 AliKink *kink = new AliKink();
4765 Int_t * usage = new Int_t[nentries];
4766 Float_t *zm = new Float_t[nentries];
4767 Float_t *z0 = new Float_t[nentries];
4768 Float_t *fim = new Float_t[nentries];
4769 Float_t *shared = new Float_t[nentries];
4770 Bool_t *circular = new Bool_t[nentries];
4771 Float_t *dca = new Float_t[nentries];
4772 //const AliESDVertex * primvertex = esd->GetVertex();
4774 // nentries = array->GetEntriesFast();
4779 for (Int_t i=0;i<nentries;i++){
4782 AliTPCseed* track = (AliTPCseed*)array->At(i);
4783 if (!track) continue;
4784 track->SetCircular(0);
4786 track->UpdatePoints();
4787 if (( track->GetPoints()[2]- track->GetPoints()[0])>5 && track->GetPoints()[3]>0.8){
4789 nclusters[i]=track->GetNumberOfClusters();
4790 alpha[i] = track->GetAlpha();
4791 new (&helixes[i]) AliHelix(*track);
4793 helixes[i].Evaluate(0,xyz);
4794 sign[i] = (track->GetC()>0) ? -1:1;
4797 if (track->GetProlongation(x,y,z)){
4799 fim[i] = alpha[i]+TMath::ATan2(y,x);
4802 zm[i] = track->GetZ();
4806 circular[i]= kFALSE;
4807 if (track->GetProlongation(0,y,z)) z0[i] = z;
4808 dca[i] = track->GetD(0,0);
4814 Int_t ncandidates =0;
4817 Double_t phase[2][2]={{0,0},{0,0}},radius[2]={0,0};
4820 // Find circling track
4822 for (Int_t i0=0;i0<nentries;i0++){
4823 AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
4824 if (!track0) continue;
4825 if (track0->GetNumberOfClusters()<40) continue;
4826 if (TMath::Abs(1./track0->GetC())>200) continue;
4827 for (Int_t i1=i0+1;i1<nentries;i1++){
4828 AliTPCseed * track1 = (AliTPCseed*)array->At(i1);
4829 if (!track1) continue;
4830 if (track1->GetNumberOfClusters()<40) continue;
4831 if ( TMath::Abs(track1->GetTgl()+track0->GetTgl())>0.1) continue;
4832 if (track0->GetBConstrain()&&track1->GetBConstrain()) continue;
4833 if (TMath::Abs(1./track1->GetC())>200) continue;
4834 if (track1->GetSigned1Pt()*track0->GetSigned1Pt()>0) continue;
4835 if (track1->GetTgl()*track0->GetTgl()>0) continue;
4836 if (TMath::Max(TMath::Abs(1./track0->GetC()),TMath::Abs(1./track1->GetC()))>190) continue;
4837 if (track0->GetBConstrain()&&track1->OneOverPt()<track0->OneOverPt()) continue; //returning - lower momenta
4838 if (track1->GetBConstrain()&&track0->OneOverPt()<track1->OneOverPt()) continue; //returning - lower momenta
4840 Float_t mindcar = TMath::Min(TMath::Abs(dca[i0]),TMath::Abs(dca[i1]));
4841 if (mindcar<5) continue;
4842 Float_t mindcaz = TMath::Min(TMath::Abs(z0[i0]-GetZ()),TMath::Abs(z0[i1]-GetZ()));
4843 if (mindcaz<5) continue;
4844 if (mindcar+mindcaz<20) continue;
4847 Float_t xc0 = helixes[i0].GetHelix(6);
4848 Float_t yc0 = helixes[i0].GetHelix(7);
4849 Float_t r0 = helixes[i0].GetHelix(8);
4850 Float_t xc1 = helixes[i1].GetHelix(6);
4851 Float_t yc1 = helixes[i1].GetHelix(7);
4852 Float_t r1 = helixes[i1].GetHelix(8);
4854 Float_t rmean = (r0+r1)*0.5;
4855 Float_t delta =TMath::Sqrt((xc1-xc0)*(xc1-xc0)+(yc1-yc0)*(yc1-yc0));
4856 //if (delta>30) continue;
4857 if (delta>rmean*0.25) continue;
4858 if (TMath::Abs(r0-r1)/rmean>0.3) continue;
4860 Int_t npoints = helixes[i0].GetRPHIintersections(helixes[i1], phase, radius,10);
4861 if (npoints==0) continue;
4862 helixes[i0].GetClosestPhases(helixes[i1], phase);
4866 Double_t hangles[3];
4867 helixes[i0].Evaluate(phase[0][0],xyz0);
4868 helixes[i1].Evaluate(phase[0][1],xyz1);
4870 helixes[i0].GetAngle(phase[0][0],helixes[i1],phase[0][1],hangles);
4871 Double_t deltah[2],deltabest;
4872 if (hangles[2]<2.8) continue;
4875 helixes[i0].ParabolicDCA(helixes[i1],phase[0][0],phase[0][1],radius[0],deltah[0],2);
4877 helixes[i0].ParabolicDCA(helixes[i1],phase[1][0],phase[1][1],radius[1],deltah[1],2);
4878 if (deltah[1]<deltah[0]) ibest=1;
4880 deltabest = TMath::Sqrt(deltah[ibest]);
4881 helixes[i0].Evaluate(phase[ibest][0],xyz0);
4882 helixes[i1].Evaluate(phase[ibest][1],xyz1);
4883 helixes[i0].GetAngle(phase[ibest][0],helixes[i1],phase[ibest][1],hangles);
4884 Double_t radiusbest = TMath::Sqrt(radius[ibest]);
4886 if (deltabest>6) continue;
4887 if (mindcar+mindcaz<40 && (hangles[2]<3.12||deltabest>3)) continue;
4888 Bool_t lsign =kFALSE;
4889 if (hangles[2]>3.06) lsign =kTRUE;
4892 circular[i0] = kTRUE;
4893 circular[i1] = kTRUE;
4894 if (track0->OneOverPt()<track1->OneOverPt()){
4895 track0->SetCircular(track0->GetCircular()+1);
4896 track1->SetCircular(track1->GetCircular()+2);
4899 track1->SetCircular(track1->GetCircular()+1);
4900 track0->SetCircular(track0->GetCircular()+2);
4903 if (lsign&&AliTPCReconstructor::StreamLevel()>1){
4905 Int_t lab0=track0->GetLabel();
4906 Int_t lab1=track1->GetLabel();
4907 TTreeSRedirector &cstream = *fDebugStreamer;
4908 cstream<<"Curling"<<
4915 "mindcar="<<mindcar<<
4916 "mindcaz="<<mindcaz<<
4919 "npoints="<<npoints<<
4920 "hangles0="<<hangles[0]<<
4921 "hangles2="<<hangles[2]<<
4926 "radius="<<radiusbest<<
4927 "deltabest="<<deltabest<<
4928 "phase0="<<phase[ibest][0]<<
4929 "phase1="<<phase[ibest][1]<<
4939 for (Int_t i =0;i<nentries;i++){
4940 if (sign[i]==0) continue;
4941 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
4948 Double_t cradius0 = 40*40;
4949 Double_t cradius1 = 270*270;
4952 Double_t cdist3=0.55;
4953 for (Int_t j =i+1;j<nentries;j++){
4955 if (sign[j]*sign[i]<1) continue;
4956 if ( (nclusters[i]+nclusters[j])>200) continue;
4957 if ( (nclusters[i]+nclusters[j])<80) continue;
4958 if ( TMath::Abs(zm[i]-zm[j])>60.) continue;
4959 if ( TMath::Abs(fim[i]-fim[j])>0.6 && TMath::Abs(fim[i]-fim[j])<5.7 ) continue;
4960 //AliTPCseed * track1 = (AliTPCseed*)array->At(j); Double_t phase[2][2],radius[2];
4961 Int_t npoints = helixes[i].GetRPHIintersections(helixes[j], phase, radius,20);
4962 if (npoints<1) continue;
4965 if (radius[0]<cradius0||radius[0]>cradius1) continue;
4968 if ( (radius[0]<cradius0||radius[0]>cradius1) && (radius[1]<cradius0||radius[1]>cradius1) ) continue;
4971 Double_t delta1=10000,delta2=10000;
4972 // cuts on the intersection radius
4973 helixes[i].LinearDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta1);
4974 if (radius[0]<20&&delta1<1) continue; //intersection at vertex
4975 if (radius[0]<10&&delta1<3) continue; //intersection at vertex
4977 helixes[i].LinearDCA(helixes[j],phase[1][0],phase[1][1],radius[1],delta2);
4978 if (radius[1]<20&&delta2<1) continue; //intersection at vertex
4979 if (radius[1]<10&&delta2<3) continue; //intersection at vertex
4982 Double_t distance1 = TMath::Min(delta1,delta2);
4983 if (distance1>cdist1) continue; // cut on DCA linear approximation
4985 npoints = helixes[i].GetRPHIintersections(helixes[j], phase, radius,20);
4986 helixes[i].ParabolicDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta1);
4987 if (radius[0]<20&&delta1<1) continue; //intersection at vertex
4988 if (radius[0]<10&&delta1<3) continue; //intersection at vertex
4991 helixes[i].ParabolicDCA(helixes[j],phase[1][0],phase[1][1],radius[1],delta2);
4992 if (radius[1]<20&&delta2<1) continue; //intersection at vertex
4993 if (radius[1]<10&&delta2<3) continue; //intersection at vertex
4995 distance1 = TMath::Min(delta1,delta2);
4998 rkink = TMath::Sqrt(radius[0]);
5001 rkink = TMath::Sqrt(radius[1]);
5003 if (distance1>cdist2) continue;
5006 AliTPCseed * track1 = (AliTPCseed*)array->At(j);
5009 Int_t row0 = GetRowNumber(rkink);
5010 if (row0<10) continue;
5011 if (row0>150) continue;
5014 Float_t dens00=-1,dens01=-1;
5015 Float_t dens10=-1,dens11=-1;
5017 Int_t found,foundable,ishared;
5018 track0->GetClusterStatistic(0,row0-5, found, foundable,ishared,kFALSE);
5019 if (foundable>5) dens00 = Float_t(found)/Float_t(foundable);
5020 track0->GetClusterStatistic(row0+5,155, found, foundable,ishared,kFALSE);
5021 if (foundable>5) dens01 = Float_t(found)/Float_t(foundable);
5023 track1->GetClusterStatistic(0,row0-5, found, foundable,ishared,kFALSE);
5024 if (foundable>10) dens10 = Float_t(found)/Float_t(foundable);
5025 track1->GetClusterStatistic(row0+5,155, found, foundable,ishared,kFALSE);
5026 if (foundable>10) dens11 = Float_t(found)/Float_t(foundable);
5028 if (dens00<dens10 && dens01<dens11) continue;
5029 if (dens00>dens10 && dens01>dens11) continue;
5030 if (TMath::Max(dens00,dens10)<0.1) continue;
5031 if (TMath::Max(dens01,dens11)<0.3) continue;
5033 if (TMath::Min(dens00,dens10)>0.6) continue;
5034 if (TMath::Min(dens01,dens11)>0.6) continue;
5037 AliTPCseed * ktrack0, *ktrack1;
5046 if (TMath::Abs(ktrack0->GetC())>5) continue; // cut on the curvature for mother particle
5047 AliExternalTrackParam paramm(*ktrack0);
5048 AliExternalTrackParam paramd(*ktrack1);
5049 if (row0>60&&ktrack1->GetReference().GetX()>90.)new (¶md) AliExternalTrackParam(ktrack1->GetReference());
5052 kink->SetMother(paramm);
5053 kink->SetDaughter(paramd);
5056 Float_t x[3] = { kink->GetPosition()[0],kink->GetPosition()[1],kink->GetPosition()[2]};
5058 fkParam->Transform0to1(x,index);
5059 fkParam->Transform1to2(x,index);
5060 row0 = GetRowNumber(x[0]);
5062 if (kink->GetR()<100) continue;
5063 if (kink->GetR()>240) continue;
5064 if (kink->GetPosition()[2]/kink->GetR()>AliTPCReconstructor::GetCtgRange()) continue; //out of fiducial volume
5065 if (kink->GetDistance()>cdist3) continue;
5066 Float_t dird = kink->GetDaughterP()[0]*kink->GetPosition()[0]+kink->GetDaughterP()[1]*kink->GetPosition()[1]; // rough direction estimate
5067 if (dird<0) continue;
5069 Float_t dirm = kink->GetMotherP()[0]*kink->GetPosition()[0]+kink->GetMotherP()[1]*kink->GetPosition()[1]; // rough direction estimate
5070 if (dirm<0) continue;
5071 Float_t mpt = TMath::Sqrt(kink->GetMotherP()[0]*kink->GetMotherP()[0]+kink->GetMotherP()[1]*kink->GetMotherP()[1]);
5072 if (mpt<0.2) continue;
5075 //for high momenta momentum not defined well in first iteration
5076 Double_t qt = TMath::Sin(kink->GetAngle(2))*ktrack1->GetP();
5077 if (qt>0.35) continue;
5080 kink->SetLabel(CookLabel(ktrack0,0.4,0,row0),0);
5081 kink->SetLabel(CookLabel(ktrack1,0.4,row0,160),1);
5083 kink->SetTPCDensity(dens00,0,0);
5084 kink->SetTPCDensity(dens01,0,1);
5085 kink->SetTPCDensity(dens10,1,0);
5086 kink->SetTPCDensity(dens11,1,1);
5087 kink->SetIndex(i,0);
5088 kink->SetIndex(j,1);
5091 kink->SetTPCDensity(dens10,0,0);
5092 kink->SetTPCDensity(dens11,0,1);
5093 kink->SetTPCDensity(dens00,1,0);
5094 kink->SetTPCDensity(dens01,1,1);
5095 kink->SetIndex(j,0);
5096 kink->SetIndex(i,1);
5099 if (mpt<1||kink->GetAngle(2)>0.1){
5100 // angle and densities not defined yet
5101 if (kink->GetTPCDensityFactor()<0.8) continue;
5102 if ((2-kink->GetTPCDensityFactor())*kink->GetDistance() >0.25) continue;
5103 if (kink->GetAngle(2)*ktrack0->GetP()<0.003) continue; //too small angle
5104 if (kink->GetAngle(2)>0.2&&kink->GetTPCDensityFactor()<1.15) continue;
5105 if (kink->GetAngle(2)>0.2&&kink->GetTPCDensity(0,1)>0.05) continue;
5107 Float_t criticalangle = track0->GetSigmaSnp2()+track0->GetSigmaTgl2();
5108 criticalangle+= track1->GetSigmaSnp2()+track1->GetSigmaTgl2();
5109 criticalangle= 3*TMath::Sqrt(criticalangle);
5110 if (criticalangle>0.02) criticalangle=0.02;
5111 if (kink->GetAngle(2)<criticalangle) continue;
5114 Int_t drow = Int_t(2.+0.5/(0.05+kink->GetAngle(2))); // overlap region defined
5115 Float_t shapesum =0;
5117 for ( Int_t row = row0-drow; row<row0+drow;row++){
5118 if (row<0) continue;
5119 if (row>155) continue;
5120 if (ktrack0->GetClusterPointer(row)){
5121 AliTPCTrackerPoint *point =ktrack0->GetTrackPoint(row);
5122 shapesum+=point->GetSigmaY()+point->GetSigmaZ();
5125 if (ktrack1->GetClusterPointer(row)){
5126 AliTPCTrackerPoint *point =ktrack1->GetTrackPoint(row);
5127 shapesum+=point->GetSigmaY()+point->GetSigmaZ();
5132 kink->SetShapeFactor(-1.);
5135 kink->SetShapeFactor(shapesum/sum);
5137 // esd->AddKink(kink);
5139 // kink->SetMother(paramm);
5140 //kink->SetDaughter(paramd);
5142 Double_t chi2P2 = paramm.GetParameter()[2]-paramd.GetParameter()[2];
5144 chi2P2/=paramm.GetCovariance()[5]+paramd.GetCovariance()[5];
5145 Double_t chi2P3 = paramm.GetParameter()[3]-paramd.GetParameter()[3];
5147 chi2P3/=paramm.GetCovariance()[9]+paramd.GetCovariance()[9];
5149 if (AliTPCReconstructor::StreamLevel()>1) {
5150 (*fDebugStreamer)<<"kinkLpt"<<
5158 if ( chi2P2+chi2P3<AliTPCReconstructor::GetRecoParam()->GetKinkAngleCutChi2(0)){
5162 kinks->AddLast(kink);
5168 // sort the kinks according quality - and refit them towards vertex
5170 Int_t nkinks = kinks->GetEntriesFast();
5171 Float_t *quality = new Float_t[nkinks];
5172 Int_t *indexes = new Int_t[nkinks];
5173 AliTPCseed *mothers = new AliTPCseed[nkinks];
5174 AliTPCseed *daughters = new AliTPCseed[nkinks];
5177 for (Int_t i=0;i<nkinks;i++){
5179 AliKink *kinkl = (AliKink*)kinks->At(i);
5181 // refit kinks towards vertex
5183 Int_t index0 = kinkl->GetIndex(0);
5184 Int_t index1 = kinkl->GetIndex(1);
5185 AliTPCseed * ktrack0 = (AliTPCseed*)array->At(index0);
5186 AliTPCseed * ktrack1 = (AliTPCseed*)array->At(index1);
5188 Int_t sumn=ktrack0->GetNumberOfClusters()+ktrack1->GetNumberOfClusters();
5190 // Refit Kink under if too small angle
5192 if (kinkl->GetAngle(2)<0.05){
5193 kinkl->SetTPCRow0(GetRowNumber(kinkl->GetR()));
5194 Int_t row0 = kinkl->GetTPCRow0();
5195 Int_t drow = Int_t(2.+0.5/(0.05+kinkl->GetAngle(2)));
5198 Int_t last = row0-drow;
5199 if (last<40) last=40;
5200 if (last<ktrack0->GetFirstPoint()+25) last = ktrack0->GetFirstPoint()+25;
5201 AliTPCseed* seed0 = ReSeed(ktrack0,last,kFALSE);
5204 Int_t first = row0+drow;
5205 if (first>130) first=130;
5206 if (first>ktrack1->GetLastPoint()-25) first = TMath::Max(ktrack1->GetLastPoint()-25,30);
5207 AliTPCseed* seed1 = ReSeed(ktrack1,first,kTRUE);
5209 if (seed0 && seed1){
5210 kinkl->SetStatus(1,8);
5211 if (RefitKink(*seed0,*seed1,*kinkl)) kinkl->SetStatus(1,9);
5212 row0 = GetRowNumber(kinkl->GetR());
5213 sumn = seed0->GetNumberOfClusters()+seed1->GetNumberOfClusters();
5214 mothers[i] = *seed0;
5215 daughters[i] = *seed1;
5218 delete kinks->RemoveAt(i);
5219 if (seed0) delete seed0;
5220 if (seed1) delete seed1;
5223 if (kinkl->GetDistance()>0.5 || kinkl->GetR()<110 || kinkl->GetR()>240) {
5224 delete kinks->RemoveAt(i);
5225 if (seed0) delete seed0;
5226 if (seed1) delete seed1;
5234 if (kinkl) quality[i] = 160*((0.1+kinkl->GetDistance())*(2.-kinkl->GetTPCDensityFactor()))/(sumn+40.); //the longest -clossest will win
5236 TMath::Sort(nkinks,quality,indexes,kFALSE);
5238 //remove double find kinks
5240 for (Int_t ikink0=1;ikink0<nkinks;ikink0++){
5241 AliKink * kink0 = (AliKink*) kinks->At(indexes[ikink0]);
5242 if (!kink0) continue;
5244 for (Int_t ikink1=0;ikink1<ikink0;ikink1++){
5245 kink0 = (AliKink*) kinks->At(indexes[ikink0]);
5246 if (!kink0) continue;
5247 AliKink * kink1 = (AliKink*) kinks->At(indexes[ikink1]);
5248 if (!kink1) continue;
5249 // if not close kink continue
5250 if (TMath::Abs(kink1->GetPosition()[2]-kink0->GetPosition()[2])>10) continue;
5251 if (TMath::Abs(kink1->GetPosition()[1]-kink0->GetPosition()[1])>10) continue;
5252 if (TMath::Abs(kink1->GetPosition()[0]-kink0->GetPosition()[0])>10) continue;
5254 AliTPCseed &mother0 = mothers[indexes[ikink0]];
5255 AliTPCseed &daughter0 = daughters[indexes[ikink0]];
5256 AliTPCseed &mother1 = mothers[indexes[ikink1]];
5257 AliTPCseed &daughter1 = daughters[indexes[ikink1]];
5258 Int_t row0 = (kink0->GetTPCRow0()+kink1->GetTPCRow0())/2;
5267 for (Int_t i=0;i<row0;i++){
5268 if (mother0.GetClusterIndex(i)>0 && mother1.GetClusterIndex(i)>0){
5271 if (mother0.GetClusterIndex(i)==mother1.GetClusterIndex(i)){
5278 for (Int_t i=row0;i<158;i++){
5279 if (daughter0.GetClusterIndex(i)>0 && daughter0.GetClusterIndex(i)>0){
5282 if (mother0.GetClusterIndex(i)==mother1.GetClusterIndex(i)){
5288 Float_t ratio = Float_t(same+1)/Float_t(both+1);
5289 Float_t ratiom = Float_t(samem+1)/Float_t(bothm+1);
5290 Float_t ratiod = Float_t(samed+1)/Float_t(bothd+1);
5291 if (ratio>0.3 && ratiom>0.5 &&ratiod>0.5) {
5292 Int_t sum0 = mother0.GetNumberOfClusters()+daughter0.GetNumberOfClusters();
5293 Int_t sum1 = mother1.GetNumberOfClusters()+daughter1.GetNumberOfClusters();
5295 shared[kink0->GetIndex(0)]= kTRUE;
5296 shared[kink0->GetIndex(1)]= kTRUE;
5297 delete kinks->RemoveAt(indexes[ikink0]);
5301 shared[kink1->GetIndex(0)]= kTRUE;
5302 shared[kink1->GetIndex(1)]= kTRUE;
5303 delete kinks->RemoveAt(indexes[ikink1]);
5310 for (Int_t i=0;i<nkinks;i++){
5311 AliKink * kinkl = (AliKink*) kinks->At(indexes[i]);
5312 if (!kinkl) continue;
5313 kinkl->SetTPCRow0(GetRowNumber(kinkl->GetR()));
5314 Int_t index0 = kinkl->GetIndex(0);
5315 Int_t index1 = kinkl->GetIndex(1);
5316 if (circular[index0]||(circular[index1]&&kinkl->GetDistance()>0.2)) continue;
5317 kinkl->SetMultiple(usage[index0],0);
5318 kinkl->SetMultiple(usage[index1],1);
5319 if (kinkl->GetMultiple()[0]+kinkl->GetMultiple()[1]>2) continue;
5320 if (kinkl->GetMultiple()[0]+kinkl->GetMultiple()[1]>0 && quality[indexes[i]]>0.2) continue;
5321 if (kinkl->GetMultiple()[0]+kinkl->GetMultiple()[1]>0 && kinkl->GetDistance()>0.2) continue;
5322 if (circular[index0]||(circular[index1]&&kinkl->GetDistance()>0.1)) continue;
5324 AliTPCseed * ktrack0 = (AliTPCseed*)array->At(index0);
5325 AliTPCseed * ktrack1 = (AliTPCseed*)array->At(index1);
5326 if (!ktrack0 || !ktrack1) continue;
5327 Int_t index = esd->AddKink(kinkl);
5330 if ( ktrack0->GetKinkIndex(0)==0 && ktrack1->GetKinkIndex(0)==0) { //best kink
5331 if (mothers[indexes[i]].GetNumberOfClusters()>20 && daughters[indexes[i]].GetNumberOfClusters()>20 && (mothers[indexes[i]].GetNumberOfClusters()+daughters[indexes[i]].GetNumberOfClusters())>100){
5332 *ktrack0 = mothers[indexes[i]];
5333 *ktrack1 = daughters[indexes[i]];
5337 ktrack0->SetKinkIndex(usage[index0],-(index+1));
5338 ktrack1->SetKinkIndex(usage[index1], (index+1));
5343 // Remove tracks corresponding to shared kink's
5345 for (Int_t i=0;i<nentries;i++){
5346 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5347 if (!track0) continue;
5348 if (track0->GetKinkIndex(0)!=0) continue;
5349 if (shared[i]) delete array->RemoveAt(i);
5354 RemoveUsed2(array,0.5,0.4,30);
5356 for (Int_t i=0;i<nentries;i++){
5357 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5358 if (!track0) continue;
5359 track0->CookdEdx(0.02,0.6);
5363 for (Int_t i=0;i<nentries;i++){
5364 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5365 if (!track0) continue;
5366 if (track0->Pt()<1.4) continue;
5367 //remove double high momenta tracks - overlapped with kink candidates
5370 for (Int_t icl=track0->GetFirstPoint();icl<track0->GetLastPoint(); icl++){
5371 if (track0->GetClusterPointer(icl)!=0){
5373 if (track0->GetClusterPointer(icl)->IsUsed(10)) ishared++;
5376 if (Float_t(ishared+1)/Float_t(all+1)>0.5) {
5377 delete array->RemoveAt(i);
5381 if (track0->GetKinkIndex(0)!=0) continue;
5382 if (track0->GetNumberOfClusters()<80) continue;
5384 AliTPCseed *pmother = new AliTPCseed();
5385 AliTPCseed *pdaughter = new AliTPCseed();
5386 AliKink *pkink = new AliKink;
5388 AliTPCseed & mother = *pmother;
5389 AliTPCseed & daughter = *pdaughter;
5390 AliKink & kinkl = *pkink;
5391 if (CheckKinkPoint(track0,mother,daughter, kinkl)){
5392 if (mother.GetNumberOfClusters()<30||daughter.GetNumberOfClusters()<20) {
5396 continue; //too short tracks
5398 if (mother.Pt()<1.4) {
5404 Int_t row0= kinkl.GetTPCRow0();
5405 if (kinkl.GetDistance()>0.5 || kinkl.GetR()<110. || kinkl.GetR()>240.) {
5412 Int_t index = esd->AddKink(&kinkl);
5413 mother.SetKinkIndex(0,-(index+1));
5414 daughter.SetKinkIndex(0,index+1);
5415 if (mother.GetNumberOfClusters()>50) {
5416 delete array->RemoveAt(i);
5417 array->AddAt(new AliTPCseed(mother),i);
5420 array->AddLast(new AliTPCseed(mother));
5422 array->AddLast(new AliTPCseed(daughter));
5423 for (Int_t icl=0;icl<row0;icl++) {
5424 if (mother.GetClusterPointer(icl)) mother.GetClusterPointer(icl)->Use(20);
5427 for (Int_t icl=row0;icl<158;icl++) {
5428 if (daughter.GetClusterPointer(icl)) daughter.GetClusterPointer(icl)->Use(20);
5437 delete [] daughters;
5459 AliInfo(Form("Ncandidates=\t%d\t%d\t%d\t%d\n",esd->GetNumberOfKinks(),ncandidates,ntracks,nall));
5464 Int_t AliTPCtrackerMI::RefitKink(AliTPCseed &mother, AliTPCseed &daughter, const AliESDkink &knk)
5467 // refit kink towards to the vertex
5470 AliKink &kink=(AliKink &)knk;
5472 Int_t row0 = GetRowNumber(kink.GetR());
5473 FollowProlongation(mother,0);
5474 mother.Reset(kFALSE);
5476 FollowProlongation(daughter,row0);
5477 daughter.Reset(kFALSE);
5478 FollowBackProlongation(daughter,158);
5479 daughter.Reset(kFALSE);
5480 Int_t first = TMath::Max(row0-20,30);
5481 Int_t last = TMath::Min(row0+20,140);
5483 const Int_t kNdiv =5;
5484 AliTPCseed param0[kNdiv]; // parameters along the track
5485 AliTPCseed param1[kNdiv]; // parameters along the track
5486 AliKink kinks[kNdiv]; // corresponding kink parameters
5489 for (Int_t irow=0; irow<kNdiv;irow++){
5490 rows[irow] = first +((last-first)*irow)/(kNdiv-1);
5492 // store parameters along the track
5494 for (Int_t irow=0;irow<kNdiv;irow++){
5495 FollowBackProlongation(mother, rows[irow]);
5496 FollowProlongation(daughter,rows[kNdiv-1-irow]);
5497 param0[irow] = mother;
5498 param1[kNdiv-1-irow] = daughter;
5502 for (Int_t irow=0; irow<kNdiv-1;irow++){
5503 if (param0[irow].GetNumberOfClusters()<kNdiv||param1[irow].GetNumberOfClusters()<kNdiv) continue;
5504 kinks[irow].SetMother(param0[irow]);
5505 kinks[irow].SetDaughter(param1[irow]);
5506 kinks[irow].Update();
5509 // choose kink with best "quality"
5511 Double_t mindist = 10000;
5512 for (Int_t irow=0;irow<kNdiv;irow++){
5513 if (param0[irow].GetNumberOfClusters()<20||param1[irow].GetNumberOfClusters()<20) continue;
5514 if (TMath::Abs(kinks[irow].GetR())>240.) continue;
5515 if (TMath::Abs(kinks[irow].GetR())<100.) continue;
5517 Float_t normdist = TMath::Abs(param0[irow].GetX()-kinks[irow].GetR())*(0.1+kink.GetDistance());
5518 normdist/= (param0[irow].GetNumberOfClusters()+param1[irow].GetNumberOfClusters()+40.);
5519 if (normdist < mindist){
5525 if (index==-1) return 0;
5528 param0[index].Reset(kTRUE);
5529 FollowProlongation(param0[index],0);
5531 mother = param0[index];
5532 daughter = param1[index]; // daughter in vertex
5534 kink.SetMother(mother);
5535 kink.SetDaughter(daughter);
5537 kink.SetTPCRow0(GetRowNumber(kink.GetR()));
5538 kink.SetTPCncls(param0[index].GetNumberOfClusters(),0);
5539 kink.SetTPCncls(param1[index].GetNumberOfClusters(),1);
5540 kink.SetLabel(CookLabel(&mother,0.4, 0,kink.GetTPCRow0()),0);
5541 kink.SetLabel(CookLabel(&daughter,0.4, kink.GetTPCRow0(),160),1);
5542 mother.SetLabel(kink.GetLabel(0));
5543 daughter.SetLabel(kink.GetLabel(1));
5549 void AliTPCtrackerMI::UpdateKinkQualityM(AliTPCseed * seed){
5551 // update Kink quality information for mother after back propagation
5553 if (seed->GetKinkIndex(0)>=0) return;
5554 for (Int_t ikink=0;ikink<3;ikink++){
5555 Int_t index = seed->GetKinkIndex(ikink);
5556 if (index>=0) break;
5557 index = TMath::Abs(index)-1;
5558 AliESDkink * kink = fEvent->GetKink(index);
5559 kink->SetTPCDensity(-1,0,0);
5560 kink->SetTPCDensity(1,0,1);
5562 Int_t row0 = kink->GetTPCRow0() - 2 - Int_t( 0.5/ (0.05+kink->GetAngle(2)));
5563 if (row0<15) row0=15;
5565 Int_t row1 = kink->GetTPCRow0() + 2 + Int_t( 0.5/ (0.05+kink->GetAngle(2)));
5566 if (row1>145) row1=145;
5568 Int_t found,foundable,shared;
5569 seed->GetClusterStatistic(0,row0, found, foundable,shared,kFALSE);
5570 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),0,0);
5571 seed->GetClusterStatistic(row1,155, found, foundable,shared,kFALSE);
5572 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),0,1);
5577 void AliTPCtrackerMI::UpdateKinkQualityD(AliTPCseed * seed){
5579 // update Kink quality information for daughter after refit
5581 if (seed->GetKinkIndex(0)<=0) return;
5582 for (Int_t ikink=0;ikink<3;ikink++){
5583 Int_t index = seed->GetKinkIndex(ikink);
5584 if (index<=0) break;
5585 index = TMath::Abs(index)-1;
5586 AliESDkink * kink = fEvent->GetKink(index);
5587 kink->SetTPCDensity(-1,1,0);
5588 kink->SetTPCDensity(-1,1,1);
5590 Int_t row0 = kink->GetTPCRow0() -2 - Int_t( 0.5/ (0.05+kink->GetAngle(2)));
5591 if (row0<15) row0=15;
5593 Int_t row1 = kink->GetTPCRow0() +2 + Int_t( 0.5/ (0.05+kink->GetAngle(2)));
5594 if (row1>145) row1=145;
5596 Int_t found,foundable,shared;
5597 seed->GetClusterStatistic(0,row0, found, foundable,shared,kFALSE);
5598 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),1,0);
5599 seed->GetClusterStatistic(row1,155, found, foundable,shared,kFALSE);
5600 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),1,1);
5606 Int_t AliTPCtrackerMI::CheckKinkPoint(AliTPCseed*seed,AliTPCseed &mother, AliTPCseed &daughter, const AliESDkink &knk)
5609 // check kink point for given track
5610 // if return value=0 kink point not found
5611 // otherwise seed0 correspond to mother particle
5612 // seed1 correspond to daughter particle
5613 // kink parameter of kink point
5614 AliKink &kink=(AliKink &)knk;
5616 Int_t middlerow = (seed->GetFirstPoint()+seed->GetLastPoint())/2;
5617 Int_t first = seed->GetFirstPoint();
5618 Int_t last = seed->GetLastPoint();
5619 if (last-first<20) return 0; // shortest length - 2*30 = 60 pad-rows
5622 AliTPCseed *seed1 = ReSeed(seed,middlerow+20, kTRUE); //middle of chamber
5623 if (!seed1) return 0;
5624 FollowProlongation(*seed1,seed->GetLastPoint()-20);
5625 seed1->Reset(kTRUE);
5626 FollowProlongation(*seed1,158);
5627 seed1->Reset(kTRUE);
5628 last = seed1->GetLastPoint();
5630 AliTPCseed *seed0 = new AliTPCseed(*seed);
5631 seed0->Reset(kFALSE);
5634 AliTPCseed param0[20]; // parameters along the track
5635 AliTPCseed param1[20]; // parameters along the track
5636 AliKink kinks[20]; // corresponding kink parameters
5638 for (Int_t irow=0; irow<20;irow++){
5639 rows[irow] = first +((last-first)*irow)/19;
5641 // store parameters along the track
5643 for (Int_t irow=0;irow<20;irow++){
5644 FollowBackProlongation(*seed0, rows[irow]);
5645 FollowProlongation(*seed1,rows[19-irow]);
5646 param0[irow] = *seed0;
5647 param1[19-irow] = *seed1;
5651 for (Int_t irow=0; irow<19;irow++){
5652 kinks[irow].SetMother(param0[irow]);
5653 kinks[irow].SetDaughter(param1[irow]);
5654 kinks[irow].Update();
5657 // choose kink with biggest change of angle
5659 Double_t maxchange= 0;
5660 for (Int_t irow=1;irow<19;irow++){
5661 if (TMath::Abs(kinks[irow].GetR())>240.) continue;
5662 if (TMath::Abs(kinks[irow].GetR())<110.) continue;
5663 Float_t quality = TMath::Abs(kinks[irow].GetAngle(2))/(3.+TMath::Abs(kinks[irow].GetR()-param0[irow].GetX()));
5664 if ( quality > maxchange){
5665 maxchange = quality;
5672 if (index<0) return 0;
5674 Int_t row0 = GetRowNumber(kinks[index].GetR()); //row 0 estimate
5675 seed0 = new AliTPCseed(param0[index]);
5676 seed1 = new AliTPCseed(param1[index]);
5677 seed0->Reset(kFALSE);
5678 seed1->Reset(kFALSE);
5679 seed0->ResetCovariance(10.);
5680 seed1->ResetCovariance(10.);
5681 FollowProlongation(*seed0,0);
5682 FollowBackProlongation(*seed1,158);
5683 mother = *seed0; // backup mother at position 0
5684 seed0->Reset(kFALSE);
5685 seed1->Reset(kFALSE);
5686 seed0->ResetCovariance(10.);
5687 seed1->ResetCovariance(10.);
5689 first = TMath::Max(row0-20,0);
5690 last = TMath::Min(row0+20,158);
5692 for (Int_t irow=0; irow<20;irow++){
5693 rows[irow] = first +((last-first)*irow)/19;
5695 // store parameters along the track
5697 for (Int_t irow=0;irow<20;irow++){
5698 FollowBackProlongation(*seed0, rows[irow]);
5699 FollowProlongation(*seed1,rows[19-irow]);
5700 param0[irow] = *seed0;
5701 param1[19-irow] = *seed1;
5705 for (Int_t irow=0; irow<19;irow++){
5706 kinks[irow].SetMother(param0[irow]);
5707 kinks[irow].SetDaughter(param1[irow]);
5708 // param0[irow].Dump();
5709 //param1[irow].Dump();
5710 kinks[irow].Update();
5713 // choose kink with biggest change of angle
5716 for (Int_t irow=0;irow<20;irow++){
5717 if (TMath::Abs(kinks[irow].GetR())>250.) continue;
5718 if (TMath::Abs(kinks[irow].GetR())<90.) continue;
5719 Float_t quality = TMath::Abs(kinks[irow].GetAngle(2))/(3.+TMath::Abs(kinks[irow].GetR()-param0[irow].GetX()));
5720 if ( quality > maxchange){
5721 maxchange = quality;
5728 if (index==-1 || param0[index].GetNumberOfClusters()+param1[index].GetNumberOfClusters()<100){
5734 // Float_t anglesigma = TMath::Sqrt(param0[index].fC22+param0[index].fC33+param1[index].fC22+param1[index].fC33);
5736 kink.SetMother(param0[index]);
5737 kink.SetDaughter(param1[index]);
5740 Double_t chi2P2 = param0[index].GetParameter()[2]-param1[index].GetParameter()[2];
5742 chi2P2/=param0[index].GetCovariance()[5]+param1[index].GetCovariance()[5];
5743 Double_t chi2P3 = param0[index].GetParameter()[3]-param1[index].GetParameter()[3];
5745 chi2P3/=param0[index].GetCovariance()[9]+param1[index].GetCovariance()[9];
5747 if (AliTPCReconstructor::StreamLevel()>1) {
5748 (*fDebugStreamer)<<"kinkHpt"<<
5751 "p0.="<<¶m0[index]<<
5752 "p1.="<<¶m1[index]<<
5756 if ( chi2P2+chi2P3<AliTPCReconstructor::GetRecoParam()->GetKinkAngleCutChi2(0)){
5763 row0 = GetRowNumber(kink.GetR());
5764 kink.SetTPCRow0(row0);
5765 kink.SetLabel(CookLabel(seed0,0.5,0,row0),0);
5766 kink.SetLabel(CookLabel(seed1,0.5,row0,158),1);
5767 kink.SetIndex(-10,0);
5768 kink.SetIndex(int(param0[index].GetNumberOfClusters()+param1[index].GetNumberOfClusters()),1);
5769 kink.SetTPCncls(param0[index].GetNumberOfClusters(),0);
5770 kink.SetTPCncls(param1[index].GetNumberOfClusters(),1);
5773 // new (&mother) AliTPCseed(param0[index]);
5774 daughter = param1[index];
5775 daughter.SetLabel(kink.GetLabel(1));
5776 param0[index].Reset(kTRUE);
5777 FollowProlongation(param0[index],0);
5778 mother = param0[index];
5779 mother.SetLabel(kink.GetLabel(0));
5780 if ( chi2P2+chi2P3<AliTPCReconstructor::GetRecoParam()->GetKinkAngleCutChi2(1)){
5792 AliTPCseed* AliTPCtrackerMI::ReSeed(AliTPCseed *t)
5795 // reseed - refit - track
5798 // Int_t last = fSectors->GetNRows()-1;
5800 if (fSectors == fOuterSec){
5801 first = TMath::Max(first, t->GetFirstPoint()-fInnerSec->GetNRows());
5805 first = t->GetFirstPoint();
5807 AliTPCseed * seed = MakeSeed(t,0.1,0.5,0.9);
5808 FollowBackProlongation(*t,fSectors->GetNRows()-1);
5810 FollowProlongation(*t,first);
5820 //_____________________________________________________________________________
5821 Int_t AliTPCtrackerMI::ReadSeeds(const TFile *inp) {
5822 //-----------------------------------------------------------------
5823 // This function reades track seeds.
5824 //-----------------------------------------------------------------
5825 TDirectory *savedir=gDirectory;
5827 TFile *in=(TFile*)inp;
5828 if (!in->IsOpen()) {
5829 cerr<<"AliTPCtrackerMI::ReadSeeds(): input file is not open !\n";
5834 TTree *seedTree=(TTree*)in->Get("Seeds");
5836 cerr<<"AliTPCtrackerMI::ReadSeeds(): ";
5837 cerr<<"can't get a tree with track seeds !\n";
5840 AliTPCtrack *seed=new AliTPCtrack;
5841 seedTree->SetBranchAddress("tracks",&seed);
5843 if (fSeeds==0) fSeeds=new TObjArray(15000);
5845 Int_t n=(Int_t)seedTree->GetEntries();
5846 for (Int_t i=0; i<n; i++) {
5847 seedTree->GetEvent(i);
5848 fSeeds->AddLast(new AliTPCseed(*seed/*,seed->GetAlpha()*/));
5857 Int_t AliTPCtrackerMI::Clusters2Tracks (AliESDEvent *const esd)
5861 if (fSeeds) DeleteSeeds();
5864 if (!fSeeds) return 1;
5866 if (AliTPCReconstructor::StreamLevel()>3) DumpClusters(0,fSeeds);
5872 //_____________________________________________________________________________
5873 Int_t AliTPCtrackerMI::Clusters2Tracks() {
5874 //-----------------------------------------------------------------
5875 // This is a track finder.
5876 //-----------------------------------------------------------------
5877 TDirectory *savedir=gDirectory;
5881 fSeeds = Tracking();
5884 Info("Clusters2Tracks","Time for tracking: \t");timer.Print();timer.Start();
5886 //activate again some tracks
5887 for (Int_t i=0; i<fSeeds->GetEntriesFast(); i++) {
5888 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
5890 Int_t nc=t.GetNumberOfClusters();
5892 delete fSeeds->RemoveAt(i);
5896 if (pt->GetRemoval()==10) {
5897 if (pt->GetDensityFirst(20)>0.8 || pt->GetDensityFirst(30)>0.8 || pt->GetDensityFirst(40)>0.7)
5898 pt->Desactivate(10); // make track again active
5900 pt->Desactivate(20);
5901 delete fSeeds->RemoveAt(i);
5906 RemoveUsed2(fSeeds,0.85,0.85,0);
5907 if (AliTPCReconstructor::GetRecoParam()->GetDoKinks()) FindKinks(fSeeds,fEvent);
5908 //FindCurling(fSeeds, fEvent,0);
5909 if (AliTPCReconstructor::StreamLevel()>5) FindMultiMC(fSeeds, fEvent,-1); // find multi found tracks
5910 RemoveUsed2(fSeeds,0.5,0.4,20);
5911 FindSplitted(fSeeds, fEvent,0); // find multi found tracks
5912 if (AliTPCReconstructor::StreamLevel()>5) FindMultiMC(fSeeds, fEvent,0); // find multi found tracks
5915 // // refit short tracks
5917 Int_t nseed=fSeeds->GetEntriesFast();
5920 for (Int_t i=0; i<nseed; i++) {
5921 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
5923 Int_t nc=t.GetNumberOfClusters();
5925 delete fSeeds->RemoveAt(i);
5928 CookLabel(pt,0.1); //For comparison only
5929 //if ((pt->IsActive() || (pt->fRemoval==10) )&& nc>50 &&pt->GetNumberOfClusters()>0.4*pt->fNFoundable){
5930 if ((pt->IsActive() || (pt->GetRemoval()==10) )){
5932 if (fDebug>0) cerr<<found<<'\r';
5936 delete fSeeds->RemoveAt(i);
5940 //RemoveOverlap(fSeeds,0.99,7,kTRUE);
5942 //RemoveUsed(fSeeds,0.9,0.9,6);
5944 nseed=fSeeds->GetEntriesFast();
5946 for (Int_t i=0; i<nseed; i++) {
5947 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
5949 Int_t nc=t.GetNumberOfClusters();
5951 delete fSeeds->RemoveAt(i);
5955 t.CookdEdx(0.02,0.6);
5956 // CheckKinkPoint(&t,0.05);
5957 //if ((pt->IsActive() || (pt->fRemoval==10) )&& nc>50 &&pt->GetNumberOfClusters()>0.4*pt->fNFoundable){
5958 if ((pt->IsActive() || (pt->GetRemoval()==10) )){
5966 delete fSeeds->RemoveAt(i);
5967 //AliTPCseed * seed1 = ReSeed(pt,0.05,0.5,1);
5969 // FollowProlongation(*seed1,0);
5970 // Int_t n = seed1->GetNumberOfClusters();
5971 // printf("fP4\t%f\t%f\n",seed1->GetC(),pt->GetC());
5972 // printf("fN\t%d\t%d\n", seed1->GetNumberOfClusters(),pt->GetNumberOfClusters());
5975 //AliTPCseed * seed2 = ReSeed(pt,0.95,0.5,0.05);
5979 SortTracks(fSeeds, 1);
5983 PrepareForBackProlongation(fSeeds,5.);
5984 PropagateBack(fSeeds);
5985 printf("Time for back propagation: \t");timer.Print();timer.Start();
5989 PrepareForProlongation(fSeeds,5.);
5990 PropagateForward2(fSeeds);
5992 printf("Time for FORWARD propagation: \t");timer.Print();timer.Start();
5993 // RemoveUsed(fSeeds,0.7,0.7,6);
5994 //RemoveOverlap(fSeeds,0.9,7,kTRUE);
5996 nseed=fSeeds->GetEntriesFast();
5998 for (Int_t i=0; i<nseed; i++) {
5999 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
6001 Int_t nc=t.GetNumberOfClusters();
6003 delete fSeeds->RemoveAt(i);
6006 t.CookdEdx(0.02,0.6);
6007 // CookLabel(pt,0.1); //For comparison only
6008 //if ((pt->IsActive() || (pt->fRemoval==10) )&& nc>50 &&pt->GetNumberOfClusters()>0.4*pt->fNFoundable){
6009 if ((pt->IsActive() || (pt->fRemoval==10) )){
6010 cerr<<found++<<'\r';
6013 delete fSeeds->RemoveAt(i);
6018 // fNTracks = found;
6020 Info("Clusters2Tracks","Time for overlap removal, track writing and dedx cooking: \t"); timer.Print();timer.Start();
6023 // cerr<<"Number of found tracks : "<<"\t"<<found<<endl;
6024 Info("Clusters2Tracks","Number of found tracks %d",found);
6026 // UnloadClusters();
6031 void AliTPCtrackerMI::Tracking(TObjArray * arr)
6034 // tracking of the seeds
6037 fSectors = fOuterSec;
6038 ParallelTracking(arr,150,63);
6039 fSectors = fOuterSec;
6040 ParallelTracking(arr,63,0);
6043 TObjArray * AliTPCtrackerMI::Tracking(Int_t seedtype, Int_t i1, Int_t i2, Float_t cuts[4], Float_t dy, Int_t dsec)
6048 TObjArray * arr = new TObjArray;
6050 fSectors = fOuterSec;
6053 for (Int_t sec=0;sec<fkNOS;sec++){
6054 if (seedtype==3) MakeSeeds3(arr,sec,i1,i2,cuts,dy, dsec);
6055 if (seedtype==4) MakeSeeds5(arr,sec,i1,i2,cuts,dy);
6056 if (seedtype==2) MakeSeeds2(arr,sec,i1,i2,cuts,dy);
6059 Info("Tracking","\nSeeding - %d\t%d\t%d\t%d\n",seedtype,i1,i2,arr->GetEntriesFast());
6071 TObjArray * AliTPCtrackerMI::Tracking()
6075 if (AliTPCReconstructor::GetRecoParam()->GetSpecialSeeding()) return TrackingSpecial();
6078 Int_t nup=fOuterSec->GetNRows()+fInnerSec->GetNRows();
6080 TObjArray * seeds = new TObjArray;
6082 Int_t fLastSeedRowSec=AliTPCReconstructor::GetRecoParam()->GetLastSeedRowSec();
6083 Int_t gapPrim = AliTPCReconstructor::GetRecoParam()->GetSeedGapPrim();
6084 Int_t gapSec = AliTPCReconstructor::GetRecoParam()->GetSeedGapSec();
6092 Float_t fnumber = 3.0;
6093 Float_t fdensity = 3.0;
6098 for (Int_t delta = 0; delta<18; delta+=gapPrim){
6102 arr = Tracking(3,nup-1-delta,nup-1-delta-gap,cuts,-1,1);
6103 SumTracks(seeds,arr);
6104 SignClusters(seeds,fnumber,fdensity);
6106 for (Int_t i=2;i<6;i+=2){
6107 // seed high pt tracks
6110 arr = Tracking(3,nup-i-delta,nup-i-delta-gap,cuts,-1,0);
6111 SumTracks(seeds,arr);
6112 SignClusters(seeds,fnumber,fdensity);
6117 // RemoveUsed(seeds,0.9,0.9,1);
6118 // UnsignClusters();
6119 // SignClusters(seeds,fnumber,fdensity);
6123 for (Int_t delta = 20; delta<120; delta+=gapPrim){
6125 // seed high pt tracks
6129 arr = Tracking(3,nup-delta,nup-delta-gap,cuts,-1);
6130 SumTracks(seeds,arr);
6131 SignClusters(seeds,fnumber,fdensity);
6136 arr = Tracking(3,nup-delta-5,nup-delta-5-gap,cuts,-1);
6137 SumTracks(seeds,arr);
6138 SignClusters(seeds,fnumber,fdensity);
6149 Info("Tracking()","\n\nPrimary seeding\t%d\n\n",seeds->GetEntriesFast());
6153 // RemoveUsed(seeds,0.75,0.75,1);
6155 //SignClusters(seeds,fnumber,fdensity);
6164 arr = Tracking(4,nup-1,nup-1-gap,cuts,-1);
6165 SumTracks(seeds,arr);
6166 SignClusters(seeds,fnumber,fdensity);
6168 arr = Tracking(4,nup-2,nup-2-gap,cuts,-1);
6169 SumTracks(seeds,arr);
6170 SignClusters(seeds,fnumber,fdensity);
6172 arr = Tracking(4,nup-3,nup-3-gap,cuts,-1);
6173 SumTracks(seeds,arr);
6174 SignClusters(seeds,fnumber,fdensity);
6176 arr = Tracking(4,nup-5,nup-5-gap,cuts,-1);
6177 SumTracks(seeds,arr);
6178 SignClusters(seeds,fnumber,fdensity);
6180 arr = Tracking(4,nup-7,nup-7-gap,cuts,-1);
6181 SumTracks(seeds,arr);
6182 SignClusters(seeds,fnumber,fdensity);
6185 arr = Tracking(4,nup-9,nup-9-gap,cuts,-1);
6186 SumTracks(seeds,arr);
6187 SignClusters(seeds,fnumber,fdensity);
6191 for (Int_t delta = 9; delta<30; delta+=gapSec){
6197 arr = Tracking(4,nup-1-delta,nup-1-delta-gap,cuts,-1);
6198 SumTracks(seeds,arr);
6199 SignClusters(seeds,fnumber,fdensity);
6201 arr = Tracking(4,nup-3-delta,nup-5-delta-gap,cuts,4);
6202 SumTracks(seeds,arr);
6203 SignClusters(seeds,fnumber,fdensity);
6216 for (Int_t delta = 30; delta<fLastSeedRowSec; delta+=gapSec){
6222 arr = Tracking(4,nup-1-delta,nup-1-delta-gap,cuts,-1);
6223 SumTracks(seeds,arr);
6224 SignClusters(seeds,fnumber,fdensity);
6226 arr = Tracking(4,nup-5-delta,nup-5-delta-gap,cuts,5 );
6227 SumTracks(seeds,arr);
6228 SignClusters(seeds,fnumber,fdensity);
6232 Info("Tracking()","\n\nSecondary seeding\t%d\n\n",seeds->GetEntriesFast());
6243 TObjArray * AliTPCtrackerMI::TrackingSpecial()
6246 // seeding adjusted for laser and cosmic tests - short tracks with big inclination angle
6247 // no primary vertex seeding tried
6251 Int_t nup=fOuterSec->GetNRows()+fInnerSec->GetNRows();
6253 TObjArray * seeds = new TObjArray;
6258 Float_t fnumber = 3.0;
6259 Float_t fdensity = 3.0;
6262 cuts[0] = AliTPCReconstructor::GetRecoParam()->GetMaxC(); // max curvature
6263 cuts[1] = 3.5; // max tan(phi) angle for seeding
6264 cuts[2] = 3.; // not used (cut on z primary vertex)
6265 cuts[3] = 3.5; // max tan(theta) angle for seeding
6267 for (Int_t delta = 0; nup-delta-gap-1>0; delta+=3){
6269 arr = Tracking(4,nup-1-delta,nup-1-delta-gap,cuts,-1);
6270 SumTracks(seeds,arr);
6271 SignClusters(seeds,fnumber,fdensity);
6275 Info("Tracking()","\n\nSecondary seeding\t%d\n\n",seeds->GetEntriesFast());
6286 void AliTPCtrackerMI::SumTracks(TObjArray *arr1,TObjArray *&arr2) const
6289 //sum tracks to common container
6290 //remove suspicious tracks
6291 Int_t nseed = arr2->GetEntriesFast();
6292 for (Int_t i=0;i<nseed;i++){
6293 AliTPCseed *pt=(AliTPCseed*)arr2->UncheckedAt(i);
6296 // remove tracks with too big curvature
6298 if (TMath::Abs(pt->GetC())>AliTPCReconstructor::GetRecoParam()->GetMaxC()){
6299 delete arr2->RemoveAt(i);
6302 // REMOVE VERY SHORT TRACKS
6303 if (pt->GetNumberOfClusters()<20){
6304 delete arr2->RemoveAt(i);
6307 // NORMAL ACTIVE TRACK
6308 if (pt->IsActive()){
6309 arr1->AddLast(arr2->RemoveAt(i));
6312 //remove not usable tracks
6313 if (pt->GetRemoval()!=10){
6314 delete arr2->RemoveAt(i);
6318 // ENABLE ONLY ENOUGH GOOD STOPPED TRACKS
6319 if (pt->GetDensityFirst(20)>0.8 || pt->GetDensityFirst(30)>0.8 || pt->GetDensityFirst(40)>0.7)
6320 arr1->AddLast(arr2->RemoveAt(i));
6322 delete arr2->RemoveAt(i);
6326 delete arr2; arr2 = 0;
6331 void AliTPCtrackerMI::ParallelTracking(TObjArray *const arr, Int_t rfirst, Int_t rlast)
6334 // try to track in parralel
6336 Int_t nseed=arr->GetEntriesFast();
6337 //prepare seeds for tracking
6338 for (Int_t i=0; i<nseed; i++) {
6339 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i), &t=*pt;
6341 if (!t.IsActive()) continue;
6342 // follow prolongation to the first layer
6343 if ( (fSectors ==fInnerSec) || (t.GetFirstPoint()-fkParam->GetNRowLow()>rfirst+1) )
6344 FollowProlongation(t, rfirst+1);
6349 for (Int_t nr=rfirst; nr>=rlast; nr--){
6350 if (nr<fInnerSec->GetNRows())
6351 fSectors = fInnerSec;
6353 fSectors = fOuterSec;
6354 // make indexes with the cluster tracks for given
6356 // find nearest cluster
6357 for (Int_t i=0; i<nseed; i++) {
6358 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i), &t=*pt;
6360 if (nr==80) pt->UpdateReference();
6361 if (!pt->IsActive()) continue;
6362 // if ( (fSectors ==fOuterSec) && (pt->fFirstPoint-fkParam->GetNRowLow())<nr) continue;
6363 if (pt->GetRelativeSector()>17) {
6366 UpdateClusters(t,nr);
6368 // prolonagate to the nearest cluster - if founded
6369 for (Int_t i=0; i<nseed; i++) {
6370 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
6372 if (!pt->IsActive()) continue;
6373 // if ((fSectors ==fOuterSec) && (pt->fFirstPoint-fkParam->GetNRowLow())<nr) continue;
6374 if (pt->GetRelativeSector()>17) {
6377 FollowToNextCluster(*pt,nr);
6382 void AliTPCtrackerMI::PrepareForBackProlongation(TObjArray *const arr,Float_t fac) const
6386 // if we use TPC track itself we have to "update" covariance
6388 Int_t nseed= arr->GetEntriesFast();
6389 for (Int_t i=0;i<nseed;i++){
6390 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
6394 //rotate to current local system at first accepted point
6395 Int_t index = pt->GetClusterIndex2(pt->GetFirstPoint());
6396 Int_t sec = (index&0xff000000)>>24;
6398 Float_t angle1 = fInnerSec->GetAlpha()*sec+fInnerSec->GetAlphaShift();
6399 if (angle1>TMath::Pi())
6400 angle1-=2.*TMath::Pi();
6401 Float_t angle2 = pt->GetAlpha();
6403 if (TMath::Abs(angle1-angle2)>0.001){
6404 pt->Rotate(angle1-angle2);
6405 //angle2 = pt->GetAlpha();
6406 //pt->fRelativeSector = pt->GetAlpha()/fInnerSec->GetAlpha();
6407 //if (pt->GetAlpha()<0)
6408 // pt->fRelativeSector+=18;
6409 //sec = pt->fRelativeSector;
6418 void AliTPCtrackerMI::PrepareForProlongation(TObjArray *const arr, Float_t fac) const
6422 // if we use TPC track itself we have to "update" covariance
6424 Int_t nseed= arr->GetEntriesFast();
6425 for (Int_t i=0;i<nseed;i++){
6426 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
6429 pt->SetFirstPoint(pt->GetLastPoint());
6437 Int_t AliTPCtrackerMI::PropagateBack(TObjArray *const arr)
6440 // make back propagation
6442 Int_t nseed= arr->GetEntriesFast();
6443 for (Int_t i=0;i<nseed;i++){
6444 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
6445 if (pt&& pt->GetKinkIndex(0)<=0) {
6446 //AliTPCseed *pt2 = new AliTPCseed(*pt);
6447 fSectors = fInnerSec;
6448 //FollowBackProlongation(*pt,fInnerSec->GetNRows()-1);
6449 //fSectors = fOuterSec;
6450 FollowBackProlongation(*pt,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1);
6451 //if (pt->GetNumberOfClusters()<(pt->fEsd->GetTPCclusters(0)) ){
6452 // Error("PropagateBack","Not prolonged track %d",pt->GetLabel());
6453 // FollowBackProlongation(*pt2,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1);
6456 if (pt&& pt->GetKinkIndex(0)>0) {
6457 AliESDkink * kink = fEvent->GetKink(pt->GetKinkIndex(0)-1);
6458 pt->SetFirstPoint(kink->GetTPCRow0());
6459 fSectors = fInnerSec;
6460 FollowBackProlongation(*pt,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1);
6468 Int_t AliTPCtrackerMI::PropagateForward2(TObjArray *const arr)
6471 // make forward propagation
6473 Int_t nseed= arr->GetEntriesFast();
6475 for (Int_t i=0;i<nseed;i++){
6476 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
6478 FollowProlongation(*pt,0);
6487 Int_t AliTPCtrackerMI::PropagateForward()
6490 // propagate track forward
6492 Int_t nseed = fSeeds->GetEntriesFast();
6493 for (Int_t i=0;i<nseed;i++){
6494 AliTPCseed *pt = (AliTPCseed*)fSeeds->UncheckedAt(i);
6496 AliTPCseed &t = *pt;
6497 Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
6498 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
6499 if (alpha < 0. ) alpha += 2.*TMath::Pi();
6500 t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN);
6504 fSectors = fOuterSec;
6505 ParallelTracking(fSeeds,fOuterSec->GetNRows()+fInnerSec->GetNRows()-1,fInnerSec->GetNRows());
6506 fSectors = fInnerSec;
6507 ParallelTracking(fSeeds,fInnerSec->GetNRows()-1,0);
6516 Int_t AliTPCtrackerMI::PropagateBack(AliTPCseed *const pt, Int_t row0, Int_t row1)
6519 // make back propagation, in between row0 and row1
6523 fSectors = fInnerSec;
6526 if (row1<fSectors->GetNRows())
6529 r1 = fSectors->GetNRows()-1;
6531 if (row0<fSectors->GetNRows()&& r1>0 )
6532 FollowBackProlongation(*pt,r1);
6533 if (row1<=fSectors->GetNRows())
6536 r1 = row1 - fSectors->GetNRows();
6537 if (r1<=0) return 0;
6538 if (r1>=fOuterSec->GetNRows()) return 0;
6539 fSectors = fOuterSec;
6540 return FollowBackProlongation(*pt,r1);
6548 void AliTPCtrackerMI::GetShape(AliTPCseed * seed, Int_t row)
6552 AliTPCClusterParam * clparam = AliTPCcalibDB::Instance()->GetClusterParam();
6553 Float_t zdrift = TMath::Abs((fkParam->GetZLength(0)-TMath::Abs(seed->GetZ())));
6554 Int_t type = (seed->GetSector() < fkParam->GetNSector()/2) ? 0: (row>126) ? 1:2;
6555 Double_t angulary = seed->GetSnp();
6557 if (TMath::Abs(angulary)>AliTPCReconstructor::GetMaxSnpTracker()) {
6558 angulary = TMath::Sign(AliTPCReconstructor::GetMaxSnpTracker(),angulary);
6561 angulary = angulary*angulary/((1.-angulary)*(1.+angulary));
6562 Double_t angularz = seed->GetTgl()*seed->GetTgl()*(1.+angulary);
6564 Double_t sigmay = clparam->GetRMS0(0,type,zdrift,TMath::Sqrt(TMath::Abs(angulary)));
6565 Double_t sigmaz = clparam->GetRMS0(1,type,zdrift,TMath::Sqrt(TMath::Abs(angularz)));
6566 seed->SetCurrentSigmaY2(sigmay*sigmay);
6567 seed->SetCurrentSigmaZ2(sigmaz*sigmaz);
6568 // Float_t sd2 = TMath::Abs((fkParam->GetZLength(0)-TMath::Abs(seed->GetZ())))*fkParam->GetDiffL()*fkParam->GetDiffL();
6569 // // Float_t padlength = fkParam->GetPadPitchLength(seed->fSector);
6570 // Float_t padlength = GetPadPitchLength(row);
6572 // Float_t sresy = (seed->GetSector() < fkParam->GetNSector()/2) ? 0.2 :0.3;
6573 // seed->SetCurrentSigmaY2(sd2+padlength*padlength*angulary/12.+sresy*sresy);
6575 // Float_t sresz = fkParam->GetZSigma();
6576 // seed->SetCurrentSigmaZ2(sd2+padlength*padlength*angularz*angularz*(1+angulary)/12.+sresz*sresz);
6578 Float_t wy = GetSigmaY(seed);
6579 Float_t wz = GetSigmaZ(seed);
6582 if (TMath::Abs(wy/seed->fCurrentSigmaY2-1)>0.0001 || TMath::Abs(wz/seed->fCurrentSigmaZ2-1)>0.0001 ){
6583 printf("problem\n");
6590 //__________________________________________________________________________
6591 void AliTPCtrackerMI::CookLabel(AliKalmanTrack *tk, Float_t wrong) const {
6592 //--------------------------------------------------------------------
6593 //This function "cooks" a track label. If label<0, this track is fake.
6594 //--------------------------------------------------------------------
6595 AliTPCseed * t = dynamic_cast<AliTPCseed*>(tk);
6597 printf("%s:%d wrong type \n",(char*)__FILE__,__LINE__);
6601 Int_t noc=t->GetNumberOfClusters();
6603 //printf("\nnot founded prolongation\n\n\n");
6609 AliTPCclusterMI *clusters[160];
6611 for (Int_t i=0;i<160;i++) {
6618 for (i=0; i<160 && current<noc; i++) {
6620 Int_t index=t->GetClusterIndex2(i);
6621 if (index<=0) continue;
6622 if (index&0x8000) continue;
6624 //clusters[current]=GetClusterMI(index);
6625 if (t->GetClusterPointer(i)){
6626 clusters[current]=t->GetClusterPointer(i);
6632 Int_t lab=123456789;
6633 for (i=0; i<noc; i++) {
6634 AliTPCclusterMI *c=clusters[i];
6636 lab=TMath::Abs(c->GetLabel(0));
6638 for (j=0; j<noc; j++) if (lb[j]==lab || mx[j]==0) break;
6644 for (i=0; i<noc; i++) if (mx[i]>max) {max=mx[i]; lab=lb[i];}
6646 for (i=0; i<noc; i++) {
6647 AliTPCclusterMI *c=clusters[i];
6649 if (TMath::Abs(c->GetLabel(1)) == lab ||
6650 TMath::Abs(c->GetLabel(2)) == lab ) max++;
6652 if (noc<=0) { lab=-1; return;}
6653 if ((1.- Float_t(max)/(noc)) > wrong) lab=-lab;
6656 Int_t tail=Int_t(0.10*noc);
6659 for (i=1; i<160&&ind<tail; i++) {
6660 // AliTPCclusterMI *c=clusters[noc-i];
6661 AliTPCclusterMI *c=clusters[i];
6663 if (lab == TMath::Abs(c->GetLabel(0)) ||
6664 lab == TMath::Abs(c->GetLabel(1)) ||
6665 lab == TMath::Abs(c->GetLabel(2))) max++;
6668 if (max < Int_t(0.5*tail)) lab=-lab;
6675 //delete[] clusters;
6679 //__________________________________________________________________________
6680 Int_t AliTPCtrackerMI::CookLabel(AliTPCseed *const t, Float_t wrong,Int_t first, Int_t last) const {
6681 //--------------------------------------------------------------------
6682 //This function "cooks" a track label. If label<0, this track is fake.
6683 //--------------------------------------------------------------------
6684 Int_t noc=t->GetNumberOfClusters();
6686 //printf("\nnot founded prolongation\n\n\n");
6692 AliTPCclusterMI *clusters[160];
6694 for (Int_t i=0;i<160;i++) {
6701 for (i=0; i<160 && current<noc; i++) {
6702 if (i<first) continue;
6703 if (i>last) continue;
6704 Int_t index=t->GetClusterIndex2(i);
6705 if (index<=0) continue;
6706 if (index&0x8000) continue;
6708 //clusters[current]=GetClusterMI(index);
6709 if (t->GetClusterPointer(i)){
6710 clusters[current]=t->GetClusterPointer(i);
6715 if (noc<5) return -1;
6716 Int_t lab=123456789;
6717 for (i=0; i<noc; i++) {
6718 AliTPCclusterMI *c=clusters[i];
6720 lab=TMath::Abs(c->GetLabel(0));
6722 for (j=0; j<noc; j++) if (lb[j]==lab || mx[j]==0) break;
6728 for (i=0; i<noc; i++) if (mx[i]>max) {max=mx[i]; lab=lb[i];}
6730 for (i=0; i<noc; i++) {
6731 AliTPCclusterMI *c=clusters[i];
6733 if (TMath::Abs(c->GetLabel(1)) == lab ||
6734 TMath::Abs(c->GetLabel(2)) == lab ) max++;
6736 if (noc<=0) { lab=-1; return -1;}
6737 if ((1.- Float_t(max)/(noc)) > wrong) lab=-lab;
6740 Int_t tail=Int_t(0.10*noc);
6743 for (i=1; i<160&&ind<tail; i++) {
6744 // AliTPCclusterMI *c=clusters[noc-i];
6745 AliTPCclusterMI *c=clusters[i];
6747 if (lab == TMath::Abs(c->GetLabel(0)) ||
6748 lab == TMath::Abs(c->GetLabel(1)) ||
6749 lab == TMath::Abs(c->GetLabel(2))) max++;
6752 if (max < Int_t(0.5*tail)) lab=-lab;
6755 // t->SetLabel(lab);
6759 //delete[] clusters;
6763 Int_t AliTPCtrackerMI::GetRowNumber(Double_t x[3]) const
6765 //return pad row number for given x vector
6766 Float_t phi = TMath::ATan2(x[1],x[0]);
6767 if(phi<0) phi=2.*TMath::Pi()+phi;
6768 // Get the local angle in the sector philoc
6769 const Float_t kRaddeg = 180/3.14159265358979312;
6770 Float_t phiangle = (Int_t (phi*kRaddeg/20.) + 0.5)*20./kRaddeg;
6771 Double_t localx = x[0]*TMath::Cos(phiangle)-x[1]*TMath::Sin(phiangle);
6772 return GetRowNumber(localx);
6777 void AliTPCtrackerMI::MakeBitmaps(AliTPCseed *t)
6779 //-----------------------------------------------------------------------
6780 // Fill the cluster and sharing bitmaps of the track
6781 //-----------------------------------------------------------------------
6783 Int_t firstpoint = 0;
6784 Int_t lastpoint = 159;
6785 AliTPCTrackerPoint *point;
6786 AliTPCclusterMI *cluster;
6788 for (int iter=firstpoint; iter<lastpoint; iter++) {
6789 // Change to cluster pointers to see if we have a cluster at given padrow
6790 cluster = t->GetClusterPointer(iter);
6792 t->SetClusterMapBit(iter, kTRUE);
6793 point = t->GetTrackPoint(iter);
6794 if (point->IsShared())
6795 t->SetSharedMapBit(iter,kTRUE);
6797 t->SetSharedMapBit(iter, kFALSE);
6800 t->SetClusterMapBit(iter, kFALSE);
6801 t->SetSharedMapBit(iter, kFALSE);
6806 Bool_t AliTPCtrackerMI::IsFindable(AliTPCseed & track){
6808 // return flag if there is findable cluster at given position
6811 Float_t z = track.GetZ();
6813 if (TMath::Abs(z)<(AliTPCReconstructor::GetCtgRange()*track.GetX()+kDeltaZ) &&
6814 TMath::Abs(z)<fkParam->GetZLength(0) &&
6815 (TMath::Abs(track.GetSnp())<AliTPCReconstructor::GetMaxSnpTracker()))
6821 void AliTPCtrackerMI::AddCovariance(AliTPCseed * seed){
6823 // Adding systematic error
6824 // !!!! the systematic error for element 4 is in 1/cm not in pt
6826 const Double_t *param = AliTPCReconstructor::GetRecoParam()->GetSystematicError();
6827 Double_t *covarIn= (Double_t*)seed->GetCovariance();
6829 for (Int_t i=0;i<15;i++) covar[i]=0;
6835 covar[0] = param[0]*param[0];
6836 covar[2] = param[1]*param[1];
6837 covar[5] = param[2]*param[2];
6838 covar[9] = param[3]*param[3];
6839 Double_t facC = AliTracker::GetBz()*kB2C;
6840 covar[14]= param[4]*param[4]*facC*facC;
6842 covar[1]=TMath::Sqrt((covar[0]*covar[2]))*covarIn[1]/TMath::Sqrt((covarIn[0]*covarIn[2]));
6844 covar[3]=TMath::Sqrt((covar[0]*covar[5]))*covarIn[3]/TMath::Sqrt((covarIn[0]*covarIn[5]));
6845 covar[4]=TMath::Sqrt((covar[2]*covar[5]))*covarIn[4]/TMath::Sqrt((covarIn[2]*covarIn[5]));
6847 covar[6]=TMath::Sqrt((covar[0]*covar[9]))*covarIn[6]/TMath::Sqrt((covarIn[0]*covarIn[9]));
6848 covar[7]=TMath::Sqrt((covar[2]*covar[9]))*covarIn[7]/TMath::Sqrt((covarIn[2]*covarIn[9]));
6849 covar[8]=TMath::Sqrt((covar[5]*covar[9]))*covarIn[8]/TMath::Sqrt((covarIn[5]*covarIn[9]));
6851 covar[10]=TMath::Sqrt((covar[0]*covar[14]))*covarIn[10]/TMath::Sqrt((covarIn[0]*covarIn[14]));
6852 covar[11]=TMath::Sqrt((covar[2]*covar[14]))*covarIn[11]/TMath::Sqrt((covarIn[2]*covarIn[14]));
6853 covar[12]=TMath::Sqrt((covar[5]*covar[14]))*covarIn[12]/TMath::Sqrt((covarIn[5]*covarIn[14]));
6854 covar[13]=TMath::Sqrt((covar[9]*covar[14]))*covarIn[13]/TMath::Sqrt((covarIn[9]*covarIn[14]));
6856 seed->AddCovariance(covar);