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 //------------------------------------------------------------------
460 for (Int_t irow=0; irow<200; irow++){
467 AliTPCtrackerMI & AliTPCtrackerMI::operator=(const AliTPCtrackerMI& /*r*/)
469 //------------------------------
471 //--------------------------------------------------------------
474 //_____________________________________________________________________________
475 AliTPCtrackerMI::~AliTPCtrackerMI() {
476 //------------------------------------------------------------------
477 // TPC tracker destructor
478 //------------------------------------------------------------------
485 if (fDebugStreamer) delete fDebugStreamer;
489 void AliTPCtrackerMI::FillESD(const TObjArray* arr)
493 //fill esds using updated tracks
496 // write tracks to the event
497 // store index of the track
498 Int_t nseed=arr->GetEntriesFast();
499 //FindKinks(arr,fEvent);
500 for (Int_t i=0; i<nseed; i++) {
501 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
505 if (AliTPCReconstructor::StreamLevel()>1) {
506 (*fDebugStreamer)<<"Track0"<<
510 // pt->PropagateTo(fkParam->GetInnerRadiusLow());
511 if (pt->GetKinkIndex(0)<=0){ //don't propagate daughter tracks
512 pt->PropagateTo(fkParam->GetInnerRadiusLow());
515 if (( pt->GetPoints()[2]- pt->GetPoints()[0])>5 && pt->GetPoints()[3]>0.8){
517 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
518 iotrack.SetTPCPoints(pt->GetPoints());
519 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
520 iotrack.SetV0Indexes(pt->GetV0Indexes());
521 // iotrack.SetTPCpid(pt->fTPCr);
522 //iotrack.SetTPCindex(i);
523 fEvent->AddTrack(&iotrack);
527 if ( (pt->GetNumberOfClusters()>70)&& (Float_t(pt->GetNumberOfClusters())/Float_t(pt->GetNFoundable()))>0.55) {
529 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
530 iotrack.SetTPCPoints(pt->GetPoints());
531 //iotrack.SetTPCindex(i);
532 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
533 iotrack.SetV0Indexes(pt->GetV0Indexes());
534 // iotrack.SetTPCpid(pt->fTPCr);
535 fEvent->AddTrack(&iotrack);
539 // short tracks - maybe decays
541 if ( (pt->GetNumberOfClusters()>30) && (Float_t(pt->GetNumberOfClusters())/Float_t(pt->GetNFoundable()))>0.70) {
542 Int_t found,foundable,shared;
543 pt->GetClusterStatistic(0,60,found, foundable,shared,kFALSE);
544 if ( (found>20) && (pt->GetNShared()/float(pt->GetNumberOfClusters())<0.2)){
546 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
547 //iotrack.SetTPCindex(i);
548 iotrack.SetTPCPoints(pt->GetPoints());
549 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
550 iotrack.SetV0Indexes(pt->GetV0Indexes());
551 //iotrack.SetTPCpid(pt->fTPCr);
552 fEvent->AddTrack(&iotrack);
557 if ( (pt->GetNumberOfClusters()>20) && (Float_t(pt->GetNumberOfClusters())/Float_t(pt->GetNFoundable()))>0.8) {
558 Int_t found,foundable,shared;
559 pt->GetClusterStatistic(0,60,found, foundable,shared,kFALSE);
560 if (found<20) continue;
561 if (pt->GetNShared()/float(pt->GetNumberOfClusters())>0.2) continue;
564 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
565 iotrack.SetTPCPoints(pt->GetPoints());
566 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
567 iotrack.SetV0Indexes(pt->GetV0Indexes());
568 //iotrack.SetTPCpid(pt->fTPCr);
569 //iotrack.SetTPCindex(i);
570 fEvent->AddTrack(&iotrack);
573 // short tracks - secondaties
575 if ( (pt->GetNumberOfClusters()>30) ) {
576 Int_t found,foundable,shared;
577 pt->GetClusterStatistic(128,158,found, foundable,shared,kFALSE);
578 if ( (found>20) && (pt->GetNShared()/float(pt->GetNumberOfClusters())<0.2) &&float(found)/float(foundable)>0.8){
580 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
581 iotrack.SetTPCPoints(pt->GetPoints());
582 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
583 iotrack.SetV0Indexes(pt->GetV0Indexes());
584 //iotrack.SetTPCpid(pt->fTPCr);
585 //iotrack.SetTPCindex(i);
586 fEvent->AddTrack(&iotrack);
591 if ( (pt->GetNumberOfClusters()>15)) {
592 Int_t found,foundable,shared;
593 pt->GetClusterStatistic(138,158,found, foundable,shared,kFALSE);
594 if (found<15) continue;
595 if (foundable<=0) continue;
596 if (pt->GetNShared()/float(pt->GetNumberOfClusters())>0.2) continue;
597 if (float(found)/float(foundable)<0.8) continue;
600 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
601 iotrack.SetTPCPoints(pt->GetPoints());
602 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
603 iotrack.SetV0Indexes(pt->GetV0Indexes());
604 // iotrack.SetTPCpid(pt->fTPCr);
605 //iotrack.SetTPCindex(i);
606 fEvent->AddTrack(&iotrack);
610 // >> account for suppressed tracks in the kink indices (RS)
611 int nESDtracks = fEvent->GetNumberOfTracks();
612 for (int it=nESDtracks;it--;) {
613 AliESDtrack* esdTr = fEvent->GetTrack(it);
614 if (!esdTr || !esdTr->GetKinkIndex(0)) continue;
615 for (int ik=0;ik<3;ik++) {
617 if (!(knkId=esdTr->GetKinkIndex(ik))) break; // no more kinks for this track
618 AliESDkink* kink = fEvent->GetKink(TMath::Abs(knkId)-1);
620 AliError(Form("ESDTrack%d refers to non-existing kink %d",it,TMath::Abs(knkId)-1));
623 kink->SetIndex(it, knkId<0 ? 0:1); // update track index of the kink: mother at 0, daughter at 1
626 // << account for suppressed tracks in the kink indices (RS)
627 AliInfo(Form("Number of filled ESDs-\t%d\n",fEvent->GetNumberOfTracks()));
635 Double_t AliTPCtrackerMI::ErrY2(AliTPCseed* seed, const AliTPCclusterMI * cl){
638 // Use calibrated cluster error from OCDB
640 AliTPCClusterParam * clparam = AliTPCcalibDB::Instance()->GetClusterParam();
642 Float_t z = TMath::Abs(fkParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
643 Int_t ctype = cl->GetType();
644 Int_t type = (cl->GetRow()<63) ? 0: (cl->GetRow()>126) ? 1:2;
645 Double_t angle = seed->GetSnp()*seed->GetSnp();
646 angle = TMath::Sqrt(TMath::Abs(angle/(1.-angle)));
647 Double_t erry2 = clparam->GetError0Par(0,type, z,angle);
649 erry2+=0.5; // edge cluster
652 seed->SetErrorY2(erry2);
656 //calculate look-up table at the beginning
657 // static Bool_t ginit = kFALSE;
658 // static Float_t gnoise1,gnoise2,gnoise3;
659 // static Float_t ggg1[10000];
660 // static Float_t ggg2[10000];
661 // static Float_t ggg3[10000];
662 // static Float_t glandau1[10000];
663 // static Float_t glandau2[10000];
664 // static Float_t glandau3[10000];
666 // static Float_t gcor01[500];
667 // static Float_t gcor02[500];
668 // static Float_t gcorp[500];
672 // if (ginit==kFALSE){
673 // for (Int_t i=1;i<500;i++){
674 // Float_t rsigma = float(i)/100.;
675 // gcor02[i] = TMath::Max(0.78 +TMath::Exp(7.4*(rsigma-1.2)),0.6);
676 // gcor01[i] = TMath::Max(0.72 +TMath::Exp(3.36*(rsigma-1.2)),0.6);
677 // gcorp[i] = TMath::Max(TMath::Power((rsigma+0.5),1.5),1.2);
681 // for (Int_t i=3;i<10000;i++){
685 // Float_t amp = float(i);
686 // Float_t padlength =0.75;
687 // gnoise1 = 0.0004/padlength;
688 // Float_t nel = 0.268*amp;
689 // Float_t nprim = 0.155*amp;
690 // ggg1[i] = fkParam->GetDiffT()*fkParam->GetDiffT()*(2+0.001*nel/(padlength*padlength))/nel;
691 // glandau1[i] = (2.+0.12*nprim)*0.5* (2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
692 // if (glandau1[i]>1) glandau1[i]=1;
693 // glandau1[i]*=padlength*padlength/12.;
697 // gnoise2 = 0.0004/padlength;
699 // nprim = 0.133*amp;
700 // ggg2[i] = fkParam->GetDiffT()*fkParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
701 // glandau2[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
702 // if (glandau2[i]>1) glandau2[i]=1;
703 // glandau2[i]*=padlength*padlength/12.;
708 // gnoise3 = 0.0004/padlength;
710 // nprim = 0.133*amp;
711 // ggg3[i] = fkParam->GetDiffT()*fkParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
712 // glandau3[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
713 // if (glandau3[i]>1) glandau3[i]=1;
714 // glandau3[i]*=padlength*padlength/12.;
722 // Int_t amp = int(TMath::Abs(cl->GetQ()));
724 // seed->SetErrorY2(1.);
728 // Float_t z = TMath::Abs(fkParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
729 // Int_t ctype = cl->GetType();
730 // Float_t padlength= GetPadPitchLength(seed->GetRow());
731 // Double_t angle2 = seed->GetSnp()*seed->GetSnp();
732 // angle2 = angle2/(1-angle2);
734 // //cluster "quality"
735 // Int_t rsigmay = int(100.*cl->GetSigmaY2()/(seed->GetCurrentSigmaY2()));
738 // if (fSectors==fInnerSec){
739 // snoise2 = gnoise1;
740 // res = ggg1[amp]*z+glandau1[amp]*angle2;
741 // if (ctype==0) res *= gcor01[rsigmay];
744 // res*= gcorp[rsigmay];
748 // if (padlength<1.1){
749 // snoise2 = gnoise2;
750 // res = ggg2[amp]*z+glandau2[amp]*angle2;
751 // if (ctype==0) res *= gcor02[rsigmay];
754 // res*= gcorp[rsigmay];
758 // snoise2 = gnoise3;
759 // res = ggg3[amp]*z+glandau3[amp]*angle2;
760 // if (ctype==0) res *= gcor02[rsigmay];
763 // res*= gcorp[rsigmay];
770 // res*=2.4; // overestimate error 2 times
774 // if (res<2*snoise2)
777 // seed->SetErrorY2(res);
785 Double_t AliTPCtrackerMI::ErrZ2(AliTPCseed* seed, const AliTPCclusterMI * cl){
788 // Use calibrated cluster error from OCDB
790 AliTPCClusterParam * clparam = AliTPCcalibDB::Instance()->GetClusterParam();
792 Float_t z = TMath::Abs(fkParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
793 Int_t ctype = cl->GetType();
794 Int_t type = (cl->GetRow()<63) ? 0: (cl->GetRow()>126) ? 1:2;
796 Double_t angle2 = seed->GetSnp()*seed->GetSnp();
797 angle2 = seed->GetTgl()*seed->GetTgl()*(1+angle2/(1-angle2));
798 Double_t angle = TMath::Sqrt(TMath::Abs(angle2));
799 Double_t errz2 = clparam->GetError0Par(1,type, z,angle);
801 errz2+=0.5; // edge cluster
804 seed->SetErrorZ2(errz2);
810 // //seed->SetErrorY2(0.1);
812 // //calculate look-up table at the beginning
813 // static Bool_t ginit = kFALSE;
814 // static Float_t gnoise1,gnoise2,gnoise3;
815 // static Float_t ggg1[10000];
816 // static Float_t ggg2[10000];
817 // static Float_t ggg3[10000];
818 // static Float_t glandau1[10000];
819 // static Float_t glandau2[10000];
820 // static Float_t glandau3[10000];
822 // static Float_t gcor01[1000];
823 // static Float_t gcor02[1000];
824 // static Float_t gcorp[1000];
828 // if (ginit==kFALSE){
829 // for (Int_t i=1;i<1000;i++){
830 // Float_t rsigma = float(i)/100.;
831 // gcor02[i] = TMath::Max(0.81 +TMath::Exp(6.8*(rsigma-1.2)),0.6);
832 // gcor01[i] = TMath::Max(0.72 +TMath::Exp(2.04*(rsigma-1.2)),0.6);
833 // gcorp[i] = TMath::Max(TMath::Power((rsigma+0.5),1.5),1.2);
837 // for (Int_t i=3;i<10000;i++){
841 // Float_t amp = float(i);
842 // Float_t padlength =0.75;
843 // gnoise1 = 0.0004/padlength;
844 // Float_t nel = 0.268*amp;
845 // Float_t nprim = 0.155*amp;
846 // ggg1[i] = fkParam->GetDiffT()*fkParam->GetDiffT()*(2+0.001*nel/(padlength*padlength))/nel;
847 // glandau1[i] = (2.+0.12*nprim)*0.5* (2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
848 // if (glandau1[i]>1) glandau1[i]=1;
849 // glandau1[i]*=padlength*padlength/12.;
853 // gnoise2 = 0.0004/padlength;
855 // nprim = 0.133*amp;
856 // ggg2[i] = fkParam->GetDiffT()*fkParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
857 // glandau2[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
858 // if (glandau2[i]>1) glandau2[i]=1;
859 // glandau2[i]*=padlength*padlength/12.;
864 // gnoise3 = 0.0004/padlength;
866 // nprim = 0.133*amp;
867 // ggg3[i] = fkParam->GetDiffT()*fkParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
868 // glandau3[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
869 // if (glandau3[i]>1) glandau3[i]=1;
870 // glandau3[i]*=padlength*padlength/12.;
878 // Int_t amp = int(TMath::Abs(cl->GetQ()));
880 // seed->SetErrorY2(1.);
884 // Float_t z = TMath::Abs(fkParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
885 // Int_t ctype = cl->GetType();
886 // Float_t padlength= GetPadPitchLength(seed->GetRow());
888 // Double_t angle2 = seed->GetSnp()*seed->GetSnp();
889 // // if (angle2<0.6) angle2 = 0.6;
890 // angle2 = seed->GetTgl()*seed->GetTgl()*(1+angle2/(1-angle2));
892 // //cluster "quality"
893 // Int_t rsigmaz = int(100.*cl->GetSigmaZ2()/(seed->GetCurrentSigmaZ2()));
896 // if (fSectors==fInnerSec){
897 // snoise2 = gnoise1;
898 // res = ggg1[amp]*z+glandau1[amp]*angle2;
899 // if (ctype==0) res *= gcor01[rsigmaz];
902 // res*= gcorp[rsigmaz];
906 // if (padlength<1.1){
907 // snoise2 = gnoise2;
908 // res = ggg2[amp]*z+glandau2[amp]*angle2;
909 // if (ctype==0) res *= gcor02[rsigmaz];
912 // res*= gcorp[rsigmaz];
916 // snoise2 = gnoise3;
917 // res = ggg3[amp]*z+glandau3[amp]*angle2;
918 // if (ctype==0) res *= gcor02[rsigmaz];
921 // res*= gcorp[rsigmaz];
930 // if ((ctype<0) &&<70){
935 // if (res<2*snoise2)
937 // if (res>3) res =3;
938 // seed->SetErrorZ2(res);
946 void AliTPCtrackerMI::RotateToLocal(AliTPCseed *seed)
948 //rotate to track "local coordinata
949 Float_t x = seed->GetX();
950 Float_t y = seed->GetY();
951 Float_t ymax = x*TMath::Tan(0.5*fSectors->GetAlpha());
954 seed->SetRelativeSector((seed->GetRelativeSector()+1) % fN);
955 if (!seed->Rotate(fSectors->GetAlpha()))
957 } else if (y <-ymax) {
958 seed->SetRelativeSector((seed->GetRelativeSector()-1+fN) % fN);
959 if (!seed->Rotate(-fSectors->GetAlpha()))
967 //_____________________________________________________________________________
968 Double_t AliTPCtrackerMI::F1old(Double_t x1,Double_t y1,
969 Double_t x2,Double_t y2,
970 Double_t x3,Double_t y3) const
972 //-----------------------------------------------------------------
973 // Initial approximation of the track curvature
974 //-----------------------------------------------------------------
975 Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
976 Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
977 (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
978 Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
979 (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
981 Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
982 if ( xr*xr+yr*yr<=0.00000000000001) return 100;
983 return -xr*yr/sqrt(xr*xr+yr*yr);
988 //_____________________________________________________________________________
989 Double_t AliTPCtrackerMI::F1(Double_t x1,Double_t y1,
990 Double_t x2,Double_t y2,
991 Double_t x3,Double_t y3) const
993 //-----------------------------------------------------------------
994 // Initial approximation of the track curvature
995 //-----------------------------------------------------------------
1001 Double_t det = x3*y2-x2*y3;
1002 if (TMath::Abs(det)<1e-10){
1006 Double_t u = 0.5* (x2*(x2-x3)+y2*(y2-y3))/det;
1007 Double_t x0 = x3*0.5-y3*u;
1008 Double_t y0 = y3*0.5+x3*u;
1009 Double_t c2 = 1/TMath::Sqrt(x0*x0+y0*y0);
1015 Double_t AliTPCtrackerMI::F2(Double_t x1,Double_t y1,
1016 Double_t x2,Double_t y2,
1017 Double_t x3,Double_t y3) const
1019 //-----------------------------------------------------------------
1020 // Initial approximation of the track curvature
1021 //-----------------------------------------------------------------
1027 Double_t det = x3*y2-x2*y3;
1028 if (TMath::Abs(det)<1e-10) {
1032 Double_t u = 0.5* (x2*(x2-x3)+y2*(y2-y3))/det;
1033 Double_t x0 = x3*0.5-y3*u;
1034 Double_t y0 = y3*0.5+x3*u;
1035 Double_t c2 = 1/TMath::Sqrt(x0*x0+y0*y0);
1044 //_____________________________________________________________________________
1045 Double_t AliTPCtrackerMI::F2old(Double_t x1,Double_t y1,
1046 Double_t x2,Double_t y2,
1047 Double_t x3,Double_t y3) const
1049 //-----------------------------------------------------------------
1050 // Initial approximation of the track curvature times center of curvature
1051 //-----------------------------------------------------------------
1052 Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
1053 Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
1054 (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
1055 Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
1056 (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
1058 Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
1060 return -a/(d*y1-b)*xr/sqrt(xr*xr+yr*yr);
1063 //_____________________________________________________________________________
1064 Double_t AliTPCtrackerMI::F3(Double_t x1,Double_t y1,
1065 Double_t x2,Double_t y2,
1066 Double_t z1,Double_t z2) const
1068 //-----------------------------------------------------------------
1069 // Initial approximation of the tangent of the track dip angle
1070 //-----------------------------------------------------------------
1071 return (z1 - z2)/sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
1075 Double_t AliTPCtrackerMI::F3n(Double_t x1,Double_t y1,
1076 Double_t x2,Double_t y2,
1077 Double_t z1,Double_t z2, Double_t c) const
1079 //-----------------------------------------------------------------
1080 // Initial approximation of the tangent of the track dip angle
1081 //-----------------------------------------------------------------
1085 //angle1 = (z1-z2)*c/(TMath::ASin(c*x1-ni)-TMath::ASin(c*x2-ni));
1087 Double_t d = TMath::Sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
1088 if (TMath::Abs(d*c*0.5)>1) return 0;
1089 // Double_t angle2 = TMath::ASin(d*c*0.5);
1090 // Double_t angle2 = AliTPCFastMath::FastAsin(d*c*0.5);
1091 Double_t angle2 = (d*c*0.5>0.1)? TMath::ASin(d*c*0.5): AliTPCFastMath::FastAsin(d*c*0.5);
1093 angle2 = (z1-z2)*c/(angle2*2.);
1097 Bool_t AliTPCtrackerMI::GetProlongation(Double_t x1, Double_t x2, Double_t x[5], Double_t &y, Double_t &z) const
1098 {//-----------------------------------------------------------------
1099 // This function find proloncation of a track to a reference plane x=x2.
1100 //-----------------------------------------------------------------
1104 if (TMath::Abs(x[4]*x1 - x[2]) >= 0.999) {
1108 Double_t c1=x[4]*x1 - x[2], r1=TMath::Sqrt((1.-c1)*(1.+c1));
1109 Double_t c2=x[4]*x2 - x[2], r2=TMath::Sqrt((1.-c2)*(1.+c2));
1113 Double_t dy = dx*(c1+c2)/(r1+r2);
1116 Double_t delta = x[4]*dx*(c1+c2)/(c1*r2 + c2*r1);
1118 if (TMath::Abs(delta)>0.01){
1119 dz = x[3]*TMath::ASin(delta)/x[4];
1121 dz = x[3]*AliTPCFastMath::FastAsin(delta)/x[4];
1124 //dz = x[3]*AliTPCFastMath::FastAsin(delta)/x[4];
1132 Int_t AliTPCtrackerMI::LoadClusters (TTree *const tree)
1137 return LoadClusters();
1141 Int_t AliTPCtrackerMI::LoadClusters(const TObjArray *arr)
1144 // load clusters to the memory
1145 AliTPCClustersRow *clrow = new AliTPCClustersRow("AliTPCclusterMI");
1146 Int_t lower = arr->LowerBound();
1147 Int_t entries = arr->GetEntriesFast();
1149 for (Int_t i=lower; i<entries; i++) {
1150 clrow = (AliTPCClustersRow*) arr->At(i);
1151 if(!clrow) continue;
1152 if(!clrow->GetArray()) continue;
1156 fkParam->AdjustSectorRow(clrow->GetID(),sec,row);
1158 for (Int_t icl=0; icl<clrow->GetArray()->GetEntriesFast(); icl++){
1159 Transform((AliTPCclusterMI*)(clrow->GetArray()->At(icl)));
1162 if (clrow->GetArray()->GetEntriesFast()<=0) continue;
1163 AliTPCtrackerRow * tpcrow=0;
1166 tpcrow = &(fInnerSec[sec%fkNIS][row]);
1170 tpcrow = &(fOuterSec[(sec-fkNIS*2)%fkNOS][row]);
1171 left = (sec-fkNIS*2)/fkNOS;
1174 tpcrow->SetN1(clrow->GetArray()->GetEntriesFast());
1175 tpcrow->SetClusters1(new AliTPCclusterMI[tpcrow->GetN1()]);
1176 for (Int_t j=0;j<tpcrow->GetN1();++j)
1177 tpcrow->SetCluster1(j, *(AliTPCclusterMI*)(clrow->GetArray()->At(j)));
1180 tpcrow->SetN2(clrow->GetArray()->GetEntriesFast());
1181 tpcrow->SetClusters2(new AliTPCclusterMI[tpcrow->GetN2()]);
1182 for (Int_t j=0;j<tpcrow->GetN2();++j)
1183 tpcrow->SetCluster2(j,*(AliTPCclusterMI*)(clrow->GetArray()->At(j)));
1185 clrow->GetArray()->Clear("C");
1194 Int_t AliTPCtrackerMI::LoadClusters(const TClonesArray *arr)
1197 // load clusters to the memory from one
1200 AliTPCclusterMI *clust=0;
1201 Int_t count[72][96] = { {0} , {0} };
1203 // loop over clusters
1204 for (Int_t icl=0; icl<arr->GetEntriesFast(); icl++) {
1205 clust = (AliTPCclusterMI*)arr->At(icl);
1206 if(!clust) continue;
1207 //printf("cluster: det %d, row %d \n", clust->GetDetector(),clust->GetRow());
1209 // transform clusters
1212 // count clusters per pad row
1213 count[clust->GetDetector()][clust->GetRow()]++;
1216 // insert clusters to sectors
1217 for (Int_t icl=0; icl<arr->GetEntriesFast(); icl++) {
1218 clust = (AliTPCclusterMI*)arr->At(icl);
1219 if(!clust) continue;
1221 Int_t sec = clust->GetDetector();
1222 Int_t row = clust->GetRow();
1224 // filter overlapping pad rows needed by HLT
1225 if(sec<fkNIS*2) { //IROCs
1226 if(row == 30) continue;
1229 if(row == 27 || row == 76) continue;
1235 fInnerSec[sec%fkNIS].InsertCluster(clust, count[sec][row], fkParam);
1238 left = (sec-fkNIS*2)/fkNOS;
1239 fOuterSec[(sec-fkNIS*2)%fkNOS].InsertCluster(clust, count[sec][row], fkParam);
1243 // Load functions must be called behind LoadCluster(TClonesArray*)
1245 //LoadOuterSectors();
1246 //LoadInnerSectors();
1252 Int_t AliTPCtrackerMI::LoadClusters()
1255 // load clusters to the memory
1256 AliTPCClustersRow *clrow= new AliTPCClustersRow("AliTPCclusterMI");
1258 // TTree * tree = fClustersArray.GetTree();
1260 TTree * tree = fInput;
1261 TBranch * br = tree->GetBranch("Segment");
1262 br->SetAddress(&clrow);
1264 Int_t j=Int_t(tree->GetEntries());
1265 for (Int_t i=0; i<j; i++) {
1269 fkParam->AdjustSectorRow(clrow->GetID(),sec,row);
1270 for (Int_t icl=0; icl<clrow->GetArray()->GetEntriesFast(); icl++){
1271 Transform((AliTPCclusterMI*)(clrow->GetArray()->At(icl)));
1274 AliTPCtrackerRow * tpcrow=0;
1277 tpcrow = &(fInnerSec[sec%fkNIS][row]);
1281 tpcrow = &(fOuterSec[(sec-fkNIS*2)%fkNOS][row]);
1282 left = (sec-fkNIS*2)/fkNOS;
1285 tpcrow->SetN1(clrow->GetArray()->GetEntriesFast());
1286 tpcrow->SetClusters1(new AliTPCclusterMI[tpcrow->GetN1()]);
1287 for (Int_t k=0;k<tpcrow->GetN1();++k)
1288 tpcrow->SetCluster1(k, *(AliTPCclusterMI*)(clrow->GetArray()->At(k)));
1291 tpcrow->SetN2(clrow->GetArray()->GetEntriesFast());
1292 tpcrow->SetClusters2(new AliTPCclusterMI[tpcrow->GetN2()]);
1293 for (Int_t k=0;k<tpcrow->GetN2();++k)
1294 tpcrow->SetCluster2(k,*(AliTPCclusterMI*)(clrow->GetArray()->At(k)));
1305 void AliTPCtrackerMI::UnloadClusters()
1308 // unload clusters from the memory
1310 Int_t nrows = fOuterSec->GetNRows();
1311 for (Int_t sec = 0;sec<fkNOS;sec++)
1312 for (Int_t row = 0;row<nrows;row++){
1313 AliTPCtrackerRow* tpcrow = &(fOuterSec[sec%fkNOS][row]);
1315 // if (tpcrow->fClusters1) delete []tpcrow->fClusters1;
1316 // if (tpcrow->fClusters2) delete []tpcrow->fClusters2;
1318 tpcrow->ResetClusters();
1321 nrows = fInnerSec->GetNRows();
1322 for (Int_t sec = 0;sec<fkNIS;sec++)
1323 for (Int_t row = 0;row<nrows;row++){
1324 AliTPCtrackerRow* tpcrow = &(fInnerSec[sec%fkNIS][row]);
1326 // if (tpcrow->fClusters1) delete []tpcrow->fClusters1;
1327 //if (tpcrow->fClusters2) delete []tpcrow->fClusters2;
1329 tpcrow->ResetClusters();
1335 void AliTPCtrackerMI::FillClusterArray(TObjArray* array) const{
1337 // Filling cluster to the array - For visualization purposes
1340 nrows = fOuterSec->GetNRows();
1341 for (Int_t sec = 0;sec<fkNOS;sec++)
1342 for (Int_t row = 0;row<nrows;row++){
1343 AliTPCtrackerRow* tpcrow = &(fOuterSec[sec%fkNOS][row]);
1344 if (!tpcrow) continue;
1345 for (Int_t icl = 0;icl<tpcrow->GetN();icl++){
1346 array->AddLast((TObject*)((*tpcrow)[icl]));
1349 nrows = fInnerSec->GetNRows();
1350 for (Int_t sec = 0;sec<fkNIS;sec++)
1351 for (Int_t row = 0;row<nrows;row++){
1352 AliTPCtrackerRow* tpcrow = &(fInnerSec[sec%fkNIS][row]);
1353 if (!tpcrow) continue;
1354 for (Int_t icl = 0;icl<tpcrow->GetN();icl++){
1355 array->AddLast((TObject*)(*tpcrow)[icl]);
1361 void AliTPCtrackerMI::Transform(AliTPCclusterMI * cluster){
1365 AliTPCcalibDB * calibDB = AliTPCcalibDB::Instance();
1366 AliTPCTransform *transform = calibDB->GetTransform() ;
1368 AliFatal("Tranformations not in calibDB");
1371 transform->SetCurrentRecoParam((AliTPCRecoParam*)AliTPCReconstructor::GetRecoParam());
1372 Double_t x[3]={cluster->GetRow(),cluster->GetPad(),cluster->GetTimeBin()};
1373 Int_t i[1]={cluster->GetDetector()};
1374 transform->Transform(x,i,0,1);
1375 // if (cluster->GetDetector()%36>17){
1380 // in debug mode check the transformation
1382 if (AliTPCReconstructor::StreamLevel()>2) {
1384 cluster->GetGlobalXYZ(gx);
1385 Int_t event = (fEvent==NULL)? 0: fEvent->GetEventNumberInFile();
1386 TTreeSRedirector &cstream = *fDebugStreamer;
1387 cstream<<"Transform"<<
1398 cluster->SetX(x[0]);
1399 cluster->SetY(x[1]);
1400 cluster->SetZ(x[2]);
1405 //if (!fkParam->IsGeoRead()) fkParam->ReadGeoMatrices();
1406 if (AliTPCReconstructor::GetRecoParam()->GetUseSectorAlignment() && (!calibDB->HasAlignmentOCDB())){
1407 TGeoHMatrix *mat = fkParam->GetClusterMatrix(cluster->GetDetector());
1409 Double_t pos[3]= {cluster->GetX(),cluster->GetY(),cluster->GetZ()};
1410 Double_t posC[3]={cluster->GetX(),cluster->GetY(),cluster->GetZ()};
1411 if (mat) mat->LocalToMaster(pos,posC);
1413 // chack Loading of Geo matrices from GeoManager - TEMPORARY FIX
1415 cluster->SetX(posC[0]);
1416 cluster->SetY(posC[1]);
1417 cluster->SetZ(posC[2]);
1421 //_____________________________________________________________________________
1422 Int_t AliTPCtrackerMI::LoadOuterSectors() {
1423 //-----------------------------------------------------------------
1424 // This function fills outer TPC sectors with clusters.
1425 //-----------------------------------------------------------------
1426 Int_t nrows = fOuterSec->GetNRows();
1428 for (Int_t sec = 0;sec<fkNOS;sec++)
1429 for (Int_t row = 0;row<nrows;row++){
1430 AliTPCtrackerRow* tpcrow = &(fOuterSec[sec%fkNOS][row]);
1431 Int_t sec2 = sec+2*fkNIS;
1433 Int_t ncl = tpcrow->GetN1();
1435 AliTPCclusterMI *c= (tpcrow->GetCluster1(ncl));
1436 index=(((sec2<<8)+row)<<16)+ncl;
1437 tpcrow->InsertCluster(c,index);
1440 ncl = tpcrow->GetN2();
1442 AliTPCclusterMI *c= (tpcrow->GetCluster2(ncl));
1443 index=((((sec2+fkNOS)<<8)+row)<<16)+ncl;
1444 tpcrow->InsertCluster(c,index);
1447 // write indexes for fast acces
1449 for (Int_t i=0;i<510;i++)
1450 tpcrow->SetFastCluster(i,-1);
1451 for (Int_t i=0;i<tpcrow->GetN();i++){
1452 Int_t zi = Int_t((*tpcrow)[i]->GetZ()+255.);
1453 tpcrow->SetFastCluster(zi,i); // write index
1456 for (Int_t i=0;i<510;i++){
1457 if (tpcrow->GetFastCluster(i)<0)
1458 tpcrow->SetFastCluster(i,last);
1460 last = tpcrow->GetFastCluster(i);
1469 //_____________________________________________________________________________
1470 Int_t AliTPCtrackerMI::LoadInnerSectors() {
1471 //-----------------------------------------------------------------
1472 // This function fills inner TPC sectors with clusters.
1473 //-----------------------------------------------------------------
1474 Int_t nrows = fInnerSec->GetNRows();
1476 for (Int_t sec = 0;sec<fkNIS;sec++)
1477 for (Int_t row = 0;row<nrows;row++){
1478 AliTPCtrackerRow* tpcrow = &(fInnerSec[sec%fkNIS][row]);
1481 Int_t ncl = tpcrow->GetN1();
1483 AliTPCclusterMI *c= (tpcrow->GetCluster1(ncl));
1484 index=(((sec<<8)+row)<<16)+ncl;
1485 tpcrow->InsertCluster(c,index);
1488 ncl = tpcrow->GetN2();
1490 AliTPCclusterMI *c= (tpcrow->GetCluster2(ncl));
1491 index=((((sec+fkNIS)<<8)+row)<<16)+ncl;
1492 tpcrow->InsertCluster(c,index);
1495 // write indexes for fast acces
1497 for (Int_t i=0;i<510;i++)
1498 tpcrow->SetFastCluster(i,-1);
1499 for (Int_t i=0;i<tpcrow->GetN();i++){
1500 Int_t zi = Int_t((*tpcrow)[i]->GetZ()+255.);
1501 tpcrow->SetFastCluster(zi,i); // write index
1504 for (Int_t i=0;i<510;i++){
1505 if (tpcrow->GetFastCluster(i)<0)
1506 tpcrow->SetFastCluster(i,last);
1508 last = tpcrow->GetFastCluster(i);
1520 //_________________________________________________________________________
1521 AliTPCclusterMI *AliTPCtrackerMI::GetClusterMI(Int_t index) const {
1522 //--------------------------------------------------------------------
1523 // Return pointer to a given cluster
1524 //--------------------------------------------------------------------
1525 if (index<0) return 0; // no cluster
1526 Int_t sec=(index&0xff000000)>>24;
1527 Int_t row=(index&0x00ff0000)>>16;
1528 Int_t ncl=(index&0x00007fff)>>00;
1530 const AliTPCtrackerRow * tpcrow=0;
1531 AliTPCclusterMI * clrow =0;
1533 if (sec<0 || sec>=fkNIS*4) {
1534 AliWarning(Form("Wrong sector %d",sec));
1539 AliTPCtrackerSector& tracksec = fInnerSec[sec%fkNIS];
1540 if (tracksec.GetNRows()<=row) return 0;
1541 tpcrow = &(tracksec[row]);
1542 if (tpcrow==0) return 0;
1545 if (tpcrow->GetN1()<=ncl) return 0;
1546 clrow = tpcrow->GetClusters1();
1549 if (tpcrow->GetN2()<=ncl) return 0;
1550 clrow = tpcrow->GetClusters2();
1554 AliTPCtrackerSector& tracksec = fOuterSec[(sec-fkNIS*2)%fkNOS];
1555 if (tracksec.GetNRows()<=row) return 0;
1556 tpcrow = &(tracksec[row]);
1557 if (tpcrow==0) return 0;
1559 if (sec-2*fkNIS<fkNOS) {
1560 if (tpcrow->GetN1()<=ncl) return 0;
1561 clrow = tpcrow->GetClusters1();
1564 if (tpcrow->GetN2()<=ncl) return 0;
1565 clrow = tpcrow->GetClusters2();
1569 return &(clrow[ncl]);
1575 Int_t AliTPCtrackerMI::FollowToNext(AliTPCseed& t, Int_t nr) {
1576 //-----------------------------------------------------------------
1577 // This function tries to find a track prolongation to next pad row
1578 //-----------------------------------------------------------------
1580 Double_t x= GetXrow(nr), ymax=GetMaxY(nr);
1583 AliTPCclusterMI *cl=0;
1584 Int_t tpcindex= t.GetClusterIndex2(nr);
1586 // update current shape info every 5 pad-row
1587 // if ( (nr%5==0) || t.GetNumberOfClusters()<2 || (t.fCurrentSigmaY2<0.0001) ){
1591 if (fIteration>0 && tpcindex>=-1){ //if we have already clusters
1593 if (tpcindex==-1) return 0; //track in dead zone
1595 cl = t.GetClusterPointer(nr);
1596 if ( (cl==0) ) cl = GetClusterMI(tpcindex);
1597 t.SetCurrentClusterIndex1(tpcindex);
1600 Int_t relativesector = ((tpcindex&0xff000000)>>24)%18; // if previously accepted cluster in different sector
1601 Float_t angle = relativesector*fSectors->GetAlpha()+fSectors->GetAlphaShift();
1603 if (angle<-TMath::Pi()) angle += 2*TMath::Pi();
1604 if (angle>=TMath::Pi()) angle -= 2*TMath::Pi();
1606 if (TMath::Abs(angle-t.GetAlpha())>0.001){
1607 Double_t rotation = angle-t.GetAlpha();
1608 t.SetRelativeSector(relativesector);
1609 if (!t.Rotate(rotation)) return 0;
1611 if (!t.PropagateTo(x)) return 0;
1613 t.SetCurrentCluster(cl);
1615 Int_t accept = AcceptCluster(&t,t.GetCurrentCluster());
1616 if ((tpcindex&0x8000)==0) accept =0;
1618 //if founded cluster is acceptible
1619 if (cl->IsUsed(11)) { // id cluster is shared inrease uncertainty
1620 t.SetErrorY2(t.GetErrorY2()+0.03);
1621 t.SetErrorZ2(t.GetErrorZ2()+0.03);
1622 t.SetErrorY2(t.GetErrorY2()*3);
1623 t.SetErrorZ2(t.GetErrorZ2()*3);
1625 t.SetNFoundable(t.GetNFoundable()+1);
1626 UpdateTrack(&t,accept);
1631 if (TMath::Abs(t.GetSnp())>AliTPCReconstructor::GetMaxSnpTracker()) return 0; // cut on angle
1632 if (fIteration>1 && IsFindable(t)){
1633 // not look for new cluster during refitting
1634 t.SetNFoundable(t.GetNFoundable()+1);
1639 // if (TMath::Abs(t.GetSnp())>0.95 || TMath::Abs(x*t.GetC()-t.GetEta())>0.95) return 0;// patch 28 fev 06
1640 Double_t y=t.GetYat(x);
1641 if (TMath::Abs(y)>ymax){
1643 t.SetRelativeSector((t.GetRelativeSector()+1) % fN);
1644 if (!t.Rotate(fSectors->GetAlpha()))
1646 } else if (y <-ymax) {
1647 t.SetRelativeSector((t.GetRelativeSector()-1+fN) % fN);
1648 if (!t.Rotate(-fSectors->GetAlpha()))
1654 if (!t.PropagateTo(x)) {
1655 if (fIteration==0) t.SetRemoval(10);
1659 Double_t z=t.GetZ();
1662 if (!IsActive(t.GetRelativeSector(),nr)) {
1664 t.SetClusterIndex2(nr,-1);
1667 //AliInfo(Form("A - Sector%d phi %f - alpha %f", t.fRelativeSector,y/x, t.GetAlpha()));
1668 Bool_t isActive = IsActive(t.GetRelativeSector(),nr);
1669 Bool_t isActive2 = (nr>=fInnerSec->GetNRows()) ? fOuterSec[t.GetRelativeSector()][nr-fInnerSec->GetNRows()].GetN()>0:fInnerSec[t.GetRelativeSector()][nr].GetN()>0;
1671 if (!isActive || !isActive2) return 0;
1673 const AliTPCtrackerRow &krow=GetRow(t.GetRelativeSector(),nr);
1674 if ( (t.GetSigmaY2()<0) || t.GetSigmaZ2()<0) return 0;
1676 Double_t roadz = 1.;
1678 if (TMath::Abs(TMath::Abs(y)-ymax)<krow.GetDeadZone()){
1680 t.SetClusterIndex2(nr,-1);
1686 // if (TMath::Abs(z)<(AliTPCReconstructor::GetCtgRange()*x+10) && TMath::Abs(z)<fkParam->GetZLength(0) && (TMath::Abs(t.GetSnp())<AliTPCReconstructor::GetMaxSnpTracker()))
1687 t.SetNFoundable(t.GetNFoundable()+1);
1693 // cl = krow.FindNearest2(y+10.,z,roady,roadz,index);
1694 cl = krow.FindNearest2(y,z,roady,roadz,index);
1695 if (cl) t.SetCurrentClusterIndex1(krow.GetIndex(index));
1698 t.SetCurrentCluster(cl);
1700 if (fIteration==2&&cl->IsUsed(10)) return 0;
1701 Int_t accept = AcceptCluster(&t,t.GetCurrentCluster());
1702 if (fIteration==2&&cl->IsUsed(11)) {
1703 t.SetErrorY2(t.GetErrorY2()+0.03);
1704 t.SetErrorZ2(t.GetErrorZ2()+0.03);
1705 t.SetErrorY2(t.GetErrorY2()*3);
1706 t.SetErrorZ2(t.GetErrorZ2()*3);
1709 if (t.fCurrentCluster->IsUsed(10)){
1714 if (t.fNShared>0.7*t.GetNumberOfClusters()) {
1720 if (accept<3) UpdateTrack(&t,accept);
1723 if ( fIteration==0 && t.GetNFoundable()*0.5 > t.GetNumberOfClusters()) t.SetRemoval(10);
1731 //_________________________________________________________________________
1732 Bool_t AliTPCtrackerMI::GetTrackPoint(Int_t index, AliTrackPoint &p ) const
1734 // Get track space point by index
1735 // return false in case the cluster doesn't exist
1736 AliTPCclusterMI *cl = GetClusterMI(index);
1737 if (!cl) return kFALSE;
1738 Int_t sector = (index&0xff000000)>>24;
1739 // Int_t row = (index&0x00ff0000)>>16;
1741 // xyz[0] = fkParam->GetPadRowRadii(sector,row);
1742 xyz[0] = cl->GetX();
1743 xyz[1] = cl->GetY();
1744 xyz[2] = cl->GetZ();
1746 fkParam->AdjustCosSin(sector,cos,sin);
1747 Float_t x = cos*xyz[0]-sin*xyz[1];
1748 Float_t y = cos*xyz[1]+sin*xyz[0];
1750 Float_t sigmaY2 = 0.027*cl->GetSigmaY2();
1751 if (sector < fkParam->GetNInnerSector()) sigmaY2 *= 2.07;
1752 Float_t sigmaZ2 = 0.066*cl->GetSigmaZ2();
1753 if (sector < fkParam->GetNInnerSector()) sigmaZ2 *= 1.77;
1754 cov[0] = sin*sin*sigmaY2;
1755 cov[1] = -sin*cos*sigmaY2;
1757 cov[3] = cos*cos*sigmaY2;
1760 p.SetXYZ(x,y,xyz[2],cov);
1761 AliGeomManager::ELayerID iLayer;
1763 if (sector < fkParam->GetNInnerSector()) {
1764 iLayer = AliGeomManager::kTPC1;
1768 iLayer = AliGeomManager::kTPC2;
1769 idet = sector - fkParam->GetNInnerSector();
1771 UShort_t volid = AliGeomManager::LayerToVolUID(iLayer,idet);
1772 p.SetVolumeID(volid);
1778 Int_t AliTPCtrackerMI::UpdateClusters(AliTPCseed& t, Int_t nr) {
1779 //-----------------------------------------------------------------
1780 // This function tries to find a track prolongation to next pad row
1781 //-----------------------------------------------------------------
1782 t.SetCurrentCluster(0);
1783 t.SetCurrentClusterIndex1(0);
1785 Double_t xt=t.GetX();
1786 Int_t row = GetRowNumber(xt)-1;
1787 Double_t ymax= GetMaxY(nr);
1789 if (row < nr) return 1; // don't prolongate if not information until now -
1790 // if (TMath::Abs(t.GetSnp())>0.9 && t.GetNumberOfClusters()>40. && fIteration!=2) {
1792 // return 0; // not prolongate strongly inclined tracks
1794 // if (TMath::Abs(t.GetSnp())>0.95) {
1796 // return 0; // not prolongate strongly inclined tracks
1797 // }// patch 28 fev 06
1799 Double_t x= GetXrow(nr);
1801 //t.PropagateTo(x+0.02);
1802 //t.PropagateTo(x+0.01);
1803 if (!t.PropagateTo(x)){
1810 if (TMath::Abs(y)>ymax){
1812 t.SetRelativeSector((t.GetRelativeSector()+1) % fN);
1813 if (!t.Rotate(fSectors->GetAlpha()))
1815 } else if (y <-ymax) {
1816 t.SetRelativeSector((t.GetRelativeSector()-1+fN) % fN);
1817 if (!t.Rotate(-fSectors->GetAlpha()))
1820 // if (!t.PropagateTo(x)){
1827 if (TMath::Abs(t.GetSnp())>AliTPCReconstructor::GetMaxSnpTracker()) return 0;
1829 if (!IsActive(t.GetRelativeSector(),nr)) {
1831 t.SetClusterIndex2(nr,-1);
1834 //AliInfo(Form("A - Sector%d phi %f - alpha %f", t.fRelativeSector,y/x, t.GetAlpha()));
1836 AliTPCtrackerRow &krow=GetRow(t.GetRelativeSector(),nr);
1838 if (TMath::Abs(TMath::Abs(y)-ymax)<krow.GetDeadZone()){
1840 t.SetClusterIndex2(nr,-1);
1846 // if (TMath::Abs(t.GetZ())<(AliTPCReconstructor::GetCtgRange()*t.GetX()+10) && (TMath::Abs(t.GetSnp())<AliTPCReconstructor::GetMaxSnpTracker()))
1847 if (IsFindable(t)) t.SetNFoundable(t.GetNFoundable()+1);
1853 if ( (nr%2==0) || t.GetNumberOfClusters()<2){
1854 // t.fCurrentSigmaY = GetSigmaY(&t);
1855 //t.fCurrentSigmaZ = GetSigmaZ(&t);
1859 AliTPCclusterMI *cl=0;
1862 Double_t roady = 1.;
1863 Double_t roadz = 1.;
1867 index = t.GetClusterIndex2(nr);
1868 if ( (index>0) && (index&0x8000)==0){
1869 cl = t.GetClusterPointer(nr);
1870 if ( (cl==0) && (index>0)) cl = GetClusterMI(index);
1871 t.SetCurrentClusterIndex1(index);
1873 t.SetCurrentCluster(cl);
1879 // if (index<0) return 0;
1880 UInt_t uindex = TMath::Abs(index);
1883 //cl = krow.FindNearest2(y+10,z,roady,roadz,uindex);
1884 cl = krow.FindNearest2(y,z,roady,roadz,uindex);
1887 if (cl) t.SetCurrentClusterIndex1(krow.GetIndex(uindex));
1888 t.SetCurrentCluster(cl);
1894 Int_t AliTPCtrackerMI::FollowToNextCluster(AliTPCseed & t, Int_t nr) {
1895 //-----------------------------------------------------------------
1896 // This function tries to find a track prolongation to next pad row
1897 //-----------------------------------------------------------------
1899 //update error according neighborhoud
1901 if (t.GetCurrentCluster()) {
1903 Int_t accept = AcceptCluster(&t,t.GetCurrentCluster());
1905 if (t.GetCurrentCluster()->IsUsed(10)){
1910 t.SetNShared(t.GetNShared()+1);
1911 if (t.GetNShared()>0.7*t.GetNumberOfClusters()) {
1916 if (fIteration>0) accept = 0;
1917 if (accept<3) UpdateTrack(&t,accept);
1921 if ( ( (t.GetSigmaY2()+t.GetSigmaZ2())>0.16)&& t.GetNumberOfClusters()>18) t.SetRemoval(10);
1922 if ( t.GetChi2()/t.GetNumberOfClusters()>6 &&t.GetNumberOfClusters()>18) t.SetRemoval(10);
1924 if (( (t.GetNFoundable()*0.5 > t.GetNumberOfClusters()) || t.GetNoCluster()>15)) t.SetRemoval(10);
1932 //_____________________________________________________________________________
1933 Int_t AliTPCtrackerMI::FollowProlongation(AliTPCseed& t, Int_t rf, Int_t step) {
1934 //-----------------------------------------------------------------
1935 // This function tries to find a track prolongation.
1936 //-----------------------------------------------------------------
1937 Double_t xt=t.GetX();
1939 Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
1940 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
1941 if (alpha < 0. ) alpha += 2.*TMath::Pi();
1943 t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN);
1945 Int_t first = GetRowNumber(xt)-1;
1946 for (Int_t nr= first; nr>=rf; nr-=step) {
1948 if (t.GetKinkIndexes()[0]>0){
1949 for (Int_t i=0;i<3;i++){
1950 Int_t index = t.GetKinkIndexes()[i];
1951 if (index==0) break;
1952 if (index<0) continue;
1954 AliKink * kink = (AliKink*)fEvent->GetKink(index-1);
1956 printf("PROBLEM\n");
1959 Int_t kinkrow = kink->GetTPCRow0()+2+Int_t(0.5/(0.05+kink->GetAngle(2)));
1961 AliExternalTrackParam paramd(t);
1962 kink->SetDaughter(paramd);
1963 kink->SetStatus(2,5);
1970 if (nr==80) t.UpdateReference();
1971 if (nr<fInnerSec->GetNRows())
1972 fSectors = fInnerSec;
1974 fSectors = fOuterSec;
1975 if (FollowToNext(t,nr)==0)
1988 Int_t AliTPCtrackerMI::FollowBackProlongation(AliTPCseed& t, Int_t rf) {
1989 //-----------------------------------------------------------------
1990 // This function tries to find a track prolongation.
1991 //-----------------------------------------------------------------
1993 Double_t xt=t.GetX();
1994 Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
1995 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
1996 if (alpha < 0. ) alpha += 2.*TMath::Pi();
1997 t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN);
1999 Int_t first = t.GetFirstPoint();
2000 if (first<GetRowNumber(xt)+1) first = GetRowNumber(xt)+1;
2002 if (first<0) first=0;
2003 for (Int_t nr=first; nr<=rf; nr++) {
2004 // if ( (TMath::Abs(t.GetSnp())>0.95)) break;//patch 28 fev 06
2005 if (t.GetKinkIndexes()[0]<0){
2006 for (Int_t i=0;i<3;i++){
2007 Int_t index = t.GetKinkIndexes()[i];
2008 if (index==0) break;
2009 if (index>0) continue;
2010 index = TMath::Abs(index);
2011 AliKink * kink = (AliKink*)fEvent->GetKink(index-1);
2013 printf("PROBLEM\n");
2016 Int_t kinkrow = kink->GetTPCRow0()-2-Int_t(0.5/(0.05+kink->GetAngle(2)));
2018 AliExternalTrackParam paramm(t);
2019 kink->SetMother(paramm);
2020 kink->SetStatus(2,1);
2027 if (nr<fInnerSec->GetNRows())
2028 fSectors = fInnerSec;
2030 fSectors = fOuterSec;
2041 Float_t AliTPCtrackerMI::OverlapFactor(AliTPCseed * s1, AliTPCseed * s2, Int_t &sum1, Int_t & sum2)
2049 Float_t dz2 =(s1->GetZ() - s2->GetZ());
2052 Float_t dy2 =TMath::Abs((s1->GetY() - s2->GetY()));
2054 Float_t distance = TMath::Sqrt(dz2+dy2);
2055 if (distance>4.) return 0; // if there are far away - not overlap - to reduce combinatorics
2058 Int_t firstpoint = TMath::Min(s1->GetFirstPoint(),s2->GetFirstPoint());
2059 Int_t lastpoint = TMath::Max(s1->GetLastPoint(),s2->GetLastPoint());
2064 if (firstpoint>lastpoint) {
2065 firstpoint =lastpoint;
2070 for (Int_t i=firstpoint-1;i<lastpoint+1;i++){
2071 if (s1->GetClusterIndex2(i)>0) sum1++;
2072 if (s2->GetClusterIndex2(i)>0) sum2++;
2073 if (s1->GetClusterIndex2(i)==s2->GetClusterIndex2(i) && s1->GetClusterIndex2(i)>0) {
2077 if (sum<5) return 0;
2079 Float_t summin = TMath::Min(sum1+1,sum2+1);
2080 Float_t ratio = (sum+1)/Float_t(summin);
2084 void AliTPCtrackerMI::SignShared(AliTPCseed * s1, AliTPCseed * s2)
2088 Float_t thetaCut = 0.2;//+10.*TMath::Sqrt(s1->GetSigmaTglZ()+ s2->GetSigmaTglZ());
2089 if (TMath::Abs(s1->GetTgl()-s2->GetTgl())>thetaCut) return;
2090 Float_t minCl = TMath::Min(s1->GetNumberOfClusters(),s2->GetNumberOfClusters());
2091 Int_t cutN0 = TMath::Max(5,TMath::Nint(0.1*minCl));
2096 //Int_t firstpoint = TMath::Max(s1->GetFirstPoint(),s2->GetFirstPoint());
2097 //Int_t lastpoint = TMath::Min(s1->GetLastPoint(),s2->GetLastPoint());
2098 Int_t firstpoint = 0;
2099 Int_t lastpoint = 160;
2101 // if (firstpoint>=lastpoint-5) return;;
2103 for (Int_t i=firstpoint;i<lastpoint;i++){
2104 // if ( (s1->GetClusterIndex2(i)&0xFFFF8FFF)==(s2->GetClusterIndex2(i)&0xFFFF8FFF) && s1->GetClusterIndex2(i)>0) {
2105 if ( (s1->GetClusterIndex2(i))==(s2->GetClusterIndex2(i)) && s1->GetClusterIndex2(i)>0) {
2107 s1->SetSharedMapBit(i, kTRUE);
2108 s2->SetSharedMapBit(i, kTRUE);
2110 if (s1->GetClusterIndex2(i)>0)
2111 s1->SetClusterMapBit(i, kTRUE);
2113 if (sumshared>cutN0){
2116 for (Int_t i=firstpoint;i<lastpoint;i++){
2117 // if ( (s1->GetClusterIndex2(i)&0xFFFF8FFF)==(s2->GetClusterIndex2(i)&0xFFFF8FFF) && s1->GetClusterIndex2(i)>0) {
2118 if ( (s1->GetClusterIndex2(i))==(s2->GetClusterIndex2(i)) && s1->GetClusterIndex2(i)>0) {
2119 AliTPCTrackerPoint *p1 = s1->GetTrackPoint(i);
2120 AliTPCTrackerPoint *p2 = s2->GetTrackPoint(i);;
2121 if (s1->IsActive()&&s2->IsActive()){
2122 p1->SetShared(kTRUE);
2123 p2->SetShared(kTRUE);
2129 if (sumshared>cutN0){
2130 for (Int_t i=0;i<4;i++){
2131 if (s1->GetOverlapLabel(3*i)==0){
2132 s1->SetOverlapLabel(3*i, s2->GetLabel());
2133 s1->SetOverlapLabel(3*i+1,sumshared);
2134 s1->SetOverlapLabel(3*i+2,s2->GetUniqueID());
2138 for (Int_t i=0;i<4;i++){
2139 if (s2->GetOverlapLabel(3*i)==0){
2140 s2->SetOverlapLabel(3*i, s1->GetLabel());
2141 s2->SetOverlapLabel(3*i+1,sumshared);
2142 s2->SetOverlapLabel(3*i+2,s1->GetUniqueID());
2149 void AliTPCtrackerMI::SignShared(TObjArray * arr)
2152 //sort trackss according sectors
2154 for (Int_t i=0; i<arr->GetEntriesFast(); i++) {
2155 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2157 //if (pt) RotateToLocal(pt);
2161 arr->Sort(); // sorting according relative sectors
2162 arr->Expand(arr->GetEntries());
2165 Int_t nseed=arr->GetEntriesFast();
2166 for (Int_t i=0; i<nseed; i++) {
2167 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2169 for (Int_t j=0;j<12;j++){
2170 pt->SetOverlapLabel(j,0);
2173 for (Int_t i=0; i<nseed; i++) {
2174 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2176 if (pt->GetRemoval()>10) continue;
2177 for (Int_t j=i+1; j<nseed; j++){
2178 AliTPCseed *pt2=(AliTPCseed*)arr->UncheckedAt(j);
2179 if (TMath::Abs(pt->GetRelativeSector()-pt2->GetRelativeSector())>1) continue;
2181 if (pt2->GetRemoval()<=10) {
2182 //if ( TMath::Abs(pt->GetRelativeSector()-pt2->GetRelativeSector())>0) break;
2190 void AliTPCtrackerMI::SortTracks(TObjArray * arr, Int_t mode) const
2193 //sort tracks in array according mode criteria
2194 Int_t nseed = arr->GetEntriesFast();
2195 for (Int_t i=0; i<nseed; i++) {
2196 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2207 void AliTPCtrackerMI::RemoveUsed2(TObjArray * arr, Float_t factor1, Float_t factor2, Int_t minimal)
2210 // Loop over all tracks and remove overlaped tracks (with lower quality)
2212 // 1. Unsign clusters
2213 // 2. Sort tracks according quality
2214 // Quality is defined by the number of cluster between first and last points
2216 // 3. Loop over tracks - decreasing quality order
2217 // a.) remove - If the fraction of shared cluster less than factor (1- n or 2)
2218 // b.) remove - If the minimal number of clusters less than minimal and not ITS
2219 // c.) if track accepted - sign clusters
2221 //Called in - AliTPCtrackerMI::Clusters2Tracks()
2222 // - AliTPCtrackerMI::PropagateBack()
2223 // - AliTPCtrackerMI::RefitInward()
2226 // factor1 - factor for constrained
2227 // factor2 - for non constrained tracks
2228 // if (Float_t(shared+1)/Float_t(found+1)>factor) - DELETE
2232 Int_t nseed = arr->GetEntriesFast();
2233 Float_t * quality = new Float_t[nseed];
2234 Int_t * indexes = new Int_t[nseed];
2238 for (Int_t i=0; i<nseed; i++) {
2239 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2244 pt->UpdatePoints(); //select first last max dens points
2245 Float_t * points = pt->GetPoints();
2246 if (points[3]<0.8) quality[i] =-1;
2247 quality[i] = (points[2]-points[0])+pt->GetNumberOfClusters();
2248 //prefer high momenta tracks if overlaps
2249 quality[i] *= TMath::Sqrt(TMath::Abs(pt->Pt())+0.5);
2251 TMath::Sort(nseed,quality,indexes);
2254 for (Int_t itrack=0; itrack<nseed; itrack++) {
2255 Int_t trackindex = indexes[itrack];
2256 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(trackindex);
2259 if (quality[trackindex]<0){
2260 delete arr->RemoveAt(trackindex);
2265 Int_t first = Int_t(pt->GetPoints()[0]);
2266 Int_t last = Int_t(pt->GetPoints()[2]);
2267 Double_t factor = (pt->GetBConstrain()) ? factor1: factor2;
2269 Int_t found,foundable,shared;
2270 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
2271 // pt->GetClusterStatistic(0,160, found, foundable,shared,kFALSE);
2272 Bool_t itsgold =kFALSE;
2275 if (pt->GetESD()->GetITSclusters(dummy)>4) itsgold= kTRUE;
2279 if (Float_t(shared+1)/Float_t(found+1)>factor){
2280 if (pt->GetKinkIndexes()[0]!=0) continue; //don't remove tracks - part of the kinks
2281 if( AliTPCReconstructor::StreamLevel()>3){
2282 TTreeSRedirector &cstream = *fDebugStreamer;
2283 cstream<<"RemoveUsed"<<
2284 "iter="<<fIteration<<
2288 delete arr->RemoveAt(trackindex);
2291 if (pt->GetNumberOfClusters()<50&&(found-0.5*shared)<minimal){ //remove short tracks
2292 if (pt->GetKinkIndexes()[0]!=0) continue; //don't remove tracks - part of the kinks
2293 if( AliTPCReconstructor::StreamLevel()>3){
2294 TTreeSRedirector &cstream = *fDebugStreamer;
2295 cstream<<"RemoveShort"<<
2296 "iter="<<fIteration<<
2300 delete arr->RemoveAt(trackindex);
2306 //if (sharedfactor>0.4) continue;
2307 if (pt->GetKinkIndexes()[0]>0) continue;
2308 //Remove tracks with undefined properties - seems
2309 if (pt->GetSigmaY2()<kAlmost0) continue; // ? what is the origin ?
2311 for (Int_t i=first; i<last; i++) {
2312 Int_t index=pt->GetClusterIndex2(i);
2313 // if (index<0 || index&0x8000 ) continue;
2314 if (index<0 || index&0x8000 ) continue;
2315 AliTPCclusterMI *c= pt->GetClusterPointer(i);
2322 Info("RemoveUsed","\n*****\nNumber of good tracks after shared removal\t%d\n",fNtracks);
2328 void AliTPCtrackerMI::DumpClusters(Int_t iter, TObjArray *trackArray)
2331 // Dump clusters after reco
2332 // signed and unsigned cluster can be visualized
2333 // 1. Unsign all cluster
2334 // 2. Sign all used clusters
2337 Int_t nseed = trackArray->GetEntries();
2338 for (Int_t i=0; i<nseed; i++){
2339 AliTPCseed *pt=(AliTPCseed*)trackArray->UncheckedAt(i);
2343 Bool_t isKink=pt->GetKinkIndex(0)!=0;
2344 for (Int_t j=0; j<160; ++j) {
2345 Int_t index=pt->GetClusterIndex2(j);
2346 if (index<0) continue;
2347 AliTPCclusterMI *c= pt->GetClusterPointer(j);
2349 if (isKink) c->Use(100); // kink
2350 c->Use(10); // by default usage 10
2355 for (Int_t sec=0;sec<fkNIS;sec++){
2356 for (Int_t row=0;row<fInnerSec->GetNRows();row++){
2357 AliTPCclusterMI *cl = fInnerSec[sec][row].GetClusters1();
2358 for (Int_t icl =0;icl< fInnerSec[sec][row].GetN1();icl++){
2359 Float_t gx[3]; cl[icl].GetGlobalXYZ(gx);
2360 (*fDebugStreamer)<<"clDump"<<
2368 cl = fInnerSec[sec][row].GetClusters2();
2369 for (Int_t icl =0;icl< fInnerSec[sec][row].GetN2();icl++){
2370 Float_t gx[3]; cl[icl].GetGlobalXYZ(gx);
2371 (*fDebugStreamer)<<"clDump"<<
2382 for (Int_t sec=0;sec<fkNOS;sec++){
2383 for (Int_t row=0;row<fOuterSec->GetNRows();row++){
2384 AliTPCclusterMI *cl = fOuterSec[sec][row].GetClusters1();
2385 for (Int_t icl =0;icl< fOuterSec[sec][row].GetN1();icl++){
2386 Float_t gx[3]; cl[icl].GetGlobalXYZ(gx);
2387 (*fDebugStreamer)<<"clDump"<<
2395 cl = fOuterSec[sec][row].GetClusters2();
2396 for (Int_t icl =0;icl< fOuterSec[sec][row].GetN2();icl++){
2397 Float_t gx[3]; cl[icl].GetGlobalXYZ(gx);
2398 (*fDebugStreamer)<<"clDump"<<
2410 void AliTPCtrackerMI::UnsignClusters()
2413 // loop over all clusters and unsign them
2416 for (Int_t sec=0;sec<fkNIS;sec++){
2417 for (Int_t row=0;row<fInnerSec->GetNRows();row++){
2418 AliTPCclusterMI *cl = fInnerSec[sec][row].GetClusters1();
2419 for (Int_t icl =0;icl< fInnerSec[sec][row].GetN1();icl++)
2420 // if (cl[icl].IsUsed(10))
2422 cl = fInnerSec[sec][row].GetClusters2();
2423 for (Int_t icl =0;icl< fInnerSec[sec][row].GetN2();icl++)
2424 //if (cl[icl].IsUsed(10))
2429 for (Int_t sec=0;sec<fkNOS;sec++){
2430 for (Int_t row=0;row<fOuterSec->GetNRows();row++){
2431 AliTPCclusterMI *cl = fOuterSec[sec][row].GetClusters1();
2432 for (Int_t icl =0;icl< fOuterSec[sec][row].GetN1();icl++)
2433 //if (cl[icl].IsUsed(10))
2435 cl = fOuterSec[sec][row].GetClusters2();
2436 for (Int_t icl =0;icl< fOuterSec[sec][row].GetN2();icl++)
2437 //if (cl[icl].IsUsed(10))
2446 void AliTPCtrackerMI::SignClusters(const TObjArray * arr, Float_t fnumber, Float_t fdensity)
2449 //sign clusters to be "used"
2451 // snumber and sdensity sign number of sigmas - bellow mean value to be accepted
2452 // loop over "primaries"
2466 Int_t nseed = arr->GetEntriesFast();
2467 for (Int_t i=0; i<nseed; i++) {
2468 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2472 if (!(pt->IsActive())) continue;
2473 Float_t dens = pt->GetNumberOfClusters()/Float_t(pt->GetNFoundable());
2474 if ( (dens>0.7) && (pt->GetNumberOfClusters()>70)){
2476 sumdens2+= dens*dens;
2477 sumn += pt->GetNumberOfClusters();
2478 sumn2 += pt->GetNumberOfClusters()*pt->GetNumberOfClusters();
2479 Float_t chi2 = pt->GetChi2()/pt->GetNumberOfClusters();
2482 sumchi2 +=chi2*chi2;
2487 Float_t mdensity = 0.9;
2488 Float_t meann = 130;
2489 Float_t meanchi = 1;
2490 Float_t sdensity = 0.1;
2491 Float_t smeann = 10;
2492 Float_t smeanchi =0.4;
2496 mdensity = sumdens/sum;
2498 meanchi = sumchi/sum;
2500 sdensity = sumdens2/sum-mdensity*mdensity;
2502 sdensity = TMath::Sqrt(sdensity);
2506 smeann = sumn2/sum-meann*meann;
2508 smeann = TMath::Sqrt(smeann);
2512 smeanchi = sumchi2/sum - meanchi*meanchi;
2514 smeanchi = TMath::Sqrt(smeanchi);
2520 //REMOVE SHORT DELTAS or tracks going out of sensitive volume of TPC
2522 for (Int_t i=0; i<nseed; i++) {
2523 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2527 if (pt->GetBSigned()) continue;
2528 if (pt->GetBConstrain()) continue;
2529 //if (!(pt->IsActive())) continue;
2531 Int_t found,foundable,shared;
2532 pt->GetClusterStatistic(0,160,found, foundable,shared);
2533 if (shared/float(found)>0.3) {
2534 if (shared/float(found)>0.9 ){
2535 //delete arr->RemoveAt(i);
2540 Bool_t isok =kFALSE;
2541 if ( (pt->GetNShared()/pt->GetNumberOfClusters()<0.5) &&pt->GetNumberOfClusters()>60)
2543 if ((TMath::Abs(1/pt->GetC())<100.) && (pt->GetNShared()/pt->GetNumberOfClusters()<0.7))
2545 if (TMath::Abs(pt->GetZ()/pt->GetX())>1.1)
2547 if ( (TMath::Abs(pt->GetSnp()>0.7) && pt->GetD(0,0)>60.))
2551 for (Int_t j=0; j<160; ++j) {
2552 Int_t index=pt->GetClusterIndex2(j);
2553 if (index<0) continue;
2554 AliTPCclusterMI *c= pt->GetClusterPointer(j);
2556 //if (!(c->IsUsed(10))) c->Use();
2563 Double_t maxchi = meanchi+2.*smeanchi;
2565 for (Int_t i=0; i<nseed; i++) {
2566 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2570 //if (!(pt->IsActive())) continue;
2571 if (pt->GetBSigned()) continue;
2572 Double_t chi = pt->GetChi2()/pt->GetNumberOfClusters();
2573 if (chi>maxchi) continue;
2576 Float_t dens = pt->GetNumberOfClusters()/Float_t(pt->GetNFoundable());
2578 //sign only tracks with enoug big density at the beginning
2580 if ((pt->GetDensityFirst(40)<0.75) && pt->GetNumberOfClusters()<meann) continue;
2583 Double_t mindens = TMath::Max(double(mdensity-sdensity*fdensity*bfactor),0.65);
2584 Double_t minn = TMath::Max(Int_t(meann-fnumber*smeann*bfactor),50);
2586 // if (pt->fBConstrain) mindens = TMath::Max(mdensity-sdensity*fdensity*bfactor,0.65);
2587 if ( (pt->GetRemoval()==10) && (pt->GetSnp()>0.8)&&(dens>mindens))
2590 if ((dens>mindens && pt->GetNumberOfClusters()>minn) && chi<maxchi ){
2591 //Int_t noc=pt->GetNumberOfClusters();
2592 pt->SetBSigned(kTRUE);
2593 for (Int_t j=0; j<160; ++j) {
2595 Int_t index=pt->GetClusterIndex2(j);
2596 if (index<0) continue;
2597 AliTPCclusterMI *c= pt->GetClusterPointer(j);
2599 // if (!(c->IsUsed(10))) c->Use();
2604 // gLastCheck = nseed;
2612 void AliTPCtrackerMI::StopNotActive(const TObjArray * arr, Int_t row0, Float_t th0, Float_t th1, Float_t th2) const
2614 // stop not active tracks
2615 // take th1 as threshold for number of founded to number of foundable on last 10 active rows
2616 // take th2 as threshold for number of founded to number of foundable on last 20 active rows
2617 Int_t nseed = arr->GetEntriesFast();
2619 for (Int_t i=0; i<nseed; i++) {
2620 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2624 if (!(pt->IsActive())) continue;
2625 StopNotActive(pt,row0,th0, th1,th2);
2631 void AliTPCtrackerMI::StopNotActive(AliTPCseed * seed, Int_t row0, Float_t th0, Float_t th1,
2634 // stop not active tracks
2635 // take th1 as threshold for number of founded to number of foundable on last 10 active rows
2636 // take th2 as threshold for number of founded to number of foundable on last 20 active rows
2639 Int_t foundable = 0;
2640 Int_t maxindex = seed->GetLastPoint(); //last foundable row
2641 if (seed->GetNFoundable()*th0 > seed->GetNumberOfClusters()) {
2642 seed->Desactivate(10) ;
2646 for (Int_t i=row0; i<maxindex; i++){
2647 Int_t index = seed->GetClusterIndex2(i);
2648 if (index!=-1) foundable++;
2650 if (foundable<=30) sumgood1++;
2651 if (foundable<=50) {
2658 if (foundable>=30.){
2659 if (sumgood1<(th1*30.)) seed->Desactivate(10);
2662 if (sumgood2<(th2*50.)) seed->Desactivate(10);
2666 Int_t AliTPCtrackerMI::RefitInward(AliESDEvent *event)
2669 // back propagation of ESD tracks
2672 if (!event) return 0;
2673 const Int_t kMaxFriendTracks=2000;
2675 // extract correction object for multiplicity dependence of dEdx
2676 TObjArray * gainCalibArray = AliTPCcalibDB::Instance()->GetTimeGainSplinesRun(event->GetRunNumber());
2678 AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ;
2680 AliFatal("Tranformations not in RefitInward");
2683 transform->SetCurrentRecoParam((AliTPCRecoParam*)AliTPCReconstructor::GetRecoParam());
2684 const AliTPCRecoParam * recoParam = AliTPCcalibDB::Instance()->GetTransform()->GetCurrentRecoParam();
2685 Int_t nContribut = event->GetNumberOfTracks();
2686 TGraphErrors * graphMultDependenceDeDx = 0x0;
2687 if (recoParam && recoParam->GetUseMultiplicityCorrectionDedx() && gainCalibArray) {
2688 if (recoParam->GetUseTotCharge()) {
2689 graphMultDependenceDeDx = (TGraphErrors *) gainCalibArray->FindObject("TGRAPHERRORS_MEANQTOT_MULTIPLICITYDEPENDENCE_BEAM_ALL");
2691 graphMultDependenceDeDx = (TGraphErrors *) gainCalibArray->FindObject("TGRAPHERRORS_MEANQMAX_MULTIPLICITYDEPENDENCE_BEAM_ALL");
2697 //PrepareForProlongation(fSeeds,1);
2698 PropagateForward2(fSeeds);
2699 RemoveUsed2(fSeeds,0.4,0.4,20);
2701 TObjArray arraySeed(fSeeds->GetEntries());
2702 for (Int_t i=0;i<fSeeds->GetEntries();i++) {
2703 arraySeed.AddAt(fSeeds->At(i),i);
2705 SignShared(&arraySeed);
2706 // FindCurling(fSeeds, event,2); // find multi found tracks
2707 FindSplitted(fSeeds, event,2); // find multi found tracks
2708 if (AliTPCReconstructor::StreamLevel()>5) FindMultiMC(fSeeds, fEvent,2); // find multi found tracks
2711 Int_t nseed = fSeeds->GetEntriesFast();
2712 for (Int_t i=0;i<nseed;i++){
2713 AliTPCseed * seed = (AliTPCseed*) fSeeds->UncheckedAt(i);
2714 if (!seed) continue;
2715 if (seed->GetKinkIndex(0)>0) UpdateKinkQualityD(seed); // update quality informations for kinks
2716 AliESDtrack *esd=event->GetTrack(i);
2718 if (seed->GetNumberOfClusters()<60 && seed->GetNumberOfClusters()<(esd->GetTPCclusters(0) -5)*0.8){
2719 AliExternalTrackParam paramIn;
2720 AliExternalTrackParam paramOut;
2721 Int_t ncl = seed->RefitTrack(seed,¶mIn,¶mOut);
2722 if (AliTPCReconstructor::StreamLevel()>2) {
2723 (*fDebugStreamer)<<"RecoverIn"<<
2727 "pout.="<<¶mOut<<
2732 seed->Set(paramIn.GetX(),paramIn.GetAlpha(),paramIn.GetParameter(),paramIn.GetCovariance());
2733 seed->SetNumberOfClusters(ncl);
2737 seed->PropagateTo(fkParam->GetInnerRadiusLow());
2738 seed->UpdatePoints();
2739 AddCovariance(seed);
2741 seed->CookdEdx(0.02,0.6);
2742 CookLabel(seed,0.1); //For comparison only
2743 esd->SetTPCClusterMap(seed->GetClusterMap());
2744 esd->SetTPCSharedMap(seed->GetSharedMap());
2746 if (AliTPCReconstructor::StreamLevel()>1 && seed!=0) {
2747 TTreeSRedirector &cstream = *fDebugStreamer;
2754 if (seed->GetNumberOfClusters()>15){
2755 esd->UpdateTrackParams(seed,AliESDtrack::kTPCrefit);
2756 esd->SetTPCPoints(seed->GetPoints());
2757 esd->SetTPCPointsF(seed->GetNFoundable());
2758 Int_t ndedx = seed->GetNCDEDX(0);
2759 Float_t sdedx = seed->GetSDEDX(0);
2760 Float_t dedx = seed->GetdEdx();
2761 // apply mutliplicity dependent dEdx correction if available
2762 if (graphMultDependenceDeDx) {
2763 Double_t corrGain = AliTPCcalibDButil::EvalGraphConst(graphMultDependenceDeDx, nContribut);
2764 dedx += (1 - corrGain)*50.; // MIP is normalized to 50
2766 esd->SetTPCsignal(dedx, sdedx, ndedx);
2768 // add seed to the esd track in Calib level
2770 Bool_t storeFriend = gRandom->Rndm()<(kMaxFriendTracks)/Float_t(nseed);
2771 if (AliTPCReconstructor::StreamLevel()>0 &&storeFriend){
2772 AliTPCseed * seedCopy = new AliTPCseed(*seed, kTRUE);
2773 esd->AddCalibObject(seedCopy);
2778 //printf("problem\n");
2781 //FindKinks(fSeeds,event);
2782 if (AliTPCReconstructor::StreamLevel()>3) DumpClusters(2,fSeeds);
2783 Info("RefitInward","Number of refitted tracks %d",ntracks);
2788 Int_t AliTPCtrackerMI::PropagateBack(AliESDEvent *event)
2791 // back propagation of ESD tracks
2793 if (!event) return 0;
2797 PropagateBack(fSeeds);
2798 RemoveUsed2(fSeeds,0.4,0.4,20);
2799 //FindCurling(fSeeds, fEvent,1);
2800 FindSplitted(fSeeds, event,1); // find multi found tracks
2801 if (AliTPCReconstructor::StreamLevel()>5) FindMultiMC(fSeeds, fEvent,1); // find multi found tracks
2804 Int_t nseed = fSeeds->GetEntriesFast();
2806 for (Int_t i=0;i<nseed;i++){
2807 AliTPCseed * seed = (AliTPCseed*) fSeeds->UncheckedAt(i);
2808 if (!seed) continue;
2809 if (seed->GetKinkIndex(0)<0) UpdateKinkQualityM(seed); // update quality informations for kinks
2810 seed->UpdatePoints();
2811 AddCovariance(seed);
2812 AliESDtrack *esd=event->GetTrack(i);
2813 if (!esd) continue; //never happen
2814 if (seed->GetNumberOfClusters()<60 && seed->GetNumberOfClusters()<(esd->GetTPCclusters(0) -5)*0.8){
2815 AliExternalTrackParam paramIn;
2816 AliExternalTrackParam paramOut;
2817 Int_t ncl = seed->RefitTrack(seed,¶mIn,¶mOut);
2818 if (AliTPCReconstructor::StreamLevel()>2) {
2819 (*fDebugStreamer)<<"RecoverBack"<<
2823 "pout.="<<¶mOut<<
2828 seed->Set(paramOut.GetX(),paramOut.GetAlpha(),paramOut.GetParameter(),paramOut.GetCovariance());
2829 seed->SetNumberOfClusters(ncl);
2832 seed->CookdEdx(0.02,0.6);
2833 CookLabel(seed,0.1); //For comparison only
2834 if (seed->GetNumberOfClusters()>15){
2835 esd->UpdateTrackParams(seed,AliESDtrack::kTPCout);
2836 esd->SetTPCPoints(seed->GetPoints());
2837 esd->SetTPCPointsF(seed->GetNFoundable());
2838 Int_t ndedx = seed->GetNCDEDX(0);
2839 Float_t sdedx = seed->GetSDEDX(0);
2840 Float_t dedx = seed->GetdEdx();
2841 esd->SetTPCsignal(dedx, sdedx, ndedx);
2843 Int_t eventnumber = event->GetEventNumberInFile();// patch 28 fev 06
2844 // This is most likely NOT the event number you'd like to use. It has nothing to do with the 'real' event number
2845 if (AliTPCReconstructor::StreamLevel()>1 && esd) {
2846 (*fDebugStreamer)<<"Cback"<<
2849 "EventNrInFile="<<eventnumber<<
2854 if (AliTPCReconstructor::StreamLevel()>3) DumpClusters(1,fSeeds);
2855 //FindKinks(fSeeds,event);
2856 Info("PropagateBack","Number of back propagated tracks %d",ntracks);
2863 void AliTPCtrackerMI::DeleteSeeds()
2868 Int_t nseed = fSeeds->GetEntriesFast();
2869 for (Int_t i=0;i<nseed;i++){
2870 AliTPCseed * seed = (AliTPCseed*)fSeeds->At(i);
2871 if (seed) delete fSeeds->RemoveAt(i);
2878 void AliTPCtrackerMI::ReadSeeds(AliESDEvent *const event, Int_t direction)
2881 //read seeds from the event
2883 Int_t nentr=event->GetNumberOfTracks();
2885 Info("ReadSeeds", "Number of ESD tracks: %d\n", nentr);
2890 fSeeds = new TObjArray(nentr);
2894 for (Int_t i=0; i<nentr; i++) {
2895 AliESDtrack *esd=event->GetTrack(i);
2896 ULong_t status=esd->GetStatus();
2897 if (!(status&AliESDtrack::kTPCin)) continue;
2898 AliTPCtrack t(*esd);
2899 t.SetNumberOfClusters(0);
2900 // AliTPCseed *seed = new AliTPCseed(t,t.GetAlpha());
2901 AliTPCseed *seed = new AliTPCseed(t/*,t.GetAlpha()*/);
2902 seed->SetUniqueID(esd->GetID());
2903 AddCovariance(seed); //add systematic ucertainty
2904 for (Int_t ikink=0;ikink<3;ikink++) {
2905 Int_t index = esd->GetKinkIndex(ikink);
2906 seed->GetKinkIndexes()[ikink] = index;
2907 if (index==0) continue;
2908 index = TMath::Abs(index);
2909 AliESDkink * kink = fEvent->GetKink(index-1);
2910 if (kink&&esd->GetKinkIndex(ikink)<0){
2911 if ((status & AliESDtrack::kTRDrefit) != 0) kink->SetStatus(1,2);
2912 if ((status & AliESDtrack::kITSout) != 0) kink->SetStatus(1,0);
2914 if (kink&&esd->GetKinkIndex(ikink)>0){
2915 if ((status & AliESDtrack::kTRDrefit) != 0) kink->SetStatus(1,6);
2916 if ((status & AliESDtrack::kITSout) != 0) kink->SetStatus(1,4);
2920 if (((status&AliESDtrack::kITSout)==0)&&(direction==1)) seed->ResetCovariance(10.);
2921 if ( direction ==2 &&(status & AliESDtrack::kTRDrefit) == 0 ) seed->ResetCovariance(10.);
2922 //if ( direction ==2 && ((status & AliESDtrack::kTPCout) == 0) ) {
2923 // fSeeds->AddAt(0,i);
2927 if ( direction ==2 &&(status & AliESDtrack::kTRDrefit) > 0 ) {
2928 Double_t par0[5],par1[5],alpha,x;
2929 esd->GetInnerExternalParameters(alpha,x,par0);
2930 esd->GetExternalParameters(x,par1);
2931 Double_t delta1 = TMath::Abs(par0[4]-par1[4])/(0.000000001+TMath::Abs(par0[4]+par1[4]));
2932 Double_t delta2 = TMath::Abs(par0[3]-par1[3]);
2934 if (esd->GetTRDncls()>0) trdchi2 = esd->GetTRDchi2()/esd->GetTRDncls();
2935 //reset covariance if suspicious
2936 if ( (delta1>0.1) || (delta2>0.006) ||trdchi2>7.)
2937 seed->ResetCovariance(10.);
2942 // rotate to the local coordinate system
2944 fSectors=fInnerSec; fN=fkNIS;
2945 Double_t alpha=seed->GetAlpha() - fSectors->GetAlphaShift();
2946 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
2947 if (alpha < 0. ) alpha += 2.*TMath::Pi();
2948 Int_t ns=Int_t(alpha/fSectors->GetAlpha())%fN;
2949 alpha =ns*fSectors->GetAlpha() + fSectors->GetAlphaShift();
2950 if (alpha<-TMath::Pi()) alpha += 2*TMath::Pi();
2951 if (alpha>=TMath::Pi()) alpha -= 2*TMath::Pi();
2952 alpha-=seed->GetAlpha();
2953 if (!seed->Rotate(alpha)) {
2959 if (esd->GetKinkIndex(0)<=0){
2960 for (Int_t irow=0;irow<160;irow++){
2961 Int_t index = seed->GetClusterIndex2(irow);
2964 AliTPCclusterMI * cl = GetClusterMI(index);
2965 seed->SetClusterPointer(irow,cl);
2967 if ((index & 0x8000)==0){
2968 cl->Use(10); // accepted cluster
2970 cl->Use(6); // close cluster not accepted
2973 Info("ReadSeeds","Not found cluster");
2978 fSeeds->AddAt(seed,i);
2984 //_____________________________________________________________________________
2985 void AliTPCtrackerMI::MakeSeeds3(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, Float_t cuts[4],
2986 Float_t deltay, Int_t ddsec) {
2987 //-----------------------------------------------------------------
2988 // This function creates track seeds.
2989 // SEEDING WITH VERTEX CONSTRAIN
2990 //-----------------------------------------------------------------
2991 // cuts[0] - fP4 cut
2992 // cuts[1] - tan(phi) cut
2993 // cuts[2] - zvertex cut
2994 // cuts[3] - fP3 cut
3002 Double_t x[5], c[15];
3003 // Int_t di = i1-i2;
3005 AliTPCseed * seed = new AliTPCseed();
3006 Double_t alpha=fSectors->GetAlpha(), shift=fSectors->GetAlphaShift();
3007 Double_t cs=cos(alpha), sn=sin(alpha);
3009 // Double_t x1 =fOuterSec->GetX(i1);
3010 //Double_t xx2=fOuterSec->GetX(i2);
3012 Double_t x1 =GetXrow(i1);
3013 Double_t xx2=GetXrow(i2);
3015 Double_t x3=GetX(), y3=GetY(), z3=GetZ();
3017 Int_t imiddle = (i2+i1)/2; //middle pad row index
3018 Double_t xm = GetXrow(imiddle); // radius of middle pad-row
3019 const AliTPCtrackerRow& krm=GetRow(sec,imiddle); //middle pad -row
3023 const AliTPCtrackerRow& kr1=GetRow(ns,i1);
3024 Double_t ymax = GetMaxY(i1)-kr1.GetDeadZone()-1.5;
3025 Double_t ymaxm = GetMaxY(imiddle)-kr1.GetDeadZone()-1.5;
3028 // change cut on curvature if it can't reach this layer
3029 // maximal curvature set to reach it
3030 Double_t dvertexmax = TMath::Sqrt((x1-x3)*(x1-x3)+(ymax+5-y3)*(ymax+5-y3));
3031 if (dvertexmax*0.5*cuts[0]>0.85){
3032 cuts[0] = 0.85/(dvertexmax*0.5+1.);
3034 Double_t r2min = 1/(cuts[0]*cuts[0]); //minimal square of radius given by cut
3037 if (deltay>0) ddsec = 0;
3038 // loop over clusters
3039 for (Int_t is=0; is < kr1; is++) {
3041 if (kr1[is]->IsUsed(10)) continue;
3042 Double_t y1=kr1[is]->GetY(), z1=kr1[is]->GetZ();
3043 //if (TMath::Abs(y1)>ymax) continue;
3045 if (deltay>0 && TMath::Abs(ymax-TMath::Abs(y1))> deltay ) continue; // seed only at the edge
3047 // find possible directions
3048 Float_t anglez = (z1-z3)/(x1-x3);
3049 Float_t extraz = z1 - anglez*(x1-xx2); // extrapolated z
3052 //find rotation angles relative to line given by vertex and point 1
3053 Double_t dvertex2 = (x1-x3)*(x1-x3)+(y1-y3)*(y1-y3);
3054 Double_t dvertex = TMath::Sqrt(dvertex2);
3055 Double_t angle13 = TMath::ATan((y1-y3)/(x1-x3));
3056 Double_t cs13 = cos(-angle13), sn13 = sin(-angle13);
3059 // loop over 2 sectors
3065 Double_t dddz1=0; // direction of delta inclination in z axis
3072 for (Int_t dsec = dsec1; dsec<=dsec2;dsec++){
3073 Int_t sec2 = sec + dsec;
3075 // AliTPCtrackerRow& kr2 = fOuterSec[(sec2+fkNOS)%fkNOS][i2];
3076 //AliTPCtrackerRow& kr2m = fOuterSec[(sec2+fkNOS)%fkNOS][imiddle];
3077 AliTPCtrackerRow& kr2 = GetRow((sec2+fkNOS)%fkNOS,i2);
3078 AliTPCtrackerRow& kr2m = GetRow((sec2+fkNOS)%fkNOS,imiddle);
3079 Int_t index1 = TMath::Max(kr2.Find(extraz-0.6-dddz1*TMath::Abs(z1)*0.05)-1,0);
3080 Int_t index2 = TMath::Min(kr2.Find(extraz+0.6+dddz2*TMath::Abs(z1)*0.05)+1,kr2);
3082 // rotation angles to p1-p3
3083 Double_t cs13r = cos(-angle13+dsec*alpha)/dvertex, sn13r = sin(-angle13+dsec*alpha)/dvertex;
3084 Double_t x2, y2, z2;
3086 // Double_t dymax = maxangle*TMath::Abs(x1-xx2);
3089 Double_t dxx0 = (xx2-x3)*cs13r;
3090 Double_t dyy0 = (xx2-x3)*sn13r;
3091 for (Int_t js=index1; js < index2; js++) {
3092 const AliTPCclusterMI *kcl = kr2[js];
3093 if (kcl->IsUsed(10)) continue;
3095 //calcutate parameters
3097 Double_t yy0 = dyy0 +(kcl->GetY()-y3)*cs13r;
3099 if (TMath::Abs(yy0)<0.000001) continue;
3100 Double_t xx0 = dxx0 -(kcl->GetY()-y3)*sn13r;
3101 Double_t y0 = 0.5*(xx0*xx0+yy0*yy0-xx0)/yy0;
3102 Double_t r02 = (0.25+y0*y0)*dvertex2;
3103 //curvature (radius) cut
3104 if (r02<r2min) continue;
3108 Double_t c0 = 1/TMath::Sqrt(r02);
3112 //Double_t dfi0 = 2.*TMath::ASin(dvertex*c0*0.5);
3113 //Double_t dfi1 = 2.*TMath::ASin(TMath::Sqrt(yy0*yy0+(1-xx0)*(1-xx0))*dvertex*c0*0.5);
3114 Double_t dfi0 = 2.*AliTPCFastMath::FastAsin(dvertex*c0*0.5);
3115 Double_t dfi1 = 2.*AliTPCFastMath::FastAsin(TMath::Sqrt(yy0*yy0+(1-xx0)*(1-xx0))*dvertex*c0*0.5);
3118 Double_t z0 = kcl->GetZ();
3119 Double_t zzzz2 = z1-(z1-z3)*dfi1/dfi0;
3120 if (TMath::Abs(zzzz2-z0)>0.5) continue;
3123 Double_t dip = (z1-z0)*c0/dfi1;
3124 Double_t x0 = (0.5*cs13+y0*sn13)*dvertex*c0;
3135 x2= xx2*cs-y2*sn*dsec;
3136 y2=+xx2*sn*dsec+y2*cs;
3146 // do we have cluster at the middle ?
3148 GetProlongation(x1,xm,x,ym,zm);
3150 AliTPCclusterMI * cm=0;
3151 if (TMath::Abs(ym)-ymaxm<0){
3152 cm = krm.FindNearest2(ym,zm,1.0,0.6,dummy);
3153 if ((!cm) || (cm->IsUsed(10))) {
3158 // rotate y1 to system 0
3159 // get state vector in rotated system
3160 Double_t yr1 = (-0.5*sn13+y0*cs13)*dvertex*c0;
3161 Double_t xr2 = x0*cs+yr1*sn*dsec;
3162 Double_t xr[5]={kcl->GetY(),kcl->GetZ(), xr2, dip, c0};
3164 GetProlongation(xx2,xm,xr,ym,zm);
3165 if (TMath::Abs(ym)-ymaxm<0){
3166 cm = kr2m.FindNearest2(ym,zm,1.0,0.6,dummy);
3167 if ((!cm) || (cm->IsUsed(10))) {
3177 dym = ym - cm->GetY();
3178 dzm = zm - cm->GetZ();
3185 Double_t sy1=kr1[is]->GetSigmaY2()*2., sz1=kr1[is]->GetSigmaZ2()*2.;
3186 Double_t sy2=kcl->GetSigmaY2()*2., sz2=kcl->GetSigmaZ2()*2.;
3187 //Double_t sy3=400*3./12., sy=0.1, sz=0.1;
3188 Double_t sy3=25000*x[4]*x[4]+0.1, sy=0.1, sz=0.1;
3189 //Double_t sy3=25000*x[4]*x[4]*60+0.5, sy=0.1, sz=0.1;
3191 Double_t f40=(F1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
3192 Double_t f42=(F1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
3193 Double_t f43=(F1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
3194 Double_t f20=(F2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
3195 Double_t f22=(F2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
3196 Double_t f23=(F2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
3198 Double_t f30=(F3(x1,y1+sy,x2,y2,z1,z2)-x[3])/sy;
3199 Double_t f31=(F3(x1,y1,x2,y2,z1+sz,z2)-x[3])/sz;
3200 Double_t f32=(F3(x1,y1,x2,y2+sy,z1,z2)-x[3])/sy;
3201 Double_t f34=(F3(x1,y1,x2,y2,z1,z2+sz)-x[3])/sz;
3205 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3206 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3207 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3208 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3209 c[13]=f30*sy1*f40+f32*sy2*f42;
3210 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3212 // if (!BuildSeed(kr1[is],kcl,0,x1,x2,x3,x,c)) continue;
3214 UInt_t index=kr1.GetIndex(is);
3215 seed->~AliTPCseed(); // this does not set the pointer to 0...
3216 AliTPCseed *track=new(seed) AliTPCseed(x1, ns*alpha+shift, x, c, index);
3218 track->SetIsSeeding(kTRUE);
3219 track->SetSeed1(i1);
3220 track->SetSeed2(i2);
3221 track->SetSeedType(3);
3225 FollowProlongation(*track, (i1+i2)/2,1);
3226 Int_t foundable,found,shared;
3227 track->GetClusterStatistic((i1+i2)/2,i1, found, foundable, shared, kTRUE);
3228 if ((found<0.55*foundable) || shared>0.5*found || (track->GetSigmaY2()+track->GetSigmaZ2())>0.5){
3230 seed->~AliTPCseed();
3236 FollowProlongation(*track, i2,1);
3240 track->SetBConstrain(1);
3241 // track->fLastPoint = i1+fInnerSec->GetNRows(); // first cluster in track position
3242 track->SetLastPoint(i1); // first cluster in track position
3243 track->SetFirstPoint(track->GetLastPoint());
3245 if (track->GetNumberOfClusters()<(i1-i2)*0.5 ||
3246 track->GetNumberOfClusters() < track->GetNFoundable()*0.6 ||
3247 track->GetNShared()>0.4*track->GetNumberOfClusters() ) {
3249 seed->~AliTPCseed();
3253 // Z VERTEX CONDITION
3254 Double_t zv, bz=GetBz();
3255 if ( !track->GetZAt(0.,bz,zv) ) continue;
3256 if (TMath::Abs(zv-z3)>cuts[2]) {
3257 FollowProlongation(*track, TMath::Max(i2-20,0));
3258 if ( !track->GetZAt(0.,bz,zv) ) continue;
3259 if (TMath::Abs(zv-z3)>cuts[2]){
3260 FollowProlongation(*track, TMath::Max(i2-40,0));
3261 if ( !track->GetZAt(0.,bz,zv) ) continue;
3262 if (TMath::Abs(zv-z3)>cuts[2] &&(track->GetNumberOfClusters() > track->GetNFoundable()*0.7)){
3263 // make seed without constrain
3264 AliTPCseed * track2 = MakeSeed(track,0.2,0.5,1.);
3265 FollowProlongation(*track2, i2,1);
3266 track2->SetBConstrain(kFALSE);
3267 track2->SetSeedType(1);
3268 arr->AddLast(track2);
3270 seed->~AliTPCseed();
3275 seed->~AliTPCseed();
3282 track->SetSeedType(0);
3283 arr->AddLast(track);
3284 seed = new AliTPCseed;
3286 // don't consider other combinations
3287 if (track->GetNumberOfClusters() > track->GetNFoundable()*0.8)
3293 Info("MakeSeeds3","\nSeeding statistic:\t%d\t%d\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin,nout1,nout2);
3299 void AliTPCtrackerMI::MakeSeeds5(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, Float_t cuts[4],
3304 //-----------------------------------------------------------------
3305 // This function creates track seeds.
3306 //-----------------------------------------------------------------
3307 // cuts[0] - fP4 cut
3308 // cuts[1] - tan(phi) cut
3309 // cuts[2] - zvertex cut
3310 // cuts[3] - fP3 cut
3320 Double_t x[5], c[15];
3322 // make temporary seed
3323 AliTPCseed * seed = new AliTPCseed;
3324 Double_t alpha=fOuterSec->GetAlpha(), shift=fOuterSec->GetAlphaShift();
3325 // Double_t cs=cos(alpha), sn=sin(alpha);
3330 Double_t x1 = GetXrow(i1-1);
3331 const AliTPCtrackerRow& kr1=GetRow(sec,i1-1);
3332 Double_t y1max = GetMaxY(i1-1)-kr1.GetDeadZone()-1.5;
3334 Double_t x1p = GetXrow(i1);
3335 const AliTPCtrackerRow& kr1p=GetRow(sec,i1);
3337 Double_t x1m = GetXrow(i1-2);
3338 const AliTPCtrackerRow& kr1m=GetRow(sec,i1-2);
3341 //last 3 padrow for seeding
3342 AliTPCtrackerRow& kr3 = GetRow((sec+fkNOS)%fkNOS,i1-7);
3343 Double_t x3 = GetXrow(i1-7);
3344 // Double_t y3max= GetMaxY(i1-7)-kr3.fDeadZone-1.5;
3346 AliTPCtrackerRow& kr3p = GetRow((sec+fkNOS)%fkNOS,i1-6);
3347 Double_t x3p = GetXrow(i1-6);
3349 AliTPCtrackerRow& kr3m = GetRow((sec+fkNOS)%fkNOS,i1-8);
3350 Double_t x3m = GetXrow(i1-8);
3355 Int_t im = i1-4; //middle pad row index
3356 Double_t xm = GetXrow(im); // radius of middle pad-row
3357 const AliTPCtrackerRow& krm=GetRow(sec,im); //middle pad -row
3358 // Double_t ymmax = GetMaxY(im)-kr1.fDeadZone-1.5;
3361 Double_t deltax = x1-x3;
3362 Double_t dymax = deltax*cuts[1];
3363 Double_t dzmax = deltax*cuts[3];
3365 // loop over clusters
3366 for (Int_t is=0; is < kr1; is++) {
3368 if (kr1[is]->IsUsed(10)) continue;
3369 Double_t y1=kr1[is]->GetY(), z1=kr1[is]->GetZ();
3371 if (deltay>0 && TMath::Abs(y1max-TMath::Abs(y1))> deltay ) continue; // seed only at the edge
3373 Int_t index1 = TMath::Max(kr3.Find(z1-dzmax)-1,0);
3374 Int_t index2 = TMath::Min(kr3.Find(z1+dzmax)+1,kr3);
3380 for (Int_t js=index1; js < index2; js++) {
3381 const AliTPCclusterMI *kcl = kr3[js];
3382 if (kcl->IsUsed(10)) continue;
3384 // apply angular cuts
3385 if (TMath::Abs(y1-y3)>dymax) continue;
3388 if (TMath::Abs(z1-z3)>dzmax) continue;
3390 Double_t angley = (y1-y3)/(x1-x3);
3391 Double_t anglez = (z1-z3)/(x1-x3);
3393 Double_t erry = TMath::Abs(angley)*(x1-x1m)*0.5+0.5;
3394 Double_t errz = TMath::Abs(anglez)*(x1-x1m)*0.5+0.5;
3396 Double_t yyym = angley*(xm-x1)+y1;
3397 Double_t zzzm = anglez*(xm-x1)+z1;
3399 const AliTPCclusterMI *kcm = krm.FindNearest2(yyym,zzzm,erry,errz,index);
3401 if (kcm->IsUsed(10)) continue;
3403 erry = TMath::Abs(angley)*(x1-x1m)*0.4+0.5;
3404 errz = TMath::Abs(anglez)*(x1-x1m)*0.4+0.5;
3411 // look around first
3412 const AliTPCclusterMI *kc1m = kr1m.FindNearest2(angley*(x1m-x1)+y1,
3418 if (kc1m->IsUsed(10)) used++;
3420 const AliTPCclusterMI *kc1p = kr1p.FindNearest2(angley*(x1p-x1)+y1,
3426 if (kc1p->IsUsed(10)) used++;
3428 if (used>1) continue;
3429 if (found<1) continue;
3433 const AliTPCclusterMI *kc3m = kr3m.FindNearest2(angley*(x3m-x3)+y3,
3439 if (kc3m->IsUsed(10)) used++;
3443 const AliTPCclusterMI *kc3p = kr3p.FindNearest2(angley*(x3p-x3)+y3,
3449 if (kc3p->IsUsed(10)) used++;
3453 if (used>1) continue;
3454 if (found<3) continue;
3464 x[4]=F1(x1,y1,x2,y2,x3,y3);
3465 //if (TMath::Abs(x[4]) >= cuts[0]) continue;
3468 x[2]=F2(x1,y1,x2,y2,x3,y3);
3471 x[3]=F3n(x1,y1,x2,y2,z1,z2,x[4]);
3472 //if (TMath::Abs(x[3]) > cuts[3]) continue;
3476 Double_t sy1=0.1, sz1=0.1;
3477 Double_t sy2=0.1, sz2=0.1;
3478 Double_t sy3=0.1, sy=0.1, sz=0.1;
3480 Double_t f40=(F1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
3481 Double_t f42=(F1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
3482 Double_t f43=(F1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
3483 Double_t f20=(F2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
3484 Double_t f22=(F2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
3485 Double_t f23=(F2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
3487 Double_t f30=(F3(x1,y1+sy,x2,y2,z1,z2)-x[3])/sy;
3488 Double_t f31=(F3(x1,y1,x2,y2,z1+sz,z2)-x[3])/sz;
3489 Double_t f32=(F3(x1,y1,x2,y2+sy,z1,z2)-x[3])/sy;
3490 Double_t f34=(F3(x1,y1,x2,y2,z1,z2+sz)-x[3])/sz;
3494 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3495 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3496 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3497 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3498 c[13]=f30*sy1*f40+f32*sy2*f42;
3499 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3501 // if (!BuildSeed(kr1[is],kcl,0,x1,x2,x3,x,c)) continue;
3503 index=kr1.GetIndex(is);
3504 seed->~AliTPCseed();
3505 AliTPCseed *track=new(seed) AliTPCseed(x1, sec*alpha+shift, x, c, index);
3507 track->SetIsSeeding(kTRUE);
3510 FollowProlongation(*track, i1-7,1);
3511 if (track->GetNumberOfClusters() < track->GetNFoundable()*0.75 ||
3512 track->GetNShared()>0.6*track->GetNumberOfClusters() || ( track->GetSigmaY2()+ track->GetSigmaZ2())>0.6){
3514 seed->~AliTPCseed();
3520 FollowProlongation(*track, i2,1);
3521 track->SetBConstrain(0);
3522 track->SetLastPoint(i1+fInnerSec->GetNRows()); // first cluster in track position
3523 track->SetFirstPoint(track->GetLastPoint());
3525 if (track->GetNumberOfClusters()<(i1-i2)*0.5 ||
3526 track->GetNumberOfClusters()<track->GetNFoundable()*0.7 ||
3527 track->GetNShared()>2. || track->GetChi2()/track->GetNumberOfClusters()>6 || ( track->GetSigmaY2()+ track->GetSigmaZ2())>0.5 ) {
3529 seed->~AliTPCseed();
3534 FollowProlongation(*track, TMath::Max(i2-10,0),1);
3535 AliTPCseed * track2 = MakeSeed(track,0.2,0.5,0.9);
3536 FollowProlongation(*track2, i2,1);
3537 track2->SetBConstrain(kFALSE);
3538 track2->SetSeedType(4);
3539 arr->AddLast(track2);
3541 seed->~AliTPCseed();
3545 //arr->AddLast(track);
3546 //seed = new AliTPCseed;
3552 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);
3558 //_____________________________________________________________________________
3559 void AliTPCtrackerMI::MakeSeeds2(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, Float_t */*cuts[4]*/,
3560 Float_t deltay, Bool_t /*bconstrain*/) {
3561 //-----------------------------------------------------------------
3562 // This function creates track seeds - without vertex constraint
3563 //-----------------------------------------------------------------
3564 // cuts[0] - fP4 cut - not applied
3565 // cuts[1] - tan(phi) cut
3566 // cuts[2] - zvertex cut - not applied
3567 // cuts[3] - fP3 cut
3577 Double_t alpha=fOuterSec->GetAlpha(), shift=fOuterSec->GetAlphaShift();
3578 // Double_t cs=cos(alpha), sn=sin(alpha);
3579 Int_t row0 = (i1+i2)/2;
3580 Int_t drow = (i1-i2)/2;
3581 const AliTPCtrackerRow& kr0=fSectors[sec][row0];
3582 AliTPCtrackerRow * kr=0;
3584 AliTPCpolyTrack polytrack;
3585 Int_t nclusters=fSectors[sec][row0];
3586 AliTPCseed * seed = new AliTPCseed;
3591 for (Int_t is=0; is < nclusters; is++) { //LOOP over clusters
3593 Int_t nfoundable =0;
3594 for (Int_t iter =1; iter<2; iter++){ //iterations
3595 const AliTPCtrackerRow& krm=fSectors[sec][row0-iter];
3596 const AliTPCtrackerRow& krp=fSectors[sec][row0+iter];
3597 const AliTPCclusterMI * cl= kr0[is];
3599 if (cl->IsUsed(10)) {
3605 Double_t x = kr0.GetX();
3606 // Initialization of the polytrack
3611 Double_t y0= cl->GetY();
3612 Double_t z0= cl->GetZ();
3616 Double_t ymax = fSectors->GetMaxY(row0)-kr0.GetDeadZone()-1.5;
3617 if (deltay>0 && TMath::Abs(ymax-TMath::Abs(y0))> deltay ) continue; // seed only at the edge
3619 erry = (0.5)*cl->GetSigmaY2()/TMath::Sqrt(cl->GetQ())*6;
3620 errz = (0.5)*cl->GetSigmaZ2()/TMath::Sqrt(cl->GetQ())*6;
3621 polytrack.AddPoint(x,y0,z0,erry, errz);
3624 if (cl->IsUsed(10)) sumused++;
3627 Float_t roady = (5*TMath::Sqrt(cl->GetSigmaY2()+0.2)+1.)*iter;
3628 Float_t roadz = (5*TMath::Sqrt(cl->GetSigmaZ2()+0.2)+1.)*iter;
3631 AliTPCclusterMI * cl1 = krm.FindNearest(y0,z0,roady,roadz);
3632 if (cl1 && TMath::Abs(ymax-TMath::Abs(y0))) {
3633 erry = (0.5)*cl1->GetSigmaY2()/TMath::Sqrt(cl1->GetQ())*3;
3634 errz = (0.5)*cl1->GetSigmaZ2()/TMath::Sqrt(cl1->GetQ())*3;
3635 if (cl1->IsUsed(10)) sumused++;
3636 polytrack.AddPoint(x,cl1->GetY(),cl1->GetZ(),erry,errz);
3640 AliTPCclusterMI * cl2 = krp.FindNearest(y0,z0,roady,roadz);
3642 erry = (0.5)*cl2->GetSigmaY2()/TMath::Sqrt(cl2->GetQ())*3;
3643 errz = (0.5)*cl2->GetSigmaZ2()/TMath::Sqrt(cl2->GetQ())*3;
3644 if (cl2->IsUsed(10)) sumused++;
3645 polytrack.AddPoint(x,cl2->GetY(),cl2->GetZ(),erry,errz);
3648 if (sumused>0) continue;
3650 polytrack.UpdateParameters();
3656 nfoundable = polytrack.GetN();
3657 nfound = nfoundable;
3659 for (Int_t ddrow = iter+1; ddrow<drow;ddrow++){
3660 Float_t maxdist = 0.8*(1.+3./(ddrow));
3661 for (Int_t delta = -1;delta<=1;delta+=2){
3662 Int_t row = row0+ddrow*delta;
3663 kr = &(fSectors[sec][row]);
3664 Double_t xn = kr->GetX();
3665 Double_t ymax1 = fSectors->GetMaxY(row)-kr->GetDeadZone()-1.5;
3666 polytrack.GetFitPoint(xn,yn,zn);
3667 if (TMath::Abs(yn)>ymax1) continue;
3669 AliTPCclusterMI * cln = kr->FindNearest(yn,zn,roady,roadz);
3671 Float_t dist = TMath::Sqrt( (yn-cln->GetY())*(yn-cln->GetY())+(zn-cln->GetZ())*(zn-cln->GetZ()));
3674 erry = (dist+0.3)*cln->GetSigmaY2()/TMath::Sqrt(cln->GetQ())*(1.+1./(ddrow));
3675 errz = (dist+0.3)*cln->GetSigmaZ2()/TMath::Sqrt(cln->GetQ())*(1.+1./(ddrow));
3676 if (cln->IsUsed(10)) {
3677 // printf("used\n");
3685 polytrack.AddPoint(xn,cln->GetY(),cln->GetZ(),erry, errz);
3690 if ( (sumused>3) || (sumused>0.5*nfound) || (nfound<0.6*nfoundable)) break;
3691 polytrack.UpdateParameters();
3694 if ( (sumused>3) || (sumused>0.5*nfound)) {
3695 //printf("sumused %d\n",sumused);
3700 polytrack.GetFitDerivation(kr0.GetX(),dy,dz);
3701 AliTPCpolyTrack track2;
3703 polytrack.Refit(track2,0.5+TMath::Abs(dy)*0.3,0.4+TMath::Abs(dz)*0.3);
3704 if (track2.GetN()<0.5*nfoundable) continue;
3707 if ((nfound>0.6*nfoundable) &&( nfoundable>0.4*(i1-i2))) {
3709 // test seed with and without constrain
3710 for (Int_t constrain=0; constrain<=0;constrain++){
3711 // add polytrack candidate
3713 Double_t x[5], c[15];
3714 Double_t x1,x2,x3,y1,y2,y3,z1,z2,z3;
3715 track2.GetBoundaries(x3,x1);
3717 track2.GetFitPoint(x1,y1,z1);
3718 track2.GetFitPoint(x2,y2,z2);
3719 track2.GetFitPoint(x3,y3,z3);
3721 //is track pointing to the vertex ?
3724 polytrack.GetFitPoint(x0,y0,z0);
3737 x[4]=F1(x1,y1,x2,y2,x3,y3);
3739 // if (TMath::Abs(x[4]) >= cuts[0]) continue; //
3740 x[2]=F2(x1,y1,x2,y2,x3,y3);
3742 //if (TMath::Abs(x[4]*x1-x[2]) >= cuts[1]) continue;
3743 //x[3]=F3(x1,y1,x2,y2,z1,z2);
3744 x[3]=F3n(x1,y1,x3,y3,z1,z3,x[4]);
3745 //if (TMath::Abs(x[3]) > cuts[3]) continue;
3748 Double_t sy =0.1, sz =0.1;
3749 Double_t sy1=0.02, sz1=0.02;
3750 Double_t sy2=0.02, sz2=0.02;
3754 sy3=25000*x[4]*x[4]+0.1, sy=0.1, sz=0.1;
3757 Double_t f40=(F1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
3758 Double_t f42=(F1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
3759 Double_t f43=(F1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
3760 Double_t f20=(F2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
3761 Double_t f22=(F2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
3762 Double_t f23=(F2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
3764 Double_t f30=(F3(x1,y1+sy,x3,y3,z1,z3)-x[3])/sy;
3765 Double_t f31=(F3(x1,y1,x3,y3,z1+sz,z3)-x[3])/sz;
3766 Double_t f32=(F3(x1,y1,x3,y3+sy,z1,z3)-x[3])/sy;
3767 Double_t f34=(F3(x1,y1,x3,y3,z1,z3+sz)-x[3])/sz;
3772 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3773 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3774 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3775 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3776 c[13]=f30*sy1*f40+f32*sy2*f42;
3777 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3779 //Int_t row1 = fSectors->GetRowNumber(x1);
3780 Int_t row1 = GetRowNumber(x1);
3784 seed->~AliTPCseed();
3785 AliTPCseed *track=new(seed) AliTPCseed(x1,sec*alpha+shift,x,c,index);
3786 track->SetIsSeeding(kTRUE);
3787 Int_t rc=FollowProlongation(*track, i2);
3788 if (constrain) track->SetBConstrain(1);
3790 track->SetBConstrain(0);
3791 track->SetLastPoint(row1+fInnerSec->GetNRows()); // first cluster in track position
3792 track->SetFirstPoint(track->GetLastPoint());
3794 if (rc==0 || track->GetNumberOfClusters()<(i1-i2)*0.5 ||
3795 track->GetNumberOfClusters() < track->GetNFoundable()*0.6 ||
3796 track->GetNShared()>0.4*track->GetNumberOfClusters()) {
3799 seed->~AliTPCseed();
3802 arr->AddLast(track);
3803 seed = new AliTPCseed;
3807 } // if accepted seed
3810 Info("MakeSeeds2","\nSeeding statiistic:\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin3);
3816 AliTPCseed *AliTPCtrackerMI::MakeSeed(AliTPCseed *const track, Float_t r0, Float_t r1, Float_t r2)
3820 //reseed using track points
3821 Int_t p0 = int(r0*track->GetNumberOfClusters()); // point 0
3822 Int_t p1 = int(r1*track->GetNumberOfClusters());
3823 Int_t p2 = int(r2*track->GetNumberOfClusters()); // last point
3825 Double_t x0[3],x1[3],x2[3];
3826 for (Int_t i=0;i<3;i++){
3832 // find track position at given ratio of the length
3833 Int_t sec0=0, sec1=0, sec2=0;
3836 for (Int_t i=0;i<160;i++){
3837 if (track->GetClusterPointer(i)){
3839 AliTPCTrackerPoint *trpoint =track->GetTrackPoint(i);
3840 if ( (index<p0) || x0[0]<0 ){
3841 if (trpoint->GetX()>1){
3842 clindex = track->GetClusterIndex2(i);
3844 x0[0] = trpoint->GetX();
3845 x0[1] = trpoint->GetY();
3846 x0[2] = trpoint->GetZ();
3847 sec0 = ((clindex&0xff000000)>>24)%18;
3852 if ( (index<p1) &&(trpoint->GetX()>1)){
3853 clindex = track->GetClusterIndex2(i);
3855 x1[0] = trpoint->GetX();
3856 x1[1] = trpoint->GetY();
3857 x1[2] = trpoint->GetZ();
3858 sec1 = ((clindex&0xff000000)>>24)%18;
3861 if ( (index<p2) &&(trpoint->GetX()>1)){
3862 clindex = track->GetClusterIndex2(i);
3864 x2[0] = trpoint->GetX();
3865 x2[1] = trpoint->GetY();
3866 x2[2] = trpoint->GetZ();
3867 sec2 = ((clindex&0xff000000)>>24)%18;
3874 Double_t alpha, cs,sn, xx2,yy2;
3876 alpha = (sec1-sec2)*fSectors->GetAlpha();
3877 cs = TMath::Cos(alpha);
3878 sn = TMath::Sin(alpha);
3879 xx2= x1[0]*cs-x1[1]*sn;
3880 yy2= x1[0]*sn+x1[1]*cs;
3884 alpha = (sec0-sec2)*fSectors->GetAlpha();
3885 cs = TMath::Cos(alpha);
3886 sn = TMath::Sin(alpha);
3887 xx2= x0[0]*cs-x0[1]*sn;
3888 yy2= x0[0]*sn+x0[1]*cs;
3894 Double_t x[5],c[15];
3898 x[4]=F1(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]);
3899 // if (x[4]>1) return 0;
3900 x[2]=F2(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]);
3901 x[3]=F3n(x2[0],x2[1],x0[0],x0[1],x2[2],x0[2],x[4]);
3902 //if (TMath::Abs(x[3]) > 2.2) return 0;
3903 //if (TMath::Abs(x[2]) > 1.99) return 0;
3905 Double_t sy =0.1, sz =0.1;
3907 Double_t sy1=0.02+track->GetSigmaY2(), sz1=0.02+track->GetSigmaZ2();
3908 Double_t sy2=0.01+track->GetSigmaY2(), sz2=0.01+track->GetSigmaZ2();
3909 Double_t sy3=0.01+track->GetSigmaY2();
3911 Double_t f40=(F1(x2[0],x2[1]+sy,x1[0],x1[1],x0[0],x0[1])-x[4])/sy;
3912 Double_t f42=(F1(x2[0],x2[1],x1[0],x1[1]+sy,x0[0],x0[1])-x[4])/sy;
3913 Double_t f43=(F1(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]+sy)-x[4])/sy;
3914 Double_t f20=(F2(x2[0],x2[1]+sy,x1[0],x1[1],x0[0],x0[1])-x[2])/sy;
3915 Double_t f22=(F2(x2[0],x2[1],x1[0],x1[1]+sy,x0[0],x0[1])-x[2])/sy;
3916 Double_t f23=(F2(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]+sy)-x[2])/sy;
3918 Double_t f30=(F3(x2[0],x2[1]+sy,x0[0],x0[1],x2[2],x0[2])-x[3])/sy;
3919 Double_t f31=(F3(x2[0],x2[1],x0[0],x0[1],x2[2]+sz,x0[2])-x[3])/sz;
3920 Double_t f32=(F3(x2[0],x2[1],x0[0],x0[1]+sy,x2[2],x0[2])-x[3])/sy;
3921 Double_t f34=(F3(x2[0],x2[1],x0[0],x0[1],x2[2],x0[2]+sz)-x[3])/sz;
3926 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3927 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3928 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3929 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3930 c[13]=f30*sy1*f40+f32*sy2*f42;
3931 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3933 // Int_t row1 = fSectors->GetRowNumber(x2[0]);
3934 AliTPCseed *seed=new AliTPCseed(x2[0], sec2*fSectors->GetAlpha()+fSectors->GetAlphaShift(), x, c, 0);
3935 // Double_t y0,z0,y1,z1, y2,z2;
3936 //seed->GetProlongation(x0[0],y0,z0);
3937 // seed->GetProlongation(x1[0],y1,z1);
3938 //seed->GetProlongation(x2[0],y2,z2);
3940 seed->SetLastPoint(pp2);
3941 seed->SetFirstPoint(pp2);
3948 AliTPCseed *AliTPCtrackerMI::ReSeed(const AliTPCseed *track, Float_t r0, Float_t r1, Float_t r2)
3952 //reseed using founded clusters
3954 // Find the number of clusters
3955 Int_t nclusters = 0;
3956 for (Int_t irow=0;irow<160;irow++){
3957 if (track->GetClusterIndex(irow)>0) nclusters++;
3961 ipos[0] = TMath::Max(int(r0*nclusters),0); // point 0 cluster
3962 ipos[1] = TMath::Min(int(r1*nclusters),nclusters-1); //
3963 ipos[2] = TMath::Min(int(r2*nclusters),nclusters-1); // last point
3966 Double_t xyz[3][3]={{0}};
3967 Int_t row[3]={0},sec[3]={0,0,0};
3969 // find track row position at given ratio of the length
3971 for (Int_t irow=0;irow<160;irow++){
3972 if (track->GetClusterIndex2(irow)<0) continue;
3974 for (Int_t ipoint=0;ipoint<3;ipoint++){
3975 if (index<=ipos[ipoint]) row[ipoint] = irow;
3979 //Get cluster and sector position
3980 for (Int_t ipoint=0;ipoint<3;ipoint++){
3981 Int_t clindex = track->GetClusterIndex2(row[ipoint]);
3982 AliTPCclusterMI * cl = GetClusterMI(clindex);
3985 // AliTPCclusterMI * cl = GetClusterMI(clindex);
3988 sec[ipoint] = ((clindex&0xff000000)>>24)%18;
3989 xyz[ipoint][0] = GetXrow(row[ipoint]);
3990 xyz[ipoint][1] = cl->GetY();
3991 xyz[ipoint][2] = cl->GetZ();
3995 // Calculate seed state vector and covariance matrix
3997 Double_t alpha, cs,sn, xx2,yy2;
3999 alpha = (sec[1]-sec[2])*fSectors->GetAlpha();
4000 cs = TMath::Cos(alpha);
4001 sn = TMath::Sin(alpha);
4002 xx2= xyz[1][0]*cs-xyz[1][1]*sn;
4003 yy2= xyz[1][0]*sn+xyz[1][1]*cs;
4007 alpha = (sec[0]-sec[2])*fSectors->GetAlpha();
4008 cs = TMath::Cos(alpha);
4009 sn = TMath::Sin(alpha);
4010 xx2= xyz[0][0]*cs-xyz[0][1]*sn;
4011 yy2= xyz[0][0]*sn+xyz[0][1]*cs;
4017 Double_t x[5],c[15];
4021 x[4]=F1(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
4022 x[2]=F2(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
4023 x[3]=F3n(xyz[2][0],xyz[2][1],xyz[0][0],xyz[0][1],xyz[2][2],xyz[0][2],x[4]);
4025 Double_t sy =0.1, sz =0.1;
4027 Double_t sy1=0.2, sz1=0.2;
4028 Double_t sy2=0.2, sz2=0.2;
4031 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;
4032 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;
4033 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;
4034 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;
4035 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;
4036 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;
4038 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;
4039 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;
4040 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;
4041 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;
4046 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
4047 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
4048 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
4049 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
4050 c[13]=f30*sy1*f40+f32*sy2*f42;
4051 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
4053 // Int_t row1 = fSectors->GetRowNumber(xyz[2][0]);
4054 AliTPCseed *seed=new AliTPCseed(xyz[2][0], sec[2]*fSectors->GetAlpha()+fSectors->GetAlphaShift(), x, c, 0);
4055 seed->SetLastPoint(row[2]);
4056 seed->SetFirstPoint(row[2]);
4061 AliTPCseed *AliTPCtrackerMI::ReSeed(AliTPCseed *track,Int_t r0, Bool_t forward)
4065 //reseed using founded clusters
4068 Int_t row[3]={0,0,0};
4069 Int_t sec[3]={0,0,0};
4071 // forward direction
4073 for (Int_t irow=r0;irow<160;irow++){
4074 if (track->GetClusterIndex(irow)>0){
4079 for (Int_t irow=160;irow>r0;irow--){
4080 if (track->GetClusterIndex(irow)>0){
4085 for (Int_t irow=row[2]-15;irow>row[0];irow--){
4086 if (track->GetClusterIndex(irow)>0){
4094 for (Int_t irow=0;irow<r0;irow++){
4095 if (track->GetClusterIndex(irow)>0){
4100 for (Int_t irow=r0;irow>0;irow--){
4101 if (track->GetClusterIndex(irow)>0){
4106 for (Int_t irow=row[2]-15;irow>row[0];irow--){
4107 if (track->GetClusterIndex(irow)>0){
4114 if ((row[2]-row[0])<20) return 0;
4115 if (row[1]==0) return 0;
4118 //Get cluster and sector position
4119 for (Int_t ipoint=0;ipoint<3;ipoint++){
4120 Int_t clindex = track->GetClusterIndex2(row[ipoint]);
4121 AliTPCclusterMI * cl = GetClusterMI(clindex);
4124 // AliTPCclusterMI * cl = GetClusterMI(clindex);
4127 sec[ipoint] = ((clindex&0xff000000)>>24)%18;
4128 xyz[ipoint][0] = GetXrow(row[ipoint]);
4129 AliTPCTrackerPoint * point = track->GetTrackPoint(row[ipoint]);
4130 if (point&&ipoint<2){
4132 xyz[ipoint][1] = point->GetY();
4133 xyz[ipoint][2] = point->GetZ();
4136 xyz[ipoint][1] = cl->GetY();
4137 xyz[ipoint][2] = cl->GetZ();
4144 // Calculate seed state vector and covariance matrix
4146 Double_t alpha, cs,sn, xx2,yy2;
4148 alpha = (sec[1]-sec[2])*fSectors->GetAlpha();
4149 cs = TMath::Cos(alpha);
4150 sn = TMath::Sin(alpha);
4151 xx2= xyz[1][0]*cs-xyz[1][1]*sn;
4152 yy2= xyz[1][0]*sn+xyz[1][1]*cs;
4156 alpha = (sec[0]-sec[2])*fSectors->GetAlpha();
4157 cs = TMath::Cos(alpha);
4158 sn = TMath::Sin(alpha);
4159 xx2= xyz[0][0]*cs-xyz[0][1]*sn;
4160 yy2= xyz[0][0]*sn+xyz[0][1]*cs;
4166 Double_t x[5],c[15];
4170 x[4]=F1(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
4171 x[2]=F2(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
4172 x[3]=F3n(xyz[2][0],xyz[2][1],xyz[0][0],xyz[0][1],xyz[2][2],xyz[0][2],x[4]);
4174 Double_t sy =0.1, sz =0.1;
4176 Double_t sy1=0.2, sz1=0.2;
4177 Double_t sy2=0.2, sz2=0.2;
4180 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;
4181 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;
4182 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;
4183 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;
4184 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;
4185 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;
4187 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;
4188 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;
4189 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;
4190 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;
4195 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
4196 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
4197 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
4198 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
4199 c[13]=f30*sy1*f40+f32*sy2*f42;
4200 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
4202 // Int_t row1 = fSectors->GetRowNumber(xyz[2][0]);
4203 AliTPCseed *seed=new AliTPCseed(xyz[2][0], sec[2]*fSectors->GetAlpha()+fSectors->GetAlphaShift(), x, c, 0);
4204 seed->SetLastPoint(row[2]);
4205 seed->SetFirstPoint(row[2]);
4206 for (Int_t i=row[0];i<row[2];i++){
4207 seed->SetClusterIndex(i, track->GetClusterIndex(i));
4215 void AliTPCtrackerMI::FindMultiMC(const TObjArray * array, AliESDEvent */*esd*/, Int_t iter)
4218 // find multi tracks - THIS FUNCTION IS ONLY FOR DEBUG PURPOSES
4220 // Use AliTPCReconstructor::StreamLevel()>2 if you want to tune parameters - cuts
4222 // Two reasons to have multiple find tracks
4223 // 1. Curling tracks can be find more than once
4224 // 2. Splitted tracks
4225 // a.) Multiple seeding to increase tracking efficiency - (~ 100% reached)
4226 // b.) Edge effect on the sector boundaries
4229 // Algorithm done in 2 phases - because of CPU consumption
4230 // it is n^2 algorithm - for lead-lead 20000x20000 combination are investigated
4232 // Algorihm for curling tracks sign:
4233 // 1 phase -makes a very rough fast cuts to minimize combinatorics
4234 // a.) opposite sign
4235 // b.) one of the tracks - not pointing to the primary vertex -
4236 // c.) delta tan(theta)
4238 // 2 phase - calculates DCA between tracks - time consument
4243 // General cuts - for splitted tracks and for curling tracks
4245 const Float_t kMaxdPhi = 0.2; // maximal distance in phi
4247 // Curling tracks cuts
4252 Int_t nentries = array->GetEntriesFast();
4253 AliHelix *helixes = new AliHelix[nentries];
4254 Float_t *xm = new Float_t[nentries];
4255 Float_t *dz0 = new Float_t[nentries];
4256 Float_t *dz1 = new Float_t[nentries];
4262 // Find track COG in x direction - point with best defined parameters
4264 for (Int_t i=0;i<nentries;i++){
4265 AliTPCseed* track = (AliTPCseed*)array->At(i);
4266 if (!track) continue;
4267 track->SetCircular(0);
4268 new (&helixes[i]) AliHelix(*track);
4272 track->GetDZ(GetX(),GetY(),GetZ(),GetBz(),dz);
4275 for (Int_t icl=0; icl<160; icl++){
4276 AliTPCclusterMI * cl = track->GetClusterPointer(icl);
4282 if (ncl>0) xm[i]/=Float_t(ncl);
4285 for (Int_t i0=0;i0<nentries;i0++){
4286 AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
4287 if (!track0) continue;
4288 Float_t xc0 = helixes[i0].GetHelix(6);
4289 Float_t yc0 = helixes[i0].GetHelix(7);
4290 Float_t r0 = helixes[i0].GetHelix(8);
4291 Float_t rc0 = TMath::Sqrt(xc0*xc0+yc0*yc0);
4292 Float_t fi0 = TMath::ATan2(yc0,xc0);
4294 for (Int_t i1=i0+1;i1<nentries;i1++){
4295 AliTPCseed * track1 = (AliTPCseed*)array->At(i1);
4296 if (!track1) continue;
4297 Int_t lab0=track0->GetLabel();
4298 Int_t lab1=track1->GetLabel();
4299 if (TMath::Abs(lab0)!=TMath::Abs(lab1)) continue;
4301 Float_t xc1 = helixes[i1].GetHelix(6);
4302 Float_t yc1 = helixes[i1].GetHelix(7);
4303 Float_t r1 = helixes[i1].GetHelix(8);
4304 Float_t rc1 = TMath::Sqrt(xc1*xc1+yc1*yc1);
4305 Float_t fi1 = TMath::ATan2(yc1,xc1);
4307 Float_t dfi = fi0-fi1;
4310 if (dfi>1.5*TMath::Pi()) dfi-=TMath::Pi(); // take care about edge effect
4311 if (dfi<-1.5*TMath::Pi()) dfi+=TMath::Pi(); //
4312 if (TMath::Abs(dfi)>kMaxdPhi&&helixes[i0].GetHelix(4)*helixes[i1].GetHelix(4)<0){
4314 // if short tracks with undefined sign
4315 fi1 = -TMath::ATan2(yc1,-xc1);
4318 Float_t dtheta = TMath::Abs(track0->GetTgl()-track1->GetTgl())<TMath::Abs(track0->GetTgl()+track1->GetTgl())? track0->GetTgl()-track1->GetTgl():track0->GetTgl()+track1->GetTgl();
4321 // debug stream to tune "fast cuts"
4323 Double_t dist[3]; // distance at X
4324 Double_t mdist[3]={0,0,0}; // mean distance X+-40cm
4325 track0->GetDistance(track1,0.5*(xm[i0]+xm[i1])-40.,dist,AliTracker::GetBz());
4326 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[i]);
4327 track0->GetDistance(track1,0.5*(xm[i0]+xm[i1])+40.,dist,AliTracker::GetBz());
4328 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[i]);
4329 track0->GetDistance(track1,0.5*(xm[i0]+xm[i1]),dist,AliTracker::GetBz());
4330 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[i]);
4331 for (Int_t i=0;i<3;i++) mdist[i]*=0.33333;
4335 for (Int_t icl=0; icl<160; icl++){
4336 AliTPCclusterMI * cl0 = track0->GetClusterPointer(icl);
4337 AliTPCclusterMI * cl1 = track1->GetClusterPointer(icl);
4340 if (cl0==cl1) sums++;
4344 if (AliTPCReconstructor::StreamLevel()>5) {
4345 TTreeSRedirector &cstream = *fDebugStreamer;
4350 "Tr0.="<<track0<< // seed0
4351 "Tr1.="<<track1<< // seed1
4352 "h0.="<<&helixes[i0]<<
4353 "h1.="<<&helixes[i1]<<
4355 "sum="<<sum<< //the sum of rows with cl in both
4356 "sums="<<sums<< //the sum of shared clusters
4357 "xm0="<<xm[i0]<< // the center of track
4358 "xm1="<<xm[i1]<< // the x center of track
4359 // General cut variables
4360 "dfi="<<dfi<< // distance in fi angle
4361 "dtheta="<<dtheta<< // distance int theta angle
4367 "dist0="<<dist[0]<< //distance x
4368 "dist1="<<dist[1]<< //distance y
4369 "dist2="<<dist[2]<< //distance z
4370 "mdist0="<<mdist[0]<< //distance x
4371 "mdist1="<<mdist[1]<< //distance y
4372 "mdist2="<<mdist[2]<< //distance z
4388 if (AliTPCReconstructor::StreamLevel()>1) {
4389 AliInfo("Time for curling tracks removal DEBUGGING MC");
4396 void AliTPCtrackerMI::FindSplitted(TObjArray * array, AliESDEvent */*esd*/, Int_t /*iter*/){
4398 // Find Splitted tracks and remove the one with worst quality
4399 // Corresponding debug streamer to tune selections - "Splitted2"
4401 // 0. Sort tracks according quility
4402 // 1. Propagate the tracks to the reference radius
4403 // 2. Double_t loop to select close tracks (only to speed up process)
4404 // 3. Calculate cluster overlap ratio - and remove the track if bigger than a threshold
4405 // 4. Delete temporary parameters
4407 const Double_t xref=GetXrow(63); // reference radius -IROC/OROC boundary
4409 const Double_t kCutP1=10; // delta Z cut 10 cm
4410 const Double_t kCutP2=0.15; // delta snp(fi) cut 0.15
4411 const Double_t kCutP3=0.15; // delta tgl(theta) cut 0.15
4412 const Double_t kCutAlpha=0.15; // delta alpha cut
4413 Int_t firstpoint = 0;
4414 Int_t lastpoint = 160;
4416 Int_t nentries = array->GetEntriesFast();
4417 AliExternalTrackParam *params = new AliExternalTrackParam[nentries];
4423 //0. Sort tracks according quality
4424 //1. Propagate the ext. param to reference radius
4425 Int_t nseed = array->GetEntriesFast();
4426 if (nseed<=0) return;
4427 Float_t * quality = new Float_t[nseed];
4428 Int_t * indexes = new Int_t[nseed];
4429 for (Int_t i=0; i<nseed; i++) {
4430 AliTPCseed *pt=(AliTPCseed*)array->UncheckedAt(i);
4435 pt->UpdatePoints(); //select first last max dens points
4436 Float_t * points = pt->GetPoints();
4437 if (points[3]<0.8) quality[i] =-1;
4438 quality[i] = (points[2]-points[0])+pt->GetNumberOfClusters();
4439 //prefer high momenta tracks if overlaps
4440 quality[i] *= TMath::Sqrt(TMath::Abs(pt->Pt())+0.5);
4442 AliTracker::PropagateTrackToBxByBz(&(params[i]),xref,pt->GetMass(),5.,kTRUE);
4443 AliTracker::PropagateTrackToBxByBz(&(params[i]),xref,pt->GetMass(),1.,kTRUE);
4445 TMath::Sort(nseed,quality,indexes);
4447 // 3. Loop over pair of tracks
4449 for (Int_t i0=0; i0<nseed; i0++) {
4450 Int_t index0=indexes[i0];
4451 if (!(array->UncheckedAt(index0))) continue;
4452 AliTPCseed *s1 = (AliTPCseed*)array->UncheckedAt(index0);
4453 if (!s1->IsActive()) continue;
4454 AliExternalTrackParam &par0=params[index0];
4455 for (Int_t i1=i0+1; i1<nseed; i1++) {
4456 Int_t index1=indexes[i1];
4457 if (!(array->UncheckedAt(index1))) continue;
4458 AliTPCseed *s2 = (AliTPCseed*)array->UncheckedAt(index1);
4459 if (!s2->IsActive()) continue;
4460 if (s2->GetKinkIndexes()[0]!=0)
4461 if (s2->GetKinkIndexes()[0] == -s1->GetKinkIndexes()[0]) continue;
4462 AliExternalTrackParam &par1=params[index1];
4463 if (TMath::Abs(par0.GetParameter()[3]-par1.GetParameter()[3])>kCutP3) continue;
4464 if (TMath::Abs(par0.GetParameter()[1]-par1.GetParameter()[1])>kCutP1) continue;
4465 if (TMath::Abs(par0.GetParameter()[2]-par1.GetParameter()[2])>kCutP2) continue;
4466 Double_t dAlpha= TMath::Abs(par0.GetAlpha()-par1.GetAlpha());
4467 if (dAlpha>TMath::Pi()) dAlpha-=TMath::Pi();
4468 if (TMath::Abs(dAlpha)>kCutAlpha) continue;
4473 Int_t firstShared=lastpoint, lastShared=firstpoint;
4474 Int_t firstRow=lastpoint, lastRow=firstpoint;
4476 for (Int_t i=firstpoint;i<lastpoint;i++){
4477 if (s1->GetClusterIndex2(i)>0) nall0++;
4478 if (s2->GetClusterIndex2(i)>0) nall1++;
4479 if (s1->GetClusterIndex2(i)>0 && s2->GetClusterIndex2(i)>0) {
4480 if (i<firstRow) firstRow=i;
4481 if (i>lastRow) lastRow=i;
4483 if ( (s1->GetClusterIndex2(i))==(s2->GetClusterIndex2(i)) && s1->GetClusterIndex2(i)>0) {
4484 if (i<firstShared) firstShared=i;
4485 if (i>lastShared) lastShared=i;
4489 Double_t ratio0 = Float_t(sumShared)/Float_t(TMath::Min(nall0+1,nall1+1));
4490 Double_t ratio1 = Float_t(sumShared)/Float_t(TMath::Max(nall0+1,nall1+1));
4492 if( AliTPCReconstructor::StreamLevel()>1){
4493 TTreeSRedirector &cstream = *fDebugStreamer;
4494 Int_t n0=s1->GetNumberOfClusters();
4495 Int_t n1=s2->GetNumberOfClusters();
4496 Int_t n0F=s1->GetNFoundable();
4497 Int_t n1F=s2->GetNFoundable();
4498 Int_t lab0=s1->GetLabel();
4499 Int_t lab1=s2->GetLabel();
4501 cstream<<"Splitted2"<<
4502 "iter="<<fIteration<<
4503 "lab0="<<lab0<< // MC label if exist
4504 "lab1="<<lab1<< // MC label if exist
4507 "ratio0="<<ratio0<< // shared ratio
4508 "ratio1="<<ratio1<< // shared ratio
4509 "p0.="<<&par0<< // track parameters
4511 "s0.="<<s1<< // full seed
4513 "n0="<<n0<< // number of clusters track 0
4514 "n1="<<n1<< // number of clusters track 1
4515 "nall0="<<nall0<< // number of clusters track 0
4516 "nall1="<<nall1<< // number of clusters track 1
4517 "n0F="<<n0F<< // number of findable
4518 "n1F="<<n1F<< // number of findable
4519 "shared="<<sumShared<< // number of shared clusters
4520 "firstS="<<firstShared<< // first and the last shared row
4521 "lastS="<<lastShared<<
4522 "firstRow="<<firstRow<< // first and the last row with cluster
4523 "lastRow="<<lastRow<< //
4527 // remove track with lower quality
4529 if (ratio0>AliTPCReconstructor::GetRecoParam()->GetCutSharedClusters(0) ||
4530 ratio1>AliTPCReconstructor::GetRecoParam()->GetCutSharedClusters(1)){
4534 delete array->RemoveAt(index1);
4539 // 4. Delete temporary array
4549 void AliTPCtrackerMI::FindCurling(const TObjArray * array, AliESDEvent */*esd*/, Int_t iter)
4552 // find Curling tracks
4553 // Use AliTPCReconstructor::StreamLevel()>1 if you want to tune parameters - cuts
4556 // Algorithm done in 2 phases - because of CPU consumption
4557 // it is n^2 algorithm - for lead-lead 20000x20000 combination are investigated
4558 // see detal in MC part what can be used to cut
4562 const Float_t kMaxC = 400; // maximal curvature to of the track
4563 const Float_t kMaxdTheta = 0.15; // maximal distance in theta
4564 const Float_t kMaxdPhi = 0.15; // maximal distance in phi
4565 const Float_t kPtRatio = 0.3; // ratio between pt
4566 const Float_t kMinDCAR = 2.; // distance to the primary vertex in r - see cpipe cut
4569 // Curling tracks cuts
4572 const Float_t kMaxDeltaRMax = 40; // distance in outer radius
4573 const Float_t kMaxDeltaRMin = 5.; // distance in lower radius - see cpipe cut
4574 const Float_t kMinAngle = 2.9; // angle between tracks
4575 const Float_t kMaxDist = 5; // biggest distance
4577 // The cuts can be tuned using the "MC information stored in Multi tree ==> see FindMultiMC
4580 TCut csign("csign","Tr0.fP[4]*Tr1.fP[4]<0"); //opposite sign
4581 TCut cmax("cmax","abs(Tr0.GetC())>1/400");
4582 TCut cda("cda","sqrt(dtheta^2+dfi^2)<0.15");
4583 TCut ccratio("ccratio","abs((Tr0.fP[4]+Tr1.fP[4])/(abs(Tr0.fP[4])+abs(Tr1.fP[4])))<0.3");
4584 TCut cpipe("cpipe", "min(abs(r0-rc0),abs(r1-rc1))>5");
4586 TCut cdrmax("cdrmax","abs(abs(rc0+r0)-abs(rc1+r1))<40")
4587 TCut cdrmin("cdrmin","abs(abs(rc0+r0)-abs(rc1+r1))<10")
4589 Multi->Draw("dfi","iter==0"+csign+cmax+cda+ccratio); ~94% of curling tracks fulfill
4590 Multi->Draw("min(abs(r0-rc0),abs(r1-rc1))","iter==0&&abs(lab1)==abs(lab0)"+csign+cmax+cda+ccratio+cpipe+cdrmin+cdrmax); //80%
4592 Curling2->Draw("dfi","iter==0&&abs(lab0)==abs(lab1)"+csign+cmax+cdtheta+cdfi+ccratio)
4598 Int_t nentries = array->GetEntriesFast();
4599 AliHelix *helixes = new AliHelix[nentries];
4600 for (Int_t i=0;i<nentries;i++){
4601 AliTPCseed* track = (AliTPCseed*)array->At(i);
4602 if (!track) continue;
4603 track->SetCircular(0);
4604 new (&helixes[i]) AliHelix(*track);
4610 Double_t phase[2][2]={{0,0},{0,0}},radius[2]={0,0};
4616 for (Int_t i0=0;i0<nentries;i0++){
4617 AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
4618 if (!track0) continue;
4619 if (TMath::Abs(track0->GetC())<1/kMaxC) continue;
4620 Float_t xc0 = helixes[i0].GetHelix(6);
4621 Float_t yc0 = helixes[i0].GetHelix(7);
4622 Float_t r0 = helixes[i0].GetHelix(8);
4623 Float_t rc0 = TMath::Sqrt(xc0*xc0+yc0*yc0);
4624 Float_t fi0 = TMath::ATan2(yc0,xc0);
4626 for (Int_t i1=i0+1;i1<nentries;i1++){
4627 AliTPCseed * track1 = (AliTPCseed*)array->At(i1);
4628 if (!track1) continue;
4629 if (TMath::Abs(track1->GetC())<1/kMaxC) continue;
4630 Float_t xc1 = helixes[i1].GetHelix(6);
4631 Float_t yc1 = helixes[i1].GetHelix(7);
4632 Float_t r1 = helixes[i1].GetHelix(8);
4633 Float_t rc1 = TMath::Sqrt(xc1*xc1+yc1*yc1);
4634 Float_t fi1 = TMath::ATan2(yc1,xc1);
4636 Float_t dfi = fi0-fi1;
4639 if (dfi>1.5*TMath::Pi()) dfi-=TMath::Pi(); // take care about edge effect
4640 if (dfi<-1.5*TMath::Pi()) dfi+=TMath::Pi(); //
4641 Float_t dtheta = TMath::Abs(track0->GetTgl()-track1->GetTgl())<TMath::Abs(track0->GetTgl()+track1->GetTgl())? track0->GetTgl()-track1->GetTgl():track0->GetTgl()+track1->GetTgl();
4645 if (track0->GetBConstrain()&&track1->GetBConstrain()) continue; // not constrained
4646 if (track1->GetSigned1Pt()*track0->GetSigned1Pt()>0) continue; // not the same sign
4647 if ( TMath::Abs(track1->GetTgl()+track0->GetTgl())>kMaxdTheta) continue; //distance in the Theta
4648 if ( TMath::Abs(dfi)>kMaxdPhi) continue; //distance in phi
4649 if ( TMath::Sqrt(dfi*dfi+dtheta*dtheta)>kMaxdPhi) continue; //common angular offset
4651 Float_t pt0 = track0->GetSignedPt();
4652 Float_t pt1 = track1->GetSignedPt();
4653 if ((TMath::Abs(pt0+pt1)/(TMath::Abs(pt0)+TMath::Abs(pt1)))>kPtRatio) continue;
4654 if ((iter==1) && TMath::Abs(TMath::Abs(rc0+r0)-TMath::Abs(rc1+r1))>kMaxDeltaRMax) continue;
4655 if ((iter!=1) &&TMath::Abs(TMath::Abs(rc0-r0)-TMath::Abs(rc1-r1))>kMaxDeltaRMin) continue;
4656 if (TMath::Min(TMath::Abs(rc0-r0),TMath::Abs(rc1-r1))<kMinDCAR) continue;
4659 // Now find closest approach
4663 Int_t npoints = helixes[i0].GetRPHIintersections(helixes[i1], phase, radius,10);
4664 if (npoints==0) continue;
4665 helixes[i0].GetClosestPhases(helixes[i1], phase);
4669 Double_t hangles[3];
4670 helixes[i0].Evaluate(phase[0][0],xyz0);
4671 helixes[i1].Evaluate(phase[0][1],xyz1);
4673 helixes[i0].GetAngle(phase[0][0],helixes[i1],phase[0][1],hangles);
4674 Double_t deltah[2],deltabest;
4675 if (TMath::Abs(hangles[2])<kMinAngle) continue;
4679 helixes[i0].ParabolicDCA(helixes[i1],phase[0][0],phase[0][1],radius[0],deltah[0],2);
4681 helixes[i0].ParabolicDCA(helixes[i1],phase[1][0],phase[1][1],radius[1],deltah[1],2);
4682 if (deltah[1]<deltah[0]) ibest=1;
4684 deltabest = TMath::Sqrt(deltah[ibest]);
4685 helixes[i0].Evaluate(phase[ibest][0],xyz0);
4686 helixes[i1].Evaluate(phase[ibest][1],xyz1);
4687 helixes[i0].GetAngle(phase[ibest][0],helixes[i1],phase[ibest][1],hangles);
4688 Double_t radiusbest = TMath::Sqrt(radius[ibest]);
4690 if (deltabest>kMaxDist) continue;
4691 // if (mindcar+mindcaz<40 && (TMath::Abs(hangles[2])<kMinAngle ||deltabest>3)) continue;
4692 Bool_t sign =kFALSE;
4693 if (hangles[2]>kMinAngle) sign =kTRUE;
4696 // circular[i0] = kTRUE;
4697 // circular[i1] = kTRUE;
4698 if (track0->OneOverPt()<track1->OneOverPt()){
4699 track0->SetCircular(track0->GetCircular()+1);
4700 track1->SetCircular(track1->GetCircular()+2);
4703 track1->SetCircular(track1->GetCircular()+1);
4704 track0->SetCircular(track0->GetCircular()+2);
4707 if (AliTPCReconstructor::StreamLevel()>2){
4709 //debug stream to tune "fine" cuts
4710 Int_t lab0=track0->GetLabel();
4711 Int_t lab1=track1->GetLabel();
4712 TTreeSRedirector &cstream = *fDebugStreamer;
4713 cstream<<"Curling2"<<
4729 "npoints="<<npoints<<
4730 "hangles0="<<hangles[0]<<
4731 "hangles1="<<hangles[1]<<
4732 "hangles2="<<hangles[2]<<
4735 "radius="<<radiusbest<<
4736 "deltabest="<<deltabest<<
4737 "phase0="<<phase[ibest][0]<<
4738 "phase1="<<phase[ibest][1]<<
4746 if (AliTPCReconstructor::StreamLevel()>1) {
4747 AliInfo("Time for curling tracks removal");
4756 void AliTPCtrackerMI::FindKinks(TObjArray * array, AliESDEvent *esd)
4763 TObjArray *kinks= new TObjArray(10000);
4764 // TObjArray *v0s= new TObjArray(10000);
4765 Int_t nentries = array->GetEntriesFast();
4766 AliHelix *helixes = new AliHelix[nentries];
4767 Int_t *sign = new Int_t[nentries];
4768 Int_t *nclusters = new Int_t[nentries];
4769 Float_t *alpha = new Float_t[nentries];
4770 AliKink *kink = new AliKink();
4771 Int_t * usage = new Int_t[nentries];
4772 Float_t *zm = new Float_t[nentries];
4773 Float_t *z0 = new Float_t[nentries];
4774 Float_t *fim = new Float_t[nentries];
4775 Float_t *shared = new Float_t[nentries];
4776 Bool_t *circular = new Bool_t[nentries];
4777 Float_t *dca = new Float_t[nentries];
4778 //const AliESDVertex * primvertex = esd->GetVertex();
4780 // nentries = array->GetEntriesFast();
4785 for (Int_t i=0;i<nentries;i++){
4788 AliTPCseed* track = (AliTPCseed*)array->At(i);
4789 if (!track) continue;
4790 track->SetCircular(0);
4792 track->UpdatePoints();
4793 if (( track->GetPoints()[2]- track->GetPoints()[0])>5 && track->GetPoints()[3]>0.8){
4795 nclusters[i]=track->GetNumberOfClusters();
4796 alpha[i] = track->GetAlpha();
4797 new (&helixes[i]) AliHelix(*track);
4799 helixes[i].Evaluate(0,xyz);
4800 sign[i] = (track->GetC()>0) ? -1:1;
4803 if (track->GetProlongation(x,y,z)){
4805 fim[i] = alpha[i]+TMath::ATan2(y,x);
4808 zm[i] = track->GetZ();
4812 circular[i]= kFALSE;
4813 if (track->GetProlongation(0,y,z)) z0[i] = z;
4814 dca[i] = track->GetD(0,0);
4820 Int_t ncandidates =0;
4823 Double_t phase[2][2]={{0,0},{0,0}},radius[2]={0,0};
4826 // Find circling track
4828 for (Int_t i0=0;i0<nentries;i0++){
4829 AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
4830 if (!track0) continue;
4831 if (track0->GetNumberOfClusters()<40) continue;
4832 if (TMath::Abs(1./track0->GetC())>200) continue;
4833 for (Int_t i1=i0+1;i1<nentries;i1++){
4834 AliTPCseed * track1 = (AliTPCseed*)array->At(i1);
4835 if (!track1) continue;
4836 if (track1->GetNumberOfClusters()<40) continue;
4837 if ( TMath::Abs(track1->GetTgl()+track0->GetTgl())>0.1) continue;
4838 if (track0->GetBConstrain()&&track1->GetBConstrain()) continue;
4839 if (TMath::Abs(1./track1->GetC())>200) continue;
4840 if (track1->GetSigned1Pt()*track0->GetSigned1Pt()>0) continue;
4841 if (track1->GetTgl()*track0->GetTgl()>0) continue;
4842 if (TMath::Max(TMath::Abs(1./track0->GetC()),TMath::Abs(1./track1->GetC()))>190) continue;
4843 if (track0->GetBConstrain()&&track1->OneOverPt()<track0->OneOverPt()) continue; //returning - lower momenta
4844 if (track1->GetBConstrain()&&track0->OneOverPt()<track1->OneOverPt()) continue; //returning - lower momenta
4846 Float_t mindcar = TMath::Min(TMath::Abs(dca[i0]),TMath::Abs(dca[i1]));
4847 if (mindcar<5) continue;
4848 Float_t mindcaz = TMath::Min(TMath::Abs(z0[i0]-GetZ()),TMath::Abs(z0[i1]-GetZ()));
4849 if (mindcaz<5) continue;
4850 if (mindcar+mindcaz<20) continue;
4853 Float_t xc0 = helixes[i0].GetHelix(6);
4854 Float_t yc0 = helixes[i0].GetHelix(7);
4855 Float_t r0 = helixes[i0].GetHelix(8);
4856 Float_t xc1 = helixes[i1].GetHelix(6);
4857 Float_t yc1 = helixes[i1].GetHelix(7);
4858 Float_t r1 = helixes[i1].GetHelix(8);
4860 Float_t rmean = (r0+r1)*0.5;
4861 Float_t delta =TMath::Sqrt((xc1-xc0)*(xc1-xc0)+(yc1-yc0)*(yc1-yc0));
4862 //if (delta>30) continue;
4863 if (delta>rmean*0.25) continue;
4864 if (TMath::Abs(r0-r1)/rmean>0.3) continue;
4866 Int_t npoints = helixes[i0].GetRPHIintersections(helixes[i1], phase, radius,10);
4867 if (npoints==0) continue;
4868 helixes[i0].GetClosestPhases(helixes[i1], phase);
4872 Double_t hangles[3];
4873 helixes[i0].Evaluate(phase[0][0],xyz0);
4874 helixes[i1].Evaluate(phase[0][1],xyz1);
4876 helixes[i0].GetAngle(phase[0][0],helixes[i1],phase[0][1],hangles);
4877 Double_t deltah[2],deltabest;
4878 if (hangles[2]<2.8) continue;
4881 helixes[i0].ParabolicDCA(helixes[i1],phase[0][0],phase[0][1],radius[0],deltah[0],2);
4883 helixes[i0].ParabolicDCA(helixes[i1],phase[1][0],phase[1][1],radius[1],deltah[1],2);
4884 if (deltah[1]<deltah[0]) ibest=1;
4886 deltabest = TMath::Sqrt(deltah[ibest]);
4887 helixes[i0].Evaluate(phase[ibest][0],xyz0);
4888 helixes[i1].Evaluate(phase[ibest][1],xyz1);
4889 helixes[i0].GetAngle(phase[ibest][0],helixes[i1],phase[ibest][1],hangles);
4890 Double_t radiusbest = TMath::Sqrt(radius[ibest]);
4892 if (deltabest>6) continue;
4893 if (mindcar+mindcaz<40 && (hangles[2]<3.12||deltabest>3)) continue;
4894 Bool_t lsign =kFALSE;
4895 if (hangles[2]>3.06) lsign =kTRUE;
4898 circular[i0] = kTRUE;
4899 circular[i1] = kTRUE;
4900 if (track0->OneOverPt()<track1->OneOverPt()){
4901 track0->SetCircular(track0->GetCircular()+1);
4902 track1->SetCircular(track1->GetCircular()+2);
4905 track1->SetCircular(track1->GetCircular()+1);
4906 track0->SetCircular(track0->GetCircular()+2);
4909 if (lsign&&AliTPCReconstructor::StreamLevel()>1){
4911 Int_t lab0=track0->GetLabel();
4912 Int_t lab1=track1->GetLabel();
4913 TTreeSRedirector &cstream = *fDebugStreamer;
4914 cstream<<"Curling"<<
4921 "mindcar="<<mindcar<<
4922 "mindcaz="<<mindcaz<<
4925 "npoints="<<npoints<<
4926 "hangles0="<<hangles[0]<<
4927 "hangles2="<<hangles[2]<<
4932 "radius="<<radiusbest<<
4933 "deltabest="<<deltabest<<
4934 "phase0="<<phase[ibest][0]<<
4935 "phase1="<<phase[ibest][1]<<
4945 for (Int_t i =0;i<nentries;i++){
4946 if (sign[i]==0) continue;
4947 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
4954 Double_t cradius0 = 40*40;
4955 Double_t cradius1 = 270*270;
4958 Double_t cdist3=0.55;
4959 for (Int_t j =i+1;j<nentries;j++){
4961 if (sign[j]*sign[i]<1) continue;
4962 if ( (nclusters[i]+nclusters[j])>200) continue;
4963 if ( (nclusters[i]+nclusters[j])<80) continue;
4964 if ( TMath::Abs(zm[i]-zm[j])>60.) continue;
4965 if ( TMath::Abs(fim[i]-fim[j])>0.6 && TMath::Abs(fim[i]-fim[j])<5.7 ) continue;
4966 //AliTPCseed * track1 = (AliTPCseed*)array->At(j); Double_t phase[2][2],radius[2];
4967 Int_t npoints = helixes[i].GetRPHIintersections(helixes[j], phase, radius,20);
4968 if (npoints<1) continue;
4971 if (radius[0]<cradius0||radius[0]>cradius1) continue;
4974 if ( (radius[0]<cradius0||radius[0]>cradius1) && (radius[1]<cradius0||radius[1]>cradius1) ) continue;
4977 Double_t delta1=10000,delta2=10000;
4978 // cuts on the intersection radius
4979 helixes[i].LinearDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta1);
4980 if (radius[0]<20&&delta1<1) continue; //intersection at vertex
4981 if (radius[0]<10&&delta1<3) continue; //intersection at vertex
4983 helixes[i].LinearDCA(helixes[j],phase[1][0],phase[1][1],radius[1],delta2);
4984 if (radius[1]<20&&delta2<1) continue; //intersection at vertex
4985 if (radius[1]<10&&delta2<3) continue; //intersection at vertex
4988 Double_t distance1 = TMath::Min(delta1,delta2);
4989 if (distance1>cdist1) continue; // cut on DCA linear approximation
4991 npoints = helixes[i].GetRPHIintersections(helixes[j], phase, radius,20);
4992 helixes[i].ParabolicDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta1);
4993 if (radius[0]<20&&delta1<1) continue; //intersection at vertex
4994 if (radius[0]<10&&delta1<3) continue; //intersection at vertex
4997 helixes[i].ParabolicDCA(helixes[j],phase[1][0],phase[1][1],radius[1],delta2);
4998 if (radius[1]<20&&delta2<1) continue; //intersection at vertex
4999 if (radius[1]<10&&delta2<3) continue; //intersection at vertex
5001 distance1 = TMath::Min(delta1,delta2);
5004 rkink = TMath::Sqrt(radius[0]);
5007 rkink = TMath::Sqrt(radius[1]);
5009 if (distance1>cdist2) continue;
5012 AliTPCseed * track1 = (AliTPCseed*)array->At(j);
5015 Int_t row0 = GetRowNumber(rkink);
5016 if (row0<10) continue;
5017 if (row0>150) continue;
5020 Float_t dens00=-1,dens01=-1;
5021 Float_t dens10=-1,dens11=-1;
5023 Int_t found,foundable,ishared;
5024 track0->GetClusterStatistic(0,row0-5, found, foundable,ishared,kFALSE);
5025 if (foundable>5) dens00 = Float_t(found)/Float_t(foundable);
5026 track0->GetClusterStatistic(row0+5,155, found, foundable,ishared,kFALSE);
5027 if (foundable>5) dens01 = Float_t(found)/Float_t(foundable);
5029 track1->GetClusterStatistic(0,row0-5, found, foundable,ishared,kFALSE);
5030 if (foundable>10) dens10 = Float_t(found)/Float_t(foundable);
5031 track1->GetClusterStatistic(row0+5,155, found, foundable,ishared,kFALSE);
5032 if (foundable>10) dens11 = Float_t(found)/Float_t(foundable);
5034 if (dens00<dens10 && dens01<dens11) continue;
5035 if (dens00>dens10 && dens01>dens11) continue;
5036 if (TMath::Max(dens00,dens10)<0.1) continue;
5037 if (TMath::Max(dens01,dens11)<0.3) continue;
5039 if (TMath::Min(dens00,dens10)>0.6) continue;
5040 if (TMath::Min(dens01,dens11)>0.6) continue;
5043 AliTPCseed * ktrack0, *ktrack1;
5052 if (TMath::Abs(ktrack0->GetC())>5) continue; // cut on the curvature for mother particle
5053 AliExternalTrackParam paramm(*ktrack0);
5054 AliExternalTrackParam paramd(*ktrack1);
5055 if (row0>60&&ktrack1->GetReference().GetX()>90.)new (¶md) AliExternalTrackParam(ktrack1->GetReference());
5058 kink->SetMother(paramm);
5059 kink->SetDaughter(paramd);
5062 Float_t x[3] = { kink->GetPosition()[0],kink->GetPosition()[1],kink->GetPosition()[2]};
5064 fkParam->Transform0to1(x,index);
5065 fkParam->Transform1to2(x,index);
5066 row0 = GetRowNumber(x[0]);
5068 if (kink->GetR()<100) continue;
5069 if (kink->GetR()>240) continue;
5070 if (kink->GetPosition()[2]/kink->GetR()>AliTPCReconstructor::GetCtgRange()) continue; //out of fiducial volume
5071 if (kink->GetDistance()>cdist3) continue;
5072 Float_t dird = kink->GetDaughterP()[0]*kink->GetPosition()[0]+kink->GetDaughterP()[1]*kink->GetPosition()[1]; // rough direction estimate
5073 if (dird<0) continue;
5075 Float_t dirm = kink->GetMotherP()[0]*kink->GetPosition()[0]+kink->GetMotherP()[1]*kink->GetPosition()[1]; // rough direction estimate
5076 if (dirm<0) continue;
5077 Float_t mpt = TMath::Sqrt(kink->GetMotherP()[0]*kink->GetMotherP()[0]+kink->GetMotherP()[1]*kink->GetMotherP()[1]);
5078 if (mpt<0.2) continue;
5081 //for high momenta momentum not defined well in first iteration
5082 Double_t qt = TMath::Sin(kink->GetAngle(2))*ktrack1->GetP();
5083 if (qt>0.35) continue;
5086 kink->SetLabel(CookLabel(ktrack0,0.4,0,row0),0);
5087 kink->SetLabel(CookLabel(ktrack1,0.4,row0,160),1);
5089 kink->SetTPCDensity(dens00,0,0);
5090 kink->SetTPCDensity(dens01,0,1);
5091 kink->SetTPCDensity(dens10,1,0);
5092 kink->SetTPCDensity(dens11,1,1);
5093 kink->SetIndex(i,0);
5094 kink->SetIndex(j,1);
5097 kink->SetTPCDensity(dens10,0,0);
5098 kink->SetTPCDensity(dens11,0,1);
5099 kink->SetTPCDensity(dens00,1,0);
5100 kink->SetTPCDensity(dens01,1,1);
5101 kink->SetIndex(j,0);
5102 kink->SetIndex(i,1);
5105 if (mpt<1||kink->GetAngle(2)>0.1){
5106 // angle and densities not defined yet
5107 if (kink->GetTPCDensityFactor()<0.8) continue;
5108 if ((2-kink->GetTPCDensityFactor())*kink->GetDistance() >0.25) continue;
5109 if (kink->GetAngle(2)*ktrack0->GetP()<0.003) continue; //too small angle
5110 if (kink->GetAngle(2)>0.2&&kink->GetTPCDensityFactor()<1.15) continue;
5111 if (kink->GetAngle(2)>0.2&&kink->GetTPCDensity(0,1)>0.05) continue;
5113 Float_t criticalangle = track0->GetSigmaSnp2()+track0->GetSigmaTgl2();
5114 criticalangle+= track1->GetSigmaSnp2()+track1->GetSigmaTgl2();
5115 criticalangle= 3*TMath::Sqrt(criticalangle);
5116 if (criticalangle>0.02) criticalangle=0.02;
5117 if (kink->GetAngle(2)<criticalangle) continue;
5120 Int_t drow = Int_t(2.+0.5/(0.05+kink->GetAngle(2))); // overlap region defined
5121 Float_t shapesum =0;
5123 for ( Int_t row = row0-drow; row<row0+drow;row++){
5124 if (row<0) continue;
5125 if (row>155) continue;
5126 if (ktrack0->GetClusterPointer(row)){
5127 AliTPCTrackerPoint *point =ktrack0->GetTrackPoint(row);
5128 shapesum+=point->GetSigmaY()+point->GetSigmaZ();
5131 if (ktrack1->GetClusterPointer(row)){
5132 AliTPCTrackerPoint *point =ktrack1->GetTrackPoint(row);
5133 shapesum+=point->GetSigmaY()+point->GetSigmaZ();
5138 kink->SetShapeFactor(-1.);
5141 kink->SetShapeFactor(shapesum/sum);
5143 // esd->AddKink(kink);
5145 // kink->SetMother(paramm);
5146 //kink->SetDaughter(paramd);
5148 Double_t chi2P2 = paramm.GetParameter()[2]-paramd.GetParameter()[2];
5150 chi2P2/=paramm.GetCovariance()[5]+paramd.GetCovariance()[5];
5151 Double_t chi2P3 = paramm.GetParameter()[3]-paramd.GetParameter()[3];
5153 chi2P3/=paramm.GetCovariance()[9]+paramd.GetCovariance()[9];
5155 if (AliTPCReconstructor::StreamLevel()>1) {
5156 (*fDebugStreamer)<<"kinkLpt"<<
5164 if ( chi2P2+chi2P3<AliTPCReconstructor::GetRecoParam()->GetKinkAngleCutChi2(0)){
5168 kinks->AddLast(kink);
5174 // sort the kinks according quality - and refit them towards vertex
5176 Int_t nkinks = kinks->GetEntriesFast();
5177 Float_t *quality = new Float_t[nkinks];
5178 Int_t *indexes = new Int_t[nkinks];
5179 AliTPCseed *mothers = new AliTPCseed[nkinks];
5180 AliTPCseed *daughters = new AliTPCseed[nkinks];
5183 for (Int_t i=0;i<nkinks;i++){
5185 AliKink *kinkl = (AliKink*)kinks->At(i);
5187 // refit kinks towards vertex
5189 Int_t index0 = kinkl->GetIndex(0);
5190 Int_t index1 = kinkl->GetIndex(1);
5191 AliTPCseed * ktrack0 = (AliTPCseed*)array->At(index0);
5192 AliTPCseed * ktrack1 = (AliTPCseed*)array->At(index1);
5194 Int_t sumn=ktrack0->GetNumberOfClusters()+ktrack1->GetNumberOfClusters();
5196 // Refit Kink under if too small angle
5198 if (kinkl->GetAngle(2)<0.05){
5199 kinkl->SetTPCRow0(GetRowNumber(kinkl->GetR()));
5200 Int_t row0 = kinkl->GetTPCRow0();
5201 Int_t drow = Int_t(2.+0.5/(0.05+kinkl->GetAngle(2)));
5204 Int_t last = row0-drow;
5205 if (last<40) last=40;
5206 if (last<ktrack0->GetFirstPoint()+25) last = ktrack0->GetFirstPoint()+25;
5207 AliTPCseed* seed0 = ReSeed(ktrack0,last,kFALSE);
5210 Int_t first = row0+drow;
5211 if (first>130) first=130;
5212 if (first>ktrack1->GetLastPoint()-25) first = TMath::Max(ktrack1->GetLastPoint()-25,30);
5213 AliTPCseed* seed1 = ReSeed(ktrack1,first,kTRUE);
5215 if (seed0 && seed1){
5216 kinkl->SetStatus(1,8);
5217 if (RefitKink(*seed0,*seed1,*kinkl)) kinkl->SetStatus(1,9);
5218 row0 = GetRowNumber(kinkl->GetR());
5219 sumn = seed0->GetNumberOfClusters()+seed1->GetNumberOfClusters();
5220 mothers[i] = *seed0;
5221 daughters[i] = *seed1;
5224 delete kinks->RemoveAt(i);
5225 if (seed0) delete seed0;
5226 if (seed1) delete seed1;
5229 if (kinkl->GetDistance()>0.5 || kinkl->GetR()<110 || kinkl->GetR()>240) {
5230 delete kinks->RemoveAt(i);
5231 if (seed0) delete seed0;
5232 if (seed1) delete seed1;
5240 if (kinkl) quality[i] = 160*((0.1+kinkl->GetDistance())*(2.-kinkl->GetTPCDensityFactor()))/(sumn+40.); //the longest -clossest will win
5242 TMath::Sort(nkinks,quality,indexes,kFALSE);
5244 //remove double find kinks
5246 for (Int_t ikink0=1;ikink0<nkinks;ikink0++){
5247 AliKink * kink0 = (AliKink*) kinks->At(indexes[ikink0]);
5248 if (!kink0) continue;
5250 for (Int_t ikink1=0;ikink1<ikink0;ikink1++){
5251 kink0 = (AliKink*) kinks->At(indexes[ikink0]);
5252 if (!kink0) continue;
5253 AliKink * kink1 = (AliKink*) kinks->At(indexes[ikink1]);
5254 if (!kink1) continue;
5255 // if not close kink continue
5256 if (TMath::Abs(kink1->GetPosition()[2]-kink0->GetPosition()[2])>10) continue;
5257 if (TMath::Abs(kink1->GetPosition()[1]-kink0->GetPosition()[1])>10) continue;
5258 if (TMath::Abs(kink1->GetPosition()[0]-kink0->GetPosition()[0])>10) continue;
5260 AliTPCseed &mother0 = mothers[indexes[ikink0]];
5261 AliTPCseed &daughter0 = daughters[indexes[ikink0]];
5262 AliTPCseed &mother1 = mothers[indexes[ikink1]];
5263 AliTPCseed &daughter1 = daughters[indexes[ikink1]];
5264 Int_t row0 = (kink0->GetTPCRow0()+kink1->GetTPCRow0())/2;
5273 for (Int_t i=0;i<row0;i++){
5274 if (mother0.GetClusterIndex(i)>0 && mother1.GetClusterIndex(i)>0){
5277 if (mother0.GetClusterIndex(i)==mother1.GetClusterIndex(i)){
5284 for (Int_t i=row0;i<158;i++){
5285 if (daughter0.GetClusterIndex(i)>0 && daughter0.GetClusterIndex(i)>0){
5288 if (mother0.GetClusterIndex(i)==mother1.GetClusterIndex(i)){
5294 Float_t ratio = Float_t(same+1)/Float_t(both+1);
5295 Float_t ratiom = Float_t(samem+1)/Float_t(bothm+1);
5296 Float_t ratiod = Float_t(samed+1)/Float_t(bothd+1);
5297 if (ratio>0.3 && ratiom>0.5 &&ratiod>0.5) {
5298 Int_t sum0 = mother0.GetNumberOfClusters()+daughter0.GetNumberOfClusters();
5299 Int_t sum1 = mother1.GetNumberOfClusters()+daughter1.GetNumberOfClusters();
5301 shared[kink0->GetIndex(0)]= kTRUE;
5302 shared[kink0->GetIndex(1)]= kTRUE;
5303 delete kinks->RemoveAt(indexes[ikink0]);
5307 shared[kink1->GetIndex(0)]= kTRUE;
5308 shared[kink1->GetIndex(1)]= kTRUE;
5309 delete kinks->RemoveAt(indexes[ikink1]);
5316 for (Int_t i=0;i<nkinks;i++){
5317 AliKink * kinkl = (AliKink*) kinks->At(indexes[i]);
5318 if (!kinkl) continue;
5319 kinkl->SetTPCRow0(GetRowNumber(kinkl->GetR()));
5320 Int_t index0 = kinkl->GetIndex(0);
5321 Int_t index1 = kinkl->GetIndex(1);
5322 if (circular[index0]||(circular[index1]&&kinkl->GetDistance()>0.2)) continue;
5323 kinkl->SetMultiple(usage[index0],0);
5324 kinkl->SetMultiple(usage[index1],1);
5325 if (kinkl->GetMultiple()[0]+kinkl->GetMultiple()[1]>2) continue;
5326 if (kinkl->GetMultiple()[0]+kinkl->GetMultiple()[1]>0 && quality[indexes[i]]>0.2) continue;
5327 if (kinkl->GetMultiple()[0]+kinkl->GetMultiple()[1]>0 && kinkl->GetDistance()>0.2) continue;
5328 if (circular[index0]||(circular[index1]&&kinkl->GetDistance()>0.1)) continue;
5330 AliTPCseed * ktrack0 = (AliTPCseed*)array->At(index0);
5331 AliTPCseed * ktrack1 = (AliTPCseed*)array->At(index1);
5332 if (!ktrack0 || !ktrack1) continue;
5333 Int_t index = esd->AddKink(kinkl);
5336 if ( ktrack0->GetKinkIndex(0)==0 && ktrack1->GetKinkIndex(0)==0) { //best kink
5337 if (mothers[indexes[i]].GetNumberOfClusters()>20 && daughters[indexes[i]].GetNumberOfClusters()>20 && (mothers[indexes[i]].GetNumberOfClusters()+daughters[indexes[i]].GetNumberOfClusters())>100){
5338 *ktrack0 = mothers[indexes[i]];
5339 *ktrack1 = daughters[indexes[i]];
5343 ktrack0->SetKinkIndex(usage[index0],-(index+1));
5344 ktrack1->SetKinkIndex(usage[index1], (index+1));
5349 // Remove tracks corresponding to shared kink's
5351 for (Int_t i=0;i<nentries;i++){
5352 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5353 if (!track0) continue;
5354 if (track0->GetKinkIndex(0)!=0) continue;
5355 if (shared[i]) delete array->RemoveAt(i);
5360 RemoveUsed2(array,0.5,0.4,30);
5362 for (Int_t i=0;i<nentries;i++){
5363 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5364 if (!track0) continue;
5365 track0->CookdEdx(0.02,0.6);
5369 for (Int_t i=0;i<nentries;i++){
5370 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5371 if (!track0) continue;
5372 if (track0->Pt()<1.4) continue;
5373 //remove double high momenta tracks - overlapped with kink candidates
5376 for (Int_t icl=track0->GetFirstPoint();icl<track0->GetLastPoint(); icl++){
5377 if (track0->GetClusterPointer(icl)!=0){
5379 if (track0->GetClusterPointer(icl)->IsUsed(10)) ishared++;
5382 if (Float_t(ishared+1)/Float_t(all+1)>0.5) {
5383 delete array->RemoveAt(i);
5387 if (track0->GetKinkIndex(0)!=0) continue;
5388 if (track0->GetNumberOfClusters()<80) continue;
5390 AliTPCseed *pmother = new AliTPCseed();
5391 AliTPCseed *pdaughter = new AliTPCseed();
5392 AliKink *pkink = new AliKink;
5394 AliTPCseed & mother = *pmother;
5395 AliTPCseed & daughter = *pdaughter;
5396 AliKink & kinkl = *pkink;
5397 if (CheckKinkPoint(track0,mother,daughter, kinkl)){
5398 if (mother.GetNumberOfClusters()<30||daughter.GetNumberOfClusters()<20) {
5402 continue; //too short tracks
5404 if (mother.Pt()<1.4) {
5410 Int_t row0= kinkl.GetTPCRow0();
5411 if (kinkl.GetDistance()>0.5 || kinkl.GetR()<110. || kinkl.GetR()>240.) {
5418 Int_t index = esd->AddKink(&kinkl);
5419 mother.SetKinkIndex(0,-(index+1));
5420 daughter.SetKinkIndex(0,index+1);
5421 if (mother.GetNumberOfClusters()>50) {
5422 delete array->RemoveAt(i);
5423 array->AddAt(new AliTPCseed(mother),i);
5426 array->AddLast(new AliTPCseed(mother));
5428 array->AddLast(new AliTPCseed(daughter));
5429 for (Int_t icl=0;icl<row0;icl++) {
5430 if (mother.GetClusterPointer(icl)) mother.GetClusterPointer(icl)->Use(20);
5433 for (Int_t icl=row0;icl<158;icl++) {
5434 if (daughter.GetClusterPointer(icl)) daughter.GetClusterPointer(icl)->Use(20);
5443 delete [] daughters;
5465 AliInfo(Form("Ncandidates=\t%d\t%d\t%d\t%d\n",esd->GetNumberOfKinks(),ncandidates,ntracks,nall));
5470 Int_t AliTPCtrackerMI::RefitKink(AliTPCseed &mother, AliTPCseed &daughter, const AliESDkink &knk)
5473 // refit kink towards to the vertex
5476 AliKink &kink=(AliKink &)knk;
5478 Int_t row0 = GetRowNumber(kink.GetR());
5479 FollowProlongation(mother,0);
5480 mother.Reset(kFALSE);
5482 FollowProlongation(daughter,row0);
5483 daughter.Reset(kFALSE);
5484 FollowBackProlongation(daughter,158);
5485 daughter.Reset(kFALSE);
5486 Int_t first = TMath::Max(row0-20,30);
5487 Int_t last = TMath::Min(row0+20,140);
5489 const Int_t kNdiv =5;
5490 AliTPCseed param0[kNdiv]; // parameters along the track
5491 AliTPCseed param1[kNdiv]; // parameters along the track
5492 AliKink kinks[kNdiv]; // corresponding kink parameters
5495 for (Int_t irow=0; irow<kNdiv;irow++){
5496 rows[irow] = first +((last-first)*irow)/(kNdiv-1);
5498 // store parameters along the track
5500 for (Int_t irow=0;irow<kNdiv;irow++){
5501 FollowBackProlongation(mother, rows[irow]);
5502 FollowProlongation(daughter,rows[kNdiv-1-irow]);
5503 param0[irow] = mother;
5504 param1[kNdiv-1-irow] = daughter;
5508 for (Int_t irow=0; irow<kNdiv-1;irow++){
5509 if (param0[irow].GetNumberOfClusters()<kNdiv||param1[irow].GetNumberOfClusters()<kNdiv) continue;
5510 kinks[irow].SetMother(param0[irow]);
5511 kinks[irow].SetDaughter(param1[irow]);
5512 kinks[irow].Update();
5515 // choose kink with best "quality"
5517 Double_t mindist = 10000;
5518 for (Int_t irow=0;irow<kNdiv;irow++){
5519 if (param0[irow].GetNumberOfClusters()<20||param1[irow].GetNumberOfClusters()<20) continue;
5520 if (TMath::Abs(kinks[irow].GetR())>240.) continue;
5521 if (TMath::Abs(kinks[irow].GetR())<100.) continue;
5523 Float_t normdist = TMath::Abs(param0[irow].GetX()-kinks[irow].GetR())*(0.1+kink.GetDistance());
5524 normdist/= (param0[irow].GetNumberOfClusters()+param1[irow].GetNumberOfClusters()+40.);
5525 if (normdist < mindist){
5531 if (index==-1) return 0;
5534 param0[index].Reset(kTRUE);
5535 FollowProlongation(param0[index],0);
5537 mother = param0[index];
5538 daughter = param1[index]; // daughter in vertex
5540 kink.SetMother(mother);
5541 kink.SetDaughter(daughter);
5543 kink.SetTPCRow0(GetRowNumber(kink.GetR()));
5544 kink.SetTPCncls(param0[index].GetNumberOfClusters(),0);
5545 kink.SetTPCncls(param1[index].GetNumberOfClusters(),1);
5546 kink.SetLabel(CookLabel(&mother,0.4, 0,kink.GetTPCRow0()),0);
5547 kink.SetLabel(CookLabel(&daughter,0.4, kink.GetTPCRow0(),160),1);
5548 mother.SetLabel(kink.GetLabel(0));
5549 daughter.SetLabel(kink.GetLabel(1));
5555 void AliTPCtrackerMI::UpdateKinkQualityM(AliTPCseed * seed){
5557 // update Kink quality information for mother after back propagation
5559 if (seed->GetKinkIndex(0)>=0) return;
5560 for (Int_t ikink=0;ikink<3;ikink++){
5561 Int_t index = seed->GetKinkIndex(ikink);
5562 if (index>=0) break;
5563 index = TMath::Abs(index)-1;
5564 AliESDkink * kink = fEvent->GetKink(index);
5565 kink->SetTPCDensity(-1,0,0);
5566 kink->SetTPCDensity(1,0,1);
5568 Int_t row0 = kink->GetTPCRow0() - 2 - Int_t( 0.5/ (0.05+kink->GetAngle(2)));
5569 if (row0<15) row0=15;
5571 Int_t row1 = kink->GetTPCRow0() + 2 + Int_t( 0.5/ (0.05+kink->GetAngle(2)));
5572 if (row1>145) row1=145;
5574 Int_t found,foundable,shared;
5575 seed->GetClusterStatistic(0,row0, found, foundable,shared,kFALSE);
5576 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),0,0);
5577 seed->GetClusterStatistic(row1,155, found, foundable,shared,kFALSE);
5578 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),0,1);
5583 void AliTPCtrackerMI::UpdateKinkQualityD(AliTPCseed * seed){
5585 // update Kink quality information for daughter after refit
5587 if (seed->GetKinkIndex(0)<=0) return;
5588 for (Int_t ikink=0;ikink<3;ikink++){
5589 Int_t index = seed->GetKinkIndex(ikink);
5590 if (index<=0) break;
5591 index = TMath::Abs(index)-1;
5592 AliESDkink * kink = fEvent->GetKink(index);
5593 kink->SetTPCDensity(-1,1,0);
5594 kink->SetTPCDensity(-1,1,1);
5596 Int_t row0 = kink->GetTPCRow0() -2 - Int_t( 0.5/ (0.05+kink->GetAngle(2)));
5597 if (row0<15) row0=15;
5599 Int_t row1 = kink->GetTPCRow0() +2 + Int_t( 0.5/ (0.05+kink->GetAngle(2)));
5600 if (row1>145) row1=145;
5602 Int_t found,foundable,shared;
5603 seed->GetClusterStatistic(0,row0, found, foundable,shared,kFALSE);
5604 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),1,0);
5605 seed->GetClusterStatistic(row1,155, found, foundable,shared,kFALSE);
5606 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),1,1);
5612 Int_t AliTPCtrackerMI::CheckKinkPoint(AliTPCseed*seed,AliTPCseed &mother, AliTPCseed &daughter, const AliESDkink &knk)
5615 // check kink point for given track
5616 // if return value=0 kink point not found
5617 // otherwise seed0 correspond to mother particle
5618 // seed1 correspond to daughter particle
5619 // kink parameter of kink point
5620 AliKink &kink=(AliKink &)knk;
5622 Int_t middlerow = (seed->GetFirstPoint()+seed->GetLastPoint())/2;
5623 Int_t first = seed->GetFirstPoint();
5624 Int_t last = seed->GetLastPoint();
5625 if (last-first<20) return 0; // shortest length - 2*30 = 60 pad-rows
5628 AliTPCseed *seed1 = ReSeed(seed,middlerow+20, kTRUE); //middle of chamber
5629 if (!seed1) return 0;
5630 FollowProlongation(*seed1,seed->GetLastPoint()-20);
5631 seed1->Reset(kTRUE);
5632 FollowProlongation(*seed1,158);
5633 seed1->Reset(kTRUE);
5634 last = seed1->GetLastPoint();
5636 AliTPCseed *seed0 = new AliTPCseed(*seed);
5637 seed0->Reset(kFALSE);
5640 AliTPCseed param0[20]; // parameters along the track
5641 AliTPCseed param1[20]; // parameters along the track
5642 AliKink kinks[20]; // corresponding kink parameters
5644 for (Int_t irow=0; irow<20;irow++){
5645 rows[irow] = first +((last-first)*irow)/19;
5647 // store parameters along the track
5649 for (Int_t irow=0;irow<20;irow++){
5650 FollowBackProlongation(*seed0, rows[irow]);
5651 FollowProlongation(*seed1,rows[19-irow]);
5652 param0[irow] = *seed0;
5653 param1[19-irow] = *seed1;
5657 for (Int_t irow=0; irow<19;irow++){
5658 kinks[irow].SetMother(param0[irow]);
5659 kinks[irow].SetDaughter(param1[irow]);
5660 kinks[irow].Update();
5663 // choose kink with biggest change of angle
5665 Double_t maxchange= 0;
5666 for (Int_t irow=1;irow<19;irow++){
5667 if (TMath::Abs(kinks[irow].GetR())>240.) continue;
5668 if (TMath::Abs(kinks[irow].GetR())<110.) continue;
5669 Float_t quality = TMath::Abs(kinks[irow].GetAngle(2))/(3.+TMath::Abs(kinks[irow].GetR()-param0[irow].GetX()));
5670 if ( quality > maxchange){
5671 maxchange = quality;
5678 if (index<0) return 0;
5680 Int_t row0 = GetRowNumber(kinks[index].GetR()); //row 0 estimate
5681 seed0 = new AliTPCseed(param0[index]);
5682 seed1 = new AliTPCseed(param1[index]);
5683 seed0->Reset(kFALSE);
5684 seed1->Reset(kFALSE);
5685 seed0->ResetCovariance(10.);
5686 seed1->ResetCovariance(10.);
5687 FollowProlongation(*seed0,0);
5688 FollowBackProlongation(*seed1,158);
5689 mother = *seed0; // backup mother at position 0
5690 seed0->Reset(kFALSE);
5691 seed1->Reset(kFALSE);
5692 seed0->ResetCovariance(10.);
5693 seed1->ResetCovariance(10.);
5695 first = TMath::Max(row0-20,0);
5696 last = TMath::Min(row0+20,158);
5698 for (Int_t irow=0; irow<20;irow++){
5699 rows[irow] = first +((last-first)*irow)/19;
5701 // store parameters along the track
5703 for (Int_t irow=0;irow<20;irow++){
5704 FollowBackProlongation(*seed0, rows[irow]);
5705 FollowProlongation(*seed1,rows[19-irow]);
5706 param0[irow] = *seed0;
5707 param1[19-irow] = *seed1;
5711 for (Int_t irow=0; irow<19;irow++){
5712 kinks[irow].SetMother(param0[irow]);
5713 kinks[irow].SetDaughter(param1[irow]);
5714 // param0[irow].Dump();
5715 //param1[irow].Dump();
5716 kinks[irow].Update();
5719 // choose kink with biggest change of angle
5722 for (Int_t irow=0;irow<20;irow++){
5723 if (TMath::Abs(kinks[irow].GetR())>250.) continue;
5724 if (TMath::Abs(kinks[irow].GetR())<90.) continue;
5725 Float_t quality = TMath::Abs(kinks[irow].GetAngle(2))/(3.+TMath::Abs(kinks[irow].GetR()-param0[irow].GetX()));
5726 if ( quality > maxchange){
5727 maxchange = quality;
5734 if (index==-1 || param0[index].GetNumberOfClusters()+param1[index].GetNumberOfClusters()<100){
5740 // Float_t anglesigma = TMath::Sqrt(param0[index].fC22+param0[index].fC33+param1[index].fC22+param1[index].fC33);
5742 kink.SetMother(param0[index]);
5743 kink.SetDaughter(param1[index]);
5746 Double_t chi2P2 = param0[index].GetParameter()[2]-param1[index].GetParameter()[2];
5748 chi2P2/=param0[index].GetCovariance()[5]+param1[index].GetCovariance()[5];
5749 Double_t chi2P3 = param0[index].GetParameter()[3]-param1[index].GetParameter()[3];
5751 chi2P3/=param0[index].GetCovariance()[9]+param1[index].GetCovariance()[9];
5753 if (AliTPCReconstructor::StreamLevel()>1) {
5754 (*fDebugStreamer)<<"kinkHpt"<<
5757 "p0.="<<¶m0[index]<<
5758 "p1.="<<¶m1[index]<<
5762 if ( chi2P2+chi2P3<AliTPCReconstructor::GetRecoParam()->GetKinkAngleCutChi2(0)){
5769 row0 = GetRowNumber(kink.GetR());
5770 kink.SetTPCRow0(row0);
5771 kink.SetLabel(CookLabel(seed0,0.5,0,row0),0);
5772 kink.SetLabel(CookLabel(seed1,0.5,row0,158),1);
5773 kink.SetIndex(-10,0);
5774 kink.SetIndex(int(param0[index].GetNumberOfClusters()+param1[index].GetNumberOfClusters()),1);
5775 kink.SetTPCncls(param0[index].GetNumberOfClusters(),0);
5776 kink.SetTPCncls(param1[index].GetNumberOfClusters(),1);
5779 // new (&mother) AliTPCseed(param0[index]);
5780 daughter = param1[index];
5781 daughter.SetLabel(kink.GetLabel(1));
5782 param0[index].Reset(kTRUE);
5783 FollowProlongation(param0[index],0);
5784 mother = param0[index];
5785 mother.SetLabel(kink.GetLabel(0));
5786 if ( chi2P2+chi2P3<AliTPCReconstructor::GetRecoParam()->GetKinkAngleCutChi2(1)){
5798 AliTPCseed* AliTPCtrackerMI::ReSeed(AliTPCseed *t)
5801 // reseed - refit - track
5804 // Int_t last = fSectors->GetNRows()-1;
5806 if (fSectors == fOuterSec){
5807 first = TMath::Max(first, t->GetFirstPoint()-fInnerSec->GetNRows());
5811 first = t->GetFirstPoint();
5813 AliTPCseed * seed = MakeSeed(t,0.1,0.5,0.9);
5814 FollowBackProlongation(*t,fSectors->GetNRows()-1);
5816 FollowProlongation(*t,first);
5826 //_____________________________________________________________________________
5827 Int_t AliTPCtrackerMI::ReadSeeds(const TFile *inp) {
5828 //-----------------------------------------------------------------
5829 // This function reades track seeds.
5830 //-----------------------------------------------------------------
5831 TDirectory *savedir=gDirectory;
5833 TFile *in=(TFile*)inp;
5834 if (!in->IsOpen()) {
5835 cerr<<"AliTPCtrackerMI::ReadSeeds(): input file is not open !\n";
5840 TTree *seedTree=(TTree*)in->Get("Seeds");
5842 cerr<<"AliTPCtrackerMI::ReadSeeds(): ";
5843 cerr<<"can't get a tree with track seeds !\n";
5846 AliTPCtrack *seed=new AliTPCtrack;
5847 seedTree->SetBranchAddress("tracks",&seed);
5849 if (fSeeds==0) fSeeds=new TObjArray(15000);
5851 Int_t n=(Int_t)seedTree->GetEntries();
5852 for (Int_t i=0; i<n; i++) {
5853 seedTree->GetEvent(i);
5854 fSeeds->AddLast(new AliTPCseed(*seed/*,seed->GetAlpha()*/));
5863 Int_t AliTPCtrackerMI::Clusters2Tracks (AliESDEvent *const esd)
5867 if (fSeeds) DeleteSeeds();
5870 AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ;
5871 transform->SetCurrentTimeStamp( esd->GetTimeStamp());
5872 transform->SetCurrentRun(esd->GetRunNumber());
5875 if (!fSeeds) return 1;
5877 if (AliTPCReconstructor::StreamLevel()>3) DumpClusters(0,fSeeds);
5883 //_____________________________________________________________________________
5884 Int_t AliTPCtrackerMI::Clusters2Tracks() {
5885 //-----------------------------------------------------------------
5886 // This is a track finder.
5887 //-----------------------------------------------------------------
5888 TDirectory *savedir=gDirectory;
5892 fSeeds = Tracking();
5895 Info("Clusters2Tracks","Time for tracking: \t");timer.Print();timer.Start();
5897 //activate again some tracks
5898 for (Int_t i=0; i<fSeeds->GetEntriesFast(); i++) {
5899 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
5901 Int_t nc=t.GetNumberOfClusters();
5903 delete fSeeds->RemoveAt(i);
5907 if (pt->GetRemoval()==10) {
5908 if (pt->GetDensityFirst(20)>0.8 || pt->GetDensityFirst(30)>0.8 || pt->GetDensityFirst(40)>0.7)
5909 pt->Desactivate(10); // make track again active
5911 pt->Desactivate(20);
5912 delete fSeeds->RemoveAt(i);
5917 RemoveUsed2(fSeeds,0.85,0.85,0);
5918 if (AliTPCReconstructor::GetRecoParam()->GetDoKinks()) FindKinks(fSeeds,fEvent);
5919 //FindCurling(fSeeds, fEvent,0);
5920 if (AliTPCReconstructor::StreamLevel()>5) FindMultiMC(fSeeds, fEvent,-1); // find multi found tracks
5921 RemoveUsed2(fSeeds,0.5,0.4,20);
5922 FindSplitted(fSeeds, fEvent,0); // find multi found tracks
5923 if (AliTPCReconstructor::StreamLevel()>5) FindMultiMC(fSeeds, fEvent,0); // find multi found tracks
5926 // // refit short tracks
5928 Int_t nseed=fSeeds->GetEntriesFast();
5931 for (Int_t i=0; i<nseed; i++) {
5932 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
5934 Int_t nc=t.GetNumberOfClusters();
5936 delete fSeeds->RemoveAt(i);
5939 CookLabel(pt,0.1); //For comparison only
5940 //if ((pt->IsActive() || (pt->fRemoval==10) )&& nc>50 &&pt->GetNumberOfClusters()>0.4*pt->fNFoundable){
5941 if ((pt->IsActive() || (pt->GetRemoval()==10) )){
5943 if (fDebug>0) cerr<<found<<'\r';
5947 delete fSeeds->RemoveAt(i);
5951 //RemoveOverlap(fSeeds,0.99,7,kTRUE);
5953 //RemoveUsed(fSeeds,0.9,0.9,6);
5955 nseed=fSeeds->GetEntriesFast();
5957 for (Int_t i=0; i<nseed; i++) {
5958 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
5960 Int_t nc=t.GetNumberOfClusters();
5962 delete fSeeds->RemoveAt(i);
5966 t.CookdEdx(0.02,0.6);
5967 // CheckKinkPoint(&t,0.05);
5968 //if ((pt->IsActive() || (pt->fRemoval==10) )&& nc>50 &&pt->GetNumberOfClusters()>0.4*pt->fNFoundable){
5969 if ((pt->IsActive() || (pt->GetRemoval()==10) )){
5977 delete fSeeds->RemoveAt(i);
5978 //AliTPCseed * seed1 = ReSeed(pt,0.05,0.5,1);
5980 // FollowProlongation(*seed1,0);
5981 // Int_t n = seed1->GetNumberOfClusters();
5982 // printf("fP4\t%f\t%f\n",seed1->GetC(),pt->GetC());
5983 // printf("fN\t%d\t%d\n", seed1->GetNumberOfClusters(),pt->GetNumberOfClusters());
5986 //AliTPCseed * seed2 = ReSeed(pt,0.95,0.5,0.05);
5990 SortTracks(fSeeds, 1);
5994 PrepareForBackProlongation(fSeeds,5.);
5995 PropagateBack(fSeeds);
5996 printf("Time for back propagation: \t");timer.Print();timer.Start();
6000 PrepareForProlongation(fSeeds,5.);
6001 PropagateForward2(fSeeds);
6003 printf("Time for FORWARD propagation: \t");timer.Print();timer.Start();
6004 // RemoveUsed(fSeeds,0.7,0.7,6);
6005 //RemoveOverlap(fSeeds,0.9,7,kTRUE);
6007 nseed=fSeeds->GetEntriesFast();
6009 for (Int_t i=0; i<nseed; i++) {
6010 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
6012 Int_t nc=t.GetNumberOfClusters();
6014 delete fSeeds->RemoveAt(i);
6017 t.CookdEdx(0.02,0.6);
6018 // CookLabel(pt,0.1); //For comparison only
6019 //if ((pt->IsActive() || (pt->fRemoval==10) )&& nc>50 &&pt->GetNumberOfClusters()>0.4*pt->fNFoundable){
6020 if ((pt->IsActive() || (pt->fRemoval==10) )){
6021 cerr<<found++<<'\r';
6024 delete fSeeds->RemoveAt(i);
6029 // fNTracks = found;
6031 Info("Clusters2Tracks","Time for overlap removal, track writing and dedx cooking: \t"); timer.Print();timer.Start();
6034 // cerr<<"Number of found tracks : "<<"\t"<<found<<endl;
6035 Info("Clusters2Tracks","Number of found tracks %d",found);
6037 // UnloadClusters();
6042 void AliTPCtrackerMI::Tracking(TObjArray * arr)
6045 // tracking of the seeds
6048 fSectors = fOuterSec;
6049 ParallelTracking(arr,150,63);
6050 fSectors = fOuterSec;
6051 ParallelTracking(arr,63,0);
6054 TObjArray * AliTPCtrackerMI::Tracking(Int_t seedtype, Int_t i1, Int_t i2, Float_t cuts[4], Float_t dy, Int_t dsec)
6059 TObjArray * arr = new TObjArray;
6061 fSectors = fOuterSec;
6064 for (Int_t sec=0;sec<fkNOS;sec++){
6065 if (seedtype==3) MakeSeeds3(arr,sec,i1,i2,cuts,dy, dsec);
6066 if (seedtype==4) MakeSeeds5(arr,sec,i1,i2,cuts,dy);
6067 if (seedtype==2) MakeSeeds2(arr,sec,i1,i2,cuts,dy);
6070 Info("Tracking","\nSeeding - %d\t%d\t%d\t%d\n",seedtype,i1,i2,arr->GetEntriesFast());
6082 TObjArray * AliTPCtrackerMI::Tracking()
6086 if (AliTPCReconstructor::GetRecoParam()->GetSpecialSeeding()) return TrackingSpecial();
6089 Int_t nup=fOuterSec->GetNRows()+fInnerSec->GetNRows();
6091 TObjArray * seeds = new TObjArray;
6093 Int_t fLastSeedRowSec=AliTPCReconstructor::GetRecoParam()->GetLastSeedRowSec();
6094 Int_t gapPrim = AliTPCReconstructor::GetRecoParam()->GetSeedGapPrim();
6095 Int_t gapSec = AliTPCReconstructor::GetRecoParam()->GetSeedGapSec();
6103 Float_t fnumber = 3.0;
6104 Float_t fdensity = 3.0;
6109 for (Int_t delta = 0; delta<18; delta+=gapPrim){
6113 arr = Tracking(3,nup-1-delta,nup-1-delta-gap,cuts,-1,1);
6114 SumTracks(seeds,arr);
6115 SignClusters(seeds,fnumber,fdensity);
6117 for (Int_t i=2;i<6;i+=2){
6118 // seed high pt tracks
6121 arr = Tracking(3,nup-i-delta,nup-i-delta-gap,cuts,-1,0);
6122 SumTracks(seeds,arr);
6123 SignClusters(seeds,fnumber,fdensity);
6128 // RemoveUsed(seeds,0.9,0.9,1);
6129 // UnsignClusters();
6130 // SignClusters(seeds,fnumber,fdensity);
6134 for (Int_t delta = 20; delta<120; delta+=gapPrim){
6136 // seed high pt tracks
6140 arr = Tracking(3,nup-delta,nup-delta-gap,cuts,-1);
6141 SumTracks(seeds,arr);
6142 SignClusters(seeds,fnumber,fdensity);
6147 arr = Tracking(3,nup-delta-5,nup-delta-5-gap,cuts,-1);
6148 SumTracks(seeds,arr);
6149 SignClusters(seeds,fnumber,fdensity);
6160 Info("Tracking()","\n\nPrimary seeding\t%d\n\n",seeds->GetEntriesFast());
6164 // RemoveUsed(seeds,0.75,0.75,1);
6166 //SignClusters(seeds,fnumber,fdensity);
6175 arr = Tracking(4,nup-1,nup-1-gap,cuts,-1);
6176 SumTracks(seeds,arr);
6177 SignClusters(seeds,fnumber,fdensity);
6179 arr = Tracking(4,nup-2,nup-2-gap,cuts,-1);
6180 SumTracks(seeds,arr);
6181 SignClusters(seeds,fnumber,fdensity);
6183 arr = Tracking(4,nup-3,nup-3-gap,cuts,-1);
6184 SumTracks(seeds,arr);
6185 SignClusters(seeds,fnumber,fdensity);
6187 arr = Tracking(4,nup-5,nup-5-gap,cuts,-1);
6188 SumTracks(seeds,arr);
6189 SignClusters(seeds,fnumber,fdensity);
6191 arr = Tracking(4,nup-7,nup-7-gap,cuts,-1);
6192 SumTracks(seeds,arr);
6193 SignClusters(seeds,fnumber,fdensity);
6196 arr = Tracking(4,nup-9,nup-9-gap,cuts,-1);
6197 SumTracks(seeds,arr);
6198 SignClusters(seeds,fnumber,fdensity);
6202 for (Int_t delta = 9; delta<30; delta+=gapSec){
6208 arr = Tracking(4,nup-1-delta,nup-1-delta-gap,cuts,-1);
6209 SumTracks(seeds,arr);
6210 SignClusters(seeds,fnumber,fdensity);
6212 arr = Tracking(4,nup-3-delta,nup-5-delta-gap,cuts,4);
6213 SumTracks(seeds,arr);
6214 SignClusters(seeds,fnumber,fdensity);
6227 for (Int_t delta = 30; delta<fLastSeedRowSec; delta+=gapSec){
6233 arr = Tracking(4,nup-1-delta,nup-1-delta-gap,cuts,-1);
6234 SumTracks(seeds,arr);
6235 SignClusters(seeds,fnumber,fdensity);
6237 arr = Tracking(4,nup-5-delta,nup-5-delta-gap,cuts,5 );
6238 SumTracks(seeds,arr);
6239 SignClusters(seeds,fnumber,fdensity);
6243 Info("Tracking()","\n\nSecondary seeding\t%d\n\n",seeds->GetEntriesFast());
6254 TObjArray * AliTPCtrackerMI::TrackingSpecial()
6257 // seeding adjusted for laser and cosmic tests - short tracks with big inclination angle
6258 // no primary vertex seeding tried
6262 Int_t nup=fOuterSec->GetNRows()+fInnerSec->GetNRows();
6264 TObjArray * seeds = new TObjArray;
6269 Float_t fnumber = 3.0;
6270 Float_t fdensity = 3.0;
6273 cuts[0] = AliTPCReconstructor::GetRecoParam()->GetMaxC(); // max curvature
6274 cuts[1] = 3.5; // max tan(phi) angle for seeding
6275 cuts[2] = 3.; // not used (cut on z primary vertex)
6276 cuts[3] = 3.5; // max tan(theta) angle for seeding
6278 for (Int_t delta = 0; nup-delta-gap-1>0; delta+=3){
6280 arr = Tracking(4,nup-1-delta,nup-1-delta-gap,cuts,-1);
6281 SumTracks(seeds,arr);
6282 SignClusters(seeds,fnumber,fdensity);
6286 Info("Tracking()","\n\nSecondary seeding\t%d\n\n",seeds->GetEntriesFast());
6297 void AliTPCtrackerMI::SumTracks(TObjArray *arr1,TObjArray *&arr2) const
6300 //sum tracks to common container
6301 //remove suspicious tracks
6302 Int_t nseed = arr2->GetEntriesFast();
6303 for (Int_t i=0;i<nseed;i++){
6304 AliTPCseed *pt=(AliTPCseed*)arr2->UncheckedAt(i);
6307 // remove tracks with too big curvature
6309 if (TMath::Abs(pt->GetC())>AliTPCReconstructor::GetRecoParam()->GetMaxC()){
6310 delete arr2->RemoveAt(i);
6313 // REMOVE VERY SHORT TRACKS
6314 if (pt->GetNumberOfClusters()<20){
6315 delete arr2->RemoveAt(i);
6318 // NORMAL ACTIVE TRACK
6319 if (pt->IsActive()){
6320 arr1->AddLast(arr2->RemoveAt(i));
6323 //remove not usable tracks
6324 if (pt->GetRemoval()!=10){
6325 delete arr2->RemoveAt(i);
6329 // ENABLE ONLY ENOUGH GOOD STOPPED TRACKS
6330 if (pt->GetDensityFirst(20)>0.8 || pt->GetDensityFirst(30)>0.8 || pt->GetDensityFirst(40)>0.7)
6331 arr1->AddLast(arr2->RemoveAt(i));
6333 delete arr2->RemoveAt(i);
6337 delete arr2; arr2 = 0;
6342 void AliTPCtrackerMI::ParallelTracking(TObjArray *const arr, Int_t rfirst, Int_t rlast)
6345 // try to track in parralel
6347 Int_t nseed=arr->GetEntriesFast();
6348 //prepare seeds for tracking
6349 for (Int_t i=0; i<nseed; i++) {
6350 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i), &t=*pt;
6352 if (!t.IsActive()) continue;
6353 // follow prolongation to the first layer
6354 if ( (fSectors ==fInnerSec) || (t.GetFirstPoint()-fkParam->GetNRowLow()>rfirst+1) )
6355 FollowProlongation(t, rfirst+1);
6360 for (Int_t nr=rfirst; nr>=rlast; nr--){
6361 if (nr<fInnerSec->GetNRows())
6362 fSectors = fInnerSec;
6364 fSectors = fOuterSec;
6365 // make indexes with the cluster tracks for given
6367 // find nearest cluster
6368 for (Int_t i=0; i<nseed; i++) {
6369 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i), &t=*pt;
6371 if (nr==80) pt->UpdateReference();
6372 if (!pt->IsActive()) continue;
6373 // if ( (fSectors ==fOuterSec) && (pt->fFirstPoint-fkParam->GetNRowLow())<nr) continue;
6374 if (pt->GetRelativeSector()>17) {
6377 UpdateClusters(t,nr);
6379 // prolonagate to the nearest cluster - if founded
6380 for (Int_t i=0; i<nseed; i++) {
6381 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
6383 if (!pt->IsActive()) continue;
6384 // if ((fSectors ==fOuterSec) && (pt->fFirstPoint-fkParam->GetNRowLow())<nr) continue;
6385 if (pt->GetRelativeSector()>17) {
6388 FollowToNextCluster(*pt,nr);
6393 void AliTPCtrackerMI::PrepareForBackProlongation(TObjArray *const arr,Float_t fac) const
6397 // if we use TPC track itself we have to "update" covariance
6399 Int_t nseed= arr->GetEntriesFast();
6400 for (Int_t i=0;i<nseed;i++){
6401 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
6405 //rotate to current local system at first accepted point
6406 Int_t index = pt->GetClusterIndex2(pt->GetFirstPoint());
6407 Int_t sec = (index&0xff000000)>>24;
6409 Float_t angle1 = fInnerSec->GetAlpha()*sec+fInnerSec->GetAlphaShift();
6410 if (angle1>TMath::Pi())
6411 angle1-=2.*TMath::Pi();
6412 Float_t angle2 = pt->GetAlpha();
6414 if (TMath::Abs(angle1-angle2)>0.001){
6415 if (!pt->Rotate(angle1-angle2)) return;
6416 //angle2 = pt->GetAlpha();
6417 //pt->fRelativeSector = pt->GetAlpha()/fInnerSec->GetAlpha();
6418 //if (pt->GetAlpha()<0)
6419 // pt->fRelativeSector+=18;
6420 //sec = pt->fRelativeSector;
6429 void AliTPCtrackerMI::PrepareForProlongation(TObjArray *const arr, Float_t fac) const
6433 // if we use TPC track itself we have to "update" covariance
6435 Int_t nseed= arr->GetEntriesFast();
6436 for (Int_t i=0;i<nseed;i++){
6437 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
6440 pt->SetFirstPoint(pt->GetLastPoint());
6448 Int_t AliTPCtrackerMI::PropagateBack(TObjArray *const arr)
6451 // make back propagation
6453 Int_t nseed= arr->GetEntriesFast();
6454 for (Int_t i=0;i<nseed;i++){
6455 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
6456 if (pt&& pt->GetKinkIndex(0)<=0) {
6457 //AliTPCseed *pt2 = new AliTPCseed(*pt);
6458 fSectors = fInnerSec;
6459 //FollowBackProlongation(*pt,fInnerSec->GetNRows()-1);
6460 //fSectors = fOuterSec;
6461 FollowBackProlongation(*pt,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1);
6462 //if (pt->GetNumberOfClusters()<(pt->fEsd->GetTPCclusters(0)) ){
6463 // Error("PropagateBack","Not prolonged track %d",pt->GetLabel());
6464 // FollowBackProlongation(*pt2,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1);
6467 if (pt&& pt->GetKinkIndex(0)>0) {
6468 AliESDkink * kink = fEvent->GetKink(pt->GetKinkIndex(0)-1);
6469 pt->SetFirstPoint(kink->GetTPCRow0());
6470 fSectors = fInnerSec;
6471 FollowBackProlongation(*pt,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1);
6479 Int_t AliTPCtrackerMI::PropagateForward2(TObjArray *const arr)
6482 // make forward propagation
6484 Int_t nseed= arr->GetEntriesFast();
6486 for (Int_t i=0;i<nseed;i++){
6487 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
6489 FollowProlongation(*pt,0);
6498 Int_t AliTPCtrackerMI::PropagateForward()
6501 // propagate track forward
6503 Int_t nseed = fSeeds->GetEntriesFast();
6504 for (Int_t i=0;i<nseed;i++){
6505 AliTPCseed *pt = (AliTPCseed*)fSeeds->UncheckedAt(i);
6507 AliTPCseed &t = *pt;
6508 Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
6509 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
6510 if (alpha < 0. ) alpha += 2.*TMath::Pi();
6511 t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN);
6515 fSectors = fOuterSec;
6516 ParallelTracking(fSeeds,fOuterSec->GetNRows()+fInnerSec->GetNRows()-1,fInnerSec->GetNRows());
6517 fSectors = fInnerSec;
6518 ParallelTracking(fSeeds,fInnerSec->GetNRows()-1,0);
6527 Int_t AliTPCtrackerMI::PropagateBack(AliTPCseed *const pt, Int_t row0, Int_t row1)
6530 // make back propagation, in between row0 and row1
6534 fSectors = fInnerSec;
6537 if (row1<fSectors->GetNRows())
6540 r1 = fSectors->GetNRows()-1;
6542 if (row0<fSectors->GetNRows()&& r1>0 )
6543 FollowBackProlongation(*pt,r1);
6544 if (row1<=fSectors->GetNRows())
6547 r1 = row1 - fSectors->GetNRows();
6548 if (r1<=0) return 0;
6549 if (r1>=fOuterSec->GetNRows()) return 0;
6550 fSectors = fOuterSec;
6551 return FollowBackProlongation(*pt,r1);
6559 void AliTPCtrackerMI::GetShape(AliTPCseed * seed, Int_t row)
6563 AliTPCClusterParam * clparam = AliTPCcalibDB::Instance()->GetClusterParam();
6564 Float_t zdrift = TMath::Abs((fkParam->GetZLength(0)-TMath::Abs(seed->GetZ())));
6565 Int_t type = (seed->GetSector() < fkParam->GetNSector()/2) ? 0: (row>126) ? 1:2;
6566 Double_t angulary = seed->GetSnp();
6568 if (TMath::Abs(angulary)>AliTPCReconstructor::GetMaxSnpTracker()) {
6569 angulary = TMath::Sign(AliTPCReconstructor::GetMaxSnpTracker(),angulary);
6572 angulary = angulary*angulary/((1.-angulary)*(1.+angulary));
6573 Double_t angularz = seed->GetTgl()*seed->GetTgl()*(1.+angulary);
6575 Double_t sigmay = clparam->GetRMS0(0,type,zdrift,TMath::Sqrt(TMath::Abs(angulary)));
6576 Double_t sigmaz = clparam->GetRMS0(1,type,zdrift,TMath::Sqrt(TMath::Abs(angularz)));
6577 seed->SetCurrentSigmaY2(sigmay*sigmay);
6578 seed->SetCurrentSigmaZ2(sigmaz*sigmaz);
6579 // Float_t sd2 = TMath::Abs((fkParam->GetZLength(0)-TMath::Abs(seed->GetZ())))*fkParam->GetDiffL()*fkParam->GetDiffL();
6580 // // Float_t padlength = fkParam->GetPadPitchLength(seed->fSector);
6581 // Float_t padlength = GetPadPitchLength(row);
6583 // Float_t sresy = (seed->GetSector() < fkParam->GetNSector()/2) ? 0.2 :0.3;
6584 // seed->SetCurrentSigmaY2(sd2+padlength*padlength*angulary/12.+sresy*sresy);
6586 // Float_t sresz = fkParam->GetZSigma();
6587 // seed->SetCurrentSigmaZ2(sd2+padlength*padlength*angularz*angularz*(1+angulary)/12.+sresz*sresz);
6589 Float_t wy = GetSigmaY(seed);
6590 Float_t wz = GetSigmaZ(seed);
6593 if (TMath::Abs(wy/seed->fCurrentSigmaY2-1)>0.0001 || TMath::Abs(wz/seed->fCurrentSigmaZ2-1)>0.0001 ){
6594 printf("problem\n");
6601 //__________________________________________________________________________
6602 void AliTPCtrackerMI::CookLabel(AliKalmanTrack *tk, Float_t wrong) const {
6603 //--------------------------------------------------------------------
6604 //This function "cooks" a track label. If label<0, this track is fake.
6605 //--------------------------------------------------------------------
6606 AliTPCseed * t = dynamic_cast<AliTPCseed*>(tk);
6608 printf("%s:%d wrong type \n",(char*)__FILE__,__LINE__);
6612 Int_t noc=t->GetNumberOfClusters();
6614 //printf("\nnot founded prolongation\n\n\n");
6620 AliTPCclusterMI *clusters[160];
6622 for (Int_t i=0;i<160;i++) {
6629 for (i=0; i<160 && current<noc; i++) {
6631 Int_t index=t->GetClusterIndex2(i);
6632 if (index<=0) continue;
6633 if (index&0x8000) continue;
6635 //clusters[current]=GetClusterMI(index);
6636 if (t->GetClusterPointer(i)){
6637 clusters[current]=t->GetClusterPointer(i);
6643 Int_t lab=123456789;
6644 for (i=0; i<noc; i++) {
6645 AliTPCclusterMI *c=clusters[i];
6647 lab=TMath::Abs(c->GetLabel(0));
6649 for (j=0; j<noc; j++) if (lb[j]==lab || mx[j]==0) break;
6655 for (i=0; i<noc; i++) if (mx[i]>max) {max=mx[i]; lab=lb[i];}
6657 for (i=0; i<noc; i++) {
6658 AliTPCclusterMI *c=clusters[i];
6660 if (TMath::Abs(c->GetLabel(1)) == lab ||
6661 TMath::Abs(c->GetLabel(2)) == lab ) max++;
6663 if (noc<=0) { lab=-1; return;}
6664 if ((1.- Float_t(max)/(noc)) > wrong) lab=-lab;
6667 Int_t tail=Int_t(0.10*noc);
6670 for (i=1; i<160&&ind<tail; i++) {
6671 // AliTPCclusterMI *c=clusters[noc-i];
6672 AliTPCclusterMI *c=clusters[i];
6674 if (lab == TMath::Abs(c->GetLabel(0)) ||
6675 lab == TMath::Abs(c->GetLabel(1)) ||
6676 lab == TMath::Abs(c->GetLabel(2))) max++;
6679 if (max < Int_t(0.5*tail)) lab=-lab;
6686 //delete[] clusters;
6690 //__________________________________________________________________________
6691 Int_t AliTPCtrackerMI::CookLabel(AliTPCseed *const t, Float_t wrong,Int_t first, Int_t last) const {
6692 //--------------------------------------------------------------------
6693 //This function "cooks" a track label. If label<0, this track is fake.
6694 //--------------------------------------------------------------------
6695 Int_t noc=t->GetNumberOfClusters();
6697 //printf("\nnot founded prolongation\n\n\n");
6703 AliTPCclusterMI *clusters[160];
6705 for (Int_t i=0;i<160;i++) {
6712 for (i=0; i<160 && current<noc; i++) {
6713 if (i<first) continue;
6714 if (i>last) continue;
6715 Int_t index=t->GetClusterIndex2(i);
6716 if (index<=0) continue;
6717 if (index&0x8000) continue;
6719 //clusters[current]=GetClusterMI(index);
6720 if (t->GetClusterPointer(i)){
6721 clusters[current]=t->GetClusterPointer(i);
6726 //if (noc<5) return -1;
6727 Int_t lab=123456789;
6728 for (i=0; i<noc; i++) {
6729 AliTPCclusterMI *c=clusters[i];
6731 lab=TMath::Abs(c->GetLabel(0));
6733 for (j=0; j<noc; j++) if (lb[j]==lab || mx[j]==0) break;
6739 for (i=0; i<noc; i++) if (mx[i]>max) {max=mx[i]; lab=lb[i];}
6741 for (i=0; i<noc; i++) {
6742 AliTPCclusterMI *c=clusters[i];
6744 if (TMath::Abs(c->GetLabel(1)) == lab ||
6745 TMath::Abs(c->GetLabel(2)) == lab ) max++;
6747 if (noc<=0) { lab=-1; return -1;}
6748 if ((1.- Float_t(max)/(noc)) > wrong) lab=-lab;
6751 Int_t tail=Int_t(0.10*noc);
6754 for (i=1; i<160&&ind<tail; i++) {
6755 // AliTPCclusterMI *c=clusters[noc-i];
6756 AliTPCclusterMI *c=clusters[i];
6758 if (lab == TMath::Abs(c->GetLabel(0)) ||
6759 lab == TMath::Abs(c->GetLabel(1)) ||
6760 lab == TMath::Abs(c->GetLabel(2))) max++;
6763 if (max < Int_t(0.5*tail)) lab=-lab;
6766 // t->SetLabel(lab);
6770 //delete[] clusters;
6774 Int_t AliTPCtrackerMI::GetRowNumber(Double_t x[3]) const
6776 //return pad row number for given x vector
6777 Float_t phi = TMath::ATan2(x[1],x[0]);
6778 if(phi<0) phi=2.*TMath::Pi()+phi;
6779 // Get the local angle in the sector philoc
6780 const Float_t kRaddeg = 180/3.14159265358979312;
6781 Float_t phiangle = (Int_t (phi*kRaddeg/20.) + 0.5)*20./kRaddeg;
6782 Double_t localx = x[0]*TMath::Cos(phiangle)-x[1]*TMath::Sin(phiangle);
6783 return GetRowNumber(localx);
6788 void AliTPCtrackerMI::MakeBitmaps(AliTPCseed *t)
6790 //-----------------------------------------------------------------------
6791 // Fill the cluster and sharing bitmaps of the track
6792 //-----------------------------------------------------------------------
6794 Int_t firstpoint = 0;
6795 Int_t lastpoint = 159;
6796 AliTPCTrackerPoint *point;
6797 AliTPCclusterMI *cluster;
6799 for (int iter=firstpoint; iter<lastpoint; iter++) {
6800 // Change to cluster pointers to see if we have a cluster at given padrow
6801 cluster = t->GetClusterPointer(iter);
6803 t->SetClusterMapBit(iter, kTRUE);
6804 point = t->GetTrackPoint(iter);
6805 if (point->IsShared())
6806 t->SetSharedMapBit(iter,kTRUE);
6808 t->SetSharedMapBit(iter, kFALSE);
6811 t->SetClusterMapBit(iter, kFALSE);
6812 t->SetSharedMapBit(iter, kFALSE);
6817 Bool_t AliTPCtrackerMI::IsFindable(AliTPCseed & track){
6819 // return flag if there is findable cluster at given position
6822 Float_t z = track.GetZ();
6824 if (TMath::Abs(z)<(AliTPCReconstructor::GetCtgRange()*track.GetX()+kDeltaZ) &&
6825 TMath::Abs(z)<fkParam->GetZLength(0) &&
6826 (TMath::Abs(track.GetSnp())<AliTPCReconstructor::GetMaxSnpTracker()))
6832 void AliTPCtrackerMI::AddCovariance(AliTPCseed * seed){
6834 // Adding systematic error estimate to the covariance matrix
6835 // !!!! the systematic error for element 4 is in 1/cm not in pt
6837 const Double_t *param = AliTPCReconstructor::GetRecoParam()->GetSystematicError();
6839 // use only the diagonal part if not specified otherwise
6840 if (!AliTPCReconstructor::GetRecoParam()->GetUseSystematicCorrelation()) return AddCovarianceAdd(seed);
6842 Double_t *covarS= (Double_t*)seed->GetCovariance();
6843 Double_t factor[5]={1,1,1,1,1};
6844 Double_t facC = AliTracker::GetBz()*kB2C;
6845 factor[0]= TMath::Sqrt(TMath::Abs((covarS[0] + param[0]*param[0])/covarS[0]));
6846 factor[1]= TMath::Sqrt(TMath::Abs((covarS[2] + param[1]*param[1])/covarS[2]));
6847 factor[2]= TMath::Sqrt(TMath::Abs((covarS[5] + param[2]*param[2])/covarS[5]));
6848 factor[3]= TMath::Sqrt(TMath::Abs((covarS[9] + param[3]*param[3])/covarS[9]));
6849 factor[4]= TMath::Sqrt(TMath::Abs((covarS[14] + facC*facC*param[4]*param[4])/covarS[14]));
6855 for (Int_t i=0; i<5; i++){
6856 for (Int_t j=i; j<5; j++){
6857 Int_t index=seed->GetIndex(i,j);
6858 covarS[index]*=factor[i]*factor[j];
6864 void AliTPCtrackerMI::AddCovarianceAdd(AliTPCseed * seed){
6866 // Adding systematic error - as additive factor without correlation
6868 const Double_t *param = AliTPCReconstructor::GetRecoParam()->GetSystematicError();
6869 Double_t *covarIn= (Double_t*)seed->GetCovariance();
6871 for (Int_t i=0;i<15;i++) covar[i]=0;
6877 covar[0] = param[0]*param[0];
6878 covar[2] = param[1]*param[1];
6879 covar[5] = param[2]*param[2];
6880 covar[9] = param[3]*param[3];
6881 Double_t facC = AliTracker::GetBz()*kB2C;
6882 covar[14]= param[4]*param[4]*facC*facC;
6884 covar[1]=TMath::Sqrt((covar[0]*covar[2]))*covarIn[1]/TMath::Sqrt((covarIn[0]*covarIn[2]));
6886 covar[3]=TMath::Sqrt((covar[0]*covar[5]))*covarIn[3]/TMath::Sqrt((covarIn[0]*covarIn[5]));
6887 covar[4]=TMath::Sqrt((covar[2]*covar[5]))*covarIn[4]/TMath::Sqrt((covarIn[2]*covarIn[5]));
6889 covar[6]=TMath::Sqrt((covar[0]*covar[9]))*covarIn[6]/TMath::Sqrt((covarIn[0]*covarIn[9]));
6890 covar[7]=TMath::Sqrt((covar[2]*covar[9]))*covarIn[7]/TMath::Sqrt((covarIn[2]*covarIn[9]));
6891 covar[8]=TMath::Sqrt((covar[5]*covar[9]))*covarIn[8]/TMath::Sqrt((covarIn[5]*covarIn[9]));
6893 covar[10]=TMath::Sqrt((covar[0]*covar[14]))*covarIn[10]/TMath::Sqrt((covarIn[0]*covarIn[14]));
6894 covar[11]=TMath::Sqrt((covar[2]*covar[14]))*covarIn[11]/TMath::Sqrt((covarIn[2]*covarIn[14]));
6895 covar[12]=TMath::Sqrt((covar[5]*covar[14]))*covarIn[12]/TMath::Sqrt((covarIn[5]*covarIn[14]));
6896 covar[13]=TMath::Sqrt((covar[9]*covar[14]))*covarIn[13]/TMath::Sqrt((covarIn[9]*covarIn[14]));
6898 seed->AddCovariance(covar);