1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
17 //-------------------------------------------------------
18 // Implementation of the TPC tracker
20 // Origin: Marian Ivanov Marian.Ivanov@cern.ch
22 // AliTPC parallel tracker
24 // The track fitting is based on Kalaman filtering approach
26 // The track finding steps:
27 // 1. Seeding - with and without vertex constraint
28 // - seeding with vertex constain done at first n^2 proble
29 // - seeding without vertex constraint n^3 problem
30 // 2. Tracking - follow prolongation road - find cluster - update kalman track
32 // The seeding and tracking is repeated several times, in different seeding region.
33 // This approach enables to find the track which cannot be seeded in some region of TPC
34 // This can happen because of low momenta (track do not reach outer radius), or track is currently in the ded region between sectors, or the track is for the moment overlapped with other track (seed quality is poor) ...
36 // With this approach we reach almost 100 % efficiency also for high occupancy events.
37 // (If the seeding efficiency in a region is about 90 % than with logical or of several
38 // regions we will reach 100% (in theory - supposing independence)
40 // Repeating several seeding - tracking procedures some of the tracks can be find
43 // The procedures to remove multi find tacks are impremented:
44 // RemoveUsed2 - fast procedure n problem -
45 // Algorithm - Sorting tracks according quality
46 // remove tracks with some shared fraction
47 // Sharing in respect to all tacks
48 // Signing clusters in gold region
49 // FindSplitted - slower algorithm n^2
50 // Sort the tracks according quality
51 // Loop over pair of tracks
52 // If overlap with other track bigger than threshold - remove track
54 // FindCurling - Finds the pair of tracks which are curling
55 // - About 10% of tracks can be find with this procedure
56 // The combinatorial background is too big to be used in High
57 // multiplicity environment
58 // - n^2 problem - Slow procedure - currently it is disabled because of
61 // The number of splitted tracks can be reduced disabling the sharing of the cluster.
62 // tpcRecoParam-> SetClusterSharing(kFALSE);
63 // IT IS HIGHLY non recomended to use it in high flux enviroonment
64 // Even using this switch some tracks can be found more than once
65 // (because of multiple seeding and low quality tracks which will not cross full chamber)
68 // The tracker itself can be debugged - the information about tracks can be stored in several // phases of the reconstruction
69 // To enable storage of the TPC tracks in the ESD friend track
70 // use AliTPCReconstructor::SetStreamLevel(n); where nis bigger 0
72 // The debug level - different procedure produce tree for numerical debugging
73 // To enable them set AliTPCReconstructor::SetStreamLevel(n); where nis bigger 1
77 // Adding systematic errors to the covariance:
79 // The systematic errors due to the misalignment and miscalibration are added to the covariance matrix
80 // of the tracks (not to the clusters as they are dependent):
81 // The parameters form AliTPCRecoParam are used AliTPCRecoParam::GetSystematicError
82 // The systematic errors are expressed there in RMS - position (cm), angle (rad), curvature (1/cm)
83 // The default values are 0.
85 // The sytematic errors are added to the covariance matrix in following places:
87 // 1. During fisrt itteration - AliTPCtrackerMI::FillESD
88 // 2. Second iteration -
89 // 2.a ITS->TPC - AliTPCtrackerMI::ReadSeeds
90 // 2.b TPC->TRD - AliTPCtrackerMI::PropagateBack
91 // 3. Third iteration -
92 // 3.a TRD->TPC - AliTPCtrackerMI::ReadSeeds
93 // 3.b TPC->ITS - AliTPCtrackerMI::RefitInward
95 // There are several places in the code which can be numerically debuged
96 // This code is keeped in order to enable code development and to check the calibration implementtion
98 // 1. ErrParam stream (Log level 9) - dump information about
100 // 2.a) cluster error estimate
101 // 3.a) cluster shape estimate
104 //-------------------------------------------------------
109 #include "Riostream.h"
110 #include <TClonesArray.h>
112 #include <TObjArray.h>
115 #include "AliComplexCluster.h"
116 #include "AliESDEvent.h"
117 #include "AliESDtrack.h"
118 #include "AliESDVertex.h"
121 #include "AliHelix.h"
122 #include "AliRunLoader.h"
123 #include "AliTPCClustersRow.h"
124 #include "AliTPCParam.h"
125 #include "AliTPCReconstructor.h"
126 #include "AliTPCpolyTrack.h"
127 #include "AliTPCreco.h"
128 #include "AliTPCseed.h"
130 #include "AliTPCtrackerSector.h"
131 #include "AliTPCtrackerMI.h"
132 #include "TStopwatch.h"
133 #include "AliTPCReconstructor.h"
135 #include "TTreeStream.h"
136 #include "AliAlignObj.h"
137 #include "AliTrackPointArray.h"
139 #include "AliTPCcalibDB.h"
140 #include "AliTPCTransform.h"
141 #include "AliTPCClusterParam.h"
145 ClassImp(AliTPCtrackerMI)
149 class AliTPCFastMath {
152 static Double_t FastAsin(Double_t x);
154 static Double_t fgFastAsin[20000]; //lookup table for fast asin computation
157 Double_t AliTPCFastMath::fgFastAsin[20000];
158 AliTPCFastMath gAliTPCFastMath; // needed to fill the LUT
160 AliTPCFastMath::AliTPCFastMath(){
162 // initialized lookup table;
163 for (Int_t i=0;i<10000;i++){
164 fgFastAsin[2*i] = TMath::ASin(i/10000.);
165 fgFastAsin[2*i+1] = (TMath::ASin((i+1)/10000.)-fgFastAsin[2*i]);
169 Double_t AliTPCFastMath::FastAsin(Double_t x){
171 // return asin using lookup table
173 Int_t index = int(x*10000);
174 return fgFastAsin[2*index]+(x*10000.-index)*fgFastAsin[2*index+1];
177 Int_t index = int(x*10000);
178 return -(fgFastAsin[2*index]+(x*10000.-index)*fgFastAsin[2*index+1]);
180 //__________________________________________________________________
181 AliTPCtrackerMI::AliTPCtrackerMI()
203 // default constructor
206 //_____________________________________________________________________
210 Int_t AliTPCtrackerMI::UpdateTrack(AliTPCseed * track, Int_t accept){
212 //update track information using current cluster - track->fCurrentCluster
215 AliTPCclusterMI* c =track->GetCurrentCluster();
216 if (accept>0) track->SetCurrentClusterIndex1(track->GetCurrentClusterIndex1() | 0x8000); //sign not accepted clusters
218 UInt_t i = track->GetCurrentClusterIndex1();
220 Int_t sec=(i&0xff000000)>>24;
221 //Int_t row = (i&0x00ff0000)>>16;
222 track->SetRow((i&0x00ff0000)>>16);
223 track->SetSector(sec);
224 // Int_t index = i&0xFFFF;
225 if (sec>=fParam->GetNInnerSector()) track->SetRow(track->GetRow()+fParam->GetNRowLow());
226 track->SetClusterIndex2(track->GetRow(), i);
227 //track->fFirstPoint = row;
228 //if ( track->fLastPoint<row) track->fLastPoint =row;
229 // if (track->fRow<0 || track->fRow>160) {
230 // printf("problem\n");
232 if (track->GetFirstPoint()>track->GetRow())
233 track->SetFirstPoint(track->GetRow());
234 if (track->GetLastPoint()<track->GetRow())
235 track->SetLastPoint(track->GetRow());
238 track->SetClusterPointer(track->GetRow(),c);
241 Double_t angle2 = track->GetSnp()*track->GetSnp();
243 //SET NEW Track Point
245 if (angle2<1) //PH sometimes angle2 is very big. To be investigated...
247 angle2 = TMath::Sqrt(angle2/(1-angle2));
248 AliTPCTrackerPoint &point =*(track->GetTrackPoint(track->GetRow()));
250 point.SetSigmaY(c->GetSigmaY2()/track->GetCurrentSigmaY2());
251 point.SetSigmaZ(c->GetSigmaZ2()/track->GetCurrentSigmaZ2());
252 point.SetErrY(sqrt(track->GetErrorY2()));
253 point.SetErrZ(sqrt(track->GetErrorZ2()));
255 point.SetX(track->GetX());
256 point.SetY(track->GetY());
257 point.SetZ(track->GetZ());
258 point.SetAngleY(angle2);
259 point.SetAngleZ(track->GetTgl());
260 if (point.IsShared()){
261 track->SetErrorY2(track->GetErrorY2()*4);
262 track->SetErrorZ2(track->GetErrorZ2()*4);
266 Double_t chi2 = track->GetPredictedChi2(track->GetCurrentCluster());
268 // track->SetErrorY2(track->GetErrorY2()*1.3);
269 // track->SetErrorY2(track->GetErrorY2()+0.01);
270 // track->SetErrorZ2(track->GetErrorZ2()*1.3);
271 // track->SetErrorZ2(track->GetErrorZ2()+0.005);
273 if (accept>0) return 0;
274 if (track->GetNumberOfClusters()%20==0){
275 // if (track->fHelixIn){
276 // TClonesArray & larr = *(track->fHelixIn);
277 // Int_t ihelix = larr.GetEntriesFast();
278 // new(larr[ihelix]) AliHelix(*track) ;
281 track->SetNoCluster(0);
282 return track->Update(c,chi2,i);
287 Int_t AliTPCtrackerMI::AcceptCluster(AliTPCseed * seed, AliTPCclusterMI * cluster)
290 // decide according desired precision to accept given
291 // cluster for tracking
292 Double_t sy2=ErrY2(seed,cluster);
293 Double_t sz2=ErrZ2(seed,cluster);
295 Double_t sdistancey2 = sy2+seed->GetSigmaY2();
296 Double_t sdistancez2 = sz2+seed->GetSigmaZ2();
298 Double_t rdistancey2 = (seed->GetCurrentCluster()->GetY()-seed->GetY())*
299 (seed->GetCurrentCluster()->GetY()-seed->GetY())/sdistancey2;
300 Double_t rdistancez2 = (seed->GetCurrentCluster()->GetZ()-seed->GetZ())*
301 (seed->GetCurrentCluster()->GetZ()-seed->GetZ())/sdistancez2;
303 Double_t rdistance2 = rdistancey2+rdistancez2;
306 if (AliTPCReconstructor::StreamLevel()>5 && seed->GetNumberOfClusters()>20) {
307 Float_t rmsy2 = seed->GetCurrentSigmaY2();
308 Float_t rmsz2 = seed->GetCurrentSigmaZ2();
309 Float_t rmsy2p30 = seed->GetCMeanSigmaY2p30();
310 Float_t rmsz2p30 = seed->GetCMeanSigmaZ2p30();
311 Float_t rmsy2p30R = seed->GetCMeanSigmaY2p30R();
312 Float_t rmsz2p30R = seed->GetCMeanSigmaZ2p30R();
314 AliExternalTrackParam param(*seed);
315 static TVectorD gcl(3),gtr(3);
317 param.GetXYZ(gcl.GetMatrixArray());
318 cluster->GetGlobalXYZ(gclf);
319 gcl[0]=gclf[0]; gcl[1]=gclf[1]; gcl[2]=gclf[2];
321 if (AliTPCReconstructor::StreamLevel()>0) {
322 (*fDebugStreamer)<<"ErrParam"<<
331 "rmsy2p30="<<rmsy2p30<<
332 "rmsz2p30="<<rmsz2p30<<
333 "rmsy2p30R="<<rmsy2p30R<<
334 "rmsz2p30R="<<rmsz2p30R<<
339 if (rdistance2>16) return 3;
342 if ((rdistancey2>9. || rdistancez2>9.) && cluster->GetType()==0)
343 return 2; //suspisiouce - will be changed
345 if ((rdistancey2>6.25 || rdistancez2>6.25) && cluster->GetType()>0)
346 // strict cut on overlaped cluster
347 return 2; //suspisiouce - will be changed
349 if ( (rdistancey2>1. || rdistancez2>6.25 )
350 && cluster->GetType()<0){
351 seed->SetNFoundable(seed->GetNFoundable()-1);
360 //_____________________________________________________________________________
361 AliTPCtrackerMI::AliTPCtrackerMI(const AliTPCParam *par):
363 fkNIS(par->GetNInnerSector()/2),
365 fkNOS(par->GetNOuterSector()/2),
382 //---------------------------------------------------------------------
383 // The main TPC tracker constructor
384 //---------------------------------------------------------------------
385 fInnerSec=new AliTPCtrackerSector[fkNIS];
386 fOuterSec=new AliTPCtrackerSector[fkNOS];
389 for (i=0; i<fkNIS; i++) fInnerSec[i].Setup(par,0);
390 for (i=0; i<fkNOS; i++) fOuterSec[i].Setup(par,1);
393 Int_t nrowlow = par->GetNRowLow();
394 Int_t nrowup = par->GetNRowUp();
397 for (i=0;i<nrowlow;i++){
398 fXRow[i] = par->GetPadRowRadiiLow(i);
399 fPadLength[i]= par->GetPadPitchLength(0,i);
400 fYMax[i] = fXRow[i]*TMath::Tan(0.5*par->GetInnerAngle());
404 for (i=0;i<nrowup;i++){
405 fXRow[i+nrowlow] = par->GetPadRowRadiiUp(i);
406 fPadLength[i+nrowlow] = par->GetPadPitchLength(60,i);
407 fYMax[i+nrowlow] = fXRow[i+nrowlow]*TMath::Tan(0.5*par->GetOuterAngle());
410 if (AliTPCReconstructor::StreamLevel()>0) {
411 fDebugStreamer = new TTreeSRedirector("TPCdebug.root");
414 //________________________________________________________________________
415 AliTPCtrackerMI::AliTPCtrackerMI(const AliTPCtrackerMI &t):
436 //------------------------------------
437 // dummy copy constructor
438 //------------------------------------------------------------------
441 AliTPCtrackerMI & AliTPCtrackerMI::operator=(const AliTPCtrackerMI& /*r*/){
442 //------------------------------
444 //--------------------------------------------------------------
447 //_____________________________________________________________________________
448 AliTPCtrackerMI::~AliTPCtrackerMI() {
449 //------------------------------------------------------------------
450 // TPC tracker destructor
451 //------------------------------------------------------------------
458 if (fDebugStreamer) delete fDebugStreamer;
462 void AliTPCtrackerMI::FillESD(TObjArray* arr)
466 //fill esds using updated tracks
468 // write tracks to the event
469 // store index of the track
470 Int_t nseed=arr->GetEntriesFast();
471 //FindKinks(arr,fEvent);
472 for (Int_t i=0; i<nseed; i++) {
473 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
477 if (AliTPCReconstructor::StreamLevel()>1) {
478 (*fDebugStreamer)<<"Track0"<<
482 // pt->PropagateTo(fParam->GetInnerRadiusLow());
483 if (pt->GetKinkIndex(0)<=0){ //don't propagate daughter tracks
484 pt->PropagateTo(fParam->GetInnerRadiusLow());
487 if (( pt->GetPoints()[2]- pt->GetPoints()[0])>5 && pt->GetPoints()[3]>0.8){
489 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
490 iotrack.SetTPCPoints(pt->GetPoints());
491 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
492 iotrack.SetV0Indexes(pt->GetV0Indexes());
493 // iotrack.SetTPCpid(pt->fTPCr);
494 //iotrack.SetTPCindex(i);
495 fEvent->AddTrack(&iotrack);
499 if ( (pt->GetNumberOfClusters()>70)&& (Float_t(pt->GetNumberOfClusters())/Float_t(pt->GetNFoundable()))>0.55) {
501 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
502 iotrack.SetTPCPoints(pt->GetPoints());
503 //iotrack.SetTPCindex(i);
504 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
505 iotrack.SetV0Indexes(pt->GetV0Indexes());
506 // iotrack.SetTPCpid(pt->fTPCr);
507 fEvent->AddTrack(&iotrack);
511 // short tracks - maybe decays
513 if ( (pt->GetNumberOfClusters()>30) && (Float_t(pt->GetNumberOfClusters())/Float_t(pt->GetNFoundable()))>0.70) {
514 Int_t found,foundable,shared;
515 pt->GetClusterStatistic(0,60,found, foundable,shared,kFALSE);
516 if ( (found>20) && (pt->GetNShared()/float(pt->GetNumberOfClusters())<0.2)){
518 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
519 //iotrack.SetTPCindex(i);
520 iotrack.SetTPCPoints(pt->GetPoints());
521 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
522 iotrack.SetV0Indexes(pt->GetV0Indexes());
523 //iotrack.SetTPCpid(pt->fTPCr);
524 fEvent->AddTrack(&iotrack);
529 if ( (pt->GetNumberOfClusters()>20) && (Float_t(pt->GetNumberOfClusters())/Float_t(pt->GetNFoundable()))>0.8) {
530 Int_t found,foundable,shared;
531 pt->GetClusterStatistic(0,60,found, foundable,shared,kFALSE);
532 if (found<20) continue;
533 if (pt->GetNShared()/float(pt->GetNumberOfClusters())>0.2) continue;
536 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
537 iotrack.SetTPCPoints(pt->GetPoints());
538 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
539 iotrack.SetV0Indexes(pt->GetV0Indexes());
540 //iotrack.SetTPCpid(pt->fTPCr);
541 //iotrack.SetTPCindex(i);
542 fEvent->AddTrack(&iotrack);
545 // short tracks - secondaties
547 if ( (pt->GetNumberOfClusters()>30) ) {
548 Int_t found,foundable,shared;
549 pt->GetClusterStatistic(128,158,found, foundable,shared,kFALSE);
550 if ( (found>20) && (pt->GetNShared()/float(pt->GetNumberOfClusters())<0.2) &&float(found)/float(foundable)>0.8){
552 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
553 iotrack.SetTPCPoints(pt->GetPoints());
554 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
555 iotrack.SetV0Indexes(pt->GetV0Indexes());
556 //iotrack.SetTPCpid(pt->fTPCr);
557 //iotrack.SetTPCindex(i);
558 fEvent->AddTrack(&iotrack);
563 if ( (pt->GetNumberOfClusters()>15)) {
564 Int_t found,foundable,shared;
565 pt->GetClusterStatistic(138,158,found, foundable,shared,kFALSE);
566 if (found<15) continue;
567 if (foundable<=0) continue;
568 if (pt->GetNShared()/float(pt->GetNumberOfClusters())>0.2) continue;
569 if (float(found)/float(foundable)<0.8) continue;
572 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
573 iotrack.SetTPCPoints(pt->GetPoints());
574 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
575 iotrack.SetV0Indexes(pt->GetV0Indexes());
576 // iotrack.SetTPCpid(pt->fTPCr);
577 //iotrack.SetTPCindex(i);
578 fEvent->AddTrack(&iotrack);
583 printf("Number of filled ESDs-\t%d\n",fEvent->GetNumberOfTracks());
590 Double_t AliTPCtrackerMI::ErrY2(AliTPCseed* seed, AliTPCclusterMI * cl){
593 // Use calibrated cluster error from OCDB
595 AliTPCClusterParam * clparam = AliTPCcalibDB::Instance()->GetClusterParam();
597 Float_t z = TMath::Abs(fParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
598 Int_t ctype = cl->GetType();
599 Int_t type = (cl->GetRow()<63) ? 0: (cl->GetRow()>126) ? 1:2;
600 Double_t angle = seed->GetSnp()*seed->GetSnp();
601 angle = TMath::Sqrt(TMath::Abs(angle/(1.-angle)));
602 Double_t erry2 = clparam->GetError0Par(0,type, z,angle);
604 erry2+=0.5; // edge cluster
607 seed->SetErrorY2(erry2);
611 //calculate look-up table at the beginning
612 // static Bool_t ginit = kFALSE;
613 // static Float_t gnoise1,gnoise2,gnoise3;
614 // static Float_t ggg1[10000];
615 // static Float_t ggg2[10000];
616 // static Float_t ggg3[10000];
617 // static Float_t glandau1[10000];
618 // static Float_t glandau2[10000];
619 // static Float_t glandau3[10000];
621 // static Float_t gcor01[500];
622 // static Float_t gcor02[500];
623 // static Float_t gcorp[500];
627 // if (ginit==kFALSE){
628 // for (Int_t i=1;i<500;i++){
629 // Float_t rsigma = float(i)/100.;
630 // gcor02[i] = TMath::Max(0.78 +TMath::Exp(7.4*(rsigma-1.2)),0.6);
631 // gcor01[i] = TMath::Max(0.72 +TMath::Exp(3.36*(rsigma-1.2)),0.6);
632 // gcorp[i] = TMath::Max(TMath::Power((rsigma+0.5),1.5),1.2);
636 // for (Int_t i=3;i<10000;i++){
640 // Float_t amp = float(i);
641 // Float_t padlength =0.75;
642 // gnoise1 = 0.0004/padlength;
643 // Float_t nel = 0.268*amp;
644 // Float_t nprim = 0.155*amp;
645 // ggg1[i] = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.001*nel/(padlength*padlength))/nel;
646 // glandau1[i] = (2.+0.12*nprim)*0.5* (2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
647 // if (glandau1[i]>1) glandau1[i]=1;
648 // glandau1[i]*=padlength*padlength/12.;
652 // gnoise2 = 0.0004/padlength;
654 // nprim = 0.133*amp;
655 // ggg2[i] = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
656 // glandau2[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
657 // if (glandau2[i]>1) glandau2[i]=1;
658 // glandau2[i]*=padlength*padlength/12.;
663 // gnoise3 = 0.0004/padlength;
665 // nprim = 0.133*amp;
666 // ggg3[i] = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
667 // glandau3[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
668 // if (glandau3[i]>1) glandau3[i]=1;
669 // glandau3[i]*=padlength*padlength/12.;
677 // Int_t amp = int(TMath::Abs(cl->GetQ()));
679 // seed->SetErrorY2(1.);
683 // Float_t z = TMath::Abs(fParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
684 // Int_t ctype = cl->GetType();
685 // Float_t padlength= GetPadPitchLength(seed->GetRow());
686 // Double_t angle2 = seed->GetSnp()*seed->GetSnp();
687 // angle2 = angle2/(1-angle2);
689 // //cluster "quality"
690 // Int_t rsigmay = int(100.*cl->GetSigmaY2()/(seed->GetCurrentSigmaY2()));
693 // if (fSectors==fInnerSec){
694 // snoise2 = gnoise1;
695 // res = ggg1[amp]*z+glandau1[amp]*angle2;
696 // if (ctype==0) res *= gcor01[rsigmay];
699 // res*= gcorp[rsigmay];
703 // if (padlength<1.1){
704 // snoise2 = gnoise2;
705 // res = ggg2[amp]*z+glandau2[amp]*angle2;
706 // if (ctype==0) res *= gcor02[rsigmay];
709 // res*= gcorp[rsigmay];
713 // snoise2 = gnoise3;
714 // res = ggg3[amp]*z+glandau3[amp]*angle2;
715 // if (ctype==0) res *= gcor02[rsigmay];
718 // res*= gcorp[rsigmay];
725 // res*=2.4; // overestimate error 2 times
729 // if (res<2*snoise2)
732 // seed->SetErrorY2(res);
740 Double_t AliTPCtrackerMI::ErrZ2(AliTPCseed* seed, AliTPCclusterMI * cl){
743 // Use calibrated cluster error from OCDB
745 AliTPCClusterParam * clparam = AliTPCcalibDB::Instance()->GetClusterParam();
747 Float_t z = TMath::Abs(fParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
748 Int_t ctype = cl->GetType();
749 Int_t type = (cl->GetRow()<63) ? 0: (cl->GetRow()>126) ? 1:2;
751 Double_t angle2 = seed->GetSnp()*seed->GetSnp();
752 angle2 = seed->GetTgl()*seed->GetTgl()*(1+angle2/(1-angle2));
753 Double_t angle = TMath::Sqrt(TMath::Abs(angle2));
754 Double_t errz2 = clparam->GetError0Par(1,type, z,angle);
756 errz2+=0.5; // edge cluster
759 seed->SetErrorZ2(errz2);
765 // //seed->SetErrorY2(0.1);
767 // //calculate look-up table at the beginning
768 // static Bool_t ginit = kFALSE;
769 // static Float_t gnoise1,gnoise2,gnoise3;
770 // static Float_t ggg1[10000];
771 // static Float_t ggg2[10000];
772 // static Float_t ggg3[10000];
773 // static Float_t glandau1[10000];
774 // static Float_t glandau2[10000];
775 // static Float_t glandau3[10000];
777 // static Float_t gcor01[1000];
778 // static Float_t gcor02[1000];
779 // static Float_t gcorp[1000];
783 // if (ginit==kFALSE){
784 // for (Int_t i=1;i<1000;i++){
785 // Float_t rsigma = float(i)/100.;
786 // gcor02[i] = TMath::Max(0.81 +TMath::Exp(6.8*(rsigma-1.2)),0.6);
787 // gcor01[i] = TMath::Max(0.72 +TMath::Exp(2.04*(rsigma-1.2)),0.6);
788 // gcorp[i] = TMath::Max(TMath::Power((rsigma+0.5),1.5),1.2);
792 // for (Int_t i=3;i<10000;i++){
796 // Float_t amp = float(i);
797 // Float_t padlength =0.75;
798 // gnoise1 = 0.0004/padlength;
799 // Float_t nel = 0.268*amp;
800 // Float_t nprim = 0.155*amp;
801 // ggg1[i] = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.001*nel/(padlength*padlength))/nel;
802 // glandau1[i] = (2.+0.12*nprim)*0.5* (2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
803 // if (glandau1[i]>1) glandau1[i]=1;
804 // glandau1[i]*=padlength*padlength/12.;
808 // gnoise2 = 0.0004/padlength;
810 // nprim = 0.133*amp;
811 // ggg2[i] = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
812 // glandau2[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
813 // if (glandau2[i]>1) glandau2[i]=1;
814 // glandau2[i]*=padlength*padlength/12.;
819 // gnoise3 = 0.0004/padlength;
821 // nprim = 0.133*amp;
822 // ggg3[i] = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
823 // glandau3[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
824 // if (glandau3[i]>1) glandau3[i]=1;
825 // glandau3[i]*=padlength*padlength/12.;
833 // Int_t amp = int(TMath::Abs(cl->GetQ()));
835 // seed->SetErrorY2(1.);
839 // Float_t z = TMath::Abs(fParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
840 // Int_t ctype = cl->GetType();
841 // Float_t padlength= GetPadPitchLength(seed->GetRow());
843 // Double_t angle2 = seed->GetSnp()*seed->GetSnp();
844 // // if (angle2<0.6) angle2 = 0.6;
845 // angle2 = seed->GetTgl()*seed->GetTgl()*(1+angle2/(1-angle2));
847 // //cluster "quality"
848 // Int_t rsigmaz = int(100.*cl->GetSigmaZ2()/(seed->GetCurrentSigmaZ2()));
851 // if (fSectors==fInnerSec){
852 // snoise2 = gnoise1;
853 // res = ggg1[amp]*z+glandau1[amp]*angle2;
854 // if (ctype==0) res *= gcor01[rsigmaz];
857 // res*= gcorp[rsigmaz];
861 // if (padlength<1.1){
862 // snoise2 = gnoise2;
863 // res = ggg2[amp]*z+glandau2[amp]*angle2;
864 // if (ctype==0) res *= gcor02[rsigmaz];
867 // res*= gcorp[rsigmaz];
871 // snoise2 = gnoise3;
872 // res = ggg3[amp]*z+glandau3[amp]*angle2;
873 // if (ctype==0) res *= gcor02[rsigmaz];
876 // res*= gcorp[rsigmaz];
885 // if ((ctype<0) &&<70){
890 // if (res<2*snoise2)
892 // if (res>3) res =3;
893 // seed->SetErrorZ2(res);
901 void AliTPCtrackerMI::RotateToLocal(AliTPCseed *seed)
903 //rotate to track "local coordinata
904 Float_t x = seed->GetX();
905 Float_t y = seed->GetY();
906 Float_t ymax = x*TMath::Tan(0.5*fSectors->GetAlpha());
909 seed->SetRelativeSector((seed->GetRelativeSector()+1) % fN);
910 if (!seed->Rotate(fSectors->GetAlpha()))
912 } else if (y <-ymax) {
913 seed->SetRelativeSector((seed->GetRelativeSector()-1+fN) % fN);
914 if (!seed->Rotate(-fSectors->GetAlpha()))
922 //_____________________________________________________________________________
923 Double_t AliTPCtrackerMI::F1old(Double_t x1,Double_t y1,
924 Double_t x2,Double_t y2,
925 Double_t x3,Double_t y3)
927 //-----------------------------------------------------------------
928 // Initial approximation of the track curvature
929 //-----------------------------------------------------------------
930 Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
931 Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
932 (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
933 Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
934 (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
936 Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
937 if ( xr*xr+yr*yr<=0.00000000000001) return 100;
938 return -xr*yr/sqrt(xr*xr+yr*yr);
943 //_____________________________________________________________________________
944 Double_t AliTPCtrackerMI::F1(Double_t x1,Double_t y1,
945 Double_t x2,Double_t y2,
946 Double_t x3,Double_t y3)
948 //-----------------------------------------------------------------
949 // Initial approximation of the track curvature
950 //-----------------------------------------------------------------
956 Double_t det = x3*y2-x2*y3;
961 Double_t u = 0.5* (x2*(x2-x3)+y2*(y2-y3))/det;
962 Double_t x0 = x3*0.5-y3*u;
963 Double_t y0 = y3*0.5+x3*u;
964 Double_t c2 = 1/TMath::Sqrt(x0*x0+y0*y0);
970 Double_t AliTPCtrackerMI::F2(Double_t x1,Double_t y1,
971 Double_t x2,Double_t y2,
972 Double_t x3,Double_t y3)
974 //-----------------------------------------------------------------
975 // Initial approximation of the track curvature
976 //-----------------------------------------------------------------
982 Double_t det = x3*y2-x2*y3;
987 Double_t u = 0.5* (x2*(x2-x3)+y2*(y2-y3))/det;
988 Double_t x0 = x3*0.5-y3*u;
989 Double_t y0 = y3*0.5+x3*u;
990 Double_t c2 = 1/TMath::Sqrt(x0*x0+y0*y0);
999 //_____________________________________________________________________________
1000 Double_t AliTPCtrackerMI::F2old(Double_t x1,Double_t y1,
1001 Double_t x2,Double_t y2,
1002 Double_t x3,Double_t y3)
1004 //-----------------------------------------------------------------
1005 // Initial approximation of the track curvature times center of curvature
1006 //-----------------------------------------------------------------
1007 Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
1008 Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
1009 (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
1010 Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
1011 (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
1013 Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
1015 return -a/(d*y1-b)*xr/sqrt(xr*xr+yr*yr);
1018 //_____________________________________________________________________________
1019 Double_t AliTPCtrackerMI::F3(Double_t x1,Double_t y1,
1020 Double_t x2,Double_t y2,
1021 Double_t z1,Double_t z2)
1023 //-----------------------------------------------------------------
1024 // Initial approximation of the tangent of the track dip angle
1025 //-----------------------------------------------------------------
1026 return (z1 - z2)/sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
1030 Double_t AliTPCtrackerMI::F3n(Double_t x1,Double_t y1,
1031 Double_t x2,Double_t y2,
1032 Double_t z1,Double_t z2, Double_t c)
1034 //-----------------------------------------------------------------
1035 // Initial approximation of the tangent of the track dip angle
1036 //-----------------------------------------------------------------
1040 //angle1 = (z1-z2)*c/(TMath::ASin(c*x1-ni)-TMath::ASin(c*x2-ni));
1042 Double_t d = TMath::Sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
1043 if (TMath::Abs(d*c*0.5)>1) return 0;
1044 // Double_t angle2 = TMath::ASin(d*c*0.5);
1045 // Double_t angle2 = AliTPCFastMath::FastAsin(d*c*0.5);
1046 Double_t angle2 = (d*c*0.5>0.1)? TMath::ASin(d*c*0.5): AliTPCFastMath::FastAsin(d*c*0.5);
1048 angle2 = (z1-z2)*c/(angle2*2.);
1052 Bool_t AliTPCtrackerMI::GetProlongation(Double_t x1, Double_t x2, Double_t x[5], Double_t &y, Double_t &z)
1053 {//-----------------------------------------------------------------
1054 // This function find proloncation of a track to a reference plane x=x2.
1055 //-----------------------------------------------------------------
1059 if (TMath::Abs(x[4]*x1 - x[2]) >= 0.999) {
1063 Double_t c1=x[4]*x1 - x[2], r1=sqrt((1.-c1)*(1.+c1));
1064 Double_t c2=x[4]*x2 - x[2], r2=sqrt((1.-c2)*(1.+c2));
1068 Double_t dy = dx*(c1+c2)/(r1+r2);
1071 Double_t delta = x[4]*dx*(c1+c2)/(c1*r2 + c2*r1);
1073 if (TMath::Abs(delta)>0.01){
1074 dz = x[3]*TMath::ASin(delta)/x[4];
1076 dz = x[3]*AliTPCFastMath::FastAsin(delta)/x[4];
1079 //dz = x[3]*AliTPCFastMath::FastAsin(delta)/x[4];
1087 Int_t AliTPCtrackerMI::LoadClusters (TTree *tree)
1092 return LoadClusters();
1096 Int_t AliTPCtrackerMI::LoadClusters(TObjArray *arr)
1099 // load clusters to the memory
1100 AliTPCClustersRow *clrow = 0x0;
1101 Int_t lower = arr->LowerBound();
1102 Int_t entries = arr->GetEntriesFast();
1104 for (Int_t i=lower; i<entries; i++) {
1105 clrow = (AliTPCClustersRow*) arr->At(i);
1106 if(!clrow) continue;
1107 if(!clrow->GetArray()) continue;
1111 fParam->AdjustSectorRow(clrow->GetID(),sec,row);
1113 for (Int_t icl=0; icl<clrow->GetArray()->GetEntriesFast(); icl++){
1114 Transform((AliTPCclusterMI*)(clrow->GetArray()->At(icl)));
1117 if (clrow->GetArray()->GetEntriesFast()<=0) continue;
1118 AliTPCtrackerRow * tpcrow=0;
1121 tpcrow = &(fInnerSec[sec%fkNIS][row]);
1125 tpcrow = &(fOuterSec[(sec-fkNIS*2)%fkNOS][row]);
1126 left = (sec-fkNIS*2)/fkNOS;
1129 tpcrow->SetN1(clrow->GetArray()->GetEntriesFast());
1130 tpcrow->SetClusters1(new AliTPCclusterMI[tpcrow->GetN1()]);
1131 for (Int_t j=0;j<tpcrow->GetN1();++j)
1132 tpcrow->SetCluster1(j, *(AliTPCclusterMI*)(clrow->GetArray()->At(j)));
1135 tpcrow->SetN2(clrow->GetArray()->GetEntriesFast());
1136 tpcrow->SetClusters2(new AliTPCclusterMI[tpcrow->GetN2()]);
1137 for (Int_t j=0;j<tpcrow->GetN2();++j)
1138 tpcrow->SetCluster2(j,*(AliTPCclusterMI*)(clrow->GetArray()->At(j)));
1148 Int_t AliTPCtrackerMI::LoadClusters(TClonesArray *arr)
1151 // load clusters to the memory from one
1154 AliTPCclusterMI *clust=0;
1155 Int_t count[72][96] = { {0} , {0} };
1157 // loop over clusters
1158 for (Int_t icl=0; icl<arr->GetEntriesFast(); icl++) {
1159 clust = (AliTPCclusterMI*)arr->At(icl);
1160 if(!clust) continue;
1161 //printf("cluster: det %d, row %d \n", clust->GetDetector(),clust->GetRow());
1163 // transform clusters
1166 // count clusters per pad row
1167 count[clust->GetDetector()][clust->GetRow()]++;
1170 // insert clusters to sectors
1171 for (Int_t icl=0; icl<arr->GetEntriesFast(); icl++) {
1172 clust = (AliTPCclusterMI*)arr->At(icl);
1173 if(!clust) continue;
1175 Int_t sec = clust->GetDetector();
1176 Int_t row = clust->GetRow();
1178 // filter overlapping pad rows needed by HLT
1179 if(sec<fkNIS*2) { //IROCs
1180 if(row == 30) continue;
1183 if(row == 27 || row == 76) continue;
1189 fInnerSec[sec%fkNIS].InsertCluster(clust, count[sec][row], fParam);
1192 left = (sec-fkNIS*2)/fkNOS;
1193 fOuterSec[(sec-fkNIS*2)%fkNOS].InsertCluster(clust, count[sec][row], fParam);
1197 // Load functions must be called behind LoadCluster(TClonesArray*)
1199 //LoadOuterSectors();
1200 //LoadInnerSectors();
1206 Int_t AliTPCtrackerMI::LoadClusters()
1209 // load clusters to the memory
1210 AliTPCClustersRow *clrow= new AliTPCClustersRow;
1211 clrow->SetClass("AliTPCclusterMI");
1213 clrow->GetArray()->ExpandCreateFast(10000);
1215 // TTree * tree = fClustersArray.GetTree();
1217 TTree * tree = fInput;
1218 TBranch * br = tree->GetBranch("Segment");
1219 br->SetAddress(&clrow);
1221 Int_t j=Int_t(tree->GetEntries());
1222 for (Int_t i=0; i<j; i++) {
1226 fParam->AdjustSectorRow(clrow->GetID(),sec,row);
1227 for (Int_t icl=0; icl<clrow->GetArray()->GetEntriesFast(); icl++){
1228 Transform((AliTPCclusterMI*)(clrow->GetArray()->At(icl)));
1231 AliTPCtrackerRow * tpcrow=0;
1234 tpcrow = &(fInnerSec[sec%fkNIS][row]);
1238 tpcrow = &(fOuterSec[(sec-fkNIS*2)%fkNOS][row]);
1239 left = (sec-fkNIS*2)/fkNOS;
1242 tpcrow->SetN1(clrow->GetArray()->GetEntriesFast());
1243 tpcrow->SetClusters1(new AliTPCclusterMI[tpcrow->GetN1()]);
1244 for (Int_t k=0;k<tpcrow->GetN1();++k)
1245 tpcrow->SetCluster1(k, *(AliTPCclusterMI*)(clrow->GetArray()->At(k)));
1248 tpcrow->SetN2(clrow->GetArray()->GetEntriesFast());
1249 tpcrow->SetClusters2(new AliTPCclusterMI[tpcrow->GetN2()]);
1250 for (Int_t k=0;k<tpcrow->GetN2();++k)
1251 tpcrow->SetCluster2(k,*(AliTPCclusterMI*)(clrow->GetArray()->At(k)));
1262 void AliTPCtrackerMI::UnloadClusters()
1265 // unload clusters from the memory
1267 Int_t nrows = fOuterSec->GetNRows();
1268 for (Int_t sec = 0;sec<fkNOS;sec++)
1269 for (Int_t row = 0;row<nrows;row++){
1270 AliTPCtrackerRow* tpcrow = &(fOuterSec[sec%fkNOS][row]);
1272 // if (tpcrow->fClusters1) delete []tpcrow->fClusters1;
1273 // if (tpcrow->fClusters2) delete []tpcrow->fClusters2;
1275 tpcrow->ResetClusters();
1278 nrows = fInnerSec->GetNRows();
1279 for (Int_t sec = 0;sec<fkNIS;sec++)
1280 for (Int_t row = 0;row<nrows;row++){
1281 AliTPCtrackerRow* tpcrow = &(fInnerSec[sec%fkNIS][row]);
1283 // if (tpcrow->fClusters1) delete []tpcrow->fClusters1;
1284 //if (tpcrow->fClusters2) delete []tpcrow->fClusters2;
1286 tpcrow->ResetClusters();
1292 void AliTPCtrackerMI::FillClusterArray(TObjArray* array) const{
1294 // Filling cluster to the array - For visualization purposes
1297 nrows = fOuterSec->GetNRows();
1298 for (Int_t sec = 0;sec<fkNOS;sec++)
1299 for (Int_t row = 0;row<nrows;row++){
1300 AliTPCtrackerRow* tpcrow = &(fOuterSec[sec%fkNOS][row]);
1301 if (!tpcrow) continue;
1302 for (Int_t icl = 0;icl<tpcrow->GetN();icl++){
1303 array->AddLast((TObject*)((*tpcrow)[icl]));
1306 nrows = fInnerSec->GetNRows();
1307 for (Int_t sec = 0;sec<fkNIS;sec++)
1308 for (Int_t row = 0;row<nrows;row++){
1309 AliTPCtrackerRow* tpcrow = &(fInnerSec[sec%fkNIS][row]);
1310 if (!tpcrow) continue;
1311 for (Int_t icl = 0;icl<tpcrow->GetN();icl++){
1312 array->AddLast((TObject*)(*tpcrow)[icl]);
1318 void AliTPCtrackerMI::Transform(AliTPCclusterMI * cluster){
1322 AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ;
1324 AliFatal("Tranformations not in calibDB");
1326 Double_t x[3]={cluster->GetRow(),cluster->GetPad(),cluster->GetTimeBin()};
1327 Int_t i[1]={cluster->GetDetector()};
1328 transform->Transform(x,i,0,1);
1329 // if (cluster->GetDetector()%36>17){
1334 // in debug mode check the transformation
1336 if (AliTPCReconstructor::StreamLevel()>1) {
1338 cluster->GetGlobalXYZ(gx);
1339 Int_t event = (fEvent==NULL)? 0: fEvent->GetEventNumberInFile();
1340 TTreeSRedirector &cstream = *fDebugStreamer;
1341 cstream<<"Transform"<<
1352 cluster->SetX(x[0]);
1353 cluster->SetY(x[1]);
1354 cluster->SetZ(x[2]);
1360 //if (!fParam->IsGeoRead()) fParam->ReadGeoMatrices();
1361 TGeoHMatrix *mat = fParam->GetClusterMatrix(cluster->GetDetector());
1363 Double_t pos[3]= {cluster->GetX(),cluster->GetY(),cluster->GetZ()};
1364 Double_t posC[3]={cluster->GetX(),cluster->GetY(),cluster->GetZ()};
1365 if (mat) mat->LocalToMaster(pos,posC);
1367 // chack Loading of Geo matrices from GeoManager - TEMPORARY FIX
1369 cluster->SetX(posC[0]);
1370 cluster->SetY(posC[1]);
1371 cluster->SetZ(posC[2]);
1374 //_____________________________________________________________________________
1375 Int_t AliTPCtrackerMI::LoadOuterSectors() {
1376 //-----------------------------------------------------------------
1377 // This function fills outer TPC sectors with clusters.
1378 //-----------------------------------------------------------------
1379 Int_t nrows = fOuterSec->GetNRows();
1381 for (Int_t sec = 0;sec<fkNOS;sec++)
1382 for (Int_t row = 0;row<nrows;row++){
1383 AliTPCtrackerRow* tpcrow = &(fOuterSec[sec%fkNOS][row]);
1384 Int_t sec2 = sec+2*fkNIS;
1386 Int_t ncl = tpcrow->GetN1();
1388 AliTPCclusterMI *c= (tpcrow->GetCluster1(ncl));
1389 index=(((sec2<<8)+row)<<16)+ncl;
1390 tpcrow->InsertCluster(c,index);
1393 ncl = tpcrow->GetN2();
1395 AliTPCclusterMI *c= (tpcrow->GetCluster2(ncl));
1396 index=((((sec2+fkNOS)<<8)+row)<<16)+ncl;
1397 tpcrow->InsertCluster(c,index);
1400 // write indexes for fast acces
1402 for (Int_t i=0;i<510;i++)
1403 tpcrow->SetFastCluster(i,-1);
1404 for (Int_t i=0;i<tpcrow->GetN();i++){
1405 Int_t zi = Int_t((*tpcrow)[i]->GetZ()+255.);
1406 tpcrow->SetFastCluster(zi,i); // write index
1409 for (Int_t i=0;i<510;i++){
1410 if (tpcrow->GetFastCluster(i)<0)
1411 tpcrow->SetFastCluster(i,last);
1413 last = tpcrow->GetFastCluster(i);
1422 //_____________________________________________________________________________
1423 Int_t AliTPCtrackerMI::LoadInnerSectors() {
1424 //-----------------------------------------------------------------
1425 // This function fills inner TPC sectors with clusters.
1426 //-----------------------------------------------------------------
1427 Int_t nrows = fInnerSec->GetNRows();
1429 for (Int_t sec = 0;sec<fkNIS;sec++)
1430 for (Int_t row = 0;row<nrows;row++){
1431 AliTPCtrackerRow* tpcrow = &(fInnerSec[sec%fkNIS][row]);
1434 Int_t ncl = tpcrow->GetN1();
1436 AliTPCclusterMI *c= (tpcrow->GetCluster1(ncl));
1437 index=(((sec<<8)+row)<<16)+ncl;
1438 tpcrow->InsertCluster(c,index);
1441 ncl = tpcrow->GetN2();
1443 AliTPCclusterMI *c= (tpcrow->GetCluster2(ncl));
1444 index=((((sec+fkNIS)<<8)+row)<<16)+ncl;
1445 tpcrow->InsertCluster(c,index);
1448 // write indexes for fast acces
1450 for (Int_t i=0;i<510;i++)
1451 tpcrow->SetFastCluster(i,-1);
1452 for (Int_t i=0;i<tpcrow->GetN();i++){
1453 Int_t zi = Int_t((*tpcrow)[i]->GetZ()+255.);
1454 tpcrow->SetFastCluster(zi,i); // write index
1457 for (Int_t i=0;i<510;i++){
1458 if (tpcrow->GetFastCluster(i)<0)
1459 tpcrow->SetFastCluster(i,last);
1461 last = tpcrow->GetFastCluster(i);
1473 //_________________________________________________________________________
1474 AliTPCclusterMI *AliTPCtrackerMI::GetClusterMI(Int_t index) const {
1475 //--------------------------------------------------------------------
1476 // Return pointer to a given cluster
1477 //--------------------------------------------------------------------
1478 if (index<0) return 0; // no cluster
1479 Int_t sec=(index&0xff000000)>>24;
1480 Int_t row=(index&0x00ff0000)>>16;
1481 Int_t ncl=(index&0x00007fff)>>00;
1483 const AliTPCtrackerRow * tpcrow=0;
1484 AliTPCclusterMI * clrow =0;
1486 if (sec<0 || sec>=fkNIS*4) {
1487 AliWarning(Form("Wrong sector %d",sec));
1492 AliTPCtrackerSector& tracksec = fInnerSec[sec%fkNIS];
1493 if (tracksec.GetNRows()<=row) return 0;
1494 tpcrow = &(tracksec[row]);
1495 if (tpcrow==0) return 0;
1498 if (tpcrow->GetN1()<=ncl) return 0;
1499 clrow = tpcrow->GetClusters1();
1502 if (tpcrow->GetN2()<=ncl) return 0;
1503 clrow = tpcrow->GetClusters2();
1507 AliTPCtrackerSector& tracksec = fOuterSec[(sec-fkNIS*2)%fkNOS];
1508 if (tracksec.GetNRows()<=row) return 0;
1509 tpcrow = &(tracksec[row]);
1510 if (tpcrow==0) return 0;
1512 if (sec-2*fkNIS<fkNOS) {
1513 if (tpcrow->GetN1()<=ncl) return 0;
1514 clrow = tpcrow->GetClusters1();
1517 if (tpcrow->GetN2()<=ncl) return 0;
1518 clrow = tpcrow->GetClusters2();
1522 return &(clrow[ncl]);
1528 Int_t AliTPCtrackerMI::FollowToNext(AliTPCseed& t, Int_t nr) {
1529 //-----------------------------------------------------------------
1530 // This function tries to find a track prolongation to next pad row
1531 //-----------------------------------------------------------------
1533 Double_t x= GetXrow(nr), ymax=GetMaxY(nr);
1536 AliTPCclusterMI *cl=0;
1537 Int_t tpcindex= t.GetClusterIndex2(nr);
1539 // update current shape info every 5 pad-row
1540 // if ( (nr%5==0) || t.GetNumberOfClusters()<2 || (t.fCurrentSigmaY2<0.0001) ){
1544 if (fIteration>0 && tpcindex>=-1){ //if we have already clusters
1546 if (tpcindex==-1) return 0; //track in dead zone
1548 cl = t.GetClusterPointer(nr);
1549 if ( (cl==0) ) cl = GetClusterMI(tpcindex);
1550 t.SetCurrentClusterIndex1(tpcindex);
1553 Int_t relativesector = ((tpcindex&0xff000000)>>24)%18; // if previously accepted cluster in different sector
1554 Float_t angle = relativesector*fSectors->GetAlpha()+fSectors->GetAlphaShift();
1556 if (angle<-TMath::Pi()) angle += 2*TMath::Pi();
1557 if (angle>=TMath::Pi()) angle -= 2*TMath::Pi();
1559 if (TMath::Abs(angle-t.GetAlpha())>0.001){
1560 Double_t rotation = angle-t.GetAlpha();
1561 t.SetRelativeSector(relativesector);
1562 if (!t.Rotate(rotation)) return 0;
1564 if (!t.PropagateTo(x)) return 0;
1566 t.SetCurrentCluster(cl);
1568 Int_t accept = AcceptCluster(&t,t.GetCurrentCluster());
1569 if ((tpcindex&0x8000)==0) accept =0;
1571 //if founded cluster is acceptible
1572 if (cl->IsUsed(11)) { // id cluster is shared inrease uncertainty
1573 t.SetErrorY2(t.GetErrorY2()+0.03);
1574 t.SetErrorZ2(t.GetErrorZ2()+0.03);
1575 t.SetErrorY2(t.GetErrorY2()*3);
1576 t.SetErrorZ2(t.GetErrorZ2()*3);
1578 t.SetNFoundable(t.GetNFoundable()+1);
1579 UpdateTrack(&t,accept);
1584 if (TMath::Abs(t.GetSnp())>AliTPCReconstructor::GetMaxSnpTracker()) return 0; // cut on angle
1585 if (fIteration>1 && IsFindable(t)){
1586 // not look for new cluster during refitting
1587 t.SetNFoundable(t.GetNFoundable()+1);
1592 // if (TMath::Abs(t.GetSnp())>0.95 || TMath::Abs(x*t.GetC()-t.GetEta())>0.95) return 0;// patch 28 fev 06
1593 Double_t y=t.GetYat(x);
1594 if (TMath::Abs(y)>ymax){
1596 t.SetRelativeSector((t.GetRelativeSector()+1) % fN);
1597 if (!t.Rotate(fSectors->GetAlpha()))
1599 } else if (y <-ymax) {
1600 t.SetRelativeSector((t.GetRelativeSector()-1+fN) % fN);
1601 if (!t.Rotate(-fSectors->GetAlpha()))
1607 if (!t.PropagateTo(x)) {
1608 if (fIteration==0) t.SetRemoval(10);
1612 Double_t z=t.GetZ();
1615 if (!IsActive(t.GetRelativeSector(),nr)) {
1617 t.SetClusterIndex2(nr,-1);
1620 //AliInfo(Form("A - Sector%d phi %f - alpha %f", t.fRelativeSector,y/x, t.GetAlpha()));
1621 Bool_t isActive = IsActive(t.GetRelativeSector(),nr);
1622 Bool_t isActive2 = (nr>=fInnerSec->GetNRows()) ? fOuterSec[t.GetRelativeSector()][nr-fInnerSec->GetNRows()].GetN()>0:fInnerSec[t.GetRelativeSector()][nr].GetN()>0;
1624 if (!isActive || !isActive2) return 0;
1626 const AliTPCtrackerRow &krow=GetRow(t.GetRelativeSector(),nr);
1627 if ( (t.GetSigmaY2()<0) || t.GetSigmaZ2()<0) return 0;
1629 Double_t roadz = 1.;
1631 if (TMath::Abs(TMath::Abs(y)-ymax)<krow.GetDeadZone()){
1633 t.SetClusterIndex2(nr,-1);
1639 // if (TMath::Abs(z)<(AliTPCReconstructor::GetCtgRange()*x+10) && TMath::Abs(z)<fParam->GetZLength(0) && (TMath::Abs(t.GetSnp())<AliTPCReconstructor::GetMaxSnpTracker()))
1640 t.SetNFoundable(t.GetNFoundable()+1);
1646 // cl = krow.FindNearest2(y+10.,z,roady,roadz,index);
1647 cl = krow.FindNearest2(y,z,roady,roadz,index);
1648 if (cl) t.SetCurrentClusterIndex1(krow.GetIndex(index));
1651 t.SetCurrentCluster(cl);
1653 if (fIteration==2&&cl->IsUsed(10)) return 0;
1654 Int_t accept = AcceptCluster(&t,t.GetCurrentCluster());
1655 if (fIteration==2&&cl->IsUsed(11)) {
1656 t.SetErrorY2(t.GetErrorY2()+0.03);
1657 t.SetErrorZ2(t.GetErrorZ2()+0.03);
1658 t.SetErrorY2(t.GetErrorY2()*3);
1659 t.SetErrorZ2(t.GetErrorZ2()*3);
1662 if (t.fCurrentCluster->IsUsed(10)){
1667 if (t.fNShared>0.7*t.GetNumberOfClusters()) {
1673 if (accept<3) UpdateTrack(&t,accept);
1676 if ( fIteration==0 && t.GetNFoundable()*0.5 > t.GetNumberOfClusters()) t.SetRemoval(10);
1684 //_________________________________________________________________________
1685 Bool_t AliTPCtrackerMI::GetTrackPoint(Int_t index, AliTrackPoint &p ) const
1687 // Get track space point by index
1688 // return false in case the cluster doesn't exist
1689 AliTPCclusterMI *cl = GetClusterMI(index);
1690 if (!cl) return kFALSE;
1691 Int_t sector = (index&0xff000000)>>24;
1692 // Int_t row = (index&0x00ff0000)>>16;
1694 // xyz[0] = fParam->GetPadRowRadii(sector,row);
1695 xyz[0] = cl->GetX();
1696 xyz[1] = cl->GetY();
1697 xyz[2] = cl->GetZ();
1699 fParam->AdjustCosSin(sector,cos,sin);
1700 Float_t x = cos*xyz[0]-sin*xyz[1];
1701 Float_t y = cos*xyz[1]+sin*xyz[0];
1703 Float_t sigmaY2 = 0.027*cl->GetSigmaY2();
1704 if (sector < fParam->GetNInnerSector()) sigmaY2 *= 2.07;
1705 Float_t sigmaZ2 = 0.066*cl->GetSigmaZ2();
1706 if (sector < fParam->GetNInnerSector()) sigmaZ2 *= 1.77;
1707 cov[0] = sin*sin*sigmaY2;
1708 cov[1] = -sin*cos*sigmaY2;
1710 cov[3] = cos*cos*sigmaY2;
1713 p.SetXYZ(x,y,xyz[2],cov);
1714 AliGeomManager::ELayerID iLayer;
1716 if (sector < fParam->GetNInnerSector()) {
1717 iLayer = AliGeomManager::kTPC1;
1721 iLayer = AliGeomManager::kTPC2;
1722 idet = sector - fParam->GetNInnerSector();
1724 UShort_t volid = AliGeomManager::LayerToVolUID(iLayer,idet);
1725 p.SetVolumeID(volid);
1731 Int_t AliTPCtrackerMI::UpdateClusters(AliTPCseed& t, Int_t nr) {
1732 //-----------------------------------------------------------------
1733 // This function tries to find a track prolongation to next pad row
1734 //-----------------------------------------------------------------
1735 t.SetCurrentCluster(0);
1736 t.SetCurrentClusterIndex1(0);
1738 Double_t xt=t.GetX();
1739 Int_t row = GetRowNumber(xt)-1;
1740 Double_t ymax= GetMaxY(nr);
1742 if (row < nr) return 1; // don't prolongate if not information until now -
1743 // if (TMath::Abs(t.GetSnp())>0.9 && t.GetNumberOfClusters()>40. && fIteration!=2) {
1745 // return 0; // not prolongate strongly inclined tracks
1747 // if (TMath::Abs(t.GetSnp())>0.95) {
1749 // return 0; // not prolongate strongly inclined tracks
1750 // }// patch 28 fev 06
1752 Double_t x= GetXrow(nr);
1754 //t.PropagateTo(x+0.02);
1755 //t.PropagateTo(x+0.01);
1756 if (!t.PropagateTo(x)){
1763 if (TMath::Abs(y)>ymax){
1765 t.SetRelativeSector((t.GetRelativeSector()+1) % fN);
1766 if (!t.Rotate(fSectors->GetAlpha()))
1768 } else if (y <-ymax) {
1769 t.SetRelativeSector((t.GetRelativeSector()-1+fN) % fN);
1770 if (!t.Rotate(-fSectors->GetAlpha()))
1773 // if (!t.PropagateTo(x)){
1780 if (TMath::Abs(t.GetSnp())>AliTPCReconstructor::GetMaxSnpTracker()) return 0;
1782 if (!IsActive(t.GetRelativeSector(),nr)) {
1784 t.SetClusterIndex2(nr,-1);
1787 //AliInfo(Form("A - Sector%d phi %f - alpha %f", t.fRelativeSector,y/x, t.GetAlpha()));
1789 AliTPCtrackerRow &krow=GetRow(t.GetRelativeSector(),nr);
1791 if (TMath::Abs(TMath::Abs(y)-ymax)<krow.GetDeadZone()){
1793 t.SetClusterIndex2(nr,-1);
1799 // if (TMath::Abs(t.GetZ())<(AliTPCReconstructor::GetCtgRange()*t.GetX()+10) && (TMath::Abs(t.GetSnp())<AliTPCReconstructor::GetMaxSnpTracker()))
1800 if (IsFindable(t)) t.SetNFoundable(t.GetNFoundable()+1);
1806 if ( (nr%2==0) || t.GetNumberOfClusters()<2){
1807 // t.fCurrentSigmaY = GetSigmaY(&t);
1808 //t.fCurrentSigmaZ = GetSigmaZ(&t);
1812 AliTPCclusterMI *cl=0;
1815 Double_t roady = 1.;
1816 Double_t roadz = 1.;
1820 index = t.GetClusterIndex2(nr);
1821 if ( (index>0) && (index&0x8000)==0){
1822 cl = t.GetClusterPointer(nr);
1823 if ( (cl==0) && (index>0)) cl = GetClusterMI(index);
1824 t.SetCurrentClusterIndex1(index);
1826 t.SetCurrentCluster(cl);
1832 // if (index<0) return 0;
1833 UInt_t uindex = TMath::Abs(index);
1836 //cl = krow.FindNearest2(y+10,z,roady,roadz,uindex);
1837 cl = krow.FindNearest2(y,z,roady,roadz,uindex);
1840 if (cl) t.SetCurrentClusterIndex1(krow.GetIndex(uindex));
1841 t.SetCurrentCluster(cl);
1847 Int_t AliTPCtrackerMI::FollowToNextCluster(AliTPCseed & t, Int_t nr) {
1848 //-----------------------------------------------------------------
1849 // This function tries to find a track prolongation to next pad row
1850 //-----------------------------------------------------------------
1852 //update error according neighborhoud
1854 if (t.GetCurrentCluster()) {
1856 Int_t accept = AcceptCluster(&t,t.GetCurrentCluster());
1858 if (t.GetCurrentCluster()->IsUsed(10)){
1863 t.SetNShared(t.GetNShared()+1);
1864 if (t.GetNShared()>0.7*t.GetNumberOfClusters()) {
1869 if (fIteration>0) accept = 0;
1870 if (accept<3) UpdateTrack(&t,accept);
1874 if ( ( (t.GetSigmaY2()+t.GetSigmaZ2())>0.16)&& t.GetNumberOfClusters()>18) t.SetRemoval(10);
1875 if ( t.GetChi2()/t.GetNumberOfClusters()>6 &&t.GetNumberOfClusters()>18) t.SetRemoval(10);
1877 if (( (t.GetNFoundable()*0.5 > t.GetNumberOfClusters()) || t.GetNoCluster()>15)) t.SetRemoval(10);
1885 //_____________________________________________________________________________
1886 Int_t AliTPCtrackerMI::FollowProlongation(AliTPCseed& t, Int_t rf, Int_t step) {
1887 //-----------------------------------------------------------------
1888 // This function tries to find a track prolongation.
1889 //-----------------------------------------------------------------
1890 Double_t xt=t.GetX();
1892 Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
1893 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
1894 if (alpha < 0. ) alpha += 2.*TMath::Pi();
1896 t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN);
1898 Int_t first = GetRowNumber(xt)-1;
1899 for (Int_t nr= first; nr>=rf; nr-=step) {
1901 if (t.GetKinkIndexes()[0]>0){
1902 for (Int_t i=0;i<3;i++){
1903 Int_t index = t.GetKinkIndexes()[i];
1904 if (index==0) break;
1905 if (index<0) continue;
1907 AliKink * kink = (AliKink*)fEvent->GetKink(index-1);
1909 printf("PROBLEM\n");
1912 Int_t kinkrow = kink->GetTPCRow0()+2+Int_t(0.5/(0.05+kink->GetAngle(2)));
1914 AliExternalTrackParam paramd(t);
1915 kink->SetDaughter(paramd);
1916 kink->SetStatus(2,5);
1923 if (nr==80) t.UpdateReference();
1924 if (nr<fInnerSec->GetNRows())
1925 fSectors = fInnerSec;
1927 fSectors = fOuterSec;
1928 if (FollowToNext(t,nr)==0)
1941 Int_t AliTPCtrackerMI::FollowBackProlongation(AliTPCseed& t, Int_t rf) {
1942 //-----------------------------------------------------------------
1943 // This function tries to find a track prolongation.
1944 //-----------------------------------------------------------------
1946 Double_t xt=t.GetX();
1947 Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
1948 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
1949 if (alpha < 0. ) alpha += 2.*TMath::Pi();
1950 t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN);
1952 Int_t first = t.GetFirstPoint();
1953 if (first<GetRowNumber(xt)+1) first = GetRowNumber(xt)+1;
1955 if (first<0) first=0;
1956 for (Int_t nr=first; nr<=rf; nr++) {
1957 // if ( (TMath::Abs(t.GetSnp())>0.95)) break;//patch 28 fev 06
1958 if (t.GetKinkIndexes()[0]<0){
1959 for (Int_t i=0;i<3;i++){
1960 Int_t index = t.GetKinkIndexes()[i];
1961 if (index==0) break;
1962 if (index>0) continue;
1963 index = TMath::Abs(index);
1964 AliKink * kink = (AliKink*)fEvent->GetKink(index-1);
1966 printf("PROBLEM\n");
1969 Int_t kinkrow = kink->GetTPCRow0()-2-Int_t(0.5/(0.05+kink->GetAngle(2)));
1971 AliExternalTrackParam paramm(t);
1972 kink->SetMother(paramm);
1973 kink->SetStatus(2,1);
1980 if (nr<fInnerSec->GetNRows())
1981 fSectors = fInnerSec;
1983 fSectors = fOuterSec;
1994 Float_t AliTPCtrackerMI::OverlapFactor(AliTPCseed * s1, AliTPCseed * s2, Int_t &sum1, Int_t & sum2)
2002 Float_t dz2 =(s1->GetZ() - s2->GetZ());
2005 Float_t dy2 =TMath::Abs((s1->GetY() - s2->GetY()));
2007 Float_t distance = TMath::Sqrt(dz2+dy2);
2008 if (distance>4.) return 0; // if there are far away - not overlap - to reduce combinatorics
2011 Int_t firstpoint = TMath::Min(s1->GetFirstPoint(),s2->GetFirstPoint());
2012 Int_t lastpoint = TMath::Max(s1->GetLastPoint(),s2->GetLastPoint());
2017 if (firstpoint>lastpoint) {
2018 firstpoint =lastpoint;
2023 for (Int_t i=firstpoint-1;i<lastpoint+1;i++){
2024 if (s1->GetClusterIndex2(i)>0) sum1++;
2025 if (s2->GetClusterIndex2(i)>0) sum2++;
2026 if (s1->GetClusterIndex2(i)==s2->GetClusterIndex2(i) && s1->GetClusterIndex2(i)>0) {
2030 if (sum<5) return 0;
2032 Float_t summin = TMath::Min(sum1+1,sum2+1);
2033 Float_t ratio = (sum+1)/Float_t(summin);
2037 void AliTPCtrackerMI::SignShared(AliTPCseed * s1, AliTPCseed * s2)
2041 Float_t thetaCut = 0.2;//+10.*TMath::Sqrt(s1->GetSigmaTglZ()+ s2->GetSigmaTglZ());
2042 if (TMath::Abs(s1->GetTgl()-s2->GetTgl())>thetaCut) return;
2043 Float_t minCl = TMath::Min(s1->GetNumberOfClusters(),s2->GetNumberOfClusters());
2044 Int_t cutN0 = TMath::Max(5,TMath::Nint(0.1*minCl));
2049 //Int_t firstpoint = TMath::Max(s1->GetFirstPoint(),s2->GetFirstPoint());
2050 //Int_t lastpoint = TMath::Min(s1->GetLastPoint(),s2->GetLastPoint());
2051 Int_t firstpoint = 0;
2052 Int_t lastpoint = 160;
2054 // if (firstpoint>=lastpoint-5) return;;
2056 for (Int_t i=firstpoint;i<lastpoint;i++){
2057 // if ( (s1->GetClusterIndex2(i)&0xFFFF8FFF)==(s2->GetClusterIndex2(i)&0xFFFF8FFF) && s1->GetClusterIndex2(i)>0) {
2058 if ( (s1->GetClusterIndex2(i))==(s2->GetClusterIndex2(i)) && s1->GetClusterIndex2(i)>0) {
2060 s1->SetSharedMapBit(i, kTRUE);
2061 s2->SetSharedMapBit(i, kTRUE);
2063 if (s1->GetClusterIndex2(i)>0)
2064 s1->SetClusterMapBit(i, kTRUE);
2066 if (sumshared>cutN0){
2069 for (Int_t i=firstpoint;i<lastpoint;i++){
2070 // if ( (s1->GetClusterIndex2(i)&0xFFFF8FFF)==(s2->GetClusterIndex2(i)&0xFFFF8FFF) && s1->GetClusterIndex2(i)>0) {
2071 if ( (s1->GetClusterIndex2(i))==(s2->GetClusterIndex2(i)) && s1->GetClusterIndex2(i)>0) {
2072 AliTPCTrackerPoint *p1 = s1->GetTrackPoint(i);
2073 AliTPCTrackerPoint *p2 = s2->GetTrackPoint(i);;
2074 if (s1->IsActive()&&s2->IsActive()){
2075 p1->SetShared(kTRUE);
2076 p2->SetShared(kTRUE);
2082 if (sumshared>cutN0){
2083 for (Int_t i=0;i<4;i++){
2084 if (s1->GetOverlapLabel(3*i)==0){
2085 s1->SetOverlapLabel(3*i, s2->GetLabel());
2086 s1->SetOverlapLabel(3*i+1,sumshared);
2087 s1->SetOverlapLabel(3*i+2,s2->GetUniqueID());
2091 for (Int_t i=0;i<4;i++){
2092 if (s2->GetOverlapLabel(3*i)==0){
2093 s2->SetOverlapLabel(3*i, s1->GetLabel());
2094 s2->SetOverlapLabel(3*i+1,sumshared);
2095 s2->SetOverlapLabel(3*i+2,s1->GetUniqueID());
2102 void AliTPCtrackerMI::SignShared(TObjArray * arr)
2105 //sort trackss according sectors
2107 for (Int_t i=0; i<arr->GetEntriesFast(); i++) {
2108 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2110 //if (pt) RotateToLocal(pt);
2114 arr->Sort(); // sorting according relative sectors
2115 arr->Expand(arr->GetEntries());
2118 Int_t nseed=arr->GetEntriesFast();
2119 for (Int_t i=0; i<nseed; i++) {
2120 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2122 for (Int_t j=0;j<=12;j++){
2123 pt->SetOverlapLabel(j,0);
2126 for (Int_t i=0; i<nseed; i++) {
2127 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2129 if (pt->GetRemoval()>10) continue;
2130 for (Int_t j=i+1; j<nseed; j++){
2131 AliTPCseed *pt2=(AliTPCseed*)arr->UncheckedAt(j);
2132 if (TMath::Abs(pt->GetRelativeSector()-pt2->GetRelativeSector())>1) continue;
2134 if (pt2->GetRemoval()<=10) {
2135 //if ( TMath::Abs(pt->GetRelativeSector()-pt2->GetRelativeSector())>0) break;
2143 void AliTPCtrackerMI::SortTracks(TObjArray * arr, Int_t mode) const
2146 //sort tracks in array according mode criteria
2147 Int_t nseed = arr->GetEntriesFast();
2148 for (Int_t i=0; i<nseed; i++) {
2149 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2160 void AliTPCtrackerMI::RemoveUsed2(TObjArray * arr, Float_t factor1, Float_t factor2, Int_t minimal)
2163 // Loop over all tracks and remove overlaped tracks (with lower quality)
2165 // 1. Unsign clusters
2166 // 2. Sort tracks according quality
2167 // Quality is defined by the number of cluster between first and last points
2169 // 3. Loop over tracks - decreasing quality order
2170 // a.) remove - If the fraction of shared cluster less than factor (1- n or 2)
2171 // b.) remove - If the minimal number of clusters less than minimal and not ITS
2172 // c.) if track accepted - sign clusters
2174 //Called in - AliTPCtrackerMI::Clusters2Tracks()
2175 // - AliTPCtrackerMI::PropagateBack()
2176 // - AliTPCtrackerMI::RefitInward()
2179 // factor1 - factor for constrained
2180 // factor2 - for non constrained tracks
2181 // if (Float_t(shared+1)/Float_t(found+1)>factor) - DELETE
2185 Int_t nseed = arr->GetEntriesFast();
2186 Float_t * quality = new Float_t[nseed];
2187 Int_t * indexes = new Int_t[nseed];
2191 for (Int_t i=0; i<nseed; i++) {
2192 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2197 pt->UpdatePoints(); //select first last max dens points
2198 Float_t * points = pt->GetPoints();
2199 if (points[3]<0.8) quality[i] =-1;
2200 quality[i] = (points[2]-points[0])+pt->GetNumberOfClusters();
2201 //prefer high momenta tracks if overlaps
2202 quality[i] *= TMath::Sqrt(TMath::Abs(pt->Pt())+0.5);
2204 TMath::Sort(nseed,quality,indexes);
2207 for (Int_t itrack=0; itrack<nseed; itrack++) {
2208 Int_t trackindex = indexes[itrack];
2209 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(trackindex);
2212 if (quality[trackindex]<0){
2214 delete arr->RemoveAt(trackindex);
2217 arr->RemoveAt(trackindex);
2223 Int_t first = Int_t(pt->GetPoints()[0]);
2224 Int_t last = Int_t(pt->GetPoints()[2]);
2225 Double_t factor = (pt->GetBConstrain()) ? factor1: factor2;
2227 Int_t found,foundable,shared;
2228 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
2229 // pt->GetClusterStatistic(0,160, found, foundable,shared,kFALSE);
2230 Bool_t itsgold =kFALSE;
2233 if (pt->GetESD()->GetITSclusters(dummy)>4) itsgold= kTRUE;
2237 if (Float_t(shared+1)/Float_t(found+1)>factor){
2238 if (pt->GetKinkIndexes()[0]!=0) continue; //don't remove tracks - part of the kinks
2239 if( AliTPCReconstructor::StreamLevel()>15){
2240 TTreeSRedirector &cstream = *fDebugStreamer;
2241 cstream<<"RemoveUsed"<<
2242 "iter="<<fIteration<<
2246 delete arr->RemoveAt(trackindex);
2249 if (pt->GetNumberOfClusters()<50&&(found-0.5*shared)<minimal){ //remove short tracks
2250 if (pt->GetKinkIndexes()[0]!=0) continue; //don't remove tracks - part of the kinks
2251 if( AliTPCReconstructor::StreamLevel()>15){
2252 TTreeSRedirector &cstream = *fDebugStreamer;
2253 cstream<<"RemoveShort"<<
2254 "iter="<<fIteration<<
2258 delete arr->RemoveAt(trackindex);
2264 //if (sharedfactor>0.4) continue;
2265 if (pt->GetKinkIndexes()[0]>0) continue;
2266 //Remove tracks with undefined properties - seems
2267 if (pt->GetSigmaY2()<kAlmost0) continue; // ? what is the origin ?
2269 for (Int_t i=first; i<last; i++) {
2270 Int_t index=pt->GetClusterIndex2(i);
2271 // if (index<0 || index&0x8000 ) continue;
2272 if (index<0 || index&0x8000 ) continue;
2273 AliTPCclusterMI *c= pt->GetClusterPointer(i);
2280 Info("RemoveUsed","\n*****\nNumber of good tracks after shared removal\t%d\n",fNtracks);
2286 void AliTPCtrackerMI::UnsignClusters()
2289 // loop over all clusters and unsign them
2292 for (Int_t sec=0;sec<fkNIS;sec++){
2293 for (Int_t row=0;row<fInnerSec->GetNRows();row++){
2294 AliTPCclusterMI *cl = fInnerSec[sec][row].GetClusters1();
2295 for (Int_t icl =0;icl< fInnerSec[sec][row].GetN1();icl++)
2296 // if (cl[icl].IsUsed(10))
2298 cl = fInnerSec[sec][row].GetClusters2();
2299 for (Int_t icl =0;icl< fInnerSec[sec][row].GetN2();icl++)
2300 //if (cl[icl].IsUsed(10))
2305 for (Int_t sec=0;sec<fkNOS;sec++){
2306 for (Int_t row=0;row<fOuterSec->GetNRows();row++){
2307 AliTPCclusterMI *cl = fOuterSec[sec][row].GetClusters1();
2308 for (Int_t icl =0;icl< fOuterSec[sec][row].GetN1();icl++)
2309 //if (cl[icl].IsUsed(10))
2311 cl = fOuterSec[sec][row].GetClusters2();
2312 for (Int_t icl =0;icl< fOuterSec[sec][row].GetN2();icl++)
2313 //if (cl[icl].IsUsed(10))
2322 void AliTPCtrackerMI::SignClusters(TObjArray * arr, Float_t fnumber, Float_t fdensity)
2325 //sign clusters to be "used"
2327 // snumber and sdensity sign number of sigmas - bellow mean value to be accepted
2328 // loop over "primaries"
2342 Int_t nseed = arr->GetEntriesFast();
2343 for (Int_t i=0; i<nseed; i++) {
2344 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2348 if (!(pt->IsActive())) continue;
2349 Float_t dens = pt->GetNumberOfClusters()/Float_t(pt->GetNFoundable());
2350 if ( (dens>0.7) && (pt->GetNumberOfClusters()>70)){
2352 sumdens2+= dens*dens;
2353 sumn += pt->GetNumberOfClusters();
2354 sumn2 += pt->GetNumberOfClusters()*pt->GetNumberOfClusters();
2355 Float_t chi2 = pt->GetChi2()/pt->GetNumberOfClusters();
2358 sumchi2 +=chi2*chi2;
2363 Float_t mdensity = 0.9;
2364 Float_t meann = 130;
2365 Float_t meanchi = 1;
2366 Float_t sdensity = 0.1;
2367 Float_t smeann = 10;
2368 Float_t smeanchi =0.4;
2372 mdensity = sumdens/sum;
2374 meanchi = sumchi/sum;
2376 sdensity = sumdens2/sum-mdensity*mdensity;
2378 sdensity = TMath::Sqrt(sdensity);
2382 smeann = sumn2/sum-meann*meann;
2384 smeann = TMath::Sqrt(smeann);
2388 smeanchi = sumchi2/sum - meanchi*meanchi;
2390 smeanchi = TMath::Sqrt(smeanchi);
2396 //REMOVE SHORT DELTAS or tracks going out of sensitive volume of TPC
2398 for (Int_t i=0; i<nseed; i++) {
2399 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2403 if (pt->GetBSigned()) continue;
2404 if (pt->GetBConstrain()) continue;
2405 //if (!(pt->IsActive())) continue;
2407 Int_t found,foundable,shared;
2408 pt->GetClusterStatistic(0,160,found, foundable,shared);
2409 if (shared/float(found)>0.3) {
2410 if (shared/float(found)>0.9 ){
2411 //delete arr->RemoveAt(i);
2416 Bool_t isok =kFALSE;
2417 if ( (pt->GetNShared()/pt->GetNumberOfClusters()<0.5) &&pt->GetNumberOfClusters()>60)
2419 if ((TMath::Abs(1/pt->GetC())<100.) && (pt->GetNShared()/pt->GetNumberOfClusters()<0.7))
2421 if (TMath::Abs(pt->GetZ()/pt->GetX())>1.1)
2423 if ( (TMath::Abs(pt->GetSnp()>0.7) && pt->GetD(0,0)>60.))
2427 for (Int_t j=0; j<160; ++j) {
2428 Int_t index=pt->GetClusterIndex2(j);
2429 if (index<0) continue;
2430 AliTPCclusterMI *c= pt->GetClusterPointer(j);
2432 //if (!(c->IsUsed(10))) c->Use();
2439 Double_t maxchi = meanchi+2.*smeanchi;
2441 for (Int_t i=0; i<nseed; i++) {
2442 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2446 //if (!(pt->IsActive())) continue;
2447 if (pt->GetBSigned()) continue;
2448 Double_t chi = pt->GetChi2()/pt->GetNumberOfClusters();
2449 if (chi>maxchi) continue;
2452 Float_t dens = pt->GetNumberOfClusters()/Float_t(pt->GetNFoundable());
2454 //sign only tracks with enoug big density at the beginning
2456 if ((pt->GetDensityFirst(40)<0.75) && pt->GetNumberOfClusters()<meann) continue;
2459 Double_t mindens = TMath::Max(double(mdensity-sdensity*fdensity*bfactor),0.65);
2460 Double_t minn = TMath::Max(Int_t(meann-fnumber*smeann*bfactor),50);
2462 // if (pt->fBConstrain) mindens = TMath::Max(mdensity-sdensity*fdensity*bfactor,0.65);
2463 if ( (pt->GetRemoval()==10) && (pt->GetSnp()>0.8)&&(dens>mindens))
2466 if ((dens>mindens && pt->GetNumberOfClusters()>minn) && chi<maxchi ){
2467 //Int_t noc=pt->GetNumberOfClusters();
2468 pt->SetBSigned(kTRUE);
2469 for (Int_t j=0; j<160; ++j) {
2471 Int_t index=pt->GetClusterIndex2(j);
2472 if (index<0) continue;
2473 AliTPCclusterMI *c= pt->GetClusterPointer(j);
2475 // if (!(c->IsUsed(10))) c->Use();
2480 // gLastCheck = nseed;
2488 void AliTPCtrackerMI::StopNotActive(TObjArray * arr, Int_t row0, Float_t th0, Float_t th1, Float_t th2) const
2490 // stop not active tracks
2491 // take th1 as threshold for number of founded to number of foundable on last 10 active rows
2492 // take th2 as threshold for number of founded to number of foundable on last 20 active rows
2493 Int_t nseed = arr->GetEntriesFast();
2495 for (Int_t i=0; i<nseed; i++) {
2496 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2500 if (!(pt->IsActive())) continue;
2501 StopNotActive(pt,row0,th0, th1,th2);
2507 void AliTPCtrackerMI::StopNotActive(AliTPCseed * seed, Int_t row0, Float_t th0, Float_t th1,
2510 // stop not active tracks
2511 // take th1 as threshold for number of founded to number of foundable on last 10 active rows
2512 // take th2 as threshold for number of founded to number of foundable on last 20 active rows
2515 Int_t foundable = 0;
2516 Int_t maxindex = seed->GetLastPoint(); //last foundable row
2517 if (seed->GetNFoundable()*th0 > seed->GetNumberOfClusters()) {
2518 seed->Desactivate(10) ;
2522 for (Int_t i=row0; i<maxindex; i++){
2523 Int_t index = seed->GetClusterIndex2(i);
2524 if (index!=-1) foundable++;
2526 if (foundable<=30) sumgood1++;
2527 if (foundable<=50) {
2534 if (foundable>=30.){
2535 if (sumgood1<(th1*30.)) seed->Desactivate(10);
2538 if (sumgood2<(th2*50.)) seed->Desactivate(10);
2542 Int_t AliTPCtrackerMI::RefitInward(AliESDEvent *event)
2545 // back propagation of ESD tracks
2548 const Int_t kMaxFriendTracks=2000;
2552 //PrepareForProlongation(fSeeds,1);
2553 PropagateForward2(fSeeds);
2554 RemoveUsed2(fSeeds,0.4,0.4,20);
2556 TObjArray arraySeed(fSeeds->GetEntries());
2557 for (Int_t i=0;i<fSeeds->GetEntries();i++) {
2558 arraySeed.AddAt(fSeeds->At(i),i);
2560 SignShared(&arraySeed);
2561 // FindCurling(fSeeds, event,2); // find multi found tracks
2562 FindSplitted(fSeeds, event,2); // find multi found tracks
2563 if (AliTPCReconstructor::StreamLevel()>2) FindMultiMC(fSeeds, fEvent,2); // find multi found tracks
2566 Int_t nseed = fSeeds->GetEntriesFast();
2567 for (Int_t i=0;i<nseed;i++){
2568 AliTPCseed * seed = (AliTPCseed*) fSeeds->UncheckedAt(i);
2569 if (!seed) continue;
2570 if (seed->GetKinkIndex(0)>0) UpdateKinkQualityD(seed); // update quality informations for kinks
2571 AliESDtrack *esd=event->GetTrack(i);
2573 if (seed->GetNumberOfClusters()<60 && seed->GetNumberOfClusters()<(esd->GetTPCclusters(0) -5)*0.8){
2574 AliExternalTrackParam paramIn;
2575 AliExternalTrackParam paramOut;
2576 Int_t ncl = seed->RefitTrack(seed,¶mIn,¶mOut);
2577 if (AliTPCReconstructor::StreamLevel()>0) {
2578 (*fDebugStreamer)<<"RecoverIn"<<
2582 "pout.="<<¶mOut<<
2587 seed->Set(paramIn.GetX(),paramIn.GetAlpha(),paramIn.GetParameter(),paramIn.GetCovariance());
2588 seed->SetNumberOfClusters(ncl);
2592 seed->PropagateTo(fParam->GetInnerRadiusLow());
2593 seed->UpdatePoints();
2594 AddCovariance(seed);
2596 seed->CookdEdx(0.02,0.6);
2597 CookLabel(seed,0.1); //For comparison only
2598 esd->SetTPCClusterMap(seed->GetClusterMap());
2599 esd->SetTPCSharedMap(seed->GetSharedMap());
2601 if (AliTPCReconstructor::StreamLevel()>1 && seed!=0&&esd!=0) {
2602 TTreeSRedirector &cstream = *fDebugStreamer;
2609 if (seed->GetNumberOfClusters()>15){
2610 esd->UpdateTrackParams(seed,AliESDtrack::kTPCrefit);
2611 esd->SetTPCPoints(seed->GetPoints());
2612 esd->SetTPCPointsF(seed->GetNFoundable());
2613 Int_t ndedx = seed->GetNCDEDX(0)+seed->GetNCDEDX(1)+seed->GetNCDEDX(2)+seed->GetNCDEDX(3);
2614 Float_t sdedx = (seed->GetSDEDX(0)+seed->GetSDEDX(1)+seed->GetSDEDX(2)+seed->GetSDEDX(3))*0.25;
2615 Float_t dedx = seed->GetdEdx();
2616 esd->SetTPCsignal(dedx, sdedx, ndedx);
2618 // add seed to the esd track in Calib level
2620 Bool_t storeFriend = gRandom->Rndm()<(kMaxFriendTracks)/Float_t(nseed);
2621 if (AliTPCReconstructor::StreamLevel()>0 &&storeFriend){
2622 AliTPCseed * seedCopy = new AliTPCseed(*seed, kTRUE);
2623 esd->AddCalibObject(seedCopy);
2628 //printf("problem\n");
2631 //FindKinks(fSeeds,event);
2632 Info("RefitInward","Number of refitted tracks %d",ntracks);
2637 Int_t AliTPCtrackerMI::PropagateBack(AliESDEvent *event)
2640 // back propagation of ESD tracks
2646 PropagateBack(fSeeds);
2647 RemoveUsed2(fSeeds,0.4,0.4,20);
2648 //FindCurling(fSeeds, fEvent,1);
2649 FindSplitted(fSeeds, event,1); // find multi found tracks
2650 if (AliTPCReconstructor::StreamLevel()>2) FindMultiMC(fSeeds, fEvent,1); // find multi found tracks
2653 Int_t nseed = fSeeds->GetEntriesFast();
2655 for (Int_t i=0;i<nseed;i++){
2656 AliTPCseed * seed = (AliTPCseed*) fSeeds->UncheckedAt(i);
2657 if (!seed) continue;
2658 if (seed->GetKinkIndex(0)<0) UpdateKinkQualityM(seed); // update quality informations for kinks
2659 seed->UpdatePoints();
2660 AddCovariance(seed);
2661 AliESDtrack *esd=event->GetTrack(i);
2662 if (seed->GetNumberOfClusters()<60 && seed->GetNumberOfClusters()<(esd->GetTPCclusters(0) -5)*0.8){
2663 AliExternalTrackParam paramIn;
2664 AliExternalTrackParam paramOut;
2665 Int_t ncl = seed->RefitTrack(seed,¶mIn,¶mOut);
2666 if (AliTPCReconstructor::StreamLevel()>0) {
2667 (*fDebugStreamer)<<"RecoverBack"<<
2671 "pout.="<<¶mOut<<
2676 seed->Set(paramOut.GetX(),paramOut.GetAlpha(),paramOut.GetParameter(),paramOut.GetCovariance());
2677 seed->SetNumberOfClusters(ncl);
2680 seed->CookdEdx(0.02,0.6);
2681 CookLabel(seed,0.1); //For comparison only
2682 if (seed->GetNumberOfClusters()>15){
2683 esd->UpdateTrackParams(seed,AliESDtrack::kTPCout);
2684 esd->SetTPCPoints(seed->GetPoints());
2685 esd->SetTPCPointsF(seed->GetNFoundable());
2686 Int_t ndedx = seed->GetNCDEDX(0)+seed->GetNCDEDX(1)+seed->GetNCDEDX(2)+seed->GetNCDEDX(3);
2687 Float_t sdedx = (seed->GetSDEDX(0)+seed->GetSDEDX(1)+seed->GetSDEDX(2)+seed->GetSDEDX(3))*0.25;
2688 Float_t dedx = seed->GetdEdx();
2689 esd->SetTPCsignal(dedx, sdedx, ndedx);
2691 Int_t eventnumber = event->GetEventNumberInFile();// patch 28 fev 06
2692 // This is most likely NOT the event number you'd like to use. It has nothing to do with the 'real' event number
2693 if (AliTPCReconstructor::StreamLevel()>1 && esd) {
2694 (*fDebugStreamer)<<"Cback"<<
2697 "EventNrInFile="<<eventnumber<<
2702 //FindKinks(fSeeds,event);
2703 Info("PropagateBack","Number of back propagated tracks %d",ntracks);
2710 void AliTPCtrackerMI::DeleteSeeds()
2715 Int_t nseed = fSeeds->GetEntriesFast();
2716 for (Int_t i=0;i<nseed;i++){
2717 AliTPCseed * seed = (AliTPCseed*)fSeeds->At(i);
2718 if (seed) delete fSeeds->RemoveAt(i);
2725 void AliTPCtrackerMI::ReadSeeds(AliESDEvent *event, Int_t direction)
2728 //read seeds from the event
2730 Int_t nentr=event->GetNumberOfTracks();
2732 Info("ReadSeeds", "Number of ESD tracks: %d\n", nentr);
2737 fSeeds = new TObjArray(nentr);
2741 for (Int_t i=0; i<nentr; i++) {
2742 AliESDtrack *esd=event->GetTrack(i);
2743 ULong_t status=esd->GetStatus();
2744 if (!(status&AliESDtrack::kTPCin)) continue;
2745 AliTPCtrack t(*esd);
2746 t.SetNumberOfClusters(0);
2747 // AliTPCseed *seed = new AliTPCseed(t,t.GetAlpha());
2748 AliTPCseed *seed = new AliTPCseed(t/*,t.GetAlpha()*/);
2749 seed->SetUniqueID(esd->GetID());
2750 AddCovariance(seed); //add systematic ucertainty
2751 for (Int_t ikink=0;ikink<3;ikink++) {
2752 Int_t index = esd->GetKinkIndex(ikink);
2753 seed->GetKinkIndexes()[ikink] = index;
2754 if (index==0) continue;
2755 index = TMath::Abs(index);
2756 AliESDkink * kink = fEvent->GetKink(index-1);
2757 if (kink&&esd->GetKinkIndex(ikink)<0){
2758 if ((status & AliESDtrack::kTRDrefit) != 0) kink->SetStatus(1,2);
2759 if ((status & AliESDtrack::kITSout) != 0) kink->SetStatus(1,0);
2761 if (kink&&esd->GetKinkIndex(ikink)>0){
2762 if ((status & AliESDtrack::kTRDrefit) != 0) kink->SetStatus(1,6);
2763 if ((status & AliESDtrack::kITSout) != 0) kink->SetStatus(1,4);
2767 if (((status&AliESDtrack::kITSout)==0)&&(direction==1)) seed->ResetCovariance(10.);
2768 if ( direction ==2 &&(status & AliESDtrack::kTRDrefit) == 0 ) seed->ResetCovariance(10.);
2769 //if ( direction ==2 && ((status & AliESDtrack::kTPCout) == 0) ) {
2770 // fSeeds->AddAt(0,i);
2774 if ( direction ==2 &&(status & AliESDtrack::kTRDrefit) > 0 ) {
2775 Double_t par0[5],par1[5],alpha,x;
2776 esd->GetInnerExternalParameters(alpha,x,par0);
2777 esd->GetExternalParameters(x,par1);
2778 Double_t delta1 = TMath::Abs(par0[4]-par1[4])/(0.000000001+TMath::Abs(par0[4]+par1[4]));
2779 Double_t delta2 = TMath::Abs(par0[3]-par1[3]);
2781 if (esd->GetTRDncls()>0) trdchi2 = esd->GetTRDchi2()/esd->GetTRDncls();
2782 //reset covariance if suspicious
2783 if ( (delta1>0.1) || (delta2>0.006) ||trdchi2>7.)
2784 seed->ResetCovariance(10.);
2789 // rotate to the local coordinate system
2791 fSectors=fInnerSec; fN=fkNIS;
2792 Double_t alpha=seed->GetAlpha() - fSectors->GetAlphaShift();
2793 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
2794 if (alpha < 0. ) alpha += 2.*TMath::Pi();
2795 Int_t ns=Int_t(alpha/fSectors->GetAlpha())%fN;
2796 alpha =ns*fSectors->GetAlpha() + fSectors->GetAlphaShift();
2797 if (alpha<-TMath::Pi()) alpha += 2*TMath::Pi();
2798 if (alpha>=TMath::Pi()) alpha -= 2*TMath::Pi();
2799 alpha-=seed->GetAlpha();
2800 if (!seed->Rotate(alpha)) {
2806 if (esd->GetKinkIndex(0)<=0){
2807 for (Int_t irow=0;irow<160;irow++){
2808 Int_t index = seed->GetClusterIndex2(irow);
2811 AliTPCclusterMI * cl = GetClusterMI(index);
2812 seed->SetClusterPointer(irow,cl);
2814 if ((index & 0x8000)==0){
2815 cl->Use(10); // accepted cluster
2817 cl->Use(6); // close cluster not accepted
2820 Info("ReadSeeds","Not found cluster");
2825 fSeeds->AddAt(seed,i);
2831 //_____________________________________________________________________________
2832 void AliTPCtrackerMI::MakeSeeds3(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, Float_t cuts[4],
2833 Float_t deltay, Int_t ddsec) {
2834 //-----------------------------------------------------------------
2835 // This function creates track seeds.
2836 // SEEDING WITH VERTEX CONSTRAIN
2837 //-----------------------------------------------------------------
2838 // cuts[0] - fP4 cut
2839 // cuts[1] - tan(phi) cut
2840 // cuts[2] - zvertex cut
2841 // cuts[3] - fP3 cut
2849 Double_t x[5], c[15];
2850 // Int_t di = i1-i2;
2852 Double_t alpha=fSectors->GetAlpha(), shift=fSectors->GetAlphaShift();
2853 Double_t cs=cos(alpha), sn=sin(alpha);
2855 Double_t x1 =GetXrow(i1);
2856 Double_t xx2=GetXrow(i2);
2858 Double_t x3=GetX(), y3=GetY(), z3=GetZ();
2860 Int_t imiddle = (i2+i1)/2; //middle pad row index
2861 Double_t xm = GetXrow(imiddle); // radius of middle pad-row
2862 const AliTPCtrackerRow& krm=GetRow(sec,imiddle); //middle pad -row
2866 const AliTPCtrackerRow& kr1=GetRow(ns,i1);
2867 Double_t ymax = GetMaxY(i1)-kr1.GetDeadZone()-1.5;
2868 Double_t ymaxm = GetMaxY(imiddle)-kr1.GetDeadZone()-1.5;
2871 // change cut on curvature if it can't reach this layer
2872 // maximal curvature set to reach it
2873 Double_t dvertexmax = TMath::Sqrt((x1-x3)*(x1-x3)+(ymax+5-y3)*(ymax+5-y3));
2874 if (dvertexmax*0.5*cuts[0]>0.85){
2875 cuts[0] = 0.85/(dvertexmax*0.5+1.);
2877 Double_t r2min = 1/(cuts[0]*cuts[0]); //minimal square of radius given by cut
2880 if (deltay>0) ddsec = 0;
2881 // loop over clusters
2882 for (Int_t is=0; is < kr1; is++) {
2884 if (kr1[is]->IsUsed(10)) continue;
2885 Double_t y1=kr1[is]->GetY(), z1=kr1[is]->GetZ();
2886 //if (TMath::Abs(y1)>ymax) continue;
2888 if (deltay>0 && TMath::Abs(ymax-TMath::Abs(y1))> deltay ) continue; // seed only at the edge
2890 // find possible directions
2891 Float_t anglez = (z1-z3)/(x1-x3);
2892 Float_t extraz = z1 - anglez*(x1-xx2); // extrapolated z
2895 //find rotation angles relative to line given by vertex and point 1
2896 Double_t dvertex2 = (x1-x3)*(x1-x3)+(y1-y3)*(y1-y3);
2897 Double_t dvertex = TMath::Sqrt(dvertex2);
2898 Double_t angle13 = TMath::ATan((y1-y3)/(x1-x3));
2899 Double_t cs13 = cos(-angle13), sn13 = sin(-angle13);
2902 // loop over 2 sectors
2908 Double_t dddz1=0; // direction of delta inclination in z axis
2915 for (Int_t dsec = dsec1; dsec<=dsec2;dsec++){
2916 Int_t sec2 = sec + dsec;
2918 // AliTPCtrackerRow& kr2 = fOuterSec[(sec2+fkNOS)%fkNOS][i2];
2919 //AliTPCtrackerRow& kr2m = fOuterSec[(sec2+fkNOS)%fkNOS][imiddle];
2920 AliTPCtrackerRow& kr2 = GetRow((sec2+fkNOS)%fkNOS,i2);
2921 AliTPCtrackerRow& kr2m = GetRow((sec2+fkNOS)%fkNOS,imiddle);
2922 Int_t index1 = TMath::Max(kr2.Find(extraz-0.6-dddz1*TMath::Abs(z1)*0.05)-1,0);
2923 Int_t index2 = TMath::Min(kr2.Find(extraz+0.6+dddz2*TMath::Abs(z1)*0.05)+1,kr2);
2925 // rotation angles to p1-p3
2926 Double_t cs13r = cos(-angle13+dsec*alpha)/dvertex, sn13r = sin(-angle13+dsec*alpha)/dvertex;
2927 Double_t x2, y2, z2;
2929 // Double_t dymax = maxangle*TMath::Abs(x1-xx2);
2932 Double_t dxx0 = (xx2-x3)*cs13r;
2933 Double_t dyy0 = (xx2-x3)*sn13r;
2934 for (Int_t js=index1; js < index2; js++) {
2935 const AliTPCclusterMI *kcl = kr2[js];
2936 if (kcl->IsUsed(10)) continue;
2938 //calcutate parameters
2940 Double_t yy0 = dyy0 +(kcl->GetY()-y3)*cs13r;
2942 if (TMath::Abs(yy0)<0.000001) continue;
2943 Double_t xx0 = dxx0 -(kcl->GetY()-y3)*sn13r;
2944 Double_t y0 = 0.5*(xx0*xx0+yy0*yy0-xx0)/yy0;
2945 Double_t r02 = (0.25+y0*y0)*dvertex2;
2946 //curvature (radius) cut
2947 if (r02<r2min) continue;
2951 Double_t c0 = 1/TMath::Sqrt(r02);
2955 //Double_t dfi0 = 2.*TMath::ASin(dvertex*c0*0.5);
2956 //Double_t dfi1 = 2.*TMath::ASin(TMath::Sqrt(yy0*yy0+(1-xx0)*(1-xx0))*dvertex*c0*0.5);
2957 Double_t dfi0 = 2.*AliTPCFastMath::FastAsin(dvertex*c0*0.5);
2958 Double_t dfi1 = 2.*AliTPCFastMath::FastAsin(TMath::Sqrt(yy0*yy0+(1-xx0)*(1-xx0))*dvertex*c0*0.5);
2961 Double_t z0 = kcl->GetZ();
2962 Double_t zzzz2 = z1-(z1-z3)*dfi1/dfi0;
2963 if (TMath::Abs(zzzz2-z0)>0.5) continue;
2966 Double_t dip = (z1-z0)*c0/dfi1;
2967 Double_t x0 = (0.5*cs13+y0*sn13)*dvertex*c0;
2978 x2= xx2*cs-y2*sn*dsec;
2979 y2=+xx2*sn*dsec+y2*cs;
2989 // do we have cluster at the middle ?
2991 GetProlongation(x1,xm,x,ym,zm);
2993 AliTPCclusterMI * cm=0;
2994 if (TMath::Abs(ym)-ymaxm<0){
2995 cm = krm.FindNearest2(ym,zm,1.0,0.6,dummy);
2996 if ((!cm) || (cm->IsUsed(10))) {
3001 // rotate y1 to system 0
3002 // get state vector in rotated system
3003 Double_t yr1 = (-0.5*sn13+y0*cs13)*dvertex*c0;
3004 Double_t xr2 = x0*cs+yr1*sn*dsec;
3005 Double_t xr[5]={kcl->GetY(),kcl->GetZ(), xr2, dip, c0};
3007 GetProlongation(xx2,xm,xr,ym,zm);
3008 if (TMath::Abs(ym)-ymaxm<0){
3009 cm = kr2m.FindNearest2(ym,zm,1.0,0.6,dummy);
3010 if ((!cm) || (cm->IsUsed(10))) {
3020 dym = ym - cm->GetY();
3021 dzm = zm - cm->GetZ();
3028 Double_t sy1=kr1[is]->GetSigmaY2()*2., sz1=kr1[is]->GetSigmaZ2()*2.;
3029 Double_t sy2=kcl->GetSigmaY2()*2., sz2=kcl->GetSigmaZ2()*2.;
3030 //Double_t sy3=400*3./12., sy=0.1, sz=0.1;
3031 Double_t sy3=25000*x[4]*x[4]+0.1, sy=0.1, sz=0.1;
3032 //Double_t sy3=25000*x[4]*x[4]*60+0.5, sy=0.1, sz=0.1;
3034 Double_t f40=(F1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
3035 Double_t f42=(F1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
3036 Double_t f43=(F1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
3037 Double_t f20=(F2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
3038 Double_t f22=(F2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
3039 Double_t f23=(F2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
3041 Double_t f30=(F3(x1,y1+sy,x2,y2,z1,z2)-x[3])/sy;
3042 Double_t f31=(F3(x1,y1,x2,y2,z1+sz,z2)-x[3])/sz;
3043 Double_t f32=(F3(x1,y1,x2,y2+sy,z1,z2)-x[3])/sy;
3044 Double_t f34=(F3(x1,y1,x2,y2,z1,z2+sz)-x[3])/sz;
3048 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3049 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3050 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3051 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3052 c[13]=f30*sy1*f40+f32*sy2*f42;
3053 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3056 UInt_t index=kr1.GetIndex(is);
3057 AliTPCseed track(x1, ns*alpha+shift, x, c, index);
3059 track.SetIsSeeding(kTRUE);
3062 track.SetSeedType(3);
3066 FollowProlongation(track, (i1+i2)/2,1);
3067 Int_t foundable,found,shared;
3068 track.GetClusterStatistic((i1+i2)/2,i1, found, foundable, shared, kTRUE);
3069 if ((found<0.55*foundable) || shared>0.5*found || (track.GetSigmaY2()+track.GetSigmaZ2())>0.5){
3075 FollowProlongation(track, i2,1);
3079 track.SetBConstrain(1);
3080 // track->fLastPoint = i1+fInnerSec->GetNRows(); // first cluster in track position
3081 track.SetLastPoint(i1); // first cluster in track position
3082 track.SetFirstPoint(track.GetLastPoint());
3084 if (track.GetNumberOfClusters()<(i1-i2)*0.5 ||
3085 track.GetNumberOfClusters() < track.GetNFoundable()*0.6 ||
3086 track.GetNShared()>0.4*track.GetNumberOfClusters() ) {
3090 // Z VERTEX CONDITION
3091 Double_t zv, bz=GetBz();
3092 if ( !track.GetZAt(0.,bz,zv) ) continue;
3093 if (TMath::Abs(zv-z3)>cuts[2]) {
3094 FollowProlongation(track, TMath::Max(i2-20,0));
3095 if ( track.GetZAt(0.,bz,zv) ) continue;
3096 if (TMath::Abs(zv-z3)>cuts[2]){
3097 // If track do not point to the primary vertex, but sufficientlu long
3098 //try to refit it without constrain
3101 FollowProlongation(track, TMath::Max(i2-40,0));
3102 if ( !track.GetZAt(0.,bz,zv) ) continue;
3103 if (TMath::Abs(zv-z3)>cuts[2] &&(track.GetNumberOfClusters() > track.GetNFoundable()*0.7)){
3104 // make seed without constrain
3105 AliTPCseed * track2 = MakeSeed(&track,0.2,0.5,1.);
3106 FollowProlongation(*track2, i2,1);
3107 track2->SetBConstrain(kFALSE);
3108 track2->SetSeedType(1);
3109 arr->AddLast(track2->Clone());
3114 //seed->~AliTPCseed();
3121 track.SetSeedType(0);
3122 arr->AddLast(track.Clone());
3123 //seed = new AliTPCseed;
3125 // don't consider other combinations
3126 if (track.GetNumberOfClusters() > track.GetNFoundable()*0.8)
3132 Info("MakeSeeds3","\nSeeding statistic:\t%d\t%d\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin,nout1,nout2);
3138 void AliTPCtrackerMI::MakeSeeds5(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, Float_t cuts[4],
3143 //-----------------------------------------------------------------
3144 // This function creates track seeds.
3145 //-----------------------------------------------------------------
3146 // cuts[0] - fP4 cut
3147 // cuts[1] - tan(phi) cut
3148 // cuts[2] - zvertex cut
3149 // cuts[3] - fP3 cut
3159 Double_t x[5], c[15];
3161 // make temporary seed
3162 AliTPCseed * seed = new AliTPCseed;
3163 Double_t alpha=fOuterSec->GetAlpha(), shift=fOuterSec->GetAlphaShift();
3164 // Double_t cs=cos(alpha), sn=sin(alpha);
3169 Double_t x1 = GetXrow(i1-1);
3170 const AliTPCtrackerRow& kr1=GetRow(sec,i1-1);
3171 Double_t y1max = GetMaxY(i1-1)-kr1.GetDeadZone()-1.5;
3173 Double_t x1p = GetXrow(i1);
3174 const AliTPCtrackerRow& kr1p=GetRow(sec,i1);
3176 Double_t x1m = GetXrow(i1-2);
3177 const AliTPCtrackerRow& kr1m=GetRow(sec,i1-2);
3180 //last 3 padrow for seeding
3181 AliTPCtrackerRow& kr3 = GetRow((sec+fkNOS)%fkNOS,i1-7);
3182 Double_t x3 = GetXrow(i1-7);
3183 // Double_t y3max= GetMaxY(i1-7)-kr3.fDeadZone-1.5;
3185 AliTPCtrackerRow& kr3p = GetRow((sec+fkNOS)%fkNOS,i1-6);
3186 Double_t x3p = GetXrow(i1-6);
3188 AliTPCtrackerRow& kr3m = GetRow((sec+fkNOS)%fkNOS,i1-8);
3189 Double_t x3m = GetXrow(i1-8);
3194 Int_t im = i1-4; //middle pad row index
3195 Double_t xm = GetXrow(im); // radius of middle pad-row
3196 const AliTPCtrackerRow& krm=GetRow(sec,im); //middle pad -row
3197 // Double_t ymmax = GetMaxY(im)-kr1.fDeadZone-1.5;
3200 Double_t deltax = x1-x3;
3201 Double_t dymax = deltax*cuts[1];
3202 Double_t dzmax = deltax*cuts[3];
3204 // loop over clusters
3205 for (Int_t is=0; is < kr1; is++) {
3207 if (kr1[is]->IsUsed(10)) continue;
3208 Double_t y1=kr1[is]->GetY(), z1=kr1[is]->GetZ();
3210 if (deltay>0 && TMath::Abs(y1max-TMath::Abs(y1))> deltay ) continue; // seed only at the edge
3212 Int_t index1 = TMath::Max(kr3.Find(z1-dzmax)-1,0);
3213 Int_t index2 = TMath::Min(kr3.Find(z1+dzmax)+1,kr3);
3219 for (Int_t js=index1; js < index2; js++) {
3220 const AliTPCclusterMI *kcl = kr3[js];
3221 if (kcl->IsUsed(10)) continue;
3223 // apply angular cuts
3224 if (TMath::Abs(y1-y3)>dymax) continue;
3227 if (TMath::Abs(z1-z3)>dzmax) continue;
3229 Double_t angley = (y1-y3)/(x1-x3);
3230 Double_t anglez = (z1-z3)/(x1-x3);
3232 Double_t erry = TMath::Abs(angley)*(x1-x1m)*0.5+0.5;
3233 Double_t errz = TMath::Abs(anglez)*(x1-x1m)*0.5+0.5;
3235 Double_t yyym = angley*(xm-x1)+y1;
3236 Double_t zzzm = anglez*(xm-x1)+z1;
3238 const AliTPCclusterMI *kcm = krm.FindNearest2(yyym,zzzm,erry,errz,index);
3240 if (kcm->IsUsed(10)) continue;
3242 erry = TMath::Abs(angley)*(x1-x1m)*0.4+0.5;
3243 errz = TMath::Abs(anglez)*(x1-x1m)*0.4+0.5;
3250 // look around first
3251 const AliTPCclusterMI *kc1m = kr1m.FindNearest2(angley*(x1m-x1)+y1,
3257 if (kc1m->IsUsed(10)) used++;
3259 const AliTPCclusterMI *kc1p = kr1p.FindNearest2(angley*(x1p-x1)+y1,
3265 if (kc1p->IsUsed(10)) used++;
3267 if (used>1) continue;
3268 if (found<1) continue;
3272 const AliTPCclusterMI *kc3m = kr3m.FindNearest2(angley*(x3m-x3)+y3,
3278 if (kc3m->IsUsed(10)) used++;
3282 const AliTPCclusterMI *kc3p = kr3p.FindNearest2(angley*(x3p-x3)+y3,
3288 if (kc3p->IsUsed(10)) used++;
3292 if (used>1) continue;
3293 if (found<3) continue;
3303 x[4]=F1(x1,y1,x2,y2,x3,y3);
3304 //if (TMath::Abs(x[4]) >= cuts[0]) continue;
3307 x[2]=F2(x1,y1,x2,y2,x3,y3);
3310 x[3]=F3n(x1,y1,x2,y2,z1,z2,x[4]);
3311 //if (TMath::Abs(x[3]) > cuts[3]) continue;
3315 Double_t sy1=0.1, sz1=0.1;
3316 Double_t sy2=0.1, sz2=0.1;
3317 Double_t sy3=0.1, sy=0.1, sz=0.1;
3319 Double_t f40=(F1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
3320 Double_t f42=(F1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
3321 Double_t f43=(F1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
3322 Double_t f20=(F2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
3323 Double_t f22=(F2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
3324 Double_t f23=(F2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
3326 Double_t f30=(F3(x1,y1+sy,x2,y2,z1,z2)-x[3])/sy;
3327 Double_t f31=(F3(x1,y1,x2,y2,z1+sz,z2)-x[3])/sz;
3328 Double_t f32=(F3(x1,y1,x2,y2+sy,z1,z2)-x[3])/sy;
3329 Double_t f34=(F3(x1,y1,x2,y2,z1,z2+sz)-x[3])/sz;
3333 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3334 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3335 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3336 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3337 c[13]=f30*sy1*f40+f32*sy2*f42;
3338 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3340 // if (!BuildSeed(kr1[is],kcl,0,x1,x2,x3,x,c)) continue;
3342 index=kr1.GetIndex(is);
3343 seed->~AliTPCseed();
3344 AliTPCseed *track=new(seed) AliTPCseed(x1, sec*alpha+shift, x, c, index);
3346 track->SetIsSeeding(kTRUE);
3349 FollowProlongation(*track, i1-7,1);
3350 if (track->GetNumberOfClusters() < track->GetNFoundable()*0.75 ||
3351 track->GetNShared()>0.6*track->GetNumberOfClusters() || ( track->GetSigmaY2()+ track->GetSigmaZ2())>0.6){
3353 seed->~AliTPCseed();
3359 FollowProlongation(*track, i2,1);
3360 track->SetBConstrain(0);
3361 track->SetLastPoint(i1+fInnerSec->GetNRows()); // first cluster in track position
3362 track->SetFirstPoint(track->GetLastPoint());
3364 if (track->GetNumberOfClusters()<(i1-i2)*0.5 ||
3365 track->GetNumberOfClusters()<track->GetNFoundable()*0.7 ||
3366 track->GetNShared()>2. || track->GetChi2()/track->GetNumberOfClusters()>6 || ( track->GetSigmaY2()+ track->GetSigmaZ2())>0.5 ) {
3368 seed->~AliTPCseed();
3373 FollowProlongation(*track, TMath::Max(i2-10,0),1);
3374 AliTPCseed * track2 = MakeSeed(track,0.2,0.5,0.9);
3375 FollowProlongation(*track2, i2,1);
3376 track2->SetBConstrain(kFALSE);
3377 track2->SetSeedType(4);
3378 arr->AddLast(track2);
3380 seed->~AliTPCseed();
3384 //arr->AddLast(track);
3385 //seed = new AliTPCseed;
3391 Info("MakeSeeds5","\nSeeding statiistic:\t%d\t%d\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin,nout1,nout2,nout3);
3397 //_____________________________________________________________________________
3398 void AliTPCtrackerMI::MakeSeeds2(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, Float_t */*cuts[4]*/,
3399 Float_t deltay, Bool_t /*bconstrain*/) {
3400 //-----------------------------------------------------------------
3401 // This function creates track seeds - without vertex constraint
3402 //-----------------------------------------------------------------
3403 // cuts[0] - fP4 cut - not applied
3404 // cuts[1] - tan(phi) cut
3405 // cuts[2] - zvertex cut - not applied
3406 // cuts[3] - fP3 cut
3416 Double_t alpha=fOuterSec->GetAlpha(), shift=fOuterSec->GetAlphaShift();
3417 // Double_t cs=cos(alpha), sn=sin(alpha);
3418 Int_t row0 = (i1+i2)/2;
3419 Int_t drow = (i1-i2)/2;
3420 const AliTPCtrackerRow& kr0=fSectors[sec][row0];
3421 AliTPCtrackerRow * kr=0;
3423 AliTPCpolyTrack polytrack;
3424 Int_t nclusters=fSectors[sec][row0];
3425 AliTPCseed * seed = new AliTPCseed;
3430 for (Int_t is=0; is < nclusters; is++) { //LOOP over clusters
3432 Int_t nfoundable =0;
3433 for (Int_t iter =1; iter<2; iter++){ //iterations
3434 const AliTPCtrackerRow& krm=fSectors[sec][row0-iter];
3435 const AliTPCtrackerRow& krp=fSectors[sec][row0+iter];
3436 const AliTPCclusterMI * cl= kr0[is];
3438 if (cl->IsUsed(10)) {
3444 Double_t x = kr0.GetX();
3445 // Initialization of the polytrack
3450 Double_t y0= cl->GetY();
3451 Double_t z0= cl->GetZ();
3455 Double_t ymax = fSectors->GetMaxY(row0)-kr0.GetDeadZone()-1.5;
3456 if (deltay>0 && TMath::Abs(ymax-TMath::Abs(y0))> deltay ) continue; // seed only at the edge
3458 erry = (0.5)*cl->GetSigmaY2()/TMath::Sqrt(cl->GetQ())*6;
3459 errz = (0.5)*cl->GetSigmaZ2()/TMath::Sqrt(cl->GetQ())*6;
3460 polytrack.AddPoint(x,y0,z0,erry, errz);
3463 if (cl->IsUsed(10)) sumused++;
3466 Float_t roady = (5*TMath::Sqrt(cl->GetSigmaY2()+0.2)+1.)*iter;
3467 Float_t roadz = (5*TMath::Sqrt(cl->GetSigmaZ2()+0.2)+1.)*iter;
3470 AliTPCclusterMI * cl1 = krm.FindNearest(y0,z0,roady,roadz);
3471 if (cl1 && TMath::Abs(ymax-TMath::Abs(y0))) {
3472 erry = (0.5)*cl1->GetSigmaY2()/TMath::Sqrt(cl1->GetQ())*3;
3473 errz = (0.5)*cl1->GetSigmaZ2()/TMath::Sqrt(cl1->GetQ())*3;
3474 if (cl1->IsUsed(10)) sumused++;
3475 polytrack.AddPoint(x,cl1->GetY(),cl1->GetZ(),erry,errz);
3479 AliTPCclusterMI * cl2 = krp.FindNearest(y0,z0,roady,roadz);
3481 erry = (0.5)*cl2->GetSigmaY2()/TMath::Sqrt(cl2->GetQ())*3;
3482 errz = (0.5)*cl2->GetSigmaZ2()/TMath::Sqrt(cl2->GetQ())*3;
3483 if (cl2->IsUsed(10)) sumused++;
3484 polytrack.AddPoint(x,cl2->GetY(),cl2->GetZ(),erry,errz);
3487 if (sumused>0) continue;
3489 polytrack.UpdateParameters();
3495 nfoundable = polytrack.GetN();
3496 nfound = nfoundable;
3498 for (Int_t ddrow = iter+1; ddrow<drow;ddrow++){
3499 Float_t maxdist = 0.8*(1.+3./(ddrow));
3500 for (Int_t delta = -1;delta<=1;delta+=2){
3501 Int_t row = row0+ddrow*delta;
3502 kr = &(fSectors[sec][row]);
3503 Double_t xn = kr->GetX();
3504 Double_t ymax1 = fSectors->GetMaxY(row)-kr->GetDeadZone()-1.5;
3505 polytrack.GetFitPoint(xn,yn,zn);
3506 if (TMath::Abs(yn)>ymax1) continue;
3508 AliTPCclusterMI * cln = kr->FindNearest(yn,zn,roady,roadz);
3510 Float_t dist = TMath::Sqrt( (yn-cln->GetY())*(yn-cln->GetY())+(zn-cln->GetZ())*(zn-cln->GetZ()));
3513 erry = (dist+0.3)*cln->GetSigmaY2()/TMath::Sqrt(cln->GetQ())*(1.+1./(ddrow));
3514 errz = (dist+0.3)*cln->GetSigmaZ2()/TMath::Sqrt(cln->GetQ())*(1.+1./(ddrow));
3515 if (cln->IsUsed(10)) {
3516 // printf("used\n");
3524 polytrack.AddPoint(xn,cln->GetY(),cln->GetZ(),erry, errz);
3529 if ( (sumused>3) || (sumused>0.5*nfound) || (nfound<0.6*nfoundable)) break;
3530 polytrack.UpdateParameters();
3533 if ( (sumused>3) || (sumused>0.5*nfound)) {
3534 //printf("sumused %d\n",sumused);
3539 polytrack.GetFitDerivation(kr0.GetX(),dy,dz);
3540 AliTPCpolyTrack track2;
3542 polytrack.Refit(track2,0.5+TMath::Abs(dy)*0.3,0.4+TMath::Abs(dz)*0.3);
3543 if (track2.GetN()<0.5*nfoundable) continue;
3546 if ((nfound>0.6*nfoundable) &&( nfoundable>0.4*(i1-i2))) {
3548 // test seed with and without constrain
3549 for (Int_t constrain=0; constrain<=0;constrain++){
3550 // add polytrack candidate
3552 Double_t x[5], c[15];
3553 Double_t x1,x2,x3,y1,y2,y3,z1,z2,z3;
3554 track2.GetBoundaries(x3,x1);
3556 track2.GetFitPoint(x1,y1,z1);
3557 track2.GetFitPoint(x2,y2,z2);
3558 track2.GetFitPoint(x3,y3,z3);
3560 //is track pointing to the vertex ?
3563 polytrack.GetFitPoint(x0,y0,z0);
3576 x[4]=F1(x1,y1,x2,y2,x3,y3);
3578 // if (TMath::Abs(x[4]) >= cuts[0]) continue; //
3579 x[2]=F2(x1,y1,x2,y2,x3,y3);
3581 //if (TMath::Abs(x[4]*x1-x[2]) >= cuts[1]) continue;
3582 //x[3]=F3(x1,y1,x2,y2,z1,z2);
3583 x[3]=F3n(x1,y1,x3,y3,z1,z3,x[4]);
3584 //if (TMath::Abs(x[3]) > cuts[3]) continue;
3587 Double_t sy =0.1, sz =0.1;
3588 Double_t sy1=0.02, sz1=0.02;
3589 Double_t sy2=0.02, sz2=0.02;
3593 sy3=25000*x[4]*x[4]+0.1, sy=0.1, sz=0.1;
3596 Double_t f40=(F1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
3597 Double_t f42=(F1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
3598 Double_t f43=(F1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
3599 Double_t f20=(F2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
3600 Double_t f22=(F2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
3601 Double_t f23=(F2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
3603 Double_t f30=(F3(x1,y1+sy,x3,y3,z1,z3)-x[3])/sy;
3604 Double_t f31=(F3(x1,y1,x3,y3,z1+sz,z3)-x[3])/sz;
3605 Double_t f32=(F3(x1,y1,x3,y3+sy,z1,z3)-x[3])/sy;
3606 Double_t f34=(F3(x1,y1,x3,y3,z1,z3+sz)-x[3])/sz;
3611 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3612 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3613 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3614 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3615 c[13]=f30*sy1*f40+f32*sy2*f42;
3616 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3618 //Int_t row1 = fSectors->GetRowNumber(x1);
3619 Int_t row1 = GetRowNumber(x1);
3623 seed->~AliTPCseed();
3624 AliTPCseed *track=new(seed) AliTPCseed(x1,sec*alpha+shift,x,c,index);
3625 track->SetIsSeeding(kTRUE);
3626 Int_t rc=FollowProlongation(*track, i2);
3627 if (constrain) track->SetBConstrain(1);
3629 track->SetBConstrain(0);
3630 track->SetLastPoint(row1+fInnerSec->GetNRows()); // first cluster in track position
3631 track->SetFirstPoint(track->GetLastPoint());
3633 if (rc==0 || track->GetNumberOfClusters()<(i1-i2)*0.5 ||
3634 track->GetNumberOfClusters() < track->GetNFoundable()*0.6 ||
3635 track->GetNShared()>0.4*track->GetNumberOfClusters()) {
3638 seed->~AliTPCseed();
3641 arr->AddLast(track);
3642 seed = new AliTPCseed;
3646 } // if accepted seed
3649 Info("MakeSeeds2","\nSeeding statiistic:\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin3);
3655 AliTPCseed *AliTPCtrackerMI::MakeSeed(AliTPCseed *track, Float_t r0, Float_t r1, Float_t r2)
3659 //reseed using track points
3660 Int_t p0 = int(r0*track->GetNumberOfClusters()); // point 0
3661 Int_t p1 = int(r1*track->GetNumberOfClusters());
3662 Int_t p2 = int(r2*track->GetNumberOfClusters()); // last point
3664 Double_t x0[3],x1[3],x2[3];
3665 for (Int_t i=0;i<3;i++){
3671 // find track position at given ratio of the length
3672 Int_t sec0=0, sec1=0, sec2=0;
3675 for (Int_t i=0;i<160;i++){
3676 if (track->GetClusterPointer(i)){
3678 AliTPCTrackerPoint *trpoint =track->GetTrackPoint(i);
3679 if ( (index<p0) || x0[0]<0 ){
3680 if (trpoint->GetX()>1){
3681 clindex = track->GetClusterIndex2(i);
3683 x0[0] = trpoint->GetX();
3684 x0[1] = trpoint->GetY();
3685 x0[2] = trpoint->GetZ();
3686 sec0 = ((clindex&0xff000000)>>24)%18;
3691 if ( (index<p1) &&(trpoint->GetX()>1)){
3692 clindex = track->GetClusterIndex2(i);
3694 x1[0] = trpoint->GetX();
3695 x1[1] = trpoint->GetY();
3696 x1[2] = trpoint->GetZ();
3697 sec1 = ((clindex&0xff000000)>>24)%18;
3700 if ( (index<p2) &&(trpoint->GetX()>1)){
3701 clindex = track->GetClusterIndex2(i);
3703 x2[0] = trpoint->GetX();
3704 x2[1] = trpoint->GetY();
3705 x2[2] = trpoint->GetZ();
3706 sec2 = ((clindex&0xff000000)>>24)%18;
3713 Double_t alpha, cs,sn, xx2,yy2;
3715 alpha = (sec1-sec2)*fSectors->GetAlpha();
3716 cs = TMath::Cos(alpha);
3717 sn = TMath::Sin(alpha);
3718 xx2= x1[0]*cs-x1[1]*sn;
3719 yy2= x1[0]*sn+x1[1]*cs;
3723 alpha = (sec0-sec2)*fSectors->GetAlpha();
3724 cs = TMath::Cos(alpha);
3725 sn = TMath::Sin(alpha);
3726 xx2= x0[0]*cs-x0[1]*sn;
3727 yy2= x0[0]*sn+x0[1]*cs;
3733 Double_t x[5],c[15];
3737 x[4]=F1(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]);
3738 // if (x[4]>1) return 0;
3739 x[2]=F2(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]);
3740 x[3]=F3n(x2[0],x2[1],x0[0],x0[1],x2[2],x0[2],x[4]);
3741 //if (TMath::Abs(x[3]) > 2.2) return 0;
3742 //if (TMath::Abs(x[2]) > 1.99) return 0;
3744 Double_t sy =0.1, sz =0.1;
3746 Double_t sy1=0.02+track->GetSigmaY2(), sz1=0.02+track->GetSigmaZ2();
3747 Double_t sy2=0.01+track->GetSigmaY2(), sz2=0.01+track->GetSigmaZ2();
3748 Double_t sy3=0.01+track->GetSigmaY2();
3750 Double_t f40=(F1(x2[0],x2[1]+sy,x1[0],x1[1],x0[0],x0[1])-x[4])/sy;
3751 Double_t f42=(F1(x2[0],x2[1],x1[0],x1[1]+sy,x0[0],x0[1])-x[4])/sy;
3752 Double_t f43=(F1(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]+sy)-x[4])/sy;
3753 Double_t f20=(F2(x2[0],x2[1]+sy,x1[0],x1[1],x0[0],x0[1])-x[2])/sy;
3754 Double_t f22=(F2(x2[0],x2[1],x1[0],x1[1]+sy,x0[0],x0[1])-x[2])/sy;
3755 Double_t f23=(F2(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]+sy)-x[2])/sy;
3757 Double_t f30=(F3(x2[0],x2[1]+sy,x0[0],x0[1],x2[2],x0[2])-x[3])/sy;
3758 Double_t f31=(F3(x2[0],x2[1],x0[0],x0[1],x2[2]+sz,x0[2])-x[3])/sz;
3759 Double_t f32=(F3(x2[0],x2[1],x0[0],x0[1]+sy,x2[2],x0[2])-x[3])/sy;
3760 Double_t f34=(F3(x2[0],x2[1],x0[0],x0[1],x2[2],x0[2]+sz)-x[3])/sz;
3765 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3766 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3767 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3768 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3769 c[13]=f30*sy1*f40+f32*sy2*f42;
3770 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3772 // Int_t row1 = fSectors->GetRowNumber(x2[0]);
3773 AliTPCseed *seed=new AliTPCseed(x2[0], sec2*fSectors->GetAlpha()+fSectors->GetAlphaShift(), x, c, 0);
3774 // Double_t y0,z0,y1,z1, y2,z2;
3775 //seed->GetProlongation(x0[0],y0,z0);
3776 // seed->GetProlongation(x1[0],y1,z1);
3777 //seed->GetProlongation(x2[0],y2,z2);
3779 seed->SetLastPoint(pp2);
3780 seed->SetFirstPoint(pp2);
3787 AliTPCseed *AliTPCtrackerMI::ReSeed(AliTPCseed *track, Float_t r0, Float_t r1, Float_t r2)
3791 //reseed using founded clusters
3793 // Find the number of clusters
3794 Int_t nclusters = 0;
3795 for (Int_t irow=0;irow<160;irow++){
3796 if (track->GetClusterIndex(irow)>0) nclusters++;
3800 ipos[0] = TMath::Max(int(r0*nclusters),0); // point 0 cluster
3801 ipos[1] = TMath::Min(int(r1*nclusters),nclusters-1); //
3802 ipos[2] = TMath::Min(int(r2*nclusters),nclusters-1); // last point
3806 Int_t row[3],sec[3]={0,0,0};
3808 // find track row position at given ratio of the length
3810 for (Int_t irow=0;irow<160;irow++){
3811 if (track->GetClusterIndex2(irow)<0) continue;
3813 for (Int_t ipoint=0;ipoint<3;ipoint++){
3814 if (index<=ipos[ipoint]) row[ipoint] = irow;
3818 //Get cluster and sector position
3819 for (Int_t ipoint=0;ipoint<3;ipoint++){
3820 Int_t clindex = track->GetClusterIndex2(row[ipoint]);
3821 AliTPCclusterMI * cl = GetClusterMI(clindex);
3824 // AliTPCclusterMI * cl = GetClusterMI(clindex);
3827 sec[ipoint] = ((clindex&0xff000000)>>24)%18;
3828 xyz[ipoint][0] = GetXrow(row[ipoint]);
3829 xyz[ipoint][1] = cl->GetY();
3830 xyz[ipoint][2] = cl->GetZ();
3834 // Calculate seed state vector and covariance matrix
3836 Double_t alpha, cs,sn, xx2,yy2;
3838 alpha = (sec[1]-sec[2])*fSectors->GetAlpha();
3839 cs = TMath::Cos(alpha);
3840 sn = TMath::Sin(alpha);
3841 xx2= xyz[1][0]*cs-xyz[1][1]*sn;
3842 yy2= xyz[1][0]*sn+xyz[1][1]*cs;
3846 alpha = (sec[0]-sec[2])*fSectors->GetAlpha();
3847 cs = TMath::Cos(alpha);
3848 sn = TMath::Sin(alpha);
3849 xx2= xyz[0][0]*cs-xyz[0][1]*sn;
3850 yy2= xyz[0][0]*sn+xyz[0][1]*cs;
3856 Double_t x[5],c[15];
3860 x[4]=F1(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
3861 x[2]=F2(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
3862 x[3]=F3n(xyz[2][0],xyz[2][1],xyz[0][0],xyz[0][1],xyz[2][2],xyz[0][2],x[4]);
3864 Double_t sy =0.1, sz =0.1;
3866 Double_t sy1=0.2, sz1=0.2;
3867 Double_t sy2=0.2, sz2=0.2;
3870 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;
3871 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;
3872 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;
3873 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;
3874 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;
3875 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;
3877 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;
3878 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;
3879 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;
3880 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;
3885 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3886 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3887 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3888 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3889 c[13]=f30*sy1*f40+f32*sy2*f42;
3890 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3892 // Int_t row1 = fSectors->GetRowNumber(xyz[2][0]);
3893 AliTPCseed *seed=new AliTPCseed(xyz[2][0], sec[2]*fSectors->GetAlpha()+fSectors->GetAlphaShift(), x, c, 0);
3894 seed->SetLastPoint(row[2]);
3895 seed->SetFirstPoint(row[2]);
3900 AliTPCseed *AliTPCtrackerMI::ReSeed(AliTPCseed *track,Int_t r0, Bool_t forward)
3904 //reseed using founded clusters
3907 Int_t row[3]={0,0,0};
3908 Int_t sec[3]={0,0,0};
3910 // forward direction
3912 for (Int_t irow=r0;irow<160;irow++){
3913 if (track->GetClusterIndex(irow)>0){
3918 for (Int_t irow=160;irow>r0;irow--){
3919 if (track->GetClusterIndex(irow)>0){
3924 for (Int_t irow=row[2]-15;irow>row[0];irow--){
3925 if (track->GetClusterIndex(irow)>0){
3933 for (Int_t irow=0;irow<r0;irow++){
3934 if (track->GetClusterIndex(irow)>0){
3939 for (Int_t irow=r0;irow>0;irow--){
3940 if (track->GetClusterIndex(irow)>0){
3945 for (Int_t irow=row[2]-15;irow>row[0];irow--){
3946 if (track->GetClusterIndex(irow)>0){
3953 if ((row[2]-row[0])<20) return 0;
3954 if (row[1]==0) return 0;
3957 //Get cluster and sector position
3958 for (Int_t ipoint=0;ipoint<3;ipoint++){
3959 Int_t clindex = track->GetClusterIndex2(row[ipoint]);
3960 AliTPCclusterMI * cl = GetClusterMI(clindex);
3963 // AliTPCclusterMI * cl = GetClusterMI(clindex);
3966 sec[ipoint] = ((clindex&0xff000000)>>24)%18;
3967 xyz[ipoint][0] = GetXrow(row[ipoint]);
3968 AliTPCTrackerPoint * point = track->GetTrackPoint(row[ipoint]);
3969 if (point&&ipoint<2){
3971 xyz[ipoint][1] = point->GetY();
3972 xyz[ipoint][2] = point->GetZ();
3975 xyz[ipoint][1] = cl->GetY();
3976 xyz[ipoint][2] = cl->GetZ();
3983 // Calculate seed state vector and covariance matrix
3985 Double_t alpha, cs,sn, xx2,yy2;
3987 alpha = (sec[1]-sec[2])*fSectors->GetAlpha();
3988 cs = TMath::Cos(alpha);
3989 sn = TMath::Sin(alpha);
3990 xx2= xyz[1][0]*cs-xyz[1][1]*sn;
3991 yy2= xyz[1][0]*sn+xyz[1][1]*cs;
3995 alpha = (sec[0]-sec[2])*fSectors->GetAlpha();
3996 cs = TMath::Cos(alpha);
3997 sn = TMath::Sin(alpha);
3998 xx2= xyz[0][0]*cs-xyz[0][1]*sn;
3999 yy2= xyz[0][0]*sn+xyz[0][1]*cs;
4005 Double_t x[5],c[15];
4009 x[4]=F1(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
4010 x[2]=F2(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
4011 x[3]=F3n(xyz[2][0],xyz[2][1],xyz[0][0],xyz[0][1],xyz[2][2],xyz[0][2],x[4]);
4013 Double_t sy =0.1, sz =0.1;
4015 Double_t sy1=0.2, sz1=0.2;
4016 Double_t sy2=0.2, sz2=0.2;
4019 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;
4020 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;
4021 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;
4022 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;
4023 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;
4024 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;
4026 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;
4027 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;
4028 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;
4029 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;
4034 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
4035 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
4036 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
4037 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
4038 c[13]=f30*sy1*f40+f32*sy2*f42;
4039 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
4041 // Int_t row1 = fSectors->GetRowNumber(xyz[2][0]);
4042 AliTPCseed *seed=new AliTPCseed(xyz[2][0], sec[2]*fSectors->GetAlpha()+fSectors->GetAlphaShift(), x, c, 0);
4043 seed->SetLastPoint(row[2]);
4044 seed->SetFirstPoint(row[2]);
4045 for (Int_t i=row[0];i<row[2];i++){
4046 seed->SetClusterIndex(i, track->GetClusterIndex(i));
4054 void AliTPCtrackerMI::FindMultiMC(TObjArray * array, AliESDEvent */*esd*/, Int_t iter)
4057 // find multi tracks - THIS FUNCTION IS ONLY FOR DEBUG PURPOSES
4059 // Use AliTPCReconstructor::StreamLevel()>2 if you want to tune parameters - cuts
4061 // Two reasons to have multiple find tracks
4062 // 1. Curling tracks can be find more than once
4063 // 2. Splitted tracks
4064 // a.) Multiple seeding to increase tracking efficiency - (~ 100% reached)
4065 // b.) Edge effect on the sector boundaries
4068 // Algorithm done in 2 phases - because of CPU consumption
4069 // it is n^2 algorithm - for lead-lead 20000x20000 combination are investigated
4071 // Algorihm for curling tracks sign:
4072 // 1 phase -makes a very rough fast cuts to minimize combinatorics
4073 // a.) opposite sign
4074 // b.) one of the tracks - not pointing to the primary vertex -
4075 // c.) delta tan(theta)
4077 // 2 phase - calculates DCA between tracks - time consument
4082 // General cuts - for splitted tracks and for curling tracks
4084 const Float_t kMaxdPhi = 0.2; // maximal distance in phi
4086 // Curling tracks cuts
4091 Int_t nentries = array->GetEntriesFast();
4092 AliHelix *helixes = new AliHelix[nentries];
4093 Float_t *xm = new Float_t[nentries];
4094 Float_t *dz0 = new Float_t[nentries];
4095 Float_t *dz1 = new Float_t[nentries];
4101 // Find track COG in x direction - point with best defined parameters
4103 for (Int_t i=0;i<nentries;i++){
4104 AliTPCseed* track = (AliTPCseed*)array->At(i);
4105 if (!track) continue;
4106 track->SetCircular(0);
4107 new (&helixes[i]) AliHelix(*track);
4111 track->GetDZ(GetX(),GetY(),GetZ(),GetBz(),dz);
4114 for (Int_t icl=0; icl<160; icl++){
4115 AliTPCclusterMI * cl = track->GetClusterPointer(icl);
4121 if (ncl>0) xm[i]/=Float_t(ncl);
4124 for (Int_t i0=0;i0<nentries;i0++){
4125 AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
4126 if (!track0) continue;
4127 Float_t xc0 = helixes[i0].GetHelix(6);
4128 Float_t yc0 = helixes[i0].GetHelix(7);
4129 Float_t r0 = helixes[i0].GetHelix(8);
4130 Float_t rc0 = TMath::Sqrt(xc0*xc0+yc0*yc0);
4131 Float_t fi0 = TMath::ATan2(yc0,xc0);
4133 for (Int_t i1=i0+1;i1<nentries;i1++){
4134 AliTPCseed * track1 = (AliTPCseed*)array->At(i1);
4135 if (!track1) continue;
4136 Int_t lab0=track0->GetLabel();
4137 Int_t lab1=track1->GetLabel();
4138 if (TMath::Abs(lab0)!=TMath::Abs(lab1)) continue;
4140 Float_t xc1 = helixes[i1].GetHelix(6);
4141 Float_t yc1 = helixes[i1].GetHelix(7);
4142 Float_t r1 = helixes[i1].GetHelix(8);
4143 Float_t rc1 = TMath::Sqrt(xc1*xc1+yc1*yc1);
4144 Float_t fi1 = TMath::ATan2(yc1,xc1);
4146 Float_t dfi = fi0-fi1;
4149 if (dfi>1.5*TMath::Pi()) dfi-=TMath::Pi(); // take care about edge effect
4150 if (dfi<-1.5*TMath::Pi()) dfi+=TMath::Pi(); //
4151 if (TMath::Abs(dfi)>kMaxdPhi&&helixes[i0].GetHelix(4)*helixes[i1].GetHelix(4)<0){
4153 // if short tracks with undefined sign
4154 fi1 = -TMath::ATan2(yc1,-xc1);
4157 Float_t dtheta = TMath::Abs(track0->GetTgl()-track1->GetTgl())<TMath::Abs(track0->GetTgl()+track1->GetTgl())? track0->GetTgl()-track1->GetTgl():track0->GetTgl()+track1->GetTgl();
4160 // debug stream to tune "fast cuts"
4162 Double_t dist[3]; // distance at X
4163 Double_t mdist[3]={0,0,0}; // mean distance X+-40cm
4164 track0->GetDistance(track1,0.5*(xm[i0]+xm[i1])-40.,dist,AliTracker::GetBz());
4165 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[i]);
4166 track0->GetDistance(track1,0.5*(xm[i0]+xm[i1])+40.,dist,AliTracker::GetBz());
4167 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[i]);
4168 track0->GetDistance(track1,0.5*(xm[i0]+xm[i1]),dist,AliTracker::GetBz());
4169 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[i]);
4170 for (Int_t i=0;i<3;i++) mdist[i]*=0.33333;
4174 for (Int_t icl=0; icl<160; icl++){
4175 AliTPCclusterMI * cl0 = track0->GetClusterPointer(icl);
4176 AliTPCclusterMI * cl1 = track1->GetClusterPointer(icl);
4179 if (cl0==cl1) sums++;
4183 if (AliTPCReconstructor::StreamLevel()>0) {
4184 TTreeSRedirector &cstream = *fDebugStreamer;
4189 "Tr0.="<<track0<< // seed0
4190 "Tr1.="<<track1<< // seed1
4191 "h0.="<<&helixes[i0]<<
4192 "h1.="<<&helixes[i1]<<
4194 "sum="<<sum<< //the sum of rows with cl in both
4195 "sums="<<sums<< //the sum of shared clusters
4196 "xm0="<<xm[i0]<< // the center of track
4197 "xm1="<<xm[i1]<< // the x center of track
4198 // General cut variables
4199 "dfi="<<dfi<< // distance in fi angle
4200 "dtheta="<<dtheta<< // distance int theta angle
4206 "dist0="<<dist[0]<< //distance x
4207 "dist1="<<dist[1]<< //distance y
4208 "dist2="<<dist[2]<< //distance z
4209 "mdist0="<<mdist[0]<< //distance x
4210 "mdist1="<<mdist[1]<< //distance y
4211 "mdist2="<<mdist[2]<< //distance z
4225 if (AliTPCReconstructor::StreamLevel()>1) {
4226 AliInfo("Time for curling tracks removal DEBUGGING MC");
4232 void AliTPCtrackerMI::FindSplitted(TObjArray * array, AliESDEvent */*esd*/, Int_t iter)
4236 // Two reasons to have multiple find tracks
4237 // 1. Curling tracks can be find more than once
4238 // 2. Splitted tracks
4239 // a.) Multiple seeding to increase tracking efficiency - (~ 100% reached)
4240 // b.) Edge effect on the sector boundaries
4242 // This function tries to find tracks closed in the parametric space
4244 // cut logic if distance is bigger than cut continue - Do Nothing
4245 const Float_t kMaxdTheta = 0.05; // maximal distance in theta
4246 const Float_t kMaxdPhi = 0.05; // maximal deistance in phi
4247 const Float_t kdelta = 40.; //delta r to calculate track distance
4249 // const Float_t kMaxDist0 = 1.; // maximal distance 0
4250 //const Float_t kMaxDist1 = 0.3; // maximal distance 1 - cut if track in separate rows
4253 TCut csec("csec","abs(Tr0.fRelativeSector-Tr1.fRelativeSector)<2");
4254 TCut cdtheta("cdtheta","abs(dtheta)<0.05");
4259 Int_t nentries = array->GetEntriesFast();
4260 AliHelix *helixes = new AliHelix[nentries];
4261 Float_t *xm = new Float_t[nentries];
4267 //Sort tracks according quality
4269 Int_t nseed = array->GetEntriesFast();
4270 Float_t * quality = new Float_t[nseed];
4271 Int_t * indexes = new Int_t[nseed];
4272 for (Int_t i=0; i<nseed; i++) {
4273 AliTPCseed *pt=(AliTPCseed*)array->UncheckedAt(i);
4278 pt->UpdatePoints(); //select first last max dens points
4279 Float_t * points = pt->GetPoints();
4280 if (points[3]<0.8) quality[i] =-1;
4281 quality[i] = (points[2]-points[0])+pt->GetNumberOfClusters();
4282 //prefer high momenta tracks if overlaps
4283 quality[i] *= TMath::Sqrt(TMath::Abs(pt->Pt())+0.5);
4285 TMath::Sort(nseed,quality,indexes);
4289 // Find track COG in x direction - point with best defined parameters
4291 for (Int_t i=0;i<nentries;i++){
4292 AliTPCseed* track = (AliTPCseed*)array->At(i);
4293 if (!track) continue;
4294 track->SetCircular(0);
4295 new (&helixes[i]) AliHelix(*track);
4298 for (Int_t icl=0; icl<160; icl++){
4299 AliTPCclusterMI * cl = track->GetClusterPointer(icl);
4305 if (ncl>0) xm[i]/=Float_t(ncl);
4308 for (Int_t is0=0;is0<nentries;is0++){
4309 Int_t i0 = indexes[is0];
4310 AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
4311 if (!track0) continue;
4312 if (track0->GetKinkIndexes()[0]!=0) continue;
4313 Float_t xc0 = helixes[i0].GetHelix(6);
4314 Float_t yc0 = helixes[i0].GetHelix(7);
4315 Float_t fi0 = TMath::ATan2(yc0,xc0);
4317 for (Int_t is1=is0+1;is1<nentries;is1++){
4318 Int_t i1 = indexes[is1];
4319 AliTPCseed * track1 = (AliTPCseed*)array->At(i1);
4320 if (!track1) continue;
4322 if (TMath::Abs(track0->GetRelativeSector()-track1->GetRelativeSector())>1) continue;
4323 if (track1->GetKinkIndexes()[0]>0 &&track0->GetKinkIndexes()[0]<0) continue;
4324 if (track1->GetKinkIndexes()[0]!=0) continue;
4326 Float_t dtheta = TMath::Abs(track0->GetTgl()-track1->GetTgl())<TMath::Abs(track0->GetTgl()+track1->GetTgl())? track0->GetTgl()-track1->GetTgl():track0->GetTgl()+track1->GetTgl();
4327 if (TMath::Abs(dtheta)>kMaxdTheta) continue;
4329 Float_t xc1 = helixes[i1].GetHelix(6);
4330 Float_t yc1 = helixes[i1].GetHelix(7);
4331 Float_t fi1 = TMath::ATan2(yc1,xc1);
4333 Float_t dfi = fi0-fi1;
4334 if (dfi>1.5*TMath::Pi()) dfi-=TMath::Pi(); // take care about edge effect
4335 if (dfi<-1.5*TMath::Pi()) dfi+=TMath::Pi(); //
4336 if (TMath::Abs(dfi)>kMaxdPhi&&helixes[i0].GetHelix(4)*helixes[i1].GetHelix(4)<0){
4338 // if short tracks with undefined sign
4339 fi1 = -TMath::ATan2(yc1,-xc1);
4342 if (TMath::Abs(dfi)>kMaxdPhi) continue;
4349 for (Int_t icl=0; icl<160; icl++){
4350 Int_t index0=track0->GetClusterIndex2(icl);
4351 Int_t index1=track1->GetClusterIndex2(icl);
4352 Bool_t used0 = (index0>0 && !(index0&0x8000));
4353 Bool_t used1 = (index1>0 && !(index1&0x8000));
4355 if (used0) sum0++; // used cluster0
4356 if (used1) sum1++; // used clusters1
4357 if (used0&&used1) sum++;
4358 if (index0==index1 && used0 && used1) sums++;
4362 if (sums<10) continue;
4363 if (sum<40) continue;
4364 if (sums/Float_t(TMath::Min(sum0,sum1))<0.5) continue;
4366 Double_t dist[5][4]; // distance at X
4367 Double_t mdist[4]={0,0,0,0}; // mean distance on range +-delta
4371 track0->GetDistance(track1,xm[i0],dist[0],AliTracker::GetBz());
4372 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[0][i]);
4373 track0->GetDistance(track1,xm[i1],dist[1],AliTracker::GetBz());
4374 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[1][i]);
4376 track0->GetDistance(track1,TMath::Min(xm[i1],xm[i0])-kdelta,dist[2],AliTracker::GetBz());
4377 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[2][i]);
4378 track0->GetDistance(track1,TMath::Max(xm[i1],xm[i0])+kdelta,dist[3],AliTracker::GetBz());
4379 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[3][i]);
4381 track0->GetDistance(track1,(xm[i1]+xm[i0])*0.5,dist[4],AliTracker::GetBz());
4382 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[4][i]);
4383 for (Int_t i=0;i<3;i++) mdist[i]*=0.2;
4386 Int_t lab0=track0->GetLabel();
4387 Int_t lab1=track1->GetLabel();
4388 if( AliTPCReconstructor::StreamLevel()>5){
4389 TTreeSRedirector &cstream = *fDebugStreamer;
4390 cstream<<"Splitted"<<
4394 "Tr0.="<<track0<< // seed0
4395 "Tr1.="<<track1<< // seed1
4396 "h0.="<<&helixes[i0]<<
4397 "h1.="<<&helixes[i1]<<
4399 "sum="<<sum<< //the sum of rows with cl in both
4400 "sum0="<<sum0<< //the sum of rows with cl in 0 track
4401 "sum1="<<sum1<< //the sum of rows with cl in 1 track
4402 "sums="<<sums<< //the sum of shared clusters
4403 "xm0="<<xm[i0]<< // the center of track
4404 "xm1="<<xm[i1]<< // the x center of track
4405 // General cut variables
4406 "dfi="<<dfi<< // distance in fi angle
4407 "dtheta="<<dtheta<< // distance int theta angle
4410 "dist0="<<dist[4][0]<< //distance x
4411 "dist1="<<dist[4][1]<< //distance y
4412 "dist2="<<dist[4][2]<< //distance z
4413 "mdist0="<<mdist[0]<< //distance x
4414 "mdist1="<<mdist[1]<< //distance y
4415 "mdist2="<<mdist[2]<< //distance z
4418 delete array->RemoveAt(i1);
4425 AliInfo("Time for splitted tracks removal");
4431 void AliTPCtrackerMI::FindCurling(TObjArray * array, AliESDEvent */*esd*/, Int_t iter)
4434 // find Curling tracks
4435 // Use AliTPCReconstructor::StreamLevel()>1 if you want to tune parameters - cuts
4438 // Algorithm done in 2 phases - because of CPU consumption
4439 // it is n^2 algorithm - for lead-lead 20000x20000 combination are investigated
4440 // see detal in MC part what can be used to cut
4444 const Float_t kMaxC = 400; // maximal curvature to of the track
4445 const Float_t kMaxdTheta = 0.15; // maximal distance in theta
4446 const Float_t kMaxdPhi = 0.15; // maximal distance in phi
4447 const Float_t kPtRatio = 0.3; // ratio between pt
4448 const Float_t kMinDCAR = 2.; // distance to the primary vertex in r - see cpipe cut
4451 // Curling tracks cuts
4454 const Float_t kMaxDeltaRMax = 40; // distance in outer radius
4455 const Float_t kMaxDeltaRMin = 5.; // distance in lower radius - see cpipe cut
4456 const Float_t kMinAngle = 2.9; // angle between tracks
4457 const Float_t kMaxDist = 5; // biggest distance
4459 // The cuts can be tuned using the "MC information stored in Multi tree ==> see FindMultiMC
4462 TCut csign("csign","Tr0.fP[4]*Tr1.fP[4]<0"); //opposite sign
4463 TCut cmax("cmax","abs(Tr0.GetC())>1/400");
4464 TCut cda("cda","sqrt(dtheta^2+dfi^2)<0.15");
4465 TCut ccratio("ccratio","abs((Tr0.fP[4]+Tr1.fP[4])/(abs(Tr0.fP[4])+abs(Tr1.fP[4])))<0.3");
4466 TCut cpipe("cpipe", "min(abs(r0-rc0),abs(r1-rc1))>5");
4468 TCut cdrmax("cdrmax","abs(abs(rc0+r0)-abs(rc1+r1))<40")
4469 TCut cdrmin("cdrmin","abs(abs(rc0+r0)-abs(rc1+r1))<10")
4471 Multi->Draw("dfi","iter==0"+csign+cmax+cda+ccratio); ~94% of curling tracks fulfill
4472 Multi->Draw("min(abs(r0-rc0),abs(r1-rc1))","iter==0&&abs(lab1)==abs(lab0)"+csign+cmax+cda+ccratio+cpipe+cdrmin+cdrmax); //80%
4474 Curling2->Draw("dfi","iter==0&&abs(lab0)==abs(lab1)"+csign+cmax+cdtheta+cdfi+ccratio)
4480 Int_t nentries = array->GetEntriesFast();
4481 AliHelix *helixes = new AliHelix[nentries];
4482 for (Int_t i=0;i<nentries;i++){
4483 AliTPCseed* track = (AliTPCseed*)array->At(i);
4484 if (!track) continue;
4485 track->SetCircular(0);
4486 new (&helixes[i]) AliHelix(*track);
4492 Double_t phase[2][2],radius[2];
4497 for (Int_t i0=0;i0<nentries;i0++){
4498 AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
4499 if (!track0) continue;
4500 if (TMath::Abs(track0->GetC())<1/kMaxC) continue;
4501 Float_t xc0 = helixes[i0].GetHelix(6);
4502 Float_t yc0 = helixes[i0].GetHelix(7);
4503 Float_t r0 = helixes[i0].GetHelix(8);
4504 Float_t rc0 = TMath::Sqrt(xc0*xc0+yc0*yc0);
4505 Float_t fi0 = TMath::ATan2(yc0,xc0);
4507 for (Int_t i1=i0+1;i1<nentries;i1++){
4508 AliTPCseed * track1 = (AliTPCseed*)array->At(i1);
4509 if (!track1) continue;
4510 if (TMath::Abs(track1->GetC())<1/kMaxC) continue;
4511 Float_t xc1 = helixes[i1].GetHelix(6);
4512 Float_t yc1 = helixes[i1].GetHelix(7);
4513 Float_t r1 = helixes[i1].GetHelix(8);
4514 Float_t rc1 = TMath::Sqrt(xc1*xc1+yc1*yc1);
4515 Float_t fi1 = TMath::ATan2(yc1,xc1);
4517 Float_t dfi = fi0-fi1;
4520 if (dfi>1.5*TMath::Pi()) dfi-=TMath::Pi(); // take care about edge effect
4521 if (dfi<-1.5*TMath::Pi()) dfi+=TMath::Pi(); //
4522 Float_t dtheta = TMath::Abs(track0->GetTgl()-track1->GetTgl())<TMath::Abs(track0->GetTgl()+track1->GetTgl())? track0->GetTgl()-track1->GetTgl():track0->GetTgl()+track1->GetTgl();
4526 if (track0->GetBConstrain()&&track1->GetBConstrain()) continue; // not constrained
4527 if (track1->GetSigned1Pt()*track0->GetSigned1Pt()>0) continue; // not the same sign
4528 if ( TMath::Abs(track1->GetTgl()+track0->GetTgl())>kMaxdTheta) continue; //distance in the Theta
4529 if ( TMath::Abs(dfi)>kMaxdPhi) continue; //distance in phi
4530 if ( TMath::Sqrt(dfi*dfi+dtheta*dtheta)>kMaxdPhi) continue; //common angular offset
4532 Float_t pt0 = track0->GetSignedPt();
4533 Float_t pt1 = track1->GetSignedPt();
4534 if ((TMath::Abs(pt0+pt1)/(TMath::Abs(pt0)+TMath::Abs(pt1)))>kPtRatio) continue;
4535 if ((iter==1) && TMath::Abs(TMath::Abs(rc0+r0)-TMath::Abs(rc1+r1))>kMaxDeltaRMax) continue;
4536 if ((iter!=1) &&TMath::Abs(TMath::Abs(rc0-r0)-TMath::Abs(rc1-r1))>kMaxDeltaRMin) continue;
4537 if (TMath::Min(TMath::Abs(rc0-r0),TMath::Abs(rc1-r1))<kMinDCAR) continue;
4540 // Now find closest approach
4544 Int_t npoints = helixes[i0].GetRPHIintersections(helixes[i1], phase, radius,10);
4545 if (npoints==0) continue;
4546 helixes[i0].GetClosestPhases(helixes[i1], phase);
4550 Double_t hangles[3];
4551 helixes[i0].Evaluate(phase[0][0],xyz0);
4552 helixes[i1].Evaluate(phase[0][1],xyz1);
4554 helixes[i0].GetAngle(phase[0][0],helixes[i1],phase[0][1],hangles);
4555 Double_t deltah[2],deltabest;
4556 if (TMath::Abs(hangles[2])<kMinAngle) continue;
4560 helixes[i0].ParabolicDCA(helixes[i1],phase[0][0],phase[0][1],radius[0],deltah[0],2);
4562 helixes[i0].ParabolicDCA(helixes[i1],phase[1][0],phase[1][1],radius[1],deltah[1],2);
4563 if (deltah[1]<deltah[0]) ibest=1;
4565 deltabest = TMath::Sqrt(deltah[ibest]);
4566 helixes[i0].Evaluate(phase[ibest][0],xyz0);
4567 helixes[i1].Evaluate(phase[ibest][1],xyz1);
4568 helixes[i0].GetAngle(phase[ibest][0],helixes[i1],phase[ibest][1],hangles);
4569 Double_t radiusbest = TMath::Sqrt(radius[ibest]);
4571 if (deltabest>kMaxDist) continue;
4572 // if (mindcar+mindcaz<40 && (TMath::Abs(hangles[2])<kMinAngle ||deltabest>3)) continue;
4573 Bool_t sign =kFALSE;
4574 if (hangles[2]>kMinAngle) sign =kTRUE;
4577 // circular[i0] = kTRUE;
4578 // circular[i1] = kTRUE;
4579 if (track0->OneOverPt()<track1->OneOverPt()){
4580 track0->SetCircular(track0->GetCircular()+1);
4581 track1->SetCircular(track1->GetCircular()+2);
4584 track1->SetCircular(track1->GetCircular()+1);
4585 track0->SetCircular(track0->GetCircular()+2);
4588 if (AliTPCReconstructor::StreamLevel()>1){
4590 //debug stream to tune "fine" cuts
4591 Int_t lab0=track0->GetLabel();
4592 Int_t lab1=track1->GetLabel();
4593 TTreeSRedirector &cstream = *fDebugStreamer;
4594 cstream<<"Curling2"<<
4610 "npoints="<<npoints<<
4611 "hangles0="<<hangles[0]<<
4612 "hangles1="<<hangles[1]<<
4613 "hangles2="<<hangles[2]<<
4616 "radius="<<radiusbest<<
4617 "deltabest="<<deltabest<<
4618 "phase0="<<phase[ibest][0]<<
4619 "phase1="<<phase[ibest][1]<<
4627 if (AliTPCReconstructor::StreamLevel()>1) {
4628 AliInfo("Time for curling tracks removal");
4637 void AliTPCtrackerMI::FindKinks(TObjArray * array, AliESDEvent *esd)
4644 TObjArray *kinks= new TObjArray(10000);
4645 // TObjArray *v0s= new TObjArray(10000);
4646 Int_t nentries = array->GetEntriesFast();
4647 AliHelix *helixes = new AliHelix[nentries];
4648 Int_t *sign = new Int_t[nentries];
4649 Int_t *nclusters = new Int_t[nentries];
4650 Float_t *alpha = new Float_t[nentries];
4651 AliKink *kink = new AliKink();
4652 Int_t * usage = new Int_t[nentries];
4653 Float_t *zm = new Float_t[nentries];
4654 Float_t *z0 = new Float_t[nentries];
4655 Float_t *fim = new Float_t[nentries];
4656 Float_t *shared = new Float_t[nentries];
4657 Bool_t *circular = new Bool_t[nentries];
4658 Float_t *dca = new Float_t[nentries];
4659 //const AliESDVertex * primvertex = esd->GetVertex();
4661 // nentries = array->GetEntriesFast();
4666 for (Int_t i=0;i<nentries;i++){
4669 AliTPCseed* track = (AliTPCseed*)array->At(i);
4670 if (!track) continue;
4671 track->SetCircular(0);
4673 track->UpdatePoints();
4674 if (( track->GetPoints()[2]- track->GetPoints()[0])>5 && track->GetPoints()[3]>0.8){
4676 nclusters[i]=track->GetNumberOfClusters();
4677 alpha[i] = track->GetAlpha();
4678 new (&helixes[i]) AliHelix(*track);
4680 helixes[i].Evaluate(0,xyz);
4681 sign[i] = (track->GetC()>0) ? -1:1;
4684 if (track->GetProlongation(x,y,z)){
4686 fim[i] = alpha[i]+TMath::ATan2(y,x);
4689 zm[i] = track->GetZ();
4693 circular[i]= kFALSE;
4694 if (track->GetProlongation(0,y,z)) z0[i] = z;
4695 dca[i] = track->GetD(0,0);
4701 Int_t ncandidates =0;
4704 Double_t phase[2][2],radius[2];
4707 // Find circling track
4709 for (Int_t i0=0;i0<nentries;i0++){
4710 AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
4711 if (!track0) continue;
4712 if (track0->GetNumberOfClusters()<40) continue;
4713 if (TMath::Abs(1./track0->GetC())>200) continue;
4714 for (Int_t i1=i0+1;i1<nentries;i1++){
4715 AliTPCseed * track1 = (AliTPCseed*)array->At(i1);
4716 if (!track1) continue;
4717 if (track1->GetNumberOfClusters()<40) continue;
4718 if ( TMath::Abs(track1->GetTgl()+track0->GetTgl())>0.1) continue;
4719 if (track0->GetBConstrain()&&track1->GetBConstrain()) continue;
4720 if (TMath::Abs(1./track1->GetC())>200) continue;
4721 if (track1->GetSigned1Pt()*track0->GetSigned1Pt()>0) continue;
4722 if (track1->GetTgl()*track0->GetTgl()>0) continue;
4723 if (TMath::Max(TMath::Abs(1./track0->GetC()),TMath::Abs(1./track1->GetC()))>190) continue;
4724 if (track0->GetBConstrain()&&track1->OneOverPt()<track0->OneOverPt()) continue; //returning - lower momenta
4725 if (track1->GetBConstrain()&&track0->OneOverPt()<track1->OneOverPt()) continue; //returning - lower momenta
4727 Float_t mindcar = TMath::Min(TMath::Abs(dca[i0]),TMath::Abs(dca[i1]));
4728 if (mindcar<5) continue;
4729 Float_t mindcaz = TMath::Min(TMath::Abs(z0[i0]-GetZ()),TMath::Abs(z0[i1]-GetZ()));
4730 if (mindcaz<5) continue;
4731 if (mindcar+mindcaz<20) continue;
4734 Float_t xc0 = helixes[i0].GetHelix(6);
4735 Float_t yc0 = helixes[i0].GetHelix(7);
4736 Float_t r0 = helixes[i0].GetHelix(8);
4737 Float_t xc1 = helixes[i1].GetHelix(6);
4738 Float_t yc1 = helixes[i1].GetHelix(7);
4739 Float_t r1 = helixes[i1].GetHelix(8);
4741 Float_t rmean = (r0+r1)*0.5;
4742 Float_t delta =TMath::Sqrt((xc1-xc0)*(xc1-xc0)+(yc1-yc0)*(yc1-yc0));
4743 //if (delta>30) continue;
4744 if (delta>rmean*0.25) continue;
4745 if (TMath::Abs(r0-r1)/rmean>0.3) continue;
4747 Int_t npoints = helixes[i0].GetRPHIintersections(helixes[i1], phase, radius,10);
4748 if (npoints==0) continue;
4749 helixes[i0].GetClosestPhases(helixes[i1], phase);
4753 Double_t hangles[3];
4754 helixes[i0].Evaluate(phase[0][0],xyz0);
4755 helixes[i1].Evaluate(phase[0][1],xyz1);
4757 helixes[i0].GetAngle(phase[0][0],helixes[i1],phase[0][1],hangles);
4758 Double_t deltah[2],deltabest;
4759 if (hangles[2]<2.8) continue;
4762 helixes[i0].ParabolicDCA(helixes[i1],phase[0][0],phase[0][1],radius[0],deltah[0],2);
4764 helixes[i0].ParabolicDCA(helixes[i1],phase[1][0],phase[1][1],radius[1],deltah[1],2);
4765 if (deltah[1]<deltah[0]) ibest=1;
4767 deltabest = TMath::Sqrt(deltah[ibest]);
4768 helixes[i0].Evaluate(phase[ibest][0],xyz0);
4769 helixes[i1].Evaluate(phase[ibest][1],xyz1);
4770 helixes[i0].GetAngle(phase[ibest][0],helixes[i1],phase[ibest][1],hangles);
4771 Double_t radiusbest = TMath::Sqrt(radius[ibest]);
4773 if (deltabest>6) continue;
4774 if (mindcar+mindcaz<40 && (hangles[2]<3.12||deltabest>3)) continue;
4775 Bool_t lsign =kFALSE;
4776 if (hangles[2]>3.06) lsign =kTRUE;
4779 circular[i0] = kTRUE;
4780 circular[i1] = kTRUE;
4781 if (track0->OneOverPt()<track1->OneOverPt()){
4782 track0->SetCircular(track0->GetCircular()+1);
4783 track1->SetCircular(track1->GetCircular()+2);
4786 track1->SetCircular(track1->GetCircular()+1);
4787 track0->SetCircular(track0->GetCircular()+2);
4790 if (lsign&&AliTPCReconstructor::StreamLevel()>1){
4792 Int_t lab0=track0->GetLabel();
4793 Int_t lab1=track1->GetLabel();
4794 TTreeSRedirector &cstream = *fDebugStreamer;
4795 cstream<<"Curling"<<
4802 "mindcar="<<mindcar<<
4803 "mindcaz="<<mindcaz<<
4806 "npoints="<<npoints<<
4807 "hangles0="<<hangles[0]<<
4808 "hangles2="<<hangles[2]<<
4813 "radius="<<radiusbest<<
4814 "deltabest="<<deltabest<<
4815 "phase0="<<phase[ibest][0]<<
4816 "phase1="<<phase[ibest][1]<<
4826 for (Int_t i =0;i<nentries;i++){
4827 if (sign[i]==0) continue;
4828 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
4835 Double_t cradius0 = 40*40;
4836 Double_t cradius1 = 270*270;
4839 Double_t cdist3=0.55;
4840 for (Int_t j =i+1;j<nentries;j++){
4842 if (sign[j]*sign[i]<1) continue;
4843 if ( (nclusters[i]+nclusters[j])>200) continue;
4844 if ( (nclusters[i]+nclusters[j])<80) continue;
4845 if ( TMath::Abs(zm[i]-zm[j])>60.) continue;
4846 if ( TMath::Abs(fim[i]-fim[j])>0.6 && TMath::Abs(fim[i]-fim[j])<5.7 ) continue;
4847 //AliTPCseed * track1 = (AliTPCseed*)array->At(j); Double_t phase[2][2],radius[2];
4848 Int_t npoints = helixes[i].GetRPHIintersections(helixes[j], phase, radius,20);
4849 if (npoints<1) continue;
4852 if (radius[0]<cradius0||radius[0]>cradius1) continue;
4855 if ( (radius[0]<cradius0||radius[0]>cradius1) && (radius[1]<cradius0||radius[1]>cradius1) ) continue;
4858 Double_t delta1=10000,delta2=10000;
4859 // cuts on the intersection radius
4860 helixes[i].LinearDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta1);
4861 if (radius[0]<20&&delta1<1) continue; //intersection at vertex
4862 if (radius[0]<10&&delta1<3) continue; //intersection at vertex
4864 helixes[i].LinearDCA(helixes[j],phase[1][0],phase[1][1],radius[1],delta2);
4865 if (radius[1]<20&&delta2<1) continue; //intersection at vertex
4866 if (radius[1]<10&&delta2<3) continue; //intersection at vertex
4869 Double_t distance1 = TMath::Min(delta1,delta2);
4870 if (distance1>cdist1) continue; // cut on DCA linear approximation
4872 npoints = helixes[i].GetRPHIintersections(helixes[j], phase, radius,20);
4873 helixes[i].ParabolicDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta1);
4874 if (radius[0]<20&&delta1<1) continue; //intersection at vertex
4875 if (radius[0]<10&&delta1<3) continue; //intersection at vertex
4878 helixes[i].ParabolicDCA(helixes[j],phase[1][0],phase[1][1],radius[1],delta2);
4879 if (radius[1]<20&&delta2<1) continue; //intersection at vertex
4880 if (radius[1]<10&&delta2<3) continue; //intersection at vertex
4882 distance1 = TMath::Min(delta1,delta2);
4885 rkink = TMath::Sqrt(radius[0]);
4888 rkink = TMath::Sqrt(radius[1]);
4890 if (distance1>cdist2) continue;
4893 AliTPCseed * track1 = (AliTPCseed*)array->At(j);
4896 Int_t row0 = GetRowNumber(rkink);
4897 if (row0<10) continue;
4898 if (row0>150) continue;
4901 Float_t dens00=-1,dens01=-1;
4902 Float_t dens10=-1,dens11=-1;
4904 Int_t found,foundable,ishared;
4905 track0->GetClusterStatistic(0,row0-5, found, foundable,ishared,kFALSE);
4906 if (foundable>5) dens00 = Float_t(found)/Float_t(foundable);
4907 track0->GetClusterStatistic(row0+5,155, found, foundable,ishared,kFALSE);
4908 if (foundable>5) dens01 = Float_t(found)/Float_t(foundable);
4910 track1->GetClusterStatistic(0,row0-5, found, foundable,ishared,kFALSE);
4911 if (foundable>10) dens10 = Float_t(found)/Float_t(foundable);
4912 track1->GetClusterStatistic(row0+5,155, found, foundable,ishared,kFALSE);
4913 if (foundable>10) dens11 = Float_t(found)/Float_t(foundable);
4915 if (dens00<dens10 && dens01<dens11) continue;
4916 if (dens00>dens10 && dens01>dens11) continue;
4917 if (TMath::Max(dens00,dens10)<0.1) continue;
4918 if (TMath::Max(dens01,dens11)<0.3) continue;
4920 if (TMath::Min(dens00,dens10)>0.6) continue;
4921 if (TMath::Min(dens01,dens11)>0.6) continue;
4924 AliTPCseed * ktrack0, *ktrack1;
4933 if (TMath::Abs(ktrack0->GetC())>5) continue; // cut on the curvature for mother particle
4934 AliExternalTrackParam paramm(*ktrack0);
4935 AliExternalTrackParam paramd(*ktrack1);
4936 if (row0>60&&ktrack1->GetReference().GetX()>90.)new (¶md) AliExternalTrackParam(ktrack1->GetReference());
4939 kink->SetMother(paramm);
4940 kink->SetDaughter(paramd);
4943 Float_t x[3] = { kink->GetPosition()[0],kink->GetPosition()[1],kink->GetPosition()[2]};
4945 fParam->Transform0to1(x,index);
4946 fParam->Transform1to2(x,index);
4947 row0 = GetRowNumber(x[0]);
4949 if (kink->GetR()<100) continue;
4950 if (kink->GetR()>240) continue;
4951 if (kink->GetPosition()[2]/kink->GetR()>AliTPCReconstructor::GetCtgRange()) continue; //out of fiducial volume
4952 if (kink->GetDistance()>cdist3) continue;
4953 Float_t dird = kink->GetDaughterP()[0]*kink->GetPosition()[0]+kink->GetDaughterP()[1]*kink->GetPosition()[1]; // rough direction estimate
4954 if (dird<0) continue;
4956 Float_t dirm = kink->GetMotherP()[0]*kink->GetPosition()[0]+kink->GetMotherP()[1]*kink->GetPosition()[1]; // rough direction estimate
4957 if (dirm<0) continue;
4958 Float_t mpt = TMath::Sqrt(kink->GetMotherP()[0]*kink->GetMotherP()[0]+kink->GetMotherP()[1]*kink->GetMotherP()[1]);
4959 if (mpt<0.2) continue;
4962 //for high momenta momentum not defined well in first iteration
4963 Double_t qt = TMath::Sin(kink->GetAngle(2))*ktrack1->GetP();
4964 if (qt>0.35) continue;
4967 kink->SetLabel(CookLabel(ktrack0,0.4,0,row0),0);
4968 kink->SetLabel(CookLabel(ktrack1,0.4,row0,160),1);
4970 kink->SetTPCDensity(dens00,0,0);
4971 kink->SetTPCDensity(dens01,0,1);
4972 kink->SetTPCDensity(dens10,1,0);
4973 kink->SetTPCDensity(dens11,1,1);
4974 kink->SetIndex(i,0);
4975 kink->SetIndex(j,1);
4978 kink->SetTPCDensity(dens10,0,0);
4979 kink->SetTPCDensity(dens11,0,1);
4980 kink->SetTPCDensity(dens00,1,0);
4981 kink->SetTPCDensity(dens01,1,1);
4982 kink->SetIndex(j,0);
4983 kink->SetIndex(i,1);
4986 if (mpt<1||kink->GetAngle(2)>0.1){
4987 // angle and densities not defined yet
4988 if (kink->GetTPCDensityFactor()<0.8) continue;
4989 if ((2-kink->GetTPCDensityFactor())*kink->GetDistance() >0.25) continue;
4990 if (kink->GetAngle(2)*ktrack0->GetP()<0.003) continue; //too small angle
4991 if (kink->GetAngle(2)>0.2&&kink->GetTPCDensityFactor()<1.15) continue;
4992 if (kink->GetAngle(2)>0.2&&kink->GetTPCDensity(0,1)>0.05) continue;
4994 Float_t criticalangle = track0->GetSigmaSnp2()+track0->GetSigmaTgl2();
4995 criticalangle+= track1->GetSigmaSnp2()+track1->GetSigmaTgl2();
4996 criticalangle= 3*TMath::Sqrt(criticalangle);
4997 if (criticalangle>0.02) criticalangle=0.02;
4998 if (kink->GetAngle(2)<criticalangle) continue;
5001 Int_t drow = Int_t(2.+0.5/(0.05+kink->GetAngle(2))); // overlap region defined
5002 Float_t shapesum =0;
5004 for ( Int_t row = row0-drow; row<row0+drow;row++){
5005 if (row<0) continue;
5006 if (row>155) continue;
5007 if (ktrack0->GetClusterPointer(row)){
5008 AliTPCTrackerPoint *point =ktrack0->GetTrackPoint(row);
5009 shapesum+=point->GetSigmaY()+point->GetSigmaZ();
5012 if (ktrack1->GetClusterPointer(row)){
5013 AliTPCTrackerPoint *point =ktrack1->GetTrackPoint(row);
5014 shapesum+=point->GetSigmaY()+point->GetSigmaZ();
5019 kink->SetShapeFactor(-1.);
5022 kink->SetShapeFactor(shapesum/sum);
5024 // esd->AddKink(kink);
5025 kinks->AddLast(kink);
5031 // sort the kinks according quality - and refit them towards vertex
5033 Int_t nkinks = kinks->GetEntriesFast();
5034 Float_t *quality = new Float_t[nkinks];
5035 Int_t *indexes = new Int_t[nkinks];
5036 AliTPCseed *mothers = new AliTPCseed[nkinks];
5037 AliTPCseed *daughters = new AliTPCseed[nkinks];
5040 for (Int_t i=0;i<nkinks;i++){
5042 AliKink *kinkl = (AliKink*)kinks->At(i);
5044 // refit kinks towards vertex
5046 Int_t index0 = kinkl->GetIndex(0);
5047 Int_t index1 = kinkl->GetIndex(1);
5048 AliTPCseed * ktrack0 = (AliTPCseed*)array->At(index0);
5049 AliTPCseed * ktrack1 = (AliTPCseed*)array->At(index1);
5051 Int_t sumn=ktrack0->GetNumberOfClusters()+ktrack1->GetNumberOfClusters();
5053 // Refit Kink under if too small angle
5055 if (kinkl->GetAngle(2)<0.05){
5056 kinkl->SetTPCRow0(GetRowNumber(kinkl->GetR()));
5057 Int_t row0 = kinkl->GetTPCRow0();
5058 Int_t drow = Int_t(2.+0.5/(0.05+kinkl->GetAngle(2)));
5061 Int_t last = row0-drow;
5062 if (last<40) last=40;
5063 if (last<ktrack0->GetFirstPoint()+25) last = ktrack0->GetFirstPoint()+25;
5064 AliTPCseed* seed0 = ReSeed(ktrack0,last,kFALSE);
5067 Int_t first = row0+drow;
5068 if (first>130) first=130;
5069 if (first>ktrack1->GetLastPoint()-25) first = TMath::Max(ktrack1->GetLastPoint()-25,30);
5070 AliTPCseed* seed1 = ReSeed(ktrack1,first,kTRUE);
5072 if (seed0 && seed1){
5073 kinkl->SetStatus(1,8);
5074 if (RefitKink(*seed0,*seed1,*kinkl)) kinkl->SetStatus(1,9);
5075 row0 = GetRowNumber(kinkl->GetR());
5076 sumn = seed0->GetNumberOfClusters()+seed1->GetNumberOfClusters();
5077 mothers[i] = *seed0;
5078 daughters[i] = *seed1;
5081 delete kinks->RemoveAt(i);
5082 if (seed0) delete seed0;
5083 if (seed1) delete seed1;
5086 if (kinkl->GetDistance()>0.5 || kinkl->GetR()<110 || kinkl->GetR()>240) {
5087 delete kinks->RemoveAt(i);
5088 if (seed0) delete seed0;
5089 if (seed1) delete seed1;
5097 if (kinkl) quality[i] = 160*((0.1+kinkl->GetDistance())*(2.-kinkl->GetTPCDensityFactor()))/(sumn+40.); //the longest -clossest will win
5099 TMath::Sort(nkinks,quality,indexes,kFALSE);
5101 //remove double find kinks
5103 for (Int_t ikink0=1;ikink0<nkinks;ikink0++){
5104 AliKink * kink0 = (AliKink*) kinks->At(indexes[ikink0]);
5105 if (!kink0) continue;
5107 for (Int_t ikink1=0;ikink1<ikink0;ikink1++){
5108 if (!kink0) continue;
5109 AliKink * kink1 = (AliKink*) kinks->At(indexes[ikink1]);
5110 if (!kink1) continue;
5111 // if not close kink continue
5112 if (TMath::Abs(kink1->GetPosition()[2]-kink0->GetPosition()[2])>10) continue;
5113 if (TMath::Abs(kink1->GetPosition()[1]-kink0->GetPosition()[1])>10) continue;
5114 if (TMath::Abs(kink1->GetPosition()[0]-kink0->GetPosition()[0])>10) continue;
5116 AliTPCseed &mother0 = mothers[indexes[ikink0]];
5117 AliTPCseed &daughter0 = daughters[indexes[ikink0]];
5118 AliTPCseed &mother1 = mothers[indexes[ikink1]];
5119 AliTPCseed &daughter1 = daughters[indexes[ikink1]];
5120 Int_t row0 = (kink0->GetTPCRow0()+kink1->GetTPCRow0())/2;
5129 for (Int_t i=0;i<row0;i++){
5130 if (mother0.GetClusterIndex(i)>0 && mother1.GetClusterIndex(i)>0){
5133 if (mother0.GetClusterIndex(i)==mother1.GetClusterIndex(i)){
5140 for (Int_t i=row0;i<158;i++){
5141 if (daughter0.GetClusterIndex(i)>0 && daughter0.GetClusterIndex(i)>0){
5144 if (mother0.GetClusterIndex(i)==mother1.GetClusterIndex(i)){
5150 Float_t ratio = Float_t(same+1)/Float_t(both+1);
5151 Float_t ratiom = Float_t(samem+1)/Float_t(bothm+1);
5152 Float_t ratiod = Float_t(samed+1)/Float_t(bothd+1);
5153 if (ratio>0.3 && ratiom>0.5 &&ratiod>0.5) {
5154 Int_t sum0 = mother0.GetNumberOfClusters()+daughter0.GetNumberOfClusters();
5155 Int_t sum1 = mother1.GetNumberOfClusters()+daughter1.GetNumberOfClusters();
5157 shared[kink0->GetIndex(0)]= kTRUE;
5158 shared[kink0->GetIndex(1)]= kTRUE;
5159 delete kinks->RemoveAt(indexes[ikink0]);
5162 shared[kink1->GetIndex(0)]= kTRUE;
5163 shared[kink1->GetIndex(1)]= kTRUE;
5164 delete kinks->RemoveAt(indexes[ikink1]);
5171 for (Int_t i=0;i<nkinks;i++){
5172 AliKink * kinkl = (AliKink*) kinks->At(indexes[i]);
5173 if (!kinkl) continue;
5174 kinkl->SetTPCRow0(GetRowNumber(kinkl->GetR()));
5175 Int_t index0 = kinkl->GetIndex(0);
5176 Int_t index1 = kinkl->GetIndex(1);
5177 if (circular[index0]||(circular[index1]&&kinkl->GetDistance()>0.2)) continue;
5178 kinkl->SetMultiple(usage[index0],0);
5179 kinkl->SetMultiple(usage[index1],1);
5180 if (kinkl->GetMultiple()[0]+kinkl->GetMultiple()[1]>2) continue;
5181 if (kinkl->GetMultiple()[0]+kinkl->GetMultiple()[1]>0 && quality[indexes[i]]>0.2) continue;
5182 if (kinkl->GetMultiple()[0]+kinkl->GetMultiple()[1]>0 && kinkl->GetDistance()>0.2) continue;
5183 if (circular[index0]||(circular[index1]&&kinkl->GetDistance()>0.1)) continue;
5185 AliTPCseed * ktrack0 = (AliTPCseed*)array->At(index0);
5186 AliTPCseed * ktrack1 = (AliTPCseed*)array->At(index1);
5187 if (!ktrack0 || !ktrack1) continue;
5188 Int_t index = esd->AddKink(kinkl);
5191 if ( ktrack0->GetKinkIndex(0)==0 && ktrack1->GetKinkIndex(0)==0) { //best kink
5192 if (mothers[indexes[i]].GetNumberOfClusters()>20 && daughters[indexes[i]].GetNumberOfClusters()>20 && (mothers[indexes[i]].GetNumberOfClusters()+daughters[indexes[i]].GetNumberOfClusters())>100){
5193 *ktrack0 = mothers[indexes[i]];
5194 *ktrack1 = daughters[indexes[i]];
5198 ktrack0->SetKinkIndex(usage[index0],-(index+1));
5199 ktrack1->SetKinkIndex(usage[index1], (index+1));
5204 // Remove tracks corresponding to shared kink's
5206 for (Int_t i=0;i<nentries;i++){
5207 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5208 if (!track0) continue;
5209 if (track0->GetKinkIndex(0)!=0) continue;
5210 if (shared[i]) delete array->RemoveAt(i);
5215 RemoveUsed2(array,0.5,0.4,30);
5217 for (Int_t i=0;i<nentries;i++){
5218 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5219 if (!track0) continue;
5220 track0->CookdEdx(0.02,0.6);
5224 for (Int_t i=0;i<nentries;i++){
5225 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5226 if (!track0) continue;
5227 if (track0->Pt()<1.4) continue;
5228 //remove double high momenta tracks - overlapped with kink candidates
5231 for (Int_t icl=track0->GetFirstPoint();icl<track0->GetLastPoint(); icl++){
5232 if (track0->GetClusterPointer(icl)!=0){
5234 if (track0->GetClusterPointer(icl)->IsUsed(10)) ishared++;
5237 if (Float_t(ishared+1)/Float_t(all+1)>0.5) {
5238 delete array->RemoveAt(i);
5242 if (track0->GetKinkIndex(0)!=0) continue;
5243 if (track0->GetNumberOfClusters()<80) continue;
5245 AliTPCseed *pmother = new AliTPCseed();
5246 AliTPCseed *pdaughter = new AliTPCseed();
5247 AliKink *pkink = new AliKink;
5249 AliTPCseed & mother = *pmother;
5250 AliTPCseed & daughter = *pdaughter;
5251 AliKink & kinkl = *pkink;
5252 if (CheckKinkPoint(track0,mother,daughter, kinkl)){
5253 if (mother.GetNumberOfClusters()<30||daughter.GetNumberOfClusters()<20) {
5257 continue; //too short tracks
5259 if (mother.Pt()<1.4) {
5265 Int_t row0= kinkl.GetTPCRow0();
5266 if (kinkl.GetDistance()>0.5 || kinkl.GetR()<110. || kinkl.GetR()>240.) {
5273 Int_t index = esd->AddKink(&kinkl);
5274 mother.SetKinkIndex(0,-(index+1));
5275 daughter.SetKinkIndex(0,index+1);
5276 if (mother.GetNumberOfClusters()>50) {
5277 delete array->RemoveAt(i);
5278 array->AddAt(new AliTPCseed(mother),i);
5281 array->AddLast(new AliTPCseed(mother));
5283 array->AddLast(new AliTPCseed(daughter));
5284 for (Int_t icl=0;icl<row0;icl++) {
5285 if (mother.GetClusterPointer(icl)) mother.GetClusterPointer(icl)->Use(20);
5288 for (Int_t icl=row0;icl<158;icl++) {
5289 if (daughter.GetClusterPointer(icl)) daughter.GetClusterPointer(icl)->Use(20);
5298 delete [] daughters;
5320 printf("Ncandidates=\t%d\t%d\t%d\t%d\n",esd->GetNumberOfKinks(),ncandidates,ntracks,nall);
5324 void AliTPCtrackerMI::FindV0s(TObjArray * array, AliESDEvent *esd)
5330 TObjArray *tpcv0s = new TObjArray(100000);
5331 Int_t nentries = array->GetEntriesFast();
5332 AliHelix *helixes = new AliHelix[nentries];
5333 Int_t *sign = new Int_t[nentries];
5334 Float_t *alpha = new Float_t[nentries];
5335 Float_t *z0 = new Float_t[nentries];
5336 Float_t *dca = new Float_t[nentries];
5337 Float_t *sdcar = new Float_t[nentries];
5338 Float_t *cdcar = new Float_t[nentries];
5339 Float_t *pulldcar = new Float_t[nentries];
5340 Float_t *pulldcaz = new Float_t[nentries];
5341 Float_t *pulldca = new Float_t[nentries];
5342 Bool_t *isPrim = new Bool_t[nentries];
5343 const AliESDVertex * primvertex = esd->GetVertex();
5344 Double_t zvertex = primvertex->GetZv();
5346 // nentries = array->GetEntriesFast();
5348 for (Int_t i=0;i<nentries;i++){
5351 AliTPCseed* track = (AliTPCseed*)array->At(i);
5352 if (!track) continue;
5353 track->GetV0Indexes()[0] = 0; //rest v0 indexes
5354 track->GetV0Indexes()[1] = 0; //rest v0 indexes
5355 track->GetV0Indexes()[2] = 0; //rest v0 indexes
5357 alpha[i] = track->GetAlpha();
5358 new (&helixes[i]) AliHelix(*track);
5360 helixes[i].Evaluate(0,xyz);
5361 sign[i] = (track->GetC()>0) ? -1:1;
5365 if (track->GetProlongation(0,y,z)) z0[i] = z;
5366 dca[i] = track->GetD(0,0);
5368 // dca error parrameterezation + pulls
5370 sdcar[i] = TMath::Sqrt(0.150*0.150+(100*track->GetC())*(100*track->GetC()));
5371 if (TMath::Abs(track->GetTgl())>1) sdcar[i]*=2.5;
5372 cdcar[i] = TMath::Exp((TMath::Abs(track->GetC())-0.0106)*525.3);
5373 pulldcar[i] = (dca[i]-cdcar[i])/sdcar[i];
5374 pulldcaz[i] = (z0[i]-zvertex)/sdcar[i];
5375 pulldca[i] = TMath::Sqrt(pulldcar[i]*pulldcar[i]+pulldcaz[i]*pulldcaz[i]);
5376 if (track->TPCrPID(1)+track->TPCrPID(2)+track->TPCrPID(3)>0.5) {
5377 if (pulldca[i]<3.) isPrim[i]=kTRUE; //pion, muon and Kaon 3 sigma cut
5379 if (track->TPCrPID(4)>0.5) {
5380 if (pulldca[i]<0.5) isPrim[i]=kTRUE; //proton 0.5 sigma cut
5382 if (track->TPCrPID(0)>0.4) {
5383 isPrim[i]=kFALSE; //electron no sigma cut
5390 Int_t ncandidates =0;
5393 Double_t phase[2][2],radius[2];
5399 Float_t fprimvertex[3]={GetX(),GetY(),GetZ()};
5401 Double_t cradius0 = 10*10;
5402 Double_t cradius1 = 200*200;
5405 Double_t cpointAngle = 0.95;
5407 Double_t delta[2]={10000,10000};
5408 for (Int_t i =0;i<nentries;i++){
5409 if (sign[i]==0) continue;
5410 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5411 if (!track0) continue;
5412 if (AliTPCReconstructor::StreamLevel()>1){
5413 TTreeSRedirector &cstream = *fDebugStreamer;
5418 "zvertex="<<zvertex<<
5419 "sdcar0="<<sdcar[i]<<
5420 "cdcar0="<<cdcar[i]<<
5421 "pulldcar0="<<pulldcar[i]<<
5422 "pulldcaz0="<<pulldcaz[i]<<
5423 "pulldca0="<<pulldca[i]<<
5424 "isPrim="<<isPrim[i]<<
5428 if (track0->GetSigned1Pt()<0) continue;
5429 if (track0->GetKinkIndex(0)>0||isPrim[i]) continue; //daughter kink
5431 if (TMath::Abs(helixes[i].GetHelix(4))<0.000000001) continue;
5436 for (Int_t j =0;j<nentries;j++){
5437 AliTPCseed * track1 = (AliTPCseed*)array->At(j);
5438 if (!track1) continue;
5439 if (track1->GetKinkIndex(0)>0 || isPrim[j]) continue; //daughter kink
5440 if (sign[j]*sign[i]>0) continue;
5441 if (TMath::Abs(helixes[j].GetHelix(4))<0.000001) continue;
5442 if (track0->GetCircular()+track1->GetCircular()>1) continue; //circling -returning track
5445 // DCA to prim vertex cut
5451 Int_t npoints = helixes[i].GetRPHIintersections(helixes[j], phase, radius,cdist2);
5452 if (npoints<1) continue;
5456 if (radius[0]<cradius0||radius[0]>cradius1) continue;
5457 helixes[i].LinearDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta[0]);
5458 if (delta[0]>cdist1) continue;
5461 if (TMath::Max(radius[0],radius[1])<cradius0|| TMath::Min(radius[0],radius[1])>cradius1) continue;
5462 helixes[i].LinearDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta[0]);
5463 helixes[i].LinearDCA(helixes[j],phase[1][0],phase[1][1],radius[1],delta[1]);
5464 if (delta[1]<delta[0]) iclosest=1;
5465 if (delta[iclosest]>cdist1) continue;
5467 helixes[i].ParabolicDCA(helixes[j],phase[iclosest][0],phase[iclosest][1],radius[iclosest],delta[iclosest]);
5468 if (radius[iclosest]<cradius0 || radius[iclosest]>cradius1 || delta[iclosest]>cdist1) continue;
5470 Double_t pointAngle = helixes[i].GetPointAngle(helixes[j],phase[iclosest],fprimvertex);
5471 if (pointAngle<cpointAngle) continue;
5473 Bool_t isGamma = kFALSE;
5474 vertex.SetParamP(*track0); //track0 - plus
5475 vertex.SetParamN(*track1); //track1 - minus
5476 vertex.Update(fprimvertex);
5477 if (track0->TPCrPID(0)>0.3&&track1->TPCrPID(0)>0.3&&vertex.GetAnglep()[2]<0.15) isGamma=kTRUE; // gamma conversion candidate
5478 Double_t pointAngle2 = vertex.GetV0CosineOfPointingAngle();
5480 if (vertex.GetV0CosineOfPointingAngle()<cpointAngle && (!isGamma)) continue;// point angle cut
5481 //Bo: if (vertex.GetDist2()>2&&(!isGamma)) continue; // point angle cut
5482 if (vertex.GetDcaV0Daughters()>2&&(!isGamma)) continue;//Bo: // point angle cut
5483 Float_t sigmae = 0.15*0.15;
5484 if (vertex.GetRr()<80)
5485 sigmae += (sdcar[i]*sdcar[i]+sdcar[j]*sdcar[j])*(1.-vertex.GetRr()/80.)*(1.-vertex.GetRr()/80.);
5486 sigmae+= TMath::Sqrt(sigmae);
5487 //Bo: if (vertex.GetDist2()/sigmae>3.&&(!isGamma)) continue;
5488 if (vertex.GetDcaV0Daughters()/sigmae>3.&&(!isGamma)) continue;
5489 Float_t densb0=0,densb1=0,densa0=0,densa1=0;
5490 Int_t row0 = GetRowNumber(vertex.GetRr());
5492 //Bo: if (vertex.GetDist2()>0.2) continue;
5493 if (vertex.GetDcaV0Daughters()>0.2) continue;
5494 densb0 = track0->Density2(0,row0-5);
5495 densb1 = track1->Density2(0,row0-5);
5496 if (densb0>0.3|| densb1>0.3) continue; //clusters before vertex
5497 densa0 = track0->Density2(row0+5,row0+40);
5498 densa1 = track1->Density2(row0+5,row0+40);
5499 if ((densa0<0.4|| densa1<0.4)&&(!isGamma)) continue; //missing clusters after vertex
5502 densa0 = track0->Density2(0,40); //cluster density
5503 densa1 = track1->Density2(0,40); //cluster density
5504 if ((vertex.GetRr()<80&&densa0+densa1<1.)&&(!isGamma)) continue;
5506 //Bo: vertex.SetLab(0,track0->GetLabel());
5507 //Bo: vertex.SetLab(1,track1->GetLabel());
5508 vertex.SetChi2After((densa0+densa1)*0.5);
5509 vertex.SetChi2Before((densb0+densb1)*0.5);
5510 vertex.SetIndex(0,i);
5511 vertex.SetIndex(1,j);
5512 //Bo: vertex.SetStatus(1); // TPC v0 candidate
5513 vertex.SetOnFlyStatus(2);//Bo: // TPC v0 candidate
5514 //Bo: vertex.SetRp(track0->TPCrPIDs());
5515 //Bo: vertex.SetRm(track1->TPCrPIDs());
5516 tpcv0s->AddLast(new AliESDv0(vertex));
5519 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
5520 Double_t radiusm= (delta[0]<delta[1])? TMath::Sqrt(radius[0]):TMath::Sqrt(radius[1]);
5521 Double_t deltam= (delta[0]<delta[1])? TMath::Sqrt(delta[0]):TMath::Sqrt(delta[1]);
5522 if (AliTPCReconstructor::StreamLevel()>1) {
5523 Int_t lab0=track0->GetLabel();
5524 Int_t lab1=track1->GetLabel();
5525 Char_t c0=track0->GetCircular();
5526 Char_t c1=track1->GetCircular();
5527 TTreeSRedirector &cstream = *fDebugStreamer;
5530 "vertex.="<<&vertex<<
5533 "Helix0.="<<&helixes[i]<<
5536 "Helix1.="<<&helixes[j]<<
5537 "pointAngle="<<pointAngle<<
5538 "pointAngle2="<<pointAngle2<<
5543 "zvertex="<<zvertex<<
5546 "npoints="<<npoints<<
5547 "radius0="<<radius[0]<<
5548 "delta0="<<delta[0]<<
5549 "radius1="<<radius[1]<<
5550 "delta1="<<delta[1]<<
5551 "radiusm="<<radiusm<<
5553 "sdcar0="<<sdcar[i]<<
5554 "sdcar1="<<sdcar[j]<<
5555 "cdcar0="<<cdcar[i]<<
5556 "cdcar1="<<cdcar[j]<<
5557 "pulldcar0="<<pulldcar[i]<<
5558 "pulldcar1="<<pulldcar[j]<<
5559 "pulldcaz0="<<pulldcaz[i]<<
5560 "pulldcaz1="<<pulldcaz[j]<<
5561 "pulldca0="<<pulldca[i]<<
5562 "pulldca1="<<pulldca[j]<<
5572 Float_t *quality = new Float_t[ncandidates];
5573 Int_t *indexes = new Int_t[ncandidates];
5575 for (Int_t i=0;i<ncandidates;i++){
5577 AliESDv0 *v0 = (AliESDv0*)tpcv0s->At(i);
5578 quality[i] = 1./(1.00001-v0->GetV0CosineOfPointingAngle()); //base point angle
5579 // quality[i] /= (0.5+v0->GetDist2());
5580 // quality[i] *= v0->GetChi2After(); //density factor
5582 Int_t index0 = v0->GetIndex(0);
5583 Int_t index1 = v0->GetIndex(1);
5584 //Bo: Double_t minpulldca = TMath::Min(2.+pulldca[v0->GetIndex(0)],(2.+pulldca[v0->GetIndex(1)]) ); //pull
5585 Double_t minpulldca = TMath::Min(2.+pulldca[index0],(2.+pulldca[index1]) );//Bo:
5589 AliTPCseed * track0 = (AliTPCseed*)array->At(index0);
5590 AliTPCseed * track1 = (AliTPCseed*)array->At(index1);
5591 if (track0->TPCrPID(0)>0.3&&track1->TPCrPID(0)>0.3&&v0->GetAnglep()[2]<0.15) quality[i]+=1000000; // gamma conversion candidate
5592 if (track0->TPCrPID(4)>0.9||(track1->TPCrPID(4)>0.9&&minpulldca>4)) quality[i]*=10; // lambda candidate candidate
5595 TMath::Sort(ncandidates,quality,indexes,kTRUE);
5598 for (Int_t i=0;i<ncandidates;i++){
5599 AliESDv0 * v0 = (AliESDv0*)tpcv0s->At(indexes[i]);
5601 Int_t index0 = v0->GetIndex(0);
5602 Int_t index1 = v0->GetIndex(1);
5603 AliTPCseed * track0 = (AliTPCseed*)array->At(index0);
5604 AliTPCseed * track1 = (AliTPCseed*)array->At(index1);
5605 if (!track0||!track1) {
5609 Bool_t accept =kTRUE; //default accept
5610 Int_t *v0indexes0 = track0->GetV0Indexes();
5611 Int_t *v0indexes1 = track1->GetV0Indexes();
5613 Int_t order0 = (v0indexes0[0]!=0) ? 1:0;
5614 Int_t order1 = (v0indexes1[0]!=0) ? 1:0;
5615 if (v0indexes0[1]!=0) order0 =2;
5616 if (v0indexes1[1]!=0) order1 =2;
5618 if (v0indexes0[2]!=0) {order0=3; accept=kFALSE;}
5619 if (v0indexes0[2]!=0) {order1=3; accept=kFALSE;}
5621 AliESDv0 * v02 = v0;
5623 //Bo: v0->SetOrder(0,order0);
5624 //Bo: v0->SetOrder(1,order1);
5625 //Bo: v0->SetOrder(1,order0+order1);
5626 v0->SetOnFlyStatus(kTRUE);
5627 Int_t index = esd->AddV0(v0);
5628 v02 = esd->GetV0(index);
5629 v0indexes0[order0]=index;
5630 v0indexes1[order1]=index;
5634 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
5635 if (AliTPCReconstructor::StreamLevel()>1) {
5636 Int_t lab0=track0->GetLabel();
5637 Int_t lab1=track1->GetLabel();
5638 TTreeSRedirector &cstream = *fDebugStreamer;
5647 "dca0="<<dca[index0]<<
5648 "dca1="<<dca[index1]<<
5652 "quality="<<quality[i]<<
5653 "pulldca0="<<pulldca[index0]<<
5654 "pulldca1="<<pulldca[index1]<<
5678 printf("TPC V0 finder : naccepted\t%d\tncandidates\t%d\tntracks\t%d\tnall\t%d\n",naccepted,ncandidates,ntracks,nall);
5682 Int_t AliTPCtrackerMI::RefitKink(AliTPCseed &mother, AliTPCseed &daughter, AliESDkink &knk)
5685 // refit kink towards to the vertex
5688 AliKink &kink=(AliKink &)knk;
5690 Int_t row0 = GetRowNumber(kink.GetR());
5691 FollowProlongation(mother,0);
5692 mother.Reset(kFALSE);
5694 FollowProlongation(daughter,row0);
5695 daughter.Reset(kFALSE);
5696 FollowBackProlongation(daughter,158);
5697 daughter.Reset(kFALSE);
5698 Int_t first = TMath::Max(row0-20,30);
5699 Int_t last = TMath::Min(row0+20,140);
5701 const Int_t kNdiv =5;
5702 AliTPCseed param0[kNdiv]; // parameters along the track
5703 AliTPCseed param1[kNdiv]; // parameters along the track
5704 AliKink kinks[kNdiv]; // corresponding kink parameters
5707 for (Int_t irow=0; irow<kNdiv;irow++){
5708 rows[irow] = first +((last-first)*irow)/(kNdiv-1);
5710 // store parameters along the track
5712 for (Int_t irow=0;irow<kNdiv;irow++){
5713 FollowBackProlongation(mother, rows[irow]);
5714 FollowProlongation(daughter,rows[kNdiv-1-irow]);
5715 param0[irow] = mother;
5716 param1[kNdiv-1-irow] = daughter;
5720 for (Int_t irow=0; irow<kNdiv-1;irow++){
5721 if (param0[irow].GetNumberOfClusters()<kNdiv||param1[irow].GetNumberOfClusters()<kNdiv) continue;
5722 kinks[irow].SetMother(param0[irow]);
5723 kinks[irow].SetDaughter(param1[irow]);
5724 kinks[irow].Update();
5727 // choose kink with best "quality"
5729 Double_t mindist = 10000;
5730 for (Int_t irow=0;irow<kNdiv;irow++){
5731 if (param0[irow].GetNumberOfClusters()<20||param1[irow].GetNumberOfClusters()<20) continue;
5732 if (TMath::Abs(kinks[irow].GetR())>240.) continue;
5733 if (TMath::Abs(kinks[irow].GetR())<100.) continue;
5735 Float_t normdist = TMath::Abs(param0[irow].GetX()-kinks[irow].GetR())*(0.1+kink.GetDistance());
5736 normdist/= (param0[irow].GetNumberOfClusters()+param1[irow].GetNumberOfClusters()+40.);
5737 if (normdist < mindist){
5743 if (index==-1) return 0;
5746 param0[index].Reset(kTRUE);
5747 FollowProlongation(param0[index],0);
5749 mother = param0[index];
5750 daughter = param1[index]; // daughter in vertex
5752 kink.SetMother(mother);
5753 kink.SetDaughter(daughter);
5755 kink.SetTPCRow0(GetRowNumber(kink.GetR()));
5756 kink.SetTPCncls(param0[index].GetNumberOfClusters(),0);
5757 kink.SetTPCncls(param1[index].GetNumberOfClusters(),1);
5758 kink.SetLabel(CookLabel(&mother,0.4, 0,kink.GetTPCRow0()),0);
5759 kink.SetLabel(CookLabel(&daughter,0.4, kink.GetTPCRow0(),160),1);
5760 mother.SetLabel(kink.GetLabel(0));
5761 daughter.SetLabel(kink.GetLabel(1));
5767 void AliTPCtrackerMI::UpdateKinkQualityM(AliTPCseed * seed){
5769 // update Kink quality information for mother after back propagation
5771 if (seed->GetKinkIndex(0)>=0) return;
5772 for (Int_t ikink=0;ikink<3;ikink++){
5773 Int_t index = seed->GetKinkIndex(ikink);
5774 if (index>=0) break;
5775 index = TMath::Abs(index)-1;
5776 AliESDkink * kink = fEvent->GetKink(index);
5777 kink->SetTPCDensity(-1,0,0);
5778 kink->SetTPCDensity(1,0,1);
5780 Int_t row0 = kink->GetTPCRow0() - 2 - Int_t( 0.5/ (0.05+kink->GetAngle(2)));
5781 if (row0<15) row0=15;
5783 Int_t row1 = kink->GetTPCRow0() + 2 + Int_t( 0.5/ (0.05+kink->GetAngle(2)));
5784 if (row1>145) row1=145;
5786 Int_t found,foundable,shared;
5787 seed->GetClusterStatistic(0,row0, found, foundable,shared,kFALSE);
5788 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),0,0);
5789 seed->GetClusterStatistic(row1,155, found, foundable,shared,kFALSE);
5790 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),0,1);
5795 void AliTPCtrackerMI::UpdateKinkQualityD(AliTPCseed * seed){
5797 // update Kink quality information for daughter after refit
5799 if (seed->GetKinkIndex(0)<=0) return;
5800 for (Int_t ikink=0;ikink<3;ikink++){
5801 Int_t index = seed->GetKinkIndex(ikink);
5802 if (index<=0) break;
5803 index = TMath::Abs(index)-1;
5804 AliESDkink * kink = fEvent->GetKink(index);
5805 kink->SetTPCDensity(-1,1,0);
5806 kink->SetTPCDensity(-1,1,1);
5808 Int_t row0 = kink->GetTPCRow0() -2 - Int_t( 0.5/ (0.05+kink->GetAngle(2)));
5809 if (row0<15) row0=15;
5811 Int_t row1 = kink->GetTPCRow0() +2 + Int_t( 0.5/ (0.05+kink->GetAngle(2)));
5812 if (row1>145) row1=145;
5814 Int_t found,foundable,shared;
5815 seed->GetClusterStatistic(0,row0, found, foundable,shared,kFALSE);
5816 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),1,0);
5817 seed->GetClusterStatistic(row1,155, found, foundable,shared,kFALSE);
5818 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),1,1);
5824 Int_t AliTPCtrackerMI::CheckKinkPoint(AliTPCseed*seed,AliTPCseed &mother, AliTPCseed &daughter, AliESDkink &knk)
5827 // check kink point for given track
5828 // if return value=0 kink point not found
5829 // otherwise seed0 correspond to mother particle
5830 // seed1 correspond to daughter particle
5831 // kink parameter of kink point
5832 AliKink &kink=(AliKink &)knk;
5834 Int_t middlerow = (seed->GetFirstPoint()+seed->GetLastPoint())/2;
5835 Int_t first = seed->GetFirstPoint();
5836 Int_t last = seed->GetLastPoint();
5837 if (last-first<20) return 0; // shortest length - 2*30 = 60 pad-rows
5840 AliTPCseed *seed1 = ReSeed(seed,middlerow+20, kTRUE); //middle of chamber
5841 if (!seed1) return 0;
5842 FollowProlongation(*seed1,seed->GetLastPoint()-20);
5843 seed1->Reset(kTRUE);
5844 FollowProlongation(*seed1,158);
5845 seed1->Reset(kTRUE);
5846 last = seed1->GetLastPoint();
5848 AliTPCseed *seed0 = new AliTPCseed(*seed);
5849 seed0->Reset(kFALSE);
5852 AliTPCseed param0[20]; // parameters along the track
5853 AliTPCseed param1[20]; // parameters along the track
5854 AliKink kinks[20]; // corresponding kink parameters
5856 for (Int_t irow=0; irow<20;irow++){
5857 rows[irow] = first +((last-first)*irow)/19;
5859 // store parameters along the track
5861 for (Int_t irow=0;irow<20;irow++){
5862 FollowBackProlongation(*seed0, rows[irow]);
5863 FollowProlongation(*seed1,rows[19-irow]);
5864 param0[irow] = *seed0;
5865 param1[19-irow] = *seed1;
5869 for (Int_t irow=0; irow<19;irow++){
5870 kinks[irow].SetMother(param0[irow]);
5871 kinks[irow].SetDaughter(param1[irow]);
5872 kinks[irow].Update();
5875 // choose kink with biggest change of angle
5877 Double_t maxchange= 0;
5878 for (Int_t irow=1;irow<19;irow++){
5879 if (TMath::Abs(kinks[irow].GetR())>240.) continue;
5880 if (TMath::Abs(kinks[irow].GetR())<110.) continue;
5881 Float_t quality = TMath::Abs(kinks[irow].GetAngle(2))/(3.+TMath::Abs(kinks[irow].GetR()-param0[irow].GetX()));
5882 if ( quality > maxchange){
5883 maxchange = quality;
5890 if (index<0) return 0;
5892 Int_t row0 = GetRowNumber(kinks[index].GetR()); //row 0 estimate
5893 seed0 = new AliTPCseed(param0[index]);
5894 seed1 = new AliTPCseed(param1[index]);
5895 seed0->Reset(kFALSE);
5896 seed1->Reset(kFALSE);
5897 seed0->ResetCovariance(10.);
5898 seed1->ResetCovariance(10.);
5899 FollowProlongation(*seed0,0);
5900 FollowBackProlongation(*seed1,158);
5901 mother = *seed0; // backup mother at position 0
5902 seed0->Reset(kFALSE);
5903 seed1->Reset(kFALSE);
5904 seed0->ResetCovariance(10.);
5905 seed1->ResetCovariance(10.);
5907 first = TMath::Max(row0-20,0);
5908 last = TMath::Min(row0+20,158);
5910 for (Int_t irow=0; irow<20;irow++){
5911 rows[irow] = first +((last-first)*irow)/19;
5913 // store parameters along the track
5915 for (Int_t irow=0;irow<20;irow++){
5916 FollowBackProlongation(*seed0, rows[irow]);
5917 FollowProlongation(*seed1,rows[19-irow]);
5918 param0[irow] = *seed0;
5919 param1[19-irow] = *seed1;
5923 for (Int_t irow=0; irow<19;irow++){
5924 kinks[irow].SetMother(param0[irow]);
5925 kinks[irow].SetDaughter(param1[irow]);
5926 // param0[irow].Dump();
5927 //param1[irow].Dump();
5928 kinks[irow].Update();
5931 // choose kink with biggest change of angle
5934 for (Int_t irow=0;irow<20;irow++){
5935 if (TMath::Abs(kinks[irow].GetR())>250.) continue;
5936 if (TMath::Abs(kinks[irow].GetR())<90.) continue;
5937 Float_t quality = TMath::Abs(kinks[irow].GetAngle(2))/(3.+TMath::Abs(kinks[irow].GetR()-param0[irow].GetX()));
5938 if ( quality > maxchange){
5939 maxchange = quality;
5946 if (index==-1 || param0[index].GetNumberOfClusters()+param1[index].GetNumberOfClusters()<100){
5951 // Float_t anglesigma = TMath::Sqrt(param0[index].fC22+param0[index].fC33+param1[index].fC22+param1[index].fC33);
5953 kink.SetMother(param0[index]);
5954 kink.SetDaughter(param1[index]);
5956 row0 = GetRowNumber(kink.GetR());
5957 kink.SetTPCRow0(row0);
5958 kink.SetLabel(CookLabel(seed0,0.5,0,row0),0);
5959 kink.SetLabel(CookLabel(seed1,0.5,row0,158),1);
5960 kink.SetIndex(-10,0);
5961 kink.SetIndex(int(param0[index].GetNumberOfClusters()+param1[index].GetNumberOfClusters()),1);
5962 kink.SetTPCncls(param0[index].GetNumberOfClusters(),0);
5963 kink.SetTPCncls(param1[index].GetNumberOfClusters(),1);
5966 // new (&mother) AliTPCseed(param0[index]);
5967 daughter = param1[index];
5968 daughter.SetLabel(kink.GetLabel(1));
5969 param0[index].Reset(kTRUE);
5970 FollowProlongation(param0[index],0);
5971 mother = param0[index];
5972 mother.SetLabel(kink.GetLabel(0));
5982 AliTPCseed* AliTPCtrackerMI::ReSeed(AliTPCseed *t)
5985 // reseed - refit - track
5988 // Int_t last = fSectors->GetNRows()-1;
5990 if (fSectors == fOuterSec){
5991 first = TMath::Max(first, t->GetFirstPoint()-fInnerSec->GetNRows());
5995 first = t->GetFirstPoint();
5997 AliTPCseed * seed = MakeSeed(t,0.1,0.5,0.9);
5998 FollowBackProlongation(*t,fSectors->GetNRows()-1);
6000 FollowProlongation(*t,first);
6010 //_____________________________________________________________________________
6011 Int_t AliTPCtrackerMI::ReadSeeds(const TFile *inp) {
6012 //-----------------------------------------------------------------
6013 // This function reades track seeds.
6014 //-----------------------------------------------------------------
6015 TDirectory *savedir=gDirectory;
6017 TFile *in=(TFile*)inp;
6018 if (!in->IsOpen()) {
6019 cerr<<"AliTPCtrackerMI::ReadSeeds(): input file is not open !\n";
6024 TTree *seedTree=(TTree*)in->Get("Seeds");
6026 cerr<<"AliTPCtrackerMI::ReadSeeds(): ";
6027 cerr<<"can't get a tree with track seeds !\n";
6030 AliTPCtrack *seed=new AliTPCtrack;
6031 seedTree->SetBranchAddress("tracks",&seed);
6033 if (fSeeds==0) fSeeds=new TObjArray(15000);
6035 Int_t n=(Int_t)seedTree->GetEntries();
6036 for (Int_t i=0; i<n; i++) {
6037 seedTree->GetEvent(i);
6038 fSeeds->AddLast(new AliTPCseed(*seed/*,seed->GetAlpha()*/));
6047 Int_t AliTPCtrackerMI::Clusters2Tracks (AliESDEvent *esd)
6050 if (fSeeds) DeleteSeeds();
6053 if (!fSeeds) return 1;
6060 //_____________________________________________________________________________
6061 Int_t AliTPCtrackerMI::Clusters2Tracks() {
6062 //-----------------------------------------------------------------
6063 // This is a track finder.
6064 //-----------------------------------------------------------------
6065 TDirectory *savedir=gDirectory;
6069 fSeeds = Tracking();
6072 Info("Clusters2Tracks","Time for tracking: \t");timer.Print();timer.Start();
6074 //activate again some tracks
6075 for (Int_t i=0; i<fSeeds->GetEntriesFast(); i++) {
6076 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
6078 Int_t nc=t.GetNumberOfClusters();
6080 delete fSeeds->RemoveAt(i);
6084 if (pt->GetRemoval()==10) {
6085 if (pt->GetDensityFirst(20)>0.8 || pt->GetDensityFirst(30)>0.8 || pt->GetDensityFirst(40)>0.7)
6086 pt->Desactivate(10); // make track again active
6088 pt->Desactivate(20);
6089 delete fSeeds->RemoveAt(i);
6094 RemoveUsed2(fSeeds,0.85,0.85,0);
6095 if (AliTPCReconstructor::GetRecoParam()->GetDoKinks()) FindKinks(fSeeds,fEvent);
6096 //FindCurling(fSeeds, fEvent,0);
6097 if (AliTPCReconstructor::StreamLevel()>2) FindMultiMC(fSeeds, fEvent,-1); // find multi found tracks
6098 RemoveUsed2(fSeeds,0.5,0.4,20);
6099 FindSplitted(fSeeds, fEvent,0); // find multi found tracks
6100 if (AliTPCReconstructor::StreamLevel()>2) FindMultiMC(fSeeds, fEvent,0); // find multi found tracks
6103 // // refit short tracks
6105 Int_t nseed=fSeeds->GetEntriesFast();
6108 for (Int_t i=0; i<nseed; i++) {
6109 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
6111 Int_t nc=t.GetNumberOfClusters();
6113 delete fSeeds->RemoveAt(i);
6116 CookLabel(pt,0.1); //For comparison only
6117 //if ((pt->IsActive() || (pt->fRemoval==10) )&& nc>50 &&pt->GetNumberOfClusters()>0.4*pt->fNFoundable){
6118 if ((pt->IsActive() || (pt->GetRemoval()==10) )){
6120 if (fDebug>0) cerr<<found<<'\r';
6124 delete fSeeds->RemoveAt(i);
6128 //RemoveOverlap(fSeeds,0.99,7,kTRUE);
6130 //RemoveUsed(fSeeds,0.9,0.9,6);
6132 nseed=fSeeds->GetEntriesFast();
6134 for (Int_t i=0; i<nseed; i++) {
6135 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
6137 Int_t nc=t.GetNumberOfClusters();
6139 delete fSeeds->RemoveAt(i);
6143 t.CookdEdx(0.02,0.6);
6144 // CheckKinkPoint(&t,0.05);
6145 //if ((pt->IsActive() || (pt->fRemoval==10) )&& nc>50 &&pt->GetNumberOfClusters()>0.4*pt->fNFoundable){
6146 if ((pt->IsActive() || (pt->GetRemoval()==10) )){
6154 delete fSeeds->RemoveAt(i);
6155 //AliTPCseed * seed1 = ReSeed(pt,0.05,0.5,1);
6157 // FollowProlongation(*seed1,0);
6158 // Int_t n = seed1->GetNumberOfClusters();
6159 // printf("fP4\t%f\t%f\n",seed1->GetC(),pt->GetC());
6160 // printf("fN\t%d\t%d\n", seed1->GetNumberOfClusters(),pt->GetNumberOfClusters());
6163 //AliTPCseed * seed2 = ReSeed(pt,0.95,0.5,0.05);
6167 SortTracks(fSeeds, 1);
6171 PrepareForBackProlongation(fSeeds,5.);
6172 PropagateBack(fSeeds);
6173 printf("Time for back propagation: \t");timer.Print();timer.Start();
6177 PrepareForProlongation(fSeeds,5.);
6178 PropagateForward2(fSeeds);
6180 printf("Time for FORWARD propagation: \t");timer.Print();timer.Start();
6181 // RemoveUsed(fSeeds,0.7,0.7,6);
6182 //RemoveOverlap(fSeeds,0.9,7,kTRUE);
6184 nseed=fSeeds->GetEntriesFast();
6186 for (Int_t i=0; i<nseed; i++) {
6187 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
6189 Int_t nc=t.GetNumberOfClusters();
6191 delete fSeeds->RemoveAt(i);
6194 t.CookdEdx(0.02,0.6);
6195 // CookLabel(pt,0.1); //For comparison only
6196 //if ((pt->IsActive() || (pt->fRemoval==10) )&& nc>50 &&pt->GetNumberOfClusters()>0.4*pt->fNFoundable){
6197 if ((pt->IsActive() || (pt->fRemoval==10) )){
6198 cerr<<found++<<'\r';
6201 delete fSeeds->RemoveAt(i);
6206 // fNTracks = found;
6208 Info("Clusters2Tracks","Time for overlap removal, track writing and dedx cooking: \t"); timer.Print();timer.Start();
6211 // cerr<<"Number of found tracks : "<<"\t"<<found<<endl;
6212 Info("Clusters2Tracks","Number of found tracks %d",found);
6214 // UnloadClusters();
6219 void AliTPCtrackerMI::Tracking(TObjArray * arr)
6222 // tracking of the seeds
6225 fSectors = fOuterSec;
6226 ParallelTracking(arr,150,63);
6227 fSectors = fOuterSec;
6228 ParallelTracking(arr,63,0);
6231 TObjArray * AliTPCtrackerMI::Tracking(Int_t seedtype, Int_t i1, Int_t i2, Float_t cuts[4], Float_t dy, Int_t dsec)
6236 TObjArray * arr = new TObjArray;
6238 fSectors = fOuterSec;
6241 for (Int_t sec=0;sec<fkNOS;sec++){
6242 if (seedtype==3) MakeSeeds3(arr,sec,i1,i2,cuts,dy, dsec);
6243 if (seedtype==4) MakeSeeds5(arr,sec,i1,i2,cuts,dy);
6244 if (seedtype==2) MakeSeeds2(arr,sec,i1,i2,cuts,dy);
6247 Info("Tracking","\nSeeding - %d\t%d\t%d\t%d\n",seedtype,i1,i2,arr->GetEntriesFast());
6259 TObjArray * AliTPCtrackerMI::Tracking()
6263 if (AliTPCReconstructor::GetRecoParam()->GetSpecialSeeding()) return TrackingSpecial();
6266 Int_t nup=fOuterSec->GetNRows()+fInnerSec->GetNRows();
6268 TObjArray * seeds = new TObjArray;
6277 Float_t fnumber = 3.0;
6278 Float_t fdensity = 3.0;
6283 for (Int_t delta = 0; delta<18; delta+=6){
6287 arr = Tracking(3,nup-1-delta,nup-1-delta-gap,cuts,-1,1);
6288 SumTracks(seeds,arr);
6289 SignClusters(seeds,fnumber,fdensity);
6291 for (Int_t i=2;i<6;i+=2){
6292 // seed high pt tracks
6295 arr = Tracking(3,nup-i-delta,nup-i-delta-gap,cuts,-1,0);
6296 SumTracks(seeds,arr);
6297 SignClusters(seeds,fnumber,fdensity);
6302 // RemoveUsed(seeds,0.9,0.9,1);
6303 // UnsignClusters();
6304 // SignClusters(seeds,fnumber,fdensity);
6308 for (Int_t delta = 20; delta<120; delta+=10){
6310 // seed high pt tracks
6314 arr = Tracking(3,nup-delta,nup-delta-gap,cuts,-1);
6315 SumTracks(seeds,arr);
6316 SignClusters(seeds,fnumber,fdensity);
6321 arr = Tracking(3,nup-delta-5,nup-delta-5-gap,cuts,-1);
6322 SumTracks(seeds,arr);
6323 SignClusters(seeds,fnumber,fdensity);
6334 Info("Tracking()","\n\nPrimary seeding\t%d\n\n",seeds->GetEntriesFast());
6338 // RemoveUsed(seeds,0.75,0.75,1);
6340 //SignClusters(seeds,fnumber,fdensity);
6349 arr = Tracking(4,nup-1,nup-1-gap,cuts,-1);
6350 SumTracks(seeds,arr);
6351 SignClusters(seeds,fnumber,fdensity);
6353 arr = Tracking(4,nup-2,nup-2-gap,cuts,-1);
6354 SumTracks(seeds,arr);
6355 SignClusters(seeds,fnumber,fdensity);
6357 arr = Tracking(4,nup-3,nup-3-gap,cuts,-1);
6358 SumTracks(seeds,arr);
6359 SignClusters(seeds,fnumber,fdensity);
6363 for (Int_t delta = 3; delta<30; delta+=5){
6369 arr = Tracking(4,nup-1-delta,nup-1-delta-gap,cuts,-1);
6370 SumTracks(seeds,arr);
6371 SignClusters(seeds,fnumber,fdensity);
6373 arr = Tracking(4,nup-3-delta,nup-5-delta-gap,cuts,4);
6374 SumTracks(seeds,arr);
6375 SignClusters(seeds,fnumber,fdensity);
6386 Int_t fLastSeedRowSec=AliTPCReconstructor::GetRecoParam()->GetLastSeedRowSec();
6389 for (Int_t delta = 30; delta<fLastSeedRowSec; delta+=10){
6395 arr = Tracking(4,nup-1-delta,nup-1-delta-gap,cuts,-1);
6396 SumTracks(seeds,arr);
6397 SignClusters(seeds,fnumber,fdensity);
6399 arr = Tracking(4,nup-5-delta,nup-5-delta-gap,cuts,5 );
6400 SumTracks(seeds,arr);
6401 SignClusters(seeds,fnumber,fdensity);
6405 Info("Tracking()","\n\nSecondary seeding\t%d\n\n",seeds->GetEntriesFast());
6416 TObjArray * AliTPCtrackerMI::TrackingSpecial()
6419 // seeding adjusted for laser and cosmic tests - short tracks with big inclination angle
6420 // no primary vertex seeding tried
6424 Int_t nup=fOuterSec->GetNRows()+fInnerSec->GetNRows();
6426 TObjArray * seeds = new TObjArray;
6431 Float_t fnumber = 3.0;
6432 Float_t fdensity = 3.0;
6435 cuts[0] = AliTPCReconstructor::GetRecoParam()->GetMaxC(); // max curvature
6436 cuts[1] = 3.5; // max tan(phi) angle for seeding
6437 cuts[2] = 3.; // not used (cut on z primary vertex)
6438 cuts[3] = 3.5; // max tan(theta) angle for seeding
6440 for (Int_t delta = 0; nup-delta-gap-1>0; delta+=3){
6442 arr = Tracking(4,nup-1-delta,nup-1-delta-gap,cuts,-1);
6443 SumTracks(seeds,arr);
6444 SignClusters(seeds,fnumber,fdensity);
6448 Info("Tracking()","\n\nSecondary seeding\t%d\n\n",seeds->GetEntriesFast());
6459 void AliTPCtrackerMI::SumTracks(TObjArray *arr1,TObjArray *&arr2) const
6462 //sum tracks to common container
6463 //remove suspicious tracks
6464 Int_t nseed = arr2->GetEntriesFast();
6465 for (Int_t i=0;i<nseed;i++){
6466 AliTPCseed *pt=(AliTPCseed*)arr2->UncheckedAt(i);
6469 // remove tracks with too big curvature
6471 if (TMath::Abs(pt->GetC())>AliTPCReconstructor::GetRecoParam()->GetMaxC()){
6472 delete arr2->RemoveAt(i);
6475 // REMOVE VERY SHORT TRACKS
6476 if (pt->GetNumberOfClusters()<20){
6477 delete arr2->RemoveAt(i);
6480 // NORMAL ACTIVE TRACK
6481 if (pt->IsActive()){
6482 arr1->AddLast(arr2->RemoveAt(i));
6485 //remove not usable tracks
6486 if (pt->GetRemoval()!=10){
6487 delete arr2->RemoveAt(i);
6491 // ENABLE ONLY ENOUGH GOOD STOPPED TRACKS
6492 if (pt->GetDensityFirst(20)>0.8 || pt->GetDensityFirst(30)>0.8 || pt->GetDensityFirst(40)>0.7)
6493 arr1->AddLast(arr2->RemoveAt(i));
6495 delete arr2->RemoveAt(i);
6499 delete arr2; arr2 = 0;
6504 void AliTPCtrackerMI::ParallelTracking(TObjArray * arr, Int_t rfirst, Int_t rlast)
6507 // try to track in parralel
6509 Int_t nseed=arr->GetEntriesFast();
6510 //prepare seeds for tracking
6511 for (Int_t i=0; i<nseed; i++) {
6512 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i), &t=*pt;
6514 if (!t.IsActive()) continue;
6515 // follow prolongation to the first layer
6516 if ( (fSectors ==fInnerSec) || (t.GetFirstPoint()-fParam->GetNRowLow()>rfirst+1) )
6517 FollowProlongation(t, rfirst+1);
6522 for (Int_t nr=rfirst; nr>=rlast; nr--){
6523 if (nr<fInnerSec->GetNRows())
6524 fSectors = fInnerSec;
6526 fSectors = fOuterSec;
6527 // make indexes with the cluster tracks for given
6529 // find nearest cluster
6530 for (Int_t i=0; i<nseed; i++) {
6531 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i), &t=*pt;
6533 if (nr==80) pt->UpdateReference();
6534 if (!pt->IsActive()) continue;
6535 // if ( (fSectors ==fOuterSec) && (pt->fFirstPoint-fParam->GetNRowLow())<nr) continue;
6536 if (pt->GetRelativeSector()>17) {
6539 UpdateClusters(t,nr);
6541 // prolonagate to the nearest cluster - if founded
6542 for (Int_t i=0; i<nseed; i++) {
6543 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
6545 if (!pt->IsActive()) continue;
6546 // if ((fSectors ==fOuterSec) && (pt->fFirstPoint-fParam->GetNRowLow())<nr) continue;
6547 if (pt->GetRelativeSector()>17) {
6550 FollowToNextCluster(*pt,nr);
6555 void AliTPCtrackerMI::PrepareForBackProlongation(TObjArray * arr,Float_t fac) const
6559 // if we use TPC track itself we have to "update" covariance
6561 Int_t nseed= arr->GetEntriesFast();
6562 for (Int_t i=0;i<nseed;i++){
6563 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
6567 //rotate to current local system at first accepted point
6568 Int_t index = pt->GetClusterIndex2(pt->GetFirstPoint());
6569 Int_t sec = (index&0xff000000)>>24;
6571 Float_t angle1 = fInnerSec->GetAlpha()*sec+fInnerSec->GetAlphaShift();
6572 if (angle1>TMath::Pi())
6573 angle1-=2.*TMath::Pi();
6574 Float_t angle2 = pt->GetAlpha();
6576 if (TMath::Abs(angle1-angle2)>0.001){
6577 pt->Rotate(angle1-angle2);
6578 //angle2 = pt->GetAlpha();
6579 //pt->fRelativeSector = pt->GetAlpha()/fInnerSec->GetAlpha();
6580 //if (pt->GetAlpha()<0)
6581 // pt->fRelativeSector+=18;
6582 //sec = pt->fRelativeSector;
6591 void AliTPCtrackerMI::PrepareForProlongation(TObjArray * arr, Float_t fac) const
6595 // if we use TPC track itself we have to "update" covariance
6597 Int_t nseed= arr->GetEntriesFast();
6598 for (Int_t i=0;i<nseed;i++){
6599 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
6602 pt->SetFirstPoint(pt->GetLastPoint());
6610 Int_t AliTPCtrackerMI::PropagateBack(TObjArray * arr)
6613 // make back propagation
6615 Int_t nseed= arr->GetEntriesFast();
6616 for (Int_t i=0;i<nseed;i++){
6617 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
6618 if (pt&& pt->GetKinkIndex(0)<=0) {
6619 //AliTPCseed *pt2 = new AliTPCseed(*pt);
6620 fSectors = fInnerSec;
6621 //FollowBackProlongation(*pt,fInnerSec->GetNRows()-1);
6622 //fSectors = fOuterSec;
6623 FollowBackProlongation(*pt,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1);
6624 //if (pt->GetNumberOfClusters()<(pt->fEsd->GetTPCclusters(0)) ){
6625 // Error("PropagateBack","Not prolonged track %d",pt->GetLabel());
6626 // FollowBackProlongation(*pt2,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1);
6629 if (pt&& pt->GetKinkIndex(0)>0) {
6630 AliESDkink * kink = fEvent->GetKink(pt->GetKinkIndex(0)-1);
6631 pt->SetFirstPoint(kink->GetTPCRow0());
6632 fSectors = fInnerSec;
6633 FollowBackProlongation(*pt,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1);
6641 Int_t AliTPCtrackerMI::PropagateForward2(TObjArray * arr)
6644 // make forward propagation
6646 Int_t nseed= arr->GetEntriesFast();
6648 for (Int_t i=0;i<nseed;i++){
6649 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
6651 FollowProlongation(*pt,0);
6660 Int_t AliTPCtrackerMI::PropagateForward()
6663 // propagate track forward
6665 Int_t nseed = fSeeds->GetEntriesFast();
6666 for (Int_t i=0;i<nseed;i++){
6667 AliTPCseed *pt = (AliTPCseed*)fSeeds->UncheckedAt(i);
6669 AliTPCseed &t = *pt;
6670 Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
6671 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
6672 if (alpha < 0. ) alpha += 2.*TMath::Pi();
6673 t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN);
6677 fSectors = fOuterSec;
6678 ParallelTracking(fSeeds,fOuterSec->GetNRows()+fInnerSec->GetNRows()-1,fInnerSec->GetNRows());
6679 fSectors = fInnerSec;
6680 ParallelTracking(fSeeds,fInnerSec->GetNRows()-1,0);
6689 Int_t AliTPCtrackerMI::PropagateBack(AliTPCseed * pt, Int_t row0, Int_t row1)
6692 // make back propagation, in between row0 and row1
6696 fSectors = fInnerSec;
6699 if (row1<fSectors->GetNRows())
6702 r1 = fSectors->GetNRows()-1;
6704 if (row0<fSectors->GetNRows()&& r1>0 )
6705 FollowBackProlongation(*pt,r1);
6706 if (row1<=fSectors->GetNRows())
6709 r1 = row1 - fSectors->GetNRows();
6710 if (r1<=0) return 0;
6711 if (r1>=fOuterSec->GetNRows()) return 0;
6712 fSectors = fOuterSec;
6713 return FollowBackProlongation(*pt,r1);
6721 void AliTPCtrackerMI::GetShape(AliTPCseed * seed, Int_t row)
6725 AliTPCClusterParam * clparam = AliTPCcalibDB::Instance()->GetClusterParam();
6726 Float_t zdrift = TMath::Abs((fParam->GetZLength(0)-TMath::Abs(seed->GetZ())));
6727 Int_t type = (seed->GetSector() < fParam->GetNSector()/2) ? 0: (row>126) ? 1:2;
6728 Double_t angulary = seed->GetSnp();
6729 angulary = angulary*angulary/(1.-angulary*angulary);
6730 Double_t angularz = seed->GetTgl()*seed->GetTgl()*(1.+angulary);
6732 Double_t sigmay = clparam->GetRMS0(0,type,zdrift,TMath::Sqrt(TMath::Abs(angulary)));
6733 Double_t sigmaz = clparam->GetRMS0(1,type,zdrift,TMath::Sqrt(TMath::Abs(angularz)));
6734 seed->SetCurrentSigmaY2(sigmay*sigmay);
6735 seed->SetCurrentSigmaZ2(sigmaz*sigmaz);
6736 // Float_t sd2 = TMath::Abs((fParam->GetZLength(0)-TMath::Abs(seed->GetZ())))*fParam->GetDiffL()*fParam->GetDiffL();
6737 // // Float_t padlength = fParam->GetPadPitchLength(seed->fSector);
6738 // Float_t padlength = GetPadPitchLength(row);
6740 // Float_t sresy = (seed->GetSector() < fParam->GetNSector()/2) ? 0.2 :0.3;
6741 // seed->SetCurrentSigmaY2(sd2+padlength*padlength*angulary/12.+sresy*sresy);
6743 // Float_t sresz = fParam->GetZSigma();
6744 // seed->SetCurrentSigmaZ2(sd2+padlength*padlength*angularz*angularz*(1+angulary)/12.+sresz*sresz);
6746 Float_t wy = GetSigmaY(seed);
6747 Float_t wz = GetSigmaZ(seed);
6750 if (TMath::Abs(wy/seed->fCurrentSigmaY2-1)>0.0001 || TMath::Abs(wz/seed->fCurrentSigmaZ2-1)>0.0001 ){
6751 printf("problem\n");
6758 //__________________________________________________________________________
6759 void AliTPCtrackerMI::CookLabel(AliKalmanTrack *tk, Float_t wrong) const {
6760 //--------------------------------------------------------------------
6761 //This function "cooks" a track label. If label<0, this track is fake.
6762 //--------------------------------------------------------------------
6763 AliTPCseed * t = dynamic_cast<AliTPCseed*>(tk);
6765 printf("%s:%d wrong type \n",(char*)__FILE__,__LINE__);
6769 Int_t noc=t->GetNumberOfClusters();
6771 //printf("\nnot founded prolongation\n\n\n");
6777 AliTPCclusterMI *clusters[160];
6779 for (Int_t i=0;i<160;i++) {
6786 for (i=0; i<160 && current<noc; i++) {
6788 Int_t index=t->GetClusterIndex2(i);
6789 if (index<=0) continue;
6790 if (index&0x8000) continue;
6792 //clusters[current]=GetClusterMI(index);
6793 if (t->GetClusterPointer(i)){
6794 clusters[current]=t->GetClusterPointer(i);
6800 Int_t lab=123456789;
6801 for (i=0; i<noc; i++) {
6802 AliTPCclusterMI *c=clusters[i];
6804 lab=TMath::Abs(c->GetLabel(0));
6806 for (j=0; j<noc; j++) if (lb[j]==lab || mx[j]==0) break;
6812 for (i=0; i<noc; i++) if (mx[i]>max) {max=mx[i]; lab=lb[i];}
6814 for (i=0; i<noc; i++) {
6815 AliTPCclusterMI *c=clusters[i];
6817 if (TMath::Abs(c->GetLabel(1)) == lab ||
6818 TMath::Abs(c->GetLabel(2)) == lab ) max++;
6821 if ((1.- Float_t(max)/noc) > wrong) lab=-lab;
6824 Int_t tail=Int_t(0.10*noc);
6827 for (i=1; i<=160&&ind<tail; i++) {
6828 // AliTPCclusterMI *c=clusters[noc-i];
6829 AliTPCclusterMI *c=clusters[i];
6831 if (lab == TMath::Abs(c->GetLabel(0)) ||
6832 lab == TMath::Abs(c->GetLabel(1)) ||
6833 lab == TMath::Abs(c->GetLabel(2))) max++;
6836 if (max < Int_t(0.5*tail)) lab=-lab;
6843 //delete[] clusters;
6847 //__________________________________________________________________________
6848 Int_t AliTPCtrackerMI::CookLabel(AliTPCseed *t, Float_t wrong,Int_t first, Int_t last) const {
6849 //--------------------------------------------------------------------
6850 //This function "cooks" a track label. If label<0, this track is fake.
6851 //--------------------------------------------------------------------
6852 Int_t noc=t->GetNumberOfClusters();
6854 //printf("\nnot founded prolongation\n\n\n");
6860 AliTPCclusterMI *clusters[160];
6862 for (Int_t i=0;i<160;i++) {
6869 for (i=0; i<160 && current<noc; i++) {
6870 if (i<first) continue;
6871 if (i>last) continue;
6872 Int_t index=t->GetClusterIndex2(i);
6873 if (index<=0) continue;
6874 if (index&0x8000) continue;
6876 //clusters[current]=GetClusterMI(index);
6877 if (t->GetClusterPointer(i)){
6878 clusters[current]=t->GetClusterPointer(i);
6883 if (noc<5) return -1;
6884 Int_t lab=123456789;
6885 for (i=0; i<noc; i++) {
6886 AliTPCclusterMI *c=clusters[i];
6888 lab=TMath::Abs(c->GetLabel(0));
6890 for (j=0; j<noc; j++) if (lb[j]==lab || mx[j]==0) break;
6896 for (i=0; i<noc; i++) if (mx[i]>max) {max=mx[i]; lab=lb[i];}
6898 for (i=0; i<noc; i++) {
6899 AliTPCclusterMI *c=clusters[i];
6901 if (TMath::Abs(c->GetLabel(1)) == lab ||
6902 TMath::Abs(c->GetLabel(2)) == lab ) max++;
6905 if ((1.- Float_t(max)/noc) > wrong) lab=-lab;
6908 Int_t tail=Int_t(0.10*noc);
6911 for (i=1; i<=160&&ind<tail; i++) {
6912 // AliTPCclusterMI *c=clusters[noc-i];
6913 AliTPCclusterMI *c=clusters[i];
6915 if (lab == TMath::Abs(c->GetLabel(0)) ||
6916 lab == TMath::Abs(c->GetLabel(1)) ||
6917 lab == TMath::Abs(c->GetLabel(2))) max++;
6920 if (max < Int_t(0.5*tail)) lab=-lab;
6923 // t->SetLabel(lab);
6927 //delete[] clusters;
6931 Int_t AliTPCtrackerMI::GetRowNumber(Double_t x[3]) const
6933 //return pad row number for given x vector
6934 Float_t phi = TMath::ATan2(x[1],x[0]);
6935 if(phi<0) phi=2.*TMath::Pi()+phi;
6936 // Get the local angle in the sector philoc
6937 const Float_t kRaddeg = 180/3.14159265358979312;
6938 Float_t phiangle = (Int_t (phi*kRaddeg/20.) + 0.5)*20./kRaddeg;
6939 Double_t localx = x[0]*TMath::Cos(phiangle)-x[1]*TMath::Sin(phiangle);
6940 return GetRowNumber(localx);
6945 void AliTPCtrackerMI::MakeBitmaps(AliTPCseed *t)
6947 //-----------------------------------------------------------------------
6948 // Fill the cluster and sharing bitmaps of the track
6949 //-----------------------------------------------------------------------
6951 Int_t firstpoint = 0;
6952 Int_t lastpoint = 159;
6953 AliTPCTrackerPoint *point;
6955 for (int iter=firstpoint; iter<lastpoint; iter++) {
6956 point = t->GetTrackPoint(iter);
6958 t->SetClusterMapBit(iter, kTRUE);
6959 if (point->IsShared())
6960 t->SetSharedMapBit(iter,kTRUE);
6962 t->SetSharedMapBit(iter, kFALSE);
6965 t->SetClusterMapBit(iter, kFALSE);
6966 t->SetSharedMapBit(iter, kFALSE);
6971 Bool_t AliTPCtrackerMI::IsFindable(AliTPCseed & track){
6973 // return flag if there is findable cluster at given position
6976 Float_t z = track.GetZ();
6978 if (TMath::Abs(z)<(AliTPCReconstructor::GetCtgRange()*track.GetX()+kDeltaZ) &&
6979 TMath::Abs(z)<fParam->GetZLength(0) &&
6980 (TMath::Abs(track.GetSnp())<AliTPCReconstructor::GetMaxSnpTracker()))
6986 void AliTPCtrackerMI::AddCovariance(AliTPCseed * seed){
6988 // Adding systematic error
6989 // !!!! the systematic error for element 4 is in 1/cm not in pt
6991 const Double_t *param = AliTPCReconstructor::GetRecoParam()->GetSystematicError();
6993 for (Int_t i=0;i<15;i++) covar[i]=0;
6999 covar[0] = param[0]*param[0];
7000 covar[2] = param[1]*param[1];
7001 covar[5] = param[2]*param[2];
7002 covar[9] = param[3]*param[3];
7003 Double_t facC = AliTracker::GetBz()*kB2C;
7004 covar[14]= param[4]*param[4]*facC*facC;
7005 seed->AddCovariance(covar);