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); where nis bigger 0
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 (Log level 9) - dump information about
100 // 2.a) cluster error estimate
101 // 3.a) cluster shape estimate
104 //-------------------------------------------------------
109 #include "Riostream.h"
110 #include <TClonesArray.h>
112 #include <TObjArray.h>
115 #include "AliComplexCluster.h"
116 #include "AliESDEvent.h"
117 #include "AliESDtrack.h"
118 #include "AliESDVertex.h"
121 #include "AliHelix.h"
122 #include "AliRunLoader.h"
123 #include "AliTPCClustersRow.h"
124 #include "AliTPCParam.h"
125 #include "AliTPCReconstructor.h"
126 #include "AliTPCpolyTrack.h"
127 #include "AliTPCreco.h"
128 #include "AliTPCseed.h"
130 #include "AliTPCtrackerSector.h"
131 #include "AliTPCtrackerMI.h"
132 #include "TStopwatch.h"
133 #include "AliTPCReconstructor.h"
135 #include "TTreeStream.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>=fParam->GetNInnerSector()) track->SetRow(track->GetRow()+fParam->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
292 Double_t sy2=ErrY2(seed,cluster);
293 Double_t sz2=ErrZ2(seed,cluster);
295 Double_t sdistancey2 = sy2+seed->GetSigmaY2();
296 Double_t sdistancez2 = sz2+seed->GetSigmaZ2();
298 Double_t rdistancey2 = (seed->GetCurrentCluster()->GetY()-seed->GetY())*
299 (seed->GetCurrentCluster()->GetY()-seed->GetY())/sdistancey2;
300 Double_t rdistancez2 = (seed->GetCurrentCluster()->GetZ()-seed->GetZ())*
301 (seed->GetCurrentCluster()->GetZ()-seed->GetZ())/sdistancez2;
303 Double_t rdistance2 = rdistancey2+rdistancez2;
306 if (AliTPCReconstructor::StreamLevel()>5 && seed->GetNumberOfClusters()>20) {
307 Float_t rmsy2 = seed->GetCurrentSigmaY2();
308 Float_t rmsz2 = seed->GetCurrentSigmaZ2();
309 Float_t rmsy2p30 = seed->GetCMeanSigmaY2p30();
310 Float_t rmsz2p30 = seed->GetCMeanSigmaZ2p30();
311 Float_t rmsy2p30R = seed->GetCMeanSigmaY2p30R();
312 Float_t rmsz2p30R = seed->GetCMeanSigmaZ2p30R();
314 AliExternalTrackParam param(*seed);
315 static TVectorD gcl(3),gtr(3);
317 param.GetXYZ(gcl.GetMatrixArray());
318 cluster->GetGlobalXYZ(gclf);
319 gcl[0]=gclf[0]; gcl[1]=gclf[1]; gcl[2]=gclf[2];
320 (*fDebugStreamer)<<"ErrParam"<<
329 "rmsy2p30="<<rmsy2p30<<
330 "rmsz2p30="<<rmsz2p30<<
331 "rmsy2p30R="<<rmsy2p30R<<
332 "rmsz2p30R="<<rmsz2p30R<<
336 if (rdistance2>16) return 3;
339 if ((rdistancey2>9. || rdistancez2>9.) && cluster->GetType()==0)
340 return 2; //suspisiouce - will be changed
342 if ((rdistancey2>6.25 || rdistancez2>6.25) && cluster->GetType()>0)
343 // strict cut on overlaped cluster
344 return 2; //suspisiouce - will be changed
346 if ( (rdistancey2>1. || rdistancez2>6.25 )
347 && cluster->GetType()<0){
348 seed->SetNFoundable(seed->GetNFoundable()-1);
357 //_____________________________________________________________________________
358 AliTPCtrackerMI::AliTPCtrackerMI(const AliTPCParam *par):
360 fkNIS(par->GetNInnerSector()/2),
362 fkNOS(par->GetNOuterSector()/2),
379 //---------------------------------------------------------------------
380 // The main TPC tracker constructor
381 //---------------------------------------------------------------------
382 fInnerSec=new AliTPCtrackerSector[fkNIS];
383 fOuterSec=new AliTPCtrackerSector[fkNOS];
386 for (i=0; i<fkNIS; i++) fInnerSec[i].Setup(par,0);
387 for (i=0; i<fkNOS; i++) fOuterSec[i].Setup(par,1);
390 Int_t nrowlow = par->GetNRowLow();
391 Int_t nrowup = par->GetNRowUp();
394 for (Int_t i=0;i<nrowlow;i++){
395 fXRow[i] = par->GetPadRowRadiiLow(i);
396 fPadLength[i]= par->GetPadPitchLength(0,i);
397 fYMax[i] = fXRow[i]*TMath::Tan(0.5*par->GetInnerAngle());
401 for (Int_t i=0;i<nrowup;i++){
402 fXRow[i+nrowlow] = par->GetPadRowRadiiUp(i);
403 fPadLength[i+nrowlow] = par->GetPadPitchLength(60,i);
404 fYMax[i+nrowlow] = fXRow[i+nrowlow]*TMath::Tan(0.5*par->GetOuterAngle());
407 fDebugStreamer = new TTreeSRedirector("TPCdebug.root");
409 //________________________________________________________________________
410 AliTPCtrackerMI::AliTPCtrackerMI(const AliTPCtrackerMI &t):
431 //------------------------------------
432 // dummy copy constructor
433 //------------------------------------------------------------------
436 AliTPCtrackerMI & AliTPCtrackerMI::operator=(const AliTPCtrackerMI& /*r*/){
437 //------------------------------
439 //--------------------------------------------------------------
442 //_____________________________________________________________________________
443 AliTPCtrackerMI::~AliTPCtrackerMI() {
444 //------------------------------------------------------------------
445 // TPC tracker destructor
446 //------------------------------------------------------------------
453 if (fDebugStreamer) delete fDebugStreamer;
457 void AliTPCtrackerMI::FillESD(TObjArray* arr)
461 //fill esds using updated tracks
463 // write tracks to the event
464 // store index of the track
465 Int_t nseed=arr->GetEntriesFast();
466 //FindKinks(arr,fEvent);
467 for (Int_t i=0; i<nseed; i++) {
468 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
472 if (AliTPCReconstructor::StreamLevel()>1) {
473 (*fDebugStreamer)<<"Track0"<<
477 // pt->PropagateTo(fParam->GetInnerRadiusLow());
478 if (pt->GetKinkIndex(0)<=0){ //don't propagate daughter tracks
479 pt->PropagateTo(fParam->GetInnerRadiusLow());
482 if (( pt->GetPoints()[2]- pt->GetPoints()[0])>5 && pt->GetPoints()[3]>0.8){
484 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
485 iotrack.SetTPCPoints(pt->GetPoints());
486 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
487 iotrack.SetV0Indexes(pt->GetV0Indexes());
488 // iotrack.SetTPCpid(pt->fTPCr);
489 //iotrack.SetTPCindex(i);
490 fEvent->AddTrack(&iotrack);
494 if ( (pt->GetNumberOfClusters()>70)&& (Float_t(pt->GetNumberOfClusters())/Float_t(pt->GetNFoundable()))>0.55) {
496 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
497 iotrack.SetTPCPoints(pt->GetPoints());
498 //iotrack.SetTPCindex(i);
499 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
500 iotrack.SetV0Indexes(pt->GetV0Indexes());
501 // iotrack.SetTPCpid(pt->fTPCr);
502 fEvent->AddTrack(&iotrack);
506 // short tracks - maybe decays
508 if ( (pt->GetNumberOfClusters()>30) && (Float_t(pt->GetNumberOfClusters())/Float_t(pt->GetNFoundable()))>0.70) {
509 Int_t found,foundable,shared;
510 pt->GetClusterStatistic(0,60,found, foundable,shared,kFALSE);
511 if ( (found>20) && (pt->GetNShared()/float(pt->GetNumberOfClusters())<0.2)){
513 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
514 //iotrack.SetTPCindex(i);
515 iotrack.SetTPCPoints(pt->GetPoints());
516 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
517 iotrack.SetV0Indexes(pt->GetV0Indexes());
518 //iotrack.SetTPCpid(pt->fTPCr);
519 fEvent->AddTrack(&iotrack);
524 if ( (pt->GetNumberOfClusters()>20) && (Float_t(pt->GetNumberOfClusters())/Float_t(pt->GetNFoundable()))>0.8) {
525 Int_t found,foundable,shared;
526 pt->GetClusterStatistic(0,60,found, foundable,shared,kFALSE);
527 if (found<20) continue;
528 if (pt->GetNShared()/float(pt->GetNumberOfClusters())>0.2) continue;
531 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
532 iotrack.SetTPCPoints(pt->GetPoints());
533 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
534 iotrack.SetV0Indexes(pt->GetV0Indexes());
535 //iotrack.SetTPCpid(pt->fTPCr);
536 //iotrack.SetTPCindex(i);
537 fEvent->AddTrack(&iotrack);
540 // short tracks - secondaties
542 if ( (pt->GetNumberOfClusters()>30) ) {
543 Int_t found,foundable,shared;
544 pt->GetClusterStatistic(128,158,found, foundable,shared,kFALSE);
545 if ( (found>20) && (pt->GetNShared()/float(pt->GetNumberOfClusters())<0.2) &&float(found)/float(foundable)>0.8){
547 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
548 iotrack.SetTPCPoints(pt->GetPoints());
549 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
550 iotrack.SetV0Indexes(pt->GetV0Indexes());
551 //iotrack.SetTPCpid(pt->fTPCr);
552 //iotrack.SetTPCindex(i);
553 fEvent->AddTrack(&iotrack);
558 if ( (pt->GetNumberOfClusters()>15)) {
559 Int_t found,foundable,shared;
560 pt->GetClusterStatistic(138,158,found, foundable,shared,kFALSE);
561 if (found<15) continue;
562 if (foundable<=0) continue;
563 if (pt->GetNShared()/float(pt->GetNumberOfClusters())>0.2) continue;
564 if (float(found)/float(foundable)<0.8) continue;
567 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
568 iotrack.SetTPCPoints(pt->GetPoints());
569 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
570 iotrack.SetV0Indexes(pt->GetV0Indexes());
571 // iotrack.SetTPCpid(pt->fTPCr);
572 //iotrack.SetTPCindex(i);
573 fEvent->AddTrack(&iotrack);
578 printf("Number of filled ESDs-\t%d\n",fEvent->GetNumberOfTracks());
585 Double_t AliTPCtrackerMI::ErrY2(AliTPCseed* seed, AliTPCclusterMI * cl){
588 // Use calibrated cluster error from OCDB
590 AliTPCClusterParam * clparam = AliTPCcalibDB::Instance()->GetClusterParam();
592 Float_t z = TMath::Abs(fParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
593 Int_t ctype = cl->GetType();
594 Int_t type = (cl->GetRow()<63) ? 0: (cl->GetRow()>126) ? 1:2;
595 Double_t angle = seed->GetSnp()*seed->GetSnp();
596 angle = TMath::Sqrt(TMath::Abs(angle/(1.-angle)));
597 Double_t erry2 = clparam->GetError0Par(0,type, z,angle);
599 erry2+=0.5; // edge cluster
602 seed->SetErrorY2(erry2);
606 //calculate look-up table at the beginning
607 // static Bool_t ginit = kFALSE;
608 // static Float_t gnoise1,gnoise2,gnoise3;
609 // static Float_t ggg1[10000];
610 // static Float_t ggg2[10000];
611 // static Float_t ggg3[10000];
612 // static Float_t glandau1[10000];
613 // static Float_t glandau2[10000];
614 // static Float_t glandau3[10000];
616 // static Float_t gcor01[500];
617 // static Float_t gcor02[500];
618 // static Float_t gcorp[500];
622 // if (ginit==kFALSE){
623 // for (Int_t i=1;i<500;i++){
624 // Float_t rsigma = float(i)/100.;
625 // gcor02[i] = TMath::Max(0.78 +TMath::Exp(7.4*(rsigma-1.2)),0.6);
626 // gcor01[i] = TMath::Max(0.72 +TMath::Exp(3.36*(rsigma-1.2)),0.6);
627 // gcorp[i] = TMath::Max(TMath::Power((rsigma+0.5),1.5),1.2);
631 // for (Int_t i=3;i<10000;i++){
635 // Float_t amp = float(i);
636 // Float_t padlength =0.75;
637 // gnoise1 = 0.0004/padlength;
638 // Float_t nel = 0.268*amp;
639 // Float_t nprim = 0.155*amp;
640 // ggg1[i] = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.001*nel/(padlength*padlength))/nel;
641 // glandau1[i] = (2.+0.12*nprim)*0.5* (2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
642 // if (glandau1[i]>1) glandau1[i]=1;
643 // glandau1[i]*=padlength*padlength/12.;
647 // gnoise2 = 0.0004/padlength;
649 // nprim = 0.133*amp;
650 // ggg2[i] = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
651 // glandau2[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
652 // if (glandau2[i]>1) glandau2[i]=1;
653 // glandau2[i]*=padlength*padlength/12.;
658 // gnoise3 = 0.0004/padlength;
660 // nprim = 0.133*amp;
661 // ggg3[i] = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
662 // glandau3[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
663 // if (glandau3[i]>1) glandau3[i]=1;
664 // glandau3[i]*=padlength*padlength/12.;
672 // Int_t amp = int(TMath::Abs(cl->GetQ()));
674 // seed->SetErrorY2(1.);
678 // Float_t z = TMath::Abs(fParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
679 // Int_t ctype = cl->GetType();
680 // Float_t padlength= GetPadPitchLength(seed->GetRow());
681 // Double_t angle2 = seed->GetSnp()*seed->GetSnp();
682 // angle2 = angle2/(1-angle2);
684 // //cluster "quality"
685 // Int_t rsigmay = int(100.*cl->GetSigmaY2()/(seed->GetCurrentSigmaY2()));
688 // if (fSectors==fInnerSec){
689 // snoise2 = gnoise1;
690 // res = ggg1[amp]*z+glandau1[amp]*angle2;
691 // if (ctype==0) res *= gcor01[rsigmay];
694 // res*= gcorp[rsigmay];
698 // if (padlength<1.1){
699 // snoise2 = gnoise2;
700 // res = ggg2[amp]*z+glandau2[amp]*angle2;
701 // if (ctype==0) res *= gcor02[rsigmay];
704 // res*= gcorp[rsigmay];
708 // snoise2 = gnoise3;
709 // res = ggg3[amp]*z+glandau3[amp]*angle2;
710 // if (ctype==0) res *= gcor02[rsigmay];
713 // res*= gcorp[rsigmay];
720 // res*=2.4; // overestimate error 2 times
724 // if (res<2*snoise2)
727 // seed->SetErrorY2(res);
735 Double_t AliTPCtrackerMI::ErrZ2(AliTPCseed* seed, AliTPCclusterMI * cl){
738 // Use calibrated cluster error from OCDB
740 AliTPCClusterParam * clparam = AliTPCcalibDB::Instance()->GetClusterParam();
742 Float_t z = TMath::Abs(fParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
743 Int_t ctype = cl->GetType();
744 Int_t type = (cl->GetRow()<63) ? 0: (cl->GetRow()>126) ? 1:2;
746 Double_t angle2 = seed->GetSnp()*seed->GetSnp();
747 angle2 = seed->GetTgl()*seed->GetTgl()*(1+angle2/(1-angle2));
748 Double_t angle = TMath::Sqrt(TMath::Abs(angle2));
749 Double_t errz2 = clparam->GetError0Par(1,type, z,angle);
751 errz2+=0.5; // edge cluster
754 seed->SetErrorZ2(errz2);
760 // //seed->SetErrorY2(0.1);
762 // //calculate look-up table at the beginning
763 // static Bool_t ginit = kFALSE;
764 // static Float_t gnoise1,gnoise2,gnoise3;
765 // static Float_t ggg1[10000];
766 // static Float_t ggg2[10000];
767 // static Float_t ggg3[10000];
768 // static Float_t glandau1[10000];
769 // static Float_t glandau2[10000];
770 // static Float_t glandau3[10000];
772 // static Float_t gcor01[1000];
773 // static Float_t gcor02[1000];
774 // static Float_t gcorp[1000];
778 // if (ginit==kFALSE){
779 // for (Int_t i=1;i<1000;i++){
780 // Float_t rsigma = float(i)/100.;
781 // gcor02[i] = TMath::Max(0.81 +TMath::Exp(6.8*(rsigma-1.2)),0.6);
782 // gcor01[i] = TMath::Max(0.72 +TMath::Exp(2.04*(rsigma-1.2)),0.6);
783 // gcorp[i] = TMath::Max(TMath::Power((rsigma+0.5),1.5),1.2);
787 // for (Int_t i=3;i<10000;i++){
791 // Float_t amp = float(i);
792 // Float_t padlength =0.75;
793 // gnoise1 = 0.0004/padlength;
794 // Float_t nel = 0.268*amp;
795 // Float_t nprim = 0.155*amp;
796 // ggg1[i] = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.001*nel/(padlength*padlength))/nel;
797 // glandau1[i] = (2.+0.12*nprim)*0.5* (2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
798 // if (glandau1[i]>1) glandau1[i]=1;
799 // glandau1[i]*=padlength*padlength/12.;
803 // gnoise2 = 0.0004/padlength;
805 // nprim = 0.133*amp;
806 // ggg2[i] = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
807 // glandau2[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
808 // if (glandau2[i]>1) glandau2[i]=1;
809 // glandau2[i]*=padlength*padlength/12.;
814 // gnoise3 = 0.0004/padlength;
816 // nprim = 0.133*amp;
817 // ggg3[i] = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
818 // glandau3[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
819 // if (glandau3[i]>1) glandau3[i]=1;
820 // glandau3[i]*=padlength*padlength/12.;
828 // Int_t amp = int(TMath::Abs(cl->GetQ()));
830 // seed->SetErrorY2(1.);
834 // Float_t z = TMath::Abs(fParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
835 // Int_t ctype = cl->GetType();
836 // Float_t padlength= GetPadPitchLength(seed->GetRow());
838 // Double_t angle2 = seed->GetSnp()*seed->GetSnp();
839 // // if (angle2<0.6) angle2 = 0.6;
840 // angle2 = seed->GetTgl()*seed->GetTgl()*(1+angle2/(1-angle2));
842 // //cluster "quality"
843 // Int_t rsigmaz = int(100.*cl->GetSigmaZ2()/(seed->GetCurrentSigmaZ2()));
846 // if (fSectors==fInnerSec){
847 // snoise2 = gnoise1;
848 // res = ggg1[amp]*z+glandau1[amp]*angle2;
849 // if (ctype==0) res *= gcor01[rsigmaz];
852 // res*= gcorp[rsigmaz];
856 // if (padlength<1.1){
857 // snoise2 = gnoise2;
858 // res = ggg2[amp]*z+glandau2[amp]*angle2;
859 // if (ctype==0) res *= gcor02[rsigmaz];
862 // res*= gcorp[rsigmaz];
866 // snoise2 = gnoise3;
867 // res = ggg3[amp]*z+glandau3[amp]*angle2;
868 // if (ctype==0) res *= gcor02[rsigmaz];
871 // res*= gcorp[rsigmaz];
880 // if ((ctype<0) &&<70){
885 // if (res<2*snoise2)
887 // if (res>3) res =3;
888 // seed->SetErrorZ2(res);
896 void AliTPCtrackerMI::RotateToLocal(AliTPCseed *seed)
898 //rotate to track "local coordinata
899 Float_t x = seed->GetX();
900 Float_t y = seed->GetY();
901 Float_t ymax = x*TMath::Tan(0.5*fSectors->GetAlpha());
904 seed->SetRelativeSector((seed->GetRelativeSector()+1) % fN);
905 if (!seed->Rotate(fSectors->GetAlpha()))
907 } else if (y <-ymax) {
908 seed->SetRelativeSector((seed->GetRelativeSector()-1+fN) % fN);
909 if (!seed->Rotate(-fSectors->GetAlpha()))
917 //_____________________________________________________________________________
918 Double_t AliTPCtrackerMI::F1old(Double_t x1,Double_t y1,
919 Double_t x2,Double_t y2,
920 Double_t x3,Double_t y3)
922 //-----------------------------------------------------------------
923 // Initial approximation of the track curvature
924 //-----------------------------------------------------------------
925 Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
926 Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
927 (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
928 Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
929 (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
931 Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
932 if ( xr*xr+yr*yr<=0.00000000000001) return 100;
933 return -xr*yr/sqrt(xr*xr+yr*yr);
938 //_____________________________________________________________________________
939 Double_t AliTPCtrackerMI::F1(Double_t x1,Double_t y1,
940 Double_t x2,Double_t y2,
941 Double_t x3,Double_t y3)
943 //-----------------------------------------------------------------
944 // Initial approximation of the track curvature
945 //-----------------------------------------------------------------
951 Double_t det = x3*y2-x2*y3;
956 Double_t u = 0.5* (x2*(x2-x3)+y2*(y2-y3))/det;
957 Double_t x0 = x3*0.5-y3*u;
958 Double_t y0 = y3*0.5+x3*u;
959 Double_t c2 = 1/TMath::Sqrt(x0*x0+y0*y0);
965 Double_t AliTPCtrackerMI::F2(Double_t x1,Double_t y1,
966 Double_t x2,Double_t y2,
967 Double_t x3,Double_t y3)
969 //-----------------------------------------------------------------
970 // Initial approximation of the track curvature
971 //-----------------------------------------------------------------
977 Double_t det = x3*y2-x2*y3;
982 Double_t u = 0.5* (x2*(x2-x3)+y2*(y2-y3))/det;
983 Double_t x0 = x3*0.5-y3*u;
984 Double_t y0 = y3*0.5+x3*u;
985 Double_t c2 = 1/TMath::Sqrt(x0*x0+y0*y0);
994 //_____________________________________________________________________________
995 Double_t AliTPCtrackerMI::F2old(Double_t x1,Double_t y1,
996 Double_t x2,Double_t y2,
997 Double_t x3,Double_t y3)
999 //-----------------------------------------------------------------
1000 // Initial approximation of the track curvature times center of curvature
1001 //-----------------------------------------------------------------
1002 Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
1003 Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
1004 (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
1005 Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
1006 (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
1008 Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
1010 return -a/(d*y1-b)*xr/sqrt(xr*xr+yr*yr);
1013 //_____________________________________________________________________________
1014 Double_t AliTPCtrackerMI::F3(Double_t x1,Double_t y1,
1015 Double_t x2,Double_t y2,
1016 Double_t z1,Double_t z2)
1018 //-----------------------------------------------------------------
1019 // Initial approximation of the tangent of the track dip angle
1020 //-----------------------------------------------------------------
1021 return (z1 - z2)/sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
1025 Double_t AliTPCtrackerMI::F3n(Double_t x1,Double_t y1,
1026 Double_t x2,Double_t y2,
1027 Double_t z1,Double_t z2, Double_t c)
1029 //-----------------------------------------------------------------
1030 // Initial approximation of the tangent of the track dip angle
1031 //-----------------------------------------------------------------
1035 //angle1 = (z1-z2)*c/(TMath::ASin(c*x1-ni)-TMath::ASin(c*x2-ni));
1037 Double_t d = TMath::Sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
1038 if (TMath::Abs(d*c*0.5)>1) return 0;
1039 // Double_t angle2 = TMath::ASin(d*c*0.5);
1040 // Double_t angle2 = AliTPCFastMath::FastAsin(d*c*0.5);
1041 Double_t angle2 = (d*c*0.5>0.1)? TMath::ASin(d*c*0.5): AliTPCFastMath::FastAsin(d*c*0.5);
1043 angle2 = (z1-z2)*c/(angle2*2.);
1047 Bool_t AliTPCtrackerMI::GetProlongation(Double_t x1, Double_t x2, Double_t x[5], Double_t &y, Double_t &z)
1048 {//-----------------------------------------------------------------
1049 // This function find proloncation of a track to a reference plane x=x2.
1050 //-----------------------------------------------------------------
1054 if (TMath::Abs(x[4]*x1 - x[2]) >= 0.999) {
1058 Double_t c1=x[4]*x1 - x[2], r1=sqrt(1.- c1*c1);
1059 Double_t c2=x[4]*x2 - x[2], r2=sqrt(1.- c2*c2);
1063 Double_t dy = dx*(c1+c2)/(r1+r2);
1066 Double_t delta = x[4]*dx*(c1+c2)/(c1*r2 + c2*r1);
1068 if (TMath::Abs(delta)>0.01){
1069 dz = x[3]*TMath::ASin(delta)/x[4];
1071 dz = x[3]*AliTPCFastMath::FastAsin(delta)/x[4];
1074 //dz = x[3]*AliTPCFastMath::FastAsin(delta)/x[4];
1082 Int_t AliTPCtrackerMI::LoadClusters (TTree *tree)
1087 return LoadClusters();
1091 Int_t AliTPCtrackerMI::LoadClusters(TObjArray *arr)
1094 // load clusters to the memory
1095 AliTPCClustersRow *clrow = 0x0;
1096 Int_t lower = arr->LowerBound();
1097 Int_t entries = arr->GetEntriesFast();
1099 for (Int_t i=lower; i<entries; i++) {
1100 clrow = (AliTPCClustersRow*) arr->At(i);
1101 if(!clrow) continue;
1102 if(!clrow->GetArray()) continue;
1106 fParam->AdjustSectorRow(clrow->GetID(),sec,row);
1108 for (Int_t icl=0; icl<clrow->GetArray()->GetEntriesFast(); icl++){
1109 Transform((AliTPCclusterMI*)(clrow->GetArray()->At(icl)));
1112 if (clrow->GetArray()->GetEntriesFast()<=0) continue;
1113 AliTPCtrackerRow * tpcrow=0;
1116 tpcrow = &(fInnerSec[sec%fkNIS][row]);
1120 tpcrow = &(fOuterSec[(sec-fkNIS*2)%fkNOS][row]);
1121 left = (sec-fkNIS*2)/fkNOS;
1124 tpcrow->SetN1(clrow->GetArray()->GetEntriesFast());
1125 tpcrow->SetClusters1(new AliTPCclusterMI[tpcrow->GetN1()]);
1126 for (Int_t i=0;i<tpcrow->GetN1();i++)
1127 tpcrow->SetCluster1(i, *(AliTPCclusterMI*)(clrow->GetArray()->At(i)));
1130 tpcrow->SetN2(clrow->GetArray()->GetEntriesFast());
1131 tpcrow->SetClusters2(new AliTPCclusterMI[tpcrow->GetN2()]);
1132 for (Int_t i=0;i<tpcrow->GetN2();i++)
1133 tpcrow->SetCluster2(i,*(AliTPCclusterMI*)(clrow->GetArray()->At(i)));
1143 Int_t AliTPCtrackerMI::LoadClusters(TClonesArray *arr)
1146 // load clusters to the memory from one
1149 AliTPCclusterMI *clust=0;
1150 Int_t count[72][96] = { {0} , {0} };
1152 // loop over clusters
1153 for (Int_t icl=0; icl<arr->GetEntriesFast(); icl++) {
1154 clust = (AliTPCclusterMI*)arr->At(icl);
1155 if(!clust) continue;
1156 //printf("cluster: det %d, row %d \n", clust->GetDetector(),clust->GetRow());
1158 // transform clusters
1161 // count clusters per pad row
1162 count[clust->GetDetector()][clust->GetRow()]++;
1165 // insert clusters to sectors
1166 for (Int_t icl=0; icl<arr->GetEntriesFast(); icl++) {
1167 clust = (AliTPCclusterMI*)arr->At(icl);
1168 if(!clust) continue;
1170 Int_t sec = clust->GetDetector();
1171 Int_t row = clust->GetRow();
1173 // filter overlapping pad rows needed by HLT
1174 if(sec<fkNIS*2) { //IROCs
1175 if(row == 30) continue;
1178 if(row == 27 || row == 76) continue;
1184 fInnerSec[sec%fkNIS].InsertCluster(clust, count[sec][row], fParam);
1187 left = (sec-fkNIS*2)/fkNOS;
1188 fOuterSec[(sec-fkNIS*2)%fkNOS].InsertCluster(clust, count[sec][row], fParam);
1192 // Load functions must be called behind LoadCluster(TClonesArray*)
1194 //LoadOuterSectors();
1195 //LoadInnerSectors();
1201 Int_t AliTPCtrackerMI::LoadClusters()
1204 // load clusters to the memory
1205 AliTPCClustersRow *clrow= new AliTPCClustersRow;
1206 clrow->SetClass("AliTPCclusterMI");
1208 clrow->GetArray()->ExpandCreateFast(10000);
1210 // TTree * tree = fClustersArray.GetTree();
1212 TTree * tree = fInput;
1213 TBranch * br = tree->GetBranch("Segment");
1214 br->SetAddress(&clrow);
1216 Int_t j=Int_t(tree->GetEntries());
1217 for (Int_t i=0; i<j; i++) {
1221 fParam->AdjustSectorRow(clrow->GetID(),sec,row);
1222 for (Int_t icl=0; icl<clrow->GetArray()->GetEntriesFast(); icl++){
1223 Transform((AliTPCclusterMI*)(clrow->GetArray()->At(icl)));
1226 AliTPCtrackerRow * tpcrow=0;
1229 tpcrow = &(fInnerSec[sec%fkNIS][row]);
1233 tpcrow = &(fOuterSec[(sec-fkNIS*2)%fkNOS][row]);
1234 left = (sec-fkNIS*2)/fkNOS;
1237 tpcrow->SetN1(clrow->GetArray()->GetEntriesFast());
1238 tpcrow->SetClusters1(new AliTPCclusterMI[tpcrow->GetN1()]);
1239 for (Int_t i=0;i<tpcrow->GetN1();i++)
1240 tpcrow->SetCluster1(i, *(AliTPCclusterMI*)(clrow->GetArray()->At(i)));
1243 tpcrow->SetN2(clrow->GetArray()->GetEntriesFast());
1244 tpcrow->SetClusters2(new AliTPCclusterMI[tpcrow->GetN2()]);
1245 for (Int_t i=0;i<tpcrow->GetN2();i++)
1246 tpcrow->SetCluster2(i,*(AliTPCclusterMI*)(clrow->GetArray()->At(i)));
1257 void AliTPCtrackerMI::UnloadClusters()
1260 // unload clusters from the memory
1262 Int_t nrows = fOuterSec->GetNRows();
1263 for (Int_t sec = 0;sec<fkNOS;sec++)
1264 for (Int_t row = 0;row<nrows;row++){
1265 AliTPCtrackerRow* tpcrow = &(fOuterSec[sec%fkNOS][row]);
1267 // if (tpcrow->fClusters1) delete []tpcrow->fClusters1;
1268 // if (tpcrow->fClusters2) delete []tpcrow->fClusters2;
1270 tpcrow->ResetClusters();
1273 nrows = fInnerSec->GetNRows();
1274 for (Int_t sec = 0;sec<fkNIS;sec++)
1275 for (Int_t row = 0;row<nrows;row++){
1276 AliTPCtrackerRow* tpcrow = &(fInnerSec[sec%fkNIS][row]);
1278 // if (tpcrow->fClusters1) delete []tpcrow->fClusters1;
1279 //if (tpcrow->fClusters2) delete []tpcrow->fClusters2;
1281 tpcrow->ResetClusters();
1287 void AliTPCtrackerMI::FillClusterArray(TObjArray* array) const{
1289 // Filling cluster to the array - For visualization purposes
1292 nrows = fOuterSec->GetNRows();
1293 for (Int_t sec = 0;sec<fkNOS;sec++)
1294 for (Int_t row = 0;row<nrows;row++){
1295 AliTPCtrackerRow* tpcrow = &(fOuterSec[sec%fkNOS][row]);
1296 if (!tpcrow) continue;
1297 for (Int_t icl = 0;icl<tpcrow->GetN();icl++){
1298 array->AddLast((TObject*)((*tpcrow)[icl]));
1301 nrows = fInnerSec->GetNRows();
1302 for (Int_t sec = 0;sec<fkNIS;sec++)
1303 for (Int_t row = 0;row<nrows;row++){
1304 AliTPCtrackerRow* tpcrow = &(fInnerSec[sec%fkNIS][row]);
1305 if (!tpcrow) continue;
1306 for (Int_t icl = 0;icl<tpcrow->GetN();icl++){
1307 array->AddLast((TObject*)(*tpcrow)[icl]);
1313 void AliTPCtrackerMI::Transform(AliTPCclusterMI * cluster){
1317 AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ;
1319 AliFatal("Tranformations not in calibDB");
1321 Double_t x[3]={cluster->GetRow(),cluster->GetPad(),cluster->GetTimeBin()};
1322 Int_t i[1]={cluster->GetDetector()};
1323 transform->Transform(x,i,0,1);
1324 if (cluster->GetDetector()%36>17){
1329 // in debug mode check the transformation
1331 if (AliTPCReconstructor::StreamLevel()>1) {
1333 cluster->GetGlobalXYZ(gx);
1334 Int_t event = (fEvent==NULL)? 0: fEvent->GetEventNumberInFile();
1335 TTreeSRedirector &cstream = *fDebugStreamer;
1336 cstream<<"Transform"<<
1347 cluster->SetX(x[0]);
1348 cluster->SetY(x[1]);
1349 cluster->SetZ(x[2]);
1355 //if (!fParam->IsGeoRead()) fParam->ReadGeoMatrices();
1356 TGeoHMatrix *mat = fParam->GetClusterMatrix(cluster->GetDetector());
1358 Double_t pos[3]= {cluster->GetX(),cluster->GetY(),cluster->GetZ()};
1359 Double_t posC[3]={cluster->GetX(),cluster->GetY(),cluster->GetZ()};
1360 if (mat) mat->LocalToMaster(pos,posC);
1362 // chack Loading of Geo matrices from GeoManager - TEMPORARY FIX
1364 cluster->SetX(posC[0]);
1365 cluster->SetY(posC[1]);
1366 cluster->SetZ(posC[2]);
1369 //_____________________________________________________________________________
1370 Int_t AliTPCtrackerMI::LoadOuterSectors() {
1371 //-----------------------------------------------------------------
1372 // This function fills outer TPC sectors with clusters.
1373 //-----------------------------------------------------------------
1374 Int_t nrows = fOuterSec->GetNRows();
1376 for (Int_t sec = 0;sec<fkNOS;sec++)
1377 for (Int_t row = 0;row<nrows;row++){
1378 AliTPCtrackerRow* tpcrow = &(fOuterSec[sec%fkNOS][row]);
1379 Int_t sec2 = sec+2*fkNIS;
1381 Int_t ncl = tpcrow->GetN1();
1383 AliTPCclusterMI *c= (tpcrow->GetCluster1(ncl));
1384 index=(((sec2<<8)+row)<<16)+ncl;
1385 tpcrow->InsertCluster(c,index);
1388 ncl = tpcrow->GetN2();
1390 AliTPCclusterMI *c= (tpcrow->GetCluster2(ncl));
1391 index=((((sec2+fkNOS)<<8)+row)<<16)+ncl;
1392 tpcrow->InsertCluster(c,index);
1395 // write indexes for fast acces
1397 for (Int_t i=0;i<510;i++)
1398 tpcrow->SetFastCluster(i,-1);
1399 for (Int_t i=0;i<tpcrow->GetN();i++){
1400 Int_t zi = Int_t((*tpcrow)[i]->GetZ()+255.);
1401 tpcrow->SetFastCluster(zi,i); // write index
1404 for (Int_t i=0;i<510;i++){
1405 if (tpcrow->GetFastCluster(i)<0)
1406 tpcrow->SetFastCluster(i,last);
1408 last = tpcrow->GetFastCluster(i);
1417 //_____________________________________________________________________________
1418 Int_t AliTPCtrackerMI::LoadInnerSectors() {
1419 //-----------------------------------------------------------------
1420 // This function fills inner TPC sectors with clusters.
1421 //-----------------------------------------------------------------
1422 Int_t nrows = fInnerSec->GetNRows();
1424 for (Int_t sec = 0;sec<fkNIS;sec++)
1425 for (Int_t row = 0;row<nrows;row++){
1426 AliTPCtrackerRow* tpcrow = &(fInnerSec[sec%fkNIS][row]);
1429 Int_t ncl = tpcrow->GetN1();
1431 AliTPCclusterMI *c= (tpcrow->GetCluster1(ncl));
1432 index=(((sec<<8)+row)<<16)+ncl;
1433 tpcrow->InsertCluster(c,index);
1436 ncl = tpcrow->GetN2();
1438 AliTPCclusterMI *c= (tpcrow->GetCluster2(ncl));
1439 index=((((sec+fkNIS)<<8)+row)<<16)+ncl;
1440 tpcrow->InsertCluster(c,index);
1443 // write indexes for fast acces
1445 for (Int_t i=0;i<510;i++)
1446 tpcrow->SetFastCluster(i,-1);
1447 for (Int_t i=0;i<tpcrow->GetN();i++){
1448 Int_t zi = Int_t((*tpcrow)[i]->GetZ()+255.);
1449 tpcrow->SetFastCluster(zi,i); // write index
1452 for (Int_t i=0;i<510;i++){
1453 if (tpcrow->GetFastCluster(i)<0)
1454 tpcrow->SetFastCluster(i,last);
1456 last = tpcrow->GetFastCluster(i);
1468 //_________________________________________________________________________
1469 AliTPCclusterMI *AliTPCtrackerMI::GetClusterMI(Int_t index) const {
1470 //--------------------------------------------------------------------
1471 // Return pointer to a given cluster
1472 //--------------------------------------------------------------------
1473 if (index<0) return 0; // no cluster
1474 Int_t sec=(index&0xff000000)>>24;
1475 Int_t row=(index&0x00ff0000)>>16;
1476 Int_t ncl=(index&0x00007fff)>>00;
1478 const AliTPCtrackerRow * tpcrow=0;
1479 AliTPCclusterMI * clrow =0;
1481 if (sec<0 || sec>=fkNIS*4) {
1482 AliWarning(Form("Wrong sector %d",sec));
1487 tpcrow = &(fInnerSec[sec%fkNIS][row]);
1488 if (tpcrow==0) return 0;
1491 if (tpcrow->GetN1()<=ncl) return 0;
1492 clrow = tpcrow->GetClusters1();
1495 if (tpcrow->GetN2()<=ncl) return 0;
1496 clrow = tpcrow->GetClusters2();
1500 tpcrow = &(fOuterSec[(sec-fkNIS*2)%fkNOS][row]);
1501 if (tpcrow==0) return 0;
1503 if (sec-2*fkNIS<fkNOS) {
1504 if (tpcrow->GetN1()<=ncl) return 0;
1505 clrow = tpcrow->GetClusters1();
1508 if (tpcrow->GetN2()<=ncl) return 0;
1509 clrow = tpcrow->GetClusters2();
1513 return &(clrow[ncl]);
1519 Int_t AliTPCtrackerMI::FollowToNext(AliTPCseed& t, Int_t nr) {
1520 //-----------------------------------------------------------------
1521 // This function tries to find a track prolongation to next pad row
1522 //-----------------------------------------------------------------
1524 Double_t x= GetXrow(nr), ymax=GetMaxY(nr);
1525 AliTPCclusterMI *cl=0;
1526 Int_t tpcindex= t.GetClusterIndex2(nr);
1528 // update current shape info every 5 pad-row
1529 // if ( (nr%5==0) || t.GetNumberOfClusters()<2 || (t.fCurrentSigmaY2<0.0001) ){
1533 if (fIteration>0 && tpcindex>=-1){ //if we have already clusters
1535 if (tpcindex==-1) return 0; //track in dead zone
1537 cl = t.GetClusterPointer(nr);
1538 if ( (cl==0) ) cl = GetClusterMI(tpcindex);
1539 t.SetCurrentClusterIndex1(tpcindex);
1542 Int_t relativesector = ((tpcindex&0xff000000)>>24)%18; // if previously accepted cluster in different sector
1543 Float_t angle = relativesector*fSectors->GetAlpha()+fSectors->GetAlphaShift();
1545 if (angle<-TMath::Pi()) angle += 2*TMath::Pi();
1546 if (angle>=TMath::Pi()) angle -= 2*TMath::Pi();
1548 if (TMath::Abs(angle-t.GetAlpha())>0.001){
1549 Double_t rotation = angle-t.GetAlpha();
1550 t.SetRelativeSector(relativesector);
1551 if (!t.Rotate(rotation)) return 0;
1553 if (!t.PropagateTo(x)) return 0;
1555 t.SetCurrentCluster(cl);
1557 Int_t accept = AcceptCluster(&t,t.GetCurrentCluster());
1558 if ((tpcindex&0x8000)==0) accept =0;
1560 //if founded cluster is acceptible
1561 if (cl->IsUsed(11)) { // id cluster is shared inrease uncertainty
1562 t.SetErrorY2(t.GetErrorY2()+0.03);
1563 t.SetErrorZ2(t.GetErrorZ2()+0.03);
1564 t.SetErrorY2(t.GetErrorY2()*3);
1565 t.SetErrorZ2(t.GetErrorZ2()*3);
1567 t.SetNFoundable(t.GetNFoundable()+1);
1568 UpdateTrack(&t,accept);
1573 if (TMath::Abs(t.GetSnp())>AliTPCReconstructor::GetMaxSnpTracker()) return 0; // cut on angle
1575 // not look for new cluster during refitting
1576 t.SetNFoundable(t.GetNFoundable()+1);
1581 // if (TMath::Abs(t.GetSnp())>0.95 || TMath::Abs(x*t.GetC()-t.GetEta())>0.95) return 0;// patch 28 fev 06
1582 Double_t y=t.GetYat(x);
1583 if (TMath::Abs(y)>ymax){
1585 t.SetRelativeSector((t.GetRelativeSector()+1) % fN);
1586 if (!t.Rotate(fSectors->GetAlpha()))
1588 } else if (y <-ymax) {
1589 t.SetRelativeSector((t.GetRelativeSector()-1+fN) % fN);
1590 if (!t.Rotate(-fSectors->GetAlpha()))
1596 if (!t.PropagateTo(x)) {
1597 if (fIteration==0) t.SetRemoval(10);
1601 Double_t z=t.GetZ();
1604 if (!IsActive(t.GetRelativeSector(),nr)) {
1606 t.SetClusterIndex2(nr,-1);
1609 //AliInfo(Form("A - Sector%d phi %f - alpha %f", t.fRelativeSector,y/x, t.GetAlpha()));
1610 Bool_t isActive = IsActive(t.GetRelativeSector(),nr);
1611 Bool_t isActive2 = (nr>=fInnerSec->GetNRows()) ? fOuterSec[t.GetRelativeSector()][nr-fInnerSec->GetNRows()].GetN()>0:fInnerSec[t.GetRelativeSector()][nr].GetN()>0;
1613 if (!isActive || !isActive2) return 0;
1615 const AliTPCtrackerRow &krow=GetRow(t.GetRelativeSector(),nr);
1616 if ( (t.GetSigmaY2()<0) || t.GetSigmaZ2()<0) return 0;
1618 Double_t roadz = 1.;
1620 if (TMath::Abs(TMath::Abs(y)-ymax)<krow.GetDeadZone()){
1622 t.SetClusterIndex2(nr,-1);
1627 if (TMath::Abs(z)<(AliTPCReconstructor::GetCtgRange()*x+10) && TMath::Abs(z)<fParam->GetZLength(0) && (TMath::Abs(t.GetSnp())<AliTPCReconstructor::GetMaxSnpTracker()))
1628 t.SetNFoundable(t.GetNFoundable()+1);
1634 // cl = krow.FindNearest2(y+10.,z,roady,roadz,index);
1635 cl = krow.FindNearest2(y,z,roady,roadz,index);
1636 if (cl) t.SetCurrentClusterIndex1(krow.GetIndex(index));
1639 t.SetCurrentCluster(cl);
1641 if (fIteration==2&&cl->IsUsed(10)) return 0;
1642 Int_t accept = AcceptCluster(&t,t.GetCurrentCluster());
1643 if (fIteration==2&&cl->IsUsed(11)) {
1644 t.SetErrorY2(t.GetErrorY2()+0.03);
1645 t.SetErrorZ2(t.GetErrorZ2()+0.03);
1646 t.SetErrorY2(t.GetErrorY2()*3);
1647 t.SetErrorZ2(t.GetErrorZ2()*3);
1650 if (t.fCurrentCluster->IsUsed(10)){
1655 if (t.fNShared>0.7*t.GetNumberOfClusters()) {
1661 if (accept<3) UpdateTrack(&t,accept);
1664 if ( fIteration==0 && t.GetNFoundable()*0.5 > t.GetNumberOfClusters()) t.SetRemoval(10);
1670 Int_t AliTPCtrackerMI::FollowToNextFast(AliTPCseed& t, Int_t nr) {
1671 //-----------------------------------------------------------------
1672 // This function tries to find a track prolongation to next pad row
1673 //-----------------------------------------------------------------
1675 Double_t x= GetXrow(nr), ymax=GetMaxY(nr);
1677 if (!t.GetProlongation(x,y,z)) {
1683 if (TMath::Abs(y)>ymax){
1686 t.SetRelativeSector((t.GetRelativeSector()+1) % fN);
1687 if (!t.Rotate(fSectors->GetAlpha()))
1689 } else if (y <-ymax) {
1690 t.SetRelativeSector((t.GetRelativeSector()-1+fN) % fN);
1691 if (!t.Rotate(-fSectors->GetAlpha()))
1694 if (!t.PropagateTo(x)) {
1697 t.GetProlongation(x,y,z);
1700 // update current shape info every 2 pad-row
1701 if ( (nr%2==0) || t.GetNumberOfClusters()<2 || (t.GetCurrentSigmaY2()<0.0001) ){
1702 // t.fCurrentSigmaY = GetSigmaY(&t);
1703 //t.fCurrentSigmaZ = GetSigmaZ(&t);
1707 AliTPCclusterMI *cl=0;
1712 const AliTPCtrackerRow &krow=GetRow(t.GetRelativeSector(),nr);
1713 if ( (t.GetSigmaY2()<0) || t.GetSigmaZ2()<0) return 0;
1715 Double_t roadz = 1.;
1718 if (TMath::Abs(TMath::Abs(y)-ymax)<krow.GetDeadZone()){
1720 t.SetClusterIndex2(row,-1);
1725 if (TMath::Abs(z)>(AliTPCReconstructor::GetCtgRange()*x+10)) t.SetClusterIndex2(row,-1);
1729 if ((cl==0)&&(krow)) {
1730 // cl = krow.FindNearest2(y+10,z,roady,roadz,index);
1731 cl = krow.FindNearest2(y,z,roady,roadz,index);
1733 if (cl) t.SetCurrentClusterIndex1(krow.GetIndex(index));
1737 t.SetCurrentCluster(cl);
1738 // Int_t accept = AcceptCluster(&t,t.fCurrentCluster);
1740 t.SetClusterIndex2(row,index);
1741 t.SetClusterPointer(row, cl);
1749 //_________________________________________________________________________
1750 Bool_t AliTPCtrackerMI::GetTrackPoint(Int_t index, AliTrackPoint &p ) const
1752 // Get track space point by index
1753 // return false in case the cluster doesn't exist
1754 AliTPCclusterMI *cl = GetClusterMI(index);
1755 if (!cl) return kFALSE;
1756 Int_t sector = (index&0xff000000)>>24;
1757 // Int_t row = (index&0x00ff0000)>>16;
1759 // xyz[0] = fParam->GetPadRowRadii(sector,row);
1760 xyz[0] = cl->GetX();
1761 xyz[1] = cl->GetY();
1762 xyz[2] = cl->GetZ();
1764 fParam->AdjustCosSin(sector,cos,sin);
1765 Float_t x = cos*xyz[0]-sin*xyz[1];
1766 Float_t y = cos*xyz[1]+sin*xyz[0];
1768 Float_t sigmaY2 = 0.027*cl->GetSigmaY2();
1769 if (sector < fParam->GetNInnerSector()) sigmaY2 *= 2.07;
1770 Float_t sigmaZ2 = 0.066*cl->GetSigmaZ2();
1771 if (sector < fParam->GetNInnerSector()) sigmaZ2 *= 1.77;
1772 cov[0] = sin*sin*sigmaY2;
1773 cov[1] = -sin*cos*sigmaY2;
1775 cov[3] = cos*cos*sigmaY2;
1778 p.SetXYZ(x,y,xyz[2],cov);
1779 AliGeomManager::ELayerID iLayer;
1781 if (sector < fParam->GetNInnerSector()) {
1782 iLayer = AliGeomManager::kTPC1;
1786 iLayer = AliGeomManager::kTPC2;
1787 idet = sector - fParam->GetNInnerSector();
1789 UShort_t volid = AliGeomManager::LayerToVolUID(iLayer,idet);
1790 p.SetVolumeID(volid);
1796 Int_t AliTPCtrackerMI::UpdateClusters(AliTPCseed& t, Int_t nr) {
1797 //-----------------------------------------------------------------
1798 // This function tries to find a track prolongation to next pad row
1799 //-----------------------------------------------------------------
1800 t.SetCurrentCluster(0);
1801 t.SetCurrentClusterIndex1(0);
1803 Double_t xt=t.GetX();
1804 Int_t row = GetRowNumber(xt)-1;
1805 Double_t ymax= GetMaxY(nr);
1807 if (row < nr) return 1; // don't prolongate if not information until now -
1808 // if (TMath::Abs(t.GetSnp())>0.9 && t.GetNumberOfClusters()>40. && fIteration!=2) {
1810 // return 0; // not prolongate strongly inclined tracks
1812 // if (TMath::Abs(t.GetSnp())>0.95) {
1814 // return 0; // not prolongate strongly inclined tracks
1815 // }// patch 28 fev 06
1817 Double_t x= GetXrow(nr);
1819 //t.PropagateTo(x+0.02);
1820 //t.PropagateTo(x+0.01);
1821 if (!t.PropagateTo(x)){
1828 if (TMath::Abs(y)>ymax){
1830 t.SetRelativeSector((t.GetRelativeSector()+1) % fN);
1831 if (!t.Rotate(fSectors->GetAlpha()))
1833 } else if (y <-ymax) {
1834 t.SetRelativeSector((t.GetRelativeSector()-1+fN) % fN);
1835 if (!t.Rotate(-fSectors->GetAlpha()))
1838 // if (!t.PropagateTo(x)){
1845 if (TMath::Abs(t.GetSnp())>AliTPCReconstructor::GetMaxSnpTracker()) return 0;
1847 if (!IsActive(t.GetRelativeSector(),nr)) {
1849 t.SetClusterIndex2(nr,-1);
1852 //AliInfo(Form("A - Sector%d phi %f - alpha %f", t.fRelativeSector,y/x, t.GetAlpha()));
1854 AliTPCtrackerRow &krow=GetRow(t.GetRelativeSector(),nr);
1856 if (TMath::Abs(TMath::Abs(y)-ymax)<krow.GetDeadZone()){
1858 t.SetClusterIndex2(nr,-1);
1863 if (TMath::Abs(t.GetZ())<(AliTPCReconstructor::GetCtgRange()*t.GetX()+10) && (TMath::Abs(t.GetSnp())<AliTPCReconstructor::GetMaxSnpTracker())) t.SetNFoundable(t.GetNFoundable()+1);
1869 if ( (nr%2==0) || t.GetNumberOfClusters()<2){
1870 // t.fCurrentSigmaY = GetSigmaY(&t);
1871 //t.fCurrentSigmaZ = GetSigmaZ(&t);
1875 AliTPCclusterMI *cl=0;
1878 Double_t roady = 1.;
1879 Double_t roadz = 1.;
1883 index = t.GetClusterIndex2(nr);
1884 if ( (index>0) && (index&0x8000)==0){
1885 cl = t.GetClusterPointer(nr);
1886 if ( (cl==0) && (index>0)) cl = GetClusterMI(index);
1887 t.SetCurrentClusterIndex1(index);
1889 t.SetCurrentCluster(cl);
1895 // if (index<0) return 0;
1896 UInt_t uindex = TMath::Abs(index);
1899 //cl = krow.FindNearest2(y+10,z,roady,roadz,uindex);
1900 cl = krow.FindNearest2(y,z,roady,roadz,uindex);
1903 if (cl) t.SetCurrentClusterIndex1(krow.GetIndex(uindex));
1904 t.SetCurrentCluster(cl);
1910 Int_t AliTPCtrackerMI::FollowToNextCluster(AliTPCseed & t, Int_t nr) {
1911 //-----------------------------------------------------------------
1912 // This function tries to find a track prolongation to next pad row
1913 //-----------------------------------------------------------------
1915 //update error according neighborhoud
1917 if (t.GetCurrentCluster()) {
1919 Int_t accept = AcceptCluster(&t,t.GetCurrentCluster());
1921 if (t.GetCurrentCluster()->IsUsed(10)){
1926 t.SetNShared(t.GetNShared()+1);
1927 if (t.GetNShared()>0.7*t.GetNumberOfClusters()) {
1932 if (fIteration>0) accept = 0;
1933 if (accept<3) UpdateTrack(&t,accept);
1937 if ( ( (t.GetSigmaY2()+t.GetSigmaZ2())>0.16)&& t.GetNumberOfClusters()>18) t.SetRemoval(10);
1938 if ( t.GetChi2()/t.GetNumberOfClusters()>6 &&t.GetNumberOfClusters()>18) t.SetRemoval(10);
1940 if (( (t.GetNFoundable()*0.5 > t.GetNumberOfClusters()) || t.GetNoCluster()>15)) t.SetRemoval(10);
1948 //_____________________________________________________________________________
1949 Int_t AliTPCtrackerMI::FollowProlongation(AliTPCseed& t, Int_t rf, Int_t step) {
1950 //-----------------------------------------------------------------
1951 // This function tries to find a track prolongation.
1952 //-----------------------------------------------------------------
1953 Double_t xt=t.GetX();
1955 Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
1956 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
1957 if (alpha < 0. ) alpha += 2.*TMath::Pi();
1959 t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN);
1961 Int_t first = GetRowNumber(xt)-1;
1962 for (Int_t nr= first; nr>=rf; nr-=step) {
1964 if (t.GetKinkIndexes()[0]>0){
1965 for (Int_t i=0;i<3;i++){
1966 Int_t index = t.GetKinkIndexes()[i];
1967 if (index==0) break;
1968 if (index<0) continue;
1970 AliKink * kink = (AliKink*)fEvent->GetKink(index-1);
1972 printf("PROBLEM\n");
1975 Int_t kinkrow = kink->GetTPCRow0()+2+Int_t(0.5/(0.05+kink->GetAngle(2)));
1977 AliExternalTrackParam paramd(t);
1978 kink->SetDaughter(paramd);
1979 kink->SetStatus(2,5);
1986 if (nr==80) t.UpdateReference();
1987 if (nr<fInnerSec->GetNRows())
1988 fSectors = fInnerSec;
1990 fSectors = fOuterSec;
1991 if (FollowToNext(t,nr)==0)
2000 //_____________________________________________________________________________
2001 Int_t AliTPCtrackerMI::FollowProlongationFast(AliTPCseed& t, Int_t rf, Int_t step) {
2002 //-----------------------------------------------------------------
2003 // This function tries to find a track prolongation.
2004 //-----------------------------------------------------------------
2005 Double_t xt=t.GetX();
2007 Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
2008 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
2009 if (alpha < 0. ) alpha += 2.*TMath::Pi();
2010 t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN);
2012 for (Int_t nr=GetRowNumber(xt)-1; nr>=rf; nr-=step) {
2014 if (FollowToNextFast(t,nr)==0)
2015 if (!t.IsActive()) return 0;
2025 Int_t AliTPCtrackerMI::FollowBackProlongation(AliTPCseed& t, Int_t rf) {
2026 //-----------------------------------------------------------------
2027 // This function tries to find a track prolongation.
2028 //-----------------------------------------------------------------
2030 Double_t xt=t.GetX();
2031 Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
2032 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
2033 if (alpha < 0. ) alpha += 2.*TMath::Pi();
2034 t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN);
2036 Int_t first = t.GetFirstPoint();
2037 if (first<GetRowNumber(xt)+1) first = GetRowNumber(xt)+1;
2039 if (first<0) first=0;
2040 for (Int_t nr=first; nr<=rf; nr++) {
2041 // if ( (TMath::Abs(t.GetSnp())>0.95)) break;//patch 28 fev 06
2042 if (t.GetKinkIndexes()[0]<0){
2043 for (Int_t i=0;i<3;i++){
2044 Int_t index = t.GetKinkIndexes()[i];
2045 if (index==0) break;
2046 if (index>0) continue;
2047 index = TMath::Abs(index);
2048 AliKink * kink = (AliKink*)fEvent->GetKink(index-1);
2050 printf("PROBLEM\n");
2053 Int_t kinkrow = kink->GetTPCRow0()-2-Int_t(0.5/(0.05+kink->GetAngle(2)));
2055 AliExternalTrackParam paramm(t);
2056 kink->SetMother(paramm);
2057 kink->SetStatus(2,1);
2064 if (nr<fInnerSec->GetNRows())
2065 fSectors = fInnerSec;
2067 fSectors = fOuterSec;
2078 Float_t AliTPCtrackerMI::OverlapFactor(AliTPCseed * s1, AliTPCseed * s2, Int_t &sum1, Int_t & sum2)
2086 Float_t dz2 =(s1->GetZ() - s2->GetZ());
2089 Float_t dy2 =TMath::Abs((s1->GetY() - s2->GetY()));
2091 Float_t distance = TMath::Sqrt(dz2+dy2);
2092 if (distance>4.) return 0; // if there are far away - not overlap - to reduce combinatorics
2095 Int_t firstpoint = TMath::Min(s1->GetFirstPoint(),s2->GetFirstPoint());
2096 Int_t lastpoint = TMath::Max(s1->GetLastPoint(),s2->GetLastPoint());
2101 if (firstpoint>lastpoint) {
2102 firstpoint =lastpoint;
2107 for (Int_t i=firstpoint-1;i<lastpoint+1;i++){
2108 if (s1->GetClusterIndex2(i)>0) sum1++;
2109 if (s2->GetClusterIndex2(i)>0) sum2++;
2110 if (s1->GetClusterIndex2(i)==s2->GetClusterIndex2(i) && s1->GetClusterIndex2(i)>0) {
2114 if (sum<5) return 0;
2116 Float_t summin = TMath::Min(sum1+1,sum2+1);
2117 Float_t ratio = (sum+1)/Float_t(summin);
2121 void AliTPCtrackerMI::SignShared(AliTPCseed * s1, AliTPCseed * s2)
2125 Float_t thetaCut = 0.2;//+10.*TMath::Sqrt(s1->GetSigmaTglZ()+ s2->GetSigmaTglZ());
2126 if (TMath::Abs(s1->GetTgl()-s2->GetTgl())>thetaCut) return;
2127 Float_t minCl = TMath::Min(s1->GetNumberOfClusters(),s2->GetNumberOfClusters());
2128 Int_t cutN0 = TMath::Max(5,TMath::Nint(0.1*minCl));
2133 //Int_t firstpoint = TMath::Max(s1->GetFirstPoint(),s2->GetFirstPoint());
2134 //Int_t lastpoint = TMath::Min(s1->GetLastPoint(),s2->GetLastPoint());
2135 Int_t firstpoint = 0;
2136 Int_t lastpoint = 160;
2138 // if (firstpoint>=lastpoint-5) return;;
2140 for (Int_t i=firstpoint;i<lastpoint;i++){
2141 // if ( (s1->GetClusterIndex2(i)&0xFFFF8FFF)==(s2->GetClusterIndex2(i)&0xFFFF8FFF) && s1->GetClusterIndex2(i)>0) {
2142 if ( (s1->GetClusterIndex2(i))==(s2->GetClusterIndex2(i)) && s1->GetClusterIndex2(i)>0) {
2144 s1->SetSharedMapBit(i, kTRUE);
2145 s2->SetSharedMapBit(i, kTRUE);
2147 if (s1->GetClusterIndex2(i)>0)
2148 s1->SetClusterMapBit(i, kTRUE);
2150 if (sumshared>cutN0){
2153 for (Int_t i=firstpoint;i<lastpoint;i++){
2154 // if ( (s1->GetClusterIndex2(i)&0xFFFF8FFF)==(s2->GetClusterIndex2(i)&0xFFFF8FFF) && s1->GetClusterIndex2(i)>0) {
2155 if ( (s1->GetClusterIndex2(i))==(s2->GetClusterIndex2(i)) && s1->GetClusterIndex2(i)>0) {
2156 AliTPCTrackerPoint *p1 = s1->GetTrackPoint(i);
2157 AliTPCTrackerPoint *p2 = s2->GetTrackPoint(i);;
2158 if (s1->IsActive()&&s2->IsActive()){
2159 p1->SetShared(kTRUE);
2160 p2->SetShared(kTRUE);
2166 if (sumshared>cutN0){
2167 for (Int_t i=0;i<4;i++){
2168 if (s1->GetOverlapLabel(3*i)==0){
2169 s1->SetOverlapLabel(3*i, s2->GetLabel());
2170 s1->SetOverlapLabel(3*i+1,sumshared);
2171 s1->SetOverlapLabel(3*i+2,s2->GetUniqueID());
2175 for (Int_t i=0;i<4;i++){
2176 if (s2->GetOverlapLabel(3*i)==0){
2177 s2->SetOverlapLabel(3*i, s1->GetLabel());
2178 s2->SetOverlapLabel(3*i+1,sumshared);
2179 s2->SetOverlapLabel(3*i+2,s1->GetUniqueID());
2186 void AliTPCtrackerMI::SignShared(TObjArray * arr)
2189 //sort trackss according sectors
2191 for (Int_t i=0; i<arr->GetEntriesFast(); i++) {
2192 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2194 //if (pt) RotateToLocal(pt);
2198 arr->Sort(); // sorting according relative sectors
2199 arr->Expand(arr->GetEntries());
2202 Int_t nseed=arr->GetEntriesFast();
2203 for (Int_t i=0; i<nseed; i++) {
2204 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2206 for (Int_t j=0;j<=12;j++){
2207 pt->SetOverlapLabel(j,0);
2210 for (Int_t i=0; i<nseed; i++) {
2211 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2213 if (pt->GetRemoval()>10) continue;
2214 for (Int_t j=i+1; j<nseed; j++){
2215 AliTPCseed *pt2=(AliTPCseed*)arr->UncheckedAt(j);
2216 if (TMath::Abs(pt->GetRelativeSector()-pt2->GetRelativeSector())>1) continue;
2218 if (pt2->GetRemoval()<=10) {
2219 //if ( TMath::Abs(pt->GetRelativeSector()-pt2->GetRelativeSector())>0) break;
2227 void AliTPCtrackerMI::SortTracks(TObjArray * arr, Int_t mode) const
2230 //sort tracks in array according mode criteria
2231 Int_t nseed = arr->GetEntriesFast();
2232 for (Int_t i=0; i<nseed; i++) {
2233 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2244 void AliTPCtrackerMI::RemoveUsed2(TObjArray * arr, Float_t factor1, Float_t factor2, Int_t minimal)
2247 // Loop over all tracks and remove overlaped tracks (with lower quality)
2249 // 1. Unsign clusters
2250 // 2. Sort tracks according quality
2251 // Quality is defined by the number of cluster between first and last points
2253 // 3. Loop over tracks - decreasing quality order
2254 // a.) remove - If the fraction of shared cluster less than factor (1- n or 2)
2255 // b.) remove - If the minimal number of clusters less than minimal and not ITS
2256 // c.) if track accepted - sign clusters
2258 //Called in - AliTPCtrackerMI::Clusters2Tracks()
2259 // - AliTPCtrackerMI::PropagateBack()
2260 // - AliTPCtrackerMI::RefitInward()
2266 Int_t nseed = arr->GetEntriesFast();
2267 Float_t * quality = new Float_t[nseed];
2268 Int_t * indexes = new Int_t[nseed];
2272 for (Int_t i=0; i<nseed; i++) {
2273 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2278 pt->UpdatePoints(); //select first last max dens points
2279 Float_t * points = pt->GetPoints();
2280 if (points[3]<0.8) quality[i] =-1;
2281 quality[i] = (points[2]-points[0])+pt->GetNumberOfClusters();
2282 //prefer high momenta tracks if overlaps
2283 quality[i] *= TMath::Sqrt(TMath::Abs(pt->Pt())+0.5);
2285 TMath::Sort(nseed,quality,indexes);
2288 for (Int_t itrack=0; itrack<nseed; itrack++) {
2289 Int_t trackindex = indexes[itrack];
2290 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(trackindex);
2293 if (quality[trackindex]<0){
2295 delete arr->RemoveAt(trackindex);
2298 arr->RemoveAt(trackindex);
2304 Int_t first = Int_t(pt->GetPoints()[0]);
2305 Int_t last = Int_t(pt->GetPoints()[2]);
2306 Double_t factor = (pt->GetBConstrain()) ? factor1: factor2;
2308 Int_t found,foundable,shared;
2309 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
2310 // pt->GetClusterStatistic(0,160, found, foundable,shared,kFALSE);
2311 Bool_t itsgold =kFALSE;
2314 if (pt->GetESD()->GetITSclusters(dummy)>4) itsgold= kTRUE;
2318 if (Float_t(shared+1)/Float_t(found+1)>factor){
2319 if (pt->GetKinkIndexes()[0]!=0) continue; //don't remove tracks - part of the kinks
2320 delete arr->RemoveAt(trackindex);
2323 if (pt->GetNumberOfClusters()<50&&(found-0.5*shared)<minimal){ //remove short tracks
2324 if (pt->GetKinkIndexes()[0]!=0) continue; //don't remove tracks - part of the kinks
2325 delete arr->RemoveAt(trackindex);
2331 //if (sharedfactor>0.4) continue;
2332 if (pt->GetKinkIndexes()[0]>0) continue;
2333 //Remove tracks with undefined properties - seems
2334 if (pt->GetSigmaY2()<kAlmost0) continue; // ? what is the origin ?
2336 for (Int_t i=first; i<last; i++) {
2337 Int_t index=pt->GetClusterIndex2(i);
2338 // if (index<0 || index&0x8000 ) continue;
2339 if (index<0 || index&0x8000 ) continue;
2340 AliTPCclusterMI *c= pt->GetClusterPointer(i);
2347 Info("RemoveUsed","\n*****\nNumber of good tracks after shared removal\t%d\n",fNtracks);
2353 void AliTPCtrackerMI::UnsignClusters()
2356 // loop over all clusters and unsign them
2359 for (Int_t sec=0;sec<fkNIS;sec++){
2360 for (Int_t row=0;row<fInnerSec->GetNRows();row++){
2361 AliTPCclusterMI *cl = fInnerSec[sec][row].GetClusters1();
2362 for (Int_t icl =0;icl< fInnerSec[sec][row].GetN1();icl++)
2363 // if (cl[icl].IsUsed(10))
2365 cl = fInnerSec[sec][row].GetClusters2();
2366 for (Int_t icl =0;icl< fInnerSec[sec][row].GetN2();icl++)
2367 //if (cl[icl].IsUsed(10))
2372 for (Int_t sec=0;sec<fkNOS;sec++){
2373 for (Int_t row=0;row<fOuterSec->GetNRows();row++){
2374 AliTPCclusterMI *cl = fOuterSec[sec][row].GetClusters1();
2375 for (Int_t icl =0;icl< fOuterSec[sec][row].GetN1();icl++)
2376 //if (cl[icl].IsUsed(10))
2378 cl = fOuterSec[sec][row].GetClusters2();
2379 for (Int_t icl =0;icl< fOuterSec[sec][row].GetN2();icl++)
2380 //if (cl[icl].IsUsed(10))
2389 void AliTPCtrackerMI::SignClusters(TObjArray * arr, Float_t fnumber, Float_t fdensity)
2392 //sign clusters to be "used"
2394 // snumber and sdensity sign number of sigmas - bellow mean value to be accepted
2395 // loop over "primaries"
2409 Int_t nseed = arr->GetEntriesFast();
2410 for (Int_t i=0; i<nseed; i++) {
2411 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2415 if (!(pt->IsActive())) continue;
2416 Float_t dens = pt->GetNumberOfClusters()/Float_t(pt->GetNFoundable());
2417 if ( (dens>0.7) && (pt->GetNumberOfClusters()>70)){
2419 sumdens2+= dens*dens;
2420 sumn += pt->GetNumberOfClusters();
2421 sumn2 += pt->GetNumberOfClusters()*pt->GetNumberOfClusters();
2422 Float_t chi2 = pt->GetChi2()/pt->GetNumberOfClusters();
2425 sumchi2 +=chi2*chi2;
2430 Float_t mdensity = 0.9;
2431 Float_t meann = 130;
2432 Float_t meanchi = 1;
2433 Float_t sdensity = 0.1;
2434 Float_t smeann = 10;
2435 Float_t smeanchi =0.4;
2439 mdensity = sumdens/sum;
2441 meanchi = sumchi/sum;
2443 sdensity = sumdens2/sum-mdensity*mdensity;
2445 sdensity = TMath::Sqrt(sdensity);
2449 smeann = sumn2/sum-meann*meann;
2451 smeann = TMath::Sqrt(smeann);
2455 smeanchi = sumchi2/sum - meanchi*meanchi;
2457 smeanchi = TMath::Sqrt(smeanchi);
2463 //REMOVE SHORT DELTAS or tracks going out of sensitive volume of TPC
2465 for (Int_t i=0; i<nseed; i++) {
2466 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2470 if (pt->GetBSigned()) continue;
2471 if (pt->GetBConstrain()) continue;
2472 //if (!(pt->IsActive())) continue;
2474 Int_t found,foundable,shared;
2475 pt->GetClusterStatistic(0,160,found, foundable,shared);
2476 if (shared/float(found)>0.3) {
2477 if (shared/float(found)>0.9 ){
2478 //delete arr->RemoveAt(i);
2483 Bool_t isok =kFALSE;
2484 if ( (pt->GetNShared()/pt->GetNumberOfClusters()<0.5) &&pt->GetNumberOfClusters()>60)
2486 if ((TMath::Abs(1/pt->GetC())<100.) && (pt->GetNShared()/pt->GetNumberOfClusters()<0.7))
2488 if (TMath::Abs(pt->GetZ()/pt->GetX())>1.1)
2490 if ( (TMath::Abs(pt->GetSnp()>0.7) && pt->GetD(0,0)>60.))
2494 for (Int_t i=0; i<160; i++) {
2495 Int_t index=pt->GetClusterIndex2(i);
2496 if (index<0) continue;
2497 AliTPCclusterMI *c= pt->GetClusterPointer(i);
2499 //if (!(c->IsUsed(10))) c->Use();
2506 Double_t maxchi = meanchi+2.*smeanchi;
2508 for (Int_t i=0; i<nseed; i++) {
2509 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2513 //if (!(pt->IsActive())) continue;
2514 if (pt->GetBSigned()) continue;
2515 Double_t chi = pt->GetChi2()/pt->GetNumberOfClusters();
2516 if (chi>maxchi) continue;
2519 Float_t dens = pt->GetNumberOfClusters()/Float_t(pt->GetNFoundable());
2521 //sign only tracks with enoug big density at the beginning
2523 if ((pt->GetDensityFirst(40)<0.75) && pt->GetNumberOfClusters()<meann) continue;
2526 Double_t mindens = TMath::Max(double(mdensity-sdensity*fdensity*bfactor),0.65);
2527 Double_t minn = TMath::Max(Int_t(meann-fnumber*smeann*bfactor),50);
2529 // if (pt->fBConstrain) mindens = TMath::Max(mdensity-sdensity*fdensity*bfactor,0.65);
2530 if ( (pt->GetRemoval()==10) && (pt->GetSnp()>0.8)&&(dens>mindens))
2533 if ((dens>mindens && pt->GetNumberOfClusters()>minn) && chi<maxchi ){
2534 //Int_t noc=pt->GetNumberOfClusters();
2535 pt->SetBSigned(kTRUE);
2536 for (Int_t i=0; i<160; i++) {
2538 Int_t index=pt->GetClusterIndex2(i);
2539 if (index<0) continue;
2540 AliTPCclusterMI *c= pt->GetClusterPointer(i);
2542 // if (!(c->IsUsed(10))) c->Use();
2547 // gLastCheck = nseed;
2555 void AliTPCtrackerMI::StopNotActive(TObjArray * arr, Int_t row0, Float_t th0, Float_t th1, Float_t th2) const
2557 // stop not active tracks
2558 // take th1 as threshold for number of founded to number of foundable on last 10 active rows
2559 // take th2 as threshold for number of founded to number of foundable on last 20 active rows
2560 Int_t nseed = arr->GetEntriesFast();
2562 for (Int_t i=0; i<nseed; i++) {
2563 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2567 if (!(pt->IsActive())) continue;
2568 StopNotActive(pt,row0,th0, th1,th2);
2574 void AliTPCtrackerMI::StopNotActive(AliTPCseed * seed, Int_t row0, Float_t th0, Float_t th1,
2577 // stop not active tracks
2578 // take th1 as threshold for number of founded to number of foundable on last 10 active rows
2579 // take th2 as threshold for number of founded to number of foundable on last 20 active rows
2582 Int_t foundable = 0;
2583 Int_t maxindex = seed->GetLastPoint(); //last foundable row
2584 if (seed->GetNFoundable()*th0 > seed->GetNumberOfClusters()) {
2585 seed->Desactivate(10) ;
2589 for (Int_t i=row0; i<maxindex; i++){
2590 Int_t index = seed->GetClusterIndex2(i);
2591 if (index!=-1) foundable++;
2593 if (foundable<=30) sumgood1++;
2594 if (foundable<=50) {
2601 if (foundable>=30.){
2602 if (sumgood1<(th1*30.)) seed->Desactivate(10);
2605 if (sumgood2<(th2*50.)) seed->Desactivate(10);
2609 Int_t AliTPCtrackerMI::RefitInward(AliESDEvent *event)
2612 // back propagation of ESD tracks
2615 const Int_t kMaxFriendTracks=2000;
2619 //PrepareForProlongation(fSeeds,1);
2620 PropagateForward2(fSeeds);
2621 RemoveUsed2(fSeeds,0.4,0.4,20);
2623 TObjArray arraySeed(fSeeds->GetEntries());
2624 for (Int_t i=0;i<fSeeds->GetEntries();i++) {
2625 arraySeed.AddAt(fSeeds->At(i),i);
2627 SignShared(&arraySeed);
2628 // FindCurling(fSeeds, event,2); // find multi found tracks
2629 FindSplitted(fSeeds, event,2); // find multi found tracks
2632 Int_t nseed = fSeeds->GetEntriesFast();
2633 for (Int_t i=0;i<nseed;i++){
2634 AliTPCseed * seed = (AliTPCseed*) fSeeds->UncheckedAt(i);
2635 if (!seed) continue;
2636 if (seed->GetKinkIndex(0)>0) UpdateKinkQualityD(seed); // update quality informations for kinks
2638 seed->PropagateTo(fParam->GetInnerRadiusLow());
2639 seed->UpdatePoints();
2640 AddCovariance(seed);
2642 AliESDtrack *esd=event->GetTrack(i);
2643 seed->CookdEdx(0.02,0.6);
2644 CookLabel(seed,0.1); //For comparison only
2645 esd->SetTPCClusterMap(seed->GetClusterMap());
2646 esd->SetTPCSharedMap(seed->GetSharedMap());
2648 if (AliTPCReconstructor::StreamLevel()>1 && seed!=0&&esd!=0) {
2649 TTreeSRedirector &cstream = *fDebugStreamer;
2656 if (seed->GetNumberOfClusters()>15){
2657 esd->UpdateTrackParams(seed,AliESDtrack::kTPCrefit);
2658 esd->SetTPCPoints(seed->GetPoints());
2659 esd->SetTPCPointsF(seed->GetNFoundable());
2660 Int_t ndedx = seed->GetNCDEDX(0)+seed->GetNCDEDX(1)+seed->GetNCDEDX(2)+seed->GetNCDEDX(3);
2661 Float_t sdedx = (seed->GetSDEDX(0)+seed->GetSDEDX(1)+seed->GetSDEDX(2)+seed->GetSDEDX(3))*0.25;
2662 Float_t dedx = seed->GetdEdx();
2663 esd->SetTPCsignal(dedx, sdedx, ndedx);
2665 // add seed to the esd track in Calib level
2667 Bool_t storeFriend = gRandom->Rndm()<(kMaxFriendTracks)/Float_t(nseed);
2668 if (AliTPCReconstructor::StreamLevel()>0 &&storeFriend){
2669 AliTPCseed * seedCopy = new AliTPCseed(*seed, kTRUE);
2670 esd->AddCalibObject(seedCopy);
2675 //printf("problem\n");
2678 //FindKinks(fSeeds,event);
2679 Info("RefitInward","Number of refitted tracks %d",ntracks);
2684 Int_t AliTPCtrackerMI::PropagateBack(AliESDEvent *event)
2687 // back propagation of ESD tracks
2693 PropagateBack(fSeeds);
2694 RemoveUsed2(fSeeds,0.4,0.4,20);
2695 //FindCurling(fSeeds, fEvent,1);
2696 FindSplitted(fSeeds, event,1); // find multi found tracks
2699 Int_t nseed = fSeeds->GetEntriesFast();
2701 for (Int_t i=0;i<nseed;i++){
2702 AliTPCseed * seed = (AliTPCseed*) fSeeds->UncheckedAt(i);
2703 if (!seed) continue;
2704 if (seed->GetKinkIndex(0)<0) UpdateKinkQualityM(seed); // update quality informations for kinks
2705 seed->UpdatePoints();
2706 AddCovariance(seed);
2707 AliESDtrack *esd=event->GetTrack(i);
2708 seed->CookdEdx(0.02,0.6);
2709 CookLabel(seed,0.1); //For comparison only
2710 if (seed->GetNumberOfClusters()>15){
2711 esd->UpdateTrackParams(seed,AliESDtrack::kTPCout);
2712 esd->SetTPCPoints(seed->GetPoints());
2713 esd->SetTPCPointsF(seed->GetNFoundable());
2714 Int_t ndedx = seed->GetNCDEDX(0)+seed->GetNCDEDX(1)+seed->GetNCDEDX(2)+seed->GetNCDEDX(3);
2715 Float_t sdedx = (seed->GetSDEDX(0)+seed->GetSDEDX(1)+seed->GetSDEDX(2)+seed->GetSDEDX(3))*0.25;
2716 Float_t dedx = seed->GetdEdx();
2717 esd->SetTPCsignal(dedx, sdedx, ndedx);
2719 Int_t eventnumber = event->GetEventNumberInFile();// patch 28 fev 06
2720 // This is most likely NOT the event number you'd like to use. It has nothing to do with the 'real' event number
2721 if (AliTPCReconstructor::StreamLevel()>1 && esd) {
2722 (*fDebugStreamer)<<"Cback"<<
2725 "EventNrInFile="<<eventnumber<<
2730 //FindKinks(fSeeds,event);
2731 Info("PropagateBack","Number of back propagated tracks %d",ntracks);
2738 void AliTPCtrackerMI::DeleteSeeds()
2743 Int_t nseed = fSeeds->GetEntriesFast();
2744 for (Int_t i=0;i<nseed;i++){
2745 AliTPCseed * seed = (AliTPCseed*)fSeeds->At(i);
2746 if (seed) delete fSeeds->RemoveAt(i);
2753 void AliTPCtrackerMI::ReadSeeds(AliESDEvent *event, Int_t direction)
2756 //read seeds from the event
2758 Int_t nentr=event->GetNumberOfTracks();
2760 Info("ReadSeeds", "Number of ESD tracks: %d\n", nentr);
2765 fSeeds = new TObjArray(nentr);
2769 for (Int_t i=0; i<nentr; i++) {
2770 AliESDtrack *esd=event->GetTrack(i);
2771 ULong_t status=esd->GetStatus();
2772 if (!(status&AliESDtrack::kTPCin)) continue;
2773 AliTPCtrack t(*esd);
2774 t.SetNumberOfClusters(0);
2775 // AliTPCseed *seed = new AliTPCseed(t,t.GetAlpha());
2776 AliTPCseed *seed = new AliTPCseed(t/*,t.GetAlpha()*/);
2777 seed->SetUniqueID(esd->GetID());
2778 AddCovariance(seed); //add systematic ucertainty
2779 for (Int_t ikink=0;ikink<3;ikink++) {
2780 Int_t index = esd->GetKinkIndex(ikink);
2781 seed->GetKinkIndexes()[ikink] = index;
2782 if (index==0) continue;
2783 index = TMath::Abs(index);
2784 AliESDkink * kink = fEvent->GetKink(index-1);
2785 if (kink&&esd->GetKinkIndex(ikink)<0){
2786 if ((status & AliESDtrack::kTRDrefit) != 0) kink->SetStatus(1,2);
2787 if ((status & AliESDtrack::kITSout) != 0) kink->SetStatus(1,0);
2789 if (kink&&esd->GetKinkIndex(ikink)>0){
2790 if ((status & AliESDtrack::kTRDrefit) != 0) kink->SetStatus(1,6);
2791 if ((status & AliESDtrack::kITSout) != 0) kink->SetStatus(1,4);
2795 if (((status&AliESDtrack::kITSout)==0)&&(direction==1)) seed->ResetCovariance(10.);
2796 if ( direction ==2 &&(status & AliESDtrack::kTRDrefit) == 0 ) seed->ResetCovariance(10.);
2797 if ( direction ==2 && ((status & AliESDtrack::kTPCout) == 0) ) {
2802 if ( direction ==2 &&(status & AliESDtrack::kTRDrefit) > 0 ) {
2803 Double_t par0[5],par1[5],alpha,x;
2804 esd->GetInnerExternalParameters(alpha,x,par0);
2805 esd->GetExternalParameters(x,par1);
2806 Double_t delta1 = TMath::Abs(par0[4]-par1[4])/(0.000000001+TMath::Abs(par0[4]+par1[4]));
2807 Double_t delta2 = TMath::Abs(par0[3]-par1[3]);
2809 if (esd->GetTRDncls()>0) trdchi2 = esd->GetTRDchi2()/esd->GetTRDncls();
2810 //reset covariance if suspicious
2811 if ( (delta1>0.1) || (delta2>0.006) ||trdchi2>7.)
2812 seed->ResetCovariance(10.);
2817 // rotate to the local coordinate system
2819 fSectors=fInnerSec; fN=fkNIS;
2820 Double_t alpha=seed->GetAlpha() - fSectors->GetAlphaShift();
2821 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
2822 if (alpha < 0. ) alpha += 2.*TMath::Pi();
2823 Int_t ns=Int_t(alpha/fSectors->GetAlpha())%fN;
2824 alpha =ns*fSectors->GetAlpha() + fSectors->GetAlphaShift();
2825 if (alpha<-TMath::Pi()) alpha += 2*TMath::Pi();
2826 if (alpha>=TMath::Pi()) alpha -= 2*TMath::Pi();
2827 alpha-=seed->GetAlpha();
2828 if (!seed->Rotate(alpha)) {
2834 if (esd->GetKinkIndex(0)<=0){
2835 for (Int_t irow=0;irow<160;irow++){
2836 Int_t index = seed->GetClusterIndex2(irow);
2839 AliTPCclusterMI * cl = GetClusterMI(index);
2840 seed->SetClusterPointer(irow,cl);
2842 if ((index & 0x8000)==0){
2843 cl->Use(10); // accepted cluster
2845 cl->Use(6); // close cluster not accepted
2848 Info("ReadSeeds","Not found cluster");
2853 fSeeds->AddAt(seed,i);
2859 //_____________________________________________________________________________
2860 void AliTPCtrackerMI::MakeSeeds3(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, Float_t cuts[4],
2861 Float_t deltay, Int_t ddsec) {
2862 //-----------------------------------------------------------------
2863 // This function creates track seeds.
2864 // SEEDING WITH VERTEX CONSTRAIN
2865 //-----------------------------------------------------------------
2866 // cuts[0] - fP4 cut
2867 // cuts[1] - tan(phi) cut
2868 // cuts[2] - zvertex cut
2869 // cuts[3] - fP3 cut
2877 Double_t x[5], c[15];
2878 // Int_t di = i1-i2;
2880 AliTPCseed * seed = new AliTPCseed();
2881 Double_t alpha=fSectors->GetAlpha(), shift=fSectors->GetAlphaShift();
2882 Double_t cs=cos(alpha), sn=sin(alpha);
2884 // Double_t x1 =fOuterSec->GetX(i1);
2885 //Double_t xx2=fOuterSec->GetX(i2);
2887 Double_t x1 =GetXrow(i1);
2888 Double_t xx2=GetXrow(i2);
2890 Double_t x3=GetX(), y3=GetY(), z3=GetZ();
2892 Int_t imiddle = (i2+i1)/2; //middle pad row index
2893 Double_t xm = GetXrow(imiddle); // radius of middle pad-row
2894 const AliTPCtrackerRow& krm=GetRow(sec,imiddle); //middle pad -row
2898 const AliTPCtrackerRow& kr1=GetRow(ns,i1);
2899 Double_t ymax = GetMaxY(i1)-kr1.GetDeadZone()-1.5;
2900 Double_t ymaxm = GetMaxY(imiddle)-kr1.GetDeadZone()-1.5;
2903 // change cut on curvature if it can't reach this layer
2904 // maximal curvature set to reach it
2905 Double_t dvertexmax = TMath::Sqrt((x1-x3)*(x1-x3)+(ymax+5-y3)*(ymax+5-y3));
2906 if (dvertexmax*0.5*cuts[0]>0.85){
2907 cuts[0] = 0.85/(dvertexmax*0.5+1.);
2909 Double_t r2min = 1/(cuts[0]*cuts[0]); //minimal square of radius given by cut
2912 if (deltay>0) ddsec = 0;
2913 // loop over clusters
2914 for (Int_t is=0; is < kr1; is++) {
2916 if (kr1[is]->IsUsed(10)) continue;
2917 Double_t y1=kr1[is]->GetY(), z1=kr1[is]->GetZ();
2918 //if (TMath::Abs(y1)>ymax) continue;
2920 if (deltay>0 && TMath::Abs(ymax-TMath::Abs(y1))> deltay ) continue; // seed only at the edge
2922 // find possible directions
2923 Float_t anglez = (z1-z3)/(x1-x3);
2924 Float_t extraz = z1 - anglez*(x1-xx2); // extrapolated z
2927 //find rotation angles relative to line given by vertex and point 1
2928 Double_t dvertex2 = (x1-x3)*(x1-x3)+(y1-y3)*(y1-y3);
2929 Double_t dvertex = TMath::Sqrt(dvertex2);
2930 Double_t angle13 = TMath::ATan((y1-y3)/(x1-x3));
2931 Double_t cs13 = cos(-angle13), sn13 = sin(-angle13);
2934 // loop over 2 sectors
2940 Double_t dddz1=0; // direction of delta inclination in z axis
2947 for (Int_t dsec = dsec1; dsec<=dsec2;dsec++){
2948 Int_t sec2 = sec + dsec;
2950 // AliTPCtrackerRow& kr2 = fOuterSec[(sec2+fkNOS)%fkNOS][i2];
2951 //AliTPCtrackerRow& kr2m = fOuterSec[(sec2+fkNOS)%fkNOS][imiddle];
2952 AliTPCtrackerRow& kr2 = GetRow((sec2+fkNOS)%fkNOS,i2);
2953 AliTPCtrackerRow& kr2m = GetRow((sec2+fkNOS)%fkNOS,imiddle);
2954 Int_t index1 = TMath::Max(kr2.Find(extraz-0.6-dddz1*TMath::Abs(z1)*0.05)-1,0);
2955 Int_t index2 = TMath::Min(kr2.Find(extraz+0.6+dddz2*TMath::Abs(z1)*0.05)+1,kr2);
2957 // rotation angles to p1-p3
2958 Double_t cs13r = cos(-angle13+dsec*alpha)/dvertex, sn13r = sin(-angle13+dsec*alpha)/dvertex;
2959 Double_t x2, y2, z2;
2961 // Double_t dymax = maxangle*TMath::Abs(x1-xx2);
2964 Double_t dxx0 = (xx2-x3)*cs13r;
2965 Double_t dyy0 = (xx2-x3)*sn13r;
2966 for (Int_t js=index1; js < index2; js++) {
2967 const AliTPCclusterMI *kcl = kr2[js];
2968 if (kcl->IsUsed(10)) continue;
2970 //calcutate parameters
2972 Double_t yy0 = dyy0 +(kcl->GetY()-y3)*cs13r;
2974 if (TMath::Abs(yy0)<0.000001) continue;
2975 Double_t xx0 = dxx0 -(kcl->GetY()-y3)*sn13r;
2976 Double_t y0 = 0.5*(xx0*xx0+yy0*yy0-xx0)/yy0;
2977 Double_t r02 = (0.25+y0*y0)*dvertex2;
2978 //curvature (radius) cut
2979 if (r02<r2min) continue;
2983 Double_t c0 = 1/TMath::Sqrt(r02);
2987 //Double_t dfi0 = 2.*TMath::ASin(dvertex*c0*0.5);
2988 //Double_t dfi1 = 2.*TMath::ASin(TMath::Sqrt(yy0*yy0+(1-xx0)*(1-xx0))*dvertex*c0*0.5);
2989 Double_t dfi0 = 2.*AliTPCFastMath::FastAsin(dvertex*c0*0.5);
2990 Double_t dfi1 = 2.*AliTPCFastMath::FastAsin(TMath::Sqrt(yy0*yy0+(1-xx0)*(1-xx0))*dvertex*c0*0.5);
2993 Double_t z0 = kcl->GetZ();
2994 Double_t zzzz2 = z1-(z1-z3)*dfi1/dfi0;
2995 if (TMath::Abs(zzzz2-z0)>0.5) continue;
2998 Double_t dip = (z1-z0)*c0/dfi1;
2999 Double_t x0 = (0.5*cs13+y0*sn13)*dvertex*c0;
3010 x2= xx2*cs-y2*sn*dsec;
3011 y2=+xx2*sn*dsec+y2*cs;
3021 // do we have cluster at the middle ?
3023 GetProlongation(x1,xm,x,ym,zm);
3025 AliTPCclusterMI * cm=0;
3026 if (TMath::Abs(ym)-ymaxm<0){
3027 cm = krm.FindNearest2(ym,zm,1.0,0.6,dummy);
3028 if ((!cm) || (cm->IsUsed(10))) {
3033 // rotate y1 to system 0
3034 // get state vector in rotated system
3035 Double_t yr1 = (-0.5*sn13+y0*cs13)*dvertex*c0;
3036 Double_t xr2 = x0*cs+yr1*sn*dsec;
3037 Double_t xr[5]={kcl->GetY(),kcl->GetZ(), xr2, dip, c0};
3039 GetProlongation(xx2,xm,xr,ym,zm);
3040 if (TMath::Abs(ym)-ymaxm<0){
3041 cm = kr2m.FindNearest2(ym,zm,1.0,0.6,dummy);
3042 if ((!cm) || (cm->IsUsed(10))) {
3052 dym = ym - cm->GetY();
3053 dzm = zm - cm->GetZ();
3060 Double_t sy1=kr1[is]->GetSigmaY2()*2., sz1=kr1[is]->GetSigmaZ2()*2.;
3061 Double_t sy2=kcl->GetSigmaY2()*2., sz2=kcl->GetSigmaZ2()*2.;
3062 //Double_t sy3=400*3./12., sy=0.1, sz=0.1;
3063 Double_t sy3=25000*x[4]*x[4]+0.1, sy=0.1, sz=0.1;
3064 //Double_t sy3=25000*x[4]*x[4]*60+0.5, sy=0.1, sz=0.1;
3066 Double_t f40=(F1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
3067 Double_t f42=(F1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
3068 Double_t f43=(F1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
3069 Double_t f20=(F2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
3070 Double_t f22=(F2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
3071 Double_t f23=(F2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
3073 Double_t f30=(F3(x1,y1+sy,x2,y2,z1,z2)-x[3])/sy;
3074 Double_t f31=(F3(x1,y1,x2,y2,z1+sz,z2)-x[3])/sz;
3075 Double_t f32=(F3(x1,y1,x2,y2+sy,z1,z2)-x[3])/sy;
3076 Double_t f34=(F3(x1,y1,x2,y2,z1,z2+sz)-x[3])/sz;
3080 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3081 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3082 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3083 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3084 c[13]=f30*sy1*f40+f32*sy2*f42;
3085 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3087 // if (!BuildSeed(kr1[is],kcl,0,x1,x2,x3,x,c)) continue;
3089 UInt_t index=kr1.GetIndex(is);
3090 seed->~AliTPCseed(); // this does not set the pointer to 0...
3091 AliTPCseed *track=new(seed) AliTPCseed(x1, ns*alpha+shift, x, c, index);
3093 track->SetIsSeeding(kTRUE);
3094 track->SetSeed1(i1);
3095 track->SetSeed2(i2);
3096 track->SetSeedType(3);
3100 FollowProlongation(*track, (i1+i2)/2,1);
3101 Int_t foundable,found,shared;
3102 track->GetClusterStatistic((i1+i2)/2,i1, found, foundable, shared, kTRUE);
3103 if ((found<0.55*foundable) || shared>0.5*found || (track->GetSigmaY2()+track->GetSigmaZ2())>0.5){
3105 seed->~AliTPCseed();
3111 FollowProlongation(*track, i2,1);
3115 track->SetBConstrain(1);
3116 // track->fLastPoint = i1+fInnerSec->GetNRows(); // first cluster in track position
3117 track->SetLastPoint(i1); // first cluster in track position
3118 track->SetFirstPoint(track->GetLastPoint());
3120 if (track->GetNumberOfClusters()<(i1-i2)*0.5 ||
3121 track->GetNumberOfClusters() < track->GetNFoundable()*0.6 ||
3122 track->GetNShared()>0.4*track->GetNumberOfClusters() ) {
3124 seed->~AliTPCseed();
3128 // Z VERTEX CONDITION
3129 Double_t zv, bz=GetBz();
3130 if ( !track->GetZAt(0.,bz,zv) ) continue;
3131 if (TMath::Abs(zv-z3)>cuts[2]) {
3132 FollowProlongation(*track, TMath::Max(i2-20,0));
3133 if ( !track->GetZAt(0.,bz,zv) ) continue;
3134 if (TMath::Abs(zv-z3)>cuts[2]){
3135 FollowProlongation(*track, TMath::Max(i2-40,0));
3136 if ( !track->GetZAt(0.,bz,zv) ) continue;
3137 if (TMath::Abs(zv-z3)>cuts[2] &&(track->GetNumberOfClusters() > track->GetNFoundable()*0.7)){
3138 // make seed without constrain
3139 AliTPCseed * track2 = MakeSeed(track,0.2,0.5,1.);
3140 FollowProlongation(*track2, i2,1);
3141 track2->SetBConstrain(kFALSE);
3142 track2->SetSeedType(1);
3143 arr->AddLast(track2);
3145 seed->~AliTPCseed();
3150 seed->~AliTPCseed();
3157 track->SetSeedType(0);
3158 arr->AddLast(track);
3159 seed = new AliTPCseed;
3161 // don't consider other combinations
3162 if (track->GetNumberOfClusters() > track->GetNFoundable()*0.8)
3168 Info("MakeSeeds3","\nSeeding statistic:\t%d\t%d\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin,nout1,nout2);
3174 void AliTPCtrackerMI::MakeSeeds5(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, Float_t cuts[4],
3179 //-----------------------------------------------------------------
3180 // This function creates track seeds.
3181 //-----------------------------------------------------------------
3182 // cuts[0] - fP4 cut
3183 // cuts[1] - tan(phi) cut
3184 // cuts[2] - zvertex cut
3185 // cuts[3] - fP3 cut
3195 Double_t x[5], c[15];
3197 // make temporary seed
3198 AliTPCseed * seed = new AliTPCseed;
3199 Double_t alpha=fOuterSec->GetAlpha(), shift=fOuterSec->GetAlphaShift();
3200 // Double_t cs=cos(alpha), sn=sin(alpha);
3205 Double_t x1 = GetXrow(i1-1);
3206 const AliTPCtrackerRow& kr1=GetRow(sec,i1-1);
3207 Double_t y1max = GetMaxY(i1-1)-kr1.GetDeadZone()-1.5;
3209 Double_t x1p = GetXrow(i1);
3210 const AliTPCtrackerRow& kr1p=GetRow(sec,i1);
3212 Double_t x1m = GetXrow(i1-2);
3213 const AliTPCtrackerRow& kr1m=GetRow(sec,i1-2);
3216 //last 3 padrow for seeding
3217 AliTPCtrackerRow& kr3 = GetRow((sec+fkNOS)%fkNOS,i1-7);
3218 Double_t x3 = GetXrow(i1-7);
3219 // Double_t y3max= GetMaxY(i1-7)-kr3.fDeadZone-1.5;
3221 AliTPCtrackerRow& kr3p = GetRow((sec+fkNOS)%fkNOS,i1-6);
3222 Double_t x3p = GetXrow(i1-6);
3224 AliTPCtrackerRow& kr3m = GetRow((sec+fkNOS)%fkNOS,i1-8);
3225 Double_t x3m = GetXrow(i1-8);
3230 Int_t im = i1-4; //middle pad row index
3231 Double_t xm = GetXrow(im); // radius of middle pad-row
3232 const AliTPCtrackerRow& krm=GetRow(sec,im); //middle pad -row
3233 // Double_t ymmax = GetMaxY(im)-kr1.fDeadZone-1.5;
3236 Double_t deltax = x1-x3;
3237 Double_t dymax = deltax*cuts[1];
3238 Double_t dzmax = deltax*cuts[3];
3240 // loop over clusters
3241 for (Int_t is=0; is < kr1; is++) {
3243 if (kr1[is]->IsUsed(10)) continue;
3244 Double_t y1=kr1[is]->GetY(), z1=kr1[is]->GetZ();
3246 if (deltay>0 && TMath::Abs(y1max-TMath::Abs(y1))> deltay ) continue; // seed only at the edge
3248 Int_t index1 = TMath::Max(kr3.Find(z1-dzmax)-1,0);
3249 Int_t index2 = TMath::Min(kr3.Find(z1+dzmax)+1,kr3);
3255 for (Int_t js=index1; js < index2; js++) {
3256 const AliTPCclusterMI *kcl = kr3[js];
3257 if (kcl->IsUsed(10)) continue;
3259 // apply angular cuts
3260 if (TMath::Abs(y1-y3)>dymax) continue;
3263 if (TMath::Abs(z1-z3)>dzmax) continue;
3265 Double_t angley = (y1-y3)/(x1-x3);
3266 Double_t anglez = (z1-z3)/(x1-x3);
3268 Double_t erry = TMath::Abs(angley)*(x1-x1m)*0.5+0.5;
3269 Double_t errz = TMath::Abs(anglez)*(x1-x1m)*0.5+0.5;
3271 Double_t yyym = angley*(xm-x1)+y1;
3272 Double_t zzzm = anglez*(xm-x1)+z1;
3274 const AliTPCclusterMI *kcm = krm.FindNearest2(yyym,zzzm,erry,errz,index);
3276 if (kcm->IsUsed(10)) continue;
3278 erry = TMath::Abs(angley)*(x1-x1m)*0.4+0.5;
3279 errz = TMath::Abs(anglez)*(x1-x1m)*0.4+0.5;
3286 // look around first
3287 const AliTPCclusterMI *kc1m = kr1m.FindNearest2(angley*(x1m-x1)+y1,
3293 if (kc1m->IsUsed(10)) used++;
3295 const AliTPCclusterMI *kc1p = kr1p.FindNearest2(angley*(x1p-x1)+y1,
3301 if (kc1p->IsUsed(10)) used++;
3303 if (used>1) continue;
3304 if (found<1) continue;
3308 const AliTPCclusterMI *kc3m = kr3m.FindNearest2(angley*(x3m-x3)+y3,
3314 if (kc3m->IsUsed(10)) used++;
3318 const AliTPCclusterMI *kc3p = kr3p.FindNearest2(angley*(x3p-x3)+y3,
3324 if (kc3p->IsUsed(10)) used++;
3328 if (used>1) continue;
3329 if (found<3) continue;
3339 x[4]=F1(x1,y1,x2,y2,x3,y3);
3340 //if (TMath::Abs(x[4]) >= cuts[0]) continue;
3343 x[2]=F2(x1,y1,x2,y2,x3,y3);
3346 x[3]=F3n(x1,y1,x2,y2,z1,z2,x[4]);
3347 //if (TMath::Abs(x[3]) > cuts[3]) continue;
3351 Double_t sy1=0.1, sz1=0.1;
3352 Double_t sy2=0.1, sz2=0.1;
3353 Double_t sy3=0.1, sy=0.1, sz=0.1;
3355 Double_t f40=(F1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
3356 Double_t f42=(F1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
3357 Double_t f43=(F1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
3358 Double_t f20=(F2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
3359 Double_t f22=(F2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
3360 Double_t f23=(F2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
3362 Double_t f30=(F3(x1,y1+sy,x2,y2,z1,z2)-x[3])/sy;
3363 Double_t f31=(F3(x1,y1,x2,y2,z1+sz,z2)-x[3])/sz;
3364 Double_t f32=(F3(x1,y1,x2,y2+sy,z1,z2)-x[3])/sy;
3365 Double_t f34=(F3(x1,y1,x2,y2,z1,z2+sz)-x[3])/sz;
3369 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3370 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3371 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3372 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3373 c[13]=f30*sy1*f40+f32*sy2*f42;
3374 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3376 // if (!BuildSeed(kr1[is],kcl,0,x1,x2,x3,x,c)) continue;
3378 UInt_t index=kr1.GetIndex(is);
3379 seed->~AliTPCseed();
3380 AliTPCseed *track=new(seed) AliTPCseed(x1, sec*alpha+shift, x, c, index);
3382 track->SetIsSeeding(kTRUE);
3385 FollowProlongation(*track, i1-7,1);
3386 if (track->GetNumberOfClusters() < track->GetNFoundable()*0.75 ||
3387 track->GetNShared()>0.6*track->GetNumberOfClusters() || ( track->GetSigmaY2()+ track->GetSigmaZ2())>0.6){
3389 seed->~AliTPCseed();
3395 FollowProlongation(*track, i2,1);
3396 track->SetBConstrain(0);
3397 track->SetLastPoint(i1+fInnerSec->GetNRows()); // first cluster in track position
3398 track->SetFirstPoint(track->GetLastPoint());
3400 if (track->GetNumberOfClusters()<(i1-i2)*0.5 ||
3401 track->GetNumberOfClusters()<track->GetNFoundable()*0.7 ||
3402 track->GetNShared()>2. || track->GetChi2()/track->GetNumberOfClusters()>6 || ( track->GetSigmaY2()+ track->GetSigmaZ2())>0.5 ) {
3404 seed->~AliTPCseed();
3409 FollowProlongation(*track, TMath::Max(i2-10,0),1);
3410 AliTPCseed * track2 = MakeSeed(track,0.2,0.5,0.9);
3411 FollowProlongation(*track2, i2,1);
3412 track2->SetBConstrain(kFALSE);
3413 track2->SetSeedType(4);
3414 arr->AddLast(track2);
3416 seed->~AliTPCseed();
3420 //arr->AddLast(track);
3421 //seed = new AliTPCseed;
3427 Info("MakeSeeds5","\nSeeding statiistic:\t%d\t%d\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin,nout1,nout2,nout3);
3433 //_____________________________________________________________________________
3434 void AliTPCtrackerMI::MakeSeeds2(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, Float_t */*cuts[4]*/,
3435 Float_t deltay, Bool_t /*bconstrain*/) {
3436 //-----------------------------------------------------------------
3437 // This function creates track seeds - without vertex constraint
3438 //-----------------------------------------------------------------
3439 // cuts[0] - fP4 cut - not applied
3440 // cuts[1] - tan(phi) cut
3441 // cuts[2] - zvertex cut - not applied
3442 // cuts[3] - fP3 cut
3452 Double_t alpha=fOuterSec->GetAlpha(), shift=fOuterSec->GetAlphaShift();
3453 // Double_t cs=cos(alpha), sn=sin(alpha);
3454 Int_t row0 = (i1+i2)/2;
3455 Int_t drow = (i1-i2)/2;
3456 const AliTPCtrackerRow& kr0=fSectors[sec][row0];
3457 AliTPCtrackerRow * kr=0;
3459 AliTPCpolyTrack polytrack;
3460 Int_t nclusters=fSectors[sec][row0];
3461 AliTPCseed * seed = new AliTPCseed;
3466 for (Int_t is=0; is < nclusters; is++) { //LOOP over clusters
3468 Int_t nfoundable =0;
3469 for (Int_t iter =1; iter<2; iter++){ //iterations
3470 const AliTPCtrackerRow& krm=fSectors[sec][row0-iter];
3471 const AliTPCtrackerRow& krp=fSectors[sec][row0+iter];
3472 const AliTPCclusterMI * cl= kr0[is];
3474 if (cl->IsUsed(10)) {
3480 Double_t x = kr0.GetX();
3481 // Initialization of the polytrack
3486 Double_t y0= cl->GetY();
3487 Double_t z0= cl->GetZ();
3491 Double_t ymax = fSectors->GetMaxY(row0)-kr0.GetDeadZone()-1.5;
3492 if (deltay>0 && TMath::Abs(ymax-TMath::Abs(y0))> deltay ) continue; // seed only at the edge
3494 erry = (0.5)*cl->GetSigmaY2()/TMath::Sqrt(cl->GetQ())*6;
3495 errz = (0.5)*cl->GetSigmaZ2()/TMath::Sqrt(cl->GetQ())*6;
3496 polytrack.AddPoint(x,y0,z0,erry, errz);
3499 if (cl->IsUsed(10)) sumused++;
3502 Float_t roady = (5*TMath::Sqrt(cl->GetSigmaY2()+0.2)+1.)*iter;
3503 Float_t roadz = (5*TMath::Sqrt(cl->GetSigmaZ2()+0.2)+1.)*iter;
3506 AliTPCclusterMI * cl1 = krm.FindNearest(y0,z0,roady,roadz);
3507 if (cl1 && TMath::Abs(ymax-TMath::Abs(y0))) {
3508 erry = (0.5)*cl1->GetSigmaY2()/TMath::Sqrt(cl1->GetQ())*3;
3509 errz = (0.5)*cl1->GetSigmaZ2()/TMath::Sqrt(cl1->GetQ())*3;
3510 if (cl1->IsUsed(10)) sumused++;
3511 polytrack.AddPoint(x,cl1->GetY(),cl1->GetZ(),erry,errz);
3515 AliTPCclusterMI * cl2 = krp.FindNearest(y0,z0,roady,roadz);
3517 erry = (0.5)*cl2->GetSigmaY2()/TMath::Sqrt(cl2->GetQ())*3;
3518 errz = (0.5)*cl2->GetSigmaZ2()/TMath::Sqrt(cl2->GetQ())*3;
3519 if (cl2->IsUsed(10)) sumused++;
3520 polytrack.AddPoint(x,cl2->GetY(),cl2->GetZ(),erry,errz);
3523 if (sumused>0) continue;
3525 polytrack.UpdateParameters();
3531 nfoundable = polytrack.GetN();
3532 nfound = nfoundable;
3534 for (Int_t ddrow = iter+1; ddrow<drow;ddrow++){
3535 Float_t maxdist = 0.8*(1.+3./(ddrow));
3536 for (Int_t delta = -1;delta<=1;delta+=2){
3537 Int_t row = row0+ddrow*delta;
3538 kr = &(fSectors[sec][row]);
3539 Double_t xn = kr->GetX();
3540 Double_t ymax = fSectors->GetMaxY(row)-kr->GetDeadZone()-1.5;
3541 polytrack.GetFitPoint(xn,yn,zn);
3542 if (TMath::Abs(yn)>ymax) continue;
3544 AliTPCclusterMI * cln = kr->FindNearest(yn,zn,roady,roadz);
3546 Float_t dist = TMath::Sqrt( (yn-cln->GetY())*(yn-cln->GetY())+(zn-cln->GetZ())*(zn-cln->GetZ()));
3549 erry = (dist+0.3)*cln->GetSigmaY2()/TMath::Sqrt(cln->GetQ())*(1.+1./(ddrow));
3550 errz = (dist+0.3)*cln->GetSigmaZ2()/TMath::Sqrt(cln->GetQ())*(1.+1./(ddrow));
3551 if (cln->IsUsed(10)) {
3552 // printf("used\n");
3560 polytrack.AddPoint(xn,cln->GetY(),cln->GetZ(),erry, errz);
3565 if ( (sumused>3) || (sumused>0.5*nfound) || (nfound<0.6*nfoundable)) break;
3566 polytrack.UpdateParameters();
3569 if ( (sumused>3) || (sumused>0.5*nfound)) {
3570 //printf("sumused %d\n",sumused);
3575 polytrack.GetFitDerivation(kr0.GetX(),dy,dz);
3576 AliTPCpolyTrack track2;
3578 polytrack.Refit(track2,0.5+TMath::Abs(dy)*0.3,0.4+TMath::Abs(dz)*0.3);
3579 if (track2.GetN()<0.5*nfoundable) continue;
3582 if ((nfound>0.6*nfoundable) &&( nfoundable>0.4*(i1-i2))) {
3584 // test seed with and without constrain
3585 for (Int_t constrain=0; constrain<=0;constrain++){
3586 // add polytrack candidate
3588 Double_t x[5], c[15];
3589 Double_t x1,x2,x3,y1,y2,y3,z1,z2,z3;
3590 track2.GetBoundaries(x3,x1);
3592 track2.GetFitPoint(x1,y1,z1);
3593 track2.GetFitPoint(x2,y2,z2);
3594 track2.GetFitPoint(x3,y3,z3);
3596 //is track pointing to the vertex ?
3599 polytrack.GetFitPoint(x0,y0,z0);
3612 x[4]=F1(x1,y1,x2,y2,x3,y3);
3614 // if (TMath::Abs(x[4]) >= cuts[0]) continue; //
3615 x[2]=F2(x1,y1,x2,y2,x3,y3);
3617 //if (TMath::Abs(x[4]*x1-x[2]) >= cuts[1]) continue;
3618 //x[3]=F3(x1,y1,x2,y2,z1,z2);
3619 x[3]=F3n(x1,y1,x3,y3,z1,z3,x[4]);
3620 //if (TMath::Abs(x[3]) > cuts[3]) continue;
3623 Double_t sy =0.1, sz =0.1;
3624 Double_t sy1=0.02, sz1=0.02;
3625 Double_t sy2=0.02, sz2=0.02;
3629 sy3=25000*x[4]*x[4]+0.1, sy=0.1, sz=0.1;
3632 Double_t f40=(F1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
3633 Double_t f42=(F1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
3634 Double_t f43=(F1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
3635 Double_t f20=(F2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
3636 Double_t f22=(F2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
3637 Double_t f23=(F2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
3639 Double_t f30=(F3(x1,y1+sy,x3,y3,z1,z3)-x[3])/sy;
3640 Double_t f31=(F3(x1,y1,x3,y3,z1+sz,z3)-x[3])/sz;
3641 Double_t f32=(F3(x1,y1,x3,y3+sy,z1,z3)-x[3])/sy;
3642 Double_t f34=(F3(x1,y1,x3,y3,z1,z3+sz)-x[3])/sz;
3647 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3648 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3649 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3650 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3651 c[13]=f30*sy1*f40+f32*sy2*f42;
3652 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3654 //Int_t row1 = fSectors->GetRowNumber(x1);
3655 Int_t row1 = GetRowNumber(x1);
3659 seed->~AliTPCseed();
3660 AliTPCseed *track=new(seed) AliTPCseed(x1,sec*alpha+shift,x,c,index);
3661 track->SetIsSeeding(kTRUE);
3662 Int_t rc=FollowProlongation(*track, i2);
3663 if (constrain) track->SetBConstrain(1);
3665 track->SetBConstrain(0);
3666 track->SetLastPoint(row1+fInnerSec->GetNRows()); // first cluster in track position
3667 track->SetFirstPoint(track->GetLastPoint());
3669 if (rc==0 || track->GetNumberOfClusters()<(i1-i2)*0.5 ||
3670 track->GetNumberOfClusters() < track->GetNFoundable()*0.6 ||
3671 track->GetNShared()>0.4*track->GetNumberOfClusters()) {
3674 seed->~AliTPCseed();
3677 arr->AddLast(track);
3678 seed = new AliTPCseed;
3682 } // if accepted seed
3685 Info("MakeSeeds2","\nSeeding statiistic:\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin3);
3691 AliTPCseed *AliTPCtrackerMI::MakeSeed(AliTPCseed *track, Float_t r0, Float_t r1, Float_t r2)
3695 //reseed using track points
3696 Int_t p0 = int(r0*track->GetNumberOfClusters()); // point 0
3697 Int_t p1 = int(r1*track->GetNumberOfClusters());
3698 Int_t p2 = int(r2*track->GetNumberOfClusters()); // last point
3700 Double_t x0[3],x1[3],x2[3];
3701 for (Int_t i=0;i<3;i++){
3707 // find track position at given ratio of the length
3708 Int_t sec0=0, sec1=0, sec2=0;
3711 for (Int_t i=0;i<160;i++){
3712 if (track->GetClusterPointer(i)){
3714 AliTPCTrackerPoint *trpoint =track->GetTrackPoint(i);
3715 if ( (index<p0) || x0[0]<0 ){
3716 if (trpoint->GetX()>1){
3717 clindex = track->GetClusterIndex2(i);
3719 x0[0] = trpoint->GetX();
3720 x0[1] = trpoint->GetY();
3721 x0[2] = trpoint->GetZ();
3722 sec0 = ((clindex&0xff000000)>>24)%18;
3727 if ( (index<p1) &&(trpoint->GetX()>1)){
3728 clindex = track->GetClusterIndex2(i);
3730 x1[0] = trpoint->GetX();
3731 x1[1] = trpoint->GetY();
3732 x1[2] = trpoint->GetZ();
3733 sec1 = ((clindex&0xff000000)>>24)%18;
3736 if ( (index<p2) &&(trpoint->GetX()>1)){
3737 clindex = track->GetClusterIndex2(i);
3739 x2[0] = trpoint->GetX();
3740 x2[1] = trpoint->GetY();
3741 x2[2] = trpoint->GetZ();
3742 sec2 = ((clindex&0xff000000)>>24)%18;
3749 Double_t alpha, cs,sn, xx2,yy2;
3751 alpha = (sec1-sec2)*fSectors->GetAlpha();
3752 cs = TMath::Cos(alpha);
3753 sn = TMath::Sin(alpha);
3754 xx2= x1[0]*cs-x1[1]*sn;
3755 yy2= x1[0]*sn+x1[1]*cs;
3759 alpha = (sec0-sec2)*fSectors->GetAlpha();
3760 cs = TMath::Cos(alpha);
3761 sn = TMath::Sin(alpha);
3762 xx2= x0[0]*cs-x0[1]*sn;
3763 yy2= x0[0]*sn+x0[1]*cs;
3769 Double_t x[5],c[15];
3773 x[4]=F1(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]);
3774 // if (x[4]>1) return 0;
3775 x[2]=F2(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]);
3776 x[3]=F3n(x2[0],x2[1],x0[0],x0[1],x2[2],x0[2],x[4]);
3777 //if (TMath::Abs(x[3]) > 2.2) return 0;
3778 //if (TMath::Abs(x[2]) > 1.99) return 0;
3780 Double_t sy =0.1, sz =0.1;
3782 Double_t sy1=0.02+track->GetSigmaY2(), sz1=0.02+track->GetSigmaZ2();
3783 Double_t sy2=0.01+track->GetSigmaY2(), sz2=0.01+track->GetSigmaZ2();
3784 Double_t sy3=0.01+track->GetSigmaY2();
3786 Double_t f40=(F1(x2[0],x2[1]+sy,x1[0],x1[1],x0[0],x0[1])-x[4])/sy;
3787 Double_t f42=(F1(x2[0],x2[1],x1[0],x1[1]+sy,x0[0],x0[1])-x[4])/sy;
3788 Double_t f43=(F1(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]+sy)-x[4])/sy;
3789 Double_t f20=(F2(x2[0],x2[1]+sy,x1[0],x1[1],x0[0],x0[1])-x[2])/sy;
3790 Double_t f22=(F2(x2[0],x2[1],x1[0],x1[1]+sy,x0[0],x0[1])-x[2])/sy;
3791 Double_t f23=(F2(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]+sy)-x[2])/sy;
3793 Double_t f30=(F3(x2[0],x2[1]+sy,x0[0],x0[1],x2[2],x0[2])-x[3])/sy;
3794 Double_t f31=(F3(x2[0],x2[1],x0[0],x0[1],x2[2]+sz,x0[2])-x[3])/sz;
3795 Double_t f32=(F3(x2[0],x2[1],x0[0],x0[1]+sy,x2[2],x0[2])-x[3])/sy;
3796 Double_t f34=(F3(x2[0],x2[1],x0[0],x0[1],x2[2],x0[2]+sz)-x[3])/sz;
3801 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3802 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3803 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3804 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3805 c[13]=f30*sy1*f40+f32*sy2*f42;
3806 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3808 // Int_t row1 = fSectors->GetRowNumber(x2[0]);
3809 AliTPCseed *seed=new AliTPCseed(x2[0], sec2*fSectors->GetAlpha()+fSectors->GetAlphaShift(), x, c, 0);
3810 // Double_t y0,z0,y1,z1, y2,z2;
3811 //seed->GetProlongation(x0[0],y0,z0);
3812 // seed->GetProlongation(x1[0],y1,z1);
3813 //seed->GetProlongation(x2[0],y2,z2);
3815 seed->SetLastPoint(pp2);
3816 seed->SetFirstPoint(pp2);
3823 AliTPCseed *AliTPCtrackerMI::ReSeed(AliTPCseed *track, Float_t r0, Float_t r1, Float_t r2)
3827 //reseed using founded clusters
3829 // Find the number of clusters
3830 Int_t nclusters = 0;
3831 for (Int_t irow=0;irow<160;irow++){
3832 if (track->GetClusterIndex(irow)>0) nclusters++;
3836 ipos[0] = TMath::Max(int(r0*nclusters),0); // point 0 cluster
3837 ipos[1] = TMath::Min(int(r1*nclusters),nclusters-1); //
3838 ipos[2] = TMath::Min(int(r2*nclusters),nclusters-1); // last point
3842 Int_t row[3],sec[3]={0,0,0};
3844 // find track row position at given ratio of the length
3846 for (Int_t irow=0;irow<160;irow++){
3847 if (track->GetClusterIndex2(irow)<0) continue;
3849 for (Int_t ipoint=0;ipoint<3;ipoint++){
3850 if (index<=ipos[ipoint]) row[ipoint] = irow;
3854 //Get cluster and sector position
3855 for (Int_t ipoint=0;ipoint<3;ipoint++){
3856 Int_t clindex = track->GetClusterIndex2(row[ipoint]);
3857 AliTPCclusterMI * cl = GetClusterMI(clindex);
3860 // AliTPCclusterMI * cl = GetClusterMI(clindex);
3863 sec[ipoint] = ((clindex&0xff000000)>>24)%18;
3864 xyz[ipoint][0] = GetXrow(row[ipoint]);
3865 xyz[ipoint][1] = cl->GetY();
3866 xyz[ipoint][2] = cl->GetZ();
3870 // Calculate seed state vector and covariance matrix
3872 Double_t alpha, cs,sn, xx2,yy2;
3874 alpha = (sec[1]-sec[2])*fSectors->GetAlpha();
3875 cs = TMath::Cos(alpha);
3876 sn = TMath::Sin(alpha);
3877 xx2= xyz[1][0]*cs-xyz[1][1]*sn;
3878 yy2= xyz[1][0]*sn+xyz[1][1]*cs;
3882 alpha = (sec[0]-sec[2])*fSectors->GetAlpha();
3883 cs = TMath::Cos(alpha);
3884 sn = TMath::Sin(alpha);
3885 xx2= xyz[0][0]*cs-xyz[0][1]*sn;
3886 yy2= xyz[0][0]*sn+xyz[0][1]*cs;
3892 Double_t x[5],c[15];
3896 x[4]=F1(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
3897 x[2]=F2(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
3898 x[3]=F3n(xyz[2][0],xyz[2][1],xyz[0][0],xyz[0][1],xyz[2][2],xyz[0][2],x[4]);
3900 Double_t sy =0.1, sz =0.1;
3902 Double_t sy1=0.2, sz1=0.2;
3903 Double_t sy2=0.2, sz2=0.2;
3906 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;
3907 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;
3908 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;
3909 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;
3910 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;
3911 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;
3913 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;
3914 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;
3915 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;
3916 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;
3921 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3922 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3923 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3924 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3925 c[13]=f30*sy1*f40+f32*sy2*f42;
3926 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3928 // Int_t row1 = fSectors->GetRowNumber(xyz[2][0]);
3929 AliTPCseed *seed=new AliTPCseed(xyz[2][0], sec[2]*fSectors->GetAlpha()+fSectors->GetAlphaShift(), x, c, 0);
3930 seed->SetLastPoint(row[2]);
3931 seed->SetFirstPoint(row[2]);
3936 AliTPCseed *AliTPCtrackerMI::ReSeed(AliTPCseed *track,Int_t r0, Bool_t forward)
3940 //reseed using founded clusters
3943 Int_t row[3]={0,0,0};
3944 Int_t sec[3]={0,0,0};
3946 // forward direction
3948 for (Int_t irow=r0;irow<160;irow++){
3949 if (track->GetClusterIndex(irow)>0){
3954 for (Int_t irow=160;irow>r0;irow--){
3955 if (track->GetClusterIndex(irow)>0){
3960 for (Int_t irow=row[2]-15;irow>row[0];irow--){
3961 if (track->GetClusterIndex(irow)>0){
3969 for (Int_t irow=0;irow<r0;irow++){
3970 if (track->GetClusterIndex(irow)>0){
3975 for (Int_t irow=r0;irow>0;irow--){
3976 if (track->GetClusterIndex(irow)>0){
3981 for (Int_t irow=row[2]-15;irow>row[0];irow--){
3982 if (track->GetClusterIndex(irow)>0){
3989 if ((row[2]-row[0])<20) return 0;
3990 if (row[1]==0) return 0;
3993 //Get cluster and sector position
3994 for (Int_t ipoint=0;ipoint<3;ipoint++){
3995 Int_t clindex = track->GetClusterIndex2(row[ipoint]);
3996 AliTPCclusterMI * cl = GetClusterMI(clindex);
3999 // AliTPCclusterMI * cl = GetClusterMI(clindex);
4002 sec[ipoint] = ((clindex&0xff000000)>>24)%18;
4003 xyz[ipoint][0] = GetXrow(row[ipoint]);
4004 AliTPCTrackerPoint * point = track->GetTrackPoint(row[ipoint]);
4005 if (point&&ipoint<2){
4007 xyz[ipoint][1] = point->GetY();
4008 xyz[ipoint][2] = point->GetZ();
4011 xyz[ipoint][1] = cl->GetY();
4012 xyz[ipoint][2] = cl->GetZ();
4019 // Calculate seed state vector and covariance matrix
4021 Double_t alpha, cs,sn, xx2,yy2;
4023 alpha = (sec[1]-sec[2])*fSectors->GetAlpha();
4024 cs = TMath::Cos(alpha);
4025 sn = TMath::Sin(alpha);
4026 xx2= xyz[1][0]*cs-xyz[1][1]*sn;
4027 yy2= xyz[1][0]*sn+xyz[1][1]*cs;
4031 alpha = (sec[0]-sec[2])*fSectors->GetAlpha();
4032 cs = TMath::Cos(alpha);
4033 sn = TMath::Sin(alpha);
4034 xx2= xyz[0][0]*cs-xyz[0][1]*sn;
4035 yy2= xyz[0][0]*sn+xyz[0][1]*cs;
4041 Double_t x[5],c[15];
4045 x[4]=F1(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
4046 x[2]=F2(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
4047 x[3]=F3n(xyz[2][0],xyz[2][1],xyz[0][0],xyz[0][1],xyz[2][2],xyz[0][2],x[4]);
4049 Double_t sy =0.1, sz =0.1;
4051 Double_t sy1=0.2, sz1=0.2;
4052 Double_t sy2=0.2, sz2=0.2;
4055 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;
4056 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;
4057 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;
4058 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;
4059 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;
4060 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;
4062 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;
4063 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;
4064 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;
4065 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;
4070 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
4071 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
4072 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
4073 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
4074 c[13]=f30*sy1*f40+f32*sy2*f42;
4075 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
4077 // Int_t row1 = fSectors->GetRowNumber(xyz[2][0]);
4078 AliTPCseed *seed=new AliTPCseed(xyz[2][0], sec[2]*fSectors->GetAlpha()+fSectors->GetAlphaShift(), x, c, 0);
4079 seed->SetLastPoint(row[2]);
4080 seed->SetFirstPoint(row[2]);
4081 for (Int_t i=row[0];i<row[2];i++){
4082 seed->SetClusterIndex(i, track->GetClusterIndex(i));
4090 void AliTPCtrackerMI::FindMultiMC(TObjArray * array, AliESDEvent */*esd*/, Int_t iter)
4093 // find multi tracks - THIS FUNCTION IS ONLY FOR DEBUG PURPOSES
4095 // Use AliTPCReconstructor::StreamLevel()>2 if you want to tune parameters - cuts
4097 // Two reasons to have multiple find tracks
4098 // 1. Curling tracks can be find more than once
4099 // 2. Splitted tracks
4100 // a.) Multiple seeding to increase tracking efficiency - (~ 100% reached)
4101 // b.) Edge effect on the sector boundaries
4104 // Algorithm done in 2 phases - because of CPU consumption
4105 // it is n^2 algorithm - for lead-lead 20000x20000 combination are investigated
4107 // Algorihm for curling tracks sign:
4108 // 1 phase -makes a very rough fast cuts to minimize combinatorics
4109 // a.) opposite sign
4110 // b.) one of the tracks - not pointing to the primary vertex -
4111 // c.) delta tan(theta)
4113 // 2 phase - calculates DCA between tracks - time consument
4118 // General cuts - for splitted tracks and for curling tracks
4120 const Float_t kMaxdPhi = 0.2; // maximal distance in phi
4122 // Curling tracks cuts
4127 Int_t nentries = array->GetEntriesFast();
4128 AliHelix *helixes = new AliHelix[nentries];
4129 Float_t *xm = new Float_t[nentries];
4130 Float_t *dz0 = new Float_t[nentries];
4131 Float_t *dz1 = new Float_t[nentries];
4137 // Find track COG in x direction - point with best defined parameters
4139 for (Int_t i=0;i<nentries;i++){
4140 AliTPCseed* track = (AliTPCseed*)array->At(i);
4141 if (!track) continue;
4142 track->SetCircular(0);
4143 new (&helixes[i]) AliHelix(*track);
4147 track->GetDZ(GetX(),GetY(),GetZ(),GetBz(),dz);
4150 for (Int_t icl=0; icl<160; icl++){
4151 AliTPCclusterMI * cl = track->GetClusterPointer(icl);
4157 if (ncl>0) xm[i]/=Float_t(ncl);
4159 TTreeSRedirector &cstream = *fDebugStreamer;
4161 for (Int_t i0=0;i0<nentries;i0++){
4162 AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
4163 if (!track0) continue;
4164 Float_t xc0 = helixes[i0].GetHelix(6);
4165 Float_t yc0 = helixes[i0].GetHelix(7);
4166 Float_t r0 = helixes[i0].GetHelix(8);
4167 Float_t rc0 = TMath::Sqrt(xc0*xc0+yc0*yc0);
4168 Float_t fi0 = TMath::ATan2(yc0,xc0);
4170 for (Int_t i1=i0+1;i1<nentries;i1++){
4171 AliTPCseed * track1 = (AliTPCseed*)array->At(i1);
4172 if (!track1) continue;
4173 Int_t lab0=track0->GetLabel();
4174 Int_t lab1=track1->GetLabel();
4175 if (TMath::Abs(lab0)!=TMath::Abs(lab1)) continue;
4177 Float_t xc1 = helixes[i1].GetHelix(6);
4178 Float_t yc1 = helixes[i1].GetHelix(7);
4179 Float_t r1 = helixes[i1].GetHelix(8);
4180 Float_t rc1 = TMath::Sqrt(xc1*xc1+yc1*yc1);
4181 Float_t fi1 = TMath::ATan2(yc1,xc1);
4183 Float_t dfi = fi0-fi1;
4186 if (dfi>1.5*TMath::Pi()) dfi-=TMath::Pi(); // take care about edge effect
4187 if (dfi<-1.5*TMath::Pi()) dfi+=TMath::Pi(); //
4188 if (TMath::Abs(dfi)>kMaxdPhi&&helixes[i0].GetHelix(4)*helixes[i1].GetHelix(4)<0){
4190 // if short tracks with undefined sign
4191 fi1 = -TMath::ATan2(yc1,-xc1);
4194 Float_t dtheta = TMath::Abs(track0->GetTgl()-track1->GetTgl())<TMath::Abs(track0->GetTgl()+track1->GetTgl())? track0->GetTgl()-track1->GetTgl():track0->GetTgl()+track1->GetTgl();
4197 // debug stream to tune "fast cuts"
4199 Double_t dist[3]; // distance at X
4200 Double_t mdist[3]={0,0,0}; // mean distance X+-40cm
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])+40.,dist,AliTracker::GetBz());
4204 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[i]);
4205 track0->GetDistance(track1,0.5*(xm[i0]+xm[i1]),dist,AliTracker::GetBz());
4206 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[i]);
4207 for (Int_t i=0;i<3;i++) mdist[i]*=0.33333;
4211 for (Int_t icl=0; icl<160; icl++){
4212 AliTPCclusterMI * cl0 = track0->GetClusterPointer(icl);
4213 AliTPCclusterMI * cl1 = track1->GetClusterPointer(icl);
4216 if (cl0==cl1) sums++;
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
4259 if (AliTPCReconstructor::StreamLevel()>1) {
4260 AliInfo("Time for curling tracks removal DEBUGGING MC");
4266 void AliTPCtrackerMI::FindSplitted(TObjArray * array, AliESDEvent */*esd*/, Int_t iter)
4270 // Two reasons to have multiple find tracks
4271 // 1. Curling tracks can be find more than once
4272 // 2. Splitted tracks
4273 // a.) Multiple seeding to increase tracking efficiency - (~ 100% reached)
4274 // b.) Edge effect on the sector boundaries
4276 // This function tries to find tracks closed in the parametric space
4278 // cut logic if distance is bigger than cut continue - Do Nothing
4279 const Float_t kMaxdTheta = 0.05; // maximal distance in theta
4280 const Float_t kMaxdPhi = 0.05; // maximal deistance in phi
4281 const Float_t kdelta = 40.; //delta r to calculate track distance
4283 // const Float_t kMaxDist0 = 1.; // maximal distance 0
4284 //const Float_t kMaxDist1 = 0.3; // maximal distance 1 - cut if track in separate rows
4287 TCut csec("csec","abs(Tr0.fRelativeSector-Tr1.fRelativeSector)<2");
4288 TCut cdtheta("cdtheta","abs(dtheta)<0.05");
4293 Int_t nentries = array->GetEntriesFast();
4294 AliHelix *helixes = new AliHelix[nentries];
4295 Float_t *xm = new Float_t[nentries];
4301 //Sort tracks according quality
4303 Int_t nseed = array->GetEntriesFast();
4304 Float_t * quality = new Float_t[nseed];
4305 Int_t * indexes = new Int_t[nseed];
4306 for (Int_t i=0; i<nseed; i++) {
4307 AliTPCseed *pt=(AliTPCseed*)array->UncheckedAt(i);
4312 pt->UpdatePoints(); //select first last max dens points
4313 Float_t * points = pt->GetPoints();
4314 if (points[3]<0.8) quality[i] =-1;
4315 quality[i] = (points[2]-points[0])+pt->GetNumberOfClusters();
4316 //prefer high momenta tracks if overlaps
4317 quality[i] *= TMath::Sqrt(TMath::Abs(pt->Pt())+0.5);
4319 TMath::Sort(nseed,quality,indexes);
4323 // Find track COG in x direction - point with best defined parameters
4325 for (Int_t i=0;i<nentries;i++){
4326 AliTPCseed* track = (AliTPCseed*)array->At(i);
4327 if (!track) continue;
4328 track->SetCircular(0);
4329 new (&helixes[i]) AliHelix(*track);
4332 for (Int_t icl=0; icl<160; icl++){
4333 AliTPCclusterMI * cl = track->GetClusterPointer(icl);
4339 if (ncl>0) xm[i]/=Float_t(ncl);
4341 TTreeSRedirector &cstream = *fDebugStreamer;
4343 for (Int_t is0=0;is0<nentries;is0++){
4344 Int_t i0 = indexes[is0];
4345 AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
4346 if (!track0) continue;
4347 if (track0->GetKinkIndexes()[0]!=0) continue;
4348 Float_t xc0 = helixes[i0].GetHelix(6);
4349 Float_t yc0 = helixes[i0].GetHelix(7);
4350 Float_t fi0 = TMath::ATan2(yc0,xc0);
4352 for (Int_t is1=is0+1;is1<nentries;is1++){
4353 Int_t i1 = indexes[is1];
4354 AliTPCseed * track1 = (AliTPCseed*)array->At(i1);
4355 if (!track1) continue;
4357 if (TMath::Abs(track0->GetRelativeSector()-track1->GetRelativeSector())>1) continue;
4358 if (track1->GetKinkIndexes()[0]>0 &&track0->GetKinkIndexes()[0]<0) continue;
4359 if (track1->GetKinkIndexes()[0]!=0) continue;
4361 Float_t dtheta = TMath::Abs(track0->GetTgl()-track1->GetTgl())<TMath::Abs(track0->GetTgl()+track1->GetTgl())? track0->GetTgl()-track1->GetTgl():track0->GetTgl()+track1->GetTgl();
4362 if (TMath::Abs(dtheta)>kMaxdTheta) continue;
4364 Float_t xc1 = helixes[i1].GetHelix(6);
4365 Float_t yc1 = helixes[i1].GetHelix(7);
4366 Float_t fi1 = TMath::ATan2(yc1,xc1);
4368 Float_t dfi = fi0-fi1;
4369 if (dfi>1.5*TMath::Pi()) dfi-=TMath::Pi(); // take care about edge effect
4370 if (dfi<-1.5*TMath::Pi()) dfi+=TMath::Pi(); //
4371 if (TMath::Abs(dfi)>kMaxdPhi&&helixes[i0].GetHelix(4)*helixes[i1].GetHelix(4)<0){
4373 // if short tracks with undefined sign
4374 fi1 = -TMath::ATan2(yc1,-xc1);
4377 if (TMath::Abs(dfi)>kMaxdPhi) continue;
4384 for (Int_t icl=0; icl<160; icl++){
4385 Int_t index0=track0->GetClusterIndex2(icl);
4386 Int_t index1=track1->GetClusterIndex2(icl);
4387 Bool_t used0 = (index0>0 && !(index0&0x8000));
4388 Bool_t used1 = (index1>0 && !(index1&0x8000));
4390 if (used0) sum0++; // used cluster0
4391 if (used1) sum1++; // used clusters1
4392 if (used0&&used1) sum++;
4393 if (index0==index1 && used0 && used1) sums++;
4397 if (sums<10) continue;
4398 if (sum<40) continue;
4399 if (sums/Float_t(TMath::Min(sum0,sum1))<0.5) continue;
4401 Double_t dist[5][4]; // distance at X
4402 Double_t mdist[4]={0,0,0,0}; // mean distance on range +-delta
4406 track0->GetDistance(track1,xm[i0],dist[0],AliTracker::GetBz());
4407 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[0][i]);
4408 track0->GetDistance(track1,xm[i1],dist[1],AliTracker::GetBz());
4409 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[1][i]);
4411 track0->GetDistance(track1,TMath::Min(xm[i1],xm[i0])-kdelta,dist[2],AliTracker::GetBz());
4412 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[2][i]);
4413 track0->GetDistance(track1,TMath::Max(xm[i1],xm[i0])+kdelta,dist[3],AliTracker::GetBz());
4414 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[3][i]);
4416 track0->GetDistance(track1,(xm[i1]+xm[i0])*0.5,dist[4],AliTracker::GetBz());
4417 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[4][i]);
4418 for (Int_t i=0;i<3;i++) mdist[i]*=0.2;
4421 Int_t lab0=track0->GetLabel();
4422 Int_t lab1=track1->GetLabel();
4423 cstream<<"Splitted"<<
4427 "Tr0.="<<track0<< // seed0
4428 "Tr1.="<<track1<< // seed1
4429 "h0.="<<&helixes[i0]<<
4430 "h1.="<<&helixes[i1]<<
4432 "sum="<<sum<< //the sum of rows with cl in both
4433 "sum0="<<sum0<< //the sum of rows with cl in 0 track
4434 "sum1="<<sum1<< //the sum of rows with cl in 1 track
4435 "sums="<<sums<< //the sum of shared clusters
4436 "xm0="<<xm[i0]<< // the center of track
4437 "xm1="<<xm[i1]<< // the x center of track
4438 // General cut variables
4439 "dfi="<<dfi<< // distance in fi angle
4440 "dtheta="<<dtheta<< // distance int theta angle
4443 "dist0="<<dist[4][0]<< //distance x
4444 "dist1="<<dist[4][1]<< //distance y
4445 "dist2="<<dist[4][2]<< //distance z
4446 "mdist0="<<mdist[0]<< //distance x
4447 "mdist1="<<mdist[1]<< //distance y
4448 "mdist2="<<mdist[2]<< //distance z
4451 delete array->RemoveAt(i1);
4456 AliInfo("Time for splitted tracks removal");
4462 void AliTPCtrackerMI::FindCurling(TObjArray * array, AliESDEvent */*esd*/, Int_t iter)
4465 // find Curling tracks
4466 // Use AliTPCReconstructor::StreamLevel()>1 if you want to tune parameters - cuts
4469 // Algorithm done in 2 phases - because of CPU consumption
4470 // it is n^2 algorithm - for lead-lead 20000x20000 combination are investigated
4471 // see detal in MC part what can be used to cut
4475 const Float_t kMaxC = 400; // maximal curvature to of the track
4476 const Float_t kMaxdTheta = 0.15; // maximal distance in theta
4477 const Float_t kMaxdPhi = 0.15; // maximal distance in phi
4478 const Float_t kPtRatio = 0.3; // ratio between pt
4479 const Float_t kMinDCAR = 2.; // distance to the primary vertex in r - see cpipe cut
4482 // Curling tracks cuts
4485 const Float_t kMaxDeltaRMax = 40; // distance in outer radius
4486 const Float_t kMaxDeltaRMin = 5.; // distance in lower radius - see cpipe cut
4487 const Float_t kMinAngle = 2.9; // angle between tracks
4488 const Float_t kMaxDist = 5; // biggest distance
4490 // The cuts can be tuned using the "MC information stored in Multi tree ==> see FindMultiMC
4493 TCut csign("csign","Tr0.fP[4]*Tr1.fP[4]<0"); //opposite sign
4494 TCut cmax("cmax","abs(Tr0.GetC())>1/400");
4495 TCut cda("cda","sqrt(dtheta^2+dfi^2)<0.15");
4496 TCut ccratio("ccratio","abs((Tr0.fP[4]+Tr1.fP[4])/(abs(Tr0.fP[4])+abs(Tr1.fP[4])))<0.3");
4497 TCut cpipe("cpipe", "min(abs(r0-rc0),abs(r1-rc1))>5");
4499 TCut cdrmax("cdrmax","abs(abs(rc0+r0)-abs(rc1+r1))<40")
4500 TCut cdrmin("cdrmin","abs(abs(rc0+r0)-abs(rc1+r1))<10")
4502 Multi->Draw("dfi","iter==0"+csign+cmax+cda+ccratio); ~94% of curling tracks fulfill
4503 Multi->Draw("min(abs(r0-rc0),abs(r1-rc1))","iter==0&&abs(lab1)==abs(lab0)"+csign+cmax+cda+ccratio+cpipe+cdrmin+cdrmax); //80%
4505 Curling2->Draw("dfi","iter==0&&abs(lab0)==abs(lab1)"+csign+cmax+cdtheta+cdfi+ccratio)
4511 Int_t nentries = array->GetEntriesFast();
4512 AliHelix *helixes = new AliHelix[nentries];
4513 for (Int_t i=0;i<nentries;i++){
4514 AliTPCseed* track = (AliTPCseed*)array->At(i);
4515 if (!track) continue;
4516 track->SetCircular(0);
4517 new (&helixes[i]) AliHelix(*track);
4523 Double_t phase[2][2],radius[2];
4527 TTreeSRedirector &cstream = *fDebugStreamer;
4529 for (Int_t i0=0;i0<nentries;i0++){
4530 AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
4531 if (!track0) continue;
4532 if (TMath::Abs(track0->GetC())<1/kMaxC) continue;
4533 Float_t xc0 = helixes[i0].GetHelix(6);
4534 Float_t yc0 = helixes[i0].GetHelix(7);
4535 Float_t r0 = helixes[i0].GetHelix(8);
4536 Float_t rc0 = TMath::Sqrt(xc0*xc0+yc0*yc0);
4537 Float_t fi0 = TMath::ATan2(yc0,xc0);
4539 for (Int_t i1=i0+1;i1<nentries;i1++){
4540 AliTPCseed * track1 = (AliTPCseed*)array->At(i1);
4541 if (!track1) continue;
4542 if (TMath::Abs(track1->GetC())<1/kMaxC) continue;
4543 Float_t xc1 = helixes[i1].GetHelix(6);
4544 Float_t yc1 = helixes[i1].GetHelix(7);
4545 Float_t r1 = helixes[i1].GetHelix(8);
4546 Float_t rc1 = TMath::Sqrt(xc1*xc1+yc1*yc1);
4547 Float_t fi1 = TMath::ATan2(yc1,xc1);
4549 Float_t dfi = fi0-fi1;
4552 if (dfi>1.5*TMath::Pi()) dfi-=TMath::Pi(); // take care about edge effect
4553 if (dfi<-1.5*TMath::Pi()) dfi+=TMath::Pi(); //
4554 Float_t dtheta = TMath::Abs(track0->GetTgl()-track1->GetTgl())<TMath::Abs(track0->GetTgl()+track1->GetTgl())? track0->GetTgl()-track1->GetTgl():track0->GetTgl()+track1->GetTgl();
4558 if (track0->GetBConstrain()&&track1->GetBConstrain()) continue; // not constrained
4559 if (track1->GetSigned1Pt()*track0->GetSigned1Pt()>0) continue; // not the same sign
4560 if ( TMath::Abs(track1->GetTgl()+track0->GetTgl())>kMaxdTheta) continue; //distance in the Theta
4561 if ( TMath::Abs(dfi)>kMaxdPhi) continue; //distance in phi
4562 if ( TMath::Sqrt(dfi*dfi+dtheta*dtheta)>kMaxdPhi) continue; //common angular offset
4564 Float_t pt0 = track0->GetSignedPt();
4565 Float_t pt1 = track1->GetSignedPt();
4566 if ((TMath::Abs(pt0+pt1)/(TMath::Abs(pt0)+TMath::Abs(pt1)))>kPtRatio) continue;
4567 if ((iter==1) && TMath::Abs(TMath::Abs(rc0+r0)-TMath::Abs(rc1+r1))>kMaxDeltaRMax) continue;
4568 if ((iter!=1) &&TMath::Abs(TMath::Abs(rc0-r0)-TMath::Abs(rc1-r1))>kMaxDeltaRMin) continue;
4569 if (TMath::Min(TMath::Abs(rc0-r0),TMath::Abs(rc1-r1))<kMinDCAR) continue;
4572 // Now find closest approach
4576 Int_t npoints = helixes[i0].GetRPHIintersections(helixes[i1], phase, radius,10);
4577 if (npoints==0) continue;
4578 helixes[i0].GetClosestPhases(helixes[i1], phase);
4582 Double_t hangles[3];
4583 helixes[i0].Evaluate(phase[0][0],xyz0);
4584 helixes[i1].Evaluate(phase[0][1],xyz1);
4586 helixes[i0].GetAngle(phase[0][0],helixes[i1],phase[0][1],hangles);
4587 Double_t deltah[2],deltabest;
4588 if (TMath::Abs(hangles[2])<kMinAngle) continue;
4592 helixes[i0].ParabolicDCA(helixes[i1],phase[0][0],phase[0][1],radius[0],deltah[0],2);
4594 helixes[i0].ParabolicDCA(helixes[i1],phase[1][0],phase[1][1],radius[1],deltah[1],2);
4595 if (deltah[1]<deltah[0]) ibest=1;
4597 deltabest = TMath::Sqrt(deltah[ibest]);
4598 helixes[i0].Evaluate(phase[ibest][0],xyz0);
4599 helixes[i1].Evaluate(phase[ibest][1],xyz1);
4600 helixes[i0].GetAngle(phase[ibest][0],helixes[i1],phase[ibest][1],hangles);
4601 Double_t radiusbest = TMath::Sqrt(radius[ibest]);
4603 if (deltabest>kMaxDist) continue;
4604 // if (mindcar+mindcaz<40 && (TMath::Abs(hangles[2])<kMinAngle ||deltabest>3)) continue;
4605 Bool_t sign =kFALSE;
4606 if (hangles[2]>kMinAngle) sign =kTRUE;
4609 // circular[i0] = kTRUE;
4610 // circular[i1] = kTRUE;
4611 if (track0->OneOverPt()<track1->OneOverPt()){
4612 track0->SetCircular(track0->GetCircular()+1);
4613 track1->SetCircular(track1->GetCircular()+2);
4616 track1->SetCircular(track1->GetCircular()+1);
4617 track0->SetCircular(track0->GetCircular()+2);
4620 if (AliTPCReconstructor::StreamLevel()>1){
4622 //debug stream to tune "fine" cuts
4623 Int_t lab0=track0->GetLabel();
4624 Int_t lab1=track1->GetLabel();
4625 cstream<<"Curling2"<<
4641 "npoints="<<npoints<<
4642 "hangles0="<<hangles[0]<<
4643 "hangles1="<<hangles[1]<<
4644 "hangles2="<<hangles[2]<<
4647 "radius="<<radiusbest<<
4648 "deltabest="<<deltabest<<
4649 "phase0="<<phase[ibest][0]<<
4650 "phase1="<<phase[ibest][1]<<
4658 if (AliTPCReconstructor::StreamLevel()>1) {
4659 AliInfo("Time for curling tracks removal");
4668 void AliTPCtrackerMI::FindKinks(TObjArray * array, AliESDEvent *esd)
4675 TObjArray *kinks= new TObjArray(10000);
4676 // TObjArray *v0s= new TObjArray(10000);
4677 Int_t nentries = array->GetEntriesFast();
4678 AliHelix *helixes = new AliHelix[nentries];
4679 Int_t *sign = new Int_t[nentries];
4680 Int_t *nclusters = new Int_t[nentries];
4681 Float_t *alpha = new Float_t[nentries];
4682 AliKink *kink = new AliKink();
4683 Int_t * usage = new Int_t[nentries];
4684 Float_t *zm = new Float_t[nentries];
4685 Float_t *z0 = new Float_t[nentries];
4686 Float_t *fim = new Float_t[nentries];
4687 Float_t *shared = new Float_t[nentries];
4688 Bool_t *circular = new Bool_t[nentries];
4689 Float_t *dca = new Float_t[nentries];
4690 //const AliESDVertex * primvertex = esd->GetVertex();
4692 // nentries = array->GetEntriesFast();
4697 for (Int_t i=0;i<nentries;i++){
4700 AliTPCseed* track = (AliTPCseed*)array->At(i);
4701 if (!track) continue;
4702 track->SetCircular(0);
4704 track->UpdatePoints();
4705 if (( track->GetPoints()[2]- track->GetPoints()[0])>5 && track->GetPoints()[3]>0.8){
4707 nclusters[i]=track->GetNumberOfClusters();
4708 alpha[i] = track->GetAlpha();
4709 new (&helixes[i]) AliHelix(*track);
4711 helixes[i].Evaluate(0,xyz);
4712 sign[i] = (track->GetC()>0) ? -1:1;
4715 if (track->GetProlongation(x,y,z)){
4717 fim[i] = alpha[i]+TMath::ATan2(y,x);
4720 zm[i] = track->GetZ();
4724 circular[i]= kFALSE;
4725 if (track->GetProlongation(0,y,z)) z0[i] = z;
4726 dca[i] = track->GetD(0,0);
4732 Int_t ncandidates =0;
4735 Double_t phase[2][2],radius[2];
4738 // Find circling track
4739 TTreeSRedirector &cstream = *fDebugStreamer;
4741 for (Int_t i0=0;i0<nentries;i0++){
4742 AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
4743 if (!track0) continue;
4744 if (track0->GetNumberOfClusters()<40) continue;
4745 if (TMath::Abs(1./track0->GetC())>200) continue;
4746 for (Int_t i1=i0+1;i1<nentries;i1++){
4747 AliTPCseed * track1 = (AliTPCseed*)array->At(i1);
4748 if (!track1) continue;
4749 if (track1->GetNumberOfClusters()<40) continue;
4750 if ( TMath::Abs(track1->GetTgl()+track0->GetTgl())>0.1) continue;
4751 if (track0->GetBConstrain()&&track1->GetBConstrain()) continue;
4752 if (TMath::Abs(1./track1->GetC())>200) continue;
4753 if (track1->GetSigned1Pt()*track0->GetSigned1Pt()>0) continue;
4754 if (track1->GetTgl()*track0->GetTgl()>0) continue;
4755 if (TMath::Max(TMath::Abs(1./track0->GetC()),TMath::Abs(1./track1->GetC()))>190) continue;
4756 if (track0->GetBConstrain()&&track1->OneOverPt()<track0->OneOverPt()) continue; //returning - lower momenta
4757 if (track1->GetBConstrain()&&track0->OneOverPt()<track1->OneOverPt()) continue; //returning - lower momenta
4759 Float_t mindcar = TMath::Min(TMath::Abs(dca[i0]),TMath::Abs(dca[i1]));
4760 if (mindcar<5) continue;
4761 Float_t mindcaz = TMath::Min(TMath::Abs(z0[i0]-GetZ()),TMath::Abs(z0[i1]-GetZ()));
4762 if (mindcaz<5) continue;
4763 if (mindcar+mindcaz<20) continue;
4766 Float_t xc0 = helixes[i0].GetHelix(6);
4767 Float_t yc0 = helixes[i0].GetHelix(7);
4768 Float_t r0 = helixes[i0].GetHelix(8);
4769 Float_t xc1 = helixes[i1].GetHelix(6);
4770 Float_t yc1 = helixes[i1].GetHelix(7);
4771 Float_t r1 = helixes[i1].GetHelix(8);
4773 Float_t rmean = (r0+r1)*0.5;
4774 Float_t delta =TMath::Sqrt((xc1-xc0)*(xc1-xc0)+(yc1-yc0)*(yc1-yc0));
4775 //if (delta>30) continue;
4776 if (delta>rmean*0.25) continue;
4777 if (TMath::Abs(r0-r1)/rmean>0.3) continue;
4779 Int_t npoints = helixes[i0].GetRPHIintersections(helixes[i1], phase, radius,10);
4780 if (npoints==0) continue;
4781 helixes[i0].GetClosestPhases(helixes[i1], phase);
4785 Double_t hangles[3];
4786 helixes[i0].Evaluate(phase[0][0],xyz0);
4787 helixes[i1].Evaluate(phase[0][1],xyz1);
4789 helixes[i0].GetAngle(phase[0][0],helixes[i1],phase[0][1],hangles);
4790 Double_t deltah[2],deltabest;
4791 if (hangles[2]<2.8) continue;
4794 helixes[i0].ParabolicDCA(helixes[i1],phase[0][0],phase[0][1],radius[0],deltah[0],2);
4796 helixes[i0].ParabolicDCA(helixes[i1],phase[1][0],phase[1][1],radius[1],deltah[1],2);
4797 if (deltah[1]<deltah[0]) ibest=1;
4799 deltabest = TMath::Sqrt(deltah[ibest]);
4800 helixes[i0].Evaluate(phase[ibest][0],xyz0);
4801 helixes[i1].Evaluate(phase[ibest][1],xyz1);
4802 helixes[i0].GetAngle(phase[ibest][0],helixes[i1],phase[ibest][1],hangles);
4803 Double_t radiusbest = TMath::Sqrt(radius[ibest]);
4805 if (deltabest>6) continue;
4806 if (mindcar+mindcaz<40 && (hangles[2]<3.12||deltabest>3)) continue;
4807 Bool_t sign =kFALSE;
4808 if (hangles[2]>3.06) sign =kTRUE;
4811 circular[i0] = kTRUE;
4812 circular[i1] = kTRUE;
4813 if (track0->OneOverPt()<track1->OneOverPt()){
4814 track0->SetCircular(track0->GetCircular()+1);
4815 track1->SetCircular(track1->GetCircular()+2);
4818 track1->SetCircular(track1->GetCircular()+1);
4819 track0->SetCircular(track0->GetCircular()+2);
4822 if (sign&&AliTPCReconstructor::StreamLevel()>1){
4824 Int_t lab0=track0->GetLabel();
4825 Int_t lab1=track1->GetLabel();
4826 cstream<<"Curling"<<
4833 "mindcar="<<mindcar<<
4834 "mindcaz="<<mindcaz<<
4837 "npoints="<<npoints<<
4838 "hangles0="<<hangles[0]<<
4839 "hangles2="<<hangles[2]<<
4844 "radius="<<radiusbest<<
4845 "deltabest="<<deltabest<<
4846 "phase0="<<phase[ibest][0]<<
4847 "phase1="<<phase[ibest][1]<<
4857 for (Int_t i =0;i<nentries;i++){
4858 if (sign[i]==0) continue;
4859 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
4866 Double_t cradius0 = 40*40;
4867 Double_t cradius1 = 270*270;
4870 Double_t cdist3=0.55;
4871 for (Int_t j =i+1;j<nentries;j++){
4873 if (sign[j]*sign[i]<1) continue;
4874 if ( (nclusters[i]+nclusters[j])>200) continue;
4875 if ( (nclusters[i]+nclusters[j])<80) continue;
4876 if ( TMath::Abs(zm[i]-zm[j])>60.) continue;
4877 if ( TMath::Abs(fim[i]-fim[j])>0.6 && TMath::Abs(fim[i]-fim[j])<5.7 ) continue;
4878 //AliTPCseed * track1 = (AliTPCseed*)array->At(j); Double_t phase[2][2],radius[2];
4879 Int_t npoints = helixes[i].GetRPHIintersections(helixes[j], phase, radius,20);
4880 if (npoints<1) continue;
4883 if (radius[0]<cradius0||radius[0]>cradius1) continue;
4886 if ( (radius[0]<cradius0||radius[0]>cradius1) && (radius[1]<cradius0||radius[1]>cradius1) ) continue;
4889 Double_t delta1=10000,delta2=10000;
4890 // cuts on the intersection radius
4891 helixes[i].LinearDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta1);
4892 if (radius[0]<20&&delta1<1) continue; //intersection at vertex
4893 if (radius[0]<10&&delta1<3) continue; //intersection at vertex
4895 helixes[i].LinearDCA(helixes[j],phase[1][0],phase[1][1],radius[1],delta2);
4896 if (radius[1]<20&&delta2<1) continue; //intersection at vertex
4897 if (radius[1]<10&&delta2<3) continue; //intersection at vertex
4900 Double_t distance1 = TMath::Min(delta1,delta2);
4901 if (distance1>cdist1) continue; // cut on DCA linear approximation
4903 npoints = helixes[i].GetRPHIintersections(helixes[j], phase, radius,20);
4904 helixes[i].ParabolicDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta1);
4905 if (radius[0]<20&&delta1<1) continue; //intersection at vertex
4906 if (radius[0]<10&&delta1<3) continue; //intersection at vertex
4909 helixes[i].ParabolicDCA(helixes[j],phase[1][0],phase[1][1],radius[1],delta2);
4910 if (radius[1]<20&&delta2<1) continue; //intersection at vertex
4911 if (radius[1]<10&&delta2<3) continue; //intersection at vertex
4913 distance1 = TMath::Min(delta1,delta2);
4916 rkink = TMath::Sqrt(radius[0]);
4919 rkink = TMath::Sqrt(radius[1]);
4921 if (distance1>cdist2) continue;
4924 AliTPCseed * track1 = (AliTPCseed*)array->At(j);
4927 Int_t row0 = GetRowNumber(rkink);
4928 if (row0<10) continue;
4929 if (row0>150) continue;
4932 Float_t dens00=-1,dens01=-1;
4933 Float_t dens10=-1,dens11=-1;
4935 Int_t found,foundable,shared;
4936 track0->GetClusterStatistic(0,row0-5, found, foundable,shared,kFALSE);
4937 if (foundable>5) dens00 = Float_t(found)/Float_t(foundable);
4938 track0->GetClusterStatistic(row0+5,155, found, foundable,shared,kFALSE);
4939 if (foundable>5) dens01 = Float_t(found)/Float_t(foundable);
4941 track1->GetClusterStatistic(0,row0-5, found, foundable,shared,kFALSE);
4942 if (foundable>10) dens10 = Float_t(found)/Float_t(foundable);
4943 track1->GetClusterStatistic(row0+5,155, found, foundable,shared,kFALSE);
4944 if (foundable>10) dens11 = Float_t(found)/Float_t(foundable);
4946 if (dens00<dens10 && dens01<dens11) continue;
4947 if (dens00>dens10 && dens01>dens11) continue;
4948 if (TMath::Max(dens00,dens10)<0.1) continue;
4949 if (TMath::Max(dens01,dens11)<0.3) continue;
4951 if (TMath::Min(dens00,dens10)>0.6) continue;
4952 if (TMath::Min(dens01,dens11)>0.6) continue;
4955 AliTPCseed * ktrack0, *ktrack1;
4964 if (TMath::Abs(ktrack0->GetC())>5) continue; // cut on the curvature for mother particle
4965 AliExternalTrackParam paramm(*ktrack0);
4966 AliExternalTrackParam paramd(*ktrack1);
4967 if (row0>60&&ktrack1->GetReference().GetX()>90.)new (¶md) AliExternalTrackParam(ktrack1->GetReference());
4970 kink->SetMother(paramm);
4971 kink->SetDaughter(paramd);
4974 Float_t x[3] = { kink->GetPosition()[0],kink->GetPosition()[1],kink->GetPosition()[2]};
4976 fParam->Transform0to1(x,index);
4977 fParam->Transform1to2(x,index);
4978 row0 = GetRowNumber(x[0]);
4980 if (kink->GetR()<100) continue;
4981 if (kink->GetR()>240) continue;
4982 if (kink->GetPosition()[2]/kink->GetR()>AliTPCReconstructor::GetCtgRange()) continue; //out of fiducial volume
4983 if (kink->GetDistance()>cdist3) continue;
4984 Float_t dird = kink->GetDaughterP()[0]*kink->GetPosition()[0]+kink->GetDaughterP()[1]*kink->GetPosition()[1]; // rough direction estimate
4985 if (dird<0) continue;
4987 Float_t dirm = kink->GetMotherP()[0]*kink->GetPosition()[0]+kink->GetMotherP()[1]*kink->GetPosition()[1]; // rough direction estimate
4988 if (dirm<0) continue;
4989 Float_t mpt = TMath::Sqrt(kink->GetMotherP()[0]*kink->GetMotherP()[0]+kink->GetMotherP()[1]*kink->GetMotherP()[1]);
4990 if (mpt<0.2) continue;
4993 //for high momenta momentum not defined well in first iteration
4994 Double_t qt = TMath::Sin(kink->GetAngle(2))*ktrack1->GetP();
4995 if (qt>0.35) continue;
4998 kink->SetLabel(CookLabel(ktrack0,0.4,0,row0),0);
4999 kink->SetLabel(CookLabel(ktrack1,0.4,row0,160),1);
5001 kink->SetTPCDensity(dens00,0,0);
5002 kink->SetTPCDensity(dens01,0,1);
5003 kink->SetTPCDensity(dens10,1,0);
5004 kink->SetTPCDensity(dens11,1,1);
5005 kink->SetIndex(i,0);
5006 kink->SetIndex(j,1);
5009 kink->SetTPCDensity(dens10,0,0);
5010 kink->SetTPCDensity(dens11,0,1);
5011 kink->SetTPCDensity(dens00,1,0);
5012 kink->SetTPCDensity(dens01,1,1);
5013 kink->SetIndex(j,0);
5014 kink->SetIndex(i,1);
5017 if (mpt<1||kink->GetAngle(2)>0.1){
5018 // angle and densities not defined yet
5019 if (kink->GetTPCDensityFactor()<0.8) continue;
5020 if ((2-kink->GetTPCDensityFactor())*kink->GetDistance() >0.25) continue;
5021 if (kink->GetAngle(2)*ktrack0->GetP()<0.003) continue; //too small angle
5022 if (kink->GetAngle(2)>0.2&&kink->GetTPCDensityFactor()<1.15) continue;
5023 if (kink->GetAngle(2)>0.2&&kink->GetTPCDensity(0,1)>0.05) continue;
5025 Float_t criticalangle = track0->GetSigmaSnp2()+track0->GetSigmaTgl2();
5026 criticalangle+= track1->GetSigmaSnp2()+track1->GetSigmaTgl2();
5027 criticalangle= 3*TMath::Sqrt(criticalangle);
5028 if (criticalangle>0.02) criticalangle=0.02;
5029 if (kink->GetAngle(2)<criticalangle) continue;
5032 Int_t drow = Int_t(2.+0.5/(0.05+kink->GetAngle(2))); // overlap region defined
5033 Float_t shapesum =0;
5035 for ( Int_t row = row0-drow; row<row0+drow;row++){
5036 if (row<0) continue;
5037 if (row>155) continue;
5038 if (ktrack0->GetClusterPointer(row)){
5039 AliTPCTrackerPoint *point =ktrack0->GetTrackPoint(row);
5040 shapesum+=point->GetSigmaY()+point->GetSigmaZ();
5043 if (ktrack1->GetClusterPointer(row)){
5044 AliTPCTrackerPoint *point =ktrack1->GetTrackPoint(row);
5045 shapesum+=point->GetSigmaY()+point->GetSigmaZ();
5050 kink->SetShapeFactor(-1.);
5053 kink->SetShapeFactor(shapesum/sum);
5055 // esd->AddKink(kink);
5056 kinks->AddLast(kink);
5062 // sort the kinks according quality - and refit them towards vertex
5064 Int_t nkinks = kinks->GetEntriesFast();
5065 Float_t *quality = new Float_t[nkinks];
5066 Int_t *indexes = new Int_t[nkinks];
5067 AliTPCseed *mothers = new AliTPCseed[nkinks];
5068 AliTPCseed *daughters = new AliTPCseed[nkinks];
5071 for (Int_t i=0;i<nkinks;i++){
5073 AliKink *kink = (AliKink*)kinks->At(i);
5075 // refit kinks towards vertex
5077 Int_t index0 = kink->GetIndex(0);
5078 Int_t index1 = kink->GetIndex(1);
5079 AliTPCseed * ktrack0 = (AliTPCseed*)array->At(index0);
5080 AliTPCseed * ktrack1 = (AliTPCseed*)array->At(index1);
5082 Int_t sumn=ktrack0->GetNumberOfClusters()+ktrack1->GetNumberOfClusters();
5084 // Refit Kink under if too small angle
5086 if (kink->GetAngle(2)<0.05){
5087 kink->SetTPCRow0(GetRowNumber(kink->GetR()));
5088 Int_t row0 = kink->GetTPCRow0();
5089 Int_t drow = Int_t(2.+0.5/(0.05+kink->GetAngle(2)));
5092 Int_t last = row0-drow;
5093 if (last<40) last=40;
5094 if (last<ktrack0->GetFirstPoint()+25) last = ktrack0->GetFirstPoint()+25;
5095 AliTPCseed* seed0 = ReSeed(ktrack0,last,kFALSE);
5098 Int_t first = row0+drow;
5099 if (first>130) first=130;
5100 if (first>ktrack1->GetLastPoint()-25) first = TMath::Max(ktrack1->GetLastPoint()-25,30);
5101 AliTPCseed* seed1 = ReSeed(ktrack1,first,kTRUE);
5103 if (seed0 && seed1){
5104 kink->SetStatus(1,8);
5105 if (RefitKink(*seed0,*seed1,*kink)) kink->SetStatus(1,9);
5106 row0 = GetRowNumber(kink->GetR());
5107 sumn = seed0->GetNumberOfClusters()+seed1->GetNumberOfClusters();
5108 mothers[i] = *seed0;
5109 daughters[i] = *seed1;
5112 delete kinks->RemoveAt(i);
5113 if (seed0) delete seed0;
5114 if (seed1) delete seed1;
5117 if (kink->GetDistance()>0.5 || kink->GetR()<110 || kink->GetR()>240) {
5118 delete kinks->RemoveAt(i);
5119 if (seed0) delete seed0;
5120 if (seed1) delete seed1;
5128 if (kink) quality[i] = 160*((0.1+kink->GetDistance())*(2.-kink->GetTPCDensityFactor()))/(sumn+40.); //the longest -clossest will win
5130 TMath::Sort(nkinks,quality,indexes,kFALSE);
5132 //remove double find kinks
5134 for (Int_t ikink0=1;ikink0<nkinks;ikink0++){
5135 AliKink * kink0 = (AliKink*) kinks->At(indexes[ikink0]);
5136 if (!kink0) continue;
5138 for (Int_t ikink1=0;ikink1<ikink0;ikink1++){
5139 if (!kink0) continue;
5140 AliKink * kink1 = (AliKink*) kinks->At(indexes[ikink1]);
5141 if (!kink1) continue;
5142 // if not close kink continue
5143 if (TMath::Abs(kink1->GetPosition()[2]-kink0->GetPosition()[2])>10) continue;
5144 if (TMath::Abs(kink1->GetPosition()[1]-kink0->GetPosition()[1])>10) continue;
5145 if (TMath::Abs(kink1->GetPosition()[0]-kink0->GetPosition()[0])>10) continue;
5147 AliTPCseed &mother0 = mothers[indexes[ikink0]];
5148 AliTPCseed &daughter0 = daughters[indexes[ikink0]];
5149 AliTPCseed &mother1 = mothers[indexes[ikink1]];
5150 AliTPCseed &daughter1 = daughters[indexes[ikink1]];
5151 Int_t row0 = (kink0->GetTPCRow0()+kink1->GetTPCRow0())/2;
5160 for (Int_t i=0;i<row0;i++){
5161 if (mother0.GetClusterIndex(i)>0 && mother1.GetClusterIndex(i)>0){
5164 if (mother0.GetClusterIndex(i)==mother1.GetClusterIndex(i)){
5171 for (Int_t i=row0;i<158;i++){
5172 if (daughter0.GetClusterIndex(i)>0 && daughter0.GetClusterIndex(i)>0){
5175 if (mother0.GetClusterIndex(i)==mother1.GetClusterIndex(i)){
5181 Float_t ratio = Float_t(same+1)/Float_t(both+1);
5182 Float_t ratiom = Float_t(samem+1)/Float_t(bothm+1);
5183 Float_t ratiod = Float_t(samed+1)/Float_t(bothd+1);
5184 if (ratio>0.3 && ratiom>0.5 &&ratiod>0.5) {
5185 Int_t sum0 = mother0.GetNumberOfClusters()+daughter0.GetNumberOfClusters();
5186 Int_t sum1 = mother1.GetNumberOfClusters()+daughter1.GetNumberOfClusters();
5188 shared[kink0->GetIndex(0)]= kTRUE;
5189 shared[kink0->GetIndex(1)]= kTRUE;
5190 delete kinks->RemoveAt(indexes[ikink0]);
5193 shared[kink1->GetIndex(0)]= kTRUE;
5194 shared[kink1->GetIndex(1)]= kTRUE;
5195 delete kinks->RemoveAt(indexes[ikink1]);
5202 for (Int_t i=0;i<nkinks;i++){
5203 AliKink * kink = (AliKink*) kinks->At(indexes[i]);
5204 if (!kink) continue;
5205 kink->SetTPCRow0(GetRowNumber(kink->GetR()));
5206 Int_t index0 = kink->GetIndex(0);
5207 Int_t index1 = kink->GetIndex(1);
5208 if (circular[index0]||circular[index1]&&kink->GetDistance()>0.2) continue;
5209 kink->SetMultiple(usage[index0],0);
5210 kink->SetMultiple(usage[index1],1);
5211 if (kink->GetMultiple()[0]+kink->GetMultiple()[1]>2) continue;
5212 if (kink->GetMultiple()[0]+kink->GetMultiple()[1]>0 && quality[indexes[i]]>0.2) continue;
5213 if (kink->GetMultiple()[0]+kink->GetMultiple()[1]>0 && kink->GetDistance()>0.2) continue;
5214 if (circular[index0]||circular[index1]&&kink->GetDistance()>0.1) continue;
5216 AliTPCseed * ktrack0 = (AliTPCseed*)array->At(index0);
5217 AliTPCseed * ktrack1 = (AliTPCseed*)array->At(index1);
5218 if (!ktrack0 || !ktrack1) continue;
5219 Int_t index = esd->AddKink(kink);
5222 if ( ktrack0->GetKinkIndex(0)==0 && ktrack1->GetKinkIndex(0)==0) { //best kink
5223 if (mothers[indexes[i]].GetNumberOfClusters()>20 && daughters[indexes[i]].GetNumberOfClusters()>20 && (mothers[indexes[i]].GetNumberOfClusters()+daughters[indexes[i]].GetNumberOfClusters())>100){
5224 *ktrack0 = mothers[indexes[i]];
5225 *ktrack1 = daughters[indexes[i]];
5229 ktrack0->SetKinkIndex(usage[index0],-(index+1));
5230 ktrack1->SetKinkIndex(usage[index1], (index+1));
5235 // Remove tracks corresponding to shared kink's
5237 for (Int_t i=0;i<nentries;i++){
5238 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5239 if (!track0) continue;
5240 if (track0->GetKinkIndex(0)!=0) continue;
5241 if (shared[i]) delete array->RemoveAt(i);
5246 RemoveUsed2(array,0.5,0.4,30);
5248 for (Int_t i=0;i<nentries;i++){
5249 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5250 if (!track0) continue;
5251 track0->CookdEdx(0.02,0.6);
5255 for (Int_t i=0;i<nentries;i++){
5256 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5257 if (!track0) continue;
5258 if (track0->Pt()<1.4) continue;
5259 //remove double high momenta tracks - overlapped with kink candidates
5262 for (Int_t icl=track0->GetFirstPoint();icl<track0->GetLastPoint(); icl++){
5263 if (track0->GetClusterPointer(icl)!=0){
5265 if (track0->GetClusterPointer(icl)->IsUsed(10)) shared++;
5268 if (Float_t(shared+1)/Float_t(all+1)>0.5) {
5269 delete array->RemoveAt(i);
5273 if (track0->GetKinkIndex(0)!=0) continue;
5274 if (track0->GetNumberOfClusters()<80) continue;
5276 AliTPCseed *pmother = new AliTPCseed();
5277 AliTPCseed *pdaughter = new AliTPCseed();
5278 AliKink *pkink = new AliKink;
5280 AliTPCseed & mother = *pmother;
5281 AliTPCseed & daughter = *pdaughter;
5282 AliKink & kink = *pkink;
5283 if (CheckKinkPoint(track0,mother,daughter, kink)){
5284 if (mother.GetNumberOfClusters()<30||daughter.GetNumberOfClusters()<20) {
5288 continue; //too short tracks
5290 if (mother.Pt()<1.4) {
5296 Int_t row0= kink.GetTPCRow0();
5297 if (kink.GetDistance()>0.5 || kink.GetR()<110. || kink.GetR()>240.) {
5304 Int_t index = esd->AddKink(&kink);
5305 mother.SetKinkIndex(0,-(index+1));
5306 daughter.SetKinkIndex(0,index+1);
5307 if (mother.GetNumberOfClusters()>50) {
5308 delete array->RemoveAt(i);
5309 array->AddAt(new AliTPCseed(mother),i);
5312 array->AddLast(new AliTPCseed(mother));
5314 array->AddLast(new AliTPCseed(daughter));
5315 for (Int_t icl=0;icl<row0;icl++) {
5316 if (mother.GetClusterPointer(icl)) mother.GetClusterPointer(icl)->Use(20);
5319 for (Int_t icl=row0;icl<158;icl++) {
5320 if (daughter.GetClusterPointer(icl)) daughter.GetClusterPointer(icl)->Use(20);
5329 delete [] daughters;
5351 printf("Ncandidates=\t%d\t%d\t%d\t%d\n",esd->GetNumberOfKinks(),ncandidates,ntracks,nall);
5355 void AliTPCtrackerMI::FindV0s(TObjArray * array, AliESDEvent *esd)
5361 TObjArray *tpcv0s = new TObjArray(100000);
5362 Int_t nentries = array->GetEntriesFast();
5363 AliHelix *helixes = new AliHelix[nentries];
5364 Int_t *sign = new Int_t[nentries];
5365 Float_t *alpha = new Float_t[nentries];
5366 Float_t *z0 = new Float_t[nentries];
5367 Float_t *dca = new Float_t[nentries];
5368 Float_t *sdcar = new Float_t[nentries];
5369 Float_t *cdcar = new Float_t[nentries];
5370 Float_t *pulldcar = new Float_t[nentries];
5371 Float_t *pulldcaz = new Float_t[nentries];
5372 Float_t *pulldca = new Float_t[nentries];
5373 Bool_t *isPrim = new Bool_t[nentries];
5374 const AliESDVertex * primvertex = esd->GetVertex();
5375 Double_t zvertex = primvertex->GetZv();
5377 // nentries = array->GetEntriesFast();
5379 for (Int_t i=0;i<nentries;i++){
5382 AliTPCseed* track = (AliTPCseed*)array->At(i);
5383 if (!track) continue;
5384 track->GetV0Indexes()[0] = 0; //rest v0 indexes
5385 track->GetV0Indexes()[1] = 0; //rest v0 indexes
5386 track->GetV0Indexes()[2] = 0; //rest v0 indexes
5388 alpha[i] = track->GetAlpha();
5389 new (&helixes[i]) AliHelix(*track);
5391 helixes[i].Evaluate(0,xyz);
5392 sign[i] = (track->GetC()>0) ? -1:1;
5396 if (track->GetProlongation(0,y,z)) z0[i] = z;
5397 dca[i] = track->GetD(0,0);
5399 // dca error parrameterezation + pulls
5401 sdcar[i] = TMath::Sqrt(0.150*0.150+(100*track->GetC())*(100*track->GetC()));
5402 if (TMath::Abs(track->GetTgl())>1) sdcar[i]*=2.5;
5403 cdcar[i] = TMath::Exp((TMath::Abs(track->GetC())-0.0106)*525.3);
5404 pulldcar[i] = (dca[i]-cdcar[i])/sdcar[i];
5405 pulldcaz[i] = (z0[i]-zvertex)/sdcar[i];
5406 pulldca[i] = TMath::Sqrt(pulldcar[i]*pulldcar[i]+pulldcaz[i]*pulldcaz[i]);
5407 if (track->TPCrPID(1)+track->TPCrPID(2)+track->TPCrPID(3)>0.5) {
5408 if (pulldca[i]<3.) isPrim[i]=kTRUE; //pion, muon and Kaon 3 sigma cut
5410 if (track->TPCrPID(4)>0.5) {
5411 if (pulldca[i]<0.5) isPrim[i]=kTRUE; //proton 0.5 sigma cut
5413 if (track->TPCrPID(0)>0.4) {
5414 isPrim[i]=kFALSE; //electron no sigma cut
5421 Int_t ncandidates =0;
5424 Double_t phase[2][2],radius[2];
5430 TTreeSRedirector &cstream = *fDebugStreamer;
5431 Float_t fprimvertex[3]={GetX(),GetY(),GetZ()};
5433 Double_t cradius0 = 10*10;
5434 Double_t cradius1 = 200*200;
5437 Double_t cpointAngle = 0.95;
5439 Double_t delta[2]={10000,10000};
5440 for (Int_t i =0;i<nentries;i++){
5441 if (sign[i]==0) continue;
5442 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5443 if (!track0) continue;
5444 if (AliTPCReconstructor::StreamLevel()>1){
5449 "zvertex="<<zvertex<<
5450 "sdcar0="<<sdcar[i]<<
5451 "cdcar0="<<cdcar[i]<<
5452 "pulldcar0="<<pulldcar[i]<<
5453 "pulldcaz0="<<pulldcaz[i]<<
5454 "pulldca0="<<pulldca[i]<<
5455 "isPrim="<<isPrim[i]<<
5459 if (track0->GetSigned1Pt()<0) continue;
5460 if (track0->GetKinkIndex(0)>0||isPrim[i]) continue; //daughter kink
5462 if (TMath::Abs(helixes[i].GetHelix(4))<0.000000001) continue;
5467 for (Int_t j =0;j<nentries;j++){
5468 AliTPCseed * track1 = (AliTPCseed*)array->At(j);
5469 if (!track1) continue;
5470 if (track1->GetKinkIndex(0)>0 || isPrim[j]) continue; //daughter kink
5471 if (sign[j]*sign[i]>0) continue;
5472 if (TMath::Abs(helixes[j].GetHelix(4))<0.000001) continue;
5473 if (track0->GetCircular()+track1->GetCircular()>1) continue; //circling -returning track
5476 // DCA to prim vertex cut
5482 Int_t npoints = helixes[i].GetRPHIintersections(helixes[j], phase, radius,cdist2);
5483 if (npoints<1) continue;
5487 if (radius[0]<cradius0||radius[0]>cradius1) continue;
5488 helixes[i].LinearDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta[0]);
5489 if (delta[0]>cdist1) continue;
5492 if (TMath::Max(radius[0],radius[1])<cradius0|| TMath::Min(radius[0],radius[1])>cradius1) continue;
5493 helixes[i].LinearDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta[0]);
5494 helixes[i].LinearDCA(helixes[j],phase[1][0],phase[1][1],radius[1],delta[1]);
5495 if (delta[1]<delta[0]) iclosest=1;
5496 if (delta[iclosest]>cdist1) continue;
5498 helixes[i].ParabolicDCA(helixes[j],phase[iclosest][0],phase[iclosest][1],radius[iclosest],delta[iclosest]);
5499 if (radius[iclosest]<cradius0 || radius[iclosest]>cradius1 || delta[iclosest]>cdist1) continue;
5501 Double_t pointAngle = helixes[i].GetPointAngle(helixes[j],phase[iclosest],fprimvertex);
5502 if (pointAngle<cpointAngle) continue;
5504 Bool_t isGamma = kFALSE;
5505 vertex.SetParamP(*track0); //track0 - plus
5506 vertex.SetParamN(*track1); //track1 - minus
5507 vertex.Update(fprimvertex);
5508 if (track0->TPCrPID(0)>0.3&&track1->TPCrPID(0)>0.3&&vertex.GetAnglep()[2]<0.15) isGamma=kTRUE; // gamma conversion candidate
5509 Double_t pointAngle2 = vertex.GetV0CosineOfPointingAngle();
5511 if (vertex.GetV0CosineOfPointingAngle()<cpointAngle && (!isGamma)) continue;// point angle cut
5512 //Bo: if (vertex.GetDist2()>2&&(!isGamma)) continue; // point angle cut
5513 if (vertex.GetDcaV0Daughters()>2&&(!isGamma)) continue;//Bo: // point angle cut
5514 Float_t sigmae = 0.15*0.15;
5515 if (vertex.GetRr()<80)
5516 sigmae += (sdcar[i]*sdcar[i]+sdcar[j]*sdcar[j])*(1.-vertex.GetRr()/80.)*(1.-vertex.GetRr()/80.);
5517 sigmae+= TMath::Sqrt(sigmae);
5518 //Bo: if (vertex.GetDist2()/sigmae>3.&&(!isGamma)) continue;
5519 if (vertex.GetDcaV0Daughters()/sigmae>3.&&(!isGamma)) continue;
5520 Float_t densb0=0,densb1=0,densa0=0,densa1=0;
5521 Int_t row0 = GetRowNumber(vertex.GetRr());
5523 //Bo: if (vertex.GetDist2()>0.2) continue;
5524 if (vertex.GetDcaV0Daughters()>0.2) continue;
5525 densb0 = track0->Density2(0,row0-5);
5526 densb1 = track1->Density2(0,row0-5);
5527 if (densb0>0.3|| densb1>0.3) continue; //clusters before vertex
5528 densa0 = track0->Density2(row0+5,row0+40);
5529 densa1 = track1->Density2(row0+5,row0+40);
5530 if ((densa0<0.4|| densa1<0.4)&&(!isGamma)) continue; //missing clusters after vertex
5533 densa0 = track0->Density2(0,40); //cluster density
5534 densa1 = track1->Density2(0,40); //cluster density
5535 if ((vertex.GetRr()<80&&densa0+densa1<1.)&&(!isGamma)) continue;
5537 //Bo: vertex.SetLab(0,track0->GetLabel());
5538 //Bo: vertex.SetLab(1,track1->GetLabel());
5539 vertex.SetChi2After((densa0+densa1)*0.5);
5540 vertex.SetChi2Before((densb0+densb1)*0.5);
5541 vertex.SetIndex(0,i);
5542 vertex.SetIndex(1,j);
5543 //Bo: vertex.SetStatus(1); // TPC v0 candidate
5544 vertex.SetOnFlyStatus(2);//Bo: // TPC v0 candidate
5545 //Bo: vertex.SetRp(track0->TPCrPIDs());
5546 //Bo: vertex.SetRm(track1->TPCrPIDs());
5547 tpcv0s->AddLast(new AliESDv0(vertex));
5550 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
5551 Double_t radiusm= (delta[0]<delta[1])? TMath::Sqrt(radius[0]):TMath::Sqrt(radius[1]);
5552 Double_t deltam= (delta[0]<delta[1])? TMath::Sqrt(delta[0]):TMath::Sqrt(delta[1]);
5553 if (AliTPCReconstructor::StreamLevel()>1) {
5554 Int_t lab0=track0->GetLabel();
5555 Int_t lab1=track1->GetLabel();
5556 Char_t c0=track0->GetCircular();
5557 Char_t c1=track1->GetCircular();
5560 "vertex.="<<&vertex<<
5563 "Helix0.="<<&helixes[i]<<
5566 "Helix1.="<<&helixes[j]<<
5567 "pointAngle="<<pointAngle<<
5568 "pointAngle2="<<pointAngle2<<
5573 "zvertex="<<zvertex<<
5576 "npoints="<<npoints<<
5577 "radius0="<<radius[0]<<
5578 "delta0="<<delta[0]<<
5579 "radius1="<<radius[1]<<
5580 "delta1="<<delta[1]<<
5581 "radiusm="<<radiusm<<
5583 "sdcar0="<<sdcar[i]<<
5584 "sdcar1="<<sdcar[j]<<
5585 "cdcar0="<<cdcar[i]<<
5586 "cdcar1="<<cdcar[j]<<
5587 "pulldcar0="<<pulldcar[i]<<
5588 "pulldcar1="<<pulldcar[j]<<
5589 "pulldcaz0="<<pulldcaz[i]<<
5590 "pulldcaz1="<<pulldcaz[j]<<
5591 "pulldca0="<<pulldca[i]<<
5592 "pulldca1="<<pulldca[j]<<
5602 Float_t *quality = new Float_t[ncandidates];
5603 Int_t *indexes = new Int_t[ncandidates];
5605 for (Int_t i=0;i<ncandidates;i++){
5607 AliESDv0 *v0 = (AliESDv0*)tpcv0s->At(i);
5608 quality[i] = 1./(1.00001-v0->GetV0CosineOfPointingAngle()); //base point angle
5609 // quality[i] /= (0.5+v0->GetDist2());
5610 // quality[i] *= v0->GetChi2After(); //density factor
5612 Int_t index0 = v0->GetIndex(0);
5613 Int_t index1 = v0->GetIndex(1);
5614 //Bo: Double_t minpulldca = TMath::Min(2.+pulldca[v0->GetIndex(0)],(2.+pulldca[v0->GetIndex(1)]) ); //pull
5615 Double_t minpulldca = TMath::Min(2.+pulldca[index0],(2.+pulldca[index1]) );//Bo:
5619 AliTPCseed * track0 = (AliTPCseed*)array->At(index0);
5620 AliTPCseed * track1 = (AliTPCseed*)array->At(index1);
5621 if (track0->TPCrPID(0)>0.3&&track1->TPCrPID(0)>0.3&&v0->GetAnglep()[2]<0.15) quality[i]+=1000000; // gamma conversion candidate
5622 if (track0->TPCrPID(4)>0.9||track1->TPCrPID(4)>0.9&&minpulldca>4) quality[i]*=10; // lambda candidate candidate
5625 TMath::Sort(ncandidates,quality,indexes,kTRUE);
5628 for (Int_t i=0;i<ncandidates;i++){
5629 AliESDv0 * v0 = (AliESDv0*)tpcv0s->At(indexes[i]);
5631 Int_t index0 = v0->GetIndex(0);
5632 Int_t index1 = v0->GetIndex(1);
5633 AliTPCseed * track0 = (AliTPCseed*)array->At(index0);
5634 AliTPCseed * track1 = (AliTPCseed*)array->At(index1);
5635 if (!track0||!track1) {
5639 Bool_t accept =kTRUE; //default accept
5640 Int_t *v0indexes0 = track0->GetV0Indexes();
5641 Int_t *v0indexes1 = track1->GetV0Indexes();
5643 Int_t order0 = (v0indexes0[0]!=0) ? 1:0;
5644 Int_t order1 = (v0indexes1[0]!=0) ? 1:0;
5645 if (v0indexes0[1]!=0) order0 =2;
5646 if (v0indexes1[1]!=0) order1 =2;
5648 if (v0indexes0[2]!=0) {order0=3; accept=kFALSE;}
5649 if (v0indexes0[2]!=0) {order1=3; accept=kFALSE;}
5651 AliESDv0 * v02 = v0;
5653 //Bo: v0->SetOrder(0,order0);
5654 //Bo: v0->SetOrder(1,order1);
5655 //Bo: v0->SetOrder(1,order0+order1);
5656 v0->SetOnFlyStatus(kTRUE);
5657 Int_t index = esd->AddV0(v0);
5658 v02 = esd->GetV0(index);
5659 v0indexes0[order0]=index;
5660 v0indexes1[order1]=index;
5664 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
5665 if (AliTPCReconstructor::StreamLevel()>1) {
5666 Int_t lab0=track0->GetLabel();
5667 Int_t lab1=track1->GetLabel();
5676 "dca0="<<dca[index0]<<
5677 "dca1="<<dca[index1]<<
5681 "quality="<<quality[i]<<
5682 "pulldca0="<<pulldca[index0]<<
5683 "pulldca1="<<pulldca[index1]<<
5707 printf("TPC V0 finder : naccepted\t%d\tncandidates\t%d\tntracks\t%d\tnall\t%d\n",naccepted,ncandidates,ntracks,nall);
5711 Int_t AliTPCtrackerMI::RefitKink(AliTPCseed &mother, AliTPCseed &daughter, AliESDkink &knk)
5714 // refit kink towards to the vertex
5717 AliKink &kink=(AliKink &)knk;
5719 Int_t row0 = GetRowNumber(kink.GetR());
5720 FollowProlongation(mother,0);
5721 mother.Reset(kFALSE);
5723 FollowProlongation(daughter,row0);
5724 daughter.Reset(kFALSE);
5725 FollowBackProlongation(daughter,158);
5726 daughter.Reset(kFALSE);
5727 Int_t first = TMath::Max(row0-20,30);
5728 Int_t last = TMath::Min(row0+20,140);
5730 const Int_t kNdiv =5;
5731 AliTPCseed param0[kNdiv]; // parameters along the track
5732 AliTPCseed param1[kNdiv]; // parameters along the track
5733 AliKink kinks[kNdiv]; // corresponding kink parameters
5736 for (Int_t irow=0; irow<kNdiv;irow++){
5737 rows[irow] = first +((last-first)*irow)/(kNdiv-1);
5739 // store parameters along the track
5741 for (Int_t irow=0;irow<kNdiv;irow++){
5742 FollowBackProlongation(mother, rows[irow]);
5743 FollowProlongation(daughter,rows[kNdiv-1-irow]);
5744 param0[irow] = mother;
5745 param1[kNdiv-1-irow] = daughter;
5749 for (Int_t irow=0; irow<kNdiv-1;irow++){
5750 if (param0[irow].GetNumberOfClusters()<kNdiv||param1[irow].GetNumberOfClusters()<kNdiv) continue;
5751 kinks[irow].SetMother(param0[irow]);
5752 kinks[irow].SetDaughter(param1[irow]);
5753 kinks[irow].Update();
5756 // choose kink with best "quality"
5758 Double_t mindist = 10000;
5759 for (Int_t irow=0;irow<kNdiv;irow++){
5760 if (param0[irow].GetNumberOfClusters()<20||param1[irow].GetNumberOfClusters()<20) continue;
5761 if (TMath::Abs(kinks[irow].GetR())>240.) continue;
5762 if (TMath::Abs(kinks[irow].GetR())<100.) continue;
5764 Float_t normdist = TMath::Abs(param0[irow].GetX()-kinks[irow].GetR())*(0.1+kink.GetDistance());
5765 normdist/= (param0[irow].GetNumberOfClusters()+param1[irow].GetNumberOfClusters()+40.);
5766 if (normdist < mindist){
5772 if (index==-1) return 0;
5775 param0[index].Reset(kTRUE);
5776 FollowProlongation(param0[index],0);
5778 mother = param0[index];
5779 daughter = param1[index]; // daughter in vertex
5781 kink.SetMother(mother);
5782 kink.SetDaughter(daughter);
5784 kink.SetTPCRow0(GetRowNumber(kink.GetR()));
5785 kink.SetTPCncls(param0[index].GetNumberOfClusters(),0);
5786 kink.SetTPCncls(param1[index].GetNumberOfClusters(),1);
5787 kink.SetLabel(CookLabel(&mother,0.4, 0,kink.GetTPCRow0()),0);
5788 kink.SetLabel(CookLabel(&daughter,0.4, kink.GetTPCRow0(),160),1);
5789 mother.SetLabel(kink.GetLabel(0));
5790 daughter.SetLabel(kink.GetLabel(1));
5796 void AliTPCtrackerMI::UpdateKinkQualityM(AliTPCseed * seed){
5798 // update Kink quality information for mother after back propagation
5800 if (seed->GetKinkIndex(0)>=0) return;
5801 for (Int_t ikink=0;ikink<3;ikink++){
5802 Int_t index = seed->GetKinkIndex(ikink);
5803 if (index>=0) break;
5804 index = TMath::Abs(index)-1;
5805 AliESDkink * kink = fEvent->GetKink(index);
5806 kink->SetTPCDensity(-1,0,0);
5807 kink->SetTPCDensity(1,0,1);
5809 Int_t row0 = kink->GetTPCRow0() - 2 - Int_t( 0.5/ (0.05+kink->GetAngle(2)));
5810 if (row0<15) row0=15;
5812 Int_t row1 = kink->GetTPCRow0() + 2 + Int_t( 0.5/ (0.05+kink->GetAngle(2)));
5813 if (row1>145) row1=145;
5815 Int_t found,foundable,shared;
5816 seed->GetClusterStatistic(0,row0, found, foundable,shared,kFALSE);
5817 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),0,0);
5818 seed->GetClusterStatistic(row1,155, found, foundable,shared,kFALSE);
5819 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),0,1);
5824 void AliTPCtrackerMI::UpdateKinkQualityD(AliTPCseed * seed){
5826 // update Kink quality information for daughter after refit
5828 if (seed->GetKinkIndex(0)<=0) return;
5829 for (Int_t ikink=0;ikink<3;ikink++){
5830 Int_t index = seed->GetKinkIndex(ikink);
5831 if (index<=0) break;
5832 index = TMath::Abs(index)-1;
5833 AliESDkink * kink = fEvent->GetKink(index);
5834 kink->SetTPCDensity(-1,1,0);
5835 kink->SetTPCDensity(-1,1,1);
5837 Int_t row0 = kink->GetTPCRow0() -2 - Int_t( 0.5/ (0.05+kink->GetAngle(2)));
5838 if (row0<15) row0=15;
5840 Int_t row1 = kink->GetTPCRow0() +2 + Int_t( 0.5/ (0.05+kink->GetAngle(2)));
5841 if (row1>145) row1=145;
5843 Int_t found,foundable,shared;
5844 seed->GetClusterStatistic(0,row0, found, foundable,shared,kFALSE);
5845 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),1,0);
5846 seed->GetClusterStatistic(row1,155, found, foundable,shared,kFALSE);
5847 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),1,1);
5853 Int_t AliTPCtrackerMI::CheckKinkPoint(AliTPCseed*seed,AliTPCseed &mother, AliTPCseed &daughter, AliESDkink &knk)
5856 // check kink point for given track
5857 // if return value=0 kink point not found
5858 // otherwise seed0 correspond to mother particle
5859 // seed1 correspond to daughter particle
5860 // kink parameter of kink point
5861 AliKink &kink=(AliKink &)knk;
5863 Int_t middlerow = (seed->GetFirstPoint()+seed->GetLastPoint())/2;
5864 Int_t first = seed->GetFirstPoint();
5865 Int_t last = seed->GetLastPoint();
5866 if (last-first<20) return 0; // shortest length - 2*30 = 60 pad-rows
5869 AliTPCseed *seed1 = ReSeed(seed,middlerow+20, kTRUE); //middle of chamber
5870 if (!seed1) return 0;
5871 FollowProlongation(*seed1,seed->GetLastPoint()-20);
5872 seed1->Reset(kTRUE);
5873 FollowProlongation(*seed1,158);
5874 seed1->Reset(kTRUE);
5875 last = seed1->GetLastPoint();
5877 AliTPCseed *seed0 = new AliTPCseed(*seed);
5878 seed0->Reset(kFALSE);
5881 AliTPCseed param0[20]; // parameters along the track
5882 AliTPCseed param1[20]; // parameters along the track
5883 AliKink kinks[20]; // corresponding kink parameters
5885 for (Int_t irow=0; irow<20;irow++){
5886 rows[irow] = first +((last-first)*irow)/19;
5888 // store parameters along the track
5890 for (Int_t irow=0;irow<20;irow++){
5891 FollowBackProlongation(*seed0, rows[irow]);
5892 FollowProlongation(*seed1,rows[19-irow]);
5893 param0[irow] = *seed0;
5894 param1[19-irow] = *seed1;
5898 for (Int_t irow=0; irow<19;irow++){
5899 kinks[irow].SetMother(param0[irow]);
5900 kinks[irow].SetDaughter(param1[irow]);
5901 kinks[irow].Update();
5904 // choose kink with biggest change of angle
5906 Double_t maxchange= 0;
5907 for (Int_t irow=1;irow<19;irow++){
5908 if (TMath::Abs(kinks[irow].GetR())>240.) continue;
5909 if (TMath::Abs(kinks[irow].GetR())<110.) continue;
5910 Float_t quality = TMath::Abs(kinks[irow].GetAngle(2))/(3.+TMath::Abs(kinks[irow].GetR()-param0[irow].GetX()));
5911 if ( quality > maxchange){
5912 maxchange = quality;
5919 if (index<0) return 0;
5921 Int_t row0 = GetRowNumber(kinks[index].GetR()); //row 0 estimate
5922 seed0 = new AliTPCseed(param0[index]);
5923 seed1 = new AliTPCseed(param1[index]);
5924 seed0->Reset(kFALSE);
5925 seed1->Reset(kFALSE);
5926 seed0->ResetCovariance(10.);
5927 seed1->ResetCovariance(10.);
5928 FollowProlongation(*seed0,0);
5929 FollowBackProlongation(*seed1,158);
5930 mother = *seed0; // backup mother at position 0
5931 seed0->Reset(kFALSE);
5932 seed1->Reset(kFALSE);
5933 seed0->ResetCovariance(10.);
5934 seed1->ResetCovariance(10.);
5936 first = TMath::Max(row0-20,0);
5937 last = TMath::Min(row0+20,158);
5939 for (Int_t irow=0; irow<20;irow++){
5940 rows[irow] = first +((last-first)*irow)/19;
5942 // store parameters along the track
5944 for (Int_t irow=0;irow<20;irow++){
5945 FollowBackProlongation(*seed0, rows[irow]);
5946 FollowProlongation(*seed1,rows[19-irow]);
5947 param0[irow] = *seed0;
5948 param1[19-irow] = *seed1;
5952 for (Int_t irow=0; irow<19;irow++){
5953 kinks[irow].SetMother(param0[irow]);
5954 kinks[irow].SetDaughter(param1[irow]);
5955 // param0[irow].Dump();
5956 //param1[irow].Dump();
5957 kinks[irow].Update();
5960 // choose kink with biggest change of angle
5963 for (Int_t irow=0;irow<20;irow++){
5964 if (TMath::Abs(kinks[irow].GetR())>250.) continue;
5965 if (TMath::Abs(kinks[irow].GetR())<90.) continue;
5966 Float_t quality = TMath::Abs(kinks[irow].GetAngle(2))/(3.+TMath::Abs(kinks[irow].GetR()-param0[irow].GetX()));
5967 if ( quality > maxchange){
5968 maxchange = quality;
5975 if (index==-1 || param0[index].GetNumberOfClusters()+param1[index].GetNumberOfClusters()<100){
5980 // Float_t anglesigma = TMath::Sqrt(param0[index].fC22+param0[index].fC33+param1[index].fC22+param1[index].fC33);
5982 kink.SetMother(param0[index]);
5983 kink.SetDaughter(param1[index]);
5985 row0 = GetRowNumber(kink.GetR());
5986 kink.SetTPCRow0(row0);
5987 kink.SetLabel(CookLabel(seed0,0.5,0,row0),0);
5988 kink.SetLabel(CookLabel(seed1,0.5,row0,158),1);
5989 kink.SetIndex(-10,0);
5990 kink.SetIndex(int(param0[index].GetNumberOfClusters()+param1[index].GetNumberOfClusters()),1);
5991 kink.SetTPCncls(param0[index].GetNumberOfClusters(),0);
5992 kink.SetTPCncls(param1[index].GetNumberOfClusters(),1);
5995 // new (&mother) AliTPCseed(param0[index]);
5996 daughter = param1[index];
5997 daughter.SetLabel(kink.GetLabel(1));
5998 param0[index].Reset(kTRUE);
5999 FollowProlongation(param0[index],0);
6000 mother = param0[index];
6001 mother.SetLabel(kink.GetLabel(0));
6011 AliTPCseed* AliTPCtrackerMI::ReSeed(AliTPCseed *t)
6014 // reseed - refit - track
6017 // Int_t last = fSectors->GetNRows()-1;
6019 if (fSectors == fOuterSec){
6020 first = TMath::Max(first, t->GetFirstPoint()-fInnerSec->GetNRows());
6024 first = t->GetFirstPoint();
6026 AliTPCseed * seed = MakeSeed(t,0.1,0.5,0.9);
6027 FollowBackProlongation(*t,fSectors->GetNRows()-1);
6029 FollowProlongation(*t,first);
6039 //_____________________________________________________________________________
6040 Int_t AliTPCtrackerMI::ReadSeeds(const TFile *inp) {
6041 //-----------------------------------------------------------------
6042 // This function reades track seeds.
6043 //-----------------------------------------------------------------
6044 TDirectory *savedir=gDirectory;
6046 TFile *in=(TFile*)inp;
6047 if (!in->IsOpen()) {
6048 cerr<<"AliTPCtrackerMI::ReadSeeds(): input file is not open !\n";
6053 TTree *seedTree=(TTree*)in->Get("Seeds");
6055 cerr<<"AliTPCtrackerMI::ReadSeeds(): ";
6056 cerr<<"can't get a tree with track seeds !\n";
6059 AliTPCtrack *seed=new AliTPCtrack;
6060 seedTree->SetBranchAddress("tracks",&seed);
6062 if (fSeeds==0) fSeeds=new TObjArray(15000);
6064 Int_t n=(Int_t)seedTree->GetEntries();
6065 for (Int_t i=0; i<n; i++) {
6066 seedTree->GetEvent(i);
6067 fSeeds->AddLast(new AliTPCseed(*seed/*,seed->GetAlpha()*/));
6076 Int_t AliTPCtrackerMI::Clusters2Tracks (AliESDEvent *esd)
6079 if (fSeeds) DeleteSeeds();
6082 if (!fSeeds) return 1;
6089 //_____________________________________________________________________________
6090 Int_t AliTPCtrackerMI::Clusters2Tracks() {
6091 //-----------------------------------------------------------------
6092 // This is a track finder.
6093 //-----------------------------------------------------------------
6094 TDirectory *savedir=gDirectory;
6098 fSeeds = Tracking();
6101 Info("Clusters2Tracks","Time for tracking: \t");timer.Print();timer.Start();
6103 //activate again some tracks
6104 for (Int_t i=0; i<fSeeds->GetEntriesFast(); i++) {
6105 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
6107 Int_t nc=t.GetNumberOfClusters();
6109 delete fSeeds->RemoveAt(i);
6113 if (pt->GetRemoval()==10) {
6114 if (pt->GetDensityFirst(20)>0.8 || pt->GetDensityFirst(30)>0.8 || pt->GetDensityFirst(40)>0.7)
6115 pt->Desactivate(10); // make track again active
6117 pt->Desactivate(20);
6118 delete fSeeds->RemoveAt(i);
6123 RemoveUsed2(fSeeds,0.85,0.85,0);
6124 if (AliTPCReconstructor::GetRecoParam()->GetDoKinks()) FindKinks(fSeeds,fEvent);
6125 //FindCurling(fSeeds, fEvent,0);
6126 if (AliTPCReconstructor::StreamLevel()>2) FindMultiMC(fSeeds, fEvent,0); // find multi found tracks
6127 RemoveUsed2(fSeeds,0.5,0.4,20);
6128 FindSplitted(fSeeds, fEvent,0); // find multi found tracks
6131 // // refit short tracks
6133 Int_t nseed=fSeeds->GetEntriesFast();
6136 for (Int_t i=0; i<nseed; i++) {
6137 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
6139 Int_t nc=t.GetNumberOfClusters();
6141 delete fSeeds->RemoveAt(i);
6144 CookLabel(pt,0.1); //For comparison only
6145 //if ((pt->IsActive() || (pt->fRemoval==10) )&& nc>50 &&pt->GetNumberOfClusters()>0.4*pt->fNFoundable){
6146 if ((pt->IsActive() || (pt->GetRemoval()==10) )){
6148 if (fDebug>0) cerr<<found<<'\r';
6152 delete fSeeds->RemoveAt(i);
6156 //RemoveOverlap(fSeeds,0.99,7,kTRUE);
6158 //RemoveUsed(fSeeds,0.9,0.9,6);
6160 nseed=fSeeds->GetEntriesFast();
6162 for (Int_t i=0; i<nseed; i++) {
6163 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
6165 Int_t nc=t.GetNumberOfClusters();
6167 delete fSeeds->RemoveAt(i);
6171 t.CookdEdx(0.02,0.6);
6172 // CheckKinkPoint(&t,0.05);
6173 //if ((pt->IsActive() || (pt->fRemoval==10) )&& nc>50 &&pt->GetNumberOfClusters()>0.4*pt->fNFoundable){
6174 if ((pt->IsActive() || (pt->GetRemoval()==10) )){
6182 delete fSeeds->RemoveAt(i);
6183 //AliTPCseed * seed1 = ReSeed(pt,0.05,0.5,1);
6185 // FollowProlongation(*seed1,0);
6186 // Int_t n = seed1->GetNumberOfClusters();
6187 // printf("fP4\t%f\t%f\n",seed1->GetC(),pt->GetC());
6188 // printf("fN\t%d\t%d\n", seed1->GetNumberOfClusters(),pt->GetNumberOfClusters());
6191 //AliTPCseed * seed2 = ReSeed(pt,0.95,0.5,0.05);
6195 SortTracks(fSeeds, 1);
6199 PrepareForBackProlongation(fSeeds,5.);
6200 PropagateBack(fSeeds);
6201 printf("Time for back propagation: \t");timer.Print();timer.Start();
6205 PrepareForProlongation(fSeeds,5.);
6206 PropagateForward2(fSeeds);
6208 printf("Time for FORWARD propagation: \t");timer.Print();timer.Start();
6209 // RemoveUsed(fSeeds,0.7,0.7,6);
6210 //RemoveOverlap(fSeeds,0.9,7,kTRUE);
6212 nseed=fSeeds->GetEntriesFast();
6214 for (Int_t i=0; i<nseed; i++) {
6215 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
6217 Int_t nc=t.GetNumberOfClusters();
6219 delete fSeeds->RemoveAt(i);
6222 t.CookdEdx(0.02,0.6);
6223 // CookLabel(pt,0.1); //For comparison only
6224 //if ((pt->IsActive() || (pt->fRemoval==10) )&& nc>50 &&pt->GetNumberOfClusters()>0.4*pt->fNFoundable){
6225 if ((pt->IsActive() || (pt->fRemoval==10) )){
6226 cerr<<found++<<'\r';
6229 delete fSeeds->RemoveAt(i);
6234 // fNTracks = found;
6236 Info("Clusters2Tracks","Time for overlap removal, track writing and dedx cooking: \t"); timer.Print();timer.Start();
6239 // cerr<<"Number of found tracks : "<<"\t"<<found<<endl;
6240 Info("Clusters2Tracks","Number of found tracks %d",found);
6242 // UnloadClusters();
6247 void AliTPCtrackerMI::Tracking(TObjArray * arr)
6250 // tracking of the seeds
6253 fSectors = fOuterSec;
6254 ParallelTracking(arr,150,63);
6255 fSectors = fOuterSec;
6256 ParallelTracking(arr,63,0);
6259 TObjArray * AliTPCtrackerMI::Tracking(Int_t seedtype, Int_t i1, Int_t i2, Float_t cuts[4], Float_t dy, Int_t dsec)
6264 TObjArray * arr = new TObjArray;
6266 fSectors = fOuterSec;
6269 for (Int_t sec=0;sec<fkNOS;sec++){
6270 if (seedtype==3) MakeSeeds3(arr,sec,i1,i2,cuts,dy, dsec);
6271 if (seedtype==4) MakeSeeds5(arr,sec,i1,i2,cuts,dy);
6272 if (seedtype==2) MakeSeeds2(arr,sec,i1,i2,cuts,dy);
6275 Info("Tracking","\nSeeding - %d\t%d\t%d\t%d\n",seedtype,i1,i2,arr->GetEntriesFast());
6287 TObjArray * AliTPCtrackerMI::Tracking()
6291 if (AliTPCReconstructor::GetRecoParam()->GetSpecialSeeding()) return TrackingSpecial();
6294 Int_t nup=fOuterSec->GetNRows()+fInnerSec->GetNRows();
6296 TObjArray * seeds = new TObjArray;
6305 Float_t fnumber = 3.0;
6306 Float_t fdensity = 3.0;
6311 for (Int_t delta = 0; delta<18; delta+=6){
6315 arr = Tracking(3,nup-1-delta,nup-1-delta-gap,cuts,-1,1);
6316 SumTracks(seeds,arr);
6317 SignClusters(seeds,fnumber,fdensity);
6319 for (Int_t i=2;i<6;i+=2){
6320 // seed high pt tracks
6323 arr = Tracking(3,nup-i-delta,nup-i-delta-gap,cuts,-1,0);
6324 SumTracks(seeds,arr);
6325 SignClusters(seeds,fnumber,fdensity);
6330 // RemoveUsed(seeds,0.9,0.9,1);
6331 // UnsignClusters();
6332 // SignClusters(seeds,fnumber,fdensity);
6336 for (Int_t delta = 20; delta<120; delta+=10){
6338 // seed high pt tracks
6342 arr = Tracking(3,nup-delta,nup-delta-gap,cuts,-1);
6343 SumTracks(seeds,arr);
6344 SignClusters(seeds,fnumber,fdensity);
6349 arr = Tracking(3,nup-delta-5,nup-delta-5-gap,cuts,-1);
6350 SumTracks(seeds,arr);
6351 SignClusters(seeds,fnumber,fdensity);
6362 Info("Tracking()","\n\nPrimary seeding\t%d\n\n",seeds->GetEntriesFast());
6366 // RemoveUsed(seeds,0.75,0.75,1);
6368 //SignClusters(seeds,fnumber,fdensity);
6377 arr = Tracking(4,nup-1,nup-1-gap,cuts,-1);
6378 SumTracks(seeds,arr);
6379 SignClusters(seeds,fnumber,fdensity);
6381 arr = Tracking(4,nup-2,nup-2-gap,cuts,-1);
6382 SumTracks(seeds,arr);
6383 SignClusters(seeds,fnumber,fdensity);
6385 arr = Tracking(4,nup-3,nup-3-gap,cuts,-1);
6386 SumTracks(seeds,arr);
6387 SignClusters(seeds,fnumber,fdensity);
6391 for (Int_t delta = 3; delta<30; delta+=5){
6397 arr = Tracking(4,nup-1-delta,nup-1-delta-gap,cuts,-1);
6398 SumTracks(seeds,arr);
6399 SignClusters(seeds,fnumber,fdensity);
6401 arr = Tracking(4,nup-3-delta,nup-5-delta-gap,cuts,4);
6402 SumTracks(seeds,arr);
6403 SignClusters(seeds,fnumber,fdensity);
6414 Int_t fLastSeedRowSec=AliTPCReconstructor::GetRecoParam()->GetLastSeedRowSec();
6417 for (Int_t delta = 30; delta<fLastSeedRowSec; delta+=10){
6423 arr = Tracking(4,nup-1-delta,nup-1-delta-gap,cuts,-1);
6424 SumTracks(seeds,arr);
6425 SignClusters(seeds,fnumber,fdensity);
6427 arr = Tracking(4,nup-5-delta,nup-5-delta-gap,cuts,5 );
6428 SumTracks(seeds,arr);
6429 SignClusters(seeds,fnumber,fdensity);
6433 Info("Tracking()","\n\nSecondary seeding\t%d\n\n",seeds->GetEntriesFast());
6444 TObjArray * AliTPCtrackerMI::TrackingSpecial()
6447 // seeding adjusted for laser and cosmic tests - short tracks with big inclination angle
6448 // no primary vertex seeding tried
6452 Int_t nup=fOuterSec->GetNRows()+fInnerSec->GetNRows();
6454 TObjArray * seeds = new TObjArray;
6459 Float_t fnumber = 3.0;
6460 Float_t fdensity = 3.0;
6463 cuts[0] = AliTPCReconstructor::GetRecoParam()->GetMaxC(); // max curvature
6464 cuts[1] = 3.5; // max tan(phi) angle for seeding
6465 cuts[2] = 3.; // not used (cut on z primary vertex)
6466 cuts[3] = 3.5; // max tan(theta) angle for seeding
6468 for (Int_t delta = 0; nup-delta-gap-1>0; delta+=3){
6470 arr = Tracking(4,nup-1-delta,nup-1-delta-gap,cuts,-1);
6471 SumTracks(seeds,arr);
6472 SignClusters(seeds,fnumber,fdensity);
6476 Info("Tracking()","\n\nSecondary seeding\t%d\n\n",seeds->GetEntriesFast());
6487 void AliTPCtrackerMI::SumTracks(TObjArray *arr1,TObjArray *&arr2) const
6490 //sum tracks to common container
6491 //remove suspicious tracks
6492 Int_t nseed = arr2->GetEntriesFast();
6493 for (Int_t i=0;i<nseed;i++){
6494 AliTPCseed *pt=(AliTPCseed*)arr2->UncheckedAt(i);
6497 // remove tracks with too big curvature
6499 if (TMath::Abs(pt->GetC())>AliTPCReconstructor::GetRecoParam()->GetMaxC()){
6500 delete arr2->RemoveAt(i);
6503 // REMOVE VERY SHORT TRACKS
6504 if (pt->GetNumberOfClusters()<20){
6505 delete arr2->RemoveAt(i);
6508 // NORMAL ACTIVE TRACK
6509 if (pt->IsActive()){
6510 arr1->AddLast(arr2->RemoveAt(i));
6513 //remove not usable tracks
6514 if (pt->GetRemoval()!=10){
6515 delete arr2->RemoveAt(i);
6519 // ENABLE ONLY ENOUGH GOOD STOPPED TRACKS
6520 if (pt->GetDensityFirst(20)>0.8 || pt->GetDensityFirst(30)>0.8 || pt->GetDensityFirst(40)>0.7)
6521 arr1->AddLast(arr2->RemoveAt(i));
6523 delete arr2->RemoveAt(i);
6527 delete arr2; arr2 = 0;
6532 void AliTPCtrackerMI::ParallelTracking(TObjArray * arr, Int_t rfirst, Int_t rlast)
6535 // try to track in parralel
6537 Int_t nseed=arr->GetEntriesFast();
6538 //prepare seeds for tracking
6539 for (Int_t i=0; i<nseed; i++) {
6540 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i), &t=*pt;
6542 if (!t.IsActive()) continue;
6543 // follow prolongation to the first layer
6544 if ( (fSectors ==fInnerSec) || (t.GetFirstPoint()-fParam->GetNRowLow()>rfirst+1) )
6545 FollowProlongation(t, rfirst+1);
6550 for (Int_t nr=rfirst; nr>=rlast; nr--){
6551 if (nr<fInnerSec->GetNRows())
6552 fSectors = fInnerSec;
6554 fSectors = fOuterSec;
6555 // make indexes with the cluster tracks for given
6557 // find nearest cluster
6558 for (Int_t i=0; i<nseed; i++) {
6559 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i), &t=*pt;
6561 if (nr==80) pt->UpdateReference();
6562 if (!pt->IsActive()) continue;
6563 // if ( (fSectors ==fOuterSec) && (pt->fFirstPoint-fParam->GetNRowLow())<nr) continue;
6564 if (pt->GetRelativeSector()>17) {
6567 UpdateClusters(t,nr);
6569 // prolonagate to the nearest cluster - if founded
6570 for (Int_t i=0; i<nseed; i++) {
6571 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
6573 if (!pt->IsActive()) continue;
6574 // if ((fSectors ==fOuterSec) && (pt->fFirstPoint-fParam->GetNRowLow())<nr) continue;
6575 if (pt->GetRelativeSector()>17) {
6578 FollowToNextCluster(*pt,nr);
6583 void AliTPCtrackerMI::PrepareForBackProlongation(TObjArray * arr,Float_t fac) const
6587 // if we use TPC track itself we have to "update" covariance
6589 Int_t nseed= arr->GetEntriesFast();
6590 for (Int_t i=0;i<nseed;i++){
6591 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
6595 //rotate to current local system at first accepted point
6596 Int_t index = pt->GetClusterIndex2(pt->GetFirstPoint());
6597 Int_t sec = (index&0xff000000)>>24;
6599 Float_t angle1 = fInnerSec->GetAlpha()*sec+fInnerSec->GetAlphaShift();
6600 if (angle1>TMath::Pi())
6601 angle1-=2.*TMath::Pi();
6602 Float_t angle2 = pt->GetAlpha();
6604 if (TMath::Abs(angle1-angle2)>0.001){
6605 pt->Rotate(angle1-angle2);
6606 //angle2 = pt->GetAlpha();
6607 //pt->fRelativeSector = pt->GetAlpha()/fInnerSec->GetAlpha();
6608 //if (pt->GetAlpha()<0)
6609 // pt->fRelativeSector+=18;
6610 //sec = pt->fRelativeSector;
6619 void AliTPCtrackerMI::PrepareForProlongation(TObjArray * arr, Float_t fac) const
6623 // if we use TPC track itself we have to "update" covariance
6625 Int_t nseed= arr->GetEntriesFast();
6626 for (Int_t i=0;i<nseed;i++){
6627 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
6630 pt->SetFirstPoint(pt->GetLastPoint());
6638 Int_t AliTPCtrackerMI::PropagateBack(TObjArray * arr)
6641 // make back propagation
6643 Int_t nseed= arr->GetEntriesFast();
6644 for (Int_t i=0;i<nseed;i++){
6645 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
6646 if (pt&& pt->GetKinkIndex(0)<=0) {
6647 //AliTPCseed *pt2 = new AliTPCseed(*pt);
6648 fSectors = fInnerSec;
6649 //FollowBackProlongation(*pt,fInnerSec->GetNRows()-1);
6650 //fSectors = fOuterSec;
6651 FollowBackProlongation(*pt,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1);
6652 //if (pt->GetNumberOfClusters()<(pt->fEsd->GetTPCclusters(0)) ){
6653 // Error("PropagateBack","Not prolonged track %d",pt->GetLabel());
6654 // FollowBackProlongation(*pt2,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1);
6657 if (pt&& pt->GetKinkIndex(0)>0) {
6658 AliESDkink * kink = fEvent->GetKink(pt->GetKinkIndex(0)-1);
6659 pt->SetFirstPoint(kink->GetTPCRow0());
6660 fSectors = fInnerSec;
6661 FollowBackProlongation(*pt,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1);
6669 Int_t AliTPCtrackerMI::PropagateForward2(TObjArray * arr)
6672 // make forward propagation
6674 Int_t nseed= arr->GetEntriesFast();
6676 for (Int_t i=0;i<nseed;i++){
6677 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
6679 FollowProlongation(*pt,0);
6688 Int_t AliTPCtrackerMI::PropagateForward()
6691 // propagate track forward
6693 Int_t nseed = fSeeds->GetEntriesFast();
6694 for (Int_t i=0;i<nseed;i++){
6695 AliTPCseed *pt = (AliTPCseed*)fSeeds->UncheckedAt(i);
6697 AliTPCseed &t = *pt;
6698 Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
6699 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
6700 if (alpha < 0. ) alpha += 2.*TMath::Pi();
6701 t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN);
6705 fSectors = fOuterSec;
6706 ParallelTracking(fSeeds,fOuterSec->GetNRows()+fInnerSec->GetNRows()-1,fInnerSec->GetNRows());
6707 fSectors = fInnerSec;
6708 ParallelTracking(fSeeds,fInnerSec->GetNRows()-1,0);
6717 Int_t AliTPCtrackerMI::PropagateBack(AliTPCseed * pt, Int_t row0, Int_t row1)
6720 // make back propagation, in between row0 and row1
6724 fSectors = fInnerSec;
6727 if (row1<fSectors->GetNRows())
6730 r1 = fSectors->GetNRows()-1;
6732 if (row0<fSectors->GetNRows()&& r1>0 )
6733 FollowBackProlongation(*pt,r1);
6734 if (row1<=fSectors->GetNRows())
6737 r1 = row1 - fSectors->GetNRows();
6738 if (r1<=0) return 0;
6739 if (r1>=fOuterSec->GetNRows()) return 0;
6740 fSectors = fOuterSec;
6741 return FollowBackProlongation(*pt,r1);
6749 void AliTPCtrackerMI::GetShape(AliTPCseed * seed, Int_t row)
6753 AliTPCClusterParam * clparam = AliTPCcalibDB::Instance()->GetClusterParam();
6754 Float_t zdrift = TMath::Abs((fParam->GetZLength(0)-TMath::Abs(seed->GetZ())));
6755 Int_t type = (seed->GetSector() < fParam->GetNSector()/2) ? 0: (row>126) ? 1:2;
6756 Double_t angulary = seed->GetSnp();
6757 angulary = angulary*angulary/(1.-angulary*angulary);
6758 Double_t angularz = seed->GetTgl()*seed->GetTgl()*(1.+angulary);
6760 Double_t sigmay = clparam->GetRMS0(0,type,zdrift,TMath::Sqrt(TMath::Abs(angulary)));
6761 Double_t sigmaz = clparam->GetRMS0(1,type,zdrift,TMath::Sqrt(TMath::Abs(angularz)));
6762 seed->SetCurrentSigmaY2(sigmay*sigmay);
6763 seed->SetCurrentSigmaZ2(sigmaz*sigmaz);
6764 // Float_t sd2 = TMath::Abs((fParam->GetZLength(0)-TMath::Abs(seed->GetZ())))*fParam->GetDiffL()*fParam->GetDiffL();
6765 // // Float_t padlength = fParam->GetPadPitchLength(seed->fSector);
6766 // Float_t padlength = GetPadPitchLength(row);
6768 // Float_t sresy = (seed->GetSector() < fParam->GetNSector()/2) ? 0.2 :0.3;
6769 // seed->SetCurrentSigmaY2(sd2+padlength*padlength*angulary/12.+sresy*sresy);
6771 // Float_t sresz = fParam->GetZSigma();
6772 // seed->SetCurrentSigmaZ2(sd2+padlength*padlength*angularz*angularz*(1+angulary)/12.+sresz*sresz);
6774 Float_t wy = GetSigmaY(seed);
6775 Float_t wz = GetSigmaZ(seed);
6778 if (TMath::Abs(wy/seed->fCurrentSigmaY2-1)>0.0001 || TMath::Abs(wz/seed->fCurrentSigmaZ2-1)>0.0001 ){
6779 printf("problem\n");
6786 //__________________________________________________________________________
6787 void AliTPCtrackerMI::CookLabel(AliKalmanTrack *tk, Float_t wrong) const {
6788 //--------------------------------------------------------------------
6789 //This function "cooks" a track label. If label<0, this track is fake.
6790 //--------------------------------------------------------------------
6791 AliTPCseed * t = dynamic_cast<AliTPCseed*>(tk);
6793 printf("%s:%d wrong type \n",(char*)__FILE__,__LINE__);
6797 Int_t noc=t->GetNumberOfClusters();
6799 //printf("\nnot founded prolongation\n\n\n");
6805 AliTPCclusterMI *clusters[160];
6807 for (Int_t i=0;i<160;i++) {
6814 for (i=0; i<160 && current<noc; i++) {
6816 Int_t index=t->GetClusterIndex2(i);
6817 if (index<=0) continue;
6818 if (index&0x8000) continue;
6820 //clusters[current]=GetClusterMI(index);
6821 if (t->GetClusterPointer(i)){
6822 clusters[current]=t->GetClusterPointer(i);
6828 Int_t lab=123456789;
6829 for (i=0; i<noc; i++) {
6830 AliTPCclusterMI *c=clusters[i];
6832 lab=TMath::Abs(c->GetLabel(0));
6834 for (j=0; j<noc; j++) if (lb[j]==lab || mx[j]==0) break;
6840 for (i=0; i<noc; i++) if (mx[i]>max) {max=mx[i]; lab=lb[i];}
6842 for (i=0; i<noc; i++) {
6843 AliTPCclusterMI *c=clusters[i];
6845 if (TMath::Abs(c->GetLabel(1)) == lab ||
6846 TMath::Abs(c->GetLabel(2)) == lab ) max++;
6849 if ((1.- Float_t(max)/noc) > wrong) lab=-lab;
6852 Int_t tail=Int_t(0.10*noc);
6855 for (i=1; i<=160&&ind<tail; i++) {
6856 // AliTPCclusterMI *c=clusters[noc-i];
6857 AliTPCclusterMI *c=clusters[i];
6859 if (lab == TMath::Abs(c->GetLabel(0)) ||
6860 lab == TMath::Abs(c->GetLabel(1)) ||
6861 lab == TMath::Abs(c->GetLabel(2))) max++;
6864 if (max < Int_t(0.5*tail)) lab=-lab;
6871 //delete[] clusters;
6875 //__________________________________________________________________________
6876 Int_t AliTPCtrackerMI::CookLabel(AliTPCseed *t, Float_t wrong,Int_t first, Int_t last) const {
6877 //--------------------------------------------------------------------
6878 //This function "cooks" a track label. If label<0, this track is fake.
6879 //--------------------------------------------------------------------
6880 Int_t noc=t->GetNumberOfClusters();
6882 //printf("\nnot founded prolongation\n\n\n");
6888 AliTPCclusterMI *clusters[160];
6890 for (Int_t i=0;i<160;i++) {
6897 for (i=0; i<160 && current<noc; i++) {
6898 if (i<first) continue;
6899 if (i>last) continue;
6900 Int_t index=t->GetClusterIndex2(i);
6901 if (index<=0) continue;
6902 if (index&0x8000) continue;
6904 //clusters[current]=GetClusterMI(index);
6905 if (t->GetClusterPointer(i)){
6906 clusters[current]=t->GetClusterPointer(i);
6911 if (noc<5) return -1;
6912 Int_t lab=123456789;
6913 for (i=0; i<noc; i++) {
6914 AliTPCclusterMI *c=clusters[i];
6916 lab=TMath::Abs(c->GetLabel(0));
6918 for (j=0; j<noc; j++) if (lb[j]==lab || mx[j]==0) break;
6924 for (i=0; i<noc; i++) if (mx[i]>max) {max=mx[i]; lab=lb[i];}
6926 for (i=0; i<noc; i++) {
6927 AliTPCclusterMI *c=clusters[i];
6929 if (TMath::Abs(c->GetLabel(1)) == lab ||
6930 TMath::Abs(c->GetLabel(2)) == lab ) max++;
6933 if ((1.- Float_t(max)/noc) > wrong) lab=-lab;
6936 Int_t tail=Int_t(0.10*noc);
6939 for (i=1; i<=160&&ind<tail; i++) {
6940 // AliTPCclusterMI *c=clusters[noc-i];
6941 AliTPCclusterMI *c=clusters[i];
6943 if (lab == TMath::Abs(c->GetLabel(0)) ||
6944 lab == TMath::Abs(c->GetLabel(1)) ||
6945 lab == TMath::Abs(c->GetLabel(2))) max++;
6948 if (max < Int_t(0.5*tail)) lab=-lab;
6951 // t->SetLabel(lab);
6955 //delete[] clusters;
6959 Int_t AliTPCtrackerMI::GetRowNumber(Double_t x[3]) const
6961 //return pad row number for given x vector
6962 Float_t phi = TMath::ATan2(x[1],x[0]);
6963 if(phi<0) phi=2.*TMath::Pi()+phi;
6964 // Get the local angle in the sector philoc
6965 const Float_t kRaddeg = 180/3.14159265358979312;
6966 Float_t phiangle = (Int_t (phi*kRaddeg/20.) + 0.5)*20./kRaddeg;
6967 Double_t localx = x[0]*TMath::Cos(phiangle)-x[1]*TMath::Sin(phiangle);
6968 return GetRowNumber(localx);
6973 void AliTPCtrackerMI::MakeBitmaps(AliTPCseed *t)
6975 //-----------------------------------------------------------------------
6976 // Fill the cluster and sharing bitmaps of the track
6977 //-----------------------------------------------------------------------
6979 Int_t firstpoint = 0;
6980 Int_t lastpoint = 159;
6981 AliTPCTrackerPoint *point;
6983 for (int iter=firstpoint; iter<lastpoint; iter++) {
6984 point = t->GetTrackPoint(iter);
6986 t->SetClusterMapBit(iter, kTRUE);
6987 if (point->IsShared())
6988 t->SetSharedMapBit(iter,kTRUE);
6990 t->SetSharedMapBit(iter, kFALSE);
6993 t->SetClusterMapBit(iter, kFALSE);
6994 t->SetSharedMapBit(iter, kFALSE);
7001 void AliTPCtrackerMI::AddCovariance(AliTPCseed * seed){
7003 // Adding systematic error
7004 // !!!! the systematic error for element 4 is in 1/cm not in pt
7006 const Double_t *param = AliTPCReconstructor::GetRecoParam()->GetSystematicError();
7008 for (Int_t i=0;i<15;i++) covar[i]=0;
7014 covar[0] = param[0]*param[0];
7015 covar[2] = param[1]*param[1];
7016 covar[5] = param[2]*param[2];
7017 covar[9] = param[3]*param[3];
7018 Double_t facC = AliTracker::GetBz()*kB2C;
7019 covar[14]= param[4]*param[4]*facC*facC;
7020 seed->AddCovariance(covar);