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) //sign not accepted clusters
225 track->SetCurrentClusterIndex1(track->GetCurrentClusterIndex1() | 0x8000);
226 else // unsign accpeted clusters
227 track->SetCurrentClusterIndex1(track->GetCurrentClusterIndex1() & 0xffff7fff);
228 UInt_t i = track->GetCurrentClusterIndex1();
230 Int_t sec=(i&0xff000000)>>24;
231 //Int_t row = (i&0x00ff0000)>>16;
232 track->SetRow((i&0x00ff0000)>>16);
233 track->SetSector(sec);
234 // Int_t index = i&0xFFFF;
235 if (sec>=fkParam->GetNInnerSector()) track->SetRow(track->GetRow()+fkParam->GetNRowLow());
236 track->SetClusterIndex2(track->GetRow(), i);
237 //track->fFirstPoint = row;
238 //if ( track->fLastPoint<row) track->fLastPoint =row;
239 // if (track->fRow<0 || track->fRow>160) {
240 // printf("problem\n");
242 if (track->GetFirstPoint()>track->GetRow())
243 track->SetFirstPoint(track->GetRow());
244 if (track->GetLastPoint()<track->GetRow())
245 track->SetLastPoint(track->GetRow());
248 track->SetClusterPointer(track->GetRow(),c);
251 Double_t angle2 = track->GetSnp()*track->GetSnp();
253 //SET NEW Track Point
255 if (angle2<1) //PH sometimes angle2 is very big. To be investigated...
257 angle2 = TMath::Sqrt(angle2/(1-angle2));
258 AliTPCTrackerPoint &point =*(track->GetTrackPoint(track->GetRow()));
260 point.SetSigmaY(c->GetSigmaY2()/track->GetCurrentSigmaY2());
261 point.SetSigmaZ(c->GetSigmaZ2()/track->GetCurrentSigmaZ2());
262 point.SetErrY(sqrt(track->GetErrorY2()));
263 point.SetErrZ(sqrt(track->GetErrorZ2()));
265 point.SetX(track->GetX());
266 point.SetY(track->GetY());
267 point.SetZ(track->GetZ());
268 point.SetAngleY(angle2);
269 point.SetAngleZ(track->GetTgl());
270 if (point.IsShared()){
271 track->SetErrorY2(track->GetErrorY2()*4);
272 track->SetErrorZ2(track->GetErrorZ2()*4);
276 Double_t chi2 = track->GetPredictedChi2(track->GetCurrentCluster());
278 // track->SetErrorY2(track->GetErrorY2()*1.3);
279 // track->SetErrorY2(track->GetErrorY2()+0.01);
280 // track->SetErrorZ2(track->GetErrorZ2()*1.3);
281 // track->SetErrorZ2(track->GetErrorZ2()+0.005);
283 if (accept>0) return 0;
284 if (track->GetNumberOfClusters()%20==0){
285 // if (track->fHelixIn){
286 // TClonesArray & larr = *(track->fHelixIn);
287 // Int_t ihelix = larr.GetEntriesFast();
288 // new(larr[ihelix]) AliHelix(*track) ;
291 track->SetNoCluster(0);
292 return track->Update(c,chi2,i);
297 Int_t AliTPCtrackerMI::AcceptCluster(AliTPCseed * seed, AliTPCclusterMI * cluster)
300 // decide according desired precision to accept given
301 // cluster for tracking
303 seed->GetProlongation(cluster->GetX(),yt,zt);
304 Double_t sy2=ErrY2(seed,cluster);
305 Double_t sz2=ErrZ2(seed,cluster);
307 Double_t sdistancey2 = sy2+seed->GetSigmaY2();
308 Double_t sdistancez2 = sz2+seed->GetSigmaZ2();
309 Double_t dy=seed->GetCurrentCluster()->GetY()-yt;
310 Double_t dz=seed->GetCurrentCluster()->GetZ()-zt;
311 Double_t rdistancey2 = (seed->GetCurrentCluster()->GetY()-yt)*
312 (seed->GetCurrentCluster()->GetY()-yt)/sdistancey2;
313 Double_t rdistancez2 = (seed->GetCurrentCluster()->GetZ()-zt)*
314 (seed->GetCurrentCluster()->GetZ()-zt)/sdistancez2;
316 Double_t rdistance2 = rdistancey2+rdistancez2;
319 if (AliTPCReconstructor::StreamLevel()>2 && seed->GetNumberOfClusters()>20) {
320 Float_t rmsy2 = seed->GetCurrentSigmaY2();
321 Float_t rmsz2 = seed->GetCurrentSigmaZ2();
322 Float_t rmsy2p30 = seed->GetCMeanSigmaY2p30();
323 Float_t rmsz2p30 = seed->GetCMeanSigmaZ2p30();
324 Float_t rmsy2p30R = seed->GetCMeanSigmaY2p30R();
325 Float_t rmsz2p30R = seed->GetCMeanSigmaZ2p30R();
326 AliExternalTrackParam param(*seed);
327 static TVectorD gcl(3),gtr(3);
329 param.GetXYZ(gcl.GetMatrixArray());
330 cluster->GetGlobalXYZ(gclf);
331 gcl[0]=gclf[0]; gcl[1]=gclf[1]; gcl[2]=gclf[2];
334 if (AliTPCReconstructor::StreamLevel()>2) {
335 (*fDebugStreamer)<<"ErrParam"<<
348 "rmsy2p30="<<rmsy2p30<<
349 "rmsz2p30="<<rmsz2p30<<
350 "rmsy2p30R="<<rmsy2p30R<<
351 "rmsz2p30R="<<rmsz2p30R<<
352 // normalize distance -
353 "rdisty="<<rdistancey2<<
354 "rdistz="<<rdistancez2<<
355 "rdist="<<rdistance2<< //
359 //return 0; // temporary
360 if (rdistance2>32) return 3;
363 if ((rdistancey2>9. || rdistancez2>9.) && cluster->GetType()==0)
364 return 2; //suspisiouce - will be changed
366 if ((rdistancey2>6.25 || rdistancez2>6.25) && cluster->GetType()>0)
367 // strict cut on overlaped cluster
368 return 2; //suspisiouce - will be changed
370 if ( (rdistancey2>1. || rdistancez2>6.25 )
371 && cluster->GetType()<0){
372 seed->SetNFoundable(seed->GetNFoundable()-1);
376 if (AliTPCReconstructor::GetRecoParam()->GetUseHLTClusters() == 3 ||
377 AliTPCReconstructor::GetRecoParam()->GetUseHLTClusters() == 4 ) {
378 if(!AliTPCReconstructor::GetRecoParam()->GetUseOnePadCluster())
379 if (TMath::Abs(cluster->GetSigmaY2()) < kAlmost0)
389 //_____________________________________________________________________________
390 AliTPCtrackerMI::AliTPCtrackerMI(const AliTPCParam *par):
392 fkNIS(par->GetNInnerSector()/2),
394 fkNOS(par->GetNOuterSector()/2),
411 //---------------------------------------------------------------------
412 // The main TPC tracker constructor
413 //---------------------------------------------------------------------
414 fInnerSec=new AliTPCtrackerSector[fkNIS];
415 fOuterSec=new AliTPCtrackerSector[fkNOS];
418 for (i=0; i<fkNIS; i++) fInnerSec[i].Setup(par,0);
419 for (i=0; i<fkNOS; i++) fOuterSec[i].Setup(par,1);
422 Int_t nrowlow = par->GetNRowLow();
423 Int_t nrowup = par->GetNRowUp();
426 for (i=0;i<nrowlow;i++){
427 fXRow[i] = par->GetPadRowRadiiLow(i);
428 fPadLength[i]= par->GetPadPitchLength(0,i);
429 fYMax[i] = fXRow[i]*TMath::Tan(0.5*par->GetInnerAngle());
433 for (i=0;i<nrowup;i++){
434 fXRow[i+nrowlow] = par->GetPadRowRadiiUp(i);
435 fPadLength[i+nrowlow] = par->GetPadPitchLength(60,i);
436 fYMax[i+nrowlow] = fXRow[i+nrowlow]*TMath::Tan(0.5*par->GetOuterAngle());
439 if (AliTPCReconstructor::StreamLevel()>0) {
440 fDebugStreamer = new TTreeSRedirector("TPCdebug.root");
443 //________________________________________________________________________
444 AliTPCtrackerMI::AliTPCtrackerMI(const AliTPCtrackerMI &t):
465 //------------------------------------
466 // dummy copy constructor
467 //------------------------------------------------------------------
469 for (Int_t irow=0; irow<200; irow++){
476 AliTPCtrackerMI & AliTPCtrackerMI::operator=(const AliTPCtrackerMI& /*r*/)
478 //------------------------------
480 //--------------------------------------------------------------
483 //_____________________________________________________________________________
484 AliTPCtrackerMI::~AliTPCtrackerMI() {
485 //------------------------------------------------------------------
486 // TPC tracker destructor
487 //------------------------------------------------------------------
494 if (fDebugStreamer) delete fDebugStreamer;
498 void AliTPCtrackerMI::FillESD(const TObjArray* arr)
502 //fill esds using updated tracks
505 // write tracks to the event
506 // store index of the track
507 Int_t nseed=arr->GetEntriesFast();
508 //FindKinks(arr,fEvent);
509 for (Int_t i=0; i<nseed; i++) {
510 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
514 if (AliTPCReconstructor::StreamLevel()>1) {
515 (*fDebugStreamer)<<"Track0"<<
519 // pt->PropagateTo(fkParam->GetInnerRadiusLow());
520 if (pt->GetKinkIndex(0)<=0){ //don't propagate daughter tracks
521 pt->PropagateTo(fkParam->GetInnerRadiusLow());
524 if (( pt->GetPoints()[2]- pt->GetPoints()[0])>5 && pt->GetPoints()[3]>0.8){
526 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
527 iotrack.SetTPCPoints(pt->GetPoints());
528 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
529 iotrack.SetV0Indexes(pt->GetV0Indexes());
530 // iotrack.SetTPCpid(pt->fTPCr);
531 //iotrack.SetTPCindex(i);
532 fEvent->AddTrack(&iotrack);
536 if ( (pt->GetNumberOfClusters()>70)&& (Float_t(pt->GetNumberOfClusters())/Float_t(pt->GetNFoundable()))>0.55) {
538 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
539 iotrack.SetTPCPoints(pt->GetPoints());
540 //iotrack.SetTPCindex(i);
541 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
542 iotrack.SetV0Indexes(pt->GetV0Indexes());
543 // iotrack.SetTPCpid(pt->fTPCr);
544 fEvent->AddTrack(&iotrack);
548 // short tracks - maybe decays
550 if ( (pt->GetNumberOfClusters()>30) && (Float_t(pt->GetNumberOfClusters())/Float_t(pt->GetNFoundable()))>0.70) {
551 Int_t found,foundable,shared;
552 pt->GetClusterStatistic(0,60,found, foundable,shared,kFALSE);
553 if ( (found>20) && (pt->GetNShared()/float(pt->GetNumberOfClusters())<0.2)){
555 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
556 //iotrack.SetTPCindex(i);
557 iotrack.SetTPCPoints(pt->GetPoints());
558 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
559 iotrack.SetV0Indexes(pt->GetV0Indexes());
560 //iotrack.SetTPCpid(pt->fTPCr);
561 fEvent->AddTrack(&iotrack);
566 if ( (pt->GetNumberOfClusters()>20) && (Float_t(pt->GetNumberOfClusters())/Float_t(pt->GetNFoundable()))>0.8) {
567 Int_t found,foundable,shared;
568 pt->GetClusterStatistic(0,60,found, foundable,shared,kFALSE);
569 if (found<20) continue;
570 if (pt->GetNShared()/float(pt->GetNumberOfClusters())>0.2) continue;
573 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
574 iotrack.SetTPCPoints(pt->GetPoints());
575 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
576 iotrack.SetV0Indexes(pt->GetV0Indexes());
577 //iotrack.SetTPCpid(pt->fTPCr);
578 //iotrack.SetTPCindex(i);
579 fEvent->AddTrack(&iotrack);
582 // short tracks - secondaties
584 if ( (pt->GetNumberOfClusters()>30) ) {
585 Int_t found,foundable,shared;
586 pt->GetClusterStatistic(128,158,found, foundable,shared,kFALSE);
587 if ( (found>20) && (pt->GetNShared()/float(pt->GetNumberOfClusters())<0.2) &&float(found)/float(foundable)>0.8){
589 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
590 iotrack.SetTPCPoints(pt->GetPoints());
591 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
592 iotrack.SetV0Indexes(pt->GetV0Indexes());
593 //iotrack.SetTPCpid(pt->fTPCr);
594 //iotrack.SetTPCindex(i);
595 fEvent->AddTrack(&iotrack);
600 if ( (pt->GetNumberOfClusters()>15)) {
601 Int_t found,foundable,shared;
602 pt->GetClusterStatistic(138,158,found, foundable,shared,kFALSE);
603 if (found<15) continue;
604 if (foundable<=0) continue;
605 if (pt->GetNShared()/float(pt->GetNumberOfClusters())>0.2) continue;
606 if (float(found)/float(foundable)<0.8) continue;
609 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
610 iotrack.SetTPCPoints(pt->GetPoints());
611 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
612 iotrack.SetV0Indexes(pt->GetV0Indexes());
613 // iotrack.SetTPCpid(pt->fTPCr);
614 //iotrack.SetTPCindex(i);
615 fEvent->AddTrack(&iotrack);
619 // >> account for suppressed tracks in the kink indices (RS)
620 int nESDtracks = fEvent->GetNumberOfTracks();
621 for (int it=nESDtracks;it--;) {
622 AliESDtrack* esdTr = fEvent->GetTrack(it);
623 if (!esdTr || !esdTr->GetKinkIndex(0)) continue;
624 for (int ik=0;ik<3;ik++) {
626 if (!(knkId=esdTr->GetKinkIndex(ik))) break; // no more kinks for this track
627 AliESDkink* kink = fEvent->GetKink(TMath::Abs(knkId)-1);
629 AliError(Form("ESDTrack%d refers to non-existing kink %d",it,TMath::Abs(knkId)-1));
632 kink->SetIndex(it, knkId<0 ? 0:1); // update track index of the kink: mother at 0, daughter at 1
635 // << account for suppressed tracks in the kink indices (RS)
636 AliInfo(Form("Number of filled ESDs-\t%d\n",fEvent->GetNumberOfTracks()));
644 Double_t AliTPCtrackerMI::ErrY2(AliTPCseed* seed, const AliTPCclusterMI * cl){
647 // Use calibrated cluster error from OCDB
649 AliTPCClusterParam * clparam = AliTPCcalibDB::Instance()->GetClusterParam();
651 Float_t z = TMath::Abs(fkParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
652 Int_t ctype = cl->GetType();
653 Int_t type = (cl->GetRow()<63) ? 0: (cl->GetRow()>126) ? 1:2;
654 Double_t angle = seed->GetSnp()*seed->GetSnp();
655 angle = TMath::Sqrt(TMath::Abs(angle/(1.-angle)));
656 Double_t erry2 = clparam->GetError0Par(0,type, z,angle);
658 erry2+=0.5; // edge cluster
661 seed->SetErrorY2(erry2);
665 //calculate look-up table at the beginning
666 // static Bool_t ginit = kFALSE;
667 // static Float_t gnoise1,gnoise2,gnoise3;
668 // static Float_t ggg1[10000];
669 // static Float_t ggg2[10000];
670 // static Float_t ggg3[10000];
671 // static Float_t glandau1[10000];
672 // static Float_t glandau2[10000];
673 // static Float_t glandau3[10000];
675 // static Float_t gcor01[500];
676 // static Float_t gcor02[500];
677 // static Float_t gcorp[500];
681 // if (ginit==kFALSE){
682 // for (Int_t i=1;i<500;i++){
683 // Float_t rsigma = float(i)/100.;
684 // gcor02[i] = TMath::Max(0.78 +TMath::Exp(7.4*(rsigma-1.2)),0.6);
685 // gcor01[i] = TMath::Max(0.72 +TMath::Exp(3.36*(rsigma-1.2)),0.6);
686 // gcorp[i] = TMath::Max(TMath::Power((rsigma+0.5),1.5),1.2);
690 // for (Int_t i=3;i<10000;i++){
694 // Float_t amp = float(i);
695 // Float_t padlength =0.75;
696 // gnoise1 = 0.0004/padlength;
697 // Float_t nel = 0.268*amp;
698 // Float_t nprim = 0.155*amp;
699 // ggg1[i] = fkParam->GetDiffT()*fkParam->GetDiffT()*(2+0.001*nel/(padlength*padlength))/nel;
700 // glandau1[i] = (2.+0.12*nprim)*0.5* (2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
701 // if (glandau1[i]>1) glandau1[i]=1;
702 // glandau1[i]*=padlength*padlength/12.;
706 // gnoise2 = 0.0004/padlength;
708 // nprim = 0.133*amp;
709 // ggg2[i] = fkParam->GetDiffT()*fkParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
710 // glandau2[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
711 // if (glandau2[i]>1) glandau2[i]=1;
712 // glandau2[i]*=padlength*padlength/12.;
717 // gnoise3 = 0.0004/padlength;
719 // nprim = 0.133*amp;
720 // ggg3[i] = fkParam->GetDiffT()*fkParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
721 // glandau3[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
722 // if (glandau3[i]>1) glandau3[i]=1;
723 // glandau3[i]*=padlength*padlength/12.;
731 // Int_t amp = int(TMath::Abs(cl->GetQ()));
733 // seed->SetErrorY2(1.);
737 // Float_t z = TMath::Abs(fkParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
738 // Int_t ctype = cl->GetType();
739 // Float_t padlength= GetPadPitchLength(seed->GetRow());
740 // Double_t angle2 = seed->GetSnp()*seed->GetSnp();
741 // angle2 = angle2/(1-angle2);
743 // //cluster "quality"
744 // Int_t rsigmay = int(100.*cl->GetSigmaY2()/(seed->GetCurrentSigmaY2()));
747 // if (fSectors==fInnerSec){
748 // snoise2 = gnoise1;
749 // res = ggg1[amp]*z+glandau1[amp]*angle2;
750 // if (ctype==0) res *= gcor01[rsigmay];
753 // res*= gcorp[rsigmay];
757 // if (padlength<1.1){
758 // snoise2 = gnoise2;
759 // res = ggg2[amp]*z+glandau2[amp]*angle2;
760 // if (ctype==0) res *= gcor02[rsigmay];
763 // res*= gcorp[rsigmay];
767 // snoise2 = gnoise3;
768 // res = ggg3[amp]*z+glandau3[amp]*angle2;
769 // if (ctype==0) res *= gcor02[rsigmay];
772 // res*= gcorp[rsigmay];
779 // res*=2.4; // overestimate error 2 times
783 // if (res<2*snoise2)
786 // seed->SetErrorY2(res);
794 Double_t AliTPCtrackerMI::ErrZ2(AliTPCseed* seed, const AliTPCclusterMI * cl){
797 // Use calibrated cluster error from OCDB
799 AliTPCClusterParam * clparam = AliTPCcalibDB::Instance()->GetClusterParam();
801 Float_t z = TMath::Abs(fkParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
802 Int_t ctype = cl->GetType();
803 Int_t type = (cl->GetRow()<63) ? 0: (cl->GetRow()>126) ? 1:2;
805 Double_t angle2 = seed->GetSnp()*seed->GetSnp();
806 angle2 = seed->GetTgl()*seed->GetTgl()*(1+angle2/(1-angle2));
807 Double_t angle = TMath::Sqrt(TMath::Abs(angle2));
808 Double_t errz2 = clparam->GetError0Par(1,type, z,angle);
810 errz2+=0.5; // edge cluster
813 seed->SetErrorZ2(errz2);
819 // //seed->SetErrorY2(0.1);
821 // //calculate look-up table at the beginning
822 // static Bool_t ginit = kFALSE;
823 // static Float_t gnoise1,gnoise2,gnoise3;
824 // static Float_t ggg1[10000];
825 // static Float_t ggg2[10000];
826 // static Float_t ggg3[10000];
827 // static Float_t glandau1[10000];
828 // static Float_t glandau2[10000];
829 // static Float_t glandau3[10000];
831 // static Float_t gcor01[1000];
832 // static Float_t gcor02[1000];
833 // static Float_t gcorp[1000];
837 // if (ginit==kFALSE){
838 // for (Int_t i=1;i<1000;i++){
839 // Float_t rsigma = float(i)/100.;
840 // gcor02[i] = TMath::Max(0.81 +TMath::Exp(6.8*(rsigma-1.2)),0.6);
841 // gcor01[i] = TMath::Max(0.72 +TMath::Exp(2.04*(rsigma-1.2)),0.6);
842 // gcorp[i] = TMath::Max(TMath::Power((rsigma+0.5),1.5),1.2);
846 // for (Int_t i=3;i<10000;i++){
850 // Float_t amp = float(i);
851 // Float_t padlength =0.75;
852 // gnoise1 = 0.0004/padlength;
853 // Float_t nel = 0.268*amp;
854 // Float_t nprim = 0.155*amp;
855 // ggg1[i] = fkParam->GetDiffT()*fkParam->GetDiffT()*(2+0.001*nel/(padlength*padlength))/nel;
856 // glandau1[i] = (2.+0.12*nprim)*0.5* (2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
857 // if (glandau1[i]>1) glandau1[i]=1;
858 // glandau1[i]*=padlength*padlength/12.;
862 // gnoise2 = 0.0004/padlength;
864 // nprim = 0.133*amp;
865 // ggg2[i] = fkParam->GetDiffT()*fkParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
866 // glandau2[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
867 // if (glandau2[i]>1) glandau2[i]=1;
868 // glandau2[i]*=padlength*padlength/12.;
873 // gnoise3 = 0.0004/padlength;
875 // nprim = 0.133*amp;
876 // ggg3[i] = fkParam->GetDiffT()*fkParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
877 // glandau3[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
878 // if (glandau3[i]>1) glandau3[i]=1;
879 // glandau3[i]*=padlength*padlength/12.;
887 // Int_t amp = int(TMath::Abs(cl->GetQ()));
889 // seed->SetErrorY2(1.);
893 // Float_t z = TMath::Abs(fkParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
894 // Int_t ctype = cl->GetType();
895 // Float_t padlength= GetPadPitchLength(seed->GetRow());
897 // Double_t angle2 = seed->GetSnp()*seed->GetSnp();
898 // // if (angle2<0.6) angle2 = 0.6;
899 // angle2 = seed->GetTgl()*seed->GetTgl()*(1+angle2/(1-angle2));
901 // //cluster "quality"
902 // Int_t rsigmaz = int(100.*cl->GetSigmaZ2()/(seed->GetCurrentSigmaZ2()));
905 // if (fSectors==fInnerSec){
906 // snoise2 = gnoise1;
907 // res = ggg1[amp]*z+glandau1[amp]*angle2;
908 // if (ctype==0) res *= gcor01[rsigmaz];
911 // res*= gcorp[rsigmaz];
915 // if (padlength<1.1){
916 // snoise2 = gnoise2;
917 // res = ggg2[amp]*z+glandau2[amp]*angle2;
918 // if (ctype==0) res *= gcor02[rsigmaz];
921 // res*= gcorp[rsigmaz];
925 // snoise2 = gnoise3;
926 // res = ggg3[amp]*z+glandau3[amp]*angle2;
927 // if (ctype==0) res *= gcor02[rsigmaz];
930 // res*= gcorp[rsigmaz];
939 // if ((ctype<0) &&<70){
944 // if (res<2*snoise2)
946 // if (res>3) res =3;
947 // seed->SetErrorZ2(res);
955 void AliTPCtrackerMI::RotateToLocal(AliTPCseed *seed)
957 //rotate to track "local coordinata
958 Float_t x = seed->GetX();
959 Float_t y = seed->GetY();
960 Float_t ymax = x*TMath::Tan(0.5*fSectors->GetAlpha());
963 seed->SetRelativeSector((seed->GetRelativeSector()+1) % fN);
964 if (!seed->Rotate(fSectors->GetAlpha()))
966 } else if (y <-ymax) {
967 seed->SetRelativeSector((seed->GetRelativeSector()-1+fN) % fN);
968 if (!seed->Rotate(-fSectors->GetAlpha()))
976 //_____________________________________________________________________________
977 Double_t AliTPCtrackerMI::F1old(Double_t x1,Double_t y1,
978 Double_t x2,Double_t y2,
979 Double_t x3,Double_t y3) const
981 //-----------------------------------------------------------------
982 // Initial approximation of the track curvature
983 //-----------------------------------------------------------------
984 Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
985 Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
986 (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
987 Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
988 (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
990 Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
991 if ( xr*xr+yr*yr<=0.00000000000001) return 100;
992 return -xr*yr/sqrt(xr*xr+yr*yr);
997 //_____________________________________________________________________________
998 Double_t AliTPCtrackerMI::F1(Double_t x1,Double_t y1,
999 Double_t x2,Double_t y2,
1000 Double_t x3,Double_t y3) const
1002 //-----------------------------------------------------------------
1003 // Initial approximation of the track curvature
1004 //-----------------------------------------------------------------
1010 Double_t det = x3*y2-x2*y3;
1011 if (TMath::Abs(det)<1e-10){
1015 Double_t u = 0.5* (x2*(x2-x3)+y2*(y2-y3))/det;
1016 Double_t x0 = x3*0.5-y3*u;
1017 Double_t y0 = y3*0.5+x3*u;
1018 Double_t c2 = 1/TMath::Sqrt(x0*x0+y0*y0);
1024 Double_t AliTPCtrackerMI::F2(Double_t x1,Double_t y1,
1025 Double_t x2,Double_t y2,
1026 Double_t x3,Double_t y3) const
1028 //-----------------------------------------------------------------
1029 // Initial approximation of the track curvature
1030 //-----------------------------------------------------------------
1036 Double_t det = x3*y2-x2*y3;
1037 if (TMath::Abs(det)<1e-10) {
1041 Double_t u = 0.5* (x2*(x2-x3)+y2*(y2-y3))/det;
1042 Double_t x0 = x3*0.5-y3*u;
1043 Double_t y0 = y3*0.5+x3*u;
1044 Double_t c2 = 1/TMath::Sqrt(x0*x0+y0*y0);
1053 //_____________________________________________________________________________
1054 Double_t AliTPCtrackerMI::F2old(Double_t x1,Double_t y1,
1055 Double_t x2,Double_t y2,
1056 Double_t x3,Double_t y3) const
1058 //-----------------------------------------------------------------
1059 // Initial approximation of the track curvature times center of curvature
1060 //-----------------------------------------------------------------
1061 Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
1062 Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
1063 (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
1064 Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
1065 (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
1067 Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
1069 return -a/(d*y1-b)*xr/sqrt(xr*xr+yr*yr);
1072 //_____________________________________________________________________________
1073 Double_t AliTPCtrackerMI::F3(Double_t x1,Double_t y1,
1074 Double_t x2,Double_t y2,
1075 Double_t z1,Double_t z2) const
1077 //-----------------------------------------------------------------
1078 // Initial approximation of the tangent of the track dip angle
1079 //-----------------------------------------------------------------
1080 return (z1 - z2)/sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
1084 Double_t AliTPCtrackerMI::F3n(Double_t x1,Double_t y1,
1085 Double_t x2,Double_t y2,
1086 Double_t z1,Double_t z2, Double_t c) const
1088 //-----------------------------------------------------------------
1089 // Initial approximation of the tangent of the track dip angle
1090 //-----------------------------------------------------------------
1094 //angle1 = (z1-z2)*c/(TMath::ASin(c*x1-ni)-TMath::ASin(c*x2-ni));
1096 Double_t d = TMath::Sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
1097 if (TMath::Abs(d*c*0.5)>1) return 0;
1098 // Double_t angle2 = TMath::ASin(d*c*0.5);
1099 // Double_t angle2 = AliTPCFastMath::FastAsin(d*c*0.5);
1100 Double_t angle2 = (d*c*0.5>0.1)? TMath::ASin(d*c*0.5): AliTPCFastMath::FastAsin(d*c*0.5);
1102 angle2 = (z1-z2)*c/(angle2*2.);
1106 Bool_t AliTPCtrackerMI::GetProlongation(Double_t x1, Double_t x2, Double_t x[5], Double_t &y, Double_t &z) const
1107 {//-----------------------------------------------------------------
1108 // This function find proloncation of a track to a reference plane x=x2.
1109 //-----------------------------------------------------------------
1113 if (TMath::Abs(x[4]*x1 - x[2]) >= 0.999) {
1117 Double_t c1=x[4]*x1 - x[2], r1=TMath::Sqrt((1.-c1)*(1.+c1));
1118 Double_t c2=x[4]*x2 - x[2], r2=TMath::Sqrt((1.-c2)*(1.+c2));
1122 Double_t dy = dx*(c1+c2)/(r1+r2);
1125 Double_t delta = x[4]*dx*(c1+c2)/(c1*r2 + c2*r1);
1127 if (TMath::Abs(delta)>0.01){
1128 dz = x[3]*TMath::ASin(delta)/x[4];
1130 dz = x[3]*AliTPCFastMath::FastAsin(delta)/x[4];
1133 //dz = x[3]*AliTPCFastMath::FastAsin(delta)/x[4];
1141 Int_t AliTPCtrackerMI::LoadClusters (TTree *const tree)
1146 return LoadClusters();
1150 Int_t AliTPCtrackerMI::LoadClusters(const TObjArray *arr)
1153 // load clusters to the memory
1154 AliTPCClustersRow *clrow = new AliTPCClustersRow("AliTPCclusterMI");
1155 Int_t lower = arr->LowerBound();
1156 Int_t entries = arr->GetEntriesFast();
1158 for (Int_t i=lower; i<entries; i++) {
1159 clrow = (AliTPCClustersRow*) arr->At(i);
1160 if(!clrow) continue;
1161 if(!clrow->GetArray()) continue;
1165 fkParam->AdjustSectorRow(clrow->GetID(),sec,row);
1167 for (Int_t icl=0; icl<clrow->GetArray()->GetEntriesFast(); icl++){
1168 Transform((AliTPCclusterMI*)(clrow->GetArray()->At(icl)));
1171 if (clrow->GetArray()->GetEntriesFast()<=0) continue;
1172 AliTPCtrackerRow * tpcrow=0;
1175 tpcrow = &(fInnerSec[sec%fkNIS][row]);
1179 tpcrow = &(fOuterSec[(sec-fkNIS*2)%fkNOS][row]);
1180 left = (sec-fkNIS*2)/fkNOS;
1183 tpcrow->SetN1(clrow->GetArray()->GetEntriesFast());
1184 tpcrow->SetClusters1(new AliTPCclusterMI[tpcrow->GetN1()]);
1185 for (Int_t j=0;j<tpcrow->GetN1();++j)
1186 tpcrow->SetCluster1(j, *(AliTPCclusterMI*)(clrow->GetArray()->At(j)));
1189 tpcrow->SetN2(clrow->GetArray()->GetEntriesFast());
1190 tpcrow->SetClusters2(new AliTPCclusterMI[tpcrow->GetN2()]);
1191 for (Int_t j=0;j<tpcrow->GetN2();++j)
1192 tpcrow->SetCluster2(j,*(AliTPCclusterMI*)(clrow->GetArray()->At(j)));
1194 clrow->GetArray()->Clear("C");
1203 Int_t AliTPCtrackerMI::LoadClusters(const TClonesArray *arr)
1206 // load clusters to the memory from one
1209 AliTPCclusterMI *clust=0;
1210 Int_t count[72][96] = { {0} , {0} };
1212 // loop over clusters
1213 for (Int_t icl=0; icl<arr->GetEntriesFast(); icl++) {
1214 clust = (AliTPCclusterMI*)arr->At(icl);
1215 if(!clust) continue;
1216 //printf("cluster: det %d, row %d \n", clust->GetDetector(),clust->GetRow());
1218 // transform clusters
1221 // count clusters per pad row
1222 count[clust->GetDetector()][clust->GetRow()]++;
1225 // insert clusters to sectors
1226 for (Int_t icl=0; icl<arr->GetEntriesFast(); icl++) {
1227 clust = (AliTPCclusterMI*)arr->At(icl);
1228 if(!clust) continue;
1230 Int_t sec = clust->GetDetector();
1231 Int_t row = clust->GetRow();
1233 // filter overlapping pad rows needed by HLT
1234 if(sec<fkNIS*2) { //IROCs
1235 if(row == 30) continue;
1238 if(row == 27 || row == 76) continue;
1244 fInnerSec[sec%fkNIS].InsertCluster(clust, count[sec][row], fkParam);
1247 left = (sec-fkNIS*2)/fkNOS;
1248 fOuterSec[(sec-fkNIS*2)%fkNOS].InsertCluster(clust, count[sec][row], fkParam);
1252 // Load functions must be called behind LoadCluster(TClonesArray*)
1254 //LoadOuterSectors();
1255 //LoadInnerSectors();
1261 Int_t AliTPCtrackerMI::LoadClusters()
1264 // load clusters to the memory
1265 AliTPCClustersRow *clrow= new AliTPCClustersRow("AliTPCclusterMI");
1267 // TTree * tree = fClustersArray.GetTree();
1269 TTree * tree = fInput;
1270 TBranch * br = tree->GetBranch("Segment");
1271 br->SetAddress(&clrow);
1273 Int_t j=Int_t(tree->GetEntries());
1274 for (Int_t i=0; i<j; i++) {
1278 fkParam->AdjustSectorRow(clrow->GetID(),sec,row);
1279 for (Int_t icl=0; icl<clrow->GetArray()->GetEntriesFast(); icl++){
1280 Transform((AliTPCclusterMI*)(clrow->GetArray()->At(icl)));
1283 AliTPCtrackerRow * tpcrow=0;
1286 tpcrow = &(fInnerSec[sec%fkNIS][row]);
1290 tpcrow = &(fOuterSec[(sec-fkNIS*2)%fkNOS][row]);
1291 left = (sec-fkNIS*2)/fkNOS;
1294 tpcrow->SetN1(clrow->GetArray()->GetEntriesFast());
1295 tpcrow->SetClusters1(new AliTPCclusterMI[tpcrow->GetN1()]);
1296 for (Int_t k=0;k<tpcrow->GetN1();++k)
1297 tpcrow->SetCluster1(k, *(AliTPCclusterMI*)(clrow->GetArray()->At(k)));
1300 tpcrow->SetN2(clrow->GetArray()->GetEntriesFast());
1301 tpcrow->SetClusters2(new AliTPCclusterMI[tpcrow->GetN2()]);
1302 for (Int_t k=0;k<tpcrow->GetN2();++k)
1303 tpcrow->SetCluster2(k,*(AliTPCclusterMI*)(clrow->GetArray()->At(k)));
1314 void AliTPCtrackerMI::UnloadClusters()
1317 // unload clusters from the memory
1319 Int_t nrows = fOuterSec->GetNRows();
1320 for (Int_t sec = 0;sec<fkNOS;sec++)
1321 for (Int_t row = 0;row<nrows;row++){
1322 AliTPCtrackerRow* tpcrow = &(fOuterSec[sec%fkNOS][row]);
1324 // if (tpcrow->fClusters1) delete []tpcrow->fClusters1;
1325 // if (tpcrow->fClusters2) delete []tpcrow->fClusters2;
1327 tpcrow->ResetClusters();
1330 nrows = fInnerSec->GetNRows();
1331 for (Int_t sec = 0;sec<fkNIS;sec++)
1332 for (Int_t row = 0;row<nrows;row++){
1333 AliTPCtrackerRow* tpcrow = &(fInnerSec[sec%fkNIS][row]);
1335 // if (tpcrow->fClusters1) delete []tpcrow->fClusters1;
1336 //if (tpcrow->fClusters2) delete []tpcrow->fClusters2;
1338 tpcrow->ResetClusters();
1344 void AliTPCtrackerMI::FillClusterArray(TObjArray* array) const{
1346 // Filling cluster to the array - For visualization purposes
1349 nrows = fOuterSec->GetNRows();
1350 for (Int_t sec = 0;sec<fkNOS;sec++)
1351 for (Int_t row = 0;row<nrows;row++){
1352 AliTPCtrackerRow* tpcrow = &(fOuterSec[sec%fkNOS][row]);
1353 if (!tpcrow) continue;
1354 for (Int_t icl = 0;icl<tpcrow->GetN();icl++){
1355 array->AddLast((TObject*)((*tpcrow)[icl]));
1358 nrows = fInnerSec->GetNRows();
1359 for (Int_t sec = 0;sec<fkNIS;sec++)
1360 for (Int_t row = 0;row<nrows;row++){
1361 AliTPCtrackerRow* tpcrow = &(fInnerSec[sec%fkNIS][row]);
1362 if (!tpcrow) continue;
1363 for (Int_t icl = 0;icl<tpcrow->GetN();icl++){
1364 array->AddLast((TObject*)(*tpcrow)[icl]);
1370 void AliTPCtrackerMI::Transform(AliTPCclusterMI * cluster){
1374 AliTPCcalibDB * calibDB = AliTPCcalibDB::Instance();
1375 AliTPCTransform *transform = calibDB->GetTransform() ;
1377 AliFatal("Tranformations not in calibDB");
1380 transform->SetCurrentRecoParam((AliTPCRecoParam*)AliTPCReconstructor::GetRecoParam());
1381 Double_t x[3]={cluster->GetRow(),cluster->GetPad(),cluster->GetTimeBin()};
1382 Int_t i[1]={cluster->GetDetector()};
1383 transform->Transform(x,i,0,1);
1384 // if (cluster->GetDetector()%36>17){
1389 // in debug mode check the transformation
1391 if (AliTPCReconstructor::StreamLevel()>2) {
1393 cluster->GetGlobalXYZ(gx);
1394 Int_t event = (fEvent==NULL)? 0: fEvent->GetEventNumberInFile();
1395 TTreeSRedirector &cstream = *fDebugStreamer;
1396 cstream<<"Transform"<<
1407 cluster->SetX(x[0]);
1408 cluster->SetY(x[1]);
1409 cluster->SetZ(x[2]);
1414 //if (!fkParam->IsGeoRead()) fkParam->ReadGeoMatrices();
1415 if (AliTPCReconstructor::GetRecoParam()->GetUseSectorAlignment() && (!calibDB->HasAlignmentOCDB())){
1416 TGeoHMatrix *mat = fkParam->GetClusterMatrix(cluster->GetDetector());
1418 Double_t pos[3]= {cluster->GetX(),cluster->GetY(),cluster->GetZ()};
1419 Double_t posC[3]={cluster->GetX(),cluster->GetY(),cluster->GetZ()};
1420 if (mat) mat->LocalToMaster(pos,posC);
1422 // chack Loading of Geo matrices from GeoManager - TEMPORARY FIX
1424 cluster->SetX(posC[0]);
1425 cluster->SetY(posC[1]);
1426 cluster->SetZ(posC[2]);
1430 //_____________________________________________________________________________
1431 Int_t AliTPCtrackerMI::LoadOuterSectors() {
1432 //-----------------------------------------------------------------
1433 // This function fills outer TPC sectors with clusters.
1434 //-----------------------------------------------------------------
1435 Int_t nrows = fOuterSec->GetNRows();
1437 for (Int_t sec = 0;sec<fkNOS;sec++)
1438 for (Int_t row = 0;row<nrows;row++){
1439 AliTPCtrackerRow* tpcrow = &(fOuterSec[sec%fkNOS][row]);
1440 Int_t sec2 = sec+2*fkNIS;
1442 Int_t ncl = tpcrow->GetN1();
1444 AliTPCclusterMI *c= (tpcrow->GetCluster1(ncl));
1445 index=(((sec2<<8)+row)<<16)+ncl;
1446 tpcrow->InsertCluster(c,index);
1449 ncl = tpcrow->GetN2();
1451 AliTPCclusterMI *c= (tpcrow->GetCluster2(ncl));
1452 index=((((sec2+fkNOS)<<8)+row)<<16)+ncl;
1453 tpcrow->InsertCluster(c,index);
1456 // write indexes for fast acces
1458 for (Int_t i=0;i<510;i++)
1459 tpcrow->SetFastCluster(i,-1);
1460 for (Int_t i=0;i<tpcrow->GetN();i++){
1461 Int_t zi = Int_t((*tpcrow)[i]->GetZ()+255.);
1462 tpcrow->SetFastCluster(zi,i); // write index
1465 for (Int_t i=0;i<510;i++){
1466 if (tpcrow->GetFastCluster(i)<0)
1467 tpcrow->SetFastCluster(i,last);
1469 last = tpcrow->GetFastCluster(i);
1478 //_____________________________________________________________________________
1479 Int_t AliTPCtrackerMI::LoadInnerSectors() {
1480 //-----------------------------------------------------------------
1481 // This function fills inner TPC sectors with clusters.
1482 //-----------------------------------------------------------------
1483 Int_t nrows = fInnerSec->GetNRows();
1485 for (Int_t sec = 0;sec<fkNIS;sec++)
1486 for (Int_t row = 0;row<nrows;row++){
1487 AliTPCtrackerRow* tpcrow = &(fInnerSec[sec%fkNIS][row]);
1490 Int_t ncl = tpcrow->GetN1();
1492 AliTPCclusterMI *c= (tpcrow->GetCluster1(ncl));
1493 index=(((sec<<8)+row)<<16)+ncl;
1494 tpcrow->InsertCluster(c,index);
1497 ncl = tpcrow->GetN2();
1499 AliTPCclusterMI *c= (tpcrow->GetCluster2(ncl));
1500 index=((((sec+fkNIS)<<8)+row)<<16)+ncl;
1501 tpcrow->InsertCluster(c,index);
1504 // write indexes for fast acces
1506 for (Int_t i=0;i<510;i++)
1507 tpcrow->SetFastCluster(i,-1);
1508 for (Int_t i=0;i<tpcrow->GetN();i++){
1509 Int_t zi = Int_t((*tpcrow)[i]->GetZ()+255.);
1510 tpcrow->SetFastCluster(zi,i); // write index
1513 for (Int_t i=0;i<510;i++){
1514 if (tpcrow->GetFastCluster(i)<0)
1515 tpcrow->SetFastCluster(i,last);
1517 last = tpcrow->GetFastCluster(i);
1529 //_________________________________________________________________________
1530 AliTPCclusterMI *AliTPCtrackerMI::GetClusterMI(Int_t index) const {
1531 //--------------------------------------------------------------------
1532 // Return pointer to a given cluster
1533 //--------------------------------------------------------------------
1534 if (index<0) return 0; // no cluster
1535 Int_t sec=(index&0xff000000)>>24;
1536 Int_t row=(index&0x00ff0000)>>16;
1537 Int_t ncl=(index&0x00007fff)>>00;
1539 const AliTPCtrackerRow * tpcrow=0;
1540 AliTPCclusterMI * clrow =0;
1542 if (sec<0 || sec>=fkNIS*4) {
1543 AliWarning(Form("Wrong sector %d",sec));
1548 AliTPCtrackerSector& tracksec = fInnerSec[sec%fkNIS];
1549 if (tracksec.GetNRows()<=row) return 0;
1550 tpcrow = &(tracksec[row]);
1551 if (tpcrow==0) return 0;
1554 if (tpcrow->GetN1()<=ncl) return 0;
1555 clrow = tpcrow->GetClusters1();
1558 if (tpcrow->GetN2()<=ncl) return 0;
1559 clrow = tpcrow->GetClusters2();
1563 AliTPCtrackerSector& tracksec = fOuterSec[(sec-fkNIS*2)%fkNOS];
1564 if (tracksec.GetNRows()<=row) return 0;
1565 tpcrow = &(tracksec[row]);
1566 if (tpcrow==0) return 0;
1568 if (sec-2*fkNIS<fkNOS) {
1569 if (tpcrow->GetN1()<=ncl) return 0;
1570 clrow = tpcrow->GetClusters1();
1573 if (tpcrow->GetN2()<=ncl) return 0;
1574 clrow = tpcrow->GetClusters2();
1578 return &(clrow[ncl]);
1584 Int_t AliTPCtrackerMI::FollowToNext(AliTPCseed& t, Int_t nr) {
1585 //-----------------------------------------------------------------
1586 // This function tries to find a track prolongation to next pad row
1587 //-----------------------------------------------------------------
1589 Double_t x= GetXrow(nr), ymax=GetMaxY(nr);
1592 AliTPCclusterMI *cl=0;
1593 Int_t tpcindex= t.GetClusterIndex2(nr);
1595 // update current shape info every 5 pad-row
1596 // if ( (nr%5==0) || t.GetNumberOfClusters()<2 || (t.fCurrentSigmaY2<0.0001) ){
1600 if (fIteration>0 && tpcindex>=-1){ //if we have already clusters
1602 if (tpcindex==-1) return 0; //track in dead zone
1603 if (tpcindex >= 0){ //
1604 cl = t.GetClusterPointer(nr);
1605 if ( (cl==0) ) cl = GetClusterMI(tpcindex);
1606 t.SetCurrentClusterIndex1(tpcindex);
1609 Int_t relativesector = ((tpcindex&0xff000000)>>24)%18; // if previously accepted cluster in different sector
1610 Float_t angle = relativesector*fSectors->GetAlpha()+fSectors->GetAlphaShift();
1612 if (angle<-TMath::Pi()) angle += 2*TMath::Pi();
1613 if (angle>=TMath::Pi()) angle -= 2*TMath::Pi();
1615 if (TMath::Abs(angle-t.GetAlpha())>0.001){
1616 Double_t rotation = angle-t.GetAlpha();
1617 t.SetRelativeSector(relativesector);
1618 if (!t.Rotate(rotation)) return 0;
1620 if (!t.PropagateTo(x)) return 0;
1622 t.SetCurrentCluster(cl);
1624 Int_t accept = AcceptCluster(&t,t.GetCurrentCluster());
1625 if ((tpcindex&0x8000)==0) accept =0;
1627 //if founded cluster is acceptible
1628 if (cl->IsUsed(11)) { // id cluster is shared inrease uncertainty
1629 t.SetErrorY2(t.GetErrorY2()+0.03);
1630 t.SetErrorZ2(t.GetErrorZ2()+0.03);
1631 t.SetErrorY2(t.GetErrorY2()*3);
1632 t.SetErrorZ2(t.GetErrorZ2()*3);
1634 t.SetNFoundable(t.GetNFoundable()+1);
1635 UpdateTrack(&t,accept);
1638 else { // Remove old cluster from track
1639 t.SetClusterIndex(nr, -3);
1640 t.SetClusterPointer(nr, 0);
1644 if (TMath::Abs(t.GetSnp())>AliTPCReconstructor::GetMaxSnpTracker()) return 0; // cut on angle
1645 if (fIteration>1 && IsFindable(t)){
1646 // not look for new cluster during refitting
1647 t.SetNFoundable(t.GetNFoundable()+1);
1652 // if (TMath::Abs(t.GetSnp())>0.95 || TMath::Abs(x*t.GetC()-t.GetEta())>0.95) return 0;// patch 28 fev 06
1653 if (!t.PropagateTo(x)) {
1654 if (fIteration==0) t.SetRemoval(10);
1657 Double_t y = t.GetY();
1658 if (TMath::Abs(y)>ymax){
1660 t.SetRelativeSector((t.GetRelativeSector()+1) % fN);
1661 if (!t.Rotate(fSectors->GetAlpha()))
1663 } else if (y <-ymax) {
1664 t.SetRelativeSector((t.GetRelativeSector()-1+fN) % fN);
1665 if (!t.Rotate(-fSectors->GetAlpha()))
1668 if (!t.PropagateTo(x)) {
1669 if (fIteration==0) t.SetRemoval(10);
1675 Double_t z=t.GetZ();
1678 if (!IsActive(t.GetRelativeSector(),nr)) {
1680 t.SetClusterIndex2(nr,-1);
1683 //AliInfo(Form("A - Sector%d phi %f - alpha %f", t.fRelativeSector,y/x, t.GetAlpha()));
1684 Bool_t isActive = IsActive(t.GetRelativeSector(),nr);
1685 Bool_t isActive2 = (nr>=fInnerSec->GetNRows()) ? fOuterSec[t.GetRelativeSector()][nr-fInnerSec->GetNRows()].GetN()>0:fInnerSec[t.GetRelativeSector()][nr].GetN()>0;
1687 if (!isActive || !isActive2) return 0;
1689 const AliTPCtrackerRow &krow=GetRow(t.GetRelativeSector(),nr);
1690 if ( (t.GetSigmaY2()<0) || t.GetSigmaZ2()<0) return 0;
1692 Double_t roadz = 1.;
1694 if (TMath::Abs(TMath::Abs(y)-ymax)<krow.GetDeadZone()){
1696 t.SetClusterIndex2(nr,-1);
1702 // if (TMath::Abs(z)<(AliTPCReconstructor::GetCtgRange()*x+10) && TMath::Abs(z)<fkParam->GetZLength(0) && (TMath::Abs(t.GetSnp())<AliTPCReconstructor::GetMaxSnpTracker()))
1703 t.SetNFoundable(t.GetNFoundable()+1);
1709 // cl = krow.FindNearest2(y+10.,z,roady,roadz,index);
1710 cl = krow.FindNearest2(y,z,roady,roadz,index);
1711 if (cl) t.SetCurrentClusterIndex1(krow.GetIndex(index));
1714 t.SetCurrentCluster(cl);
1716 if (fIteration==2&&cl->IsUsed(10)) return 0;
1717 Int_t accept = AcceptCluster(&t,t.GetCurrentCluster());
1718 if (fIteration==2&&cl->IsUsed(11)) {
1719 t.SetErrorY2(t.GetErrorY2()+0.03);
1720 t.SetErrorZ2(t.GetErrorZ2()+0.03);
1721 t.SetErrorY2(t.GetErrorY2()*3);
1722 t.SetErrorZ2(t.GetErrorZ2()*3);
1725 if (t.fCurrentCluster->IsUsed(10)){
1730 if (t.fNShared>0.7*t.GetNumberOfClusters()) {
1736 if (accept<3) UpdateTrack(&t,accept);
1739 if ( fIteration==0 && t.GetNFoundable()*0.5 > t.GetNumberOfClusters()) t.SetRemoval(10);
1747 //_________________________________________________________________________
1748 Bool_t AliTPCtrackerMI::GetTrackPoint(Int_t index, AliTrackPoint &p ) const
1750 // Get track space point by index
1751 // return false in case the cluster doesn't exist
1752 AliTPCclusterMI *cl = GetClusterMI(index);
1753 if (!cl) return kFALSE;
1754 Int_t sector = (index&0xff000000)>>24;
1755 // Int_t row = (index&0x00ff0000)>>16;
1757 // xyz[0] = fkParam->GetPadRowRadii(sector,row);
1758 xyz[0] = cl->GetX();
1759 xyz[1] = cl->GetY();
1760 xyz[2] = cl->GetZ();
1762 fkParam->AdjustCosSin(sector,cos,sin);
1763 Float_t x = cos*xyz[0]-sin*xyz[1];
1764 Float_t y = cos*xyz[1]+sin*xyz[0];
1766 Float_t sigmaY2 = 0.027*cl->GetSigmaY2();
1767 if (sector < fkParam->GetNInnerSector()) sigmaY2 *= 2.07;
1768 Float_t sigmaZ2 = 0.066*cl->GetSigmaZ2();
1769 if (sector < fkParam->GetNInnerSector()) sigmaZ2 *= 1.77;
1770 cov[0] = sin*sin*sigmaY2;
1771 cov[1] = -sin*cos*sigmaY2;
1773 cov[3] = cos*cos*sigmaY2;
1776 p.SetXYZ(x,y,xyz[2],cov);
1777 AliGeomManager::ELayerID iLayer;
1779 if (sector < fkParam->GetNInnerSector()) {
1780 iLayer = AliGeomManager::kTPC1;
1784 iLayer = AliGeomManager::kTPC2;
1785 idet = sector - fkParam->GetNInnerSector();
1787 UShort_t volid = AliGeomManager::LayerToVolUID(iLayer,idet);
1788 p.SetVolumeID(volid);
1794 Int_t AliTPCtrackerMI::UpdateClusters(AliTPCseed& t, Int_t nr) {
1795 //-----------------------------------------------------------------
1796 // This function tries to find a track prolongation to next pad row
1797 //-----------------------------------------------------------------
1798 t.SetCurrentCluster(0);
1799 t.SetCurrentClusterIndex1(-3);
1801 Double_t xt=t.GetX();
1802 Int_t row = GetRowNumber(xt)-1;
1803 Double_t ymax= GetMaxY(nr);
1805 if (row < nr) return 1; // don't prolongate if not information until now -
1806 // if (TMath::Abs(t.GetSnp())>0.9 && t.GetNumberOfClusters()>40. && fIteration!=2) {
1808 // return 0; // not prolongate strongly inclined tracks
1810 // if (TMath::Abs(t.GetSnp())>0.95) {
1812 // return 0; // not prolongate strongly inclined tracks
1813 // }// patch 28 fev 06
1815 Double_t x= GetXrow(nr);
1817 //t.PropagateTo(x+0.02);
1818 //t.PropagateTo(x+0.01);
1819 if (!t.PropagateTo(x)){
1826 if (TMath::Abs(y)>ymax){
1828 t.SetRelativeSector((t.GetRelativeSector()+1) % fN);
1829 if (!t.Rotate(fSectors->GetAlpha()))
1831 } else if (y <-ymax) {
1832 t.SetRelativeSector((t.GetRelativeSector()-1+fN) % fN);
1833 if (!t.Rotate(-fSectors->GetAlpha()))
1836 // if (!t.PropagateTo(x)){
1843 if (TMath::Abs(t.GetSnp())>AliTPCReconstructor::GetMaxSnpTracker()) return 0;
1845 if (!IsActive(t.GetRelativeSector(),nr)) {
1847 t.SetClusterIndex2(nr,-1);
1850 //AliInfo(Form("A - Sector%d phi %f - alpha %f", t.fRelativeSector,y/x, t.GetAlpha()));
1852 AliTPCtrackerRow &krow=GetRow(t.GetRelativeSector(),nr);
1854 if (TMath::Abs(TMath::Abs(y)-ymax)<krow.GetDeadZone()){
1856 t.SetClusterIndex2(nr,-1);
1862 // if (TMath::Abs(t.GetZ())<(AliTPCReconstructor::GetCtgRange()*t.GetX()+10) && (TMath::Abs(t.GetSnp())<AliTPCReconstructor::GetMaxSnpTracker()))
1863 if (IsFindable(t)) t.SetNFoundable(t.GetNFoundable()+1);
1869 if ( (nr%2==0) || t.GetNumberOfClusters()<2){
1870 // t.fCurrentSigmaY = GetSigmaY(&t);
1871 //t.fCurrentSigmaZ = GetSigmaZ(&t);
1875 AliTPCclusterMI *cl=0;
1878 Double_t roady = 1.;
1879 Double_t roadz = 1.;
1883 index = t.GetClusterIndex2(nr);
1884 if ( (index>0) && (index&0x8000)==0){
1885 cl = t.GetClusterPointer(nr);
1886 if ( (cl==0) && (index>0)) cl = GetClusterMI(index);
1887 t.SetCurrentClusterIndex1(index);
1889 t.SetCurrentCluster(cl);
1895 // if (index<0) return 0;
1896 UInt_t uindex = TMath::Abs(index);
1899 //cl = krow.FindNearest2(y+10,z,roady,roadz,uindex);
1900 cl = krow.FindNearest2(y,z,roady,roadz,uindex);
1903 if (cl) t.SetCurrentClusterIndex1(krow.GetIndex(uindex));
1904 t.SetCurrentCluster(cl);
1910 Int_t AliTPCtrackerMI::FollowToNextCluster(AliTPCseed & t, Int_t nr) {
1911 //-----------------------------------------------------------------
1912 // This function tries to find a track prolongation to next pad row
1913 //-----------------------------------------------------------------
1915 //update error according neighborhoud
1917 if (t.GetCurrentCluster()) {
1919 Int_t accept = AcceptCluster(&t,t.GetCurrentCluster());
1921 if (t.GetCurrentCluster()->IsUsed(10)){
1926 t.SetNShared(t.GetNShared()+1);
1927 if (t.GetNShared()>0.7*t.GetNumberOfClusters()) {
1932 if (fIteration>0) accept = 0;
1933 if (accept<3) UpdateTrack(&t,accept);
1937 if ( t.GetNumberOfClusters()>18 && ( (t.GetSigmaY2()+t.GetSigmaZ2())>0.16)) t.SetRemoval(10);
1938 if ( t.GetNumberOfClusters()>18 && t.GetChi2()/t.GetNumberOfClusters()>6 ) t.SetRemoval(10);
1940 if (( (t.GetNFoundable()*0.5 > t.GetNumberOfClusters()) || t.GetNoCluster()>15)) t.SetRemoval(10);
1948 //_____________________________________________________________________________
1949 Int_t AliTPCtrackerMI::FollowProlongation(AliTPCseed& t, Int_t rf, Int_t step, Bool_t fromSeeds) {
1950 //-----------------------------------------------------------------
1951 // This function tries to find a track prolongation.
1952 //-----------------------------------------------------------------
1953 Double_t xt=t.GetX();
1955 Double_t alpha=t.GetAlpha();
1956 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
1957 if (alpha < 0. ) alpha += 2.*TMath::Pi();
1959 t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN);
1961 Int_t first = GetRowNumber(xt);
1966 for (Int_t nr= first; nr>=rf; nr-=step) {
1968 if (t.GetKinkIndexes()[0]>0){
1969 for (Int_t i=0;i<3;i++){
1970 Int_t index = t.GetKinkIndexes()[i];
1971 if (index==0) break;
1972 if (index<0) continue;
1974 AliKink * kink = (AliKink*)fEvent->GetKink(index-1);
1976 printf("PROBLEM\n");
1979 Int_t kinkrow = kink->GetTPCRow0()+2+Int_t(0.5/(0.05+kink->GetAngle(2)));
1981 AliExternalTrackParam paramd(t);
1982 kink->SetDaughter(paramd);
1983 kink->SetStatus(2,5);
1990 if (nr==80) t.UpdateReference();
1991 if (nr<fInnerSec->GetNRows())
1992 fSectors = fInnerSec;
1994 fSectors = fOuterSec;
1995 if (FollowToNext(t,nr)==0)
2008 Int_t AliTPCtrackerMI::FollowBackProlongation(AliTPCseed& t, Int_t rf, Bool_t fromSeeds) {
2009 //-----------------------------------------------------------------
2010 // This function tries to find a track prolongation.
2011 //-----------------------------------------------------------------
2013 Double_t xt=t.GetX();
2014 Double_t alpha=t.GetAlpha();
2015 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
2016 if (alpha < 0. ) alpha += 2.*TMath::Pi();
2017 t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha())%fN);
2019 Int_t first = t.GetFirstPoint();
2020 Int_t ri = GetRowNumber(xt);
2024 if (first<ri) first = ri;
2026 if (first<0) first=0;
2027 for (Int_t nr=first; nr<=rf; nr++) {
2028 // if ( (TMath::Abs(t.GetSnp())>0.95)) break;//patch 28 fev 06
2029 if (t.GetKinkIndexes()[0]<0){
2030 for (Int_t i=0;i<3;i++){
2031 Int_t index = t.GetKinkIndexes()[i];
2032 if (index==0) break;
2033 if (index>0) continue;
2034 index = TMath::Abs(index);
2035 AliKink * kink = (AliKink*)fEvent->GetKink(index-1);
2037 printf("PROBLEM\n");
2040 Int_t kinkrow = kink->GetTPCRow0()-2-Int_t(0.5/(0.05+kink->GetAngle(2)));
2042 AliExternalTrackParam paramm(t);
2043 kink->SetMother(paramm);
2044 kink->SetStatus(2,1);
2051 if (nr<fInnerSec->GetNRows())
2052 fSectors = fInnerSec;
2054 fSectors = fOuterSec;
2065 Float_t AliTPCtrackerMI::OverlapFactor(AliTPCseed * s1, AliTPCseed * s2, Int_t &sum1, Int_t & sum2)
2067 // overlapping factor
2073 Float_t dz2 =(s1->GetZ() - s2->GetZ());
2076 Float_t dy2 =TMath::Abs((s1->GetY() - s2->GetY()));
2078 Float_t distance = TMath::Sqrt(dz2+dy2);
2079 if (distance>4.) return 0; // if there are far away - not overlap - to reduce combinatorics
2082 Int_t firstpoint = TMath::Min(s1->GetFirstPoint(),s2->GetFirstPoint());
2083 Int_t lastpoint = TMath::Max(s1->GetLastPoint(),s2->GetLastPoint());
2088 if (firstpoint>lastpoint) {
2089 firstpoint =lastpoint;
2094 for (Int_t i=firstpoint-1;i<lastpoint+1;i++){
2095 if (s1->GetClusterIndex2(i)>0) sum1++;
2096 if (s2->GetClusterIndex2(i)>0) sum2++;
2097 if (s1->GetClusterIndex2(i)==s2->GetClusterIndex2(i) && s1->GetClusterIndex2(i)>0) {
2101 if (sum<5) return 0;
2103 Float_t summin = TMath::Min(sum1+1,sum2+1);
2104 Float_t ratio = (sum+1)/Float_t(summin);
2108 void AliTPCtrackerMI::SignShared(AliTPCseed * s1, AliTPCseed * s2)
2112 Float_t thetaCut = 0.2;//+10.*TMath::Sqrt(s1->GetSigmaTglZ()+ s2->GetSigmaTglZ());
2113 if (TMath::Abs(s1->GetTgl()-s2->GetTgl())>thetaCut) return;
2114 Float_t minCl = TMath::Min(s1->GetNumberOfClusters(),s2->GetNumberOfClusters());
2115 Int_t cutN0 = TMath::Max(5,TMath::Nint(0.1*minCl));
2120 //Int_t firstpoint = TMath::Max(s1->GetFirstPoint(),s2->GetFirstPoint());
2121 //Int_t lastpoint = TMath::Min(s1->GetLastPoint(),s2->GetLastPoint());
2122 Int_t firstpoint = 0;
2123 Int_t lastpoint = 160;
2125 // if (firstpoint>=lastpoint-5) return;;
2127 for (Int_t i=firstpoint;i<lastpoint;i++){
2128 // if ( (s1->GetClusterIndex2(i)&0xFFFF8FFF)==(s2->GetClusterIndex2(i)&0xFFFF8FFF) && s1->GetClusterIndex2(i)>0) {
2129 if ( (s1->GetClusterIndex2(i))==(s2->GetClusterIndex2(i)) && s1->GetClusterIndex2(i)>0) {
2131 s1->SetSharedMapBit(i, kTRUE);
2132 s2->SetSharedMapBit(i, kTRUE);
2134 if (s1->GetClusterIndex2(i)>0)
2135 s1->SetClusterMapBit(i, kTRUE);
2137 if (sumshared>cutN0){
2140 for (Int_t i=firstpoint;i<lastpoint;i++){
2141 // if ( (s1->GetClusterIndex2(i)&0xFFFF8FFF)==(s2->GetClusterIndex2(i)&0xFFFF8FFF) && s1->GetClusterIndex2(i)>0) {
2142 if ( (s1->GetClusterIndex2(i))==(s2->GetClusterIndex2(i)) && s1->GetClusterIndex2(i)>0) {
2143 AliTPCTrackerPoint *p1 = s1->GetTrackPoint(i);
2144 AliTPCTrackerPoint *p2 = s2->GetTrackPoint(i);;
2145 if (s1->IsActive()&&s2->IsActive()){
2146 p1->SetShared(kTRUE);
2147 p2->SetShared(kTRUE);
2153 if (sumshared>cutN0){
2154 for (Int_t i=0;i<4;i++){
2155 if (s1->GetOverlapLabel(3*i)==0){
2156 s1->SetOverlapLabel(3*i, s2->GetLabel());
2157 s1->SetOverlapLabel(3*i+1,sumshared);
2158 s1->SetOverlapLabel(3*i+2,s2->GetUniqueID());
2162 for (Int_t i=0;i<4;i++){
2163 if (s2->GetOverlapLabel(3*i)==0){
2164 s2->SetOverlapLabel(3*i, s1->GetLabel());
2165 s2->SetOverlapLabel(3*i+1,sumshared);
2166 s2->SetOverlapLabel(3*i+2,s1->GetUniqueID());
2173 void AliTPCtrackerMI::SignShared(TObjArray * arr)
2176 //sort trackss according sectors
2178 for (Int_t i=0; i<arr->GetEntriesFast(); i++) {
2179 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2181 //if (pt) RotateToLocal(pt);
2185 arr->Sort(); // sorting according relative sectors
2186 arr->Expand(arr->GetEntries());
2189 Int_t nseed=arr->GetEntriesFast();
2190 for (Int_t i=0; i<nseed; i++) {
2191 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2193 for (Int_t j=0;j<12;j++){
2194 pt->SetOverlapLabel(j,0);
2197 for (Int_t i=0; i<nseed; i++) {
2198 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2200 if (pt->GetRemoval()>10) continue;
2201 for (Int_t j=i+1; j<nseed; j++){
2202 AliTPCseed *pt2=(AliTPCseed*)arr->UncheckedAt(j);
2203 if (TMath::Abs(pt->GetRelativeSector()-pt2->GetRelativeSector())>1) continue;
2205 if (pt2->GetRemoval()<=10) {
2206 //if ( TMath::Abs(pt->GetRelativeSector()-pt2->GetRelativeSector())>0) break;
2214 void AliTPCtrackerMI::SortTracks(TObjArray * arr, Int_t mode) const
2217 //sort tracks in array according mode criteria
2218 Int_t nseed = arr->GetEntriesFast();
2219 for (Int_t i=0; i<nseed; i++) {
2220 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2231 void AliTPCtrackerMI::RemoveUsed2(TObjArray * arr, Float_t factor1, Float_t factor2, Int_t minimal)
2234 // Loop over all tracks and remove overlaped tracks (with lower quality)
2236 // 1. Unsign clusters
2237 // 2. Sort tracks according quality
2238 // Quality is defined by the number of cluster between first and last points
2240 // 3. Loop over tracks - decreasing quality order
2241 // a.) remove - If the fraction of shared cluster less than factor (1- n or 2)
2242 // b.) remove - If the minimal number of clusters less than minimal and not ITS
2243 // c.) if track accepted - sign clusters
2245 //Called in - AliTPCtrackerMI::Clusters2Tracks()
2246 // - AliTPCtrackerMI::PropagateBack()
2247 // - AliTPCtrackerMI::RefitInward()
2250 // factor1 - factor for constrained
2251 // factor2 - for non constrained tracks
2252 // if (Float_t(shared+1)/Float_t(found+1)>factor) - DELETE
2256 Int_t nseed = arr->GetEntriesFast();
2257 Float_t * quality = new Float_t[nseed];
2258 Int_t * indexes = new Int_t[nseed];
2262 for (Int_t i=0; i<nseed; i++) {
2263 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2268 pt->UpdatePoints(); //select first last max dens points
2269 Float_t * points = pt->GetPoints();
2270 if (points[3]<0.8) quality[i] =-1;
2271 quality[i] = (points[2]-points[0])+pt->GetNumberOfClusters();
2272 //prefer high momenta tracks if overlaps
2273 quality[i] *= TMath::Sqrt(TMath::Abs(pt->Pt())+0.5);
2275 TMath::Sort(nseed,quality,indexes);
2278 for (Int_t itrack=0; itrack<nseed; itrack++) {
2279 Int_t trackindex = indexes[itrack];
2280 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(trackindex);
2283 if (quality[trackindex]<0){
2284 delete arr->RemoveAt(trackindex);
2289 Int_t first = Int_t(pt->GetPoints()[0]);
2290 Int_t last = Int_t(pt->GetPoints()[2]);
2291 Double_t factor = (pt->GetBConstrain()) ? factor1: factor2;
2293 Int_t found,foundable,shared;
2294 pt->GetClusterStatistic(first,last, found, foundable,shared,kFALSE); // better to get statistic in "high-dens" region do't use full track as in line bellow
2295 // pt->GetClusterStatistic(0,160, found, foundable,shared,kFALSE);
2296 Bool_t itsgold =kFALSE;
2299 if (pt->GetESD()->GetITSclusters(dummy)>4) itsgold= kTRUE;
2303 if (Float_t(shared+1)/Float_t(found+1)>factor){
2304 if (pt->GetKinkIndexes()[0]!=0) continue; //don't remove tracks - part of the kinks
2305 if( AliTPCReconstructor::StreamLevel()>3){
2306 TTreeSRedirector &cstream = *fDebugStreamer;
2307 cstream<<"RemoveUsed"<<
2308 "iter="<<fIteration<<
2312 delete arr->RemoveAt(trackindex);
2315 if (pt->GetNumberOfClusters()<50&&(found-0.5*shared)<minimal){ //remove short tracks
2316 if (pt->GetKinkIndexes()[0]!=0) continue; //don't remove tracks - part of the kinks
2317 if( AliTPCReconstructor::StreamLevel()>3){
2318 TTreeSRedirector &cstream = *fDebugStreamer;
2319 cstream<<"RemoveShort"<<
2320 "iter="<<fIteration<<
2324 delete arr->RemoveAt(trackindex);
2330 //if (sharedfactor>0.4) continue;
2331 if (pt->GetKinkIndexes()[0]>0) continue;
2332 //Remove tracks with undefined properties - seems
2333 if (pt->GetSigmaY2()<kAlmost0) continue; // ? what is the origin ?
2335 for (Int_t i=first; i<last; i++) {
2336 Int_t index=pt->GetClusterIndex2(i);
2337 // if (index<0 || index&0x8000 ) continue;
2338 if (index<0 || index&0x8000 ) continue;
2339 AliTPCclusterMI *c= pt->GetClusterPointer(i);
2346 Info("RemoveUsed","\n*****\nNumber of good tracks after shared removal\t%d\n",fNtracks);
2352 void AliTPCtrackerMI::DumpClusters(Int_t iter, TObjArray *trackArray)
2355 // Dump clusters after reco
2356 // signed and unsigned cluster can be visualized
2357 // 1. Unsign all cluster
2358 // 2. Sign all used clusters
2361 Int_t nseed = trackArray->GetEntries();
2362 for (Int_t i=0; i<nseed; i++){
2363 AliTPCseed *pt=(AliTPCseed*)trackArray->UncheckedAt(i);
2367 Bool_t isKink=pt->GetKinkIndex(0)!=0;
2368 for (Int_t j=0; j<160; ++j) {
2369 Int_t index=pt->GetClusterIndex2(j);
2370 if (index<0) continue;
2371 AliTPCclusterMI *c= pt->GetClusterPointer(j);
2373 if (isKink) c->Use(100); // kink
2374 c->Use(10); // by default usage 10
2379 for (Int_t sec=0;sec<fkNIS;sec++){
2380 for (Int_t row=0;row<fInnerSec->GetNRows();row++){
2381 AliTPCclusterMI *cl = fInnerSec[sec][row].GetClusters1();
2382 for (Int_t icl =0;icl< fInnerSec[sec][row].GetN1();icl++){
2383 Float_t gx[3]; cl[icl].GetGlobalXYZ(gx);
2384 (*fDebugStreamer)<<"clDump"<<
2392 cl = fInnerSec[sec][row].GetClusters2();
2393 for (Int_t icl =0;icl< fInnerSec[sec][row].GetN2();icl++){
2394 Float_t gx[3]; cl[icl].GetGlobalXYZ(gx);
2395 (*fDebugStreamer)<<"clDump"<<
2406 for (Int_t sec=0;sec<fkNOS;sec++){
2407 for (Int_t row=0;row<fOuterSec->GetNRows();row++){
2408 AliTPCclusterMI *cl = fOuterSec[sec][row].GetClusters1();
2409 for (Int_t icl =0;icl< fOuterSec[sec][row].GetN1();icl++){
2410 Float_t gx[3]; cl[icl].GetGlobalXYZ(gx);
2411 (*fDebugStreamer)<<"clDump"<<
2419 cl = fOuterSec[sec][row].GetClusters2();
2420 for (Int_t icl =0;icl< fOuterSec[sec][row].GetN2();icl++){
2421 Float_t gx[3]; cl[icl].GetGlobalXYZ(gx);
2422 (*fDebugStreamer)<<"clDump"<<
2434 void AliTPCtrackerMI::UnsignClusters()
2437 // loop over all clusters and unsign them
2440 for (Int_t sec=0;sec<fkNIS;sec++){
2441 for (Int_t row=0;row<fInnerSec->GetNRows();row++){
2442 AliTPCclusterMI *cl = fInnerSec[sec][row].GetClusters1();
2443 for (Int_t icl =0;icl< fInnerSec[sec][row].GetN1();icl++)
2444 // if (cl[icl].IsUsed(10))
2446 cl = fInnerSec[sec][row].GetClusters2();
2447 for (Int_t icl =0;icl< fInnerSec[sec][row].GetN2();icl++)
2448 //if (cl[icl].IsUsed(10))
2453 for (Int_t sec=0;sec<fkNOS;sec++){
2454 for (Int_t row=0;row<fOuterSec->GetNRows();row++){
2455 AliTPCclusterMI *cl = fOuterSec[sec][row].GetClusters1();
2456 for (Int_t icl =0;icl< fOuterSec[sec][row].GetN1();icl++)
2457 //if (cl[icl].IsUsed(10))
2459 cl = fOuterSec[sec][row].GetClusters2();
2460 for (Int_t icl =0;icl< fOuterSec[sec][row].GetN2();icl++)
2461 //if (cl[icl].IsUsed(10))
2470 void AliTPCtrackerMI::SignClusters(const TObjArray * arr, Float_t fnumber, Float_t fdensity)
2473 //sign clusters to be "used"
2475 // snumber and sdensity sign number of sigmas - bellow mean value to be accepted
2476 // loop over "primaries"
2490 Int_t nseed = arr->GetEntriesFast();
2491 for (Int_t i=0; i<nseed; i++) {
2492 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2496 if (!(pt->IsActive())) continue;
2497 Float_t dens = pt->GetNumberOfClusters()/Float_t(pt->GetNFoundable());
2498 if ( (dens>0.7) && (pt->GetNumberOfClusters()>70)){
2500 sumdens2+= dens*dens;
2501 sumn += pt->GetNumberOfClusters();
2502 sumn2 += pt->GetNumberOfClusters()*pt->GetNumberOfClusters();
2503 Float_t chi2 = pt->GetChi2()/pt->GetNumberOfClusters();
2506 sumchi2 +=chi2*chi2;
2511 Float_t mdensity = 0.9;
2512 Float_t meann = 130;
2513 Float_t meanchi = 1;
2514 Float_t sdensity = 0.1;
2515 Float_t smeann = 10;
2516 Float_t smeanchi =0.4;
2520 mdensity = sumdens/sum;
2522 meanchi = sumchi/sum;
2524 sdensity = sumdens2/sum-mdensity*mdensity;
2526 sdensity = TMath::Sqrt(sdensity);
2530 smeann = sumn2/sum-meann*meann;
2532 smeann = TMath::Sqrt(smeann);
2536 smeanchi = sumchi2/sum - meanchi*meanchi;
2538 smeanchi = TMath::Sqrt(smeanchi);
2544 //REMOVE SHORT DELTAS or tracks going out of sensitive volume of TPC
2546 for (Int_t i=0; i<nseed; i++) {
2547 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2551 if (pt->GetBSigned()) continue;
2552 if (pt->GetBConstrain()) continue;
2553 //if (!(pt->IsActive())) continue;
2555 Int_t found,foundable,shared;
2556 pt->GetClusterStatistic(0,160,found, foundable,shared);
2557 if (shared/float(found)>0.3) {
2558 if (shared/float(found)>0.9 ){
2559 //delete arr->RemoveAt(i);
2564 Bool_t isok =kFALSE;
2565 if ( (pt->GetNShared()/pt->GetNumberOfClusters()<0.5) &&pt->GetNumberOfClusters()>60)
2567 if ((TMath::Abs(1/pt->GetC())<100.) && (pt->GetNShared()/pt->GetNumberOfClusters()<0.7))
2569 if (TMath::Abs(pt->GetZ()/pt->GetX())>1.1)
2571 if ( (TMath::Abs(pt->GetSnp()>0.7) && pt->GetD(0,0)>60.))
2575 for (Int_t j=0; j<160; ++j) {
2576 Int_t index=pt->GetClusterIndex2(j);
2577 if (index<0) continue;
2578 AliTPCclusterMI *c= pt->GetClusterPointer(j);
2580 //if (!(c->IsUsed(10))) c->Use();
2587 Double_t maxchi = meanchi+2.*smeanchi;
2589 for (Int_t i=0; i<nseed; i++) {
2590 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2594 //if (!(pt->IsActive())) continue;
2595 if (pt->GetBSigned()) continue;
2596 Double_t chi = pt->GetChi2()/pt->GetNumberOfClusters();
2597 if (chi>maxchi) continue;
2600 Float_t dens = pt->GetNumberOfClusters()/Float_t(pt->GetNFoundable());
2602 //sign only tracks with enoug big density at the beginning
2604 if ((pt->GetDensityFirst(40)<0.75) && pt->GetNumberOfClusters()<meann) continue;
2607 Double_t mindens = TMath::Max(double(mdensity-sdensity*fdensity*bfactor),0.65);
2608 Double_t minn = TMath::Max(Int_t(meann-fnumber*smeann*bfactor),50);
2610 // if (pt->fBConstrain) mindens = TMath::Max(mdensity-sdensity*fdensity*bfactor,0.65);
2611 if ( (pt->GetRemoval()==10) && (pt->GetSnp()>0.8)&&(dens>mindens))
2614 if ((dens>mindens && pt->GetNumberOfClusters()>minn) && chi<maxchi ){
2615 //Int_t noc=pt->GetNumberOfClusters();
2616 pt->SetBSigned(kTRUE);
2617 for (Int_t j=0; j<160; ++j) {
2619 Int_t index=pt->GetClusterIndex2(j);
2620 if (index<0) continue;
2621 AliTPCclusterMI *c= pt->GetClusterPointer(j);
2623 // if (!(c->IsUsed(10))) c->Use();
2628 // gLastCheck = nseed;
2636 void AliTPCtrackerMI::StopNotActive(const TObjArray * arr, Int_t row0, Float_t th0, Float_t th1, Float_t th2) const
2638 // stop not active tracks
2639 // take th1 as threshold for number of founded to number of foundable on last 10 active rows
2640 // take th2 as threshold for number of founded to number of foundable on last 20 active rows
2641 Int_t nseed = arr->GetEntriesFast();
2643 for (Int_t i=0; i<nseed; i++) {
2644 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2648 if (!(pt->IsActive())) continue;
2649 StopNotActive(pt,row0,th0, th1,th2);
2655 void AliTPCtrackerMI::StopNotActive(AliTPCseed * seed, Int_t row0, Float_t th0, Float_t th1,
2658 // stop not active tracks
2659 // take th1 as threshold for number of founded to number of foundable on last 10 active rows
2660 // take th2 as threshold for number of founded to number of foundable on last 20 active rows
2663 Int_t foundable = 0;
2664 Int_t maxindex = seed->GetLastPoint(); //last foundable row
2665 if (seed->GetNFoundable()*th0 > seed->GetNumberOfClusters()) {
2666 seed->Desactivate(10) ;
2670 for (Int_t i=row0; i<maxindex; i++){
2671 Int_t index = seed->GetClusterIndex2(i);
2672 if (index!=-1) foundable++;
2674 if (foundable<=30) sumgood1++;
2675 if (foundable<=50) {
2682 if (foundable>=30.){
2683 if (sumgood1<(th1*30.)) seed->Desactivate(10);
2686 if (sumgood2<(th2*50.)) seed->Desactivate(10);
2690 Int_t AliTPCtrackerMI::RefitInward(AliESDEvent *event)
2693 // back propagation of ESD tracks
2696 if (!event) return 0;
2697 const Int_t kMaxFriendTracks=2000;
2699 // extract correction object for multiplicity dependence of dEdx
2700 TObjArray * gainCalibArray = AliTPCcalibDB::Instance()->GetTimeGainSplinesRun(event->GetRunNumber());
2702 AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ;
2704 AliFatal("Tranformations not in RefitInward");
2707 transform->SetCurrentRecoParam((AliTPCRecoParam*)AliTPCReconstructor::GetRecoParam());
2708 const AliTPCRecoParam * recoParam = AliTPCcalibDB::Instance()->GetTransform()->GetCurrentRecoParam();
2709 Int_t nContribut = event->GetNumberOfTracks();
2710 TGraphErrors * graphMultDependenceDeDx = 0x0;
2711 if (recoParam && recoParam->GetUseMultiplicityCorrectionDedx() && gainCalibArray) {
2712 if (recoParam->GetUseTotCharge()) {
2713 graphMultDependenceDeDx = (TGraphErrors *) gainCalibArray->FindObject("TGRAPHERRORS_MEANQTOT_MULTIPLICITYDEPENDENCE_BEAM_ALL");
2715 graphMultDependenceDeDx = (TGraphErrors *) gainCalibArray->FindObject("TGRAPHERRORS_MEANQMAX_MULTIPLICITYDEPENDENCE_BEAM_ALL");
2721 //PrepareForProlongation(fSeeds,1);
2722 PropagateForward2(fSeeds);
2723 RemoveUsed2(fSeeds,0.4,0.4,20);
2725 TObjArray arraySeed(fSeeds->GetEntries());
2726 for (Int_t i=0;i<fSeeds->GetEntries();i++) {
2727 arraySeed.AddAt(fSeeds->At(i),i);
2729 SignShared(&arraySeed);
2730 // FindCurling(fSeeds, event,2); // find multi found tracks
2731 FindSplitted(fSeeds, event,2); // find multi found tracks
2732 if (AliTPCReconstructor::StreamLevel()>5) FindMultiMC(fSeeds, fEvent,2); // find multi found tracks
2735 Int_t nseed = fSeeds->GetEntriesFast();
2736 for (Int_t i=0;i<nseed;i++){
2737 AliTPCseed * seed = (AliTPCseed*) fSeeds->UncheckedAt(i);
2738 if (!seed) continue;
2739 if (seed->GetKinkIndex(0)>0) UpdateKinkQualityD(seed); // update quality informations for kinks
2740 AliESDtrack *esd=event->GetTrack(i);
2742 if (seed->GetNumberOfClusters()<60 && seed->GetNumberOfClusters()<(esd->GetTPCclusters(0) -5)*0.8){
2743 AliExternalTrackParam paramIn;
2744 AliExternalTrackParam paramOut;
2745 Int_t ncl = seed->RefitTrack(seed,¶mIn,¶mOut);
2746 if (AliTPCReconstructor::StreamLevel()>2) {
2747 (*fDebugStreamer)<<"RecoverIn"<<
2751 "pout.="<<¶mOut<<
2756 seed->Set(paramIn.GetX(),paramIn.GetAlpha(),paramIn.GetParameter(),paramIn.GetCovariance());
2757 seed->SetNumberOfClusters(ncl);
2761 seed->PropagateTo(fkParam->GetInnerRadiusLow());
2762 seed->UpdatePoints();
2763 AddCovariance(seed);
2765 seed->CookdEdx(0.02,0.6);
2766 CookLabel(seed,0.1); //For comparison only
2767 esd->SetTPCClusterMap(seed->GetClusterMap());
2768 esd->SetTPCSharedMap(seed->GetSharedMap());
2770 if (AliTPCReconstructor::StreamLevel()>1 && seed!=0) {
2771 TTreeSRedirector &cstream = *fDebugStreamer;
2778 if (seed->GetNumberOfClusters()>15){
2779 esd->UpdateTrackParams(seed,AliESDtrack::kTPCrefit);
2780 esd->SetTPCPoints(seed->GetPoints());
2781 esd->SetTPCPointsF(seed->GetNFoundable());
2782 Int_t ndedx = seed->GetNCDEDX(0);
2783 Float_t sdedx = seed->GetSDEDX(0);
2784 Float_t dedx = seed->GetdEdx();
2785 // apply mutliplicity dependent dEdx correction if available
2786 if (graphMultDependenceDeDx) {
2787 Double_t corrGain = AliTPCcalibDButil::EvalGraphConst(graphMultDependenceDeDx, nContribut);
2788 dedx += (1 - corrGain)*50.; // MIP is normalized to 50
2790 esd->SetTPCsignal(dedx, sdedx, ndedx);
2792 // add seed to the esd track in Calib level
2794 Bool_t storeFriend = gRandom->Rndm()<(kMaxFriendTracks)/Float_t(nseed);
2795 if (AliTPCReconstructor::StreamLevel()>0 &&storeFriend){
2796 AliTPCseed * seedCopy = new AliTPCseed(*seed, kTRUE);
2797 esd->AddCalibObject(seedCopy);
2802 //printf("problem\n");
2805 //FindKinks(fSeeds,event);
2806 if (AliTPCReconstructor::StreamLevel()>3) DumpClusters(2,fSeeds);
2807 Info("RefitInward","Number of refitted tracks %d",ntracks);
2812 Int_t AliTPCtrackerMI::PropagateBack(AliESDEvent *event)
2815 // back propagation of ESD tracks
2817 if (!event) return 0;
2821 PropagateBack(fSeeds);
2822 RemoveUsed2(fSeeds,0.4,0.4,20);
2823 //FindCurling(fSeeds, fEvent,1);
2824 FindSplitted(fSeeds, event,1); // find multi found tracks
2825 if (AliTPCReconstructor::StreamLevel()>5) FindMultiMC(fSeeds, fEvent,1); // find multi found tracks
2828 Int_t nseed = fSeeds->GetEntriesFast();
2830 for (Int_t i=0;i<nseed;i++){
2831 AliTPCseed * seed = (AliTPCseed*) fSeeds->UncheckedAt(i);
2832 if (!seed) continue;
2833 if (seed->GetKinkIndex(0)<0) UpdateKinkQualityM(seed); // update quality informations for kinks
2834 seed->UpdatePoints();
2835 AddCovariance(seed);
2836 AliESDtrack *esd=event->GetTrack(i);
2837 if (!esd) continue; //never happen
2838 if (seed->GetNumberOfClusters()<60 && seed->GetNumberOfClusters()<(esd->GetTPCclusters(0) -5)*0.8){
2839 AliExternalTrackParam paramIn;
2840 AliExternalTrackParam paramOut;
2841 Int_t ncl = seed->RefitTrack(seed,¶mIn,¶mOut);
2842 if (AliTPCReconstructor::StreamLevel()>2) {
2843 (*fDebugStreamer)<<"RecoverBack"<<
2847 "pout.="<<¶mOut<<
2852 seed->Set(paramOut.GetX(),paramOut.GetAlpha(),paramOut.GetParameter(),paramOut.GetCovariance());
2853 seed->SetNumberOfClusters(ncl);
2856 seed->CookdEdx(0.02,0.6);
2857 CookLabel(seed,0.1); //For comparison only
2858 if (seed->GetNumberOfClusters()>15){
2859 esd->UpdateTrackParams(seed,AliESDtrack::kTPCout);
2860 esd->SetTPCPoints(seed->GetPoints());
2861 esd->SetTPCPointsF(seed->GetNFoundable());
2862 Int_t ndedx = seed->GetNCDEDX(0);
2863 Float_t sdedx = seed->GetSDEDX(0);
2864 Float_t dedx = seed->GetdEdx();
2865 esd->SetTPCsignal(dedx, sdedx, ndedx);
2867 Int_t eventnumber = event->GetEventNumberInFile();// patch 28 fev 06
2868 // This is most likely NOT the event number you'd like to use. It has nothing to do with the 'real' event number
2869 if (AliTPCReconstructor::StreamLevel()>1 && esd) {
2870 (*fDebugStreamer)<<"Cback"<<
2873 "EventNrInFile="<<eventnumber<<
2878 if (AliTPCReconstructor::StreamLevel()>3) DumpClusters(1,fSeeds);
2879 //FindKinks(fSeeds,event);
2880 Info("PropagateBack","Number of back propagated tracks %d",ntracks);
2887 void AliTPCtrackerMI::DeleteSeeds()
2892 Int_t nseed = fSeeds->GetEntriesFast();
2893 for (Int_t i=0;i<nseed;i++){
2894 AliTPCseed * seed = (AliTPCseed*)fSeeds->At(i);
2895 if (seed) delete fSeeds->RemoveAt(i);
2902 void AliTPCtrackerMI::ReadSeeds(const AliESDEvent *const event, Int_t direction)
2905 //read seeds from the event
2907 Int_t nentr=event->GetNumberOfTracks();
2909 Info("ReadSeeds", "Number of ESD tracks: %d\n", nentr);
2914 fSeeds = new TObjArray(nentr);
2918 for (Int_t i=0; i<nentr; i++) {
2919 AliESDtrack *esd=event->GetTrack(i);
2920 ULong_t status=esd->GetStatus();
2921 if (!(status&AliESDtrack::kTPCin)) continue;
2922 AliTPCtrack t(*esd);
2923 t.SetNumberOfClusters(0);
2924 // AliTPCseed *seed = new AliTPCseed(t,t.GetAlpha());
2925 AliTPCseed *seed = new AliTPCseed(t/*,t.GetAlpha()*/);
2926 seed->SetUniqueID(esd->GetID());
2927 AddCovariance(seed); //add systematic ucertainty
2928 for (Int_t ikink=0;ikink<3;ikink++) {
2929 Int_t index = esd->GetKinkIndex(ikink);
2930 seed->GetKinkIndexes()[ikink] = index;
2931 if (index==0) continue;
2932 index = TMath::Abs(index);
2933 AliESDkink * kink = fEvent->GetKink(index-1);
2934 if (kink&&esd->GetKinkIndex(ikink)<0){
2935 if ((status & AliESDtrack::kTRDrefit) != 0) kink->SetStatus(1,2);
2936 if ((status & AliESDtrack::kITSout) != 0) kink->SetStatus(1,0);
2938 if (kink&&esd->GetKinkIndex(ikink)>0){
2939 if ((status & AliESDtrack::kTRDrefit) != 0) kink->SetStatus(1,6);
2940 if ((status & AliESDtrack::kITSout) != 0) kink->SetStatus(1,4);
2944 if (((status&AliESDtrack::kITSout)==0)&&(direction==1)) seed->ResetCovariance(10.);
2945 if ( direction ==2 &&(status & AliESDtrack::kTRDrefit) == 0 ) seed->ResetCovariance(10.);
2946 //if ( direction ==2 && ((status & AliESDtrack::kTPCout) == 0) ) {
2947 // fSeeds->AddAt(0,i);
2951 if ( direction ==2 &&(status & AliESDtrack::kTRDrefit) > 0 ) {
2952 Double_t par0[5],par1[5],alpha,x;
2953 esd->GetInnerExternalParameters(alpha,x,par0);
2954 esd->GetExternalParameters(x,par1);
2955 Double_t delta1 = TMath::Abs(par0[4]-par1[4])/(0.000000001+TMath::Abs(par0[4]+par1[4]));
2956 Double_t delta2 = TMath::Abs(par0[3]-par1[3]);
2958 if (esd->GetTRDncls()>0) trdchi2 = esd->GetTRDchi2()/esd->GetTRDncls();
2959 //reset covariance if suspicious
2960 if ( (delta1>0.1) || (delta2>0.006) ||trdchi2>7.)
2961 seed->ResetCovariance(10.);
2966 // rotate to the local coordinate system
2968 fSectors=fInnerSec; fN=fkNIS;
2969 Double_t alpha=seed->GetAlpha();
2970 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
2971 if (alpha < 0. ) alpha += 2.*TMath::Pi();
2972 Int_t ns=Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN;
2973 alpha =ns*fSectors->GetAlpha() + fSectors->GetAlphaShift();
2974 alpha-=seed->GetAlpha();
2975 if (alpha<-TMath::Pi()) alpha += 2*TMath::Pi();
2976 if (alpha>=TMath::Pi()) alpha -= 2*TMath::Pi();
2977 if (TMath::Abs(alpha) > 0.001) { // This should not happen normally
2978 AliWarning(Form("Rotating track over %f",alpha));
2979 if (!seed->Rotate(alpha)) {
2986 if (esd->GetKinkIndex(0)<=0){
2987 for (Int_t irow=0;irow<160;irow++){
2988 Int_t index = seed->GetClusterIndex2(irow);
2991 AliTPCclusterMI * cl = GetClusterMI(index);
2992 seed->SetClusterPointer(irow,cl);
2994 if ((index & 0x8000)==0){
2995 cl->Use(10); // accepted cluster
2997 cl->Use(6); // close cluster not accepted
3000 Info("ReadSeeds","Not found cluster");
3005 fSeeds->AddAt(seed,i);
3011 //_____________________________________________________________________________
3012 void AliTPCtrackerMI::MakeSeeds3(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, Float_t cuts[4],
3013 Float_t deltay, Int_t ddsec) {
3014 //-----------------------------------------------------------------
3015 // This function creates track seeds.
3016 // SEEDING WITH VERTEX CONSTRAIN
3017 //-----------------------------------------------------------------
3018 // cuts[0] - fP4 cut
3019 // cuts[1] - tan(phi) cut
3020 // cuts[2] - zvertex cut
3021 // cuts[3] - fP3 cut
3029 Double_t x[5], c[15];
3030 // Int_t di = i1-i2;
3032 AliTPCseed * seed = new AliTPCseed();
3033 Double_t alpha=fSectors->GetAlpha(), shift=fSectors->GetAlphaShift();
3034 Double_t cs=cos(alpha), sn=sin(alpha);
3036 // Double_t x1 =fOuterSec->GetX(i1);
3037 //Double_t xx2=fOuterSec->GetX(i2);
3039 Double_t x1 =GetXrow(i1);
3040 Double_t xx2=GetXrow(i2);
3042 Double_t x3=GetX(), y3=GetY(), z3=GetZ();
3044 Int_t imiddle = (i2+i1)/2; //middle pad row index
3045 Double_t xm = GetXrow(imiddle); // radius of middle pad-row
3046 const AliTPCtrackerRow& krm=GetRow(sec,imiddle); //middle pad -row
3050 const AliTPCtrackerRow& kr1=GetRow(ns,i1);
3051 Double_t ymax = GetMaxY(i1)-kr1.GetDeadZone()-1.5;
3052 Double_t ymaxm = GetMaxY(imiddle)-kr1.GetDeadZone()-1.5;
3055 // change cut on curvature if it can't reach this layer
3056 // maximal curvature set to reach it
3057 Double_t dvertexmax = TMath::Sqrt((x1-x3)*(x1-x3)+(ymax+5-y3)*(ymax+5-y3));
3058 if (dvertexmax*0.5*cuts[0]>0.85){
3059 cuts[0] = 0.85/(dvertexmax*0.5+1.);
3061 Double_t r2min = 1/(cuts[0]*cuts[0]); //minimal square of radius given by cut
3064 if (deltay>0) ddsec = 0;
3065 // loop over clusters
3066 for (Int_t is=0; is < kr1; is++) {
3068 if (kr1[is]->IsUsed(10)) continue;
3069 Double_t y1=kr1[is]->GetY(), z1=kr1[is]->GetZ();
3070 //if (TMath::Abs(y1)>ymax) continue;
3072 if (deltay>0 && TMath::Abs(ymax-TMath::Abs(y1))> deltay ) continue; // seed only at the edge
3074 // find possible directions
3075 Float_t anglez = (z1-z3)/(x1-x3);
3076 Float_t extraz = z1 - anglez*(x1-xx2); // extrapolated z
3079 //find rotation angles relative to line given by vertex and point 1
3080 Double_t dvertex2 = (x1-x3)*(x1-x3)+(y1-y3)*(y1-y3);
3081 Double_t dvertex = TMath::Sqrt(dvertex2);
3082 Double_t angle13 = TMath::ATan((y1-y3)/(x1-x3));
3083 Double_t cs13 = cos(-angle13), sn13 = sin(-angle13);
3086 // loop over 2 sectors
3092 Double_t dddz1=0; // direction of delta inclination in z axis
3099 for (Int_t dsec = dsec1; dsec<=dsec2;dsec++){
3100 Int_t sec2 = sec + dsec;
3102 // AliTPCtrackerRow& kr2 = fOuterSec[(sec2+fkNOS)%fkNOS][i2];
3103 //AliTPCtrackerRow& kr2m = fOuterSec[(sec2+fkNOS)%fkNOS][imiddle];
3104 AliTPCtrackerRow& kr2 = GetRow((sec2+fkNOS)%fkNOS,i2);
3105 AliTPCtrackerRow& kr2m = GetRow((sec2+fkNOS)%fkNOS,imiddle);
3106 Int_t index1 = TMath::Max(kr2.Find(extraz-0.6-dddz1*TMath::Abs(z1)*0.05)-1,0);
3107 Int_t index2 = TMath::Min(kr2.Find(extraz+0.6+dddz2*TMath::Abs(z1)*0.05)+1,kr2);
3109 // rotation angles to p1-p3
3110 Double_t cs13r = cos(-angle13+dsec*alpha)/dvertex, sn13r = sin(-angle13+dsec*alpha)/dvertex;
3111 Double_t x2, y2, z2;
3113 // Double_t dymax = maxangle*TMath::Abs(x1-xx2);
3116 Double_t dxx0 = (xx2-x3)*cs13r;
3117 Double_t dyy0 = (xx2-x3)*sn13r;
3118 for (Int_t js=index1; js < index2; js++) {
3119 const AliTPCclusterMI *kcl = kr2[js];
3120 if (kcl->IsUsed(10)) continue;
3122 //calcutate parameters
3124 Double_t yy0 = dyy0 +(kcl->GetY()-y3)*cs13r;
3126 if (TMath::Abs(yy0)<0.000001) continue;
3127 Double_t xx0 = dxx0 -(kcl->GetY()-y3)*sn13r;
3128 Double_t y0 = 0.5*(xx0*xx0+yy0*yy0-xx0)/yy0;
3129 Double_t r02 = (0.25+y0*y0)*dvertex2;
3130 //curvature (radius) cut
3131 if (r02<r2min) continue;
3135 Double_t c0 = 1/TMath::Sqrt(r02);
3139 //Double_t dfi0 = 2.*TMath::ASin(dvertex*c0*0.5);
3140 //Double_t dfi1 = 2.*TMath::ASin(TMath::Sqrt(yy0*yy0+(1-xx0)*(1-xx0))*dvertex*c0*0.5);
3141 Double_t dfi0 = 2.*AliTPCFastMath::FastAsin(dvertex*c0*0.5);
3142 Double_t dfi1 = 2.*AliTPCFastMath::FastAsin(TMath::Sqrt(yy0*yy0+(1-xx0)*(1-xx0))*dvertex*c0*0.5);
3145 Double_t z0 = kcl->GetZ();
3146 Double_t zzzz2 = z1-(z1-z3)*dfi1/dfi0;
3147 if (TMath::Abs(zzzz2-z0)>0.5) continue;
3150 Double_t dip = (z1-z0)*c0/dfi1;
3151 Double_t x0 = (0.5*cs13+y0*sn13)*dvertex*c0;
3162 x2= xx2*cs-y2*sn*dsec;
3163 y2=+xx2*sn*dsec+y2*cs;
3173 // do we have cluster at the middle ?
3175 GetProlongation(x1,xm,x,ym,zm);
3177 AliTPCclusterMI * cm=0;
3178 if (TMath::Abs(ym)-ymaxm<0){
3179 cm = krm.FindNearest2(ym,zm,1.0,0.6,dummy);
3180 if ((!cm) || (cm->IsUsed(10))) {
3185 // rotate y1 to system 0
3186 // get state vector in rotated system
3187 Double_t yr1 = (-0.5*sn13+y0*cs13)*dvertex*c0;
3188 Double_t xr2 = x0*cs+yr1*sn*dsec;
3189 Double_t xr[5]={kcl->GetY(),kcl->GetZ(), xr2, dip, c0};
3191 GetProlongation(xx2,xm,xr,ym,zm);
3192 if (TMath::Abs(ym)-ymaxm<0){
3193 cm = kr2m.FindNearest2(ym,zm,1.0,0.6,dummy);
3194 if ((!cm) || (cm->IsUsed(10))) {
3204 dym = ym - cm->GetY();
3205 dzm = zm - cm->GetZ();
3212 Double_t sy1=kr1[is]->GetSigmaY2()*2., sz1=kr1[is]->GetSigmaZ2()*2.;
3213 Double_t sy2=kcl->GetSigmaY2()*2., sz2=kcl->GetSigmaZ2()*2.;
3214 //Double_t sy3=400*3./12., sy=0.1, sz=0.1;
3215 Double_t sy3=25000*x[4]*x[4]+0.1, sy=0.1, sz=0.1;
3216 //Double_t sy3=25000*x[4]*x[4]*60+0.5, sy=0.1, sz=0.1;
3218 Double_t f40=(F1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
3219 Double_t f42=(F1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
3220 Double_t f43=(F1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
3221 Double_t f20=(F2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
3222 Double_t f22=(F2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
3223 Double_t f23=(F2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
3225 Double_t f30=(F3(x1,y1+sy,x2,y2,z1,z2)-x[3])/sy;
3226 Double_t f31=(F3(x1,y1,x2,y2,z1+sz,z2)-x[3])/sz;
3227 Double_t f32=(F3(x1,y1,x2,y2+sy,z1,z2)-x[3])/sy;
3228 Double_t f34=(F3(x1,y1,x2,y2,z1,z2+sz)-x[3])/sz;
3232 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3233 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3234 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3235 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3236 c[13]=f30*sy1*f40+f32*sy2*f42;
3237 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3239 // if (!BuildSeed(kr1[is],kcl,0,x1,x2,x3,x,c)) continue;
3241 UInt_t index=kr1.GetIndex(is);
3242 seed->~AliTPCseed(); // this does not set the pointer to 0...
3243 AliTPCseed *track=new(seed) AliTPCseed(x1, ns*alpha+shift, x, c, index);
3245 track->SetIsSeeding(kTRUE);
3246 track->SetSeed1(i1);
3247 track->SetSeed2(i2);
3248 track->SetSeedType(3);
3252 FollowProlongation(*track, (i1+i2)/2,1);
3253 Int_t foundable,found,shared;
3254 track->GetClusterStatistic((i1+i2)/2,i1, found, foundable, shared, kTRUE);
3255 if ((found<0.55*foundable) || shared>0.5*found || (track->GetSigmaY2()+track->GetSigmaZ2())>0.5){
3257 seed->~AliTPCseed();
3263 FollowProlongation(*track, i2,1);
3267 track->SetBConstrain(1);
3268 // track->fLastPoint = i1+fInnerSec->GetNRows(); // first cluster in track position
3269 track->SetLastPoint(i1); // first cluster in track position
3270 track->SetFirstPoint(track->GetLastPoint());
3272 if (track->GetNumberOfClusters()<(i1-i2)*0.5 ||
3273 track->GetNumberOfClusters() < track->GetNFoundable()*0.6 ||
3274 track->GetNShared()>0.4*track->GetNumberOfClusters() ) {
3276 seed->~AliTPCseed();
3280 // Z VERTEX CONDITION
3281 Double_t zv, bz=GetBz();
3282 if ( !track->GetZAt(0.,bz,zv) ) continue;
3283 if (TMath::Abs(zv-z3)>cuts[2]) {
3284 FollowProlongation(*track, TMath::Max(i2-20,0));
3285 if ( !track->GetZAt(0.,bz,zv) ) continue;
3286 if (TMath::Abs(zv-z3)>cuts[2]){
3287 FollowProlongation(*track, TMath::Max(i2-40,0));
3288 if ( !track->GetZAt(0.,bz,zv) ) continue;
3289 if (TMath::Abs(zv-z3)>cuts[2] &&(track->GetNumberOfClusters() > track->GetNFoundable()*0.7)){
3290 // make seed without constrain
3291 AliTPCseed * track2 = MakeSeed(track,0.2,0.5,1.);
3292 FollowProlongation(*track2, i2,1);
3293 track2->SetBConstrain(kFALSE);
3294 track2->SetSeedType(1);
3295 arr->AddLast(track2);
3297 seed->~AliTPCseed();
3302 seed->~AliTPCseed();
3309 track->SetSeedType(0);
3310 arr->AddLast(track);
3311 seed = new AliTPCseed;
3313 // don't consider other combinations
3314 if (track->GetNumberOfClusters() > track->GetNFoundable()*0.8)
3320 Info("MakeSeeds3","\nSeeding statistic:\t%d\t%d\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin,nout1,nout2);
3326 void AliTPCtrackerMI::MakeSeeds5(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, Float_t cuts[4],
3331 //-----------------------------------------------------------------
3332 // This function creates track seeds.
3333 //-----------------------------------------------------------------
3334 // cuts[0] - fP4 cut
3335 // cuts[1] - tan(phi) cut
3336 // cuts[2] - zvertex cut
3337 // cuts[3] - fP3 cut
3347 Double_t x[5], c[15];
3349 // make temporary seed
3350 AliTPCseed * seed = new AliTPCseed;
3351 Double_t alpha=fOuterSec->GetAlpha(), shift=fOuterSec->GetAlphaShift();
3352 // Double_t cs=cos(alpha), sn=sin(alpha);
3357 Double_t x1 = GetXrow(i1-1);
3358 const AliTPCtrackerRow& kr1=GetRow(sec,i1-1);
3359 Double_t y1max = GetMaxY(i1-1)-kr1.GetDeadZone()-1.5;
3361 Double_t x1p = GetXrow(i1);
3362 const AliTPCtrackerRow& kr1p=GetRow(sec,i1);
3364 Double_t x1m = GetXrow(i1-2);
3365 const AliTPCtrackerRow& kr1m=GetRow(sec,i1-2);
3368 //last 3 padrow for seeding
3369 AliTPCtrackerRow& kr3 = GetRow((sec+fkNOS)%fkNOS,i1-7);
3370 Double_t x3 = GetXrow(i1-7);
3371 // Double_t y3max= GetMaxY(i1-7)-kr3.fDeadZone-1.5;
3373 AliTPCtrackerRow& kr3p = GetRow((sec+fkNOS)%fkNOS,i1-6);
3374 Double_t x3p = GetXrow(i1-6);
3376 AliTPCtrackerRow& kr3m = GetRow((sec+fkNOS)%fkNOS,i1-8);
3377 Double_t x3m = GetXrow(i1-8);
3382 Int_t im = i1-4; //middle pad row index
3383 Double_t xm = GetXrow(im); // radius of middle pad-row
3384 const AliTPCtrackerRow& krm=GetRow(sec,im); //middle pad -row
3385 // Double_t ymmax = GetMaxY(im)-kr1.fDeadZone-1.5;
3388 Double_t deltax = x1-x3;
3389 Double_t dymax = deltax*cuts[1];
3390 Double_t dzmax = deltax*cuts[3];
3392 // loop over clusters
3393 for (Int_t is=0; is < kr1; is++) {
3395 if (kr1[is]->IsUsed(10)) continue;
3396 Double_t y1=kr1[is]->GetY(), z1=kr1[is]->GetZ();
3398 if (deltay>0 && TMath::Abs(y1max-TMath::Abs(y1))> deltay ) continue; // seed only at the edge
3400 Int_t index1 = TMath::Max(kr3.Find(z1-dzmax)-1,0);
3401 Int_t index2 = TMath::Min(kr3.Find(z1+dzmax)+1,kr3);
3407 for (Int_t js=index1; js < index2; js++) {
3408 const AliTPCclusterMI *kcl = kr3[js];
3409 if (kcl->IsUsed(10)) continue;
3411 // apply angular cuts
3412 if (TMath::Abs(y1-y3)>dymax) continue;
3415 if (TMath::Abs(z1-z3)>dzmax) continue;
3417 Double_t angley = (y1-y3)/(x1-x3);
3418 Double_t anglez = (z1-z3)/(x1-x3);
3420 Double_t erry = TMath::Abs(angley)*(x1-x1m)*0.5+0.5;
3421 Double_t errz = TMath::Abs(anglez)*(x1-x1m)*0.5+0.5;
3423 Double_t yyym = angley*(xm-x1)+y1;
3424 Double_t zzzm = anglez*(xm-x1)+z1;
3426 const AliTPCclusterMI *kcm = krm.FindNearest2(yyym,zzzm,erry,errz,index);
3428 if (kcm->IsUsed(10)) continue;
3430 erry = TMath::Abs(angley)*(x1-x1m)*0.4+0.5;
3431 errz = TMath::Abs(anglez)*(x1-x1m)*0.4+0.5;
3438 // look around first
3439 const AliTPCclusterMI *kc1m = kr1m.FindNearest2(angley*(x1m-x1)+y1,
3445 if (kc1m->IsUsed(10)) used++;
3447 const AliTPCclusterMI *kc1p = kr1p.FindNearest2(angley*(x1p-x1)+y1,
3453 if (kc1p->IsUsed(10)) used++;
3455 if (used>1) continue;
3456 if (found<1) continue;
3460 const AliTPCclusterMI *kc3m = kr3m.FindNearest2(angley*(x3m-x3)+y3,
3466 if (kc3m->IsUsed(10)) used++;
3470 const AliTPCclusterMI *kc3p = kr3p.FindNearest2(angley*(x3p-x3)+y3,
3476 if (kc3p->IsUsed(10)) used++;
3480 if (used>1) continue;
3481 if (found<3) continue;
3491 x[4]=F1(x1,y1,x2,y2,x3,y3);
3492 //if (TMath::Abs(x[4]) >= cuts[0]) continue;
3495 x[2]=F2(x1,y1,x2,y2,x3,y3);
3498 x[3]=F3n(x1,y1,x2,y2,z1,z2,x[4]);
3499 //if (TMath::Abs(x[3]) > cuts[3]) continue;
3503 Double_t sy1=0.1, sz1=0.1;
3504 Double_t sy2=0.1, sz2=0.1;
3505 Double_t sy3=0.1, sy=0.1, sz=0.1;
3507 Double_t f40=(F1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
3508 Double_t f42=(F1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
3509 Double_t f43=(F1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
3510 Double_t f20=(F2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
3511 Double_t f22=(F2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
3512 Double_t f23=(F2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
3514 Double_t f30=(F3(x1,y1+sy,x2,y2,z1,z2)-x[3])/sy;
3515 Double_t f31=(F3(x1,y1,x2,y2,z1+sz,z2)-x[3])/sz;
3516 Double_t f32=(F3(x1,y1,x2,y2+sy,z1,z2)-x[3])/sy;
3517 Double_t f34=(F3(x1,y1,x2,y2,z1,z2+sz)-x[3])/sz;
3521 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3522 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3523 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3524 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3525 c[13]=f30*sy1*f40+f32*sy2*f42;
3526 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3528 // if (!BuildSeed(kr1[is],kcl,0,x1,x2,x3,x,c)) continue;
3530 index=kr1.GetIndex(is);
3531 seed->~AliTPCseed();
3532 AliTPCseed *track=new(seed) AliTPCseed(x1, sec*alpha+shift, x, c, index);
3534 track->SetIsSeeding(kTRUE);
3537 FollowProlongation(*track, i1-7,1);
3538 if (track->GetNumberOfClusters() < track->GetNFoundable()*0.75 ||
3539 track->GetNShared()>0.6*track->GetNumberOfClusters() || ( track->GetSigmaY2()+ track->GetSigmaZ2())>0.6){
3541 seed->~AliTPCseed();
3547 FollowProlongation(*track, i2,1);
3548 track->SetBConstrain(0);
3549 track->SetLastPoint(i1+fInnerSec->GetNRows()); // first cluster in track position
3550 track->SetFirstPoint(track->GetLastPoint());
3552 if (track->GetNumberOfClusters()<(i1-i2)*0.5 ||
3553 track->GetNumberOfClusters()<track->GetNFoundable()*0.7 ||
3554 track->GetNShared()>2. || track->GetChi2()/track->GetNumberOfClusters()>6 || ( track->GetSigmaY2()+ track->GetSigmaZ2())>0.5 ) {
3556 seed->~AliTPCseed();
3561 FollowProlongation(*track, TMath::Max(i2-10,0),1);
3562 AliTPCseed * track2 = MakeSeed(track,0.2,0.5,0.9);
3563 FollowProlongation(*track2, i2,1);
3564 track2->SetBConstrain(kFALSE);
3565 track2->SetSeedType(4);
3566 arr->AddLast(track2);
3568 seed->~AliTPCseed();
3572 //arr->AddLast(track);
3573 //seed = new AliTPCseed;
3579 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);
3585 //_____________________________________________________________________________
3586 void AliTPCtrackerMI::MakeSeeds2(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, Float_t */*cuts[4]*/,
3587 Float_t deltay, Bool_t /*bconstrain*/) {
3588 //-----------------------------------------------------------------
3589 // This function creates track seeds - without vertex constraint
3590 //-----------------------------------------------------------------
3591 // cuts[0] - fP4 cut - not applied
3592 // cuts[1] - tan(phi) cut
3593 // cuts[2] - zvertex cut - not applied
3594 // cuts[3] - fP3 cut
3604 Double_t alpha=fOuterSec->GetAlpha(), shift=fOuterSec->GetAlphaShift();
3605 // Double_t cs=cos(alpha), sn=sin(alpha);
3606 Int_t row0 = (i1+i2)/2;
3607 Int_t drow = (i1-i2)/2;
3608 const AliTPCtrackerRow& kr0=fSectors[sec][row0];
3609 AliTPCtrackerRow * kr=0;
3611 AliTPCpolyTrack polytrack;
3612 Int_t nclusters=fSectors[sec][row0];
3613 AliTPCseed * seed = new AliTPCseed;
3618 for (Int_t is=0; is < nclusters; is++) { //LOOP over clusters
3620 Int_t nfoundable =0;
3621 for (Int_t iter =1; iter<2; iter++){ //iterations
3622 const AliTPCtrackerRow& krm=fSectors[sec][row0-iter];
3623 const AliTPCtrackerRow& krp=fSectors[sec][row0+iter];
3624 const AliTPCclusterMI * cl= kr0[is];
3626 if (cl->IsUsed(10)) {
3632 Double_t x = kr0.GetX();
3633 // Initialization of the polytrack
3638 Double_t y0= cl->GetY();
3639 Double_t z0= cl->GetZ();
3643 Double_t ymax = fSectors->GetMaxY(row0)-kr0.GetDeadZone()-1.5;
3644 if (deltay>0 && TMath::Abs(ymax-TMath::Abs(y0))> deltay ) continue; // seed only at the edge
3646 erry = (0.5)*cl->GetSigmaY2()/TMath::Sqrt(cl->GetQ())*6;
3647 errz = (0.5)*cl->GetSigmaZ2()/TMath::Sqrt(cl->GetQ())*6;
3648 polytrack.AddPoint(x,y0,z0,erry, errz);
3651 if (cl->IsUsed(10)) sumused++;
3654 Float_t roady = (5*TMath::Sqrt(cl->GetSigmaY2()+0.2)+1.)*iter;
3655 Float_t roadz = (5*TMath::Sqrt(cl->GetSigmaZ2()+0.2)+1.)*iter;
3658 AliTPCclusterMI * cl1 = krm.FindNearest(y0,z0,roady,roadz);
3659 if (cl1 && TMath::Abs(ymax-TMath::Abs(y0))) {
3660 erry = (0.5)*cl1->GetSigmaY2()/TMath::Sqrt(cl1->GetQ())*3;
3661 errz = (0.5)*cl1->GetSigmaZ2()/TMath::Sqrt(cl1->GetQ())*3;
3662 if (cl1->IsUsed(10)) sumused++;
3663 polytrack.AddPoint(x,cl1->GetY(),cl1->GetZ(),erry,errz);
3667 AliTPCclusterMI * cl2 = krp.FindNearest(y0,z0,roady,roadz);
3669 erry = (0.5)*cl2->GetSigmaY2()/TMath::Sqrt(cl2->GetQ())*3;
3670 errz = (0.5)*cl2->GetSigmaZ2()/TMath::Sqrt(cl2->GetQ())*3;
3671 if (cl2->IsUsed(10)) sumused++;
3672 polytrack.AddPoint(x,cl2->GetY(),cl2->GetZ(),erry,errz);
3675 if (sumused>0) continue;
3677 polytrack.UpdateParameters();
3683 nfoundable = polytrack.GetN();
3684 nfound = nfoundable;
3686 for (Int_t ddrow = iter+1; ddrow<drow;ddrow++){
3687 Float_t maxdist = 0.8*(1.+3./(ddrow));
3688 for (Int_t delta = -1;delta<=1;delta+=2){
3689 Int_t row = row0+ddrow*delta;
3690 kr = &(fSectors[sec][row]);
3691 Double_t xn = kr->GetX();
3692 Double_t ymax1 = fSectors->GetMaxY(row)-kr->GetDeadZone()-1.5;
3693 polytrack.GetFitPoint(xn,yn,zn);
3694 if (TMath::Abs(yn)>ymax1) continue;
3696 AliTPCclusterMI * cln = kr->FindNearest(yn,zn,roady,roadz);
3698 Float_t dist = TMath::Sqrt( (yn-cln->GetY())*(yn-cln->GetY())+(zn-cln->GetZ())*(zn-cln->GetZ()));
3701 erry = (dist+0.3)*cln->GetSigmaY2()/TMath::Sqrt(cln->GetQ())*(1.+1./(ddrow));
3702 errz = (dist+0.3)*cln->GetSigmaZ2()/TMath::Sqrt(cln->GetQ())*(1.+1./(ddrow));
3703 if (cln->IsUsed(10)) {
3704 // printf("used\n");
3712 polytrack.AddPoint(xn,cln->GetY(),cln->GetZ(),erry, errz);
3717 if ( (sumused>3) || (sumused>0.5*nfound) || (nfound<0.6*nfoundable)) break;
3718 polytrack.UpdateParameters();
3721 if ( (sumused>3) || (sumused>0.5*nfound)) {
3722 //printf("sumused %d\n",sumused);
3727 polytrack.GetFitDerivation(kr0.GetX(),dy,dz);
3728 AliTPCpolyTrack track2;
3730 polytrack.Refit(track2,0.5+TMath::Abs(dy)*0.3,0.4+TMath::Abs(dz)*0.3);
3731 if (track2.GetN()<0.5*nfoundable) continue;
3734 if ((nfound>0.6*nfoundable) &&( nfoundable>0.4*(i1-i2))) {
3736 // test seed with and without constrain
3737 for (Int_t constrain=0; constrain<=0;constrain++){
3738 // add polytrack candidate
3740 Double_t x[5], c[15];
3741 Double_t x1,x2,x3,y1,y2,y3,z1,z2,z3;
3742 track2.GetBoundaries(x3,x1);
3744 track2.GetFitPoint(x1,y1,z1);
3745 track2.GetFitPoint(x2,y2,z2);
3746 track2.GetFitPoint(x3,y3,z3);
3748 //is track pointing to the vertex ?
3751 polytrack.GetFitPoint(x0,y0,z0);
3764 x[4]=F1(x1,y1,x2,y2,x3,y3);
3766 // if (TMath::Abs(x[4]) >= cuts[0]) continue; //
3767 x[2]=F2(x1,y1,x2,y2,x3,y3);
3769 //if (TMath::Abs(x[4]*x1-x[2]) >= cuts[1]) continue;
3770 //x[3]=F3(x1,y1,x2,y2,z1,z2);
3771 x[3]=F3n(x1,y1,x3,y3,z1,z3,x[4]);
3772 //if (TMath::Abs(x[3]) > cuts[3]) continue;
3775 Double_t sy =0.1, sz =0.1;
3776 Double_t sy1=0.02, sz1=0.02;
3777 Double_t sy2=0.02, sz2=0.02;
3781 sy3=25000*x[4]*x[4]+0.1, sy=0.1, sz=0.1;
3784 Double_t f40=(F1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
3785 Double_t f42=(F1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
3786 Double_t f43=(F1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
3787 Double_t f20=(F2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
3788 Double_t f22=(F2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
3789 Double_t f23=(F2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
3791 Double_t f30=(F3(x1,y1+sy,x3,y3,z1,z3)-x[3])/sy;
3792 Double_t f31=(F3(x1,y1,x3,y3,z1+sz,z3)-x[3])/sz;
3793 Double_t f32=(F3(x1,y1,x3,y3+sy,z1,z3)-x[3])/sy;
3794 Double_t f34=(F3(x1,y1,x3,y3,z1,z3+sz)-x[3])/sz;
3799 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3800 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3801 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3802 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3803 c[13]=f30*sy1*f40+f32*sy2*f42;
3804 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3806 //Int_t row1 = fSectors->GetRowNumber(x1);
3807 Int_t row1 = GetRowNumber(x1);
3811 seed->~AliTPCseed();
3812 AliTPCseed *track=new(seed) AliTPCseed(x1,sec*alpha+shift,x,c,index);
3813 track->SetIsSeeding(kTRUE);
3814 Int_t rc=FollowProlongation(*track, i2);
3815 if (constrain) track->SetBConstrain(1);
3817 track->SetBConstrain(0);
3818 track->SetLastPoint(row1+fInnerSec->GetNRows()); // first cluster in track position
3819 track->SetFirstPoint(track->GetLastPoint());
3821 if (rc==0 || track->GetNumberOfClusters()<(i1-i2)*0.5 ||
3822 track->GetNumberOfClusters() < track->GetNFoundable()*0.6 ||
3823 track->GetNShared()>0.4*track->GetNumberOfClusters()) {
3826 seed->~AliTPCseed();
3829 arr->AddLast(track);
3830 seed = new AliTPCseed;
3834 } // if accepted seed
3837 Info("MakeSeeds2","\nSeeding statiistic:\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin3);
3843 AliTPCseed *AliTPCtrackerMI::MakeSeed(AliTPCseed *const track, Float_t r0, Float_t r1, Float_t r2)
3847 //reseed using track points
3848 Int_t p0 = int(r0*track->GetNumberOfClusters()); // point 0
3849 Int_t p1 = int(r1*track->GetNumberOfClusters());
3850 Int_t p2 = int(r2*track->GetNumberOfClusters()); // last point
3852 Double_t x0[3],x1[3],x2[3];
3853 for (Int_t i=0;i<3;i++){
3859 // find track position at given ratio of the length
3860 Int_t sec0=0, sec1=0, sec2=0;
3863 for (Int_t i=0;i<160;i++){
3864 if (track->GetClusterPointer(i)){
3866 AliTPCTrackerPoint *trpoint =track->GetTrackPoint(i);
3867 if ( (index<p0) || x0[0]<0 ){
3868 if (trpoint->GetX()>1){
3869 clindex = track->GetClusterIndex2(i);
3871 x0[0] = trpoint->GetX();
3872 x0[1] = trpoint->GetY();
3873 x0[2] = trpoint->GetZ();
3874 sec0 = ((clindex&0xff000000)>>24)%18;
3879 if ( (index<p1) &&(trpoint->GetX()>1)){
3880 clindex = track->GetClusterIndex2(i);
3882 x1[0] = trpoint->GetX();
3883 x1[1] = trpoint->GetY();
3884 x1[2] = trpoint->GetZ();
3885 sec1 = ((clindex&0xff000000)>>24)%18;
3888 if ( (index<p2) &&(trpoint->GetX()>1)){
3889 clindex = track->GetClusterIndex2(i);
3891 x2[0] = trpoint->GetX();
3892 x2[1] = trpoint->GetY();
3893 x2[2] = trpoint->GetZ();
3894 sec2 = ((clindex&0xff000000)>>24)%18;
3901 Double_t alpha, cs,sn, xx2,yy2;
3903 alpha = (sec1-sec2)*fSectors->GetAlpha();
3904 cs = TMath::Cos(alpha);
3905 sn = TMath::Sin(alpha);
3906 xx2= x1[0]*cs-x1[1]*sn;
3907 yy2= x1[0]*sn+x1[1]*cs;
3911 alpha = (sec0-sec2)*fSectors->GetAlpha();
3912 cs = TMath::Cos(alpha);
3913 sn = TMath::Sin(alpha);
3914 xx2= x0[0]*cs-x0[1]*sn;
3915 yy2= x0[0]*sn+x0[1]*cs;
3921 Double_t x[5],c[15];
3925 x[4]=F1(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]);
3926 // if (x[4]>1) return 0;
3927 x[2]=F2(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]);
3928 x[3]=F3n(x2[0],x2[1],x0[0],x0[1],x2[2],x0[2],x[4]);
3929 //if (TMath::Abs(x[3]) > 2.2) return 0;
3930 //if (TMath::Abs(x[2]) > 1.99) return 0;
3932 Double_t sy =0.1, sz =0.1;
3934 Double_t sy1=0.02+track->GetSigmaY2(), sz1=0.02+track->GetSigmaZ2();
3935 Double_t sy2=0.01+track->GetSigmaY2(), sz2=0.01+track->GetSigmaZ2();
3936 Double_t sy3=0.01+track->GetSigmaY2();
3938 Double_t f40=(F1(x2[0],x2[1]+sy,x1[0],x1[1],x0[0],x0[1])-x[4])/sy;
3939 Double_t f42=(F1(x2[0],x2[1],x1[0],x1[1]+sy,x0[0],x0[1])-x[4])/sy;
3940 Double_t f43=(F1(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]+sy)-x[4])/sy;
3941 Double_t f20=(F2(x2[0],x2[1]+sy,x1[0],x1[1],x0[0],x0[1])-x[2])/sy;
3942 Double_t f22=(F2(x2[0],x2[1],x1[0],x1[1]+sy,x0[0],x0[1])-x[2])/sy;
3943 Double_t f23=(F2(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]+sy)-x[2])/sy;
3945 Double_t f30=(F3(x2[0],x2[1]+sy,x0[0],x0[1],x2[2],x0[2])-x[3])/sy;
3946 Double_t f31=(F3(x2[0],x2[1],x0[0],x0[1],x2[2]+sz,x0[2])-x[3])/sz;
3947 Double_t f32=(F3(x2[0],x2[1],x0[0],x0[1]+sy,x2[2],x0[2])-x[3])/sy;
3948 Double_t f34=(F3(x2[0],x2[1],x0[0],x0[1],x2[2],x0[2]+sz)-x[3])/sz;
3953 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3954 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3955 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3956 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3957 c[13]=f30*sy1*f40+f32*sy2*f42;
3958 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3960 // Int_t row1 = fSectors->GetRowNumber(x2[0]);
3961 AliTPCseed *seed=new AliTPCseed(x2[0], sec2*fSectors->GetAlpha()+fSectors->GetAlphaShift(), x, c, 0);
3962 // Double_t y0,z0,y1,z1, y2,z2;
3963 //seed->GetProlongation(x0[0],y0,z0);
3964 // seed->GetProlongation(x1[0],y1,z1);
3965 //seed->GetProlongation(x2[0],y2,z2);
3967 seed->SetLastPoint(pp2);
3968 seed->SetFirstPoint(pp2);
3975 AliTPCseed *AliTPCtrackerMI::ReSeed(const AliTPCseed *track, Float_t r0, Float_t r1, Float_t r2)
3979 //reseed using founded clusters
3981 // Find the number of clusters
3982 Int_t nclusters = 0;
3983 for (Int_t irow=0;irow<160;irow++){
3984 if (track->GetClusterIndex(irow)>0) nclusters++;
3988 ipos[0] = TMath::Max(int(r0*nclusters),0); // point 0 cluster
3989 ipos[1] = TMath::Min(int(r1*nclusters),nclusters-1); //
3990 ipos[2] = TMath::Min(int(r2*nclusters),nclusters-1); // last point
3993 Double_t xyz[3][3]={{0}};
3994 Int_t row[3]={0},sec[3]={0,0,0};
3996 // find track row position at given ratio of the length
3998 for (Int_t irow=0;irow<160;irow++){
3999 if (track->GetClusterIndex2(irow)<0) continue;
4001 for (Int_t ipoint=0;ipoint<3;ipoint++){
4002 if (index<=ipos[ipoint]) row[ipoint] = irow;
4006 //Get cluster and sector position
4007 for (Int_t ipoint=0;ipoint<3;ipoint++){
4008 Int_t clindex = track->GetClusterIndex2(row[ipoint]);
4009 AliTPCclusterMI * cl = GetClusterMI(clindex);
4012 // AliTPCclusterMI * cl = GetClusterMI(clindex);
4015 sec[ipoint] = ((clindex&0xff000000)>>24)%18;
4016 xyz[ipoint][0] = GetXrow(row[ipoint]);
4017 xyz[ipoint][1] = cl->GetY();
4018 xyz[ipoint][2] = cl->GetZ();
4022 // Calculate seed state vector and covariance matrix
4024 Double_t alpha, cs,sn, xx2,yy2;
4026 alpha = (sec[1]-sec[2])*fSectors->GetAlpha();
4027 cs = TMath::Cos(alpha);
4028 sn = TMath::Sin(alpha);
4029 xx2= xyz[1][0]*cs-xyz[1][1]*sn;
4030 yy2= xyz[1][0]*sn+xyz[1][1]*cs;
4034 alpha = (sec[0]-sec[2])*fSectors->GetAlpha();
4035 cs = TMath::Cos(alpha);
4036 sn = TMath::Sin(alpha);
4037 xx2= xyz[0][0]*cs-xyz[0][1]*sn;
4038 yy2= xyz[0][0]*sn+xyz[0][1]*cs;
4044 Double_t x[5],c[15];
4048 x[4]=F1(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
4049 x[2]=F2(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
4050 x[3]=F3n(xyz[2][0],xyz[2][1],xyz[0][0],xyz[0][1],xyz[2][2],xyz[0][2],x[4]);
4052 Double_t sy =0.1, sz =0.1;
4054 Double_t sy1=0.2, sz1=0.2;
4055 Double_t sy2=0.2, sz2=0.2;
4058 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;
4059 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;
4060 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;
4061 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;
4062 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;
4063 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;
4065 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;
4066 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;
4067 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;
4068 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;
4073 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
4074 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
4075 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
4076 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
4077 c[13]=f30*sy1*f40+f32*sy2*f42;
4078 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
4080 // Int_t row1 = fSectors->GetRowNumber(xyz[2][0]);
4081 AliTPCseed *seed=new AliTPCseed(xyz[2][0], sec[2]*fSectors->GetAlpha()+fSectors->GetAlphaShift(), x, c, 0);
4082 seed->SetLastPoint(row[2]);
4083 seed->SetFirstPoint(row[2]);
4088 AliTPCseed *AliTPCtrackerMI::ReSeed(AliTPCseed *track,Int_t r0, Bool_t forward)
4092 //reseed using founded clusters
4095 Int_t row[3]={0,0,0};
4096 Int_t sec[3]={0,0,0};
4098 // forward direction
4100 for (Int_t irow=r0;irow<160;irow++){
4101 if (track->GetClusterIndex(irow)>0){
4106 for (Int_t irow=160;irow>r0;irow--){
4107 if (track->GetClusterIndex(irow)>0){
4112 for (Int_t irow=row[2]-15;irow>row[0];irow--){
4113 if (track->GetClusterIndex(irow)>0){
4121 for (Int_t irow=0;irow<r0;irow++){
4122 if (track->GetClusterIndex(irow)>0){
4127 for (Int_t irow=r0;irow>0;irow--){
4128 if (track->GetClusterIndex(irow)>0){
4133 for (Int_t irow=row[2]-15;irow>row[0];irow--){
4134 if (track->GetClusterIndex(irow)>0){
4141 if ((row[2]-row[0])<20) return 0;
4142 if (row[1]==0) return 0;
4145 //Get cluster and sector position
4146 for (Int_t ipoint=0;ipoint<3;ipoint++){
4147 Int_t clindex = track->GetClusterIndex2(row[ipoint]);
4148 AliTPCclusterMI * cl = GetClusterMI(clindex);
4151 // AliTPCclusterMI * cl = GetClusterMI(clindex);
4154 sec[ipoint] = ((clindex&0xff000000)>>24)%18;
4155 xyz[ipoint][0] = GetXrow(row[ipoint]);
4156 AliTPCTrackerPoint * point = track->GetTrackPoint(row[ipoint]);
4157 if (point&&ipoint<2){
4159 xyz[ipoint][1] = point->GetY();
4160 xyz[ipoint][2] = point->GetZ();
4163 xyz[ipoint][1] = cl->GetY();
4164 xyz[ipoint][2] = cl->GetZ();
4171 // Calculate seed state vector and covariance matrix
4173 Double_t alpha, cs,sn, xx2,yy2;
4175 alpha = (sec[1]-sec[2])*fSectors->GetAlpha();
4176 cs = TMath::Cos(alpha);
4177 sn = TMath::Sin(alpha);
4178 xx2= xyz[1][0]*cs-xyz[1][1]*sn;
4179 yy2= xyz[1][0]*sn+xyz[1][1]*cs;
4183 alpha = (sec[0]-sec[2])*fSectors->GetAlpha();
4184 cs = TMath::Cos(alpha);
4185 sn = TMath::Sin(alpha);
4186 xx2= xyz[0][0]*cs-xyz[0][1]*sn;
4187 yy2= xyz[0][0]*sn+xyz[0][1]*cs;
4193 Double_t x[5],c[15];
4197 x[4]=F1(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
4198 x[2]=F2(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
4199 x[3]=F3n(xyz[2][0],xyz[2][1],xyz[0][0],xyz[0][1],xyz[2][2],xyz[0][2],x[4]);
4201 Double_t sy =0.1, sz =0.1;
4203 Double_t sy1=0.2, sz1=0.2;
4204 Double_t sy2=0.2, sz2=0.2;
4207 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;
4208 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;
4209 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;
4210 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;
4211 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;
4212 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;
4214 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;
4215 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;
4216 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;
4217 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;
4222 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
4223 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
4224 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
4225 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
4226 c[13]=f30*sy1*f40+f32*sy2*f42;
4227 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
4229 // Int_t row1 = fSectors->GetRowNumber(xyz[2][0]);
4230 AliTPCseed *seed=new AliTPCseed(xyz[2][0], sec[2]*fSectors->GetAlpha()+fSectors->GetAlphaShift(), x, c, 0);
4231 seed->SetLastPoint(row[2]);
4232 seed->SetFirstPoint(row[2]);
4233 for (Int_t i=row[0];i<row[2];i++){
4234 seed->SetClusterIndex(i, track->GetClusterIndex(i));
4242 void AliTPCtrackerMI::FindMultiMC(const TObjArray * array, AliESDEvent */*esd*/, Int_t iter)
4245 // find multi tracks - THIS FUNCTION IS ONLY FOR DEBUG PURPOSES
4247 // Use AliTPCReconstructor::StreamLevel()>2 if you want to tune parameters - cuts
4249 // Two reasons to have multiple find tracks
4250 // 1. Curling tracks can be find more than once
4251 // 2. Splitted tracks
4252 // a.) Multiple seeding to increase tracking efficiency - (~ 100% reached)
4253 // b.) Edge effect on the sector boundaries
4256 // Algorithm done in 2 phases - because of CPU consumption
4257 // it is n^2 algorithm - for lead-lead 20000x20000 combination are investigated
4259 // Algorihm for curling tracks sign:
4260 // 1 phase -makes a very rough fast cuts to minimize combinatorics
4261 // a.) opposite sign
4262 // b.) one of the tracks - not pointing to the primary vertex -
4263 // c.) delta tan(theta)
4265 // 2 phase - calculates DCA between tracks - time consument
4270 // General cuts - for splitted tracks and for curling tracks
4272 const Float_t kMaxdPhi = 0.2; // maximal distance in phi
4274 // Curling tracks cuts
4279 Int_t nentries = array->GetEntriesFast();
4280 AliHelix *helixes = new AliHelix[nentries];
4281 Float_t *xm = new Float_t[nentries];
4282 Float_t *dz0 = new Float_t[nentries];
4283 Float_t *dz1 = new Float_t[nentries];
4289 // Find track COG in x direction - point with best defined parameters
4291 for (Int_t i=0;i<nentries;i++){
4292 AliTPCseed* track = (AliTPCseed*)array->At(i);
4293 if (!track) continue;
4294 track->SetCircular(0);
4295 new (&helixes[i]) AliHelix(*track);
4299 track->GetDZ(GetX(),GetY(),GetZ(),GetBz(),dz);
4302 for (Int_t icl=0; icl<160; icl++){
4303 AliTPCclusterMI * cl = track->GetClusterPointer(icl);
4309 if (ncl>0) xm[i]/=Float_t(ncl);
4312 for (Int_t i0=0;i0<nentries;i0++){
4313 AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
4314 if (!track0) continue;
4315 Float_t xc0 = helixes[i0].GetHelix(6);
4316 Float_t yc0 = helixes[i0].GetHelix(7);
4317 Float_t r0 = helixes[i0].GetHelix(8);
4318 Float_t rc0 = TMath::Sqrt(xc0*xc0+yc0*yc0);
4319 Float_t fi0 = TMath::ATan2(yc0,xc0);
4321 for (Int_t i1=i0+1;i1<nentries;i1++){
4322 AliTPCseed * track1 = (AliTPCseed*)array->At(i1);
4323 if (!track1) continue;
4324 Int_t lab0=track0->GetLabel();
4325 Int_t lab1=track1->GetLabel();
4326 if (TMath::Abs(lab0)!=TMath::Abs(lab1)) continue;
4328 Float_t xc1 = helixes[i1].GetHelix(6);
4329 Float_t yc1 = helixes[i1].GetHelix(7);
4330 Float_t r1 = helixes[i1].GetHelix(8);
4331 Float_t rc1 = TMath::Sqrt(xc1*xc1+yc1*yc1);
4332 Float_t fi1 = TMath::ATan2(yc1,xc1);
4334 Float_t dfi = fi0-fi1;
4337 if (dfi>1.5*TMath::Pi()) dfi-=TMath::Pi(); // take care about edge effect
4338 if (dfi<-1.5*TMath::Pi()) dfi+=TMath::Pi(); //
4339 if (TMath::Abs(dfi)>kMaxdPhi&&helixes[i0].GetHelix(4)*helixes[i1].GetHelix(4)<0){
4341 // if short tracks with undefined sign
4342 fi1 = -TMath::ATan2(yc1,-xc1);
4345 Float_t dtheta = TMath::Abs(track0->GetTgl()-track1->GetTgl())<TMath::Abs(track0->GetTgl()+track1->GetTgl())? track0->GetTgl()-track1->GetTgl():track0->GetTgl()+track1->GetTgl();
4348 // debug stream to tune "fast cuts"
4350 Double_t dist[3]; // distance at X
4351 Double_t mdist[3]={0,0,0}; // mean distance X+-40cm
4352 track0->GetDistance(track1,0.5*(xm[i0]+xm[i1])-40.,dist,AliTracker::GetBz());
4353 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[i]);
4354 track0->GetDistance(track1,0.5*(xm[i0]+xm[i1])+40.,dist,AliTracker::GetBz());
4355 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[i]);
4356 track0->GetDistance(track1,0.5*(xm[i0]+xm[i1]),dist,AliTracker::GetBz());
4357 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[i]);
4358 for (Int_t i=0;i<3;i++) mdist[i]*=0.33333;
4362 for (Int_t icl=0; icl<160; icl++){
4363 AliTPCclusterMI * cl0 = track0->GetClusterPointer(icl);
4364 AliTPCclusterMI * cl1 = track1->GetClusterPointer(icl);
4367 if (cl0==cl1) sums++;
4371 if (AliTPCReconstructor::StreamLevel()>5) {
4372 TTreeSRedirector &cstream = *fDebugStreamer;
4377 "Tr0.="<<track0<< // seed0
4378 "Tr1.="<<track1<< // seed1
4379 "h0.="<<&helixes[i0]<<
4380 "h1.="<<&helixes[i1]<<
4382 "sum="<<sum<< //the sum of rows with cl in both
4383 "sums="<<sums<< //the sum of shared clusters
4384 "xm0="<<xm[i0]<< // the center of track
4385 "xm1="<<xm[i1]<< // the x center of track
4386 // General cut variables
4387 "dfi="<<dfi<< // distance in fi angle
4388 "dtheta="<<dtheta<< // distance int theta angle
4394 "dist0="<<dist[0]<< //distance x
4395 "dist1="<<dist[1]<< //distance y
4396 "dist2="<<dist[2]<< //distance z
4397 "mdist0="<<mdist[0]<< //distance x
4398 "mdist1="<<mdist[1]<< //distance y
4399 "mdist2="<<mdist[2]<< //distance z
4415 if (AliTPCReconstructor::StreamLevel()>1) {
4416 AliInfo("Time for curling tracks removal DEBUGGING MC");
4423 void AliTPCtrackerMI::FindSplitted(TObjArray * array, AliESDEvent */*esd*/, Int_t /*iter*/){
4425 // Find Splitted tracks and remove the one with worst quality
4426 // Corresponding debug streamer to tune selections - "Splitted2"
4428 // 0. Sort tracks according quility
4429 // 1. Propagate the tracks to the reference radius
4430 // 2. Double_t loop to select close tracks (only to speed up process)
4431 // 3. Calculate cluster overlap ratio - and remove the track if bigger than a threshold
4432 // 4. Delete temporary parameters
4434 const Double_t xref=GetXrow(63); // reference radius -IROC/OROC boundary
4436 const Double_t kCutP1=10; // delta Z cut 10 cm
4437 const Double_t kCutP2=0.15; // delta snp(fi) cut 0.15
4438 const Double_t kCutP3=0.15; // delta tgl(theta) cut 0.15
4439 const Double_t kCutAlpha=0.15; // delta alpha cut
4440 Int_t firstpoint = 0;
4441 Int_t lastpoint = 160;
4443 Int_t nentries = array->GetEntriesFast();
4444 AliExternalTrackParam *params = new AliExternalTrackParam[nentries];
4450 //0. Sort tracks according quality
4451 //1. Propagate the ext. param to reference radius
4452 Int_t nseed = array->GetEntriesFast();
4453 if (nseed<=0) return;
4454 Float_t * quality = new Float_t[nseed];
4455 Int_t * indexes = new Int_t[nseed];
4456 for (Int_t i=0; i<nseed; i++) {
4457 AliTPCseed *pt=(AliTPCseed*)array->UncheckedAt(i);
4462 pt->UpdatePoints(); //select first last max dens points
4463 Float_t * points = pt->GetPoints();
4464 if (points[3]<0.8) quality[i] =-1;
4465 quality[i] = (points[2]-points[0])+pt->GetNumberOfClusters();
4466 //prefer high momenta tracks if overlaps
4467 quality[i] *= TMath::Sqrt(TMath::Abs(pt->Pt())+0.5);
4469 AliTracker::PropagateTrackToBxByBz(&(params[i]),xref,pt->GetMass(),5.,kTRUE);
4470 AliTracker::PropagateTrackToBxByBz(&(params[i]),xref,pt->GetMass(),1.,kTRUE);
4472 TMath::Sort(nseed,quality,indexes);
4474 // 3. Loop over pair of tracks
4476 for (Int_t i0=0; i0<nseed; i0++) {
4477 Int_t index0=indexes[i0];
4478 if (!(array->UncheckedAt(index0))) continue;
4479 AliTPCseed *s1 = (AliTPCseed*)array->UncheckedAt(index0);
4480 if (!s1->IsActive()) continue;
4481 AliExternalTrackParam &par0=params[index0];
4482 for (Int_t i1=i0+1; i1<nseed; i1++) {
4483 Int_t index1=indexes[i1];
4484 if (!(array->UncheckedAt(index1))) continue;
4485 AliTPCseed *s2 = (AliTPCseed*)array->UncheckedAt(index1);
4486 if (!s2->IsActive()) continue;
4487 if (s2->GetKinkIndexes()[0]!=0)
4488 if (s2->GetKinkIndexes()[0] == -s1->GetKinkIndexes()[0]) continue;
4489 AliExternalTrackParam &par1=params[index1];
4490 if (TMath::Abs(par0.GetParameter()[3]-par1.GetParameter()[3])>kCutP3) continue;
4491 if (TMath::Abs(par0.GetParameter()[1]-par1.GetParameter()[1])>kCutP1) continue;
4492 if (TMath::Abs(par0.GetParameter()[2]-par1.GetParameter()[2])>kCutP2) continue;
4493 Double_t dAlpha= TMath::Abs(par0.GetAlpha()-par1.GetAlpha());
4494 if (dAlpha>TMath::Pi()) dAlpha-=TMath::Pi();
4495 if (TMath::Abs(dAlpha)>kCutAlpha) continue;
4500 Int_t firstShared=lastpoint, lastShared=firstpoint;
4501 Int_t firstRow=lastpoint, lastRow=firstpoint;
4503 for (Int_t i=firstpoint;i<lastpoint;i++){
4504 if (s1->GetClusterIndex2(i)>0) nall0++;
4505 if (s2->GetClusterIndex2(i)>0) nall1++;
4506 if (s1->GetClusterIndex2(i)>0 && s2->GetClusterIndex2(i)>0) {
4507 if (i<firstRow) firstRow=i;
4508 if (i>lastRow) lastRow=i;
4510 if ( (s1->GetClusterIndex2(i))==(s2->GetClusterIndex2(i)) && s1->GetClusterIndex2(i)>0) {
4511 if (i<firstShared) firstShared=i;
4512 if (i>lastShared) lastShared=i;
4516 Double_t ratio0 = Float_t(sumShared)/Float_t(TMath::Min(nall0+1,nall1+1));
4517 Double_t ratio1 = Float_t(sumShared)/Float_t(TMath::Max(nall0+1,nall1+1));
4519 if( AliTPCReconstructor::StreamLevel()>1){
4520 TTreeSRedirector &cstream = *fDebugStreamer;
4521 Int_t n0=s1->GetNumberOfClusters();
4522 Int_t n1=s2->GetNumberOfClusters();
4523 Int_t n0F=s1->GetNFoundable();
4524 Int_t n1F=s2->GetNFoundable();
4525 Int_t lab0=s1->GetLabel();
4526 Int_t lab1=s2->GetLabel();
4528 cstream<<"Splitted2"<<
4529 "iter="<<fIteration<<
4530 "lab0="<<lab0<< // MC label if exist
4531 "lab1="<<lab1<< // MC label if exist
4534 "ratio0="<<ratio0<< // shared ratio
4535 "ratio1="<<ratio1<< // shared ratio
4536 "p0.="<<&par0<< // track parameters
4538 "s0.="<<s1<< // full seed
4540 "n0="<<n0<< // number of clusters track 0
4541 "n1="<<n1<< // number of clusters track 1
4542 "nall0="<<nall0<< // number of clusters track 0
4543 "nall1="<<nall1<< // number of clusters track 1
4544 "n0F="<<n0F<< // number of findable
4545 "n1F="<<n1F<< // number of findable
4546 "shared="<<sumShared<< // number of shared clusters
4547 "firstS="<<firstShared<< // first and the last shared row
4548 "lastS="<<lastShared<<
4549 "firstRow="<<firstRow<< // first and the last row with cluster
4550 "lastRow="<<lastRow<< //
4554 // remove track with lower quality
4556 if (ratio0>AliTPCReconstructor::GetRecoParam()->GetCutSharedClusters(0) ||
4557 ratio1>AliTPCReconstructor::GetRecoParam()->GetCutSharedClusters(1)){
4561 delete array->RemoveAt(index1);
4566 // 4. Delete temporary array
4576 void AliTPCtrackerMI::FindCurling(const TObjArray * array, AliESDEvent */*esd*/, Int_t iter)
4579 // find Curling tracks
4580 // Use AliTPCReconstructor::StreamLevel()>1 if you want to tune parameters - cuts
4583 // Algorithm done in 2 phases - because of CPU consumption
4584 // it is n^2 algorithm - for lead-lead 20000x20000 combination are investigated
4585 // see detal in MC part what can be used to cut
4589 const Float_t kMaxC = 400; // maximal curvature to of the track
4590 const Float_t kMaxdTheta = 0.15; // maximal distance in theta
4591 const Float_t kMaxdPhi = 0.15; // maximal distance in phi
4592 const Float_t kPtRatio = 0.3; // ratio between pt
4593 const Float_t kMinDCAR = 2.; // distance to the primary vertex in r - see cpipe cut
4596 // Curling tracks cuts
4599 const Float_t kMaxDeltaRMax = 40; // distance in outer radius
4600 const Float_t kMaxDeltaRMin = 5.; // distance in lower radius - see cpipe cut
4601 const Float_t kMinAngle = 2.9; // angle between tracks
4602 const Float_t kMaxDist = 5; // biggest distance
4604 // The cuts can be tuned using the "MC information stored in Multi tree ==> see FindMultiMC
4607 TCut csign("csign","Tr0.fP[4]*Tr1.fP[4]<0"); //opposite sign
4608 TCut cmax("cmax","abs(Tr0.GetC())>1/400");
4609 TCut cda("cda","sqrt(dtheta^2+dfi^2)<0.15");
4610 TCut ccratio("ccratio","abs((Tr0.fP[4]+Tr1.fP[4])/(abs(Tr0.fP[4])+abs(Tr1.fP[4])))<0.3");
4611 TCut cpipe("cpipe", "min(abs(r0-rc0),abs(r1-rc1))>5");
4613 TCut cdrmax("cdrmax","abs(abs(rc0+r0)-abs(rc1+r1))<40")
4614 TCut cdrmin("cdrmin","abs(abs(rc0+r0)-abs(rc1+r1))<10")
4616 Multi->Draw("dfi","iter==0"+csign+cmax+cda+ccratio); ~94% of curling tracks fulfill
4617 Multi->Draw("min(abs(r0-rc0),abs(r1-rc1))","iter==0&&abs(lab1)==abs(lab0)"+csign+cmax+cda+ccratio+cpipe+cdrmin+cdrmax); //80%
4619 Curling2->Draw("dfi","iter==0&&abs(lab0)==abs(lab1)"+csign+cmax+cdtheta+cdfi+ccratio)
4625 Int_t nentries = array->GetEntriesFast();
4626 AliHelix *helixes = new AliHelix[nentries];
4627 for (Int_t i=0;i<nentries;i++){
4628 AliTPCseed* track = (AliTPCseed*)array->At(i);
4629 if (!track) continue;
4630 track->SetCircular(0);
4631 new (&helixes[i]) AliHelix(*track);
4637 Double_t phase[2][2]={{0,0},{0,0}},radius[2]={0,0};
4643 for (Int_t i0=0;i0<nentries;i0++){
4644 AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
4645 if (!track0) continue;
4646 if (TMath::Abs(track0->GetC())<1/kMaxC) continue;
4647 Float_t xc0 = helixes[i0].GetHelix(6);
4648 Float_t yc0 = helixes[i0].GetHelix(7);
4649 Float_t r0 = helixes[i0].GetHelix(8);
4650 Float_t rc0 = TMath::Sqrt(xc0*xc0+yc0*yc0);
4651 Float_t fi0 = TMath::ATan2(yc0,xc0);
4653 for (Int_t i1=i0+1;i1<nentries;i1++){
4654 AliTPCseed * track1 = (AliTPCseed*)array->At(i1);
4655 if (!track1) continue;
4656 if (TMath::Abs(track1->GetC())<1/kMaxC) continue;
4657 Float_t xc1 = helixes[i1].GetHelix(6);
4658 Float_t yc1 = helixes[i1].GetHelix(7);
4659 Float_t r1 = helixes[i1].GetHelix(8);
4660 Float_t rc1 = TMath::Sqrt(xc1*xc1+yc1*yc1);
4661 Float_t fi1 = TMath::ATan2(yc1,xc1);
4663 Float_t dfi = fi0-fi1;
4666 if (dfi>1.5*TMath::Pi()) dfi-=TMath::Pi(); // take care about edge effect
4667 if (dfi<-1.5*TMath::Pi()) dfi+=TMath::Pi(); //
4668 Float_t dtheta = TMath::Abs(track0->GetTgl()-track1->GetTgl())<TMath::Abs(track0->GetTgl()+track1->GetTgl())? track0->GetTgl()-track1->GetTgl():track0->GetTgl()+track1->GetTgl();
4672 if (track0->GetBConstrain()&&track1->GetBConstrain()) continue; // not constrained
4673 if (track1->GetSigned1Pt()*track0->GetSigned1Pt()>0) continue; // not the same sign
4674 if ( TMath::Abs(track1->GetTgl()+track0->GetTgl())>kMaxdTheta) continue; //distance in the Theta
4675 if ( TMath::Abs(dfi)>kMaxdPhi) continue; //distance in phi
4676 if ( TMath::Sqrt(dfi*dfi+dtheta*dtheta)>kMaxdPhi) continue; //common angular offset
4678 Float_t pt0 = track0->GetSignedPt();
4679 Float_t pt1 = track1->GetSignedPt();
4680 if ((TMath::Abs(pt0+pt1)/(TMath::Abs(pt0)+TMath::Abs(pt1)))>kPtRatio) continue;
4681 if ((iter==1) && TMath::Abs(TMath::Abs(rc0+r0)-TMath::Abs(rc1+r1))>kMaxDeltaRMax) continue;
4682 if ((iter!=1) &&TMath::Abs(TMath::Abs(rc0-r0)-TMath::Abs(rc1-r1))>kMaxDeltaRMin) continue;
4683 if (TMath::Min(TMath::Abs(rc0-r0),TMath::Abs(rc1-r1))<kMinDCAR) continue;
4686 // Now find closest approach
4690 Int_t npoints = helixes[i0].GetRPHIintersections(helixes[i1], phase, radius,10);
4691 if (npoints==0) continue;
4692 helixes[i0].GetClosestPhases(helixes[i1], phase);
4696 Double_t hangles[3];
4697 helixes[i0].Evaluate(phase[0][0],xyz0);
4698 helixes[i1].Evaluate(phase[0][1],xyz1);
4700 helixes[i0].GetAngle(phase[0][0],helixes[i1],phase[0][1],hangles);
4701 Double_t deltah[2],deltabest;
4702 if (TMath::Abs(hangles[2])<kMinAngle) continue;
4706 helixes[i0].ParabolicDCA(helixes[i1],phase[0][0],phase[0][1],radius[0],deltah[0],2);
4708 helixes[i0].ParabolicDCA(helixes[i1],phase[1][0],phase[1][1],radius[1],deltah[1],2);
4709 if (deltah[1]<deltah[0]) ibest=1;
4711 deltabest = TMath::Sqrt(deltah[ibest]);
4712 helixes[i0].Evaluate(phase[ibest][0],xyz0);
4713 helixes[i1].Evaluate(phase[ibest][1],xyz1);
4714 helixes[i0].GetAngle(phase[ibest][0],helixes[i1],phase[ibest][1],hangles);
4715 Double_t radiusbest = TMath::Sqrt(radius[ibest]);
4717 if (deltabest>kMaxDist) continue;
4718 // if (mindcar+mindcaz<40 && (TMath::Abs(hangles[2])<kMinAngle ||deltabest>3)) continue;
4719 Bool_t sign =kFALSE;
4720 if (hangles[2]>kMinAngle) sign =kTRUE;
4723 // circular[i0] = kTRUE;
4724 // circular[i1] = kTRUE;
4725 if (track0->OneOverPt()<track1->OneOverPt()){
4726 track0->SetCircular(track0->GetCircular()+1);
4727 track1->SetCircular(track1->GetCircular()+2);
4730 track1->SetCircular(track1->GetCircular()+1);
4731 track0->SetCircular(track0->GetCircular()+2);
4734 if (AliTPCReconstructor::StreamLevel()>2){
4736 //debug stream to tune "fine" cuts
4737 Int_t lab0=track0->GetLabel();
4738 Int_t lab1=track1->GetLabel();
4739 TTreeSRedirector &cstream = *fDebugStreamer;
4740 cstream<<"Curling2"<<
4756 "npoints="<<npoints<<
4757 "hangles0="<<hangles[0]<<
4758 "hangles1="<<hangles[1]<<
4759 "hangles2="<<hangles[2]<<
4762 "radius="<<radiusbest<<
4763 "deltabest="<<deltabest<<
4764 "phase0="<<phase[ibest][0]<<
4765 "phase1="<<phase[ibest][1]<<
4773 if (AliTPCReconstructor::StreamLevel()>1) {
4774 AliInfo("Time for curling tracks removal");
4783 void AliTPCtrackerMI::FindKinks(TObjArray * array, AliESDEvent *esd)
4790 TObjArray *kinks= new TObjArray(10000);
4791 // TObjArray *v0s= new TObjArray(10000);
4792 Int_t nentries = array->GetEntriesFast();
4793 AliHelix *helixes = new AliHelix[nentries];
4794 Int_t *sign = new Int_t[nentries];
4795 Int_t *nclusters = new Int_t[nentries];
4796 Float_t *alpha = new Float_t[nentries];
4797 AliKink *kink = new AliKink();
4798 Int_t * usage = new Int_t[nentries];
4799 Float_t *zm = new Float_t[nentries];
4800 Float_t *z0 = new Float_t[nentries];
4801 Float_t *fim = new Float_t[nentries];
4802 Float_t *shared = new Float_t[nentries];
4803 Bool_t *circular = new Bool_t[nentries];
4804 Float_t *dca = new Float_t[nentries];
4805 //const AliESDVertex * primvertex = esd->GetVertex();
4807 // nentries = array->GetEntriesFast();
4812 for (Int_t i=0;i<nentries;i++){
4815 AliTPCseed* track = (AliTPCseed*)array->At(i);
4816 if (!track) continue;
4817 track->SetCircular(0);
4819 track->UpdatePoints();
4820 if (( track->GetPoints()[2]- track->GetPoints()[0])>5 && track->GetPoints()[3]>0.8){
4822 nclusters[i]=track->GetNumberOfClusters();
4823 alpha[i] = track->GetAlpha();
4824 new (&helixes[i]) AliHelix(*track);
4826 helixes[i].Evaluate(0,xyz);
4827 sign[i] = (track->GetC()>0) ? -1:1;
4830 if (track->GetProlongation(x,y,z)){
4832 fim[i] = alpha[i]+TMath::ATan2(y,x);
4835 zm[i] = track->GetZ();
4839 circular[i]= kFALSE;
4840 if (track->GetProlongation(0,y,z)) z0[i] = z;
4841 dca[i] = track->GetD(0,0);
4847 Int_t ncandidates =0;
4850 Double_t phase[2][2]={{0,0},{0,0}},radius[2]={0,0};
4853 // Find circling track
4855 for (Int_t i0=0;i0<nentries;i0++){
4856 AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
4857 if (!track0) continue;
4858 if (track0->GetNumberOfClusters()<40) continue;
4859 if (TMath::Abs(1./track0->GetC())>200) continue;
4860 for (Int_t i1=i0+1;i1<nentries;i1++){
4861 AliTPCseed * track1 = (AliTPCseed*)array->At(i1);
4862 if (!track1) continue;
4863 if (track1->GetNumberOfClusters()<40) continue;
4864 if ( TMath::Abs(track1->GetTgl()+track0->GetTgl())>0.1) continue;
4865 if (track0->GetBConstrain()&&track1->GetBConstrain()) continue;
4866 if (TMath::Abs(1./track1->GetC())>200) continue;
4867 if (track1->GetSigned1Pt()*track0->GetSigned1Pt()>0) continue;
4868 if (track1->GetTgl()*track0->GetTgl()>0) continue;
4869 if (TMath::Max(TMath::Abs(1./track0->GetC()),TMath::Abs(1./track1->GetC()))>190) continue;
4870 if (track0->GetBConstrain()&&track1->OneOverPt()<track0->OneOverPt()) continue; //returning - lower momenta
4871 if (track1->GetBConstrain()&&track0->OneOverPt()<track1->OneOverPt()) continue; //returning - lower momenta
4873 Float_t mindcar = TMath::Min(TMath::Abs(dca[i0]),TMath::Abs(dca[i1]));
4874 if (mindcar<5) continue;
4875 Float_t mindcaz = TMath::Min(TMath::Abs(z0[i0]-GetZ()),TMath::Abs(z0[i1]-GetZ()));
4876 if (mindcaz<5) continue;
4877 if (mindcar+mindcaz<20) continue;
4880 Float_t xc0 = helixes[i0].GetHelix(6);
4881 Float_t yc0 = helixes[i0].GetHelix(7);
4882 Float_t r0 = helixes[i0].GetHelix(8);
4883 Float_t xc1 = helixes[i1].GetHelix(6);
4884 Float_t yc1 = helixes[i1].GetHelix(7);
4885 Float_t r1 = helixes[i1].GetHelix(8);
4887 Float_t rmean = (r0+r1)*0.5;
4888 Float_t delta =TMath::Sqrt((xc1-xc0)*(xc1-xc0)+(yc1-yc0)*(yc1-yc0));
4889 //if (delta>30) continue;
4890 if (delta>rmean*0.25) continue;
4891 if (TMath::Abs(r0-r1)/rmean>0.3) continue;
4893 Int_t npoints = helixes[i0].GetRPHIintersections(helixes[i1], phase, radius,10);
4894 if (npoints==0) continue;
4895 helixes[i0].GetClosestPhases(helixes[i1], phase);
4899 Double_t hangles[3];
4900 helixes[i0].Evaluate(phase[0][0],xyz0);
4901 helixes[i1].Evaluate(phase[0][1],xyz1);
4903 helixes[i0].GetAngle(phase[0][0],helixes[i1],phase[0][1],hangles);
4904 Double_t deltah[2],deltabest;
4905 if (hangles[2]<2.8) continue;
4908 helixes[i0].ParabolicDCA(helixes[i1],phase[0][0],phase[0][1],radius[0],deltah[0],2);
4910 helixes[i0].ParabolicDCA(helixes[i1],phase[1][0],phase[1][1],radius[1],deltah[1],2);
4911 if (deltah[1]<deltah[0]) ibest=1;
4913 deltabest = TMath::Sqrt(deltah[ibest]);
4914 helixes[i0].Evaluate(phase[ibest][0],xyz0);
4915 helixes[i1].Evaluate(phase[ibest][1],xyz1);
4916 helixes[i0].GetAngle(phase[ibest][0],helixes[i1],phase[ibest][1],hangles);
4917 Double_t radiusbest = TMath::Sqrt(radius[ibest]);
4919 if (deltabest>6) continue;
4920 if (mindcar+mindcaz<40 && (hangles[2]<3.12||deltabest>3)) continue;
4921 Bool_t lsign =kFALSE;
4922 if (hangles[2]>3.06) lsign =kTRUE;
4925 circular[i0] = kTRUE;
4926 circular[i1] = kTRUE;
4927 if (track0->OneOverPt()<track1->OneOverPt()){
4928 track0->SetCircular(track0->GetCircular()+1);
4929 track1->SetCircular(track1->GetCircular()+2);
4932 track1->SetCircular(track1->GetCircular()+1);
4933 track0->SetCircular(track0->GetCircular()+2);
4936 if (lsign&&AliTPCReconstructor::StreamLevel()>1){
4938 Int_t lab0=track0->GetLabel();
4939 Int_t lab1=track1->GetLabel();
4940 TTreeSRedirector &cstream = *fDebugStreamer;
4941 cstream<<"Curling"<<
4948 "mindcar="<<mindcar<<
4949 "mindcaz="<<mindcaz<<
4952 "npoints="<<npoints<<
4953 "hangles0="<<hangles[0]<<
4954 "hangles2="<<hangles[2]<<
4959 "radius="<<radiusbest<<
4960 "deltabest="<<deltabest<<
4961 "phase0="<<phase[ibest][0]<<
4962 "phase1="<<phase[ibest][1]<<
4972 for (Int_t i =0;i<nentries;i++){
4973 if (sign[i]==0) continue;
4974 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
4981 Double_t cradius0 = 40*40;
4982 Double_t cradius1 = 270*270;
4985 Double_t cdist3=0.55;
4986 for (Int_t j =i+1;j<nentries;j++){
4988 if (sign[j]*sign[i]<1) continue;
4989 if ( (nclusters[i]+nclusters[j])>200) continue;
4990 if ( (nclusters[i]+nclusters[j])<80) continue;
4991 if ( TMath::Abs(zm[i]-zm[j])>60.) continue;
4992 if ( TMath::Abs(fim[i]-fim[j])>0.6 && TMath::Abs(fim[i]-fim[j])<5.7 ) continue;
4993 //AliTPCseed * track1 = (AliTPCseed*)array->At(j); Double_t phase[2][2],radius[2];
4994 Int_t npoints = helixes[i].GetRPHIintersections(helixes[j], phase, radius,20);
4995 if (npoints<1) continue;
4998 if (radius[0]<cradius0||radius[0]>cradius1) continue;
5001 if ( (radius[0]<cradius0||radius[0]>cradius1) && (radius[1]<cradius0||radius[1]>cradius1) ) continue;
5004 Double_t delta1=10000,delta2=10000;
5005 // cuts on the intersection radius
5006 helixes[i].LinearDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta1);
5007 if (radius[0]<20&&delta1<1) continue; //intersection at vertex
5008 if (radius[0]<10&&delta1<3) continue; //intersection at vertex
5010 helixes[i].LinearDCA(helixes[j],phase[1][0],phase[1][1],radius[1],delta2);
5011 if (radius[1]<20&&delta2<1) continue; //intersection at vertex
5012 if (radius[1]<10&&delta2<3) continue; //intersection at vertex
5015 Double_t distance1 = TMath::Min(delta1,delta2);
5016 if (distance1>cdist1) continue; // cut on DCA linear approximation
5018 npoints = helixes[i].GetRPHIintersections(helixes[j], phase, radius,20);
5019 helixes[i].ParabolicDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta1);
5020 if (radius[0]<20&&delta1<1) continue; //intersection at vertex
5021 if (radius[0]<10&&delta1<3) continue; //intersection at vertex
5024 helixes[i].ParabolicDCA(helixes[j],phase[1][0],phase[1][1],radius[1],delta2);
5025 if (radius[1]<20&&delta2<1) continue; //intersection at vertex
5026 if (radius[1]<10&&delta2<3) continue; //intersection at vertex
5028 distance1 = TMath::Min(delta1,delta2);
5031 rkink = TMath::Sqrt(radius[0]);
5034 rkink = TMath::Sqrt(radius[1]);
5036 if (distance1>cdist2) continue;
5039 AliTPCseed * track1 = (AliTPCseed*)array->At(j);
5042 Int_t row0 = GetRowNumber(rkink);
5043 if (row0<10) continue;
5044 if (row0>150) continue;
5047 Float_t dens00=-1,dens01=-1;
5048 Float_t dens10=-1,dens11=-1;
5050 Int_t found,foundable,ishared;
5051 track0->GetClusterStatistic(0,row0-5, found, foundable,ishared,kFALSE);
5052 if (foundable>5) dens00 = Float_t(found)/Float_t(foundable);
5053 track0->GetClusterStatistic(row0+5,155, found, foundable,ishared,kFALSE);
5054 if (foundable>5) dens01 = Float_t(found)/Float_t(foundable);
5056 track1->GetClusterStatistic(0,row0-5, found, foundable,ishared,kFALSE);
5057 if (foundable>10) dens10 = Float_t(found)/Float_t(foundable);
5058 track1->GetClusterStatistic(row0+5,155, found, foundable,ishared,kFALSE);
5059 if (foundable>10) dens11 = Float_t(found)/Float_t(foundable);
5061 if (dens00<dens10 && dens01<dens11) continue;
5062 if (dens00>dens10 && dens01>dens11) continue;
5063 if (TMath::Max(dens00,dens10)<0.1) continue;
5064 if (TMath::Max(dens01,dens11)<0.3) continue;
5066 if (TMath::Min(dens00,dens10)>0.6) continue;
5067 if (TMath::Min(dens01,dens11)>0.6) continue;
5070 AliTPCseed * ktrack0, *ktrack1;
5079 if (TMath::Abs(ktrack0->GetC())>5) continue; // cut on the curvature for mother particle
5080 AliExternalTrackParam paramm(*ktrack0);
5081 AliExternalTrackParam paramd(*ktrack1);
5082 if (row0>60&&ktrack1->GetReference().GetX()>90.)new (¶md) AliExternalTrackParam(ktrack1->GetReference());
5085 kink->SetMother(paramm);
5086 kink->SetDaughter(paramd);
5089 Float_t x[3] = { kink->GetPosition()[0],kink->GetPosition()[1],kink->GetPosition()[2]};
5091 fkParam->Transform0to1(x,index);
5092 fkParam->Transform1to2(x,index);
5093 row0 = GetRowNumber(x[0]);
5095 if (kink->GetR()<100) continue;
5096 if (kink->GetR()>240) continue;
5097 if (kink->GetPosition()[2]/kink->GetR()>AliTPCReconstructor::GetCtgRange()) continue; //out of fiducial volume
5098 if (kink->GetDistance()>cdist3) continue;
5099 Float_t dird = kink->GetDaughterP()[0]*kink->GetPosition()[0]+kink->GetDaughterP()[1]*kink->GetPosition()[1]; // rough direction estimate
5100 if (dird<0) continue;
5102 Float_t dirm = kink->GetMotherP()[0]*kink->GetPosition()[0]+kink->GetMotherP()[1]*kink->GetPosition()[1]; // rough direction estimate
5103 if (dirm<0) continue;
5104 Float_t mpt = TMath::Sqrt(kink->GetMotherP()[0]*kink->GetMotherP()[0]+kink->GetMotherP()[1]*kink->GetMotherP()[1]);
5105 if (mpt<0.2) continue;
5108 //for high momenta momentum not defined well in first iteration
5109 Double_t qt = TMath::Sin(kink->GetAngle(2))*ktrack1->GetP();
5110 if (qt>0.35) continue;
5113 kink->SetLabel(CookLabel(ktrack0,0.4,0,row0),0);
5114 kink->SetLabel(CookLabel(ktrack1,0.4,row0,160),1);
5116 kink->SetTPCDensity(dens00,0,0);
5117 kink->SetTPCDensity(dens01,0,1);
5118 kink->SetTPCDensity(dens10,1,0);
5119 kink->SetTPCDensity(dens11,1,1);
5120 kink->SetIndex(i,0);
5121 kink->SetIndex(j,1);
5124 kink->SetTPCDensity(dens10,0,0);
5125 kink->SetTPCDensity(dens11,0,1);
5126 kink->SetTPCDensity(dens00,1,0);
5127 kink->SetTPCDensity(dens01,1,1);
5128 kink->SetIndex(j,0);
5129 kink->SetIndex(i,1);
5132 if (mpt<1||kink->GetAngle(2)>0.1){
5133 // angle and densities not defined yet
5134 if (kink->GetTPCDensityFactor()<0.8) continue;
5135 if ((2-kink->GetTPCDensityFactor())*kink->GetDistance() >0.25) continue;
5136 if (kink->GetAngle(2)*ktrack0->GetP()<0.003) continue; //too small angle
5137 if (kink->GetAngle(2)>0.2&&kink->GetTPCDensityFactor()<1.15) continue;
5138 if (kink->GetAngle(2)>0.2&&kink->GetTPCDensity(0,1)>0.05) continue;
5140 Float_t criticalangle = track0->GetSigmaSnp2()+track0->GetSigmaTgl2();
5141 criticalangle+= track1->GetSigmaSnp2()+track1->GetSigmaTgl2();
5142 criticalangle= 3*TMath::Sqrt(criticalangle);
5143 if (criticalangle>0.02) criticalangle=0.02;
5144 if (kink->GetAngle(2)<criticalangle) continue;
5147 Int_t drow = Int_t(2.+0.5/(0.05+kink->GetAngle(2))); // overlap region defined
5148 Float_t shapesum =0;
5150 for ( Int_t row = row0-drow; row<row0+drow;row++){
5151 if (row<0) continue;
5152 if (row>155) continue;
5153 if (ktrack0->GetClusterPointer(row)){
5154 AliTPCTrackerPoint *point =ktrack0->GetTrackPoint(row);
5155 shapesum+=point->GetSigmaY()+point->GetSigmaZ();
5158 if (ktrack1->GetClusterPointer(row)){
5159 AliTPCTrackerPoint *point =ktrack1->GetTrackPoint(row);
5160 shapesum+=point->GetSigmaY()+point->GetSigmaZ();
5165 kink->SetShapeFactor(-1.);
5168 kink->SetShapeFactor(shapesum/sum);
5170 // esd->AddKink(kink);
5172 // kink->SetMother(paramm);
5173 //kink->SetDaughter(paramd);
5175 Double_t chi2P2 = paramm.GetParameter()[2]-paramd.GetParameter()[2];
5177 chi2P2/=paramm.GetCovariance()[5]+paramd.GetCovariance()[5];
5178 Double_t chi2P3 = paramm.GetParameter()[3]-paramd.GetParameter()[3];
5180 chi2P3/=paramm.GetCovariance()[9]+paramd.GetCovariance()[9];
5182 if (AliTPCReconstructor::StreamLevel()>1) {
5183 (*fDebugStreamer)<<"kinkLpt"<<
5191 if ( chi2P2+chi2P3<AliTPCReconstructor::GetRecoParam()->GetKinkAngleCutChi2(0)){
5195 kinks->AddLast(kink);
5201 // sort the kinks according quality - and refit them towards vertex
5203 Int_t nkinks = kinks->GetEntriesFast();
5204 Float_t *quality = new Float_t[nkinks];
5205 Int_t *indexes = new Int_t[nkinks];
5206 AliTPCseed *mothers = new AliTPCseed[nkinks];
5207 AliTPCseed *daughters = new AliTPCseed[nkinks];
5210 for (Int_t i=0;i<nkinks;i++){
5212 AliKink *kinkl = (AliKink*)kinks->At(i);
5214 // refit kinks towards vertex
5216 Int_t index0 = kinkl->GetIndex(0);
5217 Int_t index1 = kinkl->GetIndex(1);
5218 AliTPCseed * ktrack0 = (AliTPCseed*)array->At(index0);
5219 AliTPCseed * ktrack1 = (AliTPCseed*)array->At(index1);
5221 Int_t sumn=ktrack0->GetNumberOfClusters()+ktrack1->GetNumberOfClusters();
5223 // Refit Kink under if too small angle
5225 if (kinkl->GetAngle(2)<0.05){
5226 kinkl->SetTPCRow0(GetRowNumber(kinkl->GetR()));
5227 Int_t row0 = kinkl->GetTPCRow0();
5228 Int_t drow = Int_t(2.+0.5/(0.05+kinkl->GetAngle(2)));
5231 Int_t last = row0-drow;
5232 if (last<40) last=40;
5233 if (last<ktrack0->GetFirstPoint()+25) last = ktrack0->GetFirstPoint()+25;
5234 AliTPCseed* seed0 = ReSeed(ktrack0,last,kFALSE);
5237 Int_t first = row0+drow;
5238 if (first>130) first=130;
5239 if (first>ktrack1->GetLastPoint()-25) first = TMath::Max(ktrack1->GetLastPoint()-25,30);
5240 AliTPCseed* seed1 = ReSeed(ktrack1,first,kTRUE);
5242 if (seed0 && seed1){
5243 kinkl->SetStatus(1,8);
5244 if (RefitKink(*seed0,*seed1,*kinkl)) kinkl->SetStatus(1,9);
5245 row0 = GetRowNumber(kinkl->GetR());
5246 sumn = seed0->GetNumberOfClusters()+seed1->GetNumberOfClusters();
5247 mothers[i] = *seed0;
5248 daughters[i] = *seed1;
5251 delete kinks->RemoveAt(i);
5252 if (seed0) delete seed0;
5253 if (seed1) delete seed1;
5256 if (kinkl->GetDistance()>0.5 || kinkl->GetR()<110 || kinkl->GetR()>240) {
5257 delete kinks->RemoveAt(i);
5258 if (seed0) delete seed0;
5259 if (seed1) delete seed1;
5267 if (kinkl) quality[i] = 160*((0.1+kinkl->GetDistance())*(2.-kinkl->GetTPCDensityFactor()))/(sumn+40.); //the longest -clossest will win
5269 TMath::Sort(nkinks,quality,indexes,kFALSE);
5271 //remove double find kinks
5273 for (Int_t ikink0=1;ikink0<nkinks;ikink0++){
5274 AliKink * kink0 = (AliKink*) kinks->At(indexes[ikink0]);
5275 if (!kink0) continue;
5277 for (Int_t ikink1=0;ikink1<ikink0;ikink1++){
5278 kink0 = (AliKink*) kinks->At(indexes[ikink0]);
5279 if (!kink0) continue;
5280 AliKink * kink1 = (AliKink*) kinks->At(indexes[ikink1]);
5281 if (!kink1) continue;
5282 // if not close kink continue
5283 if (TMath::Abs(kink1->GetPosition()[2]-kink0->GetPosition()[2])>10) continue;
5284 if (TMath::Abs(kink1->GetPosition()[1]-kink0->GetPosition()[1])>10) continue;
5285 if (TMath::Abs(kink1->GetPosition()[0]-kink0->GetPosition()[0])>10) continue;
5287 AliTPCseed &mother0 = mothers[indexes[ikink0]];
5288 AliTPCseed &daughter0 = daughters[indexes[ikink0]];
5289 AliTPCseed &mother1 = mothers[indexes[ikink1]];
5290 AliTPCseed &daughter1 = daughters[indexes[ikink1]];
5291 Int_t row0 = (kink0->GetTPCRow0()+kink1->GetTPCRow0())/2;
5300 for (Int_t i=0;i<row0;i++){
5301 if (mother0.GetClusterIndex(i)>0 && mother1.GetClusterIndex(i)>0){
5304 if (mother0.GetClusterIndex(i)==mother1.GetClusterIndex(i)){
5311 for (Int_t i=row0;i<158;i++){
5312 if (daughter0.GetClusterIndex(i)>0 && daughter0.GetClusterIndex(i)>0){
5315 if (mother0.GetClusterIndex(i)==mother1.GetClusterIndex(i)){
5321 Float_t ratio = Float_t(same+1)/Float_t(both+1);
5322 Float_t ratiom = Float_t(samem+1)/Float_t(bothm+1);
5323 Float_t ratiod = Float_t(samed+1)/Float_t(bothd+1);
5324 if (ratio>0.3 && ratiom>0.5 &&ratiod>0.5) {
5325 Int_t sum0 = mother0.GetNumberOfClusters()+daughter0.GetNumberOfClusters();
5326 Int_t sum1 = mother1.GetNumberOfClusters()+daughter1.GetNumberOfClusters();
5328 shared[kink0->GetIndex(0)]= kTRUE;
5329 shared[kink0->GetIndex(1)]= kTRUE;
5330 delete kinks->RemoveAt(indexes[ikink0]);
5334 shared[kink1->GetIndex(0)]= kTRUE;
5335 shared[kink1->GetIndex(1)]= kTRUE;
5336 delete kinks->RemoveAt(indexes[ikink1]);
5343 for (Int_t i=0;i<nkinks;i++){
5344 AliKink * kinkl = (AliKink*) kinks->At(indexes[i]);
5345 if (!kinkl) continue;
5346 kinkl->SetTPCRow0(GetRowNumber(kinkl->GetR()));
5347 Int_t index0 = kinkl->GetIndex(0);
5348 Int_t index1 = kinkl->GetIndex(1);
5349 if (circular[index0]||(circular[index1]&&kinkl->GetDistance()>0.2)) continue;
5350 kinkl->SetMultiple(usage[index0],0);
5351 kinkl->SetMultiple(usage[index1],1);
5352 if (kinkl->GetMultiple()[0]+kinkl->GetMultiple()[1]>2) continue;
5353 if (kinkl->GetMultiple()[0]+kinkl->GetMultiple()[1]>0 && quality[indexes[i]]>0.2) continue;
5354 if (kinkl->GetMultiple()[0]+kinkl->GetMultiple()[1]>0 && kinkl->GetDistance()>0.2) continue;
5355 if (circular[index0]||(circular[index1]&&kinkl->GetDistance()>0.1)) continue;
5357 AliTPCseed * ktrack0 = (AliTPCseed*)array->At(index0);
5358 AliTPCseed * ktrack1 = (AliTPCseed*)array->At(index1);
5359 if (!ktrack0 || !ktrack1) continue;
5360 Int_t index = esd->AddKink(kinkl);
5363 if ( ktrack0->GetKinkIndex(0)==0 && ktrack1->GetKinkIndex(0)==0) { //best kink
5364 if (mothers[indexes[i]].GetNumberOfClusters()>20 && daughters[indexes[i]].GetNumberOfClusters()>20 && (mothers[indexes[i]].GetNumberOfClusters()+daughters[indexes[i]].GetNumberOfClusters())>100){
5365 *ktrack0 = mothers[indexes[i]];
5366 *ktrack1 = daughters[indexes[i]];
5370 ktrack0->SetKinkIndex(usage[index0],-(index+1));
5371 ktrack1->SetKinkIndex(usage[index1], (index+1));
5376 // Remove tracks corresponding to shared kink's
5378 for (Int_t i=0;i<nentries;i++){
5379 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5380 if (!track0) continue;
5381 if (track0->GetKinkIndex(0)!=0) continue;
5382 if (shared[i]) delete array->RemoveAt(i);
5387 RemoveUsed2(array,0.5,0.4,30);
5389 for (Int_t i=0;i<nentries;i++){
5390 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5391 if (!track0) continue;
5392 track0->CookdEdx(0.02,0.6);
5396 for (Int_t i=0;i<nentries;i++){
5397 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5398 if (!track0) continue;
5399 if (track0->Pt()<1.4) continue;
5400 //remove double high momenta tracks - overlapped with kink candidates
5403 for (Int_t icl=track0->GetFirstPoint();icl<track0->GetLastPoint(); icl++){
5404 if (track0->GetClusterPointer(icl)!=0){
5406 if (track0->GetClusterPointer(icl)->IsUsed(10)) ishared++;
5409 if (Float_t(ishared+1)/Float_t(all+1)>0.5) {
5410 delete array->RemoveAt(i);
5414 if (track0->GetKinkIndex(0)!=0) continue;
5415 if (track0->GetNumberOfClusters()<80) continue;
5417 AliTPCseed *pmother = new AliTPCseed();
5418 AliTPCseed *pdaughter = new AliTPCseed();
5419 AliKink *pkink = new AliKink;
5421 AliTPCseed & mother = *pmother;
5422 AliTPCseed & daughter = *pdaughter;
5423 AliKink & kinkl = *pkink;
5424 if (CheckKinkPoint(track0,mother,daughter, kinkl)){
5425 if (mother.GetNumberOfClusters()<30||daughter.GetNumberOfClusters()<20) {
5429 continue; //too short tracks
5431 if (mother.Pt()<1.4) {
5437 Int_t row0= kinkl.GetTPCRow0();
5438 if (kinkl.GetDistance()>0.5 || kinkl.GetR()<110. || kinkl.GetR()>240.) {
5445 Int_t index = esd->AddKink(&kinkl);
5446 mother.SetKinkIndex(0,-(index+1));
5447 daughter.SetKinkIndex(0,index+1);
5448 if (mother.GetNumberOfClusters()>50) {
5449 delete array->RemoveAt(i);
5450 array->AddAt(new AliTPCseed(mother),i);
5453 array->AddLast(new AliTPCseed(mother));
5455 array->AddLast(new AliTPCseed(daughter));
5456 for (Int_t icl=0;icl<row0;icl++) {
5457 if (mother.GetClusterPointer(icl)) mother.GetClusterPointer(icl)->Use(20);
5460 for (Int_t icl=row0;icl<158;icl++) {
5461 if (daughter.GetClusterPointer(icl)) daughter.GetClusterPointer(icl)->Use(20);
5470 delete [] daughters;
5492 AliInfo(Form("Ncandidates=\t%d\t%d\t%d\t%d\n",esd->GetNumberOfKinks(),ncandidates,ntracks,nall));
5497 Int_t AliTPCtrackerMI::RefitKink(AliTPCseed &mother, AliTPCseed &daughter, const AliESDkink &knk)
5500 // refit kink towards to the vertex
5503 AliKink &kink=(AliKink &)knk;
5505 Int_t row0 = GetRowNumber(kink.GetR());
5506 FollowProlongation(mother,0);
5507 mother.Reset(kFALSE);
5509 FollowProlongation(daughter,row0);
5510 daughter.Reset(kFALSE);
5511 FollowBackProlongation(daughter,158);
5512 daughter.Reset(kFALSE);
5513 Int_t first = TMath::Max(row0-20,30);
5514 Int_t last = TMath::Min(row0+20,140);
5516 const Int_t kNdiv =5;
5517 AliTPCseed param0[kNdiv]; // parameters along the track
5518 AliTPCseed param1[kNdiv]; // parameters along the track
5519 AliKink kinks[kNdiv]; // corresponding kink parameters
5522 for (Int_t irow=0; irow<kNdiv;irow++){
5523 rows[irow] = first +((last-first)*irow)/(kNdiv-1);
5525 // store parameters along the track
5527 for (Int_t irow=0;irow<kNdiv;irow++){
5528 FollowBackProlongation(mother, rows[irow]);
5529 FollowProlongation(daughter,rows[kNdiv-1-irow]);
5530 param0[irow] = mother;
5531 param1[kNdiv-1-irow] = daughter;
5535 for (Int_t irow=0; irow<kNdiv-1;irow++){
5536 if (param0[irow].GetNumberOfClusters()<kNdiv||param1[irow].GetNumberOfClusters()<kNdiv) continue;
5537 kinks[irow].SetMother(param0[irow]);
5538 kinks[irow].SetDaughter(param1[irow]);
5539 kinks[irow].Update();
5542 // choose kink with best "quality"
5544 Double_t mindist = 10000;
5545 for (Int_t irow=0;irow<kNdiv;irow++){
5546 if (param0[irow].GetNumberOfClusters()<20||param1[irow].GetNumberOfClusters()<20) continue;
5547 if (TMath::Abs(kinks[irow].GetR())>240.) continue;
5548 if (TMath::Abs(kinks[irow].GetR())<100.) continue;
5550 Float_t normdist = TMath::Abs(param0[irow].GetX()-kinks[irow].GetR())*(0.1+kink.GetDistance());
5551 normdist/= (param0[irow].GetNumberOfClusters()+param1[irow].GetNumberOfClusters()+40.);
5552 if (normdist < mindist){
5558 if (index==-1) return 0;
5561 param0[index].Reset(kTRUE);
5562 FollowProlongation(param0[index],0);
5564 mother = param0[index];
5565 daughter = param1[index]; // daughter in vertex
5567 kink.SetMother(mother);
5568 kink.SetDaughter(daughter);
5570 kink.SetTPCRow0(GetRowNumber(kink.GetR()));
5571 kink.SetTPCncls(param0[index].GetNumberOfClusters(),0);
5572 kink.SetTPCncls(param1[index].GetNumberOfClusters(),1);
5573 kink.SetLabel(CookLabel(&mother,0.4, 0,kink.GetTPCRow0()),0);
5574 kink.SetLabel(CookLabel(&daughter,0.4, kink.GetTPCRow0(),160),1);
5575 mother.SetLabel(kink.GetLabel(0));
5576 daughter.SetLabel(kink.GetLabel(1));
5582 void AliTPCtrackerMI::UpdateKinkQualityM(AliTPCseed * seed){
5584 // update Kink quality information for mother after back propagation
5586 if (seed->GetKinkIndex(0)>=0) return;
5587 for (Int_t ikink=0;ikink<3;ikink++){
5588 Int_t index = seed->GetKinkIndex(ikink);
5589 if (index>=0) break;
5590 index = TMath::Abs(index)-1;
5591 AliESDkink * kink = fEvent->GetKink(index);
5592 kink->SetTPCDensity(-1,0,0);
5593 kink->SetTPCDensity(1,0,1);
5595 Int_t row0 = kink->GetTPCRow0() - 2 - Int_t( 0.5/ (0.05+kink->GetAngle(2)));
5596 if (row0<15) row0=15;
5598 Int_t row1 = kink->GetTPCRow0() + 2 + Int_t( 0.5/ (0.05+kink->GetAngle(2)));
5599 if (row1>145) row1=145;
5601 Int_t found,foundable,shared;
5602 seed->GetClusterStatistic(0,row0, found, foundable,shared,kFALSE);
5603 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),0,0);
5604 seed->GetClusterStatistic(row1,155, found, foundable,shared,kFALSE);
5605 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),0,1);
5610 void AliTPCtrackerMI::UpdateKinkQualityD(AliTPCseed * seed){
5612 // update Kink quality information for daughter after refit
5614 if (seed->GetKinkIndex(0)<=0) return;
5615 for (Int_t ikink=0;ikink<3;ikink++){
5616 Int_t index = seed->GetKinkIndex(ikink);
5617 if (index<=0) break;
5618 index = TMath::Abs(index)-1;
5619 AliESDkink * kink = fEvent->GetKink(index);
5620 kink->SetTPCDensity(-1,1,0);
5621 kink->SetTPCDensity(-1,1,1);
5623 Int_t row0 = kink->GetTPCRow0() -2 - Int_t( 0.5/ (0.05+kink->GetAngle(2)));
5624 if (row0<15) row0=15;
5626 Int_t row1 = kink->GetTPCRow0() +2 + Int_t( 0.5/ (0.05+kink->GetAngle(2)));
5627 if (row1>145) row1=145;
5629 Int_t found,foundable,shared;
5630 seed->GetClusterStatistic(0,row0, found, foundable,shared,kFALSE);
5631 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),1,0);
5632 seed->GetClusterStatistic(row1,155, found, foundable,shared,kFALSE);
5633 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),1,1);
5639 Int_t AliTPCtrackerMI::CheckKinkPoint(AliTPCseed*seed,AliTPCseed &mother, AliTPCseed &daughter, const AliESDkink &knk)
5642 // check kink point for given track
5643 // if return value=0 kink point not found
5644 // otherwise seed0 correspond to mother particle
5645 // seed1 correspond to daughter particle
5646 // kink parameter of kink point
5647 AliKink &kink=(AliKink &)knk;
5649 Int_t middlerow = (seed->GetFirstPoint()+seed->GetLastPoint())/2;
5650 Int_t first = seed->GetFirstPoint();
5651 Int_t last = seed->GetLastPoint();
5652 if (last-first<20) return 0; // shortest length - 2*30 = 60 pad-rows
5655 AliTPCseed *seed1 = ReSeed(seed,middlerow+20, kTRUE); //middle of chamber
5656 if (!seed1) return 0;
5657 FollowProlongation(*seed1,seed->GetLastPoint()-20);
5658 seed1->Reset(kTRUE);
5659 FollowProlongation(*seed1,158);
5660 seed1->Reset(kTRUE);
5661 last = seed1->GetLastPoint();
5663 AliTPCseed *seed0 = new AliTPCseed(*seed);
5664 seed0->Reset(kFALSE);
5667 AliTPCseed param0[20]; // parameters along the track
5668 AliTPCseed param1[20]; // parameters along the track
5669 AliKink kinks[20]; // corresponding kink parameters
5671 for (Int_t irow=0; irow<20;irow++){
5672 rows[irow] = first +((last-first)*irow)/19;
5674 // store parameters along the track
5676 for (Int_t irow=0;irow<20;irow++){
5677 FollowBackProlongation(*seed0, rows[irow]);
5678 FollowProlongation(*seed1,rows[19-irow]);
5679 param0[irow] = *seed0;
5680 param1[19-irow] = *seed1;
5684 for (Int_t irow=0; irow<19;irow++){
5685 kinks[irow].SetMother(param0[irow]);
5686 kinks[irow].SetDaughter(param1[irow]);
5687 kinks[irow].Update();
5690 // choose kink with biggest change of angle
5692 Double_t maxchange= 0;
5693 for (Int_t irow=1;irow<19;irow++){
5694 if (TMath::Abs(kinks[irow].GetR())>240.) continue;
5695 if (TMath::Abs(kinks[irow].GetR())<110.) continue;
5696 Float_t quality = TMath::Abs(kinks[irow].GetAngle(2))/(3.+TMath::Abs(kinks[irow].GetR()-param0[irow].GetX()));
5697 if ( quality > maxchange){
5698 maxchange = quality;
5705 if (index<0) return 0;
5707 Int_t row0 = GetRowNumber(kinks[index].GetR()); //row 0 estimate
5708 seed0 = new AliTPCseed(param0[index]);
5709 seed1 = new AliTPCseed(param1[index]);
5710 seed0->Reset(kFALSE);
5711 seed1->Reset(kFALSE);
5712 seed0->ResetCovariance(10.);
5713 seed1->ResetCovariance(10.);
5714 FollowProlongation(*seed0,0);
5715 FollowBackProlongation(*seed1,158);
5716 mother = *seed0; // backup mother at position 0
5717 seed0->Reset(kFALSE);
5718 seed1->Reset(kFALSE);
5719 seed0->ResetCovariance(10.);
5720 seed1->ResetCovariance(10.);
5722 first = TMath::Max(row0-20,0);
5723 last = TMath::Min(row0+20,158);
5725 for (Int_t irow=0; irow<20;irow++){
5726 rows[irow] = first +((last-first)*irow)/19;
5728 // store parameters along the track
5730 for (Int_t irow=0;irow<20;irow++){
5731 FollowBackProlongation(*seed0, rows[irow]);
5732 FollowProlongation(*seed1,rows[19-irow]);
5733 param0[irow] = *seed0;
5734 param1[19-irow] = *seed1;
5738 for (Int_t irow=0; irow<19;irow++){
5739 kinks[irow].SetMother(param0[irow]);
5740 kinks[irow].SetDaughter(param1[irow]);
5741 // param0[irow].Dump();
5742 //param1[irow].Dump();
5743 kinks[irow].Update();
5746 // choose kink with biggest change of angle
5749 for (Int_t irow=0;irow<20;irow++){
5750 if (TMath::Abs(kinks[irow].GetR())>250.) continue;
5751 if (TMath::Abs(kinks[irow].GetR())<90.) continue;
5752 Float_t quality = TMath::Abs(kinks[irow].GetAngle(2))/(3.+TMath::Abs(kinks[irow].GetR()-param0[irow].GetX()));
5753 if ( quality > maxchange){
5754 maxchange = quality;
5761 if (index==-1 || param0[index].GetNumberOfClusters()+param1[index].GetNumberOfClusters()<100){
5767 // Float_t anglesigma = TMath::Sqrt(param0[index].fC22+param0[index].fC33+param1[index].fC22+param1[index].fC33);
5769 kink.SetMother(param0[index]);
5770 kink.SetDaughter(param1[index]);
5773 Double_t chi2P2 = param0[index].GetParameter()[2]-param1[index].GetParameter()[2];
5775 chi2P2/=param0[index].GetCovariance()[5]+param1[index].GetCovariance()[5];
5776 Double_t chi2P3 = param0[index].GetParameter()[3]-param1[index].GetParameter()[3];
5778 chi2P3/=param0[index].GetCovariance()[9]+param1[index].GetCovariance()[9];
5780 if (AliTPCReconstructor::StreamLevel()>1) {
5781 (*fDebugStreamer)<<"kinkHpt"<<
5784 "p0.="<<¶m0[index]<<
5785 "p1.="<<¶m1[index]<<
5789 if ( chi2P2+chi2P3<AliTPCReconstructor::GetRecoParam()->GetKinkAngleCutChi2(0)){
5796 row0 = GetRowNumber(kink.GetR());
5797 kink.SetTPCRow0(row0);
5798 kink.SetLabel(CookLabel(seed0,0.5,0,row0),0);
5799 kink.SetLabel(CookLabel(seed1,0.5,row0,158),1);
5800 kink.SetIndex(-10,0);
5801 kink.SetIndex(int(param0[index].GetNumberOfClusters()+param1[index].GetNumberOfClusters()),1);
5802 kink.SetTPCncls(param0[index].GetNumberOfClusters(),0);
5803 kink.SetTPCncls(param1[index].GetNumberOfClusters(),1);
5806 // new (&mother) AliTPCseed(param0[index]);
5807 daughter = param1[index];
5808 daughter.SetLabel(kink.GetLabel(1));
5809 param0[index].Reset(kTRUE);
5810 FollowProlongation(param0[index],0);
5811 mother = param0[index];
5812 mother.SetLabel(kink.GetLabel(0));
5813 if ( chi2P2+chi2P3<AliTPCReconstructor::GetRecoParam()->GetKinkAngleCutChi2(1)){
5825 AliTPCseed* AliTPCtrackerMI::ReSeed(AliTPCseed *t)
5828 // reseed - refit - track
5831 // Int_t last = fSectors->GetNRows()-1;
5833 if (fSectors == fOuterSec){
5834 first = TMath::Max(first, t->GetFirstPoint()-fInnerSec->GetNRows());
5838 first = t->GetFirstPoint();
5840 AliTPCseed * seed = MakeSeed(t,0.1,0.5,0.9);
5841 FollowBackProlongation(*t,fSectors->GetNRows()-1);
5843 FollowProlongation(*t,first);
5853 //_____________________________________________________________________________
5854 Int_t AliTPCtrackerMI::ReadSeeds(const TFile *inp) {
5855 //-----------------------------------------------------------------
5856 // This function reades track seeds.
5857 //-----------------------------------------------------------------
5858 TDirectory *savedir=gDirectory;
5860 TFile *in=(TFile*)inp;
5861 if (!in->IsOpen()) {
5862 cerr<<"AliTPCtrackerMI::ReadSeeds(): input file is not open !\n";
5867 TTree *seedTree=(TTree*)in->Get("Seeds");
5869 cerr<<"AliTPCtrackerMI::ReadSeeds(): ";
5870 cerr<<"can't get a tree with track seeds !\n";
5873 AliTPCtrack *seed=new AliTPCtrack;
5874 seedTree->SetBranchAddress("tracks",&seed);
5876 if (fSeeds==0) fSeeds=new TObjArray(15000);
5878 Int_t n=(Int_t)seedTree->GetEntries();
5879 for (Int_t i=0; i<n; i++) {
5880 seedTree->GetEvent(i);
5881 fSeeds->AddLast(new AliTPCseed(*seed/*,seed->GetAlpha()*/));
5890 Int_t AliTPCtrackerMI::Clusters2Tracks (AliESDEvent *const esd)
5893 // clusters to tracks
5895 if (fSeeds) DeleteSeeds();
5898 AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ;
5899 transform->SetCurrentTimeStamp( esd->GetTimeStamp());
5900 transform->SetCurrentRun(esd->GetRunNumber());
5903 if (!fSeeds) return 1;
5905 if (AliTPCReconstructor::StreamLevel()>3) DumpClusters(0,fSeeds);
5911 //_____________________________________________________________________________
5912 Int_t AliTPCtrackerMI::Clusters2Tracks() {
5913 //-----------------------------------------------------------------
5914 // This is a track finder.
5915 //-----------------------------------------------------------------
5916 TDirectory *savedir=gDirectory;
5920 fSeeds = Tracking();
5923 Info("Clusters2Tracks","Time for tracking: \t");timer.Print();timer.Start();
5925 //activate again some tracks
5926 for (Int_t i=0; i<fSeeds->GetEntriesFast(); i++) {
5927 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
5929 Int_t nc=t.GetNumberOfClusters();
5931 delete fSeeds->RemoveAt(i);
5935 if (pt->GetRemoval()==10) {
5936 if (pt->GetDensityFirst(20)>0.8 || pt->GetDensityFirst(30)>0.8 || pt->GetDensityFirst(40)>0.7)
5937 pt->Desactivate(10); // make track again active
5939 pt->Desactivate(20);
5940 delete fSeeds->RemoveAt(i);
5945 RemoveUsed2(fSeeds,0.85,0.85,0);
5946 if (AliTPCReconstructor::GetRecoParam()->GetDoKinks()) FindKinks(fSeeds,fEvent);
5947 //FindCurling(fSeeds, fEvent,0);
5948 if (AliTPCReconstructor::StreamLevel()>5) FindMultiMC(fSeeds, fEvent,-1); // find multi found tracks
5949 RemoveUsed2(fSeeds,0.5,0.4,20);
5950 FindSplitted(fSeeds, fEvent,0); // find multi found tracks
5951 if (AliTPCReconstructor::StreamLevel()>5) FindMultiMC(fSeeds, fEvent,0); // find multi found tracks
5954 // // refit short tracks
5956 Int_t nseed=fSeeds->GetEntriesFast();
5959 for (Int_t i=0; i<nseed; i++) {
5960 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
5962 Int_t nc=t.GetNumberOfClusters();
5964 delete fSeeds->RemoveAt(i);
5967 CookLabel(pt,0.1); //For comparison only
5968 //if ((pt->IsActive() || (pt->fRemoval==10) )&& nc>50 &&pt->GetNumberOfClusters()>0.4*pt->fNFoundable){
5969 if ((pt->IsActive() || (pt->GetRemoval()==10) )){
5971 if (fDebug>0) cerr<<found<<'\r';
5975 delete fSeeds->RemoveAt(i);
5979 //RemoveOverlap(fSeeds,0.99,7,kTRUE);
5981 //RemoveUsed(fSeeds,0.9,0.9,6);
5983 nseed=fSeeds->GetEntriesFast();
5985 for (Int_t i=0; i<nseed; i++) {
5986 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
5988 Int_t nc=t.GetNumberOfClusters();
5990 delete fSeeds->RemoveAt(i);
5994 t.CookdEdx(0.02,0.6);
5995 // CheckKinkPoint(&t,0.05);
5996 //if ((pt->IsActive() || (pt->fRemoval==10) )&& nc>50 &&pt->GetNumberOfClusters()>0.4*pt->fNFoundable){
5997 if ((pt->IsActive() || (pt->GetRemoval()==10) )){
6005 delete fSeeds->RemoveAt(i);
6006 //AliTPCseed * seed1 = ReSeed(pt,0.05,0.5,1);
6008 // FollowProlongation(*seed1,0);
6009 // Int_t n = seed1->GetNumberOfClusters();
6010 // printf("fP4\t%f\t%f\n",seed1->GetC(),pt->GetC());
6011 // printf("fN\t%d\t%d\n", seed1->GetNumberOfClusters(),pt->GetNumberOfClusters());
6014 //AliTPCseed * seed2 = ReSeed(pt,0.95,0.5,0.05);
6018 SortTracks(fSeeds, 1);
6022 PrepareForBackProlongation(fSeeds,5.);
6023 PropagateBack(fSeeds);
6024 printf("Time for back propagation: \t");timer.Print();timer.Start();
6028 PrepareForProlongation(fSeeds,5.);
6029 PropagateForard2(fSeeds);
6031 printf("Time for FORWARD propagation: \t");timer.Print();timer.Start();
6032 // RemoveUsed(fSeeds,0.7,0.7,6);
6033 //RemoveOverlap(fSeeds,0.9,7,kTRUE);
6035 nseed=fSeeds->GetEntriesFast();
6037 for (Int_t i=0; i<nseed; i++) {
6038 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
6040 Int_t nc=t.GetNumberOfClusters();
6042 delete fSeeds->RemoveAt(i);
6045 t.CookdEdx(0.02,0.6);
6046 // CookLabel(pt,0.1); //For comparison only
6047 //if ((pt->IsActive() || (pt->fRemoval==10) )&& nc>50 &&pt->GetNumberOfClusters()>0.4*pt->fNFoundable){
6048 if ((pt->IsActive() || (pt->fRemoval==10) )){
6049 cerr<<found++<<'\r';
6052 delete fSeeds->RemoveAt(i);
6057 // fNTracks = found;
6059 Info("Clusters2Tracks","Time for overlap removal, track writing and dedx cooking: \t"); timer.Print();timer.Start();
6062 // cerr<<"Number of found tracks : "<<"\t"<<found<<endl;
6063 Info("Clusters2Tracks","Number of found tracks %d",found);
6065 // UnloadClusters();
6070 void AliTPCtrackerMI::Tracking(TObjArray * arr)
6073 // tracking of the seeds
6076 fSectors = fOuterSec;
6077 ParallelTracking(arr,150,63);
6078 fSectors = fOuterSec;
6079 ParallelTracking(arr,63,0);
6082 TObjArray * AliTPCtrackerMI::Tracking(Int_t seedtype, Int_t i1, Int_t i2, Float_t cuts[4], Float_t dy, Int_t dsec)
6087 TObjArray * arr = new TObjArray;
6089 fSectors = fOuterSec;
6092 for (Int_t sec=0;sec<fkNOS;sec++){
6093 if (seedtype==3) MakeSeeds3(arr,sec,i1,i2,cuts,dy, dsec);
6094 if (seedtype==4) MakeSeeds5(arr,sec,i1,i2,cuts,dy);
6095 if (seedtype==2) MakeSeeds2(arr,sec,i1,i2,cuts,dy);
6098 Info("Tracking","\nSeeding - %d\t%d\t%d\t%d\n",seedtype,i1,i2,arr->GetEntriesFast());
6110 TObjArray * AliTPCtrackerMI::Tracking()
6114 if (AliTPCReconstructor::GetRecoParam()->GetSpecialSeeding()) return TrackingSpecial();
6117 Int_t nup=fOuterSec->GetNRows()+fInnerSec->GetNRows();
6119 TObjArray * seeds = new TObjArray;
6121 Int_t fLastSeedRowSec=AliTPCReconstructor::GetRecoParam()->GetLastSeedRowSec();
6122 Int_t gapPrim = AliTPCReconstructor::GetRecoParam()->GetSeedGapPrim();
6123 Int_t gapSec = AliTPCReconstructor::GetRecoParam()->GetSeedGapSec();
6131 Float_t fnumber = 3.0;
6132 Float_t fdensity = 3.0;
6137 for (Int_t delta = 0; delta<18; delta+=gapPrim){
6141 arr = Tracking(3,nup-1-delta,nup-1-delta-gap,cuts,-1,1);
6142 SumTracks(seeds,arr);
6143 SignClusters(seeds,fnumber,fdensity);
6145 for (Int_t i=2;i<6;i+=2){
6146 // seed high pt tracks
6149 arr = Tracking(3,nup-i-delta,nup-i-delta-gap,cuts,-1,0);
6150 SumTracks(seeds,arr);
6151 SignClusters(seeds,fnumber,fdensity);
6156 // RemoveUsed(seeds,0.9,0.9,1);
6157 // UnsignClusters();
6158 // SignClusters(seeds,fnumber,fdensity);
6162 for (Int_t delta = 20; delta<120; delta+=gapPrim){
6164 // seed high pt tracks
6168 arr = Tracking(3,nup-delta,nup-delta-gap,cuts,-1);
6169 SumTracks(seeds,arr);
6170 SignClusters(seeds,fnumber,fdensity);
6175 arr = Tracking(3,nup-delta-5,nup-delta-5-gap,cuts,-1);
6176 SumTracks(seeds,arr);
6177 SignClusters(seeds,fnumber,fdensity);
6188 Info("Tracking()","\n\nPrimary seeding\t%d\n\n",seeds->GetEntriesFast());
6192 // RemoveUsed(seeds,0.75,0.75,1);
6194 //SignClusters(seeds,fnumber,fdensity);
6203 arr = Tracking(4,nup-1,nup-1-gap,cuts,-1);
6204 SumTracks(seeds,arr);
6205 SignClusters(seeds,fnumber,fdensity);
6207 arr = Tracking(4,nup-2,nup-2-gap,cuts,-1);
6208 SumTracks(seeds,arr);
6209 SignClusters(seeds,fnumber,fdensity);
6211 arr = Tracking(4,nup-3,nup-3-gap,cuts,-1);
6212 SumTracks(seeds,arr);
6213 SignClusters(seeds,fnumber,fdensity);
6215 arr = Tracking(4,nup-5,nup-5-gap,cuts,-1);
6216 SumTracks(seeds,arr);
6217 SignClusters(seeds,fnumber,fdensity);
6219 arr = Tracking(4,nup-7,nup-7-gap,cuts,-1);
6220 SumTracks(seeds,arr);
6221 SignClusters(seeds,fnumber,fdensity);
6224 arr = Tracking(4,nup-9,nup-9-gap,cuts,-1);
6225 SumTracks(seeds,arr);
6226 SignClusters(seeds,fnumber,fdensity);
6230 for (Int_t delta = 9; delta<30; delta+=gapSec){
6236 arr = Tracking(4,nup-1-delta,nup-1-delta-gap,cuts,-1);
6237 SumTracks(seeds,arr);
6238 SignClusters(seeds,fnumber,fdensity);
6240 arr = Tracking(4,nup-3-delta,nup-5-delta-gap,cuts,4);
6241 SumTracks(seeds,arr);
6242 SignClusters(seeds,fnumber,fdensity);
6255 for (Int_t delta = 30; delta<fLastSeedRowSec; delta+=gapSec){
6261 arr = Tracking(4,nup-1-delta,nup-1-delta-gap,cuts,-1);
6262 SumTracks(seeds,arr);
6263 SignClusters(seeds,fnumber,fdensity);
6265 arr = Tracking(4,nup-5-delta,nup-5-delta-gap,cuts,5 );
6266 SumTracks(seeds,arr);
6267 SignClusters(seeds,fnumber,fdensity);
6271 Info("Tracking()","\n\nSecondary seeding\t%d\n\n",seeds->GetEntriesFast());
6282 TObjArray * AliTPCtrackerMI::TrackingSpecial()
6285 // seeding adjusted for laser and cosmic tests - short tracks with big inclination angle
6286 // no primary vertex seeding tried
6290 Int_t nup=fOuterSec->GetNRows()+fInnerSec->GetNRows();
6292 TObjArray * seeds = new TObjArray;
6297 Float_t fnumber = 3.0;
6298 Float_t fdensity = 3.0;
6301 cuts[0] = AliTPCReconstructor::GetRecoParam()->GetMaxC(); // max curvature
6302 cuts[1] = 3.5; // max tan(phi) angle for seeding
6303 cuts[2] = 3.; // not used (cut on z primary vertex)
6304 cuts[3] = 3.5; // max tan(theta) angle for seeding
6306 for (Int_t delta = 0; nup-delta-gap-1>0; delta+=3){
6308 arr = Tracking(4,nup-1-delta,nup-1-delta-gap,cuts,-1);
6309 SumTracks(seeds,arr);
6310 SignClusters(seeds,fnumber,fdensity);
6314 Info("Tracking()","\n\nSecondary seeding\t%d\n\n",seeds->GetEntriesFast());
6325 void AliTPCtrackerMI::SumTracks(TObjArray *arr1,TObjArray *&arr2) const
6328 //sum tracks to common container
6329 //remove suspicious tracks
6330 Int_t nseed = arr2->GetEntriesFast();
6331 for (Int_t i=0;i<nseed;i++){
6332 AliTPCseed *pt=(AliTPCseed*)arr2->UncheckedAt(i);
6335 // remove tracks with too big curvature
6337 if (TMath::Abs(pt->GetC())>AliTPCReconstructor::GetRecoParam()->GetMaxC()){
6338 delete arr2->RemoveAt(i);
6341 // REMOVE VERY SHORT TRACKS
6342 if (pt->GetNumberOfClusters()<20){
6343 delete arr2->RemoveAt(i);
6346 // NORMAL ACTIVE TRACK
6347 if (pt->IsActive()){
6348 arr1->AddLast(arr2->RemoveAt(i));
6351 //remove not usable tracks
6352 if (pt->GetRemoval()!=10){
6353 delete arr2->RemoveAt(i);
6357 // ENABLE ONLY ENOUGH GOOD STOPPED TRACKS
6358 if (pt->GetDensityFirst(20)>0.8 || pt->GetDensityFirst(30)>0.8 || pt->GetDensityFirst(40)>0.7)
6359 arr1->AddLast(arr2->RemoveAt(i));
6361 delete arr2->RemoveAt(i);
6365 delete arr2; arr2 = 0;
6370 void AliTPCtrackerMI::ParallelTracking(TObjArray *const arr, Int_t rfirst, Int_t rlast)
6373 // try to track in parralel
6375 Int_t nseed=arr->GetEntriesFast();
6376 //prepare seeds for tracking
6377 for (Int_t i=0; i<nseed; i++) {
6378 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i), &t=*pt;
6380 if (!t.IsActive()) continue;
6381 // follow prolongation to the first layer
6382 if ( (fSectors ==fInnerSec) || (t.GetFirstPoint()-fkParam->GetNRowLow()>rfirst+1) )
6383 FollowProlongation(t, rfirst+1);
6388 for (Int_t nr=rfirst; nr>=rlast; nr--){
6389 if (nr<fInnerSec->GetNRows())
6390 fSectors = fInnerSec;
6392 fSectors = fOuterSec;
6393 // make indexes with the cluster tracks for given
6395 // find nearest cluster
6396 for (Int_t i=0; i<nseed; i++) {
6397 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i), &t=*pt;
6399 if (nr==80) pt->UpdateReference();
6400 if (!pt->IsActive()) continue;
6401 // if ( (fSectors ==fOuterSec) && (pt->fFirstPoint-fkParam->GetNRowLow())<nr) continue;
6402 if (pt->GetRelativeSector()>17) {
6405 UpdateClusters(t,nr);
6407 // prolonagate to the nearest cluster - if founded
6408 for (Int_t i=0; i<nseed; i++) {
6409 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
6411 if (!pt->IsActive()) continue;
6412 // if ((fSectors ==fOuterSec) && (pt->fFirstPoint-fkParam->GetNRowLow())<nr) continue;
6413 if (pt->GetRelativeSector()>17) {
6416 FollowToNextCluster(*pt,nr);
6421 void AliTPCtrackerMI::PrepareForBackProlongation(const TObjArray *const arr,Float_t fac) const
6425 // if we use TPC track itself we have to "update" covariance
6427 Int_t nseed= arr->GetEntriesFast();
6428 for (Int_t i=0;i<nseed;i++){
6429 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
6433 //rotate to current local system at first accepted point
6434 Int_t index = pt->GetClusterIndex2(pt->GetFirstPoint());
6435 Int_t sec = (index&0xff000000)>>24;
6437 Float_t angle1 = fInnerSec->GetAlpha()*sec+fInnerSec->GetAlphaShift();
6438 if (angle1>TMath::Pi())
6439 angle1-=2.*TMath::Pi();
6440 Float_t angle2 = pt->GetAlpha();
6442 if (TMath::Abs(angle1-angle2)>0.001){
6443 if (!pt->Rotate(angle1-angle2)) return;
6444 //angle2 = pt->GetAlpha();
6445 //pt->fRelativeSector = pt->GetAlpha()/fInnerSec->GetAlpha();
6446 //if (pt->GetAlpha()<0)
6447 // pt->fRelativeSector+=18;
6448 //sec = pt->fRelativeSector;
6457 void AliTPCtrackerMI::PrepareForProlongation(TObjArray *const arr, Float_t fac) const
6461 // if we use TPC track itself we have to "update" covariance
6463 Int_t nseed= arr->GetEntriesFast();
6464 for (Int_t i=0;i<nseed;i++){
6465 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
6468 pt->SetFirstPoint(pt->GetLastPoint());
6476 Int_t AliTPCtrackerMI::PropagateBack(const TObjArray *const arr)
6479 // make back propagation
6481 Int_t nseed= arr->GetEntriesFast();
6482 for (Int_t i=0;i<nseed;i++){
6483 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
6484 if (pt&& pt->GetKinkIndex(0)<=0) {
6485 //AliTPCseed *pt2 = new AliTPCseed(*pt);
6486 fSectors = fInnerSec;
6487 //FollowBackProlongation(*pt,fInnerSec->GetNRows()-1);
6488 //fSectors = fOuterSec;
6489 FollowBackProlongation(*pt,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1,1);
6490 //if (pt->GetNumberOfClusters()<(pt->fEsd->GetTPCclusters(0)) ){
6491 // Error("PropagateBack","Not prolonged track %d",pt->GetLabel());
6492 // FollowBackProlongation(*pt2,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1);
6495 if (pt&& pt->GetKinkIndex(0)>0) {
6496 AliESDkink * kink = fEvent->GetKink(pt->GetKinkIndex(0)-1);
6497 pt->SetFirstPoint(kink->GetTPCRow0());
6498 fSectors = fInnerSec;
6499 FollowBackProlongation(*pt,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1,1);
6507 Int_t AliTPCtrackerMI::PropagateForward2(const TObjArray *const arr)
6510 // make forward propagation
6512 Int_t nseed= arr->GetEntriesFast();
6514 for (Int_t i=0;i<nseed;i++){
6515 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
6517 FollowProlongation(*pt,0,1,1);
6526 Int_t AliTPCtrackerMI::PropagateForward()
6529 // propagate track forward
6531 Int_t nseed = fSeeds->GetEntriesFast();
6532 for (Int_t i=0;i<nseed;i++){
6533 AliTPCseed *pt = (AliTPCseed*)fSeeds->UncheckedAt(i);
6535 AliTPCseed &t = *pt;
6536 Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
6537 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
6538 if (alpha < 0. ) alpha += 2.*TMath::Pi();
6539 t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN);
6543 fSectors = fOuterSec;
6544 ParallelTracking(fSeeds,fOuterSec->GetNRows()+fInnerSec->GetNRows()-1,fInnerSec->GetNRows());
6545 fSectors = fInnerSec;
6546 ParallelTracking(fSeeds,fInnerSec->GetNRows()-1,0);
6555 Int_t AliTPCtrackerMI::PropagateBack(AliTPCseed *const pt, Int_t row0, Int_t row1)
6558 // make back propagation, in between row0 and row1
6562 fSectors = fInnerSec;
6565 if (row1<fSectors->GetNRows())
6568 r1 = fSectors->GetNRows()-1;
6570 if (row0<fSectors->GetNRows()&& r1>0 )
6571 FollowBackProlongation(*pt,r1);
6572 if (row1<=fSectors->GetNRows())
6575 r1 = row1 - fSectors->GetNRows();
6576 if (r1<=0) return 0;
6577 if (r1>=fOuterSec->GetNRows()) return 0;
6578 fSectors = fOuterSec;
6579 return FollowBackProlongation(*pt,r1);
6587 void AliTPCtrackerMI::GetShape(AliTPCseed * seed, Int_t row)
6589 // gets cluster shape
6591 AliTPCClusterParam * clparam = AliTPCcalibDB::Instance()->GetClusterParam();
6592 Float_t zdrift = TMath::Abs((fkParam->GetZLength(0)-TMath::Abs(seed->GetZ())));
6593 Int_t type = (seed->GetSector() < fkParam->GetNSector()/2) ? 0: (row>126) ? 1:2;
6594 Double_t angulary = seed->GetSnp();
6596 if (TMath::Abs(angulary)>AliTPCReconstructor::GetMaxSnpTracker()) {
6597 angulary = TMath::Sign(AliTPCReconstructor::GetMaxSnpTracker(),angulary);
6600 angulary = angulary*angulary/((1.-angulary)*(1.+angulary));
6601 Double_t angularz = seed->GetTgl()*seed->GetTgl()*(1.+angulary);
6603 Double_t sigmay = clparam->GetRMS0(0,type,zdrift,TMath::Sqrt(TMath::Abs(angulary)));
6604 Double_t sigmaz = clparam->GetRMS0(1,type,zdrift,TMath::Sqrt(TMath::Abs(angularz)));
6605 seed->SetCurrentSigmaY2(sigmay*sigmay);
6606 seed->SetCurrentSigmaZ2(sigmaz*sigmaz);
6607 // Float_t sd2 = TMath::Abs((fkParam->GetZLength(0)-TMath::Abs(seed->GetZ())))*fkParam->GetDiffL()*fkParam->GetDiffL();
6608 // // Float_t padlength = fkParam->GetPadPitchLength(seed->fSector);
6609 // Float_t padlength = GetPadPitchLength(row);
6611 // Float_t sresy = (seed->GetSector() < fkParam->GetNSector()/2) ? 0.2 :0.3;
6612 // seed->SetCurrentSigmaY2(sd2+padlength*padlength*angulary/12.+sresy*sresy);
6614 // Float_t sresz = fkParam->GetZSigma();
6615 // seed->SetCurrentSigmaZ2(sd2+padlength*padlength*angularz*angularz*(1+angulary)/12.+sresz*sresz);
6617 Float_t wy = GetSigmaY(seed);
6618 Float_t wz = GetSigmaZ(seed);
6621 if (TMath::Abs(wy/seed->fCurrentSigmaY2-1)>0.0001 || TMath::Abs(wz/seed->fCurrentSigmaZ2-1)>0.0001 ){
6622 printf("problem\n");
6629 //__________________________________________________________________________
6630 void AliTPCtrackerMI::CookLabel(AliKalmanTrack *tk, Float_t wrong) const {
6631 //--------------------------------------------------------------------
6632 //This function "cooks" a track label. If label<0, this track is fake.
6633 //--------------------------------------------------------------------
6634 AliTPCseed * t = dynamic_cast<AliTPCseed*>(tk);
6636 printf("%s:%d wrong type \n",(char*)__FILE__,__LINE__);
6640 Int_t noc=t->GetNumberOfClusters();
6642 //printf("\nnot founded prolongation\n\n\n");
6648 AliTPCclusterMI *clusters[160];
6650 for (Int_t i=0;i<160;i++) {
6657 for (i=0; i<160 && current<noc; i++) {
6659 Int_t index=t->GetClusterIndex2(i);
6660 if (index<=0) continue;
6661 if (index&0x8000) continue;
6663 //clusters[current]=GetClusterMI(index);
6664 if (t->GetClusterPointer(i)){
6665 clusters[current]=t->GetClusterPointer(i);
6671 Int_t lab=123456789;
6672 for (i=0; i<noc; i++) {
6673 AliTPCclusterMI *c=clusters[i];
6675 lab=TMath::Abs(c->GetLabel(0));
6677 for (j=0; j<noc; j++) if (lb[j]==lab || mx[j]==0) break;
6683 for (i=0; i<noc; i++) if (mx[i]>max) {max=mx[i]; lab=lb[i];}
6685 for (i=0; i<noc; i++) {
6686 AliTPCclusterMI *c=clusters[i];
6688 if (TMath::Abs(c->GetLabel(1)) == lab ||
6689 TMath::Abs(c->GetLabel(2)) == lab ) max++;
6691 if (noc<=0) { lab=-1; return;}
6692 if ((1.- Float_t(max)/(noc)) > wrong) lab=-lab;
6695 Int_t tail=Int_t(0.10*noc);
6698 for (i=1; i<160&&ind<tail; i++) {
6699 // AliTPCclusterMI *c=clusters[noc-i];
6700 AliTPCclusterMI *c=clusters[i];
6702 if (lab == TMath::Abs(c->GetLabel(0)) ||
6703 lab == TMath::Abs(c->GetLabel(1)) ||
6704 lab == TMath::Abs(c->GetLabel(2))) max++;
6707 if (max < Int_t(0.5*tail)) lab=-lab;
6714 //delete[] clusters;
6718 //__________________________________________________________________________
6719 Int_t AliTPCtrackerMI::CookLabel(AliTPCseed *const t, Float_t wrong,Int_t first, Int_t last) const {
6720 //--------------------------------------------------------------------
6721 //This function "cooks" a track label. If label<0, this track is fake.
6722 //--------------------------------------------------------------------
6723 Int_t noc=t->GetNumberOfClusters();
6725 //printf("\nnot founded prolongation\n\n\n");
6731 AliTPCclusterMI *clusters[160];
6733 for (Int_t i=0;i<160;i++) {
6740 for (i=0; i<160 && current<noc; i++) {
6741 if (i<first) continue;
6742 if (i>last) continue;
6743 Int_t index=t->GetClusterIndex2(i);
6744 if (index<=0) continue;
6745 if (index&0x8000) continue;
6747 //clusters[current]=GetClusterMI(index);
6748 if (t->GetClusterPointer(i)){
6749 clusters[current]=t->GetClusterPointer(i);
6754 //if (noc<5) return -1;
6755 Int_t lab=123456789;
6756 for (i=0; i<noc; i++) {
6757 AliTPCclusterMI *c=clusters[i];
6759 lab=TMath::Abs(c->GetLabel(0));
6761 for (j=0; j<noc; j++) if (lb[j]==lab || mx[j]==0) break;
6767 for (i=0; i<noc; i++) if (mx[i]>max) {max=mx[i]; lab=lb[i];}
6769 for (i=0; i<noc; i++) {
6770 AliTPCclusterMI *c=clusters[i];
6772 if (TMath::Abs(c->GetLabel(1)) == lab ||
6773 TMath::Abs(c->GetLabel(2)) == lab ) max++;
6775 if (noc<=0) { lab=-1; return -1;}
6776 if ((1.- Float_t(max)/(noc)) > wrong) lab=-lab;
6779 Int_t tail=Int_t(0.10*noc);
6782 for (i=1; i<160&&ind<tail; i++) {
6783 // AliTPCclusterMI *c=clusters[noc-i];
6784 AliTPCclusterMI *c=clusters[i];
6786 if (lab == TMath::Abs(c->GetLabel(0)) ||
6787 lab == TMath::Abs(c->GetLabel(1)) ||
6788 lab == TMath::Abs(c->GetLabel(2))) max++;
6791 if (max < Int_t(0.5*tail)) lab=-lab;
6794 // t->SetLabel(lab);
6798 //delete[] clusters;
6802 Int_t AliTPCtrackerMI::GetRowNumber(Double_t x[3]) const
6804 //return pad row number for given x vector
6805 Float_t phi = TMath::ATan2(x[1],x[0]);
6806 if(phi<0) phi=2.*TMath::Pi()+phi;
6807 // Get the local angle in the sector philoc
6808 const Float_t kRaddeg = 180/3.14159265358979312;
6809 Float_t phiangle = (Int_t (phi*kRaddeg/20.) + 0.5)*20./kRaddeg;
6810 Double_t localx = x[0]*TMath::Cos(phiangle)-x[1]*TMath::Sin(phiangle);
6811 return GetRowNumber(localx);
6816 void AliTPCtrackerMI::MakeBitmaps(AliTPCseed *t)
6818 //-----------------------------------------------------------------------
6819 // Fill the cluster and sharing bitmaps of the track
6820 //-----------------------------------------------------------------------
6822 Int_t firstpoint = 0;
6823 Int_t lastpoint = 159;
6824 AliTPCTrackerPoint *point;
6825 AliTPCclusterMI *cluster;
6827 for (int iter=firstpoint; iter<lastpoint; iter++) {
6828 // Change to cluster pointers to see if we have a cluster at given padrow
6829 cluster = t->GetClusterPointer(iter);
6831 t->SetClusterMapBit(iter, kTRUE);
6832 point = t->GetTrackPoint(iter);
6833 if (point->IsShared())
6834 t->SetSharedMapBit(iter,kTRUE);
6836 t->SetSharedMapBit(iter, kFALSE);
6839 t->SetClusterMapBit(iter, kFALSE);
6840 t->SetSharedMapBit(iter, kFALSE);
6845 Bool_t AliTPCtrackerMI::IsFindable(AliTPCseed & track){
6847 // return flag if there is findable cluster at given position
6850 Float_t z = track.GetZ();
6852 if (TMath::Abs(z)<(AliTPCReconstructor::GetCtgRange()*track.GetX()+kDeltaZ) &&
6853 TMath::Abs(z)<fkParam->GetZLength(0) &&
6854 (TMath::Abs(track.GetSnp())<AliTPCReconstructor::GetMaxSnpTracker()))
6860 void AliTPCtrackerMI::AddCovariance(AliTPCseed * seed){
6862 // Adding systematic error estimate to the covariance matrix
6863 // !!!! the systematic error for element 4 is in 1/cm not in pt
6865 const Double_t *param = AliTPCReconstructor::GetRecoParam()->GetSystematicError();
6867 // use only the diagonal part if not specified otherwise
6868 if (!AliTPCReconstructor::GetRecoParam()->GetUseSystematicCorrelation()) return AddCovarianceAdd(seed);
6870 Double_t *covarS= (Double_t*)seed->GetCovariance();
6871 Double_t factor[5]={1,1,1,1,1};
6872 Double_t facC = AliTracker::GetBz()*kB2C;
6873 factor[0]= TMath::Sqrt(TMath::Abs((covarS[0] + param[0]*param[0])/covarS[0]));
6874 factor[1]= TMath::Sqrt(TMath::Abs((covarS[2] + param[1]*param[1])/covarS[2]));
6875 factor[2]= TMath::Sqrt(TMath::Abs((covarS[5] + param[2]*param[2])/covarS[5]));
6876 factor[3]= TMath::Sqrt(TMath::Abs((covarS[9] + param[3]*param[3])/covarS[9]));
6877 factor[4]= TMath::Sqrt(TMath::Abs((covarS[14] + facC*facC*param[4]*param[4])/covarS[14]));
6883 for (Int_t i=0; i<5; i++){
6884 for (Int_t j=i; j<5; j++){
6885 Int_t index=seed->GetIndex(i,j);
6886 covarS[index]*=factor[i]*factor[j];
6892 void AliTPCtrackerMI::AddCovarianceAdd(AliTPCseed * seed){
6894 // Adding systematic error - as additive factor without correlation
6896 const Double_t *param = AliTPCReconstructor::GetRecoParam()->GetSystematicError();
6897 Double_t *covarIn= (Double_t*)seed->GetCovariance();
6899 for (Int_t i=0;i<15;i++) covar[i]=0;
6905 covar[0] = param[0]*param[0];
6906 covar[2] = param[1]*param[1];
6907 covar[5] = param[2]*param[2];
6908 covar[9] = param[3]*param[3];
6909 Double_t facC = AliTracker::GetBz()*kB2C;
6910 covar[14]= param[4]*param[4]*facC*facC;
6912 covar[1]=TMath::Sqrt((covar[0]*covar[2]))*covarIn[1]/TMath::Sqrt((covarIn[0]*covarIn[2]));
6914 covar[3]=TMath::Sqrt((covar[0]*covar[5]))*covarIn[3]/TMath::Sqrt((covarIn[0]*covarIn[5]));
6915 covar[4]=TMath::Sqrt((covar[2]*covar[5]))*covarIn[4]/TMath::Sqrt((covarIn[2]*covarIn[5]));
6917 covar[6]=TMath::Sqrt((covar[0]*covar[9]))*covarIn[6]/TMath::Sqrt((covarIn[0]*covarIn[9]));
6918 covar[7]=TMath::Sqrt((covar[2]*covar[9]))*covarIn[7]/TMath::Sqrt((covarIn[2]*covarIn[9]));
6919 covar[8]=TMath::Sqrt((covar[5]*covar[9]))*covarIn[8]/TMath::Sqrt((covarIn[5]*covarIn[9]));
6921 covar[10]=TMath::Sqrt((covar[0]*covar[14]))*covarIn[10]/TMath::Sqrt((covarIn[0]*covarIn[14]));
6922 covar[11]=TMath::Sqrt((covar[2]*covar[14]))*covarIn[11]/TMath::Sqrt((covarIn[2]*covarIn[14]));
6923 covar[12]=TMath::Sqrt((covar[5]*covar[14]))*covarIn[12]/TMath::Sqrt((covarIn[5]*covarIn[14]));
6924 covar[13]=TMath::Sqrt((covar[9]*covar[14]))*covarIn[13]/TMath::Sqrt((covarIn[9]*covarIn[14]));
6926 seed->AddCovariance(covar);