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];
321 if (AliTPCReconstructor::StreamLevel()>0) {
322 (*fDebugStreamer)<<"ErrParam"<<
331 "rmsy2p30="<<rmsy2p30<<
332 "rmsz2p30="<<rmsz2p30<<
333 "rmsy2p30R="<<rmsy2p30R<<
334 "rmsz2p30R="<<rmsz2p30R<<
339 if (rdistance2>16) return 3;
342 if ((rdistancey2>9. || rdistancez2>9.) && cluster->GetType()==0)
343 return 2; //suspisiouce - will be changed
345 if ((rdistancey2>6.25 || rdistancez2>6.25) && cluster->GetType()>0)
346 // strict cut on overlaped cluster
347 return 2; //suspisiouce - will be changed
349 if ( (rdistancey2>1. || rdistancez2>6.25 )
350 && cluster->GetType()<0){
351 seed->SetNFoundable(seed->GetNFoundable()-1);
360 //_____________________________________________________________________________
361 AliTPCtrackerMI::AliTPCtrackerMI(const AliTPCParam *par):
363 fkNIS(par->GetNInnerSector()/2),
365 fkNOS(par->GetNOuterSector()/2),
382 //---------------------------------------------------------------------
383 // The main TPC tracker constructor
384 //---------------------------------------------------------------------
385 fInnerSec=new AliTPCtrackerSector[fkNIS];
386 fOuterSec=new AliTPCtrackerSector[fkNOS];
389 for (i=0; i<fkNIS; i++) fInnerSec[i].Setup(par,0);
390 for (i=0; i<fkNOS; i++) fOuterSec[i].Setup(par,1);
393 Int_t nrowlow = par->GetNRowLow();
394 Int_t nrowup = par->GetNRowUp();
397 for (Int_t i=0;i<nrowlow;i++){
398 fXRow[i] = par->GetPadRowRadiiLow(i);
399 fPadLength[i]= par->GetPadPitchLength(0,i);
400 fYMax[i] = fXRow[i]*TMath::Tan(0.5*par->GetInnerAngle());
404 for (Int_t i=0;i<nrowup;i++){
405 fXRow[i+nrowlow] = par->GetPadRowRadiiUp(i);
406 fPadLength[i+nrowlow] = par->GetPadPitchLength(60,i);
407 fYMax[i+nrowlow] = fXRow[i+nrowlow]*TMath::Tan(0.5*par->GetOuterAngle());
410 if (AliTPCReconstructor::StreamLevel()>0) {
411 fDebugStreamer = new TTreeSRedirector("TPCdebug.root");
414 //________________________________________________________________________
415 AliTPCtrackerMI::AliTPCtrackerMI(const AliTPCtrackerMI &t):
436 //------------------------------------
437 // dummy copy constructor
438 //------------------------------------------------------------------
441 AliTPCtrackerMI & AliTPCtrackerMI::operator=(const AliTPCtrackerMI& /*r*/){
442 //------------------------------
444 //--------------------------------------------------------------
447 //_____________________________________________________________________________
448 AliTPCtrackerMI::~AliTPCtrackerMI() {
449 //------------------------------------------------------------------
450 // TPC tracker destructor
451 //------------------------------------------------------------------
458 if (fDebugStreamer) delete fDebugStreamer;
462 void AliTPCtrackerMI::FillESD(TObjArray* arr)
466 //fill esds using updated tracks
468 // write tracks to the event
469 // store index of the track
470 Int_t nseed=arr->GetEntriesFast();
471 //FindKinks(arr,fEvent);
472 for (Int_t i=0; i<nseed; i++) {
473 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
477 if (AliTPCReconstructor::StreamLevel()>1) {
478 (*fDebugStreamer)<<"Track0"<<
482 // pt->PropagateTo(fParam->GetInnerRadiusLow());
483 if (pt->GetKinkIndex(0)<=0){ //don't propagate daughter tracks
484 pt->PropagateTo(fParam->GetInnerRadiusLow());
487 if (( pt->GetPoints()[2]- pt->GetPoints()[0])>5 && pt->GetPoints()[3]>0.8){
489 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
490 iotrack.SetTPCPoints(pt->GetPoints());
491 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
492 iotrack.SetV0Indexes(pt->GetV0Indexes());
493 // iotrack.SetTPCpid(pt->fTPCr);
494 //iotrack.SetTPCindex(i);
495 fEvent->AddTrack(&iotrack);
499 if ( (pt->GetNumberOfClusters()>70)&& (Float_t(pt->GetNumberOfClusters())/Float_t(pt->GetNFoundable()))>0.55) {
501 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
502 iotrack.SetTPCPoints(pt->GetPoints());
503 //iotrack.SetTPCindex(i);
504 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
505 iotrack.SetV0Indexes(pt->GetV0Indexes());
506 // iotrack.SetTPCpid(pt->fTPCr);
507 fEvent->AddTrack(&iotrack);
511 // short tracks - maybe decays
513 if ( (pt->GetNumberOfClusters()>30) && (Float_t(pt->GetNumberOfClusters())/Float_t(pt->GetNFoundable()))>0.70) {
514 Int_t found,foundable,shared;
515 pt->GetClusterStatistic(0,60,found, foundable,shared,kFALSE);
516 if ( (found>20) && (pt->GetNShared()/float(pt->GetNumberOfClusters())<0.2)){
518 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
519 //iotrack.SetTPCindex(i);
520 iotrack.SetTPCPoints(pt->GetPoints());
521 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
522 iotrack.SetV0Indexes(pt->GetV0Indexes());
523 //iotrack.SetTPCpid(pt->fTPCr);
524 fEvent->AddTrack(&iotrack);
529 if ( (pt->GetNumberOfClusters()>20) && (Float_t(pt->GetNumberOfClusters())/Float_t(pt->GetNFoundable()))>0.8) {
530 Int_t found,foundable,shared;
531 pt->GetClusterStatistic(0,60,found, foundable,shared,kFALSE);
532 if (found<20) continue;
533 if (pt->GetNShared()/float(pt->GetNumberOfClusters())>0.2) continue;
536 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
537 iotrack.SetTPCPoints(pt->GetPoints());
538 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
539 iotrack.SetV0Indexes(pt->GetV0Indexes());
540 //iotrack.SetTPCpid(pt->fTPCr);
541 //iotrack.SetTPCindex(i);
542 fEvent->AddTrack(&iotrack);
545 // short tracks - secondaties
547 if ( (pt->GetNumberOfClusters()>30) ) {
548 Int_t found,foundable,shared;
549 pt->GetClusterStatistic(128,158,found, foundable,shared,kFALSE);
550 if ( (found>20) && (pt->GetNShared()/float(pt->GetNumberOfClusters())<0.2) &&float(found)/float(foundable)>0.8){
552 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
553 iotrack.SetTPCPoints(pt->GetPoints());
554 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
555 iotrack.SetV0Indexes(pt->GetV0Indexes());
556 //iotrack.SetTPCpid(pt->fTPCr);
557 //iotrack.SetTPCindex(i);
558 fEvent->AddTrack(&iotrack);
563 if ( (pt->GetNumberOfClusters()>15)) {
564 Int_t found,foundable,shared;
565 pt->GetClusterStatistic(138,158,found, foundable,shared,kFALSE);
566 if (found<15) continue;
567 if (foundable<=0) continue;
568 if (pt->GetNShared()/float(pt->GetNumberOfClusters())>0.2) continue;
569 if (float(found)/float(foundable)<0.8) continue;
572 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
573 iotrack.SetTPCPoints(pt->GetPoints());
574 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
575 iotrack.SetV0Indexes(pt->GetV0Indexes());
576 // iotrack.SetTPCpid(pt->fTPCr);
577 //iotrack.SetTPCindex(i);
578 fEvent->AddTrack(&iotrack);
583 printf("Number of filled ESDs-\t%d\n",fEvent->GetNumberOfTracks());
590 Double_t AliTPCtrackerMI::ErrY2(AliTPCseed* seed, AliTPCclusterMI * cl){
593 // Use calibrated cluster error from OCDB
595 AliTPCClusterParam * clparam = AliTPCcalibDB::Instance()->GetClusterParam();
597 Float_t z = TMath::Abs(fParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
598 Int_t ctype = cl->GetType();
599 Int_t type = (cl->GetRow()<63) ? 0: (cl->GetRow()>126) ? 1:2;
600 Double_t angle = seed->GetSnp()*seed->GetSnp();
601 angle = TMath::Sqrt(TMath::Abs(angle/(1.-angle)));
602 Double_t erry2 = clparam->GetError0Par(0,type, z,angle);
604 erry2+=0.5; // edge cluster
607 seed->SetErrorY2(erry2);
611 //calculate look-up table at the beginning
612 // static Bool_t ginit = kFALSE;
613 // static Float_t gnoise1,gnoise2,gnoise3;
614 // static Float_t ggg1[10000];
615 // static Float_t ggg2[10000];
616 // static Float_t ggg3[10000];
617 // static Float_t glandau1[10000];
618 // static Float_t glandau2[10000];
619 // static Float_t glandau3[10000];
621 // static Float_t gcor01[500];
622 // static Float_t gcor02[500];
623 // static Float_t gcorp[500];
627 // if (ginit==kFALSE){
628 // for (Int_t i=1;i<500;i++){
629 // Float_t rsigma = float(i)/100.;
630 // gcor02[i] = TMath::Max(0.78 +TMath::Exp(7.4*(rsigma-1.2)),0.6);
631 // gcor01[i] = TMath::Max(0.72 +TMath::Exp(3.36*(rsigma-1.2)),0.6);
632 // gcorp[i] = TMath::Max(TMath::Power((rsigma+0.5),1.5),1.2);
636 // for (Int_t i=3;i<10000;i++){
640 // Float_t amp = float(i);
641 // Float_t padlength =0.75;
642 // gnoise1 = 0.0004/padlength;
643 // Float_t nel = 0.268*amp;
644 // Float_t nprim = 0.155*amp;
645 // ggg1[i] = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.001*nel/(padlength*padlength))/nel;
646 // glandau1[i] = (2.+0.12*nprim)*0.5* (2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
647 // if (glandau1[i]>1) glandau1[i]=1;
648 // glandau1[i]*=padlength*padlength/12.;
652 // gnoise2 = 0.0004/padlength;
654 // nprim = 0.133*amp;
655 // ggg2[i] = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
656 // glandau2[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
657 // if (glandau2[i]>1) glandau2[i]=1;
658 // glandau2[i]*=padlength*padlength/12.;
663 // gnoise3 = 0.0004/padlength;
665 // nprim = 0.133*amp;
666 // ggg3[i] = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
667 // glandau3[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
668 // if (glandau3[i]>1) glandau3[i]=1;
669 // glandau3[i]*=padlength*padlength/12.;
677 // Int_t amp = int(TMath::Abs(cl->GetQ()));
679 // seed->SetErrorY2(1.);
683 // Float_t z = TMath::Abs(fParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
684 // Int_t ctype = cl->GetType();
685 // Float_t padlength= GetPadPitchLength(seed->GetRow());
686 // Double_t angle2 = seed->GetSnp()*seed->GetSnp();
687 // angle2 = angle2/(1-angle2);
689 // //cluster "quality"
690 // Int_t rsigmay = int(100.*cl->GetSigmaY2()/(seed->GetCurrentSigmaY2()));
693 // if (fSectors==fInnerSec){
694 // snoise2 = gnoise1;
695 // res = ggg1[amp]*z+glandau1[amp]*angle2;
696 // if (ctype==0) res *= gcor01[rsigmay];
699 // res*= gcorp[rsigmay];
703 // if (padlength<1.1){
704 // snoise2 = gnoise2;
705 // res = ggg2[amp]*z+glandau2[amp]*angle2;
706 // if (ctype==0) res *= gcor02[rsigmay];
709 // res*= gcorp[rsigmay];
713 // snoise2 = gnoise3;
714 // res = ggg3[amp]*z+glandau3[amp]*angle2;
715 // if (ctype==0) res *= gcor02[rsigmay];
718 // res*= gcorp[rsigmay];
725 // res*=2.4; // overestimate error 2 times
729 // if (res<2*snoise2)
732 // seed->SetErrorY2(res);
740 Double_t AliTPCtrackerMI::ErrZ2(AliTPCseed* seed, AliTPCclusterMI * cl){
743 // Use calibrated cluster error from OCDB
745 AliTPCClusterParam * clparam = AliTPCcalibDB::Instance()->GetClusterParam();
747 Float_t z = TMath::Abs(fParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
748 Int_t ctype = cl->GetType();
749 Int_t type = (cl->GetRow()<63) ? 0: (cl->GetRow()>126) ? 1:2;
751 Double_t angle2 = seed->GetSnp()*seed->GetSnp();
752 angle2 = seed->GetTgl()*seed->GetTgl()*(1+angle2/(1-angle2));
753 Double_t angle = TMath::Sqrt(TMath::Abs(angle2));
754 Double_t errz2 = clparam->GetError0Par(1,type, z,angle);
756 errz2+=0.5; // edge cluster
759 seed->SetErrorZ2(errz2);
765 // //seed->SetErrorY2(0.1);
767 // //calculate look-up table at the beginning
768 // static Bool_t ginit = kFALSE;
769 // static Float_t gnoise1,gnoise2,gnoise3;
770 // static Float_t ggg1[10000];
771 // static Float_t ggg2[10000];
772 // static Float_t ggg3[10000];
773 // static Float_t glandau1[10000];
774 // static Float_t glandau2[10000];
775 // static Float_t glandau3[10000];
777 // static Float_t gcor01[1000];
778 // static Float_t gcor02[1000];
779 // static Float_t gcorp[1000];
783 // if (ginit==kFALSE){
784 // for (Int_t i=1;i<1000;i++){
785 // Float_t rsigma = float(i)/100.;
786 // gcor02[i] = TMath::Max(0.81 +TMath::Exp(6.8*(rsigma-1.2)),0.6);
787 // gcor01[i] = TMath::Max(0.72 +TMath::Exp(2.04*(rsigma-1.2)),0.6);
788 // gcorp[i] = TMath::Max(TMath::Power((rsigma+0.5),1.5),1.2);
792 // for (Int_t i=3;i<10000;i++){
796 // Float_t amp = float(i);
797 // Float_t padlength =0.75;
798 // gnoise1 = 0.0004/padlength;
799 // Float_t nel = 0.268*amp;
800 // Float_t nprim = 0.155*amp;
801 // ggg1[i] = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.001*nel/(padlength*padlength))/nel;
802 // glandau1[i] = (2.+0.12*nprim)*0.5* (2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
803 // if (glandau1[i]>1) glandau1[i]=1;
804 // glandau1[i]*=padlength*padlength/12.;
808 // gnoise2 = 0.0004/padlength;
810 // nprim = 0.133*amp;
811 // ggg2[i] = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
812 // glandau2[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
813 // if (glandau2[i]>1) glandau2[i]=1;
814 // glandau2[i]*=padlength*padlength/12.;
819 // gnoise3 = 0.0004/padlength;
821 // nprim = 0.133*amp;
822 // ggg3[i] = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
823 // glandau3[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
824 // if (glandau3[i]>1) glandau3[i]=1;
825 // glandau3[i]*=padlength*padlength/12.;
833 // Int_t amp = int(TMath::Abs(cl->GetQ()));
835 // seed->SetErrorY2(1.);
839 // Float_t z = TMath::Abs(fParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
840 // Int_t ctype = cl->GetType();
841 // Float_t padlength= GetPadPitchLength(seed->GetRow());
843 // Double_t angle2 = seed->GetSnp()*seed->GetSnp();
844 // // if (angle2<0.6) angle2 = 0.6;
845 // angle2 = seed->GetTgl()*seed->GetTgl()*(1+angle2/(1-angle2));
847 // //cluster "quality"
848 // Int_t rsigmaz = int(100.*cl->GetSigmaZ2()/(seed->GetCurrentSigmaZ2()));
851 // if (fSectors==fInnerSec){
852 // snoise2 = gnoise1;
853 // res = ggg1[amp]*z+glandau1[amp]*angle2;
854 // if (ctype==0) res *= gcor01[rsigmaz];
857 // res*= gcorp[rsigmaz];
861 // if (padlength<1.1){
862 // snoise2 = gnoise2;
863 // res = ggg2[amp]*z+glandau2[amp]*angle2;
864 // if (ctype==0) res *= gcor02[rsigmaz];
867 // res*= gcorp[rsigmaz];
871 // snoise2 = gnoise3;
872 // res = ggg3[amp]*z+glandau3[amp]*angle2;
873 // if (ctype==0) res *= gcor02[rsigmaz];
876 // res*= gcorp[rsigmaz];
885 // if ((ctype<0) &&<70){
890 // if (res<2*snoise2)
892 // if (res>3) res =3;
893 // seed->SetErrorZ2(res);
901 void AliTPCtrackerMI::RotateToLocal(AliTPCseed *seed)
903 //rotate to track "local coordinata
904 Float_t x = seed->GetX();
905 Float_t y = seed->GetY();
906 Float_t ymax = x*TMath::Tan(0.5*fSectors->GetAlpha());
909 seed->SetRelativeSector((seed->GetRelativeSector()+1) % fN);
910 if (!seed->Rotate(fSectors->GetAlpha()))
912 } else if (y <-ymax) {
913 seed->SetRelativeSector((seed->GetRelativeSector()-1+fN) % fN);
914 if (!seed->Rotate(-fSectors->GetAlpha()))
922 //_____________________________________________________________________________
923 Double_t AliTPCtrackerMI::F1old(Double_t x1,Double_t y1,
924 Double_t x2,Double_t y2,
925 Double_t x3,Double_t y3)
927 //-----------------------------------------------------------------
928 // Initial approximation of the track curvature
929 //-----------------------------------------------------------------
930 Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
931 Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
932 (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
933 Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
934 (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
936 Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
937 if ( xr*xr+yr*yr<=0.00000000000001) return 100;
938 return -xr*yr/sqrt(xr*xr+yr*yr);
943 //_____________________________________________________________________________
944 Double_t AliTPCtrackerMI::F1(Double_t x1,Double_t y1,
945 Double_t x2,Double_t y2,
946 Double_t x3,Double_t y3)
948 //-----------------------------------------------------------------
949 // Initial approximation of the track curvature
950 //-----------------------------------------------------------------
956 Double_t det = x3*y2-x2*y3;
961 Double_t u = 0.5* (x2*(x2-x3)+y2*(y2-y3))/det;
962 Double_t x0 = x3*0.5-y3*u;
963 Double_t y0 = y3*0.5+x3*u;
964 Double_t c2 = 1/TMath::Sqrt(x0*x0+y0*y0);
970 Double_t AliTPCtrackerMI::F2(Double_t x1,Double_t y1,
971 Double_t x2,Double_t y2,
972 Double_t x3,Double_t y3)
974 //-----------------------------------------------------------------
975 // Initial approximation of the track curvature
976 //-----------------------------------------------------------------
982 Double_t det = x3*y2-x2*y3;
987 Double_t u = 0.5* (x2*(x2-x3)+y2*(y2-y3))/det;
988 Double_t x0 = x3*0.5-y3*u;
989 Double_t y0 = y3*0.5+x3*u;
990 Double_t c2 = 1/TMath::Sqrt(x0*x0+y0*y0);
999 //_____________________________________________________________________________
1000 Double_t AliTPCtrackerMI::F2old(Double_t x1,Double_t y1,
1001 Double_t x2,Double_t y2,
1002 Double_t x3,Double_t y3)
1004 //-----------------------------------------------------------------
1005 // Initial approximation of the track curvature times center of curvature
1006 //-----------------------------------------------------------------
1007 Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
1008 Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
1009 (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
1010 Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
1011 (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
1013 Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
1015 return -a/(d*y1-b)*xr/sqrt(xr*xr+yr*yr);
1018 //_____________________________________________________________________________
1019 Double_t AliTPCtrackerMI::F3(Double_t x1,Double_t y1,
1020 Double_t x2,Double_t y2,
1021 Double_t z1,Double_t z2)
1023 //-----------------------------------------------------------------
1024 // Initial approximation of the tangent of the track dip angle
1025 //-----------------------------------------------------------------
1026 return (z1 - z2)/sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
1030 Double_t AliTPCtrackerMI::F3n(Double_t x1,Double_t y1,
1031 Double_t x2,Double_t y2,
1032 Double_t z1,Double_t z2, Double_t c)
1034 //-----------------------------------------------------------------
1035 // Initial approximation of the tangent of the track dip angle
1036 //-----------------------------------------------------------------
1040 //angle1 = (z1-z2)*c/(TMath::ASin(c*x1-ni)-TMath::ASin(c*x2-ni));
1042 Double_t d = TMath::Sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
1043 if (TMath::Abs(d*c*0.5)>1) return 0;
1044 // Double_t angle2 = TMath::ASin(d*c*0.5);
1045 // Double_t angle2 = AliTPCFastMath::FastAsin(d*c*0.5);
1046 Double_t angle2 = (d*c*0.5>0.1)? TMath::ASin(d*c*0.5): AliTPCFastMath::FastAsin(d*c*0.5);
1048 angle2 = (z1-z2)*c/(angle2*2.);
1052 Bool_t AliTPCtrackerMI::GetProlongation(Double_t x1, Double_t x2, Double_t x[5], Double_t &y, Double_t &z)
1053 {//-----------------------------------------------------------------
1054 // This function find proloncation of a track to a reference plane x=x2.
1055 //-----------------------------------------------------------------
1059 if (TMath::Abs(x[4]*x1 - x[2]) >= 0.999) {
1063 Double_t c1=x[4]*x1 - x[2], r1=sqrt(1.- c1*c1);
1064 Double_t c2=x[4]*x2 - x[2], r2=sqrt(1.- c2*c2);
1068 Double_t dy = dx*(c1+c2)/(r1+r2);
1071 Double_t delta = x[4]*dx*(c1+c2)/(c1*r2 + c2*r1);
1073 if (TMath::Abs(delta)>0.01){
1074 dz = x[3]*TMath::ASin(delta)/x[4];
1076 dz = x[3]*AliTPCFastMath::FastAsin(delta)/x[4];
1079 //dz = x[3]*AliTPCFastMath::FastAsin(delta)/x[4];
1087 Int_t AliTPCtrackerMI::LoadClusters (TTree *tree)
1092 return LoadClusters();
1096 Int_t AliTPCtrackerMI::LoadClusters(TObjArray *arr)
1099 // load clusters to the memory
1100 AliTPCClustersRow *clrow = 0x0;
1101 Int_t lower = arr->LowerBound();
1102 Int_t entries = arr->GetEntriesFast();
1104 for (Int_t i=lower; i<entries; i++) {
1105 clrow = (AliTPCClustersRow*) arr->At(i);
1106 if(!clrow) continue;
1107 if(!clrow->GetArray()) continue;
1111 fParam->AdjustSectorRow(clrow->GetID(),sec,row);
1113 for (Int_t icl=0; icl<clrow->GetArray()->GetEntriesFast(); icl++){
1114 Transform((AliTPCclusterMI*)(clrow->GetArray()->At(icl)));
1117 if (clrow->GetArray()->GetEntriesFast()<=0) continue;
1118 AliTPCtrackerRow * tpcrow=0;
1121 tpcrow = &(fInnerSec[sec%fkNIS][row]);
1125 tpcrow = &(fOuterSec[(sec-fkNIS*2)%fkNOS][row]);
1126 left = (sec-fkNIS*2)/fkNOS;
1129 tpcrow->SetN1(clrow->GetArray()->GetEntriesFast());
1130 tpcrow->SetClusters1(new AliTPCclusterMI[tpcrow->GetN1()]);
1131 for (Int_t i=0;i<tpcrow->GetN1();i++)
1132 tpcrow->SetCluster1(i, *(AliTPCclusterMI*)(clrow->GetArray()->At(i)));
1135 tpcrow->SetN2(clrow->GetArray()->GetEntriesFast());
1136 tpcrow->SetClusters2(new AliTPCclusterMI[tpcrow->GetN2()]);
1137 for (Int_t i=0;i<tpcrow->GetN2();i++)
1138 tpcrow->SetCluster2(i,*(AliTPCclusterMI*)(clrow->GetArray()->At(i)));
1148 Int_t AliTPCtrackerMI::LoadClusters(TClonesArray *arr)
1151 // load clusters to the memory from one
1154 AliTPCclusterMI *clust=0;
1155 Int_t count[72][96] = { {0} , {0} };
1157 // loop over clusters
1158 for (Int_t icl=0; icl<arr->GetEntriesFast(); icl++) {
1159 clust = (AliTPCclusterMI*)arr->At(icl);
1160 if(!clust) continue;
1161 //printf("cluster: det %d, row %d \n", clust->GetDetector(),clust->GetRow());
1163 // transform clusters
1166 // count clusters per pad row
1167 count[clust->GetDetector()][clust->GetRow()]++;
1170 // insert clusters to sectors
1171 for (Int_t icl=0; icl<arr->GetEntriesFast(); icl++) {
1172 clust = (AliTPCclusterMI*)arr->At(icl);
1173 if(!clust) continue;
1175 Int_t sec = clust->GetDetector();
1176 Int_t row = clust->GetRow();
1178 // filter overlapping pad rows needed by HLT
1179 if(sec<fkNIS*2) { //IROCs
1180 if(row == 30) continue;
1183 if(row == 27 || row == 76) continue;
1189 fInnerSec[sec%fkNIS].InsertCluster(clust, count[sec][row], fParam);
1192 left = (sec-fkNIS*2)/fkNOS;
1193 fOuterSec[(sec-fkNIS*2)%fkNOS].InsertCluster(clust, count[sec][row], fParam);
1197 // Load functions must be called behind LoadCluster(TClonesArray*)
1199 //LoadOuterSectors();
1200 //LoadInnerSectors();
1206 Int_t AliTPCtrackerMI::LoadClusters()
1209 // load clusters to the memory
1210 AliTPCClustersRow *clrow= new AliTPCClustersRow;
1211 clrow->SetClass("AliTPCclusterMI");
1213 clrow->GetArray()->ExpandCreateFast(10000);
1215 // TTree * tree = fClustersArray.GetTree();
1217 TTree * tree = fInput;
1218 TBranch * br = tree->GetBranch("Segment");
1219 br->SetAddress(&clrow);
1221 Int_t j=Int_t(tree->GetEntries());
1222 for (Int_t i=0; i<j; i++) {
1226 fParam->AdjustSectorRow(clrow->GetID(),sec,row);
1227 for (Int_t icl=0; icl<clrow->GetArray()->GetEntriesFast(); icl++){
1228 Transform((AliTPCclusterMI*)(clrow->GetArray()->At(icl)));
1231 AliTPCtrackerRow * tpcrow=0;
1234 tpcrow = &(fInnerSec[sec%fkNIS][row]);
1238 tpcrow = &(fOuterSec[(sec-fkNIS*2)%fkNOS][row]);
1239 left = (sec-fkNIS*2)/fkNOS;
1242 tpcrow->SetN1(clrow->GetArray()->GetEntriesFast());
1243 tpcrow->SetClusters1(new AliTPCclusterMI[tpcrow->GetN1()]);
1244 for (Int_t i=0;i<tpcrow->GetN1();i++)
1245 tpcrow->SetCluster1(i, *(AliTPCclusterMI*)(clrow->GetArray()->At(i)));
1248 tpcrow->SetN2(clrow->GetArray()->GetEntriesFast());
1249 tpcrow->SetClusters2(new AliTPCclusterMI[tpcrow->GetN2()]);
1250 for (Int_t i=0;i<tpcrow->GetN2();i++)
1251 tpcrow->SetCluster2(i,*(AliTPCclusterMI*)(clrow->GetArray()->At(i)));
1262 void AliTPCtrackerMI::UnloadClusters()
1265 // unload clusters from the memory
1267 Int_t nrows = fOuterSec->GetNRows();
1268 for (Int_t sec = 0;sec<fkNOS;sec++)
1269 for (Int_t row = 0;row<nrows;row++){
1270 AliTPCtrackerRow* tpcrow = &(fOuterSec[sec%fkNOS][row]);
1272 // if (tpcrow->fClusters1) delete []tpcrow->fClusters1;
1273 // if (tpcrow->fClusters2) delete []tpcrow->fClusters2;
1275 tpcrow->ResetClusters();
1278 nrows = fInnerSec->GetNRows();
1279 for (Int_t sec = 0;sec<fkNIS;sec++)
1280 for (Int_t row = 0;row<nrows;row++){
1281 AliTPCtrackerRow* tpcrow = &(fInnerSec[sec%fkNIS][row]);
1283 // if (tpcrow->fClusters1) delete []tpcrow->fClusters1;
1284 //if (tpcrow->fClusters2) delete []tpcrow->fClusters2;
1286 tpcrow->ResetClusters();
1292 void AliTPCtrackerMI::FillClusterArray(TObjArray* array) const{
1294 // Filling cluster to the array - For visualization purposes
1297 nrows = fOuterSec->GetNRows();
1298 for (Int_t sec = 0;sec<fkNOS;sec++)
1299 for (Int_t row = 0;row<nrows;row++){
1300 AliTPCtrackerRow* tpcrow = &(fOuterSec[sec%fkNOS][row]);
1301 if (!tpcrow) continue;
1302 for (Int_t icl = 0;icl<tpcrow->GetN();icl++){
1303 array->AddLast((TObject*)((*tpcrow)[icl]));
1306 nrows = fInnerSec->GetNRows();
1307 for (Int_t sec = 0;sec<fkNIS;sec++)
1308 for (Int_t row = 0;row<nrows;row++){
1309 AliTPCtrackerRow* tpcrow = &(fInnerSec[sec%fkNIS][row]);
1310 if (!tpcrow) continue;
1311 for (Int_t icl = 0;icl<tpcrow->GetN();icl++){
1312 array->AddLast((TObject*)(*tpcrow)[icl]);
1318 void AliTPCtrackerMI::Transform(AliTPCclusterMI * cluster){
1322 AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ;
1324 AliFatal("Tranformations not in calibDB");
1326 Double_t x[3]={cluster->GetRow(),cluster->GetPad(),cluster->GetTimeBin()};
1327 Int_t i[1]={cluster->GetDetector()};
1328 transform->Transform(x,i,0,1);
1329 if (cluster->GetDetector()%36>17){
1334 // in debug mode check the transformation
1336 if (AliTPCReconstructor::StreamLevel()>1) {
1338 cluster->GetGlobalXYZ(gx);
1339 Int_t event = (fEvent==NULL)? 0: fEvent->GetEventNumberInFile();
1340 TTreeSRedirector &cstream = *fDebugStreamer;
1341 cstream<<"Transform"<<
1352 cluster->SetX(x[0]);
1353 cluster->SetY(x[1]);
1354 cluster->SetZ(x[2]);
1360 //if (!fParam->IsGeoRead()) fParam->ReadGeoMatrices();
1361 TGeoHMatrix *mat = fParam->GetClusterMatrix(cluster->GetDetector());
1363 Double_t pos[3]= {cluster->GetX(),cluster->GetY(),cluster->GetZ()};
1364 Double_t posC[3]={cluster->GetX(),cluster->GetY(),cluster->GetZ()};
1365 if (mat) mat->LocalToMaster(pos,posC);
1367 // chack Loading of Geo matrices from GeoManager - TEMPORARY FIX
1369 cluster->SetX(posC[0]);
1370 cluster->SetY(posC[1]);
1371 cluster->SetZ(posC[2]);
1374 //_____________________________________________________________________________
1375 Int_t AliTPCtrackerMI::LoadOuterSectors() {
1376 //-----------------------------------------------------------------
1377 // This function fills outer TPC sectors with clusters.
1378 //-----------------------------------------------------------------
1379 Int_t nrows = fOuterSec->GetNRows();
1381 for (Int_t sec = 0;sec<fkNOS;sec++)
1382 for (Int_t row = 0;row<nrows;row++){
1383 AliTPCtrackerRow* tpcrow = &(fOuterSec[sec%fkNOS][row]);
1384 Int_t sec2 = sec+2*fkNIS;
1386 Int_t ncl = tpcrow->GetN1();
1388 AliTPCclusterMI *c= (tpcrow->GetCluster1(ncl));
1389 index=(((sec2<<8)+row)<<16)+ncl;
1390 tpcrow->InsertCluster(c,index);
1393 ncl = tpcrow->GetN2();
1395 AliTPCclusterMI *c= (tpcrow->GetCluster2(ncl));
1396 index=((((sec2+fkNOS)<<8)+row)<<16)+ncl;
1397 tpcrow->InsertCluster(c,index);
1400 // write indexes for fast acces
1402 for (Int_t i=0;i<510;i++)
1403 tpcrow->SetFastCluster(i,-1);
1404 for (Int_t i=0;i<tpcrow->GetN();i++){
1405 Int_t zi = Int_t((*tpcrow)[i]->GetZ()+255.);
1406 tpcrow->SetFastCluster(zi,i); // write index
1409 for (Int_t i=0;i<510;i++){
1410 if (tpcrow->GetFastCluster(i)<0)
1411 tpcrow->SetFastCluster(i,last);
1413 last = tpcrow->GetFastCluster(i);
1422 //_____________________________________________________________________________
1423 Int_t AliTPCtrackerMI::LoadInnerSectors() {
1424 //-----------------------------------------------------------------
1425 // This function fills inner TPC sectors with clusters.
1426 //-----------------------------------------------------------------
1427 Int_t nrows = fInnerSec->GetNRows();
1429 for (Int_t sec = 0;sec<fkNIS;sec++)
1430 for (Int_t row = 0;row<nrows;row++){
1431 AliTPCtrackerRow* tpcrow = &(fInnerSec[sec%fkNIS][row]);
1434 Int_t ncl = tpcrow->GetN1();
1436 AliTPCclusterMI *c= (tpcrow->GetCluster1(ncl));
1437 index=(((sec<<8)+row)<<16)+ncl;
1438 tpcrow->InsertCluster(c,index);
1441 ncl = tpcrow->GetN2();
1443 AliTPCclusterMI *c= (tpcrow->GetCluster2(ncl));
1444 index=((((sec+fkNIS)<<8)+row)<<16)+ncl;
1445 tpcrow->InsertCluster(c,index);
1448 // write indexes for fast acces
1450 for (Int_t i=0;i<510;i++)
1451 tpcrow->SetFastCluster(i,-1);
1452 for (Int_t i=0;i<tpcrow->GetN();i++){
1453 Int_t zi = Int_t((*tpcrow)[i]->GetZ()+255.);
1454 tpcrow->SetFastCluster(zi,i); // write index
1457 for (Int_t i=0;i<510;i++){
1458 if (tpcrow->GetFastCluster(i)<0)
1459 tpcrow->SetFastCluster(i,last);
1461 last = tpcrow->GetFastCluster(i);
1473 //_________________________________________________________________________
1474 AliTPCclusterMI *AliTPCtrackerMI::GetClusterMI(Int_t index) const {
1475 //--------------------------------------------------------------------
1476 // Return pointer to a given cluster
1477 //--------------------------------------------------------------------
1478 if (index<0) return 0; // no cluster
1479 Int_t sec=(index&0xff000000)>>24;
1480 Int_t row=(index&0x00ff0000)>>16;
1481 Int_t ncl=(index&0x00007fff)>>00;
1483 const AliTPCtrackerRow * tpcrow=0;
1484 AliTPCclusterMI * clrow =0;
1486 if (sec<0 || sec>=fkNIS*4) {
1487 AliWarning(Form("Wrong sector %d",sec));
1492 tpcrow = &(fInnerSec[sec%fkNIS][row]);
1493 if (tpcrow==0) return 0;
1496 if (tpcrow->GetN1()<=ncl) return 0;
1497 clrow = tpcrow->GetClusters1();
1500 if (tpcrow->GetN2()<=ncl) return 0;
1501 clrow = tpcrow->GetClusters2();
1505 tpcrow = &(fOuterSec[(sec-fkNIS*2)%fkNOS][row]);
1506 if (tpcrow==0) return 0;
1508 if (sec-2*fkNIS<fkNOS) {
1509 if (tpcrow->GetN1()<=ncl) return 0;
1510 clrow = tpcrow->GetClusters1();
1513 if (tpcrow->GetN2()<=ncl) return 0;
1514 clrow = tpcrow->GetClusters2();
1518 return &(clrow[ncl]);
1524 Int_t AliTPCtrackerMI::FollowToNext(AliTPCseed& t, Int_t nr) {
1525 //-----------------------------------------------------------------
1526 // This function tries to find a track prolongation to next pad row
1527 //-----------------------------------------------------------------
1529 Double_t x= GetXrow(nr), ymax=GetMaxY(nr);
1530 AliTPCclusterMI *cl=0;
1531 Int_t tpcindex= t.GetClusterIndex2(nr);
1533 // update current shape info every 5 pad-row
1534 // if ( (nr%5==0) || t.GetNumberOfClusters()<2 || (t.fCurrentSigmaY2<0.0001) ){
1538 if (fIteration>0 && tpcindex>=-1){ //if we have already clusters
1540 if (tpcindex==-1) return 0; //track in dead zone
1542 cl = t.GetClusterPointer(nr);
1543 if ( (cl==0) ) cl = GetClusterMI(tpcindex);
1544 t.SetCurrentClusterIndex1(tpcindex);
1547 Int_t relativesector = ((tpcindex&0xff000000)>>24)%18; // if previously accepted cluster in different sector
1548 Float_t angle = relativesector*fSectors->GetAlpha()+fSectors->GetAlphaShift();
1550 if (angle<-TMath::Pi()) angle += 2*TMath::Pi();
1551 if (angle>=TMath::Pi()) angle -= 2*TMath::Pi();
1553 if (TMath::Abs(angle-t.GetAlpha())>0.001){
1554 Double_t rotation = angle-t.GetAlpha();
1555 t.SetRelativeSector(relativesector);
1556 if (!t.Rotate(rotation)) return 0;
1558 if (!t.PropagateTo(x)) return 0;
1560 t.SetCurrentCluster(cl);
1562 Int_t accept = AcceptCluster(&t,t.GetCurrentCluster());
1563 if ((tpcindex&0x8000)==0) accept =0;
1565 //if founded cluster is acceptible
1566 if (cl->IsUsed(11)) { // id cluster is shared inrease uncertainty
1567 t.SetErrorY2(t.GetErrorY2()+0.03);
1568 t.SetErrorZ2(t.GetErrorZ2()+0.03);
1569 t.SetErrorY2(t.GetErrorY2()*3);
1570 t.SetErrorZ2(t.GetErrorZ2()*3);
1572 t.SetNFoundable(t.GetNFoundable()+1);
1573 UpdateTrack(&t,accept);
1578 if (TMath::Abs(t.GetSnp())>AliTPCReconstructor::GetMaxSnpTracker()) return 0; // cut on angle
1580 // not look for new cluster during refitting
1581 t.SetNFoundable(t.GetNFoundable()+1);
1586 // if (TMath::Abs(t.GetSnp())>0.95 || TMath::Abs(x*t.GetC()-t.GetEta())>0.95) return 0;// patch 28 fev 06
1587 Double_t y=t.GetYat(x);
1588 if (TMath::Abs(y)>ymax){
1590 t.SetRelativeSector((t.GetRelativeSector()+1) % fN);
1591 if (!t.Rotate(fSectors->GetAlpha()))
1593 } else if (y <-ymax) {
1594 t.SetRelativeSector((t.GetRelativeSector()-1+fN) % fN);
1595 if (!t.Rotate(-fSectors->GetAlpha()))
1601 if (!t.PropagateTo(x)) {
1602 if (fIteration==0) t.SetRemoval(10);
1606 Double_t z=t.GetZ();
1609 if (!IsActive(t.GetRelativeSector(),nr)) {
1611 t.SetClusterIndex2(nr,-1);
1614 //AliInfo(Form("A - Sector%d phi %f - alpha %f", t.fRelativeSector,y/x, t.GetAlpha()));
1615 Bool_t isActive = IsActive(t.GetRelativeSector(),nr);
1616 Bool_t isActive2 = (nr>=fInnerSec->GetNRows()) ? fOuterSec[t.GetRelativeSector()][nr-fInnerSec->GetNRows()].GetN()>0:fInnerSec[t.GetRelativeSector()][nr].GetN()>0;
1618 if (!isActive || !isActive2) return 0;
1620 const AliTPCtrackerRow &krow=GetRow(t.GetRelativeSector(),nr);
1621 if ( (t.GetSigmaY2()<0) || t.GetSigmaZ2()<0) return 0;
1623 Double_t roadz = 1.;
1625 if (TMath::Abs(TMath::Abs(y)-ymax)<krow.GetDeadZone()){
1627 t.SetClusterIndex2(nr,-1);
1632 if (TMath::Abs(z)<(AliTPCReconstructor::GetCtgRange()*x+10) && TMath::Abs(z)<fParam->GetZLength(0) && (TMath::Abs(t.GetSnp())<AliTPCReconstructor::GetMaxSnpTracker()))
1633 t.SetNFoundable(t.GetNFoundable()+1);
1639 // cl = krow.FindNearest2(y+10.,z,roady,roadz,index);
1640 cl = krow.FindNearest2(y,z,roady,roadz,index);
1641 if (cl) t.SetCurrentClusterIndex1(krow.GetIndex(index));
1644 t.SetCurrentCluster(cl);
1646 if (fIteration==2&&cl->IsUsed(10)) return 0;
1647 Int_t accept = AcceptCluster(&t,t.GetCurrentCluster());
1648 if (fIteration==2&&cl->IsUsed(11)) {
1649 t.SetErrorY2(t.GetErrorY2()+0.03);
1650 t.SetErrorZ2(t.GetErrorZ2()+0.03);
1651 t.SetErrorY2(t.GetErrorY2()*3);
1652 t.SetErrorZ2(t.GetErrorZ2()*3);
1655 if (t.fCurrentCluster->IsUsed(10)){
1660 if (t.fNShared>0.7*t.GetNumberOfClusters()) {
1666 if (accept<3) UpdateTrack(&t,accept);
1669 if ( fIteration==0 && t.GetNFoundable()*0.5 > t.GetNumberOfClusters()) t.SetRemoval(10);
1675 Int_t AliTPCtrackerMI::FollowToNextFast(AliTPCseed& t, Int_t nr) {
1676 //-----------------------------------------------------------------
1677 // This function tries to find a track prolongation to next pad row
1678 //-----------------------------------------------------------------
1680 Double_t x= GetXrow(nr), ymax=GetMaxY(nr);
1682 if (!t.GetProlongation(x,y,z)) {
1688 if (TMath::Abs(y)>ymax){
1691 t.SetRelativeSector((t.GetRelativeSector()+1) % fN);
1692 if (!t.Rotate(fSectors->GetAlpha()))
1694 } else if (y <-ymax) {
1695 t.SetRelativeSector((t.GetRelativeSector()-1+fN) % fN);
1696 if (!t.Rotate(-fSectors->GetAlpha()))
1699 if (!t.PropagateTo(x)) {
1702 t.GetProlongation(x,y,z);
1705 // update current shape info every 2 pad-row
1706 if ( (nr%2==0) || t.GetNumberOfClusters()<2 || (t.GetCurrentSigmaY2()<0.0001) ){
1707 // t.fCurrentSigmaY = GetSigmaY(&t);
1708 //t.fCurrentSigmaZ = GetSigmaZ(&t);
1712 AliTPCclusterMI *cl=0;
1717 const AliTPCtrackerRow &krow=GetRow(t.GetRelativeSector(),nr);
1718 if ( (t.GetSigmaY2()<0) || t.GetSigmaZ2()<0) return 0;
1720 Double_t roadz = 1.;
1723 if (TMath::Abs(TMath::Abs(y)-ymax)<krow.GetDeadZone()){
1725 t.SetClusterIndex2(row,-1);
1730 if (TMath::Abs(z)>(AliTPCReconstructor::GetCtgRange()*x+10)) t.SetClusterIndex2(row,-1);
1734 if ((cl==0)&&(krow)) {
1735 // cl = krow.FindNearest2(y+10,z,roady,roadz,index);
1736 cl = krow.FindNearest2(y,z,roady,roadz,index);
1738 if (cl) t.SetCurrentClusterIndex1(krow.GetIndex(index));
1742 t.SetCurrentCluster(cl);
1743 // Int_t accept = AcceptCluster(&t,t.fCurrentCluster);
1745 t.SetClusterIndex2(row,index);
1746 t.SetClusterPointer(row, cl);
1754 //_________________________________________________________________________
1755 Bool_t AliTPCtrackerMI::GetTrackPoint(Int_t index, AliTrackPoint &p ) const
1757 // Get track space point by index
1758 // return false in case the cluster doesn't exist
1759 AliTPCclusterMI *cl = GetClusterMI(index);
1760 if (!cl) return kFALSE;
1761 Int_t sector = (index&0xff000000)>>24;
1762 // Int_t row = (index&0x00ff0000)>>16;
1764 // xyz[0] = fParam->GetPadRowRadii(sector,row);
1765 xyz[0] = cl->GetX();
1766 xyz[1] = cl->GetY();
1767 xyz[2] = cl->GetZ();
1769 fParam->AdjustCosSin(sector,cos,sin);
1770 Float_t x = cos*xyz[0]-sin*xyz[1];
1771 Float_t y = cos*xyz[1]+sin*xyz[0];
1773 Float_t sigmaY2 = 0.027*cl->GetSigmaY2();
1774 if (sector < fParam->GetNInnerSector()) sigmaY2 *= 2.07;
1775 Float_t sigmaZ2 = 0.066*cl->GetSigmaZ2();
1776 if (sector < fParam->GetNInnerSector()) sigmaZ2 *= 1.77;
1777 cov[0] = sin*sin*sigmaY2;
1778 cov[1] = -sin*cos*sigmaY2;
1780 cov[3] = cos*cos*sigmaY2;
1783 p.SetXYZ(x,y,xyz[2],cov);
1784 AliGeomManager::ELayerID iLayer;
1786 if (sector < fParam->GetNInnerSector()) {
1787 iLayer = AliGeomManager::kTPC1;
1791 iLayer = AliGeomManager::kTPC2;
1792 idet = sector - fParam->GetNInnerSector();
1794 UShort_t volid = AliGeomManager::LayerToVolUID(iLayer,idet);
1795 p.SetVolumeID(volid);
1801 Int_t AliTPCtrackerMI::UpdateClusters(AliTPCseed& t, Int_t nr) {
1802 //-----------------------------------------------------------------
1803 // This function tries to find a track prolongation to next pad row
1804 //-----------------------------------------------------------------
1805 t.SetCurrentCluster(0);
1806 t.SetCurrentClusterIndex1(0);
1808 Double_t xt=t.GetX();
1809 Int_t row = GetRowNumber(xt)-1;
1810 Double_t ymax= GetMaxY(nr);
1812 if (row < nr) return 1; // don't prolongate if not information until now -
1813 // if (TMath::Abs(t.GetSnp())>0.9 && t.GetNumberOfClusters()>40. && fIteration!=2) {
1815 // return 0; // not prolongate strongly inclined tracks
1817 // if (TMath::Abs(t.GetSnp())>0.95) {
1819 // return 0; // not prolongate strongly inclined tracks
1820 // }// patch 28 fev 06
1822 Double_t x= GetXrow(nr);
1824 //t.PropagateTo(x+0.02);
1825 //t.PropagateTo(x+0.01);
1826 if (!t.PropagateTo(x)){
1833 if (TMath::Abs(y)>ymax){
1835 t.SetRelativeSector((t.GetRelativeSector()+1) % fN);
1836 if (!t.Rotate(fSectors->GetAlpha()))
1838 } else if (y <-ymax) {
1839 t.SetRelativeSector((t.GetRelativeSector()-1+fN) % fN);
1840 if (!t.Rotate(-fSectors->GetAlpha()))
1843 // if (!t.PropagateTo(x)){
1850 if (TMath::Abs(t.GetSnp())>AliTPCReconstructor::GetMaxSnpTracker()) return 0;
1852 if (!IsActive(t.GetRelativeSector(),nr)) {
1854 t.SetClusterIndex2(nr,-1);
1857 //AliInfo(Form("A - Sector%d phi %f - alpha %f", t.fRelativeSector,y/x, t.GetAlpha()));
1859 AliTPCtrackerRow &krow=GetRow(t.GetRelativeSector(),nr);
1861 if (TMath::Abs(TMath::Abs(y)-ymax)<krow.GetDeadZone()){
1863 t.SetClusterIndex2(nr,-1);
1868 if (TMath::Abs(t.GetZ())<(AliTPCReconstructor::GetCtgRange()*t.GetX()+10) && (TMath::Abs(t.GetSnp())<AliTPCReconstructor::GetMaxSnpTracker())) t.SetNFoundable(t.GetNFoundable()+1);
1874 if ( (nr%2==0) || t.GetNumberOfClusters()<2){
1875 // t.fCurrentSigmaY = GetSigmaY(&t);
1876 //t.fCurrentSigmaZ = GetSigmaZ(&t);
1880 AliTPCclusterMI *cl=0;
1883 Double_t roady = 1.;
1884 Double_t roadz = 1.;
1888 index = t.GetClusterIndex2(nr);
1889 if ( (index>0) && (index&0x8000)==0){
1890 cl = t.GetClusterPointer(nr);
1891 if ( (cl==0) && (index>0)) cl = GetClusterMI(index);
1892 t.SetCurrentClusterIndex1(index);
1894 t.SetCurrentCluster(cl);
1900 // if (index<0) return 0;
1901 UInt_t uindex = TMath::Abs(index);
1904 //cl = krow.FindNearest2(y+10,z,roady,roadz,uindex);
1905 cl = krow.FindNearest2(y,z,roady,roadz,uindex);
1908 if (cl) t.SetCurrentClusterIndex1(krow.GetIndex(uindex));
1909 t.SetCurrentCluster(cl);
1915 Int_t AliTPCtrackerMI::FollowToNextCluster(AliTPCseed & t, Int_t nr) {
1916 //-----------------------------------------------------------------
1917 // This function tries to find a track prolongation to next pad row
1918 //-----------------------------------------------------------------
1920 //update error according neighborhoud
1922 if (t.GetCurrentCluster()) {
1924 Int_t accept = AcceptCluster(&t,t.GetCurrentCluster());
1926 if (t.GetCurrentCluster()->IsUsed(10)){
1931 t.SetNShared(t.GetNShared()+1);
1932 if (t.GetNShared()>0.7*t.GetNumberOfClusters()) {
1937 if (fIteration>0) accept = 0;
1938 if (accept<3) UpdateTrack(&t,accept);
1942 if ( ( (t.GetSigmaY2()+t.GetSigmaZ2())>0.16)&& t.GetNumberOfClusters()>18) t.SetRemoval(10);
1943 if ( t.GetChi2()/t.GetNumberOfClusters()>6 &&t.GetNumberOfClusters()>18) t.SetRemoval(10);
1945 if (( (t.GetNFoundable()*0.5 > t.GetNumberOfClusters()) || t.GetNoCluster()>15)) t.SetRemoval(10);
1953 //_____________________________________________________________________________
1954 Int_t AliTPCtrackerMI::FollowProlongation(AliTPCseed& t, Int_t rf, Int_t step) {
1955 //-----------------------------------------------------------------
1956 // This function tries to find a track prolongation.
1957 //-----------------------------------------------------------------
1958 Double_t xt=t.GetX();
1960 Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
1961 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
1962 if (alpha < 0. ) alpha += 2.*TMath::Pi();
1964 t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN);
1966 Int_t first = GetRowNumber(xt)-1;
1967 for (Int_t nr= first; nr>=rf; nr-=step) {
1969 if (t.GetKinkIndexes()[0]>0){
1970 for (Int_t i=0;i<3;i++){
1971 Int_t index = t.GetKinkIndexes()[i];
1972 if (index==0) break;
1973 if (index<0) continue;
1975 AliKink * kink = (AliKink*)fEvent->GetKink(index-1);
1977 printf("PROBLEM\n");
1980 Int_t kinkrow = kink->GetTPCRow0()+2+Int_t(0.5/(0.05+kink->GetAngle(2)));
1982 AliExternalTrackParam paramd(t);
1983 kink->SetDaughter(paramd);
1984 kink->SetStatus(2,5);
1991 if (nr==80) t.UpdateReference();
1992 if (nr<fInnerSec->GetNRows())
1993 fSectors = fInnerSec;
1995 fSectors = fOuterSec;
1996 if (FollowToNext(t,nr)==0)
2005 //_____________________________________________________________________________
2006 Int_t AliTPCtrackerMI::FollowProlongationFast(AliTPCseed& t, Int_t rf, Int_t step) {
2007 //-----------------------------------------------------------------
2008 // This function tries to find a track prolongation.
2009 //-----------------------------------------------------------------
2010 Double_t xt=t.GetX();
2012 Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
2013 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
2014 if (alpha < 0. ) alpha += 2.*TMath::Pi();
2015 t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN);
2017 for (Int_t nr=GetRowNumber(xt)-1; nr>=rf; nr-=step) {
2019 if (FollowToNextFast(t,nr)==0)
2020 if (!t.IsActive()) return 0;
2030 Int_t AliTPCtrackerMI::FollowBackProlongation(AliTPCseed& t, Int_t rf) {
2031 //-----------------------------------------------------------------
2032 // This function tries to find a track prolongation.
2033 //-----------------------------------------------------------------
2035 Double_t xt=t.GetX();
2036 Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
2037 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
2038 if (alpha < 0. ) alpha += 2.*TMath::Pi();
2039 t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN);
2041 Int_t first = t.GetFirstPoint();
2042 if (first<GetRowNumber(xt)+1) first = GetRowNumber(xt)+1;
2044 if (first<0) first=0;
2045 for (Int_t nr=first; nr<=rf; nr++) {
2046 // if ( (TMath::Abs(t.GetSnp())>0.95)) break;//patch 28 fev 06
2047 if (t.GetKinkIndexes()[0]<0){
2048 for (Int_t i=0;i<3;i++){
2049 Int_t index = t.GetKinkIndexes()[i];
2050 if (index==0) break;
2051 if (index>0) continue;
2052 index = TMath::Abs(index);
2053 AliKink * kink = (AliKink*)fEvent->GetKink(index-1);
2055 printf("PROBLEM\n");
2058 Int_t kinkrow = kink->GetTPCRow0()-2-Int_t(0.5/(0.05+kink->GetAngle(2)));
2060 AliExternalTrackParam paramm(t);
2061 kink->SetMother(paramm);
2062 kink->SetStatus(2,1);
2069 if (nr<fInnerSec->GetNRows())
2070 fSectors = fInnerSec;
2072 fSectors = fOuterSec;
2083 Float_t AliTPCtrackerMI::OverlapFactor(AliTPCseed * s1, AliTPCseed * s2, Int_t &sum1, Int_t & sum2)
2091 Float_t dz2 =(s1->GetZ() - s2->GetZ());
2094 Float_t dy2 =TMath::Abs((s1->GetY() - s2->GetY()));
2096 Float_t distance = TMath::Sqrt(dz2+dy2);
2097 if (distance>4.) return 0; // if there are far away - not overlap - to reduce combinatorics
2100 Int_t firstpoint = TMath::Min(s1->GetFirstPoint(),s2->GetFirstPoint());
2101 Int_t lastpoint = TMath::Max(s1->GetLastPoint(),s2->GetLastPoint());
2106 if (firstpoint>lastpoint) {
2107 firstpoint =lastpoint;
2112 for (Int_t i=firstpoint-1;i<lastpoint+1;i++){
2113 if (s1->GetClusterIndex2(i)>0) sum1++;
2114 if (s2->GetClusterIndex2(i)>0) sum2++;
2115 if (s1->GetClusterIndex2(i)==s2->GetClusterIndex2(i) && s1->GetClusterIndex2(i)>0) {
2119 if (sum<5) return 0;
2121 Float_t summin = TMath::Min(sum1+1,sum2+1);
2122 Float_t ratio = (sum+1)/Float_t(summin);
2126 void AliTPCtrackerMI::SignShared(AliTPCseed * s1, AliTPCseed * s2)
2130 Float_t thetaCut = 0.2;//+10.*TMath::Sqrt(s1->GetSigmaTglZ()+ s2->GetSigmaTglZ());
2131 if (TMath::Abs(s1->GetTgl()-s2->GetTgl())>thetaCut) return;
2132 Float_t minCl = TMath::Min(s1->GetNumberOfClusters(),s2->GetNumberOfClusters());
2133 Int_t cutN0 = TMath::Max(5,TMath::Nint(0.1*minCl));
2138 //Int_t firstpoint = TMath::Max(s1->GetFirstPoint(),s2->GetFirstPoint());
2139 //Int_t lastpoint = TMath::Min(s1->GetLastPoint(),s2->GetLastPoint());
2140 Int_t firstpoint = 0;
2141 Int_t lastpoint = 160;
2143 // if (firstpoint>=lastpoint-5) return;;
2145 for (Int_t i=firstpoint;i<lastpoint;i++){
2146 // if ( (s1->GetClusterIndex2(i)&0xFFFF8FFF)==(s2->GetClusterIndex2(i)&0xFFFF8FFF) && s1->GetClusterIndex2(i)>0) {
2147 if ( (s1->GetClusterIndex2(i))==(s2->GetClusterIndex2(i)) && s1->GetClusterIndex2(i)>0) {
2149 s1->SetSharedMapBit(i, kTRUE);
2150 s2->SetSharedMapBit(i, kTRUE);
2152 if (s1->GetClusterIndex2(i)>0)
2153 s1->SetClusterMapBit(i, kTRUE);
2155 if (sumshared>cutN0){
2158 for (Int_t i=firstpoint;i<lastpoint;i++){
2159 // if ( (s1->GetClusterIndex2(i)&0xFFFF8FFF)==(s2->GetClusterIndex2(i)&0xFFFF8FFF) && s1->GetClusterIndex2(i)>0) {
2160 if ( (s1->GetClusterIndex2(i))==(s2->GetClusterIndex2(i)) && s1->GetClusterIndex2(i)>0) {
2161 AliTPCTrackerPoint *p1 = s1->GetTrackPoint(i);
2162 AliTPCTrackerPoint *p2 = s2->GetTrackPoint(i);;
2163 if (s1->IsActive()&&s2->IsActive()){
2164 p1->SetShared(kTRUE);
2165 p2->SetShared(kTRUE);
2171 if (sumshared>cutN0){
2172 for (Int_t i=0;i<4;i++){
2173 if (s1->GetOverlapLabel(3*i)==0){
2174 s1->SetOverlapLabel(3*i, s2->GetLabel());
2175 s1->SetOverlapLabel(3*i+1,sumshared);
2176 s1->SetOverlapLabel(3*i+2,s2->GetUniqueID());
2180 for (Int_t i=0;i<4;i++){
2181 if (s2->GetOverlapLabel(3*i)==0){
2182 s2->SetOverlapLabel(3*i, s1->GetLabel());
2183 s2->SetOverlapLabel(3*i+1,sumshared);
2184 s2->SetOverlapLabel(3*i+2,s1->GetUniqueID());
2191 void AliTPCtrackerMI::SignShared(TObjArray * arr)
2194 //sort trackss according sectors
2196 for (Int_t i=0; i<arr->GetEntriesFast(); i++) {
2197 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2199 //if (pt) RotateToLocal(pt);
2203 arr->Sort(); // sorting according relative sectors
2204 arr->Expand(arr->GetEntries());
2207 Int_t nseed=arr->GetEntriesFast();
2208 for (Int_t i=0; i<nseed; i++) {
2209 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2211 for (Int_t j=0;j<=12;j++){
2212 pt->SetOverlapLabel(j,0);
2215 for (Int_t i=0; i<nseed; i++) {
2216 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2218 if (pt->GetRemoval()>10) continue;
2219 for (Int_t j=i+1; j<nseed; j++){
2220 AliTPCseed *pt2=(AliTPCseed*)arr->UncheckedAt(j);
2221 if (TMath::Abs(pt->GetRelativeSector()-pt2->GetRelativeSector())>1) continue;
2223 if (pt2->GetRemoval()<=10) {
2224 //if ( TMath::Abs(pt->GetRelativeSector()-pt2->GetRelativeSector())>0) break;
2232 void AliTPCtrackerMI::SortTracks(TObjArray * arr, Int_t mode) const
2235 //sort tracks in array according mode criteria
2236 Int_t nseed = arr->GetEntriesFast();
2237 for (Int_t i=0; i<nseed; i++) {
2238 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2249 void AliTPCtrackerMI::RemoveUsed2(TObjArray * arr, Float_t factor1, Float_t factor2, Int_t minimal)
2252 // Loop over all tracks and remove overlaped tracks (with lower quality)
2254 // 1. Unsign clusters
2255 // 2. Sort tracks according quality
2256 // Quality is defined by the number of cluster between first and last points
2258 // 3. Loop over tracks - decreasing quality order
2259 // a.) remove - If the fraction of shared cluster less than factor (1- n or 2)
2260 // b.) remove - If the minimal number of clusters less than minimal and not ITS
2261 // c.) if track accepted - sign clusters
2263 //Called in - AliTPCtrackerMI::Clusters2Tracks()
2264 // - AliTPCtrackerMI::PropagateBack()
2265 // - AliTPCtrackerMI::RefitInward()
2271 Int_t nseed = arr->GetEntriesFast();
2272 Float_t * quality = new Float_t[nseed];
2273 Int_t * indexes = new Int_t[nseed];
2277 for (Int_t i=0; i<nseed; i++) {
2278 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2283 pt->UpdatePoints(); //select first last max dens points
2284 Float_t * points = pt->GetPoints();
2285 if (points[3]<0.8) quality[i] =-1;
2286 quality[i] = (points[2]-points[0])+pt->GetNumberOfClusters();
2287 //prefer high momenta tracks if overlaps
2288 quality[i] *= TMath::Sqrt(TMath::Abs(pt->Pt())+0.5);
2290 TMath::Sort(nseed,quality,indexes);
2293 for (Int_t itrack=0; itrack<nseed; itrack++) {
2294 Int_t trackindex = indexes[itrack];
2295 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(trackindex);
2298 if (quality[trackindex]<0){
2300 delete arr->RemoveAt(trackindex);
2303 arr->RemoveAt(trackindex);
2309 Int_t first = Int_t(pt->GetPoints()[0]);
2310 Int_t last = Int_t(pt->GetPoints()[2]);
2311 Double_t factor = (pt->GetBConstrain()) ? factor1: factor2;
2313 Int_t found,foundable,shared;
2314 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
2315 // pt->GetClusterStatistic(0,160, found, foundable,shared,kFALSE);
2316 Bool_t itsgold =kFALSE;
2319 if (pt->GetESD()->GetITSclusters(dummy)>4) itsgold= kTRUE;
2323 if (Float_t(shared+1)/Float_t(found+1)>factor){
2324 if (pt->GetKinkIndexes()[0]!=0) continue; //don't remove tracks - part of the kinks
2325 delete arr->RemoveAt(trackindex);
2328 if (pt->GetNumberOfClusters()<50&&(found-0.5*shared)<minimal){ //remove short tracks
2329 if (pt->GetKinkIndexes()[0]!=0) continue; //don't remove tracks - part of the kinks
2330 delete arr->RemoveAt(trackindex);
2336 //if (sharedfactor>0.4) continue;
2337 if (pt->GetKinkIndexes()[0]>0) continue;
2338 //Remove tracks with undefined properties - seems
2339 if (pt->GetSigmaY2()<kAlmost0) continue; // ? what is the origin ?
2341 for (Int_t i=first; i<last; i++) {
2342 Int_t index=pt->GetClusterIndex2(i);
2343 // if (index<0 || index&0x8000 ) continue;
2344 if (index<0 || index&0x8000 ) continue;
2345 AliTPCclusterMI *c= pt->GetClusterPointer(i);
2352 Info("RemoveUsed","\n*****\nNumber of good tracks after shared removal\t%d\n",fNtracks);
2358 void AliTPCtrackerMI::UnsignClusters()
2361 // loop over all clusters and unsign them
2364 for (Int_t sec=0;sec<fkNIS;sec++){
2365 for (Int_t row=0;row<fInnerSec->GetNRows();row++){
2366 AliTPCclusterMI *cl = fInnerSec[sec][row].GetClusters1();
2367 for (Int_t icl =0;icl< fInnerSec[sec][row].GetN1();icl++)
2368 // if (cl[icl].IsUsed(10))
2370 cl = fInnerSec[sec][row].GetClusters2();
2371 for (Int_t icl =0;icl< fInnerSec[sec][row].GetN2();icl++)
2372 //if (cl[icl].IsUsed(10))
2377 for (Int_t sec=0;sec<fkNOS;sec++){
2378 for (Int_t row=0;row<fOuterSec->GetNRows();row++){
2379 AliTPCclusterMI *cl = fOuterSec[sec][row].GetClusters1();
2380 for (Int_t icl =0;icl< fOuterSec[sec][row].GetN1();icl++)
2381 //if (cl[icl].IsUsed(10))
2383 cl = fOuterSec[sec][row].GetClusters2();
2384 for (Int_t icl =0;icl< fOuterSec[sec][row].GetN2();icl++)
2385 //if (cl[icl].IsUsed(10))
2394 void AliTPCtrackerMI::SignClusters(TObjArray * arr, Float_t fnumber, Float_t fdensity)
2397 //sign clusters to be "used"
2399 // snumber and sdensity sign number of sigmas - bellow mean value to be accepted
2400 // loop over "primaries"
2414 Int_t nseed = arr->GetEntriesFast();
2415 for (Int_t i=0; i<nseed; i++) {
2416 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2420 if (!(pt->IsActive())) continue;
2421 Float_t dens = pt->GetNumberOfClusters()/Float_t(pt->GetNFoundable());
2422 if ( (dens>0.7) && (pt->GetNumberOfClusters()>70)){
2424 sumdens2+= dens*dens;
2425 sumn += pt->GetNumberOfClusters();
2426 sumn2 += pt->GetNumberOfClusters()*pt->GetNumberOfClusters();
2427 Float_t chi2 = pt->GetChi2()/pt->GetNumberOfClusters();
2430 sumchi2 +=chi2*chi2;
2435 Float_t mdensity = 0.9;
2436 Float_t meann = 130;
2437 Float_t meanchi = 1;
2438 Float_t sdensity = 0.1;
2439 Float_t smeann = 10;
2440 Float_t smeanchi =0.4;
2444 mdensity = sumdens/sum;
2446 meanchi = sumchi/sum;
2448 sdensity = sumdens2/sum-mdensity*mdensity;
2450 sdensity = TMath::Sqrt(sdensity);
2454 smeann = sumn2/sum-meann*meann;
2456 smeann = TMath::Sqrt(smeann);
2460 smeanchi = sumchi2/sum - meanchi*meanchi;
2462 smeanchi = TMath::Sqrt(smeanchi);
2468 //REMOVE SHORT DELTAS or tracks going out of sensitive volume of TPC
2470 for (Int_t i=0; i<nseed; i++) {
2471 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2475 if (pt->GetBSigned()) continue;
2476 if (pt->GetBConstrain()) continue;
2477 //if (!(pt->IsActive())) continue;
2479 Int_t found,foundable,shared;
2480 pt->GetClusterStatistic(0,160,found, foundable,shared);
2481 if (shared/float(found)>0.3) {
2482 if (shared/float(found)>0.9 ){
2483 //delete arr->RemoveAt(i);
2488 Bool_t isok =kFALSE;
2489 if ( (pt->GetNShared()/pt->GetNumberOfClusters()<0.5) &&pt->GetNumberOfClusters()>60)
2491 if ((TMath::Abs(1/pt->GetC())<100.) && (pt->GetNShared()/pt->GetNumberOfClusters()<0.7))
2493 if (TMath::Abs(pt->GetZ()/pt->GetX())>1.1)
2495 if ( (TMath::Abs(pt->GetSnp()>0.7) && pt->GetD(0,0)>60.))
2499 for (Int_t i=0; i<160; i++) {
2500 Int_t index=pt->GetClusterIndex2(i);
2501 if (index<0) continue;
2502 AliTPCclusterMI *c= pt->GetClusterPointer(i);
2504 //if (!(c->IsUsed(10))) c->Use();
2511 Double_t maxchi = meanchi+2.*smeanchi;
2513 for (Int_t i=0; i<nseed; i++) {
2514 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2518 //if (!(pt->IsActive())) continue;
2519 if (pt->GetBSigned()) continue;
2520 Double_t chi = pt->GetChi2()/pt->GetNumberOfClusters();
2521 if (chi>maxchi) continue;
2524 Float_t dens = pt->GetNumberOfClusters()/Float_t(pt->GetNFoundable());
2526 //sign only tracks with enoug big density at the beginning
2528 if ((pt->GetDensityFirst(40)<0.75) && pt->GetNumberOfClusters()<meann) continue;
2531 Double_t mindens = TMath::Max(double(mdensity-sdensity*fdensity*bfactor),0.65);
2532 Double_t minn = TMath::Max(Int_t(meann-fnumber*smeann*bfactor),50);
2534 // if (pt->fBConstrain) mindens = TMath::Max(mdensity-sdensity*fdensity*bfactor,0.65);
2535 if ( (pt->GetRemoval()==10) && (pt->GetSnp()>0.8)&&(dens>mindens))
2538 if ((dens>mindens && pt->GetNumberOfClusters()>minn) && chi<maxchi ){
2539 //Int_t noc=pt->GetNumberOfClusters();
2540 pt->SetBSigned(kTRUE);
2541 for (Int_t i=0; i<160; i++) {
2543 Int_t index=pt->GetClusterIndex2(i);
2544 if (index<0) continue;
2545 AliTPCclusterMI *c= pt->GetClusterPointer(i);
2547 // if (!(c->IsUsed(10))) c->Use();
2552 // gLastCheck = nseed;
2560 void AliTPCtrackerMI::StopNotActive(TObjArray * arr, Int_t row0, Float_t th0, Float_t th1, Float_t th2) const
2562 // stop not active tracks
2563 // take th1 as threshold for number of founded to number of foundable on last 10 active rows
2564 // take th2 as threshold for number of founded to number of foundable on last 20 active rows
2565 Int_t nseed = arr->GetEntriesFast();
2567 for (Int_t i=0; i<nseed; i++) {
2568 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2572 if (!(pt->IsActive())) continue;
2573 StopNotActive(pt,row0,th0, th1,th2);
2579 void AliTPCtrackerMI::StopNotActive(AliTPCseed * seed, Int_t row0, Float_t th0, Float_t th1,
2582 // stop not active tracks
2583 // take th1 as threshold for number of founded to number of foundable on last 10 active rows
2584 // take th2 as threshold for number of founded to number of foundable on last 20 active rows
2587 Int_t foundable = 0;
2588 Int_t maxindex = seed->GetLastPoint(); //last foundable row
2589 if (seed->GetNFoundable()*th0 > seed->GetNumberOfClusters()) {
2590 seed->Desactivate(10) ;
2594 for (Int_t i=row0; i<maxindex; i++){
2595 Int_t index = seed->GetClusterIndex2(i);
2596 if (index!=-1) foundable++;
2598 if (foundable<=30) sumgood1++;
2599 if (foundable<=50) {
2606 if (foundable>=30.){
2607 if (sumgood1<(th1*30.)) seed->Desactivate(10);
2610 if (sumgood2<(th2*50.)) seed->Desactivate(10);
2614 Int_t AliTPCtrackerMI::RefitInward(AliESDEvent *event)
2617 // back propagation of ESD tracks
2620 const Int_t kMaxFriendTracks=2000;
2624 //PrepareForProlongation(fSeeds,1);
2625 PropagateForward2(fSeeds);
2626 RemoveUsed2(fSeeds,0.4,0.4,20);
2628 TObjArray arraySeed(fSeeds->GetEntries());
2629 for (Int_t i=0;i<fSeeds->GetEntries();i++) {
2630 arraySeed.AddAt(fSeeds->At(i),i);
2632 SignShared(&arraySeed);
2633 // FindCurling(fSeeds, event,2); // find multi found tracks
2634 FindSplitted(fSeeds, event,2); // find multi found tracks
2637 Int_t nseed = fSeeds->GetEntriesFast();
2638 for (Int_t i=0;i<nseed;i++){
2639 AliTPCseed * seed = (AliTPCseed*) fSeeds->UncheckedAt(i);
2640 if (!seed) continue;
2641 if (seed->GetKinkIndex(0)>0) UpdateKinkQualityD(seed); // update quality informations for kinks
2643 seed->PropagateTo(fParam->GetInnerRadiusLow());
2644 seed->UpdatePoints();
2645 AddCovariance(seed);
2647 AliESDtrack *esd=event->GetTrack(i);
2648 seed->CookdEdx(0.02,0.6);
2649 CookLabel(seed,0.1); //For comparison only
2650 esd->SetTPCClusterMap(seed->GetClusterMap());
2651 esd->SetTPCSharedMap(seed->GetSharedMap());
2653 if (AliTPCReconstructor::StreamLevel()>1 && seed!=0&&esd!=0) {
2654 TTreeSRedirector &cstream = *fDebugStreamer;
2661 if (seed->GetNumberOfClusters()>15){
2662 esd->UpdateTrackParams(seed,AliESDtrack::kTPCrefit);
2663 esd->SetTPCPoints(seed->GetPoints());
2664 esd->SetTPCPointsF(seed->GetNFoundable());
2665 Int_t ndedx = seed->GetNCDEDX(0)+seed->GetNCDEDX(1)+seed->GetNCDEDX(2)+seed->GetNCDEDX(3);
2666 Float_t sdedx = (seed->GetSDEDX(0)+seed->GetSDEDX(1)+seed->GetSDEDX(2)+seed->GetSDEDX(3))*0.25;
2667 Float_t dedx = seed->GetdEdx();
2668 esd->SetTPCsignal(dedx, sdedx, ndedx);
2670 // add seed to the esd track in Calib level
2672 Bool_t storeFriend = gRandom->Rndm()<(kMaxFriendTracks)/Float_t(nseed);
2673 if (AliTPCReconstructor::StreamLevel()>0 &&storeFriend){
2674 AliTPCseed * seedCopy = new AliTPCseed(*seed, kTRUE);
2675 esd->AddCalibObject(seedCopy);
2680 //printf("problem\n");
2683 //FindKinks(fSeeds,event);
2684 Info("RefitInward","Number of refitted tracks %d",ntracks);
2689 Int_t AliTPCtrackerMI::PropagateBack(AliESDEvent *event)
2692 // back propagation of ESD tracks
2698 PropagateBack(fSeeds);
2699 RemoveUsed2(fSeeds,0.4,0.4,20);
2700 //FindCurling(fSeeds, fEvent,1);
2701 FindSplitted(fSeeds, event,1); // find multi found tracks
2704 Int_t nseed = fSeeds->GetEntriesFast();
2706 for (Int_t i=0;i<nseed;i++){
2707 AliTPCseed * seed = (AliTPCseed*) fSeeds->UncheckedAt(i);
2708 if (!seed) continue;
2709 if (seed->GetKinkIndex(0)<0) UpdateKinkQualityM(seed); // update quality informations for kinks
2710 seed->UpdatePoints();
2711 AddCovariance(seed);
2712 AliESDtrack *esd=event->GetTrack(i);
2713 seed->CookdEdx(0.02,0.6);
2714 CookLabel(seed,0.1); //For comparison only
2715 if (seed->GetNumberOfClusters()>15){
2716 esd->UpdateTrackParams(seed,AliESDtrack::kTPCout);
2717 esd->SetTPCPoints(seed->GetPoints());
2718 esd->SetTPCPointsF(seed->GetNFoundable());
2719 Int_t ndedx = seed->GetNCDEDX(0)+seed->GetNCDEDX(1)+seed->GetNCDEDX(2)+seed->GetNCDEDX(3);
2720 Float_t sdedx = (seed->GetSDEDX(0)+seed->GetSDEDX(1)+seed->GetSDEDX(2)+seed->GetSDEDX(3))*0.25;
2721 Float_t dedx = seed->GetdEdx();
2722 esd->SetTPCsignal(dedx, sdedx, ndedx);
2724 Int_t eventnumber = event->GetEventNumberInFile();// patch 28 fev 06
2725 // This is most likely NOT the event number you'd like to use. It has nothing to do with the 'real' event number
2726 if (AliTPCReconstructor::StreamLevel()>1 && esd) {
2727 (*fDebugStreamer)<<"Cback"<<
2730 "EventNrInFile="<<eventnumber<<
2735 //FindKinks(fSeeds,event);
2736 Info("PropagateBack","Number of back propagated tracks %d",ntracks);
2743 void AliTPCtrackerMI::DeleteSeeds()
2748 Int_t nseed = fSeeds->GetEntriesFast();
2749 for (Int_t i=0;i<nseed;i++){
2750 AliTPCseed * seed = (AliTPCseed*)fSeeds->At(i);
2751 if (seed) delete fSeeds->RemoveAt(i);
2758 void AliTPCtrackerMI::ReadSeeds(AliESDEvent *event, Int_t direction)
2761 //read seeds from the event
2763 Int_t nentr=event->GetNumberOfTracks();
2765 Info("ReadSeeds", "Number of ESD tracks: %d\n", nentr);
2770 fSeeds = new TObjArray(nentr);
2774 for (Int_t i=0; i<nentr; i++) {
2775 AliESDtrack *esd=event->GetTrack(i);
2776 ULong_t status=esd->GetStatus();
2777 if (!(status&AliESDtrack::kTPCin)) continue;
2778 AliTPCtrack t(*esd);
2779 t.SetNumberOfClusters(0);
2780 // AliTPCseed *seed = new AliTPCseed(t,t.GetAlpha());
2781 AliTPCseed *seed = new AliTPCseed(t/*,t.GetAlpha()*/);
2782 seed->SetUniqueID(esd->GetID());
2783 AddCovariance(seed); //add systematic ucertainty
2784 for (Int_t ikink=0;ikink<3;ikink++) {
2785 Int_t index = esd->GetKinkIndex(ikink);
2786 seed->GetKinkIndexes()[ikink] = index;
2787 if (index==0) continue;
2788 index = TMath::Abs(index);
2789 AliESDkink * kink = fEvent->GetKink(index-1);
2790 if (kink&&esd->GetKinkIndex(ikink)<0){
2791 if ((status & AliESDtrack::kTRDrefit) != 0) kink->SetStatus(1,2);
2792 if ((status & AliESDtrack::kITSout) != 0) kink->SetStatus(1,0);
2794 if (kink&&esd->GetKinkIndex(ikink)>0){
2795 if ((status & AliESDtrack::kTRDrefit) != 0) kink->SetStatus(1,6);
2796 if ((status & AliESDtrack::kITSout) != 0) kink->SetStatus(1,4);
2800 if (((status&AliESDtrack::kITSout)==0)&&(direction==1)) seed->ResetCovariance(10.);
2801 if ( direction ==2 &&(status & AliESDtrack::kTRDrefit) == 0 ) seed->ResetCovariance(10.);
2802 if ( direction ==2 && ((status & AliESDtrack::kTPCout) == 0) ) {
2807 if ( direction ==2 &&(status & AliESDtrack::kTRDrefit) > 0 ) {
2808 Double_t par0[5],par1[5],alpha,x;
2809 esd->GetInnerExternalParameters(alpha,x,par0);
2810 esd->GetExternalParameters(x,par1);
2811 Double_t delta1 = TMath::Abs(par0[4]-par1[4])/(0.000000001+TMath::Abs(par0[4]+par1[4]));
2812 Double_t delta2 = TMath::Abs(par0[3]-par1[3]);
2814 if (esd->GetTRDncls()>0) trdchi2 = esd->GetTRDchi2()/esd->GetTRDncls();
2815 //reset covariance if suspicious
2816 if ( (delta1>0.1) || (delta2>0.006) ||trdchi2>7.)
2817 seed->ResetCovariance(10.);
2822 // rotate to the local coordinate system
2824 fSectors=fInnerSec; fN=fkNIS;
2825 Double_t alpha=seed->GetAlpha() - fSectors->GetAlphaShift();
2826 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
2827 if (alpha < 0. ) alpha += 2.*TMath::Pi();
2828 Int_t ns=Int_t(alpha/fSectors->GetAlpha())%fN;
2829 alpha =ns*fSectors->GetAlpha() + fSectors->GetAlphaShift();
2830 if (alpha<-TMath::Pi()) alpha += 2*TMath::Pi();
2831 if (alpha>=TMath::Pi()) alpha -= 2*TMath::Pi();
2832 alpha-=seed->GetAlpha();
2833 if (!seed->Rotate(alpha)) {
2839 if (esd->GetKinkIndex(0)<=0){
2840 for (Int_t irow=0;irow<160;irow++){
2841 Int_t index = seed->GetClusterIndex2(irow);
2844 AliTPCclusterMI * cl = GetClusterMI(index);
2845 seed->SetClusterPointer(irow,cl);
2847 if ((index & 0x8000)==0){
2848 cl->Use(10); // accepted cluster
2850 cl->Use(6); // close cluster not accepted
2853 Info("ReadSeeds","Not found cluster");
2858 fSeeds->AddAt(seed,i);
2864 //_____________________________________________________________________________
2865 void AliTPCtrackerMI::MakeSeeds3(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, Float_t cuts[4],
2866 Float_t deltay, Int_t ddsec) {
2867 //-----------------------------------------------------------------
2868 // This function creates track seeds.
2869 // SEEDING WITH VERTEX CONSTRAIN
2870 //-----------------------------------------------------------------
2871 // cuts[0] - fP4 cut
2872 // cuts[1] - tan(phi) cut
2873 // cuts[2] - zvertex cut
2874 // cuts[3] - fP3 cut
2882 Double_t x[5], c[15];
2883 // Int_t di = i1-i2;
2885 AliTPCseed * seed = new AliTPCseed();
2886 Double_t alpha=fSectors->GetAlpha(), shift=fSectors->GetAlphaShift();
2887 Double_t cs=cos(alpha), sn=sin(alpha);
2889 // Double_t x1 =fOuterSec->GetX(i1);
2890 //Double_t xx2=fOuterSec->GetX(i2);
2892 Double_t x1 =GetXrow(i1);
2893 Double_t xx2=GetXrow(i2);
2895 Double_t x3=GetX(), y3=GetY(), z3=GetZ();
2897 Int_t imiddle = (i2+i1)/2; //middle pad row index
2898 Double_t xm = GetXrow(imiddle); // radius of middle pad-row
2899 const AliTPCtrackerRow& krm=GetRow(sec,imiddle); //middle pad -row
2903 const AliTPCtrackerRow& kr1=GetRow(ns,i1);
2904 Double_t ymax = GetMaxY(i1)-kr1.GetDeadZone()-1.5;
2905 Double_t ymaxm = GetMaxY(imiddle)-kr1.GetDeadZone()-1.5;
2908 // change cut on curvature if it can't reach this layer
2909 // maximal curvature set to reach it
2910 Double_t dvertexmax = TMath::Sqrt((x1-x3)*(x1-x3)+(ymax+5-y3)*(ymax+5-y3));
2911 if (dvertexmax*0.5*cuts[0]>0.85){
2912 cuts[0] = 0.85/(dvertexmax*0.5+1.);
2914 Double_t r2min = 1/(cuts[0]*cuts[0]); //minimal square of radius given by cut
2917 if (deltay>0) ddsec = 0;
2918 // loop over clusters
2919 for (Int_t is=0; is < kr1; is++) {
2921 if (kr1[is]->IsUsed(10)) continue;
2922 Double_t y1=kr1[is]->GetY(), z1=kr1[is]->GetZ();
2923 //if (TMath::Abs(y1)>ymax) continue;
2925 if (deltay>0 && TMath::Abs(ymax-TMath::Abs(y1))> deltay ) continue; // seed only at the edge
2927 // find possible directions
2928 Float_t anglez = (z1-z3)/(x1-x3);
2929 Float_t extraz = z1 - anglez*(x1-xx2); // extrapolated z
2932 //find rotation angles relative to line given by vertex and point 1
2933 Double_t dvertex2 = (x1-x3)*(x1-x3)+(y1-y3)*(y1-y3);
2934 Double_t dvertex = TMath::Sqrt(dvertex2);
2935 Double_t angle13 = TMath::ATan((y1-y3)/(x1-x3));
2936 Double_t cs13 = cos(-angle13), sn13 = sin(-angle13);
2939 // loop over 2 sectors
2945 Double_t dddz1=0; // direction of delta inclination in z axis
2952 for (Int_t dsec = dsec1; dsec<=dsec2;dsec++){
2953 Int_t sec2 = sec + dsec;
2955 // AliTPCtrackerRow& kr2 = fOuterSec[(sec2+fkNOS)%fkNOS][i2];
2956 //AliTPCtrackerRow& kr2m = fOuterSec[(sec2+fkNOS)%fkNOS][imiddle];
2957 AliTPCtrackerRow& kr2 = GetRow((sec2+fkNOS)%fkNOS,i2);
2958 AliTPCtrackerRow& kr2m = GetRow((sec2+fkNOS)%fkNOS,imiddle);
2959 Int_t index1 = TMath::Max(kr2.Find(extraz-0.6-dddz1*TMath::Abs(z1)*0.05)-1,0);
2960 Int_t index2 = TMath::Min(kr2.Find(extraz+0.6+dddz2*TMath::Abs(z1)*0.05)+1,kr2);
2962 // rotation angles to p1-p3
2963 Double_t cs13r = cos(-angle13+dsec*alpha)/dvertex, sn13r = sin(-angle13+dsec*alpha)/dvertex;
2964 Double_t x2, y2, z2;
2966 // Double_t dymax = maxangle*TMath::Abs(x1-xx2);
2969 Double_t dxx0 = (xx2-x3)*cs13r;
2970 Double_t dyy0 = (xx2-x3)*sn13r;
2971 for (Int_t js=index1; js < index2; js++) {
2972 const AliTPCclusterMI *kcl = kr2[js];
2973 if (kcl->IsUsed(10)) continue;
2975 //calcutate parameters
2977 Double_t yy0 = dyy0 +(kcl->GetY()-y3)*cs13r;
2979 if (TMath::Abs(yy0)<0.000001) continue;
2980 Double_t xx0 = dxx0 -(kcl->GetY()-y3)*sn13r;
2981 Double_t y0 = 0.5*(xx0*xx0+yy0*yy0-xx0)/yy0;
2982 Double_t r02 = (0.25+y0*y0)*dvertex2;
2983 //curvature (radius) cut
2984 if (r02<r2min) continue;
2988 Double_t c0 = 1/TMath::Sqrt(r02);
2992 //Double_t dfi0 = 2.*TMath::ASin(dvertex*c0*0.5);
2993 //Double_t dfi1 = 2.*TMath::ASin(TMath::Sqrt(yy0*yy0+(1-xx0)*(1-xx0))*dvertex*c0*0.5);
2994 Double_t dfi0 = 2.*AliTPCFastMath::FastAsin(dvertex*c0*0.5);
2995 Double_t dfi1 = 2.*AliTPCFastMath::FastAsin(TMath::Sqrt(yy0*yy0+(1-xx0)*(1-xx0))*dvertex*c0*0.5);
2998 Double_t z0 = kcl->GetZ();
2999 Double_t zzzz2 = z1-(z1-z3)*dfi1/dfi0;
3000 if (TMath::Abs(zzzz2-z0)>0.5) continue;
3003 Double_t dip = (z1-z0)*c0/dfi1;
3004 Double_t x0 = (0.5*cs13+y0*sn13)*dvertex*c0;
3015 x2= xx2*cs-y2*sn*dsec;
3016 y2=+xx2*sn*dsec+y2*cs;
3026 // do we have cluster at the middle ?
3028 GetProlongation(x1,xm,x,ym,zm);
3030 AliTPCclusterMI * cm=0;
3031 if (TMath::Abs(ym)-ymaxm<0){
3032 cm = krm.FindNearest2(ym,zm,1.0,0.6,dummy);
3033 if ((!cm) || (cm->IsUsed(10))) {
3038 // rotate y1 to system 0
3039 // get state vector in rotated system
3040 Double_t yr1 = (-0.5*sn13+y0*cs13)*dvertex*c0;
3041 Double_t xr2 = x0*cs+yr1*sn*dsec;
3042 Double_t xr[5]={kcl->GetY(),kcl->GetZ(), xr2, dip, c0};
3044 GetProlongation(xx2,xm,xr,ym,zm);
3045 if (TMath::Abs(ym)-ymaxm<0){
3046 cm = kr2m.FindNearest2(ym,zm,1.0,0.6,dummy);
3047 if ((!cm) || (cm->IsUsed(10))) {
3057 dym = ym - cm->GetY();
3058 dzm = zm - cm->GetZ();
3065 Double_t sy1=kr1[is]->GetSigmaY2()*2., sz1=kr1[is]->GetSigmaZ2()*2.;
3066 Double_t sy2=kcl->GetSigmaY2()*2., sz2=kcl->GetSigmaZ2()*2.;
3067 //Double_t sy3=400*3./12., sy=0.1, sz=0.1;
3068 Double_t sy3=25000*x[4]*x[4]+0.1, sy=0.1, sz=0.1;
3069 //Double_t sy3=25000*x[4]*x[4]*60+0.5, sy=0.1, sz=0.1;
3071 Double_t f40=(F1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
3072 Double_t f42=(F1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
3073 Double_t f43=(F1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
3074 Double_t f20=(F2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
3075 Double_t f22=(F2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
3076 Double_t f23=(F2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
3078 Double_t f30=(F3(x1,y1+sy,x2,y2,z1,z2)-x[3])/sy;
3079 Double_t f31=(F3(x1,y1,x2,y2,z1+sz,z2)-x[3])/sz;
3080 Double_t f32=(F3(x1,y1,x2,y2+sy,z1,z2)-x[3])/sy;
3081 Double_t f34=(F3(x1,y1,x2,y2,z1,z2+sz)-x[3])/sz;
3085 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3086 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3087 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3088 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3089 c[13]=f30*sy1*f40+f32*sy2*f42;
3090 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3092 // if (!BuildSeed(kr1[is],kcl,0,x1,x2,x3,x,c)) continue;
3094 UInt_t index=kr1.GetIndex(is);
3095 seed->~AliTPCseed(); // this does not set the pointer to 0...
3096 AliTPCseed *track=new(seed) AliTPCseed(x1, ns*alpha+shift, x, c, index);
3098 track->SetIsSeeding(kTRUE);
3099 track->SetSeed1(i1);
3100 track->SetSeed2(i2);
3101 track->SetSeedType(3);
3105 FollowProlongation(*track, (i1+i2)/2,1);
3106 Int_t foundable,found,shared;
3107 track->GetClusterStatistic((i1+i2)/2,i1, found, foundable, shared, kTRUE);
3108 if ((found<0.55*foundable) || shared>0.5*found || (track->GetSigmaY2()+track->GetSigmaZ2())>0.5){
3110 seed->~AliTPCseed();
3116 FollowProlongation(*track, i2,1);
3120 track->SetBConstrain(1);
3121 // track->fLastPoint = i1+fInnerSec->GetNRows(); // first cluster in track position
3122 track->SetLastPoint(i1); // first cluster in track position
3123 track->SetFirstPoint(track->GetLastPoint());
3125 if (track->GetNumberOfClusters()<(i1-i2)*0.5 ||
3126 track->GetNumberOfClusters() < track->GetNFoundable()*0.6 ||
3127 track->GetNShared()>0.4*track->GetNumberOfClusters() ) {
3129 seed->~AliTPCseed();
3133 // Z VERTEX CONDITION
3134 Double_t zv, bz=GetBz();
3135 if ( !track->GetZAt(0.,bz,zv) ) continue;
3136 if (TMath::Abs(zv-z3)>cuts[2]) {
3137 FollowProlongation(*track, TMath::Max(i2-20,0));
3138 if ( !track->GetZAt(0.,bz,zv) ) continue;
3139 if (TMath::Abs(zv-z3)>cuts[2]){
3140 FollowProlongation(*track, TMath::Max(i2-40,0));
3141 if ( !track->GetZAt(0.,bz,zv) ) continue;
3142 if (TMath::Abs(zv-z3)>cuts[2] &&(track->GetNumberOfClusters() > track->GetNFoundable()*0.7)){
3143 // make seed without constrain
3144 AliTPCseed * track2 = MakeSeed(track,0.2,0.5,1.);
3145 FollowProlongation(*track2, i2,1);
3146 track2->SetBConstrain(kFALSE);
3147 track2->SetSeedType(1);
3148 arr->AddLast(track2);
3150 seed->~AliTPCseed();
3155 seed->~AliTPCseed();
3162 track->SetSeedType(0);
3163 arr->AddLast(track);
3164 seed = new AliTPCseed;
3166 // don't consider other combinations
3167 if (track->GetNumberOfClusters() > track->GetNFoundable()*0.8)
3173 Info("MakeSeeds3","\nSeeding statistic:\t%d\t%d\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin,nout1,nout2);
3179 void AliTPCtrackerMI::MakeSeeds5(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, Float_t cuts[4],
3184 //-----------------------------------------------------------------
3185 // This function creates track seeds.
3186 //-----------------------------------------------------------------
3187 // cuts[0] - fP4 cut
3188 // cuts[1] - tan(phi) cut
3189 // cuts[2] - zvertex cut
3190 // cuts[3] - fP3 cut
3200 Double_t x[5], c[15];
3202 // make temporary seed
3203 AliTPCseed * seed = new AliTPCseed;
3204 Double_t alpha=fOuterSec->GetAlpha(), shift=fOuterSec->GetAlphaShift();
3205 // Double_t cs=cos(alpha), sn=sin(alpha);
3210 Double_t x1 = GetXrow(i1-1);
3211 const AliTPCtrackerRow& kr1=GetRow(sec,i1-1);
3212 Double_t y1max = GetMaxY(i1-1)-kr1.GetDeadZone()-1.5;
3214 Double_t x1p = GetXrow(i1);
3215 const AliTPCtrackerRow& kr1p=GetRow(sec,i1);
3217 Double_t x1m = GetXrow(i1-2);
3218 const AliTPCtrackerRow& kr1m=GetRow(sec,i1-2);
3221 //last 3 padrow for seeding
3222 AliTPCtrackerRow& kr3 = GetRow((sec+fkNOS)%fkNOS,i1-7);
3223 Double_t x3 = GetXrow(i1-7);
3224 // Double_t y3max= GetMaxY(i1-7)-kr3.fDeadZone-1.5;
3226 AliTPCtrackerRow& kr3p = GetRow((sec+fkNOS)%fkNOS,i1-6);
3227 Double_t x3p = GetXrow(i1-6);
3229 AliTPCtrackerRow& kr3m = GetRow((sec+fkNOS)%fkNOS,i1-8);
3230 Double_t x3m = GetXrow(i1-8);
3235 Int_t im = i1-4; //middle pad row index
3236 Double_t xm = GetXrow(im); // radius of middle pad-row
3237 const AliTPCtrackerRow& krm=GetRow(sec,im); //middle pad -row
3238 // Double_t ymmax = GetMaxY(im)-kr1.fDeadZone-1.5;
3241 Double_t deltax = x1-x3;
3242 Double_t dymax = deltax*cuts[1];
3243 Double_t dzmax = deltax*cuts[3];
3245 // loop over clusters
3246 for (Int_t is=0; is < kr1; is++) {
3248 if (kr1[is]->IsUsed(10)) continue;
3249 Double_t y1=kr1[is]->GetY(), z1=kr1[is]->GetZ();
3251 if (deltay>0 && TMath::Abs(y1max-TMath::Abs(y1))> deltay ) continue; // seed only at the edge
3253 Int_t index1 = TMath::Max(kr3.Find(z1-dzmax)-1,0);
3254 Int_t index2 = TMath::Min(kr3.Find(z1+dzmax)+1,kr3);
3260 for (Int_t js=index1; js < index2; js++) {
3261 const AliTPCclusterMI *kcl = kr3[js];
3262 if (kcl->IsUsed(10)) continue;
3264 // apply angular cuts
3265 if (TMath::Abs(y1-y3)>dymax) continue;
3268 if (TMath::Abs(z1-z3)>dzmax) continue;
3270 Double_t angley = (y1-y3)/(x1-x3);
3271 Double_t anglez = (z1-z3)/(x1-x3);
3273 Double_t erry = TMath::Abs(angley)*(x1-x1m)*0.5+0.5;
3274 Double_t errz = TMath::Abs(anglez)*(x1-x1m)*0.5+0.5;
3276 Double_t yyym = angley*(xm-x1)+y1;
3277 Double_t zzzm = anglez*(xm-x1)+z1;
3279 const AliTPCclusterMI *kcm = krm.FindNearest2(yyym,zzzm,erry,errz,index);
3281 if (kcm->IsUsed(10)) continue;
3283 erry = TMath::Abs(angley)*(x1-x1m)*0.4+0.5;
3284 errz = TMath::Abs(anglez)*(x1-x1m)*0.4+0.5;
3291 // look around first
3292 const AliTPCclusterMI *kc1m = kr1m.FindNearest2(angley*(x1m-x1)+y1,
3298 if (kc1m->IsUsed(10)) used++;
3300 const AliTPCclusterMI *kc1p = kr1p.FindNearest2(angley*(x1p-x1)+y1,
3306 if (kc1p->IsUsed(10)) used++;
3308 if (used>1) continue;
3309 if (found<1) continue;
3313 const AliTPCclusterMI *kc3m = kr3m.FindNearest2(angley*(x3m-x3)+y3,
3319 if (kc3m->IsUsed(10)) used++;
3323 const AliTPCclusterMI *kc3p = kr3p.FindNearest2(angley*(x3p-x3)+y3,
3329 if (kc3p->IsUsed(10)) used++;
3333 if (used>1) continue;
3334 if (found<3) continue;
3344 x[4]=F1(x1,y1,x2,y2,x3,y3);
3345 //if (TMath::Abs(x[4]) >= cuts[0]) continue;
3348 x[2]=F2(x1,y1,x2,y2,x3,y3);
3351 x[3]=F3n(x1,y1,x2,y2,z1,z2,x[4]);
3352 //if (TMath::Abs(x[3]) > cuts[3]) continue;
3356 Double_t sy1=0.1, sz1=0.1;
3357 Double_t sy2=0.1, sz2=0.1;
3358 Double_t sy3=0.1, sy=0.1, sz=0.1;
3360 Double_t f40=(F1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
3361 Double_t f42=(F1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
3362 Double_t f43=(F1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
3363 Double_t f20=(F2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
3364 Double_t f22=(F2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
3365 Double_t f23=(F2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
3367 Double_t f30=(F3(x1,y1+sy,x2,y2,z1,z2)-x[3])/sy;
3368 Double_t f31=(F3(x1,y1,x2,y2,z1+sz,z2)-x[3])/sz;
3369 Double_t f32=(F3(x1,y1,x2,y2+sy,z1,z2)-x[3])/sy;
3370 Double_t f34=(F3(x1,y1,x2,y2,z1,z2+sz)-x[3])/sz;
3374 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3375 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3376 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3377 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3378 c[13]=f30*sy1*f40+f32*sy2*f42;
3379 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3381 // if (!BuildSeed(kr1[is],kcl,0,x1,x2,x3,x,c)) continue;
3383 UInt_t index=kr1.GetIndex(is);
3384 seed->~AliTPCseed();
3385 AliTPCseed *track=new(seed) AliTPCseed(x1, sec*alpha+shift, x, c, index);
3387 track->SetIsSeeding(kTRUE);
3390 FollowProlongation(*track, i1-7,1);
3391 if (track->GetNumberOfClusters() < track->GetNFoundable()*0.75 ||
3392 track->GetNShared()>0.6*track->GetNumberOfClusters() || ( track->GetSigmaY2()+ track->GetSigmaZ2())>0.6){
3394 seed->~AliTPCseed();
3400 FollowProlongation(*track, i2,1);
3401 track->SetBConstrain(0);
3402 track->SetLastPoint(i1+fInnerSec->GetNRows()); // first cluster in track position
3403 track->SetFirstPoint(track->GetLastPoint());
3405 if (track->GetNumberOfClusters()<(i1-i2)*0.5 ||
3406 track->GetNumberOfClusters()<track->GetNFoundable()*0.7 ||
3407 track->GetNShared()>2. || track->GetChi2()/track->GetNumberOfClusters()>6 || ( track->GetSigmaY2()+ track->GetSigmaZ2())>0.5 ) {
3409 seed->~AliTPCseed();
3414 FollowProlongation(*track, TMath::Max(i2-10,0),1);
3415 AliTPCseed * track2 = MakeSeed(track,0.2,0.5,0.9);
3416 FollowProlongation(*track2, i2,1);
3417 track2->SetBConstrain(kFALSE);
3418 track2->SetSeedType(4);
3419 arr->AddLast(track2);
3421 seed->~AliTPCseed();
3425 //arr->AddLast(track);
3426 //seed = new AliTPCseed;
3432 Info("MakeSeeds5","\nSeeding statiistic:\t%d\t%d\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin,nout1,nout2,nout3);
3438 //_____________________________________________________________________________
3439 void AliTPCtrackerMI::MakeSeeds2(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, Float_t */*cuts[4]*/,
3440 Float_t deltay, Bool_t /*bconstrain*/) {
3441 //-----------------------------------------------------------------
3442 // This function creates track seeds - without vertex constraint
3443 //-----------------------------------------------------------------
3444 // cuts[0] - fP4 cut - not applied
3445 // cuts[1] - tan(phi) cut
3446 // cuts[2] - zvertex cut - not applied
3447 // cuts[3] - fP3 cut
3457 Double_t alpha=fOuterSec->GetAlpha(), shift=fOuterSec->GetAlphaShift();
3458 // Double_t cs=cos(alpha), sn=sin(alpha);
3459 Int_t row0 = (i1+i2)/2;
3460 Int_t drow = (i1-i2)/2;
3461 const AliTPCtrackerRow& kr0=fSectors[sec][row0];
3462 AliTPCtrackerRow * kr=0;
3464 AliTPCpolyTrack polytrack;
3465 Int_t nclusters=fSectors[sec][row0];
3466 AliTPCseed * seed = new AliTPCseed;
3471 for (Int_t is=0; is < nclusters; is++) { //LOOP over clusters
3473 Int_t nfoundable =0;
3474 for (Int_t iter =1; iter<2; iter++){ //iterations
3475 const AliTPCtrackerRow& krm=fSectors[sec][row0-iter];
3476 const AliTPCtrackerRow& krp=fSectors[sec][row0+iter];
3477 const AliTPCclusterMI * cl= kr0[is];
3479 if (cl->IsUsed(10)) {
3485 Double_t x = kr0.GetX();
3486 // Initialization of the polytrack
3491 Double_t y0= cl->GetY();
3492 Double_t z0= cl->GetZ();
3496 Double_t ymax = fSectors->GetMaxY(row0)-kr0.GetDeadZone()-1.5;
3497 if (deltay>0 && TMath::Abs(ymax-TMath::Abs(y0))> deltay ) continue; // seed only at the edge
3499 erry = (0.5)*cl->GetSigmaY2()/TMath::Sqrt(cl->GetQ())*6;
3500 errz = (0.5)*cl->GetSigmaZ2()/TMath::Sqrt(cl->GetQ())*6;
3501 polytrack.AddPoint(x,y0,z0,erry, errz);
3504 if (cl->IsUsed(10)) sumused++;
3507 Float_t roady = (5*TMath::Sqrt(cl->GetSigmaY2()+0.2)+1.)*iter;
3508 Float_t roadz = (5*TMath::Sqrt(cl->GetSigmaZ2()+0.2)+1.)*iter;
3511 AliTPCclusterMI * cl1 = krm.FindNearest(y0,z0,roady,roadz);
3512 if (cl1 && TMath::Abs(ymax-TMath::Abs(y0))) {
3513 erry = (0.5)*cl1->GetSigmaY2()/TMath::Sqrt(cl1->GetQ())*3;
3514 errz = (0.5)*cl1->GetSigmaZ2()/TMath::Sqrt(cl1->GetQ())*3;
3515 if (cl1->IsUsed(10)) sumused++;
3516 polytrack.AddPoint(x,cl1->GetY(),cl1->GetZ(),erry,errz);
3520 AliTPCclusterMI * cl2 = krp.FindNearest(y0,z0,roady,roadz);
3522 erry = (0.5)*cl2->GetSigmaY2()/TMath::Sqrt(cl2->GetQ())*3;
3523 errz = (0.5)*cl2->GetSigmaZ2()/TMath::Sqrt(cl2->GetQ())*3;
3524 if (cl2->IsUsed(10)) sumused++;
3525 polytrack.AddPoint(x,cl2->GetY(),cl2->GetZ(),erry,errz);
3528 if (sumused>0) continue;
3530 polytrack.UpdateParameters();
3536 nfoundable = polytrack.GetN();
3537 nfound = nfoundable;
3539 for (Int_t ddrow = iter+1; ddrow<drow;ddrow++){
3540 Float_t maxdist = 0.8*(1.+3./(ddrow));
3541 for (Int_t delta = -1;delta<=1;delta+=2){
3542 Int_t row = row0+ddrow*delta;
3543 kr = &(fSectors[sec][row]);
3544 Double_t xn = kr->GetX();
3545 Double_t ymax = fSectors->GetMaxY(row)-kr->GetDeadZone()-1.5;
3546 polytrack.GetFitPoint(xn,yn,zn);
3547 if (TMath::Abs(yn)>ymax) continue;
3549 AliTPCclusterMI * cln = kr->FindNearest(yn,zn,roady,roadz);
3551 Float_t dist = TMath::Sqrt( (yn-cln->GetY())*(yn-cln->GetY())+(zn-cln->GetZ())*(zn-cln->GetZ()));
3554 erry = (dist+0.3)*cln->GetSigmaY2()/TMath::Sqrt(cln->GetQ())*(1.+1./(ddrow));
3555 errz = (dist+0.3)*cln->GetSigmaZ2()/TMath::Sqrt(cln->GetQ())*(1.+1./(ddrow));
3556 if (cln->IsUsed(10)) {
3557 // printf("used\n");
3565 polytrack.AddPoint(xn,cln->GetY(),cln->GetZ(),erry, errz);
3570 if ( (sumused>3) || (sumused>0.5*nfound) || (nfound<0.6*nfoundable)) break;
3571 polytrack.UpdateParameters();
3574 if ( (sumused>3) || (sumused>0.5*nfound)) {
3575 //printf("sumused %d\n",sumused);
3580 polytrack.GetFitDerivation(kr0.GetX(),dy,dz);
3581 AliTPCpolyTrack track2;
3583 polytrack.Refit(track2,0.5+TMath::Abs(dy)*0.3,0.4+TMath::Abs(dz)*0.3);
3584 if (track2.GetN()<0.5*nfoundable) continue;
3587 if ((nfound>0.6*nfoundable) &&( nfoundable>0.4*(i1-i2))) {
3589 // test seed with and without constrain
3590 for (Int_t constrain=0; constrain<=0;constrain++){
3591 // add polytrack candidate
3593 Double_t x[5], c[15];
3594 Double_t x1,x2,x3,y1,y2,y3,z1,z2,z3;
3595 track2.GetBoundaries(x3,x1);
3597 track2.GetFitPoint(x1,y1,z1);
3598 track2.GetFitPoint(x2,y2,z2);
3599 track2.GetFitPoint(x3,y3,z3);
3601 //is track pointing to the vertex ?
3604 polytrack.GetFitPoint(x0,y0,z0);
3617 x[4]=F1(x1,y1,x2,y2,x3,y3);
3619 // if (TMath::Abs(x[4]) >= cuts[0]) continue; //
3620 x[2]=F2(x1,y1,x2,y2,x3,y3);
3622 //if (TMath::Abs(x[4]*x1-x[2]) >= cuts[1]) continue;
3623 //x[3]=F3(x1,y1,x2,y2,z1,z2);
3624 x[3]=F3n(x1,y1,x3,y3,z1,z3,x[4]);
3625 //if (TMath::Abs(x[3]) > cuts[3]) continue;
3628 Double_t sy =0.1, sz =0.1;
3629 Double_t sy1=0.02, sz1=0.02;
3630 Double_t sy2=0.02, sz2=0.02;
3634 sy3=25000*x[4]*x[4]+0.1, sy=0.1, sz=0.1;
3637 Double_t f40=(F1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
3638 Double_t f42=(F1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
3639 Double_t f43=(F1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
3640 Double_t f20=(F2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
3641 Double_t f22=(F2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
3642 Double_t f23=(F2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
3644 Double_t f30=(F3(x1,y1+sy,x3,y3,z1,z3)-x[3])/sy;
3645 Double_t f31=(F3(x1,y1,x3,y3,z1+sz,z3)-x[3])/sz;
3646 Double_t f32=(F3(x1,y1,x3,y3+sy,z1,z3)-x[3])/sy;
3647 Double_t f34=(F3(x1,y1,x3,y3,z1,z3+sz)-x[3])/sz;
3652 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3653 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3654 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3655 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3656 c[13]=f30*sy1*f40+f32*sy2*f42;
3657 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3659 //Int_t row1 = fSectors->GetRowNumber(x1);
3660 Int_t row1 = GetRowNumber(x1);
3664 seed->~AliTPCseed();
3665 AliTPCseed *track=new(seed) AliTPCseed(x1,sec*alpha+shift,x,c,index);
3666 track->SetIsSeeding(kTRUE);
3667 Int_t rc=FollowProlongation(*track, i2);
3668 if (constrain) track->SetBConstrain(1);
3670 track->SetBConstrain(0);
3671 track->SetLastPoint(row1+fInnerSec->GetNRows()); // first cluster in track position
3672 track->SetFirstPoint(track->GetLastPoint());
3674 if (rc==0 || track->GetNumberOfClusters()<(i1-i2)*0.5 ||
3675 track->GetNumberOfClusters() < track->GetNFoundable()*0.6 ||
3676 track->GetNShared()>0.4*track->GetNumberOfClusters()) {
3679 seed->~AliTPCseed();
3682 arr->AddLast(track);
3683 seed = new AliTPCseed;
3687 } // if accepted seed
3690 Info("MakeSeeds2","\nSeeding statiistic:\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin3);
3696 AliTPCseed *AliTPCtrackerMI::MakeSeed(AliTPCseed *track, Float_t r0, Float_t r1, Float_t r2)
3700 //reseed using track points
3701 Int_t p0 = int(r0*track->GetNumberOfClusters()); // point 0
3702 Int_t p1 = int(r1*track->GetNumberOfClusters());
3703 Int_t p2 = int(r2*track->GetNumberOfClusters()); // last point
3705 Double_t x0[3],x1[3],x2[3];
3706 for (Int_t i=0;i<3;i++){
3712 // find track position at given ratio of the length
3713 Int_t sec0=0, sec1=0, sec2=0;
3716 for (Int_t i=0;i<160;i++){
3717 if (track->GetClusterPointer(i)){
3719 AliTPCTrackerPoint *trpoint =track->GetTrackPoint(i);
3720 if ( (index<p0) || x0[0]<0 ){
3721 if (trpoint->GetX()>1){
3722 clindex = track->GetClusterIndex2(i);
3724 x0[0] = trpoint->GetX();
3725 x0[1] = trpoint->GetY();
3726 x0[2] = trpoint->GetZ();
3727 sec0 = ((clindex&0xff000000)>>24)%18;
3732 if ( (index<p1) &&(trpoint->GetX()>1)){
3733 clindex = track->GetClusterIndex2(i);
3735 x1[0] = trpoint->GetX();
3736 x1[1] = trpoint->GetY();
3737 x1[2] = trpoint->GetZ();
3738 sec1 = ((clindex&0xff000000)>>24)%18;
3741 if ( (index<p2) &&(trpoint->GetX()>1)){
3742 clindex = track->GetClusterIndex2(i);
3744 x2[0] = trpoint->GetX();
3745 x2[1] = trpoint->GetY();
3746 x2[2] = trpoint->GetZ();
3747 sec2 = ((clindex&0xff000000)>>24)%18;
3754 Double_t alpha, cs,sn, xx2,yy2;
3756 alpha = (sec1-sec2)*fSectors->GetAlpha();
3757 cs = TMath::Cos(alpha);
3758 sn = TMath::Sin(alpha);
3759 xx2= x1[0]*cs-x1[1]*sn;
3760 yy2= x1[0]*sn+x1[1]*cs;
3764 alpha = (sec0-sec2)*fSectors->GetAlpha();
3765 cs = TMath::Cos(alpha);
3766 sn = TMath::Sin(alpha);
3767 xx2= x0[0]*cs-x0[1]*sn;
3768 yy2= x0[0]*sn+x0[1]*cs;
3774 Double_t x[5],c[15];
3778 x[4]=F1(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]);
3779 // if (x[4]>1) return 0;
3780 x[2]=F2(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]);
3781 x[3]=F3n(x2[0],x2[1],x0[0],x0[1],x2[2],x0[2],x[4]);
3782 //if (TMath::Abs(x[3]) > 2.2) return 0;
3783 //if (TMath::Abs(x[2]) > 1.99) return 0;
3785 Double_t sy =0.1, sz =0.1;
3787 Double_t sy1=0.02+track->GetSigmaY2(), sz1=0.02+track->GetSigmaZ2();
3788 Double_t sy2=0.01+track->GetSigmaY2(), sz2=0.01+track->GetSigmaZ2();
3789 Double_t sy3=0.01+track->GetSigmaY2();
3791 Double_t f40=(F1(x2[0],x2[1]+sy,x1[0],x1[1],x0[0],x0[1])-x[4])/sy;
3792 Double_t f42=(F1(x2[0],x2[1],x1[0],x1[1]+sy,x0[0],x0[1])-x[4])/sy;
3793 Double_t f43=(F1(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]+sy)-x[4])/sy;
3794 Double_t f20=(F2(x2[0],x2[1]+sy,x1[0],x1[1],x0[0],x0[1])-x[2])/sy;
3795 Double_t f22=(F2(x2[0],x2[1],x1[0],x1[1]+sy,x0[0],x0[1])-x[2])/sy;
3796 Double_t f23=(F2(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]+sy)-x[2])/sy;
3798 Double_t f30=(F3(x2[0],x2[1]+sy,x0[0],x0[1],x2[2],x0[2])-x[3])/sy;
3799 Double_t f31=(F3(x2[0],x2[1],x0[0],x0[1],x2[2]+sz,x0[2])-x[3])/sz;
3800 Double_t f32=(F3(x2[0],x2[1],x0[0],x0[1]+sy,x2[2],x0[2])-x[3])/sy;
3801 Double_t f34=(F3(x2[0],x2[1],x0[0],x0[1],x2[2],x0[2]+sz)-x[3])/sz;
3806 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3807 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3808 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3809 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3810 c[13]=f30*sy1*f40+f32*sy2*f42;
3811 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3813 // Int_t row1 = fSectors->GetRowNumber(x2[0]);
3814 AliTPCseed *seed=new AliTPCseed(x2[0], sec2*fSectors->GetAlpha()+fSectors->GetAlphaShift(), x, c, 0);
3815 // Double_t y0,z0,y1,z1, y2,z2;
3816 //seed->GetProlongation(x0[0],y0,z0);
3817 // seed->GetProlongation(x1[0],y1,z1);
3818 //seed->GetProlongation(x2[0],y2,z2);
3820 seed->SetLastPoint(pp2);
3821 seed->SetFirstPoint(pp2);
3828 AliTPCseed *AliTPCtrackerMI::ReSeed(AliTPCseed *track, Float_t r0, Float_t r1, Float_t r2)
3832 //reseed using founded clusters
3834 // Find the number of clusters
3835 Int_t nclusters = 0;
3836 for (Int_t irow=0;irow<160;irow++){
3837 if (track->GetClusterIndex(irow)>0) nclusters++;
3841 ipos[0] = TMath::Max(int(r0*nclusters),0); // point 0 cluster
3842 ipos[1] = TMath::Min(int(r1*nclusters),nclusters-1); //
3843 ipos[2] = TMath::Min(int(r2*nclusters),nclusters-1); // last point
3847 Int_t row[3],sec[3]={0,0,0};
3849 // find track row position at given ratio of the length
3851 for (Int_t irow=0;irow<160;irow++){
3852 if (track->GetClusterIndex2(irow)<0) continue;
3854 for (Int_t ipoint=0;ipoint<3;ipoint++){
3855 if (index<=ipos[ipoint]) row[ipoint] = irow;
3859 //Get cluster and sector position
3860 for (Int_t ipoint=0;ipoint<3;ipoint++){
3861 Int_t clindex = track->GetClusterIndex2(row[ipoint]);
3862 AliTPCclusterMI * cl = GetClusterMI(clindex);
3865 // AliTPCclusterMI * cl = GetClusterMI(clindex);
3868 sec[ipoint] = ((clindex&0xff000000)>>24)%18;
3869 xyz[ipoint][0] = GetXrow(row[ipoint]);
3870 xyz[ipoint][1] = cl->GetY();
3871 xyz[ipoint][2] = cl->GetZ();
3875 // Calculate seed state vector and covariance matrix
3877 Double_t alpha, cs,sn, xx2,yy2;
3879 alpha = (sec[1]-sec[2])*fSectors->GetAlpha();
3880 cs = TMath::Cos(alpha);
3881 sn = TMath::Sin(alpha);
3882 xx2= xyz[1][0]*cs-xyz[1][1]*sn;
3883 yy2= xyz[1][0]*sn+xyz[1][1]*cs;
3887 alpha = (sec[0]-sec[2])*fSectors->GetAlpha();
3888 cs = TMath::Cos(alpha);
3889 sn = TMath::Sin(alpha);
3890 xx2= xyz[0][0]*cs-xyz[0][1]*sn;
3891 yy2= xyz[0][0]*sn+xyz[0][1]*cs;
3897 Double_t x[5],c[15];
3901 x[4]=F1(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
3902 x[2]=F2(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
3903 x[3]=F3n(xyz[2][0],xyz[2][1],xyz[0][0],xyz[0][1],xyz[2][2],xyz[0][2],x[4]);
3905 Double_t sy =0.1, sz =0.1;
3907 Double_t sy1=0.2, sz1=0.2;
3908 Double_t sy2=0.2, sz2=0.2;
3911 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;
3912 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;
3913 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;
3914 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;
3915 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;
3916 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;
3918 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;
3919 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;
3920 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;
3921 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;
3926 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3927 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3928 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3929 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3930 c[13]=f30*sy1*f40+f32*sy2*f42;
3931 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3933 // Int_t row1 = fSectors->GetRowNumber(xyz[2][0]);
3934 AliTPCseed *seed=new AliTPCseed(xyz[2][0], sec[2]*fSectors->GetAlpha()+fSectors->GetAlphaShift(), x, c, 0);
3935 seed->SetLastPoint(row[2]);
3936 seed->SetFirstPoint(row[2]);
3941 AliTPCseed *AliTPCtrackerMI::ReSeed(AliTPCseed *track,Int_t r0, Bool_t forward)
3945 //reseed using founded clusters
3948 Int_t row[3]={0,0,0};
3949 Int_t sec[3]={0,0,0};
3951 // forward direction
3953 for (Int_t irow=r0;irow<160;irow++){
3954 if (track->GetClusterIndex(irow)>0){
3959 for (Int_t irow=160;irow>r0;irow--){
3960 if (track->GetClusterIndex(irow)>0){
3965 for (Int_t irow=row[2]-15;irow>row[0];irow--){
3966 if (track->GetClusterIndex(irow)>0){
3974 for (Int_t irow=0;irow<r0;irow++){
3975 if (track->GetClusterIndex(irow)>0){
3980 for (Int_t irow=r0;irow>0;irow--){
3981 if (track->GetClusterIndex(irow)>0){
3986 for (Int_t irow=row[2]-15;irow>row[0];irow--){
3987 if (track->GetClusterIndex(irow)>0){
3994 if ((row[2]-row[0])<20) return 0;
3995 if (row[1]==0) return 0;
3998 //Get cluster and sector position
3999 for (Int_t ipoint=0;ipoint<3;ipoint++){
4000 Int_t clindex = track->GetClusterIndex2(row[ipoint]);
4001 AliTPCclusterMI * cl = GetClusterMI(clindex);
4004 // AliTPCclusterMI * cl = GetClusterMI(clindex);
4007 sec[ipoint] = ((clindex&0xff000000)>>24)%18;
4008 xyz[ipoint][0] = GetXrow(row[ipoint]);
4009 AliTPCTrackerPoint * point = track->GetTrackPoint(row[ipoint]);
4010 if (point&&ipoint<2){
4012 xyz[ipoint][1] = point->GetY();
4013 xyz[ipoint][2] = point->GetZ();
4016 xyz[ipoint][1] = cl->GetY();
4017 xyz[ipoint][2] = cl->GetZ();
4024 // Calculate seed state vector and covariance matrix
4026 Double_t alpha, cs,sn, xx2,yy2;
4028 alpha = (sec[1]-sec[2])*fSectors->GetAlpha();
4029 cs = TMath::Cos(alpha);
4030 sn = TMath::Sin(alpha);
4031 xx2= xyz[1][0]*cs-xyz[1][1]*sn;
4032 yy2= xyz[1][0]*sn+xyz[1][1]*cs;
4036 alpha = (sec[0]-sec[2])*fSectors->GetAlpha();
4037 cs = TMath::Cos(alpha);
4038 sn = TMath::Sin(alpha);
4039 xx2= xyz[0][0]*cs-xyz[0][1]*sn;
4040 yy2= xyz[0][0]*sn+xyz[0][1]*cs;
4046 Double_t x[5],c[15];
4050 x[4]=F1(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
4051 x[2]=F2(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
4052 x[3]=F3n(xyz[2][0],xyz[2][1],xyz[0][0],xyz[0][1],xyz[2][2],xyz[0][2],x[4]);
4054 Double_t sy =0.1, sz =0.1;
4056 Double_t sy1=0.2, sz1=0.2;
4057 Double_t sy2=0.2, sz2=0.2;
4060 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;
4061 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;
4062 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;
4063 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;
4064 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;
4065 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;
4067 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;
4068 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;
4069 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;
4070 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;
4075 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
4076 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
4077 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
4078 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
4079 c[13]=f30*sy1*f40+f32*sy2*f42;
4080 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
4082 // Int_t row1 = fSectors->GetRowNumber(xyz[2][0]);
4083 AliTPCseed *seed=new AliTPCseed(xyz[2][0], sec[2]*fSectors->GetAlpha()+fSectors->GetAlphaShift(), x, c, 0);
4084 seed->SetLastPoint(row[2]);
4085 seed->SetFirstPoint(row[2]);
4086 for (Int_t i=row[0];i<row[2];i++){
4087 seed->SetClusterIndex(i, track->GetClusterIndex(i));
4095 void AliTPCtrackerMI::FindMultiMC(TObjArray * array, AliESDEvent */*esd*/, Int_t iter)
4098 // find multi tracks - THIS FUNCTION IS ONLY FOR DEBUG PURPOSES
4100 // Use AliTPCReconstructor::StreamLevel()>2 if you want to tune parameters - cuts
4102 // Two reasons to have multiple find tracks
4103 // 1. Curling tracks can be find more than once
4104 // 2. Splitted tracks
4105 // a.) Multiple seeding to increase tracking efficiency - (~ 100% reached)
4106 // b.) Edge effect on the sector boundaries
4109 // Algorithm done in 2 phases - because of CPU consumption
4110 // it is n^2 algorithm - for lead-lead 20000x20000 combination are investigated
4112 // Algorihm for curling tracks sign:
4113 // 1 phase -makes a very rough fast cuts to minimize combinatorics
4114 // a.) opposite sign
4115 // b.) one of the tracks - not pointing to the primary vertex -
4116 // c.) delta tan(theta)
4118 // 2 phase - calculates DCA between tracks - time consument
4123 // General cuts - for splitted tracks and for curling tracks
4125 const Float_t kMaxdPhi = 0.2; // maximal distance in phi
4127 // Curling tracks cuts
4132 Int_t nentries = array->GetEntriesFast();
4133 AliHelix *helixes = new AliHelix[nentries];
4134 Float_t *xm = new Float_t[nentries];
4135 Float_t *dz0 = new Float_t[nentries];
4136 Float_t *dz1 = new Float_t[nentries];
4142 // Find track COG in x direction - point with best defined parameters
4144 for (Int_t i=0;i<nentries;i++){
4145 AliTPCseed* track = (AliTPCseed*)array->At(i);
4146 if (!track) continue;
4147 track->SetCircular(0);
4148 new (&helixes[i]) AliHelix(*track);
4152 track->GetDZ(GetX(),GetY(),GetZ(),GetBz(),dz);
4155 for (Int_t icl=0; icl<160; icl++){
4156 AliTPCclusterMI * cl = track->GetClusterPointer(icl);
4162 if (ncl>0) xm[i]/=Float_t(ncl);
4165 for (Int_t i0=0;i0<nentries;i0++){
4166 AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
4167 if (!track0) continue;
4168 Float_t xc0 = helixes[i0].GetHelix(6);
4169 Float_t yc0 = helixes[i0].GetHelix(7);
4170 Float_t r0 = helixes[i0].GetHelix(8);
4171 Float_t rc0 = TMath::Sqrt(xc0*xc0+yc0*yc0);
4172 Float_t fi0 = TMath::ATan2(yc0,xc0);
4174 for (Int_t i1=i0+1;i1<nentries;i1++){
4175 AliTPCseed * track1 = (AliTPCseed*)array->At(i1);
4176 if (!track1) continue;
4177 Int_t lab0=track0->GetLabel();
4178 Int_t lab1=track1->GetLabel();
4179 if (TMath::Abs(lab0)!=TMath::Abs(lab1)) continue;
4181 Float_t xc1 = helixes[i1].GetHelix(6);
4182 Float_t yc1 = helixes[i1].GetHelix(7);
4183 Float_t r1 = helixes[i1].GetHelix(8);
4184 Float_t rc1 = TMath::Sqrt(xc1*xc1+yc1*yc1);
4185 Float_t fi1 = TMath::ATan2(yc1,xc1);
4187 Float_t dfi = fi0-fi1;
4190 if (dfi>1.5*TMath::Pi()) dfi-=TMath::Pi(); // take care about edge effect
4191 if (dfi<-1.5*TMath::Pi()) dfi+=TMath::Pi(); //
4192 if (TMath::Abs(dfi)>kMaxdPhi&&helixes[i0].GetHelix(4)*helixes[i1].GetHelix(4)<0){
4194 // if short tracks with undefined sign
4195 fi1 = -TMath::ATan2(yc1,-xc1);
4198 Float_t dtheta = TMath::Abs(track0->GetTgl()-track1->GetTgl())<TMath::Abs(track0->GetTgl()+track1->GetTgl())? track0->GetTgl()-track1->GetTgl():track0->GetTgl()+track1->GetTgl();
4201 // debug stream to tune "fast cuts"
4203 Double_t dist[3]; // distance at X
4204 Double_t mdist[3]={0,0,0}; // mean distance X+-40cm
4205 track0->GetDistance(track1,0.5*(xm[i0]+xm[i1])-40.,dist,AliTracker::GetBz());
4206 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[i]);
4207 track0->GetDistance(track1,0.5*(xm[i0]+xm[i1])+40.,dist,AliTracker::GetBz());
4208 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[i]);
4209 track0->GetDistance(track1,0.5*(xm[i0]+xm[i1]),dist,AliTracker::GetBz());
4210 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[i]);
4211 for (Int_t i=0;i<3;i++) mdist[i]*=0.33333;
4215 for (Int_t icl=0; icl<160; icl++){
4216 AliTPCclusterMI * cl0 = track0->GetClusterPointer(icl);
4217 AliTPCclusterMI * cl1 = track1->GetClusterPointer(icl);
4220 if (cl0==cl1) sums++;
4224 if (AliTPCReconstructor::StreamLevel()>0) {
4225 TTreeSRedirector &cstream = *fDebugStreamer;
4230 "Tr0.="<<track0<< // seed0
4231 "Tr1.="<<track1<< // seed1
4232 "h0.="<<&helixes[i0]<<
4233 "h1.="<<&helixes[i1]<<
4235 "sum="<<sum<< //the sum of rows with cl in both
4236 "sums="<<sums<< //the sum of shared clusters
4237 "xm0="<<xm[i0]<< // the center of track
4238 "xm1="<<xm[i1]<< // the x center of track
4239 // General cut variables
4240 "dfi="<<dfi<< // distance in fi angle
4241 "dtheta="<<dtheta<< // distance int theta angle
4247 "dist0="<<dist[0]<< //distance x
4248 "dist1="<<dist[1]<< //distance y
4249 "dist2="<<dist[2]<< //distance z
4250 "mdist0="<<mdist[0]<< //distance x
4251 "mdist1="<<mdist[1]<< //distance y
4252 "mdist2="<<mdist[2]<< //distance z
4266 if (AliTPCReconstructor::StreamLevel()>1) {
4267 AliInfo("Time for curling tracks removal DEBUGGING MC");
4273 void AliTPCtrackerMI::FindSplitted(TObjArray * array, AliESDEvent */*esd*/, Int_t iter)
4277 // Two reasons to have multiple find tracks
4278 // 1. Curling tracks can be find more than once
4279 // 2. Splitted tracks
4280 // a.) Multiple seeding to increase tracking efficiency - (~ 100% reached)
4281 // b.) Edge effect on the sector boundaries
4283 // This function tries to find tracks closed in the parametric space
4285 // cut logic if distance is bigger than cut continue - Do Nothing
4286 const Float_t kMaxdTheta = 0.05; // maximal distance in theta
4287 const Float_t kMaxdPhi = 0.05; // maximal deistance in phi
4288 const Float_t kdelta = 40.; //delta r to calculate track distance
4290 // const Float_t kMaxDist0 = 1.; // maximal distance 0
4291 //const Float_t kMaxDist1 = 0.3; // maximal distance 1 - cut if track in separate rows
4294 TCut csec("csec","abs(Tr0.fRelativeSector-Tr1.fRelativeSector)<2");
4295 TCut cdtheta("cdtheta","abs(dtheta)<0.05");
4300 Int_t nentries = array->GetEntriesFast();
4301 AliHelix *helixes = new AliHelix[nentries];
4302 Float_t *xm = new Float_t[nentries];
4308 //Sort tracks according quality
4310 Int_t nseed = array->GetEntriesFast();
4311 Float_t * quality = new Float_t[nseed];
4312 Int_t * indexes = new Int_t[nseed];
4313 for (Int_t i=0; i<nseed; i++) {
4314 AliTPCseed *pt=(AliTPCseed*)array->UncheckedAt(i);
4319 pt->UpdatePoints(); //select first last max dens points
4320 Float_t * points = pt->GetPoints();
4321 if (points[3]<0.8) quality[i] =-1;
4322 quality[i] = (points[2]-points[0])+pt->GetNumberOfClusters();
4323 //prefer high momenta tracks if overlaps
4324 quality[i] *= TMath::Sqrt(TMath::Abs(pt->Pt())+0.5);
4326 TMath::Sort(nseed,quality,indexes);
4330 // Find track COG in x direction - point with best defined parameters
4332 for (Int_t i=0;i<nentries;i++){
4333 AliTPCseed* track = (AliTPCseed*)array->At(i);
4334 if (!track) continue;
4335 track->SetCircular(0);
4336 new (&helixes[i]) AliHelix(*track);
4339 for (Int_t icl=0; icl<160; icl++){
4340 AliTPCclusterMI * cl = track->GetClusterPointer(icl);
4346 if (ncl>0) xm[i]/=Float_t(ncl);
4349 for (Int_t is0=0;is0<nentries;is0++){
4350 Int_t i0 = indexes[is0];
4351 AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
4352 if (!track0) continue;
4353 if (track0->GetKinkIndexes()[0]!=0) continue;
4354 Float_t xc0 = helixes[i0].GetHelix(6);
4355 Float_t yc0 = helixes[i0].GetHelix(7);
4356 Float_t fi0 = TMath::ATan2(yc0,xc0);
4358 for (Int_t is1=is0+1;is1<nentries;is1++){
4359 Int_t i1 = indexes[is1];
4360 AliTPCseed * track1 = (AliTPCseed*)array->At(i1);
4361 if (!track1) continue;
4363 if (TMath::Abs(track0->GetRelativeSector()-track1->GetRelativeSector())>1) continue;
4364 if (track1->GetKinkIndexes()[0]>0 &&track0->GetKinkIndexes()[0]<0) continue;
4365 if (track1->GetKinkIndexes()[0]!=0) continue;
4367 Float_t dtheta = TMath::Abs(track0->GetTgl()-track1->GetTgl())<TMath::Abs(track0->GetTgl()+track1->GetTgl())? track0->GetTgl()-track1->GetTgl():track0->GetTgl()+track1->GetTgl();
4368 if (TMath::Abs(dtheta)>kMaxdTheta) continue;
4370 Float_t xc1 = helixes[i1].GetHelix(6);
4371 Float_t yc1 = helixes[i1].GetHelix(7);
4372 Float_t fi1 = TMath::ATan2(yc1,xc1);
4374 Float_t dfi = fi0-fi1;
4375 if (dfi>1.5*TMath::Pi()) dfi-=TMath::Pi(); // take care about edge effect
4376 if (dfi<-1.5*TMath::Pi()) dfi+=TMath::Pi(); //
4377 if (TMath::Abs(dfi)>kMaxdPhi&&helixes[i0].GetHelix(4)*helixes[i1].GetHelix(4)<0){
4379 // if short tracks with undefined sign
4380 fi1 = -TMath::ATan2(yc1,-xc1);
4383 if (TMath::Abs(dfi)>kMaxdPhi) continue;
4390 for (Int_t icl=0; icl<160; icl++){
4391 Int_t index0=track0->GetClusterIndex2(icl);
4392 Int_t index1=track1->GetClusterIndex2(icl);
4393 Bool_t used0 = (index0>0 && !(index0&0x8000));
4394 Bool_t used1 = (index1>0 && !(index1&0x8000));
4396 if (used0) sum0++; // used cluster0
4397 if (used1) sum1++; // used clusters1
4398 if (used0&&used1) sum++;
4399 if (index0==index1 && used0 && used1) sums++;
4403 if (sums<10) continue;
4404 if (sum<40) continue;
4405 if (sums/Float_t(TMath::Min(sum0,sum1))<0.5) continue;
4407 Double_t dist[5][4]; // distance at X
4408 Double_t mdist[4]={0,0,0,0}; // mean distance on range +-delta
4412 track0->GetDistance(track1,xm[i0],dist[0],AliTracker::GetBz());
4413 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[0][i]);
4414 track0->GetDistance(track1,xm[i1],dist[1],AliTracker::GetBz());
4415 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[1][i]);
4417 track0->GetDistance(track1,TMath::Min(xm[i1],xm[i0])-kdelta,dist[2],AliTracker::GetBz());
4418 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[2][i]);
4419 track0->GetDistance(track1,TMath::Max(xm[i1],xm[i0])+kdelta,dist[3],AliTracker::GetBz());
4420 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[3][i]);
4422 track0->GetDistance(track1,(xm[i1]+xm[i0])*0.5,dist[4],AliTracker::GetBz());
4423 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[4][i]);
4424 for (Int_t i=0;i<3;i++) mdist[i]*=0.2;
4427 Int_t lab0=track0->GetLabel();
4428 Int_t lab1=track1->GetLabel();
4429 if( AliTPCReconstructor::StreamLevel()>5){
4430 TTreeSRedirector &cstream = *fDebugStreamer;
4431 cstream<<"Splitted"<<
4435 "Tr0.="<<track0<< // seed0
4436 "Tr1.="<<track1<< // seed1
4437 "h0.="<<&helixes[i0]<<
4438 "h1.="<<&helixes[i1]<<
4440 "sum="<<sum<< //the sum of rows with cl in both
4441 "sum0="<<sum0<< //the sum of rows with cl in 0 track
4442 "sum1="<<sum1<< //the sum of rows with cl in 1 track
4443 "sums="<<sums<< //the sum of shared clusters
4444 "xm0="<<xm[i0]<< // the center of track
4445 "xm1="<<xm[i1]<< // the x center of track
4446 // General cut variables
4447 "dfi="<<dfi<< // distance in fi angle
4448 "dtheta="<<dtheta<< // distance int theta angle
4451 "dist0="<<dist[4][0]<< //distance x
4452 "dist1="<<dist[4][1]<< //distance y
4453 "dist2="<<dist[4][2]<< //distance z
4454 "mdist0="<<mdist[0]<< //distance x
4455 "mdist1="<<mdist[1]<< //distance y
4456 "mdist2="<<mdist[2]<< //distance z
4459 delete array->RemoveAt(i1);
4466 AliInfo("Time for splitted tracks removal");
4472 void AliTPCtrackerMI::FindCurling(TObjArray * array, AliESDEvent */*esd*/, Int_t iter)
4475 // find Curling tracks
4476 // Use AliTPCReconstructor::StreamLevel()>1 if you want to tune parameters - cuts
4479 // Algorithm done in 2 phases - because of CPU consumption
4480 // it is n^2 algorithm - for lead-lead 20000x20000 combination are investigated
4481 // see detal in MC part what can be used to cut
4485 const Float_t kMaxC = 400; // maximal curvature to of the track
4486 const Float_t kMaxdTheta = 0.15; // maximal distance in theta
4487 const Float_t kMaxdPhi = 0.15; // maximal distance in phi
4488 const Float_t kPtRatio = 0.3; // ratio between pt
4489 const Float_t kMinDCAR = 2.; // distance to the primary vertex in r - see cpipe cut
4492 // Curling tracks cuts
4495 const Float_t kMaxDeltaRMax = 40; // distance in outer radius
4496 const Float_t kMaxDeltaRMin = 5.; // distance in lower radius - see cpipe cut
4497 const Float_t kMinAngle = 2.9; // angle between tracks
4498 const Float_t kMaxDist = 5; // biggest distance
4500 // The cuts can be tuned using the "MC information stored in Multi tree ==> see FindMultiMC
4503 TCut csign("csign","Tr0.fP[4]*Tr1.fP[4]<0"); //opposite sign
4504 TCut cmax("cmax","abs(Tr0.GetC())>1/400");
4505 TCut cda("cda","sqrt(dtheta^2+dfi^2)<0.15");
4506 TCut ccratio("ccratio","abs((Tr0.fP[4]+Tr1.fP[4])/(abs(Tr0.fP[4])+abs(Tr1.fP[4])))<0.3");
4507 TCut cpipe("cpipe", "min(abs(r0-rc0),abs(r1-rc1))>5");
4509 TCut cdrmax("cdrmax","abs(abs(rc0+r0)-abs(rc1+r1))<40")
4510 TCut cdrmin("cdrmin","abs(abs(rc0+r0)-abs(rc1+r1))<10")
4512 Multi->Draw("dfi","iter==0"+csign+cmax+cda+ccratio); ~94% of curling tracks fulfill
4513 Multi->Draw("min(abs(r0-rc0),abs(r1-rc1))","iter==0&&abs(lab1)==abs(lab0)"+csign+cmax+cda+ccratio+cpipe+cdrmin+cdrmax); //80%
4515 Curling2->Draw("dfi","iter==0&&abs(lab0)==abs(lab1)"+csign+cmax+cdtheta+cdfi+ccratio)
4521 Int_t nentries = array->GetEntriesFast();
4522 AliHelix *helixes = new AliHelix[nentries];
4523 for (Int_t i=0;i<nentries;i++){
4524 AliTPCseed* track = (AliTPCseed*)array->At(i);
4525 if (!track) continue;
4526 track->SetCircular(0);
4527 new (&helixes[i]) AliHelix(*track);
4533 Double_t phase[2][2],radius[2];
4538 for (Int_t i0=0;i0<nentries;i0++){
4539 AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
4540 if (!track0) continue;
4541 if (TMath::Abs(track0->GetC())<1/kMaxC) continue;
4542 Float_t xc0 = helixes[i0].GetHelix(6);
4543 Float_t yc0 = helixes[i0].GetHelix(7);
4544 Float_t r0 = helixes[i0].GetHelix(8);
4545 Float_t rc0 = TMath::Sqrt(xc0*xc0+yc0*yc0);
4546 Float_t fi0 = TMath::ATan2(yc0,xc0);
4548 for (Int_t i1=i0+1;i1<nentries;i1++){
4549 AliTPCseed * track1 = (AliTPCseed*)array->At(i1);
4550 if (!track1) continue;
4551 if (TMath::Abs(track1->GetC())<1/kMaxC) continue;
4552 Float_t xc1 = helixes[i1].GetHelix(6);
4553 Float_t yc1 = helixes[i1].GetHelix(7);
4554 Float_t r1 = helixes[i1].GetHelix(8);
4555 Float_t rc1 = TMath::Sqrt(xc1*xc1+yc1*yc1);
4556 Float_t fi1 = TMath::ATan2(yc1,xc1);
4558 Float_t dfi = fi0-fi1;
4561 if (dfi>1.5*TMath::Pi()) dfi-=TMath::Pi(); // take care about edge effect
4562 if (dfi<-1.5*TMath::Pi()) dfi+=TMath::Pi(); //
4563 Float_t dtheta = TMath::Abs(track0->GetTgl()-track1->GetTgl())<TMath::Abs(track0->GetTgl()+track1->GetTgl())? track0->GetTgl()-track1->GetTgl():track0->GetTgl()+track1->GetTgl();
4567 if (track0->GetBConstrain()&&track1->GetBConstrain()) continue; // not constrained
4568 if (track1->GetSigned1Pt()*track0->GetSigned1Pt()>0) continue; // not the same sign
4569 if ( TMath::Abs(track1->GetTgl()+track0->GetTgl())>kMaxdTheta) continue; //distance in the Theta
4570 if ( TMath::Abs(dfi)>kMaxdPhi) continue; //distance in phi
4571 if ( TMath::Sqrt(dfi*dfi+dtheta*dtheta)>kMaxdPhi) continue; //common angular offset
4573 Float_t pt0 = track0->GetSignedPt();
4574 Float_t pt1 = track1->GetSignedPt();
4575 if ((TMath::Abs(pt0+pt1)/(TMath::Abs(pt0)+TMath::Abs(pt1)))>kPtRatio) continue;
4576 if ((iter==1) && TMath::Abs(TMath::Abs(rc0+r0)-TMath::Abs(rc1+r1))>kMaxDeltaRMax) continue;
4577 if ((iter!=1) &&TMath::Abs(TMath::Abs(rc0-r0)-TMath::Abs(rc1-r1))>kMaxDeltaRMin) continue;
4578 if (TMath::Min(TMath::Abs(rc0-r0),TMath::Abs(rc1-r1))<kMinDCAR) continue;
4581 // Now find closest approach
4585 Int_t npoints = helixes[i0].GetRPHIintersections(helixes[i1], phase, radius,10);
4586 if (npoints==0) continue;
4587 helixes[i0].GetClosestPhases(helixes[i1], phase);
4591 Double_t hangles[3];
4592 helixes[i0].Evaluate(phase[0][0],xyz0);
4593 helixes[i1].Evaluate(phase[0][1],xyz1);
4595 helixes[i0].GetAngle(phase[0][0],helixes[i1],phase[0][1],hangles);
4596 Double_t deltah[2],deltabest;
4597 if (TMath::Abs(hangles[2])<kMinAngle) continue;
4601 helixes[i0].ParabolicDCA(helixes[i1],phase[0][0],phase[0][1],radius[0],deltah[0],2);
4603 helixes[i0].ParabolicDCA(helixes[i1],phase[1][0],phase[1][1],radius[1],deltah[1],2);
4604 if (deltah[1]<deltah[0]) ibest=1;
4606 deltabest = TMath::Sqrt(deltah[ibest]);
4607 helixes[i0].Evaluate(phase[ibest][0],xyz0);
4608 helixes[i1].Evaluate(phase[ibest][1],xyz1);
4609 helixes[i0].GetAngle(phase[ibest][0],helixes[i1],phase[ibest][1],hangles);
4610 Double_t radiusbest = TMath::Sqrt(radius[ibest]);
4612 if (deltabest>kMaxDist) continue;
4613 // if (mindcar+mindcaz<40 && (TMath::Abs(hangles[2])<kMinAngle ||deltabest>3)) continue;
4614 Bool_t sign =kFALSE;
4615 if (hangles[2]>kMinAngle) sign =kTRUE;
4618 // circular[i0] = kTRUE;
4619 // circular[i1] = kTRUE;
4620 if (track0->OneOverPt()<track1->OneOverPt()){
4621 track0->SetCircular(track0->GetCircular()+1);
4622 track1->SetCircular(track1->GetCircular()+2);
4625 track1->SetCircular(track1->GetCircular()+1);
4626 track0->SetCircular(track0->GetCircular()+2);
4629 if (AliTPCReconstructor::StreamLevel()>1){
4631 //debug stream to tune "fine" cuts
4632 Int_t lab0=track0->GetLabel();
4633 Int_t lab1=track1->GetLabel();
4634 TTreeSRedirector &cstream = *fDebugStreamer;
4635 cstream<<"Curling2"<<
4651 "npoints="<<npoints<<
4652 "hangles0="<<hangles[0]<<
4653 "hangles1="<<hangles[1]<<
4654 "hangles2="<<hangles[2]<<
4657 "radius="<<radiusbest<<
4658 "deltabest="<<deltabest<<
4659 "phase0="<<phase[ibest][0]<<
4660 "phase1="<<phase[ibest][1]<<
4668 if (AliTPCReconstructor::StreamLevel()>1) {
4669 AliInfo("Time for curling tracks removal");
4678 void AliTPCtrackerMI::FindKinks(TObjArray * array, AliESDEvent *esd)
4685 TObjArray *kinks= new TObjArray(10000);
4686 // TObjArray *v0s= new TObjArray(10000);
4687 Int_t nentries = array->GetEntriesFast();
4688 AliHelix *helixes = new AliHelix[nentries];
4689 Int_t *sign = new Int_t[nentries];
4690 Int_t *nclusters = new Int_t[nentries];
4691 Float_t *alpha = new Float_t[nentries];
4692 AliKink *kink = new AliKink();
4693 Int_t * usage = new Int_t[nentries];
4694 Float_t *zm = new Float_t[nentries];
4695 Float_t *z0 = new Float_t[nentries];
4696 Float_t *fim = new Float_t[nentries];
4697 Float_t *shared = new Float_t[nentries];
4698 Bool_t *circular = new Bool_t[nentries];
4699 Float_t *dca = new Float_t[nentries];
4700 //const AliESDVertex * primvertex = esd->GetVertex();
4702 // nentries = array->GetEntriesFast();
4707 for (Int_t i=0;i<nentries;i++){
4710 AliTPCseed* track = (AliTPCseed*)array->At(i);
4711 if (!track) continue;
4712 track->SetCircular(0);
4714 track->UpdatePoints();
4715 if (( track->GetPoints()[2]- track->GetPoints()[0])>5 && track->GetPoints()[3]>0.8){
4717 nclusters[i]=track->GetNumberOfClusters();
4718 alpha[i] = track->GetAlpha();
4719 new (&helixes[i]) AliHelix(*track);
4721 helixes[i].Evaluate(0,xyz);
4722 sign[i] = (track->GetC()>0) ? -1:1;
4725 if (track->GetProlongation(x,y,z)){
4727 fim[i] = alpha[i]+TMath::ATan2(y,x);
4730 zm[i] = track->GetZ();
4734 circular[i]= kFALSE;
4735 if (track->GetProlongation(0,y,z)) z0[i] = z;
4736 dca[i] = track->GetD(0,0);
4742 Int_t ncandidates =0;
4745 Double_t phase[2][2],radius[2];
4748 // Find circling track
4750 for (Int_t i0=0;i0<nentries;i0++){
4751 AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
4752 if (!track0) continue;
4753 if (track0->GetNumberOfClusters()<40) continue;
4754 if (TMath::Abs(1./track0->GetC())>200) continue;
4755 for (Int_t i1=i0+1;i1<nentries;i1++){
4756 AliTPCseed * track1 = (AliTPCseed*)array->At(i1);
4757 if (!track1) continue;
4758 if (track1->GetNumberOfClusters()<40) continue;
4759 if ( TMath::Abs(track1->GetTgl()+track0->GetTgl())>0.1) continue;
4760 if (track0->GetBConstrain()&&track1->GetBConstrain()) continue;
4761 if (TMath::Abs(1./track1->GetC())>200) continue;
4762 if (track1->GetSigned1Pt()*track0->GetSigned1Pt()>0) continue;
4763 if (track1->GetTgl()*track0->GetTgl()>0) continue;
4764 if (TMath::Max(TMath::Abs(1./track0->GetC()),TMath::Abs(1./track1->GetC()))>190) continue;
4765 if (track0->GetBConstrain()&&track1->OneOverPt()<track0->OneOverPt()) continue; //returning - lower momenta
4766 if (track1->GetBConstrain()&&track0->OneOverPt()<track1->OneOverPt()) continue; //returning - lower momenta
4768 Float_t mindcar = TMath::Min(TMath::Abs(dca[i0]),TMath::Abs(dca[i1]));
4769 if (mindcar<5) continue;
4770 Float_t mindcaz = TMath::Min(TMath::Abs(z0[i0]-GetZ()),TMath::Abs(z0[i1]-GetZ()));
4771 if (mindcaz<5) continue;
4772 if (mindcar+mindcaz<20) continue;
4775 Float_t xc0 = helixes[i0].GetHelix(6);
4776 Float_t yc0 = helixes[i0].GetHelix(7);
4777 Float_t r0 = helixes[i0].GetHelix(8);
4778 Float_t xc1 = helixes[i1].GetHelix(6);
4779 Float_t yc1 = helixes[i1].GetHelix(7);
4780 Float_t r1 = helixes[i1].GetHelix(8);
4782 Float_t rmean = (r0+r1)*0.5;
4783 Float_t delta =TMath::Sqrt((xc1-xc0)*(xc1-xc0)+(yc1-yc0)*(yc1-yc0));
4784 //if (delta>30) continue;
4785 if (delta>rmean*0.25) continue;
4786 if (TMath::Abs(r0-r1)/rmean>0.3) continue;
4788 Int_t npoints = helixes[i0].GetRPHIintersections(helixes[i1], phase, radius,10);
4789 if (npoints==0) continue;
4790 helixes[i0].GetClosestPhases(helixes[i1], phase);
4794 Double_t hangles[3];
4795 helixes[i0].Evaluate(phase[0][0],xyz0);
4796 helixes[i1].Evaluate(phase[0][1],xyz1);
4798 helixes[i0].GetAngle(phase[0][0],helixes[i1],phase[0][1],hangles);
4799 Double_t deltah[2],deltabest;
4800 if (hangles[2]<2.8) continue;
4803 helixes[i0].ParabolicDCA(helixes[i1],phase[0][0],phase[0][1],radius[0],deltah[0],2);
4805 helixes[i0].ParabolicDCA(helixes[i1],phase[1][0],phase[1][1],radius[1],deltah[1],2);
4806 if (deltah[1]<deltah[0]) ibest=1;
4808 deltabest = TMath::Sqrt(deltah[ibest]);
4809 helixes[i0].Evaluate(phase[ibest][0],xyz0);
4810 helixes[i1].Evaluate(phase[ibest][1],xyz1);
4811 helixes[i0].GetAngle(phase[ibest][0],helixes[i1],phase[ibest][1],hangles);
4812 Double_t radiusbest = TMath::Sqrt(radius[ibest]);
4814 if (deltabest>6) continue;
4815 if (mindcar+mindcaz<40 && (hangles[2]<3.12||deltabest>3)) continue;
4816 Bool_t sign =kFALSE;
4817 if (hangles[2]>3.06) sign =kTRUE;
4820 circular[i0] = kTRUE;
4821 circular[i1] = kTRUE;
4822 if (track0->OneOverPt()<track1->OneOverPt()){
4823 track0->SetCircular(track0->GetCircular()+1);
4824 track1->SetCircular(track1->GetCircular()+2);
4827 track1->SetCircular(track1->GetCircular()+1);
4828 track0->SetCircular(track0->GetCircular()+2);
4831 if (sign&&AliTPCReconstructor::StreamLevel()>1){
4833 Int_t lab0=track0->GetLabel();
4834 Int_t lab1=track1->GetLabel();
4835 TTreeSRedirector &cstream = *fDebugStreamer;
4836 cstream<<"Curling"<<
4843 "mindcar="<<mindcar<<
4844 "mindcaz="<<mindcaz<<
4847 "npoints="<<npoints<<
4848 "hangles0="<<hangles[0]<<
4849 "hangles2="<<hangles[2]<<
4854 "radius="<<radiusbest<<
4855 "deltabest="<<deltabest<<
4856 "phase0="<<phase[ibest][0]<<
4857 "phase1="<<phase[ibest][1]<<
4867 for (Int_t i =0;i<nentries;i++){
4868 if (sign[i]==0) continue;
4869 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
4876 Double_t cradius0 = 40*40;
4877 Double_t cradius1 = 270*270;
4880 Double_t cdist3=0.55;
4881 for (Int_t j =i+1;j<nentries;j++){
4883 if (sign[j]*sign[i]<1) continue;
4884 if ( (nclusters[i]+nclusters[j])>200) continue;
4885 if ( (nclusters[i]+nclusters[j])<80) continue;
4886 if ( TMath::Abs(zm[i]-zm[j])>60.) continue;
4887 if ( TMath::Abs(fim[i]-fim[j])>0.6 && TMath::Abs(fim[i]-fim[j])<5.7 ) continue;
4888 //AliTPCseed * track1 = (AliTPCseed*)array->At(j); Double_t phase[2][2],radius[2];
4889 Int_t npoints = helixes[i].GetRPHIintersections(helixes[j], phase, radius,20);
4890 if (npoints<1) continue;
4893 if (radius[0]<cradius0||radius[0]>cradius1) continue;
4896 if ( (radius[0]<cradius0||radius[0]>cradius1) && (radius[1]<cradius0||radius[1]>cradius1) ) continue;
4899 Double_t delta1=10000,delta2=10000;
4900 // cuts on the intersection radius
4901 helixes[i].LinearDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta1);
4902 if (radius[0]<20&&delta1<1) continue; //intersection at vertex
4903 if (radius[0]<10&&delta1<3) continue; //intersection at vertex
4905 helixes[i].LinearDCA(helixes[j],phase[1][0],phase[1][1],radius[1],delta2);
4906 if (radius[1]<20&&delta2<1) continue; //intersection at vertex
4907 if (radius[1]<10&&delta2<3) continue; //intersection at vertex
4910 Double_t distance1 = TMath::Min(delta1,delta2);
4911 if (distance1>cdist1) continue; // cut on DCA linear approximation
4913 npoints = helixes[i].GetRPHIintersections(helixes[j], phase, radius,20);
4914 helixes[i].ParabolicDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta1);
4915 if (radius[0]<20&&delta1<1) continue; //intersection at vertex
4916 if (radius[0]<10&&delta1<3) continue; //intersection at vertex
4919 helixes[i].ParabolicDCA(helixes[j],phase[1][0],phase[1][1],radius[1],delta2);
4920 if (radius[1]<20&&delta2<1) continue; //intersection at vertex
4921 if (radius[1]<10&&delta2<3) continue; //intersection at vertex
4923 distance1 = TMath::Min(delta1,delta2);
4926 rkink = TMath::Sqrt(radius[0]);
4929 rkink = TMath::Sqrt(radius[1]);
4931 if (distance1>cdist2) continue;
4934 AliTPCseed * track1 = (AliTPCseed*)array->At(j);
4937 Int_t row0 = GetRowNumber(rkink);
4938 if (row0<10) continue;
4939 if (row0>150) continue;
4942 Float_t dens00=-1,dens01=-1;
4943 Float_t dens10=-1,dens11=-1;
4945 Int_t found,foundable,shared;
4946 track0->GetClusterStatistic(0,row0-5, found, foundable,shared,kFALSE);
4947 if (foundable>5) dens00 = Float_t(found)/Float_t(foundable);
4948 track0->GetClusterStatistic(row0+5,155, found, foundable,shared,kFALSE);
4949 if (foundable>5) dens01 = Float_t(found)/Float_t(foundable);
4951 track1->GetClusterStatistic(0,row0-5, found, foundable,shared,kFALSE);
4952 if (foundable>10) dens10 = Float_t(found)/Float_t(foundable);
4953 track1->GetClusterStatistic(row0+5,155, found, foundable,shared,kFALSE);
4954 if (foundable>10) dens11 = Float_t(found)/Float_t(foundable);
4956 if (dens00<dens10 && dens01<dens11) continue;
4957 if (dens00>dens10 && dens01>dens11) continue;
4958 if (TMath::Max(dens00,dens10)<0.1) continue;
4959 if (TMath::Max(dens01,dens11)<0.3) continue;
4961 if (TMath::Min(dens00,dens10)>0.6) continue;
4962 if (TMath::Min(dens01,dens11)>0.6) continue;
4965 AliTPCseed * ktrack0, *ktrack1;
4974 if (TMath::Abs(ktrack0->GetC())>5) continue; // cut on the curvature for mother particle
4975 AliExternalTrackParam paramm(*ktrack0);
4976 AliExternalTrackParam paramd(*ktrack1);
4977 if (row0>60&&ktrack1->GetReference().GetX()>90.)new (¶md) AliExternalTrackParam(ktrack1->GetReference());
4980 kink->SetMother(paramm);
4981 kink->SetDaughter(paramd);
4984 Float_t x[3] = { kink->GetPosition()[0],kink->GetPosition()[1],kink->GetPosition()[2]};
4986 fParam->Transform0to1(x,index);
4987 fParam->Transform1to2(x,index);
4988 row0 = GetRowNumber(x[0]);
4990 if (kink->GetR()<100) continue;
4991 if (kink->GetR()>240) continue;
4992 if (kink->GetPosition()[2]/kink->GetR()>AliTPCReconstructor::GetCtgRange()) continue; //out of fiducial volume
4993 if (kink->GetDistance()>cdist3) continue;
4994 Float_t dird = kink->GetDaughterP()[0]*kink->GetPosition()[0]+kink->GetDaughterP()[1]*kink->GetPosition()[1]; // rough direction estimate
4995 if (dird<0) continue;
4997 Float_t dirm = kink->GetMotherP()[0]*kink->GetPosition()[0]+kink->GetMotherP()[1]*kink->GetPosition()[1]; // rough direction estimate
4998 if (dirm<0) continue;
4999 Float_t mpt = TMath::Sqrt(kink->GetMotherP()[0]*kink->GetMotherP()[0]+kink->GetMotherP()[1]*kink->GetMotherP()[1]);
5000 if (mpt<0.2) continue;
5003 //for high momenta momentum not defined well in first iteration
5004 Double_t qt = TMath::Sin(kink->GetAngle(2))*ktrack1->GetP();
5005 if (qt>0.35) continue;
5008 kink->SetLabel(CookLabel(ktrack0,0.4,0,row0),0);
5009 kink->SetLabel(CookLabel(ktrack1,0.4,row0,160),1);
5011 kink->SetTPCDensity(dens00,0,0);
5012 kink->SetTPCDensity(dens01,0,1);
5013 kink->SetTPCDensity(dens10,1,0);
5014 kink->SetTPCDensity(dens11,1,1);
5015 kink->SetIndex(i,0);
5016 kink->SetIndex(j,1);
5019 kink->SetTPCDensity(dens10,0,0);
5020 kink->SetTPCDensity(dens11,0,1);
5021 kink->SetTPCDensity(dens00,1,0);
5022 kink->SetTPCDensity(dens01,1,1);
5023 kink->SetIndex(j,0);
5024 kink->SetIndex(i,1);
5027 if (mpt<1||kink->GetAngle(2)>0.1){
5028 // angle and densities not defined yet
5029 if (kink->GetTPCDensityFactor()<0.8) continue;
5030 if ((2-kink->GetTPCDensityFactor())*kink->GetDistance() >0.25) continue;
5031 if (kink->GetAngle(2)*ktrack0->GetP()<0.003) continue; //too small angle
5032 if (kink->GetAngle(2)>0.2&&kink->GetTPCDensityFactor()<1.15) continue;
5033 if (kink->GetAngle(2)>0.2&&kink->GetTPCDensity(0,1)>0.05) continue;
5035 Float_t criticalangle = track0->GetSigmaSnp2()+track0->GetSigmaTgl2();
5036 criticalangle+= track1->GetSigmaSnp2()+track1->GetSigmaTgl2();
5037 criticalangle= 3*TMath::Sqrt(criticalangle);
5038 if (criticalangle>0.02) criticalangle=0.02;
5039 if (kink->GetAngle(2)<criticalangle) continue;
5042 Int_t drow = Int_t(2.+0.5/(0.05+kink->GetAngle(2))); // overlap region defined
5043 Float_t shapesum =0;
5045 for ( Int_t row = row0-drow; row<row0+drow;row++){
5046 if (row<0) continue;
5047 if (row>155) continue;
5048 if (ktrack0->GetClusterPointer(row)){
5049 AliTPCTrackerPoint *point =ktrack0->GetTrackPoint(row);
5050 shapesum+=point->GetSigmaY()+point->GetSigmaZ();
5053 if (ktrack1->GetClusterPointer(row)){
5054 AliTPCTrackerPoint *point =ktrack1->GetTrackPoint(row);
5055 shapesum+=point->GetSigmaY()+point->GetSigmaZ();
5060 kink->SetShapeFactor(-1.);
5063 kink->SetShapeFactor(shapesum/sum);
5065 // esd->AddKink(kink);
5066 kinks->AddLast(kink);
5072 // sort the kinks according quality - and refit them towards vertex
5074 Int_t nkinks = kinks->GetEntriesFast();
5075 Float_t *quality = new Float_t[nkinks];
5076 Int_t *indexes = new Int_t[nkinks];
5077 AliTPCseed *mothers = new AliTPCseed[nkinks];
5078 AliTPCseed *daughters = new AliTPCseed[nkinks];
5081 for (Int_t i=0;i<nkinks;i++){
5083 AliKink *kink = (AliKink*)kinks->At(i);
5085 // refit kinks towards vertex
5087 Int_t index0 = kink->GetIndex(0);
5088 Int_t index1 = kink->GetIndex(1);
5089 AliTPCseed * ktrack0 = (AliTPCseed*)array->At(index0);
5090 AliTPCseed * ktrack1 = (AliTPCseed*)array->At(index1);
5092 Int_t sumn=ktrack0->GetNumberOfClusters()+ktrack1->GetNumberOfClusters();
5094 // Refit Kink under if too small angle
5096 if (kink->GetAngle(2)<0.05){
5097 kink->SetTPCRow0(GetRowNumber(kink->GetR()));
5098 Int_t row0 = kink->GetTPCRow0();
5099 Int_t drow = Int_t(2.+0.5/(0.05+kink->GetAngle(2)));
5102 Int_t last = row0-drow;
5103 if (last<40) last=40;
5104 if (last<ktrack0->GetFirstPoint()+25) last = ktrack0->GetFirstPoint()+25;
5105 AliTPCseed* seed0 = ReSeed(ktrack0,last,kFALSE);
5108 Int_t first = row0+drow;
5109 if (first>130) first=130;
5110 if (first>ktrack1->GetLastPoint()-25) first = TMath::Max(ktrack1->GetLastPoint()-25,30);
5111 AliTPCseed* seed1 = ReSeed(ktrack1,first,kTRUE);
5113 if (seed0 && seed1){
5114 kink->SetStatus(1,8);
5115 if (RefitKink(*seed0,*seed1,*kink)) kink->SetStatus(1,9);
5116 row0 = GetRowNumber(kink->GetR());
5117 sumn = seed0->GetNumberOfClusters()+seed1->GetNumberOfClusters();
5118 mothers[i] = *seed0;
5119 daughters[i] = *seed1;
5122 delete kinks->RemoveAt(i);
5123 if (seed0) delete seed0;
5124 if (seed1) delete seed1;
5127 if (kink->GetDistance()>0.5 || kink->GetR()<110 || kink->GetR()>240) {
5128 delete kinks->RemoveAt(i);
5129 if (seed0) delete seed0;
5130 if (seed1) delete seed1;
5138 if (kink) quality[i] = 160*((0.1+kink->GetDistance())*(2.-kink->GetTPCDensityFactor()))/(sumn+40.); //the longest -clossest will win
5140 TMath::Sort(nkinks,quality,indexes,kFALSE);
5142 //remove double find kinks
5144 for (Int_t ikink0=1;ikink0<nkinks;ikink0++){
5145 AliKink * kink0 = (AliKink*) kinks->At(indexes[ikink0]);
5146 if (!kink0) continue;
5148 for (Int_t ikink1=0;ikink1<ikink0;ikink1++){
5149 if (!kink0) continue;
5150 AliKink * kink1 = (AliKink*) kinks->At(indexes[ikink1]);
5151 if (!kink1) continue;
5152 // if not close kink continue
5153 if (TMath::Abs(kink1->GetPosition()[2]-kink0->GetPosition()[2])>10) continue;
5154 if (TMath::Abs(kink1->GetPosition()[1]-kink0->GetPosition()[1])>10) continue;
5155 if (TMath::Abs(kink1->GetPosition()[0]-kink0->GetPosition()[0])>10) continue;
5157 AliTPCseed &mother0 = mothers[indexes[ikink0]];
5158 AliTPCseed &daughter0 = daughters[indexes[ikink0]];
5159 AliTPCseed &mother1 = mothers[indexes[ikink1]];
5160 AliTPCseed &daughter1 = daughters[indexes[ikink1]];
5161 Int_t row0 = (kink0->GetTPCRow0()+kink1->GetTPCRow0())/2;
5170 for (Int_t i=0;i<row0;i++){
5171 if (mother0.GetClusterIndex(i)>0 && mother1.GetClusterIndex(i)>0){
5174 if (mother0.GetClusterIndex(i)==mother1.GetClusterIndex(i)){
5181 for (Int_t i=row0;i<158;i++){
5182 if (daughter0.GetClusterIndex(i)>0 && daughter0.GetClusterIndex(i)>0){
5185 if (mother0.GetClusterIndex(i)==mother1.GetClusterIndex(i)){
5191 Float_t ratio = Float_t(same+1)/Float_t(both+1);
5192 Float_t ratiom = Float_t(samem+1)/Float_t(bothm+1);
5193 Float_t ratiod = Float_t(samed+1)/Float_t(bothd+1);
5194 if (ratio>0.3 && ratiom>0.5 &&ratiod>0.5) {
5195 Int_t sum0 = mother0.GetNumberOfClusters()+daughter0.GetNumberOfClusters();
5196 Int_t sum1 = mother1.GetNumberOfClusters()+daughter1.GetNumberOfClusters();
5198 shared[kink0->GetIndex(0)]= kTRUE;
5199 shared[kink0->GetIndex(1)]= kTRUE;
5200 delete kinks->RemoveAt(indexes[ikink0]);
5203 shared[kink1->GetIndex(0)]= kTRUE;
5204 shared[kink1->GetIndex(1)]= kTRUE;
5205 delete kinks->RemoveAt(indexes[ikink1]);
5212 for (Int_t i=0;i<nkinks;i++){
5213 AliKink * kink = (AliKink*) kinks->At(indexes[i]);
5214 if (!kink) continue;
5215 kink->SetTPCRow0(GetRowNumber(kink->GetR()));
5216 Int_t index0 = kink->GetIndex(0);
5217 Int_t index1 = kink->GetIndex(1);
5218 if (circular[index0]||circular[index1]&&kink->GetDistance()>0.2) continue;
5219 kink->SetMultiple(usage[index0],0);
5220 kink->SetMultiple(usage[index1],1);
5221 if (kink->GetMultiple()[0]+kink->GetMultiple()[1]>2) continue;
5222 if (kink->GetMultiple()[0]+kink->GetMultiple()[1]>0 && quality[indexes[i]]>0.2) continue;
5223 if (kink->GetMultiple()[0]+kink->GetMultiple()[1]>0 && kink->GetDistance()>0.2) continue;
5224 if (circular[index0]||circular[index1]&&kink->GetDistance()>0.1) continue;
5226 AliTPCseed * ktrack0 = (AliTPCseed*)array->At(index0);
5227 AliTPCseed * ktrack1 = (AliTPCseed*)array->At(index1);
5228 if (!ktrack0 || !ktrack1) continue;
5229 Int_t index = esd->AddKink(kink);
5232 if ( ktrack0->GetKinkIndex(0)==0 && ktrack1->GetKinkIndex(0)==0) { //best kink
5233 if (mothers[indexes[i]].GetNumberOfClusters()>20 && daughters[indexes[i]].GetNumberOfClusters()>20 && (mothers[indexes[i]].GetNumberOfClusters()+daughters[indexes[i]].GetNumberOfClusters())>100){
5234 *ktrack0 = mothers[indexes[i]];
5235 *ktrack1 = daughters[indexes[i]];
5239 ktrack0->SetKinkIndex(usage[index0],-(index+1));
5240 ktrack1->SetKinkIndex(usage[index1], (index+1));
5245 // Remove tracks corresponding to shared kink's
5247 for (Int_t i=0;i<nentries;i++){
5248 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5249 if (!track0) continue;
5250 if (track0->GetKinkIndex(0)!=0) continue;
5251 if (shared[i]) delete array->RemoveAt(i);
5256 RemoveUsed2(array,0.5,0.4,30);
5258 for (Int_t i=0;i<nentries;i++){
5259 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5260 if (!track0) continue;
5261 track0->CookdEdx(0.02,0.6);
5265 for (Int_t i=0;i<nentries;i++){
5266 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5267 if (!track0) continue;
5268 if (track0->Pt()<1.4) continue;
5269 //remove double high momenta tracks - overlapped with kink candidates
5272 for (Int_t icl=track0->GetFirstPoint();icl<track0->GetLastPoint(); icl++){
5273 if (track0->GetClusterPointer(icl)!=0){
5275 if (track0->GetClusterPointer(icl)->IsUsed(10)) shared++;
5278 if (Float_t(shared+1)/Float_t(all+1)>0.5) {
5279 delete array->RemoveAt(i);
5283 if (track0->GetKinkIndex(0)!=0) continue;
5284 if (track0->GetNumberOfClusters()<80) continue;
5286 AliTPCseed *pmother = new AliTPCseed();
5287 AliTPCseed *pdaughter = new AliTPCseed();
5288 AliKink *pkink = new AliKink;
5290 AliTPCseed & mother = *pmother;
5291 AliTPCseed & daughter = *pdaughter;
5292 AliKink & kink = *pkink;
5293 if (CheckKinkPoint(track0,mother,daughter, kink)){
5294 if (mother.GetNumberOfClusters()<30||daughter.GetNumberOfClusters()<20) {
5298 continue; //too short tracks
5300 if (mother.Pt()<1.4) {
5306 Int_t row0= kink.GetTPCRow0();
5307 if (kink.GetDistance()>0.5 || kink.GetR()<110. || kink.GetR()>240.) {
5314 Int_t index = esd->AddKink(&kink);
5315 mother.SetKinkIndex(0,-(index+1));
5316 daughter.SetKinkIndex(0,index+1);
5317 if (mother.GetNumberOfClusters()>50) {
5318 delete array->RemoveAt(i);
5319 array->AddAt(new AliTPCseed(mother),i);
5322 array->AddLast(new AliTPCseed(mother));
5324 array->AddLast(new AliTPCseed(daughter));
5325 for (Int_t icl=0;icl<row0;icl++) {
5326 if (mother.GetClusterPointer(icl)) mother.GetClusterPointer(icl)->Use(20);
5329 for (Int_t icl=row0;icl<158;icl++) {
5330 if (daughter.GetClusterPointer(icl)) daughter.GetClusterPointer(icl)->Use(20);
5339 delete [] daughters;
5361 printf("Ncandidates=\t%d\t%d\t%d\t%d\n",esd->GetNumberOfKinks(),ncandidates,ntracks,nall);
5365 void AliTPCtrackerMI::FindV0s(TObjArray * array, AliESDEvent *esd)
5371 TObjArray *tpcv0s = new TObjArray(100000);
5372 Int_t nentries = array->GetEntriesFast();
5373 AliHelix *helixes = new AliHelix[nentries];
5374 Int_t *sign = new Int_t[nentries];
5375 Float_t *alpha = new Float_t[nentries];
5376 Float_t *z0 = new Float_t[nentries];
5377 Float_t *dca = new Float_t[nentries];
5378 Float_t *sdcar = new Float_t[nentries];
5379 Float_t *cdcar = new Float_t[nentries];
5380 Float_t *pulldcar = new Float_t[nentries];
5381 Float_t *pulldcaz = new Float_t[nentries];
5382 Float_t *pulldca = new Float_t[nentries];
5383 Bool_t *isPrim = new Bool_t[nentries];
5384 const AliESDVertex * primvertex = esd->GetVertex();
5385 Double_t zvertex = primvertex->GetZv();
5387 // nentries = array->GetEntriesFast();
5389 for (Int_t i=0;i<nentries;i++){
5392 AliTPCseed* track = (AliTPCseed*)array->At(i);
5393 if (!track) continue;
5394 track->GetV0Indexes()[0] = 0; //rest v0 indexes
5395 track->GetV0Indexes()[1] = 0; //rest v0 indexes
5396 track->GetV0Indexes()[2] = 0; //rest v0 indexes
5398 alpha[i] = track->GetAlpha();
5399 new (&helixes[i]) AliHelix(*track);
5401 helixes[i].Evaluate(0,xyz);
5402 sign[i] = (track->GetC()>0) ? -1:1;
5406 if (track->GetProlongation(0,y,z)) z0[i] = z;
5407 dca[i] = track->GetD(0,0);
5409 // dca error parrameterezation + pulls
5411 sdcar[i] = TMath::Sqrt(0.150*0.150+(100*track->GetC())*(100*track->GetC()));
5412 if (TMath::Abs(track->GetTgl())>1) sdcar[i]*=2.5;
5413 cdcar[i] = TMath::Exp((TMath::Abs(track->GetC())-0.0106)*525.3);
5414 pulldcar[i] = (dca[i]-cdcar[i])/sdcar[i];
5415 pulldcaz[i] = (z0[i]-zvertex)/sdcar[i];
5416 pulldca[i] = TMath::Sqrt(pulldcar[i]*pulldcar[i]+pulldcaz[i]*pulldcaz[i]);
5417 if (track->TPCrPID(1)+track->TPCrPID(2)+track->TPCrPID(3)>0.5) {
5418 if (pulldca[i]<3.) isPrim[i]=kTRUE; //pion, muon and Kaon 3 sigma cut
5420 if (track->TPCrPID(4)>0.5) {
5421 if (pulldca[i]<0.5) isPrim[i]=kTRUE; //proton 0.5 sigma cut
5423 if (track->TPCrPID(0)>0.4) {
5424 isPrim[i]=kFALSE; //electron no sigma cut
5431 Int_t ncandidates =0;
5434 Double_t phase[2][2],radius[2];
5440 Float_t fprimvertex[3]={GetX(),GetY(),GetZ()};
5442 Double_t cradius0 = 10*10;
5443 Double_t cradius1 = 200*200;
5446 Double_t cpointAngle = 0.95;
5448 Double_t delta[2]={10000,10000};
5449 for (Int_t i =0;i<nentries;i++){
5450 if (sign[i]==0) continue;
5451 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5452 if (!track0) continue;
5453 if (AliTPCReconstructor::StreamLevel()>1){
5454 TTreeSRedirector &cstream = *fDebugStreamer;
5459 "zvertex="<<zvertex<<
5460 "sdcar0="<<sdcar[i]<<
5461 "cdcar0="<<cdcar[i]<<
5462 "pulldcar0="<<pulldcar[i]<<
5463 "pulldcaz0="<<pulldcaz[i]<<
5464 "pulldca0="<<pulldca[i]<<
5465 "isPrim="<<isPrim[i]<<
5469 if (track0->GetSigned1Pt()<0) continue;
5470 if (track0->GetKinkIndex(0)>0||isPrim[i]) continue; //daughter kink
5472 if (TMath::Abs(helixes[i].GetHelix(4))<0.000000001) continue;
5477 for (Int_t j =0;j<nentries;j++){
5478 AliTPCseed * track1 = (AliTPCseed*)array->At(j);
5479 if (!track1) continue;
5480 if (track1->GetKinkIndex(0)>0 || isPrim[j]) continue; //daughter kink
5481 if (sign[j]*sign[i]>0) continue;
5482 if (TMath::Abs(helixes[j].GetHelix(4))<0.000001) continue;
5483 if (track0->GetCircular()+track1->GetCircular()>1) continue; //circling -returning track
5486 // DCA to prim vertex cut
5492 Int_t npoints = helixes[i].GetRPHIintersections(helixes[j], phase, radius,cdist2);
5493 if (npoints<1) continue;
5497 if (radius[0]<cradius0||radius[0]>cradius1) continue;
5498 helixes[i].LinearDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta[0]);
5499 if (delta[0]>cdist1) continue;
5502 if (TMath::Max(radius[0],radius[1])<cradius0|| TMath::Min(radius[0],radius[1])>cradius1) continue;
5503 helixes[i].LinearDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta[0]);
5504 helixes[i].LinearDCA(helixes[j],phase[1][0],phase[1][1],radius[1],delta[1]);
5505 if (delta[1]<delta[0]) iclosest=1;
5506 if (delta[iclosest]>cdist1) continue;
5508 helixes[i].ParabolicDCA(helixes[j],phase[iclosest][0],phase[iclosest][1],radius[iclosest],delta[iclosest]);
5509 if (radius[iclosest]<cradius0 || radius[iclosest]>cradius1 || delta[iclosest]>cdist1) continue;
5511 Double_t pointAngle = helixes[i].GetPointAngle(helixes[j],phase[iclosest],fprimvertex);
5512 if (pointAngle<cpointAngle) continue;
5514 Bool_t isGamma = kFALSE;
5515 vertex.SetParamP(*track0); //track0 - plus
5516 vertex.SetParamN(*track1); //track1 - minus
5517 vertex.Update(fprimvertex);
5518 if (track0->TPCrPID(0)>0.3&&track1->TPCrPID(0)>0.3&&vertex.GetAnglep()[2]<0.15) isGamma=kTRUE; // gamma conversion candidate
5519 Double_t pointAngle2 = vertex.GetV0CosineOfPointingAngle();
5521 if (vertex.GetV0CosineOfPointingAngle()<cpointAngle && (!isGamma)) continue;// point angle cut
5522 //Bo: if (vertex.GetDist2()>2&&(!isGamma)) continue; // point angle cut
5523 if (vertex.GetDcaV0Daughters()>2&&(!isGamma)) continue;//Bo: // point angle cut
5524 Float_t sigmae = 0.15*0.15;
5525 if (vertex.GetRr()<80)
5526 sigmae += (sdcar[i]*sdcar[i]+sdcar[j]*sdcar[j])*(1.-vertex.GetRr()/80.)*(1.-vertex.GetRr()/80.);
5527 sigmae+= TMath::Sqrt(sigmae);
5528 //Bo: if (vertex.GetDist2()/sigmae>3.&&(!isGamma)) continue;
5529 if (vertex.GetDcaV0Daughters()/sigmae>3.&&(!isGamma)) continue;
5530 Float_t densb0=0,densb1=0,densa0=0,densa1=0;
5531 Int_t row0 = GetRowNumber(vertex.GetRr());
5533 //Bo: if (vertex.GetDist2()>0.2) continue;
5534 if (vertex.GetDcaV0Daughters()>0.2) continue;
5535 densb0 = track0->Density2(0,row0-5);
5536 densb1 = track1->Density2(0,row0-5);
5537 if (densb0>0.3|| densb1>0.3) continue; //clusters before vertex
5538 densa0 = track0->Density2(row0+5,row0+40);
5539 densa1 = track1->Density2(row0+5,row0+40);
5540 if ((densa0<0.4|| densa1<0.4)&&(!isGamma)) continue; //missing clusters after vertex
5543 densa0 = track0->Density2(0,40); //cluster density
5544 densa1 = track1->Density2(0,40); //cluster density
5545 if ((vertex.GetRr()<80&&densa0+densa1<1.)&&(!isGamma)) continue;
5547 //Bo: vertex.SetLab(0,track0->GetLabel());
5548 //Bo: vertex.SetLab(1,track1->GetLabel());
5549 vertex.SetChi2After((densa0+densa1)*0.5);
5550 vertex.SetChi2Before((densb0+densb1)*0.5);
5551 vertex.SetIndex(0,i);
5552 vertex.SetIndex(1,j);
5553 //Bo: vertex.SetStatus(1); // TPC v0 candidate
5554 vertex.SetOnFlyStatus(2);//Bo: // TPC v0 candidate
5555 //Bo: vertex.SetRp(track0->TPCrPIDs());
5556 //Bo: vertex.SetRm(track1->TPCrPIDs());
5557 tpcv0s->AddLast(new AliESDv0(vertex));
5560 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
5561 Double_t radiusm= (delta[0]<delta[1])? TMath::Sqrt(radius[0]):TMath::Sqrt(radius[1]);
5562 Double_t deltam= (delta[0]<delta[1])? TMath::Sqrt(delta[0]):TMath::Sqrt(delta[1]);
5563 if (AliTPCReconstructor::StreamLevel()>1) {
5564 Int_t lab0=track0->GetLabel();
5565 Int_t lab1=track1->GetLabel();
5566 Char_t c0=track0->GetCircular();
5567 Char_t c1=track1->GetCircular();
5568 TTreeSRedirector &cstream = *fDebugStreamer;
5571 "vertex.="<<&vertex<<
5574 "Helix0.="<<&helixes[i]<<
5577 "Helix1.="<<&helixes[j]<<
5578 "pointAngle="<<pointAngle<<
5579 "pointAngle2="<<pointAngle2<<
5584 "zvertex="<<zvertex<<
5587 "npoints="<<npoints<<
5588 "radius0="<<radius[0]<<
5589 "delta0="<<delta[0]<<
5590 "radius1="<<radius[1]<<
5591 "delta1="<<delta[1]<<
5592 "radiusm="<<radiusm<<
5594 "sdcar0="<<sdcar[i]<<
5595 "sdcar1="<<sdcar[j]<<
5596 "cdcar0="<<cdcar[i]<<
5597 "cdcar1="<<cdcar[j]<<
5598 "pulldcar0="<<pulldcar[i]<<
5599 "pulldcar1="<<pulldcar[j]<<
5600 "pulldcaz0="<<pulldcaz[i]<<
5601 "pulldcaz1="<<pulldcaz[j]<<
5602 "pulldca0="<<pulldca[i]<<
5603 "pulldca1="<<pulldca[j]<<
5613 Float_t *quality = new Float_t[ncandidates];
5614 Int_t *indexes = new Int_t[ncandidates];
5616 for (Int_t i=0;i<ncandidates;i++){
5618 AliESDv0 *v0 = (AliESDv0*)tpcv0s->At(i);
5619 quality[i] = 1./(1.00001-v0->GetV0CosineOfPointingAngle()); //base point angle
5620 // quality[i] /= (0.5+v0->GetDist2());
5621 // quality[i] *= v0->GetChi2After(); //density factor
5623 Int_t index0 = v0->GetIndex(0);
5624 Int_t index1 = v0->GetIndex(1);
5625 //Bo: Double_t minpulldca = TMath::Min(2.+pulldca[v0->GetIndex(0)],(2.+pulldca[v0->GetIndex(1)]) ); //pull
5626 Double_t minpulldca = TMath::Min(2.+pulldca[index0],(2.+pulldca[index1]) );//Bo:
5630 AliTPCseed * track0 = (AliTPCseed*)array->At(index0);
5631 AliTPCseed * track1 = (AliTPCseed*)array->At(index1);
5632 if (track0->TPCrPID(0)>0.3&&track1->TPCrPID(0)>0.3&&v0->GetAnglep()[2]<0.15) quality[i]+=1000000; // gamma conversion candidate
5633 if (track0->TPCrPID(4)>0.9||track1->TPCrPID(4)>0.9&&minpulldca>4) quality[i]*=10; // lambda candidate candidate
5636 TMath::Sort(ncandidates,quality,indexes,kTRUE);
5639 for (Int_t i=0;i<ncandidates;i++){
5640 AliESDv0 * v0 = (AliESDv0*)tpcv0s->At(indexes[i]);
5642 Int_t index0 = v0->GetIndex(0);
5643 Int_t index1 = v0->GetIndex(1);
5644 AliTPCseed * track0 = (AliTPCseed*)array->At(index0);
5645 AliTPCseed * track1 = (AliTPCseed*)array->At(index1);
5646 if (!track0||!track1) {
5650 Bool_t accept =kTRUE; //default accept
5651 Int_t *v0indexes0 = track0->GetV0Indexes();
5652 Int_t *v0indexes1 = track1->GetV0Indexes();
5654 Int_t order0 = (v0indexes0[0]!=0) ? 1:0;
5655 Int_t order1 = (v0indexes1[0]!=0) ? 1:0;
5656 if (v0indexes0[1]!=0) order0 =2;
5657 if (v0indexes1[1]!=0) order1 =2;
5659 if (v0indexes0[2]!=0) {order0=3; accept=kFALSE;}
5660 if (v0indexes0[2]!=0) {order1=3; accept=kFALSE;}
5662 AliESDv0 * v02 = v0;
5664 //Bo: v0->SetOrder(0,order0);
5665 //Bo: v0->SetOrder(1,order1);
5666 //Bo: v0->SetOrder(1,order0+order1);
5667 v0->SetOnFlyStatus(kTRUE);
5668 Int_t index = esd->AddV0(v0);
5669 v02 = esd->GetV0(index);
5670 v0indexes0[order0]=index;
5671 v0indexes1[order1]=index;
5675 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
5676 if (AliTPCReconstructor::StreamLevel()>1) {
5677 Int_t lab0=track0->GetLabel();
5678 Int_t lab1=track1->GetLabel();
5679 TTreeSRedirector &cstream = *fDebugStreamer;
5688 "dca0="<<dca[index0]<<
5689 "dca1="<<dca[index1]<<
5693 "quality="<<quality[i]<<
5694 "pulldca0="<<pulldca[index0]<<
5695 "pulldca1="<<pulldca[index1]<<
5719 printf("TPC V0 finder : naccepted\t%d\tncandidates\t%d\tntracks\t%d\tnall\t%d\n",naccepted,ncandidates,ntracks,nall);
5723 Int_t AliTPCtrackerMI::RefitKink(AliTPCseed &mother, AliTPCseed &daughter, AliESDkink &knk)
5726 // refit kink towards to the vertex
5729 AliKink &kink=(AliKink &)knk;
5731 Int_t row0 = GetRowNumber(kink.GetR());
5732 FollowProlongation(mother,0);
5733 mother.Reset(kFALSE);
5735 FollowProlongation(daughter,row0);
5736 daughter.Reset(kFALSE);
5737 FollowBackProlongation(daughter,158);
5738 daughter.Reset(kFALSE);
5739 Int_t first = TMath::Max(row0-20,30);
5740 Int_t last = TMath::Min(row0+20,140);
5742 const Int_t kNdiv =5;
5743 AliTPCseed param0[kNdiv]; // parameters along the track
5744 AliTPCseed param1[kNdiv]; // parameters along the track
5745 AliKink kinks[kNdiv]; // corresponding kink parameters
5748 for (Int_t irow=0; irow<kNdiv;irow++){
5749 rows[irow] = first +((last-first)*irow)/(kNdiv-1);
5751 // store parameters along the track
5753 for (Int_t irow=0;irow<kNdiv;irow++){
5754 FollowBackProlongation(mother, rows[irow]);
5755 FollowProlongation(daughter,rows[kNdiv-1-irow]);
5756 param0[irow] = mother;
5757 param1[kNdiv-1-irow] = daughter;
5761 for (Int_t irow=0; irow<kNdiv-1;irow++){
5762 if (param0[irow].GetNumberOfClusters()<kNdiv||param1[irow].GetNumberOfClusters()<kNdiv) continue;
5763 kinks[irow].SetMother(param0[irow]);
5764 kinks[irow].SetDaughter(param1[irow]);
5765 kinks[irow].Update();
5768 // choose kink with best "quality"
5770 Double_t mindist = 10000;
5771 for (Int_t irow=0;irow<kNdiv;irow++){
5772 if (param0[irow].GetNumberOfClusters()<20||param1[irow].GetNumberOfClusters()<20) continue;
5773 if (TMath::Abs(kinks[irow].GetR())>240.) continue;
5774 if (TMath::Abs(kinks[irow].GetR())<100.) continue;
5776 Float_t normdist = TMath::Abs(param0[irow].GetX()-kinks[irow].GetR())*(0.1+kink.GetDistance());
5777 normdist/= (param0[irow].GetNumberOfClusters()+param1[irow].GetNumberOfClusters()+40.);
5778 if (normdist < mindist){
5784 if (index==-1) return 0;
5787 param0[index].Reset(kTRUE);
5788 FollowProlongation(param0[index],0);
5790 mother = param0[index];
5791 daughter = param1[index]; // daughter in vertex
5793 kink.SetMother(mother);
5794 kink.SetDaughter(daughter);
5796 kink.SetTPCRow0(GetRowNumber(kink.GetR()));
5797 kink.SetTPCncls(param0[index].GetNumberOfClusters(),0);
5798 kink.SetTPCncls(param1[index].GetNumberOfClusters(),1);
5799 kink.SetLabel(CookLabel(&mother,0.4, 0,kink.GetTPCRow0()),0);
5800 kink.SetLabel(CookLabel(&daughter,0.4, kink.GetTPCRow0(),160),1);
5801 mother.SetLabel(kink.GetLabel(0));
5802 daughter.SetLabel(kink.GetLabel(1));
5808 void AliTPCtrackerMI::UpdateKinkQualityM(AliTPCseed * seed){
5810 // update Kink quality information for mother after back propagation
5812 if (seed->GetKinkIndex(0)>=0) return;
5813 for (Int_t ikink=0;ikink<3;ikink++){
5814 Int_t index = seed->GetKinkIndex(ikink);
5815 if (index>=0) break;
5816 index = TMath::Abs(index)-1;
5817 AliESDkink * kink = fEvent->GetKink(index);
5818 kink->SetTPCDensity(-1,0,0);
5819 kink->SetTPCDensity(1,0,1);
5821 Int_t row0 = kink->GetTPCRow0() - 2 - Int_t( 0.5/ (0.05+kink->GetAngle(2)));
5822 if (row0<15) row0=15;
5824 Int_t row1 = kink->GetTPCRow0() + 2 + Int_t( 0.5/ (0.05+kink->GetAngle(2)));
5825 if (row1>145) row1=145;
5827 Int_t found,foundable,shared;
5828 seed->GetClusterStatistic(0,row0, found, foundable,shared,kFALSE);
5829 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),0,0);
5830 seed->GetClusterStatistic(row1,155, found, foundable,shared,kFALSE);
5831 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),0,1);
5836 void AliTPCtrackerMI::UpdateKinkQualityD(AliTPCseed * seed){
5838 // update Kink quality information for daughter after refit
5840 if (seed->GetKinkIndex(0)<=0) return;
5841 for (Int_t ikink=0;ikink<3;ikink++){
5842 Int_t index = seed->GetKinkIndex(ikink);
5843 if (index<=0) break;
5844 index = TMath::Abs(index)-1;
5845 AliESDkink * kink = fEvent->GetKink(index);
5846 kink->SetTPCDensity(-1,1,0);
5847 kink->SetTPCDensity(-1,1,1);
5849 Int_t row0 = kink->GetTPCRow0() -2 - Int_t( 0.5/ (0.05+kink->GetAngle(2)));
5850 if (row0<15) row0=15;
5852 Int_t row1 = kink->GetTPCRow0() +2 + Int_t( 0.5/ (0.05+kink->GetAngle(2)));
5853 if (row1>145) row1=145;
5855 Int_t found,foundable,shared;
5856 seed->GetClusterStatistic(0,row0, found, foundable,shared,kFALSE);
5857 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),1,0);
5858 seed->GetClusterStatistic(row1,155, found, foundable,shared,kFALSE);
5859 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),1,1);
5865 Int_t AliTPCtrackerMI::CheckKinkPoint(AliTPCseed*seed,AliTPCseed &mother, AliTPCseed &daughter, AliESDkink &knk)
5868 // check kink point for given track
5869 // if return value=0 kink point not found
5870 // otherwise seed0 correspond to mother particle
5871 // seed1 correspond to daughter particle
5872 // kink parameter of kink point
5873 AliKink &kink=(AliKink &)knk;
5875 Int_t middlerow = (seed->GetFirstPoint()+seed->GetLastPoint())/2;
5876 Int_t first = seed->GetFirstPoint();
5877 Int_t last = seed->GetLastPoint();
5878 if (last-first<20) return 0; // shortest length - 2*30 = 60 pad-rows
5881 AliTPCseed *seed1 = ReSeed(seed,middlerow+20, kTRUE); //middle of chamber
5882 if (!seed1) return 0;
5883 FollowProlongation(*seed1,seed->GetLastPoint()-20);
5884 seed1->Reset(kTRUE);
5885 FollowProlongation(*seed1,158);
5886 seed1->Reset(kTRUE);
5887 last = seed1->GetLastPoint();
5889 AliTPCseed *seed0 = new AliTPCseed(*seed);
5890 seed0->Reset(kFALSE);
5893 AliTPCseed param0[20]; // parameters along the track
5894 AliTPCseed param1[20]; // parameters along the track
5895 AliKink kinks[20]; // corresponding kink parameters
5897 for (Int_t irow=0; irow<20;irow++){
5898 rows[irow] = first +((last-first)*irow)/19;
5900 // store parameters along the track
5902 for (Int_t irow=0;irow<20;irow++){
5903 FollowBackProlongation(*seed0, rows[irow]);
5904 FollowProlongation(*seed1,rows[19-irow]);
5905 param0[irow] = *seed0;
5906 param1[19-irow] = *seed1;
5910 for (Int_t irow=0; irow<19;irow++){
5911 kinks[irow].SetMother(param0[irow]);
5912 kinks[irow].SetDaughter(param1[irow]);
5913 kinks[irow].Update();
5916 // choose kink with biggest change of angle
5918 Double_t maxchange= 0;
5919 for (Int_t irow=1;irow<19;irow++){
5920 if (TMath::Abs(kinks[irow].GetR())>240.) continue;
5921 if (TMath::Abs(kinks[irow].GetR())<110.) continue;
5922 Float_t quality = TMath::Abs(kinks[irow].GetAngle(2))/(3.+TMath::Abs(kinks[irow].GetR()-param0[irow].GetX()));
5923 if ( quality > maxchange){
5924 maxchange = quality;
5931 if (index<0) return 0;
5933 Int_t row0 = GetRowNumber(kinks[index].GetR()); //row 0 estimate
5934 seed0 = new AliTPCseed(param0[index]);
5935 seed1 = new AliTPCseed(param1[index]);
5936 seed0->Reset(kFALSE);
5937 seed1->Reset(kFALSE);
5938 seed0->ResetCovariance(10.);
5939 seed1->ResetCovariance(10.);
5940 FollowProlongation(*seed0,0);
5941 FollowBackProlongation(*seed1,158);
5942 mother = *seed0; // backup mother at position 0
5943 seed0->Reset(kFALSE);
5944 seed1->Reset(kFALSE);
5945 seed0->ResetCovariance(10.);
5946 seed1->ResetCovariance(10.);
5948 first = TMath::Max(row0-20,0);
5949 last = TMath::Min(row0+20,158);
5951 for (Int_t irow=0; irow<20;irow++){
5952 rows[irow] = first +((last-first)*irow)/19;
5954 // store parameters along the track
5956 for (Int_t irow=0;irow<20;irow++){
5957 FollowBackProlongation(*seed0, rows[irow]);
5958 FollowProlongation(*seed1,rows[19-irow]);
5959 param0[irow] = *seed0;
5960 param1[19-irow] = *seed1;
5964 for (Int_t irow=0; irow<19;irow++){
5965 kinks[irow].SetMother(param0[irow]);
5966 kinks[irow].SetDaughter(param1[irow]);
5967 // param0[irow].Dump();
5968 //param1[irow].Dump();
5969 kinks[irow].Update();
5972 // choose kink with biggest change of angle
5975 for (Int_t irow=0;irow<20;irow++){
5976 if (TMath::Abs(kinks[irow].GetR())>250.) continue;
5977 if (TMath::Abs(kinks[irow].GetR())<90.) continue;
5978 Float_t quality = TMath::Abs(kinks[irow].GetAngle(2))/(3.+TMath::Abs(kinks[irow].GetR()-param0[irow].GetX()));
5979 if ( quality > maxchange){
5980 maxchange = quality;
5987 if (index==-1 || param0[index].GetNumberOfClusters()+param1[index].GetNumberOfClusters()<100){
5992 // Float_t anglesigma = TMath::Sqrt(param0[index].fC22+param0[index].fC33+param1[index].fC22+param1[index].fC33);
5994 kink.SetMother(param0[index]);
5995 kink.SetDaughter(param1[index]);
5997 row0 = GetRowNumber(kink.GetR());
5998 kink.SetTPCRow0(row0);
5999 kink.SetLabel(CookLabel(seed0,0.5,0,row0),0);
6000 kink.SetLabel(CookLabel(seed1,0.5,row0,158),1);
6001 kink.SetIndex(-10,0);
6002 kink.SetIndex(int(param0[index].GetNumberOfClusters()+param1[index].GetNumberOfClusters()),1);
6003 kink.SetTPCncls(param0[index].GetNumberOfClusters(),0);
6004 kink.SetTPCncls(param1[index].GetNumberOfClusters(),1);
6007 // new (&mother) AliTPCseed(param0[index]);
6008 daughter = param1[index];
6009 daughter.SetLabel(kink.GetLabel(1));
6010 param0[index].Reset(kTRUE);
6011 FollowProlongation(param0[index],0);
6012 mother = param0[index];
6013 mother.SetLabel(kink.GetLabel(0));
6023 AliTPCseed* AliTPCtrackerMI::ReSeed(AliTPCseed *t)
6026 // reseed - refit - track
6029 // Int_t last = fSectors->GetNRows()-1;
6031 if (fSectors == fOuterSec){
6032 first = TMath::Max(first, t->GetFirstPoint()-fInnerSec->GetNRows());
6036 first = t->GetFirstPoint();
6038 AliTPCseed * seed = MakeSeed(t,0.1,0.5,0.9);
6039 FollowBackProlongation(*t,fSectors->GetNRows()-1);
6041 FollowProlongation(*t,first);
6051 //_____________________________________________________________________________
6052 Int_t AliTPCtrackerMI::ReadSeeds(const TFile *inp) {
6053 //-----------------------------------------------------------------
6054 // This function reades track seeds.
6055 //-----------------------------------------------------------------
6056 TDirectory *savedir=gDirectory;
6058 TFile *in=(TFile*)inp;
6059 if (!in->IsOpen()) {
6060 cerr<<"AliTPCtrackerMI::ReadSeeds(): input file is not open !\n";
6065 TTree *seedTree=(TTree*)in->Get("Seeds");
6067 cerr<<"AliTPCtrackerMI::ReadSeeds(): ";
6068 cerr<<"can't get a tree with track seeds !\n";
6071 AliTPCtrack *seed=new AliTPCtrack;
6072 seedTree->SetBranchAddress("tracks",&seed);
6074 if (fSeeds==0) fSeeds=new TObjArray(15000);
6076 Int_t n=(Int_t)seedTree->GetEntries();
6077 for (Int_t i=0; i<n; i++) {
6078 seedTree->GetEvent(i);
6079 fSeeds->AddLast(new AliTPCseed(*seed/*,seed->GetAlpha()*/));
6088 Int_t AliTPCtrackerMI::Clusters2Tracks (AliESDEvent *esd)
6091 if (fSeeds) DeleteSeeds();
6094 if (!fSeeds) return 1;
6101 //_____________________________________________________________________________
6102 Int_t AliTPCtrackerMI::Clusters2Tracks() {
6103 //-----------------------------------------------------------------
6104 // This is a track finder.
6105 //-----------------------------------------------------------------
6106 TDirectory *savedir=gDirectory;
6110 fSeeds = Tracking();
6113 Info("Clusters2Tracks","Time for tracking: \t");timer.Print();timer.Start();
6115 //activate again some tracks
6116 for (Int_t i=0; i<fSeeds->GetEntriesFast(); i++) {
6117 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
6119 Int_t nc=t.GetNumberOfClusters();
6121 delete fSeeds->RemoveAt(i);
6125 if (pt->GetRemoval()==10) {
6126 if (pt->GetDensityFirst(20)>0.8 || pt->GetDensityFirst(30)>0.8 || pt->GetDensityFirst(40)>0.7)
6127 pt->Desactivate(10); // make track again active
6129 pt->Desactivate(20);
6130 delete fSeeds->RemoveAt(i);
6135 RemoveUsed2(fSeeds,0.85,0.85,0);
6136 if (AliTPCReconstructor::GetRecoParam()->GetDoKinks()) FindKinks(fSeeds,fEvent);
6137 //FindCurling(fSeeds, fEvent,0);
6138 if (AliTPCReconstructor::StreamLevel()>2) FindMultiMC(fSeeds, fEvent,0); // find multi found tracks
6139 RemoveUsed2(fSeeds,0.5,0.4,20);
6140 FindSplitted(fSeeds, fEvent,0); // find multi found tracks
6143 // // refit short tracks
6145 Int_t nseed=fSeeds->GetEntriesFast();
6148 for (Int_t i=0; i<nseed; i++) {
6149 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
6151 Int_t nc=t.GetNumberOfClusters();
6153 delete fSeeds->RemoveAt(i);
6156 CookLabel(pt,0.1); //For comparison only
6157 //if ((pt->IsActive() || (pt->fRemoval==10) )&& nc>50 &&pt->GetNumberOfClusters()>0.4*pt->fNFoundable){
6158 if ((pt->IsActive() || (pt->GetRemoval()==10) )){
6160 if (fDebug>0) cerr<<found<<'\r';
6164 delete fSeeds->RemoveAt(i);
6168 //RemoveOverlap(fSeeds,0.99,7,kTRUE);
6170 //RemoveUsed(fSeeds,0.9,0.9,6);
6172 nseed=fSeeds->GetEntriesFast();
6174 for (Int_t i=0; i<nseed; i++) {
6175 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
6177 Int_t nc=t.GetNumberOfClusters();
6179 delete fSeeds->RemoveAt(i);
6183 t.CookdEdx(0.02,0.6);
6184 // CheckKinkPoint(&t,0.05);
6185 //if ((pt->IsActive() || (pt->fRemoval==10) )&& nc>50 &&pt->GetNumberOfClusters()>0.4*pt->fNFoundable){
6186 if ((pt->IsActive() || (pt->GetRemoval()==10) )){
6194 delete fSeeds->RemoveAt(i);
6195 //AliTPCseed * seed1 = ReSeed(pt,0.05,0.5,1);
6197 // FollowProlongation(*seed1,0);
6198 // Int_t n = seed1->GetNumberOfClusters();
6199 // printf("fP4\t%f\t%f\n",seed1->GetC(),pt->GetC());
6200 // printf("fN\t%d\t%d\n", seed1->GetNumberOfClusters(),pt->GetNumberOfClusters());
6203 //AliTPCseed * seed2 = ReSeed(pt,0.95,0.5,0.05);
6207 SortTracks(fSeeds, 1);
6211 PrepareForBackProlongation(fSeeds,5.);
6212 PropagateBack(fSeeds);
6213 printf("Time for back propagation: \t");timer.Print();timer.Start();
6217 PrepareForProlongation(fSeeds,5.);
6218 PropagateForward2(fSeeds);
6220 printf("Time for FORWARD propagation: \t");timer.Print();timer.Start();
6221 // RemoveUsed(fSeeds,0.7,0.7,6);
6222 //RemoveOverlap(fSeeds,0.9,7,kTRUE);
6224 nseed=fSeeds->GetEntriesFast();
6226 for (Int_t i=0; i<nseed; i++) {
6227 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
6229 Int_t nc=t.GetNumberOfClusters();
6231 delete fSeeds->RemoveAt(i);
6234 t.CookdEdx(0.02,0.6);
6235 // CookLabel(pt,0.1); //For comparison only
6236 //if ((pt->IsActive() || (pt->fRemoval==10) )&& nc>50 &&pt->GetNumberOfClusters()>0.4*pt->fNFoundable){
6237 if ((pt->IsActive() || (pt->fRemoval==10) )){
6238 cerr<<found++<<'\r';
6241 delete fSeeds->RemoveAt(i);
6246 // fNTracks = found;
6248 Info("Clusters2Tracks","Time for overlap removal, track writing and dedx cooking: \t"); timer.Print();timer.Start();
6251 // cerr<<"Number of found tracks : "<<"\t"<<found<<endl;
6252 Info("Clusters2Tracks","Number of found tracks %d",found);
6254 // UnloadClusters();
6259 void AliTPCtrackerMI::Tracking(TObjArray * arr)
6262 // tracking of the seeds
6265 fSectors = fOuterSec;
6266 ParallelTracking(arr,150,63);
6267 fSectors = fOuterSec;
6268 ParallelTracking(arr,63,0);
6271 TObjArray * AliTPCtrackerMI::Tracking(Int_t seedtype, Int_t i1, Int_t i2, Float_t cuts[4], Float_t dy, Int_t dsec)
6276 TObjArray * arr = new TObjArray;
6278 fSectors = fOuterSec;
6281 for (Int_t sec=0;sec<fkNOS;sec++){
6282 if (seedtype==3) MakeSeeds3(arr,sec,i1,i2,cuts,dy, dsec);
6283 if (seedtype==4) MakeSeeds5(arr,sec,i1,i2,cuts,dy);
6284 if (seedtype==2) MakeSeeds2(arr,sec,i1,i2,cuts,dy);
6287 Info("Tracking","\nSeeding - %d\t%d\t%d\t%d\n",seedtype,i1,i2,arr->GetEntriesFast());
6299 TObjArray * AliTPCtrackerMI::Tracking()
6303 if (AliTPCReconstructor::GetRecoParam()->GetSpecialSeeding()) return TrackingSpecial();
6306 Int_t nup=fOuterSec->GetNRows()+fInnerSec->GetNRows();
6308 TObjArray * seeds = new TObjArray;
6317 Float_t fnumber = 3.0;
6318 Float_t fdensity = 3.0;
6323 for (Int_t delta = 0; delta<18; delta+=6){
6327 arr = Tracking(3,nup-1-delta,nup-1-delta-gap,cuts,-1,1);
6328 SumTracks(seeds,arr);
6329 SignClusters(seeds,fnumber,fdensity);
6331 for (Int_t i=2;i<6;i+=2){
6332 // seed high pt tracks
6335 arr = Tracking(3,nup-i-delta,nup-i-delta-gap,cuts,-1,0);
6336 SumTracks(seeds,arr);
6337 SignClusters(seeds,fnumber,fdensity);
6342 // RemoveUsed(seeds,0.9,0.9,1);
6343 // UnsignClusters();
6344 // SignClusters(seeds,fnumber,fdensity);
6348 for (Int_t delta = 20; delta<120; delta+=10){
6350 // seed high pt tracks
6354 arr = Tracking(3,nup-delta,nup-delta-gap,cuts,-1);
6355 SumTracks(seeds,arr);
6356 SignClusters(seeds,fnumber,fdensity);
6361 arr = Tracking(3,nup-delta-5,nup-delta-5-gap,cuts,-1);
6362 SumTracks(seeds,arr);
6363 SignClusters(seeds,fnumber,fdensity);
6374 Info("Tracking()","\n\nPrimary seeding\t%d\n\n",seeds->GetEntriesFast());
6378 // RemoveUsed(seeds,0.75,0.75,1);
6380 //SignClusters(seeds,fnumber,fdensity);
6389 arr = Tracking(4,nup-1,nup-1-gap,cuts,-1);
6390 SumTracks(seeds,arr);
6391 SignClusters(seeds,fnumber,fdensity);
6393 arr = Tracking(4,nup-2,nup-2-gap,cuts,-1);
6394 SumTracks(seeds,arr);
6395 SignClusters(seeds,fnumber,fdensity);
6397 arr = Tracking(4,nup-3,nup-3-gap,cuts,-1);
6398 SumTracks(seeds,arr);
6399 SignClusters(seeds,fnumber,fdensity);
6403 for (Int_t delta = 3; delta<30; delta+=5){
6409 arr = Tracking(4,nup-1-delta,nup-1-delta-gap,cuts,-1);
6410 SumTracks(seeds,arr);
6411 SignClusters(seeds,fnumber,fdensity);
6413 arr = Tracking(4,nup-3-delta,nup-5-delta-gap,cuts,4);
6414 SumTracks(seeds,arr);
6415 SignClusters(seeds,fnumber,fdensity);
6426 Int_t fLastSeedRowSec=AliTPCReconstructor::GetRecoParam()->GetLastSeedRowSec();
6429 for (Int_t delta = 30; delta<fLastSeedRowSec; delta+=10){
6435 arr = Tracking(4,nup-1-delta,nup-1-delta-gap,cuts,-1);
6436 SumTracks(seeds,arr);
6437 SignClusters(seeds,fnumber,fdensity);
6439 arr = Tracking(4,nup-5-delta,nup-5-delta-gap,cuts,5 );
6440 SumTracks(seeds,arr);
6441 SignClusters(seeds,fnumber,fdensity);
6445 Info("Tracking()","\n\nSecondary seeding\t%d\n\n",seeds->GetEntriesFast());
6456 TObjArray * AliTPCtrackerMI::TrackingSpecial()
6459 // seeding adjusted for laser and cosmic tests - short tracks with big inclination angle
6460 // no primary vertex seeding tried
6464 Int_t nup=fOuterSec->GetNRows()+fInnerSec->GetNRows();
6466 TObjArray * seeds = new TObjArray;
6471 Float_t fnumber = 3.0;
6472 Float_t fdensity = 3.0;
6475 cuts[0] = AliTPCReconstructor::GetRecoParam()->GetMaxC(); // max curvature
6476 cuts[1] = 3.5; // max tan(phi) angle for seeding
6477 cuts[2] = 3.; // not used (cut on z primary vertex)
6478 cuts[3] = 3.5; // max tan(theta) angle for seeding
6480 for (Int_t delta = 0; nup-delta-gap-1>0; delta+=3){
6482 arr = Tracking(4,nup-1-delta,nup-1-delta-gap,cuts,-1);
6483 SumTracks(seeds,arr);
6484 SignClusters(seeds,fnumber,fdensity);
6488 Info("Tracking()","\n\nSecondary seeding\t%d\n\n",seeds->GetEntriesFast());
6499 void AliTPCtrackerMI::SumTracks(TObjArray *arr1,TObjArray *&arr2) const
6502 //sum tracks to common container
6503 //remove suspicious tracks
6504 Int_t nseed = arr2->GetEntriesFast();
6505 for (Int_t i=0;i<nseed;i++){
6506 AliTPCseed *pt=(AliTPCseed*)arr2->UncheckedAt(i);
6509 // remove tracks with too big curvature
6511 if (TMath::Abs(pt->GetC())>AliTPCReconstructor::GetRecoParam()->GetMaxC()){
6512 delete arr2->RemoveAt(i);
6515 // REMOVE VERY SHORT TRACKS
6516 if (pt->GetNumberOfClusters()<20){
6517 delete arr2->RemoveAt(i);
6520 // NORMAL ACTIVE TRACK
6521 if (pt->IsActive()){
6522 arr1->AddLast(arr2->RemoveAt(i));
6525 //remove not usable tracks
6526 if (pt->GetRemoval()!=10){
6527 delete arr2->RemoveAt(i);
6531 // ENABLE ONLY ENOUGH GOOD STOPPED TRACKS
6532 if (pt->GetDensityFirst(20)>0.8 || pt->GetDensityFirst(30)>0.8 || pt->GetDensityFirst(40)>0.7)
6533 arr1->AddLast(arr2->RemoveAt(i));
6535 delete arr2->RemoveAt(i);
6539 delete arr2; arr2 = 0;
6544 void AliTPCtrackerMI::ParallelTracking(TObjArray * arr, Int_t rfirst, Int_t rlast)
6547 // try to track in parralel
6549 Int_t nseed=arr->GetEntriesFast();
6550 //prepare seeds for tracking
6551 for (Int_t i=0; i<nseed; i++) {
6552 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i), &t=*pt;
6554 if (!t.IsActive()) continue;
6555 // follow prolongation to the first layer
6556 if ( (fSectors ==fInnerSec) || (t.GetFirstPoint()-fParam->GetNRowLow()>rfirst+1) )
6557 FollowProlongation(t, rfirst+1);
6562 for (Int_t nr=rfirst; nr>=rlast; nr--){
6563 if (nr<fInnerSec->GetNRows())
6564 fSectors = fInnerSec;
6566 fSectors = fOuterSec;
6567 // make indexes with the cluster tracks for given
6569 // find nearest cluster
6570 for (Int_t i=0; i<nseed; i++) {
6571 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i), &t=*pt;
6573 if (nr==80) pt->UpdateReference();
6574 if (!pt->IsActive()) continue;
6575 // if ( (fSectors ==fOuterSec) && (pt->fFirstPoint-fParam->GetNRowLow())<nr) continue;
6576 if (pt->GetRelativeSector()>17) {
6579 UpdateClusters(t,nr);
6581 // prolonagate to the nearest cluster - if founded
6582 for (Int_t i=0; i<nseed; i++) {
6583 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
6585 if (!pt->IsActive()) continue;
6586 // if ((fSectors ==fOuterSec) && (pt->fFirstPoint-fParam->GetNRowLow())<nr) continue;
6587 if (pt->GetRelativeSector()>17) {
6590 FollowToNextCluster(*pt,nr);
6595 void AliTPCtrackerMI::PrepareForBackProlongation(TObjArray * arr,Float_t fac) const
6599 // if we use TPC track itself we have to "update" covariance
6601 Int_t nseed= arr->GetEntriesFast();
6602 for (Int_t i=0;i<nseed;i++){
6603 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
6607 //rotate to current local system at first accepted point
6608 Int_t index = pt->GetClusterIndex2(pt->GetFirstPoint());
6609 Int_t sec = (index&0xff000000)>>24;
6611 Float_t angle1 = fInnerSec->GetAlpha()*sec+fInnerSec->GetAlphaShift();
6612 if (angle1>TMath::Pi())
6613 angle1-=2.*TMath::Pi();
6614 Float_t angle2 = pt->GetAlpha();
6616 if (TMath::Abs(angle1-angle2)>0.001){
6617 pt->Rotate(angle1-angle2);
6618 //angle2 = pt->GetAlpha();
6619 //pt->fRelativeSector = pt->GetAlpha()/fInnerSec->GetAlpha();
6620 //if (pt->GetAlpha()<0)
6621 // pt->fRelativeSector+=18;
6622 //sec = pt->fRelativeSector;
6631 void AliTPCtrackerMI::PrepareForProlongation(TObjArray * arr, Float_t fac) const
6635 // if we use TPC track itself we have to "update" covariance
6637 Int_t nseed= arr->GetEntriesFast();
6638 for (Int_t i=0;i<nseed;i++){
6639 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
6642 pt->SetFirstPoint(pt->GetLastPoint());
6650 Int_t AliTPCtrackerMI::PropagateBack(TObjArray * arr)
6653 // make back propagation
6655 Int_t nseed= arr->GetEntriesFast();
6656 for (Int_t i=0;i<nseed;i++){
6657 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
6658 if (pt&& pt->GetKinkIndex(0)<=0) {
6659 //AliTPCseed *pt2 = new AliTPCseed(*pt);
6660 fSectors = fInnerSec;
6661 //FollowBackProlongation(*pt,fInnerSec->GetNRows()-1);
6662 //fSectors = fOuterSec;
6663 FollowBackProlongation(*pt,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1);
6664 //if (pt->GetNumberOfClusters()<(pt->fEsd->GetTPCclusters(0)) ){
6665 // Error("PropagateBack","Not prolonged track %d",pt->GetLabel());
6666 // FollowBackProlongation(*pt2,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1);
6669 if (pt&& pt->GetKinkIndex(0)>0) {
6670 AliESDkink * kink = fEvent->GetKink(pt->GetKinkIndex(0)-1);
6671 pt->SetFirstPoint(kink->GetTPCRow0());
6672 fSectors = fInnerSec;
6673 FollowBackProlongation(*pt,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1);
6681 Int_t AliTPCtrackerMI::PropagateForward2(TObjArray * arr)
6684 // make forward propagation
6686 Int_t nseed= arr->GetEntriesFast();
6688 for (Int_t i=0;i<nseed;i++){
6689 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
6691 FollowProlongation(*pt,0);
6700 Int_t AliTPCtrackerMI::PropagateForward()
6703 // propagate track forward
6705 Int_t nseed = fSeeds->GetEntriesFast();
6706 for (Int_t i=0;i<nseed;i++){
6707 AliTPCseed *pt = (AliTPCseed*)fSeeds->UncheckedAt(i);
6709 AliTPCseed &t = *pt;
6710 Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
6711 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
6712 if (alpha < 0. ) alpha += 2.*TMath::Pi();
6713 t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN);
6717 fSectors = fOuterSec;
6718 ParallelTracking(fSeeds,fOuterSec->GetNRows()+fInnerSec->GetNRows()-1,fInnerSec->GetNRows());
6719 fSectors = fInnerSec;
6720 ParallelTracking(fSeeds,fInnerSec->GetNRows()-1,0);
6729 Int_t AliTPCtrackerMI::PropagateBack(AliTPCseed * pt, Int_t row0, Int_t row1)
6732 // make back propagation, in between row0 and row1
6736 fSectors = fInnerSec;
6739 if (row1<fSectors->GetNRows())
6742 r1 = fSectors->GetNRows()-1;
6744 if (row0<fSectors->GetNRows()&& r1>0 )
6745 FollowBackProlongation(*pt,r1);
6746 if (row1<=fSectors->GetNRows())
6749 r1 = row1 - fSectors->GetNRows();
6750 if (r1<=0) return 0;
6751 if (r1>=fOuterSec->GetNRows()) return 0;
6752 fSectors = fOuterSec;
6753 return FollowBackProlongation(*pt,r1);
6761 void AliTPCtrackerMI::GetShape(AliTPCseed * seed, Int_t row)
6765 AliTPCClusterParam * clparam = AliTPCcalibDB::Instance()->GetClusterParam();
6766 Float_t zdrift = TMath::Abs((fParam->GetZLength(0)-TMath::Abs(seed->GetZ())));
6767 Int_t type = (seed->GetSector() < fParam->GetNSector()/2) ? 0: (row>126) ? 1:2;
6768 Double_t angulary = seed->GetSnp();
6769 angulary = angulary*angulary/(1.-angulary*angulary);
6770 Double_t angularz = seed->GetTgl()*seed->GetTgl()*(1.+angulary);
6772 Double_t sigmay = clparam->GetRMS0(0,type,zdrift,TMath::Sqrt(TMath::Abs(angulary)));
6773 Double_t sigmaz = clparam->GetRMS0(1,type,zdrift,TMath::Sqrt(TMath::Abs(angularz)));
6774 seed->SetCurrentSigmaY2(sigmay*sigmay);
6775 seed->SetCurrentSigmaZ2(sigmaz*sigmaz);
6776 // Float_t sd2 = TMath::Abs((fParam->GetZLength(0)-TMath::Abs(seed->GetZ())))*fParam->GetDiffL()*fParam->GetDiffL();
6777 // // Float_t padlength = fParam->GetPadPitchLength(seed->fSector);
6778 // Float_t padlength = GetPadPitchLength(row);
6780 // Float_t sresy = (seed->GetSector() < fParam->GetNSector()/2) ? 0.2 :0.3;
6781 // seed->SetCurrentSigmaY2(sd2+padlength*padlength*angulary/12.+sresy*sresy);
6783 // Float_t sresz = fParam->GetZSigma();
6784 // seed->SetCurrentSigmaZ2(sd2+padlength*padlength*angularz*angularz*(1+angulary)/12.+sresz*sresz);
6786 Float_t wy = GetSigmaY(seed);
6787 Float_t wz = GetSigmaZ(seed);
6790 if (TMath::Abs(wy/seed->fCurrentSigmaY2-1)>0.0001 || TMath::Abs(wz/seed->fCurrentSigmaZ2-1)>0.0001 ){
6791 printf("problem\n");
6798 //__________________________________________________________________________
6799 void AliTPCtrackerMI::CookLabel(AliKalmanTrack *tk, Float_t wrong) const {
6800 //--------------------------------------------------------------------
6801 //This function "cooks" a track label. If label<0, this track is fake.
6802 //--------------------------------------------------------------------
6803 AliTPCseed * t = dynamic_cast<AliTPCseed*>(tk);
6805 printf("%s:%d wrong type \n",(char*)__FILE__,__LINE__);
6809 Int_t noc=t->GetNumberOfClusters();
6811 //printf("\nnot founded prolongation\n\n\n");
6817 AliTPCclusterMI *clusters[160];
6819 for (Int_t i=0;i<160;i++) {
6826 for (i=0; i<160 && current<noc; i++) {
6828 Int_t index=t->GetClusterIndex2(i);
6829 if (index<=0) continue;
6830 if (index&0x8000) continue;
6832 //clusters[current]=GetClusterMI(index);
6833 if (t->GetClusterPointer(i)){
6834 clusters[current]=t->GetClusterPointer(i);
6840 Int_t lab=123456789;
6841 for (i=0; i<noc; i++) {
6842 AliTPCclusterMI *c=clusters[i];
6844 lab=TMath::Abs(c->GetLabel(0));
6846 for (j=0; j<noc; j++) if (lb[j]==lab || mx[j]==0) break;
6852 for (i=0; i<noc; i++) if (mx[i]>max) {max=mx[i]; lab=lb[i];}
6854 for (i=0; i<noc; i++) {
6855 AliTPCclusterMI *c=clusters[i];
6857 if (TMath::Abs(c->GetLabel(1)) == lab ||
6858 TMath::Abs(c->GetLabel(2)) == lab ) max++;
6861 if ((1.- Float_t(max)/noc) > wrong) lab=-lab;
6864 Int_t tail=Int_t(0.10*noc);
6867 for (i=1; i<=160&&ind<tail; i++) {
6868 // AliTPCclusterMI *c=clusters[noc-i];
6869 AliTPCclusterMI *c=clusters[i];
6871 if (lab == TMath::Abs(c->GetLabel(0)) ||
6872 lab == TMath::Abs(c->GetLabel(1)) ||
6873 lab == TMath::Abs(c->GetLabel(2))) max++;
6876 if (max < Int_t(0.5*tail)) lab=-lab;
6883 //delete[] clusters;
6887 //__________________________________________________________________________
6888 Int_t AliTPCtrackerMI::CookLabel(AliTPCseed *t, Float_t wrong,Int_t first, Int_t last) const {
6889 //--------------------------------------------------------------------
6890 //This function "cooks" a track label. If label<0, this track is fake.
6891 //--------------------------------------------------------------------
6892 Int_t noc=t->GetNumberOfClusters();
6894 //printf("\nnot founded prolongation\n\n\n");
6900 AliTPCclusterMI *clusters[160];
6902 for (Int_t i=0;i<160;i++) {
6909 for (i=0; i<160 && current<noc; i++) {
6910 if (i<first) continue;
6911 if (i>last) continue;
6912 Int_t index=t->GetClusterIndex2(i);
6913 if (index<=0) continue;
6914 if (index&0x8000) continue;
6916 //clusters[current]=GetClusterMI(index);
6917 if (t->GetClusterPointer(i)){
6918 clusters[current]=t->GetClusterPointer(i);
6923 if (noc<5) return -1;
6924 Int_t lab=123456789;
6925 for (i=0; i<noc; i++) {
6926 AliTPCclusterMI *c=clusters[i];
6928 lab=TMath::Abs(c->GetLabel(0));
6930 for (j=0; j<noc; j++) if (lb[j]==lab || mx[j]==0) break;
6936 for (i=0; i<noc; i++) if (mx[i]>max) {max=mx[i]; lab=lb[i];}
6938 for (i=0; i<noc; i++) {
6939 AliTPCclusterMI *c=clusters[i];
6941 if (TMath::Abs(c->GetLabel(1)) == lab ||
6942 TMath::Abs(c->GetLabel(2)) == lab ) max++;
6945 if ((1.- Float_t(max)/noc) > wrong) lab=-lab;
6948 Int_t tail=Int_t(0.10*noc);
6951 for (i=1; i<=160&&ind<tail; i++) {
6952 // AliTPCclusterMI *c=clusters[noc-i];
6953 AliTPCclusterMI *c=clusters[i];
6955 if (lab == TMath::Abs(c->GetLabel(0)) ||
6956 lab == TMath::Abs(c->GetLabel(1)) ||
6957 lab == TMath::Abs(c->GetLabel(2))) max++;
6960 if (max < Int_t(0.5*tail)) lab=-lab;
6963 // t->SetLabel(lab);
6967 //delete[] clusters;
6971 Int_t AliTPCtrackerMI::GetRowNumber(Double_t x[3]) const
6973 //return pad row number for given x vector
6974 Float_t phi = TMath::ATan2(x[1],x[0]);
6975 if(phi<0) phi=2.*TMath::Pi()+phi;
6976 // Get the local angle in the sector philoc
6977 const Float_t kRaddeg = 180/3.14159265358979312;
6978 Float_t phiangle = (Int_t (phi*kRaddeg/20.) + 0.5)*20./kRaddeg;
6979 Double_t localx = x[0]*TMath::Cos(phiangle)-x[1]*TMath::Sin(phiangle);
6980 return GetRowNumber(localx);
6985 void AliTPCtrackerMI::MakeBitmaps(AliTPCseed *t)
6987 //-----------------------------------------------------------------------
6988 // Fill the cluster and sharing bitmaps of the track
6989 //-----------------------------------------------------------------------
6991 Int_t firstpoint = 0;
6992 Int_t lastpoint = 159;
6993 AliTPCTrackerPoint *point;
6995 for (int iter=firstpoint; iter<lastpoint; iter++) {
6996 point = t->GetTrackPoint(iter);
6998 t->SetClusterMapBit(iter, kTRUE);
6999 if (point->IsShared())
7000 t->SetSharedMapBit(iter,kTRUE);
7002 t->SetSharedMapBit(iter, kFALSE);
7005 t->SetClusterMapBit(iter, kFALSE);
7006 t->SetSharedMapBit(iter, kFALSE);
7013 void AliTPCtrackerMI::AddCovariance(AliTPCseed * seed){
7015 // Adding systematic error
7016 // !!!! the systematic error for element 4 is in 1/cm not in pt
7018 const Double_t *param = AliTPCReconstructor::GetRecoParam()->GetSystematicError();
7020 for (Int_t i=0;i<15;i++) covar[i]=0;
7026 covar[0] = param[0]*param[0];
7027 covar[2] = param[1]*param[1];
7028 covar[5] = param[2]*param[2];
7029 covar[9] = param[3]*param[3];
7030 Double_t facC = AliTracker::GetBz()*kB2C;
7031 covar[14]= param[4]*param[4]*facC*facC;
7032 seed->AddCovariance(covar);