1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
17 //-------------------------------------------------------
18 // Implementation of the TPC tracker
20 // Origin: Marian Ivanov Marian.Ivanov@cern.ch
22 // AliTPC parallel tracker
24 // The track fitting is based on Kalaman filtering approach
26 // The track finding steps:
27 // 1. Seeding - with and without vertex constraint
28 // - seeding with vertex constain done at first n^2 proble
29 // - seeding without vertex constraint n^3 problem
30 // 2. Tracking - follow prolongation road - find cluster - update kalman track
32 // The seeding and tracking is repeated several times, in different seeding region.
33 // This approach enables to find the track which cannot be seeded in some region of TPC
34 // This can happen because of low momenta (track do not reach outer radius), or track is currently in the ded region between sectors, or the track is for the moment overlapped with other track (seed quality is poor) ...
36 // With this approach we reach almost 100 % efficiency also for high occupancy events.
37 // (If the seeding efficiency in a region is about 90 % than with logical or of several
38 // regions we will reach 100% (in theory - supposing independence)
40 // Repeating several seeding - tracking procedures some of the tracks can be find
43 // The procedures to remove multi find tacks are impremented:
44 // RemoveUsed2 - fast procedure n problem -
45 // Algorithm - Sorting tracks according quality
46 // remove tracks with some shared fraction
47 // Sharing in respect to all tacks
48 // Signing clusters in gold region
49 // FindSplitted - slower algorithm n^2
50 // Sort the tracks according quality
51 // Loop over pair of tracks
52 // If overlap with other track bigger than threshold - remove track
54 // FindCurling - Finds the pair of tracks which are curling
55 // - About 10% of tracks can be find with this procedure
56 // The combinatorial background is too big to be used in High
57 // multiplicity environment
58 // - n^2 problem - Slow procedure - currently it is disabled because of
61 // The number of splitted tracks can be reduced disabling the sharing of the cluster.
62 // tpcRecoParam-> SetClusterSharing(kFALSE);
63 // IT IS HIGHLY non recomended to use it in high flux enviroonment
64 // Even using this switch some tracks can be found more than once
65 // (because of multiple seeding and low quality tracks which will not cross full chamber)
68 // The tracker itself can be debugged - the information about tracks can be stored in several // phases of the reconstruction
69 // To enable storage of the TPC tracks in the ESD friend track
70 // use AliTPCReconstructor::SetStreamLevel(n); where nis bigger 0
72 // The debug level - different procedure produce tree for numerical debugging
73 // To enable them set AliTPCReconstructor::SetStreamLevel(n); where nis bigger 1
77 // Adding systematic errors to the covariance:
79 // The systematic errors due to the misalignment and miscalibration are added to the covariance matrix
80 // of the tracks (not to the clusters as they are dependent):
81 // The parameters form AliTPCRecoParam are used AliTPCRecoParam::GetSystematicError
82 // The systematic errors are expressed there in RMS - position (cm), angle (rad), curvature (1/cm)
83 // The default values are 0.
85 // The sytematic errors are added to the covariance matrix in following places:
87 // 1. During fisrt itteration - AliTPCtrackerMI::FillESD
88 // 2. Second iteration -
89 // 2.a ITS->TPC - AliTPCtrackerMI::ReadSeeds
90 // 2.b TPC->TRD - AliTPCtrackerMI::PropagateBack
91 // 3. Third iteration -
92 // 3.a TRD->TPC - AliTPCtrackerMI::ReadSeeds
93 // 3.b TPC->ITS - AliTPCtrackerMI::RefitInward
95 // There are several places in the code which can be numerically debuged
96 // This code is keeped in order to enable code development and to check the calibration implementtion
98 // 1. ErrParam stream (Log level 9) - dump information about
100 // 2.a) cluster error estimate
101 // 3.a) cluster shape estimate
104 //-------------------------------------------------------
109 #include "Riostream.h"
110 #include <TClonesArray.h>
112 #include <TObjArray.h>
115 #include "AliComplexCluster.h"
116 #include "AliESDEvent.h"
117 #include "AliESDtrack.h"
118 #include "AliESDVertex.h"
121 #include "AliHelix.h"
122 #include "AliRunLoader.h"
123 #include "AliTPCClustersRow.h"
124 #include "AliTPCParam.h"
125 #include "AliTPCReconstructor.h"
126 #include "AliTPCpolyTrack.h"
127 #include "AliTPCreco.h"
128 #include "AliTPCseed.h"
130 #include "AliTPCtrackerSector.h"
131 #include "AliTPCtrackerMI.h"
132 #include "TStopwatch.h"
133 #include "AliTPCReconstructor.h"
135 #include "TTreeStream.h"
136 #include "AliAlignObj.h"
137 #include "AliTrackPointArray.h"
139 #include "AliTPCcalibDB.h"
140 #include "AliTPCTransform.h"
141 #include "AliTPCClusterParam.h"
145 ClassImp(AliTPCtrackerMI)
149 class AliTPCFastMath {
152 static Double_t FastAsin(Double_t x);
154 static Double_t fgFastAsin[20000]; //lookup table for fast asin computation
157 Double_t AliTPCFastMath::fgFastAsin[20000];
158 AliTPCFastMath gAliTPCFastMath; // needed to fill the LUT
160 AliTPCFastMath::AliTPCFastMath(){
162 // initialized lookup table;
163 for (Int_t i=0;i<10000;i++){
164 fgFastAsin[2*i] = TMath::ASin(i/10000.);
165 fgFastAsin[2*i+1] = (TMath::ASin((i+1)/10000.)-fgFastAsin[2*i]);
169 Double_t AliTPCFastMath::FastAsin(Double_t x){
171 // return asin using lookup table
173 Int_t index = int(x*10000);
174 return fgFastAsin[2*index]+(x*10000.-index)*fgFastAsin[2*index+1];
177 Int_t index = int(x*10000);
178 return -(fgFastAsin[2*index]+(x*10000.-index)*fgFastAsin[2*index+1]);
180 //__________________________________________________________________
181 AliTPCtrackerMI::AliTPCtrackerMI()
203 // default constructor
206 //_____________________________________________________________________
210 Int_t AliTPCtrackerMI::UpdateTrack(AliTPCseed * track, Int_t accept){
212 //update track information using current cluster - track->fCurrentCluster
215 AliTPCclusterMI* c =track->GetCurrentCluster();
216 if (accept>0) track->SetCurrentClusterIndex1(track->GetCurrentClusterIndex1() | 0x8000); //sign not accepted clusters
218 UInt_t i = track->GetCurrentClusterIndex1();
220 Int_t sec=(i&0xff000000)>>24;
221 //Int_t row = (i&0x00ff0000)>>16;
222 track->SetRow((i&0x00ff0000)>>16);
223 track->SetSector(sec);
224 // Int_t index = i&0xFFFF;
225 if (sec>=fParam->GetNInnerSector()) track->SetRow(track->GetRow()+fParam->GetNRowLow());
226 track->SetClusterIndex2(track->GetRow(), i);
227 //track->fFirstPoint = row;
228 //if ( track->fLastPoint<row) track->fLastPoint =row;
229 // if (track->fRow<0 || track->fRow>160) {
230 // printf("problem\n");
232 if (track->GetFirstPoint()>track->GetRow())
233 track->SetFirstPoint(track->GetRow());
234 if (track->GetLastPoint()<track->GetRow())
235 track->SetLastPoint(track->GetRow());
238 track->SetClusterPointer(track->GetRow(),c);
241 Double_t angle2 = track->GetSnp()*track->GetSnp();
243 //SET NEW Track Point
245 if (angle2<1) //PH sometimes angle2 is very big. To be investigated...
247 angle2 = TMath::Sqrt(angle2/(1-angle2));
248 AliTPCTrackerPoint &point =*(track->GetTrackPoint(track->GetRow()));
250 point.SetSigmaY(c->GetSigmaY2()/track->GetCurrentSigmaY2());
251 point.SetSigmaZ(c->GetSigmaZ2()/track->GetCurrentSigmaZ2());
252 point.SetErrY(sqrt(track->GetErrorY2()));
253 point.SetErrZ(sqrt(track->GetErrorZ2()));
255 point.SetX(track->GetX());
256 point.SetY(track->GetY());
257 point.SetZ(track->GetZ());
258 point.SetAngleY(angle2);
259 point.SetAngleZ(track->GetTgl());
260 if (point.IsShared()){
261 track->SetErrorY2(track->GetErrorY2()*4);
262 track->SetErrorZ2(track->GetErrorZ2()*4);
266 Double_t chi2 = track->GetPredictedChi2(track->GetCurrentCluster());
268 // track->SetErrorY2(track->GetErrorY2()*1.3);
269 // track->SetErrorY2(track->GetErrorY2()+0.01);
270 // track->SetErrorZ2(track->GetErrorZ2()*1.3);
271 // track->SetErrorZ2(track->GetErrorZ2()+0.005);
273 if (accept>0) return 0;
274 if (track->GetNumberOfClusters()%20==0){
275 // if (track->fHelixIn){
276 // TClonesArray & larr = *(track->fHelixIn);
277 // Int_t ihelix = larr.GetEntriesFast();
278 // new(larr[ihelix]) AliHelix(*track) ;
281 track->SetNoCluster(0);
282 return track->Update(c,chi2,i);
287 Int_t AliTPCtrackerMI::AcceptCluster(AliTPCseed * seed, AliTPCclusterMI * cluster)
290 // decide according desired precision to accept given
291 // cluster for tracking
292 Double_t sy2=ErrY2(seed,cluster);
293 Double_t sz2=ErrZ2(seed,cluster);
295 Double_t sdistancey2 = sy2+seed->GetSigmaY2();
296 Double_t sdistancez2 = sz2+seed->GetSigmaZ2();
298 Double_t rdistancey2 = (seed->GetCurrentCluster()->GetY()-seed->GetY())*
299 (seed->GetCurrentCluster()->GetY()-seed->GetY())/sdistancey2;
300 Double_t rdistancez2 = (seed->GetCurrentCluster()->GetZ()-seed->GetZ())*
301 (seed->GetCurrentCluster()->GetZ()-seed->GetZ())/sdistancez2;
303 Double_t rdistance2 = rdistancey2+rdistancez2;
306 if (AliTPCReconstructor::StreamLevel()>5 && seed->GetNumberOfClusters()>20) {
307 Float_t rmsy2 = seed->GetCurrentSigmaY2();
308 Float_t rmsz2 = seed->GetCurrentSigmaZ2();
309 Float_t rmsy2p30 = seed->GetCMeanSigmaY2p30();
310 Float_t rmsz2p30 = seed->GetCMeanSigmaZ2p30();
311 Float_t rmsy2p30R = seed->GetCMeanSigmaY2p30R();
312 Float_t rmsz2p30R = seed->GetCMeanSigmaZ2p30R();
314 AliExternalTrackParam param(*seed);
315 static TVectorD gcl(3),gtr(3);
317 param.GetXYZ(gcl.GetMatrixArray());
318 cluster->GetGlobalXYZ(gclf);
319 gcl[0]=gclf[0]; gcl[1]=gclf[1]; gcl[2]=gclf[2];
320 (*fDebugStreamer)<<"ErrParam"<<
329 "rmsy2p30="<<rmsy2p30<<
330 "rmsz2p30="<<rmsz2p30<<
331 "rmsy2p30R="<<rmsy2p30R<<
332 "rmsz2p30R="<<rmsz2p30R<<
336 if (rdistance2>16) return 3;
339 if ((rdistancey2>9. || rdistancez2>9.) && cluster->GetType()==0)
340 return 2; //suspisiouce - will be changed
342 if ((rdistancey2>6.25 || rdistancez2>6.25) && cluster->GetType()>0)
343 // strict cut on overlaped cluster
344 return 2; //suspisiouce - will be changed
346 if ( (rdistancey2>1. || rdistancez2>6.25 )
347 && cluster->GetType()<0){
348 seed->SetNFoundable(seed->GetNFoundable()-1);
357 //_____________________________________________________________________________
358 AliTPCtrackerMI::AliTPCtrackerMI(const AliTPCParam *par):
360 fkNIS(par->GetNInnerSector()/2),
362 fkNOS(par->GetNOuterSector()/2),
379 //---------------------------------------------------------------------
380 // The main TPC tracker constructor
381 //---------------------------------------------------------------------
382 fInnerSec=new AliTPCtrackerSector[fkNIS];
383 fOuterSec=new AliTPCtrackerSector[fkNOS];
386 for (i=0; i<fkNIS; i++) fInnerSec[i].Setup(par,0);
387 for (i=0; i<fkNOS; i++) fOuterSec[i].Setup(par,1);
390 Int_t nrowlow = par->GetNRowLow();
391 Int_t nrowup = par->GetNRowUp();
394 for (Int_t i=0;i<nrowlow;i++){
395 fXRow[i] = par->GetPadRowRadiiLow(i);
396 fPadLength[i]= par->GetPadPitchLength(0,i);
397 fYMax[i] = fXRow[i]*TMath::Tan(0.5*par->GetInnerAngle());
401 for (Int_t i=0;i<nrowup;i++){
402 fXRow[i+nrowlow] = par->GetPadRowRadiiUp(i);
403 fPadLength[i+nrowlow] = par->GetPadPitchLength(60,i);
404 fYMax[i+nrowlow] = fXRow[i+nrowlow]*TMath::Tan(0.5*par->GetOuterAngle());
407 fDebugStreamer = new TTreeSRedirector("TPCdebug.root");
409 //________________________________________________________________________
410 AliTPCtrackerMI::AliTPCtrackerMI(const AliTPCtrackerMI &t):
431 //------------------------------------
432 // dummy copy constructor
433 //------------------------------------------------------------------
436 AliTPCtrackerMI & AliTPCtrackerMI::operator=(const AliTPCtrackerMI& /*r*/){
437 //------------------------------
439 //--------------------------------------------------------------
442 //_____________________________________________________________________________
443 AliTPCtrackerMI::~AliTPCtrackerMI() {
444 //------------------------------------------------------------------
445 // TPC tracker destructor
446 //------------------------------------------------------------------
453 if (fDebugStreamer) delete fDebugStreamer;
457 void AliTPCtrackerMI::FillESD(TObjArray* arr)
461 //fill esds using updated tracks
463 // write tracks to the event
464 // store index of the track
465 Int_t nseed=arr->GetEntriesFast();
466 //FindKinks(arr,fEvent);
467 for (Int_t i=0; i<nseed; i++) {
468 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
472 if (AliTPCReconstructor::StreamLevel()>1) {
473 (*fDebugStreamer)<<"Track0"<<
477 // pt->PropagateTo(fParam->GetInnerRadiusLow());
478 if (pt->GetKinkIndex(0)<=0){ //don't propagate daughter tracks
479 pt->PropagateTo(fParam->GetInnerRadiusLow());
482 if (( pt->GetPoints()[2]- pt->GetPoints()[0])>5 && pt->GetPoints()[3]>0.8){
484 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
485 iotrack.SetTPCPoints(pt->GetPoints());
486 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
487 iotrack.SetV0Indexes(pt->GetV0Indexes());
488 // iotrack.SetTPCpid(pt->fTPCr);
489 //iotrack.SetTPCindex(i);
490 fEvent->AddTrack(&iotrack);
494 if ( (pt->GetNumberOfClusters()>70)&& (Float_t(pt->GetNumberOfClusters())/Float_t(pt->GetNFoundable()))>0.55) {
496 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
497 iotrack.SetTPCPoints(pt->GetPoints());
498 //iotrack.SetTPCindex(i);
499 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
500 iotrack.SetV0Indexes(pt->GetV0Indexes());
501 // iotrack.SetTPCpid(pt->fTPCr);
502 fEvent->AddTrack(&iotrack);
506 // short tracks - maybe decays
508 if ( (pt->GetNumberOfClusters()>30) && (Float_t(pt->GetNumberOfClusters())/Float_t(pt->GetNFoundable()))>0.70) {
509 Int_t found,foundable,shared;
510 pt->GetClusterStatistic(0,60,found, foundable,shared,kFALSE);
511 if ( (found>20) && (pt->GetNShared()/float(pt->GetNumberOfClusters())<0.2)){
513 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
514 //iotrack.SetTPCindex(i);
515 iotrack.SetTPCPoints(pt->GetPoints());
516 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
517 iotrack.SetV0Indexes(pt->GetV0Indexes());
518 //iotrack.SetTPCpid(pt->fTPCr);
519 fEvent->AddTrack(&iotrack);
524 if ( (pt->GetNumberOfClusters()>20) && (Float_t(pt->GetNumberOfClusters())/Float_t(pt->GetNFoundable()))>0.8) {
525 Int_t found,foundable,shared;
526 pt->GetClusterStatistic(0,60,found, foundable,shared,kFALSE);
527 if (found<20) continue;
528 if (pt->GetNShared()/float(pt->GetNumberOfClusters())>0.2) continue;
531 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
532 iotrack.SetTPCPoints(pt->GetPoints());
533 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
534 iotrack.SetV0Indexes(pt->GetV0Indexes());
535 //iotrack.SetTPCpid(pt->fTPCr);
536 //iotrack.SetTPCindex(i);
537 fEvent->AddTrack(&iotrack);
540 // short tracks - secondaties
542 if ( (pt->GetNumberOfClusters()>30) ) {
543 Int_t found,foundable,shared;
544 pt->GetClusterStatistic(128,158,found, foundable,shared,kFALSE);
545 if ( (found>20) && (pt->GetNShared()/float(pt->GetNumberOfClusters())<0.2) &&float(found)/float(foundable)>0.8){
547 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
548 iotrack.SetTPCPoints(pt->GetPoints());
549 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
550 iotrack.SetV0Indexes(pt->GetV0Indexes());
551 //iotrack.SetTPCpid(pt->fTPCr);
552 //iotrack.SetTPCindex(i);
553 fEvent->AddTrack(&iotrack);
558 if ( (pt->GetNumberOfClusters()>15)) {
559 Int_t found,foundable,shared;
560 pt->GetClusterStatistic(138,158,found, foundable,shared,kFALSE);
561 if (found<15) continue;
562 if (foundable<=0) continue;
563 if (pt->GetNShared()/float(pt->GetNumberOfClusters())>0.2) continue;
564 if (float(found)/float(foundable)<0.8) continue;
567 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
568 iotrack.SetTPCPoints(pt->GetPoints());
569 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
570 iotrack.SetV0Indexes(pt->GetV0Indexes());
571 // iotrack.SetTPCpid(pt->fTPCr);
572 //iotrack.SetTPCindex(i);
573 fEvent->AddTrack(&iotrack);
578 printf("Number of filled ESDs-\t%d\n",fEvent->GetNumberOfTracks());
585 Double_t AliTPCtrackerMI::ErrY2(AliTPCseed* seed, AliTPCclusterMI * cl){
588 // Use calibrated cluster error from OCDB
590 AliTPCClusterParam * clparam = AliTPCcalibDB::Instance()->GetClusterParam();
592 Float_t z = TMath::Abs(fParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
593 Int_t ctype = cl->GetType();
594 Int_t type = (cl->GetRow()<63) ? 0: (cl->GetRow()>126) ? 1:2;
595 Double_t angle = seed->GetSnp()*seed->GetSnp();
596 angle = TMath::Sqrt(TMath::Abs(angle/(1.-angle)));
597 Double_t erry2 = clparam->GetError0Par(0,type, z,angle);
599 erry2+=0.5; // edge cluster
602 seed->SetErrorY2(erry2);
606 //calculate look-up table at the beginning
607 // static Bool_t ginit = kFALSE;
608 // static Float_t gnoise1,gnoise2,gnoise3;
609 // static Float_t ggg1[10000];
610 // static Float_t ggg2[10000];
611 // static Float_t ggg3[10000];
612 // static Float_t glandau1[10000];
613 // static Float_t glandau2[10000];
614 // static Float_t glandau3[10000];
616 // static Float_t gcor01[500];
617 // static Float_t gcor02[500];
618 // static Float_t gcorp[500];
622 // if (ginit==kFALSE){
623 // for (Int_t i=1;i<500;i++){
624 // Float_t rsigma = float(i)/100.;
625 // gcor02[i] = TMath::Max(0.78 +TMath::Exp(7.4*(rsigma-1.2)),0.6);
626 // gcor01[i] = TMath::Max(0.72 +TMath::Exp(3.36*(rsigma-1.2)),0.6);
627 // gcorp[i] = TMath::Max(TMath::Power((rsigma+0.5),1.5),1.2);
631 // for (Int_t i=3;i<10000;i++){
635 // Float_t amp = float(i);
636 // Float_t padlength =0.75;
637 // gnoise1 = 0.0004/padlength;
638 // Float_t nel = 0.268*amp;
639 // Float_t nprim = 0.155*amp;
640 // ggg1[i] = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.001*nel/(padlength*padlength))/nel;
641 // glandau1[i] = (2.+0.12*nprim)*0.5* (2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
642 // if (glandau1[i]>1) glandau1[i]=1;
643 // glandau1[i]*=padlength*padlength/12.;
647 // gnoise2 = 0.0004/padlength;
649 // nprim = 0.133*amp;
650 // ggg2[i] = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
651 // glandau2[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
652 // if (glandau2[i]>1) glandau2[i]=1;
653 // glandau2[i]*=padlength*padlength/12.;
658 // gnoise3 = 0.0004/padlength;
660 // nprim = 0.133*amp;
661 // ggg3[i] = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
662 // glandau3[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
663 // if (glandau3[i]>1) glandau3[i]=1;
664 // glandau3[i]*=padlength*padlength/12.;
672 // Int_t amp = int(TMath::Abs(cl->GetQ()));
674 // seed->SetErrorY2(1.);
678 // Float_t z = TMath::Abs(fParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
679 // Int_t ctype = cl->GetType();
680 // Float_t padlength= GetPadPitchLength(seed->GetRow());
681 // Double_t angle2 = seed->GetSnp()*seed->GetSnp();
682 // angle2 = angle2/(1-angle2);
684 // //cluster "quality"
685 // Int_t rsigmay = int(100.*cl->GetSigmaY2()/(seed->GetCurrentSigmaY2()));
688 // if (fSectors==fInnerSec){
689 // snoise2 = gnoise1;
690 // res = ggg1[amp]*z+glandau1[amp]*angle2;
691 // if (ctype==0) res *= gcor01[rsigmay];
694 // res*= gcorp[rsigmay];
698 // if (padlength<1.1){
699 // snoise2 = gnoise2;
700 // res = ggg2[amp]*z+glandau2[amp]*angle2;
701 // if (ctype==0) res *= gcor02[rsigmay];
704 // res*= gcorp[rsigmay];
708 // snoise2 = gnoise3;
709 // res = ggg3[amp]*z+glandau3[amp]*angle2;
710 // if (ctype==0) res *= gcor02[rsigmay];
713 // res*= gcorp[rsigmay];
720 // res*=2.4; // overestimate error 2 times
724 // if (res<2*snoise2)
727 // seed->SetErrorY2(res);
735 Double_t AliTPCtrackerMI::ErrZ2(AliTPCseed* seed, AliTPCclusterMI * cl){
738 // Use calibrated cluster error from OCDB
740 AliTPCClusterParam * clparam = AliTPCcalibDB::Instance()->GetClusterParam();
742 Float_t z = TMath::Abs(fParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
743 Int_t ctype = cl->GetType();
744 Int_t type = (cl->GetRow()<63) ? 0: (cl->GetRow()>126) ? 1:2;
746 Double_t angle2 = seed->GetSnp()*seed->GetSnp();
747 angle2 = seed->GetTgl()*seed->GetTgl()*(1+angle2/(1-angle2));
748 Double_t angle = TMath::Sqrt(TMath::Abs(angle2));
749 Double_t errz2 = clparam->GetError0Par(1,type, z,angle);
751 errz2+=0.5; // edge cluster
754 seed->SetErrorZ2(errz2);
760 // //seed->SetErrorY2(0.1);
762 // //calculate look-up table at the beginning
763 // static Bool_t ginit = kFALSE;
764 // static Float_t gnoise1,gnoise2,gnoise3;
765 // static Float_t ggg1[10000];
766 // static Float_t ggg2[10000];
767 // static Float_t ggg3[10000];
768 // static Float_t glandau1[10000];
769 // static Float_t glandau2[10000];
770 // static Float_t glandau3[10000];
772 // static Float_t gcor01[1000];
773 // static Float_t gcor02[1000];
774 // static Float_t gcorp[1000];
778 // if (ginit==kFALSE){
779 // for (Int_t i=1;i<1000;i++){
780 // Float_t rsigma = float(i)/100.;
781 // gcor02[i] = TMath::Max(0.81 +TMath::Exp(6.8*(rsigma-1.2)),0.6);
782 // gcor01[i] = TMath::Max(0.72 +TMath::Exp(2.04*(rsigma-1.2)),0.6);
783 // gcorp[i] = TMath::Max(TMath::Power((rsigma+0.5),1.5),1.2);
787 // for (Int_t i=3;i<10000;i++){
791 // Float_t amp = float(i);
792 // Float_t padlength =0.75;
793 // gnoise1 = 0.0004/padlength;
794 // Float_t nel = 0.268*amp;
795 // Float_t nprim = 0.155*amp;
796 // ggg1[i] = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.001*nel/(padlength*padlength))/nel;
797 // glandau1[i] = (2.+0.12*nprim)*0.5* (2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
798 // if (glandau1[i]>1) glandau1[i]=1;
799 // glandau1[i]*=padlength*padlength/12.;
803 // gnoise2 = 0.0004/padlength;
805 // nprim = 0.133*amp;
806 // ggg2[i] = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
807 // glandau2[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
808 // if (glandau2[i]>1) glandau2[i]=1;
809 // glandau2[i]*=padlength*padlength/12.;
814 // gnoise3 = 0.0004/padlength;
816 // nprim = 0.133*amp;
817 // ggg3[i] = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
818 // glandau3[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
819 // if (glandau3[i]>1) glandau3[i]=1;
820 // glandau3[i]*=padlength*padlength/12.;
828 // Int_t amp = int(TMath::Abs(cl->GetQ()));
830 // seed->SetErrorY2(1.);
834 // Float_t z = TMath::Abs(fParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
835 // Int_t ctype = cl->GetType();
836 // Float_t padlength= GetPadPitchLength(seed->GetRow());
838 // Double_t angle2 = seed->GetSnp()*seed->GetSnp();
839 // // if (angle2<0.6) angle2 = 0.6;
840 // angle2 = seed->GetTgl()*seed->GetTgl()*(1+angle2/(1-angle2));
842 // //cluster "quality"
843 // Int_t rsigmaz = int(100.*cl->GetSigmaZ2()/(seed->GetCurrentSigmaZ2()));
846 // if (fSectors==fInnerSec){
847 // snoise2 = gnoise1;
848 // res = ggg1[amp]*z+glandau1[amp]*angle2;
849 // if (ctype==0) res *= gcor01[rsigmaz];
852 // res*= gcorp[rsigmaz];
856 // if (padlength<1.1){
857 // snoise2 = gnoise2;
858 // res = ggg2[amp]*z+glandau2[amp]*angle2;
859 // if (ctype==0) res *= gcor02[rsigmaz];
862 // res*= gcorp[rsigmaz];
866 // snoise2 = gnoise3;
867 // res = ggg3[amp]*z+glandau3[amp]*angle2;
868 // if (ctype==0) res *= gcor02[rsigmaz];
871 // res*= gcorp[rsigmaz];
880 // if ((ctype<0) &&<70){
885 // if (res<2*snoise2)
887 // if (res>3) res =3;
888 // seed->SetErrorZ2(res);
896 void AliTPCtrackerMI::RotateToLocal(AliTPCseed *seed)
898 //rotate to track "local coordinata
899 Float_t x = seed->GetX();
900 Float_t y = seed->GetY();
901 Float_t ymax = x*TMath::Tan(0.5*fSectors->GetAlpha());
904 seed->SetRelativeSector((seed->GetRelativeSector()+1) % fN);
905 if (!seed->Rotate(fSectors->GetAlpha()))
907 } else if (y <-ymax) {
908 seed->SetRelativeSector((seed->GetRelativeSector()-1+fN) % fN);
909 if (!seed->Rotate(-fSectors->GetAlpha()))
917 //_____________________________________________________________________________
918 Double_t AliTPCtrackerMI::F1old(Double_t x1,Double_t y1,
919 Double_t x2,Double_t y2,
920 Double_t x3,Double_t y3)
922 //-----------------------------------------------------------------
923 // Initial approximation of the track curvature
924 //-----------------------------------------------------------------
925 Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
926 Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
927 (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
928 Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
929 (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
931 Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
932 if ( xr*xr+yr*yr<=0.00000000000001) return 100;
933 return -xr*yr/sqrt(xr*xr+yr*yr);
938 //_____________________________________________________________________________
939 Double_t AliTPCtrackerMI::F1(Double_t x1,Double_t y1,
940 Double_t x2,Double_t y2,
941 Double_t x3,Double_t y3)
943 //-----------------------------------------------------------------
944 // Initial approximation of the track curvature
945 //-----------------------------------------------------------------
951 Double_t det = x3*y2-x2*y3;
956 Double_t u = 0.5* (x2*(x2-x3)+y2*(y2-y3))/det;
957 Double_t x0 = x3*0.5-y3*u;
958 Double_t y0 = y3*0.5+x3*u;
959 Double_t c2 = 1/TMath::Sqrt(x0*x0+y0*y0);
965 Double_t AliTPCtrackerMI::F2(Double_t x1,Double_t y1,
966 Double_t x2,Double_t y2,
967 Double_t x3,Double_t y3)
969 //-----------------------------------------------------------------
970 // Initial approximation of the track curvature
971 //-----------------------------------------------------------------
977 Double_t det = x3*y2-x2*y3;
982 Double_t u = 0.5* (x2*(x2-x3)+y2*(y2-y3))/det;
983 Double_t x0 = x3*0.5-y3*u;
984 Double_t y0 = y3*0.5+x3*u;
985 Double_t c2 = 1/TMath::Sqrt(x0*x0+y0*y0);
994 //_____________________________________________________________________________
995 Double_t AliTPCtrackerMI::F2old(Double_t x1,Double_t y1,
996 Double_t x2,Double_t y2,
997 Double_t x3,Double_t y3)
999 //-----------------------------------------------------------------
1000 // Initial approximation of the track curvature times center of curvature
1001 //-----------------------------------------------------------------
1002 Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
1003 Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
1004 (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
1005 Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
1006 (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
1008 Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
1010 return -a/(d*y1-b)*xr/sqrt(xr*xr+yr*yr);
1013 //_____________________________________________________________________________
1014 Double_t AliTPCtrackerMI::F3(Double_t x1,Double_t y1,
1015 Double_t x2,Double_t y2,
1016 Double_t z1,Double_t z2)
1018 //-----------------------------------------------------------------
1019 // Initial approximation of the tangent of the track dip angle
1020 //-----------------------------------------------------------------
1021 return (z1 - z2)/sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
1025 Double_t AliTPCtrackerMI::F3n(Double_t x1,Double_t y1,
1026 Double_t x2,Double_t y2,
1027 Double_t z1,Double_t z2, Double_t c)
1029 //-----------------------------------------------------------------
1030 // Initial approximation of the tangent of the track dip angle
1031 //-----------------------------------------------------------------
1035 //angle1 = (z1-z2)*c/(TMath::ASin(c*x1-ni)-TMath::ASin(c*x2-ni));
1037 Double_t d = TMath::Sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
1038 if (TMath::Abs(d*c*0.5)>1) return 0;
1039 // Double_t angle2 = TMath::ASin(d*c*0.5);
1040 // Double_t angle2 = AliTPCFastMath::FastAsin(d*c*0.5);
1041 Double_t angle2 = (d*c*0.5>0.1)? TMath::ASin(d*c*0.5): AliTPCFastMath::FastAsin(d*c*0.5);
1043 angle2 = (z1-z2)*c/(angle2*2.);
1047 Bool_t AliTPCtrackerMI::GetProlongation(Double_t x1, Double_t x2, Double_t x[5], Double_t &y, Double_t &z)
1048 {//-----------------------------------------------------------------
1049 // This function find proloncation of a track to a reference plane x=x2.
1050 //-----------------------------------------------------------------
1054 if (TMath::Abs(x[4]*x1 - x[2]) >= 0.999) {
1058 Double_t c1=x[4]*x1 - x[2], r1=sqrt(1.- c1*c1);
1059 Double_t c2=x[4]*x2 - x[2], r2=sqrt(1.- c2*c2);
1063 Double_t dy = dx*(c1+c2)/(r1+r2);
1066 Double_t delta = x[4]*dx*(c1+c2)/(c1*r2 + c2*r1);
1068 if (TMath::Abs(delta)>0.01){
1069 dz = x[3]*TMath::ASin(delta)/x[4];
1071 dz = x[3]*AliTPCFastMath::FastAsin(delta)/x[4];
1074 //dz = x[3]*AliTPCFastMath::FastAsin(delta)/x[4];
1082 Int_t AliTPCtrackerMI::LoadClusters (TTree *tree)
1087 return LoadClusters();
1091 Int_t AliTPCtrackerMI::LoadClusters(TObjArray *arr)
1094 // load clusters to the memory
1095 AliTPCClustersRow *clrow = 0x0;
1096 Int_t lower = arr->LowerBound();
1097 Int_t entries = arr->GetEntriesFast();
1099 for (Int_t i=lower; i<entries; i++) {
1100 clrow = (AliTPCClustersRow*) arr->At(i);
1101 if(!clrow) continue;
1102 if(!clrow->GetArray()) continue;
1106 fParam->AdjustSectorRow(clrow->GetID(),sec,row);
1108 for (Int_t icl=0; icl<clrow->GetArray()->GetEntriesFast(); icl++){
1109 Transform((AliTPCclusterMI*)(clrow->GetArray()->At(icl)));
1112 if (clrow->GetArray()->GetEntriesFast()<=0) continue;
1113 AliTPCtrackerRow * tpcrow=0;
1116 tpcrow = &(fInnerSec[sec%fkNIS][row]);
1120 tpcrow = &(fOuterSec[(sec-fkNIS*2)%fkNOS][row]);
1121 left = (sec-fkNIS*2)/fkNOS;
1124 tpcrow->SetN1(clrow->GetArray()->GetEntriesFast());
1125 tpcrow->SetClusters1(new AliTPCclusterMI[tpcrow->GetN1()]);
1126 for (Int_t i=0;i<tpcrow->GetN1();i++)
1127 tpcrow->SetCluster1(i, *(AliTPCclusterMI*)(clrow->GetArray()->At(i)));
1130 tpcrow->SetN2(clrow->GetArray()->GetEntriesFast());
1131 tpcrow->SetClusters2(new AliTPCclusterMI[tpcrow->GetN2()]);
1132 for (Int_t i=0;i<tpcrow->GetN2();i++)
1133 tpcrow->SetCluster2(i,*(AliTPCclusterMI*)(clrow->GetArray()->At(i)));
1144 Int_t AliTPCtrackerMI::LoadClusters()
1147 // load clusters to the memory
1148 AliTPCClustersRow *clrow= new AliTPCClustersRow;
1149 clrow->SetClass("AliTPCclusterMI");
1151 clrow->GetArray()->ExpandCreateFast(10000);
1153 // TTree * tree = fClustersArray.GetTree();
1155 TTree * tree = fInput;
1156 TBranch * br = tree->GetBranch("Segment");
1157 br->SetAddress(&clrow);
1159 Int_t j=Int_t(tree->GetEntries());
1160 for (Int_t i=0; i<j; i++) {
1164 fParam->AdjustSectorRow(clrow->GetID(),sec,row);
1165 for (Int_t icl=0; icl<clrow->GetArray()->GetEntriesFast(); icl++){
1166 Transform((AliTPCclusterMI*)(clrow->GetArray()->At(icl)));
1169 AliTPCtrackerRow * tpcrow=0;
1172 tpcrow = &(fInnerSec[sec%fkNIS][row]);
1176 tpcrow = &(fOuterSec[(sec-fkNIS*2)%fkNOS][row]);
1177 left = (sec-fkNIS*2)/fkNOS;
1180 tpcrow->SetN1(clrow->GetArray()->GetEntriesFast());
1181 tpcrow->SetClusters1(new AliTPCclusterMI[tpcrow->GetN1()]);
1182 for (Int_t i=0;i<tpcrow->GetN1();i++)
1183 tpcrow->SetCluster1(i, *(AliTPCclusterMI*)(clrow->GetArray()->At(i)));
1186 tpcrow->SetN2(clrow->GetArray()->GetEntriesFast());
1187 tpcrow->SetClusters2(new AliTPCclusterMI[tpcrow->GetN2()]);
1188 for (Int_t i=0;i<tpcrow->GetN2();i++)
1189 tpcrow->SetCluster2(i,*(AliTPCclusterMI*)(clrow->GetArray()->At(i)));
1200 void AliTPCtrackerMI::UnloadClusters()
1203 // unload clusters from the memory
1205 Int_t nrows = fOuterSec->GetNRows();
1206 for (Int_t sec = 0;sec<fkNOS;sec++)
1207 for (Int_t row = 0;row<nrows;row++){
1208 AliTPCtrackerRow* tpcrow = &(fOuterSec[sec%fkNOS][row]);
1210 // if (tpcrow->fClusters1) delete []tpcrow->fClusters1;
1211 // if (tpcrow->fClusters2) delete []tpcrow->fClusters2;
1213 tpcrow->ResetClusters();
1216 nrows = fInnerSec->GetNRows();
1217 for (Int_t sec = 0;sec<fkNIS;sec++)
1218 for (Int_t row = 0;row<nrows;row++){
1219 AliTPCtrackerRow* tpcrow = &(fInnerSec[sec%fkNIS][row]);
1221 // if (tpcrow->fClusters1) delete []tpcrow->fClusters1;
1222 //if (tpcrow->fClusters2) delete []tpcrow->fClusters2;
1224 tpcrow->ResetClusters();
1230 void AliTPCtrackerMI::FillClusterArray(TObjArray* array) const{
1232 // Filling cluster to the array - For visualization purposes
1235 nrows = fOuterSec->GetNRows();
1236 for (Int_t sec = 0;sec<fkNOS;sec++)
1237 for (Int_t row = 0;row<nrows;row++){
1238 AliTPCtrackerRow* tpcrow = &(fOuterSec[sec%fkNOS][row]);
1239 if (!tpcrow) continue;
1240 for (Int_t icl = 0;icl<tpcrow->GetN();icl++){
1241 array->AddLast((TObject*)((*tpcrow)[icl]));
1244 nrows = fInnerSec->GetNRows();
1245 for (Int_t sec = 0;sec<fkNIS;sec++)
1246 for (Int_t row = 0;row<nrows;row++){
1247 AliTPCtrackerRow* tpcrow = &(fInnerSec[sec%fkNIS][row]);
1248 if (!tpcrow) continue;
1249 for (Int_t icl = 0;icl<tpcrow->GetN();icl++){
1250 array->AddLast((TObject*)(*tpcrow)[icl]);
1256 void AliTPCtrackerMI::Transform(AliTPCclusterMI * cluster){
1260 AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ;
1262 AliFatal("Tranformations not in calibDB");
1264 Double_t x[3]={cluster->GetRow(),cluster->GetPad(),cluster->GetTimeBin()};
1265 Int_t i[1]={cluster->GetDetector()};
1266 transform->Transform(x,i,0,1);
1267 if (cluster->GetDetector()%36>17){
1272 // in debug mode check the transformation
1274 if (AliTPCReconstructor::StreamLevel()>1) {
1276 cluster->GetGlobalXYZ(gx);
1277 Int_t event = (fEvent==NULL)? 0: fEvent->GetEventNumberInFile();
1278 TTreeSRedirector &cstream = *fDebugStreamer;
1279 cstream<<"Transform"<<
1290 cluster->SetX(x[0]);
1291 cluster->SetY(x[1]);
1292 cluster->SetZ(x[2]);
1298 //if (!fParam->IsGeoRead()) fParam->ReadGeoMatrices();
1299 TGeoHMatrix *mat = fParam->GetClusterMatrix(cluster->GetDetector());
1301 Double_t pos[3]= {cluster->GetX(),cluster->GetY(),cluster->GetZ()};
1302 Double_t posC[3]={cluster->GetX(),cluster->GetY(),cluster->GetZ()};
1303 if (mat) mat->LocalToMaster(pos,posC);
1305 // chack Loading of Geo matrices from GeoManager - TEMPORARY FIX
1307 cluster->SetX(posC[0]);
1308 cluster->SetY(posC[1]);
1309 cluster->SetZ(posC[2]);
1312 //_____________________________________________________________________________
1313 Int_t AliTPCtrackerMI::LoadOuterSectors() {
1314 //-----------------------------------------------------------------
1315 // This function fills outer TPC sectors with clusters.
1316 //-----------------------------------------------------------------
1317 Int_t nrows = fOuterSec->GetNRows();
1319 for (Int_t sec = 0;sec<fkNOS;sec++)
1320 for (Int_t row = 0;row<nrows;row++){
1321 AliTPCtrackerRow* tpcrow = &(fOuterSec[sec%fkNOS][row]);
1322 Int_t sec2 = sec+2*fkNIS;
1324 Int_t ncl = tpcrow->GetN1();
1326 AliTPCclusterMI *c= (tpcrow->GetCluster1(ncl));
1327 index=(((sec2<<8)+row)<<16)+ncl;
1328 tpcrow->InsertCluster(c,index);
1331 ncl = tpcrow->GetN2();
1333 AliTPCclusterMI *c= (tpcrow->GetCluster2(ncl));
1334 index=((((sec2+fkNOS)<<8)+row)<<16)+ncl;
1335 tpcrow->InsertCluster(c,index);
1338 // write indexes for fast acces
1340 for (Int_t i=0;i<510;i++)
1341 tpcrow->SetFastCluster(i,-1);
1342 for (Int_t i=0;i<tpcrow->GetN();i++){
1343 Int_t zi = Int_t((*tpcrow)[i]->GetZ()+255.);
1344 tpcrow->SetFastCluster(zi,i); // write index
1347 for (Int_t i=0;i<510;i++){
1348 if (tpcrow->GetFastCluster(i)<0)
1349 tpcrow->SetFastCluster(i,last);
1351 last = tpcrow->GetFastCluster(i);
1360 //_____________________________________________________________________________
1361 Int_t AliTPCtrackerMI::LoadInnerSectors() {
1362 //-----------------------------------------------------------------
1363 // This function fills inner TPC sectors with clusters.
1364 //-----------------------------------------------------------------
1365 Int_t nrows = fInnerSec->GetNRows();
1367 for (Int_t sec = 0;sec<fkNIS;sec++)
1368 for (Int_t row = 0;row<nrows;row++){
1369 AliTPCtrackerRow* tpcrow = &(fInnerSec[sec%fkNIS][row]);
1372 Int_t ncl = tpcrow->GetN1();
1374 AliTPCclusterMI *c= (tpcrow->GetCluster1(ncl));
1375 index=(((sec<<8)+row)<<16)+ncl;
1376 tpcrow->InsertCluster(c,index);
1379 ncl = tpcrow->GetN2();
1381 AliTPCclusterMI *c= (tpcrow->GetCluster2(ncl));
1382 index=((((sec+fkNIS)<<8)+row)<<16)+ncl;
1383 tpcrow->InsertCluster(c,index);
1386 // write indexes for fast acces
1388 for (Int_t i=0;i<510;i++)
1389 tpcrow->SetFastCluster(i,-1);
1390 for (Int_t i=0;i<tpcrow->GetN();i++){
1391 Int_t zi = Int_t((*tpcrow)[i]->GetZ()+255.);
1392 tpcrow->SetFastCluster(zi,i); // write index
1395 for (Int_t i=0;i<510;i++){
1396 if (tpcrow->GetFastCluster(i)<0)
1397 tpcrow->SetFastCluster(i,last);
1399 last = tpcrow->GetFastCluster(i);
1411 //_________________________________________________________________________
1412 AliTPCclusterMI *AliTPCtrackerMI::GetClusterMI(Int_t index) const {
1413 //--------------------------------------------------------------------
1414 // Return pointer to a given cluster
1415 //--------------------------------------------------------------------
1416 if (index<0) return 0; // no cluster
1417 Int_t sec=(index&0xff000000)>>24;
1418 Int_t row=(index&0x00ff0000)>>16;
1419 Int_t ncl=(index&0x00007fff)>>00;
1421 const AliTPCtrackerRow * tpcrow=0;
1422 AliTPCclusterMI * clrow =0;
1424 if (sec<0 || sec>=fkNIS*4) {
1425 AliWarning(Form("Wrong sector %d",sec));
1430 tpcrow = &(fInnerSec[sec%fkNIS][row]);
1431 if (tpcrow==0) return 0;
1434 if (tpcrow->GetN1()<=ncl) return 0;
1435 clrow = tpcrow->GetClusters1();
1438 if (tpcrow->GetN2()<=ncl) return 0;
1439 clrow = tpcrow->GetClusters2();
1443 tpcrow = &(fOuterSec[(sec-fkNIS*2)%fkNOS][row]);
1444 if (tpcrow==0) return 0;
1446 if (sec-2*fkNIS<fkNOS) {
1447 if (tpcrow->GetN1()<=ncl) return 0;
1448 clrow = tpcrow->GetClusters1();
1451 if (tpcrow->GetN2()<=ncl) return 0;
1452 clrow = tpcrow->GetClusters2();
1456 return &(clrow[ncl]);
1462 Int_t AliTPCtrackerMI::FollowToNext(AliTPCseed& t, Int_t nr) {
1463 //-----------------------------------------------------------------
1464 // This function tries to find a track prolongation to next pad row
1465 //-----------------------------------------------------------------
1467 Double_t x= GetXrow(nr), ymax=GetMaxY(nr);
1468 AliTPCclusterMI *cl=0;
1469 Int_t tpcindex= t.GetClusterIndex2(nr);
1471 // update current shape info every 5 pad-row
1472 // if ( (nr%5==0) || t.GetNumberOfClusters()<2 || (t.fCurrentSigmaY2<0.0001) ){
1476 if (fIteration>0 && tpcindex>=-1){ //if we have already clusters
1478 if (tpcindex==-1) return 0; //track in dead zone
1480 cl = t.GetClusterPointer(nr);
1481 if ( (cl==0) ) cl = GetClusterMI(tpcindex);
1482 t.SetCurrentClusterIndex1(tpcindex);
1485 Int_t relativesector = ((tpcindex&0xff000000)>>24)%18; // if previously accepted cluster in different sector
1486 Float_t angle = relativesector*fSectors->GetAlpha()+fSectors->GetAlphaShift();
1488 if (angle<-TMath::Pi()) angle += 2*TMath::Pi();
1489 if (angle>=TMath::Pi()) angle -= 2*TMath::Pi();
1491 if (TMath::Abs(angle-t.GetAlpha())>0.001){
1492 Double_t rotation = angle-t.GetAlpha();
1493 t.SetRelativeSector(relativesector);
1494 if (!t.Rotate(rotation)) return 0;
1496 if (!t.PropagateTo(x)) return 0;
1498 t.SetCurrentCluster(cl);
1500 Int_t accept = AcceptCluster(&t,t.GetCurrentCluster());
1501 if ((tpcindex&0x8000)==0) accept =0;
1503 //if founded cluster is acceptible
1504 if (cl->IsUsed(11)) { // id cluster is shared inrease uncertainty
1505 t.SetErrorY2(t.GetErrorY2()+0.03);
1506 t.SetErrorZ2(t.GetErrorZ2()+0.03);
1507 t.SetErrorY2(t.GetErrorY2()*3);
1508 t.SetErrorZ2(t.GetErrorZ2()*3);
1510 t.SetNFoundable(t.GetNFoundable()+1);
1511 UpdateTrack(&t,accept);
1516 if (TMath::Abs(t.GetSnp())>AliTPCReconstructor::GetMaxSnpTracker()) return 0; // cut on angle
1518 // not look for new cluster during refitting
1519 t.SetNFoundable(t.GetNFoundable()+1);
1524 // if (TMath::Abs(t.GetSnp())>0.95 || TMath::Abs(x*t.GetC()-t.GetEta())>0.95) return 0;// patch 28 fev 06
1525 Double_t y=t.GetYat(x);
1526 if (TMath::Abs(y)>ymax){
1528 t.SetRelativeSector((t.GetRelativeSector()+1) % fN);
1529 if (!t.Rotate(fSectors->GetAlpha()))
1531 } else if (y <-ymax) {
1532 t.SetRelativeSector((t.GetRelativeSector()-1+fN) % fN);
1533 if (!t.Rotate(-fSectors->GetAlpha()))
1539 if (!t.PropagateTo(x)) {
1540 if (fIteration==0) t.SetRemoval(10);
1544 Double_t z=t.GetZ();
1547 if (!IsActive(t.GetRelativeSector(),nr)) {
1549 t.SetClusterIndex2(nr,-1);
1552 //AliInfo(Form("A - Sector%d phi %f - alpha %f", t.fRelativeSector,y/x, t.GetAlpha()));
1553 Bool_t isActive = IsActive(t.GetRelativeSector(),nr);
1554 Bool_t isActive2 = (nr>=fInnerSec->GetNRows()) ? fOuterSec[t.GetRelativeSector()][nr-fInnerSec->GetNRows()].GetN()>0:fInnerSec[t.GetRelativeSector()][nr].GetN()>0;
1556 if (!isActive || !isActive2) return 0;
1558 const AliTPCtrackerRow &krow=GetRow(t.GetRelativeSector(),nr);
1559 if ( (t.GetSigmaY2()<0) || t.GetSigmaZ2()<0) return 0;
1561 Double_t roadz = 1.;
1563 if (TMath::Abs(TMath::Abs(y)-ymax)<krow.GetDeadZone()){
1565 t.SetClusterIndex2(nr,-1);
1570 if (TMath::Abs(z)<(AliTPCReconstructor::GetCtgRange()*x+10) && TMath::Abs(z)<fParam->GetZLength(0) && (TMath::Abs(t.GetSnp())<AliTPCReconstructor::GetMaxSnpTracker()))
1571 t.SetNFoundable(t.GetNFoundable()+1);
1577 // cl = krow.FindNearest2(y+10.,z,roady,roadz,index);
1578 cl = krow.FindNearest2(y,z,roady,roadz,index);
1579 if (cl) t.SetCurrentClusterIndex1(krow.GetIndex(index));
1582 t.SetCurrentCluster(cl);
1584 if (fIteration==2&&cl->IsUsed(10)) return 0;
1585 Int_t accept = AcceptCluster(&t,t.GetCurrentCluster());
1586 if (fIteration==2&&cl->IsUsed(11)) {
1587 t.SetErrorY2(t.GetErrorY2()+0.03);
1588 t.SetErrorZ2(t.GetErrorZ2()+0.03);
1589 t.SetErrorY2(t.GetErrorY2()*3);
1590 t.SetErrorZ2(t.GetErrorZ2()*3);
1593 if (t.fCurrentCluster->IsUsed(10)){
1598 if (t.fNShared>0.7*t.GetNumberOfClusters()) {
1604 if (accept<3) UpdateTrack(&t,accept);
1607 if ( fIteration==0 && t.GetNFoundable()*0.5 > t.GetNumberOfClusters()) t.SetRemoval(10);
1613 Int_t AliTPCtrackerMI::FollowToNextFast(AliTPCseed& t, Int_t nr) {
1614 //-----------------------------------------------------------------
1615 // This function tries to find a track prolongation to next pad row
1616 //-----------------------------------------------------------------
1618 Double_t x= GetXrow(nr), ymax=GetMaxY(nr);
1620 if (!t.GetProlongation(x,y,z)) {
1626 if (TMath::Abs(y)>ymax){
1629 t.SetRelativeSector((t.GetRelativeSector()+1) % fN);
1630 if (!t.Rotate(fSectors->GetAlpha()))
1632 } else if (y <-ymax) {
1633 t.SetRelativeSector((t.GetRelativeSector()-1+fN) % fN);
1634 if (!t.Rotate(-fSectors->GetAlpha()))
1637 if (!t.PropagateTo(x)) {
1640 t.GetProlongation(x,y,z);
1643 // update current shape info every 2 pad-row
1644 if ( (nr%2==0) || t.GetNumberOfClusters()<2 || (t.GetCurrentSigmaY2()<0.0001) ){
1645 // t.fCurrentSigmaY = GetSigmaY(&t);
1646 //t.fCurrentSigmaZ = GetSigmaZ(&t);
1650 AliTPCclusterMI *cl=0;
1655 const AliTPCtrackerRow &krow=GetRow(t.GetRelativeSector(),nr);
1656 if ( (t.GetSigmaY2()<0) || t.GetSigmaZ2()<0) return 0;
1658 Double_t roadz = 1.;
1661 if (TMath::Abs(TMath::Abs(y)-ymax)<krow.GetDeadZone()){
1663 t.SetClusterIndex2(row,-1);
1668 if (TMath::Abs(z)>(AliTPCReconstructor::GetCtgRange()*x+10)) t.SetClusterIndex2(row,-1);
1672 if ((cl==0)&&(krow)) {
1673 // cl = krow.FindNearest2(y+10,z,roady,roadz,index);
1674 cl = krow.FindNearest2(y,z,roady,roadz,index);
1676 if (cl) t.SetCurrentClusterIndex1(krow.GetIndex(index));
1680 t.SetCurrentCluster(cl);
1681 // Int_t accept = AcceptCluster(&t,t.fCurrentCluster);
1683 t.SetClusterIndex2(row,index);
1684 t.SetClusterPointer(row, cl);
1692 //_________________________________________________________________________
1693 Bool_t AliTPCtrackerMI::GetTrackPoint(Int_t index, AliTrackPoint &p ) const
1695 // Get track space point by index
1696 // return false in case the cluster doesn't exist
1697 AliTPCclusterMI *cl = GetClusterMI(index);
1698 if (!cl) return kFALSE;
1699 Int_t sector = (index&0xff000000)>>24;
1700 // Int_t row = (index&0x00ff0000)>>16;
1702 // xyz[0] = fParam->GetPadRowRadii(sector,row);
1703 xyz[0] = cl->GetX();
1704 xyz[1] = cl->GetY();
1705 xyz[2] = cl->GetZ();
1707 fParam->AdjustCosSin(sector,cos,sin);
1708 Float_t x = cos*xyz[0]-sin*xyz[1];
1709 Float_t y = cos*xyz[1]+sin*xyz[0];
1711 Float_t sigmaY2 = 0.027*cl->GetSigmaY2();
1712 if (sector < fParam->GetNInnerSector()) sigmaY2 *= 2.07;
1713 Float_t sigmaZ2 = 0.066*cl->GetSigmaZ2();
1714 if (sector < fParam->GetNInnerSector()) sigmaZ2 *= 1.77;
1715 cov[0] = sin*sin*sigmaY2;
1716 cov[1] = -sin*cos*sigmaY2;
1718 cov[3] = cos*cos*sigmaY2;
1721 p.SetXYZ(x,y,xyz[2],cov);
1722 AliGeomManager::ELayerID iLayer;
1724 if (sector < fParam->GetNInnerSector()) {
1725 iLayer = AliGeomManager::kTPC1;
1729 iLayer = AliGeomManager::kTPC2;
1730 idet = sector - fParam->GetNInnerSector();
1732 UShort_t volid = AliGeomManager::LayerToVolUID(iLayer,idet);
1733 p.SetVolumeID(volid);
1739 Int_t AliTPCtrackerMI::UpdateClusters(AliTPCseed& t, Int_t nr) {
1740 //-----------------------------------------------------------------
1741 // This function tries to find a track prolongation to next pad row
1742 //-----------------------------------------------------------------
1743 t.SetCurrentCluster(0);
1744 t.SetCurrentClusterIndex1(0);
1746 Double_t xt=t.GetX();
1747 Int_t row = GetRowNumber(xt)-1;
1748 Double_t ymax= GetMaxY(nr);
1750 if (row < nr) return 1; // don't prolongate if not information until now -
1751 // if (TMath::Abs(t.GetSnp())>0.9 && t.GetNumberOfClusters()>40. && fIteration!=2) {
1753 // return 0; // not prolongate strongly inclined tracks
1755 // if (TMath::Abs(t.GetSnp())>0.95) {
1757 // return 0; // not prolongate strongly inclined tracks
1758 // }// patch 28 fev 06
1760 Double_t x= GetXrow(nr);
1762 //t.PropagateTo(x+0.02);
1763 //t.PropagateTo(x+0.01);
1764 if (!t.PropagateTo(x)){
1771 if (TMath::Abs(y)>ymax){
1773 t.SetRelativeSector((t.GetRelativeSector()+1) % fN);
1774 if (!t.Rotate(fSectors->GetAlpha()))
1776 } else if (y <-ymax) {
1777 t.SetRelativeSector((t.GetRelativeSector()-1+fN) % fN);
1778 if (!t.Rotate(-fSectors->GetAlpha()))
1781 // if (!t.PropagateTo(x)){
1788 if (TMath::Abs(t.GetSnp())>AliTPCReconstructor::GetMaxSnpTracker()) return 0;
1790 if (!IsActive(t.GetRelativeSector(),nr)) {
1792 t.SetClusterIndex2(nr,-1);
1795 //AliInfo(Form("A - Sector%d phi %f - alpha %f", t.fRelativeSector,y/x, t.GetAlpha()));
1797 AliTPCtrackerRow &krow=GetRow(t.GetRelativeSector(),nr);
1799 if (TMath::Abs(TMath::Abs(y)-ymax)<krow.GetDeadZone()){
1801 t.SetClusterIndex2(nr,-1);
1806 if (TMath::Abs(t.GetZ())<(AliTPCReconstructor::GetCtgRange()*t.GetX()+10) && (TMath::Abs(t.GetSnp())<AliTPCReconstructor::GetMaxSnpTracker())) t.SetNFoundable(t.GetNFoundable()+1);
1812 if ( (nr%2==0) || t.GetNumberOfClusters()<2){
1813 // t.fCurrentSigmaY = GetSigmaY(&t);
1814 //t.fCurrentSigmaZ = GetSigmaZ(&t);
1818 AliTPCclusterMI *cl=0;
1821 Double_t roady = 1.;
1822 Double_t roadz = 1.;
1826 index = t.GetClusterIndex2(nr);
1827 if ( (index>0) && (index&0x8000)==0){
1828 cl = t.GetClusterPointer(nr);
1829 if ( (cl==0) && (index>0)) cl = GetClusterMI(index);
1830 t.SetCurrentClusterIndex1(index);
1832 t.SetCurrentCluster(cl);
1838 // if (index<0) return 0;
1839 UInt_t uindex = TMath::Abs(index);
1842 //cl = krow.FindNearest2(y+10,z,roady,roadz,uindex);
1843 cl = krow.FindNearest2(y,z,roady,roadz,uindex);
1846 if (cl) t.SetCurrentClusterIndex1(krow.GetIndex(uindex));
1847 t.SetCurrentCluster(cl);
1853 Int_t AliTPCtrackerMI::FollowToNextCluster(AliTPCseed & t, Int_t nr) {
1854 //-----------------------------------------------------------------
1855 // This function tries to find a track prolongation to next pad row
1856 //-----------------------------------------------------------------
1858 //update error according neighborhoud
1860 if (t.GetCurrentCluster()) {
1862 Int_t accept = AcceptCluster(&t,t.GetCurrentCluster());
1864 if (t.GetCurrentCluster()->IsUsed(10)){
1869 t.SetNShared(t.GetNShared()+1);
1870 if (t.GetNShared()>0.7*t.GetNumberOfClusters()) {
1875 if (fIteration>0) accept = 0;
1876 if (accept<3) UpdateTrack(&t,accept);
1880 if ( ( (t.GetSigmaY2()+t.GetSigmaZ2())>0.16)&& t.GetNumberOfClusters()>18) t.SetRemoval(10);
1881 if ( t.GetChi2()/t.GetNumberOfClusters()>6 &&t.GetNumberOfClusters()>18) t.SetRemoval(10);
1883 if (( (t.GetNFoundable()*0.5 > t.GetNumberOfClusters()) || t.GetNoCluster()>15)) t.SetRemoval(10);
1891 //_____________________________________________________________________________
1892 Int_t AliTPCtrackerMI::FollowProlongation(AliTPCseed& t, Int_t rf, Int_t step) {
1893 //-----------------------------------------------------------------
1894 // This function tries to find a track prolongation.
1895 //-----------------------------------------------------------------
1896 Double_t xt=t.GetX();
1898 Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
1899 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
1900 if (alpha < 0. ) alpha += 2.*TMath::Pi();
1902 t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN);
1904 Int_t first = GetRowNumber(xt)-1;
1905 for (Int_t nr= first; nr>=rf; nr-=step) {
1907 if (t.GetKinkIndexes()[0]>0){
1908 for (Int_t i=0;i<3;i++){
1909 Int_t index = t.GetKinkIndexes()[i];
1910 if (index==0) break;
1911 if (index<0) continue;
1913 AliKink * kink = (AliKink*)fEvent->GetKink(index-1);
1915 printf("PROBLEM\n");
1918 Int_t kinkrow = kink->GetTPCRow0()+2+Int_t(0.5/(0.05+kink->GetAngle(2)));
1920 AliExternalTrackParam paramd(t);
1921 kink->SetDaughter(paramd);
1922 kink->SetStatus(2,5);
1929 if (nr==80) t.UpdateReference();
1930 if (nr<fInnerSec->GetNRows())
1931 fSectors = fInnerSec;
1933 fSectors = fOuterSec;
1934 if (FollowToNext(t,nr)==0)
1943 //_____________________________________________________________________________
1944 Int_t AliTPCtrackerMI::FollowProlongationFast(AliTPCseed& t, Int_t rf, Int_t step) {
1945 //-----------------------------------------------------------------
1946 // This function tries to find a track prolongation.
1947 //-----------------------------------------------------------------
1948 Double_t xt=t.GetX();
1950 Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
1951 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
1952 if (alpha < 0. ) alpha += 2.*TMath::Pi();
1953 t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN);
1955 for (Int_t nr=GetRowNumber(xt)-1; nr>=rf; nr-=step) {
1957 if (FollowToNextFast(t,nr)==0)
1958 if (!t.IsActive()) return 0;
1968 Int_t AliTPCtrackerMI::FollowBackProlongation(AliTPCseed& t, Int_t rf) {
1969 //-----------------------------------------------------------------
1970 // This function tries to find a track prolongation.
1971 //-----------------------------------------------------------------
1973 Double_t xt=t.GetX();
1974 Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
1975 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
1976 if (alpha < 0. ) alpha += 2.*TMath::Pi();
1977 t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN);
1979 Int_t first = t.GetFirstPoint();
1980 if (first<GetRowNumber(xt)+1) first = GetRowNumber(xt)+1;
1982 if (first<0) first=0;
1983 for (Int_t nr=first; nr<=rf; nr++) {
1984 // if ( (TMath::Abs(t.GetSnp())>0.95)) break;//patch 28 fev 06
1985 if (t.GetKinkIndexes()[0]<0){
1986 for (Int_t i=0;i<3;i++){
1987 Int_t index = t.GetKinkIndexes()[i];
1988 if (index==0) break;
1989 if (index>0) continue;
1990 index = TMath::Abs(index);
1991 AliKink * kink = (AliKink*)fEvent->GetKink(index-1);
1993 printf("PROBLEM\n");
1996 Int_t kinkrow = kink->GetTPCRow0()-2-Int_t(0.5/(0.05+kink->GetAngle(2)));
1998 AliExternalTrackParam paramm(t);
1999 kink->SetMother(paramm);
2000 kink->SetStatus(2,1);
2007 if (nr<fInnerSec->GetNRows())
2008 fSectors = fInnerSec;
2010 fSectors = fOuterSec;
2021 Float_t AliTPCtrackerMI::OverlapFactor(AliTPCseed * s1, AliTPCseed * s2, Int_t &sum1, Int_t & sum2)
2029 Float_t dz2 =(s1->GetZ() - s2->GetZ());
2032 Float_t dy2 =TMath::Abs((s1->GetY() - s2->GetY()));
2034 Float_t distance = TMath::Sqrt(dz2+dy2);
2035 if (distance>4.) return 0; // if there are far away - not overlap - to reduce combinatorics
2038 Int_t firstpoint = TMath::Min(s1->GetFirstPoint(),s2->GetFirstPoint());
2039 Int_t lastpoint = TMath::Max(s1->GetLastPoint(),s2->GetLastPoint());
2044 if (firstpoint>lastpoint) {
2045 firstpoint =lastpoint;
2050 for (Int_t i=firstpoint-1;i<lastpoint+1;i++){
2051 if (s1->GetClusterIndex2(i)>0) sum1++;
2052 if (s2->GetClusterIndex2(i)>0) sum2++;
2053 if (s1->GetClusterIndex2(i)==s2->GetClusterIndex2(i) && s1->GetClusterIndex2(i)>0) {
2057 if (sum<5) return 0;
2059 Float_t summin = TMath::Min(sum1+1,sum2+1);
2060 Float_t ratio = (sum+1)/Float_t(summin);
2064 void AliTPCtrackerMI::SignShared(AliTPCseed * s1, AliTPCseed * s2)
2068 Float_t thetaCut = 0.2;//+10.*TMath::Sqrt(s1->GetSigmaTglZ()+ s2->GetSigmaTglZ());
2069 if (TMath::Abs(s1->GetTgl()-s2->GetTgl())>thetaCut) return;
2070 Float_t minCl = TMath::Min(s1->GetNumberOfClusters(),s2->GetNumberOfClusters());
2071 Int_t cutN0 = TMath::Max(5,TMath::Nint(0.1*minCl));
2076 //Int_t firstpoint = TMath::Max(s1->GetFirstPoint(),s2->GetFirstPoint());
2077 //Int_t lastpoint = TMath::Min(s1->GetLastPoint(),s2->GetLastPoint());
2078 Int_t firstpoint = 0;
2079 Int_t lastpoint = 160;
2081 // if (firstpoint>=lastpoint-5) return;;
2083 for (Int_t i=firstpoint;i<lastpoint;i++){
2084 // if ( (s1->GetClusterIndex2(i)&0xFFFF8FFF)==(s2->GetClusterIndex2(i)&0xFFFF8FFF) && s1->GetClusterIndex2(i)>0) {
2085 if ( (s1->GetClusterIndex2(i))==(s2->GetClusterIndex2(i)) && s1->GetClusterIndex2(i)>0) {
2087 s1->SetSharedMapBit(i, kTRUE);
2088 s2->SetSharedMapBit(i, kTRUE);
2090 if (s1->GetClusterIndex2(i)>0)
2091 s1->SetClusterMapBit(i, kTRUE);
2093 if (sumshared>cutN0){
2096 for (Int_t i=firstpoint;i<lastpoint;i++){
2097 // if ( (s1->GetClusterIndex2(i)&0xFFFF8FFF)==(s2->GetClusterIndex2(i)&0xFFFF8FFF) && s1->GetClusterIndex2(i)>0) {
2098 if ( (s1->GetClusterIndex2(i))==(s2->GetClusterIndex2(i)) && s1->GetClusterIndex2(i)>0) {
2099 AliTPCTrackerPoint *p1 = s1->GetTrackPoint(i);
2100 AliTPCTrackerPoint *p2 = s2->GetTrackPoint(i);;
2101 if (s1->IsActive()&&s2->IsActive()){
2102 p1->SetShared(kTRUE);
2103 p2->SetShared(kTRUE);
2109 if (sumshared>cutN0){
2110 for (Int_t i=0;i<4;i++){
2111 if (s1->GetOverlapLabel(3*i)==0){
2112 s1->SetOverlapLabel(3*i, s2->GetLabel());
2113 s1->SetOverlapLabel(3*i+1,sumshared);
2114 s1->SetOverlapLabel(3*i+2,s2->GetUniqueID());
2118 for (Int_t i=0;i<4;i++){
2119 if (s2->GetOverlapLabel(3*i)==0){
2120 s2->SetOverlapLabel(3*i, s1->GetLabel());
2121 s2->SetOverlapLabel(3*i+1,sumshared);
2122 s2->SetOverlapLabel(3*i+2,s1->GetUniqueID());
2129 void AliTPCtrackerMI::SignShared(TObjArray * arr)
2132 //sort trackss according sectors
2134 for (Int_t i=0; i<arr->GetEntriesFast(); i++) {
2135 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2137 //if (pt) RotateToLocal(pt);
2141 arr->Sort(); // sorting according relative sectors
2142 arr->Expand(arr->GetEntries());
2145 Int_t nseed=arr->GetEntriesFast();
2146 for (Int_t i=0; i<nseed; i++) {
2147 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2149 for (Int_t j=0;j<=12;j++){
2150 pt->SetOverlapLabel(j,0);
2153 for (Int_t i=0; i<nseed; i++) {
2154 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2156 if (pt->GetRemoval()>10) continue;
2157 for (Int_t j=i+1; j<nseed; j++){
2158 AliTPCseed *pt2=(AliTPCseed*)arr->UncheckedAt(j);
2159 if (TMath::Abs(pt->GetRelativeSector()-pt2->GetRelativeSector())>1) continue;
2161 if (pt2->GetRemoval()<=10) {
2162 //if ( TMath::Abs(pt->GetRelativeSector()-pt2->GetRelativeSector())>0) break;
2170 void AliTPCtrackerMI::SortTracks(TObjArray * arr, Int_t mode) const
2173 //sort tracks in array according mode criteria
2174 Int_t nseed = arr->GetEntriesFast();
2175 for (Int_t i=0; i<nseed; i++) {
2176 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2187 void AliTPCtrackerMI::RemoveUsed2(TObjArray * arr, Float_t factor1, Float_t factor2, Int_t minimal)
2190 // Loop over all tracks and remove overlaped tracks (with lower quality)
2192 // 1. Unsign clusters
2193 // 2. Sort tracks according quality
2194 // Quality is defined by the number of cluster between first and last points
2196 // 3. Loop over tracks - decreasing quality order
2197 // a.) remove - If the fraction of shared cluster less than factor (1- n or 2)
2198 // b.) remove - If the minimal number of clusters less than minimal and not ITS
2199 // c.) if track accepted - sign clusters
2201 //Called in - AliTPCtrackerMI::Clusters2Tracks()
2202 // - AliTPCtrackerMI::PropagateBack()
2203 // - AliTPCtrackerMI::RefitInward()
2209 Int_t nseed = arr->GetEntriesFast();
2210 Float_t * quality = new Float_t[nseed];
2211 Int_t * indexes = new Int_t[nseed];
2215 for (Int_t i=0; i<nseed; i++) {
2216 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2221 pt->UpdatePoints(); //select first last max dens points
2222 Float_t * points = pt->GetPoints();
2223 if (points[3]<0.8) quality[i] =-1;
2224 quality[i] = (points[2]-points[0])+pt->GetNumberOfClusters();
2225 //prefer high momenta tracks if overlaps
2226 quality[i] *= TMath::Sqrt(TMath::Abs(pt->Pt())+0.5);
2228 TMath::Sort(nseed,quality,indexes);
2231 for (Int_t itrack=0; itrack<nseed; itrack++) {
2232 Int_t trackindex = indexes[itrack];
2233 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(trackindex);
2236 if (quality[trackindex]<0){
2238 delete arr->RemoveAt(trackindex);
2241 arr->RemoveAt(trackindex);
2247 Int_t first = Int_t(pt->GetPoints()[0]);
2248 Int_t last = Int_t(pt->GetPoints()[2]);
2249 Double_t factor = (pt->GetBConstrain()) ? factor1: factor2;
2251 Int_t found,foundable,shared;
2252 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
2253 // pt->GetClusterStatistic(0,160, found, foundable,shared,kFALSE);
2254 Bool_t itsgold =kFALSE;
2257 if (pt->GetESD()->GetITSclusters(dummy)>4) itsgold= kTRUE;
2261 if (Float_t(shared+1)/Float_t(found+1)>factor){
2262 if (pt->GetKinkIndexes()[0]!=0) continue; //don't remove tracks - part of the kinks
2263 delete arr->RemoveAt(trackindex);
2266 if (pt->GetNumberOfClusters()<50&&(found-0.5*shared)<minimal){ //remove short tracks
2267 if (pt->GetKinkIndexes()[0]!=0) continue; //don't remove tracks - part of the kinks
2268 delete arr->RemoveAt(trackindex);
2274 //if (sharedfactor>0.4) continue;
2275 if (pt->GetKinkIndexes()[0]>0) continue;
2276 //Remove tracks with undefined properties - seems
2277 if (pt->GetSigmaY2()<kAlmost0) continue; // ? what is the origin ?
2279 for (Int_t i=first; i<last; i++) {
2280 Int_t index=pt->GetClusterIndex2(i);
2281 // if (index<0 || index&0x8000 ) continue;
2282 if (index<0 || index&0x8000 ) continue;
2283 AliTPCclusterMI *c= pt->GetClusterPointer(i);
2290 Info("RemoveUsed","\n*****\nNumber of good tracks after shared removal\t%d\n",fNtracks);
2296 void AliTPCtrackerMI::UnsignClusters()
2299 // loop over all clusters and unsign them
2302 for (Int_t sec=0;sec<fkNIS;sec++){
2303 for (Int_t row=0;row<fInnerSec->GetNRows();row++){
2304 AliTPCclusterMI *cl = fInnerSec[sec][row].GetClusters1();
2305 for (Int_t icl =0;icl< fInnerSec[sec][row].GetN1();icl++)
2306 // if (cl[icl].IsUsed(10))
2308 cl = fInnerSec[sec][row].GetClusters2();
2309 for (Int_t icl =0;icl< fInnerSec[sec][row].GetN2();icl++)
2310 //if (cl[icl].IsUsed(10))
2315 for (Int_t sec=0;sec<fkNOS;sec++){
2316 for (Int_t row=0;row<fOuterSec->GetNRows();row++){
2317 AliTPCclusterMI *cl = fOuterSec[sec][row].GetClusters1();
2318 for (Int_t icl =0;icl< fOuterSec[sec][row].GetN1();icl++)
2319 //if (cl[icl].IsUsed(10))
2321 cl = fOuterSec[sec][row].GetClusters2();
2322 for (Int_t icl =0;icl< fOuterSec[sec][row].GetN2();icl++)
2323 //if (cl[icl].IsUsed(10))
2332 void AliTPCtrackerMI::SignClusters(TObjArray * arr, Float_t fnumber, Float_t fdensity)
2335 //sign clusters to be "used"
2337 // snumber and sdensity sign number of sigmas - bellow mean value to be accepted
2338 // loop over "primaries"
2352 Int_t nseed = arr->GetEntriesFast();
2353 for (Int_t i=0; i<nseed; i++) {
2354 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2358 if (!(pt->IsActive())) continue;
2359 Float_t dens = pt->GetNumberOfClusters()/Float_t(pt->GetNFoundable());
2360 if ( (dens>0.7) && (pt->GetNumberOfClusters()>70)){
2362 sumdens2+= dens*dens;
2363 sumn += pt->GetNumberOfClusters();
2364 sumn2 += pt->GetNumberOfClusters()*pt->GetNumberOfClusters();
2365 Float_t chi2 = pt->GetChi2()/pt->GetNumberOfClusters();
2368 sumchi2 +=chi2*chi2;
2373 Float_t mdensity = 0.9;
2374 Float_t meann = 130;
2375 Float_t meanchi = 1;
2376 Float_t sdensity = 0.1;
2377 Float_t smeann = 10;
2378 Float_t smeanchi =0.4;
2382 mdensity = sumdens/sum;
2384 meanchi = sumchi/sum;
2386 sdensity = sumdens2/sum-mdensity*mdensity;
2388 sdensity = TMath::Sqrt(sdensity);
2392 smeann = sumn2/sum-meann*meann;
2394 smeann = TMath::Sqrt(smeann);
2398 smeanchi = sumchi2/sum - meanchi*meanchi;
2400 smeanchi = TMath::Sqrt(smeanchi);
2406 //REMOVE SHORT DELTAS or tracks going out of sensitive volume of TPC
2408 for (Int_t i=0; i<nseed; i++) {
2409 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2413 if (pt->GetBSigned()) continue;
2414 if (pt->GetBConstrain()) continue;
2415 //if (!(pt->IsActive())) continue;
2417 Int_t found,foundable,shared;
2418 pt->GetClusterStatistic(0,160,found, foundable,shared);
2419 if (shared/float(found)>0.3) {
2420 if (shared/float(found)>0.9 ){
2421 //delete arr->RemoveAt(i);
2426 Bool_t isok =kFALSE;
2427 if ( (pt->GetNShared()/pt->GetNumberOfClusters()<0.5) &&pt->GetNumberOfClusters()>60)
2429 if ((TMath::Abs(1/pt->GetC())<100.) && (pt->GetNShared()/pt->GetNumberOfClusters()<0.7))
2431 if (TMath::Abs(pt->GetZ()/pt->GetX())>1.1)
2433 if ( (TMath::Abs(pt->GetSnp()>0.7) && pt->GetD(0,0)>60.))
2437 for (Int_t i=0; i<160; i++) {
2438 Int_t index=pt->GetClusterIndex2(i);
2439 if (index<0) continue;
2440 AliTPCclusterMI *c= pt->GetClusterPointer(i);
2442 //if (!(c->IsUsed(10))) c->Use();
2449 Double_t maxchi = meanchi+2.*smeanchi;
2451 for (Int_t i=0; i<nseed; i++) {
2452 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2456 //if (!(pt->IsActive())) continue;
2457 if (pt->GetBSigned()) continue;
2458 Double_t chi = pt->GetChi2()/pt->GetNumberOfClusters();
2459 if (chi>maxchi) continue;
2462 Float_t dens = pt->GetNumberOfClusters()/Float_t(pt->GetNFoundable());
2464 //sign only tracks with enoug big density at the beginning
2466 if ((pt->GetDensityFirst(40)<0.75) && pt->GetNumberOfClusters()<meann) continue;
2469 Double_t mindens = TMath::Max(double(mdensity-sdensity*fdensity*bfactor),0.65);
2470 Double_t minn = TMath::Max(Int_t(meann-fnumber*smeann*bfactor),50);
2472 // if (pt->fBConstrain) mindens = TMath::Max(mdensity-sdensity*fdensity*bfactor,0.65);
2473 if ( (pt->GetRemoval()==10) && (pt->GetSnp()>0.8)&&(dens>mindens))
2476 if ((dens>mindens && pt->GetNumberOfClusters()>minn) && chi<maxchi ){
2477 //Int_t noc=pt->GetNumberOfClusters();
2478 pt->SetBSigned(kTRUE);
2479 for (Int_t i=0; i<160; i++) {
2481 Int_t index=pt->GetClusterIndex2(i);
2482 if (index<0) continue;
2483 AliTPCclusterMI *c= pt->GetClusterPointer(i);
2485 // if (!(c->IsUsed(10))) c->Use();
2490 // gLastCheck = nseed;
2498 void AliTPCtrackerMI::StopNotActive(TObjArray * arr, Int_t row0, Float_t th0, Float_t th1, Float_t th2) const
2500 // stop not active tracks
2501 // take th1 as threshold for number of founded to number of foundable on last 10 active rows
2502 // take th2 as threshold for number of founded to number of foundable on last 20 active rows
2503 Int_t nseed = arr->GetEntriesFast();
2505 for (Int_t i=0; i<nseed; i++) {
2506 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2510 if (!(pt->IsActive())) continue;
2511 StopNotActive(pt,row0,th0, th1,th2);
2517 void AliTPCtrackerMI::StopNotActive(AliTPCseed * seed, Int_t row0, Float_t th0, Float_t th1,
2520 // stop not active tracks
2521 // take th1 as threshold for number of founded to number of foundable on last 10 active rows
2522 // take th2 as threshold for number of founded to number of foundable on last 20 active rows
2525 Int_t foundable = 0;
2526 Int_t maxindex = seed->GetLastPoint(); //last foundable row
2527 if (seed->GetNFoundable()*th0 > seed->GetNumberOfClusters()) {
2528 seed->Desactivate(10) ;
2532 for (Int_t i=row0; i<maxindex; i++){
2533 Int_t index = seed->GetClusterIndex2(i);
2534 if (index!=-1) foundable++;
2536 if (foundable<=30) sumgood1++;
2537 if (foundable<=50) {
2544 if (foundable>=30.){
2545 if (sumgood1<(th1*30.)) seed->Desactivate(10);
2548 if (sumgood2<(th2*50.)) seed->Desactivate(10);
2552 Int_t AliTPCtrackerMI::RefitInward(AliESDEvent *event)
2555 // back propagation of ESD tracks
2558 const Int_t kMaxFriendTracks=2000;
2562 //PrepareForProlongation(fSeeds,1);
2563 PropagateForward2(fSeeds);
2564 RemoveUsed2(fSeeds,0.4,0.4,20);
2566 TObjArray arraySeed(fSeeds->GetEntries());
2567 for (Int_t i=0;i<fSeeds->GetEntries();i++) {
2568 arraySeed.AddAt(fSeeds->At(i),i);
2570 SignShared(&arraySeed);
2571 // FindCurling(fSeeds, event,2); // find multi found tracks
2572 FindSplitted(fSeeds, event,2); // find multi found tracks
2575 Int_t nseed = fSeeds->GetEntriesFast();
2576 for (Int_t i=0;i<nseed;i++){
2577 AliTPCseed * seed = (AliTPCseed*) fSeeds->UncheckedAt(i);
2578 if (!seed) continue;
2579 if (seed->GetKinkIndex(0)>0) UpdateKinkQualityD(seed); // update quality informations for kinks
2581 seed->PropagateTo(fParam->GetInnerRadiusLow());
2582 seed->UpdatePoints();
2583 AddCovariance(seed);
2585 AliESDtrack *esd=event->GetTrack(i);
2586 seed->CookdEdx(0.02,0.6);
2587 CookLabel(seed,0.1); //For comparison only
2588 esd->SetTPCClusterMap(seed->GetClusterMap());
2589 esd->SetTPCSharedMap(seed->GetSharedMap());
2591 if (AliTPCReconstructor::StreamLevel()>1 && seed!=0&&esd!=0) {
2592 TTreeSRedirector &cstream = *fDebugStreamer;
2599 if (seed->GetNumberOfClusters()>15){
2600 esd->UpdateTrackParams(seed,AliESDtrack::kTPCrefit);
2601 esd->SetTPCPoints(seed->GetPoints());
2602 esd->SetTPCPointsF(seed->GetNFoundable());
2603 Int_t ndedx = seed->GetNCDEDX(0)+seed->GetNCDEDX(1)+seed->GetNCDEDX(2)+seed->GetNCDEDX(3);
2604 Float_t sdedx = (seed->GetSDEDX(0)+seed->GetSDEDX(1)+seed->GetSDEDX(2)+seed->GetSDEDX(3))*0.25;
2605 Float_t dedx = seed->GetdEdx();
2606 esd->SetTPCsignal(dedx, sdedx, ndedx);
2608 // add seed to the esd track in Calib level
2610 Bool_t storeFriend = gRandom->Rndm()<(kMaxFriendTracks)/Float_t(nseed);
2611 if (AliTPCReconstructor::StreamLevel()>0 &&storeFriend){
2612 AliTPCseed * seedCopy = new AliTPCseed(*seed, kTRUE);
2613 esd->AddCalibObject(seedCopy);
2618 //printf("problem\n");
2621 //FindKinks(fSeeds,event);
2622 Info("RefitInward","Number of refitted tracks %d",ntracks);
2627 Int_t AliTPCtrackerMI::PropagateBack(AliESDEvent *event)
2630 // back propagation of ESD tracks
2636 PropagateBack(fSeeds);
2637 RemoveUsed2(fSeeds,0.4,0.4,20);
2638 //FindCurling(fSeeds, fEvent,1);
2639 FindSplitted(fSeeds, event,1); // find multi found tracks
2642 Int_t nseed = fSeeds->GetEntriesFast();
2644 for (Int_t i=0;i<nseed;i++){
2645 AliTPCseed * seed = (AliTPCseed*) fSeeds->UncheckedAt(i);
2646 if (!seed) continue;
2647 if (seed->GetKinkIndex(0)<0) UpdateKinkQualityM(seed); // update quality informations for kinks
2648 seed->UpdatePoints();
2649 AddCovariance(seed);
2650 AliESDtrack *esd=event->GetTrack(i);
2651 seed->CookdEdx(0.02,0.6);
2652 CookLabel(seed,0.1); //For comparison only
2653 if (seed->GetNumberOfClusters()>15){
2654 esd->UpdateTrackParams(seed,AliESDtrack::kTPCout);
2655 esd->SetTPCPoints(seed->GetPoints());
2656 esd->SetTPCPointsF(seed->GetNFoundable());
2657 Int_t ndedx = seed->GetNCDEDX(0)+seed->GetNCDEDX(1)+seed->GetNCDEDX(2)+seed->GetNCDEDX(3);
2658 Float_t sdedx = (seed->GetSDEDX(0)+seed->GetSDEDX(1)+seed->GetSDEDX(2)+seed->GetSDEDX(3))*0.25;
2659 Float_t dedx = seed->GetdEdx();
2660 esd->SetTPCsignal(dedx, sdedx, ndedx);
2662 Int_t eventnumber = event->GetEventNumberInFile();// patch 28 fev 06
2663 // This is most likely NOT the event number you'd like to use. It has nothing to do with the 'real' event number
2664 if (AliTPCReconstructor::StreamLevel()>1 && esd) {
2665 (*fDebugStreamer)<<"Cback"<<
2668 "EventNrInFile="<<eventnumber<<
2673 //FindKinks(fSeeds,event);
2674 Info("PropagateBack","Number of back propagated tracks %d",ntracks);
2681 void AliTPCtrackerMI::DeleteSeeds()
2686 Int_t nseed = fSeeds->GetEntriesFast();
2687 for (Int_t i=0;i<nseed;i++){
2688 AliTPCseed * seed = (AliTPCseed*)fSeeds->At(i);
2689 if (seed) delete fSeeds->RemoveAt(i);
2696 void AliTPCtrackerMI::ReadSeeds(AliESDEvent *event, Int_t direction)
2699 //read seeds from the event
2701 Int_t nentr=event->GetNumberOfTracks();
2703 Info("ReadSeeds", "Number of ESD tracks: %d\n", nentr);
2708 fSeeds = new TObjArray(nentr);
2712 for (Int_t i=0; i<nentr; i++) {
2713 AliESDtrack *esd=event->GetTrack(i);
2714 ULong_t status=esd->GetStatus();
2715 if (!(status&AliESDtrack::kTPCin)) continue;
2716 AliTPCtrack t(*esd);
2717 t.SetNumberOfClusters(0);
2718 // AliTPCseed *seed = new AliTPCseed(t,t.GetAlpha());
2719 AliTPCseed *seed = new AliTPCseed(t/*,t.GetAlpha()*/);
2720 seed->SetUniqueID(esd->GetID());
2721 AddCovariance(seed); //add systematic ucertainty
2722 for (Int_t ikink=0;ikink<3;ikink++) {
2723 Int_t index = esd->GetKinkIndex(ikink);
2724 seed->GetKinkIndexes()[ikink] = index;
2725 if (index==0) continue;
2726 index = TMath::Abs(index);
2727 AliESDkink * kink = fEvent->GetKink(index-1);
2728 if (kink&&esd->GetKinkIndex(ikink)<0){
2729 if ((status & AliESDtrack::kTRDrefit) != 0) kink->SetStatus(1,2);
2730 if ((status & AliESDtrack::kITSout) != 0) kink->SetStatus(1,0);
2732 if (kink&&esd->GetKinkIndex(ikink)>0){
2733 if ((status & AliESDtrack::kTRDrefit) != 0) kink->SetStatus(1,6);
2734 if ((status & AliESDtrack::kITSout) != 0) kink->SetStatus(1,4);
2738 if (((status&AliESDtrack::kITSout)==0)&&(direction==1)) seed->ResetCovariance(10.);
2739 if ( direction ==2 &&(status & AliESDtrack::kTRDrefit) == 0 ) seed->ResetCovariance(10.);
2740 if ( direction ==2 && ((status & AliESDtrack::kTPCout) == 0) ) {
2745 if ( direction ==2 &&(status & AliESDtrack::kTRDrefit) > 0 ) {
2746 Double_t par0[5],par1[5],alpha,x;
2747 esd->GetInnerExternalParameters(alpha,x,par0);
2748 esd->GetExternalParameters(x,par1);
2749 Double_t delta1 = TMath::Abs(par0[4]-par1[4])/(0.000000001+TMath::Abs(par0[4]+par1[4]));
2750 Double_t delta2 = TMath::Abs(par0[3]-par1[3]);
2752 if (esd->GetTRDncls()>0) trdchi2 = esd->GetTRDchi2()/esd->GetTRDncls();
2753 //reset covariance if suspicious
2754 if ( (delta1>0.1) || (delta2>0.006) ||trdchi2>7.)
2755 seed->ResetCovariance(10.);
2760 // rotate to the local coordinate system
2762 fSectors=fInnerSec; fN=fkNIS;
2763 Double_t alpha=seed->GetAlpha() - fSectors->GetAlphaShift();
2764 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
2765 if (alpha < 0. ) alpha += 2.*TMath::Pi();
2766 Int_t ns=Int_t(alpha/fSectors->GetAlpha())%fN;
2767 alpha =ns*fSectors->GetAlpha() + fSectors->GetAlphaShift();
2768 if (alpha<-TMath::Pi()) alpha += 2*TMath::Pi();
2769 if (alpha>=TMath::Pi()) alpha -= 2*TMath::Pi();
2770 alpha-=seed->GetAlpha();
2771 if (!seed->Rotate(alpha)) {
2777 if (esd->GetKinkIndex(0)<=0){
2778 for (Int_t irow=0;irow<160;irow++){
2779 Int_t index = seed->GetClusterIndex2(irow);
2782 AliTPCclusterMI * cl = GetClusterMI(index);
2783 seed->SetClusterPointer(irow,cl);
2785 if ((index & 0x8000)==0){
2786 cl->Use(10); // accepted cluster
2788 cl->Use(6); // close cluster not accepted
2791 Info("ReadSeeds","Not found cluster");
2796 fSeeds->AddAt(seed,i);
2802 //_____________________________________________________________________________
2803 void AliTPCtrackerMI::MakeSeeds3(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, Float_t cuts[4],
2804 Float_t deltay, Int_t ddsec) {
2805 //-----------------------------------------------------------------
2806 // This function creates track seeds.
2807 // SEEDING WITH VERTEX CONSTRAIN
2808 //-----------------------------------------------------------------
2809 // cuts[0] - fP4 cut
2810 // cuts[1] - tan(phi) cut
2811 // cuts[2] - zvertex cut
2812 // cuts[3] - fP3 cut
2820 Double_t x[5], c[15];
2821 // Int_t di = i1-i2;
2823 AliTPCseed * seed = new AliTPCseed();
2824 Double_t alpha=fSectors->GetAlpha(), shift=fSectors->GetAlphaShift();
2825 Double_t cs=cos(alpha), sn=sin(alpha);
2827 // Double_t x1 =fOuterSec->GetX(i1);
2828 //Double_t xx2=fOuterSec->GetX(i2);
2830 Double_t x1 =GetXrow(i1);
2831 Double_t xx2=GetXrow(i2);
2833 Double_t x3=GetX(), y3=GetY(), z3=GetZ();
2835 Int_t imiddle = (i2+i1)/2; //middle pad row index
2836 Double_t xm = GetXrow(imiddle); // radius of middle pad-row
2837 const AliTPCtrackerRow& krm=GetRow(sec,imiddle); //middle pad -row
2841 const AliTPCtrackerRow& kr1=GetRow(ns,i1);
2842 Double_t ymax = GetMaxY(i1)-kr1.GetDeadZone()-1.5;
2843 Double_t ymaxm = GetMaxY(imiddle)-kr1.GetDeadZone()-1.5;
2846 // change cut on curvature if it can't reach this layer
2847 // maximal curvature set to reach it
2848 Double_t dvertexmax = TMath::Sqrt((x1-x3)*(x1-x3)+(ymax+5-y3)*(ymax+5-y3));
2849 if (dvertexmax*0.5*cuts[0]>0.85){
2850 cuts[0] = 0.85/(dvertexmax*0.5+1.);
2852 Double_t r2min = 1/(cuts[0]*cuts[0]); //minimal square of radius given by cut
2855 if (deltay>0) ddsec = 0;
2856 // loop over clusters
2857 for (Int_t is=0; is < kr1; is++) {
2859 if (kr1[is]->IsUsed(10)) continue;
2860 Double_t y1=kr1[is]->GetY(), z1=kr1[is]->GetZ();
2861 //if (TMath::Abs(y1)>ymax) continue;
2863 if (deltay>0 && TMath::Abs(ymax-TMath::Abs(y1))> deltay ) continue; // seed only at the edge
2865 // find possible directions
2866 Float_t anglez = (z1-z3)/(x1-x3);
2867 Float_t extraz = z1 - anglez*(x1-xx2); // extrapolated z
2870 //find rotation angles relative to line given by vertex and point 1
2871 Double_t dvertex2 = (x1-x3)*(x1-x3)+(y1-y3)*(y1-y3);
2872 Double_t dvertex = TMath::Sqrt(dvertex2);
2873 Double_t angle13 = TMath::ATan((y1-y3)/(x1-x3));
2874 Double_t cs13 = cos(-angle13), sn13 = sin(-angle13);
2877 // loop over 2 sectors
2883 Double_t dddz1=0; // direction of delta inclination in z axis
2890 for (Int_t dsec = dsec1; dsec<=dsec2;dsec++){
2891 Int_t sec2 = sec + dsec;
2893 // AliTPCtrackerRow& kr2 = fOuterSec[(sec2+fkNOS)%fkNOS][i2];
2894 //AliTPCtrackerRow& kr2m = fOuterSec[(sec2+fkNOS)%fkNOS][imiddle];
2895 AliTPCtrackerRow& kr2 = GetRow((sec2+fkNOS)%fkNOS,i2);
2896 AliTPCtrackerRow& kr2m = GetRow((sec2+fkNOS)%fkNOS,imiddle);
2897 Int_t index1 = TMath::Max(kr2.Find(extraz-0.6-dddz1*TMath::Abs(z1)*0.05)-1,0);
2898 Int_t index2 = TMath::Min(kr2.Find(extraz+0.6+dddz2*TMath::Abs(z1)*0.05)+1,kr2);
2900 // rotation angles to p1-p3
2901 Double_t cs13r = cos(-angle13+dsec*alpha)/dvertex, sn13r = sin(-angle13+dsec*alpha)/dvertex;
2902 Double_t x2, y2, z2;
2904 // Double_t dymax = maxangle*TMath::Abs(x1-xx2);
2907 Double_t dxx0 = (xx2-x3)*cs13r;
2908 Double_t dyy0 = (xx2-x3)*sn13r;
2909 for (Int_t js=index1; js < index2; js++) {
2910 const AliTPCclusterMI *kcl = kr2[js];
2911 if (kcl->IsUsed(10)) continue;
2913 //calcutate parameters
2915 Double_t yy0 = dyy0 +(kcl->GetY()-y3)*cs13r;
2917 if (TMath::Abs(yy0)<0.000001) continue;
2918 Double_t xx0 = dxx0 -(kcl->GetY()-y3)*sn13r;
2919 Double_t y0 = 0.5*(xx0*xx0+yy0*yy0-xx0)/yy0;
2920 Double_t r02 = (0.25+y0*y0)*dvertex2;
2921 //curvature (radius) cut
2922 if (r02<r2min) continue;
2926 Double_t c0 = 1/TMath::Sqrt(r02);
2930 //Double_t dfi0 = 2.*TMath::ASin(dvertex*c0*0.5);
2931 //Double_t dfi1 = 2.*TMath::ASin(TMath::Sqrt(yy0*yy0+(1-xx0)*(1-xx0))*dvertex*c0*0.5);
2932 Double_t dfi0 = 2.*AliTPCFastMath::FastAsin(dvertex*c0*0.5);
2933 Double_t dfi1 = 2.*AliTPCFastMath::FastAsin(TMath::Sqrt(yy0*yy0+(1-xx0)*(1-xx0))*dvertex*c0*0.5);
2936 Double_t z0 = kcl->GetZ();
2937 Double_t zzzz2 = z1-(z1-z3)*dfi1/dfi0;
2938 if (TMath::Abs(zzzz2-z0)>0.5) continue;
2941 Double_t dip = (z1-z0)*c0/dfi1;
2942 Double_t x0 = (0.5*cs13+y0*sn13)*dvertex*c0;
2953 x2= xx2*cs-y2*sn*dsec;
2954 y2=+xx2*sn*dsec+y2*cs;
2964 // do we have cluster at the middle ?
2966 GetProlongation(x1,xm,x,ym,zm);
2968 AliTPCclusterMI * cm=0;
2969 if (TMath::Abs(ym)-ymaxm<0){
2970 cm = krm.FindNearest2(ym,zm,1.0,0.6,dummy);
2971 if ((!cm) || (cm->IsUsed(10))) {
2976 // rotate y1 to system 0
2977 // get state vector in rotated system
2978 Double_t yr1 = (-0.5*sn13+y0*cs13)*dvertex*c0;
2979 Double_t xr2 = x0*cs+yr1*sn*dsec;
2980 Double_t xr[5]={kcl->GetY(),kcl->GetZ(), xr2, dip, c0};
2982 GetProlongation(xx2,xm,xr,ym,zm);
2983 if (TMath::Abs(ym)-ymaxm<0){
2984 cm = kr2m.FindNearest2(ym,zm,1.0,0.6,dummy);
2985 if ((!cm) || (cm->IsUsed(10))) {
2995 dym = ym - cm->GetY();
2996 dzm = zm - cm->GetZ();
3003 Double_t sy1=kr1[is]->GetSigmaY2()*2., sz1=kr1[is]->GetSigmaZ2()*2.;
3004 Double_t sy2=kcl->GetSigmaY2()*2., sz2=kcl->GetSigmaZ2()*2.;
3005 //Double_t sy3=400*3./12., sy=0.1, sz=0.1;
3006 Double_t sy3=25000*x[4]*x[4]+0.1, sy=0.1, sz=0.1;
3007 //Double_t sy3=25000*x[4]*x[4]*60+0.5, sy=0.1, sz=0.1;
3009 Double_t f40=(F1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
3010 Double_t f42=(F1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
3011 Double_t f43=(F1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
3012 Double_t f20=(F2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
3013 Double_t f22=(F2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
3014 Double_t f23=(F2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
3016 Double_t f30=(F3(x1,y1+sy,x2,y2,z1,z2)-x[3])/sy;
3017 Double_t f31=(F3(x1,y1,x2,y2,z1+sz,z2)-x[3])/sz;
3018 Double_t f32=(F3(x1,y1,x2,y2+sy,z1,z2)-x[3])/sy;
3019 Double_t f34=(F3(x1,y1,x2,y2,z1,z2+sz)-x[3])/sz;
3023 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3024 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3025 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3026 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3027 c[13]=f30*sy1*f40+f32*sy2*f42;
3028 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3030 // if (!BuildSeed(kr1[is],kcl,0,x1,x2,x3,x,c)) continue;
3032 UInt_t index=kr1.GetIndex(is);
3033 seed->~AliTPCseed(); // this does not set the pointer to 0...
3034 AliTPCseed *track=new(seed) AliTPCseed(x1, ns*alpha+shift, x, c, index);
3036 track->SetIsSeeding(kTRUE);
3037 track->SetSeed1(i1);
3038 track->SetSeed2(i2);
3039 track->SetSeedType(3);
3043 FollowProlongation(*track, (i1+i2)/2,1);
3044 Int_t foundable,found,shared;
3045 track->GetClusterStatistic((i1+i2)/2,i1, found, foundable, shared, kTRUE);
3046 if ((found<0.55*foundable) || shared>0.5*found || (track->GetSigmaY2()+track->GetSigmaZ2())>0.5){
3048 seed->~AliTPCseed();
3054 FollowProlongation(*track, i2,1);
3058 track->SetBConstrain(1);
3059 // track->fLastPoint = i1+fInnerSec->GetNRows(); // first cluster in track position
3060 track->SetLastPoint(i1); // first cluster in track position
3061 track->SetFirstPoint(track->GetLastPoint());
3063 if (track->GetNumberOfClusters()<(i1-i2)*0.5 ||
3064 track->GetNumberOfClusters() < track->GetNFoundable()*0.6 ||
3065 track->GetNShared()>0.4*track->GetNumberOfClusters() ) {
3067 seed->~AliTPCseed();
3071 // Z VERTEX CONDITION
3072 Double_t zv, bz=GetBz();
3073 if ( !track->GetZAt(0.,bz,zv) ) continue;
3074 if (TMath::Abs(zv-z3)>cuts[2]) {
3075 FollowProlongation(*track, TMath::Max(i2-20,0));
3076 if ( !track->GetZAt(0.,bz,zv) ) continue;
3077 if (TMath::Abs(zv-z3)>cuts[2]){
3078 FollowProlongation(*track, TMath::Max(i2-40,0));
3079 if ( !track->GetZAt(0.,bz,zv) ) continue;
3080 if (TMath::Abs(zv-z3)>cuts[2] &&(track->GetNumberOfClusters() > track->GetNFoundable()*0.7)){
3081 // make seed without constrain
3082 AliTPCseed * track2 = MakeSeed(track,0.2,0.5,1.);
3083 FollowProlongation(*track2, i2,1);
3084 track2->SetBConstrain(kFALSE);
3085 track2->SetSeedType(1);
3086 arr->AddLast(track2);
3088 seed->~AliTPCseed();
3093 seed->~AliTPCseed();
3100 track->SetSeedType(0);
3101 arr->AddLast(track);
3102 seed = new AliTPCseed;
3104 // don't consider other combinations
3105 if (track->GetNumberOfClusters() > track->GetNFoundable()*0.8)
3111 Info("MakeSeeds3","\nSeeding statistic:\t%d\t%d\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin,nout1,nout2);
3117 void AliTPCtrackerMI::MakeSeeds5(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, Float_t cuts[4],
3122 //-----------------------------------------------------------------
3123 // This function creates track seeds.
3124 //-----------------------------------------------------------------
3125 // cuts[0] - fP4 cut
3126 // cuts[1] - tan(phi) cut
3127 // cuts[2] - zvertex cut
3128 // cuts[3] - fP3 cut
3138 Double_t x[5], c[15];
3140 // make temporary seed
3141 AliTPCseed * seed = new AliTPCseed;
3142 Double_t alpha=fOuterSec->GetAlpha(), shift=fOuterSec->GetAlphaShift();
3143 // Double_t cs=cos(alpha), sn=sin(alpha);
3148 Double_t x1 = GetXrow(i1-1);
3149 const AliTPCtrackerRow& kr1=GetRow(sec,i1-1);
3150 Double_t y1max = GetMaxY(i1-1)-kr1.GetDeadZone()-1.5;
3152 Double_t x1p = GetXrow(i1);
3153 const AliTPCtrackerRow& kr1p=GetRow(sec,i1);
3155 Double_t x1m = GetXrow(i1-2);
3156 const AliTPCtrackerRow& kr1m=GetRow(sec,i1-2);
3159 //last 3 padrow for seeding
3160 AliTPCtrackerRow& kr3 = GetRow((sec+fkNOS)%fkNOS,i1-7);
3161 Double_t x3 = GetXrow(i1-7);
3162 // Double_t y3max= GetMaxY(i1-7)-kr3.fDeadZone-1.5;
3164 AliTPCtrackerRow& kr3p = GetRow((sec+fkNOS)%fkNOS,i1-6);
3165 Double_t x3p = GetXrow(i1-6);
3167 AliTPCtrackerRow& kr3m = GetRow((sec+fkNOS)%fkNOS,i1-8);
3168 Double_t x3m = GetXrow(i1-8);
3173 Int_t im = i1-4; //middle pad row index
3174 Double_t xm = GetXrow(im); // radius of middle pad-row
3175 const AliTPCtrackerRow& krm=GetRow(sec,im); //middle pad -row
3176 // Double_t ymmax = GetMaxY(im)-kr1.fDeadZone-1.5;
3179 Double_t deltax = x1-x3;
3180 Double_t dymax = deltax*cuts[1];
3181 Double_t dzmax = deltax*cuts[3];
3183 // loop over clusters
3184 for (Int_t is=0; is < kr1; is++) {
3186 if (kr1[is]->IsUsed(10)) continue;
3187 Double_t y1=kr1[is]->GetY(), z1=kr1[is]->GetZ();
3189 if (deltay>0 && TMath::Abs(y1max-TMath::Abs(y1))> deltay ) continue; // seed only at the edge
3191 Int_t index1 = TMath::Max(kr3.Find(z1-dzmax)-1,0);
3192 Int_t index2 = TMath::Min(kr3.Find(z1+dzmax)+1,kr3);
3198 for (Int_t js=index1; js < index2; js++) {
3199 const AliTPCclusterMI *kcl = kr3[js];
3200 if (kcl->IsUsed(10)) continue;
3202 // apply angular cuts
3203 if (TMath::Abs(y1-y3)>dymax) continue;
3206 if (TMath::Abs(z1-z3)>dzmax) continue;
3208 Double_t angley = (y1-y3)/(x1-x3);
3209 Double_t anglez = (z1-z3)/(x1-x3);
3211 Double_t erry = TMath::Abs(angley)*(x1-x1m)*0.5+0.5;
3212 Double_t errz = TMath::Abs(anglez)*(x1-x1m)*0.5+0.5;
3214 Double_t yyym = angley*(xm-x1)+y1;
3215 Double_t zzzm = anglez*(xm-x1)+z1;
3217 const AliTPCclusterMI *kcm = krm.FindNearest2(yyym,zzzm,erry,errz,index);
3219 if (kcm->IsUsed(10)) continue;
3221 erry = TMath::Abs(angley)*(x1-x1m)*0.4+0.5;
3222 errz = TMath::Abs(anglez)*(x1-x1m)*0.4+0.5;
3229 // look around first
3230 const AliTPCclusterMI *kc1m = kr1m.FindNearest2(angley*(x1m-x1)+y1,
3236 if (kc1m->IsUsed(10)) used++;
3238 const AliTPCclusterMI *kc1p = kr1p.FindNearest2(angley*(x1p-x1)+y1,
3244 if (kc1p->IsUsed(10)) used++;
3246 if (used>1) continue;
3247 if (found<1) continue;
3251 const AliTPCclusterMI *kc3m = kr3m.FindNearest2(angley*(x3m-x3)+y3,
3257 if (kc3m->IsUsed(10)) used++;
3261 const AliTPCclusterMI *kc3p = kr3p.FindNearest2(angley*(x3p-x3)+y3,
3267 if (kc3p->IsUsed(10)) used++;
3271 if (used>1) continue;
3272 if (found<3) continue;
3282 x[4]=F1(x1,y1,x2,y2,x3,y3);
3283 //if (TMath::Abs(x[4]) >= cuts[0]) continue;
3286 x[2]=F2(x1,y1,x2,y2,x3,y3);
3289 x[3]=F3n(x1,y1,x2,y2,z1,z2,x[4]);
3290 //if (TMath::Abs(x[3]) > cuts[3]) continue;
3294 Double_t sy1=0.1, sz1=0.1;
3295 Double_t sy2=0.1, sz2=0.1;
3296 Double_t sy3=0.1, sy=0.1, sz=0.1;
3298 Double_t f40=(F1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
3299 Double_t f42=(F1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
3300 Double_t f43=(F1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
3301 Double_t f20=(F2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
3302 Double_t f22=(F2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
3303 Double_t f23=(F2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
3305 Double_t f30=(F3(x1,y1+sy,x2,y2,z1,z2)-x[3])/sy;
3306 Double_t f31=(F3(x1,y1,x2,y2,z1+sz,z2)-x[3])/sz;
3307 Double_t f32=(F3(x1,y1,x2,y2+sy,z1,z2)-x[3])/sy;
3308 Double_t f34=(F3(x1,y1,x2,y2,z1,z2+sz)-x[3])/sz;
3312 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3313 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3314 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3315 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3316 c[13]=f30*sy1*f40+f32*sy2*f42;
3317 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3319 // if (!BuildSeed(kr1[is],kcl,0,x1,x2,x3,x,c)) continue;
3321 UInt_t index=kr1.GetIndex(is);
3322 seed->~AliTPCseed();
3323 AliTPCseed *track=new(seed) AliTPCseed(x1, sec*alpha+shift, x, c, index);
3325 track->SetIsSeeding(kTRUE);
3328 FollowProlongation(*track, i1-7,1);
3329 if (track->GetNumberOfClusters() < track->GetNFoundable()*0.75 ||
3330 track->GetNShared()>0.6*track->GetNumberOfClusters() || ( track->GetSigmaY2()+ track->GetSigmaZ2())>0.6){
3332 seed->~AliTPCseed();
3338 FollowProlongation(*track, i2,1);
3339 track->SetBConstrain(0);
3340 track->SetLastPoint(i1+fInnerSec->GetNRows()); // first cluster in track position
3341 track->SetFirstPoint(track->GetLastPoint());
3343 if (track->GetNumberOfClusters()<(i1-i2)*0.5 ||
3344 track->GetNumberOfClusters()<track->GetNFoundable()*0.7 ||
3345 track->GetNShared()>2. || track->GetChi2()/track->GetNumberOfClusters()>6 || ( track->GetSigmaY2()+ track->GetSigmaZ2())>0.5 ) {
3347 seed->~AliTPCseed();
3352 FollowProlongation(*track, TMath::Max(i2-10,0),1);
3353 AliTPCseed * track2 = MakeSeed(track,0.2,0.5,0.9);
3354 FollowProlongation(*track2, i2,1);
3355 track2->SetBConstrain(kFALSE);
3356 track2->SetSeedType(4);
3357 arr->AddLast(track2);
3359 seed->~AliTPCseed();
3363 //arr->AddLast(track);
3364 //seed = new AliTPCseed;
3370 Info("MakeSeeds5","\nSeeding statiistic:\t%d\t%d\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin,nout1,nout2,nout3);
3376 //_____________________________________________________________________________
3377 void AliTPCtrackerMI::MakeSeeds2(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, Float_t */*cuts[4]*/,
3378 Float_t deltay, Bool_t /*bconstrain*/) {
3379 //-----------------------------------------------------------------
3380 // This function creates track seeds - without vertex constraint
3381 //-----------------------------------------------------------------
3382 // cuts[0] - fP4 cut - not applied
3383 // cuts[1] - tan(phi) cut
3384 // cuts[2] - zvertex cut - not applied
3385 // cuts[3] - fP3 cut
3395 Double_t alpha=fOuterSec->GetAlpha(), shift=fOuterSec->GetAlphaShift();
3396 // Double_t cs=cos(alpha), sn=sin(alpha);
3397 Int_t row0 = (i1+i2)/2;
3398 Int_t drow = (i1-i2)/2;
3399 const AliTPCtrackerRow& kr0=fSectors[sec][row0];
3400 AliTPCtrackerRow * kr=0;
3402 AliTPCpolyTrack polytrack;
3403 Int_t nclusters=fSectors[sec][row0];
3404 AliTPCseed * seed = new AliTPCseed;
3409 for (Int_t is=0; is < nclusters; is++) { //LOOP over clusters
3411 Int_t nfoundable =0;
3412 for (Int_t iter =1; iter<2; iter++){ //iterations
3413 const AliTPCtrackerRow& krm=fSectors[sec][row0-iter];
3414 const AliTPCtrackerRow& krp=fSectors[sec][row0+iter];
3415 const AliTPCclusterMI * cl= kr0[is];
3417 if (cl->IsUsed(10)) {
3423 Double_t x = kr0.GetX();
3424 // Initialization of the polytrack
3429 Double_t y0= cl->GetY();
3430 Double_t z0= cl->GetZ();
3434 Double_t ymax = fSectors->GetMaxY(row0)-kr0.GetDeadZone()-1.5;
3435 if (deltay>0 && TMath::Abs(ymax-TMath::Abs(y0))> deltay ) continue; // seed only at the edge
3437 erry = (0.5)*cl->GetSigmaY2()/TMath::Sqrt(cl->GetQ())*6;
3438 errz = (0.5)*cl->GetSigmaZ2()/TMath::Sqrt(cl->GetQ())*6;
3439 polytrack.AddPoint(x,y0,z0,erry, errz);
3442 if (cl->IsUsed(10)) sumused++;
3445 Float_t roady = (5*TMath::Sqrt(cl->GetSigmaY2()+0.2)+1.)*iter;
3446 Float_t roadz = (5*TMath::Sqrt(cl->GetSigmaZ2()+0.2)+1.)*iter;
3449 AliTPCclusterMI * cl1 = krm.FindNearest(y0,z0,roady,roadz);
3450 if (cl1 && TMath::Abs(ymax-TMath::Abs(y0))) {
3451 erry = (0.5)*cl1->GetSigmaY2()/TMath::Sqrt(cl1->GetQ())*3;
3452 errz = (0.5)*cl1->GetSigmaZ2()/TMath::Sqrt(cl1->GetQ())*3;
3453 if (cl1->IsUsed(10)) sumused++;
3454 polytrack.AddPoint(x,cl1->GetY(),cl1->GetZ(),erry,errz);
3458 AliTPCclusterMI * cl2 = krp.FindNearest(y0,z0,roady,roadz);
3460 erry = (0.5)*cl2->GetSigmaY2()/TMath::Sqrt(cl2->GetQ())*3;
3461 errz = (0.5)*cl2->GetSigmaZ2()/TMath::Sqrt(cl2->GetQ())*3;
3462 if (cl2->IsUsed(10)) sumused++;
3463 polytrack.AddPoint(x,cl2->GetY(),cl2->GetZ(),erry,errz);
3466 if (sumused>0) continue;
3468 polytrack.UpdateParameters();
3474 nfoundable = polytrack.GetN();
3475 nfound = nfoundable;
3477 for (Int_t ddrow = iter+1; ddrow<drow;ddrow++){
3478 Float_t maxdist = 0.8*(1.+3./(ddrow));
3479 for (Int_t delta = -1;delta<=1;delta+=2){
3480 Int_t row = row0+ddrow*delta;
3481 kr = &(fSectors[sec][row]);
3482 Double_t xn = kr->GetX();
3483 Double_t ymax = fSectors->GetMaxY(row)-kr->GetDeadZone()-1.5;
3484 polytrack.GetFitPoint(xn,yn,zn);
3485 if (TMath::Abs(yn)>ymax) continue;
3487 AliTPCclusterMI * cln = kr->FindNearest(yn,zn,roady,roadz);
3489 Float_t dist = TMath::Sqrt( (yn-cln->GetY())*(yn-cln->GetY())+(zn-cln->GetZ())*(zn-cln->GetZ()));
3492 erry = (dist+0.3)*cln->GetSigmaY2()/TMath::Sqrt(cln->GetQ())*(1.+1./(ddrow));
3493 errz = (dist+0.3)*cln->GetSigmaZ2()/TMath::Sqrt(cln->GetQ())*(1.+1./(ddrow));
3494 if (cln->IsUsed(10)) {
3495 // printf("used\n");
3503 polytrack.AddPoint(xn,cln->GetY(),cln->GetZ(),erry, errz);
3508 if ( (sumused>3) || (sumused>0.5*nfound) || (nfound<0.6*nfoundable)) break;
3509 polytrack.UpdateParameters();
3512 if ( (sumused>3) || (sumused>0.5*nfound)) {
3513 //printf("sumused %d\n",sumused);
3518 polytrack.GetFitDerivation(kr0.GetX(),dy,dz);
3519 AliTPCpolyTrack track2;
3521 polytrack.Refit(track2,0.5+TMath::Abs(dy)*0.3,0.4+TMath::Abs(dz)*0.3);
3522 if (track2.GetN()<0.5*nfoundable) continue;
3525 if ((nfound>0.6*nfoundable) &&( nfoundable>0.4*(i1-i2))) {
3527 // test seed with and without constrain
3528 for (Int_t constrain=0; constrain<=0;constrain++){
3529 // add polytrack candidate
3531 Double_t x[5], c[15];
3532 Double_t x1,x2,x3,y1,y2,y3,z1,z2,z3;
3533 track2.GetBoundaries(x3,x1);
3535 track2.GetFitPoint(x1,y1,z1);
3536 track2.GetFitPoint(x2,y2,z2);
3537 track2.GetFitPoint(x3,y3,z3);
3539 //is track pointing to the vertex ?
3542 polytrack.GetFitPoint(x0,y0,z0);
3555 x[4]=F1(x1,y1,x2,y2,x3,y3);
3557 // if (TMath::Abs(x[4]) >= cuts[0]) continue; //
3558 x[2]=F2(x1,y1,x2,y2,x3,y3);
3560 //if (TMath::Abs(x[4]*x1-x[2]) >= cuts[1]) continue;
3561 //x[3]=F3(x1,y1,x2,y2,z1,z2);
3562 x[3]=F3n(x1,y1,x3,y3,z1,z3,x[4]);
3563 //if (TMath::Abs(x[3]) > cuts[3]) continue;
3566 Double_t sy =0.1, sz =0.1;
3567 Double_t sy1=0.02, sz1=0.02;
3568 Double_t sy2=0.02, sz2=0.02;
3572 sy3=25000*x[4]*x[4]+0.1, sy=0.1, sz=0.1;
3575 Double_t f40=(F1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
3576 Double_t f42=(F1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
3577 Double_t f43=(F1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
3578 Double_t f20=(F2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
3579 Double_t f22=(F2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
3580 Double_t f23=(F2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
3582 Double_t f30=(F3(x1,y1+sy,x3,y3,z1,z3)-x[3])/sy;
3583 Double_t f31=(F3(x1,y1,x3,y3,z1+sz,z3)-x[3])/sz;
3584 Double_t f32=(F3(x1,y1,x3,y3+sy,z1,z3)-x[3])/sy;
3585 Double_t f34=(F3(x1,y1,x3,y3,z1,z3+sz)-x[3])/sz;
3590 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3591 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3592 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3593 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3594 c[13]=f30*sy1*f40+f32*sy2*f42;
3595 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3597 //Int_t row1 = fSectors->GetRowNumber(x1);
3598 Int_t row1 = GetRowNumber(x1);
3602 seed->~AliTPCseed();
3603 AliTPCseed *track=new(seed) AliTPCseed(x1,sec*alpha+shift,x,c,index);
3604 track->SetIsSeeding(kTRUE);
3605 Int_t rc=FollowProlongation(*track, i2);
3606 if (constrain) track->SetBConstrain(1);
3608 track->SetBConstrain(0);
3609 track->SetLastPoint(row1+fInnerSec->GetNRows()); // first cluster in track position
3610 track->SetFirstPoint(track->GetLastPoint());
3612 if (rc==0 || track->GetNumberOfClusters()<(i1-i2)*0.5 ||
3613 track->GetNumberOfClusters() < track->GetNFoundable()*0.6 ||
3614 track->GetNShared()>0.4*track->GetNumberOfClusters()) {
3617 seed->~AliTPCseed();
3620 arr->AddLast(track);
3621 seed = new AliTPCseed;
3625 } // if accepted seed
3628 Info("MakeSeeds2","\nSeeding statiistic:\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin3);
3634 AliTPCseed *AliTPCtrackerMI::MakeSeed(AliTPCseed *track, Float_t r0, Float_t r1, Float_t r2)
3638 //reseed using track points
3639 Int_t p0 = int(r0*track->GetNumberOfClusters()); // point 0
3640 Int_t p1 = int(r1*track->GetNumberOfClusters());
3641 Int_t p2 = int(r2*track->GetNumberOfClusters()); // last point
3643 Double_t x0[3],x1[3],x2[3];
3644 for (Int_t i=0;i<3;i++){
3650 // find track position at given ratio of the length
3651 Int_t sec0=0, sec1=0, sec2=0;
3654 for (Int_t i=0;i<160;i++){
3655 if (track->GetClusterPointer(i)){
3657 AliTPCTrackerPoint *trpoint =track->GetTrackPoint(i);
3658 if ( (index<p0) || x0[0]<0 ){
3659 if (trpoint->GetX()>1){
3660 clindex = track->GetClusterIndex2(i);
3662 x0[0] = trpoint->GetX();
3663 x0[1] = trpoint->GetY();
3664 x0[2] = trpoint->GetZ();
3665 sec0 = ((clindex&0xff000000)>>24)%18;
3670 if ( (index<p1) &&(trpoint->GetX()>1)){
3671 clindex = track->GetClusterIndex2(i);
3673 x1[0] = trpoint->GetX();
3674 x1[1] = trpoint->GetY();
3675 x1[2] = trpoint->GetZ();
3676 sec1 = ((clindex&0xff000000)>>24)%18;
3679 if ( (index<p2) &&(trpoint->GetX()>1)){
3680 clindex = track->GetClusterIndex2(i);
3682 x2[0] = trpoint->GetX();
3683 x2[1] = trpoint->GetY();
3684 x2[2] = trpoint->GetZ();
3685 sec2 = ((clindex&0xff000000)>>24)%18;
3692 Double_t alpha, cs,sn, xx2,yy2;
3694 alpha = (sec1-sec2)*fSectors->GetAlpha();
3695 cs = TMath::Cos(alpha);
3696 sn = TMath::Sin(alpha);
3697 xx2= x1[0]*cs-x1[1]*sn;
3698 yy2= x1[0]*sn+x1[1]*cs;
3702 alpha = (sec0-sec2)*fSectors->GetAlpha();
3703 cs = TMath::Cos(alpha);
3704 sn = TMath::Sin(alpha);
3705 xx2= x0[0]*cs-x0[1]*sn;
3706 yy2= x0[0]*sn+x0[1]*cs;
3712 Double_t x[5],c[15];
3716 x[4]=F1(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]);
3717 // if (x[4]>1) return 0;
3718 x[2]=F2(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]);
3719 x[3]=F3n(x2[0],x2[1],x0[0],x0[1],x2[2],x0[2],x[4]);
3720 //if (TMath::Abs(x[3]) > 2.2) return 0;
3721 //if (TMath::Abs(x[2]) > 1.99) return 0;
3723 Double_t sy =0.1, sz =0.1;
3725 Double_t sy1=0.02+track->GetSigmaY2(), sz1=0.02+track->GetSigmaZ2();
3726 Double_t sy2=0.01+track->GetSigmaY2(), sz2=0.01+track->GetSigmaZ2();
3727 Double_t sy3=0.01+track->GetSigmaY2();
3729 Double_t f40=(F1(x2[0],x2[1]+sy,x1[0],x1[1],x0[0],x0[1])-x[4])/sy;
3730 Double_t f42=(F1(x2[0],x2[1],x1[0],x1[1]+sy,x0[0],x0[1])-x[4])/sy;
3731 Double_t f43=(F1(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]+sy)-x[4])/sy;
3732 Double_t f20=(F2(x2[0],x2[1]+sy,x1[0],x1[1],x0[0],x0[1])-x[2])/sy;
3733 Double_t f22=(F2(x2[0],x2[1],x1[0],x1[1]+sy,x0[0],x0[1])-x[2])/sy;
3734 Double_t f23=(F2(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]+sy)-x[2])/sy;
3736 Double_t f30=(F3(x2[0],x2[1]+sy,x0[0],x0[1],x2[2],x0[2])-x[3])/sy;
3737 Double_t f31=(F3(x2[0],x2[1],x0[0],x0[1],x2[2]+sz,x0[2])-x[3])/sz;
3738 Double_t f32=(F3(x2[0],x2[1],x0[0],x0[1]+sy,x2[2],x0[2])-x[3])/sy;
3739 Double_t f34=(F3(x2[0],x2[1],x0[0],x0[1],x2[2],x0[2]+sz)-x[3])/sz;
3744 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3745 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3746 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3747 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3748 c[13]=f30*sy1*f40+f32*sy2*f42;
3749 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3751 // Int_t row1 = fSectors->GetRowNumber(x2[0]);
3752 AliTPCseed *seed=new AliTPCseed(x2[0], sec2*fSectors->GetAlpha()+fSectors->GetAlphaShift(), x, c, 0);
3753 // Double_t y0,z0,y1,z1, y2,z2;
3754 //seed->GetProlongation(x0[0],y0,z0);
3755 // seed->GetProlongation(x1[0],y1,z1);
3756 //seed->GetProlongation(x2[0],y2,z2);
3758 seed->SetLastPoint(pp2);
3759 seed->SetFirstPoint(pp2);
3766 AliTPCseed *AliTPCtrackerMI::ReSeed(AliTPCseed *track, Float_t r0, Float_t r1, Float_t r2)
3770 //reseed using founded clusters
3772 // Find the number of clusters
3773 Int_t nclusters = 0;
3774 for (Int_t irow=0;irow<160;irow++){
3775 if (track->GetClusterIndex(irow)>0) nclusters++;
3779 ipos[0] = TMath::Max(int(r0*nclusters),0); // point 0 cluster
3780 ipos[1] = TMath::Min(int(r1*nclusters),nclusters-1); //
3781 ipos[2] = TMath::Min(int(r2*nclusters),nclusters-1); // last point
3785 Int_t row[3],sec[3]={0,0,0};
3787 // find track row position at given ratio of the length
3789 for (Int_t irow=0;irow<160;irow++){
3790 if (track->GetClusterIndex2(irow)<0) continue;
3792 for (Int_t ipoint=0;ipoint<3;ipoint++){
3793 if (index<=ipos[ipoint]) row[ipoint] = irow;
3797 //Get cluster and sector position
3798 for (Int_t ipoint=0;ipoint<3;ipoint++){
3799 Int_t clindex = track->GetClusterIndex2(row[ipoint]);
3800 AliTPCclusterMI * cl = GetClusterMI(clindex);
3803 // AliTPCclusterMI * cl = GetClusterMI(clindex);
3806 sec[ipoint] = ((clindex&0xff000000)>>24)%18;
3807 xyz[ipoint][0] = GetXrow(row[ipoint]);
3808 xyz[ipoint][1] = cl->GetY();
3809 xyz[ipoint][2] = cl->GetZ();
3813 // Calculate seed state vector and covariance matrix
3815 Double_t alpha, cs,sn, xx2,yy2;
3817 alpha = (sec[1]-sec[2])*fSectors->GetAlpha();
3818 cs = TMath::Cos(alpha);
3819 sn = TMath::Sin(alpha);
3820 xx2= xyz[1][0]*cs-xyz[1][1]*sn;
3821 yy2= xyz[1][0]*sn+xyz[1][1]*cs;
3825 alpha = (sec[0]-sec[2])*fSectors->GetAlpha();
3826 cs = TMath::Cos(alpha);
3827 sn = TMath::Sin(alpha);
3828 xx2= xyz[0][0]*cs-xyz[0][1]*sn;
3829 yy2= xyz[0][0]*sn+xyz[0][1]*cs;
3835 Double_t x[5],c[15];
3839 x[4]=F1(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
3840 x[2]=F2(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
3841 x[3]=F3n(xyz[2][0],xyz[2][1],xyz[0][0],xyz[0][1],xyz[2][2],xyz[0][2],x[4]);
3843 Double_t sy =0.1, sz =0.1;
3845 Double_t sy1=0.2, sz1=0.2;
3846 Double_t sy2=0.2, sz2=0.2;
3849 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;
3850 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;
3851 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;
3852 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;
3853 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;
3854 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;
3856 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;
3857 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;
3858 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;
3859 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;
3864 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3865 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3866 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3867 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3868 c[13]=f30*sy1*f40+f32*sy2*f42;
3869 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3871 // Int_t row1 = fSectors->GetRowNumber(xyz[2][0]);
3872 AliTPCseed *seed=new AliTPCseed(xyz[2][0], sec[2]*fSectors->GetAlpha()+fSectors->GetAlphaShift(), x, c, 0);
3873 seed->SetLastPoint(row[2]);
3874 seed->SetFirstPoint(row[2]);
3879 AliTPCseed *AliTPCtrackerMI::ReSeed(AliTPCseed *track,Int_t r0, Bool_t forward)
3883 //reseed using founded clusters
3886 Int_t row[3]={0,0,0};
3887 Int_t sec[3]={0,0,0};
3889 // forward direction
3891 for (Int_t irow=r0;irow<160;irow++){
3892 if (track->GetClusterIndex(irow)>0){
3897 for (Int_t irow=160;irow>r0;irow--){
3898 if (track->GetClusterIndex(irow)>0){
3903 for (Int_t irow=row[2]-15;irow>row[0];irow--){
3904 if (track->GetClusterIndex(irow)>0){
3912 for (Int_t irow=0;irow<r0;irow++){
3913 if (track->GetClusterIndex(irow)>0){
3918 for (Int_t irow=r0;irow>0;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){
3932 if ((row[2]-row[0])<20) return 0;
3933 if (row[1]==0) return 0;
3936 //Get cluster and sector position
3937 for (Int_t ipoint=0;ipoint<3;ipoint++){
3938 Int_t clindex = track->GetClusterIndex2(row[ipoint]);
3939 AliTPCclusterMI * cl = GetClusterMI(clindex);
3942 // AliTPCclusterMI * cl = GetClusterMI(clindex);
3945 sec[ipoint] = ((clindex&0xff000000)>>24)%18;
3946 xyz[ipoint][0] = GetXrow(row[ipoint]);
3947 AliTPCTrackerPoint * point = track->GetTrackPoint(row[ipoint]);
3948 if (point&&ipoint<2){
3950 xyz[ipoint][1] = point->GetY();
3951 xyz[ipoint][2] = point->GetZ();
3954 xyz[ipoint][1] = cl->GetY();
3955 xyz[ipoint][2] = cl->GetZ();
3962 // Calculate seed state vector and covariance matrix
3964 Double_t alpha, cs,sn, xx2,yy2;
3966 alpha = (sec[1]-sec[2])*fSectors->GetAlpha();
3967 cs = TMath::Cos(alpha);
3968 sn = TMath::Sin(alpha);
3969 xx2= xyz[1][0]*cs-xyz[1][1]*sn;
3970 yy2= xyz[1][0]*sn+xyz[1][1]*cs;
3974 alpha = (sec[0]-sec[2])*fSectors->GetAlpha();
3975 cs = TMath::Cos(alpha);
3976 sn = TMath::Sin(alpha);
3977 xx2= xyz[0][0]*cs-xyz[0][1]*sn;
3978 yy2= xyz[0][0]*sn+xyz[0][1]*cs;
3984 Double_t x[5],c[15];
3988 x[4]=F1(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
3989 x[2]=F2(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
3990 x[3]=F3n(xyz[2][0],xyz[2][1],xyz[0][0],xyz[0][1],xyz[2][2],xyz[0][2],x[4]);
3992 Double_t sy =0.1, sz =0.1;
3994 Double_t sy1=0.2, sz1=0.2;
3995 Double_t sy2=0.2, sz2=0.2;
3998 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;
3999 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;
4000 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;
4001 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;
4002 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;
4003 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;
4005 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;
4006 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;
4007 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;
4008 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;
4013 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
4014 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
4015 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
4016 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
4017 c[13]=f30*sy1*f40+f32*sy2*f42;
4018 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
4020 // Int_t row1 = fSectors->GetRowNumber(xyz[2][0]);
4021 AliTPCseed *seed=new AliTPCseed(xyz[2][0], sec[2]*fSectors->GetAlpha()+fSectors->GetAlphaShift(), x, c, 0);
4022 seed->SetLastPoint(row[2]);
4023 seed->SetFirstPoint(row[2]);
4024 for (Int_t i=row[0];i<row[2];i++){
4025 seed->SetClusterIndex(i, track->GetClusterIndex(i));
4033 void AliTPCtrackerMI::FindMultiMC(TObjArray * array, AliESDEvent */*esd*/, Int_t iter)
4036 // find multi tracks - THIS FUNCTION IS ONLY FOR DEBUG PURPOSES
4038 // Use AliTPCReconstructor::StreamLevel()>2 if you want to tune parameters - cuts
4040 // Two reasons to have multiple find tracks
4041 // 1. Curling tracks can be find more than once
4042 // 2. Splitted tracks
4043 // a.) Multiple seeding to increase tracking efficiency - (~ 100% reached)
4044 // b.) Edge effect on the sector boundaries
4047 // Algorithm done in 2 phases - because of CPU consumption
4048 // it is n^2 algorithm - for lead-lead 20000x20000 combination are investigated
4050 // Algorihm for curling tracks sign:
4051 // 1 phase -makes a very rough fast cuts to minimize combinatorics
4052 // a.) opposite sign
4053 // b.) one of the tracks - not pointing to the primary vertex -
4054 // c.) delta tan(theta)
4056 // 2 phase - calculates DCA between tracks - time consument
4061 // General cuts - for splitted tracks and for curling tracks
4063 const Float_t kMaxdPhi = 0.2; // maximal distance in phi
4065 // Curling tracks cuts
4070 Int_t nentries = array->GetEntriesFast();
4071 AliHelix *helixes = new AliHelix[nentries];
4072 Float_t *xm = new Float_t[nentries];
4073 Float_t *dz0 = new Float_t[nentries];
4074 Float_t *dz1 = new Float_t[nentries];
4080 // Find track COG in x direction - point with best defined parameters
4082 for (Int_t i=0;i<nentries;i++){
4083 AliTPCseed* track = (AliTPCseed*)array->At(i);
4084 if (!track) continue;
4085 track->SetCircular(0);
4086 new (&helixes[i]) AliHelix(*track);
4090 track->GetDZ(GetX(),GetY(),GetZ(),GetBz(),dz);
4093 for (Int_t icl=0; icl<160; icl++){
4094 AliTPCclusterMI * cl = track->GetClusterPointer(icl);
4100 if (ncl>0) xm[i]/=Float_t(ncl);
4102 TTreeSRedirector &cstream = *fDebugStreamer;
4104 for (Int_t i0=0;i0<nentries;i0++){
4105 AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
4106 if (!track0) continue;
4107 Float_t xc0 = helixes[i0].GetHelix(6);
4108 Float_t yc0 = helixes[i0].GetHelix(7);
4109 Float_t r0 = helixes[i0].GetHelix(8);
4110 Float_t rc0 = TMath::Sqrt(xc0*xc0+yc0*yc0);
4111 Float_t fi0 = TMath::ATan2(yc0,xc0);
4113 for (Int_t i1=i0+1;i1<nentries;i1++){
4114 AliTPCseed * track1 = (AliTPCseed*)array->At(i1);
4115 if (!track1) continue;
4116 Int_t lab0=track0->GetLabel();
4117 Int_t lab1=track1->GetLabel();
4118 if (TMath::Abs(lab0)!=TMath::Abs(lab1)) continue;
4120 Float_t xc1 = helixes[i1].GetHelix(6);
4121 Float_t yc1 = helixes[i1].GetHelix(7);
4122 Float_t r1 = helixes[i1].GetHelix(8);
4123 Float_t rc1 = TMath::Sqrt(xc1*xc1+yc1*yc1);
4124 Float_t fi1 = TMath::ATan2(yc1,xc1);
4126 Float_t dfi = fi0-fi1;
4129 if (dfi>1.5*TMath::Pi()) dfi-=TMath::Pi(); // take care about edge effect
4130 if (dfi<-1.5*TMath::Pi()) dfi+=TMath::Pi(); //
4131 if (TMath::Abs(dfi)>kMaxdPhi&&helixes[i0].GetHelix(4)*helixes[i1].GetHelix(4)<0){
4133 // if short tracks with undefined sign
4134 fi1 = -TMath::ATan2(yc1,-xc1);
4137 Float_t dtheta = TMath::Abs(track0->GetTgl()-track1->GetTgl())<TMath::Abs(track0->GetTgl()+track1->GetTgl())? track0->GetTgl()-track1->GetTgl():track0->GetTgl()+track1->GetTgl();
4140 // debug stream to tune "fast cuts"
4142 Double_t dist[3]; // distance at X
4143 Double_t mdist[3]={0,0,0}; // mean distance X+-40cm
4144 track0->GetDistance(track1,0.5*(xm[i0]+xm[i1])-40.,dist,AliTracker::GetBz());
4145 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[i]);
4146 track0->GetDistance(track1,0.5*(xm[i0]+xm[i1])+40.,dist,AliTracker::GetBz());
4147 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[i]);
4148 track0->GetDistance(track1,0.5*(xm[i0]+xm[i1]),dist,AliTracker::GetBz());
4149 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[i]);
4150 for (Int_t i=0;i<3;i++) mdist[i]*=0.33333;
4154 for (Int_t icl=0; icl<160; icl++){
4155 AliTPCclusterMI * cl0 = track0->GetClusterPointer(icl);
4156 AliTPCclusterMI * cl1 = track1->GetClusterPointer(icl);
4159 if (cl0==cl1) sums++;
4167 "Tr0.="<<track0<< // seed0
4168 "Tr1.="<<track1<< // seed1
4169 "h0.="<<&helixes[i0]<<
4170 "h1.="<<&helixes[i1]<<
4172 "sum="<<sum<< //the sum of rows with cl in both
4173 "sums="<<sums<< //the sum of shared clusters
4174 "xm0="<<xm[i0]<< // the center of track
4175 "xm1="<<xm[i1]<< // the x center of track
4176 // General cut variables
4177 "dfi="<<dfi<< // distance in fi angle
4178 "dtheta="<<dtheta<< // distance int theta angle
4184 "dist0="<<dist[0]<< //distance x
4185 "dist1="<<dist[1]<< //distance y
4186 "dist2="<<dist[2]<< //distance z
4187 "mdist0="<<mdist[0]<< //distance x
4188 "mdist1="<<mdist[1]<< //distance y
4189 "mdist2="<<mdist[2]<< //distance z
4202 if (AliTPCReconstructor::StreamLevel()>1) {
4203 AliInfo("Time for curling tracks removal DEBUGGING MC");
4209 void AliTPCtrackerMI::FindSplitted(TObjArray * array, AliESDEvent */*esd*/, Int_t iter)
4213 // Two reasons to have multiple find tracks
4214 // 1. Curling tracks can be find more than once
4215 // 2. Splitted tracks
4216 // a.) Multiple seeding to increase tracking efficiency - (~ 100% reached)
4217 // b.) Edge effect on the sector boundaries
4219 // This function tries to find tracks closed in the parametric space
4221 // cut logic if distance is bigger than cut continue - Do Nothing
4222 const Float_t kMaxdTheta = 0.05; // maximal distance in theta
4223 const Float_t kMaxdPhi = 0.05; // maximal deistance in phi
4224 const Float_t kdelta = 40.; //delta r to calculate track distance
4226 // const Float_t kMaxDist0 = 1.; // maximal distance 0
4227 //const Float_t kMaxDist1 = 0.3; // maximal distance 1 - cut if track in separate rows
4230 TCut csec("csec","abs(Tr0.fRelativeSector-Tr1.fRelativeSector)<2");
4231 TCut cdtheta("cdtheta","abs(dtheta)<0.05");
4236 Int_t nentries = array->GetEntriesFast();
4237 AliHelix *helixes = new AliHelix[nentries];
4238 Float_t *xm = new Float_t[nentries];
4244 //Sort tracks according quality
4246 Int_t nseed = array->GetEntriesFast();
4247 Float_t * quality = new Float_t[nseed];
4248 Int_t * indexes = new Int_t[nseed];
4249 for (Int_t i=0; i<nseed; i++) {
4250 AliTPCseed *pt=(AliTPCseed*)array->UncheckedAt(i);
4255 pt->UpdatePoints(); //select first last max dens points
4256 Float_t * points = pt->GetPoints();
4257 if (points[3]<0.8) quality[i] =-1;
4258 quality[i] = (points[2]-points[0])+pt->GetNumberOfClusters();
4259 //prefer high momenta tracks if overlaps
4260 quality[i] *= TMath::Sqrt(TMath::Abs(pt->Pt())+0.5);
4262 TMath::Sort(nseed,quality,indexes);
4266 // Find track COG in x direction - point with best defined parameters
4268 for (Int_t i=0;i<nentries;i++){
4269 AliTPCseed* track = (AliTPCseed*)array->At(i);
4270 if (!track) continue;
4271 track->SetCircular(0);
4272 new (&helixes[i]) AliHelix(*track);
4275 for (Int_t icl=0; icl<160; icl++){
4276 AliTPCclusterMI * cl = track->GetClusterPointer(icl);
4282 if (ncl>0) xm[i]/=Float_t(ncl);
4284 TTreeSRedirector &cstream = *fDebugStreamer;
4286 for (Int_t is0=0;is0<nentries;is0++){
4287 Int_t i0 = indexes[is0];
4288 AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
4289 if (!track0) continue;
4290 if (track0->GetKinkIndexes()[0]!=0) continue;
4291 Float_t xc0 = helixes[i0].GetHelix(6);
4292 Float_t yc0 = helixes[i0].GetHelix(7);
4293 Float_t fi0 = TMath::ATan2(yc0,xc0);
4295 for (Int_t is1=is0+1;is1<nentries;is1++){
4296 Int_t i1 = indexes[is1];
4297 AliTPCseed * track1 = (AliTPCseed*)array->At(i1);
4298 if (!track1) continue;
4300 if (TMath::Abs(track0->GetRelativeSector()-track1->GetRelativeSector())>1) continue;
4301 if (track1->GetKinkIndexes()[0]>0 &&track0->GetKinkIndexes()[0]<0) continue;
4302 if (track1->GetKinkIndexes()[0]!=0) continue;
4304 Float_t dtheta = TMath::Abs(track0->GetTgl()-track1->GetTgl())<TMath::Abs(track0->GetTgl()+track1->GetTgl())? track0->GetTgl()-track1->GetTgl():track0->GetTgl()+track1->GetTgl();
4305 if (TMath::Abs(dtheta)>kMaxdTheta) continue;
4307 Float_t xc1 = helixes[i1].GetHelix(6);
4308 Float_t yc1 = helixes[i1].GetHelix(7);
4309 Float_t fi1 = TMath::ATan2(yc1,xc1);
4311 Float_t dfi = fi0-fi1;
4312 if (dfi>1.5*TMath::Pi()) dfi-=TMath::Pi(); // take care about edge effect
4313 if (dfi<-1.5*TMath::Pi()) dfi+=TMath::Pi(); //
4314 if (TMath::Abs(dfi)>kMaxdPhi&&helixes[i0].GetHelix(4)*helixes[i1].GetHelix(4)<0){
4316 // if short tracks with undefined sign
4317 fi1 = -TMath::ATan2(yc1,-xc1);
4320 if (TMath::Abs(dfi)>kMaxdPhi) continue;
4327 for (Int_t icl=0; icl<160; icl++){
4328 Int_t index0=track0->GetClusterIndex2(icl);
4329 Int_t index1=track1->GetClusterIndex2(icl);
4330 Bool_t used0 = (index0>0 && !(index0&0x8000));
4331 Bool_t used1 = (index1>0 && !(index1&0x8000));
4333 if (used0) sum0++; // used cluster0
4334 if (used1) sum1++; // used clusters1
4335 if (used0&&used1) sum++;
4336 if (index0==index1 && used0 && used1) sums++;
4340 if (sums<10) continue;
4341 if (sum<40) continue;
4342 if (sums/Float_t(TMath::Min(sum0,sum1))<0.5) continue;
4344 Double_t dist[5][4]; // distance at X
4345 Double_t mdist[4]={0,0,0,0}; // mean distance on range +-delta
4349 track0->GetDistance(track1,xm[i0],dist[0],AliTracker::GetBz());
4350 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[0][i]);
4351 track0->GetDistance(track1,xm[i1],dist[1],AliTracker::GetBz());
4352 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[1][i]);
4354 track0->GetDistance(track1,TMath::Min(xm[i1],xm[i0])-kdelta,dist[2],AliTracker::GetBz());
4355 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[2][i]);
4356 track0->GetDistance(track1,TMath::Max(xm[i1],xm[i0])+kdelta,dist[3],AliTracker::GetBz());
4357 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[3][i]);
4359 track0->GetDistance(track1,(xm[i1]+xm[i0])*0.5,dist[4],AliTracker::GetBz());
4360 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[4][i]);
4361 for (Int_t i=0;i<3;i++) mdist[i]*=0.2;
4364 Int_t lab0=track0->GetLabel();
4365 Int_t lab1=track1->GetLabel();
4366 cstream<<"Splitted"<<
4370 "Tr0.="<<track0<< // seed0
4371 "Tr1.="<<track1<< // seed1
4372 "h0.="<<&helixes[i0]<<
4373 "h1.="<<&helixes[i1]<<
4375 "sum="<<sum<< //the sum of rows with cl in both
4376 "sum0="<<sum0<< //the sum of rows with cl in 0 track
4377 "sum1="<<sum1<< //the sum of rows with cl in 1 track
4378 "sums="<<sums<< //the sum of shared clusters
4379 "xm0="<<xm[i0]<< // the center of track
4380 "xm1="<<xm[i1]<< // the x center of track
4381 // General cut variables
4382 "dfi="<<dfi<< // distance in fi angle
4383 "dtheta="<<dtheta<< // distance int theta angle
4386 "dist0="<<dist[4][0]<< //distance x
4387 "dist1="<<dist[4][1]<< //distance y
4388 "dist2="<<dist[4][2]<< //distance z
4389 "mdist0="<<mdist[0]<< //distance x
4390 "mdist1="<<mdist[1]<< //distance y
4391 "mdist2="<<mdist[2]<< //distance z
4394 delete array->RemoveAt(i1);
4399 AliInfo("Time for splitted tracks removal");
4405 void AliTPCtrackerMI::FindCurling(TObjArray * array, AliESDEvent */*esd*/, Int_t iter)
4408 // find Curling tracks
4409 // Use AliTPCReconstructor::StreamLevel()>1 if you want to tune parameters - cuts
4412 // Algorithm done in 2 phases - because of CPU consumption
4413 // it is n^2 algorithm - for lead-lead 20000x20000 combination are investigated
4414 // see detal in MC part what can be used to cut
4418 const Float_t kMaxC = 400; // maximal curvature to of the track
4419 const Float_t kMaxdTheta = 0.15; // maximal distance in theta
4420 const Float_t kMaxdPhi = 0.15; // maximal distance in phi
4421 const Float_t kPtRatio = 0.3; // ratio between pt
4422 const Float_t kMinDCAR = 2.; // distance to the primary vertex in r - see cpipe cut
4425 // Curling tracks cuts
4428 const Float_t kMaxDeltaRMax = 40; // distance in outer radius
4429 const Float_t kMaxDeltaRMin = 5.; // distance in lower radius - see cpipe cut
4430 const Float_t kMinAngle = 2.9; // angle between tracks
4431 const Float_t kMaxDist = 5; // biggest distance
4433 // The cuts can be tuned using the "MC information stored in Multi tree ==> see FindMultiMC
4436 TCut csign("csign","Tr0.fP[4]*Tr1.fP[4]<0"); //opposite sign
4437 TCut cmax("cmax","abs(Tr0.GetC())>1/400");
4438 TCut cda("cda","sqrt(dtheta^2+dfi^2)<0.15");
4439 TCut ccratio("ccratio","abs((Tr0.fP[4]+Tr1.fP[4])/(abs(Tr0.fP[4])+abs(Tr1.fP[4])))<0.3");
4440 TCut cpipe("cpipe", "min(abs(r0-rc0),abs(r1-rc1))>5");
4442 TCut cdrmax("cdrmax","abs(abs(rc0+r0)-abs(rc1+r1))<40")
4443 TCut cdrmin("cdrmin","abs(abs(rc0+r0)-abs(rc1+r1))<10")
4445 Multi->Draw("dfi","iter==0"+csign+cmax+cda+ccratio); ~94% of curling tracks fulfill
4446 Multi->Draw("min(abs(r0-rc0),abs(r1-rc1))","iter==0&&abs(lab1)==abs(lab0)"+csign+cmax+cda+ccratio+cpipe+cdrmin+cdrmax); //80%
4448 Curling2->Draw("dfi","iter==0&&abs(lab0)==abs(lab1)"+csign+cmax+cdtheta+cdfi+ccratio)
4454 Int_t nentries = array->GetEntriesFast();
4455 AliHelix *helixes = new AliHelix[nentries];
4456 for (Int_t i=0;i<nentries;i++){
4457 AliTPCseed* track = (AliTPCseed*)array->At(i);
4458 if (!track) continue;
4459 track->SetCircular(0);
4460 new (&helixes[i]) AliHelix(*track);
4466 Double_t phase[2][2],radius[2];
4470 TTreeSRedirector &cstream = *fDebugStreamer;
4472 for (Int_t i0=0;i0<nentries;i0++){
4473 AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
4474 if (!track0) continue;
4475 if (TMath::Abs(track0->GetC())<1/kMaxC) continue;
4476 Float_t xc0 = helixes[i0].GetHelix(6);
4477 Float_t yc0 = helixes[i0].GetHelix(7);
4478 Float_t r0 = helixes[i0].GetHelix(8);
4479 Float_t rc0 = TMath::Sqrt(xc0*xc0+yc0*yc0);
4480 Float_t fi0 = TMath::ATan2(yc0,xc0);
4482 for (Int_t i1=i0+1;i1<nentries;i1++){
4483 AliTPCseed * track1 = (AliTPCseed*)array->At(i1);
4484 if (!track1) continue;
4485 if (TMath::Abs(track1->GetC())<1/kMaxC) continue;
4486 Float_t xc1 = helixes[i1].GetHelix(6);
4487 Float_t yc1 = helixes[i1].GetHelix(7);
4488 Float_t r1 = helixes[i1].GetHelix(8);
4489 Float_t rc1 = TMath::Sqrt(xc1*xc1+yc1*yc1);
4490 Float_t fi1 = TMath::ATan2(yc1,xc1);
4492 Float_t dfi = fi0-fi1;
4495 if (dfi>1.5*TMath::Pi()) dfi-=TMath::Pi(); // take care about edge effect
4496 if (dfi<-1.5*TMath::Pi()) dfi+=TMath::Pi(); //
4497 Float_t dtheta = TMath::Abs(track0->GetTgl()-track1->GetTgl())<TMath::Abs(track0->GetTgl()+track1->GetTgl())? track0->GetTgl()-track1->GetTgl():track0->GetTgl()+track1->GetTgl();
4501 if (track0->GetBConstrain()&&track1->GetBConstrain()) continue; // not constrained
4502 if (track1->GetSigned1Pt()*track0->GetSigned1Pt()>0) continue; // not the same sign
4503 if ( TMath::Abs(track1->GetTgl()+track0->GetTgl())>kMaxdTheta) continue; //distance in the Theta
4504 if ( TMath::Abs(dfi)>kMaxdPhi) continue; //distance in phi
4505 if ( TMath::Sqrt(dfi*dfi+dtheta*dtheta)>kMaxdPhi) continue; //common angular offset
4507 Float_t pt0 = track0->GetSignedPt();
4508 Float_t pt1 = track1->GetSignedPt();
4509 if ((TMath::Abs(pt0+pt1)/(TMath::Abs(pt0)+TMath::Abs(pt1)))>kPtRatio) continue;
4510 if ((iter==1) && TMath::Abs(TMath::Abs(rc0+r0)-TMath::Abs(rc1+r1))>kMaxDeltaRMax) continue;
4511 if ((iter!=1) &&TMath::Abs(TMath::Abs(rc0-r0)-TMath::Abs(rc1-r1))>kMaxDeltaRMin) continue;
4512 if (TMath::Min(TMath::Abs(rc0-r0),TMath::Abs(rc1-r1))<kMinDCAR) continue;
4515 // Now find closest approach
4519 Int_t npoints = helixes[i0].GetRPHIintersections(helixes[i1], phase, radius,10);
4520 if (npoints==0) continue;
4521 helixes[i0].GetClosestPhases(helixes[i1], phase);
4525 Double_t hangles[3];
4526 helixes[i0].Evaluate(phase[0][0],xyz0);
4527 helixes[i1].Evaluate(phase[0][1],xyz1);
4529 helixes[i0].GetAngle(phase[0][0],helixes[i1],phase[0][1],hangles);
4530 Double_t deltah[2],deltabest;
4531 if (TMath::Abs(hangles[2])<kMinAngle) continue;
4535 helixes[i0].ParabolicDCA(helixes[i1],phase[0][0],phase[0][1],radius[0],deltah[0],2);
4537 helixes[i0].ParabolicDCA(helixes[i1],phase[1][0],phase[1][1],radius[1],deltah[1],2);
4538 if (deltah[1]<deltah[0]) ibest=1;
4540 deltabest = TMath::Sqrt(deltah[ibest]);
4541 helixes[i0].Evaluate(phase[ibest][0],xyz0);
4542 helixes[i1].Evaluate(phase[ibest][1],xyz1);
4543 helixes[i0].GetAngle(phase[ibest][0],helixes[i1],phase[ibest][1],hangles);
4544 Double_t radiusbest = TMath::Sqrt(radius[ibest]);
4546 if (deltabest>kMaxDist) continue;
4547 // if (mindcar+mindcaz<40 && (TMath::Abs(hangles[2])<kMinAngle ||deltabest>3)) continue;
4548 Bool_t sign =kFALSE;
4549 if (hangles[2]>kMinAngle) sign =kTRUE;
4552 // circular[i0] = kTRUE;
4553 // circular[i1] = kTRUE;
4554 if (track0->OneOverPt()<track1->OneOverPt()){
4555 track0->SetCircular(track0->GetCircular()+1);
4556 track1->SetCircular(track1->GetCircular()+2);
4559 track1->SetCircular(track1->GetCircular()+1);
4560 track0->SetCircular(track0->GetCircular()+2);
4563 if (AliTPCReconstructor::StreamLevel()>1){
4565 //debug stream to tune "fine" cuts
4566 Int_t lab0=track0->GetLabel();
4567 Int_t lab1=track1->GetLabel();
4568 cstream<<"Curling2"<<
4584 "npoints="<<npoints<<
4585 "hangles0="<<hangles[0]<<
4586 "hangles1="<<hangles[1]<<
4587 "hangles2="<<hangles[2]<<
4590 "radius="<<radiusbest<<
4591 "deltabest="<<deltabest<<
4592 "phase0="<<phase[ibest][0]<<
4593 "phase1="<<phase[ibest][1]<<
4601 if (AliTPCReconstructor::StreamLevel()>1) {
4602 AliInfo("Time for curling tracks removal");
4611 void AliTPCtrackerMI::FindKinks(TObjArray * array, AliESDEvent *esd)
4618 TObjArray *kinks= new TObjArray(10000);
4619 // TObjArray *v0s= new TObjArray(10000);
4620 Int_t nentries = array->GetEntriesFast();
4621 AliHelix *helixes = new AliHelix[nentries];
4622 Int_t *sign = new Int_t[nentries];
4623 Int_t *nclusters = new Int_t[nentries];
4624 Float_t *alpha = new Float_t[nentries];
4625 AliKink *kink = new AliKink();
4626 Int_t * usage = new Int_t[nentries];
4627 Float_t *zm = new Float_t[nentries];
4628 Float_t *z0 = new Float_t[nentries];
4629 Float_t *fim = new Float_t[nentries];
4630 Float_t *shared = new Float_t[nentries];
4631 Bool_t *circular = new Bool_t[nentries];
4632 Float_t *dca = new Float_t[nentries];
4633 //const AliESDVertex * primvertex = esd->GetVertex();
4635 // nentries = array->GetEntriesFast();
4640 for (Int_t i=0;i<nentries;i++){
4643 AliTPCseed* track = (AliTPCseed*)array->At(i);
4644 if (!track) continue;
4645 track->SetCircular(0);
4647 track->UpdatePoints();
4648 if (( track->GetPoints()[2]- track->GetPoints()[0])>5 && track->GetPoints()[3]>0.8){
4650 nclusters[i]=track->GetNumberOfClusters();
4651 alpha[i] = track->GetAlpha();
4652 new (&helixes[i]) AliHelix(*track);
4654 helixes[i].Evaluate(0,xyz);
4655 sign[i] = (track->GetC()>0) ? -1:1;
4658 if (track->GetProlongation(x,y,z)){
4660 fim[i] = alpha[i]+TMath::ATan2(y,x);
4663 zm[i] = track->GetZ();
4667 circular[i]= kFALSE;
4668 if (track->GetProlongation(0,y,z)) z0[i] = z;
4669 dca[i] = track->GetD(0,0);
4675 Int_t ncandidates =0;
4678 Double_t phase[2][2],radius[2];
4681 // Find circling track
4682 TTreeSRedirector &cstream = *fDebugStreamer;
4684 for (Int_t i0=0;i0<nentries;i0++){
4685 AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
4686 if (!track0) continue;
4687 if (track0->GetNumberOfClusters()<40) continue;
4688 if (TMath::Abs(1./track0->GetC())>200) continue;
4689 for (Int_t i1=i0+1;i1<nentries;i1++){
4690 AliTPCseed * track1 = (AliTPCseed*)array->At(i1);
4691 if (!track1) continue;
4692 if (track1->GetNumberOfClusters()<40) continue;
4693 if ( TMath::Abs(track1->GetTgl()+track0->GetTgl())>0.1) continue;
4694 if (track0->GetBConstrain()&&track1->GetBConstrain()) continue;
4695 if (TMath::Abs(1./track1->GetC())>200) continue;
4696 if (track1->GetSigned1Pt()*track0->GetSigned1Pt()>0) continue;
4697 if (track1->GetTgl()*track0->GetTgl()>0) continue;
4698 if (TMath::Max(TMath::Abs(1./track0->GetC()),TMath::Abs(1./track1->GetC()))>190) continue;
4699 if (track0->GetBConstrain()&&track1->OneOverPt()<track0->OneOverPt()) continue; //returning - lower momenta
4700 if (track1->GetBConstrain()&&track0->OneOverPt()<track1->OneOverPt()) continue; //returning - lower momenta
4702 Float_t mindcar = TMath::Min(TMath::Abs(dca[i0]),TMath::Abs(dca[i1]));
4703 if (mindcar<5) continue;
4704 Float_t mindcaz = TMath::Min(TMath::Abs(z0[i0]-GetZ()),TMath::Abs(z0[i1]-GetZ()));
4705 if (mindcaz<5) continue;
4706 if (mindcar+mindcaz<20) continue;
4709 Float_t xc0 = helixes[i0].GetHelix(6);
4710 Float_t yc0 = helixes[i0].GetHelix(7);
4711 Float_t r0 = helixes[i0].GetHelix(8);
4712 Float_t xc1 = helixes[i1].GetHelix(6);
4713 Float_t yc1 = helixes[i1].GetHelix(7);
4714 Float_t r1 = helixes[i1].GetHelix(8);
4716 Float_t rmean = (r0+r1)*0.5;
4717 Float_t delta =TMath::Sqrt((xc1-xc0)*(xc1-xc0)+(yc1-yc0)*(yc1-yc0));
4718 //if (delta>30) continue;
4719 if (delta>rmean*0.25) continue;
4720 if (TMath::Abs(r0-r1)/rmean>0.3) continue;
4722 Int_t npoints = helixes[i0].GetRPHIintersections(helixes[i1], phase, radius,10);
4723 if (npoints==0) continue;
4724 helixes[i0].GetClosestPhases(helixes[i1], phase);
4728 Double_t hangles[3];
4729 helixes[i0].Evaluate(phase[0][0],xyz0);
4730 helixes[i1].Evaluate(phase[0][1],xyz1);
4732 helixes[i0].GetAngle(phase[0][0],helixes[i1],phase[0][1],hangles);
4733 Double_t deltah[2],deltabest;
4734 if (hangles[2]<2.8) continue;
4737 helixes[i0].ParabolicDCA(helixes[i1],phase[0][0],phase[0][1],radius[0],deltah[0],2);
4739 helixes[i0].ParabolicDCA(helixes[i1],phase[1][0],phase[1][1],radius[1],deltah[1],2);
4740 if (deltah[1]<deltah[0]) ibest=1;
4742 deltabest = TMath::Sqrt(deltah[ibest]);
4743 helixes[i0].Evaluate(phase[ibest][0],xyz0);
4744 helixes[i1].Evaluate(phase[ibest][1],xyz1);
4745 helixes[i0].GetAngle(phase[ibest][0],helixes[i1],phase[ibest][1],hangles);
4746 Double_t radiusbest = TMath::Sqrt(radius[ibest]);
4748 if (deltabest>6) continue;
4749 if (mindcar+mindcaz<40 && (hangles[2]<3.12||deltabest>3)) continue;
4750 Bool_t sign =kFALSE;
4751 if (hangles[2]>3.06) sign =kTRUE;
4754 circular[i0] = kTRUE;
4755 circular[i1] = kTRUE;
4756 if (track0->OneOverPt()<track1->OneOverPt()){
4757 track0->SetCircular(track0->GetCircular()+1);
4758 track1->SetCircular(track1->GetCircular()+2);
4761 track1->SetCircular(track1->GetCircular()+1);
4762 track0->SetCircular(track0->GetCircular()+2);
4765 if (sign&&AliTPCReconstructor::StreamLevel()>1){
4767 Int_t lab0=track0->GetLabel();
4768 Int_t lab1=track1->GetLabel();
4769 cstream<<"Curling"<<
4776 "mindcar="<<mindcar<<
4777 "mindcaz="<<mindcaz<<
4780 "npoints="<<npoints<<
4781 "hangles0="<<hangles[0]<<
4782 "hangles2="<<hangles[2]<<
4787 "radius="<<radiusbest<<
4788 "deltabest="<<deltabest<<
4789 "phase0="<<phase[ibest][0]<<
4790 "phase1="<<phase[ibest][1]<<
4800 for (Int_t i =0;i<nentries;i++){
4801 if (sign[i]==0) continue;
4802 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
4809 Double_t cradius0 = 40*40;
4810 Double_t cradius1 = 270*270;
4813 Double_t cdist3=0.55;
4814 for (Int_t j =i+1;j<nentries;j++){
4816 if (sign[j]*sign[i]<1) continue;
4817 if ( (nclusters[i]+nclusters[j])>200) continue;
4818 if ( (nclusters[i]+nclusters[j])<80) continue;
4819 if ( TMath::Abs(zm[i]-zm[j])>60.) continue;
4820 if ( TMath::Abs(fim[i]-fim[j])>0.6 && TMath::Abs(fim[i]-fim[j])<5.7 ) continue;
4821 //AliTPCseed * track1 = (AliTPCseed*)array->At(j); Double_t phase[2][2],radius[2];
4822 Int_t npoints = helixes[i].GetRPHIintersections(helixes[j], phase, radius,20);
4823 if (npoints<1) continue;
4826 if (radius[0]<cradius0||radius[0]>cradius1) continue;
4829 if ( (radius[0]<cradius0||radius[0]>cradius1) && (radius[1]<cradius0||radius[1]>cradius1) ) continue;
4832 Double_t delta1=10000,delta2=10000;
4833 // cuts on the intersection radius
4834 helixes[i].LinearDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta1);
4835 if (radius[0]<20&&delta1<1) continue; //intersection at vertex
4836 if (radius[0]<10&&delta1<3) continue; //intersection at vertex
4838 helixes[i].LinearDCA(helixes[j],phase[1][0],phase[1][1],radius[1],delta2);
4839 if (radius[1]<20&&delta2<1) continue; //intersection at vertex
4840 if (radius[1]<10&&delta2<3) continue; //intersection at vertex
4843 Double_t distance1 = TMath::Min(delta1,delta2);
4844 if (distance1>cdist1) continue; // cut on DCA linear approximation
4846 npoints = helixes[i].GetRPHIintersections(helixes[j], phase, radius,20);
4847 helixes[i].ParabolicDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta1);
4848 if (radius[0]<20&&delta1<1) continue; //intersection at vertex
4849 if (radius[0]<10&&delta1<3) continue; //intersection at vertex
4852 helixes[i].ParabolicDCA(helixes[j],phase[1][0],phase[1][1],radius[1],delta2);
4853 if (radius[1]<20&&delta2<1) continue; //intersection at vertex
4854 if (radius[1]<10&&delta2<3) continue; //intersection at vertex
4856 distance1 = TMath::Min(delta1,delta2);
4859 rkink = TMath::Sqrt(radius[0]);
4862 rkink = TMath::Sqrt(radius[1]);
4864 if (distance1>cdist2) continue;
4867 AliTPCseed * track1 = (AliTPCseed*)array->At(j);
4870 Int_t row0 = GetRowNumber(rkink);
4871 if (row0<10) continue;
4872 if (row0>150) continue;
4875 Float_t dens00=-1,dens01=-1;
4876 Float_t dens10=-1,dens11=-1;
4878 Int_t found,foundable,shared;
4879 track0->GetClusterStatistic(0,row0-5, found, foundable,shared,kFALSE);
4880 if (foundable>5) dens00 = Float_t(found)/Float_t(foundable);
4881 track0->GetClusterStatistic(row0+5,155, found, foundable,shared,kFALSE);
4882 if (foundable>5) dens01 = Float_t(found)/Float_t(foundable);
4884 track1->GetClusterStatistic(0,row0-5, found, foundable,shared,kFALSE);
4885 if (foundable>10) dens10 = Float_t(found)/Float_t(foundable);
4886 track1->GetClusterStatistic(row0+5,155, found, foundable,shared,kFALSE);
4887 if (foundable>10) dens11 = Float_t(found)/Float_t(foundable);
4889 if (dens00<dens10 && dens01<dens11) continue;
4890 if (dens00>dens10 && dens01>dens11) continue;
4891 if (TMath::Max(dens00,dens10)<0.1) continue;
4892 if (TMath::Max(dens01,dens11)<0.3) continue;
4894 if (TMath::Min(dens00,dens10)>0.6) continue;
4895 if (TMath::Min(dens01,dens11)>0.6) continue;
4898 AliTPCseed * ktrack0, *ktrack1;
4907 if (TMath::Abs(ktrack0->GetC())>5) continue; // cut on the curvature for mother particle
4908 AliExternalTrackParam paramm(*ktrack0);
4909 AliExternalTrackParam paramd(*ktrack1);
4910 if (row0>60&&ktrack1->GetReference().GetX()>90.)new (¶md) AliExternalTrackParam(ktrack1->GetReference());
4913 kink->SetMother(paramm);
4914 kink->SetDaughter(paramd);
4917 Float_t x[3] = { kink->GetPosition()[0],kink->GetPosition()[1],kink->GetPosition()[2]};
4919 fParam->Transform0to1(x,index);
4920 fParam->Transform1to2(x,index);
4921 row0 = GetRowNumber(x[0]);
4923 if (kink->GetR()<100) continue;
4924 if (kink->GetR()>240) continue;
4925 if (kink->GetPosition()[2]/kink->GetR()>AliTPCReconstructor::GetCtgRange()) continue; //out of fiducial volume
4926 if (kink->GetDistance()>cdist3) continue;
4927 Float_t dird = kink->GetDaughterP()[0]*kink->GetPosition()[0]+kink->GetDaughterP()[1]*kink->GetPosition()[1]; // rough direction estimate
4928 if (dird<0) continue;
4930 Float_t dirm = kink->GetMotherP()[0]*kink->GetPosition()[0]+kink->GetMotherP()[1]*kink->GetPosition()[1]; // rough direction estimate
4931 if (dirm<0) continue;
4932 Float_t mpt = TMath::Sqrt(kink->GetMotherP()[0]*kink->GetMotherP()[0]+kink->GetMotherP()[1]*kink->GetMotherP()[1]);
4933 if (mpt<0.2) continue;
4936 //for high momenta momentum not defined well in first iteration
4937 Double_t qt = TMath::Sin(kink->GetAngle(2))*ktrack1->GetP();
4938 if (qt>0.35) continue;
4941 kink->SetLabel(CookLabel(ktrack0,0.4,0,row0),0);
4942 kink->SetLabel(CookLabel(ktrack1,0.4,row0,160),1);
4944 kink->SetTPCDensity(dens00,0,0);
4945 kink->SetTPCDensity(dens01,0,1);
4946 kink->SetTPCDensity(dens10,1,0);
4947 kink->SetTPCDensity(dens11,1,1);
4948 kink->SetIndex(i,0);
4949 kink->SetIndex(j,1);
4952 kink->SetTPCDensity(dens10,0,0);
4953 kink->SetTPCDensity(dens11,0,1);
4954 kink->SetTPCDensity(dens00,1,0);
4955 kink->SetTPCDensity(dens01,1,1);
4956 kink->SetIndex(j,0);
4957 kink->SetIndex(i,1);
4960 if (mpt<1||kink->GetAngle(2)>0.1){
4961 // angle and densities not defined yet
4962 if (kink->GetTPCDensityFactor()<0.8) continue;
4963 if ((2-kink->GetTPCDensityFactor())*kink->GetDistance() >0.25) continue;
4964 if (kink->GetAngle(2)*ktrack0->GetP()<0.003) continue; //too small angle
4965 if (kink->GetAngle(2)>0.2&&kink->GetTPCDensityFactor()<1.15) continue;
4966 if (kink->GetAngle(2)>0.2&&kink->GetTPCDensity(0,1)>0.05) continue;
4968 Float_t criticalangle = track0->GetSigmaSnp2()+track0->GetSigmaTgl2();
4969 criticalangle+= track1->GetSigmaSnp2()+track1->GetSigmaTgl2();
4970 criticalangle= 3*TMath::Sqrt(criticalangle);
4971 if (criticalangle>0.02) criticalangle=0.02;
4972 if (kink->GetAngle(2)<criticalangle) continue;
4975 Int_t drow = Int_t(2.+0.5/(0.05+kink->GetAngle(2))); // overlap region defined
4976 Float_t shapesum =0;
4978 for ( Int_t row = row0-drow; row<row0+drow;row++){
4979 if (row<0) continue;
4980 if (row>155) continue;
4981 if (ktrack0->GetClusterPointer(row)){
4982 AliTPCTrackerPoint *point =ktrack0->GetTrackPoint(row);
4983 shapesum+=point->GetSigmaY()+point->GetSigmaZ();
4986 if (ktrack1->GetClusterPointer(row)){
4987 AliTPCTrackerPoint *point =ktrack1->GetTrackPoint(row);
4988 shapesum+=point->GetSigmaY()+point->GetSigmaZ();
4993 kink->SetShapeFactor(-1.);
4996 kink->SetShapeFactor(shapesum/sum);
4998 // esd->AddKink(kink);
4999 kinks->AddLast(kink);
5005 // sort the kinks according quality - and refit them towards vertex
5007 Int_t nkinks = kinks->GetEntriesFast();
5008 Float_t *quality = new Float_t[nkinks];
5009 Int_t *indexes = new Int_t[nkinks];
5010 AliTPCseed *mothers = new AliTPCseed[nkinks];
5011 AliTPCseed *daughters = new AliTPCseed[nkinks];
5014 for (Int_t i=0;i<nkinks;i++){
5016 AliKink *kink = (AliKink*)kinks->At(i);
5018 // refit kinks towards vertex
5020 Int_t index0 = kink->GetIndex(0);
5021 Int_t index1 = kink->GetIndex(1);
5022 AliTPCseed * ktrack0 = (AliTPCseed*)array->At(index0);
5023 AliTPCseed * ktrack1 = (AliTPCseed*)array->At(index1);
5025 Int_t sumn=ktrack0->GetNumberOfClusters()+ktrack1->GetNumberOfClusters();
5027 // Refit Kink under if too small angle
5029 if (kink->GetAngle(2)<0.05){
5030 kink->SetTPCRow0(GetRowNumber(kink->GetR()));
5031 Int_t row0 = kink->GetTPCRow0();
5032 Int_t drow = Int_t(2.+0.5/(0.05+kink->GetAngle(2)));
5035 Int_t last = row0-drow;
5036 if (last<40) last=40;
5037 if (last<ktrack0->GetFirstPoint()+25) last = ktrack0->GetFirstPoint()+25;
5038 AliTPCseed* seed0 = ReSeed(ktrack0,last,kFALSE);
5041 Int_t first = row0+drow;
5042 if (first>130) first=130;
5043 if (first>ktrack1->GetLastPoint()-25) first = TMath::Max(ktrack1->GetLastPoint()-25,30);
5044 AliTPCseed* seed1 = ReSeed(ktrack1,first,kTRUE);
5046 if (seed0 && seed1){
5047 kink->SetStatus(1,8);
5048 if (RefitKink(*seed0,*seed1,*kink)) kink->SetStatus(1,9);
5049 row0 = GetRowNumber(kink->GetR());
5050 sumn = seed0->GetNumberOfClusters()+seed1->GetNumberOfClusters();
5051 mothers[i] = *seed0;
5052 daughters[i] = *seed1;
5055 delete kinks->RemoveAt(i);
5056 if (seed0) delete seed0;
5057 if (seed1) delete seed1;
5060 if (kink->GetDistance()>0.5 || kink->GetR()<110 || kink->GetR()>240) {
5061 delete kinks->RemoveAt(i);
5062 if (seed0) delete seed0;
5063 if (seed1) delete seed1;
5071 if (kink) quality[i] = 160*((0.1+kink->GetDistance())*(2.-kink->GetTPCDensityFactor()))/(sumn+40.); //the longest -clossest will win
5073 TMath::Sort(nkinks,quality,indexes,kFALSE);
5075 //remove double find kinks
5077 for (Int_t ikink0=1;ikink0<nkinks;ikink0++){
5078 AliKink * kink0 = (AliKink*) kinks->At(indexes[ikink0]);
5079 if (!kink0) continue;
5081 for (Int_t ikink1=0;ikink1<ikink0;ikink1++){
5082 if (!kink0) continue;
5083 AliKink * kink1 = (AliKink*) kinks->At(indexes[ikink1]);
5084 if (!kink1) continue;
5085 // if not close kink continue
5086 if (TMath::Abs(kink1->GetPosition()[2]-kink0->GetPosition()[2])>10) continue;
5087 if (TMath::Abs(kink1->GetPosition()[1]-kink0->GetPosition()[1])>10) continue;
5088 if (TMath::Abs(kink1->GetPosition()[0]-kink0->GetPosition()[0])>10) continue;
5090 AliTPCseed &mother0 = mothers[indexes[ikink0]];
5091 AliTPCseed &daughter0 = daughters[indexes[ikink0]];
5092 AliTPCseed &mother1 = mothers[indexes[ikink1]];
5093 AliTPCseed &daughter1 = daughters[indexes[ikink1]];
5094 Int_t row0 = (kink0->GetTPCRow0()+kink1->GetTPCRow0())/2;
5103 for (Int_t i=0;i<row0;i++){
5104 if (mother0.GetClusterIndex(i)>0 && mother1.GetClusterIndex(i)>0){
5107 if (mother0.GetClusterIndex(i)==mother1.GetClusterIndex(i)){
5114 for (Int_t i=row0;i<158;i++){
5115 if (daughter0.GetClusterIndex(i)>0 && daughter0.GetClusterIndex(i)>0){
5118 if (mother0.GetClusterIndex(i)==mother1.GetClusterIndex(i)){
5124 Float_t ratio = Float_t(same+1)/Float_t(both+1);
5125 Float_t ratiom = Float_t(samem+1)/Float_t(bothm+1);
5126 Float_t ratiod = Float_t(samed+1)/Float_t(bothd+1);
5127 if (ratio>0.3 && ratiom>0.5 &&ratiod>0.5) {
5128 Int_t sum0 = mother0.GetNumberOfClusters()+daughter0.GetNumberOfClusters();
5129 Int_t sum1 = mother1.GetNumberOfClusters()+daughter1.GetNumberOfClusters();
5131 shared[kink0->GetIndex(0)]= kTRUE;
5132 shared[kink0->GetIndex(1)]= kTRUE;
5133 delete kinks->RemoveAt(indexes[ikink0]);
5136 shared[kink1->GetIndex(0)]= kTRUE;
5137 shared[kink1->GetIndex(1)]= kTRUE;
5138 delete kinks->RemoveAt(indexes[ikink1]);
5145 for (Int_t i=0;i<nkinks;i++){
5146 AliKink * kink = (AliKink*) kinks->At(indexes[i]);
5147 if (!kink) continue;
5148 kink->SetTPCRow0(GetRowNumber(kink->GetR()));
5149 Int_t index0 = kink->GetIndex(0);
5150 Int_t index1 = kink->GetIndex(1);
5151 if (circular[index0]||circular[index1]&&kink->GetDistance()>0.2) continue;
5152 kink->SetMultiple(usage[index0],0);
5153 kink->SetMultiple(usage[index1],1);
5154 if (kink->GetMultiple()[0]+kink->GetMultiple()[1]>2) continue;
5155 if (kink->GetMultiple()[0]+kink->GetMultiple()[1]>0 && quality[indexes[i]]>0.2) continue;
5156 if (kink->GetMultiple()[0]+kink->GetMultiple()[1]>0 && kink->GetDistance()>0.2) continue;
5157 if (circular[index0]||circular[index1]&&kink->GetDistance()>0.1) continue;
5159 AliTPCseed * ktrack0 = (AliTPCseed*)array->At(index0);
5160 AliTPCseed * ktrack1 = (AliTPCseed*)array->At(index1);
5161 if (!ktrack0 || !ktrack1) continue;
5162 Int_t index = esd->AddKink(kink);
5165 if ( ktrack0->GetKinkIndex(0)==0 && ktrack1->GetKinkIndex(0)==0) { //best kink
5166 if (mothers[indexes[i]].GetNumberOfClusters()>20 && daughters[indexes[i]].GetNumberOfClusters()>20 && (mothers[indexes[i]].GetNumberOfClusters()+daughters[indexes[i]].GetNumberOfClusters())>100){
5167 *ktrack0 = mothers[indexes[i]];
5168 *ktrack1 = daughters[indexes[i]];
5172 ktrack0->SetKinkIndex(usage[index0],-(index+1));
5173 ktrack1->SetKinkIndex(usage[index1], (index+1));
5178 // Remove tracks corresponding to shared kink's
5180 for (Int_t i=0;i<nentries;i++){
5181 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5182 if (!track0) continue;
5183 if (track0->GetKinkIndex(0)!=0) continue;
5184 if (shared[i]) delete array->RemoveAt(i);
5189 RemoveUsed2(array,0.5,0.4,30);
5191 for (Int_t i=0;i<nentries;i++){
5192 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5193 if (!track0) continue;
5194 track0->CookdEdx(0.02,0.6);
5198 for (Int_t i=0;i<nentries;i++){
5199 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5200 if (!track0) continue;
5201 if (track0->Pt()<1.4) continue;
5202 //remove double high momenta tracks - overlapped with kink candidates
5205 for (Int_t icl=track0->GetFirstPoint();icl<track0->GetLastPoint(); icl++){
5206 if (track0->GetClusterPointer(icl)!=0){
5208 if (track0->GetClusterPointer(icl)->IsUsed(10)) shared++;
5211 if (Float_t(shared+1)/Float_t(all+1)>0.5) {
5212 delete array->RemoveAt(i);
5216 if (track0->GetKinkIndex(0)!=0) continue;
5217 if (track0->GetNumberOfClusters()<80) continue;
5219 AliTPCseed *pmother = new AliTPCseed();
5220 AliTPCseed *pdaughter = new AliTPCseed();
5221 AliKink *pkink = new AliKink;
5223 AliTPCseed & mother = *pmother;
5224 AliTPCseed & daughter = *pdaughter;
5225 AliKink & kink = *pkink;
5226 if (CheckKinkPoint(track0,mother,daughter, kink)){
5227 if (mother.GetNumberOfClusters()<30||daughter.GetNumberOfClusters()<20) {
5231 continue; //too short tracks
5233 if (mother.Pt()<1.4) {
5239 Int_t row0= kink.GetTPCRow0();
5240 if (kink.GetDistance()>0.5 || kink.GetR()<110. || kink.GetR()>240.) {
5247 Int_t index = esd->AddKink(&kink);
5248 mother.SetKinkIndex(0,-(index+1));
5249 daughter.SetKinkIndex(0,index+1);
5250 if (mother.GetNumberOfClusters()>50) {
5251 delete array->RemoveAt(i);
5252 array->AddAt(new AliTPCseed(mother),i);
5255 array->AddLast(new AliTPCseed(mother));
5257 array->AddLast(new AliTPCseed(daughter));
5258 for (Int_t icl=0;icl<row0;icl++) {
5259 if (mother.GetClusterPointer(icl)) mother.GetClusterPointer(icl)->Use(20);
5262 for (Int_t icl=row0;icl<158;icl++) {
5263 if (daughter.GetClusterPointer(icl)) daughter.GetClusterPointer(icl)->Use(20);
5272 delete [] daughters;
5294 printf("Ncandidates=\t%d\t%d\t%d\t%d\n",esd->GetNumberOfKinks(),ncandidates,ntracks,nall);
5298 void AliTPCtrackerMI::FindV0s(TObjArray * array, AliESDEvent *esd)
5304 TObjArray *tpcv0s = new TObjArray(100000);
5305 Int_t nentries = array->GetEntriesFast();
5306 AliHelix *helixes = new AliHelix[nentries];
5307 Int_t *sign = new Int_t[nentries];
5308 Float_t *alpha = new Float_t[nentries];
5309 Float_t *z0 = new Float_t[nentries];
5310 Float_t *dca = new Float_t[nentries];
5311 Float_t *sdcar = new Float_t[nentries];
5312 Float_t *cdcar = new Float_t[nentries];
5313 Float_t *pulldcar = new Float_t[nentries];
5314 Float_t *pulldcaz = new Float_t[nentries];
5315 Float_t *pulldca = new Float_t[nentries];
5316 Bool_t *isPrim = new Bool_t[nentries];
5317 const AliESDVertex * primvertex = esd->GetVertex();
5318 Double_t zvertex = primvertex->GetZv();
5320 // nentries = array->GetEntriesFast();
5322 for (Int_t i=0;i<nentries;i++){
5325 AliTPCseed* track = (AliTPCseed*)array->At(i);
5326 if (!track) continue;
5327 track->GetV0Indexes()[0] = 0; //rest v0 indexes
5328 track->GetV0Indexes()[1] = 0; //rest v0 indexes
5329 track->GetV0Indexes()[2] = 0; //rest v0 indexes
5331 alpha[i] = track->GetAlpha();
5332 new (&helixes[i]) AliHelix(*track);
5334 helixes[i].Evaluate(0,xyz);
5335 sign[i] = (track->GetC()>0) ? -1:1;
5339 if (track->GetProlongation(0,y,z)) z0[i] = z;
5340 dca[i] = track->GetD(0,0);
5342 // dca error parrameterezation + pulls
5344 sdcar[i] = TMath::Sqrt(0.150*0.150+(100*track->GetC())*(100*track->GetC()));
5345 if (TMath::Abs(track->GetTgl())>1) sdcar[i]*=2.5;
5346 cdcar[i] = TMath::Exp((TMath::Abs(track->GetC())-0.0106)*525.3);
5347 pulldcar[i] = (dca[i]-cdcar[i])/sdcar[i];
5348 pulldcaz[i] = (z0[i]-zvertex)/sdcar[i];
5349 pulldca[i] = TMath::Sqrt(pulldcar[i]*pulldcar[i]+pulldcaz[i]*pulldcaz[i]);
5350 if (track->TPCrPID(1)+track->TPCrPID(2)+track->TPCrPID(3)>0.5) {
5351 if (pulldca[i]<3.) isPrim[i]=kTRUE; //pion, muon and Kaon 3 sigma cut
5353 if (track->TPCrPID(4)>0.5) {
5354 if (pulldca[i]<0.5) isPrim[i]=kTRUE; //proton 0.5 sigma cut
5356 if (track->TPCrPID(0)>0.4) {
5357 isPrim[i]=kFALSE; //electron no sigma cut
5364 Int_t ncandidates =0;
5367 Double_t phase[2][2],radius[2];
5373 TTreeSRedirector &cstream = *fDebugStreamer;
5374 Float_t fprimvertex[3]={GetX(),GetY(),GetZ()};
5376 Double_t cradius0 = 10*10;
5377 Double_t cradius1 = 200*200;
5380 Double_t cpointAngle = 0.95;
5382 Double_t delta[2]={10000,10000};
5383 for (Int_t i =0;i<nentries;i++){
5384 if (sign[i]==0) continue;
5385 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5386 if (!track0) continue;
5387 if (AliTPCReconstructor::StreamLevel()>1){
5392 "zvertex="<<zvertex<<
5393 "sdcar0="<<sdcar[i]<<
5394 "cdcar0="<<cdcar[i]<<
5395 "pulldcar0="<<pulldcar[i]<<
5396 "pulldcaz0="<<pulldcaz[i]<<
5397 "pulldca0="<<pulldca[i]<<
5398 "isPrim="<<isPrim[i]<<
5402 if (track0->GetSigned1Pt()<0) continue;
5403 if (track0->GetKinkIndex(0)>0||isPrim[i]) continue; //daughter kink
5405 if (TMath::Abs(helixes[i].GetHelix(4))<0.000000001) continue;
5410 for (Int_t j =0;j<nentries;j++){
5411 AliTPCseed * track1 = (AliTPCseed*)array->At(j);
5412 if (!track1) continue;
5413 if (track1->GetKinkIndex(0)>0 || isPrim[j]) continue; //daughter kink
5414 if (sign[j]*sign[i]>0) continue;
5415 if (TMath::Abs(helixes[j].GetHelix(4))<0.000001) continue;
5416 if (track0->GetCircular()+track1->GetCircular()>1) continue; //circling -returning track
5419 // DCA to prim vertex cut
5425 Int_t npoints = helixes[i].GetRPHIintersections(helixes[j], phase, radius,cdist2);
5426 if (npoints<1) continue;
5430 if (radius[0]<cradius0||radius[0]>cradius1) continue;
5431 helixes[i].LinearDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta[0]);
5432 if (delta[0]>cdist1) continue;
5435 if (TMath::Max(radius[0],radius[1])<cradius0|| TMath::Min(radius[0],radius[1])>cradius1) continue;
5436 helixes[i].LinearDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta[0]);
5437 helixes[i].LinearDCA(helixes[j],phase[1][0],phase[1][1],radius[1],delta[1]);
5438 if (delta[1]<delta[0]) iclosest=1;
5439 if (delta[iclosest]>cdist1) continue;
5441 helixes[i].ParabolicDCA(helixes[j],phase[iclosest][0],phase[iclosest][1],radius[iclosest],delta[iclosest]);
5442 if (radius[iclosest]<cradius0 || radius[iclosest]>cradius1 || delta[iclosest]>cdist1) continue;
5444 Double_t pointAngle = helixes[i].GetPointAngle(helixes[j],phase[iclosest],fprimvertex);
5445 if (pointAngle<cpointAngle) continue;
5447 Bool_t isGamma = kFALSE;
5448 vertex.SetParamP(*track0); //track0 - plus
5449 vertex.SetParamN(*track1); //track1 - minus
5450 vertex.Update(fprimvertex);
5451 if (track0->TPCrPID(0)>0.3&&track1->TPCrPID(0)>0.3&&vertex.GetAnglep()[2]<0.15) isGamma=kTRUE; // gamma conversion candidate
5452 Double_t pointAngle2 = vertex.GetV0CosineOfPointingAngle();
5454 if (vertex.GetV0CosineOfPointingAngle()<cpointAngle && (!isGamma)) continue;// point angle cut
5455 //Bo: if (vertex.GetDist2()>2&&(!isGamma)) continue; // point angle cut
5456 if (vertex.GetDcaV0Daughters()>2&&(!isGamma)) continue;//Bo: // point angle cut
5457 Float_t sigmae = 0.15*0.15;
5458 if (vertex.GetRr()<80)
5459 sigmae += (sdcar[i]*sdcar[i]+sdcar[j]*sdcar[j])*(1.-vertex.GetRr()/80.)*(1.-vertex.GetRr()/80.);
5460 sigmae+= TMath::Sqrt(sigmae);
5461 //Bo: if (vertex.GetDist2()/sigmae>3.&&(!isGamma)) continue;
5462 if (vertex.GetDcaV0Daughters()/sigmae>3.&&(!isGamma)) continue;
5463 Float_t densb0=0,densb1=0,densa0=0,densa1=0;
5464 Int_t row0 = GetRowNumber(vertex.GetRr());
5466 //Bo: if (vertex.GetDist2()>0.2) continue;
5467 if (vertex.GetDcaV0Daughters()>0.2) continue;
5468 densb0 = track0->Density2(0,row0-5);
5469 densb1 = track1->Density2(0,row0-5);
5470 if (densb0>0.3|| densb1>0.3) continue; //clusters before vertex
5471 densa0 = track0->Density2(row0+5,row0+40);
5472 densa1 = track1->Density2(row0+5,row0+40);
5473 if ((densa0<0.4|| densa1<0.4)&&(!isGamma)) continue; //missing clusters after vertex
5476 densa0 = track0->Density2(0,40); //cluster density
5477 densa1 = track1->Density2(0,40); //cluster density
5478 if ((vertex.GetRr()<80&&densa0+densa1<1.)&&(!isGamma)) continue;
5480 //Bo: vertex.SetLab(0,track0->GetLabel());
5481 //Bo: vertex.SetLab(1,track1->GetLabel());
5482 vertex.SetChi2After((densa0+densa1)*0.5);
5483 vertex.SetChi2Before((densb0+densb1)*0.5);
5484 vertex.SetIndex(0,i);
5485 vertex.SetIndex(1,j);
5486 //Bo: vertex.SetStatus(1); // TPC v0 candidate
5487 vertex.SetOnFlyStatus(2);//Bo: // TPC v0 candidate
5488 //Bo: vertex.SetRp(track0->TPCrPIDs());
5489 //Bo: vertex.SetRm(track1->TPCrPIDs());
5490 tpcv0s->AddLast(new AliESDv0(vertex));
5493 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
5494 Double_t radiusm= (delta[0]<delta[1])? TMath::Sqrt(radius[0]):TMath::Sqrt(radius[1]);
5495 Double_t deltam= (delta[0]<delta[1])? TMath::Sqrt(delta[0]):TMath::Sqrt(delta[1]);
5496 if (AliTPCReconstructor::StreamLevel()>1) {
5497 Int_t lab0=track0->GetLabel();
5498 Int_t lab1=track1->GetLabel();
5499 Char_t c0=track0->GetCircular();
5500 Char_t c1=track1->GetCircular();
5503 "vertex.="<<&vertex<<
5506 "Helix0.="<<&helixes[i]<<
5509 "Helix1.="<<&helixes[j]<<
5510 "pointAngle="<<pointAngle<<
5511 "pointAngle2="<<pointAngle2<<
5516 "zvertex="<<zvertex<<
5519 "npoints="<<npoints<<
5520 "radius0="<<radius[0]<<
5521 "delta0="<<delta[0]<<
5522 "radius1="<<radius[1]<<
5523 "delta1="<<delta[1]<<
5524 "radiusm="<<radiusm<<
5526 "sdcar0="<<sdcar[i]<<
5527 "sdcar1="<<sdcar[j]<<
5528 "cdcar0="<<cdcar[i]<<
5529 "cdcar1="<<cdcar[j]<<
5530 "pulldcar0="<<pulldcar[i]<<
5531 "pulldcar1="<<pulldcar[j]<<
5532 "pulldcaz0="<<pulldcaz[i]<<
5533 "pulldcaz1="<<pulldcaz[j]<<
5534 "pulldca0="<<pulldca[i]<<
5535 "pulldca1="<<pulldca[j]<<
5545 Float_t *quality = new Float_t[ncandidates];
5546 Int_t *indexes = new Int_t[ncandidates];
5548 for (Int_t i=0;i<ncandidates;i++){
5550 AliESDv0 *v0 = (AliESDv0*)tpcv0s->At(i);
5551 quality[i] = 1./(1.00001-v0->GetV0CosineOfPointingAngle()); //base point angle
5552 // quality[i] /= (0.5+v0->GetDist2());
5553 // quality[i] *= v0->GetChi2After(); //density factor
5555 Int_t index0 = v0->GetIndex(0);
5556 Int_t index1 = v0->GetIndex(1);
5557 //Bo: Double_t minpulldca = TMath::Min(2.+pulldca[v0->GetIndex(0)],(2.+pulldca[v0->GetIndex(1)]) ); //pull
5558 Double_t minpulldca = TMath::Min(2.+pulldca[index0],(2.+pulldca[index1]) );//Bo:
5562 AliTPCseed * track0 = (AliTPCseed*)array->At(index0);
5563 AliTPCseed * track1 = (AliTPCseed*)array->At(index1);
5564 if (track0->TPCrPID(0)>0.3&&track1->TPCrPID(0)>0.3&&v0->GetAnglep()[2]<0.15) quality[i]+=1000000; // gamma conversion candidate
5565 if (track0->TPCrPID(4)>0.9||track1->TPCrPID(4)>0.9&&minpulldca>4) quality[i]*=10; // lambda candidate candidate
5568 TMath::Sort(ncandidates,quality,indexes,kTRUE);
5571 for (Int_t i=0;i<ncandidates;i++){
5572 AliESDv0 * v0 = (AliESDv0*)tpcv0s->At(indexes[i]);
5574 Int_t index0 = v0->GetIndex(0);
5575 Int_t index1 = v0->GetIndex(1);
5576 AliTPCseed * track0 = (AliTPCseed*)array->At(index0);
5577 AliTPCseed * track1 = (AliTPCseed*)array->At(index1);
5578 if (!track0||!track1) {
5582 Bool_t accept =kTRUE; //default accept
5583 Int_t *v0indexes0 = track0->GetV0Indexes();
5584 Int_t *v0indexes1 = track1->GetV0Indexes();
5586 Int_t order0 = (v0indexes0[0]!=0) ? 1:0;
5587 Int_t order1 = (v0indexes1[0]!=0) ? 1:0;
5588 if (v0indexes0[1]!=0) order0 =2;
5589 if (v0indexes1[1]!=0) order1 =2;
5591 if (v0indexes0[2]!=0) {order0=3; accept=kFALSE;}
5592 if (v0indexes0[2]!=0) {order1=3; accept=kFALSE;}
5594 AliESDv0 * v02 = v0;
5596 //Bo: v0->SetOrder(0,order0);
5597 //Bo: v0->SetOrder(1,order1);
5598 //Bo: v0->SetOrder(1,order0+order1);
5599 v0->SetOnFlyStatus(kTRUE);
5600 Int_t index = esd->AddV0(v0);
5601 v02 = esd->GetV0(index);
5602 v0indexes0[order0]=index;
5603 v0indexes1[order1]=index;
5607 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
5608 if (AliTPCReconstructor::StreamLevel()>1) {
5609 Int_t lab0=track0->GetLabel();
5610 Int_t lab1=track1->GetLabel();
5619 "dca0="<<dca[index0]<<
5620 "dca1="<<dca[index1]<<
5624 "quality="<<quality[i]<<
5625 "pulldca0="<<pulldca[index0]<<
5626 "pulldca1="<<pulldca[index1]<<
5650 printf("TPC V0 finder : naccepted\t%d\tncandidates\t%d\tntracks\t%d\tnall\t%d\n",naccepted,ncandidates,ntracks,nall);
5654 Int_t AliTPCtrackerMI::RefitKink(AliTPCseed &mother, AliTPCseed &daughter, AliESDkink &knk)
5657 // refit kink towards to the vertex
5660 AliKink &kink=(AliKink &)knk;
5662 Int_t row0 = GetRowNumber(kink.GetR());
5663 FollowProlongation(mother,0);
5664 mother.Reset(kFALSE);
5666 FollowProlongation(daughter,row0);
5667 daughter.Reset(kFALSE);
5668 FollowBackProlongation(daughter,158);
5669 daughter.Reset(kFALSE);
5670 Int_t first = TMath::Max(row0-20,30);
5671 Int_t last = TMath::Min(row0+20,140);
5673 const Int_t kNdiv =5;
5674 AliTPCseed param0[kNdiv]; // parameters along the track
5675 AliTPCseed param1[kNdiv]; // parameters along the track
5676 AliKink kinks[kNdiv]; // corresponding kink parameters
5679 for (Int_t irow=0; irow<kNdiv;irow++){
5680 rows[irow] = first +((last-first)*irow)/(kNdiv-1);
5682 // store parameters along the track
5684 for (Int_t irow=0;irow<kNdiv;irow++){
5685 FollowBackProlongation(mother, rows[irow]);
5686 FollowProlongation(daughter,rows[kNdiv-1-irow]);
5687 param0[irow] = mother;
5688 param1[kNdiv-1-irow] = daughter;
5692 for (Int_t irow=0; irow<kNdiv-1;irow++){
5693 if (param0[irow].GetNumberOfClusters()<kNdiv||param1[irow].GetNumberOfClusters()<kNdiv) continue;
5694 kinks[irow].SetMother(param0[irow]);
5695 kinks[irow].SetDaughter(param1[irow]);
5696 kinks[irow].Update();
5699 // choose kink with best "quality"
5701 Double_t mindist = 10000;
5702 for (Int_t irow=0;irow<kNdiv;irow++){
5703 if (param0[irow].GetNumberOfClusters()<20||param1[irow].GetNumberOfClusters()<20) continue;
5704 if (TMath::Abs(kinks[irow].GetR())>240.) continue;
5705 if (TMath::Abs(kinks[irow].GetR())<100.) continue;
5707 Float_t normdist = TMath::Abs(param0[irow].GetX()-kinks[irow].GetR())*(0.1+kink.GetDistance());
5708 normdist/= (param0[irow].GetNumberOfClusters()+param1[irow].GetNumberOfClusters()+40.);
5709 if (normdist < mindist){
5715 if (index==-1) return 0;
5718 param0[index].Reset(kTRUE);
5719 FollowProlongation(param0[index],0);
5721 mother = param0[index];
5722 daughter = param1[index]; // daughter in vertex
5724 kink.SetMother(mother);
5725 kink.SetDaughter(daughter);
5727 kink.SetTPCRow0(GetRowNumber(kink.GetR()));
5728 kink.SetTPCncls(param0[index].GetNumberOfClusters(),0);
5729 kink.SetTPCncls(param1[index].GetNumberOfClusters(),1);
5730 kink.SetLabel(CookLabel(&mother,0.4, 0,kink.GetTPCRow0()),0);
5731 kink.SetLabel(CookLabel(&daughter,0.4, kink.GetTPCRow0(),160),1);
5732 mother.SetLabel(kink.GetLabel(0));
5733 daughter.SetLabel(kink.GetLabel(1));
5739 void AliTPCtrackerMI::UpdateKinkQualityM(AliTPCseed * seed){
5741 // update Kink quality information for mother after back propagation
5743 if (seed->GetKinkIndex(0)>=0) return;
5744 for (Int_t ikink=0;ikink<3;ikink++){
5745 Int_t index = seed->GetKinkIndex(ikink);
5746 if (index>=0) break;
5747 index = TMath::Abs(index)-1;
5748 AliESDkink * kink = fEvent->GetKink(index);
5749 kink->SetTPCDensity(-1,0,0);
5750 kink->SetTPCDensity(1,0,1);
5752 Int_t row0 = kink->GetTPCRow0() - 2 - Int_t( 0.5/ (0.05+kink->GetAngle(2)));
5753 if (row0<15) row0=15;
5755 Int_t row1 = kink->GetTPCRow0() + 2 + Int_t( 0.5/ (0.05+kink->GetAngle(2)));
5756 if (row1>145) row1=145;
5758 Int_t found,foundable,shared;
5759 seed->GetClusterStatistic(0,row0, found, foundable,shared,kFALSE);
5760 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),0,0);
5761 seed->GetClusterStatistic(row1,155, found, foundable,shared,kFALSE);
5762 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),0,1);
5767 void AliTPCtrackerMI::UpdateKinkQualityD(AliTPCseed * seed){
5769 // update Kink quality information for daughter after refit
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,1,0);
5778 kink->SetTPCDensity(-1,1,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),1,0);
5789 seed->GetClusterStatistic(row1,155, found, foundable,shared,kFALSE);
5790 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),1,1);
5796 Int_t AliTPCtrackerMI::CheckKinkPoint(AliTPCseed*seed,AliTPCseed &mother, AliTPCseed &daughter, AliESDkink &knk)
5799 // check kink point for given track
5800 // if return value=0 kink point not found
5801 // otherwise seed0 correspond to mother particle
5802 // seed1 correspond to daughter particle
5803 // kink parameter of kink point
5804 AliKink &kink=(AliKink &)knk;
5806 Int_t middlerow = (seed->GetFirstPoint()+seed->GetLastPoint())/2;
5807 Int_t first = seed->GetFirstPoint();
5808 Int_t last = seed->GetLastPoint();
5809 if (last-first<20) return 0; // shortest length - 2*30 = 60 pad-rows
5812 AliTPCseed *seed1 = ReSeed(seed,middlerow+20, kTRUE); //middle of chamber
5813 if (!seed1) return 0;
5814 FollowProlongation(*seed1,seed->GetLastPoint()-20);
5815 seed1->Reset(kTRUE);
5816 FollowProlongation(*seed1,158);
5817 seed1->Reset(kTRUE);
5818 last = seed1->GetLastPoint();
5820 AliTPCseed *seed0 = new AliTPCseed(*seed);
5821 seed0->Reset(kFALSE);
5824 AliTPCseed param0[20]; // parameters along the track
5825 AliTPCseed param1[20]; // parameters along the track
5826 AliKink kinks[20]; // corresponding kink parameters
5828 for (Int_t irow=0; irow<20;irow++){
5829 rows[irow] = first +((last-first)*irow)/19;
5831 // store parameters along the track
5833 for (Int_t irow=0;irow<20;irow++){
5834 FollowBackProlongation(*seed0, rows[irow]);
5835 FollowProlongation(*seed1,rows[19-irow]);
5836 param0[irow] = *seed0;
5837 param1[19-irow] = *seed1;
5841 for (Int_t irow=0; irow<19;irow++){
5842 kinks[irow].SetMother(param0[irow]);
5843 kinks[irow].SetDaughter(param1[irow]);
5844 kinks[irow].Update();
5847 // choose kink with biggest change of angle
5849 Double_t maxchange= 0;
5850 for (Int_t irow=1;irow<19;irow++){
5851 if (TMath::Abs(kinks[irow].GetR())>240.) continue;
5852 if (TMath::Abs(kinks[irow].GetR())<110.) continue;
5853 Float_t quality = TMath::Abs(kinks[irow].GetAngle(2))/(3.+TMath::Abs(kinks[irow].GetR()-param0[irow].GetX()));
5854 if ( quality > maxchange){
5855 maxchange = quality;
5862 if (index<0) return 0;
5864 Int_t row0 = GetRowNumber(kinks[index].GetR()); //row 0 estimate
5865 seed0 = new AliTPCseed(param0[index]);
5866 seed1 = new AliTPCseed(param1[index]);
5867 seed0->Reset(kFALSE);
5868 seed1->Reset(kFALSE);
5869 seed0->ResetCovariance(10.);
5870 seed1->ResetCovariance(10.);
5871 FollowProlongation(*seed0,0);
5872 FollowBackProlongation(*seed1,158);
5873 mother = *seed0; // backup mother at position 0
5874 seed0->Reset(kFALSE);
5875 seed1->Reset(kFALSE);
5876 seed0->ResetCovariance(10.);
5877 seed1->ResetCovariance(10.);
5879 first = TMath::Max(row0-20,0);
5880 last = TMath::Min(row0+20,158);
5882 for (Int_t irow=0; irow<20;irow++){
5883 rows[irow] = first +((last-first)*irow)/19;
5885 // store parameters along the track
5887 for (Int_t irow=0;irow<20;irow++){
5888 FollowBackProlongation(*seed0, rows[irow]);
5889 FollowProlongation(*seed1,rows[19-irow]);
5890 param0[irow] = *seed0;
5891 param1[19-irow] = *seed1;
5895 for (Int_t irow=0; irow<19;irow++){
5896 kinks[irow].SetMother(param0[irow]);
5897 kinks[irow].SetDaughter(param1[irow]);
5898 // param0[irow].Dump();
5899 //param1[irow].Dump();
5900 kinks[irow].Update();
5903 // choose kink with biggest change of angle
5906 for (Int_t irow=0;irow<20;irow++){
5907 if (TMath::Abs(kinks[irow].GetR())>250.) continue;
5908 if (TMath::Abs(kinks[irow].GetR())<90.) continue;
5909 Float_t quality = TMath::Abs(kinks[irow].GetAngle(2))/(3.+TMath::Abs(kinks[irow].GetR()-param0[irow].GetX()));
5910 if ( quality > maxchange){
5911 maxchange = quality;
5918 if (index==-1 || param0[index].GetNumberOfClusters()+param1[index].GetNumberOfClusters()<100){
5923 // Float_t anglesigma = TMath::Sqrt(param0[index].fC22+param0[index].fC33+param1[index].fC22+param1[index].fC33);
5925 kink.SetMother(param0[index]);
5926 kink.SetDaughter(param1[index]);
5928 row0 = GetRowNumber(kink.GetR());
5929 kink.SetTPCRow0(row0);
5930 kink.SetLabel(CookLabel(seed0,0.5,0,row0),0);
5931 kink.SetLabel(CookLabel(seed1,0.5,row0,158),1);
5932 kink.SetIndex(-10,0);
5933 kink.SetIndex(int(param0[index].GetNumberOfClusters()+param1[index].GetNumberOfClusters()),1);
5934 kink.SetTPCncls(param0[index].GetNumberOfClusters(),0);
5935 kink.SetTPCncls(param1[index].GetNumberOfClusters(),1);
5938 // new (&mother) AliTPCseed(param0[index]);
5939 daughter = param1[index];
5940 daughter.SetLabel(kink.GetLabel(1));
5941 param0[index].Reset(kTRUE);
5942 FollowProlongation(param0[index],0);
5943 mother = param0[index];
5944 mother.SetLabel(kink.GetLabel(0));
5954 AliTPCseed* AliTPCtrackerMI::ReSeed(AliTPCseed *t)
5957 // reseed - refit - track
5960 // Int_t last = fSectors->GetNRows()-1;
5962 if (fSectors == fOuterSec){
5963 first = TMath::Max(first, t->GetFirstPoint()-fInnerSec->GetNRows());
5967 first = t->GetFirstPoint();
5969 AliTPCseed * seed = MakeSeed(t,0.1,0.5,0.9);
5970 FollowBackProlongation(*t,fSectors->GetNRows()-1);
5972 FollowProlongation(*t,first);
5982 //_____________________________________________________________________________
5983 Int_t AliTPCtrackerMI::ReadSeeds(const TFile *inp) {
5984 //-----------------------------------------------------------------
5985 // This function reades track seeds.
5986 //-----------------------------------------------------------------
5987 TDirectory *savedir=gDirectory;
5989 TFile *in=(TFile*)inp;
5990 if (!in->IsOpen()) {
5991 cerr<<"AliTPCtrackerMI::ReadSeeds(): input file is not open !\n";
5996 TTree *seedTree=(TTree*)in->Get("Seeds");
5998 cerr<<"AliTPCtrackerMI::ReadSeeds(): ";
5999 cerr<<"can't get a tree with track seeds !\n";
6002 AliTPCtrack *seed=new AliTPCtrack;
6003 seedTree->SetBranchAddress("tracks",&seed);
6005 if (fSeeds==0) fSeeds=new TObjArray(15000);
6007 Int_t n=(Int_t)seedTree->GetEntries();
6008 for (Int_t i=0; i<n; i++) {
6009 seedTree->GetEvent(i);
6010 fSeeds->AddLast(new AliTPCseed(*seed/*,seed->GetAlpha()*/));
6019 Int_t AliTPCtrackerMI::Clusters2Tracks (AliESDEvent *esd)
6022 if (fSeeds) DeleteSeeds();
6025 if (!fSeeds) return 1;
6032 //_____________________________________________________________________________
6033 Int_t AliTPCtrackerMI::Clusters2Tracks() {
6034 //-----------------------------------------------------------------
6035 // This is a track finder.
6036 //-----------------------------------------------------------------
6037 TDirectory *savedir=gDirectory;
6041 fSeeds = Tracking();
6044 Info("Clusters2Tracks","Time for tracking: \t");timer.Print();timer.Start();
6046 //activate again some tracks
6047 for (Int_t i=0; i<fSeeds->GetEntriesFast(); i++) {
6048 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
6050 Int_t nc=t.GetNumberOfClusters();
6052 delete fSeeds->RemoveAt(i);
6056 if (pt->GetRemoval()==10) {
6057 if (pt->GetDensityFirst(20)>0.8 || pt->GetDensityFirst(30)>0.8 || pt->GetDensityFirst(40)>0.7)
6058 pt->Desactivate(10); // make track again active
6060 pt->Desactivate(20);
6061 delete fSeeds->RemoveAt(i);
6066 RemoveUsed2(fSeeds,0.85,0.85,0);
6067 if (AliTPCReconstructor::GetRecoParam()->GetDoKinks()) FindKinks(fSeeds,fEvent);
6068 //FindCurling(fSeeds, fEvent,0);
6069 if (AliTPCReconstructor::StreamLevel()>2) FindMultiMC(fSeeds, fEvent,0); // find multi found tracks
6070 RemoveUsed2(fSeeds,0.5,0.4,20);
6071 FindSplitted(fSeeds, fEvent,0); // find multi found tracks
6074 // // refit short tracks
6076 Int_t nseed=fSeeds->GetEntriesFast();
6079 for (Int_t i=0; i<nseed; i++) {
6080 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
6082 Int_t nc=t.GetNumberOfClusters();
6084 delete fSeeds->RemoveAt(i);
6087 CookLabel(pt,0.1); //For comparison only
6088 //if ((pt->IsActive() || (pt->fRemoval==10) )&& nc>50 &&pt->GetNumberOfClusters()>0.4*pt->fNFoundable){
6089 if ((pt->IsActive() || (pt->GetRemoval()==10) )){
6091 if (fDebug>0) cerr<<found<<'\r';
6095 delete fSeeds->RemoveAt(i);
6099 //RemoveOverlap(fSeeds,0.99,7,kTRUE);
6101 //RemoveUsed(fSeeds,0.9,0.9,6);
6103 nseed=fSeeds->GetEntriesFast();
6105 for (Int_t i=0; i<nseed; i++) {
6106 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
6108 Int_t nc=t.GetNumberOfClusters();
6110 delete fSeeds->RemoveAt(i);
6114 t.CookdEdx(0.02,0.6);
6115 // CheckKinkPoint(&t,0.05);
6116 //if ((pt->IsActive() || (pt->fRemoval==10) )&& nc>50 &&pt->GetNumberOfClusters()>0.4*pt->fNFoundable){
6117 if ((pt->IsActive() || (pt->GetRemoval()==10) )){
6125 delete fSeeds->RemoveAt(i);
6126 //AliTPCseed * seed1 = ReSeed(pt,0.05,0.5,1);
6128 // FollowProlongation(*seed1,0);
6129 // Int_t n = seed1->GetNumberOfClusters();
6130 // printf("fP4\t%f\t%f\n",seed1->GetC(),pt->GetC());
6131 // printf("fN\t%d\t%d\n", seed1->GetNumberOfClusters(),pt->GetNumberOfClusters());
6134 //AliTPCseed * seed2 = ReSeed(pt,0.95,0.5,0.05);
6138 SortTracks(fSeeds, 1);
6142 PrepareForBackProlongation(fSeeds,5.);
6143 PropagateBack(fSeeds);
6144 printf("Time for back propagation: \t");timer.Print();timer.Start();
6148 PrepareForProlongation(fSeeds,5.);
6149 PropagateForward2(fSeeds);
6151 printf("Time for FORWARD propagation: \t");timer.Print();timer.Start();
6152 // RemoveUsed(fSeeds,0.7,0.7,6);
6153 //RemoveOverlap(fSeeds,0.9,7,kTRUE);
6155 nseed=fSeeds->GetEntriesFast();
6157 for (Int_t i=0; i<nseed; i++) {
6158 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
6160 Int_t nc=t.GetNumberOfClusters();
6162 delete fSeeds->RemoveAt(i);
6165 t.CookdEdx(0.02,0.6);
6166 // CookLabel(pt,0.1); //For comparison only
6167 //if ((pt->IsActive() || (pt->fRemoval==10) )&& nc>50 &&pt->GetNumberOfClusters()>0.4*pt->fNFoundable){
6168 if ((pt->IsActive() || (pt->fRemoval==10) )){
6169 cerr<<found++<<'\r';
6172 delete fSeeds->RemoveAt(i);
6177 // fNTracks = found;
6179 Info("Clusters2Tracks","Time for overlap removal, track writing and dedx cooking: \t"); timer.Print();timer.Start();
6182 // cerr<<"Number of found tracks : "<<"\t"<<found<<endl;
6183 Info("Clusters2Tracks","Number of found tracks %d",found);
6185 // UnloadClusters();
6190 void AliTPCtrackerMI::Tracking(TObjArray * arr)
6193 // tracking of the seeds
6196 fSectors = fOuterSec;
6197 ParallelTracking(arr,150,63);
6198 fSectors = fOuterSec;
6199 ParallelTracking(arr,63,0);
6202 TObjArray * AliTPCtrackerMI::Tracking(Int_t seedtype, Int_t i1, Int_t i2, Float_t cuts[4], Float_t dy, Int_t dsec)
6207 TObjArray * arr = new TObjArray;
6209 fSectors = fOuterSec;
6212 for (Int_t sec=0;sec<fkNOS;sec++){
6213 if (seedtype==3) MakeSeeds3(arr,sec,i1,i2,cuts,dy, dsec);
6214 if (seedtype==4) MakeSeeds5(arr,sec,i1,i2,cuts,dy);
6215 if (seedtype==2) MakeSeeds2(arr,sec,i1,i2,cuts,dy);
6218 Info("Tracking","\nSeeding - %d\t%d\t%d\t%d\n",seedtype,i1,i2,arr->GetEntriesFast());
6230 TObjArray * AliTPCtrackerMI::Tracking()
6234 if (AliTPCReconstructor::GetRecoParam()->GetSpecialSeeding()) return TrackingSpecial();
6237 Int_t nup=fOuterSec->GetNRows()+fInnerSec->GetNRows();
6239 TObjArray * seeds = new TObjArray;
6248 Float_t fnumber = 3.0;
6249 Float_t fdensity = 3.0;
6254 for (Int_t delta = 0; delta<18; delta+=6){
6258 arr = Tracking(3,nup-1-delta,nup-1-delta-gap,cuts,-1,1);
6259 SumTracks(seeds,arr);
6260 SignClusters(seeds,fnumber,fdensity);
6262 for (Int_t i=2;i<6;i+=2){
6263 // seed high pt tracks
6266 arr = Tracking(3,nup-i-delta,nup-i-delta-gap,cuts,-1,0);
6267 SumTracks(seeds,arr);
6268 SignClusters(seeds,fnumber,fdensity);
6273 // RemoveUsed(seeds,0.9,0.9,1);
6274 // UnsignClusters();
6275 // SignClusters(seeds,fnumber,fdensity);
6279 for (Int_t delta = 20; delta<120; delta+=10){
6281 // seed high pt tracks
6285 arr = Tracking(3,nup-delta,nup-delta-gap,cuts,-1);
6286 SumTracks(seeds,arr);
6287 SignClusters(seeds,fnumber,fdensity);
6292 arr = Tracking(3,nup-delta-5,nup-delta-5-gap,cuts,-1);
6293 SumTracks(seeds,arr);
6294 SignClusters(seeds,fnumber,fdensity);
6305 Info("Tracking()","\n\nPrimary seeding\t%d\n\n",seeds->GetEntriesFast());
6309 // RemoveUsed(seeds,0.75,0.75,1);
6311 //SignClusters(seeds,fnumber,fdensity);
6320 arr = Tracking(4,nup-1,nup-1-gap,cuts,-1);
6321 SumTracks(seeds,arr);
6322 SignClusters(seeds,fnumber,fdensity);
6324 arr = Tracking(4,nup-2,nup-2-gap,cuts,-1);
6325 SumTracks(seeds,arr);
6326 SignClusters(seeds,fnumber,fdensity);
6328 arr = Tracking(4,nup-3,nup-3-gap,cuts,-1);
6329 SumTracks(seeds,arr);
6330 SignClusters(seeds,fnumber,fdensity);
6334 for (Int_t delta = 3; delta<30; delta+=5){
6340 arr = Tracking(4,nup-1-delta,nup-1-delta-gap,cuts,-1);
6341 SumTracks(seeds,arr);
6342 SignClusters(seeds,fnumber,fdensity);
6344 arr = Tracking(4,nup-3-delta,nup-5-delta-gap,cuts,4);
6345 SumTracks(seeds,arr);
6346 SignClusters(seeds,fnumber,fdensity);
6357 Int_t fLastSeedRowSec=AliTPCReconstructor::GetRecoParam()->GetLastSeedRowSec();
6360 for (Int_t delta = 30; delta<fLastSeedRowSec; delta+=10){
6366 arr = Tracking(4,nup-1-delta,nup-1-delta-gap,cuts,-1);
6367 SumTracks(seeds,arr);
6368 SignClusters(seeds,fnumber,fdensity);
6370 arr = Tracking(4,nup-5-delta,nup-5-delta-gap,cuts,5 );
6371 SumTracks(seeds,arr);
6372 SignClusters(seeds,fnumber,fdensity);
6376 Info("Tracking()","\n\nSecondary seeding\t%d\n\n",seeds->GetEntriesFast());
6387 TObjArray * AliTPCtrackerMI::TrackingSpecial()
6390 // seeding adjusted for laser and cosmic tests - short tracks with big inclination angle
6391 // no primary vertex seeding tried
6395 Int_t nup=fOuterSec->GetNRows()+fInnerSec->GetNRows();
6397 TObjArray * seeds = new TObjArray;
6402 Float_t fnumber = 3.0;
6403 Float_t fdensity = 3.0;
6406 cuts[0] = AliTPCReconstructor::GetRecoParam()->GetMaxC(); // max curvature
6407 cuts[1] = 3.5; // max tan(phi) angle for seeding
6408 cuts[2] = 3.; // not used (cut on z primary vertex)
6409 cuts[3] = 3.5; // max tan(theta) angle for seeding
6411 for (Int_t delta = 0; nup-delta-gap-1>0; delta+=3){
6413 arr = Tracking(4,nup-1-delta,nup-1-delta-gap,cuts,-1);
6414 SumTracks(seeds,arr);
6415 SignClusters(seeds,fnumber,fdensity);
6419 Info("Tracking()","\n\nSecondary seeding\t%d\n\n",seeds->GetEntriesFast());
6430 void AliTPCtrackerMI::SumTracks(TObjArray *arr1,TObjArray *&arr2) const
6433 //sum tracks to common container
6434 //remove suspicious tracks
6435 Int_t nseed = arr2->GetEntriesFast();
6436 for (Int_t i=0;i<nseed;i++){
6437 AliTPCseed *pt=(AliTPCseed*)arr2->UncheckedAt(i);
6440 // remove tracks with too big curvature
6442 if (TMath::Abs(pt->GetC())>AliTPCReconstructor::GetRecoParam()->GetMaxC()){
6443 delete arr2->RemoveAt(i);
6446 // REMOVE VERY SHORT TRACKS
6447 if (pt->GetNumberOfClusters()<20){
6448 delete arr2->RemoveAt(i);
6451 // NORMAL ACTIVE TRACK
6452 if (pt->IsActive()){
6453 arr1->AddLast(arr2->RemoveAt(i));
6456 //remove not usable tracks
6457 if (pt->GetRemoval()!=10){
6458 delete arr2->RemoveAt(i);
6462 // ENABLE ONLY ENOUGH GOOD STOPPED TRACKS
6463 if (pt->GetDensityFirst(20)>0.8 || pt->GetDensityFirst(30)>0.8 || pt->GetDensityFirst(40)>0.7)
6464 arr1->AddLast(arr2->RemoveAt(i));
6466 delete arr2->RemoveAt(i);
6470 delete arr2; arr2 = 0;
6475 void AliTPCtrackerMI::ParallelTracking(TObjArray * arr, Int_t rfirst, Int_t rlast)
6478 // try to track in parralel
6480 Int_t nseed=arr->GetEntriesFast();
6481 //prepare seeds for tracking
6482 for (Int_t i=0; i<nseed; i++) {
6483 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i), &t=*pt;
6485 if (!t.IsActive()) continue;
6486 // follow prolongation to the first layer
6487 if ( (fSectors ==fInnerSec) || (t.GetFirstPoint()-fParam->GetNRowLow()>rfirst+1) )
6488 FollowProlongation(t, rfirst+1);
6493 for (Int_t nr=rfirst; nr>=rlast; nr--){
6494 if (nr<fInnerSec->GetNRows())
6495 fSectors = fInnerSec;
6497 fSectors = fOuterSec;
6498 // make indexes with the cluster tracks for given
6500 // find nearest cluster
6501 for (Int_t i=0; i<nseed; i++) {
6502 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i), &t=*pt;
6504 if (nr==80) pt->UpdateReference();
6505 if (!pt->IsActive()) continue;
6506 // if ( (fSectors ==fOuterSec) && (pt->fFirstPoint-fParam->GetNRowLow())<nr) continue;
6507 if (pt->GetRelativeSector()>17) {
6510 UpdateClusters(t,nr);
6512 // prolonagate to the nearest cluster - if founded
6513 for (Int_t i=0; i<nseed; i++) {
6514 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
6516 if (!pt->IsActive()) continue;
6517 // if ((fSectors ==fOuterSec) && (pt->fFirstPoint-fParam->GetNRowLow())<nr) continue;
6518 if (pt->GetRelativeSector()>17) {
6521 FollowToNextCluster(*pt,nr);
6526 void AliTPCtrackerMI::PrepareForBackProlongation(TObjArray * arr,Float_t fac) const
6530 // if we use TPC track itself we have to "update" covariance
6532 Int_t nseed= arr->GetEntriesFast();
6533 for (Int_t i=0;i<nseed;i++){
6534 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
6538 //rotate to current local system at first accepted point
6539 Int_t index = pt->GetClusterIndex2(pt->GetFirstPoint());
6540 Int_t sec = (index&0xff000000)>>24;
6542 Float_t angle1 = fInnerSec->GetAlpha()*sec+fInnerSec->GetAlphaShift();
6543 if (angle1>TMath::Pi())
6544 angle1-=2.*TMath::Pi();
6545 Float_t angle2 = pt->GetAlpha();
6547 if (TMath::Abs(angle1-angle2)>0.001){
6548 pt->Rotate(angle1-angle2);
6549 //angle2 = pt->GetAlpha();
6550 //pt->fRelativeSector = pt->GetAlpha()/fInnerSec->GetAlpha();
6551 //if (pt->GetAlpha()<0)
6552 // pt->fRelativeSector+=18;
6553 //sec = pt->fRelativeSector;
6562 void AliTPCtrackerMI::PrepareForProlongation(TObjArray * arr, Float_t fac) const
6566 // if we use TPC track itself we have to "update" covariance
6568 Int_t nseed= arr->GetEntriesFast();
6569 for (Int_t i=0;i<nseed;i++){
6570 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
6573 pt->SetFirstPoint(pt->GetLastPoint());
6581 Int_t AliTPCtrackerMI::PropagateBack(TObjArray * arr)
6584 // make back propagation
6586 Int_t nseed= arr->GetEntriesFast();
6587 for (Int_t i=0;i<nseed;i++){
6588 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
6589 if (pt&& pt->GetKinkIndex(0)<=0) {
6590 //AliTPCseed *pt2 = new AliTPCseed(*pt);
6591 fSectors = fInnerSec;
6592 //FollowBackProlongation(*pt,fInnerSec->GetNRows()-1);
6593 //fSectors = fOuterSec;
6594 FollowBackProlongation(*pt,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1);
6595 //if (pt->GetNumberOfClusters()<(pt->fEsd->GetTPCclusters(0)) ){
6596 // Error("PropagateBack","Not prolonged track %d",pt->GetLabel());
6597 // FollowBackProlongation(*pt2,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1);
6600 if (pt&& pt->GetKinkIndex(0)>0) {
6601 AliESDkink * kink = fEvent->GetKink(pt->GetKinkIndex(0)-1);
6602 pt->SetFirstPoint(kink->GetTPCRow0());
6603 fSectors = fInnerSec;
6604 FollowBackProlongation(*pt,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1);
6612 Int_t AliTPCtrackerMI::PropagateForward2(TObjArray * arr)
6615 // make forward propagation
6617 Int_t nseed= arr->GetEntriesFast();
6619 for (Int_t i=0;i<nseed;i++){
6620 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
6622 FollowProlongation(*pt,0);
6631 Int_t AliTPCtrackerMI::PropagateForward()
6634 // propagate track forward
6636 Int_t nseed = fSeeds->GetEntriesFast();
6637 for (Int_t i=0;i<nseed;i++){
6638 AliTPCseed *pt = (AliTPCseed*)fSeeds->UncheckedAt(i);
6640 AliTPCseed &t = *pt;
6641 Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
6642 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
6643 if (alpha < 0. ) alpha += 2.*TMath::Pi();
6644 t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN);
6648 fSectors = fOuterSec;
6649 ParallelTracking(fSeeds,fOuterSec->GetNRows()+fInnerSec->GetNRows()-1,fInnerSec->GetNRows());
6650 fSectors = fInnerSec;
6651 ParallelTracking(fSeeds,fInnerSec->GetNRows()-1,0);
6660 Int_t AliTPCtrackerMI::PropagateBack(AliTPCseed * pt, Int_t row0, Int_t row1)
6663 // make back propagation, in between row0 and row1
6667 fSectors = fInnerSec;
6670 if (row1<fSectors->GetNRows())
6673 r1 = fSectors->GetNRows()-1;
6675 if (row0<fSectors->GetNRows()&& r1>0 )
6676 FollowBackProlongation(*pt,r1);
6677 if (row1<=fSectors->GetNRows())
6680 r1 = row1 - fSectors->GetNRows();
6681 if (r1<=0) return 0;
6682 if (r1>=fOuterSec->GetNRows()) return 0;
6683 fSectors = fOuterSec;
6684 return FollowBackProlongation(*pt,r1);
6692 void AliTPCtrackerMI::GetShape(AliTPCseed * seed, Int_t row)
6696 AliTPCClusterParam * clparam = AliTPCcalibDB::Instance()->GetClusterParam();
6697 Float_t zdrift = TMath::Abs((fParam->GetZLength(0)-TMath::Abs(seed->GetZ())));
6698 Int_t type = (seed->GetSector() < fParam->GetNSector()/2) ? 0: (row>126) ? 1:2;
6699 Double_t angulary = seed->GetSnp();
6700 angulary = angulary*angulary/(1.-angulary*angulary);
6701 Double_t angularz = seed->GetTgl()*seed->GetTgl()*(1.+angulary);
6703 Double_t sigmay = clparam->GetRMS0(0,type,zdrift,TMath::Sqrt(TMath::Abs(angulary)));
6704 Double_t sigmaz = clparam->GetRMS0(1,type,zdrift,TMath::Sqrt(TMath::Abs(angularz)));
6705 seed->SetCurrentSigmaY2(sigmay*sigmay);
6706 seed->SetCurrentSigmaZ2(sigmaz*sigmaz);
6707 // Float_t sd2 = TMath::Abs((fParam->GetZLength(0)-TMath::Abs(seed->GetZ())))*fParam->GetDiffL()*fParam->GetDiffL();
6708 // // Float_t padlength = fParam->GetPadPitchLength(seed->fSector);
6709 // Float_t padlength = GetPadPitchLength(row);
6711 // Float_t sresy = (seed->GetSector() < fParam->GetNSector()/2) ? 0.2 :0.3;
6712 // seed->SetCurrentSigmaY2(sd2+padlength*padlength*angulary/12.+sresy*sresy);
6714 // Float_t sresz = fParam->GetZSigma();
6715 // seed->SetCurrentSigmaZ2(sd2+padlength*padlength*angularz*angularz*(1+angulary)/12.+sresz*sresz);
6717 Float_t wy = GetSigmaY(seed);
6718 Float_t wz = GetSigmaZ(seed);
6721 if (TMath::Abs(wy/seed->fCurrentSigmaY2-1)>0.0001 || TMath::Abs(wz/seed->fCurrentSigmaZ2-1)>0.0001 ){
6722 printf("problem\n");
6729 //__________________________________________________________________________
6730 void AliTPCtrackerMI::CookLabel(AliKalmanTrack *tk, Float_t wrong) const {
6731 //--------------------------------------------------------------------
6732 //This function "cooks" a track label. If label<0, this track is fake.
6733 //--------------------------------------------------------------------
6734 AliTPCseed * t = dynamic_cast<AliTPCseed*>(tk);
6736 printf("%s:%d wrong type \n",(char*)__FILE__,__LINE__);
6740 Int_t noc=t->GetNumberOfClusters();
6742 //printf("\nnot founded prolongation\n\n\n");
6748 AliTPCclusterMI *clusters[160];
6750 for (Int_t i=0;i<160;i++) {
6757 for (i=0; i<160 && current<noc; i++) {
6759 Int_t index=t->GetClusterIndex2(i);
6760 if (index<=0) continue;
6761 if (index&0x8000) continue;
6763 //clusters[current]=GetClusterMI(index);
6764 if (t->GetClusterPointer(i)){
6765 clusters[current]=t->GetClusterPointer(i);
6771 Int_t lab=123456789;
6772 for (i=0; i<noc; i++) {
6773 AliTPCclusterMI *c=clusters[i];
6775 lab=TMath::Abs(c->GetLabel(0));
6777 for (j=0; j<noc; j++) if (lb[j]==lab || mx[j]==0) break;
6783 for (i=0; i<noc; i++) if (mx[i]>max) {max=mx[i]; lab=lb[i];}
6785 for (i=0; i<noc; i++) {
6786 AliTPCclusterMI *c=clusters[i];
6788 if (TMath::Abs(c->GetLabel(1)) == lab ||
6789 TMath::Abs(c->GetLabel(2)) == lab ) max++;
6792 if ((1.- Float_t(max)/noc) > wrong) lab=-lab;
6795 Int_t tail=Int_t(0.10*noc);
6798 for (i=1; i<=160&&ind<tail; i++) {
6799 // AliTPCclusterMI *c=clusters[noc-i];
6800 AliTPCclusterMI *c=clusters[i];
6802 if (lab == TMath::Abs(c->GetLabel(0)) ||
6803 lab == TMath::Abs(c->GetLabel(1)) ||
6804 lab == TMath::Abs(c->GetLabel(2))) max++;
6807 if (max < Int_t(0.5*tail)) lab=-lab;
6814 //delete[] clusters;
6818 //__________________________________________________________________________
6819 Int_t AliTPCtrackerMI::CookLabel(AliTPCseed *t, Float_t wrong,Int_t first, Int_t last) const {
6820 //--------------------------------------------------------------------
6821 //This function "cooks" a track label. If label<0, this track is fake.
6822 //--------------------------------------------------------------------
6823 Int_t noc=t->GetNumberOfClusters();
6825 //printf("\nnot founded prolongation\n\n\n");
6831 AliTPCclusterMI *clusters[160];
6833 for (Int_t i=0;i<160;i++) {
6840 for (i=0; i<160 && current<noc; i++) {
6841 if (i<first) continue;
6842 if (i>last) continue;
6843 Int_t index=t->GetClusterIndex2(i);
6844 if (index<=0) continue;
6845 if (index&0x8000) continue;
6847 //clusters[current]=GetClusterMI(index);
6848 if (t->GetClusterPointer(i)){
6849 clusters[current]=t->GetClusterPointer(i);
6854 if (noc<5) return -1;
6855 Int_t lab=123456789;
6856 for (i=0; i<noc; i++) {
6857 AliTPCclusterMI *c=clusters[i];
6859 lab=TMath::Abs(c->GetLabel(0));
6861 for (j=0; j<noc; j++) if (lb[j]==lab || mx[j]==0) break;
6867 for (i=0; i<noc; i++) if (mx[i]>max) {max=mx[i]; lab=lb[i];}
6869 for (i=0; i<noc; i++) {
6870 AliTPCclusterMI *c=clusters[i];
6872 if (TMath::Abs(c->GetLabel(1)) == lab ||
6873 TMath::Abs(c->GetLabel(2)) == lab ) max++;
6876 if ((1.- Float_t(max)/noc) > wrong) lab=-lab;
6879 Int_t tail=Int_t(0.10*noc);
6882 for (i=1; i<=160&&ind<tail; i++) {
6883 // AliTPCclusterMI *c=clusters[noc-i];
6884 AliTPCclusterMI *c=clusters[i];
6886 if (lab == TMath::Abs(c->GetLabel(0)) ||
6887 lab == TMath::Abs(c->GetLabel(1)) ||
6888 lab == TMath::Abs(c->GetLabel(2))) max++;
6891 if (max < Int_t(0.5*tail)) lab=-lab;
6894 // t->SetLabel(lab);
6898 //delete[] clusters;
6902 Int_t AliTPCtrackerMI::GetRowNumber(Double_t x[3]) const
6904 //return pad row number for given x vector
6905 Float_t phi = TMath::ATan2(x[1],x[0]);
6906 if(phi<0) phi=2.*TMath::Pi()+phi;
6907 // Get the local angle in the sector philoc
6908 const Float_t kRaddeg = 180/3.14159265358979312;
6909 Float_t phiangle = (Int_t (phi*kRaddeg/20.) + 0.5)*20./kRaddeg;
6910 Double_t localx = x[0]*TMath::Cos(phiangle)-x[1]*TMath::Sin(phiangle);
6911 return GetRowNumber(localx);
6916 void AliTPCtrackerMI::MakeBitmaps(AliTPCseed *t)
6918 //-----------------------------------------------------------------------
6919 // Fill the cluster and sharing bitmaps of the track
6920 //-----------------------------------------------------------------------
6922 Int_t firstpoint = 0;
6923 Int_t lastpoint = 159;
6924 AliTPCTrackerPoint *point;
6926 for (int iter=firstpoint; iter<lastpoint; iter++) {
6927 point = t->GetTrackPoint(iter);
6929 t->SetClusterMapBit(iter, kTRUE);
6930 if (point->IsShared())
6931 t->SetSharedMapBit(iter,kTRUE);
6933 t->SetSharedMapBit(iter, kFALSE);
6936 t->SetClusterMapBit(iter, kFALSE);
6937 t->SetSharedMapBit(iter, kFALSE);
6944 void AliTPCtrackerMI::AddCovariance(AliTPCseed * seed){
6946 // Adding systematic error
6947 // !!!! the systematic error for element 4 is in 1/cm not in pt
6949 const Double_t *param = AliTPCReconstructor::GetRecoParam()->GetSystematicError();
6951 for (Int_t i=0;i<15;i++) covar[i]=0;
6957 covar[0] = param[0]*param[0];
6958 covar[2] = param[1]*param[1];
6959 covar[5] = param[2]*param[2];
6960 covar[9] = param[3]*param[3];
6961 Double_t facC = AliTracker::GetBz()*kB2C;
6962 covar[14]= param[4]*param[4]*facC*facC;
6963 seed->AddCovariance(covar);