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>
117 #include "AliComplexCluster.h"
118 #include "AliESDEvent.h"
119 #include "AliESDtrack.h"
120 #include "AliESDVertex.h"
123 #include "AliHelix.h"
124 #include "AliRunLoader.h"
125 #include "AliTPCClustersRow.h"
126 #include "AliTPCParam.h"
127 #include "AliTPCReconstructor.h"
128 #include "AliTPCpolyTrack.h"
129 #include "AliTPCreco.h"
130 #include "AliTPCseed.h"
132 #include "AliTPCtrackerSector.h"
133 #include "AliTPCtrackerMI.h"
134 #include "TStopwatch.h"
135 #include "AliTPCReconstructor.h"
136 #include "AliAlignObj.h"
137 #include "AliTrackPointArray.h"
139 #include "AliTPCcalibDB.h"
140 #include "AliTPCTransform.h"
141 #include "AliTPCClusterParam.h"
145 ClassImp(AliTPCtrackerMI)
149 class AliTPCFastMath {
152 static Double_t FastAsin(Double_t x);
154 static Double_t fgFastAsin[20000]; //lookup table for fast asin computation
157 Double_t AliTPCFastMath::fgFastAsin[20000];
158 AliTPCFastMath gAliTPCFastMath; // needed to fill the LUT
160 AliTPCFastMath::AliTPCFastMath(){
162 // initialized lookup table;
163 for (Int_t i=0;i<10000;i++){
164 fgFastAsin[2*i] = TMath::ASin(i/10000.);
165 fgFastAsin[2*i+1] = (TMath::ASin((i+1)/10000.)-fgFastAsin[2*i]);
169 Double_t AliTPCFastMath::FastAsin(Double_t x){
171 // return asin using lookup table
173 Int_t index = int(x*10000);
174 return fgFastAsin[2*index]+(x*10000.-index)*fgFastAsin[2*index+1];
177 Int_t index = int(x*10000);
178 return -(fgFastAsin[2*index]+(x*10000.-index)*fgFastAsin[2*index+1]);
180 //__________________________________________________________________
181 AliTPCtrackerMI::AliTPCtrackerMI()
203 // default constructor
206 //_____________________________________________________________________
210 Int_t AliTPCtrackerMI::UpdateTrack(AliTPCseed * track, Int_t accept){
212 //update track information using current cluster - track->fCurrentCluster
215 AliTPCclusterMI* c =track->GetCurrentCluster();
216 if (accept>0) track->SetCurrentClusterIndex1(track->GetCurrentClusterIndex1() | 0x8000); //sign not accepted clusters
218 UInt_t i = track->GetCurrentClusterIndex1();
220 Int_t sec=(i&0xff000000)>>24;
221 //Int_t row = (i&0x00ff0000)>>16;
222 track->SetRow((i&0x00ff0000)>>16);
223 track->SetSector(sec);
224 // Int_t index = i&0xFFFF;
225 if (sec>=fkParam->GetNInnerSector()) track->SetRow(track->GetRow()+fkParam->GetNRowLow());
226 track->SetClusterIndex2(track->GetRow(), i);
227 //track->fFirstPoint = row;
228 //if ( track->fLastPoint<row) track->fLastPoint =row;
229 // if (track->fRow<0 || track->fRow>160) {
230 // printf("problem\n");
232 if (track->GetFirstPoint()>track->GetRow())
233 track->SetFirstPoint(track->GetRow());
234 if (track->GetLastPoint()<track->GetRow())
235 track->SetLastPoint(track->GetRow());
238 track->SetClusterPointer(track->GetRow(),c);
241 Double_t angle2 = track->GetSnp()*track->GetSnp();
243 //SET NEW Track Point
245 if (angle2<1) //PH sometimes angle2 is very big. To be investigated...
247 angle2 = TMath::Sqrt(angle2/(1-angle2));
248 AliTPCTrackerPoint &point =*(track->GetTrackPoint(track->GetRow()));
250 point.SetSigmaY(c->GetSigmaY2()/track->GetCurrentSigmaY2());
251 point.SetSigmaZ(c->GetSigmaZ2()/track->GetCurrentSigmaZ2());
252 point.SetErrY(sqrt(track->GetErrorY2()));
253 point.SetErrZ(sqrt(track->GetErrorZ2()));
255 point.SetX(track->GetX());
256 point.SetY(track->GetY());
257 point.SetZ(track->GetZ());
258 point.SetAngleY(angle2);
259 point.SetAngleZ(track->GetTgl());
260 if (point.IsShared()){
261 track->SetErrorY2(track->GetErrorY2()*4);
262 track->SetErrorZ2(track->GetErrorZ2()*4);
266 Double_t chi2 = track->GetPredictedChi2(track->GetCurrentCluster());
268 // track->SetErrorY2(track->GetErrorY2()*1.3);
269 // track->SetErrorY2(track->GetErrorY2()+0.01);
270 // track->SetErrorZ2(track->GetErrorZ2()*1.3);
271 // track->SetErrorZ2(track->GetErrorZ2()+0.005);
273 if (accept>0) return 0;
274 if (track->GetNumberOfClusters()%20==0){
275 // if (track->fHelixIn){
276 // TClonesArray & larr = *(track->fHelixIn);
277 // Int_t ihelix = larr.GetEntriesFast();
278 // new(larr[ihelix]) AliHelix(*track) ;
281 track->SetNoCluster(0);
282 return track->Update(c,chi2,i);
287 Int_t AliTPCtrackerMI::AcceptCluster(AliTPCseed * seed, AliTPCclusterMI * cluster)
290 // decide according desired precision to accept given
291 // cluster for tracking
293 seed->GetProlongation(cluster->GetX(),yt,zt);
294 Double_t sy2=ErrY2(seed,cluster);
295 Double_t sz2=ErrZ2(seed,cluster);
297 Double_t sdistancey2 = sy2+seed->GetSigmaY2();
298 Double_t sdistancez2 = sz2+seed->GetSigmaZ2();
299 Double_t dy=seed->GetCurrentCluster()->GetY()-yt;
300 Double_t dz=seed->GetCurrentCluster()->GetY()-zt;
301 Double_t rdistancey2 = (seed->GetCurrentCluster()->GetY()-yt)*
302 (seed->GetCurrentCluster()->GetY()-yt)/sdistancey2;
303 Double_t rdistancez2 = (seed->GetCurrentCluster()->GetZ()-zt)*
304 (seed->GetCurrentCluster()->GetZ()-zt)/sdistancez2;
306 Double_t rdistance2 = rdistancey2+rdistancez2;
309 if (AliTPCReconstructor::StreamLevel()>2 && seed->GetNumberOfClusters()>20) {
310 Float_t rmsy2 = seed->GetCurrentSigmaY2();
311 Float_t rmsz2 = seed->GetCurrentSigmaZ2();
312 Float_t rmsy2p30 = seed->GetCMeanSigmaY2p30();
313 Float_t rmsz2p30 = seed->GetCMeanSigmaZ2p30();
314 Float_t rmsy2p30R = seed->GetCMeanSigmaY2p30R();
315 Float_t rmsz2p30R = seed->GetCMeanSigmaZ2p30R();
316 AliExternalTrackParam param(*seed);
317 static TVectorD gcl(3),gtr(3);
319 param.GetXYZ(gcl.GetMatrixArray());
320 cluster->GetGlobalXYZ(gclf);
321 gcl[0]=gclf[0]; gcl[1]=gclf[1]; gcl[2]=gclf[2];
324 if (AliTPCReconstructor::StreamLevel()>2) {
325 (*fDebugStreamer)<<"ErrParam"<<
338 "rmsy2p30="<<rmsy2p30<<
339 "rmsz2p30="<<rmsz2p30<<
340 "rmsy2p30R="<<rmsy2p30R<<
341 "rmsz2p30R="<<rmsz2p30R<<
345 //return 0; // temporary
346 if (rdistance2>32) return 3;
349 if ((rdistancey2>9. || rdistancez2>9.) && cluster->GetType()==0)
350 return 2; //suspisiouce - will be changed
352 if ((rdistancey2>6.25 || rdistancez2>6.25) && cluster->GetType()>0)
353 // strict cut on overlaped cluster
354 return 2; //suspisiouce - will be changed
356 if ( (rdistancey2>1. || rdistancez2>6.25 )
357 && cluster->GetType()<0){
358 seed->SetNFoundable(seed->GetNFoundable()-1);
368 //_____________________________________________________________________________
369 AliTPCtrackerMI::AliTPCtrackerMI(const AliTPCParam *par):
371 fkNIS(par->GetNInnerSector()/2),
373 fkNOS(par->GetNOuterSector()/2),
390 //---------------------------------------------------------------------
391 // The main TPC tracker constructor
392 //---------------------------------------------------------------------
393 fInnerSec=new AliTPCtrackerSector[fkNIS];
394 fOuterSec=new AliTPCtrackerSector[fkNOS];
397 for (i=0; i<fkNIS; i++) fInnerSec[i].Setup(par,0);
398 for (i=0; i<fkNOS; i++) fOuterSec[i].Setup(par,1);
401 Int_t nrowlow = par->GetNRowLow();
402 Int_t nrowup = par->GetNRowUp();
405 for (i=0;i<nrowlow;i++){
406 fXRow[i] = par->GetPadRowRadiiLow(i);
407 fPadLength[i]= par->GetPadPitchLength(0,i);
408 fYMax[i] = fXRow[i]*TMath::Tan(0.5*par->GetInnerAngle());
412 for (i=0;i<nrowup;i++){
413 fXRow[i+nrowlow] = par->GetPadRowRadiiUp(i);
414 fPadLength[i+nrowlow] = par->GetPadPitchLength(60,i);
415 fYMax[i+nrowlow] = fXRow[i+nrowlow]*TMath::Tan(0.5*par->GetOuterAngle());
418 if (AliTPCReconstructor::StreamLevel()>0) {
419 fDebugStreamer = new TTreeSRedirector("TPCdebug.root");
422 //________________________________________________________________________
423 AliTPCtrackerMI::AliTPCtrackerMI(const AliTPCtrackerMI &t):
444 //------------------------------------
445 // dummy copy constructor
446 //------------------------------------------------------------------
449 AliTPCtrackerMI & AliTPCtrackerMI::operator=(const AliTPCtrackerMI& /*r*/)
451 //------------------------------
453 //--------------------------------------------------------------
456 //_____________________________________________________________________________
457 AliTPCtrackerMI::~AliTPCtrackerMI() {
458 //------------------------------------------------------------------
459 // TPC tracker destructor
460 //------------------------------------------------------------------
467 if (fDebugStreamer) delete fDebugStreamer;
471 void AliTPCtrackerMI::FillESD(const TObjArray* arr)
475 //fill esds using updated tracks
477 // write tracks to the event
478 // store index of the track
479 Int_t nseed=arr->GetEntriesFast();
480 //FindKinks(arr,fEvent);
481 for (Int_t i=0; i<nseed; i++) {
482 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
486 if (AliTPCReconstructor::StreamLevel()>1) {
487 (*fDebugStreamer)<<"Track0"<<
491 // pt->PropagateTo(fkParam->GetInnerRadiusLow());
492 if (pt->GetKinkIndex(0)<=0){ //don't propagate daughter tracks
493 pt->PropagateTo(fkParam->GetInnerRadiusLow());
496 if (( pt->GetPoints()[2]- pt->GetPoints()[0])>5 && pt->GetPoints()[3]>0.8){
498 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
499 iotrack.SetTPCPoints(pt->GetPoints());
500 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
501 iotrack.SetV0Indexes(pt->GetV0Indexes());
502 // iotrack.SetTPCpid(pt->fTPCr);
503 //iotrack.SetTPCindex(i);
504 fEvent->AddTrack(&iotrack);
508 if ( (pt->GetNumberOfClusters()>70)&& (Float_t(pt->GetNumberOfClusters())/Float_t(pt->GetNFoundable()))>0.55) {
510 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
511 iotrack.SetTPCPoints(pt->GetPoints());
512 //iotrack.SetTPCindex(i);
513 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
514 iotrack.SetV0Indexes(pt->GetV0Indexes());
515 // iotrack.SetTPCpid(pt->fTPCr);
516 fEvent->AddTrack(&iotrack);
520 // short tracks - maybe decays
522 if ( (pt->GetNumberOfClusters()>30) && (Float_t(pt->GetNumberOfClusters())/Float_t(pt->GetNFoundable()))>0.70) {
523 Int_t found,foundable,shared;
524 pt->GetClusterStatistic(0,60,found, foundable,shared,kFALSE);
525 if ( (found>20) && (pt->GetNShared()/float(pt->GetNumberOfClusters())<0.2)){
527 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
528 //iotrack.SetTPCindex(i);
529 iotrack.SetTPCPoints(pt->GetPoints());
530 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
531 iotrack.SetV0Indexes(pt->GetV0Indexes());
532 //iotrack.SetTPCpid(pt->fTPCr);
533 fEvent->AddTrack(&iotrack);
538 if ( (pt->GetNumberOfClusters()>20) && (Float_t(pt->GetNumberOfClusters())/Float_t(pt->GetNFoundable()))>0.8) {
539 Int_t found,foundable,shared;
540 pt->GetClusterStatistic(0,60,found, foundable,shared,kFALSE);
541 if (found<20) continue;
542 if (pt->GetNShared()/float(pt->GetNumberOfClusters())>0.2) continue;
545 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
546 iotrack.SetTPCPoints(pt->GetPoints());
547 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
548 iotrack.SetV0Indexes(pt->GetV0Indexes());
549 //iotrack.SetTPCpid(pt->fTPCr);
550 //iotrack.SetTPCindex(i);
551 fEvent->AddTrack(&iotrack);
554 // short tracks - secondaties
556 if ( (pt->GetNumberOfClusters()>30) ) {
557 Int_t found,foundable,shared;
558 pt->GetClusterStatistic(128,158,found, foundable,shared,kFALSE);
559 if ( (found>20) && (pt->GetNShared()/float(pt->GetNumberOfClusters())<0.2) &&float(found)/float(foundable)>0.8){
561 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
562 iotrack.SetTPCPoints(pt->GetPoints());
563 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
564 iotrack.SetV0Indexes(pt->GetV0Indexes());
565 //iotrack.SetTPCpid(pt->fTPCr);
566 //iotrack.SetTPCindex(i);
567 fEvent->AddTrack(&iotrack);
572 if ( (pt->GetNumberOfClusters()>15)) {
573 Int_t found,foundable,shared;
574 pt->GetClusterStatistic(138,158,found, foundable,shared,kFALSE);
575 if (found<15) continue;
576 if (foundable<=0) continue;
577 if (pt->GetNShared()/float(pt->GetNumberOfClusters())>0.2) continue;
578 if (float(found)/float(foundable)<0.8) continue;
581 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
582 iotrack.SetTPCPoints(pt->GetPoints());
583 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
584 iotrack.SetV0Indexes(pt->GetV0Indexes());
585 // iotrack.SetTPCpid(pt->fTPCr);
586 //iotrack.SetTPCindex(i);
587 fEvent->AddTrack(&iotrack);
591 // >> account for suppressed tracks in the kink indices (RS)
592 int nESDtracks = fEvent->GetNumberOfTracks();
593 for (int it=nESDtracks;it--;) {
594 AliESDtrack* esdTr = fEvent->GetTrack(it);
595 if (!esdTr || !esdTr->GetKinkIndex(0)) continue;
596 for (int ik=0;ik<3;ik++) {
598 if (!(knkId=esdTr->GetKinkIndex(ik))) break; // no more kinks for this track
599 AliESDkink* kink = fEvent->GetKink(TMath::Abs(knkId)-1);
601 AliError(Form("ESDTrack%d refers to non-existing kink %d",it,TMath::Abs(knkId)-1));
604 kink->SetIndex(it, knkId<0 ? 0:1); // update track index of the kink: mother at 0, daughter at 1
607 // << account for suppressed tracks in the kink indices (RS)
609 AliInfo(Form("Number of filled ESDs-\t%d\n",fEvent->GetNumberOfTracks()));
616 Double_t AliTPCtrackerMI::ErrY2(AliTPCseed* seed, const AliTPCclusterMI * cl){
619 // Use calibrated cluster error from OCDB
621 AliTPCClusterParam * clparam = AliTPCcalibDB::Instance()->GetClusterParam();
623 Float_t z = TMath::Abs(fkParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
624 Int_t ctype = cl->GetType();
625 Int_t type = (cl->GetRow()<63) ? 0: (cl->GetRow()>126) ? 1:2;
626 Double_t angle = seed->GetSnp()*seed->GetSnp();
627 angle = TMath::Sqrt(TMath::Abs(angle/(1.-angle)));
628 Double_t erry2 = clparam->GetError0Par(0,type, z,angle);
630 erry2+=0.5; // edge cluster
633 seed->SetErrorY2(erry2);
637 //calculate look-up table at the beginning
638 // static Bool_t ginit = kFALSE;
639 // static Float_t gnoise1,gnoise2,gnoise3;
640 // static Float_t ggg1[10000];
641 // static Float_t ggg2[10000];
642 // static Float_t ggg3[10000];
643 // static Float_t glandau1[10000];
644 // static Float_t glandau2[10000];
645 // static Float_t glandau3[10000];
647 // static Float_t gcor01[500];
648 // static Float_t gcor02[500];
649 // static Float_t gcorp[500];
653 // if (ginit==kFALSE){
654 // for (Int_t i=1;i<500;i++){
655 // Float_t rsigma = float(i)/100.;
656 // gcor02[i] = TMath::Max(0.78 +TMath::Exp(7.4*(rsigma-1.2)),0.6);
657 // gcor01[i] = TMath::Max(0.72 +TMath::Exp(3.36*(rsigma-1.2)),0.6);
658 // gcorp[i] = TMath::Max(TMath::Power((rsigma+0.5),1.5),1.2);
662 // for (Int_t i=3;i<10000;i++){
666 // Float_t amp = float(i);
667 // Float_t padlength =0.75;
668 // gnoise1 = 0.0004/padlength;
669 // Float_t nel = 0.268*amp;
670 // Float_t nprim = 0.155*amp;
671 // ggg1[i] = fkParam->GetDiffT()*fkParam->GetDiffT()*(2+0.001*nel/(padlength*padlength))/nel;
672 // glandau1[i] = (2.+0.12*nprim)*0.5* (2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
673 // if (glandau1[i]>1) glandau1[i]=1;
674 // glandau1[i]*=padlength*padlength/12.;
678 // gnoise2 = 0.0004/padlength;
680 // nprim = 0.133*amp;
681 // ggg2[i] = fkParam->GetDiffT()*fkParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
682 // glandau2[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
683 // if (glandau2[i]>1) glandau2[i]=1;
684 // glandau2[i]*=padlength*padlength/12.;
689 // gnoise3 = 0.0004/padlength;
691 // nprim = 0.133*amp;
692 // ggg3[i] = fkParam->GetDiffT()*fkParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
693 // glandau3[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
694 // if (glandau3[i]>1) glandau3[i]=1;
695 // glandau3[i]*=padlength*padlength/12.;
703 // Int_t amp = int(TMath::Abs(cl->GetQ()));
705 // seed->SetErrorY2(1.);
709 // Float_t z = TMath::Abs(fkParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
710 // Int_t ctype = cl->GetType();
711 // Float_t padlength= GetPadPitchLength(seed->GetRow());
712 // Double_t angle2 = seed->GetSnp()*seed->GetSnp();
713 // angle2 = angle2/(1-angle2);
715 // //cluster "quality"
716 // Int_t rsigmay = int(100.*cl->GetSigmaY2()/(seed->GetCurrentSigmaY2()));
719 // if (fSectors==fInnerSec){
720 // snoise2 = gnoise1;
721 // res = ggg1[amp]*z+glandau1[amp]*angle2;
722 // if (ctype==0) res *= gcor01[rsigmay];
725 // res*= gcorp[rsigmay];
729 // if (padlength<1.1){
730 // snoise2 = gnoise2;
731 // res = ggg2[amp]*z+glandau2[amp]*angle2;
732 // if (ctype==0) res *= gcor02[rsigmay];
735 // res*= gcorp[rsigmay];
739 // snoise2 = gnoise3;
740 // res = ggg3[amp]*z+glandau3[amp]*angle2;
741 // if (ctype==0) res *= gcor02[rsigmay];
744 // res*= gcorp[rsigmay];
751 // res*=2.4; // overestimate error 2 times
755 // if (res<2*snoise2)
758 // seed->SetErrorY2(res);
766 Double_t AliTPCtrackerMI::ErrZ2(AliTPCseed* seed, const AliTPCclusterMI * cl){
769 // Use calibrated cluster error from OCDB
771 AliTPCClusterParam * clparam = AliTPCcalibDB::Instance()->GetClusterParam();
773 Float_t z = TMath::Abs(fkParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
774 Int_t ctype = cl->GetType();
775 Int_t type = (cl->GetRow()<63) ? 0: (cl->GetRow()>126) ? 1:2;
777 Double_t angle2 = seed->GetSnp()*seed->GetSnp();
778 angle2 = seed->GetTgl()*seed->GetTgl()*(1+angle2/(1-angle2));
779 Double_t angle = TMath::Sqrt(TMath::Abs(angle2));
780 Double_t errz2 = clparam->GetError0Par(1,type, z,angle);
782 errz2+=0.5; // edge cluster
785 seed->SetErrorZ2(errz2);
791 // //seed->SetErrorY2(0.1);
793 // //calculate look-up table at the beginning
794 // static Bool_t ginit = kFALSE;
795 // static Float_t gnoise1,gnoise2,gnoise3;
796 // static Float_t ggg1[10000];
797 // static Float_t ggg2[10000];
798 // static Float_t ggg3[10000];
799 // static Float_t glandau1[10000];
800 // static Float_t glandau2[10000];
801 // static Float_t glandau3[10000];
803 // static Float_t gcor01[1000];
804 // static Float_t gcor02[1000];
805 // static Float_t gcorp[1000];
809 // if (ginit==kFALSE){
810 // for (Int_t i=1;i<1000;i++){
811 // Float_t rsigma = float(i)/100.;
812 // gcor02[i] = TMath::Max(0.81 +TMath::Exp(6.8*(rsigma-1.2)),0.6);
813 // gcor01[i] = TMath::Max(0.72 +TMath::Exp(2.04*(rsigma-1.2)),0.6);
814 // gcorp[i] = TMath::Max(TMath::Power((rsigma+0.5),1.5),1.2);
818 // for (Int_t i=3;i<10000;i++){
822 // Float_t amp = float(i);
823 // Float_t padlength =0.75;
824 // gnoise1 = 0.0004/padlength;
825 // Float_t nel = 0.268*amp;
826 // Float_t nprim = 0.155*amp;
827 // ggg1[i] = fkParam->GetDiffT()*fkParam->GetDiffT()*(2+0.001*nel/(padlength*padlength))/nel;
828 // glandau1[i] = (2.+0.12*nprim)*0.5* (2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
829 // if (glandau1[i]>1) glandau1[i]=1;
830 // glandau1[i]*=padlength*padlength/12.;
834 // gnoise2 = 0.0004/padlength;
836 // nprim = 0.133*amp;
837 // ggg2[i] = fkParam->GetDiffT()*fkParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
838 // glandau2[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
839 // if (glandau2[i]>1) glandau2[i]=1;
840 // glandau2[i]*=padlength*padlength/12.;
845 // gnoise3 = 0.0004/padlength;
847 // nprim = 0.133*amp;
848 // ggg3[i] = fkParam->GetDiffT()*fkParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
849 // glandau3[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
850 // if (glandau3[i]>1) glandau3[i]=1;
851 // glandau3[i]*=padlength*padlength/12.;
859 // Int_t amp = int(TMath::Abs(cl->GetQ()));
861 // seed->SetErrorY2(1.);
865 // Float_t z = TMath::Abs(fkParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
866 // Int_t ctype = cl->GetType();
867 // Float_t padlength= GetPadPitchLength(seed->GetRow());
869 // Double_t angle2 = seed->GetSnp()*seed->GetSnp();
870 // // if (angle2<0.6) angle2 = 0.6;
871 // angle2 = seed->GetTgl()*seed->GetTgl()*(1+angle2/(1-angle2));
873 // //cluster "quality"
874 // Int_t rsigmaz = int(100.*cl->GetSigmaZ2()/(seed->GetCurrentSigmaZ2()));
877 // if (fSectors==fInnerSec){
878 // snoise2 = gnoise1;
879 // res = ggg1[amp]*z+glandau1[amp]*angle2;
880 // if (ctype==0) res *= gcor01[rsigmaz];
883 // res*= gcorp[rsigmaz];
887 // if (padlength<1.1){
888 // snoise2 = gnoise2;
889 // res = ggg2[amp]*z+glandau2[amp]*angle2;
890 // if (ctype==0) res *= gcor02[rsigmaz];
893 // res*= gcorp[rsigmaz];
897 // snoise2 = gnoise3;
898 // res = ggg3[amp]*z+glandau3[amp]*angle2;
899 // if (ctype==0) res *= gcor02[rsigmaz];
902 // res*= gcorp[rsigmaz];
911 // if ((ctype<0) &&<70){
916 // if (res<2*snoise2)
918 // if (res>3) res =3;
919 // seed->SetErrorZ2(res);
927 void AliTPCtrackerMI::RotateToLocal(AliTPCseed *seed)
929 //rotate to track "local coordinata
930 Float_t x = seed->GetX();
931 Float_t y = seed->GetY();
932 Float_t ymax = x*TMath::Tan(0.5*fSectors->GetAlpha());
935 seed->SetRelativeSector((seed->GetRelativeSector()+1) % fN);
936 if (!seed->Rotate(fSectors->GetAlpha()))
938 } else if (y <-ymax) {
939 seed->SetRelativeSector((seed->GetRelativeSector()-1+fN) % fN);
940 if (!seed->Rotate(-fSectors->GetAlpha()))
948 //_____________________________________________________________________________
949 Double_t AliTPCtrackerMI::F1old(Double_t x1,Double_t y1,
950 Double_t x2,Double_t y2,
951 Double_t x3,Double_t y3) const
953 //-----------------------------------------------------------------
954 // Initial approximation of the track curvature
955 //-----------------------------------------------------------------
956 Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
957 Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
958 (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
959 Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
960 (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
962 Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
963 if ( xr*xr+yr*yr<=0.00000000000001) return 100;
964 return -xr*yr/sqrt(xr*xr+yr*yr);
969 //_____________________________________________________________________________
970 Double_t AliTPCtrackerMI::F1(Double_t x1,Double_t y1,
971 Double_t x2,Double_t y2,
972 Double_t x3,Double_t y3) const
974 //-----------------------------------------------------------------
975 // Initial approximation of the track curvature
976 //-----------------------------------------------------------------
982 Double_t det = x3*y2-x2*y3;
983 if (TMath::Abs(det)<1e-10){
987 Double_t u = 0.5* (x2*(x2-x3)+y2*(y2-y3))/det;
988 Double_t x0 = x3*0.5-y3*u;
989 Double_t y0 = y3*0.5+x3*u;
990 Double_t c2 = 1/TMath::Sqrt(x0*x0+y0*y0);
996 Double_t AliTPCtrackerMI::F2(Double_t x1,Double_t y1,
997 Double_t x2,Double_t y2,
998 Double_t x3,Double_t y3) const
1000 //-----------------------------------------------------------------
1001 // Initial approximation of the track curvature
1002 //-----------------------------------------------------------------
1008 Double_t det = x3*y2-x2*y3;
1009 if (TMath::Abs(det)<1e-10) {
1013 Double_t u = 0.5* (x2*(x2-x3)+y2*(y2-y3))/det;
1014 Double_t x0 = x3*0.5-y3*u;
1015 Double_t y0 = y3*0.5+x3*u;
1016 Double_t c2 = 1/TMath::Sqrt(x0*x0+y0*y0);
1025 //_____________________________________________________________________________
1026 Double_t AliTPCtrackerMI::F2old(Double_t x1,Double_t y1,
1027 Double_t x2,Double_t y2,
1028 Double_t x3,Double_t y3) const
1030 //-----------------------------------------------------------------
1031 // Initial approximation of the track curvature times center of curvature
1032 //-----------------------------------------------------------------
1033 Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
1034 Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
1035 (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
1036 Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
1037 (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
1039 Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
1041 return -a/(d*y1-b)*xr/sqrt(xr*xr+yr*yr);
1044 //_____________________________________________________________________________
1045 Double_t AliTPCtrackerMI::F3(Double_t x1,Double_t y1,
1046 Double_t x2,Double_t y2,
1047 Double_t z1,Double_t z2) const
1049 //-----------------------------------------------------------------
1050 // Initial approximation of the tangent of the track dip angle
1051 //-----------------------------------------------------------------
1052 return (z1 - z2)/sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
1056 Double_t AliTPCtrackerMI::F3n(Double_t x1,Double_t y1,
1057 Double_t x2,Double_t y2,
1058 Double_t z1,Double_t z2, Double_t c) const
1060 //-----------------------------------------------------------------
1061 // Initial approximation of the tangent of the track dip angle
1062 //-----------------------------------------------------------------
1066 //angle1 = (z1-z2)*c/(TMath::ASin(c*x1-ni)-TMath::ASin(c*x2-ni));
1068 Double_t d = TMath::Sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
1069 if (TMath::Abs(d*c*0.5)>1) return 0;
1070 // Double_t angle2 = TMath::ASin(d*c*0.5);
1071 // Double_t angle2 = AliTPCFastMath::FastAsin(d*c*0.5);
1072 Double_t angle2 = (d*c*0.5>0.1)? TMath::ASin(d*c*0.5): AliTPCFastMath::FastAsin(d*c*0.5);
1074 angle2 = (z1-z2)*c/(angle2*2.);
1078 Bool_t AliTPCtrackerMI::GetProlongation(Double_t x1, Double_t x2, Double_t x[5], Double_t &y, Double_t &z) const
1079 {//-----------------------------------------------------------------
1080 // This function find proloncation of a track to a reference plane x=x2.
1081 //-----------------------------------------------------------------
1085 if (TMath::Abs(x[4]*x1 - x[2]) >= 0.999) {
1089 Double_t c1=x[4]*x1 - x[2], r1=TMath::Sqrt((1.-c1)*(1.+c1));
1090 Double_t c2=x[4]*x2 - x[2], r2=TMath::Sqrt((1.-c2)*(1.+c2));
1094 Double_t dy = dx*(c1+c2)/(r1+r2);
1097 Double_t delta = x[4]*dx*(c1+c2)/(c1*r2 + c2*r1);
1099 if (TMath::Abs(delta)>0.01){
1100 dz = x[3]*TMath::ASin(delta)/x[4];
1102 dz = x[3]*AliTPCFastMath::FastAsin(delta)/x[4];
1105 //dz = x[3]*AliTPCFastMath::FastAsin(delta)/x[4];
1113 Int_t AliTPCtrackerMI::LoadClusters (TTree *const tree)
1118 return LoadClusters();
1122 Int_t AliTPCtrackerMI::LoadClusters(const TObjArray *arr)
1125 // load clusters to the memory
1126 AliTPCClustersRow *clrow = 0x0;
1127 Int_t lower = arr->LowerBound();
1128 Int_t entries = arr->GetEntriesFast();
1130 for (Int_t i=lower; i<entries; i++) {
1131 clrow = (AliTPCClustersRow*) arr->At(i);
1132 if(!clrow) continue;
1133 if(!clrow->GetArray()) continue;
1137 fkParam->AdjustSectorRow(clrow->GetID(),sec,row);
1139 for (Int_t icl=0; icl<clrow->GetArray()->GetEntriesFast(); icl++){
1140 Transform((AliTPCclusterMI*)(clrow->GetArray()->At(icl)));
1143 if (clrow->GetArray()->GetEntriesFast()<=0) continue;
1144 AliTPCtrackerRow * tpcrow=0;
1147 tpcrow = &(fInnerSec[sec%fkNIS][row]);
1151 tpcrow = &(fOuterSec[(sec-fkNIS*2)%fkNOS][row]);
1152 left = (sec-fkNIS*2)/fkNOS;
1155 tpcrow->SetN1(clrow->GetArray()->GetEntriesFast());
1156 tpcrow->SetClusters1(new AliTPCclusterMI[tpcrow->GetN1()]);
1157 for (Int_t j=0;j<tpcrow->GetN1();++j)
1158 tpcrow->SetCluster1(j, *(AliTPCclusterMI*)(clrow->GetArray()->At(j)));
1161 tpcrow->SetN2(clrow->GetArray()->GetEntriesFast());
1162 tpcrow->SetClusters2(new AliTPCclusterMI[tpcrow->GetN2()]);
1163 for (Int_t j=0;j<tpcrow->GetN2();++j)
1164 tpcrow->SetCluster2(j,*(AliTPCclusterMI*)(clrow->GetArray()->At(j)));
1174 Int_t AliTPCtrackerMI::LoadClusters(const TClonesArray *arr)
1177 // load clusters to the memory from one
1180 AliTPCclusterMI *clust=0;
1181 Int_t count[72][96] = { {0} , {0} };
1183 // loop over clusters
1184 for (Int_t icl=0; icl<arr->GetEntriesFast(); icl++) {
1185 clust = (AliTPCclusterMI*)arr->At(icl);
1186 if(!clust) continue;
1187 //printf("cluster: det %d, row %d \n", clust->GetDetector(),clust->GetRow());
1189 // transform clusters
1192 // count clusters per pad row
1193 count[clust->GetDetector()][clust->GetRow()]++;
1196 // insert clusters to sectors
1197 for (Int_t icl=0; icl<arr->GetEntriesFast(); icl++) {
1198 clust = (AliTPCclusterMI*)arr->At(icl);
1199 if(!clust) continue;
1201 Int_t sec = clust->GetDetector();
1202 Int_t row = clust->GetRow();
1204 // filter overlapping pad rows needed by HLT
1205 if(sec<fkNIS*2) { //IROCs
1206 if(row == 30) continue;
1209 if(row == 27 || row == 76) continue;
1215 fInnerSec[sec%fkNIS].InsertCluster(clust, count[sec][row], fkParam);
1218 left = (sec-fkNIS*2)/fkNOS;
1219 fOuterSec[(sec-fkNIS*2)%fkNOS].InsertCluster(clust, count[sec][row], fkParam);
1223 // Load functions must be called behind LoadCluster(TClonesArray*)
1225 //LoadOuterSectors();
1226 //LoadInnerSectors();
1232 Int_t AliTPCtrackerMI::LoadClusters()
1235 // load clusters to the memory
1236 AliTPCClustersRow *clrow= new AliTPCClustersRow;
1237 clrow->SetClass("AliTPCclusterMI");
1239 clrow->GetArray()->ExpandCreateFast(10000);
1241 // TTree * tree = fClustersArray.GetTree();
1243 TTree * tree = fInput;
1244 TBranch * br = tree->GetBranch("Segment");
1245 br->SetAddress(&clrow);
1247 Int_t j=Int_t(tree->GetEntries());
1248 for (Int_t i=0; i<j; i++) {
1252 fkParam->AdjustSectorRow(clrow->GetID(),sec,row);
1253 for (Int_t icl=0; icl<clrow->GetArray()->GetEntriesFast(); icl++){
1254 Transform((AliTPCclusterMI*)(clrow->GetArray()->At(icl)));
1257 AliTPCtrackerRow * tpcrow=0;
1260 tpcrow = &(fInnerSec[sec%fkNIS][row]);
1264 tpcrow = &(fOuterSec[(sec-fkNIS*2)%fkNOS][row]);
1265 left = (sec-fkNIS*2)/fkNOS;
1268 tpcrow->SetN1(clrow->GetArray()->GetEntriesFast());
1269 tpcrow->SetClusters1(new AliTPCclusterMI[tpcrow->GetN1()]);
1270 for (Int_t k=0;k<tpcrow->GetN1();++k)
1271 tpcrow->SetCluster1(k, *(AliTPCclusterMI*)(clrow->GetArray()->At(k)));
1274 tpcrow->SetN2(clrow->GetArray()->GetEntriesFast());
1275 tpcrow->SetClusters2(new AliTPCclusterMI[tpcrow->GetN2()]);
1276 for (Int_t k=0;k<tpcrow->GetN2();++k)
1277 tpcrow->SetCluster2(k,*(AliTPCclusterMI*)(clrow->GetArray()->At(k)));
1288 void AliTPCtrackerMI::UnloadClusters()
1291 // unload clusters from the memory
1293 Int_t nrows = fOuterSec->GetNRows();
1294 for (Int_t sec = 0;sec<fkNOS;sec++)
1295 for (Int_t row = 0;row<nrows;row++){
1296 AliTPCtrackerRow* tpcrow = &(fOuterSec[sec%fkNOS][row]);
1298 // if (tpcrow->fClusters1) delete []tpcrow->fClusters1;
1299 // if (tpcrow->fClusters2) delete []tpcrow->fClusters2;
1301 tpcrow->ResetClusters();
1304 nrows = fInnerSec->GetNRows();
1305 for (Int_t sec = 0;sec<fkNIS;sec++)
1306 for (Int_t row = 0;row<nrows;row++){
1307 AliTPCtrackerRow* tpcrow = &(fInnerSec[sec%fkNIS][row]);
1309 // if (tpcrow->fClusters1) delete []tpcrow->fClusters1;
1310 //if (tpcrow->fClusters2) delete []tpcrow->fClusters2;
1312 tpcrow->ResetClusters();
1318 void AliTPCtrackerMI::FillClusterArray(TObjArray* array) const{
1320 // Filling cluster to the array - For visualization purposes
1323 nrows = fOuterSec->GetNRows();
1324 for (Int_t sec = 0;sec<fkNOS;sec++)
1325 for (Int_t row = 0;row<nrows;row++){
1326 AliTPCtrackerRow* tpcrow = &(fOuterSec[sec%fkNOS][row]);
1327 if (!tpcrow) continue;
1328 for (Int_t icl = 0;icl<tpcrow->GetN();icl++){
1329 array->AddLast((TObject*)((*tpcrow)[icl]));
1332 nrows = fInnerSec->GetNRows();
1333 for (Int_t sec = 0;sec<fkNIS;sec++)
1334 for (Int_t row = 0;row<nrows;row++){
1335 AliTPCtrackerRow* tpcrow = &(fInnerSec[sec%fkNIS][row]);
1336 if (!tpcrow) continue;
1337 for (Int_t icl = 0;icl<tpcrow->GetN();icl++){
1338 array->AddLast((TObject*)(*tpcrow)[icl]);
1344 void AliTPCtrackerMI::Transform(AliTPCclusterMI * cluster){
1348 AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ;
1350 AliFatal("Tranformations not in calibDB");
1352 transform->SetCurrentRecoParam((AliTPCRecoParam*)AliTPCReconstructor::GetRecoParam());
1353 Double_t x[3]={cluster->GetRow(),cluster->GetPad(),cluster->GetTimeBin()};
1354 Int_t i[1]={cluster->GetDetector()};
1355 transform->Transform(x,i,0,1);
1356 // if (cluster->GetDetector()%36>17){
1361 // in debug mode check the transformation
1363 if (AliTPCReconstructor::StreamLevel()>2) {
1365 cluster->GetGlobalXYZ(gx);
1366 Int_t event = (fEvent==NULL)? 0: fEvent->GetEventNumberInFile();
1367 TTreeSRedirector &cstream = *fDebugStreamer;
1368 cstream<<"Transform"<<
1379 cluster->SetX(x[0]);
1380 cluster->SetY(x[1]);
1381 cluster->SetZ(x[2]);
1387 //if (!fkParam->IsGeoRead()) fkParam->ReadGeoMatrices();
1388 TGeoHMatrix *mat = fkParam->GetClusterMatrix(cluster->GetDetector());
1390 Double_t pos[3]= {cluster->GetX(),cluster->GetY(),cluster->GetZ()};
1391 Double_t posC[3]={cluster->GetX(),cluster->GetY(),cluster->GetZ()};
1392 if (mat) mat->LocalToMaster(pos,posC);
1394 // chack Loading of Geo matrices from GeoManager - TEMPORARY FIX
1396 cluster->SetX(posC[0]);
1397 cluster->SetY(posC[1]);
1398 cluster->SetZ(posC[2]);
1401 //_____________________________________________________________________________
1402 Int_t AliTPCtrackerMI::LoadOuterSectors() {
1403 //-----------------------------------------------------------------
1404 // This function fills outer TPC sectors with clusters.
1405 //-----------------------------------------------------------------
1406 Int_t nrows = fOuterSec->GetNRows();
1408 for (Int_t sec = 0;sec<fkNOS;sec++)
1409 for (Int_t row = 0;row<nrows;row++){
1410 AliTPCtrackerRow* tpcrow = &(fOuterSec[sec%fkNOS][row]);
1411 Int_t sec2 = sec+2*fkNIS;
1413 Int_t ncl = tpcrow->GetN1();
1415 AliTPCclusterMI *c= (tpcrow->GetCluster1(ncl));
1416 index=(((sec2<<8)+row)<<16)+ncl;
1417 tpcrow->InsertCluster(c,index);
1420 ncl = tpcrow->GetN2();
1422 AliTPCclusterMI *c= (tpcrow->GetCluster2(ncl));
1423 index=((((sec2+fkNOS)<<8)+row)<<16)+ncl;
1424 tpcrow->InsertCluster(c,index);
1427 // write indexes for fast acces
1429 for (Int_t i=0;i<510;i++)
1430 tpcrow->SetFastCluster(i,-1);
1431 for (Int_t i=0;i<tpcrow->GetN();i++){
1432 Int_t zi = Int_t((*tpcrow)[i]->GetZ()+255.);
1433 tpcrow->SetFastCluster(zi,i); // write index
1436 for (Int_t i=0;i<510;i++){
1437 if (tpcrow->GetFastCluster(i)<0)
1438 tpcrow->SetFastCluster(i,last);
1440 last = tpcrow->GetFastCluster(i);
1449 //_____________________________________________________________________________
1450 Int_t AliTPCtrackerMI::LoadInnerSectors() {
1451 //-----------------------------------------------------------------
1452 // This function fills inner TPC sectors with clusters.
1453 //-----------------------------------------------------------------
1454 Int_t nrows = fInnerSec->GetNRows();
1456 for (Int_t sec = 0;sec<fkNIS;sec++)
1457 for (Int_t row = 0;row<nrows;row++){
1458 AliTPCtrackerRow* tpcrow = &(fInnerSec[sec%fkNIS][row]);
1461 Int_t ncl = tpcrow->GetN1();
1463 AliTPCclusterMI *c= (tpcrow->GetCluster1(ncl));
1464 index=(((sec<<8)+row)<<16)+ncl;
1465 tpcrow->InsertCluster(c,index);
1468 ncl = tpcrow->GetN2();
1470 AliTPCclusterMI *c= (tpcrow->GetCluster2(ncl));
1471 index=((((sec+fkNIS)<<8)+row)<<16)+ncl;
1472 tpcrow->InsertCluster(c,index);
1475 // write indexes for fast acces
1477 for (Int_t i=0;i<510;i++)
1478 tpcrow->SetFastCluster(i,-1);
1479 for (Int_t i=0;i<tpcrow->GetN();i++){
1480 Int_t zi = Int_t((*tpcrow)[i]->GetZ()+255.);
1481 tpcrow->SetFastCluster(zi,i); // write index
1484 for (Int_t i=0;i<510;i++){
1485 if (tpcrow->GetFastCluster(i)<0)
1486 tpcrow->SetFastCluster(i,last);
1488 last = tpcrow->GetFastCluster(i);
1500 //_________________________________________________________________________
1501 AliTPCclusterMI *AliTPCtrackerMI::GetClusterMI(Int_t index) const {
1502 //--------------------------------------------------------------------
1503 // Return pointer to a given cluster
1504 //--------------------------------------------------------------------
1505 if (index<0) return 0; // no cluster
1506 Int_t sec=(index&0xff000000)>>24;
1507 Int_t row=(index&0x00ff0000)>>16;
1508 Int_t ncl=(index&0x00007fff)>>00;
1510 const AliTPCtrackerRow * tpcrow=0;
1511 AliTPCclusterMI * clrow =0;
1513 if (sec<0 || sec>=fkNIS*4) {
1514 AliWarning(Form("Wrong sector %d",sec));
1519 AliTPCtrackerSector& tracksec = fInnerSec[sec%fkNIS];
1520 if (tracksec.GetNRows()<=row) return 0;
1521 tpcrow = &(tracksec[row]);
1522 if (tpcrow==0) return 0;
1525 if (tpcrow->GetN1()<=ncl) return 0;
1526 clrow = tpcrow->GetClusters1();
1529 if (tpcrow->GetN2()<=ncl) return 0;
1530 clrow = tpcrow->GetClusters2();
1534 AliTPCtrackerSector& tracksec = fOuterSec[(sec-fkNIS*2)%fkNOS];
1535 if (tracksec.GetNRows()<=row) return 0;
1536 tpcrow = &(tracksec[row]);
1537 if (tpcrow==0) return 0;
1539 if (sec-2*fkNIS<fkNOS) {
1540 if (tpcrow->GetN1()<=ncl) return 0;
1541 clrow = tpcrow->GetClusters1();
1544 if (tpcrow->GetN2()<=ncl) return 0;
1545 clrow = tpcrow->GetClusters2();
1549 return &(clrow[ncl]);
1555 Int_t AliTPCtrackerMI::FollowToNext(AliTPCseed& t, Int_t nr) {
1556 //-----------------------------------------------------------------
1557 // This function tries to find a track prolongation to next pad row
1558 //-----------------------------------------------------------------
1560 Double_t x= GetXrow(nr), ymax=GetMaxY(nr);
1563 AliTPCclusterMI *cl=0;
1564 Int_t tpcindex= t.GetClusterIndex2(nr);
1566 // update current shape info every 5 pad-row
1567 // if ( (nr%5==0) || t.GetNumberOfClusters()<2 || (t.fCurrentSigmaY2<0.0001) ){
1571 if (fIteration>0 && tpcindex>=-1){ //if we have already clusters
1573 if (tpcindex==-1) return 0; //track in dead zone
1575 cl = t.GetClusterPointer(nr);
1576 if ( (cl==0) ) cl = GetClusterMI(tpcindex);
1577 t.SetCurrentClusterIndex1(tpcindex);
1580 Int_t relativesector = ((tpcindex&0xff000000)>>24)%18; // if previously accepted cluster in different sector
1581 Float_t angle = relativesector*fSectors->GetAlpha()+fSectors->GetAlphaShift();
1583 if (angle<-TMath::Pi()) angle += 2*TMath::Pi();
1584 if (angle>=TMath::Pi()) angle -= 2*TMath::Pi();
1586 if (TMath::Abs(angle-t.GetAlpha())>0.001){
1587 Double_t rotation = angle-t.GetAlpha();
1588 t.SetRelativeSector(relativesector);
1589 if (!t.Rotate(rotation)) return 0;
1591 if (!t.PropagateTo(x)) return 0;
1593 t.SetCurrentCluster(cl);
1595 Int_t accept = AcceptCluster(&t,t.GetCurrentCluster());
1596 if ((tpcindex&0x8000)==0) accept =0;
1598 //if founded cluster is acceptible
1599 if (cl->IsUsed(11)) { // id cluster is shared inrease uncertainty
1600 t.SetErrorY2(t.GetErrorY2()+0.03);
1601 t.SetErrorZ2(t.GetErrorZ2()+0.03);
1602 t.SetErrorY2(t.GetErrorY2()*3);
1603 t.SetErrorZ2(t.GetErrorZ2()*3);
1605 t.SetNFoundable(t.GetNFoundable()+1);
1606 UpdateTrack(&t,accept);
1611 if (TMath::Abs(t.GetSnp())>AliTPCReconstructor::GetMaxSnpTracker()) return 0; // cut on angle
1612 if (fIteration>1 && IsFindable(t)){
1613 // not look for new cluster during refitting
1614 t.SetNFoundable(t.GetNFoundable()+1);
1619 // if (TMath::Abs(t.GetSnp())>0.95 || TMath::Abs(x*t.GetC()-t.GetEta())>0.95) return 0;// patch 28 fev 06
1620 Double_t y=t.GetYat(x);
1621 if (TMath::Abs(y)>ymax){
1623 t.SetRelativeSector((t.GetRelativeSector()+1) % fN);
1624 if (!t.Rotate(fSectors->GetAlpha()))
1626 } else if (y <-ymax) {
1627 t.SetRelativeSector((t.GetRelativeSector()-1+fN) % fN);
1628 if (!t.Rotate(-fSectors->GetAlpha()))
1634 if (!t.PropagateTo(x)) {
1635 if (fIteration==0) t.SetRemoval(10);
1639 Double_t z=t.GetZ();
1642 if (!IsActive(t.GetRelativeSector(),nr)) {
1644 t.SetClusterIndex2(nr,-1);
1647 //AliInfo(Form("A - Sector%d phi %f - alpha %f", t.fRelativeSector,y/x, t.GetAlpha()));
1648 Bool_t isActive = IsActive(t.GetRelativeSector(),nr);
1649 Bool_t isActive2 = (nr>=fInnerSec->GetNRows()) ? fOuterSec[t.GetRelativeSector()][nr-fInnerSec->GetNRows()].GetN()>0:fInnerSec[t.GetRelativeSector()][nr].GetN()>0;
1651 if (!isActive || !isActive2) return 0;
1653 const AliTPCtrackerRow &krow=GetRow(t.GetRelativeSector(),nr);
1654 if ( (t.GetSigmaY2()<0) || t.GetSigmaZ2()<0) return 0;
1656 Double_t roadz = 1.;
1658 if (TMath::Abs(TMath::Abs(y)-ymax)<krow.GetDeadZone()){
1660 t.SetClusterIndex2(nr,-1);
1666 // if (TMath::Abs(z)<(AliTPCReconstructor::GetCtgRange()*x+10) && TMath::Abs(z)<fkParam->GetZLength(0) && (TMath::Abs(t.GetSnp())<AliTPCReconstructor::GetMaxSnpTracker()))
1667 t.SetNFoundable(t.GetNFoundable()+1);
1673 // cl = krow.FindNearest2(y+10.,z,roady,roadz,index);
1674 cl = krow.FindNearest2(y,z,roady,roadz,index);
1675 if (cl) t.SetCurrentClusterIndex1(krow.GetIndex(index));
1678 t.SetCurrentCluster(cl);
1680 if (fIteration==2&&cl->IsUsed(10)) return 0;
1681 Int_t accept = AcceptCluster(&t,t.GetCurrentCluster());
1682 if (fIteration==2&&cl->IsUsed(11)) {
1683 t.SetErrorY2(t.GetErrorY2()+0.03);
1684 t.SetErrorZ2(t.GetErrorZ2()+0.03);
1685 t.SetErrorY2(t.GetErrorY2()*3);
1686 t.SetErrorZ2(t.GetErrorZ2()*3);
1689 if (t.fCurrentCluster->IsUsed(10)){
1694 if (t.fNShared>0.7*t.GetNumberOfClusters()) {
1700 if (accept<3) UpdateTrack(&t,accept);
1703 if ( fIteration==0 && t.GetNFoundable()*0.5 > t.GetNumberOfClusters()) t.SetRemoval(10);
1711 //_________________________________________________________________________
1712 Bool_t AliTPCtrackerMI::GetTrackPoint(Int_t index, AliTrackPoint &p ) const
1714 // Get track space point by index
1715 // return false in case the cluster doesn't exist
1716 AliTPCclusterMI *cl = GetClusterMI(index);
1717 if (!cl) return kFALSE;
1718 Int_t sector = (index&0xff000000)>>24;
1719 // Int_t row = (index&0x00ff0000)>>16;
1721 // xyz[0] = fkParam->GetPadRowRadii(sector,row);
1722 xyz[0] = cl->GetX();
1723 xyz[1] = cl->GetY();
1724 xyz[2] = cl->GetZ();
1726 fkParam->AdjustCosSin(sector,cos,sin);
1727 Float_t x = cos*xyz[0]-sin*xyz[1];
1728 Float_t y = cos*xyz[1]+sin*xyz[0];
1730 Float_t sigmaY2 = 0.027*cl->GetSigmaY2();
1731 if (sector < fkParam->GetNInnerSector()) sigmaY2 *= 2.07;
1732 Float_t sigmaZ2 = 0.066*cl->GetSigmaZ2();
1733 if (sector < fkParam->GetNInnerSector()) sigmaZ2 *= 1.77;
1734 cov[0] = sin*sin*sigmaY2;
1735 cov[1] = -sin*cos*sigmaY2;
1737 cov[3] = cos*cos*sigmaY2;
1740 p.SetXYZ(x,y,xyz[2],cov);
1741 AliGeomManager::ELayerID iLayer;
1743 if (sector < fkParam->GetNInnerSector()) {
1744 iLayer = AliGeomManager::kTPC1;
1748 iLayer = AliGeomManager::kTPC2;
1749 idet = sector - fkParam->GetNInnerSector();
1751 UShort_t volid = AliGeomManager::LayerToVolUID(iLayer,idet);
1752 p.SetVolumeID(volid);
1758 Int_t AliTPCtrackerMI::UpdateClusters(AliTPCseed& t, Int_t nr) {
1759 //-----------------------------------------------------------------
1760 // This function tries to find a track prolongation to next pad row
1761 //-----------------------------------------------------------------
1762 t.SetCurrentCluster(0);
1763 t.SetCurrentClusterIndex1(0);
1765 Double_t xt=t.GetX();
1766 Int_t row = GetRowNumber(xt)-1;
1767 Double_t ymax= GetMaxY(nr);
1769 if (row < nr) return 1; // don't prolongate if not information until now -
1770 // if (TMath::Abs(t.GetSnp())>0.9 && t.GetNumberOfClusters()>40. && fIteration!=2) {
1772 // return 0; // not prolongate strongly inclined tracks
1774 // if (TMath::Abs(t.GetSnp())>0.95) {
1776 // return 0; // not prolongate strongly inclined tracks
1777 // }// patch 28 fev 06
1779 Double_t x= GetXrow(nr);
1781 //t.PropagateTo(x+0.02);
1782 //t.PropagateTo(x+0.01);
1783 if (!t.PropagateTo(x)){
1790 if (TMath::Abs(y)>ymax){
1792 t.SetRelativeSector((t.GetRelativeSector()+1) % fN);
1793 if (!t.Rotate(fSectors->GetAlpha()))
1795 } else if (y <-ymax) {
1796 t.SetRelativeSector((t.GetRelativeSector()-1+fN) % fN);
1797 if (!t.Rotate(-fSectors->GetAlpha()))
1800 // if (!t.PropagateTo(x)){
1807 if (TMath::Abs(t.GetSnp())>AliTPCReconstructor::GetMaxSnpTracker()) return 0;
1809 if (!IsActive(t.GetRelativeSector(),nr)) {
1811 t.SetClusterIndex2(nr,-1);
1814 //AliInfo(Form("A - Sector%d phi %f - alpha %f", t.fRelativeSector,y/x, t.GetAlpha()));
1816 AliTPCtrackerRow &krow=GetRow(t.GetRelativeSector(),nr);
1818 if (TMath::Abs(TMath::Abs(y)-ymax)<krow.GetDeadZone()){
1820 t.SetClusterIndex2(nr,-1);
1826 // if (TMath::Abs(t.GetZ())<(AliTPCReconstructor::GetCtgRange()*t.GetX()+10) && (TMath::Abs(t.GetSnp())<AliTPCReconstructor::GetMaxSnpTracker()))
1827 if (IsFindable(t)) t.SetNFoundable(t.GetNFoundable()+1);
1833 if ( (nr%2==0) || t.GetNumberOfClusters()<2){
1834 // t.fCurrentSigmaY = GetSigmaY(&t);
1835 //t.fCurrentSigmaZ = GetSigmaZ(&t);
1839 AliTPCclusterMI *cl=0;
1842 Double_t roady = 1.;
1843 Double_t roadz = 1.;
1847 index = t.GetClusterIndex2(nr);
1848 if ( (index>0) && (index&0x8000)==0){
1849 cl = t.GetClusterPointer(nr);
1850 if ( (cl==0) && (index>0)) cl = GetClusterMI(index);
1851 t.SetCurrentClusterIndex1(index);
1853 t.SetCurrentCluster(cl);
1859 // if (index<0) return 0;
1860 UInt_t uindex = TMath::Abs(index);
1863 //cl = krow.FindNearest2(y+10,z,roady,roadz,uindex);
1864 cl = krow.FindNearest2(y,z,roady,roadz,uindex);
1867 if (cl) t.SetCurrentClusterIndex1(krow.GetIndex(uindex));
1868 t.SetCurrentCluster(cl);
1874 Int_t AliTPCtrackerMI::FollowToNextCluster(AliTPCseed & t, Int_t nr) {
1875 //-----------------------------------------------------------------
1876 // This function tries to find a track prolongation to next pad row
1877 //-----------------------------------------------------------------
1879 //update error according neighborhoud
1881 if (t.GetCurrentCluster()) {
1883 Int_t accept = AcceptCluster(&t,t.GetCurrentCluster());
1885 if (t.GetCurrentCluster()->IsUsed(10)){
1890 t.SetNShared(t.GetNShared()+1);
1891 if (t.GetNShared()>0.7*t.GetNumberOfClusters()) {
1896 if (fIteration>0) accept = 0;
1897 if (accept<3) UpdateTrack(&t,accept);
1901 if ( ( (t.GetSigmaY2()+t.GetSigmaZ2())>0.16)&& t.GetNumberOfClusters()>18) t.SetRemoval(10);
1902 if ( t.GetChi2()/t.GetNumberOfClusters()>6 &&t.GetNumberOfClusters()>18) t.SetRemoval(10);
1904 if (( (t.GetNFoundable()*0.5 > t.GetNumberOfClusters()) || t.GetNoCluster()>15)) t.SetRemoval(10);
1912 //_____________________________________________________________________________
1913 Int_t AliTPCtrackerMI::FollowProlongation(AliTPCseed& t, Int_t rf, Int_t step) {
1914 //-----------------------------------------------------------------
1915 // This function tries to find a track prolongation.
1916 //-----------------------------------------------------------------
1917 Double_t xt=t.GetX();
1919 Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
1920 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
1921 if (alpha < 0. ) alpha += 2.*TMath::Pi();
1923 t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN);
1925 Int_t first = GetRowNumber(xt)-1;
1926 for (Int_t nr= first; nr>=rf; nr-=step) {
1928 if (t.GetKinkIndexes()[0]>0){
1929 for (Int_t i=0;i<3;i++){
1930 Int_t index = t.GetKinkIndexes()[i];
1931 if (index==0) break;
1932 if (index<0) continue;
1934 AliKink * kink = (AliKink*)fEvent->GetKink(index-1);
1936 printf("PROBLEM\n");
1939 Int_t kinkrow = kink->GetTPCRow0()+2+Int_t(0.5/(0.05+kink->GetAngle(2)));
1941 AliExternalTrackParam paramd(t);
1942 kink->SetDaughter(paramd);
1943 kink->SetStatus(2,5);
1950 if (nr==80) t.UpdateReference();
1951 if (nr<fInnerSec->GetNRows())
1952 fSectors = fInnerSec;
1954 fSectors = fOuterSec;
1955 if (FollowToNext(t,nr)==0)
1968 Int_t AliTPCtrackerMI::FollowBackProlongation(AliTPCseed& t, Int_t rf) {
1969 //-----------------------------------------------------------------
1970 // This function tries to find a track prolongation.
1971 //-----------------------------------------------------------------
1973 Double_t xt=t.GetX();
1974 Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
1975 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
1976 if (alpha < 0. ) alpha += 2.*TMath::Pi();
1977 t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN);
1979 Int_t first = t.GetFirstPoint();
1980 if (first<GetRowNumber(xt)+1) first = GetRowNumber(xt)+1;
1982 if (first<0) first=0;
1983 for (Int_t nr=first; nr<=rf; nr++) {
1984 // if ( (TMath::Abs(t.GetSnp())>0.95)) break;//patch 28 fev 06
1985 if (t.GetKinkIndexes()[0]<0){
1986 for (Int_t i=0;i<3;i++){
1987 Int_t index = t.GetKinkIndexes()[i];
1988 if (index==0) break;
1989 if (index>0) continue;
1990 index = TMath::Abs(index);
1991 AliKink * kink = (AliKink*)fEvent->GetKink(index-1);
1993 printf("PROBLEM\n");
1996 Int_t kinkrow = kink->GetTPCRow0()-2-Int_t(0.5/(0.05+kink->GetAngle(2)));
1998 AliExternalTrackParam paramm(t);
1999 kink->SetMother(paramm);
2000 kink->SetStatus(2,1);
2007 if (nr<fInnerSec->GetNRows())
2008 fSectors = fInnerSec;
2010 fSectors = fOuterSec;
2021 Float_t AliTPCtrackerMI::OverlapFactor(AliTPCseed * s1, AliTPCseed * s2, Int_t &sum1, Int_t & sum2)
2029 Float_t dz2 =(s1->GetZ() - s2->GetZ());
2032 Float_t dy2 =TMath::Abs((s1->GetY() - s2->GetY()));
2034 Float_t distance = TMath::Sqrt(dz2+dy2);
2035 if (distance>4.) return 0; // if there are far away - not overlap - to reduce combinatorics
2038 Int_t firstpoint = TMath::Min(s1->GetFirstPoint(),s2->GetFirstPoint());
2039 Int_t lastpoint = TMath::Max(s1->GetLastPoint(),s2->GetLastPoint());
2044 if (firstpoint>lastpoint) {
2045 firstpoint =lastpoint;
2050 for (Int_t i=firstpoint-1;i<lastpoint+1;i++){
2051 if (s1->GetClusterIndex2(i)>0) sum1++;
2052 if (s2->GetClusterIndex2(i)>0) sum2++;
2053 if (s1->GetClusterIndex2(i)==s2->GetClusterIndex2(i) && s1->GetClusterIndex2(i)>0) {
2057 if (sum<5) return 0;
2059 Float_t summin = TMath::Min(sum1+1,sum2+1);
2060 Float_t ratio = (sum+1)/Float_t(summin);
2064 void AliTPCtrackerMI::SignShared(AliTPCseed * s1, AliTPCseed * s2)
2068 Float_t thetaCut = 0.2;//+10.*TMath::Sqrt(s1->GetSigmaTglZ()+ s2->GetSigmaTglZ());
2069 if (TMath::Abs(s1->GetTgl()-s2->GetTgl())>thetaCut) return;
2070 Float_t minCl = TMath::Min(s1->GetNumberOfClusters(),s2->GetNumberOfClusters());
2071 Int_t cutN0 = TMath::Max(5,TMath::Nint(0.1*minCl));
2076 //Int_t firstpoint = TMath::Max(s1->GetFirstPoint(),s2->GetFirstPoint());
2077 //Int_t lastpoint = TMath::Min(s1->GetLastPoint(),s2->GetLastPoint());
2078 Int_t firstpoint = 0;
2079 Int_t lastpoint = 160;
2081 // if (firstpoint>=lastpoint-5) return;;
2083 for (Int_t i=firstpoint;i<lastpoint;i++){
2084 // if ( (s1->GetClusterIndex2(i)&0xFFFF8FFF)==(s2->GetClusterIndex2(i)&0xFFFF8FFF) && s1->GetClusterIndex2(i)>0) {
2085 if ( (s1->GetClusterIndex2(i))==(s2->GetClusterIndex2(i)) && s1->GetClusterIndex2(i)>0) {
2087 s1->SetSharedMapBit(i, kTRUE);
2088 s2->SetSharedMapBit(i, kTRUE);
2090 if (s1->GetClusterIndex2(i)>0)
2091 s1->SetClusterMapBit(i, kTRUE);
2093 if (sumshared>cutN0){
2096 for (Int_t i=firstpoint;i<lastpoint;i++){
2097 // if ( (s1->GetClusterIndex2(i)&0xFFFF8FFF)==(s2->GetClusterIndex2(i)&0xFFFF8FFF) && s1->GetClusterIndex2(i)>0) {
2098 if ( (s1->GetClusterIndex2(i))==(s2->GetClusterIndex2(i)) && s1->GetClusterIndex2(i)>0) {
2099 AliTPCTrackerPoint *p1 = s1->GetTrackPoint(i);
2100 AliTPCTrackerPoint *p2 = s2->GetTrackPoint(i);;
2101 if (s1->IsActive()&&s2->IsActive()){
2102 p1->SetShared(kTRUE);
2103 p2->SetShared(kTRUE);
2109 if (sumshared>cutN0){
2110 for (Int_t i=0;i<4;i++){
2111 if (s1->GetOverlapLabel(3*i)==0){
2112 s1->SetOverlapLabel(3*i, s2->GetLabel());
2113 s1->SetOverlapLabel(3*i+1,sumshared);
2114 s1->SetOverlapLabel(3*i+2,s2->GetUniqueID());
2118 for (Int_t i=0;i<4;i++){
2119 if (s2->GetOverlapLabel(3*i)==0){
2120 s2->SetOverlapLabel(3*i, s1->GetLabel());
2121 s2->SetOverlapLabel(3*i+1,sumshared);
2122 s2->SetOverlapLabel(3*i+2,s1->GetUniqueID());
2129 void AliTPCtrackerMI::SignShared(TObjArray * arr)
2132 //sort trackss according sectors
2134 for (Int_t i=0; i<arr->GetEntriesFast(); i++) {
2135 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2137 //if (pt) RotateToLocal(pt);
2141 arr->Sort(); // sorting according relative sectors
2142 arr->Expand(arr->GetEntries());
2145 Int_t nseed=arr->GetEntriesFast();
2146 for (Int_t i=0; i<nseed; i++) {
2147 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2149 for (Int_t j=0;j<=12;j++){
2150 pt->SetOverlapLabel(j,0);
2153 for (Int_t i=0; i<nseed; i++) {
2154 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2156 if (pt->GetRemoval()>10) continue;
2157 for (Int_t j=i+1; j<nseed; j++){
2158 AliTPCseed *pt2=(AliTPCseed*)arr->UncheckedAt(j);
2159 if (TMath::Abs(pt->GetRelativeSector()-pt2->GetRelativeSector())>1) continue;
2161 if (pt2->GetRemoval()<=10) {
2162 //if ( TMath::Abs(pt->GetRelativeSector()-pt2->GetRelativeSector())>0) break;
2170 void AliTPCtrackerMI::SortTracks(TObjArray * arr, Int_t mode) const
2173 //sort tracks in array according mode criteria
2174 Int_t nseed = arr->GetEntriesFast();
2175 for (Int_t i=0; i<nseed; i++) {
2176 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2187 void AliTPCtrackerMI::RemoveUsed2(TObjArray * arr, Float_t factor1, Float_t factor2, Int_t minimal)
2190 // Loop over all tracks and remove overlaped tracks (with lower quality)
2192 // 1. Unsign clusters
2193 // 2. Sort tracks according quality
2194 // Quality is defined by the number of cluster between first and last points
2196 // 3. Loop over tracks - decreasing quality order
2197 // a.) remove - If the fraction of shared cluster less than factor (1- n or 2)
2198 // b.) remove - If the minimal number of clusters less than minimal and not ITS
2199 // c.) if track accepted - sign clusters
2201 //Called in - AliTPCtrackerMI::Clusters2Tracks()
2202 // - AliTPCtrackerMI::PropagateBack()
2203 // - AliTPCtrackerMI::RefitInward()
2206 // factor1 - factor for constrained
2207 // factor2 - for non constrained tracks
2208 // if (Float_t(shared+1)/Float_t(found+1)>factor) - DELETE
2212 Int_t nseed = arr->GetEntriesFast();
2213 Float_t * quality = new Float_t[nseed];
2214 Int_t * indexes = new Int_t[nseed];
2218 for (Int_t i=0; i<nseed; i++) {
2219 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2224 pt->UpdatePoints(); //select first last max dens points
2225 Float_t * points = pt->GetPoints();
2226 if (points[3]<0.8) quality[i] =-1;
2227 quality[i] = (points[2]-points[0])+pt->GetNumberOfClusters();
2228 //prefer high momenta tracks if overlaps
2229 quality[i] *= TMath::Sqrt(TMath::Abs(pt->Pt())+0.5);
2231 TMath::Sort(nseed,quality,indexes);
2234 for (Int_t itrack=0; itrack<nseed; itrack++) {
2235 Int_t trackindex = indexes[itrack];
2236 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(trackindex);
2239 if (quality[trackindex]<0){
2241 delete arr->RemoveAt(trackindex);
2244 arr->RemoveAt(trackindex);
2250 Int_t first = Int_t(pt->GetPoints()[0]);
2251 Int_t last = Int_t(pt->GetPoints()[2]);
2252 Double_t factor = (pt->GetBConstrain()) ? factor1: factor2;
2254 Int_t found,foundable,shared;
2255 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
2256 // pt->GetClusterStatistic(0,160, found, foundable,shared,kFALSE);
2257 Bool_t itsgold =kFALSE;
2260 if (pt->GetESD()->GetITSclusters(dummy)>4) itsgold= kTRUE;
2264 if (Float_t(shared+1)/Float_t(found+1)>factor){
2265 if (pt->GetKinkIndexes()[0]!=0) continue; //don't remove tracks - part of the kinks
2266 if( AliTPCReconstructor::StreamLevel()>3){
2267 TTreeSRedirector &cstream = *fDebugStreamer;
2268 cstream<<"RemoveUsed"<<
2269 "iter="<<fIteration<<
2273 delete arr->RemoveAt(trackindex);
2276 if (pt->GetNumberOfClusters()<50&&(found-0.5*shared)<minimal){ //remove short tracks
2277 if (pt->GetKinkIndexes()[0]!=0) continue; //don't remove tracks - part of the kinks
2278 if( AliTPCReconstructor::StreamLevel()>3){
2279 TTreeSRedirector &cstream = *fDebugStreamer;
2280 cstream<<"RemoveShort"<<
2281 "iter="<<fIteration<<
2285 delete arr->RemoveAt(trackindex);
2291 //if (sharedfactor>0.4) continue;
2292 if (pt->GetKinkIndexes()[0]>0) continue;
2293 //Remove tracks with undefined properties - seems
2294 if (pt->GetSigmaY2()<kAlmost0) continue; // ? what is the origin ?
2296 for (Int_t i=first; i<last; i++) {
2297 Int_t index=pt->GetClusterIndex2(i);
2298 // if (index<0 || index&0x8000 ) continue;
2299 if (index<0 || index&0x8000 ) continue;
2300 AliTPCclusterMI *c= pt->GetClusterPointer(i);
2307 Info("RemoveUsed","\n*****\nNumber of good tracks after shared removal\t%d\n",fNtracks);
2313 void AliTPCtrackerMI::UnsignClusters()
2316 // loop over all clusters and unsign them
2319 for (Int_t sec=0;sec<fkNIS;sec++){
2320 for (Int_t row=0;row<fInnerSec->GetNRows();row++){
2321 AliTPCclusterMI *cl = fInnerSec[sec][row].GetClusters1();
2322 for (Int_t icl =0;icl< fInnerSec[sec][row].GetN1();icl++)
2323 // if (cl[icl].IsUsed(10))
2325 cl = fInnerSec[sec][row].GetClusters2();
2326 for (Int_t icl =0;icl< fInnerSec[sec][row].GetN2();icl++)
2327 //if (cl[icl].IsUsed(10))
2332 for (Int_t sec=0;sec<fkNOS;sec++){
2333 for (Int_t row=0;row<fOuterSec->GetNRows();row++){
2334 AliTPCclusterMI *cl = fOuterSec[sec][row].GetClusters1();
2335 for (Int_t icl =0;icl< fOuterSec[sec][row].GetN1();icl++)
2336 //if (cl[icl].IsUsed(10))
2338 cl = fOuterSec[sec][row].GetClusters2();
2339 for (Int_t icl =0;icl< fOuterSec[sec][row].GetN2();icl++)
2340 //if (cl[icl].IsUsed(10))
2349 void AliTPCtrackerMI::SignClusters(const TObjArray * arr, Float_t fnumber, Float_t fdensity)
2352 //sign clusters to be "used"
2354 // snumber and sdensity sign number of sigmas - bellow mean value to be accepted
2355 // loop over "primaries"
2369 Int_t nseed = arr->GetEntriesFast();
2370 for (Int_t i=0; i<nseed; i++) {
2371 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2375 if (!(pt->IsActive())) continue;
2376 Float_t dens = pt->GetNumberOfClusters()/Float_t(pt->GetNFoundable());
2377 if ( (dens>0.7) && (pt->GetNumberOfClusters()>70)){
2379 sumdens2+= dens*dens;
2380 sumn += pt->GetNumberOfClusters();
2381 sumn2 += pt->GetNumberOfClusters()*pt->GetNumberOfClusters();
2382 Float_t chi2 = pt->GetChi2()/pt->GetNumberOfClusters();
2385 sumchi2 +=chi2*chi2;
2390 Float_t mdensity = 0.9;
2391 Float_t meann = 130;
2392 Float_t meanchi = 1;
2393 Float_t sdensity = 0.1;
2394 Float_t smeann = 10;
2395 Float_t smeanchi =0.4;
2399 mdensity = sumdens/sum;
2401 meanchi = sumchi/sum;
2403 sdensity = sumdens2/sum-mdensity*mdensity;
2405 sdensity = TMath::Sqrt(sdensity);
2409 smeann = sumn2/sum-meann*meann;
2411 smeann = TMath::Sqrt(smeann);
2415 smeanchi = sumchi2/sum - meanchi*meanchi;
2417 smeanchi = TMath::Sqrt(smeanchi);
2423 //REMOVE SHORT DELTAS or tracks going out of sensitive volume of TPC
2425 for (Int_t i=0; i<nseed; i++) {
2426 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2430 if (pt->GetBSigned()) continue;
2431 if (pt->GetBConstrain()) continue;
2432 //if (!(pt->IsActive())) continue;
2434 Int_t found,foundable,shared;
2435 pt->GetClusterStatistic(0,160,found, foundable,shared);
2436 if (shared/float(found)>0.3) {
2437 if (shared/float(found)>0.9 ){
2438 //delete arr->RemoveAt(i);
2443 Bool_t isok =kFALSE;
2444 if ( (pt->GetNShared()/pt->GetNumberOfClusters()<0.5) &&pt->GetNumberOfClusters()>60)
2446 if ((TMath::Abs(1/pt->GetC())<100.) && (pt->GetNShared()/pt->GetNumberOfClusters()<0.7))
2448 if (TMath::Abs(pt->GetZ()/pt->GetX())>1.1)
2450 if ( (TMath::Abs(pt->GetSnp()>0.7) && pt->GetD(0,0)>60.))
2454 for (Int_t j=0; j<160; ++j) {
2455 Int_t index=pt->GetClusterIndex2(j);
2456 if (index<0) continue;
2457 AliTPCclusterMI *c= pt->GetClusterPointer(j);
2459 //if (!(c->IsUsed(10))) c->Use();
2466 Double_t maxchi = meanchi+2.*smeanchi;
2468 for (Int_t i=0; i<nseed; i++) {
2469 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2473 //if (!(pt->IsActive())) continue;
2474 if (pt->GetBSigned()) continue;
2475 Double_t chi = pt->GetChi2()/pt->GetNumberOfClusters();
2476 if (chi>maxchi) continue;
2479 Float_t dens = pt->GetNumberOfClusters()/Float_t(pt->GetNFoundable());
2481 //sign only tracks with enoug big density at the beginning
2483 if ((pt->GetDensityFirst(40)<0.75) && pt->GetNumberOfClusters()<meann) continue;
2486 Double_t mindens = TMath::Max(double(mdensity-sdensity*fdensity*bfactor),0.65);
2487 Double_t minn = TMath::Max(Int_t(meann-fnumber*smeann*bfactor),50);
2489 // if (pt->fBConstrain) mindens = TMath::Max(mdensity-sdensity*fdensity*bfactor,0.65);
2490 if ( (pt->GetRemoval()==10) && (pt->GetSnp()>0.8)&&(dens>mindens))
2493 if ((dens>mindens && pt->GetNumberOfClusters()>minn) && chi<maxchi ){
2494 //Int_t noc=pt->GetNumberOfClusters();
2495 pt->SetBSigned(kTRUE);
2496 for (Int_t j=0; j<160; ++j) {
2498 Int_t index=pt->GetClusterIndex2(j);
2499 if (index<0) continue;
2500 AliTPCclusterMI *c= pt->GetClusterPointer(j);
2502 // if (!(c->IsUsed(10))) c->Use();
2507 // gLastCheck = nseed;
2515 void AliTPCtrackerMI::StopNotActive(const TObjArray * arr, Int_t row0, Float_t th0, Float_t th1, Float_t th2) const
2517 // stop not active tracks
2518 // take th1 as threshold for number of founded to number of foundable on last 10 active rows
2519 // take th2 as threshold for number of founded to number of foundable on last 20 active rows
2520 Int_t nseed = arr->GetEntriesFast();
2522 for (Int_t i=0; i<nseed; i++) {
2523 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2527 if (!(pt->IsActive())) continue;
2528 StopNotActive(pt,row0,th0, th1,th2);
2534 void AliTPCtrackerMI::StopNotActive(AliTPCseed * seed, Int_t row0, Float_t th0, Float_t th1,
2537 // stop not active tracks
2538 // take th1 as threshold for number of founded to number of foundable on last 10 active rows
2539 // take th2 as threshold for number of founded to number of foundable on last 20 active rows
2542 Int_t foundable = 0;
2543 Int_t maxindex = seed->GetLastPoint(); //last foundable row
2544 if (seed->GetNFoundable()*th0 > seed->GetNumberOfClusters()) {
2545 seed->Desactivate(10) ;
2549 for (Int_t i=row0; i<maxindex; i++){
2550 Int_t index = seed->GetClusterIndex2(i);
2551 if (index!=-1) foundable++;
2553 if (foundable<=30) sumgood1++;
2554 if (foundable<=50) {
2561 if (foundable>=30.){
2562 if (sumgood1<(th1*30.)) seed->Desactivate(10);
2565 if (sumgood2<(th2*50.)) seed->Desactivate(10);
2569 Int_t AliTPCtrackerMI::RefitInward(AliESDEvent *event)
2572 // back propagation of ESD tracks
2575 const Int_t kMaxFriendTracks=2000;
2579 //PrepareForProlongation(fSeeds,1);
2580 PropagateForward2(fSeeds);
2581 RemoveUsed2(fSeeds,0.4,0.4,20);
2583 TObjArray arraySeed(fSeeds->GetEntries());
2584 for (Int_t i=0;i<fSeeds->GetEntries();i++) {
2585 arraySeed.AddAt(fSeeds->At(i),i);
2587 SignShared(&arraySeed);
2588 // FindCurling(fSeeds, event,2); // find multi found tracks
2589 FindSplitted(fSeeds, event,2); // find multi found tracks
2590 if (AliTPCReconstructor::StreamLevel()>5) FindMultiMC(fSeeds, fEvent,2); // find multi found tracks
2593 Int_t nseed = fSeeds->GetEntriesFast();
2594 for (Int_t i=0;i<nseed;i++){
2595 AliTPCseed * seed = (AliTPCseed*) fSeeds->UncheckedAt(i);
2596 if (!seed) continue;
2597 if (seed->GetKinkIndex(0)>0) UpdateKinkQualityD(seed); // update quality informations for kinks
2598 AliESDtrack *esd=event->GetTrack(i);
2600 if (seed->GetNumberOfClusters()<60 && seed->GetNumberOfClusters()<(esd->GetTPCclusters(0) -5)*0.8){
2601 AliExternalTrackParam paramIn;
2602 AliExternalTrackParam paramOut;
2603 Int_t ncl = seed->RefitTrack(seed,¶mIn,¶mOut);
2604 if (AliTPCReconstructor::StreamLevel()>2) {
2605 (*fDebugStreamer)<<"RecoverIn"<<
2609 "pout.="<<¶mOut<<
2614 seed->Set(paramIn.GetX(),paramIn.GetAlpha(),paramIn.GetParameter(),paramIn.GetCovariance());
2615 seed->SetNumberOfClusters(ncl);
2619 seed->PropagateTo(fkParam->GetInnerRadiusLow());
2620 seed->UpdatePoints();
2621 AddCovariance(seed);
2623 seed->CookdEdx(0.02,0.6);
2624 CookLabel(seed,0.1); //For comparison only
2625 esd->SetTPCClusterMap(seed->GetClusterMap());
2626 esd->SetTPCSharedMap(seed->GetSharedMap());
2628 if (AliTPCReconstructor::StreamLevel()>1 && seed!=0&&esd!=0) {
2629 TTreeSRedirector &cstream = *fDebugStreamer;
2636 if (seed->GetNumberOfClusters()>15){
2637 esd->UpdateTrackParams(seed,AliESDtrack::kTPCrefit);
2638 esd->SetTPCPoints(seed->GetPoints());
2639 esd->SetTPCPointsF(seed->GetNFoundable());
2640 Int_t ndedx = seed->GetNCDEDX(0);
2641 Float_t sdedx = seed->GetSDEDX(0);
2642 Float_t dedx = seed->GetdEdx();
2643 esd->SetTPCsignal(dedx, sdedx, ndedx);
2645 // add seed to the esd track in Calib level
2647 Bool_t storeFriend = gRandom->Rndm()<(kMaxFriendTracks)/Float_t(nseed);
2648 if (AliTPCReconstructor::StreamLevel()>0 &&storeFriend){
2649 AliTPCseed * seedCopy = new AliTPCseed(*seed, kTRUE);
2650 esd->AddCalibObject(seedCopy);
2655 //printf("problem\n");
2658 //FindKinks(fSeeds,event);
2659 Info("RefitInward","Number of refitted tracks %d",ntracks);
2664 Int_t AliTPCtrackerMI::PropagateBack(AliESDEvent *event)
2667 // back propagation of ESD tracks
2673 PropagateBack(fSeeds);
2674 RemoveUsed2(fSeeds,0.4,0.4,20);
2675 //FindCurling(fSeeds, fEvent,1);
2676 FindSplitted(fSeeds, event,1); // find multi found tracks
2677 if (AliTPCReconstructor::StreamLevel()>5) FindMultiMC(fSeeds, fEvent,1); // find multi found tracks
2680 Int_t nseed = fSeeds->GetEntriesFast();
2682 for (Int_t i=0;i<nseed;i++){
2683 AliTPCseed * seed = (AliTPCseed*) fSeeds->UncheckedAt(i);
2684 if (!seed) continue;
2685 if (seed->GetKinkIndex(0)<0) UpdateKinkQualityM(seed); // update quality informations for kinks
2686 seed->UpdatePoints();
2687 AddCovariance(seed);
2688 AliESDtrack *esd=event->GetTrack(i);
2689 if (seed->GetNumberOfClusters()<60 && seed->GetNumberOfClusters()<(esd->GetTPCclusters(0) -5)*0.8){
2690 AliExternalTrackParam paramIn;
2691 AliExternalTrackParam paramOut;
2692 Int_t ncl = seed->RefitTrack(seed,¶mIn,¶mOut);
2693 if (AliTPCReconstructor::StreamLevel()>2) {
2694 (*fDebugStreamer)<<"RecoverBack"<<
2698 "pout.="<<¶mOut<<
2703 seed->Set(paramOut.GetX(),paramOut.GetAlpha(),paramOut.GetParameter(),paramOut.GetCovariance());
2704 seed->SetNumberOfClusters(ncl);
2707 seed->CookdEdx(0.02,0.6);
2708 CookLabel(seed,0.1); //For comparison only
2709 if (seed->GetNumberOfClusters()>15){
2710 esd->UpdateTrackParams(seed,AliESDtrack::kTPCout);
2711 esd->SetTPCPoints(seed->GetPoints());
2712 esd->SetTPCPointsF(seed->GetNFoundable());
2713 Int_t ndedx = seed->GetNCDEDX(0);
2714 Float_t sdedx = seed->GetSDEDX(0);
2715 Float_t dedx = seed->GetdEdx();
2716 esd->SetTPCsignal(dedx, sdedx, ndedx);
2718 Int_t eventnumber = event->GetEventNumberInFile();// patch 28 fev 06
2719 // This is most likely NOT the event number you'd like to use. It has nothing to do with the 'real' event number
2720 if (AliTPCReconstructor::StreamLevel()>1 && esd) {
2721 (*fDebugStreamer)<<"Cback"<<
2724 "EventNrInFile="<<eventnumber<<
2729 //FindKinks(fSeeds,event);
2730 Info("PropagateBack","Number of back propagated tracks %d",ntracks);
2737 void AliTPCtrackerMI::DeleteSeeds()
2742 Int_t nseed = fSeeds->GetEntriesFast();
2743 for (Int_t i=0;i<nseed;i++){
2744 AliTPCseed * seed = (AliTPCseed*)fSeeds->At(i);
2745 if (seed) delete fSeeds->RemoveAt(i);
2752 void AliTPCtrackerMI::ReadSeeds(AliESDEvent *const event, Int_t direction)
2755 //read seeds from the event
2757 Int_t nentr=event->GetNumberOfTracks();
2759 Info("ReadSeeds", "Number of ESD tracks: %d\n", nentr);
2764 fSeeds = new TObjArray(nentr);
2768 for (Int_t i=0; i<nentr; i++) {
2769 AliESDtrack *esd=event->GetTrack(i);
2770 ULong_t status=esd->GetStatus();
2771 if (!(status&AliESDtrack::kTPCin)) continue;
2772 AliTPCtrack t(*esd);
2773 t.SetNumberOfClusters(0);
2774 // AliTPCseed *seed = new AliTPCseed(t,t.GetAlpha());
2775 AliTPCseed *seed = new AliTPCseed(t/*,t.GetAlpha()*/);
2776 seed->SetUniqueID(esd->GetID());
2777 AddCovariance(seed); //add systematic ucertainty
2778 for (Int_t ikink=0;ikink<3;ikink++) {
2779 Int_t index = esd->GetKinkIndex(ikink);
2780 seed->GetKinkIndexes()[ikink] = index;
2781 if (index==0) continue;
2782 index = TMath::Abs(index);
2783 AliESDkink * kink = fEvent->GetKink(index-1);
2784 if (kink&&esd->GetKinkIndex(ikink)<0){
2785 if ((status & AliESDtrack::kTRDrefit) != 0) kink->SetStatus(1,2);
2786 if ((status & AliESDtrack::kITSout) != 0) kink->SetStatus(1,0);
2788 if (kink&&esd->GetKinkIndex(ikink)>0){
2789 if ((status & AliESDtrack::kTRDrefit) != 0) kink->SetStatus(1,6);
2790 if ((status & AliESDtrack::kITSout) != 0) kink->SetStatus(1,4);
2794 if (((status&AliESDtrack::kITSout)==0)&&(direction==1)) seed->ResetCovariance(10.);
2795 if ( direction ==2 &&(status & AliESDtrack::kTRDrefit) == 0 ) seed->ResetCovariance(10.);
2796 //if ( direction ==2 && ((status & AliESDtrack::kTPCout) == 0) ) {
2797 // fSeeds->AddAt(0,i);
2801 if ( direction ==2 &&(status & AliESDtrack::kTRDrefit) > 0 ) {
2802 Double_t par0[5],par1[5],alpha,x;
2803 esd->GetInnerExternalParameters(alpha,x,par0);
2804 esd->GetExternalParameters(x,par1);
2805 Double_t delta1 = TMath::Abs(par0[4]-par1[4])/(0.000000001+TMath::Abs(par0[4]+par1[4]));
2806 Double_t delta2 = TMath::Abs(par0[3]-par1[3]);
2808 if (esd->GetTRDncls()>0) trdchi2 = esd->GetTRDchi2()/esd->GetTRDncls();
2809 //reset covariance if suspicious
2810 if ( (delta1>0.1) || (delta2>0.006) ||trdchi2>7.)
2811 seed->ResetCovariance(10.);
2816 // rotate to the local coordinate system
2818 fSectors=fInnerSec; fN=fkNIS;
2819 Double_t alpha=seed->GetAlpha() - fSectors->GetAlphaShift();
2820 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
2821 if (alpha < 0. ) alpha += 2.*TMath::Pi();
2822 Int_t ns=Int_t(alpha/fSectors->GetAlpha())%fN;
2823 alpha =ns*fSectors->GetAlpha() + fSectors->GetAlphaShift();
2824 if (alpha<-TMath::Pi()) alpha += 2*TMath::Pi();
2825 if (alpha>=TMath::Pi()) alpha -= 2*TMath::Pi();
2826 alpha-=seed->GetAlpha();
2827 if (!seed->Rotate(alpha)) {
2833 if (esd->GetKinkIndex(0)<=0){
2834 for (Int_t irow=0;irow<160;irow++){
2835 Int_t index = seed->GetClusterIndex2(irow);
2838 AliTPCclusterMI * cl = GetClusterMI(index);
2839 seed->SetClusterPointer(irow,cl);
2841 if ((index & 0x8000)==0){
2842 cl->Use(10); // accepted cluster
2844 cl->Use(6); // close cluster not accepted
2847 Info("ReadSeeds","Not found cluster");
2852 fSeeds->AddAt(seed,i);
2858 //_____________________________________________________________________________
2859 void AliTPCtrackerMI::MakeSeeds3(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, Float_t cuts[4],
2860 Float_t deltay, Int_t ddsec) {
2861 //-----------------------------------------------------------------
2862 // This function creates track seeds.
2863 // SEEDING WITH VERTEX CONSTRAIN
2864 //-----------------------------------------------------------------
2865 // cuts[0] - fP4 cut
2866 // cuts[1] - tan(phi) cut
2867 // cuts[2] - zvertex cut
2868 // cuts[3] - fP3 cut
2876 Double_t x[5], c[15];
2877 // Int_t di = i1-i2;
2879 AliTPCseed * seed = new AliTPCseed();
2880 Double_t alpha=fSectors->GetAlpha(), shift=fSectors->GetAlphaShift();
2881 Double_t cs=cos(alpha), sn=sin(alpha);
2883 // Double_t x1 =fOuterSec->GetX(i1);
2884 //Double_t xx2=fOuterSec->GetX(i2);
2886 Double_t x1 =GetXrow(i1);
2887 Double_t xx2=GetXrow(i2);
2889 Double_t x3=GetX(), y3=GetY(), z3=GetZ();
2891 Int_t imiddle = (i2+i1)/2; //middle pad row index
2892 Double_t xm = GetXrow(imiddle); // radius of middle pad-row
2893 const AliTPCtrackerRow& krm=GetRow(sec,imiddle); //middle pad -row
2897 const AliTPCtrackerRow& kr1=GetRow(ns,i1);
2898 Double_t ymax = GetMaxY(i1)-kr1.GetDeadZone()-1.5;
2899 Double_t ymaxm = GetMaxY(imiddle)-kr1.GetDeadZone()-1.5;
2902 // change cut on curvature if it can't reach this layer
2903 // maximal curvature set to reach it
2904 Double_t dvertexmax = TMath::Sqrt((x1-x3)*(x1-x3)+(ymax+5-y3)*(ymax+5-y3));
2905 if (dvertexmax*0.5*cuts[0]>0.85){
2906 cuts[0] = 0.85/(dvertexmax*0.5+1.);
2908 Double_t r2min = 1/(cuts[0]*cuts[0]); //minimal square of radius given by cut
2911 if (deltay>0) ddsec = 0;
2912 // loop over clusters
2913 for (Int_t is=0; is < kr1; is++) {
2915 if (kr1[is]->IsUsed(10)) continue;
2916 Double_t y1=kr1[is]->GetY(), z1=kr1[is]->GetZ();
2917 //if (TMath::Abs(y1)>ymax) continue;
2919 if (deltay>0 && TMath::Abs(ymax-TMath::Abs(y1))> deltay ) continue; // seed only at the edge
2921 // find possible directions
2922 Float_t anglez = (z1-z3)/(x1-x3);
2923 Float_t extraz = z1 - anglez*(x1-xx2); // extrapolated z
2926 //find rotation angles relative to line given by vertex and point 1
2927 Double_t dvertex2 = (x1-x3)*(x1-x3)+(y1-y3)*(y1-y3);
2928 Double_t dvertex = TMath::Sqrt(dvertex2);
2929 Double_t angle13 = TMath::ATan((y1-y3)/(x1-x3));
2930 Double_t cs13 = cos(-angle13), sn13 = sin(-angle13);
2933 // loop over 2 sectors
2939 Double_t dddz1=0; // direction of delta inclination in z axis
2946 for (Int_t dsec = dsec1; dsec<=dsec2;dsec++){
2947 Int_t sec2 = sec + dsec;
2949 // AliTPCtrackerRow& kr2 = fOuterSec[(sec2+fkNOS)%fkNOS][i2];
2950 //AliTPCtrackerRow& kr2m = fOuterSec[(sec2+fkNOS)%fkNOS][imiddle];
2951 AliTPCtrackerRow& kr2 = GetRow((sec2+fkNOS)%fkNOS,i2);
2952 AliTPCtrackerRow& kr2m = GetRow((sec2+fkNOS)%fkNOS,imiddle);
2953 Int_t index1 = TMath::Max(kr2.Find(extraz-0.6-dddz1*TMath::Abs(z1)*0.05)-1,0);
2954 Int_t index2 = TMath::Min(kr2.Find(extraz+0.6+dddz2*TMath::Abs(z1)*0.05)+1,kr2);
2956 // rotation angles to p1-p3
2957 Double_t cs13r = cos(-angle13+dsec*alpha)/dvertex, sn13r = sin(-angle13+dsec*alpha)/dvertex;
2958 Double_t x2, y2, z2;
2960 // Double_t dymax = maxangle*TMath::Abs(x1-xx2);
2963 Double_t dxx0 = (xx2-x3)*cs13r;
2964 Double_t dyy0 = (xx2-x3)*sn13r;
2965 for (Int_t js=index1; js < index2; js++) {
2966 const AliTPCclusterMI *kcl = kr2[js];
2967 if (kcl->IsUsed(10)) continue;
2969 //calcutate parameters
2971 Double_t yy0 = dyy0 +(kcl->GetY()-y3)*cs13r;
2973 if (TMath::Abs(yy0)<0.000001) continue;
2974 Double_t xx0 = dxx0 -(kcl->GetY()-y3)*sn13r;
2975 Double_t y0 = 0.5*(xx0*xx0+yy0*yy0-xx0)/yy0;
2976 Double_t r02 = (0.25+y0*y0)*dvertex2;
2977 //curvature (radius) cut
2978 if (r02<r2min) continue;
2982 Double_t c0 = 1/TMath::Sqrt(r02);
2986 //Double_t dfi0 = 2.*TMath::ASin(dvertex*c0*0.5);
2987 //Double_t dfi1 = 2.*TMath::ASin(TMath::Sqrt(yy0*yy0+(1-xx0)*(1-xx0))*dvertex*c0*0.5);
2988 Double_t dfi0 = 2.*AliTPCFastMath::FastAsin(dvertex*c0*0.5);
2989 Double_t dfi1 = 2.*AliTPCFastMath::FastAsin(TMath::Sqrt(yy0*yy0+(1-xx0)*(1-xx0))*dvertex*c0*0.5);
2992 Double_t z0 = kcl->GetZ();
2993 Double_t zzzz2 = z1-(z1-z3)*dfi1/dfi0;
2994 if (TMath::Abs(zzzz2-z0)>0.5) continue;
2997 Double_t dip = (z1-z0)*c0/dfi1;
2998 Double_t x0 = (0.5*cs13+y0*sn13)*dvertex*c0;
3009 x2= xx2*cs-y2*sn*dsec;
3010 y2=+xx2*sn*dsec+y2*cs;
3020 // do we have cluster at the middle ?
3022 GetProlongation(x1,xm,x,ym,zm);
3024 AliTPCclusterMI * cm=0;
3025 if (TMath::Abs(ym)-ymaxm<0){
3026 cm = krm.FindNearest2(ym,zm,1.0,0.6,dummy);
3027 if ((!cm) || (cm->IsUsed(10))) {
3032 // rotate y1 to system 0
3033 // get state vector in rotated system
3034 Double_t yr1 = (-0.5*sn13+y0*cs13)*dvertex*c0;
3035 Double_t xr2 = x0*cs+yr1*sn*dsec;
3036 Double_t xr[5]={kcl->GetY(),kcl->GetZ(), xr2, dip, c0};
3038 GetProlongation(xx2,xm,xr,ym,zm);
3039 if (TMath::Abs(ym)-ymaxm<0){
3040 cm = kr2m.FindNearest2(ym,zm,1.0,0.6,dummy);
3041 if ((!cm) || (cm->IsUsed(10))) {
3051 dym = ym - cm->GetY();
3052 dzm = zm - cm->GetZ();
3059 Double_t sy1=kr1[is]->GetSigmaY2()*2., sz1=kr1[is]->GetSigmaZ2()*2.;
3060 Double_t sy2=kcl->GetSigmaY2()*2., sz2=kcl->GetSigmaZ2()*2.;
3061 //Double_t sy3=400*3./12., sy=0.1, sz=0.1;
3062 Double_t sy3=25000*x[4]*x[4]+0.1, sy=0.1, sz=0.1;
3063 //Double_t sy3=25000*x[4]*x[4]*60+0.5, sy=0.1, sz=0.1;
3065 Double_t f40=(F1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
3066 Double_t f42=(F1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
3067 Double_t f43=(F1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
3068 Double_t f20=(F2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
3069 Double_t f22=(F2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
3070 Double_t f23=(F2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
3072 Double_t f30=(F3(x1,y1+sy,x2,y2,z1,z2)-x[3])/sy;
3073 Double_t f31=(F3(x1,y1,x2,y2,z1+sz,z2)-x[3])/sz;
3074 Double_t f32=(F3(x1,y1,x2,y2+sy,z1,z2)-x[3])/sy;
3075 Double_t f34=(F3(x1,y1,x2,y2,z1,z2+sz)-x[3])/sz;
3079 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3080 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3081 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3082 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3083 c[13]=f30*sy1*f40+f32*sy2*f42;
3084 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3086 // if (!BuildSeed(kr1[is],kcl,0,x1,x2,x3,x,c)) continue;
3088 UInt_t index=kr1.GetIndex(is);
3089 seed->~AliTPCseed(); // this does not set the pointer to 0...
3090 AliTPCseed *track=new(seed) AliTPCseed(x1, ns*alpha+shift, x, c, index);
3092 track->SetIsSeeding(kTRUE);
3093 track->SetSeed1(i1);
3094 track->SetSeed2(i2);
3095 track->SetSeedType(3);
3099 FollowProlongation(*track, (i1+i2)/2,1);
3100 Int_t foundable,found,shared;
3101 track->GetClusterStatistic((i1+i2)/2,i1, found, foundable, shared, kTRUE);
3102 if ((found<0.55*foundable) || shared>0.5*found || (track->GetSigmaY2()+track->GetSigmaZ2())>0.5){
3104 seed->~AliTPCseed();
3110 FollowProlongation(*track, i2,1);
3114 track->SetBConstrain(1);
3115 // track->fLastPoint = i1+fInnerSec->GetNRows(); // first cluster in track position
3116 track->SetLastPoint(i1); // first cluster in track position
3117 track->SetFirstPoint(track->GetLastPoint());
3119 if (track->GetNumberOfClusters()<(i1-i2)*0.5 ||
3120 track->GetNumberOfClusters() < track->GetNFoundable()*0.6 ||
3121 track->GetNShared()>0.4*track->GetNumberOfClusters() ) {
3123 seed->~AliTPCseed();
3127 // Z VERTEX CONDITION
3128 Double_t zv, bz=GetBz();
3129 if ( !track->GetZAt(0.,bz,zv) ) continue;
3130 if (TMath::Abs(zv-z3)>cuts[2]) {
3131 FollowProlongation(*track, TMath::Max(i2-20,0));
3132 if ( !track->GetZAt(0.,bz,zv) ) continue;
3133 if (TMath::Abs(zv-z3)>cuts[2]){
3134 FollowProlongation(*track, TMath::Max(i2-40,0));
3135 if ( !track->GetZAt(0.,bz,zv) ) continue;
3136 if (TMath::Abs(zv-z3)>cuts[2] &&(track->GetNumberOfClusters() > track->GetNFoundable()*0.7)){
3137 // make seed without constrain
3138 AliTPCseed * track2 = MakeSeed(track,0.2,0.5,1.);
3139 FollowProlongation(*track2, i2,1);
3140 track2->SetBConstrain(kFALSE);
3141 track2->SetSeedType(1);
3142 arr->AddLast(track2);
3144 seed->~AliTPCseed();
3149 seed->~AliTPCseed();
3156 track->SetSeedType(0);
3157 arr->AddLast(track);
3158 seed = new AliTPCseed;
3160 // don't consider other combinations
3161 if (track->GetNumberOfClusters() > track->GetNFoundable()*0.8)
3167 Info("MakeSeeds3","\nSeeding statistic:\t%d\t%d\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin,nout1,nout2);
3173 void AliTPCtrackerMI::MakeSeeds5(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, Float_t cuts[4],
3178 //-----------------------------------------------------------------
3179 // This function creates track seeds.
3180 //-----------------------------------------------------------------
3181 // cuts[0] - fP4 cut
3182 // cuts[1] - tan(phi) cut
3183 // cuts[2] - zvertex cut
3184 // cuts[3] - fP3 cut
3194 Double_t x[5], c[15];
3196 // make temporary seed
3197 AliTPCseed * seed = new AliTPCseed;
3198 Double_t alpha=fOuterSec->GetAlpha(), shift=fOuterSec->GetAlphaShift();
3199 // Double_t cs=cos(alpha), sn=sin(alpha);
3204 Double_t x1 = GetXrow(i1-1);
3205 const AliTPCtrackerRow& kr1=GetRow(sec,i1-1);
3206 Double_t y1max = GetMaxY(i1-1)-kr1.GetDeadZone()-1.5;
3208 Double_t x1p = GetXrow(i1);
3209 const AliTPCtrackerRow& kr1p=GetRow(sec,i1);
3211 Double_t x1m = GetXrow(i1-2);
3212 const AliTPCtrackerRow& kr1m=GetRow(sec,i1-2);
3215 //last 3 padrow for seeding
3216 AliTPCtrackerRow& kr3 = GetRow((sec+fkNOS)%fkNOS,i1-7);
3217 Double_t x3 = GetXrow(i1-7);
3218 // Double_t y3max= GetMaxY(i1-7)-kr3.fDeadZone-1.5;
3220 AliTPCtrackerRow& kr3p = GetRow((sec+fkNOS)%fkNOS,i1-6);
3221 Double_t x3p = GetXrow(i1-6);
3223 AliTPCtrackerRow& kr3m = GetRow((sec+fkNOS)%fkNOS,i1-8);
3224 Double_t x3m = GetXrow(i1-8);
3229 Int_t im = i1-4; //middle pad row index
3230 Double_t xm = GetXrow(im); // radius of middle pad-row
3231 const AliTPCtrackerRow& krm=GetRow(sec,im); //middle pad -row
3232 // Double_t ymmax = GetMaxY(im)-kr1.fDeadZone-1.5;
3235 Double_t deltax = x1-x3;
3236 Double_t dymax = deltax*cuts[1];
3237 Double_t dzmax = deltax*cuts[3];
3239 // loop over clusters
3240 for (Int_t is=0; is < kr1; is++) {
3242 if (kr1[is]->IsUsed(10)) continue;
3243 Double_t y1=kr1[is]->GetY(), z1=kr1[is]->GetZ();
3245 if (deltay>0 && TMath::Abs(y1max-TMath::Abs(y1))> deltay ) continue; // seed only at the edge
3247 Int_t index1 = TMath::Max(kr3.Find(z1-dzmax)-1,0);
3248 Int_t index2 = TMath::Min(kr3.Find(z1+dzmax)+1,kr3);
3254 for (Int_t js=index1; js < index2; js++) {
3255 const AliTPCclusterMI *kcl = kr3[js];
3256 if (kcl->IsUsed(10)) continue;
3258 // apply angular cuts
3259 if (TMath::Abs(y1-y3)>dymax) continue;
3262 if (TMath::Abs(z1-z3)>dzmax) continue;
3264 Double_t angley = (y1-y3)/(x1-x3);
3265 Double_t anglez = (z1-z3)/(x1-x3);
3267 Double_t erry = TMath::Abs(angley)*(x1-x1m)*0.5+0.5;
3268 Double_t errz = TMath::Abs(anglez)*(x1-x1m)*0.5+0.5;
3270 Double_t yyym = angley*(xm-x1)+y1;
3271 Double_t zzzm = anglez*(xm-x1)+z1;
3273 const AliTPCclusterMI *kcm = krm.FindNearest2(yyym,zzzm,erry,errz,index);
3275 if (kcm->IsUsed(10)) continue;
3277 erry = TMath::Abs(angley)*(x1-x1m)*0.4+0.5;
3278 errz = TMath::Abs(anglez)*(x1-x1m)*0.4+0.5;
3285 // look around first
3286 const AliTPCclusterMI *kc1m = kr1m.FindNearest2(angley*(x1m-x1)+y1,
3292 if (kc1m->IsUsed(10)) used++;
3294 const AliTPCclusterMI *kc1p = kr1p.FindNearest2(angley*(x1p-x1)+y1,
3300 if (kc1p->IsUsed(10)) used++;
3302 if (used>1) continue;
3303 if (found<1) continue;
3307 const AliTPCclusterMI *kc3m = kr3m.FindNearest2(angley*(x3m-x3)+y3,
3313 if (kc3m->IsUsed(10)) used++;
3317 const AliTPCclusterMI *kc3p = kr3p.FindNearest2(angley*(x3p-x3)+y3,
3323 if (kc3p->IsUsed(10)) used++;
3327 if (used>1) continue;
3328 if (found<3) continue;
3338 x[4]=F1(x1,y1,x2,y2,x3,y3);
3339 //if (TMath::Abs(x[4]) >= cuts[0]) continue;
3342 x[2]=F2(x1,y1,x2,y2,x3,y3);
3345 x[3]=F3n(x1,y1,x2,y2,z1,z2,x[4]);
3346 //if (TMath::Abs(x[3]) > cuts[3]) continue;
3350 Double_t sy1=0.1, sz1=0.1;
3351 Double_t sy2=0.1, sz2=0.1;
3352 Double_t sy3=0.1, sy=0.1, sz=0.1;
3354 Double_t f40=(F1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
3355 Double_t f42=(F1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
3356 Double_t f43=(F1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
3357 Double_t f20=(F2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
3358 Double_t f22=(F2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
3359 Double_t f23=(F2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
3361 Double_t f30=(F3(x1,y1+sy,x2,y2,z1,z2)-x[3])/sy;
3362 Double_t f31=(F3(x1,y1,x2,y2,z1+sz,z2)-x[3])/sz;
3363 Double_t f32=(F3(x1,y1,x2,y2+sy,z1,z2)-x[3])/sy;
3364 Double_t f34=(F3(x1,y1,x2,y2,z1,z2+sz)-x[3])/sz;
3368 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3369 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3370 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3371 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3372 c[13]=f30*sy1*f40+f32*sy2*f42;
3373 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3375 // if (!BuildSeed(kr1[is],kcl,0,x1,x2,x3,x,c)) continue;
3377 index=kr1.GetIndex(is);
3378 seed->~AliTPCseed();
3379 AliTPCseed *track=new(seed) AliTPCseed(x1, sec*alpha+shift, x, c, index);
3381 track->SetIsSeeding(kTRUE);
3384 FollowProlongation(*track, i1-7,1);
3385 if (track->GetNumberOfClusters() < track->GetNFoundable()*0.75 ||
3386 track->GetNShared()>0.6*track->GetNumberOfClusters() || ( track->GetSigmaY2()+ track->GetSigmaZ2())>0.6){
3388 seed->~AliTPCseed();
3394 FollowProlongation(*track, i2,1);
3395 track->SetBConstrain(0);
3396 track->SetLastPoint(i1+fInnerSec->GetNRows()); // first cluster in track position
3397 track->SetFirstPoint(track->GetLastPoint());
3399 if (track->GetNumberOfClusters()<(i1-i2)*0.5 ||
3400 track->GetNumberOfClusters()<track->GetNFoundable()*0.7 ||
3401 track->GetNShared()>2. || track->GetChi2()/track->GetNumberOfClusters()>6 || ( track->GetSigmaY2()+ track->GetSigmaZ2())>0.5 ) {
3403 seed->~AliTPCseed();
3408 FollowProlongation(*track, TMath::Max(i2-10,0),1);
3409 AliTPCseed * track2 = MakeSeed(track,0.2,0.5,0.9);
3410 FollowProlongation(*track2, i2,1);
3411 track2->SetBConstrain(kFALSE);
3412 track2->SetSeedType(4);
3413 arr->AddLast(track2);
3415 seed->~AliTPCseed();
3419 //arr->AddLast(track);
3420 //seed = new AliTPCseed;
3426 Info("MakeSeeds5","\nSeeding statiistic:\t%d\t%d\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin,nout1,nout2,nout3);
3432 //_____________________________________________________________________________
3433 void AliTPCtrackerMI::MakeSeeds2(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, Float_t */*cuts[4]*/,
3434 Float_t deltay, Bool_t /*bconstrain*/) {
3435 //-----------------------------------------------------------------
3436 // This function creates track seeds - without vertex constraint
3437 //-----------------------------------------------------------------
3438 // cuts[0] - fP4 cut - not applied
3439 // cuts[1] - tan(phi) cut
3440 // cuts[2] - zvertex cut - not applied
3441 // cuts[3] - fP3 cut
3451 Double_t alpha=fOuterSec->GetAlpha(), shift=fOuterSec->GetAlphaShift();
3452 // Double_t cs=cos(alpha), sn=sin(alpha);
3453 Int_t row0 = (i1+i2)/2;
3454 Int_t drow = (i1-i2)/2;
3455 const AliTPCtrackerRow& kr0=fSectors[sec][row0];
3456 AliTPCtrackerRow * kr=0;
3458 AliTPCpolyTrack polytrack;
3459 Int_t nclusters=fSectors[sec][row0];
3460 AliTPCseed * seed = new AliTPCseed;
3465 for (Int_t is=0; is < nclusters; is++) { //LOOP over clusters
3467 Int_t nfoundable =0;
3468 for (Int_t iter =1; iter<2; iter++){ //iterations
3469 const AliTPCtrackerRow& krm=fSectors[sec][row0-iter];
3470 const AliTPCtrackerRow& krp=fSectors[sec][row0+iter];
3471 const AliTPCclusterMI * cl= kr0[is];
3473 if (cl->IsUsed(10)) {
3479 Double_t x = kr0.GetX();
3480 // Initialization of the polytrack
3485 Double_t y0= cl->GetY();
3486 Double_t z0= cl->GetZ();
3490 Double_t ymax = fSectors->GetMaxY(row0)-kr0.GetDeadZone()-1.5;
3491 if (deltay>0 && TMath::Abs(ymax-TMath::Abs(y0))> deltay ) continue; // seed only at the edge
3493 erry = (0.5)*cl->GetSigmaY2()/TMath::Sqrt(cl->GetQ())*6;
3494 errz = (0.5)*cl->GetSigmaZ2()/TMath::Sqrt(cl->GetQ())*6;
3495 polytrack.AddPoint(x,y0,z0,erry, errz);
3498 if (cl->IsUsed(10)) sumused++;
3501 Float_t roady = (5*TMath::Sqrt(cl->GetSigmaY2()+0.2)+1.)*iter;
3502 Float_t roadz = (5*TMath::Sqrt(cl->GetSigmaZ2()+0.2)+1.)*iter;
3505 AliTPCclusterMI * cl1 = krm.FindNearest(y0,z0,roady,roadz);
3506 if (cl1 && TMath::Abs(ymax-TMath::Abs(y0))) {
3507 erry = (0.5)*cl1->GetSigmaY2()/TMath::Sqrt(cl1->GetQ())*3;
3508 errz = (0.5)*cl1->GetSigmaZ2()/TMath::Sqrt(cl1->GetQ())*3;
3509 if (cl1->IsUsed(10)) sumused++;
3510 polytrack.AddPoint(x,cl1->GetY(),cl1->GetZ(),erry,errz);
3514 AliTPCclusterMI * cl2 = krp.FindNearest(y0,z0,roady,roadz);
3516 erry = (0.5)*cl2->GetSigmaY2()/TMath::Sqrt(cl2->GetQ())*3;
3517 errz = (0.5)*cl2->GetSigmaZ2()/TMath::Sqrt(cl2->GetQ())*3;
3518 if (cl2->IsUsed(10)) sumused++;
3519 polytrack.AddPoint(x,cl2->GetY(),cl2->GetZ(),erry,errz);
3522 if (sumused>0) continue;
3524 polytrack.UpdateParameters();
3530 nfoundable = polytrack.GetN();
3531 nfound = nfoundable;
3533 for (Int_t ddrow = iter+1; ddrow<drow;ddrow++){
3534 Float_t maxdist = 0.8*(1.+3./(ddrow));
3535 for (Int_t delta = -1;delta<=1;delta+=2){
3536 Int_t row = row0+ddrow*delta;
3537 kr = &(fSectors[sec][row]);
3538 Double_t xn = kr->GetX();
3539 Double_t ymax1 = fSectors->GetMaxY(row)-kr->GetDeadZone()-1.5;
3540 polytrack.GetFitPoint(xn,yn,zn);
3541 if (TMath::Abs(yn)>ymax1) continue;
3543 AliTPCclusterMI * cln = kr->FindNearest(yn,zn,roady,roadz);
3545 Float_t dist = TMath::Sqrt( (yn-cln->GetY())*(yn-cln->GetY())+(zn-cln->GetZ())*(zn-cln->GetZ()));
3548 erry = (dist+0.3)*cln->GetSigmaY2()/TMath::Sqrt(cln->GetQ())*(1.+1./(ddrow));
3549 errz = (dist+0.3)*cln->GetSigmaZ2()/TMath::Sqrt(cln->GetQ())*(1.+1./(ddrow));
3550 if (cln->IsUsed(10)) {
3551 // printf("used\n");
3559 polytrack.AddPoint(xn,cln->GetY(),cln->GetZ(),erry, errz);
3564 if ( (sumused>3) || (sumused>0.5*nfound) || (nfound<0.6*nfoundable)) break;
3565 polytrack.UpdateParameters();
3568 if ( (sumused>3) || (sumused>0.5*nfound)) {
3569 //printf("sumused %d\n",sumused);
3574 polytrack.GetFitDerivation(kr0.GetX(),dy,dz);
3575 AliTPCpolyTrack track2;
3577 polytrack.Refit(track2,0.5+TMath::Abs(dy)*0.3,0.4+TMath::Abs(dz)*0.3);
3578 if (track2.GetN()<0.5*nfoundable) continue;
3581 if ((nfound>0.6*nfoundable) &&( nfoundable>0.4*(i1-i2))) {
3583 // test seed with and without constrain
3584 for (Int_t constrain=0; constrain<=0;constrain++){
3585 // add polytrack candidate
3587 Double_t x[5], c[15];
3588 Double_t x1,x2,x3,y1,y2,y3,z1,z2,z3;
3589 track2.GetBoundaries(x3,x1);
3591 track2.GetFitPoint(x1,y1,z1);
3592 track2.GetFitPoint(x2,y2,z2);
3593 track2.GetFitPoint(x3,y3,z3);
3595 //is track pointing to the vertex ?
3598 polytrack.GetFitPoint(x0,y0,z0);
3611 x[4]=F1(x1,y1,x2,y2,x3,y3);
3613 // if (TMath::Abs(x[4]) >= cuts[0]) continue; //
3614 x[2]=F2(x1,y1,x2,y2,x3,y3);
3616 //if (TMath::Abs(x[4]*x1-x[2]) >= cuts[1]) continue;
3617 //x[3]=F3(x1,y1,x2,y2,z1,z2);
3618 x[3]=F3n(x1,y1,x3,y3,z1,z3,x[4]);
3619 //if (TMath::Abs(x[3]) > cuts[3]) continue;
3622 Double_t sy =0.1, sz =0.1;
3623 Double_t sy1=0.02, sz1=0.02;
3624 Double_t sy2=0.02, sz2=0.02;
3628 sy3=25000*x[4]*x[4]+0.1, sy=0.1, sz=0.1;
3631 Double_t f40=(F1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
3632 Double_t f42=(F1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
3633 Double_t f43=(F1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
3634 Double_t f20=(F2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
3635 Double_t f22=(F2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
3636 Double_t f23=(F2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
3638 Double_t f30=(F3(x1,y1+sy,x3,y3,z1,z3)-x[3])/sy;
3639 Double_t f31=(F3(x1,y1,x3,y3,z1+sz,z3)-x[3])/sz;
3640 Double_t f32=(F3(x1,y1,x3,y3+sy,z1,z3)-x[3])/sy;
3641 Double_t f34=(F3(x1,y1,x3,y3,z1,z3+sz)-x[3])/sz;
3646 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3647 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3648 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3649 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3650 c[13]=f30*sy1*f40+f32*sy2*f42;
3651 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3653 //Int_t row1 = fSectors->GetRowNumber(x1);
3654 Int_t row1 = GetRowNumber(x1);
3658 seed->~AliTPCseed();
3659 AliTPCseed *track=new(seed) AliTPCseed(x1,sec*alpha+shift,x,c,index);
3660 track->SetIsSeeding(kTRUE);
3661 Int_t rc=FollowProlongation(*track, i2);
3662 if (constrain) track->SetBConstrain(1);
3664 track->SetBConstrain(0);
3665 track->SetLastPoint(row1+fInnerSec->GetNRows()); // first cluster in track position
3666 track->SetFirstPoint(track->GetLastPoint());
3668 if (rc==0 || track->GetNumberOfClusters()<(i1-i2)*0.5 ||
3669 track->GetNumberOfClusters() < track->GetNFoundable()*0.6 ||
3670 track->GetNShared()>0.4*track->GetNumberOfClusters()) {
3673 seed->~AliTPCseed();
3676 arr->AddLast(track);
3677 seed = new AliTPCseed;
3681 } // if accepted seed
3684 Info("MakeSeeds2","\nSeeding statiistic:\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin3);
3690 AliTPCseed *AliTPCtrackerMI::MakeSeed(AliTPCseed *const track, Float_t r0, Float_t r1, Float_t r2)
3694 //reseed using track points
3695 Int_t p0 = int(r0*track->GetNumberOfClusters()); // point 0
3696 Int_t p1 = int(r1*track->GetNumberOfClusters());
3697 Int_t p2 = int(r2*track->GetNumberOfClusters()); // last point
3699 Double_t x0[3],x1[3],x2[3];
3700 for (Int_t i=0;i<3;i++){
3706 // find track position at given ratio of the length
3707 Int_t sec0=0, sec1=0, sec2=0;
3710 for (Int_t i=0;i<160;i++){
3711 if (track->GetClusterPointer(i)){
3713 AliTPCTrackerPoint *trpoint =track->GetTrackPoint(i);
3714 if ( (index<p0) || x0[0]<0 ){
3715 if (trpoint->GetX()>1){
3716 clindex = track->GetClusterIndex2(i);
3718 x0[0] = trpoint->GetX();
3719 x0[1] = trpoint->GetY();
3720 x0[2] = trpoint->GetZ();
3721 sec0 = ((clindex&0xff000000)>>24)%18;
3726 if ( (index<p1) &&(trpoint->GetX()>1)){
3727 clindex = track->GetClusterIndex2(i);
3729 x1[0] = trpoint->GetX();
3730 x1[1] = trpoint->GetY();
3731 x1[2] = trpoint->GetZ();
3732 sec1 = ((clindex&0xff000000)>>24)%18;
3735 if ( (index<p2) &&(trpoint->GetX()>1)){
3736 clindex = track->GetClusterIndex2(i);
3738 x2[0] = trpoint->GetX();
3739 x2[1] = trpoint->GetY();
3740 x2[2] = trpoint->GetZ();
3741 sec2 = ((clindex&0xff000000)>>24)%18;
3748 Double_t alpha, cs,sn, xx2,yy2;
3750 alpha = (sec1-sec2)*fSectors->GetAlpha();
3751 cs = TMath::Cos(alpha);
3752 sn = TMath::Sin(alpha);
3753 xx2= x1[0]*cs-x1[1]*sn;
3754 yy2= x1[0]*sn+x1[1]*cs;
3758 alpha = (sec0-sec2)*fSectors->GetAlpha();
3759 cs = TMath::Cos(alpha);
3760 sn = TMath::Sin(alpha);
3761 xx2= x0[0]*cs-x0[1]*sn;
3762 yy2= x0[0]*sn+x0[1]*cs;
3768 Double_t x[5],c[15];
3772 x[4]=F1(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]);
3773 // if (x[4]>1) return 0;
3774 x[2]=F2(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]);
3775 x[3]=F3n(x2[0],x2[1],x0[0],x0[1],x2[2],x0[2],x[4]);
3776 //if (TMath::Abs(x[3]) > 2.2) return 0;
3777 //if (TMath::Abs(x[2]) > 1.99) return 0;
3779 Double_t sy =0.1, sz =0.1;
3781 Double_t sy1=0.02+track->GetSigmaY2(), sz1=0.02+track->GetSigmaZ2();
3782 Double_t sy2=0.01+track->GetSigmaY2(), sz2=0.01+track->GetSigmaZ2();
3783 Double_t sy3=0.01+track->GetSigmaY2();
3785 Double_t f40=(F1(x2[0],x2[1]+sy,x1[0],x1[1],x0[0],x0[1])-x[4])/sy;
3786 Double_t f42=(F1(x2[0],x2[1],x1[0],x1[1]+sy,x0[0],x0[1])-x[4])/sy;
3787 Double_t f43=(F1(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]+sy)-x[4])/sy;
3788 Double_t f20=(F2(x2[0],x2[1]+sy,x1[0],x1[1],x0[0],x0[1])-x[2])/sy;
3789 Double_t f22=(F2(x2[0],x2[1],x1[0],x1[1]+sy,x0[0],x0[1])-x[2])/sy;
3790 Double_t f23=(F2(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]+sy)-x[2])/sy;
3792 Double_t f30=(F3(x2[0],x2[1]+sy,x0[0],x0[1],x2[2],x0[2])-x[3])/sy;
3793 Double_t f31=(F3(x2[0],x2[1],x0[0],x0[1],x2[2]+sz,x0[2])-x[3])/sz;
3794 Double_t f32=(F3(x2[0],x2[1],x0[0],x0[1]+sy,x2[2],x0[2])-x[3])/sy;
3795 Double_t f34=(F3(x2[0],x2[1],x0[0],x0[1],x2[2],x0[2]+sz)-x[3])/sz;
3800 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3801 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3802 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3803 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3804 c[13]=f30*sy1*f40+f32*sy2*f42;
3805 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3807 // Int_t row1 = fSectors->GetRowNumber(x2[0]);
3808 AliTPCseed *seed=new AliTPCseed(x2[0], sec2*fSectors->GetAlpha()+fSectors->GetAlphaShift(), x, c, 0);
3809 // Double_t y0,z0,y1,z1, y2,z2;
3810 //seed->GetProlongation(x0[0],y0,z0);
3811 // seed->GetProlongation(x1[0],y1,z1);
3812 //seed->GetProlongation(x2[0],y2,z2);
3814 seed->SetLastPoint(pp2);
3815 seed->SetFirstPoint(pp2);
3822 AliTPCseed *AliTPCtrackerMI::ReSeed(const AliTPCseed *track, Float_t r0, Float_t r1, Float_t r2)
3826 //reseed using founded clusters
3828 // Find the number of clusters
3829 Int_t nclusters = 0;
3830 for (Int_t irow=0;irow<160;irow++){
3831 if (track->GetClusterIndex(irow)>0) nclusters++;
3835 ipos[0] = TMath::Max(int(r0*nclusters),0); // point 0 cluster
3836 ipos[1] = TMath::Min(int(r1*nclusters),nclusters-1); //
3837 ipos[2] = TMath::Min(int(r2*nclusters),nclusters-1); // last point
3841 Int_t row[3],sec[3]={0,0,0};
3843 // find track row position at given ratio of the length
3845 for (Int_t irow=0;irow<160;irow++){
3846 if (track->GetClusterIndex2(irow)<0) continue;
3848 for (Int_t ipoint=0;ipoint<3;ipoint++){
3849 if (index<=ipos[ipoint]) row[ipoint] = irow;
3853 //Get cluster and sector position
3854 for (Int_t ipoint=0;ipoint<3;ipoint++){
3855 Int_t clindex = track->GetClusterIndex2(row[ipoint]);
3856 AliTPCclusterMI * cl = GetClusterMI(clindex);
3859 // AliTPCclusterMI * cl = GetClusterMI(clindex);
3862 sec[ipoint] = ((clindex&0xff000000)>>24)%18;
3863 xyz[ipoint][0] = GetXrow(row[ipoint]);
3864 xyz[ipoint][1] = cl->GetY();
3865 xyz[ipoint][2] = cl->GetZ();
3869 // Calculate seed state vector and covariance matrix
3871 Double_t alpha, cs,sn, xx2,yy2;
3873 alpha = (sec[1]-sec[2])*fSectors->GetAlpha();
3874 cs = TMath::Cos(alpha);
3875 sn = TMath::Sin(alpha);
3876 xx2= xyz[1][0]*cs-xyz[1][1]*sn;
3877 yy2= xyz[1][0]*sn+xyz[1][1]*cs;
3881 alpha = (sec[0]-sec[2])*fSectors->GetAlpha();
3882 cs = TMath::Cos(alpha);
3883 sn = TMath::Sin(alpha);
3884 xx2= xyz[0][0]*cs-xyz[0][1]*sn;
3885 yy2= xyz[0][0]*sn+xyz[0][1]*cs;
3891 Double_t x[5],c[15];
3895 x[4]=F1(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
3896 x[2]=F2(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
3897 x[3]=F3n(xyz[2][0],xyz[2][1],xyz[0][0],xyz[0][1],xyz[2][2],xyz[0][2],x[4]);
3899 Double_t sy =0.1, sz =0.1;
3901 Double_t sy1=0.2, sz1=0.2;
3902 Double_t sy2=0.2, sz2=0.2;
3905 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;
3906 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;
3907 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;
3908 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;
3909 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;
3910 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;
3912 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;
3913 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;
3914 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;
3915 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;
3920 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3921 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3922 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3923 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3924 c[13]=f30*sy1*f40+f32*sy2*f42;
3925 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3927 // Int_t row1 = fSectors->GetRowNumber(xyz[2][0]);
3928 AliTPCseed *seed=new AliTPCseed(xyz[2][0], sec[2]*fSectors->GetAlpha()+fSectors->GetAlphaShift(), x, c, 0);
3929 seed->SetLastPoint(row[2]);
3930 seed->SetFirstPoint(row[2]);
3935 AliTPCseed *AliTPCtrackerMI::ReSeed(AliTPCseed *track,Int_t r0, Bool_t forward)
3939 //reseed using founded clusters
3942 Int_t row[3]={0,0,0};
3943 Int_t sec[3]={0,0,0};
3945 // forward direction
3947 for (Int_t irow=r0;irow<160;irow++){
3948 if (track->GetClusterIndex(irow)>0){
3953 for (Int_t irow=160;irow>r0;irow--){
3954 if (track->GetClusterIndex(irow)>0){
3959 for (Int_t irow=row[2]-15;irow>row[0];irow--){
3960 if (track->GetClusterIndex(irow)>0){
3968 for (Int_t irow=0;irow<r0;irow++){
3969 if (track->GetClusterIndex(irow)>0){
3974 for (Int_t irow=r0;irow>0;irow--){
3975 if (track->GetClusterIndex(irow)>0){
3980 for (Int_t irow=row[2]-15;irow>row[0];irow--){
3981 if (track->GetClusterIndex(irow)>0){
3988 if ((row[2]-row[0])<20) return 0;
3989 if (row[1]==0) return 0;
3992 //Get cluster and sector position
3993 for (Int_t ipoint=0;ipoint<3;ipoint++){
3994 Int_t clindex = track->GetClusterIndex2(row[ipoint]);
3995 AliTPCclusterMI * cl = GetClusterMI(clindex);
3998 // AliTPCclusterMI * cl = GetClusterMI(clindex);
4001 sec[ipoint] = ((clindex&0xff000000)>>24)%18;
4002 xyz[ipoint][0] = GetXrow(row[ipoint]);
4003 AliTPCTrackerPoint * point = track->GetTrackPoint(row[ipoint]);
4004 if (point&&ipoint<2){
4006 xyz[ipoint][1] = point->GetY();
4007 xyz[ipoint][2] = point->GetZ();
4010 xyz[ipoint][1] = cl->GetY();
4011 xyz[ipoint][2] = cl->GetZ();
4018 // Calculate seed state vector and covariance matrix
4020 Double_t alpha, cs,sn, xx2,yy2;
4022 alpha = (sec[1]-sec[2])*fSectors->GetAlpha();
4023 cs = TMath::Cos(alpha);
4024 sn = TMath::Sin(alpha);
4025 xx2= xyz[1][0]*cs-xyz[1][1]*sn;
4026 yy2= xyz[1][0]*sn+xyz[1][1]*cs;
4030 alpha = (sec[0]-sec[2])*fSectors->GetAlpha();
4031 cs = TMath::Cos(alpha);
4032 sn = TMath::Sin(alpha);
4033 xx2= xyz[0][0]*cs-xyz[0][1]*sn;
4034 yy2= xyz[0][0]*sn+xyz[0][1]*cs;
4040 Double_t x[5],c[15];
4044 x[4]=F1(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
4045 x[2]=F2(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
4046 x[3]=F3n(xyz[2][0],xyz[2][1],xyz[0][0],xyz[0][1],xyz[2][2],xyz[0][2],x[4]);
4048 Double_t sy =0.1, sz =0.1;
4050 Double_t sy1=0.2, sz1=0.2;
4051 Double_t sy2=0.2, sz2=0.2;
4054 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;
4055 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;
4056 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;
4057 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;
4058 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;
4059 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;
4061 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;
4062 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;
4063 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;
4064 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;
4069 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
4070 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
4071 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
4072 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
4073 c[13]=f30*sy1*f40+f32*sy2*f42;
4074 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
4076 // Int_t row1 = fSectors->GetRowNumber(xyz[2][0]);
4077 AliTPCseed *seed=new AliTPCseed(xyz[2][0], sec[2]*fSectors->GetAlpha()+fSectors->GetAlphaShift(), x, c, 0);
4078 seed->SetLastPoint(row[2]);
4079 seed->SetFirstPoint(row[2]);
4080 for (Int_t i=row[0];i<row[2];i++){
4081 seed->SetClusterIndex(i, track->GetClusterIndex(i));
4089 void AliTPCtrackerMI::FindMultiMC(const TObjArray * array, AliESDEvent */*esd*/, Int_t iter)
4092 // find multi tracks - THIS FUNCTION IS ONLY FOR DEBUG PURPOSES
4094 // Use AliTPCReconstructor::StreamLevel()>2 if you want to tune parameters - cuts
4096 // Two reasons to have multiple find tracks
4097 // 1. Curling tracks can be find more than once
4098 // 2. Splitted tracks
4099 // a.) Multiple seeding to increase tracking efficiency - (~ 100% reached)
4100 // b.) Edge effect on the sector boundaries
4103 // Algorithm done in 2 phases - because of CPU consumption
4104 // it is n^2 algorithm - for lead-lead 20000x20000 combination are investigated
4106 // Algorihm for curling tracks sign:
4107 // 1 phase -makes a very rough fast cuts to minimize combinatorics
4108 // a.) opposite sign
4109 // b.) one of the tracks - not pointing to the primary vertex -
4110 // c.) delta tan(theta)
4112 // 2 phase - calculates DCA between tracks - time consument
4117 // General cuts - for splitted tracks and for curling tracks
4119 const Float_t kMaxdPhi = 0.2; // maximal distance in phi
4121 // Curling tracks cuts
4126 Int_t nentries = array->GetEntriesFast();
4127 AliHelix *helixes = new AliHelix[nentries];
4128 Float_t *xm = new Float_t[nentries];
4129 Float_t *dz0 = new Float_t[nentries];
4130 Float_t *dz1 = new Float_t[nentries];
4136 // Find track COG in x direction - point with best defined parameters
4138 for (Int_t i=0;i<nentries;i++){
4139 AliTPCseed* track = (AliTPCseed*)array->At(i);
4140 if (!track) continue;
4141 track->SetCircular(0);
4142 new (&helixes[i]) AliHelix(*track);
4146 track->GetDZ(GetX(),GetY(),GetZ(),GetBz(),dz);
4149 for (Int_t icl=0; icl<160; icl++){
4150 AliTPCclusterMI * cl = track->GetClusterPointer(icl);
4156 if (ncl>0) xm[i]/=Float_t(ncl);
4159 for (Int_t i0=0;i0<nentries;i0++){
4160 AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
4161 if (!track0) continue;
4162 Float_t xc0 = helixes[i0].GetHelix(6);
4163 Float_t yc0 = helixes[i0].GetHelix(7);
4164 Float_t r0 = helixes[i0].GetHelix(8);
4165 Float_t rc0 = TMath::Sqrt(xc0*xc0+yc0*yc0);
4166 Float_t fi0 = TMath::ATan2(yc0,xc0);
4168 for (Int_t i1=i0+1;i1<nentries;i1++){
4169 AliTPCseed * track1 = (AliTPCseed*)array->At(i1);
4170 if (!track1) continue;
4171 Int_t lab0=track0->GetLabel();
4172 Int_t lab1=track1->GetLabel();
4173 if (TMath::Abs(lab0)!=TMath::Abs(lab1)) continue;
4175 Float_t xc1 = helixes[i1].GetHelix(6);
4176 Float_t yc1 = helixes[i1].GetHelix(7);
4177 Float_t r1 = helixes[i1].GetHelix(8);
4178 Float_t rc1 = TMath::Sqrt(xc1*xc1+yc1*yc1);
4179 Float_t fi1 = TMath::ATan2(yc1,xc1);
4181 Float_t dfi = fi0-fi1;
4184 if (dfi>1.5*TMath::Pi()) dfi-=TMath::Pi(); // take care about edge effect
4185 if (dfi<-1.5*TMath::Pi()) dfi+=TMath::Pi(); //
4186 if (TMath::Abs(dfi)>kMaxdPhi&&helixes[i0].GetHelix(4)*helixes[i1].GetHelix(4)<0){
4188 // if short tracks with undefined sign
4189 fi1 = -TMath::ATan2(yc1,-xc1);
4192 Float_t dtheta = TMath::Abs(track0->GetTgl()-track1->GetTgl())<TMath::Abs(track0->GetTgl()+track1->GetTgl())? track0->GetTgl()-track1->GetTgl():track0->GetTgl()+track1->GetTgl();
4195 // debug stream to tune "fast cuts"
4197 Double_t dist[3]; // distance at X
4198 Double_t mdist[3]={0,0,0}; // mean distance X+-40cm
4199 track0->GetDistance(track1,0.5*(xm[i0]+xm[i1])-40.,dist,AliTracker::GetBz());
4200 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[i]);
4201 track0->GetDistance(track1,0.5*(xm[i0]+xm[i1])+40.,dist,AliTracker::GetBz());
4202 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[i]);
4203 track0->GetDistance(track1,0.5*(xm[i0]+xm[i1]),dist,AliTracker::GetBz());
4204 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[i]);
4205 for (Int_t i=0;i<3;i++) mdist[i]*=0.33333;
4209 for (Int_t icl=0; icl<160; icl++){
4210 AliTPCclusterMI * cl0 = track0->GetClusterPointer(icl);
4211 AliTPCclusterMI * cl1 = track1->GetClusterPointer(icl);
4214 if (cl0==cl1) sums++;
4218 if (AliTPCReconstructor::StreamLevel()>5) {
4219 TTreeSRedirector &cstream = *fDebugStreamer;
4224 "Tr0.="<<track0<< // seed0
4225 "Tr1.="<<track1<< // seed1
4226 "h0.="<<&helixes[i0]<<
4227 "h1.="<<&helixes[i1]<<
4229 "sum="<<sum<< //the sum of rows with cl in both
4230 "sums="<<sums<< //the sum of shared clusters
4231 "xm0="<<xm[i0]<< // the center of track
4232 "xm1="<<xm[i1]<< // the x center of track
4233 // General cut variables
4234 "dfi="<<dfi<< // distance in fi angle
4235 "dtheta="<<dtheta<< // distance int theta angle
4241 "dist0="<<dist[0]<< //distance x
4242 "dist1="<<dist[1]<< //distance y
4243 "dist2="<<dist[2]<< //distance z
4244 "mdist0="<<mdist[0]<< //distance x
4245 "mdist1="<<mdist[1]<< //distance y
4246 "mdist2="<<mdist[2]<< //distance z
4260 if (AliTPCReconstructor::StreamLevel()>1) {
4261 AliInfo("Time for curling tracks removal DEBUGGING MC");
4268 void AliTPCtrackerMI::FindSplitted(TObjArray * array, AliESDEvent */*esd*/, Int_t /*iter*/){
4270 // Find Splitted tracks and remove the one with worst quality
4271 // Corresponding debug streamer to tune selections - "Splitted2"
4273 // 0. Sort tracks according quility
4274 // 1. Propagate the tracks to the reference radius
4275 // 2. Double_t loop to select close tracks (only to speed up process)
4276 // 3. Calculate cluster overlap ratio - and remove the track if bigger than a threshold
4277 // 4. Delete temporary parameters
4279 const Double_t xref=GetXrow(63); // reference radius -IROC/OROC boundary
4281 const Double_t kCutP1=10; // delta Z cut 10 cm
4282 const Double_t kCutP2=0.15; // delta snp(fi) cut 0.15
4283 const Double_t kCutP3=0.15; // delta tgl(theta) cut 0.15
4284 const Double_t kCutAlpha=0.15; // delta alpha cut
4285 Int_t firstpoint = 0;
4286 Int_t lastpoint = 160;
4288 Int_t nentries = array->GetEntriesFast();
4289 AliExternalTrackParam *params = new AliExternalTrackParam[nentries];
4295 //0. Sort tracks according quality
4296 //1. Propagate the ext. param to reference radius
4297 Int_t nseed = array->GetEntriesFast();
4298 if (nseed<=0) return;
4299 Float_t * quality = new Float_t[nseed];
4300 Int_t * indexes = new Int_t[nseed];
4301 for (Int_t i=0; i<nseed; i++) {
4302 AliTPCseed *pt=(AliTPCseed*)array->UncheckedAt(i);
4307 pt->UpdatePoints(); //select first last max dens points
4308 Float_t * points = pt->GetPoints();
4309 if (points[3]<0.8) quality[i] =-1;
4310 quality[i] = (points[2]-points[0])+pt->GetNumberOfClusters();
4311 //prefer high momenta tracks if overlaps
4312 quality[i] *= TMath::Sqrt(TMath::Abs(pt->Pt())+0.5);
4314 AliTracker::PropagateTrackToBxByBz(&(params[i]),xref,pt->GetMass(),5.,kTRUE);
4315 AliTracker::PropagateTrackToBxByBz(&(params[i]),xref,pt->GetMass(),1.,kTRUE);
4317 TMath::Sort(nseed,quality,indexes);
4319 // 3. Loop over pair of tracks
4321 for (Int_t i0=0; i0<nseed; i0++) {
4322 Int_t index0=indexes[i0];
4323 if (!(array->UncheckedAt(index0))) continue;
4324 AliTPCseed *s1 = (AliTPCseed*)array->UncheckedAt(index0);
4325 if (!s1->IsActive()) continue;
4326 AliExternalTrackParam &par0=params[index0];
4327 for (Int_t i1=i0+1; i1<nseed; i1++) {
4328 Int_t index1=indexes[i1];
4329 if (!(array->UncheckedAt(index1))) continue;
4330 AliTPCseed *s2 = (AliTPCseed*)array->UncheckedAt(index1);
4331 if (!s2->IsActive()) continue;
4332 if (s2->GetKinkIndexes()[0]!=0)
4333 if (s2->GetKinkIndexes()[0] == -s1->GetKinkIndexes()[0]) continue;
4334 AliExternalTrackParam &par1=params[index1];
4335 if (TMath::Abs(par0.GetParameter()[3]-par1.GetParameter()[3])>kCutP3) continue;
4336 if (TMath::Abs(par0.GetParameter()[1]-par1.GetParameter()[1])>kCutP1) continue;
4337 if (TMath::Abs(par0.GetParameter()[2]-par1.GetParameter()[2])>kCutP2) continue;
4338 Double_t dAlpha= TMath::Abs(par0.GetAlpha()-par1.GetAlpha());
4339 if (dAlpha>TMath::Pi()) dAlpha-=TMath::Pi();
4340 if (TMath::Abs(dAlpha)>kCutAlpha) continue;
4345 Int_t firstShared=lastpoint, lastShared=firstpoint;
4346 Int_t firstRow=lastpoint, lastRow=firstpoint;
4348 for (Int_t i=firstpoint;i<lastpoint;i++){
4349 if (s1->GetClusterIndex2(i)>0) nall0++;
4350 if (s2->GetClusterIndex2(i)>0) nall1++;
4351 if (s1->GetClusterIndex2(i)>0 && s2->GetClusterIndex2(i)>0) {
4352 if (i<firstRow) firstRow=i;
4353 if (i>lastRow) lastRow=i;
4355 if ( (s1->GetClusterIndex2(i))==(s2->GetClusterIndex2(i)) && s1->GetClusterIndex2(i)>0) {
4356 if (i<firstShared) firstShared=i;
4357 if (i>lastShared) lastShared=i;
4361 Double_t ratio0 = Float_t(sumShared)/Float_t(TMath::Min(nall0+1,nall1+1));
4362 Double_t ratio1 = Float_t(sumShared)/Float_t(TMath::Max(nall0+1,nall1+1));
4364 if( AliTPCReconstructor::StreamLevel()>1){
4365 TTreeSRedirector &cstream = *fDebugStreamer;
4366 Int_t n0=s1->GetNumberOfClusters();
4367 Int_t n1=s2->GetNumberOfClusters();
4368 Int_t n0F=s1->GetNFoundable();
4369 Int_t n1F=s2->GetNFoundable();
4370 Int_t lab0=s1->GetLabel();
4371 Int_t lab1=s2->GetLabel();
4373 cstream<<"Splitted2"<<
4374 "iter="<<fIteration<<
4375 "lab0="<<lab0<< // MC label if exist
4376 "lab1="<<lab1<< // MC label if exist
4379 "ratio0="<<ratio0<< // shared ratio
4380 "ratio1="<<ratio1<< // shared ratio
4381 "p0.="<<&par0<< // track parameters
4383 "s0.="<<s1<< // full seed
4385 "n0="<<n0<< // number of clusters track 0
4386 "n1="<<n1<< // number of clusters track 1
4387 "nall0="<<nall0<< // number of clusters track 0
4388 "nall1="<<nall1<< // number of clusters track 1
4389 "n0F="<<n0F<< // number of findable
4390 "n1F="<<n1F<< // number of findable
4391 "shared="<<sumShared<< // number of shared clusters
4392 "firstS="<<firstShared<< // first and the last shared row
4393 "lastS="<<lastShared<<
4394 "firstRow="<<firstRow<< // first and the last row with cluster
4395 "lastRow="<<lastRow<< //
4399 // remove track with lower quality
4401 if (ratio0>AliTPCReconstructor::GetRecoParam()->GetCutSharedClusters(0) ||
4402 ratio1>AliTPCReconstructor::GetRecoParam()->GetCutSharedClusters(1)){
4406 delete array->RemoveAt(index1);
4411 // 4. Delete temporary array
4421 void AliTPCtrackerMI::FindCurling(const TObjArray * array, AliESDEvent */*esd*/, Int_t iter)
4424 // find Curling tracks
4425 // Use AliTPCReconstructor::StreamLevel()>1 if you want to tune parameters - cuts
4428 // Algorithm done in 2 phases - because of CPU consumption
4429 // it is n^2 algorithm - for lead-lead 20000x20000 combination are investigated
4430 // see detal in MC part what can be used to cut
4434 const Float_t kMaxC = 400; // maximal curvature to of the track
4435 const Float_t kMaxdTheta = 0.15; // maximal distance in theta
4436 const Float_t kMaxdPhi = 0.15; // maximal distance in phi
4437 const Float_t kPtRatio = 0.3; // ratio between pt
4438 const Float_t kMinDCAR = 2.; // distance to the primary vertex in r - see cpipe cut
4441 // Curling tracks cuts
4444 const Float_t kMaxDeltaRMax = 40; // distance in outer radius
4445 const Float_t kMaxDeltaRMin = 5.; // distance in lower radius - see cpipe cut
4446 const Float_t kMinAngle = 2.9; // angle between tracks
4447 const Float_t kMaxDist = 5; // biggest distance
4449 // The cuts can be tuned using the "MC information stored in Multi tree ==> see FindMultiMC
4452 TCut csign("csign","Tr0.fP[4]*Tr1.fP[4]<0"); //opposite sign
4453 TCut cmax("cmax","abs(Tr0.GetC())>1/400");
4454 TCut cda("cda","sqrt(dtheta^2+dfi^2)<0.15");
4455 TCut ccratio("ccratio","abs((Tr0.fP[4]+Tr1.fP[4])/(abs(Tr0.fP[4])+abs(Tr1.fP[4])))<0.3");
4456 TCut cpipe("cpipe", "min(abs(r0-rc0),abs(r1-rc1))>5");
4458 TCut cdrmax("cdrmax","abs(abs(rc0+r0)-abs(rc1+r1))<40")
4459 TCut cdrmin("cdrmin","abs(abs(rc0+r0)-abs(rc1+r1))<10")
4461 Multi->Draw("dfi","iter==0"+csign+cmax+cda+ccratio); ~94% of curling tracks fulfill
4462 Multi->Draw("min(abs(r0-rc0),abs(r1-rc1))","iter==0&&abs(lab1)==abs(lab0)"+csign+cmax+cda+ccratio+cpipe+cdrmin+cdrmax); //80%
4464 Curling2->Draw("dfi","iter==0&&abs(lab0)==abs(lab1)"+csign+cmax+cdtheta+cdfi+ccratio)
4470 Int_t nentries = array->GetEntriesFast();
4471 AliHelix *helixes = new AliHelix[nentries];
4472 for (Int_t i=0;i<nentries;i++){
4473 AliTPCseed* track = (AliTPCseed*)array->At(i);
4474 if (!track) continue;
4475 track->SetCircular(0);
4476 new (&helixes[i]) AliHelix(*track);
4482 Double_t phase[2][2],radius[2];
4487 for (Int_t i0=0;i0<nentries;i0++){
4488 AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
4489 if (!track0) continue;
4490 if (TMath::Abs(track0->GetC())<1/kMaxC) continue;
4491 Float_t xc0 = helixes[i0].GetHelix(6);
4492 Float_t yc0 = helixes[i0].GetHelix(7);
4493 Float_t r0 = helixes[i0].GetHelix(8);
4494 Float_t rc0 = TMath::Sqrt(xc0*xc0+yc0*yc0);
4495 Float_t fi0 = TMath::ATan2(yc0,xc0);
4497 for (Int_t i1=i0+1;i1<nentries;i1++){
4498 AliTPCseed * track1 = (AliTPCseed*)array->At(i1);
4499 if (!track1) continue;
4500 if (TMath::Abs(track1->GetC())<1/kMaxC) continue;
4501 Float_t xc1 = helixes[i1].GetHelix(6);
4502 Float_t yc1 = helixes[i1].GetHelix(7);
4503 Float_t r1 = helixes[i1].GetHelix(8);
4504 Float_t rc1 = TMath::Sqrt(xc1*xc1+yc1*yc1);
4505 Float_t fi1 = TMath::ATan2(yc1,xc1);
4507 Float_t dfi = fi0-fi1;
4510 if (dfi>1.5*TMath::Pi()) dfi-=TMath::Pi(); // take care about edge effect
4511 if (dfi<-1.5*TMath::Pi()) dfi+=TMath::Pi(); //
4512 Float_t dtheta = TMath::Abs(track0->GetTgl()-track1->GetTgl())<TMath::Abs(track0->GetTgl()+track1->GetTgl())? track0->GetTgl()-track1->GetTgl():track0->GetTgl()+track1->GetTgl();
4516 if (track0->GetBConstrain()&&track1->GetBConstrain()) continue; // not constrained
4517 if (track1->GetSigned1Pt()*track0->GetSigned1Pt()>0) continue; // not the same sign
4518 if ( TMath::Abs(track1->GetTgl()+track0->GetTgl())>kMaxdTheta) continue; //distance in the Theta
4519 if ( TMath::Abs(dfi)>kMaxdPhi) continue; //distance in phi
4520 if ( TMath::Sqrt(dfi*dfi+dtheta*dtheta)>kMaxdPhi) continue; //common angular offset
4522 Float_t pt0 = track0->GetSignedPt();
4523 Float_t pt1 = track1->GetSignedPt();
4524 if ((TMath::Abs(pt0+pt1)/(TMath::Abs(pt0)+TMath::Abs(pt1)))>kPtRatio) continue;
4525 if ((iter==1) && TMath::Abs(TMath::Abs(rc0+r0)-TMath::Abs(rc1+r1))>kMaxDeltaRMax) continue;
4526 if ((iter!=1) &&TMath::Abs(TMath::Abs(rc0-r0)-TMath::Abs(rc1-r1))>kMaxDeltaRMin) continue;
4527 if (TMath::Min(TMath::Abs(rc0-r0),TMath::Abs(rc1-r1))<kMinDCAR) continue;
4530 // Now find closest approach
4534 Int_t npoints = helixes[i0].GetRPHIintersections(helixes[i1], phase, radius,10);
4535 if (npoints==0) continue;
4536 helixes[i0].GetClosestPhases(helixes[i1], phase);
4540 Double_t hangles[3];
4541 helixes[i0].Evaluate(phase[0][0],xyz0);
4542 helixes[i1].Evaluate(phase[0][1],xyz1);
4544 helixes[i0].GetAngle(phase[0][0],helixes[i1],phase[0][1],hangles);
4545 Double_t deltah[2],deltabest;
4546 if (TMath::Abs(hangles[2])<kMinAngle) continue;
4550 helixes[i0].ParabolicDCA(helixes[i1],phase[0][0],phase[0][1],radius[0],deltah[0],2);
4552 helixes[i0].ParabolicDCA(helixes[i1],phase[1][0],phase[1][1],radius[1],deltah[1],2);
4553 if (deltah[1]<deltah[0]) ibest=1;
4555 deltabest = TMath::Sqrt(deltah[ibest]);
4556 helixes[i0].Evaluate(phase[ibest][0],xyz0);
4557 helixes[i1].Evaluate(phase[ibest][1],xyz1);
4558 helixes[i0].GetAngle(phase[ibest][0],helixes[i1],phase[ibest][1],hangles);
4559 Double_t radiusbest = TMath::Sqrt(radius[ibest]);
4561 if (deltabest>kMaxDist) continue;
4562 // if (mindcar+mindcaz<40 && (TMath::Abs(hangles[2])<kMinAngle ||deltabest>3)) continue;
4563 Bool_t sign =kFALSE;
4564 if (hangles[2]>kMinAngle) sign =kTRUE;
4567 // circular[i0] = kTRUE;
4568 // circular[i1] = kTRUE;
4569 if (track0->OneOverPt()<track1->OneOverPt()){
4570 track0->SetCircular(track0->GetCircular()+1);
4571 track1->SetCircular(track1->GetCircular()+2);
4574 track1->SetCircular(track1->GetCircular()+1);
4575 track0->SetCircular(track0->GetCircular()+2);
4578 if (AliTPCReconstructor::StreamLevel()>2){
4580 //debug stream to tune "fine" cuts
4581 Int_t lab0=track0->GetLabel();
4582 Int_t lab1=track1->GetLabel();
4583 TTreeSRedirector &cstream = *fDebugStreamer;
4584 cstream<<"Curling2"<<
4600 "npoints="<<npoints<<
4601 "hangles0="<<hangles[0]<<
4602 "hangles1="<<hangles[1]<<
4603 "hangles2="<<hangles[2]<<
4606 "radius="<<radiusbest<<
4607 "deltabest="<<deltabest<<
4608 "phase0="<<phase[ibest][0]<<
4609 "phase1="<<phase[ibest][1]<<
4617 if (AliTPCReconstructor::StreamLevel()>1) {
4618 AliInfo("Time for curling tracks removal");
4627 void AliTPCtrackerMI::FindKinks(TObjArray * array, AliESDEvent *esd)
4634 TObjArray *kinks= new TObjArray(10000);
4635 // TObjArray *v0s= new TObjArray(10000);
4636 Int_t nentries = array->GetEntriesFast();
4637 AliHelix *helixes = new AliHelix[nentries];
4638 Int_t *sign = new Int_t[nentries];
4639 Int_t *nclusters = new Int_t[nentries];
4640 Float_t *alpha = new Float_t[nentries];
4641 AliKink *kink = new AliKink();
4642 Int_t * usage = new Int_t[nentries];
4643 Float_t *zm = new Float_t[nentries];
4644 Float_t *z0 = new Float_t[nentries];
4645 Float_t *fim = new Float_t[nentries];
4646 Float_t *shared = new Float_t[nentries];
4647 Bool_t *circular = new Bool_t[nentries];
4648 Float_t *dca = new Float_t[nentries];
4649 //const AliESDVertex * primvertex = esd->GetVertex();
4651 // nentries = array->GetEntriesFast();
4656 for (Int_t i=0;i<nentries;i++){
4659 AliTPCseed* track = (AliTPCseed*)array->At(i);
4660 if (!track) continue;
4661 track->SetCircular(0);
4663 track->UpdatePoints();
4664 if (( track->GetPoints()[2]- track->GetPoints()[0])>5 && track->GetPoints()[3]>0.8){
4666 nclusters[i]=track->GetNumberOfClusters();
4667 alpha[i] = track->GetAlpha();
4668 new (&helixes[i]) AliHelix(*track);
4670 helixes[i].Evaluate(0,xyz);
4671 sign[i] = (track->GetC()>0) ? -1:1;
4674 if (track->GetProlongation(x,y,z)){
4676 fim[i] = alpha[i]+TMath::ATan2(y,x);
4679 zm[i] = track->GetZ();
4683 circular[i]= kFALSE;
4684 if (track->GetProlongation(0,y,z)) z0[i] = z;
4685 dca[i] = track->GetD(0,0);
4691 Int_t ncandidates =0;
4694 Double_t phase[2][2],radius[2];
4697 // Find circling track
4699 for (Int_t i0=0;i0<nentries;i0++){
4700 AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
4701 if (!track0) continue;
4702 if (track0->GetNumberOfClusters()<40) continue;
4703 if (TMath::Abs(1./track0->GetC())>200) continue;
4704 for (Int_t i1=i0+1;i1<nentries;i1++){
4705 AliTPCseed * track1 = (AliTPCseed*)array->At(i1);
4706 if (!track1) continue;
4707 if (track1->GetNumberOfClusters()<40) continue;
4708 if ( TMath::Abs(track1->GetTgl()+track0->GetTgl())>0.1) continue;
4709 if (track0->GetBConstrain()&&track1->GetBConstrain()) continue;
4710 if (TMath::Abs(1./track1->GetC())>200) continue;
4711 if (track1->GetSigned1Pt()*track0->GetSigned1Pt()>0) continue;
4712 if (track1->GetTgl()*track0->GetTgl()>0) continue;
4713 if (TMath::Max(TMath::Abs(1./track0->GetC()),TMath::Abs(1./track1->GetC()))>190) continue;
4714 if (track0->GetBConstrain()&&track1->OneOverPt()<track0->OneOverPt()) continue; //returning - lower momenta
4715 if (track1->GetBConstrain()&&track0->OneOverPt()<track1->OneOverPt()) continue; //returning - lower momenta
4717 Float_t mindcar = TMath::Min(TMath::Abs(dca[i0]),TMath::Abs(dca[i1]));
4718 if (mindcar<5) continue;
4719 Float_t mindcaz = TMath::Min(TMath::Abs(z0[i0]-GetZ()),TMath::Abs(z0[i1]-GetZ()));
4720 if (mindcaz<5) continue;
4721 if (mindcar+mindcaz<20) continue;
4724 Float_t xc0 = helixes[i0].GetHelix(6);
4725 Float_t yc0 = helixes[i0].GetHelix(7);
4726 Float_t r0 = helixes[i0].GetHelix(8);
4727 Float_t xc1 = helixes[i1].GetHelix(6);
4728 Float_t yc1 = helixes[i1].GetHelix(7);
4729 Float_t r1 = helixes[i1].GetHelix(8);
4731 Float_t rmean = (r0+r1)*0.5;
4732 Float_t delta =TMath::Sqrt((xc1-xc0)*(xc1-xc0)+(yc1-yc0)*(yc1-yc0));
4733 //if (delta>30) continue;
4734 if (delta>rmean*0.25) continue;
4735 if (TMath::Abs(r0-r1)/rmean>0.3) continue;
4737 Int_t npoints = helixes[i0].GetRPHIintersections(helixes[i1], phase, radius,10);
4738 if (npoints==0) continue;
4739 helixes[i0].GetClosestPhases(helixes[i1], phase);
4743 Double_t hangles[3];
4744 helixes[i0].Evaluate(phase[0][0],xyz0);
4745 helixes[i1].Evaluate(phase[0][1],xyz1);
4747 helixes[i0].GetAngle(phase[0][0],helixes[i1],phase[0][1],hangles);
4748 Double_t deltah[2],deltabest;
4749 if (hangles[2]<2.8) continue;
4752 helixes[i0].ParabolicDCA(helixes[i1],phase[0][0],phase[0][1],radius[0],deltah[0],2);
4754 helixes[i0].ParabolicDCA(helixes[i1],phase[1][0],phase[1][1],radius[1],deltah[1],2);
4755 if (deltah[1]<deltah[0]) ibest=1;
4757 deltabest = TMath::Sqrt(deltah[ibest]);
4758 helixes[i0].Evaluate(phase[ibest][0],xyz0);
4759 helixes[i1].Evaluate(phase[ibest][1],xyz1);
4760 helixes[i0].GetAngle(phase[ibest][0],helixes[i1],phase[ibest][1],hangles);
4761 Double_t radiusbest = TMath::Sqrt(radius[ibest]);
4763 if (deltabest>6) continue;
4764 if (mindcar+mindcaz<40 && (hangles[2]<3.12||deltabest>3)) continue;
4765 Bool_t lsign =kFALSE;
4766 if (hangles[2]>3.06) lsign =kTRUE;
4769 circular[i0] = kTRUE;
4770 circular[i1] = kTRUE;
4771 if (track0->OneOverPt()<track1->OneOverPt()){
4772 track0->SetCircular(track0->GetCircular()+1);
4773 track1->SetCircular(track1->GetCircular()+2);
4776 track1->SetCircular(track1->GetCircular()+1);
4777 track0->SetCircular(track0->GetCircular()+2);
4780 if (lsign&&AliTPCReconstructor::StreamLevel()>1){
4782 Int_t lab0=track0->GetLabel();
4783 Int_t lab1=track1->GetLabel();
4784 TTreeSRedirector &cstream = *fDebugStreamer;
4785 cstream<<"Curling"<<
4792 "mindcar="<<mindcar<<
4793 "mindcaz="<<mindcaz<<
4796 "npoints="<<npoints<<
4797 "hangles0="<<hangles[0]<<
4798 "hangles2="<<hangles[2]<<
4803 "radius="<<radiusbest<<
4804 "deltabest="<<deltabest<<
4805 "phase0="<<phase[ibest][0]<<
4806 "phase1="<<phase[ibest][1]<<
4816 for (Int_t i =0;i<nentries;i++){
4817 if (sign[i]==0) continue;
4818 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
4825 Double_t cradius0 = 40*40;
4826 Double_t cradius1 = 270*270;
4829 Double_t cdist3=0.55;
4830 for (Int_t j =i+1;j<nentries;j++){
4832 if (sign[j]*sign[i]<1) continue;
4833 if ( (nclusters[i]+nclusters[j])>200) continue;
4834 if ( (nclusters[i]+nclusters[j])<80) continue;
4835 if ( TMath::Abs(zm[i]-zm[j])>60.) continue;
4836 if ( TMath::Abs(fim[i]-fim[j])>0.6 && TMath::Abs(fim[i]-fim[j])<5.7 ) continue;
4837 //AliTPCseed * track1 = (AliTPCseed*)array->At(j); Double_t phase[2][2],radius[2];
4838 Int_t npoints = helixes[i].GetRPHIintersections(helixes[j], phase, radius,20);
4839 if (npoints<1) continue;
4842 if (radius[0]<cradius0||radius[0]>cradius1) continue;
4845 if ( (radius[0]<cradius0||radius[0]>cradius1) && (radius[1]<cradius0||radius[1]>cradius1) ) continue;
4848 Double_t delta1=10000,delta2=10000;
4849 // cuts on the intersection radius
4850 helixes[i].LinearDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta1);
4851 if (radius[0]<20&&delta1<1) continue; //intersection at vertex
4852 if (radius[0]<10&&delta1<3) continue; //intersection at vertex
4854 helixes[i].LinearDCA(helixes[j],phase[1][0],phase[1][1],radius[1],delta2);
4855 if (radius[1]<20&&delta2<1) continue; //intersection at vertex
4856 if (radius[1]<10&&delta2<3) continue; //intersection at vertex
4859 Double_t distance1 = TMath::Min(delta1,delta2);
4860 if (distance1>cdist1) continue; // cut on DCA linear approximation
4862 npoints = helixes[i].GetRPHIintersections(helixes[j], phase, radius,20);
4863 helixes[i].ParabolicDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta1);
4864 if (radius[0]<20&&delta1<1) continue; //intersection at vertex
4865 if (radius[0]<10&&delta1<3) continue; //intersection at vertex
4868 helixes[i].ParabolicDCA(helixes[j],phase[1][0],phase[1][1],radius[1],delta2);
4869 if (radius[1]<20&&delta2<1) continue; //intersection at vertex
4870 if (radius[1]<10&&delta2<3) continue; //intersection at vertex
4872 distance1 = TMath::Min(delta1,delta2);
4875 rkink = TMath::Sqrt(radius[0]);
4878 rkink = TMath::Sqrt(radius[1]);
4880 if (distance1>cdist2) continue;
4883 AliTPCseed * track1 = (AliTPCseed*)array->At(j);
4886 Int_t row0 = GetRowNumber(rkink);
4887 if (row0<10) continue;
4888 if (row0>150) continue;
4891 Float_t dens00=-1,dens01=-1;
4892 Float_t dens10=-1,dens11=-1;
4894 Int_t found,foundable,ishared;
4895 track0->GetClusterStatistic(0,row0-5, found, foundable,ishared,kFALSE);
4896 if (foundable>5) dens00 = Float_t(found)/Float_t(foundable);
4897 track0->GetClusterStatistic(row0+5,155, found, foundable,ishared,kFALSE);
4898 if (foundable>5) dens01 = Float_t(found)/Float_t(foundable);
4900 track1->GetClusterStatistic(0,row0-5, found, foundable,ishared,kFALSE);
4901 if (foundable>10) dens10 = Float_t(found)/Float_t(foundable);
4902 track1->GetClusterStatistic(row0+5,155, found, foundable,ishared,kFALSE);
4903 if (foundable>10) dens11 = Float_t(found)/Float_t(foundable);
4905 if (dens00<dens10 && dens01<dens11) continue;
4906 if (dens00>dens10 && dens01>dens11) continue;
4907 if (TMath::Max(dens00,dens10)<0.1) continue;
4908 if (TMath::Max(dens01,dens11)<0.3) continue;
4910 if (TMath::Min(dens00,dens10)>0.6) continue;
4911 if (TMath::Min(dens01,dens11)>0.6) continue;
4914 AliTPCseed * ktrack0, *ktrack1;
4923 if (TMath::Abs(ktrack0->GetC())>5) continue; // cut on the curvature for mother particle
4924 AliExternalTrackParam paramm(*ktrack0);
4925 AliExternalTrackParam paramd(*ktrack1);
4926 if (row0>60&&ktrack1->GetReference().GetX()>90.)new (¶md) AliExternalTrackParam(ktrack1->GetReference());
4929 kink->SetMother(paramm);
4930 kink->SetDaughter(paramd);
4933 Float_t x[3] = { kink->GetPosition()[0],kink->GetPosition()[1],kink->GetPosition()[2]};
4935 fkParam->Transform0to1(x,index);
4936 fkParam->Transform1to2(x,index);
4937 row0 = GetRowNumber(x[0]);
4939 if (kink->GetR()<100) continue;
4940 if (kink->GetR()>240) continue;
4941 if (kink->GetPosition()[2]/kink->GetR()>AliTPCReconstructor::GetCtgRange()) continue; //out of fiducial volume
4942 if (kink->GetDistance()>cdist3) continue;
4943 Float_t dird = kink->GetDaughterP()[0]*kink->GetPosition()[0]+kink->GetDaughterP()[1]*kink->GetPosition()[1]; // rough direction estimate
4944 if (dird<0) continue;
4946 Float_t dirm = kink->GetMotherP()[0]*kink->GetPosition()[0]+kink->GetMotherP()[1]*kink->GetPosition()[1]; // rough direction estimate
4947 if (dirm<0) continue;
4948 Float_t mpt = TMath::Sqrt(kink->GetMotherP()[0]*kink->GetMotherP()[0]+kink->GetMotherP()[1]*kink->GetMotherP()[1]);
4949 if (mpt<0.2) continue;
4952 //for high momenta momentum not defined well in first iteration
4953 Double_t qt = TMath::Sin(kink->GetAngle(2))*ktrack1->GetP();
4954 if (qt>0.35) continue;
4957 kink->SetLabel(CookLabel(ktrack0,0.4,0,row0),0);
4958 kink->SetLabel(CookLabel(ktrack1,0.4,row0,160),1);
4960 kink->SetTPCDensity(dens00,0,0);
4961 kink->SetTPCDensity(dens01,0,1);
4962 kink->SetTPCDensity(dens10,1,0);
4963 kink->SetTPCDensity(dens11,1,1);
4964 kink->SetIndex(i,0);
4965 kink->SetIndex(j,1);
4968 kink->SetTPCDensity(dens10,0,0);
4969 kink->SetTPCDensity(dens11,0,1);
4970 kink->SetTPCDensity(dens00,1,0);
4971 kink->SetTPCDensity(dens01,1,1);
4972 kink->SetIndex(j,0);
4973 kink->SetIndex(i,1);
4976 if (mpt<1||kink->GetAngle(2)>0.1){
4977 // angle and densities not defined yet
4978 if (kink->GetTPCDensityFactor()<0.8) continue;
4979 if ((2-kink->GetTPCDensityFactor())*kink->GetDistance() >0.25) continue;
4980 if (kink->GetAngle(2)*ktrack0->GetP()<0.003) continue; //too small angle
4981 if (kink->GetAngle(2)>0.2&&kink->GetTPCDensityFactor()<1.15) continue;
4982 if (kink->GetAngle(2)>0.2&&kink->GetTPCDensity(0,1)>0.05) continue;
4984 Float_t criticalangle = track0->GetSigmaSnp2()+track0->GetSigmaTgl2();
4985 criticalangle+= track1->GetSigmaSnp2()+track1->GetSigmaTgl2();
4986 criticalangle= 3*TMath::Sqrt(criticalangle);
4987 if (criticalangle>0.02) criticalangle=0.02;
4988 if (kink->GetAngle(2)<criticalangle) continue;
4991 Int_t drow = Int_t(2.+0.5/(0.05+kink->GetAngle(2))); // overlap region defined
4992 Float_t shapesum =0;
4994 for ( Int_t row = row0-drow; row<row0+drow;row++){
4995 if (row<0) continue;
4996 if (row>155) continue;
4997 if (ktrack0->GetClusterPointer(row)){
4998 AliTPCTrackerPoint *point =ktrack0->GetTrackPoint(row);
4999 shapesum+=point->GetSigmaY()+point->GetSigmaZ();
5002 if (ktrack1->GetClusterPointer(row)){
5003 AliTPCTrackerPoint *point =ktrack1->GetTrackPoint(row);
5004 shapesum+=point->GetSigmaY()+point->GetSigmaZ();
5009 kink->SetShapeFactor(-1.);
5012 kink->SetShapeFactor(shapesum/sum);
5014 // esd->AddKink(kink);
5016 // kink->SetMother(paramm);
5017 //kink->SetDaughter(paramd);
5019 Double_t chi2P2 = paramm.GetParameter()[2]-paramd.GetParameter()[2];
5021 chi2P2/=paramm.GetCovariance()[5]+paramd.GetCovariance()[5];
5022 Double_t chi2P3 = paramm.GetParameter()[3]-paramd.GetParameter()[3];
5024 chi2P3/=paramm.GetCovariance()[9]+paramd.GetCovariance()[9];
5026 if (AliTPCReconstructor::StreamLevel()>1) {
5027 (*fDebugStreamer)<<"kinkLpt"<<
5035 if ( chi2P2+chi2P3<AliTPCReconstructor::GetRecoParam()->GetKinkAngleCutChi2(0)){
5039 kinks->AddLast(kink);
5045 // sort the kinks according quality - and refit them towards vertex
5047 Int_t nkinks = kinks->GetEntriesFast();
5048 Float_t *quality = new Float_t[nkinks];
5049 Int_t *indexes = new Int_t[nkinks];
5050 AliTPCseed *mothers = new AliTPCseed[nkinks];
5051 AliTPCseed *daughters = new AliTPCseed[nkinks];
5054 for (Int_t i=0;i<nkinks;i++){
5056 AliKink *kinkl = (AliKink*)kinks->At(i);
5058 // refit kinks towards vertex
5060 Int_t index0 = kinkl->GetIndex(0);
5061 Int_t index1 = kinkl->GetIndex(1);
5062 AliTPCseed * ktrack0 = (AliTPCseed*)array->At(index0);
5063 AliTPCseed * ktrack1 = (AliTPCseed*)array->At(index1);
5065 Int_t sumn=ktrack0->GetNumberOfClusters()+ktrack1->GetNumberOfClusters();
5067 // Refit Kink under if too small angle
5069 if (kinkl->GetAngle(2)<0.05){
5070 kinkl->SetTPCRow0(GetRowNumber(kinkl->GetR()));
5071 Int_t row0 = kinkl->GetTPCRow0();
5072 Int_t drow = Int_t(2.+0.5/(0.05+kinkl->GetAngle(2)));
5075 Int_t last = row0-drow;
5076 if (last<40) last=40;
5077 if (last<ktrack0->GetFirstPoint()+25) last = ktrack0->GetFirstPoint()+25;
5078 AliTPCseed* seed0 = ReSeed(ktrack0,last,kFALSE);
5081 Int_t first = row0+drow;
5082 if (first>130) first=130;
5083 if (first>ktrack1->GetLastPoint()-25) first = TMath::Max(ktrack1->GetLastPoint()-25,30);
5084 AliTPCseed* seed1 = ReSeed(ktrack1,first,kTRUE);
5086 if (seed0 && seed1){
5087 kinkl->SetStatus(1,8);
5088 if (RefitKink(*seed0,*seed1,*kinkl)) kinkl->SetStatus(1,9);
5089 row0 = GetRowNumber(kinkl->GetR());
5090 sumn = seed0->GetNumberOfClusters()+seed1->GetNumberOfClusters();
5091 mothers[i] = *seed0;
5092 daughters[i] = *seed1;
5095 delete kinks->RemoveAt(i);
5096 if (seed0) delete seed0;
5097 if (seed1) delete seed1;
5100 if (kinkl->GetDistance()>0.5 || kinkl->GetR()<110 || kinkl->GetR()>240) {
5101 delete kinks->RemoveAt(i);
5102 if (seed0) delete seed0;
5103 if (seed1) delete seed1;
5111 if (kinkl) quality[i] = 160*((0.1+kinkl->GetDistance())*(2.-kinkl->GetTPCDensityFactor()))/(sumn+40.); //the longest -clossest will win
5113 TMath::Sort(nkinks,quality,indexes,kFALSE);
5115 //remove double find kinks
5117 for (Int_t ikink0=1;ikink0<nkinks;ikink0++){
5118 AliKink * kink0 = (AliKink*) kinks->At(indexes[ikink0]);
5119 if (!kink0) continue;
5121 for (Int_t ikink1=0;ikink1<ikink0;ikink1++){
5122 if (!kink0) continue;
5123 AliKink * kink1 = (AliKink*) kinks->At(indexes[ikink1]);
5124 if (!kink1) continue;
5125 // if not close kink continue
5126 if (TMath::Abs(kink1->GetPosition()[2]-kink0->GetPosition()[2])>10) continue;
5127 if (TMath::Abs(kink1->GetPosition()[1]-kink0->GetPosition()[1])>10) continue;
5128 if (TMath::Abs(kink1->GetPosition()[0]-kink0->GetPosition()[0])>10) continue;
5130 AliTPCseed &mother0 = mothers[indexes[ikink0]];
5131 AliTPCseed &daughter0 = daughters[indexes[ikink0]];
5132 AliTPCseed &mother1 = mothers[indexes[ikink1]];
5133 AliTPCseed &daughter1 = daughters[indexes[ikink1]];
5134 Int_t row0 = (kink0->GetTPCRow0()+kink1->GetTPCRow0())/2;
5143 for (Int_t i=0;i<row0;i++){
5144 if (mother0.GetClusterIndex(i)>0 && mother1.GetClusterIndex(i)>0){
5147 if (mother0.GetClusterIndex(i)==mother1.GetClusterIndex(i)){
5154 for (Int_t i=row0;i<158;i++){
5155 if (daughter0.GetClusterIndex(i)>0 && daughter0.GetClusterIndex(i)>0){
5158 if (mother0.GetClusterIndex(i)==mother1.GetClusterIndex(i)){
5164 Float_t ratio = Float_t(same+1)/Float_t(both+1);
5165 Float_t ratiom = Float_t(samem+1)/Float_t(bothm+1);
5166 Float_t ratiod = Float_t(samed+1)/Float_t(bothd+1);
5167 if (ratio>0.3 && ratiom>0.5 &&ratiod>0.5) {
5168 Int_t sum0 = mother0.GetNumberOfClusters()+daughter0.GetNumberOfClusters();
5169 Int_t sum1 = mother1.GetNumberOfClusters()+daughter1.GetNumberOfClusters();
5171 shared[kink0->GetIndex(0)]= kTRUE;
5172 shared[kink0->GetIndex(1)]= kTRUE;
5173 delete kinks->RemoveAt(indexes[ikink0]);
5176 shared[kink1->GetIndex(0)]= kTRUE;
5177 shared[kink1->GetIndex(1)]= kTRUE;
5178 delete kinks->RemoveAt(indexes[ikink1]);
5185 for (Int_t i=0;i<nkinks;i++){
5186 AliKink * kinkl = (AliKink*) kinks->At(indexes[i]);
5187 if (!kinkl) continue;
5188 kinkl->SetTPCRow0(GetRowNumber(kinkl->GetR()));
5189 Int_t index0 = kinkl->GetIndex(0);
5190 Int_t index1 = kinkl->GetIndex(1);
5191 if (circular[index0]||(circular[index1]&&kinkl->GetDistance()>0.2)) continue;
5192 kinkl->SetMultiple(usage[index0],0);
5193 kinkl->SetMultiple(usage[index1],1);
5194 if (kinkl->GetMultiple()[0]+kinkl->GetMultiple()[1]>2) continue;
5195 if (kinkl->GetMultiple()[0]+kinkl->GetMultiple()[1]>0 && quality[indexes[i]]>0.2) continue;
5196 if (kinkl->GetMultiple()[0]+kinkl->GetMultiple()[1]>0 && kinkl->GetDistance()>0.2) continue;
5197 if (circular[index0]||(circular[index1]&&kinkl->GetDistance()>0.1)) continue;
5199 AliTPCseed * ktrack0 = (AliTPCseed*)array->At(index0);
5200 AliTPCseed * ktrack1 = (AliTPCseed*)array->At(index1);
5201 if (!ktrack0 || !ktrack1) continue;
5202 Int_t index = esd->AddKink(kinkl);
5205 if ( ktrack0->GetKinkIndex(0)==0 && ktrack1->GetKinkIndex(0)==0) { //best kink
5206 if (mothers[indexes[i]].GetNumberOfClusters()>20 && daughters[indexes[i]].GetNumberOfClusters()>20 && (mothers[indexes[i]].GetNumberOfClusters()+daughters[indexes[i]].GetNumberOfClusters())>100){
5207 *ktrack0 = mothers[indexes[i]];
5208 *ktrack1 = daughters[indexes[i]];
5212 ktrack0->SetKinkIndex(usage[index0],-(index+1));
5213 ktrack1->SetKinkIndex(usage[index1], (index+1));
5218 // Remove tracks corresponding to shared kink's
5220 for (Int_t i=0;i<nentries;i++){
5221 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5222 if (!track0) continue;
5223 if (track0->GetKinkIndex(0)!=0) continue;
5224 if (shared[i]) delete array->RemoveAt(i);
5229 RemoveUsed2(array,0.5,0.4,30);
5231 for (Int_t i=0;i<nentries;i++){
5232 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5233 if (!track0) continue;
5234 track0->CookdEdx(0.02,0.6);
5238 for (Int_t i=0;i<nentries;i++){
5239 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5240 if (!track0) continue;
5241 if (track0->Pt()<1.4) continue;
5242 //remove double high momenta tracks - overlapped with kink candidates
5245 for (Int_t icl=track0->GetFirstPoint();icl<track0->GetLastPoint(); icl++){
5246 if (track0->GetClusterPointer(icl)!=0){
5248 if (track0->GetClusterPointer(icl)->IsUsed(10)) ishared++;
5251 if (Float_t(ishared+1)/Float_t(all+1)>0.5) {
5252 delete array->RemoveAt(i);
5256 if (track0->GetKinkIndex(0)!=0) continue;
5257 if (track0->GetNumberOfClusters()<80) continue;
5259 AliTPCseed *pmother = new AliTPCseed();
5260 AliTPCseed *pdaughter = new AliTPCseed();
5261 AliKink *pkink = new AliKink;
5263 AliTPCseed & mother = *pmother;
5264 AliTPCseed & daughter = *pdaughter;
5265 AliKink & kinkl = *pkink;
5266 if (CheckKinkPoint(track0,mother,daughter, kinkl)){
5267 if (mother.GetNumberOfClusters()<30||daughter.GetNumberOfClusters()<20) {
5271 continue; //too short tracks
5273 if (mother.Pt()<1.4) {
5279 Int_t row0= kinkl.GetTPCRow0();
5280 if (kinkl.GetDistance()>0.5 || kinkl.GetR()<110. || kinkl.GetR()>240.) {
5287 Int_t index = esd->AddKink(&kinkl);
5288 mother.SetKinkIndex(0,-(index+1));
5289 daughter.SetKinkIndex(0,index+1);
5290 if (mother.GetNumberOfClusters()>50) {
5291 delete array->RemoveAt(i);
5292 array->AddAt(new AliTPCseed(mother),i);
5295 array->AddLast(new AliTPCseed(mother));
5297 array->AddLast(new AliTPCseed(daughter));
5298 for (Int_t icl=0;icl<row0;icl++) {
5299 if (mother.GetClusterPointer(icl)) mother.GetClusterPointer(icl)->Use(20);
5302 for (Int_t icl=row0;icl<158;icl++) {
5303 if (daughter.GetClusterPointer(icl)) daughter.GetClusterPointer(icl)->Use(20);
5312 delete [] daughters;
5334 AliInfo(Form("Ncandidates=\t%d\t%d\t%d\t%d\n",esd->GetNumberOfKinks(),ncandidates,ntracks,nall));
5338 void AliTPCtrackerMI::FindV0s(const TObjArray * array, AliESDEvent *const esd)
5344 TObjArray *tpcv0s = new TObjArray(100000);
5345 Int_t nentries = array->GetEntriesFast();
5346 AliHelix *helixes = new AliHelix[nentries];
5347 Int_t *sign = new Int_t[nentries];
5348 Float_t *alpha = new Float_t[nentries];
5349 Float_t *z0 = new Float_t[nentries];
5350 Float_t *dca = new Float_t[nentries];
5351 Float_t *sdcar = new Float_t[nentries];
5352 Float_t *cdcar = new Float_t[nentries];
5353 Float_t *pulldcar = new Float_t[nentries];
5354 Float_t *pulldcaz = new Float_t[nentries];
5355 Float_t *pulldca = new Float_t[nentries];
5356 Bool_t *isPrim = new Bool_t[nentries];
5357 const AliESDVertex * primvertex = esd->GetVertex();
5358 Double_t zvertex = primvertex->GetZv();
5360 // nentries = array->GetEntriesFast();
5362 for (Int_t i=0;i<nentries;i++){
5365 AliTPCseed* track = (AliTPCseed*)array->At(i);
5366 if (!track) continue;
5367 track->GetV0Indexes()[0] = 0; //rest v0 indexes
5368 track->GetV0Indexes()[1] = 0; //rest v0 indexes
5369 track->GetV0Indexes()[2] = 0; //rest v0 indexes
5371 alpha[i] = track->GetAlpha();
5372 new (&helixes[i]) AliHelix(*track);
5374 helixes[i].Evaluate(0,xyz);
5375 sign[i] = (track->GetC()>0) ? -1:1;
5379 if (track->GetProlongation(0,y,z)) z0[i] = z;
5380 dca[i] = track->GetD(0,0);
5382 // dca error parrameterezation + pulls
5384 sdcar[i] = TMath::Sqrt(0.150*0.150+(100*track->GetC())*(100*track->GetC()));
5385 if (TMath::Abs(track->GetTgl())>1) sdcar[i]*=2.5;
5386 cdcar[i] = TMath::Exp((TMath::Abs(track->GetC())-0.0106)*525.3);
5387 pulldcar[i] = (dca[i]-cdcar[i])/sdcar[i];
5388 pulldcaz[i] = (z0[i]-zvertex)/sdcar[i];
5389 pulldca[i] = TMath::Sqrt(pulldcar[i]*pulldcar[i]+pulldcaz[i]*pulldcaz[i]);
5390 if (track->TPCrPID(1)+track->TPCrPID(2)+track->TPCrPID(3)>0.5) {
5391 if (pulldca[i]<3.) isPrim[i]=kTRUE; //pion, muon and Kaon 3 sigma cut
5393 if (track->TPCrPID(4)>0.5) {
5394 if (pulldca[i]<0.5) isPrim[i]=kTRUE; //proton 0.5 sigma cut
5396 if (track->TPCrPID(0)>0.4) {
5397 isPrim[i]=kFALSE; //electron no sigma cut
5404 Int_t ncandidates =0;
5407 Double_t phase[2][2],radius[2];
5413 Float_t fprimvertex[3]={GetX(),GetY(),GetZ()};
5415 Double_t cradius0 = 10*10;
5416 Double_t cradius1 = 200*200;
5419 Double_t cpointAngle = 0.95;
5421 Double_t delta[2]={10000,10000};
5422 for (Int_t i =0;i<nentries;i++){
5423 if (sign[i]==0) continue;
5424 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5425 if (!track0) continue;
5426 if (AliTPCReconstructor::StreamLevel()>1){
5427 TTreeSRedirector &cstream = *fDebugStreamer;
5432 "zvertex="<<zvertex<<
5433 "sdcar0="<<sdcar[i]<<
5434 "cdcar0="<<cdcar[i]<<
5435 "pulldcar0="<<pulldcar[i]<<
5436 "pulldcaz0="<<pulldcaz[i]<<
5437 "pulldca0="<<pulldca[i]<<
5438 "isPrim="<<isPrim[i]<<
5442 if (track0->GetSigned1Pt()<0) continue;
5443 if (track0->GetKinkIndex(0)>0||isPrim[i]) continue; //daughter kink
5445 if (TMath::Abs(helixes[i].GetHelix(4))<0.000000001) continue;
5450 for (Int_t j =0;j<nentries;j++){
5451 AliTPCseed * track1 = (AliTPCseed*)array->At(j);
5452 if (!track1) continue;
5453 if (track1->GetKinkIndex(0)>0 || isPrim[j]) continue; //daughter kink
5454 if (sign[j]*sign[i]>0) continue;
5455 if (TMath::Abs(helixes[j].GetHelix(4))<0.000001) continue;
5456 if (track0->GetCircular()+track1->GetCircular()>1) continue; //circling -returning track
5459 // DCA to prim vertex cut
5465 Int_t npoints = helixes[i].GetRPHIintersections(helixes[j], phase, radius,cdist2);
5466 if (npoints<1) continue;
5470 if (radius[0]<cradius0||radius[0]>cradius1) continue;
5471 helixes[i].LinearDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta[0]);
5472 if (delta[0]>cdist1) continue;
5475 if (TMath::Max(radius[0],radius[1])<cradius0|| TMath::Min(radius[0],radius[1])>cradius1) continue;
5476 helixes[i].LinearDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta[0]);
5477 helixes[i].LinearDCA(helixes[j],phase[1][0],phase[1][1],radius[1],delta[1]);
5478 if (delta[1]<delta[0]) iclosest=1;
5479 if (delta[iclosest]>cdist1) continue;
5481 helixes[i].ParabolicDCA(helixes[j],phase[iclosest][0],phase[iclosest][1],radius[iclosest],delta[iclosest]);
5482 if (radius[iclosest]<cradius0 || radius[iclosest]>cradius1 || delta[iclosest]>cdist1) continue;
5484 Double_t pointAngle = helixes[i].GetPointAngle(helixes[j],phase[iclosest],fprimvertex);
5485 if (pointAngle<cpointAngle) continue;
5487 Bool_t isGamma = kFALSE;
5488 vertex.SetParamP(*track0); //track0 - plus
5489 vertex.SetParamN(*track1); //track1 - minus
5490 vertex.Update(fprimvertex);
5491 if (track0->TPCrPID(0)>0.3&&track1->TPCrPID(0)>0.3&&vertex.GetAnglep()[2]<0.15) isGamma=kTRUE; // gamma conversion candidate
5492 Double_t pointAngle2 = vertex.GetV0CosineOfPointingAngle();
5494 if (vertex.GetV0CosineOfPointingAngle()<cpointAngle && (!isGamma)) continue;// point angle cut
5495 //Bo: if (vertex.GetDist2()>2&&(!isGamma)) continue; // point angle cut
5496 if (vertex.GetDcaV0Daughters()>2&&(!isGamma)) continue;//Bo: // point angle cut
5497 Float_t sigmae = 0.15*0.15;
5498 if (vertex.GetRr()<80)
5499 sigmae += (sdcar[i]*sdcar[i]+sdcar[j]*sdcar[j])*(1.-vertex.GetRr()/80.)*(1.-vertex.GetRr()/80.);
5500 sigmae+= TMath::Sqrt(sigmae);
5501 //Bo: if (vertex.GetDist2()/sigmae>3.&&(!isGamma)) continue;
5502 if (vertex.GetDcaV0Daughters()/sigmae>3.&&(!isGamma)) continue;
5503 Float_t densb0=0,densb1=0,densa0=0,densa1=0;
5504 Int_t row0 = GetRowNumber(vertex.GetRr());
5506 //Bo: if (vertex.GetDist2()>0.2) continue;
5507 if (vertex.GetDcaV0Daughters()>0.2) continue;
5508 densb0 = track0->Density2(0,row0-5);
5509 densb1 = track1->Density2(0,row0-5);
5510 if (densb0>0.3|| densb1>0.3) continue; //clusters before vertex
5511 densa0 = track0->Density2(row0+5,row0+40);
5512 densa1 = track1->Density2(row0+5,row0+40);
5513 if ((densa0<0.4|| densa1<0.4)&&(!isGamma)) continue; //missing clusters after vertex
5516 densa0 = track0->Density2(0,40); //cluster density
5517 densa1 = track1->Density2(0,40); //cluster density
5518 if ((vertex.GetRr()<80&&densa0+densa1<1.)&&(!isGamma)) continue;
5520 //Bo: vertex.SetLab(0,track0->GetLabel());
5521 //Bo: vertex.SetLab(1,track1->GetLabel());
5522 vertex.SetChi2After((densa0+densa1)*0.5);
5523 vertex.SetChi2Before((densb0+densb1)*0.5);
5524 vertex.SetIndex(0,i);
5525 vertex.SetIndex(1,j);
5526 //Bo: vertex.SetStatus(1); // TPC v0 candidate
5527 vertex.SetOnFlyStatus(2);//Bo: // TPC v0 candidate
5528 //Bo: vertex.SetRp(track0->TPCrPIDs());
5529 //Bo: vertex.SetRm(track1->TPCrPIDs());
5530 tpcv0s->AddLast(new AliESDv0(vertex));
5533 Int_t eventNr = esd->GetEventNumberInFile(); // This is most likely NOT the event number you'd like to use. It has nothing to do with the 'real' event number
5534 Double_t radiusm= (delta[0]<delta[1])? TMath::Sqrt(radius[0]):TMath::Sqrt(radius[1]);
5535 Double_t deltam= (delta[0]<delta[1])? TMath::Sqrt(delta[0]):TMath::Sqrt(delta[1]);
5536 if (AliTPCReconstructor::StreamLevel()>1) {
5537 Int_t lab0=track0->GetLabel();
5538 Int_t lab1=track1->GetLabel();
5539 Char_t c0=track0->GetCircular();
5540 Char_t c1=track1->GetCircular();
5541 TTreeSRedirector &cstream = *fDebugStreamer;
5544 "vertex.="<<&vertex<<
5547 "Helix0.="<<&helixes[i]<<
5550 "Helix1.="<<&helixes[j]<<
5551 "pointAngle="<<pointAngle<<
5552 "pointAngle2="<<pointAngle2<<
5557 "zvertex="<<zvertex<<
5560 "npoints="<<npoints<<
5561 "radius0="<<radius[0]<<
5562 "delta0="<<delta[0]<<
5563 "radius1="<<radius[1]<<
5564 "delta1="<<delta[1]<<
5565 "radiusm="<<radiusm<<
5567 "sdcar0="<<sdcar[i]<<
5568 "sdcar1="<<sdcar[j]<<
5569 "cdcar0="<<cdcar[i]<<
5570 "cdcar1="<<cdcar[j]<<
5571 "pulldcar0="<<pulldcar[i]<<
5572 "pulldcar1="<<pulldcar[j]<<
5573 "pulldcaz0="<<pulldcaz[i]<<
5574 "pulldcaz1="<<pulldcaz[j]<<
5575 "pulldca0="<<pulldca[i]<<
5576 "pulldca1="<<pulldca[j]<<
5586 Float_t *quality = new Float_t[ncandidates];
5587 Int_t *indexes = new Int_t[ncandidates];
5589 for (Int_t i=0;i<ncandidates;i++){
5591 AliESDv0 *v0 = (AliESDv0*)tpcv0s->At(i);
5592 quality[i] = 1./(1.00001-v0->GetV0CosineOfPointingAngle()); //base point angle
5593 // quality[i] /= (0.5+v0->GetDist2());
5594 // quality[i] *= v0->GetChi2After(); //density factor
5596 Int_t index0 = v0->GetIndex(0);
5597 Int_t index1 = v0->GetIndex(1);
5598 //Bo: Double_t minpulldca = TMath::Min(2.+pulldca[v0->GetIndex(0)],(2.+pulldca[v0->GetIndex(1)]) ); //pull
5599 Double_t minpulldca = TMath::Min(2.+pulldca[index0],(2.+pulldca[index1]) );//Bo:
5603 AliTPCseed * track0 = (AliTPCseed*)array->At(index0);
5604 AliTPCseed * track1 = (AliTPCseed*)array->At(index1);
5605 if (track0->TPCrPID(0)>0.3&&track1->TPCrPID(0)>0.3&&v0->GetAnglep()[2]<0.15) quality[i]+=1000000; // gamma conversion candidate
5606 if (track0->TPCrPID(4)>0.9||(track1->TPCrPID(4)>0.9&&minpulldca>4)) quality[i]*=10; // lambda candidate candidate
5609 TMath::Sort(ncandidates,quality,indexes,kTRUE);
5612 for (Int_t i=0;i<ncandidates;i++){
5613 AliESDv0 * v0 = (AliESDv0*)tpcv0s->At(indexes[i]);
5615 Int_t index0 = v0->GetIndex(0);
5616 Int_t index1 = v0->GetIndex(1);
5617 AliTPCseed * track0 = (AliTPCseed*)array->At(index0);
5618 AliTPCseed * track1 = (AliTPCseed*)array->At(index1);
5619 if (!track0||!track1) {
5623 Bool_t accept =kTRUE; //default accept
5624 Int_t *v0indexes0 = track0->GetV0Indexes();
5625 Int_t *v0indexes1 = track1->GetV0Indexes();
5627 Int_t order0 = (v0indexes0[0]!=0) ? 1:0;
5628 Int_t order1 = (v0indexes1[0]!=0) ? 1:0;
5629 if (v0indexes0[1]!=0) order0 =2;
5630 if (v0indexes1[1]!=0) order1 =2;
5632 if (v0indexes0[2]!=0) {order0=3; accept=kFALSE;}
5633 if (v0indexes0[2]!=0) {order1=3; accept=kFALSE;}
5635 AliESDv0 * v02 = v0;
5637 //Bo: v0->SetOrder(0,order0);
5638 //Bo: v0->SetOrder(1,order1);
5639 //Bo: v0->SetOrder(1,order0+order1);
5640 v0->SetOnFlyStatus(kTRUE);
5641 Int_t index = esd->AddV0(v0);
5642 v02 = esd->GetV0(index);
5643 v0indexes0[order0]=index;
5644 v0indexes1[order1]=index;
5648 Int_t eventNr = esd->GetEventNumberInFile(); // This is most likely NOT the event number you'd like to use. It has nothing to do with the 'real' event number
5649 if (AliTPCReconstructor::StreamLevel()>1) {
5650 Int_t lab0=track0->GetLabel();
5651 Int_t lab1=track1->GetLabel();
5652 TTreeSRedirector &cstream = *fDebugStreamer;
5661 "dca0="<<dca[index0]<<
5662 "dca1="<<dca[index1]<<
5666 "quality="<<quality[i]<<
5667 "pulldca0="<<pulldca[index0]<<
5668 "pulldca1="<<pulldca[index1]<<
5692 printf("TPC V0 finder : naccepted\t%d\tncandidates\t%d\tntracks\t%d\tnall\t%d\n",naccepted,ncandidates,ntracks,nall);
5696 Int_t AliTPCtrackerMI::RefitKink(AliTPCseed &mother, AliTPCseed &daughter, const AliESDkink &knk)
5699 // refit kink towards to the vertex
5702 AliKink &kink=(AliKink &)knk;
5704 Int_t row0 = GetRowNumber(kink.GetR());
5705 FollowProlongation(mother,0);
5706 mother.Reset(kFALSE);
5708 FollowProlongation(daughter,row0);
5709 daughter.Reset(kFALSE);
5710 FollowBackProlongation(daughter,158);
5711 daughter.Reset(kFALSE);
5712 Int_t first = TMath::Max(row0-20,30);
5713 Int_t last = TMath::Min(row0+20,140);
5715 const Int_t kNdiv =5;
5716 AliTPCseed param0[kNdiv]; // parameters along the track
5717 AliTPCseed param1[kNdiv]; // parameters along the track
5718 AliKink kinks[kNdiv]; // corresponding kink parameters
5721 for (Int_t irow=0; irow<kNdiv;irow++){
5722 rows[irow] = first +((last-first)*irow)/(kNdiv-1);
5724 // store parameters along the track
5726 for (Int_t irow=0;irow<kNdiv;irow++){
5727 FollowBackProlongation(mother, rows[irow]);
5728 FollowProlongation(daughter,rows[kNdiv-1-irow]);
5729 param0[irow] = mother;
5730 param1[kNdiv-1-irow] = daughter;
5734 for (Int_t irow=0; irow<kNdiv-1;irow++){
5735 if (param0[irow].GetNumberOfClusters()<kNdiv||param1[irow].GetNumberOfClusters()<kNdiv) continue;
5736 kinks[irow].SetMother(param0[irow]);
5737 kinks[irow].SetDaughter(param1[irow]);
5738 kinks[irow].Update();
5741 // choose kink with best "quality"
5743 Double_t mindist = 10000;
5744 for (Int_t irow=0;irow<kNdiv;irow++){
5745 if (param0[irow].GetNumberOfClusters()<20||param1[irow].GetNumberOfClusters()<20) continue;
5746 if (TMath::Abs(kinks[irow].GetR())>240.) continue;
5747 if (TMath::Abs(kinks[irow].GetR())<100.) continue;
5749 Float_t normdist = TMath::Abs(param0[irow].GetX()-kinks[irow].GetR())*(0.1+kink.GetDistance());
5750 normdist/= (param0[irow].GetNumberOfClusters()+param1[irow].GetNumberOfClusters()+40.);
5751 if (normdist < mindist){
5757 if (index==-1) return 0;
5760 param0[index].Reset(kTRUE);
5761 FollowProlongation(param0[index],0);
5763 mother = param0[index];
5764 daughter = param1[index]; // daughter in vertex
5766 kink.SetMother(mother);
5767 kink.SetDaughter(daughter);
5769 kink.SetTPCRow0(GetRowNumber(kink.GetR()));
5770 kink.SetTPCncls(param0[index].GetNumberOfClusters(),0);
5771 kink.SetTPCncls(param1[index].GetNumberOfClusters(),1);
5772 kink.SetLabel(CookLabel(&mother,0.4, 0,kink.GetTPCRow0()),0);
5773 kink.SetLabel(CookLabel(&daughter,0.4, kink.GetTPCRow0(),160),1);
5774 mother.SetLabel(kink.GetLabel(0));
5775 daughter.SetLabel(kink.GetLabel(1));
5781 void AliTPCtrackerMI::UpdateKinkQualityM(AliTPCseed * seed){
5783 // update Kink quality information for mother after back propagation
5785 if (seed->GetKinkIndex(0)>=0) return;
5786 for (Int_t ikink=0;ikink<3;ikink++){
5787 Int_t index = seed->GetKinkIndex(ikink);
5788 if (index>=0) break;
5789 index = TMath::Abs(index)-1;
5790 AliESDkink * kink = fEvent->GetKink(index);
5791 kink->SetTPCDensity(-1,0,0);
5792 kink->SetTPCDensity(1,0,1);
5794 Int_t row0 = kink->GetTPCRow0() - 2 - Int_t( 0.5/ (0.05+kink->GetAngle(2)));
5795 if (row0<15) row0=15;
5797 Int_t row1 = kink->GetTPCRow0() + 2 + Int_t( 0.5/ (0.05+kink->GetAngle(2)));
5798 if (row1>145) row1=145;
5800 Int_t found,foundable,shared;
5801 seed->GetClusterStatistic(0,row0, found, foundable,shared,kFALSE);
5802 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),0,0);
5803 seed->GetClusterStatistic(row1,155, found, foundable,shared,kFALSE);
5804 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),0,1);
5809 void AliTPCtrackerMI::UpdateKinkQualityD(AliTPCseed * seed){
5811 // update Kink quality information for daughter after refit
5813 if (seed->GetKinkIndex(0)<=0) return;
5814 for (Int_t ikink=0;ikink<3;ikink++){
5815 Int_t index = seed->GetKinkIndex(ikink);
5816 if (index<=0) break;
5817 index = TMath::Abs(index)-1;
5818 AliESDkink * kink = fEvent->GetKink(index);
5819 kink->SetTPCDensity(-1,1,0);
5820 kink->SetTPCDensity(-1,1,1);
5822 Int_t row0 = kink->GetTPCRow0() -2 - Int_t( 0.5/ (0.05+kink->GetAngle(2)));
5823 if (row0<15) row0=15;
5825 Int_t row1 = kink->GetTPCRow0() +2 + Int_t( 0.5/ (0.05+kink->GetAngle(2)));
5826 if (row1>145) row1=145;
5828 Int_t found,foundable,shared;
5829 seed->GetClusterStatistic(0,row0, found, foundable,shared,kFALSE);
5830 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),1,0);
5831 seed->GetClusterStatistic(row1,155, found, foundable,shared,kFALSE);
5832 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),1,1);
5838 Int_t AliTPCtrackerMI::CheckKinkPoint(AliTPCseed*seed,AliTPCseed &mother, AliTPCseed &daughter, const AliESDkink &knk)
5841 // check kink point for given track
5842 // if return value=0 kink point not found
5843 // otherwise seed0 correspond to mother particle
5844 // seed1 correspond to daughter particle
5845 // kink parameter of kink point
5846 AliKink &kink=(AliKink &)knk;
5848 Int_t middlerow = (seed->GetFirstPoint()+seed->GetLastPoint())/2;
5849 Int_t first = seed->GetFirstPoint();
5850 Int_t last = seed->GetLastPoint();
5851 if (last-first<20) return 0; // shortest length - 2*30 = 60 pad-rows
5854 AliTPCseed *seed1 = ReSeed(seed,middlerow+20, kTRUE); //middle of chamber
5855 if (!seed1) return 0;
5856 FollowProlongation(*seed1,seed->GetLastPoint()-20);
5857 seed1->Reset(kTRUE);
5858 FollowProlongation(*seed1,158);
5859 seed1->Reset(kTRUE);
5860 last = seed1->GetLastPoint();
5862 AliTPCseed *seed0 = new AliTPCseed(*seed);
5863 seed0->Reset(kFALSE);
5866 AliTPCseed param0[20]; // parameters along the track
5867 AliTPCseed param1[20]; // parameters along the track
5868 AliKink kinks[20]; // corresponding kink parameters
5870 for (Int_t irow=0; irow<20;irow++){
5871 rows[irow] = first +((last-first)*irow)/19;
5873 // store parameters along the track
5875 for (Int_t irow=0;irow<20;irow++){
5876 FollowBackProlongation(*seed0, rows[irow]);
5877 FollowProlongation(*seed1,rows[19-irow]);
5878 param0[irow] = *seed0;
5879 param1[19-irow] = *seed1;
5883 for (Int_t irow=0; irow<19;irow++){
5884 kinks[irow].SetMother(param0[irow]);
5885 kinks[irow].SetDaughter(param1[irow]);
5886 kinks[irow].Update();
5889 // choose kink with biggest change of angle
5891 Double_t maxchange= 0;
5892 for (Int_t irow=1;irow<19;irow++){
5893 if (TMath::Abs(kinks[irow].GetR())>240.) continue;
5894 if (TMath::Abs(kinks[irow].GetR())<110.) continue;
5895 Float_t quality = TMath::Abs(kinks[irow].GetAngle(2))/(3.+TMath::Abs(kinks[irow].GetR()-param0[irow].GetX()));
5896 if ( quality > maxchange){
5897 maxchange = quality;
5904 if (index<0) return 0;
5906 Int_t row0 = GetRowNumber(kinks[index].GetR()); //row 0 estimate
5907 seed0 = new AliTPCseed(param0[index]);
5908 seed1 = new AliTPCseed(param1[index]);
5909 seed0->Reset(kFALSE);
5910 seed1->Reset(kFALSE);
5911 seed0->ResetCovariance(10.);
5912 seed1->ResetCovariance(10.);
5913 FollowProlongation(*seed0,0);
5914 FollowBackProlongation(*seed1,158);
5915 mother = *seed0; // backup mother at position 0
5916 seed0->Reset(kFALSE);
5917 seed1->Reset(kFALSE);
5918 seed0->ResetCovariance(10.);
5919 seed1->ResetCovariance(10.);
5921 first = TMath::Max(row0-20,0);
5922 last = TMath::Min(row0+20,158);
5924 for (Int_t irow=0; irow<20;irow++){
5925 rows[irow] = first +((last-first)*irow)/19;
5927 // store parameters along the track
5929 for (Int_t irow=0;irow<20;irow++){
5930 FollowBackProlongation(*seed0, rows[irow]);
5931 FollowProlongation(*seed1,rows[19-irow]);
5932 param0[irow] = *seed0;
5933 param1[19-irow] = *seed1;
5937 for (Int_t irow=0; irow<19;irow++){
5938 kinks[irow].SetMother(param0[irow]);
5939 kinks[irow].SetDaughter(param1[irow]);
5940 // param0[irow].Dump();
5941 //param1[irow].Dump();
5942 kinks[irow].Update();
5945 // choose kink with biggest change of angle
5948 for (Int_t irow=0;irow<20;irow++){
5949 if (TMath::Abs(kinks[irow].GetR())>250.) continue;
5950 if (TMath::Abs(kinks[irow].GetR())<90.) continue;
5951 Float_t quality = TMath::Abs(kinks[irow].GetAngle(2))/(3.+TMath::Abs(kinks[irow].GetR()-param0[irow].GetX()));
5952 if ( quality > maxchange){
5953 maxchange = quality;
5960 if (index==-1 || param0[index].GetNumberOfClusters()+param1[index].GetNumberOfClusters()<100){
5966 // Float_t anglesigma = TMath::Sqrt(param0[index].fC22+param0[index].fC33+param1[index].fC22+param1[index].fC33);
5968 kink.SetMother(param0[index]);
5969 kink.SetDaughter(param1[index]);
5972 Double_t chi2P2 = param0[index].GetParameter()[2]-param1[index].GetParameter()[2];
5974 chi2P2/=param0[index].GetCovariance()[5]+param1[index].GetCovariance()[5];
5975 Double_t chi2P3 = param0[index].GetParameter()[3]-param1[index].GetParameter()[3];
5977 chi2P3/=param0[index].GetCovariance()[9]+param1[index].GetCovariance()[9];
5979 if (AliTPCReconstructor::StreamLevel()>1) {
5980 (*fDebugStreamer)<<"kinkHpt"<<
5983 "p0.="<<¶m0[index]<<
5984 "p1.="<<¶m1[index]<<
5988 if ( chi2P2+chi2P3<AliTPCReconstructor::GetRecoParam()->GetKinkAngleCutChi2(0)){
5995 row0 = GetRowNumber(kink.GetR());
5996 kink.SetTPCRow0(row0);
5997 kink.SetLabel(CookLabel(seed0,0.5,0,row0),0);
5998 kink.SetLabel(CookLabel(seed1,0.5,row0,158),1);
5999 kink.SetIndex(-10,0);
6000 kink.SetIndex(int(param0[index].GetNumberOfClusters()+param1[index].GetNumberOfClusters()),1);
6001 kink.SetTPCncls(param0[index].GetNumberOfClusters(),0);
6002 kink.SetTPCncls(param1[index].GetNumberOfClusters(),1);
6005 // new (&mother) AliTPCseed(param0[index]);
6006 daughter = param1[index];
6007 daughter.SetLabel(kink.GetLabel(1));
6008 param0[index].Reset(kTRUE);
6009 FollowProlongation(param0[index],0);
6010 mother = param0[index];
6011 mother.SetLabel(kink.GetLabel(0));
6012 if ( chi2P2+chi2P3<AliTPCReconstructor::GetRecoParam()->GetKinkAngleCutChi2(1)){
6024 AliTPCseed* AliTPCtrackerMI::ReSeed(AliTPCseed *t)
6027 // reseed - refit - track
6030 // Int_t last = fSectors->GetNRows()-1;
6032 if (fSectors == fOuterSec){
6033 first = TMath::Max(first, t->GetFirstPoint()-fInnerSec->GetNRows());
6037 first = t->GetFirstPoint();
6039 AliTPCseed * seed = MakeSeed(t,0.1,0.5,0.9);
6040 FollowBackProlongation(*t,fSectors->GetNRows()-1);
6042 FollowProlongation(*t,first);
6052 //_____________________________________________________________________________
6053 Int_t AliTPCtrackerMI::ReadSeeds(const TFile *inp) {
6054 //-----------------------------------------------------------------
6055 // This function reades track seeds.
6056 //-----------------------------------------------------------------
6057 TDirectory *savedir=gDirectory;
6059 TFile *in=(TFile*)inp;
6060 if (!in->IsOpen()) {
6061 cerr<<"AliTPCtrackerMI::ReadSeeds(): input file is not open !\n";
6066 TTree *seedTree=(TTree*)in->Get("Seeds");
6068 cerr<<"AliTPCtrackerMI::ReadSeeds(): ";
6069 cerr<<"can't get a tree with track seeds !\n";
6072 AliTPCtrack *seed=new AliTPCtrack;
6073 seedTree->SetBranchAddress("tracks",&seed);
6075 if (fSeeds==0) fSeeds=new TObjArray(15000);
6077 Int_t n=(Int_t)seedTree->GetEntries();
6078 for (Int_t i=0; i<n; i++) {
6079 seedTree->GetEvent(i);
6080 fSeeds->AddLast(new AliTPCseed(*seed/*,seed->GetAlpha()*/));
6089 Int_t AliTPCtrackerMI::Clusters2Tracks (AliESDEvent *const esd)
6092 if (fSeeds) DeleteSeeds();
6095 if (!fSeeds) return 1;
6102 //_____________________________________________________________________________
6103 Int_t AliTPCtrackerMI::Clusters2Tracks() {
6104 //-----------------------------------------------------------------
6105 // This is a track finder.
6106 //-----------------------------------------------------------------
6107 TDirectory *savedir=gDirectory;
6111 fSeeds = Tracking();
6114 Info("Clusters2Tracks","Time for tracking: \t");timer.Print();timer.Start();
6116 //activate again some tracks
6117 for (Int_t i=0; i<fSeeds->GetEntriesFast(); i++) {
6118 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
6120 Int_t nc=t.GetNumberOfClusters();
6122 delete fSeeds->RemoveAt(i);
6126 if (pt->GetRemoval()==10) {
6127 if (pt->GetDensityFirst(20)>0.8 || pt->GetDensityFirst(30)>0.8 || pt->GetDensityFirst(40)>0.7)
6128 pt->Desactivate(10); // make track again active
6130 pt->Desactivate(20);
6131 delete fSeeds->RemoveAt(i);
6136 RemoveUsed2(fSeeds,0.85,0.85,0);
6137 if (AliTPCReconstructor::GetRecoParam()->GetDoKinks()) FindKinks(fSeeds,fEvent);
6138 //FindCurling(fSeeds, fEvent,0);
6139 if (AliTPCReconstructor::StreamLevel()>5) FindMultiMC(fSeeds, fEvent,-1); // find multi found tracks
6140 RemoveUsed2(fSeeds,0.5,0.4,20);
6141 FindSplitted(fSeeds, fEvent,0); // find multi found tracks
6142 if (AliTPCReconstructor::StreamLevel()>5) FindMultiMC(fSeeds, fEvent,0); // find multi found tracks
6145 // // refit short tracks
6147 Int_t nseed=fSeeds->GetEntriesFast();
6150 for (Int_t i=0; i<nseed; i++) {
6151 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
6153 Int_t nc=t.GetNumberOfClusters();
6155 delete fSeeds->RemoveAt(i);
6158 CookLabel(pt,0.1); //For comparison only
6159 //if ((pt->IsActive() || (pt->fRemoval==10) )&& nc>50 &&pt->GetNumberOfClusters()>0.4*pt->fNFoundable){
6160 if ((pt->IsActive() || (pt->GetRemoval()==10) )){
6162 if (fDebug>0) cerr<<found<<'\r';
6166 delete fSeeds->RemoveAt(i);
6170 //RemoveOverlap(fSeeds,0.99,7,kTRUE);
6172 //RemoveUsed(fSeeds,0.9,0.9,6);
6174 nseed=fSeeds->GetEntriesFast();
6176 for (Int_t i=0; i<nseed; i++) {
6177 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
6179 Int_t nc=t.GetNumberOfClusters();
6181 delete fSeeds->RemoveAt(i);
6185 t.CookdEdx(0.02,0.6);
6186 // CheckKinkPoint(&t,0.05);
6187 //if ((pt->IsActive() || (pt->fRemoval==10) )&& nc>50 &&pt->GetNumberOfClusters()>0.4*pt->fNFoundable){
6188 if ((pt->IsActive() || (pt->GetRemoval()==10) )){
6196 delete fSeeds->RemoveAt(i);
6197 //AliTPCseed * seed1 = ReSeed(pt,0.05,0.5,1);
6199 // FollowProlongation(*seed1,0);
6200 // Int_t n = seed1->GetNumberOfClusters();
6201 // printf("fP4\t%f\t%f\n",seed1->GetC(),pt->GetC());
6202 // printf("fN\t%d\t%d\n", seed1->GetNumberOfClusters(),pt->GetNumberOfClusters());
6205 //AliTPCseed * seed2 = ReSeed(pt,0.95,0.5,0.05);
6209 SortTracks(fSeeds, 1);
6213 PrepareForBackProlongation(fSeeds,5.);
6214 PropagateBack(fSeeds);
6215 printf("Time for back propagation: \t");timer.Print();timer.Start();
6219 PrepareForProlongation(fSeeds,5.);
6220 PropagateForward2(fSeeds);
6222 printf("Time for FORWARD propagation: \t");timer.Print();timer.Start();
6223 // RemoveUsed(fSeeds,0.7,0.7,6);
6224 //RemoveOverlap(fSeeds,0.9,7,kTRUE);
6226 nseed=fSeeds->GetEntriesFast();
6228 for (Int_t i=0; i<nseed; i++) {
6229 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
6231 Int_t nc=t.GetNumberOfClusters();
6233 delete fSeeds->RemoveAt(i);
6236 t.CookdEdx(0.02,0.6);
6237 // CookLabel(pt,0.1); //For comparison only
6238 //if ((pt->IsActive() || (pt->fRemoval==10) )&& nc>50 &&pt->GetNumberOfClusters()>0.4*pt->fNFoundable){
6239 if ((pt->IsActive() || (pt->fRemoval==10) )){
6240 cerr<<found++<<'\r';
6243 delete fSeeds->RemoveAt(i);
6248 // fNTracks = found;
6250 Info("Clusters2Tracks","Time for overlap removal, track writing and dedx cooking: \t"); timer.Print();timer.Start();
6253 // cerr<<"Number of found tracks : "<<"\t"<<found<<endl;
6254 Info("Clusters2Tracks","Number of found tracks %d",found);
6256 // UnloadClusters();
6261 void AliTPCtrackerMI::Tracking(TObjArray * arr)
6264 // tracking of the seeds
6267 fSectors = fOuterSec;
6268 ParallelTracking(arr,150,63);
6269 fSectors = fOuterSec;
6270 ParallelTracking(arr,63,0);
6273 TObjArray * AliTPCtrackerMI::Tracking(Int_t seedtype, Int_t i1, Int_t i2, Float_t cuts[4], Float_t dy, Int_t dsec)
6278 TObjArray * arr = new TObjArray;
6280 fSectors = fOuterSec;
6283 for (Int_t sec=0;sec<fkNOS;sec++){
6284 if (seedtype==3) MakeSeeds3(arr,sec,i1,i2,cuts,dy, dsec);
6285 if (seedtype==4) MakeSeeds5(arr,sec,i1,i2,cuts,dy);
6286 if (seedtype==2) MakeSeeds2(arr,sec,i1,i2,cuts,dy);
6289 Info("Tracking","\nSeeding - %d\t%d\t%d\t%d\n",seedtype,i1,i2,arr->GetEntriesFast());
6301 TObjArray * AliTPCtrackerMI::Tracking()
6305 if (AliTPCReconstructor::GetRecoParam()->GetSpecialSeeding()) return TrackingSpecial();
6308 Int_t nup=fOuterSec->GetNRows()+fInnerSec->GetNRows();
6310 TObjArray * seeds = new TObjArray;
6319 Float_t fnumber = 3.0;
6320 Float_t fdensity = 3.0;
6325 for (Int_t delta = 0; delta<18; delta+=6){
6329 arr = Tracking(3,nup-1-delta,nup-1-delta-gap,cuts,-1,1);
6330 SumTracks(seeds,arr);
6331 SignClusters(seeds,fnumber,fdensity);
6333 for (Int_t i=2;i<6;i+=2){
6334 // seed high pt tracks
6337 arr = Tracking(3,nup-i-delta,nup-i-delta-gap,cuts,-1,0);
6338 SumTracks(seeds,arr);
6339 SignClusters(seeds,fnumber,fdensity);
6344 // RemoveUsed(seeds,0.9,0.9,1);
6345 // UnsignClusters();
6346 // SignClusters(seeds,fnumber,fdensity);
6350 for (Int_t delta = 20; delta<120; delta+=10){
6352 // seed high pt tracks
6356 arr = Tracking(3,nup-delta,nup-delta-gap,cuts,-1);
6357 SumTracks(seeds,arr);
6358 SignClusters(seeds,fnumber,fdensity);
6363 arr = Tracking(3,nup-delta-5,nup-delta-5-gap,cuts,-1);
6364 SumTracks(seeds,arr);
6365 SignClusters(seeds,fnumber,fdensity);
6376 Info("Tracking()","\n\nPrimary seeding\t%d\n\n",seeds->GetEntriesFast());
6380 // RemoveUsed(seeds,0.75,0.75,1);
6382 //SignClusters(seeds,fnumber,fdensity);
6391 arr = Tracking(4,nup-1,nup-1-gap,cuts,-1);
6392 SumTracks(seeds,arr);
6393 SignClusters(seeds,fnumber,fdensity);
6395 arr = Tracking(4,nup-2,nup-2-gap,cuts,-1);
6396 SumTracks(seeds,arr);
6397 SignClusters(seeds,fnumber,fdensity);
6399 arr = Tracking(4,nup-3,nup-3-gap,cuts,-1);
6400 SumTracks(seeds,arr);
6401 SignClusters(seeds,fnumber,fdensity);
6405 for (Int_t delta = 3; delta<30; delta+=5){
6411 arr = Tracking(4,nup-1-delta,nup-1-delta-gap,cuts,-1);
6412 SumTracks(seeds,arr);
6413 SignClusters(seeds,fnumber,fdensity);
6415 arr = Tracking(4,nup-3-delta,nup-5-delta-gap,cuts,4);
6416 SumTracks(seeds,arr);
6417 SignClusters(seeds,fnumber,fdensity);
6428 Int_t fLastSeedRowSec=AliTPCReconstructor::GetRecoParam()->GetLastSeedRowSec();
6431 for (Int_t delta = 30; delta<fLastSeedRowSec; delta+=10){
6437 arr = Tracking(4,nup-1-delta,nup-1-delta-gap,cuts,-1);
6438 SumTracks(seeds,arr);
6439 SignClusters(seeds,fnumber,fdensity);
6441 arr = Tracking(4,nup-5-delta,nup-5-delta-gap,cuts,5 );
6442 SumTracks(seeds,arr);
6443 SignClusters(seeds,fnumber,fdensity);
6447 Info("Tracking()","\n\nSecondary seeding\t%d\n\n",seeds->GetEntriesFast());
6458 TObjArray * AliTPCtrackerMI::TrackingSpecial()
6461 // seeding adjusted for laser and cosmic tests - short tracks with big inclination angle
6462 // no primary vertex seeding tried
6466 Int_t nup=fOuterSec->GetNRows()+fInnerSec->GetNRows();
6468 TObjArray * seeds = new TObjArray;
6473 Float_t fnumber = 3.0;
6474 Float_t fdensity = 3.0;
6477 cuts[0] = AliTPCReconstructor::GetRecoParam()->GetMaxC(); // max curvature
6478 cuts[1] = 3.5; // max tan(phi) angle for seeding
6479 cuts[2] = 3.; // not used (cut on z primary vertex)
6480 cuts[3] = 3.5; // max tan(theta) angle for seeding
6482 for (Int_t delta = 0; nup-delta-gap-1>0; delta+=3){
6484 arr = Tracking(4,nup-1-delta,nup-1-delta-gap,cuts,-1);
6485 SumTracks(seeds,arr);
6486 SignClusters(seeds,fnumber,fdensity);
6490 Info("Tracking()","\n\nSecondary seeding\t%d\n\n",seeds->GetEntriesFast());
6501 void AliTPCtrackerMI::SumTracks(TObjArray *arr1,TObjArray *&arr2) const
6504 //sum tracks to common container
6505 //remove suspicious tracks
6506 Int_t nseed = arr2->GetEntriesFast();
6507 for (Int_t i=0;i<nseed;i++){
6508 AliTPCseed *pt=(AliTPCseed*)arr2->UncheckedAt(i);
6511 // remove tracks with too big curvature
6513 if (TMath::Abs(pt->GetC())>AliTPCReconstructor::GetRecoParam()->GetMaxC()){
6514 delete arr2->RemoveAt(i);
6517 // REMOVE VERY SHORT TRACKS
6518 if (pt->GetNumberOfClusters()<20){
6519 delete arr2->RemoveAt(i);
6522 // NORMAL ACTIVE TRACK
6523 if (pt->IsActive()){
6524 arr1->AddLast(arr2->RemoveAt(i));
6527 //remove not usable tracks
6528 if (pt->GetRemoval()!=10){
6529 delete arr2->RemoveAt(i);
6533 // ENABLE ONLY ENOUGH GOOD STOPPED TRACKS
6534 if (pt->GetDensityFirst(20)>0.8 || pt->GetDensityFirst(30)>0.8 || pt->GetDensityFirst(40)>0.7)
6535 arr1->AddLast(arr2->RemoveAt(i));
6537 delete arr2->RemoveAt(i);
6541 delete arr2; arr2 = 0;
6546 void AliTPCtrackerMI::ParallelTracking(TObjArray *const arr, Int_t rfirst, Int_t rlast)
6549 // try to track in parralel
6551 Int_t nseed=arr->GetEntriesFast();
6552 //prepare seeds for tracking
6553 for (Int_t i=0; i<nseed; i++) {
6554 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i), &t=*pt;
6556 if (!t.IsActive()) continue;
6557 // follow prolongation to the first layer
6558 if ( (fSectors ==fInnerSec) || (t.GetFirstPoint()-fkParam->GetNRowLow()>rfirst+1) )
6559 FollowProlongation(t, rfirst+1);
6564 for (Int_t nr=rfirst; nr>=rlast; nr--){
6565 if (nr<fInnerSec->GetNRows())
6566 fSectors = fInnerSec;
6568 fSectors = fOuterSec;
6569 // make indexes with the cluster tracks for given
6571 // find nearest cluster
6572 for (Int_t i=0; i<nseed; i++) {
6573 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i), &t=*pt;
6575 if (nr==80) pt->UpdateReference();
6576 if (!pt->IsActive()) continue;
6577 // if ( (fSectors ==fOuterSec) && (pt->fFirstPoint-fkParam->GetNRowLow())<nr) continue;
6578 if (pt->GetRelativeSector()>17) {
6581 UpdateClusters(t,nr);
6583 // prolonagate to the nearest cluster - if founded
6584 for (Int_t i=0; i<nseed; i++) {
6585 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
6587 if (!pt->IsActive()) continue;
6588 // if ((fSectors ==fOuterSec) && (pt->fFirstPoint-fkParam->GetNRowLow())<nr) continue;
6589 if (pt->GetRelativeSector()>17) {
6592 FollowToNextCluster(*pt,nr);
6597 void AliTPCtrackerMI::PrepareForBackProlongation(TObjArray *const arr,Float_t fac) const
6601 // if we use TPC track itself we have to "update" covariance
6603 Int_t nseed= arr->GetEntriesFast();
6604 for (Int_t i=0;i<nseed;i++){
6605 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
6609 //rotate to current local system at first accepted point
6610 Int_t index = pt->GetClusterIndex2(pt->GetFirstPoint());
6611 Int_t sec = (index&0xff000000)>>24;
6613 Float_t angle1 = fInnerSec->GetAlpha()*sec+fInnerSec->GetAlphaShift();
6614 if (angle1>TMath::Pi())
6615 angle1-=2.*TMath::Pi();
6616 Float_t angle2 = pt->GetAlpha();
6618 if (TMath::Abs(angle1-angle2)>0.001){
6619 pt->Rotate(angle1-angle2);
6620 //angle2 = pt->GetAlpha();
6621 //pt->fRelativeSector = pt->GetAlpha()/fInnerSec->GetAlpha();
6622 //if (pt->GetAlpha()<0)
6623 // pt->fRelativeSector+=18;
6624 //sec = pt->fRelativeSector;
6633 void AliTPCtrackerMI::PrepareForProlongation(TObjArray *const arr, Float_t fac) const
6637 // if we use TPC track itself we have to "update" covariance
6639 Int_t nseed= arr->GetEntriesFast();
6640 for (Int_t i=0;i<nseed;i++){
6641 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
6644 pt->SetFirstPoint(pt->GetLastPoint());
6652 Int_t AliTPCtrackerMI::PropagateBack(TObjArray *const arr)
6655 // make back propagation
6657 Int_t nseed= arr->GetEntriesFast();
6658 for (Int_t i=0;i<nseed;i++){
6659 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
6660 if (pt&& pt->GetKinkIndex(0)<=0) {
6661 //AliTPCseed *pt2 = new AliTPCseed(*pt);
6662 fSectors = fInnerSec;
6663 //FollowBackProlongation(*pt,fInnerSec->GetNRows()-1);
6664 //fSectors = fOuterSec;
6665 FollowBackProlongation(*pt,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1);
6666 //if (pt->GetNumberOfClusters()<(pt->fEsd->GetTPCclusters(0)) ){
6667 // Error("PropagateBack","Not prolonged track %d",pt->GetLabel());
6668 // FollowBackProlongation(*pt2,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1);
6671 if (pt&& pt->GetKinkIndex(0)>0) {
6672 AliESDkink * kink = fEvent->GetKink(pt->GetKinkIndex(0)-1);
6673 pt->SetFirstPoint(kink->GetTPCRow0());
6674 fSectors = fInnerSec;
6675 FollowBackProlongation(*pt,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1);
6683 Int_t AliTPCtrackerMI::PropagateForward2(TObjArray *const arr)
6686 // make forward propagation
6688 Int_t nseed= arr->GetEntriesFast();
6690 for (Int_t i=0;i<nseed;i++){
6691 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
6693 FollowProlongation(*pt,0);
6702 Int_t AliTPCtrackerMI::PropagateForward()
6705 // propagate track forward
6707 Int_t nseed = fSeeds->GetEntriesFast();
6708 for (Int_t i=0;i<nseed;i++){
6709 AliTPCseed *pt = (AliTPCseed*)fSeeds->UncheckedAt(i);
6711 AliTPCseed &t = *pt;
6712 Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
6713 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
6714 if (alpha < 0. ) alpha += 2.*TMath::Pi();
6715 t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN);
6719 fSectors = fOuterSec;
6720 ParallelTracking(fSeeds,fOuterSec->GetNRows()+fInnerSec->GetNRows()-1,fInnerSec->GetNRows());
6721 fSectors = fInnerSec;
6722 ParallelTracking(fSeeds,fInnerSec->GetNRows()-1,0);
6731 Int_t AliTPCtrackerMI::PropagateBack(AliTPCseed *const pt, Int_t row0, Int_t row1)
6734 // make back propagation, in between row0 and row1
6738 fSectors = fInnerSec;
6741 if (row1<fSectors->GetNRows())
6744 r1 = fSectors->GetNRows()-1;
6746 if (row0<fSectors->GetNRows()&& r1>0 )
6747 FollowBackProlongation(*pt,r1);
6748 if (row1<=fSectors->GetNRows())
6751 r1 = row1 - fSectors->GetNRows();
6752 if (r1<=0) return 0;
6753 if (r1>=fOuterSec->GetNRows()) return 0;
6754 fSectors = fOuterSec;
6755 return FollowBackProlongation(*pt,r1);
6763 void AliTPCtrackerMI::GetShape(AliTPCseed * seed, Int_t row)
6767 AliTPCClusterParam * clparam = AliTPCcalibDB::Instance()->GetClusterParam();
6768 Float_t zdrift = TMath::Abs((fkParam->GetZLength(0)-TMath::Abs(seed->GetZ())));
6769 Int_t type = (seed->GetSector() < fkParam->GetNSector()/2) ? 0: (row>126) ? 1:2;
6770 Double_t angulary = seed->GetSnp();
6771 angulary = angulary*angulary/((1.-angulary)*(1.+angulary));
6772 Double_t angularz = seed->GetTgl()*seed->GetTgl()*(1.+angulary);
6774 Double_t sigmay = clparam->GetRMS0(0,type,zdrift,TMath::Sqrt(TMath::Abs(angulary)));
6775 Double_t sigmaz = clparam->GetRMS0(1,type,zdrift,TMath::Sqrt(TMath::Abs(angularz)));
6776 seed->SetCurrentSigmaY2(sigmay*sigmay);
6777 seed->SetCurrentSigmaZ2(sigmaz*sigmaz);
6778 // Float_t sd2 = TMath::Abs((fkParam->GetZLength(0)-TMath::Abs(seed->GetZ())))*fkParam->GetDiffL()*fkParam->GetDiffL();
6779 // // Float_t padlength = fkParam->GetPadPitchLength(seed->fSector);
6780 // Float_t padlength = GetPadPitchLength(row);
6782 // Float_t sresy = (seed->GetSector() < fkParam->GetNSector()/2) ? 0.2 :0.3;
6783 // seed->SetCurrentSigmaY2(sd2+padlength*padlength*angulary/12.+sresy*sresy);
6785 // Float_t sresz = fkParam->GetZSigma();
6786 // seed->SetCurrentSigmaZ2(sd2+padlength*padlength*angularz*angularz*(1+angulary)/12.+sresz*sresz);
6788 Float_t wy = GetSigmaY(seed);
6789 Float_t wz = GetSigmaZ(seed);
6792 if (TMath::Abs(wy/seed->fCurrentSigmaY2-1)>0.0001 || TMath::Abs(wz/seed->fCurrentSigmaZ2-1)>0.0001 ){
6793 printf("problem\n");
6800 //__________________________________________________________________________
6801 void AliTPCtrackerMI::CookLabel(AliKalmanTrack *tk, Float_t wrong) const {
6802 //--------------------------------------------------------------------
6803 //This function "cooks" a track label. If label<0, this track is fake.
6804 //--------------------------------------------------------------------
6805 AliTPCseed * t = dynamic_cast<AliTPCseed*>(tk);
6807 printf("%s:%d wrong type \n",(char*)__FILE__,__LINE__);
6811 Int_t noc=t->GetNumberOfClusters();
6813 //printf("\nnot founded prolongation\n\n\n");
6819 AliTPCclusterMI *clusters[160];
6821 for (Int_t i=0;i<160;i++) {
6828 for (i=0; i<160 && current<noc; i++) {
6830 Int_t index=t->GetClusterIndex2(i);
6831 if (index<=0) continue;
6832 if (index&0x8000) continue;
6834 //clusters[current]=GetClusterMI(index);
6835 if (t->GetClusterPointer(i)){
6836 clusters[current]=t->GetClusterPointer(i);
6842 Int_t lab=123456789;
6843 for (i=0; i<noc; i++) {
6844 AliTPCclusterMI *c=clusters[i];
6846 lab=TMath::Abs(c->GetLabel(0));
6848 for (j=0; j<noc; j++) if (lb[j]==lab || mx[j]==0) break;
6854 for (i=0; i<noc; i++) if (mx[i]>max) {max=mx[i]; lab=lb[i];}
6856 for (i=0; i<noc; i++) {
6857 AliTPCclusterMI *c=clusters[i];
6859 if (TMath::Abs(c->GetLabel(1)) == lab ||
6860 TMath::Abs(c->GetLabel(2)) == lab ) max++;
6863 if ((1.- Float_t(max)/noc) > wrong) lab=-lab;
6866 Int_t tail=Int_t(0.10*noc);
6869 for (i=1; i<=160&&ind<tail; i++) {
6870 // AliTPCclusterMI *c=clusters[noc-i];
6871 AliTPCclusterMI *c=clusters[i];
6873 if (lab == TMath::Abs(c->GetLabel(0)) ||
6874 lab == TMath::Abs(c->GetLabel(1)) ||
6875 lab == TMath::Abs(c->GetLabel(2))) max++;
6878 if (max < Int_t(0.5*tail)) lab=-lab;
6885 //delete[] clusters;
6889 //__________________________________________________________________________
6890 Int_t AliTPCtrackerMI::CookLabel(AliTPCseed *const t, Float_t wrong,Int_t first, Int_t last) const {
6891 //--------------------------------------------------------------------
6892 //This function "cooks" a track label. If label<0, this track is fake.
6893 //--------------------------------------------------------------------
6894 Int_t noc=t->GetNumberOfClusters();
6896 //printf("\nnot founded prolongation\n\n\n");
6902 AliTPCclusterMI *clusters[160];
6904 for (Int_t i=0;i<160;i++) {
6911 for (i=0; i<160 && current<noc; i++) {
6912 if (i<first) continue;
6913 if (i>last) continue;
6914 Int_t index=t->GetClusterIndex2(i);
6915 if (index<=0) continue;
6916 if (index&0x8000) continue;
6918 //clusters[current]=GetClusterMI(index);
6919 if (t->GetClusterPointer(i)){
6920 clusters[current]=t->GetClusterPointer(i);
6925 if (noc<5) return -1;
6926 Int_t lab=123456789;
6927 for (i=0; i<noc; i++) {
6928 AliTPCclusterMI *c=clusters[i];
6930 lab=TMath::Abs(c->GetLabel(0));
6932 for (j=0; j<noc; j++) if (lb[j]==lab || mx[j]==0) break;
6938 for (i=0; i<noc; i++) if (mx[i]>max) {max=mx[i]; lab=lb[i];}
6940 for (i=0; i<noc; i++) {
6941 AliTPCclusterMI *c=clusters[i];
6943 if (TMath::Abs(c->GetLabel(1)) == lab ||
6944 TMath::Abs(c->GetLabel(2)) == lab ) max++;
6947 if ((1.- Float_t(max)/noc) > wrong) lab=-lab;
6950 Int_t tail=Int_t(0.10*noc);
6953 for (i=1; i<=160&&ind<tail; i++) {
6954 // AliTPCclusterMI *c=clusters[noc-i];
6955 AliTPCclusterMI *c=clusters[i];
6957 if (lab == TMath::Abs(c->GetLabel(0)) ||
6958 lab == TMath::Abs(c->GetLabel(1)) ||
6959 lab == TMath::Abs(c->GetLabel(2))) max++;
6962 if (max < Int_t(0.5*tail)) lab=-lab;
6965 // t->SetLabel(lab);
6969 //delete[] clusters;
6973 Int_t AliTPCtrackerMI::GetRowNumber(Double_t x[3]) const
6975 //return pad row number for given x vector
6976 Float_t phi = TMath::ATan2(x[1],x[0]);
6977 if(phi<0) phi=2.*TMath::Pi()+phi;
6978 // Get the local angle in the sector philoc
6979 const Float_t kRaddeg = 180/3.14159265358979312;
6980 Float_t phiangle = (Int_t (phi*kRaddeg/20.) + 0.5)*20./kRaddeg;
6981 Double_t localx = x[0]*TMath::Cos(phiangle)-x[1]*TMath::Sin(phiangle);
6982 return GetRowNumber(localx);
6987 void AliTPCtrackerMI::MakeBitmaps(AliTPCseed *t)
6989 //-----------------------------------------------------------------------
6990 // Fill the cluster and sharing bitmaps of the track
6991 //-----------------------------------------------------------------------
6993 Int_t firstpoint = 0;
6994 Int_t lastpoint = 159;
6995 AliTPCTrackerPoint *point;
6996 AliTPCclusterMI *cluster;
6998 for (int iter=firstpoint; iter<lastpoint; iter++) {
6999 // Change to cluster pointers to see if we have a cluster at given padrow
7000 cluster = t->GetClusterPointer(iter);
7002 t->SetClusterMapBit(iter, kTRUE);
7003 point = t->GetTrackPoint(iter);
7004 if (point->IsShared())
7005 t->SetSharedMapBit(iter,kTRUE);
7007 t->SetSharedMapBit(iter, kFALSE);
7010 t->SetClusterMapBit(iter, kFALSE);
7011 t->SetSharedMapBit(iter, kFALSE);
7016 Bool_t AliTPCtrackerMI::IsFindable(AliTPCseed & track){
7018 // return flag if there is findable cluster at given position
7021 Float_t z = track.GetZ();
7023 if (TMath::Abs(z)<(AliTPCReconstructor::GetCtgRange()*track.GetX()+kDeltaZ) &&
7024 TMath::Abs(z)<fkParam->GetZLength(0) &&
7025 (TMath::Abs(track.GetSnp())<AliTPCReconstructor::GetMaxSnpTracker()))
7031 void AliTPCtrackerMI::AddCovariance(AliTPCseed * seed){
7033 // Adding systematic error
7034 // !!!! the systematic error for element 4 is in 1/cm not in pt
7036 const Double_t *param = AliTPCReconstructor::GetRecoParam()->GetSystematicError();
7038 for (Int_t i=0;i<15;i++) covar[i]=0;
7044 covar[0] = param[0]*param[0];
7045 covar[2] = param[1]*param[1];
7046 covar[5] = param[2]*param[2];
7047 covar[9] = param[3]*param[3];
7048 Double_t facC = AliTracker::GetBz()*kB2C;
7049 covar[14]= param[4]*param[4]*facC*facC;
7050 seed->AddCovariance(covar);