1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
17 //-------------------------------------------------------
18 // Implementation of the TPC tracker
20 // Origin: Marian Ivanov Marian.Ivanov@cern.ch
22 // AliTPC parallel tracker
24 // The track fitting is based on Kalaman filtering approach
26 // The track finding steps:
27 // 1. Seeding - with and without vertex constraint
28 // - seeding with vertex constain done at first n^2 proble
29 // - seeding without vertex constraint n^3 problem
30 // 2. Tracking - follow prolongation road - find cluster - update kalman track
32 // The seeding and tracking is repeated several times, in different seeding region.
33 // This approach enables to find the track which cannot be seeded in some region of TPC
34 // This can happen because of low momenta (track do not reach outer radius), or track is currently in the ded region between sectors, or the track is for the moment overlapped with other track (seed quality is poor) ...
36 // With this approach we reach almost 100 % efficiency also for high occupancy events.
37 // (If the seeding efficiency in a region is about 90 % than with logical or of several
38 // regions we will reach 100% (in theory - supposing independence)
40 // Repeating several seeding - tracking procedures some of the tracks can be find
43 // The procedures to remove multi find tacks are impremented:
44 // RemoveUsed2 - fast procedure n problem -
45 // Algorithm - Sorting tracks according quality
46 // remove tracks with some shared fraction
47 // Sharing in respect to all tacks
48 // Signing clusters in gold region
49 // FindSplitted - slower algorithm n^2
50 // Sort the tracks according quality
51 // Loop over pair of tracks
52 // If overlap with other track bigger than threshold - remove track
54 // FindCurling - Finds the pair of tracks which are curling
55 // - About 10% of tracks can be find with this procedure
56 // The combinatorial background is too big to be used in High
57 // multiplicity environment
58 // - n^2 problem - Slow procedure - currently it is disabled because of
61 // The number of splitted tracks can be reduced disabling the sharing of the cluster.
62 // tpcRecoParam-> SetClusterSharing(kFALSE);
63 // IT IS HIGHLY non recomended to use it in high flux enviroonment
64 // Even using this switch some tracks can be found more than once
65 // (because of multiple seeding and low quality tracks which will not cross full chamber)
68 // The tracker itself can be debugged - the information about tracks can be stored in several // phases of the reconstruction
69 // To enable storage of the TPC tracks in the ESD friend track
70 // use AliTPCReconstructor::SetStreamLevel(n); where nis bigger 0
72 // The debug level - different procedure produce tree for numerical debugging
73 // To enable them set AliTPCReconstructor::SetStreamLevel(n); where nis bigger 1
77 // Adding systematic errors to the covariance:
79 // The systematic errors due to the misalignment and miscalibration are added to the covariance matrix
80 // of the tracks (not to the clusters as they are dependent):
81 // The parameters form AliTPCRecoParam are used AliTPCRecoParam::GetSystematicError
82 // The systematic errors are expressed there in RMS - position (cm), angle (rad), curvature (1/cm)
83 // The default values are 0.
85 // The sytematic errors are added to the covariance matrix in following places:
87 // 1. During fisrt itteration - AliTPCtrackerMI::FillESD
88 // 2. Second iteration -
89 // 2.a ITS->TPC - AliTPCtrackerMI::ReadSeeds
90 // 2.b TPC->TRD - AliTPCtrackerMI::PropagateBack
91 // 3. Third iteration -
92 // 3.a TRD->TPC - AliTPCtrackerMI::ReadSeeds
93 // 3.b TPC->ITS - AliTPCtrackerMI::RefitInward
95 // There are several places in the code which can be numerically debuged
96 // This code is keeped in order to enable code development and to check the calibration implementtion
98 // 1. ErrParam stream (Log level 9) - dump information about
100 // 2.a) cluster error estimate
101 // 3.a) cluster shape estimate
104 //-------------------------------------------------------
109 #include "Riostream.h"
110 #include <TClonesArray.h>
112 #include <TObjArray.h>
115 #include "AliComplexCluster.h"
116 #include "AliESDEvent.h"
117 #include "AliESDtrack.h"
118 #include "AliESDVertex.h"
121 #include "AliHelix.h"
122 #include "AliRunLoader.h"
123 #include "AliTPCClustersRow.h"
124 #include "AliTPCParam.h"
125 #include "AliTPCReconstructor.h"
126 #include "AliTPCpolyTrack.h"
127 #include "AliTPCreco.h"
128 #include "AliTPCseed.h"
130 #include "AliTPCtrackerSector.h"
131 #include "AliTPCtrackerMI.h"
132 #include "TStopwatch.h"
133 #include "AliTPCReconstructor.h"
135 #include "TTreeStream.h"
136 #include "AliAlignObj.h"
137 #include "AliTrackPointArray.h"
139 #include "AliTPCcalibDB.h"
140 #include "AliTPCTransform.h"
141 #include "AliTPCClusterParam.h"
145 ClassImp(AliTPCtrackerMI)
149 class AliTPCFastMath {
152 static Double_t FastAsin(Double_t x);
154 static Double_t fgFastAsin[20000]; //lookup table for fast asin computation
157 Double_t AliTPCFastMath::fgFastAsin[20000];
158 AliTPCFastMath gAliTPCFastMath; // needed to fill the LUT
160 AliTPCFastMath::AliTPCFastMath(){
162 // initialized lookup table;
163 for (Int_t i=0;i<10000;i++){
164 fgFastAsin[2*i] = TMath::ASin(i/10000.);
165 fgFastAsin[2*i+1] = (TMath::ASin((i+1)/10000.)-fgFastAsin[2*i]);
169 Double_t AliTPCFastMath::FastAsin(Double_t x){
171 // return asin using lookup table
173 Int_t index = int(x*10000);
174 return fgFastAsin[2*index]+(x*10000.-index)*fgFastAsin[2*index+1];
177 Int_t index = int(x*10000);
178 return -(fgFastAsin[2*index]+(x*10000.-index)*fgFastAsin[2*index+1]);
180 //__________________________________________________________________
181 AliTPCtrackerMI::AliTPCtrackerMI()
203 // default constructor
206 //_____________________________________________________________________
210 Int_t AliTPCtrackerMI::UpdateTrack(AliTPCseed * track, Int_t accept){
212 //update track information using current cluster - track->fCurrentCluster
215 AliTPCclusterMI* c =track->GetCurrentCluster();
216 if (accept>0) track->SetCurrentClusterIndex1(track->GetCurrentClusterIndex1() | 0x8000); //sign not accepted clusters
218 UInt_t i = track->GetCurrentClusterIndex1();
220 Int_t sec=(i&0xff000000)>>24;
221 //Int_t row = (i&0x00ff0000)>>16;
222 track->SetRow((i&0x00ff0000)>>16);
223 track->SetSector(sec);
224 // Int_t index = i&0xFFFF;
225 if (sec>=fParam->GetNInnerSector()) track->SetRow(track->GetRow()+fParam->GetNRowLow());
226 track->SetClusterIndex2(track->GetRow(), i);
227 //track->fFirstPoint = row;
228 //if ( track->fLastPoint<row) track->fLastPoint =row;
229 // if (track->fRow<0 || track->fRow>160) {
230 // printf("problem\n");
232 if (track->GetFirstPoint()>track->GetRow())
233 track->SetFirstPoint(track->GetRow());
234 if (track->GetLastPoint()<track->GetRow())
235 track->SetLastPoint(track->GetRow());
238 track->SetClusterPointer(track->GetRow(),c);
241 Double_t angle2 = track->GetSnp()*track->GetSnp();
243 //SET NEW Track Point
245 if (angle2<1) //PH sometimes angle2 is very big. To be investigated...
247 angle2 = TMath::Sqrt(angle2/(1-angle2));
248 AliTPCTrackerPoint &point =*(track->GetTrackPoint(track->GetRow()));
250 point.SetSigmaY(c->GetSigmaY2()/track->GetCurrentSigmaY2());
251 point.SetSigmaZ(c->GetSigmaZ2()/track->GetCurrentSigmaZ2());
252 point.SetErrY(sqrt(track->GetErrorY2()));
253 point.SetErrZ(sqrt(track->GetErrorZ2()));
255 point.SetX(track->GetX());
256 point.SetY(track->GetY());
257 point.SetZ(track->GetZ());
258 point.SetAngleY(angle2);
259 point.SetAngleZ(track->GetTgl());
260 if (point.IsShared()){
261 track->SetErrorY2(track->GetErrorY2()*4);
262 track->SetErrorZ2(track->GetErrorZ2()*4);
266 Double_t chi2 = track->GetPredictedChi2(track->GetCurrentCluster());
268 // track->SetErrorY2(track->GetErrorY2()*1.3);
269 // track->SetErrorY2(track->GetErrorY2()+0.01);
270 // track->SetErrorZ2(track->GetErrorZ2()*1.3);
271 // track->SetErrorZ2(track->GetErrorZ2()+0.005);
273 if (accept>0) return 0;
274 if (track->GetNumberOfClusters()%20==0){
275 // if (track->fHelixIn){
276 // TClonesArray & larr = *(track->fHelixIn);
277 // Int_t ihelix = larr.GetEntriesFast();
278 // new(larr[ihelix]) AliHelix(*track) ;
281 track->SetNoCluster(0);
282 return track->Update(c,chi2,i);
287 Int_t AliTPCtrackerMI::AcceptCluster(AliTPCseed * seed, AliTPCclusterMI * cluster)
290 // decide according desired precision to accept given
291 // cluster for tracking
292 Double_t sy2=ErrY2(seed,cluster);
293 Double_t sz2=ErrZ2(seed,cluster);
295 Double_t sdistancey2 = sy2+seed->GetSigmaY2();
296 Double_t sdistancez2 = sz2+seed->GetSigmaZ2();
298 Double_t rdistancey2 = (seed->GetCurrentCluster()->GetY()-seed->GetY())*
299 (seed->GetCurrentCluster()->GetY()-seed->GetY())/sdistancey2;
300 Double_t rdistancez2 = (seed->GetCurrentCluster()->GetZ()-seed->GetZ())*
301 (seed->GetCurrentCluster()->GetZ()-seed->GetZ())/sdistancez2;
303 Double_t rdistance2 = rdistancey2+rdistancez2;
306 if (AliTPCReconstructor::StreamLevel()>5 && seed->GetNumberOfClusters()>20) {
307 Float_t rmsy2 = seed->GetCurrentSigmaY2();
308 Float_t rmsz2 = seed->GetCurrentSigmaZ2();
309 Float_t rmsy2p30 = seed->GetCMeanSigmaY2p30();
310 Float_t rmsz2p30 = seed->GetCMeanSigmaZ2p30();
311 Float_t rmsy2p30R = seed->GetCMeanSigmaY2p30R();
312 Float_t rmsz2p30R = seed->GetCMeanSigmaZ2p30R();
314 AliExternalTrackParam param(*seed);
315 static TVectorD gcl(3),gtr(3);
317 param.GetXYZ(gcl.GetMatrixArray());
318 cluster->GetGlobalXYZ(gclf);
319 gcl[0]=gclf[0]; gcl[1]=gclf[1]; gcl[2]=gclf[2];
321 if (AliTPCReconstructor::StreamLevel()>0) {
322 (*fDebugStreamer)<<"ErrParam"<<
331 "rmsy2p30="<<rmsy2p30<<
332 "rmsz2p30="<<rmsz2p30<<
333 "rmsy2p30R="<<rmsy2p30R<<
334 "rmsz2p30R="<<rmsz2p30R<<
339 if (rdistance2>16) return 3;
342 if ((rdistancey2>9. || rdistancez2>9.) && cluster->GetType()==0)
343 return 2; //suspisiouce - will be changed
345 if ((rdistancey2>6.25 || rdistancez2>6.25) && cluster->GetType()>0)
346 // strict cut on overlaped cluster
347 return 2; //suspisiouce - will be changed
349 if ( (rdistancey2>1. || rdistancez2>6.25 )
350 && cluster->GetType()<0){
351 seed->SetNFoundable(seed->GetNFoundable()-1);
360 //_____________________________________________________________________________
361 AliTPCtrackerMI::AliTPCtrackerMI(const AliTPCParam *par):
363 fkNIS(par->GetNInnerSector()/2),
365 fkNOS(par->GetNOuterSector()/2),
382 //---------------------------------------------------------------------
383 // The main TPC tracker constructor
384 //---------------------------------------------------------------------
385 fInnerSec=new AliTPCtrackerSector[fkNIS];
386 fOuterSec=new AliTPCtrackerSector[fkNOS];
389 for (i=0; i<fkNIS; i++) fInnerSec[i].Setup(par,0);
390 for (i=0; i<fkNOS; i++) fOuterSec[i].Setup(par,1);
393 Int_t nrowlow = par->GetNRowLow();
394 Int_t nrowup = par->GetNRowUp();
397 for (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*/){
442 //------------------------------
444 //--------------------------------------------------------------
447 //_____________________________________________________________________________
448 AliTPCtrackerMI::~AliTPCtrackerMI() {
449 //------------------------------------------------------------------
450 // TPC tracker destructor
451 //------------------------------------------------------------------
458 if (fDebugStreamer) delete fDebugStreamer;
462 void AliTPCtrackerMI::FillESD(TObjArray* arr)
466 //fill esds using updated tracks
468 // write tracks to the event
469 // store index of the track
470 Int_t nseed=arr->GetEntriesFast();
471 //FindKinks(arr,fEvent);
472 for (Int_t i=0; i<nseed; i++) {
473 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
477 if (AliTPCReconstructor::StreamLevel()>1) {
478 (*fDebugStreamer)<<"Track0"<<
482 // pt->PropagateTo(fParam->GetInnerRadiusLow());
483 if (pt->GetKinkIndex(0)<=0){ //don't propagate daughter tracks
484 pt->PropagateTo(fParam->GetInnerRadiusLow());
487 if (( pt->GetPoints()[2]- pt->GetPoints()[0])>5 && pt->GetPoints()[3]>0.8){
489 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
490 iotrack.SetTPCPoints(pt->GetPoints());
491 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
492 iotrack.SetV0Indexes(pt->GetV0Indexes());
493 // iotrack.SetTPCpid(pt->fTPCr);
494 //iotrack.SetTPCindex(i);
495 fEvent->AddTrack(&iotrack);
499 if ( (pt->GetNumberOfClusters()>70)&& (Float_t(pt->GetNumberOfClusters())/Float_t(pt->GetNFoundable()))>0.55) {
501 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
502 iotrack.SetTPCPoints(pt->GetPoints());
503 //iotrack.SetTPCindex(i);
504 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
505 iotrack.SetV0Indexes(pt->GetV0Indexes());
506 // iotrack.SetTPCpid(pt->fTPCr);
507 fEvent->AddTrack(&iotrack);
511 // short tracks - maybe decays
513 if ( (pt->GetNumberOfClusters()>30) && (Float_t(pt->GetNumberOfClusters())/Float_t(pt->GetNFoundable()))>0.70) {
514 Int_t found,foundable,shared;
515 pt->GetClusterStatistic(0,60,found, foundable,shared,kFALSE);
516 if ( (found>20) && (pt->GetNShared()/float(pt->GetNumberOfClusters())<0.2)){
518 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
519 //iotrack.SetTPCindex(i);
520 iotrack.SetTPCPoints(pt->GetPoints());
521 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
522 iotrack.SetV0Indexes(pt->GetV0Indexes());
523 //iotrack.SetTPCpid(pt->fTPCr);
524 fEvent->AddTrack(&iotrack);
529 if ( (pt->GetNumberOfClusters()>20) && (Float_t(pt->GetNumberOfClusters())/Float_t(pt->GetNFoundable()))>0.8) {
530 Int_t found,foundable,shared;
531 pt->GetClusterStatistic(0,60,found, foundable,shared,kFALSE);
532 if (found<20) continue;
533 if (pt->GetNShared()/float(pt->GetNumberOfClusters())>0.2) continue;
536 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
537 iotrack.SetTPCPoints(pt->GetPoints());
538 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
539 iotrack.SetV0Indexes(pt->GetV0Indexes());
540 //iotrack.SetTPCpid(pt->fTPCr);
541 //iotrack.SetTPCindex(i);
542 fEvent->AddTrack(&iotrack);
545 // short tracks - secondaties
547 if ( (pt->GetNumberOfClusters()>30) ) {
548 Int_t found,foundable,shared;
549 pt->GetClusterStatistic(128,158,found, foundable,shared,kFALSE);
550 if ( (found>20) && (pt->GetNShared()/float(pt->GetNumberOfClusters())<0.2) &&float(found)/float(foundable)>0.8){
552 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
553 iotrack.SetTPCPoints(pt->GetPoints());
554 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
555 iotrack.SetV0Indexes(pt->GetV0Indexes());
556 //iotrack.SetTPCpid(pt->fTPCr);
557 //iotrack.SetTPCindex(i);
558 fEvent->AddTrack(&iotrack);
563 if ( (pt->GetNumberOfClusters()>15)) {
564 Int_t found,foundable,shared;
565 pt->GetClusterStatistic(138,158,found, foundable,shared,kFALSE);
566 if (found<15) continue;
567 if (foundable<=0) continue;
568 if (pt->GetNShared()/float(pt->GetNumberOfClusters())>0.2) continue;
569 if (float(found)/float(foundable)<0.8) continue;
572 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
573 iotrack.SetTPCPoints(pt->GetPoints());
574 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
575 iotrack.SetV0Indexes(pt->GetV0Indexes());
576 // iotrack.SetTPCpid(pt->fTPCr);
577 //iotrack.SetTPCindex(i);
578 fEvent->AddTrack(&iotrack);
583 printf("Number of filled ESDs-\t%d\n",fEvent->GetNumberOfTracks());
590 Double_t AliTPCtrackerMI::ErrY2(AliTPCseed* seed, AliTPCclusterMI * cl){
593 // Use calibrated cluster error from OCDB
595 AliTPCClusterParam * clparam = AliTPCcalibDB::Instance()->GetClusterParam();
597 Float_t z = TMath::Abs(fParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
598 Int_t ctype = cl->GetType();
599 Int_t type = (cl->GetRow()<63) ? 0: (cl->GetRow()>126) ? 1:2;
600 Double_t angle = seed->GetSnp()*seed->GetSnp();
601 angle = TMath::Sqrt(TMath::Abs(angle/(1.-angle)));
602 Double_t erry2 = clparam->GetError0Par(0,type, z,angle);
604 erry2+=0.5; // edge cluster
607 seed->SetErrorY2(erry2);
611 //calculate look-up table at the beginning
612 // static Bool_t ginit = kFALSE;
613 // static Float_t gnoise1,gnoise2,gnoise3;
614 // static Float_t ggg1[10000];
615 // static Float_t ggg2[10000];
616 // static Float_t ggg3[10000];
617 // static Float_t glandau1[10000];
618 // static Float_t glandau2[10000];
619 // static Float_t glandau3[10000];
621 // static Float_t gcor01[500];
622 // static Float_t gcor02[500];
623 // static Float_t gcorp[500];
627 // if (ginit==kFALSE){
628 // for (Int_t i=1;i<500;i++){
629 // Float_t rsigma = float(i)/100.;
630 // gcor02[i] = TMath::Max(0.78 +TMath::Exp(7.4*(rsigma-1.2)),0.6);
631 // gcor01[i] = TMath::Max(0.72 +TMath::Exp(3.36*(rsigma-1.2)),0.6);
632 // gcorp[i] = TMath::Max(TMath::Power((rsigma+0.5),1.5),1.2);
636 // for (Int_t i=3;i<10000;i++){
640 // Float_t amp = float(i);
641 // Float_t padlength =0.75;
642 // gnoise1 = 0.0004/padlength;
643 // Float_t nel = 0.268*amp;
644 // Float_t nprim = 0.155*amp;
645 // ggg1[i] = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.001*nel/(padlength*padlength))/nel;
646 // glandau1[i] = (2.+0.12*nprim)*0.5* (2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
647 // if (glandau1[i]>1) glandau1[i]=1;
648 // glandau1[i]*=padlength*padlength/12.;
652 // gnoise2 = 0.0004/padlength;
654 // nprim = 0.133*amp;
655 // ggg2[i] = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
656 // glandau2[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
657 // if (glandau2[i]>1) glandau2[i]=1;
658 // glandau2[i]*=padlength*padlength/12.;
663 // gnoise3 = 0.0004/padlength;
665 // nprim = 0.133*amp;
666 // ggg3[i] = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
667 // glandau3[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
668 // if (glandau3[i]>1) glandau3[i]=1;
669 // glandau3[i]*=padlength*padlength/12.;
677 // Int_t amp = int(TMath::Abs(cl->GetQ()));
679 // seed->SetErrorY2(1.);
683 // Float_t z = TMath::Abs(fParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
684 // Int_t ctype = cl->GetType();
685 // Float_t padlength= GetPadPitchLength(seed->GetRow());
686 // Double_t angle2 = seed->GetSnp()*seed->GetSnp();
687 // angle2 = angle2/(1-angle2);
689 // //cluster "quality"
690 // Int_t rsigmay = int(100.*cl->GetSigmaY2()/(seed->GetCurrentSigmaY2()));
693 // if (fSectors==fInnerSec){
694 // snoise2 = gnoise1;
695 // res = ggg1[amp]*z+glandau1[amp]*angle2;
696 // if (ctype==0) res *= gcor01[rsigmay];
699 // res*= gcorp[rsigmay];
703 // if (padlength<1.1){
704 // snoise2 = gnoise2;
705 // res = ggg2[amp]*z+glandau2[amp]*angle2;
706 // if (ctype==0) res *= gcor02[rsigmay];
709 // res*= gcorp[rsigmay];
713 // snoise2 = gnoise3;
714 // res = ggg3[amp]*z+glandau3[amp]*angle2;
715 // if (ctype==0) res *= gcor02[rsigmay];
718 // res*= gcorp[rsigmay];
725 // res*=2.4; // overestimate error 2 times
729 // if (res<2*snoise2)
732 // seed->SetErrorY2(res);
740 Double_t AliTPCtrackerMI::ErrZ2(AliTPCseed* seed, AliTPCclusterMI * cl){
743 // Use calibrated cluster error from OCDB
745 AliTPCClusterParam * clparam = AliTPCcalibDB::Instance()->GetClusterParam();
747 Float_t z = TMath::Abs(fParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
748 Int_t ctype = cl->GetType();
749 Int_t type = (cl->GetRow()<63) ? 0: (cl->GetRow()>126) ? 1:2;
751 Double_t angle2 = seed->GetSnp()*seed->GetSnp();
752 angle2 = seed->GetTgl()*seed->GetTgl()*(1+angle2/(1-angle2));
753 Double_t angle = TMath::Sqrt(TMath::Abs(angle2));
754 Double_t errz2 = clparam->GetError0Par(1,type, z,angle);
756 errz2+=0.5; // edge cluster
759 seed->SetErrorZ2(errz2);
765 // //seed->SetErrorY2(0.1);
767 // //calculate look-up table at the beginning
768 // static Bool_t ginit = kFALSE;
769 // static Float_t gnoise1,gnoise2,gnoise3;
770 // static Float_t ggg1[10000];
771 // static Float_t ggg2[10000];
772 // static Float_t ggg3[10000];
773 // static Float_t glandau1[10000];
774 // static Float_t glandau2[10000];
775 // static Float_t glandau3[10000];
777 // static Float_t gcor01[1000];
778 // static Float_t gcor02[1000];
779 // static Float_t gcorp[1000];
783 // if (ginit==kFALSE){
784 // for (Int_t i=1;i<1000;i++){
785 // Float_t rsigma = float(i)/100.;
786 // gcor02[i] = TMath::Max(0.81 +TMath::Exp(6.8*(rsigma-1.2)),0.6);
787 // gcor01[i] = TMath::Max(0.72 +TMath::Exp(2.04*(rsigma-1.2)),0.6);
788 // gcorp[i] = TMath::Max(TMath::Power((rsigma+0.5),1.5),1.2);
792 // for (Int_t i=3;i<10000;i++){
796 // Float_t amp = float(i);
797 // Float_t padlength =0.75;
798 // gnoise1 = 0.0004/padlength;
799 // Float_t nel = 0.268*amp;
800 // Float_t nprim = 0.155*amp;
801 // ggg1[i] = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.001*nel/(padlength*padlength))/nel;
802 // glandau1[i] = (2.+0.12*nprim)*0.5* (2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
803 // if (glandau1[i]>1) glandau1[i]=1;
804 // glandau1[i]*=padlength*padlength/12.;
808 // gnoise2 = 0.0004/padlength;
810 // nprim = 0.133*amp;
811 // ggg2[i] = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
812 // glandau2[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
813 // if (glandau2[i]>1) glandau2[i]=1;
814 // glandau2[i]*=padlength*padlength/12.;
819 // gnoise3 = 0.0004/padlength;
821 // nprim = 0.133*amp;
822 // ggg3[i] = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
823 // glandau3[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
824 // if (glandau3[i]>1) glandau3[i]=1;
825 // glandau3[i]*=padlength*padlength/12.;
833 // Int_t amp = int(TMath::Abs(cl->GetQ()));
835 // seed->SetErrorY2(1.);
839 // Float_t z = TMath::Abs(fParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
840 // Int_t ctype = cl->GetType();
841 // Float_t padlength= GetPadPitchLength(seed->GetRow());
843 // Double_t angle2 = seed->GetSnp()*seed->GetSnp();
844 // // if (angle2<0.6) angle2 = 0.6;
845 // angle2 = seed->GetTgl()*seed->GetTgl()*(1+angle2/(1-angle2));
847 // //cluster "quality"
848 // Int_t rsigmaz = int(100.*cl->GetSigmaZ2()/(seed->GetCurrentSigmaZ2()));
851 // if (fSectors==fInnerSec){
852 // snoise2 = gnoise1;
853 // res = ggg1[amp]*z+glandau1[amp]*angle2;
854 // if (ctype==0) res *= gcor01[rsigmaz];
857 // res*= gcorp[rsigmaz];
861 // if (padlength<1.1){
862 // snoise2 = gnoise2;
863 // res = ggg2[amp]*z+glandau2[amp]*angle2;
864 // if (ctype==0) res *= gcor02[rsigmaz];
867 // res*= gcorp[rsigmaz];
871 // snoise2 = gnoise3;
872 // res = ggg3[amp]*z+glandau3[amp]*angle2;
873 // if (ctype==0) res *= gcor02[rsigmaz];
876 // res*= gcorp[rsigmaz];
885 // if ((ctype<0) &&<70){
890 // if (res<2*snoise2)
892 // if (res>3) res =3;
893 // seed->SetErrorZ2(res);
901 void AliTPCtrackerMI::RotateToLocal(AliTPCseed *seed)
903 //rotate to track "local coordinata
904 Float_t x = seed->GetX();
905 Float_t y = seed->GetY();
906 Float_t ymax = x*TMath::Tan(0.5*fSectors->GetAlpha());
909 seed->SetRelativeSector((seed->GetRelativeSector()+1) % fN);
910 if (!seed->Rotate(fSectors->GetAlpha()))
912 } else if (y <-ymax) {
913 seed->SetRelativeSector((seed->GetRelativeSector()-1+fN) % fN);
914 if (!seed->Rotate(-fSectors->GetAlpha()))
922 //_____________________________________________________________________________
923 Double_t AliTPCtrackerMI::F1old(Double_t x1,Double_t y1,
924 Double_t x2,Double_t y2,
925 Double_t x3,Double_t y3)
927 //-----------------------------------------------------------------
928 // Initial approximation of the track curvature
929 //-----------------------------------------------------------------
930 Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
931 Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
932 (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
933 Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
934 (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
936 Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
937 if ( xr*xr+yr*yr<=0.00000000000001) return 100;
938 return -xr*yr/sqrt(xr*xr+yr*yr);
943 //_____________________________________________________________________________
944 Double_t AliTPCtrackerMI::F1(Double_t x1,Double_t y1,
945 Double_t x2,Double_t y2,
946 Double_t x3,Double_t y3)
948 //-----------------------------------------------------------------
949 // Initial approximation of the track curvature
950 //-----------------------------------------------------------------
956 Double_t det = x3*y2-x2*y3;
961 Double_t u = 0.5* (x2*(x2-x3)+y2*(y2-y3))/det;
962 Double_t x0 = x3*0.5-y3*u;
963 Double_t y0 = y3*0.5+x3*u;
964 Double_t c2 = 1/TMath::Sqrt(x0*x0+y0*y0);
970 Double_t AliTPCtrackerMI::F2(Double_t x1,Double_t y1,
971 Double_t x2,Double_t y2,
972 Double_t x3,Double_t y3)
974 //-----------------------------------------------------------------
975 // Initial approximation of the track curvature
976 //-----------------------------------------------------------------
982 Double_t det = x3*y2-x2*y3;
987 Double_t u = 0.5* (x2*(x2-x3)+y2*(y2-y3))/det;
988 Double_t x0 = x3*0.5-y3*u;
989 Double_t y0 = y3*0.5+x3*u;
990 Double_t c2 = 1/TMath::Sqrt(x0*x0+y0*y0);
999 //_____________________________________________________________________________
1000 Double_t AliTPCtrackerMI::F2old(Double_t x1,Double_t y1,
1001 Double_t x2,Double_t y2,
1002 Double_t x3,Double_t y3)
1004 //-----------------------------------------------------------------
1005 // Initial approximation of the track curvature times center of curvature
1006 //-----------------------------------------------------------------
1007 Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
1008 Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
1009 (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
1010 Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
1011 (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
1013 Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
1015 return -a/(d*y1-b)*xr/sqrt(xr*xr+yr*yr);
1018 //_____________________________________________________________________________
1019 Double_t AliTPCtrackerMI::F3(Double_t x1,Double_t y1,
1020 Double_t x2,Double_t y2,
1021 Double_t z1,Double_t z2)
1023 //-----------------------------------------------------------------
1024 // Initial approximation of the tangent of the track dip angle
1025 //-----------------------------------------------------------------
1026 return (z1 - z2)/sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
1030 Double_t AliTPCtrackerMI::F3n(Double_t x1,Double_t y1,
1031 Double_t x2,Double_t y2,
1032 Double_t z1,Double_t z2, Double_t c)
1034 //-----------------------------------------------------------------
1035 // Initial approximation of the tangent of the track dip angle
1036 //-----------------------------------------------------------------
1040 //angle1 = (z1-z2)*c/(TMath::ASin(c*x1-ni)-TMath::ASin(c*x2-ni));
1042 Double_t d = TMath::Sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
1043 if (TMath::Abs(d*c*0.5)>1) return 0;
1044 // Double_t angle2 = TMath::ASin(d*c*0.5);
1045 // Double_t angle2 = AliTPCFastMath::FastAsin(d*c*0.5);
1046 Double_t angle2 = (d*c*0.5>0.1)? TMath::ASin(d*c*0.5): AliTPCFastMath::FastAsin(d*c*0.5);
1048 angle2 = (z1-z2)*c/(angle2*2.);
1052 Bool_t AliTPCtrackerMI::GetProlongation(Double_t x1, Double_t x2, Double_t x[5], Double_t &y, Double_t &z)
1053 {//-----------------------------------------------------------------
1054 // This function find proloncation of a track to a reference plane x=x2.
1055 //-----------------------------------------------------------------
1059 if (TMath::Abs(x[4]*x1 - x[2]) >= 0.999) {
1063 Double_t c1=x[4]*x1 - x[2], r1=sqrt(1.- c1*c1);
1064 Double_t c2=x[4]*x2 - x[2], r2=sqrt(1.- c2*c2);
1068 Double_t dy = dx*(c1+c2)/(r1+r2);
1071 Double_t delta = x[4]*dx*(c1+c2)/(c1*r2 + c2*r1);
1073 if (TMath::Abs(delta)>0.01){
1074 dz = x[3]*TMath::ASin(delta)/x[4];
1076 dz = x[3]*AliTPCFastMath::FastAsin(delta)/x[4];
1079 //dz = x[3]*AliTPCFastMath::FastAsin(delta)/x[4];
1087 Int_t AliTPCtrackerMI::LoadClusters (TTree *tree)
1092 return LoadClusters();
1096 Int_t AliTPCtrackerMI::LoadClusters(TObjArray *arr)
1099 // load clusters to the memory
1100 AliTPCClustersRow *clrow = 0x0;
1101 Int_t lower = arr->LowerBound();
1102 Int_t entries = arr->GetEntriesFast();
1104 for (Int_t i=lower; i<entries; i++) {
1105 clrow = (AliTPCClustersRow*) arr->At(i);
1106 if(!clrow) continue;
1107 if(!clrow->GetArray()) continue;
1111 fParam->AdjustSectorRow(clrow->GetID(),sec,row);
1113 for (Int_t icl=0; icl<clrow->GetArray()->GetEntriesFast(); icl++){
1114 Transform((AliTPCclusterMI*)(clrow->GetArray()->At(icl)));
1117 if (clrow->GetArray()->GetEntriesFast()<=0) continue;
1118 AliTPCtrackerRow * tpcrow=0;
1121 tpcrow = &(fInnerSec[sec%fkNIS][row]);
1125 tpcrow = &(fOuterSec[(sec-fkNIS*2)%fkNOS][row]);
1126 left = (sec-fkNIS*2)/fkNOS;
1129 tpcrow->SetN1(clrow->GetArray()->GetEntriesFast());
1130 tpcrow->SetClusters1(new AliTPCclusterMI[tpcrow->GetN1()]);
1131 for (Int_t j=0;j<tpcrow->GetN1();++j)
1132 tpcrow->SetCluster1(j, *(AliTPCclusterMI*)(clrow->GetArray()->At(j)));
1135 tpcrow->SetN2(clrow->GetArray()->GetEntriesFast());
1136 tpcrow->SetClusters2(new AliTPCclusterMI[tpcrow->GetN2()]);
1137 for (Int_t j=0;j<tpcrow->GetN2();++j)
1138 tpcrow->SetCluster2(j,*(AliTPCclusterMI*)(clrow->GetArray()->At(j)));
1148 Int_t AliTPCtrackerMI::LoadClusters(TClonesArray *arr)
1151 // load clusters to the memory from one
1154 AliTPCclusterMI *clust=0;
1155 Int_t count[72][96] = { {0} , {0} };
1157 // loop over clusters
1158 for (Int_t icl=0; icl<arr->GetEntriesFast(); icl++) {
1159 clust = (AliTPCclusterMI*)arr->At(icl);
1160 if(!clust) continue;
1161 //printf("cluster: det %d, row %d \n", clust->GetDetector(),clust->GetRow());
1163 // transform clusters
1166 // count clusters per pad row
1167 count[clust->GetDetector()][clust->GetRow()]++;
1170 // insert clusters to sectors
1171 for (Int_t icl=0; icl<arr->GetEntriesFast(); icl++) {
1172 clust = (AliTPCclusterMI*)arr->At(icl);
1173 if(!clust) continue;
1175 Int_t sec = clust->GetDetector();
1176 Int_t row = clust->GetRow();
1178 // filter overlapping pad rows needed by HLT
1179 if(sec<fkNIS*2) { //IROCs
1180 if(row == 30) continue;
1183 if(row == 27 || row == 76) continue;
1189 fInnerSec[sec%fkNIS].InsertCluster(clust, count[sec][row], fParam);
1192 left = (sec-fkNIS*2)/fkNOS;
1193 fOuterSec[(sec-fkNIS*2)%fkNOS].InsertCluster(clust, count[sec][row], fParam);
1197 // Load functions must be called behind LoadCluster(TClonesArray*)
1199 //LoadOuterSectors();
1200 //LoadInnerSectors();
1206 Int_t AliTPCtrackerMI::LoadClusters()
1209 // load clusters to the memory
1210 AliTPCClustersRow *clrow= new AliTPCClustersRow;
1211 clrow->SetClass("AliTPCclusterMI");
1213 clrow->GetArray()->ExpandCreateFast(10000);
1215 // TTree * tree = fClustersArray.GetTree();
1217 TTree * tree = fInput;
1218 TBranch * br = tree->GetBranch("Segment");
1219 br->SetAddress(&clrow);
1221 Int_t j=Int_t(tree->GetEntries());
1222 for (Int_t i=0; i<j; i++) {
1226 fParam->AdjustSectorRow(clrow->GetID(),sec,row);
1227 for (Int_t icl=0; icl<clrow->GetArray()->GetEntriesFast(); icl++){
1228 Transform((AliTPCclusterMI*)(clrow->GetArray()->At(icl)));
1231 AliTPCtrackerRow * tpcrow=0;
1234 tpcrow = &(fInnerSec[sec%fkNIS][row]);
1238 tpcrow = &(fOuterSec[(sec-fkNIS*2)%fkNOS][row]);
1239 left = (sec-fkNIS*2)/fkNOS;
1242 tpcrow->SetN1(clrow->GetArray()->GetEntriesFast());
1243 tpcrow->SetClusters1(new AliTPCclusterMI[tpcrow->GetN1()]);
1244 for (Int_t k=0;k<tpcrow->GetN1();++k)
1245 tpcrow->SetCluster1(k, *(AliTPCclusterMI*)(clrow->GetArray()->At(k)));
1248 tpcrow->SetN2(clrow->GetArray()->GetEntriesFast());
1249 tpcrow->SetClusters2(new AliTPCclusterMI[tpcrow->GetN2()]);
1250 for (Int_t k=0;k<tpcrow->GetN2();++k)
1251 tpcrow->SetCluster2(k,*(AliTPCclusterMI*)(clrow->GetArray()->At(k)));
1262 void AliTPCtrackerMI::UnloadClusters()
1265 // unload clusters from the memory
1267 Int_t nrows = fOuterSec->GetNRows();
1268 for (Int_t sec = 0;sec<fkNOS;sec++)
1269 for (Int_t row = 0;row<nrows;row++){
1270 AliTPCtrackerRow* tpcrow = &(fOuterSec[sec%fkNOS][row]);
1272 // if (tpcrow->fClusters1) delete []tpcrow->fClusters1;
1273 // if (tpcrow->fClusters2) delete []tpcrow->fClusters2;
1275 tpcrow->ResetClusters();
1278 nrows = fInnerSec->GetNRows();
1279 for (Int_t sec = 0;sec<fkNIS;sec++)
1280 for (Int_t row = 0;row<nrows;row++){
1281 AliTPCtrackerRow* tpcrow = &(fInnerSec[sec%fkNIS][row]);
1283 // if (tpcrow->fClusters1) delete []tpcrow->fClusters1;
1284 //if (tpcrow->fClusters2) delete []tpcrow->fClusters2;
1286 tpcrow->ResetClusters();
1292 void AliTPCtrackerMI::FillClusterArray(TObjArray* array) const{
1294 // Filling cluster to the array - For visualization purposes
1297 nrows = fOuterSec->GetNRows();
1298 for (Int_t sec = 0;sec<fkNOS;sec++)
1299 for (Int_t row = 0;row<nrows;row++){
1300 AliTPCtrackerRow* tpcrow = &(fOuterSec[sec%fkNOS][row]);
1301 if (!tpcrow) continue;
1302 for (Int_t icl = 0;icl<tpcrow->GetN();icl++){
1303 array->AddLast((TObject*)((*tpcrow)[icl]));
1306 nrows = fInnerSec->GetNRows();
1307 for (Int_t sec = 0;sec<fkNIS;sec++)
1308 for (Int_t row = 0;row<nrows;row++){
1309 AliTPCtrackerRow* tpcrow = &(fInnerSec[sec%fkNIS][row]);
1310 if (!tpcrow) continue;
1311 for (Int_t icl = 0;icl<tpcrow->GetN();icl++){
1312 array->AddLast((TObject*)(*tpcrow)[icl]);
1318 void AliTPCtrackerMI::Transform(AliTPCclusterMI * cluster){
1322 AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ;
1324 AliFatal("Tranformations not in calibDB");
1326 Double_t x[3]={cluster->GetRow(),cluster->GetPad(),cluster->GetTimeBin()};
1327 Int_t i[1]={cluster->GetDetector()};
1328 transform->Transform(x,i,0,1);
1329 // if (cluster->GetDetector()%36>17){
1334 // in debug mode check the transformation
1336 if (AliTPCReconstructor::StreamLevel()>1) {
1338 cluster->GetGlobalXYZ(gx);
1339 Int_t event = (fEvent==NULL)? 0: fEvent->GetEventNumberInFile();
1340 TTreeSRedirector &cstream = *fDebugStreamer;
1341 cstream<<"Transform"<<
1352 cluster->SetX(x[0]);
1353 cluster->SetY(x[1]);
1354 cluster->SetZ(x[2]);
1360 //if (!fParam->IsGeoRead()) fParam->ReadGeoMatrices();
1361 TGeoHMatrix *mat = fParam->GetClusterMatrix(cluster->GetDetector());
1363 Double_t pos[3]= {cluster->GetX(),cluster->GetY(),cluster->GetZ()};
1364 Double_t posC[3]={cluster->GetX(),cluster->GetY(),cluster->GetZ()};
1365 if (mat) mat->LocalToMaster(pos,posC);
1367 // chack Loading of Geo matrices from GeoManager - TEMPORARY FIX
1369 cluster->SetX(posC[0]);
1370 cluster->SetY(posC[1]);
1371 cluster->SetZ(posC[2]);
1374 //_____________________________________________________________________________
1375 Int_t AliTPCtrackerMI::LoadOuterSectors() {
1376 //-----------------------------------------------------------------
1377 // This function fills outer TPC sectors with clusters.
1378 //-----------------------------------------------------------------
1379 Int_t nrows = fOuterSec->GetNRows();
1381 for (Int_t sec = 0;sec<fkNOS;sec++)
1382 for (Int_t row = 0;row<nrows;row++){
1383 AliTPCtrackerRow* tpcrow = &(fOuterSec[sec%fkNOS][row]);
1384 Int_t sec2 = sec+2*fkNIS;
1386 Int_t ncl = tpcrow->GetN1();
1388 AliTPCclusterMI *c= (tpcrow->GetCluster1(ncl));
1389 index=(((sec2<<8)+row)<<16)+ncl;
1390 tpcrow->InsertCluster(c,index);
1393 ncl = tpcrow->GetN2();
1395 AliTPCclusterMI *c= (tpcrow->GetCluster2(ncl));
1396 index=((((sec2+fkNOS)<<8)+row)<<16)+ncl;
1397 tpcrow->InsertCluster(c,index);
1400 // write indexes for fast acces
1402 for (Int_t i=0;i<510;i++)
1403 tpcrow->SetFastCluster(i,-1);
1404 for (Int_t i=0;i<tpcrow->GetN();i++){
1405 Int_t zi = Int_t((*tpcrow)[i]->GetZ()+255.);
1406 tpcrow->SetFastCluster(zi,i); // write index
1409 for (Int_t i=0;i<510;i++){
1410 if (tpcrow->GetFastCluster(i)<0)
1411 tpcrow->SetFastCluster(i,last);
1413 last = tpcrow->GetFastCluster(i);
1422 //_____________________________________________________________________________
1423 Int_t AliTPCtrackerMI::LoadInnerSectors() {
1424 //-----------------------------------------------------------------
1425 // This function fills inner TPC sectors with clusters.
1426 //-----------------------------------------------------------------
1427 Int_t nrows = fInnerSec->GetNRows();
1429 for (Int_t sec = 0;sec<fkNIS;sec++)
1430 for (Int_t row = 0;row<nrows;row++){
1431 AliTPCtrackerRow* tpcrow = &(fInnerSec[sec%fkNIS][row]);
1434 Int_t ncl = tpcrow->GetN1();
1436 AliTPCclusterMI *c= (tpcrow->GetCluster1(ncl));
1437 index=(((sec<<8)+row)<<16)+ncl;
1438 tpcrow->InsertCluster(c,index);
1441 ncl = tpcrow->GetN2();
1443 AliTPCclusterMI *c= (tpcrow->GetCluster2(ncl));
1444 index=((((sec+fkNIS)<<8)+row)<<16)+ncl;
1445 tpcrow->InsertCluster(c,index);
1448 // write indexes for fast acces
1450 for (Int_t i=0;i<510;i++)
1451 tpcrow->SetFastCluster(i,-1);
1452 for (Int_t i=0;i<tpcrow->GetN();i++){
1453 Int_t zi = Int_t((*tpcrow)[i]->GetZ()+255.);
1454 tpcrow->SetFastCluster(zi,i); // write index
1457 for (Int_t i=0;i<510;i++){
1458 if (tpcrow->GetFastCluster(i)<0)
1459 tpcrow->SetFastCluster(i,last);
1461 last = tpcrow->GetFastCluster(i);
1473 //_________________________________________________________________________
1474 AliTPCclusterMI *AliTPCtrackerMI::GetClusterMI(Int_t index) const {
1475 //--------------------------------------------------------------------
1476 // Return pointer to a given cluster
1477 //--------------------------------------------------------------------
1478 if (index<0) return 0; // no cluster
1479 Int_t sec=(index&0xff000000)>>24;
1480 Int_t row=(index&0x00ff0000)>>16;
1481 Int_t ncl=(index&0x00007fff)>>00;
1483 const AliTPCtrackerRow * tpcrow=0;
1484 AliTPCclusterMI * clrow =0;
1486 if (sec<0 || sec>=fkNIS*4) {
1487 AliWarning(Form("Wrong sector %d",sec));
1492 AliTPCtrackerSector& tracksec = fInnerSec[sec%fkNIS];
1493 if (tracksec.GetNRows()<=row) return 0;
1494 tpcrow = &(tracksec[row]);
1495 if (tpcrow==0) return 0;
1498 if (tpcrow->GetN1()<=ncl) return 0;
1499 clrow = tpcrow->GetClusters1();
1502 if (tpcrow->GetN2()<=ncl) return 0;
1503 clrow = tpcrow->GetClusters2();
1507 AliTPCtrackerSector& tracksec = fOuterSec[(sec-fkNIS*2)%fkNOS];
1508 if (tracksec.GetNRows()<=row) return 0;
1509 tpcrow = &(tracksec[row]);
1510 if (tpcrow==0) return 0;
1512 if (sec-2*fkNIS<fkNOS) {
1513 if (tpcrow->GetN1()<=ncl) return 0;
1514 clrow = tpcrow->GetClusters1();
1517 if (tpcrow->GetN2()<=ncl) return 0;
1518 clrow = tpcrow->GetClusters2();
1522 return &(clrow[ncl]);
1528 Int_t AliTPCtrackerMI::FollowToNext(AliTPCseed& t, Int_t nr) {
1529 //-----------------------------------------------------------------
1530 // This function tries to find a track prolongation to next pad row
1531 //-----------------------------------------------------------------
1533 Double_t x= GetXrow(nr), ymax=GetMaxY(nr);
1534 AliTPCclusterMI *cl=0;
1535 Int_t tpcindex= t.GetClusterIndex2(nr);
1537 // update current shape info every 5 pad-row
1538 // if ( (nr%5==0) || t.GetNumberOfClusters()<2 || (t.fCurrentSigmaY2<0.0001) ){
1542 if (fIteration>0 && tpcindex>=-1){ //if we have already clusters
1544 if (tpcindex==-1) return 0; //track in dead zone
1546 cl = t.GetClusterPointer(nr);
1547 if ( (cl==0) ) cl = GetClusterMI(tpcindex);
1548 t.SetCurrentClusterIndex1(tpcindex);
1551 Int_t relativesector = ((tpcindex&0xff000000)>>24)%18; // if previously accepted cluster in different sector
1552 Float_t angle = relativesector*fSectors->GetAlpha()+fSectors->GetAlphaShift();
1554 if (angle<-TMath::Pi()) angle += 2*TMath::Pi();
1555 if (angle>=TMath::Pi()) angle -= 2*TMath::Pi();
1557 if (TMath::Abs(angle-t.GetAlpha())>0.001){
1558 Double_t rotation = angle-t.GetAlpha();
1559 t.SetRelativeSector(relativesector);
1560 if (!t.Rotate(rotation)) return 0;
1562 if (!t.PropagateTo(x)) return 0;
1564 t.SetCurrentCluster(cl);
1566 Int_t accept = AcceptCluster(&t,t.GetCurrentCluster());
1567 if ((tpcindex&0x8000)==0) accept =0;
1569 //if founded cluster is acceptible
1570 if (cl->IsUsed(11)) { // id cluster is shared inrease uncertainty
1571 t.SetErrorY2(t.GetErrorY2()+0.03);
1572 t.SetErrorZ2(t.GetErrorZ2()+0.03);
1573 t.SetErrorY2(t.GetErrorY2()*3);
1574 t.SetErrorZ2(t.GetErrorZ2()*3);
1576 t.SetNFoundable(t.GetNFoundable()+1);
1577 UpdateTrack(&t,accept);
1582 if (TMath::Abs(t.GetSnp())>AliTPCReconstructor::GetMaxSnpTracker()) return 0; // cut on angle
1584 // not look for new cluster during refitting
1585 t.SetNFoundable(t.GetNFoundable()+1);
1590 // if (TMath::Abs(t.GetSnp())>0.95 || TMath::Abs(x*t.GetC()-t.GetEta())>0.95) return 0;// patch 28 fev 06
1591 Double_t y=t.GetYat(x);
1592 if (TMath::Abs(y)>ymax){
1594 t.SetRelativeSector((t.GetRelativeSector()+1) % fN);
1595 if (!t.Rotate(fSectors->GetAlpha()))
1597 } else if (y <-ymax) {
1598 t.SetRelativeSector((t.GetRelativeSector()-1+fN) % fN);
1599 if (!t.Rotate(-fSectors->GetAlpha()))
1605 if (!t.PropagateTo(x)) {
1606 if (fIteration==0) t.SetRemoval(10);
1610 Double_t z=t.GetZ();
1613 if (!IsActive(t.GetRelativeSector(),nr)) {
1615 t.SetClusterIndex2(nr,-1);
1618 //AliInfo(Form("A - Sector%d phi %f - alpha %f", t.fRelativeSector,y/x, t.GetAlpha()));
1619 Bool_t isActive = IsActive(t.GetRelativeSector(),nr);
1620 Bool_t isActive2 = (nr>=fInnerSec->GetNRows()) ? fOuterSec[t.GetRelativeSector()][nr-fInnerSec->GetNRows()].GetN()>0:fInnerSec[t.GetRelativeSector()][nr].GetN()>0;
1622 if (!isActive || !isActive2) return 0;
1624 const AliTPCtrackerRow &krow=GetRow(t.GetRelativeSector(),nr);
1625 if ( (t.GetSigmaY2()<0) || t.GetSigmaZ2()<0) return 0;
1627 Double_t roadz = 1.;
1629 if (TMath::Abs(TMath::Abs(y)-ymax)<krow.GetDeadZone()){
1631 t.SetClusterIndex2(nr,-1);
1636 if (TMath::Abs(z)<(AliTPCReconstructor::GetCtgRange()*x+10) && TMath::Abs(z)<fParam->GetZLength(0) && (TMath::Abs(t.GetSnp())<AliTPCReconstructor::GetMaxSnpTracker()))
1637 t.SetNFoundable(t.GetNFoundable()+1);
1643 // cl = krow.FindNearest2(y+10.,z,roady,roadz,index);
1644 cl = krow.FindNearest2(y,z,roady,roadz,index);
1645 if (cl) t.SetCurrentClusterIndex1(krow.GetIndex(index));
1648 t.SetCurrentCluster(cl);
1650 if (fIteration==2&&cl->IsUsed(10)) return 0;
1651 Int_t accept = AcceptCluster(&t,t.GetCurrentCluster());
1652 if (fIteration==2&&cl->IsUsed(11)) {
1653 t.SetErrorY2(t.GetErrorY2()+0.03);
1654 t.SetErrorZ2(t.GetErrorZ2()+0.03);
1655 t.SetErrorY2(t.GetErrorY2()*3);
1656 t.SetErrorZ2(t.GetErrorZ2()*3);
1659 if (t.fCurrentCluster->IsUsed(10)){
1664 if (t.fNShared>0.7*t.GetNumberOfClusters()) {
1670 if (accept<3) UpdateTrack(&t,accept);
1673 if ( fIteration==0 && t.GetNFoundable()*0.5 > t.GetNumberOfClusters()) t.SetRemoval(10);
1679 Int_t AliTPCtrackerMI::FollowToNextFast(AliTPCseed& t, Int_t nr) {
1680 //-----------------------------------------------------------------
1681 // This function tries to find a track prolongation to next pad row
1682 //-----------------------------------------------------------------
1684 Double_t x= GetXrow(nr), ymax=GetMaxY(nr);
1686 if (!t.GetProlongation(x,y,z)) {
1692 if (TMath::Abs(y)>ymax){
1695 t.SetRelativeSector((t.GetRelativeSector()+1) % fN);
1696 if (!t.Rotate(fSectors->GetAlpha()))
1698 } else if (y <-ymax) {
1699 t.SetRelativeSector((t.GetRelativeSector()-1+fN) % fN);
1700 if (!t.Rotate(-fSectors->GetAlpha()))
1703 if (!t.PropagateTo(x)) {
1706 t.GetProlongation(x,y,z);
1709 // update current shape info every 2 pad-row
1710 if ( (nr%2==0) || t.GetNumberOfClusters()<2 || (t.GetCurrentSigmaY2()<0.0001) ){
1711 // t.fCurrentSigmaY = GetSigmaY(&t);
1712 //t.fCurrentSigmaZ = GetSigmaZ(&t);
1716 AliTPCclusterMI *cl=0;
1721 const AliTPCtrackerRow &krow=GetRow(t.GetRelativeSector(),nr);
1722 if ( (t.GetSigmaY2()<0) || t.GetSigmaZ2()<0) return 0;
1724 Double_t roadz = 1.;
1727 if (TMath::Abs(TMath::Abs(y)-ymax)<krow.GetDeadZone()){
1729 t.SetClusterIndex2(row,-1);
1734 if (TMath::Abs(z)>(AliTPCReconstructor::GetCtgRange()*x+10)) t.SetClusterIndex2(row,-1);
1738 if ((cl==0)&&(krow)) {
1739 // cl = krow.FindNearest2(y+10,z,roady,roadz,index);
1740 cl = krow.FindNearest2(y,z,roady,roadz,index);
1742 if (cl) t.SetCurrentClusterIndex1(krow.GetIndex(index));
1746 t.SetCurrentCluster(cl);
1747 // Int_t accept = AcceptCluster(&t,t.fCurrentCluster);
1749 t.SetClusterIndex2(row,index);
1750 t.SetClusterPointer(row, cl);
1758 //_________________________________________________________________________
1759 Bool_t AliTPCtrackerMI::GetTrackPoint(Int_t index, AliTrackPoint &p ) const
1761 // Get track space point by index
1762 // return false in case the cluster doesn't exist
1763 AliTPCclusterMI *cl = GetClusterMI(index);
1764 if (!cl) return kFALSE;
1765 Int_t sector = (index&0xff000000)>>24;
1766 // Int_t row = (index&0x00ff0000)>>16;
1768 // xyz[0] = fParam->GetPadRowRadii(sector,row);
1769 xyz[0] = cl->GetX();
1770 xyz[1] = cl->GetY();
1771 xyz[2] = cl->GetZ();
1773 fParam->AdjustCosSin(sector,cos,sin);
1774 Float_t x = cos*xyz[0]-sin*xyz[1];
1775 Float_t y = cos*xyz[1]+sin*xyz[0];
1777 Float_t sigmaY2 = 0.027*cl->GetSigmaY2();
1778 if (sector < fParam->GetNInnerSector()) sigmaY2 *= 2.07;
1779 Float_t sigmaZ2 = 0.066*cl->GetSigmaZ2();
1780 if (sector < fParam->GetNInnerSector()) sigmaZ2 *= 1.77;
1781 cov[0] = sin*sin*sigmaY2;
1782 cov[1] = -sin*cos*sigmaY2;
1784 cov[3] = cos*cos*sigmaY2;
1787 p.SetXYZ(x,y,xyz[2],cov);
1788 AliGeomManager::ELayerID iLayer;
1790 if (sector < fParam->GetNInnerSector()) {
1791 iLayer = AliGeomManager::kTPC1;
1795 iLayer = AliGeomManager::kTPC2;
1796 idet = sector - fParam->GetNInnerSector();
1798 UShort_t volid = AliGeomManager::LayerToVolUID(iLayer,idet);
1799 p.SetVolumeID(volid);
1805 Int_t AliTPCtrackerMI::UpdateClusters(AliTPCseed& t, Int_t nr) {
1806 //-----------------------------------------------------------------
1807 // This function tries to find a track prolongation to next pad row
1808 //-----------------------------------------------------------------
1809 t.SetCurrentCluster(0);
1810 t.SetCurrentClusterIndex1(0);
1812 Double_t xt=t.GetX();
1813 Int_t row = GetRowNumber(xt)-1;
1814 Double_t ymax= GetMaxY(nr);
1816 if (row < nr) return 1; // don't prolongate if not information until now -
1817 // if (TMath::Abs(t.GetSnp())>0.9 && t.GetNumberOfClusters()>40. && fIteration!=2) {
1819 // return 0; // not prolongate strongly inclined tracks
1821 // if (TMath::Abs(t.GetSnp())>0.95) {
1823 // return 0; // not prolongate strongly inclined tracks
1824 // }// patch 28 fev 06
1826 Double_t x= GetXrow(nr);
1828 //t.PropagateTo(x+0.02);
1829 //t.PropagateTo(x+0.01);
1830 if (!t.PropagateTo(x)){
1837 if (TMath::Abs(y)>ymax){
1839 t.SetRelativeSector((t.GetRelativeSector()+1) % fN);
1840 if (!t.Rotate(fSectors->GetAlpha()))
1842 } else if (y <-ymax) {
1843 t.SetRelativeSector((t.GetRelativeSector()-1+fN) % fN);
1844 if (!t.Rotate(-fSectors->GetAlpha()))
1847 // if (!t.PropagateTo(x)){
1854 if (TMath::Abs(t.GetSnp())>AliTPCReconstructor::GetMaxSnpTracker()) return 0;
1856 if (!IsActive(t.GetRelativeSector(),nr)) {
1858 t.SetClusterIndex2(nr,-1);
1861 //AliInfo(Form("A - Sector%d phi %f - alpha %f", t.fRelativeSector,y/x, t.GetAlpha()));
1863 AliTPCtrackerRow &krow=GetRow(t.GetRelativeSector(),nr);
1865 if (TMath::Abs(TMath::Abs(y)-ymax)<krow.GetDeadZone()){
1867 t.SetClusterIndex2(nr,-1);
1872 if (TMath::Abs(t.GetZ())<(AliTPCReconstructor::GetCtgRange()*t.GetX()+10) && (TMath::Abs(t.GetSnp())<AliTPCReconstructor::GetMaxSnpTracker())) t.SetNFoundable(t.GetNFoundable()+1);
1878 if ( (nr%2==0) || t.GetNumberOfClusters()<2){
1879 // t.fCurrentSigmaY = GetSigmaY(&t);
1880 //t.fCurrentSigmaZ = GetSigmaZ(&t);
1884 AliTPCclusterMI *cl=0;
1887 Double_t roady = 1.;
1888 Double_t roadz = 1.;
1892 index = t.GetClusterIndex2(nr);
1893 if ( (index>0) && (index&0x8000)==0){
1894 cl = t.GetClusterPointer(nr);
1895 if ( (cl==0) && (index>0)) cl = GetClusterMI(index);
1896 t.SetCurrentClusterIndex1(index);
1898 t.SetCurrentCluster(cl);
1904 // if (index<0) return 0;
1905 UInt_t uindex = TMath::Abs(index);
1908 //cl = krow.FindNearest2(y+10,z,roady,roadz,uindex);
1909 cl = krow.FindNearest2(y,z,roady,roadz,uindex);
1912 if (cl) t.SetCurrentClusterIndex1(krow.GetIndex(uindex));
1913 t.SetCurrentCluster(cl);
1919 Int_t AliTPCtrackerMI::FollowToNextCluster(AliTPCseed & t, Int_t nr) {
1920 //-----------------------------------------------------------------
1921 // This function tries to find a track prolongation to next pad row
1922 //-----------------------------------------------------------------
1924 //update error according neighborhoud
1926 if (t.GetCurrentCluster()) {
1928 Int_t accept = AcceptCluster(&t,t.GetCurrentCluster());
1930 if (t.GetCurrentCluster()->IsUsed(10)){
1935 t.SetNShared(t.GetNShared()+1);
1936 if (t.GetNShared()>0.7*t.GetNumberOfClusters()) {
1941 if (fIteration>0) accept = 0;
1942 if (accept<3) UpdateTrack(&t,accept);
1946 if ( ( (t.GetSigmaY2()+t.GetSigmaZ2())>0.16)&& t.GetNumberOfClusters()>18) t.SetRemoval(10);
1947 if ( t.GetChi2()/t.GetNumberOfClusters()>6 &&t.GetNumberOfClusters()>18) t.SetRemoval(10);
1949 if (( (t.GetNFoundable()*0.5 > t.GetNumberOfClusters()) || t.GetNoCluster()>15)) t.SetRemoval(10);
1957 //_____________________________________________________________________________
1958 Int_t AliTPCtrackerMI::FollowProlongation(AliTPCseed& t, Int_t rf, Int_t step) {
1959 //-----------------------------------------------------------------
1960 // This function tries to find a track prolongation.
1961 //-----------------------------------------------------------------
1962 Double_t xt=t.GetX();
1964 Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
1965 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
1966 if (alpha < 0. ) alpha += 2.*TMath::Pi();
1968 t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN);
1970 Int_t first = GetRowNumber(xt)-1;
1971 for (Int_t nr= first; nr>=rf; nr-=step) {
1973 if (t.GetKinkIndexes()[0]>0){
1974 for (Int_t i=0;i<3;i++){
1975 Int_t index = t.GetKinkIndexes()[i];
1976 if (index==0) break;
1977 if (index<0) continue;
1979 AliKink * kink = (AliKink*)fEvent->GetKink(index-1);
1981 printf("PROBLEM\n");
1984 Int_t kinkrow = kink->GetTPCRow0()+2+Int_t(0.5/(0.05+kink->GetAngle(2)));
1986 AliExternalTrackParam paramd(t);
1987 kink->SetDaughter(paramd);
1988 kink->SetStatus(2,5);
1995 if (nr==80) t.UpdateReference();
1996 if (nr<fInnerSec->GetNRows())
1997 fSectors = fInnerSec;
1999 fSectors = fOuterSec;
2000 if (FollowToNext(t,nr)==0)
2009 //_____________________________________________________________________________
2010 Int_t AliTPCtrackerMI::FollowProlongationFast(AliTPCseed& t, Int_t rf, Int_t step) {
2011 //-----------------------------------------------------------------
2012 // This function tries to find a track prolongation.
2013 //-----------------------------------------------------------------
2014 Double_t xt=t.GetX();
2016 Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
2017 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
2018 if (alpha < 0. ) alpha += 2.*TMath::Pi();
2019 t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN);
2021 for (Int_t nr=GetRowNumber(xt)-1; nr>=rf; nr-=step) {
2023 if (FollowToNextFast(t,nr)==0)
2024 if (!t.IsActive()) return 0;
2034 Int_t AliTPCtrackerMI::FollowBackProlongation(AliTPCseed& t, Int_t rf) {
2035 //-----------------------------------------------------------------
2036 // This function tries to find a track prolongation.
2037 //-----------------------------------------------------------------
2039 Double_t xt=t.GetX();
2040 Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
2041 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
2042 if (alpha < 0. ) alpha += 2.*TMath::Pi();
2043 t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN);
2045 Int_t first = t.GetFirstPoint();
2046 if (first<GetRowNumber(xt)+1) first = GetRowNumber(xt)+1;
2048 if (first<0) first=0;
2049 for (Int_t nr=first; nr<=rf; nr++) {
2050 // if ( (TMath::Abs(t.GetSnp())>0.95)) break;//patch 28 fev 06
2051 if (t.GetKinkIndexes()[0]<0){
2052 for (Int_t i=0;i<3;i++){
2053 Int_t index = t.GetKinkIndexes()[i];
2054 if (index==0) break;
2055 if (index>0) continue;
2056 index = TMath::Abs(index);
2057 AliKink * kink = (AliKink*)fEvent->GetKink(index-1);
2059 printf("PROBLEM\n");
2062 Int_t kinkrow = kink->GetTPCRow0()-2-Int_t(0.5/(0.05+kink->GetAngle(2)));
2064 AliExternalTrackParam paramm(t);
2065 kink->SetMother(paramm);
2066 kink->SetStatus(2,1);
2073 if (nr<fInnerSec->GetNRows())
2074 fSectors = fInnerSec;
2076 fSectors = fOuterSec;
2087 Float_t AliTPCtrackerMI::OverlapFactor(AliTPCseed * s1, AliTPCseed * s2, Int_t &sum1, Int_t & sum2)
2095 Float_t dz2 =(s1->GetZ() - s2->GetZ());
2098 Float_t dy2 =TMath::Abs((s1->GetY() - s2->GetY()));
2100 Float_t distance = TMath::Sqrt(dz2+dy2);
2101 if (distance>4.) return 0; // if there are far away - not overlap - to reduce combinatorics
2104 Int_t firstpoint = TMath::Min(s1->GetFirstPoint(),s2->GetFirstPoint());
2105 Int_t lastpoint = TMath::Max(s1->GetLastPoint(),s2->GetLastPoint());
2110 if (firstpoint>lastpoint) {
2111 firstpoint =lastpoint;
2116 for (Int_t i=firstpoint-1;i<lastpoint+1;i++){
2117 if (s1->GetClusterIndex2(i)>0) sum1++;
2118 if (s2->GetClusterIndex2(i)>0) sum2++;
2119 if (s1->GetClusterIndex2(i)==s2->GetClusterIndex2(i) && s1->GetClusterIndex2(i)>0) {
2123 if (sum<5) return 0;
2125 Float_t summin = TMath::Min(sum1+1,sum2+1);
2126 Float_t ratio = (sum+1)/Float_t(summin);
2130 void AliTPCtrackerMI::SignShared(AliTPCseed * s1, AliTPCseed * s2)
2134 Float_t thetaCut = 0.2;//+10.*TMath::Sqrt(s1->GetSigmaTglZ()+ s2->GetSigmaTglZ());
2135 if (TMath::Abs(s1->GetTgl()-s2->GetTgl())>thetaCut) return;
2136 Float_t minCl = TMath::Min(s1->GetNumberOfClusters(),s2->GetNumberOfClusters());
2137 Int_t cutN0 = TMath::Max(5,TMath::Nint(0.1*minCl));
2142 //Int_t firstpoint = TMath::Max(s1->GetFirstPoint(),s2->GetFirstPoint());
2143 //Int_t lastpoint = TMath::Min(s1->GetLastPoint(),s2->GetLastPoint());
2144 Int_t firstpoint = 0;
2145 Int_t lastpoint = 160;
2147 // if (firstpoint>=lastpoint-5) return;;
2149 for (Int_t i=firstpoint;i<lastpoint;i++){
2150 // if ( (s1->GetClusterIndex2(i)&0xFFFF8FFF)==(s2->GetClusterIndex2(i)&0xFFFF8FFF) && s1->GetClusterIndex2(i)>0) {
2151 if ( (s1->GetClusterIndex2(i))==(s2->GetClusterIndex2(i)) && s1->GetClusterIndex2(i)>0) {
2153 s1->SetSharedMapBit(i, kTRUE);
2154 s2->SetSharedMapBit(i, kTRUE);
2156 if (s1->GetClusterIndex2(i)>0)
2157 s1->SetClusterMapBit(i, kTRUE);
2159 if (sumshared>cutN0){
2162 for (Int_t i=firstpoint;i<lastpoint;i++){
2163 // if ( (s1->GetClusterIndex2(i)&0xFFFF8FFF)==(s2->GetClusterIndex2(i)&0xFFFF8FFF) && s1->GetClusterIndex2(i)>0) {
2164 if ( (s1->GetClusterIndex2(i))==(s2->GetClusterIndex2(i)) && s1->GetClusterIndex2(i)>0) {
2165 AliTPCTrackerPoint *p1 = s1->GetTrackPoint(i);
2166 AliTPCTrackerPoint *p2 = s2->GetTrackPoint(i);;
2167 if (s1->IsActive()&&s2->IsActive()){
2168 p1->SetShared(kTRUE);
2169 p2->SetShared(kTRUE);
2175 if (sumshared>cutN0){
2176 for (Int_t i=0;i<4;i++){
2177 if (s1->GetOverlapLabel(3*i)==0){
2178 s1->SetOverlapLabel(3*i, s2->GetLabel());
2179 s1->SetOverlapLabel(3*i+1,sumshared);
2180 s1->SetOverlapLabel(3*i+2,s2->GetUniqueID());
2184 for (Int_t i=0;i<4;i++){
2185 if (s2->GetOverlapLabel(3*i)==0){
2186 s2->SetOverlapLabel(3*i, s1->GetLabel());
2187 s2->SetOverlapLabel(3*i+1,sumshared);
2188 s2->SetOverlapLabel(3*i+2,s1->GetUniqueID());
2195 void AliTPCtrackerMI::SignShared(TObjArray * arr)
2198 //sort trackss according sectors
2200 for (Int_t i=0; i<arr->GetEntriesFast(); i++) {
2201 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2203 //if (pt) RotateToLocal(pt);
2207 arr->Sort(); // sorting according relative sectors
2208 arr->Expand(arr->GetEntries());
2211 Int_t nseed=arr->GetEntriesFast();
2212 for (Int_t i=0; i<nseed; i++) {
2213 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2215 for (Int_t j=0;j<=12;j++){
2216 pt->SetOverlapLabel(j,0);
2219 for (Int_t i=0; i<nseed; i++) {
2220 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2222 if (pt->GetRemoval()>10) continue;
2223 for (Int_t j=i+1; j<nseed; j++){
2224 AliTPCseed *pt2=(AliTPCseed*)arr->UncheckedAt(j);
2225 if (TMath::Abs(pt->GetRelativeSector()-pt2->GetRelativeSector())>1) continue;
2227 if (pt2->GetRemoval()<=10) {
2228 //if ( TMath::Abs(pt->GetRelativeSector()-pt2->GetRelativeSector())>0) break;
2236 void AliTPCtrackerMI::SortTracks(TObjArray * arr, Int_t mode) const
2239 //sort tracks in array according mode criteria
2240 Int_t nseed = arr->GetEntriesFast();
2241 for (Int_t i=0; i<nseed; i++) {
2242 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2253 void AliTPCtrackerMI::RemoveUsed2(TObjArray * arr, Float_t factor1, Float_t factor2, Int_t minimal)
2256 // Loop over all tracks and remove overlaped tracks (with lower quality)
2258 // 1. Unsign clusters
2259 // 2. Sort tracks according quality
2260 // Quality is defined by the number of cluster between first and last points
2262 // 3. Loop over tracks - decreasing quality order
2263 // a.) remove - If the fraction of shared cluster less than factor (1- n or 2)
2264 // b.) remove - If the minimal number of clusters less than minimal and not ITS
2265 // c.) if track accepted - sign clusters
2267 //Called in - AliTPCtrackerMI::Clusters2Tracks()
2268 // - AliTPCtrackerMI::PropagateBack()
2269 // - AliTPCtrackerMI::RefitInward()
2272 // factor1 - factor for constrained
2273 // factor2 - for non constrained tracks
2274 // if (Float_t(shared+1)/Float_t(found+1)>factor) - DELETE
2278 Int_t nseed = arr->GetEntriesFast();
2279 Float_t * quality = new Float_t[nseed];
2280 Int_t * indexes = new Int_t[nseed];
2284 for (Int_t i=0; i<nseed; i++) {
2285 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2290 pt->UpdatePoints(); //select first last max dens points
2291 Float_t * points = pt->GetPoints();
2292 if (points[3]<0.8) quality[i] =-1;
2293 quality[i] = (points[2]-points[0])+pt->GetNumberOfClusters();
2294 //prefer high momenta tracks if overlaps
2295 quality[i] *= TMath::Sqrt(TMath::Abs(pt->Pt())+0.5);
2297 TMath::Sort(nseed,quality,indexes);
2300 for (Int_t itrack=0; itrack<nseed; itrack++) {
2301 Int_t trackindex = indexes[itrack];
2302 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(trackindex);
2305 if (quality[trackindex]<0){
2307 delete arr->RemoveAt(trackindex);
2310 arr->RemoveAt(trackindex);
2316 Int_t first = Int_t(pt->GetPoints()[0]);
2317 Int_t last = Int_t(pt->GetPoints()[2]);
2318 Double_t factor = (pt->GetBConstrain()) ? factor1: factor2;
2320 Int_t found,foundable,shared;
2321 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
2322 // pt->GetClusterStatistic(0,160, found, foundable,shared,kFALSE);
2323 Bool_t itsgold =kFALSE;
2326 if (pt->GetESD()->GetITSclusters(dummy)>4) itsgold= kTRUE;
2330 if (Float_t(shared+1)/Float_t(found+1)>factor){
2331 if (pt->GetKinkIndexes()[0]!=0) continue; //don't remove tracks - part of the kinks
2332 if( AliTPCReconstructor::StreamLevel()>15){
2333 TTreeSRedirector &cstream = *fDebugStreamer;
2334 cstream<<"RemoveUsed"<<
2335 "iter="<<fIteration<<
2339 delete arr->RemoveAt(trackindex);
2342 if (pt->GetNumberOfClusters()<50&&(found-0.5*shared)<minimal){ //remove short tracks
2343 if (pt->GetKinkIndexes()[0]!=0) continue; //don't remove tracks - part of the kinks
2344 if( AliTPCReconstructor::StreamLevel()>15){
2345 TTreeSRedirector &cstream = *fDebugStreamer;
2346 cstream<<"RemoveShort"<<
2347 "iter="<<fIteration<<
2351 delete arr->RemoveAt(trackindex);
2357 //if (sharedfactor>0.4) continue;
2358 if (pt->GetKinkIndexes()[0]>0) continue;
2359 //Remove tracks with undefined properties - seems
2360 if (pt->GetSigmaY2()<kAlmost0) continue; // ? what is the origin ?
2362 for (Int_t i=first; i<last; i++) {
2363 Int_t index=pt->GetClusterIndex2(i);
2364 // if (index<0 || index&0x8000 ) continue;
2365 if (index<0 || index&0x8000 ) continue;
2366 AliTPCclusterMI *c= pt->GetClusterPointer(i);
2373 Info("RemoveUsed","\n*****\nNumber of good tracks after shared removal\t%d\n",fNtracks);
2379 void AliTPCtrackerMI::UnsignClusters()
2382 // loop over all clusters and unsign them
2385 for (Int_t sec=0;sec<fkNIS;sec++){
2386 for (Int_t row=0;row<fInnerSec->GetNRows();row++){
2387 AliTPCclusterMI *cl = fInnerSec[sec][row].GetClusters1();
2388 for (Int_t icl =0;icl< fInnerSec[sec][row].GetN1();icl++)
2389 // if (cl[icl].IsUsed(10))
2391 cl = fInnerSec[sec][row].GetClusters2();
2392 for (Int_t icl =0;icl< fInnerSec[sec][row].GetN2();icl++)
2393 //if (cl[icl].IsUsed(10))
2398 for (Int_t sec=0;sec<fkNOS;sec++){
2399 for (Int_t row=0;row<fOuterSec->GetNRows();row++){
2400 AliTPCclusterMI *cl = fOuterSec[sec][row].GetClusters1();
2401 for (Int_t icl =0;icl< fOuterSec[sec][row].GetN1();icl++)
2402 //if (cl[icl].IsUsed(10))
2404 cl = fOuterSec[sec][row].GetClusters2();
2405 for (Int_t icl =0;icl< fOuterSec[sec][row].GetN2();icl++)
2406 //if (cl[icl].IsUsed(10))
2415 void AliTPCtrackerMI::SignClusters(TObjArray * arr, Float_t fnumber, Float_t fdensity)
2418 //sign clusters to be "used"
2420 // snumber and sdensity sign number of sigmas - bellow mean value to be accepted
2421 // loop over "primaries"
2435 Int_t nseed = arr->GetEntriesFast();
2436 for (Int_t i=0; i<nseed; i++) {
2437 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2441 if (!(pt->IsActive())) continue;
2442 Float_t dens = pt->GetNumberOfClusters()/Float_t(pt->GetNFoundable());
2443 if ( (dens>0.7) && (pt->GetNumberOfClusters()>70)){
2445 sumdens2+= dens*dens;
2446 sumn += pt->GetNumberOfClusters();
2447 sumn2 += pt->GetNumberOfClusters()*pt->GetNumberOfClusters();
2448 Float_t chi2 = pt->GetChi2()/pt->GetNumberOfClusters();
2451 sumchi2 +=chi2*chi2;
2456 Float_t mdensity = 0.9;
2457 Float_t meann = 130;
2458 Float_t meanchi = 1;
2459 Float_t sdensity = 0.1;
2460 Float_t smeann = 10;
2461 Float_t smeanchi =0.4;
2465 mdensity = sumdens/sum;
2467 meanchi = sumchi/sum;
2469 sdensity = sumdens2/sum-mdensity*mdensity;
2471 sdensity = TMath::Sqrt(sdensity);
2475 smeann = sumn2/sum-meann*meann;
2477 smeann = TMath::Sqrt(smeann);
2481 smeanchi = sumchi2/sum - meanchi*meanchi;
2483 smeanchi = TMath::Sqrt(smeanchi);
2489 //REMOVE SHORT DELTAS or tracks going out of sensitive volume of TPC
2491 for (Int_t i=0; i<nseed; i++) {
2492 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2496 if (pt->GetBSigned()) continue;
2497 if (pt->GetBConstrain()) continue;
2498 //if (!(pt->IsActive())) continue;
2500 Int_t found,foundable,shared;
2501 pt->GetClusterStatistic(0,160,found, foundable,shared);
2502 if (shared/float(found)>0.3) {
2503 if (shared/float(found)>0.9 ){
2504 //delete arr->RemoveAt(i);
2509 Bool_t isok =kFALSE;
2510 if ( (pt->GetNShared()/pt->GetNumberOfClusters()<0.5) &&pt->GetNumberOfClusters()>60)
2512 if ((TMath::Abs(1/pt->GetC())<100.) && (pt->GetNShared()/pt->GetNumberOfClusters()<0.7))
2514 if (TMath::Abs(pt->GetZ()/pt->GetX())>1.1)
2516 if ( (TMath::Abs(pt->GetSnp()>0.7) && pt->GetD(0,0)>60.))
2520 for (Int_t j=0; j<160; ++j) {
2521 Int_t index=pt->GetClusterIndex2(j);
2522 if (index<0) continue;
2523 AliTPCclusterMI *c= pt->GetClusterPointer(j);
2525 //if (!(c->IsUsed(10))) c->Use();
2532 Double_t maxchi = meanchi+2.*smeanchi;
2534 for (Int_t i=0; i<nseed; i++) {
2535 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2539 //if (!(pt->IsActive())) continue;
2540 if (pt->GetBSigned()) continue;
2541 Double_t chi = pt->GetChi2()/pt->GetNumberOfClusters();
2542 if (chi>maxchi) continue;
2545 Float_t dens = pt->GetNumberOfClusters()/Float_t(pt->GetNFoundable());
2547 //sign only tracks with enoug big density at the beginning
2549 if ((pt->GetDensityFirst(40)<0.75) && pt->GetNumberOfClusters()<meann) continue;
2552 Double_t mindens = TMath::Max(double(mdensity-sdensity*fdensity*bfactor),0.65);
2553 Double_t minn = TMath::Max(Int_t(meann-fnumber*smeann*bfactor),50);
2555 // if (pt->fBConstrain) mindens = TMath::Max(mdensity-sdensity*fdensity*bfactor,0.65);
2556 if ( (pt->GetRemoval()==10) && (pt->GetSnp()>0.8)&&(dens>mindens))
2559 if ((dens>mindens && pt->GetNumberOfClusters()>minn) && chi<maxchi ){
2560 //Int_t noc=pt->GetNumberOfClusters();
2561 pt->SetBSigned(kTRUE);
2562 for (Int_t j=0; j<160; ++j) {
2564 Int_t index=pt->GetClusterIndex2(j);
2565 if (index<0) continue;
2566 AliTPCclusterMI *c= pt->GetClusterPointer(j);
2568 // if (!(c->IsUsed(10))) c->Use();
2573 // gLastCheck = nseed;
2581 void AliTPCtrackerMI::StopNotActive(TObjArray * arr, Int_t row0, Float_t th0, Float_t th1, Float_t th2) const
2583 // stop not active tracks
2584 // take th1 as threshold for number of founded to number of foundable on last 10 active rows
2585 // take th2 as threshold for number of founded to number of foundable on last 20 active rows
2586 Int_t nseed = arr->GetEntriesFast();
2588 for (Int_t i=0; i<nseed; i++) {
2589 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2593 if (!(pt->IsActive())) continue;
2594 StopNotActive(pt,row0,th0, th1,th2);
2600 void AliTPCtrackerMI::StopNotActive(AliTPCseed * seed, Int_t row0, Float_t th0, Float_t th1,
2603 // stop not active tracks
2604 // take th1 as threshold for number of founded to number of foundable on last 10 active rows
2605 // take th2 as threshold for number of founded to number of foundable on last 20 active rows
2608 Int_t foundable = 0;
2609 Int_t maxindex = seed->GetLastPoint(); //last foundable row
2610 if (seed->GetNFoundable()*th0 > seed->GetNumberOfClusters()) {
2611 seed->Desactivate(10) ;
2615 for (Int_t i=row0; i<maxindex; i++){
2616 Int_t index = seed->GetClusterIndex2(i);
2617 if (index!=-1) foundable++;
2619 if (foundable<=30) sumgood1++;
2620 if (foundable<=50) {
2627 if (foundable>=30.){
2628 if (sumgood1<(th1*30.)) seed->Desactivate(10);
2631 if (sumgood2<(th2*50.)) seed->Desactivate(10);
2635 Int_t AliTPCtrackerMI::RefitInward(AliESDEvent *event)
2638 // back propagation of ESD tracks
2641 const Int_t kMaxFriendTracks=2000;
2645 //PrepareForProlongation(fSeeds,1);
2646 PropagateForward2(fSeeds);
2647 RemoveUsed2(fSeeds,0.4,0.4,20);
2649 TObjArray arraySeed(fSeeds->GetEntries());
2650 for (Int_t i=0;i<fSeeds->GetEntries();i++) {
2651 arraySeed.AddAt(fSeeds->At(i),i);
2653 SignShared(&arraySeed);
2654 // FindCurling(fSeeds, event,2); // find multi found tracks
2655 FindSplitted(fSeeds, event,2); // find multi found tracks
2658 Int_t nseed = fSeeds->GetEntriesFast();
2659 for (Int_t i=0;i<nseed;i++){
2660 AliTPCseed * seed = (AliTPCseed*) fSeeds->UncheckedAt(i);
2661 if (!seed) continue;
2662 if (seed->GetKinkIndex(0)>0) UpdateKinkQualityD(seed); // update quality informations for kinks
2663 AliESDtrack *esd=event->GetTrack(i);
2665 if (seed->GetNumberOfClusters()<60 && seed->GetNumberOfClusters()<(esd->GetTPCclusters(0) -5)*0.8){
2666 AliExternalTrackParam paramIn;
2667 AliExternalTrackParam paramOut;
2668 Int_t ncl = seed->RefitTrack(seed,¶mIn,¶mOut);
2669 if (AliTPCReconstructor::StreamLevel()>0) {
2670 (*fDebugStreamer)<<"RecoverIn"<<
2674 "pout.="<<¶mOut<<
2679 seed->Set(paramIn.GetX(),paramIn.GetAlpha(),paramIn.GetParameter(),paramIn.GetCovariance());
2680 seed->SetNumberOfClusters(ncl);
2684 seed->PropagateTo(fParam->GetInnerRadiusLow());
2685 seed->UpdatePoints();
2686 AddCovariance(seed);
2688 seed->CookdEdx(0.02,0.6);
2689 CookLabel(seed,0.1); //For comparison only
2690 esd->SetTPCClusterMap(seed->GetClusterMap());
2691 esd->SetTPCSharedMap(seed->GetSharedMap());
2693 if (AliTPCReconstructor::StreamLevel()>1 && seed!=0&&esd!=0) {
2694 TTreeSRedirector &cstream = *fDebugStreamer;
2701 if (seed->GetNumberOfClusters()>15){
2702 esd->UpdateTrackParams(seed,AliESDtrack::kTPCrefit);
2703 esd->SetTPCPoints(seed->GetPoints());
2704 esd->SetTPCPointsF(seed->GetNFoundable());
2705 Int_t ndedx = seed->GetNCDEDX(0)+seed->GetNCDEDX(1)+seed->GetNCDEDX(2)+seed->GetNCDEDX(3);
2706 Float_t sdedx = (seed->GetSDEDX(0)+seed->GetSDEDX(1)+seed->GetSDEDX(2)+seed->GetSDEDX(3))*0.25;
2707 Float_t dedx = seed->GetdEdx();
2708 esd->SetTPCsignal(dedx, sdedx, ndedx);
2710 // add seed to the esd track in Calib level
2712 Bool_t storeFriend = gRandom->Rndm()<(kMaxFriendTracks)/Float_t(nseed);
2713 if (AliTPCReconstructor::StreamLevel()>0 &&storeFriend){
2714 AliTPCseed * seedCopy = new AliTPCseed(*seed, kTRUE);
2715 esd->AddCalibObject(seedCopy);
2720 //printf("problem\n");
2723 //FindKinks(fSeeds,event);
2724 Info("RefitInward","Number of refitted tracks %d",ntracks);
2729 Int_t AliTPCtrackerMI::PropagateBack(AliESDEvent *event)
2732 // back propagation of ESD tracks
2738 PropagateBack(fSeeds);
2739 RemoveUsed2(fSeeds,0.4,0.4,20);
2740 //FindCurling(fSeeds, fEvent,1);
2741 FindSplitted(fSeeds, event,1); // find multi found tracks
2744 Int_t nseed = fSeeds->GetEntriesFast();
2746 for (Int_t i=0;i<nseed;i++){
2747 AliTPCseed * seed = (AliTPCseed*) fSeeds->UncheckedAt(i);
2748 if (!seed) continue;
2749 if (seed->GetKinkIndex(0)<0) UpdateKinkQualityM(seed); // update quality informations for kinks
2750 seed->UpdatePoints();
2751 AddCovariance(seed);
2752 AliESDtrack *esd=event->GetTrack(i);
2753 if (seed->GetNumberOfClusters()<60 && seed->GetNumberOfClusters()<(esd->GetTPCclusters(0) -5)*0.8){
2754 AliExternalTrackParam paramIn;
2755 AliExternalTrackParam paramOut;
2756 Int_t ncl = seed->RefitTrack(seed,¶mIn,¶mOut);
2757 if (AliTPCReconstructor::StreamLevel()>0) {
2758 (*fDebugStreamer)<<"RecoverBack"<<
2762 "pout.="<<¶mOut<<
2767 seed->Set(paramOut.GetX(),paramOut.GetAlpha(),paramOut.GetParameter(),paramOut.GetCovariance());
2768 seed->SetNumberOfClusters(ncl);
2771 seed->CookdEdx(0.02,0.6);
2772 CookLabel(seed,0.1); //For comparison only
2773 if (seed->GetNumberOfClusters()>15){
2774 esd->UpdateTrackParams(seed,AliESDtrack::kTPCout);
2775 esd->SetTPCPoints(seed->GetPoints());
2776 esd->SetTPCPointsF(seed->GetNFoundable());
2777 Int_t ndedx = seed->GetNCDEDX(0)+seed->GetNCDEDX(1)+seed->GetNCDEDX(2)+seed->GetNCDEDX(3);
2778 Float_t sdedx = (seed->GetSDEDX(0)+seed->GetSDEDX(1)+seed->GetSDEDX(2)+seed->GetSDEDX(3))*0.25;
2779 Float_t dedx = seed->GetdEdx();
2780 esd->SetTPCsignal(dedx, sdedx, ndedx);
2782 Int_t eventnumber = event->GetEventNumberInFile();// patch 28 fev 06
2783 // This is most likely NOT the event number you'd like to use. It has nothing to do with the 'real' event number
2784 if (AliTPCReconstructor::StreamLevel()>1 && esd) {
2785 (*fDebugStreamer)<<"Cback"<<
2788 "EventNrInFile="<<eventnumber<<
2793 //FindKinks(fSeeds,event);
2794 Info("PropagateBack","Number of back propagated tracks %d",ntracks);
2801 void AliTPCtrackerMI::DeleteSeeds()
2806 Int_t nseed = fSeeds->GetEntriesFast();
2807 for (Int_t i=0;i<nseed;i++){
2808 AliTPCseed * seed = (AliTPCseed*)fSeeds->At(i);
2809 if (seed) delete fSeeds->RemoveAt(i);
2816 void AliTPCtrackerMI::ReadSeeds(AliESDEvent *event, Int_t direction)
2819 //read seeds from the event
2821 Int_t nentr=event->GetNumberOfTracks();
2823 Info("ReadSeeds", "Number of ESD tracks: %d\n", nentr);
2828 fSeeds = new TObjArray(nentr);
2832 for (Int_t i=0; i<nentr; i++) {
2833 AliESDtrack *esd=event->GetTrack(i);
2834 ULong_t status=esd->GetStatus();
2835 if (!(status&AliESDtrack::kTPCin)) continue;
2836 AliTPCtrack t(*esd);
2837 t.SetNumberOfClusters(0);
2838 // AliTPCseed *seed = new AliTPCseed(t,t.GetAlpha());
2839 AliTPCseed *seed = new AliTPCseed(t/*,t.GetAlpha()*/);
2840 seed->SetUniqueID(esd->GetID());
2841 AddCovariance(seed); //add systematic ucertainty
2842 for (Int_t ikink=0;ikink<3;ikink++) {
2843 Int_t index = esd->GetKinkIndex(ikink);
2844 seed->GetKinkIndexes()[ikink] = index;
2845 if (index==0) continue;
2846 index = TMath::Abs(index);
2847 AliESDkink * kink = fEvent->GetKink(index-1);
2848 if (kink&&esd->GetKinkIndex(ikink)<0){
2849 if ((status & AliESDtrack::kTRDrefit) != 0) kink->SetStatus(1,2);
2850 if ((status & AliESDtrack::kITSout) != 0) kink->SetStatus(1,0);
2852 if (kink&&esd->GetKinkIndex(ikink)>0){
2853 if ((status & AliESDtrack::kTRDrefit) != 0) kink->SetStatus(1,6);
2854 if ((status & AliESDtrack::kITSout) != 0) kink->SetStatus(1,4);
2858 if (((status&AliESDtrack::kITSout)==0)&&(direction==1)) seed->ResetCovariance(10.);
2859 if ( direction ==2 &&(status & AliESDtrack::kTRDrefit) == 0 ) seed->ResetCovariance(10.);
2860 //if ( direction ==2 && ((status & AliESDtrack::kTPCout) == 0) ) {
2861 // fSeeds->AddAt(0,i);
2865 if ( direction ==2 &&(status & AliESDtrack::kTRDrefit) > 0 ) {
2866 Double_t par0[5],par1[5],alpha,x;
2867 esd->GetInnerExternalParameters(alpha,x,par0);
2868 esd->GetExternalParameters(x,par1);
2869 Double_t delta1 = TMath::Abs(par0[4]-par1[4])/(0.000000001+TMath::Abs(par0[4]+par1[4]));
2870 Double_t delta2 = TMath::Abs(par0[3]-par1[3]);
2872 if (esd->GetTRDncls()>0) trdchi2 = esd->GetTRDchi2()/esd->GetTRDncls();
2873 //reset covariance if suspicious
2874 if ( (delta1>0.1) || (delta2>0.006) ||trdchi2>7.)
2875 seed->ResetCovariance(10.);
2880 // rotate to the local coordinate system
2882 fSectors=fInnerSec; fN=fkNIS;
2883 Double_t alpha=seed->GetAlpha() - fSectors->GetAlphaShift();
2884 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
2885 if (alpha < 0. ) alpha += 2.*TMath::Pi();
2886 Int_t ns=Int_t(alpha/fSectors->GetAlpha())%fN;
2887 alpha =ns*fSectors->GetAlpha() + fSectors->GetAlphaShift();
2888 if (alpha<-TMath::Pi()) alpha += 2*TMath::Pi();
2889 if (alpha>=TMath::Pi()) alpha -= 2*TMath::Pi();
2890 alpha-=seed->GetAlpha();
2891 if (!seed->Rotate(alpha)) {
2897 if (esd->GetKinkIndex(0)<=0){
2898 for (Int_t irow=0;irow<160;irow++){
2899 Int_t index = seed->GetClusterIndex2(irow);
2902 AliTPCclusterMI * cl = GetClusterMI(index);
2903 seed->SetClusterPointer(irow,cl);
2905 if ((index & 0x8000)==0){
2906 cl->Use(10); // accepted cluster
2908 cl->Use(6); // close cluster not accepted
2911 Info("ReadSeeds","Not found cluster");
2916 fSeeds->AddAt(seed,i);
2922 //_____________________________________________________________________________
2923 void AliTPCtrackerMI::MakeSeeds3(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, Float_t cuts[4],
2924 Float_t deltay, Int_t ddsec) {
2925 //-----------------------------------------------------------------
2926 // This function creates track seeds.
2927 // SEEDING WITH VERTEX CONSTRAIN
2928 //-----------------------------------------------------------------
2929 // cuts[0] - fP4 cut
2930 // cuts[1] - tan(phi) cut
2931 // cuts[2] - zvertex cut
2932 // cuts[3] - fP3 cut
2940 Double_t x[5], c[15];
2941 // Int_t di = i1-i2;
2943 AliTPCseed * seed = new AliTPCseed();
2944 Double_t alpha=fSectors->GetAlpha(), shift=fSectors->GetAlphaShift();
2945 Double_t cs=cos(alpha), sn=sin(alpha);
2947 // Double_t x1 =fOuterSec->GetX(i1);
2948 //Double_t xx2=fOuterSec->GetX(i2);
2950 Double_t x1 =GetXrow(i1);
2951 Double_t xx2=GetXrow(i2);
2953 Double_t x3=GetX(), y3=GetY(), z3=GetZ();
2955 Int_t imiddle = (i2+i1)/2; //middle pad row index
2956 Double_t xm = GetXrow(imiddle); // radius of middle pad-row
2957 const AliTPCtrackerRow& krm=GetRow(sec,imiddle); //middle pad -row
2961 const AliTPCtrackerRow& kr1=GetRow(ns,i1);
2962 Double_t ymax = GetMaxY(i1)-kr1.GetDeadZone()-1.5;
2963 Double_t ymaxm = GetMaxY(imiddle)-kr1.GetDeadZone()-1.5;
2966 // change cut on curvature if it can't reach this layer
2967 // maximal curvature set to reach it
2968 Double_t dvertexmax = TMath::Sqrt((x1-x3)*(x1-x3)+(ymax+5-y3)*(ymax+5-y3));
2969 if (dvertexmax*0.5*cuts[0]>0.85){
2970 cuts[0] = 0.85/(dvertexmax*0.5+1.);
2972 Double_t r2min = 1/(cuts[0]*cuts[0]); //minimal square of radius given by cut
2975 if (deltay>0) ddsec = 0;
2976 // loop over clusters
2977 for (Int_t is=0; is < kr1; is++) {
2979 if (kr1[is]->IsUsed(10)) continue;
2980 Double_t y1=kr1[is]->GetY(), z1=kr1[is]->GetZ();
2981 //if (TMath::Abs(y1)>ymax) continue;
2983 if (deltay>0 && TMath::Abs(ymax-TMath::Abs(y1))> deltay ) continue; // seed only at the edge
2985 // find possible directions
2986 Float_t anglez = (z1-z3)/(x1-x3);
2987 Float_t extraz = z1 - anglez*(x1-xx2); // extrapolated z
2990 //find rotation angles relative to line given by vertex and point 1
2991 Double_t dvertex2 = (x1-x3)*(x1-x3)+(y1-y3)*(y1-y3);
2992 Double_t dvertex = TMath::Sqrt(dvertex2);
2993 Double_t angle13 = TMath::ATan((y1-y3)/(x1-x3));
2994 Double_t cs13 = cos(-angle13), sn13 = sin(-angle13);
2997 // loop over 2 sectors
3003 Double_t dddz1=0; // direction of delta inclination in z axis
3010 for (Int_t dsec = dsec1; dsec<=dsec2;dsec++){
3011 Int_t sec2 = sec + dsec;
3013 // AliTPCtrackerRow& kr2 = fOuterSec[(sec2+fkNOS)%fkNOS][i2];
3014 //AliTPCtrackerRow& kr2m = fOuterSec[(sec2+fkNOS)%fkNOS][imiddle];
3015 AliTPCtrackerRow& kr2 = GetRow((sec2+fkNOS)%fkNOS,i2);
3016 AliTPCtrackerRow& kr2m = GetRow((sec2+fkNOS)%fkNOS,imiddle);
3017 Int_t index1 = TMath::Max(kr2.Find(extraz-0.6-dddz1*TMath::Abs(z1)*0.05)-1,0);
3018 Int_t index2 = TMath::Min(kr2.Find(extraz+0.6+dddz2*TMath::Abs(z1)*0.05)+1,kr2);
3020 // rotation angles to p1-p3
3021 Double_t cs13r = cos(-angle13+dsec*alpha)/dvertex, sn13r = sin(-angle13+dsec*alpha)/dvertex;
3022 Double_t x2, y2, z2;
3024 // Double_t dymax = maxangle*TMath::Abs(x1-xx2);
3027 Double_t dxx0 = (xx2-x3)*cs13r;
3028 Double_t dyy0 = (xx2-x3)*sn13r;
3029 for (Int_t js=index1; js < index2; js++) {
3030 const AliTPCclusterMI *kcl = kr2[js];
3031 if (kcl->IsUsed(10)) continue;
3033 //calcutate parameters
3035 Double_t yy0 = dyy0 +(kcl->GetY()-y3)*cs13r;
3037 if (TMath::Abs(yy0)<0.000001) continue;
3038 Double_t xx0 = dxx0 -(kcl->GetY()-y3)*sn13r;
3039 Double_t y0 = 0.5*(xx0*xx0+yy0*yy0-xx0)/yy0;
3040 Double_t r02 = (0.25+y0*y0)*dvertex2;
3041 //curvature (radius) cut
3042 if (r02<r2min) continue;
3046 Double_t c0 = 1/TMath::Sqrt(r02);
3050 //Double_t dfi0 = 2.*TMath::ASin(dvertex*c0*0.5);
3051 //Double_t dfi1 = 2.*TMath::ASin(TMath::Sqrt(yy0*yy0+(1-xx0)*(1-xx0))*dvertex*c0*0.5);
3052 Double_t dfi0 = 2.*AliTPCFastMath::FastAsin(dvertex*c0*0.5);
3053 Double_t dfi1 = 2.*AliTPCFastMath::FastAsin(TMath::Sqrt(yy0*yy0+(1-xx0)*(1-xx0))*dvertex*c0*0.5);
3056 Double_t z0 = kcl->GetZ();
3057 Double_t zzzz2 = z1-(z1-z3)*dfi1/dfi0;
3058 if (TMath::Abs(zzzz2-z0)>0.5) continue;
3061 Double_t dip = (z1-z0)*c0/dfi1;
3062 Double_t x0 = (0.5*cs13+y0*sn13)*dvertex*c0;
3073 x2= xx2*cs-y2*sn*dsec;
3074 y2=+xx2*sn*dsec+y2*cs;
3084 // do we have cluster at the middle ?
3086 GetProlongation(x1,xm,x,ym,zm);
3088 AliTPCclusterMI * cm=0;
3089 if (TMath::Abs(ym)-ymaxm<0){
3090 cm = krm.FindNearest2(ym,zm,1.0,0.6,dummy);
3091 if ((!cm) || (cm->IsUsed(10))) {
3096 // rotate y1 to system 0
3097 // get state vector in rotated system
3098 Double_t yr1 = (-0.5*sn13+y0*cs13)*dvertex*c0;
3099 Double_t xr2 = x0*cs+yr1*sn*dsec;
3100 Double_t xr[5]={kcl->GetY(),kcl->GetZ(), xr2, dip, c0};
3102 GetProlongation(xx2,xm,xr,ym,zm);
3103 if (TMath::Abs(ym)-ymaxm<0){
3104 cm = kr2m.FindNearest2(ym,zm,1.0,0.6,dummy);
3105 if ((!cm) || (cm->IsUsed(10))) {
3115 dym = ym - cm->GetY();
3116 dzm = zm - cm->GetZ();
3123 Double_t sy1=kr1[is]->GetSigmaY2()*2., sz1=kr1[is]->GetSigmaZ2()*2.;
3124 Double_t sy2=kcl->GetSigmaY2()*2., sz2=kcl->GetSigmaZ2()*2.;
3125 //Double_t sy3=400*3./12., sy=0.1, sz=0.1;
3126 Double_t sy3=25000*x[4]*x[4]+0.1, sy=0.1, sz=0.1;
3127 //Double_t sy3=25000*x[4]*x[4]*60+0.5, sy=0.1, sz=0.1;
3129 Double_t f40=(F1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
3130 Double_t f42=(F1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
3131 Double_t f43=(F1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
3132 Double_t f20=(F2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
3133 Double_t f22=(F2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
3134 Double_t f23=(F2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
3136 Double_t f30=(F3(x1,y1+sy,x2,y2,z1,z2)-x[3])/sy;
3137 Double_t f31=(F3(x1,y1,x2,y2,z1+sz,z2)-x[3])/sz;
3138 Double_t f32=(F3(x1,y1,x2,y2+sy,z1,z2)-x[3])/sy;
3139 Double_t f34=(F3(x1,y1,x2,y2,z1,z2+sz)-x[3])/sz;
3143 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3144 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3145 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3146 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3147 c[13]=f30*sy1*f40+f32*sy2*f42;
3148 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3150 // if (!BuildSeed(kr1[is],kcl,0,x1,x2,x3,x,c)) continue;
3152 UInt_t index=kr1.GetIndex(is);
3153 seed->~AliTPCseed(); // this does not set the pointer to 0...
3154 AliTPCseed *track=new(seed) AliTPCseed(x1, ns*alpha+shift, x, c, index);
3156 track->SetIsSeeding(kTRUE);
3157 track->SetSeed1(i1);
3158 track->SetSeed2(i2);
3159 track->SetSeedType(3);
3163 FollowProlongation(*track, (i1+i2)/2,1);
3164 Int_t foundable,found,shared;
3165 track->GetClusterStatistic((i1+i2)/2,i1, found, foundable, shared, kTRUE);
3166 if ((found<0.55*foundable) || shared>0.5*found || (track->GetSigmaY2()+track->GetSigmaZ2())>0.5){
3168 seed->~AliTPCseed();
3174 FollowProlongation(*track, i2,1);
3178 track->SetBConstrain(1);
3179 // track->fLastPoint = i1+fInnerSec->GetNRows(); // first cluster in track position
3180 track->SetLastPoint(i1); // first cluster in track position
3181 track->SetFirstPoint(track->GetLastPoint());
3183 if (track->GetNumberOfClusters()<(i1-i2)*0.5 ||
3184 track->GetNumberOfClusters() < track->GetNFoundable()*0.6 ||
3185 track->GetNShared()>0.4*track->GetNumberOfClusters() ) {
3187 seed->~AliTPCseed();
3191 // Z VERTEX CONDITION
3192 Double_t zv, bz=GetBz();
3193 if ( !track->GetZAt(0.,bz,zv) ) continue;
3194 if (TMath::Abs(zv-z3)>cuts[2]) {
3195 FollowProlongation(*track, TMath::Max(i2-20,0));
3196 if ( !track->GetZAt(0.,bz,zv) ) continue;
3197 if (TMath::Abs(zv-z3)>cuts[2]){
3198 FollowProlongation(*track, TMath::Max(i2-40,0));
3199 if ( !track->GetZAt(0.,bz,zv) ) continue;
3200 if (TMath::Abs(zv-z3)>cuts[2] &&(track->GetNumberOfClusters() > track->GetNFoundable()*0.7)){
3201 // make seed without constrain
3202 AliTPCseed * track2 = MakeSeed(track,0.2,0.5,1.);
3203 FollowProlongation(*track2, i2,1);
3204 track2->SetBConstrain(kFALSE);
3205 track2->SetSeedType(1);
3206 arr->AddLast(track2);
3208 seed->~AliTPCseed();
3213 seed->~AliTPCseed();
3220 track->SetSeedType(0);
3221 arr->AddLast(track);
3222 seed = new AliTPCseed;
3224 // don't consider other combinations
3225 if (track->GetNumberOfClusters() > track->GetNFoundable()*0.8)
3231 Info("MakeSeeds3","\nSeeding statistic:\t%d\t%d\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin,nout1,nout2);
3237 void AliTPCtrackerMI::MakeSeeds5(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, Float_t cuts[4],
3242 //-----------------------------------------------------------------
3243 // This function creates track seeds.
3244 //-----------------------------------------------------------------
3245 // cuts[0] - fP4 cut
3246 // cuts[1] - tan(phi) cut
3247 // cuts[2] - zvertex cut
3248 // cuts[3] - fP3 cut
3258 Double_t x[5], c[15];
3260 // make temporary seed
3261 AliTPCseed * seed = new AliTPCseed;
3262 Double_t alpha=fOuterSec->GetAlpha(), shift=fOuterSec->GetAlphaShift();
3263 // Double_t cs=cos(alpha), sn=sin(alpha);
3268 Double_t x1 = GetXrow(i1-1);
3269 const AliTPCtrackerRow& kr1=GetRow(sec,i1-1);
3270 Double_t y1max = GetMaxY(i1-1)-kr1.GetDeadZone()-1.5;
3272 Double_t x1p = GetXrow(i1);
3273 const AliTPCtrackerRow& kr1p=GetRow(sec,i1);
3275 Double_t x1m = GetXrow(i1-2);
3276 const AliTPCtrackerRow& kr1m=GetRow(sec,i1-2);
3279 //last 3 padrow for seeding
3280 AliTPCtrackerRow& kr3 = GetRow((sec+fkNOS)%fkNOS,i1-7);
3281 Double_t x3 = GetXrow(i1-7);
3282 // Double_t y3max= GetMaxY(i1-7)-kr3.fDeadZone-1.5;
3284 AliTPCtrackerRow& kr3p = GetRow((sec+fkNOS)%fkNOS,i1-6);
3285 Double_t x3p = GetXrow(i1-6);
3287 AliTPCtrackerRow& kr3m = GetRow((sec+fkNOS)%fkNOS,i1-8);
3288 Double_t x3m = GetXrow(i1-8);
3293 Int_t im = i1-4; //middle pad row index
3294 Double_t xm = GetXrow(im); // radius of middle pad-row
3295 const AliTPCtrackerRow& krm=GetRow(sec,im); //middle pad -row
3296 // Double_t ymmax = GetMaxY(im)-kr1.fDeadZone-1.5;
3299 Double_t deltax = x1-x3;
3300 Double_t dymax = deltax*cuts[1];
3301 Double_t dzmax = deltax*cuts[3];
3303 // loop over clusters
3304 for (Int_t is=0; is < kr1; is++) {
3306 if (kr1[is]->IsUsed(10)) continue;
3307 Double_t y1=kr1[is]->GetY(), z1=kr1[is]->GetZ();
3309 if (deltay>0 && TMath::Abs(y1max-TMath::Abs(y1))> deltay ) continue; // seed only at the edge
3311 Int_t index1 = TMath::Max(kr3.Find(z1-dzmax)-1,0);
3312 Int_t index2 = TMath::Min(kr3.Find(z1+dzmax)+1,kr3);
3318 for (Int_t js=index1; js < index2; js++) {
3319 const AliTPCclusterMI *kcl = kr3[js];
3320 if (kcl->IsUsed(10)) continue;
3322 // apply angular cuts
3323 if (TMath::Abs(y1-y3)>dymax) continue;
3326 if (TMath::Abs(z1-z3)>dzmax) continue;
3328 Double_t angley = (y1-y3)/(x1-x3);
3329 Double_t anglez = (z1-z3)/(x1-x3);
3331 Double_t erry = TMath::Abs(angley)*(x1-x1m)*0.5+0.5;
3332 Double_t errz = TMath::Abs(anglez)*(x1-x1m)*0.5+0.5;
3334 Double_t yyym = angley*(xm-x1)+y1;
3335 Double_t zzzm = anglez*(xm-x1)+z1;
3337 const AliTPCclusterMI *kcm = krm.FindNearest2(yyym,zzzm,erry,errz,index);
3339 if (kcm->IsUsed(10)) continue;
3341 erry = TMath::Abs(angley)*(x1-x1m)*0.4+0.5;
3342 errz = TMath::Abs(anglez)*(x1-x1m)*0.4+0.5;
3349 // look around first
3350 const AliTPCclusterMI *kc1m = kr1m.FindNearest2(angley*(x1m-x1)+y1,
3356 if (kc1m->IsUsed(10)) used++;
3358 const AliTPCclusterMI *kc1p = kr1p.FindNearest2(angley*(x1p-x1)+y1,
3364 if (kc1p->IsUsed(10)) used++;
3366 if (used>1) continue;
3367 if (found<1) continue;
3371 const AliTPCclusterMI *kc3m = kr3m.FindNearest2(angley*(x3m-x3)+y3,
3377 if (kc3m->IsUsed(10)) used++;
3381 const AliTPCclusterMI *kc3p = kr3p.FindNearest2(angley*(x3p-x3)+y3,
3387 if (kc3p->IsUsed(10)) used++;
3391 if (used>1) continue;
3392 if (found<3) continue;
3402 x[4]=F1(x1,y1,x2,y2,x3,y3);
3403 //if (TMath::Abs(x[4]) >= cuts[0]) continue;
3406 x[2]=F2(x1,y1,x2,y2,x3,y3);
3409 x[3]=F3n(x1,y1,x2,y2,z1,z2,x[4]);
3410 //if (TMath::Abs(x[3]) > cuts[3]) continue;
3414 Double_t sy1=0.1, sz1=0.1;
3415 Double_t sy2=0.1, sz2=0.1;
3416 Double_t sy3=0.1, sy=0.1, sz=0.1;
3418 Double_t f40=(F1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
3419 Double_t f42=(F1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
3420 Double_t f43=(F1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
3421 Double_t f20=(F2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
3422 Double_t f22=(F2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
3423 Double_t f23=(F2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
3425 Double_t f30=(F3(x1,y1+sy,x2,y2,z1,z2)-x[3])/sy;
3426 Double_t f31=(F3(x1,y1,x2,y2,z1+sz,z2)-x[3])/sz;
3427 Double_t f32=(F3(x1,y1,x2,y2+sy,z1,z2)-x[3])/sy;
3428 Double_t f34=(F3(x1,y1,x2,y2,z1,z2+sz)-x[3])/sz;
3432 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3433 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3434 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3435 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3436 c[13]=f30*sy1*f40+f32*sy2*f42;
3437 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3439 // if (!BuildSeed(kr1[is],kcl,0,x1,x2,x3,x,c)) continue;
3441 index=kr1.GetIndex(is);
3442 seed->~AliTPCseed();
3443 AliTPCseed *track=new(seed) AliTPCseed(x1, sec*alpha+shift, x, c, index);
3445 track->SetIsSeeding(kTRUE);
3448 FollowProlongation(*track, i1-7,1);
3449 if (track->GetNumberOfClusters() < track->GetNFoundable()*0.75 ||
3450 track->GetNShared()>0.6*track->GetNumberOfClusters() || ( track->GetSigmaY2()+ track->GetSigmaZ2())>0.6){
3452 seed->~AliTPCseed();
3458 FollowProlongation(*track, i2,1);
3459 track->SetBConstrain(0);
3460 track->SetLastPoint(i1+fInnerSec->GetNRows()); // first cluster in track position
3461 track->SetFirstPoint(track->GetLastPoint());
3463 if (track->GetNumberOfClusters()<(i1-i2)*0.5 ||
3464 track->GetNumberOfClusters()<track->GetNFoundable()*0.7 ||
3465 track->GetNShared()>2. || track->GetChi2()/track->GetNumberOfClusters()>6 || ( track->GetSigmaY2()+ track->GetSigmaZ2())>0.5 ) {
3467 seed->~AliTPCseed();
3472 FollowProlongation(*track, TMath::Max(i2-10,0),1);
3473 AliTPCseed * track2 = MakeSeed(track,0.2,0.5,0.9);
3474 FollowProlongation(*track2, i2,1);
3475 track2->SetBConstrain(kFALSE);
3476 track2->SetSeedType(4);
3477 arr->AddLast(track2);
3479 seed->~AliTPCseed();
3483 //arr->AddLast(track);
3484 //seed = new AliTPCseed;
3490 Info("MakeSeeds5","\nSeeding statiistic:\t%d\t%d\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin,nout1,nout2,nout3);
3496 //_____________________________________________________________________________
3497 void AliTPCtrackerMI::MakeSeeds2(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, Float_t */*cuts[4]*/,
3498 Float_t deltay, Bool_t /*bconstrain*/) {
3499 //-----------------------------------------------------------------
3500 // This function creates track seeds - without vertex constraint
3501 //-----------------------------------------------------------------
3502 // cuts[0] - fP4 cut - not applied
3503 // cuts[1] - tan(phi) cut
3504 // cuts[2] - zvertex cut - not applied
3505 // cuts[3] - fP3 cut
3515 Double_t alpha=fOuterSec->GetAlpha(), shift=fOuterSec->GetAlphaShift();
3516 // Double_t cs=cos(alpha), sn=sin(alpha);
3517 Int_t row0 = (i1+i2)/2;
3518 Int_t drow = (i1-i2)/2;
3519 const AliTPCtrackerRow& kr0=fSectors[sec][row0];
3520 AliTPCtrackerRow * kr=0;
3522 AliTPCpolyTrack polytrack;
3523 Int_t nclusters=fSectors[sec][row0];
3524 AliTPCseed * seed = new AliTPCseed;
3529 for (Int_t is=0; is < nclusters; is++) { //LOOP over clusters
3531 Int_t nfoundable =0;
3532 for (Int_t iter =1; iter<2; iter++){ //iterations
3533 const AliTPCtrackerRow& krm=fSectors[sec][row0-iter];
3534 const AliTPCtrackerRow& krp=fSectors[sec][row0+iter];
3535 const AliTPCclusterMI * cl= kr0[is];
3537 if (cl->IsUsed(10)) {
3543 Double_t x = kr0.GetX();
3544 // Initialization of the polytrack
3549 Double_t y0= cl->GetY();
3550 Double_t z0= cl->GetZ();
3554 Double_t ymax = fSectors->GetMaxY(row0)-kr0.GetDeadZone()-1.5;
3555 if (deltay>0 && TMath::Abs(ymax-TMath::Abs(y0))> deltay ) continue; // seed only at the edge
3557 erry = (0.5)*cl->GetSigmaY2()/TMath::Sqrt(cl->GetQ())*6;
3558 errz = (0.5)*cl->GetSigmaZ2()/TMath::Sqrt(cl->GetQ())*6;
3559 polytrack.AddPoint(x,y0,z0,erry, errz);
3562 if (cl->IsUsed(10)) sumused++;
3565 Float_t roady = (5*TMath::Sqrt(cl->GetSigmaY2()+0.2)+1.)*iter;
3566 Float_t roadz = (5*TMath::Sqrt(cl->GetSigmaZ2()+0.2)+1.)*iter;
3569 AliTPCclusterMI * cl1 = krm.FindNearest(y0,z0,roady,roadz);
3570 if (cl1 && TMath::Abs(ymax-TMath::Abs(y0))) {
3571 erry = (0.5)*cl1->GetSigmaY2()/TMath::Sqrt(cl1->GetQ())*3;
3572 errz = (0.5)*cl1->GetSigmaZ2()/TMath::Sqrt(cl1->GetQ())*3;
3573 if (cl1->IsUsed(10)) sumused++;
3574 polytrack.AddPoint(x,cl1->GetY(),cl1->GetZ(),erry,errz);
3578 AliTPCclusterMI * cl2 = krp.FindNearest(y0,z0,roady,roadz);
3580 erry = (0.5)*cl2->GetSigmaY2()/TMath::Sqrt(cl2->GetQ())*3;
3581 errz = (0.5)*cl2->GetSigmaZ2()/TMath::Sqrt(cl2->GetQ())*3;
3582 if (cl2->IsUsed(10)) sumused++;
3583 polytrack.AddPoint(x,cl2->GetY(),cl2->GetZ(),erry,errz);
3586 if (sumused>0) continue;
3588 polytrack.UpdateParameters();
3594 nfoundable = polytrack.GetN();
3595 nfound = nfoundable;
3597 for (Int_t ddrow = iter+1; ddrow<drow;ddrow++){
3598 Float_t maxdist = 0.8*(1.+3./(ddrow));
3599 for (Int_t delta = -1;delta<=1;delta+=2){
3600 Int_t row = row0+ddrow*delta;
3601 kr = &(fSectors[sec][row]);
3602 Double_t xn = kr->GetX();
3603 Double_t ymax1 = fSectors->GetMaxY(row)-kr->GetDeadZone()-1.5;
3604 polytrack.GetFitPoint(xn,yn,zn);
3605 if (TMath::Abs(yn)>ymax1) continue;
3607 AliTPCclusterMI * cln = kr->FindNearest(yn,zn,roady,roadz);
3609 Float_t dist = TMath::Sqrt( (yn-cln->GetY())*(yn-cln->GetY())+(zn-cln->GetZ())*(zn-cln->GetZ()));
3612 erry = (dist+0.3)*cln->GetSigmaY2()/TMath::Sqrt(cln->GetQ())*(1.+1./(ddrow));
3613 errz = (dist+0.3)*cln->GetSigmaZ2()/TMath::Sqrt(cln->GetQ())*(1.+1./(ddrow));
3614 if (cln->IsUsed(10)) {
3615 // printf("used\n");
3623 polytrack.AddPoint(xn,cln->GetY(),cln->GetZ(),erry, errz);
3628 if ( (sumused>3) || (sumused>0.5*nfound) || (nfound<0.6*nfoundable)) break;
3629 polytrack.UpdateParameters();
3632 if ( (sumused>3) || (sumused>0.5*nfound)) {
3633 //printf("sumused %d\n",sumused);
3638 polytrack.GetFitDerivation(kr0.GetX(),dy,dz);
3639 AliTPCpolyTrack track2;
3641 polytrack.Refit(track2,0.5+TMath::Abs(dy)*0.3,0.4+TMath::Abs(dz)*0.3);
3642 if (track2.GetN()<0.5*nfoundable) continue;
3645 if ((nfound>0.6*nfoundable) &&( nfoundable>0.4*(i1-i2))) {
3647 // test seed with and without constrain
3648 for (Int_t constrain=0; constrain<=0;constrain++){
3649 // add polytrack candidate
3651 Double_t x[5], c[15];
3652 Double_t x1,x2,x3,y1,y2,y3,z1,z2,z3;
3653 track2.GetBoundaries(x3,x1);
3655 track2.GetFitPoint(x1,y1,z1);
3656 track2.GetFitPoint(x2,y2,z2);
3657 track2.GetFitPoint(x3,y3,z3);
3659 //is track pointing to the vertex ?
3662 polytrack.GetFitPoint(x0,y0,z0);
3675 x[4]=F1(x1,y1,x2,y2,x3,y3);
3677 // if (TMath::Abs(x[4]) >= cuts[0]) continue; //
3678 x[2]=F2(x1,y1,x2,y2,x3,y3);
3680 //if (TMath::Abs(x[4]*x1-x[2]) >= cuts[1]) continue;
3681 //x[3]=F3(x1,y1,x2,y2,z1,z2);
3682 x[3]=F3n(x1,y1,x3,y3,z1,z3,x[4]);
3683 //if (TMath::Abs(x[3]) > cuts[3]) continue;
3686 Double_t sy =0.1, sz =0.1;
3687 Double_t sy1=0.02, sz1=0.02;
3688 Double_t sy2=0.02, sz2=0.02;
3692 sy3=25000*x[4]*x[4]+0.1, sy=0.1, sz=0.1;
3695 Double_t f40=(F1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
3696 Double_t f42=(F1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
3697 Double_t f43=(F1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
3698 Double_t f20=(F2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
3699 Double_t f22=(F2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
3700 Double_t f23=(F2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
3702 Double_t f30=(F3(x1,y1+sy,x3,y3,z1,z3)-x[3])/sy;
3703 Double_t f31=(F3(x1,y1,x3,y3,z1+sz,z3)-x[3])/sz;
3704 Double_t f32=(F3(x1,y1,x3,y3+sy,z1,z3)-x[3])/sy;
3705 Double_t f34=(F3(x1,y1,x3,y3,z1,z3+sz)-x[3])/sz;
3710 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3711 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3712 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3713 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3714 c[13]=f30*sy1*f40+f32*sy2*f42;
3715 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3717 //Int_t row1 = fSectors->GetRowNumber(x1);
3718 Int_t row1 = GetRowNumber(x1);
3722 seed->~AliTPCseed();
3723 AliTPCseed *track=new(seed) AliTPCseed(x1,sec*alpha+shift,x,c,index);
3724 track->SetIsSeeding(kTRUE);
3725 Int_t rc=FollowProlongation(*track, i2);
3726 if (constrain) track->SetBConstrain(1);
3728 track->SetBConstrain(0);
3729 track->SetLastPoint(row1+fInnerSec->GetNRows()); // first cluster in track position
3730 track->SetFirstPoint(track->GetLastPoint());
3732 if (rc==0 || track->GetNumberOfClusters()<(i1-i2)*0.5 ||
3733 track->GetNumberOfClusters() < track->GetNFoundable()*0.6 ||
3734 track->GetNShared()>0.4*track->GetNumberOfClusters()) {
3737 seed->~AliTPCseed();
3740 arr->AddLast(track);
3741 seed = new AliTPCseed;
3745 } // if accepted seed
3748 Info("MakeSeeds2","\nSeeding statiistic:\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin3);
3754 AliTPCseed *AliTPCtrackerMI::MakeSeed(AliTPCseed *track, Float_t r0, Float_t r1, Float_t r2)
3758 //reseed using track points
3759 Int_t p0 = int(r0*track->GetNumberOfClusters()); // point 0
3760 Int_t p1 = int(r1*track->GetNumberOfClusters());
3761 Int_t p2 = int(r2*track->GetNumberOfClusters()); // last point
3763 Double_t x0[3],x1[3],x2[3];
3764 for (Int_t i=0;i<3;i++){
3770 // find track position at given ratio of the length
3771 Int_t sec0=0, sec1=0, sec2=0;
3774 for (Int_t i=0;i<160;i++){
3775 if (track->GetClusterPointer(i)){
3777 AliTPCTrackerPoint *trpoint =track->GetTrackPoint(i);
3778 if ( (index<p0) || x0[0]<0 ){
3779 if (trpoint->GetX()>1){
3780 clindex = track->GetClusterIndex2(i);
3782 x0[0] = trpoint->GetX();
3783 x0[1] = trpoint->GetY();
3784 x0[2] = trpoint->GetZ();
3785 sec0 = ((clindex&0xff000000)>>24)%18;
3790 if ( (index<p1) &&(trpoint->GetX()>1)){
3791 clindex = track->GetClusterIndex2(i);
3793 x1[0] = trpoint->GetX();
3794 x1[1] = trpoint->GetY();
3795 x1[2] = trpoint->GetZ();
3796 sec1 = ((clindex&0xff000000)>>24)%18;
3799 if ( (index<p2) &&(trpoint->GetX()>1)){
3800 clindex = track->GetClusterIndex2(i);
3802 x2[0] = trpoint->GetX();
3803 x2[1] = trpoint->GetY();
3804 x2[2] = trpoint->GetZ();
3805 sec2 = ((clindex&0xff000000)>>24)%18;
3812 Double_t alpha, cs,sn, xx2,yy2;
3814 alpha = (sec1-sec2)*fSectors->GetAlpha();
3815 cs = TMath::Cos(alpha);
3816 sn = TMath::Sin(alpha);
3817 xx2= x1[0]*cs-x1[1]*sn;
3818 yy2= x1[0]*sn+x1[1]*cs;
3822 alpha = (sec0-sec2)*fSectors->GetAlpha();
3823 cs = TMath::Cos(alpha);
3824 sn = TMath::Sin(alpha);
3825 xx2= x0[0]*cs-x0[1]*sn;
3826 yy2= x0[0]*sn+x0[1]*cs;
3832 Double_t x[5],c[15];
3836 x[4]=F1(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]);
3837 // if (x[4]>1) return 0;
3838 x[2]=F2(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]);
3839 x[3]=F3n(x2[0],x2[1],x0[0],x0[1],x2[2],x0[2],x[4]);
3840 //if (TMath::Abs(x[3]) > 2.2) return 0;
3841 //if (TMath::Abs(x[2]) > 1.99) return 0;
3843 Double_t sy =0.1, sz =0.1;
3845 Double_t sy1=0.02+track->GetSigmaY2(), sz1=0.02+track->GetSigmaZ2();
3846 Double_t sy2=0.01+track->GetSigmaY2(), sz2=0.01+track->GetSigmaZ2();
3847 Double_t sy3=0.01+track->GetSigmaY2();
3849 Double_t f40=(F1(x2[0],x2[1]+sy,x1[0],x1[1],x0[0],x0[1])-x[4])/sy;
3850 Double_t f42=(F1(x2[0],x2[1],x1[0],x1[1]+sy,x0[0],x0[1])-x[4])/sy;
3851 Double_t f43=(F1(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]+sy)-x[4])/sy;
3852 Double_t f20=(F2(x2[0],x2[1]+sy,x1[0],x1[1],x0[0],x0[1])-x[2])/sy;
3853 Double_t f22=(F2(x2[0],x2[1],x1[0],x1[1]+sy,x0[0],x0[1])-x[2])/sy;
3854 Double_t f23=(F2(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]+sy)-x[2])/sy;
3856 Double_t f30=(F3(x2[0],x2[1]+sy,x0[0],x0[1],x2[2],x0[2])-x[3])/sy;
3857 Double_t f31=(F3(x2[0],x2[1],x0[0],x0[1],x2[2]+sz,x0[2])-x[3])/sz;
3858 Double_t f32=(F3(x2[0],x2[1],x0[0],x0[1]+sy,x2[2],x0[2])-x[3])/sy;
3859 Double_t f34=(F3(x2[0],x2[1],x0[0],x0[1],x2[2],x0[2]+sz)-x[3])/sz;
3864 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3865 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3866 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3867 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3868 c[13]=f30*sy1*f40+f32*sy2*f42;
3869 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3871 // Int_t row1 = fSectors->GetRowNumber(x2[0]);
3872 AliTPCseed *seed=new AliTPCseed(x2[0], sec2*fSectors->GetAlpha()+fSectors->GetAlphaShift(), x, c, 0);
3873 // Double_t y0,z0,y1,z1, y2,z2;
3874 //seed->GetProlongation(x0[0],y0,z0);
3875 // seed->GetProlongation(x1[0],y1,z1);
3876 //seed->GetProlongation(x2[0],y2,z2);
3878 seed->SetLastPoint(pp2);
3879 seed->SetFirstPoint(pp2);
3886 AliTPCseed *AliTPCtrackerMI::ReSeed(AliTPCseed *track, Float_t r0, Float_t r1, Float_t r2)
3890 //reseed using founded clusters
3892 // Find the number of clusters
3893 Int_t nclusters = 0;
3894 for (Int_t irow=0;irow<160;irow++){
3895 if (track->GetClusterIndex(irow)>0) nclusters++;
3899 ipos[0] = TMath::Max(int(r0*nclusters),0); // point 0 cluster
3900 ipos[1] = TMath::Min(int(r1*nclusters),nclusters-1); //
3901 ipos[2] = TMath::Min(int(r2*nclusters),nclusters-1); // last point
3905 Int_t row[3],sec[3]={0,0,0};
3907 // find track row position at given ratio of the length
3909 for (Int_t irow=0;irow<160;irow++){
3910 if (track->GetClusterIndex2(irow)<0) continue;
3912 for (Int_t ipoint=0;ipoint<3;ipoint++){
3913 if (index<=ipos[ipoint]) row[ipoint] = irow;
3917 //Get cluster and sector position
3918 for (Int_t ipoint=0;ipoint<3;ipoint++){
3919 Int_t clindex = track->GetClusterIndex2(row[ipoint]);
3920 AliTPCclusterMI * cl = GetClusterMI(clindex);
3923 // AliTPCclusterMI * cl = GetClusterMI(clindex);
3926 sec[ipoint] = ((clindex&0xff000000)>>24)%18;
3927 xyz[ipoint][0] = GetXrow(row[ipoint]);
3928 xyz[ipoint][1] = cl->GetY();
3929 xyz[ipoint][2] = cl->GetZ();
3933 // Calculate seed state vector and covariance matrix
3935 Double_t alpha, cs,sn, xx2,yy2;
3937 alpha = (sec[1]-sec[2])*fSectors->GetAlpha();
3938 cs = TMath::Cos(alpha);
3939 sn = TMath::Sin(alpha);
3940 xx2= xyz[1][0]*cs-xyz[1][1]*sn;
3941 yy2= xyz[1][0]*sn+xyz[1][1]*cs;
3945 alpha = (sec[0]-sec[2])*fSectors->GetAlpha();
3946 cs = TMath::Cos(alpha);
3947 sn = TMath::Sin(alpha);
3948 xx2= xyz[0][0]*cs-xyz[0][1]*sn;
3949 yy2= xyz[0][0]*sn+xyz[0][1]*cs;
3955 Double_t x[5],c[15];
3959 x[4]=F1(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
3960 x[2]=F2(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
3961 x[3]=F3n(xyz[2][0],xyz[2][1],xyz[0][0],xyz[0][1],xyz[2][2],xyz[0][2],x[4]);
3963 Double_t sy =0.1, sz =0.1;
3965 Double_t sy1=0.2, sz1=0.2;
3966 Double_t sy2=0.2, sz2=0.2;
3969 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;
3970 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;
3971 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;
3972 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;
3973 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;
3974 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;
3976 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;
3977 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;
3978 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;
3979 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;
3984 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3985 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3986 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3987 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3988 c[13]=f30*sy1*f40+f32*sy2*f42;
3989 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3991 // Int_t row1 = fSectors->GetRowNumber(xyz[2][0]);
3992 AliTPCseed *seed=new AliTPCseed(xyz[2][0], sec[2]*fSectors->GetAlpha()+fSectors->GetAlphaShift(), x, c, 0);
3993 seed->SetLastPoint(row[2]);
3994 seed->SetFirstPoint(row[2]);
3999 AliTPCseed *AliTPCtrackerMI::ReSeed(AliTPCseed *track,Int_t r0, Bool_t forward)
4003 //reseed using founded clusters
4006 Int_t row[3]={0,0,0};
4007 Int_t sec[3]={0,0,0};
4009 // forward direction
4011 for (Int_t irow=r0;irow<160;irow++){
4012 if (track->GetClusterIndex(irow)>0){
4017 for (Int_t irow=160;irow>r0;irow--){
4018 if (track->GetClusterIndex(irow)>0){
4023 for (Int_t irow=row[2]-15;irow>row[0];irow--){
4024 if (track->GetClusterIndex(irow)>0){
4032 for (Int_t irow=0;irow<r0;irow++){
4033 if (track->GetClusterIndex(irow)>0){
4038 for (Int_t irow=r0;irow>0;irow--){
4039 if (track->GetClusterIndex(irow)>0){
4044 for (Int_t irow=row[2]-15;irow>row[0];irow--){
4045 if (track->GetClusterIndex(irow)>0){
4052 if ((row[2]-row[0])<20) return 0;
4053 if (row[1]==0) return 0;
4056 //Get cluster and sector position
4057 for (Int_t ipoint=0;ipoint<3;ipoint++){
4058 Int_t clindex = track->GetClusterIndex2(row[ipoint]);
4059 AliTPCclusterMI * cl = GetClusterMI(clindex);
4062 // AliTPCclusterMI * cl = GetClusterMI(clindex);
4065 sec[ipoint] = ((clindex&0xff000000)>>24)%18;
4066 xyz[ipoint][0] = GetXrow(row[ipoint]);
4067 AliTPCTrackerPoint * point = track->GetTrackPoint(row[ipoint]);
4068 if (point&&ipoint<2){
4070 xyz[ipoint][1] = point->GetY();
4071 xyz[ipoint][2] = point->GetZ();
4074 xyz[ipoint][1] = cl->GetY();
4075 xyz[ipoint][2] = cl->GetZ();
4082 // Calculate seed state vector and covariance matrix
4084 Double_t alpha, cs,sn, xx2,yy2;
4086 alpha = (sec[1]-sec[2])*fSectors->GetAlpha();
4087 cs = TMath::Cos(alpha);
4088 sn = TMath::Sin(alpha);
4089 xx2= xyz[1][0]*cs-xyz[1][1]*sn;
4090 yy2= xyz[1][0]*sn+xyz[1][1]*cs;
4094 alpha = (sec[0]-sec[2])*fSectors->GetAlpha();
4095 cs = TMath::Cos(alpha);
4096 sn = TMath::Sin(alpha);
4097 xx2= xyz[0][0]*cs-xyz[0][1]*sn;
4098 yy2= xyz[0][0]*sn+xyz[0][1]*cs;
4104 Double_t x[5],c[15];
4108 x[4]=F1(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
4109 x[2]=F2(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
4110 x[3]=F3n(xyz[2][0],xyz[2][1],xyz[0][0],xyz[0][1],xyz[2][2],xyz[0][2],x[4]);
4112 Double_t sy =0.1, sz =0.1;
4114 Double_t sy1=0.2, sz1=0.2;
4115 Double_t sy2=0.2, sz2=0.2;
4118 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;
4119 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;
4120 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;
4121 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;
4122 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;
4123 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;
4125 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;
4126 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;
4127 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;
4128 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;
4133 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
4134 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
4135 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
4136 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
4137 c[13]=f30*sy1*f40+f32*sy2*f42;
4138 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
4140 // Int_t row1 = fSectors->GetRowNumber(xyz[2][0]);
4141 AliTPCseed *seed=new AliTPCseed(xyz[2][0], sec[2]*fSectors->GetAlpha()+fSectors->GetAlphaShift(), x, c, 0);
4142 seed->SetLastPoint(row[2]);
4143 seed->SetFirstPoint(row[2]);
4144 for (Int_t i=row[0];i<row[2];i++){
4145 seed->SetClusterIndex(i, track->GetClusterIndex(i));
4153 void AliTPCtrackerMI::FindMultiMC(TObjArray * array, AliESDEvent */*esd*/, Int_t iter)
4156 // find multi tracks - THIS FUNCTION IS ONLY FOR DEBUG PURPOSES
4158 // Use AliTPCReconstructor::StreamLevel()>2 if you want to tune parameters - cuts
4160 // Two reasons to have multiple find tracks
4161 // 1. Curling tracks can be find more than once
4162 // 2. Splitted tracks
4163 // a.) Multiple seeding to increase tracking efficiency - (~ 100% reached)
4164 // b.) Edge effect on the sector boundaries
4167 // Algorithm done in 2 phases - because of CPU consumption
4168 // it is n^2 algorithm - for lead-lead 20000x20000 combination are investigated
4170 // Algorihm for curling tracks sign:
4171 // 1 phase -makes a very rough fast cuts to minimize combinatorics
4172 // a.) opposite sign
4173 // b.) one of the tracks - not pointing to the primary vertex -
4174 // c.) delta tan(theta)
4176 // 2 phase - calculates DCA between tracks - time consument
4181 // General cuts - for splitted tracks and for curling tracks
4183 const Float_t kMaxdPhi = 0.2; // maximal distance in phi
4185 // Curling tracks cuts
4190 Int_t nentries = array->GetEntriesFast();
4191 AliHelix *helixes = new AliHelix[nentries];
4192 Float_t *xm = new Float_t[nentries];
4193 Float_t *dz0 = new Float_t[nentries];
4194 Float_t *dz1 = new Float_t[nentries];
4200 // Find track COG in x direction - point with best defined parameters
4202 for (Int_t i=0;i<nentries;i++){
4203 AliTPCseed* track = (AliTPCseed*)array->At(i);
4204 if (!track) continue;
4205 track->SetCircular(0);
4206 new (&helixes[i]) AliHelix(*track);
4210 track->GetDZ(GetX(),GetY(),GetZ(),GetBz(),dz);
4213 for (Int_t icl=0; icl<160; icl++){
4214 AliTPCclusterMI * cl = track->GetClusterPointer(icl);
4220 if (ncl>0) xm[i]/=Float_t(ncl);
4223 for (Int_t i0=0;i0<nentries;i0++){
4224 AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
4225 if (!track0) continue;
4226 Float_t xc0 = helixes[i0].GetHelix(6);
4227 Float_t yc0 = helixes[i0].GetHelix(7);
4228 Float_t r0 = helixes[i0].GetHelix(8);
4229 Float_t rc0 = TMath::Sqrt(xc0*xc0+yc0*yc0);
4230 Float_t fi0 = TMath::ATan2(yc0,xc0);
4232 for (Int_t i1=i0+1;i1<nentries;i1++){
4233 AliTPCseed * track1 = (AliTPCseed*)array->At(i1);
4234 if (!track1) continue;
4235 Int_t lab0=track0->GetLabel();
4236 Int_t lab1=track1->GetLabel();
4237 if (TMath::Abs(lab0)!=TMath::Abs(lab1)) continue;
4239 Float_t xc1 = helixes[i1].GetHelix(6);
4240 Float_t yc1 = helixes[i1].GetHelix(7);
4241 Float_t r1 = helixes[i1].GetHelix(8);
4242 Float_t rc1 = TMath::Sqrt(xc1*xc1+yc1*yc1);
4243 Float_t fi1 = TMath::ATan2(yc1,xc1);
4245 Float_t dfi = fi0-fi1;
4248 if (dfi>1.5*TMath::Pi()) dfi-=TMath::Pi(); // take care about edge effect
4249 if (dfi<-1.5*TMath::Pi()) dfi+=TMath::Pi(); //
4250 if (TMath::Abs(dfi)>kMaxdPhi&&helixes[i0].GetHelix(4)*helixes[i1].GetHelix(4)<0){
4252 // if short tracks with undefined sign
4253 fi1 = -TMath::ATan2(yc1,-xc1);
4256 Float_t dtheta = TMath::Abs(track0->GetTgl()-track1->GetTgl())<TMath::Abs(track0->GetTgl()+track1->GetTgl())? track0->GetTgl()-track1->GetTgl():track0->GetTgl()+track1->GetTgl();
4259 // debug stream to tune "fast cuts"
4261 Double_t dist[3]; // distance at X
4262 Double_t mdist[3]={0,0,0}; // mean distance X+-40cm
4263 track0->GetDistance(track1,0.5*(xm[i0]+xm[i1])-40.,dist,AliTracker::GetBz());
4264 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[i]);
4265 track0->GetDistance(track1,0.5*(xm[i0]+xm[i1])+40.,dist,AliTracker::GetBz());
4266 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[i]);
4267 track0->GetDistance(track1,0.5*(xm[i0]+xm[i1]),dist,AliTracker::GetBz());
4268 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[i]);
4269 for (Int_t i=0;i<3;i++) mdist[i]*=0.33333;
4273 for (Int_t icl=0; icl<160; icl++){
4274 AliTPCclusterMI * cl0 = track0->GetClusterPointer(icl);
4275 AliTPCclusterMI * cl1 = track1->GetClusterPointer(icl);
4278 if (cl0==cl1) sums++;
4282 if (AliTPCReconstructor::StreamLevel()>0) {
4283 TTreeSRedirector &cstream = *fDebugStreamer;
4288 "Tr0.="<<track0<< // seed0
4289 "Tr1.="<<track1<< // seed1
4290 "h0.="<<&helixes[i0]<<
4291 "h1.="<<&helixes[i1]<<
4293 "sum="<<sum<< //the sum of rows with cl in both
4294 "sums="<<sums<< //the sum of shared clusters
4295 "xm0="<<xm[i0]<< // the center of track
4296 "xm1="<<xm[i1]<< // the x center of track
4297 // General cut variables
4298 "dfi="<<dfi<< // distance in fi angle
4299 "dtheta="<<dtheta<< // distance int theta angle
4305 "dist0="<<dist[0]<< //distance x
4306 "dist1="<<dist[1]<< //distance y
4307 "dist2="<<dist[2]<< //distance z
4308 "mdist0="<<mdist[0]<< //distance x
4309 "mdist1="<<mdist[1]<< //distance y
4310 "mdist2="<<mdist[2]<< //distance z
4324 if (AliTPCReconstructor::StreamLevel()>1) {
4325 AliInfo("Time for curling tracks removal DEBUGGING MC");
4331 void AliTPCtrackerMI::FindSplitted(TObjArray * array, AliESDEvent */*esd*/, Int_t iter)
4335 // Two reasons to have multiple find tracks
4336 // 1. Curling tracks can be find more than once
4337 // 2. Splitted tracks
4338 // a.) Multiple seeding to increase tracking efficiency - (~ 100% reached)
4339 // b.) Edge effect on the sector boundaries
4341 // This function tries to find tracks closed in the parametric space
4343 // cut logic if distance is bigger than cut continue - Do Nothing
4344 const Float_t kMaxdTheta = 0.05; // maximal distance in theta
4345 const Float_t kMaxdPhi = 0.05; // maximal deistance in phi
4346 const Float_t kdelta = 40.; //delta r to calculate track distance
4348 // const Float_t kMaxDist0 = 1.; // maximal distance 0
4349 //const Float_t kMaxDist1 = 0.3; // maximal distance 1 - cut if track in separate rows
4352 TCut csec("csec","abs(Tr0.fRelativeSector-Tr1.fRelativeSector)<2");
4353 TCut cdtheta("cdtheta","abs(dtheta)<0.05");
4358 Int_t nentries = array->GetEntriesFast();
4359 AliHelix *helixes = new AliHelix[nentries];
4360 Float_t *xm = new Float_t[nentries];
4366 //Sort tracks according quality
4368 Int_t nseed = array->GetEntriesFast();
4369 Float_t * quality = new Float_t[nseed];
4370 Int_t * indexes = new Int_t[nseed];
4371 for (Int_t i=0; i<nseed; i++) {
4372 AliTPCseed *pt=(AliTPCseed*)array->UncheckedAt(i);
4377 pt->UpdatePoints(); //select first last max dens points
4378 Float_t * points = pt->GetPoints();
4379 if (points[3]<0.8) quality[i] =-1;
4380 quality[i] = (points[2]-points[0])+pt->GetNumberOfClusters();
4381 //prefer high momenta tracks if overlaps
4382 quality[i] *= TMath::Sqrt(TMath::Abs(pt->Pt())+0.5);
4384 TMath::Sort(nseed,quality,indexes);
4388 // Find track COG in x direction - point with best defined parameters
4390 for (Int_t i=0;i<nentries;i++){
4391 AliTPCseed* track = (AliTPCseed*)array->At(i);
4392 if (!track) continue;
4393 track->SetCircular(0);
4394 new (&helixes[i]) AliHelix(*track);
4397 for (Int_t icl=0; icl<160; icl++){
4398 AliTPCclusterMI * cl = track->GetClusterPointer(icl);
4404 if (ncl>0) xm[i]/=Float_t(ncl);
4407 for (Int_t is0=0;is0<nentries;is0++){
4408 Int_t i0 = indexes[is0];
4409 AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
4410 if (!track0) continue;
4411 if (track0->GetKinkIndexes()[0]!=0) continue;
4412 Float_t xc0 = helixes[i0].GetHelix(6);
4413 Float_t yc0 = helixes[i0].GetHelix(7);
4414 Float_t fi0 = TMath::ATan2(yc0,xc0);
4416 for (Int_t is1=is0+1;is1<nentries;is1++){
4417 Int_t i1 = indexes[is1];
4418 AliTPCseed * track1 = (AliTPCseed*)array->At(i1);
4419 if (!track1) continue;
4421 if (TMath::Abs(track0->GetRelativeSector()-track1->GetRelativeSector())>1) continue;
4422 if (track1->GetKinkIndexes()[0]>0 &&track0->GetKinkIndexes()[0]<0) continue;
4423 if (track1->GetKinkIndexes()[0]!=0) continue;
4425 Float_t dtheta = TMath::Abs(track0->GetTgl()-track1->GetTgl())<TMath::Abs(track0->GetTgl()+track1->GetTgl())? track0->GetTgl()-track1->GetTgl():track0->GetTgl()+track1->GetTgl();
4426 if (TMath::Abs(dtheta)>kMaxdTheta) continue;
4428 Float_t xc1 = helixes[i1].GetHelix(6);
4429 Float_t yc1 = helixes[i1].GetHelix(7);
4430 Float_t fi1 = TMath::ATan2(yc1,xc1);
4432 Float_t dfi = fi0-fi1;
4433 if (dfi>1.5*TMath::Pi()) dfi-=TMath::Pi(); // take care about edge effect
4434 if (dfi<-1.5*TMath::Pi()) dfi+=TMath::Pi(); //
4435 if (TMath::Abs(dfi)>kMaxdPhi&&helixes[i0].GetHelix(4)*helixes[i1].GetHelix(4)<0){
4437 // if short tracks with undefined sign
4438 fi1 = -TMath::ATan2(yc1,-xc1);
4441 if (TMath::Abs(dfi)>kMaxdPhi) continue;
4448 for (Int_t icl=0; icl<160; icl++){
4449 Int_t index0=track0->GetClusterIndex2(icl);
4450 Int_t index1=track1->GetClusterIndex2(icl);
4451 Bool_t used0 = (index0>0 && !(index0&0x8000));
4452 Bool_t used1 = (index1>0 && !(index1&0x8000));
4454 if (used0) sum0++; // used cluster0
4455 if (used1) sum1++; // used clusters1
4456 if (used0&&used1) sum++;
4457 if (index0==index1 && used0 && used1) sums++;
4461 if (sums<10) continue;
4462 if (sum<40) continue;
4463 if (sums/Float_t(TMath::Min(sum0,sum1))<0.5) continue;
4465 Double_t dist[5][4]; // distance at X
4466 Double_t mdist[4]={0,0,0,0}; // mean distance on range +-delta
4470 track0->GetDistance(track1,xm[i0],dist[0],AliTracker::GetBz());
4471 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[0][i]);
4472 track0->GetDistance(track1,xm[i1],dist[1],AliTracker::GetBz());
4473 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[1][i]);
4475 track0->GetDistance(track1,TMath::Min(xm[i1],xm[i0])-kdelta,dist[2],AliTracker::GetBz());
4476 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[2][i]);
4477 track0->GetDistance(track1,TMath::Max(xm[i1],xm[i0])+kdelta,dist[3],AliTracker::GetBz());
4478 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[3][i]);
4480 track0->GetDistance(track1,(xm[i1]+xm[i0])*0.5,dist[4],AliTracker::GetBz());
4481 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[4][i]);
4482 for (Int_t i=0;i<3;i++) mdist[i]*=0.2;
4485 Int_t lab0=track0->GetLabel();
4486 Int_t lab1=track1->GetLabel();
4487 if( AliTPCReconstructor::StreamLevel()>5){
4488 TTreeSRedirector &cstream = *fDebugStreamer;
4489 cstream<<"Splitted"<<
4493 "Tr0.="<<track0<< // seed0
4494 "Tr1.="<<track1<< // seed1
4495 "h0.="<<&helixes[i0]<<
4496 "h1.="<<&helixes[i1]<<
4498 "sum="<<sum<< //the sum of rows with cl in both
4499 "sum0="<<sum0<< //the sum of rows with cl in 0 track
4500 "sum1="<<sum1<< //the sum of rows with cl in 1 track
4501 "sums="<<sums<< //the sum of shared clusters
4502 "xm0="<<xm[i0]<< // the center of track
4503 "xm1="<<xm[i1]<< // the x center of track
4504 // General cut variables
4505 "dfi="<<dfi<< // distance in fi angle
4506 "dtheta="<<dtheta<< // distance int theta angle
4509 "dist0="<<dist[4][0]<< //distance x
4510 "dist1="<<dist[4][1]<< //distance y
4511 "dist2="<<dist[4][2]<< //distance z
4512 "mdist0="<<mdist[0]<< //distance x
4513 "mdist1="<<mdist[1]<< //distance y
4514 "mdist2="<<mdist[2]<< //distance z
4517 delete array->RemoveAt(i1);
4524 AliInfo("Time for splitted tracks removal");
4530 void AliTPCtrackerMI::FindCurling(TObjArray * array, AliESDEvent */*esd*/, Int_t iter)
4533 // find Curling tracks
4534 // Use AliTPCReconstructor::StreamLevel()>1 if you want to tune parameters - cuts
4537 // Algorithm done in 2 phases - because of CPU consumption
4538 // it is n^2 algorithm - for lead-lead 20000x20000 combination are investigated
4539 // see detal in MC part what can be used to cut
4543 const Float_t kMaxC = 400; // maximal curvature to of the track
4544 const Float_t kMaxdTheta = 0.15; // maximal distance in theta
4545 const Float_t kMaxdPhi = 0.15; // maximal distance in phi
4546 const Float_t kPtRatio = 0.3; // ratio between pt
4547 const Float_t kMinDCAR = 2.; // distance to the primary vertex in r - see cpipe cut
4550 // Curling tracks cuts
4553 const Float_t kMaxDeltaRMax = 40; // distance in outer radius
4554 const Float_t kMaxDeltaRMin = 5.; // distance in lower radius - see cpipe cut
4555 const Float_t kMinAngle = 2.9; // angle between tracks
4556 const Float_t kMaxDist = 5; // biggest distance
4558 // The cuts can be tuned using the "MC information stored in Multi tree ==> see FindMultiMC
4561 TCut csign("csign","Tr0.fP[4]*Tr1.fP[4]<0"); //opposite sign
4562 TCut cmax("cmax","abs(Tr0.GetC())>1/400");
4563 TCut cda("cda","sqrt(dtheta^2+dfi^2)<0.15");
4564 TCut ccratio("ccratio","abs((Tr0.fP[4]+Tr1.fP[4])/(abs(Tr0.fP[4])+abs(Tr1.fP[4])))<0.3");
4565 TCut cpipe("cpipe", "min(abs(r0-rc0),abs(r1-rc1))>5");
4567 TCut cdrmax("cdrmax","abs(abs(rc0+r0)-abs(rc1+r1))<40")
4568 TCut cdrmin("cdrmin","abs(abs(rc0+r0)-abs(rc1+r1))<10")
4570 Multi->Draw("dfi","iter==0"+csign+cmax+cda+ccratio); ~94% of curling tracks fulfill
4571 Multi->Draw("min(abs(r0-rc0),abs(r1-rc1))","iter==0&&abs(lab1)==abs(lab0)"+csign+cmax+cda+ccratio+cpipe+cdrmin+cdrmax); //80%
4573 Curling2->Draw("dfi","iter==0&&abs(lab0)==abs(lab1)"+csign+cmax+cdtheta+cdfi+ccratio)
4579 Int_t nentries = array->GetEntriesFast();
4580 AliHelix *helixes = new AliHelix[nentries];
4581 for (Int_t i=0;i<nentries;i++){
4582 AliTPCseed* track = (AliTPCseed*)array->At(i);
4583 if (!track) continue;
4584 track->SetCircular(0);
4585 new (&helixes[i]) AliHelix(*track);
4591 Double_t phase[2][2],radius[2];
4596 for (Int_t i0=0;i0<nentries;i0++){
4597 AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
4598 if (!track0) continue;
4599 if (TMath::Abs(track0->GetC())<1/kMaxC) continue;
4600 Float_t xc0 = helixes[i0].GetHelix(6);
4601 Float_t yc0 = helixes[i0].GetHelix(7);
4602 Float_t r0 = helixes[i0].GetHelix(8);
4603 Float_t rc0 = TMath::Sqrt(xc0*xc0+yc0*yc0);
4604 Float_t fi0 = TMath::ATan2(yc0,xc0);
4606 for (Int_t i1=i0+1;i1<nentries;i1++){
4607 AliTPCseed * track1 = (AliTPCseed*)array->At(i1);
4608 if (!track1) continue;
4609 if (TMath::Abs(track1->GetC())<1/kMaxC) continue;
4610 Float_t xc1 = helixes[i1].GetHelix(6);
4611 Float_t yc1 = helixes[i1].GetHelix(7);
4612 Float_t r1 = helixes[i1].GetHelix(8);
4613 Float_t rc1 = TMath::Sqrt(xc1*xc1+yc1*yc1);
4614 Float_t fi1 = TMath::ATan2(yc1,xc1);
4616 Float_t dfi = fi0-fi1;
4619 if (dfi>1.5*TMath::Pi()) dfi-=TMath::Pi(); // take care about edge effect
4620 if (dfi<-1.5*TMath::Pi()) dfi+=TMath::Pi(); //
4621 Float_t dtheta = TMath::Abs(track0->GetTgl()-track1->GetTgl())<TMath::Abs(track0->GetTgl()+track1->GetTgl())? track0->GetTgl()-track1->GetTgl():track0->GetTgl()+track1->GetTgl();
4625 if (track0->GetBConstrain()&&track1->GetBConstrain()) continue; // not constrained
4626 if (track1->GetSigned1Pt()*track0->GetSigned1Pt()>0) continue; // not the same sign
4627 if ( TMath::Abs(track1->GetTgl()+track0->GetTgl())>kMaxdTheta) continue; //distance in the Theta
4628 if ( TMath::Abs(dfi)>kMaxdPhi) continue; //distance in phi
4629 if ( TMath::Sqrt(dfi*dfi+dtheta*dtheta)>kMaxdPhi) continue; //common angular offset
4631 Float_t pt0 = track0->GetSignedPt();
4632 Float_t pt1 = track1->GetSignedPt();
4633 if ((TMath::Abs(pt0+pt1)/(TMath::Abs(pt0)+TMath::Abs(pt1)))>kPtRatio) continue;
4634 if ((iter==1) && TMath::Abs(TMath::Abs(rc0+r0)-TMath::Abs(rc1+r1))>kMaxDeltaRMax) continue;
4635 if ((iter!=1) &&TMath::Abs(TMath::Abs(rc0-r0)-TMath::Abs(rc1-r1))>kMaxDeltaRMin) continue;
4636 if (TMath::Min(TMath::Abs(rc0-r0),TMath::Abs(rc1-r1))<kMinDCAR) continue;
4639 // Now find closest approach
4643 Int_t npoints = helixes[i0].GetRPHIintersections(helixes[i1], phase, radius,10);
4644 if (npoints==0) continue;
4645 helixes[i0].GetClosestPhases(helixes[i1], phase);
4649 Double_t hangles[3];
4650 helixes[i0].Evaluate(phase[0][0],xyz0);
4651 helixes[i1].Evaluate(phase[0][1],xyz1);
4653 helixes[i0].GetAngle(phase[0][0],helixes[i1],phase[0][1],hangles);
4654 Double_t deltah[2],deltabest;
4655 if (TMath::Abs(hangles[2])<kMinAngle) continue;
4659 helixes[i0].ParabolicDCA(helixes[i1],phase[0][0],phase[0][1],radius[0],deltah[0],2);
4661 helixes[i0].ParabolicDCA(helixes[i1],phase[1][0],phase[1][1],radius[1],deltah[1],2);
4662 if (deltah[1]<deltah[0]) ibest=1;
4664 deltabest = TMath::Sqrt(deltah[ibest]);
4665 helixes[i0].Evaluate(phase[ibest][0],xyz0);
4666 helixes[i1].Evaluate(phase[ibest][1],xyz1);
4667 helixes[i0].GetAngle(phase[ibest][0],helixes[i1],phase[ibest][1],hangles);
4668 Double_t radiusbest = TMath::Sqrt(radius[ibest]);
4670 if (deltabest>kMaxDist) continue;
4671 // if (mindcar+mindcaz<40 && (TMath::Abs(hangles[2])<kMinAngle ||deltabest>3)) continue;
4672 Bool_t sign =kFALSE;
4673 if (hangles[2]>kMinAngle) sign =kTRUE;
4676 // circular[i0] = kTRUE;
4677 // circular[i1] = kTRUE;
4678 if (track0->OneOverPt()<track1->OneOverPt()){
4679 track0->SetCircular(track0->GetCircular()+1);
4680 track1->SetCircular(track1->GetCircular()+2);
4683 track1->SetCircular(track1->GetCircular()+1);
4684 track0->SetCircular(track0->GetCircular()+2);
4687 if (AliTPCReconstructor::StreamLevel()>1){
4689 //debug stream to tune "fine" cuts
4690 Int_t lab0=track0->GetLabel();
4691 Int_t lab1=track1->GetLabel();
4692 TTreeSRedirector &cstream = *fDebugStreamer;
4693 cstream<<"Curling2"<<
4709 "npoints="<<npoints<<
4710 "hangles0="<<hangles[0]<<
4711 "hangles1="<<hangles[1]<<
4712 "hangles2="<<hangles[2]<<
4715 "radius="<<radiusbest<<
4716 "deltabest="<<deltabest<<
4717 "phase0="<<phase[ibest][0]<<
4718 "phase1="<<phase[ibest][1]<<
4726 if (AliTPCReconstructor::StreamLevel()>1) {
4727 AliInfo("Time for curling tracks removal");
4736 void AliTPCtrackerMI::FindKinks(TObjArray * array, AliESDEvent *esd)
4743 TObjArray *kinks= new TObjArray(10000);
4744 // TObjArray *v0s= new TObjArray(10000);
4745 Int_t nentries = array->GetEntriesFast();
4746 AliHelix *helixes = new AliHelix[nentries];
4747 Int_t *sign = new Int_t[nentries];
4748 Int_t *nclusters = new Int_t[nentries];
4749 Float_t *alpha = new Float_t[nentries];
4750 AliKink *kink = new AliKink();
4751 Int_t * usage = new Int_t[nentries];
4752 Float_t *zm = new Float_t[nentries];
4753 Float_t *z0 = new Float_t[nentries];
4754 Float_t *fim = new Float_t[nentries];
4755 Float_t *shared = new Float_t[nentries];
4756 Bool_t *circular = new Bool_t[nentries];
4757 Float_t *dca = new Float_t[nentries];
4758 //const AliESDVertex * primvertex = esd->GetVertex();
4760 // nentries = array->GetEntriesFast();
4765 for (Int_t i=0;i<nentries;i++){
4768 AliTPCseed* track = (AliTPCseed*)array->At(i);
4769 if (!track) continue;
4770 track->SetCircular(0);
4772 track->UpdatePoints();
4773 if (( track->GetPoints()[2]- track->GetPoints()[0])>5 && track->GetPoints()[3]>0.8){
4775 nclusters[i]=track->GetNumberOfClusters();
4776 alpha[i] = track->GetAlpha();
4777 new (&helixes[i]) AliHelix(*track);
4779 helixes[i].Evaluate(0,xyz);
4780 sign[i] = (track->GetC()>0) ? -1:1;
4783 if (track->GetProlongation(x,y,z)){
4785 fim[i] = alpha[i]+TMath::ATan2(y,x);
4788 zm[i] = track->GetZ();
4792 circular[i]= kFALSE;
4793 if (track->GetProlongation(0,y,z)) z0[i] = z;
4794 dca[i] = track->GetD(0,0);
4800 Int_t ncandidates =0;
4803 Double_t phase[2][2],radius[2];
4806 // Find circling track
4808 for (Int_t i0=0;i0<nentries;i0++){
4809 AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
4810 if (!track0) continue;
4811 if (track0->GetNumberOfClusters()<40) continue;
4812 if (TMath::Abs(1./track0->GetC())>200) continue;
4813 for (Int_t i1=i0+1;i1<nentries;i1++){
4814 AliTPCseed * track1 = (AliTPCseed*)array->At(i1);
4815 if (!track1) continue;
4816 if (track1->GetNumberOfClusters()<40) continue;
4817 if ( TMath::Abs(track1->GetTgl()+track0->GetTgl())>0.1) continue;
4818 if (track0->GetBConstrain()&&track1->GetBConstrain()) continue;
4819 if (TMath::Abs(1./track1->GetC())>200) continue;
4820 if (track1->GetSigned1Pt()*track0->GetSigned1Pt()>0) continue;
4821 if (track1->GetTgl()*track0->GetTgl()>0) continue;
4822 if (TMath::Max(TMath::Abs(1./track0->GetC()),TMath::Abs(1./track1->GetC()))>190) continue;
4823 if (track0->GetBConstrain()&&track1->OneOverPt()<track0->OneOverPt()) continue; //returning - lower momenta
4824 if (track1->GetBConstrain()&&track0->OneOverPt()<track1->OneOverPt()) continue; //returning - lower momenta
4826 Float_t mindcar = TMath::Min(TMath::Abs(dca[i0]),TMath::Abs(dca[i1]));
4827 if (mindcar<5) continue;
4828 Float_t mindcaz = TMath::Min(TMath::Abs(z0[i0]-GetZ()),TMath::Abs(z0[i1]-GetZ()));
4829 if (mindcaz<5) continue;
4830 if (mindcar+mindcaz<20) continue;
4833 Float_t xc0 = helixes[i0].GetHelix(6);
4834 Float_t yc0 = helixes[i0].GetHelix(7);
4835 Float_t r0 = helixes[i0].GetHelix(8);
4836 Float_t xc1 = helixes[i1].GetHelix(6);
4837 Float_t yc1 = helixes[i1].GetHelix(7);
4838 Float_t r1 = helixes[i1].GetHelix(8);
4840 Float_t rmean = (r0+r1)*0.5;
4841 Float_t delta =TMath::Sqrt((xc1-xc0)*(xc1-xc0)+(yc1-yc0)*(yc1-yc0));
4842 //if (delta>30) continue;
4843 if (delta>rmean*0.25) continue;
4844 if (TMath::Abs(r0-r1)/rmean>0.3) continue;
4846 Int_t npoints = helixes[i0].GetRPHIintersections(helixes[i1], phase, radius,10);
4847 if (npoints==0) continue;
4848 helixes[i0].GetClosestPhases(helixes[i1], phase);
4852 Double_t hangles[3];
4853 helixes[i0].Evaluate(phase[0][0],xyz0);
4854 helixes[i1].Evaluate(phase[0][1],xyz1);
4856 helixes[i0].GetAngle(phase[0][0],helixes[i1],phase[0][1],hangles);
4857 Double_t deltah[2],deltabest;
4858 if (hangles[2]<2.8) continue;
4861 helixes[i0].ParabolicDCA(helixes[i1],phase[0][0],phase[0][1],radius[0],deltah[0],2);
4863 helixes[i0].ParabolicDCA(helixes[i1],phase[1][0],phase[1][1],radius[1],deltah[1],2);
4864 if (deltah[1]<deltah[0]) ibest=1;
4866 deltabest = TMath::Sqrt(deltah[ibest]);
4867 helixes[i0].Evaluate(phase[ibest][0],xyz0);
4868 helixes[i1].Evaluate(phase[ibest][1],xyz1);
4869 helixes[i0].GetAngle(phase[ibest][0],helixes[i1],phase[ibest][1],hangles);
4870 Double_t radiusbest = TMath::Sqrt(radius[ibest]);
4872 if (deltabest>6) continue;
4873 if (mindcar+mindcaz<40 && (hangles[2]<3.12||deltabest>3)) continue;
4874 Bool_t lsign =kFALSE;
4875 if (hangles[2]>3.06) lsign =kTRUE;
4878 circular[i0] = kTRUE;
4879 circular[i1] = kTRUE;
4880 if (track0->OneOverPt()<track1->OneOverPt()){
4881 track0->SetCircular(track0->GetCircular()+1);
4882 track1->SetCircular(track1->GetCircular()+2);
4885 track1->SetCircular(track1->GetCircular()+1);
4886 track0->SetCircular(track0->GetCircular()+2);
4889 if (lsign&&AliTPCReconstructor::StreamLevel()>1){
4891 Int_t lab0=track0->GetLabel();
4892 Int_t lab1=track1->GetLabel();
4893 TTreeSRedirector &cstream = *fDebugStreamer;
4894 cstream<<"Curling"<<
4901 "mindcar="<<mindcar<<
4902 "mindcaz="<<mindcaz<<
4905 "npoints="<<npoints<<
4906 "hangles0="<<hangles[0]<<
4907 "hangles2="<<hangles[2]<<
4912 "radius="<<radiusbest<<
4913 "deltabest="<<deltabest<<
4914 "phase0="<<phase[ibest][0]<<
4915 "phase1="<<phase[ibest][1]<<
4925 for (Int_t i =0;i<nentries;i++){
4926 if (sign[i]==0) continue;
4927 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
4934 Double_t cradius0 = 40*40;
4935 Double_t cradius1 = 270*270;
4938 Double_t cdist3=0.55;
4939 for (Int_t j =i+1;j<nentries;j++){
4941 if (sign[j]*sign[i]<1) continue;
4942 if ( (nclusters[i]+nclusters[j])>200) continue;
4943 if ( (nclusters[i]+nclusters[j])<80) continue;
4944 if ( TMath::Abs(zm[i]-zm[j])>60.) continue;
4945 if ( TMath::Abs(fim[i]-fim[j])>0.6 && TMath::Abs(fim[i]-fim[j])<5.7 ) continue;
4946 //AliTPCseed * track1 = (AliTPCseed*)array->At(j); Double_t phase[2][2],radius[2];
4947 Int_t npoints = helixes[i].GetRPHIintersections(helixes[j], phase, radius,20);
4948 if (npoints<1) continue;
4951 if (radius[0]<cradius0||radius[0]>cradius1) continue;
4954 if ( (radius[0]<cradius0||radius[0]>cradius1) && (radius[1]<cradius0||radius[1]>cradius1) ) continue;
4957 Double_t delta1=10000,delta2=10000;
4958 // cuts on the intersection radius
4959 helixes[i].LinearDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta1);
4960 if (radius[0]<20&&delta1<1) continue; //intersection at vertex
4961 if (radius[0]<10&&delta1<3) continue; //intersection at vertex
4963 helixes[i].LinearDCA(helixes[j],phase[1][0],phase[1][1],radius[1],delta2);
4964 if (radius[1]<20&&delta2<1) continue; //intersection at vertex
4965 if (radius[1]<10&&delta2<3) continue; //intersection at vertex
4968 Double_t distance1 = TMath::Min(delta1,delta2);
4969 if (distance1>cdist1) continue; // cut on DCA linear approximation
4971 npoints = helixes[i].GetRPHIintersections(helixes[j], phase, radius,20);
4972 helixes[i].ParabolicDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta1);
4973 if (radius[0]<20&&delta1<1) continue; //intersection at vertex
4974 if (radius[0]<10&&delta1<3) continue; //intersection at vertex
4977 helixes[i].ParabolicDCA(helixes[j],phase[1][0],phase[1][1],radius[1],delta2);
4978 if (radius[1]<20&&delta2<1) continue; //intersection at vertex
4979 if (radius[1]<10&&delta2<3) continue; //intersection at vertex
4981 distance1 = TMath::Min(delta1,delta2);
4984 rkink = TMath::Sqrt(radius[0]);
4987 rkink = TMath::Sqrt(radius[1]);
4989 if (distance1>cdist2) continue;
4992 AliTPCseed * track1 = (AliTPCseed*)array->At(j);
4995 Int_t row0 = GetRowNumber(rkink);
4996 if (row0<10) continue;
4997 if (row0>150) continue;
5000 Float_t dens00=-1,dens01=-1;
5001 Float_t dens10=-1,dens11=-1;
5003 Int_t found,foundable,ishared;
5004 track0->GetClusterStatistic(0,row0-5, found, foundable,ishared,kFALSE);
5005 if (foundable>5) dens00 = Float_t(found)/Float_t(foundable);
5006 track0->GetClusterStatistic(row0+5,155, found, foundable,ishared,kFALSE);
5007 if (foundable>5) dens01 = Float_t(found)/Float_t(foundable);
5009 track1->GetClusterStatistic(0,row0-5, found, foundable,ishared,kFALSE);
5010 if (foundable>10) dens10 = Float_t(found)/Float_t(foundable);
5011 track1->GetClusterStatistic(row0+5,155, found, foundable,ishared,kFALSE);
5012 if (foundable>10) dens11 = Float_t(found)/Float_t(foundable);
5014 if (dens00<dens10 && dens01<dens11) continue;
5015 if (dens00>dens10 && dens01>dens11) continue;
5016 if (TMath::Max(dens00,dens10)<0.1) continue;
5017 if (TMath::Max(dens01,dens11)<0.3) continue;
5019 if (TMath::Min(dens00,dens10)>0.6) continue;
5020 if (TMath::Min(dens01,dens11)>0.6) continue;
5023 AliTPCseed * ktrack0, *ktrack1;
5032 if (TMath::Abs(ktrack0->GetC())>5) continue; // cut on the curvature for mother particle
5033 AliExternalTrackParam paramm(*ktrack0);
5034 AliExternalTrackParam paramd(*ktrack1);
5035 if (row0>60&&ktrack1->GetReference().GetX()>90.)new (¶md) AliExternalTrackParam(ktrack1->GetReference());
5038 kink->SetMother(paramm);
5039 kink->SetDaughter(paramd);
5042 Float_t x[3] = { kink->GetPosition()[0],kink->GetPosition()[1],kink->GetPosition()[2]};
5044 fParam->Transform0to1(x,index);
5045 fParam->Transform1to2(x,index);
5046 row0 = GetRowNumber(x[0]);
5048 if (kink->GetR()<100) continue;
5049 if (kink->GetR()>240) continue;
5050 if (kink->GetPosition()[2]/kink->GetR()>AliTPCReconstructor::GetCtgRange()) continue; //out of fiducial volume
5051 if (kink->GetDistance()>cdist3) continue;
5052 Float_t dird = kink->GetDaughterP()[0]*kink->GetPosition()[0]+kink->GetDaughterP()[1]*kink->GetPosition()[1]; // rough direction estimate
5053 if (dird<0) continue;
5055 Float_t dirm = kink->GetMotherP()[0]*kink->GetPosition()[0]+kink->GetMotherP()[1]*kink->GetPosition()[1]; // rough direction estimate
5056 if (dirm<0) continue;
5057 Float_t mpt = TMath::Sqrt(kink->GetMotherP()[0]*kink->GetMotherP()[0]+kink->GetMotherP()[1]*kink->GetMotherP()[1]);
5058 if (mpt<0.2) continue;
5061 //for high momenta momentum not defined well in first iteration
5062 Double_t qt = TMath::Sin(kink->GetAngle(2))*ktrack1->GetP();
5063 if (qt>0.35) continue;
5066 kink->SetLabel(CookLabel(ktrack0,0.4,0,row0),0);
5067 kink->SetLabel(CookLabel(ktrack1,0.4,row0,160),1);
5069 kink->SetTPCDensity(dens00,0,0);
5070 kink->SetTPCDensity(dens01,0,1);
5071 kink->SetTPCDensity(dens10,1,0);
5072 kink->SetTPCDensity(dens11,1,1);
5073 kink->SetIndex(i,0);
5074 kink->SetIndex(j,1);
5077 kink->SetTPCDensity(dens10,0,0);
5078 kink->SetTPCDensity(dens11,0,1);
5079 kink->SetTPCDensity(dens00,1,0);
5080 kink->SetTPCDensity(dens01,1,1);
5081 kink->SetIndex(j,0);
5082 kink->SetIndex(i,1);
5085 if (mpt<1||kink->GetAngle(2)>0.1){
5086 // angle and densities not defined yet
5087 if (kink->GetTPCDensityFactor()<0.8) continue;
5088 if ((2-kink->GetTPCDensityFactor())*kink->GetDistance() >0.25) continue;
5089 if (kink->GetAngle(2)*ktrack0->GetP()<0.003) continue; //too small angle
5090 if (kink->GetAngle(2)>0.2&&kink->GetTPCDensityFactor()<1.15) continue;
5091 if (kink->GetAngle(2)>0.2&&kink->GetTPCDensity(0,1)>0.05) continue;
5093 Float_t criticalangle = track0->GetSigmaSnp2()+track0->GetSigmaTgl2();
5094 criticalangle+= track1->GetSigmaSnp2()+track1->GetSigmaTgl2();
5095 criticalangle= 3*TMath::Sqrt(criticalangle);
5096 if (criticalangle>0.02) criticalangle=0.02;
5097 if (kink->GetAngle(2)<criticalangle) continue;
5100 Int_t drow = Int_t(2.+0.5/(0.05+kink->GetAngle(2))); // overlap region defined
5101 Float_t shapesum =0;
5103 for ( Int_t row = row0-drow; row<row0+drow;row++){
5104 if (row<0) continue;
5105 if (row>155) continue;
5106 if (ktrack0->GetClusterPointer(row)){
5107 AliTPCTrackerPoint *point =ktrack0->GetTrackPoint(row);
5108 shapesum+=point->GetSigmaY()+point->GetSigmaZ();
5111 if (ktrack1->GetClusterPointer(row)){
5112 AliTPCTrackerPoint *point =ktrack1->GetTrackPoint(row);
5113 shapesum+=point->GetSigmaY()+point->GetSigmaZ();
5118 kink->SetShapeFactor(-1.);
5121 kink->SetShapeFactor(shapesum/sum);
5123 // esd->AddKink(kink);
5124 kinks->AddLast(kink);
5130 // sort the kinks according quality - and refit them towards vertex
5132 Int_t nkinks = kinks->GetEntriesFast();
5133 Float_t *quality = new Float_t[nkinks];
5134 Int_t *indexes = new Int_t[nkinks];
5135 AliTPCseed *mothers = new AliTPCseed[nkinks];
5136 AliTPCseed *daughters = new AliTPCseed[nkinks];
5139 for (Int_t i=0;i<nkinks;i++){
5141 AliKink *kinkl = (AliKink*)kinks->At(i);
5143 // refit kinks towards vertex
5145 Int_t index0 = kinkl->GetIndex(0);
5146 Int_t index1 = kinkl->GetIndex(1);
5147 AliTPCseed * ktrack0 = (AliTPCseed*)array->At(index0);
5148 AliTPCseed * ktrack1 = (AliTPCseed*)array->At(index1);
5150 Int_t sumn=ktrack0->GetNumberOfClusters()+ktrack1->GetNumberOfClusters();
5152 // Refit Kink under if too small angle
5154 if (kinkl->GetAngle(2)<0.05){
5155 kinkl->SetTPCRow0(GetRowNumber(kinkl->GetR()));
5156 Int_t row0 = kinkl->GetTPCRow0();
5157 Int_t drow = Int_t(2.+0.5/(0.05+kinkl->GetAngle(2)));
5160 Int_t last = row0-drow;
5161 if (last<40) last=40;
5162 if (last<ktrack0->GetFirstPoint()+25) last = ktrack0->GetFirstPoint()+25;
5163 AliTPCseed* seed0 = ReSeed(ktrack0,last,kFALSE);
5166 Int_t first = row0+drow;
5167 if (first>130) first=130;
5168 if (first>ktrack1->GetLastPoint()-25) first = TMath::Max(ktrack1->GetLastPoint()-25,30);
5169 AliTPCseed* seed1 = ReSeed(ktrack1,first,kTRUE);
5171 if (seed0 && seed1){
5172 kinkl->SetStatus(1,8);
5173 if (RefitKink(*seed0,*seed1,*kinkl)) kinkl->SetStatus(1,9);
5174 row0 = GetRowNumber(kinkl->GetR());
5175 sumn = seed0->GetNumberOfClusters()+seed1->GetNumberOfClusters();
5176 mothers[i] = *seed0;
5177 daughters[i] = *seed1;
5180 delete kinks->RemoveAt(i);
5181 if (seed0) delete seed0;
5182 if (seed1) delete seed1;
5185 if (kinkl->GetDistance()>0.5 || kinkl->GetR()<110 || kinkl->GetR()>240) {
5186 delete kinks->RemoveAt(i);
5187 if (seed0) delete seed0;
5188 if (seed1) delete seed1;
5196 if (kinkl) quality[i] = 160*((0.1+kinkl->GetDistance())*(2.-kinkl->GetTPCDensityFactor()))/(sumn+40.); //the longest -clossest will win
5198 TMath::Sort(nkinks,quality,indexes,kFALSE);
5200 //remove double find kinks
5202 for (Int_t ikink0=1;ikink0<nkinks;ikink0++){
5203 AliKink * kink0 = (AliKink*) kinks->At(indexes[ikink0]);
5204 if (!kink0) continue;
5206 for (Int_t ikink1=0;ikink1<ikink0;ikink1++){
5207 if (!kink0) continue;
5208 AliKink * kink1 = (AliKink*) kinks->At(indexes[ikink1]);
5209 if (!kink1) continue;
5210 // if not close kink continue
5211 if (TMath::Abs(kink1->GetPosition()[2]-kink0->GetPosition()[2])>10) continue;
5212 if (TMath::Abs(kink1->GetPosition()[1]-kink0->GetPosition()[1])>10) continue;
5213 if (TMath::Abs(kink1->GetPosition()[0]-kink0->GetPosition()[0])>10) continue;
5215 AliTPCseed &mother0 = mothers[indexes[ikink0]];
5216 AliTPCseed &daughter0 = daughters[indexes[ikink0]];
5217 AliTPCseed &mother1 = mothers[indexes[ikink1]];
5218 AliTPCseed &daughter1 = daughters[indexes[ikink1]];
5219 Int_t row0 = (kink0->GetTPCRow0()+kink1->GetTPCRow0())/2;
5228 for (Int_t i=0;i<row0;i++){
5229 if (mother0.GetClusterIndex(i)>0 && mother1.GetClusterIndex(i)>0){
5232 if (mother0.GetClusterIndex(i)==mother1.GetClusterIndex(i)){
5239 for (Int_t i=row0;i<158;i++){
5240 if (daughter0.GetClusterIndex(i)>0 && daughter0.GetClusterIndex(i)>0){
5243 if (mother0.GetClusterIndex(i)==mother1.GetClusterIndex(i)){
5249 Float_t ratio = Float_t(same+1)/Float_t(both+1);
5250 Float_t ratiom = Float_t(samem+1)/Float_t(bothm+1);
5251 Float_t ratiod = Float_t(samed+1)/Float_t(bothd+1);
5252 if (ratio>0.3 && ratiom>0.5 &&ratiod>0.5) {
5253 Int_t sum0 = mother0.GetNumberOfClusters()+daughter0.GetNumberOfClusters();
5254 Int_t sum1 = mother1.GetNumberOfClusters()+daughter1.GetNumberOfClusters();
5256 shared[kink0->GetIndex(0)]= kTRUE;
5257 shared[kink0->GetIndex(1)]= kTRUE;
5258 delete kinks->RemoveAt(indexes[ikink0]);
5261 shared[kink1->GetIndex(0)]= kTRUE;
5262 shared[kink1->GetIndex(1)]= kTRUE;
5263 delete kinks->RemoveAt(indexes[ikink1]);
5270 for (Int_t i=0;i<nkinks;i++){
5271 AliKink * kinkl = (AliKink*) kinks->At(indexes[i]);
5272 if (!kinkl) continue;
5273 kinkl->SetTPCRow0(GetRowNumber(kinkl->GetR()));
5274 Int_t index0 = kinkl->GetIndex(0);
5275 Int_t index1 = kinkl->GetIndex(1);
5276 if (circular[index0]||(circular[index1]&&kinkl->GetDistance()>0.2)) continue;
5277 kinkl->SetMultiple(usage[index0],0);
5278 kinkl->SetMultiple(usage[index1],1);
5279 if (kinkl->GetMultiple()[0]+kinkl->GetMultiple()[1]>2) continue;
5280 if (kinkl->GetMultiple()[0]+kinkl->GetMultiple()[1]>0 && quality[indexes[i]]>0.2) continue;
5281 if (kinkl->GetMultiple()[0]+kinkl->GetMultiple()[1]>0 && kinkl->GetDistance()>0.2) continue;
5282 if (circular[index0]||(circular[index1]&&kinkl->GetDistance()>0.1)) continue;
5284 AliTPCseed * ktrack0 = (AliTPCseed*)array->At(index0);
5285 AliTPCseed * ktrack1 = (AliTPCseed*)array->At(index1);
5286 if (!ktrack0 || !ktrack1) continue;
5287 Int_t index = esd->AddKink(kinkl);
5290 if ( ktrack0->GetKinkIndex(0)==0 && ktrack1->GetKinkIndex(0)==0) { //best kink
5291 if (mothers[indexes[i]].GetNumberOfClusters()>20 && daughters[indexes[i]].GetNumberOfClusters()>20 && (mothers[indexes[i]].GetNumberOfClusters()+daughters[indexes[i]].GetNumberOfClusters())>100){
5292 *ktrack0 = mothers[indexes[i]];
5293 *ktrack1 = daughters[indexes[i]];
5297 ktrack0->SetKinkIndex(usage[index0],-(index+1));
5298 ktrack1->SetKinkIndex(usage[index1], (index+1));
5303 // Remove tracks corresponding to shared kink's
5305 for (Int_t i=0;i<nentries;i++){
5306 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5307 if (!track0) continue;
5308 if (track0->GetKinkIndex(0)!=0) continue;
5309 if (shared[i]) delete array->RemoveAt(i);
5314 RemoveUsed2(array,0.5,0.4,30);
5316 for (Int_t i=0;i<nentries;i++){
5317 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5318 if (!track0) continue;
5319 track0->CookdEdx(0.02,0.6);
5323 for (Int_t i=0;i<nentries;i++){
5324 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5325 if (!track0) continue;
5326 if (track0->Pt()<1.4) continue;
5327 //remove double high momenta tracks - overlapped with kink candidates
5330 for (Int_t icl=track0->GetFirstPoint();icl<track0->GetLastPoint(); icl++){
5331 if (track0->GetClusterPointer(icl)!=0){
5333 if (track0->GetClusterPointer(icl)->IsUsed(10)) ishared++;
5336 if (Float_t(ishared+1)/Float_t(all+1)>0.5) {
5337 delete array->RemoveAt(i);
5341 if (track0->GetKinkIndex(0)!=0) continue;
5342 if (track0->GetNumberOfClusters()<80) continue;
5344 AliTPCseed *pmother = new AliTPCseed();
5345 AliTPCseed *pdaughter = new AliTPCseed();
5346 AliKink *pkink = new AliKink;
5348 AliTPCseed & mother = *pmother;
5349 AliTPCseed & daughter = *pdaughter;
5350 AliKink & kinkl = *pkink;
5351 if (CheckKinkPoint(track0,mother,daughter, kinkl)){
5352 if (mother.GetNumberOfClusters()<30||daughter.GetNumberOfClusters()<20) {
5356 continue; //too short tracks
5358 if (mother.Pt()<1.4) {
5364 Int_t row0= kinkl.GetTPCRow0();
5365 if (kinkl.GetDistance()>0.5 || kinkl.GetR()<110. || kinkl.GetR()>240.) {
5372 Int_t index = esd->AddKink(&kinkl);
5373 mother.SetKinkIndex(0,-(index+1));
5374 daughter.SetKinkIndex(0,index+1);
5375 if (mother.GetNumberOfClusters()>50) {
5376 delete array->RemoveAt(i);
5377 array->AddAt(new AliTPCseed(mother),i);
5380 array->AddLast(new AliTPCseed(mother));
5382 array->AddLast(new AliTPCseed(daughter));
5383 for (Int_t icl=0;icl<row0;icl++) {
5384 if (mother.GetClusterPointer(icl)) mother.GetClusterPointer(icl)->Use(20);
5387 for (Int_t icl=row0;icl<158;icl++) {
5388 if (daughter.GetClusterPointer(icl)) daughter.GetClusterPointer(icl)->Use(20);
5397 delete [] daughters;
5419 printf("Ncandidates=\t%d\t%d\t%d\t%d\n",esd->GetNumberOfKinks(),ncandidates,ntracks,nall);
5423 void AliTPCtrackerMI::FindV0s(TObjArray * array, AliESDEvent *esd)
5429 TObjArray *tpcv0s = new TObjArray(100000);
5430 Int_t nentries = array->GetEntriesFast();
5431 AliHelix *helixes = new AliHelix[nentries];
5432 Int_t *sign = new Int_t[nentries];
5433 Float_t *alpha = new Float_t[nentries];
5434 Float_t *z0 = new Float_t[nentries];
5435 Float_t *dca = new Float_t[nentries];
5436 Float_t *sdcar = new Float_t[nentries];
5437 Float_t *cdcar = new Float_t[nentries];
5438 Float_t *pulldcar = new Float_t[nentries];
5439 Float_t *pulldcaz = new Float_t[nentries];
5440 Float_t *pulldca = new Float_t[nentries];
5441 Bool_t *isPrim = new Bool_t[nentries];
5442 const AliESDVertex * primvertex = esd->GetVertex();
5443 Double_t zvertex = primvertex->GetZv();
5445 // nentries = array->GetEntriesFast();
5447 for (Int_t i=0;i<nentries;i++){
5450 AliTPCseed* track = (AliTPCseed*)array->At(i);
5451 if (!track) continue;
5452 track->GetV0Indexes()[0] = 0; //rest v0 indexes
5453 track->GetV0Indexes()[1] = 0; //rest v0 indexes
5454 track->GetV0Indexes()[2] = 0; //rest v0 indexes
5456 alpha[i] = track->GetAlpha();
5457 new (&helixes[i]) AliHelix(*track);
5459 helixes[i].Evaluate(0,xyz);
5460 sign[i] = (track->GetC()>0) ? -1:1;
5464 if (track->GetProlongation(0,y,z)) z0[i] = z;
5465 dca[i] = track->GetD(0,0);
5467 // dca error parrameterezation + pulls
5469 sdcar[i] = TMath::Sqrt(0.150*0.150+(100*track->GetC())*(100*track->GetC()));
5470 if (TMath::Abs(track->GetTgl())>1) sdcar[i]*=2.5;
5471 cdcar[i] = TMath::Exp((TMath::Abs(track->GetC())-0.0106)*525.3);
5472 pulldcar[i] = (dca[i]-cdcar[i])/sdcar[i];
5473 pulldcaz[i] = (z0[i]-zvertex)/sdcar[i];
5474 pulldca[i] = TMath::Sqrt(pulldcar[i]*pulldcar[i]+pulldcaz[i]*pulldcaz[i]);
5475 if (track->TPCrPID(1)+track->TPCrPID(2)+track->TPCrPID(3)>0.5) {
5476 if (pulldca[i]<3.) isPrim[i]=kTRUE; //pion, muon and Kaon 3 sigma cut
5478 if (track->TPCrPID(4)>0.5) {
5479 if (pulldca[i]<0.5) isPrim[i]=kTRUE; //proton 0.5 sigma cut
5481 if (track->TPCrPID(0)>0.4) {
5482 isPrim[i]=kFALSE; //electron no sigma cut
5489 Int_t ncandidates =0;
5492 Double_t phase[2][2],radius[2];
5498 Float_t fprimvertex[3]={GetX(),GetY(),GetZ()};
5500 Double_t cradius0 = 10*10;
5501 Double_t cradius1 = 200*200;
5504 Double_t cpointAngle = 0.95;
5506 Double_t delta[2]={10000,10000};
5507 for (Int_t i =0;i<nentries;i++){
5508 if (sign[i]==0) continue;
5509 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5510 if (!track0) continue;
5511 if (AliTPCReconstructor::StreamLevel()>1){
5512 TTreeSRedirector &cstream = *fDebugStreamer;
5517 "zvertex="<<zvertex<<
5518 "sdcar0="<<sdcar[i]<<
5519 "cdcar0="<<cdcar[i]<<
5520 "pulldcar0="<<pulldcar[i]<<
5521 "pulldcaz0="<<pulldcaz[i]<<
5522 "pulldca0="<<pulldca[i]<<
5523 "isPrim="<<isPrim[i]<<
5527 if (track0->GetSigned1Pt()<0) continue;
5528 if (track0->GetKinkIndex(0)>0||isPrim[i]) continue; //daughter kink
5530 if (TMath::Abs(helixes[i].GetHelix(4))<0.000000001) continue;
5535 for (Int_t j =0;j<nentries;j++){
5536 AliTPCseed * track1 = (AliTPCseed*)array->At(j);
5537 if (!track1) continue;
5538 if (track1->GetKinkIndex(0)>0 || isPrim[j]) continue; //daughter kink
5539 if (sign[j]*sign[i]>0) continue;
5540 if (TMath::Abs(helixes[j].GetHelix(4))<0.000001) continue;
5541 if (track0->GetCircular()+track1->GetCircular()>1) continue; //circling -returning track
5544 // DCA to prim vertex cut
5550 Int_t npoints = helixes[i].GetRPHIintersections(helixes[j], phase, radius,cdist2);
5551 if (npoints<1) continue;
5555 if (radius[0]<cradius0||radius[0]>cradius1) continue;
5556 helixes[i].LinearDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta[0]);
5557 if (delta[0]>cdist1) continue;
5560 if (TMath::Max(radius[0],radius[1])<cradius0|| TMath::Min(radius[0],radius[1])>cradius1) continue;
5561 helixes[i].LinearDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta[0]);
5562 helixes[i].LinearDCA(helixes[j],phase[1][0],phase[1][1],radius[1],delta[1]);
5563 if (delta[1]<delta[0]) iclosest=1;
5564 if (delta[iclosest]>cdist1) continue;
5566 helixes[i].ParabolicDCA(helixes[j],phase[iclosest][0],phase[iclosest][1],radius[iclosest],delta[iclosest]);
5567 if (radius[iclosest]<cradius0 || radius[iclosest]>cradius1 || delta[iclosest]>cdist1) continue;
5569 Double_t pointAngle = helixes[i].GetPointAngle(helixes[j],phase[iclosest],fprimvertex);
5570 if (pointAngle<cpointAngle) continue;
5572 Bool_t isGamma = kFALSE;
5573 vertex.SetParamP(*track0); //track0 - plus
5574 vertex.SetParamN(*track1); //track1 - minus
5575 vertex.Update(fprimvertex);
5576 if (track0->TPCrPID(0)>0.3&&track1->TPCrPID(0)>0.3&&vertex.GetAnglep()[2]<0.15) isGamma=kTRUE; // gamma conversion candidate
5577 Double_t pointAngle2 = vertex.GetV0CosineOfPointingAngle();
5579 if (vertex.GetV0CosineOfPointingAngle()<cpointAngle && (!isGamma)) continue;// point angle cut
5580 //Bo: if (vertex.GetDist2()>2&&(!isGamma)) continue; // point angle cut
5581 if (vertex.GetDcaV0Daughters()>2&&(!isGamma)) continue;//Bo: // point angle cut
5582 Float_t sigmae = 0.15*0.15;
5583 if (vertex.GetRr()<80)
5584 sigmae += (sdcar[i]*sdcar[i]+sdcar[j]*sdcar[j])*(1.-vertex.GetRr()/80.)*(1.-vertex.GetRr()/80.);
5585 sigmae+= TMath::Sqrt(sigmae);
5586 //Bo: if (vertex.GetDist2()/sigmae>3.&&(!isGamma)) continue;
5587 if (vertex.GetDcaV0Daughters()/sigmae>3.&&(!isGamma)) continue;
5588 Float_t densb0=0,densb1=0,densa0=0,densa1=0;
5589 Int_t row0 = GetRowNumber(vertex.GetRr());
5591 //Bo: if (vertex.GetDist2()>0.2) continue;
5592 if (vertex.GetDcaV0Daughters()>0.2) continue;
5593 densb0 = track0->Density2(0,row0-5);
5594 densb1 = track1->Density2(0,row0-5);
5595 if (densb0>0.3|| densb1>0.3) continue; //clusters before vertex
5596 densa0 = track0->Density2(row0+5,row0+40);
5597 densa1 = track1->Density2(row0+5,row0+40);
5598 if ((densa0<0.4|| densa1<0.4)&&(!isGamma)) continue; //missing clusters after vertex
5601 densa0 = track0->Density2(0,40); //cluster density
5602 densa1 = track1->Density2(0,40); //cluster density
5603 if ((vertex.GetRr()<80&&densa0+densa1<1.)&&(!isGamma)) continue;
5605 //Bo: vertex.SetLab(0,track0->GetLabel());
5606 //Bo: vertex.SetLab(1,track1->GetLabel());
5607 vertex.SetChi2After((densa0+densa1)*0.5);
5608 vertex.SetChi2Before((densb0+densb1)*0.5);
5609 vertex.SetIndex(0,i);
5610 vertex.SetIndex(1,j);
5611 //Bo: vertex.SetStatus(1); // TPC v0 candidate
5612 vertex.SetOnFlyStatus(2);//Bo: // TPC v0 candidate
5613 //Bo: vertex.SetRp(track0->TPCrPIDs());
5614 //Bo: vertex.SetRm(track1->TPCrPIDs());
5615 tpcv0s->AddLast(new AliESDv0(vertex));
5618 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
5619 Double_t radiusm= (delta[0]<delta[1])? TMath::Sqrt(radius[0]):TMath::Sqrt(radius[1]);
5620 Double_t deltam= (delta[0]<delta[1])? TMath::Sqrt(delta[0]):TMath::Sqrt(delta[1]);
5621 if (AliTPCReconstructor::StreamLevel()>1) {
5622 Int_t lab0=track0->GetLabel();
5623 Int_t lab1=track1->GetLabel();
5624 Char_t c0=track0->GetCircular();
5625 Char_t c1=track1->GetCircular();
5626 TTreeSRedirector &cstream = *fDebugStreamer;
5629 "vertex.="<<&vertex<<
5632 "Helix0.="<<&helixes[i]<<
5635 "Helix1.="<<&helixes[j]<<
5636 "pointAngle="<<pointAngle<<
5637 "pointAngle2="<<pointAngle2<<
5642 "zvertex="<<zvertex<<
5645 "npoints="<<npoints<<
5646 "radius0="<<radius[0]<<
5647 "delta0="<<delta[0]<<
5648 "radius1="<<radius[1]<<
5649 "delta1="<<delta[1]<<
5650 "radiusm="<<radiusm<<
5652 "sdcar0="<<sdcar[i]<<
5653 "sdcar1="<<sdcar[j]<<
5654 "cdcar0="<<cdcar[i]<<
5655 "cdcar1="<<cdcar[j]<<
5656 "pulldcar0="<<pulldcar[i]<<
5657 "pulldcar1="<<pulldcar[j]<<
5658 "pulldcaz0="<<pulldcaz[i]<<
5659 "pulldcaz1="<<pulldcaz[j]<<
5660 "pulldca0="<<pulldca[i]<<
5661 "pulldca1="<<pulldca[j]<<
5671 Float_t *quality = new Float_t[ncandidates];
5672 Int_t *indexes = new Int_t[ncandidates];
5674 for (Int_t i=0;i<ncandidates;i++){
5676 AliESDv0 *v0 = (AliESDv0*)tpcv0s->At(i);
5677 quality[i] = 1./(1.00001-v0->GetV0CosineOfPointingAngle()); //base point angle
5678 // quality[i] /= (0.5+v0->GetDist2());
5679 // quality[i] *= v0->GetChi2After(); //density factor
5681 Int_t index0 = v0->GetIndex(0);
5682 Int_t index1 = v0->GetIndex(1);
5683 //Bo: Double_t minpulldca = TMath::Min(2.+pulldca[v0->GetIndex(0)],(2.+pulldca[v0->GetIndex(1)]) ); //pull
5684 Double_t minpulldca = TMath::Min(2.+pulldca[index0],(2.+pulldca[index1]) );//Bo:
5688 AliTPCseed * track0 = (AliTPCseed*)array->At(index0);
5689 AliTPCseed * track1 = (AliTPCseed*)array->At(index1);
5690 if (track0->TPCrPID(0)>0.3&&track1->TPCrPID(0)>0.3&&v0->GetAnglep()[2]<0.15) quality[i]+=1000000; // gamma conversion candidate
5691 if (track0->TPCrPID(4)>0.9||(track1->TPCrPID(4)>0.9&&minpulldca>4)) quality[i]*=10; // lambda candidate candidate
5694 TMath::Sort(ncandidates,quality,indexes,kTRUE);
5697 for (Int_t i=0;i<ncandidates;i++){
5698 AliESDv0 * v0 = (AliESDv0*)tpcv0s->At(indexes[i]);
5700 Int_t index0 = v0->GetIndex(0);
5701 Int_t index1 = v0->GetIndex(1);
5702 AliTPCseed * track0 = (AliTPCseed*)array->At(index0);
5703 AliTPCseed * track1 = (AliTPCseed*)array->At(index1);
5704 if (!track0||!track1) {
5708 Bool_t accept =kTRUE; //default accept
5709 Int_t *v0indexes0 = track0->GetV0Indexes();
5710 Int_t *v0indexes1 = track1->GetV0Indexes();
5712 Int_t order0 = (v0indexes0[0]!=0) ? 1:0;
5713 Int_t order1 = (v0indexes1[0]!=0) ? 1:0;
5714 if (v0indexes0[1]!=0) order0 =2;
5715 if (v0indexes1[1]!=0) order1 =2;
5717 if (v0indexes0[2]!=0) {order0=3; accept=kFALSE;}
5718 if (v0indexes0[2]!=0) {order1=3; accept=kFALSE;}
5720 AliESDv0 * v02 = v0;
5722 //Bo: v0->SetOrder(0,order0);
5723 //Bo: v0->SetOrder(1,order1);
5724 //Bo: v0->SetOrder(1,order0+order1);
5725 v0->SetOnFlyStatus(kTRUE);
5726 Int_t index = esd->AddV0(v0);
5727 v02 = esd->GetV0(index);
5728 v0indexes0[order0]=index;
5729 v0indexes1[order1]=index;
5733 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
5734 if (AliTPCReconstructor::StreamLevel()>1) {
5735 Int_t lab0=track0->GetLabel();
5736 Int_t lab1=track1->GetLabel();
5737 TTreeSRedirector &cstream = *fDebugStreamer;
5746 "dca0="<<dca[index0]<<
5747 "dca1="<<dca[index1]<<
5751 "quality="<<quality[i]<<
5752 "pulldca0="<<pulldca[index0]<<
5753 "pulldca1="<<pulldca[index1]<<
5777 printf("TPC V0 finder : naccepted\t%d\tncandidates\t%d\tntracks\t%d\tnall\t%d\n",naccepted,ncandidates,ntracks,nall);
5781 Int_t AliTPCtrackerMI::RefitKink(AliTPCseed &mother, AliTPCseed &daughter, AliESDkink &knk)
5784 // refit kink towards to the vertex
5787 AliKink &kink=(AliKink &)knk;
5789 Int_t row0 = GetRowNumber(kink.GetR());
5790 FollowProlongation(mother,0);
5791 mother.Reset(kFALSE);
5793 FollowProlongation(daughter,row0);
5794 daughter.Reset(kFALSE);
5795 FollowBackProlongation(daughter,158);
5796 daughter.Reset(kFALSE);
5797 Int_t first = TMath::Max(row0-20,30);
5798 Int_t last = TMath::Min(row0+20,140);
5800 const Int_t kNdiv =5;
5801 AliTPCseed param0[kNdiv]; // parameters along the track
5802 AliTPCseed param1[kNdiv]; // parameters along the track
5803 AliKink kinks[kNdiv]; // corresponding kink parameters
5806 for (Int_t irow=0; irow<kNdiv;irow++){
5807 rows[irow] = first +((last-first)*irow)/(kNdiv-1);
5809 // store parameters along the track
5811 for (Int_t irow=0;irow<kNdiv;irow++){
5812 FollowBackProlongation(mother, rows[irow]);
5813 FollowProlongation(daughter,rows[kNdiv-1-irow]);
5814 param0[irow] = mother;
5815 param1[kNdiv-1-irow] = daughter;
5819 for (Int_t irow=0; irow<kNdiv-1;irow++){
5820 if (param0[irow].GetNumberOfClusters()<kNdiv||param1[irow].GetNumberOfClusters()<kNdiv) continue;
5821 kinks[irow].SetMother(param0[irow]);
5822 kinks[irow].SetDaughter(param1[irow]);
5823 kinks[irow].Update();
5826 // choose kink with best "quality"
5828 Double_t mindist = 10000;
5829 for (Int_t irow=0;irow<kNdiv;irow++){
5830 if (param0[irow].GetNumberOfClusters()<20||param1[irow].GetNumberOfClusters()<20) continue;
5831 if (TMath::Abs(kinks[irow].GetR())>240.) continue;
5832 if (TMath::Abs(kinks[irow].GetR())<100.) continue;
5834 Float_t normdist = TMath::Abs(param0[irow].GetX()-kinks[irow].GetR())*(0.1+kink.GetDistance());
5835 normdist/= (param0[irow].GetNumberOfClusters()+param1[irow].GetNumberOfClusters()+40.);
5836 if (normdist < mindist){
5842 if (index==-1) return 0;
5845 param0[index].Reset(kTRUE);
5846 FollowProlongation(param0[index],0);
5848 mother = param0[index];
5849 daughter = param1[index]; // daughter in vertex
5851 kink.SetMother(mother);
5852 kink.SetDaughter(daughter);
5854 kink.SetTPCRow0(GetRowNumber(kink.GetR()));
5855 kink.SetTPCncls(param0[index].GetNumberOfClusters(),0);
5856 kink.SetTPCncls(param1[index].GetNumberOfClusters(),1);
5857 kink.SetLabel(CookLabel(&mother,0.4, 0,kink.GetTPCRow0()),0);
5858 kink.SetLabel(CookLabel(&daughter,0.4, kink.GetTPCRow0(),160),1);
5859 mother.SetLabel(kink.GetLabel(0));
5860 daughter.SetLabel(kink.GetLabel(1));
5866 void AliTPCtrackerMI::UpdateKinkQualityM(AliTPCseed * seed){
5868 // update Kink quality information for mother after back propagation
5870 if (seed->GetKinkIndex(0)>=0) return;
5871 for (Int_t ikink=0;ikink<3;ikink++){
5872 Int_t index = seed->GetKinkIndex(ikink);
5873 if (index>=0) break;
5874 index = TMath::Abs(index)-1;
5875 AliESDkink * kink = fEvent->GetKink(index);
5876 kink->SetTPCDensity(-1,0,0);
5877 kink->SetTPCDensity(1,0,1);
5879 Int_t row0 = kink->GetTPCRow0() - 2 - Int_t( 0.5/ (0.05+kink->GetAngle(2)));
5880 if (row0<15) row0=15;
5882 Int_t row1 = kink->GetTPCRow0() + 2 + Int_t( 0.5/ (0.05+kink->GetAngle(2)));
5883 if (row1>145) row1=145;
5885 Int_t found,foundable,shared;
5886 seed->GetClusterStatistic(0,row0, found, foundable,shared,kFALSE);
5887 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),0,0);
5888 seed->GetClusterStatistic(row1,155, found, foundable,shared,kFALSE);
5889 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),0,1);
5894 void AliTPCtrackerMI::UpdateKinkQualityD(AliTPCseed * seed){
5896 // update Kink quality information for daughter after refit
5898 if (seed->GetKinkIndex(0)<=0) return;
5899 for (Int_t ikink=0;ikink<3;ikink++){
5900 Int_t index = seed->GetKinkIndex(ikink);
5901 if (index<=0) break;
5902 index = TMath::Abs(index)-1;
5903 AliESDkink * kink = fEvent->GetKink(index);
5904 kink->SetTPCDensity(-1,1,0);
5905 kink->SetTPCDensity(-1,1,1);
5907 Int_t row0 = kink->GetTPCRow0() -2 - Int_t( 0.5/ (0.05+kink->GetAngle(2)));
5908 if (row0<15) row0=15;
5910 Int_t row1 = kink->GetTPCRow0() +2 + Int_t( 0.5/ (0.05+kink->GetAngle(2)));
5911 if (row1>145) row1=145;
5913 Int_t found,foundable,shared;
5914 seed->GetClusterStatistic(0,row0, found, foundable,shared,kFALSE);
5915 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),1,0);
5916 seed->GetClusterStatistic(row1,155, found, foundable,shared,kFALSE);
5917 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),1,1);
5923 Int_t AliTPCtrackerMI::CheckKinkPoint(AliTPCseed*seed,AliTPCseed &mother, AliTPCseed &daughter, AliESDkink &knk)
5926 // check kink point for given track
5927 // if return value=0 kink point not found
5928 // otherwise seed0 correspond to mother particle
5929 // seed1 correspond to daughter particle
5930 // kink parameter of kink point
5931 AliKink &kink=(AliKink &)knk;
5933 Int_t middlerow = (seed->GetFirstPoint()+seed->GetLastPoint())/2;
5934 Int_t first = seed->GetFirstPoint();
5935 Int_t last = seed->GetLastPoint();
5936 if (last-first<20) return 0; // shortest length - 2*30 = 60 pad-rows
5939 AliTPCseed *seed1 = ReSeed(seed,middlerow+20, kTRUE); //middle of chamber
5940 if (!seed1) return 0;
5941 FollowProlongation(*seed1,seed->GetLastPoint()-20);
5942 seed1->Reset(kTRUE);
5943 FollowProlongation(*seed1,158);
5944 seed1->Reset(kTRUE);
5945 last = seed1->GetLastPoint();
5947 AliTPCseed *seed0 = new AliTPCseed(*seed);
5948 seed0->Reset(kFALSE);
5951 AliTPCseed param0[20]; // parameters along the track
5952 AliTPCseed param1[20]; // parameters along the track
5953 AliKink kinks[20]; // corresponding kink parameters
5955 for (Int_t irow=0; irow<20;irow++){
5956 rows[irow] = first +((last-first)*irow)/19;
5958 // store parameters along the track
5960 for (Int_t irow=0;irow<20;irow++){
5961 FollowBackProlongation(*seed0, rows[irow]);
5962 FollowProlongation(*seed1,rows[19-irow]);
5963 param0[irow] = *seed0;
5964 param1[19-irow] = *seed1;
5968 for (Int_t irow=0; irow<19;irow++){
5969 kinks[irow].SetMother(param0[irow]);
5970 kinks[irow].SetDaughter(param1[irow]);
5971 kinks[irow].Update();
5974 // choose kink with biggest change of angle
5976 Double_t maxchange= 0;
5977 for (Int_t irow=1;irow<19;irow++){
5978 if (TMath::Abs(kinks[irow].GetR())>240.) continue;
5979 if (TMath::Abs(kinks[irow].GetR())<110.) continue;
5980 Float_t quality = TMath::Abs(kinks[irow].GetAngle(2))/(3.+TMath::Abs(kinks[irow].GetR()-param0[irow].GetX()));
5981 if ( quality > maxchange){
5982 maxchange = quality;
5989 if (index<0) return 0;
5991 Int_t row0 = GetRowNumber(kinks[index].GetR()); //row 0 estimate
5992 seed0 = new AliTPCseed(param0[index]);
5993 seed1 = new AliTPCseed(param1[index]);
5994 seed0->Reset(kFALSE);
5995 seed1->Reset(kFALSE);
5996 seed0->ResetCovariance(10.);
5997 seed1->ResetCovariance(10.);
5998 FollowProlongation(*seed0,0);
5999 FollowBackProlongation(*seed1,158);
6000 mother = *seed0; // backup mother at position 0
6001 seed0->Reset(kFALSE);
6002 seed1->Reset(kFALSE);
6003 seed0->ResetCovariance(10.);
6004 seed1->ResetCovariance(10.);
6006 first = TMath::Max(row0-20,0);
6007 last = TMath::Min(row0+20,158);
6009 for (Int_t irow=0; irow<20;irow++){
6010 rows[irow] = first +((last-first)*irow)/19;
6012 // store parameters along the track
6014 for (Int_t irow=0;irow<20;irow++){
6015 FollowBackProlongation(*seed0, rows[irow]);
6016 FollowProlongation(*seed1,rows[19-irow]);
6017 param0[irow] = *seed0;
6018 param1[19-irow] = *seed1;
6022 for (Int_t irow=0; irow<19;irow++){
6023 kinks[irow].SetMother(param0[irow]);
6024 kinks[irow].SetDaughter(param1[irow]);
6025 // param0[irow].Dump();
6026 //param1[irow].Dump();
6027 kinks[irow].Update();
6030 // choose kink with biggest change of angle
6033 for (Int_t irow=0;irow<20;irow++){
6034 if (TMath::Abs(kinks[irow].GetR())>250.) continue;
6035 if (TMath::Abs(kinks[irow].GetR())<90.) continue;
6036 Float_t quality = TMath::Abs(kinks[irow].GetAngle(2))/(3.+TMath::Abs(kinks[irow].GetR()-param0[irow].GetX()));
6037 if ( quality > maxchange){
6038 maxchange = quality;
6045 if (index==-1 || param0[index].GetNumberOfClusters()+param1[index].GetNumberOfClusters()<100){
6050 // Float_t anglesigma = TMath::Sqrt(param0[index].fC22+param0[index].fC33+param1[index].fC22+param1[index].fC33);
6052 kink.SetMother(param0[index]);
6053 kink.SetDaughter(param1[index]);
6055 row0 = GetRowNumber(kink.GetR());
6056 kink.SetTPCRow0(row0);
6057 kink.SetLabel(CookLabel(seed0,0.5,0,row0),0);
6058 kink.SetLabel(CookLabel(seed1,0.5,row0,158),1);
6059 kink.SetIndex(-10,0);
6060 kink.SetIndex(int(param0[index].GetNumberOfClusters()+param1[index].GetNumberOfClusters()),1);
6061 kink.SetTPCncls(param0[index].GetNumberOfClusters(),0);
6062 kink.SetTPCncls(param1[index].GetNumberOfClusters(),1);
6065 // new (&mother) AliTPCseed(param0[index]);
6066 daughter = param1[index];
6067 daughter.SetLabel(kink.GetLabel(1));
6068 param0[index].Reset(kTRUE);
6069 FollowProlongation(param0[index],0);
6070 mother = param0[index];
6071 mother.SetLabel(kink.GetLabel(0));
6081 AliTPCseed* AliTPCtrackerMI::ReSeed(AliTPCseed *t)
6084 // reseed - refit - track
6087 // Int_t last = fSectors->GetNRows()-1;
6089 if (fSectors == fOuterSec){
6090 first = TMath::Max(first, t->GetFirstPoint()-fInnerSec->GetNRows());
6094 first = t->GetFirstPoint();
6096 AliTPCseed * seed = MakeSeed(t,0.1,0.5,0.9);
6097 FollowBackProlongation(*t,fSectors->GetNRows()-1);
6099 FollowProlongation(*t,first);
6109 //_____________________________________________________________________________
6110 Int_t AliTPCtrackerMI::ReadSeeds(const TFile *inp) {
6111 //-----------------------------------------------------------------
6112 // This function reades track seeds.
6113 //-----------------------------------------------------------------
6114 TDirectory *savedir=gDirectory;
6116 TFile *in=(TFile*)inp;
6117 if (!in->IsOpen()) {
6118 cerr<<"AliTPCtrackerMI::ReadSeeds(): input file is not open !\n";
6123 TTree *seedTree=(TTree*)in->Get("Seeds");
6125 cerr<<"AliTPCtrackerMI::ReadSeeds(): ";
6126 cerr<<"can't get a tree with track seeds !\n";
6129 AliTPCtrack *seed=new AliTPCtrack;
6130 seedTree->SetBranchAddress("tracks",&seed);
6132 if (fSeeds==0) fSeeds=new TObjArray(15000);
6134 Int_t n=(Int_t)seedTree->GetEntries();
6135 for (Int_t i=0; i<n; i++) {
6136 seedTree->GetEvent(i);
6137 fSeeds->AddLast(new AliTPCseed(*seed/*,seed->GetAlpha()*/));
6146 Int_t AliTPCtrackerMI::Clusters2Tracks (AliESDEvent *esd)
6149 if (fSeeds) DeleteSeeds();
6152 if (!fSeeds) return 1;
6159 //_____________________________________________________________________________
6160 Int_t AliTPCtrackerMI::Clusters2Tracks() {
6161 //-----------------------------------------------------------------
6162 // This is a track finder.
6163 //-----------------------------------------------------------------
6164 TDirectory *savedir=gDirectory;
6168 fSeeds = Tracking();
6171 Info("Clusters2Tracks","Time for tracking: \t");timer.Print();timer.Start();
6173 //activate again some tracks
6174 for (Int_t i=0; i<fSeeds->GetEntriesFast(); i++) {
6175 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
6177 Int_t nc=t.GetNumberOfClusters();
6179 delete fSeeds->RemoveAt(i);
6183 if (pt->GetRemoval()==10) {
6184 if (pt->GetDensityFirst(20)>0.8 || pt->GetDensityFirst(30)>0.8 || pt->GetDensityFirst(40)>0.7)
6185 pt->Desactivate(10); // make track again active
6187 pt->Desactivate(20);
6188 delete fSeeds->RemoveAt(i);
6193 RemoveUsed2(fSeeds,0.85,0.85,0);
6194 if (AliTPCReconstructor::GetRecoParam()->GetDoKinks()) FindKinks(fSeeds,fEvent);
6195 //FindCurling(fSeeds, fEvent,0);
6196 if (AliTPCReconstructor::StreamLevel()>2) FindMultiMC(fSeeds, fEvent,0); // find multi found tracks
6197 RemoveUsed2(fSeeds,0.5,0.4,20);
6198 FindSplitted(fSeeds, fEvent,0); // find multi found tracks
6201 // // refit short tracks
6203 Int_t nseed=fSeeds->GetEntriesFast();
6206 for (Int_t i=0; i<nseed; i++) {
6207 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
6209 Int_t nc=t.GetNumberOfClusters();
6211 delete fSeeds->RemoveAt(i);
6214 CookLabel(pt,0.1); //For comparison only
6215 //if ((pt->IsActive() || (pt->fRemoval==10) )&& nc>50 &&pt->GetNumberOfClusters()>0.4*pt->fNFoundable){
6216 if ((pt->IsActive() || (pt->GetRemoval()==10) )){
6218 if (fDebug>0) cerr<<found<<'\r';
6222 delete fSeeds->RemoveAt(i);
6226 //RemoveOverlap(fSeeds,0.99,7,kTRUE);
6228 //RemoveUsed(fSeeds,0.9,0.9,6);
6230 nseed=fSeeds->GetEntriesFast();
6232 for (Int_t i=0; i<nseed; i++) {
6233 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
6235 Int_t nc=t.GetNumberOfClusters();
6237 delete fSeeds->RemoveAt(i);
6241 t.CookdEdx(0.02,0.6);
6242 // CheckKinkPoint(&t,0.05);
6243 //if ((pt->IsActive() || (pt->fRemoval==10) )&& nc>50 &&pt->GetNumberOfClusters()>0.4*pt->fNFoundable){
6244 if ((pt->IsActive() || (pt->GetRemoval()==10) )){
6252 delete fSeeds->RemoveAt(i);
6253 //AliTPCseed * seed1 = ReSeed(pt,0.05,0.5,1);
6255 // FollowProlongation(*seed1,0);
6256 // Int_t n = seed1->GetNumberOfClusters();
6257 // printf("fP4\t%f\t%f\n",seed1->GetC(),pt->GetC());
6258 // printf("fN\t%d\t%d\n", seed1->GetNumberOfClusters(),pt->GetNumberOfClusters());
6261 //AliTPCseed * seed2 = ReSeed(pt,0.95,0.5,0.05);
6265 SortTracks(fSeeds, 1);
6269 PrepareForBackProlongation(fSeeds,5.);
6270 PropagateBack(fSeeds);
6271 printf("Time for back propagation: \t");timer.Print();timer.Start();
6275 PrepareForProlongation(fSeeds,5.);
6276 PropagateForward2(fSeeds);
6278 printf("Time for FORWARD propagation: \t");timer.Print();timer.Start();
6279 // RemoveUsed(fSeeds,0.7,0.7,6);
6280 //RemoveOverlap(fSeeds,0.9,7,kTRUE);
6282 nseed=fSeeds->GetEntriesFast();
6284 for (Int_t i=0; i<nseed; i++) {
6285 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
6287 Int_t nc=t.GetNumberOfClusters();
6289 delete fSeeds->RemoveAt(i);
6292 t.CookdEdx(0.02,0.6);
6293 // CookLabel(pt,0.1); //For comparison only
6294 //if ((pt->IsActive() || (pt->fRemoval==10) )&& nc>50 &&pt->GetNumberOfClusters()>0.4*pt->fNFoundable){
6295 if ((pt->IsActive() || (pt->fRemoval==10) )){
6296 cerr<<found++<<'\r';
6299 delete fSeeds->RemoveAt(i);
6304 // fNTracks = found;
6306 Info("Clusters2Tracks","Time for overlap removal, track writing and dedx cooking: \t"); timer.Print();timer.Start();
6309 // cerr<<"Number of found tracks : "<<"\t"<<found<<endl;
6310 Info("Clusters2Tracks","Number of found tracks %d",found);
6312 // UnloadClusters();
6317 void AliTPCtrackerMI::Tracking(TObjArray * arr)
6320 // tracking of the seeds
6323 fSectors = fOuterSec;
6324 ParallelTracking(arr,150,63);
6325 fSectors = fOuterSec;
6326 ParallelTracking(arr,63,0);
6329 TObjArray * AliTPCtrackerMI::Tracking(Int_t seedtype, Int_t i1, Int_t i2, Float_t cuts[4], Float_t dy, Int_t dsec)
6334 TObjArray * arr = new TObjArray;
6336 fSectors = fOuterSec;
6339 for (Int_t sec=0;sec<fkNOS;sec++){
6340 if (seedtype==3) MakeSeeds3(arr,sec,i1,i2,cuts,dy, dsec);
6341 if (seedtype==4) MakeSeeds5(arr,sec,i1,i2,cuts,dy);
6342 if (seedtype==2) MakeSeeds2(arr,sec,i1,i2,cuts,dy);
6345 Info("Tracking","\nSeeding - %d\t%d\t%d\t%d\n",seedtype,i1,i2,arr->GetEntriesFast());
6357 TObjArray * AliTPCtrackerMI::Tracking()
6361 if (AliTPCReconstructor::GetRecoParam()->GetSpecialSeeding()) return TrackingSpecial();
6364 Int_t nup=fOuterSec->GetNRows()+fInnerSec->GetNRows();
6366 TObjArray * seeds = new TObjArray;
6375 Float_t fnumber = 3.0;
6376 Float_t fdensity = 3.0;
6381 for (Int_t delta = 0; delta<18; delta+=6){
6385 arr = Tracking(3,nup-1-delta,nup-1-delta-gap,cuts,-1,1);
6386 SumTracks(seeds,arr);
6387 SignClusters(seeds,fnumber,fdensity);
6389 for (Int_t i=2;i<6;i+=2){
6390 // seed high pt tracks
6393 arr = Tracking(3,nup-i-delta,nup-i-delta-gap,cuts,-1,0);
6394 SumTracks(seeds,arr);
6395 SignClusters(seeds,fnumber,fdensity);
6400 // RemoveUsed(seeds,0.9,0.9,1);
6401 // UnsignClusters();
6402 // SignClusters(seeds,fnumber,fdensity);
6406 for (Int_t delta = 20; delta<120; delta+=10){
6408 // seed high pt tracks
6412 arr = Tracking(3,nup-delta,nup-delta-gap,cuts,-1);
6413 SumTracks(seeds,arr);
6414 SignClusters(seeds,fnumber,fdensity);
6419 arr = Tracking(3,nup-delta-5,nup-delta-5-gap,cuts,-1);
6420 SumTracks(seeds,arr);
6421 SignClusters(seeds,fnumber,fdensity);
6432 Info("Tracking()","\n\nPrimary seeding\t%d\n\n",seeds->GetEntriesFast());
6436 // RemoveUsed(seeds,0.75,0.75,1);
6438 //SignClusters(seeds,fnumber,fdensity);
6447 arr = Tracking(4,nup-1,nup-1-gap,cuts,-1);
6448 SumTracks(seeds,arr);
6449 SignClusters(seeds,fnumber,fdensity);
6451 arr = Tracking(4,nup-2,nup-2-gap,cuts,-1);
6452 SumTracks(seeds,arr);
6453 SignClusters(seeds,fnumber,fdensity);
6455 arr = Tracking(4,nup-3,nup-3-gap,cuts,-1);
6456 SumTracks(seeds,arr);
6457 SignClusters(seeds,fnumber,fdensity);
6461 for (Int_t delta = 3; delta<30; delta+=5){
6467 arr = Tracking(4,nup-1-delta,nup-1-delta-gap,cuts,-1);
6468 SumTracks(seeds,arr);
6469 SignClusters(seeds,fnumber,fdensity);
6471 arr = Tracking(4,nup-3-delta,nup-5-delta-gap,cuts,4);
6472 SumTracks(seeds,arr);
6473 SignClusters(seeds,fnumber,fdensity);
6484 Int_t fLastSeedRowSec=AliTPCReconstructor::GetRecoParam()->GetLastSeedRowSec();
6487 for (Int_t delta = 30; delta<fLastSeedRowSec; delta+=10){
6493 arr = Tracking(4,nup-1-delta,nup-1-delta-gap,cuts,-1);
6494 SumTracks(seeds,arr);
6495 SignClusters(seeds,fnumber,fdensity);
6497 arr = Tracking(4,nup-5-delta,nup-5-delta-gap,cuts,5 );
6498 SumTracks(seeds,arr);
6499 SignClusters(seeds,fnumber,fdensity);
6503 Info("Tracking()","\n\nSecondary seeding\t%d\n\n",seeds->GetEntriesFast());
6514 TObjArray * AliTPCtrackerMI::TrackingSpecial()
6517 // seeding adjusted for laser and cosmic tests - short tracks with big inclination angle
6518 // no primary vertex seeding tried
6522 Int_t nup=fOuterSec->GetNRows()+fInnerSec->GetNRows();
6524 TObjArray * seeds = new TObjArray;
6529 Float_t fnumber = 3.0;
6530 Float_t fdensity = 3.0;
6533 cuts[0] = AliTPCReconstructor::GetRecoParam()->GetMaxC(); // max curvature
6534 cuts[1] = 3.5; // max tan(phi) angle for seeding
6535 cuts[2] = 3.; // not used (cut on z primary vertex)
6536 cuts[3] = 3.5; // max tan(theta) angle for seeding
6538 for (Int_t delta = 0; nup-delta-gap-1>0; delta+=3){
6540 arr = Tracking(4,nup-1-delta,nup-1-delta-gap,cuts,-1);
6541 SumTracks(seeds,arr);
6542 SignClusters(seeds,fnumber,fdensity);
6546 Info("Tracking()","\n\nSecondary seeding\t%d\n\n",seeds->GetEntriesFast());
6557 void AliTPCtrackerMI::SumTracks(TObjArray *arr1,TObjArray *&arr2) const
6560 //sum tracks to common container
6561 //remove suspicious tracks
6562 Int_t nseed = arr2->GetEntriesFast();
6563 for (Int_t i=0;i<nseed;i++){
6564 AliTPCseed *pt=(AliTPCseed*)arr2->UncheckedAt(i);
6567 // remove tracks with too big curvature
6569 if (TMath::Abs(pt->GetC())>AliTPCReconstructor::GetRecoParam()->GetMaxC()){
6570 delete arr2->RemoveAt(i);
6573 // REMOVE VERY SHORT TRACKS
6574 if (pt->GetNumberOfClusters()<20){
6575 delete arr2->RemoveAt(i);
6578 // NORMAL ACTIVE TRACK
6579 if (pt->IsActive()){
6580 arr1->AddLast(arr2->RemoveAt(i));
6583 //remove not usable tracks
6584 if (pt->GetRemoval()!=10){
6585 delete arr2->RemoveAt(i);
6589 // ENABLE ONLY ENOUGH GOOD STOPPED TRACKS
6590 if (pt->GetDensityFirst(20)>0.8 || pt->GetDensityFirst(30)>0.8 || pt->GetDensityFirst(40)>0.7)
6591 arr1->AddLast(arr2->RemoveAt(i));
6593 delete arr2->RemoveAt(i);
6597 delete arr2; arr2 = 0;
6602 void AliTPCtrackerMI::ParallelTracking(TObjArray * arr, Int_t rfirst, Int_t rlast)
6605 // try to track in parralel
6607 Int_t nseed=arr->GetEntriesFast();
6608 //prepare seeds for tracking
6609 for (Int_t i=0; i<nseed; i++) {
6610 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i), &t=*pt;
6612 if (!t.IsActive()) continue;
6613 // follow prolongation to the first layer
6614 if ( (fSectors ==fInnerSec) || (t.GetFirstPoint()-fParam->GetNRowLow()>rfirst+1) )
6615 FollowProlongation(t, rfirst+1);
6620 for (Int_t nr=rfirst; nr>=rlast; nr--){
6621 if (nr<fInnerSec->GetNRows())
6622 fSectors = fInnerSec;
6624 fSectors = fOuterSec;
6625 // make indexes with the cluster tracks for given
6627 // find nearest cluster
6628 for (Int_t i=0; i<nseed; i++) {
6629 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i), &t=*pt;
6631 if (nr==80) pt->UpdateReference();
6632 if (!pt->IsActive()) continue;
6633 // if ( (fSectors ==fOuterSec) && (pt->fFirstPoint-fParam->GetNRowLow())<nr) continue;
6634 if (pt->GetRelativeSector()>17) {
6637 UpdateClusters(t,nr);
6639 // prolonagate to the nearest cluster - if founded
6640 for (Int_t i=0; i<nseed; i++) {
6641 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
6643 if (!pt->IsActive()) continue;
6644 // if ((fSectors ==fOuterSec) && (pt->fFirstPoint-fParam->GetNRowLow())<nr) continue;
6645 if (pt->GetRelativeSector()>17) {
6648 FollowToNextCluster(*pt,nr);
6653 void AliTPCtrackerMI::PrepareForBackProlongation(TObjArray * arr,Float_t fac) const
6657 // if we use TPC track itself we have to "update" covariance
6659 Int_t nseed= arr->GetEntriesFast();
6660 for (Int_t i=0;i<nseed;i++){
6661 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
6665 //rotate to current local system at first accepted point
6666 Int_t index = pt->GetClusterIndex2(pt->GetFirstPoint());
6667 Int_t sec = (index&0xff000000)>>24;
6669 Float_t angle1 = fInnerSec->GetAlpha()*sec+fInnerSec->GetAlphaShift();
6670 if (angle1>TMath::Pi())
6671 angle1-=2.*TMath::Pi();
6672 Float_t angle2 = pt->GetAlpha();
6674 if (TMath::Abs(angle1-angle2)>0.001){
6675 pt->Rotate(angle1-angle2);
6676 //angle2 = pt->GetAlpha();
6677 //pt->fRelativeSector = pt->GetAlpha()/fInnerSec->GetAlpha();
6678 //if (pt->GetAlpha()<0)
6679 // pt->fRelativeSector+=18;
6680 //sec = pt->fRelativeSector;
6689 void AliTPCtrackerMI::PrepareForProlongation(TObjArray * arr, Float_t fac) const
6693 // if we use TPC track itself we have to "update" covariance
6695 Int_t nseed= arr->GetEntriesFast();
6696 for (Int_t i=0;i<nseed;i++){
6697 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
6700 pt->SetFirstPoint(pt->GetLastPoint());
6708 Int_t AliTPCtrackerMI::PropagateBack(TObjArray * arr)
6711 // make back propagation
6713 Int_t nseed= arr->GetEntriesFast();
6714 for (Int_t i=0;i<nseed;i++){
6715 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
6716 if (pt&& pt->GetKinkIndex(0)<=0) {
6717 //AliTPCseed *pt2 = new AliTPCseed(*pt);
6718 fSectors = fInnerSec;
6719 //FollowBackProlongation(*pt,fInnerSec->GetNRows()-1);
6720 //fSectors = fOuterSec;
6721 FollowBackProlongation(*pt,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1);
6722 //if (pt->GetNumberOfClusters()<(pt->fEsd->GetTPCclusters(0)) ){
6723 // Error("PropagateBack","Not prolonged track %d",pt->GetLabel());
6724 // FollowBackProlongation(*pt2,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1);
6727 if (pt&& pt->GetKinkIndex(0)>0) {
6728 AliESDkink * kink = fEvent->GetKink(pt->GetKinkIndex(0)-1);
6729 pt->SetFirstPoint(kink->GetTPCRow0());
6730 fSectors = fInnerSec;
6731 FollowBackProlongation(*pt,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1);
6739 Int_t AliTPCtrackerMI::PropagateForward2(TObjArray * arr)
6742 // make forward propagation
6744 Int_t nseed= arr->GetEntriesFast();
6746 for (Int_t i=0;i<nseed;i++){
6747 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
6749 FollowProlongation(*pt,0);
6758 Int_t AliTPCtrackerMI::PropagateForward()
6761 // propagate track forward
6763 Int_t nseed = fSeeds->GetEntriesFast();
6764 for (Int_t i=0;i<nseed;i++){
6765 AliTPCseed *pt = (AliTPCseed*)fSeeds->UncheckedAt(i);
6767 AliTPCseed &t = *pt;
6768 Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
6769 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
6770 if (alpha < 0. ) alpha += 2.*TMath::Pi();
6771 t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN);
6775 fSectors = fOuterSec;
6776 ParallelTracking(fSeeds,fOuterSec->GetNRows()+fInnerSec->GetNRows()-1,fInnerSec->GetNRows());
6777 fSectors = fInnerSec;
6778 ParallelTracking(fSeeds,fInnerSec->GetNRows()-1,0);
6787 Int_t AliTPCtrackerMI::PropagateBack(AliTPCseed * pt, Int_t row0, Int_t row1)
6790 // make back propagation, in between row0 and row1
6794 fSectors = fInnerSec;
6797 if (row1<fSectors->GetNRows())
6800 r1 = fSectors->GetNRows()-1;
6802 if (row0<fSectors->GetNRows()&& r1>0 )
6803 FollowBackProlongation(*pt,r1);
6804 if (row1<=fSectors->GetNRows())
6807 r1 = row1 - fSectors->GetNRows();
6808 if (r1<=0) return 0;
6809 if (r1>=fOuterSec->GetNRows()) return 0;
6810 fSectors = fOuterSec;
6811 return FollowBackProlongation(*pt,r1);
6819 void AliTPCtrackerMI::GetShape(AliTPCseed * seed, Int_t row)
6823 AliTPCClusterParam * clparam = AliTPCcalibDB::Instance()->GetClusterParam();
6824 Float_t zdrift = TMath::Abs((fParam->GetZLength(0)-TMath::Abs(seed->GetZ())));
6825 Int_t type = (seed->GetSector() < fParam->GetNSector()/2) ? 0: (row>126) ? 1:2;
6826 Double_t angulary = seed->GetSnp();
6827 angulary = angulary*angulary/(1.-angulary*angulary);
6828 Double_t angularz = seed->GetTgl()*seed->GetTgl()*(1.+angulary);
6830 Double_t sigmay = clparam->GetRMS0(0,type,zdrift,TMath::Sqrt(TMath::Abs(angulary)));
6831 Double_t sigmaz = clparam->GetRMS0(1,type,zdrift,TMath::Sqrt(TMath::Abs(angularz)));
6832 seed->SetCurrentSigmaY2(sigmay*sigmay);
6833 seed->SetCurrentSigmaZ2(sigmaz*sigmaz);
6834 // Float_t sd2 = TMath::Abs((fParam->GetZLength(0)-TMath::Abs(seed->GetZ())))*fParam->GetDiffL()*fParam->GetDiffL();
6835 // // Float_t padlength = fParam->GetPadPitchLength(seed->fSector);
6836 // Float_t padlength = GetPadPitchLength(row);
6838 // Float_t sresy = (seed->GetSector() < fParam->GetNSector()/2) ? 0.2 :0.3;
6839 // seed->SetCurrentSigmaY2(sd2+padlength*padlength*angulary/12.+sresy*sresy);
6841 // Float_t sresz = fParam->GetZSigma();
6842 // seed->SetCurrentSigmaZ2(sd2+padlength*padlength*angularz*angularz*(1+angulary)/12.+sresz*sresz);
6844 Float_t wy = GetSigmaY(seed);
6845 Float_t wz = GetSigmaZ(seed);
6848 if (TMath::Abs(wy/seed->fCurrentSigmaY2-1)>0.0001 || TMath::Abs(wz/seed->fCurrentSigmaZ2-1)>0.0001 ){
6849 printf("problem\n");
6856 //__________________________________________________________________________
6857 void AliTPCtrackerMI::CookLabel(AliKalmanTrack *tk, Float_t wrong) const {
6858 //--------------------------------------------------------------------
6859 //This function "cooks" a track label. If label<0, this track is fake.
6860 //--------------------------------------------------------------------
6861 AliTPCseed * t = dynamic_cast<AliTPCseed*>(tk);
6863 printf("%s:%d wrong type \n",(char*)__FILE__,__LINE__);
6867 Int_t noc=t->GetNumberOfClusters();
6869 //printf("\nnot founded prolongation\n\n\n");
6875 AliTPCclusterMI *clusters[160];
6877 for (Int_t i=0;i<160;i++) {
6884 for (i=0; i<160 && current<noc; i++) {
6886 Int_t index=t->GetClusterIndex2(i);
6887 if (index<=0) continue;
6888 if (index&0x8000) continue;
6890 //clusters[current]=GetClusterMI(index);
6891 if (t->GetClusterPointer(i)){
6892 clusters[current]=t->GetClusterPointer(i);
6898 Int_t lab=123456789;
6899 for (i=0; i<noc; i++) {
6900 AliTPCclusterMI *c=clusters[i];
6902 lab=TMath::Abs(c->GetLabel(0));
6904 for (j=0; j<noc; j++) if (lb[j]==lab || mx[j]==0) break;
6910 for (i=0; i<noc; i++) if (mx[i]>max) {max=mx[i]; lab=lb[i];}
6912 for (i=0; i<noc; i++) {
6913 AliTPCclusterMI *c=clusters[i];
6915 if (TMath::Abs(c->GetLabel(1)) == lab ||
6916 TMath::Abs(c->GetLabel(2)) == lab ) max++;
6919 if ((1.- Float_t(max)/noc) > wrong) lab=-lab;
6922 Int_t tail=Int_t(0.10*noc);
6925 for (i=1; i<=160&&ind<tail; i++) {
6926 // AliTPCclusterMI *c=clusters[noc-i];
6927 AliTPCclusterMI *c=clusters[i];
6929 if (lab == TMath::Abs(c->GetLabel(0)) ||
6930 lab == TMath::Abs(c->GetLabel(1)) ||
6931 lab == TMath::Abs(c->GetLabel(2))) max++;
6934 if (max < Int_t(0.5*tail)) lab=-lab;
6941 //delete[] clusters;
6945 //__________________________________________________________________________
6946 Int_t AliTPCtrackerMI::CookLabel(AliTPCseed *t, Float_t wrong,Int_t first, Int_t last) const {
6947 //--------------------------------------------------------------------
6948 //This function "cooks" a track label. If label<0, this track is fake.
6949 //--------------------------------------------------------------------
6950 Int_t noc=t->GetNumberOfClusters();
6952 //printf("\nnot founded prolongation\n\n\n");
6958 AliTPCclusterMI *clusters[160];
6960 for (Int_t i=0;i<160;i++) {
6967 for (i=0; i<160 && current<noc; i++) {
6968 if (i<first) continue;
6969 if (i>last) continue;
6970 Int_t index=t->GetClusterIndex2(i);
6971 if (index<=0) continue;
6972 if (index&0x8000) continue;
6974 //clusters[current]=GetClusterMI(index);
6975 if (t->GetClusterPointer(i)){
6976 clusters[current]=t->GetClusterPointer(i);
6981 if (noc<5) return -1;
6982 Int_t lab=123456789;
6983 for (i=0; i<noc; i++) {
6984 AliTPCclusterMI *c=clusters[i];
6986 lab=TMath::Abs(c->GetLabel(0));
6988 for (j=0; j<noc; j++) if (lb[j]==lab || mx[j]==0) break;
6994 for (i=0; i<noc; i++) if (mx[i]>max) {max=mx[i]; lab=lb[i];}
6996 for (i=0; i<noc; i++) {
6997 AliTPCclusterMI *c=clusters[i];
6999 if (TMath::Abs(c->GetLabel(1)) == lab ||
7000 TMath::Abs(c->GetLabel(2)) == lab ) max++;
7003 if ((1.- Float_t(max)/noc) > wrong) lab=-lab;
7006 Int_t tail=Int_t(0.10*noc);
7009 for (i=1; i<=160&&ind<tail; i++) {
7010 // AliTPCclusterMI *c=clusters[noc-i];
7011 AliTPCclusterMI *c=clusters[i];
7013 if (lab == TMath::Abs(c->GetLabel(0)) ||
7014 lab == TMath::Abs(c->GetLabel(1)) ||
7015 lab == TMath::Abs(c->GetLabel(2))) max++;
7018 if (max < Int_t(0.5*tail)) lab=-lab;
7021 // t->SetLabel(lab);
7025 //delete[] clusters;
7029 Int_t AliTPCtrackerMI::GetRowNumber(Double_t x[3]) const
7031 //return pad row number for given x vector
7032 Float_t phi = TMath::ATan2(x[1],x[0]);
7033 if(phi<0) phi=2.*TMath::Pi()+phi;
7034 // Get the local angle in the sector philoc
7035 const Float_t kRaddeg = 180/3.14159265358979312;
7036 Float_t phiangle = (Int_t (phi*kRaddeg/20.) + 0.5)*20./kRaddeg;
7037 Double_t localx = x[0]*TMath::Cos(phiangle)-x[1]*TMath::Sin(phiangle);
7038 return GetRowNumber(localx);
7043 void AliTPCtrackerMI::MakeBitmaps(AliTPCseed *t)
7045 //-----------------------------------------------------------------------
7046 // Fill the cluster and sharing bitmaps of the track
7047 //-----------------------------------------------------------------------
7049 Int_t firstpoint = 0;
7050 Int_t lastpoint = 159;
7051 AliTPCTrackerPoint *point;
7053 for (int iter=firstpoint; iter<lastpoint; iter++) {
7054 point = t->GetTrackPoint(iter);
7056 t->SetClusterMapBit(iter, kTRUE);
7057 if (point->IsShared())
7058 t->SetSharedMapBit(iter,kTRUE);
7060 t->SetSharedMapBit(iter, kFALSE);
7063 t->SetClusterMapBit(iter, kFALSE);
7064 t->SetSharedMapBit(iter, kFALSE);
7071 void AliTPCtrackerMI::AddCovariance(AliTPCseed * seed){
7073 // Adding systematic error
7074 // !!!! the systematic error for element 4 is in 1/cm not in pt
7076 const Double_t *param = AliTPCReconstructor::GetRecoParam()->GetSystematicError();
7078 for (Int_t i=0;i<15;i++) covar[i]=0;
7084 covar[0] = param[0]*param[0];
7085 covar[2] = param[1]*param[1];
7086 covar[5] = param[2]*param[2];
7087 covar[9] = param[3]*param[3];
7088 Double_t facC = AliTracker::GetBz()*kB2C;
7089 covar[14]= param[4]*param[4]*facC*facC;
7090 seed->AddCovariance(covar);