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 if( AliTPCReconstructor::StreamLevel()>5){
4424 cstream<<"Splitted"<<
4428 "Tr0.="<<track0<< // seed0
4429 "Tr1.="<<track1<< // seed1
4430 "h0.="<<&helixes[i0]<<
4431 "h1.="<<&helixes[i1]<<
4433 "sum="<<sum<< //the sum of rows with cl in both
4434 "sum0="<<sum0<< //the sum of rows with cl in 0 track
4435 "sum1="<<sum1<< //the sum of rows with cl in 1 track
4436 "sums="<<sums<< //the sum of shared clusters
4437 "xm0="<<xm[i0]<< // the center of track
4438 "xm1="<<xm[i1]<< // the x center of track
4439 // General cut variables
4440 "dfi="<<dfi<< // distance in fi angle
4441 "dtheta="<<dtheta<< // distance int theta angle
4444 "dist0="<<dist[4][0]<< //distance x
4445 "dist1="<<dist[4][1]<< //distance y
4446 "dist2="<<dist[4][2]<< //distance z
4447 "mdist0="<<mdist[0]<< //distance x
4448 "mdist1="<<mdist[1]<< //distance y
4449 "mdist2="<<mdist[2]<< //distance z
4452 delete array->RemoveAt(i1);
4459 AliInfo("Time for splitted tracks removal");
4465 void AliTPCtrackerMI::FindCurling(TObjArray * array, AliESDEvent */*esd*/, Int_t iter)
4468 // find Curling tracks
4469 // Use AliTPCReconstructor::StreamLevel()>1 if you want to tune parameters - cuts
4472 // Algorithm done in 2 phases - because of CPU consumption
4473 // it is n^2 algorithm - for lead-lead 20000x20000 combination are investigated
4474 // see detal in MC part what can be used to cut
4478 const Float_t kMaxC = 400; // maximal curvature to of the track
4479 const Float_t kMaxdTheta = 0.15; // maximal distance in theta
4480 const Float_t kMaxdPhi = 0.15; // maximal distance in phi
4481 const Float_t kPtRatio = 0.3; // ratio between pt
4482 const Float_t kMinDCAR = 2.; // distance to the primary vertex in r - see cpipe cut
4485 // Curling tracks cuts
4488 const Float_t kMaxDeltaRMax = 40; // distance in outer radius
4489 const Float_t kMaxDeltaRMin = 5.; // distance in lower radius - see cpipe cut
4490 const Float_t kMinAngle = 2.9; // angle between tracks
4491 const Float_t kMaxDist = 5; // biggest distance
4493 // The cuts can be tuned using the "MC information stored in Multi tree ==> see FindMultiMC
4496 TCut csign("csign","Tr0.fP[4]*Tr1.fP[4]<0"); //opposite sign
4497 TCut cmax("cmax","abs(Tr0.GetC())>1/400");
4498 TCut cda("cda","sqrt(dtheta^2+dfi^2)<0.15");
4499 TCut ccratio("ccratio","abs((Tr0.fP[4]+Tr1.fP[4])/(abs(Tr0.fP[4])+abs(Tr1.fP[4])))<0.3");
4500 TCut cpipe("cpipe", "min(abs(r0-rc0),abs(r1-rc1))>5");
4502 TCut cdrmax("cdrmax","abs(abs(rc0+r0)-abs(rc1+r1))<40")
4503 TCut cdrmin("cdrmin","abs(abs(rc0+r0)-abs(rc1+r1))<10")
4505 Multi->Draw("dfi","iter==0"+csign+cmax+cda+ccratio); ~94% of curling tracks fulfill
4506 Multi->Draw("min(abs(r0-rc0),abs(r1-rc1))","iter==0&&abs(lab1)==abs(lab0)"+csign+cmax+cda+ccratio+cpipe+cdrmin+cdrmax); //80%
4508 Curling2->Draw("dfi","iter==0&&abs(lab0)==abs(lab1)"+csign+cmax+cdtheta+cdfi+ccratio)
4514 Int_t nentries = array->GetEntriesFast();
4515 AliHelix *helixes = new AliHelix[nentries];
4516 for (Int_t i=0;i<nentries;i++){
4517 AliTPCseed* track = (AliTPCseed*)array->At(i);
4518 if (!track) continue;
4519 track->SetCircular(0);
4520 new (&helixes[i]) AliHelix(*track);
4526 Double_t phase[2][2],radius[2];
4530 TTreeSRedirector &cstream = *fDebugStreamer;
4532 for (Int_t i0=0;i0<nentries;i0++){
4533 AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
4534 if (!track0) continue;
4535 if (TMath::Abs(track0->GetC())<1/kMaxC) continue;
4536 Float_t xc0 = helixes[i0].GetHelix(6);
4537 Float_t yc0 = helixes[i0].GetHelix(7);
4538 Float_t r0 = helixes[i0].GetHelix(8);
4539 Float_t rc0 = TMath::Sqrt(xc0*xc0+yc0*yc0);
4540 Float_t fi0 = TMath::ATan2(yc0,xc0);
4542 for (Int_t i1=i0+1;i1<nentries;i1++){
4543 AliTPCseed * track1 = (AliTPCseed*)array->At(i1);
4544 if (!track1) continue;
4545 if (TMath::Abs(track1->GetC())<1/kMaxC) continue;
4546 Float_t xc1 = helixes[i1].GetHelix(6);
4547 Float_t yc1 = helixes[i1].GetHelix(7);
4548 Float_t r1 = helixes[i1].GetHelix(8);
4549 Float_t rc1 = TMath::Sqrt(xc1*xc1+yc1*yc1);
4550 Float_t fi1 = TMath::ATan2(yc1,xc1);
4552 Float_t dfi = fi0-fi1;
4555 if (dfi>1.5*TMath::Pi()) dfi-=TMath::Pi(); // take care about edge effect
4556 if (dfi<-1.5*TMath::Pi()) dfi+=TMath::Pi(); //
4557 Float_t dtheta = TMath::Abs(track0->GetTgl()-track1->GetTgl())<TMath::Abs(track0->GetTgl()+track1->GetTgl())? track0->GetTgl()-track1->GetTgl():track0->GetTgl()+track1->GetTgl();
4561 if (track0->GetBConstrain()&&track1->GetBConstrain()) continue; // not constrained
4562 if (track1->GetSigned1Pt()*track0->GetSigned1Pt()>0) continue; // not the same sign
4563 if ( TMath::Abs(track1->GetTgl()+track0->GetTgl())>kMaxdTheta) continue; //distance in the Theta
4564 if ( TMath::Abs(dfi)>kMaxdPhi) continue; //distance in phi
4565 if ( TMath::Sqrt(dfi*dfi+dtheta*dtheta)>kMaxdPhi) continue; //common angular offset
4567 Float_t pt0 = track0->GetSignedPt();
4568 Float_t pt1 = track1->GetSignedPt();
4569 if ((TMath::Abs(pt0+pt1)/(TMath::Abs(pt0)+TMath::Abs(pt1)))>kPtRatio) continue;
4570 if ((iter==1) && TMath::Abs(TMath::Abs(rc0+r0)-TMath::Abs(rc1+r1))>kMaxDeltaRMax) continue;
4571 if ((iter!=1) &&TMath::Abs(TMath::Abs(rc0-r0)-TMath::Abs(rc1-r1))>kMaxDeltaRMin) continue;
4572 if (TMath::Min(TMath::Abs(rc0-r0),TMath::Abs(rc1-r1))<kMinDCAR) continue;
4575 // Now find closest approach
4579 Int_t npoints = helixes[i0].GetRPHIintersections(helixes[i1], phase, radius,10);
4580 if (npoints==0) continue;
4581 helixes[i0].GetClosestPhases(helixes[i1], phase);
4585 Double_t hangles[3];
4586 helixes[i0].Evaluate(phase[0][0],xyz0);
4587 helixes[i1].Evaluate(phase[0][1],xyz1);
4589 helixes[i0].GetAngle(phase[0][0],helixes[i1],phase[0][1],hangles);
4590 Double_t deltah[2],deltabest;
4591 if (TMath::Abs(hangles[2])<kMinAngle) continue;
4595 helixes[i0].ParabolicDCA(helixes[i1],phase[0][0],phase[0][1],radius[0],deltah[0],2);
4597 helixes[i0].ParabolicDCA(helixes[i1],phase[1][0],phase[1][1],radius[1],deltah[1],2);
4598 if (deltah[1]<deltah[0]) ibest=1;
4600 deltabest = TMath::Sqrt(deltah[ibest]);
4601 helixes[i0].Evaluate(phase[ibest][0],xyz0);
4602 helixes[i1].Evaluate(phase[ibest][1],xyz1);
4603 helixes[i0].GetAngle(phase[ibest][0],helixes[i1],phase[ibest][1],hangles);
4604 Double_t radiusbest = TMath::Sqrt(radius[ibest]);
4606 if (deltabest>kMaxDist) continue;
4607 // if (mindcar+mindcaz<40 && (TMath::Abs(hangles[2])<kMinAngle ||deltabest>3)) continue;
4608 Bool_t sign =kFALSE;
4609 if (hangles[2]>kMinAngle) sign =kTRUE;
4612 // circular[i0] = kTRUE;
4613 // circular[i1] = kTRUE;
4614 if (track0->OneOverPt()<track1->OneOverPt()){
4615 track0->SetCircular(track0->GetCircular()+1);
4616 track1->SetCircular(track1->GetCircular()+2);
4619 track1->SetCircular(track1->GetCircular()+1);
4620 track0->SetCircular(track0->GetCircular()+2);
4623 if (AliTPCReconstructor::StreamLevel()>1){
4625 //debug stream to tune "fine" cuts
4626 Int_t lab0=track0->GetLabel();
4627 Int_t lab1=track1->GetLabel();
4628 cstream<<"Curling2"<<
4644 "npoints="<<npoints<<
4645 "hangles0="<<hangles[0]<<
4646 "hangles1="<<hangles[1]<<
4647 "hangles2="<<hangles[2]<<
4650 "radius="<<radiusbest<<
4651 "deltabest="<<deltabest<<
4652 "phase0="<<phase[ibest][0]<<
4653 "phase1="<<phase[ibest][1]<<
4661 if (AliTPCReconstructor::StreamLevel()>1) {
4662 AliInfo("Time for curling tracks removal");
4671 void AliTPCtrackerMI::FindKinks(TObjArray * array, AliESDEvent *esd)
4678 TObjArray *kinks= new TObjArray(10000);
4679 // TObjArray *v0s= new TObjArray(10000);
4680 Int_t nentries = array->GetEntriesFast();
4681 AliHelix *helixes = new AliHelix[nentries];
4682 Int_t *sign = new Int_t[nentries];
4683 Int_t *nclusters = new Int_t[nentries];
4684 Float_t *alpha = new Float_t[nentries];
4685 AliKink *kink = new AliKink();
4686 Int_t * usage = new Int_t[nentries];
4687 Float_t *zm = new Float_t[nentries];
4688 Float_t *z0 = new Float_t[nentries];
4689 Float_t *fim = new Float_t[nentries];
4690 Float_t *shared = new Float_t[nentries];
4691 Bool_t *circular = new Bool_t[nentries];
4692 Float_t *dca = new Float_t[nentries];
4693 //const AliESDVertex * primvertex = esd->GetVertex();
4695 // nentries = array->GetEntriesFast();
4700 for (Int_t i=0;i<nentries;i++){
4703 AliTPCseed* track = (AliTPCseed*)array->At(i);
4704 if (!track) continue;
4705 track->SetCircular(0);
4707 track->UpdatePoints();
4708 if (( track->GetPoints()[2]- track->GetPoints()[0])>5 && track->GetPoints()[3]>0.8){
4710 nclusters[i]=track->GetNumberOfClusters();
4711 alpha[i] = track->GetAlpha();
4712 new (&helixes[i]) AliHelix(*track);
4714 helixes[i].Evaluate(0,xyz);
4715 sign[i] = (track->GetC()>0) ? -1:1;
4718 if (track->GetProlongation(x,y,z)){
4720 fim[i] = alpha[i]+TMath::ATan2(y,x);
4723 zm[i] = track->GetZ();
4727 circular[i]= kFALSE;
4728 if (track->GetProlongation(0,y,z)) z0[i] = z;
4729 dca[i] = track->GetD(0,0);
4735 Int_t ncandidates =0;
4738 Double_t phase[2][2],radius[2];
4741 // Find circling track
4742 TTreeSRedirector &cstream = *fDebugStreamer;
4744 for (Int_t i0=0;i0<nentries;i0++){
4745 AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
4746 if (!track0) continue;
4747 if (track0->GetNumberOfClusters()<40) continue;
4748 if (TMath::Abs(1./track0->GetC())>200) continue;
4749 for (Int_t i1=i0+1;i1<nentries;i1++){
4750 AliTPCseed * track1 = (AliTPCseed*)array->At(i1);
4751 if (!track1) continue;
4752 if (track1->GetNumberOfClusters()<40) continue;
4753 if ( TMath::Abs(track1->GetTgl()+track0->GetTgl())>0.1) continue;
4754 if (track0->GetBConstrain()&&track1->GetBConstrain()) continue;
4755 if (TMath::Abs(1./track1->GetC())>200) continue;
4756 if (track1->GetSigned1Pt()*track0->GetSigned1Pt()>0) continue;
4757 if (track1->GetTgl()*track0->GetTgl()>0) continue;
4758 if (TMath::Max(TMath::Abs(1./track0->GetC()),TMath::Abs(1./track1->GetC()))>190) continue;
4759 if (track0->GetBConstrain()&&track1->OneOverPt()<track0->OneOverPt()) continue; //returning - lower momenta
4760 if (track1->GetBConstrain()&&track0->OneOverPt()<track1->OneOverPt()) continue; //returning - lower momenta
4762 Float_t mindcar = TMath::Min(TMath::Abs(dca[i0]),TMath::Abs(dca[i1]));
4763 if (mindcar<5) continue;
4764 Float_t mindcaz = TMath::Min(TMath::Abs(z0[i0]-GetZ()),TMath::Abs(z0[i1]-GetZ()));
4765 if (mindcaz<5) continue;
4766 if (mindcar+mindcaz<20) continue;
4769 Float_t xc0 = helixes[i0].GetHelix(6);
4770 Float_t yc0 = helixes[i0].GetHelix(7);
4771 Float_t r0 = helixes[i0].GetHelix(8);
4772 Float_t xc1 = helixes[i1].GetHelix(6);
4773 Float_t yc1 = helixes[i1].GetHelix(7);
4774 Float_t r1 = helixes[i1].GetHelix(8);
4776 Float_t rmean = (r0+r1)*0.5;
4777 Float_t delta =TMath::Sqrt((xc1-xc0)*(xc1-xc0)+(yc1-yc0)*(yc1-yc0));
4778 //if (delta>30) continue;
4779 if (delta>rmean*0.25) continue;
4780 if (TMath::Abs(r0-r1)/rmean>0.3) continue;
4782 Int_t npoints = helixes[i0].GetRPHIintersections(helixes[i1], phase, radius,10);
4783 if (npoints==0) continue;
4784 helixes[i0].GetClosestPhases(helixes[i1], phase);
4788 Double_t hangles[3];
4789 helixes[i0].Evaluate(phase[0][0],xyz0);
4790 helixes[i1].Evaluate(phase[0][1],xyz1);
4792 helixes[i0].GetAngle(phase[0][0],helixes[i1],phase[0][1],hangles);
4793 Double_t deltah[2],deltabest;
4794 if (hangles[2]<2.8) continue;
4797 helixes[i0].ParabolicDCA(helixes[i1],phase[0][0],phase[0][1],radius[0],deltah[0],2);
4799 helixes[i0].ParabolicDCA(helixes[i1],phase[1][0],phase[1][1],radius[1],deltah[1],2);
4800 if (deltah[1]<deltah[0]) ibest=1;
4802 deltabest = TMath::Sqrt(deltah[ibest]);
4803 helixes[i0].Evaluate(phase[ibest][0],xyz0);
4804 helixes[i1].Evaluate(phase[ibest][1],xyz1);
4805 helixes[i0].GetAngle(phase[ibest][0],helixes[i1],phase[ibest][1],hangles);
4806 Double_t radiusbest = TMath::Sqrt(radius[ibest]);
4808 if (deltabest>6) continue;
4809 if (mindcar+mindcaz<40 && (hangles[2]<3.12||deltabest>3)) continue;
4810 Bool_t sign =kFALSE;
4811 if (hangles[2]>3.06) sign =kTRUE;
4814 circular[i0] = kTRUE;
4815 circular[i1] = kTRUE;
4816 if (track0->OneOverPt()<track1->OneOverPt()){
4817 track0->SetCircular(track0->GetCircular()+1);
4818 track1->SetCircular(track1->GetCircular()+2);
4821 track1->SetCircular(track1->GetCircular()+1);
4822 track0->SetCircular(track0->GetCircular()+2);
4825 if (sign&&AliTPCReconstructor::StreamLevel()>1){
4827 Int_t lab0=track0->GetLabel();
4828 Int_t lab1=track1->GetLabel();
4829 cstream<<"Curling"<<
4836 "mindcar="<<mindcar<<
4837 "mindcaz="<<mindcaz<<
4840 "npoints="<<npoints<<
4841 "hangles0="<<hangles[0]<<
4842 "hangles2="<<hangles[2]<<
4847 "radius="<<radiusbest<<
4848 "deltabest="<<deltabest<<
4849 "phase0="<<phase[ibest][0]<<
4850 "phase1="<<phase[ibest][1]<<
4860 for (Int_t i =0;i<nentries;i++){
4861 if (sign[i]==0) continue;
4862 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
4869 Double_t cradius0 = 40*40;
4870 Double_t cradius1 = 270*270;
4873 Double_t cdist3=0.55;
4874 for (Int_t j =i+1;j<nentries;j++){
4876 if (sign[j]*sign[i]<1) continue;
4877 if ( (nclusters[i]+nclusters[j])>200) continue;
4878 if ( (nclusters[i]+nclusters[j])<80) continue;
4879 if ( TMath::Abs(zm[i]-zm[j])>60.) continue;
4880 if ( TMath::Abs(fim[i]-fim[j])>0.6 && TMath::Abs(fim[i]-fim[j])<5.7 ) continue;
4881 //AliTPCseed * track1 = (AliTPCseed*)array->At(j); Double_t phase[2][2],radius[2];
4882 Int_t npoints = helixes[i].GetRPHIintersections(helixes[j], phase, radius,20);
4883 if (npoints<1) continue;
4886 if (radius[0]<cradius0||radius[0]>cradius1) continue;
4889 if ( (radius[0]<cradius0||radius[0]>cradius1) && (radius[1]<cradius0||radius[1]>cradius1) ) continue;
4892 Double_t delta1=10000,delta2=10000;
4893 // cuts on the intersection radius
4894 helixes[i].LinearDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta1);
4895 if (radius[0]<20&&delta1<1) continue; //intersection at vertex
4896 if (radius[0]<10&&delta1<3) continue; //intersection at vertex
4898 helixes[i].LinearDCA(helixes[j],phase[1][0],phase[1][1],radius[1],delta2);
4899 if (radius[1]<20&&delta2<1) continue; //intersection at vertex
4900 if (radius[1]<10&&delta2<3) continue; //intersection at vertex
4903 Double_t distance1 = TMath::Min(delta1,delta2);
4904 if (distance1>cdist1) continue; // cut on DCA linear approximation
4906 npoints = helixes[i].GetRPHIintersections(helixes[j], phase, radius,20);
4907 helixes[i].ParabolicDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta1);
4908 if (radius[0]<20&&delta1<1) continue; //intersection at vertex
4909 if (radius[0]<10&&delta1<3) continue; //intersection at vertex
4912 helixes[i].ParabolicDCA(helixes[j],phase[1][0],phase[1][1],radius[1],delta2);
4913 if (radius[1]<20&&delta2<1) continue; //intersection at vertex
4914 if (radius[1]<10&&delta2<3) continue; //intersection at vertex
4916 distance1 = TMath::Min(delta1,delta2);
4919 rkink = TMath::Sqrt(radius[0]);
4922 rkink = TMath::Sqrt(radius[1]);
4924 if (distance1>cdist2) continue;
4927 AliTPCseed * track1 = (AliTPCseed*)array->At(j);
4930 Int_t row0 = GetRowNumber(rkink);
4931 if (row0<10) continue;
4932 if (row0>150) continue;
4935 Float_t dens00=-1,dens01=-1;
4936 Float_t dens10=-1,dens11=-1;
4938 Int_t found,foundable,shared;
4939 track0->GetClusterStatistic(0,row0-5, found, foundable,shared,kFALSE);
4940 if (foundable>5) dens00 = Float_t(found)/Float_t(foundable);
4941 track0->GetClusterStatistic(row0+5,155, found, foundable,shared,kFALSE);
4942 if (foundable>5) dens01 = Float_t(found)/Float_t(foundable);
4944 track1->GetClusterStatistic(0,row0-5, found, foundable,shared,kFALSE);
4945 if (foundable>10) dens10 = Float_t(found)/Float_t(foundable);
4946 track1->GetClusterStatistic(row0+5,155, found, foundable,shared,kFALSE);
4947 if (foundable>10) dens11 = Float_t(found)/Float_t(foundable);
4949 if (dens00<dens10 && dens01<dens11) continue;
4950 if (dens00>dens10 && dens01>dens11) continue;
4951 if (TMath::Max(dens00,dens10)<0.1) continue;
4952 if (TMath::Max(dens01,dens11)<0.3) continue;
4954 if (TMath::Min(dens00,dens10)>0.6) continue;
4955 if (TMath::Min(dens01,dens11)>0.6) continue;
4958 AliTPCseed * ktrack0, *ktrack1;
4967 if (TMath::Abs(ktrack0->GetC())>5) continue; // cut on the curvature for mother particle
4968 AliExternalTrackParam paramm(*ktrack0);
4969 AliExternalTrackParam paramd(*ktrack1);
4970 if (row0>60&&ktrack1->GetReference().GetX()>90.)new (¶md) AliExternalTrackParam(ktrack1->GetReference());
4973 kink->SetMother(paramm);
4974 kink->SetDaughter(paramd);
4977 Float_t x[3] = { kink->GetPosition()[0],kink->GetPosition()[1],kink->GetPosition()[2]};
4979 fParam->Transform0to1(x,index);
4980 fParam->Transform1to2(x,index);
4981 row0 = GetRowNumber(x[0]);
4983 if (kink->GetR()<100) continue;
4984 if (kink->GetR()>240) continue;
4985 if (kink->GetPosition()[2]/kink->GetR()>AliTPCReconstructor::GetCtgRange()) continue; //out of fiducial volume
4986 if (kink->GetDistance()>cdist3) continue;
4987 Float_t dird = kink->GetDaughterP()[0]*kink->GetPosition()[0]+kink->GetDaughterP()[1]*kink->GetPosition()[1]; // rough direction estimate
4988 if (dird<0) continue;
4990 Float_t dirm = kink->GetMotherP()[0]*kink->GetPosition()[0]+kink->GetMotherP()[1]*kink->GetPosition()[1]; // rough direction estimate
4991 if (dirm<0) continue;
4992 Float_t mpt = TMath::Sqrt(kink->GetMotherP()[0]*kink->GetMotherP()[0]+kink->GetMotherP()[1]*kink->GetMotherP()[1]);
4993 if (mpt<0.2) continue;
4996 //for high momenta momentum not defined well in first iteration
4997 Double_t qt = TMath::Sin(kink->GetAngle(2))*ktrack1->GetP();
4998 if (qt>0.35) continue;
5001 kink->SetLabel(CookLabel(ktrack0,0.4,0,row0),0);
5002 kink->SetLabel(CookLabel(ktrack1,0.4,row0,160),1);
5004 kink->SetTPCDensity(dens00,0,0);
5005 kink->SetTPCDensity(dens01,0,1);
5006 kink->SetTPCDensity(dens10,1,0);
5007 kink->SetTPCDensity(dens11,1,1);
5008 kink->SetIndex(i,0);
5009 kink->SetIndex(j,1);
5012 kink->SetTPCDensity(dens10,0,0);
5013 kink->SetTPCDensity(dens11,0,1);
5014 kink->SetTPCDensity(dens00,1,0);
5015 kink->SetTPCDensity(dens01,1,1);
5016 kink->SetIndex(j,0);
5017 kink->SetIndex(i,1);
5020 if (mpt<1||kink->GetAngle(2)>0.1){
5021 // angle and densities not defined yet
5022 if (kink->GetTPCDensityFactor()<0.8) continue;
5023 if ((2-kink->GetTPCDensityFactor())*kink->GetDistance() >0.25) continue;
5024 if (kink->GetAngle(2)*ktrack0->GetP()<0.003) continue; //too small angle
5025 if (kink->GetAngle(2)>0.2&&kink->GetTPCDensityFactor()<1.15) continue;
5026 if (kink->GetAngle(2)>0.2&&kink->GetTPCDensity(0,1)>0.05) continue;
5028 Float_t criticalangle = track0->GetSigmaSnp2()+track0->GetSigmaTgl2();
5029 criticalangle+= track1->GetSigmaSnp2()+track1->GetSigmaTgl2();
5030 criticalangle= 3*TMath::Sqrt(criticalangle);
5031 if (criticalangle>0.02) criticalangle=0.02;
5032 if (kink->GetAngle(2)<criticalangle) continue;
5035 Int_t drow = Int_t(2.+0.5/(0.05+kink->GetAngle(2))); // overlap region defined
5036 Float_t shapesum =0;
5038 for ( Int_t row = row0-drow; row<row0+drow;row++){
5039 if (row<0) continue;
5040 if (row>155) continue;
5041 if (ktrack0->GetClusterPointer(row)){
5042 AliTPCTrackerPoint *point =ktrack0->GetTrackPoint(row);
5043 shapesum+=point->GetSigmaY()+point->GetSigmaZ();
5046 if (ktrack1->GetClusterPointer(row)){
5047 AliTPCTrackerPoint *point =ktrack1->GetTrackPoint(row);
5048 shapesum+=point->GetSigmaY()+point->GetSigmaZ();
5053 kink->SetShapeFactor(-1.);
5056 kink->SetShapeFactor(shapesum/sum);
5058 // esd->AddKink(kink);
5059 kinks->AddLast(kink);
5065 // sort the kinks according quality - and refit them towards vertex
5067 Int_t nkinks = kinks->GetEntriesFast();
5068 Float_t *quality = new Float_t[nkinks];
5069 Int_t *indexes = new Int_t[nkinks];
5070 AliTPCseed *mothers = new AliTPCseed[nkinks];
5071 AliTPCseed *daughters = new AliTPCseed[nkinks];
5074 for (Int_t i=0;i<nkinks;i++){
5076 AliKink *kink = (AliKink*)kinks->At(i);
5078 // refit kinks towards vertex
5080 Int_t index0 = kink->GetIndex(0);
5081 Int_t index1 = kink->GetIndex(1);
5082 AliTPCseed * ktrack0 = (AliTPCseed*)array->At(index0);
5083 AliTPCseed * ktrack1 = (AliTPCseed*)array->At(index1);
5085 Int_t sumn=ktrack0->GetNumberOfClusters()+ktrack1->GetNumberOfClusters();
5087 // Refit Kink under if too small angle
5089 if (kink->GetAngle(2)<0.05){
5090 kink->SetTPCRow0(GetRowNumber(kink->GetR()));
5091 Int_t row0 = kink->GetTPCRow0();
5092 Int_t drow = Int_t(2.+0.5/(0.05+kink->GetAngle(2)));
5095 Int_t last = row0-drow;
5096 if (last<40) last=40;
5097 if (last<ktrack0->GetFirstPoint()+25) last = ktrack0->GetFirstPoint()+25;
5098 AliTPCseed* seed0 = ReSeed(ktrack0,last,kFALSE);
5101 Int_t first = row0+drow;
5102 if (first>130) first=130;
5103 if (first>ktrack1->GetLastPoint()-25) first = TMath::Max(ktrack1->GetLastPoint()-25,30);
5104 AliTPCseed* seed1 = ReSeed(ktrack1,first,kTRUE);
5106 if (seed0 && seed1){
5107 kink->SetStatus(1,8);
5108 if (RefitKink(*seed0,*seed1,*kink)) kink->SetStatus(1,9);
5109 row0 = GetRowNumber(kink->GetR());
5110 sumn = seed0->GetNumberOfClusters()+seed1->GetNumberOfClusters();
5111 mothers[i] = *seed0;
5112 daughters[i] = *seed1;
5115 delete kinks->RemoveAt(i);
5116 if (seed0) delete seed0;
5117 if (seed1) delete seed1;
5120 if (kink->GetDistance()>0.5 || kink->GetR()<110 || kink->GetR()>240) {
5121 delete kinks->RemoveAt(i);
5122 if (seed0) delete seed0;
5123 if (seed1) delete seed1;
5131 if (kink) quality[i] = 160*((0.1+kink->GetDistance())*(2.-kink->GetTPCDensityFactor()))/(sumn+40.); //the longest -clossest will win
5133 TMath::Sort(nkinks,quality,indexes,kFALSE);
5135 //remove double find kinks
5137 for (Int_t ikink0=1;ikink0<nkinks;ikink0++){
5138 AliKink * kink0 = (AliKink*) kinks->At(indexes[ikink0]);
5139 if (!kink0) continue;
5141 for (Int_t ikink1=0;ikink1<ikink0;ikink1++){
5142 if (!kink0) continue;
5143 AliKink * kink1 = (AliKink*) kinks->At(indexes[ikink1]);
5144 if (!kink1) continue;
5145 // if not close kink continue
5146 if (TMath::Abs(kink1->GetPosition()[2]-kink0->GetPosition()[2])>10) continue;
5147 if (TMath::Abs(kink1->GetPosition()[1]-kink0->GetPosition()[1])>10) continue;
5148 if (TMath::Abs(kink1->GetPosition()[0]-kink0->GetPosition()[0])>10) continue;
5150 AliTPCseed &mother0 = mothers[indexes[ikink0]];
5151 AliTPCseed &daughter0 = daughters[indexes[ikink0]];
5152 AliTPCseed &mother1 = mothers[indexes[ikink1]];
5153 AliTPCseed &daughter1 = daughters[indexes[ikink1]];
5154 Int_t row0 = (kink0->GetTPCRow0()+kink1->GetTPCRow0())/2;
5163 for (Int_t i=0;i<row0;i++){
5164 if (mother0.GetClusterIndex(i)>0 && mother1.GetClusterIndex(i)>0){
5167 if (mother0.GetClusterIndex(i)==mother1.GetClusterIndex(i)){
5174 for (Int_t i=row0;i<158;i++){
5175 if (daughter0.GetClusterIndex(i)>0 && daughter0.GetClusterIndex(i)>0){
5178 if (mother0.GetClusterIndex(i)==mother1.GetClusterIndex(i)){
5184 Float_t ratio = Float_t(same+1)/Float_t(both+1);
5185 Float_t ratiom = Float_t(samem+1)/Float_t(bothm+1);
5186 Float_t ratiod = Float_t(samed+1)/Float_t(bothd+1);
5187 if (ratio>0.3 && ratiom>0.5 &&ratiod>0.5) {
5188 Int_t sum0 = mother0.GetNumberOfClusters()+daughter0.GetNumberOfClusters();
5189 Int_t sum1 = mother1.GetNumberOfClusters()+daughter1.GetNumberOfClusters();
5191 shared[kink0->GetIndex(0)]= kTRUE;
5192 shared[kink0->GetIndex(1)]= kTRUE;
5193 delete kinks->RemoveAt(indexes[ikink0]);
5196 shared[kink1->GetIndex(0)]= kTRUE;
5197 shared[kink1->GetIndex(1)]= kTRUE;
5198 delete kinks->RemoveAt(indexes[ikink1]);
5205 for (Int_t i=0;i<nkinks;i++){
5206 AliKink * kink = (AliKink*) kinks->At(indexes[i]);
5207 if (!kink) continue;
5208 kink->SetTPCRow0(GetRowNumber(kink->GetR()));
5209 Int_t index0 = kink->GetIndex(0);
5210 Int_t index1 = kink->GetIndex(1);
5211 if (circular[index0]||circular[index1]&&kink->GetDistance()>0.2) continue;
5212 kink->SetMultiple(usage[index0],0);
5213 kink->SetMultiple(usage[index1],1);
5214 if (kink->GetMultiple()[0]+kink->GetMultiple()[1]>2) continue;
5215 if (kink->GetMultiple()[0]+kink->GetMultiple()[1]>0 && quality[indexes[i]]>0.2) continue;
5216 if (kink->GetMultiple()[0]+kink->GetMultiple()[1]>0 && kink->GetDistance()>0.2) continue;
5217 if (circular[index0]||circular[index1]&&kink->GetDistance()>0.1) continue;
5219 AliTPCseed * ktrack0 = (AliTPCseed*)array->At(index0);
5220 AliTPCseed * ktrack1 = (AliTPCseed*)array->At(index1);
5221 if (!ktrack0 || !ktrack1) continue;
5222 Int_t index = esd->AddKink(kink);
5225 if ( ktrack0->GetKinkIndex(0)==0 && ktrack1->GetKinkIndex(0)==0) { //best kink
5226 if (mothers[indexes[i]].GetNumberOfClusters()>20 && daughters[indexes[i]].GetNumberOfClusters()>20 && (mothers[indexes[i]].GetNumberOfClusters()+daughters[indexes[i]].GetNumberOfClusters())>100){
5227 *ktrack0 = mothers[indexes[i]];
5228 *ktrack1 = daughters[indexes[i]];
5232 ktrack0->SetKinkIndex(usage[index0],-(index+1));
5233 ktrack1->SetKinkIndex(usage[index1], (index+1));
5238 // Remove tracks corresponding to shared kink's
5240 for (Int_t i=0;i<nentries;i++){
5241 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5242 if (!track0) continue;
5243 if (track0->GetKinkIndex(0)!=0) continue;
5244 if (shared[i]) delete array->RemoveAt(i);
5249 RemoveUsed2(array,0.5,0.4,30);
5251 for (Int_t i=0;i<nentries;i++){
5252 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5253 if (!track0) continue;
5254 track0->CookdEdx(0.02,0.6);
5258 for (Int_t i=0;i<nentries;i++){
5259 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5260 if (!track0) continue;
5261 if (track0->Pt()<1.4) continue;
5262 //remove double high momenta tracks - overlapped with kink candidates
5265 for (Int_t icl=track0->GetFirstPoint();icl<track0->GetLastPoint(); icl++){
5266 if (track0->GetClusterPointer(icl)!=0){
5268 if (track0->GetClusterPointer(icl)->IsUsed(10)) shared++;
5271 if (Float_t(shared+1)/Float_t(all+1)>0.5) {
5272 delete array->RemoveAt(i);
5276 if (track0->GetKinkIndex(0)!=0) continue;
5277 if (track0->GetNumberOfClusters()<80) continue;
5279 AliTPCseed *pmother = new AliTPCseed();
5280 AliTPCseed *pdaughter = new AliTPCseed();
5281 AliKink *pkink = new AliKink;
5283 AliTPCseed & mother = *pmother;
5284 AliTPCseed & daughter = *pdaughter;
5285 AliKink & kink = *pkink;
5286 if (CheckKinkPoint(track0,mother,daughter, kink)){
5287 if (mother.GetNumberOfClusters()<30||daughter.GetNumberOfClusters()<20) {
5291 continue; //too short tracks
5293 if (mother.Pt()<1.4) {
5299 Int_t row0= kink.GetTPCRow0();
5300 if (kink.GetDistance()>0.5 || kink.GetR()<110. || kink.GetR()>240.) {
5307 Int_t index = esd->AddKink(&kink);
5308 mother.SetKinkIndex(0,-(index+1));
5309 daughter.SetKinkIndex(0,index+1);
5310 if (mother.GetNumberOfClusters()>50) {
5311 delete array->RemoveAt(i);
5312 array->AddAt(new AliTPCseed(mother),i);
5315 array->AddLast(new AliTPCseed(mother));
5317 array->AddLast(new AliTPCseed(daughter));
5318 for (Int_t icl=0;icl<row0;icl++) {
5319 if (mother.GetClusterPointer(icl)) mother.GetClusterPointer(icl)->Use(20);
5322 for (Int_t icl=row0;icl<158;icl++) {
5323 if (daughter.GetClusterPointer(icl)) daughter.GetClusterPointer(icl)->Use(20);
5332 delete [] daughters;
5354 printf("Ncandidates=\t%d\t%d\t%d\t%d\n",esd->GetNumberOfKinks(),ncandidates,ntracks,nall);
5358 void AliTPCtrackerMI::FindV0s(TObjArray * array, AliESDEvent *esd)
5364 TObjArray *tpcv0s = new TObjArray(100000);
5365 Int_t nentries = array->GetEntriesFast();
5366 AliHelix *helixes = new AliHelix[nentries];
5367 Int_t *sign = new Int_t[nentries];
5368 Float_t *alpha = new Float_t[nentries];
5369 Float_t *z0 = new Float_t[nentries];
5370 Float_t *dca = new Float_t[nentries];
5371 Float_t *sdcar = new Float_t[nentries];
5372 Float_t *cdcar = new Float_t[nentries];
5373 Float_t *pulldcar = new Float_t[nentries];
5374 Float_t *pulldcaz = new Float_t[nentries];
5375 Float_t *pulldca = new Float_t[nentries];
5376 Bool_t *isPrim = new Bool_t[nentries];
5377 const AliESDVertex * primvertex = esd->GetVertex();
5378 Double_t zvertex = primvertex->GetZv();
5380 // nentries = array->GetEntriesFast();
5382 for (Int_t i=0;i<nentries;i++){
5385 AliTPCseed* track = (AliTPCseed*)array->At(i);
5386 if (!track) continue;
5387 track->GetV0Indexes()[0] = 0; //rest v0 indexes
5388 track->GetV0Indexes()[1] = 0; //rest v0 indexes
5389 track->GetV0Indexes()[2] = 0; //rest v0 indexes
5391 alpha[i] = track->GetAlpha();
5392 new (&helixes[i]) AliHelix(*track);
5394 helixes[i].Evaluate(0,xyz);
5395 sign[i] = (track->GetC()>0) ? -1:1;
5399 if (track->GetProlongation(0,y,z)) z0[i] = z;
5400 dca[i] = track->GetD(0,0);
5402 // dca error parrameterezation + pulls
5404 sdcar[i] = TMath::Sqrt(0.150*0.150+(100*track->GetC())*(100*track->GetC()));
5405 if (TMath::Abs(track->GetTgl())>1) sdcar[i]*=2.5;
5406 cdcar[i] = TMath::Exp((TMath::Abs(track->GetC())-0.0106)*525.3);
5407 pulldcar[i] = (dca[i]-cdcar[i])/sdcar[i];
5408 pulldcaz[i] = (z0[i]-zvertex)/sdcar[i];
5409 pulldca[i] = TMath::Sqrt(pulldcar[i]*pulldcar[i]+pulldcaz[i]*pulldcaz[i]);
5410 if (track->TPCrPID(1)+track->TPCrPID(2)+track->TPCrPID(3)>0.5) {
5411 if (pulldca[i]<3.) isPrim[i]=kTRUE; //pion, muon and Kaon 3 sigma cut
5413 if (track->TPCrPID(4)>0.5) {
5414 if (pulldca[i]<0.5) isPrim[i]=kTRUE; //proton 0.5 sigma cut
5416 if (track->TPCrPID(0)>0.4) {
5417 isPrim[i]=kFALSE; //electron no sigma cut
5424 Int_t ncandidates =0;
5427 Double_t phase[2][2],radius[2];
5433 TTreeSRedirector &cstream = *fDebugStreamer;
5434 Float_t fprimvertex[3]={GetX(),GetY(),GetZ()};
5436 Double_t cradius0 = 10*10;
5437 Double_t cradius1 = 200*200;
5440 Double_t cpointAngle = 0.95;
5442 Double_t delta[2]={10000,10000};
5443 for (Int_t i =0;i<nentries;i++){
5444 if (sign[i]==0) continue;
5445 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5446 if (!track0) continue;
5447 if (AliTPCReconstructor::StreamLevel()>1){
5452 "zvertex="<<zvertex<<
5453 "sdcar0="<<sdcar[i]<<
5454 "cdcar0="<<cdcar[i]<<
5455 "pulldcar0="<<pulldcar[i]<<
5456 "pulldcaz0="<<pulldcaz[i]<<
5457 "pulldca0="<<pulldca[i]<<
5458 "isPrim="<<isPrim[i]<<
5462 if (track0->GetSigned1Pt()<0) continue;
5463 if (track0->GetKinkIndex(0)>0||isPrim[i]) continue; //daughter kink
5465 if (TMath::Abs(helixes[i].GetHelix(4))<0.000000001) continue;
5470 for (Int_t j =0;j<nentries;j++){
5471 AliTPCseed * track1 = (AliTPCseed*)array->At(j);
5472 if (!track1) continue;
5473 if (track1->GetKinkIndex(0)>0 || isPrim[j]) continue; //daughter kink
5474 if (sign[j]*sign[i]>0) continue;
5475 if (TMath::Abs(helixes[j].GetHelix(4))<0.000001) continue;
5476 if (track0->GetCircular()+track1->GetCircular()>1) continue; //circling -returning track
5479 // DCA to prim vertex cut
5485 Int_t npoints = helixes[i].GetRPHIintersections(helixes[j], phase, radius,cdist2);
5486 if (npoints<1) continue;
5490 if (radius[0]<cradius0||radius[0]>cradius1) continue;
5491 helixes[i].LinearDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta[0]);
5492 if (delta[0]>cdist1) continue;
5495 if (TMath::Max(radius[0],radius[1])<cradius0|| TMath::Min(radius[0],radius[1])>cradius1) continue;
5496 helixes[i].LinearDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta[0]);
5497 helixes[i].LinearDCA(helixes[j],phase[1][0],phase[1][1],radius[1],delta[1]);
5498 if (delta[1]<delta[0]) iclosest=1;
5499 if (delta[iclosest]>cdist1) continue;
5501 helixes[i].ParabolicDCA(helixes[j],phase[iclosest][0],phase[iclosest][1],radius[iclosest],delta[iclosest]);
5502 if (radius[iclosest]<cradius0 || radius[iclosest]>cradius1 || delta[iclosest]>cdist1) continue;
5504 Double_t pointAngle = helixes[i].GetPointAngle(helixes[j],phase[iclosest],fprimvertex);
5505 if (pointAngle<cpointAngle) continue;
5507 Bool_t isGamma = kFALSE;
5508 vertex.SetParamP(*track0); //track0 - plus
5509 vertex.SetParamN(*track1); //track1 - minus
5510 vertex.Update(fprimvertex);
5511 if (track0->TPCrPID(0)>0.3&&track1->TPCrPID(0)>0.3&&vertex.GetAnglep()[2]<0.15) isGamma=kTRUE; // gamma conversion candidate
5512 Double_t pointAngle2 = vertex.GetV0CosineOfPointingAngle();
5514 if (vertex.GetV0CosineOfPointingAngle()<cpointAngle && (!isGamma)) continue;// point angle cut
5515 //Bo: if (vertex.GetDist2()>2&&(!isGamma)) continue; // point angle cut
5516 if (vertex.GetDcaV0Daughters()>2&&(!isGamma)) continue;//Bo: // point angle cut
5517 Float_t sigmae = 0.15*0.15;
5518 if (vertex.GetRr()<80)
5519 sigmae += (sdcar[i]*sdcar[i]+sdcar[j]*sdcar[j])*(1.-vertex.GetRr()/80.)*(1.-vertex.GetRr()/80.);
5520 sigmae+= TMath::Sqrt(sigmae);
5521 //Bo: if (vertex.GetDist2()/sigmae>3.&&(!isGamma)) continue;
5522 if (vertex.GetDcaV0Daughters()/sigmae>3.&&(!isGamma)) continue;
5523 Float_t densb0=0,densb1=0,densa0=0,densa1=0;
5524 Int_t row0 = GetRowNumber(vertex.GetRr());
5526 //Bo: if (vertex.GetDist2()>0.2) continue;
5527 if (vertex.GetDcaV0Daughters()>0.2) continue;
5528 densb0 = track0->Density2(0,row0-5);
5529 densb1 = track1->Density2(0,row0-5);
5530 if (densb0>0.3|| densb1>0.3) continue; //clusters before vertex
5531 densa0 = track0->Density2(row0+5,row0+40);
5532 densa1 = track1->Density2(row0+5,row0+40);
5533 if ((densa0<0.4|| densa1<0.4)&&(!isGamma)) continue; //missing clusters after vertex
5536 densa0 = track0->Density2(0,40); //cluster density
5537 densa1 = track1->Density2(0,40); //cluster density
5538 if ((vertex.GetRr()<80&&densa0+densa1<1.)&&(!isGamma)) continue;
5540 //Bo: vertex.SetLab(0,track0->GetLabel());
5541 //Bo: vertex.SetLab(1,track1->GetLabel());
5542 vertex.SetChi2After((densa0+densa1)*0.5);
5543 vertex.SetChi2Before((densb0+densb1)*0.5);
5544 vertex.SetIndex(0,i);
5545 vertex.SetIndex(1,j);
5546 //Bo: vertex.SetStatus(1); // TPC v0 candidate
5547 vertex.SetOnFlyStatus(2);//Bo: // TPC v0 candidate
5548 //Bo: vertex.SetRp(track0->TPCrPIDs());
5549 //Bo: vertex.SetRm(track1->TPCrPIDs());
5550 tpcv0s->AddLast(new AliESDv0(vertex));
5553 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
5554 Double_t radiusm= (delta[0]<delta[1])? TMath::Sqrt(radius[0]):TMath::Sqrt(radius[1]);
5555 Double_t deltam= (delta[0]<delta[1])? TMath::Sqrt(delta[0]):TMath::Sqrt(delta[1]);
5556 if (AliTPCReconstructor::StreamLevel()>1) {
5557 Int_t lab0=track0->GetLabel();
5558 Int_t lab1=track1->GetLabel();
5559 Char_t c0=track0->GetCircular();
5560 Char_t c1=track1->GetCircular();
5563 "vertex.="<<&vertex<<
5566 "Helix0.="<<&helixes[i]<<
5569 "Helix1.="<<&helixes[j]<<
5570 "pointAngle="<<pointAngle<<
5571 "pointAngle2="<<pointAngle2<<
5576 "zvertex="<<zvertex<<
5579 "npoints="<<npoints<<
5580 "radius0="<<radius[0]<<
5581 "delta0="<<delta[0]<<
5582 "radius1="<<radius[1]<<
5583 "delta1="<<delta[1]<<
5584 "radiusm="<<radiusm<<
5586 "sdcar0="<<sdcar[i]<<
5587 "sdcar1="<<sdcar[j]<<
5588 "cdcar0="<<cdcar[i]<<
5589 "cdcar1="<<cdcar[j]<<
5590 "pulldcar0="<<pulldcar[i]<<
5591 "pulldcar1="<<pulldcar[j]<<
5592 "pulldcaz0="<<pulldcaz[i]<<
5593 "pulldcaz1="<<pulldcaz[j]<<
5594 "pulldca0="<<pulldca[i]<<
5595 "pulldca1="<<pulldca[j]<<
5605 Float_t *quality = new Float_t[ncandidates];
5606 Int_t *indexes = new Int_t[ncandidates];
5608 for (Int_t i=0;i<ncandidates;i++){
5610 AliESDv0 *v0 = (AliESDv0*)tpcv0s->At(i);
5611 quality[i] = 1./(1.00001-v0->GetV0CosineOfPointingAngle()); //base point angle
5612 // quality[i] /= (0.5+v0->GetDist2());
5613 // quality[i] *= v0->GetChi2After(); //density factor
5615 Int_t index0 = v0->GetIndex(0);
5616 Int_t index1 = v0->GetIndex(1);
5617 //Bo: Double_t minpulldca = TMath::Min(2.+pulldca[v0->GetIndex(0)],(2.+pulldca[v0->GetIndex(1)]) ); //pull
5618 Double_t minpulldca = TMath::Min(2.+pulldca[index0],(2.+pulldca[index1]) );//Bo:
5622 AliTPCseed * track0 = (AliTPCseed*)array->At(index0);
5623 AliTPCseed * track1 = (AliTPCseed*)array->At(index1);
5624 if (track0->TPCrPID(0)>0.3&&track1->TPCrPID(0)>0.3&&v0->GetAnglep()[2]<0.15) quality[i]+=1000000; // gamma conversion candidate
5625 if (track0->TPCrPID(4)>0.9||track1->TPCrPID(4)>0.9&&minpulldca>4) quality[i]*=10; // lambda candidate candidate
5628 TMath::Sort(ncandidates,quality,indexes,kTRUE);
5631 for (Int_t i=0;i<ncandidates;i++){
5632 AliESDv0 * v0 = (AliESDv0*)tpcv0s->At(indexes[i]);
5634 Int_t index0 = v0->GetIndex(0);
5635 Int_t index1 = v0->GetIndex(1);
5636 AliTPCseed * track0 = (AliTPCseed*)array->At(index0);
5637 AliTPCseed * track1 = (AliTPCseed*)array->At(index1);
5638 if (!track0||!track1) {
5642 Bool_t accept =kTRUE; //default accept
5643 Int_t *v0indexes0 = track0->GetV0Indexes();
5644 Int_t *v0indexes1 = track1->GetV0Indexes();
5646 Int_t order0 = (v0indexes0[0]!=0) ? 1:0;
5647 Int_t order1 = (v0indexes1[0]!=0) ? 1:0;
5648 if (v0indexes0[1]!=0) order0 =2;
5649 if (v0indexes1[1]!=0) order1 =2;
5651 if (v0indexes0[2]!=0) {order0=3; accept=kFALSE;}
5652 if (v0indexes0[2]!=0) {order1=3; accept=kFALSE;}
5654 AliESDv0 * v02 = v0;
5656 //Bo: v0->SetOrder(0,order0);
5657 //Bo: v0->SetOrder(1,order1);
5658 //Bo: v0->SetOrder(1,order0+order1);
5659 v0->SetOnFlyStatus(kTRUE);
5660 Int_t index = esd->AddV0(v0);
5661 v02 = esd->GetV0(index);
5662 v0indexes0[order0]=index;
5663 v0indexes1[order1]=index;
5667 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
5668 if (AliTPCReconstructor::StreamLevel()>1) {
5669 Int_t lab0=track0->GetLabel();
5670 Int_t lab1=track1->GetLabel();
5679 "dca0="<<dca[index0]<<
5680 "dca1="<<dca[index1]<<
5684 "quality="<<quality[i]<<
5685 "pulldca0="<<pulldca[index0]<<
5686 "pulldca1="<<pulldca[index1]<<
5710 printf("TPC V0 finder : naccepted\t%d\tncandidates\t%d\tntracks\t%d\tnall\t%d\n",naccepted,ncandidates,ntracks,nall);
5714 Int_t AliTPCtrackerMI::RefitKink(AliTPCseed &mother, AliTPCseed &daughter, AliESDkink &knk)
5717 // refit kink towards to the vertex
5720 AliKink &kink=(AliKink &)knk;
5722 Int_t row0 = GetRowNumber(kink.GetR());
5723 FollowProlongation(mother,0);
5724 mother.Reset(kFALSE);
5726 FollowProlongation(daughter,row0);
5727 daughter.Reset(kFALSE);
5728 FollowBackProlongation(daughter,158);
5729 daughter.Reset(kFALSE);
5730 Int_t first = TMath::Max(row0-20,30);
5731 Int_t last = TMath::Min(row0+20,140);
5733 const Int_t kNdiv =5;
5734 AliTPCseed param0[kNdiv]; // parameters along the track
5735 AliTPCseed param1[kNdiv]; // parameters along the track
5736 AliKink kinks[kNdiv]; // corresponding kink parameters
5739 for (Int_t irow=0; irow<kNdiv;irow++){
5740 rows[irow] = first +((last-first)*irow)/(kNdiv-1);
5742 // store parameters along the track
5744 for (Int_t irow=0;irow<kNdiv;irow++){
5745 FollowBackProlongation(mother, rows[irow]);
5746 FollowProlongation(daughter,rows[kNdiv-1-irow]);
5747 param0[irow] = mother;
5748 param1[kNdiv-1-irow] = daughter;
5752 for (Int_t irow=0; irow<kNdiv-1;irow++){
5753 if (param0[irow].GetNumberOfClusters()<kNdiv||param1[irow].GetNumberOfClusters()<kNdiv) continue;
5754 kinks[irow].SetMother(param0[irow]);
5755 kinks[irow].SetDaughter(param1[irow]);
5756 kinks[irow].Update();
5759 // choose kink with best "quality"
5761 Double_t mindist = 10000;
5762 for (Int_t irow=0;irow<kNdiv;irow++){
5763 if (param0[irow].GetNumberOfClusters()<20||param1[irow].GetNumberOfClusters()<20) continue;
5764 if (TMath::Abs(kinks[irow].GetR())>240.) continue;
5765 if (TMath::Abs(kinks[irow].GetR())<100.) continue;
5767 Float_t normdist = TMath::Abs(param0[irow].GetX()-kinks[irow].GetR())*(0.1+kink.GetDistance());
5768 normdist/= (param0[irow].GetNumberOfClusters()+param1[irow].GetNumberOfClusters()+40.);
5769 if (normdist < mindist){
5775 if (index==-1) return 0;
5778 param0[index].Reset(kTRUE);
5779 FollowProlongation(param0[index],0);
5781 mother = param0[index];
5782 daughter = param1[index]; // daughter in vertex
5784 kink.SetMother(mother);
5785 kink.SetDaughter(daughter);
5787 kink.SetTPCRow0(GetRowNumber(kink.GetR()));
5788 kink.SetTPCncls(param0[index].GetNumberOfClusters(),0);
5789 kink.SetTPCncls(param1[index].GetNumberOfClusters(),1);
5790 kink.SetLabel(CookLabel(&mother,0.4, 0,kink.GetTPCRow0()),0);
5791 kink.SetLabel(CookLabel(&daughter,0.4, kink.GetTPCRow0(),160),1);
5792 mother.SetLabel(kink.GetLabel(0));
5793 daughter.SetLabel(kink.GetLabel(1));
5799 void AliTPCtrackerMI::UpdateKinkQualityM(AliTPCseed * seed){
5801 // update Kink quality information for mother after back propagation
5803 if (seed->GetKinkIndex(0)>=0) return;
5804 for (Int_t ikink=0;ikink<3;ikink++){
5805 Int_t index = seed->GetKinkIndex(ikink);
5806 if (index>=0) break;
5807 index = TMath::Abs(index)-1;
5808 AliESDkink * kink = fEvent->GetKink(index);
5809 kink->SetTPCDensity(-1,0,0);
5810 kink->SetTPCDensity(1,0,1);
5812 Int_t row0 = kink->GetTPCRow0() - 2 - Int_t( 0.5/ (0.05+kink->GetAngle(2)));
5813 if (row0<15) row0=15;
5815 Int_t row1 = kink->GetTPCRow0() + 2 + Int_t( 0.5/ (0.05+kink->GetAngle(2)));
5816 if (row1>145) row1=145;
5818 Int_t found,foundable,shared;
5819 seed->GetClusterStatistic(0,row0, found, foundable,shared,kFALSE);
5820 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),0,0);
5821 seed->GetClusterStatistic(row1,155, found, foundable,shared,kFALSE);
5822 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),0,1);
5827 void AliTPCtrackerMI::UpdateKinkQualityD(AliTPCseed * seed){
5829 // update Kink quality information for daughter after refit
5831 if (seed->GetKinkIndex(0)<=0) return;
5832 for (Int_t ikink=0;ikink<3;ikink++){
5833 Int_t index = seed->GetKinkIndex(ikink);
5834 if (index<=0) break;
5835 index = TMath::Abs(index)-1;
5836 AliESDkink * kink = fEvent->GetKink(index);
5837 kink->SetTPCDensity(-1,1,0);
5838 kink->SetTPCDensity(-1,1,1);
5840 Int_t row0 = kink->GetTPCRow0() -2 - Int_t( 0.5/ (0.05+kink->GetAngle(2)));
5841 if (row0<15) row0=15;
5843 Int_t row1 = kink->GetTPCRow0() +2 + Int_t( 0.5/ (0.05+kink->GetAngle(2)));
5844 if (row1>145) row1=145;
5846 Int_t found,foundable,shared;
5847 seed->GetClusterStatistic(0,row0, found, foundable,shared,kFALSE);
5848 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),1,0);
5849 seed->GetClusterStatistic(row1,155, found, foundable,shared,kFALSE);
5850 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),1,1);
5856 Int_t AliTPCtrackerMI::CheckKinkPoint(AliTPCseed*seed,AliTPCseed &mother, AliTPCseed &daughter, AliESDkink &knk)
5859 // check kink point for given track
5860 // if return value=0 kink point not found
5861 // otherwise seed0 correspond to mother particle
5862 // seed1 correspond to daughter particle
5863 // kink parameter of kink point
5864 AliKink &kink=(AliKink &)knk;
5866 Int_t middlerow = (seed->GetFirstPoint()+seed->GetLastPoint())/2;
5867 Int_t first = seed->GetFirstPoint();
5868 Int_t last = seed->GetLastPoint();
5869 if (last-first<20) return 0; // shortest length - 2*30 = 60 pad-rows
5872 AliTPCseed *seed1 = ReSeed(seed,middlerow+20, kTRUE); //middle of chamber
5873 if (!seed1) return 0;
5874 FollowProlongation(*seed1,seed->GetLastPoint()-20);
5875 seed1->Reset(kTRUE);
5876 FollowProlongation(*seed1,158);
5877 seed1->Reset(kTRUE);
5878 last = seed1->GetLastPoint();
5880 AliTPCseed *seed0 = new AliTPCseed(*seed);
5881 seed0->Reset(kFALSE);
5884 AliTPCseed param0[20]; // parameters along the track
5885 AliTPCseed param1[20]; // parameters along the track
5886 AliKink kinks[20]; // corresponding kink parameters
5888 for (Int_t irow=0; irow<20;irow++){
5889 rows[irow] = first +((last-first)*irow)/19;
5891 // store parameters along the track
5893 for (Int_t irow=0;irow<20;irow++){
5894 FollowBackProlongation(*seed0, rows[irow]);
5895 FollowProlongation(*seed1,rows[19-irow]);
5896 param0[irow] = *seed0;
5897 param1[19-irow] = *seed1;
5901 for (Int_t irow=0; irow<19;irow++){
5902 kinks[irow].SetMother(param0[irow]);
5903 kinks[irow].SetDaughter(param1[irow]);
5904 kinks[irow].Update();
5907 // choose kink with biggest change of angle
5909 Double_t maxchange= 0;
5910 for (Int_t irow=1;irow<19;irow++){
5911 if (TMath::Abs(kinks[irow].GetR())>240.) continue;
5912 if (TMath::Abs(kinks[irow].GetR())<110.) continue;
5913 Float_t quality = TMath::Abs(kinks[irow].GetAngle(2))/(3.+TMath::Abs(kinks[irow].GetR()-param0[irow].GetX()));
5914 if ( quality > maxchange){
5915 maxchange = quality;
5922 if (index<0) return 0;
5924 Int_t row0 = GetRowNumber(kinks[index].GetR()); //row 0 estimate
5925 seed0 = new AliTPCseed(param0[index]);
5926 seed1 = new AliTPCseed(param1[index]);
5927 seed0->Reset(kFALSE);
5928 seed1->Reset(kFALSE);
5929 seed0->ResetCovariance(10.);
5930 seed1->ResetCovariance(10.);
5931 FollowProlongation(*seed0,0);
5932 FollowBackProlongation(*seed1,158);
5933 mother = *seed0; // backup mother at position 0
5934 seed0->Reset(kFALSE);
5935 seed1->Reset(kFALSE);
5936 seed0->ResetCovariance(10.);
5937 seed1->ResetCovariance(10.);
5939 first = TMath::Max(row0-20,0);
5940 last = TMath::Min(row0+20,158);
5942 for (Int_t irow=0; irow<20;irow++){
5943 rows[irow] = first +((last-first)*irow)/19;
5945 // store parameters along the track
5947 for (Int_t irow=0;irow<20;irow++){
5948 FollowBackProlongation(*seed0, rows[irow]);
5949 FollowProlongation(*seed1,rows[19-irow]);
5950 param0[irow] = *seed0;
5951 param1[19-irow] = *seed1;
5955 for (Int_t irow=0; irow<19;irow++){
5956 kinks[irow].SetMother(param0[irow]);
5957 kinks[irow].SetDaughter(param1[irow]);
5958 // param0[irow].Dump();
5959 //param1[irow].Dump();
5960 kinks[irow].Update();
5963 // choose kink with biggest change of angle
5966 for (Int_t irow=0;irow<20;irow++){
5967 if (TMath::Abs(kinks[irow].GetR())>250.) continue;
5968 if (TMath::Abs(kinks[irow].GetR())<90.) continue;
5969 Float_t quality = TMath::Abs(kinks[irow].GetAngle(2))/(3.+TMath::Abs(kinks[irow].GetR()-param0[irow].GetX()));
5970 if ( quality > maxchange){
5971 maxchange = quality;
5978 if (index==-1 || param0[index].GetNumberOfClusters()+param1[index].GetNumberOfClusters()<100){
5983 // Float_t anglesigma = TMath::Sqrt(param0[index].fC22+param0[index].fC33+param1[index].fC22+param1[index].fC33);
5985 kink.SetMother(param0[index]);
5986 kink.SetDaughter(param1[index]);
5988 row0 = GetRowNumber(kink.GetR());
5989 kink.SetTPCRow0(row0);
5990 kink.SetLabel(CookLabel(seed0,0.5,0,row0),0);
5991 kink.SetLabel(CookLabel(seed1,0.5,row0,158),1);
5992 kink.SetIndex(-10,0);
5993 kink.SetIndex(int(param0[index].GetNumberOfClusters()+param1[index].GetNumberOfClusters()),1);
5994 kink.SetTPCncls(param0[index].GetNumberOfClusters(),0);
5995 kink.SetTPCncls(param1[index].GetNumberOfClusters(),1);
5998 // new (&mother) AliTPCseed(param0[index]);
5999 daughter = param1[index];
6000 daughter.SetLabel(kink.GetLabel(1));
6001 param0[index].Reset(kTRUE);
6002 FollowProlongation(param0[index],0);
6003 mother = param0[index];
6004 mother.SetLabel(kink.GetLabel(0));
6014 AliTPCseed* AliTPCtrackerMI::ReSeed(AliTPCseed *t)
6017 // reseed - refit - track
6020 // Int_t last = fSectors->GetNRows()-1;
6022 if (fSectors == fOuterSec){
6023 first = TMath::Max(first, t->GetFirstPoint()-fInnerSec->GetNRows());
6027 first = t->GetFirstPoint();
6029 AliTPCseed * seed = MakeSeed(t,0.1,0.5,0.9);
6030 FollowBackProlongation(*t,fSectors->GetNRows()-1);
6032 FollowProlongation(*t,first);
6042 //_____________________________________________________________________________
6043 Int_t AliTPCtrackerMI::ReadSeeds(const TFile *inp) {
6044 //-----------------------------------------------------------------
6045 // This function reades track seeds.
6046 //-----------------------------------------------------------------
6047 TDirectory *savedir=gDirectory;
6049 TFile *in=(TFile*)inp;
6050 if (!in->IsOpen()) {
6051 cerr<<"AliTPCtrackerMI::ReadSeeds(): input file is not open !\n";
6056 TTree *seedTree=(TTree*)in->Get("Seeds");
6058 cerr<<"AliTPCtrackerMI::ReadSeeds(): ";
6059 cerr<<"can't get a tree with track seeds !\n";
6062 AliTPCtrack *seed=new AliTPCtrack;
6063 seedTree->SetBranchAddress("tracks",&seed);
6065 if (fSeeds==0) fSeeds=new TObjArray(15000);
6067 Int_t n=(Int_t)seedTree->GetEntries();
6068 for (Int_t i=0; i<n; i++) {
6069 seedTree->GetEvent(i);
6070 fSeeds->AddLast(new AliTPCseed(*seed/*,seed->GetAlpha()*/));
6079 Int_t AliTPCtrackerMI::Clusters2Tracks (AliESDEvent *esd)
6082 if (fSeeds) DeleteSeeds();
6085 if (!fSeeds) return 1;
6092 //_____________________________________________________________________________
6093 Int_t AliTPCtrackerMI::Clusters2Tracks() {
6094 //-----------------------------------------------------------------
6095 // This is a track finder.
6096 //-----------------------------------------------------------------
6097 TDirectory *savedir=gDirectory;
6101 fSeeds = Tracking();
6104 Info("Clusters2Tracks","Time for tracking: \t");timer.Print();timer.Start();
6106 //activate again some tracks
6107 for (Int_t i=0; i<fSeeds->GetEntriesFast(); i++) {
6108 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
6110 Int_t nc=t.GetNumberOfClusters();
6112 delete fSeeds->RemoveAt(i);
6116 if (pt->GetRemoval()==10) {
6117 if (pt->GetDensityFirst(20)>0.8 || pt->GetDensityFirst(30)>0.8 || pt->GetDensityFirst(40)>0.7)
6118 pt->Desactivate(10); // make track again active
6120 pt->Desactivate(20);
6121 delete fSeeds->RemoveAt(i);
6126 RemoveUsed2(fSeeds,0.85,0.85,0);
6127 if (AliTPCReconstructor::GetRecoParam()->GetDoKinks()) FindKinks(fSeeds,fEvent);
6128 //FindCurling(fSeeds, fEvent,0);
6129 if (AliTPCReconstructor::StreamLevel()>2) FindMultiMC(fSeeds, fEvent,0); // find multi found tracks
6130 RemoveUsed2(fSeeds,0.5,0.4,20);
6131 FindSplitted(fSeeds, fEvent,0); // find multi found tracks
6134 // // refit short tracks
6136 Int_t nseed=fSeeds->GetEntriesFast();
6139 for (Int_t i=0; i<nseed; i++) {
6140 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
6142 Int_t nc=t.GetNumberOfClusters();
6144 delete fSeeds->RemoveAt(i);
6147 CookLabel(pt,0.1); //For comparison only
6148 //if ((pt->IsActive() || (pt->fRemoval==10) )&& nc>50 &&pt->GetNumberOfClusters()>0.4*pt->fNFoundable){
6149 if ((pt->IsActive() || (pt->GetRemoval()==10) )){
6151 if (fDebug>0) cerr<<found<<'\r';
6155 delete fSeeds->RemoveAt(i);
6159 //RemoveOverlap(fSeeds,0.99,7,kTRUE);
6161 //RemoveUsed(fSeeds,0.9,0.9,6);
6163 nseed=fSeeds->GetEntriesFast();
6165 for (Int_t i=0; i<nseed; i++) {
6166 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
6168 Int_t nc=t.GetNumberOfClusters();
6170 delete fSeeds->RemoveAt(i);
6174 t.CookdEdx(0.02,0.6);
6175 // CheckKinkPoint(&t,0.05);
6176 //if ((pt->IsActive() || (pt->fRemoval==10) )&& nc>50 &&pt->GetNumberOfClusters()>0.4*pt->fNFoundable){
6177 if ((pt->IsActive() || (pt->GetRemoval()==10) )){
6185 delete fSeeds->RemoveAt(i);
6186 //AliTPCseed * seed1 = ReSeed(pt,0.05,0.5,1);
6188 // FollowProlongation(*seed1,0);
6189 // Int_t n = seed1->GetNumberOfClusters();
6190 // printf("fP4\t%f\t%f\n",seed1->GetC(),pt->GetC());
6191 // printf("fN\t%d\t%d\n", seed1->GetNumberOfClusters(),pt->GetNumberOfClusters());
6194 //AliTPCseed * seed2 = ReSeed(pt,0.95,0.5,0.05);
6198 SortTracks(fSeeds, 1);
6202 PrepareForBackProlongation(fSeeds,5.);
6203 PropagateBack(fSeeds);
6204 printf("Time for back propagation: \t");timer.Print();timer.Start();
6208 PrepareForProlongation(fSeeds,5.);
6209 PropagateForward2(fSeeds);
6211 printf("Time for FORWARD propagation: \t");timer.Print();timer.Start();
6212 // RemoveUsed(fSeeds,0.7,0.7,6);
6213 //RemoveOverlap(fSeeds,0.9,7,kTRUE);
6215 nseed=fSeeds->GetEntriesFast();
6217 for (Int_t i=0; i<nseed; i++) {
6218 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
6220 Int_t nc=t.GetNumberOfClusters();
6222 delete fSeeds->RemoveAt(i);
6225 t.CookdEdx(0.02,0.6);
6226 // CookLabel(pt,0.1); //For comparison only
6227 //if ((pt->IsActive() || (pt->fRemoval==10) )&& nc>50 &&pt->GetNumberOfClusters()>0.4*pt->fNFoundable){
6228 if ((pt->IsActive() || (pt->fRemoval==10) )){
6229 cerr<<found++<<'\r';
6232 delete fSeeds->RemoveAt(i);
6237 // fNTracks = found;
6239 Info("Clusters2Tracks","Time for overlap removal, track writing and dedx cooking: \t"); timer.Print();timer.Start();
6242 // cerr<<"Number of found tracks : "<<"\t"<<found<<endl;
6243 Info("Clusters2Tracks","Number of found tracks %d",found);
6245 // UnloadClusters();
6250 void AliTPCtrackerMI::Tracking(TObjArray * arr)
6253 // tracking of the seeds
6256 fSectors = fOuterSec;
6257 ParallelTracking(arr,150,63);
6258 fSectors = fOuterSec;
6259 ParallelTracking(arr,63,0);
6262 TObjArray * AliTPCtrackerMI::Tracking(Int_t seedtype, Int_t i1, Int_t i2, Float_t cuts[4], Float_t dy, Int_t dsec)
6267 TObjArray * arr = new TObjArray;
6269 fSectors = fOuterSec;
6272 for (Int_t sec=0;sec<fkNOS;sec++){
6273 if (seedtype==3) MakeSeeds3(arr,sec,i1,i2,cuts,dy, dsec);
6274 if (seedtype==4) MakeSeeds5(arr,sec,i1,i2,cuts,dy);
6275 if (seedtype==2) MakeSeeds2(arr,sec,i1,i2,cuts,dy);
6278 Info("Tracking","\nSeeding - %d\t%d\t%d\t%d\n",seedtype,i1,i2,arr->GetEntriesFast());
6290 TObjArray * AliTPCtrackerMI::Tracking()
6294 if (AliTPCReconstructor::GetRecoParam()->GetSpecialSeeding()) return TrackingSpecial();
6297 Int_t nup=fOuterSec->GetNRows()+fInnerSec->GetNRows();
6299 TObjArray * seeds = new TObjArray;
6308 Float_t fnumber = 3.0;
6309 Float_t fdensity = 3.0;
6314 for (Int_t delta = 0; delta<18; delta+=6){
6318 arr = Tracking(3,nup-1-delta,nup-1-delta-gap,cuts,-1,1);
6319 SumTracks(seeds,arr);
6320 SignClusters(seeds,fnumber,fdensity);
6322 for (Int_t i=2;i<6;i+=2){
6323 // seed high pt tracks
6326 arr = Tracking(3,nup-i-delta,nup-i-delta-gap,cuts,-1,0);
6327 SumTracks(seeds,arr);
6328 SignClusters(seeds,fnumber,fdensity);
6333 // RemoveUsed(seeds,0.9,0.9,1);
6334 // UnsignClusters();
6335 // SignClusters(seeds,fnumber,fdensity);
6339 for (Int_t delta = 20; delta<120; delta+=10){
6341 // seed high pt tracks
6345 arr = Tracking(3,nup-delta,nup-delta-gap,cuts,-1);
6346 SumTracks(seeds,arr);
6347 SignClusters(seeds,fnumber,fdensity);
6352 arr = Tracking(3,nup-delta-5,nup-delta-5-gap,cuts,-1);
6353 SumTracks(seeds,arr);
6354 SignClusters(seeds,fnumber,fdensity);
6365 Info("Tracking()","\n\nPrimary seeding\t%d\n\n",seeds->GetEntriesFast());
6369 // RemoveUsed(seeds,0.75,0.75,1);
6371 //SignClusters(seeds,fnumber,fdensity);
6380 arr = Tracking(4,nup-1,nup-1-gap,cuts,-1);
6381 SumTracks(seeds,arr);
6382 SignClusters(seeds,fnumber,fdensity);
6384 arr = Tracking(4,nup-2,nup-2-gap,cuts,-1);
6385 SumTracks(seeds,arr);
6386 SignClusters(seeds,fnumber,fdensity);
6388 arr = Tracking(4,nup-3,nup-3-gap,cuts,-1);
6389 SumTracks(seeds,arr);
6390 SignClusters(seeds,fnumber,fdensity);
6394 for (Int_t delta = 3; delta<30; delta+=5){
6400 arr = Tracking(4,nup-1-delta,nup-1-delta-gap,cuts,-1);
6401 SumTracks(seeds,arr);
6402 SignClusters(seeds,fnumber,fdensity);
6404 arr = Tracking(4,nup-3-delta,nup-5-delta-gap,cuts,4);
6405 SumTracks(seeds,arr);
6406 SignClusters(seeds,fnumber,fdensity);
6417 Int_t fLastSeedRowSec=AliTPCReconstructor::GetRecoParam()->GetLastSeedRowSec();
6420 for (Int_t delta = 30; delta<fLastSeedRowSec; delta+=10){
6426 arr = Tracking(4,nup-1-delta,nup-1-delta-gap,cuts,-1);
6427 SumTracks(seeds,arr);
6428 SignClusters(seeds,fnumber,fdensity);
6430 arr = Tracking(4,nup-5-delta,nup-5-delta-gap,cuts,5 );
6431 SumTracks(seeds,arr);
6432 SignClusters(seeds,fnumber,fdensity);
6436 Info("Tracking()","\n\nSecondary seeding\t%d\n\n",seeds->GetEntriesFast());
6447 TObjArray * AliTPCtrackerMI::TrackingSpecial()
6450 // seeding adjusted for laser and cosmic tests - short tracks with big inclination angle
6451 // no primary vertex seeding tried
6455 Int_t nup=fOuterSec->GetNRows()+fInnerSec->GetNRows();
6457 TObjArray * seeds = new TObjArray;
6462 Float_t fnumber = 3.0;
6463 Float_t fdensity = 3.0;
6466 cuts[0] = AliTPCReconstructor::GetRecoParam()->GetMaxC(); // max curvature
6467 cuts[1] = 3.5; // max tan(phi) angle for seeding
6468 cuts[2] = 3.; // not used (cut on z primary vertex)
6469 cuts[3] = 3.5; // max tan(theta) angle for seeding
6471 for (Int_t delta = 0; nup-delta-gap-1>0; delta+=3){
6473 arr = Tracking(4,nup-1-delta,nup-1-delta-gap,cuts,-1);
6474 SumTracks(seeds,arr);
6475 SignClusters(seeds,fnumber,fdensity);
6479 Info("Tracking()","\n\nSecondary seeding\t%d\n\n",seeds->GetEntriesFast());
6490 void AliTPCtrackerMI::SumTracks(TObjArray *arr1,TObjArray *&arr2) const
6493 //sum tracks to common container
6494 //remove suspicious tracks
6495 Int_t nseed = arr2->GetEntriesFast();
6496 for (Int_t i=0;i<nseed;i++){
6497 AliTPCseed *pt=(AliTPCseed*)arr2->UncheckedAt(i);
6500 // remove tracks with too big curvature
6502 if (TMath::Abs(pt->GetC())>AliTPCReconstructor::GetRecoParam()->GetMaxC()){
6503 delete arr2->RemoveAt(i);
6506 // REMOVE VERY SHORT TRACKS
6507 if (pt->GetNumberOfClusters()<20){
6508 delete arr2->RemoveAt(i);
6511 // NORMAL ACTIVE TRACK
6512 if (pt->IsActive()){
6513 arr1->AddLast(arr2->RemoveAt(i));
6516 //remove not usable tracks
6517 if (pt->GetRemoval()!=10){
6518 delete arr2->RemoveAt(i);
6522 // ENABLE ONLY ENOUGH GOOD STOPPED TRACKS
6523 if (pt->GetDensityFirst(20)>0.8 || pt->GetDensityFirst(30)>0.8 || pt->GetDensityFirst(40)>0.7)
6524 arr1->AddLast(arr2->RemoveAt(i));
6526 delete arr2->RemoveAt(i);
6530 delete arr2; arr2 = 0;
6535 void AliTPCtrackerMI::ParallelTracking(TObjArray * arr, Int_t rfirst, Int_t rlast)
6538 // try to track in parralel
6540 Int_t nseed=arr->GetEntriesFast();
6541 //prepare seeds for tracking
6542 for (Int_t i=0; i<nseed; i++) {
6543 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i), &t=*pt;
6545 if (!t.IsActive()) continue;
6546 // follow prolongation to the first layer
6547 if ( (fSectors ==fInnerSec) || (t.GetFirstPoint()-fParam->GetNRowLow()>rfirst+1) )
6548 FollowProlongation(t, rfirst+1);
6553 for (Int_t nr=rfirst; nr>=rlast; nr--){
6554 if (nr<fInnerSec->GetNRows())
6555 fSectors = fInnerSec;
6557 fSectors = fOuterSec;
6558 // make indexes with the cluster tracks for given
6560 // find nearest cluster
6561 for (Int_t i=0; i<nseed; i++) {
6562 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i), &t=*pt;
6564 if (nr==80) pt->UpdateReference();
6565 if (!pt->IsActive()) continue;
6566 // if ( (fSectors ==fOuterSec) && (pt->fFirstPoint-fParam->GetNRowLow())<nr) continue;
6567 if (pt->GetRelativeSector()>17) {
6570 UpdateClusters(t,nr);
6572 // prolonagate to the nearest cluster - if founded
6573 for (Int_t i=0; i<nseed; i++) {
6574 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
6576 if (!pt->IsActive()) continue;
6577 // if ((fSectors ==fOuterSec) && (pt->fFirstPoint-fParam->GetNRowLow())<nr) continue;
6578 if (pt->GetRelativeSector()>17) {
6581 FollowToNextCluster(*pt,nr);
6586 void AliTPCtrackerMI::PrepareForBackProlongation(TObjArray * arr,Float_t fac) const
6590 // if we use TPC track itself we have to "update" covariance
6592 Int_t nseed= arr->GetEntriesFast();
6593 for (Int_t i=0;i<nseed;i++){
6594 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
6598 //rotate to current local system at first accepted point
6599 Int_t index = pt->GetClusterIndex2(pt->GetFirstPoint());
6600 Int_t sec = (index&0xff000000)>>24;
6602 Float_t angle1 = fInnerSec->GetAlpha()*sec+fInnerSec->GetAlphaShift();
6603 if (angle1>TMath::Pi())
6604 angle1-=2.*TMath::Pi();
6605 Float_t angle2 = pt->GetAlpha();
6607 if (TMath::Abs(angle1-angle2)>0.001){
6608 pt->Rotate(angle1-angle2);
6609 //angle2 = pt->GetAlpha();
6610 //pt->fRelativeSector = pt->GetAlpha()/fInnerSec->GetAlpha();
6611 //if (pt->GetAlpha()<0)
6612 // pt->fRelativeSector+=18;
6613 //sec = pt->fRelativeSector;
6622 void AliTPCtrackerMI::PrepareForProlongation(TObjArray * arr, Float_t fac) const
6626 // if we use TPC track itself we have to "update" covariance
6628 Int_t nseed= arr->GetEntriesFast();
6629 for (Int_t i=0;i<nseed;i++){
6630 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
6633 pt->SetFirstPoint(pt->GetLastPoint());
6641 Int_t AliTPCtrackerMI::PropagateBack(TObjArray * arr)
6644 // make back propagation
6646 Int_t nseed= arr->GetEntriesFast();
6647 for (Int_t i=0;i<nseed;i++){
6648 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
6649 if (pt&& pt->GetKinkIndex(0)<=0) {
6650 //AliTPCseed *pt2 = new AliTPCseed(*pt);
6651 fSectors = fInnerSec;
6652 //FollowBackProlongation(*pt,fInnerSec->GetNRows()-1);
6653 //fSectors = fOuterSec;
6654 FollowBackProlongation(*pt,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1);
6655 //if (pt->GetNumberOfClusters()<(pt->fEsd->GetTPCclusters(0)) ){
6656 // Error("PropagateBack","Not prolonged track %d",pt->GetLabel());
6657 // FollowBackProlongation(*pt2,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1);
6660 if (pt&& pt->GetKinkIndex(0)>0) {
6661 AliESDkink * kink = fEvent->GetKink(pt->GetKinkIndex(0)-1);
6662 pt->SetFirstPoint(kink->GetTPCRow0());
6663 fSectors = fInnerSec;
6664 FollowBackProlongation(*pt,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1);
6672 Int_t AliTPCtrackerMI::PropagateForward2(TObjArray * arr)
6675 // make forward propagation
6677 Int_t nseed= arr->GetEntriesFast();
6679 for (Int_t i=0;i<nseed;i++){
6680 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
6682 FollowProlongation(*pt,0);
6691 Int_t AliTPCtrackerMI::PropagateForward()
6694 // propagate track forward
6696 Int_t nseed = fSeeds->GetEntriesFast();
6697 for (Int_t i=0;i<nseed;i++){
6698 AliTPCseed *pt = (AliTPCseed*)fSeeds->UncheckedAt(i);
6700 AliTPCseed &t = *pt;
6701 Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
6702 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
6703 if (alpha < 0. ) alpha += 2.*TMath::Pi();
6704 t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN);
6708 fSectors = fOuterSec;
6709 ParallelTracking(fSeeds,fOuterSec->GetNRows()+fInnerSec->GetNRows()-1,fInnerSec->GetNRows());
6710 fSectors = fInnerSec;
6711 ParallelTracking(fSeeds,fInnerSec->GetNRows()-1,0);
6720 Int_t AliTPCtrackerMI::PropagateBack(AliTPCseed * pt, Int_t row0, Int_t row1)
6723 // make back propagation, in between row0 and row1
6727 fSectors = fInnerSec;
6730 if (row1<fSectors->GetNRows())
6733 r1 = fSectors->GetNRows()-1;
6735 if (row0<fSectors->GetNRows()&& r1>0 )
6736 FollowBackProlongation(*pt,r1);
6737 if (row1<=fSectors->GetNRows())
6740 r1 = row1 - fSectors->GetNRows();
6741 if (r1<=0) return 0;
6742 if (r1>=fOuterSec->GetNRows()) return 0;
6743 fSectors = fOuterSec;
6744 return FollowBackProlongation(*pt,r1);
6752 void AliTPCtrackerMI::GetShape(AliTPCseed * seed, Int_t row)
6756 AliTPCClusterParam * clparam = AliTPCcalibDB::Instance()->GetClusterParam();
6757 Float_t zdrift = TMath::Abs((fParam->GetZLength(0)-TMath::Abs(seed->GetZ())));
6758 Int_t type = (seed->GetSector() < fParam->GetNSector()/2) ? 0: (row>126) ? 1:2;
6759 Double_t angulary = seed->GetSnp();
6760 angulary = angulary*angulary/(1.-angulary*angulary);
6761 Double_t angularz = seed->GetTgl()*seed->GetTgl()*(1.+angulary);
6763 Double_t sigmay = clparam->GetRMS0(0,type,zdrift,TMath::Sqrt(TMath::Abs(angulary)));
6764 Double_t sigmaz = clparam->GetRMS0(1,type,zdrift,TMath::Sqrt(TMath::Abs(angularz)));
6765 seed->SetCurrentSigmaY2(sigmay*sigmay);
6766 seed->SetCurrentSigmaZ2(sigmaz*sigmaz);
6767 // Float_t sd2 = TMath::Abs((fParam->GetZLength(0)-TMath::Abs(seed->GetZ())))*fParam->GetDiffL()*fParam->GetDiffL();
6768 // // Float_t padlength = fParam->GetPadPitchLength(seed->fSector);
6769 // Float_t padlength = GetPadPitchLength(row);
6771 // Float_t sresy = (seed->GetSector() < fParam->GetNSector()/2) ? 0.2 :0.3;
6772 // seed->SetCurrentSigmaY2(sd2+padlength*padlength*angulary/12.+sresy*sresy);
6774 // Float_t sresz = fParam->GetZSigma();
6775 // seed->SetCurrentSigmaZ2(sd2+padlength*padlength*angularz*angularz*(1+angulary)/12.+sresz*sresz);
6777 Float_t wy = GetSigmaY(seed);
6778 Float_t wz = GetSigmaZ(seed);
6781 if (TMath::Abs(wy/seed->fCurrentSigmaY2-1)>0.0001 || TMath::Abs(wz/seed->fCurrentSigmaZ2-1)>0.0001 ){
6782 printf("problem\n");
6789 //__________________________________________________________________________
6790 void AliTPCtrackerMI::CookLabel(AliKalmanTrack *tk, Float_t wrong) const {
6791 //--------------------------------------------------------------------
6792 //This function "cooks" a track label. If label<0, this track is fake.
6793 //--------------------------------------------------------------------
6794 AliTPCseed * t = dynamic_cast<AliTPCseed*>(tk);
6796 printf("%s:%d wrong type \n",(char*)__FILE__,__LINE__);
6800 Int_t noc=t->GetNumberOfClusters();
6802 //printf("\nnot founded prolongation\n\n\n");
6808 AliTPCclusterMI *clusters[160];
6810 for (Int_t i=0;i<160;i++) {
6817 for (i=0; i<160 && current<noc; i++) {
6819 Int_t index=t->GetClusterIndex2(i);
6820 if (index<=0) continue;
6821 if (index&0x8000) continue;
6823 //clusters[current]=GetClusterMI(index);
6824 if (t->GetClusterPointer(i)){
6825 clusters[current]=t->GetClusterPointer(i);
6831 Int_t lab=123456789;
6832 for (i=0; i<noc; i++) {
6833 AliTPCclusterMI *c=clusters[i];
6835 lab=TMath::Abs(c->GetLabel(0));
6837 for (j=0; j<noc; j++) if (lb[j]==lab || mx[j]==0) break;
6843 for (i=0; i<noc; i++) if (mx[i]>max) {max=mx[i]; lab=lb[i];}
6845 for (i=0; i<noc; i++) {
6846 AliTPCclusterMI *c=clusters[i];
6848 if (TMath::Abs(c->GetLabel(1)) == lab ||
6849 TMath::Abs(c->GetLabel(2)) == lab ) max++;
6852 if ((1.- Float_t(max)/noc) > wrong) lab=-lab;
6855 Int_t tail=Int_t(0.10*noc);
6858 for (i=1; i<=160&&ind<tail; i++) {
6859 // AliTPCclusterMI *c=clusters[noc-i];
6860 AliTPCclusterMI *c=clusters[i];
6862 if (lab == TMath::Abs(c->GetLabel(0)) ||
6863 lab == TMath::Abs(c->GetLabel(1)) ||
6864 lab == TMath::Abs(c->GetLabel(2))) max++;
6867 if (max < Int_t(0.5*tail)) lab=-lab;
6874 //delete[] clusters;
6878 //__________________________________________________________________________
6879 Int_t AliTPCtrackerMI::CookLabel(AliTPCseed *t, Float_t wrong,Int_t first, Int_t last) const {
6880 //--------------------------------------------------------------------
6881 //This function "cooks" a track label. If label<0, this track is fake.
6882 //--------------------------------------------------------------------
6883 Int_t noc=t->GetNumberOfClusters();
6885 //printf("\nnot founded prolongation\n\n\n");
6891 AliTPCclusterMI *clusters[160];
6893 for (Int_t i=0;i<160;i++) {
6900 for (i=0; i<160 && current<noc; i++) {
6901 if (i<first) continue;
6902 if (i>last) continue;
6903 Int_t index=t->GetClusterIndex2(i);
6904 if (index<=0) continue;
6905 if (index&0x8000) continue;
6907 //clusters[current]=GetClusterMI(index);
6908 if (t->GetClusterPointer(i)){
6909 clusters[current]=t->GetClusterPointer(i);
6914 if (noc<5) return -1;
6915 Int_t lab=123456789;
6916 for (i=0; i<noc; i++) {
6917 AliTPCclusterMI *c=clusters[i];
6919 lab=TMath::Abs(c->GetLabel(0));
6921 for (j=0; j<noc; j++) if (lb[j]==lab || mx[j]==0) break;
6927 for (i=0; i<noc; i++) if (mx[i]>max) {max=mx[i]; lab=lb[i];}
6929 for (i=0; i<noc; i++) {
6930 AliTPCclusterMI *c=clusters[i];
6932 if (TMath::Abs(c->GetLabel(1)) == lab ||
6933 TMath::Abs(c->GetLabel(2)) == lab ) max++;
6936 if ((1.- Float_t(max)/noc) > wrong) lab=-lab;
6939 Int_t tail=Int_t(0.10*noc);
6942 for (i=1; i<=160&&ind<tail; i++) {
6943 // AliTPCclusterMI *c=clusters[noc-i];
6944 AliTPCclusterMI *c=clusters[i];
6946 if (lab == TMath::Abs(c->GetLabel(0)) ||
6947 lab == TMath::Abs(c->GetLabel(1)) ||
6948 lab == TMath::Abs(c->GetLabel(2))) max++;
6951 if (max < Int_t(0.5*tail)) lab=-lab;
6954 // t->SetLabel(lab);
6958 //delete[] clusters;
6962 Int_t AliTPCtrackerMI::GetRowNumber(Double_t x[3]) const
6964 //return pad row number for given x vector
6965 Float_t phi = TMath::ATan2(x[1],x[0]);
6966 if(phi<0) phi=2.*TMath::Pi()+phi;
6967 // Get the local angle in the sector philoc
6968 const Float_t kRaddeg = 180/3.14159265358979312;
6969 Float_t phiangle = (Int_t (phi*kRaddeg/20.) + 0.5)*20./kRaddeg;
6970 Double_t localx = x[0]*TMath::Cos(phiangle)-x[1]*TMath::Sin(phiangle);
6971 return GetRowNumber(localx);
6976 void AliTPCtrackerMI::MakeBitmaps(AliTPCseed *t)
6978 //-----------------------------------------------------------------------
6979 // Fill the cluster and sharing bitmaps of the track
6980 //-----------------------------------------------------------------------
6982 Int_t firstpoint = 0;
6983 Int_t lastpoint = 159;
6984 AliTPCTrackerPoint *point;
6986 for (int iter=firstpoint; iter<lastpoint; iter++) {
6987 point = t->GetTrackPoint(iter);
6989 t->SetClusterMapBit(iter, kTRUE);
6990 if (point->IsShared())
6991 t->SetSharedMapBit(iter,kTRUE);
6993 t->SetSharedMapBit(iter, kFALSE);
6996 t->SetClusterMapBit(iter, kFALSE);
6997 t->SetSharedMapBit(iter, kFALSE);
7004 void AliTPCtrackerMI::AddCovariance(AliTPCseed * seed){
7006 // Adding systematic error
7007 // !!!! the systematic error for element 4 is in 1/cm not in pt
7009 const Double_t *param = AliTPCReconstructor::GetRecoParam()->GetSystematicError();
7011 for (Int_t i=0;i<15;i++) covar[i]=0;
7017 covar[0] = param[0]*param[0];
7018 covar[2] = param[1]*param[1];
7019 covar[5] = param[2]*param[2];
7020 covar[9] = param[3]*param[3];
7021 Double_t facC = AliTracker::GetBz()*kB2C;
7022 covar[14]= param[4]*param[4]*facC*facC;
7023 seed->AddCovariance(covar);