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);
374 if (AliTPCReconstructor::GetRecoParam()->GetUseHLTClusters() == 3 ||
375 AliTPCReconstructor::GetRecoParam()->GetUseHLTClusters() == 4 ) {
376 if(!AliTPCReconstructor::GetRecoParam()->GetUseOnePadCluster())
377 if (TMath::Abs(cluster->GetSigmaY2()) < kAlmost0)
387 //_____________________________________________________________________________
388 AliTPCtrackerMI::AliTPCtrackerMI(const AliTPCParam *par):
390 fkNIS(par->GetNInnerSector()/2),
392 fkNOS(par->GetNOuterSector()/2),
409 //---------------------------------------------------------------------
410 // The main TPC tracker constructor
411 //---------------------------------------------------------------------
412 fInnerSec=new AliTPCtrackerSector[fkNIS];
413 fOuterSec=new AliTPCtrackerSector[fkNOS];
416 for (i=0; i<fkNIS; i++) fInnerSec[i].Setup(par,0);
417 for (i=0; i<fkNOS; i++) fOuterSec[i].Setup(par,1);
420 Int_t nrowlow = par->GetNRowLow();
421 Int_t nrowup = par->GetNRowUp();
424 for (i=0;i<nrowlow;i++){
425 fXRow[i] = par->GetPadRowRadiiLow(i);
426 fPadLength[i]= par->GetPadPitchLength(0,i);
427 fYMax[i] = fXRow[i]*TMath::Tan(0.5*par->GetInnerAngle());
431 for (i=0;i<nrowup;i++){
432 fXRow[i+nrowlow] = par->GetPadRowRadiiUp(i);
433 fPadLength[i+nrowlow] = par->GetPadPitchLength(60,i);
434 fYMax[i+nrowlow] = fXRow[i+nrowlow]*TMath::Tan(0.5*par->GetOuterAngle());
437 if (AliTPCReconstructor::StreamLevel()>0) {
438 fDebugStreamer = new TTreeSRedirector("TPCdebug.root");
441 //________________________________________________________________________
442 AliTPCtrackerMI::AliTPCtrackerMI(const AliTPCtrackerMI &t):
463 //------------------------------------
464 // dummy copy constructor
465 //------------------------------------------------------------------
467 for (Int_t irow=0; irow<200; irow++){
474 AliTPCtrackerMI & AliTPCtrackerMI::operator=(const AliTPCtrackerMI& /*r*/)
476 //------------------------------
478 //--------------------------------------------------------------
481 //_____________________________________________________________________________
482 AliTPCtrackerMI::~AliTPCtrackerMI() {
483 //------------------------------------------------------------------
484 // TPC tracker destructor
485 //------------------------------------------------------------------
492 if (fDebugStreamer) delete fDebugStreamer;
496 void AliTPCtrackerMI::FillESD(const TObjArray* arr)
500 //fill esds using updated tracks
503 // write tracks to the event
504 // store index of the track
505 Int_t nseed=arr->GetEntriesFast();
506 //FindKinks(arr,fEvent);
507 for (Int_t i=0; i<nseed; i++) {
508 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
512 if (AliTPCReconstructor::StreamLevel()>1) {
513 (*fDebugStreamer)<<"Track0"<<
517 // pt->PropagateTo(fkParam->GetInnerRadiusLow());
518 if (pt->GetKinkIndex(0)<=0){ //don't propagate daughter tracks
519 pt->PropagateTo(fkParam->GetInnerRadiusLow());
522 if (( pt->GetPoints()[2]- pt->GetPoints()[0])>5 && pt->GetPoints()[3]>0.8){
524 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
525 iotrack.SetTPCPoints(pt->GetPoints());
526 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
527 iotrack.SetV0Indexes(pt->GetV0Indexes());
528 // iotrack.SetTPCpid(pt->fTPCr);
529 //iotrack.SetTPCindex(i);
530 fEvent->AddTrack(&iotrack);
534 if ( (pt->GetNumberOfClusters()>70)&& (Float_t(pt->GetNumberOfClusters())/Float_t(pt->GetNFoundable()))>0.55) {
536 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
537 iotrack.SetTPCPoints(pt->GetPoints());
538 //iotrack.SetTPCindex(i);
539 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
540 iotrack.SetV0Indexes(pt->GetV0Indexes());
541 // iotrack.SetTPCpid(pt->fTPCr);
542 fEvent->AddTrack(&iotrack);
546 // short tracks - maybe decays
548 if ( (pt->GetNumberOfClusters()>30) && (Float_t(pt->GetNumberOfClusters())/Float_t(pt->GetNFoundable()))>0.70) {
549 Int_t found,foundable,shared;
550 pt->GetClusterStatistic(0,60,found, foundable,shared,kFALSE);
551 if ( (found>20) && (pt->GetNShared()/float(pt->GetNumberOfClusters())<0.2)){
553 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
554 //iotrack.SetTPCindex(i);
555 iotrack.SetTPCPoints(pt->GetPoints());
556 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
557 iotrack.SetV0Indexes(pt->GetV0Indexes());
558 //iotrack.SetTPCpid(pt->fTPCr);
559 fEvent->AddTrack(&iotrack);
564 if ( (pt->GetNumberOfClusters()>20) && (Float_t(pt->GetNumberOfClusters())/Float_t(pt->GetNFoundable()))>0.8) {
565 Int_t found,foundable,shared;
566 pt->GetClusterStatistic(0,60,found, foundable,shared,kFALSE);
567 if (found<20) continue;
568 if (pt->GetNShared()/float(pt->GetNumberOfClusters())>0.2) continue;
571 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
572 iotrack.SetTPCPoints(pt->GetPoints());
573 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
574 iotrack.SetV0Indexes(pt->GetV0Indexes());
575 //iotrack.SetTPCpid(pt->fTPCr);
576 //iotrack.SetTPCindex(i);
577 fEvent->AddTrack(&iotrack);
580 // short tracks - secondaties
582 if ( (pt->GetNumberOfClusters()>30) ) {
583 Int_t found,foundable,shared;
584 pt->GetClusterStatistic(128,158,found, foundable,shared,kFALSE);
585 if ( (found>20) && (pt->GetNShared()/float(pt->GetNumberOfClusters())<0.2) &&float(found)/float(foundable)>0.8){
587 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
588 iotrack.SetTPCPoints(pt->GetPoints());
589 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
590 iotrack.SetV0Indexes(pt->GetV0Indexes());
591 //iotrack.SetTPCpid(pt->fTPCr);
592 //iotrack.SetTPCindex(i);
593 fEvent->AddTrack(&iotrack);
598 if ( (pt->GetNumberOfClusters()>15)) {
599 Int_t found,foundable,shared;
600 pt->GetClusterStatistic(138,158,found, foundable,shared,kFALSE);
601 if (found<15) continue;
602 if (foundable<=0) continue;
603 if (pt->GetNShared()/float(pt->GetNumberOfClusters())>0.2) continue;
604 if (float(found)/float(foundable)<0.8) continue;
607 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
608 iotrack.SetTPCPoints(pt->GetPoints());
609 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
610 iotrack.SetV0Indexes(pt->GetV0Indexes());
611 // iotrack.SetTPCpid(pt->fTPCr);
612 //iotrack.SetTPCindex(i);
613 fEvent->AddTrack(&iotrack);
617 // >> account for suppressed tracks in the kink indices (RS)
618 int nESDtracks = fEvent->GetNumberOfTracks();
619 for (int it=nESDtracks;it--;) {
620 AliESDtrack* esdTr = fEvent->GetTrack(it);
621 if (!esdTr || !esdTr->GetKinkIndex(0)) continue;
622 for (int ik=0;ik<3;ik++) {
624 if (!(knkId=esdTr->GetKinkIndex(ik))) break; // no more kinks for this track
625 AliESDkink* kink = fEvent->GetKink(TMath::Abs(knkId)-1);
627 AliError(Form("ESDTrack%d refers to non-existing kink %d",it,TMath::Abs(knkId)-1));
630 kink->SetIndex(it, knkId<0 ? 0:1); // update track index of the kink: mother at 0, daughter at 1
633 // << account for suppressed tracks in the kink indices (RS)
634 AliInfo(Form("Number of filled ESDs-\t%d\n",fEvent->GetNumberOfTracks()));
642 Double_t AliTPCtrackerMI::ErrY2(AliTPCseed* seed, const AliTPCclusterMI * cl){
645 // Use calibrated cluster error from OCDB
647 AliTPCClusterParam * clparam = AliTPCcalibDB::Instance()->GetClusterParam();
649 Float_t z = TMath::Abs(fkParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
650 Int_t ctype = cl->GetType();
651 Int_t type = (cl->GetRow()<63) ? 0: (cl->GetRow()>126) ? 1:2;
652 Double_t angle = seed->GetSnp()*seed->GetSnp();
653 angle = TMath::Sqrt(TMath::Abs(angle/(1.-angle)));
654 Double_t erry2 = clparam->GetError0Par(0,type, z,angle);
656 erry2+=0.5; // edge cluster
659 seed->SetErrorY2(erry2);
663 //calculate look-up table at the beginning
664 // static Bool_t ginit = kFALSE;
665 // static Float_t gnoise1,gnoise2,gnoise3;
666 // static Float_t ggg1[10000];
667 // static Float_t ggg2[10000];
668 // static Float_t ggg3[10000];
669 // static Float_t glandau1[10000];
670 // static Float_t glandau2[10000];
671 // static Float_t glandau3[10000];
673 // static Float_t gcor01[500];
674 // static Float_t gcor02[500];
675 // static Float_t gcorp[500];
679 // if (ginit==kFALSE){
680 // for (Int_t i=1;i<500;i++){
681 // Float_t rsigma = float(i)/100.;
682 // gcor02[i] = TMath::Max(0.78 +TMath::Exp(7.4*(rsigma-1.2)),0.6);
683 // gcor01[i] = TMath::Max(0.72 +TMath::Exp(3.36*(rsigma-1.2)),0.6);
684 // gcorp[i] = TMath::Max(TMath::Power((rsigma+0.5),1.5),1.2);
688 // for (Int_t i=3;i<10000;i++){
692 // Float_t amp = float(i);
693 // Float_t padlength =0.75;
694 // gnoise1 = 0.0004/padlength;
695 // Float_t nel = 0.268*amp;
696 // Float_t nprim = 0.155*amp;
697 // ggg1[i] = fkParam->GetDiffT()*fkParam->GetDiffT()*(2+0.001*nel/(padlength*padlength))/nel;
698 // glandau1[i] = (2.+0.12*nprim)*0.5* (2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
699 // if (glandau1[i]>1) glandau1[i]=1;
700 // glandau1[i]*=padlength*padlength/12.;
704 // gnoise2 = 0.0004/padlength;
706 // nprim = 0.133*amp;
707 // ggg2[i] = fkParam->GetDiffT()*fkParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
708 // glandau2[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
709 // if (glandau2[i]>1) glandau2[i]=1;
710 // glandau2[i]*=padlength*padlength/12.;
715 // gnoise3 = 0.0004/padlength;
717 // nprim = 0.133*amp;
718 // ggg3[i] = fkParam->GetDiffT()*fkParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
719 // glandau3[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
720 // if (glandau3[i]>1) glandau3[i]=1;
721 // glandau3[i]*=padlength*padlength/12.;
729 // Int_t amp = int(TMath::Abs(cl->GetQ()));
731 // seed->SetErrorY2(1.);
735 // Float_t z = TMath::Abs(fkParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
736 // Int_t ctype = cl->GetType();
737 // Float_t padlength= GetPadPitchLength(seed->GetRow());
738 // Double_t angle2 = seed->GetSnp()*seed->GetSnp();
739 // angle2 = angle2/(1-angle2);
741 // //cluster "quality"
742 // Int_t rsigmay = int(100.*cl->GetSigmaY2()/(seed->GetCurrentSigmaY2()));
745 // if (fSectors==fInnerSec){
746 // snoise2 = gnoise1;
747 // res = ggg1[amp]*z+glandau1[amp]*angle2;
748 // if (ctype==0) res *= gcor01[rsigmay];
751 // res*= gcorp[rsigmay];
755 // if (padlength<1.1){
756 // snoise2 = gnoise2;
757 // res = ggg2[amp]*z+glandau2[amp]*angle2;
758 // if (ctype==0) res *= gcor02[rsigmay];
761 // res*= gcorp[rsigmay];
765 // snoise2 = gnoise3;
766 // res = ggg3[amp]*z+glandau3[amp]*angle2;
767 // if (ctype==0) res *= gcor02[rsigmay];
770 // res*= gcorp[rsigmay];
777 // res*=2.4; // overestimate error 2 times
781 // if (res<2*snoise2)
784 // seed->SetErrorY2(res);
792 Double_t AliTPCtrackerMI::ErrZ2(AliTPCseed* seed, const AliTPCclusterMI * cl){
795 // Use calibrated cluster error from OCDB
797 AliTPCClusterParam * clparam = AliTPCcalibDB::Instance()->GetClusterParam();
799 Float_t z = TMath::Abs(fkParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
800 Int_t ctype = cl->GetType();
801 Int_t type = (cl->GetRow()<63) ? 0: (cl->GetRow()>126) ? 1:2;
803 Double_t angle2 = seed->GetSnp()*seed->GetSnp();
804 angle2 = seed->GetTgl()*seed->GetTgl()*(1+angle2/(1-angle2));
805 Double_t angle = TMath::Sqrt(TMath::Abs(angle2));
806 Double_t errz2 = clparam->GetError0Par(1,type, z,angle);
808 errz2+=0.5; // edge cluster
811 seed->SetErrorZ2(errz2);
817 // //seed->SetErrorY2(0.1);
819 // //calculate look-up table at the beginning
820 // static Bool_t ginit = kFALSE;
821 // static Float_t gnoise1,gnoise2,gnoise3;
822 // static Float_t ggg1[10000];
823 // static Float_t ggg2[10000];
824 // static Float_t ggg3[10000];
825 // static Float_t glandau1[10000];
826 // static Float_t glandau2[10000];
827 // static Float_t glandau3[10000];
829 // static Float_t gcor01[1000];
830 // static Float_t gcor02[1000];
831 // static Float_t gcorp[1000];
835 // if (ginit==kFALSE){
836 // for (Int_t i=1;i<1000;i++){
837 // Float_t rsigma = float(i)/100.;
838 // gcor02[i] = TMath::Max(0.81 +TMath::Exp(6.8*(rsigma-1.2)),0.6);
839 // gcor01[i] = TMath::Max(0.72 +TMath::Exp(2.04*(rsigma-1.2)),0.6);
840 // gcorp[i] = TMath::Max(TMath::Power((rsigma+0.5),1.5),1.2);
844 // for (Int_t i=3;i<10000;i++){
848 // Float_t amp = float(i);
849 // Float_t padlength =0.75;
850 // gnoise1 = 0.0004/padlength;
851 // Float_t nel = 0.268*amp;
852 // Float_t nprim = 0.155*amp;
853 // ggg1[i] = fkParam->GetDiffT()*fkParam->GetDiffT()*(2+0.001*nel/(padlength*padlength))/nel;
854 // glandau1[i] = (2.+0.12*nprim)*0.5* (2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
855 // if (glandau1[i]>1) glandau1[i]=1;
856 // glandau1[i]*=padlength*padlength/12.;
860 // gnoise2 = 0.0004/padlength;
862 // nprim = 0.133*amp;
863 // ggg2[i] = fkParam->GetDiffT()*fkParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
864 // glandau2[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
865 // if (glandau2[i]>1) glandau2[i]=1;
866 // glandau2[i]*=padlength*padlength/12.;
871 // gnoise3 = 0.0004/padlength;
873 // nprim = 0.133*amp;
874 // ggg3[i] = fkParam->GetDiffT()*fkParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
875 // glandau3[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
876 // if (glandau3[i]>1) glandau3[i]=1;
877 // glandau3[i]*=padlength*padlength/12.;
885 // Int_t amp = int(TMath::Abs(cl->GetQ()));
887 // seed->SetErrorY2(1.);
891 // Float_t z = TMath::Abs(fkParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
892 // Int_t ctype = cl->GetType();
893 // Float_t padlength= GetPadPitchLength(seed->GetRow());
895 // Double_t angle2 = seed->GetSnp()*seed->GetSnp();
896 // // if (angle2<0.6) angle2 = 0.6;
897 // angle2 = seed->GetTgl()*seed->GetTgl()*(1+angle2/(1-angle2));
899 // //cluster "quality"
900 // Int_t rsigmaz = int(100.*cl->GetSigmaZ2()/(seed->GetCurrentSigmaZ2()));
903 // if (fSectors==fInnerSec){
904 // snoise2 = gnoise1;
905 // res = ggg1[amp]*z+glandau1[amp]*angle2;
906 // if (ctype==0) res *= gcor01[rsigmaz];
909 // res*= gcorp[rsigmaz];
913 // if (padlength<1.1){
914 // snoise2 = gnoise2;
915 // res = ggg2[amp]*z+glandau2[amp]*angle2;
916 // if (ctype==0) res *= gcor02[rsigmaz];
919 // res*= gcorp[rsigmaz];
923 // snoise2 = gnoise3;
924 // res = ggg3[amp]*z+glandau3[amp]*angle2;
925 // if (ctype==0) res *= gcor02[rsigmaz];
928 // res*= gcorp[rsigmaz];
937 // if ((ctype<0) &&<70){
942 // if (res<2*snoise2)
944 // if (res>3) res =3;
945 // seed->SetErrorZ2(res);
953 void AliTPCtrackerMI::RotateToLocal(AliTPCseed *seed)
955 //rotate to track "local coordinata
956 Float_t x = seed->GetX();
957 Float_t y = seed->GetY();
958 Float_t ymax = x*TMath::Tan(0.5*fSectors->GetAlpha());
961 seed->SetRelativeSector((seed->GetRelativeSector()+1) % fN);
962 if (!seed->Rotate(fSectors->GetAlpha()))
964 } else if (y <-ymax) {
965 seed->SetRelativeSector((seed->GetRelativeSector()-1+fN) % fN);
966 if (!seed->Rotate(-fSectors->GetAlpha()))
974 //_____________________________________________________________________________
975 Double_t AliTPCtrackerMI::F1old(Double_t x1,Double_t y1,
976 Double_t x2,Double_t y2,
977 Double_t x3,Double_t y3) const
979 //-----------------------------------------------------------------
980 // Initial approximation of the track curvature
981 //-----------------------------------------------------------------
982 Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
983 Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
984 (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
985 Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
986 (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
988 Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
989 if ( xr*xr+yr*yr<=0.00000000000001) return 100;
990 return -xr*yr/sqrt(xr*xr+yr*yr);
995 //_____________________________________________________________________________
996 Double_t AliTPCtrackerMI::F1(Double_t x1,Double_t y1,
997 Double_t x2,Double_t y2,
998 Double_t x3,Double_t y3) const
1000 //-----------------------------------------------------------------
1001 // Initial approximation of the track curvature
1002 //-----------------------------------------------------------------
1008 Double_t det = x3*y2-x2*y3;
1009 if (TMath::Abs(det)<1e-10){
1013 Double_t u = 0.5* (x2*(x2-x3)+y2*(y2-y3))/det;
1014 Double_t x0 = x3*0.5-y3*u;
1015 Double_t y0 = y3*0.5+x3*u;
1016 Double_t c2 = 1/TMath::Sqrt(x0*x0+y0*y0);
1022 Double_t AliTPCtrackerMI::F2(Double_t x1,Double_t y1,
1023 Double_t x2,Double_t y2,
1024 Double_t x3,Double_t y3) const
1026 //-----------------------------------------------------------------
1027 // Initial approximation of the track curvature
1028 //-----------------------------------------------------------------
1034 Double_t det = x3*y2-x2*y3;
1035 if (TMath::Abs(det)<1e-10) {
1039 Double_t u = 0.5* (x2*(x2-x3)+y2*(y2-y3))/det;
1040 Double_t x0 = x3*0.5-y3*u;
1041 Double_t y0 = y3*0.5+x3*u;
1042 Double_t c2 = 1/TMath::Sqrt(x0*x0+y0*y0);
1051 //_____________________________________________________________________________
1052 Double_t AliTPCtrackerMI::F2old(Double_t x1,Double_t y1,
1053 Double_t x2,Double_t y2,
1054 Double_t x3,Double_t y3) const
1056 //-----------------------------------------------------------------
1057 // Initial approximation of the track curvature times center of curvature
1058 //-----------------------------------------------------------------
1059 Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
1060 Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
1061 (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
1062 Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
1063 (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
1065 Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
1067 return -a/(d*y1-b)*xr/sqrt(xr*xr+yr*yr);
1070 //_____________________________________________________________________________
1071 Double_t AliTPCtrackerMI::F3(Double_t x1,Double_t y1,
1072 Double_t x2,Double_t y2,
1073 Double_t z1,Double_t z2) const
1075 //-----------------------------------------------------------------
1076 // Initial approximation of the tangent of the track dip angle
1077 //-----------------------------------------------------------------
1078 return (z1 - z2)/sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
1082 Double_t AliTPCtrackerMI::F3n(Double_t x1,Double_t y1,
1083 Double_t x2,Double_t y2,
1084 Double_t z1,Double_t z2, Double_t c) const
1086 //-----------------------------------------------------------------
1087 // Initial approximation of the tangent of the track dip angle
1088 //-----------------------------------------------------------------
1092 //angle1 = (z1-z2)*c/(TMath::ASin(c*x1-ni)-TMath::ASin(c*x2-ni));
1094 Double_t d = TMath::Sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
1095 if (TMath::Abs(d*c*0.5)>1) return 0;
1096 // Double_t angle2 = TMath::ASin(d*c*0.5);
1097 // Double_t angle2 = AliTPCFastMath::FastAsin(d*c*0.5);
1098 Double_t angle2 = (d*c*0.5>0.1)? TMath::ASin(d*c*0.5): AliTPCFastMath::FastAsin(d*c*0.5);
1100 angle2 = (z1-z2)*c/(angle2*2.);
1104 Bool_t AliTPCtrackerMI::GetProlongation(Double_t x1, Double_t x2, Double_t x[5], Double_t &y, Double_t &z) const
1105 {//-----------------------------------------------------------------
1106 // This function find proloncation of a track to a reference plane x=x2.
1107 //-----------------------------------------------------------------
1111 if (TMath::Abs(x[4]*x1 - x[2]) >= 0.999) {
1115 Double_t c1=x[4]*x1 - x[2], r1=TMath::Sqrt((1.-c1)*(1.+c1));
1116 Double_t c2=x[4]*x2 - x[2], r2=TMath::Sqrt((1.-c2)*(1.+c2));
1120 Double_t dy = dx*(c1+c2)/(r1+r2);
1123 Double_t delta = x[4]*dx*(c1+c2)/(c1*r2 + c2*r1);
1125 if (TMath::Abs(delta)>0.01){
1126 dz = x[3]*TMath::ASin(delta)/x[4];
1128 dz = x[3]*AliTPCFastMath::FastAsin(delta)/x[4];
1131 //dz = x[3]*AliTPCFastMath::FastAsin(delta)/x[4];
1139 Int_t AliTPCtrackerMI::LoadClusters (TTree *const tree)
1144 return LoadClusters();
1148 Int_t AliTPCtrackerMI::LoadClusters(const TObjArray *arr)
1151 // load clusters to the memory
1152 AliTPCClustersRow *clrow = new AliTPCClustersRow("AliTPCclusterMI");
1153 Int_t lower = arr->LowerBound();
1154 Int_t entries = arr->GetEntriesFast();
1156 for (Int_t i=lower; i<entries; i++) {
1157 clrow = (AliTPCClustersRow*) arr->At(i);
1158 if(!clrow) continue;
1159 if(!clrow->GetArray()) continue;
1163 fkParam->AdjustSectorRow(clrow->GetID(),sec,row);
1165 for (Int_t icl=0; icl<clrow->GetArray()->GetEntriesFast(); icl++){
1166 Transform((AliTPCclusterMI*)(clrow->GetArray()->At(icl)));
1169 if (clrow->GetArray()->GetEntriesFast()<=0) continue;
1170 AliTPCtrackerRow * tpcrow=0;
1173 tpcrow = &(fInnerSec[sec%fkNIS][row]);
1177 tpcrow = &(fOuterSec[(sec-fkNIS*2)%fkNOS][row]);
1178 left = (sec-fkNIS*2)/fkNOS;
1181 tpcrow->SetN1(clrow->GetArray()->GetEntriesFast());
1182 tpcrow->SetClusters1(new AliTPCclusterMI[tpcrow->GetN1()]);
1183 for (Int_t j=0;j<tpcrow->GetN1();++j)
1184 tpcrow->SetCluster1(j, *(AliTPCclusterMI*)(clrow->GetArray()->At(j)));
1187 tpcrow->SetN2(clrow->GetArray()->GetEntriesFast());
1188 tpcrow->SetClusters2(new AliTPCclusterMI[tpcrow->GetN2()]);
1189 for (Int_t j=0;j<tpcrow->GetN2();++j)
1190 tpcrow->SetCluster2(j,*(AliTPCclusterMI*)(clrow->GetArray()->At(j)));
1192 clrow->GetArray()->Clear("C");
1201 Int_t AliTPCtrackerMI::LoadClusters(const TClonesArray *arr)
1204 // load clusters to the memory from one
1207 AliTPCclusterMI *clust=0;
1208 Int_t count[72][96] = { {0} , {0} };
1210 // loop over clusters
1211 for (Int_t icl=0; icl<arr->GetEntriesFast(); icl++) {
1212 clust = (AliTPCclusterMI*)arr->At(icl);
1213 if(!clust) continue;
1214 //printf("cluster: det %d, row %d \n", clust->GetDetector(),clust->GetRow());
1216 // transform clusters
1219 // count clusters per pad row
1220 count[clust->GetDetector()][clust->GetRow()]++;
1223 // insert clusters to sectors
1224 for (Int_t icl=0; icl<arr->GetEntriesFast(); icl++) {
1225 clust = (AliTPCclusterMI*)arr->At(icl);
1226 if(!clust) continue;
1228 Int_t sec = clust->GetDetector();
1229 Int_t row = clust->GetRow();
1231 // filter overlapping pad rows needed by HLT
1232 if(sec<fkNIS*2) { //IROCs
1233 if(row == 30) continue;
1236 if(row == 27 || row == 76) continue;
1242 fInnerSec[sec%fkNIS].InsertCluster(clust, count[sec][row], fkParam);
1245 left = (sec-fkNIS*2)/fkNOS;
1246 fOuterSec[(sec-fkNIS*2)%fkNOS].InsertCluster(clust, count[sec][row], fkParam);
1250 // Load functions must be called behind LoadCluster(TClonesArray*)
1252 //LoadOuterSectors();
1253 //LoadInnerSectors();
1259 Int_t AliTPCtrackerMI::LoadClusters()
1262 // load clusters to the memory
1263 AliTPCClustersRow *clrow= new AliTPCClustersRow("AliTPCclusterMI");
1265 // TTree * tree = fClustersArray.GetTree();
1267 TTree * tree = fInput;
1268 TBranch * br = tree->GetBranch("Segment");
1269 br->SetAddress(&clrow);
1271 Int_t j=Int_t(tree->GetEntries());
1272 for (Int_t i=0; i<j; i++) {
1276 fkParam->AdjustSectorRow(clrow->GetID(),sec,row);
1277 for (Int_t icl=0; icl<clrow->GetArray()->GetEntriesFast(); icl++){
1278 Transform((AliTPCclusterMI*)(clrow->GetArray()->At(icl)));
1281 AliTPCtrackerRow * tpcrow=0;
1284 tpcrow = &(fInnerSec[sec%fkNIS][row]);
1288 tpcrow = &(fOuterSec[(sec-fkNIS*2)%fkNOS][row]);
1289 left = (sec-fkNIS*2)/fkNOS;
1292 tpcrow->SetN1(clrow->GetArray()->GetEntriesFast());
1293 tpcrow->SetClusters1(new AliTPCclusterMI[tpcrow->GetN1()]);
1294 for (Int_t k=0;k<tpcrow->GetN1();++k)
1295 tpcrow->SetCluster1(k, *(AliTPCclusterMI*)(clrow->GetArray()->At(k)));
1298 tpcrow->SetN2(clrow->GetArray()->GetEntriesFast());
1299 tpcrow->SetClusters2(new AliTPCclusterMI[tpcrow->GetN2()]);
1300 for (Int_t k=0;k<tpcrow->GetN2();++k)
1301 tpcrow->SetCluster2(k,*(AliTPCclusterMI*)(clrow->GetArray()->At(k)));
1312 void AliTPCtrackerMI::UnloadClusters()
1315 // unload clusters from the memory
1317 Int_t nrows = fOuterSec->GetNRows();
1318 for (Int_t sec = 0;sec<fkNOS;sec++)
1319 for (Int_t row = 0;row<nrows;row++){
1320 AliTPCtrackerRow* tpcrow = &(fOuterSec[sec%fkNOS][row]);
1322 // if (tpcrow->fClusters1) delete []tpcrow->fClusters1;
1323 // if (tpcrow->fClusters2) delete []tpcrow->fClusters2;
1325 tpcrow->ResetClusters();
1328 nrows = fInnerSec->GetNRows();
1329 for (Int_t sec = 0;sec<fkNIS;sec++)
1330 for (Int_t row = 0;row<nrows;row++){
1331 AliTPCtrackerRow* tpcrow = &(fInnerSec[sec%fkNIS][row]);
1333 // if (tpcrow->fClusters1) delete []tpcrow->fClusters1;
1334 //if (tpcrow->fClusters2) delete []tpcrow->fClusters2;
1336 tpcrow->ResetClusters();
1342 void AliTPCtrackerMI::FillClusterArray(TObjArray* array) const{
1344 // Filling cluster to the array - For visualization purposes
1347 nrows = fOuterSec->GetNRows();
1348 for (Int_t sec = 0;sec<fkNOS;sec++)
1349 for (Int_t row = 0;row<nrows;row++){
1350 AliTPCtrackerRow* tpcrow = &(fOuterSec[sec%fkNOS][row]);
1351 if (!tpcrow) continue;
1352 for (Int_t icl = 0;icl<tpcrow->GetN();icl++){
1353 array->AddLast((TObject*)((*tpcrow)[icl]));
1356 nrows = fInnerSec->GetNRows();
1357 for (Int_t sec = 0;sec<fkNIS;sec++)
1358 for (Int_t row = 0;row<nrows;row++){
1359 AliTPCtrackerRow* tpcrow = &(fInnerSec[sec%fkNIS][row]);
1360 if (!tpcrow) continue;
1361 for (Int_t icl = 0;icl<tpcrow->GetN();icl++){
1362 array->AddLast((TObject*)(*tpcrow)[icl]);
1368 void AliTPCtrackerMI::Transform(AliTPCclusterMI * cluster){
1372 AliTPCcalibDB * calibDB = AliTPCcalibDB::Instance();
1373 AliTPCTransform *transform = calibDB->GetTransform() ;
1375 AliFatal("Tranformations not in calibDB");
1378 transform->SetCurrentRecoParam((AliTPCRecoParam*)AliTPCReconstructor::GetRecoParam());
1379 Double_t x[3]={cluster->GetRow(),cluster->GetPad(),cluster->GetTimeBin()};
1380 Int_t i[1]={cluster->GetDetector()};
1381 transform->Transform(x,i,0,1);
1382 // if (cluster->GetDetector()%36>17){
1387 // in debug mode check the transformation
1389 if (AliTPCReconstructor::StreamLevel()>2) {
1391 cluster->GetGlobalXYZ(gx);
1392 Int_t event = (fEvent==NULL)? 0: fEvent->GetEventNumberInFile();
1393 TTreeSRedirector &cstream = *fDebugStreamer;
1394 cstream<<"Transform"<<
1405 cluster->SetX(x[0]);
1406 cluster->SetY(x[1]);
1407 cluster->SetZ(x[2]);
1412 //if (!fkParam->IsGeoRead()) fkParam->ReadGeoMatrices();
1413 if (AliTPCReconstructor::GetRecoParam()->GetUseSectorAlignment() && (!calibDB->HasAlignmentOCDB())){
1414 TGeoHMatrix *mat = fkParam->GetClusterMatrix(cluster->GetDetector());
1416 Double_t pos[3]= {cluster->GetX(),cluster->GetY(),cluster->GetZ()};
1417 Double_t posC[3]={cluster->GetX(),cluster->GetY(),cluster->GetZ()};
1418 if (mat) mat->LocalToMaster(pos,posC);
1420 // chack Loading of Geo matrices from GeoManager - TEMPORARY FIX
1422 cluster->SetX(posC[0]);
1423 cluster->SetY(posC[1]);
1424 cluster->SetZ(posC[2]);
1428 //_____________________________________________________________________________
1429 Int_t AliTPCtrackerMI::LoadOuterSectors() {
1430 //-----------------------------------------------------------------
1431 // This function fills outer TPC sectors with clusters.
1432 //-----------------------------------------------------------------
1433 Int_t nrows = fOuterSec->GetNRows();
1435 for (Int_t sec = 0;sec<fkNOS;sec++)
1436 for (Int_t row = 0;row<nrows;row++){
1437 AliTPCtrackerRow* tpcrow = &(fOuterSec[sec%fkNOS][row]);
1438 Int_t sec2 = sec+2*fkNIS;
1440 Int_t ncl = tpcrow->GetN1();
1442 AliTPCclusterMI *c= (tpcrow->GetCluster1(ncl));
1443 index=(((sec2<<8)+row)<<16)+ncl;
1444 tpcrow->InsertCluster(c,index);
1447 ncl = tpcrow->GetN2();
1449 AliTPCclusterMI *c= (tpcrow->GetCluster2(ncl));
1450 index=((((sec2+fkNOS)<<8)+row)<<16)+ncl;
1451 tpcrow->InsertCluster(c,index);
1454 // write indexes for fast acces
1456 for (Int_t i=0;i<510;i++)
1457 tpcrow->SetFastCluster(i,-1);
1458 for (Int_t i=0;i<tpcrow->GetN();i++){
1459 Int_t zi = Int_t((*tpcrow)[i]->GetZ()+255.);
1460 tpcrow->SetFastCluster(zi,i); // write index
1463 for (Int_t i=0;i<510;i++){
1464 if (tpcrow->GetFastCluster(i)<0)
1465 tpcrow->SetFastCluster(i,last);
1467 last = tpcrow->GetFastCluster(i);
1476 //_____________________________________________________________________________
1477 Int_t AliTPCtrackerMI::LoadInnerSectors() {
1478 //-----------------------------------------------------------------
1479 // This function fills inner TPC sectors with clusters.
1480 //-----------------------------------------------------------------
1481 Int_t nrows = fInnerSec->GetNRows();
1483 for (Int_t sec = 0;sec<fkNIS;sec++)
1484 for (Int_t row = 0;row<nrows;row++){
1485 AliTPCtrackerRow* tpcrow = &(fInnerSec[sec%fkNIS][row]);
1488 Int_t ncl = tpcrow->GetN1();
1490 AliTPCclusterMI *c= (tpcrow->GetCluster1(ncl));
1491 index=(((sec<<8)+row)<<16)+ncl;
1492 tpcrow->InsertCluster(c,index);
1495 ncl = tpcrow->GetN2();
1497 AliTPCclusterMI *c= (tpcrow->GetCluster2(ncl));
1498 index=((((sec+fkNIS)<<8)+row)<<16)+ncl;
1499 tpcrow->InsertCluster(c,index);
1502 // write indexes for fast acces
1504 for (Int_t i=0;i<510;i++)
1505 tpcrow->SetFastCluster(i,-1);
1506 for (Int_t i=0;i<tpcrow->GetN();i++){
1507 Int_t zi = Int_t((*tpcrow)[i]->GetZ()+255.);
1508 tpcrow->SetFastCluster(zi,i); // write index
1511 for (Int_t i=0;i<510;i++){
1512 if (tpcrow->GetFastCluster(i)<0)
1513 tpcrow->SetFastCluster(i,last);
1515 last = tpcrow->GetFastCluster(i);
1527 //_________________________________________________________________________
1528 AliTPCclusterMI *AliTPCtrackerMI::GetClusterMI(Int_t index) const {
1529 //--------------------------------------------------------------------
1530 // Return pointer to a given cluster
1531 //--------------------------------------------------------------------
1532 if (index<0) return 0; // no cluster
1533 Int_t sec=(index&0xff000000)>>24;
1534 Int_t row=(index&0x00ff0000)>>16;
1535 Int_t ncl=(index&0x00007fff)>>00;
1537 const AliTPCtrackerRow * tpcrow=0;
1538 AliTPCclusterMI * clrow =0;
1540 if (sec<0 || sec>=fkNIS*4) {
1541 AliWarning(Form("Wrong sector %d",sec));
1546 AliTPCtrackerSector& tracksec = fInnerSec[sec%fkNIS];
1547 if (tracksec.GetNRows()<=row) return 0;
1548 tpcrow = &(tracksec[row]);
1549 if (tpcrow==0) return 0;
1552 if (tpcrow->GetN1()<=ncl) return 0;
1553 clrow = tpcrow->GetClusters1();
1556 if (tpcrow->GetN2()<=ncl) return 0;
1557 clrow = tpcrow->GetClusters2();
1561 AliTPCtrackerSector& tracksec = fOuterSec[(sec-fkNIS*2)%fkNOS];
1562 if (tracksec.GetNRows()<=row) return 0;
1563 tpcrow = &(tracksec[row]);
1564 if (tpcrow==0) return 0;
1566 if (sec-2*fkNIS<fkNOS) {
1567 if (tpcrow->GetN1()<=ncl) return 0;
1568 clrow = tpcrow->GetClusters1();
1571 if (tpcrow->GetN2()<=ncl) return 0;
1572 clrow = tpcrow->GetClusters2();
1576 return &(clrow[ncl]);
1582 Int_t AliTPCtrackerMI::FollowToNext(AliTPCseed& t, Int_t nr) {
1583 //-----------------------------------------------------------------
1584 // This function tries to find a track prolongation to next pad row
1585 //-----------------------------------------------------------------
1587 Double_t x= GetXrow(nr), ymax=GetMaxY(nr);
1590 AliTPCclusterMI *cl=0;
1591 Int_t tpcindex= t.GetClusterIndex2(nr);
1593 // update current shape info every 5 pad-row
1594 // if ( (nr%5==0) || t.GetNumberOfClusters()<2 || (t.fCurrentSigmaY2<0.0001) ){
1598 if (fIteration>0 && tpcindex>=-1){ //if we have already clusters
1600 if (tpcindex==-1) return 0; //track in dead zone
1602 cl = t.GetClusterPointer(nr);
1603 if ( (cl==0) ) cl = GetClusterMI(tpcindex);
1604 t.SetCurrentClusterIndex1(tpcindex);
1607 Int_t relativesector = ((tpcindex&0xff000000)>>24)%18; // if previously accepted cluster in different sector
1608 Float_t angle = relativesector*fSectors->GetAlpha()+fSectors->GetAlphaShift();
1610 if (angle<-TMath::Pi()) angle += 2*TMath::Pi();
1611 if (angle>=TMath::Pi()) angle -= 2*TMath::Pi();
1613 if (TMath::Abs(angle-t.GetAlpha())>0.001){
1614 Double_t rotation = angle-t.GetAlpha();
1615 t.SetRelativeSector(relativesector);
1616 if (!t.Rotate(rotation)) return 0;
1618 if (!t.PropagateTo(x)) return 0;
1620 t.SetCurrentCluster(cl);
1622 Int_t accept = AcceptCluster(&t,t.GetCurrentCluster());
1623 if ((tpcindex&0x8000)==0) accept =0;
1625 //if founded cluster is acceptible
1626 if (cl->IsUsed(11)) { // id cluster is shared inrease uncertainty
1627 t.SetErrorY2(t.GetErrorY2()+0.03);
1628 t.SetErrorZ2(t.GetErrorZ2()+0.03);
1629 t.SetErrorY2(t.GetErrorY2()*3);
1630 t.SetErrorZ2(t.GetErrorZ2()*3);
1632 t.SetNFoundable(t.GetNFoundable()+1);
1633 UpdateTrack(&t,accept);
1638 if (TMath::Abs(t.GetSnp())>AliTPCReconstructor::GetMaxSnpTracker()) return 0; // cut on angle
1639 if (fIteration>1 && IsFindable(t)){
1640 // not look for new cluster during refitting
1641 t.SetNFoundable(t.GetNFoundable()+1);
1646 // if (TMath::Abs(t.GetSnp())>0.95 || TMath::Abs(x*t.GetC()-t.GetEta())>0.95) return 0;// patch 28 fev 06
1647 Double_t y=t.GetYat(x);
1648 if (TMath::Abs(y)>ymax){
1650 t.SetRelativeSector((t.GetRelativeSector()+1) % fN);
1651 if (!t.Rotate(fSectors->GetAlpha()))
1653 } else if (y <-ymax) {
1654 t.SetRelativeSector((t.GetRelativeSector()-1+fN) % fN);
1655 if (!t.Rotate(-fSectors->GetAlpha()))
1661 if (!t.PropagateTo(x)) {
1662 if (fIteration==0) t.SetRemoval(10);
1666 Double_t z=t.GetZ();
1669 if (!IsActive(t.GetRelativeSector(),nr)) {
1671 t.SetClusterIndex2(nr,-1);
1674 //AliInfo(Form("A - Sector%d phi %f - alpha %f", t.fRelativeSector,y/x, t.GetAlpha()));
1675 Bool_t isActive = IsActive(t.GetRelativeSector(),nr);
1676 Bool_t isActive2 = (nr>=fInnerSec->GetNRows()) ? fOuterSec[t.GetRelativeSector()][nr-fInnerSec->GetNRows()].GetN()>0:fInnerSec[t.GetRelativeSector()][nr].GetN()>0;
1678 if (!isActive || !isActive2) return 0;
1680 const AliTPCtrackerRow &krow=GetRow(t.GetRelativeSector(),nr);
1681 if ( (t.GetSigmaY2()<0) || t.GetSigmaZ2()<0) return 0;
1683 Double_t roadz = 1.;
1685 if (TMath::Abs(TMath::Abs(y)-ymax)<krow.GetDeadZone()){
1687 t.SetClusterIndex2(nr,-1);
1693 // if (TMath::Abs(z)<(AliTPCReconstructor::GetCtgRange()*x+10) && TMath::Abs(z)<fkParam->GetZLength(0) && (TMath::Abs(t.GetSnp())<AliTPCReconstructor::GetMaxSnpTracker()))
1694 t.SetNFoundable(t.GetNFoundable()+1);
1700 // cl = krow.FindNearest2(y+10.,z,roady,roadz,index);
1701 cl = krow.FindNearest2(y,z,roady,roadz,index);
1702 if (cl) t.SetCurrentClusterIndex1(krow.GetIndex(index));
1705 t.SetCurrentCluster(cl);
1707 if (fIteration==2&&cl->IsUsed(10)) return 0;
1708 Int_t accept = AcceptCluster(&t,t.GetCurrentCluster());
1709 if (fIteration==2&&cl->IsUsed(11)) {
1710 t.SetErrorY2(t.GetErrorY2()+0.03);
1711 t.SetErrorZ2(t.GetErrorZ2()+0.03);
1712 t.SetErrorY2(t.GetErrorY2()*3);
1713 t.SetErrorZ2(t.GetErrorZ2()*3);
1716 if (t.fCurrentCluster->IsUsed(10)){
1721 if (t.fNShared>0.7*t.GetNumberOfClusters()) {
1727 if (accept<3) UpdateTrack(&t,accept);
1730 if ( fIteration==0 && t.GetNFoundable()*0.5 > t.GetNumberOfClusters()) t.SetRemoval(10);
1738 //_________________________________________________________________________
1739 Bool_t AliTPCtrackerMI::GetTrackPoint(Int_t index, AliTrackPoint &p ) const
1741 // Get track space point by index
1742 // return false in case the cluster doesn't exist
1743 AliTPCclusterMI *cl = GetClusterMI(index);
1744 if (!cl) return kFALSE;
1745 Int_t sector = (index&0xff000000)>>24;
1746 // Int_t row = (index&0x00ff0000)>>16;
1748 // xyz[0] = fkParam->GetPadRowRadii(sector,row);
1749 xyz[0] = cl->GetX();
1750 xyz[1] = cl->GetY();
1751 xyz[2] = cl->GetZ();
1753 fkParam->AdjustCosSin(sector,cos,sin);
1754 Float_t x = cos*xyz[0]-sin*xyz[1];
1755 Float_t y = cos*xyz[1]+sin*xyz[0];
1757 Float_t sigmaY2 = 0.027*cl->GetSigmaY2();
1758 if (sector < fkParam->GetNInnerSector()) sigmaY2 *= 2.07;
1759 Float_t sigmaZ2 = 0.066*cl->GetSigmaZ2();
1760 if (sector < fkParam->GetNInnerSector()) sigmaZ2 *= 1.77;
1761 cov[0] = sin*sin*sigmaY2;
1762 cov[1] = -sin*cos*sigmaY2;
1764 cov[3] = cos*cos*sigmaY2;
1767 p.SetXYZ(x,y,xyz[2],cov);
1768 AliGeomManager::ELayerID iLayer;
1770 if (sector < fkParam->GetNInnerSector()) {
1771 iLayer = AliGeomManager::kTPC1;
1775 iLayer = AliGeomManager::kTPC2;
1776 idet = sector - fkParam->GetNInnerSector();
1778 UShort_t volid = AliGeomManager::LayerToVolUID(iLayer,idet);
1779 p.SetVolumeID(volid);
1785 Int_t AliTPCtrackerMI::UpdateClusters(AliTPCseed& t, Int_t nr) {
1786 //-----------------------------------------------------------------
1787 // This function tries to find a track prolongation to next pad row
1788 //-----------------------------------------------------------------
1789 t.SetCurrentCluster(0);
1790 t.SetCurrentClusterIndex1(0);
1792 Double_t xt=t.GetX();
1793 Int_t row = GetRowNumber(xt)-1;
1794 Double_t ymax= GetMaxY(nr);
1796 if (row < nr) return 1; // don't prolongate if not information until now -
1797 // if (TMath::Abs(t.GetSnp())>0.9 && t.GetNumberOfClusters()>40. && fIteration!=2) {
1799 // return 0; // not prolongate strongly inclined tracks
1801 // if (TMath::Abs(t.GetSnp())>0.95) {
1803 // return 0; // not prolongate strongly inclined tracks
1804 // }// patch 28 fev 06
1806 Double_t x= GetXrow(nr);
1808 //t.PropagateTo(x+0.02);
1809 //t.PropagateTo(x+0.01);
1810 if (!t.PropagateTo(x)){
1817 if (TMath::Abs(y)>ymax){
1819 t.SetRelativeSector((t.GetRelativeSector()+1) % fN);
1820 if (!t.Rotate(fSectors->GetAlpha()))
1822 } else if (y <-ymax) {
1823 t.SetRelativeSector((t.GetRelativeSector()-1+fN) % fN);
1824 if (!t.Rotate(-fSectors->GetAlpha()))
1827 // if (!t.PropagateTo(x)){
1834 if (TMath::Abs(t.GetSnp())>AliTPCReconstructor::GetMaxSnpTracker()) return 0;
1836 if (!IsActive(t.GetRelativeSector(),nr)) {
1838 t.SetClusterIndex2(nr,-1);
1841 //AliInfo(Form("A - Sector%d phi %f - alpha %f", t.fRelativeSector,y/x, t.GetAlpha()));
1843 AliTPCtrackerRow &krow=GetRow(t.GetRelativeSector(),nr);
1845 if (TMath::Abs(TMath::Abs(y)-ymax)<krow.GetDeadZone()){
1847 t.SetClusterIndex2(nr,-1);
1853 // if (TMath::Abs(t.GetZ())<(AliTPCReconstructor::GetCtgRange()*t.GetX()+10) && (TMath::Abs(t.GetSnp())<AliTPCReconstructor::GetMaxSnpTracker()))
1854 if (IsFindable(t)) t.SetNFoundable(t.GetNFoundable()+1);
1860 if ( (nr%2==0) || t.GetNumberOfClusters()<2){
1861 // t.fCurrentSigmaY = GetSigmaY(&t);
1862 //t.fCurrentSigmaZ = GetSigmaZ(&t);
1866 AliTPCclusterMI *cl=0;
1869 Double_t roady = 1.;
1870 Double_t roadz = 1.;
1874 index = t.GetClusterIndex2(nr);
1875 if ( (index>0) && (index&0x8000)==0){
1876 cl = t.GetClusterPointer(nr);
1877 if ( (cl==0) && (index>0)) cl = GetClusterMI(index);
1878 t.SetCurrentClusterIndex1(index);
1880 t.SetCurrentCluster(cl);
1886 // if (index<0) return 0;
1887 UInt_t uindex = TMath::Abs(index);
1890 //cl = krow.FindNearest2(y+10,z,roady,roadz,uindex);
1891 cl = krow.FindNearest2(y,z,roady,roadz,uindex);
1894 if (cl) t.SetCurrentClusterIndex1(krow.GetIndex(uindex));
1895 t.SetCurrentCluster(cl);
1901 Int_t AliTPCtrackerMI::FollowToNextCluster(AliTPCseed & t, Int_t nr) {
1902 //-----------------------------------------------------------------
1903 // This function tries to find a track prolongation to next pad row
1904 //-----------------------------------------------------------------
1906 //update error according neighborhoud
1908 if (t.GetCurrentCluster()) {
1910 Int_t accept = AcceptCluster(&t,t.GetCurrentCluster());
1912 if (t.GetCurrentCluster()->IsUsed(10)){
1917 t.SetNShared(t.GetNShared()+1);
1918 if (t.GetNShared()>0.7*t.GetNumberOfClusters()) {
1923 if (fIteration>0) accept = 0;
1924 if (accept<3) UpdateTrack(&t,accept);
1928 if ( ( (t.GetSigmaY2()+t.GetSigmaZ2())>0.16)&& t.GetNumberOfClusters()>18) t.SetRemoval(10);
1929 if ( t.GetChi2()/t.GetNumberOfClusters()>6 &&t.GetNumberOfClusters()>18) t.SetRemoval(10);
1931 if (( (t.GetNFoundable()*0.5 > t.GetNumberOfClusters()) || t.GetNoCluster()>15)) t.SetRemoval(10);
1939 //_____________________________________________________________________________
1940 Int_t AliTPCtrackerMI::FollowProlongation(AliTPCseed& t, Int_t rf, Int_t step) {
1941 //-----------------------------------------------------------------
1942 // This function tries to find a track prolongation.
1943 //-----------------------------------------------------------------
1944 Double_t xt=t.GetX();
1946 Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
1947 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
1948 if (alpha < 0. ) alpha += 2.*TMath::Pi();
1950 t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN);
1952 Int_t first = GetRowNumber(xt)-1;
1953 for (Int_t nr= first; nr>=rf; nr-=step) {
1955 if (t.GetKinkIndexes()[0]>0){
1956 for (Int_t i=0;i<3;i++){
1957 Int_t index = t.GetKinkIndexes()[i];
1958 if (index==0) break;
1959 if (index<0) continue;
1961 AliKink * kink = (AliKink*)fEvent->GetKink(index-1);
1963 printf("PROBLEM\n");
1966 Int_t kinkrow = kink->GetTPCRow0()+2+Int_t(0.5/(0.05+kink->GetAngle(2)));
1968 AliExternalTrackParam paramd(t);
1969 kink->SetDaughter(paramd);
1970 kink->SetStatus(2,5);
1977 if (nr==80) t.UpdateReference();
1978 if (nr<fInnerSec->GetNRows())
1979 fSectors = fInnerSec;
1981 fSectors = fOuterSec;
1982 if (FollowToNext(t,nr)==0)
1995 Int_t AliTPCtrackerMI::FollowBackProlongation(AliTPCseed& t, Int_t rf) {
1996 //-----------------------------------------------------------------
1997 // This function tries to find a track prolongation.
1998 //-----------------------------------------------------------------
2000 Double_t xt=t.GetX();
2001 Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
2002 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
2003 if (alpha < 0. ) alpha += 2.*TMath::Pi();
2004 t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN);
2006 Int_t first = t.GetFirstPoint();
2007 if (first<GetRowNumber(xt)+1) first = GetRowNumber(xt)+1;
2009 if (first<0) first=0;
2010 for (Int_t nr=first; nr<=rf; nr++) {
2011 // if ( (TMath::Abs(t.GetSnp())>0.95)) break;//patch 28 fev 06
2012 if (t.GetKinkIndexes()[0]<0){
2013 for (Int_t i=0;i<3;i++){
2014 Int_t index = t.GetKinkIndexes()[i];
2015 if (index==0) break;
2016 if (index>0) continue;
2017 index = TMath::Abs(index);
2018 AliKink * kink = (AliKink*)fEvent->GetKink(index-1);
2020 printf("PROBLEM\n");
2023 Int_t kinkrow = kink->GetTPCRow0()-2-Int_t(0.5/(0.05+kink->GetAngle(2)));
2025 AliExternalTrackParam paramm(t);
2026 kink->SetMother(paramm);
2027 kink->SetStatus(2,1);
2034 if (nr<fInnerSec->GetNRows())
2035 fSectors = fInnerSec;
2037 fSectors = fOuterSec;
2048 Float_t AliTPCtrackerMI::OverlapFactor(AliTPCseed * s1, AliTPCseed * s2, Int_t &sum1, Int_t & sum2)
2050 // overlapping factor
2056 Float_t dz2 =(s1->GetZ() - s2->GetZ());
2059 Float_t dy2 =TMath::Abs((s1->GetY() - s2->GetY()));
2061 Float_t distance = TMath::Sqrt(dz2+dy2);
2062 if (distance>4.) return 0; // if there are far away - not overlap - to reduce combinatorics
2065 Int_t firstpoint = TMath::Min(s1->GetFirstPoint(),s2->GetFirstPoint());
2066 Int_t lastpoint = TMath::Max(s1->GetLastPoint(),s2->GetLastPoint());
2071 if (firstpoint>lastpoint) {
2072 firstpoint =lastpoint;
2077 for (Int_t i=firstpoint-1;i<lastpoint+1;i++){
2078 if (s1->GetClusterIndex2(i)>0) sum1++;
2079 if (s2->GetClusterIndex2(i)>0) sum2++;
2080 if (s1->GetClusterIndex2(i)==s2->GetClusterIndex2(i) && s1->GetClusterIndex2(i)>0) {
2084 if (sum<5) return 0;
2086 Float_t summin = TMath::Min(sum1+1,sum2+1);
2087 Float_t ratio = (sum+1)/Float_t(summin);
2091 void AliTPCtrackerMI::SignShared(AliTPCseed * s1, AliTPCseed * s2)
2095 Float_t thetaCut = 0.2;//+10.*TMath::Sqrt(s1->GetSigmaTglZ()+ s2->GetSigmaTglZ());
2096 if (TMath::Abs(s1->GetTgl()-s2->GetTgl())>thetaCut) return;
2097 Float_t minCl = TMath::Min(s1->GetNumberOfClusters(),s2->GetNumberOfClusters());
2098 Int_t cutN0 = TMath::Max(5,TMath::Nint(0.1*minCl));
2103 //Int_t firstpoint = TMath::Max(s1->GetFirstPoint(),s2->GetFirstPoint());
2104 //Int_t lastpoint = TMath::Min(s1->GetLastPoint(),s2->GetLastPoint());
2105 Int_t firstpoint = 0;
2106 Int_t lastpoint = 160;
2108 // if (firstpoint>=lastpoint-5) return;;
2110 for (Int_t i=firstpoint;i<lastpoint;i++){
2111 // if ( (s1->GetClusterIndex2(i)&0xFFFF8FFF)==(s2->GetClusterIndex2(i)&0xFFFF8FFF) && s1->GetClusterIndex2(i)>0) {
2112 if ( (s1->GetClusterIndex2(i))==(s2->GetClusterIndex2(i)) && s1->GetClusterIndex2(i)>0) {
2114 s1->SetSharedMapBit(i, kTRUE);
2115 s2->SetSharedMapBit(i, kTRUE);
2117 if (s1->GetClusterIndex2(i)>0)
2118 s1->SetClusterMapBit(i, kTRUE);
2120 if (sumshared>cutN0){
2123 for (Int_t i=firstpoint;i<lastpoint;i++){
2124 // if ( (s1->GetClusterIndex2(i)&0xFFFF8FFF)==(s2->GetClusterIndex2(i)&0xFFFF8FFF) && s1->GetClusterIndex2(i)>0) {
2125 if ( (s1->GetClusterIndex2(i))==(s2->GetClusterIndex2(i)) && s1->GetClusterIndex2(i)>0) {
2126 AliTPCTrackerPoint *p1 = s1->GetTrackPoint(i);
2127 AliTPCTrackerPoint *p2 = s2->GetTrackPoint(i);;
2128 if (s1->IsActive()&&s2->IsActive()){
2129 p1->SetShared(kTRUE);
2130 p2->SetShared(kTRUE);
2136 if (sumshared>cutN0){
2137 for (Int_t i=0;i<4;i++){
2138 if (s1->GetOverlapLabel(3*i)==0){
2139 s1->SetOverlapLabel(3*i, s2->GetLabel());
2140 s1->SetOverlapLabel(3*i+1,sumshared);
2141 s1->SetOverlapLabel(3*i+2,s2->GetUniqueID());
2145 for (Int_t i=0;i<4;i++){
2146 if (s2->GetOverlapLabel(3*i)==0){
2147 s2->SetOverlapLabel(3*i, s1->GetLabel());
2148 s2->SetOverlapLabel(3*i+1,sumshared);
2149 s2->SetOverlapLabel(3*i+2,s1->GetUniqueID());
2156 void AliTPCtrackerMI::SignShared(TObjArray * arr)
2159 //sort trackss according sectors
2161 for (Int_t i=0; i<arr->GetEntriesFast(); i++) {
2162 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2164 //if (pt) RotateToLocal(pt);
2168 arr->Sort(); // sorting according relative sectors
2169 arr->Expand(arr->GetEntries());
2172 Int_t nseed=arr->GetEntriesFast();
2173 for (Int_t i=0; i<nseed; i++) {
2174 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2176 for (Int_t j=0;j<12;j++){
2177 pt->SetOverlapLabel(j,0);
2180 for (Int_t i=0; i<nseed; i++) {
2181 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2183 if (pt->GetRemoval()>10) continue;
2184 for (Int_t j=i+1; j<nseed; j++){
2185 AliTPCseed *pt2=(AliTPCseed*)arr->UncheckedAt(j);
2186 if (TMath::Abs(pt->GetRelativeSector()-pt2->GetRelativeSector())>1) continue;
2188 if (pt2->GetRemoval()<=10) {
2189 //if ( TMath::Abs(pt->GetRelativeSector()-pt2->GetRelativeSector())>0) break;
2197 void AliTPCtrackerMI::SortTracks(TObjArray * arr, Int_t mode) const
2200 //sort tracks in array according mode criteria
2201 Int_t nseed = arr->GetEntriesFast();
2202 for (Int_t i=0; i<nseed; i++) {
2203 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2214 void AliTPCtrackerMI::RemoveUsed2(TObjArray * arr, Float_t factor1, Float_t factor2, Int_t minimal)
2217 // Loop over all tracks and remove overlaped tracks (with lower quality)
2219 // 1. Unsign clusters
2220 // 2. Sort tracks according quality
2221 // Quality is defined by the number of cluster between first and last points
2223 // 3. Loop over tracks - decreasing quality order
2224 // a.) remove - If the fraction of shared cluster less than factor (1- n or 2)
2225 // b.) remove - If the minimal number of clusters less than minimal and not ITS
2226 // c.) if track accepted - sign clusters
2228 //Called in - AliTPCtrackerMI::Clusters2Tracks()
2229 // - AliTPCtrackerMI::PropagateBack()
2230 // - AliTPCtrackerMI::RefitInward()
2233 // factor1 - factor for constrained
2234 // factor2 - for non constrained tracks
2235 // if (Float_t(shared+1)/Float_t(found+1)>factor) - DELETE
2239 Int_t nseed = arr->GetEntriesFast();
2240 Float_t * quality = new Float_t[nseed];
2241 Int_t * indexes = new Int_t[nseed];
2245 for (Int_t i=0; i<nseed; i++) {
2246 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2251 pt->UpdatePoints(); //select first last max dens points
2252 Float_t * points = pt->GetPoints();
2253 if (points[3]<0.8) quality[i] =-1;
2254 quality[i] = (points[2]-points[0])+pt->GetNumberOfClusters();
2255 //prefer high momenta tracks if overlaps
2256 quality[i] *= TMath::Sqrt(TMath::Abs(pt->Pt())+0.5);
2258 TMath::Sort(nseed,quality,indexes);
2261 for (Int_t itrack=0; itrack<nseed; itrack++) {
2262 Int_t trackindex = indexes[itrack];
2263 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(trackindex);
2266 if (quality[trackindex]<0){
2267 delete arr->RemoveAt(trackindex);
2272 Int_t first = Int_t(pt->GetPoints()[0]);
2273 Int_t last = Int_t(pt->GetPoints()[2]);
2274 Double_t factor = (pt->GetBConstrain()) ? factor1: factor2;
2276 Int_t found,foundable,shared;
2277 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
2278 // pt->GetClusterStatistic(0,160, found, foundable,shared,kFALSE);
2279 Bool_t itsgold =kFALSE;
2282 if (pt->GetESD()->GetITSclusters(dummy)>4) itsgold= kTRUE;
2286 if (Float_t(shared+1)/Float_t(found+1)>factor){
2287 if (pt->GetKinkIndexes()[0]!=0) continue; //don't remove tracks - part of the kinks
2288 if( AliTPCReconstructor::StreamLevel()>3){
2289 TTreeSRedirector &cstream = *fDebugStreamer;
2290 cstream<<"RemoveUsed"<<
2291 "iter="<<fIteration<<
2295 delete arr->RemoveAt(trackindex);
2298 if (pt->GetNumberOfClusters()<50&&(found-0.5*shared)<minimal){ //remove short tracks
2299 if (pt->GetKinkIndexes()[0]!=0) continue; //don't remove tracks - part of the kinks
2300 if( AliTPCReconstructor::StreamLevel()>3){
2301 TTreeSRedirector &cstream = *fDebugStreamer;
2302 cstream<<"RemoveShort"<<
2303 "iter="<<fIteration<<
2307 delete arr->RemoveAt(trackindex);
2313 //if (sharedfactor>0.4) continue;
2314 if (pt->GetKinkIndexes()[0]>0) continue;
2315 //Remove tracks with undefined properties - seems
2316 if (pt->GetSigmaY2()<kAlmost0) continue; // ? what is the origin ?
2318 for (Int_t i=first; i<last; i++) {
2319 Int_t index=pt->GetClusterIndex2(i);
2320 // if (index<0 || index&0x8000 ) continue;
2321 if (index<0 || index&0x8000 ) continue;
2322 AliTPCclusterMI *c= pt->GetClusterPointer(i);
2329 Info("RemoveUsed","\n*****\nNumber of good tracks after shared removal\t%d\n",fNtracks);
2335 void AliTPCtrackerMI::DumpClusters(Int_t iter, TObjArray *trackArray)
2338 // Dump clusters after reco
2339 // signed and unsigned cluster can be visualized
2340 // 1. Unsign all cluster
2341 // 2. Sign all used clusters
2344 Int_t nseed = trackArray->GetEntries();
2345 for (Int_t i=0; i<nseed; i++){
2346 AliTPCseed *pt=(AliTPCseed*)trackArray->UncheckedAt(i);
2350 Bool_t isKink=pt->GetKinkIndex(0)!=0;
2351 for (Int_t j=0; j<160; ++j) {
2352 Int_t index=pt->GetClusterIndex2(j);
2353 if (index<0) continue;
2354 AliTPCclusterMI *c= pt->GetClusterPointer(j);
2356 if (isKink) c->Use(100); // kink
2357 c->Use(10); // by default usage 10
2362 for (Int_t sec=0;sec<fkNIS;sec++){
2363 for (Int_t row=0;row<fInnerSec->GetNRows();row++){
2364 AliTPCclusterMI *cl = fInnerSec[sec][row].GetClusters1();
2365 for (Int_t icl =0;icl< fInnerSec[sec][row].GetN1();icl++){
2366 Float_t gx[3]; cl[icl].GetGlobalXYZ(gx);
2367 (*fDebugStreamer)<<"clDump"<<
2375 cl = fInnerSec[sec][row].GetClusters2();
2376 for (Int_t icl =0;icl< fInnerSec[sec][row].GetN2();icl++){
2377 Float_t gx[3]; cl[icl].GetGlobalXYZ(gx);
2378 (*fDebugStreamer)<<"clDump"<<
2389 for (Int_t sec=0;sec<fkNOS;sec++){
2390 for (Int_t row=0;row<fOuterSec->GetNRows();row++){
2391 AliTPCclusterMI *cl = fOuterSec[sec][row].GetClusters1();
2392 for (Int_t icl =0;icl< fOuterSec[sec][row].GetN1();icl++){
2393 Float_t gx[3]; cl[icl].GetGlobalXYZ(gx);
2394 (*fDebugStreamer)<<"clDump"<<
2402 cl = fOuterSec[sec][row].GetClusters2();
2403 for (Int_t icl =0;icl< fOuterSec[sec][row].GetN2();icl++){
2404 Float_t gx[3]; cl[icl].GetGlobalXYZ(gx);
2405 (*fDebugStreamer)<<"clDump"<<
2417 void AliTPCtrackerMI::UnsignClusters()
2420 // loop over all clusters and unsign them
2423 for (Int_t sec=0;sec<fkNIS;sec++){
2424 for (Int_t row=0;row<fInnerSec->GetNRows();row++){
2425 AliTPCclusterMI *cl = fInnerSec[sec][row].GetClusters1();
2426 for (Int_t icl =0;icl< fInnerSec[sec][row].GetN1();icl++)
2427 // if (cl[icl].IsUsed(10))
2429 cl = fInnerSec[sec][row].GetClusters2();
2430 for (Int_t icl =0;icl< fInnerSec[sec][row].GetN2();icl++)
2431 //if (cl[icl].IsUsed(10))
2436 for (Int_t sec=0;sec<fkNOS;sec++){
2437 for (Int_t row=0;row<fOuterSec->GetNRows();row++){
2438 AliTPCclusterMI *cl = fOuterSec[sec][row].GetClusters1();
2439 for (Int_t icl =0;icl< fOuterSec[sec][row].GetN1();icl++)
2440 //if (cl[icl].IsUsed(10))
2442 cl = fOuterSec[sec][row].GetClusters2();
2443 for (Int_t icl =0;icl< fOuterSec[sec][row].GetN2();icl++)
2444 //if (cl[icl].IsUsed(10))
2453 void AliTPCtrackerMI::SignClusters(const TObjArray * arr, Float_t fnumber, Float_t fdensity)
2456 //sign clusters to be "used"
2458 // snumber and sdensity sign number of sigmas - bellow mean value to be accepted
2459 // loop over "primaries"
2473 Int_t nseed = arr->GetEntriesFast();
2474 for (Int_t i=0; i<nseed; i++) {
2475 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2479 if (!(pt->IsActive())) continue;
2480 Float_t dens = pt->GetNumberOfClusters()/Float_t(pt->GetNFoundable());
2481 if ( (dens>0.7) && (pt->GetNumberOfClusters()>70)){
2483 sumdens2+= dens*dens;
2484 sumn += pt->GetNumberOfClusters();
2485 sumn2 += pt->GetNumberOfClusters()*pt->GetNumberOfClusters();
2486 Float_t chi2 = pt->GetChi2()/pt->GetNumberOfClusters();
2489 sumchi2 +=chi2*chi2;
2494 Float_t mdensity = 0.9;
2495 Float_t meann = 130;
2496 Float_t meanchi = 1;
2497 Float_t sdensity = 0.1;
2498 Float_t smeann = 10;
2499 Float_t smeanchi =0.4;
2503 mdensity = sumdens/sum;
2505 meanchi = sumchi/sum;
2507 sdensity = sumdens2/sum-mdensity*mdensity;
2509 sdensity = TMath::Sqrt(sdensity);
2513 smeann = sumn2/sum-meann*meann;
2515 smeann = TMath::Sqrt(smeann);
2519 smeanchi = sumchi2/sum - meanchi*meanchi;
2521 smeanchi = TMath::Sqrt(smeanchi);
2527 //REMOVE SHORT DELTAS or tracks going out of sensitive volume of TPC
2529 for (Int_t i=0; i<nseed; i++) {
2530 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2534 if (pt->GetBSigned()) continue;
2535 if (pt->GetBConstrain()) continue;
2536 //if (!(pt->IsActive())) continue;
2538 Int_t found,foundable,shared;
2539 pt->GetClusterStatistic(0,160,found, foundable,shared);
2540 if (shared/float(found)>0.3) {
2541 if (shared/float(found)>0.9 ){
2542 //delete arr->RemoveAt(i);
2547 Bool_t isok =kFALSE;
2548 if ( (pt->GetNShared()/pt->GetNumberOfClusters()<0.5) &&pt->GetNumberOfClusters()>60)
2550 if ((TMath::Abs(1/pt->GetC())<100.) && (pt->GetNShared()/pt->GetNumberOfClusters()<0.7))
2552 if (TMath::Abs(pt->GetZ()/pt->GetX())>1.1)
2554 if ( (TMath::Abs(pt->GetSnp()>0.7) && pt->GetD(0,0)>60.))
2558 for (Int_t j=0; j<160; ++j) {
2559 Int_t index=pt->GetClusterIndex2(j);
2560 if (index<0) continue;
2561 AliTPCclusterMI *c= pt->GetClusterPointer(j);
2563 //if (!(c->IsUsed(10))) c->Use();
2570 Double_t maxchi = meanchi+2.*smeanchi;
2572 for (Int_t i=0; i<nseed; i++) {
2573 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2577 //if (!(pt->IsActive())) continue;
2578 if (pt->GetBSigned()) continue;
2579 Double_t chi = pt->GetChi2()/pt->GetNumberOfClusters();
2580 if (chi>maxchi) continue;
2583 Float_t dens = pt->GetNumberOfClusters()/Float_t(pt->GetNFoundable());
2585 //sign only tracks with enoug big density at the beginning
2587 if ((pt->GetDensityFirst(40)<0.75) && pt->GetNumberOfClusters()<meann) continue;
2590 Double_t mindens = TMath::Max(double(mdensity-sdensity*fdensity*bfactor),0.65);
2591 Double_t minn = TMath::Max(Int_t(meann-fnumber*smeann*bfactor),50);
2593 // if (pt->fBConstrain) mindens = TMath::Max(mdensity-sdensity*fdensity*bfactor,0.65);
2594 if ( (pt->GetRemoval()==10) && (pt->GetSnp()>0.8)&&(dens>mindens))
2597 if ((dens>mindens && pt->GetNumberOfClusters()>minn) && chi<maxchi ){
2598 //Int_t noc=pt->GetNumberOfClusters();
2599 pt->SetBSigned(kTRUE);
2600 for (Int_t j=0; j<160; ++j) {
2602 Int_t index=pt->GetClusterIndex2(j);
2603 if (index<0) continue;
2604 AliTPCclusterMI *c= pt->GetClusterPointer(j);
2606 // if (!(c->IsUsed(10))) c->Use();
2611 // gLastCheck = nseed;
2619 void AliTPCtrackerMI::StopNotActive(const TObjArray * arr, Int_t row0, Float_t th0, Float_t th1, Float_t th2) const
2621 // stop not active tracks
2622 // take th1 as threshold for number of founded to number of foundable on last 10 active rows
2623 // take th2 as threshold for number of founded to number of foundable on last 20 active rows
2624 Int_t nseed = arr->GetEntriesFast();
2626 for (Int_t i=0; i<nseed; i++) {
2627 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2631 if (!(pt->IsActive())) continue;
2632 StopNotActive(pt,row0,th0, th1,th2);
2638 void AliTPCtrackerMI::StopNotActive(AliTPCseed * seed, Int_t row0, Float_t th0, Float_t th1,
2641 // stop not active tracks
2642 // take th1 as threshold for number of founded to number of foundable on last 10 active rows
2643 // take th2 as threshold for number of founded to number of foundable on last 20 active rows
2646 Int_t foundable = 0;
2647 Int_t maxindex = seed->GetLastPoint(); //last foundable row
2648 if (seed->GetNFoundable()*th0 > seed->GetNumberOfClusters()) {
2649 seed->Desactivate(10) ;
2653 for (Int_t i=row0; i<maxindex; i++){
2654 Int_t index = seed->GetClusterIndex2(i);
2655 if (index!=-1) foundable++;
2657 if (foundable<=30) sumgood1++;
2658 if (foundable<=50) {
2665 if (foundable>=30.){
2666 if (sumgood1<(th1*30.)) seed->Desactivate(10);
2669 if (sumgood2<(th2*50.)) seed->Desactivate(10);
2673 Int_t AliTPCtrackerMI::RefitInward(AliESDEvent *event)
2676 // back propagation of ESD tracks
2679 if (!event) return 0;
2680 const Int_t kMaxFriendTracks=2000;
2682 // extract correction object for multiplicity dependence of dEdx
2683 TObjArray * gainCalibArray = AliTPCcalibDB::Instance()->GetTimeGainSplinesRun(event->GetRunNumber());
2685 AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ;
2687 AliFatal("Tranformations not in RefitInward");
2690 transform->SetCurrentRecoParam((AliTPCRecoParam*)AliTPCReconstructor::GetRecoParam());
2691 const AliTPCRecoParam * recoParam = AliTPCcalibDB::Instance()->GetTransform()->GetCurrentRecoParam();
2692 Int_t nContribut = event->GetNumberOfTracks();
2693 TGraphErrors * graphMultDependenceDeDx = 0x0;
2694 if (recoParam && recoParam->GetUseMultiplicityCorrectionDedx() && gainCalibArray) {
2695 if (recoParam->GetUseTotCharge()) {
2696 graphMultDependenceDeDx = (TGraphErrors *) gainCalibArray->FindObject("TGRAPHERRORS_MEANQTOT_MULTIPLICITYDEPENDENCE_BEAM_ALL");
2698 graphMultDependenceDeDx = (TGraphErrors *) gainCalibArray->FindObject("TGRAPHERRORS_MEANQMAX_MULTIPLICITYDEPENDENCE_BEAM_ALL");
2704 //PrepareForProlongation(fSeeds,1);
2705 PropagateForward2(fSeeds);
2706 RemoveUsed2(fSeeds,0.4,0.4,20);
2708 TObjArray arraySeed(fSeeds->GetEntries());
2709 for (Int_t i=0;i<fSeeds->GetEntries();i++) {
2710 arraySeed.AddAt(fSeeds->At(i),i);
2712 SignShared(&arraySeed);
2713 // FindCurling(fSeeds, event,2); // find multi found tracks
2714 FindSplitted(fSeeds, event,2); // find multi found tracks
2715 if (AliTPCReconstructor::StreamLevel()>5) FindMultiMC(fSeeds, fEvent,2); // find multi found tracks
2718 Int_t nseed = fSeeds->GetEntriesFast();
2719 for (Int_t i=0;i<nseed;i++){
2720 AliTPCseed * seed = (AliTPCseed*) fSeeds->UncheckedAt(i);
2721 if (!seed) continue;
2722 if (seed->GetKinkIndex(0)>0) UpdateKinkQualityD(seed); // update quality informations for kinks
2723 AliESDtrack *esd=event->GetTrack(i);
2725 if (seed->GetNumberOfClusters()<60 && seed->GetNumberOfClusters()<(esd->GetTPCclusters(0) -5)*0.8){
2726 AliExternalTrackParam paramIn;
2727 AliExternalTrackParam paramOut;
2728 Int_t ncl = seed->RefitTrack(seed,¶mIn,¶mOut);
2729 if (AliTPCReconstructor::StreamLevel()>2) {
2730 (*fDebugStreamer)<<"RecoverIn"<<
2734 "pout.="<<¶mOut<<
2739 seed->Set(paramIn.GetX(),paramIn.GetAlpha(),paramIn.GetParameter(),paramIn.GetCovariance());
2740 seed->SetNumberOfClusters(ncl);
2744 seed->PropagateTo(fkParam->GetInnerRadiusLow());
2745 seed->UpdatePoints();
2746 AddCovariance(seed);
2748 seed->CookdEdx(0.02,0.6);
2749 CookLabel(seed,0.1); //For comparison only
2750 esd->SetTPCClusterMap(seed->GetClusterMap());
2751 esd->SetTPCSharedMap(seed->GetSharedMap());
2753 if (AliTPCReconstructor::StreamLevel()>1 && seed!=0) {
2754 TTreeSRedirector &cstream = *fDebugStreamer;
2761 if (seed->GetNumberOfClusters()>15){
2762 esd->UpdateTrackParams(seed,AliESDtrack::kTPCrefit);
2763 esd->SetTPCPoints(seed->GetPoints());
2764 esd->SetTPCPointsF(seed->GetNFoundable());
2765 Int_t ndedx = seed->GetNCDEDX(0);
2766 Float_t sdedx = seed->GetSDEDX(0);
2767 Float_t dedx = seed->GetdEdx();
2768 // apply mutliplicity dependent dEdx correction if available
2769 if (graphMultDependenceDeDx) {
2770 Double_t corrGain = AliTPCcalibDButil::EvalGraphConst(graphMultDependenceDeDx, nContribut);
2771 dedx += (1 - corrGain)*50.; // MIP is normalized to 50
2773 esd->SetTPCsignal(dedx, sdedx, ndedx);
2775 // add seed to the esd track in Calib level
2777 Bool_t storeFriend = gRandom->Rndm()<(kMaxFriendTracks)/Float_t(nseed);
2778 if (AliTPCReconstructor::StreamLevel()>0 &&storeFriend){
2779 AliTPCseed * seedCopy = new AliTPCseed(*seed, kTRUE);
2780 esd->AddCalibObject(seedCopy);
2785 //printf("problem\n");
2788 //FindKinks(fSeeds,event);
2789 if (AliTPCReconstructor::StreamLevel()>3) DumpClusters(2,fSeeds);
2790 Info("RefitInward","Number of refitted tracks %d",ntracks);
2795 Int_t AliTPCtrackerMI::PropagateBack(AliESDEvent *event)
2798 // back propagation of ESD tracks
2800 if (!event) return 0;
2804 PropagateBack(fSeeds);
2805 RemoveUsed2(fSeeds,0.4,0.4,20);
2806 //FindCurling(fSeeds, fEvent,1);
2807 FindSplitted(fSeeds, event,1); // find multi found tracks
2808 if (AliTPCReconstructor::StreamLevel()>5) FindMultiMC(fSeeds, fEvent,1); // find multi found tracks
2811 Int_t nseed = fSeeds->GetEntriesFast();
2813 for (Int_t i=0;i<nseed;i++){
2814 AliTPCseed * seed = (AliTPCseed*) fSeeds->UncheckedAt(i);
2815 if (!seed) continue;
2816 if (seed->GetKinkIndex(0)<0) UpdateKinkQualityM(seed); // update quality informations for kinks
2817 seed->UpdatePoints();
2818 AddCovariance(seed);
2819 AliESDtrack *esd=event->GetTrack(i);
2820 if (!esd) continue; //never happen
2821 if (seed->GetNumberOfClusters()<60 && seed->GetNumberOfClusters()<(esd->GetTPCclusters(0) -5)*0.8){
2822 AliExternalTrackParam paramIn;
2823 AliExternalTrackParam paramOut;
2824 Int_t ncl = seed->RefitTrack(seed,¶mIn,¶mOut);
2825 if (AliTPCReconstructor::StreamLevel()>2) {
2826 (*fDebugStreamer)<<"RecoverBack"<<
2830 "pout.="<<¶mOut<<
2835 seed->Set(paramOut.GetX(),paramOut.GetAlpha(),paramOut.GetParameter(),paramOut.GetCovariance());
2836 seed->SetNumberOfClusters(ncl);
2839 seed->CookdEdx(0.02,0.6);
2840 CookLabel(seed,0.1); //For comparison only
2841 if (seed->GetNumberOfClusters()>15){
2842 esd->UpdateTrackParams(seed,AliESDtrack::kTPCout);
2843 esd->SetTPCPoints(seed->GetPoints());
2844 esd->SetTPCPointsF(seed->GetNFoundable());
2845 Int_t ndedx = seed->GetNCDEDX(0);
2846 Float_t sdedx = seed->GetSDEDX(0);
2847 Float_t dedx = seed->GetdEdx();
2848 esd->SetTPCsignal(dedx, sdedx, ndedx);
2850 Int_t eventnumber = event->GetEventNumberInFile();// patch 28 fev 06
2851 // This is most likely NOT the event number you'd like to use. It has nothing to do with the 'real' event number
2852 if (AliTPCReconstructor::StreamLevel()>1 && esd) {
2853 (*fDebugStreamer)<<"Cback"<<
2856 "EventNrInFile="<<eventnumber<<
2861 if (AliTPCReconstructor::StreamLevel()>3) DumpClusters(1,fSeeds);
2862 //FindKinks(fSeeds,event);
2863 Info("PropagateBack","Number of back propagated tracks %d",ntracks);
2870 void AliTPCtrackerMI::DeleteSeeds()
2875 Int_t nseed = fSeeds->GetEntriesFast();
2876 for (Int_t i=0;i<nseed;i++){
2877 AliTPCseed * seed = (AliTPCseed*)fSeeds->At(i);
2878 if (seed) delete fSeeds->RemoveAt(i);
2885 void AliTPCtrackerMI::ReadSeeds(const AliESDEvent *const event, Int_t direction)
2888 //read seeds from the event
2890 Int_t nentr=event->GetNumberOfTracks();
2892 Info("ReadSeeds", "Number of ESD tracks: %d\n", nentr);
2897 fSeeds = new TObjArray(nentr);
2901 for (Int_t i=0; i<nentr; i++) {
2902 AliESDtrack *esd=event->GetTrack(i);
2903 ULong_t status=esd->GetStatus();
2904 if (!(status&AliESDtrack::kTPCin)) continue;
2905 AliTPCtrack t(*esd);
2906 t.SetNumberOfClusters(0);
2907 // AliTPCseed *seed = new AliTPCseed(t,t.GetAlpha());
2908 AliTPCseed *seed = new AliTPCseed(t/*,t.GetAlpha()*/);
2909 seed->SetUniqueID(esd->GetID());
2910 AddCovariance(seed); //add systematic ucertainty
2911 for (Int_t ikink=0;ikink<3;ikink++) {
2912 Int_t index = esd->GetKinkIndex(ikink);
2913 seed->GetKinkIndexes()[ikink] = index;
2914 if (index==0) continue;
2915 index = TMath::Abs(index);
2916 AliESDkink * kink = fEvent->GetKink(index-1);
2917 if (kink&&esd->GetKinkIndex(ikink)<0){
2918 if ((status & AliESDtrack::kTRDrefit) != 0) kink->SetStatus(1,2);
2919 if ((status & AliESDtrack::kITSout) != 0) kink->SetStatus(1,0);
2921 if (kink&&esd->GetKinkIndex(ikink)>0){
2922 if ((status & AliESDtrack::kTRDrefit) != 0) kink->SetStatus(1,6);
2923 if ((status & AliESDtrack::kITSout) != 0) kink->SetStatus(1,4);
2927 if (((status&AliESDtrack::kITSout)==0)&&(direction==1)) seed->ResetCovariance(10.);
2928 if ( direction ==2 &&(status & AliESDtrack::kTRDrefit) == 0 ) seed->ResetCovariance(10.);
2929 //if ( direction ==2 && ((status & AliESDtrack::kTPCout) == 0) ) {
2930 // fSeeds->AddAt(0,i);
2934 if ( direction ==2 &&(status & AliESDtrack::kTRDrefit) > 0 ) {
2935 Double_t par0[5],par1[5],alpha,x;
2936 esd->GetInnerExternalParameters(alpha,x,par0);
2937 esd->GetExternalParameters(x,par1);
2938 Double_t delta1 = TMath::Abs(par0[4]-par1[4])/(0.000000001+TMath::Abs(par0[4]+par1[4]));
2939 Double_t delta2 = TMath::Abs(par0[3]-par1[3]);
2941 if (esd->GetTRDncls()>0) trdchi2 = esd->GetTRDchi2()/esd->GetTRDncls();
2942 //reset covariance if suspicious
2943 if ( (delta1>0.1) || (delta2>0.006) ||trdchi2>7.)
2944 seed->ResetCovariance(10.);
2949 // rotate to the local coordinate system
2951 fSectors=fInnerSec; fN=fkNIS;
2952 Double_t alpha=seed->GetAlpha() - fSectors->GetAlphaShift();
2953 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
2954 if (alpha < 0. ) alpha += 2.*TMath::Pi();
2955 Int_t ns=Int_t(alpha/fSectors->GetAlpha())%fN;
2956 alpha =ns*fSectors->GetAlpha() + fSectors->GetAlphaShift();
2957 if (alpha<-TMath::Pi()) alpha += 2*TMath::Pi();
2958 if (alpha>=TMath::Pi()) alpha -= 2*TMath::Pi();
2959 alpha-=seed->GetAlpha();
2960 if (!seed->Rotate(alpha)) {
2966 if (esd->GetKinkIndex(0)<=0){
2967 for (Int_t irow=0;irow<160;irow++){
2968 Int_t index = seed->GetClusterIndex2(irow);
2971 AliTPCclusterMI * cl = GetClusterMI(index);
2972 seed->SetClusterPointer(irow,cl);
2974 if ((index & 0x8000)==0){
2975 cl->Use(10); // accepted cluster
2977 cl->Use(6); // close cluster not accepted
2980 Info("ReadSeeds","Not found cluster");
2985 fSeeds->AddAt(seed,i);
2991 //_____________________________________________________________________________
2992 void AliTPCtrackerMI::MakeSeeds3(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, Float_t cuts[4],
2993 Float_t deltay, Int_t ddsec) {
2994 //-----------------------------------------------------------------
2995 // This function creates track seeds.
2996 // SEEDING WITH VERTEX CONSTRAIN
2997 //-----------------------------------------------------------------
2998 // cuts[0] - fP4 cut
2999 // cuts[1] - tan(phi) cut
3000 // cuts[2] - zvertex cut
3001 // cuts[3] - fP3 cut
3009 Double_t x[5], c[15];
3010 // Int_t di = i1-i2;
3012 AliTPCseed * seed = new AliTPCseed();
3013 Double_t alpha=fSectors->GetAlpha(), shift=fSectors->GetAlphaShift();
3014 Double_t cs=cos(alpha), sn=sin(alpha);
3016 // Double_t x1 =fOuterSec->GetX(i1);
3017 //Double_t xx2=fOuterSec->GetX(i2);
3019 Double_t x1 =GetXrow(i1);
3020 Double_t xx2=GetXrow(i2);
3022 Double_t x3=GetX(), y3=GetY(), z3=GetZ();
3024 Int_t imiddle = (i2+i1)/2; //middle pad row index
3025 Double_t xm = GetXrow(imiddle); // radius of middle pad-row
3026 const AliTPCtrackerRow& krm=GetRow(sec,imiddle); //middle pad -row
3030 const AliTPCtrackerRow& kr1=GetRow(ns,i1);
3031 Double_t ymax = GetMaxY(i1)-kr1.GetDeadZone()-1.5;
3032 Double_t ymaxm = GetMaxY(imiddle)-kr1.GetDeadZone()-1.5;
3035 // change cut on curvature if it can't reach this layer
3036 // maximal curvature set to reach it
3037 Double_t dvertexmax = TMath::Sqrt((x1-x3)*(x1-x3)+(ymax+5-y3)*(ymax+5-y3));
3038 if (dvertexmax*0.5*cuts[0]>0.85){
3039 cuts[0] = 0.85/(dvertexmax*0.5+1.);
3041 Double_t r2min = 1/(cuts[0]*cuts[0]); //minimal square of radius given by cut
3044 if (deltay>0) ddsec = 0;
3045 // loop over clusters
3046 for (Int_t is=0; is < kr1; is++) {
3048 if (kr1[is]->IsUsed(10)) continue;
3049 Double_t y1=kr1[is]->GetY(), z1=kr1[is]->GetZ();
3050 //if (TMath::Abs(y1)>ymax) continue;
3052 if (deltay>0 && TMath::Abs(ymax-TMath::Abs(y1))> deltay ) continue; // seed only at the edge
3054 // find possible directions
3055 Float_t anglez = (z1-z3)/(x1-x3);
3056 Float_t extraz = z1 - anglez*(x1-xx2); // extrapolated z
3059 //find rotation angles relative to line given by vertex and point 1
3060 Double_t dvertex2 = (x1-x3)*(x1-x3)+(y1-y3)*(y1-y3);
3061 Double_t dvertex = TMath::Sqrt(dvertex2);
3062 Double_t angle13 = TMath::ATan((y1-y3)/(x1-x3));
3063 Double_t cs13 = cos(-angle13), sn13 = sin(-angle13);
3066 // loop over 2 sectors
3072 Double_t dddz1=0; // direction of delta inclination in z axis
3079 for (Int_t dsec = dsec1; dsec<=dsec2;dsec++){
3080 Int_t sec2 = sec + dsec;
3082 // AliTPCtrackerRow& kr2 = fOuterSec[(sec2+fkNOS)%fkNOS][i2];
3083 //AliTPCtrackerRow& kr2m = fOuterSec[(sec2+fkNOS)%fkNOS][imiddle];
3084 AliTPCtrackerRow& kr2 = GetRow((sec2+fkNOS)%fkNOS,i2);
3085 AliTPCtrackerRow& kr2m = GetRow((sec2+fkNOS)%fkNOS,imiddle);
3086 Int_t index1 = TMath::Max(kr2.Find(extraz-0.6-dddz1*TMath::Abs(z1)*0.05)-1,0);
3087 Int_t index2 = TMath::Min(kr2.Find(extraz+0.6+dddz2*TMath::Abs(z1)*0.05)+1,kr2);
3089 // rotation angles to p1-p3
3090 Double_t cs13r = cos(-angle13+dsec*alpha)/dvertex, sn13r = sin(-angle13+dsec*alpha)/dvertex;
3091 Double_t x2, y2, z2;
3093 // Double_t dymax = maxangle*TMath::Abs(x1-xx2);
3096 Double_t dxx0 = (xx2-x3)*cs13r;
3097 Double_t dyy0 = (xx2-x3)*sn13r;
3098 for (Int_t js=index1; js < index2; js++) {
3099 const AliTPCclusterMI *kcl = kr2[js];
3100 if (kcl->IsUsed(10)) continue;
3102 //calcutate parameters
3104 Double_t yy0 = dyy0 +(kcl->GetY()-y3)*cs13r;
3106 if (TMath::Abs(yy0)<0.000001) continue;
3107 Double_t xx0 = dxx0 -(kcl->GetY()-y3)*sn13r;
3108 Double_t y0 = 0.5*(xx0*xx0+yy0*yy0-xx0)/yy0;
3109 Double_t r02 = (0.25+y0*y0)*dvertex2;
3110 //curvature (radius) cut
3111 if (r02<r2min) continue;
3115 Double_t c0 = 1/TMath::Sqrt(r02);
3119 //Double_t dfi0 = 2.*TMath::ASin(dvertex*c0*0.5);
3120 //Double_t dfi1 = 2.*TMath::ASin(TMath::Sqrt(yy0*yy0+(1-xx0)*(1-xx0))*dvertex*c0*0.5);
3121 Double_t dfi0 = 2.*AliTPCFastMath::FastAsin(dvertex*c0*0.5);
3122 Double_t dfi1 = 2.*AliTPCFastMath::FastAsin(TMath::Sqrt(yy0*yy0+(1-xx0)*(1-xx0))*dvertex*c0*0.5);
3125 Double_t z0 = kcl->GetZ();
3126 Double_t zzzz2 = z1-(z1-z3)*dfi1/dfi0;
3127 if (TMath::Abs(zzzz2-z0)>0.5) continue;
3130 Double_t dip = (z1-z0)*c0/dfi1;
3131 Double_t x0 = (0.5*cs13+y0*sn13)*dvertex*c0;
3142 x2= xx2*cs-y2*sn*dsec;
3143 y2=+xx2*sn*dsec+y2*cs;
3153 // do we have cluster at the middle ?
3155 GetProlongation(x1,xm,x,ym,zm);
3157 AliTPCclusterMI * cm=0;
3158 if (TMath::Abs(ym)-ymaxm<0){
3159 cm = krm.FindNearest2(ym,zm,1.0,0.6,dummy);
3160 if ((!cm) || (cm->IsUsed(10))) {
3165 // rotate y1 to system 0
3166 // get state vector in rotated system
3167 Double_t yr1 = (-0.5*sn13+y0*cs13)*dvertex*c0;
3168 Double_t xr2 = x0*cs+yr1*sn*dsec;
3169 Double_t xr[5]={kcl->GetY(),kcl->GetZ(), xr2, dip, c0};
3171 GetProlongation(xx2,xm,xr,ym,zm);
3172 if (TMath::Abs(ym)-ymaxm<0){
3173 cm = kr2m.FindNearest2(ym,zm,1.0,0.6,dummy);
3174 if ((!cm) || (cm->IsUsed(10))) {
3184 dym = ym - cm->GetY();
3185 dzm = zm - cm->GetZ();
3192 Double_t sy1=kr1[is]->GetSigmaY2()*2., sz1=kr1[is]->GetSigmaZ2()*2.;
3193 Double_t sy2=kcl->GetSigmaY2()*2., sz2=kcl->GetSigmaZ2()*2.;
3194 //Double_t sy3=400*3./12., sy=0.1, sz=0.1;
3195 Double_t sy3=25000*x[4]*x[4]+0.1, sy=0.1, sz=0.1;
3196 //Double_t sy3=25000*x[4]*x[4]*60+0.5, sy=0.1, sz=0.1;
3198 Double_t f40=(F1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
3199 Double_t f42=(F1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
3200 Double_t f43=(F1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
3201 Double_t f20=(F2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
3202 Double_t f22=(F2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
3203 Double_t f23=(F2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
3205 Double_t f30=(F3(x1,y1+sy,x2,y2,z1,z2)-x[3])/sy;
3206 Double_t f31=(F3(x1,y1,x2,y2,z1+sz,z2)-x[3])/sz;
3207 Double_t f32=(F3(x1,y1,x2,y2+sy,z1,z2)-x[3])/sy;
3208 Double_t f34=(F3(x1,y1,x2,y2,z1,z2+sz)-x[3])/sz;
3212 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3213 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3214 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3215 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3216 c[13]=f30*sy1*f40+f32*sy2*f42;
3217 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3219 // if (!BuildSeed(kr1[is],kcl,0,x1,x2,x3,x,c)) continue;
3221 UInt_t index=kr1.GetIndex(is);
3222 seed->~AliTPCseed(); // this does not set the pointer to 0...
3223 AliTPCseed *track=new(seed) AliTPCseed(x1, ns*alpha+shift, x, c, index);
3225 track->SetIsSeeding(kTRUE);
3226 track->SetSeed1(i1);
3227 track->SetSeed2(i2);
3228 track->SetSeedType(3);
3232 FollowProlongation(*track, (i1+i2)/2,1);
3233 Int_t foundable,found,shared;
3234 track->GetClusterStatistic((i1+i2)/2,i1, found, foundable, shared, kTRUE);
3235 if ((found<0.55*foundable) || shared>0.5*found || (track->GetSigmaY2()+track->GetSigmaZ2())>0.5){
3237 seed->~AliTPCseed();
3243 FollowProlongation(*track, i2,1);
3247 track->SetBConstrain(1);
3248 // track->fLastPoint = i1+fInnerSec->GetNRows(); // first cluster in track position
3249 track->SetLastPoint(i1); // first cluster in track position
3250 track->SetFirstPoint(track->GetLastPoint());
3252 if (track->GetNumberOfClusters()<(i1-i2)*0.5 ||
3253 track->GetNumberOfClusters() < track->GetNFoundable()*0.6 ||
3254 track->GetNShared()>0.4*track->GetNumberOfClusters() ) {
3256 seed->~AliTPCseed();
3260 // Z VERTEX CONDITION
3261 Double_t zv, bz=GetBz();
3262 if ( !track->GetZAt(0.,bz,zv) ) continue;
3263 if (TMath::Abs(zv-z3)>cuts[2]) {
3264 FollowProlongation(*track, TMath::Max(i2-20,0));
3265 if ( !track->GetZAt(0.,bz,zv) ) continue;
3266 if (TMath::Abs(zv-z3)>cuts[2]){
3267 FollowProlongation(*track, TMath::Max(i2-40,0));
3268 if ( !track->GetZAt(0.,bz,zv) ) continue;
3269 if (TMath::Abs(zv-z3)>cuts[2] &&(track->GetNumberOfClusters() > track->GetNFoundable()*0.7)){
3270 // make seed without constrain
3271 AliTPCseed * track2 = MakeSeed(track,0.2,0.5,1.);
3272 FollowProlongation(*track2, i2,1);
3273 track2->SetBConstrain(kFALSE);
3274 track2->SetSeedType(1);
3275 arr->AddLast(track2);
3277 seed->~AliTPCseed();
3282 seed->~AliTPCseed();
3289 track->SetSeedType(0);
3290 arr->AddLast(track);
3291 seed = new AliTPCseed;
3293 // don't consider other combinations
3294 if (track->GetNumberOfClusters() > track->GetNFoundable()*0.8)
3300 Info("MakeSeeds3","\nSeeding statistic:\t%d\t%d\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin,nout1,nout2);
3306 void AliTPCtrackerMI::MakeSeeds5(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, Float_t cuts[4],
3311 //-----------------------------------------------------------------
3312 // This function creates track seeds.
3313 //-----------------------------------------------------------------
3314 // cuts[0] - fP4 cut
3315 // cuts[1] - tan(phi) cut
3316 // cuts[2] - zvertex cut
3317 // cuts[3] - fP3 cut
3327 Double_t x[5], c[15];
3329 // make temporary seed
3330 AliTPCseed * seed = new AliTPCseed;
3331 Double_t alpha=fOuterSec->GetAlpha(), shift=fOuterSec->GetAlphaShift();
3332 // Double_t cs=cos(alpha), sn=sin(alpha);
3337 Double_t x1 = GetXrow(i1-1);
3338 const AliTPCtrackerRow& kr1=GetRow(sec,i1-1);
3339 Double_t y1max = GetMaxY(i1-1)-kr1.GetDeadZone()-1.5;
3341 Double_t x1p = GetXrow(i1);
3342 const AliTPCtrackerRow& kr1p=GetRow(sec,i1);
3344 Double_t x1m = GetXrow(i1-2);
3345 const AliTPCtrackerRow& kr1m=GetRow(sec,i1-2);
3348 //last 3 padrow for seeding
3349 AliTPCtrackerRow& kr3 = GetRow((sec+fkNOS)%fkNOS,i1-7);
3350 Double_t x3 = GetXrow(i1-7);
3351 // Double_t y3max= GetMaxY(i1-7)-kr3.fDeadZone-1.5;
3353 AliTPCtrackerRow& kr3p = GetRow((sec+fkNOS)%fkNOS,i1-6);
3354 Double_t x3p = GetXrow(i1-6);
3356 AliTPCtrackerRow& kr3m = GetRow((sec+fkNOS)%fkNOS,i1-8);
3357 Double_t x3m = GetXrow(i1-8);
3362 Int_t im = i1-4; //middle pad row index
3363 Double_t xm = GetXrow(im); // radius of middle pad-row
3364 const AliTPCtrackerRow& krm=GetRow(sec,im); //middle pad -row
3365 // Double_t ymmax = GetMaxY(im)-kr1.fDeadZone-1.5;
3368 Double_t deltax = x1-x3;
3369 Double_t dymax = deltax*cuts[1];
3370 Double_t dzmax = deltax*cuts[3];
3372 // loop over clusters
3373 for (Int_t is=0; is < kr1; is++) {
3375 if (kr1[is]->IsUsed(10)) continue;
3376 Double_t y1=kr1[is]->GetY(), z1=kr1[is]->GetZ();
3378 if (deltay>0 && TMath::Abs(y1max-TMath::Abs(y1))> deltay ) continue; // seed only at the edge
3380 Int_t index1 = TMath::Max(kr3.Find(z1-dzmax)-1,0);
3381 Int_t index2 = TMath::Min(kr3.Find(z1+dzmax)+1,kr3);
3387 for (Int_t js=index1; js < index2; js++) {
3388 const AliTPCclusterMI *kcl = kr3[js];
3389 if (kcl->IsUsed(10)) continue;
3391 // apply angular cuts
3392 if (TMath::Abs(y1-y3)>dymax) continue;
3395 if (TMath::Abs(z1-z3)>dzmax) continue;
3397 Double_t angley = (y1-y3)/(x1-x3);
3398 Double_t anglez = (z1-z3)/(x1-x3);
3400 Double_t erry = TMath::Abs(angley)*(x1-x1m)*0.5+0.5;
3401 Double_t errz = TMath::Abs(anglez)*(x1-x1m)*0.5+0.5;
3403 Double_t yyym = angley*(xm-x1)+y1;
3404 Double_t zzzm = anglez*(xm-x1)+z1;
3406 const AliTPCclusterMI *kcm = krm.FindNearest2(yyym,zzzm,erry,errz,index);
3408 if (kcm->IsUsed(10)) continue;
3410 erry = TMath::Abs(angley)*(x1-x1m)*0.4+0.5;
3411 errz = TMath::Abs(anglez)*(x1-x1m)*0.4+0.5;
3418 // look around first
3419 const AliTPCclusterMI *kc1m = kr1m.FindNearest2(angley*(x1m-x1)+y1,
3425 if (kc1m->IsUsed(10)) used++;
3427 const AliTPCclusterMI *kc1p = kr1p.FindNearest2(angley*(x1p-x1)+y1,
3433 if (kc1p->IsUsed(10)) used++;
3435 if (used>1) continue;
3436 if (found<1) continue;
3440 const AliTPCclusterMI *kc3m = kr3m.FindNearest2(angley*(x3m-x3)+y3,
3446 if (kc3m->IsUsed(10)) used++;
3450 const AliTPCclusterMI *kc3p = kr3p.FindNearest2(angley*(x3p-x3)+y3,
3456 if (kc3p->IsUsed(10)) used++;
3460 if (used>1) continue;
3461 if (found<3) continue;
3471 x[4]=F1(x1,y1,x2,y2,x3,y3);
3472 //if (TMath::Abs(x[4]) >= cuts[0]) continue;
3475 x[2]=F2(x1,y1,x2,y2,x3,y3);
3478 x[3]=F3n(x1,y1,x2,y2,z1,z2,x[4]);
3479 //if (TMath::Abs(x[3]) > cuts[3]) continue;
3483 Double_t sy1=0.1, sz1=0.1;
3484 Double_t sy2=0.1, sz2=0.1;
3485 Double_t sy3=0.1, sy=0.1, sz=0.1;
3487 Double_t f40=(F1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
3488 Double_t f42=(F1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
3489 Double_t f43=(F1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
3490 Double_t f20=(F2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
3491 Double_t f22=(F2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
3492 Double_t f23=(F2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
3494 Double_t f30=(F3(x1,y1+sy,x2,y2,z1,z2)-x[3])/sy;
3495 Double_t f31=(F3(x1,y1,x2,y2,z1+sz,z2)-x[3])/sz;
3496 Double_t f32=(F3(x1,y1,x2,y2+sy,z1,z2)-x[3])/sy;
3497 Double_t f34=(F3(x1,y1,x2,y2,z1,z2+sz)-x[3])/sz;
3501 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3502 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3503 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3504 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3505 c[13]=f30*sy1*f40+f32*sy2*f42;
3506 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3508 // if (!BuildSeed(kr1[is],kcl,0,x1,x2,x3,x,c)) continue;
3510 index=kr1.GetIndex(is);
3511 seed->~AliTPCseed();
3512 AliTPCseed *track=new(seed) AliTPCseed(x1, sec*alpha+shift, x, c, index);
3514 track->SetIsSeeding(kTRUE);
3517 FollowProlongation(*track, i1-7,1);
3518 if (track->GetNumberOfClusters() < track->GetNFoundable()*0.75 ||
3519 track->GetNShared()>0.6*track->GetNumberOfClusters() || ( track->GetSigmaY2()+ track->GetSigmaZ2())>0.6){
3521 seed->~AliTPCseed();
3527 FollowProlongation(*track, i2,1);
3528 track->SetBConstrain(0);
3529 track->SetLastPoint(i1+fInnerSec->GetNRows()); // first cluster in track position
3530 track->SetFirstPoint(track->GetLastPoint());
3532 if (track->GetNumberOfClusters()<(i1-i2)*0.5 ||
3533 track->GetNumberOfClusters()<track->GetNFoundable()*0.7 ||
3534 track->GetNShared()>2. || track->GetChi2()/track->GetNumberOfClusters()>6 || ( track->GetSigmaY2()+ track->GetSigmaZ2())>0.5 ) {
3536 seed->~AliTPCseed();
3541 FollowProlongation(*track, TMath::Max(i2-10,0),1);
3542 AliTPCseed * track2 = MakeSeed(track,0.2,0.5,0.9);
3543 FollowProlongation(*track2, i2,1);
3544 track2->SetBConstrain(kFALSE);
3545 track2->SetSeedType(4);
3546 arr->AddLast(track2);
3548 seed->~AliTPCseed();
3552 //arr->AddLast(track);
3553 //seed = new AliTPCseed;
3559 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);
3565 //_____________________________________________________________________________
3566 void AliTPCtrackerMI::MakeSeeds2(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, Float_t */*cuts[4]*/,
3567 Float_t deltay, Bool_t /*bconstrain*/) {
3568 //-----------------------------------------------------------------
3569 // This function creates track seeds - without vertex constraint
3570 //-----------------------------------------------------------------
3571 // cuts[0] - fP4 cut - not applied
3572 // cuts[1] - tan(phi) cut
3573 // cuts[2] - zvertex cut - not applied
3574 // cuts[3] - fP3 cut
3584 Double_t alpha=fOuterSec->GetAlpha(), shift=fOuterSec->GetAlphaShift();
3585 // Double_t cs=cos(alpha), sn=sin(alpha);
3586 Int_t row0 = (i1+i2)/2;
3587 Int_t drow = (i1-i2)/2;
3588 const AliTPCtrackerRow& kr0=fSectors[sec][row0];
3589 AliTPCtrackerRow * kr=0;
3591 AliTPCpolyTrack polytrack;
3592 Int_t nclusters=fSectors[sec][row0];
3593 AliTPCseed * seed = new AliTPCseed;
3598 for (Int_t is=0; is < nclusters; is++) { //LOOP over clusters
3600 Int_t nfoundable =0;
3601 for (Int_t iter =1; iter<2; iter++){ //iterations
3602 const AliTPCtrackerRow& krm=fSectors[sec][row0-iter];
3603 const AliTPCtrackerRow& krp=fSectors[sec][row0+iter];
3604 const AliTPCclusterMI * cl= kr0[is];
3606 if (cl->IsUsed(10)) {
3612 Double_t x = kr0.GetX();
3613 // Initialization of the polytrack
3618 Double_t y0= cl->GetY();
3619 Double_t z0= cl->GetZ();
3623 Double_t ymax = fSectors->GetMaxY(row0)-kr0.GetDeadZone()-1.5;
3624 if (deltay>0 && TMath::Abs(ymax-TMath::Abs(y0))> deltay ) continue; // seed only at the edge
3626 erry = (0.5)*cl->GetSigmaY2()/TMath::Sqrt(cl->GetQ())*6;
3627 errz = (0.5)*cl->GetSigmaZ2()/TMath::Sqrt(cl->GetQ())*6;
3628 polytrack.AddPoint(x,y0,z0,erry, errz);
3631 if (cl->IsUsed(10)) sumused++;
3634 Float_t roady = (5*TMath::Sqrt(cl->GetSigmaY2()+0.2)+1.)*iter;
3635 Float_t roadz = (5*TMath::Sqrt(cl->GetSigmaZ2()+0.2)+1.)*iter;
3638 AliTPCclusterMI * cl1 = krm.FindNearest(y0,z0,roady,roadz);
3639 if (cl1 && TMath::Abs(ymax-TMath::Abs(y0))) {
3640 erry = (0.5)*cl1->GetSigmaY2()/TMath::Sqrt(cl1->GetQ())*3;
3641 errz = (0.5)*cl1->GetSigmaZ2()/TMath::Sqrt(cl1->GetQ())*3;
3642 if (cl1->IsUsed(10)) sumused++;
3643 polytrack.AddPoint(x,cl1->GetY(),cl1->GetZ(),erry,errz);
3647 AliTPCclusterMI * cl2 = krp.FindNearest(y0,z0,roady,roadz);
3649 erry = (0.5)*cl2->GetSigmaY2()/TMath::Sqrt(cl2->GetQ())*3;
3650 errz = (0.5)*cl2->GetSigmaZ2()/TMath::Sqrt(cl2->GetQ())*3;
3651 if (cl2->IsUsed(10)) sumused++;
3652 polytrack.AddPoint(x,cl2->GetY(),cl2->GetZ(),erry,errz);
3655 if (sumused>0) continue;
3657 polytrack.UpdateParameters();
3663 nfoundable = polytrack.GetN();
3664 nfound = nfoundable;
3666 for (Int_t ddrow = iter+1; ddrow<drow;ddrow++){
3667 Float_t maxdist = 0.8*(1.+3./(ddrow));
3668 for (Int_t delta = -1;delta<=1;delta+=2){
3669 Int_t row = row0+ddrow*delta;
3670 kr = &(fSectors[sec][row]);
3671 Double_t xn = kr->GetX();
3672 Double_t ymax1 = fSectors->GetMaxY(row)-kr->GetDeadZone()-1.5;
3673 polytrack.GetFitPoint(xn,yn,zn);
3674 if (TMath::Abs(yn)>ymax1) continue;
3676 AliTPCclusterMI * cln = kr->FindNearest(yn,zn,roady,roadz);
3678 Float_t dist = TMath::Sqrt( (yn-cln->GetY())*(yn-cln->GetY())+(zn-cln->GetZ())*(zn-cln->GetZ()));
3681 erry = (dist+0.3)*cln->GetSigmaY2()/TMath::Sqrt(cln->GetQ())*(1.+1./(ddrow));
3682 errz = (dist+0.3)*cln->GetSigmaZ2()/TMath::Sqrt(cln->GetQ())*(1.+1./(ddrow));
3683 if (cln->IsUsed(10)) {
3684 // printf("used\n");
3692 polytrack.AddPoint(xn,cln->GetY(),cln->GetZ(),erry, errz);
3697 if ( (sumused>3) || (sumused>0.5*nfound) || (nfound<0.6*nfoundable)) break;
3698 polytrack.UpdateParameters();
3701 if ( (sumused>3) || (sumused>0.5*nfound)) {
3702 //printf("sumused %d\n",sumused);
3707 polytrack.GetFitDerivation(kr0.GetX(),dy,dz);
3708 AliTPCpolyTrack track2;
3710 polytrack.Refit(track2,0.5+TMath::Abs(dy)*0.3,0.4+TMath::Abs(dz)*0.3);
3711 if (track2.GetN()<0.5*nfoundable) continue;
3714 if ((nfound>0.6*nfoundable) &&( nfoundable>0.4*(i1-i2))) {
3716 // test seed with and without constrain
3717 for (Int_t constrain=0; constrain<=0;constrain++){
3718 // add polytrack candidate
3720 Double_t x[5], c[15];
3721 Double_t x1,x2,x3,y1,y2,y3,z1,z2,z3;
3722 track2.GetBoundaries(x3,x1);
3724 track2.GetFitPoint(x1,y1,z1);
3725 track2.GetFitPoint(x2,y2,z2);
3726 track2.GetFitPoint(x3,y3,z3);
3728 //is track pointing to the vertex ?
3731 polytrack.GetFitPoint(x0,y0,z0);
3744 x[4]=F1(x1,y1,x2,y2,x3,y3);
3746 // if (TMath::Abs(x[4]) >= cuts[0]) continue; //
3747 x[2]=F2(x1,y1,x2,y2,x3,y3);
3749 //if (TMath::Abs(x[4]*x1-x[2]) >= cuts[1]) continue;
3750 //x[3]=F3(x1,y1,x2,y2,z1,z2);
3751 x[3]=F3n(x1,y1,x3,y3,z1,z3,x[4]);
3752 //if (TMath::Abs(x[3]) > cuts[3]) continue;
3755 Double_t sy =0.1, sz =0.1;
3756 Double_t sy1=0.02, sz1=0.02;
3757 Double_t sy2=0.02, sz2=0.02;
3761 sy3=25000*x[4]*x[4]+0.1, sy=0.1, sz=0.1;
3764 Double_t f40=(F1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
3765 Double_t f42=(F1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
3766 Double_t f43=(F1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
3767 Double_t f20=(F2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
3768 Double_t f22=(F2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
3769 Double_t f23=(F2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
3771 Double_t f30=(F3(x1,y1+sy,x3,y3,z1,z3)-x[3])/sy;
3772 Double_t f31=(F3(x1,y1,x3,y3,z1+sz,z3)-x[3])/sz;
3773 Double_t f32=(F3(x1,y1,x3,y3+sy,z1,z3)-x[3])/sy;
3774 Double_t f34=(F3(x1,y1,x3,y3,z1,z3+sz)-x[3])/sz;
3779 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3780 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3781 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3782 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3783 c[13]=f30*sy1*f40+f32*sy2*f42;
3784 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3786 //Int_t row1 = fSectors->GetRowNumber(x1);
3787 Int_t row1 = GetRowNumber(x1);
3791 seed->~AliTPCseed();
3792 AliTPCseed *track=new(seed) AliTPCseed(x1,sec*alpha+shift,x,c,index);
3793 track->SetIsSeeding(kTRUE);
3794 Int_t rc=FollowProlongation(*track, i2);
3795 if (constrain) track->SetBConstrain(1);
3797 track->SetBConstrain(0);
3798 track->SetLastPoint(row1+fInnerSec->GetNRows()); // first cluster in track position
3799 track->SetFirstPoint(track->GetLastPoint());
3801 if (rc==0 || track->GetNumberOfClusters()<(i1-i2)*0.5 ||
3802 track->GetNumberOfClusters() < track->GetNFoundable()*0.6 ||
3803 track->GetNShared()>0.4*track->GetNumberOfClusters()) {
3806 seed->~AliTPCseed();
3809 arr->AddLast(track);
3810 seed = new AliTPCseed;
3814 } // if accepted seed
3817 Info("MakeSeeds2","\nSeeding statiistic:\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin3);
3823 AliTPCseed *AliTPCtrackerMI::MakeSeed(AliTPCseed *const track, Float_t r0, Float_t r1, Float_t r2)
3827 //reseed using track points
3828 Int_t p0 = int(r0*track->GetNumberOfClusters()); // point 0
3829 Int_t p1 = int(r1*track->GetNumberOfClusters());
3830 Int_t p2 = int(r2*track->GetNumberOfClusters()); // last point
3832 Double_t x0[3],x1[3],x2[3];
3833 for (Int_t i=0;i<3;i++){
3839 // find track position at given ratio of the length
3840 Int_t sec0=0, sec1=0, sec2=0;
3843 for (Int_t i=0;i<160;i++){
3844 if (track->GetClusterPointer(i)){
3846 AliTPCTrackerPoint *trpoint =track->GetTrackPoint(i);
3847 if ( (index<p0) || x0[0]<0 ){
3848 if (trpoint->GetX()>1){
3849 clindex = track->GetClusterIndex2(i);
3851 x0[0] = trpoint->GetX();
3852 x0[1] = trpoint->GetY();
3853 x0[2] = trpoint->GetZ();
3854 sec0 = ((clindex&0xff000000)>>24)%18;
3859 if ( (index<p1) &&(trpoint->GetX()>1)){
3860 clindex = track->GetClusterIndex2(i);
3862 x1[0] = trpoint->GetX();
3863 x1[1] = trpoint->GetY();
3864 x1[2] = trpoint->GetZ();
3865 sec1 = ((clindex&0xff000000)>>24)%18;
3868 if ( (index<p2) &&(trpoint->GetX()>1)){
3869 clindex = track->GetClusterIndex2(i);
3871 x2[0] = trpoint->GetX();
3872 x2[1] = trpoint->GetY();
3873 x2[2] = trpoint->GetZ();
3874 sec2 = ((clindex&0xff000000)>>24)%18;
3881 Double_t alpha, cs,sn, xx2,yy2;
3883 alpha = (sec1-sec2)*fSectors->GetAlpha();
3884 cs = TMath::Cos(alpha);
3885 sn = TMath::Sin(alpha);
3886 xx2= x1[0]*cs-x1[1]*sn;
3887 yy2= x1[0]*sn+x1[1]*cs;
3891 alpha = (sec0-sec2)*fSectors->GetAlpha();
3892 cs = TMath::Cos(alpha);
3893 sn = TMath::Sin(alpha);
3894 xx2= x0[0]*cs-x0[1]*sn;
3895 yy2= x0[0]*sn+x0[1]*cs;
3901 Double_t x[5],c[15];
3905 x[4]=F1(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]);
3906 // if (x[4]>1) return 0;
3907 x[2]=F2(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]);
3908 x[3]=F3n(x2[0],x2[1],x0[0],x0[1],x2[2],x0[2],x[4]);
3909 //if (TMath::Abs(x[3]) > 2.2) return 0;
3910 //if (TMath::Abs(x[2]) > 1.99) return 0;
3912 Double_t sy =0.1, sz =0.1;
3914 Double_t sy1=0.02+track->GetSigmaY2(), sz1=0.02+track->GetSigmaZ2();
3915 Double_t sy2=0.01+track->GetSigmaY2(), sz2=0.01+track->GetSigmaZ2();
3916 Double_t sy3=0.01+track->GetSigmaY2();
3918 Double_t f40=(F1(x2[0],x2[1]+sy,x1[0],x1[1],x0[0],x0[1])-x[4])/sy;
3919 Double_t f42=(F1(x2[0],x2[1],x1[0],x1[1]+sy,x0[0],x0[1])-x[4])/sy;
3920 Double_t f43=(F1(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]+sy)-x[4])/sy;
3921 Double_t f20=(F2(x2[0],x2[1]+sy,x1[0],x1[1],x0[0],x0[1])-x[2])/sy;
3922 Double_t f22=(F2(x2[0],x2[1],x1[0],x1[1]+sy,x0[0],x0[1])-x[2])/sy;
3923 Double_t f23=(F2(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]+sy)-x[2])/sy;
3925 Double_t f30=(F3(x2[0],x2[1]+sy,x0[0],x0[1],x2[2],x0[2])-x[3])/sy;
3926 Double_t f31=(F3(x2[0],x2[1],x0[0],x0[1],x2[2]+sz,x0[2])-x[3])/sz;
3927 Double_t f32=(F3(x2[0],x2[1],x0[0],x0[1]+sy,x2[2],x0[2])-x[3])/sy;
3928 Double_t f34=(F3(x2[0],x2[1],x0[0],x0[1],x2[2],x0[2]+sz)-x[3])/sz;
3933 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3934 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3935 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3936 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3937 c[13]=f30*sy1*f40+f32*sy2*f42;
3938 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3940 // Int_t row1 = fSectors->GetRowNumber(x2[0]);
3941 AliTPCseed *seed=new AliTPCseed(x2[0], sec2*fSectors->GetAlpha()+fSectors->GetAlphaShift(), x, c, 0);
3942 // Double_t y0,z0,y1,z1, y2,z2;
3943 //seed->GetProlongation(x0[0],y0,z0);
3944 // seed->GetProlongation(x1[0],y1,z1);
3945 //seed->GetProlongation(x2[0],y2,z2);
3947 seed->SetLastPoint(pp2);
3948 seed->SetFirstPoint(pp2);
3955 AliTPCseed *AliTPCtrackerMI::ReSeed(const AliTPCseed *track, Float_t r0, Float_t r1, Float_t r2)
3959 //reseed using founded clusters
3961 // Find the number of clusters
3962 Int_t nclusters = 0;
3963 for (Int_t irow=0;irow<160;irow++){
3964 if (track->GetClusterIndex(irow)>0) nclusters++;
3968 ipos[0] = TMath::Max(int(r0*nclusters),0); // point 0 cluster
3969 ipos[1] = TMath::Min(int(r1*nclusters),nclusters-1); //
3970 ipos[2] = TMath::Min(int(r2*nclusters),nclusters-1); // last point
3973 Double_t xyz[3][3]={{0}};
3974 Int_t row[3]={0},sec[3]={0,0,0};
3976 // find track row position at given ratio of the length
3978 for (Int_t irow=0;irow<160;irow++){
3979 if (track->GetClusterIndex2(irow)<0) continue;
3981 for (Int_t ipoint=0;ipoint<3;ipoint++){
3982 if (index<=ipos[ipoint]) row[ipoint] = irow;
3986 //Get cluster and sector position
3987 for (Int_t ipoint=0;ipoint<3;ipoint++){
3988 Int_t clindex = track->GetClusterIndex2(row[ipoint]);
3989 AliTPCclusterMI * cl = GetClusterMI(clindex);
3992 // AliTPCclusterMI * cl = GetClusterMI(clindex);
3995 sec[ipoint] = ((clindex&0xff000000)>>24)%18;
3996 xyz[ipoint][0] = GetXrow(row[ipoint]);
3997 xyz[ipoint][1] = cl->GetY();
3998 xyz[ipoint][2] = cl->GetZ();
4002 // Calculate seed state vector and covariance matrix
4004 Double_t alpha, cs,sn, xx2,yy2;
4006 alpha = (sec[1]-sec[2])*fSectors->GetAlpha();
4007 cs = TMath::Cos(alpha);
4008 sn = TMath::Sin(alpha);
4009 xx2= xyz[1][0]*cs-xyz[1][1]*sn;
4010 yy2= xyz[1][0]*sn+xyz[1][1]*cs;
4014 alpha = (sec[0]-sec[2])*fSectors->GetAlpha();
4015 cs = TMath::Cos(alpha);
4016 sn = TMath::Sin(alpha);
4017 xx2= xyz[0][0]*cs-xyz[0][1]*sn;
4018 yy2= xyz[0][0]*sn+xyz[0][1]*cs;
4024 Double_t x[5],c[15];
4028 x[4]=F1(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
4029 x[2]=F2(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
4030 x[3]=F3n(xyz[2][0],xyz[2][1],xyz[0][0],xyz[0][1],xyz[2][2],xyz[0][2],x[4]);
4032 Double_t sy =0.1, sz =0.1;
4034 Double_t sy1=0.2, sz1=0.2;
4035 Double_t sy2=0.2, sz2=0.2;
4038 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;
4039 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;
4040 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;
4041 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;
4042 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;
4043 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;
4045 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;
4046 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;
4047 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;
4048 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;
4053 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
4054 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
4055 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
4056 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
4057 c[13]=f30*sy1*f40+f32*sy2*f42;
4058 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
4060 // Int_t row1 = fSectors->GetRowNumber(xyz[2][0]);
4061 AliTPCseed *seed=new AliTPCseed(xyz[2][0], sec[2]*fSectors->GetAlpha()+fSectors->GetAlphaShift(), x, c, 0);
4062 seed->SetLastPoint(row[2]);
4063 seed->SetFirstPoint(row[2]);
4068 AliTPCseed *AliTPCtrackerMI::ReSeed(AliTPCseed *track,Int_t r0, Bool_t forward)
4072 //reseed using founded clusters
4075 Int_t row[3]={0,0,0};
4076 Int_t sec[3]={0,0,0};
4078 // forward direction
4080 for (Int_t irow=r0;irow<160;irow++){
4081 if (track->GetClusterIndex(irow)>0){
4086 for (Int_t irow=160;irow>r0;irow--){
4087 if (track->GetClusterIndex(irow)>0){
4092 for (Int_t irow=row[2]-15;irow>row[0];irow--){
4093 if (track->GetClusterIndex(irow)>0){
4101 for (Int_t irow=0;irow<r0;irow++){
4102 if (track->GetClusterIndex(irow)>0){
4107 for (Int_t irow=r0;irow>0;irow--){
4108 if (track->GetClusterIndex(irow)>0){
4113 for (Int_t irow=row[2]-15;irow>row[0];irow--){
4114 if (track->GetClusterIndex(irow)>0){
4121 if ((row[2]-row[0])<20) return 0;
4122 if (row[1]==0) return 0;
4125 //Get cluster and sector position
4126 for (Int_t ipoint=0;ipoint<3;ipoint++){
4127 Int_t clindex = track->GetClusterIndex2(row[ipoint]);
4128 AliTPCclusterMI * cl = GetClusterMI(clindex);
4131 // AliTPCclusterMI * cl = GetClusterMI(clindex);
4134 sec[ipoint] = ((clindex&0xff000000)>>24)%18;
4135 xyz[ipoint][0] = GetXrow(row[ipoint]);
4136 AliTPCTrackerPoint * point = track->GetTrackPoint(row[ipoint]);
4137 if (point&&ipoint<2){
4139 xyz[ipoint][1] = point->GetY();
4140 xyz[ipoint][2] = point->GetZ();
4143 xyz[ipoint][1] = cl->GetY();
4144 xyz[ipoint][2] = cl->GetZ();
4151 // Calculate seed state vector and covariance matrix
4153 Double_t alpha, cs,sn, xx2,yy2;
4155 alpha = (sec[1]-sec[2])*fSectors->GetAlpha();
4156 cs = TMath::Cos(alpha);
4157 sn = TMath::Sin(alpha);
4158 xx2= xyz[1][0]*cs-xyz[1][1]*sn;
4159 yy2= xyz[1][0]*sn+xyz[1][1]*cs;
4163 alpha = (sec[0]-sec[2])*fSectors->GetAlpha();
4164 cs = TMath::Cos(alpha);
4165 sn = TMath::Sin(alpha);
4166 xx2= xyz[0][0]*cs-xyz[0][1]*sn;
4167 yy2= xyz[0][0]*sn+xyz[0][1]*cs;
4173 Double_t x[5],c[15];
4177 x[4]=F1(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
4178 x[2]=F2(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
4179 x[3]=F3n(xyz[2][0],xyz[2][1],xyz[0][0],xyz[0][1],xyz[2][2],xyz[0][2],x[4]);
4181 Double_t sy =0.1, sz =0.1;
4183 Double_t sy1=0.2, sz1=0.2;
4184 Double_t sy2=0.2, sz2=0.2;
4187 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;
4188 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;
4189 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;
4190 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;
4191 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;
4192 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;
4194 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;
4195 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;
4196 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;
4197 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;
4202 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
4203 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
4204 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
4205 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
4206 c[13]=f30*sy1*f40+f32*sy2*f42;
4207 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
4209 // Int_t row1 = fSectors->GetRowNumber(xyz[2][0]);
4210 AliTPCseed *seed=new AliTPCseed(xyz[2][0], sec[2]*fSectors->GetAlpha()+fSectors->GetAlphaShift(), x, c, 0);
4211 seed->SetLastPoint(row[2]);
4212 seed->SetFirstPoint(row[2]);
4213 for (Int_t i=row[0];i<row[2];i++){
4214 seed->SetClusterIndex(i, track->GetClusterIndex(i));
4222 void AliTPCtrackerMI::FindMultiMC(const TObjArray * array, AliESDEvent */*esd*/, Int_t iter)
4225 // find multi tracks - THIS FUNCTION IS ONLY FOR DEBUG PURPOSES
4227 // Use AliTPCReconstructor::StreamLevel()>2 if you want to tune parameters - cuts
4229 // Two reasons to have multiple find tracks
4230 // 1. Curling tracks can be find more than once
4231 // 2. Splitted tracks
4232 // a.) Multiple seeding to increase tracking efficiency - (~ 100% reached)
4233 // b.) Edge effect on the sector boundaries
4236 // Algorithm done in 2 phases - because of CPU consumption
4237 // it is n^2 algorithm - for lead-lead 20000x20000 combination are investigated
4239 // Algorihm for curling tracks sign:
4240 // 1 phase -makes a very rough fast cuts to minimize combinatorics
4241 // a.) opposite sign
4242 // b.) one of the tracks - not pointing to the primary vertex -
4243 // c.) delta tan(theta)
4245 // 2 phase - calculates DCA between tracks - time consument
4250 // General cuts - for splitted tracks and for curling tracks
4252 const Float_t kMaxdPhi = 0.2; // maximal distance in phi
4254 // Curling tracks cuts
4259 Int_t nentries = array->GetEntriesFast();
4260 AliHelix *helixes = new AliHelix[nentries];
4261 Float_t *xm = new Float_t[nentries];
4262 Float_t *dz0 = new Float_t[nentries];
4263 Float_t *dz1 = new Float_t[nentries];
4269 // Find track COG in x direction - point with best defined parameters
4271 for (Int_t i=0;i<nentries;i++){
4272 AliTPCseed* track = (AliTPCseed*)array->At(i);
4273 if (!track) continue;
4274 track->SetCircular(0);
4275 new (&helixes[i]) AliHelix(*track);
4279 track->GetDZ(GetX(),GetY(),GetZ(),GetBz(),dz);
4282 for (Int_t icl=0; icl<160; icl++){
4283 AliTPCclusterMI * cl = track->GetClusterPointer(icl);
4289 if (ncl>0) xm[i]/=Float_t(ncl);
4292 for (Int_t i0=0;i0<nentries;i0++){
4293 AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
4294 if (!track0) continue;
4295 Float_t xc0 = helixes[i0].GetHelix(6);
4296 Float_t yc0 = helixes[i0].GetHelix(7);
4297 Float_t r0 = helixes[i0].GetHelix(8);
4298 Float_t rc0 = TMath::Sqrt(xc0*xc0+yc0*yc0);
4299 Float_t fi0 = TMath::ATan2(yc0,xc0);
4301 for (Int_t i1=i0+1;i1<nentries;i1++){
4302 AliTPCseed * track1 = (AliTPCseed*)array->At(i1);
4303 if (!track1) continue;
4304 Int_t lab0=track0->GetLabel();
4305 Int_t lab1=track1->GetLabel();
4306 if (TMath::Abs(lab0)!=TMath::Abs(lab1)) continue;
4308 Float_t xc1 = helixes[i1].GetHelix(6);
4309 Float_t yc1 = helixes[i1].GetHelix(7);
4310 Float_t r1 = helixes[i1].GetHelix(8);
4311 Float_t rc1 = TMath::Sqrt(xc1*xc1+yc1*yc1);
4312 Float_t fi1 = TMath::ATan2(yc1,xc1);
4314 Float_t dfi = fi0-fi1;
4317 if (dfi>1.5*TMath::Pi()) dfi-=TMath::Pi(); // take care about edge effect
4318 if (dfi<-1.5*TMath::Pi()) dfi+=TMath::Pi(); //
4319 if (TMath::Abs(dfi)>kMaxdPhi&&helixes[i0].GetHelix(4)*helixes[i1].GetHelix(4)<0){
4321 // if short tracks with undefined sign
4322 fi1 = -TMath::ATan2(yc1,-xc1);
4325 Float_t dtheta = TMath::Abs(track0->GetTgl()-track1->GetTgl())<TMath::Abs(track0->GetTgl()+track1->GetTgl())? track0->GetTgl()-track1->GetTgl():track0->GetTgl()+track1->GetTgl();
4328 // debug stream to tune "fast cuts"
4330 Double_t dist[3]; // distance at X
4331 Double_t mdist[3]={0,0,0}; // mean distance X+-40cm
4332 track0->GetDistance(track1,0.5*(xm[i0]+xm[i1])-40.,dist,AliTracker::GetBz());
4333 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[i]);
4334 track0->GetDistance(track1,0.5*(xm[i0]+xm[i1])+40.,dist,AliTracker::GetBz());
4335 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[i]);
4336 track0->GetDistance(track1,0.5*(xm[i0]+xm[i1]),dist,AliTracker::GetBz());
4337 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[i]);
4338 for (Int_t i=0;i<3;i++) mdist[i]*=0.33333;
4342 for (Int_t icl=0; icl<160; icl++){
4343 AliTPCclusterMI * cl0 = track0->GetClusterPointer(icl);
4344 AliTPCclusterMI * cl1 = track1->GetClusterPointer(icl);
4347 if (cl0==cl1) sums++;
4351 if (AliTPCReconstructor::StreamLevel()>5) {
4352 TTreeSRedirector &cstream = *fDebugStreamer;
4357 "Tr0.="<<track0<< // seed0
4358 "Tr1.="<<track1<< // seed1
4359 "h0.="<<&helixes[i0]<<
4360 "h1.="<<&helixes[i1]<<
4362 "sum="<<sum<< //the sum of rows with cl in both
4363 "sums="<<sums<< //the sum of shared clusters
4364 "xm0="<<xm[i0]<< // the center of track
4365 "xm1="<<xm[i1]<< // the x center of track
4366 // General cut variables
4367 "dfi="<<dfi<< // distance in fi angle
4368 "dtheta="<<dtheta<< // distance int theta angle
4374 "dist0="<<dist[0]<< //distance x
4375 "dist1="<<dist[1]<< //distance y
4376 "dist2="<<dist[2]<< //distance z
4377 "mdist0="<<mdist[0]<< //distance x
4378 "mdist1="<<mdist[1]<< //distance y
4379 "mdist2="<<mdist[2]<< //distance z
4395 if (AliTPCReconstructor::StreamLevel()>1) {
4396 AliInfo("Time for curling tracks removal DEBUGGING MC");
4403 void AliTPCtrackerMI::FindSplitted(TObjArray * array, AliESDEvent */*esd*/, Int_t /*iter*/){
4405 // Find Splitted tracks and remove the one with worst quality
4406 // Corresponding debug streamer to tune selections - "Splitted2"
4408 // 0. Sort tracks according quility
4409 // 1. Propagate the tracks to the reference radius
4410 // 2. Double_t loop to select close tracks (only to speed up process)
4411 // 3. Calculate cluster overlap ratio - and remove the track if bigger than a threshold
4412 // 4. Delete temporary parameters
4414 const Double_t xref=GetXrow(63); // reference radius -IROC/OROC boundary
4416 const Double_t kCutP1=10; // delta Z cut 10 cm
4417 const Double_t kCutP2=0.15; // delta snp(fi) cut 0.15
4418 const Double_t kCutP3=0.15; // delta tgl(theta) cut 0.15
4419 const Double_t kCutAlpha=0.15; // delta alpha cut
4420 Int_t firstpoint = 0;
4421 Int_t lastpoint = 160;
4423 Int_t nentries = array->GetEntriesFast();
4424 AliExternalTrackParam *params = new AliExternalTrackParam[nentries];
4430 //0. Sort tracks according quality
4431 //1. Propagate the ext. param to reference radius
4432 Int_t nseed = array->GetEntriesFast();
4433 if (nseed<=0) return;
4434 Float_t * quality = new Float_t[nseed];
4435 Int_t * indexes = new Int_t[nseed];
4436 for (Int_t i=0; i<nseed; i++) {
4437 AliTPCseed *pt=(AliTPCseed*)array->UncheckedAt(i);
4442 pt->UpdatePoints(); //select first last max dens points
4443 Float_t * points = pt->GetPoints();
4444 if (points[3]<0.8) quality[i] =-1;
4445 quality[i] = (points[2]-points[0])+pt->GetNumberOfClusters();
4446 //prefer high momenta tracks if overlaps
4447 quality[i] *= TMath::Sqrt(TMath::Abs(pt->Pt())+0.5);
4449 AliTracker::PropagateTrackToBxByBz(&(params[i]),xref,pt->GetMass(),5.,kTRUE);
4450 AliTracker::PropagateTrackToBxByBz(&(params[i]),xref,pt->GetMass(),1.,kTRUE);
4452 TMath::Sort(nseed,quality,indexes);
4454 // 3. Loop over pair of tracks
4456 for (Int_t i0=0; i0<nseed; i0++) {
4457 Int_t index0=indexes[i0];
4458 if (!(array->UncheckedAt(index0))) continue;
4459 AliTPCseed *s1 = (AliTPCseed*)array->UncheckedAt(index0);
4460 if (!s1->IsActive()) continue;
4461 AliExternalTrackParam &par0=params[index0];
4462 for (Int_t i1=i0+1; i1<nseed; i1++) {
4463 Int_t index1=indexes[i1];
4464 if (!(array->UncheckedAt(index1))) continue;
4465 AliTPCseed *s2 = (AliTPCseed*)array->UncheckedAt(index1);
4466 if (!s2->IsActive()) continue;
4467 if (s2->GetKinkIndexes()[0]!=0)
4468 if (s2->GetKinkIndexes()[0] == -s1->GetKinkIndexes()[0]) continue;
4469 AliExternalTrackParam &par1=params[index1];
4470 if (TMath::Abs(par0.GetParameter()[3]-par1.GetParameter()[3])>kCutP3) continue;
4471 if (TMath::Abs(par0.GetParameter()[1]-par1.GetParameter()[1])>kCutP1) continue;
4472 if (TMath::Abs(par0.GetParameter()[2]-par1.GetParameter()[2])>kCutP2) continue;
4473 Double_t dAlpha= TMath::Abs(par0.GetAlpha()-par1.GetAlpha());
4474 if (dAlpha>TMath::Pi()) dAlpha-=TMath::Pi();
4475 if (TMath::Abs(dAlpha)>kCutAlpha) continue;
4480 Int_t firstShared=lastpoint, lastShared=firstpoint;
4481 Int_t firstRow=lastpoint, lastRow=firstpoint;
4483 for (Int_t i=firstpoint;i<lastpoint;i++){
4484 if (s1->GetClusterIndex2(i)>0) nall0++;
4485 if (s2->GetClusterIndex2(i)>0) nall1++;
4486 if (s1->GetClusterIndex2(i)>0 && s2->GetClusterIndex2(i)>0) {
4487 if (i<firstRow) firstRow=i;
4488 if (i>lastRow) lastRow=i;
4490 if ( (s1->GetClusterIndex2(i))==(s2->GetClusterIndex2(i)) && s1->GetClusterIndex2(i)>0) {
4491 if (i<firstShared) firstShared=i;
4492 if (i>lastShared) lastShared=i;
4496 Double_t ratio0 = Float_t(sumShared)/Float_t(TMath::Min(nall0+1,nall1+1));
4497 Double_t ratio1 = Float_t(sumShared)/Float_t(TMath::Max(nall0+1,nall1+1));
4499 if( AliTPCReconstructor::StreamLevel()>1){
4500 TTreeSRedirector &cstream = *fDebugStreamer;
4501 Int_t n0=s1->GetNumberOfClusters();
4502 Int_t n1=s2->GetNumberOfClusters();
4503 Int_t n0F=s1->GetNFoundable();
4504 Int_t n1F=s2->GetNFoundable();
4505 Int_t lab0=s1->GetLabel();
4506 Int_t lab1=s2->GetLabel();
4508 cstream<<"Splitted2"<<
4509 "iter="<<fIteration<<
4510 "lab0="<<lab0<< // MC label if exist
4511 "lab1="<<lab1<< // MC label if exist
4514 "ratio0="<<ratio0<< // shared ratio
4515 "ratio1="<<ratio1<< // shared ratio
4516 "p0.="<<&par0<< // track parameters
4518 "s0.="<<s1<< // full seed
4520 "n0="<<n0<< // number of clusters track 0
4521 "n1="<<n1<< // number of clusters track 1
4522 "nall0="<<nall0<< // number of clusters track 0
4523 "nall1="<<nall1<< // number of clusters track 1
4524 "n0F="<<n0F<< // number of findable
4525 "n1F="<<n1F<< // number of findable
4526 "shared="<<sumShared<< // number of shared clusters
4527 "firstS="<<firstShared<< // first and the last shared row
4528 "lastS="<<lastShared<<
4529 "firstRow="<<firstRow<< // first and the last row with cluster
4530 "lastRow="<<lastRow<< //
4534 // remove track with lower quality
4536 if (ratio0>AliTPCReconstructor::GetRecoParam()->GetCutSharedClusters(0) ||
4537 ratio1>AliTPCReconstructor::GetRecoParam()->GetCutSharedClusters(1)){
4541 delete array->RemoveAt(index1);
4546 // 4. Delete temporary array
4556 void AliTPCtrackerMI::FindCurling(const TObjArray * array, AliESDEvent */*esd*/, Int_t iter)
4559 // find Curling tracks
4560 // Use AliTPCReconstructor::StreamLevel()>1 if you want to tune parameters - cuts
4563 // Algorithm done in 2 phases - because of CPU consumption
4564 // it is n^2 algorithm - for lead-lead 20000x20000 combination are investigated
4565 // see detal in MC part what can be used to cut
4569 const Float_t kMaxC = 400; // maximal curvature to of the track
4570 const Float_t kMaxdTheta = 0.15; // maximal distance in theta
4571 const Float_t kMaxdPhi = 0.15; // maximal distance in phi
4572 const Float_t kPtRatio = 0.3; // ratio between pt
4573 const Float_t kMinDCAR = 2.; // distance to the primary vertex in r - see cpipe cut
4576 // Curling tracks cuts
4579 const Float_t kMaxDeltaRMax = 40; // distance in outer radius
4580 const Float_t kMaxDeltaRMin = 5.; // distance in lower radius - see cpipe cut
4581 const Float_t kMinAngle = 2.9; // angle between tracks
4582 const Float_t kMaxDist = 5; // biggest distance
4584 // The cuts can be tuned using the "MC information stored in Multi tree ==> see FindMultiMC
4587 TCut csign("csign","Tr0.fP[4]*Tr1.fP[4]<0"); //opposite sign
4588 TCut cmax("cmax","abs(Tr0.GetC())>1/400");
4589 TCut cda("cda","sqrt(dtheta^2+dfi^2)<0.15");
4590 TCut ccratio("ccratio","abs((Tr0.fP[4]+Tr1.fP[4])/(abs(Tr0.fP[4])+abs(Tr1.fP[4])))<0.3");
4591 TCut cpipe("cpipe", "min(abs(r0-rc0),abs(r1-rc1))>5");
4593 TCut cdrmax("cdrmax","abs(abs(rc0+r0)-abs(rc1+r1))<40")
4594 TCut cdrmin("cdrmin","abs(abs(rc0+r0)-abs(rc1+r1))<10")
4596 Multi->Draw("dfi","iter==0"+csign+cmax+cda+ccratio); ~94% of curling tracks fulfill
4597 Multi->Draw("min(abs(r0-rc0),abs(r1-rc1))","iter==0&&abs(lab1)==abs(lab0)"+csign+cmax+cda+ccratio+cpipe+cdrmin+cdrmax); //80%
4599 Curling2->Draw("dfi","iter==0&&abs(lab0)==abs(lab1)"+csign+cmax+cdtheta+cdfi+ccratio)
4605 Int_t nentries = array->GetEntriesFast();
4606 AliHelix *helixes = new AliHelix[nentries];
4607 for (Int_t i=0;i<nentries;i++){
4608 AliTPCseed* track = (AliTPCseed*)array->At(i);
4609 if (!track) continue;
4610 track->SetCircular(0);
4611 new (&helixes[i]) AliHelix(*track);
4617 Double_t phase[2][2]={{0,0},{0,0}},radius[2]={0,0};
4623 for (Int_t i0=0;i0<nentries;i0++){
4624 AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
4625 if (!track0) continue;
4626 if (TMath::Abs(track0->GetC())<1/kMaxC) continue;
4627 Float_t xc0 = helixes[i0].GetHelix(6);
4628 Float_t yc0 = helixes[i0].GetHelix(7);
4629 Float_t r0 = helixes[i0].GetHelix(8);
4630 Float_t rc0 = TMath::Sqrt(xc0*xc0+yc0*yc0);
4631 Float_t fi0 = TMath::ATan2(yc0,xc0);
4633 for (Int_t i1=i0+1;i1<nentries;i1++){
4634 AliTPCseed * track1 = (AliTPCseed*)array->At(i1);
4635 if (!track1) continue;
4636 if (TMath::Abs(track1->GetC())<1/kMaxC) continue;
4637 Float_t xc1 = helixes[i1].GetHelix(6);
4638 Float_t yc1 = helixes[i1].GetHelix(7);
4639 Float_t r1 = helixes[i1].GetHelix(8);
4640 Float_t rc1 = TMath::Sqrt(xc1*xc1+yc1*yc1);
4641 Float_t fi1 = TMath::ATan2(yc1,xc1);
4643 Float_t dfi = fi0-fi1;
4646 if (dfi>1.5*TMath::Pi()) dfi-=TMath::Pi(); // take care about edge effect
4647 if (dfi<-1.5*TMath::Pi()) dfi+=TMath::Pi(); //
4648 Float_t dtheta = TMath::Abs(track0->GetTgl()-track1->GetTgl())<TMath::Abs(track0->GetTgl()+track1->GetTgl())? track0->GetTgl()-track1->GetTgl():track0->GetTgl()+track1->GetTgl();
4652 if (track0->GetBConstrain()&&track1->GetBConstrain()) continue; // not constrained
4653 if (track1->GetSigned1Pt()*track0->GetSigned1Pt()>0) continue; // not the same sign
4654 if ( TMath::Abs(track1->GetTgl()+track0->GetTgl())>kMaxdTheta) continue; //distance in the Theta
4655 if ( TMath::Abs(dfi)>kMaxdPhi) continue; //distance in phi
4656 if ( TMath::Sqrt(dfi*dfi+dtheta*dtheta)>kMaxdPhi) continue; //common angular offset
4658 Float_t pt0 = track0->GetSignedPt();
4659 Float_t pt1 = track1->GetSignedPt();
4660 if ((TMath::Abs(pt0+pt1)/(TMath::Abs(pt0)+TMath::Abs(pt1)))>kPtRatio) continue;
4661 if ((iter==1) && TMath::Abs(TMath::Abs(rc0+r0)-TMath::Abs(rc1+r1))>kMaxDeltaRMax) continue;
4662 if ((iter!=1) &&TMath::Abs(TMath::Abs(rc0-r0)-TMath::Abs(rc1-r1))>kMaxDeltaRMin) continue;
4663 if (TMath::Min(TMath::Abs(rc0-r0),TMath::Abs(rc1-r1))<kMinDCAR) continue;
4666 // Now find closest approach
4670 Int_t npoints = helixes[i0].GetRPHIintersections(helixes[i1], phase, radius,10);
4671 if (npoints==0) continue;
4672 helixes[i0].GetClosestPhases(helixes[i1], phase);
4676 Double_t hangles[3];
4677 helixes[i0].Evaluate(phase[0][0],xyz0);
4678 helixes[i1].Evaluate(phase[0][1],xyz1);
4680 helixes[i0].GetAngle(phase[0][0],helixes[i1],phase[0][1],hangles);
4681 Double_t deltah[2],deltabest;
4682 if (TMath::Abs(hangles[2])<kMinAngle) continue;
4686 helixes[i0].ParabolicDCA(helixes[i1],phase[0][0],phase[0][1],radius[0],deltah[0],2);
4688 helixes[i0].ParabolicDCA(helixes[i1],phase[1][0],phase[1][1],radius[1],deltah[1],2);
4689 if (deltah[1]<deltah[0]) ibest=1;
4691 deltabest = TMath::Sqrt(deltah[ibest]);
4692 helixes[i0].Evaluate(phase[ibest][0],xyz0);
4693 helixes[i1].Evaluate(phase[ibest][1],xyz1);
4694 helixes[i0].GetAngle(phase[ibest][0],helixes[i1],phase[ibest][1],hangles);
4695 Double_t radiusbest = TMath::Sqrt(radius[ibest]);
4697 if (deltabest>kMaxDist) continue;
4698 // if (mindcar+mindcaz<40 && (TMath::Abs(hangles[2])<kMinAngle ||deltabest>3)) continue;
4699 Bool_t sign =kFALSE;
4700 if (hangles[2]>kMinAngle) sign =kTRUE;
4703 // circular[i0] = kTRUE;
4704 // circular[i1] = kTRUE;
4705 if (track0->OneOverPt()<track1->OneOverPt()){
4706 track0->SetCircular(track0->GetCircular()+1);
4707 track1->SetCircular(track1->GetCircular()+2);
4710 track1->SetCircular(track1->GetCircular()+1);
4711 track0->SetCircular(track0->GetCircular()+2);
4714 if (AliTPCReconstructor::StreamLevel()>2){
4716 //debug stream to tune "fine" cuts
4717 Int_t lab0=track0->GetLabel();
4718 Int_t lab1=track1->GetLabel();
4719 TTreeSRedirector &cstream = *fDebugStreamer;
4720 cstream<<"Curling2"<<
4736 "npoints="<<npoints<<
4737 "hangles0="<<hangles[0]<<
4738 "hangles1="<<hangles[1]<<
4739 "hangles2="<<hangles[2]<<
4742 "radius="<<radiusbest<<
4743 "deltabest="<<deltabest<<
4744 "phase0="<<phase[ibest][0]<<
4745 "phase1="<<phase[ibest][1]<<
4753 if (AliTPCReconstructor::StreamLevel()>1) {
4754 AliInfo("Time for curling tracks removal");
4763 void AliTPCtrackerMI::FindKinks(TObjArray * array, AliESDEvent *esd)
4770 TObjArray *kinks= new TObjArray(10000);
4771 // TObjArray *v0s= new TObjArray(10000);
4772 Int_t nentries = array->GetEntriesFast();
4773 AliHelix *helixes = new AliHelix[nentries];
4774 Int_t *sign = new Int_t[nentries];
4775 Int_t *nclusters = new Int_t[nentries];
4776 Float_t *alpha = new Float_t[nentries];
4777 AliKink *kink = new AliKink();
4778 Int_t * usage = new Int_t[nentries];
4779 Float_t *zm = new Float_t[nentries];
4780 Float_t *z0 = new Float_t[nentries];
4781 Float_t *fim = new Float_t[nentries];
4782 Float_t *shared = new Float_t[nentries];
4783 Bool_t *circular = new Bool_t[nentries];
4784 Float_t *dca = new Float_t[nentries];
4785 //const AliESDVertex * primvertex = esd->GetVertex();
4787 // nentries = array->GetEntriesFast();
4792 for (Int_t i=0;i<nentries;i++){
4795 AliTPCseed* track = (AliTPCseed*)array->At(i);
4796 if (!track) continue;
4797 track->SetCircular(0);
4799 track->UpdatePoints();
4800 if (( track->GetPoints()[2]- track->GetPoints()[0])>5 && track->GetPoints()[3]>0.8){
4802 nclusters[i]=track->GetNumberOfClusters();
4803 alpha[i] = track->GetAlpha();
4804 new (&helixes[i]) AliHelix(*track);
4806 helixes[i].Evaluate(0,xyz);
4807 sign[i] = (track->GetC()>0) ? -1:1;
4810 if (track->GetProlongation(x,y,z)){
4812 fim[i] = alpha[i]+TMath::ATan2(y,x);
4815 zm[i] = track->GetZ();
4819 circular[i]= kFALSE;
4820 if (track->GetProlongation(0,y,z)) z0[i] = z;
4821 dca[i] = track->GetD(0,0);
4827 Int_t ncandidates =0;
4830 Double_t phase[2][2]={{0,0},{0,0}},radius[2]={0,0};
4833 // Find circling track
4835 for (Int_t i0=0;i0<nentries;i0++){
4836 AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
4837 if (!track0) continue;
4838 if (track0->GetNumberOfClusters()<40) continue;
4839 if (TMath::Abs(1./track0->GetC())>200) continue;
4840 for (Int_t i1=i0+1;i1<nentries;i1++){
4841 AliTPCseed * track1 = (AliTPCseed*)array->At(i1);
4842 if (!track1) continue;
4843 if (track1->GetNumberOfClusters()<40) continue;
4844 if ( TMath::Abs(track1->GetTgl()+track0->GetTgl())>0.1) continue;
4845 if (track0->GetBConstrain()&&track1->GetBConstrain()) continue;
4846 if (TMath::Abs(1./track1->GetC())>200) continue;
4847 if (track1->GetSigned1Pt()*track0->GetSigned1Pt()>0) continue;
4848 if (track1->GetTgl()*track0->GetTgl()>0) continue;
4849 if (TMath::Max(TMath::Abs(1./track0->GetC()),TMath::Abs(1./track1->GetC()))>190) continue;
4850 if (track0->GetBConstrain()&&track1->OneOverPt()<track0->OneOverPt()) continue; //returning - lower momenta
4851 if (track1->GetBConstrain()&&track0->OneOverPt()<track1->OneOverPt()) continue; //returning - lower momenta
4853 Float_t mindcar = TMath::Min(TMath::Abs(dca[i0]),TMath::Abs(dca[i1]));
4854 if (mindcar<5) continue;
4855 Float_t mindcaz = TMath::Min(TMath::Abs(z0[i0]-GetZ()),TMath::Abs(z0[i1]-GetZ()));
4856 if (mindcaz<5) continue;
4857 if (mindcar+mindcaz<20) continue;
4860 Float_t xc0 = helixes[i0].GetHelix(6);
4861 Float_t yc0 = helixes[i0].GetHelix(7);
4862 Float_t r0 = helixes[i0].GetHelix(8);
4863 Float_t xc1 = helixes[i1].GetHelix(6);
4864 Float_t yc1 = helixes[i1].GetHelix(7);
4865 Float_t r1 = helixes[i1].GetHelix(8);
4867 Float_t rmean = (r0+r1)*0.5;
4868 Float_t delta =TMath::Sqrt((xc1-xc0)*(xc1-xc0)+(yc1-yc0)*(yc1-yc0));
4869 //if (delta>30) continue;
4870 if (delta>rmean*0.25) continue;
4871 if (TMath::Abs(r0-r1)/rmean>0.3) continue;
4873 Int_t npoints = helixes[i0].GetRPHIintersections(helixes[i1], phase, radius,10);
4874 if (npoints==0) continue;
4875 helixes[i0].GetClosestPhases(helixes[i1], phase);
4879 Double_t hangles[3];
4880 helixes[i0].Evaluate(phase[0][0],xyz0);
4881 helixes[i1].Evaluate(phase[0][1],xyz1);
4883 helixes[i0].GetAngle(phase[0][0],helixes[i1],phase[0][1],hangles);
4884 Double_t deltah[2],deltabest;
4885 if (hangles[2]<2.8) continue;
4888 helixes[i0].ParabolicDCA(helixes[i1],phase[0][0],phase[0][1],radius[0],deltah[0],2);
4890 helixes[i0].ParabolicDCA(helixes[i1],phase[1][0],phase[1][1],radius[1],deltah[1],2);
4891 if (deltah[1]<deltah[0]) ibest=1;
4893 deltabest = TMath::Sqrt(deltah[ibest]);
4894 helixes[i0].Evaluate(phase[ibest][0],xyz0);
4895 helixes[i1].Evaluate(phase[ibest][1],xyz1);
4896 helixes[i0].GetAngle(phase[ibest][0],helixes[i1],phase[ibest][1],hangles);
4897 Double_t radiusbest = TMath::Sqrt(radius[ibest]);
4899 if (deltabest>6) continue;
4900 if (mindcar+mindcaz<40 && (hangles[2]<3.12||deltabest>3)) continue;
4901 Bool_t lsign =kFALSE;
4902 if (hangles[2]>3.06) lsign =kTRUE;
4905 circular[i0] = kTRUE;
4906 circular[i1] = kTRUE;
4907 if (track0->OneOverPt()<track1->OneOverPt()){
4908 track0->SetCircular(track0->GetCircular()+1);
4909 track1->SetCircular(track1->GetCircular()+2);
4912 track1->SetCircular(track1->GetCircular()+1);
4913 track0->SetCircular(track0->GetCircular()+2);
4916 if (lsign&&AliTPCReconstructor::StreamLevel()>1){
4918 Int_t lab0=track0->GetLabel();
4919 Int_t lab1=track1->GetLabel();
4920 TTreeSRedirector &cstream = *fDebugStreamer;
4921 cstream<<"Curling"<<
4928 "mindcar="<<mindcar<<
4929 "mindcaz="<<mindcaz<<
4932 "npoints="<<npoints<<
4933 "hangles0="<<hangles[0]<<
4934 "hangles2="<<hangles[2]<<
4939 "radius="<<radiusbest<<
4940 "deltabest="<<deltabest<<
4941 "phase0="<<phase[ibest][0]<<
4942 "phase1="<<phase[ibest][1]<<
4952 for (Int_t i =0;i<nentries;i++){
4953 if (sign[i]==0) continue;
4954 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
4961 Double_t cradius0 = 40*40;
4962 Double_t cradius1 = 270*270;
4965 Double_t cdist3=0.55;
4966 for (Int_t j =i+1;j<nentries;j++){
4968 if (sign[j]*sign[i]<1) continue;
4969 if ( (nclusters[i]+nclusters[j])>200) continue;
4970 if ( (nclusters[i]+nclusters[j])<80) continue;
4971 if ( TMath::Abs(zm[i]-zm[j])>60.) continue;
4972 if ( TMath::Abs(fim[i]-fim[j])>0.6 && TMath::Abs(fim[i]-fim[j])<5.7 ) continue;
4973 //AliTPCseed * track1 = (AliTPCseed*)array->At(j); Double_t phase[2][2],radius[2];
4974 Int_t npoints = helixes[i].GetRPHIintersections(helixes[j], phase, radius,20);
4975 if (npoints<1) continue;
4978 if (radius[0]<cradius0||radius[0]>cradius1) continue;
4981 if ( (radius[0]<cradius0||radius[0]>cradius1) && (radius[1]<cradius0||radius[1]>cradius1) ) continue;
4984 Double_t delta1=10000,delta2=10000;
4985 // cuts on the intersection radius
4986 helixes[i].LinearDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta1);
4987 if (radius[0]<20&&delta1<1) continue; //intersection at vertex
4988 if (radius[0]<10&&delta1<3) continue; //intersection at vertex
4990 helixes[i].LinearDCA(helixes[j],phase[1][0],phase[1][1],radius[1],delta2);
4991 if (radius[1]<20&&delta2<1) continue; //intersection at vertex
4992 if (radius[1]<10&&delta2<3) continue; //intersection at vertex
4995 Double_t distance1 = TMath::Min(delta1,delta2);
4996 if (distance1>cdist1) continue; // cut on DCA linear approximation
4998 npoints = helixes[i].GetRPHIintersections(helixes[j], phase, radius,20);
4999 helixes[i].ParabolicDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta1);
5000 if (radius[0]<20&&delta1<1) continue; //intersection at vertex
5001 if (radius[0]<10&&delta1<3) continue; //intersection at vertex
5004 helixes[i].ParabolicDCA(helixes[j],phase[1][0],phase[1][1],radius[1],delta2);
5005 if (radius[1]<20&&delta2<1) continue; //intersection at vertex
5006 if (radius[1]<10&&delta2<3) continue; //intersection at vertex
5008 distance1 = TMath::Min(delta1,delta2);
5011 rkink = TMath::Sqrt(radius[0]);
5014 rkink = TMath::Sqrt(radius[1]);
5016 if (distance1>cdist2) continue;
5019 AliTPCseed * track1 = (AliTPCseed*)array->At(j);
5022 Int_t row0 = GetRowNumber(rkink);
5023 if (row0<10) continue;
5024 if (row0>150) continue;
5027 Float_t dens00=-1,dens01=-1;
5028 Float_t dens10=-1,dens11=-1;
5030 Int_t found,foundable,ishared;
5031 track0->GetClusterStatistic(0,row0-5, found, foundable,ishared,kFALSE);
5032 if (foundable>5) dens00 = Float_t(found)/Float_t(foundable);
5033 track0->GetClusterStatistic(row0+5,155, found, foundable,ishared,kFALSE);
5034 if (foundable>5) dens01 = Float_t(found)/Float_t(foundable);
5036 track1->GetClusterStatistic(0,row0-5, found, foundable,ishared,kFALSE);
5037 if (foundable>10) dens10 = Float_t(found)/Float_t(foundable);
5038 track1->GetClusterStatistic(row0+5,155, found, foundable,ishared,kFALSE);
5039 if (foundable>10) dens11 = Float_t(found)/Float_t(foundable);
5041 if (dens00<dens10 && dens01<dens11) continue;
5042 if (dens00>dens10 && dens01>dens11) continue;
5043 if (TMath::Max(dens00,dens10)<0.1) continue;
5044 if (TMath::Max(dens01,dens11)<0.3) continue;
5046 if (TMath::Min(dens00,dens10)>0.6) continue;
5047 if (TMath::Min(dens01,dens11)>0.6) continue;
5050 AliTPCseed * ktrack0, *ktrack1;
5059 if (TMath::Abs(ktrack0->GetC())>5) continue; // cut on the curvature for mother particle
5060 AliExternalTrackParam paramm(*ktrack0);
5061 AliExternalTrackParam paramd(*ktrack1);
5062 if (row0>60&&ktrack1->GetReference().GetX()>90.)new (¶md) AliExternalTrackParam(ktrack1->GetReference());
5065 kink->SetMother(paramm);
5066 kink->SetDaughter(paramd);
5069 Float_t x[3] = { kink->GetPosition()[0],kink->GetPosition()[1],kink->GetPosition()[2]};
5071 fkParam->Transform0to1(x,index);
5072 fkParam->Transform1to2(x,index);
5073 row0 = GetRowNumber(x[0]);
5075 if (kink->GetR()<100) continue;
5076 if (kink->GetR()>240) continue;
5077 if (kink->GetPosition()[2]/kink->GetR()>AliTPCReconstructor::GetCtgRange()) continue; //out of fiducial volume
5078 if (kink->GetDistance()>cdist3) continue;
5079 Float_t dird = kink->GetDaughterP()[0]*kink->GetPosition()[0]+kink->GetDaughterP()[1]*kink->GetPosition()[1]; // rough direction estimate
5080 if (dird<0) continue;
5082 Float_t dirm = kink->GetMotherP()[0]*kink->GetPosition()[0]+kink->GetMotherP()[1]*kink->GetPosition()[1]; // rough direction estimate
5083 if (dirm<0) continue;
5084 Float_t mpt = TMath::Sqrt(kink->GetMotherP()[0]*kink->GetMotherP()[0]+kink->GetMotherP()[1]*kink->GetMotherP()[1]);
5085 if (mpt<0.2) continue;
5088 //for high momenta momentum not defined well in first iteration
5089 Double_t qt = TMath::Sin(kink->GetAngle(2))*ktrack1->GetP();
5090 if (qt>0.35) continue;
5093 kink->SetLabel(CookLabel(ktrack0,0.4,0,row0),0);
5094 kink->SetLabel(CookLabel(ktrack1,0.4,row0,160),1);
5096 kink->SetTPCDensity(dens00,0,0);
5097 kink->SetTPCDensity(dens01,0,1);
5098 kink->SetTPCDensity(dens10,1,0);
5099 kink->SetTPCDensity(dens11,1,1);
5100 kink->SetIndex(i,0);
5101 kink->SetIndex(j,1);
5104 kink->SetTPCDensity(dens10,0,0);
5105 kink->SetTPCDensity(dens11,0,1);
5106 kink->SetTPCDensity(dens00,1,0);
5107 kink->SetTPCDensity(dens01,1,1);
5108 kink->SetIndex(j,0);
5109 kink->SetIndex(i,1);
5112 if (mpt<1||kink->GetAngle(2)>0.1){
5113 // angle and densities not defined yet
5114 if (kink->GetTPCDensityFactor()<0.8) continue;
5115 if ((2-kink->GetTPCDensityFactor())*kink->GetDistance() >0.25) continue;
5116 if (kink->GetAngle(2)*ktrack0->GetP()<0.003) continue; //too small angle
5117 if (kink->GetAngle(2)>0.2&&kink->GetTPCDensityFactor()<1.15) continue;
5118 if (kink->GetAngle(2)>0.2&&kink->GetTPCDensity(0,1)>0.05) continue;
5120 Float_t criticalangle = track0->GetSigmaSnp2()+track0->GetSigmaTgl2();
5121 criticalangle+= track1->GetSigmaSnp2()+track1->GetSigmaTgl2();
5122 criticalangle= 3*TMath::Sqrt(criticalangle);
5123 if (criticalangle>0.02) criticalangle=0.02;
5124 if (kink->GetAngle(2)<criticalangle) continue;
5127 Int_t drow = Int_t(2.+0.5/(0.05+kink->GetAngle(2))); // overlap region defined
5128 Float_t shapesum =0;
5130 for ( Int_t row = row0-drow; row<row0+drow;row++){
5131 if (row<0) continue;
5132 if (row>155) continue;
5133 if (ktrack0->GetClusterPointer(row)){
5134 AliTPCTrackerPoint *point =ktrack0->GetTrackPoint(row);
5135 shapesum+=point->GetSigmaY()+point->GetSigmaZ();
5138 if (ktrack1->GetClusterPointer(row)){
5139 AliTPCTrackerPoint *point =ktrack1->GetTrackPoint(row);
5140 shapesum+=point->GetSigmaY()+point->GetSigmaZ();
5145 kink->SetShapeFactor(-1.);
5148 kink->SetShapeFactor(shapesum/sum);
5150 // esd->AddKink(kink);
5152 // kink->SetMother(paramm);
5153 //kink->SetDaughter(paramd);
5155 Double_t chi2P2 = paramm.GetParameter()[2]-paramd.GetParameter()[2];
5157 chi2P2/=paramm.GetCovariance()[5]+paramd.GetCovariance()[5];
5158 Double_t chi2P3 = paramm.GetParameter()[3]-paramd.GetParameter()[3];
5160 chi2P3/=paramm.GetCovariance()[9]+paramd.GetCovariance()[9];
5162 if (AliTPCReconstructor::StreamLevel()>1) {
5163 (*fDebugStreamer)<<"kinkLpt"<<
5171 if ( chi2P2+chi2P3<AliTPCReconstructor::GetRecoParam()->GetKinkAngleCutChi2(0)){
5175 kinks->AddLast(kink);
5181 // sort the kinks according quality - and refit them towards vertex
5183 Int_t nkinks = kinks->GetEntriesFast();
5184 Float_t *quality = new Float_t[nkinks];
5185 Int_t *indexes = new Int_t[nkinks];
5186 AliTPCseed *mothers = new AliTPCseed[nkinks];
5187 AliTPCseed *daughters = new AliTPCseed[nkinks];
5190 for (Int_t i=0;i<nkinks;i++){
5192 AliKink *kinkl = (AliKink*)kinks->At(i);
5194 // refit kinks towards vertex
5196 Int_t index0 = kinkl->GetIndex(0);
5197 Int_t index1 = kinkl->GetIndex(1);
5198 AliTPCseed * ktrack0 = (AliTPCseed*)array->At(index0);
5199 AliTPCseed * ktrack1 = (AliTPCseed*)array->At(index1);
5201 Int_t sumn=ktrack0->GetNumberOfClusters()+ktrack1->GetNumberOfClusters();
5203 // Refit Kink under if too small angle
5205 if (kinkl->GetAngle(2)<0.05){
5206 kinkl->SetTPCRow0(GetRowNumber(kinkl->GetR()));
5207 Int_t row0 = kinkl->GetTPCRow0();
5208 Int_t drow = Int_t(2.+0.5/(0.05+kinkl->GetAngle(2)));
5211 Int_t last = row0-drow;
5212 if (last<40) last=40;
5213 if (last<ktrack0->GetFirstPoint()+25) last = ktrack0->GetFirstPoint()+25;
5214 AliTPCseed* seed0 = ReSeed(ktrack0,last,kFALSE);
5217 Int_t first = row0+drow;
5218 if (first>130) first=130;
5219 if (first>ktrack1->GetLastPoint()-25) first = TMath::Max(ktrack1->GetLastPoint()-25,30);
5220 AliTPCseed* seed1 = ReSeed(ktrack1,first,kTRUE);
5222 if (seed0 && seed1){
5223 kinkl->SetStatus(1,8);
5224 if (RefitKink(*seed0,*seed1,*kinkl)) kinkl->SetStatus(1,9);
5225 row0 = GetRowNumber(kinkl->GetR());
5226 sumn = seed0->GetNumberOfClusters()+seed1->GetNumberOfClusters();
5227 mothers[i] = *seed0;
5228 daughters[i] = *seed1;
5231 delete kinks->RemoveAt(i);
5232 if (seed0) delete seed0;
5233 if (seed1) delete seed1;
5236 if (kinkl->GetDistance()>0.5 || kinkl->GetR()<110 || kinkl->GetR()>240) {
5237 delete kinks->RemoveAt(i);
5238 if (seed0) delete seed0;
5239 if (seed1) delete seed1;
5247 if (kinkl) quality[i] = 160*((0.1+kinkl->GetDistance())*(2.-kinkl->GetTPCDensityFactor()))/(sumn+40.); //the longest -clossest will win
5249 TMath::Sort(nkinks,quality,indexes,kFALSE);
5251 //remove double find kinks
5253 for (Int_t ikink0=1;ikink0<nkinks;ikink0++){
5254 AliKink * kink0 = (AliKink*) kinks->At(indexes[ikink0]);
5255 if (!kink0) continue;
5257 for (Int_t ikink1=0;ikink1<ikink0;ikink1++){
5258 kink0 = (AliKink*) kinks->At(indexes[ikink0]);
5259 if (!kink0) continue;
5260 AliKink * kink1 = (AliKink*) kinks->At(indexes[ikink1]);
5261 if (!kink1) continue;
5262 // if not close kink continue
5263 if (TMath::Abs(kink1->GetPosition()[2]-kink0->GetPosition()[2])>10) continue;
5264 if (TMath::Abs(kink1->GetPosition()[1]-kink0->GetPosition()[1])>10) continue;
5265 if (TMath::Abs(kink1->GetPosition()[0]-kink0->GetPosition()[0])>10) continue;
5267 AliTPCseed &mother0 = mothers[indexes[ikink0]];
5268 AliTPCseed &daughter0 = daughters[indexes[ikink0]];
5269 AliTPCseed &mother1 = mothers[indexes[ikink1]];
5270 AliTPCseed &daughter1 = daughters[indexes[ikink1]];
5271 Int_t row0 = (kink0->GetTPCRow0()+kink1->GetTPCRow0())/2;
5280 for (Int_t i=0;i<row0;i++){
5281 if (mother0.GetClusterIndex(i)>0 && mother1.GetClusterIndex(i)>0){
5284 if (mother0.GetClusterIndex(i)==mother1.GetClusterIndex(i)){
5291 for (Int_t i=row0;i<158;i++){
5292 if (daughter0.GetClusterIndex(i)>0 && daughter0.GetClusterIndex(i)>0){
5295 if (mother0.GetClusterIndex(i)==mother1.GetClusterIndex(i)){
5301 Float_t ratio = Float_t(same+1)/Float_t(both+1);
5302 Float_t ratiom = Float_t(samem+1)/Float_t(bothm+1);
5303 Float_t ratiod = Float_t(samed+1)/Float_t(bothd+1);
5304 if (ratio>0.3 && ratiom>0.5 &&ratiod>0.5) {
5305 Int_t sum0 = mother0.GetNumberOfClusters()+daughter0.GetNumberOfClusters();
5306 Int_t sum1 = mother1.GetNumberOfClusters()+daughter1.GetNumberOfClusters();
5308 shared[kink0->GetIndex(0)]= kTRUE;
5309 shared[kink0->GetIndex(1)]= kTRUE;
5310 delete kinks->RemoveAt(indexes[ikink0]);
5314 shared[kink1->GetIndex(0)]= kTRUE;
5315 shared[kink1->GetIndex(1)]= kTRUE;
5316 delete kinks->RemoveAt(indexes[ikink1]);
5323 for (Int_t i=0;i<nkinks;i++){
5324 AliKink * kinkl = (AliKink*) kinks->At(indexes[i]);
5325 if (!kinkl) continue;
5326 kinkl->SetTPCRow0(GetRowNumber(kinkl->GetR()));
5327 Int_t index0 = kinkl->GetIndex(0);
5328 Int_t index1 = kinkl->GetIndex(1);
5329 if (circular[index0]||(circular[index1]&&kinkl->GetDistance()>0.2)) continue;
5330 kinkl->SetMultiple(usage[index0],0);
5331 kinkl->SetMultiple(usage[index1],1);
5332 if (kinkl->GetMultiple()[0]+kinkl->GetMultiple()[1]>2) continue;
5333 if (kinkl->GetMultiple()[0]+kinkl->GetMultiple()[1]>0 && quality[indexes[i]]>0.2) continue;
5334 if (kinkl->GetMultiple()[0]+kinkl->GetMultiple()[1]>0 && kinkl->GetDistance()>0.2) continue;
5335 if (circular[index0]||(circular[index1]&&kinkl->GetDistance()>0.1)) continue;
5337 AliTPCseed * ktrack0 = (AliTPCseed*)array->At(index0);
5338 AliTPCseed * ktrack1 = (AliTPCseed*)array->At(index1);
5339 if (!ktrack0 || !ktrack1) continue;
5340 Int_t index = esd->AddKink(kinkl);
5343 if ( ktrack0->GetKinkIndex(0)==0 && ktrack1->GetKinkIndex(0)==0) { //best kink
5344 if (mothers[indexes[i]].GetNumberOfClusters()>20 && daughters[indexes[i]].GetNumberOfClusters()>20 && (mothers[indexes[i]].GetNumberOfClusters()+daughters[indexes[i]].GetNumberOfClusters())>100){
5345 *ktrack0 = mothers[indexes[i]];
5346 *ktrack1 = daughters[indexes[i]];
5350 ktrack0->SetKinkIndex(usage[index0],-(index+1));
5351 ktrack1->SetKinkIndex(usage[index1], (index+1));
5356 // Remove tracks corresponding to shared kink's
5358 for (Int_t i=0;i<nentries;i++){
5359 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5360 if (!track0) continue;
5361 if (track0->GetKinkIndex(0)!=0) continue;
5362 if (shared[i]) delete array->RemoveAt(i);
5367 RemoveUsed2(array,0.5,0.4,30);
5369 for (Int_t i=0;i<nentries;i++){
5370 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5371 if (!track0) continue;
5372 track0->CookdEdx(0.02,0.6);
5376 for (Int_t i=0;i<nentries;i++){
5377 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5378 if (!track0) continue;
5379 if (track0->Pt()<1.4) continue;
5380 //remove double high momenta tracks - overlapped with kink candidates
5383 for (Int_t icl=track0->GetFirstPoint();icl<track0->GetLastPoint(); icl++){
5384 if (track0->GetClusterPointer(icl)!=0){
5386 if (track0->GetClusterPointer(icl)->IsUsed(10)) ishared++;
5389 if (Float_t(ishared+1)/Float_t(all+1)>0.5) {
5390 delete array->RemoveAt(i);
5394 if (track0->GetKinkIndex(0)!=0) continue;
5395 if (track0->GetNumberOfClusters()<80) continue;
5397 AliTPCseed *pmother = new AliTPCseed();
5398 AliTPCseed *pdaughter = new AliTPCseed();
5399 AliKink *pkink = new AliKink;
5401 AliTPCseed & mother = *pmother;
5402 AliTPCseed & daughter = *pdaughter;
5403 AliKink & kinkl = *pkink;
5404 if (CheckKinkPoint(track0,mother,daughter, kinkl)){
5405 if (mother.GetNumberOfClusters()<30||daughter.GetNumberOfClusters()<20) {
5409 continue; //too short tracks
5411 if (mother.Pt()<1.4) {
5417 Int_t row0= kinkl.GetTPCRow0();
5418 if (kinkl.GetDistance()>0.5 || kinkl.GetR()<110. || kinkl.GetR()>240.) {
5425 Int_t index = esd->AddKink(&kinkl);
5426 mother.SetKinkIndex(0,-(index+1));
5427 daughter.SetKinkIndex(0,index+1);
5428 if (mother.GetNumberOfClusters()>50) {
5429 delete array->RemoveAt(i);
5430 array->AddAt(new AliTPCseed(mother),i);
5433 array->AddLast(new AliTPCseed(mother));
5435 array->AddLast(new AliTPCseed(daughter));
5436 for (Int_t icl=0;icl<row0;icl++) {
5437 if (mother.GetClusterPointer(icl)) mother.GetClusterPointer(icl)->Use(20);
5440 for (Int_t icl=row0;icl<158;icl++) {
5441 if (daughter.GetClusterPointer(icl)) daughter.GetClusterPointer(icl)->Use(20);
5450 delete [] daughters;
5472 AliInfo(Form("Ncandidates=\t%d\t%d\t%d\t%d\n",esd->GetNumberOfKinks(),ncandidates,ntracks,nall));
5477 Int_t AliTPCtrackerMI::RefitKink(AliTPCseed &mother, AliTPCseed &daughter, const AliESDkink &knk)
5480 // refit kink towards to the vertex
5483 AliKink &kink=(AliKink &)knk;
5485 Int_t row0 = GetRowNumber(kink.GetR());
5486 FollowProlongation(mother,0);
5487 mother.Reset(kFALSE);
5489 FollowProlongation(daughter,row0);
5490 daughter.Reset(kFALSE);
5491 FollowBackProlongation(daughter,158);
5492 daughter.Reset(kFALSE);
5493 Int_t first = TMath::Max(row0-20,30);
5494 Int_t last = TMath::Min(row0+20,140);
5496 const Int_t kNdiv =5;
5497 AliTPCseed param0[kNdiv]; // parameters along the track
5498 AliTPCseed param1[kNdiv]; // parameters along the track
5499 AliKink kinks[kNdiv]; // corresponding kink parameters
5502 for (Int_t irow=0; irow<kNdiv;irow++){
5503 rows[irow] = first +((last-first)*irow)/(kNdiv-1);
5505 // store parameters along the track
5507 for (Int_t irow=0;irow<kNdiv;irow++){
5508 FollowBackProlongation(mother, rows[irow]);
5509 FollowProlongation(daughter,rows[kNdiv-1-irow]);
5510 param0[irow] = mother;
5511 param1[kNdiv-1-irow] = daughter;
5515 for (Int_t irow=0; irow<kNdiv-1;irow++){
5516 if (param0[irow].GetNumberOfClusters()<kNdiv||param1[irow].GetNumberOfClusters()<kNdiv) continue;
5517 kinks[irow].SetMother(param0[irow]);
5518 kinks[irow].SetDaughter(param1[irow]);
5519 kinks[irow].Update();
5522 // choose kink with best "quality"
5524 Double_t mindist = 10000;
5525 for (Int_t irow=0;irow<kNdiv;irow++){
5526 if (param0[irow].GetNumberOfClusters()<20||param1[irow].GetNumberOfClusters()<20) continue;
5527 if (TMath::Abs(kinks[irow].GetR())>240.) continue;
5528 if (TMath::Abs(kinks[irow].GetR())<100.) continue;
5530 Float_t normdist = TMath::Abs(param0[irow].GetX()-kinks[irow].GetR())*(0.1+kink.GetDistance());
5531 normdist/= (param0[irow].GetNumberOfClusters()+param1[irow].GetNumberOfClusters()+40.);
5532 if (normdist < mindist){
5538 if (index==-1) return 0;
5541 param0[index].Reset(kTRUE);
5542 FollowProlongation(param0[index],0);
5544 mother = param0[index];
5545 daughter = param1[index]; // daughter in vertex
5547 kink.SetMother(mother);
5548 kink.SetDaughter(daughter);
5550 kink.SetTPCRow0(GetRowNumber(kink.GetR()));
5551 kink.SetTPCncls(param0[index].GetNumberOfClusters(),0);
5552 kink.SetTPCncls(param1[index].GetNumberOfClusters(),1);
5553 kink.SetLabel(CookLabel(&mother,0.4, 0,kink.GetTPCRow0()),0);
5554 kink.SetLabel(CookLabel(&daughter,0.4, kink.GetTPCRow0(),160),1);
5555 mother.SetLabel(kink.GetLabel(0));
5556 daughter.SetLabel(kink.GetLabel(1));
5562 void AliTPCtrackerMI::UpdateKinkQualityM(AliTPCseed * seed){
5564 // update Kink quality information for mother after back propagation
5566 if (seed->GetKinkIndex(0)>=0) return;
5567 for (Int_t ikink=0;ikink<3;ikink++){
5568 Int_t index = seed->GetKinkIndex(ikink);
5569 if (index>=0) break;
5570 index = TMath::Abs(index)-1;
5571 AliESDkink * kink = fEvent->GetKink(index);
5572 kink->SetTPCDensity(-1,0,0);
5573 kink->SetTPCDensity(1,0,1);
5575 Int_t row0 = kink->GetTPCRow0() - 2 - Int_t( 0.5/ (0.05+kink->GetAngle(2)));
5576 if (row0<15) row0=15;
5578 Int_t row1 = kink->GetTPCRow0() + 2 + Int_t( 0.5/ (0.05+kink->GetAngle(2)));
5579 if (row1>145) row1=145;
5581 Int_t found,foundable,shared;
5582 seed->GetClusterStatistic(0,row0, found, foundable,shared,kFALSE);
5583 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),0,0);
5584 seed->GetClusterStatistic(row1,155, found, foundable,shared,kFALSE);
5585 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),0,1);
5590 void AliTPCtrackerMI::UpdateKinkQualityD(AliTPCseed * seed){
5592 // update Kink quality information for daughter after refit
5594 if (seed->GetKinkIndex(0)<=0) return;
5595 for (Int_t ikink=0;ikink<3;ikink++){
5596 Int_t index = seed->GetKinkIndex(ikink);
5597 if (index<=0) break;
5598 index = TMath::Abs(index)-1;
5599 AliESDkink * kink = fEvent->GetKink(index);
5600 kink->SetTPCDensity(-1,1,0);
5601 kink->SetTPCDensity(-1,1,1);
5603 Int_t row0 = kink->GetTPCRow0() -2 - Int_t( 0.5/ (0.05+kink->GetAngle(2)));
5604 if (row0<15) row0=15;
5606 Int_t row1 = kink->GetTPCRow0() +2 + Int_t( 0.5/ (0.05+kink->GetAngle(2)));
5607 if (row1>145) row1=145;
5609 Int_t found,foundable,shared;
5610 seed->GetClusterStatistic(0,row0, found, foundable,shared,kFALSE);
5611 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),1,0);
5612 seed->GetClusterStatistic(row1,155, found, foundable,shared,kFALSE);
5613 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),1,1);
5619 Int_t AliTPCtrackerMI::CheckKinkPoint(AliTPCseed*seed,AliTPCseed &mother, AliTPCseed &daughter, const AliESDkink &knk)
5622 // check kink point for given track
5623 // if return value=0 kink point not found
5624 // otherwise seed0 correspond to mother particle
5625 // seed1 correspond to daughter particle
5626 // kink parameter of kink point
5627 AliKink &kink=(AliKink &)knk;
5629 Int_t middlerow = (seed->GetFirstPoint()+seed->GetLastPoint())/2;
5630 Int_t first = seed->GetFirstPoint();
5631 Int_t last = seed->GetLastPoint();
5632 if (last-first<20) return 0; // shortest length - 2*30 = 60 pad-rows
5635 AliTPCseed *seed1 = ReSeed(seed,middlerow+20, kTRUE); //middle of chamber
5636 if (!seed1) return 0;
5637 FollowProlongation(*seed1,seed->GetLastPoint()-20);
5638 seed1->Reset(kTRUE);
5639 FollowProlongation(*seed1,158);
5640 seed1->Reset(kTRUE);
5641 last = seed1->GetLastPoint();
5643 AliTPCseed *seed0 = new AliTPCseed(*seed);
5644 seed0->Reset(kFALSE);
5647 AliTPCseed param0[20]; // parameters along the track
5648 AliTPCseed param1[20]; // parameters along the track
5649 AliKink kinks[20]; // corresponding kink parameters
5651 for (Int_t irow=0; irow<20;irow++){
5652 rows[irow] = first +((last-first)*irow)/19;
5654 // store parameters along the track
5656 for (Int_t irow=0;irow<20;irow++){
5657 FollowBackProlongation(*seed0, rows[irow]);
5658 FollowProlongation(*seed1,rows[19-irow]);
5659 param0[irow] = *seed0;
5660 param1[19-irow] = *seed1;
5664 for (Int_t irow=0; irow<19;irow++){
5665 kinks[irow].SetMother(param0[irow]);
5666 kinks[irow].SetDaughter(param1[irow]);
5667 kinks[irow].Update();
5670 // choose kink with biggest change of angle
5672 Double_t maxchange= 0;
5673 for (Int_t irow=1;irow<19;irow++){
5674 if (TMath::Abs(kinks[irow].GetR())>240.) continue;
5675 if (TMath::Abs(kinks[irow].GetR())<110.) continue;
5676 Float_t quality = TMath::Abs(kinks[irow].GetAngle(2))/(3.+TMath::Abs(kinks[irow].GetR()-param0[irow].GetX()));
5677 if ( quality > maxchange){
5678 maxchange = quality;
5685 if (index<0) return 0;
5687 Int_t row0 = GetRowNumber(kinks[index].GetR()); //row 0 estimate
5688 seed0 = new AliTPCseed(param0[index]);
5689 seed1 = new AliTPCseed(param1[index]);
5690 seed0->Reset(kFALSE);
5691 seed1->Reset(kFALSE);
5692 seed0->ResetCovariance(10.);
5693 seed1->ResetCovariance(10.);
5694 FollowProlongation(*seed0,0);
5695 FollowBackProlongation(*seed1,158);
5696 mother = *seed0; // backup mother at position 0
5697 seed0->Reset(kFALSE);
5698 seed1->Reset(kFALSE);
5699 seed0->ResetCovariance(10.);
5700 seed1->ResetCovariance(10.);
5702 first = TMath::Max(row0-20,0);
5703 last = TMath::Min(row0+20,158);
5705 for (Int_t irow=0; irow<20;irow++){
5706 rows[irow] = first +((last-first)*irow)/19;
5708 // store parameters along the track
5710 for (Int_t irow=0;irow<20;irow++){
5711 FollowBackProlongation(*seed0, rows[irow]);
5712 FollowProlongation(*seed1,rows[19-irow]);
5713 param0[irow] = *seed0;
5714 param1[19-irow] = *seed1;
5718 for (Int_t irow=0; irow<19;irow++){
5719 kinks[irow].SetMother(param0[irow]);
5720 kinks[irow].SetDaughter(param1[irow]);
5721 // param0[irow].Dump();
5722 //param1[irow].Dump();
5723 kinks[irow].Update();
5726 // choose kink with biggest change of angle
5729 for (Int_t irow=0;irow<20;irow++){
5730 if (TMath::Abs(kinks[irow].GetR())>250.) continue;
5731 if (TMath::Abs(kinks[irow].GetR())<90.) continue;
5732 Float_t quality = TMath::Abs(kinks[irow].GetAngle(2))/(3.+TMath::Abs(kinks[irow].GetR()-param0[irow].GetX()));
5733 if ( quality > maxchange){
5734 maxchange = quality;
5741 if (index==-1 || param0[index].GetNumberOfClusters()+param1[index].GetNumberOfClusters()<100){
5747 // Float_t anglesigma = TMath::Sqrt(param0[index].fC22+param0[index].fC33+param1[index].fC22+param1[index].fC33);
5749 kink.SetMother(param0[index]);
5750 kink.SetDaughter(param1[index]);
5753 Double_t chi2P2 = param0[index].GetParameter()[2]-param1[index].GetParameter()[2];
5755 chi2P2/=param0[index].GetCovariance()[5]+param1[index].GetCovariance()[5];
5756 Double_t chi2P3 = param0[index].GetParameter()[3]-param1[index].GetParameter()[3];
5758 chi2P3/=param0[index].GetCovariance()[9]+param1[index].GetCovariance()[9];
5760 if (AliTPCReconstructor::StreamLevel()>1) {
5761 (*fDebugStreamer)<<"kinkHpt"<<
5764 "p0.="<<¶m0[index]<<
5765 "p1.="<<¶m1[index]<<
5769 if ( chi2P2+chi2P3<AliTPCReconstructor::GetRecoParam()->GetKinkAngleCutChi2(0)){
5776 row0 = GetRowNumber(kink.GetR());
5777 kink.SetTPCRow0(row0);
5778 kink.SetLabel(CookLabel(seed0,0.5,0,row0),0);
5779 kink.SetLabel(CookLabel(seed1,0.5,row0,158),1);
5780 kink.SetIndex(-10,0);
5781 kink.SetIndex(int(param0[index].GetNumberOfClusters()+param1[index].GetNumberOfClusters()),1);
5782 kink.SetTPCncls(param0[index].GetNumberOfClusters(),0);
5783 kink.SetTPCncls(param1[index].GetNumberOfClusters(),1);
5786 // new (&mother) AliTPCseed(param0[index]);
5787 daughter = param1[index];
5788 daughter.SetLabel(kink.GetLabel(1));
5789 param0[index].Reset(kTRUE);
5790 FollowProlongation(param0[index],0);
5791 mother = param0[index];
5792 mother.SetLabel(kink.GetLabel(0));
5793 if ( chi2P2+chi2P3<AliTPCReconstructor::GetRecoParam()->GetKinkAngleCutChi2(1)){
5805 AliTPCseed* AliTPCtrackerMI::ReSeed(AliTPCseed *t)
5808 // reseed - refit - track
5811 // Int_t last = fSectors->GetNRows()-1;
5813 if (fSectors == fOuterSec){
5814 first = TMath::Max(first, t->GetFirstPoint()-fInnerSec->GetNRows());
5818 first = t->GetFirstPoint();
5820 AliTPCseed * seed = MakeSeed(t,0.1,0.5,0.9);
5821 FollowBackProlongation(*t,fSectors->GetNRows()-1);
5823 FollowProlongation(*t,first);
5833 //_____________________________________________________________________________
5834 Int_t AliTPCtrackerMI::ReadSeeds(const TFile *inp) {
5835 //-----------------------------------------------------------------
5836 // This function reades track seeds.
5837 //-----------------------------------------------------------------
5838 TDirectory *savedir=gDirectory;
5840 TFile *in=(TFile*)inp;
5841 if (!in->IsOpen()) {
5842 cerr<<"AliTPCtrackerMI::ReadSeeds(): input file is not open !\n";
5847 TTree *seedTree=(TTree*)in->Get("Seeds");
5849 cerr<<"AliTPCtrackerMI::ReadSeeds(): ";
5850 cerr<<"can't get a tree with track seeds !\n";
5853 AliTPCtrack *seed=new AliTPCtrack;
5854 seedTree->SetBranchAddress("tracks",&seed);
5856 if (fSeeds==0) fSeeds=new TObjArray(15000);
5858 Int_t n=(Int_t)seedTree->GetEntries();
5859 for (Int_t i=0; i<n; i++) {
5860 seedTree->GetEvent(i);
5861 fSeeds->AddLast(new AliTPCseed(*seed/*,seed->GetAlpha()*/));
5870 Int_t AliTPCtrackerMI::Clusters2Tracks (AliESDEvent *const esd)
5873 // clusters to tracks
5875 if (fSeeds) DeleteSeeds();
5878 AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ;
5879 transform->SetCurrentTimeStamp( esd->GetTimeStamp());
5880 transform->SetCurrentRun(esd->GetRunNumber());
5883 if (!fSeeds) return 1;
5885 if (AliTPCReconstructor::StreamLevel()>3) DumpClusters(0,fSeeds);
5891 //_____________________________________________________________________________
5892 Int_t AliTPCtrackerMI::Clusters2Tracks() {
5893 //-----------------------------------------------------------------
5894 // This is a track finder.
5895 //-----------------------------------------------------------------
5896 TDirectory *savedir=gDirectory;
5900 fSeeds = Tracking();
5903 Info("Clusters2Tracks","Time for tracking: \t");timer.Print();timer.Start();
5905 //activate again some tracks
5906 for (Int_t i=0; i<fSeeds->GetEntriesFast(); i++) {
5907 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
5909 Int_t nc=t.GetNumberOfClusters();
5911 delete fSeeds->RemoveAt(i);
5915 if (pt->GetRemoval()==10) {
5916 if (pt->GetDensityFirst(20)>0.8 || pt->GetDensityFirst(30)>0.8 || pt->GetDensityFirst(40)>0.7)
5917 pt->Desactivate(10); // make track again active
5919 pt->Desactivate(20);
5920 delete fSeeds->RemoveAt(i);
5925 RemoveUsed2(fSeeds,0.85,0.85,0);
5926 if (AliTPCReconstructor::GetRecoParam()->GetDoKinks()) FindKinks(fSeeds,fEvent);
5927 //FindCurling(fSeeds, fEvent,0);
5928 if (AliTPCReconstructor::StreamLevel()>5) FindMultiMC(fSeeds, fEvent,-1); // find multi found tracks
5929 RemoveUsed2(fSeeds,0.5,0.4,20);
5930 FindSplitted(fSeeds, fEvent,0); // find multi found tracks
5931 if (AliTPCReconstructor::StreamLevel()>5) FindMultiMC(fSeeds, fEvent,0); // find multi found tracks
5934 // // refit short tracks
5936 Int_t nseed=fSeeds->GetEntriesFast();
5939 for (Int_t i=0; i<nseed; i++) {
5940 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
5942 Int_t nc=t.GetNumberOfClusters();
5944 delete fSeeds->RemoveAt(i);
5947 CookLabel(pt,0.1); //For comparison only
5948 //if ((pt->IsActive() || (pt->fRemoval==10) )&& nc>50 &&pt->GetNumberOfClusters()>0.4*pt->fNFoundable){
5949 if ((pt->IsActive() || (pt->GetRemoval()==10) )){
5951 if (fDebug>0) cerr<<found<<'\r';
5955 delete fSeeds->RemoveAt(i);
5959 //RemoveOverlap(fSeeds,0.99,7,kTRUE);
5961 //RemoveUsed(fSeeds,0.9,0.9,6);
5963 nseed=fSeeds->GetEntriesFast();
5965 for (Int_t i=0; i<nseed; i++) {
5966 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
5968 Int_t nc=t.GetNumberOfClusters();
5970 delete fSeeds->RemoveAt(i);
5974 t.CookdEdx(0.02,0.6);
5975 // CheckKinkPoint(&t,0.05);
5976 //if ((pt->IsActive() || (pt->fRemoval==10) )&& nc>50 &&pt->GetNumberOfClusters()>0.4*pt->fNFoundable){
5977 if ((pt->IsActive() || (pt->GetRemoval()==10) )){
5985 delete fSeeds->RemoveAt(i);
5986 //AliTPCseed * seed1 = ReSeed(pt,0.05,0.5,1);
5988 // FollowProlongation(*seed1,0);
5989 // Int_t n = seed1->GetNumberOfClusters();
5990 // printf("fP4\t%f\t%f\n",seed1->GetC(),pt->GetC());
5991 // printf("fN\t%d\t%d\n", seed1->GetNumberOfClusters(),pt->GetNumberOfClusters());
5994 //AliTPCseed * seed2 = ReSeed(pt,0.95,0.5,0.05);
5998 SortTracks(fSeeds, 1);
6002 PrepareForBackProlongation(fSeeds,5.);
6003 PropagateBack(fSeeds);
6004 printf("Time for back propagation: \t");timer.Print();timer.Start();
6008 PrepareForProlongation(fSeeds,5.);
6009 PropagateForward2(fSeeds);
6011 printf("Time for FORWARD propagation: \t");timer.Print();timer.Start();
6012 // RemoveUsed(fSeeds,0.7,0.7,6);
6013 //RemoveOverlap(fSeeds,0.9,7,kTRUE);
6015 nseed=fSeeds->GetEntriesFast();
6017 for (Int_t i=0; i<nseed; i++) {
6018 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
6020 Int_t nc=t.GetNumberOfClusters();
6022 delete fSeeds->RemoveAt(i);
6025 t.CookdEdx(0.02,0.6);
6026 // CookLabel(pt,0.1); //For comparison only
6027 //if ((pt->IsActive() || (pt->fRemoval==10) )&& nc>50 &&pt->GetNumberOfClusters()>0.4*pt->fNFoundable){
6028 if ((pt->IsActive() || (pt->fRemoval==10) )){
6029 cerr<<found++<<'\r';
6032 delete fSeeds->RemoveAt(i);
6037 // fNTracks = found;
6039 Info("Clusters2Tracks","Time for overlap removal, track writing and dedx cooking: \t"); timer.Print();timer.Start();
6042 // cerr<<"Number of found tracks : "<<"\t"<<found<<endl;
6043 Info("Clusters2Tracks","Number of found tracks %d",found);
6045 // UnloadClusters();
6050 void AliTPCtrackerMI::Tracking(TObjArray * arr)
6053 // tracking of the seeds
6056 fSectors = fOuterSec;
6057 ParallelTracking(arr,150,63);
6058 fSectors = fOuterSec;
6059 ParallelTracking(arr,63,0);
6062 TObjArray * AliTPCtrackerMI::Tracking(Int_t seedtype, Int_t i1, Int_t i2, Float_t cuts[4], Float_t dy, Int_t dsec)
6067 TObjArray * arr = new TObjArray;
6069 fSectors = fOuterSec;
6072 for (Int_t sec=0;sec<fkNOS;sec++){
6073 if (seedtype==3) MakeSeeds3(arr,sec,i1,i2,cuts,dy, dsec);
6074 if (seedtype==4) MakeSeeds5(arr,sec,i1,i2,cuts,dy);
6075 if (seedtype==2) MakeSeeds2(arr,sec,i1,i2,cuts,dy);
6078 Info("Tracking","\nSeeding - %d\t%d\t%d\t%d\n",seedtype,i1,i2,arr->GetEntriesFast());
6090 TObjArray * AliTPCtrackerMI::Tracking()
6094 if (AliTPCReconstructor::GetRecoParam()->GetSpecialSeeding()) return TrackingSpecial();
6097 Int_t nup=fOuterSec->GetNRows()+fInnerSec->GetNRows();
6099 TObjArray * seeds = new TObjArray;
6101 Int_t fLastSeedRowSec=AliTPCReconstructor::GetRecoParam()->GetLastSeedRowSec();
6102 Int_t gapPrim = AliTPCReconstructor::GetRecoParam()->GetSeedGapPrim();
6103 Int_t gapSec = AliTPCReconstructor::GetRecoParam()->GetSeedGapSec();
6111 Float_t fnumber = 3.0;
6112 Float_t fdensity = 3.0;
6117 for (Int_t delta = 0; delta<18; delta+=gapPrim){
6121 arr = Tracking(3,nup-1-delta,nup-1-delta-gap,cuts,-1,1);
6122 SumTracks(seeds,arr);
6123 SignClusters(seeds,fnumber,fdensity);
6125 for (Int_t i=2;i<6;i+=2){
6126 // seed high pt tracks
6129 arr = Tracking(3,nup-i-delta,nup-i-delta-gap,cuts,-1,0);
6130 SumTracks(seeds,arr);
6131 SignClusters(seeds,fnumber,fdensity);
6136 // RemoveUsed(seeds,0.9,0.9,1);
6137 // UnsignClusters();
6138 // SignClusters(seeds,fnumber,fdensity);
6142 for (Int_t delta = 20; delta<120; delta+=gapPrim){
6144 // seed high pt tracks
6148 arr = Tracking(3,nup-delta,nup-delta-gap,cuts,-1);
6149 SumTracks(seeds,arr);
6150 SignClusters(seeds,fnumber,fdensity);
6155 arr = Tracking(3,nup-delta-5,nup-delta-5-gap,cuts,-1);
6156 SumTracks(seeds,arr);
6157 SignClusters(seeds,fnumber,fdensity);
6168 Info("Tracking()","\n\nPrimary seeding\t%d\n\n",seeds->GetEntriesFast());
6172 // RemoveUsed(seeds,0.75,0.75,1);
6174 //SignClusters(seeds,fnumber,fdensity);
6183 arr = Tracking(4,nup-1,nup-1-gap,cuts,-1);
6184 SumTracks(seeds,arr);
6185 SignClusters(seeds,fnumber,fdensity);
6187 arr = Tracking(4,nup-2,nup-2-gap,cuts,-1);
6188 SumTracks(seeds,arr);
6189 SignClusters(seeds,fnumber,fdensity);
6191 arr = Tracking(4,nup-3,nup-3-gap,cuts,-1);
6192 SumTracks(seeds,arr);
6193 SignClusters(seeds,fnumber,fdensity);
6195 arr = Tracking(4,nup-5,nup-5-gap,cuts,-1);
6196 SumTracks(seeds,arr);
6197 SignClusters(seeds,fnumber,fdensity);
6199 arr = Tracking(4,nup-7,nup-7-gap,cuts,-1);
6200 SumTracks(seeds,arr);
6201 SignClusters(seeds,fnumber,fdensity);
6204 arr = Tracking(4,nup-9,nup-9-gap,cuts,-1);
6205 SumTracks(seeds,arr);
6206 SignClusters(seeds,fnumber,fdensity);
6210 for (Int_t delta = 9; delta<30; delta+=gapSec){
6216 arr = Tracking(4,nup-1-delta,nup-1-delta-gap,cuts,-1);
6217 SumTracks(seeds,arr);
6218 SignClusters(seeds,fnumber,fdensity);
6220 arr = Tracking(4,nup-3-delta,nup-5-delta-gap,cuts,4);
6221 SumTracks(seeds,arr);
6222 SignClusters(seeds,fnumber,fdensity);
6235 for (Int_t delta = 30; delta<fLastSeedRowSec; delta+=gapSec){
6241 arr = Tracking(4,nup-1-delta,nup-1-delta-gap,cuts,-1);
6242 SumTracks(seeds,arr);
6243 SignClusters(seeds,fnumber,fdensity);
6245 arr = Tracking(4,nup-5-delta,nup-5-delta-gap,cuts,5 );
6246 SumTracks(seeds,arr);
6247 SignClusters(seeds,fnumber,fdensity);
6251 Info("Tracking()","\n\nSecondary seeding\t%d\n\n",seeds->GetEntriesFast());
6262 TObjArray * AliTPCtrackerMI::TrackingSpecial()
6265 // seeding adjusted for laser and cosmic tests - short tracks with big inclination angle
6266 // no primary vertex seeding tried
6270 Int_t nup=fOuterSec->GetNRows()+fInnerSec->GetNRows();
6272 TObjArray * seeds = new TObjArray;
6277 Float_t fnumber = 3.0;
6278 Float_t fdensity = 3.0;
6281 cuts[0] = AliTPCReconstructor::GetRecoParam()->GetMaxC(); // max curvature
6282 cuts[1] = 3.5; // max tan(phi) angle for seeding
6283 cuts[2] = 3.; // not used (cut on z primary vertex)
6284 cuts[3] = 3.5; // max tan(theta) angle for seeding
6286 for (Int_t delta = 0; nup-delta-gap-1>0; delta+=3){
6288 arr = Tracking(4,nup-1-delta,nup-1-delta-gap,cuts,-1);
6289 SumTracks(seeds,arr);
6290 SignClusters(seeds,fnumber,fdensity);
6294 Info("Tracking()","\n\nSecondary seeding\t%d\n\n",seeds->GetEntriesFast());
6305 void AliTPCtrackerMI::SumTracks(TObjArray *arr1,TObjArray *&arr2) const
6308 //sum tracks to common container
6309 //remove suspicious tracks
6310 Int_t nseed = arr2->GetEntriesFast();
6311 for (Int_t i=0;i<nseed;i++){
6312 AliTPCseed *pt=(AliTPCseed*)arr2->UncheckedAt(i);
6315 // remove tracks with too big curvature
6317 if (TMath::Abs(pt->GetC())>AliTPCReconstructor::GetRecoParam()->GetMaxC()){
6318 delete arr2->RemoveAt(i);
6321 // REMOVE VERY SHORT TRACKS
6322 if (pt->GetNumberOfClusters()<20){
6323 delete arr2->RemoveAt(i);
6326 // NORMAL ACTIVE TRACK
6327 if (pt->IsActive()){
6328 arr1->AddLast(arr2->RemoveAt(i));
6331 //remove not usable tracks
6332 if (pt->GetRemoval()!=10){
6333 delete arr2->RemoveAt(i);
6337 // ENABLE ONLY ENOUGH GOOD STOPPED TRACKS
6338 if (pt->GetDensityFirst(20)>0.8 || pt->GetDensityFirst(30)>0.8 || pt->GetDensityFirst(40)>0.7)
6339 arr1->AddLast(arr2->RemoveAt(i));
6341 delete arr2->RemoveAt(i);
6345 delete arr2; arr2 = 0;
6350 void AliTPCtrackerMI::ParallelTracking(TObjArray *const arr, Int_t rfirst, Int_t rlast)
6353 // try to track in parralel
6355 Int_t nseed=arr->GetEntriesFast();
6356 //prepare seeds for tracking
6357 for (Int_t i=0; i<nseed; i++) {
6358 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i), &t=*pt;
6360 if (!t.IsActive()) continue;
6361 // follow prolongation to the first layer
6362 if ( (fSectors ==fInnerSec) || (t.GetFirstPoint()-fkParam->GetNRowLow()>rfirst+1) )
6363 FollowProlongation(t, rfirst+1);
6368 for (Int_t nr=rfirst; nr>=rlast; nr--){
6369 if (nr<fInnerSec->GetNRows())
6370 fSectors = fInnerSec;
6372 fSectors = fOuterSec;
6373 // make indexes with the cluster tracks for given
6375 // find nearest cluster
6376 for (Int_t i=0; i<nseed; i++) {
6377 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i), &t=*pt;
6379 if (nr==80) pt->UpdateReference();
6380 if (!pt->IsActive()) continue;
6381 // if ( (fSectors ==fOuterSec) && (pt->fFirstPoint-fkParam->GetNRowLow())<nr) continue;
6382 if (pt->GetRelativeSector()>17) {
6385 UpdateClusters(t,nr);
6387 // prolonagate to the nearest cluster - if founded
6388 for (Int_t i=0; i<nseed; i++) {
6389 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
6391 if (!pt->IsActive()) continue;
6392 // if ((fSectors ==fOuterSec) && (pt->fFirstPoint-fkParam->GetNRowLow())<nr) continue;
6393 if (pt->GetRelativeSector()>17) {
6396 FollowToNextCluster(*pt,nr);
6401 void AliTPCtrackerMI::PrepareForBackProlongation(const TObjArray *const arr,Float_t fac) const
6405 // if we use TPC track itself we have to "update" covariance
6407 Int_t nseed= arr->GetEntriesFast();
6408 for (Int_t i=0;i<nseed;i++){
6409 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
6413 //rotate to current local system at first accepted point
6414 Int_t index = pt->GetClusterIndex2(pt->GetFirstPoint());
6415 Int_t sec = (index&0xff000000)>>24;
6417 Float_t angle1 = fInnerSec->GetAlpha()*sec+fInnerSec->GetAlphaShift();
6418 if (angle1>TMath::Pi())
6419 angle1-=2.*TMath::Pi();
6420 Float_t angle2 = pt->GetAlpha();
6422 if (TMath::Abs(angle1-angle2)>0.001){
6423 if (!pt->Rotate(angle1-angle2)) return;
6424 //angle2 = pt->GetAlpha();
6425 //pt->fRelativeSector = pt->GetAlpha()/fInnerSec->GetAlpha();
6426 //if (pt->GetAlpha()<0)
6427 // pt->fRelativeSector+=18;
6428 //sec = pt->fRelativeSector;
6437 void AliTPCtrackerMI::PrepareForProlongation(TObjArray *const arr, Float_t fac) const
6441 // if we use TPC track itself we have to "update" covariance
6443 Int_t nseed= arr->GetEntriesFast();
6444 for (Int_t i=0;i<nseed;i++){
6445 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
6448 pt->SetFirstPoint(pt->GetLastPoint());
6456 Int_t AliTPCtrackerMI::PropagateBack(const TObjArray *const arr)
6459 // make back propagation
6461 Int_t nseed= arr->GetEntriesFast();
6462 for (Int_t i=0;i<nseed;i++){
6463 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
6464 if (pt&& pt->GetKinkIndex(0)<=0) {
6465 //AliTPCseed *pt2 = new AliTPCseed(*pt);
6466 fSectors = fInnerSec;
6467 //FollowBackProlongation(*pt,fInnerSec->GetNRows()-1);
6468 //fSectors = fOuterSec;
6469 FollowBackProlongation(*pt,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1);
6470 //if (pt->GetNumberOfClusters()<(pt->fEsd->GetTPCclusters(0)) ){
6471 // Error("PropagateBack","Not prolonged track %d",pt->GetLabel());
6472 // FollowBackProlongation(*pt2,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1);
6475 if (pt&& pt->GetKinkIndex(0)>0) {
6476 AliESDkink * kink = fEvent->GetKink(pt->GetKinkIndex(0)-1);
6477 pt->SetFirstPoint(kink->GetTPCRow0());
6478 fSectors = fInnerSec;
6479 FollowBackProlongation(*pt,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1);
6487 Int_t AliTPCtrackerMI::PropagateForward2(const TObjArray *const arr)
6490 // make forward propagation
6492 Int_t nseed= arr->GetEntriesFast();
6494 for (Int_t i=0;i<nseed;i++){
6495 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
6497 FollowProlongation(*pt,0);
6506 Int_t AliTPCtrackerMI::PropagateForward()
6509 // propagate track forward
6511 Int_t nseed = fSeeds->GetEntriesFast();
6512 for (Int_t i=0;i<nseed;i++){
6513 AliTPCseed *pt = (AliTPCseed*)fSeeds->UncheckedAt(i);
6515 AliTPCseed &t = *pt;
6516 Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
6517 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
6518 if (alpha < 0. ) alpha += 2.*TMath::Pi();
6519 t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN);
6523 fSectors = fOuterSec;
6524 ParallelTracking(fSeeds,fOuterSec->GetNRows()+fInnerSec->GetNRows()-1,fInnerSec->GetNRows());
6525 fSectors = fInnerSec;
6526 ParallelTracking(fSeeds,fInnerSec->GetNRows()-1,0);
6535 Int_t AliTPCtrackerMI::PropagateBack(AliTPCseed *const pt, Int_t row0, Int_t row1)
6538 // make back propagation, in between row0 and row1
6542 fSectors = fInnerSec;
6545 if (row1<fSectors->GetNRows())
6548 r1 = fSectors->GetNRows()-1;
6550 if (row0<fSectors->GetNRows()&& r1>0 )
6551 FollowBackProlongation(*pt,r1);
6552 if (row1<=fSectors->GetNRows())
6555 r1 = row1 - fSectors->GetNRows();
6556 if (r1<=0) return 0;
6557 if (r1>=fOuterSec->GetNRows()) return 0;
6558 fSectors = fOuterSec;
6559 return FollowBackProlongation(*pt,r1);
6567 void AliTPCtrackerMI::GetShape(AliTPCseed * seed, Int_t row)
6569 // gets cluster shape
6571 AliTPCClusterParam * clparam = AliTPCcalibDB::Instance()->GetClusterParam();
6572 Float_t zdrift = TMath::Abs((fkParam->GetZLength(0)-TMath::Abs(seed->GetZ())));
6573 Int_t type = (seed->GetSector() < fkParam->GetNSector()/2) ? 0: (row>126) ? 1:2;
6574 Double_t angulary = seed->GetSnp();
6576 if (TMath::Abs(angulary)>AliTPCReconstructor::GetMaxSnpTracker()) {
6577 angulary = TMath::Sign(AliTPCReconstructor::GetMaxSnpTracker(),angulary);
6580 angulary = angulary*angulary/((1.-angulary)*(1.+angulary));
6581 Double_t angularz = seed->GetTgl()*seed->GetTgl()*(1.+angulary);
6583 Double_t sigmay = clparam->GetRMS0(0,type,zdrift,TMath::Sqrt(TMath::Abs(angulary)));
6584 Double_t sigmaz = clparam->GetRMS0(1,type,zdrift,TMath::Sqrt(TMath::Abs(angularz)));
6585 seed->SetCurrentSigmaY2(sigmay*sigmay);
6586 seed->SetCurrentSigmaZ2(sigmaz*sigmaz);
6587 // Float_t sd2 = TMath::Abs((fkParam->GetZLength(0)-TMath::Abs(seed->GetZ())))*fkParam->GetDiffL()*fkParam->GetDiffL();
6588 // // Float_t padlength = fkParam->GetPadPitchLength(seed->fSector);
6589 // Float_t padlength = GetPadPitchLength(row);
6591 // Float_t sresy = (seed->GetSector() < fkParam->GetNSector()/2) ? 0.2 :0.3;
6592 // seed->SetCurrentSigmaY2(sd2+padlength*padlength*angulary/12.+sresy*sresy);
6594 // Float_t sresz = fkParam->GetZSigma();
6595 // seed->SetCurrentSigmaZ2(sd2+padlength*padlength*angularz*angularz*(1+angulary)/12.+sresz*sresz);
6597 Float_t wy = GetSigmaY(seed);
6598 Float_t wz = GetSigmaZ(seed);
6601 if (TMath::Abs(wy/seed->fCurrentSigmaY2-1)>0.0001 || TMath::Abs(wz/seed->fCurrentSigmaZ2-1)>0.0001 ){
6602 printf("problem\n");
6609 //__________________________________________________________________________
6610 void AliTPCtrackerMI::CookLabel(AliKalmanTrack *tk, Float_t wrong) const {
6611 //--------------------------------------------------------------------
6612 //This function "cooks" a track label. If label<0, this track is fake.
6613 //--------------------------------------------------------------------
6614 AliTPCseed * t = dynamic_cast<AliTPCseed*>(tk);
6616 printf("%s:%d wrong type \n",(char*)__FILE__,__LINE__);
6620 Int_t noc=t->GetNumberOfClusters();
6622 //printf("\nnot founded prolongation\n\n\n");
6628 AliTPCclusterMI *clusters[160];
6630 for (Int_t i=0;i<160;i++) {
6637 for (i=0; i<160 && current<noc; i++) {
6639 Int_t index=t->GetClusterIndex2(i);
6640 if (index<=0) continue;
6641 if (index&0x8000) continue;
6643 //clusters[current]=GetClusterMI(index);
6644 if (t->GetClusterPointer(i)){
6645 clusters[current]=t->GetClusterPointer(i);
6651 Int_t lab=123456789;
6652 for (i=0; i<noc; i++) {
6653 AliTPCclusterMI *c=clusters[i];
6655 lab=TMath::Abs(c->GetLabel(0));
6657 for (j=0; j<noc; j++) if (lb[j]==lab || mx[j]==0) break;
6663 for (i=0; i<noc; i++) if (mx[i]>max) {max=mx[i]; lab=lb[i];}
6665 for (i=0; i<noc; i++) {
6666 AliTPCclusterMI *c=clusters[i];
6668 if (TMath::Abs(c->GetLabel(1)) == lab ||
6669 TMath::Abs(c->GetLabel(2)) == lab ) max++;
6671 if (noc<=0) { lab=-1; return;}
6672 if ((1.- Float_t(max)/(noc)) > wrong) lab=-lab;
6675 Int_t tail=Int_t(0.10*noc);
6678 for (i=1; i<160&&ind<tail; i++) {
6679 // AliTPCclusterMI *c=clusters[noc-i];
6680 AliTPCclusterMI *c=clusters[i];
6682 if (lab == TMath::Abs(c->GetLabel(0)) ||
6683 lab == TMath::Abs(c->GetLabel(1)) ||
6684 lab == TMath::Abs(c->GetLabel(2))) max++;
6687 if (max < Int_t(0.5*tail)) lab=-lab;
6694 //delete[] clusters;
6698 //__________________________________________________________________________
6699 Int_t AliTPCtrackerMI::CookLabel(AliTPCseed *const t, Float_t wrong,Int_t first, Int_t last) const {
6700 //--------------------------------------------------------------------
6701 //This function "cooks" a track label. If label<0, this track is fake.
6702 //--------------------------------------------------------------------
6703 Int_t noc=t->GetNumberOfClusters();
6705 //printf("\nnot founded prolongation\n\n\n");
6711 AliTPCclusterMI *clusters[160];
6713 for (Int_t i=0;i<160;i++) {
6720 for (i=0; i<160 && current<noc; i++) {
6721 if (i<first) continue;
6722 if (i>last) continue;
6723 Int_t index=t->GetClusterIndex2(i);
6724 if (index<=0) continue;
6725 if (index&0x8000) continue;
6727 //clusters[current]=GetClusterMI(index);
6728 if (t->GetClusterPointer(i)){
6729 clusters[current]=t->GetClusterPointer(i);
6734 //if (noc<5) return -1;
6735 Int_t lab=123456789;
6736 for (i=0; i<noc; i++) {
6737 AliTPCclusterMI *c=clusters[i];
6739 lab=TMath::Abs(c->GetLabel(0));
6741 for (j=0; j<noc; j++) if (lb[j]==lab || mx[j]==0) break;
6747 for (i=0; i<noc; i++) if (mx[i]>max) {max=mx[i]; lab=lb[i];}
6749 for (i=0; i<noc; i++) {
6750 AliTPCclusterMI *c=clusters[i];
6752 if (TMath::Abs(c->GetLabel(1)) == lab ||
6753 TMath::Abs(c->GetLabel(2)) == lab ) max++;
6755 if (noc<=0) { lab=-1; return -1;}
6756 if ((1.- Float_t(max)/(noc)) > wrong) lab=-lab;
6759 Int_t tail=Int_t(0.10*noc);
6762 for (i=1; i<160&&ind<tail; i++) {
6763 // AliTPCclusterMI *c=clusters[noc-i];
6764 AliTPCclusterMI *c=clusters[i];
6766 if (lab == TMath::Abs(c->GetLabel(0)) ||
6767 lab == TMath::Abs(c->GetLabel(1)) ||
6768 lab == TMath::Abs(c->GetLabel(2))) max++;
6771 if (max < Int_t(0.5*tail)) lab=-lab;
6774 // t->SetLabel(lab);
6778 //delete[] clusters;
6782 Int_t AliTPCtrackerMI::GetRowNumber(Double_t x[3]) const
6784 //return pad row number for given x vector
6785 Float_t phi = TMath::ATan2(x[1],x[0]);
6786 if(phi<0) phi=2.*TMath::Pi()+phi;
6787 // Get the local angle in the sector philoc
6788 const Float_t kRaddeg = 180/3.14159265358979312;
6789 Float_t phiangle = (Int_t (phi*kRaddeg/20.) + 0.5)*20./kRaddeg;
6790 Double_t localx = x[0]*TMath::Cos(phiangle)-x[1]*TMath::Sin(phiangle);
6791 return GetRowNumber(localx);
6796 void AliTPCtrackerMI::MakeBitmaps(AliTPCseed *t)
6798 //-----------------------------------------------------------------------
6799 // Fill the cluster and sharing bitmaps of the track
6800 //-----------------------------------------------------------------------
6802 Int_t firstpoint = 0;
6803 Int_t lastpoint = 159;
6804 AliTPCTrackerPoint *point;
6805 AliTPCclusterMI *cluster;
6807 for (int iter=firstpoint; iter<lastpoint; iter++) {
6808 // Change to cluster pointers to see if we have a cluster at given padrow
6809 cluster = t->GetClusterPointer(iter);
6811 t->SetClusterMapBit(iter, kTRUE);
6812 point = t->GetTrackPoint(iter);
6813 if (point->IsShared())
6814 t->SetSharedMapBit(iter,kTRUE);
6816 t->SetSharedMapBit(iter, kFALSE);
6819 t->SetClusterMapBit(iter, kFALSE);
6820 t->SetSharedMapBit(iter, kFALSE);
6825 Bool_t AliTPCtrackerMI::IsFindable(AliTPCseed & track){
6827 // return flag if there is findable cluster at given position
6830 Float_t z = track.GetZ();
6832 if (TMath::Abs(z)<(AliTPCReconstructor::GetCtgRange()*track.GetX()+kDeltaZ) &&
6833 TMath::Abs(z)<fkParam->GetZLength(0) &&
6834 (TMath::Abs(track.GetSnp())<AliTPCReconstructor::GetMaxSnpTracker()))
6840 void AliTPCtrackerMI::AddCovariance(AliTPCseed * seed){
6842 // Adding systematic error estimate to the covariance matrix
6843 // !!!! the systematic error for element 4 is in 1/cm not in pt
6845 const Double_t *param = AliTPCReconstructor::GetRecoParam()->GetSystematicError();
6847 // use only the diagonal part if not specified otherwise
6848 if (!AliTPCReconstructor::GetRecoParam()->GetUseSystematicCorrelation()) return AddCovarianceAdd(seed);
6850 Double_t *covarS= (Double_t*)seed->GetCovariance();
6851 Double_t factor[5]={1,1,1,1,1};
6852 Double_t facC = AliTracker::GetBz()*kB2C;
6853 factor[0]= TMath::Sqrt(TMath::Abs((covarS[0] + param[0]*param[0])/covarS[0]));
6854 factor[1]= TMath::Sqrt(TMath::Abs((covarS[2] + param[1]*param[1])/covarS[2]));
6855 factor[2]= TMath::Sqrt(TMath::Abs((covarS[5] + param[2]*param[2])/covarS[5]));
6856 factor[3]= TMath::Sqrt(TMath::Abs((covarS[9] + param[3]*param[3])/covarS[9]));
6857 factor[4]= TMath::Sqrt(TMath::Abs((covarS[14] + facC*facC*param[4]*param[4])/covarS[14]));
6863 for (Int_t i=0; i<5; i++){
6864 for (Int_t j=i; j<5; j++){
6865 Int_t index=seed->GetIndex(i,j);
6866 covarS[index]*=factor[i]*factor[j];
6872 void AliTPCtrackerMI::AddCovarianceAdd(AliTPCseed * seed){
6874 // Adding systematic error - as additive factor without correlation
6876 const Double_t *param = AliTPCReconstructor::GetRecoParam()->GetSystematicError();
6877 Double_t *covarIn= (Double_t*)seed->GetCovariance();
6879 for (Int_t i=0;i<15;i++) covar[i]=0;
6885 covar[0] = param[0]*param[0];
6886 covar[2] = param[1]*param[1];
6887 covar[5] = param[2]*param[2];
6888 covar[9] = param[3]*param[3];
6889 Double_t facC = AliTracker::GetBz()*kB2C;
6890 covar[14]= param[4]*param[4]*facC*facC;
6892 covar[1]=TMath::Sqrt((covar[0]*covar[2]))*covarIn[1]/TMath::Sqrt((covarIn[0]*covarIn[2]));
6894 covar[3]=TMath::Sqrt((covar[0]*covar[5]))*covarIn[3]/TMath::Sqrt((covarIn[0]*covarIn[5]));
6895 covar[4]=TMath::Sqrt((covar[2]*covar[5]))*covarIn[4]/TMath::Sqrt((covarIn[2]*covarIn[5]));
6897 covar[6]=TMath::Sqrt((covar[0]*covar[9]))*covarIn[6]/TMath::Sqrt((covarIn[0]*covarIn[9]));
6898 covar[7]=TMath::Sqrt((covar[2]*covar[9]))*covarIn[7]/TMath::Sqrt((covarIn[2]*covarIn[9]));
6899 covar[8]=TMath::Sqrt((covar[5]*covar[9]))*covarIn[8]/TMath::Sqrt((covarIn[5]*covarIn[9]));
6901 covar[10]=TMath::Sqrt((covar[0]*covar[14]))*covarIn[10]/TMath::Sqrt((covarIn[0]*covarIn[14]));
6902 covar[11]=TMath::Sqrt((covar[2]*covar[14]))*covarIn[11]/TMath::Sqrt((covarIn[2]*covarIn[14]));
6903 covar[12]=TMath::Sqrt((covar[5]*covar[14]))*covarIn[12]/TMath::Sqrt((covarIn[5]*covarIn[14]));
6904 covar[13]=TMath::Sqrt((covar[9]*covar[14]))*covarIn[13]/TMath::Sqrt((covarIn[9]*covarIn[14]));
6906 seed->AddCovariance(covar);