1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
17 //-------------------------------------------------------
18 // Implementation of the TPC tracker
20 // Origin: Marian Ivanov Marian.Ivanov@cern.ch
22 // AliTPC parallel tracker
24 // The track fitting is based on Kalaman filtering approach
26 // The track finding steps:
27 // 1. Seeding - with and without vertex constraint
28 // - seeding with vertex constain done at first n^2 proble
29 // - seeding without vertex constraint n^3 problem
30 // 2. Tracking - follow prolongation road - find cluster - update kalman track
32 // The seeding and tracking is repeated several times, in different seeding region.
33 // This approach enables to find the track which cannot be seeded in some region of TPC
34 // This can happen because of low momenta (track do not reach outer radius), or track is currently in the ded region between sectors, or the track is for the moment overlapped with other track (seed quality is poor) ...
36 // With this approach we reach almost 100 % efficiency also for high occupancy events.
37 // (If the seeding efficiency in a region is about 90 % than with logical or of several
38 // regions we will reach 100% (in theory - supposing independence)
40 // Repeating several seeding - tracking procedures some of the tracks can be find
43 // The procedures to remove multi find tacks are impremented:
44 // RemoveUsed2 - fast procedure n problem -
45 // Algorithm - Sorting tracks according quality
46 // remove tracks with some shared fraction
47 // Sharing in respect to all tacks
48 // Signing clusters in gold region
49 // FindSplitted - slower algorithm n^2
50 // Sort the tracks according quality
51 // Loop over pair of tracks
52 // If overlap with other track bigger than threshold - remove track
54 // FindCurling - Finds the pair of tracks which are curling
55 // - About 10% of tracks can be find with this procedure
56 // The combinatorial background is too big to be used in High
57 // multiplicity environment
58 // - n^2 problem - Slow procedure - currently it is disabled because of
61 // The number of splitted tracks can be reduced disabling the sharing of the cluster.
62 // tpcRecoParam-> SetClusterSharing(kFALSE);
63 // IT IS HIGHLY non recomended to use it in high flux enviroonment
64 // Even using this switch some tracks can be found more than once
65 // (because of multiple seeding and low quality tracks which will not cross full chamber)
68 // The tracker itself can be debugged - the information about tracks can be stored in several // phases of the reconstruction
69 // To enable storage of the TPC tracks in the ESD friend track
70 // use AliTPCReconstructor::SetStreamLevel(n); where nis bigger 0
72 // The debug level - different procedure produce tree for numerical debugging
73 // To enable them set AliTPCReconstructor::SetStreamLevel(n); where nis bigger 1
77 // Adding systematic errors to the covariance:
79 // The systematic errors due to the misalignment and miscalibration are added to the covariance matrix
80 // of the tracks (not to the clusters as they are dependent):
81 // The parameters form AliTPCRecoParam are used AliTPCRecoParam::GetSystematicError
82 // The systematic errors are expressed there in RMS - position (cm), angle (rad), curvature (1/cm)
83 // The default values are 0.
85 // The sytematic errors are added to the covariance matrix in following places:
87 // 1. During fisrt itteration - AliTPCtrackerMI::FillESD
88 // 2. Second iteration -
89 // 2.a ITS->TPC - AliTPCtrackerMI::ReadSeeds
90 // 2.b TPC->TRD - AliTPCtrackerMI::PropagateBack
91 // 3. Third iteration -
92 // 3.a TRD->TPC - AliTPCtrackerMI::ReadSeeds
93 // 3.b TPC->ITS - AliTPCtrackerMI::RefitInward
95 // There are several places in the code which can be numerically debuged
96 // This code is keeped in order to enable code development and to check the calibration implementtion
98 // 1. ErrParam stream (Log level 9) - dump information about
100 // 2.a) cluster error estimate
101 // 3.a) cluster shape estimate
104 //-------------------------------------------------------
109 #include "Riostream.h"
110 #include <TClonesArray.h>
112 #include <TObjArray.h>
115 #include "AliComplexCluster.h"
116 #include "AliESDEvent.h"
117 #include "AliESDtrack.h"
118 #include "AliESDVertex.h"
121 #include "AliHelix.h"
122 #include "AliRunLoader.h"
123 #include "AliTPCClustersRow.h"
124 #include "AliTPCParam.h"
125 #include "AliTPCReconstructor.h"
126 #include "AliTPCpolyTrack.h"
127 #include "AliTPCreco.h"
128 #include "AliTPCseed.h"
130 #include "AliTPCtrackerSector.h"
131 #include "AliTPCtrackerMI.h"
132 #include "TStopwatch.h"
133 #include "AliTPCReconstructor.h"
135 #include "TTreeStream.h"
136 #include "AliAlignObj.h"
137 #include "AliTrackPointArray.h"
139 #include "AliTPCcalibDB.h"
140 #include "AliTPCTransform.h"
141 #include "AliTPCClusterParam.h"
145 ClassImp(AliTPCtrackerMI)
149 class AliTPCFastMath {
152 static Double_t FastAsin(Double_t x);
154 static Double_t fgFastAsin[20000]; //lookup table for fast asin computation
157 Double_t AliTPCFastMath::fgFastAsin[20000];
158 AliTPCFastMath gAliTPCFastMath; // needed to fill the LUT
160 AliTPCFastMath::AliTPCFastMath(){
162 // initialized lookup table;
163 for (Int_t i=0;i<10000;i++){
164 fgFastAsin[2*i] = TMath::ASin(i/10000.);
165 fgFastAsin[2*i+1] = (TMath::ASin((i+1)/10000.)-fgFastAsin[2*i]);
169 Double_t AliTPCFastMath::FastAsin(Double_t x){
171 // return asin using lookup table
173 Int_t index = int(x*10000);
174 return fgFastAsin[2*index]+(x*10000.-index)*fgFastAsin[2*index+1];
177 Int_t index = int(x*10000);
178 return -(fgFastAsin[2*index]+(x*10000.-index)*fgFastAsin[2*index+1]);
180 //__________________________________________________________________
181 AliTPCtrackerMI::AliTPCtrackerMI()
203 // default constructor
206 //_____________________________________________________________________
210 Int_t AliTPCtrackerMI::UpdateTrack(AliTPCseed * track, Int_t accept){
212 //update track information using current cluster - track->fCurrentCluster
215 AliTPCclusterMI* c =track->GetCurrentCluster();
216 if (accept>0) track->SetCurrentClusterIndex1(track->GetCurrentClusterIndex1() | 0x8000); //sign not accepted clusters
218 UInt_t i = track->GetCurrentClusterIndex1();
220 Int_t sec=(i&0xff000000)>>24;
221 //Int_t row = (i&0x00ff0000)>>16;
222 track->SetRow((i&0x00ff0000)>>16);
223 track->SetSector(sec);
224 // Int_t index = i&0xFFFF;
225 if (sec>=fParam->GetNInnerSector()) track->SetRow(track->GetRow()+fParam->GetNRowLow());
226 track->SetClusterIndex2(track->GetRow(), i);
227 //track->fFirstPoint = row;
228 //if ( track->fLastPoint<row) track->fLastPoint =row;
229 // if (track->fRow<0 || track->fRow>160) {
230 // printf("problem\n");
232 if (track->GetFirstPoint()>track->GetRow())
233 track->SetFirstPoint(track->GetRow());
234 if (track->GetLastPoint()<track->GetRow())
235 track->SetLastPoint(track->GetRow());
238 track->SetClusterPointer(track->GetRow(),c);
241 Double_t angle2 = track->GetSnp()*track->GetSnp();
243 //SET NEW Track Point
245 if (angle2<1) //PH sometimes angle2 is very big. To be investigated...
247 angle2 = TMath::Sqrt(angle2/(1-angle2));
248 AliTPCTrackerPoint &point =*(track->GetTrackPoint(track->GetRow()));
250 point.SetSigmaY(c->GetSigmaY2()/track->GetCurrentSigmaY2());
251 point.SetSigmaZ(c->GetSigmaZ2()/track->GetCurrentSigmaZ2());
252 point.SetErrY(sqrt(track->GetErrorY2()));
253 point.SetErrZ(sqrt(track->GetErrorZ2()));
255 point.SetX(track->GetX());
256 point.SetY(track->GetY());
257 point.SetZ(track->GetZ());
258 point.SetAngleY(angle2);
259 point.SetAngleZ(track->GetTgl());
260 if (point.IsShared()){
261 track->SetErrorY2(track->GetErrorY2()*4);
262 track->SetErrorZ2(track->GetErrorZ2()*4);
266 Double_t chi2 = track->GetPredictedChi2(track->GetCurrentCluster());
268 // track->SetErrorY2(track->GetErrorY2()*1.3);
269 // track->SetErrorY2(track->GetErrorY2()+0.01);
270 // track->SetErrorZ2(track->GetErrorZ2()*1.3);
271 // track->SetErrorZ2(track->GetErrorZ2()+0.005);
273 if (accept>0) return 0;
274 if (track->GetNumberOfClusters()%20==0){
275 // if (track->fHelixIn){
276 // TClonesArray & larr = *(track->fHelixIn);
277 // Int_t ihelix = larr.GetEntriesFast();
278 // new(larr[ihelix]) AliHelix(*track) ;
281 track->SetNoCluster(0);
282 return track->Update(c,chi2,i);
287 Int_t AliTPCtrackerMI::AcceptCluster(AliTPCseed * seed, AliTPCclusterMI * cluster)
290 // decide according desired precision to accept given
291 // cluster for tracking
292 Double_t sy2=ErrY2(seed,cluster);
293 Double_t sz2=ErrZ2(seed,cluster);
295 Double_t sdistancey2 = sy2+seed->GetSigmaY2();
296 Double_t sdistancez2 = sz2+seed->GetSigmaZ2();
298 Double_t rdistancey2 = (seed->GetCurrentCluster()->GetY()-seed->GetY())*
299 (seed->GetCurrentCluster()->GetY()-seed->GetY())/sdistancey2;
300 Double_t rdistancez2 = (seed->GetCurrentCluster()->GetZ()-seed->GetZ())*
301 (seed->GetCurrentCluster()->GetZ()-seed->GetZ())/sdistancez2;
303 Double_t rdistance2 = rdistancey2+rdistancez2;
306 if (AliTPCReconstructor::StreamLevel()>5 && seed->GetNumberOfClusters()>20) {
307 Float_t rmsy2 = seed->GetCurrentSigmaY2();
308 Float_t rmsz2 = seed->GetCurrentSigmaZ2();
309 Float_t rmsy2p30 = seed->GetCMeanSigmaY2p30();
310 Float_t rmsz2p30 = seed->GetCMeanSigmaZ2p30();
311 Float_t rmsy2p30R = seed->GetCMeanSigmaY2p30R();
312 Float_t rmsz2p30R = seed->GetCMeanSigmaZ2p30R();
314 AliExternalTrackParam param(*seed);
315 static TVectorD gcl(3),gtr(3);
317 param.GetXYZ(gcl.GetMatrixArray());
318 cluster->GetGlobalXYZ(gclf);
319 gcl[0]=gclf[0]; gcl[1]=gclf[1]; gcl[2]=gclf[2];
320 (*fDebugStreamer)<<"ErrParam"<<
329 "rmsy2p30="<<rmsy2p30<<
330 "rmsz2p30="<<rmsz2p30<<
331 "rmsy2p30R="<<rmsy2p30R<<
332 "rmsz2p30R="<<rmsz2p30R<<
336 if (rdistance2>16) return 3;
339 if ((rdistancey2>9. || rdistancez2>9.) && cluster->GetType()==0)
340 return 2; //suspisiouce - will be changed
342 if ((rdistancey2>6.25 || rdistancez2>6.25) && cluster->GetType()>0)
343 // strict cut on overlaped cluster
344 return 2; //suspisiouce - will be changed
346 if ( (rdistancey2>1. || rdistancez2>6.25 )
347 && cluster->GetType()<0){
348 seed->SetNFoundable(seed->GetNFoundable()-1);
357 //_____________________________________________________________________________
358 AliTPCtrackerMI::AliTPCtrackerMI(const AliTPCParam *par):
360 fkNIS(par->GetNInnerSector()/2),
362 fkNOS(par->GetNOuterSector()/2),
379 //---------------------------------------------------------------------
380 // The main TPC tracker constructor
381 //---------------------------------------------------------------------
382 fInnerSec=new AliTPCtrackerSector[fkNIS];
383 fOuterSec=new AliTPCtrackerSector[fkNOS];
386 for (i=0; i<fkNIS; i++) fInnerSec[i].Setup(par,0);
387 for (i=0; i<fkNOS; i++) fOuterSec[i].Setup(par,1);
390 Int_t nrowlow = par->GetNRowLow();
391 Int_t nrowup = par->GetNRowUp();
394 for (Int_t i=0;i<nrowlow;i++){
395 fXRow[i] = par->GetPadRowRadiiLow(i);
396 fPadLength[i]= par->GetPadPitchLength(0,i);
397 fYMax[i] = fXRow[i]*TMath::Tan(0.5*par->GetInnerAngle());
401 for (Int_t i=0;i<nrowup;i++){
402 fXRow[i+nrowlow] = par->GetPadRowRadiiUp(i);
403 fPadLength[i+nrowlow] = par->GetPadPitchLength(60,i);
404 fYMax[i+nrowlow] = fXRow[i+nrowlow]*TMath::Tan(0.5*par->GetOuterAngle());
407 fDebugStreamer = new TTreeSRedirector("TPCdebug.root");
409 //________________________________________________________________________
410 AliTPCtrackerMI::AliTPCtrackerMI(const AliTPCtrackerMI &t):
431 //------------------------------------
432 // dummy copy constructor
433 //------------------------------------------------------------------
436 AliTPCtrackerMI & AliTPCtrackerMI::operator=(const AliTPCtrackerMI& /*r*/){
437 //------------------------------
439 //--------------------------------------------------------------
442 //_____________________________________________________________________________
443 AliTPCtrackerMI::~AliTPCtrackerMI() {
444 //------------------------------------------------------------------
445 // TPC tracker destructor
446 //------------------------------------------------------------------
453 if (fDebugStreamer) delete fDebugStreamer;
457 void AliTPCtrackerMI::FillESD(TObjArray* arr)
461 //fill esds using updated tracks
463 // write tracks to the event
464 // store index of the track
465 Int_t nseed=arr->GetEntriesFast();
466 //FindKinks(arr,fEvent);
467 for (Int_t i=0; i<nseed; i++) {
468 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
472 if (AliTPCReconstructor::StreamLevel()>1) {
473 (*fDebugStreamer)<<"Track0"<<
477 // pt->PropagateTo(fParam->GetInnerRadiusLow());
478 if (pt->GetKinkIndex(0)<=0){ //don't propagate daughter tracks
479 pt->PropagateTo(fParam->GetInnerRadiusLow());
482 if (( pt->GetPoints()[2]- pt->GetPoints()[0])>5 && pt->GetPoints()[3]>0.8){
484 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
485 iotrack.SetTPCPoints(pt->GetPoints());
486 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
487 iotrack.SetV0Indexes(pt->GetV0Indexes());
488 // iotrack.SetTPCpid(pt->fTPCr);
489 //iotrack.SetTPCindex(i);
490 fEvent->AddTrack(&iotrack);
494 if ( (pt->GetNumberOfClusters()>70)&& (Float_t(pt->GetNumberOfClusters())/Float_t(pt->GetNFoundable()))>0.55) {
496 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
497 iotrack.SetTPCPoints(pt->GetPoints());
498 //iotrack.SetTPCindex(i);
499 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
500 iotrack.SetV0Indexes(pt->GetV0Indexes());
501 // iotrack.SetTPCpid(pt->fTPCr);
502 fEvent->AddTrack(&iotrack);
506 // short tracks - maybe decays
508 if ( (pt->GetNumberOfClusters()>30) && (Float_t(pt->GetNumberOfClusters())/Float_t(pt->GetNFoundable()))>0.70) {
509 Int_t found,foundable,shared;
510 pt->GetClusterStatistic(0,60,found, foundable,shared,kFALSE);
511 if ( (found>20) && (pt->GetNShared()/float(pt->GetNumberOfClusters())<0.2)){
513 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
514 //iotrack.SetTPCindex(i);
515 iotrack.SetTPCPoints(pt->GetPoints());
516 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
517 iotrack.SetV0Indexes(pt->GetV0Indexes());
518 //iotrack.SetTPCpid(pt->fTPCr);
519 fEvent->AddTrack(&iotrack);
524 if ( (pt->GetNumberOfClusters()>20) && (Float_t(pt->GetNumberOfClusters())/Float_t(pt->GetNFoundable()))>0.8) {
525 Int_t found,foundable,shared;
526 pt->GetClusterStatistic(0,60,found, foundable,shared,kFALSE);
527 if (found<20) continue;
528 if (pt->GetNShared()/float(pt->GetNumberOfClusters())>0.2) continue;
531 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
532 iotrack.SetTPCPoints(pt->GetPoints());
533 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
534 iotrack.SetV0Indexes(pt->GetV0Indexes());
535 //iotrack.SetTPCpid(pt->fTPCr);
536 //iotrack.SetTPCindex(i);
537 fEvent->AddTrack(&iotrack);
540 // short tracks - secondaties
542 if ( (pt->GetNumberOfClusters()>30) ) {
543 Int_t found,foundable,shared;
544 pt->GetClusterStatistic(128,158,found, foundable,shared,kFALSE);
545 if ( (found>20) && (pt->GetNShared()/float(pt->GetNumberOfClusters())<0.2) &&float(found)/float(foundable)>0.8){
547 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
548 iotrack.SetTPCPoints(pt->GetPoints());
549 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
550 iotrack.SetV0Indexes(pt->GetV0Indexes());
551 //iotrack.SetTPCpid(pt->fTPCr);
552 //iotrack.SetTPCindex(i);
553 fEvent->AddTrack(&iotrack);
558 if ( (pt->GetNumberOfClusters()>15)) {
559 Int_t found,foundable,shared;
560 pt->GetClusterStatistic(138,158,found, foundable,shared,kFALSE);
561 if (found<15) continue;
562 if (foundable<=0) continue;
563 if (pt->GetNShared()/float(pt->GetNumberOfClusters())>0.2) continue;
564 if (float(found)/float(foundable)<0.8) continue;
567 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
568 iotrack.SetTPCPoints(pt->GetPoints());
569 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
570 iotrack.SetV0Indexes(pt->GetV0Indexes());
571 // iotrack.SetTPCpid(pt->fTPCr);
572 //iotrack.SetTPCindex(i);
573 fEvent->AddTrack(&iotrack);
578 printf("Number of filled ESDs-\t%d\n",fEvent->GetNumberOfTracks());
585 Double_t AliTPCtrackerMI::ErrY2(AliTPCseed* seed, AliTPCclusterMI * cl){
588 // Use calibrated cluster error from OCDB
590 AliTPCClusterParam * clparam = AliTPCcalibDB::Instance()->GetClusterParam();
592 Float_t z = TMath::Abs(fParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
593 Int_t ctype = cl->GetType();
594 Int_t type = (cl->GetRow()<63) ? 0: (cl->GetRow()>126) ? 1:2;
595 Double_t angle = seed->GetSnp()*seed->GetSnp();
596 angle = TMath::Sqrt(TMath::Abs(angle/(1.-angle)));
597 Double_t erry2 = clparam->GetError0Par(0,type, z,angle);
599 erry2+=0.5; // edge cluster
602 seed->SetErrorY2(erry2);
606 //calculate look-up table at the beginning
607 // static Bool_t ginit = kFALSE;
608 // static Float_t gnoise1,gnoise2,gnoise3;
609 // static Float_t ggg1[10000];
610 // static Float_t ggg2[10000];
611 // static Float_t ggg3[10000];
612 // static Float_t glandau1[10000];
613 // static Float_t glandau2[10000];
614 // static Float_t glandau3[10000];
616 // static Float_t gcor01[500];
617 // static Float_t gcor02[500];
618 // static Float_t gcorp[500];
622 // if (ginit==kFALSE){
623 // for (Int_t i=1;i<500;i++){
624 // Float_t rsigma = float(i)/100.;
625 // gcor02[i] = TMath::Max(0.78 +TMath::Exp(7.4*(rsigma-1.2)),0.6);
626 // gcor01[i] = TMath::Max(0.72 +TMath::Exp(3.36*(rsigma-1.2)),0.6);
627 // gcorp[i] = TMath::Max(TMath::Power((rsigma+0.5),1.5),1.2);
631 // for (Int_t i=3;i<10000;i++){
635 // Float_t amp = float(i);
636 // Float_t padlength =0.75;
637 // gnoise1 = 0.0004/padlength;
638 // Float_t nel = 0.268*amp;
639 // Float_t nprim = 0.155*amp;
640 // ggg1[i] = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.001*nel/(padlength*padlength))/nel;
641 // glandau1[i] = (2.+0.12*nprim)*0.5* (2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
642 // if (glandau1[i]>1) glandau1[i]=1;
643 // glandau1[i]*=padlength*padlength/12.;
647 // gnoise2 = 0.0004/padlength;
649 // nprim = 0.133*amp;
650 // ggg2[i] = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
651 // glandau2[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
652 // if (glandau2[i]>1) glandau2[i]=1;
653 // glandau2[i]*=padlength*padlength/12.;
658 // gnoise3 = 0.0004/padlength;
660 // nprim = 0.133*amp;
661 // ggg3[i] = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
662 // glandau3[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
663 // if (glandau3[i]>1) glandau3[i]=1;
664 // glandau3[i]*=padlength*padlength/12.;
672 // Int_t amp = int(TMath::Abs(cl->GetQ()));
674 // seed->SetErrorY2(1.);
678 // Float_t z = TMath::Abs(fParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
679 // Int_t ctype = cl->GetType();
680 // Float_t padlength= GetPadPitchLength(seed->GetRow());
681 // Double_t angle2 = seed->GetSnp()*seed->GetSnp();
682 // angle2 = angle2/(1-angle2);
684 // //cluster "quality"
685 // Int_t rsigmay = int(100.*cl->GetSigmaY2()/(seed->GetCurrentSigmaY2()));
688 // if (fSectors==fInnerSec){
689 // snoise2 = gnoise1;
690 // res = ggg1[amp]*z+glandau1[amp]*angle2;
691 // if (ctype==0) res *= gcor01[rsigmay];
694 // res*= gcorp[rsigmay];
698 // if (padlength<1.1){
699 // snoise2 = gnoise2;
700 // res = ggg2[amp]*z+glandau2[amp]*angle2;
701 // if (ctype==0) res *= gcor02[rsigmay];
704 // res*= gcorp[rsigmay];
708 // snoise2 = gnoise3;
709 // res = ggg3[amp]*z+glandau3[amp]*angle2;
710 // if (ctype==0) res *= gcor02[rsigmay];
713 // res*= gcorp[rsigmay];
720 // res*=2.4; // overestimate error 2 times
724 // if (res<2*snoise2)
727 // seed->SetErrorY2(res);
735 Double_t AliTPCtrackerMI::ErrZ2(AliTPCseed* seed, AliTPCclusterMI * cl){
738 // Use calibrated cluster error from OCDB
740 AliTPCClusterParam * clparam = AliTPCcalibDB::Instance()->GetClusterParam();
742 Float_t z = TMath::Abs(fParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
743 Int_t ctype = cl->GetType();
744 Int_t type = (cl->GetRow()<63) ? 0: (cl->GetRow()>126) ? 1:2;
746 Double_t angle2 = seed->GetSnp()*seed->GetSnp();
747 angle2 = seed->GetTgl()*seed->GetTgl()*(1+angle2/(1-angle2));
748 Double_t angle = TMath::Sqrt(TMath::Abs(angle2));
749 Double_t errz2 = clparam->GetError0Par(1,type, z,angle);
751 errz2+=0.5; // edge cluster
754 seed->SetErrorZ2(errz2);
760 // //seed->SetErrorY2(0.1);
762 // //calculate look-up table at the beginning
763 // static Bool_t ginit = kFALSE;
764 // static Float_t gnoise1,gnoise2,gnoise3;
765 // static Float_t ggg1[10000];
766 // static Float_t ggg2[10000];
767 // static Float_t ggg3[10000];
768 // static Float_t glandau1[10000];
769 // static Float_t glandau2[10000];
770 // static Float_t glandau3[10000];
772 // static Float_t gcor01[1000];
773 // static Float_t gcor02[1000];
774 // static Float_t gcorp[1000];
778 // if (ginit==kFALSE){
779 // for (Int_t i=1;i<1000;i++){
780 // Float_t rsigma = float(i)/100.;
781 // gcor02[i] = TMath::Max(0.81 +TMath::Exp(6.8*(rsigma-1.2)),0.6);
782 // gcor01[i] = TMath::Max(0.72 +TMath::Exp(2.04*(rsigma-1.2)),0.6);
783 // gcorp[i] = TMath::Max(TMath::Power((rsigma+0.5),1.5),1.2);
787 // for (Int_t i=3;i<10000;i++){
791 // Float_t amp = float(i);
792 // Float_t padlength =0.75;
793 // gnoise1 = 0.0004/padlength;
794 // Float_t nel = 0.268*amp;
795 // Float_t nprim = 0.155*amp;
796 // ggg1[i] = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.001*nel/(padlength*padlength))/nel;
797 // glandau1[i] = (2.+0.12*nprim)*0.5* (2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
798 // if (glandau1[i]>1) glandau1[i]=1;
799 // glandau1[i]*=padlength*padlength/12.;
803 // gnoise2 = 0.0004/padlength;
805 // nprim = 0.133*amp;
806 // ggg2[i] = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
807 // glandau2[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
808 // if (glandau2[i]>1) glandau2[i]=1;
809 // glandau2[i]*=padlength*padlength/12.;
814 // gnoise3 = 0.0004/padlength;
816 // nprim = 0.133*amp;
817 // ggg3[i] = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
818 // glandau3[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
819 // if (glandau3[i]>1) glandau3[i]=1;
820 // glandau3[i]*=padlength*padlength/12.;
828 // Int_t amp = int(TMath::Abs(cl->GetQ()));
830 // seed->SetErrorY2(1.);
834 // Float_t z = TMath::Abs(fParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
835 // Int_t ctype = cl->GetType();
836 // Float_t padlength= GetPadPitchLength(seed->GetRow());
838 // Double_t angle2 = seed->GetSnp()*seed->GetSnp();
839 // // if (angle2<0.6) angle2 = 0.6;
840 // angle2 = seed->GetTgl()*seed->GetTgl()*(1+angle2/(1-angle2));
842 // //cluster "quality"
843 // Int_t rsigmaz = int(100.*cl->GetSigmaZ2()/(seed->GetCurrentSigmaZ2()));
846 // if (fSectors==fInnerSec){
847 // snoise2 = gnoise1;
848 // res = ggg1[amp]*z+glandau1[amp]*angle2;
849 // if (ctype==0) res *= gcor01[rsigmaz];
852 // res*= gcorp[rsigmaz];
856 // if (padlength<1.1){
857 // snoise2 = gnoise2;
858 // res = ggg2[amp]*z+glandau2[amp]*angle2;
859 // if (ctype==0) res *= gcor02[rsigmaz];
862 // res*= gcorp[rsigmaz];
866 // snoise2 = gnoise3;
867 // res = ggg3[amp]*z+glandau3[amp]*angle2;
868 // if (ctype==0) res *= gcor02[rsigmaz];
871 // res*= gcorp[rsigmaz];
880 // if ((ctype<0) &&<70){
885 // if (res<2*snoise2)
887 // if (res>3) res =3;
888 // seed->SetErrorZ2(res);
896 void AliTPCtrackerMI::RotateToLocal(AliTPCseed *seed)
898 //rotate to track "local coordinata
899 Float_t x = seed->GetX();
900 Float_t y = seed->GetY();
901 Float_t ymax = x*TMath::Tan(0.5*fSectors->GetAlpha());
904 seed->SetRelativeSector((seed->GetRelativeSector()+1) % fN);
905 if (!seed->Rotate(fSectors->GetAlpha()))
907 } else if (y <-ymax) {
908 seed->SetRelativeSector((seed->GetRelativeSector()-1+fN) % fN);
909 if (!seed->Rotate(-fSectors->GetAlpha()))
917 //_____________________________________________________________________________
918 Double_t AliTPCtrackerMI::F1old(Double_t x1,Double_t y1,
919 Double_t x2,Double_t y2,
920 Double_t x3,Double_t y3)
922 //-----------------------------------------------------------------
923 // Initial approximation of the track curvature
924 //-----------------------------------------------------------------
925 Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
926 Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
927 (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
928 Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
929 (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
931 Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
932 if ( xr*xr+yr*yr<=0.00000000000001) return 100;
933 return -xr*yr/sqrt(xr*xr+yr*yr);
938 //_____________________________________________________________________________
939 Double_t AliTPCtrackerMI::F1(Double_t x1,Double_t y1,
940 Double_t x2,Double_t y2,
941 Double_t x3,Double_t y3)
943 //-----------------------------------------------------------------
944 // Initial approximation of the track curvature
945 //-----------------------------------------------------------------
951 Double_t det = x3*y2-x2*y3;
956 Double_t u = 0.5* (x2*(x2-x3)+y2*(y2-y3))/det;
957 Double_t x0 = x3*0.5-y3*u;
958 Double_t y0 = y3*0.5+x3*u;
959 Double_t c2 = 1/TMath::Sqrt(x0*x0+y0*y0);
965 Double_t AliTPCtrackerMI::F2(Double_t x1,Double_t y1,
966 Double_t x2,Double_t y2,
967 Double_t x3,Double_t y3)
969 //-----------------------------------------------------------------
970 // Initial approximation of the track curvature
971 //-----------------------------------------------------------------
977 Double_t det = x3*y2-x2*y3;
982 Double_t u = 0.5* (x2*(x2-x3)+y2*(y2-y3))/det;
983 Double_t x0 = x3*0.5-y3*u;
984 Double_t y0 = y3*0.5+x3*u;
985 Double_t c2 = 1/TMath::Sqrt(x0*x0+y0*y0);
994 //_____________________________________________________________________________
995 Double_t AliTPCtrackerMI::F2old(Double_t x1,Double_t y1,
996 Double_t x2,Double_t y2,
997 Double_t x3,Double_t y3)
999 //-----------------------------------------------------------------
1000 // Initial approximation of the track curvature times center of curvature
1001 //-----------------------------------------------------------------
1002 Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
1003 Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
1004 (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
1005 Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
1006 (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
1008 Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
1010 return -a/(d*y1-b)*xr/sqrt(xr*xr+yr*yr);
1013 //_____________________________________________________________________________
1014 Double_t AliTPCtrackerMI::F3(Double_t x1,Double_t y1,
1015 Double_t x2,Double_t y2,
1016 Double_t z1,Double_t z2)
1018 //-----------------------------------------------------------------
1019 // Initial approximation of the tangent of the track dip angle
1020 //-----------------------------------------------------------------
1021 return (z1 - z2)/sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
1025 Double_t AliTPCtrackerMI::F3n(Double_t x1,Double_t y1,
1026 Double_t x2,Double_t y2,
1027 Double_t z1,Double_t z2, Double_t c)
1029 //-----------------------------------------------------------------
1030 // Initial approximation of the tangent of the track dip angle
1031 //-----------------------------------------------------------------
1035 //angle1 = (z1-z2)*c/(TMath::ASin(c*x1-ni)-TMath::ASin(c*x2-ni));
1037 Double_t d = TMath::Sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
1038 if (TMath::Abs(d*c*0.5)>1) return 0;
1039 // Double_t angle2 = TMath::ASin(d*c*0.5);
1040 // Double_t angle2 = AliTPCFastMath::FastAsin(d*c*0.5);
1041 Double_t angle2 = (d*c*0.5>0.1)? TMath::ASin(d*c*0.5): AliTPCFastMath::FastAsin(d*c*0.5);
1043 angle2 = (z1-z2)*c/(angle2*2.);
1047 Bool_t AliTPCtrackerMI::GetProlongation(Double_t x1, Double_t x2, Double_t x[5], Double_t &y, Double_t &z)
1048 {//-----------------------------------------------------------------
1049 // This function find proloncation of a track to a reference plane x=x2.
1050 //-----------------------------------------------------------------
1054 if (TMath::Abs(x[4]*x1 - x[2]) >= 0.999) {
1058 Double_t c1=x[4]*x1 - x[2], r1=sqrt(1.- c1*c1);
1059 Double_t c2=x[4]*x2 - x[2], r2=sqrt(1.- c2*c2);
1063 Double_t dy = dx*(c1+c2)/(r1+r2);
1066 Double_t delta = x[4]*dx*(c1+c2)/(c1*r2 + c2*r1);
1068 if (TMath::Abs(delta)>0.01){
1069 dz = x[3]*TMath::ASin(delta)/x[4];
1071 dz = x[3]*AliTPCFastMath::FastAsin(delta)/x[4];
1074 //dz = x[3]*AliTPCFastMath::FastAsin(delta)/x[4];
1082 Int_t AliTPCtrackerMI::LoadClusters (TTree *tree)
1087 return LoadClusters();
1091 Int_t AliTPCtrackerMI::LoadClusters(TObjArray *arr)
1094 // load clusters to the memory
1095 AliTPCClustersRow *clrow = 0x0;
1096 Int_t lower = arr->LowerBound();
1097 Int_t entries = arr->GetEntriesFast();
1099 for (Int_t i=lower; i<entries; i++) {
1100 clrow = (AliTPCClustersRow*) arr->At(i);
1101 if(!clrow) continue;
1102 if(!clrow->GetArray()) continue;
1106 fParam->AdjustSectorRow(clrow->GetID(),sec,row);
1108 for (Int_t icl=0; icl<clrow->GetArray()->GetEntriesFast(); icl++){
1109 Transform((AliTPCclusterMI*)(clrow->GetArray()->At(icl)));
1112 if (clrow->GetArray()->GetEntriesFast()<=0) continue;
1113 AliTPCtrackerRow * tpcrow=0;
1116 tpcrow = &(fInnerSec[sec%fkNIS][row]);
1120 tpcrow = &(fOuterSec[(sec-fkNIS*2)%fkNOS][row]);
1121 left = (sec-fkNIS*2)/fkNOS;
1124 tpcrow->SetN1(clrow->GetArray()->GetEntriesFast());
1125 tpcrow->SetClusters1(new AliTPCclusterMI[tpcrow->GetN1()]);
1126 for (Int_t i=0;i<tpcrow->GetN1();i++)
1127 tpcrow->SetCluster1(i, *(AliTPCclusterMI*)(clrow->GetArray()->At(i)));
1130 tpcrow->SetN2(clrow->GetArray()->GetEntriesFast());
1131 tpcrow->SetClusters2(new AliTPCclusterMI[tpcrow->GetN2()]);
1132 for (Int_t i=0;i<tpcrow->GetN2();i++)
1133 tpcrow->SetCluster2(i,*(AliTPCclusterMI*)(clrow->GetArray()->At(i)));
1144 Int_t AliTPCtrackerMI::LoadClusters()
1147 // load clusters to the memory
1148 AliTPCClustersRow *clrow= new AliTPCClustersRow;
1149 clrow->SetClass("AliTPCclusterMI");
1151 clrow->GetArray()->ExpandCreateFast(10000);
1153 // TTree * tree = fClustersArray.GetTree();
1155 TTree * tree = fInput;
1156 TBranch * br = tree->GetBranch("Segment");
1157 br->SetAddress(&clrow);
1159 Int_t j=Int_t(tree->GetEntries());
1160 for (Int_t i=0; i<j; i++) {
1164 fParam->AdjustSectorRow(clrow->GetID(),sec,row);
1165 for (Int_t icl=0; icl<clrow->GetArray()->GetEntriesFast(); icl++){
1166 Transform((AliTPCclusterMI*)(clrow->GetArray()->At(icl)));
1169 AliTPCtrackerRow * tpcrow=0;
1172 tpcrow = &(fInnerSec[sec%fkNIS][row]);
1176 tpcrow = &(fOuterSec[(sec-fkNIS*2)%fkNOS][row]);
1177 left = (sec-fkNIS*2)/fkNOS;
1180 tpcrow->SetN1(clrow->GetArray()->GetEntriesFast());
1181 tpcrow->SetClusters1(new AliTPCclusterMI[tpcrow->GetN1()]);
1182 for (Int_t i=0;i<tpcrow->GetN1();i++)
1183 tpcrow->SetCluster1(i, *(AliTPCclusterMI*)(clrow->GetArray()->At(i)));
1186 tpcrow->SetN2(clrow->GetArray()->GetEntriesFast());
1187 tpcrow->SetClusters2(new AliTPCclusterMI[tpcrow->GetN2()]);
1188 for (Int_t i=0;i<tpcrow->GetN2();i++)
1189 tpcrow->SetCluster2(i,*(AliTPCclusterMI*)(clrow->GetArray()->At(i)));
1200 void AliTPCtrackerMI::UnloadClusters()
1203 // unload clusters from the memory
1205 Int_t nrows = fOuterSec->GetNRows();
1206 for (Int_t sec = 0;sec<fkNOS;sec++)
1207 for (Int_t row = 0;row<nrows;row++){
1208 AliTPCtrackerRow* tpcrow = &(fOuterSec[sec%fkNOS][row]);
1210 // if (tpcrow->fClusters1) delete []tpcrow->fClusters1;
1211 // if (tpcrow->fClusters2) delete []tpcrow->fClusters2;
1213 tpcrow->ResetClusters();
1216 nrows = fInnerSec->GetNRows();
1217 for (Int_t sec = 0;sec<fkNIS;sec++)
1218 for (Int_t row = 0;row<nrows;row++){
1219 AliTPCtrackerRow* tpcrow = &(fInnerSec[sec%fkNIS][row]);
1221 // if (tpcrow->fClusters1) delete []tpcrow->fClusters1;
1222 //if (tpcrow->fClusters2) delete []tpcrow->fClusters2;
1224 tpcrow->ResetClusters();
1230 void AliTPCtrackerMI::FillClusterArray(TObjArray* array) const{
1232 // Filling cluster to the array - For visualization purposes
1235 nrows = fOuterSec->GetNRows();
1236 for (Int_t sec = 0;sec<fkNOS;sec++)
1237 for (Int_t row = 0;row<nrows;row++){
1238 AliTPCtrackerRow* tpcrow = &(fOuterSec[sec%fkNOS][row]);
1239 if (!tpcrow) continue;
1240 for (Int_t icl = 0;icl<tpcrow->GetN();icl++){
1241 array->AddLast((TObject*)((*tpcrow)[icl]));
1244 nrows = fInnerSec->GetNRows();
1245 for (Int_t sec = 0;sec<fkNIS;sec++)
1246 for (Int_t row = 0;row<nrows;row++){
1247 AliTPCtrackerRow* tpcrow = &(fInnerSec[sec%fkNIS][row]);
1248 if (!tpcrow) continue;
1249 for (Int_t icl = 0;icl<tpcrow->GetN();icl++){
1250 array->AddLast((TObject*)(*tpcrow)[icl]);
1256 void AliTPCtrackerMI::Transform(AliTPCclusterMI * cluster){
1260 AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ;
1262 AliFatal("Tranformations not in calibDB");
1264 Double_t x[3]={cluster->GetRow(),cluster->GetPad(),cluster->GetTimeBin()};
1265 Int_t i[1]={cluster->GetDetector()};
1266 transform->Transform(x,i,0,1);
1267 if (!AliTPCReconstructor::GetRecoParam()->GetBYMirror()){
1268 if (cluster->GetDetector()%36>17){
1274 // in debug mode check the transformation
1276 if (AliTPCReconstructor::StreamLevel()>1) {
1278 cluster->GetGlobalXYZ(gx);
1279 Int_t event = (fEvent==NULL)? 0: fEvent->GetEventNumberInFile();
1280 TTreeSRedirector &cstream = *fDebugStreamer;
1281 cstream<<"Transform"<<
1292 cluster->SetX(x[0]);
1293 cluster->SetY(x[1]);
1294 cluster->SetZ(x[2]);
1300 //if (!fParam->IsGeoRead()) fParam->ReadGeoMatrices();
1301 TGeoHMatrix *mat = fParam->GetClusterMatrix(cluster->GetDetector());
1303 Double_t pos[3]= {cluster->GetX(),cluster->GetY(),cluster->GetZ()};
1304 Double_t posC[3]={cluster->GetX(),cluster->GetY(),cluster->GetZ()};
1305 if (mat) mat->LocalToMaster(pos,posC);
1307 // chack Loading of Geo matrices from GeoManager - TEMPORARY FIX
1309 cluster->SetX(posC[0]);
1310 cluster->SetY(posC[1]);
1311 cluster->SetZ(posC[2]);
1314 //_____________________________________________________________________________
1315 Int_t AliTPCtrackerMI::LoadOuterSectors() {
1316 //-----------------------------------------------------------------
1317 // This function fills outer TPC sectors with clusters.
1318 //-----------------------------------------------------------------
1319 Int_t nrows = fOuterSec->GetNRows();
1321 for (Int_t sec = 0;sec<fkNOS;sec++)
1322 for (Int_t row = 0;row<nrows;row++){
1323 AliTPCtrackerRow* tpcrow = &(fOuterSec[sec%fkNOS][row]);
1324 Int_t sec2 = sec+2*fkNIS;
1326 Int_t ncl = tpcrow->GetN1();
1328 AliTPCclusterMI *c= (tpcrow->GetCluster1(ncl));
1329 index=(((sec2<<8)+row)<<16)+ncl;
1330 tpcrow->InsertCluster(c,index);
1333 ncl = tpcrow->GetN2();
1335 AliTPCclusterMI *c= (tpcrow->GetCluster2(ncl));
1336 index=((((sec2+fkNOS)<<8)+row)<<16)+ncl;
1337 tpcrow->InsertCluster(c,index);
1340 // write indexes for fast acces
1342 for (Int_t i=0;i<510;i++)
1343 tpcrow->SetFastCluster(i,-1);
1344 for (Int_t i=0;i<tpcrow->GetN();i++){
1345 Int_t zi = Int_t((*tpcrow)[i]->GetZ()+255.);
1346 tpcrow->SetFastCluster(zi,i); // write index
1349 for (Int_t i=0;i<510;i++){
1350 if (tpcrow->GetFastCluster(i)<0)
1351 tpcrow->SetFastCluster(i,last);
1353 last = tpcrow->GetFastCluster(i);
1362 //_____________________________________________________________________________
1363 Int_t AliTPCtrackerMI::LoadInnerSectors() {
1364 //-----------------------------------------------------------------
1365 // This function fills inner TPC sectors with clusters.
1366 //-----------------------------------------------------------------
1367 Int_t nrows = fInnerSec->GetNRows();
1369 for (Int_t sec = 0;sec<fkNIS;sec++)
1370 for (Int_t row = 0;row<nrows;row++){
1371 AliTPCtrackerRow* tpcrow = &(fInnerSec[sec%fkNIS][row]);
1374 Int_t ncl = tpcrow->GetN1();
1376 AliTPCclusterMI *c= (tpcrow->GetCluster1(ncl));
1377 index=(((sec<<8)+row)<<16)+ncl;
1378 tpcrow->InsertCluster(c,index);
1381 ncl = tpcrow->GetN2();
1383 AliTPCclusterMI *c= (tpcrow->GetCluster2(ncl));
1384 index=((((sec+fkNIS)<<8)+row)<<16)+ncl;
1385 tpcrow->InsertCluster(c,index);
1388 // write indexes for fast acces
1390 for (Int_t i=0;i<510;i++)
1391 tpcrow->SetFastCluster(i,-1);
1392 for (Int_t i=0;i<tpcrow->GetN();i++){
1393 Int_t zi = Int_t((*tpcrow)[i]->GetZ()+255.);
1394 tpcrow->SetFastCluster(zi,i); // write index
1397 for (Int_t i=0;i<510;i++){
1398 if (tpcrow->GetFastCluster(i)<0)
1399 tpcrow->SetFastCluster(i,last);
1401 last = tpcrow->GetFastCluster(i);
1413 //_________________________________________________________________________
1414 AliTPCclusterMI *AliTPCtrackerMI::GetClusterMI(Int_t index) const {
1415 //--------------------------------------------------------------------
1416 // Return pointer to a given cluster
1417 //--------------------------------------------------------------------
1418 if (index<0) return 0; // no cluster
1419 Int_t sec=(index&0xff000000)>>24;
1420 Int_t row=(index&0x00ff0000)>>16;
1421 Int_t ncl=(index&0x00007fff)>>00;
1423 const AliTPCtrackerRow * tpcrow=0;
1424 AliTPCclusterMI * clrow =0;
1426 if (sec<0 || sec>=fkNIS*4) {
1427 AliWarning(Form("Wrong sector %d",sec));
1432 tpcrow = &(fInnerSec[sec%fkNIS][row]);
1433 if (tpcrow==0) return 0;
1436 if (tpcrow->GetN1()<=ncl) return 0;
1437 clrow = tpcrow->GetClusters1();
1440 if (tpcrow->GetN2()<=ncl) return 0;
1441 clrow = tpcrow->GetClusters2();
1445 tpcrow = &(fOuterSec[(sec-fkNIS*2)%fkNOS][row]);
1446 if (tpcrow==0) return 0;
1448 if (sec-2*fkNIS<fkNOS) {
1449 if (tpcrow->GetN1()<=ncl) return 0;
1450 clrow = tpcrow->GetClusters1();
1453 if (tpcrow->GetN2()<=ncl) return 0;
1454 clrow = tpcrow->GetClusters2();
1458 return &(clrow[ncl]);
1464 Int_t AliTPCtrackerMI::FollowToNext(AliTPCseed& t, Int_t nr) {
1465 //-----------------------------------------------------------------
1466 // This function tries to find a track prolongation to next pad row
1467 //-----------------------------------------------------------------
1469 Double_t x= GetXrow(nr), ymax=GetMaxY(nr);
1470 AliTPCclusterMI *cl=0;
1471 Int_t tpcindex= t.GetClusterIndex2(nr);
1473 // update current shape info every 5 pad-row
1474 // if ( (nr%5==0) || t.GetNumberOfClusters()<2 || (t.fCurrentSigmaY2<0.0001) ){
1478 if (fIteration>0 && tpcindex>=-1){ //if we have already clusters
1480 if (tpcindex==-1) return 0; //track in dead zone
1482 cl = t.GetClusterPointer(nr);
1483 if ( (cl==0) ) cl = GetClusterMI(tpcindex);
1484 t.SetCurrentClusterIndex1(tpcindex);
1487 Int_t relativesector = ((tpcindex&0xff000000)>>24)%18; // if previously accepted cluster in different sector
1488 Float_t angle = relativesector*fSectors->GetAlpha()+fSectors->GetAlphaShift();
1490 if (angle<-TMath::Pi()) angle += 2*TMath::Pi();
1491 if (angle>=TMath::Pi()) angle -= 2*TMath::Pi();
1493 if (TMath::Abs(angle-t.GetAlpha())>0.001){
1494 Double_t rotation = angle-t.GetAlpha();
1495 t.SetRelativeSector(relativesector);
1496 if (!t.Rotate(rotation)) return 0;
1498 if (!t.PropagateTo(x)) return 0;
1500 t.SetCurrentCluster(cl);
1502 Int_t accept = AcceptCluster(&t,t.GetCurrentCluster());
1503 if ((tpcindex&0x8000)==0) accept =0;
1505 //if founded cluster is acceptible
1506 if (cl->IsUsed(11)) { // id cluster is shared inrease uncertainty
1507 t.SetErrorY2(t.GetErrorY2()+0.03);
1508 t.SetErrorZ2(t.GetErrorZ2()+0.03);
1509 t.SetErrorY2(t.GetErrorY2()*3);
1510 t.SetErrorZ2(t.GetErrorZ2()*3);
1512 t.SetNFoundable(t.GetNFoundable()+1);
1513 UpdateTrack(&t,accept);
1518 if (TMath::Abs(t.GetSnp())>AliTPCReconstructor::GetMaxSnpTracker()) return 0; // cut on angle
1520 // not look for new cluster during refitting
1521 t.SetNFoundable(t.GetNFoundable()+1);
1526 // if (TMath::Abs(t.GetSnp())>0.95 || TMath::Abs(x*t.GetC()-t.GetEta())>0.95) return 0;// patch 28 fev 06
1527 Double_t y=t.GetYat(x);
1528 if (TMath::Abs(y)>ymax){
1530 t.SetRelativeSector((t.GetRelativeSector()+1) % fN);
1531 if (!t.Rotate(fSectors->GetAlpha()))
1533 } else if (y <-ymax) {
1534 t.SetRelativeSector((t.GetRelativeSector()-1+fN) % fN);
1535 if (!t.Rotate(-fSectors->GetAlpha()))
1541 if (!t.PropagateTo(x)) {
1542 if (fIteration==0) t.SetRemoval(10);
1546 Double_t z=t.GetZ();
1549 if (!IsActive(t.GetRelativeSector(),nr)) {
1551 t.SetClusterIndex2(nr,-1);
1554 //AliInfo(Form("A - Sector%d phi %f - alpha %f", t.fRelativeSector,y/x, t.GetAlpha()));
1555 Bool_t isActive = IsActive(t.GetRelativeSector(),nr);
1556 Bool_t isActive2 = (nr>=fInnerSec->GetNRows()) ? fOuterSec[t.GetRelativeSector()][nr-fInnerSec->GetNRows()].GetN()>0:fInnerSec[t.GetRelativeSector()][nr].GetN()>0;
1558 if (!isActive || !isActive2) return 0;
1560 const AliTPCtrackerRow &krow=GetRow(t.GetRelativeSector(),nr);
1561 if ( (t.GetSigmaY2()<0) || t.GetSigmaZ2()<0) return 0;
1563 Double_t roadz = 1.;
1565 if (TMath::Abs(TMath::Abs(y)-ymax)<krow.GetDeadZone()){
1567 t.SetClusterIndex2(nr,-1);
1572 if (TMath::Abs(z)<(AliTPCReconstructor::GetCtgRange()*x+10) && TMath::Abs(z)<fParam->GetZLength(0) && (TMath::Abs(t.GetSnp())<AliTPCReconstructor::GetMaxSnpTracker()))
1573 t.SetNFoundable(t.GetNFoundable()+1);
1579 // cl = krow.FindNearest2(y+10.,z,roady,roadz,index);
1580 cl = krow.FindNearest2(y,z,roady,roadz,index);
1581 if (cl) t.SetCurrentClusterIndex1(krow.GetIndex(index));
1584 t.SetCurrentCluster(cl);
1586 if (fIteration==2&&cl->IsUsed(10)) return 0;
1587 Int_t accept = AcceptCluster(&t,t.GetCurrentCluster());
1588 if (fIteration==2&&cl->IsUsed(11)) {
1589 t.SetErrorY2(t.GetErrorY2()+0.03);
1590 t.SetErrorZ2(t.GetErrorZ2()+0.03);
1591 t.SetErrorY2(t.GetErrorY2()*3);
1592 t.SetErrorZ2(t.GetErrorZ2()*3);
1595 if (t.fCurrentCluster->IsUsed(10)){
1600 if (t.fNShared>0.7*t.GetNumberOfClusters()) {
1606 if (accept<3) UpdateTrack(&t,accept);
1609 if ( fIteration==0 && t.GetNFoundable()*0.5 > t.GetNumberOfClusters()) t.SetRemoval(10);
1615 Int_t AliTPCtrackerMI::FollowToNextFast(AliTPCseed& t, Int_t nr) {
1616 //-----------------------------------------------------------------
1617 // This function tries to find a track prolongation to next pad row
1618 //-----------------------------------------------------------------
1620 Double_t x= GetXrow(nr), ymax=GetMaxY(nr);
1622 if (!t.GetProlongation(x,y,z)) {
1628 if (TMath::Abs(y)>ymax){
1631 t.SetRelativeSector((t.GetRelativeSector()+1) % fN);
1632 if (!t.Rotate(fSectors->GetAlpha()))
1634 } else if (y <-ymax) {
1635 t.SetRelativeSector((t.GetRelativeSector()-1+fN) % fN);
1636 if (!t.Rotate(-fSectors->GetAlpha()))
1639 if (!t.PropagateTo(x)) {
1642 t.GetProlongation(x,y,z);
1645 // update current shape info every 2 pad-row
1646 if ( (nr%2==0) || t.GetNumberOfClusters()<2 || (t.GetCurrentSigmaY2()<0.0001) ){
1647 // t.fCurrentSigmaY = GetSigmaY(&t);
1648 //t.fCurrentSigmaZ = GetSigmaZ(&t);
1652 AliTPCclusterMI *cl=0;
1657 const AliTPCtrackerRow &krow=GetRow(t.GetRelativeSector(),nr);
1658 if ( (t.GetSigmaY2()<0) || t.GetSigmaZ2()<0) return 0;
1660 Double_t roadz = 1.;
1663 if (TMath::Abs(TMath::Abs(y)-ymax)<krow.GetDeadZone()){
1665 t.SetClusterIndex2(row,-1);
1670 if (TMath::Abs(z)>(AliTPCReconstructor::GetCtgRange()*x+10)) t.SetClusterIndex2(row,-1);
1674 if ((cl==0)&&(krow)) {
1675 // cl = krow.FindNearest2(y+10,z,roady,roadz,index);
1676 cl = krow.FindNearest2(y,z,roady,roadz,index);
1678 if (cl) t.SetCurrentClusterIndex1(krow.GetIndex(index));
1682 t.SetCurrentCluster(cl);
1683 // Int_t accept = AcceptCluster(&t,t.fCurrentCluster);
1685 t.SetClusterIndex2(row,index);
1686 t.SetClusterPointer(row, cl);
1694 //_________________________________________________________________________
1695 Bool_t AliTPCtrackerMI::GetTrackPoint(Int_t index, AliTrackPoint &p ) const
1697 // Get track space point by index
1698 // return false in case the cluster doesn't exist
1699 AliTPCclusterMI *cl = GetClusterMI(index);
1700 if (!cl) return kFALSE;
1701 Int_t sector = (index&0xff000000)>>24;
1702 // Int_t row = (index&0x00ff0000)>>16;
1704 // xyz[0] = fParam->GetPadRowRadii(sector,row);
1705 xyz[0] = cl->GetX();
1706 xyz[1] = cl->GetY();
1707 xyz[2] = cl->GetZ();
1709 fParam->AdjustCosSin(sector,cos,sin);
1710 Float_t x = cos*xyz[0]-sin*xyz[1];
1711 Float_t y = cos*xyz[1]+sin*xyz[0];
1713 Float_t sigmaY2 = 0.027*cl->GetSigmaY2();
1714 if (sector < fParam->GetNInnerSector()) sigmaY2 *= 2.07;
1715 Float_t sigmaZ2 = 0.066*cl->GetSigmaZ2();
1716 if (sector < fParam->GetNInnerSector()) sigmaZ2 *= 1.77;
1717 cov[0] = sin*sin*sigmaY2;
1718 cov[1] = -sin*cos*sigmaY2;
1720 cov[3] = cos*cos*sigmaY2;
1723 p.SetXYZ(x,y,xyz[2],cov);
1724 AliGeomManager::ELayerID iLayer;
1726 if (sector < fParam->GetNInnerSector()) {
1727 iLayer = AliGeomManager::kTPC1;
1731 iLayer = AliGeomManager::kTPC2;
1732 idet = sector - fParam->GetNInnerSector();
1734 UShort_t volid = AliGeomManager::LayerToVolUID(iLayer,idet);
1735 p.SetVolumeID(volid);
1741 Int_t AliTPCtrackerMI::UpdateClusters(AliTPCseed& t, Int_t nr) {
1742 //-----------------------------------------------------------------
1743 // This function tries to find a track prolongation to next pad row
1744 //-----------------------------------------------------------------
1745 t.SetCurrentCluster(0);
1746 t.SetCurrentClusterIndex1(0);
1748 Double_t xt=t.GetX();
1749 Int_t row = GetRowNumber(xt)-1;
1750 Double_t ymax= GetMaxY(nr);
1752 if (row < nr) return 1; // don't prolongate if not information until now -
1753 // if (TMath::Abs(t.GetSnp())>0.9 && t.GetNumberOfClusters()>40. && fIteration!=2) {
1755 // return 0; // not prolongate strongly inclined tracks
1757 // if (TMath::Abs(t.GetSnp())>0.95) {
1759 // return 0; // not prolongate strongly inclined tracks
1760 // }// patch 28 fev 06
1762 Double_t x= GetXrow(nr);
1764 //t.PropagateTo(x+0.02);
1765 //t.PropagateTo(x+0.01);
1766 if (!t.PropagateTo(x)){
1773 if (TMath::Abs(y)>ymax){
1775 t.SetRelativeSector((t.GetRelativeSector()+1) % fN);
1776 if (!t.Rotate(fSectors->GetAlpha()))
1778 } else if (y <-ymax) {
1779 t.SetRelativeSector((t.GetRelativeSector()-1+fN) % fN);
1780 if (!t.Rotate(-fSectors->GetAlpha()))
1783 // if (!t.PropagateTo(x)){
1790 if (TMath::Abs(t.GetSnp())>AliTPCReconstructor::GetMaxSnpTracker()) return 0;
1792 if (!IsActive(t.GetRelativeSector(),nr)) {
1794 t.SetClusterIndex2(nr,-1);
1797 //AliInfo(Form("A - Sector%d phi %f - alpha %f", t.fRelativeSector,y/x, t.GetAlpha()));
1799 AliTPCtrackerRow &krow=GetRow(t.GetRelativeSector(),nr);
1801 if (TMath::Abs(TMath::Abs(y)-ymax)<krow.GetDeadZone()){
1803 t.SetClusterIndex2(nr,-1);
1808 if (TMath::Abs(t.GetZ())<(AliTPCReconstructor::GetCtgRange()*t.GetX()+10) && (TMath::Abs(t.GetSnp())<AliTPCReconstructor::GetMaxSnpTracker())) t.SetNFoundable(t.GetNFoundable()+1);
1814 if ( (nr%2==0) || t.GetNumberOfClusters()<2){
1815 // t.fCurrentSigmaY = GetSigmaY(&t);
1816 //t.fCurrentSigmaZ = GetSigmaZ(&t);
1820 AliTPCclusterMI *cl=0;
1823 Double_t roady = 1.;
1824 Double_t roadz = 1.;
1828 index = t.GetClusterIndex2(nr);
1829 if ( (index>0) && (index&0x8000)==0){
1830 cl = t.GetClusterPointer(nr);
1831 if ( (cl==0) && (index>0)) cl = GetClusterMI(index);
1832 t.SetCurrentClusterIndex1(index);
1834 t.SetCurrentCluster(cl);
1840 // if (index<0) return 0;
1841 UInt_t uindex = TMath::Abs(index);
1844 //cl = krow.FindNearest2(y+10,z,roady,roadz,uindex);
1845 cl = krow.FindNearest2(y,z,roady,roadz,uindex);
1848 if (cl) t.SetCurrentClusterIndex1(krow.GetIndex(uindex));
1849 t.SetCurrentCluster(cl);
1855 Int_t AliTPCtrackerMI::FollowToNextCluster(AliTPCseed & t, Int_t nr) {
1856 //-----------------------------------------------------------------
1857 // This function tries to find a track prolongation to next pad row
1858 //-----------------------------------------------------------------
1860 //update error according neighborhoud
1862 if (t.GetCurrentCluster()) {
1864 Int_t accept = AcceptCluster(&t,t.GetCurrentCluster());
1866 if (t.GetCurrentCluster()->IsUsed(10)){
1871 t.SetNShared(t.GetNShared()+1);
1872 if (t.GetNShared()>0.7*t.GetNumberOfClusters()) {
1877 if (fIteration>0) accept = 0;
1878 if (accept<3) UpdateTrack(&t,accept);
1882 if ( ( (t.GetSigmaY2()+t.GetSigmaZ2())>0.16)&& t.GetNumberOfClusters()>18) t.SetRemoval(10);
1883 if ( t.GetChi2()/t.GetNumberOfClusters()>6 &&t.GetNumberOfClusters()>18) t.SetRemoval(10);
1885 if (( (t.GetNFoundable()*0.5 > t.GetNumberOfClusters()) || t.GetNoCluster()>15)) t.SetRemoval(10);
1893 //_____________________________________________________________________________
1894 Int_t AliTPCtrackerMI::FollowProlongation(AliTPCseed& t, Int_t rf, Int_t step) {
1895 //-----------------------------------------------------------------
1896 // This function tries to find a track prolongation.
1897 //-----------------------------------------------------------------
1898 Double_t xt=t.GetX();
1900 Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
1901 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
1902 if (alpha < 0. ) alpha += 2.*TMath::Pi();
1904 t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN);
1906 Int_t first = GetRowNumber(xt)-1;
1907 for (Int_t nr= first; nr>=rf; nr-=step) {
1909 if (t.GetKinkIndexes()[0]>0){
1910 for (Int_t i=0;i<3;i++){
1911 Int_t index = t.GetKinkIndexes()[i];
1912 if (index==0) break;
1913 if (index<0) continue;
1915 AliKink * kink = (AliKink*)fEvent->GetKink(index-1);
1917 printf("PROBLEM\n");
1920 Int_t kinkrow = kink->GetTPCRow0()+2+Int_t(0.5/(0.05+kink->GetAngle(2)));
1922 AliExternalTrackParam paramd(t);
1923 kink->SetDaughter(paramd);
1924 kink->SetStatus(2,5);
1931 if (nr==80) t.UpdateReference();
1932 if (nr<fInnerSec->GetNRows())
1933 fSectors = fInnerSec;
1935 fSectors = fOuterSec;
1936 if (FollowToNext(t,nr)==0)
1945 //_____________________________________________________________________________
1946 Int_t AliTPCtrackerMI::FollowProlongationFast(AliTPCseed& t, Int_t rf, Int_t step) {
1947 //-----------------------------------------------------------------
1948 // This function tries to find a track prolongation.
1949 //-----------------------------------------------------------------
1950 Double_t xt=t.GetX();
1952 Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
1953 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
1954 if (alpha < 0. ) alpha += 2.*TMath::Pi();
1955 t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN);
1957 for (Int_t nr=GetRowNumber(xt)-1; nr>=rf; nr-=step) {
1959 if (FollowToNextFast(t,nr)==0)
1960 if (!t.IsActive()) return 0;
1970 Int_t AliTPCtrackerMI::FollowBackProlongation(AliTPCseed& t, Int_t rf) {
1971 //-----------------------------------------------------------------
1972 // This function tries to find a track prolongation.
1973 //-----------------------------------------------------------------
1975 Double_t xt=t.GetX();
1976 Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
1977 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
1978 if (alpha < 0. ) alpha += 2.*TMath::Pi();
1979 t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN);
1981 Int_t first = t.GetFirstPoint();
1982 if (first<GetRowNumber(xt)+1) first = GetRowNumber(xt)+1;
1984 if (first<0) first=0;
1985 for (Int_t nr=first; nr<=rf; nr++) {
1986 // if ( (TMath::Abs(t.GetSnp())>0.95)) break;//patch 28 fev 06
1987 if (t.GetKinkIndexes()[0]<0){
1988 for (Int_t i=0;i<3;i++){
1989 Int_t index = t.GetKinkIndexes()[i];
1990 if (index==0) break;
1991 if (index>0) continue;
1992 index = TMath::Abs(index);
1993 AliKink * kink = (AliKink*)fEvent->GetKink(index-1);
1995 printf("PROBLEM\n");
1998 Int_t kinkrow = kink->GetTPCRow0()-2-Int_t(0.5/(0.05+kink->GetAngle(2)));
2000 AliExternalTrackParam paramm(t);
2001 kink->SetMother(paramm);
2002 kink->SetStatus(2,1);
2009 if (nr<fInnerSec->GetNRows())
2010 fSectors = fInnerSec;
2012 fSectors = fOuterSec;
2023 Float_t AliTPCtrackerMI::OverlapFactor(AliTPCseed * s1, AliTPCseed * s2, Int_t &sum1, Int_t & sum2)
2031 Float_t dz2 =(s1->GetZ() - s2->GetZ());
2034 Float_t dy2 =TMath::Abs((s1->GetY() - s2->GetY()));
2036 Float_t distance = TMath::Sqrt(dz2+dy2);
2037 if (distance>4.) return 0; // if there are far away - not overlap - to reduce combinatorics
2040 Int_t firstpoint = TMath::Min(s1->GetFirstPoint(),s2->GetFirstPoint());
2041 Int_t lastpoint = TMath::Max(s1->GetLastPoint(),s2->GetLastPoint());
2046 if (firstpoint>lastpoint) {
2047 firstpoint =lastpoint;
2052 for (Int_t i=firstpoint-1;i<lastpoint+1;i++){
2053 if (s1->GetClusterIndex2(i)>0) sum1++;
2054 if (s2->GetClusterIndex2(i)>0) sum2++;
2055 if (s1->GetClusterIndex2(i)==s2->GetClusterIndex2(i) && s1->GetClusterIndex2(i)>0) {
2059 if (sum<5) return 0;
2061 Float_t summin = TMath::Min(sum1+1,sum2+1);
2062 Float_t ratio = (sum+1)/Float_t(summin);
2066 void AliTPCtrackerMI::SignShared(AliTPCseed * s1, AliTPCseed * s2)
2070 Float_t thetaCut = 0.2;//+10.*TMath::Sqrt(s1->GetSigmaTglZ()+ s2->GetSigmaTglZ());
2071 if (TMath::Abs(s1->GetTgl()-s2->GetTgl())>thetaCut) return;
2072 Float_t minCl = TMath::Min(s1->GetNumberOfClusters(),s2->GetNumberOfClusters());
2073 Int_t cutN0 = TMath::Max(5,TMath::Nint(0.1*minCl));
2078 //Int_t firstpoint = TMath::Max(s1->GetFirstPoint(),s2->GetFirstPoint());
2079 //Int_t lastpoint = TMath::Min(s1->GetLastPoint(),s2->GetLastPoint());
2080 Int_t firstpoint = 0;
2081 Int_t lastpoint = 160;
2083 // if (firstpoint>=lastpoint-5) return;;
2085 for (Int_t i=firstpoint;i<lastpoint;i++){
2086 // if ( (s1->GetClusterIndex2(i)&0xFFFF8FFF)==(s2->GetClusterIndex2(i)&0xFFFF8FFF) && s1->GetClusterIndex2(i)>0) {
2087 if ( (s1->GetClusterIndex2(i))==(s2->GetClusterIndex2(i)) && s1->GetClusterIndex2(i)>0) {
2089 s1->SetSharedMapBit(i, kTRUE);
2090 s2->SetSharedMapBit(i, kTRUE);
2092 if (s1->GetClusterIndex2(i)>0)
2093 s1->SetClusterMapBit(i, kTRUE);
2095 if (sumshared>cutN0){
2098 for (Int_t i=firstpoint;i<lastpoint;i++){
2099 // if ( (s1->GetClusterIndex2(i)&0xFFFF8FFF)==(s2->GetClusterIndex2(i)&0xFFFF8FFF) && s1->GetClusterIndex2(i)>0) {
2100 if ( (s1->GetClusterIndex2(i))==(s2->GetClusterIndex2(i)) && s1->GetClusterIndex2(i)>0) {
2101 AliTPCTrackerPoint *p1 = s1->GetTrackPoint(i);
2102 AliTPCTrackerPoint *p2 = s2->GetTrackPoint(i);;
2103 if (s1->IsActive()&&s2->IsActive()){
2104 p1->SetShared(kTRUE);
2105 p2->SetShared(kTRUE);
2111 if (sumshared>cutN0){
2112 for (Int_t i=0;i<4;i++){
2113 if (s1->GetOverlapLabel(3*i)==0){
2114 s1->SetOverlapLabel(3*i, s2->GetLabel());
2115 s1->SetOverlapLabel(3*i+1,sumshared);
2116 s1->SetOverlapLabel(3*i+2,s2->GetUniqueID());
2120 for (Int_t i=0;i<4;i++){
2121 if (s2->GetOverlapLabel(3*i)==0){
2122 s2->SetOverlapLabel(3*i, s1->GetLabel());
2123 s2->SetOverlapLabel(3*i+1,sumshared);
2124 s2->SetOverlapLabel(3*i+2,s1->GetUniqueID());
2131 void AliTPCtrackerMI::SignShared(TObjArray * arr)
2134 //sort trackss according sectors
2136 for (Int_t i=0; i<arr->GetEntriesFast(); i++) {
2137 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2139 //if (pt) RotateToLocal(pt);
2143 arr->Sort(); // sorting according relative sectors
2144 arr->Expand(arr->GetEntries());
2147 Int_t nseed=arr->GetEntriesFast();
2148 for (Int_t i=0; i<nseed; i++) {
2149 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2151 for (Int_t j=0;j<=12;j++){
2152 pt->SetOverlapLabel(j,0);
2155 for (Int_t i=0; i<nseed; i++) {
2156 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2158 if (pt->GetRemoval()>10) continue;
2159 for (Int_t j=i+1; j<nseed; j++){
2160 AliTPCseed *pt2=(AliTPCseed*)arr->UncheckedAt(j);
2161 if (TMath::Abs(pt->GetRelativeSector()-pt2->GetRelativeSector())>1) continue;
2163 if (pt2->GetRemoval()<=10) {
2164 //if ( TMath::Abs(pt->GetRelativeSector()-pt2->GetRelativeSector())>0) break;
2172 void AliTPCtrackerMI::SortTracks(TObjArray * arr, Int_t mode) const
2175 //sort tracks in array according mode criteria
2176 Int_t nseed = arr->GetEntriesFast();
2177 for (Int_t i=0; i<nseed; i++) {
2178 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2189 void AliTPCtrackerMI::RemoveUsed2(TObjArray * arr, Float_t factor1, Float_t factor2, Int_t minimal)
2192 // Loop over all tracks and remove overlaped tracks (with lower quality)
2194 // 1. Unsign clusters
2195 // 2. Sort tracks according quality
2196 // Quality is defined by the number of cluster between first and last points
2198 // 3. Loop over tracks - decreasing quality order
2199 // a.) remove - If the fraction of shared cluster less than factor (1- n or 2)
2200 // b.) remove - If the minimal number of clusters less than minimal and not ITS
2201 // c.) if track accepted - sign clusters
2203 //Called in - AliTPCtrackerMI::Clusters2Tracks()
2204 // - AliTPCtrackerMI::PropagateBack()
2205 // - AliTPCtrackerMI::RefitInward()
2211 Int_t nseed = arr->GetEntriesFast();
2212 Float_t * quality = new Float_t[nseed];
2213 Int_t * indexes = new Int_t[nseed];
2217 for (Int_t i=0; i<nseed; i++) {
2218 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2223 pt->UpdatePoints(); //select first last max dens points
2224 Float_t * points = pt->GetPoints();
2225 if (points[3]<0.8) quality[i] =-1;
2226 quality[i] = (points[2]-points[0])+pt->GetNumberOfClusters();
2227 //prefer high momenta tracks if overlaps
2228 quality[i] *= TMath::Sqrt(TMath::Abs(pt->Pt())+0.5);
2230 TMath::Sort(nseed,quality,indexes);
2233 for (Int_t itrack=0; itrack<nseed; itrack++) {
2234 Int_t trackindex = indexes[itrack];
2235 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(trackindex);
2238 if (quality[trackindex]<0){
2240 delete arr->RemoveAt(trackindex);
2243 arr->RemoveAt(trackindex);
2249 Int_t first = Int_t(pt->GetPoints()[0]);
2250 Int_t last = Int_t(pt->GetPoints()[2]);
2251 Double_t factor = (pt->GetBConstrain()) ? factor1: factor2;
2253 Int_t found,foundable,shared;
2254 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
2255 // pt->GetClusterStatistic(0,160, found, foundable,shared,kFALSE);
2256 Bool_t itsgold =kFALSE;
2259 if (pt->GetESD()->GetITSclusters(dummy)>4) itsgold= kTRUE;
2263 if (Float_t(shared+1)/Float_t(found+1)>factor){
2264 if (pt->GetKinkIndexes()[0]!=0) continue; //don't remove tracks - part of the kinks
2265 delete arr->RemoveAt(trackindex);
2268 if (pt->GetNumberOfClusters()<50&&(found-0.5*shared)<minimal){ //remove short tracks
2269 if (pt->GetKinkIndexes()[0]!=0) continue; //don't remove tracks - part of the kinks
2270 delete arr->RemoveAt(trackindex);
2276 //if (sharedfactor>0.4) continue;
2277 if (pt->GetKinkIndexes()[0]>0) continue;
2278 //Remove tracks with undefined properties - seems
2279 if (pt->GetSigmaY2()<kAlmost0) continue; // ? what is the origin ?
2281 for (Int_t i=first; i<last; i++) {
2282 Int_t index=pt->GetClusterIndex2(i);
2283 // if (index<0 || index&0x8000 ) continue;
2284 if (index<0 || index&0x8000 ) continue;
2285 AliTPCclusterMI *c= pt->GetClusterPointer(i);
2292 Info("RemoveUsed","\n*****\nNumber of good tracks after shared removal\t%d\n",fNtracks);
2298 void AliTPCtrackerMI::UnsignClusters()
2301 // loop over all clusters and unsign them
2304 for (Int_t sec=0;sec<fkNIS;sec++){
2305 for (Int_t row=0;row<fInnerSec->GetNRows();row++){
2306 AliTPCclusterMI *cl = fInnerSec[sec][row].GetClusters1();
2307 for (Int_t icl =0;icl< fInnerSec[sec][row].GetN1();icl++)
2308 // if (cl[icl].IsUsed(10))
2310 cl = fInnerSec[sec][row].GetClusters2();
2311 for (Int_t icl =0;icl< fInnerSec[sec][row].GetN2();icl++)
2312 //if (cl[icl].IsUsed(10))
2317 for (Int_t sec=0;sec<fkNOS;sec++){
2318 for (Int_t row=0;row<fOuterSec->GetNRows();row++){
2319 AliTPCclusterMI *cl = fOuterSec[sec][row].GetClusters1();
2320 for (Int_t icl =0;icl< fOuterSec[sec][row].GetN1();icl++)
2321 //if (cl[icl].IsUsed(10))
2323 cl = fOuterSec[sec][row].GetClusters2();
2324 for (Int_t icl =0;icl< fOuterSec[sec][row].GetN2();icl++)
2325 //if (cl[icl].IsUsed(10))
2334 void AliTPCtrackerMI::SignClusters(TObjArray * arr, Float_t fnumber, Float_t fdensity)
2337 //sign clusters to be "used"
2339 // snumber and sdensity sign number of sigmas - bellow mean value to be accepted
2340 // loop over "primaries"
2354 Int_t nseed = arr->GetEntriesFast();
2355 for (Int_t i=0; i<nseed; i++) {
2356 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2360 if (!(pt->IsActive())) continue;
2361 Float_t dens = pt->GetNumberOfClusters()/Float_t(pt->GetNFoundable());
2362 if ( (dens>0.7) && (pt->GetNumberOfClusters()>70)){
2364 sumdens2+= dens*dens;
2365 sumn += pt->GetNumberOfClusters();
2366 sumn2 += pt->GetNumberOfClusters()*pt->GetNumberOfClusters();
2367 Float_t chi2 = pt->GetChi2()/pt->GetNumberOfClusters();
2370 sumchi2 +=chi2*chi2;
2375 Float_t mdensity = 0.9;
2376 Float_t meann = 130;
2377 Float_t meanchi = 1;
2378 Float_t sdensity = 0.1;
2379 Float_t smeann = 10;
2380 Float_t smeanchi =0.4;
2384 mdensity = sumdens/sum;
2386 meanchi = sumchi/sum;
2388 sdensity = sumdens2/sum-mdensity*mdensity;
2390 sdensity = TMath::Sqrt(sdensity);
2394 smeann = sumn2/sum-meann*meann;
2396 smeann = TMath::Sqrt(smeann);
2400 smeanchi = sumchi2/sum - meanchi*meanchi;
2402 smeanchi = TMath::Sqrt(smeanchi);
2408 //REMOVE SHORT DELTAS or tracks going out of sensitive volume of TPC
2410 for (Int_t i=0; i<nseed; i++) {
2411 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2415 if (pt->GetBSigned()) continue;
2416 if (pt->GetBConstrain()) continue;
2417 //if (!(pt->IsActive())) continue;
2419 Int_t found,foundable,shared;
2420 pt->GetClusterStatistic(0,160,found, foundable,shared);
2421 if (shared/float(found)>0.3) {
2422 if (shared/float(found)>0.9 ){
2423 //delete arr->RemoveAt(i);
2428 Bool_t isok =kFALSE;
2429 if ( (pt->GetNShared()/pt->GetNumberOfClusters()<0.5) &&pt->GetNumberOfClusters()>60)
2431 if ((TMath::Abs(1/pt->GetC())<100.) && (pt->GetNShared()/pt->GetNumberOfClusters()<0.7))
2433 if (TMath::Abs(pt->GetZ()/pt->GetX())>1.1)
2435 if ( (TMath::Abs(pt->GetSnp()>0.7) && pt->GetD(0,0)>60.))
2439 for (Int_t i=0; i<160; i++) {
2440 Int_t index=pt->GetClusterIndex2(i);
2441 if (index<0) continue;
2442 AliTPCclusterMI *c= pt->GetClusterPointer(i);
2444 //if (!(c->IsUsed(10))) c->Use();
2451 Double_t maxchi = meanchi+2.*smeanchi;
2453 for (Int_t i=0; i<nseed; i++) {
2454 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2458 //if (!(pt->IsActive())) continue;
2459 if (pt->GetBSigned()) continue;
2460 Double_t chi = pt->GetChi2()/pt->GetNumberOfClusters();
2461 if (chi>maxchi) continue;
2464 Float_t dens = pt->GetNumberOfClusters()/Float_t(pt->GetNFoundable());
2466 //sign only tracks with enoug big density at the beginning
2468 if ((pt->GetDensityFirst(40)<0.75) && pt->GetNumberOfClusters()<meann) continue;
2471 Double_t mindens = TMath::Max(double(mdensity-sdensity*fdensity*bfactor),0.65);
2472 Double_t minn = TMath::Max(Int_t(meann-fnumber*smeann*bfactor),50);
2474 // if (pt->fBConstrain) mindens = TMath::Max(mdensity-sdensity*fdensity*bfactor,0.65);
2475 if ( (pt->GetRemoval()==10) && (pt->GetSnp()>0.8)&&(dens>mindens))
2478 if ((dens>mindens && pt->GetNumberOfClusters()>minn) && chi<maxchi ){
2479 //Int_t noc=pt->GetNumberOfClusters();
2480 pt->SetBSigned(kTRUE);
2481 for (Int_t i=0; i<160; i++) {
2483 Int_t index=pt->GetClusterIndex2(i);
2484 if (index<0) continue;
2485 AliTPCclusterMI *c= pt->GetClusterPointer(i);
2487 // if (!(c->IsUsed(10))) c->Use();
2492 // gLastCheck = nseed;
2500 void AliTPCtrackerMI::StopNotActive(TObjArray * arr, Int_t row0, Float_t th0, Float_t th1, Float_t th2) const
2502 // stop not active tracks
2503 // take th1 as threshold for number of founded to number of foundable on last 10 active rows
2504 // take th2 as threshold for number of founded to number of foundable on last 20 active rows
2505 Int_t nseed = arr->GetEntriesFast();
2507 for (Int_t i=0; i<nseed; i++) {
2508 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2512 if (!(pt->IsActive())) continue;
2513 StopNotActive(pt,row0,th0, th1,th2);
2519 void AliTPCtrackerMI::StopNotActive(AliTPCseed * seed, Int_t row0, Float_t th0, Float_t th1,
2522 // stop not active tracks
2523 // take th1 as threshold for number of founded to number of foundable on last 10 active rows
2524 // take th2 as threshold for number of founded to number of foundable on last 20 active rows
2527 Int_t foundable = 0;
2528 Int_t maxindex = seed->GetLastPoint(); //last foundable row
2529 if (seed->GetNFoundable()*th0 > seed->GetNumberOfClusters()) {
2530 seed->Desactivate(10) ;
2534 for (Int_t i=row0; i<maxindex; i++){
2535 Int_t index = seed->GetClusterIndex2(i);
2536 if (index!=-1) foundable++;
2538 if (foundable<=30) sumgood1++;
2539 if (foundable<=50) {
2546 if (foundable>=30.){
2547 if (sumgood1<(th1*30.)) seed->Desactivate(10);
2550 if (sumgood2<(th2*50.)) seed->Desactivate(10);
2554 Int_t AliTPCtrackerMI::RefitInward(AliESDEvent *event)
2557 // back propagation of ESD tracks
2560 const Int_t kMaxFriendTracks=2000;
2564 //PrepareForProlongation(fSeeds,1);
2565 PropagateForward2(fSeeds);
2566 RemoveUsed2(fSeeds,0.4,0.4,20);
2568 TObjArray arraySeed(fSeeds->GetEntries());
2569 for (Int_t i=0;i<fSeeds->GetEntries();i++) {
2570 arraySeed.AddAt(fSeeds->At(i),i);
2572 SignShared(&arraySeed);
2573 // FindCurling(fSeeds, event,2); // find multi found tracks
2574 FindSplitted(fSeeds, event,2); // find multi found tracks
2577 Int_t nseed = fSeeds->GetEntriesFast();
2578 for (Int_t i=0;i<nseed;i++){
2579 AliTPCseed * seed = (AliTPCseed*) fSeeds->UncheckedAt(i);
2580 if (!seed) continue;
2581 if (seed->GetKinkIndex(0)>0) UpdateKinkQualityD(seed); // update quality informations for kinks
2583 seed->PropagateTo(fParam->GetInnerRadiusLow());
2584 seed->UpdatePoints();
2585 AddCovariance(seed);
2587 AliESDtrack *esd=event->GetTrack(i);
2588 seed->CookdEdx(0.02,0.6);
2589 CookLabel(seed,0.1); //For comparison only
2590 esd->SetTPCClusterMap(seed->GetClusterMap());
2591 esd->SetTPCSharedMap(seed->GetSharedMap());
2593 if (AliTPCReconstructor::StreamLevel()>1 && seed!=0&&esd!=0) {
2594 TTreeSRedirector &cstream = *fDebugStreamer;
2601 if (seed->GetNumberOfClusters()>15){
2602 esd->UpdateTrackParams(seed,AliESDtrack::kTPCrefit);
2603 esd->SetTPCPoints(seed->GetPoints());
2604 esd->SetTPCPointsF(seed->GetNFoundable());
2605 Int_t ndedx = seed->GetNCDEDX(0)+seed->GetNCDEDX(1)+seed->GetNCDEDX(2)+seed->GetNCDEDX(3);
2606 Float_t sdedx = (seed->GetSDEDX(0)+seed->GetSDEDX(1)+seed->GetSDEDX(2)+seed->GetSDEDX(3))*0.25;
2607 Float_t dedx = seed->GetdEdx();
2608 esd->SetTPCsignal(dedx, sdedx, ndedx);
2610 // add seed to the esd track in Calib level
2612 Bool_t storeFriend = gRandom->Rndm()<(kMaxFriendTracks)/Float_t(nseed);
2613 if (AliTPCReconstructor::StreamLevel()>0 &&storeFriend){
2614 AliTPCseed * seedCopy = new AliTPCseed(*seed, kTRUE);
2615 esd->AddCalibObject(seedCopy);
2620 //printf("problem\n");
2623 //FindKinks(fSeeds,event);
2624 Info("RefitInward","Number of refitted tracks %d",ntracks);
2629 Int_t AliTPCtrackerMI::PropagateBack(AliESDEvent *event)
2632 // back propagation of ESD tracks
2638 PropagateBack(fSeeds);
2639 RemoveUsed2(fSeeds,0.4,0.4,20);
2640 //FindCurling(fSeeds, fEvent,1);
2641 FindSplitted(fSeeds, event,1); // find multi found tracks
2644 Int_t nseed = fSeeds->GetEntriesFast();
2646 for (Int_t i=0;i<nseed;i++){
2647 AliTPCseed * seed = (AliTPCseed*) fSeeds->UncheckedAt(i);
2648 if (!seed) continue;
2649 if (seed->GetKinkIndex(0)<0) UpdateKinkQualityM(seed); // update quality informations for kinks
2650 seed->UpdatePoints();
2651 AddCovariance(seed);
2652 AliESDtrack *esd=event->GetTrack(i);
2653 seed->CookdEdx(0.02,0.6);
2654 CookLabel(seed,0.1); //For comparison only
2655 if (seed->GetNumberOfClusters()>15){
2656 esd->UpdateTrackParams(seed,AliESDtrack::kTPCout);
2657 esd->SetTPCPoints(seed->GetPoints());
2658 esd->SetTPCPointsF(seed->GetNFoundable());
2659 Int_t ndedx = seed->GetNCDEDX(0)+seed->GetNCDEDX(1)+seed->GetNCDEDX(2)+seed->GetNCDEDX(3);
2660 Float_t sdedx = (seed->GetSDEDX(0)+seed->GetSDEDX(1)+seed->GetSDEDX(2)+seed->GetSDEDX(3))*0.25;
2661 Float_t dedx = seed->GetdEdx();
2662 esd->SetTPCsignal(dedx, sdedx, ndedx);
2664 Int_t eventnumber = event->GetEventNumberInFile();// patch 28 fev 06
2665 // This is most likely NOT the event number you'd like to use. It has nothing to do with the 'real' event number
2666 if (AliTPCReconstructor::StreamLevel()>1 && esd) {
2667 (*fDebugStreamer)<<"Cback"<<
2670 "EventNrInFile="<<eventnumber<<
2675 //FindKinks(fSeeds,event);
2676 Info("PropagateBack","Number of back propagated tracks %d",ntracks);
2683 void AliTPCtrackerMI::DeleteSeeds()
2688 Int_t nseed = fSeeds->GetEntriesFast();
2689 for (Int_t i=0;i<nseed;i++){
2690 AliTPCseed * seed = (AliTPCseed*)fSeeds->At(i);
2691 if (seed) delete fSeeds->RemoveAt(i);
2698 void AliTPCtrackerMI::ReadSeeds(AliESDEvent *event, Int_t direction)
2701 //read seeds from the event
2703 Int_t nentr=event->GetNumberOfTracks();
2705 Info("ReadSeeds", "Number of ESD tracks: %d\n", nentr);
2710 fSeeds = new TObjArray(nentr);
2714 for (Int_t i=0; i<nentr; i++) {
2715 AliESDtrack *esd=event->GetTrack(i);
2716 ULong_t status=esd->GetStatus();
2717 if (!(status&AliESDtrack::kTPCin)) continue;
2718 AliTPCtrack t(*esd);
2719 t.SetNumberOfClusters(0);
2720 // AliTPCseed *seed = new AliTPCseed(t,t.GetAlpha());
2721 AliTPCseed *seed = new AliTPCseed(t/*,t.GetAlpha()*/);
2722 seed->SetUniqueID(esd->GetID());
2723 AddCovariance(seed); //add systematic ucertainty
2724 for (Int_t ikink=0;ikink<3;ikink++) {
2725 Int_t index = esd->GetKinkIndex(ikink);
2726 seed->GetKinkIndexes()[ikink] = index;
2727 if (index==0) continue;
2728 index = TMath::Abs(index);
2729 AliESDkink * kink = fEvent->GetKink(index-1);
2730 if (kink&&esd->GetKinkIndex(ikink)<0){
2731 if ((status & AliESDtrack::kTRDrefit) != 0) kink->SetStatus(1,2);
2732 if ((status & AliESDtrack::kITSout) != 0) kink->SetStatus(1,0);
2734 if (kink&&esd->GetKinkIndex(ikink)>0){
2735 if ((status & AliESDtrack::kTRDrefit) != 0) kink->SetStatus(1,6);
2736 if ((status & AliESDtrack::kITSout) != 0) kink->SetStatus(1,4);
2740 if (((status&AliESDtrack::kITSout)==0)&&(direction==1)) seed->ResetCovariance(10.);
2741 if ( direction ==2 &&(status & AliESDtrack::kTRDrefit) == 0 ) seed->ResetCovariance(10.);
2742 if ( direction ==2 && ((status & AliESDtrack::kTPCout) == 0) ) {
2747 if ( direction ==2 &&(status & AliESDtrack::kTRDrefit) > 0 ) {
2748 Double_t par0[5],par1[5],alpha,x;
2749 esd->GetInnerExternalParameters(alpha,x,par0);
2750 esd->GetExternalParameters(x,par1);
2751 Double_t delta1 = TMath::Abs(par0[4]-par1[4])/(0.000000001+TMath::Abs(par0[4]+par1[4]));
2752 Double_t delta2 = TMath::Abs(par0[3]-par1[3]);
2754 if (esd->GetTRDncls()>0) trdchi2 = esd->GetTRDchi2()/esd->GetTRDncls();
2755 //reset covariance if suspicious
2756 if ( (delta1>0.1) || (delta2>0.006) ||trdchi2>7.)
2757 seed->ResetCovariance(10.);
2762 // rotate to the local coordinate system
2764 fSectors=fInnerSec; fN=fkNIS;
2765 Double_t alpha=seed->GetAlpha() - fSectors->GetAlphaShift();
2766 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
2767 if (alpha < 0. ) alpha += 2.*TMath::Pi();
2768 Int_t ns=Int_t(alpha/fSectors->GetAlpha())%fN;
2769 alpha =ns*fSectors->GetAlpha() + fSectors->GetAlphaShift();
2770 if (alpha<-TMath::Pi()) alpha += 2*TMath::Pi();
2771 if (alpha>=TMath::Pi()) alpha -= 2*TMath::Pi();
2772 alpha-=seed->GetAlpha();
2773 if (!seed->Rotate(alpha)) {
2779 if (esd->GetKinkIndex(0)<=0){
2780 for (Int_t irow=0;irow<160;irow++){
2781 Int_t index = seed->GetClusterIndex2(irow);
2784 AliTPCclusterMI * cl = GetClusterMI(index);
2785 seed->SetClusterPointer(irow,cl);
2787 if ((index & 0x8000)==0){
2788 cl->Use(10); // accepted cluster
2790 cl->Use(6); // close cluster not accepted
2793 Info("ReadSeeds","Not found cluster");
2798 fSeeds->AddAt(seed,i);
2804 //_____________________________________________________________________________
2805 void AliTPCtrackerMI::MakeSeeds3(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, Float_t cuts[4],
2806 Float_t deltay, Int_t ddsec) {
2807 //-----------------------------------------------------------------
2808 // This function creates track seeds.
2809 // SEEDING WITH VERTEX CONSTRAIN
2810 //-----------------------------------------------------------------
2811 // cuts[0] - fP4 cut
2812 // cuts[1] - tan(phi) cut
2813 // cuts[2] - zvertex cut
2814 // cuts[3] - fP3 cut
2822 Double_t x[5], c[15];
2823 // Int_t di = i1-i2;
2825 AliTPCseed * seed = new AliTPCseed();
2826 Double_t alpha=fSectors->GetAlpha(), shift=fSectors->GetAlphaShift();
2827 Double_t cs=cos(alpha), sn=sin(alpha);
2829 // Double_t x1 =fOuterSec->GetX(i1);
2830 //Double_t xx2=fOuterSec->GetX(i2);
2832 Double_t x1 =GetXrow(i1);
2833 Double_t xx2=GetXrow(i2);
2835 Double_t x3=GetX(), y3=GetY(), z3=GetZ();
2837 Int_t imiddle = (i2+i1)/2; //middle pad row index
2838 Double_t xm = GetXrow(imiddle); // radius of middle pad-row
2839 const AliTPCtrackerRow& krm=GetRow(sec,imiddle); //middle pad -row
2843 const AliTPCtrackerRow& kr1=GetRow(ns,i1);
2844 Double_t ymax = GetMaxY(i1)-kr1.GetDeadZone()-1.5;
2845 Double_t ymaxm = GetMaxY(imiddle)-kr1.GetDeadZone()-1.5;
2848 // change cut on curvature if it can't reach this layer
2849 // maximal curvature set to reach it
2850 Double_t dvertexmax = TMath::Sqrt((x1-x3)*(x1-x3)+(ymax+5-y3)*(ymax+5-y3));
2851 if (dvertexmax*0.5*cuts[0]>0.85){
2852 cuts[0] = 0.85/(dvertexmax*0.5+1.);
2854 Double_t r2min = 1/(cuts[0]*cuts[0]); //minimal square of radius given by cut
2857 if (deltay>0) ddsec = 0;
2858 // loop over clusters
2859 for (Int_t is=0; is < kr1; is++) {
2861 if (kr1[is]->IsUsed(10)) continue;
2862 Double_t y1=kr1[is]->GetY(), z1=kr1[is]->GetZ();
2863 //if (TMath::Abs(y1)>ymax) continue;
2865 if (deltay>0 && TMath::Abs(ymax-TMath::Abs(y1))> deltay ) continue; // seed only at the edge
2867 // find possible directions
2868 Float_t anglez = (z1-z3)/(x1-x3);
2869 Float_t extraz = z1 - anglez*(x1-xx2); // extrapolated z
2872 //find rotation angles relative to line given by vertex and point 1
2873 Double_t dvertex2 = (x1-x3)*(x1-x3)+(y1-y3)*(y1-y3);
2874 Double_t dvertex = TMath::Sqrt(dvertex2);
2875 Double_t angle13 = TMath::ATan((y1-y3)/(x1-x3));
2876 Double_t cs13 = cos(-angle13), sn13 = sin(-angle13);
2879 // loop over 2 sectors
2885 Double_t dddz1=0; // direction of delta inclination in z axis
2892 for (Int_t dsec = dsec1; dsec<=dsec2;dsec++){
2893 Int_t sec2 = sec + dsec;
2895 // AliTPCtrackerRow& kr2 = fOuterSec[(sec2+fkNOS)%fkNOS][i2];
2896 //AliTPCtrackerRow& kr2m = fOuterSec[(sec2+fkNOS)%fkNOS][imiddle];
2897 AliTPCtrackerRow& kr2 = GetRow((sec2+fkNOS)%fkNOS,i2);
2898 AliTPCtrackerRow& kr2m = GetRow((sec2+fkNOS)%fkNOS,imiddle);
2899 Int_t index1 = TMath::Max(kr2.Find(extraz-0.6-dddz1*TMath::Abs(z1)*0.05)-1,0);
2900 Int_t index2 = TMath::Min(kr2.Find(extraz+0.6+dddz2*TMath::Abs(z1)*0.05)+1,kr2);
2902 // rotation angles to p1-p3
2903 Double_t cs13r = cos(-angle13+dsec*alpha)/dvertex, sn13r = sin(-angle13+dsec*alpha)/dvertex;
2904 Double_t x2, y2, z2;
2906 // Double_t dymax = maxangle*TMath::Abs(x1-xx2);
2909 Double_t dxx0 = (xx2-x3)*cs13r;
2910 Double_t dyy0 = (xx2-x3)*sn13r;
2911 for (Int_t js=index1; js < index2; js++) {
2912 const AliTPCclusterMI *kcl = kr2[js];
2913 if (kcl->IsUsed(10)) continue;
2915 //calcutate parameters
2917 Double_t yy0 = dyy0 +(kcl->GetY()-y3)*cs13r;
2919 if (TMath::Abs(yy0)<0.000001) continue;
2920 Double_t xx0 = dxx0 -(kcl->GetY()-y3)*sn13r;
2921 Double_t y0 = 0.5*(xx0*xx0+yy0*yy0-xx0)/yy0;
2922 Double_t r02 = (0.25+y0*y0)*dvertex2;
2923 //curvature (radius) cut
2924 if (r02<r2min) continue;
2928 Double_t c0 = 1/TMath::Sqrt(r02);
2932 //Double_t dfi0 = 2.*TMath::ASin(dvertex*c0*0.5);
2933 //Double_t dfi1 = 2.*TMath::ASin(TMath::Sqrt(yy0*yy0+(1-xx0)*(1-xx0))*dvertex*c0*0.5);
2934 Double_t dfi0 = 2.*AliTPCFastMath::FastAsin(dvertex*c0*0.5);
2935 Double_t dfi1 = 2.*AliTPCFastMath::FastAsin(TMath::Sqrt(yy0*yy0+(1-xx0)*(1-xx0))*dvertex*c0*0.5);
2938 Double_t z0 = kcl->GetZ();
2939 Double_t zzzz2 = z1-(z1-z3)*dfi1/dfi0;
2940 if (TMath::Abs(zzzz2-z0)>0.5) continue;
2943 Double_t dip = (z1-z0)*c0/dfi1;
2944 Double_t x0 = (0.5*cs13+y0*sn13)*dvertex*c0;
2955 x2= xx2*cs-y2*sn*dsec;
2956 y2=+xx2*sn*dsec+y2*cs;
2966 // do we have cluster at the middle ?
2968 GetProlongation(x1,xm,x,ym,zm);
2970 AliTPCclusterMI * cm=0;
2971 if (TMath::Abs(ym)-ymaxm<0){
2972 cm = krm.FindNearest2(ym,zm,1.0,0.6,dummy);
2973 if ((!cm) || (cm->IsUsed(10))) {
2978 // rotate y1 to system 0
2979 // get state vector in rotated system
2980 Double_t yr1 = (-0.5*sn13+y0*cs13)*dvertex*c0;
2981 Double_t xr2 = x0*cs+yr1*sn*dsec;
2982 Double_t xr[5]={kcl->GetY(),kcl->GetZ(), xr2, dip, c0};
2984 GetProlongation(xx2,xm,xr,ym,zm);
2985 if (TMath::Abs(ym)-ymaxm<0){
2986 cm = kr2m.FindNearest2(ym,zm,1.0,0.6,dummy);
2987 if ((!cm) || (cm->IsUsed(10))) {
2997 dym = ym - cm->GetY();
2998 dzm = zm - cm->GetZ();
3005 Double_t sy1=kr1[is]->GetSigmaY2()*2., sz1=kr1[is]->GetSigmaZ2()*2.;
3006 Double_t sy2=kcl->GetSigmaY2()*2., sz2=kcl->GetSigmaZ2()*2.;
3007 //Double_t sy3=400*3./12., sy=0.1, sz=0.1;
3008 Double_t sy3=25000*x[4]*x[4]+0.1, sy=0.1, sz=0.1;
3009 //Double_t sy3=25000*x[4]*x[4]*60+0.5, sy=0.1, sz=0.1;
3011 Double_t f40=(F1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
3012 Double_t f42=(F1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
3013 Double_t f43=(F1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
3014 Double_t f20=(F2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
3015 Double_t f22=(F2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
3016 Double_t f23=(F2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
3018 Double_t f30=(F3(x1,y1+sy,x2,y2,z1,z2)-x[3])/sy;
3019 Double_t f31=(F3(x1,y1,x2,y2,z1+sz,z2)-x[3])/sz;
3020 Double_t f32=(F3(x1,y1,x2,y2+sy,z1,z2)-x[3])/sy;
3021 Double_t f34=(F3(x1,y1,x2,y2,z1,z2+sz)-x[3])/sz;
3025 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3026 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3027 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3028 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3029 c[13]=f30*sy1*f40+f32*sy2*f42;
3030 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3032 // if (!BuildSeed(kr1[is],kcl,0,x1,x2,x3,x,c)) continue;
3034 UInt_t index=kr1.GetIndex(is);
3035 seed->~AliTPCseed(); // this does not set the pointer to 0...
3036 AliTPCseed *track=new(seed) AliTPCseed(x1, ns*alpha+shift, x, c, index);
3038 track->SetIsSeeding(kTRUE);
3039 track->SetSeed1(i1);
3040 track->SetSeed2(i2);
3041 track->SetSeedType(3);
3045 FollowProlongation(*track, (i1+i2)/2,1);
3046 Int_t foundable,found,shared;
3047 track->GetClusterStatistic((i1+i2)/2,i1, found, foundable, shared, kTRUE);
3048 if ((found<0.55*foundable) || shared>0.5*found || (track->GetSigmaY2()+track->GetSigmaZ2())>0.5){
3050 seed->~AliTPCseed();
3056 FollowProlongation(*track, i2,1);
3060 track->SetBConstrain(1);
3061 // track->fLastPoint = i1+fInnerSec->GetNRows(); // first cluster in track position
3062 track->SetLastPoint(i1); // first cluster in track position
3063 track->SetFirstPoint(track->GetLastPoint());
3065 if (track->GetNumberOfClusters()<(i1-i2)*0.5 ||
3066 track->GetNumberOfClusters() < track->GetNFoundable()*0.6 ||
3067 track->GetNShared()>0.4*track->GetNumberOfClusters() ) {
3069 seed->~AliTPCseed();
3073 // Z VERTEX CONDITION
3074 Double_t zv, bz=GetBz();
3075 if ( !track->GetZAt(0.,bz,zv) ) continue;
3076 if (TMath::Abs(zv-z3)>cuts[2]) {
3077 FollowProlongation(*track, TMath::Max(i2-20,0));
3078 if ( !track->GetZAt(0.,bz,zv) ) continue;
3079 if (TMath::Abs(zv-z3)>cuts[2]){
3080 FollowProlongation(*track, TMath::Max(i2-40,0));
3081 if ( !track->GetZAt(0.,bz,zv) ) continue;
3082 if (TMath::Abs(zv-z3)>cuts[2] &&(track->GetNumberOfClusters() > track->GetNFoundable()*0.7)){
3083 // make seed without constrain
3084 AliTPCseed * track2 = MakeSeed(track,0.2,0.5,1.);
3085 FollowProlongation(*track2, i2,1);
3086 track2->SetBConstrain(kFALSE);
3087 track2->SetSeedType(1);
3088 arr->AddLast(track2);
3090 seed->~AliTPCseed();
3095 seed->~AliTPCseed();
3102 track->SetSeedType(0);
3103 arr->AddLast(track);
3104 seed = new AliTPCseed;
3106 // don't consider other combinations
3107 if (track->GetNumberOfClusters() > track->GetNFoundable()*0.8)
3113 Info("MakeSeeds3","\nSeeding statistic:\t%d\t%d\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin,nout1,nout2);
3119 void AliTPCtrackerMI::MakeSeeds5(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, Float_t cuts[4],
3124 //-----------------------------------------------------------------
3125 // This function creates track seeds.
3126 //-----------------------------------------------------------------
3127 // cuts[0] - fP4 cut
3128 // cuts[1] - tan(phi) cut
3129 // cuts[2] - zvertex cut
3130 // cuts[3] - fP3 cut
3140 Double_t x[5], c[15];
3142 // make temporary seed
3143 AliTPCseed * seed = new AliTPCseed;
3144 Double_t alpha=fOuterSec->GetAlpha(), shift=fOuterSec->GetAlphaShift();
3145 // Double_t cs=cos(alpha), sn=sin(alpha);
3150 Double_t x1 = GetXrow(i1-1);
3151 const AliTPCtrackerRow& kr1=GetRow(sec,i1-1);
3152 Double_t y1max = GetMaxY(i1-1)-kr1.GetDeadZone()-1.5;
3154 Double_t x1p = GetXrow(i1);
3155 const AliTPCtrackerRow& kr1p=GetRow(sec,i1);
3157 Double_t x1m = GetXrow(i1-2);
3158 const AliTPCtrackerRow& kr1m=GetRow(sec,i1-2);
3161 //last 3 padrow for seeding
3162 AliTPCtrackerRow& kr3 = GetRow((sec+fkNOS)%fkNOS,i1-7);
3163 Double_t x3 = GetXrow(i1-7);
3164 // Double_t y3max= GetMaxY(i1-7)-kr3.fDeadZone-1.5;
3166 AliTPCtrackerRow& kr3p = GetRow((sec+fkNOS)%fkNOS,i1-6);
3167 Double_t x3p = GetXrow(i1-6);
3169 AliTPCtrackerRow& kr3m = GetRow((sec+fkNOS)%fkNOS,i1-8);
3170 Double_t x3m = GetXrow(i1-8);
3175 Int_t im = i1-4; //middle pad row index
3176 Double_t xm = GetXrow(im); // radius of middle pad-row
3177 const AliTPCtrackerRow& krm=GetRow(sec,im); //middle pad -row
3178 // Double_t ymmax = GetMaxY(im)-kr1.fDeadZone-1.5;
3181 Double_t deltax = x1-x3;
3182 Double_t dymax = deltax*cuts[1];
3183 Double_t dzmax = deltax*cuts[3];
3185 // loop over clusters
3186 for (Int_t is=0; is < kr1; is++) {
3188 if (kr1[is]->IsUsed(10)) continue;
3189 Double_t y1=kr1[is]->GetY(), z1=kr1[is]->GetZ();
3191 if (deltay>0 && TMath::Abs(y1max-TMath::Abs(y1))> deltay ) continue; // seed only at the edge
3193 Int_t index1 = TMath::Max(kr3.Find(z1-dzmax)-1,0);
3194 Int_t index2 = TMath::Min(kr3.Find(z1+dzmax)+1,kr3);
3200 for (Int_t js=index1; js < index2; js++) {
3201 const AliTPCclusterMI *kcl = kr3[js];
3202 if (kcl->IsUsed(10)) continue;
3204 // apply angular cuts
3205 if (TMath::Abs(y1-y3)>dymax) continue;
3208 if (TMath::Abs(z1-z3)>dzmax) continue;
3210 Double_t angley = (y1-y3)/(x1-x3);
3211 Double_t anglez = (z1-z3)/(x1-x3);
3213 Double_t erry = TMath::Abs(angley)*(x1-x1m)*0.5+0.5;
3214 Double_t errz = TMath::Abs(anglez)*(x1-x1m)*0.5+0.5;
3216 Double_t yyym = angley*(xm-x1)+y1;
3217 Double_t zzzm = anglez*(xm-x1)+z1;
3219 const AliTPCclusterMI *kcm = krm.FindNearest2(yyym,zzzm,erry,errz,index);
3221 if (kcm->IsUsed(10)) continue;
3223 erry = TMath::Abs(angley)*(x1-x1m)*0.4+0.5;
3224 errz = TMath::Abs(anglez)*(x1-x1m)*0.4+0.5;
3231 // look around first
3232 const AliTPCclusterMI *kc1m = kr1m.FindNearest2(angley*(x1m-x1)+y1,
3238 if (kc1m->IsUsed(10)) used++;
3240 const AliTPCclusterMI *kc1p = kr1p.FindNearest2(angley*(x1p-x1)+y1,
3246 if (kc1p->IsUsed(10)) used++;
3248 if (used>1) continue;
3249 if (found<1) continue;
3253 const AliTPCclusterMI *kc3m = kr3m.FindNearest2(angley*(x3m-x3)+y3,
3259 if (kc3m->IsUsed(10)) used++;
3263 const AliTPCclusterMI *kc3p = kr3p.FindNearest2(angley*(x3p-x3)+y3,
3269 if (kc3p->IsUsed(10)) used++;
3273 if (used>1) continue;
3274 if (found<3) continue;
3284 x[4]=F1(x1,y1,x2,y2,x3,y3);
3285 //if (TMath::Abs(x[4]) >= cuts[0]) continue;
3288 x[2]=F2(x1,y1,x2,y2,x3,y3);
3291 x[3]=F3n(x1,y1,x2,y2,z1,z2,x[4]);
3292 //if (TMath::Abs(x[3]) > cuts[3]) continue;
3296 Double_t sy1=0.1, sz1=0.1;
3297 Double_t sy2=0.1, sz2=0.1;
3298 Double_t sy3=0.1, sy=0.1, sz=0.1;
3300 Double_t f40=(F1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
3301 Double_t f42=(F1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
3302 Double_t f43=(F1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
3303 Double_t f20=(F2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
3304 Double_t f22=(F2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
3305 Double_t f23=(F2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
3307 Double_t f30=(F3(x1,y1+sy,x2,y2,z1,z2)-x[3])/sy;
3308 Double_t f31=(F3(x1,y1,x2,y2,z1+sz,z2)-x[3])/sz;
3309 Double_t f32=(F3(x1,y1,x2,y2+sy,z1,z2)-x[3])/sy;
3310 Double_t f34=(F3(x1,y1,x2,y2,z1,z2+sz)-x[3])/sz;
3314 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3315 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3316 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3317 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3318 c[13]=f30*sy1*f40+f32*sy2*f42;
3319 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3321 // if (!BuildSeed(kr1[is],kcl,0,x1,x2,x3,x,c)) continue;
3323 UInt_t index=kr1.GetIndex(is);
3324 seed->~AliTPCseed();
3325 AliTPCseed *track=new(seed) AliTPCseed(x1, sec*alpha+shift, x, c, index);
3327 track->SetIsSeeding(kTRUE);
3330 FollowProlongation(*track, i1-7,1);
3331 if (track->GetNumberOfClusters() < track->GetNFoundable()*0.75 ||
3332 track->GetNShared()>0.6*track->GetNumberOfClusters() || ( track->GetSigmaY2()+ track->GetSigmaZ2())>0.6){
3334 seed->~AliTPCseed();
3340 FollowProlongation(*track, i2,1);
3341 track->SetBConstrain(0);
3342 track->SetLastPoint(i1+fInnerSec->GetNRows()); // first cluster in track position
3343 track->SetFirstPoint(track->GetLastPoint());
3345 if (track->GetNumberOfClusters()<(i1-i2)*0.5 ||
3346 track->GetNumberOfClusters()<track->GetNFoundable()*0.7 ||
3347 track->GetNShared()>2. || track->GetChi2()/track->GetNumberOfClusters()>6 || ( track->GetSigmaY2()+ track->GetSigmaZ2())>0.5 ) {
3349 seed->~AliTPCseed();
3354 FollowProlongation(*track, TMath::Max(i2-10,0),1);
3355 AliTPCseed * track2 = MakeSeed(track,0.2,0.5,0.9);
3356 FollowProlongation(*track2, i2,1);
3357 track2->SetBConstrain(kFALSE);
3358 track2->SetSeedType(4);
3359 arr->AddLast(track2);
3361 seed->~AliTPCseed();
3365 //arr->AddLast(track);
3366 //seed = new AliTPCseed;
3372 Info("MakeSeeds5","\nSeeding statiistic:\t%d\t%d\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin,nout1,nout2,nout3);
3378 //_____________________________________________________________________________
3379 void AliTPCtrackerMI::MakeSeeds2(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, Float_t */*cuts[4]*/,
3380 Float_t deltay, Bool_t /*bconstrain*/) {
3381 //-----------------------------------------------------------------
3382 // This function creates track seeds - without vertex constraint
3383 //-----------------------------------------------------------------
3384 // cuts[0] - fP4 cut - not applied
3385 // cuts[1] - tan(phi) cut
3386 // cuts[2] - zvertex cut - not applied
3387 // cuts[3] - fP3 cut
3397 Double_t alpha=fOuterSec->GetAlpha(), shift=fOuterSec->GetAlphaShift();
3398 // Double_t cs=cos(alpha), sn=sin(alpha);
3399 Int_t row0 = (i1+i2)/2;
3400 Int_t drow = (i1-i2)/2;
3401 const AliTPCtrackerRow& kr0=fSectors[sec][row0];
3402 AliTPCtrackerRow * kr=0;
3404 AliTPCpolyTrack polytrack;
3405 Int_t nclusters=fSectors[sec][row0];
3406 AliTPCseed * seed = new AliTPCseed;
3411 for (Int_t is=0; is < nclusters; is++) { //LOOP over clusters
3413 Int_t nfoundable =0;
3414 for (Int_t iter =1; iter<2; iter++){ //iterations
3415 const AliTPCtrackerRow& krm=fSectors[sec][row0-iter];
3416 const AliTPCtrackerRow& krp=fSectors[sec][row0+iter];
3417 const AliTPCclusterMI * cl= kr0[is];
3419 if (cl->IsUsed(10)) {
3425 Double_t x = kr0.GetX();
3426 // Initialization of the polytrack
3431 Double_t y0= cl->GetY();
3432 Double_t z0= cl->GetZ();
3436 Double_t ymax = fSectors->GetMaxY(row0)-kr0.GetDeadZone()-1.5;
3437 if (deltay>0 && TMath::Abs(ymax-TMath::Abs(y0))> deltay ) continue; // seed only at the edge
3439 erry = (0.5)*cl->GetSigmaY2()/TMath::Sqrt(cl->GetQ())*6;
3440 errz = (0.5)*cl->GetSigmaZ2()/TMath::Sqrt(cl->GetQ())*6;
3441 polytrack.AddPoint(x,y0,z0,erry, errz);
3444 if (cl->IsUsed(10)) sumused++;
3447 Float_t roady = (5*TMath::Sqrt(cl->GetSigmaY2()+0.2)+1.)*iter;
3448 Float_t roadz = (5*TMath::Sqrt(cl->GetSigmaZ2()+0.2)+1.)*iter;
3451 AliTPCclusterMI * cl1 = krm.FindNearest(y0,z0,roady,roadz);
3452 if (cl1 && TMath::Abs(ymax-TMath::Abs(y0))) {
3453 erry = (0.5)*cl1->GetSigmaY2()/TMath::Sqrt(cl1->GetQ())*3;
3454 errz = (0.5)*cl1->GetSigmaZ2()/TMath::Sqrt(cl1->GetQ())*3;
3455 if (cl1->IsUsed(10)) sumused++;
3456 polytrack.AddPoint(x,cl1->GetY(),cl1->GetZ(),erry,errz);
3460 AliTPCclusterMI * cl2 = krp.FindNearest(y0,z0,roady,roadz);
3462 erry = (0.5)*cl2->GetSigmaY2()/TMath::Sqrt(cl2->GetQ())*3;
3463 errz = (0.5)*cl2->GetSigmaZ2()/TMath::Sqrt(cl2->GetQ())*3;
3464 if (cl2->IsUsed(10)) sumused++;
3465 polytrack.AddPoint(x,cl2->GetY(),cl2->GetZ(),erry,errz);
3468 if (sumused>0) continue;
3470 polytrack.UpdateParameters();
3476 nfoundable = polytrack.GetN();
3477 nfound = nfoundable;
3479 for (Int_t ddrow = iter+1; ddrow<drow;ddrow++){
3480 Float_t maxdist = 0.8*(1.+3./(ddrow));
3481 for (Int_t delta = -1;delta<=1;delta+=2){
3482 Int_t row = row0+ddrow*delta;
3483 kr = &(fSectors[sec][row]);
3484 Double_t xn = kr->GetX();
3485 Double_t ymax = fSectors->GetMaxY(row)-kr->GetDeadZone()-1.5;
3486 polytrack.GetFitPoint(xn,yn,zn);
3487 if (TMath::Abs(yn)>ymax) continue;
3489 AliTPCclusterMI * cln = kr->FindNearest(yn,zn,roady,roadz);
3491 Float_t dist = TMath::Sqrt( (yn-cln->GetY())*(yn-cln->GetY())+(zn-cln->GetZ())*(zn-cln->GetZ()));
3494 erry = (dist+0.3)*cln->GetSigmaY2()/TMath::Sqrt(cln->GetQ())*(1.+1./(ddrow));
3495 errz = (dist+0.3)*cln->GetSigmaZ2()/TMath::Sqrt(cln->GetQ())*(1.+1./(ddrow));
3496 if (cln->IsUsed(10)) {
3497 // printf("used\n");
3505 polytrack.AddPoint(xn,cln->GetY(),cln->GetZ(),erry, errz);
3510 if ( (sumused>3) || (sumused>0.5*nfound) || (nfound<0.6*nfoundable)) break;
3511 polytrack.UpdateParameters();
3514 if ( (sumused>3) || (sumused>0.5*nfound)) {
3515 //printf("sumused %d\n",sumused);
3520 polytrack.GetFitDerivation(kr0.GetX(),dy,dz);
3521 AliTPCpolyTrack track2;
3523 polytrack.Refit(track2,0.5+TMath::Abs(dy)*0.3,0.4+TMath::Abs(dz)*0.3);
3524 if (track2.GetN()<0.5*nfoundable) continue;
3527 if ((nfound>0.6*nfoundable) &&( nfoundable>0.4*(i1-i2))) {
3529 // test seed with and without constrain
3530 for (Int_t constrain=0; constrain<=0;constrain++){
3531 // add polytrack candidate
3533 Double_t x[5], c[15];
3534 Double_t x1,x2,x3,y1,y2,y3,z1,z2,z3;
3535 track2.GetBoundaries(x3,x1);
3537 track2.GetFitPoint(x1,y1,z1);
3538 track2.GetFitPoint(x2,y2,z2);
3539 track2.GetFitPoint(x3,y3,z3);
3541 //is track pointing to the vertex ?
3544 polytrack.GetFitPoint(x0,y0,z0);
3557 x[4]=F1(x1,y1,x2,y2,x3,y3);
3559 // if (TMath::Abs(x[4]) >= cuts[0]) continue; //
3560 x[2]=F2(x1,y1,x2,y2,x3,y3);
3562 //if (TMath::Abs(x[4]*x1-x[2]) >= cuts[1]) continue;
3563 //x[3]=F3(x1,y1,x2,y2,z1,z2);
3564 x[3]=F3n(x1,y1,x3,y3,z1,z3,x[4]);
3565 //if (TMath::Abs(x[3]) > cuts[3]) continue;
3568 Double_t sy =0.1, sz =0.1;
3569 Double_t sy1=0.02, sz1=0.02;
3570 Double_t sy2=0.02, sz2=0.02;
3574 sy3=25000*x[4]*x[4]+0.1, sy=0.1, sz=0.1;
3577 Double_t f40=(F1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
3578 Double_t f42=(F1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
3579 Double_t f43=(F1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
3580 Double_t f20=(F2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
3581 Double_t f22=(F2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
3582 Double_t f23=(F2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
3584 Double_t f30=(F3(x1,y1+sy,x3,y3,z1,z3)-x[3])/sy;
3585 Double_t f31=(F3(x1,y1,x3,y3,z1+sz,z3)-x[3])/sz;
3586 Double_t f32=(F3(x1,y1,x3,y3+sy,z1,z3)-x[3])/sy;
3587 Double_t f34=(F3(x1,y1,x3,y3,z1,z3+sz)-x[3])/sz;
3592 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3593 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3594 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3595 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3596 c[13]=f30*sy1*f40+f32*sy2*f42;
3597 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3599 //Int_t row1 = fSectors->GetRowNumber(x1);
3600 Int_t row1 = GetRowNumber(x1);
3604 seed->~AliTPCseed();
3605 AliTPCseed *track=new(seed) AliTPCseed(x1,sec*alpha+shift,x,c,index);
3606 track->SetIsSeeding(kTRUE);
3607 Int_t rc=FollowProlongation(*track, i2);
3608 if (constrain) track->SetBConstrain(1);
3610 track->SetBConstrain(0);
3611 track->SetLastPoint(row1+fInnerSec->GetNRows()); // first cluster in track position
3612 track->SetFirstPoint(track->GetLastPoint());
3614 if (rc==0 || track->GetNumberOfClusters()<(i1-i2)*0.5 ||
3615 track->GetNumberOfClusters() < track->GetNFoundable()*0.6 ||
3616 track->GetNShared()>0.4*track->GetNumberOfClusters()) {
3619 seed->~AliTPCseed();
3622 arr->AddLast(track);
3623 seed = new AliTPCseed;
3627 } // if accepted seed
3630 Info("MakeSeeds2","\nSeeding statiistic:\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin3);
3636 AliTPCseed *AliTPCtrackerMI::MakeSeed(AliTPCseed *track, Float_t r0, Float_t r1, Float_t r2)
3640 //reseed using track points
3641 Int_t p0 = int(r0*track->GetNumberOfClusters()); // point 0
3642 Int_t p1 = int(r1*track->GetNumberOfClusters());
3643 Int_t p2 = int(r2*track->GetNumberOfClusters()); // last point
3645 Double_t x0[3],x1[3],x2[3];
3646 for (Int_t i=0;i<3;i++){
3652 // find track position at given ratio of the length
3653 Int_t sec0=0, sec1=0, sec2=0;
3656 for (Int_t i=0;i<160;i++){
3657 if (track->GetClusterPointer(i)){
3659 AliTPCTrackerPoint *trpoint =track->GetTrackPoint(i);
3660 if ( (index<p0) || x0[0]<0 ){
3661 if (trpoint->GetX()>1){
3662 clindex = track->GetClusterIndex2(i);
3664 x0[0] = trpoint->GetX();
3665 x0[1] = trpoint->GetY();
3666 x0[2] = trpoint->GetZ();
3667 sec0 = ((clindex&0xff000000)>>24)%18;
3672 if ( (index<p1) &&(trpoint->GetX()>1)){
3673 clindex = track->GetClusterIndex2(i);
3675 x1[0] = trpoint->GetX();
3676 x1[1] = trpoint->GetY();
3677 x1[2] = trpoint->GetZ();
3678 sec1 = ((clindex&0xff000000)>>24)%18;
3681 if ( (index<p2) &&(trpoint->GetX()>1)){
3682 clindex = track->GetClusterIndex2(i);
3684 x2[0] = trpoint->GetX();
3685 x2[1] = trpoint->GetY();
3686 x2[2] = trpoint->GetZ();
3687 sec2 = ((clindex&0xff000000)>>24)%18;
3694 Double_t alpha, cs,sn, xx2,yy2;
3696 alpha = (sec1-sec2)*fSectors->GetAlpha();
3697 cs = TMath::Cos(alpha);
3698 sn = TMath::Sin(alpha);
3699 xx2= x1[0]*cs-x1[1]*sn;
3700 yy2= x1[0]*sn+x1[1]*cs;
3704 alpha = (sec0-sec2)*fSectors->GetAlpha();
3705 cs = TMath::Cos(alpha);
3706 sn = TMath::Sin(alpha);
3707 xx2= x0[0]*cs-x0[1]*sn;
3708 yy2= x0[0]*sn+x0[1]*cs;
3714 Double_t x[5],c[15];
3718 x[4]=F1(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]);
3719 // if (x[4]>1) return 0;
3720 x[2]=F2(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]);
3721 x[3]=F3n(x2[0],x2[1],x0[0],x0[1],x2[2],x0[2],x[4]);
3722 //if (TMath::Abs(x[3]) > 2.2) return 0;
3723 //if (TMath::Abs(x[2]) > 1.99) return 0;
3725 Double_t sy =0.1, sz =0.1;
3727 Double_t sy1=0.02+track->GetSigmaY2(), sz1=0.02+track->GetSigmaZ2();
3728 Double_t sy2=0.01+track->GetSigmaY2(), sz2=0.01+track->GetSigmaZ2();
3729 Double_t sy3=0.01+track->GetSigmaY2();
3731 Double_t f40=(F1(x2[0],x2[1]+sy,x1[0],x1[1],x0[0],x0[1])-x[4])/sy;
3732 Double_t f42=(F1(x2[0],x2[1],x1[0],x1[1]+sy,x0[0],x0[1])-x[4])/sy;
3733 Double_t f43=(F1(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]+sy)-x[4])/sy;
3734 Double_t f20=(F2(x2[0],x2[1]+sy,x1[0],x1[1],x0[0],x0[1])-x[2])/sy;
3735 Double_t f22=(F2(x2[0],x2[1],x1[0],x1[1]+sy,x0[0],x0[1])-x[2])/sy;
3736 Double_t f23=(F2(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]+sy)-x[2])/sy;
3738 Double_t f30=(F3(x2[0],x2[1]+sy,x0[0],x0[1],x2[2],x0[2])-x[3])/sy;
3739 Double_t f31=(F3(x2[0],x2[1],x0[0],x0[1],x2[2]+sz,x0[2])-x[3])/sz;
3740 Double_t f32=(F3(x2[0],x2[1],x0[0],x0[1]+sy,x2[2],x0[2])-x[3])/sy;
3741 Double_t f34=(F3(x2[0],x2[1],x0[0],x0[1],x2[2],x0[2]+sz)-x[3])/sz;
3746 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3747 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3748 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3749 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3750 c[13]=f30*sy1*f40+f32*sy2*f42;
3751 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3753 // Int_t row1 = fSectors->GetRowNumber(x2[0]);
3754 AliTPCseed *seed=new AliTPCseed(x2[0], sec2*fSectors->GetAlpha()+fSectors->GetAlphaShift(), x, c, 0);
3755 // Double_t y0,z0,y1,z1, y2,z2;
3756 //seed->GetProlongation(x0[0],y0,z0);
3757 // seed->GetProlongation(x1[0],y1,z1);
3758 //seed->GetProlongation(x2[0],y2,z2);
3760 seed->SetLastPoint(pp2);
3761 seed->SetFirstPoint(pp2);
3768 AliTPCseed *AliTPCtrackerMI::ReSeed(AliTPCseed *track, Float_t r0, Float_t r1, Float_t r2)
3772 //reseed using founded clusters
3774 // Find the number of clusters
3775 Int_t nclusters = 0;
3776 for (Int_t irow=0;irow<160;irow++){
3777 if (track->GetClusterIndex(irow)>0) nclusters++;
3781 ipos[0] = TMath::Max(int(r0*nclusters),0); // point 0 cluster
3782 ipos[1] = TMath::Min(int(r1*nclusters),nclusters-1); //
3783 ipos[2] = TMath::Min(int(r2*nclusters),nclusters-1); // last point
3787 Int_t row[3],sec[3]={0,0,0};
3789 // find track row position at given ratio of the length
3791 for (Int_t irow=0;irow<160;irow++){
3792 if (track->GetClusterIndex2(irow)<0) continue;
3794 for (Int_t ipoint=0;ipoint<3;ipoint++){
3795 if (index<=ipos[ipoint]) row[ipoint] = irow;
3799 //Get cluster and sector position
3800 for (Int_t ipoint=0;ipoint<3;ipoint++){
3801 Int_t clindex = track->GetClusterIndex2(row[ipoint]);
3802 AliTPCclusterMI * cl = GetClusterMI(clindex);
3805 // AliTPCclusterMI * cl = GetClusterMI(clindex);
3808 sec[ipoint] = ((clindex&0xff000000)>>24)%18;
3809 xyz[ipoint][0] = GetXrow(row[ipoint]);
3810 xyz[ipoint][1] = cl->GetY();
3811 xyz[ipoint][2] = cl->GetZ();
3815 // Calculate seed state vector and covariance matrix
3817 Double_t alpha, cs,sn, xx2,yy2;
3819 alpha = (sec[1]-sec[2])*fSectors->GetAlpha();
3820 cs = TMath::Cos(alpha);
3821 sn = TMath::Sin(alpha);
3822 xx2= xyz[1][0]*cs-xyz[1][1]*sn;
3823 yy2= xyz[1][0]*sn+xyz[1][1]*cs;
3827 alpha = (sec[0]-sec[2])*fSectors->GetAlpha();
3828 cs = TMath::Cos(alpha);
3829 sn = TMath::Sin(alpha);
3830 xx2= xyz[0][0]*cs-xyz[0][1]*sn;
3831 yy2= xyz[0][0]*sn+xyz[0][1]*cs;
3837 Double_t x[5],c[15];
3841 x[4]=F1(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
3842 x[2]=F2(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
3843 x[3]=F3n(xyz[2][0],xyz[2][1],xyz[0][0],xyz[0][1],xyz[2][2],xyz[0][2],x[4]);
3845 Double_t sy =0.1, sz =0.1;
3847 Double_t sy1=0.2, sz1=0.2;
3848 Double_t sy2=0.2, sz2=0.2;
3851 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;
3852 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;
3853 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;
3854 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;
3855 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;
3856 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;
3858 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;
3859 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;
3860 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;
3861 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;
3866 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3867 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3868 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3869 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3870 c[13]=f30*sy1*f40+f32*sy2*f42;
3871 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3873 // Int_t row1 = fSectors->GetRowNumber(xyz[2][0]);
3874 AliTPCseed *seed=new AliTPCseed(xyz[2][0], sec[2]*fSectors->GetAlpha()+fSectors->GetAlphaShift(), x, c, 0);
3875 seed->SetLastPoint(row[2]);
3876 seed->SetFirstPoint(row[2]);
3881 AliTPCseed *AliTPCtrackerMI::ReSeed(AliTPCseed *track,Int_t r0, Bool_t forward)
3885 //reseed using founded clusters
3888 Int_t row[3]={0,0,0};
3889 Int_t sec[3]={0,0,0};
3891 // forward direction
3893 for (Int_t irow=r0;irow<160;irow++){
3894 if (track->GetClusterIndex(irow)>0){
3899 for (Int_t irow=160;irow>r0;irow--){
3900 if (track->GetClusterIndex(irow)>0){
3905 for (Int_t irow=row[2]-15;irow>row[0];irow--){
3906 if (track->GetClusterIndex(irow)>0){
3914 for (Int_t irow=0;irow<r0;irow++){
3915 if (track->GetClusterIndex(irow)>0){
3920 for (Int_t irow=r0;irow>0;irow--){
3921 if (track->GetClusterIndex(irow)>0){
3926 for (Int_t irow=row[2]-15;irow>row[0];irow--){
3927 if (track->GetClusterIndex(irow)>0){
3934 if ((row[2]-row[0])<20) return 0;
3935 if (row[1]==0) return 0;
3938 //Get cluster and sector position
3939 for (Int_t ipoint=0;ipoint<3;ipoint++){
3940 Int_t clindex = track->GetClusterIndex2(row[ipoint]);
3941 AliTPCclusterMI * cl = GetClusterMI(clindex);
3944 // AliTPCclusterMI * cl = GetClusterMI(clindex);
3947 sec[ipoint] = ((clindex&0xff000000)>>24)%18;
3948 xyz[ipoint][0] = GetXrow(row[ipoint]);
3949 AliTPCTrackerPoint * point = track->GetTrackPoint(row[ipoint]);
3950 if (point&&ipoint<2){
3952 xyz[ipoint][1] = point->GetY();
3953 xyz[ipoint][2] = point->GetZ();
3956 xyz[ipoint][1] = cl->GetY();
3957 xyz[ipoint][2] = cl->GetZ();
3964 // Calculate seed state vector and covariance matrix
3966 Double_t alpha, cs,sn, xx2,yy2;
3968 alpha = (sec[1]-sec[2])*fSectors->GetAlpha();
3969 cs = TMath::Cos(alpha);
3970 sn = TMath::Sin(alpha);
3971 xx2= xyz[1][0]*cs-xyz[1][1]*sn;
3972 yy2= xyz[1][0]*sn+xyz[1][1]*cs;
3976 alpha = (sec[0]-sec[2])*fSectors->GetAlpha();
3977 cs = TMath::Cos(alpha);
3978 sn = TMath::Sin(alpha);
3979 xx2= xyz[0][0]*cs-xyz[0][1]*sn;
3980 yy2= xyz[0][0]*sn+xyz[0][1]*cs;
3986 Double_t x[5],c[15];
3990 x[4]=F1(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
3991 x[2]=F2(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
3992 x[3]=F3n(xyz[2][0],xyz[2][1],xyz[0][0],xyz[0][1],xyz[2][2],xyz[0][2],x[4]);
3994 Double_t sy =0.1, sz =0.1;
3996 Double_t sy1=0.2, sz1=0.2;
3997 Double_t sy2=0.2, sz2=0.2;
4000 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;
4001 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;
4002 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;
4003 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;
4004 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;
4005 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;
4007 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;
4008 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;
4009 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;
4010 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;
4015 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
4016 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
4017 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
4018 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
4019 c[13]=f30*sy1*f40+f32*sy2*f42;
4020 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
4022 // Int_t row1 = fSectors->GetRowNumber(xyz[2][0]);
4023 AliTPCseed *seed=new AliTPCseed(xyz[2][0], sec[2]*fSectors->GetAlpha()+fSectors->GetAlphaShift(), x, c, 0);
4024 seed->SetLastPoint(row[2]);
4025 seed->SetFirstPoint(row[2]);
4026 for (Int_t i=row[0];i<row[2];i++){
4027 seed->SetClusterIndex(i, track->GetClusterIndex(i));
4035 void AliTPCtrackerMI::FindMultiMC(TObjArray * array, AliESDEvent */*esd*/, Int_t iter)
4038 // find multi tracks - THIS FUNCTION IS ONLY FOR DEBUG PURPOSES
4040 // Use AliTPCReconstructor::StreamLevel()>2 if you want to tune parameters - cuts
4042 // Two reasons to have multiple find tracks
4043 // 1. Curling tracks can be find more than once
4044 // 2. Splitted tracks
4045 // a.) Multiple seeding to increase tracking efficiency - (~ 100% reached)
4046 // b.) Edge effect on the sector boundaries
4049 // Algorithm done in 2 phases - because of CPU consumption
4050 // it is n^2 algorithm - for lead-lead 20000x20000 combination are investigated
4052 // Algorihm for curling tracks sign:
4053 // 1 phase -makes a very rough fast cuts to minimize combinatorics
4054 // a.) opposite sign
4055 // b.) one of the tracks - not pointing to the primary vertex -
4056 // c.) delta tan(theta)
4058 // 2 phase - calculates DCA between tracks - time consument
4063 // General cuts - for splitted tracks and for curling tracks
4065 const Float_t kMaxdPhi = 0.2; // maximal distance in phi
4067 // Curling tracks cuts
4072 Int_t nentries = array->GetEntriesFast();
4073 AliHelix *helixes = new AliHelix[nentries];
4074 Float_t *xm = new Float_t[nentries];
4075 Float_t *dz0 = new Float_t[nentries];
4076 Float_t *dz1 = new Float_t[nentries];
4082 // Find track COG in x direction - point with best defined parameters
4084 for (Int_t i=0;i<nentries;i++){
4085 AliTPCseed* track = (AliTPCseed*)array->At(i);
4086 if (!track) continue;
4087 track->SetCircular(0);
4088 new (&helixes[i]) AliHelix(*track);
4092 track->GetDZ(GetX(),GetY(),GetZ(),GetBz(),dz);
4095 for (Int_t icl=0; icl<160; icl++){
4096 AliTPCclusterMI * cl = track->GetClusterPointer(icl);
4102 if (ncl>0) xm[i]/=Float_t(ncl);
4104 TTreeSRedirector &cstream = *fDebugStreamer;
4106 for (Int_t i0=0;i0<nentries;i0++){
4107 AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
4108 if (!track0) continue;
4109 Float_t xc0 = helixes[i0].GetHelix(6);
4110 Float_t yc0 = helixes[i0].GetHelix(7);
4111 Float_t r0 = helixes[i0].GetHelix(8);
4112 Float_t rc0 = TMath::Sqrt(xc0*xc0+yc0*yc0);
4113 Float_t fi0 = TMath::ATan2(yc0,xc0);
4115 for (Int_t i1=i0+1;i1<nentries;i1++){
4116 AliTPCseed * track1 = (AliTPCseed*)array->At(i1);
4117 if (!track1) continue;
4118 Int_t lab0=track0->GetLabel();
4119 Int_t lab1=track1->GetLabel();
4120 if (TMath::Abs(lab0)!=TMath::Abs(lab1)) continue;
4122 Float_t xc1 = helixes[i1].GetHelix(6);
4123 Float_t yc1 = helixes[i1].GetHelix(7);
4124 Float_t r1 = helixes[i1].GetHelix(8);
4125 Float_t rc1 = TMath::Sqrt(xc1*xc1+yc1*yc1);
4126 Float_t fi1 = TMath::ATan2(yc1,xc1);
4128 Float_t dfi = fi0-fi1;
4131 if (dfi>1.5*TMath::Pi()) dfi-=TMath::Pi(); // take care about edge effect
4132 if (dfi<-1.5*TMath::Pi()) dfi+=TMath::Pi(); //
4133 if (TMath::Abs(dfi)>kMaxdPhi&&helixes[i0].GetHelix(4)*helixes[i1].GetHelix(4)<0){
4135 // if short tracks with undefined sign
4136 fi1 = -TMath::ATan2(yc1,-xc1);
4139 Float_t dtheta = TMath::Abs(track0->GetTgl()-track1->GetTgl())<TMath::Abs(track0->GetTgl()+track1->GetTgl())? track0->GetTgl()-track1->GetTgl():track0->GetTgl()+track1->GetTgl();
4142 // debug stream to tune "fast cuts"
4144 Double_t dist[3]; // distance at X
4145 Double_t mdist[3]={0,0,0}; // mean distance X+-40cm
4146 track0->GetDistance(track1,0.5*(xm[i0]+xm[i1])-40.,dist,AliTracker::GetBz());
4147 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[i]);
4148 track0->GetDistance(track1,0.5*(xm[i0]+xm[i1])+40.,dist,AliTracker::GetBz());
4149 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[i]);
4150 track0->GetDistance(track1,0.5*(xm[i0]+xm[i1]),dist,AliTracker::GetBz());
4151 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[i]);
4152 for (Int_t i=0;i<3;i++) mdist[i]*=0.33333;
4156 for (Int_t icl=0; icl<160; icl++){
4157 AliTPCclusterMI * cl0 = track0->GetClusterPointer(icl);
4158 AliTPCclusterMI * cl1 = track1->GetClusterPointer(icl);
4161 if (cl0==cl1) sums++;
4169 "Tr0.="<<track0<< // seed0
4170 "Tr1.="<<track1<< // seed1
4171 "h0.="<<&helixes[i0]<<
4172 "h1.="<<&helixes[i1]<<
4174 "sum="<<sum<< //the sum of rows with cl in both
4175 "sums="<<sums<< //the sum of shared clusters
4176 "xm0="<<xm[i0]<< // the center of track
4177 "xm1="<<xm[i1]<< // the x center of track
4178 // General cut variables
4179 "dfi="<<dfi<< // distance in fi angle
4180 "dtheta="<<dtheta<< // distance int theta angle
4186 "dist0="<<dist[0]<< //distance x
4187 "dist1="<<dist[1]<< //distance y
4188 "dist2="<<dist[2]<< //distance z
4189 "mdist0="<<mdist[0]<< //distance x
4190 "mdist1="<<mdist[1]<< //distance y
4191 "mdist2="<<mdist[2]<< //distance z
4204 if (AliTPCReconstructor::StreamLevel()>1) {
4205 AliInfo("Time for curling tracks removal DEBUGGING MC");
4211 void AliTPCtrackerMI::FindSplitted(TObjArray * array, AliESDEvent */*esd*/, Int_t iter)
4215 // Two reasons to have multiple find tracks
4216 // 1. Curling tracks can be find more than once
4217 // 2. Splitted tracks
4218 // a.) Multiple seeding to increase tracking efficiency - (~ 100% reached)
4219 // b.) Edge effect on the sector boundaries
4221 // This function tries to find tracks closed in the parametric space
4223 // cut logic if distance is bigger than cut continue - Do Nothing
4224 const Float_t kMaxdTheta = 0.05; // maximal distance in theta
4225 const Float_t kMaxdPhi = 0.05; // maximal deistance in phi
4226 const Float_t kdelta = 40.; //delta r to calculate track distance
4228 // const Float_t kMaxDist0 = 1.; // maximal distance 0
4229 //const Float_t kMaxDist1 = 0.3; // maximal distance 1 - cut if track in separate rows
4232 TCut csec("csec","abs(Tr0.fRelativeSector-Tr1.fRelativeSector)<2");
4233 TCut cdtheta("cdtheta","abs(dtheta)<0.05");
4238 Int_t nentries = array->GetEntriesFast();
4239 AliHelix *helixes = new AliHelix[nentries];
4240 Float_t *xm = new Float_t[nentries];
4246 //Sort tracks according quality
4248 Int_t nseed = array->GetEntriesFast();
4249 Float_t * quality = new Float_t[nseed];
4250 Int_t * indexes = new Int_t[nseed];
4251 for (Int_t i=0; i<nseed; i++) {
4252 AliTPCseed *pt=(AliTPCseed*)array->UncheckedAt(i);
4257 pt->UpdatePoints(); //select first last max dens points
4258 Float_t * points = pt->GetPoints();
4259 if (points[3]<0.8) quality[i] =-1;
4260 quality[i] = (points[2]-points[0])+pt->GetNumberOfClusters();
4261 //prefer high momenta tracks if overlaps
4262 quality[i] *= TMath::Sqrt(TMath::Abs(pt->Pt())+0.5);
4264 TMath::Sort(nseed,quality,indexes);
4268 // Find track COG in x direction - point with best defined parameters
4270 for (Int_t i=0;i<nentries;i++){
4271 AliTPCseed* track = (AliTPCseed*)array->At(i);
4272 if (!track) continue;
4273 track->SetCircular(0);
4274 new (&helixes[i]) AliHelix(*track);
4277 for (Int_t icl=0; icl<160; icl++){
4278 AliTPCclusterMI * cl = track->GetClusterPointer(icl);
4284 if (ncl>0) xm[i]/=Float_t(ncl);
4286 TTreeSRedirector &cstream = *fDebugStreamer;
4288 for (Int_t is0=0;is0<nentries;is0++){
4289 Int_t i0 = indexes[is0];
4290 AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
4291 if (!track0) continue;
4292 if (track0->GetKinkIndexes()[0]!=0) continue;
4293 Float_t xc0 = helixes[i0].GetHelix(6);
4294 Float_t yc0 = helixes[i0].GetHelix(7);
4295 Float_t fi0 = TMath::ATan2(yc0,xc0);
4297 for (Int_t is1=is0+1;is1<nentries;is1++){
4298 Int_t i1 = indexes[is1];
4299 AliTPCseed * track1 = (AliTPCseed*)array->At(i1);
4300 if (!track1) continue;
4302 if (TMath::Abs(track0->GetRelativeSector()-track1->GetRelativeSector())>1) continue;
4303 if (track1->GetKinkIndexes()[0]>0 &&track0->GetKinkIndexes()[0]<0) continue;
4304 if (track1->GetKinkIndexes()[0]!=0) continue;
4306 Float_t dtheta = TMath::Abs(track0->GetTgl()-track1->GetTgl())<TMath::Abs(track0->GetTgl()+track1->GetTgl())? track0->GetTgl()-track1->GetTgl():track0->GetTgl()+track1->GetTgl();
4307 if (TMath::Abs(dtheta)>kMaxdTheta) continue;
4309 Float_t xc1 = helixes[i1].GetHelix(6);
4310 Float_t yc1 = helixes[i1].GetHelix(7);
4311 Float_t fi1 = TMath::ATan2(yc1,xc1);
4313 Float_t dfi = fi0-fi1;
4314 if (dfi>1.5*TMath::Pi()) dfi-=TMath::Pi(); // take care about edge effect
4315 if (dfi<-1.5*TMath::Pi()) dfi+=TMath::Pi(); //
4316 if (TMath::Abs(dfi)>kMaxdPhi&&helixes[i0].GetHelix(4)*helixes[i1].GetHelix(4)<0){
4318 // if short tracks with undefined sign
4319 fi1 = -TMath::ATan2(yc1,-xc1);
4322 if (TMath::Abs(dfi)>kMaxdPhi) continue;
4329 for (Int_t icl=0; icl<160; icl++){
4330 Int_t index0=track0->GetClusterIndex2(icl);
4331 Int_t index1=track1->GetClusterIndex2(icl);
4332 Bool_t used0 = (index0>0 && !(index0&0x8000));
4333 Bool_t used1 = (index1>0 && !(index1&0x8000));
4335 if (used0) sum0++; // used cluster0
4336 if (used1) sum1++; // used clusters1
4337 if (used0&&used1) sum++;
4338 if (index0==index1 && used0 && used1) sums++;
4342 if (sums<10) continue;
4343 if (sum<40) continue;
4344 if (sums/Float_t(TMath::Min(sum0,sum1))<0.5) continue;
4346 Double_t dist[5][4]; // distance at X
4347 Double_t mdist[4]={0,0,0,0}; // mean distance on range +-delta
4351 track0->GetDistance(track1,xm[i0],dist[0],AliTracker::GetBz());
4352 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[0][i]);
4353 track0->GetDistance(track1,xm[i1],dist[1],AliTracker::GetBz());
4354 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[1][i]);
4356 track0->GetDistance(track1,TMath::Min(xm[i1],xm[i0])-kdelta,dist[2],AliTracker::GetBz());
4357 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[2][i]);
4358 track0->GetDistance(track1,TMath::Max(xm[i1],xm[i0])+kdelta,dist[3],AliTracker::GetBz());
4359 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[3][i]);
4361 track0->GetDistance(track1,(xm[i1]+xm[i0])*0.5,dist[4],AliTracker::GetBz());
4362 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[4][i]);
4363 for (Int_t i=0;i<3;i++) mdist[i]*=0.2;
4366 Int_t lab0=track0->GetLabel();
4367 Int_t lab1=track1->GetLabel();
4368 cstream<<"Splitted"<<
4372 "Tr0.="<<track0<< // seed0
4373 "Tr1.="<<track1<< // seed1
4374 "h0.="<<&helixes[i0]<<
4375 "h1.="<<&helixes[i1]<<
4377 "sum="<<sum<< //the sum of rows with cl in both
4378 "sum0="<<sum0<< //the sum of rows with cl in 0 track
4379 "sum1="<<sum1<< //the sum of rows with cl in 1 track
4380 "sums="<<sums<< //the sum of shared clusters
4381 "xm0="<<xm[i0]<< // the center of track
4382 "xm1="<<xm[i1]<< // the x center of track
4383 // General cut variables
4384 "dfi="<<dfi<< // distance in fi angle
4385 "dtheta="<<dtheta<< // distance int theta angle
4388 "dist0="<<dist[4][0]<< //distance x
4389 "dist1="<<dist[4][1]<< //distance y
4390 "dist2="<<dist[4][2]<< //distance z
4391 "mdist0="<<mdist[0]<< //distance x
4392 "mdist1="<<mdist[1]<< //distance y
4393 "mdist2="<<mdist[2]<< //distance z
4396 delete array->RemoveAt(i1);
4401 AliInfo("Time for splitted tracks removal");
4407 void AliTPCtrackerMI::FindCurling(TObjArray * array, AliESDEvent */*esd*/, Int_t iter)
4410 // find Curling tracks
4411 // Use AliTPCReconstructor::StreamLevel()>1 if you want to tune parameters - cuts
4414 // Algorithm done in 2 phases - because of CPU consumption
4415 // it is n^2 algorithm - for lead-lead 20000x20000 combination are investigated
4416 // see detal in MC part what can be used to cut
4420 const Float_t kMaxC = 400; // maximal curvature to of the track
4421 const Float_t kMaxdTheta = 0.15; // maximal distance in theta
4422 const Float_t kMaxdPhi = 0.15; // maximal distance in phi
4423 const Float_t kPtRatio = 0.3; // ratio between pt
4424 const Float_t kMinDCAR = 2.; // distance to the primary vertex in r - see cpipe cut
4427 // Curling tracks cuts
4430 const Float_t kMaxDeltaRMax = 40; // distance in outer radius
4431 const Float_t kMaxDeltaRMin = 5.; // distance in lower radius - see cpipe cut
4432 const Float_t kMinAngle = 2.9; // angle between tracks
4433 const Float_t kMaxDist = 5; // biggest distance
4435 // The cuts can be tuned using the "MC information stored in Multi tree ==> see FindMultiMC
4438 TCut csign("csign","Tr0.fP[4]*Tr1.fP[4]<0"); //opposite sign
4439 TCut cmax("cmax","abs(Tr0.GetC())>1/400");
4440 TCut cda("cda","sqrt(dtheta^2+dfi^2)<0.15");
4441 TCut ccratio("ccratio","abs((Tr0.fP[4]+Tr1.fP[4])/(abs(Tr0.fP[4])+abs(Tr1.fP[4])))<0.3");
4442 TCut cpipe("cpipe", "min(abs(r0-rc0),abs(r1-rc1))>5");
4444 TCut cdrmax("cdrmax","abs(abs(rc0+r0)-abs(rc1+r1))<40")
4445 TCut cdrmin("cdrmin","abs(abs(rc0+r0)-abs(rc1+r1))<10")
4447 Multi->Draw("dfi","iter==0"+csign+cmax+cda+ccratio); ~94% of curling tracks fulfill
4448 Multi->Draw("min(abs(r0-rc0),abs(r1-rc1))","iter==0&&abs(lab1)==abs(lab0)"+csign+cmax+cda+ccratio+cpipe+cdrmin+cdrmax); //80%
4450 Curling2->Draw("dfi","iter==0&&abs(lab0)==abs(lab1)"+csign+cmax+cdtheta+cdfi+ccratio)
4456 Int_t nentries = array->GetEntriesFast();
4457 AliHelix *helixes = new AliHelix[nentries];
4458 for (Int_t i=0;i<nentries;i++){
4459 AliTPCseed* track = (AliTPCseed*)array->At(i);
4460 if (!track) continue;
4461 track->SetCircular(0);
4462 new (&helixes[i]) AliHelix(*track);
4468 Double_t phase[2][2],radius[2];
4472 TTreeSRedirector &cstream = *fDebugStreamer;
4474 for (Int_t i0=0;i0<nentries;i0++){
4475 AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
4476 if (!track0) continue;
4477 if (TMath::Abs(track0->GetC())<1/kMaxC) continue;
4478 Float_t xc0 = helixes[i0].GetHelix(6);
4479 Float_t yc0 = helixes[i0].GetHelix(7);
4480 Float_t r0 = helixes[i0].GetHelix(8);
4481 Float_t rc0 = TMath::Sqrt(xc0*xc0+yc0*yc0);
4482 Float_t fi0 = TMath::ATan2(yc0,xc0);
4484 for (Int_t i1=i0+1;i1<nentries;i1++){
4485 AliTPCseed * track1 = (AliTPCseed*)array->At(i1);
4486 if (!track1) continue;
4487 if (TMath::Abs(track1->GetC())<1/kMaxC) continue;
4488 Float_t xc1 = helixes[i1].GetHelix(6);
4489 Float_t yc1 = helixes[i1].GetHelix(7);
4490 Float_t r1 = helixes[i1].GetHelix(8);
4491 Float_t rc1 = TMath::Sqrt(xc1*xc1+yc1*yc1);
4492 Float_t fi1 = TMath::ATan2(yc1,xc1);
4494 Float_t dfi = fi0-fi1;
4497 if (dfi>1.5*TMath::Pi()) dfi-=TMath::Pi(); // take care about edge effect
4498 if (dfi<-1.5*TMath::Pi()) dfi+=TMath::Pi(); //
4499 Float_t dtheta = TMath::Abs(track0->GetTgl()-track1->GetTgl())<TMath::Abs(track0->GetTgl()+track1->GetTgl())? track0->GetTgl()-track1->GetTgl():track0->GetTgl()+track1->GetTgl();
4503 if (track0->GetBConstrain()&&track1->GetBConstrain()) continue; // not constrained
4504 if (track1->GetSigned1Pt()*track0->GetSigned1Pt()>0) continue; // not the same sign
4505 if ( TMath::Abs(track1->GetTgl()+track0->GetTgl())>kMaxdTheta) continue; //distance in the Theta
4506 if ( TMath::Abs(dfi)>kMaxdPhi) continue; //distance in phi
4507 if ( TMath::Sqrt(dfi*dfi+dtheta*dtheta)>kMaxdPhi) continue; //common angular offset
4509 Float_t pt0 = track0->GetSignedPt();
4510 Float_t pt1 = track1->GetSignedPt();
4511 if ((TMath::Abs(pt0+pt1)/(TMath::Abs(pt0)+TMath::Abs(pt1)))>kPtRatio) continue;
4512 if ((iter==1) && TMath::Abs(TMath::Abs(rc0+r0)-TMath::Abs(rc1+r1))>kMaxDeltaRMax) continue;
4513 if ((iter!=1) &&TMath::Abs(TMath::Abs(rc0-r0)-TMath::Abs(rc1-r1))>kMaxDeltaRMin) continue;
4514 if (TMath::Min(TMath::Abs(rc0-r0),TMath::Abs(rc1-r1))<kMinDCAR) continue;
4517 // Now find closest approach
4521 Int_t npoints = helixes[i0].GetRPHIintersections(helixes[i1], phase, radius,10);
4522 if (npoints==0) continue;
4523 helixes[i0].GetClosestPhases(helixes[i1], phase);
4527 Double_t hangles[3];
4528 helixes[i0].Evaluate(phase[0][0],xyz0);
4529 helixes[i1].Evaluate(phase[0][1],xyz1);
4531 helixes[i0].GetAngle(phase[0][0],helixes[i1],phase[0][1],hangles);
4532 Double_t deltah[2],deltabest;
4533 if (TMath::Abs(hangles[2])<kMinAngle) continue;
4537 helixes[i0].ParabolicDCA(helixes[i1],phase[0][0],phase[0][1],radius[0],deltah[0],2);
4539 helixes[i0].ParabolicDCA(helixes[i1],phase[1][0],phase[1][1],radius[1],deltah[1],2);
4540 if (deltah[1]<deltah[0]) ibest=1;
4542 deltabest = TMath::Sqrt(deltah[ibest]);
4543 helixes[i0].Evaluate(phase[ibest][0],xyz0);
4544 helixes[i1].Evaluate(phase[ibest][1],xyz1);
4545 helixes[i0].GetAngle(phase[ibest][0],helixes[i1],phase[ibest][1],hangles);
4546 Double_t radiusbest = TMath::Sqrt(radius[ibest]);
4548 if (deltabest>kMaxDist) continue;
4549 // if (mindcar+mindcaz<40 && (TMath::Abs(hangles[2])<kMinAngle ||deltabest>3)) continue;
4550 Bool_t sign =kFALSE;
4551 if (hangles[2]>kMinAngle) sign =kTRUE;
4554 // circular[i0] = kTRUE;
4555 // circular[i1] = kTRUE;
4556 if (track0->OneOverPt()<track1->OneOverPt()){
4557 track0->SetCircular(track0->GetCircular()+1);
4558 track1->SetCircular(track1->GetCircular()+2);
4561 track1->SetCircular(track1->GetCircular()+1);
4562 track0->SetCircular(track0->GetCircular()+2);
4565 if (AliTPCReconstructor::StreamLevel()>1){
4567 //debug stream to tune "fine" cuts
4568 Int_t lab0=track0->GetLabel();
4569 Int_t lab1=track1->GetLabel();
4570 cstream<<"Curling2"<<
4586 "npoints="<<npoints<<
4587 "hangles0="<<hangles[0]<<
4588 "hangles1="<<hangles[1]<<
4589 "hangles2="<<hangles[2]<<
4592 "radius="<<radiusbest<<
4593 "deltabest="<<deltabest<<
4594 "phase0="<<phase[ibest][0]<<
4595 "phase1="<<phase[ibest][1]<<
4603 if (AliTPCReconstructor::StreamLevel()>1) {
4604 AliInfo("Time for curling tracks removal");
4613 void AliTPCtrackerMI::FindKinks(TObjArray * array, AliESDEvent *esd)
4620 TObjArray *kinks= new TObjArray(10000);
4621 // TObjArray *v0s= new TObjArray(10000);
4622 Int_t nentries = array->GetEntriesFast();
4623 AliHelix *helixes = new AliHelix[nentries];
4624 Int_t *sign = new Int_t[nentries];
4625 Int_t *nclusters = new Int_t[nentries];
4626 Float_t *alpha = new Float_t[nentries];
4627 AliKink *kink = new AliKink();
4628 Int_t * usage = new Int_t[nentries];
4629 Float_t *zm = new Float_t[nentries];
4630 Float_t *z0 = new Float_t[nentries];
4631 Float_t *fim = new Float_t[nentries];
4632 Float_t *shared = new Float_t[nentries];
4633 Bool_t *circular = new Bool_t[nentries];
4634 Float_t *dca = new Float_t[nentries];
4635 //const AliESDVertex * primvertex = esd->GetVertex();
4637 // nentries = array->GetEntriesFast();
4642 for (Int_t i=0;i<nentries;i++){
4645 AliTPCseed* track = (AliTPCseed*)array->At(i);
4646 if (!track) continue;
4647 track->SetCircular(0);
4649 track->UpdatePoints();
4650 if (( track->GetPoints()[2]- track->GetPoints()[0])>5 && track->GetPoints()[3]>0.8){
4652 nclusters[i]=track->GetNumberOfClusters();
4653 alpha[i] = track->GetAlpha();
4654 new (&helixes[i]) AliHelix(*track);
4656 helixes[i].Evaluate(0,xyz);
4657 sign[i] = (track->GetC()>0) ? -1:1;
4660 if (track->GetProlongation(x,y,z)){
4662 fim[i] = alpha[i]+TMath::ATan2(y,x);
4665 zm[i] = track->GetZ();
4669 circular[i]= kFALSE;
4670 if (track->GetProlongation(0,y,z)) z0[i] = z;
4671 dca[i] = track->GetD(0,0);
4677 Int_t ncandidates =0;
4680 Double_t phase[2][2],radius[2];
4683 // Find circling track
4684 TTreeSRedirector &cstream = *fDebugStreamer;
4686 for (Int_t i0=0;i0<nentries;i0++){
4687 AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
4688 if (!track0) continue;
4689 if (track0->GetNumberOfClusters()<40) continue;
4690 if (TMath::Abs(1./track0->GetC())>200) continue;
4691 for (Int_t i1=i0+1;i1<nentries;i1++){
4692 AliTPCseed * track1 = (AliTPCseed*)array->At(i1);
4693 if (!track1) continue;
4694 if (track1->GetNumberOfClusters()<40) continue;
4695 if ( TMath::Abs(track1->GetTgl()+track0->GetTgl())>0.1) continue;
4696 if (track0->GetBConstrain()&&track1->GetBConstrain()) continue;
4697 if (TMath::Abs(1./track1->GetC())>200) continue;
4698 if (track1->GetSigned1Pt()*track0->GetSigned1Pt()>0) continue;
4699 if (track1->GetTgl()*track0->GetTgl()>0) continue;
4700 if (TMath::Max(TMath::Abs(1./track0->GetC()),TMath::Abs(1./track1->GetC()))>190) continue;
4701 if (track0->GetBConstrain()&&track1->OneOverPt()<track0->OneOverPt()) continue; //returning - lower momenta
4702 if (track1->GetBConstrain()&&track0->OneOverPt()<track1->OneOverPt()) continue; //returning - lower momenta
4704 Float_t mindcar = TMath::Min(TMath::Abs(dca[i0]),TMath::Abs(dca[i1]));
4705 if (mindcar<5) continue;
4706 Float_t mindcaz = TMath::Min(TMath::Abs(z0[i0]-GetZ()),TMath::Abs(z0[i1]-GetZ()));
4707 if (mindcaz<5) continue;
4708 if (mindcar+mindcaz<20) continue;
4711 Float_t xc0 = helixes[i0].GetHelix(6);
4712 Float_t yc0 = helixes[i0].GetHelix(7);
4713 Float_t r0 = helixes[i0].GetHelix(8);
4714 Float_t xc1 = helixes[i1].GetHelix(6);
4715 Float_t yc1 = helixes[i1].GetHelix(7);
4716 Float_t r1 = helixes[i1].GetHelix(8);
4718 Float_t rmean = (r0+r1)*0.5;
4719 Float_t delta =TMath::Sqrt((xc1-xc0)*(xc1-xc0)+(yc1-yc0)*(yc1-yc0));
4720 //if (delta>30) continue;
4721 if (delta>rmean*0.25) continue;
4722 if (TMath::Abs(r0-r1)/rmean>0.3) continue;
4724 Int_t npoints = helixes[i0].GetRPHIintersections(helixes[i1], phase, radius,10);
4725 if (npoints==0) continue;
4726 helixes[i0].GetClosestPhases(helixes[i1], phase);
4730 Double_t hangles[3];
4731 helixes[i0].Evaluate(phase[0][0],xyz0);
4732 helixes[i1].Evaluate(phase[0][1],xyz1);
4734 helixes[i0].GetAngle(phase[0][0],helixes[i1],phase[0][1],hangles);
4735 Double_t deltah[2],deltabest;
4736 if (hangles[2]<2.8) continue;
4739 helixes[i0].ParabolicDCA(helixes[i1],phase[0][0],phase[0][1],radius[0],deltah[0],2);
4741 helixes[i0].ParabolicDCA(helixes[i1],phase[1][0],phase[1][1],radius[1],deltah[1],2);
4742 if (deltah[1]<deltah[0]) ibest=1;
4744 deltabest = TMath::Sqrt(deltah[ibest]);
4745 helixes[i0].Evaluate(phase[ibest][0],xyz0);
4746 helixes[i1].Evaluate(phase[ibest][1],xyz1);
4747 helixes[i0].GetAngle(phase[ibest][0],helixes[i1],phase[ibest][1],hangles);
4748 Double_t radiusbest = TMath::Sqrt(radius[ibest]);
4750 if (deltabest>6) continue;
4751 if (mindcar+mindcaz<40 && (hangles[2]<3.12||deltabest>3)) continue;
4752 Bool_t sign =kFALSE;
4753 if (hangles[2]>3.06) sign =kTRUE;
4756 circular[i0] = kTRUE;
4757 circular[i1] = kTRUE;
4758 if (track0->OneOverPt()<track1->OneOverPt()){
4759 track0->SetCircular(track0->GetCircular()+1);
4760 track1->SetCircular(track1->GetCircular()+2);
4763 track1->SetCircular(track1->GetCircular()+1);
4764 track0->SetCircular(track0->GetCircular()+2);
4767 if (sign&&AliTPCReconstructor::StreamLevel()>1){
4769 Int_t lab0=track0->GetLabel();
4770 Int_t lab1=track1->GetLabel();
4771 cstream<<"Curling"<<
4778 "mindcar="<<mindcar<<
4779 "mindcaz="<<mindcaz<<
4782 "npoints="<<npoints<<
4783 "hangles0="<<hangles[0]<<
4784 "hangles2="<<hangles[2]<<
4789 "radius="<<radiusbest<<
4790 "deltabest="<<deltabest<<
4791 "phase0="<<phase[ibest][0]<<
4792 "phase1="<<phase[ibest][1]<<
4802 for (Int_t i =0;i<nentries;i++){
4803 if (sign[i]==0) continue;
4804 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
4811 Double_t cradius0 = 40*40;
4812 Double_t cradius1 = 270*270;
4815 Double_t cdist3=0.55;
4816 for (Int_t j =i+1;j<nentries;j++){
4818 if (sign[j]*sign[i]<1) continue;
4819 if ( (nclusters[i]+nclusters[j])>200) continue;
4820 if ( (nclusters[i]+nclusters[j])<80) continue;
4821 if ( TMath::Abs(zm[i]-zm[j])>60.) continue;
4822 if ( TMath::Abs(fim[i]-fim[j])>0.6 && TMath::Abs(fim[i]-fim[j])<5.7 ) continue;
4823 //AliTPCseed * track1 = (AliTPCseed*)array->At(j); Double_t phase[2][2],radius[2];
4824 Int_t npoints = helixes[i].GetRPHIintersections(helixes[j], phase, radius,20);
4825 if (npoints<1) continue;
4828 if (radius[0]<cradius0||radius[0]>cradius1) continue;
4831 if ( (radius[0]<cradius0||radius[0]>cradius1) && (radius[1]<cradius0||radius[1]>cradius1) ) continue;
4834 Double_t delta1=10000,delta2=10000;
4835 // cuts on the intersection radius
4836 helixes[i].LinearDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta1);
4837 if (radius[0]<20&&delta1<1) continue; //intersection at vertex
4838 if (radius[0]<10&&delta1<3) continue; //intersection at vertex
4840 helixes[i].LinearDCA(helixes[j],phase[1][0],phase[1][1],radius[1],delta2);
4841 if (radius[1]<20&&delta2<1) continue; //intersection at vertex
4842 if (radius[1]<10&&delta2<3) continue; //intersection at vertex
4845 Double_t distance1 = TMath::Min(delta1,delta2);
4846 if (distance1>cdist1) continue; // cut on DCA linear approximation
4848 npoints = helixes[i].GetRPHIintersections(helixes[j], phase, radius,20);
4849 helixes[i].ParabolicDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta1);
4850 if (radius[0]<20&&delta1<1) continue; //intersection at vertex
4851 if (radius[0]<10&&delta1<3) continue; //intersection at vertex
4854 helixes[i].ParabolicDCA(helixes[j],phase[1][0],phase[1][1],radius[1],delta2);
4855 if (radius[1]<20&&delta2<1) continue; //intersection at vertex
4856 if (radius[1]<10&&delta2<3) continue; //intersection at vertex
4858 distance1 = TMath::Min(delta1,delta2);
4861 rkink = TMath::Sqrt(radius[0]);
4864 rkink = TMath::Sqrt(radius[1]);
4866 if (distance1>cdist2) continue;
4869 AliTPCseed * track1 = (AliTPCseed*)array->At(j);
4872 Int_t row0 = GetRowNumber(rkink);
4873 if (row0<10) continue;
4874 if (row0>150) continue;
4877 Float_t dens00=-1,dens01=-1;
4878 Float_t dens10=-1,dens11=-1;
4880 Int_t found,foundable,shared;
4881 track0->GetClusterStatistic(0,row0-5, found, foundable,shared,kFALSE);
4882 if (foundable>5) dens00 = Float_t(found)/Float_t(foundable);
4883 track0->GetClusterStatistic(row0+5,155, found, foundable,shared,kFALSE);
4884 if (foundable>5) dens01 = Float_t(found)/Float_t(foundable);
4886 track1->GetClusterStatistic(0,row0-5, found, foundable,shared,kFALSE);
4887 if (foundable>10) dens10 = Float_t(found)/Float_t(foundable);
4888 track1->GetClusterStatistic(row0+5,155, found, foundable,shared,kFALSE);
4889 if (foundable>10) dens11 = Float_t(found)/Float_t(foundable);
4891 if (dens00<dens10 && dens01<dens11) continue;
4892 if (dens00>dens10 && dens01>dens11) continue;
4893 if (TMath::Max(dens00,dens10)<0.1) continue;
4894 if (TMath::Max(dens01,dens11)<0.3) continue;
4896 if (TMath::Min(dens00,dens10)>0.6) continue;
4897 if (TMath::Min(dens01,dens11)>0.6) continue;
4900 AliTPCseed * ktrack0, *ktrack1;
4909 if (TMath::Abs(ktrack0->GetC())>5) continue; // cut on the curvature for mother particle
4910 AliExternalTrackParam paramm(*ktrack0);
4911 AliExternalTrackParam paramd(*ktrack1);
4912 if (row0>60&&ktrack1->GetReference().GetX()>90.)new (¶md) AliExternalTrackParam(ktrack1->GetReference());
4915 kink->SetMother(paramm);
4916 kink->SetDaughter(paramd);
4919 Float_t x[3] = { kink->GetPosition()[0],kink->GetPosition()[1],kink->GetPosition()[2]};
4921 fParam->Transform0to1(x,index);
4922 fParam->Transform1to2(x,index);
4923 row0 = GetRowNumber(x[0]);
4925 if (kink->GetR()<100) continue;
4926 if (kink->GetR()>240) continue;
4927 if (kink->GetPosition()[2]/kink->GetR()>AliTPCReconstructor::GetCtgRange()) continue; //out of fiducial volume
4928 if (kink->GetDistance()>cdist3) continue;
4929 Float_t dird = kink->GetDaughterP()[0]*kink->GetPosition()[0]+kink->GetDaughterP()[1]*kink->GetPosition()[1]; // rough direction estimate
4930 if (dird<0) continue;
4932 Float_t dirm = kink->GetMotherP()[0]*kink->GetPosition()[0]+kink->GetMotherP()[1]*kink->GetPosition()[1]; // rough direction estimate
4933 if (dirm<0) continue;
4934 Float_t mpt = TMath::Sqrt(kink->GetMotherP()[0]*kink->GetMotherP()[0]+kink->GetMotherP()[1]*kink->GetMotherP()[1]);
4935 if (mpt<0.2) continue;
4938 //for high momenta momentum not defined well in first iteration
4939 Double_t qt = TMath::Sin(kink->GetAngle(2))*ktrack1->GetP();
4940 if (qt>0.35) continue;
4943 kink->SetLabel(CookLabel(ktrack0,0.4,0,row0),0);
4944 kink->SetLabel(CookLabel(ktrack1,0.4,row0,160),1);
4946 kink->SetTPCDensity(dens00,0,0);
4947 kink->SetTPCDensity(dens01,0,1);
4948 kink->SetTPCDensity(dens10,1,0);
4949 kink->SetTPCDensity(dens11,1,1);
4950 kink->SetIndex(i,0);
4951 kink->SetIndex(j,1);
4954 kink->SetTPCDensity(dens10,0,0);
4955 kink->SetTPCDensity(dens11,0,1);
4956 kink->SetTPCDensity(dens00,1,0);
4957 kink->SetTPCDensity(dens01,1,1);
4958 kink->SetIndex(j,0);
4959 kink->SetIndex(i,1);
4962 if (mpt<1||kink->GetAngle(2)>0.1){
4963 // angle and densities not defined yet
4964 if (kink->GetTPCDensityFactor()<0.8) continue;
4965 if ((2-kink->GetTPCDensityFactor())*kink->GetDistance() >0.25) continue;
4966 if (kink->GetAngle(2)*ktrack0->GetP()<0.003) continue; //too small angle
4967 if (kink->GetAngle(2)>0.2&&kink->GetTPCDensityFactor()<1.15) continue;
4968 if (kink->GetAngle(2)>0.2&&kink->GetTPCDensity(0,1)>0.05) continue;
4970 Float_t criticalangle = track0->GetSigmaSnp2()+track0->GetSigmaTgl2();
4971 criticalangle+= track1->GetSigmaSnp2()+track1->GetSigmaTgl2();
4972 criticalangle= 3*TMath::Sqrt(criticalangle);
4973 if (criticalangle>0.02) criticalangle=0.02;
4974 if (kink->GetAngle(2)<criticalangle) continue;
4977 Int_t drow = Int_t(2.+0.5/(0.05+kink->GetAngle(2))); // overlap region defined
4978 Float_t shapesum =0;
4980 for ( Int_t row = row0-drow; row<row0+drow;row++){
4981 if (row<0) continue;
4982 if (row>155) continue;
4983 if (ktrack0->GetClusterPointer(row)){
4984 AliTPCTrackerPoint *point =ktrack0->GetTrackPoint(row);
4985 shapesum+=point->GetSigmaY()+point->GetSigmaZ();
4988 if (ktrack1->GetClusterPointer(row)){
4989 AliTPCTrackerPoint *point =ktrack1->GetTrackPoint(row);
4990 shapesum+=point->GetSigmaY()+point->GetSigmaZ();
4995 kink->SetShapeFactor(-1.);
4998 kink->SetShapeFactor(shapesum/sum);
5000 // esd->AddKink(kink);
5001 kinks->AddLast(kink);
5007 // sort the kinks according quality - and refit them towards vertex
5009 Int_t nkinks = kinks->GetEntriesFast();
5010 Float_t *quality = new Float_t[nkinks];
5011 Int_t *indexes = new Int_t[nkinks];
5012 AliTPCseed *mothers = new AliTPCseed[nkinks];
5013 AliTPCseed *daughters = new AliTPCseed[nkinks];
5016 for (Int_t i=0;i<nkinks;i++){
5018 AliKink *kink = (AliKink*)kinks->At(i);
5020 // refit kinks towards vertex
5022 Int_t index0 = kink->GetIndex(0);
5023 Int_t index1 = kink->GetIndex(1);
5024 AliTPCseed * ktrack0 = (AliTPCseed*)array->At(index0);
5025 AliTPCseed * ktrack1 = (AliTPCseed*)array->At(index1);
5027 Int_t sumn=ktrack0->GetNumberOfClusters()+ktrack1->GetNumberOfClusters();
5029 // Refit Kink under if too small angle
5031 if (kink->GetAngle(2)<0.05){
5032 kink->SetTPCRow0(GetRowNumber(kink->GetR()));
5033 Int_t row0 = kink->GetTPCRow0();
5034 Int_t drow = Int_t(2.+0.5/(0.05+kink->GetAngle(2)));
5037 Int_t last = row0-drow;
5038 if (last<40) last=40;
5039 if (last<ktrack0->GetFirstPoint()+25) last = ktrack0->GetFirstPoint()+25;
5040 AliTPCseed* seed0 = ReSeed(ktrack0,last,kFALSE);
5043 Int_t first = row0+drow;
5044 if (first>130) first=130;
5045 if (first>ktrack1->GetLastPoint()-25) first = TMath::Max(ktrack1->GetLastPoint()-25,30);
5046 AliTPCseed* seed1 = ReSeed(ktrack1,first,kTRUE);
5048 if (seed0 && seed1){
5049 kink->SetStatus(1,8);
5050 if (RefitKink(*seed0,*seed1,*kink)) kink->SetStatus(1,9);
5051 row0 = GetRowNumber(kink->GetR());
5052 sumn = seed0->GetNumberOfClusters()+seed1->GetNumberOfClusters();
5053 mothers[i] = *seed0;
5054 daughters[i] = *seed1;
5057 delete kinks->RemoveAt(i);
5058 if (seed0) delete seed0;
5059 if (seed1) delete seed1;
5062 if (kink->GetDistance()>0.5 || kink->GetR()<110 || kink->GetR()>240) {
5063 delete kinks->RemoveAt(i);
5064 if (seed0) delete seed0;
5065 if (seed1) delete seed1;
5073 if (kink) quality[i] = 160*((0.1+kink->GetDistance())*(2.-kink->GetTPCDensityFactor()))/(sumn+40.); //the longest -clossest will win
5075 TMath::Sort(nkinks,quality,indexes,kFALSE);
5077 //remove double find kinks
5079 for (Int_t ikink0=1;ikink0<nkinks;ikink0++){
5080 AliKink * kink0 = (AliKink*) kinks->At(indexes[ikink0]);
5081 if (!kink0) continue;
5083 for (Int_t ikink1=0;ikink1<ikink0;ikink1++){
5084 if (!kink0) continue;
5085 AliKink * kink1 = (AliKink*) kinks->At(indexes[ikink1]);
5086 if (!kink1) continue;
5087 // if not close kink continue
5088 if (TMath::Abs(kink1->GetPosition()[2]-kink0->GetPosition()[2])>10) continue;
5089 if (TMath::Abs(kink1->GetPosition()[1]-kink0->GetPosition()[1])>10) continue;
5090 if (TMath::Abs(kink1->GetPosition()[0]-kink0->GetPosition()[0])>10) continue;
5092 AliTPCseed &mother0 = mothers[indexes[ikink0]];
5093 AliTPCseed &daughter0 = daughters[indexes[ikink0]];
5094 AliTPCseed &mother1 = mothers[indexes[ikink1]];
5095 AliTPCseed &daughter1 = daughters[indexes[ikink1]];
5096 Int_t row0 = (kink0->GetTPCRow0()+kink1->GetTPCRow0())/2;
5105 for (Int_t i=0;i<row0;i++){
5106 if (mother0.GetClusterIndex(i)>0 && mother1.GetClusterIndex(i)>0){
5109 if (mother0.GetClusterIndex(i)==mother1.GetClusterIndex(i)){
5116 for (Int_t i=row0;i<158;i++){
5117 if (daughter0.GetClusterIndex(i)>0 && daughter0.GetClusterIndex(i)>0){
5120 if (mother0.GetClusterIndex(i)==mother1.GetClusterIndex(i)){
5126 Float_t ratio = Float_t(same+1)/Float_t(both+1);
5127 Float_t ratiom = Float_t(samem+1)/Float_t(bothm+1);
5128 Float_t ratiod = Float_t(samed+1)/Float_t(bothd+1);
5129 if (ratio>0.3 && ratiom>0.5 &&ratiod>0.5) {
5130 Int_t sum0 = mother0.GetNumberOfClusters()+daughter0.GetNumberOfClusters();
5131 Int_t sum1 = mother1.GetNumberOfClusters()+daughter1.GetNumberOfClusters();
5133 shared[kink0->GetIndex(0)]= kTRUE;
5134 shared[kink0->GetIndex(1)]= kTRUE;
5135 delete kinks->RemoveAt(indexes[ikink0]);
5138 shared[kink1->GetIndex(0)]= kTRUE;
5139 shared[kink1->GetIndex(1)]= kTRUE;
5140 delete kinks->RemoveAt(indexes[ikink1]);
5147 for (Int_t i=0;i<nkinks;i++){
5148 AliKink * kink = (AliKink*) kinks->At(indexes[i]);
5149 if (!kink) continue;
5150 kink->SetTPCRow0(GetRowNumber(kink->GetR()));
5151 Int_t index0 = kink->GetIndex(0);
5152 Int_t index1 = kink->GetIndex(1);
5153 if (circular[index0]||circular[index1]&&kink->GetDistance()>0.2) continue;
5154 kink->SetMultiple(usage[index0],0);
5155 kink->SetMultiple(usage[index1],1);
5156 if (kink->GetMultiple()[0]+kink->GetMultiple()[1]>2) continue;
5157 if (kink->GetMultiple()[0]+kink->GetMultiple()[1]>0 && quality[indexes[i]]>0.2) continue;
5158 if (kink->GetMultiple()[0]+kink->GetMultiple()[1]>0 && kink->GetDistance()>0.2) continue;
5159 if (circular[index0]||circular[index1]&&kink->GetDistance()>0.1) continue;
5161 AliTPCseed * ktrack0 = (AliTPCseed*)array->At(index0);
5162 AliTPCseed * ktrack1 = (AliTPCseed*)array->At(index1);
5163 if (!ktrack0 || !ktrack1) continue;
5164 Int_t index = esd->AddKink(kink);
5167 if ( ktrack0->GetKinkIndex(0)==0 && ktrack1->GetKinkIndex(0)==0) { //best kink
5168 if (mothers[indexes[i]].GetNumberOfClusters()>20 && daughters[indexes[i]].GetNumberOfClusters()>20 && (mothers[indexes[i]].GetNumberOfClusters()+daughters[indexes[i]].GetNumberOfClusters())>100){
5169 *ktrack0 = mothers[indexes[i]];
5170 *ktrack1 = daughters[indexes[i]];
5174 ktrack0->SetKinkIndex(usage[index0],-(index+1));
5175 ktrack1->SetKinkIndex(usage[index1], (index+1));
5180 // Remove tracks corresponding to shared kink's
5182 for (Int_t i=0;i<nentries;i++){
5183 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5184 if (!track0) continue;
5185 if (track0->GetKinkIndex(0)!=0) continue;
5186 if (shared[i]) delete array->RemoveAt(i);
5191 RemoveUsed2(array,0.5,0.4,30);
5193 for (Int_t i=0;i<nentries;i++){
5194 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5195 if (!track0) continue;
5196 track0->CookdEdx(0.02,0.6);
5200 for (Int_t i=0;i<nentries;i++){
5201 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5202 if (!track0) continue;
5203 if (track0->Pt()<1.4) continue;
5204 //remove double high momenta tracks - overlapped with kink candidates
5207 for (Int_t icl=track0->GetFirstPoint();icl<track0->GetLastPoint(); icl++){
5208 if (track0->GetClusterPointer(icl)!=0){
5210 if (track0->GetClusterPointer(icl)->IsUsed(10)) shared++;
5213 if (Float_t(shared+1)/Float_t(all+1)>0.5) {
5214 delete array->RemoveAt(i);
5218 if (track0->GetKinkIndex(0)!=0) continue;
5219 if (track0->GetNumberOfClusters()<80) continue;
5221 AliTPCseed *pmother = new AliTPCseed();
5222 AliTPCseed *pdaughter = new AliTPCseed();
5223 AliKink *pkink = new AliKink;
5225 AliTPCseed & mother = *pmother;
5226 AliTPCseed & daughter = *pdaughter;
5227 AliKink & kink = *pkink;
5228 if (CheckKinkPoint(track0,mother,daughter, kink)){
5229 if (mother.GetNumberOfClusters()<30||daughter.GetNumberOfClusters()<20) {
5233 continue; //too short tracks
5235 if (mother.Pt()<1.4) {
5241 Int_t row0= kink.GetTPCRow0();
5242 if (kink.GetDistance()>0.5 || kink.GetR()<110. || kink.GetR()>240.) {
5249 Int_t index = esd->AddKink(&kink);
5250 mother.SetKinkIndex(0,-(index+1));
5251 daughter.SetKinkIndex(0,index+1);
5252 if (mother.GetNumberOfClusters()>50) {
5253 delete array->RemoveAt(i);
5254 array->AddAt(new AliTPCseed(mother),i);
5257 array->AddLast(new AliTPCseed(mother));
5259 array->AddLast(new AliTPCseed(daughter));
5260 for (Int_t icl=0;icl<row0;icl++) {
5261 if (mother.GetClusterPointer(icl)) mother.GetClusterPointer(icl)->Use(20);
5264 for (Int_t icl=row0;icl<158;icl++) {
5265 if (daughter.GetClusterPointer(icl)) daughter.GetClusterPointer(icl)->Use(20);
5274 delete [] daughters;
5296 printf("Ncandidates=\t%d\t%d\t%d\t%d\n",esd->GetNumberOfKinks(),ncandidates,ntracks,nall);
5300 void AliTPCtrackerMI::FindV0s(TObjArray * array, AliESDEvent *esd)
5306 TObjArray *tpcv0s = new TObjArray(100000);
5307 Int_t nentries = array->GetEntriesFast();
5308 AliHelix *helixes = new AliHelix[nentries];
5309 Int_t *sign = new Int_t[nentries];
5310 Float_t *alpha = new Float_t[nentries];
5311 Float_t *z0 = new Float_t[nentries];
5312 Float_t *dca = new Float_t[nentries];
5313 Float_t *sdcar = new Float_t[nentries];
5314 Float_t *cdcar = new Float_t[nentries];
5315 Float_t *pulldcar = new Float_t[nentries];
5316 Float_t *pulldcaz = new Float_t[nentries];
5317 Float_t *pulldca = new Float_t[nentries];
5318 Bool_t *isPrim = new Bool_t[nentries];
5319 const AliESDVertex * primvertex = esd->GetVertex();
5320 Double_t zvertex = primvertex->GetZv();
5322 // nentries = array->GetEntriesFast();
5324 for (Int_t i=0;i<nentries;i++){
5327 AliTPCseed* track = (AliTPCseed*)array->At(i);
5328 if (!track) continue;
5329 track->GetV0Indexes()[0] = 0; //rest v0 indexes
5330 track->GetV0Indexes()[1] = 0; //rest v0 indexes
5331 track->GetV0Indexes()[2] = 0; //rest v0 indexes
5333 alpha[i] = track->GetAlpha();
5334 new (&helixes[i]) AliHelix(*track);
5336 helixes[i].Evaluate(0,xyz);
5337 sign[i] = (track->GetC()>0) ? -1:1;
5341 if (track->GetProlongation(0,y,z)) z0[i] = z;
5342 dca[i] = track->GetD(0,0);
5344 // dca error parrameterezation + pulls
5346 sdcar[i] = TMath::Sqrt(0.150*0.150+(100*track->GetC())*(100*track->GetC()));
5347 if (TMath::Abs(track->GetTgl())>1) sdcar[i]*=2.5;
5348 cdcar[i] = TMath::Exp((TMath::Abs(track->GetC())-0.0106)*525.3);
5349 pulldcar[i] = (dca[i]-cdcar[i])/sdcar[i];
5350 pulldcaz[i] = (z0[i]-zvertex)/sdcar[i];
5351 pulldca[i] = TMath::Sqrt(pulldcar[i]*pulldcar[i]+pulldcaz[i]*pulldcaz[i]);
5352 if (track->TPCrPID(1)+track->TPCrPID(2)+track->TPCrPID(3)>0.5) {
5353 if (pulldca[i]<3.) isPrim[i]=kTRUE; //pion, muon and Kaon 3 sigma cut
5355 if (track->TPCrPID(4)>0.5) {
5356 if (pulldca[i]<0.5) isPrim[i]=kTRUE; //proton 0.5 sigma cut
5358 if (track->TPCrPID(0)>0.4) {
5359 isPrim[i]=kFALSE; //electron no sigma cut
5366 Int_t ncandidates =0;
5369 Double_t phase[2][2],radius[2];
5375 TTreeSRedirector &cstream = *fDebugStreamer;
5376 Float_t fprimvertex[3]={GetX(),GetY(),GetZ()};
5378 Double_t cradius0 = 10*10;
5379 Double_t cradius1 = 200*200;
5382 Double_t cpointAngle = 0.95;
5384 Double_t delta[2]={10000,10000};
5385 for (Int_t i =0;i<nentries;i++){
5386 if (sign[i]==0) continue;
5387 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5388 if (!track0) continue;
5389 if (AliTPCReconstructor::StreamLevel()>1){
5394 "zvertex="<<zvertex<<
5395 "sdcar0="<<sdcar[i]<<
5396 "cdcar0="<<cdcar[i]<<
5397 "pulldcar0="<<pulldcar[i]<<
5398 "pulldcaz0="<<pulldcaz[i]<<
5399 "pulldca0="<<pulldca[i]<<
5400 "isPrim="<<isPrim[i]<<
5404 if (track0->GetSigned1Pt()<0) continue;
5405 if (track0->GetKinkIndex(0)>0||isPrim[i]) continue; //daughter kink
5407 if (TMath::Abs(helixes[i].GetHelix(4))<0.000000001) continue;
5412 for (Int_t j =0;j<nentries;j++){
5413 AliTPCseed * track1 = (AliTPCseed*)array->At(j);
5414 if (!track1) continue;
5415 if (track1->GetKinkIndex(0)>0 || isPrim[j]) continue; //daughter kink
5416 if (sign[j]*sign[i]>0) continue;
5417 if (TMath::Abs(helixes[j].GetHelix(4))<0.000001) continue;
5418 if (track0->GetCircular()+track1->GetCircular()>1) continue; //circling -returning track
5421 // DCA to prim vertex cut
5427 Int_t npoints = helixes[i].GetRPHIintersections(helixes[j], phase, radius,cdist2);
5428 if (npoints<1) continue;
5432 if (radius[0]<cradius0||radius[0]>cradius1) continue;
5433 helixes[i].LinearDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta[0]);
5434 if (delta[0]>cdist1) continue;
5437 if (TMath::Max(radius[0],radius[1])<cradius0|| TMath::Min(radius[0],radius[1])>cradius1) continue;
5438 helixes[i].LinearDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta[0]);
5439 helixes[i].LinearDCA(helixes[j],phase[1][0],phase[1][1],radius[1],delta[1]);
5440 if (delta[1]<delta[0]) iclosest=1;
5441 if (delta[iclosest]>cdist1) continue;
5443 helixes[i].ParabolicDCA(helixes[j],phase[iclosest][0],phase[iclosest][1],radius[iclosest],delta[iclosest]);
5444 if (radius[iclosest]<cradius0 || radius[iclosest]>cradius1 || delta[iclosest]>cdist1) continue;
5446 Double_t pointAngle = helixes[i].GetPointAngle(helixes[j],phase[iclosest],fprimvertex);
5447 if (pointAngle<cpointAngle) continue;
5449 Bool_t isGamma = kFALSE;
5450 vertex.SetParamP(*track0); //track0 - plus
5451 vertex.SetParamN(*track1); //track1 - minus
5452 vertex.Update(fprimvertex);
5453 if (track0->TPCrPID(0)>0.3&&track1->TPCrPID(0)>0.3&&vertex.GetAnglep()[2]<0.15) isGamma=kTRUE; // gamma conversion candidate
5454 Double_t pointAngle2 = vertex.GetV0CosineOfPointingAngle();
5456 if (vertex.GetV0CosineOfPointingAngle()<cpointAngle && (!isGamma)) continue;// point angle cut
5457 //Bo: if (vertex.GetDist2()>2&&(!isGamma)) continue; // point angle cut
5458 if (vertex.GetDcaV0Daughters()>2&&(!isGamma)) continue;//Bo: // point angle cut
5459 Float_t sigmae = 0.15*0.15;
5460 if (vertex.GetRr()<80)
5461 sigmae += (sdcar[i]*sdcar[i]+sdcar[j]*sdcar[j])*(1.-vertex.GetRr()/80.)*(1.-vertex.GetRr()/80.);
5462 sigmae+= TMath::Sqrt(sigmae);
5463 //Bo: if (vertex.GetDist2()/sigmae>3.&&(!isGamma)) continue;
5464 if (vertex.GetDcaV0Daughters()/sigmae>3.&&(!isGamma)) continue;
5465 Float_t densb0=0,densb1=0,densa0=0,densa1=0;
5466 Int_t row0 = GetRowNumber(vertex.GetRr());
5468 //Bo: if (vertex.GetDist2()>0.2) continue;
5469 if (vertex.GetDcaV0Daughters()>0.2) continue;
5470 densb0 = track0->Density2(0,row0-5);
5471 densb1 = track1->Density2(0,row0-5);
5472 if (densb0>0.3|| densb1>0.3) continue; //clusters before vertex
5473 densa0 = track0->Density2(row0+5,row0+40);
5474 densa1 = track1->Density2(row0+5,row0+40);
5475 if ((densa0<0.4|| densa1<0.4)&&(!isGamma)) continue; //missing clusters after vertex
5478 densa0 = track0->Density2(0,40); //cluster density
5479 densa1 = track1->Density2(0,40); //cluster density
5480 if ((vertex.GetRr()<80&&densa0+densa1<1.)&&(!isGamma)) continue;
5482 //Bo: vertex.SetLab(0,track0->GetLabel());
5483 //Bo: vertex.SetLab(1,track1->GetLabel());
5484 vertex.SetChi2After((densa0+densa1)*0.5);
5485 vertex.SetChi2Before((densb0+densb1)*0.5);
5486 vertex.SetIndex(0,i);
5487 vertex.SetIndex(1,j);
5488 //Bo: vertex.SetStatus(1); // TPC v0 candidate
5489 vertex.SetOnFlyStatus(2);//Bo: // TPC v0 candidate
5490 //Bo: vertex.SetRp(track0->TPCrPIDs());
5491 //Bo: vertex.SetRm(track1->TPCrPIDs());
5492 tpcv0s->AddLast(new AliESDv0(vertex));
5495 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
5496 Double_t radiusm= (delta[0]<delta[1])? TMath::Sqrt(radius[0]):TMath::Sqrt(radius[1]);
5497 Double_t deltam= (delta[0]<delta[1])? TMath::Sqrt(delta[0]):TMath::Sqrt(delta[1]);
5498 if (AliTPCReconstructor::StreamLevel()>1) {
5499 Int_t lab0=track0->GetLabel();
5500 Int_t lab1=track1->GetLabel();
5501 Char_t c0=track0->GetCircular();
5502 Char_t c1=track1->GetCircular();
5505 "vertex.="<<&vertex<<
5508 "Helix0.="<<&helixes[i]<<
5511 "Helix1.="<<&helixes[j]<<
5512 "pointAngle="<<pointAngle<<
5513 "pointAngle2="<<pointAngle2<<
5518 "zvertex="<<zvertex<<
5521 "npoints="<<npoints<<
5522 "radius0="<<radius[0]<<
5523 "delta0="<<delta[0]<<
5524 "radius1="<<radius[1]<<
5525 "delta1="<<delta[1]<<
5526 "radiusm="<<radiusm<<
5528 "sdcar0="<<sdcar[i]<<
5529 "sdcar1="<<sdcar[j]<<
5530 "cdcar0="<<cdcar[i]<<
5531 "cdcar1="<<cdcar[j]<<
5532 "pulldcar0="<<pulldcar[i]<<
5533 "pulldcar1="<<pulldcar[j]<<
5534 "pulldcaz0="<<pulldcaz[i]<<
5535 "pulldcaz1="<<pulldcaz[j]<<
5536 "pulldca0="<<pulldca[i]<<
5537 "pulldca1="<<pulldca[j]<<
5547 Float_t *quality = new Float_t[ncandidates];
5548 Int_t *indexes = new Int_t[ncandidates];
5550 for (Int_t i=0;i<ncandidates;i++){
5552 AliESDv0 *v0 = (AliESDv0*)tpcv0s->At(i);
5553 quality[i] = 1./(1.00001-v0->GetV0CosineOfPointingAngle()); //base point angle
5554 // quality[i] /= (0.5+v0->GetDist2());
5555 // quality[i] *= v0->GetChi2After(); //density factor
5557 Int_t index0 = v0->GetIndex(0);
5558 Int_t index1 = v0->GetIndex(1);
5559 //Bo: Double_t minpulldca = TMath::Min(2.+pulldca[v0->GetIndex(0)],(2.+pulldca[v0->GetIndex(1)]) ); //pull
5560 Double_t minpulldca = TMath::Min(2.+pulldca[index0],(2.+pulldca[index1]) );//Bo:
5564 AliTPCseed * track0 = (AliTPCseed*)array->At(index0);
5565 AliTPCseed * track1 = (AliTPCseed*)array->At(index1);
5566 if (track0->TPCrPID(0)>0.3&&track1->TPCrPID(0)>0.3&&v0->GetAnglep()[2]<0.15) quality[i]+=1000000; // gamma conversion candidate
5567 if (track0->TPCrPID(4)>0.9||track1->TPCrPID(4)>0.9&&minpulldca>4) quality[i]*=10; // lambda candidate candidate
5570 TMath::Sort(ncandidates,quality,indexes,kTRUE);
5573 for (Int_t i=0;i<ncandidates;i++){
5574 AliESDv0 * v0 = (AliESDv0*)tpcv0s->At(indexes[i]);
5576 Int_t index0 = v0->GetIndex(0);
5577 Int_t index1 = v0->GetIndex(1);
5578 AliTPCseed * track0 = (AliTPCseed*)array->At(index0);
5579 AliTPCseed * track1 = (AliTPCseed*)array->At(index1);
5580 if (!track0||!track1) {
5584 Bool_t accept =kTRUE; //default accept
5585 Int_t *v0indexes0 = track0->GetV0Indexes();
5586 Int_t *v0indexes1 = track1->GetV0Indexes();
5588 Int_t order0 = (v0indexes0[0]!=0) ? 1:0;
5589 Int_t order1 = (v0indexes1[0]!=0) ? 1:0;
5590 if (v0indexes0[1]!=0) order0 =2;
5591 if (v0indexes1[1]!=0) order1 =2;
5593 if (v0indexes0[2]!=0) {order0=3; accept=kFALSE;}
5594 if (v0indexes0[2]!=0) {order1=3; accept=kFALSE;}
5596 AliESDv0 * v02 = v0;
5598 //Bo: v0->SetOrder(0,order0);
5599 //Bo: v0->SetOrder(1,order1);
5600 //Bo: v0->SetOrder(1,order0+order1);
5601 v0->SetOnFlyStatus(kTRUE);
5602 Int_t index = esd->AddV0(v0);
5603 v02 = esd->GetV0(index);
5604 v0indexes0[order0]=index;
5605 v0indexes1[order1]=index;
5609 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
5610 if (AliTPCReconstructor::StreamLevel()>1) {
5611 Int_t lab0=track0->GetLabel();
5612 Int_t lab1=track1->GetLabel();
5621 "dca0="<<dca[index0]<<
5622 "dca1="<<dca[index1]<<
5626 "quality="<<quality[i]<<
5627 "pulldca0="<<pulldca[index0]<<
5628 "pulldca1="<<pulldca[index1]<<
5652 printf("TPC V0 finder : naccepted\t%d\tncandidates\t%d\tntracks\t%d\tnall\t%d\n",naccepted,ncandidates,ntracks,nall);
5656 Int_t AliTPCtrackerMI::RefitKink(AliTPCseed &mother, AliTPCseed &daughter, AliESDkink &knk)
5659 // refit kink towards to the vertex
5662 AliKink &kink=(AliKink &)knk;
5664 Int_t row0 = GetRowNumber(kink.GetR());
5665 FollowProlongation(mother,0);
5666 mother.Reset(kFALSE);
5668 FollowProlongation(daughter,row0);
5669 daughter.Reset(kFALSE);
5670 FollowBackProlongation(daughter,158);
5671 daughter.Reset(kFALSE);
5672 Int_t first = TMath::Max(row0-20,30);
5673 Int_t last = TMath::Min(row0+20,140);
5675 const Int_t kNdiv =5;
5676 AliTPCseed param0[kNdiv]; // parameters along the track
5677 AliTPCseed param1[kNdiv]; // parameters along the track
5678 AliKink kinks[kNdiv]; // corresponding kink parameters
5681 for (Int_t irow=0; irow<kNdiv;irow++){
5682 rows[irow] = first +((last-first)*irow)/(kNdiv-1);
5684 // store parameters along the track
5686 for (Int_t irow=0;irow<kNdiv;irow++){
5687 FollowBackProlongation(mother, rows[irow]);
5688 FollowProlongation(daughter,rows[kNdiv-1-irow]);
5689 param0[irow] = mother;
5690 param1[kNdiv-1-irow] = daughter;
5694 for (Int_t irow=0; irow<kNdiv-1;irow++){
5695 if (param0[irow].GetNumberOfClusters()<kNdiv||param1[irow].GetNumberOfClusters()<kNdiv) continue;
5696 kinks[irow].SetMother(param0[irow]);
5697 kinks[irow].SetDaughter(param1[irow]);
5698 kinks[irow].Update();
5701 // choose kink with best "quality"
5703 Double_t mindist = 10000;
5704 for (Int_t irow=0;irow<kNdiv;irow++){
5705 if (param0[irow].GetNumberOfClusters()<20||param1[irow].GetNumberOfClusters()<20) continue;
5706 if (TMath::Abs(kinks[irow].GetR())>240.) continue;
5707 if (TMath::Abs(kinks[irow].GetR())<100.) continue;
5709 Float_t normdist = TMath::Abs(param0[irow].GetX()-kinks[irow].GetR())*(0.1+kink.GetDistance());
5710 normdist/= (param0[irow].GetNumberOfClusters()+param1[irow].GetNumberOfClusters()+40.);
5711 if (normdist < mindist){
5717 if (index==-1) return 0;
5720 param0[index].Reset(kTRUE);
5721 FollowProlongation(param0[index],0);
5723 mother = param0[index];
5724 daughter = param1[index]; // daughter in vertex
5726 kink.SetMother(mother);
5727 kink.SetDaughter(daughter);
5729 kink.SetTPCRow0(GetRowNumber(kink.GetR()));
5730 kink.SetTPCncls(param0[index].GetNumberOfClusters(),0);
5731 kink.SetTPCncls(param1[index].GetNumberOfClusters(),1);
5732 kink.SetLabel(CookLabel(&mother,0.4, 0,kink.GetTPCRow0()),0);
5733 kink.SetLabel(CookLabel(&daughter,0.4, kink.GetTPCRow0(),160),1);
5734 mother.SetLabel(kink.GetLabel(0));
5735 daughter.SetLabel(kink.GetLabel(1));
5741 void AliTPCtrackerMI::UpdateKinkQualityM(AliTPCseed * seed){
5743 // update Kink quality information for mother after back propagation
5745 if (seed->GetKinkIndex(0)>=0) return;
5746 for (Int_t ikink=0;ikink<3;ikink++){
5747 Int_t index = seed->GetKinkIndex(ikink);
5748 if (index>=0) break;
5749 index = TMath::Abs(index)-1;
5750 AliESDkink * kink = fEvent->GetKink(index);
5751 kink->SetTPCDensity(-1,0,0);
5752 kink->SetTPCDensity(1,0,1);
5754 Int_t row0 = kink->GetTPCRow0() - 2 - Int_t( 0.5/ (0.05+kink->GetAngle(2)));
5755 if (row0<15) row0=15;
5757 Int_t row1 = kink->GetTPCRow0() + 2 + Int_t( 0.5/ (0.05+kink->GetAngle(2)));
5758 if (row1>145) row1=145;
5760 Int_t found,foundable,shared;
5761 seed->GetClusterStatistic(0,row0, found, foundable,shared,kFALSE);
5762 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),0,0);
5763 seed->GetClusterStatistic(row1,155, found, foundable,shared,kFALSE);
5764 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),0,1);
5769 void AliTPCtrackerMI::UpdateKinkQualityD(AliTPCseed * seed){
5771 // update Kink quality information for daughter after refit
5773 if (seed->GetKinkIndex(0)<=0) return;
5774 for (Int_t ikink=0;ikink<3;ikink++){
5775 Int_t index = seed->GetKinkIndex(ikink);
5776 if (index<=0) break;
5777 index = TMath::Abs(index)-1;
5778 AliESDkink * kink = fEvent->GetKink(index);
5779 kink->SetTPCDensity(-1,1,0);
5780 kink->SetTPCDensity(-1,1,1);
5782 Int_t row0 = kink->GetTPCRow0() -2 - Int_t( 0.5/ (0.05+kink->GetAngle(2)));
5783 if (row0<15) row0=15;
5785 Int_t row1 = kink->GetTPCRow0() +2 + Int_t( 0.5/ (0.05+kink->GetAngle(2)));
5786 if (row1>145) row1=145;
5788 Int_t found,foundable,shared;
5789 seed->GetClusterStatistic(0,row0, found, foundable,shared,kFALSE);
5790 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),1,0);
5791 seed->GetClusterStatistic(row1,155, found, foundable,shared,kFALSE);
5792 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),1,1);
5798 Int_t AliTPCtrackerMI::CheckKinkPoint(AliTPCseed*seed,AliTPCseed &mother, AliTPCseed &daughter, AliESDkink &knk)
5801 // check kink point for given track
5802 // if return value=0 kink point not found
5803 // otherwise seed0 correspond to mother particle
5804 // seed1 correspond to daughter particle
5805 // kink parameter of kink point
5806 AliKink &kink=(AliKink &)knk;
5808 Int_t middlerow = (seed->GetFirstPoint()+seed->GetLastPoint())/2;
5809 Int_t first = seed->GetFirstPoint();
5810 Int_t last = seed->GetLastPoint();
5811 if (last-first<20) return 0; // shortest length - 2*30 = 60 pad-rows
5814 AliTPCseed *seed1 = ReSeed(seed,middlerow+20, kTRUE); //middle of chamber
5815 if (!seed1) return 0;
5816 FollowProlongation(*seed1,seed->GetLastPoint()-20);
5817 seed1->Reset(kTRUE);
5818 FollowProlongation(*seed1,158);
5819 seed1->Reset(kTRUE);
5820 last = seed1->GetLastPoint();
5822 AliTPCseed *seed0 = new AliTPCseed(*seed);
5823 seed0->Reset(kFALSE);
5826 AliTPCseed param0[20]; // parameters along the track
5827 AliTPCseed param1[20]; // parameters along the track
5828 AliKink kinks[20]; // corresponding kink parameters
5830 for (Int_t irow=0; irow<20;irow++){
5831 rows[irow] = first +((last-first)*irow)/19;
5833 // store parameters along the track
5835 for (Int_t irow=0;irow<20;irow++){
5836 FollowBackProlongation(*seed0, rows[irow]);
5837 FollowProlongation(*seed1,rows[19-irow]);
5838 param0[irow] = *seed0;
5839 param1[19-irow] = *seed1;
5843 for (Int_t irow=0; irow<19;irow++){
5844 kinks[irow].SetMother(param0[irow]);
5845 kinks[irow].SetDaughter(param1[irow]);
5846 kinks[irow].Update();
5849 // choose kink with biggest change of angle
5851 Double_t maxchange= 0;
5852 for (Int_t irow=1;irow<19;irow++){
5853 if (TMath::Abs(kinks[irow].GetR())>240.) continue;
5854 if (TMath::Abs(kinks[irow].GetR())<110.) continue;
5855 Float_t quality = TMath::Abs(kinks[irow].GetAngle(2))/(3.+TMath::Abs(kinks[irow].GetR()-param0[irow].GetX()));
5856 if ( quality > maxchange){
5857 maxchange = quality;
5864 if (index<0) return 0;
5866 Int_t row0 = GetRowNumber(kinks[index].GetR()); //row 0 estimate
5867 seed0 = new AliTPCseed(param0[index]);
5868 seed1 = new AliTPCseed(param1[index]);
5869 seed0->Reset(kFALSE);
5870 seed1->Reset(kFALSE);
5871 seed0->ResetCovariance(10.);
5872 seed1->ResetCovariance(10.);
5873 FollowProlongation(*seed0,0);
5874 FollowBackProlongation(*seed1,158);
5875 mother = *seed0; // backup mother at position 0
5876 seed0->Reset(kFALSE);
5877 seed1->Reset(kFALSE);
5878 seed0->ResetCovariance(10.);
5879 seed1->ResetCovariance(10.);
5881 first = TMath::Max(row0-20,0);
5882 last = TMath::Min(row0+20,158);
5884 for (Int_t irow=0; irow<20;irow++){
5885 rows[irow] = first +((last-first)*irow)/19;
5887 // store parameters along the track
5889 for (Int_t irow=0;irow<20;irow++){
5890 FollowBackProlongation(*seed0, rows[irow]);
5891 FollowProlongation(*seed1,rows[19-irow]);
5892 param0[irow] = *seed0;
5893 param1[19-irow] = *seed1;
5897 for (Int_t irow=0; irow<19;irow++){
5898 kinks[irow].SetMother(param0[irow]);
5899 kinks[irow].SetDaughter(param1[irow]);
5900 // param0[irow].Dump();
5901 //param1[irow].Dump();
5902 kinks[irow].Update();
5905 // choose kink with biggest change of angle
5908 for (Int_t irow=0;irow<20;irow++){
5909 if (TMath::Abs(kinks[irow].GetR())>250.) continue;
5910 if (TMath::Abs(kinks[irow].GetR())<90.) continue;
5911 Float_t quality = TMath::Abs(kinks[irow].GetAngle(2))/(3.+TMath::Abs(kinks[irow].GetR()-param0[irow].GetX()));
5912 if ( quality > maxchange){
5913 maxchange = quality;
5920 if (index==-1 || param0[index].GetNumberOfClusters()+param1[index].GetNumberOfClusters()<100){
5925 // Float_t anglesigma = TMath::Sqrt(param0[index].fC22+param0[index].fC33+param1[index].fC22+param1[index].fC33);
5927 kink.SetMother(param0[index]);
5928 kink.SetDaughter(param1[index]);
5930 row0 = GetRowNumber(kink.GetR());
5931 kink.SetTPCRow0(row0);
5932 kink.SetLabel(CookLabel(seed0,0.5,0,row0),0);
5933 kink.SetLabel(CookLabel(seed1,0.5,row0,158),1);
5934 kink.SetIndex(-10,0);
5935 kink.SetIndex(int(param0[index].GetNumberOfClusters()+param1[index].GetNumberOfClusters()),1);
5936 kink.SetTPCncls(param0[index].GetNumberOfClusters(),0);
5937 kink.SetTPCncls(param1[index].GetNumberOfClusters(),1);
5940 // new (&mother) AliTPCseed(param0[index]);
5941 daughter = param1[index];
5942 daughter.SetLabel(kink.GetLabel(1));
5943 param0[index].Reset(kTRUE);
5944 FollowProlongation(param0[index],0);
5945 mother = param0[index];
5946 mother.SetLabel(kink.GetLabel(0));
5956 AliTPCseed* AliTPCtrackerMI::ReSeed(AliTPCseed *t)
5959 // reseed - refit - track
5962 // Int_t last = fSectors->GetNRows()-1;
5964 if (fSectors == fOuterSec){
5965 first = TMath::Max(first, t->GetFirstPoint()-fInnerSec->GetNRows());
5969 first = t->GetFirstPoint();
5971 AliTPCseed * seed = MakeSeed(t,0.1,0.5,0.9);
5972 FollowBackProlongation(*t,fSectors->GetNRows()-1);
5974 FollowProlongation(*t,first);
5984 //_____________________________________________________________________________
5985 Int_t AliTPCtrackerMI::ReadSeeds(const TFile *inp) {
5986 //-----------------------------------------------------------------
5987 // This function reades track seeds.
5988 //-----------------------------------------------------------------
5989 TDirectory *savedir=gDirectory;
5991 TFile *in=(TFile*)inp;
5992 if (!in->IsOpen()) {
5993 cerr<<"AliTPCtrackerMI::ReadSeeds(): input file is not open !\n";
5998 TTree *seedTree=(TTree*)in->Get("Seeds");
6000 cerr<<"AliTPCtrackerMI::ReadSeeds(): ";
6001 cerr<<"can't get a tree with track seeds !\n";
6004 AliTPCtrack *seed=new AliTPCtrack;
6005 seedTree->SetBranchAddress("tracks",&seed);
6007 if (fSeeds==0) fSeeds=new TObjArray(15000);
6009 Int_t n=(Int_t)seedTree->GetEntries();
6010 for (Int_t i=0; i<n; i++) {
6011 seedTree->GetEvent(i);
6012 fSeeds->AddLast(new AliTPCseed(*seed/*,seed->GetAlpha()*/));
6021 Int_t AliTPCtrackerMI::Clusters2Tracks (AliESDEvent *esd)
6024 if (fSeeds) DeleteSeeds();
6027 if (!fSeeds) return 1;
6034 //_____________________________________________________________________________
6035 Int_t AliTPCtrackerMI::Clusters2Tracks() {
6036 //-----------------------------------------------------------------
6037 // This is a track finder.
6038 //-----------------------------------------------------------------
6039 TDirectory *savedir=gDirectory;
6043 fSeeds = Tracking();
6046 Info("Clusters2Tracks","Time for tracking: \t");timer.Print();timer.Start();
6048 //activate again some tracks
6049 for (Int_t i=0; i<fSeeds->GetEntriesFast(); i++) {
6050 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
6052 Int_t nc=t.GetNumberOfClusters();
6054 delete fSeeds->RemoveAt(i);
6058 if (pt->GetRemoval()==10) {
6059 if (pt->GetDensityFirst(20)>0.8 || pt->GetDensityFirst(30)>0.8 || pt->GetDensityFirst(40)>0.7)
6060 pt->Desactivate(10); // make track again active
6062 pt->Desactivate(20);
6063 delete fSeeds->RemoveAt(i);
6068 RemoveUsed2(fSeeds,0.85,0.85,0);
6069 if (AliTPCReconstructor::GetRecoParam()->GetDoKinks()) FindKinks(fSeeds,fEvent);
6070 //FindCurling(fSeeds, fEvent,0);
6071 if (AliTPCReconstructor::StreamLevel()>2) FindMultiMC(fSeeds, fEvent,0); // find multi found tracks
6072 RemoveUsed2(fSeeds,0.5,0.4,20);
6073 FindSplitted(fSeeds, fEvent,0); // find multi found tracks
6076 // // refit short tracks
6078 Int_t nseed=fSeeds->GetEntriesFast();
6081 for (Int_t i=0; i<nseed; i++) {
6082 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
6084 Int_t nc=t.GetNumberOfClusters();
6086 delete fSeeds->RemoveAt(i);
6089 CookLabel(pt,0.1); //For comparison only
6090 //if ((pt->IsActive() || (pt->fRemoval==10) )&& nc>50 &&pt->GetNumberOfClusters()>0.4*pt->fNFoundable){
6091 if ((pt->IsActive() || (pt->GetRemoval()==10) )){
6093 if (fDebug>0) cerr<<found<<'\r';
6097 delete fSeeds->RemoveAt(i);
6101 //RemoveOverlap(fSeeds,0.99,7,kTRUE);
6103 //RemoveUsed(fSeeds,0.9,0.9,6);
6105 nseed=fSeeds->GetEntriesFast();
6107 for (Int_t i=0; i<nseed; i++) {
6108 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
6110 Int_t nc=t.GetNumberOfClusters();
6112 delete fSeeds->RemoveAt(i);
6116 t.CookdEdx(0.02,0.6);
6117 // CheckKinkPoint(&t,0.05);
6118 //if ((pt->IsActive() || (pt->fRemoval==10) )&& nc>50 &&pt->GetNumberOfClusters()>0.4*pt->fNFoundable){
6119 if ((pt->IsActive() || (pt->GetRemoval()==10) )){
6127 delete fSeeds->RemoveAt(i);
6128 //AliTPCseed * seed1 = ReSeed(pt,0.05,0.5,1);
6130 // FollowProlongation(*seed1,0);
6131 // Int_t n = seed1->GetNumberOfClusters();
6132 // printf("fP4\t%f\t%f\n",seed1->GetC(),pt->GetC());
6133 // printf("fN\t%d\t%d\n", seed1->GetNumberOfClusters(),pt->GetNumberOfClusters());
6136 //AliTPCseed * seed2 = ReSeed(pt,0.95,0.5,0.05);
6140 SortTracks(fSeeds, 1);
6144 PrepareForBackProlongation(fSeeds,5.);
6145 PropagateBack(fSeeds);
6146 printf("Time for back propagation: \t");timer.Print();timer.Start();
6150 PrepareForProlongation(fSeeds,5.);
6151 PropagateForward2(fSeeds);
6153 printf("Time for FORWARD propagation: \t");timer.Print();timer.Start();
6154 // RemoveUsed(fSeeds,0.7,0.7,6);
6155 //RemoveOverlap(fSeeds,0.9,7,kTRUE);
6157 nseed=fSeeds->GetEntriesFast();
6159 for (Int_t i=0; i<nseed; i++) {
6160 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
6162 Int_t nc=t.GetNumberOfClusters();
6164 delete fSeeds->RemoveAt(i);
6167 t.CookdEdx(0.02,0.6);
6168 // CookLabel(pt,0.1); //For comparison only
6169 //if ((pt->IsActive() || (pt->fRemoval==10) )&& nc>50 &&pt->GetNumberOfClusters()>0.4*pt->fNFoundable){
6170 if ((pt->IsActive() || (pt->fRemoval==10) )){
6171 cerr<<found++<<'\r';
6174 delete fSeeds->RemoveAt(i);
6179 // fNTracks = found;
6181 Info("Clusters2Tracks","Time for overlap removal, track writing and dedx cooking: \t"); timer.Print();timer.Start();
6184 // cerr<<"Number of found tracks : "<<"\t"<<found<<endl;
6185 Info("Clusters2Tracks","Number of found tracks %d",found);
6187 // UnloadClusters();
6192 void AliTPCtrackerMI::Tracking(TObjArray * arr)
6195 // tracking of the seeds
6198 fSectors = fOuterSec;
6199 ParallelTracking(arr,150,63);
6200 fSectors = fOuterSec;
6201 ParallelTracking(arr,63,0);
6204 TObjArray * AliTPCtrackerMI::Tracking(Int_t seedtype, Int_t i1, Int_t i2, Float_t cuts[4], Float_t dy, Int_t dsec)
6209 TObjArray * arr = new TObjArray;
6211 fSectors = fOuterSec;
6214 for (Int_t sec=0;sec<fkNOS;sec++){
6215 if (seedtype==3) MakeSeeds3(arr,sec,i1,i2,cuts,dy, dsec);
6216 if (seedtype==4) MakeSeeds5(arr,sec,i1,i2,cuts,dy);
6217 if (seedtype==2) MakeSeeds2(arr,sec,i1,i2,cuts,dy);
6220 Info("Tracking","\nSeeding - %d\t%d\t%d\t%d\n",seedtype,i1,i2,arr->GetEntriesFast());
6232 TObjArray * AliTPCtrackerMI::Tracking()
6236 if (AliTPCReconstructor::GetRecoParam()->GetSpecialSeeding()) return TrackingSpecial();
6239 Int_t nup=fOuterSec->GetNRows()+fInnerSec->GetNRows();
6241 TObjArray * seeds = new TObjArray;
6250 Float_t fnumber = 3.0;
6251 Float_t fdensity = 3.0;
6256 for (Int_t delta = 0; delta<18; delta+=6){
6260 arr = Tracking(3,nup-1-delta,nup-1-delta-gap,cuts,-1,1);
6261 SumTracks(seeds,arr);
6262 SignClusters(seeds,fnumber,fdensity);
6264 for (Int_t i=2;i<6;i+=2){
6265 // seed high pt tracks
6268 arr = Tracking(3,nup-i-delta,nup-i-delta-gap,cuts,-1,0);
6269 SumTracks(seeds,arr);
6270 SignClusters(seeds,fnumber,fdensity);
6275 // RemoveUsed(seeds,0.9,0.9,1);
6276 // UnsignClusters();
6277 // SignClusters(seeds,fnumber,fdensity);
6281 for (Int_t delta = 20; delta<120; delta+=10){
6283 // seed high pt tracks
6287 arr = Tracking(3,nup-delta,nup-delta-gap,cuts,-1);
6288 SumTracks(seeds,arr);
6289 SignClusters(seeds,fnumber,fdensity);
6294 arr = Tracking(3,nup-delta-5,nup-delta-5-gap,cuts,-1);
6295 SumTracks(seeds,arr);
6296 SignClusters(seeds,fnumber,fdensity);
6307 Info("Tracking()","\n\nPrimary seeding\t%d\n\n",seeds->GetEntriesFast());
6311 // RemoveUsed(seeds,0.75,0.75,1);
6313 //SignClusters(seeds,fnumber,fdensity);
6322 arr = Tracking(4,nup-1,nup-1-gap,cuts,-1);
6323 SumTracks(seeds,arr);
6324 SignClusters(seeds,fnumber,fdensity);
6326 arr = Tracking(4,nup-2,nup-2-gap,cuts,-1);
6327 SumTracks(seeds,arr);
6328 SignClusters(seeds,fnumber,fdensity);
6330 arr = Tracking(4,nup-3,nup-3-gap,cuts,-1);
6331 SumTracks(seeds,arr);
6332 SignClusters(seeds,fnumber,fdensity);
6336 for (Int_t delta = 3; delta<30; delta+=5){
6342 arr = Tracking(4,nup-1-delta,nup-1-delta-gap,cuts,-1);
6343 SumTracks(seeds,arr);
6344 SignClusters(seeds,fnumber,fdensity);
6346 arr = Tracking(4,nup-3-delta,nup-5-delta-gap,cuts,4);
6347 SumTracks(seeds,arr);
6348 SignClusters(seeds,fnumber,fdensity);
6359 Int_t fLastSeedRowSec=AliTPCReconstructor::GetRecoParam()->GetLastSeedRowSec();
6362 for (Int_t delta = 30; delta<fLastSeedRowSec; delta+=10){
6368 arr = Tracking(4,nup-1-delta,nup-1-delta-gap,cuts,-1);
6369 SumTracks(seeds,arr);
6370 SignClusters(seeds,fnumber,fdensity);
6372 arr = Tracking(4,nup-5-delta,nup-5-delta-gap,cuts,5 );
6373 SumTracks(seeds,arr);
6374 SignClusters(seeds,fnumber,fdensity);
6378 Info("Tracking()","\n\nSecondary seeding\t%d\n\n",seeds->GetEntriesFast());
6389 TObjArray * AliTPCtrackerMI::TrackingSpecial()
6392 // seeding adjusted for laser and cosmic tests - short tracks with big inclination angle
6393 // no primary vertex seeding tried
6397 Int_t nup=fOuterSec->GetNRows()+fInnerSec->GetNRows();
6399 TObjArray * seeds = new TObjArray;
6404 Float_t fnumber = 3.0;
6405 Float_t fdensity = 3.0;
6408 cuts[0] = AliTPCReconstructor::GetRecoParam()->GetMaxC(); // max curvature
6409 cuts[1] = 3.5; // max tan(phi) angle for seeding
6410 cuts[2] = 3.; // not used (cut on z primary vertex)
6411 cuts[3] = 3.5; // max tan(theta) angle for seeding
6413 for (Int_t delta = 0; nup-delta-gap-1>0; delta+=3){
6415 arr = Tracking(4,nup-1-delta,nup-1-delta-gap,cuts,-1);
6416 SumTracks(seeds,arr);
6417 SignClusters(seeds,fnumber,fdensity);
6421 Info("Tracking()","\n\nSecondary seeding\t%d\n\n",seeds->GetEntriesFast());
6432 void AliTPCtrackerMI::SumTracks(TObjArray *arr1,TObjArray *&arr2) const
6435 //sum tracks to common container
6436 //remove suspicious tracks
6437 Int_t nseed = arr2->GetEntriesFast();
6438 for (Int_t i=0;i<nseed;i++){
6439 AliTPCseed *pt=(AliTPCseed*)arr2->UncheckedAt(i);
6442 // remove tracks with too big curvature
6444 if (TMath::Abs(pt->GetC())>AliTPCReconstructor::GetRecoParam()->GetMaxC()){
6445 delete arr2->RemoveAt(i);
6448 // REMOVE VERY SHORT TRACKS
6449 if (pt->GetNumberOfClusters()<20){
6450 delete arr2->RemoveAt(i);
6453 // NORMAL ACTIVE TRACK
6454 if (pt->IsActive()){
6455 arr1->AddLast(arr2->RemoveAt(i));
6458 //remove not usable tracks
6459 if (pt->GetRemoval()!=10){
6460 delete arr2->RemoveAt(i);
6464 // ENABLE ONLY ENOUGH GOOD STOPPED TRACKS
6465 if (pt->GetDensityFirst(20)>0.8 || pt->GetDensityFirst(30)>0.8 || pt->GetDensityFirst(40)>0.7)
6466 arr1->AddLast(arr2->RemoveAt(i));
6468 delete arr2->RemoveAt(i);
6472 delete arr2; arr2 = 0;
6477 void AliTPCtrackerMI::ParallelTracking(TObjArray * arr, Int_t rfirst, Int_t rlast)
6480 // try to track in parralel
6482 Int_t nseed=arr->GetEntriesFast();
6483 //prepare seeds for tracking
6484 for (Int_t i=0; i<nseed; i++) {
6485 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i), &t=*pt;
6487 if (!t.IsActive()) continue;
6488 // follow prolongation to the first layer
6489 if ( (fSectors ==fInnerSec) || (t.GetFirstPoint()-fParam->GetNRowLow()>rfirst+1) )
6490 FollowProlongation(t, rfirst+1);
6495 for (Int_t nr=rfirst; nr>=rlast; nr--){
6496 if (nr<fInnerSec->GetNRows())
6497 fSectors = fInnerSec;
6499 fSectors = fOuterSec;
6500 // make indexes with the cluster tracks for given
6502 // find nearest cluster
6503 for (Int_t i=0; i<nseed; i++) {
6504 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i), &t=*pt;
6506 if (nr==80) pt->UpdateReference();
6507 if (!pt->IsActive()) continue;
6508 // if ( (fSectors ==fOuterSec) && (pt->fFirstPoint-fParam->GetNRowLow())<nr) continue;
6509 if (pt->GetRelativeSector()>17) {
6512 UpdateClusters(t,nr);
6514 // prolonagate to the nearest cluster - if founded
6515 for (Int_t i=0; i<nseed; i++) {
6516 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
6518 if (!pt->IsActive()) continue;
6519 // if ((fSectors ==fOuterSec) && (pt->fFirstPoint-fParam->GetNRowLow())<nr) continue;
6520 if (pt->GetRelativeSector()>17) {
6523 FollowToNextCluster(*pt,nr);
6528 void AliTPCtrackerMI::PrepareForBackProlongation(TObjArray * arr,Float_t fac) const
6532 // if we use TPC track itself we have to "update" covariance
6534 Int_t nseed= arr->GetEntriesFast();
6535 for (Int_t i=0;i<nseed;i++){
6536 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
6540 //rotate to current local system at first accepted point
6541 Int_t index = pt->GetClusterIndex2(pt->GetFirstPoint());
6542 Int_t sec = (index&0xff000000)>>24;
6544 Float_t angle1 = fInnerSec->GetAlpha()*sec+fInnerSec->GetAlphaShift();
6545 if (angle1>TMath::Pi())
6546 angle1-=2.*TMath::Pi();
6547 Float_t angle2 = pt->GetAlpha();
6549 if (TMath::Abs(angle1-angle2)>0.001){
6550 pt->Rotate(angle1-angle2);
6551 //angle2 = pt->GetAlpha();
6552 //pt->fRelativeSector = pt->GetAlpha()/fInnerSec->GetAlpha();
6553 //if (pt->GetAlpha()<0)
6554 // pt->fRelativeSector+=18;
6555 //sec = pt->fRelativeSector;
6564 void AliTPCtrackerMI::PrepareForProlongation(TObjArray * arr, Float_t fac) const
6568 // if we use TPC track itself we have to "update" covariance
6570 Int_t nseed= arr->GetEntriesFast();
6571 for (Int_t i=0;i<nseed;i++){
6572 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
6575 pt->SetFirstPoint(pt->GetLastPoint());
6583 Int_t AliTPCtrackerMI::PropagateBack(TObjArray * arr)
6586 // make back propagation
6588 Int_t nseed= arr->GetEntriesFast();
6589 for (Int_t i=0;i<nseed;i++){
6590 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
6591 if (pt&& pt->GetKinkIndex(0)<=0) {
6592 //AliTPCseed *pt2 = new AliTPCseed(*pt);
6593 fSectors = fInnerSec;
6594 //FollowBackProlongation(*pt,fInnerSec->GetNRows()-1);
6595 //fSectors = fOuterSec;
6596 FollowBackProlongation(*pt,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1);
6597 //if (pt->GetNumberOfClusters()<(pt->fEsd->GetTPCclusters(0)) ){
6598 // Error("PropagateBack","Not prolonged track %d",pt->GetLabel());
6599 // FollowBackProlongation(*pt2,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1);
6602 if (pt&& pt->GetKinkIndex(0)>0) {
6603 AliESDkink * kink = fEvent->GetKink(pt->GetKinkIndex(0)-1);
6604 pt->SetFirstPoint(kink->GetTPCRow0());
6605 fSectors = fInnerSec;
6606 FollowBackProlongation(*pt,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1);
6614 Int_t AliTPCtrackerMI::PropagateForward2(TObjArray * arr)
6617 // make forward propagation
6619 Int_t nseed= arr->GetEntriesFast();
6621 for (Int_t i=0;i<nseed;i++){
6622 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
6624 FollowProlongation(*pt,0);
6633 Int_t AliTPCtrackerMI::PropagateForward()
6636 // propagate track forward
6638 Int_t nseed = fSeeds->GetEntriesFast();
6639 for (Int_t i=0;i<nseed;i++){
6640 AliTPCseed *pt = (AliTPCseed*)fSeeds->UncheckedAt(i);
6642 AliTPCseed &t = *pt;
6643 Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
6644 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
6645 if (alpha < 0. ) alpha += 2.*TMath::Pi();
6646 t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN);
6650 fSectors = fOuterSec;
6651 ParallelTracking(fSeeds,fOuterSec->GetNRows()+fInnerSec->GetNRows()-1,fInnerSec->GetNRows());
6652 fSectors = fInnerSec;
6653 ParallelTracking(fSeeds,fInnerSec->GetNRows()-1,0);
6662 Int_t AliTPCtrackerMI::PropagateBack(AliTPCseed * pt, Int_t row0, Int_t row1)
6665 // make back propagation, in between row0 and row1
6669 fSectors = fInnerSec;
6672 if (row1<fSectors->GetNRows())
6675 r1 = fSectors->GetNRows()-1;
6677 if (row0<fSectors->GetNRows()&& r1>0 )
6678 FollowBackProlongation(*pt,r1);
6679 if (row1<=fSectors->GetNRows())
6682 r1 = row1 - fSectors->GetNRows();
6683 if (r1<=0) return 0;
6684 if (r1>=fOuterSec->GetNRows()) return 0;
6685 fSectors = fOuterSec;
6686 return FollowBackProlongation(*pt,r1);
6694 void AliTPCtrackerMI::GetShape(AliTPCseed * seed, Int_t row)
6698 AliTPCClusterParam * clparam = AliTPCcalibDB::Instance()->GetClusterParam();
6699 Float_t zdrift = TMath::Abs((fParam->GetZLength(0)-TMath::Abs(seed->GetZ())));
6700 Int_t type = (seed->GetSector() < fParam->GetNSector()/2) ? 0: (row>126) ? 1:2;
6701 Double_t angulary = seed->GetSnp();
6702 angulary = angulary*angulary/(1.-angulary*angulary);
6703 Double_t angularz = seed->GetTgl()*seed->GetTgl()*(1.+angulary);
6705 Double_t sigmay = clparam->GetRMS0(0,type,zdrift,TMath::Sqrt(TMath::Abs(angulary)));
6706 Double_t sigmaz = clparam->GetRMS0(1,type,zdrift,TMath::Sqrt(TMath::Abs(angularz)));
6707 seed->SetCurrentSigmaY2(sigmay*sigmay);
6708 seed->SetCurrentSigmaZ2(sigmaz*sigmaz);
6709 // Float_t sd2 = TMath::Abs((fParam->GetZLength(0)-TMath::Abs(seed->GetZ())))*fParam->GetDiffL()*fParam->GetDiffL();
6710 // // Float_t padlength = fParam->GetPadPitchLength(seed->fSector);
6711 // Float_t padlength = GetPadPitchLength(row);
6713 // Float_t sresy = (seed->GetSector() < fParam->GetNSector()/2) ? 0.2 :0.3;
6714 // seed->SetCurrentSigmaY2(sd2+padlength*padlength*angulary/12.+sresy*sresy);
6716 // Float_t sresz = fParam->GetZSigma();
6717 // seed->SetCurrentSigmaZ2(sd2+padlength*padlength*angularz*angularz*(1+angulary)/12.+sresz*sresz);
6719 Float_t wy = GetSigmaY(seed);
6720 Float_t wz = GetSigmaZ(seed);
6723 if (TMath::Abs(wy/seed->fCurrentSigmaY2-1)>0.0001 || TMath::Abs(wz/seed->fCurrentSigmaZ2-1)>0.0001 ){
6724 printf("problem\n");
6731 //__________________________________________________________________________
6732 void AliTPCtrackerMI::CookLabel(AliKalmanTrack *tk, Float_t wrong) const {
6733 //--------------------------------------------------------------------
6734 //This function "cooks" a track label. If label<0, this track is fake.
6735 //--------------------------------------------------------------------
6736 AliTPCseed * t = dynamic_cast<AliTPCseed*>(tk);
6738 printf("%s:%d wrong type \n",(char*)__FILE__,__LINE__);
6742 Int_t noc=t->GetNumberOfClusters();
6744 //printf("\nnot founded prolongation\n\n\n");
6750 AliTPCclusterMI *clusters[160];
6752 for (Int_t i=0;i<160;i++) {
6759 for (i=0; i<160 && current<noc; i++) {
6761 Int_t index=t->GetClusterIndex2(i);
6762 if (index<=0) continue;
6763 if (index&0x8000) continue;
6765 //clusters[current]=GetClusterMI(index);
6766 if (t->GetClusterPointer(i)){
6767 clusters[current]=t->GetClusterPointer(i);
6773 Int_t lab=123456789;
6774 for (i=0; i<noc; i++) {
6775 AliTPCclusterMI *c=clusters[i];
6777 lab=TMath::Abs(c->GetLabel(0));
6779 for (j=0; j<noc; j++) if (lb[j]==lab || mx[j]==0) break;
6785 for (i=0; i<noc; i++) if (mx[i]>max) {max=mx[i]; lab=lb[i];}
6787 for (i=0; i<noc; i++) {
6788 AliTPCclusterMI *c=clusters[i];
6790 if (TMath::Abs(c->GetLabel(1)) == lab ||
6791 TMath::Abs(c->GetLabel(2)) == lab ) max++;
6794 if ((1.- Float_t(max)/noc) > wrong) lab=-lab;
6797 Int_t tail=Int_t(0.10*noc);
6800 for (i=1; i<=160&&ind<tail; i++) {
6801 // AliTPCclusterMI *c=clusters[noc-i];
6802 AliTPCclusterMI *c=clusters[i];
6804 if (lab == TMath::Abs(c->GetLabel(0)) ||
6805 lab == TMath::Abs(c->GetLabel(1)) ||
6806 lab == TMath::Abs(c->GetLabel(2))) max++;
6809 if (max < Int_t(0.5*tail)) lab=-lab;
6816 //delete[] clusters;
6820 //__________________________________________________________________________
6821 Int_t AliTPCtrackerMI::CookLabel(AliTPCseed *t, Float_t wrong,Int_t first, Int_t last) const {
6822 //--------------------------------------------------------------------
6823 //This function "cooks" a track label. If label<0, this track is fake.
6824 //--------------------------------------------------------------------
6825 Int_t noc=t->GetNumberOfClusters();
6827 //printf("\nnot founded prolongation\n\n\n");
6833 AliTPCclusterMI *clusters[160];
6835 for (Int_t i=0;i<160;i++) {
6842 for (i=0; i<160 && current<noc; i++) {
6843 if (i<first) continue;
6844 if (i>last) continue;
6845 Int_t index=t->GetClusterIndex2(i);
6846 if (index<=0) continue;
6847 if (index&0x8000) continue;
6849 //clusters[current]=GetClusterMI(index);
6850 if (t->GetClusterPointer(i)){
6851 clusters[current]=t->GetClusterPointer(i);
6856 if (noc<5) return -1;
6857 Int_t lab=123456789;
6858 for (i=0; i<noc; i++) {
6859 AliTPCclusterMI *c=clusters[i];
6861 lab=TMath::Abs(c->GetLabel(0));
6863 for (j=0; j<noc; j++) if (lb[j]==lab || mx[j]==0) break;
6869 for (i=0; i<noc; i++) if (mx[i]>max) {max=mx[i]; lab=lb[i];}
6871 for (i=0; i<noc; i++) {
6872 AliTPCclusterMI *c=clusters[i];
6874 if (TMath::Abs(c->GetLabel(1)) == lab ||
6875 TMath::Abs(c->GetLabel(2)) == lab ) max++;
6878 if ((1.- Float_t(max)/noc) > wrong) lab=-lab;
6881 Int_t tail=Int_t(0.10*noc);
6884 for (i=1; i<=160&&ind<tail; i++) {
6885 // AliTPCclusterMI *c=clusters[noc-i];
6886 AliTPCclusterMI *c=clusters[i];
6888 if (lab == TMath::Abs(c->GetLabel(0)) ||
6889 lab == TMath::Abs(c->GetLabel(1)) ||
6890 lab == TMath::Abs(c->GetLabel(2))) max++;
6893 if (max < Int_t(0.5*tail)) lab=-lab;
6896 // t->SetLabel(lab);
6900 //delete[] clusters;
6904 Int_t AliTPCtrackerMI::GetRowNumber(Double_t x[3]) const
6906 //return pad row number for given x vector
6907 Float_t phi = TMath::ATan2(x[1],x[0]);
6908 if(phi<0) phi=2.*TMath::Pi()+phi;
6909 // Get the local angle in the sector philoc
6910 const Float_t kRaddeg = 180/3.14159265358979312;
6911 Float_t phiangle = (Int_t (phi*kRaddeg/20.) + 0.5)*20./kRaddeg;
6912 Double_t localx = x[0]*TMath::Cos(phiangle)-x[1]*TMath::Sin(phiangle);
6913 return GetRowNumber(localx);
6918 void AliTPCtrackerMI::MakeBitmaps(AliTPCseed *t)
6920 //-----------------------------------------------------------------------
6921 // Fill the cluster and sharing bitmaps of the track
6922 //-----------------------------------------------------------------------
6924 Int_t firstpoint = 0;
6925 Int_t lastpoint = 159;
6926 AliTPCTrackerPoint *point;
6928 for (int iter=firstpoint; iter<lastpoint; iter++) {
6929 point = t->GetTrackPoint(iter);
6931 t->SetClusterMapBit(iter, kTRUE);
6932 if (point->IsShared())
6933 t->SetSharedMapBit(iter,kTRUE);
6935 t->SetSharedMapBit(iter, kFALSE);
6938 t->SetClusterMapBit(iter, kFALSE);
6939 t->SetSharedMapBit(iter, kFALSE);
6946 void AliTPCtrackerMI::AddCovariance(AliTPCseed * seed){
6948 // Adding systematic error
6949 // !!!! the systematic error for element 4 is in 1/cm not in pt
6951 const Double_t *param = AliTPCReconstructor::GetRecoParam()->GetSystematicError();
6953 for (Int_t i=0;i<15;i++) covar[i]=0;
6959 covar[0] = param[0]*param[0];
6960 covar[2] = param[1]*param[1];
6961 covar[5] = param[2]*param[2];
6962 covar[9] = param[3]*param[3];
6963 Double_t facC = AliTracker::GetBz()*kB2C;
6964 covar[14]= param[4]*param[4]*facC*facC;
6965 seed->AddCovariance(covar);