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>=fkParam->GetNInnerSector()) track->SetRow(track->GetRow()+fkParam->GetNRowLow());
226 track->SetClusterIndex2(track->GetRow(), i);
227 //track->fFirstPoint = row;
228 //if ( track->fLastPoint<row) track->fLastPoint =row;
229 // if (track->fRow<0 || track->fRow>160) {
230 // printf("problem\n");
232 if (track->GetFirstPoint()>track->GetRow())
233 track->SetFirstPoint(track->GetRow());
234 if (track->GetLastPoint()<track->GetRow())
235 track->SetLastPoint(track->GetRow());
238 track->SetClusterPointer(track->GetRow(),c);
241 Double_t angle2 = track->GetSnp()*track->GetSnp();
243 //SET NEW Track Point
245 if (angle2<1) //PH sometimes angle2 is very big. To be investigated...
247 angle2 = TMath::Sqrt(angle2/(1-angle2));
248 AliTPCTrackerPoint &point =*(track->GetTrackPoint(track->GetRow()));
250 point.SetSigmaY(c->GetSigmaY2()/track->GetCurrentSigmaY2());
251 point.SetSigmaZ(c->GetSigmaZ2()/track->GetCurrentSigmaZ2());
252 point.SetErrY(sqrt(track->GetErrorY2()));
253 point.SetErrZ(sqrt(track->GetErrorZ2()));
255 point.SetX(track->GetX());
256 point.SetY(track->GetY());
257 point.SetZ(track->GetZ());
258 point.SetAngleY(angle2);
259 point.SetAngleZ(track->GetTgl());
260 if (point.IsShared()){
261 track->SetErrorY2(track->GetErrorY2()*4);
262 track->SetErrorZ2(track->GetErrorZ2()*4);
266 Double_t chi2 = track->GetPredictedChi2(track->GetCurrentCluster());
268 // track->SetErrorY2(track->GetErrorY2()*1.3);
269 // track->SetErrorY2(track->GetErrorY2()+0.01);
270 // track->SetErrorZ2(track->GetErrorZ2()*1.3);
271 // track->SetErrorZ2(track->GetErrorZ2()+0.005);
273 if (accept>0) return 0;
274 if (track->GetNumberOfClusters()%20==0){
275 // if (track->fHelixIn){
276 // TClonesArray & larr = *(track->fHelixIn);
277 // Int_t ihelix = larr.GetEntriesFast();
278 // new(larr[ihelix]) AliHelix(*track) ;
281 track->SetNoCluster(0);
282 return track->Update(c,chi2,i);
287 Int_t AliTPCtrackerMI::AcceptCluster(AliTPCseed * seed, AliTPCclusterMI * cluster)
290 // decide according desired precision to accept given
291 // cluster for tracking
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 (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 (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*/)
443 //------------------------------
445 //--------------------------------------------------------------
448 //_____________________________________________________________________________
449 AliTPCtrackerMI::~AliTPCtrackerMI() {
450 //------------------------------------------------------------------
451 // TPC tracker destructor
452 //------------------------------------------------------------------
459 if (fDebugStreamer) delete fDebugStreamer;
463 void AliTPCtrackerMI::FillESD(const TObjArray* arr)
467 //fill esds using updated tracks
469 // write tracks to the event
470 // store index of the track
471 Int_t nseed=arr->GetEntriesFast();
472 //FindKinks(arr,fEvent);
473 for (Int_t i=0; i<nseed; i++) {
474 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
478 if (AliTPCReconstructor::StreamLevel()>1) {
479 (*fDebugStreamer)<<"Track0"<<
483 // pt->PropagateTo(fkParam->GetInnerRadiusLow());
484 if (pt->GetKinkIndex(0)<=0){ //don't propagate daughter tracks
485 pt->PropagateTo(fkParam->GetInnerRadiusLow());
488 if (( pt->GetPoints()[2]- pt->GetPoints()[0])>5 && pt->GetPoints()[3]>0.8){
490 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
491 iotrack.SetTPCPoints(pt->GetPoints());
492 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
493 iotrack.SetV0Indexes(pt->GetV0Indexes());
494 // iotrack.SetTPCpid(pt->fTPCr);
495 //iotrack.SetTPCindex(i);
496 fEvent->AddTrack(&iotrack);
500 if ( (pt->GetNumberOfClusters()>70)&& (Float_t(pt->GetNumberOfClusters())/Float_t(pt->GetNFoundable()))>0.55) {
502 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
503 iotrack.SetTPCPoints(pt->GetPoints());
504 //iotrack.SetTPCindex(i);
505 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
506 iotrack.SetV0Indexes(pt->GetV0Indexes());
507 // iotrack.SetTPCpid(pt->fTPCr);
508 fEvent->AddTrack(&iotrack);
512 // short tracks - maybe decays
514 if ( (pt->GetNumberOfClusters()>30) && (Float_t(pt->GetNumberOfClusters())/Float_t(pt->GetNFoundable()))>0.70) {
515 Int_t found,foundable,shared;
516 pt->GetClusterStatistic(0,60,found, foundable,shared,kFALSE);
517 if ( (found>20) && (pt->GetNShared()/float(pt->GetNumberOfClusters())<0.2)){
519 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
520 //iotrack.SetTPCindex(i);
521 iotrack.SetTPCPoints(pt->GetPoints());
522 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
523 iotrack.SetV0Indexes(pt->GetV0Indexes());
524 //iotrack.SetTPCpid(pt->fTPCr);
525 fEvent->AddTrack(&iotrack);
530 if ( (pt->GetNumberOfClusters()>20) && (Float_t(pt->GetNumberOfClusters())/Float_t(pt->GetNFoundable()))>0.8) {
531 Int_t found,foundable,shared;
532 pt->GetClusterStatistic(0,60,found, foundable,shared,kFALSE);
533 if (found<20) continue;
534 if (pt->GetNShared()/float(pt->GetNumberOfClusters())>0.2) continue;
537 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
538 iotrack.SetTPCPoints(pt->GetPoints());
539 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
540 iotrack.SetV0Indexes(pt->GetV0Indexes());
541 //iotrack.SetTPCpid(pt->fTPCr);
542 //iotrack.SetTPCindex(i);
543 fEvent->AddTrack(&iotrack);
546 // short tracks - secondaties
548 if ( (pt->GetNumberOfClusters()>30) ) {
549 Int_t found,foundable,shared;
550 pt->GetClusterStatistic(128,158,found, foundable,shared,kFALSE);
551 if ( (found>20) && (pt->GetNShared()/float(pt->GetNumberOfClusters())<0.2) &&float(found)/float(foundable)>0.8){
553 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
554 iotrack.SetTPCPoints(pt->GetPoints());
555 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
556 iotrack.SetV0Indexes(pt->GetV0Indexes());
557 //iotrack.SetTPCpid(pt->fTPCr);
558 //iotrack.SetTPCindex(i);
559 fEvent->AddTrack(&iotrack);
564 if ( (pt->GetNumberOfClusters()>15)) {
565 Int_t found,foundable,shared;
566 pt->GetClusterStatistic(138,158,found, foundable,shared,kFALSE);
567 if (found<15) continue;
568 if (foundable<=0) continue;
569 if (pt->GetNShared()/float(pt->GetNumberOfClusters())>0.2) continue;
570 if (float(found)/float(foundable)<0.8) continue;
573 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
574 iotrack.SetTPCPoints(pt->GetPoints());
575 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
576 iotrack.SetV0Indexes(pt->GetV0Indexes());
577 // iotrack.SetTPCpid(pt->fTPCr);
578 //iotrack.SetTPCindex(i);
579 fEvent->AddTrack(&iotrack);
584 printf("Number of filled ESDs-\t%d\n",fEvent->GetNumberOfTracks());
591 Double_t AliTPCtrackerMI::ErrY2(AliTPCseed* seed, const AliTPCclusterMI * cl){
594 // Use calibrated cluster error from OCDB
596 AliTPCClusterParam * clparam = AliTPCcalibDB::Instance()->GetClusterParam();
598 Float_t z = TMath::Abs(fkParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
599 Int_t ctype = cl->GetType();
600 Int_t type = (cl->GetRow()<63) ? 0: (cl->GetRow()>126) ? 1:2;
601 Double_t angle = seed->GetSnp()*seed->GetSnp();
602 angle = TMath::Sqrt(TMath::Abs(angle/(1.-angle)));
603 Double_t erry2 = clparam->GetError0Par(0,type, z,angle);
605 erry2+=0.5; // edge cluster
608 seed->SetErrorY2(erry2);
612 //calculate look-up table at the beginning
613 // static Bool_t ginit = kFALSE;
614 // static Float_t gnoise1,gnoise2,gnoise3;
615 // static Float_t ggg1[10000];
616 // static Float_t ggg2[10000];
617 // static Float_t ggg3[10000];
618 // static Float_t glandau1[10000];
619 // static Float_t glandau2[10000];
620 // static Float_t glandau3[10000];
622 // static Float_t gcor01[500];
623 // static Float_t gcor02[500];
624 // static Float_t gcorp[500];
628 // if (ginit==kFALSE){
629 // for (Int_t i=1;i<500;i++){
630 // Float_t rsigma = float(i)/100.;
631 // gcor02[i] = TMath::Max(0.78 +TMath::Exp(7.4*(rsigma-1.2)),0.6);
632 // gcor01[i] = TMath::Max(0.72 +TMath::Exp(3.36*(rsigma-1.2)),0.6);
633 // gcorp[i] = TMath::Max(TMath::Power((rsigma+0.5),1.5),1.2);
637 // for (Int_t i=3;i<10000;i++){
641 // Float_t amp = float(i);
642 // Float_t padlength =0.75;
643 // gnoise1 = 0.0004/padlength;
644 // Float_t nel = 0.268*amp;
645 // Float_t nprim = 0.155*amp;
646 // ggg1[i] = fkParam->GetDiffT()*fkParam->GetDiffT()*(2+0.001*nel/(padlength*padlength))/nel;
647 // glandau1[i] = (2.+0.12*nprim)*0.5* (2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
648 // if (glandau1[i]>1) glandau1[i]=1;
649 // glandau1[i]*=padlength*padlength/12.;
653 // gnoise2 = 0.0004/padlength;
655 // nprim = 0.133*amp;
656 // ggg2[i] = fkParam->GetDiffT()*fkParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
657 // glandau2[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
658 // if (glandau2[i]>1) glandau2[i]=1;
659 // glandau2[i]*=padlength*padlength/12.;
664 // gnoise3 = 0.0004/padlength;
666 // nprim = 0.133*amp;
667 // ggg3[i] = fkParam->GetDiffT()*fkParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
668 // glandau3[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
669 // if (glandau3[i]>1) glandau3[i]=1;
670 // glandau3[i]*=padlength*padlength/12.;
678 // Int_t amp = int(TMath::Abs(cl->GetQ()));
680 // seed->SetErrorY2(1.);
684 // Float_t z = TMath::Abs(fkParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
685 // Int_t ctype = cl->GetType();
686 // Float_t padlength= GetPadPitchLength(seed->GetRow());
687 // Double_t angle2 = seed->GetSnp()*seed->GetSnp();
688 // angle2 = angle2/(1-angle2);
690 // //cluster "quality"
691 // Int_t rsigmay = int(100.*cl->GetSigmaY2()/(seed->GetCurrentSigmaY2()));
694 // if (fSectors==fInnerSec){
695 // snoise2 = gnoise1;
696 // res = ggg1[amp]*z+glandau1[amp]*angle2;
697 // if (ctype==0) res *= gcor01[rsigmay];
700 // res*= gcorp[rsigmay];
704 // if (padlength<1.1){
705 // snoise2 = gnoise2;
706 // res = ggg2[amp]*z+glandau2[amp]*angle2;
707 // if (ctype==0) res *= gcor02[rsigmay];
710 // res*= gcorp[rsigmay];
714 // snoise2 = gnoise3;
715 // res = ggg3[amp]*z+glandau3[amp]*angle2;
716 // if (ctype==0) res *= gcor02[rsigmay];
719 // res*= gcorp[rsigmay];
726 // res*=2.4; // overestimate error 2 times
730 // if (res<2*snoise2)
733 // seed->SetErrorY2(res);
741 Double_t AliTPCtrackerMI::ErrZ2(AliTPCseed* seed, const AliTPCclusterMI * cl){
744 // Use calibrated cluster error from OCDB
746 AliTPCClusterParam * clparam = AliTPCcalibDB::Instance()->GetClusterParam();
748 Float_t z = TMath::Abs(fkParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
749 Int_t ctype = cl->GetType();
750 Int_t type = (cl->GetRow()<63) ? 0: (cl->GetRow()>126) ? 1:2;
752 Double_t angle2 = seed->GetSnp()*seed->GetSnp();
753 angle2 = seed->GetTgl()*seed->GetTgl()*(1+angle2/(1-angle2));
754 Double_t angle = TMath::Sqrt(TMath::Abs(angle2));
755 Double_t errz2 = clparam->GetError0Par(1,type, z,angle);
757 errz2+=0.5; // edge cluster
760 seed->SetErrorZ2(errz2);
766 // //seed->SetErrorY2(0.1);
768 // //calculate look-up table at the beginning
769 // static Bool_t ginit = kFALSE;
770 // static Float_t gnoise1,gnoise2,gnoise3;
771 // static Float_t ggg1[10000];
772 // static Float_t ggg2[10000];
773 // static Float_t ggg3[10000];
774 // static Float_t glandau1[10000];
775 // static Float_t glandau2[10000];
776 // static Float_t glandau3[10000];
778 // static Float_t gcor01[1000];
779 // static Float_t gcor02[1000];
780 // static Float_t gcorp[1000];
784 // if (ginit==kFALSE){
785 // for (Int_t i=1;i<1000;i++){
786 // Float_t rsigma = float(i)/100.;
787 // gcor02[i] = TMath::Max(0.81 +TMath::Exp(6.8*(rsigma-1.2)),0.6);
788 // gcor01[i] = TMath::Max(0.72 +TMath::Exp(2.04*(rsigma-1.2)),0.6);
789 // gcorp[i] = TMath::Max(TMath::Power((rsigma+0.5),1.5),1.2);
793 // for (Int_t i=3;i<10000;i++){
797 // Float_t amp = float(i);
798 // Float_t padlength =0.75;
799 // gnoise1 = 0.0004/padlength;
800 // Float_t nel = 0.268*amp;
801 // Float_t nprim = 0.155*amp;
802 // ggg1[i] = fkParam->GetDiffT()*fkParam->GetDiffT()*(2+0.001*nel/(padlength*padlength))/nel;
803 // glandau1[i] = (2.+0.12*nprim)*0.5* (2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
804 // if (glandau1[i]>1) glandau1[i]=1;
805 // glandau1[i]*=padlength*padlength/12.;
809 // gnoise2 = 0.0004/padlength;
811 // nprim = 0.133*amp;
812 // ggg2[i] = fkParam->GetDiffT()*fkParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
813 // glandau2[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
814 // if (glandau2[i]>1) glandau2[i]=1;
815 // glandau2[i]*=padlength*padlength/12.;
820 // gnoise3 = 0.0004/padlength;
822 // nprim = 0.133*amp;
823 // ggg3[i] = fkParam->GetDiffT()*fkParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
824 // glandau3[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
825 // if (glandau3[i]>1) glandau3[i]=1;
826 // glandau3[i]*=padlength*padlength/12.;
834 // Int_t amp = int(TMath::Abs(cl->GetQ()));
836 // seed->SetErrorY2(1.);
840 // Float_t z = TMath::Abs(fkParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
841 // Int_t ctype = cl->GetType();
842 // Float_t padlength= GetPadPitchLength(seed->GetRow());
844 // Double_t angle2 = seed->GetSnp()*seed->GetSnp();
845 // // if (angle2<0.6) angle2 = 0.6;
846 // angle2 = seed->GetTgl()*seed->GetTgl()*(1+angle2/(1-angle2));
848 // //cluster "quality"
849 // Int_t rsigmaz = int(100.*cl->GetSigmaZ2()/(seed->GetCurrentSigmaZ2()));
852 // if (fSectors==fInnerSec){
853 // snoise2 = gnoise1;
854 // res = ggg1[amp]*z+glandau1[amp]*angle2;
855 // if (ctype==0) res *= gcor01[rsigmaz];
858 // res*= gcorp[rsigmaz];
862 // if (padlength<1.1){
863 // snoise2 = gnoise2;
864 // res = ggg2[amp]*z+glandau2[amp]*angle2;
865 // if (ctype==0) res *= gcor02[rsigmaz];
868 // res*= gcorp[rsigmaz];
872 // snoise2 = gnoise3;
873 // res = ggg3[amp]*z+glandau3[amp]*angle2;
874 // if (ctype==0) res *= gcor02[rsigmaz];
877 // res*= gcorp[rsigmaz];
886 // if ((ctype<0) &&<70){
891 // if (res<2*snoise2)
893 // if (res>3) res =3;
894 // seed->SetErrorZ2(res);
902 void AliTPCtrackerMI::RotateToLocal(AliTPCseed *seed)
904 //rotate to track "local coordinata
905 Float_t x = seed->GetX();
906 Float_t y = seed->GetY();
907 Float_t ymax = x*TMath::Tan(0.5*fSectors->GetAlpha());
910 seed->SetRelativeSector((seed->GetRelativeSector()+1) % fN);
911 if (!seed->Rotate(fSectors->GetAlpha()))
913 } else if (y <-ymax) {
914 seed->SetRelativeSector((seed->GetRelativeSector()-1+fN) % fN);
915 if (!seed->Rotate(-fSectors->GetAlpha()))
923 //_____________________________________________________________________________
924 Double_t AliTPCtrackerMI::F1old(Double_t x1,Double_t y1,
925 Double_t x2,Double_t y2,
926 Double_t x3,Double_t y3) const
928 //-----------------------------------------------------------------
929 // Initial approximation of the track curvature
930 //-----------------------------------------------------------------
931 Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
932 Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
933 (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
934 Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
935 (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
937 Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
938 if ( xr*xr+yr*yr<=0.00000000000001) return 100;
939 return -xr*yr/sqrt(xr*xr+yr*yr);
944 //_____________________________________________________________________________
945 Double_t AliTPCtrackerMI::F1(Double_t x1,Double_t y1,
946 Double_t x2,Double_t y2,
947 Double_t x3,Double_t y3) const
949 //-----------------------------------------------------------------
950 // Initial approximation of the track curvature
951 //-----------------------------------------------------------------
957 Double_t det = x3*y2-x2*y3;
962 Double_t u = 0.5* (x2*(x2-x3)+y2*(y2-y3))/det;
963 Double_t x0 = x3*0.5-y3*u;
964 Double_t y0 = y3*0.5+x3*u;
965 Double_t c2 = 1/TMath::Sqrt(x0*x0+y0*y0);
971 Double_t AliTPCtrackerMI::F2(Double_t x1,Double_t y1,
972 Double_t x2,Double_t y2,
973 Double_t x3,Double_t y3) const
975 //-----------------------------------------------------------------
976 // Initial approximation of the track curvature
977 //-----------------------------------------------------------------
983 Double_t det = x3*y2-x2*y3;
988 Double_t u = 0.5* (x2*(x2-x3)+y2*(y2-y3))/det;
989 Double_t x0 = x3*0.5-y3*u;
990 Double_t y0 = y3*0.5+x3*u;
991 Double_t c2 = 1/TMath::Sqrt(x0*x0+y0*y0);
1000 //_____________________________________________________________________________
1001 Double_t AliTPCtrackerMI::F2old(Double_t x1,Double_t y1,
1002 Double_t x2,Double_t y2,
1003 Double_t x3,Double_t y3) const
1005 //-----------------------------------------------------------------
1006 // Initial approximation of the track curvature times center of curvature
1007 //-----------------------------------------------------------------
1008 Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
1009 Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
1010 (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
1011 Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
1012 (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
1014 Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
1016 return -a/(d*y1-b)*xr/sqrt(xr*xr+yr*yr);
1019 //_____________________________________________________________________________
1020 Double_t AliTPCtrackerMI::F3(Double_t x1,Double_t y1,
1021 Double_t x2,Double_t y2,
1022 Double_t z1,Double_t z2) const
1024 //-----------------------------------------------------------------
1025 // Initial approximation of the tangent of the track dip angle
1026 //-----------------------------------------------------------------
1027 return (z1 - z2)/sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
1031 Double_t AliTPCtrackerMI::F3n(Double_t x1,Double_t y1,
1032 Double_t x2,Double_t y2,
1033 Double_t z1,Double_t z2, Double_t c) const
1035 //-----------------------------------------------------------------
1036 // Initial approximation of the tangent of the track dip angle
1037 //-----------------------------------------------------------------
1041 //angle1 = (z1-z2)*c/(TMath::ASin(c*x1-ni)-TMath::ASin(c*x2-ni));
1043 Double_t d = TMath::Sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
1044 if (TMath::Abs(d*c*0.5)>1) return 0;
1045 // Double_t angle2 = TMath::ASin(d*c*0.5);
1046 // Double_t angle2 = AliTPCFastMath::FastAsin(d*c*0.5);
1047 Double_t angle2 = (d*c*0.5>0.1)? TMath::ASin(d*c*0.5): AliTPCFastMath::FastAsin(d*c*0.5);
1049 angle2 = (z1-z2)*c/(angle2*2.);
1053 Bool_t AliTPCtrackerMI::GetProlongation(Double_t x1, Double_t x2, Double_t x[5], Double_t &y, Double_t &z)
1054 {//-----------------------------------------------------------------
1055 // This function find proloncation of a track to a reference plane x=x2.
1056 //-----------------------------------------------------------------
1060 if (TMath::Abs(x[4]*x1 - x[2]) >= 0.999) {
1064 Double_t c1=x[4]*x1 - x[2], r1=sqrt((1.-c1)*(1.+c1));
1065 Double_t c2=x[4]*x2 - x[2], r2=sqrt((1.-c2)*(1.+c2));
1069 Double_t dy = dx*(c1+c2)/(r1+r2);
1072 Double_t delta = x[4]*dx*(c1+c2)/(c1*r2 + c2*r1);
1074 if (TMath::Abs(delta)>0.01){
1075 dz = x[3]*TMath::ASin(delta)/x[4];
1077 dz = x[3]*AliTPCFastMath::FastAsin(delta)/x[4];
1080 //dz = x[3]*AliTPCFastMath::FastAsin(delta)/x[4];
1088 Int_t AliTPCtrackerMI::LoadClusters (TTree *tree)
1093 return LoadClusters();
1097 Int_t AliTPCtrackerMI::LoadClusters(const TObjArray *arr)
1100 // load clusters to the memory
1101 AliTPCClustersRow *clrow = 0x0;
1102 Int_t lower = arr->LowerBound();
1103 Int_t entries = arr->GetEntriesFast();
1105 for (Int_t i=lower; i<entries; i++) {
1106 clrow = (AliTPCClustersRow*) arr->At(i);
1107 if(!clrow) continue;
1108 if(!clrow->GetArray()) continue;
1112 fkParam->AdjustSectorRow(clrow->GetID(),sec,row);
1114 for (Int_t icl=0; icl<clrow->GetArray()->GetEntriesFast(); icl++){
1115 Transform((AliTPCclusterMI*)(clrow->GetArray()->At(icl)));
1118 if (clrow->GetArray()->GetEntriesFast()<=0) continue;
1119 AliTPCtrackerRow * tpcrow=0;
1122 tpcrow = &(fInnerSec[sec%fkNIS][row]);
1126 tpcrow = &(fOuterSec[(sec-fkNIS*2)%fkNOS][row]);
1127 left = (sec-fkNIS*2)/fkNOS;
1130 tpcrow->SetN1(clrow->GetArray()->GetEntriesFast());
1131 tpcrow->SetClusters1(new AliTPCclusterMI[tpcrow->GetN1()]);
1132 for (Int_t j=0;j<tpcrow->GetN1();++j)
1133 tpcrow->SetCluster1(j, *(AliTPCclusterMI*)(clrow->GetArray()->At(j)));
1136 tpcrow->SetN2(clrow->GetArray()->GetEntriesFast());
1137 tpcrow->SetClusters2(new AliTPCclusterMI[tpcrow->GetN2()]);
1138 for (Int_t j=0;j<tpcrow->GetN2();++j)
1139 tpcrow->SetCluster2(j,*(AliTPCclusterMI*)(clrow->GetArray()->At(j)));
1149 Int_t AliTPCtrackerMI::LoadClusters(const TClonesArray *arr)
1152 // load clusters to the memory from one
1155 AliTPCclusterMI *clust=0;
1156 Int_t count[72][96] = { {0} , {0} };
1158 // loop over clusters
1159 for (Int_t icl=0; icl<arr->GetEntriesFast(); icl++) {
1160 clust = (AliTPCclusterMI*)arr->At(icl);
1161 if(!clust) continue;
1162 //printf("cluster: det %d, row %d \n", clust->GetDetector(),clust->GetRow());
1164 // transform clusters
1167 // count clusters per pad row
1168 count[clust->GetDetector()][clust->GetRow()]++;
1171 // insert clusters to sectors
1172 for (Int_t icl=0; icl<arr->GetEntriesFast(); icl++) {
1173 clust = (AliTPCclusterMI*)arr->At(icl);
1174 if(!clust) continue;
1176 Int_t sec = clust->GetDetector();
1177 Int_t row = clust->GetRow();
1179 // filter overlapping pad rows needed by HLT
1180 if(sec<fkNIS*2) { //IROCs
1181 if(row == 30) continue;
1184 if(row == 27 || row == 76) continue;
1190 fInnerSec[sec%fkNIS].InsertCluster(clust, count[sec][row], fkParam);
1193 left = (sec-fkNIS*2)/fkNOS;
1194 fOuterSec[(sec-fkNIS*2)%fkNOS].InsertCluster(clust, count[sec][row], fkParam);
1198 // Load functions must be called behind LoadCluster(TClonesArray*)
1200 //LoadOuterSectors();
1201 //LoadInnerSectors();
1207 Int_t AliTPCtrackerMI::LoadClusters()
1210 // load clusters to the memory
1211 AliTPCClustersRow *clrow= new AliTPCClustersRow;
1212 clrow->SetClass("AliTPCclusterMI");
1214 clrow->GetArray()->ExpandCreateFast(10000);
1216 // TTree * tree = fClustersArray.GetTree();
1218 TTree * tree = fInput;
1219 TBranch * br = tree->GetBranch("Segment");
1220 br->SetAddress(&clrow);
1222 Int_t j=Int_t(tree->GetEntries());
1223 for (Int_t i=0; i<j; i++) {
1227 fkParam->AdjustSectorRow(clrow->GetID(),sec,row);
1228 for (Int_t icl=0; icl<clrow->GetArray()->GetEntriesFast(); icl++){
1229 Transform((AliTPCclusterMI*)(clrow->GetArray()->At(icl)));
1232 AliTPCtrackerRow * tpcrow=0;
1235 tpcrow = &(fInnerSec[sec%fkNIS][row]);
1239 tpcrow = &(fOuterSec[(sec-fkNIS*2)%fkNOS][row]);
1240 left = (sec-fkNIS*2)/fkNOS;
1243 tpcrow->SetN1(clrow->GetArray()->GetEntriesFast());
1244 tpcrow->SetClusters1(new AliTPCclusterMI[tpcrow->GetN1()]);
1245 for (Int_t k=0;k<tpcrow->GetN1();++k)
1246 tpcrow->SetCluster1(k, *(AliTPCclusterMI*)(clrow->GetArray()->At(k)));
1249 tpcrow->SetN2(clrow->GetArray()->GetEntriesFast());
1250 tpcrow->SetClusters2(new AliTPCclusterMI[tpcrow->GetN2()]);
1251 for (Int_t k=0;k<tpcrow->GetN2();++k)
1252 tpcrow->SetCluster2(k,*(AliTPCclusterMI*)(clrow->GetArray()->At(k)));
1263 void AliTPCtrackerMI::UnloadClusters()
1266 // unload clusters from the memory
1268 Int_t nrows = fOuterSec->GetNRows();
1269 for (Int_t sec = 0;sec<fkNOS;sec++)
1270 for (Int_t row = 0;row<nrows;row++){
1271 AliTPCtrackerRow* tpcrow = &(fOuterSec[sec%fkNOS][row]);
1273 // if (tpcrow->fClusters1) delete []tpcrow->fClusters1;
1274 // if (tpcrow->fClusters2) delete []tpcrow->fClusters2;
1276 tpcrow->ResetClusters();
1279 nrows = fInnerSec->GetNRows();
1280 for (Int_t sec = 0;sec<fkNIS;sec++)
1281 for (Int_t row = 0;row<nrows;row++){
1282 AliTPCtrackerRow* tpcrow = &(fInnerSec[sec%fkNIS][row]);
1284 // if (tpcrow->fClusters1) delete []tpcrow->fClusters1;
1285 //if (tpcrow->fClusters2) delete []tpcrow->fClusters2;
1287 tpcrow->ResetClusters();
1293 void AliTPCtrackerMI::FillClusterArray(TObjArray* array) const{
1295 // Filling cluster to the array - For visualization purposes
1298 nrows = fOuterSec->GetNRows();
1299 for (Int_t sec = 0;sec<fkNOS;sec++)
1300 for (Int_t row = 0;row<nrows;row++){
1301 AliTPCtrackerRow* tpcrow = &(fOuterSec[sec%fkNOS][row]);
1302 if (!tpcrow) continue;
1303 for (Int_t icl = 0;icl<tpcrow->GetN();icl++){
1304 array->AddLast((TObject*)((*tpcrow)[icl]));
1307 nrows = fInnerSec->GetNRows();
1308 for (Int_t sec = 0;sec<fkNIS;sec++)
1309 for (Int_t row = 0;row<nrows;row++){
1310 AliTPCtrackerRow* tpcrow = &(fInnerSec[sec%fkNIS][row]);
1311 if (!tpcrow) continue;
1312 for (Int_t icl = 0;icl<tpcrow->GetN();icl++){
1313 array->AddLast((TObject*)(*tpcrow)[icl]);
1319 void AliTPCtrackerMI::Transform(AliTPCclusterMI * cluster){
1323 AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ;
1325 AliFatal("Tranformations not in calibDB");
1327 Double_t x[3]={cluster->GetRow(),cluster->GetPad(),cluster->GetTimeBin()};
1328 Int_t i[1]={cluster->GetDetector()};
1329 transform->Transform(x,i,0,1);
1330 // if (cluster->GetDetector()%36>17){
1335 // in debug mode check the transformation
1337 if (AliTPCReconstructor::StreamLevel()>1) {
1339 cluster->GetGlobalXYZ(gx);
1340 Int_t event = (fEvent==NULL)? 0: fEvent->GetEventNumberInFile();
1341 TTreeSRedirector &cstream = *fDebugStreamer;
1342 cstream<<"Transform"<<
1353 cluster->SetX(x[0]);
1354 cluster->SetY(x[1]);
1355 cluster->SetZ(x[2]);
1361 //if (!fkParam->IsGeoRead()) fkParam->ReadGeoMatrices();
1362 TGeoHMatrix *mat = fkParam->GetClusterMatrix(cluster->GetDetector());
1364 Double_t pos[3]= {cluster->GetX(),cluster->GetY(),cluster->GetZ()};
1365 Double_t posC[3]={cluster->GetX(),cluster->GetY(),cluster->GetZ()};
1366 if (mat) mat->LocalToMaster(pos,posC);
1368 // chack Loading of Geo matrices from GeoManager - TEMPORARY FIX
1370 cluster->SetX(posC[0]);
1371 cluster->SetY(posC[1]);
1372 cluster->SetZ(posC[2]);
1375 //_____________________________________________________________________________
1376 Int_t AliTPCtrackerMI::LoadOuterSectors() {
1377 //-----------------------------------------------------------------
1378 // This function fills outer TPC sectors with clusters.
1379 //-----------------------------------------------------------------
1380 Int_t nrows = fOuterSec->GetNRows();
1382 for (Int_t sec = 0;sec<fkNOS;sec++)
1383 for (Int_t row = 0;row<nrows;row++){
1384 AliTPCtrackerRow* tpcrow = &(fOuterSec[sec%fkNOS][row]);
1385 Int_t sec2 = sec+2*fkNIS;
1387 Int_t ncl = tpcrow->GetN1();
1389 AliTPCclusterMI *c= (tpcrow->GetCluster1(ncl));
1390 index=(((sec2<<8)+row)<<16)+ncl;
1391 tpcrow->InsertCluster(c,index);
1394 ncl = tpcrow->GetN2();
1396 AliTPCclusterMI *c= (tpcrow->GetCluster2(ncl));
1397 index=((((sec2+fkNOS)<<8)+row)<<16)+ncl;
1398 tpcrow->InsertCluster(c,index);
1401 // write indexes for fast acces
1403 for (Int_t i=0;i<510;i++)
1404 tpcrow->SetFastCluster(i,-1);
1405 for (Int_t i=0;i<tpcrow->GetN();i++){
1406 Int_t zi = Int_t((*tpcrow)[i]->GetZ()+255.);
1407 tpcrow->SetFastCluster(zi,i); // write index
1410 for (Int_t i=0;i<510;i++){
1411 if (tpcrow->GetFastCluster(i)<0)
1412 tpcrow->SetFastCluster(i,last);
1414 last = tpcrow->GetFastCluster(i);
1423 //_____________________________________________________________________________
1424 Int_t AliTPCtrackerMI::LoadInnerSectors() {
1425 //-----------------------------------------------------------------
1426 // This function fills inner TPC sectors with clusters.
1427 //-----------------------------------------------------------------
1428 Int_t nrows = fInnerSec->GetNRows();
1430 for (Int_t sec = 0;sec<fkNIS;sec++)
1431 for (Int_t row = 0;row<nrows;row++){
1432 AliTPCtrackerRow* tpcrow = &(fInnerSec[sec%fkNIS][row]);
1435 Int_t ncl = tpcrow->GetN1();
1437 AliTPCclusterMI *c= (tpcrow->GetCluster1(ncl));
1438 index=(((sec<<8)+row)<<16)+ncl;
1439 tpcrow->InsertCluster(c,index);
1442 ncl = tpcrow->GetN2();
1444 AliTPCclusterMI *c= (tpcrow->GetCluster2(ncl));
1445 index=((((sec+fkNIS)<<8)+row)<<16)+ncl;
1446 tpcrow->InsertCluster(c,index);
1449 // write indexes for fast acces
1451 for (Int_t i=0;i<510;i++)
1452 tpcrow->SetFastCluster(i,-1);
1453 for (Int_t i=0;i<tpcrow->GetN();i++){
1454 Int_t zi = Int_t((*tpcrow)[i]->GetZ()+255.);
1455 tpcrow->SetFastCluster(zi,i); // write index
1458 for (Int_t i=0;i<510;i++){
1459 if (tpcrow->GetFastCluster(i)<0)
1460 tpcrow->SetFastCluster(i,last);
1462 last = tpcrow->GetFastCluster(i);
1474 //_________________________________________________________________________
1475 AliTPCclusterMI *AliTPCtrackerMI::GetClusterMI(Int_t index) const {
1476 //--------------------------------------------------------------------
1477 // Return pointer to a given cluster
1478 //--------------------------------------------------------------------
1479 if (index<0) return 0; // no cluster
1480 Int_t sec=(index&0xff000000)>>24;
1481 Int_t row=(index&0x00ff0000)>>16;
1482 Int_t ncl=(index&0x00007fff)>>00;
1484 const AliTPCtrackerRow * tpcrow=0;
1485 AliTPCclusterMI * clrow =0;
1487 if (sec<0 || sec>=fkNIS*4) {
1488 AliWarning(Form("Wrong sector %d",sec));
1493 AliTPCtrackerSector& tracksec = fInnerSec[sec%fkNIS];
1494 if (tracksec.GetNRows()<=row) return 0;
1495 tpcrow = &(tracksec[row]);
1496 if (tpcrow==0) return 0;
1499 if (tpcrow->GetN1()<=ncl) return 0;
1500 clrow = tpcrow->GetClusters1();
1503 if (tpcrow->GetN2()<=ncl) return 0;
1504 clrow = tpcrow->GetClusters2();
1508 AliTPCtrackerSector& tracksec = fOuterSec[(sec-fkNIS*2)%fkNOS];
1509 if (tracksec.GetNRows()<=row) return 0;
1510 tpcrow = &(tracksec[row]);
1511 if (tpcrow==0) return 0;
1513 if (sec-2*fkNIS<fkNOS) {
1514 if (tpcrow->GetN1()<=ncl) return 0;
1515 clrow = tpcrow->GetClusters1();
1518 if (tpcrow->GetN2()<=ncl) return 0;
1519 clrow = tpcrow->GetClusters2();
1523 return &(clrow[ncl]);
1529 Int_t AliTPCtrackerMI::FollowToNext(AliTPCseed& t, Int_t nr) {
1530 //-----------------------------------------------------------------
1531 // This function tries to find a track prolongation to next pad row
1532 //-----------------------------------------------------------------
1534 Double_t x= GetXrow(nr), ymax=GetMaxY(nr);
1537 AliTPCclusterMI *cl=0;
1538 Int_t tpcindex= t.GetClusterIndex2(nr);
1540 // update current shape info every 5 pad-row
1541 // if ( (nr%5==0) || t.GetNumberOfClusters()<2 || (t.fCurrentSigmaY2<0.0001) ){
1545 if (fIteration>0 && tpcindex>=-1){ //if we have already clusters
1547 if (tpcindex==-1) return 0; //track in dead zone
1549 cl = t.GetClusterPointer(nr);
1550 if ( (cl==0) ) cl = GetClusterMI(tpcindex);
1551 t.SetCurrentClusterIndex1(tpcindex);
1554 Int_t relativesector = ((tpcindex&0xff000000)>>24)%18; // if previously accepted cluster in different sector
1555 Float_t angle = relativesector*fSectors->GetAlpha()+fSectors->GetAlphaShift();
1557 if (angle<-TMath::Pi()) angle += 2*TMath::Pi();
1558 if (angle>=TMath::Pi()) angle -= 2*TMath::Pi();
1560 if (TMath::Abs(angle-t.GetAlpha())>0.001){
1561 Double_t rotation = angle-t.GetAlpha();
1562 t.SetRelativeSector(relativesector);
1563 if (!t.Rotate(rotation)) return 0;
1565 if (!t.PropagateTo(x)) return 0;
1567 t.SetCurrentCluster(cl);
1569 Int_t accept = AcceptCluster(&t,t.GetCurrentCluster());
1570 if ((tpcindex&0x8000)==0) accept =0;
1572 //if founded cluster is acceptible
1573 if (cl->IsUsed(11)) { // id cluster is shared inrease uncertainty
1574 t.SetErrorY2(t.GetErrorY2()+0.03);
1575 t.SetErrorZ2(t.GetErrorZ2()+0.03);
1576 t.SetErrorY2(t.GetErrorY2()*3);
1577 t.SetErrorZ2(t.GetErrorZ2()*3);
1579 t.SetNFoundable(t.GetNFoundable()+1);
1580 UpdateTrack(&t,accept);
1585 if (TMath::Abs(t.GetSnp())>AliTPCReconstructor::GetMaxSnpTracker()) return 0; // cut on angle
1586 if (fIteration>1 && IsFindable(t)){
1587 // not look for new cluster during refitting
1588 t.SetNFoundable(t.GetNFoundable()+1);
1593 // if (TMath::Abs(t.GetSnp())>0.95 || TMath::Abs(x*t.GetC()-t.GetEta())>0.95) return 0;// patch 28 fev 06
1594 Double_t y=t.GetYat(x);
1595 if (TMath::Abs(y)>ymax){
1597 t.SetRelativeSector((t.GetRelativeSector()+1) % fN);
1598 if (!t.Rotate(fSectors->GetAlpha()))
1600 } else if (y <-ymax) {
1601 t.SetRelativeSector((t.GetRelativeSector()-1+fN) % fN);
1602 if (!t.Rotate(-fSectors->GetAlpha()))
1608 if (!t.PropagateTo(x)) {
1609 if (fIteration==0) t.SetRemoval(10);
1613 Double_t z=t.GetZ();
1616 if (!IsActive(t.GetRelativeSector(),nr)) {
1618 t.SetClusterIndex2(nr,-1);
1621 //AliInfo(Form("A - Sector%d phi %f - alpha %f", t.fRelativeSector,y/x, t.GetAlpha()));
1622 Bool_t isActive = IsActive(t.GetRelativeSector(),nr);
1623 Bool_t isActive2 = (nr>=fInnerSec->GetNRows()) ? fOuterSec[t.GetRelativeSector()][nr-fInnerSec->GetNRows()].GetN()>0:fInnerSec[t.GetRelativeSector()][nr].GetN()>0;
1625 if (!isActive || !isActive2) return 0;
1627 const AliTPCtrackerRow &krow=GetRow(t.GetRelativeSector(),nr);
1628 if ( (t.GetSigmaY2()<0) || t.GetSigmaZ2()<0) return 0;
1630 Double_t roadz = 1.;
1632 if (TMath::Abs(TMath::Abs(y)-ymax)<krow.GetDeadZone()){
1634 t.SetClusterIndex2(nr,-1);
1640 // if (TMath::Abs(z)<(AliTPCReconstructor::GetCtgRange()*x+10) && TMath::Abs(z)<fkParam->GetZLength(0) && (TMath::Abs(t.GetSnp())<AliTPCReconstructor::GetMaxSnpTracker()))
1641 t.SetNFoundable(t.GetNFoundable()+1);
1647 // cl = krow.FindNearest2(y+10.,z,roady,roadz,index);
1648 cl = krow.FindNearest2(y,z,roady,roadz,index);
1649 if (cl) t.SetCurrentClusterIndex1(krow.GetIndex(index));
1652 t.SetCurrentCluster(cl);
1654 if (fIteration==2&&cl->IsUsed(10)) return 0;
1655 Int_t accept = AcceptCluster(&t,t.GetCurrentCluster());
1656 if (fIteration==2&&cl->IsUsed(11)) {
1657 t.SetErrorY2(t.GetErrorY2()+0.03);
1658 t.SetErrorZ2(t.GetErrorZ2()+0.03);
1659 t.SetErrorY2(t.GetErrorY2()*3);
1660 t.SetErrorZ2(t.GetErrorZ2()*3);
1663 if (t.fCurrentCluster->IsUsed(10)){
1668 if (t.fNShared>0.7*t.GetNumberOfClusters()) {
1674 if (accept<3) UpdateTrack(&t,accept);
1677 if ( fIteration==0 && t.GetNFoundable()*0.5 > t.GetNumberOfClusters()) t.SetRemoval(10);
1685 //_________________________________________________________________________
1686 Bool_t AliTPCtrackerMI::GetTrackPoint(Int_t index, AliTrackPoint &p ) const
1688 // Get track space point by index
1689 // return false in case the cluster doesn't exist
1690 AliTPCclusterMI *cl = GetClusterMI(index);
1691 if (!cl) return kFALSE;
1692 Int_t sector = (index&0xff000000)>>24;
1693 // Int_t row = (index&0x00ff0000)>>16;
1695 // xyz[0] = fkParam->GetPadRowRadii(sector,row);
1696 xyz[0] = cl->GetX();
1697 xyz[1] = cl->GetY();
1698 xyz[2] = cl->GetZ();
1700 fkParam->AdjustCosSin(sector,cos,sin);
1701 Float_t x = cos*xyz[0]-sin*xyz[1];
1702 Float_t y = cos*xyz[1]+sin*xyz[0];
1704 Float_t sigmaY2 = 0.027*cl->GetSigmaY2();
1705 if (sector < fkParam->GetNInnerSector()) sigmaY2 *= 2.07;
1706 Float_t sigmaZ2 = 0.066*cl->GetSigmaZ2();
1707 if (sector < fkParam->GetNInnerSector()) sigmaZ2 *= 1.77;
1708 cov[0] = sin*sin*sigmaY2;
1709 cov[1] = -sin*cos*sigmaY2;
1711 cov[3] = cos*cos*sigmaY2;
1714 p.SetXYZ(x,y,xyz[2],cov);
1715 AliGeomManager::ELayerID iLayer;
1717 if (sector < fkParam->GetNInnerSector()) {
1718 iLayer = AliGeomManager::kTPC1;
1722 iLayer = AliGeomManager::kTPC2;
1723 idet = sector - fkParam->GetNInnerSector();
1725 UShort_t volid = AliGeomManager::LayerToVolUID(iLayer,idet);
1726 p.SetVolumeID(volid);
1732 Int_t AliTPCtrackerMI::UpdateClusters(AliTPCseed& t, Int_t nr) {
1733 //-----------------------------------------------------------------
1734 // This function tries to find a track prolongation to next pad row
1735 //-----------------------------------------------------------------
1736 t.SetCurrentCluster(0);
1737 t.SetCurrentClusterIndex1(0);
1739 Double_t xt=t.GetX();
1740 Int_t row = GetRowNumber(xt)-1;
1741 Double_t ymax= GetMaxY(nr);
1743 if (row < nr) return 1; // don't prolongate if not information until now -
1744 // if (TMath::Abs(t.GetSnp())>0.9 && t.GetNumberOfClusters()>40. && fIteration!=2) {
1746 // return 0; // not prolongate strongly inclined tracks
1748 // if (TMath::Abs(t.GetSnp())>0.95) {
1750 // return 0; // not prolongate strongly inclined tracks
1751 // }// patch 28 fev 06
1753 Double_t x= GetXrow(nr);
1755 //t.PropagateTo(x+0.02);
1756 //t.PropagateTo(x+0.01);
1757 if (!t.PropagateTo(x)){
1764 if (TMath::Abs(y)>ymax){
1766 t.SetRelativeSector((t.GetRelativeSector()+1) % fN);
1767 if (!t.Rotate(fSectors->GetAlpha()))
1769 } else if (y <-ymax) {
1770 t.SetRelativeSector((t.GetRelativeSector()-1+fN) % fN);
1771 if (!t.Rotate(-fSectors->GetAlpha()))
1774 // if (!t.PropagateTo(x)){
1781 if (TMath::Abs(t.GetSnp())>AliTPCReconstructor::GetMaxSnpTracker()) return 0;
1783 if (!IsActive(t.GetRelativeSector(),nr)) {
1785 t.SetClusterIndex2(nr,-1);
1788 //AliInfo(Form("A - Sector%d phi %f - alpha %f", t.fRelativeSector,y/x, t.GetAlpha()));
1790 AliTPCtrackerRow &krow=GetRow(t.GetRelativeSector(),nr);
1792 if (TMath::Abs(TMath::Abs(y)-ymax)<krow.GetDeadZone()){
1794 t.SetClusterIndex2(nr,-1);
1800 // if (TMath::Abs(t.GetZ())<(AliTPCReconstructor::GetCtgRange()*t.GetX()+10) && (TMath::Abs(t.GetSnp())<AliTPCReconstructor::GetMaxSnpTracker()))
1801 if (IsFindable(t)) t.SetNFoundable(t.GetNFoundable()+1);
1807 if ( (nr%2==0) || t.GetNumberOfClusters()<2){
1808 // t.fCurrentSigmaY = GetSigmaY(&t);
1809 //t.fCurrentSigmaZ = GetSigmaZ(&t);
1813 AliTPCclusterMI *cl=0;
1816 Double_t roady = 1.;
1817 Double_t roadz = 1.;
1821 index = t.GetClusterIndex2(nr);
1822 if ( (index>0) && (index&0x8000)==0){
1823 cl = t.GetClusterPointer(nr);
1824 if ( (cl==0) && (index>0)) cl = GetClusterMI(index);
1825 t.SetCurrentClusterIndex1(index);
1827 t.SetCurrentCluster(cl);
1833 // if (index<0) return 0;
1834 UInt_t uindex = TMath::Abs(index);
1837 //cl = krow.FindNearest2(y+10,z,roady,roadz,uindex);
1838 cl = krow.FindNearest2(y,z,roady,roadz,uindex);
1841 if (cl) t.SetCurrentClusterIndex1(krow.GetIndex(uindex));
1842 t.SetCurrentCluster(cl);
1848 Int_t AliTPCtrackerMI::FollowToNextCluster(AliTPCseed & t, Int_t nr) {
1849 //-----------------------------------------------------------------
1850 // This function tries to find a track prolongation to next pad row
1851 //-----------------------------------------------------------------
1853 //update error according neighborhoud
1855 if (t.GetCurrentCluster()) {
1857 Int_t accept = AcceptCluster(&t,t.GetCurrentCluster());
1859 if (t.GetCurrentCluster()->IsUsed(10)){
1864 t.SetNShared(t.GetNShared()+1);
1865 if (t.GetNShared()>0.7*t.GetNumberOfClusters()) {
1870 if (fIteration>0) accept = 0;
1871 if (accept<3) UpdateTrack(&t,accept);
1875 if ( ( (t.GetSigmaY2()+t.GetSigmaZ2())>0.16)&& t.GetNumberOfClusters()>18) t.SetRemoval(10);
1876 if ( t.GetChi2()/t.GetNumberOfClusters()>6 &&t.GetNumberOfClusters()>18) t.SetRemoval(10);
1878 if (( (t.GetNFoundable()*0.5 > t.GetNumberOfClusters()) || t.GetNoCluster()>15)) t.SetRemoval(10);
1886 //_____________________________________________________________________________
1887 Int_t AliTPCtrackerMI::FollowProlongation(AliTPCseed& t, Int_t rf, Int_t step) {
1888 //-----------------------------------------------------------------
1889 // This function tries to find a track prolongation.
1890 //-----------------------------------------------------------------
1891 Double_t xt=t.GetX();
1893 Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
1894 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
1895 if (alpha < 0. ) alpha += 2.*TMath::Pi();
1897 t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN);
1899 Int_t first = GetRowNumber(xt)-1;
1900 for (Int_t nr= first; nr>=rf; nr-=step) {
1902 if (t.GetKinkIndexes()[0]>0){
1903 for (Int_t i=0;i<3;i++){
1904 Int_t index = t.GetKinkIndexes()[i];
1905 if (index==0) break;
1906 if (index<0) continue;
1908 AliKink * kink = (AliKink*)fEvent->GetKink(index-1);
1910 printf("PROBLEM\n");
1913 Int_t kinkrow = kink->GetTPCRow0()+2+Int_t(0.5/(0.05+kink->GetAngle(2)));
1915 AliExternalTrackParam paramd(t);
1916 kink->SetDaughter(paramd);
1917 kink->SetStatus(2,5);
1924 if (nr==80) t.UpdateReference();
1925 if (nr<fInnerSec->GetNRows())
1926 fSectors = fInnerSec;
1928 fSectors = fOuterSec;
1929 if (FollowToNext(t,nr)==0)
1942 Int_t AliTPCtrackerMI::FollowBackProlongation(AliTPCseed& t, Int_t rf) {
1943 //-----------------------------------------------------------------
1944 // This function tries to find a track prolongation.
1945 //-----------------------------------------------------------------
1947 Double_t xt=t.GetX();
1948 Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
1949 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
1950 if (alpha < 0. ) alpha += 2.*TMath::Pi();
1951 t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN);
1953 Int_t first = t.GetFirstPoint();
1954 if (first<GetRowNumber(xt)+1) first = GetRowNumber(xt)+1;
1956 if (first<0) first=0;
1957 for (Int_t nr=first; nr<=rf; nr++) {
1958 // if ( (TMath::Abs(t.GetSnp())>0.95)) break;//patch 28 fev 06
1959 if (t.GetKinkIndexes()[0]<0){
1960 for (Int_t i=0;i<3;i++){
1961 Int_t index = t.GetKinkIndexes()[i];
1962 if (index==0) break;
1963 if (index>0) continue;
1964 index = TMath::Abs(index);
1965 AliKink * kink = (AliKink*)fEvent->GetKink(index-1);
1967 printf("PROBLEM\n");
1970 Int_t kinkrow = kink->GetTPCRow0()-2-Int_t(0.5/(0.05+kink->GetAngle(2)));
1972 AliExternalTrackParam paramm(t);
1973 kink->SetMother(paramm);
1974 kink->SetStatus(2,1);
1981 if (nr<fInnerSec->GetNRows())
1982 fSectors = fInnerSec;
1984 fSectors = fOuterSec;
1995 Float_t AliTPCtrackerMI::OverlapFactor(AliTPCseed * s1, AliTPCseed * s2, Int_t &sum1, Int_t & sum2)
2003 Float_t dz2 =(s1->GetZ() - s2->GetZ());
2006 Float_t dy2 =TMath::Abs((s1->GetY() - s2->GetY()));
2008 Float_t distance = TMath::Sqrt(dz2+dy2);
2009 if (distance>4.) return 0; // if there are far away - not overlap - to reduce combinatorics
2012 Int_t firstpoint = TMath::Min(s1->GetFirstPoint(),s2->GetFirstPoint());
2013 Int_t lastpoint = TMath::Max(s1->GetLastPoint(),s2->GetLastPoint());
2018 if (firstpoint>lastpoint) {
2019 firstpoint =lastpoint;
2024 for (Int_t i=firstpoint-1;i<lastpoint+1;i++){
2025 if (s1->GetClusterIndex2(i)>0) sum1++;
2026 if (s2->GetClusterIndex2(i)>0) sum2++;
2027 if (s1->GetClusterIndex2(i)==s2->GetClusterIndex2(i) && s1->GetClusterIndex2(i)>0) {
2031 if (sum<5) return 0;
2033 Float_t summin = TMath::Min(sum1+1,sum2+1);
2034 Float_t ratio = (sum+1)/Float_t(summin);
2038 void AliTPCtrackerMI::SignShared(AliTPCseed * s1, AliTPCseed * s2)
2042 Float_t thetaCut = 0.2;//+10.*TMath::Sqrt(s1->GetSigmaTglZ()+ s2->GetSigmaTglZ());
2043 if (TMath::Abs(s1->GetTgl()-s2->GetTgl())>thetaCut) return;
2044 Float_t minCl = TMath::Min(s1->GetNumberOfClusters(),s2->GetNumberOfClusters());
2045 Int_t cutN0 = TMath::Max(5,TMath::Nint(0.1*minCl));
2050 //Int_t firstpoint = TMath::Max(s1->GetFirstPoint(),s2->GetFirstPoint());
2051 //Int_t lastpoint = TMath::Min(s1->GetLastPoint(),s2->GetLastPoint());
2052 Int_t firstpoint = 0;
2053 Int_t lastpoint = 160;
2055 // if (firstpoint>=lastpoint-5) return;;
2057 for (Int_t i=firstpoint;i<lastpoint;i++){
2058 // if ( (s1->GetClusterIndex2(i)&0xFFFF8FFF)==(s2->GetClusterIndex2(i)&0xFFFF8FFF) && s1->GetClusterIndex2(i)>0) {
2059 if ( (s1->GetClusterIndex2(i))==(s2->GetClusterIndex2(i)) && s1->GetClusterIndex2(i)>0) {
2061 s1->SetSharedMapBit(i, kTRUE);
2062 s2->SetSharedMapBit(i, kTRUE);
2064 if (s1->GetClusterIndex2(i)>0)
2065 s1->SetClusterMapBit(i, kTRUE);
2067 if (sumshared>cutN0){
2070 for (Int_t i=firstpoint;i<lastpoint;i++){
2071 // if ( (s1->GetClusterIndex2(i)&0xFFFF8FFF)==(s2->GetClusterIndex2(i)&0xFFFF8FFF) && s1->GetClusterIndex2(i)>0) {
2072 if ( (s1->GetClusterIndex2(i))==(s2->GetClusterIndex2(i)) && s1->GetClusterIndex2(i)>0) {
2073 AliTPCTrackerPoint *p1 = s1->GetTrackPoint(i);
2074 AliTPCTrackerPoint *p2 = s2->GetTrackPoint(i);;
2075 if (s1->IsActive()&&s2->IsActive()){
2076 p1->SetShared(kTRUE);
2077 p2->SetShared(kTRUE);
2083 if (sumshared>cutN0){
2084 for (Int_t i=0;i<4;i++){
2085 if (s1->GetOverlapLabel(3*i)==0){
2086 s1->SetOverlapLabel(3*i, s2->GetLabel());
2087 s1->SetOverlapLabel(3*i+1,sumshared);
2088 s1->SetOverlapLabel(3*i+2,s2->GetUniqueID());
2092 for (Int_t i=0;i<4;i++){
2093 if (s2->GetOverlapLabel(3*i)==0){
2094 s2->SetOverlapLabel(3*i, s1->GetLabel());
2095 s2->SetOverlapLabel(3*i+1,sumshared);
2096 s2->SetOverlapLabel(3*i+2,s1->GetUniqueID());
2103 void AliTPCtrackerMI::SignShared(TObjArray * arr)
2106 //sort trackss according sectors
2108 for (Int_t i=0; i<arr->GetEntriesFast(); i++) {
2109 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2111 //if (pt) RotateToLocal(pt);
2115 arr->Sort(); // sorting according relative sectors
2116 arr->Expand(arr->GetEntries());
2119 Int_t nseed=arr->GetEntriesFast();
2120 for (Int_t i=0; i<nseed; i++) {
2121 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2123 for (Int_t j=0;j<=12;j++){
2124 pt->SetOverlapLabel(j,0);
2127 for (Int_t i=0; i<nseed; i++) {
2128 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2130 if (pt->GetRemoval()>10) continue;
2131 for (Int_t j=i+1; j<nseed; j++){
2132 AliTPCseed *pt2=(AliTPCseed*)arr->UncheckedAt(j);
2133 if (TMath::Abs(pt->GetRelativeSector()-pt2->GetRelativeSector())>1) continue;
2135 if (pt2->GetRemoval()<=10) {
2136 //if ( TMath::Abs(pt->GetRelativeSector()-pt2->GetRelativeSector())>0) break;
2144 void AliTPCtrackerMI::SortTracks(TObjArray * arr, Int_t mode) const
2147 //sort tracks in array according mode criteria
2148 Int_t nseed = arr->GetEntriesFast();
2149 for (Int_t i=0; i<nseed; i++) {
2150 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2161 void AliTPCtrackerMI::RemoveUsed2(TObjArray * arr, Float_t factor1, Float_t factor2, Int_t minimal)
2164 // Loop over all tracks and remove overlaped tracks (with lower quality)
2166 // 1. Unsign clusters
2167 // 2. Sort tracks according quality
2168 // Quality is defined by the number of cluster between first and last points
2170 // 3. Loop over tracks - decreasing quality order
2171 // a.) remove - If the fraction of shared cluster less than factor (1- n or 2)
2172 // b.) remove - If the minimal number of clusters less than minimal and not ITS
2173 // c.) if track accepted - sign clusters
2175 //Called in - AliTPCtrackerMI::Clusters2Tracks()
2176 // - AliTPCtrackerMI::PropagateBack()
2177 // - AliTPCtrackerMI::RefitInward()
2180 // factor1 - factor for constrained
2181 // factor2 - for non constrained tracks
2182 // if (Float_t(shared+1)/Float_t(found+1)>factor) - DELETE
2186 Int_t nseed = arr->GetEntriesFast();
2187 Float_t * quality = new Float_t[nseed];
2188 Int_t * indexes = new Int_t[nseed];
2192 for (Int_t i=0; i<nseed; i++) {
2193 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2198 pt->UpdatePoints(); //select first last max dens points
2199 Float_t * points = pt->GetPoints();
2200 if (points[3]<0.8) quality[i] =-1;
2201 quality[i] = (points[2]-points[0])+pt->GetNumberOfClusters();
2202 //prefer high momenta tracks if overlaps
2203 quality[i] *= TMath::Sqrt(TMath::Abs(pt->Pt())+0.5);
2205 TMath::Sort(nseed,quality,indexes);
2208 for (Int_t itrack=0; itrack<nseed; itrack++) {
2209 Int_t trackindex = indexes[itrack];
2210 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(trackindex);
2213 if (quality[trackindex]<0){
2215 delete arr->RemoveAt(trackindex);
2218 arr->RemoveAt(trackindex);
2224 Int_t first = Int_t(pt->GetPoints()[0]);
2225 Int_t last = Int_t(pt->GetPoints()[2]);
2226 Double_t factor = (pt->GetBConstrain()) ? factor1: factor2;
2228 Int_t found,foundable,shared;
2229 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
2230 // pt->GetClusterStatistic(0,160, found, foundable,shared,kFALSE);
2231 Bool_t itsgold =kFALSE;
2234 if (pt->GetESD()->GetITSclusters(dummy)>4) itsgold= kTRUE;
2238 if (Float_t(shared+1)/Float_t(found+1)>factor){
2239 if (pt->GetKinkIndexes()[0]!=0) continue; //don't remove tracks - part of the kinks
2240 if( AliTPCReconstructor::StreamLevel()>15){
2241 TTreeSRedirector &cstream = *fDebugStreamer;
2242 cstream<<"RemoveUsed"<<
2243 "iter="<<fIteration<<
2247 delete arr->RemoveAt(trackindex);
2250 if (pt->GetNumberOfClusters()<50&&(found-0.5*shared)<minimal){ //remove short tracks
2251 if (pt->GetKinkIndexes()[0]!=0) continue; //don't remove tracks - part of the kinks
2252 if( AliTPCReconstructor::StreamLevel()>15){
2253 TTreeSRedirector &cstream = *fDebugStreamer;
2254 cstream<<"RemoveShort"<<
2255 "iter="<<fIteration<<
2259 delete arr->RemoveAt(trackindex);
2265 //if (sharedfactor>0.4) continue;
2266 if (pt->GetKinkIndexes()[0]>0) continue;
2267 //Remove tracks with undefined properties - seems
2268 if (pt->GetSigmaY2()<kAlmost0) continue; // ? what is the origin ?
2270 for (Int_t i=first; i<last; i++) {
2271 Int_t index=pt->GetClusterIndex2(i);
2272 // if (index<0 || index&0x8000 ) continue;
2273 if (index<0 || index&0x8000 ) continue;
2274 AliTPCclusterMI *c= pt->GetClusterPointer(i);
2281 Info("RemoveUsed","\n*****\nNumber of good tracks after shared removal\t%d\n",fNtracks);
2287 void AliTPCtrackerMI::UnsignClusters()
2290 // loop over all clusters and unsign them
2293 for (Int_t sec=0;sec<fkNIS;sec++){
2294 for (Int_t row=0;row<fInnerSec->GetNRows();row++){
2295 AliTPCclusterMI *cl = fInnerSec[sec][row].GetClusters1();
2296 for (Int_t icl =0;icl< fInnerSec[sec][row].GetN1();icl++)
2297 // if (cl[icl].IsUsed(10))
2299 cl = fInnerSec[sec][row].GetClusters2();
2300 for (Int_t icl =0;icl< fInnerSec[sec][row].GetN2();icl++)
2301 //if (cl[icl].IsUsed(10))
2306 for (Int_t sec=0;sec<fkNOS;sec++){
2307 for (Int_t row=0;row<fOuterSec->GetNRows();row++){
2308 AliTPCclusterMI *cl = fOuterSec[sec][row].GetClusters1();
2309 for (Int_t icl =0;icl< fOuterSec[sec][row].GetN1();icl++)
2310 //if (cl[icl].IsUsed(10))
2312 cl = fOuterSec[sec][row].GetClusters2();
2313 for (Int_t icl =0;icl< fOuterSec[sec][row].GetN2();icl++)
2314 //if (cl[icl].IsUsed(10))
2323 void AliTPCtrackerMI::SignClusters(const TObjArray * arr, Float_t fnumber, Float_t fdensity)
2326 //sign clusters to be "used"
2328 // snumber and sdensity sign number of sigmas - bellow mean value to be accepted
2329 // loop over "primaries"
2343 Int_t nseed = arr->GetEntriesFast();
2344 for (Int_t i=0; i<nseed; i++) {
2345 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2349 if (!(pt->IsActive())) continue;
2350 Float_t dens = pt->GetNumberOfClusters()/Float_t(pt->GetNFoundable());
2351 if ( (dens>0.7) && (pt->GetNumberOfClusters()>70)){
2353 sumdens2+= dens*dens;
2354 sumn += pt->GetNumberOfClusters();
2355 sumn2 += pt->GetNumberOfClusters()*pt->GetNumberOfClusters();
2356 Float_t chi2 = pt->GetChi2()/pt->GetNumberOfClusters();
2359 sumchi2 +=chi2*chi2;
2364 Float_t mdensity = 0.9;
2365 Float_t meann = 130;
2366 Float_t meanchi = 1;
2367 Float_t sdensity = 0.1;
2368 Float_t smeann = 10;
2369 Float_t smeanchi =0.4;
2373 mdensity = sumdens/sum;
2375 meanchi = sumchi/sum;
2377 sdensity = sumdens2/sum-mdensity*mdensity;
2379 sdensity = TMath::Sqrt(sdensity);
2383 smeann = sumn2/sum-meann*meann;
2385 smeann = TMath::Sqrt(smeann);
2389 smeanchi = sumchi2/sum - meanchi*meanchi;
2391 smeanchi = TMath::Sqrt(smeanchi);
2397 //REMOVE SHORT DELTAS or tracks going out of sensitive volume of TPC
2399 for (Int_t i=0; i<nseed; i++) {
2400 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2404 if (pt->GetBSigned()) continue;
2405 if (pt->GetBConstrain()) continue;
2406 //if (!(pt->IsActive())) continue;
2408 Int_t found,foundable,shared;
2409 pt->GetClusterStatistic(0,160,found, foundable,shared);
2410 if (shared/float(found)>0.3) {
2411 if (shared/float(found)>0.9 ){
2412 //delete arr->RemoveAt(i);
2417 Bool_t isok =kFALSE;
2418 if ( (pt->GetNShared()/pt->GetNumberOfClusters()<0.5) &&pt->GetNumberOfClusters()>60)
2420 if ((TMath::Abs(1/pt->GetC())<100.) && (pt->GetNShared()/pt->GetNumberOfClusters()<0.7))
2422 if (TMath::Abs(pt->GetZ()/pt->GetX())>1.1)
2424 if ( (TMath::Abs(pt->GetSnp()>0.7) && pt->GetD(0,0)>60.))
2428 for (Int_t j=0; j<160; ++j) {
2429 Int_t index=pt->GetClusterIndex2(j);
2430 if (index<0) continue;
2431 AliTPCclusterMI *c= pt->GetClusterPointer(j);
2433 //if (!(c->IsUsed(10))) c->Use();
2440 Double_t maxchi = meanchi+2.*smeanchi;
2442 for (Int_t i=0; i<nseed; i++) {
2443 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2447 //if (!(pt->IsActive())) continue;
2448 if (pt->GetBSigned()) continue;
2449 Double_t chi = pt->GetChi2()/pt->GetNumberOfClusters();
2450 if (chi>maxchi) continue;
2453 Float_t dens = pt->GetNumberOfClusters()/Float_t(pt->GetNFoundable());
2455 //sign only tracks with enoug big density at the beginning
2457 if ((pt->GetDensityFirst(40)<0.75) && pt->GetNumberOfClusters()<meann) continue;
2460 Double_t mindens = TMath::Max(double(mdensity-sdensity*fdensity*bfactor),0.65);
2461 Double_t minn = TMath::Max(Int_t(meann-fnumber*smeann*bfactor),50);
2463 // if (pt->fBConstrain) mindens = TMath::Max(mdensity-sdensity*fdensity*bfactor,0.65);
2464 if ( (pt->GetRemoval()==10) && (pt->GetSnp()>0.8)&&(dens>mindens))
2467 if ((dens>mindens && pt->GetNumberOfClusters()>minn) && chi<maxchi ){
2468 //Int_t noc=pt->GetNumberOfClusters();
2469 pt->SetBSigned(kTRUE);
2470 for (Int_t j=0; j<160; ++j) {
2472 Int_t index=pt->GetClusterIndex2(j);
2473 if (index<0) continue;
2474 AliTPCclusterMI *c= pt->GetClusterPointer(j);
2476 // if (!(c->IsUsed(10))) c->Use();
2481 // gLastCheck = nseed;
2489 void AliTPCtrackerMI::StopNotActive(const TObjArray * arr, Int_t row0, Float_t th0, Float_t th1, Float_t th2) const
2491 // stop not active tracks
2492 // take th1 as threshold for number of founded to number of foundable on last 10 active rows
2493 // take th2 as threshold for number of founded to number of foundable on last 20 active rows
2494 Int_t nseed = arr->GetEntriesFast();
2496 for (Int_t i=0; i<nseed; i++) {
2497 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2501 if (!(pt->IsActive())) continue;
2502 StopNotActive(pt,row0,th0, th1,th2);
2508 void AliTPCtrackerMI::StopNotActive(AliTPCseed * seed, Int_t row0, Float_t th0, Float_t th1,
2511 // stop not active tracks
2512 // take th1 as threshold for number of founded to number of foundable on last 10 active rows
2513 // take th2 as threshold for number of founded to number of foundable on last 20 active rows
2516 Int_t foundable = 0;
2517 Int_t maxindex = seed->GetLastPoint(); //last foundable row
2518 if (seed->GetNFoundable()*th0 > seed->GetNumberOfClusters()) {
2519 seed->Desactivate(10) ;
2523 for (Int_t i=row0; i<maxindex; i++){
2524 Int_t index = seed->GetClusterIndex2(i);
2525 if (index!=-1) foundable++;
2527 if (foundable<=30) sumgood1++;
2528 if (foundable<=50) {
2535 if (foundable>=30.){
2536 if (sumgood1<(th1*30.)) seed->Desactivate(10);
2539 if (sumgood2<(th2*50.)) seed->Desactivate(10);
2543 Int_t AliTPCtrackerMI::RefitInward(AliESDEvent *event)
2546 // back propagation of ESD tracks
2549 const Int_t kMaxFriendTracks=2000;
2553 //PrepareForProlongation(fSeeds,1);
2554 PropagateForward2(fSeeds);
2555 RemoveUsed2(fSeeds,0.4,0.4,20);
2557 TObjArray arraySeed(fSeeds->GetEntries());
2558 for (Int_t i=0;i<fSeeds->GetEntries();i++) {
2559 arraySeed.AddAt(fSeeds->At(i),i);
2561 SignShared(&arraySeed);
2562 // FindCurling(fSeeds, event,2); // find multi found tracks
2563 FindSplitted(fSeeds, event,2); // find multi found tracks
2564 if (AliTPCReconstructor::StreamLevel()>2) FindMultiMC(fSeeds, fEvent,2); // find multi found tracks
2567 Int_t nseed = fSeeds->GetEntriesFast();
2568 for (Int_t i=0;i<nseed;i++){
2569 AliTPCseed * seed = (AliTPCseed*) fSeeds->UncheckedAt(i);
2570 if (!seed) continue;
2571 if (seed->GetKinkIndex(0)>0) UpdateKinkQualityD(seed); // update quality informations for kinks
2572 AliESDtrack *esd=event->GetTrack(i);
2574 if (seed->GetNumberOfClusters()<60 && seed->GetNumberOfClusters()<(esd->GetTPCclusters(0) -5)*0.8){
2575 AliExternalTrackParam paramIn;
2576 AliExternalTrackParam paramOut;
2577 Int_t ncl = seed->RefitTrack(seed,¶mIn,¶mOut);
2578 if (AliTPCReconstructor::StreamLevel()>0) {
2579 (*fDebugStreamer)<<"RecoverIn"<<
2583 "pout.="<<¶mOut<<
2588 seed->Set(paramIn.GetX(),paramIn.GetAlpha(),paramIn.GetParameter(),paramIn.GetCovariance());
2589 seed->SetNumberOfClusters(ncl);
2593 seed->PropagateTo(fkParam->GetInnerRadiusLow());
2594 seed->UpdatePoints();
2595 AddCovariance(seed);
2597 seed->CookdEdx(0.02,0.6);
2598 CookLabel(seed,0.1); //For comparison only
2599 esd->SetTPCClusterMap(seed->GetClusterMap());
2600 esd->SetTPCSharedMap(seed->GetSharedMap());
2602 if (AliTPCReconstructor::StreamLevel()>1 && seed!=0&&esd!=0) {
2603 TTreeSRedirector &cstream = *fDebugStreamer;
2610 if (seed->GetNumberOfClusters()>15){
2611 esd->UpdateTrackParams(seed,AliESDtrack::kTPCrefit);
2612 esd->SetTPCPoints(seed->GetPoints());
2613 esd->SetTPCPointsF(seed->GetNFoundable());
2614 Int_t ndedx = seed->GetNCDEDX(0);
2615 Float_t sdedx = seed->GetSDEDX(0);
2616 Float_t dedx = seed->GetdEdx();
2617 esd->SetTPCsignal(dedx, sdedx, ndedx);
2619 // add seed to the esd track in Calib level
2621 Bool_t storeFriend = gRandom->Rndm()<(kMaxFriendTracks)/Float_t(nseed);
2622 if (AliTPCReconstructor::StreamLevel()>0 &&storeFriend){
2623 AliTPCseed * seedCopy = new AliTPCseed(*seed, kTRUE);
2624 esd->AddCalibObject(seedCopy);
2629 //printf("problem\n");
2632 //FindKinks(fSeeds,event);
2633 Info("RefitInward","Number of refitted tracks %d",ntracks);
2638 Int_t AliTPCtrackerMI::PropagateBack(AliESDEvent *event)
2641 // back propagation of ESD tracks
2647 PropagateBack(fSeeds);
2648 RemoveUsed2(fSeeds,0.4,0.4,20);
2649 //FindCurling(fSeeds, fEvent,1);
2650 FindSplitted(fSeeds, event,1); // find multi found tracks
2651 if (AliTPCReconstructor::StreamLevel()>2) FindMultiMC(fSeeds, fEvent,1); // find multi found tracks
2654 Int_t nseed = fSeeds->GetEntriesFast();
2656 for (Int_t i=0;i<nseed;i++){
2657 AliTPCseed * seed = (AliTPCseed*) fSeeds->UncheckedAt(i);
2658 if (!seed) continue;
2659 if (seed->GetKinkIndex(0)<0) UpdateKinkQualityM(seed); // update quality informations for kinks
2660 seed->UpdatePoints();
2661 AddCovariance(seed);
2662 AliESDtrack *esd=event->GetTrack(i);
2663 if (seed->GetNumberOfClusters()<60 && seed->GetNumberOfClusters()<(esd->GetTPCclusters(0) -5)*0.8){
2664 AliExternalTrackParam paramIn;
2665 AliExternalTrackParam paramOut;
2666 Int_t ncl = seed->RefitTrack(seed,¶mIn,¶mOut);
2667 if (AliTPCReconstructor::StreamLevel()>0) {
2668 (*fDebugStreamer)<<"RecoverBack"<<
2672 "pout.="<<¶mOut<<
2677 seed->Set(paramOut.GetX(),paramOut.GetAlpha(),paramOut.GetParameter(),paramOut.GetCovariance());
2678 seed->SetNumberOfClusters(ncl);
2681 seed->CookdEdx(0.02,0.6);
2682 CookLabel(seed,0.1); //For comparison only
2683 if (seed->GetNumberOfClusters()>15){
2684 esd->UpdateTrackParams(seed,AliESDtrack::kTPCout);
2685 esd->SetTPCPoints(seed->GetPoints());
2686 esd->SetTPCPointsF(seed->GetNFoundable());
2687 Int_t ndedx = seed->GetNCDEDX(0);
2688 Float_t sdedx = seed->GetSDEDX(0);
2689 Float_t dedx = seed->GetdEdx();
2690 esd->SetTPCsignal(dedx, sdedx, ndedx);
2692 Int_t eventnumber = event->GetEventNumberInFile();// patch 28 fev 06
2693 // This is most likely NOT the event number you'd like to use. It has nothing to do with the 'real' event number
2694 if (AliTPCReconstructor::StreamLevel()>1 && esd) {
2695 (*fDebugStreamer)<<"Cback"<<
2698 "EventNrInFile="<<eventnumber<<
2703 //FindKinks(fSeeds,event);
2704 Info("PropagateBack","Number of back propagated tracks %d",ntracks);
2711 void AliTPCtrackerMI::DeleteSeeds()
2716 Int_t nseed = fSeeds->GetEntriesFast();
2717 for (Int_t i=0;i<nseed;i++){
2718 AliTPCseed * seed = (AliTPCseed*)fSeeds->At(i);
2719 if (seed) delete fSeeds->RemoveAt(i);
2726 void AliTPCtrackerMI::ReadSeeds(AliESDEvent *event, Int_t direction)
2729 //read seeds from the event
2731 Int_t nentr=event->GetNumberOfTracks();
2733 Info("ReadSeeds", "Number of ESD tracks: %d\n", nentr);
2738 fSeeds = new TObjArray(nentr);
2742 for (Int_t i=0; i<nentr; i++) {
2743 AliESDtrack *esd=event->GetTrack(i);
2744 ULong_t status=esd->GetStatus();
2745 if (!(status&AliESDtrack::kTPCin)) continue;
2746 AliTPCtrack t(*esd);
2747 t.SetNumberOfClusters(0);
2748 // AliTPCseed *seed = new AliTPCseed(t,t.GetAlpha());
2749 AliTPCseed *seed = new AliTPCseed(t/*,t.GetAlpha()*/);
2750 seed->SetUniqueID(esd->GetID());
2751 AddCovariance(seed); //add systematic ucertainty
2752 for (Int_t ikink=0;ikink<3;ikink++) {
2753 Int_t index = esd->GetKinkIndex(ikink);
2754 seed->GetKinkIndexes()[ikink] = index;
2755 if (index==0) continue;
2756 index = TMath::Abs(index);
2757 AliESDkink * kink = fEvent->GetKink(index-1);
2758 if (kink&&esd->GetKinkIndex(ikink)<0){
2759 if ((status & AliESDtrack::kTRDrefit) != 0) kink->SetStatus(1,2);
2760 if ((status & AliESDtrack::kITSout) != 0) kink->SetStatus(1,0);
2762 if (kink&&esd->GetKinkIndex(ikink)>0){
2763 if ((status & AliESDtrack::kTRDrefit) != 0) kink->SetStatus(1,6);
2764 if ((status & AliESDtrack::kITSout) != 0) kink->SetStatus(1,4);
2768 if (((status&AliESDtrack::kITSout)==0)&&(direction==1)) seed->ResetCovariance(10.);
2769 if ( direction ==2 &&(status & AliESDtrack::kTRDrefit) == 0 ) seed->ResetCovariance(10.);
2770 //if ( direction ==2 && ((status & AliESDtrack::kTPCout) == 0) ) {
2771 // fSeeds->AddAt(0,i);
2775 if ( direction ==2 &&(status & AliESDtrack::kTRDrefit) > 0 ) {
2776 Double_t par0[5],par1[5],alpha,x;
2777 esd->GetInnerExternalParameters(alpha,x,par0);
2778 esd->GetExternalParameters(x,par1);
2779 Double_t delta1 = TMath::Abs(par0[4]-par1[4])/(0.000000001+TMath::Abs(par0[4]+par1[4]));
2780 Double_t delta2 = TMath::Abs(par0[3]-par1[3]);
2782 if (esd->GetTRDncls()>0) trdchi2 = esd->GetTRDchi2()/esd->GetTRDncls();
2783 //reset covariance if suspicious
2784 if ( (delta1>0.1) || (delta2>0.006) ||trdchi2>7.)
2785 seed->ResetCovariance(10.);
2790 // rotate to the local coordinate system
2792 fSectors=fInnerSec; fN=fkNIS;
2793 Double_t alpha=seed->GetAlpha() - fSectors->GetAlphaShift();
2794 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
2795 if (alpha < 0. ) alpha += 2.*TMath::Pi();
2796 Int_t ns=Int_t(alpha/fSectors->GetAlpha())%fN;
2797 alpha =ns*fSectors->GetAlpha() + fSectors->GetAlphaShift();
2798 if (alpha<-TMath::Pi()) alpha += 2*TMath::Pi();
2799 if (alpha>=TMath::Pi()) alpha -= 2*TMath::Pi();
2800 alpha-=seed->GetAlpha();
2801 if (!seed->Rotate(alpha)) {
2807 if (esd->GetKinkIndex(0)<=0){
2808 for (Int_t irow=0;irow<160;irow++){
2809 Int_t index = seed->GetClusterIndex2(irow);
2812 AliTPCclusterMI * cl = GetClusterMI(index);
2813 seed->SetClusterPointer(irow,cl);
2815 if ((index & 0x8000)==0){
2816 cl->Use(10); // accepted cluster
2818 cl->Use(6); // close cluster not accepted
2821 Info("ReadSeeds","Not found cluster");
2826 fSeeds->AddAt(seed,i);
2832 //_____________________________________________________________________________
2833 void AliTPCtrackerMI::MakeSeeds3(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, Float_t cuts[4],
2834 Float_t deltay, Int_t ddsec) {
2835 //-----------------------------------------------------------------
2836 // This function creates track seeds.
2837 // SEEDING WITH VERTEX CONSTRAIN
2838 //-----------------------------------------------------------------
2839 // cuts[0] - fP4 cut
2840 // cuts[1] - tan(phi) cut
2841 // cuts[2] - zvertex cut
2842 // cuts[3] - fP3 cut
2850 Double_t x[5], c[15];
2851 // Int_t di = i1-i2;
2853 AliTPCseed * seed = new AliTPCseed();
2854 Double_t alpha=fSectors->GetAlpha(), shift=fSectors->GetAlphaShift();
2855 Double_t cs=cos(alpha), sn=sin(alpha);
2857 // Double_t x1 =fOuterSec->GetX(i1);
2858 //Double_t xx2=fOuterSec->GetX(i2);
2860 Double_t x1 =GetXrow(i1);
2861 Double_t xx2=GetXrow(i2);
2863 Double_t x3=GetX(), y3=GetY(), z3=GetZ();
2865 Int_t imiddle = (i2+i1)/2; //middle pad row index
2866 Double_t xm = GetXrow(imiddle); // radius of middle pad-row
2867 const AliTPCtrackerRow& krm=GetRow(sec,imiddle); //middle pad -row
2871 const AliTPCtrackerRow& kr1=GetRow(ns,i1);
2872 Double_t ymax = GetMaxY(i1)-kr1.GetDeadZone()-1.5;
2873 Double_t ymaxm = GetMaxY(imiddle)-kr1.GetDeadZone()-1.5;
2876 // change cut on curvature if it can't reach this layer
2877 // maximal curvature set to reach it
2878 Double_t dvertexmax = TMath::Sqrt((x1-x3)*(x1-x3)+(ymax+5-y3)*(ymax+5-y3));
2879 if (dvertexmax*0.5*cuts[0]>0.85){
2880 cuts[0] = 0.85/(dvertexmax*0.5+1.);
2882 Double_t r2min = 1/(cuts[0]*cuts[0]); //minimal square of radius given by cut
2885 if (deltay>0) ddsec = 0;
2886 // loop over clusters
2887 for (Int_t is=0; is < kr1; is++) {
2889 if (kr1[is]->IsUsed(10)) continue;
2890 Double_t y1=kr1[is]->GetY(), z1=kr1[is]->GetZ();
2891 //if (TMath::Abs(y1)>ymax) continue;
2893 if (deltay>0 && TMath::Abs(ymax-TMath::Abs(y1))> deltay ) continue; // seed only at the edge
2895 // find possible directions
2896 Float_t anglez = (z1-z3)/(x1-x3);
2897 Float_t extraz = z1 - anglez*(x1-xx2); // extrapolated z
2900 //find rotation angles relative to line given by vertex and point 1
2901 Double_t dvertex2 = (x1-x3)*(x1-x3)+(y1-y3)*(y1-y3);
2902 Double_t dvertex = TMath::Sqrt(dvertex2);
2903 Double_t angle13 = TMath::ATan((y1-y3)/(x1-x3));
2904 Double_t cs13 = cos(-angle13), sn13 = sin(-angle13);
2907 // loop over 2 sectors
2913 Double_t dddz1=0; // direction of delta inclination in z axis
2920 for (Int_t dsec = dsec1; dsec<=dsec2;dsec++){
2921 Int_t sec2 = sec + dsec;
2923 // AliTPCtrackerRow& kr2 = fOuterSec[(sec2+fkNOS)%fkNOS][i2];
2924 //AliTPCtrackerRow& kr2m = fOuterSec[(sec2+fkNOS)%fkNOS][imiddle];
2925 AliTPCtrackerRow& kr2 = GetRow((sec2+fkNOS)%fkNOS,i2);
2926 AliTPCtrackerRow& kr2m = GetRow((sec2+fkNOS)%fkNOS,imiddle);
2927 Int_t index1 = TMath::Max(kr2.Find(extraz-0.6-dddz1*TMath::Abs(z1)*0.05)-1,0);
2928 Int_t index2 = TMath::Min(kr2.Find(extraz+0.6+dddz2*TMath::Abs(z1)*0.05)+1,kr2);
2930 // rotation angles to p1-p3
2931 Double_t cs13r = cos(-angle13+dsec*alpha)/dvertex, sn13r = sin(-angle13+dsec*alpha)/dvertex;
2932 Double_t x2, y2, z2;
2934 // Double_t dymax = maxangle*TMath::Abs(x1-xx2);
2937 Double_t dxx0 = (xx2-x3)*cs13r;
2938 Double_t dyy0 = (xx2-x3)*sn13r;
2939 for (Int_t js=index1; js < index2; js++) {
2940 const AliTPCclusterMI *kcl = kr2[js];
2941 if (kcl->IsUsed(10)) continue;
2943 //calcutate parameters
2945 Double_t yy0 = dyy0 +(kcl->GetY()-y3)*cs13r;
2947 if (TMath::Abs(yy0)<0.000001) continue;
2948 Double_t xx0 = dxx0 -(kcl->GetY()-y3)*sn13r;
2949 Double_t y0 = 0.5*(xx0*xx0+yy0*yy0-xx0)/yy0;
2950 Double_t r02 = (0.25+y0*y0)*dvertex2;
2951 //curvature (radius) cut
2952 if (r02<r2min) continue;
2956 Double_t c0 = 1/TMath::Sqrt(r02);
2960 //Double_t dfi0 = 2.*TMath::ASin(dvertex*c0*0.5);
2961 //Double_t dfi1 = 2.*TMath::ASin(TMath::Sqrt(yy0*yy0+(1-xx0)*(1-xx0))*dvertex*c0*0.5);
2962 Double_t dfi0 = 2.*AliTPCFastMath::FastAsin(dvertex*c0*0.5);
2963 Double_t dfi1 = 2.*AliTPCFastMath::FastAsin(TMath::Sqrt(yy0*yy0+(1-xx0)*(1-xx0))*dvertex*c0*0.5);
2966 Double_t z0 = kcl->GetZ();
2967 Double_t zzzz2 = z1-(z1-z3)*dfi1/dfi0;
2968 if (TMath::Abs(zzzz2-z0)>0.5) continue;
2971 Double_t dip = (z1-z0)*c0/dfi1;
2972 Double_t x0 = (0.5*cs13+y0*sn13)*dvertex*c0;
2983 x2= xx2*cs-y2*sn*dsec;
2984 y2=+xx2*sn*dsec+y2*cs;
2994 // do we have cluster at the middle ?
2996 GetProlongation(x1,xm,x,ym,zm);
2998 AliTPCclusterMI * cm=0;
2999 if (TMath::Abs(ym)-ymaxm<0){
3000 cm = krm.FindNearest2(ym,zm,1.0,0.6,dummy);
3001 if ((!cm) || (cm->IsUsed(10))) {
3006 // rotate y1 to system 0
3007 // get state vector in rotated system
3008 Double_t yr1 = (-0.5*sn13+y0*cs13)*dvertex*c0;
3009 Double_t xr2 = x0*cs+yr1*sn*dsec;
3010 Double_t xr[5]={kcl->GetY(),kcl->GetZ(), xr2, dip, c0};
3012 GetProlongation(xx2,xm,xr,ym,zm);
3013 if (TMath::Abs(ym)-ymaxm<0){
3014 cm = kr2m.FindNearest2(ym,zm,1.0,0.6,dummy);
3015 if ((!cm) || (cm->IsUsed(10))) {
3025 dym = ym - cm->GetY();
3026 dzm = zm - cm->GetZ();
3033 Double_t sy1=kr1[is]->GetSigmaY2()*2., sz1=kr1[is]->GetSigmaZ2()*2.;
3034 Double_t sy2=kcl->GetSigmaY2()*2., sz2=kcl->GetSigmaZ2()*2.;
3035 //Double_t sy3=400*3./12., sy=0.1, sz=0.1;
3036 Double_t sy3=25000*x[4]*x[4]+0.1, sy=0.1, sz=0.1;
3037 //Double_t sy3=25000*x[4]*x[4]*60+0.5, sy=0.1, sz=0.1;
3039 Double_t f40=(F1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
3040 Double_t f42=(F1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
3041 Double_t f43=(F1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
3042 Double_t f20=(F2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
3043 Double_t f22=(F2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
3044 Double_t f23=(F2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
3046 Double_t f30=(F3(x1,y1+sy,x2,y2,z1,z2)-x[3])/sy;
3047 Double_t f31=(F3(x1,y1,x2,y2,z1+sz,z2)-x[3])/sz;
3048 Double_t f32=(F3(x1,y1,x2,y2+sy,z1,z2)-x[3])/sy;
3049 Double_t f34=(F3(x1,y1,x2,y2,z1,z2+sz)-x[3])/sz;
3053 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3054 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3055 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3056 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3057 c[13]=f30*sy1*f40+f32*sy2*f42;
3058 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3060 // if (!BuildSeed(kr1[is],kcl,0,x1,x2,x3,x,c)) continue;
3062 UInt_t index=kr1.GetIndex(is);
3063 seed->~AliTPCseed(); // this does not set the pointer to 0...
3064 AliTPCseed *track=new(seed) AliTPCseed(x1, ns*alpha+shift, x, c, index);
3066 track->SetIsSeeding(kTRUE);
3067 track->SetSeed1(i1);
3068 track->SetSeed2(i2);
3069 track->SetSeedType(3);
3073 FollowProlongation(*track, (i1+i2)/2,1);
3074 Int_t foundable,found,shared;
3075 track->GetClusterStatistic((i1+i2)/2,i1, found, foundable, shared, kTRUE);
3076 if ((found<0.55*foundable) || shared>0.5*found || (track->GetSigmaY2()+track->GetSigmaZ2())>0.5){
3078 seed->~AliTPCseed();
3084 FollowProlongation(*track, i2,1);
3088 track->SetBConstrain(1);
3089 // track->fLastPoint = i1+fInnerSec->GetNRows(); // first cluster in track position
3090 track->SetLastPoint(i1); // first cluster in track position
3091 track->SetFirstPoint(track->GetLastPoint());
3093 if (track->GetNumberOfClusters()<(i1-i2)*0.5 ||
3094 track->GetNumberOfClusters() < track->GetNFoundable()*0.6 ||
3095 track->GetNShared()>0.4*track->GetNumberOfClusters() ) {
3097 seed->~AliTPCseed();
3101 // Z VERTEX CONDITION
3102 Double_t zv, bz=GetBz();
3103 if ( !track->GetZAt(0.,bz,zv) ) continue;
3104 if (TMath::Abs(zv-z3)>cuts[2]) {
3105 FollowProlongation(*track, TMath::Max(i2-20,0));
3106 if ( !track->GetZAt(0.,bz,zv) ) continue;
3107 if (TMath::Abs(zv-z3)>cuts[2]){
3108 FollowProlongation(*track, TMath::Max(i2-40,0));
3109 if ( !track->GetZAt(0.,bz,zv) ) continue;
3110 if (TMath::Abs(zv-z3)>cuts[2] &&(track->GetNumberOfClusters() > track->GetNFoundable()*0.7)){
3111 // make seed without constrain
3112 AliTPCseed * track2 = MakeSeed(track,0.2,0.5,1.);
3113 FollowProlongation(*track2, i2,1);
3114 track2->SetBConstrain(kFALSE);
3115 track2->SetSeedType(1);
3116 arr->AddLast(track2);
3118 seed->~AliTPCseed();
3123 seed->~AliTPCseed();
3130 track->SetSeedType(0);
3131 arr->AddLast(track);
3132 seed = new AliTPCseed;
3134 // don't consider other combinations
3135 if (track->GetNumberOfClusters() > track->GetNFoundable()*0.8)
3141 Info("MakeSeeds3","\nSeeding statistic:\t%d\t%d\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin,nout1,nout2);
3147 void AliTPCtrackerMI::MakeSeeds5(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, Float_t cuts[4],
3152 //-----------------------------------------------------------------
3153 // This function creates track seeds.
3154 //-----------------------------------------------------------------
3155 // cuts[0] - fP4 cut
3156 // cuts[1] - tan(phi) cut
3157 // cuts[2] - zvertex cut
3158 // cuts[3] - fP3 cut
3168 Double_t x[5], c[15];
3170 // make temporary seed
3171 AliTPCseed * seed = new AliTPCseed;
3172 Double_t alpha=fOuterSec->GetAlpha(), shift=fOuterSec->GetAlphaShift();
3173 // Double_t cs=cos(alpha), sn=sin(alpha);
3178 Double_t x1 = GetXrow(i1-1);
3179 const AliTPCtrackerRow& kr1=GetRow(sec,i1-1);
3180 Double_t y1max = GetMaxY(i1-1)-kr1.GetDeadZone()-1.5;
3182 Double_t x1p = GetXrow(i1);
3183 const AliTPCtrackerRow& kr1p=GetRow(sec,i1);
3185 Double_t x1m = GetXrow(i1-2);
3186 const AliTPCtrackerRow& kr1m=GetRow(sec,i1-2);
3189 //last 3 padrow for seeding
3190 AliTPCtrackerRow& kr3 = GetRow((sec+fkNOS)%fkNOS,i1-7);
3191 Double_t x3 = GetXrow(i1-7);
3192 // Double_t y3max= GetMaxY(i1-7)-kr3.fDeadZone-1.5;
3194 AliTPCtrackerRow& kr3p = GetRow((sec+fkNOS)%fkNOS,i1-6);
3195 Double_t x3p = GetXrow(i1-6);
3197 AliTPCtrackerRow& kr3m = GetRow((sec+fkNOS)%fkNOS,i1-8);
3198 Double_t x3m = GetXrow(i1-8);
3203 Int_t im = i1-4; //middle pad row index
3204 Double_t xm = GetXrow(im); // radius of middle pad-row
3205 const AliTPCtrackerRow& krm=GetRow(sec,im); //middle pad -row
3206 // Double_t ymmax = GetMaxY(im)-kr1.fDeadZone-1.5;
3209 Double_t deltax = x1-x3;
3210 Double_t dymax = deltax*cuts[1];
3211 Double_t dzmax = deltax*cuts[3];
3213 // loop over clusters
3214 for (Int_t is=0; is < kr1; is++) {
3216 if (kr1[is]->IsUsed(10)) continue;
3217 Double_t y1=kr1[is]->GetY(), z1=kr1[is]->GetZ();
3219 if (deltay>0 && TMath::Abs(y1max-TMath::Abs(y1))> deltay ) continue; // seed only at the edge
3221 Int_t index1 = TMath::Max(kr3.Find(z1-dzmax)-1,0);
3222 Int_t index2 = TMath::Min(kr3.Find(z1+dzmax)+1,kr3);
3228 for (Int_t js=index1; js < index2; js++) {
3229 const AliTPCclusterMI *kcl = kr3[js];
3230 if (kcl->IsUsed(10)) continue;
3232 // apply angular cuts
3233 if (TMath::Abs(y1-y3)>dymax) continue;
3236 if (TMath::Abs(z1-z3)>dzmax) continue;
3238 Double_t angley = (y1-y3)/(x1-x3);
3239 Double_t anglez = (z1-z3)/(x1-x3);
3241 Double_t erry = TMath::Abs(angley)*(x1-x1m)*0.5+0.5;
3242 Double_t errz = TMath::Abs(anglez)*(x1-x1m)*0.5+0.5;
3244 Double_t yyym = angley*(xm-x1)+y1;
3245 Double_t zzzm = anglez*(xm-x1)+z1;
3247 const AliTPCclusterMI *kcm = krm.FindNearest2(yyym,zzzm,erry,errz,index);
3249 if (kcm->IsUsed(10)) continue;
3251 erry = TMath::Abs(angley)*(x1-x1m)*0.4+0.5;
3252 errz = TMath::Abs(anglez)*(x1-x1m)*0.4+0.5;
3259 // look around first
3260 const AliTPCclusterMI *kc1m = kr1m.FindNearest2(angley*(x1m-x1)+y1,
3266 if (kc1m->IsUsed(10)) used++;
3268 const AliTPCclusterMI *kc1p = kr1p.FindNearest2(angley*(x1p-x1)+y1,
3274 if (kc1p->IsUsed(10)) used++;
3276 if (used>1) continue;
3277 if (found<1) continue;
3281 const AliTPCclusterMI *kc3m = kr3m.FindNearest2(angley*(x3m-x3)+y3,
3287 if (kc3m->IsUsed(10)) used++;
3291 const AliTPCclusterMI *kc3p = kr3p.FindNearest2(angley*(x3p-x3)+y3,
3297 if (kc3p->IsUsed(10)) used++;
3301 if (used>1) continue;
3302 if (found<3) continue;
3312 x[4]=F1(x1,y1,x2,y2,x3,y3);
3313 //if (TMath::Abs(x[4]) >= cuts[0]) continue;
3316 x[2]=F2(x1,y1,x2,y2,x3,y3);
3319 x[3]=F3n(x1,y1,x2,y2,z1,z2,x[4]);
3320 //if (TMath::Abs(x[3]) > cuts[3]) continue;
3324 Double_t sy1=0.1, sz1=0.1;
3325 Double_t sy2=0.1, sz2=0.1;
3326 Double_t sy3=0.1, sy=0.1, sz=0.1;
3328 Double_t f40=(F1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
3329 Double_t f42=(F1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
3330 Double_t f43=(F1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
3331 Double_t f20=(F2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
3332 Double_t f22=(F2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
3333 Double_t f23=(F2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
3335 Double_t f30=(F3(x1,y1+sy,x2,y2,z1,z2)-x[3])/sy;
3336 Double_t f31=(F3(x1,y1,x2,y2,z1+sz,z2)-x[3])/sz;
3337 Double_t f32=(F3(x1,y1,x2,y2+sy,z1,z2)-x[3])/sy;
3338 Double_t f34=(F3(x1,y1,x2,y2,z1,z2+sz)-x[3])/sz;
3342 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3343 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3344 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3345 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3346 c[13]=f30*sy1*f40+f32*sy2*f42;
3347 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3349 // if (!BuildSeed(kr1[is],kcl,0,x1,x2,x3,x,c)) continue;
3351 index=kr1.GetIndex(is);
3352 seed->~AliTPCseed();
3353 AliTPCseed *track=new(seed) AliTPCseed(x1, sec*alpha+shift, x, c, index);
3355 track->SetIsSeeding(kTRUE);
3358 FollowProlongation(*track, i1-7,1);
3359 if (track->GetNumberOfClusters() < track->GetNFoundable()*0.75 ||
3360 track->GetNShared()>0.6*track->GetNumberOfClusters() || ( track->GetSigmaY2()+ track->GetSigmaZ2())>0.6){
3362 seed->~AliTPCseed();
3368 FollowProlongation(*track, i2,1);
3369 track->SetBConstrain(0);
3370 track->SetLastPoint(i1+fInnerSec->GetNRows()); // first cluster in track position
3371 track->SetFirstPoint(track->GetLastPoint());
3373 if (track->GetNumberOfClusters()<(i1-i2)*0.5 ||
3374 track->GetNumberOfClusters()<track->GetNFoundable()*0.7 ||
3375 track->GetNShared()>2. || track->GetChi2()/track->GetNumberOfClusters()>6 || ( track->GetSigmaY2()+ track->GetSigmaZ2())>0.5 ) {
3377 seed->~AliTPCseed();
3382 FollowProlongation(*track, TMath::Max(i2-10,0),1);
3383 AliTPCseed * track2 = MakeSeed(track,0.2,0.5,0.9);
3384 FollowProlongation(*track2, i2,1);
3385 track2->SetBConstrain(kFALSE);
3386 track2->SetSeedType(4);
3387 arr->AddLast(track2);
3389 seed->~AliTPCseed();
3393 //arr->AddLast(track);
3394 //seed = new AliTPCseed;
3400 Info("MakeSeeds5","\nSeeding statiistic:\t%d\t%d\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin,nout1,nout2,nout3);
3406 //_____________________________________________________________________________
3407 void AliTPCtrackerMI::MakeSeeds2(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, Float_t */*cuts[4]*/,
3408 Float_t deltay, Bool_t /*bconstrain*/) {
3409 //-----------------------------------------------------------------
3410 // This function creates track seeds - without vertex constraint
3411 //-----------------------------------------------------------------
3412 // cuts[0] - fP4 cut - not applied
3413 // cuts[1] - tan(phi) cut
3414 // cuts[2] - zvertex cut - not applied
3415 // cuts[3] - fP3 cut
3425 Double_t alpha=fOuterSec->GetAlpha(), shift=fOuterSec->GetAlphaShift();
3426 // Double_t cs=cos(alpha), sn=sin(alpha);
3427 Int_t row0 = (i1+i2)/2;
3428 Int_t drow = (i1-i2)/2;
3429 const AliTPCtrackerRow& kr0=fSectors[sec][row0];
3430 AliTPCtrackerRow * kr=0;
3432 AliTPCpolyTrack polytrack;
3433 Int_t nclusters=fSectors[sec][row0];
3434 AliTPCseed * seed = new AliTPCseed;
3439 for (Int_t is=0; is < nclusters; is++) { //LOOP over clusters
3441 Int_t nfoundable =0;
3442 for (Int_t iter =1; iter<2; iter++){ //iterations
3443 const AliTPCtrackerRow& krm=fSectors[sec][row0-iter];
3444 const AliTPCtrackerRow& krp=fSectors[sec][row0+iter];
3445 const AliTPCclusterMI * cl= kr0[is];
3447 if (cl->IsUsed(10)) {
3453 Double_t x = kr0.GetX();
3454 // Initialization of the polytrack
3459 Double_t y0= cl->GetY();
3460 Double_t z0= cl->GetZ();
3464 Double_t ymax = fSectors->GetMaxY(row0)-kr0.GetDeadZone()-1.5;
3465 if (deltay>0 && TMath::Abs(ymax-TMath::Abs(y0))> deltay ) continue; // seed only at the edge
3467 erry = (0.5)*cl->GetSigmaY2()/TMath::Sqrt(cl->GetQ())*6;
3468 errz = (0.5)*cl->GetSigmaZ2()/TMath::Sqrt(cl->GetQ())*6;
3469 polytrack.AddPoint(x,y0,z0,erry, errz);
3472 if (cl->IsUsed(10)) sumused++;
3475 Float_t roady = (5*TMath::Sqrt(cl->GetSigmaY2()+0.2)+1.)*iter;
3476 Float_t roadz = (5*TMath::Sqrt(cl->GetSigmaZ2()+0.2)+1.)*iter;
3479 AliTPCclusterMI * cl1 = krm.FindNearest(y0,z0,roady,roadz);
3480 if (cl1 && TMath::Abs(ymax-TMath::Abs(y0))) {
3481 erry = (0.5)*cl1->GetSigmaY2()/TMath::Sqrt(cl1->GetQ())*3;
3482 errz = (0.5)*cl1->GetSigmaZ2()/TMath::Sqrt(cl1->GetQ())*3;
3483 if (cl1->IsUsed(10)) sumused++;
3484 polytrack.AddPoint(x,cl1->GetY(),cl1->GetZ(),erry,errz);
3488 AliTPCclusterMI * cl2 = krp.FindNearest(y0,z0,roady,roadz);
3490 erry = (0.5)*cl2->GetSigmaY2()/TMath::Sqrt(cl2->GetQ())*3;
3491 errz = (0.5)*cl2->GetSigmaZ2()/TMath::Sqrt(cl2->GetQ())*3;
3492 if (cl2->IsUsed(10)) sumused++;
3493 polytrack.AddPoint(x,cl2->GetY(),cl2->GetZ(),erry,errz);
3496 if (sumused>0) continue;
3498 polytrack.UpdateParameters();
3504 nfoundable = polytrack.GetN();
3505 nfound = nfoundable;
3507 for (Int_t ddrow = iter+1; ddrow<drow;ddrow++){
3508 Float_t maxdist = 0.8*(1.+3./(ddrow));
3509 for (Int_t delta = -1;delta<=1;delta+=2){
3510 Int_t row = row0+ddrow*delta;
3511 kr = &(fSectors[sec][row]);
3512 Double_t xn = kr->GetX();
3513 Double_t ymax1 = fSectors->GetMaxY(row)-kr->GetDeadZone()-1.5;
3514 polytrack.GetFitPoint(xn,yn,zn);
3515 if (TMath::Abs(yn)>ymax1) continue;
3517 AliTPCclusterMI * cln = kr->FindNearest(yn,zn,roady,roadz);
3519 Float_t dist = TMath::Sqrt( (yn-cln->GetY())*(yn-cln->GetY())+(zn-cln->GetZ())*(zn-cln->GetZ()));
3522 erry = (dist+0.3)*cln->GetSigmaY2()/TMath::Sqrt(cln->GetQ())*(1.+1./(ddrow));
3523 errz = (dist+0.3)*cln->GetSigmaZ2()/TMath::Sqrt(cln->GetQ())*(1.+1./(ddrow));
3524 if (cln->IsUsed(10)) {
3525 // printf("used\n");
3533 polytrack.AddPoint(xn,cln->GetY(),cln->GetZ(),erry, errz);
3538 if ( (sumused>3) || (sumused>0.5*nfound) || (nfound<0.6*nfoundable)) break;
3539 polytrack.UpdateParameters();
3542 if ( (sumused>3) || (sumused>0.5*nfound)) {
3543 //printf("sumused %d\n",sumused);
3548 polytrack.GetFitDerivation(kr0.GetX(),dy,dz);
3549 AliTPCpolyTrack track2;
3551 polytrack.Refit(track2,0.5+TMath::Abs(dy)*0.3,0.4+TMath::Abs(dz)*0.3);
3552 if (track2.GetN()<0.5*nfoundable) continue;
3555 if ((nfound>0.6*nfoundable) &&( nfoundable>0.4*(i1-i2))) {
3557 // test seed with and without constrain
3558 for (Int_t constrain=0; constrain<=0;constrain++){
3559 // add polytrack candidate
3561 Double_t x[5], c[15];
3562 Double_t x1,x2,x3,y1,y2,y3,z1,z2,z3;
3563 track2.GetBoundaries(x3,x1);
3565 track2.GetFitPoint(x1,y1,z1);
3566 track2.GetFitPoint(x2,y2,z2);
3567 track2.GetFitPoint(x3,y3,z3);
3569 //is track pointing to the vertex ?
3572 polytrack.GetFitPoint(x0,y0,z0);
3585 x[4]=F1(x1,y1,x2,y2,x3,y3);
3587 // if (TMath::Abs(x[4]) >= cuts[0]) continue; //
3588 x[2]=F2(x1,y1,x2,y2,x3,y3);
3590 //if (TMath::Abs(x[4]*x1-x[2]) >= cuts[1]) continue;
3591 //x[3]=F3(x1,y1,x2,y2,z1,z2);
3592 x[3]=F3n(x1,y1,x3,y3,z1,z3,x[4]);
3593 //if (TMath::Abs(x[3]) > cuts[3]) continue;
3596 Double_t sy =0.1, sz =0.1;
3597 Double_t sy1=0.02, sz1=0.02;
3598 Double_t sy2=0.02, sz2=0.02;
3602 sy3=25000*x[4]*x[4]+0.1, sy=0.1, sz=0.1;
3605 Double_t f40=(F1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
3606 Double_t f42=(F1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
3607 Double_t f43=(F1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
3608 Double_t f20=(F2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
3609 Double_t f22=(F2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
3610 Double_t f23=(F2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
3612 Double_t f30=(F3(x1,y1+sy,x3,y3,z1,z3)-x[3])/sy;
3613 Double_t f31=(F3(x1,y1,x3,y3,z1+sz,z3)-x[3])/sz;
3614 Double_t f32=(F3(x1,y1,x3,y3+sy,z1,z3)-x[3])/sy;
3615 Double_t f34=(F3(x1,y1,x3,y3,z1,z3+sz)-x[3])/sz;
3620 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3621 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3622 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3623 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3624 c[13]=f30*sy1*f40+f32*sy2*f42;
3625 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3627 //Int_t row1 = fSectors->GetRowNumber(x1);
3628 Int_t row1 = GetRowNumber(x1);
3632 seed->~AliTPCseed();
3633 AliTPCseed *track=new(seed) AliTPCseed(x1,sec*alpha+shift,x,c,index);
3634 track->SetIsSeeding(kTRUE);
3635 Int_t rc=FollowProlongation(*track, i2);
3636 if (constrain) track->SetBConstrain(1);
3638 track->SetBConstrain(0);
3639 track->SetLastPoint(row1+fInnerSec->GetNRows()); // first cluster in track position
3640 track->SetFirstPoint(track->GetLastPoint());
3642 if (rc==0 || track->GetNumberOfClusters()<(i1-i2)*0.5 ||
3643 track->GetNumberOfClusters() < track->GetNFoundable()*0.6 ||
3644 track->GetNShared()>0.4*track->GetNumberOfClusters()) {
3647 seed->~AliTPCseed();
3650 arr->AddLast(track);
3651 seed = new AliTPCseed;
3655 } // if accepted seed
3658 Info("MakeSeeds2","\nSeeding statiistic:\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin3);
3664 AliTPCseed *AliTPCtrackerMI::MakeSeed(AliTPCseed *track, Float_t r0, Float_t r1, Float_t r2)
3668 //reseed using track points
3669 Int_t p0 = int(r0*track->GetNumberOfClusters()); // point 0
3670 Int_t p1 = int(r1*track->GetNumberOfClusters());
3671 Int_t p2 = int(r2*track->GetNumberOfClusters()); // last point
3673 Double_t x0[3],x1[3],x2[3];
3674 for (Int_t i=0;i<3;i++){
3680 // find track position at given ratio of the length
3681 Int_t sec0=0, sec1=0, sec2=0;
3684 for (Int_t i=0;i<160;i++){
3685 if (track->GetClusterPointer(i)){
3687 AliTPCTrackerPoint *trpoint =track->GetTrackPoint(i);
3688 if ( (index<p0) || x0[0]<0 ){
3689 if (trpoint->GetX()>1){
3690 clindex = track->GetClusterIndex2(i);
3692 x0[0] = trpoint->GetX();
3693 x0[1] = trpoint->GetY();
3694 x0[2] = trpoint->GetZ();
3695 sec0 = ((clindex&0xff000000)>>24)%18;
3700 if ( (index<p1) &&(trpoint->GetX()>1)){
3701 clindex = track->GetClusterIndex2(i);
3703 x1[0] = trpoint->GetX();
3704 x1[1] = trpoint->GetY();
3705 x1[2] = trpoint->GetZ();
3706 sec1 = ((clindex&0xff000000)>>24)%18;
3709 if ( (index<p2) &&(trpoint->GetX()>1)){
3710 clindex = track->GetClusterIndex2(i);
3712 x2[0] = trpoint->GetX();
3713 x2[1] = trpoint->GetY();
3714 x2[2] = trpoint->GetZ();
3715 sec2 = ((clindex&0xff000000)>>24)%18;
3722 Double_t alpha, cs,sn, xx2,yy2;
3724 alpha = (sec1-sec2)*fSectors->GetAlpha();
3725 cs = TMath::Cos(alpha);
3726 sn = TMath::Sin(alpha);
3727 xx2= x1[0]*cs-x1[1]*sn;
3728 yy2= x1[0]*sn+x1[1]*cs;
3732 alpha = (sec0-sec2)*fSectors->GetAlpha();
3733 cs = TMath::Cos(alpha);
3734 sn = TMath::Sin(alpha);
3735 xx2= x0[0]*cs-x0[1]*sn;
3736 yy2= x0[0]*sn+x0[1]*cs;
3742 Double_t x[5],c[15];
3746 x[4]=F1(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]);
3747 // if (x[4]>1) return 0;
3748 x[2]=F2(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]);
3749 x[3]=F3n(x2[0],x2[1],x0[0],x0[1],x2[2],x0[2],x[4]);
3750 //if (TMath::Abs(x[3]) > 2.2) return 0;
3751 //if (TMath::Abs(x[2]) > 1.99) return 0;
3753 Double_t sy =0.1, sz =0.1;
3755 Double_t sy1=0.02+track->GetSigmaY2(), sz1=0.02+track->GetSigmaZ2();
3756 Double_t sy2=0.01+track->GetSigmaY2(), sz2=0.01+track->GetSigmaZ2();
3757 Double_t sy3=0.01+track->GetSigmaY2();
3759 Double_t f40=(F1(x2[0],x2[1]+sy,x1[0],x1[1],x0[0],x0[1])-x[4])/sy;
3760 Double_t f42=(F1(x2[0],x2[1],x1[0],x1[1]+sy,x0[0],x0[1])-x[4])/sy;
3761 Double_t f43=(F1(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]+sy)-x[4])/sy;
3762 Double_t f20=(F2(x2[0],x2[1]+sy,x1[0],x1[1],x0[0],x0[1])-x[2])/sy;
3763 Double_t f22=(F2(x2[0],x2[1],x1[0],x1[1]+sy,x0[0],x0[1])-x[2])/sy;
3764 Double_t f23=(F2(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]+sy)-x[2])/sy;
3766 Double_t f30=(F3(x2[0],x2[1]+sy,x0[0],x0[1],x2[2],x0[2])-x[3])/sy;
3767 Double_t f31=(F3(x2[0],x2[1],x0[0],x0[1],x2[2]+sz,x0[2])-x[3])/sz;
3768 Double_t f32=(F3(x2[0],x2[1],x0[0],x0[1]+sy,x2[2],x0[2])-x[3])/sy;
3769 Double_t f34=(F3(x2[0],x2[1],x0[0],x0[1],x2[2],x0[2]+sz)-x[3])/sz;
3774 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3775 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3776 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3777 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3778 c[13]=f30*sy1*f40+f32*sy2*f42;
3779 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3781 // Int_t row1 = fSectors->GetRowNumber(x2[0]);
3782 AliTPCseed *seed=new AliTPCseed(x2[0], sec2*fSectors->GetAlpha()+fSectors->GetAlphaShift(), x, c, 0);
3783 // Double_t y0,z0,y1,z1, y2,z2;
3784 //seed->GetProlongation(x0[0],y0,z0);
3785 // seed->GetProlongation(x1[0],y1,z1);
3786 //seed->GetProlongation(x2[0],y2,z2);
3788 seed->SetLastPoint(pp2);
3789 seed->SetFirstPoint(pp2);
3796 AliTPCseed *AliTPCtrackerMI::ReSeed(const AliTPCseed *track, Float_t r0, Float_t r1, Float_t r2)
3800 //reseed using founded clusters
3802 // Find the number of clusters
3803 Int_t nclusters = 0;
3804 for (Int_t irow=0;irow<160;irow++){
3805 if (track->GetClusterIndex(irow)>0) nclusters++;
3809 ipos[0] = TMath::Max(int(r0*nclusters),0); // point 0 cluster
3810 ipos[1] = TMath::Min(int(r1*nclusters),nclusters-1); //
3811 ipos[2] = TMath::Min(int(r2*nclusters),nclusters-1); // last point
3815 Int_t row[3],sec[3]={0,0,0};
3817 // find track row position at given ratio of the length
3819 for (Int_t irow=0;irow<160;irow++){
3820 if (track->GetClusterIndex2(irow)<0) continue;
3822 for (Int_t ipoint=0;ipoint<3;ipoint++){
3823 if (index<=ipos[ipoint]) row[ipoint] = irow;
3827 //Get cluster and sector position
3828 for (Int_t ipoint=0;ipoint<3;ipoint++){
3829 Int_t clindex = track->GetClusterIndex2(row[ipoint]);
3830 AliTPCclusterMI * cl = GetClusterMI(clindex);
3833 // AliTPCclusterMI * cl = GetClusterMI(clindex);
3836 sec[ipoint] = ((clindex&0xff000000)>>24)%18;
3837 xyz[ipoint][0] = GetXrow(row[ipoint]);
3838 xyz[ipoint][1] = cl->GetY();
3839 xyz[ipoint][2] = cl->GetZ();
3843 // Calculate seed state vector and covariance matrix
3845 Double_t alpha, cs,sn, xx2,yy2;
3847 alpha = (sec[1]-sec[2])*fSectors->GetAlpha();
3848 cs = TMath::Cos(alpha);
3849 sn = TMath::Sin(alpha);
3850 xx2= xyz[1][0]*cs-xyz[1][1]*sn;
3851 yy2= xyz[1][0]*sn+xyz[1][1]*cs;
3855 alpha = (sec[0]-sec[2])*fSectors->GetAlpha();
3856 cs = TMath::Cos(alpha);
3857 sn = TMath::Sin(alpha);
3858 xx2= xyz[0][0]*cs-xyz[0][1]*sn;
3859 yy2= xyz[0][0]*sn+xyz[0][1]*cs;
3865 Double_t x[5],c[15];
3869 x[4]=F1(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
3870 x[2]=F2(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
3871 x[3]=F3n(xyz[2][0],xyz[2][1],xyz[0][0],xyz[0][1],xyz[2][2],xyz[0][2],x[4]);
3873 Double_t sy =0.1, sz =0.1;
3875 Double_t sy1=0.2, sz1=0.2;
3876 Double_t sy2=0.2, sz2=0.2;
3879 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;
3880 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;
3881 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;
3882 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;
3883 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;
3884 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;
3886 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;
3887 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;
3888 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;
3889 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;
3894 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3895 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3896 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3897 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3898 c[13]=f30*sy1*f40+f32*sy2*f42;
3899 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3901 // Int_t row1 = fSectors->GetRowNumber(xyz[2][0]);
3902 AliTPCseed *seed=new AliTPCseed(xyz[2][0], sec[2]*fSectors->GetAlpha()+fSectors->GetAlphaShift(), x, c, 0);
3903 seed->SetLastPoint(row[2]);
3904 seed->SetFirstPoint(row[2]);
3909 AliTPCseed *AliTPCtrackerMI::ReSeed(AliTPCseed *track,Int_t r0, Bool_t forward)
3913 //reseed using founded clusters
3916 Int_t row[3]={0,0,0};
3917 Int_t sec[3]={0,0,0};
3919 // forward direction
3921 for (Int_t irow=r0;irow<160;irow++){
3922 if (track->GetClusterIndex(irow)>0){
3927 for (Int_t irow=160;irow>r0;irow--){
3928 if (track->GetClusterIndex(irow)>0){
3933 for (Int_t irow=row[2]-15;irow>row[0];irow--){
3934 if (track->GetClusterIndex(irow)>0){
3942 for (Int_t irow=0;irow<r0;irow++){
3943 if (track->GetClusterIndex(irow)>0){
3948 for (Int_t irow=r0;irow>0;irow--){
3949 if (track->GetClusterIndex(irow)>0){
3954 for (Int_t irow=row[2]-15;irow>row[0];irow--){
3955 if (track->GetClusterIndex(irow)>0){
3962 if ((row[2]-row[0])<20) return 0;
3963 if (row[1]==0) return 0;
3966 //Get cluster and sector position
3967 for (Int_t ipoint=0;ipoint<3;ipoint++){
3968 Int_t clindex = track->GetClusterIndex2(row[ipoint]);
3969 AliTPCclusterMI * cl = GetClusterMI(clindex);
3972 // AliTPCclusterMI * cl = GetClusterMI(clindex);
3975 sec[ipoint] = ((clindex&0xff000000)>>24)%18;
3976 xyz[ipoint][0] = GetXrow(row[ipoint]);
3977 AliTPCTrackerPoint * point = track->GetTrackPoint(row[ipoint]);
3978 if (point&&ipoint<2){
3980 xyz[ipoint][1] = point->GetY();
3981 xyz[ipoint][2] = point->GetZ();
3984 xyz[ipoint][1] = cl->GetY();
3985 xyz[ipoint][2] = cl->GetZ();
3992 // Calculate seed state vector and covariance matrix
3994 Double_t alpha, cs,sn, xx2,yy2;
3996 alpha = (sec[1]-sec[2])*fSectors->GetAlpha();
3997 cs = TMath::Cos(alpha);
3998 sn = TMath::Sin(alpha);
3999 xx2= xyz[1][0]*cs-xyz[1][1]*sn;
4000 yy2= xyz[1][0]*sn+xyz[1][1]*cs;
4004 alpha = (sec[0]-sec[2])*fSectors->GetAlpha();
4005 cs = TMath::Cos(alpha);
4006 sn = TMath::Sin(alpha);
4007 xx2= xyz[0][0]*cs-xyz[0][1]*sn;
4008 yy2= xyz[0][0]*sn+xyz[0][1]*cs;
4014 Double_t x[5],c[15];
4018 x[4]=F1(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
4019 x[2]=F2(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
4020 x[3]=F3n(xyz[2][0],xyz[2][1],xyz[0][0],xyz[0][1],xyz[2][2],xyz[0][2],x[4]);
4022 Double_t sy =0.1, sz =0.1;
4024 Double_t sy1=0.2, sz1=0.2;
4025 Double_t sy2=0.2, sz2=0.2;
4028 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;
4029 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;
4030 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;
4031 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;
4032 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;
4033 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;
4035 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;
4036 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;
4037 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;
4038 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;
4043 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
4044 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
4045 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
4046 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
4047 c[13]=f30*sy1*f40+f32*sy2*f42;
4048 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
4050 // Int_t row1 = fSectors->GetRowNumber(xyz[2][0]);
4051 AliTPCseed *seed=new AliTPCseed(xyz[2][0], sec[2]*fSectors->GetAlpha()+fSectors->GetAlphaShift(), x, c, 0);
4052 seed->SetLastPoint(row[2]);
4053 seed->SetFirstPoint(row[2]);
4054 for (Int_t i=row[0];i<row[2];i++){
4055 seed->SetClusterIndex(i, track->GetClusterIndex(i));
4063 void AliTPCtrackerMI::FindMultiMC(const TObjArray * array, AliESDEvent */*esd*/, Int_t iter)
4066 // find multi tracks - THIS FUNCTION IS ONLY FOR DEBUG PURPOSES
4068 // Use AliTPCReconstructor::StreamLevel()>2 if you want to tune parameters - cuts
4070 // Two reasons to have multiple find tracks
4071 // 1. Curling tracks can be find more than once
4072 // 2. Splitted tracks
4073 // a.) Multiple seeding to increase tracking efficiency - (~ 100% reached)
4074 // b.) Edge effect on the sector boundaries
4077 // Algorithm done in 2 phases - because of CPU consumption
4078 // it is n^2 algorithm - for lead-lead 20000x20000 combination are investigated
4080 // Algorihm for curling tracks sign:
4081 // 1 phase -makes a very rough fast cuts to minimize combinatorics
4082 // a.) opposite sign
4083 // b.) one of the tracks - not pointing to the primary vertex -
4084 // c.) delta tan(theta)
4086 // 2 phase - calculates DCA between tracks - time consument
4091 // General cuts - for splitted tracks and for curling tracks
4093 const Float_t kMaxdPhi = 0.2; // maximal distance in phi
4095 // Curling tracks cuts
4100 Int_t nentries = array->GetEntriesFast();
4101 AliHelix *helixes = new AliHelix[nentries];
4102 Float_t *xm = new Float_t[nentries];
4103 Float_t *dz0 = new Float_t[nentries];
4104 Float_t *dz1 = new Float_t[nentries];
4110 // Find track COG in x direction - point with best defined parameters
4112 for (Int_t i=0;i<nentries;i++){
4113 AliTPCseed* track = (AliTPCseed*)array->At(i);
4114 if (!track) continue;
4115 track->SetCircular(0);
4116 new (&helixes[i]) AliHelix(*track);
4120 track->GetDZ(GetX(),GetY(),GetZ(),GetBz(),dz);
4123 for (Int_t icl=0; icl<160; icl++){
4124 AliTPCclusterMI * cl = track->GetClusterPointer(icl);
4130 if (ncl>0) xm[i]/=Float_t(ncl);
4133 for (Int_t i0=0;i0<nentries;i0++){
4134 AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
4135 if (!track0) continue;
4136 Float_t xc0 = helixes[i0].GetHelix(6);
4137 Float_t yc0 = helixes[i0].GetHelix(7);
4138 Float_t r0 = helixes[i0].GetHelix(8);
4139 Float_t rc0 = TMath::Sqrt(xc0*xc0+yc0*yc0);
4140 Float_t fi0 = TMath::ATan2(yc0,xc0);
4142 for (Int_t i1=i0+1;i1<nentries;i1++){
4143 AliTPCseed * track1 = (AliTPCseed*)array->At(i1);
4144 if (!track1) continue;
4145 Int_t lab0=track0->GetLabel();
4146 Int_t lab1=track1->GetLabel();
4147 if (TMath::Abs(lab0)!=TMath::Abs(lab1)) continue;
4149 Float_t xc1 = helixes[i1].GetHelix(6);
4150 Float_t yc1 = helixes[i1].GetHelix(7);
4151 Float_t r1 = helixes[i1].GetHelix(8);
4152 Float_t rc1 = TMath::Sqrt(xc1*xc1+yc1*yc1);
4153 Float_t fi1 = TMath::ATan2(yc1,xc1);
4155 Float_t dfi = fi0-fi1;
4158 if (dfi>1.5*TMath::Pi()) dfi-=TMath::Pi(); // take care about edge effect
4159 if (dfi<-1.5*TMath::Pi()) dfi+=TMath::Pi(); //
4160 if (TMath::Abs(dfi)>kMaxdPhi&&helixes[i0].GetHelix(4)*helixes[i1].GetHelix(4)<0){
4162 // if short tracks with undefined sign
4163 fi1 = -TMath::ATan2(yc1,-xc1);
4166 Float_t dtheta = TMath::Abs(track0->GetTgl()-track1->GetTgl())<TMath::Abs(track0->GetTgl()+track1->GetTgl())? track0->GetTgl()-track1->GetTgl():track0->GetTgl()+track1->GetTgl();
4169 // debug stream to tune "fast cuts"
4171 Double_t dist[3]; // distance at X
4172 Double_t mdist[3]={0,0,0}; // mean distance X+-40cm
4173 track0->GetDistance(track1,0.5*(xm[i0]+xm[i1])-40.,dist,AliTracker::GetBz());
4174 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[i]);
4175 track0->GetDistance(track1,0.5*(xm[i0]+xm[i1])+40.,dist,AliTracker::GetBz());
4176 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[i]);
4177 track0->GetDistance(track1,0.5*(xm[i0]+xm[i1]),dist,AliTracker::GetBz());
4178 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[i]);
4179 for (Int_t i=0;i<3;i++) mdist[i]*=0.33333;
4183 for (Int_t icl=0; icl<160; icl++){
4184 AliTPCclusterMI * cl0 = track0->GetClusterPointer(icl);
4185 AliTPCclusterMI * cl1 = track1->GetClusterPointer(icl);
4188 if (cl0==cl1) sums++;
4192 if (AliTPCReconstructor::StreamLevel()>0) {
4193 TTreeSRedirector &cstream = *fDebugStreamer;
4198 "Tr0.="<<track0<< // seed0
4199 "Tr1.="<<track1<< // seed1
4200 "h0.="<<&helixes[i0]<<
4201 "h1.="<<&helixes[i1]<<
4203 "sum="<<sum<< //the sum of rows with cl in both
4204 "sums="<<sums<< //the sum of shared clusters
4205 "xm0="<<xm[i0]<< // the center of track
4206 "xm1="<<xm[i1]<< // the x center of track
4207 // General cut variables
4208 "dfi="<<dfi<< // distance in fi angle
4209 "dtheta="<<dtheta<< // distance int theta angle
4215 "dist0="<<dist[0]<< //distance x
4216 "dist1="<<dist[1]<< //distance y
4217 "dist2="<<dist[2]<< //distance z
4218 "mdist0="<<mdist[0]<< //distance x
4219 "mdist1="<<mdist[1]<< //distance y
4220 "mdist2="<<mdist[2]<< //distance z
4234 if (AliTPCReconstructor::StreamLevel()>1) {
4235 AliInfo("Time for curling tracks removal DEBUGGING MC");
4241 void AliTPCtrackerMI::FindSplitted(TObjArray * array, AliESDEvent */*esd*/, Int_t iter)
4245 // Two reasons to have multiple find tracks
4246 // 1. Curling tracks can be find more than once
4247 // 2. Splitted tracks
4248 // a.) Multiple seeding to increase tracking efficiency - (~ 100% reached)
4249 // b.) Edge effect on the sector boundaries
4251 // This function tries to find tracks closed in the parametric space
4253 // cut logic if distance is bigger than cut continue - Do Nothing
4254 const Float_t kMaxdTheta = 0.05; // maximal distance in theta
4255 const Float_t kMaxdPhi = 0.05; // maximal deistance in phi
4256 const Float_t kdelta = 40.; //delta r to calculate track distance
4258 // const Float_t kMaxDist0 = 1.; // maximal distance 0
4259 //const Float_t kMaxDist1 = 0.3; // maximal distance 1 - cut if track in separate rows
4262 TCut csec("csec","abs(Tr0.fRelativeSector-Tr1.fRelativeSector)<2");
4263 TCut cdtheta("cdtheta","abs(dtheta)<0.05");
4268 Int_t nentries = array->GetEntriesFast();
4269 AliHelix *helixes = new AliHelix[nentries];
4270 Float_t *xm = new Float_t[nentries];
4276 //Sort tracks according quality
4278 Int_t nseed = array->GetEntriesFast();
4279 Float_t * quality = new Float_t[nseed];
4280 Int_t * indexes = new Int_t[nseed];
4281 for (Int_t i=0; i<nseed; i++) {
4282 AliTPCseed *pt=(AliTPCseed*)array->UncheckedAt(i);
4287 pt->UpdatePoints(); //select first last max dens points
4288 Float_t * points = pt->GetPoints();
4289 if (points[3]<0.8) quality[i] =-1;
4290 quality[i] = (points[2]-points[0])+pt->GetNumberOfClusters();
4291 //prefer high momenta tracks if overlaps
4292 quality[i] *= TMath::Sqrt(TMath::Abs(pt->Pt())+0.5);
4294 TMath::Sort(nseed,quality,indexes);
4298 // Find track COG in x direction - point with best defined parameters
4300 for (Int_t i=0;i<nentries;i++){
4301 AliTPCseed* track = (AliTPCseed*)array->At(i);
4302 if (!track) continue;
4303 track->SetCircular(0);
4304 new (&helixes[i]) AliHelix(*track);
4307 for (Int_t icl=0; icl<160; icl++){
4308 AliTPCclusterMI * cl = track->GetClusterPointer(icl);
4314 if (ncl>0) xm[i]/=Float_t(ncl);
4317 for (Int_t is0=0;is0<nentries;is0++){
4318 Int_t i0 = indexes[is0];
4319 AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
4320 if (!track0) continue;
4321 Float_t xc0 = helixes[i0].GetHelix(6);
4322 Float_t yc0 = helixes[i0].GetHelix(7);
4323 Float_t fi0 = TMath::ATan2(yc0,xc0);
4325 for (Int_t is1=is0+1;is1<nentries;is1++){
4326 Int_t i1 = indexes[is1];
4327 AliTPCseed * track1 = (AliTPCseed*)array->At(i1);
4328 if (!track1) continue;
4330 Int_t dsec = TMath::Abs((track0->GetRelativeSector()%18)-(track1->GetRelativeSector()%18)); // sector distance
4331 if (dsec>1 && dsec<17) continue;
4333 if (track1->GetKinkIndexes()[0] == -track0->GetKinkIndexes()[0]) continue;
4335 Float_t dtheta = TMath::Abs(track0->GetTgl()-track1->GetTgl())<TMath::Abs(track0->GetTgl()+track1->GetTgl())? track0->GetTgl()-track1->GetTgl():track0->GetTgl()+track1->GetTgl();
4336 if (TMath::Abs(dtheta)>kMaxdTheta) continue;
4338 Float_t xc1 = helixes[i1].GetHelix(6);
4339 Float_t yc1 = helixes[i1].GetHelix(7);
4340 Float_t fi1 = TMath::ATan2(yc1,xc1);
4342 Float_t dfi = fi0-fi1;
4343 if (dfi>TMath::Pi()) dfi-=TMath::TwoPi(); // take care about edge effect
4344 if (dfi<-TMath::Pi()) dfi+=TMath::TwoPi(); //
4345 if (TMath::Abs(dfi)>kMaxdPhi&&helixes[i0].GetHelix(4)*helixes[i1].GetHelix(4)<0){
4347 // if short tracks with undefined sign
4348 fi1 = -TMath::ATan2(yc1,-xc1);
4350 if (dfi>TMath::Pi()) dfi-=TMath::TwoPi(); // take care about edge effect
4351 if (dfi<-TMath::Pi()) dfi+=TMath::TwoPi(); //
4353 if (TMath::Abs(dfi)>kMaxdPhi) continue;
4360 for (Int_t icl=0; icl<160; icl++){
4361 Int_t index0=track0->GetClusterIndex2(icl);
4362 Int_t index1=track1->GetClusterIndex2(icl);
4363 Bool_t used0 = (index0>0 && !(index0&0x8000));
4364 Bool_t used1 = (index1>0 && !(index1&0x8000));
4366 if (used0) sum0++; // used cluster0
4367 if (used1) sum1++; // used clusters1
4368 if (used0&&used1) sum++;
4369 if (index0==index1 && used0 && used1) sums++;
4373 if (sums<10) continue;
4374 if (sum<40) continue;
4375 if (sums/Float_t(TMath::Min(sum0,sum1))<0.5) continue;
4377 Double_t dist[5][4]; // distance at X
4378 Double_t mdist[4]={0,0,0,0}; // mean distance on range +-delta
4382 track0->GetDistance(track1,xm[i0],dist[0],AliTracker::GetBz());
4383 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[0][i]);
4384 track0->GetDistance(track1,xm[i1],dist[1],AliTracker::GetBz());
4385 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[1][i]);
4387 track0->GetDistance(track1,TMath::Min(xm[i1],xm[i0])-kdelta,dist[2],AliTracker::GetBz());
4388 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[2][i]);
4389 track0->GetDistance(track1,TMath::Max(xm[i1],xm[i0])+kdelta,dist[3],AliTracker::GetBz());
4390 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[3][i]);
4392 track0->GetDistance(track1,(xm[i1]+xm[i0])*0.5,dist[4],AliTracker::GetBz());
4393 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[4][i]);
4394 for (Int_t i=0;i<3;i++) mdist[i]*=0.2;
4397 Int_t lab0=track0->GetLabel();
4398 Int_t lab1=track1->GetLabel();
4399 if( AliTPCReconstructor::StreamLevel()>5){
4400 TTreeSRedirector &cstream = *fDebugStreamer;
4401 cstream<<"Splitted"<<
4405 "Tr0.="<<track0<< // seed0
4406 "Tr1.="<<track1<< // seed1
4407 "h0.="<<&helixes[i0]<<
4408 "h1.="<<&helixes[i1]<<
4410 "sum="<<sum<< //the sum of rows with cl in both
4411 "sum0="<<sum0<< //the sum of rows with cl in 0 track
4412 "sum1="<<sum1<< //the sum of rows with cl in 1 track
4413 "sums="<<sums<< //the sum of shared clusters
4414 "xm0="<<xm[i0]<< // the center of track
4415 "xm1="<<xm[i1]<< // the x center of track
4416 // General cut variables
4417 "dfi="<<dfi<< // distance in fi angle
4418 "dtheta="<<dtheta<< // distance int theta angle
4421 "dist0="<<dist[4][0]<< //distance x
4422 "dist1="<<dist[4][1]<< //distance y
4423 "dist2="<<dist[4][2]<< //distance z
4424 "mdist0="<<mdist[0]<< //distance x
4425 "mdist1="<<mdist[1]<< //distance y
4426 "mdist2="<<mdist[2]<< //distance z
4429 delete array->RemoveAt(i1);
4436 AliInfo("Time for splitted tracks removal");
4442 void AliTPCtrackerMI::FindCurling(const TObjArray * array, AliESDEvent */*esd*/, Int_t iter)
4445 // find Curling tracks
4446 // Use AliTPCReconstructor::StreamLevel()>1 if you want to tune parameters - cuts
4449 // Algorithm done in 2 phases - because of CPU consumption
4450 // it is n^2 algorithm - for lead-lead 20000x20000 combination are investigated
4451 // see detal in MC part what can be used to cut
4455 const Float_t kMaxC = 400; // maximal curvature to of the track
4456 const Float_t kMaxdTheta = 0.15; // maximal distance in theta
4457 const Float_t kMaxdPhi = 0.15; // maximal distance in phi
4458 const Float_t kPtRatio = 0.3; // ratio between pt
4459 const Float_t kMinDCAR = 2.; // distance to the primary vertex in r - see cpipe cut
4462 // Curling tracks cuts
4465 const Float_t kMaxDeltaRMax = 40; // distance in outer radius
4466 const Float_t kMaxDeltaRMin = 5.; // distance in lower radius - see cpipe cut
4467 const Float_t kMinAngle = 2.9; // angle between tracks
4468 const Float_t kMaxDist = 5; // biggest distance
4470 // The cuts can be tuned using the "MC information stored in Multi tree ==> see FindMultiMC
4473 TCut csign("csign","Tr0.fP[4]*Tr1.fP[4]<0"); //opposite sign
4474 TCut cmax("cmax","abs(Tr0.GetC())>1/400");
4475 TCut cda("cda","sqrt(dtheta^2+dfi^2)<0.15");
4476 TCut ccratio("ccratio","abs((Tr0.fP[4]+Tr1.fP[4])/(abs(Tr0.fP[4])+abs(Tr1.fP[4])))<0.3");
4477 TCut cpipe("cpipe", "min(abs(r0-rc0),abs(r1-rc1))>5");
4479 TCut cdrmax("cdrmax","abs(abs(rc0+r0)-abs(rc1+r1))<40")
4480 TCut cdrmin("cdrmin","abs(abs(rc0+r0)-abs(rc1+r1))<10")
4482 Multi->Draw("dfi","iter==0"+csign+cmax+cda+ccratio); ~94% of curling tracks fulfill
4483 Multi->Draw("min(abs(r0-rc0),abs(r1-rc1))","iter==0&&abs(lab1)==abs(lab0)"+csign+cmax+cda+ccratio+cpipe+cdrmin+cdrmax); //80%
4485 Curling2->Draw("dfi","iter==0&&abs(lab0)==abs(lab1)"+csign+cmax+cdtheta+cdfi+ccratio)
4491 Int_t nentries = array->GetEntriesFast();
4492 AliHelix *helixes = new AliHelix[nentries];
4493 for (Int_t i=0;i<nentries;i++){
4494 AliTPCseed* track = (AliTPCseed*)array->At(i);
4495 if (!track) continue;
4496 track->SetCircular(0);
4497 new (&helixes[i]) AliHelix(*track);
4503 Double_t phase[2][2],radius[2];
4508 for (Int_t i0=0;i0<nentries;i0++){
4509 AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
4510 if (!track0) continue;
4511 if (TMath::Abs(track0->GetC())<1/kMaxC) continue;
4512 Float_t xc0 = helixes[i0].GetHelix(6);
4513 Float_t yc0 = helixes[i0].GetHelix(7);
4514 Float_t r0 = helixes[i0].GetHelix(8);
4515 Float_t rc0 = TMath::Sqrt(xc0*xc0+yc0*yc0);
4516 Float_t fi0 = TMath::ATan2(yc0,xc0);
4518 for (Int_t i1=i0+1;i1<nentries;i1++){
4519 AliTPCseed * track1 = (AliTPCseed*)array->At(i1);
4520 if (!track1) continue;
4521 if (TMath::Abs(track1->GetC())<1/kMaxC) continue;
4522 Float_t xc1 = helixes[i1].GetHelix(6);
4523 Float_t yc1 = helixes[i1].GetHelix(7);
4524 Float_t r1 = helixes[i1].GetHelix(8);
4525 Float_t rc1 = TMath::Sqrt(xc1*xc1+yc1*yc1);
4526 Float_t fi1 = TMath::ATan2(yc1,xc1);
4528 Float_t dfi = fi0-fi1;
4531 if (dfi>1.5*TMath::Pi()) dfi-=TMath::Pi(); // take care about edge effect
4532 if (dfi<-1.5*TMath::Pi()) dfi+=TMath::Pi(); //
4533 Float_t dtheta = TMath::Abs(track0->GetTgl()-track1->GetTgl())<TMath::Abs(track0->GetTgl()+track1->GetTgl())? track0->GetTgl()-track1->GetTgl():track0->GetTgl()+track1->GetTgl();
4537 if (track0->GetBConstrain()&&track1->GetBConstrain()) continue; // not constrained
4538 if (track1->GetSigned1Pt()*track0->GetSigned1Pt()>0) continue; // not the same sign
4539 if ( TMath::Abs(track1->GetTgl()+track0->GetTgl())>kMaxdTheta) continue; //distance in the Theta
4540 if ( TMath::Abs(dfi)>kMaxdPhi) continue; //distance in phi
4541 if ( TMath::Sqrt(dfi*dfi+dtheta*dtheta)>kMaxdPhi) continue; //common angular offset
4543 Float_t pt0 = track0->GetSignedPt();
4544 Float_t pt1 = track1->GetSignedPt();
4545 if ((TMath::Abs(pt0+pt1)/(TMath::Abs(pt0)+TMath::Abs(pt1)))>kPtRatio) continue;
4546 if ((iter==1) && TMath::Abs(TMath::Abs(rc0+r0)-TMath::Abs(rc1+r1))>kMaxDeltaRMax) continue;
4547 if ((iter!=1) &&TMath::Abs(TMath::Abs(rc0-r0)-TMath::Abs(rc1-r1))>kMaxDeltaRMin) continue;
4548 if (TMath::Min(TMath::Abs(rc0-r0),TMath::Abs(rc1-r1))<kMinDCAR) continue;
4551 // Now find closest approach
4555 Int_t npoints = helixes[i0].GetRPHIintersections(helixes[i1], phase, radius,10);
4556 if (npoints==0) continue;
4557 helixes[i0].GetClosestPhases(helixes[i1], phase);
4561 Double_t hangles[3];
4562 helixes[i0].Evaluate(phase[0][0],xyz0);
4563 helixes[i1].Evaluate(phase[0][1],xyz1);
4565 helixes[i0].GetAngle(phase[0][0],helixes[i1],phase[0][1],hangles);
4566 Double_t deltah[2],deltabest;
4567 if (TMath::Abs(hangles[2])<kMinAngle) continue;
4571 helixes[i0].ParabolicDCA(helixes[i1],phase[0][0],phase[0][1],radius[0],deltah[0],2);
4573 helixes[i0].ParabolicDCA(helixes[i1],phase[1][0],phase[1][1],radius[1],deltah[1],2);
4574 if (deltah[1]<deltah[0]) ibest=1;
4576 deltabest = TMath::Sqrt(deltah[ibest]);
4577 helixes[i0].Evaluate(phase[ibest][0],xyz0);
4578 helixes[i1].Evaluate(phase[ibest][1],xyz1);
4579 helixes[i0].GetAngle(phase[ibest][0],helixes[i1],phase[ibest][1],hangles);
4580 Double_t radiusbest = TMath::Sqrt(radius[ibest]);
4582 if (deltabest>kMaxDist) continue;
4583 // if (mindcar+mindcaz<40 && (TMath::Abs(hangles[2])<kMinAngle ||deltabest>3)) continue;
4584 Bool_t sign =kFALSE;
4585 if (hangles[2]>kMinAngle) sign =kTRUE;
4588 // circular[i0] = kTRUE;
4589 // circular[i1] = kTRUE;
4590 if (track0->OneOverPt()<track1->OneOverPt()){
4591 track0->SetCircular(track0->GetCircular()+1);
4592 track1->SetCircular(track1->GetCircular()+2);
4595 track1->SetCircular(track1->GetCircular()+1);
4596 track0->SetCircular(track0->GetCircular()+2);
4599 if (AliTPCReconstructor::StreamLevel()>1){
4601 //debug stream to tune "fine" cuts
4602 Int_t lab0=track0->GetLabel();
4603 Int_t lab1=track1->GetLabel();
4604 TTreeSRedirector &cstream = *fDebugStreamer;
4605 cstream<<"Curling2"<<
4621 "npoints="<<npoints<<
4622 "hangles0="<<hangles[0]<<
4623 "hangles1="<<hangles[1]<<
4624 "hangles2="<<hangles[2]<<
4627 "radius="<<radiusbest<<
4628 "deltabest="<<deltabest<<
4629 "phase0="<<phase[ibest][0]<<
4630 "phase1="<<phase[ibest][1]<<
4638 if (AliTPCReconstructor::StreamLevel()>1) {
4639 AliInfo("Time for curling tracks removal");
4648 void AliTPCtrackerMI::FindKinks(TObjArray * array, AliESDEvent *esd)
4655 TObjArray *kinks= new TObjArray(10000);
4656 // TObjArray *v0s= new TObjArray(10000);
4657 Int_t nentries = array->GetEntriesFast();
4658 AliHelix *helixes = new AliHelix[nentries];
4659 Int_t *sign = new Int_t[nentries];
4660 Int_t *nclusters = new Int_t[nentries];
4661 Float_t *alpha = new Float_t[nentries];
4662 AliKink *kink = new AliKink();
4663 Int_t * usage = new Int_t[nentries];
4664 Float_t *zm = new Float_t[nentries];
4665 Float_t *z0 = new Float_t[nentries];
4666 Float_t *fim = new Float_t[nentries];
4667 Float_t *shared = new Float_t[nentries];
4668 Bool_t *circular = new Bool_t[nentries];
4669 Float_t *dca = new Float_t[nentries];
4670 //const AliESDVertex * primvertex = esd->GetVertex();
4672 // nentries = array->GetEntriesFast();
4677 for (Int_t i=0;i<nentries;i++){
4680 AliTPCseed* track = (AliTPCseed*)array->At(i);
4681 if (!track) continue;
4682 track->SetCircular(0);
4684 track->UpdatePoints();
4685 if (( track->GetPoints()[2]- track->GetPoints()[0])>5 && track->GetPoints()[3]>0.8){
4687 nclusters[i]=track->GetNumberOfClusters();
4688 alpha[i] = track->GetAlpha();
4689 new (&helixes[i]) AliHelix(*track);
4691 helixes[i].Evaluate(0,xyz);
4692 sign[i] = (track->GetC()>0) ? -1:1;
4695 if (track->GetProlongation(x,y,z)){
4697 fim[i] = alpha[i]+TMath::ATan2(y,x);
4700 zm[i] = track->GetZ();
4704 circular[i]= kFALSE;
4705 if (track->GetProlongation(0,y,z)) z0[i] = z;
4706 dca[i] = track->GetD(0,0);
4712 Int_t ncandidates =0;
4715 Double_t phase[2][2],radius[2];
4718 // Find circling track
4720 for (Int_t i0=0;i0<nentries;i0++){
4721 AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
4722 if (!track0) continue;
4723 if (track0->GetNumberOfClusters()<40) continue;
4724 if (TMath::Abs(1./track0->GetC())>200) continue;
4725 for (Int_t i1=i0+1;i1<nentries;i1++){
4726 AliTPCseed * track1 = (AliTPCseed*)array->At(i1);
4727 if (!track1) continue;
4728 if (track1->GetNumberOfClusters()<40) continue;
4729 if ( TMath::Abs(track1->GetTgl()+track0->GetTgl())>0.1) continue;
4730 if (track0->GetBConstrain()&&track1->GetBConstrain()) continue;
4731 if (TMath::Abs(1./track1->GetC())>200) continue;
4732 if (track1->GetSigned1Pt()*track0->GetSigned1Pt()>0) continue;
4733 if (track1->GetTgl()*track0->GetTgl()>0) continue;
4734 if (TMath::Max(TMath::Abs(1./track0->GetC()),TMath::Abs(1./track1->GetC()))>190) continue;
4735 if (track0->GetBConstrain()&&track1->OneOverPt()<track0->OneOverPt()) continue; //returning - lower momenta
4736 if (track1->GetBConstrain()&&track0->OneOverPt()<track1->OneOverPt()) continue; //returning - lower momenta
4738 Float_t mindcar = TMath::Min(TMath::Abs(dca[i0]),TMath::Abs(dca[i1]));
4739 if (mindcar<5) continue;
4740 Float_t mindcaz = TMath::Min(TMath::Abs(z0[i0]-GetZ()),TMath::Abs(z0[i1]-GetZ()));
4741 if (mindcaz<5) continue;
4742 if (mindcar+mindcaz<20) continue;
4745 Float_t xc0 = helixes[i0].GetHelix(6);
4746 Float_t yc0 = helixes[i0].GetHelix(7);
4747 Float_t r0 = helixes[i0].GetHelix(8);
4748 Float_t xc1 = helixes[i1].GetHelix(6);
4749 Float_t yc1 = helixes[i1].GetHelix(7);
4750 Float_t r1 = helixes[i1].GetHelix(8);
4752 Float_t rmean = (r0+r1)*0.5;
4753 Float_t delta =TMath::Sqrt((xc1-xc0)*(xc1-xc0)+(yc1-yc0)*(yc1-yc0));
4754 //if (delta>30) continue;
4755 if (delta>rmean*0.25) continue;
4756 if (TMath::Abs(r0-r1)/rmean>0.3) continue;
4758 Int_t npoints = helixes[i0].GetRPHIintersections(helixes[i1], phase, radius,10);
4759 if (npoints==0) continue;
4760 helixes[i0].GetClosestPhases(helixes[i1], phase);
4764 Double_t hangles[3];
4765 helixes[i0].Evaluate(phase[0][0],xyz0);
4766 helixes[i1].Evaluate(phase[0][1],xyz1);
4768 helixes[i0].GetAngle(phase[0][0],helixes[i1],phase[0][1],hangles);
4769 Double_t deltah[2],deltabest;
4770 if (hangles[2]<2.8) continue;
4773 helixes[i0].ParabolicDCA(helixes[i1],phase[0][0],phase[0][1],radius[0],deltah[0],2);
4775 helixes[i0].ParabolicDCA(helixes[i1],phase[1][0],phase[1][1],radius[1],deltah[1],2);
4776 if (deltah[1]<deltah[0]) ibest=1;
4778 deltabest = TMath::Sqrt(deltah[ibest]);
4779 helixes[i0].Evaluate(phase[ibest][0],xyz0);
4780 helixes[i1].Evaluate(phase[ibest][1],xyz1);
4781 helixes[i0].GetAngle(phase[ibest][0],helixes[i1],phase[ibest][1],hangles);
4782 Double_t radiusbest = TMath::Sqrt(radius[ibest]);
4784 if (deltabest>6) continue;
4785 if (mindcar+mindcaz<40 && (hangles[2]<3.12||deltabest>3)) continue;
4786 Bool_t lsign =kFALSE;
4787 if (hangles[2]>3.06) lsign =kTRUE;
4790 circular[i0] = kTRUE;
4791 circular[i1] = kTRUE;
4792 if (track0->OneOverPt()<track1->OneOverPt()){
4793 track0->SetCircular(track0->GetCircular()+1);
4794 track1->SetCircular(track1->GetCircular()+2);
4797 track1->SetCircular(track1->GetCircular()+1);
4798 track0->SetCircular(track0->GetCircular()+2);
4801 if (lsign&&AliTPCReconstructor::StreamLevel()>1){
4803 Int_t lab0=track0->GetLabel();
4804 Int_t lab1=track1->GetLabel();
4805 TTreeSRedirector &cstream = *fDebugStreamer;
4806 cstream<<"Curling"<<
4813 "mindcar="<<mindcar<<
4814 "mindcaz="<<mindcaz<<
4817 "npoints="<<npoints<<
4818 "hangles0="<<hangles[0]<<
4819 "hangles2="<<hangles[2]<<
4824 "radius="<<radiusbest<<
4825 "deltabest="<<deltabest<<
4826 "phase0="<<phase[ibest][0]<<
4827 "phase1="<<phase[ibest][1]<<
4837 for (Int_t i =0;i<nentries;i++){
4838 if (sign[i]==0) continue;
4839 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
4846 Double_t cradius0 = 40*40;
4847 Double_t cradius1 = 270*270;
4850 Double_t cdist3=0.55;
4851 for (Int_t j =i+1;j<nentries;j++){
4853 if (sign[j]*sign[i]<1) continue;
4854 if ( (nclusters[i]+nclusters[j])>200) continue;
4855 if ( (nclusters[i]+nclusters[j])<80) continue;
4856 if ( TMath::Abs(zm[i]-zm[j])>60.) continue;
4857 if ( TMath::Abs(fim[i]-fim[j])>0.6 && TMath::Abs(fim[i]-fim[j])<5.7 ) continue;
4858 //AliTPCseed * track1 = (AliTPCseed*)array->At(j); Double_t phase[2][2],radius[2];
4859 Int_t npoints = helixes[i].GetRPHIintersections(helixes[j], phase, radius,20);
4860 if (npoints<1) continue;
4863 if (radius[0]<cradius0||radius[0]>cradius1) continue;
4866 if ( (radius[0]<cradius0||radius[0]>cradius1) && (radius[1]<cradius0||radius[1]>cradius1) ) continue;
4869 Double_t delta1=10000,delta2=10000;
4870 // cuts on the intersection radius
4871 helixes[i].LinearDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta1);
4872 if (radius[0]<20&&delta1<1) continue; //intersection at vertex
4873 if (radius[0]<10&&delta1<3) continue; //intersection at vertex
4875 helixes[i].LinearDCA(helixes[j],phase[1][0],phase[1][1],radius[1],delta2);
4876 if (radius[1]<20&&delta2<1) continue; //intersection at vertex
4877 if (radius[1]<10&&delta2<3) continue; //intersection at vertex
4880 Double_t distance1 = TMath::Min(delta1,delta2);
4881 if (distance1>cdist1) continue; // cut on DCA linear approximation
4883 npoints = helixes[i].GetRPHIintersections(helixes[j], phase, radius,20);
4884 helixes[i].ParabolicDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta1);
4885 if (radius[0]<20&&delta1<1) continue; //intersection at vertex
4886 if (radius[0]<10&&delta1<3) continue; //intersection at vertex
4889 helixes[i].ParabolicDCA(helixes[j],phase[1][0],phase[1][1],radius[1],delta2);
4890 if (radius[1]<20&&delta2<1) continue; //intersection at vertex
4891 if (radius[1]<10&&delta2<3) continue; //intersection at vertex
4893 distance1 = TMath::Min(delta1,delta2);
4896 rkink = TMath::Sqrt(radius[0]);
4899 rkink = TMath::Sqrt(radius[1]);
4901 if (distance1>cdist2) continue;
4904 AliTPCseed * track1 = (AliTPCseed*)array->At(j);
4907 Int_t row0 = GetRowNumber(rkink);
4908 if (row0<10) continue;
4909 if (row0>150) continue;
4912 Float_t dens00=-1,dens01=-1;
4913 Float_t dens10=-1,dens11=-1;
4915 Int_t found,foundable,ishared;
4916 track0->GetClusterStatistic(0,row0-5, found, foundable,ishared,kFALSE);
4917 if (foundable>5) dens00 = Float_t(found)/Float_t(foundable);
4918 track0->GetClusterStatistic(row0+5,155, found, foundable,ishared,kFALSE);
4919 if (foundable>5) dens01 = Float_t(found)/Float_t(foundable);
4921 track1->GetClusterStatistic(0,row0-5, found, foundable,ishared,kFALSE);
4922 if (foundable>10) dens10 = Float_t(found)/Float_t(foundable);
4923 track1->GetClusterStatistic(row0+5,155, found, foundable,ishared,kFALSE);
4924 if (foundable>10) dens11 = Float_t(found)/Float_t(foundable);
4926 if (dens00<dens10 && dens01<dens11) continue;
4927 if (dens00>dens10 && dens01>dens11) continue;
4928 if (TMath::Max(dens00,dens10)<0.1) continue;
4929 if (TMath::Max(dens01,dens11)<0.3) continue;
4931 if (TMath::Min(dens00,dens10)>0.6) continue;
4932 if (TMath::Min(dens01,dens11)>0.6) continue;
4935 AliTPCseed * ktrack0, *ktrack1;
4944 if (TMath::Abs(ktrack0->GetC())>5) continue; // cut on the curvature for mother particle
4945 AliExternalTrackParam paramm(*ktrack0);
4946 AliExternalTrackParam paramd(*ktrack1);
4947 if (row0>60&&ktrack1->GetReference().GetX()>90.)new (¶md) AliExternalTrackParam(ktrack1->GetReference());
4950 kink->SetMother(paramm);
4951 kink->SetDaughter(paramd);
4954 Float_t x[3] = { kink->GetPosition()[0],kink->GetPosition()[1],kink->GetPosition()[2]};
4956 fkParam->Transform0to1(x,index);
4957 fkParam->Transform1to2(x,index);
4958 row0 = GetRowNumber(x[0]);
4960 if (kink->GetR()<100) continue;
4961 if (kink->GetR()>240) continue;
4962 if (kink->GetPosition()[2]/kink->GetR()>AliTPCReconstructor::GetCtgRange()) continue; //out of fiducial volume
4963 if (kink->GetDistance()>cdist3) continue;
4964 Float_t dird = kink->GetDaughterP()[0]*kink->GetPosition()[0]+kink->GetDaughterP()[1]*kink->GetPosition()[1]; // rough direction estimate
4965 if (dird<0) continue;
4967 Float_t dirm = kink->GetMotherP()[0]*kink->GetPosition()[0]+kink->GetMotherP()[1]*kink->GetPosition()[1]; // rough direction estimate
4968 if (dirm<0) continue;
4969 Float_t mpt = TMath::Sqrt(kink->GetMotherP()[0]*kink->GetMotherP()[0]+kink->GetMotherP()[1]*kink->GetMotherP()[1]);
4970 if (mpt<0.2) continue;
4973 //for high momenta momentum not defined well in first iteration
4974 Double_t qt = TMath::Sin(kink->GetAngle(2))*ktrack1->GetP();
4975 if (qt>0.35) continue;
4978 kink->SetLabel(CookLabel(ktrack0,0.4,0,row0),0);
4979 kink->SetLabel(CookLabel(ktrack1,0.4,row0,160),1);
4981 kink->SetTPCDensity(dens00,0,0);
4982 kink->SetTPCDensity(dens01,0,1);
4983 kink->SetTPCDensity(dens10,1,0);
4984 kink->SetTPCDensity(dens11,1,1);
4985 kink->SetIndex(i,0);
4986 kink->SetIndex(j,1);
4989 kink->SetTPCDensity(dens10,0,0);
4990 kink->SetTPCDensity(dens11,0,1);
4991 kink->SetTPCDensity(dens00,1,0);
4992 kink->SetTPCDensity(dens01,1,1);
4993 kink->SetIndex(j,0);
4994 kink->SetIndex(i,1);
4997 if (mpt<1||kink->GetAngle(2)>0.1){
4998 // angle and densities not defined yet
4999 if (kink->GetTPCDensityFactor()<0.8) continue;
5000 if ((2-kink->GetTPCDensityFactor())*kink->GetDistance() >0.25) continue;
5001 if (kink->GetAngle(2)*ktrack0->GetP()<0.003) continue; //too small angle
5002 if (kink->GetAngle(2)>0.2&&kink->GetTPCDensityFactor()<1.15) continue;
5003 if (kink->GetAngle(2)>0.2&&kink->GetTPCDensity(0,1)>0.05) continue;
5005 Float_t criticalangle = track0->GetSigmaSnp2()+track0->GetSigmaTgl2();
5006 criticalangle+= track1->GetSigmaSnp2()+track1->GetSigmaTgl2();
5007 criticalangle= 3*TMath::Sqrt(criticalangle);
5008 if (criticalangle>0.02) criticalangle=0.02;
5009 if (kink->GetAngle(2)<criticalangle) continue;
5012 Int_t drow = Int_t(2.+0.5/(0.05+kink->GetAngle(2))); // overlap region defined
5013 Float_t shapesum =0;
5015 for ( Int_t row = row0-drow; row<row0+drow;row++){
5016 if (row<0) continue;
5017 if (row>155) continue;
5018 if (ktrack0->GetClusterPointer(row)){
5019 AliTPCTrackerPoint *point =ktrack0->GetTrackPoint(row);
5020 shapesum+=point->GetSigmaY()+point->GetSigmaZ();
5023 if (ktrack1->GetClusterPointer(row)){
5024 AliTPCTrackerPoint *point =ktrack1->GetTrackPoint(row);
5025 shapesum+=point->GetSigmaY()+point->GetSigmaZ();
5030 kink->SetShapeFactor(-1.);
5033 kink->SetShapeFactor(shapesum/sum);
5035 // esd->AddKink(kink);
5036 kinks->AddLast(kink);
5042 // sort the kinks according quality - and refit them towards vertex
5044 Int_t nkinks = kinks->GetEntriesFast();
5045 Float_t *quality = new Float_t[nkinks];
5046 Int_t *indexes = new Int_t[nkinks];
5047 AliTPCseed *mothers = new AliTPCseed[nkinks];
5048 AliTPCseed *daughters = new AliTPCseed[nkinks];
5051 for (Int_t i=0;i<nkinks;i++){
5053 AliKink *kinkl = (AliKink*)kinks->At(i);
5055 // refit kinks towards vertex
5057 Int_t index0 = kinkl->GetIndex(0);
5058 Int_t index1 = kinkl->GetIndex(1);
5059 AliTPCseed * ktrack0 = (AliTPCseed*)array->At(index0);
5060 AliTPCseed * ktrack1 = (AliTPCseed*)array->At(index1);
5062 Int_t sumn=ktrack0->GetNumberOfClusters()+ktrack1->GetNumberOfClusters();
5064 // Refit Kink under if too small angle
5066 if (kinkl->GetAngle(2)<0.05){
5067 kinkl->SetTPCRow0(GetRowNumber(kinkl->GetR()));
5068 Int_t row0 = kinkl->GetTPCRow0();
5069 Int_t drow = Int_t(2.+0.5/(0.05+kinkl->GetAngle(2)));
5072 Int_t last = row0-drow;
5073 if (last<40) last=40;
5074 if (last<ktrack0->GetFirstPoint()+25) last = ktrack0->GetFirstPoint()+25;
5075 AliTPCseed* seed0 = ReSeed(ktrack0,last,kFALSE);
5078 Int_t first = row0+drow;
5079 if (first>130) first=130;
5080 if (first>ktrack1->GetLastPoint()-25) first = TMath::Max(ktrack1->GetLastPoint()-25,30);
5081 AliTPCseed* seed1 = ReSeed(ktrack1,first,kTRUE);
5083 if (seed0 && seed1){
5084 kinkl->SetStatus(1,8);
5085 if (RefitKink(*seed0,*seed1,*kinkl)) kinkl->SetStatus(1,9);
5086 row0 = GetRowNumber(kinkl->GetR());
5087 sumn = seed0->GetNumberOfClusters()+seed1->GetNumberOfClusters();
5088 mothers[i] = *seed0;
5089 daughters[i] = *seed1;
5092 delete kinks->RemoveAt(i);
5093 if (seed0) delete seed0;
5094 if (seed1) delete seed1;
5097 if (kinkl->GetDistance()>0.5 || kinkl->GetR()<110 || kinkl->GetR()>240) {
5098 delete kinks->RemoveAt(i);
5099 if (seed0) delete seed0;
5100 if (seed1) delete seed1;
5108 if (kinkl) quality[i] = 160*((0.1+kinkl->GetDistance())*(2.-kinkl->GetTPCDensityFactor()))/(sumn+40.); //the longest -clossest will win
5110 TMath::Sort(nkinks,quality,indexes,kFALSE);
5112 //remove double find kinks
5114 for (Int_t ikink0=1;ikink0<nkinks;ikink0++){
5115 AliKink * kink0 = (AliKink*) kinks->At(indexes[ikink0]);
5116 if (!kink0) continue;
5118 for (Int_t ikink1=0;ikink1<ikink0;ikink1++){
5119 if (!kink0) continue;
5120 AliKink * kink1 = (AliKink*) kinks->At(indexes[ikink1]);
5121 if (!kink1) continue;
5122 // if not close kink continue
5123 if (TMath::Abs(kink1->GetPosition()[2]-kink0->GetPosition()[2])>10) continue;
5124 if (TMath::Abs(kink1->GetPosition()[1]-kink0->GetPosition()[1])>10) continue;
5125 if (TMath::Abs(kink1->GetPosition()[0]-kink0->GetPosition()[0])>10) continue;
5127 AliTPCseed &mother0 = mothers[indexes[ikink0]];
5128 AliTPCseed &daughter0 = daughters[indexes[ikink0]];
5129 AliTPCseed &mother1 = mothers[indexes[ikink1]];
5130 AliTPCseed &daughter1 = daughters[indexes[ikink1]];
5131 Int_t row0 = (kink0->GetTPCRow0()+kink1->GetTPCRow0())/2;
5140 for (Int_t i=0;i<row0;i++){
5141 if (mother0.GetClusterIndex(i)>0 && mother1.GetClusterIndex(i)>0){
5144 if (mother0.GetClusterIndex(i)==mother1.GetClusterIndex(i)){
5151 for (Int_t i=row0;i<158;i++){
5152 if (daughter0.GetClusterIndex(i)>0 && daughter0.GetClusterIndex(i)>0){
5155 if (mother0.GetClusterIndex(i)==mother1.GetClusterIndex(i)){
5161 Float_t ratio = Float_t(same+1)/Float_t(both+1);
5162 Float_t ratiom = Float_t(samem+1)/Float_t(bothm+1);
5163 Float_t ratiod = Float_t(samed+1)/Float_t(bothd+1);
5164 if (ratio>0.3 && ratiom>0.5 &&ratiod>0.5) {
5165 Int_t sum0 = mother0.GetNumberOfClusters()+daughter0.GetNumberOfClusters();
5166 Int_t sum1 = mother1.GetNumberOfClusters()+daughter1.GetNumberOfClusters();
5168 shared[kink0->GetIndex(0)]= kTRUE;
5169 shared[kink0->GetIndex(1)]= kTRUE;
5170 delete kinks->RemoveAt(indexes[ikink0]);
5173 shared[kink1->GetIndex(0)]= kTRUE;
5174 shared[kink1->GetIndex(1)]= kTRUE;
5175 delete kinks->RemoveAt(indexes[ikink1]);
5182 for (Int_t i=0;i<nkinks;i++){
5183 AliKink * kinkl = (AliKink*) kinks->At(indexes[i]);
5184 if (!kinkl) continue;
5185 kinkl->SetTPCRow0(GetRowNumber(kinkl->GetR()));
5186 Int_t index0 = kinkl->GetIndex(0);
5187 Int_t index1 = kinkl->GetIndex(1);
5188 if (circular[index0]||(circular[index1]&&kinkl->GetDistance()>0.2)) continue;
5189 kinkl->SetMultiple(usage[index0],0);
5190 kinkl->SetMultiple(usage[index1],1);
5191 if (kinkl->GetMultiple()[0]+kinkl->GetMultiple()[1]>2) continue;
5192 if (kinkl->GetMultiple()[0]+kinkl->GetMultiple()[1]>0 && quality[indexes[i]]>0.2) continue;
5193 if (kinkl->GetMultiple()[0]+kinkl->GetMultiple()[1]>0 && kinkl->GetDistance()>0.2) continue;
5194 if (circular[index0]||(circular[index1]&&kinkl->GetDistance()>0.1)) continue;
5196 AliTPCseed * ktrack0 = (AliTPCseed*)array->At(index0);
5197 AliTPCseed * ktrack1 = (AliTPCseed*)array->At(index1);
5198 if (!ktrack0 || !ktrack1) continue;
5199 Int_t index = esd->AddKink(kinkl);
5202 if ( ktrack0->GetKinkIndex(0)==0 && ktrack1->GetKinkIndex(0)==0) { //best kink
5203 if (mothers[indexes[i]].GetNumberOfClusters()>20 && daughters[indexes[i]].GetNumberOfClusters()>20 && (mothers[indexes[i]].GetNumberOfClusters()+daughters[indexes[i]].GetNumberOfClusters())>100){
5204 *ktrack0 = mothers[indexes[i]];
5205 *ktrack1 = daughters[indexes[i]];
5209 ktrack0->SetKinkIndex(usage[index0],-(index+1));
5210 ktrack1->SetKinkIndex(usage[index1], (index+1));
5215 // Remove tracks corresponding to shared kink's
5217 for (Int_t i=0;i<nentries;i++){
5218 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5219 if (!track0) continue;
5220 if (track0->GetKinkIndex(0)!=0) continue;
5221 if (shared[i]) delete array->RemoveAt(i);
5226 RemoveUsed2(array,0.5,0.4,30);
5228 for (Int_t i=0;i<nentries;i++){
5229 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5230 if (!track0) continue;
5231 track0->CookdEdx(0.02,0.6);
5235 for (Int_t i=0;i<nentries;i++){
5236 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5237 if (!track0) continue;
5238 if (track0->Pt()<1.4) continue;
5239 //remove double high momenta tracks - overlapped with kink candidates
5242 for (Int_t icl=track0->GetFirstPoint();icl<track0->GetLastPoint(); icl++){
5243 if (track0->GetClusterPointer(icl)!=0){
5245 if (track0->GetClusterPointer(icl)->IsUsed(10)) ishared++;
5248 if (Float_t(ishared+1)/Float_t(all+1)>0.5) {
5249 delete array->RemoveAt(i);
5253 if (track0->GetKinkIndex(0)!=0) continue;
5254 if (track0->GetNumberOfClusters()<80) continue;
5256 AliTPCseed *pmother = new AliTPCseed();
5257 AliTPCseed *pdaughter = new AliTPCseed();
5258 AliKink *pkink = new AliKink;
5260 AliTPCseed & mother = *pmother;
5261 AliTPCseed & daughter = *pdaughter;
5262 AliKink & kinkl = *pkink;
5263 if (CheckKinkPoint(track0,mother,daughter, kinkl)){
5264 if (mother.GetNumberOfClusters()<30||daughter.GetNumberOfClusters()<20) {
5268 continue; //too short tracks
5270 if (mother.Pt()<1.4) {
5276 Int_t row0= kinkl.GetTPCRow0();
5277 if (kinkl.GetDistance()>0.5 || kinkl.GetR()<110. || kinkl.GetR()>240.) {
5284 Int_t index = esd->AddKink(&kinkl);
5285 mother.SetKinkIndex(0,-(index+1));
5286 daughter.SetKinkIndex(0,index+1);
5287 if (mother.GetNumberOfClusters()>50) {
5288 delete array->RemoveAt(i);
5289 array->AddAt(new AliTPCseed(mother),i);
5292 array->AddLast(new AliTPCseed(mother));
5294 array->AddLast(new AliTPCseed(daughter));
5295 for (Int_t icl=0;icl<row0;icl++) {
5296 if (mother.GetClusterPointer(icl)) mother.GetClusterPointer(icl)->Use(20);
5299 for (Int_t icl=row0;icl<158;icl++) {
5300 if (daughter.GetClusterPointer(icl)) daughter.GetClusterPointer(icl)->Use(20);
5309 delete [] daughters;
5331 printf("Ncandidates=\t%d\t%d\t%d\t%d\n",esd->GetNumberOfKinks(),ncandidates,ntracks,nall);
5335 void AliTPCtrackerMI::FindV0s(const TObjArray * array, AliESDEvent *esd)
5341 TObjArray *tpcv0s = new TObjArray(100000);
5342 Int_t nentries = array->GetEntriesFast();
5343 AliHelix *helixes = new AliHelix[nentries];
5344 Int_t *sign = new Int_t[nentries];
5345 Float_t *alpha = new Float_t[nentries];
5346 Float_t *z0 = new Float_t[nentries];
5347 Float_t *dca = new Float_t[nentries];
5348 Float_t *sdcar = new Float_t[nentries];
5349 Float_t *cdcar = new Float_t[nentries];
5350 Float_t *pulldcar = new Float_t[nentries];
5351 Float_t *pulldcaz = new Float_t[nentries];
5352 Float_t *pulldca = new Float_t[nentries];
5353 Bool_t *isPrim = new Bool_t[nentries];
5354 const AliESDVertex * primvertex = esd->GetVertex();
5355 Double_t zvertex = primvertex->GetZv();
5357 // nentries = array->GetEntriesFast();
5359 for (Int_t i=0;i<nentries;i++){
5362 AliTPCseed* track = (AliTPCseed*)array->At(i);
5363 if (!track) continue;
5364 track->GetV0Indexes()[0] = 0; //rest v0 indexes
5365 track->GetV0Indexes()[1] = 0; //rest v0 indexes
5366 track->GetV0Indexes()[2] = 0; //rest v0 indexes
5368 alpha[i] = track->GetAlpha();
5369 new (&helixes[i]) AliHelix(*track);
5371 helixes[i].Evaluate(0,xyz);
5372 sign[i] = (track->GetC()>0) ? -1:1;
5376 if (track->GetProlongation(0,y,z)) z0[i] = z;
5377 dca[i] = track->GetD(0,0);
5379 // dca error parrameterezation + pulls
5381 sdcar[i] = TMath::Sqrt(0.150*0.150+(100*track->GetC())*(100*track->GetC()));
5382 if (TMath::Abs(track->GetTgl())>1) sdcar[i]*=2.5;
5383 cdcar[i] = TMath::Exp((TMath::Abs(track->GetC())-0.0106)*525.3);
5384 pulldcar[i] = (dca[i]-cdcar[i])/sdcar[i];
5385 pulldcaz[i] = (z0[i]-zvertex)/sdcar[i];
5386 pulldca[i] = TMath::Sqrt(pulldcar[i]*pulldcar[i]+pulldcaz[i]*pulldcaz[i]);
5387 if (track->TPCrPID(1)+track->TPCrPID(2)+track->TPCrPID(3)>0.5) {
5388 if (pulldca[i]<3.) isPrim[i]=kTRUE; //pion, muon and Kaon 3 sigma cut
5390 if (track->TPCrPID(4)>0.5) {
5391 if (pulldca[i]<0.5) isPrim[i]=kTRUE; //proton 0.5 sigma cut
5393 if (track->TPCrPID(0)>0.4) {
5394 isPrim[i]=kFALSE; //electron no sigma cut
5401 Int_t ncandidates =0;
5404 Double_t phase[2][2],radius[2];
5410 Float_t fprimvertex[3]={GetX(),GetY(),GetZ()};
5412 Double_t cradius0 = 10*10;
5413 Double_t cradius1 = 200*200;
5416 Double_t cpointAngle = 0.95;
5418 Double_t delta[2]={10000,10000};
5419 for (Int_t i =0;i<nentries;i++){
5420 if (sign[i]==0) continue;
5421 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5422 if (!track0) continue;
5423 if (AliTPCReconstructor::StreamLevel()>1){
5424 TTreeSRedirector &cstream = *fDebugStreamer;
5429 "zvertex="<<zvertex<<
5430 "sdcar0="<<sdcar[i]<<
5431 "cdcar0="<<cdcar[i]<<
5432 "pulldcar0="<<pulldcar[i]<<
5433 "pulldcaz0="<<pulldcaz[i]<<
5434 "pulldca0="<<pulldca[i]<<
5435 "isPrim="<<isPrim[i]<<
5439 if (track0->GetSigned1Pt()<0) continue;
5440 if (track0->GetKinkIndex(0)>0||isPrim[i]) continue; //daughter kink
5442 if (TMath::Abs(helixes[i].GetHelix(4))<0.000000001) continue;
5447 for (Int_t j =0;j<nentries;j++){
5448 AliTPCseed * track1 = (AliTPCseed*)array->At(j);
5449 if (!track1) continue;
5450 if (track1->GetKinkIndex(0)>0 || isPrim[j]) continue; //daughter kink
5451 if (sign[j]*sign[i]>0) continue;
5452 if (TMath::Abs(helixes[j].GetHelix(4))<0.000001) continue;
5453 if (track0->GetCircular()+track1->GetCircular()>1) continue; //circling -returning track
5456 // DCA to prim vertex cut
5462 Int_t npoints = helixes[i].GetRPHIintersections(helixes[j], phase, radius,cdist2);
5463 if (npoints<1) continue;
5467 if (radius[0]<cradius0||radius[0]>cradius1) continue;
5468 helixes[i].LinearDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta[0]);
5469 if (delta[0]>cdist1) continue;
5472 if (TMath::Max(radius[0],radius[1])<cradius0|| TMath::Min(radius[0],radius[1])>cradius1) continue;
5473 helixes[i].LinearDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta[0]);
5474 helixes[i].LinearDCA(helixes[j],phase[1][0],phase[1][1],radius[1],delta[1]);
5475 if (delta[1]<delta[0]) iclosest=1;
5476 if (delta[iclosest]>cdist1) continue;
5478 helixes[i].ParabolicDCA(helixes[j],phase[iclosest][0],phase[iclosest][1],radius[iclosest],delta[iclosest]);
5479 if (radius[iclosest]<cradius0 || radius[iclosest]>cradius1 || delta[iclosest]>cdist1) continue;
5481 Double_t pointAngle = helixes[i].GetPointAngle(helixes[j],phase[iclosest],fprimvertex);
5482 if (pointAngle<cpointAngle) continue;
5484 Bool_t isGamma = kFALSE;
5485 vertex.SetParamP(*track0); //track0 - plus
5486 vertex.SetParamN(*track1); //track1 - minus
5487 vertex.Update(fprimvertex);
5488 if (track0->TPCrPID(0)>0.3&&track1->TPCrPID(0)>0.3&&vertex.GetAnglep()[2]<0.15) isGamma=kTRUE; // gamma conversion candidate
5489 Double_t pointAngle2 = vertex.GetV0CosineOfPointingAngle();
5491 if (vertex.GetV0CosineOfPointingAngle()<cpointAngle && (!isGamma)) continue;// point angle cut
5492 //Bo: if (vertex.GetDist2()>2&&(!isGamma)) continue; // point angle cut
5493 if (vertex.GetDcaV0Daughters()>2&&(!isGamma)) continue;//Bo: // point angle cut
5494 Float_t sigmae = 0.15*0.15;
5495 if (vertex.GetRr()<80)
5496 sigmae += (sdcar[i]*sdcar[i]+sdcar[j]*sdcar[j])*(1.-vertex.GetRr()/80.)*(1.-vertex.GetRr()/80.);
5497 sigmae+= TMath::Sqrt(sigmae);
5498 //Bo: if (vertex.GetDist2()/sigmae>3.&&(!isGamma)) continue;
5499 if (vertex.GetDcaV0Daughters()/sigmae>3.&&(!isGamma)) continue;
5500 Float_t densb0=0,densb1=0,densa0=0,densa1=0;
5501 Int_t row0 = GetRowNumber(vertex.GetRr());
5503 //Bo: if (vertex.GetDist2()>0.2) continue;
5504 if (vertex.GetDcaV0Daughters()>0.2) continue;
5505 densb0 = track0->Density2(0,row0-5);
5506 densb1 = track1->Density2(0,row0-5);
5507 if (densb0>0.3|| densb1>0.3) continue; //clusters before vertex
5508 densa0 = track0->Density2(row0+5,row0+40);
5509 densa1 = track1->Density2(row0+5,row0+40);
5510 if ((densa0<0.4|| densa1<0.4)&&(!isGamma)) continue; //missing clusters after vertex
5513 densa0 = track0->Density2(0,40); //cluster density
5514 densa1 = track1->Density2(0,40); //cluster density
5515 if ((vertex.GetRr()<80&&densa0+densa1<1.)&&(!isGamma)) continue;
5517 //Bo: vertex.SetLab(0,track0->GetLabel());
5518 //Bo: vertex.SetLab(1,track1->GetLabel());
5519 vertex.SetChi2After((densa0+densa1)*0.5);
5520 vertex.SetChi2Before((densb0+densb1)*0.5);
5521 vertex.SetIndex(0,i);
5522 vertex.SetIndex(1,j);
5523 //Bo: vertex.SetStatus(1); // TPC v0 candidate
5524 vertex.SetOnFlyStatus(2);//Bo: // TPC v0 candidate
5525 //Bo: vertex.SetRp(track0->TPCrPIDs());
5526 //Bo: vertex.SetRm(track1->TPCrPIDs());
5527 tpcv0s->AddLast(new AliESDv0(vertex));
5530 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
5531 Double_t radiusm= (delta[0]<delta[1])? TMath::Sqrt(radius[0]):TMath::Sqrt(radius[1]);
5532 Double_t deltam= (delta[0]<delta[1])? TMath::Sqrt(delta[0]):TMath::Sqrt(delta[1]);
5533 if (AliTPCReconstructor::StreamLevel()>1) {
5534 Int_t lab0=track0->GetLabel();
5535 Int_t lab1=track1->GetLabel();
5536 Char_t c0=track0->GetCircular();
5537 Char_t c1=track1->GetCircular();
5538 TTreeSRedirector &cstream = *fDebugStreamer;
5541 "vertex.="<<&vertex<<
5544 "Helix0.="<<&helixes[i]<<
5547 "Helix1.="<<&helixes[j]<<
5548 "pointAngle="<<pointAngle<<
5549 "pointAngle2="<<pointAngle2<<
5554 "zvertex="<<zvertex<<
5557 "npoints="<<npoints<<
5558 "radius0="<<radius[0]<<
5559 "delta0="<<delta[0]<<
5560 "radius1="<<radius[1]<<
5561 "delta1="<<delta[1]<<
5562 "radiusm="<<radiusm<<
5564 "sdcar0="<<sdcar[i]<<
5565 "sdcar1="<<sdcar[j]<<
5566 "cdcar0="<<cdcar[i]<<
5567 "cdcar1="<<cdcar[j]<<
5568 "pulldcar0="<<pulldcar[i]<<
5569 "pulldcar1="<<pulldcar[j]<<
5570 "pulldcaz0="<<pulldcaz[i]<<
5571 "pulldcaz1="<<pulldcaz[j]<<
5572 "pulldca0="<<pulldca[i]<<
5573 "pulldca1="<<pulldca[j]<<
5583 Float_t *quality = new Float_t[ncandidates];
5584 Int_t *indexes = new Int_t[ncandidates];
5586 for (Int_t i=0;i<ncandidates;i++){
5588 AliESDv0 *v0 = (AliESDv0*)tpcv0s->At(i);
5589 quality[i] = 1./(1.00001-v0->GetV0CosineOfPointingAngle()); //base point angle
5590 // quality[i] /= (0.5+v0->GetDist2());
5591 // quality[i] *= v0->GetChi2After(); //density factor
5593 Int_t index0 = v0->GetIndex(0);
5594 Int_t index1 = v0->GetIndex(1);
5595 //Bo: Double_t minpulldca = TMath::Min(2.+pulldca[v0->GetIndex(0)],(2.+pulldca[v0->GetIndex(1)]) ); //pull
5596 Double_t minpulldca = TMath::Min(2.+pulldca[index0],(2.+pulldca[index1]) );//Bo:
5600 AliTPCseed * track0 = (AliTPCseed*)array->At(index0);
5601 AliTPCseed * track1 = (AliTPCseed*)array->At(index1);
5602 if (track0->TPCrPID(0)>0.3&&track1->TPCrPID(0)>0.3&&v0->GetAnglep()[2]<0.15) quality[i]+=1000000; // gamma conversion candidate
5603 if (track0->TPCrPID(4)>0.9||(track1->TPCrPID(4)>0.9&&minpulldca>4)) quality[i]*=10; // lambda candidate candidate
5606 TMath::Sort(ncandidates,quality,indexes,kTRUE);
5609 for (Int_t i=0;i<ncandidates;i++){
5610 AliESDv0 * v0 = (AliESDv0*)tpcv0s->At(indexes[i]);
5612 Int_t index0 = v0->GetIndex(0);
5613 Int_t index1 = v0->GetIndex(1);
5614 AliTPCseed * track0 = (AliTPCseed*)array->At(index0);
5615 AliTPCseed * track1 = (AliTPCseed*)array->At(index1);
5616 if (!track0||!track1) {
5620 Bool_t accept =kTRUE; //default accept
5621 Int_t *v0indexes0 = track0->GetV0Indexes();
5622 Int_t *v0indexes1 = track1->GetV0Indexes();
5624 Int_t order0 = (v0indexes0[0]!=0) ? 1:0;
5625 Int_t order1 = (v0indexes1[0]!=0) ? 1:0;
5626 if (v0indexes0[1]!=0) order0 =2;
5627 if (v0indexes1[1]!=0) order1 =2;
5629 if (v0indexes0[2]!=0) {order0=3; accept=kFALSE;}
5630 if (v0indexes0[2]!=0) {order1=3; accept=kFALSE;}
5632 AliESDv0 * v02 = v0;
5634 //Bo: v0->SetOrder(0,order0);
5635 //Bo: v0->SetOrder(1,order1);
5636 //Bo: v0->SetOrder(1,order0+order1);
5637 v0->SetOnFlyStatus(kTRUE);
5638 Int_t index = esd->AddV0(v0);
5639 v02 = esd->GetV0(index);
5640 v0indexes0[order0]=index;
5641 v0indexes1[order1]=index;
5645 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
5646 if (AliTPCReconstructor::StreamLevel()>1) {
5647 Int_t lab0=track0->GetLabel();
5648 Int_t lab1=track1->GetLabel();
5649 TTreeSRedirector &cstream = *fDebugStreamer;
5658 "dca0="<<dca[index0]<<
5659 "dca1="<<dca[index1]<<
5663 "quality="<<quality[i]<<
5664 "pulldca0="<<pulldca[index0]<<
5665 "pulldca1="<<pulldca[index1]<<
5689 printf("TPC V0 finder : naccepted\t%d\tncandidates\t%d\tntracks\t%d\tnall\t%d\n",naccepted,ncandidates,ntracks,nall);
5693 Int_t AliTPCtrackerMI::RefitKink(AliTPCseed &mother, AliTPCseed &daughter, AliESDkink &knk)
5696 // refit kink towards to the vertex
5699 AliKink &kink=(AliKink &)knk;
5701 Int_t row0 = GetRowNumber(kink.GetR());
5702 FollowProlongation(mother,0);
5703 mother.Reset(kFALSE);
5705 FollowProlongation(daughter,row0);
5706 daughter.Reset(kFALSE);
5707 FollowBackProlongation(daughter,158);
5708 daughter.Reset(kFALSE);
5709 Int_t first = TMath::Max(row0-20,30);
5710 Int_t last = TMath::Min(row0+20,140);
5712 const Int_t kNdiv =5;
5713 AliTPCseed param0[kNdiv]; // parameters along the track
5714 AliTPCseed param1[kNdiv]; // parameters along the track
5715 AliKink kinks[kNdiv]; // corresponding kink parameters
5718 for (Int_t irow=0; irow<kNdiv;irow++){
5719 rows[irow] = first +((last-first)*irow)/(kNdiv-1);
5721 // store parameters along the track
5723 for (Int_t irow=0;irow<kNdiv;irow++){
5724 FollowBackProlongation(mother, rows[irow]);
5725 FollowProlongation(daughter,rows[kNdiv-1-irow]);
5726 param0[irow] = mother;
5727 param1[kNdiv-1-irow] = daughter;
5731 for (Int_t irow=0; irow<kNdiv-1;irow++){
5732 if (param0[irow].GetNumberOfClusters()<kNdiv||param1[irow].GetNumberOfClusters()<kNdiv) continue;
5733 kinks[irow].SetMother(param0[irow]);
5734 kinks[irow].SetDaughter(param1[irow]);
5735 kinks[irow].Update();
5738 // choose kink with best "quality"
5740 Double_t mindist = 10000;
5741 for (Int_t irow=0;irow<kNdiv;irow++){
5742 if (param0[irow].GetNumberOfClusters()<20||param1[irow].GetNumberOfClusters()<20) continue;
5743 if (TMath::Abs(kinks[irow].GetR())>240.) continue;
5744 if (TMath::Abs(kinks[irow].GetR())<100.) continue;
5746 Float_t normdist = TMath::Abs(param0[irow].GetX()-kinks[irow].GetR())*(0.1+kink.GetDistance());
5747 normdist/= (param0[irow].GetNumberOfClusters()+param1[irow].GetNumberOfClusters()+40.);
5748 if (normdist < mindist){
5754 if (index==-1) return 0;
5757 param0[index].Reset(kTRUE);
5758 FollowProlongation(param0[index],0);
5760 mother = param0[index];
5761 daughter = param1[index]; // daughter in vertex
5763 kink.SetMother(mother);
5764 kink.SetDaughter(daughter);
5766 kink.SetTPCRow0(GetRowNumber(kink.GetR()));
5767 kink.SetTPCncls(param0[index].GetNumberOfClusters(),0);
5768 kink.SetTPCncls(param1[index].GetNumberOfClusters(),1);
5769 kink.SetLabel(CookLabel(&mother,0.4, 0,kink.GetTPCRow0()),0);
5770 kink.SetLabel(CookLabel(&daughter,0.4, kink.GetTPCRow0(),160),1);
5771 mother.SetLabel(kink.GetLabel(0));
5772 daughter.SetLabel(kink.GetLabel(1));
5778 void AliTPCtrackerMI::UpdateKinkQualityM(AliTPCseed * seed){
5780 // update Kink quality information for mother after back propagation
5782 if (seed->GetKinkIndex(0)>=0) return;
5783 for (Int_t ikink=0;ikink<3;ikink++){
5784 Int_t index = seed->GetKinkIndex(ikink);
5785 if (index>=0) break;
5786 index = TMath::Abs(index)-1;
5787 AliESDkink * kink = fEvent->GetKink(index);
5788 kink->SetTPCDensity(-1,0,0);
5789 kink->SetTPCDensity(1,0,1);
5791 Int_t row0 = kink->GetTPCRow0() - 2 - Int_t( 0.5/ (0.05+kink->GetAngle(2)));
5792 if (row0<15) row0=15;
5794 Int_t row1 = kink->GetTPCRow0() + 2 + Int_t( 0.5/ (0.05+kink->GetAngle(2)));
5795 if (row1>145) row1=145;
5797 Int_t found,foundable,shared;
5798 seed->GetClusterStatistic(0,row0, found, foundable,shared,kFALSE);
5799 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),0,0);
5800 seed->GetClusterStatistic(row1,155, found, foundable,shared,kFALSE);
5801 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),0,1);
5806 void AliTPCtrackerMI::UpdateKinkQualityD(AliTPCseed * seed){
5808 // update Kink quality information for daughter after refit
5810 if (seed->GetKinkIndex(0)<=0) return;
5811 for (Int_t ikink=0;ikink<3;ikink++){
5812 Int_t index = seed->GetKinkIndex(ikink);
5813 if (index<=0) break;
5814 index = TMath::Abs(index)-1;
5815 AliESDkink * kink = fEvent->GetKink(index);
5816 kink->SetTPCDensity(-1,1,0);
5817 kink->SetTPCDensity(-1,1,1);
5819 Int_t row0 = kink->GetTPCRow0() -2 - Int_t( 0.5/ (0.05+kink->GetAngle(2)));
5820 if (row0<15) row0=15;
5822 Int_t row1 = kink->GetTPCRow0() +2 + Int_t( 0.5/ (0.05+kink->GetAngle(2)));
5823 if (row1>145) row1=145;
5825 Int_t found,foundable,shared;
5826 seed->GetClusterStatistic(0,row0, found, foundable,shared,kFALSE);
5827 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),1,0);
5828 seed->GetClusterStatistic(row1,155, found, foundable,shared,kFALSE);
5829 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),1,1);
5835 Int_t AliTPCtrackerMI::CheckKinkPoint(AliTPCseed*seed,AliTPCseed &mother, AliTPCseed &daughter, AliESDkink &knk)
5838 // check kink point for given track
5839 // if return value=0 kink point not found
5840 // otherwise seed0 correspond to mother particle
5841 // seed1 correspond to daughter particle
5842 // kink parameter of kink point
5843 AliKink &kink=(AliKink &)knk;
5845 Int_t middlerow = (seed->GetFirstPoint()+seed->GetLastPoint())/2;
5846 Int_t first = seed->GetFirstPoint();
5847 Int_t last = seed->GetLastPoint();
5848 if (last-first<20) return 0; // shortest length - 2*30 = 60 pad-rows
5851 AliTPCseed *seed1 = ReSeed(seed,middlerow+20, kTRUE); //middle of chamber
5852 if (!seed1) return 0;
5853 FollowProlongation(*seed1,seed->GetLastPoint()-20);
5854 seed1->Reset(kTRUE);
5855 FollowProlongation(*seed1,158);
5856 seed1->Reset(kTRUE);
5857 last = seed1->GetLastPoint();
5859 AliTPCseed *seed0 = new AliTPCseed(*seed);
5860 seed0->Reset(kFALSE);
5863 AliTPCseed param0[20]; // parameters along the track
5864 AliTPCseed param1[20]; // parameters along the track
5865 AliKink kinks[20]; // corresponding kink parameters
5867 for (Int_t irow=0; irow<20;irow++){
5868 rows[irow] = first +((last-first)*irow)/19;
5870 // store parameters along the track
5872 for (Int_t irow=0;irow<20;irow++){
5873 FollowBackProlongation(*seed0, rows[irow]);
5874 FollowProlongation(*seed1,rows[19-irow]);
5875 param0[irow] = *seed0;
5876 param1[19-irow] = *seed1;
5880 for (Int_t irow=0; irow<19;irow++){
5881 kinks[irow].SetMother(param0[irow]);
5882 kinks[irow].SetDaughter(param1[irow]);
5883 kinks[irow].Update();
5886 // choose kink with biggest change of angle
5888 Double_t maxchange= 0;
5889 for (Int_t irow=1;irow<19;irow++){
5890 if (TMath::Abs(kinks[irow].GetR())>240.) continue;
5891 if (TMath::Abs(kinks[irow].GetR())<110.) continue;
5892 Float_t quality = TMath::Abs(kinks[irow].GetAngle(2))/(3.+TMath::Abs(kinks[irow].GetR()-param0[irow].GetX()));
5893 if ( quality > maxchange){
5894 maxchange = quality;
5901 if (index<0) return 0;
5903 Int_t row0 = GetRowNumber(kinks[index].GetR()); //row 0 estimate
5904 seed0 = new AliTPCseed(param0[index]);
5905 seed1 = new AliTPCseed(param1[index]);
5906 seed0->Reset(kFALSE);
5907 seed1->Reset(kFALSE);
5908 seed0->ResetCovariance(10.);
5909 seed1->ResetCovariance(10.);
5910 FollowProlongation(*seed0,0);
5911 FollowBackProlongation(*seed1,158);
5912 mother = *seed0; // backup mother at position 0
5913 seed0->Reset(kFALSE);
5914 seed1->Reset(kFALSE);
5915 seed0->ResetCovariance(10.);
5916 seed1->ResetCovariance(10.);
5918 first = TMath::Max(row0-20,0);
5919 last = TMath::Min(row0+20,158);
5921 for (Int_t irow=0; irow<20;irow++){
5922 rows[irow] = first +((last-first)*irow)/19;
5924 // store parameters along the track
5926 for (Int_t irow=0;irow<20;irow++){
5927 FollowBackProlongation(*seed0, rows[irow]);
5928 FollowProlongation(*seed1,rows[19-irow]);
5929 param0[irow] = *seed0;
5930 param1[19-irow] = *seed1;
5934 for (Int_t irow=0; irow<19;irow++){
5935 kinks[irow].SetMother(param0[irow]);
5936 kinks[irow].SetDaughter(param1[irow]);
5937 // param0[irow].Dump();
5938 //param1[irow].Dump();
5939 kinks[irow].Update();
5942 // choose kink with biggest change of angle
5945 for (Int_t irow=0;irow<20;irow++){
5946 if (TMath::Abs(kinks[irow].GetR())>250.) continue;
5947 if (TMath::Abs(kinks[irow].GetR())<90.) continue;
5948 Float_t quality = TMath::Abs(kinks[irow].GetAngle(2))/(3.+TMath::Abs(kinks[irow].GetR()-param0[irow].GetX()));
5949 if ( quality > maxchange){
5950 maxchange = quality;
5957 if (index==-1 || param0[index].GetNumberOfClusters()+param1[index].GetNumberOfClusters()<100){
5962 // Float_t anglesigma = TMath::Sqrt(param0[index].fC22+param0[index].fC33+param1[index].fC22+param1[index].fC33);
5964 kink.SetMother(param0[index]);
5965 kink.SetDaughter(param1[index]);
5967 row0 = GetRowNumber(kink.GetR());
5968 kink.SetTPCRow0(row0);
5969 kink.SetLabel(CookLabel(seed0,0.5,0,row0),0);
5970 kink.SetLabel(CookLabel(seed1,0.5,row0,158),1);
5971 kink.SetIndex(-10,0);
5972 kink.SetIndex(int(param0[index].GetNumberOfClusters()+param1[index].GetNumberOfClusters()),1);
5973 kink.SetTPCncls(param0[index].GetNumberOfClusters(),0);
5974 kink.SetTPCncls(param1[index].GetNumberOfClusters(),1);
5977 // new (&mother) AliTPCseed(param0[index]);
5978 daughter = param1[index];
5979 daughter.SetLabel(kink.GetLabel(1));
5980 param0[index].Reset(kTRUE);
5981 FollowProlongation(param0[index],0);
5982 mother = param0[index];
5983 mother.SetLabel(kink.GetLabel(0));
5993 AliTPCseed* AliTPCtrackerMI::ReSeed(AliTPCseed *t)
5996 // reseed - refit - track
5999 // Int_t last = fSectors->GetNRows()-1;
6001 if (fSectors == fOuterSec){
6002 first = TMath::Max(first, t->GetFirstPoint()-fInnerSec->GetNRows());
6006 first = t->GetFirstPoint();
6008 AliTPCseed * seed = MakeSeed(t,0.1,0.5,0.9);
6009 FollowBackProlongation(*t,fSectors->GetNRows()-1);
6011 FollowProlongation(*t,first);
6021 //_____________________________________________________________________________
6022 Int_t AliTPCtrackerMI::ReadSeeds(const TFile *inp) {
6023 //-----------------------------------------------------------------
6024 // This function reades track seeds.
6025 //-----------------------------------------------------------------
6026 TDirectory *savedir=gDirectory;
6028 TFile *in=(TFile*)inp;
6029 if (!in->IsOpen()) {
6030 cerr<<"AliTPCtrackerMI::ReadSeeds(): input file is not open !\n";
6035 TTree *seedTree=(TTree*)in->Get("Seeds");
6037 cerr<<"AliTPCtrackerMI::ReadSeeds(): ";
6038 cerr<<"can't get a tree with track seeds !\n";
6041 AliTPCtrack *seed=new AliTPCtrack;
6042 seedTree->SetBranchAddress("tracks",&seed);
6044 if (fSeeds==0) fSeeds=new TObjArray(15000);
6046 Int_t n=(Int_t)seedTree->GetEntries();
6047 for (Int_t i=0; i<n; i++) {
6048 seedTree->GetEvent(i);
6049 fSeeds->AddLast(new AliTPCseed(*seed/*,seed->GetAlpha()*/));
6058 Int_t AliTPCtrackerMI::Clusters2Tracks (AliESDEvent *esd)
6061 if (fSeeds) DeleteSeeds();
6064 if (!fSeeds) return 1;
6071 //_____________________________________________________________________________
6072 Int_t AliTPCtrackerMI::Clusters2Tracks() {
6073 //-----------------------------------------------------------------
6074 // This is a track finder.
6075 //-----------------------------------------------------------------
6076 TDirectory *savedir=gDirectory;
6080 fSeeds = Tracking();
6083 Info("Clusters2Tracks","Time for tracking: \t");timer.Print();timer.Start();
6085 //activate again some tracks
6086 for (Int_t i=0; i<fSeeds->GetEntriesFast(); i++) {
6087 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
6089 Int_t nc=t.GetNumberOfClusters();
6091 delete fSeeds->RemoveAt(i);
6095 if (pt->GetRemoval()==10) {
6096 if (pt->GetDensityFirst(20)>0.8 || pt->GetDensityFirst(30)>0.8 || pt->GetDensityFirst(40)>0.7)
6097 pt->Desactivate(10); // make track again active
6099 pt->Desactivate(20);
6100 delete fSeeds->RemoveAt(i);
6105 RemoveUsed2(fSeeds,0.85,0.85,0);
6106 if (AliTPCReconstructor::GetRecoParam()->GetDoKinks()) FindKinks(fSeeds,fEvent);
6107 //FindCurling(fSeeds, fEvent,0);
6108 if (AliTPCReconstructor::StreamLevel()>2) FindMultiMC(fSeeds, fEvent,-1); // find multi found tracks
6109 RemoveUsed2(fSeeds,0.5,0.4,20);
6110 FindSplitted(fSeeds, fEvent,0); // find multi found tracks
6111 if (AliTPCReconstructor::StreamLevel()>2) FindMultiMC(fSeeds, fEvent,0); // find multi found tracks
6114 // // refit short tracks
6116 Int_t nseed=fSeeds->GetEntriesFast();
6119 for (Int_t i=0; i<nseed; i++) {
6120 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
6122 Int_t nc=t.GetNumberOfClusters();
6124 delete fSeeds->RemoveAt(i);
6127 CookLabel(pt,0.1); //For comparison only
6128 //if ((pt->IsActive() || (pt->fRemoval==10) )&& nc>50 &&pt->GetNumberOfClusters()>0.4*pt->fNFoundable){
6129 if ((pt->IsActive() || (pt->GetRemoval()==10) )){
6131 if (fDebug>0) cerr<<found<<'\r';
6135 delete fSeeds->RemoveAt(i);
6139 //RemoveOverlap(fSeeds,0.99,7,kTRUE);
6141 //RemoveUsed(fSeeds,0.9,0.9,6);
6143 nseed=fSeeds->GetEntriesFast();
6145 for (Int_t i=0; i<nseed; i++) {
6146 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
6148 Int_t nc=t.GetNumberOfClusters();
6150 delete fSeeds->RemoveAt(i);
6154 t.CookdEdx(0.02,0.6);
6155 // CheckKinkPoint(&t,0.05);
6156 //if ((pt->IsActive() || (pt->fRemoval==10) )&& nc>50 &&pt->GetNumberOfClusters()>0.4*pt->fNFoundable){
6157 if ((pt->IsActive() || (pt->GetRemoval()==10) )){
6165 delete fSeeds->RemoveAt(i);
6166 //AliTPCseed * seed1 = ReSeed(pt,0.05,0.5,1);
6168 // FollowProlongation(*seed1,0);
6169 // Int_t n = seed1->GetNumberOfClusters();
6170 // printf("fP4\t%f\t%f\n",seed1->GetC(),pt->GetC());
6171 // printf("fN\t%d\t%d\n", seed1->GetNumberOfClusters(),pt->GetNumberOfClusters());
6174 //AliTPCseed * seed2 = ReSeed(pt,0.95,0.5,0.05);
6178 SortTracks(fSeeds, 1);
6182 PrepareForBackProlongation(fSeeds,5.);
6183 PropagateBack(fSeeds);
6184 printf("Time for back propagation: \t");timer.Print();timer.Start();
6188 PrepareForProlongation(fSeeds,5.);
6189 PropagateForward2(fSeeds);
6191 printf("Time for FORWARD propagation: \t");timer.Print();timer.Start();
6192 // RemoveUsed(fSeeds,0.7,0.7,6);
6193 //RemoveOverlap(fSeeds,0.9,7,kTRUE);
6195 nseed=fSeeds->GetEntriesFast();
6197 for (Int_t i=0; i<nseed; i++) {
6198 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
6200 Int_t nc=t.GetNumberOfClusters();
6202 delete fSeeds->RemoveAt(i);
6205 t.CookdEdx(0.02,0.6);
6206 // CookLabel(pt,0.1); //For comparison only
6207 //if ((pt->IsActive() || (pt->fRemoval==10) )&& nc>50 &&pt->GetNumberOfClusters()>0.4*pt->fNFoundable){
6208 if ((pt->IsActive() || (pt->fRemoval==10) )){
6209 cerr<<found++<<'\r';
6212 delete fSeeds->RemoveAt(i);
6217 // fNTracks = found;
6219 Info("Clusters2Tracks","Time for overlap removal, track writing and dedx cooking: \t"); timer.Print();timer.Start();
6222 // cerr<<"Number of found tracks : "<<"\t"<<found<<endl;
6223 Info("Clusters2Tracks","Number of found tracks %d",found);
6225 // UnloadClusters();
6230 void AliTPCtrackerMI::Tracking(TObjArray * arr)
6233 // tracking of the seeds
6236 fSectors = fOuterSec;
6237 ParallelTracking(arr,150,63);
6238 fSectors = fOuterSec;
6239 ParallelTracking(arr,63,0);
6242 TObjArray * AliTPCtrackerMI::Tracking(Int_t seedtype, Int_t i1, Int_t i2, Float_t cuts[4], Float_t dy, Int_t dsec)
6247 TObjArray * arr = new TObjArray;
6249 fSectors = fOuterSec;
6252 for (Int_t sec=0;sec<fkNOS;sec++){
6253 if (seedtype==3) MakeSeeds3(arr,sec,i1,i2,cuts,dy, dsec);
6254 if (seedtype==4) MakeSeeds5(arr,sec,i1,i2,cuts,dy);
6255 if (seedtype==2) MakeSeeds2(arr,sec,i1,i2,cuts,dy);
6258 Info("Tracking","\nSeeding - %d\t%d\t%d\t%d\n",seedtype,i1,i2,arr->GetEntriesFast());
6270 TObjArray * AliTPCtrackerMI::Tracking()
6274 if (AliTPCReconstructor::GetRecoParam()->GetSpecialSeeding()) return TrackingSpecial();
6277 Int_t nup=fOuterSec->GetNRows()+fInnerSec->GetNRows();
6279 TObjArray * seeds = new TObjArray;
6288 Float_t fnumber = 3.0;
6289 Float_t fdensity = 3.0;
6294 for (Int_t delta = 0; delta<18; delta+=6){
6298 arr = Tracking(3,nup-1-delta,nup-1-delta-gap,cuts,-1,1);
6299 SumTracks(seeds,arr);
6300 SignClusters(seeds,fnumber,fdensity);
6302 for (Int_t i=2;i<6;i+=2){
6303 // seed high pt tracks
6306 arr = Tracking(3,nup-i-delta,nup-i-delta-gap,cuts,-1,0);
6307 SumTracks(seeds,arr);
6308 SignClusters(seeds,fnumber,fdensity);
6313 // RemoveUsed(seeds,0.9,0.9,1);
6314 // UnsignClusters();
6315 // SignClusters(seeds,fnumber,fdensity);
6319 for (Int_t delta = 20; delta<120; delta+=10){
6321 // seed high pt tracks
6325 arr = Tracking(3,nup-delta,nup-delta-gap,cuts,-1);
6326 SumTracks(seeds,arr);
6327 SignClusters(seeds,fnumber,fdensity);
6332 arr = Tracking(3,nup-delta-5,nup-delta-5-gap,cuts,-1);
6333 SumTracks(seeds,arr);
6334 SignClusters(seeds,fnumber,fdensity);
6345 Info("Tracking()","\n\nPrimary seeding\t%d\n\n",seeds->GetEntriesFast());
6349 // RemoveUsed(seeds,0.75,0.75,1);
6351 //SignClusters(seeds,fnumber,fdensity);
6360 arr = Tracking(4,nup-1,nup-1-gap,cuts,-1);
6361 SumTracks(seeds,arr);
6362 SignClusters(seeds,fnumber,fdensity);
6364 arr = Tracking(4,nup-2,nup-2-gap,cuts,-1);
6365 SumTracks(seeds,arr);
6366 SignClusters(seeds,fnumber,fdensity);
6368 arr = Tracking(4,nup-3,nup-3-gap,cuts,-1);
6369 SumTracks(seeds,arr);
6370 SignClusters(seeds,fnumber,fdensity);
6374 for (Int_t delta = 3; delta<30; delta+=5){
6380 arr = Tracking(4,nup-1-delta,nup-1-delta-gap,cuts,-1);
6381 SumTracks(seeds,arr);
6382 SignClusters(seeds,fnumber,fdensity);
6384 arr = Tracking(4,nup-3-delta,nup-5-delta-gap,cuts,4);
6385 SumTracks(seeds,arr);
6386 SignClusters(seeds,fnumber,fdensity);
6397 Int_t fLastSeedRowSec=AliTPCReconstructor::GetRecoParam()->GetLastSeedRowSec();
6400 for (Int_t delta = 30; delta<fLastSeedRowSec; delta+=10){
6406 arr = Tracking(4,nup-1-delta,nup-1-delta-gap,cuts,-1);
6407 SumTracks(seeds,arr);
6408 SignClusters(seeds,fnumber,fdensity);
6410 arr = Tracking(4,nup-5-delta,nup-5-delta-gap,cuts,5 );
6411 SumTracks(seeds,arr);
6412 SignClusters(seeds,fnumber,fdensity);
6416 Info("Tracking()","\n\nSecondary seeding\t%d\n\n",seeds->GetEntriesFast());
6427 TObjArray * AliTPCtrackerMI::TrackingSpecial()
6430 // seeding adjusted for laser and cosmic tests - short tracks with big inclination angle
6431 // no primary vertex seeding tried
6435 Int_t nup=fOuterSec->GetNRows()+fInnerSec->GetNRows();
6437 TObjArray * seeds = new TObjArray;
6442 Float_t fnumber = 3.0;
6443 Float_t fdensity = 3.0;
6446 cuts[0] = AliTPCReconstructor::GetRecoParam()->GetMaxC(); // max curvature
6447 cuts[1] = 3.5; // max tan(phi) angle for seeding
6448 cuts[2] = 3.; // not used (cut on z primary vertex)
6449 cuts[3] = 3.5; // max tan(theta) angle for seeding
6451 for (Int_t delta = 0; nup-delta-gap-1>0; delta+=3){
6453 arr = Tracking(4,nup-1-delta,nup-1-delta-gap,cuts,-1);
6454 SumTracks(seeds,arr);
6455 SignClusters(seeds,fnumber,fdensity);
6459 Info("Tracking()","\n\nSecondary seeding\t%d\n\n",seeds->GetEntriesFast());
6470 void AliTPCtrackerMI::SumTracks(TObjArray *arr1,TObjArray *&arr2) const
6473 //sum tracks to common container
6474 //remove suspicious tracks
6475 Int_t nseed = arr2->GetEntriesFast();
6476 for (Int_t i=0;i<nseed;i++){
6477 AliTPCseed *pt=(AliTPCseed*)arr2->UncheckedAt(i);
6480 // remove tracks with too big curvature
6482 if (TMath::Abs(pt->GetC())>AliTPCReconstructor::GetRecoParam()->GetMaxC()){
6483 delete arr2->RemoveAt(i);
6486 // REMOVE VERY SHORT TRACKS
6487 if (pt->GetNumberOfClusters()<20){
6488 delete arr2->RemoveAt(i);
6491 // NORMAL ACTIVE TRACK
6492 if (pt->IsActive()){
6493 arr1->AddLast(arr2->RemoveAt(i));
6496 //remove not usable tracks
6497 if (pt->GetRemoval()!=10){
6498 delete arr2->RemoveAt(i);
6502 // ENABLE ONLY ENOUGH GOOD STOPPED TRACKS
6503 if (pt->GetDensityFirst(20)>0.8 || pt->GetDensityFirst(30)>0.8 || pt->GetDensityFirst(40)>0.7)
6504 arr1->AddLast(arr2->RemoveAt(i));
6506 delete arr2->RemoveAt(i);
6510 delete arr2; arr2 = 0;
6515 void AliTPCtrackerMI::ParallelTracking(TObjArray * arr, Int_t rfirst, Int_t rlast)
6518 // try to track in parralel
6520 Int_t nseed=arr->GetEntriesFast();
6521 //prepare seeds for tracking
6522 for (Int_t i=0; i<nseed; i++) {
6523 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i), &t=*pt;
6525 if (!t.IsActive()) continue;
6526 // follow prolongation to the first layer
6527 if ( (fSectors ==fInnerSec) || (t.GetFirstPoint()-fkParam->GetNRowLow()>rfirst+1) )
6528 FollowProlongation(t, rfirst+1);
6533 for (Int_t nr=rfirst; nr>=rlast; nr--){
6534 if (nr<fInnerSec->GetNRows())
6535 fSectors = fInnerSec;
6537 fSectors = fOuterSec;
6538 // make indexes with the cluster tracks for given
6540 // find nearest cluster
6541 for (Int_t i=0; i<nseed; i++) {
6542 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i), &t=*pt;
6544 if (nr==80) pt->UpdateReference();
6545 if (!pt->IsActive()) continue;
6546 // if ( (fSectors ==fOuterSec) && (pt->fFirstPoint-fkParam->GetNRowLow())<nr) continue;
6547 if (pt->GetRelativeSector()>17) {
6550 UpdateClusters(t,nr);
6552 // prolonagate to the nearest cluster - if founded
6553 for (Int_t i=0; i<nseed; i++) {
6554 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
6556 if (!pt->IsActive()) continue;
6557 // if ((fSectors ==fOuterSec) && (pt->fFirstPoint-fkParam->GetNRowLow())<nr) continue;
6558 if (pt->GetRelativeSector()>17) {
6561 FollowToNextCluster(*pt,nr);
6566 void AliTPCtrackerMI::PrepareForBackProlongation(TObjArray * arr,Float_t fac) const
6570 // if we use TPC track itself we have to "update" covariance
6572 Int_t nseed= arr->GetEntriesFast();
6573 for (Int_t i=0;i<nseed;i++){
6574 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
6578 //rotate to current local system at first accepted point
6579 Int_t index = pt->GetClusterIndex2(pt->GetFirstPoint());
6580 Int_t sec = (index&0xff000000)>>24;
6582 Float_t angle1 = fInnerSec->GetAlpha()*sec+fInnerSec->GetAlphaShift();
6583 if (angle1>TMath::Pi())
6584 angle1-=2.*TMath::Pi();
6585 Float_t angle2 = pt->GetAlpha();
6587 if (TMath::Abs(angle1-angle2)>0.001){
6588 pt->Rotate(angle1-angle2);
6589 //angle2 = pt->GetAlpha();
6590 //pt->fRelativeSector = pt->GetAlpha()/fInnerSec->GetAlpha();
6591 //if (pt->GetAlpha()<0)
6592 // pt->fRelativeSector+=18;
6593 //sec = pt->fRelativeSector;
6602 void AliTPCtrackerMI::PrepareForProlongation(TObjArray * arr, Float_t fac) const
6606 // if we use TPC track itself we have to "update" covariance
6608 Int_t nseed= arr->GetEntriesFast();
6609 for (Int_t i=0;i<nseed;i++){
6610 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
6613 pt->SetFirstPoint(pt->GetLastPoint());
6621 Int_t AliTPCtrackerMI::PropagateBack(TObjArray * arr)
6624 // make back propagation
6626 Int_t nseed= arr->GetEntriesFast();
6627 for (Int_t i=0;i<nseed;i++){
6628 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
6629 if (pt&& pt->GetKinkIndex(0)<=0) {
6630 //AliTPCseed *pt2 = new AliTPCseed(*pt);
6631 fSectors = fInnerSec;
6632 //FollowBackProlongation(*pt,fInnerSec->GetNRows()-1);
6633 //fSectors = fOuterSec;
6634 FollowBackProlongation(*pt,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1);
6635 //if (pt->GetNumberOfClusters()<(pt->fEsd->GetTPCclusters(0)) ){
6636 // Error("PropagateBack","Not prolonged track %d",pt->GetLabel());
6637 // FollowBackProlongation(*pt2,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1);
6640 if (pt&& pt->GetKinkIndex(0)>0) {
6641 AliESDkink * kink = fEvent->GetKink(pt->GetKinkIndex(0)-1);
6642 pt->SetFirstPoint(kink->GetTPCRow0());
6643 fSectors = fInnerSec;
6644 FollowBackProlongation(*pt,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1);
6652 Int_t AliTPCtrackerMI::PropagateForward2(TObjArray * arr)
6655 // make forward propagation
6657 Int_t nseed= arr->GetEntriesFast();
6659 for (Int_t i=0;i<nseed;i++){
6660 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
6662 FollowProlongation(*pt,0);
6671 Int_t AliTPCtrackerMI::PropagateForward()
6674 // propagate track forward
6676 Int_t nseed = fSeeds->GetEntriesFast();
6677 for (Int_t i=0;i<nseed;i++){
6678 AliTPCseed *pt = (AliTPCseed*)fSeeds->UncheckedAt(i);
6680 AliTPCseed &t = *pt;
6681 Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
6682 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
6683 if (alpha < 0. ) alpha += 2.*TMath::Pi();
6684 t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN);
6688 fSectors = fOuterSec;
6689 ParallelTracking(fSeeds,fOuterSec->GetNRows()+fInnerSec->GetNRows()-1,fInnerSec->GetNRows());
6690 fSectors = fInnerSec;
6691 ParallelTracking(fSeeds,fInnerSec->GetNRows()-1,0);
6700 Int_t AliTPCtrackerMI::PropagateBack(AliTPCseed * pt, Int_t row0, Int_t row1)
6703 // make back propagation, in between row0 and row1
6707 fSectors = fInnerSec;
6710 if (row1<fSectors->GetNRows())
6713 r1 = fSectors->GetNRows()-1;
6715 if (row0<fSectors->GetNRows()&& r1>0 )
6716 FollowBackProlongation(*pt,r1);
6717 if (row1<=fSectors->GetNRows())
6720 r1 = row1 - fSectors->GetNRows();
6721 if (r1<=0) return 0;
6722 if (r1>=fOuterSec->GetNRows()) return 0;
6723 fSectors = fOuterSec;
6724 return FollowBackProlongation(*pt,r1);
6732 void AliTPCtrackerMI::GetShape(AliTPCseed * seed, Int_t row)
6736 AliTPCClusterParam * clparam = AliTPCcalibDB::Instance()->GetClusterParam();
6737 Float_t zdrift = TMath::Abs((fkParam->GetZLength(0)-TMath::Abs(seed->GetZ())));
6738 Int_t type = (seed->GetSector() < fkParam->GetNSector()/2) ? 0: (row>126) ? 1:2;
6739 Double_t angulary = seed->GetSnp();
6740 angulary = angulary*angulary/(1.-angulary*angulary);
6741 Double_t angularz = seed->GetTgl()*seed->GetTgl()*(1.+angulary);
6743 Double_t sigmay = clparam->GetRMS0(0,type,zdrift,TMath::Sqrt(TMath::Abs(angulary)));
6744 Double_t sigmaz = clparam->GetRMS0(1,type,zdrift,TMath::Sqrt(TMath::Abs(angularz)));
6745 seed->SetCurrentSigmaY2(sigmay*sigmay);
6746 seed->SetCurrentSigmaZ2(sigmaz*sigmaz);
6747 // Float_t sd2 = TMath::Abs((fkParam->GetZLength(0)-TMath::Abs(seed->GetZ())))*fkParam->GetDiffL()*fkParam->GetDiffL();
6748 // // Float_t padlength = fkParam->GetPadPitchLength(seed->fSector);
6749 // Float_t padlength = GetPadPitchLength(row);
6751 // Float_t sresy = (seed->GetSector() < fkParam->GetNSector()/2) ? 0.2 :0.3;
6752 // seed->SetCurrentSigmaY2(sd2+padlength*padlength*angulary/12.+sresy*sresy);
6754 // Float_t sresz = fkParam->GetZSigma();
6755 // seed->SetCurrentSigmaZ2(sd2+padlength*padlength*angularz*angularz*(1+angulary)/12.+sresz*sresz);
6757 Float_t wy = GetSigmaY(seed);
6758 Float_t wz = GetSigmaZ(seed);
6761 if (TMath::Abs(wy/seed->fCurrentSigmaY2-1)>0.0001 || TMath::Abs(wz/seed->fCurrentSigmaZ2-1)>0.0001 ){
6762 printf("problem\n");
6769 //__________________________________________________________________________
6770 void AliTPCtrackerMI::CookLabel(AliKalmanTrack *tk, Float_t wrong) const {
6771 //--------------------------------------------------------------------
6772 //This function "cooks" a track label. If label<0, this track is fake.
6773 //--------------------------------------------------------------------
6774 AliTPCseed * t = dynamic_cast<AliTPCseed*>(tk);
6776 printf("%s:%d wrong type \n",(char*)__FILE__,__LINE__);
6780 Int_t noc=t->GetNumberOfClusters();
6782 //printf("\nnot founded prolongation\n\n\n");
6788 AliTPCclusterMI *clusters[160];
6790 for (Int_t i=0;i<160;i++) {
6797 for (i=0; i<160 && current<noc; i++) {
6799 Int_t index=t->GetClusterIndex2(i);
6800 if (index<=0) continue;
6801 if (index&0x8000) continue;
6803 //clusters[current]=GetClusterMI(index);
6804 if (t->GetClusterPointer(i)){
6805 clusters[current]=t->GetClusterPointer(i);
6811 Int_t lab=123456789;
6812 for (i=0; i<noc; i++) {
6813 AliTPCclusterMI *c=clusters[i];
6815 lab=TMath::Abs(c->GetLabel(0));
6817 for (j=0; j<noc; j++) if (lb[j]==lab || mx[j]==0) break;
6823 for (i=0; i<noc; i++) if (mx[i]>max) {max=mx[i]; lab=lb[i];}
6825 for (i=0; i<noc; i++) {
6826 AliTPCclusterMI *c=clusters[i];
6828 if (TMath::Abs(c->GetLabel(1)) == lab ||
6829 TMath::Abs(c->GetLabel(2)) == lab ) max++;
6832 if ((1.- Float_t(max)/noc) > wrong) lab=-lab;
6835 Int_t tail=Int_t(0.10*noc);
6838 for (i=1; i<=160&&ind<tail; i++) {
6839 // AliTPCclusterMI *c=clusters[noc-i];
6840 AliTPCclusterMI *c=clusters[i];
6842 if (lab == TMath::Abs(c->GetLabel(0)) ||
6843 lab == TMath::Abs(c->GetLabel(1)) ||
6844 lab == TMath::Abs(c->GetLabel(2))) max++;
6847 if (max < Int_t(0.5*tail)) lab=-lab;
6854 //delete[] clusters;
6858 //__________________________________________________________________________
6859 Int_t AliTPCtrackerMI::CookLabel(AliTPCseed *t, Float_t wrong,Int_t first, Int_t last) const {
6860 //--------------------------------------------------------------------
6861 //This function "cooks" a track label. If label<0, this track is fake.
6862 //--------------------------------------------------------------------
6863 Int_t noc=t->GetNumberOfClusters();
6865 //printf("\nnot founded prolongation\n\n\n");
6871 AliTPCclusterMI *clusters[160];
6873 for (Int_t i=0;i<160;i++) {
6880 for (i=0; i<160 && current<noc; i++) {
6881 if (i<first) continue;
6882 if (i>last) continue;
6883 Int_t index=t->GetClusterIndex2(i);
6884 if (index<=0) continue;
6885 if (index&0x8000) continue;
6887 //clusters[current]=GetClusterMI(index);
6888 if (t->GetClusterPointer(i)){
6889 clusters[current]=t->GetClusterPointer(i);
6894 if (noc<5) return -1;
6895 Int_t lab=123456789;
6896 for (i=0; i<noc; i++) {
6897 AliTPCclusterMI *c=clusters[i];
6899 lab=TMath::Abs(c->GetLabel(0));
6901 for (j=0; j<noc; j++) if (lb[j]==lab || mx[j]==0) break;
6907 for (i=0; i<noc; i++) if (mx[i]>max) {max=mx[i]; lab=lb[i];}
6909 for (i=0; i<noc; i++) {
6910 AliTPCclusterMI *c=clusters[i];
6912 if (TMath::Abs(c->GetLabel(1)) == lab ||
6913 TMath::Abs(c->GetLabel(2)) == lab ) max++;
6916 if ((1.- Float_t(max)/noc) > wrong) lab=-lab;
6919 Int_t tail=Int_t(0.10*noc);
6922 for (i=1; i<=160&&ind<tail; i++) {
6923 // AliTPCclusterMI *c=clusters[noc-i];
6924 AliTPCclusterMI *c=clusters[i];
6926 if (lab == TMath::Abs(c->GetLabel(0)) ||
6927 lab == TMath::Abs(c->GetLabel(1)) ||
6928 lab == TMath::Abs(c->GetLabel(2))) max++;
6931 if (max < Int_t(0.5*tail)) lab=-lab;
6934 // t->SetLabel(lab);
6938 //delete[] clusters;
6942 Int_t AliTPCtrackerMI::GetRowNumber(Double_t x[3]) const
6944 //return pad row number for given x vector
6945 Float_t phi = TMath::ATan2(x[1],x[0]);
6946 if(phi<0) phi=2.*TMath::Pi()+phi;
6947 // Get the local angle in the sector philoc
6948 const Float_t kRaddeg = 180/3.14159265358979312;
6949 Float_t phiangle = (Int_t (phi*kRaddeg/20.) + 0.5)*20./kRaddeg;
6950 Double_t localx = x[0]*TMath::Cos(phiangle)-x[1]*TMath::Sin(phiangle);
6951 return GetRowNumber(localx);
6956 void AliTPCtrackerMI::MakeBitmaps(AliTPCseed *t)
6958 //-----------------------------------------------------------------------
6959 // Fill the cluster and sharing bitmaps of the track
6960 //-----------------------------------------------------------------------
6962 Int_t firstpoint = 0;
6963 Int_t lastpoint = 159;
6964 AliTPCTrackerPoint *point;
6965 AliTPCclusterMI *cluster;
6967 for (int iter=firstpoint; iter<lastpoint; iter++) {
6968 // Change to cluster pointers to see if we have a cluster at given padrow
6969 cluster = t->GetClusterPointer(iter);
6971 t->SetClusterMapBit(iter, kTRUE);
6972 point = t->GetTrackPoint(iter);
6973 if (point->IsShared())
6974 t->SetSharedMapBit(iter,kTRUE);
6976 t->SetSharedMapBit(iter, kFALSE);
6979 t->SetClusterMapBit(iter, kFALSE);
6980 t->SetSharedMapBit(iter, kFALSE);
6985 Bool_t AliTPCtrackerMI::IsFindable(AliTPCseed & track){
6987 // return flag if there is findable cluster at given position
6990 Float_t z = track.GetZ();
6992 if (TMath::Abs(z)<(AliTPCReconstructor::GetCtgRange()*track.GetX()+kDeltaZ) &&
6993 TMath::Abs(z)<fkParam->GetZLength(0) &&
6994 (TMath::Abs(track.GetSnp())<AliTPCReconstructor::GetMaxSnpTracker()))
7000 void AliTPCtrackerMI::AddCovariance(AliTPCseed * seed){
7002 // Adding systematic error
7003 // !!!! the systematic error for element 4 is in 1/cm not in pt
7005 const Double_t *param = AliTPCReconstructor::GetRecoParam()->GetSystematicError();
7007 for (Int_t i=0;i<15;i++) covar[i]=0;
7013 covar[0] = param[0]*param[0];
7014 covar[2] = param[1]*param[1];
7015 covar[5] = param[2]*param[2];
7016 covar[9] = param[3]*param[3];
7017 Double_t facC = AliTracker::GetBz()*kB2C;
7018 covar[14]= param[4]*param[4]*facC*facC;
7019 seed->AddCovariance(covar);