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 = 0x0;
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)));
1178 Int_t AliTPCtrackerMI::LoadClusters(const TClonesArray *arr)
1181 // load clusters to the memory from one
1184 AliTPCclusterMI *clust=0;
1185 Int_t count[72][96] = { {0} , {0} };
1187 // loop over clusters
1188 for (Int_t icl=0; icl<arr->GetEntriesFast(); icl++) {
1189 clust = (AliTPCclusterMI*)arr->At(icl);
1190 if(!clust) continue;
1191 //printf("cluster: det %d, row %d \n", clust->GetDetector(),clust->GetRow());
1193 // transform clusters
1196 // count clusters per pad row
1197 count[clust->GetDetector()][clust->GetRow()]++;
1200 // insert clusters to sectors
1201 for (Int_t icl=0; icl<arr->GetEntriesFast(); icl++) {
1202 clust = (AliTPCclusterMI*)arr->At(icl);
1203 if(!clust) continue;
1205 Int_t sec = clust->GetDetector();
1206 Int_t row = clust->GetRow();
1208 // filter overlapping pad rows needed by HLT
1209 if(sec<fkNIS*2) { //IROCs
1210 if(row == 30) continue;
1213 if(row == 27 || row == 76) continue;
1219 fInnerSec[sec%fkNIS].InsertCluster(clust, count[sec][row], fkParam);
1222 left = (sec-fkNIS*2)/fkNOS;
1223 fOuterSec[(sec-fkNIS*2)%fkNOS].InsertCluster(clust, count[sec][row], fkParam);
1227 // Load functions must be called behind LoadCluster(TClonesArray*)
1229 //LoadOuterSectors();
1230 //LoadInnerSectors();
1236 Int_t AliTPCtrackerMI::LoadClusters()
1239 // load clusters to the memory
1240 AliTPCClustersRow *clrow= new AliTPCClustersRow;
1241 clrow->SetClass("AliTPCclusterMI");
1243 clrow->GetArray()->ExpandCreateFast(10000);
1245 // TTree * tree = fClustersArray.GetTree();
1247 TTree * tree = fInput;
1248 TBranch * br = tree->GetBranch("Segment");
1249 br->SetAddress(&clrow);
1251 Int_t j=Int_t(tree->GetEntries());
1252 for (Int_t i=0; i<j; i++) {
1256 fkParam->AdjustSectorRow(clrow->GetID(),sec,row);
1257 for (Int_t icl=0; icl<clrow->GetArray()->GetEntriesFast(); icl++){
1258 Transform((AliTPCclusterMI*)(clrow->GetArray()->At(icl)));
1261 AliTPCtrackerRow * tpcrow=0;
1264 tpcrow = &(fInnerSec[sec%fkNIS][row]);
1268 tpcrow = &(fOuterSec[(sec-fkNIS*2)%fkNOS][row]);
1269 left = (sec-fkNIS*2)/fkNOS;
1272 tpcrow->SetN1(clrow->GetArray()->GetEntriesFast());
1273 tpcrow->SetClusters1(new AliTPCclusterMI[tpcrow->GetN1()]);
1274 for (Int_t k=0;k<tpcrow->GetN1();++k)
1275 tpcrow->SetCluster1(k, *(AliTPCclusterMI*)(clrow->GetArray()->At(k)));
1278 tpcrow->SetN2(clrow->GetArray()->GetEntriesFast());
1279 tpcrow->SetClusters2(new AliTPCclusterMI[tpcrow->GetN2()]);
1280 for (Int_t k=0;k<tpcrow->GetN2();++k)
1281 tpcrow->SetCluster2(k,*(AliTPCclusterMI*)(clrow->GetArray()->At(k)));
1292 void AliTPCtrackerMI::UnloadClusters()
1295 // unload clusters from the memory
1297 Int_t nrows = fOuterSec->GetNRows();
1298 for (Int_t sec = 0;sec<fkNOS;sec++)
1299 for (Int_t row = 0;row<nrows;row++){
1300 AliTPCtrackerRow* tpcrow = &(fOuterSec[sec%fkNOS][row]);
1302 // if (tpcrow->fClusters1) delete []tpcrow->fClusters1;
1303 // if (tpcrow->fClusters2) delete []tpcrow->fClusters2;
1305 tpcrow->ResetClusters();
1308 nrows = fInnerSec->GetNRows();
1309 for (Int_t sec = 0;sec<fkNIS;sec++)
1310 for (Int_t row = 0;row<nrows;row++){
1311 AliTPCtrackerRow* tpcrow = &(fInnerSec[sec%fkNIS][row]);
1313 // if (tpcrow->fClusters1) delete []tpcrow->fClusters1;
1314 //if (tpcrow->fClusters2) delete []tpcrow->fClusters2;
1316 tpcrow->ResetClusters();
1322 void AliTPCtrackerMI::FillClusterArray(TObjArray* array) const{
1324 // Filling cluster to the array - For visualization purposes
1327 nrows = fOuterSec->GetNRows();
1328 for (Int_t sec = 0;sec<fkNOS;sec++)
1329 for (Int_t row = 0;row<nrows;row++){
1330 AliTPCtrackerRow* tpcrow = &(fOuterSec[sec%fkNOS][row]);
1331 if (!tpcrow) continue;
1332 for (Int_t icl = 0;icl<tpcrow->GetN();icl++){
1333 array->AddLast((TObject*)((*tpcrow)[icl]));
1336 nrows = fInnerSec->GetNRows();
1337 for (Int_t sec = 0;sec<fkNIS;sec++)
1338 for (Int_t row = 0;row<nrows;row++){
1339 AliTPCtrackerRow* tpcrow = &(fInnerSec[sec%fkNIS][row]);
1340 if (!tpcrow) continue;
1341 for (Int_t icl = 0;icl<tpcrow->GetN();icl++){
1342 array->AddLast((TObject*)(*tpcrow)[icl]);
1348 void AliTPCtrackerMI::Transform(AliTPCclusterMI * cluster){
1352 AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ;
1354 AliFatal("Tranformations not in calibDB");
1356 transform->SetCurrentRecoParam((AliTPCRecoParam*)AliTPCReconstructor::GetRecoParam());
1357 Double_t x[3]={cluster->GetRow(),cluster->GetPad(),cluster->GetTimeBin()};
1358 Int_t i[1]={cluster->GetDetector()};
1359 transform->Transform(x,i,0,1);
1360 // if (cluster->GetDetector()%36>17){
1365 // in debug mode check the transformation
1367 if (AliTPCReconstructor::StreamLevel()>2) {
1369 cluster->GetGlobalXYZ(gx);
1370 Int_t event = (fEvent==NULL)? 0: fEvent->GetEventNumberInFile();
1371 TTreeSRedirector &cstream = *fDebugStreamer;
1372 cstream<<"Transform"<<
1383 cluster->SetX(x[0]);
1384 cluster->SetY(x[1]);
1385 cluster->SetZ(x[2]);
1391 //if (!fkParam->IsGeoRead()) fkParam->ReadGeoMatrices();
1392 TGeoHMatrix *mat = fkParam->GetClusterMatrix(cluster->GetDetector());
1394 Double_t pos[3]= {cluster->GetX(),cluster->GetY(),cluster->GetZ()};
1395 Double_t posC[3]={cluster->GetX(),cluster->GetY(),cluster->GetZ()};
1396 if (mat) mat->LocalToMaster(pos,posC);
1398 // chack Loading of Geo matrices from GeoManager - TEMPORARY FIX
1400 cluster->SetX(posC[0]);
1401 cluster->SetY(posC[1]);
1402 cluster->SetZ(posC[2]);
1405 //_____________________________________________________________________________
1406 Int_t AliTPCtrackerMI::LoadOuterSectors() {
1407 //-----------------------------------------------------------------
1408 // This function fills outer TPC sectors with clusters.
1409 //-----------------------------------------------------------------
1410 Int_t nrows = fOuterSec->GetNRows();
1412 for (Int_t sec = 0;sec<fkNOS;sec++)
1413 for (Int_t row = 0;row<nrows;row++){
1414 AliTPCtrackerRow* tpcrow = &(fOuterSec[sec%fkNOS][row]);
1415 Int_t sec2 = sec+2*fkNIS;
1417 Int_t ncl = tpcrow->GetN1();
1419 AliTPCclusterMI *c= (tpcrow->GetCluster1(ncl));
1420 index=(((sec2<<8)+row)<<16)+ncl;
1421 tpcrow->InsertCluster(c,index);
1424 ncl = tpcrow->GetN2();
1426 AliTPCclusterMI *c= (tpcrow->GetCluster2(ncl));
1427 index=((((sec2+fkNOS)<<8)+row)<<16)+ncl;
1428 tpcrow->InsertCluster(c,index);
1431 // write indexes for fast acces
1433 for (Int_t i=0;i<510;i++)
1434 tpcrow->SetFastCluster(i,-1);
1435 for (Int_t i=0;i<tpcrow->GetN();i++){
1436 Int_t zi = Int_t((*tpcrow)[i]->GetZ()+255.);
1437 tpcrow->SetFastCluster(zi,i); // write index
1440 for (Int_t i=0;i<510;i++){
1441 if (tpcrow->GetFastCluster(i)<0)
1442 tpcrow->SetFastCluster(i,last);
1444 last = tpcrow->GetFastCluster(i);
1453 //_____________________________________________________________________________
1454 Int_t AliTPCtrackerMI::LoadInnerSectors() {
1455 //-----------------------------------------------------------------
1456 // This function fills inner TPC sectors with clusters.
1457 //-----------------------------------------------------------------
1458 Int_t nrows = fInnerSec->GetNRows();
1460 for (Int_t sec = 0;sec<fkNIS;sec++)
1461 for (Int_t row = 0;row<nrows;row++){
1462 AliTPCtrackerRow* tpcrow = &(fInnerSec[sec%fkNIS][row]);
1465 Int_t ncl = tpcrow->GetN1();
1467 AliTPCclusterMI *c= (tpcrow->GetCluster1(ncl));
1468 index=(((sec<<8)+row)<<16)+ncl;
1469 tpcrow->InsertCluster(c,index);
1472 ncl = tpcrow->GetN2();
1474 AliTPCclusterMI *c= (tpcrow->GetCluster2(ncl));
1475 index=((((sec+fkNIS)<<8)+row)<<16)+ncl;
1476 tpcrow->InsertCluster(c,index);
1479 // write indexes for fast acces
1481 for (Int_t i=0;i<510;i++)
1482 tpcrow->SetFastCluster(i,-1);
1483 for (Int_t i=0;i<tpcrow->GetN();i++){
1484 Int_t zi = Int_t((*tpcrow)[i]->GetZ()+255.);
1485 tpcrow->SetFastCluster(zi,i); // write index
1488 for (Int_t i=0;i<510;i++){
1489 if (tpcrow->GetFastCluster(i)<0)
1490 tpcrow->SetFastCluster(i,last);
1492 last = tpcrow->GetFastCluster(i);
1504 //_________________________________________________________________________
1505 AliTPCclusterMI *AliTPCtrackerMI::GetClusterMI(Int_t index) const {
1506 //--------------------------------------------------------------------
1507 // Return pointer to a given cluster
1508 //--------------------------------------------------------------------
1509 if (index<0) return 0; // no cluster
1510 Int_t sec=(index&0xff000000)>>24;
1511 Int_t row=(index&0x00ff0000)>>16;
1512 Int_t ncl=(index&0x00007fff)>>00;
1514 const AliTPCtrackerRow * tpcrow=0;
1515 AliTPCclusterMI * clrow =0;
1517 if (sec<0 || sec>=fkNIS*4) {
1518 AliWarning(Form("Wrong sector %d",sec));
1523 AliTPCtrackerSector& tracksec = fInnerSec[sec%fkNIS];
1524 if (tracksec.GetNRows()<=row) return 0;
1525 tpcrow = &(tracksec[row]);
1526 if (tpcrow==0) return 0;
1529 if (tpcrow->GetN1()<=ncl) return 0;
1530 clrow = tpcrow->GetClusters1();
1533 if (tpcrow->GetN2()<=ncl) return 0;
1534 clrow = tpcrow->GetClusters2();
1538 AliTPCtrackerSector& tracksec = fOuterSec[(sec-fkNIS*2)%fkNOS];
1539 if (tracksec.GetNRows()<=row) return 0;
1540 tpcrow = &(tracksec[row]);
1541 if (tpcrow==0) return 0;
1543 if (sec-2*fkNIS<fkNOS) {
1544 if (tpcrow->GetN1()<=ncl) return 0;
1545 clrow = tpcrow->GetClusters1();
1548 if (tpcrow->GetN2()<=ncl) return 0;
1549 clrow = tpcrow->GetClusters2();
1553 return &(clrow[ncl]);
1559 Int_t AliTPCtrackerMI::FollowToNext(AliTPCseed& t, Int_t nr) {
1560 //-----------------------------------------------------------------
1561 // This function tries to find a track prolongation to next pad row
1562 //-----------------------------------------------------------------
1564 Double_t x= GetXrow(nr), ymax=GetMaxY(nr);
1567 AliTPCclusterMI *cl=0;
1568 Int_t tpcindex= t.GetClusterIndex2(nr);
1570 // update current shape info every 5 pad-row
1571 // if ( (nr%5==0) || t.GetNumberOfClusters()<2 || (t.fCurrentSigmaY2<0.0001) ){
1575 if (fIteration>0 && tpcindex>=-1){ //if we have already clusters
1577 if (tpcindex==-1) return 0; //track in dead zone
1579 cl = t.GetClusterPointer(nr);
1580 if ( (cl==0) ) cl = GetClusterMI(tpcindex);
1581 t.SetCurrentClusterIndex1(tpcindex);
1584 Int_t relativesector = ((tpcindex&0xff000000)>>24)%18; // if previously accepted cluster in different sector
1585 Float_t angle = relativesector*fSectors->GetAlpha()+fSectors->GetAlphaShift();
1587 if (angle<-TMath::Pi()) angle += 2*TMath::Pi();
1588 if (angle>=TMath::Pi()) angle -= 2*TMath::Pi();
1590 if (TMath::Abs(angle-t.GetAlpha())>0.001){
1591 Double_t rotation = angle-t.GetAlpha();
1592 t.SetRelativeSector(relativesector);
1593 if (!t.Rotate(rotation)) return 0;
1595 if (!t.PropagateTo(x)) return 0;
1597 t.SetCurrentCluster(cl);
1599 Int_t accept = AcceptCluster(&t,t.GetCurrentCluster());
1600 if ((tpcindex&0x8000)==0) accept =0;
1602 //if founded cluster is acceptible
1603 if (cl->IsUsed(11)) { // id cluster is shared inrease uncertainty
1604 t.SetErrorY2(t.GetErrorY2()+0.03);
1605 t.SetErrorZ2(t.GetErrorZ2()+0.03);
1606 t.SetErrorY2(t.GetErrorY2()*3);
1607 t.SetErrorZ2(t.GetErrorZ2()*3);
1609 t.SetNFoundable(t.GetNFoundable()+1);
1610 UpdateTrack(&t,accept);
1615 if (TMath::Abs(t.GetSnp())>AliTPCReconstructor::GetMaxSnpTracker()) return 0; // cut on angle
1616 if (fIteration>1 && IsFindable(t)){
1617 // not look for new cluster during refitting
1618 t.SetNFoundable(t.GetNFoundable()+1);
1623 // if (TMath::Abs(t.GetSnp())>0.95 || TMath::Abs(x*t.GetC()-t.GetEta())>0.95) return 0;// patch 28 fev 06
1624 Double_t y=t.GetYat(x);
1625 if (TMath::Abs(y)>ymax){
1627 t.SetRelativeSector((t.GetRelativeSector()+1) % fN);
1628 if (!t.Rotate(fSectors->GetAlpha()))
1630 } else if (y <-ymax) {
1631 t.SetRelativeSector((t.GetRelativeSector()-1+fN) % fN);
1632 if (!t.Rotate(-fSectors->GetAlpha()))
1638 if (!t.PropagateTo(x)) {
1639 if (fIteration==0) t.SetRemoval(10);
1643 Double_t z=t.GetZ();
1646 if (!IsActive(t.GetRelativeSector(),nr)) {
1648 t.SetClusterIndex2(nr,-1);
1651 //AliInfo(Form("A - Sector%d phi %f - alpha %f", t.fRelativeSector,y/x, t.GetAlpha()));
1652 Bool_t isActive = IsActive(t.GetRelativeSector(),nr);
1653 Bool_t isActive2 = (nr>=fInnerSec->GetNRows()) ? fOuterSec[t.GetRelativeSector()][nr-fInnerSec->GetNRows()].GetN()>0:fInnerSec[t.GetRelativeSector()][nr].GetN()>0;
1655 if (!isActive || !isActive2) return 0;
1657 const AliTPCtrackerRow &krow=GetRow(t.GetRelativeSector(),nr);
1658 if ( (t.GetSigmaY2()<0) || t.GetSigmaZ2()<0) return 0;
1660 Double_t roadz = 1.;
1662 if (TMath::Abs(TMath::Abs(y)-ymax)<krow.GetDeadZone()){
1664 t.SetClusterIndex2(nr,-1);
1670 // if (TMath::Abs(z)<(AliTPCReconstructor::GetCtgRange()*x+10) && TMath::Abs(z)<fkParam->GetZLength(0) && (TMath::Abs(t.GetSnp())<AliTPCReconstructor::GetMaxSnpTracker()))
1671 t.SetNFoundable(t.GetNFoundable()+1);
1677 // cl = krow.FindNearest2(y+10.,z,roady,roadz,index);
1678 cl = krow.FindNearest2(y,z,roady,roadz,index);
1679 if (cl) t.SetCurrentClusterIndex1(krow.GetIndex(index));
1682 t.SetCurrentCluster(cl);
1684 if (fIteration==2&&cl->IsUsed(10)) return 0;
1685 Int_t accept = AcceptCluster(&t,t.GetCurrentCluster());
1686 if (fIteration==2&&cl->IsUsed(11)) {
1687 t.SetErrorY2(t.GetErrorY2()+0.03);
1688 t.SetErrorZ2(t.GetErrorZ2()+0.03);
1689 t.SetErrorY2(t.GetErrorY2()*3);
1690 t.SetErrorZ2(t.GetErrorZ2()*3);
1693 if (t.fCurrentCluster->IsUsed(10)){
1698 if (t.fNShared>0.7*t.GetNumberOfClusters()) {
1704 if (accept<3) UpdateTrack(&t,accept);
1707 if ( fIteration==0 && t.GetNFoundable()*0.5 > t.GetNumberOfClusters()) t.SetRemoval(10);
1715 //_________________________________________________________________________
1716 Bool_t AliTPCtrackerMI::GetTrackPoint(Int_t index, AliTrackPoint &p ) const
1718 // Get track space point by index
1719 // return false in case the cluster doesn't exist
1720 AliTPCclusterMI *cl = GetClusterMI(index);
1721 if (!cl) return kFALSE;
1722 Int_t sector = (index&0xff000000)>>24;
1723 // Int_t row = (index&0x00ff0000)>>16;
1725 // xyz[0] = fkParam->GetPadRowRadii(sector,row);
1726 xyz[0] = cl->GetX();
1727 xyz[1] = cl->GetY();
1728 xyz[2] = cl->GetZ();
1730 fkParam->AdjustCosSin(sector,cos,sin);
1731 Float_t x = cos*xyz[0]-sin*xyz[1];
1732 Float_t y = cos*xyz[1]+sin*xyz[0];
1734 Float_t sigmaY2 = 0.027*cl->GetSigmaY2();
1735 if (sector < fkParam->GetNInnerSector()) sigmaY2 *= 2.07;
1736 Float_t sigmaZ2 = 0.066*cl->GetSigmaZ2();
1737 if (sector < fkParam->GetNInnerSector()) sigmaZ2 *= 1.77;
1738 cov[0] = sin*sin*sigmaY2;
1739 cov[1] = -sin*cos*sigmaY2;
1741 cov[3] = cos*cos*sigmaY2;
1744 p.SetXYZ(x,y,xyz[2],cov);
1745 AliGeomManager::ELayerID iLayer;
1747 if (sector < fkParam->GetNInnerSector()) {
1748 iLayer = AliGeomManager::kTPC1;
1752 iLayer = AliGeomManager::kTPC2;
1753 idet = sector - fkParam->GetNInnerSector();
1755 UShort_t volid = AliGeomManager::LayerToVolUID(iLayer,idet);
1756 p.SetVolumeID(volid);
1762 Int_t AliTPCtrackerMI::UpdateClusters(AliTPCseed& t, Int_t nr) {
1763 //-----------------------------------------------------------------
1764 // This function tries to find a track prolongation to next pad row
1765 //-----------------------------------------------------------------
1766 t.SetCurrentCluster(0);
1767 t.SetCurrentClusterIndex1(0);
1769 Double_t xt=t.GetX();
1770 Int_t row = GetRowNumber(xt)-1;
1771 Double_t ymax= GetMaxY(nr);
1773 if (row < nr) return 1; // don't prolongate if not information until now -
1774 // if (TMath::Abs(t.GetSnp())>0.9 && t.GetNumberOfClusters()>40. && fIteration!=2) {
1776 // return 0; // not prolongate strongly inclined tracks
1778 // if (TMath::Abs(t.GetSnp())>0.95) {
1780 // return 0; // not prolongate strongly inclined tracks
1781 // }// patch 28 fev 06
1783 Double_t x= GetXrow(nr);
1785 //t.PropagateTo(x+0.02);
1786 //t.PropagateTo(x+0.01);
1787 if (!t.PropagateTo(x)){
1794 if (TMath::Abs(y)>ymax){
1796 t.SetRelativeSector((t.GetRelativeSector()+1) % fN);
1797 if (!t.Rotate(fSectors->GetAlpha()))
1799 } else if (y <-ymax) {
1800 t.SetRelativeSector((t.GetRelativeSector()-1+fN) % fN);
1801 if (!t.Rotate(-fSectors->GetAlpha()))
1804 // if (!t.PropagateTo(x)){
1811 if (TMath::Abs(t.GetSnp())>AliTPCReconstructor::GetMaxSnpTracker()) return 0;
1813 if (!IsActive(t.GetRelativeSector(),nr)) {
1815 t.SetClusterIndex2(nr,-1);
1818 //AliInfo(Form("A - Sector%d phi %f - alpha %f", t.fRelativeSector,y/x, t.GetAlpha()));
1820 AliTPCtrackerRow &krow=GetRow(t.GetRelativeSector(),nr);
1822 if (TMath::Abs(TMath::Abs(y)-ymax)<krow.GetDeadZone()){
1824 t.SetClusterIndex2(nr,-1);
1830 // if (TMath::Abs(t.GetZ())<(AliTPCReconstructor::GetCtgRange()*t.GetX()+10) && (TMath::Abs(t.GetSnp())<AliTPCReconstructor::GetMaxSnpTracker()))
1831 if (IsFindable(t)) t.SetNFoundable(t.GetNFoundable()+1);
1837 if ( (nr%2==0) || t.GetNumberOfClusters()<2){
1838 // t.fCurrentSigmaY = GetSigmaY(&t);
1839 //t.fCurrentSigmaZ = GetSigmaZ(&t);
1843 AliTPCclusterMI *cl=0;
1846 Double_t roady = 1.;
1847 Double_t roadz = 1.;
1851 index = t.GetClusterIndex2(nr);
1852 if ( (index>0) && (index&0x8000)==0){
1853 cl = t.GetClusterPointer(nr);
1854 if ( (cl==0) && (index>0)) cl = GetClusterMI(index);
1855 t.SetCurrentClusterIndex1(index);
1857 t.SetCurrentCluster(cl);
1863 // if (index<0) return 0;
1864 UInt_t uindex = TMath::Abs(index);
1867 //cl = krow.FindNearest2(y+10,z,roady,roadz,uindex);
1868 cl = krow.FindNearest2(y,z,roady,roadz,uindex);
1871 if (cl) t.SetCurrentClusterIndex1(krow.GetIndex(uindex));
1872 t.SetCurrentCluster(cl);
1878 Int_t AliTPCtrackerMI::FollowToNextCluster(AliTPCseed & t, Int_t nr) {
1879 //-----------------------------------------------------------------
1880 // This function tries to find a track prolongation to next pad row
1881 //-----------------------------------------------------------------
1883 //update error according neighborhoud
1885 if (t.GetCurrentCluster()) {
1887 Int_t accept = AcceptCluster(&t,t.GetCurrentCluster());
1889 if (t.GetCurrentCluster()->IsUsed(10)){
1894 t.SetNShared(t.GetNShared()+1);
1895 if (t.GetNShared()>0.7*t.GetNumberOfClusters()) {
1900 if (fIteration>0) accept = 0;
1901 if (accept<3) UpdateTrack(&t,accept);
1905 if ( ( (t.GetSigmaY2()+t.GetSigmaZ2())>0.16)&& t.GetNumberOfClusters()>18) t.SetRemoval(10);
1906 if ( t.GetChi2()/t.GetNumberOfClusters()>6 &&t.GetNumberOfClusters()>18) t.SetRemoval(10);
1908 if (( (t.GetNFoundable()*0.5 > t.GetNumberOfClusters()) || t.GetNoCluster()>15)) t.SetRemoval(10);
1916 //_____________________________________________________________________________
1917 Int_t AliTPCtrackerMI::FollowProlongation(AliTPCseed& t, Int_t rf, Int_t step) {
1918 //-----------------------------------------------------------------
1919 // This function tries to find a track prolongation.
1920 //-----------------------------------------------------------------
1921 Double_t xt=t.GetX();
1923 Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
1924 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
1925 if (alpha < 0. ) alpha += 2.*TMath::Pi();
1927 t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN);
1929 Int_t first = GetRowNumber(xt)-1;
1930 for (Int_t nr= first; nr>=rf; nr-=step) {
1932 if (t.GetKinkIndexes()[0]>0){
1933 for (Int_t i=0;i<3;i++){
1934 Int_t index = t.GetKinkIndexes()[i];
1935 if (index==0) break;
1936 if (index<0) continue;
1938 AliKink * kink = (AliKink*)fEvent->GetKink(index-1);
1940 printf("PROBLEM\n");
1943 Int_t kinkrow = kink->GetTPCRow0()+2+Int_t(0.5/(0.05+kink->GetAngle(2)));
1945 AliExternalTrackParam paramd(t);
1946 kink->SetDaughter(paramd);
1947 kink->SetStatus(2,5);
1954 if (nr==80) t.UpdateReference();
1955 if (nr<fInnerSec->GetNRows())
1956 fSectors = fInnerSec;
1958 fSectors = fOuterSec;
1959 if (FollowToNext(t,nr)==0)
1972 Int_t AliTPCtrackerMI::FollowBackProlongation(AliTPCseed& t, Int_t rf) {
1973 //-----------------------------------------------------------------
1974 // This function tries to find a track prolongation.
1975 //-----------------------------------------------------------------
1977 Double_t xt=t.GetX();
1978 Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
1979 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
1980 if (alpha < 0. ) alpha += 2.*TMath::Pi();
1981 t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN);
1983 Int_t first = t.GetFirstPoint();
1984 if (first<GetRowNumber(xt)+1) first = GetRowNumber(xt)+1;
1986 if (first<0) first=0;
1987 for (Int_t nr=first; nr<=rf; nr++) {
1988 // if ( (TMath::Abs(t.GetSnp())>0.95)) break;//patch 28 fev 06
1989 if (t.GetKinkIndexes()[0]<0){
1990 for (Int_t i=0;i<3;i++){
1991 Int_t index = t.GetKinkIndexes()[i];
1992 if (index==0) break;
1993 if (index>0) continue;
1994 index = TMath::Abs(index);
1995 AliKink * kink = (AliKink*)fEvent->GetKink(index-1);
1997 printf("PROBLEM\n");
2000 Int_t kinkrow = kink->GetTPCRow0()-2-Int_t(0.5/(0.05+kink->GetAngle(2)));
2002 AliExternalTrackParam paramm(t);
2003 kink->SetMother(paramm);
2004 kink->SetStatus(2,1);
2011 if (nr<fInnerSec->GetNRows())
2012 fSectors = fInnerSec;
2014 fSectors = fOuterSec;
2025 Float_t AliTPCtrackerMI::OverlapFactor(AliTPCseed * s1, AliTPCseed * s2, Int_t &sum1, Int_t & sum2)
2033 Float_t dz2 =(s1->GetZ() - s2->GetZ());
2036 Float_t dy2 =TMath::Abs((s1->GetY() - s2->GetY()));
2038 Float_t distance = TMath::Sqrt(dz2+dy2);
2039 if (distance>4.) return 0; // if there are far away - not overlap - to reduce combinatorics
2042 Int_t firstpoint = TMath::Min(s1->GetFirstPoint(),s2->GetFirstPoint());
2043 Int_t lastpoint = TMath::Max(s1->GetLastPoint(),s2->GetLastPoint());
2048 if (firstpoint>lastpoint) {
2049 firstpoint =lastpoint;
2054 for (Int_t i=firstpoint-1;i<lastpoint+1;i++){
2055 if (s1->GetClusterIndex2(i)>0) sum1++;
2056 if (s2->GetClusterIndex2(i)>0) sum2++;
2057 if (s1->GetClusterIndex2(i)==s2->GetClusterIndex2(i) && s1->GetClusterIndex2(i)>0) {
2061 if (sum<5) return 0;
2063 Float_t summin = TMath::Min(sum1+1,sum2+1);
2064 Float_t ratio = (sum+1)/Float_t(summin);
2068 void AliTPCtrackerMI::SignShared(AliTPCseed * s1, AliTPCseed * s2)
2072 Float_t thetaCut = 0.2;//+10.*TMath::Sqrt(s1->GetSigmaTglZ()+ s2->GetSigmaTglZ());
2073 if (TMath::Abs(s1->GetTgl()-s2->GetTgl())>thetaCut) return;
2074 Float_t minCl = TMath::Min(s1->GetNumberOfClusters(),s2->GetNumberOfClusters());
2075 Int_t cutN0 = TMath::Max(5,TMath::Nint(0.1*minCl));
2080 //Int_t firstpoint = TMath::Max(s1->GetFirstPoint(),s2->GetFirstPoint());
2081 //Int_t lastpoint = TMath::Min(s1->GetLastPoint(),s2->GetLastPoint());
2082 Int_t firstpoint = 0;
2083 Int_t lastpoint = 160;
2085 // if (firstpoint>=lastpoint-5) return;;
2087 for (Int_t i=firstpoint;i<lastpoint;i++){
2088 // if ( (s1->GetClusterIndex2(i)&0xFFFF8FFF)==(s2->GetClusterIndex2(i)&0xFFFF8FFF) && s1->GetClusterIndex2(i)>0) {
2089 if ( (s1->GetClusterIndex2(i))==(s2->GetClusterIndex2(i)) && s1->GetClusterIndex2(i)>0) {
2091 s1->SetSharedMapBit(i, kTRUE);
2092 s2->SetSharedMapBit(i, kTRUE);
2094 if (s1->GetClusterIndex2(i)>0)
2095 s1->SetClusterMapBit(i, kTRUE);
2097 if (sumshared>cutN0){
2100 for (Int_t i=firstpoint;i<lastpoint;i++){
2101 // if ( (s1->GetClusterIndex2(i)&0xFFFF8FFF)==(s2->GetClusterIndex2(i)&0xFFFF8FFF) && s1->GetClusterIndex2(i)>0) {
2102 if ( (s1->GetClusterIndex2(i))==(s2->GetClusterIndex2(i)) && s1->GetClusterIndex2(i)>0) {
2103 AliTPCTrackerPoint *p1 = s1->GetTrackPoint(i);
2104 AliTPCTrackerPoint *p2 = s2->GetTrackPoint(i);;
2105 if (s1->IsActive()&&s2->IsActive()){
2106 p1->SetShared(kTRUE);
2107 p2->SetShared(kTRUE);
2113 if (sumshared>cutN0){
2114 for (Int_t i=0;i<4;i++){
2115 if (s1->GetOverlapLabel(3*i)==0){
2116 s1->SetOverlapLabel(3*i, s2->GetLabel());
2117 s1->SetOverlapLabel(3*i+1,sumshared);
2118 s1->SetOverlapLabel(3*i+2,s2->GetUniqueID());
2122 for (Int_t i=0;i<4;i++){
2123 if (s2->GetOverlapLabel(3*i)==0){
2124 s2->SetOverlapLabel(3*i, s1->GetLabel());
2125 s2->SetOverlapLabel(3*i+1,sumshared);
2126 s2->SetOverlapLabel(3*i+2,s1->GetUniqueID());
2133 void AliTPCtrackerMI::SignShared(TObjArray * arr)
2136 //sort trackss according sectors
2138 for (Int_t i=0; i<arr->GetEntriesFast(); i++) {
2139 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2141 //if (pt) RotateToLocal(pt);
2145 arr->Sort(); // sorting according relative sectors
2146 arr->Expand(arr->GetEntries());
2149 Int_t nseed=arr->GetEntriesFast();
2150 for (Int_t i=0; i<nseed; i++) {
2151 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2153 for (Int_t j=0;j<=12;j++){
2154 pt->SetOverlapLabel(j,0);
2157 for (Int_t i=0; i<nseed; i++) {
2158 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2160 if (pt->GetRemoval()>10) continue;
2161 for (Int_t j=i+1; j<nseed; j++){
2162 AliTPCseed *pt2=(AliTPCseed*)arr->UncheckedAt(j);
2163 if (TMath::Abs(pt->GetRelativeSector()-pt2->GetRelativeSector())>1) continue;
2165 if (pt2->GetRemoval()<=10) {
2166 //if ( TMath::Abs(pt->GetRelativeSector()-pt2->GetRelativeSector())>0) break;
2174 void AliTPCtrackerMI::SortTracks(TObjArray * arr, Int_t mode) const
2177 //sort tracks in array according mode criteria
2178 Int_t nseed = arr->GetEntriesFast();
2179 for (Int_t i=0; i<nseed; i++) {
2180 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2191 void AliTPCtrackerMI::RemoveUsed2(TObjArray * arr, Float_t factor1, Float_t factor2, Int_t minimal)
2194 // Loop over all tracks and remove overlaped tracks (with lower quality)
2196 // 1. Unsign clusters
2197 // 2. Sort tracks according quality
2198 // Quality is defined by the number of cluster between first and last points
2200 // 3. Loop over tracks - decreasing quality order
2201 // a.) remove - If the fraction of shared cluster less than factor (1- n or 2)
2202 // b.) remove - If the minimal number of clusters less than minimal and not ITS
2203 // c.) if track accepted - sign clusters
2205 //Called in - AliTPCtrackerMI::Clusters2Tracks()
2206 // - AliTPCtrackerMI::PropagateBack()
2207 // - AliTPCtrackerMI::RefitInward()
2210 // factor1 - factor for constrained
2211 // factor2 - for non constrained tracks
2212 // if (Float_t(shared+1)/Float_t(found+1)>factor) - DELETE
2216 Int_t nseed = arr->GetEntriesFast();
2217 Float_t * quality = new Float_t[nseed];
2218 Int_t * indexes = new Int_t[nseed];
2222 for (Int_t i=0; i<nseed; i++) {
2223 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2228 pt->UpdatePoints(); //select first last max dens points
2229 Float_t * points = pt->GetPoints();
2230 if (points[3]<0.8) quality[i] =-1;
2231 quality[i] = (points[2]-points[0])+pt->GetNumberOfClusters();
2232 //prefer high momenta tracks if overlaps
2233 quality[i] *= TMath::Sqrt(TMath::Abs(pt->Pt())+0.5);
2235 TMath::Sort(nseed,quality,indexes);
2238 for (Int_t itrack=0; itrack<nseed; itrack++) {
2239 Int_t trackindex = indexes[itrack];
2240 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(trackindex);
2243 if (quality[trackindex]<0){
2245 delete arr->RemoveAt(trackindex);
2248 arr->RemoveAt(trackindex);
2254 Int_t first = Int_t(pt->GetPoints()[0]);
2255 Int_t last = Int_t(pt->GetPoints()[2]);
2256 Double_t factor = (pt->GetBConstrain()) ? factor1: factor2;
2258 Int_t found,foundable,shared;
2259 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
2260 // pt->GetClusterStatistic(0,160, found, foundable,shared,kFALSE);
2261 Bool_t itsgold =kFALSE;
2264 if (pt->GetESD()->GetITSclusters(dummy)>4) itsgold= kTRUE;
2268 if (Float_t(shared+1)/Float_t(found+1)>factor){
2269 if (pt->GetKinkIndexes()[0]!=0) continue; //don't remove tracks - part of the kinks
2270 if( AliTPCReconstructor::StreamLevel()>3){
2271 TTreeSRedirector &cstream = *fDebugStreamer;
2272 cstream<<"RemoveUsed"<<
2273 "iter="<<fIteration<<
2277 delete arr->RemoveAt(trackindex);
2280 if (pt->GetNumberOfClusters()<50&&(found-0.5*shared)<minimal){ //remove short tracks
2281 if (pt->GetKinkIndexes()[0]!=0) continue; //don't remove tracks - part of the kinks
2282 if( AliTPCReconstructor::StreamLevel()>3){
2283 TTreeSRedirector &cstream = *fDebugStreamer;
2284 cstream<<"RemoveShort"<<
2285 "iter="<<fIteration<<
2289 delete arr->RemoveAt(trackindex);
2295 //if (sharedfactor>0.4) continue;
2296 if (pt->GetKinkIndexes()[0]>0) continue;
2297 //Remove tracks with undefined properties - seems
2298 if (pt->GetSigmaY2()<kAlmost0) continue; // ? what is the origin ?
2300 for (Int_t i=first; i<last; i++) {
2301 Int_t index=pt->GetClusterIndex2(i);
2302 // if (index<0 || index&0x8000 ) continue;
2303 if (index<0 || index&0x8000 ) continue;
2304 AliTPCclusterMI *c= pt->GetClusterPointer(i);
2311 Info("RemoveUsed","\n*****\nNumber of good tracks after shared removal\t%d\n",fNtracks);
2317 void AliTPCtrackerMI::DumpClusters(Int_t iter, TObjArray *trackArray)
2320 // Dump clusters after reco
2321 // signed and unsigned cluster can be visualized
2322 // 1. Unsign all cluster
2323 // 2. Sign all used clusters
2326 Int_t nseed = trackArray->GetEntries();
2327 for (Int_t i=0; i<nseed; i++){
2328 AliTPCseed *pt=(AliTPCseed*)trackArray->UncheckedAt(i);
2332 Bool_t isKink=pt->GetKinkIndex(0)!=0;
2333 for (Int_t j=0; j<160; ++j) {
2334 Int_t index=pt->GetClusterIndex2(j);
2335 if (index<0) continue;
2336 AliTPCclusterMI *c= pt->GetClusterPointer(j);
2338 if (isKink) c->Use(100); // kink
2339 c->Use(10); // by default usage 10
2344 for (Int_t sec=0;sec<fkNIS;sec++){
2345 for (Int_t row=0;row<fInnerSec->GetNRows();row++){
2346 AliTPCclusterMI *cl = fInnerSec[sec][row].GetClusters1();
2347 for (Int_t icl =0;icl< fInnerSec[sec][row].GetN1();icl++){
2348 Float_t gx[3]; cl[icl].GetGlobalXYZ(gx);
2349 (*fDebugStreamer)<<"clDump"<<
2357 cl = fInnerSec[sec][row].GetClusters2();
2358 for (Int_t icl =0;icl< fInnerSec[sec][row].GetN2();icl++){
2359 Float_t gx[3]; cl[icl].GetGlobalXYZ(gx);
2360 (*fDebugStreamer)<<"clDump"<<
2371 for (Int_t sec=0;sec<fkNOS;sec++){
2372 for (Int_t row=0;row<fOuterSec->GetNRows();row++){
2373 AliTPCclusterMI *cl = fOuterSec[sec][row].GetClusters1();
2374 for (Int_t icl =0;icl< fOuterSec[sec][row].GetN1();icl++){
2375 Float_t gx[3]; cl[icl].GetGlobalXYZ(gx);
2376 (*fDebugStreamer)<<"clDump"<<
2384 cl = fOuterSec[sec][row].GetClusters2();
2385 for (Int_t icl =0;icl< fOuterSec[sec][row].GetN2();icl++){
2386 Float_t gx[3]; cl[icl].GetGlobalXYZ(gx);
2387 (*fDebugStreamer)<<"clDump"<<
2399 void AliTPCtrackerMI::UnsignClusters()
2402 // loop over all clusters and unsign them
2405 for (Int_t sec=0;sec<fkNIS;sec++){
2406 for (Int_t row=0;row<fInnerSec->GetNRows();row++){
2407 AliTPCclusterMI *cl = fInnerSec[sec][row].GetClusters1();
2408 for (Int_t icl =0;icl< fInnerSec[sec][row].GetN1();icl++)
2409 // if (cl[icl].IsUsed(10))
2411 cl = fInnerSec[sec][row].GetClusters2();
2412 for (Int_t icl =0;icl< fInnerSec[sec][row].GetN2();icl++)
2413 //if (cl[icl].IsUsed(10))
2418 for (Int_t sec=0;sec<fkNOS;sec++){
2419 for (Int_t row=0;row<fOuterSec->GetNRows();row++){
2420 AliTPCclusterMI *cl = fOuterSec[sec][row].GetClusters1();
2421 for (Int_t icl =0;icl< fOuterSec[sec][row].GetN1();icl++)
2422 //if (cl[icl].IsUsed(10))
2424 cl = fOuterSec[sec][row].GetClusters2();
2425 for (Int_t icl =0;icl< fOuterSec[sec][row].GetN2();icl++)
2426 //if (cl[icl].IsUsed(10))
2435 void AliTPCtrackerMI::SignClusters(const TObjArray * arr, Float_t fnumber, Float_t fdensity)
2438 //sign clusters to be "used"
2440 // snumber and sdensity sign number of sigmas - bellow mean value to be accepted
2441 // loop over "primaries"
2455 Int_t nseed = arr->GetEntriesFast();
2456 for (Int_t i=0; i<nseed; i++) {
2457 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2461 if (!(pt->IsActive())) continue;
2462 Float_t dens = pt->GetNumberOfClusters()/Float_t(pt->GetNFoundable());
2463 if ( (dens>0.7) && (pt->GetNumberOfClusters()>70)){
2465 sumdens2+= dens*dens;
2466 sumn += pt->GetNumberOfClusters();
2467 sumn2 += pt->GetNumberOfClusters()*pt->GetNumberOfClusters();
2468 Float_t chi2 = pt->GetChi2()/pt->GetNumberOfClusters();
2471 sumchi2 +=chi2*chi2;
2476 Float_t mdensity = 0.9;
2477 Float_t meann = 130;
2478 Float_t meanchi = 1;
2479 Float_t sdensity = 0.1;
2480 Float_t smeann = 10;
2481 Float_t smeanchi =0.4;
2485 mdensity = sumdens/sum;
2487 meanchi = sumchi/sum;
2489 sdensity = sumdens2/sum-mdensity*mdensity;
2491 sdensity = TMath::Sqrt(sdensity);
2495 smeann = sumn2/sum-meann*meann;
2497 smeann = TMath::Sqrt(smeann);
2501 smeanchi = sumchi2/sum - meanchi*meanchi;
2503 smeanchi = TMath::Sqrt(smeanchi);
2509 //REMOVE SHORT DELTAS or tracks going out of sensitive volume of TPC
2511 for (Int_t i=0; i<nseed; i++) {
2512 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2516 if (pt->GetBSigned()) continue;
2517 if (pt->GetBConstrain()) continue;
2518 //if (!(pt->IsActive())) continue;
2520 Int_t found,foundable,shared;
2521 pt->GetClusterStatistic(0,160,found, foundable,shared);
2522 if (shared/float(found)>0.3) {
2523 if (shared/float(found)>0.9 ){
2524 //delete arr->RemoveAt(i);
2529 Bool_t isok =kFALSE;
2530 if ( (pt->GetNShared()/pt->GetNumberOfClusters()<0.5) &&pt->GetNumberOfClusters()>60)
2532 if ((TMath::Abs(1/pt->GetC())<100.) && (pt->GetNShared()/pt->GetNumberOfClusters()<0.7))
2534 if (TMath::Abs(pt->GetZ()/pt->GetX())>1.1)
2536 if ( (TMath::Abs(pt->GetSnp()>0.7) && pt->GetD(0,0)>60.))
2540 for (Int_t j=0; j<160; ++j) {
2541 Int_t index=pt->GetClusterIndex2(j);
2542 if (index<0) continue;
2543 AliTPCclusterMI *c= pt->GetClusterPointer(j);
2545 //if (!(c->IsUsed(10))) c->Use();
2552 Double_t maxchi = meanchi+2.*smeanchi;
2554 for (Int_t i=0; i<nseed; i++) {
2555 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2559 //if (!(pt->IsActive())) continue;
2560 if (pt->GetBSigned()) continue;
2561 Double_t chi = pt->GetChi2()/pt->GetNumberOfClusters();
2562 if (chi>maxchi) continue;
2565 Float_t dens = pt->GetNumberOfClusters()/Float_t(pt->GetNFoundable());
2567 //sign only tracks with enoug big density at the beginning
2569 if ((pt->GetDensityFirst(40)<0.75) && pt->GetNumberOfClusters()<meann) continue;
2572 Double_t mindens = TMath::Max(double(mdensity-sdensity*fdensity*bfactor),0.65);
2573 Double_t minn = TMath::Max(Int_t(meann-fnumber*smeann*bfactor),50);
2575 // if (pt->fBConstrain) mindens = TMath::Max(mdensity-sdensity*fdensity*bfactor,0.65);
2576 if ( (pt->GetRemoval()==10) && (pt->GetSnp()>0.8)&&(dens>mindens))
2579 if ((dens>mindens && pt->GetNumberOfClusters()>minn) && chi<maxchi ){
2580 //Int_t noc=pt->GetNumberOfClusters();
2581 pt->SetBSigned(kTRUE);
2582 for (Int_t j=0; j<160; ++j) {
2584 Int_t index=pt->GetClusterIndex2(j);
2585 if (index<0) continue;
2586 AliTPCclusterMI *c= pt->GetClusterPointer(j);
2588 // if (!(c->IsUsed(10))) c->Use();
2593 // gLastCheck = nseed;
2601 void AliTPCtrackerMI::StopNotActive(const TObjArray * arr, Int_t row0, Float_t th0, Float_t th1, Float_t th2) const
2603 // stop not active tracks
2604 // take th1 as threshold for number of founded to number of foundable on last 10 active rows
2605 // take th2 as threshold for number of founded to number of foundable on last 20 active rows
2606 Int_t nseed = arr->GetEntriesFast();
2608 for (Int_t i=0; i<nseed; i++) {
2609 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2613 if (!(pt->IsActive())) continue;
2614 StopNotActive(pt,row0,th0, th1,th2);
2620 void AliTPCtrackerMI::StopNotActive(AliTPCseed * seed, Int_t row0, Float_t th0, Float_t th1,
2623 // stop not active tracks
2624 // take th1 as threshold for number of founded to number of foundable on last 10 active rows
2625 // take th2 as threshold for number of founded to number of foundable on last 20 active rows
2628 Int_t foundable = 0;
2629 Int_t maxindex = seed->GetLastPoint(); //last foundable row
2630 if (seed->GetNFoundable()*th0 > seed->GetNumberOfClusters()) {
2631 seed->Desactivate(10) ;
2635 for (Int_t i=row0; i<maxindex; i++){
2636 Int_t index = seed->GetClusterIndex2(i);
2637 if (index!=-1) foundable++;
2639 if (foundable<=30) sumgood1++;
2640 if (foundable<=50) {
2647 if (foundable>=30.){
2648 if (sumgood1<(th1*30.)) seed->Desactivate(10);
2651 if (sumgood2<(th2*50.)) seed->Desactivate(10);
2655 Int_t AliTPCtrackerMI::RefitInward(AliESDEvent *event)
2658 // back propagation of ESD tracks
2661 const Int_t kMaxFriendTracks=2000;
2665 //PrepareForProlongation(fSeeds,1);
2666 PropagateForward2(fSeeds);
2667 RemoveUsed2(fSeeds,0.4,0.4,20);
2669 TObjArray arraySeed(fSeeds->GetEntries());
2670 for (Int_t i=0;i<fSeeds->GetEntries();i++) {
2671 arraySeed.AddAt(fSeeds->At(i),i);
2673 SignShared(&arraySeed);
2674 // FindCurling(fSeeds, event,2); // find multi found tracks
2675 FindSplitted(fSeeds, event,2); // find multi found tracks
2676 if (AliTPCReconstructor::StreamLevel()>5) FindMultiMC(fSeeds, fEvent,2); // find multi found tracks
2679 Int_t nseed = fSeeds->GetEntriesFast();
2680 for (Int_t i=0;i<nseed;i++){
2681 AliTPCseed * seed = (AliTPCseed*) fSeeds->UncheckedAt(i);
2682 if (!seed) continue;
2683 if (seed->GetKinkIndex(0)>0) UpdateKinkQualityD(seed); // update quality informations for kinks
2684 AliESDtrack *esd=event->GetTrack(i);
2686 if (seed->GetNumberOfClusters()<60 && seed->GetNumberOfClusters()<(esd->GetTPCclusters(0) -5)*0.8){
2687 AliExternalTrackParam paramIn;
2688 AliExternalTrackParam paramOut;
2689 Int_t ncl = seed->RefitTrack(seed,¶mIn,¶mOut);
2690 if (AliTPCReconstructor::StreamLevel()>2) {
2691 (*fDebugStreamer)<<"RecoverIn"<<
2695 "pout.="<<¶mOut<<
2700 seed->Set(paramIn.GetX(),paramIn.GetAlpha(),paramIn.GetParameter(),paramIn.GetCovariance());
2701 seed->SetNumberOfClusters(ncl);
2705 seed->PropagateTo(fkParam->GetInnerRadiusLow());
2706 seed->UpdatePoints();
2707 AddCovariance(seed);
2709 seed->CookdEdx(0.02,0.6);
2710 CookLabel(seed,0.1); //For comparison only
2711 esd->SetTPCClusterMap(seed->GetClusterMap());
2712 esd->SetTPCSharedMap(seed->GetSharedMap());
2714 if (AliTPCReconstructor::StreamLevel()>1 && seed!=0&&esd!=0) {
2715 TTreeSRedirector &cstream = *fDebugStreamer;
2722 if (seed->GetNumberOfClusters()>15){
2723 esd->UpdateTrackParams(seed,AliESDtrack::kTPCrefit);
2724 esd->SetTPCPoints(seed->GetPoints());
2725 esd->SetTPCPointsF(seed->GetNFoundable());
2726 Int_t ndedx = seed->GetNCDEDX(0);
2727 Float_t sdedx = seed->GetSDEDX(0);
2728 Float_t dedx = seed->GetdEdx();
2729 esd->SetTPCsignal(dedx, sdedx, ndedx);
2731 // add seed to the esd track in Calib level
2733 Bool_t storeFriend = gRandom->Rndm()<(kMaxFriendTracks)/Float_t(nseed);
2734 if (AliTPCReconstructor::StreamLevel()>0 &&storeFriend){
2735 AliTPCseed * seedCopy = new AliTPCseed(*seed, kTRUE);
2736 esd->AddCalibObject(seedCopy);
2741 //printf("problem\n");
2744 //FindKinks(fSeeds,event);
2745 if (AliTPCReconstructor::StreamLevel()>3) DumpClusters(2,fSeeds);
2746 Info("RefitInward","Number of refitted tracks %d",ntracks);
2751 Int_t AliTPCtrackerMI::PropagateBack(AliESDEvent *event)
2754 // back propagation of ESD tracks
2760 PropagateBack(fSeeds);
2761 RemoveUsed2(fSeeds,0.4,0.4,20);
2762 //FindCurling(fSeeds, fEvent,1);
2763 FindSplitted(fSeeds, event,1); // find multi found tracks
2764 if (AliTPCReconstructor::StreamLevel()>5) FindMultiMC(fSeeds, fEvent,1); // find multi found tracks
2767 Int_t nseed = fSeeds->GetEntriesFast();
2769 for (Int_t i=0;i<nseed;i++){
2770 AliTPCseed * seed = (AliTPCseed*) fSeeds->UncheckedAt(i);
2771 if (!seed) continue;
2772 if (seed->GetKinkIndex(0)<0) UpdateKinkQualityM(seed); // update quality informations for kinks
2773 seed->UpdatePoints();
2774 AddCovariance(seed);
2775 AliESDtrack *esd=event->GetTrack(i);
2776 if (seed->GetNumberOfClusters()<60 && seed->GetNumberOfClusters()<(esd->GetTPCclusters(0) -5)*0.8){
2777 AliExternalTrackParam paramIn;
2778 AliExternalTrackParam paramOut;
2779 Int_t ncl = seed->RefitTrack(seed,¶mIn,¶mOut);
2780 if (AliTPCReconstructor::StreamLevel()>2) {
2781 (*fDebugStreamer)<<"RecoverBack"<<
2785 "pout.="<<¶mOut<<
2790 seed->Set(paramOut.GetX(),paramOut.GetAlpha(),paramOut.GetParameter(),paramOut.GetCovariance());
2791 seed->SetNumberOfClusters(ncl);
2794 seed->CookdEdx(0.02,0.6);
2795 CookLabel(seed,0.1); //For comparison only
2796 if (seed->GetNumberOfClusters()>15){
2797 esd->UpdateTrackParams(seed,AliESDtrack::kTPCout);
2798 esd->SetTPCPoints(seed->GetPoints());
2799 esd->SetTPCPointsF(seed->GetNFoundable());
2800 Int_t ndedx = seed->GetNCDEDX(0);
2801 Float_t sdedx = seed->GetSDEDX(0);
2802 Float_t dedx = seed->GetdEdx();
2803 esd->SetTPCsignal(dedx, sdedx, ndedx);
2805 Int_t eventnumber = event->GetEventNumberInFile();// patch 28 fev 06
2806 // This is most likely NOT the event number you'd like to use. It has nothing to do with the 'real' event number
2807 if (AliTPCReconstructor::StreamLevel()>1 && esd) {
2808 (*fDebugStreamer)<<"Cback"<<
2811 "EventNrInFile="<<eventnumber<<
2816 if (AliTPCReconstructor::StreamLevel()>3) DumpClusters(1,fSeeds);
2817 //FindKinks(fSeeds,event);
2818 Info("PropagateBack","Number of back propagated tracks %d",ntracks);
2825 void AliTPCtrackerMI::DeleteSeeds()
2830 Int_t nseed = fSeeds->GetEntriesFast();
2831 for (Int_t i=0;i<nseed;i++){
2832 AliTPCseed * seed = (AliTPCseed*)fSeeds->At(i);
2833 if (seed) delete fSeeds->RemoveAt(i);
2840 void AliTPCtrackerMI::ReadSeeds(AliESDEvent *const event, Int_t direction)
2843 //read seeds from the event
2845 Int_t nentr=event->GetNumberOfTracks();
2847 Info("ReadSeeds", "Number of ESD tracks: %d\n", nentr);
2852 fSeeds = new TObjArray(nentr);
2856 for (Int_t i=0; i<nentr; i++) {
2857 AliESDtrack *esd=event->GetTrack(i);
2858 ULong_t status=esd->GetStatus();
2859 if (!(status&AliESDtrack::kTPCin)) continue;
2860 AliTPCtrack t(*esd);
2861 t.SetNumberOfClusters(0);
2862 // AliTPCseed *seed = new AliTPCseed(t,t.GetAlpha());
2863 AliTPCseed *seed = new AliTPCseed(t/*,t.GetAlpha()*/);
2864 seed->SetUniqueID(esd->GetID());
2865 AddCovariance(seed); //add systematic ucertainty
2866 for (Int_t ikink=0;ikink<3;ikink++) {
2867 Int_t index = esd->GetKinkIndex(ikink);
2868 seed->GetKinkIndexes()[ikink] = index;
2869 if (index==0) continue;
2870 index = TMath::Abs(index);
2871 AliESDkink * kink = fEvent->GetKink(index-1);
2872 if (kink&&esd->GetKinkIndex(ikink)<0){
2873 if ((status & AliESDtrack::kTRDrefit) != 0) kink->SetStatus(1,2);
2874 if ((status & AliESDtrack::kITSout) != 0) kink->SetStatus(1,0);
2876 if (kink&&esd->GetKinkIndex(ikink)>0){
2877 if ((status & AliESDtrack::kTRDrefit) != 0) kink->SetStatus(1,6);
2878 if ((status & AliESDtrack::kITSout) != 0) kink->SetStatus(1,4);
2882 if (((status&AliESDtrack::kITSout)==0)&&(direction==1)) seed->ResetCovariance(10.);
2883 if ( direction ==2 &&(status & AliESDtrack::kTRDrefit) == 0 ) seed->ResetCovariance(10.);
2884 //if ( direction ==2 && ((status & AliESDtrack::kTPCout) == 0) ) {
2885 // fSeeds->AddAt(0,i);
2889 if ( direction ==2 &&(status & AliESDtrack::kTRDrefit) > 0 ) {
2890 Double_t par0[5],par1[5],alpha,x;
2891 esd->GetInnerExternalParameters(alpha,x,par0);
2892 esd->GetExternalParameters(x,par1);
2893 Double_t delta1 = TMath::Abs(par0[4]-par1[4])/(0.000000001+TMath::Abs(par0[4]+par1[4]));
2894 Double_t delta2 = TMath::Abs(par0[3]-par1[3]);
2896 if (esd->GetTRDncls()>0) trdchi2 = esd->GetTRDchi2()/esd->GetTRDncls();
2897 //reset covariance if suspicious
2898 if ( (delta1>0.1) || (delta2>0.006) ||trdchi2>7.)
2899 seed->ResetCovariance(10.);
2904 // rotate to the local coordinate system
2906 fSectors=fInnerSec; fN=fkNIS;
2907 Double_t alpha=seed->GetAlpha() - fSectors->GetAlphaShift();
2908 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
2909 if (alpha < 0. ) alpha += 2.*TMath::Pi();
2910 Int_t ns=Int_t(alpha/fSectors->GetAlpha())%fN;
2911 alpha =ns*fSectors->GetAlpha() + fSectors->GetAlphaShift();
2912 if (alpha<-TMath::Pi()) alpha += 2*TMath::Pi();
2913 if (alpha>=TMath::Pi()) alpha -= 2*TMath::Pi();
2914 alpha-=seed->GetAlpha();
2915 if (!seed->Rotate(alpha)) {
2921 if (esd->GetKinkIndex(0)<=0){
2922 for (Int_t irow=0;irow<160;irow++){
2923 Int_t index = seed->GetClusterIndex2(irow);
2926 AliTPCclusterMI * cl = GetClusterMI(index);
2927 seed->SetClusterPointer(irow,cl);
2929 if ((index & 0x8000)==0){
2930 cl->Use(10); // accepted cluster
2932 cl->Use(6); // close cluster not accepted
2935 Info("ReadSeeds","Not found cluster");
2940 fSeeds->AddAt(seed,i);
2946 //_____________________________________________________________________________
2947 void AliTPCtrackerMI::MakeSeeds3(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, Float_t cuts[4],
2948 Float_t deltay, Int_t ddsec) {
2949 //-----------------------------------------------------------------
2950 // This function creates track seeds.
2951 // SEEDING WITH VERTEX CONSTRAIN
2952 //-----------------------------------------------------------------
2953 // cuts[0] - fP4 cut
2954 // cuts[1] - tan(phi) cut
2955 // cuts[2] - zvertex cut
2956 // cuts[3] - fP3 cut
2964 Double_t x[5], c[15];
2965 // Int_t di = i1-i2;
2967 AliTPCseed * seed = new AliTPCseed();
2968 Double_t alpha=fSectors->GetAlpha(), shift=fSectors->GetAlphaShift();
2969 Double_t cs=cos(alpha), sn=sin(alpha);
2971 // Double_t x1 =fOuterSec->GetX(i1);
2972 //Double_t xx2=fOuterSec->GetX(i2);
2974 Double_t x1 =GetXrow(i1);
2975 Double_t xx2=GetXrow(i2);
2977 Double_t x3=GetX(), y3=GetY(), z3=GetZ();
2979 Int_t imiddle = (i2+i1)/2; //middle pad row index
2980 Double_t xm = GetXrow(imiddle); // radius of middle pad-row
2981 const AliTPCtrackerRow& krm=GetRow(sec,imiddle); //middle pad -row
2985 const AliTPCtrackerRow& kr1=GetRow(ns,i1);
2986 Double_t ymax = GetMaxY(i1)-kr1.GetDeadZone()-1.5;
2987 Double_t ymaxm = GetMaxY(imiddle)-kr1.GetDeadZone()-1.5;
2990 // change cut on curvature if it can't reach this layer
2991 // maximal curvature set to reach it
2992 Double_t dvertexmax = TMath::Sqrt((x1-x3)*(x1-x3)+(ymax+5-y3)*(ymax+5-y3));
2993 if (dvertexmax*0.5*cuts[0]>0.85){
2994 cuts[0] = 0.85/(dvertexmax*0.5+1.);
2996 Double_t r2min = 1/(cuts[0]*cuts[0]); //minimal square of radius given by cut
2999 if (deltay>0) ddsec = 0;
3000 // loop over clusters
3001 for (Int_t is=0; is < kr1; is++) {
3003 if (kr1[is]->IsUsed(10)) continue;
3004 Double_t y1=kr1[is]->GetY(), z1=kr1[is]->GetZ();
3005 //if (TMath::Abs(y1)>ymax) continue;
3007 if (deltay>0 && TMath::Abs(ymax-TMath::Abs(y1))> deltay ) continue; // seed only at the edge
3009 // find possible directions
3010 Float_t anglez = (z1-z3)/(x1-x3);
3011 Float_t extraz = z1 - anglez*(x1-xx2); // extrapolated z
3014 //find rotation angles relative to line given by vertex and point 1
3015 Double_t dvertex2 = (x1-x3)*(x1-x3)+(y1-y3)*(y1-y3);
3016 Double_t dvertex = TMath::Sqrt(dvertex2);
3017 Double_t angle13 = TMath::ATan((y1-y3)/(x1-x3));
3018 Double_t cs13 = cos(-angle13), sn13 = sin(-angle13);
3021 // loop over 2 sectors
3027 Double_t dddz1=0; // direction of delta inclination in z axis
3034 for (Int_t dsec = dsec1; dsec<=dsec2;dsec++){
3035 Int_t sec2 = sec + dsec;
3037 // AliTPCtrackerRow& kr2 = fOuterSec[(sec2+fkNOS)%fkNOS][i2];
3038 //AliTPCtrackerRow& kr2m = fOuterSec[(sec2+fkNOS)%fkNOS][imiddle];
3039 AliTPCtrackerRow& kr2 = GetRow((sec2+fkNOS)%fkNOS,i2);
3040 AliTPCtrackerRow& kr2m = GetRow((sec2+fkNOS)%fkNOS,imiddle);
3041 Int_t index1 = TMath::Max(kr2.Find(extraz-0.6-dddz1*TMath::Abs(z1)*0.05)-1,0);
3042 Int_t index2 = TMath::Min(kr2.Find(extraz+0.6+dddz2*TMath::Abs(z1)*0.05)+1,kr2);
3044 // rotation angles to p1-p3
3045 Double_t cs13r = cos(-angle13+dsec*alpha)/dvertex, sn13r = sin(-angle13+dsec*alpha)/dvertex;
3046 Double_t x2, y2, z2;
3048 // Double_t dymax = maxangle*TMath::Abs(x1-xx2);
3051 Double_t dxx0 = (xx2-x3)*cs13r;
3052 Double_t dyy0 = (xx2-x3)*sn13r;
3053 for (Int_t js=index1; js < index2; js++) {
3054 const AliTPCclusterMI *kcl = kr2[js];
3055 if (kcl->IsUsed(10)) continue;
3057 //calcutate parameters
3059 Double_t yy0 = dyy0 +(kcl->GetY()-y3)*cs13r;
3061 if (TMath::Abs(yy0)<0.000001) continue;
3062 Double_t xx0 = dxx0 -(kcl->GetY()-y3)*sn13r;
3063 Double_t y0 = 0.5*(xx0*xx0+yy0*yy0-xx0)/yy0;
3064 Double_t r02 = (0.25+y0*y0)*dvertex2;
3065 //curvature (radius) cut
3066 if (r02<r2min) continue;
3070 Double_t c0 = 1/TMath::Sqrt(r02);
3074 //Double_t dfi0 = 2.*TMath::ASin(dvertex*c0*0.5);
3075 //Double_t dfi1 = 2.*TMath::ASin(TMath::Sqrt(yy0*yy0+(1-xx0)*(1-xx0))*dvertex*c0*0.5);
3076 Double_t dfi0 = 2.*AliTPCFastMath::FastAsin(dvertex*c0*0.5);
3077 Double_t dfi1 = 2.*AliTPCFastMath::FastAsin(TMath::Sqrt(yy0*yy0+(1-xx0)*(1-xx0))*dvertex*c0*0.5);
3080 Double_t z0 = kcl->GetZ();
3081 Double_t zzzz2 = z1-(z1-z3)*dfi1/dfi0;
3082 if (TMath::Abs(zzzz2-z0)>0.5) continue;
3085 Double_t dip = (z1-z0)*c0/dfi1;
3086 Double_t x0 = (0.5*cs13+y0*sn13)*dvertex*c0;
3097 x2= xx2*cs-y2*sn*dsec;
3098 y2=+xx2*sn*dsec+y2*cs;
3108 // do we have cluster at the middle ?
3110 GetProlongation(x1,xm,x,ym,zm);
3112 AliTPCclusterMI * cm=0;
3113 if (TMath::Abs(ym)-ymaxm<0){
3114 cm = krm.FindNearest2(ym,zm,1.0,0.6,dummy);
3115 if ((!cm) || (cm->IsUsed(10))) {
3120 // rotate y1 to system 0
3121 // get state vector in rotated system
3122 Double_t yr1 = (-0.5*sn13+y0*cs13)*dvertex*c0;
3123 Double_t xr2 = x0*cs+yr1*sn*dsec;
3124 Double_t xr[5]={kcl->GetY(),kcl->GetZ(), xr2, dip, c0};
3126 GetProlongation(xx2,xm,xr,ym,zm);
3127 if (TMath::Abs(ym)-ymaxm<0){
3128 cm = kr2m.FindNearest2(ym,zm,1.0,0.6,dummy);
3129 if ((!cm) || (cm->IsUsed(10))) {
3139 dym = ym - cm->GetY();
3140 dzm = zm - cm->GetZ();
3147 Double_t sy1=kr1[is]->GetSigmaY2()*2., sz1=kr1[is]->GetSigmaZ2()*2.;
3148 Double_t sy2=kcl->GetSigmaY2()*2., sz2=kcl->GetSigmaZ2()*2.;
3149 //Double_t sy3=400*3./12., sy=0.1, sz=0.1;
3150 Double_t sy3=25000*x[4]*x[4]+0.1, sy=0.1, sz=0.1;
3151 //Double_t sy3=25000*x[4]*x[4]*60+0.5, sy=0.1, sz=0.1;
3153 Double_t f40=(F1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
3154 Double_t f42=(F1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
3155 Double_t f43=(F1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
3156 Double_t f20=(F2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
3157 Double_t f22=(F2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
3158 Double_t f23=(F2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
3160 Double_t f30=(F3(x1,y1+sy,x2,y2,z1,z2)-x[3])/sy;
3161 Double_t f31=(F3(x1,y1,x2,y2,z1+sz,z2)-x[3])/sz;
3162 Double_t f32=(F3(x1,y1,x2,y2+sy,z1,z2)-x[3])/sy;
3163 Double_t f34=(F3(x1,y1,x2,y2,z1,z2+sz)-x[3])/sz;
3167 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3168 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3169 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3170 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3171 c[13]=f30*sy1*f40+f32*sy2*f42;
3172 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3174 // if (!BuildSeed(kr1[is],kcl,0,x1,x2,x3,x,c)) continue;
3176 UInt_t index=kr1.GetIndex(is);
3177 seed->~AliTPCseed(); // this does not set the pointer to 0...
3178 AliTPCseed *track=new(seed) AliTPCseed(x1, ns*alpha+shift, x, c, index);
3180 track->SetIsSeeding(kTRUE);
3181 track->SetSeed1(i1);
3182 track->SetSeed2(i2);
3183 track->SetSeedType(3);
3187 FollowProlongation(*track, (i1+i2)/2,1);
3188 Int_t foundable,found,shared;
3189 track->GetClusterStatistic((i1+i2)/2,i1, found, foundable, shared, kTRUE);
3190 if ((found<0.55*foundable) || shared>0.5*found || (track->GetSigmaY2()+track->GetSigmaZ2())>0.5){
3192 seed->~AliTPCseed();
3198 FollowProlongation(*track, i2,1);
3202 track->SetBConstrain(1);
3203 // track->fLastPoint = i1+fInnerSec->GetNRows(); // first cluster in track position
3204 track->SetLastPoint(i1); // first cluster in track position
3205 track->SetFirstPoint(track->GetLastPoint());
3207 if (track->GetNumberOfClusters()<(i1-i2)*0.5 ||
3208 track->GetNumberOfClusters() < track->GetNFoundable()*0.6 ||
3209 track->GetNShared()>0.4*track->GetNumberOfClusters() ) {
3211 seed->~AliTPCseed();
3215 // Z VERTEX CONDITION
3216 Double_t zv, bz=GetBz();
3217 if ( !track->GetZAt(0.,bz,zv) ) continue;
3218 if (TMath::Abs(zv-z3)>cuts[2]) {
3219 FollowProlongation(*track, TMath::Max(i2-20,0));
3220 if ( !track->GetZAt(0.,bz,zv) ) continue;
3221 if (TMath::Abs(zv-z3)>cuts[2]){
3222 FollowProlongation(*track, TMath::Max(i2-40,0));
3223 if ( !track->GetZAt(0.,bz,zv) ) continue;
3224 if (TMath::Abs(zv-z3)>cuts[2] &&(track->GetNumberOfClusters() > track->GetNFoundable()*0.7)){
3225 // make seed without constrain
3226 AliTPCseed * track2 = MakeSeed(track,0.2,0.5,1.);
3227 FollowProlongation(*track2, i2,1);
3228 track2->SetBConstrain(kFALSE);
3229 track2->SetSeedType(1);
3230 arr->AddLast(track2);
3232 seed->~AliTPCseed();
3237 seed->~AliTPCseed();
3244 track->SetSeedType(0);
3245 arr->AddLast(track);
3246 seed = new AliTPCseed;
3248 // don't consider other combinations
3249 if (track->GetNumberOfClusters() > track->GetNFoundable()*0.8)
3255 Info("MakeSeeds3","\nSeeding statistic:\t%d\t%d\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin,nout1,nout2);
3261 void AliTPCtrackerMI::MakeSeeds5(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, Float_t cuts[4],
3266 //-----------------------------------------------------------------
3267 // This function creates track seeds.
3268 //-----------------------------------------------------------------
3269 // cuts[0] - fP4 cut
3270 // cuts[1] - tan(phi) cut
3271 // cuts[2] - zvertex cut
3272 // cuts[3] - fP3 cut
3282 Double_t x[5], c[15];
3284 // make temporary seed
3285 AliTPCseed * seed = new AliTPCseed;
3286 Double_t alpha=fOuterSec->GetAlpha(), shift=fOuterSec->GetAlphaShift();
3287 // Double_t cs=cos(alpha), sn=sin(alpha);
3292 Double_t x1 = GetXrow(i1-1);
3293 const AliTPCtrackerRow& kr1=GetRow(sec,i1-1);
3294 Double_t y1max = GetMaxY(i1-1)-kr1.GetDeadZone()-1.5;
3296 Double_t x1p = GetXrow(i1);
3297 const AliTPCtrackerRow& kr1p=GetRow(sec,i1);
3299 Double_t x1m = GetXrow(i1-2);
3300 const AliTPCtrackerRow& kr1m=GetRow(sec,i1-2);
3303 //last 3 padrow for seeding
3304 AliTPCtrackerRow& kr3 = GetRow((sec+fkNOS)%fkNOS,i1-7);
3305 Double_t x3 = GetXrow(i1-7);
3306 // Double_t y3max= GetMaxY(i1-7)-kr3.fDeadZone-1.5;
3308 AliTPCtrackerRow& kr3p = GetRow((sec+fkNOS)%fkNOS,i1-6);
3309 Double_t x3p = GetXrow(i1-6);
3311 AliTPCtrackerRow& kr3m = GetRow((sec+fkNOS)%fkNOS,i1-8);
3312 Double_t x3m = GetXrow(i1-8);
3317 Int_t im = i1-4; //middle pad row index
3318 Double_t xm = GetXrow(im); // radius of middle pad-row
3319 const AliTPCtrackerRow& krm=GetRow(sec,im); //middle pad -row
3320 // Double_t ymmax = GetMaxY(im)-kr1.fDeadZone-1.5;
3323 Double_t deltax = x1-x3;
3324 Double_t dymax = deltax*cuts[1];
3325 Double_t dzmax = deltax*cuts[3];
3327 // loop over clusters
3328 for (Int_t is=0; is < kr1; is++) {
3330 if (kr1[is]->IsUsed(10)) continue;
3331 Double_t y1=kr1[is]->GetY(), z1=kr1[is]->GetZ();
3333 if (deltay>0 && TMath::Abs(y1max-TMath::Abs(y1))> deltay ) continue; // seed only at the edge
3335 Int_t index1 = TMath::Max(kr3.Find(z1-dzmax)-1,0);
3336 Int_t index2 = TMath::Min(kr3.Find(z1+dzmax)+1,kr3);
3342 for (Int_t js=index1; js < index2; js++) {
3343 const AliTPCclusterMI *kcl = kr3[js];
3344 if (kcl->IsUsed(10)) continue;
3346 // apply angular cuts
3347 if (TMath::Abs(y1-y3)>dymax) continue;
3350 if (TMath::Abs(z1-z3)>dzmax) continue;
3352 Double_t angley = (y1-y3)/(x1-x3);
3353 Double_t anglez = (z1-z3)/(x1-x3);
3355 Double_t erry = TMath::Abs(angley)*(x1-x1m)*0.5+0.5;
3356 Double_t errz = TMath::Abs(anglez)*(x1-x1m)*0.5+0.5;
3358 Double_t yyym = angley*(xm-x1)+y1;
3359 Double_t zzzm = anglez*(xm-x1)+z1;
3361 const AliTPCclusterMI *kcm = krm.FindNearest2(yyym,zzzm,erry,errz,index);
3363 if (kcm->IsUsed(10)) continue;
3365 erry = TMath::Abs(angley)*(x1-x1m)*0.4+0.5;
3366 errz = TMath::Abs(anglez)*(x1-x1m)*0.4+0.5;
3373 // look around first
3374 const AliTPCclusterMI *kc1m = kr1m.FindNearest2(angley*(x1m-x1)+y1,
3380 if (kc1m->IsUsed(10)) used++;
3382 const AliTPCclusterMI *kc1p = kr1p.FindNearest2(angley*(x1p-x1)+y1,
3388 if (kc1p->IsUsed(10)) used++;
3390 if (used>1) continue;
3391 if (found<1) continue;
3395 const AliTPCclusterMI *kc3m = kr3m.FindNearest2(angley*(x3m-x3)+y3,
3401 if (kc3m->IsUsed(10)) used++;
3405 const AliTPCclusterMI *kc3p = kr3p.FindNearest2(angley*(x3p-x3)+y3,
3411 if (kc3p->IsUsed(10)) used++;
3415 if (used>1) continue;
3416 if (found<3) continue;
3426 x[4]=F1(x1,y1,x2,y2,x3,y3);
3427 //if (TMath::Abs(x[4]) >= cuts[0]) continue;
3430 x[2]=F2(x1,y1,x2,y2,x3,y3);
3433 x[3]=F3n(x1,y1,x2,y2,z1,z2,x[4]);
3434 //if (TMath::Abs(x[3]) > cuts[3]) continue;
3438 Double_t sy1=0.1, sz1=0.1;
3439 Double_t sy2=0.1, sz2=0.1;
3440 Double_t sy3=0.1, sy=0.1, sz=0.1;
3442 Double_t f40=(F1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
3443 Double_t f42=(F1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
3444 Double_t f43=(F1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
3445 Double_t f20=(F2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
3446 Double_t f22=(F2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
3447 Double_t f23=(F2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
3449 Double_t f30=(F3(x1,y1+sy,x2,y2,z1,z2)-x[3])/sy;
3450 Double_t f31=(F3(x1,y1,x2,y2,z1+sz,z2)-x[3])/sz;
3451 Double_t f32=(F3(x1,y1,x2,y2+sy,z1,z2)-x[3])/sy;
3452 Double_t f34=(F3(x1,y1,x2,y2,z1,z2+sz)-x[3])/sz;
3456 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3457 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3458 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3459 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3460 c[13]=f30*sy1*f40+f32*sy2*f42;
3461 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3463 // if (!BuildSeed(kr1[is],kcl,0,x1,x2,x3,x,c)) continue;
3465 index=kr1.GetIndex(is);
3466 seed->~AliTPCseed();
3467 AliTPCseed *track=new(seed) AliTPCseed(x1, sec*alpha+shift, x, c, index);
3469 track->SetIsSeeding(kTRUE);
3472 FollowProlongation(*track, i1-7,1);
3473 if (track->GetNumberOfClusters() < track->GetNFoundable()*0.75 ||
3474 track->GetNShared()>0.6*track->GetNumberOfClusters() || ( track->GetSigmaY2()+ track->GetSigmaZ2())>0.6){
3476 seed->~AliTPCseed();
3482 FollowProlongation(*track, i2,1);
3483 track->SetBConstrain(0);
3484 track->SetLastPoint(i1+fInnerSec->GetNRows()); // first cluster in track position
3485 track->SetFirstPoint(track->GetLastPoint());
3487 if (track->GetNumberOfClusters()<(i1-i2)*0.5 ||
3488 track->GetNumberOfClusters()<track->GetNFoundable()*0.7 ||
3489 track->GetNShared()>2. || track->GetChi2()/track->GetNumberOfClusters()>6 || ( track->GetSigmaY2()+ track->GetSigmaZ2())>0.5 ) {
3491 seed->~AliTPCseed();
3496 FollowProlongation(*track, TMath::Max(i2-10,0),1);
3497 AliTPCseed * track2 = MakeSeed(track,0.2,0.5,0.9);
3498 FollowProlongation(*track2, i2,1);
3499 track2->SetBConstrain(kFALSE);
3500 track2->SetSeedType(4);
3501 arr->AddLast(track2);
3503 seed->~AliTPCseed();
3507 //arr->AddLast(track);
3508 //seed = new AliTPCseed;
3514 Info("MakeSeeds5","\nSeeding statiistic:\t%d\t%d\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin,nout1,nout2,nout3);
3520 //_____________________________________________________________________________
3521 void AliTPCtrackerMI::MakeSeeds2(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, Float_t */*cuts[4]*/,
3522 Float_t deltay, Bool_t /*bconstrain*/) {
3523 //-----------------------------------------------------------------
3524 // This function creates track seeds - without vertex constraint
3525 //-----------------------------------------------------------------
3526 // cuts[0] - fP4 cut - not applied
3527 // cuts[1] - tan(phi) cut
3528 // cuts[2] - zvertex cut - not applied
3529 // cuts[3] - fP3 cut
3539 Double_t alpha=fOuterSec->GetAlpha(), shift=fOuterSec->GetAlphaShift();
3540 // Double_t cs=cos(alpha), sn=sin(alpha);
3541 Int_t row0 = (i1+i2)/2;
3542 Int_t drow = (i1-i2)/2;
3543 const AliTPCtrackerRow& kr0=fSectors[sec][row0];
3544 AliTPCtrackerRow * kr=0;
3546 AliTPCpolyTrack polytrack;
3547 Int_t nclusters=fSectors[sec][row0];
3548 AliTPCseed * seed = new AliTPCseed;
3553 for (Int_t is=0; is < nclusters; is++) { //LOOP over clusters
3555 Int_t nfoundable =0;
3556 for (Int_t iter =1; iter<2; iter++){ //iterations
3557 const AliTPCtrackerRow& krm=fSectors[sec][row0-iter];
3558 const AliTPCtrackerRow& krp=fSectors[sec][row0+iter];
3559 const AliTPCclusterMI * cl= kr0[is];
3561 if (cl->IsUsed(10)) {
3567 Double_t x = kr0.GetX();
3568 // Initialization of the polytrack
3573 Double_t y0= cl->GetY();
3574 Double_t z0= cl->GetZ();
3578 Double_t ymax = fSectors->GetMaxY(row0)-kr0.GetDeadZone()-1.5;
3579 if (deltay>0 && TMath::Abs(ymax-TMath::Abs(y0))> deltay ) continue; // seed only at the edge
3581 erry = (0.5)*cl->GetSigmaY2()/TMath::Sqrt(cl->GetQ())*6;
3582 errz = (0.5)*cl->GetSigmaZ2()/TMath::Sqrt(cl->GetQ())*6;
3583 polytrack.AddPoint(x,y0,z0,erry, errz);
3586 if (cl->IsUsed(10)) sumused++;
3589 Float_t roady = (5*TMath::Sqrt(cl->GetSigmaY2()+0.2)+1.)*iter;
3590 Float_t roadz = (5*TMath::Sqrt(cl->GetSigmaZ2()+0.2)+1.)*iter;
3593 AliTPCclusterMI * cl1 = krm.FindNearest(y0,z0,roady,roadz);
3594 if (cl1 && TMath::Abs(ymax-TMath::Abs(y0))) {
3595 erry = (0.5)*cl1->GetSigmaY2()/TMath::Sqrt(cl1->GetQ())*3;
3596 errz = (0.5)*cl1->GetSigmaZ2()/TMath::Sqrt(cl1->GetQ())*3;
3597 if (cl1->IsUsed(10)) sumused++;
3598 polytrack.AddPoint(x,cl1->GetY(),cl1->GetZ(),erry,errz);
3602 AliTPCclusterMI * cl2 = krp.FindNearest(y0,z0,roady,roadz);
3604 erry = (0.5)*cl2->GetSigmaY2()/TMath::Sqrt(cl2->GetQ())*3;
3605 errz = (0.5)*cl2->GetSigmaZ2()/TMath::Sqrt(cl2->GetQ())*3;
3606 if (cl2->IsUsed(10)) sumused++;
3607 polytrack.AddPoint(x,cl2->GetY(),cl2->GetZ(),erry,errz);
3610 if (sumused>0) continue;
3612 polytrack.UpdateParameters();
3618 nfoundable = polytrack.GetN();
3619 nfound = nfoundable;
3621 for (Int_t ddrow = iter+1; ddrow<drow;ddrow++){
3622 Float_t maxdist = 0.8*(1.+3./(ddrow));
3623 for (Int_t delta = -1;delta<=1;delta+=2){
3624 Int_t row = row0+ddrow*delta;
3625 kr = &(fSectors[sec][row]);
3626 Double_t xn = kr->GetX();
3627 Double_t ymax1 = fSectors->GetMaxY(row)-kr->GetDeadZone()-1.5;
3628 polytrack.GetFitPoint(xn,yn,zn);
3629 if (TMath::Abs(yn)>ymax1) continue;
3631 AliTPCclusterMI * cln = kr->FindNearest(yn,zn,roady,roadz);
3633 Float_t dist = TMath::Sqrt( (yn-cln->GetY())*(yn-cln->GetY())+(zn-cln->GetZ())*(zn-cln->GetZ()));
3636 erry = (dist+0.3)*cln->GetSigmaY2()/TMath::Sqrt(cln->GetQ())*(1.+1./(ddrow));
3637 errz = (dist+0.3)*cln->GetSigmaZ2()/TMath::Sqrt(cln->GetQ())*(1.+1./(ddrow));
3638 if (cln->IsUsed(10)) {
3639 // printf("used\n");
3647 polytrack.AddPoint(xn,cln->GetY(),cln->GetZ(),erry, errz);
3652 if ( (sumused>3) || (sumused>0.5*nfound) || (nfound<0.6*nfoundable)) break;
3653 polytrack.UpdateParameters();
3656 if ( (sumused>3) || (sumused>0.5*nfound)) {
3657 //printf("sumused %d\n",sumused);
3662 polytrack.GetFitDerivation(kr0.GetX(),dy,dz);
3663 AliTPCpolyTrack track2;
3665 polytrack.Refit(track2,0.5+TMath::Abs(dy)*0.3,0.4+TMath::Abs(dz)*0.3);
3666 if (track2.GetN()<0.5*nfoundable) continue;
3669 if ((nfound>0.6*nfoundable) &&( nfoundable>0.4*(i1-i2))) {
3671 // test seed with and without constrain
3672 for (Int_t constrain=0; constrain<=0;constrain++){
3673 // add polytrack candidate
3675 Double_t x[5], c[15];
3676 Double_t x1,x2,x3,y1,y2,y3,z1,z2,z3;
3677 track2.GetBoundaries(x3,x1);
3679 track2.GetFitPoint(x1,y1,z1);
3680 track2.GetFitPoint(x2,y2,z2);
3681 track2.GetFitPoint(x3,y3,z3);
3683 //is track pointing to the vertex ?
3686 polytrack.GetFitPoint(x0,y0,z0);
3699 x[4]=F1(x1,y1,x2,y2,x3,y3);
3701 // if (TMath::Abs(x[4]) >= cuts[0]) continue; //
3702 x[2]=F2(x1,y1,x2,y2,x3,y3);
3704 //if (TMath::Abs(x[4]*x1-x[2]) >= cuts[1]) continue;
3705 //x[3]=F3(x1,y1,x2,y2,z1,z2);
3706 x[3]=F3n(x1,y1,x3,y3,z1,z3,x[4]);
3707 //if (TMath::Abs(x[3]) > cuts[3]) continue;
3710 Double_t sy =0.1, sz =0.1;
3711 Double_t sy1=0.02, sz1=0.02;
3712 Double_t sy2=0.02, sz2=0.02;
3716 sy3=25000*x[4]*x[4]+0.1, sy=0.1, sz=0.1;
3719 Double_t f40=(F1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
3720 Double_t f42=(F1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
3721 Double_t f43=(F1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
3722 Double_t f20=(F2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
3723 Double_t f22=(F2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
3724 Double_t f23=(F2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
3726 Double_t f30=(F3(x1,y1+sy,x3,y3,z1,z3)-x[3])/sy;
3727 Double_t f31=(F3(x1,y1,x3,y3,z1+sz,z3)-x[3])/sz;
3728 Double_t f32=(F3(x1,y1,x3,y3+sy,z1,z3)-x[3])/sy;
3729 Double_t f34=(F3(x1,y1,x3,y3,z1,z3+sz)-x[3])/sz;
3734 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3735 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3736 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3737 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3738 c[13]=f30*sy1*f40+f32*sy2*f42;
3739 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3741 //Int_t row1 = fSectors->GetRowNumber(x1);
3742 Int_t row1 = GetRowNumber(x1);
3746 seed->~AliTPCseed();
3747 AliTPCseed *track=new(seed) AliTPCseed(x1,sec*alpha+shift,x,c,index);
3748 track->SetIsSeeding(kTRUE);
3749 Int_t rc=FollowProlongation(*track, i2);
3750 if (constrain) track->SetBConstrain(1);
3752 track->SetBConstrain(0);
3753 track->SetLastPoint(row1+fInnerSec->GetNRows()); // first cluster in track position
3754 track->SetFirstPoint(track->GetLastPoint());
3756 if (rc==0 || track->GetNumberOfClusters()<(i1-i2)*0.5 ||
3757 track->GetNumberOfClusters() < track->GetNFoundable()*0.6 ||
3758 track->GetNShared()>0.4*track->GetNumberOfClusters()) {
3761 seed->~AliTPCseed();
3764 arr->AddLast(track);
3765 seed = new AliTPCseed;
3769 } // if accepted seed
3772 Info("MakeSeeds2","\nSeeding statiistic:\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin3);
3778 AliTPCseed *AliTPCtrackerMI::MakeSeed(AliTPCseed *const track, Float_t r0, Float_t r1, Float_t r2)
3782 //reseed using track points
3783 Int_t p0 = int(r0*track->GetNumberOfClusters()); // point 0
3784 Int_t p1 = int(r1*track->GetNumberOfClusters());
3785 Int_t p2 = int(r2*track->GetNumberOfClusters()); // last point
3787 Double_t x0[3],x1[3],x2[3];
3788 for (Int_t i=0;i<3;i++){
3794 // find track position at given ratio of the length
3795 Int_t sec0=0, sec1=0, sec2=0;
3798 for (Int_t i=0;i<160;i++){
3799 if (track->GetClusterPointer(i)){
3801 AliTPCTrackerPoint *trpoint =track->GetTrackPoint(i);
3802 if ( (index<p0) || x0[0]<0 ){
3803 if (trpoint->GetX()>1){
3804 clindex = track->GetClusterIndex2(i);
3806 x0[0] = trpoint->GetX();
3807 x0[1] = trpoint->GetY();
3808 x0[2] = trpoint->GetZ();
3809 sec0 = ((clindex&0xff000000)>>24)%18;
3814 if ( (index<p1) &&(trpoint->GetX()>1)){
3815 clindex = track->GetClusterIndex2(i);
3817 x1[0] = trpoint->GetX();
3818 x1[1] = trpoint->GetY();
3819 x1[2] = trpoint->GetZ();
3820 sec1 = ((clindex&0xff000000)>>24)%18;
3823 if ( (index<p2) &&(trpoint->GetX()>1)){
3824 clindex = track->GetClusterIndex2(i);
3826 x2[0] = trpoint->GetX();
3827 x2[1] = trpoint->GetY();
3828 x2[2] = trpoint->GetZ();
3829 sec2 = ((clindex&0xff000000)>>24)%18;
3836 Double_t alpha, cs,sn, xx2,yy2;
3838 alpha = (sec1-sec2)*fSectors->GetAlpha();
3839 cs = TMath::Cos(alpha);
3840 sn = TMath::Sin(alpha);
3841 xx2= x1[0]*cs-x1[1]*sn;
3842 yy2= x1[0]*sn+x1[1]*cs;
3846 alpha = (sec0-sec2)*fSectors->GetAlpha();
3847 cs = TMath::Cos(alpha);
3848 sn = TMath::Sin(alpha);
3849 xx2= x0[0]*cs-x0[1]*sn;
3850 yy2= x0[0]*sn+x0[1]*cs;
3856 Double_t x[5],c[15];
3860 x[4]=F1(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]);
3861 // if (x[4]>1) return 0;
3862 x[2]=F2(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]);
3863 x[3]=F3n(x2[0],x2[1],x0[0],x0[1],x2[2],x0[2],x[4]);
3864 //if (TMath::Abs(x[3]) > 2.2) return 0;
3865 //if (TMath::Abs(x[2]) > 1.99) return 0;
3867 Double_t sy =0.1, sz =0.1;
3869 Double_t sy1=0.02+track->GetSigmaY2(), sz1=0.02+track->GetSigmaZ2();
3870 Double_t sy2=0.01+track->GetSigmaY2(), sz2=0.01+track->GetSigmaZ2();
3871 Double_t sy3=0.01+track->GetSigmaY2();
3873 Double_t f40=(F1(x2[0],x2[1]+sy,x1[0],x1[1],x0[0],x0[1])-x[4])/sy;
3874 Double_t f42=(F1(x2[0],x2[1],x1[0],x1[1]+sy,x0[0],x0[1])-x[4])/sy;
3875 Double_t f43=(F1(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]+sy)-x[4])/sy;
3876 Double_t f20=(F2(x2[0],x2[1]+sy,x1[0],x1[1],x0[0],x0[1])-x[2])/sy;
3877 Double_t f22=(F2(x2[0],x2[1],x1[0],x1[1]+sy,x0[0],x0[1])-x[2])/sy;
3878 Double_t f23=(F2(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]+sy)-x[2])/sy;
3880 Double_t f30=(F3(x2[0],x2[1]+sy,x0[0],x0[1],x2[2],x0[2])-x[3])/sy;
3881 Double_t f31=(F3(x2[0],x2[1],x0[0],x0[1],x2[2]+sz,x0[2])-x[3])/sz;
3882 Double_t f32=(F3(x2[0],x2[1],x0[0],x0[1]+sy,x2[2],x0[2])-x[3])/sy;
3883 Double_t f34=(F3(x2[0],x2[1],x0[0],x0[1],x2[2],x0[2]+sz)-x[3])/sz;
3888 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3889 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3890 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3891 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3892 c[13]=f30*sy1*f40+f32*sy2*f42;
3893 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3895 // Int_t row1 = fSectors->GetRowNumber(x2[0]);
3896 AliTPCseed *seed=new AliTPCseed(x2[0], sec2*fSectors->GetAlpha()+fSectors->GetAlphaShift(), x, c, 0);
3897 // Double_t y0,z0,y1,z1, y2,z2;
3898 //seed->GetProlongation(x0[0],y0,z0);
3899 // seed->GetProlongation(x1[0],y1,z1);
3900 //seed->GetProlongation(x2[0],y2,z2);
3902 seed->SetLastPoint(pp2);
3903 seed->SetFirstPoint(pp2);
3910 AliTPCseed *AliTPCtrackerMI::ReSeed(const AliTPCseed *track, Float_t r0, Float_t r1, Float_t r2)
3914 //reseed using founded clusters
3916 // Find the number of clusters
3917 Int_t nclusters = 0;
3918 for (Int_t irow=0;irow<160;irow++){
3919 if (track->GetClusterIndex(irow)>0) nclusters++;
3923 ipos[0] = TMath::Max(int(r0*nclusters),0); // point 0 cluster
3924 ipos[1] = TMath::Min(int(r1*nclusters),nclusters-1); //
3925 ipos[2] = TMath::Min(int(r2*nclusters),nclusters-1); // last point
3929 Int_t row[3],sec[3]={0,0,0};
3931 // find track row position at given ratio of the length
3933 for (Int_t irow=0;irow<160;irow++){
3934 if (track->GetClusterIndex2(irow)<0) continue;
3936 for (Int_t ipoint=0;ipoint<3;ipoint++){
3937 if (index<=ipos[ipoint]) row[ipoint] = irow;
3941 //Get cluster and sector position
3942 for (Int_t ipoint=0;ipoint<3;ipoint++){
3943 Int_t clindex = track->GetClusterIndex2(row[ipoint]);
3944 AliTPCclusterMI * cl = GetClusterMI(clindex);
3947 // AliTPCclusterMI * cl = GetClusterMI(clindex);
3950 sec[ipoint] = ((clindex&0xff000000)>>24)%18;
3951 xyz[ipoint][0] = GetXrow(row[ipoint]);
3952 xyz[ipoint][1] = cl->GetY();
3953 xyz[ipoint][2] = cl->GetZ();
3957 // Calculate seed state vector and covariance matrix
3959 Double_t alpha, cs,sn, xx2,yy2;
3961 alpha = (sec[1]-sec[2])*fSectors->GetAlpha();
3962 cs = TMath::Cos(alpha);
3963 sn = TMath::Sin(alpha);
3964 xx2= xyz[1][0]*cs-xyz[1][1]*sn;
3965 yy2= xyz[1][0]*sn+xyz[1][1]*cs;
3969 alpha = (sec[0]-sec[2])*fSectors->GetAlpha();
3970 cs = TMath::Cos(alpha);
3971 sn = TMath::Sin(alpha);
3972 xx2= xyz[0][0]*cs-xyz[0][1]*sn;
3973 yy2= xyz[0][0]*sn+xyz[0][1]*cs;
3979 Double_t x[5],c[15];
3983 x[4]=F1(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
3984 x[2]=F2(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
3985 x[3]=F3n(xyz[2][0],xyz[2][1],xyz[0][0],xyz[0][1],xyz[2][2],xyz[0][2],x[4]);
3987 Double_t sy =0.1, sz =0.1;
3989 Double_t sy1=0.2, sz1=0.2;
3990 Double_t sy2=0.2, sz2=0.2;
3993 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;
3994 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;
3995 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;
3996 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;
3997 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;
3998 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;
4000 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;
4001 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;
4002 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;
4003 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;
4008 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
4009 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
4010 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
4011 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
4012 c[13]=f30*sy1*f40+f32*sy2*f42;
4013 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
4015 // Int_t row1 = fSectors->GetRowNumber(xyz[2][0]);
4016 AliTPCseed *seed=new AliTPCseed(xyz[2][0], sec[2]*fSectors->GetAlpha()+fSectors->GetAlphaShift(), x, c, 0);
4017 seed->SetLastPoint(row[2]);
4018 seed->SetFirstPoint(row[2]);
4023 AliTPCseed *AliTPCtrackerMI::ReSeed(AliTPCseed *track,Int_t r0, Bool_t forward)
4027 //reseed using founded clusters
4030 Int_t row[3]={0,0,0};
4031 Int_t sec[3]={0,0,0};
4033 // forward direction
4035 for (Int_t irow=r0;irow<160;irow++){
4036 if (track->GetClusterIndex(irow)>0){
4041 for (Int_t irow=160;irow>r0;irow--){
4042 if (track->GetClusterIndex(irow)>0){
4047 for (Int_t irow=row[2]-15;irow>row[0];irow--){
4048 if (track->GetClusterIndex(irow)>0){
4056 for (Int_t irow=0;irow<r0;irow++){
4057 if (track->GetClusterIndex(irow)>0){
4062 for (Int_t irow=r0;irow>0;irow--){
4063 if (track->GetClusterIndex(irow)>0){
4068 for (Int_t irow=row[2]-15;irow>row[0];irow--){
4069 if (track->GetClusterIndex(irow)>0){
4076 if ((row[2]-row[0])<20) return 0;
4077 if (row[1]==0) return 0;
4080 //Get cluster and sector position
4081 for (Int_t ipoint=0;ipoint<3;ipoint++){
4082 Int_t clindex = track->GetClusterIndex2(row[ipoint]);
4083 AliTPCclusterMI * cl = GetClusterMI(clindex);
4086 // AliTPCclusterMI * cl = GetClusterMI(clindex);
4089 sec[ipoint] = ((clindex&0xff000000)>>24)%18;
4090 xyz[ipoint][0] = GetXrow(row[ipoint]);
4091 AliTPCTrackerPoint * point = track->GetTrackPoint(row[ipoint]);
4092 if (point&&ipoint<2){
4094 xyz[ipoint][1] = point->GetY();
4095 xyz[ipoint][2] = point->GetZ();
4098 xyz[ipoint][1] = cl->GetY();
4099 xyz[ipoint][2] = cl->GetZ();
4106 // Calculate seed state vector and covariance matrix
4108 Double_t alpha, cs,sn, xx2,yy2;
4110 alpha = (sec[1]-sec[2])*fSectors->GetAlpha();
4111 cs = TMath::Cos(alpha);
4112 sn = TMath::Sin(alpha);
4113 xx2= xyz[1][0]*cs-xyz[1][1]*sn;
4114 yy2= xyz[1][0]*sn+xyz[1][1]*cs;
4118 alpha = (sec[0]-sec[2])*fSectors->GetAlpha();
4119 cs = TMath::Cos(alpha);
4120 sn = TMath::Sin(alpha);
4121 xx2= xyz[0][0]*cs-xyz[0][1]*sn;
4122 yy2= xyz[0][0]*sn+xyz[0][1]*cs;
4128 Double_t x[5],c[15];
4132 x[4]=F1(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
4133 x[2]=F2(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
4134 x[3]=F3n(xyz[2][0],xyz[2][1],xyz[0][0],xyz[0][1],xyz[2][2],xyz[0][2],x[4]);
4136 Double_t sy =0.1, sz =0.1;
4138 Double_t sy1=0.2, sz1=0.2;
4139 Double_t sy2=0.2, sz2=0.2;
4142 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;
4143 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;
4144 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;
4145 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;
4146 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;
4147 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;
4149 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;
4150 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;
4151 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;
4152 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;
4157 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
4158 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
4159 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
4160 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
4161 c[13]=f30*sy1*f40+f32*sy2*f42;
4162 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
4164 // Int_t row1 = fSectors->GetRowNumber(xyz[2][0]);
4165 AliTPCseed *seed=new AliTPCseed(xyz[2][0], sec[2]*fSectors->GetAlpha()+fSectors->GetAlphaShift(), x, c, 0);
4166 seed->SetLastPoint(row[2]);
4167 seed->SetFirstPoint(row[2]);
4168 for (Int_t i=row[0];i<row[2];i++){
4169 seed->SetClusterIndex(i, track->GetClusterIndex(i));
4177 void AliTPCtrackerMI::FindMultiMC(const TObjArray * array, AliESDEvent */*esd*/, Int_t iter)
4180 // find multi tracks - THIS FUNCTION IS ONLY FOR DEBUG PURPOSES
4182 // Use AliTPCReconstructor::StreamLevel()>2 if you want to tune parameters - cuts
4184 // Two reasons to have multiple find tracks
4185 // 1. Curling tracks can be find more than once
4186 // 2. Splitted tracks
4187 // a.) Multiple seeding to increase tracking efficiency - (~ 100% reached)
4188 // b.) Edge effect on the sector boundaries
4191 // Algorithm done in 2 phases - because of CPU consumption
4192 // it is n^2 algorithm - for lead-lead 20000x20000 combination are investigated
4194 // Algorihm for curling tracks sign:
4195 // 1 phase -makes a very rough fast cuts to minimize combinatorics
4196 // a.) opposite sign
4197 // b.) one of the tracks - not pointing to the primary vertex -
4198 // c.) delta tan(theta)
4200 // 2 phase - calculates DCA between tracks - time consument
4205 // General cuts - for splitted tracks and for curling tracks
4207 const Float_t kMaxdPhi = 0.2; // maximal distance in phi
4209 // Curling tracks cuts
4214 Int_t nentries = array->GetEntriesFast();
4215 AliHelix *helixes = new AliHelix[nentries];
4216 Float_t *xm = new Float_t[nentries];
4217 Float_t *dz0 = new Float_t[nentries];
4218 Float_t *dz1 = new Float_t[nentries];
4224 // Find track COG in x direction - point with best defined parameters
4226 for (Int_t i=0;i<nentries;i++){
4227 AliTPCseed* track = (AliTPCseed*)array->At(i);
4228 if (!track) continue;
4229 track->SetCircular(0);
4230 new (&helixes[i]) AliHelix(*track);
4234 track->GetDZ(GetX(),GetY(),GetZ(),GetBz(),dz);
4237 for (Int_t icl=0; icl<160; icl++){
4238 AliTPCclusterMI * cl = track->GetClusterPointer(icl);
4244 if (ncl>0) xm[i]/=Float_t(ncl);
4247 for (Int_t i0=0;i0<nentries;i0++){
4248 AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
4249 if (!track0) continue;
4250 Float_t xc0 = helixes[i0].GetHelix(6);
4251 Float_t yc0 = helixes[i0].GetHelix(7);
4252 Float_t r0 = helixes[i0].GetHelix(8);
4253 Float_t rc0 = TMath::Sqrt(xc0*xc0+yc0*yc0);
4254 Float_t fi0 = TMath::ATan2(yc0,xc0);
4256 for (Int_t i1=i0+1;i1<nentries;i1++){
4257 AliTPCseed * track1 = (AliTPCseed*)array->At(i1);
4258 if (!track1) continue;
4259 Int_t lab0=track0->GetLabel();
4260 Int_t lab1=track1->GetLabel();
4261 if (TMath::Abs(lab0)!=TMath::Abs(lab1)) continue;
4263 Float_t xc1 = helixes[i1].GetHelix(6);
4264 Float_t yc1 = helixes[i1].GetHelix(7);
4265 Float_t r1 = helixes[i1].GetHelix(8);
4266 Float_t rc1 = TMath::Sqrt(xc1*xc1+yc1*yc1);
4267 Float_t fi1 = TMath::ATan2(yc1,xc1);
4269 Float_t dfi = fi0-fi1;
4272 if (dfi>1.5*TMath::Pi()) dfi-=TMath::Pi(); // take care about edge effect
4273 if (dfi<-1.5*TMath::Pi()) dfi+=TMath::Pi(); //
4274 if (TMath::Abs(dfi)>kMaxdPhi&&helixes[i0].GetHelix(4)*helixes[i1].GetHelix(4)<0){
4276 // if short tracks with undefined sign
4277 fi1 = -TMath::ATan2(yc1,-xc1);
4280 Float_t dtheta = TMath::Abs(track0->GetTgl()-track1->GetTgl())<TMath::Abs(track0->GetTgl()+track1->GetTgl())? track0->GetTgl()-track1->GetTgl():track0->GetTgl()+track1->GetTgl();
4283 // debug stream to tune "fast cuts"
4285 Double_t dist[3]; // distance at X
4286 Double_t mdist[3]={0,0,0}; // mean distance X+-40cm
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])+40.,dist,AliTracker::GetBz());
4290 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[i]);
4291 track0->GetDistance(track1,0.5*(xm[i0]+xm[i1]),dist,AliTracker::GetBz());
4292 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[i]);
4293 for (Int_t i=0;i<3;i++) mdist[i]*=0.33333;
4297 for (Int_t icl=0; icl<160; icl++){
4298 AliTPCclusterMI * cl0 = track0->GetClusterPointer(icl);
4299 AliTPCclusterMI * cl1 = track1->GetClusterPointer(icl);
4302 if (cl0==cl1) sums++;
4306 if (AliTPCReconstructor::StreamLevel()>5) {
4307 TTreeSRedirector &cstream = *fDebugStreamer;
4312 "Tr0.="<<track0<< // seed0
4313 "Tr1.="<<track1<< // seed1
4314 "h0.="<<&helixes[i0]<<
4315 "h1.="<<&helixes[i1]<<
4317 "sum="<<sum<< //the sum of rows with cl in both
4318 "sums="<<sums<< //the sum of shared clusters
4319 "xm0="<<xm[i0]<< // the center of track
4320 "xm1="<<xm[i1]<< // the x center of track
4321 // General cut variables
4322 "dfi="<<dfi<< // distance in fi angle
4323 "dtheta="<<dtheta<< // distance int theta angle
4329 "dist0="<<dist[0]<< //distance x
4330 "dist1="<<dist[1]<< //distance y
4331 "dist2="<<dist[2]<< //distance z
4332 "mdist0="<<mdist[0]<< //distance x
4333 "mdist1="<<mdist[1]<< //distance y
4334 "mdist2="<<mdist[2]<< //distance z
4348 if (AliTPCReconstructor::StreamLevel()>1) {
4349 AliInfo("Time for curling tracks removal DEBUGGING MC");
4356 void AliTPCtrackerMI::FindSplitted(TObjArray * array, AliESDEvent */*esd*/, Int_t /*iter*/){
4358 // Find Splitted tracks and remove the one with worst quality
4359 // Corresponding debug streamer to tune selections - "Splitted2"
4361 // 0. Sort tracks according quility
4362 // 1. Propagate the tracks to the reference radius
4363 // 2. Double_t loop to select close tracks (only to speed up process)
4364 // 3. Calculate cluster overlap ratio - and remove the track if bigger than a threshold
4365 // 4. Delete temporary parameters
4367 const Double_t xref=GetXrow(63); // reference radius -IROC/OROC boundary
4369 const Double_t kCutP1=10; // delta Z cut 10 cm
4370 const Double_t kCutP2=0.15; // delta snp(fi) cut 0.15
4371 const Double_t kCutP3=0.15; // delta tgl(theta) cut 0.15
4372 const Double_t kCutAlpha=0.15; // delta alpha cut
4373 Int_t firstpoint = 0;
4374 Int_t lastpoint = 160;
4376 Int_t nentries = array->GetEntriesFast();
4377 AliExternalTrackParam *params = new AliExternalTrackParam[nentries];
4383 //0. Sort tracks according quality
4384 //1. Propagate the ext. param to reference radius
4385 Int_t nseed = array->GetEntriesFast();
4386 if (nseed<=0) return;
4387 Float_t * quality = new Float_t[nseed];
4388 Int_t * indexes = new Int_t[nseed];
4389 for (Int_t i=0; i<nseed; i++) {
4390 AliTPCseed *pt=(AliTPCseed*)array->UncheckedAt(i);
4395 pt->UpdatePoints(); //select first last max dens points
4396 Float_t * points = pt->GetPoints();
4397 if (points[3]<0.8) quality[i] =-1;
4398 quality[i] = (points[2]-points[0])+pt->GetNumberOfClusters();
4399 //prefer high momenta tracks if overlaps
4400 quality[i] *= TMath::Sqrt(TMath::Abs(pt->Pt())+0.5);
4402 AliTracker::PropagateTrackToBxByBz(&(params[i]),xref,pt->GetMass(),5.,kTRUE);
4403 AliTracker::PropagateTrackToBxByBz(&(params[i]),xref,pt->GetMass(),1.,kTRUE);
4405 TMath::Sort(nseed,quality,indexes);
4407 // 3. Loop over pair of tracks
4409 for (Int_t i0=0; i0<nseed; i0++) {
4410 Int_t index0=indexes[i0];
4411 if (!(array->UncheckedAt(index0))) continue;
4412 AliTPCseed *s1 = (AliTPCseed*)array->UncheckedAt(index0);
4413 if (!s1->IsActive()) continue;
4414 AliExternalTrackParam &par0=params[index0];
4415 for (Int_t i1=i0+1; i1<nseed; i1++) {
4416 Int_t index1=indexes[i1];
4417 if (!(array->UncheckedAt(index1))) continue;
4418 AliTPCseed *s2 = (AliTPCseed*)array->UncheckedAt(index1);
4419 if (!s2->IsActive()) continue;
4420 if (s2->GetKinkIndexes()[0]!=0)
4421 if (s2->GetKinkIndexes()[0] == -s1->GetKinkIndexes()[0]) continue;
4422 AliExternalTrackParam &par1=params[index1];
4423 if (TMath::Abs(par0.GetParameter()[3]-par1.GetParameter()[3])>kCutP3) continue;
4424 if (TMath::Abs(par0.GetParameter()[1]-par1.GetParameter()[1])>kCutP1) continue;
4425 if (TMath::Abs(par0.GetParameter()[2]-par1.GetParameter()[2])>kCutP2) continue;
4426 Double_t dAlpha= TMath::Abs(par0.GetAlpha()-par1.GetAlpha());
4427 if (dAlpha>TMath::Pi()) dAlpha-=TMath::Pi();
4428 if (TMath::Abs(dAlpha)>kCutAlpha) continue;
4433 Int_t firstShared=lastpoint, lastShared=firstpoint;
4434 Int_t firstRow=lastpoint, lastRow=firstpoint;
4436 for (Int_t i=firstpoint;i<lastpoint;i++){
4437 if (s1->GetClusterIndex2(i)>0) nall0++;
4438 if (s2->GetClusterIndex2(i)>0) nall1++;
4439 if (s1->GetClusterIndex2(i)>0 && s2->GetClusterIndex2(i)>0) {
4440 if (i<firstRow) firstRow=i;
4441 if (i>lastRow) lastRow=i;
4443 if ( (s1->GetClusterIndex2(i))==(s2->GetClusterIndex2(i)) && s1->GetClusterIndex2(i)>0) {
4444 if (i<firstShared) firstShared=i;
4445 if (i>lastShared) lastShared=i;
4449 Double_t ratio0 = Float_t(sumShared)/Float_t(TMath::Min(nall0+1,nall1+1));
4450 Double_t ratio1 = Float_t(sumShared)/Float_t(TMath::Max(nall0+1,nall1+1));
4452 if( AliTPCReconstructor::StreamLevel()>1){
4453 TTreeSRedirector &cstream = *fDebugStreamer;
4454 Int_t n0=s1->GetNumberOfClusters();
4455 Int_t n1=s2->GetNumberOfClusters();
4456 Int_t n0F=s1->GetNFoundable();
4457 Int_t n1F=s2->GetNFoundable();
4458 Int_t lab0=s1->GetLabel();
4459 Int_t lab1=s2->GetLabel();
4461 cstream<<"Splitted2"<<
4462 "iter="<<fIteration<<
4463 "lab0="<<lab0<< // MC label if exist
4464 "lab1="<<lab1<< // MC label if exist
4467 "ratio0="<<ratio0<< // shared ratio
4468 "ratio1="<<ratio1<< // shared ratio
4469 "p0.="<<&par0<< // track parameters
4471 "s0.="<<s1<< // full seed
4473 "n0="<<n0<< // number of clusters track 0
4474 "n1="<<n1<< // number of clusters track 1
4475 "nall0="<<nall0<< // number of clusters track 0
4476 "nall1="<<nall1<< // number of clusters track 1
4477 "n0F="<<n0F<< // number of findable
4478 "n1F="<<n1F<< // number of findable
4479 "shared="<<sumShared<< // number of shared clusters
4480 "firstS="<<firstShared<< // first and the last shared row
4481 "lastS="<<lastShared<<
4482 "firstRow="<<firstRow<< // first and the last row with cluster
4483 "lastRow="<<lastRow<< //
4487 // remove track with lower quality
4489 if (ratio0>AliTPCReconstructor::GetRecoParam()->GetCutSharedClusters(0) ||
4490 ratio1>AliTPCReconstructor::GetRecoParam()->GetCutSharedClusters(1)){
4494 delete array->RemoveAt(index1);
4499 // 4. Delete temporary array
4509 void AliTPCtrackerMI::FindCurling(const TObjArray * array, AliESDEvent */*esd*/, Int_t iter)
4512 // find Curling tracks
4513 // Use AliTPCReconstructor::StreamLevel()>1 if you want to tune parameters - cuts
4516 // Algorithm done in 2 phases - because of CPU consumption
4517 // it is n^2 algorithm - for lead-lead 20000x20000 combination are investigated
4518 // see detal in MC part what can be used to cut
4522 const Float_t kMaxC = 400; // maximal curvature to of the track
4523 const Float_t kMaxdTheta = 0.15; // maximal distance in theta
4524 const Float_t kMaxdPhi = 0.15; // maximal distance in phi
4525 const Float_t kPtRatio = 0.3; // ratio between pt
4526 const Float_t kMinDCAR = 2.; // distance to the primary vertex in r - see cpipe cut
4529 // Curling tracks cuts
4532 const Float_t kMaxDeltaRMax = 40; // distance in outer radius
4533 const Float_t kMaxDeltaRMin = 5.; // distance in lower radius - see cpipe cut
4534 const Float_t kMinAngle = 2.9; // angle between tracks
4535 const Float_t kMaxDist = 5; // biggest distance
4537 // The cuts can be tuned using the "MC information stored in Multi tree ==> see FindMultiMC
4540 TCut csign("csign","Tr0.fP[4]*Tr1.fP[4]<0"); //opposite sign
4541 TCut cmax("cmax","abs(Tr0.GetC())>1/400");
4542 TCut cda("cda","sqrt(dtheta^2+dfi^2)<0.15");
4543 TCut ccratio("ccratio","abs((Tr0.fP[4]+Tr1.fP[4])/(abs(Tr0.fP[4])+abs(Tr1.fP[4])))<0.3");
4544 TCut cpipe("cpipe", "min(abs(r0-rc0),abs(r1-rc1))>5");
4546 TCut cdrmax("cdrmax","abs(abs(rc0+r0)-abs(rc1+r1))<40")
4547 TCut cdrmin("cdrmin","abs(abs(rc0+r0)-abs(rc1+r1))<10")
4549 Multi->Draw("dfi","iter==0"+csign+cmax+cda+ccratio); ~94% of curling tracks fulfill
4550 Multi->Draw("min(abs(r0-rc0),abs(r1-rc1))","iter==0&&abs(lab1)==abs(lab0)"+csign+cmax+cda+ccratio+cpipe+cdrmin+cdrmax); //80%
4552 Curling2->Draw("dfi","iter==0&&abs(lab0)==abs(lab1)"+csign+cmax+cdtheta+cdfi+ccratio)
4558 Int_t nentries = array->GetEntriesFast();
4559 AliHelix *helixes = new AliHelix[nentries];
4560 for (Int_t i=0;i<nentries;i++){
4561 AliTPCseed* track = (AliTPCseed*)array->At(i);
4562 if (!track) continue;
4563 track->SetCircular(0);
4564 new (&helixes[i]) AliHelix(*track);
4570 Double_t phase[2][2],radius[2];
4575 for (Int_t i0=0;i0<nentries;i0++){
4576 AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
4577 if (!track0) continue;
4578 if (TMath::Abs(track0->GetC())<1/kMaxC) continue;
4579 Float_t xc0 = helixes[i0].GetHelix(6);
4580 Float_t yc0 = helixes[i0].GetHelix(7);
4581 Float_t r0 = helixes[i0].GetHelix(8);
4582 Float_t rc0 = TMath::Sqrt(xc0*xc0+yc0*yc0);
4583 Float_t fi0 = TMath::ATan2(yc0,xc0);
4585 for (Int_t i1=i0+1;i1<nentries;i1++){
4586 AliTPCseed * track1 = (AliTPCseed*)array->At(i1);
4587 if (!track1) continue;
4588 if (TMath::Abs(track1->GetC())<1/kMaxC) continue;
4589 Float_t xc1 = helixes[i1].GetHelix(6);
4590 Float_t yc1 = helixes[i1].GetHelix(7);
4591 Float_t r1 = helixes[i1].GetHelix(8);
4592 Float_t rc1 = TMath::Sqrt(xc1*xc1+yc1*yc1);
4593 Float_t fi1 = TMath::ATan2(yc1,xc1);
4595 Float_t dfi = fi0-fi1;
4598 if (dfi>1.5*TMath::Pi()) dfi-=TMath::Pi(); // take care about edge effect
4599 if (dfi<-1.5*TMath::Pi()) dfi+=TMath::Pi(); //
4600 Float_t dtheta = TMath::Abs(track0->GetTgl()-track1->GetTgl())<TMath::Abs(track0->GetTgl()+track1->GetTgl())? track0->GetTgl()-track1->GetTgl():track0->GetTgl()+track1->GetTgl();
4604 if (track0->GetBConstrain()&&track1->GetBConstrain()) continue; // not constrained
4605 if (track1->GetSigned1Pt()*track0->GetSigned1Pt()>0) continue; // not the same sign
4606 if ( TMath::Abs(track1->GetTgl()+track0->GetTgl())>kMaxdTheta) continue; //distance in the Theta
4607 if ( TMath::Abs(dfi)>kMaxdPhi) continue; //distance in phi
4608 if ( TMath::Sqrt(dfi*dfi+dtheta*dtheta)>kMaxdPhi) continue; //common angular offset
4610 Float_t pt0 = track0->GetSignedPt();
4611 Float_t pt1 = track1->GetSignedPt();
4612 if ((TMath::Abs(pt0+pt1)/(TMath::Abs(pt0)+TMath::Abs(pt1)))>kPtRatio) continue;
4613 if ((iter==1) && TMath::Abs(TMath::Abs(rc0+r0)-TMath::Abs(rc1+r1))>kMaxDeltaRMax) continue;
4614 if ((iter!=1) &&TMath::Abs(TMath::Abs(rc0-r0)-TMath::Abs(rc1-r1))>kMaxDeltaRMin) continue;
4615 if (TMath::Min(TMath::Abs(rc0-r0),TMath::Abs(rc1-r1))<kMinDCAR) continue;
4618 // Now find closest approach
4622 Int_t npoints = helixes[i0].GetRPHIintersections(helixes[i1], phase, radius,10);
4623 if (npoints==0) continue;
4624 helixes[i0].GetClosestPhases(helixes[i1], phase);
4628 Double_t hangles[3];
4629 helixes[i0].Evaluate(phase[0][0],xyz0);
4630 helixes[i1].Evaluate(phase[0][1],xyz1);
4632 helixes[i0].GetAngle(phase[0][0],helixes[i1],phase[0][1],hangles);
4633 Double_t deltah[2],deltabest;
4634 if (TMath::Abs(hangles[2])<kMinAngle) continue;
4638 helixes[i0].ParabolicDCA(helixes[i1],phase[0][0],phase[0][1],radius[0],deltah[0],2);
4640 helixes[i0].ParabolicDCA(helixes[i1],phase[1][0],phase[1][1],radius[1],deltah[1],2);
4641 if (deltah[1]<deltah[0]) ibest=1;
4643 deltabest = TMath::Sqrt(deltah[ibest]);
4644 helixes[i0].Evaluate(phase[ibest][0],xyz0);
4645 helixes[i1].Evaluate(phase[ibest][1],xyz1);
4646 helixes[i0].GetAngle(phase[ibest][0],helixes[i1],phase[ibest][1],hangles);
4647 Double_t radiusbest = TMath::Sqrt(radius[ibest]);
4649 if (deltabest>kMaxDist) continue;
4650 // if (mindcar+mindcaz<40 && (TMath::Abs(hangles[2])<kMinAngle ||deltabest>3)) continue;
4651 Bool_t sign =kFALSE;
4652 if (hangles[2]>kMinAngle) sign =kTRUE;
4655 // circular[i0] = kTRUE;
4656 // circular[i1] = kTRUE;
4657 if (track0->OneOverPt()<track1->OneOverPt()){
4658 track0->SetCircular(track0->GetCircular()+1);
4659 track1->SetCircular(track1->GetCircular()+2);
4662 track1->SetCircular(track1->GetCircular()+1);
4663 track0->SetCircular(track0->GetCircular()+2);
4666 if (AliTPCReconstructor::StreamLevel()>2){
4668 //debug stream to tune "fine" cuts
4669 Int_t lab0=track0->GetLabel();
4670 Int_t lab1=track1->GetLabel();
4671 TTreeSRedirector &cstream = *fDebugStreamer;
4672 cstream<<"Curling2"<<
4688 "npoints="<<npoints<<
4689 "hangles0="<<hangles[0]<<
4690 "hangles1="<<hangles[1]<<
4691 "hangles2="<<hangles[2]<<
4694 "radius="<<radiusbest<<
4695 "deltabest="<<deltabest<<
4696 "phase0="<<phase[ibest][0]<<
4697 "phase1="<<phase[ibest][1]<<
4705 if (AliTPCReconstructor::StreamLevel()>1) {
4706 AliInfo("Time for curling tracks removal");
4715 void AliTPCtrackerMI::FindKinks(TObjArray * array, AliESDEvent *esd)
4722 TObjArray *kinks= new TObjArray(10000);
4723 // TObjArray *v0s= new TObjArray(10000);
4724 Int_t nentries = array->GetEntriesFast();
4725 AliHelix *helixes = new AliHelix[nentries];
4726 Int_t *sign = new Int_t[nentries];
4727 Int_t *nclusters = new Int_t[nentries];
4728 Float_t *alpha = new Float_t[nentries];
4729 AliKink *kink = new AliKink();
4730 Int_t * usage = new Int_t[nentries];
4731 Float_t *zm = new Float_t[nentries];
4732 Float_t *z0 = new Float_t[nentries];
4733 Float_t *fim = new Float_t[nentries];
4734 Float_t *shared = new Float_t[nentries];
4735 Bool_t *circular = new Bool_t[nentries];
4736 Float_t *dca = new Float_t[nentries];
4737 //const AliESDVertex * primvertex = esd->GetVertex();
4739 // nentries = array->GetEntriesFast();
4744 for (Int_t i=0;i<nentries;i++){
4747 AliTPCseed* track = (AliTPCseed*)array->At(i);
4748 if (!track) continue;
4749 track->SetCircular(0);
4751 track->UpdatePoints();
4752 if (( track->GetPoints()[2]- track->GetPoints()[0])>5 && track->GetPoints()[3]>0.8){
4754 nclusters[i]=track->GetNumberOfClusters();
4755 alpha[i] = track->GetAlpha();
4756 new (&helixes[i]) AliHelix(*track);
4758 helixes[i].Evaluate(0,xyz);
4759 sign[i] = (track->GetC()>0) ? -1:1;
4762 if (track->GetProlongation(x,y,z)){
4764 fim[i] = alpha[i]+TMath::ATan2(y,x);
4767 zm[i] = track->GetZ();
4771 circular[i]= kFALSE;
4772 if (track->GetProlongation(0,y,z)) z0[i] = z;
4773 dca[i] = track->GetD(0,0);
4779 Int_t ncandidates =0;
4782 Double_t phase[2][2],radius[2];
4785 // Find circling track
4787 for (Int_t i0=0;i0<nentries;i0++){
4788 AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
4789 if (!track0) continue;
4790 if (track0->GetNumberOfClusters()<40) continue;
4791 if (TMath::Abs(1./track0->GetC())>200) continue;
4792 for (Int_t i1=i0+1;i1<nentries;i1++){
4793 AliTPCseed * track1 = (AliTPCseed*)array->At(i1);
4794 if (!track1) continue;
4795 if (track1->GetNumberOfClusters()<40) continue;
4796 if ( TMath::Abs(track1->GetTgl()+track0->GetTgl())>0.1) continue;
4797 if (track0->GetBConstrain()&&track1->GetBConstrain()) continue;
4798 if (TMath::Abs(1./track1->GetC())>200) continue;
4799 if (track1->GetSigned1Pt()*track0->GetSigned1Pt()>0) continue;
4800 if (track1->GetTgl()*track0->GetTgl()>0) continue;
4801 if (TMath::Max(TMath::Abs(1./track0->GetC()),TMath::Abs(1./track1->GetC()))>190) continue;
4802 if (track0->GetBConstrain()&&track1->OneOverPt()<track0->OneOverPt()) continue; //returning - lower momenta
4803 if (track1->GetBConstrain()&&track0->OneOverPt()<track1->OneOverPt()) continue; //returning - lower momenta
4805 Float_t mindcar = TMath::Min(TMath::Abs(dca[i0]),TMath::Abs(dca[i1]));
4806 if (mindcar<5) continue;
4807 Float_t mindcaz = TMath::Min(TMath::Abs(z0[i0]-GetZ()),TMath::Abs(z0[i1]-GetZ()));
4808 if (mindcaz<5) continue;
4809 if (mindcar+mindcaz<20) continue;
4812 Float_t xc0 = helixes[i0].GetHelix(6);
4813 Float_t yc0 = helixes[i0].GetHelix(7);
4814 Float_t r0 = helixes[i0].GetHelix(8);
4815 Float_t xc1 = helixes[i1].GetHelix(6);
4816 Float_t yc1 = helixes[i1].GetHelix(7);
4817 Float_t r1 = helixes[i1].GetHelix(8);
4819 Float_t rmean = (r0+r1)*0.5;
4820 Float_t delta =TMath::Sqrt((xc1-xc0)*(xc1-xc0)+(yc1-yc0)*(yc1-yc0));
4821 //if (delta>30) continue;
4822 if (delta>rmean*0.25) continue;
4823 if (TMath::Abs(r0-r1)/rmean>0.3) continue;
4825 Int_t npoints = helixes[i0].GetRPHIintersections(helixes[i1], phase, radius,10);
4826 if (npoints==0) continue;
4827 helixes[i0].GetClosestPhases(helixes[i1], phase);
4831 Double_t hangles[3];
4832 helixes[i0].Evaluate(phase[0][0],xyz0);
4833 helixes[i1].Evaluate(phase[0][1],xyz1);
4835 helixes[i0].GetAngle(phase[0][0],helixes[i1],phase[0][1],hangles);
4836 Double_t deltah[2],deltabest;
4837 if (hangles[2]<2.8) continue;
4840 helixes[i0].ParabolicDCA(helixes[i1],phase[0][0],phase[0][1],radius[0],deltah[0],2);
4842 helixes[i0].ParabolicDCA(helixes[i1],phase[1][0],phase[1][1],radius[1],deltah[1],2);
4843 if (deltah[1]<deltah[0]) ibest=1;
4845 deltabest = TMath::Sqrt(deltah[ibest]);
4846 helixes[i0].Evaluate(phase[ibest][0],xyz0);
4847 helixes[i1].Evaluate(phase[ibest][1],xyz1);
4848 helixes[i0].GetAngle(phase[ibest][0],helixes[i1],phase[ibest][1],hangles);
4849 Double_t radiusbest = TMath::Sqrt(radius[ibest]);
4851 if (deltabest>6) continue;
4852 if (mindcar+mindcaz<40 && (hangles[2]<3.12||deltabest>3)) continue;
4853 Bool_t lsign =kFALSE;
4854 if (hangles[2]>3.06) lsign =kTRUE;
4857 circular[i0] = kTRUE;
4858 circular[i1] = kTRUE;
4859 if (track0->OneOverPt()<track1->OneOverPt()){
4860 track0->SetCircular(track0->GetCircular()+1);
4861 track1->SetCircular(track1->GetCircular()+2);
4864 track1->SetCircular(track1->GetCircular()+1);
4865 track0->SetCircular(track0->GetCircular()+2);
4868 if (lsign&&AliTPCReconstructor::StreamLevel()>1){
4870 Int_t lab0=track0->GetLabel();
4871 Int_t lab1=track1->GetLabel();
4872 TTreeSRedirector &cstream = *fDebugStreamer;
4873 cstream<<"Curling"<<
4880 "mindcar="<<mindcar<<
4881 "mindcaz="<<mindcaz<<
4884 "npoints="<<npoints<<
4885 "hangles0="<<hangles[0]<<
4886 "hangles2="<<hangles[2]<<
4891 "radius="<<radiusbest<<
4892 "deltabest="<<deltabest<<
4893 "phase0="<<phase[ibest][0]<<
4894 "phase1="<<phase[ibest][1]<<
4904 for (Int_t i =0;i<nentries;i++){
4905 if (sign[i]==0) continue;
4906 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
4913 Double_t cradius0 = 40*40;
4914 Double_t cradius1 = 270*270;
4917 Double_t cdist3=0.55;
4918 for (Int_t j =i+1;j<nentries;j++){
4920 if (sign[j]*sign[i]<1) continue;
4921 if ( (nclusters[i]+nclusters[j])>200) continue;
4922 if ( (nclusters[i]+nclusters[j])<80) continue;
4923 if ( TMath::Abs(zm[i]-zm[j])>60.) continue;
4924 if ( TMath::Abs(fim[i]-fim[j])>0.6 && TMath::Abs(fim[i]-fim[j])<5.7 ) continue;
4925 //AliTPCseed * track1 = (AliTPCseed*)array->At(j); Double_t phase[2][2],radius[2];
4926 Int_t npoints = helixes[i].GetRPHIintersections(helixes[j], phase, radius,20);
4927 if (npoints<1) continue;
4930 if (radius[0]<cradius0||radius[0]>cradius1) continue;
4933 if ( (radius[0]<cradius0||radius[0]>cradius1) && (radius[1]<cradius0||radius[1]>cradius1) ) continue;
4936 Double_t delta1=10000,delta2=10000;
4937 // cuts on the intersection radius
4938 helixes[i].LinearDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta1);
4939 if (radius[0]<20&&delta1<1) continue; //intersection at vertex
4940 if (radius[0]<10&&delta1<3) continue; //intersection at vertex
4942 helixes[i].LinearDCA(helixes[j],phase[1][0],phase[1][1],radius[1],delta2);
4943 if (radius[1]<20&&delta2<1) continue; //intersection at vertex
4944 if (radius[1]<10&&delta2<3) continue; //intersection at vertex
4947 Double_t distance1 = TMath::Min(delta1,delta2);
4948 if (distance1>cdist1) continue; // cut on DCA linear approximation
4950 npoints = helixes[i].GetRPHIintersections(helixes[j], phase, radius,20);
4951 helixes[i].ParabolicDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta1);
4952 if (radius[0]<20&&delta1<1) continue; //intersection at vertex
4953 if (radius[0]<10&&delta1<3) continue; //intersection at vertex
4956 helixes[i].ParabolicDCA(helixes[j],phase[1][0],phase[1][1],radius[1],delta2);
4957 if (radius[1]<20&&delta2<1) continue; //intersection at vertex
4958 if (radius[1]<10&&delta2<3) continue; //intersection at vertex
4960 distance1 = TMath::Min(delta1,delta2);
4963 rkink = TMath::Sqrt(radius[0]);
4966 rkink = TMath::Sqrt(radius[1]);
4968 if (distance1>cdist2) continue;
4971 AliTPCseed * track1 = (AliTPCseed*)array->At(j);
4974 Int_t row0 = GetRowNumber(rkink);
4975 if (row0<10) continue;
4976 if (row0>150) continue;
4979 Float_t dens00=-1,dens01=-1;
4980 Float_t dens10=-1,dens11=-1;
4982 Int_t found,foundable,ishared;
4983 track0->GetClusterStatistic(0,row0-5, found, foundable,ishared,kFALSE);
4984 if (foundable>5) dens00 = Float_t(found)/Float_t(foundable);
4985 track0->GetClusterStatistic(row0+5,155, found, foundable,ishared,kFALSE);
4986 if (foundable>5) dens01 = Float_t(found)/Float_t(foundable);
4988 track1->GetClusterStatistic(0,row0-5, found, foundable,ishared,kFALSE);
4989 if (foundable>10) dens10 = Float_t(found)/Float_t(foundable);
4990 track1->GetClusterStatistic(row0+5,155, found, foundable,ishared,kFALSE);
4991 if (foundable>10) dens11 = Float_t(found)/Float_t(foundable);
4993 if (dens00<dens10 && dens01<dens11) continue;
4994 if (dens00>dens10 && dens01>dens11) continue;
4995 if (TMath::Max(dens00,dens10)<0.1) continue;
4996 if (TMath::Max(dens01,dens11)<0.3) continue;
4998 if (TMath::Min(dens00,dens10)>0.6) continue;
4999 if (TMath::Min(dens01,dens11)>0.6) continue;
5002 AliTPCseed * ktrack0, *ktrack1;
5011 if (TMath::Abs(ktrack0->GetC())>5) continue; // cut on the curvature for mother particle
5012 AliExternalTrackParam paramm(*ktrack0);
5013 AliExternalTrackParam paramd(*ktrack1);
5014 if (row0>60&&ktrack1->GetReference().GetX()>90.)new (¶md) AliExternalTrackParam(ktrack1->GetReference());
5017 kink->SetMother(paramm);
5018 kink->SetDaughter(paramd);
5021 Float_t x[3] = { kink->GetPosition()[0],kink->GetPosition()[1],kink->GetPosition()[2]};
5023 fkParam->Transform0to1(x,index);
5024 fkParam->Transform1to2(x,index);
5025 row0 = GetRowNumber(x[0]);
5027 if (kink->GetR()<100) continue;
5028 if (kink->GetR()>240) continue;
5029 if (kink->GetPosition()[2]/kink->GetR()>AliTPCReconstructor::GetCtgRange()) continue; //out of fiducial volume
5030 if (kink->GetDistance()>cdist3) continue;
5031 Float_t dird = kink->GetDaughterP()[0]*kink->GetPosition()[0]+kink->GetDaughterP()[1]*kink->GetPosition()[1]; // rough direction estimate
5032 if (dird<0) continue;
5034 Float_t dirm = kink->GetMotherP()[0]*kink->GetPosition()[0]+kink->GetMotherP()[1]*kink->GetPosition()[1]; // rough direction estimate
5035 if (dirm<0) continue;
5036 Float_t mpt = TMath::Sqrt(kink->GetMotherP()[0]*kink->GetMotherP()[0]+kink->GetMotherP()[1]*kink->GetMotherP()[1]);
5037 if (mpt<0.2) continue;
5040 //for high momenta momentum not defined well in first iteration
5041 Double_t qt = TMath::Sin(kink->GetAngle(2))*ktrack1->GetP();
5042 if (qt>0.35) continue;
5045 kink->SetLabel(CookLabel(ktrack0,0.4,0,row0),0);
5046 kink->SetLabel(CookLabel(ktrack1,0.4,row0,160),1);
5048 kink->SetTPCDensity(dens00,0,0);
5049 kink->SetTPCDensity(dens01,0,1);
5050 kink->SetTPCDensity(dens10,1,0);
5051 kink->SetTPCDensity(dens11,1,1);
5052 kink->SetIndex(i,0);
5053 kink->SetIndex(j,1);
5056 kink->SetTPCDensity(dens10,0,0);
5057 kink->SetTPCDensity(dens11,0,1);
5058 kink->SetTPCDensity(dens00,1,0);
5059 kink->SetTPCDensity(dens01,1,1);
5060 kink->SetIndex(j,0);
5061 kink->SetIndex(i,1);
5064 if (mpt<1||kink->GetAngle(2)>0.1){
5065 // angle and densities not defined yet
5066 if (kink->GetTPCDensityFactor()<0.8) continue;
5067 if ((2-kink->GetTPCDensityFactor())*kink->GetDistance() >0.25) continue;
5068 if (kink->GetAngle(2)*ktrack0->GetP()<0.003) continue; //too small angle
5069 if (kink->GetAngle(2)>0.2&&kink->GetTPCDensityFactor()<1.15) continue;
5070 if (kink->GetAngle(2)>0.2&&kink->GetTPCDensity(0,1)>0.05) continue;
5072 Float_t criticalangle = track0->GetSigmaSnp2()+track0->GetSigmaTgl2();
5073 criticalangle+= track1->GetSigmaSnp2()+track1->GetSigmaTgl2();
5074 criticalangle= 3*TMath::Sqrt(criticalangle);
5075 if (criticalangle>0.02) criticalangle=0.02;
5076 if (kink->GetAngle(2)<criticalangle) continue;
5079 Int_t drow = Int_t(2.+0.5/(0.05+kink->GetAngle(2))); // overlap region defined
5080 Float_t shapesum =0;
5082 for ( Int_t row = row0-drow; row<row0+drow;row++){
5083 if (row<0) continue;
5084 if (row>155) continue;
5085 if (ktrack0->GetClusterPointer(row)){
5086 AliTPCTrackerPoint *point =ktrack0->GetTrackPoint(row);
5087 shapesum+=point->GetSigmaY()+point->GetSigmaZ();
5090 if (ktrack1->GetClusterPointer(row)){
5091 AliTPCTrackerPoint *point =ktrack1->GetTrackPoint(row);
5092 shapesum+=point->GetSigmaY()+point->GetSigmaZ();
5097 kink->SetShapeFactor(-1.);
5100 kink->SetShapeFactor(shapesum/sum);
5102 // esd->AddKink(kink);
5104 // kink->SetMother(paramm);
5105 //kink->SetDaughter(paramd);
5107 Double_t chi2P2 = paramm.GetParameter()[2]-paramd.GetParameter()[2];
5109 chi2P2/=paramm.GetCovariance()[5]+paramd.GetCovariance()[5];
5110 Double_t chi2P3 = paramm.GetParameter()[3]-paramd.GetParameter()[3];
5112 chi2P3/=paramm.GetCovariance()[9]+paramd.GetCovariance()[9];
5114 if (AliTPCReconstructor::StreamLevel()>1) {
5115 (*fDebugStreamer)<<"kinkLpt"<<
5123 if ( chi2P2+chi2P3<AliTPCReconstructor::GetRecoParam()->GetKinkAngleCutChi2(0)){
5127 kinks->AddLast(kink);
5133 // sort the kinks according quality - and refit them towards vertex
5135 Int_t nkinks = kinks->GetEntriesFast();
5136 Float_t *quality = new Float_t[nkinks];
5137 Int_t *indexes = new Int_t[nkinks];
5138 AliTPCseed *mothers = new AliTPCseed[nkinks];
5139 AliTPCseed *daughters = new AliTPCseed[nkinks];
5142 for (Int_t i=0;i<nkinks;i++){
5144 AliKink *kinkl = (AliKink*)kinks->At(i);
5146 // refit kinks towards vertex
5148 Int_t index0 = kinkl->GetIndex(0);
5149 Int_t index1 = kinkl->GetIndex(1);
5150 AliTPCseed * ktrack0 = (AliTPCseed*)array->At(index0);
5151 AliTPCseed * ktrack1 = (AliTPCseed*)array->At(index1);
5153 Int_t sumn=ktrack0->GetNumberOfClusters()+ktrack1->GetNumberOfClusters();
5155 // Refit Kink under if too small angle
5157 if (kinkl->GetAngle(2)<0.05){
5158 kinkl->SetTPCRow0(GetRowNumber(kinkl->GetR()));
5159 Int_t row0 = kinkl->GetTPCRow0();
5160 Int_t drow = Int_t(2.+0.5/(0.05+kinkl->GetAngle(2)));
5163 Int_t last = row0-drow;
5164 if (last<40) last=40;
5165 if (last<ktrack0->GetFirstPoint()+25) last = ktrack0->GetFirstPoint()+25;
5166 AliTPCseed* seed0 = ReSeed(ktrack0,last,kFALSE);
5169 Int_t first = row0+drow;
5170 if (first>130) first=130;
5171 if (first>ktrack1->GetLastPoint()-25) first = TMath::Max(ktrack1->GetLastPoint()-25,30);
5172 AliTPCseed* seed1 = ReSeed(ktrack1,first,kTRUE);
5174 if (seed0 && seed1){
5175 kinkl->SetStatus(1,8);
5176 if (RefitKink(*seed0,*seed1,*kinkl)) kinkl->SetStatus(1,9);
5177 row0 = GetRowNumber(kinkl->GetR());
5178 sumn = seed0->GetNumberOfClusters()+seed1->GetNumberOfClusters();
5179 mothers[i] = *seed0;
5180 daughters[i] = *seed1;
5183 delete kinks->RemoveAt(i);
5184 if (seed0) delete seed0;
5185 if (seed1) delete seed1;
5188 if (kinkl->GetDistance()>0.5 || kinkl->GetR()<110 || kinkl->GetR()>240) {
5189 delete kinks->RemoveAt(i);
5190 if (seed0) delete seed0;
5191 if (seed1) delete seed1;
5199 if (kinkl) quality[i] = 160*((0.1+kinkl->GetDistance())*(2.-kinkl->GetTPCDensityFactor()))/(sumn+40.); //the longest -clossest will win
5201 TMath::Sort(nkinks,quality,indexes,kFALSE);
5203 //remove double find kinks
5205 for (Int_t ikink0=1;ikink0<nkinks;ikink0++){
5206 AliKink * kink0 = (AliKink*) kinks->At(indexes[ikink0]);
5207 if (!kink0) continue;
5209 for (Int_t ikink1=0;ikink1<ikink0;ikink1++){
5210 if (!kink0) continue;
5211 AliKink * kink1 = (AliKink*) kinks->At(indexes[ikink1]);
5212 if (!kink1) continue;
5213 // if not close kink continue
5214 if (TMath::Abs(kink1->GetPosition()[2]-kink0->GetPosition()[2])>10) continue;
5215 if (TMath::Abs(kink1->GetPosition()[1]-kink0->GetPosition()[1])>10) continue;
5216 if (TMath::Abs(kink1->GetPosition()[0]-kink0->GetPosition()[0])>10) continue;
5218 AliTPCseed &mother0 = mothers[indexes[ikink0]];
5219 AliTPCseed &daughter0 = daughters[indexes[ikink0]];
5220 AliTPCseed &mother1 = mothers[indexes[ikink1]];
5221 AliTPCseed &daughter1 = daughters[indexes[ikink1]];
5222 Int_t row0 = (kink0->GetTPCRow0()+kink1->GetTPCRow0())/2;
5231 for (Int_t i=0;i<row0;i++){
5232 if (mother0.GetClusterIndex(i)>0 && mother1.GetClusterIndex(i)>0){
5235 if (mother0.GetClusterIndex(i)==mother1.GetClusterIndex(i)){
5242 for (Int_t i=row0;i<158;i++){
5243 if (daughter0.GetClusterIndex(i)>0 && daughter0.GetClusterIndex(i)>0){
5246 if (mother0.GetClusterIndex(i)==mother1.GetClusterIndex(i)){
5252 Float_t ratio = Float_t(same+1)/Float_t(both+1);
5253 Float_t ratiom = Float_t(samem+1)/Float_t(bothm+1);
5254 Float_t ratiod = Float_t(samed+1)/Float_t(bothd+1);
5255 if (ratio>0.3 && ratiom>0.5 &&ratiod>0.5) {
5256 Int_t sum0 = mother0.GetNumberOfClusters()+daughter0.GetNumberOfClusters();
5257 Int_t sum1 = mother1.GetNumberOfClusters()+daughter1.GetNumberOfClusters();
5259 shared[kink0->GetIndex(0)]= kTRUE;
5260 shared[kink0->GetIndex(1)]= kTRUE;
5261 delete kinks->RemoveAt(indexes[ikink0]);
5264 shared[kink1->GetIndex(0)]= kTRUE;
5265 shared[kink1->GetIndex(1)]= kTRUE;
5266 delete kinks->RemoveAt(indexes[ikink1]);
5273 for (Int_t i=0;i<nkinks;i++){
5274 AliKink * kinkl = (AliKink*) kinks->At(indexes[i]);
5275 if (!kinkl) continue;
5276 kinkl->SetTPCRow0(GetRowNumber(kinkl->GetR()));
5277 Int_t index0 = kinkl->GetIndex(0);
5278 Int_t index1 = kinkl->GetIndex(1);
5279 if (circular[index0]||(circular[index1]&&kinkl->GetDistance()>0.2)) continue;
5280 kinkl->SetMultiple(usage[index0],0);
5281 kinkl->SetMultiple(usage[index1],1);
5282 if (kinkl->GetMultiple()[0]+kinkl->GetMultiple()[1]>2) continue;
5283 if (kinkl->GetMultiple()[0]+kinkl->GetMultiple()[1]>0 && quality[indexes[i]]>0.2) continue;
5284 if (kinkl->GetMultiple()[0]+kinkl->GetMultiple()[1]>0 && kinkl->GetDistance()>0.2) continue;
5285 if (circular[index0]||(circular[index1]&&kinkl->GetDistance()>0.1)) continue;
5287 AliTPCseed * ktrack0 = (AliTPCseed*)array->At(index0);
5288 AliTPCseed * ktrack1 = (AliTPCseed*)array->At(index1);
5289 if (!ktrack0 || !ktrack1) continue;
5290 Int_t index = esd->AddKink(kinkl);
5293 if ( ktrack0->GetKinkIndex(0)==0 && ktrack1->GetKinkIndex(0)==0) { //best kink
5294 if (mothers[indexes[i]].GetNumberOfClusters()>20 && daughters[indexes[i]].GetNumberOfClusters()>20 && (mothers[indexes[i]].GetNumberOfClusters()+daughters[indexes[i]].GetNumberOfClusters())>100){
5295 *ktrack0 = mothers[indexes[i]];
5296 *ktrack1 = daughters[indexes[i]];
5300 ktrack0->SetKinkIndex(usage[index0],-(index+1));
5301 ktrack1->SetKinkIndex(usage[index1], (index+1));
5306 // Remove tracks corresponding to shared kink's
5308 for (Int_t i=0;i<nentries;i++){
5309 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5310 if (!track0) continue;
5311 if (track0->GetKinkIndex(0)!=0) continue;
5312 if (shared[i]) delete array->RemoveAt(i);
5317 RemoveUsed2(array,0.5,0.4,30);
5319 for (Int_t i=0;i<nentries;i++){
5320 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5321 if (!track0) continue;
5322 track0->CookdEdx(0.02,0.6);
5326 for (Int_t i=0;i<nentries;i++){
5327 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5328 if (!track0) continue;
5329 if (track0->Pt()<1.4) continue;
5330 //remove double high momenta tracks - overlapped with kink candidates
5333 for (Int_t icl=track0->GetFirstPoint();icl<track0->GetLastPoint(); icl++){
5334 if (track0->GetClusterPointer(icl)!=0){
5336 if (track0->GetClusterPointer(icl)->IsUsed(10)) ishared++;
5339 if (Float_t(ishared+1)/Float_t(all+1)>0.5) {
5340 delete array->RemoveAt(i);
5344 if (track0->GetKinkIndex(0)!=0) continue;
5345 if (track0->GetNumberOfClusters()<80) continue;
5347 AliTPCseed *pmother = new AliTPCseed();
5348 AliTPCseed *pdaughter = new AliTPCseed();
5349 AliKink *pkink = new AliKink;
5351 AliTPCseed & mother = *pmother;
5352 AliTPCseed & daughter = *pdaughter;
5353 AliKink & kinkl = *pkink;
5354 if (CheckKinkPoint(track0,mother,daughter, kinkl)){
5355 if (mother.GetNumberOfClusters()<30||daughter.GetNumberOfClusters()<20) {
5359 continue; //too short tracks
5361 if (mother.Pt()<1.4) {
5367 Int_t row0= kinkl.GetTPCRow0();
5368 if (kinkl.GetDistance()>0.5 || kinkl.GetR()<110. || kinkl.GetR()>240.) {
5375 Int_t index = esd->AddKink(&kinkl);
5376 mother.SetKinkIndex(0,-(index+1));
5377 daughter.SetKinkIndex(0,index+1);
5378 if (mother.GetNumberOfClusters()>50) {
5379 delete array->RemoveAt(i);
5380 array->AddAt(new AliTPCseed(mother),i);
5383 array->AddLast(new AliTPCseed(mother));
5385 array->AddLast(new AliTPCseed(daughter));
5386 for (Int_t icl=0;icl<row0;icl++) {
5387 if (mother.GetClusterPointer(icl)) mother.GetClusterPointer(icl)->Use(20);
5390 for (Int_t icl=row0;icl<158;icl++) {
5391 if (daughter.GetClusterPointer(icl)) daughter.GetClusterPointer(icl)->Use(20);
5400 delete [] daughters;
5422 AliInfo(Form("Ncandidates=\t%d\t%d\t%d\t%d\n",esd->GetNumberOfKinks(),ncandidates,ntracks,nall));
5426 void AliTPCtrackerMI::FindV0s(const TObjArray * array, AliESDEvent *const esd)
5432 TObjArray *tpcv0s = new TObjArray(100000);
5433 Int_t nentries = array->GetEntriesFast();
5434 AliHelix *helixes = new AliHelix[nentries];
5435 Int_t *sign = new Int_t[nentries];
5436 Float_t *alpha = new Float_t[nentries];
5437 Float_t *z0 = new Float_t[nentries];
5438 Float_t *dca = new Float_t[nentries];
5439 Float_t *sdcar = new Float_t[nentries];
5440 Float_t *cdcar = new Float_t[nentries];
5441 Float_t *pulldcar = new Float_t[nentries];
5442 Float_t *pulldcaz = new Float_t[nentries];
5443 Float_t *pulldca = new Float_t[nentries];
5444 Bool_t *isPrim = new Bool_t[nentries];
5445 const AliESDVertex * primvertex = esd->GetVertex();
5446 Double_t zvertex = primvertex->GetZv();
5448 // nentries = array->GetEntriesFast();
5450 for (Int_t i=0;i<nentries;i++){
5453 AliTPCseed* track = (AliTPCseed*)array->At(i);
5454 if (!track) continue;
5455 track->GetV0Indexes()[0] = 0; //rest v0 indexes
5456 track->GetV0Indexes()[1] = 0; //rest v0 indexes
5457 track->GetV0Indexes()[2] = 0; //rest v0 indexes
5459 alpha[i] = track->GetAlpha();
5460 new (&helixes[i]) AliHelix(*track);
5462 helixes[i].Evaluate(0,xyz);
5463 sign[i] = (track->GetC()>0) ? -1:1;
5467 if (track->GetProlongation(0,y,z)) z0[i] = z;
5468 dca[i] = track->GetD(0,0);
5470 // dca error parrameterezation + pulls
5472 sdcar[i] = TMath::Sqrt(0.150*0.150+(100*track->GetC())*(100*track->GetC()));
5473 if (TMath::Abs(track->GetTgl())>1) sdcar[i]*=2.5;
5474 cdcar[i] = TMath::Exp((TMath::Abs(track->GetC())-0.0106)*525.3);
5475 pulldcar[i] = (dca[i]-cdcar[i])/sdcar[i];
5476 pulldcaz[i] = (z0[i]-zvertex)/sdcar[i];
5477 pulldca[i] = TMath::Sqrt(pulldcar[i]*pulldcar[i]+pulldcaz[i]*pulldcaz[i]);
5478 if (track->TPCrPID(1)+track->TPCrPID(2)+track->TPCrPID(3)>0.5) {
5479 if (pulldca[i]<3.) isPrim[i]=kTRUE; //pion, muon and Kaon 3 sigma cut
5481 if (track->TPCrPID(4)>0.5) {
5482 if (pulldca[i]<0.5) isPrim[i]=kTRUE; //proton 0.5 sigma cut
5484 if (track->TPCrPID(0)>0.4) {
5485 isPrim[i]=kFALSE; //electron no sigma cut
5492 Int_t ncandidates =0;
5495 Double_t phase[2][2],radius[2];
5501 Float_t fprimvertex[3]={GetX(),GetY(),GetZ()};
5503 Double_t cradius0 = 10*10;
5504 Double_t cradius1 = 200*200;
5507 Double_t cpointAngle = 0.95;
5509 Double_t delta[2]={10000,10000};
5510 for (Int_t i =0;i<nentries;i++){
5511 if (sign[i]==0) continue;
5512 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5513 if (!track0) continue;
5514 if (AliTPCReconstructor::StreamLevel()>1){
5515 TTreeSRedirector &cstream = *fDebugStreamer;
5520 "zvertex="<<zvertex<<
5521 "sdcar0="<<sdcar[i]<<
5522 "cdcar0="<<cdcar[i]<<
5523 "pulldcar0="<<pulldcar[i]<<
5524 "pulldcaz0="<<pulldcaz[i]<<
5525 "pulldca0="<<pulldca[i]<<
5526 "isPrim="<<isPrim[i]<<
5530 if (track0->GetSigned1Pt()<0) continue;
5531 if (track0->GetKinkIndex(0)>0||isPrim[i]) continue; //daughter kink
5533 if (TMath::Abs(helixes[i].GetHelix(4))<0.000000001) continue;
5538 for (Int_t j =0;j<nentries;j++){
5539 AliTPCseed * track1 = (AliTPCseed*)array->At(j);
5540 if (!track1) continue;
5541 if (track1->GetKinkIndex(0)>0 || isPrim[j]) continue; //daughter kink
5542 if (sign[j]*sign[i]>0) continue;
5543 if (TMath::Abs(helixes[j].GetHelix(4))<0.000001) continue;
5544 if (track0->GetCircular()+track1->GetCircular()>1) continue; //circling -returning track
5547 // DCA to prim vertex cut
5553 Int_t npoints = helixes[i].GetRPHIintersections(helixes[j], phase, radius,cdist2);
5554 if (npoints<1) continue;
5558 if (radius[0]<cradius0||radius[0]>cradius1) continue;
5559 helixes[i].LinearDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta[0]);
5560 if (delta[0]>cdist1) continue;
5563 if (TMath::Max(radius[0],radius[1])<cradius0|| TMath::Min(radius[0],radius[1])>cradius1) continue;
5564 helixes[i].LinearDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta[0]);
5565 helixes[i].LinearDCA(helixes[j],phase[1][0],phase[1][1],radius[1],delta[1]);
5566 if (delta[1]<delta[0]) iclosest=1;
5567 if (delta[iclosest]>cdist1) continue;
5569 helixes[i].ParabolicDCA(helixes[j],phase[iclosest][0],phase[iclosest][1],radius[iclosest],delta[iclosest]);
5570 if (radius[iclosest]<cradius0 || radius[iclosest]>cradius1 || delta[iclosest]>cdist1) continue;
5572 Double_t pointAngle = helixes[i].GetPointAngle(helixes[j],phase[iclosest],fprimvertex);
5573 if (pointAngle<cpointAngle) continue;
5575 Bool_t isGamma = kFALSE;
5576 vertex.SetParamP(*track0); //track0 - plus
5577 vertex.SetParamN(*track1); //track1 - minus
5578 vertex.Update(fprimvertex);
5579 if (track0->TPCrPID(0)>0.3&&track1->TPCrPID(0)>0.3&&vertex.GetAnglep()[2]<0.15) isGamma=kTRUE; // gamma conversion candidate
5580 Double_t pointAngle2 = vertex.GetV0CosineOfPointingAngle();
5582 if (vertex.GetV0CosineOfPointingAngle()<cpointAngle && (!isGamma)) continue;// point angle cut
5583 //Bo: if (vertex.GetDist2()>2&&(!isGamma)) continue; // point angle cut
5584 if (vertex.GetDcaV0Daughters()>2&&(!isGamma)) continue;//Bo: // point angle cut
5585 Float_t sigmae = 0.15*0.15;
5586 if (vertex.GetRr()<80)
5587 sigmae += (sdcar[i]*sdcar[i]+sdcar[j]*sdcar[j])*(1.-vertex.GetRr()/80.)*(1.-vertex.GetRr()/80.);
5588 sigmae+= TMath::Sqrt(sigmae);
5589 //Bo: if (vertex.GetDist2()/sigmae>3.&&(!isGamma)) continue;
5590 if (vertex.GetDcaV0Daughters()/sigmae>3.&&(!isGamma)) continue;
5591 Float_t densb0=0,densb1=0,densa0=0,densa1=0;
5592 Int_t row0 = GetRowNumber(vertex.GetRr());
5594 //Bo: if (vertex.GetDist2()>0.2) continue;
5595 if (vertex.GetDcaV0Daughters()>0.2) continue;
5596 densb0 = track0->Density2(0,row0-5);
5597 densb1 = track1->Density2(0,row0-5);
5598 if (densb0>0.3|| densb1>0.3) continue; //clusters before vertex
5599 densa0 = track0->Density2(row0+5,row0+40);
5600 densa1 = track1->Density2(row0+5,row0+40);
5601 if ((densa0<0.4|| densa1<0.4)&&(!isGamma)) continue; //missing clusters after vertex
5604 densa0 = track0->Density2(0,40); //cluster density
5605 densa1 = track1->Density2(0,40); //cluster density
5606 if ((vertex.GetRr()<80&&densa0+densa1<1.)&&(!isGamma)) continue;
5608 //Bo: vertex.SetLab(0,track0->GetLabel());
5609 //Bo: vertex.SetLab(1,track1->GetLabel());
5610 vertex.SetChi2After((densa0+densa1)*0.5);
5611 vertex.SetChi2Before((densb0+densb1)*0.5);
5612 vertex.SetIndex(0,i);
5613 vertex.SetIndex(1,j);
5614 //Bo: vertex.SetStatus(1); // TPC v0 candidate
5615 vertex.SetOnFlyStatus(2);//Bo: // TPC v0 candidate
5616 //Bo: vertex.SetRp(track0->TPCrPIDs());
5617 //Bo: vertex.SetRm(track1->TPCrPIDs());
5618 tpcv0s->AddLast(new AliESDv0(vertex));
5621 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
5622 Double_t radiusm= (delta[0]<delta[1])? TMath::Sqrt(radius[0]):TMath::Sqrt(radius[1]);
5623 Double_t deltam= (delta[0]<delta[1])? TMath::Sqrt(delta[0]):TMath::Sqrt(delta[1]);
5624 if (AliTPCReconstructor::StreamLevel()>1) {
5625 Int_t lab0=track0->GetLabel();
5626 Int_t lab1=track1->GetLabel();
5627 Char_t c0=track0->GetCircular();
5628 Char_t c1=track1->GetCircular();
5629 TTreeSRedirector &cstream = *fDebugStreamer;
5632 "vertex.="<<&vertex<<
5635 "Helix0.="<<&helixes[i]<<
5638 "Helix1.="<<&helixes[j]<<
5639 "pointAngle="<<pointAngle<<
5640 "pointAngle2="<<pointAngle2<<
5645 "zvertex="<<zvertex<<
5648 "npoints="<<npoints<<
5649 "radius0="<<radius[0]<<
5650 "delta0="<<delta[0]<<
5651 "radius1="<<radius[1]<<
5652 "delta1="<<delta[1]<<
5653 "radiusm="<<radiusm<<
5655 "sdcar0="<<sdcar[i]<<
5656 "sdcar1="<<sdcar[j]<<
5657 "cdcar0="<<cdcar[i]<<
5658 "cdcar1="<<cdcar[j]<<
5659 "pulldcar0="<<pulldcar[i]<<
5660 "pulldcar1="<<pulldcar[j]<<
5661 "pulldcaz0="<<pulldcaz[i]<<
5662 "pulldcaz1="<<pulldcaz[j]<<
5663 "pulldca0="<<pulldca[i]<<
5664 "pulldca1="<<pulldca[j]<<
5674 Float_t *quality = new Float_t[ncandidates];
5675 Int_t *indexes = new Int_t[ncandidates];
5677 for (Int_t i=0;i<ncandidates;i++){
5679 AliESDv0 *v0 = (AliESDv0*)tpcv0s->At(i);
5680 quality[i] = 1./(1.00001-v0->GetV0CosineOfPointingAngle()); //base point angle
5681 // quality[i] /= (0.5+v0->GetDist2());
5682 // quality[i] *= v0->GetChi2After(); //density factor
5684 Int_t index0 = v0->GetIndex(0);
5685 Int_t index1 = v0->GetIndex(1);
5686 //Bo: Double_t minpulldca = TMath::Min(2.+pulldca[v0->GetIndex(0)],(2.+pulldca[v0->GetIndex(1)]) ); //pull
5687 Double_t minpulldca = TMath::Min(2.+pulldca[index0],(2.+pulldca[index1]) );//Bo:
5691 AliTPCseed * track0 = (AliTPCseed*)array->At(index0);
5692 AliTPCseed * track1 = (AliTPCseed*)array->At(index1);
5693 if (track0->TPCrPID(0)>0.3&&track1->TPCrPID(0)>0.3&&v0->GetAnglep()[2]<0.15) quality[i]+=1000000; // gamma conversion candidate
5694 if (track0->TPCrPID(4)>0.9||(track1->TPCrPID(4)>0.9&&minpulldca>4)) quality[i]*=10; // lambda candidate candidate
5697 TMath::Sort(ncandidates,quality,indexes,kTRUE);
5700 for (Int_t i=0;i<ncandidates;i++){
5701 AliESDv0 * v0 = (AliESDv0*)tpcv0s->At(indexes[i]);
5703 Int_t index0 = v0->GetIndex(0);
5704 Int_t index1 = v0->GetIndex(1);
5705 AliTPCseed * track0 = (AliTPCseed*)array->At(index0);
5706 AliTPCseed * track1 = (AliTPCseed*)array->At(index1);
5707 if (!track0||!track1) {
5711 Bool_t accept =kTRUE; //default accept
5712 Int_t *v0indexes0 = track0->GetV0Indexes();
5713 Int_t *v0indexes1 = track1->GetV0Indexes();
5715 Int_t order0 = (v0indexes0[0]!=0) ? 1:0;
5716 Int_t order1 = (v0indexes1[0]!=0) ? 1:0;
5717 if (v0indexes0[1]!=0) order0 =2;
5718 if (v0indexes1[1]!=0) order1 =2;
5720 if (v0indexes0[2]!=0) {order0=3; accept=kFALSE;}
5721 if (v0indexes0[2]!=0) {order1=3; accept=kFALSE;}
5723 AliESDv0 * v02 = v0;
5725 //Bo: v0->SetOrder(0,order0);
5726 //Bo: v0->SetOrder(1,order1);
5727 //Bo: v0->SetOrder(1,order0+order1);
5728 v0->SetOnFlyStatus(kTRUE);
5729 Int_t index = esd->AddV0(v0);
5730 v02 = esd->GetV0(index);
5731 v0indexes0[order0]=index;
5732 v0indexes1[order1]=index;
5736 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
5737 if (AliTPCReconstructor::StreamLevel()>1) {
5738 Int_t lab0=track0->GetLabel();
5739 Int_t lab1=track1->GetLabel();
5740 TTreeSRedirector &cstream = *fDebugStreamer;
5749 "dca0="<<dca[index0]<<
5750 "dca1="<<dca[index1]<<
5754 "quality="<<quality[i]<<
5755 "pulldca0="<<pulldca[index0]<<
5756 "pulldca1="<<pulldca[index1]<<
5780 printf("TPC V0 finder : naccepted\t%d\tncandidates\t%d\tntracks\t%d\tnall\t%d\n",naccepted,ncandidates,ntracks,nall);
5784 Int_t AliTPCtrackerMI::RefitKink(AliTPCseed &mother, AliTPCseed &daughter, const AliESDkink &knk)
5787 // refit kink towards to the vertex
5790 AliKink &kink=(AliKink &)knk;
5792 Int_t row0 = GetRowNumber(kink.GetR());
5793 FollowProlongation(mother,0);
5794 mother.Reset(kFALSE);
5796 FollowProlongation(daughter,row0);
5797 daughter.Reset(kFALSE);
5798 FollowBackProlongation(daughter,158);
5799 daughter.Reset(kFALSE);
5800 Int_t first = TMath::Max(row0-20,30);
5801 Int_t last = TMath::Min(row0+20,140);
5803 const Int_t kNdiv =5;
5804 AliTPCseed param0[kNdiv]; // parameters along the track
5805 AliTPCseed param1[kNdiv]; // parameters along the track
5806 AliKink kinks[kNdiv]; // corresponding kink parameters
5809 for (Int_t irow=0; irow<kNdiv;irow++){
5810 rows[irow] = first +((last-first)*irow)/(kNdiv-1);
5812 // store parameters along the track
5814 for (Int_t irow=0;irow<kNdiv;irow++){
5815 FollowBackProlongation(mother, rows[irow]);
5816 FollowProlongation(daughter,rows[kNdiv-1-irow]);
5817 param0[irow] = mother;
5818 param1[kNdiv-1-irow] = daughter;
5822 for (Int_t irow=0; irow<kNdiv-1;irow++){
5823 if (param0[irow].GetNumberOfClusters()<kNdiv||param1[irow].GetNumberOfClusters()<kNdiv) continue;
5824 kinks[irow].SetMother(param0[irow]);
5825 kinks[irow].SetDaughter(param1[irow]);
5826 kinks[irow].Update();
5829 // choose kink with best "quality"
5831 Double_t mindist = 10000;
5832 for (Int_t irow=0;irow<kNdiv;irow++){
5833 if (param0[irow].GetNumberOfClusters()<20||param1[irow].GetNumberOfClusters()<20) continue;
5834 if (TMath::Abs(kinks[irow].GetR())>240.) continue;
5835 if (TMath::Abs(kinks[irow].GetR())<100.) continue;
5837 Float_t normdist = TMath::Abs(param0[irow].GetX()-kinks[irow].GetR())*(0.1+kink.GetDistance());
5838 normdist/= (param0[irow].GetNumberOfClusters()+param1[irow].GetNumberOfClusters()+40.);
5839 if (normdist < mindist){
5845 if (index==-1) return 0;
5848 param0[index].Reset(kTRUE);
5849 FollowProlongation(param0[index],0);
5851 mother = param0[index];
5852 daughter = param1[index]; // daughter in vertex
5854 kink.SetMother(mother);
5855 kink.SetDaughter(daughter);
5857 kink.SetTPCRow0(GetRowNumber(kink.GetR()));
5858 kink.SetTPCncls(param0[index].GetNumberOfClusters(),0);
5859 kink.SetTPCncls(param1[index].GetNumberOfClusters(),1);
5860 kink.SetLabel(CookLabel(&mother,0.4, 0,kink.GetTPCRow0()),0);
5861 kink.SetLabel(CookLabel(&daughter,0.4, kink.GetTPCRow0(),160),1);
5862 mother.SetLabel(kink.GetLabel(0));
5863 daughter.SetLabel(kink.GetLabel(1));
5869 void AliTPCtrackerMI::UpdateKinkQualityM(AliTPCseed * seed){
5871 // update Kink quality information for mother after back propagation
5873 if (seed->GetKinkIndex(0)>=0) return;
5874 for (Int_t ikink=0;ikink<3;ikink++){
5875 Int_t index = seed->GetKinkIndex(ikink);
5876 if (index>=0) break;
5877 index = TMath::Abs(index)-1;
5878 AliESDkink * kink = fEvent->GetKink(index);
5879 kink->SetTPCDensity(-1,0,0);
5880 kink->SetTPCDensity(1,0,1);
5882 Int_t row0 = kink->GetTPCRow0() - 2 - Int_t( 0.5/ (0.05+kink->GetAngle(2)));
5883 if (row0<15) row0=15;
5885 Int_t row1 = kink->GetTPCRow0() + 2 + Int_t( 0.5/ (0.05+kink->GetAngle(2)));
5886 if (row1>145) row1=145;
5888 Int_t found,foundable,shared;
5889 seed->GetClusterStatistic(0,row0, found, foundable,shared,kFALSE);
5890 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),0,0);
5891 seed->GetClusterStatistic(row1,155, found, foundable,shared,kFALSE);
5892 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),0,1);
5897 void AliTPCtrackerMI::UpdateKinkQualityD(AliTPCseed * seed){
5899 // update Kink quality information for daughter after refit
5901 if (seed->GetKinkIndex(0)<=0) return;
5902 for (Int_t ikink=0;ikink<3;ikink++){
5903 Int_t index = seed->GetKinkIndex(ikink);
5904 if (index<=0) break;
5905 index = TMath::Abs(index)-1;
5906 AliESDkink * kink = fEvent->GetKink(index);
5907 kink->SetTPCDensity(-1,1,0);
5908 kink->SetTPCDensity(-1,1,1);
5910 Int_t row0 = kink->GetTPCRow0() -2 - Int_t( 0.5/ (0.05+kink->GetAngle(2)));
5911 if (row0<15) row0=15;
5913 Int_t row1 = kink->GetTPCRow0() +2 + Int_t( 0.5/ (0.05+kink->GetAngle(2)));
5914 if (row1>145) row1=145;
5916 Int_t found,foundable,shared;
5917 seed->GetClusterStatistic(0,row0, found, foundable,shared,kFALSE);
5918 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),1,0);
5919 seed->GetClusterStatistic(row1,155, found, foundable,shared,kFALSE);
5920 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),1,1);
5926 Int_t AliTPCtrackerMI::CheckKinkPoint(AliTPCseed*seed,AliTPCseed &mother, AliTPCseed &daughter, const AliESDkink &knk)
5929 // check kink point for given track
5930 // if return value=0 kink point not found
5931 // otherwise seed0 correspond to mother particle
5932 // seed1 correspond to daughter particle
5933 // kink parameter of kink point
5934 AliKink &kink=(AliKink &)knk;
5936 Int_t middlerow = (seed->GetFirstPoint()+seed->GetLastPoint())/2;
5937 Int_t first = seed->GetFirstPoint();
5938 Int_t last = seed->GetLastPoint();
5939 if (last-first<20) return 0; // shortest length - 2*30 = 60 pad-rows
5942 AliTPCseed *seed1 = ReSeed(seed,middlerow+20, kTRUE); //middle of chamber
5943 if (!seed1) return 0;
5944 FollowProlongation(*seed1,seed->GetLastPoint()-20);
5945 seed1->Reset(kTRUE);
5946 FollowProlongation(*seed1,158);
5947 seed1->Reset(kTRUE);
5948 last = seed1->GetLastPoint();
5950 AliTPCseed *seed0 = new AliTPCseed(*seed);
5951 seed0->Reset(kFALSE);
5954 AliTPCseed param0[20]; // parameters along the track
5955 AliTPCseed param1[20]; // parameters along the track
5956 AliKink kinks[20]; // corresponding kink parameters
5958 for (Int_t irow=0; irow<20;irow++){
5959 rows[irow] = first +((last-first)*irow)/19;
5961 // store parameters along the track
5963 for (Int_t irow=0;irow<20;irow++){
5964 FollowBackProlongation(*seed0, rows[irow]);
5965 FollowProlongation(*seed1,rows[19-irow]);
5966 param0[irow] = *seed0;
5967 param1[19-irow] = *seed1;
5971 for (Int_t irow=0; irow<19;irow++){
5972 kinks[irow].SetMother(param0[irow]);
5973 kinks[irow].SetDaughter(param1[irow]);
5974 kinks[irow].Update();
5977 // choose kink with biggest change of angle
5979 Double_t maxchange= 0;
5980 for (Int_t irow=1;irow<19;irow++){
5981 if (TMath::Abs(kinks[irow].GetR())>240.) continue;
5982 if (TMath::Abs(kinks[irow].GetR())<110.) continue;
5983 Float_t quality = TMath::Abs(kinks[irow].GetAngle(2))/(3.+TMath::Abs(kinks[irow].GetR()-param0[irow].GetX()));
5984 if ( quality > maxchange){
5985 maxchange = quality;
5992 if (index<0) return 0;
5994 Int_t row0 = GetRowNumber(kinks[index].GetR()); //row 0 estimate
5995 seed0 = new AliTPCseed(param0[index]);
5996 seed1 = new AliTPCseed(param1[index]);
5997 seed0->Reset(kFALSE);
5998 seed1->Reset(kFALSE);
5999 seed0->ResetCovariance(10.);
6000 seed1->ResetCovariance(10.);
6001 FollowProlongation(*seed0,0);
6002 FollowBackProlongation(*seed1,158);
6003 mother = *seed0; // backup mother at position 0
6004 seed0->Reset(kFALSE);
6005 seed1->Reset(kFALSE);
6006 seed0->ResetCovariance(10.);
6007 seed1->ResetCovariance(10.);
6009 first = TMath::Max(row0-20,0);
6010 last = TMath::Min(row0+20,158);
6012 for (Int_t irow=0; irow<20;irow++){
6013 rows[irow] = first +((last-first)*irow)/19;
6015 // store parameters along the track
6017 for (Int_t irow=0;irow<20;irow++){
6018 FollowBackProlongation(*seed0, rows[irow]);
6019 FollowProlongation(*seed1,rows[19-irow]);
6020 param0[irow] = *seed0;
6021 param1[19-irow] = *seed1;
6025 for (Int_t irow=0; irow<19;irow++){
6026 kinks[irow].SetMother(param0[irow]);
6027 kinks[irow].SetDaughter(param1[irow]);
6028 // param0[irow].Dump();
6029 //param1[irow].Dump();
6030 kinks[irow].Update();
6033 // choose kink with biggest change of angle
6036 for (Int_t irow=0;irow<20;irow++){
6037 if (TMath::Abs(kinks[irow].GetR())>250.) continue;
6038 if (TMath::Abs(kinks[irow].GetR())<90.) continue;
6039 Float_t quality = TMath::Abs(kinks[irow].GetAngle(2))/(3.+TMath::Abs(kinks[irow].GetR()-param0[irow].GetX()));
6040 if ( quality > maxchange){
6041 maxchange = quality;
6048 if (index==-1 || param0[index].GetNumberOfClusters()+param1[index].GetNumberOfClusters()<100){
6054 // Float_t anglesigma = TMath::Sqrt(param0[index].fC22+param0[index].fC33+param1[index].fC22+param1[index].fC33);
6056 kink.SetMother(param0[index]);
6057 kink.SetDaughter(param1[index]);
6060 Double_t chi2P2 = param0[index].GetParameter()[2]-param1[index].GetParameter()[2];
6062 chi2P2/=param0[index].GetCovariance()[5]+param1[index].GetCovariance()[5];
6063 Double_t chi2P3 = param0[index].GetParameter()[3]-param1[index].GetParameter()[3];
6065 chi2P3/=param0[index].GetCovariance()[9]+param1[index].GetCovariance()[9];
6067 if (AliTPCReconstructor::StreamLevel()>1) {
6068 (*fDebugStreamer)<<"kinkHpt"<<
6071 "p0.="<<¶m0[index]<<
6072 "p1.="<<¶m1[index]<<
6076 if ( chi2P2+chi2P3<AliTPCReconstructor::GetRecoParam()->GetKinkAngleCutChi2(0)){
6083 row0 = GetRowNumber(kink.GetR());
6084 kink.SetTPCRow0(row0);
6085 kink.SetLabel(CookLabel(seed0,0.5,0,row0),0);
6086 kink.SetLabel(CookLabel(seed1,0.5,row0,158),1);
6087 kink.SetIndex(-10,0);
6088 kink.SetIndex(int(param0[index].GetNumberOfClusters()+param1[index].GetNumberOfClusters()),1);
6089 kink.SetTPCncls(param0[index].GetNumberOfClusters(),0);
6090 kink.SetTPCncls(param1[index].GetNumberOfClusters(),1);
6093 // new (&mother) AliTPCseed(param0[index]);
6094 daughter = param1[index];
6095 daughter.SetLabel(kink.GetLabel(1));
6096 param0[index].Reset(kTRUE);
6097 FollowProlongation(param0[index],0);
6098 mother = param0[index];
6099 mother.SetLabel(kink.GetLabel(0));
6100 if ( chi2P2+chi2P3<AliTPCReconstructor::GetRecoParam()->GetKinkAngleCutChi2(1)){
6112 AliTPCseed* AliTPCtrackerMI::ReSeed(AliTPCseed *t)
6115 // reseed - refit - track
6118 // Int_t last = fSectors->GetNRows()-1;
6120 if (fSectors == fOuterSec){
6121 first = TMath::Max(first, t->GetFirstPoint()-fInnerSec->GetNRows());
6125 first = t->GetFirstPoint();
6127 AliTPCseed * seed = MakeSeed(t,0.1,0.5,0.9);
6128 FollowBackProlongation(*t,fSectors->GetNRows()-1);
6130 FollowProlongation(*t,first);
6140 //_____________________________________________________________________________
6141 Int_t AliTPCtrackerMI::ReadSeeds(const TFile *inp) {
6142 //-----------------------------------------------------------------
6143 // This function reades track seeds.
6144 //-----------------------------------------------------------------
6145 TDirectory *savedir=gDirectory;
6147 TFile *in=(TFile*)inp;
6148 if (!in->IsOpen()) {
6149 cerr<<"AliTPCtrackerMI::ReadSeeds(): input file is not open !\n";
6154 TTree *seedTree=(TTree*)in->Get("Seeds");
6156 cerr<<"AliTPCtrackerMI::ReadSeeds(): ";
6157 cerr<<"can't get a tree with track seeds !\n";
6160 AliTPCtrack *seed=new AliTPCtrack;
6161 seedTree->SetBranchAddress("tracks",&seed);
6163 if (fSeeds==0) fSeeds=new TObjArray(15000);
6165 Int_t n=(Int_t)seedTree->GetEntries();
6166 for (Int_t i=0; i<n; i++) {
6167 seedTree->GetEvent(i);
6168 fSeeds->AddLast(new AliTPCseed(*seed/*,seed->GetAlpha()*/));
6177 Int_t AliTPCtrackerMI::Clusters2Tracks (AliESDEvent *const esd)
6180 if (fSeeds) DeleteSeeds();
6183 if (!fSeeds) return 1;
6185 if (AliTPCReconstructor::StreamLevel()>3) DumpClusters(0,fSeeds);
6191 //_____________________________________________________________________________
6192 Int_t AliTPCtrackerMI::Clusters2Tracks() {
6193 //-----------------------------------------------------------------
6194 // This is a track finder.
6195 //-----------------------------------------------------------------
6196 TDirectory *savedir=gDirectory;
6200 fSeeds = Tracking();
6203 Info("Clusters2Tracks","Time for tracking: \t");timer.Print();timer.Start();
6205 //activate again some tracks
6206 for (Int_t i=0; i<fSeeds->GetEntriesFast(); i++) {
6207 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
6209 Int_t nc=t.GetNumberOfClusters();
6211 delete fSeeds->RemoveAt(i);
6215 if (pt->GetRemoval()==10) {
6216 if (pt->GetDensityFirst(20)>0.8 || pt->GetDensityFirst(30)>0.8 || pt->GetDensityFirst(40)>0.7)
6217 pt->Desactivate(10); // make track again active
6219 pt->Desactivate(20);
6220 delete fSeeds->RemoveAt(i);
6225 RemoveUsed2(fSeeds,0.85,0.85,0);
6226 if (AliTPCReconstructor::GetRecoParam()->GetDoKinks()) FindKinks(fSeeds,fEvent);
6227 //FindCurling(fSeeds, fEvent,0);
6228 if (AliTPCReconstructor::StreamLevel()>5) FindMultiMC(fSeeds, fEvent,-1); // find multi found tracks
6229 RemoveUsed2(fSeeds,0.5,0.4,20);
6230 FindSplitted(fSeeds, fEvent,0); // find multi found tracks
6231 if (AliTPCReconstructor::StreamLevel()>5) FindMultiMC(fSeeds, fEvent,0); // find multi found tracks
6234 // // refit short tracks
6236 Int_t nseed=fSeeds->GetEntriesFast();
6239 for (Int_t i=0; i<nseed; i++) {
6240 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
6242 Int_t nc=t.GetNumberOfClusters();
6244 delete fSeeds->RemoveAt(i);
6247 CookLabel(pt,0.1); //For comparison only
6248 //if ((pt->IsActive() || (pt->fRemoval==10) )&& nc>50 &&pt->GetNumberOfClusters()>0.4*pt->fNFoundable){
6249 if ((pt->IsActive() || (pt->GetRemoval()==10) )){
6251 if (fDebug>0) cerr<<found<<'\r';
6255 delete fSeeds->RemoveAt(i);
6259 //RemoveOverlap(fSeeds,0.99,7,kTRUE);
6261 //RemoveUsed(fSeeds,0.9,0.9,6);
6263 nseed=fSeeds->GetEntriesFast();
6265 for (Int_t i=0; i<nseed; i++) {
6266 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
6268 Int_t nc=t.GetNumberOfClusters();
6270 delete fSeeds->RemoveAt(i);
6274 t.CookdEdx(0.02,0.6);
6275 // CheckKinkPoint(&t,0.05);
6276 //if ((pt->IsActive() || (pt->fRemoval==10) )&& nc>50 &&pt->GetNumberOfClusters()>0.4*pt->fNFoundable){
6277 if ((pt->IsActive() || (pt->GetRemoval()==10) )){
6285 delete fSeeds->RemoveAt(i);
6286 //AliTPCseed * seed1 = ReSeed(pt,0.05,0.5,1);
6288 // FollowProlongation(*seed1,0);
6289 // Int_t n = seed1->GetNumberOfClusters();
6290 // printf("fP4\t%f\t%f\n",seed1->GetC(),pt->GetC());
6291 // printf("fN\t%d\t%d\n", seed1->GetNumberOfClusters(),pt->GetNumberOfClusters());
6294 //AliTPCseed * seed2 = ReSeed(pt,0.95,0.5,0.05);
6298 SortTracks(fSeeds, 1);
6302 PrepareForBackProlongation(fSeeds,5.);
6303 PropagateBack(fSeeds);
6304 printf("Time for back propagation: \t");timer.Print();timer.Start();
6308 PrepareForProlongation(fSeeds,5.);
6309 PropagateForward2(fSeeds);
6311 printf("Time for FORWARD propagation: \t");timer.Print();timer.Start();
6312 // RemoveUsed(fSeeds,0.7,0.7,6);
6313 //RemoveOverlap(fSeeds,0.9,7,kTRUE);
6315 nseed=fSeeds->GetEntriesFast();
6317 for (Int_t i=0; i<nseed; i++) {
6318 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
6320 Int_t nc=t.GetNumberOfClusters();
6322 delete fSeeds->RemoveAt(i);
6325 t.CookdEdx(0.02,0.6);
6326 // CookLabel(pt,0.1); //For comparison only
6327 //if ((pt->IsActive() || (pt->fRemoval==10) )&& nc>50 &&pt->GetNumberOfClusters()>0.4*pt->fNFoundable){
6328 if ((pt->IsActive() || (pt->fRemoval==10) )){
6329 cerr<<found++<<'\r';
6332 delete fSeeds->RemoveAt(i);
6337 // fNTracks = found;
6339 Info("Clusters2Tracks","Time for overlap removal, track writing and dedx cooking: \t"); timer.Print();timer.Start();
6342 // cerr<<"Number of found tracks : "<<"\t"<<found<<endl;
6343 Info("Clusters2Tracks","Number of found tracks %d",found);
6345 // UnloadClusters();
6350 void AliTPCtrackerMI::Tracking(TObjArray * arr)
6353 // tracking of the seeds
6356 fSectors = fOuterSec;
6357 ParallelTracking(arr,150,63);
6358 fSectors = fOuterSec;
6359 ParallelTracking(arr,63,0);
6362 TObjArray * AliTPCtrackerMI::Tracking(Int_t seedtype, Int_t i1, Int_t i2, Float_t cuts[4], Float_t dy, Int_t dsec)
6367 TObjArray * arr = new TObjArray;
6369 fSectors = fOuterSec;
6372 for (Int_t sec=0;sec<fkNOS;sec++){
6373 if (seedtype==3) MakeSeeds3(arr,sec,i1,i2,cuts,dy, dsec);
6374 if (seedtype==4) MakeSeeds5(arr,sec,i1,i2,cuts,dy);
6375 if (seedtype==2) MakeSeeds2(arr,sec,i1,i2,cuts,dy);
6378 Info("Tracking","\nSeeding - %d\t%d\t%d\t%d\n",seedtype,i1,i2,arr->GetEntriesFast());
6390 TObjArray * AliTPCtrackerMI::Tracking()
6394 if (AliTPCReconstructor::GetRecoParam()->GetSpecialSeeding()) return TrackingSpecial();
6397 Int_t nup=fOuterSec->GetNRows()+fInnerSec->GetNRows();
6399 TObjArray * seeds = new TObjArray;
6401 Int_t fLastSeedRowSec=AliTPCReconstructor::GetRecoParam()->GetLastSeedRowSec();
6402 Int_t gapPrim = AliTPCReconstructor::GetRecoParam()->GetSeedGapPrim();
6403 Int_t gapSec = AliTPCReconstructor::GetRecoParam()->GetSeedGapSec();
6411 Float_t fnumber = 3.0;
6412 Float_t fdensity = 3.0;
6417 for (Int_t delta = 0; delta<18; delta+=gapPrim){
6421 arr = Tracking(3,nup-1-delta,nup-1-delta-gap,cuts,-1,1);
6422 SumTracks(seeds,arr);
6423 SignClusters(seeds,fnumber,fdensity);
6425 for (Int_t i=2;i<6;i+=2){
6426 // seed high pt tracks
6429 arr = Tracking(3,nup-i-delta,nup-i-delta-gap,cuts,-1,0);
6430 SumTracks(seeds,arr);
6431 SignClusters(seeds,fnumber,fdensity);
6436 // RemoveUsed(seeds,0.9,0.9,1);
6437 // UnsignClusters();
6438 // SignClusters(seeds,fnumber,fdensity);
6442 for (Int_t delta = 20; delta<120; delta+=gapPrim){
6444 // seed high pt tracks
6448 arr = Tracking(3,nup-delta,nup-delta-gap,cuts,-1);
6449 SumTracks(seeds,arr);
6450 SignClusters(seeds,fnumber,fdensity);
6455 arr = Tracking(3,nup-delta-5,nup-delta-5-gap,cuts,-1);
6456 SumTracks(seeds,arr);
6457 SignClusters(seeds,fnumber,fdensity);
6468 Info("Tracking()","\n\nPrimary seeding\t%d\n\n",seeds->GetEntriesFast());
6472 // RemoveUsed(seeds,0.75,0.75,1);
6474 //SignClusters(seeds,fnumber,fdensity);
6483 arr = Tracking(4,nup-1,nup-1-gap,cuts,-1);
6484 SumTracks(seeds,arr);
6485 SignClusters(seeds,fnumber,fdensity);
6487 arr = Tracking(4,nup-2,nup-2-gap,cuts,-1);
6488 SumTracks(seeds,arr);
6489 SignClusters(seeds,fnumber,fdensity);
6491 arr = Tracking(4,nup-3,nup-3-gap,cuts,-1);
6492 SumTracks(seeds,arr);
6493 SignClusters(seeds,fnumber,fdensity);
6495 arr = Tracking(4,nup-5,nup-5-gap,cuts,-1);
6496 SumTracks(seeds,arr);
6497 SignClusters(seeds,fnumber,fdensity);
6499 arr = Tracking(4,nup-7,nup-7-gap,cuts,-1);
6500 SumTracks(seeds,arr);
6501 SignClusters(seeds,fnumber,fdensity);
6504 arr = Tracking(4,nup-9,nup-9-gap,cuts,-1);
6505 SumTracks(seeds,arr);
6506 SignClusters(seeds,fnumber,fdensity);
6510 for (Int_t delta = 9; delta<30; delta+=gapSec){
6516 arr = Tracking(4,nup-1-delta,nup-1-delta-gap,cuts,-1);
6517 SumTracks(seeds,arr);
6518 SignClusters(seeds,fnumber,fdensity);
6520 arr = Tracking(4,nup-3-delta,nup-5-delta-gap,cuts,4);
6521 SumTracks(seeds,arr);
6522 SignClusters(seeds,fnumber,fdensity);
6535 for (Int_t delta = 30; delta<fLastSeedRowSec; delta+=gapSec){
6541 arr = Tracking(4,nup-1-delta,nup-1-delta-gap,cuts,-1);
6542 SumTracks(seeds,arr);
6543 SignClusters(seeds,fnumber,fdensity);
6545 arr = Tracking(4,nup-5-delta,nup-5-delta-gap,cuts,5 );
6546 SumTracks(seeds,arr);
6547 SignClusters(seeds,fnumber,fdensity);
6551 Info("Tracking()","\n\nSecondary seeding\t%d\n\n",seeds->GetEntriesFast());
6562 TObjArray * AliTPCtrackerMI::TrackingSpecial()
6565 // seeding adjusted for laser and cosmic tests - short tracks with big inclination angle
6566 // no primary vertex seeding tried
6570 Int_t nup=fOuterSec->GetNRows()+fInnerSec->GetNRows();
6572 TObjArray * seeds = new TObjArray;
6577 Float_t fnumber = 3.0;
6578 Float_t fdensity = 3.0;
6581 cuts[0] = AliTPCReconstructor::GetRecoParam()->GetMaxC(); // max curvature
6582 cuts[1] = 3.5; // max tan(phi) angle for seeding
6583 cuts[2] = 3.; // not used (cut on z primary vertex)
6584 cuts[3] = 3.5; // max tan(theta) angle for seeding
6586 for (Int_t delta = 0; nup-delta-gap-1>0; delta+=3){
6588 arr = Tracking(4,nup-1-delta,nup-1-delta-gap,cuts,-1);
6589 SumTracks(seeds,arr);
6590 SignClusters(seeds,fnumber,fdensity);
6594 Info("Tracking()","\n\nSecondary seeding\t%d\n\n",seeds->GetEntriesFast());
6605 void AliTPCtrackerMI::SumTracks(TObjArray *arr1,TObjArray *&arr2) const
6608 //sum tracks to common container
6609 //remove suspicious tracks
6610 Int_t nseed = arr2->GetEntriesFast();
6611 for (Int_t i=0;i<nseed;i++){
6612 AliTPCseed *pt=(AliTPCseed*)arr2->UncheckedAt(i);
6615 // remove tracks with too big curvature
6617 if (TMath::Abs(pt->GetC())>AliTPCReconstructor::GetRecoParam()->GetMaxC()){
6618 delete arr2->RemoveAt(i);
6621 // REMOVE VERY SHORT TRACKS
6622 if (pt->GetNumberOfClusters()<20){
6623 delete arr2->RemoveAt(i);
6626 // NORMAL ACTIVE TRACK
6627 if (pt->IsActive()){
6628 arr1->AddLast(arr2->RemoveAt(i));
6631 //remove not usable tracks
6632 if (pt->GetRemoval()!=10){
6633 delete arr2->RemoveAt(i);
6637 // ENABLE ONLY ENOUGH GOOD STOPPED TRACKS
6638 if (pt->GetDensityFirst(20)>0.8 || pt->GetDensityFirst(30)>0.8 || pt->GetDensityFirst(40)>0.7)
6639 arr1->AddLast(arr2->RemoveAt(i));
6641 delete arr2->RemoveAt(i);
6645 delete arr2; arr2 = 0;
6650 void AliTPCtrackerMI::ParallelTracking(TObjArray *const arr, Int_t rfirst, Int_t rlast)
6653 // try to track in parralel
6655 Int_t nseed=arr->GetEntriesFast();
6656 //prepare seeds for tracking
6657 for (Int_t i=0; i<nseed; i++) {
6658 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i), &t=*pt;
6660 if (!t.IsActive()) continue;
6661 // follow prolongation to the first layer
6662 if ( (fSectors ==fInnerSec) || (t.GetFirstPoint()-fkParam->GetNRowLow()>rfirst+1) )
6663 FollowProlongation(t, rfirst+1);
6668 for (Int_t nr=rfirst; nr>=rlast; nr--){
6669 if (nr<fInnerSec->GetNRows())
6670 fSectors = fInnerSec;
6672 fSectors = fOuterSec;
6673 // make indexes with the cluster tracks for given
6675 // find nearest cluster
6676 for (Int_t i=0; i<nseed; i++) {
6677 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i), &t=*pt;
6679 if (nr==80) pt->UpdateReference();
6680 if (!pt->IsActive()) continue;
6681 // if ( (fSectors ==fOuterSec) && (pt->fFirstPoint-fkParam->GetNRowLow())<nr) continue;
6682 if (pt->GetRelativeSector()>17) {
6685 UpdateClusters(t,nr);
6687 // prolonagate to the nearest cluster - if founded
6688 for (Int_t i=0; i<nseed; i++) {
6689 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
6691 if (!pt->IsActive()) continue;
6692 // if ((fSectors ==fOuterSec) && (pt->fFirstPoint-fkParam->GetNRowLow())<nr) continue;
6693 if (pt->GetRelativeSector()>17) {
6696 FollowToNextCluster(*pt,nr);
6701 void AliTPCtrackerMI::PrepareForBackProlongation(TObjArray *const arr,Float_t fac) const
6705 // if we use TPC track itself we have to "update" covariance
6707 Int_t nseed= arr->GetEntriesFast();
6708 for (Int_t i=0;i<nseed;i++){
6709 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
6713 //rotate to current local system at first accepted point
6714 Int_t index = pt->GetClusterIndex2(pt->GetFirstPoint());
6715 Int_t sec = (index&0xff000000)>>24;
6717 Float_t angle1 = fInnerSec->GetAlpha()*sec+fInnerSec->GetAlphaShift();
6718 if (angle1>TMath::Pi())
6719 angle1-=2.*TMath::Pi();
6720 Float_t angle2 = pt->GetAlpha();
6722 if (TMath::Abs(angle1-angle2)>0.001){
6723 pt->Rotate(angle1-angle2);
6724 //angle2 = pt->GetAlpha();
6725 //pt->fRelativeSector = pt->GetAlpha()/fInnerSec->GetAlpha();
6726 //if (pt->GetAlpha()<0)
6727 // pt->fRelativeSector+=18;
6728 //sec = pt->fRelativeSector;
6737 void AliTPCtrackerMI::PrepareForProlongation(TObjArray *const arr, Float_t fac) const
6741 // if we use TPC track itself we have to "update" covariance
6743 Int_t nseed= arr->GetEntriesFast();
6744 for (Int_t i=0;i<nseed;i++){
6745 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
6748 pt->SetFirstPoint(pt->GetLastPoint());
6756 Int_t AliTPCtrackerMI::PropagateBack(TObjArray *const arr)
6759 // make back propagation
6761 Int_t nseed= arr->GetEntriesFast();
6762 for (Int_t i=0;i<nseed;i++){
6763 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
6764 if (pt&& pt->GetKinkIndex(0)<=0) {
6765 //AliTPCseed *pt2 = new AliTPCseed(*pt);
6766 fSectors = fInnerSec;
6767 //FollowBackProlongation(*pt,fInnerSec->GetNRows()-1);
6768 //fSectors = fOuterSec;
6769 FollowBackProlongation(*pt,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1);
6770 //if (pt->GetNumberOfClusters()<(pt->fEsd->GetTPCclusters(0)) ){
6771 // Error("PropagateBack","Not prolonged track %d",pt->GetLabel());
6772 // FollowBackProlongation(*pt2,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1);
6775 if (pt&& pt->GetKinkIndex(0)>0) {
6776 AliESDkink * kink = fEvent->GetKink(pt->GetKinkIndex(0)-1);
6777 pt->SetFirstPoint(kink->GetTPCRow0());
6778 fSectors = fInnerSec;
6779 FollowBackProlongation(*pt,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1);
6787 Int_t AliTPCtrackerMI::PropagateForward2(TObjArray *const arr)
6790 // make forward propagation
6792 Int_t nseed= arr->GetEntriesFast();
6794 for (Int_t i=0;i<nseed;i++){
6795 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
6797 FollowProlongation(*pt,0);
6806 Int_t AliTPCtrackerMI::PropagateForward()
6809 // propagate track forward
6811 Int_t nseed = fSeeds->GetEntriesFast();
6812 for (Int_t i=0;i<nseed;i++){
6813 AliTPCseed *pt = (AliTPCseed*)fSeeds->UncheckedAt(i);
6815 AliTPCseed &t = *pt;
6816 Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
6817 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
6818 if (alpha < 0. ) alpha += 2.*TMath::Pi();
6819 t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN);
6823 fSectors = fOuterSec;
6824 ParallelTracking(fSeeds,fOuterSec->GetNRows()+fInnerSec->GetNRows()-1,fInnerSec->GetNRows());
6825 fSectors = fInnerSec;
6826 ParallelTracking(fSeeds,fInnerSec->GetNRows()-1,0);
6835 Int_t AliTPCtrackerMI::PropagateBack(AliTPCseed *const pt, Int_t row0, Int_t row1)
6838 // make back propagation, in between row0 and row1
6842 fSectors = fInnerSec;
6845 if (row1<fSectors->GetNRows())
6848 r1 = fSectors->GetNRows()-1;
6850 if (row0<fSectors->GetNRows()&& r1>0 )
6851 FollowBackProlongation(*pt,r1);
6852 if (row1<=fSectors->GetNRows())
6855 r1 = row1 - fSectors->GetNRows();
6856 if (r1<=0) return 0;
6857 if (r1>=fOuterSec->GetNRows()) return 0;
6858 fSectors = fOuterSec;
6859 return FollowBackProlongation(*pt,r1);
6867 void AliTPCtrackerMI::GetShape(AliTPCseed * seed, Int_t row)
6871 AliTPCClusterParam * clparam = AliTPCcalibDB::Instance()->GetClusterParam();
6872 Float_t zdrift = TMath::Abs((fkParam->GetZLength(0)-TMath::Abs(seed->GetZ())));
6873 Int_t type = (seed->GetSector() < fkParam->GetNSector()/2) ? 0: (row>126) ? 1:2;
6874 Double_t angulary = seed->GetSnp();
6875 angulary = angulary*angulary/((1.-angulary)*(1.+angulary));
6876 Double_t angularz = seed->GetTgl()*seed->GetTgl()*(1.+angulary);
6878 Double_t sigmay = clparam->GetRMS0(0,type,zdrift,TMath::Sqrt(TMath::Abs(angulary)));
6879 Double_t sigmaz = clparam->GetRMS0(1,type,zdrift,TMath::Sqrt(TMath::Abs(angularz)));
6880 seed->SetCurrentSigmaY2(sigmay*sigmay);
6881 seed->SetCurrentSigmaZ2(sigmaz*sigmaz);
6882 // Float_t sd2 = TMath::Abs((fkParam->GetZLength(0)-TMath::Abs(seed->GetZ())))*fkParam->GetDiffL()*fkParam->GetDiffL();
6883 // // Float_t padlength = fkParam->GetPadPitchLength(seed->fSector);
6884 // Float_t padlength = GetPadPitchLength(row);
6886 // Float_t sresy = (seed->GetSector() < fkParam->GetNSector()/2) ? 0.2 :0.3;
6887 // seed->SetCurrentSigmaY2(sd2+padlength*padlength*angulary/12.+sresy*sresy);
6889 // Float_t sresz = fkParam->GetZSigma();
6890 // seed->SetCurrentSigmaZ2(sd2+padlength*padlength*angularz*angularz*(1+angulary)/12.+sresz*sresz);
6892 Float_t wy = GetSigmaY(seed);
6893 Float_t wz = GetSigmaZ(seed);
6896 if (TMath::Abs(wy/seed->fCurrentSigmaY2-1)>0.0001 || TMath::Abs(wz/seed->fCurrentSigmaZ2-1)>0.0001 ){
6897 printf("problem\n");
6904 //__________________________________________________________________________
6905 void AliTPCtrackerMI::CookLabel(AliKalmanTrack *tk, Float_t wrong) const {
6906 //--------------------------------------------------------------------
6907 //This function "cooks" a track label. If label<0, this track is fake.
6908 //--------------------------------------------------------------------
6909 AliTPCseed * t = dynamic_cast<AliTPCseed*>(tk);
6911 printf("%s:%d wrong type \n",(char*)__FILE__,__LINE__);
6915 Int_t noc=t->GetNumberOfClusters();
6917 //printf("\nnot founded prolongation\n\n\n");
6923 AliTPCclusterMI *clusters[160];
6925 for (Int_t i=0;i<160;i++) {
6932 for (i=0; i<160 && current<noc; i++) {
6934 Int_t index=t->GetClusterIndex2(i);
6935 if (index<=0) continue;
6936 if (index&0x8000) continue;
6938 //clusters[current]=GetClusterMI(index);
6939 if (t->GetClusterPointer(i)){
6940 clusters[current]=t->GetClusterPointer(i);
6946 Int_t lab=123456789;
6947 for (i=0; i<noc; i++) {
6948 AliTPCclusterMI *c=clusters[i];
6950 lab=TMath::Abs(c->GetLabel(0));
6952 for (j=0; j<noc; j++) if (lb[j]==lab || mx[j]==0) break;
6958 for (i=0; i<noc; i++) if (mx[i]>max) {max=mx[i]; lab=lb[i];}
6960 for (i=0; i<noc; i++) {
6961 AliTPCclusterMI *c=clusters[i];
6963 if (TMath::Abs(c->GetLabel(1)) == lab ||
6964 TMath::Abs(c->GetLabel(2)) == lab ) max++;
6966 if (noc<=0) { lab=-1; return;}
6967 if ((1.- Float_t(max)/(noc)) > wrong) lab=-lab;
6970 Int_t tail=Int_t(0.10*noc);
6973 for (i=1; i<=160&&ind<tail; i++) {
6974 // AliTPCclusterMI *c=clusters[noc-i];
6975 AliTPCclusterMI *c=clusters[i];
6977 if (lab == TMath::Abs(c->GetLabel(0)) ||
6978 lab == TMath::Abs(c->GetLabel(1)) ||
6979 lab == TMath::Abs(c->GetLabel(2))) max++;
6982 if (max < Int_t(0.5*tail)) lab=-lab;
6989 //delete[] clusters;
6993 //__________________________________________________________________________
6994 Int_t AliTPCtrackerMI::CookLabel(AliTPCseed *const t, Float_t wrong,Int_t first, Int_t last) const {
6995 //--------------------------------------------------------------------
6996 //This function "cooks" a track label. If label<0, this track is fake.
6997 //--------------------------------------------------------------------
6998 Int_t noc=t->GetNumberOfClusters();
7000 //printf("\nnot founded prolongation\n\n\n");
7006 AliTPCclusterMI *clusters[160];
7008 for (Int_t i=0;i<160;i++) {
7015 for (i=0; i<160 && current<noc; i++) {
7016 if (i<first) continue;
7017 if (i>last) continue;
7018 Int_t index=t->GetClusterIndex2(i);
7019 if (index<=0) continue;
7020 if (index&0x8000) continue;
7022 //clusters[current]=GetClusterMI(index);
7023 if (t->GetClusterPointer(i)){
7024 clusters[current]=t->GetClusterPointer(i);
7029 if (noc<5) return -1;
7030 Int_t lab=123456789;
7031 for (i=0; i<noc; i++) {
7032 AliTPCclusterMI *c=clusters[i];
7034 lab=TMath::Abs(c->GetLabel(0));
7036 for (j=0; j<noc; j++) if (lb[j]==lab || mx[j]==0) break;
7042 for (i=0; i<noc; i++) if (mx[i]>max) {max=mx[i]; lab=lb[i];}
7044 for (i=0; i<noc; i++) {
7045 AliTPCclusterMI *c=clusters[i];
7047 if (TMath::Abs(c->GetLabel(1)) == lab ||
7048 TMath::Abs(c->GetLabel(2)) == lab ) max++;
7050 if (noc<=0) { lab=-1; return -1;}
7051 if ((1.- Float_t(max)/(noc)) > wrong) lab=-lab;
7054 Int_t tail=Int_t(0.10*noc);
7057 for (i=1; i<=160&&ind<tail; i++) {
7058 // AliTPCclusterMI *c=clusters[noc-i];
7059 AliTPCclusterMI *c=clusters[i];
7061 if (lab == TMath::Abs(c->GetLabel(0)) ||
7062 lab == TMath::Abs(c->GetLabel(1)) ||
7063 lab == TMath::Abs(c->GetLabel(2))) max++;
7066 if (max < Int_t(0.5*tail)) lab=-lab;
7069 // t->SetLabel(lab);
7073 //delete[] clusters;
7077 Int_t AliTPCtrackerMI::GetRowNumber(Double_t x[3]) const
7079 //return pad row number for given x vector
7080 Float_t phi = TMath::ATan2(x[1],x[0]);
7081 if(phi<0) phi=2.*TMath::Pi()+phi;
7082 // Get the local angle in the sector philoc
7083 const Float_t kRaddeg = 180/3.14159265358979312;
7084 Float_t phiangle = (Int_t (phi*kRaddeg/20.) + 0.5)*20./kRaddeg;
7085 Double_t localx = x[0]*TMath::Cos(phiangle)-x[1]*TMath::Sin(phiangle);
7086 return GetRowNumber(localx);
7091 void AliTPCtrackerMI::MakeBitmaps(AliTPCseed *t)
7093 //-----------------------------------------------------------------------
7094 // Fill the cluster and sharing bitmaps of the track
7095 //-----------------------------------------------------------------------
7097 Int_t firstpoint = 0;
7098 Int_t lastpoint = 159;
7099 AliTPCTrackerPoint *point;
7100 AliTPCclusterMI *cluster;
7102 for (int iter=firstpoint; iter<lastpoint; iter++) {
7103 // Change to cluster pointers to see if we have a cluster at given padrow
7104 cluster = t->GetClusterPointer(iter);
7106 t->SetClusterMapBit(iter, kTRUE);
7107 point = t->GetTrackPoint(iter);
7108 if (point->IsShared())
7109 t->SetSharedMapBit(iter,kTRUE);
7111 t->SetSharedMapBit(iter, kFALSE);
7114 t->SetClusterMapBit(iter, kFALSE);
7115 t->SetSharedMapBit(iter, kFALSE);
7120 Bool_t AliTPCtrackerMI::IsFindable(AliTPCseed & track){
7122 // return flag if there is findable cluster at given position
7125 Float_t z = track.GetZ();
7127 if (TMath::Abs(z)<(AliTPCReconstructor::GetCtgRange()*track.GetX()+kDeltaZ) &&
7128 TMath::Abs(z)<fkParam->GetZLength(0) &&
7129 (TMath::Abs(track.GetSnp())<AliTPCReconstructor::GetMaxSnpTracker()))
7135 void AliTPCtrackerMI::AddCovariance(AliTPCseed * seed){
7137 // Adding systematic error
7138 // !!!! the systematic error for element 4 is in 1/cm not in pt
7140 const Double_t *param = AliTPCReconstructor::GetRecoParam()->GetSystematicError();
7141 Double_t *covarIn= (Double_t*)seed->GetCovariance();
7143 for (Int_t i=0;i<15;i++) covar[i]=0;
7149 covar[0] = param[0]*param[0];
7150 covar[2] = param[1]*param[1];
7151 covar[5] = param[2]*param[2];
7152 covar[9] = param[3]*param[3];
7153 Double_t facC = AliTracker::GetBz()*kB2C;
7154 covar[14]= param[4]*param[4]*facC*facC;
7156 covar[1]=TMath::Sqrt((covar[0]*covar[2]))*covarIn[1]/TMath::Sqrt((covarIn[0]*covarIn[2]));
7158 covar[3]=TMath::Sqrt((covar[0]*covar[5]))*covarIn[3]/TMath::Sqrt((covarIn[0]*covarIn[5]));
7159 covar[4]=TMath::Sqrt((covar[2]*covar[5]))*covarIn[4]/TMath::Sqrt((covarIn[2]*covarIn[5]));
7161 covar[6]=TMath::Sqrt((covar[0]*covar[9]))*covarIn[6]/TMath::Sqrt((covarIn[0]*covarIn[9]));
7162 covar[7]=TMath::Sqrt((covar[2]*covar[9]))*covarIn[7]/TMath::Sqrt((covarIn[2]*covarIn[9]));
7163 covar[8]=TMath::Sqrt((covar[5]*covar[9]))*covarIn[8]/TMath::Sqrt((covarIn[5]*covarIn[9]));
7165 covar[10]=TMath::Sqrt((covar[0]*covar[14]))*covarIn[10]/TMath::Sqrt((covarIn[0]*covarIn[14]));
7166 covar[11]=TMath::Sqrt((covar[2]*covar[14]))*covarIn[11]/TMath::Sqrt((covarIn[2]*covarIn[14]));
7167 covar[12]=TMath::Sqrt((covar[5]*covar[14]))*covarIn[12]/TMath::Sqrt((covarIn[5]*covarIn[14]));
7168 covar[13]=TMath::Sqrt((covar[9]*covar[14]))*covarIn[13]/TMath::Sqrt((covarIn[9]*covarIn[14]));
7170 seed->AddCovariance(covar);