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]);
5263 shared[kink1->GetIndex(0)]= kTRUE;
5264 shared[kink1->GetIndex(1)]= kTRUE;
5265 delete kinks->RemoveAt(indexes[ikink1]);
5272 for (Int_t i=0;i<nkinks;i++){
5273 AliKink * kinkl = (AliKink*) kinks->At(indexes[i]);
5274 if (!kinkl) continue;
5275 kinkl->SetTPCRow0(GetRowNumber(kinkl->GetR()));
5276 Int_t index0 = kinkl->GetIndex(0);
5277 Int_t index1 = kinkl->GetIndex(1);
5278 if (circular[index0]||(circular[index1]&&kinkl->GetDistance()>0.2)) continue;
5279 kinkl->SetMultiple(usage[index0],0);
5280 kinkl->SetMultiple(usage[index1],1);
5281 if (kinkl->GetMultiple()[0]+kinkl->GetMultiple()[1]>2) continue;
5282 if (kinkl->GetMultiple()[0]+kinkl->GetMultiple()[1]>0 && quality[indexes[i]]>0.2) continue;
5283 if (kinkl->GetMultiple()[0]+kinkl->GetMultiple()[1]>0 && kinkl->GetDistance()>0.2) continue;
5284 if (circular[index0]||(circular[index1]&&kinkl->GetDistance()>0.1)) continue;
5286 AliTPCseed * ktrack0 = (AliTPCseed*)array->At(index0);
5287 AliTPCseed * ktrack1 = (AliTPCseed*)array->At(index1);
5288 if (!ktrack0 || !ktrack1) continue;
5289 Int_t index = esd->AddKink(kinkl);
5292 if ( ktrack0->GetKinkIndex(0)==0 && ktrack1->GetKinkIndex(0)==0) { //best kink
5293 if (mothers[indexes[i]].GetNumberOfClusters()>20 && daughters[indexes[i]].GetNumberOfClusters()>20 && (mothers[indexes[i]].GetNumberOfClusters()+daughters[indexes[i]].GetNumberOfClusters())>100){
5294 *ktrack0 = mothers[indexes[i]];
5295 *ktrack1 = daughters[indexes[i]];
5299 ktrack0->SetKinkIndex(usage[index0],-(index+1));
5300 ktrack1->SetKinkIndex(usage[index1], (index+1));
5305 // Remove tracks corresponding to shared kink's
5307 for (Int_t i=0;i<nentries;i++){
5308 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5309 if (!track0) continue;
5310 if (track0->GetKinkIndex(0)!=0) continue;
5311 if (shared[i]) delete array->RemoveAt(i);
5316 RemoveUsed2(array,0.5,0.4,30);
5318 for (Int_t i=0;i<nentries;i++){
5319 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5320 if (!track0) continue;
5321 track0->CookdEdx(0.02,0.6);
5325 for (Int_t i=0;i<nentries;i++){
5326 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5327 if (!track0) continue;
5328 if (track0->Pt()<1.4) continue;
5329 //remove double high momenta tracks - overlapped with kink candidates
5332 for (Int_t icl=track0->GetFirstPoint();icl<track0->GetLastPoint(); icl++){
5333 if (track0->GetClusterPointer(icl)!=0){
5335 if (track0->GetClusterPointer(icl)->IsUsed(10)) ishared++;
5338 if (Float_t(ishared+1)/Float_t(all+1)>0.5) {
5339 delete array->RemoveAt(i);
5343 if (track0->GetKinkIndex(0)!=0) continue;
5344 if (track0->GetNumberOfClusters()<80) continue;
5346 AliTPCseed *pmother = new AliTPCseed();
5347 AliTPCseed *pdaughter = new AliTPCseed();
5348 AliKink *pkink = new AliKink;
5350 AliTPCseed & mother = *pmother;
5351 AliTPCseed & daughter = *pdaughter;
5352 AliKink & kinkl = *pkink;
5353 if (CheckKinkPoint(track0,mother,daughter, kinkl)){
5354 if (mother.GetNumberOfClusters()<30||daughter.GetNumberOfClusters()<20) {
5358 continue; //too short tracks
5360 if (mother.Pt()<1.4) {
5366 Int_t row0= kinkl.GetTPCRow0();
5367 if (kinkl.GetDistance()>0.5 || kinkl.GetR()<110. || kinkl.GetR()>240.) {
5374 Int_t index = esd->AddKink(&kinkl);
5375 mother.SetKinkIndex(0,-(index+1));
5376 daughter.SetKinkIndex(0,index+1);
5377 if (mother.GetNumberOfClusters()>50) {
5378 delete array->RemoveAt(i);
5379 array->AddAt(new AliTPCseed(mother),i);
5382 array->AddLast(new AliTPCseed(mother));
5384 array->AddLast(new AliTPCseed(daughter));
5385 for (Int_t icl=0;icl<row0;icl++) {
5386 if (mother.GetClusterPointer(icl)) mother.GetClusterPointer(icl)->Use(20);
5389 for (Int_t icl=row0;icl<158;icl++) {
5390 if (daughter.GetClusterPointer(icl)) daughter.GetClusterPointer(icl)->Use(20);
5399 delete [] daughters;
5421 AliInfo(Form("Ncandidates=\t%d\t%d\t%d\t%d\n",esd->GetNumberOfKinks(),ncandidates,ntracks,nall));
5425 void AliTPCtrackerMI::FindV0s(const TObjArray * array, AliESDEvent *const esd)
5431 TObjArray *tpcv0s = new TObjArray(100000);
5432 Int_t nentries = array->GetEntriesFast();
5433 AliHelix *helixes = new AliHelix[nentries];
5434 Int_t *sign = new Int_t[nentries];
5435 Float_t *alpha = new Float_t[nentries];
5436 Float_t *z0 = new Float_t[nentries];
5437 Float_t *dca = new Float_t[nentries];
5438 Float_t *sdcar = new Float_t[nentries];
5439 Float_t *cdcar = new Float_t[nentries];
5440 Float_t *pulldcar = new Float_t[nentries];
5441 Float_t *pulldcaz = new Float_t[nentries];
5442 Float_t *pulldca = new Float_t[nentries];
5443 Bool_t *isPrim = new Bool_t[nentries];
5444 const AliESDVertex * primvertex = esd->GetVertex();
5445 Double_t zvertex = primvertex->GetZv();
5447 // nentries = array->GetEntriesFast();
5449 for (Int_t i=0;i<nentries;i++){
5452 AliTPCseed* track = (AliTPCseed*)array->At(i);
5453 if (!track) continue;
5454 track->GetV0Indexes()[0] = 0; //rest v0 indexes
5455 track->GetV0Indexes()[1] = 0; //rest v0 indexes
5456 track->GetV0Indexes()[2] = 0; //rest v0 indexes
5458 alpha[i] = track->GetAlpha();
5459 new (&helixes[i]) AliHelix(*track);
5461 helixes[i].Evaluate(0,xyz);
5462 sign[i] = (track->GetC()>0) ? -1:1;
5466 if (track->GetProlongation(0,y,z)) z0[i] = z;
5467 dca[i] = track->GetD(0,0);
5469 // dca error parrameterezation + pulls
5471 sdcar[i] = TMath::Sqrt(0.150*0.150+(100*track->GetC())*(100*track->GetC()));
5472 if (TMath::Abs(track->GetTgl())>1) sdcar[i]*=2.5;
5473 cdcar[i] = TMath::Exp((TMath::Abs(track->GetC())-0.0106)*525.3);
5474 pulldcar[i] = (dca[i]-cdcar[i])/sdcar[i];
5475 pulldcaz[i] = (z0[i]-zvertex)/sdcar[i];
5476 pulldca[i] = TMath::Sqrt(pulldcar[i]*pulldcar[i]+pulldcaz[i]*pulldcaz[i]);
5477 if (track->TPCrPID(1)+track->TPCrPID(2)+track->TPCrPID(3)>0.5) {
5478 if (pulldca[i]<3.) isPrim[i]=kTRUE; //pion, muon and Kaon 3 sigma cut
5480 if (track->TPCrPID(4)>0.5) {
5481 if (pulldca[i]<0.5) isPrim[i]=kTRUE; //proton 0.5 sigma cut
5483 if (track->TPCrPID(0)>0.4) {
5484 isPrim[i]=kFALSE; //electron no sigma cut
5491 Int_t ncandidates =0;
5494 Double_t phase[2][2],radius[2];
5500 Float_t fprimvertex[3]={GetX(),GetY(),GetZ()};
5502 Double_t cradius0 = 10*10;
5503 Double_t cradius1 = 200*200;
5506 Double_t cpointAngle = 0.95;
5508 Double_t delta[2]={10000,10000};
5509 for (Int_t i =0;i<nentries;i++){
5510 if (sign[i]==0) continue;
5511 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5512 if (!track0) continue;
5513 if (AliTPCReconstructor::StreamLevel()>1){
5514 TTreeSRedirector &cstream = *fDebugStreamer;
5519 "zvertex="<<zvertex<<
5520 "sdcar0="<<sdcar[i]<<
5521 "cdcar0="<<cdcar[i]<<
5522 "pulldcar0="<<pulldcar[i]<<
5523 "pulldcaz0="<<pulldcaz[i]<<
5524 "pulldca0="<<pulldca[i]<<
5525 "isPrim="<<isPrim[i]<<
5529 if (track0->GetSigned1Pt()<0) continue;
5530 if (track0->GetKinkIndex(0)>0||isPrim[i]) continue; //daughter kink
5532 if (TMath::Abs(helixes[i].GetHelix(4))<0.000000001) continue;
5537 for (Int_t j =0;j<nentries;j++){
5538 AliTPCseed * track1 = (AliTPCseed*)array->At(j);
5539 if (!track1) continue;
5540 if (track1->GetKinkIndex(0)>0 || isPrim[j]) continue; //daughter kink
5541 if (sign[j]*sign[i]>0) continue;
5542 if (TMath::Abs(helixes[j].GetHelix(4))<0.000001) continue;
5543 if (track0->GetCircular()+track1->GetCircular()>1) continue; //circling -returning track
5546 // DCA to prim vertex cut
5552 Int_t npoints = helixes[i].GetRPHIintersections(helixes[j], phase, radius,cdist2);
5553 if (npoints<1) continue;
5557 if (radius[0]<cradius0||radius[0]>cradius1) continue;
5558 helixes[i].LinearDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta[0]);
5559 if (delta[0]>cdist1) continue;
5562 if (TMath::Max(radius[0],radius[1])<cradius0|| TMath::Min(radius[0],radius[1])>cradius1) continue;
5563 helixes[i].LinearDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta[0]);
5564 helixes[i].LinearDCA(helixes[j],phase[1][0],phase[1][1],radius[1],delta[1]);
5565 if (delta[1]<delta[0]) iclosest=1;
5566 if (delta[iclosest]>cdist1) continue;
5568 helixes[i].ParabolicDCA(helixes[j],phase[iclosest][0],phase[iclosest][1],radius[iclosest],delta[iclosest]);
5569 if (radius[iclosest]<cradius0 || radius[iclosest]>cradius1 || delta[iclosest]>cdist1) continue;
5571 Double_t pointAngle = helixes[i].GetPointAngle(helixes[j],phase[iclosest],fprimvertex);
5572 if (pointAngle<cpointAngle) continue;
5574 Bool_t isGamma = kFALSE;
5575 vertex.SetParamP(*track0); //track0 - plus
5576 vertex.SetParamN(*track1); //track1 - minus
5577 vertex.Update(fprimvertex);
5578 if (track0->TPCrPID(0)>0.3&&track1->TPCrPID(0)>0.3&&vertex.GetAnglep()[2]<0.15) isGamma=kTRUE; // gamma conversion candidate
5579 Double_t pointAngle2 = vertex.GetV0CosineOfPointingAngle();
5581 if (vertex.GetV0CosineOfPointingAngle()<cpointAngle && (!isGamma)) continue;// point angle cut
5582 //Bo: if (vertex.GetDist2()>2&&(!isGamma)) continue; // point angle cut
5583 if (vertex.GetDcaV0Daughters()>2&&(!isGamma)) continue;//Bo: // point angle cut
5584 Float_t sigmae = 0.15*0.15;
5585 if (vertex.GetRr()<80)
5586 sigmae += (sdcar[i]*sdcar[i]+sdcar[j]*sdcar[j])*(1.-vertex.GetRr()/80.)*(1.-vertex.GetRr()/80.);
5587 sigmae+= TMath::Sqrt(sigmae);
5588 //Bo: if (vertex.GetDist2()/sigmae>3.&&(!isGamma)) continue;
5589 if (vertex.GetDcaV0Daughters()/sigmae>3.&&(!isGamma)) continue;
5590 Float_t densb0=0,densb1=0,densa0=0,densa1=0;
5591 Int_t row0 = GetRowNumber(vertex.GetRr());
5593 //Bo: if (vertex.GetDist2()>0.2) continue;
5594 if (vertex.GetDcaV0Daughters()>0.2) continue;
5595 densb0 = track0->Density2(0,row0-5);
5596 densb1 = track1->Density2(0,row0-5);
5597 if (densb0>0.3|| densb1>0.3) continue; //clusters before vertex
5598 densa0 = track0->Density2(row0+5,row0+40);
5599 densa1 = track1->Density2(row0+5,row0+40);
5600 if ((densa0<0.4|| densa1<0.4)&&(!isGamma)) continue; //missing clusters after vertex
5603 densa0 = track0->Density2(0,40); //cluster density
5604 densa1 = track1->Density2(0,40); //cluster density
5605 if ((vertex.GetRr()<80&&densa0+densa1<1.)&&(!isGamma)) continue;
5607 //Bo: vertex.SetLab(0,track0->GetLabel());
5608 //Bo: vertex.SetLab(1,track1->GetLabel());
5609 vertex.SetChi2After((densa0+densa1)*0.5);
5610 vertex.SetChi2Before((densb0+densb1)*0.5);
5611 vertex.SetIndex(0,i);
5612 vertex.SetIndex(1,j);
5613 //Bo: vertex.SetStatus(1); // TPC v0 candidate
5614 vertex.SetOnFlyStatus(2);//Bo: // TPC v0 candidate
5615 //Bo: vertex.SetRp(track0->TPCrPIDs());
5616 //Bo: vertex.SetRm(track1->TPCrPIDs());
5617 tpcv0s->AddLast(new AliESDv0(vertex));
5620 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
5621 Double_t radiusm= (delta[0]<delta[1])? TMath::Sqrt(radius[0]):TMath::Sqrt(radius[1]);
5622 Double_t deltam= (delta[0]<delta[1])? TMath::Sqrt(delta[0]):TMath::Sqrt(delta[1]);
5623 if (AliTPCReconstructor::StreamLevel()>1) {
5624 Int_t lab0=track0->GetLabel();
5625 Int_t lab1=track1->GetLabel();
5626 Char_t c0=track0->GetCircular();
5627 Char_t c1=track1->GetCircular();
5628 TTreeSRedirector &cstream = *fDebugStreamer;
5631 "vertex.="<<&vertex<<
5634 "Helix0.="<<&helixes[i]<<
5637 "Helix1.="<<&helixes[j]<<
5638 "pointAngle="<<pointAngle<<
5639 "pointAngle2="<<pointAngle2<<
5644 "zvertex="<<zvertex<<
5647 "npoints="<<npoints<<
5648 "radius0="<<radius[0]<<
5649 "delta0="<<delta[0]<<
5650 "radius1="<<radius[1]<<
5651 "delta1="<<delta[1]<<
5652 "radiusm="<<radiusm<<
5654 "sdcar0="<<sdcar[i]<<
5655 "sdcar1="<<sdcar[j]<<
5656 "cdcar0="<<cdcar[i]<<
5657 "cdcar1="<<cdcar[j]<<
5658 "pulldcar0="<<pulldcar[i]<<
5659 "pulldcar1="<<pulldcar[j]<<
5660 "pulldcaz0="<<pulldcaz[i]<<
5661 "pulldcaz1="<<pulldcaz[j]<<
5662 "pulldca0="<<pulldca[i]<<
5663 "pulldca1="<<pulldca[j]<<
5673 Float_t *quality = new Float_t[ncandidates];
5674 Int_t *indexes = new Int_t[ncandidates];
5676 for (Int_t i=0;i<ncandidates;i++){
5678 AliESDv0 *v0 = (AliESDv0*)tpcv0s->At(i);
5679 quality[i] = 1./(1.00001-v0->GetV0CosineOfPointingAngle()); //base point angle
5680 // quality[i] /= (0.5+v0->GetDist2());
5681 // quality[i] *= v0->GetChi2After(); //density factor
5683 Int_t index0 = v0->GetIndex(0);
5684 Int_t index1 = v0->GetIndex(1);
5685 //Bo: Double_t minpulldca = TMath::Min(2.+pulldca[v0->GetIndex(0)],(2.+pulldca[v0->GetIndex(1)]) ); //pull
5686 Double_t minpulldca = TMath::Min(2.+pulldca[index0],(2.+pulldca[index1]) );//Bo:
5690 AliTPCseed * track0 = (AliTPCseed*)array->At(index0);
5691 AliTPCseed * track1 = (AliTPCseed*)array->At(index1);
5692 if (track0->TPCrPID(0)>0.3&&track1->TPCrPID(0)>0.3&&v0->GetAnglep()[2]<0.15) quality[i]+=1000000; // gamma conversion candidate
5693 if (track0->TPCrPID(4)>0.9||(track1->TPCrPID(4)>0.9&&minpulldca>4)) quality[i]*=10; // lambda candidate candidate
5696 TMath::Sort(ncandidates,quality,indexes,kTRUE);
5699 for (Int_t i=0;i<ncandidates;i++){
5700 AliESDv0 * v0 = (AliESDv0*)tpcv0s->At(indexes[i]);
5702 Int_t index0 = v0->GetIndex(0);
5703 Int_t index1 = v0->GetIndex(1);
5704 AliTPCseed * track0 = (AliTPCseed*)array->At(index0);
5705 AliTPCseed * track1 = (AliTPCseed*)array->At(index1);
5706 if (!track0||!track1) {
5710 Bool_t accept =kTRUE; //default accept
5711 Int_t *v0indexes0 = track0->GetV0Indexes();
5712 Int_t *v0indexes1 = track1->GetV0Indexes();
5714 Int_t order0 = (v0indexes0[0]!=0) ? 1:0;
5715 Int_t order1 = (v0indexes1[0]!=0) ? 1:0;
5716 if (v0indexes0[1]!=0) order0 =2;
5717 if (v0indexes1[1]!=0) order1 =2;
5719 if (v0indexes0[2]!=0) {order0=3; accept=kFALSE;}
5720 if (v0indexes0[2]!=0) {order1=3; accept=kFALSE;}
5722 AliESDv0 * v02 = v0;
5724 //Bo: v0->SetOrder(0,order0);
5725 //Bo: v0->SetOrder(1,order1);
5726 //Bo: v0->SetOrder(1,order0+order1);
5727 v0->SetOnFlyStatus(kTRUE);
5728 Int_t index = esd->AddV0(v0);
5729 v02 = esd->GetV0(index);
5730 v0indexes0[order0]=index;
5731 v0indexes1[order1]=index;
5735 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
5736 if (AliTPCReconstructor::StreamLevel()>1) {
5737 Int_t lab0=track0->GetLabel();
5738 Int_t lab1=track1->GetLabel();
5739 TTreeSRedirector &cstream = *fDebugStreamer;
5748 "dca0="<<dca[index0]<<
5749 "dca1="<<dca[index1]<<
5753 "quality="<<quality[i]<<
5754 "pulldca0="<<pulldca[index0]<<
5755 "pulldca1="<<pulldca[index1]<<
5779 printf("TPC V0 finder : naccepted\t%d\tncandidates\t%d\tntracks\t%d\tnall\t%d\n",naccepted,ncandidates,ntracks,nall);
5783 Int_t AliTPCtrackerMI::RefitKink(AliTPCseed &mother, AliTPCseed &daughter, const AliESDkink &knk)
5786 // refit kink towards to the vertex
5789 AliKink &kink=(AliKink &)knk;
5791 Int_t row0 = GetRowNumber(kink.GetR());
5792 FollowProlongation(mother,0);
5793 mother.Reset(kFALSE);
5795 FollowProlongation(daughter,row0);
5796 daughter.Reset(kFALSE);
5797 FollowBackProlongation(daughter,158);
5798 daughter.Reset(kFALSE);
5799 Int_t first = TMath::Max(row0-20,30);
5800 Int_t last = TMath::Min(row0+20,140);
5802 const Int_t kNdiv =5;
5803 AliTPCseed param0[kNdiv]; // parameters along the track
5804 AliTPCseed param1[kNdiv]; // parameters along the track
5805 AliKink kinks[kNdiv]; // corresponding kink parameters
5808 for (Int_t irow=0; irow<kNdiv;irow++){
5809 rows[irow] = first +((last-first)*irow)/(kNdiv-1);
5811 // store parameters along the track
5813 for (Int_t irow=0;irow<kNdiv;irow++){
5814 FollowBackProlongation(mother, rows[irow]);
5815 FollowProlongation(daughter,rows[kNdiv-1-irow]);
5816 param0[irow] = mother;
5817 param1[kNdiv-1-irow] = daughter;
5821 for (Int_t irow=0; irow<kNdiv-1;irow++){
5822 if (param0[irow].GetNumberOfClusters()<kNdiv||param1[irow].GetNumberOfClusters()<kNdiv) continue;
5823 kinks[irow].SetMother(param0[irow]);
5824 kinks[irow].SetDaughter(param1[irow]);
5825 kinks[irow].Update();
5828 // choose kink with best "quality"
5830 Double_t mindist = 10000;
5831 for (Int_t irow=0;irow<kNdiv;irow++){
5832 if (param0[irow].GetNumberOfClusters()<20||param1[irow].GetNumberOfClusters()<20) continue;
5833 if (TMath::Abs(kinks[irow].GetR())>240.) continue;
5834 if (TMath::Abs(kinks[irow].GetR())<100.) continue;
5836 Float_t normdist = TMath::Abs(param0[irow].GetX()-kinks[irow].GetR())*(0.1+kink.GetDistance());
5837 normdist/= (param0[irow].GetNumberOfClusters()+param1[irow].GetNumberOfClusters()+40.);
5838 if (normdist < mindist){
5844 if (index==-1) return 0;
5847 param0[index].Reset(kTRUE);
5848 FollowProlongation(param0[index],0);
5850 mother = param0[index];
5851 daughter = param1[index]; // daughter in vertex
5853 kink.SetMother(mother);
5854 kink.SetDaughter(daughter);
5856 kink.SetTPCRow0(GetRowNumber(kink.GetR()));
5857 kink.SetTPCncls(param0[index].GetNumberOfClusters(),0);
5858 kink.SetTPCncls(param1[index].GetNumberOfClusters(),1);
5859 kink.SetLabel(CookLabel(&mother,0.4, 0,kink.GetTPCRow0()),0);
5860 kink.SetLabel(CookLabel(&daughter,0.4, kink.GetTPCRow0(),160),1);
5861 mother.SetLabel(kink.GetLabel(0));
5862 daughter.SetLabel(kink.GetLabel(1));
5868 void AliTPCtrackerMI::UpdateKinkQualityM(AliTPCseed * seed){
5870 // update Kink quality information for mother after back propagation
5872 if (seed->GetKinkIndex(0)>=0) return;
5873 for (Int_t ikink=0;ikink<3;ikink++){
5874 Int_t index = seed->GetKinkIndex(ikink);
5875 if (index>=0) break;
5876 index = TMath::Abs(index)-1;
5877 AliESDkink * kink = fEvent->GetKink(index);
5878 kink->SetTPCDensity(-1,0,0);
5879 kink->SetTPCDensity(1,0,1);
5881 Int_t row0 = kink->GetTPCRow0() - 2 - Int_t( 0.5/ (0.05+kink->GetAngle(2)));
5882 if (row0<15) row0=15;
5884 Int_t row1 = kink->GetTPCRow0() + 2 + Int_t( 0.5/ (0.05+kink->GetAngle(2)));
5885 if (row1>145) row1=145;
5887 Int_t found,foundable,shared;
5888 seed->GetClusterStatistic(0,row0, found, foundable,shared,kFALSE);
5889 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),0,0);
5890 seed->GetClusterStatistic(row1,155, found, foundable,shared,kFALSE);
5891 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),0,1);
5896 void AliTPCtrackerMI::UpdateKinkQualityD(AliTPCseed * seed){
5898 // update Kink quality information for daughter after refit
5900 if (seed->GetKinkIndex(0)<=0) return;
5901 for (Int_t ikink=0;ikink<3;ikink++){
5902 Int_t index = seed->GetKinkIndex(ikink);
5903 if (index<=0) break;
5904 index = TMath::Abs(index)-1;
5905 AliESDkink * kink = fEvent->GetKink(index);
5906 kink->SetTPCDensity(-1,1,0);
5907 kink->SetTPCDensity(-1,1,1);
5909 Int_t row0 = kink->GetTPCRow0() -2 - Int_t( 0.5/ (0.05+kink->GetAngle(2)));
5910 if (row0<15) row0=15;
5912 Int_t row1 = kink->GetTPCRow0() +2 + Int_t( 0.5/ (0.05+kink->GetAngle(2)));
5913 if (row1>145) row1=145;
5915 Int_t found,foundable,shared;
5916 seed->GetClusterStatistic(0,row0, found, foundable,shared,kFALSE);
5917 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),1,0);
5918 seed->GetClusterStatistic(row1,155, found, foundable,shared,kFALSE);
5919 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),1,1);
5925 Int_t AliTPCtrackerMI::CheckKinkPoint(AliTPCseed*seed,AliTPCseed &mother, AliTPCseed &daughter, const AliESDkink &knk)
5928 // check kink point for given track
5929 // if return value=0 kink point not found
5930 // otherwise seed0 correspond to mother particle
5931 // seed1 correspond to daughter particle
5932 // kink parameter of kink point
5933 AliKink &kink=(AliKink &)knk;
5935 Int_t middlerow = (seed->GetFirstPoint()+seed->GetLastPoint())/2;
5936 Int_t first = seed->GetFirstPoint();
5937 Int_t last = seed->GetLastPoint();
5938 if (last-first<20) return 0; // shortest length - 2*30 = 60 pad-rows
5941 AliTPCseed *seed1 = ReSeed(seed,middlerow+20, kTRUE); //middle of chamber
5942 if (!seed1) return 0;
5943 FollowProlongation(*seed1,seed->GetLastPoint()-20);
5944 seed1->Reset(kTRUE);
5945 FollowProlongation(*seed1,158);
5946 seed1->Reset(kTRUE);
5947 last = seed1->GetLastPoint();
5949 AliTPCseed *seed0 = new AliTPCseed(*seed);
5950 seed0->Reset(kFALSE);
5953 AliTPCseed param0[20]; // parameters along the track
5954 AliTPCseed param1[20]; // parameters along the track
5955 AliKink kinks[20]; // corresponding kink parameters
5957 for (Int_t irow=0; irow<20;irow++){
5958 rows[irow] = first +((last-first)*irow)/19;
5960 // store parameters along the track
5962 for (Int_t irow=0;irow<20;irow++){
5963 FollowBackProlongation(*seed0, rows[irow]);
5964 FollowProlongation(*seed1,rows[19-irow]);
5965 param0[irow] = *seed0;
5966 param1[19-irow] = *seed1;
5970 for (Int_t irow=0; irow<19;irow++){
5971 kinks[irow].SetMother(param0[irow]);
5972 kinks[irow].SetDaughter(param1[irow]);
5973 kinks[irow].Update();
5976 // choose kink with biggest change of angle
5978 Double_t maxchange= 0;
5979 for (Int_t irow=1;irow<19;irow++){
5980 if (TMath::Abs(kinks[irow].GetR())>240.) continue;
5981 if (TMath::Abs(kinks[irow].GetR())<110.) continue;
5982 Float_t quality = TMath::Abs(kinks[irow].GetAngle(2))/(3.+TMath::Abs(kinks[irow].GetR()-param0[irow].GetX()));
5983 if ( quality > maxchange){
5984 maxchange = quality;
5991 if (index<0) return 0;
5993 Int_t row0 = GetRowNumber(kinks[index].GetR()); //row 0 estimate
5994 seed0 = new AliTPCseed(param0[index]);
5995 seed1 = new AliTPCseed(param1[index]);
5996 seed0->Reset(kFALSE);
5997 seed1->Reset(kFALSE);
5998 seed0->ResetCovariance(10.);
5999 seed1->ResetCovariance(10.);
6000 FollowProlongation(*seed0,0);
6001 FollowBackProlongation(*seed1,158);
6002 mother = *seed0; // backup mother at position 0
6003 seed0->Reset(kFALSE);
6004 seed1->Reset(kFALSE);
6005 seed0->ResetCovariance(10.);
6006 seed1->ResetCovariance(10.);
6008 first = TMath::Max(row0-20,0);
6009 last = TMath::Min(row0+20,158);
6011 for (Int_t irow=0; irow<20;irow++){
6012 rows[irow] = first +((last-first)*irow)/19;
6014 // store parameters along the track
6016 for (Int_t irow=0;irow<20;irow++){
6017 FollowBackProlongation(*seed0, rows[irow]);
6018 FollowProlongation(*seed1,rows[19-irow]);
6019 param0[irow] = *seed0;
6020 param1[19-irow] = *seed1;
6024 for (Int_t irow=0; irow<19;irow++){
6025 kinks[irow].SetMother(param0[irow]);
6026 kinks[irow].SetDaughter(param1[irow]);
6027 // param0[irow].Dump();
6028 //param1[irow].Dump();
6029 kinks[irow].Update();
6032 // choose kink with biggest change of angle
6035 for (Int_t irow=0;irow<20;irow++){
6036 if (TMath::Abs(kinks[irow].GetR())>250.) continue;
6037 if (TMath::Abs(kinks[irow].GetR())<90.) continue;
6038 Float_t quality = TMath::Abs(kinks[irow].GetAngle(2))/(3.+TMath::Abs(kinks[irow].GetR()-param0[irow].GetX()));
6039 if ( quality > maxchange){
6040 maxchange = quality;
6047 if (index==-1 || param0[index].GetNumberOfClusters()+param1[index].GetNumberOfClusters()<100){
6053 // Float_t anglesigma = TMath::Sqrt(param0[index].fC22+param0[index].fC33+param1[index].fC22+param1[index].fC33);
6055 kink.SetMother(param0[index]);
6056 kink.SetDaughter(param1[index]);
6059 Double_t chi2P2 = param0[index].GetParameter()[2]-param1[index].GetParameter()[2];
6061 chi2P2/=param0[index].GetCovariance()[5]+param1[index].GetCovariance()[5];
6062 Double_t chi2P3 = param0[index].GetParameter()[3]-param1[index].GetParameter()[3];
6064 chi2P3/=param0[index].GetCovariance()[9]+param1[index].GetCovariance()[9];
6066 if (AliTPCReconstructor::StreamLevel()>1) {
6067 (*fDebugStreamer)<<"kinkHpt"<<
6070 "p0.="<<¶m0[index]<<
6071 "p1.="<<¶m1[index]<<
6075 if ( chi2P2+chi2P3<AliTPCReconstructor::GetRecoParam()->GetKinkAngleCutChi2(0)){
6082 row0 = GetRowNumber(kink.GetR());
6083 kink.SetTPCRow0(row0);
6084 kink.SetLabel(CookLabel(seed0,0.5,0,row0),0);
6085 kink.SetLabel(CookLabel(seed1,0.5,row0,158),1);
6086 kink.SetIndex(-10,0);
6087 kink.SetIndex(int(param0[index].GetNumberOfClusters()+param1[index].GetNumberOfClusters()),1);
6088 kink.SetTPCncls(param0[index].GetNumberOfClusters(),0);
6089 kink.SetTPCncls(param1[index].GetNumberOfClusters(),1);
6092 // new (&mother) AliTPCseed(param0[index]);
6093 daughter = param1[index];
6094 daughter.SetLabel(kink.GetLabel(1));
6095 param0[index].Reset(kTRUE);
6096 FollowProlongation(param0[index],0);
6097 mother = param0[index];
6098 mother.SetLabel(kink.GetLabel(0));
6099 if ( chi2P2+chi2P3<AliTPCReconstructor::GetRecoParam()->GetKinkAngleCutChi2(1)){
6111 AliTPCseed* AliTPCtrackerMI::ReSeed(AliTPCseed *t)
6114 // reseed - refit - track
6117 // Int_t last = fSectors->GetNRows()-1;
6119 if (fSectors == fOuterSec){
6120 first = TMath::Max(first, t->GetFirstPoint()-fInnerSec->GetNRows());
6124 first = t->GetFirstPoint();
6126 AliTPCseed * seed = MakeSeed(t,0.1,0.5,0.9);
6127 FollowBackProlongation(*t,fSectors->GetNRows()-1);
6129 FollowProlongation(*t,first);
6139 //_____________________________________________________________________________
6140 Int_t AliTPCtrackerMI::ReadSeeds(const TFile *inp) {
6141 //-----------------------------------------------------------------
6142 // This function reades track seeds.
6143 //-----------------------------------------------------------------
6144 TDirectory *savedir=gDirectory;
6146 TFile *in=(TFile*)inp;
6147 if (!in->IsOpen()) {
6148 cerr<<"AliTPCtrackerMI::ReadSeeds(): input file is not open !\n";
6153 TTree *seedTree=(TTree*)in->Get("Seeds");
6155 cerr<<"AliTPCtrackerMI::ReadSeeds(): ";
6156 cerr<<"can't get a tree with track seeds !\n";
6159 AliTPCtrack *seed=new AliTPCtrack;
6160 seedTree->SetBranchAddress("tracks",&seed);
6162 if (fSeeds==0) fSeeds=new TObjArray(15000);
6164 Int_t n=(Int_t)seedTree->GetEntries();
6165 for (Int_t i=0; i<n; i++) {
6166 seedTree->GetEvent(i);
6167 fSeeds->AddLast(new AliTPCseed(*seed/*,seed->GetAlpha()*/));
6176 Int_t AliTPCtrackerMI::Clusters2Tracks (AliESDEvent *const esd)
6179 if (fSeeds) DeleteSeeds();
6182 if (!fSeeds) return 1;
6184 if (AliTPCReconstructor::StreamLevel()>3) DumpClusters(0,fSeeds);
6190 //_____________________________________________________________________________
6191 Int_t AliTPCtrackerMI::Clusters2Tracks() {
6192 //-----------------------------------------------------------------
6193 // This is a track finder.
6194 //-----------------------------------------------------------------
6195 TDirectory *savedir=gDirectory;
6199 fSeeds = Tracking();
6202 Info("Clusters2Tracks","Time for tracking: \t");timer.Print();timer.Start();
6204 //activate again some tracks
6205 for (Int_t i=0; i<fSeeds->GetEntriesFast(); i++) {
6206 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
6208 Int_t nc=t.GetNumberOfClusters();
6210 delete fSeeds->RemoveAt(i);
6214 if (pt->GetRemoval()==10) {
6215 if (pt->GetDensityFirst(20)>0.8 || pt->GetDensityFirst(30)>0.8 || pt->GetDensityFirst(40)>0.7)
6216 pt->Desactivate(10); // make track again active
6218 pt->Desactivate(20);
6219 delete fSeeds->RemoveAt(i);
6224 RemoveUsed2(fSeeds,0.85,0.85,0);
6225 if (AliTPCReconstructor::GetRecoParam()->GetDoKinks()) FindKinks(fSeeds,fEvent);
6226 //FindCurling(fSeeds, fEvent,0);
6227 if (AliTPCReconstructor::StreamLevel()>5) FindMultiMC(fSeeds, fEvent,-1); // find multi found tracks
6228 RemoveUsed2(fSeeds,0.5,0.4,20);
6229 FindSplitted(fSeeds, fEvent,0); // find multi found tracks
6230 if (AliTPCReconstructor::StreamLevel()>5) FindMultiMC(fSeeds, fEvent,0); // find multi found tracks
6233 // // refit short tracks
6235 Int_t nseed=fSeeds->GetEntriesFast();
6238 for (Int_t i=0; i<nseed; i++) {
6239 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
6241 Int_t nc=t.GetNumberOfClusters();
6243 delete fSeeds->RemoveAt(i);
6246 CookLabel(pt,0.1); //For comparison only
6247 //if ((pt->IsActive() || (pt->fRemoval==10) )&& nc>50 &&pt->GetNumberOfClusters()>0.4*pt->fNFoundable){
6248 if ((pt->IsActive() || (pt->GetRemoval()==10) )){
6250 if (fDebug>0) cerr<<found<<'\r';
6254 delete fSeeds->RemoveAt(i);
6258 //RemoveOverlap(fSeeds,0.99,7,kTRUE);
6260 //RemoveUsed(fSeeds,0.9,0.9,6);
6262 nseed=fSeeds->GetEntriesFast();
6264 for (Int_t i=0; i<nseed; i++) {
6265 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
6267 Int_t nc=t.GetNumberOfClusters();
6269 delete fSeeds->RemoveAt(i);
6273 t.CookdEdx(0.02,0.6);
6274 // CheckKinkPoint(&t,0.05);
6275 //if ((pt->IsActive() || (pt->fRemoval==10) )&& nc>50 &&pt->GetNumberOfClusters()>0.4*pt->fNFoundable){
6276 if ((pt->IsActive() || (pt->GetRemoval()==10) )){
6284 delete fSeeds->RemoveAt(i);
6285 //AliTPCseed * seed1 = ReSeed(pt,0.05,0.5,1);
6287 // FollowProlongation(*seed1,0);
6288 // Int_t n = seed1->GetNumberOfClusters();
6289 // printf("fP4\t%f\t%f\n",seed1->GetC(),pt->GetC());
6290 // printf("fN\t%d\t%d\n", seed1->GetNumberOfClusters(),pt->GetNumberOfClusters());
6293 //AliTPCseed * seed2 = ReSeed(pt,0.95,0.5,0.05);
6297 SortTracks(fSeeds, 1);
6301 PrepareForBackProlongation(fSeeds,5.);
6302 PropagateBack(fSeeds);
6303 printf("Time for back propagation: \t");timer.Print();timer.Start();
6307 PrepareForProlongation(fSeeds,5.);
6308 PropagateForward2(fSeeds);
6310 printf("Time for FORWARD propagation: \t");timer.Print();timer.Start();
6311 // RemoveUsed(fSeeds,0.7,0.7,6);
6312 //RemoveOverlap(fSeeds,0.9,7,kTRUE);
6314 nseed=fSeeds->GetEntriesFast();
6316 for (Int_t i=0; i<nseed; i++) {
6317 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
6319 Int_t nc=t.GetNumberOfClusters();
6321 delete fSeeds->RemoveAt(i);
6324 t.CookdEdx(0.02,0.6);
6325 // CookLabel(pt,0.1); //For comparison only
6326 //if ((pt->IsActive() || (pt->fRemoval==10) )&& nc>50 &&pt->GetNumberOfClusters()>0.4*pt->fNFoundable){
6327 if ((pt->IsActive() || (pt->fRemoval==10) )){
6328 cerr<<found++<<'\r';
6331 delete fSeeds->RemoveAt(i);
6336 // fNTracks = found;
6338 Info("Clusters2Tracks","Time for overlap removal, track writing and dedx cooking: \t"); timer.Print();timer.Start();
6341 // cerr<<"Number of found tracks : "<<"\t"<<found<<endl;
6342 Info("Clusters2Tracks","Number of found tracks %d",found);
6344 // UnloadClusters();
6349 void AliTPCtrackerMI::Tracking(TObjArray * arr)
6352 // tracking of the seeds
6355 fSectors = fOuterSec;
6356 ParallelTracking(arr,150,63);
6357 fSectors = fOuterSec;
6358 ParallelTracking(arr,63,0);
6361 TObjArray * AliTPCtrackerMI::Tracking(Int_t seedtype, Int_t i1, Int_t i2, Float_t cuts[4], Float_t dy, Int_t dsec)
6366 TObjArray * arr = new TObjArray;
6368 fSectors = fOuterSec;
6371 for (Int_t sec=0;sec<fkNOS;sec++){
6372 if (seedtype==3) MakeSeeds3(arr,sec,i1,i2,cuts,dy, dsec);
6373 if (seedtype==4) MakeSeeds5(arr,sec,i1,i2,cuts,dy);
6374 if (seedtype==2) MakeSeeds2(arr,sec,i1,i2,cuts,dy);
6377 Info("Tracking","\nSeeding - %d\t%d\t%d\t%d\n",seedtype,i1,i2,arr->GetEntriesFast());
6389 TObjArray * AliTPCtrackerMI::Tracking()
6393 if (AliTPCReconstructor::GetRecoParam()->GetSpecialSeeding()) return TrackingSpecial();
6396 Int_t nup=fOuterSec->GetNRows()+fInnerSec->GetNRows();
6398 TObjArray * seeds = new TObjArray;
6400 Int_t fLastSeedRowSec=AliTPCReconstructor::GetRecoParam()->GetLastSeedRowSec();
6401 Int_t gapPrim = AliTPCReconstructor::GetRecoParam()->GetSeedGapPrim();
6402 Int_t gapSec = AliTPCReconstructor::GetRecoParam()->GetSeedGapSec();
6410 Float_t fnumber = 3.0;
6411 Float_t fdensity = 3.0;
6416 for (Int_t delta = 0; delta<18; delta+=gapPrim){
6420 arr = Tracking(3,nup-1-delta,nup-1-delta-gap,cuts,-1,1);
6421 SumTracks(seeds,arr);
6422 SignClusters(seeds,fnumber,fdensity);
6424 for (Int_t i=2;i<6;i+=2){
6425 // seed high pt tracks
6428 arr = Tracking(3,nup-i-delta,nup-i-delta-gap,cuts,-1,0);
6429 SumTracks(seeds,arr);
6430 SignClusters(seeds,fnumber,fdensity);
6435 // RemoveUsed(seeds,0.9,0.9,1);
6436 // UnsignClusters();
6437 // SignClusters(seeds,fnumber,fdensity);
6441 for (Int_t delta = 20; delta<120; delta+=gapPrim){
6443 // seed high pt tracks
6447 arr = Tracking(3,nup-delta,nup-delta-gap,cuts,-1);
6448 SumTracks(seeds,arr);
6449 SignClusters(seeds,fnumber,fdensity);
6454 arr = Tracking(3,nup-delta-5,nup-delta-5-gap,cuts,-1);
6455 SumTracks(seeds,arr);
6456 SignClusters(seeds,fnumber,fdensity);
6467 Info("Tracking()","\n\nPrimary seeding\t%d\n\n",seeds->GetEntriesFast());
6471 // RemoveUsed(seeds,0.75,0.75,1);
6473 //SignClusters(seeds,fnumber,fdensity);
6482 arr = Tracking(4,nup-1,nup-1-gap,cuts,-1);
6483 SumTracks(seeds,arr);
6484 SignClusters(seeds,fnumber,fdensity);
6486 arr = Tracking(4,nup-2,nup-2-gap,cuts,-1);
6487 SumTracks(seeds,arr);
6488 SignClusters(seeds,fnumber,fdensity);
6490 arr = Tracking(4,nup-3,nup-3-gap,cuts,-1);
6491 SumTracks(seeds,arr);
6492 SignClusters(seeds,fnumber,fdensity);
6494 arr = Tracking(4,nup-5,nup-5-gap,cuts,-1);
6495 SumTracks(seeds,arr);
6496 SignClusters(seeds,fnumber,fdensity);
6498 arr = Tracking(4,nup-7,nup-7-gap,cuts,-1);
6499 SumTracks(seeds,arr);
6500 SignClusters(seeds,fnumber,fdensity);
6503 arr = Tracking(4,nup-9,nup-9-gap,cuts,-1);
6504 SumTracks(seeds,arr);
6505 SignClusters(seeds,fnumber,fdensity);
6509 for (Int_t delta = 9; delta<30; delta+=gapSec){
6515 arr = Tracking(4,nup-1-delta,nup-1-delta-gap,cuts,-1);
6516 SumTracks(seeds,arr);
6517 SignClusters(seeds,fnumber,fdensity);
6519 arr = Tracking(4,nup-3-delta,nup-5-delta-gap,cuts,4);
6520 SumTracks(seeds,arr);
6521 SignClusters(seeds,fnumber,fdensity);
6534 for (Int_t delta = 30; delta<fLastSeedRowSec; delta+=gapSec){
6540 arr = Tracking(4,nup-1-delta,nup-1-delta-gap,cuts,-1);
6541 SumTracks(seeds,arr);
6542 SignClusters(seeds,fnumber,fdensity);
6544 arr = Tracking(4,nup-5-delta,nup-5-delta-gap,cuts,5 );
6545 SumTracks(seeds,arr);
6546 SignClusters(seeds,fnumber,fdensity);
6550 Info("Tracking()","\n\nSecondary seeding\t%d\n\n",seeds->GetEntriesFast());
6561 TObjArray * AliTPCtrackerMI::TrackingSpecial()
6564 // seeding adjusted for laser and cosmic tests - short tracks with big inclination angle
6565 // no primary vertex seeding tried
6569 Int_t nup=fOuterSec->GetNRows()+fInnerSec->GetNRows();
6571 TObjArray * seeds = new TObjArray;
6576 Float_t fnumber = 3.0;
6577 Float_t fdensity = 3.0;
6580 cuts[0] = AliTPCReconstructor::GetRecoParam()->GetMaxC(); // max curvature
6581 cuts[1] = 3.5; // max tan(phi) angle for seeding
6582 cuts[2] = 3.; // not used (cut on z primary vertex)
6583 cuts[3] = 3.5; // max tan(theta) angle for seeding
6585 for (Int_t delta = 0; nup-delta-gap-1>0; delta+=3){
6587 arr = Tracking(4,nup-1-delta,nup-1-delta-gap,cuts,-1);
6588 SumTracks(seeds,arr);
6589 SignClusters(seeds,fnumber,fdensity);
6593 Info("Tracking()","\n\nSecondary seeding\t%d\n\n",seeds->GetEntriesFast());
6604 void AliTPCtrackerMI::SumTracks(TObjArray *arr1,TObjArray *&arr2) const
6607 //sum tracks to common container
6608 //remove suspicious tracks
6609 Int_t nseed = arr2->GetEntriesFast();
6610 for (Int_t i=0;i<nseed;i++){
6611 AliTPCseed *pt=(AliTPCseed*)arr2->UncheckedAt(i);
6614 // remove tracks with too big curvature
6616 if (TMath::Abs(pt->GetC())>AliTPCReconstructor::GetRecoParam()->GetMaxC()){
6617 delete arr2->RemoveAt(i);
6620 // REMOVE VERY SHORT TRACKS
6621 if (pt->GetNumberOfClusters()<20){
6622 delete arr2->RemoveAt(i);
6625 // NORMAL ACTIVE TRACK
6626 if (pt->IsActive()){
6627 arr1->AddLast(arr2->RemoveAt(i));
6630 //remove not usable tracks
6631 if (pt->GetRemoval()!=10){
6632 delete arr2->RemoveAt(i);
6636 // ENABLE ONLY ENOUGH GOOD STOPPED TRACKS
6637 if (pt->GetDensityFirst(20)>0.8 || pt->GetDensityFirst(30)>0.8 || pt->GetDensityFirst(40)>0.7)
6638 arr1->AddLast(arr2->RemoveAt(i));
6640 delete arr2->RemoveAt(i);
6644 delete arr2; arr2 = 0;
6649 void AliTPCtrackerMI::ParallelTracking(TObjArray *const arr, Int_t rfirst, Int_t rlast)
6652 // try to track in parralel
6654 Int_t nseed=arr->GetEntriesFast();
6655 //prepare seeds for tracking
6656 for (Int_t i=0; i<nseed; i++) {
6657 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i), &t=*pt;
6659 if (!t.IsActive()) continue;
6660 // follow prolongation to the first layer
6661 if ( (fSectors ==fInnerSec) || (t.GetFirstPoint()-fkParam->GetNRowLow()>rfirst+1) )
6662 FollowProlongation(t, rfirst+1);
6667 for (Int_t nr=rfirst; nr>=rlast; nr--){
6668 if (nr<fInnerSec->GetNRows())
6669 fSectors = fInnerSec;
6671 fSectors = fOuterSec;
6672 // make indexes with the cluster tracks for given
6674 // find nearest cluster
6675 for (Int_t i=0; i<nseed; i++) {
6676 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i), &t=*pt;
6678 if (nr==80) pt->UpdateReference();
6679 if (!pt->IsActive()) continue;
6680 // if ( (fSectors ==fOuterSec) && (pt->fFirstPoint-fkParam->GetNRowLow())<nr) continue;
6681 if (pt->GetRelativeSector()>17) {
6684 UpdateClusters(t,nr);
6686 // prolonagate to the nearest cluster - if founded
6687 for (Int_t i=0; i<nseed; i++) {
6688 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
6690 if (!pt->IsActive()) continue;
6691 // if ((fSectors ==fOuterSec) && (pt->fFirstPoint-fkParam->GetNRowLow())<nr) continue;
6692 if (pt->GetRelativeSector()>17) {
6695 FollowToNextCluster(*pt,nr);
6700 void AliTPCtrackerMI::PrepareForBackProlongation(TObjArray *const arr,Float_t fac) const
6704 // if we use TPC track itself we have to "update" covariance
6706 Int_t nseed= arr->GetEntriesFast();
6707 for (Int_t i=0;i<nseed;i++){
6708 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
6712 //rotate to current local system at first accepted point
6713 Int_t index = pt->GetClusterIndex2(pt->GetFirstPoint());
6714 Int_t sec = (index&0xff000000)>>24;
6716 Float_t angle1 = fInnerSec->GetAlpha()*sec+fInnerSec->GetAlphaShift();
6717 if (angle1>TMath::Pi())
6718 angle1-=2.*TMath::Pi();
6719 Float_t angle2 = pt->GetAlpha();
6721 if (TMath::Abs(angle1-angle2)>0.001){
6722 pt->Rotate(angle1-angle2);
6723 //angle2 = pt->GetAlpha();
6724 //pt->fRelativeSector = pt->GetAlpha()/fInnerSec->GetAlpha();
6725 //if (pt->GetAlpha()<0)
6726 // pt->fRelativeSector+=18;
6727 //sec = pt->fRelativeSector;
6736 void AliTPCtrackerMI::PrepareForProlongation(TObjArray *const arr, Float_t fac) const
6740 // if we use TPC track itself we have to "update" covariance
6742 Int_t nseed= arr->GetEntriesFast();
6743 for (Int_t i=0;i<nseed;i++){
6744 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
6747 pt->SetFirstPoint(pt->GetLastPoint());
6755 Int_t AliTPCtrackerMI::PropagateBack(TObjArray *const arr)
6758 // make back propagation
6760 Int_t nseed= arr->GetEntriesFast();
6761 for (Int_t i=0;i<nseed;i++){
6762 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
6763 if (pt&& pt->GetKinkIndex(0)<=0) {
6764 //AliTPCseed *pt2 = new AliTPCseed(*pt);
6765 fSectors = fInnerSec;
6766 //FollowBackProlongation(*pt,fInnerSec->GetNRows()-1);
6767 //fSectors = fOuterSec;
6768 FollowBackProlongation(*pt,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1);
6769 //if (pt->GetNumberOfClusters()<(pt->fEsd->GetTPCclusters(0)) ){
6770 // Error("PropagateBack","Not prolonged track %d",pt->GetLabel());
6771 // FollowBackProlongation(*pt2,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1);
6774 if (pt&& pt->GetKinkIndex(0)>0) {
6775 AliESDkink * kink = fEvent->GetKink(pt->GetKinkIndex(0)-1);
6776 pt->SetFirstPoint(kink->GetTPCRow0());
6777 fSectors = fInnerSec;
6778 FollowBackProlongation(*pt,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1);
6786 Int_t AliTPCtrackerMI::PropagateForward2(TObjArray *const arr)
6789 // make forward propagation
6791 Int_t nseed= arr->GetEntriesFast();
6793 for (Int_t i=0;i<nseed;i++){
6794 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
6796 FollowProlongation(*pt,0);
6805 Int_t AliTPCtrackerMI::PropagateForward()
6808 // propagate track forward
6810 Int_t nseed = fSeeds->GetEntriesFast();
6811 for (Int_t i=0;i<nseed;i++){
6812 AliTPCseed *pt = (AliTPCseed*)fSeeds->UncheckedAt(i);
6814 AliTPCseed &t = *pt;
6815 Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
6816 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
6817 if (alpha < 0. ) alpha += 2.*TMath::Pi();
6818 t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN);
6822 fSectors = fOuterSec;
6823 ParallelTracking(fSeeds,fOuterSec->GetNRows()+fInnerSec->GetNRows()-1,fInnerSec->GetNRows());
6824 fSectors = fInnerSec;
6825 ParallelTracking(fSeeds,fInnerSec->GetNRows()-1,0);
6834 Int_t AliTPCtrackerMI::PropagateBack(AliTPCseed *const pt, Int_t row0, Int_t row1)
6837 // make back propagation, in between row0 and row1
6841 fSectors = fInnerSec;
6844 if (row1<fSectors->GetNRows())
6847 r1 = fSectors->GetNRows()-1;
6849 if (row0<fSectors->GetNRows()&& r1>0 )
6850 FollowBackProlongation(*pt,r1);
6851 if (row1<=fSectors->GetNRows())
6854 r1 = row1 - fSectors->GetNRows();
6855 if (r1<=0) return 0;
6856 if (r1>=fOuterSec->GetNRows()) return 0;
6857 fSectors = fOuterSec;
6858 return FollowBackProlongation(*pt,r1);
6866 void AliTPCtrackerMI::GetShape(AliTPCseed * seed, Int_t row)
6870 AliTPCClusterParam * clparam = AliTPCcalibDB::Instance()->GetClusterParam();
6871 Float_t zdrift = TMath::Abs((fkParam->GetZLength(0)-TMath::Abs(seed->GetZ())));
6872 Int_t type = (seed->GetSector() < fkParam->GetNSector()/2) ? 0: (row>126) ? 1:2;
6873 Double_t angulary = seed->GetSnp();
6875 if (TMath::Abs(angulary)>AliTPCReconstructor::GetMaxSnpTracker()) {
6876 angulary = TMath::Sign(AliTPCReconstructor::GetMaxSnpTracker(),angulary);
6879 angulary = angulary*angulary/((1.-angulary)*(1.+angulary));
6880 Double_t angularz = seed->GetTgl()*seed->GetTgl()*(1.+angulary);
6882 Double_t sigmay = clparam->GetRMS0(0,type,zdrift,TMath::Sqrt(TMath::Abs(angulary)));
6883 Double_t sigmaz = clparam->GetRMS0(1,type,zdrift,TMath::Sqrt(TMath::Abs(angularz)));
6884 seed->SetCurrentSigmaY2(sigmay*sigmay);
6885 seed->SetCurrentSigmaZ2(sigmaz*sigmaz);
6886 // Float_t sd2 = TMath::Abs((fkParam->GetZLength(0)-TMath::Abs(seed->GetZ())))*fkParam->GetDiffL()*fkParam->GetDiffL();
6887 // // Float_t padlength = fkParam->GetPadPitchLength(seed->fSector);
6888 // Float_t padlength = GetPadPitchLength(row);
6890 // Float_t sresy = (seed->GetSector() < fkParam->GetNSector()/2) ? 0.2 :0.3;
6891 // seed->SetCurrentSigmaY2(sd2+padlength*padlength*angulary/12.+sresy*sresy);
6893 // Float_t sresz = fkParam->GetZSigma();
6894 // seed->SetCurrentSigmaZ2(sd2+padlength*padlength*angularz*angularz*(1+angulary)/12.+sresz*sresz);
6896 Float_t wy = GetSigmaY(seed);
6897 Float_t wz = GetSigmaZ(seed);
6900 if (TMath::Abs(wy/seed->fCurrentSigmaY2-1)>0.0001 || TMath::Abs(wz/seed->fCurrentSigmaZ2-1)>0.0001 ){
6901 printf("problem\n");
6908 //__________________________________________________________________________
6909 void AliTPCtrackerMI::CookLabel(AliKalmanTrack *tk, Float_t wrong) const {
6910 //--------------------------------------------------------------------
6911 //This function "cooks" a track label. If label<0, this track is fake.
6912 //--------------------------------------------------------------------
6913 AliTPCseed * t = dynamic_cast<AliTPCseed*>(tk);
6915 printf("%s:%d wrong type \n",(char*)__FILE__,__LINE__);
6919 Int_t noc=t->GetNumberOfClusters();
6921 //printf("\nnot founded prolongation\n\n\n");
6927 AliTPCclusterMI *clusters[160];
6929 for (Int_t i=0;i<160;i++) {
6936 for (i=0; i<160 && current<noc; i++) {
6938 Int_t index=t->GetClusterIndex2(i);
6939 if (index<=0) continue;
6940 if (index&0x8000) continue;
6942 //clusters[current]=GetClusterMI(index);
6943 if (t->GetClusterPointer(i)){
6944 clusters[current]=t->GetClusterPointer(i);
6950 Int_t lab=123456789;
6951 for (i=0; i<noc; i++) {
6952 AliTPCclusterMI *c=clusters[i];
6954 lab=TMath::Abs(c->GetLabel(0));
6956 for (j=0; j<noc; j++) if (lb[j]==lab || mx[j]==0) break;
6962 for (i=0; i<noc; i++) if (mx[i]>max) {max=mx[i]; lab=lb[i];}
6964 for (i=0; i<noc; i++) {
6965 AliTPCclusterMI *c=clusters[i];
6967 if (TMath::Abs(c->GetLabel(1)) == lab ||
6968 TMath::Abs(c->GetLabel(2)) == lab ) max++;
6970 if (noc<=0) { lab=-1; return;}
6971 if ((1.- Float_t(max)/(noc)) > wrong) lab=-lab;
6974 Int_t tail=Int_t(0.10*noc);
6977 for (i=1; i<=160&&ind<tail; i++) {
6978 // AliTPCclusterMI *c=clusters[noc-i];
6979 AliTPCclusterMI *c=clusters[i];
6981 if (lab == TMath::Abs(c->GetLabel(0)) ||
6982 lab == TMath::Abs(c->GetLabel(1)) ||
6983 lab == TMath::Abs(c->GetLabel(2))) max++;
6986 if (max < Int_t(0.5*tail)) lab=-lab;
6993 //delete[] clusters;
6997 //__________________________________________________________________________
6998 Int_t AliTPCtrackerMI::CookLabel(AliTPCseed *const t, Float_t wrong,Int_t first, Int_t last) const {
6999 //--------------------------------------------------------------------
7000 //This function "cooks" a track label. If label<0, this track is fake.
7001 //--------------------------------------------------------------------
7002 Int_t noc=t->GetNumberOfClusters();
7004 //printf("\nnot founded prolongation\n\n\n");
7010 AliTPCclusterMI *clusters[160];
7012 for (Int_t i=0;i<160;i++) {
7019 for (i=0; i<160 && current<noc; i++) {
7020 if (i<first) continue;
7021 if (i>last) continue;
7022 Int_t index=t->GetClusterIndex2(i);
7023 if (index<=0) continue;
7024 if (index&0x8000) continue;
7026 //clusters[current]=GetClusterMI(index);
7027 if (t->GetClusterPointer(i)){
7028 clusters[current]=t->GetClusterPointer(i);
7033 if (noc<5) return -1;
7034 Int_t lab=123456789;
7035 for (i=0; i<noc; i++) {
7036 AliTPCclusterMI *c=clusters[i];
7038 lab=TMath::Abs(c->GetLabel(0));
7040 for (j=0; j<noc; j++) if (lb[j]==lab || mx[j]==0) break;
7046 for (i=0; i<noc; i++) if (mx[i]>max) {max=mx[i]; lab=lb[i];}
7048 for (i=0; i<noc; i++) {
7049 AliTPCclusterMI *c=clusters[i];
7051 if (TMath::Abs(c->GetLabel(1)) == lab ||
7052 TMath::Abs(c->GetLabel(2)) == lab ) max++;
7054 if (noc<=0) { lab=-1; return -1;}
7055 if ((1.- Float_t(max)/(noc)) > wrong) lab=-lab;
7058 Int_t tail=Int_t(0.10*noc);
7061 for (i=1; i<=160&&ind<tail; i++) {
7062 // AliTPCclusterMI *c=clusters[noc-i];
7063 AliTPCclusterMI *c=clusters[i];
7065 if (lab == TMath::Abs(c->GetLabel(0)) ||
7066 lab == TMath::Abs(c->GetLabel(1)) ||
7067 lab == TMath::Abs(c->GetLabel(2))) max++;
7070 if (max < Int_t(0.5*tail)) lab=-lab;
7073 // t->SetLabel(lab);
7077 //delete[] clusters;
7081 Int_t AliTPCtrackerMI::GetRowNumber(Double_t x[3]) const
7083 //return pad row number for given x vector
7084 Float_t phi = TMath::ATan2(x[1],x[0]);
7085 if(phi<0) phi=2.*TMath::Pi()+phi;
7086 // Get the local angle in the sector philoc
7087 const Float_t kRaddeg = 180/3.14159265358979312;
7088 Float_t phiangle = (Int_t (phi*kRaddeg/20.) + 0.5)*20./kRaddeg;
7089 Double_t localx = x[0]*TMath::Cos(phiangle)-x[1]*TMath::Sin(phiangle);
7090 return GetRowNumber(localx);
7095 void AliTPCtrackerMI::MakeBitmaps(AliTPCseed *t)
7097 //-----------------------------------------------------------------------
7098 // Fill the cluster and sharing bitmaps of the track
7099 //-----------------------------------------------------------------------
7101 Int_t firstpoint = 0;
7102 Int_t lastpoint = 159;
7103 AliTPCTrackerPoint *point;
7104 AliTPCclusterMI *cluster;
7106 for (int iter=firstpoint; iter<lastpoint; iter++) {
7107 // Change to cluster pointers to see if we have a cluster at given padrow
7108 cluster = t->GetClusterPointer(iter);
7110 t->SetClusterMapBit(iter, kTRUE);
7111 point = t->GetTrackPoint(iter);
7112 if (point->IsShared())
7113 t->SetSharedMapBit(iter,kTRUE);
7115 t->SetSharedMapBit(iter, kFALSE);
7118 t->SetClusterMapBit(iter, kFALSE);
7119 t->SetSharedMapBit(iter, kFALSE);
7124 Bool_t AliTPCtrackerMI::IsFindable(AliTPCseed & track){
7126 // return flag if there is findable cluster at given position
7129 Float_t z = track.GetZ();
7131 if (TMath::Abs(z)<(AliTPCReconstructor::GetCtgRange()*track.GetX()+kDeltaZ) &&
7132 TMath::Abs(z)<fkParam->GetZLength(0) &&
7133 (TMath::Abs(track.GetSnp())<AliTPCReconstructor::GetMaxSnpTracker()))
7139 void AliTPCtrackerMI::AddCovariance(AliTPCseed * seed){
7141 // Adding systematic error
7142 // !!!! the systematic error for element 4 is in 1/cm not in pt
7144 const Double_t *param = AliTPCReconstructor::GetRecoParam()->GetSystematicError();
7145 Double_t *covarIn= (Double_t*)seed->GetCovariance();
7147 for (Int_t i=0;i<15;i++) covar[i]=0;
7153 covar[0] = param[0]*param[0];
7154 covar[2] = param[1]*param[1];
7155 covar[5] = param[2]*param[2];
7156 covar[9] = param[3]*param[3];
7157 Double_t facC = AliTracker::GetBz()*kB2C;
7158 covar[14]= param[4]*param[4]*facC*facC;
7160 covar[1]=TMath::Sqrt((covar[0]*covar[2]))*covarIn[1]/TMath::Sqrt((covarIn[0]*covarIn[2]));
7162 covar[3]=TMath::Sqrt((covar[0]*covar[5]))*covarIn[3]/TMath::Sqrt((covarIn[0]*covarIn[5]));
7163 covar[4]=TMath::Sqrt((covar[2]*covar[5]))*covarIn[4]/TMath::Sqrt((covarIn[2]*covarIn[5]));
7165 covar[6]=TMath::Sqrt((covar[0]*covar[9]))*covarIn[6]/TMath::Sqrt((covarIn[0]*covarIn[9]));
7166 covar[7]=TMath::Sqrt((covar[2]*covar[9]))*covarIn[7]/TMath::Sqrt((covarIn[2]*covarIn[9]));
7167 covar[8]=TMath::Sqrt((covar[5]*covar[9]))*covarIn[8]/TMath::Sqrt((covarIn[5]*covarIn[9]));
7169 covar[10]=TMath::Sqrt((covar[0]*covar[14]))*covarIn[10]/TMath::Sqrt((covarIn[0]*covarIn[14]));
7170 covar[11]=TMath::Sqrt((covar[2]*covar[14]))*covarIn[11]/TMath::Sqrt((covarIn[2]*covarIn[14]));
7171 covar[12]=TMath::Sqrt((covar[5]*covar[14]))*covarIn[12]/TMath::Sqrt((covarIn[5]*covarIn[14]));
7172 covar[13]=TMath::Sqrt((covar[9]*covar[14]))*covarIn[13]/TMath::Sqrt((covarIn[9]*covarIn[14]));
7174 seed->AddCovariance(covar);