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 transform->SetCurrentRecoParam((AliTPCRecoParam*)AliTPCReconstructor::GetRecoParam());
1328 Double_t x[3]={cluster->GetRow(),cluster->GetPad(),cluster->GetTimeBin()};
1329 Int_t i[1]={cluster->GetDetector()};
1330 transform->Transform(x,i,0,1);
1331 // if (cluster->GetDetector()%36>17){
1336 // in debug mode check the transformation
1338 if (AliTPCReconstructor::StreamLevel()>1) {
1340 cluster->GetGlobalXYZ(gx);
1341 Int_t event = (fEvent==NULL)? 0: fEvent->GetEventNumberInFile();
1342 TTreeSRedirector &cstream = *fDebugStreamer;
1343 cstream<<"Transform"<<
1354 cluster->SetX(x[0]);
1355 cluster->SetY(x[1]);
1356 cluster->SetZ(x[2]);
1362 //if (!fkParam->IsGeoRead()) fkParam->ReadGeoMatrices();
1363 TGeoHMatrix *mat = fkParam->GetClusterMatrix(cluster->GetDetector());
1365 Double_t pos[3]= {cluster->GetX(),cluster->GetY(),cluster->GetZ()};
1366 Double_t posC[3]={cluster->GetX(),cluster->GetY(),cluster->GetZ()};
1367 if (mat) mat->LocalToMaster(pos,posC);
1369 // chack Loading of Geo matrices from GeoManager - TEMPORARY FIX
1371 cluster->SetX(posC[0]);
1372 cluster->SetY(posC[1]);
1373 cluster->SetZ(posC[2]);
1376 //_____________________________________________________________________________
1377 Int_t AliTPCtrackerMI::LoadOuterSectors() {
1378 //-----------------------------------------------------------------
1379 // This function fills outer TPC sectors with clusters.
1380 //-----------------------------------------------------------------
1381 Int_t nrows = fOuterSec->GetNRows();
1383 for (Int_t sec = 0;sec<fkNOS;sec++)
1384 for (Int_t row = 0;row<nrows;row++){
1385 AliTPCtrackerRow* tpcrow = &(fOuterSec[sec%fkNOS][row]);
1386 Int_t sec2 = sec+2*fkNIS;
1388 Int_t ncl = tpcrow->GetN1();
1390 AliTPCclusterMI *c= (tpcrow->GetCluster1(ncl));
1391 index=(((sec2<<8)+row)<<16)+ncl;
1392 tpcrow->InsertCluster(c,index);
1395 ncl = tpcrow->GetN2();
1397 AliTPCclusterMI *c= (tpcrow->GetCluster2(ncl));
1398 index=((((sec2+fkNOS)<<8)+row)<<16)+ncl;
1399 tpcrow->InsertCluster(c,index);
1402 // write indexes for fast acces
1404 for (Int_t i=0;i<510;i++)
1405 tpcrow->SetFastCluster(i,-1);
1406 for (Int_t i=0;i<tpcrow->GetN();i++){
1407 Int_t zi = Int_t((*tpcrow)[i]->GetZ()+255.);
1408 tpcrow->SetFastCluster(zi,i); // write index
1411 for (Int_t i=0;i<510;i++){
1412 if (tpcrow->GetFastCluster(i)<0)
1413 tpcrow->SetFastCluster(i,last);
1415 last = tpcrow->GetFastCluster(i);
1424 //_____________________________________________________________________________
1425 Int_t AliTPCtrackerMI::LoadInnerSectors() {
1426 //-----------------------------------------------------------------
1427 // This function fills inner TPC sectors with clusters.
1428 //-----------------------------------------------------------------
1429 Int_t nrows = fInnerSec->GetNRows();
1431 for (Int_t sec = 0;sec<fkNIS;sec++)
1432 for (Int_t row = 0;row<nrows;row++){
1433 AliTPCtrackerRow* tpcrow = &(fInnerSec[sec%fkNIS][row]);
1436 Int_t ncl = tpcrow->GetN1();
1438 AliTPCclusterMI *c= (tpcrow->GetCluster1(ncl));
1439 index=(((sec<<8)+row)<<16)+ncl;
1440 tpcrow->InsertCluster(c,index);
1443 ncl = tpcrow->GetN2();
1445 AliTPCclusterMI *c= (tpcrow->GetCluster2(ncl));
1446 index=((((sec+fkNIS)<<8)+row)<<16)+ncl;
1447 tpcrow->InsertCluster(c,index);
1450 // write indexes for fast acces
1452 for (Int_t i=0;i<510;i++)
1453 tpcrow->SetFastCluster(i,-1);
1454 for (Int_t i=0;i<tpcrow->GetN();i++){
1455 Int_t zi = Int_t((*tpcrow)[i]->GetZ()+255.);
1456 tpcrow->SetFastCluster(zi,i); // write index
1459 for (Int_t i=0;i<510;i++){
1460 if (tpcrow->GetFastCluster(i)<0)
1461 tpcrow->SetFastCluster(i,last);
1463 last = tpcrow->GetFastCluster(i);
1475 //_________________________________________________________________________
1476 AliTPCclusterMI *AliTPCtrackerMI::GetClusterMI(Int_t index) const {
1477 //--------------------------------------------------------------------
1478 // Return pointer to a given cluster
1479 //--------------------------------------------------------------------
1480 if (index<0) return 0; // no cluster
1481 Int_t sec=(index&0xff000000)>>24;
1482 Int_t row=(index&0x00ff0000)>>16;
1483 Int_t ncl=(index&0x00007fff)>>00;
1485 const AliTPCtrackerRow * tpcrow=0;
1486 AliTPCclusterMI * clrow =0;
1488 if (sec<0 || sec>=fkNIS*4) {
1489 AliWarning(Form("Wrong sector %d",sec));
1494 AliTPCtrackerSector& tracksec = fInnerSec[sec%fkNIS];
1495 if (tracksec.GetNRows()<=row) return 0;
1496 tpcrow = &(tracksec[row]);
1497 if (tpcrow==0) return 0;
1500 if (tpcrow->GetN1()<=ncl) return 0;
1501 clrow = tpcrow->GetClusters1();
1504 if (tpcrow->GetN2()<=ncl) return 0;
1505 clrow = tpcrow->GetClusters2();
1509 AliTPCtrackerSector& tracksec = fOuterSec[(sec-fkNIS*2)%fkNOS];
1510 if (tracksec.GetNRows()<=row) return 0;
1511 tpcrow = &(tracksec[row]);
1512 if (tpcrow==0) return 0;
1514 if (sec-2*fkNIS<fkNOS) {
1515 if (tpcrow->GetN1()<=ncl) return 0;
1516 clrow = tpcrow->GetClusters1();
1519 if (tpcrow->GetN2()<=ncl) return 0;
1520 clrow = tpcrow->GetClusters2();
1524 return &(clrow[ncl]);
1530 Int_t AliTPCtrackerMI::FollowToNext(AliTPCseed& t, Int_t nr) {
1531 //-----------------------------------------------------------------
1532 // This function tries to find a track prolongation to next pad row
1533 //-----------------------------------------------------------------
1535 Double_t x= GetXrow(nr), ymax=GetMaxY(nr);
1538 AliTPCclusterMI *cl=0;
1539 Int_t tpcindex= t.GetClusterIndex2(nr);
1541 // update current shape info every 5 pad-row
1542 // if ( (nr%5==0) || t.GetNumberOfClusters()<2 || (t.fCurrentSigmaY2<0.0001) ){
1546 if (fIteration>0 && tpcindex>=-1){ //if we have already clusters
1548 if (tpcindex==-1) return 0; //track in dead zone
1550 cl = t.GetClusterPointer(nr);
1551 if ( (cl==0) ) cl = GetClusterMI(tpcindex);
1552 t.SetCurrentClusterIndex1(tpcindex);
1555 Int_t relativesector = ((tpcindex&0xff000000)>>24)%18; // if previously accepted cluster in different sector
1556 Float_t angle = relativesector*fSectors->GetAlpha()+fSectors->GetAlphaShift();
1558 if (angle<-TMath::Pi()) angle += 2*TMath::Pi();
1559 if (angle>=TMath::Pi()) angle -= 2*TMath::Pi();
1561 if (TMath::Abs(angle-t.GetAlpha())>0.001){
1562 Double_t rotation = angle-t.GetAlpha();
1563 t.SetRelativeSector(relativesector);
1564 if (!t.Rotate(rotation)) return 0;
1566 if (!t.PropagateTo(x)) return 0;
1568 t.SetCurrentCluster(cl);
1570 Int_t accept = AcceptCluster(&t,t.GetCurrentCluster());
1571 if ((tpcindex&0x8000)==0) accept =0;
1573 //if founded cluster is acceptible
1574 if (cl->IsUsed(11)) { // id cluster is shared inrease uncertainty
1575 t.SetErrorY2(t.GetErrorY2()+0.03);
1576 t.SetErrorZ2(t.GetErrorZ2()+0.03);
1577 t.SetErrorY2(t.GetErrorY2()*3);
1578 t.SetErrorZ2(t.GetErrorZ2()*3);
1580 t.SetNFoundable(t.GetNFoundable()+1);
1581 UpdateTrack(&t,accept);
1586 if (TMath::Abs(t.GetSnp())>AliTPCReconstructor::GetMaxSnpTracker()) return 0; // cut on angle
1587 if (fIteration>1 && IsFindable(t)){
1588 // not look for new cluster during refitting
1589 t.SetNFoundable(t.GetNFoundable()+1);
1594 // if (TMath::Abs(t.GetSnp())>0.95 || TMath::Abs(x*t.GetC()-t.GetEta())>0.95) return 0;// patch 28 fev 06
1595 Double_t y=t.GetYat(x);
1596 if (TMath::Abs(y)>ymax){
1598 t.SetRelativeSector((t.GetRelativeSector()+1) % fN);
1599 if (!t.Rotate(fSectors->GetAlpha()))
1601 } else if (y <-ymax) {
1602 t.SetRelativeSector((t.GetRelativeSector()-1+fN) % fN);
1603 if (!t.Rotate(-fSectors->GetAlpha()))
1609 if (!t.PropagateTo(x)) {
1610 if (fIteration==0) t.SetRemoval(10);
1614 Double_t z=t.GetZ();
1617 if (!IsActive(t.GetRelativeSector(),nr)) {
1619 t.SetClusterIndex2(nr,-1);
1622 //AliInfo(Form("A - Sector%d phi %f - alpha %f", t.fRelativeSector,y/x, t.GetAlpha()));
1623 Bool_t isActive = IsActive(t.GetRelativeSector(),nr);
1624 Bool_t isActive2 = (nr>=fInnerSec->GetNRows()) ? fOuterSec[t.GetRelativeSector()][nr-fInnerSec->GetNRows()].GetN()>0:fInnerSec[t.GetRelativeSector()][nr].GetN()>0;
1626 if (!isActive || !isActive2) return 0;
1628 const AliTPCtrackerRow &krow=GetRow(t.GetRelativeSector(),nr);
1629 if ( (t.GetSigmaY2()<0) || t.GetSigmaZ2()<0) return 0;
1631 Double_t roadz = 1.;
1633 if (TMath::Abs(TMath::Abs(y)-ymax)<krow.GetDeadZone()){
1635 t.SetClusterIndex2(nr,-1);
1641 // if (TMath::Abs(z)<(AliTPCReconstructor::GetCtgRange()*x+10) && TMath::Abs(z)<fkParam->GetZLength(0) && (TMath::Abs(t.GetSnp())<AliTPCReconstructor::GetMaxSnpTracker()))
1642 t.SetNFoundable(t.GetNFoundable()+1);
1648 // cl = krow.FindNearest2(y+10.,z,roady,roadz,index);
1649 cl = krow.FindNearest2(y,z,roady,roadz,index);
1650 if (cl) t.SetCurrentClusterIndex1(krow.GetIndex(index));
1653 t.SetCurrentCluster(cl);
1655 if (fIteration==2&&cl->IsUsed(10)) return 0;
1656 Int_t accept = AcceptCluster(&t,t.GetCurrentCluster());
1657 if (fIteration==2&&cl->IsUsed(11)) {
1658 t.SetErrorY2(t.GetErrorY2()+0.03);
1659 t.SetErrorZ2(t.GetErrorZ2()+0.03);
1660 t.SetErrorY2(t.GetErrorY2()*3);
1661 t.SetErrorZ2(t.GetErrorZ2()*3);
1664 if (t.fCurrentCluster->IsUsed(10)){
1669 if (t.fNShared>0.7*t.GetNumberOfClusters()) {
1675 if (accept<3) UpdateTrack(&t,accept);
1678 if ( fIteration==0 && t.GetNFoundable()*0.5 > t.GetNumberOfClusters()) t.SetRemoval(10);
1686 //_________________________________________________________________________
1687 Bool_t AliTPCtrackerMI::GetTrackPoint(Int_t index, AliTrackPoint &p ) const
1689 // Get track space point by index
1690 // return false in case the cluster doesn't exist
1691 AliTPCclusterMI *cl = GetClusterMI(index);
1692 if (!cl) return kFALSE;
1693 Int_t sector = (index&0xff000000)>>24;
1694 // Int_t row = (index&0x00ff0000)>>16;
1696 // xyz[0] = fkParam->GetPadRowRadii(sector,row);
1697 xyz[0] = cl->GetX();
1698 xyz[1] = cl->GetY();
1699 xyz[2] = cl->GetZ();
1701 fkParam->AdjustCosSin(sector,cos,sin);
1702 Float_t x = cos*xyz[0]-sin*xyz[1];
1703 Float_t y = cos*xyz[1]+sin*xyz[0];
1705 Float_t sigmaY2 = 0.027*cl->GetSigmaY2();
1706 if (sector < fkParam->GetNInnerSector()) sigmaY2 *= 2.07;
1707 Float_t sigmaZ2 = 0.066*cl->GetSigmaZ2();
1708 if (sector < fkParam->GetNInnerSector()) sigmaZ2 *= 1.77;
1709 cov[0] = sin*sin*sigmaY2;
1710 cov[1] = -sin*cos*sigmaY2;
1712 cov[3] = cos*cos*sigmaY2;
1715 p.SetXYZ(x,y,xyz[2],cov);
1716 AliGeomManager::ELayerID iLayer;
1718 if (sector < fkParam->GetNInnerSector()) {
1719 iLayer = AliGeomManager::kTPC1;
1723 iLayer = AliGeomManager::kTPC2;
1724 idet = sector - fkParam->GetNInnerSector();
1726 UShort_t volid = AliGeomManager::LayerToVolUID(iLayer,idet);
1727 p.SetVolumeID(volid);
1733 Int_t AliTPCtrackerMI::UpdateClusters(AliTPCseed& t, Int_t nr) {
1734 //-----------------------------------------------------------------
1735 // This function tries to find a track prolongation to next pad row
1736 //-----------------------------------------------------------------
1737 t.SetCurrentCluster(0);
1738 t.SetCurrentClusterIndex1(0);
1740 Double_t xt=t.GetX();
1741 Int_t row = GetRowNumber(xt)-1;
1742 Double_t ymax= GetMaxY(nr);
1744 if (row < nr) return 1; // don't prolongate if not information until now -
1745 // if (TMath::Abs(t.GetSnp())>0.9 && t.GetNumberOfClusters()>40. && fIteration!=2) {
1747 // return 0; // not prolongate strongly inclined tracks
1749 // if (TMath::Abs(t.GetSnp())>0.95) {
1751 // return 0; // not prolongate strongly inclined tracks
1752 // }// patch 28 fev 06
1754 Double_t x= GetXrow(nr);
1756 //t.PropagateTo(x+0.02);
1757 //t.PropagateTo(x+0.01);
1758 if (!t.PropagateTo(x)){
1765 if (TMath::Abs(y)>ymax){
1767 t.SetRelativeSector((t.GetRelativeSector()+1) % fN);
1768 if (!t.Rotate(fSectors->GetAlpha()))
1770 } else if (y <-ymax) {
1771 t.SetRelativeSector((t.GetRelativeSector()-1+fN) % fN);
1772 if (!t.Rotate(-fSectors->GetAlpha()))
1775 // if (!t.PropagateTo(x)){
1782 if (TMath::Abs(t.GetSnp())>AliTPCReconstructor::GetMaxSnpTracker()) return 0;
1784 if (!IsActive(t.GetRelativeSector(),nr)) {
1786 t.SetClusterIndex2(nr,-1);
1789 //AliInfo(Form("A - Sector%d phi %f - alpha %f", t.fRelativeSector,y/x, t.GetAlpha()));
1791 AliTPCtrackerRow &krow=GetRow(t.GetRelativeSector(),nr);
1793 if (TMath::Abs(TMath::Abs(y)-ymax)<krow.GetDeadZone()){
1795 t.SetClusterIndex2(nr,-1);
1801 // if (TMath::Abs(t.GetZ())<(AliTPCReconstructor::GetCtgRange()*t.GetX()+10) && (TMath::Abs(t.GetSnp())<AliTPCReconstructor::GetMaxSnpTracker()))
1802 if (IsFindable(t)) t.SetNFoundable(t.GetNFoundable()+1);
1808 if ( (nr%2==0) || t.GetNumberOfClusters()<2){
1809 // t.fCurrentSigmaY = GetSigmaY(&t);
1810 //t.fCurrentSigmaZ = GetSigmaZ(&t);
1814 AliTPCclusterMI *cl=0;
1817 Double_t roady = 1.;
1818 Double_t roadz = 1.;
1822 index = t.GetClusterIndex2(nr);
1823 if ( (index>0) && (index&0x8000)==0){
1824 cl = t.GetClusterPointer(nr);
1825 if ( (cl==0) && (index>0)) cl = GetClusterMI(index);
1826 t.SetCurrentClusterIndex1(index);
1828 t.SetCurrentCluster(cl);
1834 // if (index<0) return 0;
1835 UInt_t uindex = TMath::Abs(index);
1838 //cl = krow.FindNearest2(y+10,z,roady,roadz,uindex);
1839 cl = krow.FindNearest2(y,z,roady,roadz,uindex);
1842 if (cl) t.SetCurrentClusterIndex1(krow.GetIndex(uindex));
1843 t.SetCurrentCluster(cl);
1849 Int_t AliTPCtrackerMI::FollowToNextCluster(AliTPCseed & t, Int_t nr) {
1850 //-----------------------------------------------------------------
1851 // This function tries to find a track prolongation to next pad row
1852 //-----------------------------------------------------------------
1854 //update error according neighborhoud
1856 if (t.GetCurrentCluster()) {
1858 Int_t accept = AcceptCluster(&t,t.GetCurrentCluster());
1860 if (t.GetCurrentCluster()->IsUsed(10)){
1865 t.SetNShared(t.GetNShared()+1);
1866 if (t.GetNShared()>0.7*t.GetNumberOfClusters()) {
1871 if (fIteration>0) accept = 0;
1872 if (accept<3) UpdateTrack(&t,accept);
1876 if ( ( (t.GetSigmaY2()+t.GetSigmaZ2())>0.16)&& t.GetNumberOfClusters()>18) t.SetRemoval(10);
1877 if ( t.GetChi2()/t.GetNumberOfClusters()>6 &&t.GetNumberOfClusters()>18) t.SetRemoval(10);
1879 if (( (t.GetNFoundable()*0.5 > t.GetNumberOfClusters()) || t.GetNoCluster()>15)) t.SetRemoval(10);
1887 //_____________________________________________________________________________
1888 Int_t AliTPCtrackerMI::FollowProlongation(AliTPCseed& t, Int_t rf, Int_t step) {
1889 //-----------------------------------------------------------------
1890 // This function tries to find a track prolongation.
1891 //-----------------------------------------------------------------
1892 Double_t xt=t.GetX();
1894 Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
1895 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
1896 if (alpha < 0. ) alpha += 2.*TMath::Pi();
1898 t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN);
1900 Int_t first = GetRowNumber(xt)-1;
1901 for (Int_t nr= first; nr>=rf; nr-=step) {
1903 if (t.GetKinkIndexes()[0]>0){
1904 for (Int_t i=0;i<3;i++){
1905 Int_t index = t.GetKinkIndexes()[i];
1906 if (index==0) break;
1907 if (index<0) continue;
1909 AliKink * kink = (AliKink*)fEvent->GetKink(index-1);
1911 printf("PROBLEM\n");
1914 Int_t kinkrow = kink->GetTPCRow0()+2+Int_t(0.5/(0.05+kink->GetAngle(2)));
1916 AliExternalTrackParam paramd(t);
1917 kink->SetDaughter(paramd);
1918 kink->SetStatus(2,5);
1925 if (nr==80) t.UpdateReference();
1926 if (nr<fInnerSec->GetNRows())
1927 fSectors = fInnerSec;
1929 fSectors = fOuterSec;
1930 if (FollowToNext(t,nr)==0)
1943 Int_t AliTPCtrackerMI::FollowBackProlongation(AliTPCseed& t, Int_t rf) {
1944 //-----------------------------------------------------------------
1945 // This function tries to find a track prolongation.
1946 //-----------------------------------------------------------------
1948 Double_t xt=t.GetX();
1949 Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
1950 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
1951 if (alpha < 0. ) alpha += 2.*TMath::Pi();
1952 t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN);
1954 Int_t first = t.GetFirstPoint();
1955 if (first<GetRowNumber(xt)+1) first = GetRowNumber(xt)+1;
1957 if (first<0) first=0;
1958 for (Int_t nr=first; nr<=rf; nr++) {
1959 // if ( (TMath::Abs(t.GetSnp())>0.95)) break;//patch 28 fev 06
1960 if (t.GetKinkIndexes()[0]<0){
1961 for (Int_t i=0;i<3;i++){
1962 Int_t index = t.GetKinkIndexes()[i];
1963 if (index==0) break;
1964 if (index>0) continue;
1965 index = TMath::Abs(index);
1966 AliKink * kink = (AliKink*)fEvent->GetKink(index-1);
1968 printf("PROBLEM\n");
1971 Int_t kinkrow = kink->GetTPCRow0()-2-Int_t(0.5/(0.05+kink->GetAngle(2)));
1973 AliExternalTrackParam paramm(t);
1974 kink->SetMother(paramm);
1975 kink->SetStatus(2,1);
1982 if (nr<fInnerSec->GetNRows())
1983 fSectors = fInnerSec;
1985 fSectors = fOuterSec;
1996 Float_t AliTPCtrackerMI::OverlapFactor(AliTPCseed * s1, AliTPCseed * s2, Int_t &sum1, Int_t & sum2)
2004 Float_t dz2 =(s1->GetZ() - s2->GetZ());
2007 Float_t dy2 =TMath::Abs((s1->GetY() - s2->GetY()));
2009 Float_t distance = TMath::Sqrt(dz2+dy2);
2010 if (distance>4.) return 0; // if there are far away - not overlap - to reduce combinatorics
2013 Int_t firstpoint = TMath::Min(s1->GetFirstPoint(),s2->GetFirstPoint());
2014 Int_t lastpoint = TMath::Max(s1->GetLastPoint(),s2->GetLastPoint());
2019 if (firstpoint>lastpoint) {
2020 firstpoint =lastpoint;
2025 for (Int_t i=firstpoint-1;i<lastpoint+1;i++){
2026 if (s1->GetClusterIndex2(i)>0) sum1++;
2027 if (s2->GetClusterIndex2(i)>0) sum2++;
2028 if (s1->GetClusterIndex2(i)==s2->GetClusterIndex2(i) && s1->GetClusterIndex2(i)>0) {
2032 if (sum<5) return 0;
2034 Float_t summin = TMath::Min(sum1+1,sum2+1);
2035 Float_t ratio = (sum+1)/Float_t(summin);
2039 void AliTPCtrackerMI::SignShared(AliTPCseed * s1, AliTPCseed * s2)
2043 Float_t thetaCut = 0.2;//+10.*TMath::Sqrt(s1->GetSigmaTglZ()+ s2->GetSigmaTglZ());
2044 if (TMath::Abs(s1->GetTgl()-s2->GetTgl())>thetaCut) return;
2045 Float_t minCl = TMath::Min(s1->GetNumberOfClusters(),s2->GetNumberOfClusters());
2046 Int_t cutN0 = TMath::Max(5,TMath::Nint(0.1*minCl));
2051 //Int_t firstpoint = TMath::Max(s1->GetFirstPoint(),s2->GetFirstPoint());
2052 //Int_t lastpoint = TMath::Min(s1->GetLastPoint(),s2->GetLastPoint());
2053 Int_t firstpoint = 0;
2054 Int_t lastpoint = 160;
2056 // if (firstpoint>=lastpoint-5) return;;
2058 for (Int_t i=firstpoint;i<lastpoint;i++){
2059 // if ( (s1->GetClusterIndex2(i)&0xFFFF8FFF)==(s2->GetClusterIndex2(i)&0xFFFF8FFF) && s1->GetClusterIndex2(i)>0) {
2060 if ( (s1->GetClusterIndex2(i))==(s2->GetClusterIndex2(i)) && s1->GetClusterIndex2(i)>0) {
2062 s1->SetSharedMapBit(i, kTRUE);
2063 s2->SetSharedMapBit(i, kTRUE);
2065 if (s1->GetClusterIndex2(i)>0)
2066 s1->SetClusterMapBit(i, kTRUE);
2068 if (sumshared>cutN0){
2071 for (Int_t i=firstpoint;i<lastpoint;i++){
2072 // if ( (s1->GetClusterIndex2(i)&0xFFFF8FFF)==(s2->GetClusterIndex2(i)&0xFFFF8FFF) && s1->GetClusterIndex2(i)>0) {
2073 if ( (s1->GetClusterIndex2(i))==(s2->GetClusterIndex2(i)) && s1->GetClusterIndex2(i)>0) {
2074 AliTPCTrackerPoint *p1 = s1->GetTrackPoint(i);
2075 AliTPCTrackerPoint *p2 = s2->GetTrackPoint(i);;
2076 if (s1->IsActive()&&s2->IsActive()){
2077 p1->SetShared(kTRUE);
2078 p2->SetShared(kTRUE);
2084 if (sumshared>cutN0){
2085 for (Int_t i=0;i<4;i++){
2086 if (s1->GetOverlapLabel(3*i)==0){
2087 s1->SetOverlapLabel(3*i, s2->GetLabel());
2088 s1->SetOverlapLabel(3*i+1,sumshared);
2089 s1->SetOverlapLabel(3*i+2,s2->GetUniqueID());
2093 for (Int_t i=0;i<4;i++){
2094 if (s2->GetOverlapLabel(3*i)==0){
2095 s2->SetOverlapLabel(3*i, s1->GetLabel());
2096 s2->SetOverlapLabel(3*i+1,sumshared);
2097 s2->SetOverlapLabel(3*i+2,s1->GetUniqueID());
2104 void AliTPCtrackerMI::SignShared(TObjArray * arr)
2107 //sort trackss according sectors
2109 for (Int_t i=0; i<arr->GetEntriesFast(); i++) {
2110 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2112 //if (pt) RotateToLocal(pt);
2116 arr->Sort(); // sorting according relative sectors
2117 arr->Expand(arr->GetEntries());
2120 Int_t nseed=arr->GetEntriesFast();
2121 for (Int_t i=0; i<nseed; i++) {
2122 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2124 for (Int_t j=0;j<=12;j++){
2125 pt->SetOverlapLabel(j,0);
2128 for (Int_t i=0; i<nseed; i++) {
2129 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2131 if (pt->GetRemoval()>10) continue;
2132 for (Int_t j=i+1; j<nseed; j++){
2133 AliTPCseed *pt2=(AliTPCseed*)arr->UncheckedAt(j);
2134 if (TMath::Abs(pt->GetRelativeSector()-pt2->GetRelativeSector())>1) continue;
2136 if (pt2->GetRemoval()<=10) {
2137 //if ( TMath::Abs(pt->GetRelativeSector()-pt2->GetRelativeSector())>0) break;
2145 void AliTPCtrackerMI::SortTracks(TObjArray * arr, Int_t mode) const
2148 //sort tracks in array according mode criteria
2149 Int_t nseed = arr->GetEntriesFast();
2150 for (Int_t i=0; i<nseed; i++) {
2151 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2162 void AliTPCtrackerMI::RemoveUsed2(TObjArray * arr, Float_t factor1, Float_t factor2, Int_t minimal)
2165 // Loop over all tracks and remove overlaped tracks (with lower quality)
2167 // 1. Unsign clusters
2168 // 2. Sort tracks according quality
2169 // Quality is defined by the number of cluster between first and last points
2171 // 3. Loop over tracks - decreasing quality order
2172 // a.) remove - If the fraction of shared cluster less than factor (1- n or 2)
2173 // b.) remove - If the minimal number of clusters less than minimal and not ITS
2174 // c.) if track accepted - sign clusters
2176 //Called in - AliTPCtrackerMI::Clusters2Tracks()
2177 // - AliTPCtrackerMI::PropagateBack()
2178 // - AliTPCtrackerMI::RefitInward()
2181 // factor1 - factor for constrained
2182 // factor2 - for non constrained tracks
2183 // if (Float_t(shared+1)/Float_t(found+1)>factor) - DELETE
2187 Int_t nseed = arr->GetEntriesFast();
2188 Float_t * quality = new Float_t[nseed];
2189 Int_t * indexes = new Int_t[nseed];
2193 for (Int_t i=0; i<nseed; i++) {
2194 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2199 pt->UpdatePoints(); //select first last max dens points
2200 Float_t * points = pt->GetPoints();
2201 if (points[3]<0.8) quality[i] =-1;
2202 quality[i] = (points[2]-points[0])+pt->GetNumberOfClusters();
2203 //prefer high momenta tracks if overlaps
2204 quality[i] *= TMath::Sqrt(TMath::Abs(pt->Pt())+0.5);
2206 TMath::Sort(nseed,quality,indexes);
2209 for (Int_t itrack=0; itrack<nseed; itrack++) {
2210 Int_t trackindex = indexes[itrack];
2211 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(trackindex);
2214 if (quality[trackindex]<0){
2216 delete arr->RemoveAt(trackindex);
2219 arr->RemoveAt(trackindex);
2225 Int_t first = Int_t(pt->GetPoints()[0]);
2226 Int_t last = Int_t(pt->GetPoints()[2]);
2227 Double_t factor = (pt->GetBConstrain()) ? factor1: factor2;
2229 Int_t found,foundable,shared;
2230 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
2231 // pt->GetClusterStatistic(0,160, found, foundable,shared,kFALSE);
2232 Bool_t itsgold =kFALSE;
2235 if (pt->GetESD()->GetITSclusters(dummy)>4) itsgold= kTRUE;
2239 if (Float_t(shared+1)/Float_t(found+1)>factor){
2240 if (pt->GetKinkIndexes()[0]!=0) continue; //don't remove tracks - part of the kinks
2241 if( AliTPCReconstructor::StreamLevel()>15){
2242 TTreeSRedirector &cstream = *fDebugStreamer;
2243 cstream<<"RemoveUsed"<<
2244 "iter="<<fIteration<<
2248 delete arr->RemoveAt(trackindex);
2251 if (pt->GetNumberOfClusters()<50&&(found-0.5*shared)<minimal){ //remove short tracks
2252 if (pt->GetKinkIndexes()[0]!=0) continue; //don't remove tracks - part of the kinks
2253 if( AliTPCReconstructor::StreamLevel()>15){
2254 TTreeSRedirector &cstream = *fDebugStreamer;
2255 cstream<<"RemoveShort"<<
2256 "iter="<<fIteration<<
2260 delete arr->RemoveAt(trackindex);
2266 //if (sharedfactor>0.4) continue;
2267 if (pt->GetKinkIndexes()[0]>0) continue;
2268 //Remove tracks with undefined properties - seems
2269 if (pt->GetSigmaY2()<kAlmost0) continue; // ? what is the origin ?
2271 for (Int_t i=first; i<last; i++) {
2272 Int_t index=pt->GetClusterIndex2(i);
2273 // if (index<0 || index&0x8000 ) continue;
2274 if (index<0 || index&0x8000 ) continue;
2275 AliTPCclusterMI *c= pt->GetClusterPointer(i);
2282 Info("RemoveUsed","\n*****\nNumber of good tracks after shared removal\t%d\n",fNtracks);
2288 void AliTPCtrackerMI::UnsignClusters()
2291 // loop over all clusters and unsign them
2294 for (Int_t sec=0;sec<fkNIS;sec++){
2295 for (Int_t row=0;row<fInnerSec->GetNRows();row++){
2296 AliTPCclusterMI *cl = fInnerSec[sec][row].GetClusters1();
2297 for (Int_t icl =0;icl< fInnerSec[sec][row].GetN1();icl++)
2298 // if (cl[icl].IsUsed(10))
2300 cl = fInnerSec[sec][row].GetClusters2();
2301 for (Int_t icl =0;icl< fInnerSec[sec][row].GetN2();icl++)
2302 //if (cl[icl].IsUsed(10))
2307 for (Int_t sec=0;sec<fkNOS;sec++){
2308 for (Int_t row=0;row<fOuterSec->GetNRows();row++){
2309 AliTPCclusterMI *cl = fOuterSec[sec][row].GetClusters1();
2310 for (Int_t icl =0;icl< fOuterSec[sec][row].GetN1();icl++)
2311 //if (cl[icl].IsUsed(10))
2313 cl = fOuterSec[sec][row].GetClusters2();
2314 for (Int_t icl =0;icl< fOuterSec[sec][row].GetN2();icl++)
2315 //if (cl[icl].IsUsed(10))
2324 void AliTPCtrackerMI::SignClusters(const TObjArray * arr, Float_t fnumber, Float_t fdensity)
2327 //sign clusters to be "used"
2329 // snumber and sdensity sign number of sigmas - bellow mean value to be accepted
2330 // loop over "primaries"
2344 Int_t nseed = arr->GetEntriesFast();
2345 for (Int_t i=0; i<nseed; i++) {
2346 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2350 if (!(pt->IsActive())) continue;
2351 Float_t dens = pt->GetNumberOfClusters()/Float_t(pt->GetNFoundable());
2352 if ( (dens>0.7) && (pt->GetNumberOfClusters()>70)){
2354 sumdens2+= dens*dens;
2355 sumn += pt->GetNumberOfClusters();
2356 sumn2 += pt->GetNumberOfClusters()*pt->GetNumberOfClusters();
2357 Float_t chi2 = pt->GetChi2()/pt->GetNumberOfClusters();
2360 sumchi2 +=chi2*chi2;
2365 Float_t mdensity = 0.9;
2366 Float_t meann = 130;
2367 Float_t meanchi = 1;
2368 Float_t sdensity = 0.1;
2369 Float_t smeann = 10;
2370 Float_t smeanchi =0.4;
2374 mdensity = sumdens/sum;
2376 meanchi = sumchi/sum;
2378 sdensity = sumdens2/sum-mdensity*mdensity;
2380 sdensity = TMath::Sqrt(sdensity);
2384 smeann = sumn2/sum-meann*meann;
2386 smeann = TMath::Sqrt(smeann);
2390 smeanchi = sumchi2/sum - meanchi*meanchi;
2392 smeanchi = TMath::Sqrt(smeanchi);
2398 //REMOVE SHORT DELTAS or tracks going out of sensitive volume of TPC
2400 for (Int_t i=0; i<nseed; i++) {
2401 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2405 if (pt->GetBSigned()) continue;
2406 if (pt->GetBConstrain()) continue;
2407 //if (!(pt->IsActive())) continue;
2409 Int_t found,foundable,shared;
2410 pt->GetClusterStatistic(0,160,found, foundable,shared);
2411 if (shared/float(found)>0.3) {
2412 if (shared/float(found)>0.9 ){
2413 //delete arr->RemoveAt(i);
2418 Bool_t isok =kFALSE;
2419 if ( (pt->GetNShared()/pt->GetNumberOfClusters()<0.5) &&pt->GetNumberOfClusters()>60)
2421 if ((TMath::Abs(1/pt->GetC())<100.) && (pt->GetNShared()/pt->GetNumberOfClusters()<0.7))
2423 if (TMath::Abs(pt->GetZ()/pt->GetX())>1.1)
2425 if ( (TMath::Abs(pt->GetSnp()>0.7) && pt->GetD(0,0)>60.))
2429 for (Int_t j=0; j<160; ++j) {
2430 Int_t index=pt->GetClusterIndex2(j);
2431 if (index<0) continue;
2432 AliTPCclusterMI *c= pt->GetClusterPointer(j);
2434 //if (!(c->IsUsed(10))) c->Use();
2441 Double_t maxchi = meanchi+2.*smeanchi;
2443 for (Int_t i=0; i<nseed; i++) {
2444 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2448 //if (!(pt->IsActive())) continue;
2449 if (pt->GetBSigned()) continue;
2450 Double_t chi = pt->GetChi2()/pt->GetNumberOfClusters();
2451 if (chi>maxchi) continue;
2454 Float_t dens = pt->GetNumberOfClusters()/Float_t(pt->GetNFoundable());
2456 //sign only tracks with enoug big density at the beginning
2458 if ((pt->GetDensityFirst(40)<0.75) && pt->GetNumberOfClusters()<meann) continue;
2461 Double_t mindens = TMath::Max(double(mdensity-sdensity*fdensity*bfactor),0.65);
2462 Double_t minn = TMath::Max(Int_t(meann-fnumber*smeann*bfactor),50);
2464 // if (pt->fBConstrain) mindens = TMath::Max(mdensity-sdensity*fdensity*bfactor,0.65);
2465 if ( (pt->GetRemoval()==10) && (pt->GetSnp()>0.8)&&(dens>mindens))
2468 if ((dens>mindens && pt->GetNumberOfClusters()>minn) && chi<maxchi ){
2469 //Int_t noc=pt->GetNumberOfClusters();
2470 pt->SetBSigned(kTRUE);
2471 for (Int_t j=0; j<160; ++j) {
2473 Int_t index=pt->GetClusterIndex2(j);
2474 if (index<0) continue;
2475 AliTPCclusterMI *c= pt->GetClusterPointer(j);
2477 // if (!(c->IsUsed(10))) c->Use();
2482 // gLastCheck = nseed;
2490 void AliTPCtrackerMI::StopNotActive(const TObjArray * arr, Int_t row0, Float_t th0, Float_t th1, Float_t th2) const
2492 // stop not active tracks
2493 // take th1 as threshold for number of founded to number of foundable on last 10 active rows
2494 // take th2 as threshold for number of founded to number of foundable on last 20 active rows
2495 Int_t nseed = arr->GetEntriesFast();
2497 for (Int_t i=0; i<nseed; i++) {
2498 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2502 if (!(pt->IsActive())) continue;
2503 StopNotActive(pt,row0,th0, th1,th2);
2509 void AliTPCtrackerMI::StopNotActive(AliTPCseed * seed, Int_t row0, Float_t th0, Float_t th1,
2512 // stop not active tracks
2513 // take th1 as threshold for number of founded to number of foundable on last 10 active rows
2514 // take th2 as threshold for number of founded to number of foundable on last 20 active rows
2517 Int_t foundable = 0;
2518 Int_t maxindex = seed->GetLastPoint(); //last foundable row
2519 if (seed->GetNFoundable()*th0 > seed->GetNumberOfClusters()) {
2520 seed->Desactivate(10) ;
2524 for (Int_t i=row0; i<maxindex; i++){
2525 Int_t index = seed->GetClusterIndex2(i);
2526 if (index!=-1) foundable++;
2528 if (foundable<=30) sumgood1++;
2529 if (foundable<=50) {
2536 if (foundable>=30.){
2537 if (sumgood1<(th1*30.)) seed->Desactivate(10);
2540 if (sumgood2<(th2*50.)) seed->Desactivate(10);
2544 Int_t AliTPCtrackerMI::RefitInward(AliESDEvent *event)
2547 // back propagation of ESD tracks
2550 const Int_t kMaxFriendTracks=2000;
2554 //PrepareForProlongation(fSeeds,1);
2555 PropagateForward2(fSeeds);
2556 RemoveUsed2(fSeeds,0.4,0.4,20);
2558 TObjArray arraySeed(fSeeds->GetEntries());
2559 for (Int_t i=0;i<fSeeds->GetEntries();i++) {
2560 arraySeed.AddAt(fSeeds->At(i),i);
2562 SignShared(&arraySeed);
2563 // FindCurling(fSeeds, event,2); // find multi found tracks
2564 FindSplitted(fSeeds, event,2); // find multi found tracks
2565 if (AliTPCReconstructor::StreamLevel()>2) FindMultiMC(fSeeds, fEvent,2); // find multi found tracks
2568 Int_t nseed = fSeeds->GetEntriesFast();
2569 for (Int_t i=0;i<nseed;i++){
2570 AliTPCseed * seed = (AliTPCseed*) fSeeds->UncheckedAt(i);
2571 if (!seed) continue;
2572 if (seed->GetKinkIndex(0)>0) UpdateKinkQualityD(seed); // update quality informations for kinks
2573 AliESDtrack *esd=event->GetTrack(i);
2575 if (seed->GetNumberOfClusters()<60 && seed->GetNumberOfClusters()<(esd->GetTPCclusters(0) -5)*0.8){
2576 AliExternalTrackParam paramIn;
2577 AliExternalTrackParam paramOut;
2578 Int_t ncl = seed->RefitTrack(seed,¶mIn,¶mOut);
2579 if (AliTPCReconstructor::StreamLevel()>0) {
2580 (*fDebugStreamer)<<"RecoverIn"<<
2584 "pout.="<<¶mOut<<
2589 seed->Set(paramIn.GetX(),paramIn.GetAlpha(),paramIn.GetParameter(),paramIn.GetCovariance());
2590 seed->SetNumberOfClusters(ncl);
2594 seed->PropagateTo(fkParam->GetInnerRadiusLow());
2595 seed->UpdatePoints();
2596 AddCovariance(seed);
2598 seed->CookdEdx(0.02,0.6);
2599 CookLabel(seed,0.1); //For comparison only
2600 esd->SetTPCClusterMap(seed->GetClusterMap());
2601 esd->SetTPCSharedMap(seed->GetSharedMap());
2603 if (AliTPCReconstructor::StreamLevel()>1 && seed!=0&&esd!=0) {
2604 TTreeSRedirector &cstream = *fDebugStreamer;
2611 if (seed->GetNumberOfClusters()>15){
2612 esd->UpdateTrackParams(seed,AliESDtrack::kTPCrefit);
2613 esd->SetTPCPoints(seed->GetPoints());
2614 esd->SetTPCPointsF(seed->GetNFoundable());
2615 Int_t ndedx = seed->GetNCDEDX(0);
2616 Float_t sdedx = seed->GetSDEDX(0);
2617 Float_t dedx = seed->GetdEdx();
2618 esd->SetTPCsignal(dedx, sdedx, ndedx);
2620 // add seed to the esd track in Calib level
2622 Bool_t storeFriend = gRandom->Rndm()<(kMaxFriendTracks)/Float_t(nseed);
2623 if (AliTPCReconstructor::StreamLevel()>0 &&storeFriend){
2624 AliTPCseed * seedCopy = new AliTPCseed(*seed, kTRUE);
2625 esd->AddCalibObject(seedCopy);
2630 //printf("problem\n");
2633 //FindKinks(fSeeds,event);
2634 Info("RefitInward","Number of refitted tracks %d",ntracks);
2639 Int_t AliTPCtrackerMI::PropagateBack(AliESDEvent *event)
2642 // back propagation of ESD tracks
2648 PropagateBack(fSeeds);
2649 RemoveUsed2(fSeeds,0.4,0.4,20);
2650 //FindCurling(fSeeds, fEvent,1);
2651 FindSplitted(fSeeds, event,1); // find multi found tracks
2652 if (AliTPCReconstructor::StreamLevel()>2) FindMultiMC(fSeeds, fEvent,1); // find multi found tracks
2655 Int_t nseed = fSeeds->GetEntriesFast();
2657 for (Int_t i=0;i<nseed;i++){
2658 AliTPCseed * seed = (AliTPCseed*) fSeeds->UncheckedAt(i);
2659 if (!seed) continue;
2660 if (seed->GetKinkIndex(0)<0) UpdateKinkQualityM(seed); // update quality informations for kinks
2661 seed->UpdatePoints();
2662 AddCovariance(seed);
2663 AliESDtrack *esd=event->GetTrack(i);
2664 if (seed->GetNumberOfClusters()<60 && seed->GetNumberOfClusters()<(esd->GetTPCclusters(0) -5)*0.8){
2665 AliExternalTrackParam paramIn;
2666 AliExternalTrackParam paramOut;
2667 Int_t ncl = seed->RefitTrack(seed,¶mIn,¶mOut);
2668 if (AliTPCReconstructor::StreamLevel()>0) {
2669 (*fDebugStreamer)<<"RecoverBack"<<
2673 "pout.="<<¶mOut<<
2678 seed->Set(paramOut.GetX(),paramOut.GetAlpha(),paramOut.GetParameter(),paramOut.GetCovariance());
2679 seed->SetNumberOfClusters(ncl);
2682 seed->CookdEdx(0.02,0.6);
2683 CookLabel(seed,0.1); //For comparison only
2684 if (seed->GetNumberOfClusters()>15){
2685 esd->UpdateTrackParams(seed,AliESDtrack::kTPCout);
2686 esd->SetTPCPoints(seed->GetPoints());
2687 esd->SetTPCPointsF(seed->GetNFoundable());
2688 Int_t ndedx = seed->GetNCDEDX(0);
2689 Float_t sdedx = seed->GetSDEDX(0);
2690 Float_t dedx = seed->GetdEdx();
2691 esd->SetTPCsignal(dedx, sdedx, ndedx);
2693 Int_t eventnumber = event->GetEventNumberInFile();// patch 28 fev 06
2694 // This is most likely NOT the event number you'd like to use. It has nothing to do with the 'real' event number
2695 if (AliTPCReconstructor::StreamLevel()>1 && esd) {
2696 (*fDebugStreamer)<<"Cback"<<
2699 "EventNrInFile="<<eventnumber<<
2704 //FindKinks(fSeeds,event);
2705 Info("PropagateBack","Number of back propagated tracks %d",ntracks);
2712 void AliTPCtrackerMI::DeleteSeeds()
2717 Int_t nseed = fSeeds->GetEntriesFast();
2718 for (Int_t i=0;i<nseed;i++){
2719 AliTPCseed * seed = (AliTPCseed*)fSeeds->At(i);
2720 if (seed) delete fSeeds->RemoveAt(i);
2727 void AliTPCtrackerMI::ReadSeeds(AliESDEvent *event, Int_t direction)
2730 //read seeds from the event
2732 Int_t nentr=event->GetNumberOfTracks();
2734 Info("ReadSeeds", "Number of ESD tracks: %d\n", nentr);
2739 fSeeds = new TObjArray(nentr);
2743 for (Int_t i=0; i<nentr; i++) {
2744 AliESDtrack *esd=event->GetTrack(i);
2745 ULong_t status=esd->GetStatus();
2746 if (!(status&AliESDtrack::kTPCin)) continue;
2747 AliTPCtrack t(*esd);
2748 t.SetNumberOfClusters(0);
2749 // AliTPCseed *seed = new AliTPCseed(t,t.GetAlpha());
2750 AliTPCseed *seed = new AliTPCseed(t/*,t.GetAlpha()*/);
2751 seed->SetUniqueID(esd->GetID());
2752 AddCovariance(seed); //add systematic ucertainty
2753 for (Int_t ikink=0;ikink<3;ikink++) {
2754 Int_t index = esd->GetKinkIndex(ikink);
2755 seed->GetKinkIndexes()[ikink] = index;
2756 if (index==0) continue;
2757 index = TMath::Abs(index);
2758 AliESDkink * kink = fEvent->GetKink(index-1);
2759 if (kink&&esd->GetKinkIndex(ikink)<0){
2760 if ((status & AliESDtrack::kTRDrefit) != 0) kink->SetStatus(1,2);
2761 if ((status & AliESDtrack::kITSout) != 0) kink->SetStatus(1,0);
2763 if (kink&&esd->GetKinkIndex(ikink)>0){
2764 if ((status & AliESDtrack::kTRDrefit) != 0) kink->SetStatus(1,6);
2765 if ((status & AliESDtrack::kITSout) != 0) kink->SetStatus(1,4);
2769 if (((status&AliESDtrack::kITSout)==0)&&(direction==1)) seed->ResetCovariance(10.);
2770 if ( direction ==2 &&(status & AliESDtrack::kTRDrefit) == 0 ) seed->ResetCovariance(10.);
2771 //if ( direction ==2 && ((status & AliESDtrack::kTPCout) == 0) ) {
2772 // fSeeds->AddAt(0,i);
2776 if ( direction ==2 &&(status & AliESDtrack::kTRDrefit) > 0 ) {
2777 Double_t par0[5],par1[5],alpha,x;
2778 esd->GetInnerExternalParameters(alpha,x,par0);
2779 esd->GetExternalParameters(x,par1);
2780 Double_t delta1 = TMath::Abs(par0[4]-par1[4])/(0.000000001+TMath::Abs(par0[4]+par1[4]));
2781 Double_t delta2 = TMath::Abs(par0[3]-par1[3]);
2783 if (esd->GetTRDncls()>0) trdchi2 = esd->GetTRDchi2()/esd->GetTRDncls();
2784 //reset covariance if suspicious
2785 if ( (delta1>0.1) || (delta2>0.006) ||trdchi2>7.)
2786 seed->ResetCovariance(10.);
2791 // rotate to the local coordinate system
2793 fSectors=fInnerSec; fN=fkNIS;
2794 Double_t alpha=seed->GetAlpha() - fSectors->GetAlphaShift();
2795 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
2796 if (alpha < 0. ) alpha += 2.*TMath::Pi();
2797 Int_t ns=Int_t(alpha/fSectors->GetAlpha())%fN;
2798 alpha =ns*fSectors->GetAlpha() + fSectors->GetAlphaShift();
2799 if (alpha<-TMath::Pi()) alpha += 2*TMath::Pi();
2800 if (alpha>=TMath::Pi()) alpha -= 2*TMath::Pi();
2801 alpha-=seed->GetAlpha();
2802 if (!seed->Rotate(alpha)) {
2808 if (esd->GetKinkIndex(0)<=0){
2809 for (Int_t irow=0;irow<160;irow++){
2810 Int_t index = seed->GetClusterIndex2(irow);
2813 AliTPCclusterMI * cl = GetClusterMI(index);
2814 seed->SetClusterPointer(irow,cl);
2816 if ((index & 0x8000)==0){
2817 cl->Use(10); // accepted cluster
2819 cl->Use(6); // close cluster not accepted
2822 Info("ReadSeeds","Not found cluster");
2827 fSeeds->AddAt(seed,i);
2833 //_____________________________________________________________________________
2834 void AliTPCtrackerMI::MakeSeeds3(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, Float_t cuts[4],
2835 Float_t deltay, Int_t ddsec) {
2836 //-----------------------------------------------------------------
2837 // This function creates track seeds.
2838 // SEEDING WITH VERTEX CONSTRAIN
2839 //-----------------------------------------------------------------
2840 // cuts[0] - fP4 cut
2841 // cuts[1] - tan(phi) cut
2842 // cuts[2] - zvertex cut
2843 // cuts[3] - fP3 cut
2851 Double_t x[5], c[15];
2852 // Int_t di = i1-i2;
2854 AliTPCseed * seed = new AliTPCseed();
2855 Double_t alpha=fSectors->GetAlpha(), shift=fSectors->GetAlphaShift();
2856 Double_t cs=cos(alpha), sn=sin(alpha);
2858 // Double_t x1 =fOuterSec->GetX(i1);
2859 //Double_t xx2=fOuterSec->GetX(i2);
2861 Double_t x1 =GetXrow(i1);
2862 Double_t xx2=GetXrow(i2);
2864 Double_t x3=GetX(), y3=GetY(), z3=GetZ();
2866 Int_t imiddle = (i2+i1)/2; //middle pad row index
2867 Double_t xm = GetXrow(imiddle); // radius of middle pad-row
2868 const AliTPCtrackerRow& krm=GetRow(sec,imiddle); //middle pad -row
2872 const AliTPCtrackerRow& kr1=GetRow(ns,i1);
2873 Double_t ymax = GetMaxY(i1)-kr1.GetDeadZone()-1.5;
2874 Double_t ymaxm = GetMaxY(imiddle)-kr1.GetDeadZone()-1.5;
2877 // change cut on curvature if it can't reach this layer
2878 // maximal curvature set to reach it
2879 Double_t dvertexmax = TMath::Sqrt((x1-x3)*(x1-x3)+(ymax+5-y3)*(ymax+5-y3));
2880 if (dvertexmax*0.5*cuts[0]>0.85){
2881 cuts[0] = 0.85/(dvertexmax*0.5+1.);
2883 Double_t r2min = 1/(cuts[0]*cuts[0]); //minimal square of radius given by cut
2886 if (deltay>0) ddsec = 0;
2887 // loop over clusters
2888 for (Int_t is=0; is < kr1; is++) {
2890 if (kr1[is]->IsUsed(10)) continue;
2891 Double_t y1=kr1[is]->GetY(), z1=kr1[is]->GetZ();
2892 //if (TMath::Abs(y1)>ymax) continue;
2894 if (deltay>0 && TMath::Abs(ymax-TMath::Abs(y1))> deltay ) continue; // seed only at the edge
2896 // find possible directions
2897 Float_t anglez = (z1-z3)/(x1-x3);
2898 Float_t extraz = z1 - anglez*(x1-xx2); // extrapolated z
2901 //find rotation angles relative to line given by vertex and point 1
2902 Double_t dvertex2 = (x1-x3)*(x1-x3)+(y1-y3)*(y1-y3);
2903 Double_t dvertex = TMath::Sqrt(dvertex2);
2904 Double_t angle13 = TMath::ATan((y1-y3)/(x1-x3));
2905 Double_t cs13 = cos(-angle13), sn13 = sin(-angle13);
2908 // loop over 2 sectors
2914 Double_t dddz1=0; // direction of delta inclination in z axis
2921 for (Int_t dsec = dsec1; dsec<=dsec2;dsec++){
2922 Int_t sec2 = sec + dsec;
2924 // AliTPCtrackerRow& kr2 = fOuterSec[(sec2+fkNOS)%fkNOS][i2];
2925 //AliTPCtrackerRow& kr2m = fOuterSec[(sec2+fkNOS)%fkNOS][imiddle];
2926 AliTPCtrackerRow& kr2 = GetRow((sec2+fkNOS)%fkNOS,i2);
2927 AliTPCtrackerRow& kr2m = GetRow((sec2+fkNOS)%fkNOS,imiddle);
2928 Int_t index1 = TMath::Max(kr2.Find(extraz-0.6-dddz1*TMath::Abs(z1)*0.05)-1,0);
2929 Int_t index2 = TMath::Min(kr2.Find(extraz+0.6+dddz2*TMath::Abs(z1)*0.05)+1,kr2);
2931 // rotation angles to p1-p3
2932 Double_t cs13r = cos(-angle13+dsec*alpha)/dvertex, sn13r = sin(-angle13+dsec*alpha)/dvertex;
2933 Double_t x2, y2, z2;
2935 // Double_t dymax = maxangle*TMath::Abs(x1-xx2);
2938 Double_t dxx0 = (xx2-x3)*cs13r;
2939 Double_t dyy0 = (xx2-x3)*sn13r;
2940 for (Int_t js=index1; js < index2; js++) {
2941 const AliTPCclusterMI *kcl = kr2[js];
2942 if (kcl->IsUsed(10)) continue;
2944 //calcutate parameters
2946 Double_t yy0 = dyy0 +(kcl->GetY()-y3)*cs13r;
2948 if (TMath::Abs(yy0)<0.000001) continue;
2949 Double_t xx0 = dxx0 -(kcl->GetY()-y3)*sn13r;
2950 Double_t y0 = 0.5*(xx0*xx0+yy0*yy0-xx0)/yy0;
2951 Double_t r02 = (0.25+y0*y0)*dvertex2;
2952 //curvature (radius) cut
2953 if (r02<r2min) continue;
2957 Double_t c0 = 1/TMath::Sqrt(r02);
2961 //Double_t dfi0 = 2.*TMath::ASin(dvertex*c0*0.5);
2962 //Double_t dfi1 = 2.*TMath::ASin(TMath::Sqrt(yy0*yy0+(1-xx0)*(1-xx0))*dvertex*c0*0.5);
2963 Double_t dfi0 = 2.*AliTPCFastMath::FastAsin(dvertex*c0*0.5);
2964 Double_t dfi1 = 2.*AliTPCFastMath::FastAsin(TMath::Sqrt(yy0*yy0+(1-xx0)*(1-xx0))*dvertex*c0*0.5);
2967 Double_t z0 = kcl->GetZ();
2968 Double_t zzzz2 = z1-(z1-z3)*dfi1/dfi0;
2969 if (TMath::Abs(zzzz2-z0)>0.5) continue;
2972 Double_t dip = (z1-z0)*c0/dfi1;
2973 Double_t x0 = (0.5*cs13+y0*sn13)*dvertex*c0;
2984 x2= xx2*cs-y2*sn*dsec;
2985 y2=+xx2*sn*dsec+y2*cs;
2995 // do we have cluster at the middle ?
2997 GetProlongation(x1,xm,x,ym,zm);
2999 AliTPCclusterMI * cm=0;
3000 if (TMath::Abs(ym)-ymaxm<0){
3001 cm = krm.FindNearest2(ym,zm,1.0,0.6,dummy);
3002 if ((!cm) || (cm->IsUsed(10))) {
3007 // rotate y1 to system 0
3008 // get state vector in rotated system
3009 Double_t yr1 = (-0.5*sn13+y0*cs13)*dvertex*c0;
3010 Double_t xr2 = x0*cs+yr1*sn*dsec;
3011 Double_t xr[5]={kcl->GetY(),kcl->GetZ(), xr2, dip, c0};
3013 GetProlongation(xx2,xm,xr,ym,zm);
3014 if (TMath::Abs(ym)-ymaxm<0){
3015 cm = kr2m.FindNearest2(ym,zm,1.0,0.6,dummy);
3016 if ((!cm) || (cm->IsUsed(10))) {
3026 dym = ym - cm->GetY();
3027 dzm = zm - cm->GetZ();
3034 Double_t sy1=kr1[is]->GetSigmaY2()*2., sz1=kr1[is]->GetSigmaZ2()*2.;
3035 Double_t sy2=kcl->GetSigmaY2()*2., sz2=kcl->GetSigmaZ2()*2.;
3036 //Double_t sy3=400*3./12., sy=0.1, sz=0.1;
3037 Double_t sy3=25000*x[4]*x[4]+0.1, sy=0.1, sz=0.1;
3038 //Double_t sy3=25000*x[4]*x[4]*60+0.5, sy=0.1, sz=0.1;
3040 Double_t f40=(F1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
3041 Double_t f42=(F1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
3042 Double_t f43=(F1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
3043 Double_t f20=(F2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
3044 Double_t f22=(F2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
3045 Double_t f23=(F2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
3047 Double_t f30=(F3(x1,y1+sy,x2,y2,z1,z2)-x[3])/sy;
3048 Double_t f31=(F3(x1,y1,x2,y2,z1+sz,z2)-x[3])/sz;
3049 Double_t f32=(F3(x1,y1,x2,y2+sy,z1,z2)-x[3])/sy;
3050 Double_t f34=(F3(x1,y1,x2,y2,z1,z2+sz)-x[3])/sz;
3054 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3055 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3056 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3057 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3058 c[13]=f30*sy1*f40+f32*sy2*f42;
3059 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3061 // if (!BuildSeed(kr1[is],kcl,0,x1,x2,x3,x,c)) continue;
3063 UInt_t index=kr1.GetIndex(is);
3064 seed->~AliTPCseed(); // this does not set the pointer to 0...
3065 AliTPCseed *track=new(seed) AliTPCseed(x1, ns*alpha+shift, x, c, index);
3067 track->SetIsSeeding(kTRUE);
3068 track->SetSeed1(i1);
3069 track->SetSeed2(i2);
3070 track->SetSeedType(3);
3074 FollowProlongation(*track, (i1+i2)/2,1);
3075 Int_t foundable,found,shared;
3076 track->GetClusterStatistic((i1+i2)/2,i1, found, foundable, shared, kTRUE);
3077 if ((found<0.55*foundable) || shared>0.5*found || (track->GetSigmaY2()+track->GetSigmaZ2())>0.5){
3079 seed->~AliTPCseed();
3085 FollowProlongation(*track, i2,1);
3089 track->SetBConstrain(1);
3090 // track->fLastPoint = i1+fInnerSec->GetNRows(); // first cluster in track position
3091 track->SetLastPoint(i1); // first cluster in track position
3092 track->SetFirstPoint(track->GetLastPoint());
3094 if (track->GetNumberOfClusters()<(i1-i2)*0.5 ||
3095 track->GetNumberOfClusters() < track->GetNFoundable()*0.6 ||
3096 track->GetNShared()>0.4*track->GetNumberOfClusters() ) {
3098 seed->~AliTPCseed();
3102 // Z VERTEX CONDITION
3103 Double_t zv, bz=GetBz();
3104 if ( !track->GetZAt(0.,bz,zv) ) continue;
3105 if (TMath::Abs(zv-z3)>cuts[2]) {
3106 FollowProlongation(*track, TMath::Max(i2-20,0));
3107 if ( !track->GetZAt(0.,bz,zv) ) continue;
3108 if (TMath::Abs(zv-z3)>cuts[2]){
3109 FollowProlongation(*track, TMath::Max(i2-40,0));
3110 if ( !track->GetZAt(0.,bz,zv) ) continue;
3111 if (TMath::Abs(zv-z3)>cuts[2] &&(track->GetNumberOfClusters() > track->GetNFoundable()*0.7)){
3112 // make seed without constrain
3113 AliTPCseed * track2 = MakeSeed(track,0.2,0.5,1.);
3114 FollowProlongation(*track2, i2,1);
3115 track2->SetBConstrain(kFALSE);
3116 track2->SetSeedType(1);
3117 arr->AddLast(track2);
3119 seed->~AliTPCseed();
3124 seed->~AliTPCseed();
3131 track->SetSeedType(0);
3132 arr->AddLast(track);
3133 seed = new AliTPCseed;
3135 // don't consider other combinations
3136 if (track->GetNumberOfClusters() > track->GetNFoundable()*0.8)
3142 Info("MakeSeeds3","\nSeeding statistic:\t%d\t%d\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin,nout1,nout2);
3148 void AliTPCtrackerMI::MakeSeeds5(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, Float_t cuts[4],
3153 //-----------------------------------------------------------------
3154 // This function creates track seeds.
3155 //-----------------------------------------------------------------
3156 // cuts[0] - fP4 cut
3157 // cuts[1] - tan(phi) cut
3158 // cuts[2] - zvertex cut
3159 // cuts[3] - fP3 cut
3169 Double_t x[5], c[15];
3171 // make temporary seed
3172 AliTPCseed * seed = new AliTPCseed;
3173 Double_t alpha=fOuterSec->GetAlpha(), shift=fOuterSec->GetAlphaShift();
3174 // Double_t cs=cos(alpha), sn=sin(alpha);
3179 Double_t x1 = GetXrow(i1-1);
3180 const AliTPCtrackerRow& kr1=GetRow(sec,i1-1);
3181 Double_t y1max = GetMaxY(i1-1)-kr1.GetDeadZone()-1.5;
3183 Double_t x1p = GetXrow(i1);
3184 const AliTPCtrackerRow& kr1p=GetRow(sec,i1);
3186 Double_t x1m = GetXrow(i1-2);
3187 const AliTPCtrackerRow& kr1m=GetRow(sec,i1-2);
3190 //last 3 padrow for seeding
3191 AliTPCtrackerRow& kr3 = GetRow((sec+fkNOS)%fkNOS,i1-7);
3192 Double_t x3 = GetXrow(i1-7);
3193 // Double_t y3max= GetMaxY(i1-7)-kr3.fDeadZone-1.5;
3195 AliTPCtrackerRow& kr3p = GetRow((sec+fkNOS)%fkNOS,i1-6);
3196 Double_t x3p = GetXrow(i1-6);
3198 AliTPCtrackerRow& kr3m = GetRow((sec+fkNOS)%fkNOS,i1-8);
3199 Double_t x3m = GetXrow(i1-8);
3204 Int_t im = i1-4; //middle pad row index
3205 Double_t xm = GetXrow(im); // radius of middle pad-row
3206 const AliTPCtrackerRow& krm=GetRow(sec,im); //middle pad -row
3207 // Double_t ymmax = GetMaxY(im)-kr1.fDeadZone-1.5;
3210 Double_t deltax = x1-x3;
3211 Double_t dymax = deltax*cuts[1];
3212 Double_t dzmax = deltax*cuts[3];
3214 // loop over clusters
3215 for (Int_t is=0; is < kr1; is++) {
3217 if (kr1[is]->IsUsed(10)) continue;
3218 Double_t y1=kr1[is]->GetY(), z1=kr1[is]->GetZ();
3220 if (deltay>0 && TMath::Abs(y1max-TMath::Abs(y1))> deltay ) continue; // seed only at the edge
3222 Int_t index1 = TMath::Max(kr3.Find(z1-dzmax)-1,0);
3223 Int_t index2 = TMath::Min(kr3.Find(z1+dzmax)+1,kr3);
3229 for (Int_t js=index1; js < index2; js++) {
3230 const AliTPCclusterMI *kcl = kr3[js];
3231 if (kcl->IsUsed(10)) continue;
3233 // apply angular cuts
3234 if (TMath::Abs(y1-y3)>dymax) continue;
3237 if (TMath::Abs(z1-z3)>dzmax) continue;
3239 Double_t angley = (y1-y3)/(x1-x3);
3240 Double_t anglez = (z1-z3)/(x1-x3);
3242 Double_t erry = TMath::Abs(angley)*(x1-x1m)*0.5+0.5;
3243 Double_t errz = TMath::Abs(anglez)*(x1-x1m)*0.5+0.5;
3245 Double_t yyym = angley*(xm-x1)+y1;
3246 Double_t zzzm = anglez*(xm-x1)+z1;
3248 const AliTPCclusterMI *kcm = krm.FindNearest2(yyym,zzzm,erry,errz,index);
3250 if (kcm->IsUsed(10)) continue;
3252 erry = TMath::Abs(angley)*(x1-x1m)*0.4+0.5;
3253 errz = TMath::Abs(anglez)*(x1-x1m)*0.4+0.5;
3260 // look around first
3261 const AliTPCclusterMI *kc1m = kr1m.FindNearest2(angley*(x1m-x1)+y1,
3267 if (kc1m->IsUsed(10)) used++;
3269 const AliTPCclusterMI *kc1p = kr1p.FindNearest2(angley*(x1p-x1)+y1,
3275 if (kc1p->IsUsed(10)) used++;
3277 if (used>1) continue;
3278 if (found<1) continue;
3282 const AliTPCclusterMI *kc3m = kr3m.FindNearest2(angley*(x3m-x3)+y3,
3288 if (kc3m->IsUsed(10)) used++;
3292 const AliTPCclusterMI *kc3p = kr3p.FindNearest2(angley*(x3p-x3)+y3,
3298 if (kc3p->IsUsed(10)) used++;
3302 if (used>1) continue;
3303 if (found<3) continue;
3313 x[4]=F1(x1,y1,x2,y2,x3,y3);
3314 //if (TMath::Abs(x[4]) >= cuts[0]) continue;
3317 x[2]=F2(x1,y1,x2,y2,x3,y3);
3320 x[3]=F3n(x1,y1,x2,y2,z1,z2,x[4]);
3321 //if (TMath::Abs(x[3]) > cuts[3]) continue;
3325 Double_t sy1=0.1, sz1=0.1;
3326 Double_t sy2=0.1, sz2=0.1;
3327 Double_t sy3=0.1, sy=0.1, sz=0.1;
3329 Double_t f40=(F1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
3330 Double_t f42=(F1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
3331 Double_t f43=(F1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
3332 Double_t f20=(F2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
3333 Double_t f22=(F2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
3334 Double_t f23=(F2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
3336 Double_t f30=(F3(x1,y1+sy,x2,y2,z1,z2)-x[3])/sy;
3337 Double_t f31=(F3(x1,y1,x2,y2,z1+sz,z2)-x[3])/sz;
3338 Double_t f32=(F3(x1,y1,x2,y2+sy,z1,z2)-x[3])/sy;
3339 Double_t f34=(F3(x1,y1,x2,y2,z1,z2+sz)-x[3])/sz;
3343 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3344 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3345 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3346 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3347 c[13]=f30*sy1*f40+f32*sy2*f42;
3348 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3350 // if (!BuildSeed(kr1[is],kcl,0,x1,x2,x3,x,c)) continue;
3352 index=kr1.GetIndex(is);
3353 seed->~AliTPCseed();
3354 AliTPCseed *track=new(seed) AliTPCseed(x1, sec*alpha+shift, x, c, index);
3356 track->SetIsSeeding(kTRUE);
3359 FollowProlongation(*track, i1-7,1);
3360 if (track->GetNumberOfClusters() < track->GetNFoundable()*0.75 ||
3361 track->GetNShared()>0.6*track->GetNumberOfClusters() || ( track->GetSigmaY2()+ track->GetSigmaZ2())>0.6){
3363 seed->~AliTPCseed();
3369 FollowProlongation(*track, i2,1);
3370 track->SetBConstrain(0);
3371 track->SetLastPoint(i1+fInnerSec->GetNRows()); // first cluster in track position
3372 track->SetFirstPoint(track->GetLastPoint());
3374 if (track->GetNumberOfClusters()<(i1-i2)*0.5 ||
3375 track->GetNumberOfClusters()<track->GetNFoundable()*0.7 ||
3376 track->GetNShared()>2. || track->GetChi2()/track->GetNumberOfClusters()>6 || ( track->GetSigmaY2()+ track->GetSigmaZ2())>0.5 ) {
3378 seed->~AliTPCseed();
3383 FollowProlongation(*track, TMath::Max(i2-10,0),1);
3384 AliTPCseed * track2 = MakeSeed(track,0.2,0.5,0.9);
3385 FollowProlongation(*track2, i2,1);
3386 track2->SetBConstrain(kFALSE);
3387 track2->SetSeedType(4);
3388 arr->AddLast(track2);
3390 seed->~AliTPCseed();
3394 //arr->AddLast(track);
3395 //seed = new AliTPCseed;
3401 Info("MakeSeeds5","\nSeeding statiistic:\t%d\t%d\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin,nout1,nout2,nout3);
3407 //_____________________________________________________________________________
3408 void AliTPCtrackerMI::MakeSeeds2(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, Float_t */*cuts[4]*/,
3409 Float_t deltay, Bool_t /*bconstrain*/) {
3410 //-----------------------------------------------------------------
3411 // This function creates track seeds - without vertex constraint
3412 //-----------------------------------------------------------------
3413 // cuts[0] - fP4 cut - not applied
3414 // cuts[1] - tan(phi) cut
3415 // cuts[2] - zvertex cut - not applied
3416 // cuts[3] - fP3 cut
3426 Double_t alpha=fOuterSec->GetAlpha(), shift=fOuterSec->GetAlphaShift();
3427 // Double_t cs=cos(alpha), sn=sin(alpha);
3428 Int_t row0 = (i1+i2)/2;
3429 Int_t drow = (i1-i2)/2;
3430 const AliTPCtrackerRow& kr0=fSectors[sec][row0];
3431 AliTPCtrackerRow * kr=0;
3433 AliTPCpolyTrack polytrack;
3434 Int_t nclusters=fSectors[sec][row0];
3435 AliTPCseed * seed = new AliTPCseed;
3440 for (Int_t is=0; is < nclusters; is++) { //LOOP over clusters
3442 Int_t nfoundable =0;
3443 for (Int_t iter =1; iter<2; iter++){ //iterations
3444 const AliTPCtrackerRow& krm=fSectors[sec][row0-iter];
3445 const AliTPCtrackerRow& krp=fSectors[sec][row0+iter];
3446 const AliTPCclusterMI * cl= kr0[is];
3448 if (cl->IsUsed(10)) {
3454 Double_t x = kr0.GetX();
3455 // Initialization of the polytrack
3460 Double_t y0= cl->GetY();
3461 Double_t z0= cl->GetZ();
3465 Double_t ymax = fSectors->GetMaxY(row0)-kr0.GetDeadZone()-1.5;
3466 if (deltay>0 && TMath::Abs(ymax-TMath::Abs(y0))> deltay ) continue; // seed only at the edge
3468 erry = (0.5)*cl->GetSigmaY2()/TMath::Sqrt(cl->GetQ())*6;
3469 errz = (0.5)*cl->GetSigmaZ2()/TMath::Sqrt(cl->GetQ())*6;
3470 polytrack.AddPoint(x,y0,z0,erry, errz);
3473 if (cl->IsUsed(10)) sumused++;
3476 Float_t roady = (5*TMath::Sqrt(cl->GetSigmaY2()+0.2)+1.)*iter;
3477 Float_t roadz = (5*TMath::Sqrt(cl->GetSigmaZ2()+0.2)+1.)*iter;
3480 AliTPCclusterMI * cl1 = krm.FindNearest(y0,z0,roady,roadz);
3481 if (cl1 && TMath::Abs(ymax-TMath::Abs(y0))) {
3482 erry = (0.5)*cl1->GetSigmaY2()/TMath::Sqrt(cl1->GetQ())*3;
3483 errz = (0.5)*cl1->GetSigmaZ2()/TMath::Sqrt(cl1->GetQ())*3;
3484 if (cl1->IsUsed(10)) sumused++;
3485 polytrack.AddPoint(x,cl1->GetY(),cl1->GetZ(),erry,errz);
3489 AliTPCclusterMI * cl2 = krp.FindNearest(y0,z0,roady,roadz);
3491 erry = (0.5)*cl2->GetSigmaY2()/TMath::Sqrt(cl2->GetQ())*3;
3492 errz = (0.5)*cl2->GetSigmaZ2()/TMath::Sqrt(cl2->GetQ())*3;
3493 if (cl2->IsUsed(10)) sumused++;
3494 polytrack.AddPoint(x,cl2->GetY(),cl2->GetZ(),erry,errz);
3497 if (sumused>0) continue;
3499 polytrack.UpdateParameters();
3505 nfoundable = polytrack.GetN();
3506 nfound = nfoundable;
3508 for (Int_t ddrow = iter+1; ddrow<drow;ddrow++){
3509 Float_t maxdist = 0.8*(1.+3./(ddrow));
3510 for (Int_t delta = -1;delta<=1;delta+=2){
3511 Int_t row = row0+ddrow*delta;
3512 kr = &(fSectors[sec][row]);
3513 Double_t xn = kr->GetX();
3514 Double_t ymax1 = fSectors->GetMaxY(row)-kr->GetDeadZone()-1.5;
3515 polytrack.GetFitPoint(xn,yn,zn);
3516 if (TMath::Abs(yn)>ymax1) continue;
3518 AliTPCclusterMI * cln = kr->FindNearest(yn,zn,roady,roadz);
3520 Float_t dist = TMath::Sqrt( (yn-cln->GetY())*(yn-cln->GetY())+(zn-cln->GetZ())*(zn-cln->GetZ()));
3523 erry = (dist+0.3)*cln->GetSigmaY2()/TMath::Sqrt(cln->GetQ())*(1.+1./(ddrow));
3524 errz = (dist+0.3)*cln->GetSigmaZ2()/TMath::Sqrt(cln->GetQ())*(1.+1./(ddrow));
3525 if (cln->IsUsed(10)) {
3526 // printf("used\n");
3534 polytrack.AddPoint(xn,cln->GetY(),cln->GetZ(),erry, errz);
3539 if ( (sumused>3) || (sumused>0.5*nfound) || (nfound<0.6*nfoundable)) break;
3540 polytrack.UpdateParameters();
3543 if ( (sumused>3) || (sumused>0.5*nfound)) {
3544 //printf("sumused %d\n",sumused);
3549 polytrack.GetFitDerivation(kr0.GetX(),dy,dz);
3550 AliTPCpolyTrack track2;
3552 polytrack.Refit(track2,0.5+TMath::Abs(dy)*0.3,0.4+TMath::Abs(dz)*0.3);
3553 if (track2.GetN()<0.5*nfoundable) continue;
3556 if ((nfound>0.6*nfoundable) &&( nfoundable>0.4*(i1-i2))) {
3558 // test seed with and without constrain
3559 for (Int_t constrain=0; constrain<=0;constrain++){
3560 // add polytrack candidate
3562 Double_t x[5], c[15];
3563 Double_t x1,x2,x3,y1,y2,y3,z1,z2,z3;
3564 track2.GetBoundaries(x3,x1);
3566 track2.GetFitPoint(x1,y1,z1);
3567 track2.GetFitPoint(x2,y2,z2);
3568 track2.GetFitPoint(x3,y3,z3);
3570 //is track pointing to the vertex ?
3573 polytrack.GetFitPoint(x0,y0,z0);
3586 x[4]=F1(x1,y1,x2,y2,x3,y3);
3588 // if (TMath::Abs(x[4]) >= cuts[0]) continue; //
3589 x[2]=F2(x1,y1,x2,y2,x3,y3);
3591 //if (TMath::Abs(x[4]*x1-x[2]) >= cuts[1]) continue;
3592 //x[3]=F3(x1,y1,x2,y2,z1,z2);
3593 x[3]=F3n(x1,y1,x3,y3,z1,z3,x[4]);
3594 //if (TMath::Abs(x[3]) > cuts[3]) continue;
3597 Double_t sy =0.1, sz =0.1;
3598 Double_t sy1=0.02, sz1=0.02;
3599 Double_t sy2=0.02, sz2=0.02;
3603 sy3=25000*x[4]*x[4]+0.1, sy=0.1, sz=0.1;
3606 Double_t f40=(F1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
3607 Double_t f42=(F1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
3608 Double_t f43=(F1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
3609 Double_t f20=(F2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
3610 Double_t f22=(F2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
3611 Double_t f23=(F2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
3613 Double_t f30=(F3(x1,y1+sy,x3,y3,z1,z3)-x[3])/sy;
3614 Double_t f31=(F3(x1,y1,x3,y3,z1+sz,z3)-x[3])/sz;
3615 Double_t f32=(F3(x1,y1,x3,y3+sy,z1,z3)-x[3])/sy;
3616 Double_t f34=(F3(x1,y1,x3,y3,z1,z3+sz)-x[3])/sz;
3621 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3622 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3623 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3624 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3625 c[13]=f30*sy1*f40+f32*sy2*f42;
3626 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3628 //Int_t row1 = fSectors->GetRowNumber(x1);
3629 Int_t row1 = GetRowNumber(x1);
3633 seed->~AliTPCseed();
3634 AliTPCseed *track=new(seed) AliTPCseed(x1,sec*alpha+shift,x,c,index);
3635 track->SetIsSeeding(kTRUE);
3636 Int_t rc=FollowProlongation(*track, i2);
3637 if (constrain) track->SetBConstrain(1);
3639 track->SetBConstrain(0);
3640 track->SetLastPoint(row1+fInnerSec->GetNRows()); // first cluster in track position
3641 track->SetFirstPoint(track->GetLastPoint());
3643 if (rc==0 || track->GetNumberOfClusters()<(i1-i2)*0.5 ||
3644 track->GetNumberOfClusters() < track->GetNFoundable()*0.6 ||
3645 track->GetNShared()>0.4*track->GetNumberOfClusters()) {
3648 seed->~AliTPCseed();
3651 arr->AddLast(track);
3652 seed = new AliTPCseed;
3656 } // if accepted seed
3659 Info("MakeSeeds2","\nSeeding statiistic:\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin3);
3665 AliTPCseed *AliTPCtrackerMI::MakeSeed(AliTPCseed *track, Float_t r0, Float_t r1, Float_t r2)
3669 //reseed using track points
3670 Int_t p0 = int(r0*track->GetNumberOfClusters()); // point 0
3671 Int_t p1 = int(r1*track->GetNumberOfClusters());
3672 Int_t p2 = int(r2*track->GetNumberOfClusters()); // last point
3674 Double_t x0[3],x1[3],x2[3];
3675 for (Int_t i=0;i<3;i++){
3681 // find track position at given ratio of the length
3682 Int_t sec0=0, sec1=0, sec2=0;
3685 for (Int_t i=0;i<160;i++){
3686 if (track->GetClusterPointer(i)){
3688 AliTPCTrackerPoint *trpoint =track->GetTrackPoint(i);
3689 if ( (index<p0) || x0[0]<0 ){
3690 if (trpoint->GetX()>1){
3691 clindex = track->GetClusterIndex2(i);
3693 x0[0] = trpoint->GetX();
3694 x0[1] = trpoint->GetY();
3695 x0[2] = trpoint->GetZ();
3696 sec0 = ((clindex&0xff000000)>>24)%18;
3701 if ( (index<p1) &&(trpoint->GetX()>1)){
3702 clindex = track->GetClusterIndex2(i);
3704 x1[0] = trpoint->GetX();
3705 x1[1] = trpoint->GetY();
3706 x1[2] = trpoint->GetZ();
3707 sec1 = ((clindex&0xff000000)>>24)%18;
3710 if ( (index<p2) &&(trpoint->GetX()>1)){
3711 clindex = track->GetClusterIndex2(i);
3713 x2[0] = trpoint->GetX();
3714 x2[1] = trpoint->GetY();
3715 x2[2] = trpoint->GetZ();
3716 sec2 = ((clindex&0xff000000)>>24)%18;
3723 Double_t alpha, cs,sn, xx2,yy2;
3725 alpha = (sec1-sec2)*fSectors->GetAlpha();
3726 cs = TMath::Cos(alpha);
3727 sn = TMath::Sin(alpha);
3728 xx2= x1[0]*cs-x1[1]*sn;
3729 yy2= x1[0]*sn+x1[1]*cs;
3733 alpha = (sec0-sec2)*fSectors->GetAlpha();
3734 cs = TMath::Cos(alpha);
3735 sn = TMath::Sin(alpha);
3736 xx2= x0[0]*cs-x0[1]*sn;
3737 yy2= x0[0]*sn+x0[1]*cs;
3743 Double_t x[5],c[15];
3747 x[4]=F1(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]);
3748 // if (x[4]>1) return 0;
3749 x[2]=F2(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]);
3750 x[3]=F3n(x2[0],x2[1],x0[0],x0[1],x2[2],x0[2],x[4]);
3751 //if (TMath::Abs(x[3]) > 2.2) return 0;
3752 //if (TMath::Abs(x[2]) > 1.99) return 0;
3754 Double_t sy =0.1, sz =0.1;
3756 Double_t sy1=0.02+track->GetSigmaY2(), sz1=0.02+track->GetSigmaZ2();
3757 Double_t sy2=0.01+track->GetSigmaY2(), sz2=0.01+track->GetSigmaZ2();
3758 Double_t sy3=0.01+track->GetSigmaY2();
3760 Double_t f40=(F1(x2[0],x2[1]+sy,x1[0],x1[1],x0[0],x0[1])-x[4])/sy;
3761 Double_t f42=(F1(x2[0],x2[1],x1[0],x1[1]+sy,x0[0],x0[1])-x[4])/sy;
3762 Double_t f43=(F1(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]+sy)-x[4])/sy;
3763 Double_t f20=(F2(x2[0],x2[1]+sy,x1[0],x1[1],x0[0],x0[1])-x[2])/sy;
3764 Double_t f22=(F2(x2[0],x2[1],x1[0],x1[1]+sy,x0[0],x0[1])-x[2])/sy;
3765 Double_t f23=(F2(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]+sy)-x[2])/sy;
3767 Double_t f30=(F3(x2[0],x2[1]+sy,x0[0],x0[1],x2[2],x0[2])-x[3])/sy;
3768 Double_t f31=(F3(x2[0],x2[1],x0[0],x0[1],x2[2]+sz,x0[2])-x[3])/sz;
3769 Double_t f32=(F3(x2[0],x2[1],x0[0],x0[1]+sy,x2[2],x0[2])-x[3])/sy;
3770 Double_t f34=(F3(x2[0],x2[1],x0[0],x0[1],x2[2],x0[2]+sz)-x[3])/sz;
3775 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3776 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3777 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3778 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3779 c[13]=f30*sy1*f40+f32*sy2*f42;
3780 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3782 // Int_t row1 = fSectors->GetRowNumber(x2[0]);
3783 AliTPCseed *seed=new AliTPCseed(x2[0], sec2*fSectors->GetAlpha()+fSectors->GetAlphaShift(), x, c, 0);
3784 // Double_t y0,z0,y1,z1, y2,z2;
3785 //seed->GetProlongation(x0[0],y0,z0);
3786 // seed->GetProlongation(x1[0],y1,z1);
3787 //seed->GetProlongation(x2[0],y2,z2);
3789 seed->SetLastPoint(pp2);
3790 seed->SetFirstPoint(pp2);
3797 AliTPCseed *AliTPCtrackerMI::ReSeed(const AliTPCseed *track, Float_t r0, Float_t r1, Float_t r2)
3801 //reseed using founded clusters
3803 // Find the number of clusters
3804 Int_t nclusters = 0;
3805 for (Int_t irow=0;irow<160;irow++){
3806 if (track->GetClusterIndex(irow)>0) nclusters++;
3810 ipos[0] = TMath::Max(int(r0*nclusters),0); // point 0 cluster
3811 ipos[1] = TMath::Min(int(r1*nclusters),nclusters-1); //
3812 ipos[2] = TMath::Min(int(r2*nclusters),nclusters-1); // last point
3816 Int_t row[3],sec[3]={0,0,0};
3818 // find track row position at given ratio of the length
3820 for (Int_t irow=0;irow<160;irow++){
3821 if (track->GetClusterIndex2(irow)<0) continue;
3823 for (Int_t ipoint=0;ipoint<3;ipoint++){
3824 if (index<=ipos[ipoint]) row[ipoint] = irow;
3828 //Get cluster and sector position
3829 for (Int_t ipoint=0;ipoint<3;ipoint++){
3830 Int_t clindex = track->GetClusterIndex2(row[ipoint]);
3831 AliTPCclusterMI * cl = GetClusterMI(clindex);
3834 // AliTPCclusterMI * cl = GetClusterMI(clindex);
3837 sec[ipoint] = ((clindex&0xff000000)>>24)%18;
3838 xyz[ipoint][0] = GetXrow(row[ipoint]);
3839 xyz[ipoint][1] = cl->GetY();
3840 xyz[ipoint][2] = cl->GetZ();
3844 // Calculate seed state vector and covariance matrix
3846 Double_t alpha, cs,sn, xx2,yy2;
3848 alpha = (sec[1]-sec[2])*fSectors->GetAlpha();
3849 cs = TMath::Cos(alpha);
3850 sn = TMath::Sin(alpha);
3851 xx2= xyz[1][0]*cs-xyz[1][1]*sn;
3852 yy2= xyz[1][0]*sn+xyz[1][1]*cs;
3856 alpha = (sec[0]-sec[2])*fSectors->GetAlpha();
3857 cs = TMath::Cos(alpha);
3858 sn = TMath::Sin(alpha);
3859 xx2= xyz[0][0]*cs-xyz[0][1]*sn;
3860 yy2= xyz[0][0]*sn+xyz[0][1]*cs;
3866 Double_t x[5],c[15];
3870 x[4]=F1(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
3871 x[2]=F2(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
3872 x[3]=F3n(xyz[2][0],xyz[2][1],xyz[0][0],xyz[0][1],xyz[2][2],xyz[0][2],x[4]);
3874 Double_t sy =0.1, sz =0.1;
3876 Double_t sy1=0.2, sz1=0.2;
3877 Double_t sy2=0.2, sz2=0.2;
3880 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;
3881 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;
3882 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;
3883 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;
3884 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;
3885 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;
3887 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;
3888 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;
3889 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;
3890 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;
3895 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3896 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3897 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3898 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3899 c[13]=f30*sy1*f40+f32*sy2*f42;
3900 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3902 // Int_t row1 = fSectors->GetRowNumber(xyz[2][0]);
3903 AliTPCseed *seed=new AliTPCseed(xyz[2][0], sec[2]*fSectors->GetAlpha()+fSectors->GetAlphaShift(), x, c, 0);
3904 seed->SetLastPoint(row[2]);
3905 seed->SetFirstPoint(row[2]);
3910 AliTPCseed *AliTPCtrackerMI::ReSeed(AliTPCseed *track,Int_t r0, Bool_t forward)
3914 //reseed using founded clusters
3917 Int_t row[3]={0,0,0};
3918 Int_t sec[3]={0,0,0};
3920 // forward direction
3922 for (Int_t irow=r0;irow<160;irow++){
3923 if (track->GetClusterIndex(irow)>0){
3928 for (Int_t irow=160;irow>r0;irow--){
3929 if (track->GetClusterIndex(irow)>0){
3934 for (Int_t irow=row[2]-15;irow>row[0];irow--){
3935 if (track->GetClusterIndex(irow)>0){
3943 for (Int_t irow=0;irow<r0;irow++){
3944 if (track->GetClusterIndex(irow)>0){
3949 for (Int_t irow=r0;irow>0;irow--){
3950 if (track->GetClusterIndex(irow)>0){
3955 for (Int_t irow=row[2]-15;irow>row[0];irow--){
3956 if (track->GetClusterIndex(irow)>0){
3963 if ((row[2]-row[0])<20) return 0;
3964 if (row[1]==0) return 0;
3967 //Get cluster and sector position
3968 for (Int_t ipoint=0;ipoint<3;ipoint++){
3969 Int_t clindex = track->GetClusterIndex2(row[ipoint]);
3970 AliTPCclusterMI * cl = GetClusterMI(clindex);
3973 // AliTPCclusterMI * cl = GetClusterMI(clindex);
3976 sec[ipoint] = ((clindex&0xff000000)>>24)%18;
3977 xyz[ipoint][0] = GetXrow(row[ipoint]);
3978 AliTPCTrackerPoint * point = track->GetTrackPoint(row[ipoint]);
3979 if (point&&ipoint<2){
3981 xyz[ipoint][1] = point->GetY();
3982 xyz[ipoint][2] = point->GetZ();
3985 xyz[ipoint][1] = cl->GetY();
3986 xyz[ipoint][2] = cl->GetZ();
3993 // Calculate seed state vector and covariance matrix
3995 Double_t alpha, cs,sn, xx2,yy2;
3997 alpha = (sec[1]-sec[2])*fSectors->GetAlpha();
3998 cs = TMath::Cos(alpha);
3999 sn = TMath::Sin(alpha);
4000 xx2= xyz[1][0]*cs-xyz[1][1]*sn;
4001 yy2= xyz[1][0]*sn+xyz[1][1]*cs;
4005 alpha = (sec[0]-sec[2])*fSectors->GetAlpha();
4006 cs = TMath::Cos(alpha);
4007 sn = TMath::Sin(alpha);
4008 xx2= xyz[0][0]*cs-xyz[0][1]*sn;
4009 yy2= xyz[0][0]*sn+xyz[0][1]*cs;
4015 Double_t x[5],c[15];
4019 x[4]=F1(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
4020 x[2]=F2(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
4021 x[3]=F3n(xyz[2][0],xyz[2][1],xyz[0][0],xyz[0][1],xyz[2][2],xyz[0][2],x[4]);
4023 Double_t sy =0.1, sz =0.1;
4025 Double_t sy1=0.2, sz1=0.2;
4026 Double_t sy2=0.2, sz2=0.2;
4029 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;
4030 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;
4031 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;
4032 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;
4033 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;
4034 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;
4036 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;
4037 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;
4038 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;
4039 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;
4044 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
4045 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
4046 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
4047 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
4048 c[13]=f30*sy1*f40+f32*sy2*f42;
4049 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
4051 // Int_t row1 = fSectors->GetRowNumber(xyz[2][0]);
4052 AliTPCseed *seed=new AliTPCseed(xyz[2][0], sec[2]*fSectors->GetAlpha()+fSectors->GetAlphaShift(), x, c, 0);
4053 seed->SetLastPoint(row[2]);
4054 seed->SetFirstPoint(row[2]);
4055 for (Int_t i=row[0];i<row[2];i++){
4056 seed->SetClusterIndex(i, track->GetClusterIndex(i));
4064 void AliTPCtrackerMI::FindMultiMC(const TObjArray * array, AliESDEvent */*esd*/, Int_t iter)
4067 // find multi tracks - THIS FUNCTION IS ONLY FOR DEBUG PURPOSES
4069 // Use AliTPCReconstructor::StreamLevel()>2 if you want to tune parameters - cuts
4071 // Two reasons to have multiple find tracks
4072 // 1. Curling tracks can be find more than once
4073 // 2. Splitted tracks
4074 // a.) Multiple seeding to increase tracking efficiency - (~ 100% reached)
4075 // b.) Edge effect on the sector boundaries
4078 // Algorithm done in 2 phases - because of CPU consumption
4079 // it is n^2 algorithm - for lead-lead 20000x20000 combination are investigated
4081 // Algorihm for curling tracks sign:
4082 // 1 phase -makes a very rough fast cuts to minimize combinatorics
4083 // a.) opposite sign
4084 // b.) one of the tracks - not pointing to the primary vertex -
4085 // c.) delta tan(theta)
4087 // 2 phase - calculates DCA between tracks - time consument
4092 // General cuts - for splitted tracks and for curling tracks
4094 const Float_t kMaxdPhi = 0.2; // maximal distance in phi
4096 // Curling tracks cuts
4101 Int_t nentries = array->GetEntriesFast();
4102 AliHelix *helixes = new AliHelix[nentries];
4103 Float_t *xm = new Float_t[nentries];
4104 Float_t *dz0 = new Float_t[nentries];
4105 Float_t *dz1 = new Float_t[nentries];
4111 // Find track COG in x direction - point with best defined parameters
4113 for (Int_t i=0;i<nentries;i++){
4114 AliTPCseed* track = (AliTPCseed*)array->At(i);
4115 if (!track) continue;
4116 track->SetCircular(0);
4117 new (&helixes[i]) AliHelix(*track);
4121 track->GetDZ(GetX(),GetY(),GetZ(),GetBz(),dz);
4124 for (Int_t icl=0; icl<160; icl++){
4125 AliTPCclusterMI * cl = track->GetClusterPointer(icl);
4131 if (ncl>0) xm[i]/=Float_t(ncl);
4134 for (Int_t i0=0;i0<nentries;i0++){
4135 AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
4136 if (!track0) continue;
4137 Float_t xc0 = helixes[i0].GetHelix(6);
4138 Float_t yc0 = helixes[i0].GetHelix(7);
4139 Float_t r0 = helixes[i0].GetHelix(8);
4140 Float_t rc0 = TMath::Sqrt(xc0*xc0+yc0*yc0);
4141 Float_t fi0 = TMath::ATan2(yc0,xc0);
4143 for (Int_t i1=i0+1;i1<nentries;i1++){
4144 AliTPCseed * track1 = (AliTPCseed*)array->At(i1);
4145 if (!track1) continue;
4146 Int_t lab0=track0->GetLabel();
4147 Int_t lab1=track1->GetLabel();
4148 if (TMath::Abs(lab0)!=TMath::Abs(lab1)) continue;
4150 Float_t xc1 = helixes[i1].GetHelix(6);
4151 Float_t yc1 = helixes[i1].GetHelix(7);
4152 Float_t r1 = helixes[i1].GetHelix(8);
4153 Float_t rc1 = TMath::Sqrt(xc1*xc1+yc1*yc1);
4154 Float_t fi1 = TMath::ATan2(yc1,xc1);
4156 Float_t dfi = fi0-fi1;
4159 if (dfi>1.5*TMath::Pi()) dfi-=TMath::Pi(); // take care about edge effect
4160 if (dfi<-1.5*TMath::Pi()) dfi+=TMath::Pi(); //
4161 if (TMath::Abs(dfi)>kMaxdPhi&&helixes[i0].GetHelix(4)*helixes[i1].GetHelix(4)<0){
4163 // if short tracks with undefined sign
4164 fi1 = -TMath::ATan2(yc1,-xc1);
4167 Float_t dtheta = TMath::Abs(track0->GetTgl()-track1->GetTgl())<TMath::Abs(track0->GetTgl()+track1->GetTgl())? track0->GetTgl()-track1->GetTgl():track0->GetTgl()+track1->GetTgl();
4170 // debug stream to tune "fast cuts"
4172 Double_t dist[3]; // distance at X
4173 Double_t mdist[3]={0,0,0}; // mean distance X+-40cm
4174 track0->GetDistance(track1,0.5*(xm[i0]+xm[i1])-40.,dist,AliTracker::GetBz());
4175 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[i]);
4176 track0->GetDistance(track1,0.5*(xm[i0]+xm[i1])+40.,dist,AliTracker::GetBz());
4177 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[i]);
4178 track0->GetDistance(track1,0.5*(xm[i0]+xm[i1]),dist,AliTracker::GetBz());
4179 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[i]);
4180 for (Int_t i=0;i<3;i++) mdist[i]*=0.33333;
4184 for (Int_t icl=0; icl<160; icl++){
4185 AliTPCclusterMI * cl0 = track0->GetClusterPointer(icl);
4186 AliTPCclusterMI * cl1 = track1->GetClusterPointer(icl);
4189 if (cl0==cl1) sums++;
4193 if (AliTPCReconstructor::StreamLevel()>0) {
4194 TTreeSRedirector &cstream = *fDebugStreamer;
4199 "Tr0.="<<track0<< // seed0
4200 "Tr1.="<<track1<< // seed1
4201 "h0.="<<&helixes[i0]<<
4202 "h1.="<<&helixes[i1]<<
4204 "sum="<<sum<< //the sum of rows with cl in both
4205 "sums="<<sums<< //the sum of shared clusters
4206 "xm0="<<xm[i0]<< // the center of track
4207 "xm1="<<xm[i1]<< // the x center of track
4208 // General cut variables
4209 "dfi="<<dfi<< // distance in fi angle
4210 "dtheta="<<dtheta<< // distance int theta angle
4216 "dist0="<<dist[0]<< //distance x
4217 "dist1="<<dist[1]<< //distance y
4218 "dist2="<<dist[2]<< //distance z
4219 "mdist0="<<mdist[0]<< //distance x
4220 "mdist1="<<mdist[1]<< //distance y
4221 "mdist2="<<mdist[2]<< //distance z
4235 if (AliTPCReconstructor::StreamLevel()>1) {
4236 AliInfo("Time for curling tracks removal DEBUGGING MC");
4242 void AliTPCtrackerMI::FindSplitted(TObjArray * array, AliESDEvent */*esd*/, Int_t iter)
4246 // Two reasons to have multiple find tracks
4247 // 1. Curling tracks can be find more than once
4248 // 2. Splitted tracks
4249 // a.) Multiple seeding to increase tracking efficiency - (~ 100% reached)
4250 // b.) Edge effect on the sector boundaries
4252 // This function tries to find tracks closed in the parametric space
4254 // cut logic if distance is bigger than cut continue - Do Nothing
4255 const Float_t kMaxdTheta = 0.05; // maximal distance in theta
4256 const Float_t kMaxdPhi = 0.05; // maximal deistance in phi
4257 const Float_t kdelta = 40.; //delta r to calculate track distance
4259 // const Float_t kMaxDist0 = 1.; // maximal distance 0
4260 //const Float_t kMaxDist1 = 0.3; // maximal distance 1 - cut if track in separate rows
4263 TCut csec("csec","abs(Tr0.fRelativeSector-Tr1.fRelativeSector)<2");
4264 TCut cdtheta("cdtheta","abs(dtheta)<0.05");
4269 Int_t nentries = array->GetEntriesFast();
4270 AliHelix *helixes = new AliHelix[nentries];
4271 Float_t *xm = new Float_t[nentries];
4277 //Sort tracks according quality
4279 Int_t nseed = array->GetEntriesFast();
4280 Float_t * quality = new Float_t[nseed];
4281 Int_t * indexes = new Int_t[nseed];
4282 for (Int_t i=0; i<nseed; i++) {
4283 AliTPCseed *pt=(AliTPCseed*)array->UncheckedAt(i);
4288 pt->UpdatePoints(); //select first last max dens points
4289 Float_t * points = pt->GetPoints();
4290 if (points[3]<0.8) quality[i] =-1;
4291 quality[i] = (points[2]-points[0])+pt->GetNumberOfClusters();
4292 //prefer high momenta tracks if overlaps
4293 quality[i] *= TMath::Sqrt(TMath::Abs(pt->Pt())+0.5);
4295 TMath::Sort(nseed,quality,indexes);
4299 // Find track COG in x direction - point with best defined parameters
4301 for (Int_t i=0;i<nentries;i++){
4302 AliTPCseed* track = (AliTPCseed*)array->At(i);
4303 if (!track) continue;
4304 track->SetCircular(0);
4305 new (&helixes[i]) AliHelix(*track);
4308 for (Int_t icl=0; icl<160; icl++){
4309 AliTPCclusterMI * cl = track->GetClusterPointer(icl);
4315 if (ncl>0) xm[i]/=Float_t(ncl);
4318 for (Int_t is0=0;is0<nentries;is0++){
4319 Int_t i0 = indexes[is0];
4320 AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
4321 if (!track0) continue;
4322 Float_t xc0 = helixes[i0].GetHelix(6);
4323 Float_t yc0 = helixes[i0].GetHelix(7);
4324 Float_t fi0 = TMath::ATan2(yc0,xc0);
4326 for (Int_t is1=is0+1;is1<nentries;is1++){
4327 Int_t i1 = indexes[is1];
4328 AliTPCseed * track1 = (AliTPCseed*)array->At(i1);
4329 if (!track1) continue;
4331 Int_t dsec = TMath::Abs((track0->GetRelativeSector()%18)-(track1->GetRelativeSector()%18)); // sector distance
4332 if (dsec>1 && dsec<17) continue;
4334 if (track1->GetKinkIndexes()[0] == -track0->GetKinkIndexes()[0]) continue;
4336 Float_t dtheta = TMath::Abs(track0->GetTgl()-track1->GetTgl())<TMath::Abs(track0->GetTgl()+track1->GetTgl())? track0->GetTgl()-track1->GetTgl():track0->GetTgl()+track1->GetTgl();
4337 if (TMath::Abs(dtheta)>kMaxdTheta) continue;
4339 Float_t xc1 = helixes[i1].GetHelix(6);
4340 Float_t yc1 = helixes[i1].GetHelix(7);
4341 Float_t fi1 = TMath::ATan2(yc1,xc1);
4343 Float_t dfi = fi0-fi1;
4344 if (dfi>TMath::Pi()) dfi-=TMath::TwoPi(); // take care about edge effect
4345 if (dfi<-TMath::Pi()) dfi+=TMath::TwoPi(); //
4346 if (TMath::Abs(dfi)>kMaxdPhi&&helixes[i0].GetHelix(4)*helixes[i1].GetHelix(4)<0){
4348 // if short tracks with undefined sign
4349 fi1 = -TMath::ATan2(yc1,-xc1);
4351 if (dfi>TMath::Pi()) dfi-=TMath::TwoPi(); // take care about edge effect
4352 if (dfi<-TMath::Pi()) dfi+=TMath::TwoPi(); //
4354 if (TMath::Abs(dfi)>kMaxdPhi) continue;
4361 for (Int_t icl=0; icl<160; icl++){
4362 Int_t index0=track0->GetClusterIndex2(icl);
4363 Int_t index1=track1->GetClusterIndex2(icl);
4364 Bool_t used0 = (index0>0 && !(index0&0x8000));
4365 Bool_t used1 = (index1>0 && !(index1&0x8000));
4367 if (used0) sum0++; // used cluster0
4368 if (used1) sum1++; // used clusters1
4369 if (used0&&used1) sum++;
4370 if (index0==index1 && used0 && used1) sums++;
4374 if (sums<10) continue;
4375 if (sum<40) continue;
4376 if (sums/Float_t(TMath::Min(sum0,sum1))<0.5) continue;
4378 Double_t dist[5][4]; // distance at X
4379 Double_t mdist[4]={0,0,0,0}; // mean distance on range +-delta
4383 track0->GetDistance(track1,xm[i0],dist[0],AliTracker::GetBz());
4384 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[0][i]);
4385 track0->GetDistance(track1,xm[i1],dist[1],AliTracker::GetBz());
4386 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[1][i]);
4388 track0->GetDistance(track1,TMath::Min(xm[i1],xm[i0])-kdelta,dist[2],AliTracker::GetBz());
4389 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[2][i]);
4390 track0->GetDistance(track1,TMath::Max(xm[i1],xm[i0])+kdelta,dist[3],AliTracker::GetBz());
4391 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[3][i]);
4393 track0->GetDistance(track1,(xm[i1]+xm[i0])*0.5,dist[4],AliTracker::GetBz());
4394 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[4][i]);
4395 for (Int_t i=0;i<3;i++) mdist[i]*=0.2;
4398 Int_t lab0=track0->GetLabel();
4399 Int_t lab1=track1->GetLabel();
4400 if( AliTPCReconstructor::StreamLevel()>5){
4401 TTreeSRedirector &cstream = *fDebugStreamer;
4402 cstream<<"Splitted"<<
4406 "Tr0.="<<track0<< // seed0
4407 "Tr1.="<<track1<< // seed1
4408 "h0.="<<&helixes[i0]<<
4409 "h1.="<<&helixes[i1]<<
4411 "sum="<<sum<< //the sum of rows with cl in both
4412 "sum0="<<sum0<< //the sum of rows with cl in 0 track
4413 "sum1="<<sum1<< //the sum of rows with cl in 1 track
4414 "sums="<<sums<< //the sum of shared clusters
4415 "xm0="<<xm[i0]<< // the center of track
4416 "xm1="<<xm[i1]<< // the x center of track
4417 // General cut variables
4418 "dfi="<<dfi<< // distance in fi angle
4419 "dtheta="<<dtheta<< // distance int theta angle
4422 "dist0="<<dist[4][0]<< //distance x
4423 "dist1="<<dist[4][1]<< //distance y
4424 "dist2="<<dist[4][2]<< //distance z
4425 "mdist0="<<mdist[0]<< //distance x
4426 "mdist1="<<mdist[1]<< //distance y
4427 "mdist2="<<mdist[2]<< //distance z
4430 delete array->RemoveAt(i1);
4437 AliInfo("Time for splitted tracks removal");
4443 void AliTPCtrackerMI::FindCurling(const TObjArray * array, AliESDEvent */*esd*/, Int_t iter)
4446 // find Curling tracks
4447 // Use AliTPCReconstructor::StreamLevel()>1 if you want to tune parameters - cuts
4450 // Algorithm done in 2 phases - because of CPU consumption
4451 // it is n^2 algorithm - for lead-lead 20000x20000 combination are investigated
4452 // see detal in MC part what can be used to cut
4456 const Float_t kMaxC = 400; // maximal curvature to of the track
4457 const Float_t kMaxdTheta = 0.15; // maximal distance in theta
4458 const Float_t kMaxdPhi = 0.15; // maximal distance in phi
4459 const Float_t kPtRatio = 0.3; // ratio between pt
4460 const Float_t kMinDCAR = 2.; // distance to the primary vertex in r - see cpipe cut
4463 // Curling tracks cuts
4466 const Float_t kMaxDeltaRMax = 40; // distance in outer radius
4467 const Float_t kMaxDeltaRMin = 5.; // distance in lower radius - see cpipe cut
4468 const Float_t kMinAngle = 2.9; // angle between tracks
4469 const Float_t kMaxDist = 5; // biggest distance
4471 // The cuts can be tuned using the "MC information stored in Multi tree ==> see FindMultiMC
4474 TCut csign("csign","Tr0.fP[4]*Tr1.fP[4]<0"); //opposite sign
4475 TCut cmax("cmax","abs(Tr0.GetC())>1/400");
4476 TCut cda("cda","sqrt(dtheta^2+dfi^2)<0.15");
4477 TCut ccratio("ccratio","abs((Tr0.fP[4]+Tr1.fP[4])/(abs(Tr0.fP[4])+abs(Tr1.fP[4])))<0.3");
4478 TCut cpipe("cpipe", "min(abs(r0-rc0),abs(r1-rc1))>5");
4480 TCut cdrmax("cdrmax","abs(abs(rc0+r0)-abs(rc1+r1))<40")
4481 TCut cdrmin("cdrmin","abs(abs(rc0+r0)-abs(rc1+r1))<10")
4483 Multi->Draw("dfi","iter==0"+csign+cmax+cda+ccratio); ~94% of curling tracks fulfill
4484 Multi->Draw("min(abs(r0-rc0),abs(r1-rc1))","iter==0&&abs(lab1)==abs(lab0)"+csign+cmax+cda+ccratio+cpipe+cdrmin+cdrmax); //80%
4486 Curling2->Draw("dfi","iter==0&&abs(lab0)==abs(lab1)"+csign+cmax+cdtheta+cdfi+ccratio)
4492 Int_t nentries = array->GetEntriesFast();
4493 AliHelix *helixes = new AliHelix[nentries];
4494 for (Int_t i=0;i<nentries;i++){
4495 AliTPCseed* track = (AliTPCseed*)array->At(i);
4496 if (!track) continue;
4497 track->SetCircular(0);
4498 new (&helixes[i]) AliHelix(*track);
4504 Double_t phase[2][2],radius[2];
4509 for (Int_t i0=0;i0<nentries;i0++){
4510 AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
4511 if (!track0) continue;
4512 if (TMath::Abs(track0->GetC())<1/kMaxC) continue;
4513 Float_t xc0 = helixes[i0].GetHelix(6);
4514 Float_t yc0 = helixes[i0].GetHelix(7);
4515 Float_t r0 = helixes[i0].GetHelix(8);
4516 Float_t rc0 = TMath::Sqrt(xc0*xc0+yc0*yc0);
4517 Float_t fi0 = TMath::ATan2(yc0,xc0);
4519 for (Int_t i1=i0+1;i1<nentries;i1++){
4520 AliTPCseed * track1 = (AliTPCseed*)array->At(i1);
4521 if (!track1) continue;
4522 if (TMath::Abs(track1->GetC())<1/kMaxC) continue;
4523 Float_t xc1 = helixes[i1].GetHelix(6);
4524 Float_t yc1 = helixes[i1].GetHelix(7);
4525 Float_t r1 = helixes[i1].GetHelix(8);
4526 Float_t rc1 = TMath::Sqrt(xc1*xc1+yc1*yc1);
4527 Float_t fi1 = TMath::ATan2(yc1,xc1);
4529 Float_t dfi = fi0-fi1;
4532 if (dfi>1.5*TMath::Pi()) dfi-=TMath::Pi(); // take care about edge effect
4533 if (dfi<-1.5*TMath::Pi()) dfi+=TMath::Pi(); //
4534 Float_t dtheta = TMath::Abs(track0->GetTgl()-track1->GetTgl())<TMath::Abs(track0->GetTgl()+track1->GetTgl())? track0->GetTgl()-track1->GetTgl():track0->GetTgl()+track1->GetTgl();
4538 if (track0->GetBConstrain()&&track1->GetBConstrain()) continue; // not constrained
4539 if (track1->GetSigned1Pt()*track0->GetSigned1Pt()>0) continue; // not the same sign
4540 if ( TMath::Abs(track1->GetTgl()+track0->GetTgl())>kMaxdTheta) continue; //distance in the Theta
4541 if ( TMath::Abs(dfi)>kMaxdPhi) continue; //distance in phi
4542 if ( TMath::Sqrt(dfi*dfi+dtheta*dtheta)>kMaxdPhi) continue; //common angular offset
4544 Float_t pt0 = track0->GetSignedPt();
4545 Float_t pt1 = track1->GetSignedPt();
4546 if ((TMath::Abs(pt0+pt1)/(TMath::Abs(pt0)+TMath::Abs(pt1)))>kPtRatio) continue;
4547 if ((iter==1) && TMath::Abs(TMath::Abs(rc0+r0)-TMath::Abs(rc1+r1))>kMaxDeltaRMax) continue;
4548 if ((iter!=1) &&TMath::Abs(TMath::Abs(rc0-r0)-TMath::Abs(rc1-r1))>kMaxDeltaRMin) continue;
4549 if (TMath::Min(TMath::Abs(rc0-r0),TMath::Abs(rc1-r1))<kMinDCAR) continue;
4552 // Now find closest approach
4556 Int_t npoints = helixes[i0].GetRPHIintersections(helixes[i1], phase, radius,10);
4557 if (npoints==0) continue;
4558 helixes[i0].GetClosestPhases(helixes[i1], phase);
4562 Double_t hangles[3];
4563 helixes[i0].Evaluate(phase[0][0],xyz0);
4564 helixes[i1].Evaluate(phase[0][1],xyz1);
4566 helixes[i0].GetAngle(phase[0][0],helixes[i1],phase[0][1],hangles);
4567 Double_t deltah[2],deltabest;
4568 if (TMath::Abs(hangles[2])<kMinAngle) continue;
4572 helixes[i0].ParabolicDCA(helixes[i1],phase[0][0],phase[0][1],radius[0],deltah[0],2);
4574 helixes[i0].ParabolicDCA(helixes[i1],phase[1][0],phase[1][1],radius[1],deltah[1],2);
4575 if (deltah[1]<deltah[0]) ibest=1;
4577 deltabest = TMath::Sqrt(deltah[ibest]);
4578 helixes[i0].Evaluate(phase[ibest][0],xyz0);
4579 helixes[i1].Evaluate(phase[ibest][1],xyz1);
4580 helixes[i0].GetAngle(phase[ibest][0],helixes[i1],phase[ibest][1],hangles);
4581 Double_t radiusbest = TMath::Sqrt(radius[ibest]);
4583 if (deltabest>kMaxDist) continue;
4584 // if (mindcar+mindcaz<40 && (TMath::Abs(hangles[2])<kMinAngle ||deltabest>3)) continue;
4585 Bool_t sign =kFALSE;
4586 if (hangles[2]>kMinAngle) sign =kTRUE;
4589 // circular[i0] = kTRUE;
4590 // circular[i1] = kTRUE;
4591 if (track0->OneOverPt()<track1->OneOverPt()){
4592 track0->SetCircular(track0->GetCircular()+1);
4593 track1->SetCircular(track1->GetCircular()+2);
4596 track1->SetCircular(track1->GetCircular()+1);
4597 track0->SetCircular(track0->GetCircular()+2);
4600 if (AliTPCReconstructor::StreamLevel()>1){
4602 //debug stream to tune "fine" cuts
4603 Int_t lab0=track0->GetLabel();
4604 Int_t lab1=track1->GetLabel();
4605 TTreeSRedirector &cstream = *fDebugStreamer;
4606 cstream<<"Curling2"<<
4622 "npoints="<<npoints<<
4623 "hangles0="<<hangles[0]<<
4624 "hangles1="<<hangles[1]<<
4625 "hangles2="<<hangles[2]<<
4628 "radius="<<radiusbest<<
4629 "deltabest="<<deltabest<<
4630 "phase0="<<phase[ibest][0]<<
4631 "phase1="<<phase[ibest][1]<<
4639 if (AliTPCReconstructor::StreamLevel()>1) {
4640 AliInfo("Time for curling tracks removal");
4649 void AliTPCtrackerMI::FindKinks(TObjArray * array, AliESDEvent *esd)
4656 TObjArray *kinks= new TObjArray(10000);
4657 // TObjArray *v0s= new TObjArray(10000);
4658 Int_t nentries = array->GetEntriesFast();
4659 AliHelix *helixes = new AliHelix[nentries];
4660 Int_t *sign = new Int_t[nentries];
4661 Int_t *nclusters = new Int_t[nentries];
4662 Float_t *alpha = new Float_t[nentries];
4663 AliKink *kink = new AliKink();
4664 Int_t * usage = new Int_t[nentries];
4665 Float_t *zm = new Float_t[nentries];
4666 Float_t *z0 = new Float_t[nentries];
4667 Float_t *fim = new Float_t[nentries];
4668 Float_t *shared = new Float_t[nentries];
4669 Bool_t *circular = new Bool_t[nentries];
4670 Float_t *dca = new Float_t[nentries];
4671 //const AliESDVertex * primvertex = esd->GetVertex();
4673 // nentries = array->GetEntriesFast();
4678 for (Int_t i=0;i<nentries;i++){
4681 AliTPCseed* track = (AliTPCseed*)array->At(i);
4682 if (!track) continue;
4683 track->SetCircular(0);
4685 track->UpdatePoints();
4686 if (( track->GetPoints()[2]- track->GetPoints()[0])>5 && track->GetPoints()[3]>0.8){
4688 nclusters[i]=track->GetNumberOfClusters();
4689 alpha[i] = track->GetAlpha();
4690 new (&helixes[i]) AliHelix(*track);
4692 helixes[i].Evaluate(0,xyz);
4693 sign[i] = (track->GetC()>0) ? -1:1;
4696 if (track->GetProlongation(x,y,z)){
4698 fim[i] = alpha[i]+TMath::ATan2(y,x);
4701 zm[i] = track->GetZ();
4705 circular[i]= kFALSE;
4706 if (track->GetProlongation(0,y,z)) z0[i] = z;
4707 dca[i] = track->GetD(0,0);
4713 Int_t ncandidates =0;
4716 Double_t phase[2][2],radius[2];
4719 // Find circling track
4721 for (Int_t i0=0;i0<nentries;i0++){
4722 AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
4723 if (!track0) continue;
4724 if (track0->GetNumberOfClusters()<40) continue;
4725 if (TMath::Abs(1./track0->GetC())>200) continue;
4726 for (Int_t i1=i0+1;i1<nentries;i1++){
4727 AliTPCseed * track1 = (AliTPCseed*)array->At(i1);
4728 if (!track1) continue;
4729 if (track1->GetNumberOfClusters()<40) continue;
4730 if ( TMath::Abs(track1->GetTgl()+track0->GetTgl())>0.1) continue;
4731 if (track0->GetBConstrain()&&track1->GetBConstrain()) continue;
4732 if (TMath::Abs(1./track1->GetC())>200) continue;
4733 if (track1->GetSigned1Pt()*track0->GetSigned1Pt()>0) continue;
4734 if (track1->GetTgl()*track0->GetTgl()>0) continue;
4735 if (TMath::Max(TMath::Abs(1./track0->GetC()),TMath::Abs(1./track1->GetC()))>190) continue;
4736 if (track0->GetBConstrain()&&track1->OneOverPt()<track0->OneOverPt()) continue; //returning - lower momenta
4737 if (track1->GetBConstrain()&&track0->OneOverPt()<track1->OneOverPt()) continue; //returning - lower momenta
4739 Float_t mindcar = TMath::Min(TMath::Abs(dca[i0]),TMath::Abs(dca[i1]));
4740 if (mindcar<5) continue;
4741 Float_t mindcaz = TMath::Min(TMath::Abs(z0[i0]-GetZ()),TMath::Abs(z0[i1]-GetZ()));
4742 if (mindcaz<5) continue;
4743 if (mindcar+mindcaz<20) continue;
4746 Float_t xc0 = helixes[i0].GetHelix(6);
4747 Float_t yc0 = helixes[i0].GetHelix(7);
4748 Float_t r0 = helixes[i0].GetHelix(8);
4749 Float_t xc1 = helixes[i1].GetHelix(6);
4750 Float_t yc1 = helixes[i1].GetHelix(7);
4751 Float_t r1 = helixes[i1].GetHelix(8);
4753 Float_t rmean = (r0+r1)*0.5;
4754 Float_t delta =TMath::Sqrt((xc1-xc0)*(xc1-xc0)+(yc1-yc0)*(yc1-yc0));
4755 //if (delta>30) continue;
4756 if (delta>rmean*0.25) continue;
4757 if (TMath::Abs(r0-r1)/rmean>0.3) continue;
4759 Int_t npoints = helixes[i0].GetRPHIintersections(helixes[i1], phase, radius,10);
4760 if (npoints==0) continue;
4761 helixes[i0].GetClosestPhases(helixes[i1], phase);
4765 Double_t hangles[3];
4766 helixes[i0].Evaluate(phase[0][0],xyz0);
4767 helixes[i1].Evaluate(phase[0][1],xyz1);
4769 helixes[i0].GetAngle(phase[0][0],helixes[i1],phase[0][1],hangles);
4770 Double_t deltah[2],deltabest;
4771 if (hangles[2]<2.8) continue;
4774 helixes[i0].ParabolicDCA(helixes[i1],phase[0][0],phase[0][1],radius[0],deltah[0],2);
4776 helixes[i0].ParabolicDCA(helixes[i1],phase[1][0],phase[1][1],radius[1],deltah[1],2);
4777 if (deltah[1]<deltah[0]) ibest=1;
4779 deltabest = TMath::Sqrt(deltah[ibest]);
4780 helixes[i0].Evaluate(phase[ibest][0],xyz0);
4781 helixes[i1].Evaluate(phase[ibest][1],xyz1);
4782 helixes[i0].GetAngle(phase[ibest][0],helixes[i1],phase[ibest][1],hangles);
4783 Double_t radiusbest = TMath::Sqrt(radius[ibest]);
4785 if (deltabest>6) continue;
4786 if (mindcar+mindcaz<40 && (hangles[2]<3.12||deltabest>3)) continue;
4787 Bool_t lsign =kFALSE;
4788 if (hangles[2]>3.06) lsign =kTRUE;
4791 circular[i0] = kTRUE;
4792 circular[i1] = kTRUE;
4793 if (track0->OneOverPt()<track1->OneOverPt()){
4794 track0->SetCircular(track0->GetCircular()+1);
4795 track1->SetCircular(track1->GetCircular()+2);
4798 track1->SetCircular(track1->GetCircular()+1);
4799 track0->SetCircular(track0->GetCircular()+2);
4802 if (lsign&&AliTPCReconstructor::StreamLevel()>1){
4804 Int_t lab0=track0->GetLabel();
4805 Int_t lab1=track1->GetLabel();
4806 TTreeSRedirector &cstream = *fDebugStreamer;
4807 cstream<<"Curling"<<
4814 "mindcar="<<mindcar<<
4815 "mindcaz="<<mindcaz<<
4818 "npoints="<<npoints<<
4819 "hangles0="<<hangles[0]<<
4820 "hangles2="<<hangles[2]<<
4825 "radius="<<radiusbest<<
4826 "deltabest="<<deltabest<<
4827 "phase0="<<phase[ibest][0]<<
4828 "phase1="<<phase[ibest][1]<<
4838 for (Int_t i =0;i<nentries;i++){
4839 if (sign[i]==0) continue;
4840 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
4847 Double_t cradius0 = 40*40;
4848 Double_t cradius1 = 270*270;
4851 Double_t cdist3=0.55;
4852 for (Int_t j =i+1;j<nentries;j++){
4854 if (sign[j]*sign[i]<1) continue;
4855 if ( (nclusters[i]+nclusters[j])>200) continue;
4856 if ( (nclusters[i]+nclusters[j])<80) continue;
4857 if ( TMath::Abs(zm[i]-zm[j])>60.) continue;
4858 if ( TMath::Abs(fim[i]-fim[j])>0.6 && TMath::Abs(fim[i]-fim[j])<5.7 ) continue;
4859 //AliTPCseed * track1 = (AliTPCseed*)array->At(j); Double_t phase[2][2],radius[2];
4860 Int_t npoints = helixes[i].GetRPHIintersections(helixes[j], phase, radius,20);
4861 if (npoints<1) continue;
4864 if (radius[0]<cradius0||radius[0]>cradius1) continue;
4867 if ( (radius[0]<cradius0||radius[0]>cradius1) && (radius[1]<cradius0||radius[1]>cradius1) ) continue;
4870 Double_t delta1=10000,delta2=10000;
4871 // cuts on the intersection radius
4872 helixes[i].LinearDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta1);
4873 if (radius[0]<20&&delta1<1) continue; //intersection at vertex
4874 if (radius[0]<10&&delta1<3) continue; //intersection at vertex
4876 helixes[i].LinearDCA(helixes[j],phase[1][0],phase[1][1],radius[1],delta2);
4877 if (radius[1]<20&&delta2<1) continue; //intersection at vertex
4878 if (radius[1]<10&&delta2<3) continue; //intersection at vertex
4881 Double_t distance1 = TMath::Min(delta1,delta2);
4882 if (distance1>cdist1) continue; // cut on DCA linear approximation
4884 npoints = helixes[i].GetRPHIintersections(helixes[j], phase, radius,20);
4885 helixes[i].ParabolicDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta1);
4886 if (radius[0]<20&&delta1<1) continue; //intersection at vertex
4887 if (radius[0]<10&&delta1<3) continue; //intersection at vertex
4890 helixes[i].ParabolicDCA(helixes[j],phase[1][0],phase[1][1],radius[1],delta2);
4891 if (radius[1]<20&&delta2<1) continue; //intersection at vertex
4892 if (radius[1]<10&&delta2<3) continue; //intersection at vertex
4894 distance1 = TMath::Min(delta1,delta2);
4897 rkink = TMath::Sqrt(radius[0]);
4900 rkink = TMath::Sqrt(radius[1]);
4902 if (distance1>cdist2) continue;
4905 AliTPCseed * track1 = (AliTPCseed*)array->At(j);
4908 Int_t row0 = GetRowNumber(rkink);
4909 if (row0<10) continue;
4910 if (row0>150) continue;
4913 Float_t dens00=-1,dens01=-1;
4914 Float_t dens10=-1,dens11=-1;
4916 Int_t found,foundable,ishared;
4917 track0->GetClusterStatistic(0,row0-5, found, foundable,ishared,kFALSE);
4918 if (foundable>5) dens00 = Float_t(found)/Float_t(foundable);
4919 track0->GetClusterStatistic(row0+5,155, found, foundable,ishared,kFALSE);
4920 if (foundable>5) dens01 = Float_t(found)/Float_t(foundable);
4922 track1->GetClusterStatistic(0,row0-5, found, foundable,ishared,kFALSE);
4923 if (foundable>10) dens10 = Float_t(found)/Float_t(foundable);
4924 track1->GetClusterStatistic(row0+5,155, found, foundable,ishared,kFALSE);
4925 if (foundable>10) dens11 = Float_t(found)/Float_t(foundable);
4927 if (dens00<dens10 && dens01<dens11) continue;
4928 if (dens00>dens10 && dens01>dens11) continue;
4929 if (TMath::Max(dens00,dens10)<0.1) continue;
4930 if (TMath::Max(dens01,dens11)<0.3) continue;
4932 if (TMath::Min(dens00,dens10)>0.6) continue;
4933 if (TMath::Min(dens01,dens11)>0.6) continue;
4936 AliTPCseed * ktrack0, *ktrack1;
4945 if (TMath::Abs(ktrack0->GetC())>5) continue; // cut on the curvature for mother particle
4946 AliExternalTrackParam paramm(*ktrack0);
4947 AliExternalTrackParam paramd(*ktrack1);
4948 if (row0>60&&ktrack1->GetReference().GetX()>90.)new (¶md) AliExternalTrackParam(ktrack1->GetReference());
4951 kink->SetMother(paramm);
4952 kink->SetDaughter(paramd);
4955 Float_t x[3] = { kink->GetPosition()[0],kink->GetPosition()[1],kink->GetPosition()[2]};
4957 fkParam->Transform0to1(x,index);
4958 fkParam->Transform1to2(x,index);
4959 row0 = GetRowNumber(x[0]);
4961 if (kink->GetR()<100) continue;
4962 if (kink->GetR()>240) continue;
4963 if (kink->GetPosition()[2]/kink->GetR()>AliTPCReconstructor::GetCtgRange()) continue; //out of fiducial volume
4964 if (kink->GetDistance()>cdist3) continue;
4965 Float_t dird = kink->GetDaughterP()[0]*kink->GetPosition()[0]+kink->GetDaughterP()[1]*kink->GetPosition()[1]; // rough direction estimate
4966 if (dird<0) continue;
4968 Float_t dirm = kink->GetMotherP()[0]*kink->GetPosition()[0]+kink->GetMotherP()[1]*kink->GetPosition()[1]; // rough direction estimate
4969 if (dirm<0) continue;
4970 Float_t mpt = TMath::Sqrt(kink->GetMotherP()[0]*kink->GetMotherP()[0]+kink->GetMotherP()[1]*kink->GetMotherP()[1]);
4971 if (mpt<0.2) continue;
4974 //for high momenta momentum not defined well in first iteration
4975 Double_t qt = TMath::Sin(kink->GetAngle(2))*ktrack1->GetP();
4976 if (qt>0.35) continue;
4979 kink->SetLabel(CookLabel(ktrack0,0.4,0,row0),0);
4980 kink->SetLabel(CookLabel(ktrack1,0.4,row0,160),1);
4982 kink->SetTPCDensity(dens00,0,0);
4983 kink->SetTPCDensity(dens01,0,1);
4984 kink->SetTPCDensity(dens10,1,0);
4985 kink->SetTPCDensity(dens11,1,1);
4986 kink->SetIndex(i,0);
4987 kink->SetIndex(j,1);
4990 kink->SetTPCDensity(dens10,0,0);
4991 kink->SetTPCDensity(dens11,0,1);
4992 kink->SetTPCDensity(dens00,1,0);
4993 kink->SetTPCDensity(dens01,1,1);
4994 kink->SetIndex(j,0);
4995 kink->SetIndex(i,1);
4998 if (mpt<1||kink->GetAngle(2)>0.1){
4999 // angle and densities not defined yet
5000 if (kink->GetTPCDensityFactor()<0.8) continue;
5001 if ((2-kink->GetTPCDensityFactor())*kink->GetDistance() >0.25) continue;
5002 if (kink->GetAngle(2)*ktrack0->GetP()<0.003) continue; //too small angle
5003 if (kink->GetAngle(2)>0.2&&kink->GetTPCDensityFactor()<1.15) continue;
5004 if (kink->GetAngle(2)>0.2&&kink->GetTPCDensity(0,1)>0.05) continue;
5006 Float_t criticalangle = track0->GetSigmaSnp2()+track0->GetSigmaTgl2();
5007 criticalangle+= track1->GetSigmaSnp2()+track1->GetSigmaTgl2();
5008 criticalangle= 3*TMath::Sqrt(criticalangle);
5009 if (criticalangle>0.02) criticalangle=0.02;
5010 if (kink->GetAngle(2)<criticalangle) continue;
5013 Int_t drow = Int_t(2.+0.5/(0.05+kink->GetAngle(2))); // overlap region defined
5014 Float_t shapesum =0;
5016 for ( Int_t row = row0-drow; row<row0+drow;row++){
5017 if (row<0) continue;
5018 if (row>155) continue;
5019 if (ktrack0->GetClusterPointer(row)){
5020 AliTPCTrackerPoint *point =ktrack0->GetTrackPoint(row);
5021 shapesum+=point->GetSigmaY()+point->GetSigmaZ();
5024 if (ktrack1->GetClusterPointer(row)){
5025 AliTPCTrackerPoint *point =ktrack1->GetTrackPoint(row);
5026 shapesum+=point->GetSigmaY()+point->GetSigmaZ();
5031 kink->SetShapeFactor(-1.);
5034 kink->SetShapeFactor(shapesum/sum);
5036 // esd->AddKink(kink);
5037 kinks->AddLast(kink);
5043 // sort the kinks according quality - and refit them towards vertex
5045 Int_t nkinks = kinks->GetEntriesFast();
5046 Float_t *quality = new Float_t[nkinks];
5047 Int_t *indexes = new Int_t[nkinks];
5048 AliTPCseed *mothers = new AliTPCseed[nkinks];
5049 AliTPCseed *daughters = new AliTPCseed[nkinks];
5052 for (Int_t i=0;i<nkinks;i++){
5054 AliKink *kinkl = (AliKink*)kinks->At(i);
5056 // refit kinks towards vertex
5058 Int_t index0 = kinkl->GetIndex(0);
5059 Int_t index1 = kinkl->GetIndex(1);
5060 AliTPCseed * ktrack0 = (AliTPCseed*)array->At(index0);
5061 AliTPCseed * ktrack1 = (AliTPCseed*)array->At(index1);
5063 Int_t sumn=ktrack0->GetNumberOfClusters()+ktrack1->GetNumberOfClusters();
5065 // Refit Kink under if too small angle
5067 if (kinkl->GetAngle(2)<0.05){
5068 kinkl->SetTPCRow0(GetRowNumber(kinkl->GetR()));
5069 Int_t row0 = kinkl->GetTPCRow0();
5070 Int_t drow = Int_t(2.+0.5/(0.05+kinkl->GetAngle(2)));
5073 Int_t last = row0-drow;
5074 if (last<40) last=40;
5075 if (last<ktrack0->GetFirstPoint()+25) last = ktrack0->GetFirstPoint()+25;
5076 AliTPCseed* seed0 = ReSeed(ktrack0,last,kFALSE);
5079 Int_t first = row0+drow;
5080 if (first>130) first=130;
5081 if (first>ktrack1->GetLastPoint()-25) first = TMath::Max(ktrack1->GetLastPoint()-25,30);
5082 AliTPCseed* seed1 = ReSeed(ktrack1,first,kTRUE);
5084 if (seed0 && seed1){
5085 kinkl->SetStatus(1,8);
5086 if (RefitKink(*seed0,*seed1,*kinkl)) kinkl->SetStatus(1,9);
5087 row0 = GetRowNumber(kinkl->GetR());
5088 sumn = seed0->GetNumberOfClusters()+seed1->GetNumberOfClusters();
5089 mothers[i] = *seed0;
5090 daughters[i] = *seed1;
5093 delete kinks->RemoveAt(i);
5094 if (seed0) delete seed0;
5095 if (seed1) delete seed1;
5098 if (kinkl->GetDistance()>0.5 || kinkl->GetR()<110 || kinkl->GetR()>240) {
5099 delete kinks->RemoveAt(i);
5100 if (seed0) delete seed0;
5101 if (seed1) delete seed1;
5109 if (kinkl) quality[i] = 160*((0.1+kinkl->GetDistance())*(2.-kinkl->GetTPCDensityFactor()))/(sumn+40.); //the longest -clossest will win
5111 TMath::Sort(nkinks,quality,indexes,kFALSE);
5113 //remove double find kinks
5115 for (Int_t ikink0=1;ikink0<nkinks;ikink0++){
5116 AliKink * kink0 = (AliKink*) kinks->At(indexes[ikink0]);
5117 if (!kink0) continue;
5119 for (Int_t ikink1=0;ikink1<ikink0;ikink1++){
5120 if (!kink0) continue;
5121 AliKink * kink1 = (AliKink*) kinks->At(indexes[ikink1]);
5122 if (!kink1) continue;
5123 // if not close kink continue
5124 if (TMath::Abs(kink1->GetPosition()[2]-kink0->GetPosition()[2])>10) continue;
5125 if (TMath::Abs(kink1->GetPosition()[1]-kink0->GetPosition()[1])>10) continue;
5126 if (TMath::Abs(kink1->GetPosition()[0]-kink0->GetPosition()[0])>10) continue;
5128 AliTPCseed &mother0 = mothers[indexes[ikink0]];
5129 AliTPCseed &daughter0 = daughters[indexes[ikink0]];
5130 AliTPCseed &mother1 = mothers[indexes[ikink1]];
5131 AliTPCseed &daughter1 = daughters[indexes[ikink1]];
5132 Int_t row0 = (kink0->GetTPCRow0()+kink1->GetTPCRow0())/2;
5141 for (Int_t i=0;i<row0;i++){
5142 if (mother0.GetClusterIndex(i)>0 && mother1.GetClusterIndex(i)>0){
5145 if (mother0.GetClusterIndex(i)==mother1.GetClusterIndex(i)){
5152 for (Int_t i=row0;i<158;i++){
5153 if (daughter0.GetClusterIndex(i)>0 && daughter0.GetClusterIndex(i)>0){
5156 if (mother0.GetClusterIndex(i)==mother1.GetClusterIndex(i)){
5162 Float_t ratio = Float_t(same+1)/Float_t(both+1);
5163 Float_t ratiom = Float_t(samem+1)/Float_t(bothm+1);
5164 Float_t ratiod = Float_t(samed+1)/Float_t(bothd+1);
5165 if (ratio>0.3 && ratiom>0.5 &&ratiod>0.5) {
5166 Int_t sum0 = mother0.GetNumberOfClusters()+daughter0.GetNumberOfClusters();
5167 Int_t sum1 = mother1.GetNumberOfClusters()+daughter1.GetNumberOfClusters();
5169 shared[kink0->GetIndex(0)]= kTRUE;
5170 shared[kink0->GetIndex(1)]= kTRUE;
5171 delete kinks->RemoveAt(indexes[ikink0]);
5174 shared[kink1->GetIndex(0)]= kTRUE;
5175 shared[kink1->GetIndex(1)]= kTRUE;
5176 delete kinks->RemoveAt(indexes[ikink1]);
5183 for (Int_t i=0;i<nkinks;i++){
5184 AliKink * kinkl = (AliKink*) kinks->At(indexes[i]);
5185 if (!kinkl) continue;
5186 kinkl->SetTPCRow0(GetRowNumber(kinkl->GetR()));
5187 Int_t index0 = kinkl->GetIndex(0);
5188 Int_t index1 = kinkl->GetIndex(1);
5189 if (circular[index0]||(circular[index1]&&kinkl->GetDistance()>0.2)) continue;
5190 kinkl->SetMultiple(usage[index0],0);
5191 kinkl->SetMultiple(usage[index1],1);
5192 if (kinkl->GetMultiple()[0]+kinkl->GetMultiple()[1]>2) continue;
5193 if (kinkl->GetMultiple()[0]+kinkl->GetMultiple()[1]>0 && quality[indexes[i]]>0.2) continue;
5194 if (kinkl->GetMultiple()[0]+kinkl->GetMultiple()[1]>0 && kinkl->GetDistance()>0.2) continue;
5195 if (circular[index0]||(circular[index1]&&kinkl->GetDistance()>0.1)) continue;
5197 AliTPCseed * ktrack0 = (AliTPCseed*)array->At(index0);
5198 AliTPCseed * ktrack1 = (AliTPCseed*)array->At(index1);
5199 if (!ktrack0 || !ktrack1) continue;
5200 Int_t index = esd->AddKink(kinkl);
5203 if ( ktrack0->GetKinkIndex(0)==0 && ktrack1->GetKinkIndex(0)==0) { //best kink
5204 if (mothers[indexes[i]].GetNumberOfClusters()>20 && daughters[indexes[i]].GetNumberOfClusters()>20 && (mothers[indexes[i]].GetNumberOfClusters()+daughters[indexes[i]].GetNumberOfClusters())>100){
5205 *ktrack0 = mothers[indexes[i]];
5206 *ktrack1 = daughters[indexes[i]];
5210 ktrack0->SetKinkIndex(usage[index0],-(index+1));
5211 ktrack1->SetKinkIndex(usage[index1], (index+1));
5216 // Remove tracks corresponding to shared kink's
5218 for (Int_t i=0;i<nentries;i++){
5219 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5220 if (!track0) continue;
5221 if (track0->GetKinkIndex(0)!=0) continue;
5222 if (shared[i]) delete array->RemoveAt(i);
5227 RemoveUsed2(array,0.5,0.4,30);
5229 for (Int_t i=0;i<nentries;i++){
5230 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5231 if (!track0) continue;
5232 track0->CookdEdx(0.02,0.6);
5236 for (Int_t i=0;i<nentries;i++){
5237 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5238 if (!track0) continue;
5239 if (track0->Pt()<1.4) continue;
5240 //remove double high momenta tracks - overlapped with kink candidates
5243 for (Int_t icl=track0->GetFirstPoint();icl<track0->GetLastPoint(); icl++){
5244 if (track0->GetClusterPointer(icl)!=0){
5246 if (track0->GetClusterPointer(icl)->IsUsed(10)) ishared++;
5249 if (Float_t(ishared+1)/Float_t(all+1)>0.5) {
5250 delete array->RemoveAt(i);
5254 if (track0->GetKinkIndex(0)!=0) continue;
5255 if (track0->GetNumberOfClusters()<80) continue;
5257 AliTPCseed *pmother = new AliTPCseed();
5258 AliTPCseed *pdaughter = new AliTPCseed();
5259 AliKink *pkink = new AliKink;
5261 AliTPCseed & mother = *pmother;
5262 AliTPCseed & daughter = *pdaughter;
5263 AliKink & kinkl = *pkink;
5264 if (CheckKinkPoint(track0,mother,daughter, kinkl)){
5265 if (mother.GetNumberOfClusters()<30||daughter.GetNumberOfClusters()<20) {
5269 continue; //too short tracks
5271 if (mother.Pt()<1.4) {
5277 Int_t row0= kinkl.GetTPCRow0();
5278 if (kinkl.GetDistance()>0.5 || kinkl.GetR()<110. || kinkl.GetR()>240.) {
5285 Int_t index = esd->AddKink(&kinkl);
5286 mother.SetKinkIndex(0,-(index+1));
5287 daughter.SetKinkIndex(0,index+1);
5288 if (mother.GetNumberOfClusters()>50) {
5289 delete array->RemoveAt(i);
5290 array->AddAt(new AliTPCseed(mother),i);
5293 array->AddLast(new AliTPCseed(mother));
5295 array->AddLast(new AliTPCseed(daughter));
5296 for (Int_t icl=0;icl<row0;icl++) {
5297 if (mother.GetClusterPointer(icl)) mother.GetClusterPointer(icl)->Use(20);
5300 for (Int_t icl=row0;icl<158;icl++) {
5301 if (daughter.GetClusterPointer(icl)) daughter.GetClusterPointer(icl)->Use(20);
5310 delete [] daughters;
5332 printf("Ncandidates=\t%d\t%d\t%d\t%d\n",esd->GetNumberOfKinks(),ncandidates,ntracks,nall);
5336 void AliTPCtrackerMI::FindV0s(const TObjArray * array, AliESDEvent *esd)
5342 TObjArray *tpcv0s = new TObjArray(100000);
5343 Int_t nentries = array->GetEntriesFast();
5344 AliHelix *helixes = new AliHelix[nentries];
5345 Int_t *sign = new Int_t[nentries];
5346 Float_t *alpha = new Float_t[nentries];
5347 Float_t *z0 = new Float_t[nentries];
5348 Float_t *dca = new Float_t[nentries];
5349 Float_t *sdcar = new Float_t[nentries];
5350 Float_t *cdcar = new Float_t[nentries];
5351 Float_t *pulldcar = new Float_t[nentries];
5352 Float_t *pulldcaz = new Float_t[nentries];
5353 Float_t *pulldca = new Float_t[nentries];
5354 Bool_t *isPrim = new Bool_t[nentries];
5355 const AliESDVertex * primvertex = esd->GetVertex();
5356 Double_t zvertex = primvertex->GetZv();
5358 // nentries = array->GetEntriesFast();
5360 for (Int_t i=0;i<nentries;i++){
5363 AliTPCseed* track = (AliTPCseed*)array->At(i);
5364 if (!track) continue;
5365 track->GetV0Indexes()[0] = 0; //rest v0 indexes
5366 track->GetV0Indexes()[1] = 0; //rest v0 indexes
5367 track->GetV0Indexes()[2] = 0; //rest v0 indexes
5369 alpha[i] = track->GetAlpha();
5370 new (&helixes[i]) AliHelix(*track);
5372 helixes[i].Evaluate(0,xyz);
5373 sign[i] = (track->GetC()>0) ? -1:1;
5377 if (track->GetProlongation(0,y,z)) z0[i] = z;
5378 dca[i] = track->GetD(0,0);
5380 // dca error parrameterezation + pulls
5382 sdcar[i] = TMath::Sqrt(0.150*0.150+(100*track->GetC())*(100*track->GetC()));
5383 if (TMath::Abs(track->GetTgl())>1) sdcar[i]*=2.5;
5384 cdcar[i] = TMath::Exp((TMath::Abs(track->GetC())-0.0106)*525.3);
5385 pulldcar[i] = (dca[i]-cdcar[i])/sdcar[i];
5386 pulldcaz[i] = (z0[i]-zvertex)/sdcar[i];
5387 pulldca[i] = TMath::Sqrt(pulldcar[i]*pulldcar[i]+pulldcaz[i]*pulldcaz[i]);
5388 if (track->TPCrPID(1)+track->TPCrPID(2)+track->TPCrPID(3)>0.5) {
5389 if (pulldca[i]<3.) isPrim[i]=kTRUE; //pion, muon and Kaon 3 sigma cut
5391 if (track->TPCrPID(4)>0.5) {
5392 if (pulldca[i]<0.5) isPrim[i]=kTRUE; //proton 0.5 sigma cut
5394 if (track->TPCrPID(0)>0.4) {
5395 isPrim[i]=kFALSE; //electron no sigma cut
5402 Int_t ncandidates =0;
5405 Double_t phase[2][2],radius[2];
5411 Float_t fprimvertex[3]={GetX(),GetY(),GetZ()};
5413 Double_t cradius0 = 10*10;
5414 Double_t cradius1 = 200*200;
5417 Double_t cpointAngle = 0.95;
5419 Double_t delta[2]={10000,10000};
5420 for (Int_t i =0;i<nentries;i++){
5421 if (sign[i]==0) continue;
5422 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5423 if (!track0) continue;
5424 if (AliTPCReconstructor::StreamLevel()>1){
5425 TTreeSRedirector &cstream = *fDebugStreamer;
5430 "zvertex="<<zvertex<<
5431 "sdcar0="<<sdcar[i]<<
5432 "cdcar0="<<cdcar[i]<<
5433 "pulldcar0="<<pulldcar[i]<<
5434 "pulldcaz0="<<pulldcaz[i]<<
5435 "pulldca0="<<pulldca[i]<<
5436 "isPrim="<<isPrim[i]<<
5440 if (track0->GetSigned1Pt()<0) continue;
5441 if (track0->GetKinkIndex(0)>0||isPrim[i]) continue; //daughter kink
5443 if (TMath::Abs(helixes[i].GetHelix(4))<0.000000001) continue;
5448 for (Int_t j =0;j<nentries;j++){
5449 AliTPCseed * track1 = (AliTPCseed*)array->At(j);
5450 if (!track1) continue;
5451 if (track1->GetKinkIndex(0)>0 || isPrim[j]) continue; //daughter kink
5452 if (sign[j]*sign[i]>0) continue;
5453 if (TMath::Abs(helixes[j].GetHelix(4))<0.000001) continue;
5454 if (track0->GetCircular()+track1->GetCircular()>1) continue; //circling -returning track
5457 // DCA to prim vertex cut
5463 Int_t npoints = helixes[i].GetRPHIintersections(helixes[j], phase, radius,cdist2);
5464 if (npoints<1) continue;
5468 if (radius[0]<cradius0||radius[0]>cradius1) continue;
5469 helixes[i].LinearDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta[0]);
5470 if (delta[0]>cdist1) continue;
5473 if (TMath::Max(radius[0],radius[1])<cradius0|| TMath::Min(radius[0],radius[1])>cradius1) continue;
5474 helixes[i].LinearDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta[0]);
5475 helixes[i].LinearDCA(helixes[j],phase[1][0],phase[1][1],radius[1],delta[1]);
5476 if (delta[1]<delta[0]) iclosest=1;
5477 if (delta[iclosest]>cdist1) continue;
5479 helixes[i].ParabolicDCA(helixes[j],phase[iclosest][0],phase[iclosest][1],radius[iclosest],delta[iclosest]);
5480 if (radius[iclosest]<cradius0 || radius[iclosest]>cradius1 || delta[iclosest]>cdist1) continue;
5482 Double_t pointAngle = helixes[i].GetPointAngle(helixes[j],phase[iclosest],fprimvertex);
5483 if (pointAngle<cpointAngle) continue;
5485 Bool_t isGamma = kFALSE;
5486 vertex.SetParamP(*track0); //track0 - plus
5487 vertex.SetParamN(*track1); //track1 - minus
5488 vertex.Update(fprimvertex);
5489 if (track0->TPCrPID(0)>0.3&&track1->TPCrPID(0)>0.3&&vertex.GetAnglep()[2]<0.15) isGamma=kTRUE; // gamma conversion candidate
5490 Double_t pointAngle2 = vertex.GetV0CosineOfPointingAngle();
5492 if (vertex.GetV0CosineOfPointingAngle()<cpointAngle && (!isGamma)) continue;// point angle cut
5493 //Bo: if (vertex.GetDist2()>2&&(!isGamma)) continue; // point angle cut
5494 if (vertex.GetDcaV0Daughters()>2&&(!isGamma)) continue;//Bo: // point angle cut
5495 Float_t sigmae = 0.15*0.15;
5496 if (vertex.GetRr()<80)
5497 sigmae += (sdcar[i]*sdcar[i]+sdcar[j]*sdcar[j])*(1.-vertex.GetRr()/80.)*(1.-vertex.GetRr()/80.);
5498 sigmae+= TMath::Sqrt(sigmae);
5499 //Bo: if (vertex.GetDist2()/sigmae>3.&&(!isGamma)) continue;
5500 if (vertex.GetDcaV0Daughters()/sigmae>3.&&(!isGamma)) continue;
5501 Float_t densb0=0,densb1=0,densa0=0,densa1=0;
5502 Int_t row0 = GetRowNumber(vertex.GetRr());
5504 //Bo: if (vertex.GetDist2()>0.2) continue;
5505 if (vertex.GetDcaV0Daughters()>0.2) continue;
5506 densb0 = track0->Density2(0,row0-5);
5507 densb1 = track1->Density2(0,row0-5);
5508 if (densb0>0.3|| densb1>0.3) continue; //clusters before vertex
5509 densa0 = track0->Density2(row0+5,row0+40);
5510 densa1 = track1->Density2(row0+5,row0+40);
5511 if ((densa0<0.4|| densa1<0.4)&&(!isGamma)) continue; //missing clusters after vertex
5514 densa0 = track0->Density2(0,40); //cluster density
5515 densa1 = track1->Density2(0,40); //cluster density
5516 if ((vertex.GetRr()<80&&densa0+densa1<1.)&&(!isGamma)) continue;
5518 //Bo: vertex.SetLab(0,track0->GetLabel());
5519 //Bo: vertex.SetLab(1,track1->GetLabel());
5520 vertex.SetChi2After((densa0+densa1)*0.5);
5521 vertex.SetChi2Before((densb0+densb1)*0.5);
5522 vertex.SetIndex(0,i);
5523 vertex.SetIndex(1,j);
5524 //Bo: vertex.SetStatus(1); // TPC v0 candidate
5525 vertex.SetOnFlyStatus(2);//Bo: // TPC v0 candidate
5526 //Bo: vertex.SetRp(track0->TPCrPIDs());
5527 //Bo: vertex.SetRm(track1->TPCrPIDs());
5528 tpcv0s->AddLast(new AliESDv0(vertex));
5531 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
5532 Double_t radiusm= (delta[0]<delta[1])? TMath::Sqrt(radius[0]):TMath::Sqrt(radius[1]);
5533 Double_t deltam= (delta[0]<delta[1])? TMath::Sqrt(delta[0]):TMath::Sqrt(delta[1]);
5534 if (AliTPCReconstructor::StreamLevel()>1) {
5535 Int_t lab0=track0->GetLabel();
5536 Int_t lab1=track1->GetLabel();
5537 Char_t c0=track0->GetCircular();
5538 Char_t c1=track1->GetCircular();
5539 TTreeSRedirector &cstream = *fDebugStreamer;
5542 "vertex.="<<&vertex<<
5545 "Helix0.="<<&helixes[i]<<
5548 "Helix1.="<<&helixes[j]<<
5549 "pointAngle="<<pointAngle<<
5550 "pointAngle2="<<pointAngle2<<
5555 "zvertex="<<zvertex<<
5558 "npoints="<<npoints<<
5559 "radius0="<<radius[0]<<
5560 "delta0="<<delta[0]<<
5561 "radius1="<<radius[1]<<
5562 "delta1="<<delta[1]<<
5563 "radiusm="<<radiusm<<
5565 "sdcar0="<<sdcar[i]<<
5566 "sdcar1="<<sdcar[j]<<
5567 "cdcar0="<<cdcar[i]<<
5568 "cdcar1="<<cdcar[j]<<
5569 "pulldcar0="<<pulldcar[i]<<
5570 "pulldcar1="<<pulldcar[j]<<
5571 "pulldcaz0="<<pulldcaz[i]<<
5572 "pulldcaz1="<<pulldcaz[j]<<
5573 "pulldca0="<<pulldca[i]<<
5574 "pulldca1="<<pulldca[j]<<
5584 Float_t *quality = new Float_t[ncandidates];
5585 Int_t *indexes = new Int_t[ncandidates];
5587 for (Int_t i=0;i<ncandidates;i++){
5589 AliESDv0 *v0 = (AliESDv0*)tpcv0s->At(i);
5590 quality[i] = 1./(1.00001-v0->GetV0CosineOfPointingAngle()); //base point angle
5591 // quality[i] /= (0.5+v0->GetDist2());
5592 // quality[i] *= v0->GetChi2After(); //density factor
5594 Int_t index0 = v0->GetIndex(0);
5595 Int_t index1 = v0->GetIndex(1);
5596 //Bo: Double_t minpulldca = TMath::Min(2.+pulldca[v0->GetIndex(0)],(2.+pulldca[v0->GetIndex(1)]) ); //pull
5597 Double_t minpulldca = TMath::Min(2.+pulldca[index0],(2.+pulldca[index1]) );//Bo:
5601 AliTPCseed * track0 = (AliTPCseed*)array->At(index0);
5602 AliTPCseed * track1 = (AliTPCseed*)array->At(index1);
5603 if (track0->TPCrPID(0)>0.3&&track1->TPCrPID(0)>0.3&&v0->GetAnglep()[2]<0.15) quality[i]+=1000000; // gamma conversion candidate
5604 if (track0->TPCrPID(4)>0.9||(track1->TPCrPID(4)>0.9&&minpulldca>4)) quality[i]*=10; // lambda candidate candidate
5607 TMath::Sort(ncandidates,quality,indexes,kTRUE);
5610 for (Int_t i=0;i<ncandidates;i++){
5611 AliESDv0 * v0 = (AliESDv0*)tpcv0s->At(indexes[i]);
5613 Int_t index0 = v0->GetIndex(0);
5614 Int_t index1 = v0->GetIndex(1);
5615 AliTPCseed * track0 = (AliTPCseed*)array->At(index0);
5616 AliTPCseed * track1 = (AliTPCseed*)array->At(index1);
5617 if (!track0||!track1) {
5621 Bool_t accept =kTRUE; //default accept
5622 Int_t *v0indexes0 = track0->GetV0Indexes();
5623 Int_t *v0indexes1 = track1->GetV0Indexes();
5625 Int_t order0 = (v0indexes0[0]!=0) ? 1:0;
5626 Int_t order1 = (v0indexes1[0]!=0) ? 1:0;
5627 if (v0indexes0[1]!=0) order0 =2;
5628 if (v0indexes1[1]!=0) order1 =2;
5630 if (v0indexes0[2]!=0) {order0=3; accept=kFALSE;}
5631 if (v0indexes0[2]!=0) {order1=3; accept=kFALSE;}
5633 AliESDv0 * v02 = v0;
5635 //Bo: v0->SetOrder(0,order0);
5636 //Bo: v0->SetOrder(1,order1);
5637 //Bo: v0->SetOrder(1,order0+order1);
5638 v0->SetOnFlyStatus(kTRUE);
5639 Int_t index = esd->AddV0(v0);
5640 v02 = esd->GetV0(index);
5641 v0indexes0[order0]=index;
5642 v0indexes1[order1]=index;
5646 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
5647 if (AliTPCReconstructor::StreamLevel()>1) {
5648 Int_t lab0=track0->GetLabel();
5649 Int_t lab1=track1->GetLabel();
5650 TTreeSRedirector &cstream = *fDebugStreamer;
5659 "dca0="<<dca[index0]<<
5660 "dca1="<<dca[index1]<<
5664 "quality="<<quality[i]<<
5665 "pulldca0="<<pulldca[index0]<<
5666 "pulldca1="<<pulldca[index1]<<
5690 printf("TPC V0 finder : naccepted\t%d\tncandidates\t%d\tntracks\t%d\tnall\t%d\n",naccepted,ncandidates,ntracks,nall);
5694 Int_t AliTPCtrackerMI::RefitKink(AliTPCseed &mother, AliTPCseed &daughter, AliESDkink &knk)
5697 // refit kink towards to the vertex
5700 AliKink &kink=(AliKink &)knk;
5702 Int_t row0 = GetRowNumber(kink.GetR());
5703 FollowProlongation(mother,0);
5704 mother.Reset(kFALSE);
5706 FollowProlongation(daughter,row0);
5707 daughter.Reset(kFALSE);
5708 FollowBackProlongation(daughter,158);
5709 daughter.Reset(kFALSE);
5710 Int_t first = TMath::Max(row0-20,30);
5711 Int_t last = TMath::Min(row0+20,140);
5713 const Int_t kNdiv =5;
5714 AliTPCseed param0[kNdiv]; // parameters along the track
5715 AliTPCseed param1[kNdiv]; // parameters along the track
5716 AliKink kinks[kNdiv]; // corresponding kink parameters
5719 for (Int_t irow=0; irow<kNdiv;irow++){
5720 rows[irow] = first +((last-first)*irow)/(kNdiv-1);
5722 // store parameters along the track
5724 for (Int_t irow=0;irow<kNdiv;irow++){
5725 FollowBackProlongation(mother, rows[irow]);
5726 FollowProlongation(daughter,rows[kNdiv-1-irow]);
5727 param0[irow] = mother;
5728 param1[kNdiv-1-irow] = daughter;
5732 for (Int_t irow=0; irow<kNdiv-1;irow++){
5733 if (param0[irow].GetNumberOfClusters()<kNdiv||param1[irow].GetNumberOfClusters()<kNdiv) continue;
5734 kinks[irow].SetMother(param0[irow]);
5735 kinks[irow].SetDaughter(param1[irow]);
5736 kinks[irow].Update();
5739 // choose kink with best "quality"
5741 Double_t mindist = 10000;
5742 for (Int_t irow=0;irow<kNdiv;irow++){
5743 if (param0[irow].GetNumberOfClusters()<20||param1[irow].GetNumberOfClusters()<20) continue;
5744 if (TMath::Abs(kinks[irow].GetR())>240.) continue;
5745 if (TMath::Abs(kinks[irow].GetR())<100.) continue;
5747 Float_t normdist = TMath::Abs(param0[irow].GetX()-kinks[irow].GetR())*(0.1+kink.GetDistance());
5748 normdist/= (param0[irow].GetNumberOfClusters()+param1[irow].GetNumberOfClusters()+40.);
5749 if (normdist < mindist){
5755 if (index==-1) return 0;
5758 param0[index].Reset(kTRUE);
5759 FollowProlongation(param0[index],0);
5761 mother = param0[index];
5762 daughter = param1[index]; // daughter in vertex
5764 kink.SetMother(mother);
5765 kink.SetDaughter(daughter);
5767 kink.SetTPCRow0(GetRowNumber(kink.GetR()));
5768 kink.SetTPCncls(param0[index].GetNumberOfClusters(),0);
5769 kink.SetTPCncls(param1[index].GetNumberOfClusters(),1);
5770 kink.SetLabel(CookLabel(&mother,0.4, 0,kink.GetTPCRow0()),0);
5771 kink.SetLabel(CookLabel(&daughter,0.4, kink.GetTPCRow0(),160),1);
5772 mother.SetLabel(kink.GetLabel(0));
5773 daughter.SetLabel(kink.GetLabel(1));
5779 void AliTPCtrackerMI::UpdateKinkQualityM(AliTPCseed * seed){
5781 // update Kink quality information for mother after back propagation
5783 if (seed->GetKinkIndex(0)>=0) return;
5784 for (Int_t ikink=0;ikink<3;ikink++){
5785 Int_t index = seed->GetKinkIndex(ikink);
5786 if (index>=0) break;
5787 index = TMath::Abs(index)-1;
5788 AliESDkink * kink = fEvent->GetKink(index);
5789 kink->SetTPCDensity(-1,0,0);
5790 kink->SetTPCDensity(1,0,1);
5792 Int_t row0 = kink->GetTPCRow0() - 2 - Int_t( 0.5/ (0.05+kink->GetAngle(2)));
5793 if (row0<15) row0=15;
5795 Int_t row1 = kink->GetTPCRow0() + 2 + Int_t( 0.5/ (0.05+kink->GetAngle(2)));
5796 if (row1>145) row1=145;
5798 Int_t found,foundable,shared;
5799 seed->GetClusterStatistic(0,row0, found, foundable,shared,kFALSE);
5800 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),0,0);
5801 seed->GetClusterStatistic(row1,155, found, foundable,shared,kFALSE);
5802 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),0,1);
5807 void AliTPCtrackerMI::UpdateKinkQualityD(AliTPCseed * seed){
5809 // update Kink quality information for daughter after refit
5811 if (seed->GetKinkIndex(0)<=0) return;
5812 for (Int_t ikink=0;ikink<3;ikink++){
5813 Int_t index = seed->GetKinkIndex(ikink);
5814 if (index<=0) break;
5815 index = TMath::Abs(index)-1;
5816 AliESDkink * kink = fEvent->GetKink(index);
5817 kink->SetTPCDensity(-1,1,0);
5818 kink->SetTPCDensity(-1,1,1);
5820 Int_t row0 = kink->GetTPCRow0() -2 - Int_t( 0.5/ (0.05+kink->GetAngle(2)));
5821 if (row0<15) row0=15;
5823 Int_t row1 = kink->GetTPCRow0() +2 + Int_t( 0.5/ (0.05+kink->GetAngle(2)));
5824 if (row1>145) row1=145;
5826 Int_t found,foundable,shared;
5827 seed->GetClusterStatistic(0,row0, found, foundable,shared,kFALSE);
5828 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),1,0);
5829 seed->GetClusterStatistic(row1,155, found, foundable,shared,kFALSE);
5830 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),1,1);
5836 Int_t AliTPCtrackerMI::CheckKinkPoint(AliTPCseed*seed,AliTPCseed &mother, AliTPCseed &daughter, AliESDkink &knk)
5839 // check kink point for given track
5840 // if return value=0 kink point not found
5841 // otherwise seed0 correspond to mother particle
5842 // seed1 correspond to daughter particle
5843 // kink parameter of kink point
5844 AliKink &kink=(AliKink &)knk;
5846 Int_t middlerow = (seed->GetFirstPoint()+seed->GetLastPoint())/2;
5847 Int_t first = seed->GetFirstPoint();
5848 Int_t last = seed->GetLastPoint();
5849 if (last-first<20) return 0; // shortest length - 2*30 = 60 pad-rows
5852 AliTPCseed *seed1 = ReSeed(seed,middlerow+20, kTRUE); //middle of chamber
5853 if (!seed1) return 0;
5854 FollowProlongation(*seed1,seed->GetLastPoint()-20);
5855 seed1->Reset(kTRUE);
5856 FollowProlongation(*seed1,158);
5857 seed1->Reset(kTRUE);
5858 last = seed1->GetLastPoint();
5860 AliTPCseed *seed0 = new AliTPCseed(*seed);
5861 seed0->Reset(kFALSE);
5864 AliTPCseed param0[20]; // parameters along the track
5865 AliTPCseed param1[20]; // parameters along the track
5866 AliKink kinks[20]; // corresponding kink parameters
5868 for (Int_t irow=0; irow<20;irow++){
5869 rows[irow] = first +((last-first)*irow)/19;
5871 // store parameters along the track
5873 for (Int_t irow=0;irow<20;irow++){
5874 FollowBackProlongation(*seed0, rows[irow]);
5875 FollowProlongation(*seed1,rows[19-irow]);
5876 param0[irow] = *seed0;
5877 param1[19-irow] = *seed1;
5881 for (Int_t irow=0; irow<19;irow++){
5882 kinks[irow].SetMother(param0[irow]);
5883 kinks[irow].SetDaughter(param1[irow]);
5884 kinks[irow].Update();
5887 // choose kink with biggest change of angle
5889 Double_t maxchange= 0;
5890 for (Int_t irow=1;irow<19;irow++){
5891 if (TMath::Abs(kinks[irow].GetR())>240.) continue;
5892 if (TMath::Abs(kinks[irow].GetR())<110.) continue;
5893 Float_t quality = TMath::Abs(kinks[irow].GetAngle(2))/(3.+TMath::Abs(kinks[irow].GetR()-param0[irow].GetX()));
5894 if ( quality > maxchange){
5895 maxchange = quality;
5902 if (index<0) return 0;
5904 Int_t row0 = GetRowNumber(kinks[index].GetR()); //row 0 estimate
5905 seed0 = new AliTPCseed(param0[index]);
5906 seed1 = new AliTPCseed(param1[index]);
5907 seed0->Reset(kFALSE);
5908 seed1->Reset(kFALSE);
5909 seed0->ResetCovariance(10.);
5910 seed1->ResetCovariance(10.);
5911 FollowProlongation(*seed0,0);
5912 FollowBackProlongation(*seed1,158);
5913 mother = *seed0; // backup mother at position 0
5914 seed0->Reset(kFALSE);
5915 seed1->Reset(kFALSE);
5916 seed0->ResetCovariance(10.);
5917 seed1->ResetCovariance(10.);
5919 first = TMath::Max(row0-20,0);
5920 last = TMath::Min(row0+20,158);
5922 for (Int_t irow=0; irow<20;irow++){
5923 rows[irow] = first +((last-first)*irow)/19;
5925 // store parameters along the track
5927 for (Int_t irow=0;irow<20;irow++){
5928 FollowBackProlongation(*seed0, rows[irow]);
5929 FollowProlongation(*seed1,rows[19-irow]);
5930 param0[irow] = *seed0;
5931 param1[19-irow] = *seed1;
5935 for (Int_t irow=0; irow<19;irow++){
5936 kinks[irow].SetMother(param0[irow]);
5937 kinks[irow].SetDaughter(param1[irow]);
5938 // param0[irow].Dump();
5939 //param1[irow].Dump();
5940 kinks[irow].Update();
5943 // choose kink with biggest change of angle
5946 for (Int_t irow=0;irow<20;irow++){
5947 if (TMath::Abs(kinks[irow].GetR())>250.) continue;
5948 if (TMath::Abs(kinks[irow].GetR())<90.) continue;
5949 Float_t quality = TMath::Abs(kinks[irow].GetAngle(2))/(3.+TMath::Abs(kinks[irow].GetR()-param0[irow].GetX()));
5950 if ( quality > maxchange){
5951 maxchange = quality;
5958 if (index==-1 || param0[index].GetNumberOfClusters()+param1[index].GetNumberOfClusters()<100){
5963 // Float_t anglesigma = TMath::Sqrt(param0[index].fC22+param0[index].fC33+param1[index].fC22+param1[index].fC33);
5965 kink.SetMother(param0[index]);
5966 kink.SetDaughter(param1[index]);
5968 row0 = GetRowNumber(kink.GetR());
5969 kink.SetTPCRow0(row0);
5970 kink.SetLabel(CookLabel(seed0,0.5,0,row0),0);
5971 kink.SetLabel(CookLabel(seed1,0.5,row0,158),1);
5972 kink.SetIndex(-10,0);
5973 kink.SetIndex(int(param0[index].GetNumberOfClusters()+param1[index].GetNumberOfClusters()),1);
5974 kink.SetTPCncls(param0[index].GetNumberOfClusters(),0);
5975 kink.SetTPCncls(param1[index].GetNumberOfClusters(),1);
5978 // new (&mother) AliTPCseed(param0[index]);
5979 daughter = param1[index];
5980 daughter.SetLabel(kink.GetLabel(1));
5981 param0[index].Reset(kTRUE);
5982 FollowProlongation(param0[index],0);
5983 mother = param0[index];
5984 mother.SetLabel(kink.GetLabel(0));
5994 AliTPCseed* AliTPCtrackerMI::ReSeed(AliTPCseed *t)
5997 // reseed - refit - track
6000 // Int_t last = fSectors->GetNRows()-1;
6002 if (fSectors == fOuterSec){
6003 first = TMath::Max(first, t->GetFirstPoint()-fInnerSec->GetNRows());
6007 first = t->GetFirstPoint();
6009 AliTPCseed * seed = MakeSeed(t,0.1,0.5,0.9);
6010 FollowBackProlongation(*t,fSectors->GetNRows()-1);
6012 FollowProlongation(*t,first);
6022 //_____________________________________________________________________________
6023 Int_t AliTPCtrackerMI::ReadSeeds(const TFile *inp) {
6024 //-----------------------------------------------------------------
6025 // This function reades track seeds.
6026 //-----------------------------------------------------------------
6027 TDirectory *savedir=gDirectory;
6029 TFile *in=(TFile*)inp;
6030 if (!in->IsOpen()) {
6031 cerr<<"AliTPCtrackerMI::ReadSeeds(): input file is not open !\n";
6036 TTree *seedTree=(TTree*)in->Get("Seeds");
6038 cerr<<"AliTPCtrackerMI::ReadSeeds(): ";
6039 cerr<<"can't get a tree with track seeds !\n";
6042 AliTPCtrack *seed=new AliTPCtrack;
6043 seedTree->SetBranchAddress("tracks",&seed);
6045 if (fSeeds==0) fSeeds=new TObjArray(15000);
6047 Int_t n=(Int_t)seedTree->GetEntries();
6048 for (Int_t i=0; i<n; i++) {
6049 seedTree->GetEvent(i);
6050 fSeeds->AddLast(new AliTPCseed(*seed/*,seed->GetAlpha()*/));
6059 Int_t AliTPCtrackerMI::Clusters2Tracks (AliESDEvent *esd)
6062 if (fSeeds) DeleteSeeds();
6065 if (!fSeeds) return 1;
6072 //_____________________________________________________________________________
6073 Int_t AliTPCtrackerMI::Clusters2Tracks() {
6074 //-----------------------------------------------------------------
6075 // This is a track finder.
6076 //-----------------------------------------------------------------
6077 TDirectory *savedir=gDirectory;
6081 fSeeds = Tracking();
6084 Info("Clusters2Tracks","Time for tracking: \t");timer.Print();timer.Start();
6086 //activate again some tracks
6087 for (Int_t i=0; i<fSeeds->GetEntriesFast(); i++) {
6088 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
6090 Int_t nc=t.GetNumberOfClusters();
6092 delete fSeeds->RemoveAt(i);
6096 if (pt->GetRemoval()==10) {
6097 if (pt->GetDensityFirst(20)>0.8 || pt->GetDensityFirst(30)>0.8 || pt->GetDensityFirst(40)>0.7)
6098 pt->Desactivate(10); // make track again active
6100 pt->Desactivate(20);
6101 delete fSeeds->RemoveAt(i);
6106 RemoveUsed2(fSeeds,0.85,0.85,0);
6107 if (AliTPCReconstructor::GetRecoParam()->GetDoKinks()) FindKinks(fSeeds,fEvent);
6108 //FindCurling(fSeeds, fEvent,0);
6109 if (AliTPCReconstructor::StreamLevel()>2) FindMultiMC(fSeeds, fEvent,-1); // find multi found tracks
6110 RemoveUsed2(fSeeds,0.5,0.4,20);
6111 FindSplitted(fSeeds, fEvent,0); // find multi found tracks
6112 if (AliTPCReconstructor::StreamLevel()>2) FindMultiMC(fSeeds, fEvent,0); // find multi found tracks
6115 // // refit short tracks
6117 Int_t nseed=fSeeds->GetEntriesFast();
6120 for (Int_t i=0; i<nseed; i++) {
6121 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
6123 Int_t nc=t.GetNumberOfClusters();
6125 delete fSeeds->RemoveAt(i);
6128 CookLabel(pt,0.1); //For comparison only
6129 //if ((pt->IsActive() || (pt->fRemoval==10) )&& nc>50 &&pt->GetNumberOfClusters()>0.4*pt->fNFoundable){
6130 if ((pt->IsActive() || (pt->GetRemoval()==10) )){
6132 if (fDebug>0) cerr<<found<<'\r';
6136 delete fSeeds->RemoveAt(i);
6140 //RemoveOverlap(fSeeds,0.99,7,kTRUE);
6142 //RemoveUsed(fSeeds,0.9,0.9,6);
6144 nseed=fSeeds->GetEntriesFast();
6146 for (Int_t i=0; i<nseed; i++) {
6147 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
6149 Int_t nc=t.GetNumberOfClusters();
6151 delete fSeeds->RemoveAt(i);
6155 t.CookdEdx(0.02,0.6);
6156 // CheckKinkPoint(&t,0.05);
6157 //if ((pt->IsActive() || (pt->fRemoval==10) )&& nc>50 &&pt->GetNumberOfClusters()>0.4*pt->fNFoundable){
6158 if ((pt->IsActive() || (pt->GetRemoval()==10) )){
6166 delete fSeeds->RemoveAt(i);
6167 //AliTPCseed * seed1 = ReSeed(pt,0.05,0.5,1);
6169 // FollowProlongation(*seed1,0);
6170 // Int_t n = seed1->GetNumberOfClusters();
6171 // printf("fP4\t%f\t%f\n",seed1->GetC(),pt->GetC());
6172 // printf("fN\t%d\t%d\n", seed1->GetNumberOfClusters(),pt->GetNumberOfClusters());
6175 //AliTPCseed * seed2 = ReSeed(pt,0.95,0.5,0.05);
6179 SortTracks(fSeeds, 1);
6183 PrepareForBackProlongation(fSeeds,5.);
6184 PropagateBack(fSeeds);
6185 printf("Time for back propagation: \t");timer.Print();timer.Start();
6189 PrepareForProlongation(fSeeds,5.);
6190 PropagateForward2(fSeeds);
6192 printf("Time for FORWARD propagation: \t");timer.Print();timer.Start();
6193 // RemoveUsed(fSeeds,0.7,0.7,6);
6194 //RemoveOverlap(fSeeds,0.9,7,kTRUE);
6196 nseed=fSeeds->GetEntriesFast();
6198 for (Int_t i=0; i<nseed; i++) {
6199 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
6201 Int_t nc=t.GetNumberOfClusters();
6203 delete fSeeds->RemoveAt(i);
6206 t.CookdEdx(0.02,0.6);
6207 // CookLabel(pt,0.1); //For comparison only
6208 //if ((pt->IsActive() || (pt->fRemoval==10) )&& nc>50 &&pt->GetNumberOfClusters()>0.4*pt->fNFoundable){
6209 if ((pt->IsActive() || (pt->fRemoval==10) )){
6210 cerr<<found++<<'\r';
6213 delete fSeeds->RemoveAt(i);
6218 // fNTracks = found;
6220 Info("Clusters2Tracks","Time for overlap removal, track writing and dedx cooking: \t"); timer.Print();timer.Start();
6223 // cerr<<"Number of found tracks : "<<"\t"<<found<<endl;
6224 Info("Clusters2Tracks","Number of found tracks %d",found);
6226 // UnloadClusters();
6231 void AliTPCtrackerMI::Tracking(TObjArray * arr)
6234 // tracking of the seeds
6237 fSectors = fOuterSec;
6238 ParallelTracking(arr,150,63);
6239 fSectors = fOuterSec;
6240 ParallelTracking(arr,63,0);
6243 TObjArray * AliTPCtrackerMI::Tracking(Int_t seedtype, Int_t i1, Int_t i2, Float_t cuts[4], Float_t dy, Int_t dsec)
6248 TObjArray * arr = new TObjArray;
6250 fSectors = fOuterSec;
6253 for (Int_t sec=0;sec<fkNOS;sec++){
6254 if (seedtype==3) MakeSeeds3(arr,sec,i1,i2,cuts,dy, dsec);
6255 if (seedtype==4) MakeSeeds5(arr,sec,i1,i2,cuts,dy);
6256 if (seedtype==2) MakeSeeds2(arr,sec,i1,i2,cuts,dy);
6259 Info("Tracking","\nSeeding - %d\t%d\t%d\t%d\n",seedtype,i1,i2,arr->GetEntriesFast());
6271 TObjArray * AliTPCtrackerMI::Tracking()
6275 if (AliTPCReconstructor::GetRecoParam()->GetSpecialSeeding()) return TrackingSpecial();
6278 Int_t nup=fOuterSec->GetNRows()+fInnerSec->GetNRows();
6280 TObjArray * seeds = new TObjArray;
6289 Float_t fnumber = 3.0;
6290 Float_t fdensity = 3.0;
6295 for (Int_t delta = 0; delta<18; delta+=6){
6299 arr = Tracking(3,nup-1-delta,nup-1-delta-gap,cuts,-1,1);
6300 SumTracks(seeds,arr);
6301 SignClusters(seeds,fnumber,fdensity);
6303 for (Int_t i=2;i<6;i+=2){
6304 // seed high pt tracks
6307 arr = Tracking(3,nup-i-delta,nup-i-delta-gap,cuts,-1,0);
6308 SumTracks(seeds,arr);
6309 SignClusters(seeds,fnumber,fdensity);
6314 // RemoveUsed(seeds,0.9,0.9,1);
6315 // UnsignClusters();
6316 // SignClusters(seeds,fnumber,fdensity);
6320 for (Int_t delta = 20; delta<120; delta+=10){
6322 // seed high pt tracks
6326 arr = Tracking(3,nup-delta,nup-delta-gap,cuts,-1);
6327 SumTracks(seeds,arr);
6328 SignClusters(seeds,fnumber,fdensity);
6333 arr = Tracking(3,nup-delta-5,nup-delta-5-gap,cuts,-1);
6334 SumTracks(seeds,arr);
6335 SignClusters(seeds,fnumber,fdensity);
6346 Info("Tracking()","\n\nPrimary seeding\t%d\n\n",seeds->GetEntriesFast());
6350 // RemoveUsed(seeds,0.75,0.75,1);
6352 //SignClusters(seeds,fnumber,fdensity);
6361 arr = Tracking(4,nup-1,nup-1-gap,cuts,-1);
6362 SumTracks(seeds,arr);
6363 SignClusters(seeds,fnumber,fdensity);
6365 arr = Tracking(4,nup-2,nup-2-gap,cuts,-1);
6366 SumTracks(seeds,arr);
6367 SignClusters(seeds,fnumber,fdensity);
6369 arr = Tracking(4,nup-3,nup-3-gap,cuts,-1);
6370 SumTracks(seeds,arr);
6371 SignClusters(seeds,fnumber,fdensity);
6375 for (Int_t delta = 3; delta<30; delta+=5){
6381 arr = Tracking(4,nup-1-delta,nup-1-delta-gap,cuts,-1);
6382 SumTracks(seeds,arr);
6383 SignClusters(seeds,fnumber,fdensity);
6385 arr = Tracking(4,nup-3-delta,nup-5-delta-gap,cuts,4);
6386 SumTracks(seeds,arr);
6387 SignClusters(seeds,fnumber,fdensity);
6398 Int_t fLastSeedRowSec=AliTPCReconstructor::GetRecoParam()->GetLastSeedRowSec();
6401 for (Int_t delta = 30; delta<fLastSeedRowSec; delta+=10){
6407 arr = Tracking(4,nup-1-delta,nup-1-delta-gap,cuts,-1);
6408 SumTracks(seeds,arr);
6409 SignClusters(seeds,fnumber,fdensity);
6411 arr = Tracking(4,nup-5-delta,nup-5-delta-gap,cuts,5 );
6412 SumTracks(seeds,arr);
6413 SignClusters(seeds,fnumber,fdensity);
6417 Info("Tracking()","\n\nSecondary seeding\t%d\n\n",seeds->GetEntriesFast());
6428 TObjArray * AliTPCtrackerMI::TrackingSpecial()
6431 // seeding adjusted for laser and cosmic tests - short tracks with big inclination angle
6432 // no primary vertex seeding tried
6436 Int_t nup=fOuterSec->GetNRows()+fInnerSec->GetNRows();
6438 TObjArray * seeds = new TObjArray;
6443 Float_t fnumber = 3.0;
6444 Float_t fdensity = 3.0;
6447 cuts[0] = AliTPCReconstructor::GetRecoParam()->GetMaxC(); // max curvature
6448 cuts[1] = 3.5; // max tan(phi) angle for seeding
6449 cuts[2] = 3.; // not used (cut on z primary vertex)
6450 cuts[3] = 3.5; // max tan(theta) angle for seeding
6452 for (Int_t delta = 0; nup-delta-gap-1>0; delta+=3){
6454 arr = Tracking(4,nup-1-delta,nup-1-delta-gap,cuts,-1);
6455 SumTracks(seeds,arr);
6456 SignClusters(seeds,fnumber,fdensity);
6460 Info("Tracking()","\n\nSecondary seeding\t%d\n\n",seeds->GetEntriesFast());
6471 void AliTPCtrackerMI::SumTracks(TObjArray *arr1,TObjArray *&arr2) const
6474 //sum tracks to common container
6475 //remove suspicious tracks
6476 Int_t nseed = arr2->GetEntriesFast();
6477 for (Int_t i=0;i<nseed;i++){
6478 AliTPCseed *pt=(AliTPCseed*)arr2->UncheckedAt(i);
6481 // remove tracks with too big curvature
6483 if (TMath::Abs(pt->GetC())>AliTPCReconstructor::GetRecoParam()->GetMaxC()){
6484 delete arr2->RemoveAt(i);
6487 // REMOVE VERY SHORT TRACKS
6488 if (pt->GetNumberOfClusters()<20){
6489 delete arr2->RemoveAt(i);
6492 // NORMAL ACTIVE TRACK
6493 if (pt->IsActive()){
6494 arr1->AddLast(arr2->RemoveAt(i));
6497 //remove not usable tracks
6498 if (pt->GetRemoval()!=10){
6499 delete arr2->RemoveAt(i);
6503 // ENABLE ONLY ENOUGH GOOD STOPPED TRACKS
6504 if (pt->GetDensityFirst(20)>0.8 || pt->GetDensityFirst(30)>0.8 || pt->GetDensityFirst(40)>0.7)
6505 arr1->AddLast(arr2->RemoveAt(i));
6507 delete arr2->RemoveAt(i);
6511 delete arr2; arr2 = 0;
6516 void AliTPCtrackerMI::ParallelTracking(TObjArray * arr, Int_t rfirst, Int_t rlast)
6519 // try to track in parralel
6521 Int_t nseed=arr->GetEntriesFast();
6522 //prepare seeds for tracking
6523 for (Int_t i=0; i<nseed; i++) {
6524 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i), &t=*pt;
6526 if (!t.IsActive()) continue;
6527 // follow prolongation to the first layer
6528 if ( (fSectors ==fInnerSec) || (t.GetFirstPoint()-fkParam->GetNRowLow()>rfirst+1) )
6529 FollowProlongation(t, rfirst+1);
6534 for (Int_t nr=rfirst; nr>=rlast; nr--){
6535 if (nr<fInnerSec->GetNRows())
6536 fSectors = fInnerSec;
6538 fSectors = fOuterSec;
6539 // make indexes with the cluster tracks for given
6541 // find nearest cluster
6542 for (Int_t i=0; i<nseed; i++) {
6543 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i), &t=*pt;
6545 if (nr==80) pt->UpdateReference();
6546 if (!pt->IsActive()) continue;
6547 // if ( (fSectors ==fOuterSec) && (pt->fFirstPoint-fkParam->GetNRowLow())<nr) continue;
6548 if (pt->GetRelativeSector()>17) {
6551 UpdateClusters(t,nr);
6553 // prolonagate to the nearest cluster - if founded
6554 for (Int_t i=0; i<nseed; i++) {
6555 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
6557 if (!pt->IsActive()) continue;
6558 // if ((fSectors ==fOuterSec) && (pt->fFirstPoint-fkParam->GetNRowLow())<nr) continue;
6559 if (pt->GetRelativeSector()>17) {
6562 FollowToNextCluster(*pt,nr);
6567 void AliTPCtrackerMI::PrepareForBackProlongation(TObjArray * arr,Float_t fac) const
6571 // if we use TPC track itself we have to "update" covariance
6573 Int_t nseed= arr->GetEntriesFast();
6574 for (Int_t i=0;i<nseed;i++){
6575 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
6579 //rotate to current local system at first accepted point
6580 Int_t index = pt->GetClusterIndex2(pt->GetFirstPoint());
6581 Int_t sec = (index&0xff000000)>>24;
6583 Float_t angle1 = fInnerSec->GetAlpha()*sec+fInnerSec->GetAlphaShift();
6584 if (angle1>TMath::Pi())
6585 angle1-=2.*TMath::Pi();
6586 Float_t angle2 = pt->GetAlpha();
6588 if (TMath::Abs(angle1-angle2)>0.001){
6589 pt->Rotate(angle1-angle2);
6590 //angle2 = pt->GetAlpha();
6591 //pt->fRelativeSector = pt->GetAlpha()/fInnerSec->GetAlpha();
6592 //if (pt->GetAlpha()<0)
6593 // pt->fRelativeSector+=18;
6594 //sec = pt->fRelativeSector;
6603 void AliTPCtrackerMI::PrepareForProlongation(TObjArray * arr, Float_t fac) const
6607 // if we use TPC track itself we have to "update" covariance
6609 Int_t nseed= arr->GetEntriesFast();
6610 for (Int_t i=0;i<nseed;i++){
6611 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
6614 pt->SetFirstPoint(pt->GetLastPoint());
6622 Int_t AliTPCtrackerMI::PropagateBack(TObjArray * arr)
6625 // make back propagation
6627 Int_t nseed= arr->GetEntriesFast();
6628 for (Int_t i=0;i<nseed;i++){
6629 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
6630 if (pt&& pt->GetKinkIndex(0)<=0) {
6631 //AliTPCseed *pt2 = new AliTPCseed(*pt);
6632 fSectors = fInnerSec;
6633 //FollowBackProlongation(*pt,fInnerSec->GetNRows()-1);
6634 //fSectors = fOuterSec;
6635 FollowBackProlongation(*pt,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1);
6636 //if (pt->GetNumberOfClusters()<(pt->fEsd->GetTPCclusters(0)) ){
6637 // Error("PropagateBack","Not prolonged track %d",pt->GetLabel());
6638 // FollowBackProlongation(*pt2,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1);
6641 if (pt&& pt->GetKinkIndex(0)>0) {
6642 AliESDkink * kink = fEvent->GetKink(pt->GetKinkIndex(0)-1);
6643 pt->SetFirstPoint(kink->GetTPCRow0());
6644 fSectors = fInnerSec;
6645 FollowBackProlongation(*pt,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1);
6653 Int_t AliTPCtrackerMI::PropagateForward2(TObjArray * arr)
6656 // make forward propagation
6658 Int_t nseed= arr->GetEntriesFast();
6660 for (Int_t i=0;i<nseed;i++){
6661 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
6663 FollowProlongation(*pt,0);
6672 Int_t AliTPCtrackerMI::PropagateForward()
6675 // propagate track forward
6677 Int_t nseed = fSeeds->GetEntriesFast();
6678 for (Int_t i=0;i<nseed;i++){
6679 AliTPCseed *pt = (AliTPCseed*)fSeeds->UncheckedAt(i);
6681 AliTPCseed &t = *pt;
6682 Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
6683 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
6684 if (alpha < 0. ) alpha += 2.*TMath::Pi();
6685 t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN);
6689 fSectors = fOuterSec;
6690 ParallelTracking(fSeeds,fOuterSec->GetNRows()+fInnerSec->GetNRows()-1,fInnerSec->GetNRows());
6691 fSectors = fInnerSec;
6692 ParallelTracking(fSeeds,fInnerSec->GetNRows()-1,0);
6701 Int_t AliTPCtrackerMI::PropagateBack(AliTPCseed * pt, Int_t row0, Int_t row1)
6704 // make back propagation, in between row0 and row1
6708 fSectors = fInnerSec;
6711 if (row1<fSectors->GetNRows())
6714 r1 = fSectors->GetNRows()-1;
6716 if (row0<fSectors->GetNRows()&& r1>0 )
6717 FollowBackProlongation(*pt,r1);
6718 if (row1<=fSectors->GetNRows())
6721 r1 = row1 - fSectors->GetNRows();
6722 if (r1<=0) return 0;
6723 if (r1>=fOuterSec->GetNRows()) return 0;
6724 fSectors = fOuterSec;
6725 return FollowBackProlongation(*pt,r1);
6733 void AliTPCtrackerMI::GetShape(AliTPCseed * seed, Int_t row)
6737 AliTPCClusterParam * clparam = AliTPCcalibDB::Instance()->GetClusterParam();
6738 Float_t zdrift = TMath::Abs((fkParam->GetZLength(0)-TMath::Abs(seed->GetZ())));
6739 Int_t type = (seed->GetSector() < fkParam->GetNSector()/2) ? 0: (row>126) ? 1:2;
6740 Double_t angulary = seed->GetSnp();
6741 angulary = angulary*angulary/(1.-angulary*angulary);
6742 Double_t angularz = seed->GetTgl()*seed->GetTgl()*(1.+angulary);
6744 Double_t sigmay = clparam->GetRMS0(0,type,zdrift,TMath::Sqrt(TMath::Abs(angulary)));
6745 Double_t sigmaz = clparam->GetRMS0(1,type,zdrift,TMath::Sqrt(TMath::Abs(angularz)));
6746 seed->SetCurrentSigmaY2(sigmay*sigmay);
6747 seed->SetCurrentSigmaZ2(sigmaz*sigmaz);
6748 // Float_t sd2 = TMath::Abs((fkParam->GetZLength(0)-TMath::Abs(seed->GetZ())))*fkParam->GetDiffL()*fkParam->GetDiffL();
6749 // // Float_t padlength = fkParam->GetPadPitchLength(seed->fSector);
6750 // Float_t padlength = GetPadPitchLength(row);
6752 // Float_t sresy = (seed->GetSector() < fkParam->GetNSector()/2) ? 0.2 :0.3;
6753 // seed->SetCurrentSigmaY2(sd2+padlength*padlength*angulary/12.+sresy*sresy);
6755 // Float_t sresz = fkParam->GetZSigma();
6756 // seed->SetCurrentSigmaZ2(sd2+padlength*padlength*angularz*angularz*(1+angulary)/12.+sresz*sresz);
6758 Float_t wy = GetSigmaY(seed);
6759 Float_t wz = GetSigmaZ(seed);
6762 if (TMath::Abs(wy/seed->fCurrentSigmaY2-1)>0.0001 || TMath::Abs(wz/seed->fCurrentSigmaZ2-1)>0.0001 ){
6763 printf("problem\n");
6770 //__________________________________________________________________________
6771 void AliTPCtrackerMI::CookLabel(AliKalmanTrack *tk, Float_t wrong) const {
6772 //--------------------------------------------------------------------
6773 //This function "cooks" a track label. If label<0, this track is fake.
6774 //--------------------------------------------------------------------
6775 AliTPCseed * t = dynamic_cast<AliTPCseed*>(tk);
6777 printf("%s:%d wrong type \n",(char*)__FILE__,__LINE__);
6781 Int_t noc=t->GetNumberOfClusters();
6783 //printf("\nnot founded prolongation\n\n\n");
6789 AliTPCclusterMI *clusters[160];
6791 for (Int_t i=0;i<160;i++) {
6798 for (i=0; i<160 && current<noc; i++) {
6800 Int_t index=t->GetClusterIndex2(i);
6801 if (index<=0) continue;
6802 if (index&0x8000) continue;
6804 //clusters[current]=GetClusterMI(index);
6805 if (t->GetClusterPointer(i)){
6806 clusters[current]=t->GetClusterPointer(i);
6812 Int_t lab=123456789;
6813 for (i=0; i<noc; i++) {
6814 AliTPCclusterMI *c=clusters[i];
6816 lab=TMath::Abs(c->GetLabel(0));
6818 for (j=0; j<noc; j++) if (lb[j]==lab || mx[j]==0) break;
6824 for (i=0; i<noc; i++) if (mx[i]>max) {max=mx[i]; lab=lb[i];}
6826 for (i=0; i<noc; i++) {
6827 AliTPCclusterMI *c=clusters[i];
6829 if (TMath::Abs(c->GetLabel(1)) == lab ||
6830 TMath::Abs(c->GetLabel(2)) == lab ) max++;
6833 if ((1.- Float_t(max)/noc) > wrong) lab=-lab;
6836 Int_t tail=Int_t(0.10*noc);
6839 for (i=1; i<=160&&ind<tail; i++) {
6840 // AliTPCclusterMI *c=clusters[noc-i];
6841 AliTPCclusterMI *c=clusters[i];
6843 if (lab == TMath::Abs(c->GetLabel(0)) ||
6844 lab == TMath::Abs(c->GetLabel(1)) ||
6845 lab == TMath::Abs(c->GetLabel(2))) max++;
6848 if (max < Int_t(0.5*tail)) lab=-lab;
6855 //delete[] clusters;
6859 //__________________________________________________________________________
6860 Int_t AliTPCtrackerMI::CookLabel(AliTPCseed *t, Float_t wrong,Int_t first, Int_t last) const {
6861 //--------------------------------------------------------------------
6862 //This function "cooks" a track label. If label<0, this track is fake.
6863 //--------------------------------------------------------------------
6864 Int_t noc=t->GetNumberOfClusters();
6866 //printf("\nnot founded prolongation\n\n\n");
6872 AliTPCclusterMI *clusters[160];
6874 for (Int_t i=0;i<160;i++) {
6881 for (i=0; i<160 && current<noc; i++) {
6882 if (i<first) continue;
6883 if (i>last) continue;
6884 Int_t index=t->GetClusterIndex2(i);
6885 if (index<=0) continue;
6886 if (index&0x8000) continue;
6888 //clusters[current]=GetClusterMI(index);
6889 if (t->GetClusterPointer(i)){
6890 clusters[current]=t->GetClusterPointer(i);
6895 if (noc<5) return -1;
6896 Int_t lab=123456789;
6897 for (i=0; i<noc; i++) {
6898 AliTPCclusterMI *c=clusters[i];
6900 lab=TMath::Abs(c->GetLabel(0));
6902 for (j=0; j<noc; j++) if (lb[j]==lab || mx[j]==0) break;
6908 for (i=0; i<noc; i++) if (mx[i]>max) {max=mx[i]; lab=lb[i];}
6910 for (i=0; i<noc; i++) {
6911 AliTPCclusterMI *c=clusters[i];
6913 if (TMath::Abs(c->GetLabel(1)) == lab ||
6914 TMath::Abs(c->GetLabel(2)) == lab ) max++;
6917 if ((1.- Float_t(max)/noc) > wrong) lab=-lab;
6920 Int_t tail=Int_t(0.10*noc);
6923 for (i=1; i<=160&&ind<tail; i++) {
6924 // AliTPCclusterMI *c=clusters[noc-i];
6925 AliTPCclusterMI *c=clusters[i];
6927 if (lab == TMath::Abs(c->GetLabel(0)) ||
6928 lab == TMath::Abs(c->GetLabel(1)) ||
6929 lab == TMath::Abs(c->GetLabel(2))) max++;
6932 if (max < Int_t(0.5*tail)) lab=-lab;
6935 // t->SetLabel(lab);
6939 //delete[] clusters;
6943 Int_t AliTPCtrackerMI::GetRowNumber(Double_t x[3]) const
6945 //return pad row number for given x vector
6946 Float_t phi = TMath::ATan2(x[1],x[0]);
6947 if(phi<0) phi=2.*TMath::Pi()+phi;
6948 // Get the local angle in the sector philoc
6949 const Float_t kRaddeg = 180/3.14159265358979312;
6950 Float_t phiangle = (Int_t (phi*kRaddeg/20.) + 0.5)*20./kRaddeg;
6951 Double_t localx = x[0]*TMath::Cos(phiangle)-x[1]*TMath::Sin(phiangle);
6952 return GetRowNumber(localx);
6957 void AliTPCtrackerMI::MakeBitmaps(AliTPCseed *t)
6959 //-----------------------------------------------------------------------
6960 // Fill the cluster and sharing bitmaps of the track
6961 //-----------------------------------------------------------------------
6963 Int_t firstpoint = 0;
6964 Int_t lastpoint = 159;
6965 AliTPCTrackerPoint *point;
6966 AliTPCclusterMI *cluster;
6968 for (int iter=firstpoint; iter<lastpoint; iter++) {
6969 // Change to cluster pointers to see if we have a cluster at given padrow
6970 cluster = t->GetClusterPointer(iter);
6972 t->SetClusterMapBit(iter, kTRUE);
6973 point = t->GetTrackPoint(iter);
6974 if (point->IsShared())
6975 t->SetSharedMapBit(iter,kTRUE);
6977 t->SetSharedMapBit(iter, kFALSE);
6980 t->SetClusterMapBit(iter, kFALSE);
6981 t->SetSharedMapBit(iter, kFALSE);
6986 Bool_t AliTPCtrackerMI::IsFindable(AliTPCseed & track){
6988 // return flag if there is findable cluster at given position
6991 Float_t z = track.GetZ();
6993 if (TMath::Abs(z)<(AliTPCReconstructor::GetCtgRange()*track.GetX()+kDeltaZ) &&
6994 TMath::Abs(z)<fkParam->GetZLength(0) &&
6995 (TMath::Abs(track.GetSnp())<AliTPCReconstructor::GetMaxSnpTracker()))
7001 void AliTPCtrackerMI::AddCovariance(AliTPCseed * seed){
7003 // Adding systematic error
7004 // !!!! the systematic error for element 4 is in 1/cm not in pt
7006 const Double_t *param = AliTPCReconstructor::GetRecoParam()->GetSystematicError();
7008 for (Int_t i=0;i<15;i++) covar[i]=0;
7014 covar[0] = param[0]*param[0];
7015 covar[2] = param[1]*param[1];
7016 covar[5] = param[2]*param[2];
7017 covar[9] = param[3]*param[3];
7018 Double_t facC = AliTracker::GetBz()*kB2C;
7019 covar[14]= param[4]*param[4]*facC*facC;
7020 seed->AddCovariance(covar);