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);
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 - dump information about
100 // 2.a) cluster error estimate
101 // 3.a) cluster shape estimate
104 // Debug streamer levels:
106 //-------------------------------------------------------
111 #include "Riostream.h"
112 #include <TClonesArray.h>
114 #include <TObjArray.h>
117 #include "AliComplexCluster.h"
118 #include "AliESDEvent.h"
119 #include "AliESDtrack.h"
120 #include "AliESDVertex.h"
123 #include "AliHelix.h"
124 #include "AliRunLoader.h"
125 #include "AliTPCClustersRow.h"
126 #include "AliTPCParam.h"
127 #include "AliTPCReconstructor.h"
128 #include "AliTPCpolyTrack.h"
129 #include "AliTPCreco.h"
130 #include "AliTPCseed.h"
132 #include "AliTPCtrackerSector.h"
133 #include "AliTPCtrackerMI.h"
134 #include "TStopwatch.h"
135 #include "AliTPCReconstructor.h"
136 #include "AliAlignObj.h"
137 #include "AliTrackPointArray.h"
139 #include "AliTPCcalibDB.h"
140 #include "AliTPCTransform.h"
141 #include "AliTPCClusterParam.h"
145 ClassImp(AliTPCtrackerMI)
149 class AliTPCFastMath {
152 static Double_t FastAsin(Double_t x);
154 static Double_t fgFastAsin[20000]; //lookup table for fast asin computation
157 Double_t AliTPCFastMath::fgFastAsin[20000];
158 AliTPCFastMath gAliTPCFastMath; // needed to fill the LUT
160 AliTPCFastMath::AliTPCFastMath(){
162 // initialized lookup table;
163 for (Int_t i=0;i<10000;i++){
164 fgFastAsin[2*i] = TMath::ASin(i/10000.);
165 fgFastAsin[2*i+1] = (TMath::ASin((i+1)/10000.)-fgFastAsin[2*i]);
169 Double_t AliTPCFastMath::FastAsin(Double_t x){
171 // return asin using lookup table
173 Int_t index = int(x*10000);
174 return fgFastAsin[2*index]+(x*10000.-index)*fgFastAsin[2*index+1];
177 Int_t index = int(x*10000);
178 return -(fgFastAsin[2*index]+(x*10000.-index)*fgFastAsin[2*index+1]);
180 //__________________________________________________________________
181 AliTPCtrackerMI::AliTPCtrackerMI()
203 // default constructor
206 //_____________________________________________________________________
210 Int_t AliTPCtrackerMI::UpdateTrack(AliTPCseed * track, Int_t accept){
212 //update track information using current cluster - track->fCurrentCluster
215 AliTPCclusterMI* c =track->GetCurrentCluster();
216 if (accept>0) track->SetCurrentClusterIndex1(track->GetCurrentClusterIndex1() | 0x8000); //sign not accepted clusters
218 UInt_t i = track->GetCurrentClusterIndex1();
220 Int_t sec=(i&0xff000000)>>24;
221 //Int_t row = (i&0x00ff0000)>>16;
222 track->SetRow((i&0x00ff0000)>>16);
223 track->SetSector(sec);
224 // Int_t index = i&0xFFFF;
225 if (sec>=fkParam->GetNInnerSector()) track->SetRow(track->GetRow()+fkParam->GetNRowLow());
226 track->SetClusterIndex2(track->GetRow(), i);
227 //track->fFirstPoint = row;
228 //if ( track->fLastPoint<row) track->fLastPoint =row;
229 // if (track->fRow<0 || track->fRow>160) {
230 // printf("problem\n");
232 if (track->GetFirstPoint()>track->GetRow())
233 track->SetFirstPoint(track->GetRow());
234 if (track->GetLastPoint()<track->GetRow())
235 track->SetLastPoint(track->GetRow());
238 track->SetClusterPointer(track->GetRow(),c);
241 Double_t angle2 = track->GetSnp()*track->GetSnp();
243 //SET NEW Track Point
245 if (angle2<1) //PH sometimes angle2 is very big. To be investigated...
247 angle2 = TMath::Sqrt(angle2/(1-angle2));
248 AliTPCTrackerPoint &point =*(track->GetTrackPoint(track->GetRow()));
250 point.SetSigmaY(c->GetSigmaY2()/track->GetCurrentSigmaY2());
251 point.SetSigmaZ(c->GetSigmaZ2()/track->GetCurrentSigmaZ2());
252 point.SetErrY(sqrt(track->GetErrorY2()));
253 point.SetErrZ(sqrt(track->GetErrorZ2()));
255 point.SetX(track->GetX());
256 point.SetY(track->GetY());
257 point.SetZ(track->GetZ());
258 point.SetAngleY(angle2);
259 point.SetAngleZ(track->GetTgl());
260 if (point.IsShared()){
261 track->SetErrorY2(track->GetErrorY2()*4);
262 track->SetErrorZ2(track->GetErrorZ2()*4);
266 Double_t chi2 = track->GetPredictedChi2(track->GetCurrentCluster());
268 // track->SetErrorY2(track->GetErrorY2()*1.3);
269 // track->SetErrorY2(track->GetErrorY2()+0.01);
270 // track->SetErrorZ2(track->GetErrorZ2()*1.3);
271 // track->SetErrorZ2(track->GetErrorZ2()+0.005);
273 if (accept>0) return 0;
274 if (track->GetNumberOfClusters()%20==0){
275 // if (track->fHelixIn){
276 // TClonesArray & larr = *(track->fHelixIn);
277 // Int_t ihelix = larr.GetEntriesFast();
278 // new(larr[ihelix]) AliHelix(*track) ;
281 track->SetNoCluster(0);
282 return track->Update(c,chi2,i);
287 Int_t AliTPCtrackerMI::AcceptCluster(AliTPCseed * seed, AliTPCclusterMI * cluster)
290 // decide according desired precision to accept given
291 // cluster for tracking
293 seed->GetProlongation(cluster->GetX(),yt,zt);
294 Double_t sy2=ErrY2(seed,cluster);
295 Double_t sz2=ErrZ2(seed,cluster);
297 Double_t sdistancey2 = sy2+seed->GetSigmaY2();
298 Double_t sdistancez2 = sz2+seed->GetSigmaZ2();
299 Double_t dy=seed->GetCurrentCluster()->GetY()-yt;
300 Double_t dz=seed->GetCurrentCluster()->GetZ()-zt;
301 Double_t rdistancey2 = (seed->GetCurrentCluster()->GetY()-yt)*
302 (seed->GetCurrentCluster()->GetY()-yt)/sdistancey2;
303 Double_t rdistancez2 = (seed->GetCurrentCluster()->GetZ()-zt)*
304 (seed->GetCurrentCluster()->GetZ()-zt)/sdistancez2;
306 Double_t rdistance2 = rdistancey2+rdistancez2;
309 if (AliTPCReconstructor::StreamLevel()>2 && seed->GetNumberOfClusters()>20) {
310 Float_t rmsy2 = seed->GetCurrentSigmaY2();
311 Float_t rmsz2 = seed->GetCurrentSigmaZ2();
312 Float_t rmsy2p30 = seed->GetCMeanSigmaY2p30();
313 Float_t rmsz2p30 = seed->GetCMeanSigmaZ2p30();
314 Float_t rmsy2p30R = seed->GetCMeanSigmaY2p30R();
315 Float_t rmsz2p30R = seed->GetCMeanSigmaZ2p30R();
316 AliExternalTrackParam param(*seed);
317 static TVectorD gcl(3),gtr(3);
319 param.GetXYZ(gcl.GetMatrixArray());
320 cluster->GetGlobalXYZ(gclf);
321 gcl[0]=gclf[0]; gcl[1]=gclf[1]; gcl[2]=gclf[2];
324 if (AliTPCReconstructor::StreamLevel()>2) {
325 (*fDebugStreamer)<<"ErrParam"<<
338 "rmsy2p30="<<rmsy2p30<<
339 "rmsz2p30="<<rmsz2p30<<
340 "rmsy2p30R="<<rmsy2p30R<<
341 "rmsz2p30R="<<rmsz2p30R<<
342 // normalize distance -
343 "rdisty="<<rdistancey2<<
344 "rdistz="<<rdistancez2<<
345 "rdist="<<rdistance2<< //
349 //return 0; // temporary
350 if (rdistance2>32) return 3;
353 if ((rdistancey2>9. || rdistancez2>9.) && cluster->GetType()==0)
354 return 2; //suspisiouce - will be changed
356 if ((rdistancey2>6.25 || rdistancez2>6.25) && cluster->GetType()>0)
357 // strict cut on overlaped cluster
358 return 2; //suspisiouce - will be changed
360 if ( (rdistancey2>1. || rdistancez2>6.25 )
361 && cluster->GetType()<0){
362 seed->SetNFoundable(seed->GetNFoundable()-1);
372 //_____________________________________________________________________________
373 AliTPCtrackerMI::AliTPCtrackerMI(const AliTPCParam *par):
375 fkNIS(par->GetNInnerSector()/2),
377 fkNOS(par->GetNOuterSector()/2),
394 //---------------------------------------------------------------------
395 // The main TPC tracker constructor
396 //---------------------------------------------------------------------
397 fInnerSec=new AliTPCtrackerSector[fkNIS];
398 fOuterSec=new AliTPCtrackerSector[fkNOS];
401 for (i=0; i<fkNIS; i++) fInnerSec[i].Setup(par,0);
402 for (i=0; i<fkNOS; i++) fOuterSec[i].Setup(par,1);
405 Int_t nrowlow = par->GetNRowLow();
406 Int_t nrowup = par->GetNRowUp();
409 for (i=0;i<nrowlow;i++){
410 fXRow[i] = par->GetPadRowRadiiLow(i);
411 fPadLength[i]= par->GetPadPitchLength(0,i);
412 fYMax[i] = fXRow[i]*TMath::Tan(0.5*par->GetInnerAngle());
416 for (i=0;i<nrowup;i++){
417 fXRow[i+nrowlow] = par->GetPadRowRadiiUp(i);
418 fPadLength[i+nrowlow] = par->GetPadPitchLength(60,i);
419 fYMax[i+nrowlow] = fXRow[i+nrowlow]*TMath::Tan(0.5*par->GetOuterAngle());
422 if (AliTPCReconstructor::StreamLevel()>0) {
423 fDebugStreamer = new TTreeSRedirector("TPCdebug.root");
426 //________________________________________________________________________
427 AliTPCtrackerMI::AliTPCtrackerMI(const AliTPCtrackerMI &t):
448 //------------------------------------
449 // dummy copy constructor
450 //------------------------------------------------------------------
453 AliTPCtrackerMI & AliTPCtrackerMI::operator=(const AliTPCtrackerMI& /*r*/)
455 //------------------------------
457 //--------------------------------------------------------------
460 //_____________________________________________________________________________
461 AliTPCtrackerMI::~AliTPCtrackerMI() {
462 //------------------------------------------------------------------
463 // TPC tracker destructor
464 //------------------------------------------------------------------
471 if (fDebugStreamer) delete fDebugStreamer;
475 void AliTPCtrackerMI::FillESD(const TObjArray* arr)
479 //fill esds using updated tracks
481 // write tracks to the event
482 // store index of the track
483 Int_t nseed=arr->GetEntriesFast();
484 //FindKinks(arr,fEvent);
485 for (Int_t i=0; i<nseed; i++) {
486 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
490 if (AliTPCReconstructor::StreamLevel()>1) {
491 (*fDebugStreamer)<<"Track0"<<
495 // pt->PropagateTo(fkParam->GetInnerRadiusLow());
496 if (pt->GetKinkIndex(0)<=0){ //don't propagate daughter tracks
497 pt->PropagateTo(fkParam->GetInnerRadiusLow());
500 if (( pt->GetPoints()[2]- pt->GetPoints()[0])>5 && pt->GetPoints()[3]>0.8){
502 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
503 iotrack.SetTPCPoints(pt->GetPoints());
504 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
505 iotrack.SetV0Indexes(pt->GetV0Indexes());
506 // iotrack.SetTPCpid(pt->fTPCr);
507 //iotrack.SetTPCindex(i);
508 fEvent->AddTrack(&iotrack);
512 if ( (pt->GetNumberOfClusters()>70)&& (Float_t(pt->GetNumberOfClusters())/Float_t(pt->GetNFoundable()))>0.55) {
514 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
515 iotrack.SetTPCPoints(pt->GetPoints());
516 //iotrack.SetTPCindex(i);
517 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
518 iotrack.SetV0Indexes(pt->GetV0Indexes());
519 // iotrack.SetTPCpid(pt->fTPCr);
520 fEvent->AddTrack(&iotrack);
524 // short tracks - maybe decays
526 if ( (pt->GetNumberOfClusters()>30) && (Float_t(pt->GetNumberOfClusters())/Float_t(pt->GetNFoundable()))>0.70) {
527 Int_t found,foundable,shared;
528 pt->GetClusterStatistic(0,60,found, foundable,shared,kFALSE);
529 if ( (found>20) && (pt->GetNShared()/float(pt->GetNumberOfClusters())<0.2)){
531 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
532 //iotrack.SetTPCindex(i);
533 iotrack.SetTPCPoints(pt->GetPoints());
534 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
535 iotrack.SetV0Indexes(pt->GetV0Indexes());
536 //iotrack.SetTPCpid(pt->fTPCr);
537 fEvent->AddTrack(&iotrack);
542 if ( (pt->GetNumberOfClusters()>20) && (Float_t(pt->GetNumberOfClusters())/Float_t(pt->GetNFoundable()))>0.8) {
543 Int_t found,foundable,shared;
544 pt->GetClusterStatistic(0,60,found, foundable,shared,kFALSE);
545 if (found<20) continue;
546 if (pt->GetNShared()/float(pt->GetNumberOfClusters())>0.2) continue;
549 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
550 iotrack.SetTPCPoints(pt->GetPoints());
551 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
552 iotrack.SetV0Indexes(pt->GetV0Indexes());
553 //iotrack.SetTPCpid(pt->fTPCr);
554 //iotrack.SetTPCindex(i);
555 fEvent->AddTrack(&iotrack);
558 // short tracks - secondaties
560 if ( (pt->GetNumberOfClusters()>30) ) {
561 Int_t found,foundable,shared;
562 pt->GetClusterStatistic(128,158,found, foundable,shared,kFALSE);
563 if ( (found>20) && (pt->GetNShared()/float(pt->GetNumberOfClusters())<0.2) &&float(found)/float(foundable)>0.8){
565 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
566 iotrack.SetTPCPoints(pt->GetPoints());
567 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
568 iotrack.SetV0Indexes(pt->GetV0Indexes());
569 //iotrack.SetTPCpid(pt->fTPCr);
570 //iotrack.SetTPCindex(i);
571 fEvent->AddTrack(&iotrack);
576 if ( (pt->GetNumberOfClusters()>15)) {
577 Int_t found,foundable,shared;
578 pt->GetClusterStatistic(138,158,found, foundable,shared,kFALSE);
579 if (found<15) continue;
580 if (foundable<=0) continue;
581 if (pt->GetNShared()/float(pt->GetNumberOfClusters())>0.2) continue;
582 if (float(found)/float(foundable)<0.8) continue;
585 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
586 iotrack.SetTPCPoints(pt->GetPoints());
587 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
588 iotrack.SetV0Indexes(pt->GetV0Indexes());
589 // iotrack.SetTPCpid(pt->fTPCr);
590 //iotrack.SetTPCindex(i);
591 fEvent->AddTrack(&iotrack);
595 // >> account for suppressed tracks in the kink indices (RS)
596 int nESDtracks = fEvent->GetNumberOfTracks();
597 for (int it=nESDtracks;it--;) {
598 AliESDtrack* esdTr = fEvent->GetTrack(it);
599 if (!esdTr || !esdTr->GetKinkIndex(0)) continue;
600 for (int ik=0;ik<3;ik++) {
602 if (!(knkId=esdTr->GetKinkIndex(ik))) break; // no more kinks for this track
603 AliESDkink* kink = fEvent->GetKink(TMath::Abs(knkId)-1);
605 AliError(Form("ESDTrack%d refers to non-existing kink %d",it,TMath::Abs(knkId)-1));
608 kink->SetIndex(it, knkId<0 ? 0:1); // update track index of the kink: mother at 0, daughter at 1
611 // << account for suppressed tracks in the kink indices (RS)
613 AliInfo(Form("Number of filled ESDs-\t%d\n",fEvent->GetNumberOfTracks()));
620 Double_t AliTPCtrackerMI::ErrY2(AliTPCseed* seed, const AliTPCclusterMI * cl){
623 // Use calibrated cluster error from OCDB
625 AliTPCClusterParam * clparam = AliTPCcalibDB::Instance()->GetClusterParam();
627 Float_t z = TMath::Abs(fkParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
628 Int_t ctype = cl->GetType();
629 Int_t type = (cl->GetRow()<63) ? 0: (cl->GetRow()>126) ? 1:2;
630 Double_t angle = seed->GetSnp()*seed->GetSnp();
631 angle = TMath::Sqrt(TMath::Abs(angle/(1.-angle)));
632 Double_t erry2 = clparam->GetError0Par(0,type, z,angle);
634 erry2+=0.5; // edge cluster
637 seed->SetErrorY2(erry2);
641 //calculate look-up table at the beginning
642 // static Bool_t ginit = kFALSE;
643 // static Float_t gnoise1,gnoise2,gnoise3;
644 // static Float_t ggg1[10000];
645 // static Float_t ggg2[10000];
646 // static Float_t ggg3[10000];
647 // static Float_t glandau1[10000];
648 // static Float_t glandau2[10000];
649 // static Float_t glandau3[10000];
651 // static Float_t gcor01[500];
652 // static Float_t gcor02[500];
653 // static Float_t gcorp[500];
657 // if (ginit==kFALSE){
658 // for (Int_t i=1;i<500;i++){
659 // Float_t rsigma = float(i)/100.;
660 // gcor02[i] = TMath::Max(0.78 +TMath::Exp(7.4*(rsigma-1.2)),0.6);
661 // gcor01[i] = TMath::Max(0.72 +TMath::Exp(3.36*(rsigma-1.2)),0.6);
662 // gcorp[i] = TMath::Max(TMath::Power((rsigma+0.5),1.5),1.2);
666 // for (Int_t i=3;i<10000;i++){
670 // Float_t amp = float(i);
671 // Float_t padlength =0.75;
672 // gnoise1 = 0.0004/padlength;
673 // Float_t nel = 0.268*amp;
674 // Float_t nprim = 0.155*amp;
675 // ggg1[i] = fkParam->GetDiffT()*fkParam->GetDiffT()*(2+0.001*nel/(padlength*padlength))/nel;
676 // glandau1[i] = (2.+0.12*nprim)*0.5* (2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
677 // if (glandau1[i]>1) glandau1[i]=1;
678 // glandau1[i]*=padlength*padlength/12.;
682 // gnoise2 = 0.0004/padlength;
684 // nprim = 0.133*amp;
685 // ggg2[i] = fkParam->GetDiffT()*fkParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
686 // glandau2[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
687 // if (glandau2[i]>1) glandau2[i]=1;
688 // glandau2[i]*=padlength*padlength/12.;
693 // gnoise3 = 0.0004/padlength;
695 // nprim = 0.133*amp;
696 // ggg3[i] = fkParam->GetDiffT()*fkParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
697 // glandau3[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
698 // if (glandau3[i]>1) glandau3[i]=1;
699 // glandau3[i]*=padlength*padlength/12.;
707 // Int_t amp = int(TMath::Abs(cl->GetQ()));
709 // seed->SetErrorY2(1.);
713 // Float_t z = TMath::Abs(fkParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
714 // Int_t ctype = cl->GetType();
715 // Float_t padlength= GetPadPitchLength(seed->GetRow());
716 // Double_t angle2 = seed->GetSnp()*seed->GetSnp();
717 // angle2 = angle2/(1-angle2);
719 // //cluster "quality"
720 // Int_t rsigmay = int(100.*cl->GetSigmaY2()/(seed->GetCurrentSigmaY2()));
723 // if (fSectors==fInnerSec){
724 // snoise2 = gnoise1;
725 // res = ggg1[amp]*z+glandau1[amp]*angle2;
726 // if (ctype==0) res *= gcor01[rsigmay];
729 // res*= gcorp[rsigmay];
733 // if (padlength<1.1){
734 // snoise2 = gnoise2;
735 // res = ggg2[amp]*z+glandau2[amp]*angle2;
736 // if (ctype==0) res *= gcor02[rsigmay];
739 // res*= gcorp[rsigmay];
743 // snoise2 = gnoise3;
744 // res = ggg3[amp]*z+glandau3[amp]*angle2;
745 // if (ctype==0) res *= gcor02[rsigmay];
748 // res*= gcorp[rsigmay];
755 // res*=2.4; // overestimate error 2 times
759 // if (res<2*snoise2)
762 // seed->SetErrorY2(res);
770 Double_t AliTPCtrackerMI::ErrZ2(AliTPCseed* seed, const AliTPCclusterMI * cl){
773 // Use calibrated cluster error from OCDB
775 AliTPCClusterParam * clparam = AliTPCcalibDB::Instance()->GetClusterParam();
777 Float_t z = TMath::Abs(fkParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
778 Int_t ctype = cl->GetType();
779 Int_t type = (cl->GetRow()<63) ? 0: (cl->GetRow()>126) ? 1:2;
781 Double_t angle2 = seed->GetSnp()*seed->GetSnp();
782 angle2 = seed->GetTgl()*seed->GetTgl()*(1+angle2/(1-angle2));
783 Double_t angle = TMath::Sqrt(TMath::Abs(angle2));
784 Double_t errz2 = clparam->GetError0Par(1,type, z,angle);
786 errz2+=0.5; // edge cluster
789 seed->SetErrorZ2(errz2);
795 // //seed->SetErrorY2(0.1);
797 // //calculate look-up table at the beginning
798 // static Bool_t ginit = kFALSE;
799 // static Float_t gnoise1,gnoise2,gnoise3;
800 // static Float_t ggg1[10000];
801 // static Float_t ggg2[10000];
802 // static Float_t ggg3[10000];
803 // static Float_t glandau1[10000];
804 // static Float_t glandau2[10000];
805 // static Float_t glandau3[10000];
807 // static Float_t gcor01[1000];
808 // static Float_t gcor02[1000];
809 // static Float_t gcorp[1000];
813 // if (ginit==kFALSE){
814 // for (Int_t i=1;i<1000;i++){
815 // Float_t rsigma = float(i)/100.;
816 // gcor02[i] = TMath::Max(0.81 +TMath::Exp(6.8*(rsigma-1.2)),0.6);
817 // gcor01[i] = TMath::Max(0.72 +TMath::Exp(2.04*(rsigma-1.2)),0.6);
818 // gcorp[i] = TMath::Max(TMath::Power((rsigma+0.5),1.5),1.2);
822 // for (Int_t i=3;i<10000;i++){
826 // Float_t amp = float(i);
827 // Float_t padlength =0.75;
828 // gnoise1 = 0.0004/padlength;
829 // Float_t nel = 0.268*amp;
830 // Float_t nprim = 0.155*amp;
831 // ggg1[i] = fkParam->GetDiffT()*fkParam->GetDiffT()*(2+0.001*nel/(padlength*padlength))/nel;
832 // glandau1[i] = (2.+0.12*nprim)*0.5* (2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
833 // if (glandau1[i]>1) glandau1[i]=1;
834 // glandau1[i]*=padlength*padlength/12.;
838 // gnoise2 = 0.0004/padlength;
840 // nprim = 0.133*amp;
841 // ggg2[i] = fkParam->GetDiffT()*fkParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
842 // glandau2[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
843 // if (glandau2[i]>1) glandau2[i]=1;
844 // glandau2[i]*=padlength*padlength/12.;
849 // gnoise3 = 0.0004/padlength;
851 // nprim = 0.133*amp;
852 // ggg3[i] = fkParam->GetDiffT()*fkParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
853 // glandau3[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
854 // if (glandau3[i]>1) glandau3[i]=1;
855 // glandau3[i]*=padlength*padlength/12.;
863 // Int_t amp = int(TMath::Abs(cl->GetQ()));
865 // seed->SetErrorY2(1.);
869 // Float_t z = TMath::Abs(fkParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
870 // Int_t ctype = cl->GetType();
871 // Float_t padlength= GetPadPitchLength(seed->GetRow());
873 // Double_t angle2 = seed->GetSnp()*seed->GetSnp();
874 // // if (angle2<0.6) angle2 = 0.6;
875 // angle2 = seed->GetTgl()*seed->GetTgl()*(1+angle2/(1-angle2));
877 // //cluster "quality"
878 // Int_t rsigmaz = int(100.*cl->GetSigmaZ2()/(seed->GetCurrentSigmaZ2()));
881 // if (fSectors==fInnerSec){
882 // snoise2 = gnoise1;
883 // res = ggg1[amp]*z+glandau1[amp]*angle2;
884 // if (ctype==0) res *= gcor01[rsigmaz];
887 // res*= gcorp[rsigmaz];
891 // if (padlength<1.1){
892 // snoise2 = gnoise2;
893 // res = ggg2[amp]*z+glandau2[amp]*angle2;
894 // if (ctype==0) res *= gcor02[rsigmaz];
897 // res*= gcorp[rsigmaz];
901 // snoise2 = gnoise3;
902 // res = ggg3[amp]*z+glandau3[amp]*angle2;
903 // if (ctype==0) res *= gcor02[rsigmaz];
906 // res*= gcorp[rsigmaz];
915 // if ((ctype<0) &&<70){
920 // if (res<2*snoise2)
922 // if (res>3) res =3;
923 // seed->SetErrorZ2(res);
931 void AliTPCtrackerMI::RotateToLocal(AliTPCseed *seed)
933 //rotate to track "local coordinata
934 Float_t x = seed->GetX();
935 Float_t y = seed->GetY();
936 Float_t ymax = x*TMath::Tan(0.5*fSectors->GetAlpha());
939 seed->SetRelativeSector((seed->GetRelativeSector()+1) % fN);
940 if (!seed->Rotate(fSectors->GetAlpha()))
942 } else if (y <-ymax) {
943 seed->SetRelativeSector((seed->GetRelativeSector()-1+fN) % fN);
944 if (!seed->Rotate(-fSectors->GetAlpha()))
952 //_____________________________________________________________________________
953 Double_t AliTPCtrackerMI::F1old(Double_t x1,Double_t y1,
954 Double_t x2,Double_t y2,
955 Double_t x3,Double_t y3) const
957 //-----------------------------------------------------------------
958 // Initial approximation of the track curvature
959 //-----------------------------------------------------------------
960 Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
961 Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
962 (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
963 Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
964 (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
966 Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
967 if ( xr*xr+yr*yr<=0.00000000000001) return 100;
968 return -xr*yr/sqrt(xr*xr+yr*yr);
973 //_____________________________________________________________________________
974 Double_t AliTPCtrackerMI::F1(Double_t x1,Double_t y1,
975 Double_t x2,Double_t y2,
976 Double_t x3,Double_t y3) const
978 //-----------------------------------------------------------------
979 // Initial approximation of the track curvature
980 //-----------------------------------------------------------------
986 Double_t det = x3*y2-x2*y3;
987 if (TMath::Abs(det)<1e-10){
991 Double_t u = 0.5* (x2*(x2-x3)+y2*(y2-y3))/det;
992 Double_t x0 = x3*0.5-y3*u;
993 Double_t y0 = y3*0.5+x3*u;
994 Double_t c2 = 1/TMath::Sqrt(x0*x0+y0*y0);
1000 Double_t AliTPCtrackerMI::F2(Double_t x1,Double_t y1,
1001 Double_t x2,Double_t y2,
1002 Double_t x3,Double_t y3) const
1004 //-----------------------------------------------------------------
1005 // Initial approximation of the track curvature
1006 //-----------------------------------------------------------------
1012 Double_t det = x3*y2-x2*y3;
1013 if (TMath::Abs(det)<1e-10) {
1017 Double_t u = 0.5* (x2*(x2-x3)+y2*(y2-y3))/det;
1018 Double_t x0 = x3*0.5-y3*u;
1019 Double_t y0 = y3*0.5+x3*u;
1020 Double_t c2 = 1/TMath::Sqrt(x0*x0+y0*y0);
1029 //_____________________________________________________________________________
1030 Double_t AliTPCtrackerMI::F2old(Double_t x1,Double_t y1,
1031 Double_t x2,Double_t y2,
1032 Double_t x3,Double_t y3) const
1034 //-----------------------------------------------------------------
1035 // Initial approximation of the track curvature times center of curvature
1036 //-----------------------------------------------------------------
1037 Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
1038 Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
1039 (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
1040 Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
1041 (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
1043 Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
1045 return -a/(d*y1-b)*xr/sqrt(xr*xr+yr*yr);
1048 //_____________________________________________________________________________
1049 Double_t AliTPCtrackerMI::F3(Double_t x1,Double_t y1,
1050 Double_t x2,Double_t y2,
1051 Double_t z1,Double_t z2) const
1053 //-----------------------------------------------------------------
1054 // Initial approximation of the tangent of the track dip angle
1055 //-----------------------------------------------------------------
1056 return (z1 - z2)/sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
1060 Double_t AliTPCtrackerMI::F3n(Double_t x1,Double_t y1,
1061 Double_t x2,Double_t y2,
1062 Double_t z1,Double_t z2, Double_t c) const
1064 //-----------------------------------------------------------------
1065 // Initial approximation of the tangent of the track dip angle
1066 //-----------------------------------------------------------------
1070 //angle1 = (z1-z2)*c/(TMath::ASin(c*x1-ni)-TMath::ASin(c*x2-ni));
1072 Double_t d = TMath::Sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
1073 if (TMath::Abs(d*c*0.5)>1) return 0;
1074 // Double_t angle2 = TMath::ASin(d*c*0.5);
1075 // Double_t angle2 = AliTPCFastMath::FastAsin(d*c*0.5);
1076 Double_t angle2 = (d*c*0.5>0.1)? TMath::ASin(d*c*0.5): AliTPCFastMath::FastAsin(d*c*0.5);
1078 angle2 = (z1-z2)*c/(angle2*2.);
1082 Bool_t AliTPCtrackerMI::GetProlongation(Double_t x1, Double_t x2, Double_t x[5], Double_t &y, Double_t &z) const
1083 {//-----------------------------------------------------------------
1084 // This function find proloncation of a track to a reference plane x=x2.
1085 //-----------------------------------------------------------------
1089 if (TMath::Abs(x[4]*x1 - x[2]) >= 0.999) {
1093 Double_t c1=x[4]*x1 - x[2], r1=TMath::Sqrt((1.-c1)*(1.+c1));
1094 Double_t c2=x[4]*x2 - x[2], r2=TMath::Sqrt((1.-c2)*(1.+c2));
1098 Double_t dy = dx*(c1+c2)/(r1+r2);
1101 Double_t delta = x[4]*dx*(c1+c2)/(c1*r2 + c2*r1);
1103 if (TMath::Abs(delta)>0.01){
1104 dz = x[3]*TMath::ASin(delta)/x[4];
1106 dz = x[3]*AliTPCFastMath::FastAsin(delta)/x[4];
1109 //dz = x[3]*AliTPCFastMath::FastAsin(delta)/x[4];
1117 Int_t AliTPCtrackerMI::LoadClusters (TTree *const tree)
1122 return LoadClusters();
1126 Int_t AliTPCtrackerMI::LoadClusters(const TObjArray *arr)
1129 // load clusters to the memory
1130 AliTPCClustersRow *clrow = new AliTPCClustersRow("AliTPCclusterMI");
1131 Int_t lower = arr->LowerBound();
1132 Int_t entries = arr->GetEntriesFast();
1134 for (Int_t i=lower; i<entries; i++) {
1135 clrow = (AliTPCClustersRow*) arr->At(i);
1136 if(!clrow) continue;
1137 if(!clrow->GetArray()) continue;
1141 fkParam->AdjustSectorRow(clrow->GetID(),sec,row);
1143 for (Int_t icl=0; icl<clrow->GetArray()->GetEntriesFast(); icl++){
1144 Transform((AliTPCclusterMI*)(clrow->GetArray()->At(icl)));
1147 if (clrow->GetArray()->GetEntriesFast()<=0) continue;
1148 AliTPCtrackerRow * tpcrow=0;
1151 tpcrow = &(fInnerSec[sec%fkNIS][row]);
1155 tpcrow = &(fOuterSec[(sec-fkNIS*2)%fkNOS][row]);
1156 left = (sec-fkNIS*2)/fkNOS;
1159 tpcrow->SetN1(clrow->GetArray()->GetEntriesFast());
1160 tpcrow->SetClusters1(new AliTPCclusterMI[tpcrow->GetN1()]);
1161 for (Int_t j=0;j<tpcrow->GetN1();++j)
1162 tpcrow->SetCluster1(j, *(AliTPCclusterMI*)(clrow->GetArray()->At(j)));
1165 tpcrow->SetN2(clrow->GetArray()->GetEntriesFast());
1166 tpcrow->SetClusters2(new AliTPCclusterMI[tpcrow->GetN2()]);
1167 for (Int_t j=0;j<tpcrow->GetN2();++j)
1168 tpcrow->SetCluster2(j,*(AliTPCclusterMI*)(clrow->GetArray()->At(j)));
1170 clrow->GetArray()->Clear("C");
1179 Int_t AliTPCtrackerMI::LoadClusters(const TClonesArray *arr)
1182 // load clusters to the memory from one
1185 AliTPCclusterMI *clust=0;
1186 Int_t count[72][96] = { {0} , {0} };
1188 // loop over clusters
1189 for (Int_t icl=0; icl<arr->GetEntriesFast(); icl++) {
1190 clust = (AliTPCclusterMI*)arr->At(icl);
1191 if(!clust) continue;
1192 //printf("cluster: det %d, row %d \n", clust->GetDetector(),clust->GetRow());
1194 // transform clusters
1197 // count clusters per pad row
1198 count[clust->GetDetector()][clust->GetRow()]++;
1201 // insert clusters to sectors
1202 for (Int_t icl=0; icl<arr->GetEntriesFast(); icl++) {
1203 clust = (AliTPCclusterMI*)arr->At(icl);
1204 if(!clust) continue;
1206 Int_t sec = clust->GetDetector();
1207 Int_t row = clust->GetRow();
1209 // filter overlapping pad rows needed by HLT
1210 if(sec<fkNIS*2) { //IROCs
1211 if(row == 30) continue;
1214 if(row == 27 || row == 76) continue;
1220 fInnerSec[sec%fkNIS].InsertCluster(clust, count[sec][row], fkParam);
1223 left = (sec-fkNIS*2)/fkNOS;
1224 fOuterSec[(sec-fkNIS*2)%fkNOS].InsertCluster(clust, count[sec][row], fkParam);
1228 // Load functions must be called behind LoadCluster(TClonesArray*)
1230 //LoadOuterSectors();
1231 //LoadInnerSectors();
1237 Int_t AliTPCtrackerMI::LoadClusters()
1240 // load clusters to the memory
1241 AliTPCClustersRow *clrow= new AliTPCClustersRow("AliTPCclusterMI");
1243 // TTree * tree = fClustersArray.GetTree();
1245 TTree * tree = fInput;
1246 TBranch * br = tree->GetBranch("Segment");
1247 br->SetAddress(&clrow);
1249 Int_t j=Int_t(tree->GetEntries());
1250 for (Int_t i=0; i<j; i++) {
1254 fkParam->AdjustSectorRow(clrow->GetID(),sec,row);
1255 for (Int_t icl=0; icl<clrow->GetArray()->GetEntriesFast(); icl++){
1256 Transform((AliTPCclusterMI*)(clrow->GetArray()->At(icl)));
1259 AliTPCtrackerRow * tpcrow=0;
1262 tpcrow = &(fInnerSec[sec%fkNIS][row]);
1266 tpcrow = &(fOuterSec[(sec-fkNIS*2)%fkNOS][row]);
1267 left = (sec-fkNIS*2)/fkNOS;
1270 tpcrow->SetN1(clrow->GetArray()->GetEntriesFast());
1271 tpcrow->SetClusters1(new AliTPCclusterMI[tpcrow->GetN1()]);
1272 for (Int_t k=0;k<tpcrow->GetN1();++k)
1273 tpcrow->SetCluster1(k, *(AliTPCclusterMI*)(clrow->GetArray()->At(k)));
1276 tpcrow->SetN2(clrow->GetArray()->GetEntriesFast());
1277 tpcrow->SetClusters2(new AliTPCclusterMI[tpcrow->GetN2()]);
1278 for (Int_t k=0;k<tpcrow->GetN2();++k)
1279 tpcrow->SetCluster2(k,*(AliTPCclusterMI*)(clrow->GetArray()->At(k)));
1290 void AliTPCtrackerMI::UnloadClusters()
1293 // unload clusters from the memory
1295 Int_t nrows = fOuterSec->GetNRows();
1296 for (Int_t sec = 0;sec<fkNOS;sec++)
1297 for (Int_t row = 0;row<nrows;row++){
1298 AliTPCtrackerRow* tpcrow = &(fOuterSec[sec%fkNOS][row]);
1300 // if (tpcrow->fClusters1) delete []tpcrow->fClusters1;
1301 // if (tpcrow->fClusters2) delete []tpcrow->fClusters2;
1303 tpcrow->ResetClusters();
1306 nrows = fInnerSec->GetNRows();
1307 for (Int_t sec = 0;sec<fkNIS;sec++)
1308 for (Int_t row = 0;row<nrows;row++){
1309 AliTPCtrackerRow* tpcrow = &(fInnerSec[sec%fkNIS][row]);
1311 // if (tpcrow->fClusters1) delete []tpcrow->fClusters1;
1312 //if (tpcrow->fClusters2) delete []tpcrow->fClusters2;
1314 tpcrow->ResetClusters();
1320 void AliTPCtrackerMI::FillClusterArray(TObjArray* array) const{
1322 // Filling cluster to the array - For visualization purposes
1325 nrows = fOuterSec->GetNRows();
1326 for (Int_t sec = 0;sec<fkNOS;sec++)
1327 for (Int_t row = 0;row<nrows;row++){
1328 AliTPCtrackerRow* tpcrow = &(fOuterSec[sec%fkNOS][row]);
1329 if (!tpcrow) continue;
1330 for (Int_t icl = 0;icl<tpcrow->GetN();icl++){
1331 array->AddLast((TObject*)((*tpcrow)[icl]));
1334 nrows = fInnerSec->GetNRows();
1335 for (Int_t sec = 0;sec<fkNIS;sec++)
1336 for (Int_t row = 0;row<nrows;row++){
1337 AliTPCtrackerRow* tpcrow = &(fInnerSec[sec%fkNIS][row]);
1338 if (!tpcrow) continue;
1339 for (Int_t icl = 0;icl<tpcrow->GetN();icl++){
1340 array->AddLast((TObject*)(*tpcrow)[icl]);
1346 void AliTPCtrackerMI::Transform(AliTPCclusterMI * cluster){
1350 AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ;
1352 AliFatal("Tranformations not in calibDB");
1354 transform->SetCurrentRecoParam((AliTPCRecoParam*)AliTPCReconstructor::GetRecoParam());
1355 Double_t x[3]={cluster->GetRow(),cluster->GetPad(),cluster->GetTimeBin()};
1356 Int_t i[1]={cluster->GetDetector()};
1357 transform->Transform(x,i,0,1);
1358 // if (cluster->GetDetector()%36>17){
1363 // in debug mode check the transformation
1365 if (AliTPCReconstructor::StreamLevel()>2) {
1367 cluster->GetGlobalXYZ(gx);
1368 Int_t event = (fEvent==NULL)? 0: fEvent->GetEventNumberInFile();
1369 TTreeSRedirector &cstream = *fDebugStreamer;
1370 cstream<<"Transform"<<
1381 cluster->SetX(x[0]);
1382 cluster->SetY(x[1]);
1383 cluster->SetZ(x[2]);
1389 //if (!fkParam->IsGeoRead()) fkParam->ReadGeoMatrices();
1390 TGeoHMatrix *mat = fkParam->GetClusterMatrix(cluster->GetDetector());
1392 Double_t pos[3]= {cluster->GetX(),cluster->GetY(),cluster->GetZ()};
1393 Double_t posC[3]={cluster->GetX(),cluster->GetY(),cluster->GetZ()};
1394 if (mat) mat->LocalToMaster(pos,posC);
1396 // chack Loading of Geo matrices from GeoManager - TEMPORARY FIX
1398 cluster->SetX(posC[0]);
1399 cluster->SetY(posC[1]);
1400 cluster->SetZ(posC[2]);
1403 //_____________________________________________________________________________
1404 Int_t AliTPCtrackerMI::LoadOuterSectors() {
1405 //-----------------------------------------------------------------
1406 // This function fills outer TPC sectors with clusters.
1407 //-----------------------------------------------------------------
1408 Int_t nrows = fOuterSec->GetNRows();
1410 for (Int_t sec = 0;sec<fkNOS;sec++)
1411 for (Int_t row = 0;row<nrows;row++){
1412 AliTPCtrackerRow* tpcrow = &(fOuterSec[sec%fkNOS][row]);
1413 Int_t sec2 = sec+2*fkNIS;
1415 Int_t ncl = tpcrow->GetN1();
1417 AliTPCclusterMI *c= (tpcrow->GetCluster1(ncl));
1418 index=(((sec2<<8)+row)<<16)+ncl;
1419 tpcrow->InsertCluster(c,index);
1422 ncl = tpcrow->GetN2();
1424 AliTPCclusterMI *c= (tpcrow->GetCluster2(ncl));
1425 index=((((sec2+fkNOS)<<8)+row)<<16)+ncl;
1426 tpcrow->InsertCluster(c,index);
1429 // write indexes for fast acces
1431 for (Int_t i=0;i<510;i++)
1432 tpcrow->SetFastCluster(i,-1);
1433 for (Int_t i=0;i<tpcrow->GetN();i++){
1434 Int_t zi = Int_t((*tpcrow)[i]->GetZ()+255.);
1435 tpcrow->SetFastCluster(zi,i); // write index
1438 for (Int_t i=0;i<510;i++){
1439 if (tpcrow->GetFastCluster(i)<0)
1440 tpcrow->SetFastCluster(i,last);
1442 last = tpcrow->GetFastCluster(i);
1451 //_____________________________________________________________________________
1452 Int_t AliTPCtrackerMI::LoadInnerSectors() {
1453 //-----------------------------------------------------------------
1454 // This function fills inner TPC sectors with clusters.
1455 //-----------------------------------------------------------------
1456 Int_t nrows = fInnerSec->GetNRows();
1458 for (Int_t sec = 0;sec<fkNIS;sec++)
1459 for (Int_t row = 0;row<nrows;row++){
1460 AliTPCtrackerRow* tpcrow = &(fInnerSec[sec%fkNIS][row]);
1463 Int_t ncl = tpcrow->GetN1();
1465 AliTPCclusterMI *c= (tpcrow->GetCluster1(ncl));
1466 index=(((sec<<8)+row)<<16)+ncl;
1467 tpcrow->InsertCluster(c,index);
1470 ncl = tpcrow->GetN2();
1472 AliTPCclusterMI *c= (tpcrow->GetCluster2(ncl));
1473 index=((((sec+fkNIS)<<8)+row)<<16)+ncl;
1474 tpcrow->InsertCluster(c,index);
1477 // write indexes for fast acces
1479 for (Int_t i=0;i<510;i++)
1480 tpcrow->SetFastCluster(i,-1);
1481 for (Int_t i=0;i<tpcrow->GetN();i++){
1482 Int_t zi = Int_t((*tpcrow)[i]->GetZ()+255.);
1483 tpcrow->SetFastCluster(zi,i); // write index
1486 for (Int_t i=0;i<510;i++){
1487 if (tpcrow->GetFastCluster(i)<0)
1488 tpcrow->SetFastCluster(i,last);
1490 last = tpcrow->GetFastCluster(i);
1502 //_________________________________________________________________________
1503 AliTPCclusterMI *AliTPCtrackerMI::GetClusterMI(Int_t index) const {
1504 //--------------------------------------------------------------------
1505 // Return pointer to a given cluster
1506 //--------------------------------------------------------------------
1507 if (index<0) return 0; // no cluster
1508 Int_t sec=(index&0xff000000)>>24;
1509 Int_t row=(index&0x00ff0000)>>16;
1510 Int_t ncl=(index&0x00007fff)>>00;
1512 const AliTPCtrackerRow * tpcrow=0;
1513 AliTPCclusterMI * clrow =0;
1515 if (sec<0 || sec>=fkNIS*4) {
1516 AliWarning(Form("Wrong sector %d",sec));
1521 AliTPCtrackerSector& tracksec = fInnerSec[sec%fkNIS];
1522 if (tracksec.GetNRows()<=row) return 0;
1523 tpcrow = &(tracksec[row]);
1524 if (tpcrow==0) return 0;
1527 if (tpcrow->GetN1()<=ncl) return 0;
1528 clrow = tpcrow->GetClusters1();
1531 if (tpcrow->GetN2()<=ncl) return 0;
1532 clrow = tpcrow->GetClusters2();
1536 AliTPCtrackerSector& tracksec = fOuterSec[(sec-fkNIS*2)%fkNOS];
1537 if (tracksec.GetNRows()<=row) return 0;
1538 tpcrow = &(tracksec[row]);
1539 if (tpcrow==0) return 0;
1541 if (sec-2*fkNIS<fkNOS) {
1542 if (tpcrow->GetN1()<=ncl) return 0;
1543 clrow = tpcrow->GetClusters1();
1546 if (tpcrow->GetN2()<=ncl) return 0;
1547 clrow = tpcrow->GetClusters2();
1551 return &(clrow[ncl]);
1557 Int_t AliTPCtrackerMI::FollowToNext(AliTPCseed& t, Int_t nr) {
1558 //-----------------------------------------------------------------
1559 // This function tries to find a track prolongation to next pad row
1560 //-----------------------------------------------------------------
1562 Double_t x= GetXrow(nr), ymax=GetMaxY(nr);
1565 AliTPCclusterMI *cl=0;
1566 Int_t tpcindex= t.GetClusterIndex2(nr);
1568 // update current shape info every 5 pad-row
1569 // if ( (nr%5==0) || t.GetNumberOfClusters()<2 || (t.fCurrentSigmaY2<0.0001) ){
1573 if (fIteration>0 && tpcindex>=-1){ //if we have already clusters
1575 if (tpcindex==-1) return 0; //track in dead zone
1577 cl = t.GetClusterPointer(nr);
1578 if ( (cl==0) ) cl = GetClusterMI(tpcindex);
1579 t.SetCurrentClusterIndex1(tpcindex);
1582 Int_t relativesector = ((tpcindex&0xff000000)>>24)%18; // if previously accepted cluster in different sector
1583 Float_t angle = relativesector*fSectors->GetAlpha()+fSectors->GetAlphaShift();
1585 if (angle<-TMath::Pi()) angle += 2*TMath::Pi();
1586 if (angle>=TMath::Pi()) angle -= 2*TMath::Pi();
1588 if (TMath::Abs(angle-t.GetAlpha())>0.001){
1589 Double_t rotation = angle-t.GetAlpha();
1590 t.SetRelativeSector(relativesector);
1591 if (!t.Rotate(rotation)) return 0;
1593 if (!t.PropagateTo(x)) return 0;
1595 t.SetCurrentCluster(cl);
1597 Int_t accept = AcceptCluster(&t,t.GetCurrentCluster());
1598 if ((tpcindex&0x8000)==0) accept =0;
1600 //if founded cluster is acceptible
1601 if (cl->IsUsed(11)) { // id cluster is shared inrease uncertainty
1602 t.SetErrorY2(t.GetErrorY2()+0.03);
1603 t.SetErrorZ2(t.GetErrorZ2()+0.03);
1604 t.SetErrorY2(t.GetErrorY2()*3);
1605 t.SetErrorZ2(t.GetErrorZ2()*3);
1607 t.SetNFoundable(t.GetNFoundable()+1);
1608 UpdateTrack(&t,accept);
1613 if (TMath::Abs(t.GetSnp())>AliTPCReconstructor::GetMaxSnpTracker()) return 0; // cut on angle
1614 if (fIteration>1 && IsFindable(t)){
1615 // not look for new cluster during refitting
1616 t.SetNFoundable(t.GetNFoundable()+1);
1621 // if (TMath::Abs(t.GetSnp())>0.95 || TMath::Abs(x*t.GetC()-t.GetEta())>0.95) return 0;// patch 28 fev 06
1622 Double_t y=t.GetYat(x);
1623 if (TMath::Abs(y)>ymax){
1625 t.SetRelativeSector((t.GetRelativeSector()+1) % fN);
1626 if (!t.Rotate(fSectors->GetAlpha()))
1628 } else if (y <-ymax) {
1629 t.SetRelativeSector((t.GetRelativeSector()-1+fN) % fN);
1630 if (!t.Rotate(-fSectors->GetAlpha()))
1636 if (!t.PropagateTo(x)) {
1637 if (fIteration==0) t.SetRemoval(10);
1641 Double_t z=t.GetZ();
1644 if (!IsActive(t.GetRelativeSector(),nr)) {
1646 t.SetClusterIndex2(nr,-1);
1649 //AliInfo(Form("A - Sector%d phi %f - alpha %f", t.fRelativeSector,y/x, t.GetAlpha()));
1650 Bool_t isActive = IsActive(t.GetRelativeSector(),nr);
1651 Bool_t isActive2 = (nr>=fInnerSec->GetNRows()) ? fOuterSec[t.GetRelativeSector()][nr-fInnerSec->GetNRows()].GetN()>0:fInnerSec[t.GetRelativeSector()][nr].GetN()>0;
1653 if (!isActive || !isActive2) return 0;
1655 const AliTPCtrackerRow &krow=GetRow(t.GetRelativeSector(),nr);
1656 if ( (t.GetSigmaY2()<0) || t.GetSigmaZ2()<0) return 0;
1658 Double_t roadz = 1.;
1660 if (TMath::Abs(TMath::Abs(y)-ymax)<krow.GetDeadZone()){
1662 t.SetClusterIndex2(nr,-1);
1668 // if (TMath::Abs(z)<(AliTPCReconstructor::GetCtgRange()*x+10) && TMath::Abs(z)<fkParam->GetZLength(0) && (TMath::Abs(t.GetSnp())<AliTPCReconstructor::GetMaxSnpTracker()))
1669 t.SetNFoundable(t.GetNFoundable()+1);
1675 // cl = krow.FindNearest2(y+10.,z,roady,roadz,index);
1676 cl = krow.FindNearest2(y,z,roady,roadz,index);
1677 if (cl) t.SetCurrentClusterIndex1(krow.GetIndex(index));
1680 t.SetCurrentCluster(cl);
1682 if (fIteration==2&&cl->IsUsed(10)) return 0;
1683 Int_t accept = AcceptCluster(&t,t.GetCurrentCluster());
1684 if (fIteration==2&&cl->IsUsed(11)) {
1685 t.SetErrorY2(t.GetErrorY2()+0.03);
1686 t.SetErrorZ2(t.GetErrorZ2()+0.03);
1687 t.SetErrorY2(t.GetErrorY2()*3);
1688 t.SetErrorZ2(t.GetErrorZ2()*3);
1691 if (t.fCurrentCluster->IsUsed(10)){
1696 if (t.fNShared>0.7*t.GetNumberOfClusters()) {
1702 if (accept<3) UpdateTrack(&t,accept);
1705 if ( fIteration==0 && t.GetNFoundable()*0.5 > t.GetNumberOfClusters()) t.SetRemoval(10);
1713 //_________________________________________________________________________
1714 Bool_t AliTPCtrackerMI::GetTrackPoint(Int_t index, AliTrackPoint &p ) const
1716 // Get track space point by index
1717 // return false in case the cluster doesn't exist
1718 AliTPCclusterMI *cl = GetClusterMI(index);
1719 if (!cl) return kFALSE;
1720 Int_t sector = (index&0xff000000)>>24;
1721 // Int_t row = (index&0x00ff0000)>>16;
1723 // xyz[0] = fkParam->GetPadRowRadii(sector,row);
1724 xyz[0] = cl->GetX();
1725 xyz[1] = cl->GetY();
1726 xyz[2] = cl->GetZ();
1728 fkParam->AdjustCosSin(sector,cos,sin);
1729 Float_t x = cos*xyz[0]-sin*xyz[1];
1730 Float_t y = cos*xyz[1]+sin*xyz[0];
1732 Float_t sigmaY2 = 0.027*cl->GetSigmaY2();
1733 if (sector < fkParam->GetNInnerSector()) sigmaY2 *= 2.07;
1734 Float_t sigmaZ2 = 0.066*cl->GetSigmaZ2();
1735 if (sector < fkParam->GetNInnerSector()) sigmaZ2 *= 1.77;
1736 cov[0] = sin*sin*sigmaY2;
1737 cov[1] = -sin*cos*sigmaY2;
1739 cov[3] = cos*cos*sigmaY2;
1742 p.SetXYZ(x,y,xyz[2],cov);
1743 AliGeomManager::ELayerID iLayer;
1745 if (sector < fkParam->GetNInnerSector()) {
1746 iLayer = AliGeomManager::kTPC1;
1750 iLayer = AliGeomManager::kTPC2;
1751 idet = sector - fkParam->GetNInnerSector();
1753 UShort_t volid = AliGeomManager::LayerToVolUID(iLayer,idet);
1754 p.SetVolumeID(volid);
1760 Int_t AliTPCtrackerMI::UpdateClusters(AliTPCseed& t, Int_t nr) {
1761 //-----------------------------------------------------------------
1762 // This function tries to find a track prolongation to next pad row
1763 //-----------------------------------------------------------------
1764 t.SetCurrentCluster(0);
1765 t.SetCurrentClusterIndex1(0);
1767 Double_t xt=t.GetX();
1768 Int_t row = GetRowNumber(xt)-1;
1769 Double_t ymax= GetMaxY(nr);
1771 if (row < nr) return 1; // don't prolongate if not information until now -
1772 // if (TMath::Abs(t.GetSnp())>0.9 && t.GetNumberOfClusters()>40. && fIteration!=2) {
1774 // return 0; // not prolongate strongly inclined tracks
1776 // if (TMath::Abs(t.GetSnp())>0.95) {
1778 // return 0; // not prolongate strongly inclined tracks
1779 // }// patch 28 fev 06
1781 Double_t x= GetXrow(nr);
1783 //t.PropagateTo(x+0.02);
1784 //t.PropagateTo(x+0.01);
1785 if (!t.PropagateTo(x)){
1792 if (TMath::Abs(y)>ymax){
1794 t.SetRelativeSector((t.GetRelativeSector()+1) % fN);
1795 if (!t.Rotate(fSectors->GetAlpha()))
1797 } else if (y <-ymax) {
1798 t.SetRelativeSector((t.GetRelativeSector()-1+fN) % fN);
1799 if (!t.Rotate(-fSectors->GetAlpha()))
1802 // if (!t.PropagateTo(x)){
1809 if (TMath::Abs(t.GetSnp())>AliTPCReconstructor::GetMaxSnpTracker()) return 0;
1811 if (!IsActive(t.GetRelativeSector(),nr)) {
1813 t.SetClusterIndex2(nr,-1);
1816 //AliInfo(Form("A - Sector%d phi %f - alpha %f", t.fRelativeSector,y/x, t.GetAlpha()));
1818 AliTPCtrackerRow &krow=GetRow(t.GetRelativeSector(),nr);
1820 if (TMath::Abs(TMath::Abs(y)-ymax)<krow.GetDeadZone()){
1822 t.SetClusterIndex2(nr,-1);
1828 // if (TMath::Abs(t.GetZ())<(AliTPCReconstructor::GetCtgRange()*t.GetX()+10) && (TMath::Abs(t.GetSnp())<AliTPCReconstructor::GetMaxSnpTracker()))
1829 if (IsFindable(t)) t.SetNFoundable(t.GetNFoundable()+1);
1835 if ( (nr%2==0) || t.GetNumberOfClusters()<2){
1836 // t.fCurrentSigmaY = GetSigmaY(&t);
1837 //t.fCurrentSigmaZ = GetSigmaZ(&t);
1841 AliTPCclusterMI *cl=0;
1844 Double_t roady = 1.;
1845 Double_t roadz = 1.;
1849 index = t.GetClusterIndex2(nr);
1850 if ( (index>0) && (index&0x8000)==0){
1851 cl = t.GetClusterPointer(nr);
1852 if ( (cl==0) && (index>0)) cl = GetClusterMI(index);
1853 t.SetCurrentClusterIndex1(index);
1855 t.SetCurrentCluster(cl);
1861 // if (index<0) return 0;
1862 UInt_t uindex = TMath::Abs(index);
1865 //cl = krow.FindNearest2(y+10,z,roady,roadz,uindex);
1866 cl = krow.FindNearest2(y,z,roady,roadz,uindex);
1869 if (cl) t.SetCurrentClusterIndex1(krow.GetIndex(uindex));
1870 t.SetCurrentCluster(cl);
1876 Int_t AliTPCtrackerMI::FollowToNextCluster(AliTPCseed & t, Int_t nr) {
1877 //-----------------------------------------------------------------
1878 // This function tries to find a track prolongation to next pad row
1879 //-----------------------------------------------------------------
1881 //update error according neighborhoud
1883 if (t.GetCurrentCluster()) {
1885 Int_t accept = AcceptCluster(&t,t.GetCurrentCluster());
1887 if (t.GetCurrentCluster()->IsUsed(10)){
1892 t.SetNShared(t.GetNShared()+1);
1893 if (t.GetNShared()>0.7*t.GetNumberOfClusters()) {
1898 if (fIteration>0) accept = 0;
1899 if (accept<3) UpdateTrack(&t,accept);
1903 if ( ( (t.GetSigmaY2()+t.GetSigmaZ2())>0.16)&& t.GetNumberOfClusters()>18) t.SetRemoval(10);
1904 if ( t.GetChi2()/t.GetNumberOfClusters()>6 &&t.GetNumberOfClusters()>18) t.SetRemoval(10);
1906 if (( (t.GetNFoundable()*0.5 > t.GetNumberOfClusters()) || t.GetNoCluster()>15)) t.SetRemoval(10);
1914 //_____________________________________________________________________________
1915 Int_t AliTPCtrackerMI::FollowProlongation(AliTPCseed& t, Int_t rf, Int_t step) {
1916 //-----------------------------------------------------------------
1917 // This function tries to find a track prolongation.
1918 //-----------------------------------------------------------------
1919 Double_t xt=t.GetX();
1921 Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
1922 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
1923 if (alpha < 0. ) alpha += 2.*TMath::Pi();
1925 t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN);
1927 Int_t first = GetRowNumber(xt)-1;
1928 for (Int_t nr= first; nr>=rf; nr-=step) {
1930 if (t.GetKinkIndexes()[0]>0){
1931 for (Int_t i=0;i<3;i++){
1932 Int_t index = t.GetKinkIndexes()[i];
1933 if (index==0) break;
1934 if (index<0) continue;
1936 AliKink * kink = (AliKink*)fEvent->GetKink(index-1);
1938 printf("PROBLEM\n");
1941 Int_t kinkrow = kink->GetTPCRow0()+2+Int_t(0.5/(0.05+kink->GetAngle(2)));
1943 AliExternalTrackParam paramd(t);
1944 kink->SetDaughter(paramd);
1945 kink->SetStatus(2,5);
1952 if (nr==80) t.UpdateReference();
1953 if (nr<fInnerSec->GetNRows())
1954 fSectors = fInnerSec;
1956 fSectors = fOuterSec;
1957 if (FollowToNext(t,nr)==0)
1970 Int_t AliTPCtrackerMI::FollowBackProlongation(AliTPCseed& t, Int_t rf) {
1971 //-----------------------------------------------------------------
1972 // This function tries to find a track prolongation.
1973 //-----------------------------------------------------------------
1975 Double_t xt=t.GetX();
1976 Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
1977 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
1978 if (alpha < 0. ) alpha += 2.*TMath::Pi();
1979 t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN);
1981 Int_t first = t.GetFirstPoint();
1982 if (first<GetRowNumber(xt)+1) first = GetRowNumber(xt)+1;
1984 if (first<0) first=0;
1985 for (Int_t nr=first; nr<=rf; nr++) {
1986 // if ( (TMath::Abs(t.GetSnp())>0.95)) break;//patch 28 fev 06
1987 if (t.GetKinkIndexes()[0]<0){
1988 for (Int_t i=0;i<3;i++){
1989 Int_t index = t.GetKinkIndexes()[i];
1990 if (index==0) break;
1991 if (index>0) continue;
1992 index = TMath::Abs(index);
1993 AliKink * kink = (AliKink*)fEvent->GetKink(index-1);
1995 printf("PROBLEM\n");
1998 Int_t kinkrow = kink->GetTPCRow0()-2-Int_t(0.5/(0.05+kink->GetAngle(2)));
2000 AliExternalTrackParam paramm(t);
2001 kink->SetMother(paramm);
2002 kink->SetStatus(2,1);
2009 if (nr<fInnerSec->GetNRows())
2010 fSectors = fInnerSec;
2012 fSectors = fOuterSec;
2023 Float_t AliTPCtrackerMI::OverlapFactor(AliTPCseed * s1, AliTPCseed * s2, Int_t &sum1, Int_t & sum2)
2031 Float_t dz2 =(s1->GetZ() - s2->GetZ());
2034 Float_t dy2 =TMath::Abs((s1->GetY() - s2->GetY()));
2036 Float_t distance = TMath::Sqrt(dz2+dy2);
2037 if (distance>4.) return 0; // if there are far away - not overlap - to reduce combinatorics
2040 Int_t firstpoint = TMath::Min(s1->GetFirstPoint(),s2->GetFirstPoint());
2041 Int_t lastpoint = TMath::Max(s1->GetLastPoint(),s2->GetLastPoint());
2046 if (firstpoint>lastpoint) {
2047 firstpoint =lastpoint;
2052 for (Int_t i=firstpoint-1;i<lastpoint+1;i++){
2053 if (s1->GetClusterIndex2(i)>0) sum1++;
2054 if (s2->GetClusterIndex2(i)>0) sum2++;
2055 if (s1->GetClusterIndex2(i)==s2->GetClusterIndex2(i) && s1->GetClusterIndex2(i)>0) {
2059 if (sum<5) return 0;
2061 Float_t summin = TMath::Min(sum1+1,sum2+1);
2062 Float_t ratio = (sum+1)/Float_t(summin);
2066 void AliTPCtrackerMI::SignShared(AliTPCseed * s1, AliTPCseed * s2)
2070 Float_t thetaCut = 0.2;//+10.*TMath::Sqrt(s1->GetSigmaTglZ()+ s2->GetSigmaTglZ());
2071 if (TMath::Abs(s1->GetTgl()-s2->GetTgl())>thetaCut) return;
2072 Float_t minCl = TMath::Min(s1->GetNumberOfClusters(),s2->GetNumberOfClusters());
2073 Int_t cutN0 = TMath::Max(5,TMath::Nint(0.1*minCl));
2078 //Int_t firstpoint = TMath::Max(s1->GetFirstPoint(),s2->GetFirstPoint());
2079 //Int_t lastpoint = TMath::Min(s1->GetLastPoint(),s2->GetLastPoint());
2080 Int_t firstpoint = 0;
2081 Int_t lastpoint = 160;
2083 // if (firstpoint>=lastpoint-5) return;;
2085 for (Int_t i=firstpoint;i<lastpoint;i++){
2086 // if ( (s1->GetClusterIndex2(i)&0xFFFF8FFF)==(s2->GetClusterIndex2(i)&0xFFFF8FFF) && s1->GetClusterIndex2(i)>0) {
2087 if ( (s1->GetClusterIndex2(i))==(s2->GetClusterIndex2(i)) && s1->GetClusterIndex2(i)>0) {
2089 s1->SetSharedMapBit(i, kTRUE);
2090 s2->SetSharedMapBit(i, kTRUE);
2092 if (s1->GetClusterIndex2(i)>0)
2093 s1->SetClusterMapBit(i, kTRUE);
2095 if (sumshared>cutN0){
2098 for (Int_t i=firstpoint;i<lastpoint;i++){
2099 // if ( (s1->GetClusterIndex2(i)&0xFFFF8FFF)==(s2->GetClusterIndex2(i)&0xFFFF8FFF) && s1->GetClusterIndex2(i)>0) {
2100 if ( (s1->GetClusterIndex2(i))==(s2->GetClusterIndex2(i)) && s1->GetClusterIndex2(i)>0) {
2101 AliTPCTrackerPoint *p1 = s1->GetTrackPoint(i);
2102 AliTPCTrackerPoint *p2 = s2->GetTrackPoint(i);;
2103 if (s1->IsActive()&&s2->IsActive()){
2104 p1->SetShared(kTRUE);
2105 p2->SetShared(kTRUE);
2111 if (sumshared>cutN0){
2112 for (Int_t i=0;i<4;i++){
2113 if (s1->GetOverlapLabel(3*i)==0){
2114 s1->SetOverlapLabel(3*i, s2->GetLabel());
2115 s1->SetOverlapLabel(3*i+1,sumshared);
2116 s1->SetOverlapLabel(3*i+2,s2->GetUniqueID());
2120 for (Int_t i=0;i<4;i++){
2121 if (s2->GetOverlapLabel(3*i)==0){
2122 s2->SetOverlapLabel(3*i, s1->GetLabel());
2123 s2->SetOverlapLabel(3*i+1,sumshared);
2124 s2->SetOverlapLabel(3*i+2,s1->GetUniqueID());
2131 void AliTPCtrackerMI::SignShared(TObjArray * arr)
2134 //sort trackss according sectors
2136 for (Int_t i=0; i<arr->GetEntriesFast(); i++) {
2137 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2139 //if (pt) RotateToLocal(pt);
2143 arr->Sort(); // sorting according relative sectors
2144 arr->Expand(arr->GetEntries());
2147 Int_t nseed=arr->GetEntriesFast();
2148 for (Int_t i=0; i<nseed; i++) {
2149 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2151 for (Int_t j=0;j<=12;j++){
2152 pt->SetOverlapLabel(j,0);
2155 for (Int_t i=0; i<nseed; i++) {
2156 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2158 if (pt->GetRemoval()>10) continue;
2159 for (Int_t j=i+1; j<nseed; j++){
2160 AliTPCseed *pt2=(AliTPCseed*)arr->UncheckedAt(j);
2161 if (TMath::Abs(pt->GetRelativeSector()-pt2->GetRelativeSector())>1) continue;
2163 if (pt2->GetRemoval()<=10) {
2164 //if ( TMath::Abs(pt->GetRelativeSector()-pt2->GetRelativeSector())>0) break;
2172 void AliTPCtrackerMI::SortTracks(TObjArray * arr, Int_t mode) const
2175 //sort tracks in array according mode criteria
2176 Int_t nseed = arr->GetEntriesFast();
2177 for (Int_t i=0; i<nseed; i++) {
2178 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2189 void AliTPCtrackerMI::RemoveUsed2(TObjArray * arr, Float_t factor1, Float_t factor2, Int_t minimal)
2192 // Loop over all tracks and remove overlaped tracks (with lower quality)
2194 // 1. Unsign clusters
2195 // 2. Sort tracks according quality
2196 // Quality is defined by the number of cluster between first and last points
2198 // 3. Loop over tracks - decreasing quality order
2199 // a.) remove - If the fraction of shared cluster less than factor (1- n or 2)
2200 // b.) remove - If the minimal number of clusters less than minimal and not ITS
2201 // c.) if track accepted - sign clusters
2203 //Called in - AliTPCtrackerMI::Clusters2Tracks()
2204 // - AliTPCtrackerMI::PropagateBack()
2205 // - AliTPCtrackerMI::RefitInward()
2208 // factor1 - factor for constrained
2209 // factor2 - for non constrained tracks
2210 // if (Float_t(shared+1)/Float_t(found+1)>factor) - DELETE
2214 Int_t nseed = arr->GetEntriesFast();
2215 Float_t * quality = new Float_t[nseed];
2216 Int_t * indexes = new Int_t[nseed];
2220 for (Int_t i=0; i<nseed; i++) {
2221 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2226 pt->UpdatePoints(); //select first last max dens points
2227 Float_t * points = pt->GetPoints();
2228 if (points[3]<0.8) quality[i] =-1;
2229 quality[i] = (points[2]-points[0])+pt->GetNumberOfClusters();
2230 //prefer high momenta tracks if overlaps
2231 quality[i] *= TMath::Sqrt(TMath::Abs(pt->Pt())+0.5);
2233 TMath::Sort(nseed,quality,indexes);
2236 for (Int_t itrack=0; itrack<nseed; itrack++) {
2237 Int_t trackindex = indexes[itrack];
2238 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(trackindex);
2241 if (quality[trackindex]<0){
2243 delete arr->RemoveAt(trackindex);
2246 arr->RemoveAt(trackindex);
2252 Int_t first = Int_t(pt->GetPoints()[0]);
2253 Int_t last = Int_t(pt->GetPoints()[2]);
2254 Double_t factor = (pt->GetBConstrain()) ? factor1: factor2;
2256 Int_t found,foundable,shared;
2257 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
2258 // pt->GetClusterStatistic(0,160, found, foundable,shared,kFALSE);
2259 Bool_t itsgold =kFALSE;
2262 if (pt->GetESD()->GetITSclusters(dummy)>4) itsgold= kTRUE;
2266 if (Float_t(shared+1)/Float_t(found+1)>factor){
2267 if (pt->GetKinkIndexes()[0]!=0) continue; //don't remove tracks - part of the kinks
2268 if( AliTPCReconstructor::StreamLevel()>3){
2269 TTreeSRedirector &cstream = *fDebugStreamer;
2270 cstream<<"RemoveUsed"<<
2271 "iter="<<fIteration<<
2275 delete arr->RemoveAt(trackindex);
2278 if (pt->GetNumberOfClusters()<50&&(found-0.5*shared)<minimal){ //remove short tracks
2279 if (pt->GetKinkIndexes()[0]!=0) continue; //don't remove tracks - part of the kinks
2280 if( AliTPCReconstructor::StreamLevel()>3){
2281 TTreeSRedirector &cstream = *fDebugStreamer;
2282 cstream<<"RemoveShort"<<
2283 "iter="<<fIteration<<
2287 delete arr->RemoveAt(trackindex);
2293 //if (sharedfactor>0.4) continue;
2294 if (pt->GetKinkIndexes()[0]>0) continue;
2295 //Remove tracks with undefined properties - seems
2296 if (pt->GetSigmaY2()<kAlmost0) continue; // ? what is the origin ?
2298 for (Int_t i=first; i<last; i++) {
2299 Int_t index=pt->GetClusterIndex2(i);
2300 // if (index<0 || index&0x8000 ) continue;
2301 if (index<0 || index&0x8000 ) continue;
2302 AliTPCclusterMI *c= pt->GetClusterPointer(i);
2309 Info("RemoveUsed","\n*****\nNumber of good tracks after shared removal\t%d\n",fNtracks);
2315 void AliTPCtrackerMI::DumpClusters(Int_t iter, TObjArray *trackArray)
2318 // Dump clusters after reco
2319 // signed and unsigned cluster can be visualized
2320 // 1. Unsign all cluster
2321 // 2. Sign all used clusters
2324 Int_t nseed = trackArray->GetEntries();
2325 for (Int_t i=0; i<nseed; i++){
2326 AliTPCseed *pt=(AliTPCseed*)trackArray->UncheckedAt(i);
2330 Bool_t isKink=pt->GetKinkIndex(0)!=0;
2331 for (Int_t j=0; j<160; ++j) {
2332 Int_t index=pt->GetClusterIndex2(j);
2333 if (index<0) continue;
2334 AliTPCclusterMI *c= pt->GetClusterPointer(j);
2336 if (isKink) c->Use(100); // kink
2337 c->Use(10); // by default usage 10
2342 for (Int_t sec=0;sec<fkNIS;sec++){
2343 for (Int_t row=0;row<fInnerSec->GetNRows();row++){
2344 AliTPCclusterMI *cl = fInnerSec[sec][row].GetClusters1();
2345 for (Int_t icl =0;icl< fInnerSec[sec][row].GetN1();icl++){
2346 Float_t gx[3]; cl[icl].GetGlobalXYZ(gx);
2347 (*fDebugStreamer)<<"clDump"<<
2355 cl = fInnerSec[sec][row].GetClusters2();
2356 for (Int_t icl =0;icl< fInnerSec[sec][row].GetN2();icl++){
2357 Float_t gx[3]; cl[icl].GetGlobalXYZ(gx);
2358 (*fDebugStreamer)<<"clDump"<<
2369 for (Int_t sec=0;sec<fkNOS;sec++){
2370 for (Int_t row=0;row<fOuterSec->GetNRows();row++){
2371 AliTPCclusterMI *cl = fOuterSec[sec][row].GetClusters1();
2372 for (Int_t icl =0;icl< fOuterSec[sec][row].GetN1();icl++){
2373 Float_t gx[3]; cl[icl].GetGlobalXYZ(gx);
2374 (*fDebugStreamer)<<"clDump"<<
2382 cl = fOuterSec[sec][row].GetClusters2();
2383 for (Int_t icl =0;icl< fOuterSec[sec][row].GetN2();icl++){
2384 Float_t gx[3]; cl[icl].GetGlobalXYZ(gx);
2385 (*fDebugStreamer)<<"clDump"<<
2397 void AliTPCtrackerMI::UnsignClusters()
2400 // loop over all clusters and unsign them
2403 for (Int_t sec=0;sec<fkNIS;sec++){
2404 for (Int_t row=0;row<fInnerSec->GetNRows();row++){
2405 AliTPCclusterMI *cl = fInnerSec[sec][row].GetClusters1();
2406 for (Int_t icl =0;icl< fInnerSec[sec][row].GetN1();icl++)
2407 // if (cl[icl].IsUsed(10))
2409 cl = fInnerSec[sec][row].GetClusters2();
2410 for (Int_t icl =0;icl< fInnerSec[sec][row].GetN2();icl++)
2411 //if (cl[icl].IsUsed(10))
2416 for (Int_t sec=0;sec<fkNOS;sec++){
2417 for (Int_t row=0;row<fOuterSec->GetNRows();row++){
2418 AliTPCclusterMI *cl = fOuterSec[sec][row].GetClusters1();
2419 for (Int_t icl =0;icl< fOuterSec[sec][row].GetN1();icl++)
2420 //if (cl[icl].IsUsed(10))
2422 cl = fOuterSec[sec][row].GetClusters2();
2423 for (Int_t icl =0;icl< fOuterSec[sec][row].GetN2();icl++)
2424 //if (cl[icl].IsUsed(10))
2433 void AliTPCtrackerMI::SignClusters(const TObjArray * arr, Float_t fnumber, Float_t fdensity)
2436 //sign clusters to be "used"
2438 // snumber and sdensity sign number of sigmas - bellow mean value to be accepted
2439 // loop over "primaries"
2453 Int_t nseed = arr->GetEntriesFast();
2454 for (Int_t i=0; i<nseed; i++) {
2455 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2459 if (!(pt->IsActive())) continue;
2460 Float_t dens = pt->GetNumberOfClusters()/Float_t(pt->GetNFoundable());
2461 if ( (dens>0.7) && (pt->GetNumberOfClusters()>70)){
2463 sumdens2+= dens*dens;
2464 sumn += pt->GetNumberOfClusters();
2465 sumn2 += pt->GetNumberOfClusters()*pt->GetNumberOfClusters();
2466 Float_t chi2 = pt->GetChi2()/pt->GetNumberOfClusters();
2469 sumchi2 +=chi2*chi2;
2474 Float_t mdensity = 0.9;
2475 Float_t meann = 130;
2476 Float_t meanchi = 1;
2477 Float_t sdensity = 0.1;
2478 Float_t smeann = 10;
2479 Float_t smeanchi =0.4;
2483 mdensity = sumdens/sum;
2485 meanchi = sumchi/sum;
2487 sdensity = sumdens2/sum-mdensity*mdensity;
2489 sdensity = TMath::Sqrt(sdensity);
2493 smeann = sumn2/sum-meann*meann;
2495 smeann = TMath::Sqrt(smeann);
2499 smeanchi = sumchi2/sum - meanchi*meanchi;
2501 smeanchi = TMath::Sqrt(smeanchi);
2507 //REMOVE SHORT DELTAS or tracks going out of sensitive volume of TPC
2509 for (Int_t i=0; i<nseed; i++) {
2510 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2514 if (pt->GetBSigned()) continue;
2515 if (pt->GetBConstrain()) continue;
2516 //if (!(pt->IsActive())) continue;
2518 Int_t found,foundable,shared;
2519 pt->GetClusterStatistic(0,160,found, foundable,shared);
2520 if (shared/float(found)>0.3) {
2521 if (shared/float(found)>0.9 ){
2522 //delete arr->RemoveAt(i);
2527 Bool_t isok =kFALSE;
2528 if ( (pt->GetNShared()/pt->GetNumberOfClusters()<0.5) &&pt->GetNumberOfClusters()>60)
2530 if ((TMath::Abs(1/pt->GetC())<100.) && (pt->GetNShared()/pt->GetNumberOfClusters()<0.7))
2532 if (TMath::Abs(pt->GetZ()/pt->GetX())>1.1)
2534 if ( (TMath::Abs(pt->GetSnp()>0.7) && pt->GetD(0,0)>60.))
2538 for (Int_t j=0; j<160; ++j) {
2539 Int_t index=pt->GetClusterIndex2(j);
2540 if (index<0) continue;
2541 AliTPCclusterMI *c= pt->GetClusterPointer(j);
2543 //if (!(c->IsUsed(10))) c->Use();
2550 Double_t maxchi = meanchi+2.*smeanchi;
2552 for (Int_t i=0; i<nseed; i++) {
2553 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2557 //if (!(pt->IsActive())) continue;
2558 if (pt->GetBSigned()) continue;
2559 Double_t chi = pt->GetChi2()/pt->GetNumberOfClusters();
2560 if (chi>maxchi) continue;
2563 Float_t dens = pt->GetNumberOfClusters()/Float_t(pt->GetNFoundable());
2565 //sign only tracks with enoug big density at the beginning
2567 if ((pt->GetDensityFirst(40)<0.75) && pt->GetNumberOfClusters()<meann) continue;
2570 Double_t mindens = TMath::Max(double(mdensity-sdensity*fdensity*bfactor),0.65);
2571 Double_t minn = TMath::Max(Int_t(meann-fnumber*smeann*bfactor),50);
2573 // if (pt->fBConstrain) mindens = TMath::Max(mdensity-sdensity*fdensity*bfactor,0.65);
2574 if ( (pt->GetRemoval()==10) && (pt->GetSnp()>0.8)&&(dens>mindens))
2577 if ((dens>mindens && pt->GetNumberOfClusters()>minn) && chi<maxchi ){
2578 //Int_t noc=pt->GetNumberOfClusters();
2579 pt->SetBSigned(kTRUE);
2580 for (Int_t j=0; j<160; ++j) {
2582 Int_t index=pt->GetClusterIndex2(j);
2583 if (index<0) continue;
2584 AliTPCclusterMI *c= pt->GetClusterPointer(j);
2586 // if (!(c->IsUsed(10))) c->Use();
2591 // gLastCheck = nseed;
2599 void AliTPCtrackerMI::StopNotActive(const TObjArray * arr, Int_t row0, Float_t th0, Float_t th1, Float_t th2) const
2601 // stop not active tracks
2602 // take th1 as threshold for number of founded to number of foundable on last 10 active rows
2603 // take th2 as threshold for number of founded to number of foundable on last 20 active rows
2604 Int_t nseed = arr->GetEntriesFast();
2606 for (Int_t i=0; i<nseed; i++) {
2607 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2611 if (!(pt->IsActive())) continue;
2612 StopNotActive(pt,row0,th0, th1,th2);
2618 void AliTPCtrackerMI::StopNotActive(AliTPCseed * seed, Int_t row0, Float_t th0, Float_t th1,
2621 // stop not active tracks
2622 // take th1 as threshold for number of founded to number of foundable on last 10 active rows
2623 // take th2 as threshold for number of founded to number of foundable on last 20 active rows
2626 Int_t foundable = 0;
2627 Int_t maxindex = seed->GetLastPoint(); //last foundable row
2628 if (seed->GetNFoundable()*th0 > seed->GetNumberOfClusters()) {
2629 seed->Desactivate(10) ;
2633 for (Int_t i=row0; i<maxindex; i++){
2634 Int_t index = seed->GetClusterIndex2(i);
2635 if (index!=-1) foundable++;
2637 if (foundable<=30) sumgood1++;
2638 if (foundable<=50) {
2645 if (foundable>=30.){
2646 if (sumgood1<(th1*30.)) seed->Desactivate(10);
2649 if (sumgood2<(th2*50.)) seed->Desactivate(10);
2653 Int_t AliTPCtrackerMI::RefitInward(AliESDEvent *event)
2656 // back propagation of ESD tracks
2659 const Int_t kMaxFriendTracks=2000;
2663 //PrepareForProlongation(fSeeds,1);
2664 PropagateForward2(fSeeds);
2665 RemoveUsed2(fSeeds,0.4,0.4,20);
2667 TObjArray arraySeed(fSeeds->GetEntries());
2668 for (Int_t i=0;i<fSeeds->GetEntries();i++) {
2669 arraySeed.AddAt(fSeeds->At(i),i);
2671 SignShared(&arraySeed);
2672 // FindCurling(fSeeds, event,2); // find multi found tracks
2673 FindSplitted(fSeeds, event,2); // find multi found tracks
2674 if (AliTPCReconstructor::StreamLevel()>5) FindMultiMC(fSeeds, fEvent,2); // find multi found tracks
2677 Int_t nseed = fSeeds->GetEntriesFast();
2678 for (Int_t i=0;i<nseed;i++){
2679 AliTPCseed * seed = (AliTPCseed*) fSeeds->UncheckedAt(i);
2680 if (!seed) continue;
2681 if (seed->GetKinkIndex(0)>0) UpdateKinkQualityD(seed); // update quality informations for kinks
2682 AliESDtrack *esd=event->GetTrack(i);
2684 if (seed->GetNumberOfClusters()<60 && seed->GetNumberOfClusters()<(esd->GetTPCclusters(0) -5)*0.8){
2685 AliExternalTrackParam paramIn;
2686 AliExternalTrackParam paramOut;
2687 Int_t ncl = seed->RefitTrack(seed,¶mIn,¶mOut);
2688 if (AliTPCReconstructor::StreamLevel()>2) {
2689 (*fDebugStreamer)<<"RecoverIn"<<
2693 "pout.="<<¶mOut<<
2698 seed->Set(paramIn.GetX(),paramIn.GetAlpha(),paramIn.GetParameter(),paramIn.GetCovariance());
2699 seed->SetNumberOfClusters(ncl);
2703 seed->PropagateTo(fkParam->GetInnerRadiusLow());
2704 seed->UpdatePoints();
2705 AddCovariance(seed);
2707 seed->CookdEdx(0.02,0.6);
2708 CookLabel(seed,0.1); //For comparison only
2709 esd->SetTPCClusterMap(seed->GetClusterMap());
2710 esd->SetTPCSharedMap(seed->GetSharedMap());
2712 if (AliTPCReconstructor::StreamLevel()>1 && seed!=0&&esd!=0) {
2713 TTreeSRedirector &cstream = *fDebugStreamer;
2720 if (seed->GetNumberOfClusters()>15){
2721 esd->UpdateTrackParams(seed,AliESDtrack::kTPCrefit);
2722 esd->SetTPCPoints(seed->GetPoints());
2723 esd->SetTPCPointsF(seed->GetNFoundable());
2724 Int_t ndedx = seed->GetNCDEDX(0);
2725 Float_t sdedx = seed->GetSDEDX(0);
2726 Float_t dedx = seed->GetdEdx();
2727 esd->SetTPCsignal(dedx, sdedx, ndedx);
2729 // add seed to the esd track in Calib level
2731 Bool_t storeFriend = gRandom->Rndm()<(kMaxFriendTracks)/Float_t(nseed);
2732 if (AliTPCReconstructor::StreamLevel()>0 &&storeFriend){
2733 AliTPCseed * seedCopy = new AliTPCseed(*seed, kTRUE);
2734 esd->AddCalibObject(seedCopy);
2739 //printf("problem\n");
2742 //FindKinks(fSeeds,event);
2743 if (AliTPCReconstructor::StreamLevel()>3) DumpClusters(2,fSeeds);
2744 Info("RefitInward","Number of refitted tracks %d",ntracks);
2749 Int_t AliTPCtrackerMI::PropagateBack(AliESDEvent *event)
2752 // back propagation of ESD tracks
2758 PropagateBack(fSeeds);
2759 RemoveUsed2(fSeeds,0.4,0.4,20);
2760 //FindCurling(fSeeds, fEvent,1);
2761 FindSplitted(fSeeds, event,1); // find multi found tracks
2762 if (AliTPCReconstructor::StreamLevel()>5) FindMultiMC(fSeeds, fEvent,1); // find multi found tracks
2765 Int_t nseed = fSeeds->GetEntriesFast();
2767 for (Int_t i=0;i<nseed;i++){
2768 AliTPCseed * seed = (AliTPCseed*) fSeeds->UncheckedAt(i);
2769 if (!seed) continue;
2770 if (seed->GetKinkIndex(0)<0) UpdateKinkQualityM(seed); // update quality informations for kinks
2771 seed->UpdatePoints();
2772 AddCovariance(seed);
2773 AliESDtrack *esd=event->GetTrack(i);
2774 if (seed->GetNumberOfClusters()<60 && seed->GetNumberOfClusters()<(esd->GetTPCclusters(0) -5)*0.8){
2775 AliExternalTrackParam paramIn;
2776 AliExternalTrackParam paramOut;
2777 Int_t ncl = seed->RefitTrack(seed,¶mIn,¶mOut);
2778 if (AliTPCReconstructor::StreamLevel()>2) {
2779 (*fDebugStreamer)<<"RecoverBack"<<
2783 "pout.="<<¶mOut<<
2788 seed->Set(paramOut.GetX(),paramOut.GetAlpha(),paramOut.GetParameter(),paramOut.GetCovariance());
2789 seed->SetNumberOfClusters(ncl);
2792 seed->CookdEdx(0.02,0.6);
2793 CookLabel(seed,0.1); //For comparison only
2794 if (seed->GetNumberOfClusters()>15){
2795 esd->UpdateTrackParams(seed,AliESDtrack::kTPCout);
2796 esd->SetTPCPoints(seed->GetPoints());
2797 esd->SetTPCPointsF(seed->GetNFoundable());
2798 Int_t ndedx = seed->GetNCDEDX(0);
2799 Float_t sdedx = seed->GetSDEDX(0);
2800 Float_t dedx = seed->GetdEdx();
2801 esd->SetTPCsignal(dedx, sdedx, ndedx);
2803 Int_t eventnumber = event->GetEventNumberInFile();// patch 28 fev 06
2804 // This is most likely NOT the event number you'd like to use. It has nothing to do with the 'real' event number
2805 if (AliTPCReconstructor::StreamLevel()>1 && esd) {
2806 (*fDebugStreamer)<<"Cback"<<
2809 "EventNrInFile="<<eventnumber<<
2814 if (AliTPCReconstructor::StreamLevel()>3) DumpClusters(1,fSeeds);
2815 //FindKinks(fSeeds,event);
2816 Info("PropagateBack","Number of back propagated tracks %d",ntracks);
2823 void AliTPCtrackerMI::DeleteSeeds()
2828 Int_t nseed = fSeeds->GetEntriesFast();
2829 for (Int_t i=0;i<nseed;i++){
2830 AliTPCseed * seed = (AliTPCseed*)fSeeds->At(i);
2831 if (seed) delete fSeeds->RemoveAt(i);
2838 void AliTPCtrackerMI::ReadSeeds(AliESDEvent *const event, Int_t direction)
2841 //read seeds from the event
2843 Int_t nentr=event->GetNumberOfTracks();
2845 Info("ReadSeeds", "Number of ESD tracks: %d\n", nentr);
2850 fSeeds = new TObjArray(nentr);
2854 for (Int_t i=0; i<nentr; i++) {
2855 AliESDtrack *esd=event->GetTrack(i);
2856 ULong_t status=esd->GetStatus();
2857 if (!(status&AliESDtrack::kTPCin)) continue;
2858 AliTPCtrack t(*esd);
2859 t.SetNumberOfClusters(0);
2860 // AliTPCseed *seed = new AliTPCseed(t,t.GetAlpha());
2861 AliTPCseed *seed = new AliTPCseed(t/*,t.GetAlpha()*/);
2862 seed->SetUniqueID(esd->GetID());
2863 AddCovariance(seed); //add systematic ucertainty
2864 for (Int_t ikink=0;ikink<3;ikink++) {
2865 Int_t index = esd->GetKinkIndex(ikink);
2866 seed->GetKinkIndexes()[ikink] = index;
2867 if (index==0) continue;
2868 index = TMath::Abs(index);
2869 AliESDkink * kink = fEvent->GetKink(index-1);
2870 if (kink&&esd->GetKinkIndex(ikink)<0){
2871 if ((status & AliESDtrack::kTRDrefit) != 0) kink->SetStatus(1,2);
2872 if ((status & AliESDtrack::kITSout) != 0) kink->SetStatus(1,0);
2874 if (kink&&esd->GetKinkIndex(ikink)>0){
2875 if ((status & AliESDtrack::kTRDrefit) != 0) kink->SetStatus(1,6);
2876 if ((status & AliESDtrack::kITSout) != 0) kink->SetStatus(1,4);
2880 if (((status&AliESDtrack::kITSout)==0)&&(direction==1)) seed->ResetCovariance(10.);
2881 if ( direction ==2 &&(status & AliESDtrack::kTRDrefit) == 0 ) seed->ResetCovariance(10.);
2882 //if ( direction ==2 && ((status & AliESDtrack::kTPCout) == 0) ) {
2883 // fSeeds->AddAt(0,i);
2887 if ( direction ==2 &&(status & AliESDtrack::kTRDrefit) > 0 ) {
2888 Double_t par0[5],par1[5],alpha,x;
2889 esd->GetInnerExternalParameters(alpha,x,par0);
2890 esd->GetExternalParameters(x,par1);
2891 Double_t delta1 = TMath::Abs(par0[4]-par1[4])/(0.000000001+TMath::Abs(par0[4]+par1[4]));
2892 Double_t delta2 = TMath::Abs(par0[3]-par1[3]);
2894 if (esd->GetTRDncls()>0) trdchi2 = esd->GetTRDchi2()/esd->GetTRDncls();
2895 //reset covariance if suspicious
2896 if ( (delta1>0.1) || (delta2>0.006) ||trdchi2>7.)
2897 seed->ResetCovariance(10.);
2902 // rotate to the local coordinate system
2904 fSectors=fInnerSec; fN=fkNIS;
2905 Double_t alpha=seed->GetAlpha() - fSectors->GetAlphaShift();
2906 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
2907 if (alpha < 0. ) alpha += 2.*TMath::Pi();
2908 Int_t ns=Int_t(alpha/fSectors->GetAlpha())%fN;
2909 alpha =ns*fSectors->GetAlpha() + fSectors->GetAlphaShift();
2910 if (alpha<-TMath::Pi()) alpha += 2*TMath::Pi();
2911 if (alpha>=TMath::Pi()) alpha -= 2*TMath::Pi();
2912 alpha-=seed->GetAlpha();
2913 if (!seed->Rotate(alpha)) {
2919 if (esd->GetKinkIndex(0)<=0){
2920 for (Int_t irow=0;irow<160;irow++){
2921 Int_t index = seed->GetClusterIndex2(irow);
2924 AliTPCclusterMI * cl = GetClusterMI(index);
2925 seed->SetClusterPointer(irow,cl);
2927 if ((index & 0x8000)==0){
2928 cl->Use(10); // accepted cluster
2930 cl->Use(6); // close cluster not accepted
2933 Info("ReadSeeds","Not found cluster");
2938 fSeeds->AddAt(seed,i);
2944 //_____________________________________________________________________________
2945 void AliTPCtrackerMI::MakeSeeds3(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, Float_t cuts[4],
2946 Float_t deltay, Int_t ddsec) {
2947 //-----------------------------------------------------------------
2948 // This function creates track seeds.
2949 // SEEDING WITH VERTEX CONSTRAIN
2950 //-----------------------------------------------------------------
2951 // cuts[0] - fP4 cut
2952 // cuts[1] - tan(phi) cut
2953 // cuts[2] - zvertex cut
2954 // cuts[3] - fP3 cut
2962 Double_t x[5], c[15];
2963 // Int_t di = i1-i2;
2965 AliTPCseed * seed = new AliTPCseed();
2966 Double_t alpha=fSectors->GetAlpha(), shift=fSectors->GetAlphaShift();
2967 Double_t cs=cos(alpha), sn=sin(alpha);
2969 // Double_t x1 =fOuterSec->GetX(i1);
2970 //Double_t xx2=fOuterSec->GetX(i2);
2972 Double_t x1 =GetXrow(i1);
2973 Double_t xx2=GetXrow(i2);
2975 Double_t x3=GetX(), y3=GetY(), z3=GetZ();
2977 Int_t imiddle = (i2+i1)/2; //middle pad row index
2978 Double_t xm = GetXrow(imiddle); // radius of middle pad-row
2979 const AliTPCtrackerRow& krm=GetRow(sec,imiddle); //middle pad -row
2983 const AliTPCtrackerRow& kr1=GetRow(ns,i1);
2984 Double_t ymax = GetMaxY(i1)-kr1.GetDeadZone()-1.5;
2985 Double_t ymaxm = GetMaxY(imiddle)-kr1.GetDeadZone()-1.5;
2988 // change cut on curvature if it can't reach this layer
2989 // maximal curvature set to reach it
2990 Double_t dvertexmax = TMath::Sqrt((x1-x3)*(x1-x3)+(ymax+5-y3)*(ymax+5-y3));
2991 if (dvertexmax*0.5*cuts[0]>0.85){
2992 cuts[0] = 0.85/(dvertexmax*0.5+1.);
2994 Double_t r2min = 1/(cuts[0]*cuts[0]); //minimal square of radius given by cut
2997 if (deltay>0) ddsec = 0;
2998 // loop over clusters
2999 for (Int_t is=0; is < kr1; is++) {
3001 if (kr1[is]->IsUsed(10)) continue;
3002 Double_t y1=kr1[is]->GetY(), z1=kr1[is]->GetZ();
3003 //if (TMath::Abs(y1)>ymax) continue;
3005 if (deltay>0 && TMath::Abs(ymax-TMath::Abs(y1))> deltay ) continue; // seed only at the edge
3007 // find possible directions
3008 Float_t anglez = (z1-z3)/(x1-x3);
3009 Float_t extraz = z1 - anglez*(x1-xx2); // extrapolated z
3012 //find rotation angles relative to line given by vertex and point 1
3013 Double_t dvertex2 = (x1-x3)*(x1-x3)+(y1-y3)*(y1-y3);
3014 Double_t dvertex = TMath::Sqrt(dvertex2);
3015 Double_t angle13 = TMath::ATan((y1-y3)/(x1-x3));
3016 Double_t cs13 = cos(-angle13), sn13 = sin(-angle13);
3019 // loop over 2 sectors
3025 Double_t dddz1=0; // direction of delta inclination in z axis
3032 for (Int_t dsec = dsec1; dsec<=dsec2;dsec++){
3033 Int_t sec2 = sec + dsec;
3035 // AliTPCtrackerRow& kr2 = fOuterSec[(sec2+fkNOS)%fkNOS][i2];
3036 //AliTPCtrackerRow& kr2m = fOuterSec[(sec2+fkNOS)%fkNOS][imiddle];
3037 AliTPCtrackerRow& kr2 = GetRow((sec2+fkNOS)%fkNOS,i2);
3038 AliTPCtrackerRow& kr2m = GetRow((sec2+fkNOS)%fkNOS,imiddle);
3039 Int_t index1 = TMath::Max(kr2.Find(extraz-0.6-dddz1*TMath::Abs(z1)*0.05)-1,0);
3040 Int_t index2 = TMath::Min(kr2.Find(extraz+0.6+dddz2*TMath::Abs(z1)*0.05)+1,kr2);
3042 // rotation angles to p1-p3
3043 Double_t cs13r = cos(-angle13+dsec*alpha)/dvertex, sn13r = sin(-angle13+dsec*alpha)/dvertex;
3044 Double_t x2, y2, z2;
3046 // Double_t dymax = maxangle*TMath::Abs(x1-xx2);
3049 Double_t dxx0 = (xx2-x3)*cs13r;
3050 Double_t dyy0 = (xx2-x3)*sn13r;
3051 for (Int_t js=index1; js < index2; js++) {
3052 const AliTPCclusterMI *kcl = kr2[js];
3053 if (kcl->IsUsed(10)) continue;
3055 //calcutate parameters
3057 Double_t yy0 = dyy0 +(kcl->GetY()-y3)*cs13r;
3059 if (TMath::Abs(yy0)<0.000001) continue;
3060 Double_t xx0 = dxx0 -(kcl->GetY()-y3)*sn13r;
3061 Double_t y0 = 0.5*(xx0*xx0+yy0*yy0-xx0)/yy0;
3062 Double_t r02 = (0.25+y0*y0)*dvertex2;
3063 //curvature (radius) cut
3064 if (r02<r2min) continue;
3068 Double_t c0 = 1/TMath::Sqrt(r02);
3072 //Double_t dfi0 = 2.*TMath::ASin(dvertex*c0*0.5);
3073 //Double_t dfi1 = 2.*TMath::ASin(TMath::Sqrt(yy0*yy0+(1-xx0)*(1-xx0))*dvertex*c0*0.5);
3074 Double_t dfi0 = 2.*AliTPCFastMath::FastAsin(dvertex*c0*0.5);
3075 Double_t dfi1 = 2.*AliTPCFastMath::FastAsin(TMath::Sqrt(yy0*yy0+(1-xx0)*(1-xx0))*dvertex*c0*0.5);
3078 Double_t z0 = kcl->GetZ();
3079 Double_t zzzz2 = z1-(z1-z3)*dfi1/dfi0;
3080 if (TMath::Abs(zzzz2-z0)>0.5) continue;
3083 Double_t dip = (z1-z0)*c0/dfi1;
3084 Double_t x0 = (0.5*cs13+y0*sn13)*dvertex*c0;
3095 x2= xx2*cs-y2*sn*dsec;
3096 y2=+xx2*sn*dsec+y2*cs;
3106 // do we have cluster at the middle ?
3108 GetProlongation(x1,xm,x,ym,zm);
3110 AliTPCclusterMI * cm=0;
3111 if (TMath::Abs(ym)-ymaxm<0){
3112 cm = krm.FindNearest2(ym,zm,1.0,0.6,dummy);
3113 if ((!cm) || (cm->IsUsed(10))) {
3118 // rotate y1 to system 0
3119 // get state vector in rotated system
3120 Double_t yr1 = (-0.5*sn13+y0*cs13)*dvertex*c0;
3121 Double_t xr2 = x0*cs+yr1*sn*dsec;
3122 Double_t xr[5]={kcl->GetY(),kcl->GetZ(), xr2, dip, c0};
3124 GetProlongation(xx2,xm,xr,ym,zm);
3125 if (TMath::Abs(ym)-ymaxm<0){
3126 cm = kr2m.FindNearest2(ym,zm,1.0,0.6,dummy);
3127 if ((!cm) || (cm->IsUsed(10))) {
3137 dym = ym - cm->GetY();
3138 dzm = zm - cm->GetZ();
3145 Double_t sy1=kr1[is]->GetSigmaY2()*2., sz1=kr1[is]->GetSigmaZ2()*2.;
3146 Double_t sy2=kcl->GetSigmaY2()*2., sz2=kcl->GetSigmaZ2()*2.;
3147 //Double_t sy3=400*3./12., sy=0.1, sz=0.1;
3148 Double_t sy3=25000*x[4]*x[4]+0.1, sy=0.1, sz=0.1;
3149 //Double_t sy3=25000*x[4]*x[4]*60+0.5, sy=0.1, sz=0.1;
3151 Double_t f40=(F1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
3152 Double_t f42=(F1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
3153 Double_t f43=(F1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
3154 Double_t f20=(F2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
3155 Double_t f22=(F2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
3156 Double_t f23=(F2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
3158 Double_t f30=(F3(x1,y1+sy,x2,y2,z1,z2)-x[3])/sy;
3159 Double_t f31=(F3(x1,y1,x2,y2,z1+sz,z2)-x[3])/sz;
3160 Double_t f32=(F3(x1,y1,x2,y2+sy,z1,z2)-x[3])/sy;
3161 Double_t f34=(F3(x1,y1,x2,y2,z1,z2+sz)-x[3])/sz;
3165 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3166 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3167 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3168 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3169 c[13]=f30*sy1*f40+f32*sy2*f42;
3170 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3172 // if (!BuildSeed(kr1[is],kcl,0,x1,x2,x3,x,c)) continue;
3174 UInt_t index=kr1.GetIndex(is);
3175 seed->~AliTPCseed(); // this does not set the pointer to 0...
3176 AliTPCseed *track=new(seed) AliTPCseed(x1, ns*alpha+shift, x, c, index);
3178 track->SetIsSeeding(kTRUE);
3179 track->SetSeed1(i1);
3180 track->SetSeed2(i2);
3181 track->SetSeedType(3);
3185 FollowProlongation(*track, (i1+i2)/2,1);
3186 Int_t foundable,found,shared;
3187 track->GetClusterStatistic((i1+i2)/2,i1, found, foundable, shared, kTRUE);
3188 if ((found<0.55*foundable) || shared>0.5*found || (track->GetSigmaY2()+track->GetSigmaZ2())>0.5){
3190 seed->~AliTPCseed();
3196 FollowProlongation(*track, i2,1);
3200 track->SetBConstrain(1);
3201 // track->fLastPoint = i1+fInnerSec->GetNRows(); // first cluster in track position
3202 track->SetLastPoint(i1); // first cluster in track position
3203 track->SetFirstPoint(track->GetLastPoint());
3205 if (track->GetNumberOfClusters()<(i1-i2)*0.5 ||
3206 track->GetNumberOfClusters() < track->GetNFoundable()*0.6 ||
3207 track->GetNShared()>0.4*track->GetNumberOfClusters() ) {
3209 seed->~AliTPCseed();
3213 // Z VERTEX CONDITION
3214 Double_t zv, bz=GetBz();
3215 if ( !track->GetZAt(0.,bz,zv) ) continue;
3216 if (TMath::Abs(zv-z3)>cuts[2]) {
3217 FollowProlongation(*track, TMath::Max(i2-20,0));
3218 if ( !track->GetZAt(0.,bz,zv) ) continue;
3219 if (TMath::Abs(zv-z3)>cuts[2]){
3220 FollowProlongation(*track, TMath::Max(i2-40,0));
3221 if ( !track->GetZAt(0.,bz,zv) ) continue;
3222 if (TMath::Abs(zv-z3)>cuts[2] &&(track->GetNumberOfClusters() > track->GetNFoundable()*0.7)){
3223 // make seed without constrain
3224 AliTPCseed * track2 = MakeSeed(track,0.2,0.5,1.);
3225 FollowProlongation(*track2, i2,1);
3226 track2->SetBConstrain(kFALSE);
3227 track2->SetSeedType(1);
3228 arr->AddLast(track2);
3230 seed->~AliTPCseed();
3235 seed->~AliTPCseed();
3242 track->SetSeedType(0);
3243 arr->AddLast(track);
3244 seed = new AliTPCseed;
3246 // don't consider other combinations
3247 if (track->GetNumberOfClusters() > track->GetNFoundable()*0.8)
3253 Info("MakeSeeds3","\nSeeding statistic:\t%d\t%d\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin,nout1,nout2);
3259 void AliTPCtrackerMI::MakeSeeds5(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, Float_t cuts[4],
3264 //-----------------------------------------------------------------
3265 // This function creates track seeds.
3266 //-----------------------------------------------------------------
3267 // cuts[0] - fP4 cut
3268 // cuts[1] - tan(phi) cut
3269 // cuts[2] - zvertex cut
3270 // cuts[3] - fP3 cut
3280 Double_t x[5], c[15];
3282 // make temporary seed
3283 AliTPCseed * seed = new AliTPCseed;
3284 Double_t alpha=fOuterSec->GetAlpha(), shift=fOuterSec->GetAlphaShift();
3285 // Double_t cs=cos(alpha), sn=sin(alpha);
3290 Double_t x1 = GetXrow(i1-1);
3291 const AliTPCtrackerRow& kr1=GetRow(sec,i1-1);
3292 Double_t y1max = GetMaxY(i1-1)-kr1.GetDeadZone()-1.5;
3294 Double_t x1p = GetXrow(i1);
3295 const AliTPCtrackerRow& kr1p=GetRow(sec,i1);
3297 Double_t x1m = GetXrow(i1-2);
3298 const AliTPCtrackerRow& kr1m=GetRow(sec,i1-2);
3301 //last 3 padrow for seeding
3302 AliTPCtrackerRow& kr3 = GetRow((sec+fkNOS)%fkNOS,i1-7);
3303 Double_t x3 = GetXrow(i1-7);
3304 // Double_t y3max= GetMaxY(i1-7)-kr3.fDeadZone-1.5;
3306 AliTPCtrackerRow& kr3p = GetRow((sec+fkNOS)%fkNOS,i1-6);
3307 Double_t x3p = GetXrow(i1-6);
3309 AliTPCtrackerRow& kr3m = GetRow((sec+fkNOS)%fkNOS,i1-8);
3310 Double_t x3m = GetXrow(i1-8);
3315 Int_t im = i1-4; //middle pad row index
3316 Double_t xm = GetXrow(im); // radius of middle pad-row
3317 const AliTPCtrackerRow& krm=GetRow(sec,im); //middle pad -row
3318 // Double_t ymmax = GetMaxY(im)-kr1.fDeadZone-1.5;
3321 Double_t deltax = x1-x3;
3322 Double_t dymax = deltax*cuts[1];
3323 Double_t dzmax = deltax*cuts[3];
3325 // loop over clusters
3326 for (Int_t is=0; is < kr1; is++) {
3328 if (kr1[is]->IsUsed(10)) continue;
3329 Double_t y1=kr1[is]->GetY(), z1=kr1[is]->GetZ();
3331 if (deltay>0 && TMath::Abs(y1max-TMath::Abs(y1))> deltay ) continue; // seed only at the edge
3333 Int_t index1 = TMath::Max(kr3.Find(z1-dzmax)-1,0);
3334 Int_t index2 = TMath::Min(kr3.Find(z1+dzmax)+1,kr3);
3340 for (Int_t js=index1; js < index2; js++) {
3341 const AliTPCclusterMI *kcl = kr3[js];
3342 if (kcl->IsUsed(10)) continue;
3344 // apply angular cuts
3345 if (TMath::Abs(y1-y3)>dymax) continue;
3348 if (TMath::Abs(z1-z3)>dzmax) continue;
3350 Double_t angley = (y1-y3)/(x1-x3);
3351 Double_t anglez = (z1-z3)/(x1-x3);
3353 Double_t erry = TMath::Abs(angley)*(x1-x1m)*0.5+0.5;
3354 Double_t errz = TMath::Abs(anglez)*(x1-x1m)*0.5+0.5;
3356 Double_t yyym = angley*(xm-x1)+y1;
3357 Double_t zzzm = anglez*(xm-x1)+z1;
3359 const AliTPCclusterMI *kcm = krm.FindNearest2(yyym,zzzm,erry,errz,index);
3361 if (kcm->IsUsed(10)) continue;
3363 erry = TMath::Abs(angley)*(x1-x1m)*0.4+0.5;
3364 errz = TMath::Abs(anglez)*(x1-x1m)*0.4+0.5;
3371 // look around first
3372 const AliTPCclusterMI *kc1m = kr1m.FindNearest2(angley*(x1m-x1)+y1,
3378 if (kc1m->IsUsed(10)) used++;
3380 const AliTPCclusterMI *kc1p = kr1p.FindNearest2(angley*(x1p-x1)+y1,
3386 if (kc1p->IsUsed(10)) used++;
3388 if (used>1) continue;
3389 if (found<1) continue;
3393 const AliTPCclusterMI *kc3m = kr3m.FindNearest2(angley*(x3m-x3)+y3,
3399 if (kc3m->IsUsed(10)) used++;
3403 const AliTPCclusterMI *kc3p = kr3p.FindNearest2(angley*(x3p-x3)+y3,
3409 if (kc3p->IsUsed(10)) used++;
3413 if (used>1) continue;
3414 if (found<3) continue;
3424 x[4]=F1(x1,y1,x2,y2,x3,y3);
3425 //if (TMath::Abs(x[4]) >= cuts[0]) continue;
3428 x[2]=F2(x1,y1,x2,y2,x3,y3);
3431 x[3]=F3n(x1,y1,x2,y2,z1,z2,x[4]);
3432 //if (TMath::Abs(x[3]) > cuts[3]) continue;
3436 Double_t sy1=0.1, sz1=0.1;
3437 Double_t sy2=0.1, sz2=0.1;
3438 Double_t sy3=0.1, sy=0.1, sz=0.1;
3440 Double_t f40=(F1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
3441 Double_t f42=(F1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
3442 Double_t f43=(F1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
3443 Double_t f20=(F2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
3444 Double_t f22=(F2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
3445 Double_t f23=(F2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
3447 Double_t f30=(F3(x1,y1+sy,x2,y2,z1,z2)-x[3])/sy;
3448 Double_t f31=(F3(x1,y1,x2,y2,z1+sz,z2)-x[3])/sz;
3449 Double_t f32=(F3(x1,y1,x2,y2+sy,z1,z2)-x[3])/sy;
3450 Double_t f34=(F3(x1,y1,x2,y2,z1,z2+sz)-x[3])/sz;
3454 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3455 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3456 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3457 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3458 c[13]=f30*sy1*f40+f32*sy2*f42;
3459 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3461 // if (!BuildSeed(kr1[is],kcl,0,x1,x2,x3,x,c)) continue;
3463 index=kr1.GetIndex(is);
3464 seed->~AliTPCseed();
3465 AliTPCseed *track=new(seed) AliTPCseed(x1, sec*alpha+shift, x, c, index);
3467 track->SetIsSeeding(kTRUE);
3470 FollowProlongation(*track, i1-7,1);
3471 if (track->GetNumberOfClusters() < track->GetNFoundable()*0.75 ||
3472 track->GetNShared()>0.6*track->GetNumberOfClusters() || ( track->GetSigmaY2()+ track->GetSigmaZ2())>0.6){
3474 seed->~AliTPCseed();
3480 FollowProlongation(*track, i2,1);
3481 track->SetBConstrain(0);
3482 track->SetLastPoint(i1+fInnerSec->GetNRows()); // first cluster in track position
3483 track->SetFirstPoint(track->GetLastPoint());
3485 if (track->GetNumberOfClusters()<(i1-i2)*0.5 ||
3486 track->GetNumberOfClusters()<track->GetNFoundable()*0.7 ||
3487 track->GetNShared()>2. || track->GetChi2()/track->GetNumberOfClusters()>6 || ( track->GetSigmaY2()+ track->GetSigmaZ2())>0.5 ) {
3489 seed->~AliTPCseed();
3494 FollowProlongation(*track, TMath::Max(i2-10,0),1);
3495 AliTPCseed * track2 = MakeSeed(track,0.2,0.5,0.9);
3496 FollowProlongation(*track2, i2,1);
3497 track2->SetBConstrain(kFALSE);
3498 track2->SetSeedType(4);
3499 arr->AddLast(track2);
3501 seed->~AliTPCseed();
3505 //arr->AddLast(track);
3506 //seed = new AliTPCseed;
3512 Info("MakeSeeds5","\nSeeding statiistic:\t%d\t%d\t%d\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin,nout1,nout2,nout3);
3518 //_____________________________________________________________________________
3519 void AliTPCtrackerMI::MakeSeeds2(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, Float_t */*cuts[4]*/,
3520 Float_t deltay, Bool_t /*bconstrain*/) {
3521 //-----------------------------------------------------------------
3522 // This function creates track seeds - without vertex constraint
3523 //-----------------------------------------------------------------
3524 // cuts[0] - fP4 cut - not applied
3525 // cuts[1] - tan(phi) cut
3526 // cuts[2] - zvertex cut - not applied
3527 // cuts[3] - fP3 cut
3537 Double_t alpha=fOuterSec->GetAlpha(), shift=fOuterSec->GetAlphaShift();
3538 // Double_t cs=cos(alpha), sn=sin(alpha);
3539 Int_t row0 = (i1+i2)/2;
3540 Int_t drow = (i1-i2)/2;
3541 const AliTPCtrackerRow& kr0=fSectors[sec][row0];
3542 AliTPCtrackerRow * kr=0;
3544 AliTPCpolyTrack polytrack;
3545 Int_t nclusters=fSectors[sec][row0];
3546 AliTPCseed * seed = new AliTPCseed;
3551 for (Int_t is=0; is < nclusters; is++) { //LOOP over clusters
3553 Int_t nfoundable =0;
3554 for (Int_t iter =1; iter<2; iter++){ //iterations
3555 const AliTPCtrackerRow& krm=fSectors[sec][row0-iter];
3556 const AliTPCtrackerRow& krp=fSectors[sec][row0+iter];
3557 const AliTPCclusterMI * cl= kr0[is];
3559 if (cl->IsUsed(10)) {
3565 Double_t x = kr0.GetX();
3566 // Initialization of the polytrack
3571 Double_t y0= cl->GetY();
3572 Double_t z0= cl->GetZ();
3576 Double_t ymax = fSectors->GetMaxY(row0)-kr0.GetDeadZone()-1.5;
3577 if (deltay>0 && TMath::Abs(ymax-TMath::Abs(y0))> deltay ) continue; // seed only at the edge
3579 erry = (0.5)*cl->GetSigmaY2()/TMath::Sqrt(cl->GetQ())*6;
3580 errz = (0.5)*cl->GetSigmaZ2()/TMath::Sqrt(cl->GetQ())*6;
3581 polytrack.AddPoint(x,y0,z0,erry, errz);
3584 if (cl->IsUsed(10)) sumused++;
3587 Float_t roady = (5*TMath::Sqrt(cl->GetSigmaY2()+0.2)+1.)*iter;
3588 Float_t roadz = (5*TMath::Sqrt(cl->GetSigmaZ2()+0.2)+1.)*iter;
3591 AliTPCclusterMI * cl1 = krm.FindNearest(y0,z0,roady,roadz);
3592 if (cl1 && TMath::Abs(ymax-TMath::Abs(y0))) {
3593 erry = (0.5)*cl1->GetSigmaY2()/TMath::Sqrt(cl1->GetQ())*3;
3594 errz = (0.5)*cl1->GetSigmaZ2()/TMath::Sqrt(cl1->GetQ())*3;
3595 if (cl1->IsUsed(10)) sumused++;
3596 polytrack.AddPoint(x,cl1->GetY(),cl1->GetZ(),erry,errz);
3600 AliTPCclusterMI * cl2 = krp.FindNearest(y0,z0,roady,roadz);
3602 erry = (0.5)*cl2->GetSigmaY2()/TMath::Sqrt(cl2->GetQ())*3;
3603 errz = (0.5)*cl2->GetSigmaZ2()/TMath::Sqrt(cl2->GetQ())*3;
3604 if (cl2->IsUsed(10)) sumused++;
3605 polytrack.AddPoint(x,cl2->GetY(),cl2->GetZ(),erry,errz);
3608 if (sumused>0) continue;
3610 polytrack.UpdateParameters();
3616 nfoundable = polytrack.GetN();
3617 nfound = nfoundable;
3619 for (Int_t ddrow = iter+1; ddrow<drow;ddrow++){
3620 Float_t maxdist = 0.8*(1.+3./(ddrow));
3621 for (Int_t delta = -1;delta<=1;delta+=2){
3622 Int_t row = row0+ddrow*delta;
3623 kr = &(fSectors[sec][row]);
3624 Double_t xn = kr->GetX();
3625 Double_t ymax1 = fSectors->GetMaxY(row)-kr->GetDeadZone()-1.5;
3626 polytrack.GetFitPoint(xn,yn,zn);
3627 if (TMath::Abs(yn)>ymax1) continue;
3629 AliTPCclusterMI * cln = kr->FindNearest(yn,zn,roady,roadz);
3631 Float_t dist = TMath::Sqrt( (yn-cln->GetY())*(yn-cln->GetY())+(zn-cln->GetZ())*(zn-cln->GetZ()));
3634 erry = (dist+0.3)*cln->GetSigmaY2()/TMath::Sqrt(cln->GetQ())*(1.+1./(ddrow));
3635 errz = (dist+0.3)*cln->GetSigmaZ2()/TMath::Sqrt(cln->GetQ())*(1.+1./(ddrow));
3636 if (cln->IsUsed(10)) {
3637 // printf("used\n");
3645 polytrack.AddPoint(xn,cln->GetY(),cln->GetZ(),erry, errz);
3650 if ( (sumused>3) || (sumused>0.5*nfound) || (nfound<0.6*nfoundable)) break;
3651 polytrack.UpdateParameters();
3654 if ( (sumused>3) || (sumused>0.5*nfound)) {
3655 //printf("sumused %d\n",sumused);
3660 polytrack.GetFitDerivation(kr0.GetX(),dy,dz);
3661 AliTPCpolyTrack track2;
3663 polytrack.Refit(track2,0.5+TMath::Abs(dy)*0.3,0.4+TMath::Abs(dz)*0.3);
3664 if (track2.GetN()<0.5*nfoundable) continue;
3667 if ((nfound>0.6*nfoundable) &&( nfoundable>0.4*(i1-i2))) {
3669 // test seed with and without constrain
3670 for (Int_t constrain=0; constrain<=0;constrain++){
3671 // add polytrack candidate
3673 Double_t x[5], c[15];
3674 Double_t x1,x2,x3,y1,y2,y3,z1,z2,z3;
3675 track2.GetBoundaries(x3,x1);
3677 track2.GetFitPoint(x1,y1,z1);
3678 track2.GetFitPoint(x2,y2,z2);
3679 track2.GetFitPoint(x3,y3,z3);
3681 //is track pointing to the vertex ?
3684 polytrack.GetFitPoint(x0,y0,z0);
3697 x[4]=F1(x1,y1,x2,y2,x3,y3);
3699 // if (TMath::Abs(x[4]) >= cuts[0]) continue; //
3700 x[2]=F2(x1,y1,x2,y2,x3,y3);
3702 //if (TMath::Abs(x[4]*x1-x[2]) >= cuts[1]) continue;
3703 //x[3]=F3(x1,y1,x2,y2,z1,z2);
3704 x[3]=F3n(x1,y1,x3,y3,z1,z3,x[4]);
3705 //if (TMath::Abs(x[3]) > cuts[3]) continue;
3708 Double_t sy =0.1, sz =0.1;
3709 Double_t sy1=0.02, sz1=0.02;
3710 Double_t sy2=0.02, sz2=0.02;
3714 sy3=25000*x[4]*x[4]+0.1, sy=0.1, sz=0.1;
3717 Double_t f40=(F1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
3718 Double_t f42=(F1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
3719 Double_t f43=(F1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
3720 Double_t f20=(F2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
3721 Double_t f22=(F2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
3722 Double_t f23=(F2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
3724 Double_t f30=(F3(x1,y1+sy,x3,y3,z1,z3)-x[3])/sy;
3725 Double_t f31=(F3(x1,y1,x3,y3,z1+sz,z3)-x[3])/sz;
3726 Double_t f32=(F3(x1,y1,x3,y3+sy,z1,z3)-x[3])/sy;
3727 Double_t f34=(F3(x1,y1,x3,y3,z1,z3+sz)-x[3])/sz;
3732 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3733 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3734 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3735 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3736 c[13]=f30*sy1*f40+f32*sy2*f42;
3737 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3739 //Int_t row1 = fSectors->GetRowNumber(x1);
3740 Int_t row1 = GetRowNumber(x1);
3744 seed->~AliTPCseed();
3745 AliTPCseed *track=new(seed) AliTPCseed(x1,sec*alpha+shift,x,c,index);
3746 track->SetIsSeeding(kTRUE);
3747 Int_t rc=FollowProlongation(*track, i2);
3748 if (constrain) track->SetBConstrain(1);
3750 track->SetBConstrain(0);
3751 track->SetLastPoint(row1+fInnerSec->GetNRows()); // first cluster in track position
3752 track->SetFirstPoint(track->GetLastPoint());
3754 if (rc==0 || track->GetNumberOfClusters()<(i1-i2)*0.5 ||
3755 track->GetNumberOfClusters() < track->GetNFoundable()*0.6 ||
3756 track->GetNShared()>0.4*track->GetNumberOfClusters()) {
3759 seed->~AliTPCseed();
3762 arr->AddLast(track);
3763 seed = new AliTPCseed;
3767 } // if accepted seed
3770 Info("MakeSeeds2","\nSeeding statiistic:\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin3);
3776 AliTPCseed *AliTPCtrackerMI::MakeSeed(AliTPCseed *const track, Float_t r0, Float_t r1, Float_t r2)
3780 //reseed using track points
3781 Int_t p0 = int(r0*track->GetNumberOfClusters()); // point 0
3782 Int_t p1 = int(r1*track->GetNumberOfClusters());
3783 Int_t p2 = int(r2*track->GetNumberOfClusters()); // last point
3785 Double_t x0[3],x1[3],x2[3];
3786 for (Int_t i=0;i<3;i++){
3792 // find track position at given ratio of the length
3793 Int_t sec0=0, sec1=0, sec2=0;
3796 for (Int_t i=0;i<160;i++){
3797 if (track->GetClusterPointer(i)){
3799 AliTPCTrackerPoint *trpoint =track->GetTrackPoint(i);
3800 if ( (index<p0) || x0[0]<0 ){
3801 if (trpoint->GetX()>1){
3802 clindex = track->GetClusterIndex2(i);
3804 x0[0] = trpoint->GetX();
3805 x0[1] = trpoint->GetY();
3806 x0[2] = trpoint->GetZ();
3807 sec0 = ((clindex&0xff000000)>>24)%18;
3812 if ( (index<p1) &&(trpoint->GetX()>1)){
3813 clindex = track->GetClusterIndex2(i);
3815 x1[0] = trpoint->GetX();
3816 x1[1] = trpoint->GetY();
3817 x1[2] = trpoint->GetZ();
3818 sec1 = ((clindex&0xff000000)>>24)%18;
3821 if ( (index<p2) &&(trpoint->GetX()>1)){
3822 clindex = track->GetClusterIndex2(i);
3824 x2[0] = trpoint->GetX();
3825 x2[1] = trpoint->GetY();
3826 x2[2] = trpoint->GetZ();
3827 sec2 = ((clindex&0xff000000)>>24)%18;
3834 Double_t alpha, cs,sn, xx2,yy2;
3836 alpha = (sec1-sec2)*fSectors->GetAlpha();
3837 cs = TMath::Cos(alpha);
3838 sn = TMath::Sin(alpha);
3839 xx2= x1[0]*cs-x1[1]*sn;
3840 yy2= x1[0]*sn+x1[1]*cs;
3844 alpha = (sec0-sec2)*fSectors->GetAlpha();
3845 cs = TMath::Cos(alpha);
3846 sn = TMath::Sin(alpha);
3847 xx2= x0[0]*cs-x0[1]*sn;
3848 yy2= x0[0]*sn+x0[1]*cs;
3854 Double_t x[5],c[15];
3858 x[4]=F1(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]);
3859 // if (x[4]>1) return 0;
3860 x[2]=F2(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]);
3861 x[3]=F3n(x2[0],x2[1],x0[0],x0[1],x2[2],x0[2],x[4]);
3862 //if (TMath::Abs(x[3]) > 2.2) return 0;
3863 //if (TMath::Abs(x[2]) > 1.99) return 0;
3865 Double_t sy =0.1, sz =0.1;
3867 Double_t sy1=0.02+track->GetSigmaY2(), sz1=0.02+track->GetSigmaZ2();
3868 Double_t sy2=0.01+track->GetSigmaY2(), sz2=0.01+track->GetSigmaZ2();
3869 Double_t sy3=0.01+track->GetSigmaY2();
3871 Double_t f40=(F1(x2[0],x2[1]+sy,x1[0],x1[1],x0[0],x0[1])-x[4])/sy;
3872 Double_t f42=(F1(x2[0],x2[1],x1[0],x1[1]+sy,x0[0],x0[1])-x[4])/sy;
3873 Double_t f43=(F1(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]+sy)-x[4])/sy;
3874 Double_t f20=(F2(x2[0],x2[1]+sy,x1[0],x1[1],x0[0],x0[1])-x[2])/sy;
3875 Double_t f22=(F2(x2[0],x2[1],x1[0],x1[1]+sy,x0[0],x0[1])-x[2])/sy;
3876 Double_t f23=(F2(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]+sy)-x[2])/sy;
3878 Double_t f30=(F3(x2[0],x2[1]+sy,x0[0],x0[1],x2[2],x0[2])-x[3])/sy;
3879 Double_t f31=(F3(x2[0],x2[1],x0[0],x0[1],x2[2]+sz,x0[2])-x[3])/sz;
3880 Double_t f32=(F3(x2[0],x2[1],x0[0],x0[1]+sy,x2[2],x0[2])-x[3])/sy;
3881 Double_t f34=(F3(x2[0],x2[1],x0[0],x0[1],x2[2],x0[2]+sz)-x[3])/sz;
3886 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3887 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3888 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3889 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3890 c[13]=f30*sy1*f40+f32*sy2*f42;
3891 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3893 // Int_t row1 = fSectors->GetRowNumber(x2[0]);
3894 AliTPCseed *seed=new AliTPCseed(x2[0], sec2*fSectors->GetAlpha()+fSectors->GetAlphaShift(), x, c, 0);
3895 // Double_t y0,z0,y1,z1, y2,z2;
3896 //seed->GetProlongation(x0[0],y0,z0);
3897 // seed->GetProlongation(x1[0],y1,z1);
3898 //seed->GetProlongation(x2[0],y2,z2);
3900 seed->SetLastPoint(pp2);
3901 seed->SetFirstPoint(pp2);
3908 AliTPCseed *AliTPCtrackerMI::ReSeed(const AliTPCseed *track, Float_t r0, Float_t r1, Float_t r2)
3912 //reseed using founded clusters
3914 // Find the number of clusters
3915 Int_t nclusters = 0;
3916 for (Int_t irow=0;irow<160;irow++){
3917 if (track->GetClusterIndex(irow)>0) nclusters++;
3921 ipos[0] = TMath::Max(int(r0*nclusters),0); // point 0 cluster
3922 ipos[1] = TMath::Min(int(r1*nclusters),nclusters-1); //
3923 ipos[2] = TMath::Min(int(r2*nclusters),nclusters-1); // last point
3927 Int_t row[3],sec[3]={0,0,0};
3929 // find track row position at given ratio of the length
3931 for (Int_t irow=0;irow<160;irow++){
3932 if (track->GetClusterIndex2(irow)<0) continue;
3934 for (Int_t ipoint=0;ipoint<3;ipoint++){
3935 if (index<=ipos[ipoint]) row[ipoint] = irow;
3939 //Get cluster and sector position
3940 for (Int_t ipoint=0;ipoint<3;ipoint++){
3941 Int_t clindex = track->GetClusterIndex2(row[ipoint]);
3942 AliTPCclusterMI * cl = GetClusterMI(clindex);
3945 // AliTPCclusterMI * cl = GetClusterMI(clindex);
3948 sec[ipoint] = ((clindex&0xff000000)>>24)%18;
3949 xyz[ipoint][0] = GetXrow(row[ipoint]);
3950 xyz[ipoint][1] = cl->GetY();
3951 xyz[ipoint][2] = cl->GetZ();
3955 // Calculate seed state vector and covariance matrix
3957 Double_t alpha, cs,sn, xx2,yy2;
3959 alpha = (sec[1]-sec[2])*fSectors->GetAlpha();
3960 cs = TMath::Cos(alpha);
3961 sn = TMath::Sin(alpha);
3962 xx2= xyz[1][0]*cs-xyz[1][1]*sn;
3963 yy2= xyz[1][0]*sn+xyz[1][1]*cs;
3967 alpha = (sec[0]-sec[2])*fSectors->GetAlpha();
3968 cs = TMath::Cos(alpha);
3969 sn = TMath::Sin(alpha);
3970 xx2= xyz[0][0]*cs-xyz[0][1]*sn;
3971 yy2= xyz[0][0]*sn+xyz[0][1]*cs;
3977 Double_t x[5],c[15];
3981 x[4]=F1(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
3982 x[2]=F2(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
3983 x[3]=F3n(xyz[2][0],xyz[2][1],xyz[0][0],xyz[0][1],xyz[2][2],xyz[0][2],x[4]);
3985 Double_t sy =0.1, sz =0.1;
3987 Double_t sy1=0.2, sz1=0.2;
3988 Double_t sy2=0.2, sz2=0.2;
3991 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;
3992 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;
3993 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;
3994 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;
3995 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;
3996 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;
3998 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;
3999 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;
4000 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;
4001 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;
4006 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
4007 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
4008 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
4009 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
4010 c[13]=f30*sy1*f40+f32*sy2*f42;
4011 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
4013 // Int_t row1 = fSectors->GetRowNumber(xyz[2][0]);
4014 AliTPCseed *seed=new AliTPCseed(xyz[2][0], sec[2]*fSectors->GetAlpha()+fSectors->GetAlphaShift(), x, c, 0);
4015 seed->SetLastPoint(row[2]);
4016 seed->SetFirstPoint(row[2]);
4021 AliTPCseed *AliTPCtrackerMI::ReSeed(AliTPCseed *track,Int_t r0, Bool_t forward)
4025 //reseed using founded clusters
4028 Int_t row[3]={0,0,0};
4029 Int_t sec[3]={0,0,0};
4031 // forward direction
4033 for (Int_t irow=r0;irow<160;irow++){
4034 if (track->GetClusterIndex(irow)>0){
4039 for (Int_t irow=160;irow>r0;irow--){
4040 if (track->GetClusterIndex(irow)>0){
4045 for (Int_t irow=row[2]-15;irow>row[0];irow--){
4046 if (track->GetClusterIndex(irow)>0){
4054 for (Int_t irow=0;irow<r0;irow++){
4055 if (track->GetClusterIndex(irow)>0){
4060 for (Int_t irow=r0;irow>0;irow--){
4061 if (track->GetClusterIndex(irow)>0){
4066 for (Int_t irow=row[2]-15;irow>row[0];irow--){
4067 if (track->GetClusterIndex(irow)>0){
4074 if ((row[2]-row[0])<20) return 0;
4075 if (row[1]==0) return 0;
4078 //Get cluster and sector position
4079 for (Int_t ipoint=0;ipoint<3;ipoint++){
4080 Int_t clindex = track->GetClusterIndex2(row[ipoint]);
4081 AliTPCclusterMI * cl = GetClusterMI(clindex);
4084 // AliTPCclusterMI * cl = GetClusterMI(clindex);
4087 sec[ipoint] = ((clindex&0xff000000)>>24)%18;
4088 xyz[ipoint][0] = GetXrow(row[ipoint]);
4089 AliTPCTrackerPoint * point = track->GetTrackPoint(row[ipoint]);
4090 if (point&&ipoint<2){
4092 xyz[ipoint][1] = point->GetY();
4093 xyz[ipoint][2] = point->GetZ();
4096 xyz[ipoint][1] = cl->GetY();
4097 xyz[ipoint][2] = cl->GetZ();
4104 // Calculate seed state vector and covariance matrix
4106 Double_t alpha, cs,sn, xx2,yy2;
4108 alpha = (sec[1]-sec[2])*fSectors->GetAlpha();
4109 cs = TMath::Cos(alpha);
4110 sn = TMath::Sin(alpha);
4111 xx2= xyz[1][0]*cs-xyz[1][1]*sn;
4112 yy2= xyz[1][0]*sn+xyz[1][1]*cs;
4116 alpha = (sec[0]-sec[2])*fSectors->GetAlpha();
4117 cs = TMath::Cos(alpha);
4118 sn = TMath::Sin(alpha);
4119 xx2= xyz[0][0]*cs-xyz[0][1]*sn;
4120 yy2= xyz[0][0]*sn+xyz[0][1]*cs;
4126 Double_t x[5],c[15];
4130 x[4]=F1(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
4131 x[2]=F2(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
4132 x[3]=F3n(xyz[2][0],xyz[2][1],xyz[0][0],xyz[0][1],xyz[2][2],xyz[0][2],x[4]);
4134 Double_t sy =0.1, sz =0.1;
4136 Double_t sy1=0.2, sz1=0.2;
4137 Double_t sy2=0.2, sz2=0.2;
4140 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;
4141 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;
4142 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;
4143 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;
4144 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;
4145 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;
4147 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;
4148 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;
4149 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;
4150 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;
4155 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
4156 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
4157 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
4158 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
4159 c[13]=f30*sy1*f40+f32*sy2*f42;
4160 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
4162 // Int_t row1 = fSectors->GetRowNumber(xyz[2][0]);
4163 AliTPCseed *seed=new AliTPCseed(xyz[2][0], sec[2]*fSectors->GetAlpha()+fSectors->GetAlphaShift(), x, c, 0);
4164 seed->SetLastPoint(row[2]);
4165 seed->SetFirstPoint(row[2]);
4166 for (Int_t i=row[0];i<row[2];i++){
4167 seed->SetClusterIndex(i, track->GetClusterIndex(i));
4175 void AliTPCtrackerMI::FindMultiMC(const TObjArray * array, AliESDEvent */*esd*/, Int_t iter)
4178 // find multi tracks - THIS FUNCTION IS ONLY FOR DEBUG PURPOSES
4180 // Use AliTPCReconstructor::StreamLevel()>2 if you want to tune parameters - cuts
4182 // Two reasons to have multiple find tracks
4183 // 1. Curling tracks can be find more than once
4184 // 2. Splitted tracks
4185 // a.) Multiple seeding to increase tracking efficiency - (~ 100% reached)
4186 // b.) Edge effect on the sector boundaries
4189 // Algorithm done in 2 phases - because of CPU consumption
4190 // it is n^2 algorithm - for lead-lead 20000x20000 combination are investigated
4192 // Algorihm for curling tracks sign:
4193 // 1 phase -makes a very rough fast cuts to minimize combinatorics
4194 // a.) opposite sign
4195 // b.) one of the tracks - not pointing to the primary vertex -
4196 // c.) delta tan(theta)
4198 // 2 phase - calculates DCA between tracks - time consument
4203 // General cuts - for splitted tracks and for curling tracks
4205 const Float_t kMaxdPhi = 0.2; // maximal distance in phi
4207 // Curling tracks cuts
4212 Int_t nentries = array->GetEntriesFast();
4213 AliHelix *helixes = new AliHelix[nentries];
4214 Float_t *xm = new Float_t[nentries];
4215 Float_t *dz0 = new Float_t[nentries];
4216 Float_t *dz1 = new Float_t[nentries];
4222 // Find track COG in x direction - point with best defined parameters
4224 for (Int_t i=0;i<nentries;i++){
4225 AliTPCseed* track = (AliTPCseed*)array->At(i);
4226 if (!track) continue;
4227 track->SetCircular(0);
4228 new (&helixes[i]) AliHelix(*track);
4232 track->GetDZ(GetX(),GetY(),GetZ(),GetBz(),dz);
4235 for (Int_t icl=0; icl<160; icl++){
4236 AliTPCclusterMI * cl = track->GetClusterPointer(icl);
4242 if (ncl>0) xm[i]/=Float_t(ncl);
4245 for (Int_t i0=0;i0<nentries;i0++){
4246 AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
4247 if (!track0) continue;
4248 Float_t xc0 = helixes[i0].GetHelix(6);
4249 Float_t yc0 = helixes[i0].GetHelix(7);
4250 Float_t r0 = helixes[i0].GetHelix(8);
4251 Float_t rc0 = TMath::Sqrt(xc0*xc0+yc0*yc0);
4252 Float_t fi0 = TMath::ATan2(yc0,xc0);
4254 for (Int_t i1=i0+1;i1<nentries;i1++){
4255 AliTPCseed * track1 = (AliTPCseed*)array->At(i1);
4256 if (!track1) continue;
4257 Int_t lab0=track0->GetLabel();
4258 Int_t lab1=track1->GetLabel();
4259 if (TMath::Abs(lab0)!=TMath::Abs(lab1)) continue;
4261 Float_t xc1 = helixes[i1].GetHelix(6);
4262 Float_t yc1 = helixes[i1].GetHelix(7);
4263 Float_t r1 = helixes[i1].GetHelix(8);
4264 Float_t rc1 = TMath::Sqrt(xc1*xc1+yc1*yc1);
4265 Float_t fi1 = TMath::ATan2(yc1,xc1);
4267 Float_t dfi = fi0-fi1;
4270 if (dfi>1.5*TMath::Pi()) dfi-=TMath::Pi(); // take care about edge effect
4271 if (dfi<-1.5*TMath::Pi()) dfi+=TMath::Pi(); //
4272 if (TMath::Abs(dfi)>kMaxdPhi&&helixes[i0].GetHelix(4)*helixes[i1].GetHelix(4)<0){
4274 // if short tracks with undefined sign
4275 fi1 = -TMath::ATan2(yc1,-xc1);
4278 Float_t dtheta = TMath::Abs(track0->GetTgl()-track1->GetTgl())<TMath::Abs(track0->GetTgl()+track1->GetTgl())? track0->GetTgl()-track1->GetTgl():track0->GetTgl()+track1->GetTgl();
4281 // debug stream to tune "fast cuts"
4283 Double_t dist[3]; // distance at X
4284 Double_t mdist[3]={0,0,0}; // mean distance X+-40cm
4285 track0->GetDistance(track1,0.5*(xm[i0]+xm[i1])-40.,dist,AliTracker::GetBz());
4286 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[i]);
4287 track0->GetDistance(track1,0.5*(xm[i0]+xm[i1])+40.,dist,AliTracker::GetBz());
4288 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[i]);
4289 track0->GetDistance(track1,0.5*(xm[i0]+xm[i1]),dist,AliTracker::GetBz());
4290 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[i]);
4291 for (Int_t i=0;i<3;i++) mdist[i]*=0.33333;
4295 for (Int_t icl=0; icl<160; icl++){
4296 AliTPCclusterMI * cl0 = track0->GetClusterPointer(icl);
4297 AliTPCclusterMI * cl1 = track1->GetClusterPointer(icl);
4300 if (cl0==cl1) sums++;
4304 if (AliTPCReconstructor::StreamLevel()>5) {
4305 TTreeSRedirector &cstream = *fDebugStreamer;
4310 "Tr0.="<<track0<< // seed0
4311 "Tr1.="<<track1<< // seed1
4312 "h0.="<<&helixes[i0]<<
4313 "h1.="<<&helixes[i1]<<
4315 "sum="<<sum<< //the sum of rows with cl in both
4316 "sums="<<sums<< //the sum of shared clusters
4317 "xm0="<<xm[i0]<< // the center of track
4318 "xm1="<<xm[i1]<< // the x center of track
4319 // General cut variables
4320 "dfi="<<dfi<< // distance in fi angle
4321 "dtheta="<<dtheta<< // distance int theta angle
4327 "dist0="<<dist[0]<< //distance x
4328 "dist1="<<dist[1]<< //distance y
4329 "dist2="<<dist[2]<< //distance z
4330 "mdist0="<<mdist[0]<< //distance x
4331 "mdist1="<<mdist[1]<< //distance y
4332 "mdist2="<<mdist[2]<< //distance z
4346 if (AliTPCReconstructor::StreamLevel()>1) {
4347 AliInfo("Time for curling tracks removal DEBUGGING MC");
4354 void AliTPCtrackerMI::FindSplitted(TObjArray * array, AliESDEvent */*esd*/, Int_t /*iter*/){
4356 // Find Splitted tracks and remove the one with worst quality
4357 // Corresponding debug streamer to tune selections - "Splitted2"
4359 // 0. Sort tracks according quility
4360 // 1. Propagate the tracks to the reference radius
4361 // 2. Double_t loop to select close tracks (only to speed up process)
4362 // 3. Calculate cluster overlap ratio - and remove the track if bigger than a threshold
4363 // 4. Delete temporary parameters
4365 const Double_t xref=GetXrow(63); // reference radius -IROC/OROC boundary
4367 const Double_t kCutP1=10; // delta Z cut 10 cm
4368 const Double_t kCutP2=0.15; // delta snp(fi) cut 0.15
4369 const Double_t kCutP3=0.15; // delta tgl(theta) cut 0.15
4370 const Double_t kCutAlpha=0.15; // delta alpha cut
4371 Int_t firstpoint = 0;
4372 Int_t lastpoint = 160;
4374 Int_t nentries = array->GetEntriesFast();
4375 AliExternalTrackParam *params = new AliExternalTrackParam[nentries];
4381 //0. Sort tracks according quality
4382 //1. Propagate the ext. param to reference radius
4383 Int_t nseed = array->GetEntriesFast();
4384 if (nseed<=0) return;
4385 Float_t * quality = new Float_t[nseed];
4386 Int_t * indexes = new Int_t[nseed];
4387 for (Int_t i=0; i<nseed; i++) {
4388 AliTPCseed *pt=(AliTPCseed*)array->UncheckedAt(i);
4393 pt->UpdatePoints(); //select first last max dens points
4394 Float_t * points = pt->GetPoints();
4395 if (points[3]<0.8) quality[i] =-1;
4396 quality[i] = (points[2]-points[0])+pt->GetNumberOfClusters();
4397 //prefer high momenta tracks if overlaps
4398 quality[i] *= TMath::Sqrt(TMath::Abs(pt->Pt())+0.5);
4400 AliTracker::PropagateTrackToBxByBz(&(params[i]),xref,pt->GetMass(),5.,kTRUE);
4401 AliTracker::PropagateTrackToBxByBz(&(params[i]),xref,pt->GetMass(),1.,kTRUE);
4403 TMath::Sort(nseed,quality,indexes);
4405 // 3. Loop over pair of tracks
4407 for (Int_t i0=0; i0<nseed; i0++) {
4408 Int_t index0=indexes[i0];
4409 if (!(array->UncheckedAt(index0))) continue;
4410 AliTPCseed *s1 = (AliTPCseed*)array->UncheckedAt(index0);
4411 if (!s1->IsActive()) continue;
4412 AliExternalTrackParam &par0=params[index0];
4413 for (Int_t i1=i0+1; i1<nseed; i1++) {
4414 Int_t index1=indexes[i1];
4415 if (!(array->UncheckedAt(index1))) continue;
4416 AliTPCseed *s2 = (AliTPCseed*)array->UncheckedAt(index1);
4417 if (!s2->IsActive()) continue;
4418 if (s2->GetKinkIndexes()[0]!=0)
4419 if (s2->GetKinkIndexes()[0] == -s1->GetKinkIndexes()[0]) continue;
4420 AliExternalTrackParam &par1=params[index1];
4421 if (TMath::Abs(par0.GetParameter()[3]-par1.GetParameter()[3])>kCutP3) continue;
4422 if (TMath::Abs(par0.GetParameter()[1]-par1.GetParameter()[1])>kCutP1) continue;
4423 if (TMath::Abs(par0.GetParameter()[2]-par1.GetParameter()[2])>kCutP2) continue;
4424 Double_t dAlpha= TMath::Abs(par0.GetAlpha()-par1.GetAlpha());
4425 if (dAlpha>TMath::Pi()) dAlpha-=TMath::Pi();
4426 if (TMath::Abs(dAlpha)>kCutAlpha) continue;
4431 Int_t firstShared=lastpoint, lastShared=firstpoint;
4432 Int_t firstRow=lastpoint, lastRow=firstpoint;
4434 for (Int_t i=firstpoint;i<lastpoint;i++){
4435 if (s1->GetClusterIndex2(i)>0) nall0++;
4436 if (s2->GetClusterIndex2(i)>0) nall1++;
4437 if (s1->GetClusterIndex2(i)>0 && s2->GetClusterIndex2(i)>0) {
4438 if (i<firstRow) firstRow=i;
4439 if (i>lastRow) lastRow=i;
4441 if ( (s1->GetClusterIndex2(i))==(s2->GetClusterIndex2(i)) && s1->GetClusterIndex2(i)>0) {
4442 if (i<firstShared) firstShared=i;
4443 if (i>lastShared) lastShared=i;
4447 Double_t ratio0 = Float_t(sumShared)/Float_t(TMath::Min(nall0+1,nall1+1));
4448 Double_t ratio1 = Float_t(sumShared)/Float_t(TMath::Max(nall0+1,nall1+1));
4450 if( AliTPCReconstructor::StreamLevel()>1){
4451 TTreeSRedirector &cstream = *fDebugStreamer;
4452 Int_t n0=s1->GetNumberOfClusters();
4453 Int_t n1=s2->GetNumberOfClusters();
4454 Int_t n0F=s1->GetNFoundable();
4455 Int_t n1F=s2->GetNFoundable();
4456 Int_t lab0=s1->GetLabel();
4457 Int_t lab1=s2->GetLabel();
4459 cstream<<"Splitted2"<<
4460 "iter="<<fIteration<<
4461 "lab0="<<lab0<< // MC label if exist
4462 "lab1="<<lab1<< // MC label if exist
4465 "ratio0="<<ratio0<< // shared ratio
4466 "ratio1="<<ratio1<< // shared ratio
4467 "p0.="<<&par0<< // track parameters
4469 "s0.="<<s1<< // full seed
4471 "n0="<<n0<< // number of clusters track 0
4472 "n1="<<n1<< // number of clusters track 1
4473 "nall0="<<nall0<< // number of clusters track 0
4474 "nall1="<<nall1<< // number of clusters track 1
4475 "n0F="<<n0F<< // number of findable
4476 "n1F="<<n1F<< // number of findable
4477 "shared="<<sumShared<< // number of shared clusters
4478 "firstS="<<firstShared<< // first and the last shared row
4479 "lastS="<<lastShared<<
4480 "firstRow="<<firstRow<< // first and the last row with cluster
4481 "lastRow="<<lastRow<< //
4485 // remove track with lower quality
4487 if (ratio0>AliTPCReconstructor::GetRecoParam()->GetCutSharedClusters(0) ||
4488 ratio1>AliTPCReconstructor::GetRecoParam()->GetCutSharedClusters(1)){
4492 delete array->RemoveAt(index1);
4497 // 4. Delete temporary array
4507 void AliTPCtrackerMI::FindCurling(const TObjArray * array, AliESDEvent */*esd*/, Int_t iter)
4510 // find Curling tracks
4511 // Use AliTPCReconstructor::StreamLevel()>1 if you want to tune parameters - cuts
4514 // Algorithm done in 2 phases - because of CPU consumption
4515 // it is n^2 algorithm - for lead-lead 20000x20000 combination are investigated
4516 // see detal in MC part what can be used to cut
4520 const Float_t kMaxC = 400; // maximal curvature to of the track
4521 const Float_t kMaxdTheta = 0.15; // maximal distance in theta
4522 const Float_t kMaxdPhi = 0.15; // maximal distance in phi
4523 const Float_t kPtRatio = 0.3; // ratio between pt
4524 const Float_t kMinDCAR = 2.; // distance to the primary vertex in r - see cpipe cut
4527 // Curling tracks cuts
4530 const Float_t kMaxDeltaRMax = 40; // distance in outer radius
4531 const Float_t kMaxDeltaRMin = 5.; // distance in lower radius - see cpipe cut
4532 const Float_t kMinAngle = 2.9; // angle between tracks
4533 const Float_t kMaxDist = 5; // biggest distance
4535 // The cuts can be tuned using the "MC information stored in Multi tree ==> see FindMultiMC
4538 TCut csign("csign","Tr0.fP[4]*Tr1.fP[4]<0"); //opposite sign
4539 TCut cmax("cmax","abs(Tr0.GetC())>1/400");
4540 TCut cda("cda","sqrt(dtheta^2+dfi^2)<0.15");
4541 TCut ccratio("ccratio","abs((Tr0.fP[4]+Tr1.fP[4])/(abs(Tr0.fP[4])+abs(Tr1.fP[4])))<0.3");
4542 TCut cpipe("cpipe", "min(abs(r0-rc0),abs(r1-rc1))>5");
4544 TCut cdrmax("cdrmax","abs(abs(rc0+r0)-abs(rc1+r1))<40")
4545 TCut cdrmin("cdrmin","abs(abs(rc0+r0)-abs(rc1+r1))<10")
4547 Multi->Draw("dfi","iter==0"+csign+cmax+cda+ccratio); ~94% of curling tracks fulfill
4548 Multi->Draw("min(abs(r0-rc0),abs(r1-rc1))","iter==0&&abs(lab1)==abs(lab0)"+csign+cmax+cda+ccratio+cpipe+cdrmin+cdrmax); //80%
4550 Curling2->Draw("dfi","iter==0&&abs(lab0)==abs(lab1)"+csign+cmax+cdtheta+cdfi+ccratio)
4556 Int_t nentries = array->GetEntriesFast();
4557 AliHelix *helixes = new AliHelix[nentries];
4558 for (Int_t i=0;i<nentries;i++){
4559 AliTPCseed* track = (AliTPCseed*)array->At(i);
4560 if (!track) continue;
4561 track->SetCircular(0);
4562 new (&helixes[i]) AliHelix(*track);
4568 Double_t phase[2][2],radius[2];
4573 for (Int_t i0=0;i0<nentries;i0++){
4574 AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
4575 if (!track0) continue;
4576 if (TMath::Abs(track0->GetC())<1/kMaxC) continue;
4577 Float_t xc0 = helixes[i0].GetHelix(6);
4578 Float_t yc0 = helixes[i0].GetHelix(7);
4579 Float_t r0 = helixes[i0].GetHelix(8);
4580 Float_t rc0 = TMath::Sqrt(xc0*xc0+yc0*yc0);
4581 Float_t fi0 = TMath::ATan2(yc0,xc0);
4583 for (Int_t i1=i0+1;i1<nentries;i1++){
4584 AliTPCseed * track1 = (AliTPCseed*)array->At(i1);
4585 if (!track1) continue;
4586 if (TMath::Abs(track1->GetC())<1/kMaxC) continue;
4587 Float_t xc1 = helixes[i1].GetHelix(6);
4588 Float_t yc1 = helixes[i1].GetHelix(7);
4589 Float_t r1 = helixes[i1].GetHelix(8);
4590 Float_t rc1 = TMath::Sqrt(xc1*xc1+yc1*yc1);
4591 Float_t fi1 = TMath::ATan2(yc1,xc1);
4593 Float_t dfi = fi0-fi1;
4596 if (dfi>1.5*TMath::Pi()) dfi-=TMath::Pi(); // take care about edge effect
4597 if (dfi<-1.5*TMath::Pi()) dfi+=TMath::Pi(); //
4598 Float_t dtheta = TMath::Abs(track0->GetTgl()-track1->GetTgl())<TMath::Abs(track0->GetTgl()+track1->GetTgl())? track0->GetTgl()-track1->GetTgl():track0->GetTgl()+track1->GetTgl();
4602 if (track0->GetBConstrain()&&track1->GetBConstrain()) continue; // not constrained
4603 if (track1->GetSigned1Pt()*track0->GetSigned1Pt()>0) continue; // not the same sign
4604 if ( TMath::Abs(track1->GetTgl()+track0->GetTgl())>kMaxdTheta) continue; //distance in the Theta
4605 if ( TMath::Abs(dfi)>kMaxdPhi) continue; //distance in phi
4606 if ( TMath::Sqrt(dfi*dfi+dtheta*dtheta)>kMaxdPhi) continue; //common angular offset
4608 Float_t pt0 = track0->GetSignedPt();
4609 Float_t pt1 = track1->GetSignedPt();
4610 if ((TMath::Abs(pt0+pt1)/(TMath::Abs(pt0)+TMath::Abs(pt1)))>kPtRatio) continue;
4611 if ((iter==1) && TMath::Abs(TMath::Abs(rc0+r0)-TMath::Abs(rc1+r1))>kMaxDeltaRMax) continue;
4612 if ((iter!=1) &&TMath::Abs(TMath::Abs(rc0-r0)-TMath::Abs(rc1-r1))>kMaxDeltaRMin) continue;
4613 if (TMath::Min(TMath::Abs(rc0-r0),TMath::Abs(rc1-r1))<kMinDCAR) continue;
4616 // Now find closest approach
4620 Int_t npoints = helixes[i0].GetRPHIintersections(helixes[i1], phase, radius,10);
4621 if (npoints==0) continue;
4622 helixes[i0].GetClosestPhases(helixes[i1], phase);
4626 Double_t hangles[3];
4627 helixes[i0].Evaluate(phase[0][0],xyz0);
4628 helixes[i1].Evaluate(phase[0][1],xyz1);
4630 helixes[i0].GetAngle(phase[0][0],helixes[i1],phase[0][1],hangles);
4631 Double_t deltah[2],deltabest;
4632 if (TMath::Abs(hangles[2])<kMinAngle) continue;
4636 helixes[i0].ParabolicDCA(helixes[i1],phase[0][0],phase[0][1],radius[0],deltah[0],2);
4638 helixes[i0].ParabolicDCA(helixes[i1],phase[1][0],phase[1][1],radius[1],deltah[1],2);
4639 if (deltah[1]<deltah[0]) ibest=1;
4641 deltabest = TMath::Sqrt(deltah[ibest]);
4642 helixes[i0].Evaluate(phase[ibest][0],xyz0);
4643 helixes[i1].Evaluate(phase[ibest][1],xyz1);
4644 helixes[i0].GetAngle(phase[ibest][0],helixes[i1],phase[ibest][1],hangles);
4645 Double_t radiusbest = TMath::Sqrt(radius[ibest]);
4647 if (deltabest>kMaxDist) continue;
4648 // if (mindcar+mindcaz<40 && (TMath::Abs(hangles[2])<kMinAngle ||deltabest>3)) continue;
4649 Bool_t sign =kFALSE;
4650 if (hangles[2]>kMinAngle) sign =kTRUE;
4653 // circular[i0] = kTRUE;
4654 // circular[i1] = kTRUE;
4655 if (track0->OneOverPt()<track1->OneOverPt()){
4656 track0->SetCircular(track0->GetCircular()+1);
4657 track1->SetCircular(track1->GetCircular()+2);
4660 track1->SetCircular(track1->GetCircular()+1);
4661 track0->SetCircular(track0->GetCircular()+2);
4664 if (AliTPCReconstructor::StreamLevel()>2){
4666 //debug stream to tune "fine" cuts
4667 Int_t lab0=track0->GetLabel();
4668 Int_t lab1=track1->GetLabel();
4669 TTreeSRedirector &cstream = *fDebugStreamer;
4670 cstream<<"Curling2"<<
4686 "npoints="<<npoints<<
4687 "hangles0="<<hangles[0]<<
4688 "hangles1="<<hangles[1]<<
4689 "hangles2="<<hangles[2]<<
4692 "radius="<<radiusbest<<
4693 "deltabest="<<deltabest<<
4694 "phase0="<<phase[ibest][0]<<
4695 "phase1="<<phase[ibest][1]<<
4703 if (AliTPCReconstructor::StreamLevel()>1) {
4704 AliInfo("Time for curling tracks removal");
4713 void AliTPCtrackerMI::FindKinks(TObjArray * array, AliESDEvent *esd)
4720 TObjArray *kinks= new TObjArray(10000);
4721 // TObjArray *v0s= new TObjArray(10000);
4722 Int_t nentries = array->GetEntriesFast();
4723 AliHelix *helixes = new AliHelix[nentries];
4724 Int_t *sign = new Int_t[nentries];
4725 Int_t *nclusters = new Int_t[nentries];
4726 Float_t *alpha = new Float_t[nentries];
4727 AliKink *kink = new AliKink();
4728 Int_t * usage = new Int_t[nentries];
4729 Float_t *zm = new Float_t[nentries];
4730 Float_t *z0 = new Float_t[nentries];
4731 Float_t *fim = new Float_t[nentries];
4732 Float_t *shared = new Float_t[nentries];
4733 Bool_t *circular = new Bool_t[nentries];
4734 Float_t *dca = new Float_t[nentries];
4735 //const AliESDVertex * primvertex = esd->GetVertex();
4737 // nentries = array->GetEntriesFast();
4742 for (Int_t i=0;i<nentries;i++){
4745 AliTPCseed* track = (AliTPCseed*)array->At(i);
4746 if (!track) continue;
4747 track->SetCircular(0);
4749 track->UpdatePoints();
4750 if (( track->GetPoints()[2]- track->GetPoints()[0])>5 && track->GetPoints()[3]>0.8){
4752 nclusters[i]=track->GetNumberOfClusters();
4753 alpha[i] = track->GetAlpha();
4754 new (&helixes[i]) AliHelix(*track);
4756 helixes[i].Evaluate(0,xyz);
4757 sign[i] = (track->GetC()>0) ? -1:1;
4760 if (track->GetProlongation(x,y,z)){
4762 fim[i] = alpha[i]+TMath::ATan2(y,x);
4765 zm[i] = track->GetZ();
4769 circular[i]= kFALSE;
4770 if (track->GetProlongation(0,y,z)) z0[i] = z;
4771 dca[i] = track->GetD(0,0);
4777 Int_t ncandidates =0;
4780 Double_t phase[2][2],radius[2];
4783 // Find circling track
4785 for (Int_t i0=0;i0<nentries;i0++){
4786 AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
4787 if (!track0) continue;
4788 if (track0->GetNumberOfClusters()<40) continue;
4789 if (TMath::Abs(1./track0->GetC())>200) continue;
4790 for (Int_t i1=i0+1;i1<nentries;i1++){
4791 AliTPCseed * track1 = (AliTPCseed*)array->At(i1);
4792 if (!track1) continue;
4793 if (track1->GetNumberOfClusters()<40) continue;
4794 if ( TMath::Abs(track1->GetTgl()+track0->GetTgl())>0.1) continue;
4795 if (track0->GetBConstrain()&&track1->GetBConstrain()) continue;
4796 if (TMath::Abs(1./track1->GetC())>200) continue;
4797 if (track1->GetSigned1Pt()*track0->GetSigned1Pt()>0) continue;
4798 if (track1->GetTgl()*track0->GetTgl()>0) continue;
4799 if (TMath::Max(TMath::Abs(1./track0->GetC()),TMath::Abs(1./track1->GetC()))>190) continue;
4800 if (track0->GetBConstrain()&&track1->OneOverPt()<track0->OneOverPt()) continue; //returning - lower momenta
4801 if (track1->GetBConstrain()&&track0->OneOverPt()<track1->OneOverPt()) continue; //returning - lower momenta
4803 Float_t mindcar = TMath::Min(TMath::Abs(dca[i0]),TMath::Abs(dca[i1]));
4804 if (mindcar<5) continue;
4805 Float_t mindcaz = TMath::Min(TMath::Abs(z0[i0]-GetZ()),TMath::Abs(z0[i1]-GetZ()));
4806 if (mindcaz<5) continue;
4807 if (mindcar+mindcaz<20) continue;
4810 Float_t xc0 = helixes[i0].GetHelix(6);
4811 Float_t yc0 = helixes[i0].GetHelix(7);
4812 Float_t r0 = helixes[i0].GetHelix(8);
4813 Float_t xc1 = helixes[i1].GetHelix(6);
4814 Float_t yc1 = helixes[i1].GetHelix(7);
4815 Float_t r1 = helixes[i1].GetHelix(8);
4817 Float_t rmean = (r0+r1)*0.5;
4818 Float_t delta =TMath::Sqrt((xc1-xc0)*(xc1-xc0)+(yc1-yc0)*(yc1-yc0));
4819 //if (delta>30) continue;
4820 if (delta>rmean*0.25) continue;
4821 if (TMath::Abs(r0-r1)/rmean>0.3) continue;
4823 Int_t npoints = helixes[i0].GetRPHIintersections(helixes[i1], phase, radius,10);
4824 if (npoints==0) continue;
4825 helixes[i0].GetClosestPhases(helixes[i1], phase);
4829 Double_t hangles[3];
4830 helixes[i0].Evaluate(phase[0][0],xyz0);
4831 helixes[i1].Evaluate(phase[0][1],xyz1);
4833 helixes[i0].GetAngle(phase[0][0],helixes[i1],phase[0][1],hangles);
4834 Double_t deltah[2],deltabest;
4835 if (hangles[2]<2.8) continue;
4838 helixes[i0].ParabolicDCA(helixes[i1],phase[0][0],phase[0][1],radius[0],deltah[0],2);
4840 helixes[i0].ParabolicDCA(helixes[i1],phase[1][0],phase[1][1],radius[1],deltah[1],2);
4841 if (deltah[1]<deltah[0]) ibest=1;
4843 deltabest = TMath::Sqrt(deltah[ibest]);
4844 helixes[i0].Evaluate(phase[ibest][0],xyz0);
4845 helixes[i1].Evaluate(phase[ibest][1],xyz1);
4846 helixes[i0].GetAngle(phase[ibest][0],helixes[i1],phase[ibest][1],hangles);
4847 Double_t radiusbest = TMath::Sqrt(radius[ibest]);
4849 if (deltabest>6) continue;
4850 if (mindcar+mindcaz<40 && (hangles[2]<3.12||deltabest>3)) continue;
4851 Bool_t lsign =kFALSE;
4852 if (hangles[2]>3.06) lsign =kTRUE;
4855 circular[i0] = kTRUE;
4856 circular[i1] = kTRUE;
4857 if (track0->OneOverPt()<track1->OneOverPt()){
4858 track0->SetCircular(track0->GetCircular()+1);
4859 track1->SetCircular(track1->GetCircular()+2);
4862 track1->SetCircular(track1->GetCircular()+1);
4863 track0->SetCircular(track0->GetCircular()+2);
4866 if (lsign&&AliTPCReconstructor::StreamLevel()>1){
4868 Int_t lab0=track0->GetLabel();
4869 Int_t lab1=track1->GetLabel();
4870 TTreeSRedirector &cstream = *fDebugStreamer;
4871 cstream<<"Curling"<<
4878 "mindcar="<<mindcar<<
4879 "mindcaz="<<mindcaz<<
4882 "npoints="<<npoints<<
4883 "hangles0="<<hangles[0]<<
4884 "hangles2="<<hangles[2]<<
4889 "radius="<<radiusbest<<
4890 "deltabest="<<deltabest<<
4891 "phase0="<<phase[ibest][0]<<
4892 "phase1="<<phase[ibest][1]<<
4902 for (Int_t i =0;i<nentries;i++){
4903 if (sign[i]==0) continue;
4904 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
4911 Double_t cradius0 = 40*40;
4912 Double_t cradius1 = 270*270;
4915 Double_t cdist3=0.55;
4916 for (Int_t j =i+1;j<nentries;j++){
4918 if (sign[j]*sign[i]<1) continue;
4919 if ( (nclusters[i]+nclusters[j])>200) continue;
4920 if ( (nclusters[i]+nclusters[j])<80) continue;
4921 if ( TMath::Abs(zm[i]-zm[j])>60.) continue;
4922 if ( TMath::Abs(fim[i]-fim[j])>0.6 && TMath::Abs(fim[i]-fim[j])<5.7 ) continue;
4923 //AliTPCseed * track1 = (AliTPCseed*)array->At(j); Double_t phase[2][2],radius[2];
4924 Int_t npoints = helixes[i].GetRPHIintersections(helixes[j], phase, radius,20);
4925 if (npoints<1) continue;
4928 if (radius[0]<cradius0||radius[0]>cradius1) continue;
4931 if ( (radius[0]<cradius0||radius[0]>cradius1) && (radius[1]<cradius0||radius[1]>cradius1) ) continue;
4934 Double_t delta1=10000,delta2=10000;
4935 // cuts on the intersection radius
4936 helixes[i].LinearDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta1);
4937 if (radius[0]<20&&delta1<1) continue; //intersection at vertex
4938 if (radius[0]<10&&delta1<3) continue; //intersection at vertex
4940 helixes[i].LinearDCA(helixes[j],phase[1][0],phase[1][1],radius[1],delta2);
4941 if (radius[1]<20&&delta2<1) continue; //intersection at vertex
4942 if (radius[1]<10&&delta2<3) continue; //intersection at vertex
4945 Double_t distance1 = TMath::Min(delta1,delta2);
4946 if (distance1>cdist1) continue; // cut on DCA linear approximation
4948 npoints = helixes[i].GetRPHIintersections(helixes[j], phase, radius,20);
4949 helixes[i].ParabolicDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta1);
4950 if (radius[0]<20&&delta1<1) continue; //intersection at vertex
4951 if (radius[0]<10&&delta1<3) continue; //intersection at vertex
4954 helixes[i].ParabolicDCA(helixes[j],phase[1][0],phase[1][1],radius[1],delta2);
4955 if (radius[1]<20&&delta2<1) continue; //intersection at vertex
4956 if (radius[1]<10&&delta2<3) continue; //intersection at vertex
4958 distance1 = TMath::Min(delta1,delta2);
4961 rkink = TMath::Sqrt(radius[0]);
4964 rkink = TMath::Sqrt(radius[1]);
4966 if (distance1>cdist2) continue;
4969 AliTPCseed * track1 = (AliTPCseed*)array->At(j);
4972 Int_t row0 = GetRowNumber(rkink);
4973 if (row0<10) continue;
4974 if (row0>150) continue;
4977 Float_t dens00=-1,dens01=-1;
4978 Float_t dens10=-1,dens11=-1;
4980 Int_t found,foundable,ishared;
4981 track0->GetClusterStatistic(0,row0-5, found, foundable,ishared,kFALSE);
4982 if (foundable>5) dens00 = Float_t(found)/Float_t(foundable);
4983 track0->GetClusterStatistic(row0+5,155, found, foundable,ishared,kFALSE);
4984 if (foundable>5) dens01 = Float_t(found)/Float_t(foundable);
4986 track1->GetClusterStatistic(0,row0-5, found, foundable,ishared,kFALSE);
4987 if (foundable>10) dens10 = Float_t(found)/Float_t(foundable);
4988 track1->GetClusterStatistic(row0+5,155, found, foundable,ishared,kFALSE);
4989 if (foundable>10) dens11 = Float_t(found)/Float_t(foundable);
4991 if (dens00<dens10 && dens01<dens11) continue;
4992 if (dens00>dens10 && dens01>dens11) continue;
4993 if (TMath::Max(dens00,dens10)<0.1) continue;
4994 if (TMath::Max(dens01,dens11)<0.3) continue;
4996 if (TMath::Min(dens00,dens10)>0.6) continue;
4997 if (TMath::Min(dens01,dens11)>0.6) continue;
5000 AliTPCseed * ktrack0, *ktrack1;
5009 if (TMath::Abs(ktrack0->GetC())>5) continue; // cut on the curvature for mother particle
5010 AliExternalTrackParam paramm(*ktrack0);
5011 AliExternalTrackParam paramd(*ktrack1);
5012 if (row0>60&&ktrack1->GetReference().GetX()>90.)new (¶md) AliExternalTrackParam(ktrack1->GetReference());
5015 kink->SetMother(paramm);
5016 kink->SetDaughter(paramd);
5019 Float_t x[3] = { kink->GetPosition()[0],kink->GetPosition()[1],kink->GetPosition()[2]};
5021 fkParam->Transform0to1(x,index);
5022 fkParam->Transform1to2(x,index);
5023 row0 = GetRowNumber(x[0]);
5025 if (kink->GetR()<100) continue;
5026 if (kink->GetR()>240) continue;
5027 if (kink->GetPosition()[2]/kink->GetR()>AliTPCReconstructor::GetCtgRange()) continue; //out of fiducial volume
5028 if (kink->GetDistance()>cdist3) continue;
5029 Float_t dird = kink->GetDaughterP()[0]*kink->GetPosition()[0]+kink->GetDaughterP()[1]*kink->GetPosition()[1]; // rough direction estimate
5030 if (dird<0) continue;
5032 Float_t dirm = kink->GetMotherP()[0]*kink->GetPosition()[0]+kink->GetMotherP()[1]*kink->GetPosition()[1]; // rough direction estimate
5033 if (dirm<0) continue;
5034 Float_t mpt = TMath::Sqrt(kink->GetMotherP()[0]*kink->GetMotherP()[0]+kink->GetMotherP()[1]*kink->GetMotherP()[1]);
5035 if (mpt<0.2) continue;
5038 //for high momenta momentum not defined well in first iteration
5039 Double_t qt = TMath::Sin(kink->GetAngle(2))*ktrack1->GetP();
5040 if (qt>0.35) continue;
5043 kink->SetLabel(CookLabel(ktrack0,0.4,0,row0),0);
5044 kink->SetLabel(CookLabel(ktrack1,0.4,row0,160),1);
5046 kink->SetTPCDensity(dens00,0,0);
5047 kink->SetTPCDensity(dens01,0,1);
5048 kink->SetTPCDensity(dens10,1,0);
5049 kink->SetTPCDensity(dens11,1,1);
5050 kink->SetIndex(i,0);
5051 kink->SetIndex(j,1);
5054 kink->SetTPCDensity(dens10,0,0);
5055 kink->SetTPCDensity(dens11,0,1);
5056 kink->SetTPCDensity(dens00,1,0);
5057 kink->SetTPCDensity(dens01,1,1);
5058 kink->SetIndex(j,0);
5059 kink->SetIndex(i,1);
5062 if (mpt<1||kink->GetAngle(2)>0.1){
5063 // angle and densities not defined yet
5064 if (kink->GetTPCDensityFactor()<0.8) continue;
5065 if ((2-kink->GetTPCDensityFactor())*kink->GetDistance() >0.25) continue;
5066 if (kink->GetAngle(2)*ktrack0->GetP()<0.003) continue; //too small angle
5067 if (kink->GetAngle(2)>0.2&&kink->GetTPCDensityFactor()<1.15) continue;
5068 if (kink->GetAngle(2)>0.2&&kink->GetTPCDensity(0,1)>0.05) continue;
5070 Float_t criticalangle = track0->GetSigmaSnp2()+track0->GetSigmaTgl2();
5071 criticalangle+= track1->GetSigmaSnp2()+track1->GetSigmaTgl2();
5072 criticalangle= 3*TMath::Sqrt(criticalangle);
5073 if (criticalangle>0.02) criticalangle=0.02;
5074 if (kink->GetAngle(2)<criticalangle) continue;
5077 Int_t drow = Int_t(2.+0.5/(0.05+kink->GetAngle(2))); // overlap region defined
5078 Float_t shapesum =0;
5080 for ( Int_t row = row0-drow; row<row0+drow;row++){
5081 if (row<0) continue;
5082 if (row>155) continue;
5083 if (ktrack0->GetClusterPointer(row)){
5084 AliTPCTrackerPoint *point =ktrack0->GetTrackPoint(row);
5085 shapesum+=point->GetSigmaY()+point->GetSigmaZ();
5088 if (ktrack1->GetClusterPointer(row)){
5089 AliTPCTrackerPoint *point =ktrack1->GetTrackPoint(row);
5090 shapesum+=point->GetSigmaY()+point->GetSigmaZ();
5095 kink->SetShapeFactor(-1.);
5098 kink->SetShapeFactor(shapesum/sum);
5100 // esd->AddKink(kink);
5102 // kink->SetMother(paramm);
5103 //kink->SetDaughter(paramd);
5105 Double_t chi2P2 = paramm.GetParameter()[2]-paramd.GetParameter()[2];
5107 chi2P2/=paramm.GetCovariance()[5]+paramd.GetCovariance()[5];
5108 Double_t chi2P3 = paramm.GetParameter()[3]-paramd.GetParameter()[3];
5110 chi2P3/=paramm.GetCovariance()[9]+paramd.GetCovariance()[9];
5112 if (AliTPCReconstructor::StreamLevel()>1) {
5113 (*fDebugStreamer)<<"kinkLpt"<<
5121 if ( chi2P2+chi2P3<AliTPCReconstructor::GetRecoParam()->GetKinkAngleCutChi2(0)){
5125 kinks->AddLast(kink);
5131 // sort the kinks according quality - and refit them towards vertex
5133 Int_t nkinks = kinks->GetEntriesFast();
5134 Float_t *quality = new Float_t[nkinks];
5135 Int_t *indexes = new Int_t[nkinks];
5136 AliTPCseed *mothers = new AliTPCseed[nkinks];
5137 AliTPCseed *daughters = new AliTPCseed[nkinks];
5140 for (Int_t i=0;i<nkinks;i++){
5142 AliKink *kinkl = (AliKink*)kinks->At(i);
5144 // refit kinks towards vertex
5146 Int_t index0 = kinkl->GetIndex(0);
5147 Int_t index1 = kinkl->GetIndex(1);
5148 AliTPCseed * ktrack0 = (AliTPCseed*)array->At(index0);
5149 AliTPCseed * ktrack1 = (AliTPCseed*)array->At(index1);
5151 Int_t sumn=ktrack0->GetNumberOfClusters()+ktrack1->GetNumberOfClusters();
5153 // Refit Kink under if too small angle
5155 if (kinkl->GetAngle(2)<0.05){
5156 kinkl->SetTPCRow0(GetRowNumber(kinkl->GetR()));
5157 Int_t row0 = kinkl->GetTPCRow0();
5158 Int_t drow = Int_t(2.+0.5/(0.05+kinkl->GetAngle(2)));
5161 Int_t last = row0-drow;
5162 if (last<40) last=40;
5163 if (last<ktrack0->GetFirstPoint()+25) last = ktrack0->GetFirstPoint()+25;
5164 AliTPCseed* seed0 = ReSeed(ktrack0,last,kFALSE);
5167 Int_t first = row0+drow;
5168 if (first>130) first=130;
5169 if (first>ktrack1->GetLastPoint()-25) first = TMath::Max(ktrack1->GetLastPoint()-25,30);
5170 AliTPCseed* seed1 = ReSeed(ktrack1,first,kTRUE);
5172 if (seed0 && seed1){
5173 kinkl->SetStatus(1,8);
5174 if (RefitKink(*seed0,*seed1,*kinkl)) kinkl->SetStatus(1,9);
5175 row0 = GetRowNumber(kinkl->GetR());
5176 sumn = seed0->GetNumberOfClusters()+seed1->GetNumberOfClusters();
5177 mothers[i] = *seed0;
5178 daughters[i] = *seed1;
5181 delete kinks->RemoveAt(i);
5182 if (seed0) delete seed0;
5183 if (seed1) delete seed1;
5186 if (kinkl->GetDistance()>0.5 || kinkl->GetR()<110 || kinkl->GetR()>240) {
5187 delete kinks->RemoveAt(i);
5188 if (seed0) delete seed0;
5189 if (seed1) delete seed1;
5197 if (kinkl) quality[i] = 160*((0.1+kinkl->GetDistance())*(2.-kinkl->GetTPCDensityFactor()))/(sumn+40.); //the longest -clossest will win
5199 TMath::Sort(nkinks,quality,indexes,kFALSE);
5201 //remove double find kinks
5203 for (Int_t ikink0=1;ikink0<nkinks;ikink0++){
5204 AliKink * kink0 = (AliKink*) kinks->At(indexes[ikink0]);
5205 if (!kink0) continue;
5207 for (Int_t ikink1=0;ikink1<ikink0;ikink1++){
5208 if (!kink0) continue;
5209 AliKink * kink1 = (AliKink*) kinks->At(indexes[ikink1]);
5210 if (!kink1) continue;
5211 // if not close kink continue
5212 if (TMath::Abs(kink1->GetPosition()[2]-kink0->GetPosition()[2])>10) continue;
5213 if (TMath::Abs(kink1->GetPosition()[1]-kink0->GetPosition()[1])>10) continue;
5214 if (TMath::Abs(kink1->GetPosition()[0]-kink0->GetPosition()[0])>10) continue;
5216 AliTPCseed &mother0 = mothers[indexes[ikink0]];
5217 AliTPCseed &daughter0 = daughters[indexes[ikink0]];
5218 AliTPCseed &mother1 = mothers[indexes[ikink1]];
5219 AliTPCseed &daughter1 = daughters[indexes[ikink1]];
5220 Int_t row0 = (kink0->GetTPCRow0()+kink1->GetTPCRow0())/2;
5229 for (Int_t i=0;i<row0;i++){
5230 if (mother0.GetClusterIndex(i)>0 && mother1.GetClusterIndex(i)>0){
5233 if (mother0.GetClusterIndex(i)==mother1.GetClusterIndex(i)){
5240 for (Int_t i=row0;i<158;i++){
5241 if (daughter0.GetClusterIndex(i)>0 && daughter0.GetClusterIndex(i)>0){
5244 if (mother0.GetClusterIndex(i)==mother1.GetClusterIndex(i)){
5250 Float_t ratio = Float_t(same+1)/Float_t(both+1);
5251 Float_t ratiom = Float_t(samem+1)/Float_t(bothm+1);
5252 Float_t ratiod = Float_t(samed+1)/Float_t(bothd+1);
5253 if (ratio>0.3 && ratiom>0.5 &&ratiod>0.5) {
5254 Int_t sum0 = mother0.GetNumberOfClusters()+daughter0.GetNumberOfClusters();
5255 Int_t sum1 = mother1.GetNumberOfClusters()+daughter1.GetNumberOfClusters();
5257 shared[kink0->GetIndex(0)]= kTRUE;
5258 shared[kink0->GetIndex(1)]= kTRUE;
5259 delete kinks->RemoveAt(indexes[ikink0]);
5262 shared[kink1->GetIndex(0)]= kTRUE;
5263 shared[kink1->GetIndex(1)]= kTRUE;
5264 delete kinks->RemoveAt(indexes[ikink1]);
5271 for (Int_t i=0;i<nkinks;i++){
5272 AliKink * kinkl = (AliKink*) kinks->At(indexes[i]);
5273 if (!kinkl) continue;
5274 kinkl->SetTPCRow0(GetRowNumber(kinkl->GetR()));
5275 Int_t index0 = kinkl->GetIndex(0);
5276 Int_t index1 = kinkl->GetIndex(1);
5277 if (circular[index0]||(circular[index1]&&kinkl->GetDistance()>0.2)) continue;
5278 kinkl->SetMultiple(usage[index0],0);
5279 kinkl->SetMultiple(usage[index1],1);
5280 if (kinkl->GetMultiple()[0]+kinkl->GetMultiple()[1]>2) continue;
5281 if (kinkl->GetMultiple()[0]+kinkl->GetMultiple()[1]>0 && quality[indexes[i]]>0.2) continue;
5282 if (kinkl->GetMultiple()[0]+kinkl->GetMultiple()[1]>0 && kinkl->GetDistance()>0.2) continue;
5283 if (circular[index0]||(circular[index1]&&kinkl->GetDistance()>0.1)) continue;
5285 AliTPCseed * ktrack0 = (AliTPCseed*)array->At(index0);
5286 AliTPCseed * ktrack1 = (AliTPCseed*)array->At(index1);
5287 if (!ktrack0 || !ktrack1) continue;
5288 Int_t index = esd->AddKink(kinkl);
5291 if ( ktrack0->GetKinkIndex(0)==0 && ktrack1->GetKinkIndex(0)==0) { //best kink
5292 if (mothers[indexes[i]].GetNumberOfClusters()>20 && daughters[indexes[i]].GetNumberOfClusters()>20 && (mothers[indexes[i]].GetNumberOfClusters()+daughters[indexes[i]].GetNumberOfClusters())>100){
5293 *ktrack0 = mothers[indexes[i]];
5294 *ktrack1 = daughters[indexes[i]];
5298 ktrack0->SetKinkIndex(usage[index0],-(index+1));
5299 ktrack1->SetKinkIndex(usage[index1], (index+1));
5304 // Remove tracks corresponding to shared kink's
5306 for (Int_t i=0;i<nentries;i++){
5307 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5308 if (!track0) continue;
5309 if (track0->GetKinkIndex(0)!=0) continue;
5310 if (shared[i]) delete array->RemoveAt(i);
5315 RemoveUsed2(array,0.5,0.4,30);
5317 for (Int_t i=0;i<nentries;i++){
5318 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5319 if (!track0) continue;
5320 track0->CookdEdx(0.02,0.6);
5324 for (Int_t i=0;i<nentries;i++){
5325 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5326 if (!track0) continue;
5327 if (track0->Pt()<1.4) continue;
5328 //remove double high momenta tracks - overlapped with kink candidates
5331 for (Int_t icl=track0->GetFirstPoint();icl<track0->GetLastPoint(); icl++){
5332 if (track0->GetClusterPointer(icl)!=0){
5334 if (track0->GetClusterPointer(icl)->IsUsed(10)) ishared++;
5337 if (Float_t(ishared+1)/Float_t(all+1)>0.5) {
5338 delete array->RemoveAt(i);
5342 if (track0->GetKinkIndex(0)!=0) continue;
5343 if (track0->GetNumberOfClusters()<80) continue;
5345 AliTPCseed *pmother = new AliTPCseed();
5346 AliTPCseed *pdaughter = new AliTPCseed();
5347 AliKink *pkink = new AliKink;
5349 AliTPCseed & mother = *pmother;
5350 AliTPCseed & daughter = *pdaughter;
5351 AliKink & kinkl = *pkink;
5352 if (CheckKinkPoint(track0,mother,daughter, kinkl)){
5353 if (mother.GetNumberOfClusters()<30||daughter.GetNumberOfClusters()<20) {
5357 continue; //too short tracks
5359 if (mother.Pt()<1.4) {
5365 Int_t row0= kinkl.GetTPCRow0();
5366 if (kinkl.GetDistance()>0.5 || kinkl.GetR()<110. || kinkl.GetR()>240.) {
5373 Int_t index = esd->AddKink(&kinkl);
5374 mother.SetKinkIndex(0,-(index+1));
5375 daughter.SetKinkIndex(0,index+1);
5376 if (mother.GetNumberOfClusters()>50) {
5377 delete array->RemoveAt(i);
5378 array->AddAt(new AliTPCseed(mother),i);
5381 array->AddLast(new AliTPCseed(mother));
5383 array->AddLast(new AliTPCseed(daughter));
5384 for (Int_t icl=0;icl<row0;icl++) {
5385 if (mother.GetClusterPointer(icl)) mother.GetClusterPointer(icl)->Use(20);
5388 for (Int_t icl=row0;icl<158;icl++) {
5389 if (daughter.GetClusterPointer(icl)) daughter.GetClusterPointer(icl)->Use(20);
5398 delete [] daughters;
5420 AliInfo(Form("Ncandidates=\t%d\t%d\t%d\t%d\n",esd->GetNumberOfKinks(),ncandidates,ntracks,nall));
5424 void AliTPCtrackerMI::FindV0s(const TObjArray * array, AliESDEvent *const esd)
5430 TObjArray *tpcv0s = new TObjArray(100000);
5431 Int_t nentries = array->GetEntriesFast();
5432 AliHelix *helixes = new AliHelix[nentries];
5433 Int_t *sign = new Int_t[nentries];
5434 Float_t *alpha = new Float_t[nentries];
5435 Float_t *z0 = new Float_t[nentries];
5436 Float_t *dca = new Float_t[nentries];
5437 Float_t *sdcar = new Float_t[nentries];
5438 Float_t *cdcar = new Float_t[nentries];
5439 Float_t *pulldcar = new Float_t[nentries];
5440 Float_t *pulldcaz = new Float_t[nentries];
5441 Float_t *pulldca = new Float_t[nentries];
5442 Bool_t *isPrim = new Bool_t[nentries];
5443 const AliESDVertex * primvertex = esd->GetVertex();
5444 Double_t zvertex = primvertex->GetZv();
5446 // nentries = array->GetEntriesFast();
5448 for (Int_t i=0;i<nentries;i++){
5451 AliTPCseed* track = (AliTPCseed*)array->At(i);
5452 if (!track) continue;
5453 track->GetV0Indexes()[0] = 0; //rest v0 indexes
5454 track->GetV0Indexes()[1] = 0; //rest v0 indexes
5455 track->GetV0Indexes()[2] = 0; //rest v0 indexes
5457 alpha[i] = track->GetAlpha();
5458 new (&helixes[i]) AliHelix(*track);
5460 helixes[i].Evaluate(0,xyz);
5461 sign[i] = (track->GetC()>0) ? -1:1;
5465 if (track->GetProlongation(0,y,z)) z0[i] = z;
5466 dca[i] = track->GetD(0,0);
5468 // dca error parrameterezation + pulls
5470 sdcar[i] = TMath::Sqrt(0.150*0.150+(100*track->GetC())*(100*track->GetC()));
5471 if (TMath::Abs(track->GetTgl())>1) sdcar[i]*=2.5;
5472 cdcar[i] = TMath::Exp((TMath::Abs(track->GetC())-0.0106)*525.3);
5473 pulldcar[i] = (dca[i]-cdcar[i])/sdcar[i];
5474 pulldcaz[i] = (z0[i]-zvertex)/sdcar[i];
5475 pulldca[i] = TMath::Sqrt(pulldcar[i]*pulldcar[i]+pulldcaz[i]*pulldcaz[i]);
5476 if (track->TPCrPID(1)+track->TPCrPID(2)+track->TPCrPID(3)>0.5) {
5477 if (pulldca[i]<3.) isPrim[i]=kTRUE; //pion, muon and Kaon 3 sigma cut
5479 if (track->TPCrPID(4)>0.5) {
5480 if (pulldca[i]<0.5) isPrim[i]=kTRUE; //proton 0.5 sigma cut
5482 if (track->TPCrPID(0)>0.4) {
5483 isPrim[i]=kFALSE; //electron no sigma cut
5490 Int_t ncandidates =0;
5493 Double_t phase[2][2],radius[2];
5499 Float_t fprimvertex[3]={GetX(),GetY(),GetZ()};
5501 Double_t cradius0 = 10*10;
5502 Double_t cradius1 = 200*200;
5505 Double_t cpointAngle = 0.95;
5507 Double_t delta[2]={10000,10000};
5508 for (Int_t i =0;i<nentries;i++){
5509 if (sign[i]==0) continue;
5510 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5511 if (!track0) continue;
5512 if (AliTPCReconstructor::StreamLevel()>1){
5513 TTreeSRedirector &cstream = *fDebugStreamer;
5518 "zvertex="<<zvertex<<
5519 "sdcar0="<<sdcar[i]<<
5520 "cdcar0="<<cdcar[i]<<
5521 "pulldcar0="<<pulldcar[i]<<
5522 "pulldcaz0="<<pulldcaz[i]<<
5523 "pulldca0="<<pulldca[i]<<
5524 "isPrim="<<isPrim[i]<<
5528 if (track0->GetSigned1Pt()<0) continue;
5529 if (track0->GetKinkIndex(0)>0||isPrim[i]) continue; //daughter kink
5531 if (TMath::Abs(helixes[i].GetHelix(4))<0.000000001) continue;
5536 for (Int_t j =0;j<nentries;j++){
5537 AliTPCseed * track1 = (AliTPCseed*)array->At(j);
5538 if (!track1) continue;
5539 if (track1->GetKinkIndex(0)>0 || isPrim[j]) continue; //daughter kink
5540 if (sign[j]*sign[i]>0) continue;
5541 if (TMath::Abs(helixes[j].GetHelix(4))<0.000001) continue;
5542 if (track0->GetCircular()+track1->GetCircular()>1) continue; //circling -returning track
5545 // DCA to prim vertex cut
5551 Int_t npoints = helixes[i].GetRPHIintersections(helixes[j], phase, radius,cdist2);
5552 if (npoints<1) continue;
5556 if (radius[0]<cradius0||radius[0]>cradius1) continue;
5557 helixes[i].LinearDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta[0]);
5558 if (delta[0]>cdist1) continue;
5561 if (TMath::Max(radius[0],radius[1])<cradius0|| TMath::Min(radius[0],radius[1])>cradius1) continue;
5562 helixes[i].LinearDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta[0]);
5563 helixes[i].LinearDCA(helixes[j],phase[1][0],phase[1][1],radius[1],delta[1]);
5564 if (delta[1]<delta[0]) iclosest=1;
5565 if (delta[iclosest]>cdist1) continue;
5567 helixes[i].ParabolicDCA(helixes[j],phase[iclosest][0],phase[iclosest][1],radius[iclosest],delta[iclosest]);
5568 if (radius[iclosest]<cradius0 || radius[iclosest]>cradius1 || delta[iclosest]>cdist1) continue;
5570 Double_t pointAngle = helixes[i].GetPointAngle(helixes[j],phase[iclosest],fprimvertex);
5571 if (pointAngle<cpointAngle) continue;
5573 Bool_t isGamma = kFALSE;
5574 vertex.SetParamP(*track0); //track0 - plus
5575 vertex.SetParamN(*track1); //track1 - minus
5576 vertex.Update(fprimvertex);
5577 if (track0->TPCrPID(0)>0.3&&track1->TPCrPID(0)>0.3&&vertex.GetAnglep()[2]<0.15) isGamma=kTRUE; // gamma conversion candidate
5578 Double_t pointAngle2 = vertex.GetV0CosineOfPointingAngle();
5580 if (vertex.GetV0CosineOfPointingAngle()<cpointAngle && (!isGamma)) continue;// point angle cut
5581 //Bo: if (vertex.GetDist2()>2&&(!isGamma)) continue; // point angle cut
5582 if (vertex.GetDcaV0Daughters()>2&&(!isGamma)) continue;//Bo: // point angle cut
5583 Float_t sigmae = 0.15*0.15;
5584 if (vertex.GetRr()<80)
5585 sigmae += (sdcar[i]*sdcar[i]+sdcar[j]*sdcar[j])*(1.-vertex.GetRr()/80.)*(1.-vertex.GetRr()/80.);
5586 sigmae+= TMath::Sqrt(sigmae);
5587 //Bo: if (vertex.GetDist2()/sigmae>3.&&(!isGamma)) continue;
5588 if (vertex.GetDcaV0Daughters()/sigmae>3.&&(!isGamma)) continue;
5589 Float_t densb0=0,densb1=0,densa0=0,densa1=0;
5590 Int_t row0 = GetRowNumber(vertex.GetRr());
5592 //Bo: if (vertex.GetDist2()>0.2) continue;
5593 if (vertex.GetDcaV0Daughters()>0.2) continue;
5594 densb0 = track0->Density2(0,row0-5);
5595 densb1 = track1->Density2(0,row0-5);
5596 if (densb0>0.3|| densb1>0.3) continue; //clusters before vertex
5597 densa0 = track0->Density2(row0+5,row0+40);
5598 densa1 = track1->Density2(row0+5,row0+40);
5599 if ((densa0<0.4|| densa1<0.4)&&(!isGamma)) continue; //missing clusters after vertex
5602 densa0 = track0->Density2(0,40); //cluster density
5603 densa1 = track1->Density2(0,40); //cluster density
5604 if ((vertex.GetRr()<80&&densa0+densa1<1.)&&(!isGamma)) continue;
5606 //Bo: vertex.SetLab(0,track0->GetLabel());
5607 //Bo: vertex.SetLab(1,track1->GetLabel());
5608 vertex.SetChi2After((densa0+densa1)*0.5);
5609 vertex.SetChi2Before((densb0+densb1)*0.5);
5610 vertex.SetIndex(0,i);
5611 vertex.SetIndex(1,j);
5612 //Bo: vertex.SetStatus(1); // TPC v0 candidate
5613 vertex.SetOnFlyStatus(2);//Bo: // TPC v0 candidate
5614 //Bo: vertex.SetRp(track0->TPCrPIDs());
5615 //Bo: vertex.SetRm(track1->TPCrPIDs());
5616 tpcv0s->AddLast(new AliESDv0(vertex));
5619 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
5620 Double_t radiusm= (delta[0]<delta[1])? TMath::Sqrt(radius[0]):TMath::Sqrt(radius[1]);
5621 Double_t deltam= (delta[0]<delta[1])? TMath::Sqrt(delta[0]):TMath::Sqrt(delta[1]);
5622 if (AliTPCReconstructor::StreamLevel()>1) {
5623 Int_t lab0=track0->GetLabel();
5624 Int_t lab1=track1->GetLabel();
5625 Char_t c0=track0->GetCircular();
5626 Char_t c1=track1->GetCircular();
5627 TTreeSRedirector &cstream = *fDebugStreamer;
5630 "vertex.="<<&vertex<<
5633 "Helix0.="<<&helixes[i]<<
5636 "Helix1.="<<&helixes[j]<<
5637 "pointAngle="<<pointAngle<<
5638 "pointAngle2="<<pointAngle2<<
5643 "zvertex="<<zvertex<<
5646 "npoints="<<npoints<<
5647 "radius0="<<radius[0]<<
5648 "delta0="<<delta[0]<<
5649 "radius1="<<radius[1]<<
5650 "delta1="<<delta[1]<<
5651 "radiusm="<<radiusm<<
5653 "sdcar0="<<sdcar[i]<<
5654 "sdcar1="<<sdcar[j]<<
5655 "cdcar0="<<cdcar[i]<<
5656 "cdcar1="<<cdcar[j]<<
5657 "pulldcar0="<<pulldcar[i]<<
5658 "pulldcar1="<<pulldcar[j]<<
5659 "pulldcaz0="<<pulldcaz[i]<<
5660 "pulldcaz1="<<pulldcaz[j]<<
5661 "pulldca0="<<pulldca[i]<<
5662 "pulldca1="<<pulldca[j]<<
5672 Float_t *quality = new Float_t[ncandidates];
5673 Int_t *indexes = new Int_t[ncandidates];
5675 for (Int_t i=0;i<ncandidates;i++){
5677 AliESDv0 *v0 = (AliESDv0*)tpcv0s->At(i);
5678 quality[i] = 1./(1.00001-v0->GetV0CosineOfPointingAngle()); //base point angle
5679 // quality[i] /= (0.5+v0->GetDist2());
5680 // quality[i] *= v0->GetChi2After(); //density factor
5682 Int_t index0 = v0->GetIndex(0);
5683 Int_t index1 = v0->GetIndex(1);
5684 //Bo: Double_t minpulldca = TMath::Min(2.+pulldca[v0->GetIndex(0)],(2.+pulldca[v0->GetIndex(1)]) ); //pull
5685 Double_t minpulldca = TMath::Min(2.+pulldca[index0],(2.+pulldca[index1]) );//Bo:
5689 AliTPCseed * track0 = (AliTPCseed*)array->At(index0);
5690 AliTPCseed * track1 = (AliTPCseed*)array->At(index1);
5691 if (track0->TPCrPID(0)>0.3&&track1->TPCrPID(0)>0.3&&v0->GetAnglep()[2]<0.15) quality[i]+=1000000; // gamma conversion candidate
5692 if (track0->TPCrPID(4)>0.9||(track1->TPCrPID(4)>0.9&&minpulldca>4)) quality[i]*=10; // lambda candidate candidate
5695 TMath::Sort(ncandidates,quality,indexes,kTRUE);
5698 for (Int_t i=0;i<ncandidates;i++){
5699 AliESDv0 * v0 = (AliESDv0*)tpcv0s->At(indexes[i]);
5701 Int_t index0 = v0->GetIndex(0);
5702 Int_t index1 = v0->GetIndex(1);
5703 AliTPCseed * track0 = (AliTPCseed*)array->At(index0);
5704 AliTPCseed * track1 = (AliTPCseed*)array->At(index1);
5705 if (!track0||!track1) {
5709 Bool_t accept =kTRUE; //default accept
5710 Int_t *v0indexes0 = track0->GetV0Indexes();
5711 Int_t *v0indexes1 = track1->GetV0Indexes();
5713 Int_t order0 = (v0indexes0[0]!=0) ? 1:0;
5714 Int_t order1 = (v0indexes1[0]!=0) ? 1:0;
5715 if (v0indexes0[1]!=0) order0 =2;
5716 if (v0indexes1[1]!=0) order1 =2;
5718 if (v0indexes0[2]!=0) {order0=3; accept=kFALSE;}
5719 if (v0indexes0[2]!=0) {order1=3; accept=kFALSE;}
5721 AliESDv0 * v02 = v0;
5723 //Bo: v0->SetOrder(0,order0);
5724 //Bo: v0->SetOrder(1,order1);
5725 //Bo: v0->SetOrder(1,order0+order1);
5726 v0->SetOnFlyStatus(kTRUE);
5727 Int_t index = esd->AddV0(v0);
5728 v02 = esd->GetV0(index);
5729 v0indexes0[order0]=index;
5730 v0indexes1[order1]=index;
5734 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
5735 if (AliTPCReconstructor::StreamLevel()>1) {
5736 Int_t lab0=track0->GetLabel();
5737 Int_t lab1=track1->GetLabel();
5738 TTreeSRedirector &cstream = *fDebugStreamer;
5747 "dca0="<<dca[index0]<<
5748 "dca1="<<dca[index1]<<
5752 "quality="<<quality[i]<<
5753 "pulldca0="<<pulldca[index0]<<
5754 "pulldca1="<<pulldca[index1]<<
5778 printf("TPC V0 finder : naccepted\t%d\tncandidates\t%d\tntracks\t%d\tnall\t%d\n",naccepted,ncandidates,ntracks,nall);
5782 Int_t AliTPCtrackerMI::RefitKink(AliTPCseed &mother, AliTPCseed &daughter, const AliESDkink &knk)
5785 // refit kink towards to the vertex
5788 AliKink &kink=(AliKink &)knk;
5790 Int_t row0 = GetRowNumber(kink.GetR());
5791 FollowProlongation(mother,0);
5792 mother.Reset(kFALSE);
5794 FollowProlongation(daughter,row0);
5795 daughter.Reset(kFALSE);
5796 FollowBackProlongation(daughter,158);
5797 daughter.Reset(kFALSE);
5798 Int_t first = TMath::Max(row0-20,30);
5799 Int_t last = TMath::Min(row0+20,140);
5801 const Int_t kNdiv =5;
5802 AliTPCseed param0[kNdiv]; // parameters along the track
5803 AliTPCseed param1[kNdiv]; // parameters along the track
5804 AliKink kinks[kNdiv]; // corresponding kink parameters
5807 for (Int_t irow=0; irow<kNdiv;irow++){
5808 rows[irow] = first +((last-first)*irow)/(kNdiv-1);
5810 // store parameters along the track
5812 for (Int_t irow=0;irow<kNdiv;irow++){
5813 FollowBackProlongation(mother, rows[irow]);
5814 FollowProlongation(daughter,rows[kNdiv-1-irow]);
5815 param0[irow] = mother;
5816 param1[kNdiv-1-irow] = daughter;
5820 for (Int_t irow=0; irow<kNdiv-1;irow++){
5821 if (param0[irow].GetNumberOfClusters()<kNdiv||param1[irow].GetNumberOfClusters()<kNdiv) continue;
5822 kinks[irow].SetMother(param0[irow]);
5823 kinks[irow].SetDaughter(param1[irow]);
5824 kinks[irow].Update();
5827 // choose kink with best "quality"
5829 Double_t mindist = 10000;
5830 for (Int_t irow=0;irow<kNdiv;irow++){
5831 if (param0[irow].GetNumberOfClusters()<20||param1[irow].GetNumberOfClusters()<20) continue;
5832 if (TMath::Abs(kinks[irow].GetR())>240.) continue;
5833 if (TMath::Abs(kinks[irow].GetR())<100.) continue;
5835 Float_t normdist = TMath::Abs(param0[irow].GetX()-kinks[irow].GetR())*(0.1+kink.GetDistance());
5836 normdist/= (param0[irow].GetNumberOfClusters()+param1[irow].GetNumberOfClusters()+40.);
5837 if (normdist < mindist){
5843 if (index==-1) return 0;
5846 param0[index].Reset(kTRUE);
5847 FollowProlongation(param0[index],0);
5849 mother = param0[index];
5850 daughter = param1[index]; // daughter in vertex
5852 kink.SetMother(mother);
5853 kink.SetDaughter(daughter);
5855 kink.SetTPCRow0(GetRowNumber(kink.GetR()));
5856 kink.SetTPCncls(param0[index].GetNumberOfClusters(),0);
5857 kink.SetTPCncls(param1[index].GetNumberOfClusters(),1);
5858 kink.SetLabel(CookLabel(&mother,0.4, 0,kink.GetTPCRow0()),0);
5859 kink.SetLabel(CookLabel(&daughter,0.4, kink.GetTPCRow0(),160),1);
5860 mother.SetLabel(kink.GetLabel(0));
5861 daughter.SetLabel(kink.GetLabel(1));
5867 void AliTPCtrackerMI::UpdateKinkQualityM(AliTPCseed * seed){
5869 // update Kink quality information for mother after back propagation
5871 if (seed->GetKinkIndex(0)>=0) return;
5872 for (Int_t ikink=0;ikink<3;ikink++){
5873 Int_t index = seed->GetKinkIndex(ikink);
5874 if (index>=0) break;
5875 index = TMath::Abs(index)-1;
5876 AliESDkink * kink = fEvent->GetKink(index);
5877 kink->SetTPCDensity(-1,0,0);
5878 kink->SetTPCDensity(1,0,1);
5880 Int_t row0 = kink->GetTPCRow0() - 2 - Int_t( 0.5/ (0.05+kink->GetAngle(2)));
5881 if (row0<15) row0=15;
5883 Int_t row1 = kink->GetTPCRow0() + 2 + Int_t( 0.5/ (0.05+kink->GetAngle(2)));
5884 if (row1>145) row1=145;
5886 Int_t found,foundable,shared;
5887 seed->GetClusterStatistic(0,row0, found, foundable,shared,kFALSE);
5888 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),0,0);
5889 seed->GetClusterStatistic(row1,155, found, foundable,shared,kFALSE);
5890 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),0,1);
5895 void AliTPCtrackerMI::UpdateKinkQualityD(AliTPCseed * seed){
5897 // update Kink quality information for daughter after refit
5899 if (seed->GetKinkIndex(0)<=0) return;
5900 for (Int_t ikink=0;ikink<3;ikink++){
5901 Int_t index = seed->GetKinkIndex(ikink);
5902 if (index<=0) break;
5903 index = TMath::Abs(index)-1;
5904 AliESDkink * kink = fEvent->GetKink(index);
5905 kink->SetTPCDensity(-1,1,0);
5906 kink->SetTPCDensity(-1,1,1);
5908 Int_t row0 = kink->GetTPCRow0() -2 - Int_t( 0.5/ (0.05+kink->GetAngle(2)));
5909 if (row0<15) row0=15;
5911 Int_t row1 = kink->GetTPCRow0() +2 + Int_t( 0.5/ (0.05+kink->GetAngle(2)));
5912 if (row1>145) row1=145;
5914 Int_t found,foundable,shared;
5915 seed->GetClusterStatistic(0,row0, found, foundable,shared,kFALSE);
5916 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),1,0);
5917 seed->GetClusterStatistic(row1,155, found, foundable,shared,kFALSE);
5918 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),1,1);
5924 Int_t AliTPCtrackerMI::CheckKinkPoint(AliTPCseed*seed,AliTPCseed &mother, AliTPCseed &daughter, const AliESDkink &knk)
5927 // check kink point for given track
5928 // if return value=0 kink point not found
5929 // otherwise seed0 correspond to mother particle
5930 // seed1 correspond to daughter particle
5931 // kink parameter of kink point
5932 AliKink &kink=(AliKink &)knk;
5934 Int_t middlerow = (seed->GetFirstPoint()+seed->GetLastPoint())/2;
5935 Int_t first = seed->GetFirstPoint();
5936 Int_t last = seed->GetLastPoint();
5937 if (last-first<20) return 0; // shortest length - 2*30 = 60 pad-rows
5940 AliTPCseed *seed1 = ReSeed(seed,middlerow+20, kTRUE); //middle of chamber
5941 if (!seed1) return 0;
5942 FollowProlongation(*seed1,seed->GetLastPoint()-20);
5943 seed1->Reset(kTRUE);
5944 FollowProlongation(*seed1,158);
5945 seed1->Reset(kTRUE);
5946 last = seed1->GetLastPoint();
5948 AliTPCseed *seed0 = new AliTPCseed(*seed);
5949 seed0->Reset(kFALSE);
5952 AliTPCseed param0[20]; // parameters along the track
5953 AliTPCseed param1[20]; // parameters along the track
5954 AliKink kinks[20]; // corresponding kink parameters
5956 for (Int_t irow=0; irow<20;irow++){
5957 rows[irow] = first +((last-first)*irow)/19;
5959 // store parameters along the track
5961 for (Int_t irow=0;irow<20;irow++){
5962 FollowBackProlongation(*seed0, rows[irow]);
5963 FollowProlongation(*seed1,rows[19-irow]);
5964 param0[irow] = *seed0;
5965 param1[19-irow] = *seed1;
5969 for (Int_t irow=0; irow<19;irow++){
5970 kinks[irow].SetMother(param0[irow]);
5971 kinks[irow].SetDaughter(param1[irow]);
5972 kinks[irow].Update();
5975 // choose kink with biggest change of angle
5977 Double_t maxchange= 0;
5978 for (Int_t irow=1;irow<19;irow++){
5979 if (TMath::Abs(kinks[irow].GetR())>240.) continue;
5980 if (TMath::Abs(kinks[irow].GetR())<110.) continue;
5981 Float_t quality = TMath::Abs(kinks[irow].GetAngle(2))/(3.+TMath::Abs(kinks[irow].GetR()-param0[irow].GetX()));
5982 if ( quality > maxchange){
5983 maxchange = quality;
5990 if (index<0) return 0;
5992 Int_t row0 = GetRowNumber(kinks[index].GetR()); //row 0 estimate
5993 seed0 = new AliTPCseed(param0[index]);
5994 seed1 = new AliTPCseed(param1[index]);
5995 seed0->Reset(kFALSE);
5996 seed1->Reset(kFALSE);
5997 seed0->ResetCovariance(10.);
5998 seed1->ResetCovariance(10.);
5999 FollowProlongation(*seed0,0);
6000 FollowBackProlongation(*seed1,158);
6001 mother = *seed0; // backup mother at position 0
6002 seed0->Reset(kFALSE);
6003 seed1->Reset(kFALSE);
6004 seed0->ResetCovariance(10.);
6005 seed1->ResetCovariance(10.);
6007 first = TMath::Max(row0-20,0);
6008 last = TMath::Min(row0+20,158);
6010 for (Int_t irow=0; irow<20;irow++){
6011 rows[irow] = first +((last-first)*irow)/19;
6013 // store parameters along the track
6015 for (Int_t irow=0;irow<20;irow++){
6016 FollowBackProlongation(*seed0, rows[irow]);
6017 FollowProlongation(*seed1,rows[19-irow]);
6018 param0[irow] = *seed0;
6019 param1[19-irow] = *seed1;
6023 for (Int_t irow=0; irow<19;irow++){
6024 kinks[irow].SetMother(param0[irow]);
6025 kinks[irow].SetDaughter(param1[irow]);
6026 // param0[irow].Dump();
6027 //param1[irow].Dump();
6028 kinks[irow].Update();
6031 // choose kink with biggest change of angle
6034 for (Int_t irow=0;irow<20;irow++){
6035 if (TMath::Abs(kinks[irow].GetR())>250.) continue;
6036 if (TMath::Abs(kinks[irow].GetR())<90.) continue;
6037 Float_t quality = TMath::Abs(kinks[irow].GetAngle(2))/(3.+TMath::Abs(kinks[irow].GetR()-param0[irow].GetX()));
6038 if ( quality > maxchange){
6039 maxchange = quality;
6046 if (index==-1 || param0[index].GetNumberOfClusters()+param1[index].GetNumberOfClusters()<100){
6052 // Float_t anglesigma = TMath::Sqrt(param0[index].fC22+param0[index].fC33+param1[index].fC22+param1[index].fC33);
6054 kink.SetMother(param0[index]);
6055 kink.SetDaughter(param1[index]);
6058 Double_t chi2P2 = param0[index].GetParameter()[2]-param1[index].GetParameter()[2];
6060 chi2P2/=param0[index].GetCovariance()[5]+param1[index].GetCovariance()[5];
6061 Double_t chi2P3 = param0[index].GetParameter()[3]-param1[index].GetParameter()[3];
6063 chi2P3/=param0[index].GetCovariance()[9]+param1[index].GetCovariance()[9];
6065 if (AliTPCReconstructor::StreamLevel()>1) {
6066 (*fDebugStreamer)<<"kinkHpt"<<
6069 "p0.="<<¶m0[index]<<
6070 "p1.="<<¶m1[index]<<
6074 if ( chi2P2+chi2P3<AliTPCReconstructor::GetRecoParam()->GetKinkAngleCutChi2(0)){
6081 row0 = GetRowNumber(kink.GetR());
6082 kink.SetTPCRow0(row0);
6083 kink.SetLabel(CookLabel(seed0,0.5,0,row0),0);
6084 kink.SetLabel(CookLabel(seed1,0.5,row0,158),1);
6085 kink.SetIndex(-10,0);
6086 kink.SetIndex(int(param0[index].GetNumberOfClusters()+param1[index].GetNumberOfClusters()),1);
6087 kink.SetTPCncls(param0[index].GetNumberOfClusters(),0);
6088 kink.SetTPCncls(param1[index].GetNumberOfClusters(),1);
6091 // new (&mother) AliTPCseed(param0[index]);
6092 daughter = param1[index];
6093 daughter.SetLabel(kink.GetLabel(1));
6094 param0[index].Reset(kTRUE);
6095 FollowProlongation(param0[index],0);
6096 mother = param0[index];
6097 mother.SetLabel(kink.GetLabel(0));
6098 if ( chi2P2+chi2P3<AliTPCReconstructor::GetRecoParam()->GetKinkAngleCutChi2(1)){
6110 AliTPCseed* AliTPCtrackerMI::ReSeed(AliTPCseed *t)
6113 // reseed - refit - track
6116 // Int_t last = fSectors->GetNRows()-1;
6118 if (fSectors == fOuterSec){
6119 first = TMath::Max(first, t->GetFirstPoint()-fInnerSec->GetNRows());
6123 first = t->GetFirstPoint();
6125 AliTPCseed * seed = MakeSeed(t,0.1,0.5,0.9);
6126 FollowBackProlongation(*t,fSectors->GetNRows()-1);
6128 FollowProlongation(*t,first);
6138 //_____________________________________________________________________________
6139 Int_t AliTPCtrackerMI::ReadSeeds(const TFile *inp) {
6140 //-----------------------------------------------------------------
6141 // This function reades track seeds.
6142 //-----------------------------------------------------------------
6143 TDirectory *savedir=gDirectory;
6145 TFile *in=(TFile*)inp;
6146 if (!in->IsOpen()) {
6147 cerr<<"AliTPCtrackerMI::ReadSeeds(): input file is not open !\n";
6152 TTree *seedTree=(TTree*)in->Get("Seeds");
6154 cerr<<"AliTPCtrackerMI::ReadSeeds(): ";
6155 cerr<<"can't get a tree with track seeds !\n";
6158 AliTPCtrack *seed=new AliTPCtrack;
6159 seedTree->SetBranchAddress("tracks",&seed);
6161 if (fSeeds==0) fSeeds=new TObjArray(15000);
6163 Int_t n=(Int_t)seedTree->GetEntries();
6164 for (Int_t i=0; i<n; i++) {
6165 seedTree->GetEvent(i);
6166 fSeeds->AddLast(new AliTPCseed(*seed/*,seed->GetAlpha()*/));
6175 Int_t AliTPCtrackerMI::Clusters2Tracks (AliESDEvent *const esd)
6178 if (fSeeds) DeleteSeeds();
6181 if (!fSeeds) return 1;
6183 if (AliTPCReconstructor::StreamLevel()>3) DumpClusters(0,fSeeds);
6189 //_____________________________________________________________________________
6190 Int_t AliTPCtrackerMI::Clusters2Tracks() {
6191 //-----------------------------------------------------------------
6192 // This is a track finder.
6193 //-----------------------------------------------------------------
6194 TDirectory *savedir=gDirectory;
6198 fSeeds = Tracking();
6201 Info("Clusters2Tracks","Time for tracking: \t");timer.Print();timer.Start();
6203 //activate again some tracks
6204 for (Int_t i=0; i<fSeeds->GetEntriesFast(); i++) {
6205 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
6207 Int_t nc=t.GetNumberOfClusters();
6209 delete fSeeds->RemoveAt(i);
6213 if (pt->GetRemoval()==10) {
6214 if (pt->GetDensityFirst(20)>0.8 || pt->GetDensityFirst(30)>0.8 || pt->GetDensityFirst(40)>0.7)
6215 pt->Desactivate(10); // make track again active
6217 pt->Desactivate(20);
6218 delete fSeeds->RemoveAt(i);
6223 RemoveUsed2(fSeeds,0.85,0.85,0);
6224 if (AliTPCReconstructor::GetRecoParam()->GetDoKinks()) FindKinks(fSeeds,fEvent);
6225 //FindCurling(fSeeds, fEvent,0);
6226 if (AliTPCReconstructor::StreamLevel()>5) FindMultiMC(fSeeds, fEvent,-1); // find multi found tracks
6227 RemoveUsed2(fSeeds,0.5,0.4,20);
6228 FindSplitted(fSeeds, fEvent,0); // find multi found tracks
6229 if (AliTPCReconstructor::StreamLevel()>5) FindMultiMC(fSeeds, fEvent,0); // find multi found tracks
6232 // // refit short tracks
6234 Int_t nseed=fSeeds->GetEntriesFast();
6237 for (Int_t i=0; i<nseed; i++) {
6238 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
6240 Int_t nc=t.GetNumberOfClusters();
6242 delete fSeeds->RemoveAt(i);
6245 CookLabel(pt,0.1); //For comparison only
6246 //if ((pt->IsActive() || (pt->fRemoval==10) )&& nc>50 &&pt->GetNumberOfClusters()>0.4*pt->fNFoundable){
6247 if ((pt->IsActive() || (pt->GetRemoval()==10) )){
6249 if (fDebug>0) cerr<<found<<'\r';
6253 delete fSeeds->RemoveAt(i);
6257 //RemoveOverlap(fSeeds,0.99,7,kTRUE);
6259 //RemoveUsed(fSeeds,0.9,0.9,6);
6261 nseed=fSeeds->GetEntriesFast();
6263 for (Int_t i=0; i<nseed; i++) {
6264 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
6266 Int_t nc=t.GetNumberOfClusters();
6268 delete fSeeds->RemoveAt(i);
6272 t.CookdEdx(0.02,0.6);
6273 // CheckKinkPoint(&t,0.05);
6274 //if ((pt->IsActive() || (pt->fRemoval==10) )&& nc>50 &&pt->GetNumberOfClusters()>0.4*pt->fNFoundable){
6275 if ((pt->IsActive() || (pt->GetRemoval()==10) )){
6283 delete fSeeds->RemoveAt(i);
6284 //AliTPCseed * seed1 = ReSeed(pt,0.05,0.5,1);
6286 // FollowProlongation(*seed1,0);
6287 // Int_t n = seed1->GetNumberOfClusters();
6288 // printf("fP4\t%f\t%f\n",seed1->GetC(),pt->GetC());
6289 // printf("fN\t%d\t%d\n", seed1->GetNumberOfClusters(),pt->GetNumberOfClusters());
6292 //AliTPCseed * seed2 = ReSeed(pt,0.95,0.5,0.05);
6296 SortTracks(fSeeds, 1);
6300 PrepareForBackProlongation(fSeeds,5.);
6301 PropagateBack(fSeeds);
6302 printf("Time for back propagation: \t");timer.Print();timer.Start();
6306 PrepareForProlongation(fSeeds,5.);
6307 PropagateForward2(fSeeds);
6309 printf("Time for FORWARD propagation: \t");timer.Print();timer.Start();
6310 // RemoveUsed(fSeeds,0.7,0.7,6);
6311 //RemoveOverlap(fSeeds,0.9,7,kTRUE);
6313 nseed=fSeeds->GetEntriesFast();
6315 for (Int_t i=0; i<nseed; i++) {
6316 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
6318 Int_t nc=t.GetNumberOfClusters();
6320 delete fSeeds->RemoveAt(i);
6323 t.CookdEdx(0.02,0.6);
6324 // CookLabel(pt,0.1); //For comparison only
6325 //if ((pt->IsActive() || (pt->fRemoval==10) )&& nc>50 &&pt->GetNumberOfClusters()>0.4*pt->fNFoundable){
6326 if ((pt->IsActive() || (pt->fRemoval==10) )){
6327 cerr<<found++<<'\r';
6330 delete fSeeds->RemoveAt(i);
6335 // fNTracks = found;
6337 Info("Clusters2Tracks","Time for overlap removal, track writing and dedx cooking: \t"); timer.Print();timer.Start();
6340 // cerr<<"Number of found tracks : "<<"\t"<<found<<endl;
6341 Info("Clusters2Tracks","Number of found tracks %d",found);
6343 // UnloadClusters();
6348 void AliTPCtrackerMI::Tracking(TObjArray * arr)
6351 // tracking of the seeds
6354 fSectors = fOuterSec;
6355 ParallelTracking(arr,150,63);
6356 fSectors = fOuterSec;
6357 ParallelTracking(arr,63,0);
6360 TObjArray * AliTPCtrackerMI::Tracking(Int_t seedtype, Int_t i1, Int_t i2, Float_t cuts[4], Float_t dy, Int_t dsec)
6365 TObjArray * arr = new TObjArray;
6367 fSectors = fOuterSec;
6370 for (Int_t sec=0;sec<fkNOS;sec++){
6371 if (seedtype==3) MakeSeeds3(arr,sec,i1,i2,cuts,dy, dsec);
6372 if (seedtype==4) MakeSeeds5(arr,sec,i1,i2,cuts,dy);
6373 if (seedtype==2) MakeSeeds2(arr,sec,i1,i2,cuts,dy);
6376 Info("Tracking","\nSeeding - %d\t%d\t%d\t%d\n",seedtype,i1,i2,arr->GetEntriesFast());
6388 TObjArray * AliTPCtrackerMI::Tracking()
6392 if (AliTPCReconstructor::GetRecoParam()->GetSpecialSeeding()) return TrackingSpecial();
6395 Int_t nup=fOuterSec->GetNRows()+fInnerSec->GetNRows();
6397 TObjArray * seeds = new TObjArray;
6399 Int_t fLastSeedRowSec=AliTPCReconstructor::GetRecoParam()->GetLastSeedRowSec();
6400 Int_t gapPrim = AliTPCReconstructor::GetRecoParam()->GetSeedGapPrim();
6401 Int_t gapSec = AliTPCReconstructor::GetRecoParam()->GetSeedGapSec();
6409 Float_t fnumber = 3.0;
6410 Float_t fdensity = 3.0;
6415 for (Int_t delta = 0; delta<18; delta+=gapPrim){
6419 arr = Tracking(3,nup-1-delta,nup-1-delta-gap,cuts,-1,1);
6420 SumTracks(seeds,arr);
6421 SignClusters(seeds,fnumber,fdensity);
6423 for (Int_t i=2;i<6;i+=2){
6424 // seed high pt tracks
6427 arr = Tracking(3,nup-i-delta,nup-i-delta-gap,cuts,-1,0);
6428 SumTracks(seeds,arr);
6429 SignClusters(seeds,fnumber,fdensity);
6434 // RemoveUsed(seeds,0.9,0.9,1);
6435 // UnsignClusters();
6436 // SignClusters(seeds,fnumber,fdensity);
6440 for (Int_t delta = 20; delta<120; delta+=gapPrim){
6442 // seed high pt tracks
6446 arr = Tracking(3,nup-delta,nup-delta-gap,cuts,-1);
6447 SumTracks(seeds,arr);
6448 SignClusters(seeds,fnumber,fdensity);
6453 arr = Tracking(3,nup-delta-5,nup-delta-5-gap,cuts,-1);
6454 SumTracks(seeds,arr);
6455 SignClusters(seeds,fnumber,fdensity);
6466 Info("Tracking()","\n\nPrimary seeding\t%d\n\n",seeds->GetEntriesFast());
6470 // RemoveUsed(seeds,0.75,0.75,1);
6472 //SignClusters(seeds,fnumber,fdensity);
6481 arr = Tracking(4,nup-1,nup-1-gap,cuts,-1);
6482 SumTracks(seeds,arr);
6483 SignClusters(seeds,fnumber,fdensity);
6485 arr = Tracking(4,nup-2,nup-2-gap,cuts,-1);
6486 SumTracks(seeds,arr);
6487 SignClusters(seeds,fnumber,fdensity);
6489 arr = Tracking(4,nup-3,nup-3-gap,cuts,-1);
6490 SumTracks(seeds,arr);
6491 SignClusters(seeds,fnumber,fdensity);
6493 arr = Tracking(4,nup-5,nup-5-gap,cuts,-1);
6494 SumTracks(seeds,arr);
6495 SignClusters(seeds,fnumber,fdensity);
6497 arr = Tracking(4,nup-7,nup-7-gap,cuts,-1);
6498 SumTracks(seeds,arr);
6499 SignClusters(seeds,fnumber,fdensity);
6502 arr = Tracking(4,nup-9,nup-9-gap,cuts,-1);
6503 SumTracks(seeds,arr);
6504 SignClusters(seeds,fnumber,fdensity);
6508 for (Int_t delta = 9; delta<30; delta+=gapSec){
6514 arr = Tracking(4,nup-1-delta,nup-1-delta-gap,cuts,-1);
6515 SumTracks(seeds,arr);
6516 SignClusters(seeds,fnumber,fdensity);
6518 arr = Tracking(4,nup-3-delta,nup-5-delta-gap,cuts,4);
6519 SumTracks(seeds,arr);
6520 SignClusters(seeds,fnumber,fdensity);
6533 for (Int_t delta = 30; delta<fLastSeedRowSec; delta+=gapSec){
6539 arr = Tracking(4,nup-1-delta,nup-1-delta-gap,cuts,-1);
6540 SumTracks(seeds,arr);
6541 SignClusters(seeds,fnumber,fdensity);
6543 arr = Tracking(4,nup-5-delta,nup-5-delta-gap,cuts,5 );
6544 SumTracks(seeds,arr);
6545 SignClusters(seeds,fnumber,fdensity);
6549 Info("Tracking()","\n\nSecondary seeding\t%d\n\n",seeds->GetEntriesFast());
6560 TObjArray * AliTPCtrackerMI::TrackingSpecial()
6563 // seeding adjusted for laser and cosmic tests - short tracks with big inclination angle
6564 // no primary vertex seeding tried
6568 Int_t nup=fOuterSec->GetNRows()+fInnerSec->GetNRows();
6570 TObjArray * seeds = new TObjArray;
6575 Float_t fnumber = 3.0;
6576 Float_t fdensity = 3.0;
6579 cuts[0] = AliTPCReconstructor::GetRecoParam()->GetMaxC(); // max curvature
6580 cuts[1] = 3.5; // max tan(phi) angle for seeding
6581 cuts[2] = 3.; // not used (cut on z primary vertex)
6582 cuts[3] = 3.5; // max tan(theta) angle for seeding
6584 for (Int_t delta = 0; nup-delta-gap-1>0; delta+=3){
6586 arr = Tracking(4,nup-1-delta,nup-1-delta-gap,cuts,-1);
6587 SumTracks(seeds,arr);
6588 SignClusters(seeds,fnumber,fdensity);
6592 Info("Tracking()","\n\nSecondary seeding\t%d\n\n",seeds->GetEntriesFast());
6603 void AliTPCtrackerMI::SumTracks(TObjArray *arr1,TObjArray *&arr2) const
6606 //sum tracks to common container
6607 //remove suspicious tracks
6608 Int_t nseed = arr2->GetEntriesFast();
6609 for (Int_t i=0;i<nseed;i++){
6610 AliTPCseed *pt=(AliTPCseed*)arr2->UncheckedAt(i);
6613 // remove tracks with too big curvature
6615 if (TMath::Abs(pt->GetC())>AliTPCReconstructor::GetRecoParam()->GetMaxC()){
6616 delete arr2->RemoveAt(i);
6619 // REMOVE VERY SHORT TRACKS
6620 if (pt->GetNumberOfClusters()<20){
6621 delete arr2->RemoveAt(i);
6624 // NORMAL ACTIVE TRACK
6625 if (pt->IsActive()){
6626 arr1->AddLast(arr2->RemoveAt(i));
6629 //remove not usable tracks
6630 if (pt->GetRemoval()!=10){
6631 delete arr2->RemoveAt(i);
6635 // ENABLE ONLY ENOUGH GOOD STOPPED TRACKS
6636 if (pt->GetDensityFirst(20)>0.8 || pt->GetDensityFirst(30)>0.8 || pt->GetDensityFirst(40)>0.7)
6637 arr1->AddLast(arr2->RemoveAt(i));
6639 delete arr2->RemoveAt(i);
6643 delete arr2; arr2 = 0;
6648 void AliTPCtrackerMI::ParallelTracking(TObjArray *const arr, Int_t rfirst, Int_t rlast)
6651 // try to track in parralel
6653 Int_t nseed=arr->GetEntriesFast();
6654 //prepare seeds for tracking
6655 for (Int_t i=0; i<nseed; i++) {
6656 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i), &t=*pt;
6658 if (!t.IsActive()) continue;
6659 // follow prolongation to the first layer
6660 if ( (fSectors ==fInnerSec) || (t.GetFirstPoint()-fkParam->GetNRowLow()>rfirst+1) )
6661 FollowProlongation(t, rfirst+1);
6666 for (Int_t nr=rfirst; nr>=rlast; nr--){
6667 if (nr<fInnerSec->GetNRows())
6668 fSectors = fInnerSec;
6670 fSectors = fOuterSec;
6671 // make indexes with the cluster tracks for given
6673 // find nearest cluster
6674 for (Int_t i=0; i<nseed; i++) {
6675 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i), &t=*pt;
6677 if (nr==80) pt->UpdateReference();
6678 if (!pt->IsActive()) continue;
6679 // if ( (fSectors ==fOuterSec) && (pt->fFirstPoint-fkParam->GetNRowLow())<nr) continue;
6680 if (pt->GetRelativeSector()>17) {
6683 UpdateClusters(t,nr);
6685 // prolonagate to the nearest cluster - if founded
6686 for (Int_t i=0; i<nseed; i++) {
6687 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
6689 if (!pt->IsActive()) continue;
6690 // if ((fSectors ==fOuterSec) && (pt->fFirstPoint-fkParam->GetNRowLow())<nr) continue;
6691 if (pt->GetRelativeSector()>17) {
6694 FollowToNextCluster(*pt,nr);
6699 void AliTPCtrackerMI::PrepareForBackProlongation(TObjArray *const arr,Float_t fac) const
6703 // if we use TPC track itself we have to "update" covariance
6705 Int_t nseed= arr->GetEntriesFast();
6706 for (Int_t i=0;i<nseed;i++){
6707 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
6711 //rotate to current local system at first accepted point
6712 Int_t index = pt->GetClusterIndex2(pt->GetFirstPoint());
6713 Int_t sec = (index&0xff000000)>>24;
6715 Float_t angle1 = fInnerSec->GetAlpha()*sec+fInnerSec->GetAlphaShift();
6716 if (angle1>TMath::Pi())
6717 angle1-=2.*TMath::Pi();
6718 Float_t angle2 = pt->GetAlpha();
6720 if (TMath::Abs(angle1-angle2)>0.001){
6721 pt->Rotate(angle1-angle2);
6722 //angle2 = pt->GetAlpha();
6723 //pt->fRelativeSector = pt->GetAlpha()/fInnerSec->GetAlpha();
6724 //if (pt->GetAlpha()<0)
6725 // pt->fRelativeSector+=18;
6726 //sec = pt->fRelativeSector;
6735 void AliTPCtrackerMI::PrepareForProlongation(TObjArray *const arr, Float_t fac) const
6739 // if we use TPC track itself we have to "update" covariance
6741 Int_t nseed= arr->GetEntriesFast();
6742 for (Int_t i=0;i<nseed;i++){
6743 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
6746 pt->SetFirstPoint(pt->GetLastPoint());
6754 Int_t AliTPCtrackerMI::PropagateBack(TObjArray *const arr)
6757 // make back propagation
6759 Int_t nseed= arr->GetEntriesFast();
6760 for (Int_t i=0;i<nseed;i++){
6761 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
6762 if (pt&& pt->GetKinkIndex(0)<=0) {
6763 //AliTPCseed *pt2 = new AliTPCseed(*pt);
6764 fSectors = fInnerSec;
6765 //FollowBackProlongation(*pt,fInnerSec->GetNRows()-1);
6766 //fSectors = fOuterSec;
6767 FollowBackProlongation(*pt,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1);
6768 //if (pt->GetNumberOfClusters()<(pt->fEsd->GetTPCclusters(0)) ){
6769 // Error("PropagateBack","Not prolonged track %d",pt->GetLabel());
6770 // FollowBackProlongation(*pt2,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1);
6773 if (pt&& pt->GetKinkIndex(0)>0) {
6774 AliESDkink * kink = fEvent->GetKink(pt->GetKinkIndex(0)-1);
6775 pt->SetFirstPoint(kink->GetTPCRow0());
6776 fSectors = fInnerSec;
6777 FollowBackProlongation(*pt,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1);
6785 Int_t AliTPCtrackerMI::PropagateForward2(TObjArray *const arr)
6788 // make forward propagation
6790 Int_t nseed= arr->GetEntriesFast();
6792 for (Int_t i=0;i<nseed;i++){
6793 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
6795 FollowProlongation(*pt,0);
6804 Int_t AliTPCtrackerMI::PropagateForward()
6807 // propagate track forward
6809 Int_t nseed = fSeeds->GetEntriesFast();
6810 for (Int_t i=0;i<nseed;i++){
6811 AliTPCseed *pt = (AliTPCseed*)fSeeds->UncheckedAt(i);
6813 AliTPCseed &t = *pt;
6814 Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
6815 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
6816 if (alpha < 0. ) alpha += 2.*TMath::Pi();
6817 t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN);
6821 fSectors = fOuterSec;
6822 ParallelTracking(fSeeds,fOuterSec->GetNRows()+fInnerSec->GetNRows()-1,fInnerSec->GetNRows());
6823 fSectors = fInnerSec;
6824 ParallelTracking(fSeeds,fInnerSec->GetNRows()-1,0);
6833 Int_t AliTPCtrackerMI::PropagateBack(AliTPCseed *const pt, Int_t row0, Int_t row1)
6836 // make back propagation, in between row0 and row1
6840 fSectors = fInnerSec;
6843 if (row1<fSectors->GetNRows())
6846 r1 = fSectors->GetNRows()-1;
6848 if (row0<fSectors->GetNRows()&& r1>0 )
6849 FollowBackProlongation(*pt,r1);
6850 if (row1<=fSectors->GetNRows())
6853 r1 = row1 - fSectors->GetNRows();
6854 if (r1<=0) return 0;
6855 if (r1>=fOuterSec->GetNRows()) return 0;
6856 fSectors = fOuterSec;
6857 return FollowBackProlongation(*pt,r1);
6865 void AliTPCtrackerMI::GetShape(AliTPCseed * seed, Int_t row)
6869 AliTPCClusterParam * clparam = AliTPCcalibDB::Instance()->GetClusterParam();
6870 Float_t zdrift = TMath::Abs((fkParam->GetZLength(0)-TMath::Abs(seed->GetZ())));
6871 Int_t type = (seed->GetSector() < fkParam->GetNSector()/2) ? 0: (row>126) ? 1:2;
6872 Double_t angulary = seed->GetSnp();
6873 angulary = angulary*angulary/((1.-angulary)*(1.+angulary));
6874 Double_t angularz = seed->GetTgl()*seed->GetTgl()*(1.+angulary);
6876 Double_t sigmay = clparam->GetRMS0(0,type,zdrift,TMath::Sqrt(TMath::Abs(angulary)));
6877 Double_t sigmaz = clparam->GetRMS0(1,type,zdrift,TMath::Sqrt(TMath::Abs(angularz)));
6878 seed->SetCurrentSigmaY2(sigmay*sigmay);
6879 seed->SetCurrentSigmaZ2(sigmaz*sigmaz);
6880 // Float_t sd2 = TMath::Abs((fkParam->GetZLength(0)-TMath::Abs(seed->GetZ())))*fkParam->GetDiffL()*fkParam->GetDiffL();
6881 // // Float_t padlength = fkParam->GetPadPitchLength(seed->fSector);
6882 // Float_t padlength = GetPadPitchLength(row);
6884 // Float_t sresy = (seed->GetSector() < fkParam->GetNSector()/2) ? 0.2 :0.3;
6885 // seed->SetCurrentSigmaY2(sd2+padlength*padlength*angulary/12.+sresy*sresy);
6887 // Float_t sresz = fkParam->GetZSigma();
6888 // seed->SetCurrentSigmaZ2(sd2+padlength*padlength*angularz*angularz*(1+angulary)/12.+sresz*sresz);
6890 Float_t wy = GetSigmaY(seed);
6891 Float_t wz = GetSigmaZ(seed);
6894 if (TMath::Abs(wy/seed->fCurrentSigmaY2-1)>0.0001 || TMath::Abs(wz/seed->fCurrentSigmaZ2-1)>0.0001 ){
6895 printf("problem\n");
6902 //__________________________________________________________________________
6903 void AliTPCtrackerMI::CookLabel(AliKalmanTrack *tk, Float_t wrong) const {
6904 //--------------------------------------------------------------------
6905 //This function "cooks" a track label. If label<0, this track is fake.
6906 //--------------------------------------------------------------------
6907 AliTPCseed * t = dynamic_cast<AliTPCseed*>(tk);
6909 printf("%s:%d wrong type \n",(char*)__FILE__,__LINE__);
6913 Int_t noc=t->GetNumberOfClusters();
6915 //printf("\nnot founded prolongation\n\n\n");
6921 AliTPCclusterMI *clusters[160];
6923 for (Int_t i=0;i<160;i++) {
6930 for (i=0; i<160 && current<noc; i++) {
6932 Int_t index=t->GetClusterIndex2(i);
6933 if (index<=0) continue;
6934 if (index&0x8000) continue;
6936 //clusters[current]=GetClusterMI(index);
6937 if (t->GetClusterPointer(i)){
6938 clusters[current]=t->GetClusterPointer(i);
6944 Int_t lab=123456789;
6945 for (i=0; i<noc; i++) {
6946 AliTPCclusterMI *c=clusters[i];
6948 lab=TMath::Abs(c->GetLabel(0));
6950 for (j=0; j<noc; j++) if (lb[j]==lab || mx[j]==0) break;
6956 for (i=0; i<noc; i++) if (mx[i]>max) {max=mx[i]; lab=lb[i];}
6958 for (i=0; i<noc; i++) {
6959 AliTPCclusterMI *c=clusters[i];
6961 if (TMath::Abs(c->GetLabel(1)) == lab ||
6962 TMath::Abs(c->GetLabel(2)) == lab ) max++;
6964 if (noc<=0) { lab=-1; return;}
6965 if ((1.- Float_t(max)/(noc)) > wrong) lab=-lab;
6968 Int_t tail=Int_t(0.10*noc);
6971 for (i=1; i<=160&&ind<tail; i++) {
6972 // AliTPCclusterMI *c=clusters[noc-i];
6973 AliTPCclusterMI *c=clusters[i];
6975 if (lab == TMath::Abs(c->GetLabel(0)) ||
6976 lab == TMath::Abs(c->GetLabel(1)) ||
6977 lab == TMath::Abs(c->GetLabel(2))) max++;
6980 if (max < Int_t(0.5*tail)) lab=-lab;
6987 //delete[] clusters;
6991 //__________________________________________________________________________
6992 Int_t AliTPCtrackerMI::CookLabel(AliTPCseed *const t, Float_t wrong,Int_t first, Int_t last) const {
6993 //--------------------------------------------------------------------
6994 //This function "cooks" a track label. If label<0, this track is fake.
6995 //--------------------------------------------------------------------
6996 Int_t noc=t->GetNumberOfClusters();
6998 //printf("\nnot founded prolongation\n\n\n");
7004 AliTPCclusterMI *clusters[160];
7006 for (Int_t i=0;i<160;i++) {
7013 for (i=0; i<160 && current<noc; i++) {
7014 if (i<first) continue;
7015 if (i>last) continue;
7016 Int_t index=t->GetClusterIndex2(i);
7017 if (index<=0) continue;
7018 if (index&0x8000) continue;
7020 //clusters[current]=GetClusterMI(index);
7021 if (t->GetClusterPointer(i)){
7022 clusters[current]=t->GetClusterPointer(i);
7027 if (noc<5) return -1;
7028 Int_t lab=123456789;
7029 for (i=0; i<noc; i++) {
7030 AliTPCclusterMI *c=clusters[i];
7032 lab=TMath::Abs(c->GetLabel(0));
7034 for (j=0; j<noc; j++) if (lb[j]==lab || mx[j]==0) break;
7040 for (i=0; i<noc; i++) if (mx[i]>max) {max=mx[i]; lab=lb[i];}
7042 for (i=0; i<noc; i++) {
7043 AliTPCclusterMI *c=clusters[i];
7045 if (TMath::Abs(c->GetLabel(1)) == lab ||
7046 TMath::Abs(c->GetLabel(2)) == lab ) max++;
7048 if (noc<=0) { lab=-1; return -1;}
7049 if ((1.- Float_t(max)/(noc)) > wrong) lab=-lab;
7052 Int_t tail=Int_t(0.10*noc);
7055 for (i=1; i<=160&&ind<tail; i++) {
7056 // AliTPCclusterMI *c=clusters[noc-i];
7057 AliTPCclusterMI *c=clusters[i];
7059 if (lab == TMath::Abs(c->GetLabel(0)) ||
7060 lab == TMath::Abs(c->GetLabel(1)) ||
7061 lab == TMath::Abs(c->GetLabel(2))) max++;
7064 if (max < Int_t(0.5*tail)) lab=-lab;
7067 // t->SetLabel(lab);
7071 //delete[] clusters;
7075 Int_t AliTPCtrackerMI::GetRowNumber(Double_t x[3]) const
7077 //return pad row number for given x vector
7078 Float_t phi = TMath::ATan2(x[1],x[0]);
7079 if(phi<0) phi=2.*TMath::Pi()+phi;
7080 // Get the local angle in the sector philoc
7081 const Float_t kRaddeg = 180/3.14159265358979312;
7082 Float_t phiangle = (Int_t (phi*kRaddeg/20.) + 0.5)*20./kRaddeg;
7083 Double_t localx = x[0]*TMath::Cos(phiangle)-x[1]*TMath::Sin(phiangle);
7084 return GetRowNumber(localx);
7089 void AliTPCtrackerMI::MakeBitmaps(AliTPCseed *t)
7091 //-----------------------------------------------------------------------
7092 // Fill the cluster and sharing bitmaps of the track
7093 //-----------------------------------------------------------------------
7095 Int_t firstpoint = 0;
7096 Int_t lastpoint = 159;
7097 AliTPCTrackerPoint *point;
7098 AliTPCclusterMI *cluster;
7100 for (int iter=firstpoint; iter<lastpoint; iter++) {
7101 // Change to cluster pointers to see if we have a cluster at given padrow
7102 cluster = t->GetClusterPointer(iter);
7104 t->SetClusterMapBit(iter, kTRUE);
7105 point = t->GetTrackPoint(iter);
7106 if (point->IsShared())
7107 t->SetSharedMapBit(iter,kTRUE);
7109 t->SetSharedMapBit(iter, kFALSE);
7112 t->SetClusterMapBit(iter, kFALSE);
7113 t->SetSharedMapBit(iter, kFALSE);
7118 Bool_t AliTPCtrackerMI::IsFindable(AliTPCseed & track){
7120 // return flag if there is findable cluster at given position
7123 Float_t z = track.GetZ();
7125 if (TMath::Abs(z)<(AliTPCReconstructor::GetCtgRange()*track.GetX()+kDeltaZ) &&
7126 TMath::Abs(z)<fkParam->GetZLength(0) &&
7127 (TMath::Abs(track.GetSnp())<AliTPCReconstructor::GetMaxSnpTracker()))
7133 void AliTPCtrackerMI::AddCovariance(AliTPCseed * seed){
7135 // Adding systematic error
7136 // !!!! the systematic error for element 4 is in 1/cm not in pt
7138 const Double_t *param = AliTPCReconstructor::GetRecoParam()->GetSystematicError();
7139 Double_t *covarIn= (Double_t*)seed->GetCovariance();
7141 for (Int_t i=0;i<15;i++) covar[i]=0;
7147 covar[0] = param[0]*param[0];
7148 covar[2] = param[1]*param[1];
7149 covar[5] = param[2]*param[2];
7150 covar[9] = param[3]*param[3];
7151 Double_t facC = AliTracker::GetBz()*kB2C;
7152 covar[14]= param[4]*param[4]*facC*facC;
7154 covar[1]=TMath::Sqrt((covar[0]*covar[2]))*covarIn[1]/TMath::Sqrt((covarIn[0]*covarIn[2]));
7156 covar[3]=TMath::Sqrt((covar[0]*covar[5]))*covarIn[3]/TMath::Sqrt((covarIn[0]*covarIn[5]));
7157 covar[4]=TMath::Sqrt((covar[2]*covar[5]))*covarIn[4]/TMath::Sqrt((covarIn[2]*covarIn[5]));
7159 covar[6]=TMath::Sqrt((covar[0]*covar[9]))*covarIn[6]/TMath::Sqrt((covarIn[0]*covarIn[9]));
7160 covar[7]=TMath::Sqrt((covar[2]*covar[9]))*covarIn[7]/TMath::Sqrt((covarIn[2]*covarIn[9]));
7161 covar[8]=TMath::Sqrt((covar[5]*covar[9]))*covarIn[8]/TMath::Sqrt((covarIn[5]*covarIn[9]));
7163 covar[10]=TMath::Sqrt((covar[0]*covar[14]))*covarIn[10]/TMath::Sqrt((covarIn[0]*covarIn[14]));
7164 covar[11]=TMath::Sqrt((covar[2]*covar[14]))*covarIn[11]/TMath::Sqrt((covarIn[2]*covarIn[14]));
7165 covar[12]=TMath::Sqrt((covar[5]*covar[14]))*covarIn[12]/TMath::Sqrt((covarIn[5]*covarIn[14]));
7166 covar[13]=TMath::Sqrt((covar[9]*covar[14]))*covarIn[13]/TMath::Sqrt((covarIn[9]*covarIn[14]));
7168 seed->AddCovariance(covar);