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>
116 #include <TGraphErrors.h>
118 #include "AliComplexCluster.h"
119 #include "AliESDEvent.h"
120 #include "AliESDtrack.h"
121 #include "AliESDVertex.h"
124 #include "AliHelix.h"
125 #include "AliRunLoader.h"
126 #include "AliTPCClustersRow.h"
127 #include "AliTPCParam.h"
128 #include "AliTPCReconstructor.h"
129 #include "AliTPCpolyTrack.h"
130 #include "AliTPCreco.h"
131 #include "AliTPCseed.h"
133 #include "AliTPCtrackerSector.h"
134 #include "AliTPCtrackerMI.h"
135 #include "TStopwatch.h"
136 #include "AliTPCReconstructor.h"
137 #include "AliAlignObj.h"
138 #include "AliTrackPointArray.h"
140 #include "AliTPCcalibDB.h"
141 #include "AliTPCcalibDButil.h"
142 #include "AliTPCTransform.h"
143 #include "AliTPCClusterParam.h"
144 #include "AliTPCdEdxInfo.h"
145 #include "AliDCSSensorArray.h"
146 #include "AliDCSSensor.h"
151 ClassImp(AliTPCtrackerMI)
155 class AliTPCFastMath {
158 static Double_t FastAsin(Double_t x);
160 static Double_t fgFastAsin[20000]; //lookup table for fast asin computation
163 Double_t AliTPCFastMath::fgFastAsin[20000];
164 AliTPCFastMath gAliTPCFastMath; // needed to fill the LUT
166 AliTPCFastMath::AliTPCFastMath(){
168 // initialized lookup table;
169 for (Int_t i=0;i<10000;i++){
170 fgFastAsin[2*i] = TMath::ASin(i/10000.);
171 fgFastAsin[2*i+1] = (TMath::ASin((i+1)/10000.)-fgFastAsin[2*i]);
175 Double_t AliTPCFastMath::FastAsin(Double_t x){
177 // return asin using lookup table
179 Int_t index = int(x*10000);
180 return fgFastAsin[2*index]+(x*10000.-index)*fgFastAsin[2*index+1];
183 Int_t index = int(x*10000);
184 return -(fgFastAsin[2*index]+(x*10000.-index)*fgFastAsin[2*index+1]);
186 //__________________________________________________________________
187 AliTPCtrackerMI::AliTPCtrackerMI()
214 // default constructor
216 for (Int_t irow=0; irow<200; irow++){
223 //_____________________________________________________________________
227 Int_t AliTPCtrackerMI::UpdateTrack(AliTPCseed * track, Int_t accept){
229 //update track information using current cluster - track->fCurrentCluster
232 AliTPCclusterMI* c =track->GetCurrentCluster();
233 if (accept > 0) //sign not accepted clusters
234 track->SetCurrentClusterIndex1(track->GetCurrentClusterIndex1() | 0x8000);
235 else // unsign accpeted clusters
236 track->SetCurrentClusterIndex1(track->GetCurrentClusterIndex1() & 0xffff7fff);
237 UInt_t i = track->GetCurrentClusterIndex1();
239 Int_t sec=(i&0xff000000)>>24;
240 //Int_t row = (i&0x00ff0000)>>16;
241 track->SetRow((i&0x00ff0000)>>16);
242 track->SetSector(sec);
243 // Int_t index = i&0xFFFF;
244 if (sec>=fkParam->GetNInnerSector()) track->SetRow(track->GetRow()+fkParam->GetNRowLow());
245 track->SetClusterIndex2(track->GetRow(), i);
246 //track->fFirstPoint = row;
247 //if ( track->fLastPoint<row) track->fLastPoint =row;
248 // if (track->fRow<0 || track->fRow>160) {
249 // printf("problem\n");
251 if (track->GetFirstPoint()>track->GetRow())
252 track->SetFirstPoint(track->GetRow());
253 if (track->GetLastPoint()<track->GetRow())
254 track->SetLastPoint(track->GetRow());
257 track->SetClusterPointer(track->GetRow(),c);
260 Double_t angle2 = track->GetSnp()*track->GetSnp();
262 //SET NEW Track Point
264 if (angle2<1) //PH sometimes angle2 is very big. To be investigated...
266 angle2 = TMath::Sqrt(angle2/(1-angle2));
267 AliTPCTrackerPoint &point =*(track->GetTrackPoint(track->GetRow()));
269 point.SetSigmaY(c->GetSigmaY2()/track->GetCurrentSigmaY2());
270 point.SetSigmaZ(c->GetSigmaZ2()/track->GetCurrentSigmaZ2());
271 point.SetErrY(sqrt(track->GetErrorY2()));
272 point.SetErrZ(sqrt(track->GetErrorZ2()));
274 point.SetX(track->GetX());
275 point.SetY(track->GetY());
276 point.SetZ(track->GetZ());
277 point.SetAngleY(angle2);
278 point.SetAngleZ(track->GetTgl());
279 if (point.IsShared()){
280 track->SetErrorY2(track->GetErrorY2()*4);
281 track->SetErrorZ2(track->GetErrorZ2()*4);
285 Double_t chi2 = track->GetPredictedChi2(track->GetCurrentCluster());
287 // track->SetErrorY2(track->GetErrorY2()*1.3);
288 // track->SetErrorY2(track->GetErrorY2()+0.01);
289 // track->SetErrorZ2(track->GetErrorZ2()*1.3);
290 // track->SetErrorZ2(track->GetErrorZ2()+0.005);
292 if (accept>0) return 0;
293 if (track->GetNumberOfClusters()%20==0){
294 // if (track->fHelixIn){
295 // TClonesArray & larr = *(track->fHelixIn);
296 // Int_t ihelix = larr.GetEntriesFast();
297 // new(larr[ihelix]) AliHelix(*track) ;
300 track->SetNoCluster(0);
301 return track->Update(c,chi2,i);
306 Int_t AliTPCtrackerMI::AcceptCluster(AliTPCseed * seed, AliTPCclusterMI * cluster)
309 // decide according desired precision to accept given
310 // cluster for tracking
312 seed->GetProlongation(cluster->GetX(),yt,zt);
313 Double_t sy2=ErrY2(seed,cluster);
314 Double_t sz2=ErrZ2(seed,cluster);
316 Double_t sdistancey2 = sy2+seed->GetSigmaY2();
317 Double_t sdistancez2 = sz2+seed->GetSigmaZ2();
318 Double_t dy=seed->GetCurrentCluster()->GetY()-yt;
319 Double_t dz=seed->GetCurrentCluster()->GetZ()-zt;
320 Double_t rdistancey2 = (seed->GetCurrentCluster()->GetY()-yt)*
321 (seed->GetCurrentCluster()->GetY()-yt)/sdistancey2;
322 Double_t rdistancez2 = (seed->GetCurrentCluster()->GetZ()-zt)*
323 (seed->GetCurrentCluster()->GetZ()-zt)/sdistancez2;
325 Double_t rdistance2 = rdistancey2+rdistancez2;
328 if (AliTPCReconstructor::StreamLevel()>2 && seed->GetNumberOfClusters()>20) {
329 Float_t rmsy2 = seed->GetCurrentSigmaY2();
330 Float_t rmsz2 = seed->GetCurrentSigmaZ2();
331 Float_t rmsy2p30 = seed->GetCMeanSigmaY2p30();
332 Float_t rmsz2p30 = seed->GetCMeanSigmaZ2p30();
333 Float_t rmsy2p30R = seed->GetCMeanSigmaY2p30R();
334 Float_t rmsz2p30R = seed->GetCMeanSigmaZ2p30R();
335 AliExternalTrackParam param(*seed);
336 static TVectorD gcl(3),gtr(3);
338 param.GetXYZ(gcl.GetMatrixArray());
339 cluster->GetGlobalXYZ(gclf);
340 gcl[0]=gclf[0]; gcl[1]=gclf[1]; gcl[2]=gclf[2];
343 if (AliTPCReconstructor::StreamLevel()>2) {
344 (*fDebugStreamer)<<"ErrParam"<<
357 "rmsy2p30="<<rmsy2p30<<
358 "rmsz2p30="<<rmsz2p30<<
359 "rmsy2p30R="<<rmsy2p30R<<
360 "rmsz2p30R="<<rmsz2p30R<<
361 // normalize distance -
362 "rdisty="<<rdistancey2<<
363 "rdistz="<<rdistancez2<<
364 "rdist="<<rdistance2<< //
368 //return 0; // temporary
369 if (rdistance2>32) return 3;
372 if ((rdistancey2>9. || rdistancez2>9.) && cluster->GetType()==0)
373 return 2; //suspisiouce - will be changed
375 if ((rdistancey2>6.25 || rdistancez2>6.25) && cluster->GetType()>0)
376 // strict cut on overlaped cluster
377 return 2; //suspisiouce - will be changed
379 if ( (rdistancey2>1. || rdistancez2>6.25 )
380 && cluster->GetType()<0){
381 seed->SetNFoundable(seed->GetNFoundable()-1);
385 if (fUseHLTClusters == 3 || fUseHLTClusters == 4) {
387 if(!AliTPCReconstructor::GetRecoParam()->GetUseHLTOnePadCluster()) {
388 if (TMath::Abs(cluster->GetSigmaY2()) < kAlmost0)
401 //_____________________________________________________________________________
402 AliTPCtrackerMI::AliTPCtrackerMI(const AliTPCParam *par):
404 fkNIS(par->GetNInnerSector()/2),
406 fkNOS(par->GetNOuterSector()/2),
428 //---------------------------------------------------------------------
429 // The main TPC tracker constructor
430 //---------------------------------------------------------------------
431 fInnerSec=new AliTPCtrackerSector[fkNIS];
432 fOuterSec=new AliTPCtrackerSector[fkNOS];
435 for (i=0; i<fkNIS; i++) fInnerSec[i].Setup(par,0);
436 for (i=0; i<fkNOS; i++) fOuterSec[i].Setup(par,1);
439 Int_t nrowlow = par->GetNRowLow();
440 Int_t nrowup = par->GetNRowUp();
443 for (i=0;i<nrowlow;i++){
444 fXRow[i] = par->GetPadRowRadiiLow(i);
445 fPadLength[i]= par->GetPadPitchLength(0,i);
446 fYMax[i] = fXRow[i]*TMath::Tan(0.5*par->GetInnerAngle());
450 for (i=0;i<nrowup;i++){
451 fXRow[i+nrowlow] = par->GetPadRowRadiiUp(i);
452 fPadLength[i+nrowlow] = par->GetPadPitchLength(60,i);
453 fYMax[i+nrowlow] = fXRow[i+nrowlow]*TMath::Tan(0.5*par->GetOuterAngle());
456 if (AliTPCReconstructor::StreamLevel()>0) {
457 fDebugStreamer = new TTreeSRedirector("TPCdebug.root");
460 fSeedsPool = new TClonesArray("AliTPCseed",1000);
462 //________________________________________________________________________
463 AliTPCtrackerMI::AliTPCtrackerMI(const AliTPCtrackerMI &t):
489 //------------------------------------
490 // dummy copy constructor
491 //------------------------------------------------------------------
493 for (Int_t irow=0; irow<200; irow++){
500 AliTPCtrackerMI & AliTPCtrackerMI::operator=(const AliTPCtrackerMI& /*r*/)
502 //------------------------------
504 //--------------------------------------------------------------
507 //_____________________________________________________________________________
508 AliTPCtrackerMI::~AliTPCtrackerMI() {
509 //------------------------------------------------------------------
510 // TPC tracker destructor
511 //------------------------------------------------------------------
518 if (fDebugStreamer) delete fDebugStreamer;
519 if (fSeedsPool) delete fSeedsPool;
523 void AliTPCtrackerMI::FillESD(const TObjArray* arr)
527 //fill esds using updated tracks
530 // write tracks to the event
531 // store index of the track
532 Int_t nseed=arr->GetEntriesFast();
533 //FindKinks(arr,fEvent);
534 for (Int_t i=0; i<nseed; i++) {
535 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
539 if (AliTPCReconstructor::StreamLevel()>1) {
540 (*fDebugStreamer)<<"Track0"<<
544 // pt->PropagateTo(fkParam->GetInnerRadiusLow());
545 if (pt->GetKinkIndex(0)<=0){ //don't propagate daughter tracks
546 pt->PropagateTo(fkParam->GetInnerRadiusLow());
549 if (( pt->GetPoints()[2]- pt->GetPoints()[0])>5 && pt->GetPoints()[3]>0.8){
551 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
552 iotrack.SetTPCPoints(pt->GetPoints());
553 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
554 iotrack.SetV0Indexes(pt->GetV0Indexes());
555 // iotrack.SetTPCpid(pt->fTPCr);
556 //iotrack.SetTPCindex(i);
557 fEvent->AddTrack(&iotrack);
561 if ( (pt->GetNumberOfClusters()>70)&& (Float_t(pt->GetNumberOfClusters())/Float_t(pt->GetNFoundable()))>0.55) {
563 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
564 iotrack.SetTPCPoints(pt->GetPoints());
565 //iotrack.SetTPCindex(i);
566 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
567 iotrack.SetV0Indexes(pt->GetV0Indexes());
568 // iotrack.SetTPCpid(pt->fTPCr);
569 fEvent->AddTrack(&iotrack);
573 // short tracks - maybe decays
575 if ( (pt->GetNumberOfClusters()>30) && (Float_t(pt->GetNumberOfClusters())/Float_t(pt->GetNFoundable()))>0.70) {
576 Int_t found,foundable,shared;
577 pt->GetClusterStatistic(0,60,found, foundable,shared,kFALSE);
578 if ( (found>20) && (pt->GetNShared()/float(pt->GetNumberOfClusters())<0.2)){
580 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
581 //iotrack.SetTPCindex(i);
582 iotrack.SetTPCPoints(pt->GetPoints());
583 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
584 iotrack.SetV0Indexes(pt->GetV0Indexes());
585 //iotrack.SetTPCpid(pt->fTPCr);
586 fEvent->AddTrack(&iotrack);
591 if ( (pt->GetNumberOfClusters()>20) && (Float_t(pt->GetNumberOfClusters())/Float_t(pt->GetNFoundable()))>0.8) {
592 Int_t found,foundable,shared;
593 pt->GetClusterStatistic(0,60,found, foundable,shared,kFALSE);
594 if (found<20) continue;
595 if (pt->GetNShared()/float(pt->GetNumberOfClusters())>0.2) continue;
598 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
599 iotrack.SetTPCPoints(pt->GetPoints());
600 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
601 iotrack.SetV0Indexes(pt->GetV0Indexes());
602 //iotrack.SetTPCpid(pt->fTPCr);
603 //iotrack.SetTPCindex(i);
604 fEvent->AddTrack(&iotrack);
607 // short tracks - secondaties
609 if ( (pt->GetNumberOfClusters()>30) ) {
610 Int_t found,foundable,shared;
611 pt->GetClusterStatistic(128,158,found, foundable,shared,kFALSE);
612 if ( (found>20) && (pt->GetNShared()/float(pt->GetNumberOfClusters())<0.2) &&float(found)/float(foundable)>0.8){
614 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
615 iotrack.SetTPCPoints(pt->GetPoints());
616 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
617 iotrack.SetV0Indexes(pt->GetV0Indexes());
618 //iotrack.SetTPCpid(pt->fTPCr);
619 //iotrack.SetTPCindex(i);
620 fEvent->AddTrack(&iotrack);
625 if ( (pt->GetNumberOfClusters()>15)) {
626 Int_t found,foundable,shared;
627 pt->GetClusterStatistic(138,158,found, foundable,shared,kFALSE);
628 if (found<15) continue;
629 if (foundable<=0) continue;
630 if (pt->GetNShared()/float(pt->GetNumberOfClusters())>0.2) continue;
631 if (float(found)/float(foundable)<0.8) continue;
634 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
635 iotrack.SetTPCPoints(pt->GetPoints());
636 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
637 iotrack.SetV0Indexes(pt->GetV0Indexes());
638 // iotrack.SetTPCpid(pt->fTPCr);
639 //iotrack.SetTPCindex(i);
640 fEvent->AddTrack(&iotrack);
644 // >> account for suppressed tracks in the kink indices (RS)
645 int nESDtracks = fEvent->GetNumberOfTracks();
646 for (int it=nESDtracks;it--;) {
647 AliESDtrack* esdTr = fEvent->GetTrack(it);
648 if (!esdTr || !esdTr->GetKinkIndex(0)) continue;
649 for (int ik=0;ik<3;ik++) {
651 if (!(knkId=esdTr->GetKinkIndex(ik))) break; // no more kinks for this track
652 AliESDkink* kink = fEvent->GetKink(TMath::Abs(knkId)-1);
654 AliError(Form("ESDTrack%d refers to non-existing kink %d",it,TMath::Abs(knkId)-1));
657 kink->SetIndex(it, knkId<0 ? 0:1); // update track index of the kink: mother at 0, daughter at 1
660 // << account for suppressed tracks in the kink indices (RS)
661 AliInfo(Form("Number of filled ESDs-\t%d\n",fEvent->GetNumberOfTracks()));
669 Double_t AliTPCtrackerMI::ErrY2(AliTPCseed* seed, const AliTPCclusterMI * cl){
672 // Use calibrated cluster error from OCDB
674 AliTPCClusterParam * clparam = AliTPCcalibDB::Instance()->GetClusterParam();
676 Float_t z = TMath::Abs(fkParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
677 Int_t ctype = cl->GetType();
678 Int_t type = (cl->GetRow()<63) ? 0: (cl->GetRow()>126) ? 1:2;
679 Double_t angle = seed->GetSnp()*seed->GetSnp();
680 angle = TMath::Sqrt(TMath::Abs(angle/(1.-angle)));
681 Double_t erry2 = clparam->GetError0Par(0,type, z,angle);
683 erry2+=0.5; // edge cluster
686 seed->SetErrorY2(erry2);
690 //calculate look-up table at the beginning
691 // static Bool_t ginit = kFALSE;
692 // static Float_t gnoise1,gnoise2,gnoise3;
693 // static Float_t ggg1[10000];
694 // static Float_t ggg2[10000];
695 // static Float_t ggg3[10000];
696 // static Float_t glandau1[10000];
697 // static Float_t glandau2[10000];
698 // static Float_t glandau3[10000];
700 // static Float_t gcor01[500];
701 // static Float_t gcor02[500];
702 // static Float_t gcorp[500];
706 // if (ginit==kFALSE){
707 // for (Int_t i=1;i<500;i++){
708 // Float_t rsigma = float(i)/100.;
709 // gcor02[i] = TMath::Max(0.78 +TMath::Exp(7.4*(rsigma-1.2)),0.6);
710 // gcor01[i] = TMath::Max(0.72 +TMath::Exp(3.36*(rsigma-1.2)),0.6);
711 // gcorp[i] = TMath::Max(TMath::Power((rsigma+0.5),1.5),1.2);
715 // for (Int_t i=3;i<10000;i++){
719 // Float_t amp = float(i);
720 // Float_t padlength =0.75;
721 // gnoise1 = 0.0004/padlength;
722 // Float_t nel = 0.268*amp;
723 // Float_t nprim = 0.155*amp;
724 // ggg1[i] = fkParam->GetDiffT()*fkParam->GetDiffT()*(2+0.001*nel/(padlength*padlength))/nel;
725 // glandau1[i] = (2.+0.12*nprim)*0.5* (2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
726 // if (glandau1[i]>1) glandau1[i]=1;
727 // glandau1[i]*=padlength*padlength/12.;
731 // gnoise2 = 0.0004/padlength;
733 // nprim = 0.133*amp;
734 // ggg2[i] = fkParam->GetDiffT()*fkParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
735 // glandau2[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
736 // if (glandau2[i]>1) glandau2[i]=1;
737 // glandau2[i]*=padlength*padlength/12.;
742 // gnoise3 = 0.0004/padlength;
744 // nprim = 0.133*amp;
745 // ggg3[i] = fkParam->GetDiffT()*fkParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
746 // glandau3[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
747 // if (glandau3[i]>1) glandau3[i]=1;
748 // glandau3[i]*=padlength*padlength/12.;
756 // Int_t amp = int(TMath::Abs(cl->GetQ()));
758 // seed->SetErrorY2(1.);
762 // Float_t z = TMath::Abs(fkParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
763 // Int_t ctype = cl->GetType();
764 // Float_t padlength= GetPadPitchLength(seed->GetRow());
765 // Double_t angle2 = seed->GetSnp()*seed->GetSnp();
766 // angle2 = angle2/(1-angle2);
768 // //cluster "quality"
769 // Int_t rsigmay = int(100.*cl->GetSigmaY2()/(seed->GetCurrentSigmaY2()));
772 // if (fSectors==fInnerSec){
773 // snoise2 = gnoise1;
774 // res = ggg1[amp]*z+glandau1[amp]*angle2;
775 // if (ctype==0) res *= gcor01[rsigmay];
778 // res*= gcorp[rsigmay];
782 // if (padlength<1.1){
783 // snoise2 = gnoise2;
784 // res = ggg2[amp]*z+glandau2[amp]*angle2;
785 // if (ctype==0) res *= gcor02[rsigmay];
788 // res*= gcorp[rsigmay];
792 // snoise2 = gnoise3;
793 // res = ggg3[amp]*z+glandau3[amp]*angle2;
794 // if (ctype==0) res *= gcor02[rsigmay];
797 // res*= gcorp[rsigmay];
804 // res*=2.4; // overestimate error 2 times
808 // if (res<2*snoise2)
811 // seed->SetErrorY2(res);
819 Double_t AliTPCtrackerMI::ErrZ2(AliTPCseed* seed, const AliTPCclusterMI * cl){
822 // Use calibrated cluster error from OCDB
824 AliTPCClusterParam * clparam = AliTPCcalibDB::Instance()->GetClusterParam();
826 Float_t z = TMath::Abs(fkParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
827 Int_t ctype = cl->GetType();
828 Int_t type = (cl->GetRow()<63) ? 0: (cl->GetRow()>126) ? 1:2;
830 Double_t angle2 = seed->GetSnp()*seed->GetSnp();
831 angle2 = seed->GetTgl()*seed->GetTgl()*(1+angle2/(1-angle2));
832 Double_t angle = TMath::Sqrt(TMath::Abs(angle2));
833 Double_t errz2 = clparam->GetError0Par(1,type, z,angle);
835 errz2+=0.5; // edge cluster
838 seed->SetErrorZ2(errz2);
844 // //seed->SetErrorY2(0.1);
846 // //calculate look-up table at the beginning
847 // static Bool_t ginit = kFALSE;
848 // static Float_t gnoise1,gnoise2,gnoise3;
849 // static Float_t ggg1[10000];
850 // static Float_t ggg2[10000];
851 // static Float_t ggg3[10000];
852 // static Float_t glandau1[10000];
853 // static Float_t glandau2[10000];
854 // static Float_t glandau3[10000];
856 // static Float_t gcor01[1000];
857 // static Float_t gcor02[1000];
858 // static Float_t gcorp[1000];
862 // if (ginit==kFALSE){
863 // for (Int_t i=1;i<1000;i++){
864 // Float_t rsigma = float(i)/100.;
865 // gcor02[i] = TMath::Max(0.81 +TMath::Exp(6.8*(rsigma-1.2)),0.6);
866 // gcor01[i] = TMath::Max(0.72 +TMath::Exp(2.04*(rsigma-1.2)),0.6);
867 // gcorp[i] = TMath::Max(TMath::Power((rsigma+0.5),1.5),1.2);
871 // for (Int_t i=3;i<10000;i++){
875 // Float_t amp = float(i);
876 // Float_t padlength =0.75;
877 // gnoise1 = 0.0004/padlength;
878 // Float_t nel = 0.268*amp;
879 // Float_t nprim = 0.155*amp;
880 // ggg1[i] = fkParam->GetDiffT()*fkParam->GetDiffT()*(2+0.001*nel/(padlength*padlength))/nel;
881 // glandau1[i] = (2.+0.12*nprim)*0.5* (2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
882 // if (glandau1[i]>1) glandau1[i]=1;
883 // glandau1[i]*=padlength*padlength/12.;
887 // gnoise2 = 0.0004/padlength;
889 // nprim = 0.133*amp;
890 // ggg2[i] = fkParam->GetDiffT()*fkParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
891 // glandau2[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
892 // if (glandau2[i]>1) glandau2[i]=1;
893 // glandau2[i]*=padlength*padlength/12.;
898 // gnoise3 = 0.0004/padlength;
900 // nprim = 0.133*amp;
901 // ggg3[i] = fkParam->GetDiffT()*fkParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
902 // glandau3[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
903 // if (glandau3[i]>1) glandau3[i]=1;
904 // glandau3[i]*=padlength*padlength/12.;
912 // Int_t amp = int(TMath::Abs(cl->GetQ()));
914 // seed->SetErrorY2(1.);
918 // Float_t z = TMath::Abs(fkParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
919 // Int_t ctype = cl->GetType();
920 // Float_t padlength= GetPadPitchLength(seed->GetRow());
922 // Double_t angle2 = seed->GetSnp()*seed->GetSnp();
923 // // if (angle2<0.6) angle2 = 0.6;
924 // angle2 = seed->GetTgl()*seed->GetTgl()*(1+angle2/(1-angle2));
926 // //cluster "quality"
927 // Int_t rsigmaz = int(100.*cl->GetSigmaZ2()/(seed->GetCurrentSigmaZ2()));
930 // if (fSectors==fInnerSec){
931 // snoise2 = gnoise1;
932 // res = ggg1[amp]*z+glandau1[amp]*angle2;
933 // if (ctype==0) res *= gcor01[rsigmaz];
936 // res*= gcorp[rsigmaz];
940 // if (padlength<1.1){
941 // snoise2 = gnoise2;
942 // res = ggg2[amp]*z+glandau2[amp]*angle2;
943 // if (ctype==0) res *= gcor02[rsigmaz];
946 // res*= gcorp[rsigmaz];
950 // snoise2 = gnoise3;
951 // res = ggg3[amp]*z+glandau3[amp]*angle2;
952 // if (ctype==0) res *= gcor02[rsigmaz];
955 // res*= gcorp[rsigmaz];
964 // if ((ctype<0) &&<70){
969 // if (res<2*snoise2)
971 // if (res>3) res =3;
972 // seed->SetErrorZ2(res);
980 void AliTPCtrackerMI::RotateToLocal(AliTPCseed *seed)
982 //rotate to track "local coordinata
983 Float_t x = seed->GetX();
984 Float_t y = seed->GetY();
985 Float_t ymax = x*TMath::Tan(0.5*fSectors->GetAlpha());
988 seed->SetRelativeSector((seed->GetRelativeSector()+1) % fN);
989 if (!seed->Rotate(fSectors->GetAlpha()))
991 } else if (y <-ymax) {
992 seed->SetRelativeSector((seed->GetRelativeSector()-1+fN) % fN);
993 if (!seed->Rotate(-fSectors->GetAlpha()))
1001 //_____________________________________________________________________________
1002 Double_t AliTPCtrackerMI::F1old(Double_t x1,Double_t y1,
1003 Double_t x2,Double_t y2,
1004 Double_t x3,Double_t y3) const
1006 //-----------------------------------------------------------------
1007 // Initial approximation of the track curvature
1008 //-----------------------------------------------------------------
1009 Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
1010 Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
1011 (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
1012 Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
1013 (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
1015 Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
1016 if ( xr*xr+yr*yr<=0.00000000000001) return 100;
1017 return -xr*yr/sqrt(xr*xr+yr*yr);
1022 //_____________________________________________________________________________
1023 Double_t AliTPCtrackerMI::F1(Double_t x1,Double_t y1,
1024 Double_t x2,Double_t y2,
1025 Double_t x3,Double_t y3) const
1027 //-----------------------------------------------------------------
1028 // Initial approximation of the track curvature
1029 //-----------------------------------------------------------------
1035 Double_t det = x3*y2-x2*y3;
1036 if (TMath::Abs(det)<1e-10){
1040 Double_t u = 0.5* (x2*(x2-x3)+y2*(y2-y3))/det;
1041 Double_t x0 = x3*0.5-y3*u;
1042 Double_t y0 = y3*0.5+x3*u;
1043 Double_t c2 = 1/TMath::Sqrt(x0*x0+y0*y0);
1049 Double_t AliTPCtrackerMI::F2(Double_t x1,Double_t y1,
1050 Double_t x2,Double_t y2,
1051 Double_t x3,Double_t y3) const
1053 //-----------------------------------------------------------------
1054 // Initial approximation of the track curvature
1055 //-----------------------------------------------------------------
1061 Double_t det = x3*y2-x2*y3;
1062 if (TMath::Abs(det)<1e-10) {
1066 Double_t u = 0.5* (x2*(x2-x3)+y2*(y2-y3))/det;
1067 Double_t x0 = x3*0.5-y3*u;
1068 Double_t y0 = y3*0.5+x3*u;
1069 Double_t c2 = 1/TMath::Sqrt(x0*x0+y0*y0);
1078 //_____________________________________________________________________________
1079 Double_t AliTPCtrackerMI::F2old(Double_t x1,Double_t y1,
1080 Double_t x2,Double_t y2,
1081 Double_t x3,Double_t y3) const
1083 //-----------------------------------------------------------------
1084 // Initial approximation of the track curvature times center of curvature
1085 //-----------------------------------------------------------------
1086 Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
1087 Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
1088 (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
1089 Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
1090 (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
1092 Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
1094 return -a/(d*y1-b)*xr/sqrt(xr*xr+yr*yr);
1097 //_____________________________________________________________________________
1098 Double_t AliTPCtrackerMI::F3(Double_t x1,Double_t y1,
1099 Double_t x2,Double_t y2,
1100 Double_t z1,Double_t z2) const
1102 //-----------------------------------------------------------------
1103 // Initial approximation of the tangent of the track dip angle
1104 //-----------------------------------------------------------------
1105 return (z1 - z2)/sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
1109 Double_t AliTPCtrackerMI::F3n(Double_t x1,Double_t y1,
1110 Double_t x2,Double_t y2,
1111 Double_t z1,Double_t z2, Double_t c) const
1113 //-----------------------------------------------------------------
1114 // Initial approximation of the tangent of the track dip angle
1115 //-----------------------------------------------------------------
1119 //angle1 = (z1-z2)*c/(TMath::ASin(c*x1-ni)-TMath::ASin(c*x2-ni));
1121 Double_t d = TMath::Sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
1122 if (TMath::Abs(d*c*0.5)>1) return 0;
1123 // Double_t angle2 = TMath::ASin(d*c*0.5);
1124 // Double_t angle2 = AliTPCFastMath::FastAsin(d*c*0.5);
1125 Double_t angle2 = (d*c*0.5>0.1)? TMath::ASin(d*c*0.5): AliTPCFastMath::FastAsin(d*c*0.5);
1127 angle2 = (z1-z2)*c/(angle2*2.);
1131 Bool_t AliTPCtrackerMI::GetProlongation(Double_t x1, Double_t x2, Double_t x[5], Double_t &y, Double_t &z) const
1132 {//-----------------------------------------------------------------
1133 // This function find proloncation of a track to a reference plane x=x2.
1134 //-----------------------------------------------------------------
1138 if (TMath::Abs(x[4]*x1 - x[2]) >= 0.999) {
1142 Double_t c1=x[4]*x1 - x[2], r1=TMath::Sqrt((1.-c1)*(1.+c1));
1143 Double_t c2=x[4]*x2 - x[2], r2=TMath::Sqrt((1.-c2)*(1.+c2));
1147 Double_t dy = dx*(c1+c2)/(r1+r2);
1150 Double_t delta = x[4]*dx*(c1+c2)/(c1*r2 + c2*r1);
1152 if (TMath::Abs(delta)>0.01){
1153 dz = x[3]*TMath::ASin(delta)/x[4];
1155 dz = x[3]*AliTPCFastMath::FastAsin(delta)/x[4];
1158 //dz = x[3]*AliTPCFastMath::FastAsin(delta)/x[4];
1166 Int_t AliTPCtrackerMI::LoadClusters (TTree *const tree)
1171 return LoadClusters();
1175 Int_t AliTPCtrackerMI::LoadClusters(const TObjArray *arr)
1178 // load clusters to the memory
1179 AliTPCClustersRow *clrow = new AliTPCClustersRow("AliTPCclusterMI");
1180 Int_t lower = arr->LowerBound();
1181 Int_t entries = arr->GetEntriesFast();
1183 for (Int_t i=lower; i<entries; i++) {
1184 clrow = (AliTPCClustersRow*) arr->At(i);
1185 if(!clrow) continue;
1186 if(!clrow->GetArray()) continue;
1190 fkParam->AdjustSectorRow(clrow->GetID(),sec,row);
1192 for (Int_t icl=0; icl<clrow->GetArray()->GetEntriesFast(); icl++){
1193 Transform((AliTPCclusterMI*)(clrow->GetArray()->At(icl)));
1196 if (clrow->GetArray()->GetEntriesFast()<=0) continue;
1197 AliTPCtrackerRow * tpcrow=0;
1200 tpcrow = &(fInnerSec[sec%fkNIS][row]);
1204 tpcrow = &(fOuterSec[(sec-fkNIS*2)%fkNOS][row]);
1205 left = (sec-fkNIS*2)/fkNOS;
1208 tpcrow->SetN1(clrow->GetArray()->GetEntriesFast());
1209 for (Int_t j=0;j<tpcrow->GetN1();++j)
1210 tpcrow->SetCluster1(j, *(AliTPCclusterMI*)(clrow->GetArray()->At(j)));
1213 tpcrow->SetN2(clrow->GetArray()->GetEntriesFast());
1214 for (Int_t j=0;j<tpcrow->GetN2();++j)
1215 tpcrow->SetCluster2(j, *(AliTPCclusterMI*)(clrow->GetArray()->At(j)));
1217 clrow->GetArray()->Clear("C");
1226 Int_t AliTPCtrackerMI::LoadClusters(const TClonesArray *arr)
1229 // load clusters to the memory from one
1232 AliTPCclusterMI *clust=0;
1233 Int_t count[72][96] = { {0} , {0} };
1235 // loop over clusters
1236 for (Int_t icl=0; icl<arr->GetEntriesFast(); icl++) {
1237 clust = (AliTPCclusterMI*)arr->At(icl);
1238 if(!clust) continue;
1239 //printf("cluster: det %d, row %d \n", clust->GetDetector(),clust->GetRow());
1241 // transform clusters
1244 // count clusters per pad row
1245 count[clust->GetDetector()][clust->GetRow()]++;
1248 // insert clusters to sectors
1249 for (Int_t icl=0; icl<arr->GetEntriesFast(); icl++) {
1250 clust = (AliTPCclusterMI*)arr->At(icl);
1251 if(!clust) continue;
1253 Int_t sec = clust->GetDetector();
1254 Int_t row = clust->GetRow();
1256 // filter overlapping pad rows needed by HLT
1257 if(sec<fkNIS*2) { //IROCs
1258 if(row == 30) continue;
1261 if(row == 27 || row == 76) continue;
1267 fInnerSec[sec%fkNIS].InsertCluster(clust, count[sec][row], fkParam);
1270 left = (sec-fkNIS*2)/fkNOS;
1271 fOuterSec[(sec-fkNIS*2)%fkNOS].InsertCluster(clust, count[sec][row], fkParam);
1275 // Load functions must be called behind LoadCluster(TClonesArray*)
1277 //LoadOuterSectors();
1278 //LoadInnerSectors();
1284 Int_t AliTPCtrackerMI::LoadClusters()
1287 // load clusters to the memory
1288 static AliTPCClustersRow *clrow= new AliTPCClustersRow("AliTPCclusterMI");
1290 // TTree * tree = fClustersArray.GetTree();
1292 TTree * tree = fInput;
1293 TBranch * br = tree->GetBranch("Segment");
1294 br->SetAddress(&clrow);
1296 // Conversion of pad, row coordinates in local tracking coords.
1297 // Could be skipped here; is already done in clusterfinder
1299 Int_t j=Int_t(tree->GetEntries());
1300 for (Int_t i=0; i<j; i++) {
1304 fkParam->AdjustSectorRow(clrow->GetID(),sec,row);
1305 for (Int_t icl=0; icl<clrow->GetArray()->GetEntriesFast(); icl++){
1306 Transform((AliTPCclusterMI*)(clrow->GetArray()->At(icl)));
1309 AliTPCtrackerRow * tpcrow=0;
1312 tpcrow = &(fInnerSec[sec%fkNIS][row]);
1316 tpcrow = &(fOuterSec[(sec-fkNIS*2)%fkNOS][row]);
1317 left = (sec-fkNIS*2)/fkNOS;
1320 tpcrow->SetN1(clrow->GetArray()->GetEntriesFast());
1321 for (Int_t k=0;k<tpcrow->GetN1();++k)
1322 tpcrow->SetCluster1(k, *(AliTPCclusterMI*)(clrow->GetArray()->At(k)));
1325 tpcrow->SetN2(clrow->GetArray()->GetEntriesFast());
1326 for (Int_t k=0;k<tpcrow->GetN2();++k)
1327 tpcrow->SetCluster2(k, *(AliTPCclusterMI*)(clrow->GetArray()->At(k)));
1338 void AliTPCtrackerMI::UnloadClusters()
1341 // unload clusters from the memory
1343 Int_t nrows = fOuterSec->GetNRows();
1344 for (Int_t sec = 0;sec<fkNOS;sec++)
1345 for (Int_t row = 0;row<nrows;row++){
1346 AliTPCtrackerRow* tpcrow = &(fOuterSec[sec%fkNOS][row]);
1348 // if (tpcrow->fClusters1) delete []tpcrow->fClusters1;
1349 // if (tpcrow->fClusters2) delete []tpcrow->fClusters2;
1351 tpcrow->ResetClusters();
1354 nrows = fInnerSec->GetNRows();
1355 for (Int_t sec = 0;sec<fkNIS;sec++)
1356 for (Int_t row = 0;row<nrows;row++){
1357 AliTPCtrackerRow* tpcrow = &(fInnerSec[sec%fkNIS][row]);
1359 // if (tpcrow->fClusters1) delete []tpcrow->fClusters1;
1360 //if (tpcrow->fClusters2) delete []tpcrow->fClusters2;
1362 tpcrow->ResetClusters();
1368 void AliTPCtrackerMI::FillClusterArray(TObjArray* array) const{
1370 // Filling cluster to the array - For visualization purposes
1373 nrows = fOuterSec->GetNRows();
1374 for (Int_t sec = 0;sec<fkNOS;sec++)
1375 for (Int_t row = 0;row<nrows;row++){
1376 AliTPCtrackerRow* tpcrow = &(fOuterSec[sec%fkNOS][row]);
1377 if (!tpcrow) continue;
1378 for (Int_t icl = 0;icl<tpcrow->GetN();icl++){
1379 array->AddLast((TObject*)((*tpcrow)[icl]));
1382 nrows = fInnerSec->GetNRows();
1383 for (Int_t sec = 0;sec<fkNIS;sec++)
1384 for (Int_t row = 0;row<nrows;row++){
1385 AliTPCtrackerRow* tpcrow = &(fInnerSec[sec%fkNIS][row]);
1386 if (!tpcrow) continue;
1387 for (Int_t icl = 0;icl<tpcrow->GetN();icl++){
1388 array->AddLast((TObject*)(*tpcrow)[icl]);
1394 void AliTPCtrackerMI::Transform(AliTPCclusterMI * cluster){
1398 AliTPCcalibDB * calibDB = AliTPCcalibDB::Instance();
1399 AliTPCTransform *transform = calibDB->GetTransform() ;
1401 AliFatal("Tranformations not in calibDB");
1404 transform->SetCurrentRecoParam((AliTPCRecoParam*)AliTPCReconstructor::GetRecoParam());
1405 Double_t x[3]={cluster->GetRow(),cluster->GetPad(),cluster->GetTimeBin()};
1406 Int_t i[1]={cluster->GetDetector()};
1407 transform->Transform(x,i,0,1);
1408 // if (cluster->GetDetector()%36>17){
1413 // in debug mode check the transformation
1415 if (AliTPCReconstructor::StreamLevel()>2) {
1417 cluster->GetGlobalXYZ(gx);
1418 Int_t event = (fEvent==NULL)? 0: fEvent->GetEventNumberInFile();
1419 TTreeSRedirector &cstream = *fDebugStreamer;
1420 cstream<<"Transform"<<
1431 cluster->SetX(x[0]);
1432 cluster->SetY(x[1]);
1433 cluster->SetZ(x[2]);
1438 //if (!fkParam->IsGeoRead()) fkParam->ReadGeoMatrices();
1439 if (AliTPCReconstructor::GetRecoParam()->GetUseSectorAlignment() && (!calibDB->HasAlignmentOCDB())){
1440 TGeoHMatrix *mat = fkParam->GetClusterMatrix(cluster->GetDetector());
1442 Double_t pos[3]= {cluster->GetX(),cluster->GetY(),cluster->GetZ()};
1443 Double_t posC[3]={cluster->GetX(),cluster->GetY(),cluster->GetZ()};
1444 if (mat) mat->LocalToMaster(pos,posC);
1446 // chack Loading of Geo matrices from GeoManager - TEMPORARY FIX
1448 cluster->SetX(posC[0]);
1449 cluster->SetY(posC[1]);
1450 cluster->SetZ(posC[2]);
1454 //_____________________________________________________________________________
1455 Int_t AliTPCtrackerMI::LoadOuterSectors() {
1456 //-----------------------------------------------------------------
1457 // This function fills outer TPC sectors with clusters.
1458 //-----------------------------------------------------------------
1459 Int_t nrows = fOuterSec->GetNRows();
1461 for (Int_t sec = 0;sec<fkNOS;sec++)
1462 for (Int_t row = 0;row<nrows;row++){
1463 AliTPCtrackerRow* tpcrow = &(fOuterSec[sec%fkNOS][row]);
1464 Int_t sec2 = sec+2*fkNIS;
1466 Int_t ncl = tpcrow->GetN1();
1468 AliTPCclusterMI *c= (tpcrow->GetCluster1(ncl));
1469 index=(((sec2<<8)+row)<<16)+ncl;
1470 tpcrow->InsertCluster(c,index);
1473 ncl = tpcrow->GetN2();
1475 AliTPCclusterMI *c= (tpcrow->GetCluster2(ncl));
1476 index=((((sec2+fkNOS)<<8)+row)<<16)+ncl;
1477 tpcrow->InsertCluster(c,index);
1480 // write indexes for fast acces
1482 for (Int_t i=0;i<510;i++)
1483 tpcrow->SetFastCluster(i,-1);
1484 for (Int_t i=0;i<tpcrow->GetN();i++){
1485 Int_t zi = Int_t((*tpcrow)[i]->GetZ()+255.);
1486 tpcrow->SetFastCluster(zi,i); // write index
1489 for (Int_t i=0;i<510;i++){
1490 if (tpcrow->GetFastCluster(i)<0)
1491 tpcrow->SetFastCluster(i,last);
1493 last = tpcrow->GetFastCluster(i);
1502 //_____________________________________________________________________________
1503 Int_t AliTPCtrackerMI::LoadInnerSectors() {
1504 //-----------------------------------------------------------------
1505 // This function fills inner TPC sectors with clusters.
1506 //-----------------------------------------------------------------
1507 Int_t nrows = fInnerSec->GetNRows();
1509 for (Int_t sec = 0;sec<fkNIS;sec++)
1510 for (Int_t row = 0;row<nrows;row++){
1511 AliTPCtrackerRow* tpcrow = &(fInnerSec[sec%fkNIS][row]);
1514 Int_t ncl = tpcrow->GetN1();
1516 AliTPCclusterMI *c= (tpcrow->GetCluster1(ncl));
1517 index=(((sec<<8)+row)<<16)+ncl;
1518 tpcrow->InsertCluster(c,index);
1521 ncl = tpcrow->GetN2();
1523 AliTPCclusterMI *c= (tpcrow->GetCluster2(ncl));
1524 index=((((sec+fkNIS)<<8)+row)<<16)+ncl;
1525 tpcrow->InsertCluster(c,index);
1528 // write indexes for fast acces
1530 for (Int_t i=0;i<510;i++)
1531 tpcrow->SetFastCluster(i,-1);
1532 for (Int_t i=0;i<tpcrow->GetN();i++){
1533 Int_t zi = Int_t((*tpcrow)[i]->GetZ()+255.);
1534 tpcrow->SetFastCluster(zi,i); // write index
1537 for (Int_t i=0;i<510;i++){
1538 if (tpcrow->GetFastCluster(i)<0)
1539 tpcrow->SetFastCluster(i,last);
1541 last = tpcrow->GetFastCluster(i);
1553 //_________________________________________________________________________
1554 AliTPCclusterMI *AliTPCtrackerMI::GetClusterMI(Int_t index) const {
1555 //--------------------------------------------------------------------
1556 // Return pointer to a given cluster
1557 //--------------------------------------------------------------------
1558 if (index<0) return 0; // no cluster
1559 Int_t sec=(index&0xff000000)>>24;
1560 Int_t row=(index&0x00ff0000)>>16;
1561 Int_t ncl=(index&0x00007fff)>>00;
1563 const AliTPCtrackerRow * tpcrow=0;
1564 TClonesArray * clrow =0;
1566 if (sec<0 || sec>=fkNIS*4) {
1567 AliWarning(Form("Wrong sector %d",sec));
1572 AliTPCtrackerSector& tracksec = fInnerSec[sec%fkNIS];
1573 if (tracksec.GetNRows()<=row) return 0;
1574 tpcrow = &(tracksec[row]);
1575 if (tpcrow==0) return 0;
1578 if (tpcrow->GetN1()<=ncl) return 0;
1579 clrow = tpcrow->GetClusters1();
1582 if (tpcrow->GetN2()<=ncl) return 0;
1583 clrow = tpcrow->GetClusters2();
1587 AliTPCtrackerSector& tracksec = fOuterSec[(sec-fkNIS*2)%fkNOS];
1588 if (tracksec.GetNRows()<=row) return 0;
1589 tpcrow = &(tracksec[row]);
1590 if (tpcrow==0) return 0;
1592 if (sec-2*fkNIS<fkNOS) {
1593 if (tpcrow->GetN1()<=ncl) return 0;
1594 clrow = tpcrow->GetClusters1();
1597 if (tpcrow->GetN2()<=ncl) return 0;
1598 clrow = tpcrow->GetClusters2();
1602 return (AliTPCclusterMI*)clrow->At(ncl);
1608 Int_t AliTPCtrackerMI::FollowToNext(AliTPCseed& t, Int_t nr) {
1609 //-----------------------------------------------------------------
1610 // This function tries to find a track prolongation to next pad row
1611 //-----------------------------------------------------------------
1613 Double_t x= GetXrow(nr), ymax=GetMaxY(nr);
1616 AliTPCclusterMI *cl=0;
1617 Int_t tpcindex= t.GetClusterIndex2(nr);
1619 // update current shape info every 5 pad-row
1620 // if ( (nr%5==0) || t.GetNumberOfClusters()<2 || (t.fCurrentSigmaY2<0.0001) ){
1624 if (fIteration>0 && tpcindex>=-1){ //if we have already clusters
1626 if (tpcindex==-1) return 0; //track in dead zone
1627 if (tpcindex >= 0){ //
1628 cl = t.GetClusterPointer(nr);
1629 //if (cl==0) cl = GetClusterMI(tpcindex);
1630 if (!cl) cl = GetClusterMI(tpcindex);
1631 t.SetCurrentClusterIndex1(tpcindex);
1634 Int_t relativesector = ((tpcindex&0xff000000)>>24)%18; // if previously accepted cluster in different sector
1635 Float_t angle = relativesector*fSectors->GetAlpha()+fSectors->GetAlphaShift();
1637 if (angle<-TMath::Pi()) angle += 2*TMath::Pi();
1638 if (angle>=TMath::Pi()) angle -= 2*TMath::Pi();
1640 if (TMath::Abs(angle-t.GetAlpha())>0.001){
1641 Double_t rotation = angle-t.GetAlpha();
1642 t.SetRelativeSector(relativesector);
1643 if (!t.Rotate(rotation)) return 0;
1645 if (!t.PropagateTo(x)) return 0;
1647 t.SetCurrentCluster(cl);
1649 Int_t accept = AcceptCluster(&t,t.GetCurrentCluster());
1650 if ((tpcindex&0x8000)==0) accept =0;
1652 //if founded cluster is acceptible
1653 if (cl->IsUsed(11)) { // id cluster is shared inrease uncertainty
1654 t.SetErrorY2(t.GetErrorY2()+0.03);
1655 t.SetErrorZ2(t.GetErrorZ2()+0.03);
1656 t.SetErrorY2(t.GetErrorY2()*3);
1657 t.SetErrorZ2(t.GetErrorZ2()*3);
1659 t.SetNFoundable(t.GetNFoundable()+1);
1660 UpdateTrack(&t,accept);
1663 else { // Remove old cluster from track
1664 t.SetClusterIndex(nr, -3);
1665 t.SetClusterPointer(nr, 0);
1669 if (TMath::Abs(t.GetSnp())>AliTPCReconstructor::GetMaxSnpTracker()) return 0; // cut on angle
1670 if (fIteration>1 && IsFindable(t)){
1671 // not look for new cluster during refitting
1672 t.SetNFoundable(t.GetNFoundable()+1);
1677 // if (TMath::Abs(t.GetSnp())>0.95 || TMath::Abs(x*t.GetC()-t.GetEta())>0.95) return 0;// patch 28 fev 06
1678 if (!t.PropagateTo(x)) {
1679 if (fIteration==0) t.SetRemoval(10);
1682 Double_t y = t.GetY();
1683 if (TMath::Abs(y)>ymax){
1685 t.SetRelativeSector((t.GetRelativeSector()+1) % fN);
1686 if (!t.Rotate(fSectors->GetAlpha()))
1688 } else if (y <-ymax) {
1689 t.SetRelativeSector((t.GetRelativeSector()-1+fN) % fN);
1690 if (!t.Rotate(-fSectors->GetAlpha()))
1693 if (!t.PropagateTo(x)) {
1694 if (fIteration==0) t.SetRemoval(10);
1700 Double_t z=t.GetZ();
1703 if (!IsActive(t.GetRelativeSector(),nr)) {
1705 t.SetClusterIndex2(nr,-1);
1708 //AliInfo(Form("A - Sector%d phi %f - alpha %f", t.fRelativeSector,y/x, t.GetAlpha()));
1709 Bool_t isActive = IsActive(t.GetRelativeSector(),nr);
1710 Bool_t isActive2 = (nr>=fInnerSec->GetNRows()) ? fOuterSec[t.GetRelativeSector()][nr-fInnerSec->GetNRows()].GetN()>0:fInnerSec[t.GetRelativeSector()][nr].GetN()>0;
1712 if (!isActive || !isActive2) return 0;
1714 const AliTPCtrackerRow &krow=GetRow(t.GetRelativeSector(),nr);
1715 if ( (t.GetSigmaY2()<0) || t.GetSigmaZ2()<0) return 0;
1717 Double_t roadz = 1.;
1719 if (TMath::Abs(TMath::Abs(y)-ymax)<krow.GetDeadZone()){
1721 t.SetClusterIndex2(nr,-1);
1727 // if (TMath::Abs(z)<(AliTPCReconstructor::GetCtgRange()*x+10) && TMath::Abs(z)<fkParam->GetZLength(0) && (TMath::Abs(t.GetSnp())<AliTPCReconstructor::GetMaxSnpTracker()))
1728 t.SetNFoundable(t.GetNFoundable()+1);
1734 // cl = krow.FindNearest2(y+10.,z,roady,roadz,index);
1735 cl = krow.FindNearest2(y,z,roady,roadz,index);
1736 if (cl) t.SetCurrentClusterIndex1(krow.GetIndex(index));
1739 t.SetCurrentCluster(cl);
1741 if (fIteration==2&&cl->IsUsed(10)) return 0;
1742 Int_t accept = AcceptCluster(&t,t.GetCurrentCluster());
1743 if (fIteration==2&&cl->IsUsed(11)) {
1744 t.SetErrorY2(t.GetErrorY2()+0.03);
1745 t.SetErrorZ2(t.GetErrorZ2()+0.03);
1746 t.SetErrorY2(t.GetErrorY2()*3);
1747 t.SetErrorZ2(t.GetErrorZ2()*3);
1750 if (t.fCurrentCluster->IsUsed(10)){
1755 if (t.fNShared>0.7*t.GetNumberOfClusters()) {
1761 if (accept<3) UpdateTrack(&t,accept);
1764 if ( fIteration==0 && t.GetNFoundable()*0.5 > t.GetNumberOfClusters()) t.SetRemoval(10);
1772 //_________________________________________________________________________
1773 Bool_t AliTPCtrackerMI::GetTrackPoint(Int_t index, AliTrackPoint &p ) const
1775 // Get track space point by index
1776 // return false in case the cluster doesn't exist
1777 AliTPCclusterMI *cl = GetClusterMI(index);
1778 if (!cl) return kFALSE;
1779 Int_t sector = (index&0xff000000)>>24;
1780 // Int_t row = (index&0x00ff0000)>>16;
1782 // xyz[0] = fkParam->GetPadRowRadii(sector,row);
1783 xyz[0] = cl->GetX();
1784 xyz[1] = cl->GetY();
1785 xyz[2] = cl->GetZ();
1787 fkParam->AdjustCosSin(sector,cos,sin);
1788 Float_t x = cos*xyz[0]-sin*xyz[1];
1789 Float_t y = cos*xyz[1]+sin*xyz[0];
1791 Float_t sigmaY2 = 0.027*cl->GetSigmaY2();
1792 if (sector < fkParam->GetNInnerSector()) sigmaY2 *= 2.07;
1793 Float_t sigmaZ2 = 0.066*cl->GetSigmaZ2();
1794 if (sector < fkParam->GetNInnerSector()) sigmaZ2 *= 1.77;
1795 cov[0] = sin*sin*sigmaY2;
1796 cov[1] = -sin*cos*sigmaY2;
1798 cov[3] = cos*cos*sigmaY2;
1801 p.SetXYZ(x,y,xyz[2],cov);
1802 AliGeomManager::ELayerID iLayer;
1804 if (sector < fkParam->GetNInnerSector()) {
1805 iLayer = AliGeomManager::kTPC1;
1809 iLayer = AliGeomManager::kTPC2;
1810 idet = sector - fkParam->GetNInnerSector();
1812 UShort_t volid = AliGeomManager::LayerToVolUID(iLayer,idet);
1813 p.SetVolumeID(volid);
1819 Int_t AliTPCtrackerMI::UpdateClusters(AliTPCseed& t, Int_t nr) {
1820 //-----------------------------------------------------------------
1821 // This function tries to find a track prolongation to next pad row
1822 //-----------------------------------------------------------------
1823 t.SetCurrentCluster(0);
1824 t.SetCurrentClusterIndex1(-3);
1826 Double_t xt=t.GetX();
1827 Int_t row = GetRowNumber(xt)-1;
1828 Double_t ymax= GetMaxY(nr);
1830 if (row < nr) return 1; // don't prolongate if not information until now -
1831 // if (TMath::Abs(t.GetSnp())>0.9 && t.GetNumberOfClusters()>40. && fIteration!=2) {
1833 // return 0; // not prolongate strongly inclined tracks
1835 // if (TMath::Abs(t.GetSnp())>0.95) {
1837 // return 0; // not prolongate strongly inclined tracks
1838 // }// patch 28 fev 06
1840 Double_t x= GetXrow(nr);
1842 //t.PropagateTo(x+0.02);
1843 //t.PropagateTo(x+0.01);
1844 if (!t.PropagateTo(x)){
1851 if (TMath::Abs(y)>ymax){
1853 t.SetRelativeSector((t.GetRelativeSector()+1) % fN);
1854 if (!t.Rotate(fSectors->GetAlpha()))
1856 } else if (y <-ymax) {
1857 t.SetRelativeSector((t.GetRelativeSector()-1+fN) % fN);
1858 if (!t.Rotate(-fSectors->GetAlpha()))
1861 // if (!t.PropagateTo(x)){
1868 if (TMath::Abs(t.GetSnp())>AliTPCReconstructor::GetMaxSnpTracker()) return 0;
1870 if (!IsActive(t.GetRelativeSector(),nr)) {
1872 t.SetClusterIndex2(nr,-1);
1875 //AliInfo(Form("A - Sector%d phi %f - alpha %f", t.fRelativeSector,y/x, t.GetAlpha()));
1877 AliTPCtrackerRow &krow=GetRow(t.GetRelativeSector(),nr);
1879 if (TMath::Abs(TMath::Abs(y)-ymax)<krow.GetDeadZone()){
1881 t.SetClusterIndex2(nr,-1);
1887 // if (TMath::Abs(t.GetZ())<(AliTPCReconstructor::GetCtgRange()*t.GetX()+10) && (TMath::Abs(t.GetSnp())<AliTPCReconstructor::GetMaxSnpTracker()))
1888 if (IsFindable(t)) t.SetNFoundable(t.GetNFoundable()+1);
1894 if ( (nr%2==0) || t.GetNumberOfClusters()<2){
1895 // t.fCurrentSigmaY = GetSigmaY(&t);
1896 //t.fCurrentSigmaZ = GetSigmaZ(&t);
1900 AliTPCclusterMI *cl=0;
1903 Double_t roady = 1.;
1904 Double_t roadz = 1.;
1908 index = t.GetClusterIndex2(nr);
1909 if ( (index >= 0) && (index&0x8000)==0){
1910 cl = t.GetClusterPointer(nr);
1911 if ( (cl==0) && (index >= 0)) cl = GetClusterMI(index);
1912 t.SetCurrentClusterIndex1(index);
1914 t.SetCurrentCluster(cl);
1920 // if (index<0) return 0;
1921 UInt_t uindex = TMath::Abs(index);
1924 //cl = krow.FindNearest2(y+10,z,roady,roadz,uindex);
1925 cl = krow.FindNearest2(y,z,roady,roadz,uindex);
1928 if (cl) t.SetCurrentClusterIndex1(krow.GetIndex(uindex));
1929 t.SetCurrentCluster(cl);
1935 Int_t AliTPCtrackerMI::FollowToNextCluster(AliTPCseed & t, Int_t nr) {
1936 //-----------------------------------------------------------------
1937 // This function tries to find a track prolongation to next pad row
1938 //-----------------------------------------------------------------
1940 //update error according neighborhoud
1942 if (t.GetCurrentCluster()) {
1944 Int_t accept = AcceptCluster(&t,t.GetCurrentCluster());
1946 if (t.GetCurrentCluster()->IsUsed(10)){
1951 t.SetNShared(t.GetNShared()+1);
1952 if (t.GetNShared()>0.7*t.GetNumberOfClusters()) {
1957 if (fIteration>0) accept = 0;
1958 if (accept<3) UpdateTrack(&t,accept);
1962 if ( t.GetNumberOfClusters()>18 && ( (t.GetSigmaY2()+t.GetSigmaZ2())>0.16)) t.SetRemoval(10);
1963 if ( t.GetNumberOfClusters()>18 && t.GetChi2()/t.GetNumberOfClusters()>6 ) t.SetRemoval(10);
1965 if (( (t.GetNFoundable()*0.5 > t.GetNumberOfClusters()) || t.GetNoCluster()>15)) t.SetRemoval(10);
1973 //_____________________________________________________________________________
1974 Int_t AliTPCtrackerMI::FollowProlongation(AliTPCseed& t, Int_t rf, Int_t step, Bool_t fromSeeds) {
1975 //-----------------------------------------------------------------
1976 // This function tries to find a track prolongation.
1977 //-----------------------------------------------------------------
1978 Double_t xt=t.GetX();
1980 Double_t alpha=t.GetAlpha();
1981 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
1982 if (alpha < 0. ) alpha += 2.*TMath::Pi();
1984 t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN);
1986 Int_t first = GetRowNumber(xt);
1991 for (Int_t nr= first; nr>=rf; nr-=step) {
1993 if (t.GetKinkIndexes()[0]>0){
1994 for (Int_t i=0;i<3;i++){
1995 Int_t index = t.GetKinkIndexes()[i];
1996 if (index==0) break;
1997 if (index<0) continue;
1999 AliKink * kink = (AliKink*)fEvent->GetKink(index-1);
2001 printf("PROBLEM\n");
2004 Int_t kinkrow = kink->GetTPCRow0()+2+Int_t(0.5/(0.05+kink->GetAngle(2)));
2006 AliExternalTrackParam paramd(t);
2007 kink->SetDaughter(paramd);
2008 kink->SetStatus(2,5);
2015 if (nr==80) t.UpdateReference();
2016 if (nr<fInnerSec->GetNRows())
2017 fSectors = fInnerSec;
2019 fSectors = fOuterSec;
2020 if (FollowToNext(t,nr)==0)
2033 Int_t AliTPCtrackerMI::FollowBackProlongation(AliTPCseed& t, Int_t rf, Bool_t fromSeeds) {
2034 //-----------------------------------------------------------------
2035 // This function tries to find a track prolongation.
2036 //-----------------------------------------------------------------
2038 Double_t xt=t.GetX();
2039 Double_t alpha=t.GetAlpha();
2040 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
2041 if (alpha < 0. ) alpha += 2.*TMath::Pi();
2042 t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN);
2044 Int_t first = t.GetFirstPoint();
2045 Int_t ri = GetRowNumber(xt);
2049 if (first<ri) first = ri;
2051 if (first<0) first=0;
2052 for (Int_t nr=first; nr<=rf; nr++) {
2053 // if ( (TMath::Abs(t.GetSnp())>0.95)) break;//patch 28 fev 06
2054 if (t.GetKinkIndexes()[0]<0){
2055 for (Int_t i=0;i<3;i++){
2056 Int_t index = t.GetKinkIndexes()[i];
2057 if (index==0) break;
2058 if (index>0) continue;
2059 index = TMath::Abs(index);
2060 AliKink * kink = (AliKink*)fEvent->GetKink(index-1);
2062 printf("PROBLEM\n");
2065 Int_t kinkrow = kink->GetTPCRow0()-2-Int_t(0.5/(0.05+kink->GetAngle(2)));
2067 AliExternalTrackParam paramm(t);
2068 kink->SetMother(paramm);
2069 kink->SetStatus(2,1);
2076 if (nr<fInnerSec->GetNRows())
2077 fSectors = fInnerSec;
2079 fSectors = fOuterSec;
2090 Float_t AliTPCtrackerMI::OverlapFactor(AliTPCseed * s1, AliTPCseed * s2, Int_t &sum1, Int_t & sum2)
2092 // overlapping factor
2098 Float_t dz2 =(s1->GetZ() - s2->GetZ());
2101 Float_t dy2 =TMath::Abs((s1->GetY() - s2->GetY()));
2103 Float_t distance = TMath::Sqrt(dz2+dy2);
2104 if (distance>4.) return 0; // if there are far away - not overlap - to reduce combinatorics
2107 Int_t firstpoint = TMath::Min(s1->GetFirstPoint(),s2->GetFirstPoint());
2108 Int_t lastpoint = TMath::Max(s1->GetLastPoint(),s2->GetLastPoint());
2113 if (firstpoint>lastpoint) {
2114 firstpoint =lastpoint;
2119 for (Int_t i=firstpoint-1;i<lastpoint+1;i++){
2120 if (s1->GetClusterIndex2(i)>0) sum1++;
2121 if (s2->GetClusterIndex2(i)>0) sum2++;
2122 if (s1->GetClusterIndex2(i)==s2->GetClusterIndex2(i) && s1->GetClusterIndex2(i)>0) {
2126 if (sum<5) return 0;
2128 Float_t summin = TMath::Min(sum1+1,sum2+1);
2129 Float_t ratio = (sum+1)/Float_t(summin);
2133 void AliTPCtrackerMI::SignShared(AliTPCseed * s1, AliTPCseed * s2)
2137 Float_t thetaCut = 0.2;//+10.*TMath::Sqrt(s1->GetSigmaTglZ()+ s2->GetSigmaTglZ());
2138 if (TMath::Abs(s1->GetTgl()-s2->GetTgl())>thetaCut) return;
2139 Float_t minCl = TMath::Min(s1->GetNumberOfClusters(),s2->GetNumberOfClusters());
2140 Int_t cutN0 = TMath::Max(5,TMath::Nint(0.1*minCl));
2145 //Int_t firstpoint = TMath::Max(s1->GetFirstPoint(),s2->GetFirstPoint());
2146 //Int_t lastpoint = TMath::Min(s1->GetLastPoint(),s2->GetLastPoint());
2147 Int_t firstpoint = 0;
2148 Int_t lastpoint = 160;
2150 // if (firstpoint>=lastpoint-5) return;;
2152 for (Int_t i=firstpoint;i<lastpoint;i++){
2153 // if ( (s1->GetClusterIndex2(i)&0xFFFF8FFF)==(s2->GetClusterIndex2(i)&0xFFFF8FFF) && s1->GetClusterIndex2(i)>0) {
2154 if ( (s1->GetClusterIndex2(i))==(s2->GetClusterIndex2(i)) && s1->GetClusterIndex2(i)>0) {
2158 if (sumshared>cutN0){
2161 for (Int_t i=firstpoint;i<lastpoint;i++){
2162 // if ( (s1->GetClusterIndex2(i)&0xFFFF8FFF)==(s2->GetClusterIndex2(i)&0xFFFF8FFF) && s1->GetClusterIndex2(i)>0) {
2163 if ( (s1->GetClusterIndex2(i))==(s2->GetClusterIndex2(i)) && s1->GetClusterIndex2(i)>0) {
2164 AliTPCTrackerPoint *p1 = s1->GetTrackPoint(i);
2165 AliTPCTrackerPoint *p2 = s2->GetTrackPoint(i);;
2166 if (s1->IsActive()&&s2->IsActive()){
2167 p1->SetShared(kTRUE);
2168 p2->SetShared(kTRUE);
2174 if (sumshared>cutN0){
2175 for (Int_t i=0;i<4;i++){
2176 if (s1->GetOverlapLabel(3*i)==0){
2177 s1->SetOverlapLabel(3*i, s2->GetLabel());
2178 s1->SetOverlapLabel(3*i+1,sumshared);
2179 s1->SetOverlapLabel(3*i+2,s2->GetUniqueID());
2183 for (Int_t i=0;i<4;i++){
2184 if (s2->GetOverlapLabel(3*i)==0){
2185 s2->SetOverlapLabel(3*i, s1->GetLabel());
2186 s2->SetOverlapLabel(3*i+1,sumshared);
2187 s2->SetOverlapLabel(3*i+2,s1->GetUniqueID());
2194 void AliTPCtrackerMI::SignShared(TObjArray * arr)
2197 //sort trackss according sectors
2199 for (Int_t i=0; i<arr->GetEntriesFast(); i++) {
2200 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2202 //if (pt) RotateToLocal(pt);
2206 arr->Sort(); // sorting according relative sectors
2207 arr->Expand(arr->GetEntries());
2210 Int_t nseed=arr->GetEntriesFast();
2211 for (Int_t i=0; i<nseed; i++) {
2212 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2214 for (Int_t j=0;j<12;j++){
2215 pt->SetOverlapLabel(j,0);
2218 for (Int_t i=0; i<nseed; i++) {
2219 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2221 if (pt->GetRemoval()>10) continue;
2222 for (Int_t j=i+1; j<nseed; j++){
2223 AliTPCseed *pt2=(AliTPCseed*)arr->UncheckedAt(j);
2224 if (TMath::Abs(pt->GetRelativeSector()-pt2->GetRelativeSector())>1) continue;
2226 if (pt2->GetRemoval()<=10) {
2227 //if ( TMath::Abs(pt->GetRelativeSector()-pt2->GetRelativeSector())>0) break;
2235 void AliTPCtrackerMI::SortTracks(TObjArray * arr, Int_t mode) const
2238 //sort tracks in array according mode criteria
2239 Int_t nseed = arr->GetEntriesFast();
2240 for (Int_t i=0; i<nseed; i++) {
2241 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2252 void AliTPCtrackerMI::RemoveUsed2(TObjArray * arr, Float_t factor1, Float_t factor2, Int_t minimal)
2255 // Loop over all tracks and remove overlaped tracks (with lower quality)
2257 // 1. Unsign clusters
2258 // 2. Sort tracks according quality
2259 // Quality is defined by the number of cluster between first and last points
2261 // 3. Loop over tracks - decreasing quality order
2262 // a.) remove - If the fraction of shared cluster less than factor (1- n or 2)
2263 // b.) remove - If the minimal number of clusters less than minimal and not ITS
2264 // c.) if track accepted - sign clusters
2266 //Called in - AliTPCtrackerMI::Clusters2Tracks()
2267 // - AliTPCtrackerMI::PropagateBack()
2268 // - AliTPCtrackerMI::RefitInward()
2271 // factor1 - factor for constrained
2272 // factor2 - for non constrained tracks
2273 // if (Float_t(shared+1)/Float_t(found+1)>factor) - DELETE
2277 Int_t nseed = arr->GetEntriesFast();
2278 Float_t * quality = new Float_t[nseed];
2279 Int_t * indexes = new Int_t[nseed];
2283 for (Int_t i=0; i<nseed; i++) {
2284 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2289 pt->UpdatePoints(); //select first last max dens points
2290 Float_t * points = pt->GetPoints();
2291 if (points[3]<0.8) quality[i] =-1;
2292 quality[i] = (points[2]-points[0])+pt->GetNumberOfClusters();
2293 //prefer high momenta tracks if overlaps
2294 quality[i] *= TMath::Sqrt(TMath::Abs(pt->Pt())+0.5);
2296 TMath::Sort(nseed,quality,indexes);
2299 for (Int_t itrack=0; itrack<nseed; itrack++) {
2300 Int_t trackindex = indexes[itrack];
2301 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(trackindex);
2304 if (quality[trackindex]<0){
2305 MarkSeedFree( arr->RemoveAt(trackindex) );
2310 Int_t first = Int_t(pt->GetPoints()[0]);
2311 Int_t last = Int_t(pt->GetPoints()[2]);
2312 Double_t factor = (pt->GetBConstrain()) ? factor1: factor2;
2314 Int_t found,foundable,shared;
2315 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
2316 // pt->GetClusterStatistic(0,160, found, foundable,shared,kFALSE);
2317 Bool_t itsgold =kFALSE;
2320 if (pt->GetESD()->GetITSclusters(dummy)>4) itsgold= kTRUE;
2324 if (Float_t(shared+1)/Float_t(found+1)>factor){
2325 if (pt->GetKinkIndexes()[0]!=0) continue; //don't remove tracks - part of the kinks
2326 if( AliTPCReconstructor::StreamLevel()>3){
2327 TTreeSRedirector &cstream = *fDebugStreamer;
2328 cstream<<"RemoveUsed"<<
2329 "iter="<<fIteration<<
2333 MarkSeedFree( arr->RemoveAt(trackindex) );
2336 if (pt->GetNumberOfClusters()<50&&(found-0.5*shared)<minimal){ //remove short tracks
2337 if (pt->GetKinkIndexes()[0]!=0) continue; //don't remove tracks - part of the kinks
2338 if( AliTPCReconstructor::StreamLevel()>3){
2339 TTreeSRedirector &cstream = *fDebugStreamer;
2340 cstream<<"RemoveShort"<<
2341 "iter="<<fIteration<<
2345 MarkSeedFree( arr->RemoveAt(trackindex) );
2351 //if (sharedfactor>0.4) continue;
2352 if (pt->GetKinkIndexes()[0]>0) continue;
2353 //Remove tracks with undefined properties - seems
2354 if (pt->GetSigmaY2()<kAlmost0) continue; // ? what is the origin ?
2356 for (Int_t i=first; i<last; i++) {
2357 Int_t index=pt->GetClusterIndex2(i);
2358 // if (index<0 || index&0x8000 ) continue;
2359 if (index<0 || index&0x8000 ) continue;
2360 AliTPCclusterMI *c= pt->GetClusterPointer(i);
2367 Info("RemoveUsed","\n*****\nNumber of good tracks after shared removal\t%d\n",fNtracks);
2373 void AliTPCtrackerMI::DumpClusters(Int_t iter, TObjArray *trackArray)
2376 // Dump clusters after reco
2377 // signed and unsigned cluster can be visualized
2378 // 1. Unsign all cluster
2379 // 2. Sign all used clusters
2382 Int_t nseed = trackArray->GetEntries();
2383 for (Int_t i=0; i<nseed; i++){
2384 AliTPCseed *pt=(AliTPCseed*)trackArray->UncheckedAt(i);
2388 Bool_t isKink=pt->GetKinkIndex(0)!=0;
2389 for (Int_t j=0; j<160; ++j) {
2390 Int_t index=pt->GetClusterIndex2(j);
2391 if (index<0) continue;
2392 AliTPCclusterMI *c= pt->GetClusterPointer(j);
2394 if (isKink) c->Use(100); // kink
2395 c->Use(10); // by default usage 10
2400 for (Int_t sec=0;sec<fkNIS;sec++){
2401 for (Int_t row=0;row<fInnerSec->GetNRows();row++){
2402 TClonesArray *cla = fInnerSec[sec][row].GetClusters1();
2403 for (Int_t icl =0;icl< fInnerSec[sec][row].GetN1();icl++){
2404 AliTPCclusterMI* cl = (AliTPCclusterMI*)cla->At(icl);
2405 Float_t gx[3]; cl->GetGlobalXYZ(gx);
2406 (*fDebugStreamer)<<"clDump"<<
2414 cla = fInnerSec[sec][row].GetClusters2();
2415 for (Int_t icl =0;icl< fInnerSec[sec][row].GetN2();icl++){
2416 AliTPCclusterMI* cl = (AliTPCclusterMI*)cla->At(icl);
2417 Float_t gx[3]; cl->GetGlobalXYZ(gx);
2418 (*fDebugStreamer)<<"clDump"<<
2429 for (Int_t sec=0;sec<fkNOS;sec++){
2430 for (Int_t row=0;row<fOuterSec->GetNRows();row++){
2431 TClonesArray *cla = fOuterSec[sec][row].GetClusters1();
2432 for (Int_t icl =0;icl< fOuterSec[sec][row].GetN1();icl++){
2434 AliTPCclusterMI* cl = (AliTPCclusterMI*) cla->At(icl);
2435 cl->GetGlobalXYZ(gx);
2436 (*fDebugStreamer)<<"clDump"<<
2444 cla = fOuterSec[sec][row].GetClusters2();
2445 for (Int_t icl =0;icl< fOuterSec[sec][row].GetN2();icl++){
2447 AliTPCclusterMI* cl = (AliTPCclusterMI*) cla->At(icl);
2448 cl->GetGlobalXYZ(gx);
2449 (*fDebugStreamer)<<"clDump"<<
2461 void AliTPCtrackerMI::UnsignClusters()
2464 // loop over all clusters and unsign them
2467 for (Int_t sec=0;sec<fkNIS;sec++){
2468 for (Int_t row=0;row<fInnerSec->GetNRows();row++){
2469 TClonesArray *cla = fInnerSec[sec][row].GetClusters1();
2470 for (Int_t icl =0;icl< fInnerSec[sec][row].GetN1();icl++)
2471 // if (cl[icl].IsUsed(10))
2472 ((AliTPCclusterMI*) cla->At(icl))->Use(-1);
2473 cla = fInnerSec[sec][row].GetClusters2();
2474 for (Int_t icl =0;icl< fInnerSec[sec][row].GetN2();icl++)
2475 //if (cl[icl].IsUsed(10))
2476 ((AliTPCclusterMI*) cla->At(icl))->Use(-1);
2480 for (Int_t sec=0;sec<fkNOS;sec++){
2481 for (Int_t row=0;row<fOuterSec->GetNRows();row++){
2482 TClonesArray *cla = fOuterSec[sec][row].GetClusters1();
2483 for (Int_t icl =0;icl< fOuterSec[sec][row].GetN1();icl++)
2484 //if (cl[icl].IsUsed(10))
2485 ((AliTPCclusterMI*) cla->At(icl))->Use(-1);
2486 cla = fOuterSec[sec][row].GetClusters2();
2487 for (Int_t icl =0;icl< fOuterSec[sec][row].GetN2();icl++)
2488 //if (cl[icl].IsUsed(10))
2489 ((AliTPCclusterMI*) cla->At(icl))->Use(-1);
2497 void AliTPCtrackerMI::SignClusters(const TObjArray * arr, Float_t fnumber, Float_t fdensity)
2500 //sign clusters to be "used"
2502 // snumber and sdensity sign number of sigmas - bellow mean value to be accepted
2503 // loop over "primaries"
2517 Int_t nseed = arr->GetEntriesFast();
2518 for (Int_t i=0; i<nseed; i++) {
2519 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2523 if (!(pt->IsActive())) continue;
2524 Float_t dens = pt->GetNumberOfClusters()/Float_t(pt->GetNFoundable());
2525 if ( (dens>0.7) && (pt->GetNumberOfClusters()>70)){
2527 sumdens2+= dens*dens;
2528 sumn += pt->GetNumberOfClusters();
2529 sumn2 += pt->GetNumberOfClusters()*pt->GetNumberOfClusters();
2530 Float_t chi2 = pt->GetChi2()/pt->GetNumberOfClusters();
2533 sumchi2 +=chi2*chi2;
2538 Float_t mdensity = 0.9;
2539 Float_t meann = 130;
2540 Float_t meanchi = 1;
2541 Float_t sdensity = 0.1;
2542 Float_t smeann = 10;
2543 Float_t smeanchi =0.4;
2547 mdensity = sumdens/sum;
2549 meanchi = sumchi/sum;
2551 sdensity = sumdens2/sum-mdensity*mdensity;
2553 sdensity = TMath::Sqrt(sdensity);
2557 smeann = sumn2/sum-meann*meann;
2559 smeann = TMath::Sqrt(smeann);
2563 smeanchi = sumchi2/sum - meanchi*meanchi;
2565 smeanchi = TMath::Sqrt(smeanchi);
2571 //REMOVE SHORT DELTAS or tracks going out of sensitive volume of TPC
2573 for (Int_t i=0; i<nseed; i++) {
2574 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2578 if (pt->GetBSigned()) continue;
2579 if (pt->GetBConstrain()) continue;
2580 //if (!(pt->IsActive())) continue;
2582 Int_t found,foundable,shared;
2583 pt->GetClusterStatistic(0,160,found, foundable,shared);
2584 if (shared/float(found)>0.3) {
2585 if (shared/float(found)>0.9 ){
2586 //MarkSeedFree( arr->RemoveAt(i) );
2591 Bool_t isok =kFALSE;
2592 if ( (pt->GetNShared()/pt->GetNumberOfClusters()<0.5) &&pt->GetNumberOfClusters()>60)
2594 if ((TMath::Abs(1/pt->GetC())<100.) && (pt->GetNShared()/pt->GetNumberOfClusters()<0.7))
2596 if (TMath::Abs(pt->GetZ()/pt->GetX())>1.1)
2598 if ( (TMath::Abs(pt->GetSnp()>0.7) && pt->GetD(0,0)>60.))
2602 for (Int_t j=0; j<160; ++j) {
2603 Int_t index=pt->GetClusterIndex2(j);
2604 if (index<0) continue;
2605 AliTPCclusterMI *c= pt->GetClusterPointer(j);
2607 //if (!(c->IsUsed(10))) c->Use();
2614 Double_t maxchi = meanchi+2.*smeanchi;
2616 for (Int_t i=0; i<nseed; i++) {
2617 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2621 //if (!(pt->IsActive())) continue;
2622 if (pt->GetBSigned()) continue;
2623 Double_t chi = pt->GetChi2()/pt->GetNumberOfClusters();
2624 if (chi>maxchi) continue;
2627 Float_t dens = pt->GetNumberOfClusters()/Float_t(pt->GetNFoundable());
2629 //sign only tracks with enoug big density at the beginning
2631 if ((pt->GetDensityFirst(40)<0.75) && pt->GetNumberOfClusters()<meann) continue;
2634 Double_t mindens = TMath::Max(double(mdensity-sdensity*fdensity*bfactor),0.65);
2635 Double_t minn = TMath::Max(Int_t(meann-fnumber*smeann*bfactor),50);
2637 // if (pt->fBConstrain) mindens = TMath::Max(mdensity-sdensity*fdensity*bfactor,0.65);
2638 if ( (pt->GetRemoval()==10) && (pt->GetSnp()>0.8)&&(dens>mindens))
2641 if ((dens>mindens && pt->GetNumberOfClusters()>minn) && chi<maxchi ){
2642 //Int_t noc=pt->GetNumberOfClusters();
2643 pt->SetBSigned(kTRUE);
2644 for (Int_t j=0; j<160; ++j) {
2646 Int_t index=pt->GetClusterIndex2(j);
2647 if (index<0) continue;
2648 AliTPCclusterMI *c= pt->GetClusterPointer(j);
2650 // if (!(c->IsUsed(10))) c->Use();
2655 // gLastCheck = nseed;
2664 Int_t AliTPCtrackerMI::RefitInward(AliESDEvent *event)
2667 // back propagation of ESD tracks
2670 if (!event) return 0;
2671 const Int_t kMaxFriendTracks=2000;
2673 // extract correction object for multiplicity dependence of dEdx
2674 TObjArray * gainCalibArray = AliTPCcalibDB::Instance()->GetTimeGainSplinesRun(event->GetRunNumber());
2676 AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ;
2678 AliFatal("Tranformations not in RefitInward");
2681 transform->SetCurrentRecoParam((AliTPCRecoParam*)AliTPCReconstructor::GetRecoParam());
2682 const AliTPCRecoParam * recoParam = AliTPCcalibDB::Instance()->GetTransform()->GetCurrentRecoParam();
2683 Int_t nContribut = event->GetNumberOfTracks();
2684 TGraphErrors * graphMultDependenceDeDx = 0x0;
2685 if (recoParam && recoParam->GetUseMultiplicityCorrectionDedx() && gainCalibArray) {
2686 if (recoParam->GetUseTotCharge()) {
2687 graphMultDependenceDeDx = (TGraphErrors *) gainCalibArray->FindObject("TGRAPHERRORS_MEANQTOT_MULTIPLICITYDEPENDENCE_BEAM_ALL");
2689 graphMultDependenceDeDx = (TGraphErrors *) gainCalibArray->FindObject("TGRAPHERRORS_MEANQMAX_MULTIPLICITYDEPENDENCE_BEAM_ALL");
2695 //PrepareForProlongation(fSeeds,1);
2696 PropagateForward2(fSeeds);
2697 RemoveUsed2(fSeeds,0.4,0.4,20);
2699 TObjArray arraySeed(fSeeds->GetEntries());
2700 for (Int_t i=0;i<fSeeds->GetEntries();i++) {
2701 arraySeed.AddAt(fSeeds->At(i),i);
2703 SignShared(&arraySeed);
2704 // FindCurling(fSeeds, event,2); // find multi found tracks
2705 FindSplitted(fSeeds, event,2); // find multi found tracks
2706 if (AliTPCReconstructor::StreamLevel()>5) FindMultiMC(fSeeds, fEvent,2); // find multi found tracks
2709 Int_t nseed = fSeeds->GetEntriesFast();
2710 for (Int_t i=0;i<nseed;i++){
2711 AliTPCseed * seed = (AliTPCseed*) fSeeds->UncheckedAt(i);
2712 if (!seed) continue;
2713 if (seed->GetKinkIndex(0)>0) UpdateKinkQualityD(seed); // update quality informations for kinks
2714 AliESDtrack *esd=event->GetTrack(i);
2716 if (seed->GetNumberOfClusters()<60 && seed->GetNumberOfClusters()<(esd->GetTPCclusters(0) -5)*0.8){
2717 AliExternalTrackParam paramIn;
2718 AliExternalTrackParam paramOut;
2719 Int_t ncl = seed->RefitTrack(seed,¶mIn,¶mOut);
2720 if (AliTPCReconstructor::StreamLevel()>2) {
2721 (*fDebugStreamer)<<"RecoverIn"<<
2725 "pout.="<<¶mOut<<
2730 seed->Set(paramIn.GetX(),paramIn.GetAlpha(),paramIn.GetParameter(),paramIn.GetCovariance());
2731 seed->SetNumberOfClusters(ncl);
2735 seed->PropagateTo(fkParam->GetInnerRadiusLow());
2736 seed->UpdatePoints();
2737 AddCovariance(seed);
2738 MakeESDBitmaps(seed, esd);
2739 seed->CookdEdx(0.02,0.6);
2740 CookLabel(seed,0.1); //For comparison only
2742 if (AliTPCReconstructor::StreamLevel()>1 && seed!=0) {
2743 TTreeSRedirector &cstream = *fDebugStreamer;
2750 if (seed->GetNumberOfClusters()>15){
2751 esd->UpdateTrackParams(seed,AliESDtrack::kTPCrefit);
2752 esd->SetTPCPoints(seed->GetPoints());
2753 esd->SetTPCPointsF(seed->GetNFoundable());
2754 Int_t ndedx = seed->GetNCDEDX(0);
2755 Float_t sdedx = seed->GetSDEDX(0);
2756 Float_t dedx = seed->GetdEdx();
2757 // apply mutliplicity dependent dEdx correction if available
2758 if (graphMultDependenceDeDx) {
2759 Double_t corrGain = AliTPCcalibDButil::EvalGraphConst(graphMultDependenceDeDx, nContribut);
2760 dedx += (1 - corrGain)*50.; // MIP is normalized to 50
2762 esd->SetTPCsignal(dedx, sdedx, ndedx);
2764 // fill new dEdx information
2766 Double32_t signal[4];
2770 for(Int_t iarr=0;iarr<3;iarr++) {
2771 signal[iarr] = seed->GetDEDXregion(iarr+1);
2772 ncl[iarr] = seed->GetNCDEDX(iarr+1);
2773 nrows[iarr] = seed->GetNCDEDXInclThres(iarr+1);
2775 signal[3] = seed->GetDEDXregion(4);
2777 AliTPCdEdxInfo * infoTpcPid = new AliTPCdEdxInfo();
2778 infoTpcPid->SetTPCSignalRegionInfo(signal, ncl, nrows);
2779 esd->SetTPCdEdxInfo(infoTpcPid);
2781 // add seed to the esd track in Calib level
2783 Bool_t storeFriend = gRandom->Rndm()<(kMaxFriendTracks)/Float_t(nseed);
2784 if (AliTPCReconstructor::StreamLevel()>0 &&storeFriend){
2785 // RS: this is the only place where the seed is created not in the pool,
2786 // since it should belong to ESDevent
2787 AliTPCseed * seedCopy = new AliTPCseed(*seed, kTRUE);
2788 esd->AddCalibObject(seedCopy);
2793 //printf("problem\n");
2796 //FindKinks(fSeeds,event);
2797 if (AliTPCReconstructor::StreamLevel()>3) DumpClusters(2,fSeeds);
2798 Info("RefitInward","Number of refitted tracks %d",ntracks);
2803 Int_t AliTPCtrackerMI::PropagateBack(AliESDEvent *event)
2806 // back propagation of ESD tracks
2808 if (!event) return 0;
2812 PropagateBack(fSeeds);
2813 RemoveUsed2(fSeeds,0.4,0.4,20);
2814 //FindCurling(fSeeds, fEvent,1);
2815 FindSplitted(fSeeds, event,1); // find multi found tracks
2816 if (AliTPCReconstructor::StreamLevel()>5) FindMultiMC(fSeeds, fEvent,1); // find multi found tracks
2819 Int_t nseed = fSeeds->GetEntriesFast();
2821 for (Int_t i=0;i<nseed;i++){
2822 AliTPCseed * seed = (AliTPCseed*) fSeeds->UncheckedAt(i);
2823 if (!seed) continue;
2824 if (seed->GetKinkIndex(0)<0) UpdateKinkQualityM(seed); // update quality informations for kinks
2825 seed->UpdatePoints();
2826 AddCovariance(seed);
2827 AliESDtrack *esd=event->GetTrack(i);
2828 if (!esd) continue; //never happen
2829 if (seed->GetNumberOfClusters()<60 && seed->GetNumberOfClusters()<(esd->GetTPCclusters(0) -5)*0.8){
2830 AliExternalTrackParam paramIn;
2831 AliExternalTrackParam paramOut;
2832 Int_t ncl = seed->RefitTrack(seed,¶mIn,¶mOut);
2833 if (AliTPCReconstructor::StreamLevel()>2) {
2834 (*fDebugStreamer)<<"RecoverBack"<<
2838 "pout.="<<¶mOut<<
2843 seed->Set(paramOut.GetX(),paramOut.GetAlpha(),paramOut.GetParameter(),paramOut.GetCovariance());
2844 seed->SetNumberOfClusters(ncl);
2847 seed->CookdEdx(0.02,0.6);
2848 CookLabel(seed,0.1); //For comparison only
2849 if (seed->GetNumberOfClusters()>15){
2850 esd->UpdateTrackParams(seed,AliESDtrack::kTPCout);
2851 esd->SetTPCPoints(seed->GetPoints());
2852 esd->SetTPCPointsF(seed->GetNFoundable());
2853 Int_t ndedx = seed->GetNCDEDX(0);
2854 Float_t sdedx = seed->GetSDEDX(0);
2855 Float_t dedx = seed->GetdEdx();
2856 esd->SetTPCsignal(dedx, sdedx, ndedx);
2858 Int_t eventnumber = event->GetEventNumberInFile();// patch 28 fev 06
2859 // This is most likely NOT the event number you'd like to use. It has nothing to do with the 'real' event number
2860 if (AliTPCReconstructor::StreamLevel()>1 && esd) {
2861 (*fDebugStreamer)<<"Cback"<<
2864 "EventNrInFile="<<eventnumber<<
2869 if (AliTPCReconstructor::StreamLevel()>3) DumpClusters(1,fSeeds);
2870 //FindKinks(fSeeds,event);
2871 Info("PropagateBack","Number of back propagated tracks %d",ntracks);
2878 Int_t AliTPCtrackerMI::PostProcess(AliESDEvent *event)
2881 // Post process events
2883 if (!event) return 0;
2886 // Set TPC event status
2889 // event affected by HV dip
2891 if(IsTPCHVDipEvent(event)) {
2892 event->ResetDetectorStatus(AliDAQ::kTPC);
2895 //printf("Status %d \n", event->IsDetectorOn(AliDAQ::kTPC));
2901 void AliTPCtrackerMI::DeleteSeeds()
2910 void AliTPCtrackerMI::ReadSeeds(const AliESDEvent *const event, Int_t direction)
2913 //read seeds from the event
2915 Int_t nentr=event->GetNumberOfTracks();
2917 Info("ReadSeeds", "Number of ESD tracks: %d\n", nentr);
2922 fSeeds = new TObjArray(nentr);
2926 for (Int_t i=0; i<nentr; i++) {
2927 AliESDtrack *esd=event->GetTrack(i);
2928 ULong_t status=esd->GetStatus();
2929 if (!(status&AliESDtrack::kTPCin)) continue;
2930 AliTPCtrack t(*esd);
2931 t.SetNumberOfClusters(0);
2932 // AliTPCseed *seed = new AliTPCseed(t,t.GetAlpha());
2933 AliTPCseed *seed = new( NextFreeSeed() ) AliTPCseed(t/*,t.GetAlpha()*/);
2934 seed->SetPoolID(fLastSeedID);
2935 seed->SetUniqueID(esd->GetID());
2936 AddCovariance(seed); //add systematic ucertainty
2937 for (Int_t ikink=0;ikink<3;ikink++) {
2938 Int_t index = esd->GetKinkIndex(ikink);
2939 seed->GetKinkIndexes()[ikink] = index;
2940 if (index==0) continue;
2941 index = TMath::Abs(index);
2942 AliESDkink * kink = fEvent->GetKink(index-1);
2943 if (kink&&esd->GetKinkIndex(ikink)<0){
2944 if ((status & AliESDtrack::kTRDrefit) != 0) kink->SetStatus(1,2);
2945 if ((status & AliESDtrack::kITSout) != 0) kink->SetStatus(1,0);
2947 if (kink&&esd->GetKinkIndex(ikink)>0){
2948 if ((status & AliESDtrack::kTRDrefit) != 0) kink->SetStatus(1,6);
2949 if ((status & AliESDtrack::kITSout) != 0) kink->SetStatus(1,4);
2953 if (((status&AliESDtrack::kITSout)==0)&&(direction==1)) seed->ResetCovariance(10.);
2954 if ( direction ==2 &&(status & AliESDtrack::kTRDrefit) == 0 ) seed->ResetCovariance(10.);
2955 //if ( direction ==2 && ((status & AliESDtrack::kTPCout) == 0) ) {
2956 // fSeeds->AddAt(0,i);
2957 // MarkSeedFree( seed );
2960 if ( direction ==2 &&(status & AliESDtrack::kTRDrefit) > 0 ) {
2961 Double_t par0[5],par1[5],alpha,x;
2962 esd->GetInnerExternalParameters(alpha,x,par0);
2963 esd->GetExternalParameters(x,par1);
2964 Double_t delta1 = TMath::Abs(par0[4]-par1[4])/(0.000000001+TMath::Abs(par0[4]+par1[4]));
2965 Double_t delta2 = TMath::Abs(par0[3]-par1[3]);
2967 if (esd->GetTRDncls()>0) trdchi2 = esd->GetTRDchi2()/esd->GetTRDncls();
2968 //reset covariance if suspicious
2969 if ( (delta1>0.1) || (delta2>0.006) ||trdchi2>7.)
2970 seed->ResetCovariance(10.);
2975 // rotate to the local coordinate system
2977 fSectors=fInnerSec; fN=fkNIS;
2978 Double_t alpha=seed->GetAlpha();
2979 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
2980 if (alpha < 0. ) alpha += 2.*TMath::Pi();
2981 Int_t ns=Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN;
2982 alpha =ns*fSectors->GetAlpha() + fSectors->GetAlphaShift();
2983 alpha-=seed->GetAlpha();
2984 if (alpha<-TMath::Pi()) alpha += 2*TMath::Pi();
2985 if (alpha>=TMath::Pi()) alpha -= 2*TMath::Pi();
2986 if (TMath::Abs(alpha) > 0.001) { // This should not happen normally
2987 AliWarning(Form("Rotating track over %f",alpha));
2988 if (!seed->Rotate(alpha)) {
2989 MarkSeedFree( seed );
2995 if (esd->GetKinkIndex(0)<=0){
2996 for (Int_t irow=0;irow<160;irow++){
2997 Int_t index = seed->GetClusterIndex2(irow);
3000 AliTPCclusterMI * cl = GetClusterMI(index);
3001 seed->SetClusterPointer(irow,cl);
3003 if ((index & 0x8000)==0){
3004 cl->Use(10); // accepted cluster
3006 cl->Use(6); // close cluster not accepted
3009 Info("ReadSeeds","Not found cluster");
3014 fSeeds->AddAt(seed,i);
3020 //_____________________________________________________________________________
3021 void AliTPCtrackerMI::MakeSeeds3(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, Float_t cuts[4],
3022 Float_t deltay, Int_t ddsec) {
3023 //-----------------------------------------------------------------
3024 // This function creates track seeds.
3025 // SEEDING WITH VERTEX CONSTRAIN
3026 //-----------------------------------------------------------------
3027 // cuts[0] - fP4 cut
3028 // cuts[1] - tan(phi) cut
3029 // cuts[2] - zvertex cut
3030 // cuts[3] - fP3 cut
3038 Double_t x[5], c[15];
3039 // Int_t di = i1-i2;
3041 AliTPCseed * seed = new( NextFreeSeed() ) AliTPCseed();
3042 seed->SetPoolID(fLastSeedID);
3043 Double_t alpha=fSectors->GetAlpha(), shift=fSectors->GetAlphaShift();
3044 Double_t cs=cos(alpha), sn=sin(alpha);
3046 // Double_t x1 =fOuterSec->GetX(i1);
3047 //Double_t xx2=fOuterSec->GetX(i2);
3049 Double_t x1 =GetXrow(i1);
3050 Double_t xx2=GetXrow(i2);
3052 Double_t x3=GetX(), y3=GetY(), z3=GetZ();
3054 Int_t imiddle = (i2+i1)/2; //middle pad row index
3055 Double_t xm = GetXrow(imiddle); // radius of middle pad-row
3056 const AliTPCtrackerRow& krm=GetRow(sec,imiddle); //middle pad -row
3060 const AliTPCtrackerRow& kr1=GetRow(ns,i1);
3061 Double_t ymax = GetMaxY(i1)-kr1.GetDeadZone()-1.5;
3062 Double_t ymaxm = GetMaxY(imiddle)-kr1.GetDeadZone()-1.5;
3065 // change cut on curvature if it can't reach this layer
3066 // maximal curvature set to reach it
3067 Double_t dvertexmax = TMath::Sqrt((x1-x3)*(x1-x3)+(ymax+5-y3)*(ymax+5-y3));
3068 if (dvertexmax*0.5*cuts[0]>0.85){
3069 cuts[0] = 0.85/(dvertexmax*0.5+1.);
3071 Double_t r2min = 1/(cuts[0]*cuts[0]); //minimal square of radius given by cut
3074 if (deltay>0) ddsec = 0;
3075 // loop over clusters
3076 for (Int_t is=0; is < kr1; is++) {
3078 if (kr1[is]->IsUsed(10)) continue;
3079 Double_t y1=kr1[is]->GetY(), z1=kr1[is]->GetZ();
3080 //if (TMath::Abs(y1)>ymax) continue;
3082 if (deltay>0 && TMath::Abs(ymax-TMath::Abs(y1))> deltay ) continue; // seed only at the edge
3084 // find possible directions
3085 Float_t anglez = (z1-z3)/(x1-x3);
3086 Float_t extraz = z1 - anglez*(x1-xx2); // extrapolated z
3089 //find rotation angles relative to line given by vertex and point 1
3090 Double_t dvertex2 = (x1-x3)*(x1-x3)+(y1-y3)*(y1-y3);
3091 Double_t dvertex = TMath::Sqrt(dvertex2);
3092 Double_t angle13 = TMath::ATan((y1-y3)/(x1-x3));
3093 Double_t cs13 = cos(-angle13), sn13 = sin(-angle13);
3096 // loop over 2 sectors
3102 Double_t dddz1=0; // direction of delta inclination in z axis
3109 for (Int_t dsec = dsec1; dsec<=dsec2;dsec++){
3110 Int_t sec2 = sec + dsec;
3112 // AliTPCtrackerRow& kr2 = fOuterSec[(sec2+fkNOS)%fkNOS][i2];
3113 //AliTPCtrackerRow& kr2m = fOuterSec[(sec2+fkNOS)%fkNOS][imiddle];
3114 AliTPCtrackerRow& kr2 = GetRow((sec2+fkNOS)%fkNOS,i2);
3115 AliTPCtrackerRow& kr2m = GetRow((sec2+fkNOS)%fkNOS,imiddle);
3116 Int_t index1 = TMath::Max(kr2.Find(extraz-0.6-dddz1*TMath::Abs(z1)*0.05)-1,0);
3117 Int_t index2 = TMath::Min(kr2.Find(extraz+0.6+dddz2*TMath::Abs(z1)*0.05)+1,kr2);
3119 // rotation angles to p1-p3
3120 Double_t cs13r = cos(-angle13+dsec*alpha)/dvertex, sn13r = sin(-angle13+dsec*alpha)/dvertex;
3121 Double_t x2, y2, z2;
3123 // Double_t dymax = maxangle*TMath::Abs(x1-xx2);
3126 Double_t dxx0 = (xx2-x3)*cs13r;
3127 Double_t dyy0 = (xx2-x3)*sn13r;
3128 for (Int_t js=index1; js < index2; js++) {
3129 const AliTPCclusterMI *kcl = kr2[js];
3130 if (kcl->IsUsed(10)) continue;
3132 //calcutate parameters
3134 Double_t yy0 = dyy0 +(kcl->GetY()-y3)*cs13r;
3136 if (TMath::Abs(yy0)<0.000001) continue;
3137 Double_t xx0 = dxx0 -(kcl->GetY()-y3)*sn13r;
3138 Double_t y0 = 0.5*(xx0*xx0+yy0*yy0-xx0)/yy0;
3139 Double_t r02 = (0.25+y0*y0)*dvertex2;
3140 //curvature (radius) cut
3141 if (r02<r2min) continue;
3145 Double_t c0 = 1/TMath::Sqrt(r02);
3149 //Double_t dfi0 = 2.*TMath::ASin(dvertex*c0*0.5);
3150 //Double_t dfi1 = 2.*TMath::ASin(TMath::Sqrt(yy0*yy0+(1-xx0)*(1-xx0))*dvertex*c0*0.5);
3151 Double_t dfi0 = 2.*AliTPCFastMath::FastAsin(dvertex*c0*0.5);
3152 Double_t dfi1 = 2.*AliTPCFastMath::FastAsin(TMath::Sqrt(yy0*yy0+(1-xx0)*(1-xx0))*dvertex*c0*0.5);
3155 Double_t z0 = kcl->GetZ();
3156 Double_t zzzz2 = z1-(z1-z3)*dfi1/dfi0;
3157 if (TMath::Abs(zzzz2-z0)>0.5) continue;
3160 Double_t dip = (z1-z0)*c0/dfi1;
3161 Double_t x0 = (0.5*cs13+y0*sn13)*dvertex*c0;
3172 x2= xx2*cs-y2*sn*dsec;
3173 y2=+xx2*sn*dsec+y2*cs;
3183 // do we have cluster at the middle ?
3185 GetProlongation(x1,xm,x,ym,zm);
3187 AliTPCclusterMI * cm=0;
3188 if (TMath::Abs(ym)-ymaxm<0){
3189 cm = krm.FindNearest2(ym,zm,1.0,0.6,dummy);
3190 if ((!cm) || (cm->IsUsed(10))) {
3195 // rotate y1 to system 0
3196 // get state vector in rotated system
3197 Double_t yr1 = (-0.5*sn13+y0*cs13)*dvertex*c0;
3198 Double_t xr2 = x0*cs+yr1*sn*dsec;
3199 Double_t xr[5]={kcl->GetY(),kcl->GetZ(), xr2, dip, c0};
3201 GetProlongation(xx2,xm,xr,ym,zm);
3202 if (TMath::Abs(ym)-ymaxm<0){
3203 cm = kr2m.FindNearest2(ym,zm,1.0,0.6,dummy);
3204 if ((!cm) || (cm->IsUsed(10))) {
3214 dym = ym - cm->GetY();
3215 dzm = zm - cm->GetZ();
3222 Double_t sy1=kr1[is]->GetSigmaY2()*2., sz1=kr1[is]->GetSigmaZ2()*2.;
3223 Double_t sy2=kcl->GetSigmaY2()*2., sz2=kcl->GetSigmaZ2()*2.;
3224 //Double_t sy3=400*3./12., sy=0.1, sz=0.1;
3225 Double_t sy3=25000*x[4]*x[4]+0.1, sy=0.1, sz=0.1;
3226 //Double_t sy3=25000*x[4]*x[4]*60+0.5, sy=0.1, sz=0.1;
3228 Double_t f40=(F1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
3229 Double_t f42=(F1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
3230 Double_t f43=(F1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
3231 Double_t f20=(F2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
3232 Double_t f22=(F2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
3233 Double_t f23=(F2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
3235 Double_t f30=(F3(x1,y1+sy,x2,y2,z1,z2)-x[3])/sy;
3236 Double_t f31=(F3(x1,y1,x2,y2,z1+sz,z2)-x[3])/sz;
3237 Double_t f32=(F3(x1,y1,x2,y2+sy,z1,z2)-x[3])/sy;
3238 Double_t f34=(F3(x1,y1,x2,y2,z1,z2+sz)-x[3])/sz;
3242 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3243 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3244 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3245 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3246 c[13]=f30*sy1*f40+f32*sy2*f42;
3247 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3249 // if (!BuildSeed(kr1[is],kcl,0,x1,x2,x3,x,c)) continue;
3251 UInt_t index=kr1.GetIndex(is);
3252 if (seed) {MarkSeedFree(seed); seed = 0;}
3253 AliTPCseed *track = seed = new( NextFreeSeed() ) AliTPCseed(x1, ns*alpha+shift, x, c, index);
3254 seed->SetPoolID(fLastSeedID);
3255 track->SetIsSeeding(kTRUE);
3256 track->SetSeed1(i1);
3257 track->SetSeed2(i2);
3258 track->SetSeedType(3);
3262 FollowProlongation(*track, (i1+i2)/2,1);
3263 Int_t foundable,found,shared;
3264 track->GetClusterStatistic((i1+i2)/2,i1, found, foundable, shared, kTRUE);
3265 if ((found<0.55*foundable) || shared>0.5*found || (track->GetSigmaY2()+track->GetSigmaZ2())>0.5){
3266 MarkSeedFree(seed); seed = 0;
3272 FollowProlongation(*track, i2,1);
3276 track->SetBConstrain(1);
3277 // track->fLastPoint = i1+fInnerSec->GetNRows(); // first cluster in track position
3278 track->SetLastPoint(i1); // first cluster in track position
3279 track->SetFirstPoint(track->GetLastPoint());
3281 if (track->GetNumberOfClusters()<(i1-i2)*0.5 ||
3282 track->GetNumberOfClusters() < track->GetNFoundable()*0.6 ||
3283 track->GetNShared()>0.4*track->GetNumberOfClusters() ) {
3284 MarkSeedFree(seed); seed = 0;
3288 // Z VERTEX CONDITION
3289 Double_t zv, bz=GetBz();
3290 if ( !track->GetZAt(0.,bz,zv) ) continue;
3291 if (TMath::Abs(zv-z3)>cuts[2]) {
3292 FollowProlongation(*track, TMath::Max(i2-20,0));
3293 if ( !track->GetZAt(0.,bz,zv) ) continue;
3294 if (TMath::Abs(zv-z3)>cuts[2]){
3295 FollowProlongation(*track, TMath::Max(i2-40,0));
3296 if ( !track->GetZAt(0.,bz,zv) ) continue;
3297 if (TMath::Abs(zv-z3)>cuts[2] &&(track->GetNumberOfClusters() > track->GetNFoundable()*0.7)){
3298 // make seed without constrain
3299 AliTPCseed * track2 = MakeSeed(track,0.2,0.5,1.);
3300 FollowProlongation(*track2, i2,1);
3301 track2->SetBConstrain(kFALSE);
3302 track2->SetSeedType(1);
3303 arr->AddLast(track2);
3304 MarkSeedFree( seed ); seed = 0;
3308 MarkSeedFree( seed ); seed = 0;
3315 track->SetSeedType(0);
3316 arr->AddLast(track); // note, track is seed, don't free the seed
3317 seed = new( NextFreeSeed() ) AliTPCseed;
3318 seed->SetPoolID(fLastSeedID);
3320 // don't consider other combinations
3321 if (track->GetNumberOfClusters() > track->GetNFoundable()*0.8)
3327 Info("MakeSeeds3","\nSeeding statistic:\t%d\t%d\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin,nout1,nout2);
3329 if (seed) MarkSeedFree( seed );
3333 void AliTPCtrackerMI::MakeSeeds5(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, Float_t cuts[4],
3338 //-----------------------------------------------------------------
3339 // This function creates track seeds.
3340 //-----------------------------------------------------------------
3341 // cuts[0] - fP4 cut
3342 // cuts[1] - tan(phi) cut
3343 // cuts[2] - zvertex cut
3344 // cuts[3] - fP3 cut
3354 Double_t x[5], c[15];
3356 // make temporary seed
3357 AliTPCseed * seed = new( NextFreeSeed() ) AliTPCseed;
3358 seed->SetPoolID(fLastSeedID);
3359 Double_t alpha=fOuterSec->GetAlpha(), shift=fOuterSec->GetAlphaShift();
3360 // Double_t cs=cos(alpha), sn=sin(alpha);
3365 Double_t x1 = GetXrow(i1-1);
3366 const AliTPCtrackerRow& kr1=GetRow(sec,i1-1);
3367 Double_t y1max = GetMaxY(i1-1)-kr1.GetDeadZone()-1.5;
3369 Double_t x1p = GetXrow(i1);
3370 const AliTPCtrackerRow& kr1p=GetRow(sec,i1);
3372 Double_t x1m = GetXrow(i1-2);
3373 const AliTPCtrackerRow& kr1m=GetRow(sec,i1-2);
3376 //last 3 padrow for seeding
3377 AliTPCtrackerRow& kr3 = GetRow((sec+fkNOS)%fkNOS,i1-7);
3378 Double_t x3 = GetXrow(i1-7);
3379 // Double_t y3max= GetMaxY(i1-7)-kr3.fDeadZone-1.5;
3381 AliTPCtrackerRow& kr3p = GetRow((sec+fkNOS)%fkNOS,i1-6);
3382 Double_t x3p = GetXrow(i1-6);
3384 AliTPCtrackerRow& kr3m = GetRow((sec+fkNOS)%fkNOS,i1-8);
3385 Double_t x3m = GetXrow(i1-8);
3390 Int_t im = i1-4; //middle pad row index
3391 Double_t xm = GetXrow(im); // radius of middle pad-row
3392 const AliTPCtrackerRow& krm=GetRow(sec,im); //middle pad -row
3393 // Double_t ymmax = GetMaxY(im)-kr1.fDeadZone-1.5;
3396 Double_t deltax = x1-x3;
3397 Double_t dymax = deltax*cuts[1];
3398 Double_t dzmax = deltax*cuts[3];
3400 // loop over clusters
3401 for (Int_t is=0; is < kr1; is++) {
3403 if (kr1[is]->IsUsed(10)) continue;
3404 Double_t y1=kr1[is]->GetY(), z1=kr1[is]->GetZ();
3406 if (deltay>0 && TMath::Abs(y1max-TMath::Abs(y1))> deltay ) continue; // seed only at the edge
3408 Int_t index1 = TMath::Max(kr3.Find(z1-dzmax)-1,0);
3409 Int_t index2 = TMath::Min(kr3.Find(z1+dzmax)+1,kr3);
3415 for (Int_t js=index1; js < index2; js++) {
3416 const AliTPCclusterMI *kcl = kr3[js];
3417 if (kcl->IsUsed(10)) continue;
3419 // apply angular cuts
3420 if (TMath::Abs(y1-y3)>dymax) continue;
3423 if (TMath::Abs(z1-z3)>dzmax) continue;
3425 Double_t angley = (y1-y3)/(x1-x3);
3426 Double_t anglez = (z1-z3)/(x1-x3);
3428 Double_t erry = TMath::Abs(angley)*(x1-x1m)*0.5+0.5;
3429 Double_t errz = TMath::Abs(anglez)*(x1-x1m)*0.5+0.5;
3431 Double_t yyym = angley*(xm-x1)+y1;
3432 Double_t zzzm = anglez*(xm-x1)+z1;
3434 const AliTPCclusterMI *kcm = krm.FindNearest2(yyym,zzzm,erry,errz,index);
3436 if (kcm->IsUsed(10)) continue;
3438 erry = TMath::Abs(angley)*(x1-x1m)*0.4+0.5;
3439 errz = TMath::Abs(anglez)*(x1-x1m)*0.4+0.5;
3446 // look around first
3447 const AliTPCclusterMI *kc1m = kr1m.FindNearest2(angley*(x1m-x1)+y1,
3453 if (kc1m->IsUsed(10)) used++;
3455 const AliTPCclusterMI *kc1p = kr1p.FindNearest2(angley*(x1p-x1)+y1,
3461 if (kc1p->IsUsed(10)) used++;
3463 if (used>1) continue;
3464 if (found<1) continue;
3468 const AliTPCclusterMI *kc3m = kr3m.FindNearest2(angley*(x3m-x3)+y3,
3474 if (kc3m->IsUsed(10)) used++;
3478 const AliTPCclusterMI *kc3p = kr3p.FindNearest2(angley*(x3p-x3)+y3,
3484 if (kc3p->IsUsed(10)) used++;
3488 if (used>1) continue;
3489 if (found<3) continue;
3499 x[4]=F1(x1,y1,x2,y2,x3,y3);
3500 //if (TMath::Abs(x[4]) >= cuts[0]) continue;
3503 x[2]=F2(x1,y1,x2,y2,x3,y3);
3506 x[3]=F3n(x1,y1,x2,y2,z1,z2,x[4]);
3507 //if (TMath::Abs(x[3]) > cuts[3]) continue;
3511 Double_t sy1=0.1, sz1=0.1;
3512 Double_t sy2=0.1, sz2=0.1;
3513 Double_t sy3=0.1, sy=0.1, sz=0.1;
3515 Double_t f40=(F1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
3516 Double_t f42=(F1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
3517 Double_t f43=(F1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
3518 Double_t f20=(F2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
3519 Double_t f22=(F2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
3520 Double_t f23=(F2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
3522 Double_t f30=(F3(x1,y1+sy,x2,y2,z1,z2)-x[3])/sy;
3523 Double_t f31=(F3(x1,y1,x2,y2,z1+sz,z2)-x[3])/sz;
3524 Double_t f32=(F3(x1,y1,x2,y2+sy,z1,z2)-x[3])/sy;
3525 Double_t f34=(F3(x1,y1,x2,y2,z1,z2+sz)-x[3])/sz;
3529 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3530 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3531 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3532 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3533 c[13]=f30*sy1*f40+f32*sy2*f42;
3534 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3536 // if (!BuildSeed(kr1[is],kcl,0,x1,x2,x3,x,c)) continue;
3538 index=kr1.GetIndex(is);
3539 if (seed) {MarkSeedFree( seed ); seed = 0;}
3540 AliTPCseed *track = seed = new( NextFreeSeed() ) AliTPCseed(x1, sec*alpha+shift, x, c, index);
3541 seed->SetPoolID(fLastSeedID);
3543 track->SetIsSeeding(kTRUE);
3546 FollowProlongation(*track, i1-7,1);
3547 if (track->GetNumberOfClusters() < track->GetNFoundable()*0.75 ||
3548 track->GetNShared()>0.6*track->GetNumberOfClusters() || ( track->GetSigmaY2()+ track->GetSigmaZ2())>0.6){
3549 MarkSeedFree( seed ); seed = 0;
3555 FollowProlongation(*track, i2,1);
3556 track->SetBConstrain(0);
3557 track->SetLastPoint(i1+fInnerSec->GetNRows()); // first cluster in track position
3558 track->SetFirstPoint(track->GetLastPoint());
3560 if (track->GetNumberOfClusters()<(i1-i2)*0.5 ||
3561 track->GetNumberOfClusters()<track->GetNFoundable()*0.7 ||
3562 track->GetNShared()>2. || track->GetChi2()/track->GetNumberOfClusters()>6 || ( track->GetSigmaY2()+ track->GetSigmaZ2())>0.5 ) {
3563 MarkSeedFree( seed ); seed = 0;
3568 FollowProlongation(*track, TMath::Max(i2-10,0),1);
3569 AliTPCseed * track2 = MakeSeed(track,0.2,0.5,0.9);
3570 FollowProlongation(*track2, i2,1);
3571 track2->SetBConstrain(kFALSE);
3572 track2->SetSeedType(4);
3573 arr->AddLast(track2);
3574 MarkSeedFree( seed ); seed = 0;
3578 //arr->AddLast(track);
3579 //seed = new AliTPCseed;
3585 Info("MakeSeeds5","\nSeeding statiistic:\t%d\t%d\t%d\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin,nout1,nout2,nout3);
3587 if (seed) MarkSeedFree(seed);
3591 //_____________________________________________________________________________
3592 void AliTPCtrackerMI::MakeSeeds2(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, Float_t */*cuts[4]*/,
3593 Float_t deltay, Bool_t /*bconstrain*/) {
3594 //-----------------------------------------------------------------
3595 // This function creates track seeds - without vertex constraint
3596 //-----------------------------------------------------------------
3597 // cuts[0] - fP4 cut - not applied
3598 // cuts[1] - tan(phi) cut
3599 // cuts[2] - zvertex cut - not applied
3600 // cuts[3] - fP3 cut
3610 Double_t alpha=fOuterSec->GetAlpha(), shift=fOuterSec->GetAlphaShift();
3611 // Double_t cs=cos(alpha), sn=sin(alpha);
3612 Int_t row0 = (i1+i2)/2;
3613 Int_t drow = (i1-i2)/2;
3614 const AliTPCtrackerRow& kr0=fSectors[sec][row0];
3615 AliTPCtrackerRow * kr=0;
3617 AliTPCpolyTrack polytrack;
3618 Int_t nclusters=fSectors[sec][row0];
3619 AliTPCseed * seed = new( NextFreeSeed() ) AliTPCseed;
3620 seed->SetPoolID(fLastSeedID);
3625 for (Int_t is=0; is < nclusters; is++) { //LOOP over clusters
3627 Int_t nfoundable =0;
3628 for (Int_t iter =1; iter<2; iter++){ //iterations
3629 const AliTPCtrackerRow& krm=fSectors[sec][row0-iter];
3630 const AliTPCtrackerRow& krp=fSectors[sec][row0+iter];
3631 const AliTPCclusterMI * cl= kr0[is];
3633 if (cl->IsUsed(10)) {
3639 Double_t x = kr0.GetX();
3640 // Initialization of the polytrack
3645 Double_t y0= cl->GetY();
3646 Double_t z0= cl->GetZ();
3650 Double_t ymax = fSectors->GetMaxY(row0)-kr0.GetDeadZone()-1.5;
3651 if (deltay>0 && TMath::Abs(ymax-TMath::Abs(y0))> deltay ) continue; // seed only at the edge
3653 erry = (0.5)*cl->GetSigmaY2()/TMath::Sqrt(cl->GetQ())*6;
3654 errz = (0.5)*cl->GetSigmaZ2()/TMath::Sqrt(cl->GetQ())*6;
3655 polytrack.AddPoint(x,y0,z0,erry, errz);
3658 if (cl->IsUsed(10)) sumused++;
3661 Float_t roady = (5*TMath::Sqrt(cl->GetSigmaY2()+0.2)+1.)*iter;
3662 Float_t roadz = (5*TMath::Sqrt(cl->GetSigmaZ2()+0.2)+1.)*iter;
3665 AliTPCclusterMI * cl1 = krm.FindNearest(y0,z0,roady,roadz);
3666 if (cl1 && TMath::Abs(ymax-TMath::Abs(y0))) {
3667 erry = (0.5)*cl1->GetSigmaY2()/TMath::Sqrt(cl1->GetQ())*3;
3668 errz = (0.5)*cl1->GetSigmaZ2()/TMath::Sqrt(cl1->GetQ())*3;
3669 if (cl1->IsUsed(10)) sumused++;
3670 polytrack.AddPoint(x,cl1->GetY(),cl1->GetZ(),erry,errz);
3674 AliTPCclusterMI * cl2 = krp.FindNearest(y0,z0,roady,roadz);
3676 erry = (0.5)*cl2->GetSigmaY2()/TMath::Sqrt(cl2->GetQ())*3;
3677 errz = (0.5)*cl2->GetSigmaZ2()/TMath::Sqrt(cl2->GetQ())*3;
3678 if (cl2->IsUsed(10)) sumused++;
3679 polytrack.AddPoint(x,cl2->GetY(),cl2->GetZ(),erry,errz);
3682 if (sumused>0) continue;
3684 polytrack.UpdateParameters();
3690 nfoundable = polytrack.GetN();
3691 nfound = nfoundable;
3693 for (Int_t ddrow = iter+1; ddrow<drow;ddrow++){
3694 Float_t maxdist = 0.8*(1.+3./(ddrow));
3695 for (Int_t delta = -1;delta<=1;delta+=2){
3696 Int_t row = row0+ddrow*delta;
3697 kr = &(fSectors[sec][row]);
3698 Double_t xn = kr->GetX();
3699 Double_t ymax1 = fSectors->GetMaxY(row)-kr->GetDeadZone()-1.5;
3700 polytrack.GetFitPoint(xn,yn,zn);
3701 if (TMath::Abs(yn)>ymax1) continue;
3703 AliTPCclusterMI * cln = kr->FindNearest(yn,zn,roady,roadz);
3705 Float_t dist = TMath::Sqrt( (yn-cln->GetY())*(yn-cln->GetY())+(zn-cln->GetZ())*(zn-cln->GetZ()));
3708 erry = (dist+0.3)*cln->GetSigmaY2()/TMath::Sqrt(cln->GetQ())*(1.+1./(ddrow));
3709 errz = (dist+0.3)*cln->GetSigmaZ2()/TMath::Sqrt(cln->GetQ())*(1.+1./(ddrow));
3710 if (cln->IsUsed(10)) {
3711 // printf("used\n");
3719 polytrack.AddPoint(xn,cln->GetY(),cln->GetZ(),erry, errz);
3724 if ( (sumused>3) || (sumused>0.5*nfound) || (nfound<0.6*nfoundable)) break;
3725 polytrack.UpdateParameters();
3728 if ( (sumused>3) || (sumused>0.5*nfound)) {
3729 //printf("sumused %d\n",sumused);
3734 polytrack.GetFitDerivation(kr0.GetX(),dy,dz);
3735 AliTPCpolyTrack track2;
3737 polytrack.Refit(track2,0.5+TMath::Abs(dy)*0.3,0.4+TMath::Abs(dz)*0.3);
3738 if (track2.GetN()<0.5*nfoundable) continue;
3741 if ((nfound>0.6*nfoundable) &&( nfoundable>0.4*(i1-i2))) {
3743 // test seed with and without constrain
3744 for (Int_t constrain=0; constrain<=0;constrain++){
3745 // add polytrack candidate
3747 Double_t x[5], c[15];
3748 Double_t x1,x2,x3,y1,y2,y3,z1,z2,z3;
3749 track2.GetBoundaries(x3,x1);
3751 track2.GetFitPoint(x1,y1,z1);
3752 track2.GetFitPoint(x2,y2,z2);
3753 track2.GetFitPoint(x3,y3,z3);
3755 //is track pointing to the vertex ?
3758 polytrack.GetFitPoint(x0,y0,z0);
3771 x[4]=F1(x1,y1,x2,y2,x3,y3);
3773 // if (TMath::Abs(x[4]) >= cuts[0]) continue; //
3774 x[2]=F2(x1,y1,x2,y2,x3,y3);
3776 //if (TMath::Abs(x[4]*x1-x[2]) >= cuts[1]) continue;
3777 //x[3]=F3(x1,y1,x2,y2,z1,z2);
3778 x[3]=F3n(x1,y1,x3,y3,z1,z3,x[4]);
3779 //if (TMath::Abs(x[3]) > cuts[3]) continue;
3782 Double_t sy =0.1, sz =0.1;
3783 Double_t sy1=0.02, sz1=0.02;
3784 Double_t sy2=0.02, sz2=0.02;
3788 sy3=25000*x[4]*x[4]+0.1, sy=0.1, sz=0.1;
3791 Double_t f40=(F1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
3792 Double_t f42=(F1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
3793 Double_t f43=(F1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
3794 Double_t f20=(F2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
3795 Double_t f22=(F2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
3796 Double_t f23=(F2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
3798 Double_t f30=(F3(x1,y1+sy,x3,y3,z1,z3)-x[3])/sy;
3799 Double_t f31=(F3(x1,y1,x3,y3,z1+sz,z3)-x[3])/sz;
3800 Double_t f32=(F3(x1,y1,x3,y3+sy,z1,z3)-x[3])/sy;
3801 Double_t f34=(F3(x1,y1,x3,y3,z1,z3+sz)-x[3])/sz;
3806 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3807 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3808 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3809 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3810 c[13]=f30*sy1*f40+f32*sy2*f42;
3811 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3813 //Int_t row1 = fSectors->GetRowNumber(x1);
3814 Int_t row1 = GetRowNumber(x1);
3818 if (seed) {MarkSeedFree( seed ); seed = 0;}
3819 AliTPCseed *track = seed = new( NextFreeSeed() ) AliTPCseed(x1,sec*alpha+shift,x,c,index);
3820 seed->SetPoolID(fLastSeedID);
3821 track->SetIsSeeding(kTRUE);
3822 Int_t rc=FollowProlongation(*track, i2);
3823 if (constrain) track->SetBConstrain(1);
3825 track->SetBConstrain(0);
3826 track->SetLastPoint(row1+fInnerSec->GetNRows()); // first cluster in track position
3827 track->SetFirstPoint(track->GetLastPoint());
3829 if (rc==0 || track->GetNumberOfClusters()<(i1-i2)*0.5 ||
3830 track->GetNumberOfClusters() < track->GetNFoundable()*0.6 ||
3831 track->GetNShared()>0.4*track->GetNumberOfClusters()) {
3832 MarkSeedFree( seed ); seed = 0;
3835 arr->AddLast(track); // track IS seed, don't free seed
3836 seed = new( NextFreeSeed() ) AliTPCseed;
3837 seed->SetPoolID(fLastSeedID);
3841 } // if accepted seed
3844 Info("MakeSeeds2","\nSeeding statiistic:\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin3);
3846 if (seed) MarkSeedFree( seed );
3850 AliTPCseed *AliTPCtrackerMI::MakeSeed(AliTPCseed *const track, Float_t r0, Float_t r1, Float_t r2)
3854 //reseed using track points
3855 Int_t p0 = int(r0*track->GetNumberOfClusters()); // point 0
3856 Int_t p1 = int(r1*track->GetNumberOfClusters());
3857 Int_t p2 = int(r2*track->GetNumberOfClusters()); // last point
3859 Double_t x0[3],x1[3],x2[3];
3860 for (Int_t i=0;i<3;i++){
3866 // find track position at given ratio of the length
3867 Int_t sec0=0, sec1=0, sec2=0;
3870 for (Int_t i=0;i<160;i++){
3871 if (track->GetClusterPointer(i)){
3873 AliTPCTrackerPoint *trpoint =track->GetTrackPoint(i);
3874 if ( (index<p0) || x0[0]<0 ){
3875 if (trpoint->GetX()>1){
3876 clindex = track->GetClusterIndex2(i);
3878 x0[0] = trpoint->GetX();
3879 x0[1] = trpoint->GetY();
3880 x0[2] = trpoint->GetZ();
3881 sec0 = ((clindex&0xff000000)>>24)%18;
3886 if ( (index<p1) &&(trpoint->GetX()>1)){
3887 clindex = track->GetClusterIndex2(i);
3889 x1[0] = trpoint->GetX();
3890 x1[1] = trpoint->GetY();
3891 x1[2] = trpoint->GetZ();
3892 sec1 = ((clindex&0xff000000)>>24)%18;
3895 if ( (index<p2) &&(trpoint->GetX()>1)){
3896 clindex = track->GetClusterIndex2(i);
3898 x2[0] = trpoint->GetX();
3899 x2[1] = trpoint->GetY();
3900 x2[2] = trpoint->GetZ();
3901 sec2 = ((clindex&0xff000000)>>24)%18;
3908 Double_t alpha, cs,sn, xx2,yy2;
3910 alpha = (sec1-sec2)*fSectors->GetAlpha();
3911 cs = TMath::Cos(alpha);
3912 sn = TMath::Sin(alpha);
3913 xx2= x1[0]*cs-x1[1]*sn;
3914 yy2= x1[0]*sn+x1[1]*cs;
3918 alpha = (sec0-sec2)*fSectors->GetAlpha();
3919 cs = TMath::Cos(alpha);
3920 sn = TMath::Sin(alpha);
3921 xx2= x0[0]*cs-x0[1]*sn;
3922 yy2= x0[0]*sn+x0[1]*cs;
3928 Double_t x[5],c[15];
3932 x[4]=F1(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]);
3933 // if (x[4]>1) return 0;
3934 x[2]=F2(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]);
3935 x[3]=F3n(x2[0],x2[1],x0[0],x0[1],x2[2],x0[2],x[4]);
3936 //if (TMath::Abs(x[3]) > 2.2) return 0;
3937 //if (TMath::Abs(x[2]) > 1.99) return 0;
3939 Double_t sy =0.1, sz =0.1;
3941 Double_t sy1=0.02+track->GetSigmaY2(), sz1=0.02+track->GetSigmaZ2();
3942 Double_t sy2=0.01+track->GetSigmaY2(), sz2=0.01+track->GetSigmaZ2();
3943 Double_t sy3=0.01+track->GetSigmaY2();
3945 Double_t f40=(F1(x2[0],x2[1]+sy,x1[0],x1[1],x0[0],x0[1])-x[4])/sy;
3946 Double_t f42=(F1(x2[0],x2[1],x1[0],x1[1]+sy,x0[0],x0[1])-x[4])/sy;
3947 Double_t f43=(F1(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]+sy)-x[4])/sy;
3948 Double_t f20=(F2(x2[0],x2[1]+sy,x1[0],x1[1],x0[0],x0[1])-x[2])/sy;
3949 Double_t f22=(F2(x2[0],x2[1],x1[0],x1[1]+sy,x0[0],x0[1])-x[2])/sy;
3950 Double_t f23=(F2(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]+sy)-x[2])/sy;
3952 Double_t f30=(F3(x2[0],x2[1]+sy,x0[0],x0[1],x2[2],x0[2])-x[3])/sy;
3953 Double_t f31=(F3(x2[0],x2[1],x0[0],x0[1],x2[2]+sz,x0[2])-x[3])/sz;
3954 Double_t f32=(F3(x2[0],x2[1],x0[0],x0[1]+sy,x2[2],x0[2])-x[3])/sy;
3955 Double_t f34=(F3(x2[0],x2[1],x0[0],x0[1],x2[2],x0[2]+sz)-x[3])/sz;
3960 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3961 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3962 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3963 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3964 c[13]=f30*sy1*f40+f32*sy2*f42;
3965 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3967 // Int_t row1 = fSectors->GetRowNumber(x2[0]);
3968 AliTPCseed *seed = new( NextFreeSeed() ) AliTPCseed(x2[0], sec2*fSectors->GetAlpha()+fSectors->GetAlphaShift(), x, c, 0);
3969 seed->SetPoolID(fLastSeedID);
3970 // Double_t y0,z0,y1,z1, y2,z2;
3971 //seed->GetProlongation(x0[0],y0,z0);
3972 // seed->GetProlongation(x1[0],y1,z1);
3973 //seed->GetProlongation(x2[0],y2,z2);
3975 seed->SetLastPoint(pp2);
3976 seed->SetFirstPoint(pp2);
3983 AliTPCseed *AliTPCtrackerMI::ReSeed(const AliTPCseed *track, Float_t r0, Float_t r1, Float_t r2)
3987 //reseed using founded clusters
3989 // Find the number of clusters
3990 Int_t nclusters = 0;
3991 for (Int_t irow=0;irow<160;irow++){
3992 if (track->GetClusterIndex(irow)>0) nclusters++;
3996 ipos[0] = TMath::Max(int(r0*nclusters),0); // point 0 cluster
3997 ipos[1] = TMath::Min(int(r1*nclusters),nclusters-1); //
3998 ipos[2] = TMath::Min(int(r2*nclusters),nclusters-1); // last point
4001 Double_t xyz[3][3]={{0}};
4002 Int_t row[3]={0},sec[3]={0,0,0};
4004 // find track row position at given ratio of the length
4006 for (Int_t irow=0;irow<160;irow++){
4007 if (track->GetClusterIndex2(irow)<0) continue;
4009 for (Int_t ipoint=0;ipoint<3;ipoint++){
4010 if (index<=ipos[ipoint]) row[ipoint] = irow;
4014 //Get cluster and sector position
4015 for (Int_t ipoint=0;ipoint<3;ipoint++){
4016 Int_t clindex = track->GetClusterIndex2(row[ipoint]);
4017 AliTPCclusterMI * cl = GetClusterMI(clindex);
4020 // AliTPCclusterMI * cl = GetClusterMI(clindex);
4023 sec[ipoint] = ((clindex&0xff000000)>>24)%18;
4024 xyz[ipoint][0] = GetXrow(row[ipoint]);
4025 xyz[ipoint][1] = cl->GetY();
4026 xyz[ipoint][2] = cl->GetZ();
4030 // Calculate seed state vector and covariance matrix
4032 Double_t alpha, cs,sn, xx2,yy2;
4034 alpha = (sec[1]-sec[2])*fSectors->GetAlpha();
4035 cs = TMath::Cos(alpha);
4036 sn = TMath::Sin(alpha);
4037 xx2= xyz[1][0]*cs-xyz[1][1]*sn;
4038 yy2= xyz[1][0]*sn+xyz[1][1]*cs;
4042 alpha = (sec[0]-sec[2])*fSectors->GetAlpha();
4043 cs = TMath::Cos(alpha);
4044 sn = TMath::Sin(alpha);
4045 xx2= xyz[0][0]*cs-xyz[0][1]*sn;
4046 yy2= xyz[0][0]*sn+xyz[0][1]*cs;
4052 Double_t x[5],c[15];
4056 x[4]=F1(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
4057 x[2]=F2(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
4058 x[3]=F3n(xyz[2][0],xyz[2][1],xyz[0][0],xyz[0][1],xyz[2][2],xyz[0][2],x[4]);
4060 Double_t sy =0.1, sz =0.1;
4062 Double_t sy1=0.2, sz1=0.2;
4063 Double_t sy2=0.2, sz2=0.2;
4066 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;
4067 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;
4068 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;
4069 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;
4070 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;
4071 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;
4073 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;
4074 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;
4075 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;
4076 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;
4081 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
4082 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
4083 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
4084 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
4085 c[13]=f30*sy1*f40+f32*sy2*f42;
4086 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
4088 // Int_t row1 = fSectors->GetRowNumber(xyz[2][0]);
4089 AliTPCseed *seed=new( NextFreeSeed() ) AliTPCseed(xyz[2][0], sec[2]*fSectors->GetAlpha()+fSectors->GetAlphaShift(), x, c, 0);
4090 seed->SetPoolID(fLastSeedID);
4091 seed->SetLastPoint(row[2]);
4092 seed->SetFirstPoint(row[2]);
4097 AliTPCseed *AliTPCtrackerMI::ReSeed(AliTPCseed *track,Int_t r0, Bool_t forward)
4101 //reseed using founded clusters
4104 Int_t row[3]={0,0,0};
4105 Int_t sec[3]={0,0,0};
4107 // forward direction
4109 for (Int_t irow=r0;irow<160;irow++){
4110 if (track->GetClusterIndex(irow)>0){
4115 for (Int_t irow=160;irow>r0;irow--){
4116 if (track->GetClusterIndex(irow)>0){
4121 for (Int_t irow=row[2]-15;irow>row[0];irow--){
4122 if (track->GetClusterIndex(irow)>0){
4130 for (Int_t irow=0;irow<r0;irow++){
4131 if (track->GetClusterIndex(irow)>0){
4136 for (Int_t irow=r0;irow>0;irow--){
4137 if (track->GetClusterIndex(irow)>0){
4142 for (Int_t irow=row[2]-15;irow>row[0];irow--){
4143 if (track->GetClusterIndex(irow)>0){
4150 if ((row[2]-row[0])<20) return 0;
4151 if (row[1]==0) return 0;
4154 //Get cluster and sector position
4155 for (Int_t ipoint=0;ipoint<3;ipoint++){
4156 Int_t clindex = track->GetClusterIndex2(row[ipoint]);
4157 AliTPCclusterMI * cl = GetClusterMI(clindex);
4160 // AliTPCclusterMI * cl = GetClusterMI(clindex);
4163 sec[ipoint] = ((clindex&0xff000000)>>24)%18;
4164 xyz[ipoint][0] = GetXrow(row[ipoint]);
4165 AliTPCTrackerPoint * point = track->GetTrackPoint(row[ipoint]);
4166 if (point&&ipoint<2){
4168 xyz[ipoint][1] = point->GetY();
4169 xyz[ipoint][2] = point->GetZ();
4172 xyz[ipoint][1] = cl->GetY();
4173 xyz[ipoint][2] = cl->GetZ();
4180 // Calculate seed state vector and covariance matrix
4182 Double_t alpha, cs,sn, xx2,yy2;
4184 alpha = (sec[1]-sec[2])*fSectors->GetAlpha();
4185 cs = TMath::Cos(alpha);
4186 sn = TMath::Sin(alpha);
4187 xx2= xyz[1][0]*cs-xyz[1][1]*sn;
4188 yy2= xyz[1][0]*sn+xyz[1][1]*cs;
4192 alpha = (sec[0]-sec[2])*fSectors->GetAlpha();
4193 cs = TMath::Cos(alpha);
4194 sn = TMath::Sin(alpha);
4195 xx2= xyz[0][0]*cs-xyz[0][1]*sn;
4196 yy2= xyz[0][0]*sn+xyz[0][1]*cs;
4202 Double_t x[5],c[15];
4206 x[4]=F1(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
4207 x[2]=F2(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
4208 x[3]=F3n(xyz[2][0],xyz[2][1],xyz[0][0],xyz[0][1],xyz[2][2],xyz[0][2],x[4]);
4210 Double_t sy =0.1, sz =0.1;
4212 Double_t sy1=0.2, sz1=0.2;
4213 Double_t sy2=0.2, sz2=0.2;
4216 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;
4217 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;
4218 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;
4219 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;
4220 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;
4221 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;
4223 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;
4224 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;
4225 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;
4226 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;
4231 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
4232 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
4233 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
4234 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
4235 c[13]=f30*sy1*f40+f32*sy2*f42;
4236 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
4238 // Int_t row1 = fSectors->GetRowNumber(xyz[2][0]);
4239 AliTPCseed *seed=new( NextFreeSeed() ) AliTPCseed(xyz[2][0], sec[2]*fSectors->GetAlpha()+fSectors->GetAlphaShift(), x, c, 0);
4240 seed->SetPoolID(fLastSeedID);
4241 seed->SetLastPoint(row[2]);
4242 seed->SetFirstPoint(row[2]);
4243 for (Int_t i=row[0];i<row[2];i++){
4244 seed->SetClusterIndex(i, track->GetClusterIndex(i));
4252 void AliTPCtrackerMI::FindMultiMC(const TObjArray * array, AliESDEvent */*esd*/, Int_t iter)
4255 // find multi tracks - THIS FUNCTION IS ONLY FOR DEBUG PURPOSES
4257 // Use AliTPCReconstructor::StreamLevel()>2 if you want to tune parameters - cuts
4259 // Two reasons to have multiple find tracks
4260 // 1. Curling tracks can be find more than once
4261 // 2. Splitted tracks
4262 // a.) Multiple seeding to increase tracking efficiency - (~ 100% reached)
4263 // b.) Edge effect on the sector boundaries
4266 // Algorithm done in 2 phases - because of CPU consumption
4267 // it is n^2 algorithm - for lead-lead 20000x20000 combination are investigated
4269 // Algorihm for curling tracks sign:
4270 // 1 phase -makes a very rough fast cuts to minimize combinatorics
4271 // a.) opposite sign
4272 // b.) one of the tracks - not pointing to the primary vertex -
4273 // c.) delta tan(theta)
4275 // 2 phase - calculates DCA between tracks - time consument
4280 // General cuts - for splitted tracks and for curling tracks
4282 const Float_t kMaxdPhi = 0.2; // maximal distance in phi
4284 // Curling tracks cuts
4289 Int_t nentries = array->GetEntriesFast();
4290 AliHelix *helixes = new AliHelix[nentries];
4291 Float_t *xm = new Float_t[nentries];
4292 Float_t *dz0 = new Float_t[nentries];
4293 Float_t *dz1 = new Float_t[nentries];
4299 // Find track COG in x direction - point with best defined parameters
4301 for (Int_t i=0;i<nentries;i++){
4302 AliTPCseed* track = (AliTPCseed*)array->At(i);
4303 if (!track) continue;
4304 track->SetCircular(0);
4305 new (&helixes[i]) AliHelix(*track);
4309 track->GetDZ(GetX(),GetY(),GetZ(),GetBz(),dz);
4312 for (Int_t icl=0; icl<160; icl++){
4313 AliTPCclusterMI * cl = track->GetClusterPointer(icl);
4319 if (ncl>0) xm[i]/=Float_t(ncl);
4322 for (Int_t i0=0;i0<nentries;i0++){
4323 AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
4324 if (!track0) continue;
4325 Float_t xc0 = helixes[i0].GetHelix(6);
4326 Float_t yc0 = helixes[i0].GetHelix(7);
4327 Float_t r0 = helixes[i0].GetHelix(8);
4328 Float_t rc0 = TMath::Sqrt(xc0*xc0+yc0*yc0);
4329 Float_t fi0 = TMath::ATan2(yc0,xc0);
4331 for (Int_t i1=i0+1;i1<nentries;i1++){
4332 AliTPCseed * track1 = (AliTPCseed*)array->At(i1);
4333 if (!track1) continue;
4334 Int_t lab0=track0->GetLabel();
4335 Int_t lab1=track1->GetLabel();
4336 if (TMath::Abs(lab0)!=TMath::Abs(lab1)) continue;
4338 Float_t xc1 = helixes[i1].GetHelix(6);
4339 Float_t yc1 = helixes[i1].GetHelix(7);
4340 Float_t r1 = helixes[i1].GetHelix(8);
4341 Float_t rc1 = TMath::Sqrt(xc1*xc1+yc1*yc1);
4342 Float_t fi1 = TMath::ATan2(yc1,xc1);
4344 Float_t dfi = fi0-fi1;
4347 if (dfi>1.5*TMath::Pi()) dfi-=TMath::Pi(); // take care about edge effect
4348 if (dfi<-1.5*TMath::Pi()) dfi+=TMath::Pi(); //
4349 if (TMath::Abs(dfi)>kMaxdPhi&&helixes[i0].GetHelix(4)*helixes[i1].GetHelix(4)<0){
4351 // if short tracks with undefined sign
4352 fi1 = -TMath::ATan2(yc1,-xc1);
4355 Float_t dtheta = TMath::Abs(track0->GetTgl()-track1->GetTgl())<TMath::Abs(track0->GetTgl()+track1->GetTgl())? track0->GetTgl()-track1->GetTgl():track0->GetTgl()+track1->GetTgl();
4358 // debug stream to tune "fast cuts"
4360 Double_t dist[3]; // distance at X
4361 Double_t mdist[3]={0,0,0}; // mean distance X+-40cm
4362 track0->GetDistance(track1,0.5*(xm[i0]+xm[i1])-40.,dist,AliTracker::GetBz());
4363 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[i]);
4364 track0->GetDistance(track1,0.5*(xm[i0]+xm[i1])+40.,dist,AliTracker::GetBz());
4365 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[i]);
4366 track0->GetDistance(track1,0.5*(xm[i0]+xm[i1]),dist,AliTracker::GetBz());
4367 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[i]);
4368 for (Int_t i=0;i<3;i++) mdist[i]*=0.33333;
4372 for (Int_t icl=0; icl<160; icl++){
4373 AliTPCclusterMI * cl0 = track0->GetClusterPointer(icl);
4374 AliTPCclusterMI * cl1 = track1->GetClusterPointer(icl);
4377 if (cl0==cl1) sums++;
4381 if (AliTPCReconstructor::StreamLevel()>5) {
4382 TTreeSRedirector &cstream = *fDebugStreamer;
4387 "Tr0.="<<track0<< // seed0
4388 "Tr1.="<<track1<< // seed1
4389 "h0.="<<&helixes[i0]<<
4390 "h1.="<<&helixes[i1]<<
4392 "sum="<<sum<< //the sum of rows with cl in both
4393 "sums="<<sums<< //the sum of shared clusters
4394 "xm0="<<xm[i0]<< // the center of track
4395 "xm1="<<xm[i1]<< // the x center of track
4396 // General cut variables
4397 "dfi="<<dfi<< // distance in fi angle
4398 "dtheta="<<dtheta<< // distance int theta angle
4404 "dist0="<<dist[0]<< //distance x
4405 "dist1="<<dist[1]<< //distance y
4406 "dist2="<<dist[2]<< //distance z
4407 "mdist0="<<mdist[0]<< //distance x
4408 "mdist1="<<mdist[1]<< //distance y
4409 "mdist2="<<mdist[2]<< //distance z
4425 if (AliTPCReconstructor::StreamLevel()>1) {
4426 AliInfo("Time for curling tracks removal DEBUGGING MC");
4433 void AliTPCtrackerMI::FindSplitted(TObjArray * array, AliESDEvent */*esd*/, Int_t /*iter*/){
4435 // Find Splitted tracks and remove the one with worst quality
4436 // Corresponding debug streamer to tune selections - "Splitted2"
4438 // 0. Sort tracks according quility
4439 // 1. Propagate the tracks to the reference radius
4440 // 2. Double_t loop to select close tracks (only to speed up process)
4441 // 3. Calculate cluster overlap ratio - and remove the track if bigger than a threshold
4442 // 4. Delete temporary parameters
4444 const Double_t xref=GetXrow(63); // reference radius -IROC/OROC boundary
4446 const Double_t kCutP1=10; // delta Z cut 10 cm
4447 const Double_t kCutP2=0.15; // delta snp(fi) cut 0.15
4448 const Double_t kCutP3=0.15; // delta tgl(theta) cut 0.15
4449 const Double_t kCutAlpha=0.15; // delta alpha cut
4450 Int_t firstpoint = 0;
4451 Int_t lastpoint = 160;
4453 Int_t nentries = array->GetEntriesFast();
4454 AliExternalTrackParam *params = new AliExternalTrackParam[nentries];
4460 //0. Sort tracks according quality
4461 //1. Propagate the ext. param to reference radius
4462 Int_t nseed = array->GetEntriesFast();
4463 if (nseed<=0) return;
4464 Float_t * quality = new Float_t[nseed];
4465 Int_t * indexes = new Int_t[nseed];
4466 for (Int_t i=0; i<nseed; i++) {
4467 AliTPCseed *pt=(AliTPCseed*)array->UncheckedAt(i);
4472 pt->UpdatePoints(); //select first last max dens points
4473 Float_t * points = pt->GetPoints();
4474 if (points[3]<0.8) quality[i] =-1;
4475 quality[i] = (points[2]-points[0])+pt->GetNumberOfClusters();
4476 //prefer high momenta tracks if overlaps
4477 quality[i] *= TMath::Sqrt(TMath::Abs(pt->Pt())+0.5);
4479 AliTracker::PropagateTrackToBxByBz(&(params[i]),xref,pt->GetMass(),5.,kTRUE);
4480 AliTracker::PropagateTrackToBxByBz(&(params[i]),xref,pt->GetMass(),1.,kTRUE);
4482 TMath::Sort(nseed,quality,indexes);
4484 // 3. Loop over pair of tracks
4486 for (Int_t i0=0; i0<nseed; i0++) {
4487 Int_t index0=indexes[i0];
4488 if (!(array->UncheckedAt(index0))) continue;
4489 AliTPCseed *s1 = (AliTPCseed*)array->UncheckedAt(index0);
4490 if (!s1->IsActive()) continue;
4491 AliExternalTrackParam &par0=params[index0];
4492 for (Int_t i1=i0+1; i1<nseed; i1++) {
4493 Int_t index1=indexes[i1];
4494 if (!(array->UncheckedAt(index1))) continue;
4495 AliTPCseed *s2 = (AliTPCseed*)array->UncheckedAt(index1);
4496 if (!s2->IsActive()) continue;
4497 if (s2->GetKinkIndexes()[0]!=0)
4498 if (s2->GetKinkIndexes()[0] == -s1->GetKinkIndexes()[0]) continue;
4499 AliExternalTrackParam &par1=params[index1];
4500 if (TMath::Abs(par0.GetParameter()[3]-par1.GetParameter()[3])>kCutP3) continue;
4501 if (TMath::Abs(par0.GetParameter()[1]-par1.GetParameter()[1])>kCutP1) continue;
4502 if (TMath::Abs(par0.GetParameter()[2]-par1.GetParameter()[2])>kCutP2) continue;
4503 Double_t dAlpha= TMath::Abs(par0.GetAlpha()-par1.GetAlpha());
4504 if (dAlpha>TMath::Pi()) dAlpha-=TMath::Pi();
4505 if (TMath::Abs(dAlpha)>kCutAlpha) continue;
4510 Int_t firstShared=lastpoint, lastShared=firstpoint;
4511 Int_t firstRow=lastpoint, lastRow=firstpoint;
4513 for (Int_t i=firstpoint;i<lastpoint;i++){
4514 if (s1->GetClusterIndex2(i)>0) nall0++;
4515 if (s2->GetClusterIndex2(i)>0) nall1++;
4516 if (s1->GetClusterIndex2(i)>0 && s2->GetClusterIndex2(i)>0) {
4517 if (i<firstRow) firstRow=i;
4518 if (i>lastRow) lastRow=i;
4520 if ( (s1->GetClusterIndex2(i))==(s2->GetClusterIndex2(i)) && s1->GetClusterIndex2(i)>0) {
4521 if (i<firstShared) firstShared=i;
4522 if (i>lastShared) lastShared=i;
4526 Double_t ratio0 = Float_t(sumShared)/Float_t(TMath::Min(nall0+1,nall1+1));
4527 Double_t ratio1 = Float_t(sumShared)/Float_t(TMath::Max(nall0+1,nall1+1));
4529 if( AliTPCReconstructor::StreamLevel()>1){
4530 TTreeSRedirector &cstream = *fDebugStreamer;
4531 Int_t n0=s1->GetNumberOfClusters();
4532 Int_t n1=s2->GetNumberOfClusters();
4533 Int_t n0F=s1->GetNFoundable();
4534 Int_t n1F=s2->GetNFoundable();
4535 Int_t lab0=s1->GetLabel();
4536 Int_t lab1=s2->GetLabel();
4538 cstream<<"Splitted2"<<
4539 "iter="<<fIteration<<
4540 "lab0="<<lab0<< // MC label if exist
4541 "lab1="<<lab1<< // MC label if exist
4544 "ratio0="<<ratio0<< // shared ratio
4545 "ratio1="<<ratio1<< // shared ratio
4546 "p0.="<<&par0<< // track parameters
4548 "s0.="<<s1<< // full seed
4550 "n0="<<n0<< // number of clusters track 0
4551 "n1="<<n1<< // number of clusters track 1
4552 "nall0="<<nall0<< // number of clusters track 0
4553 "nall1="<<nall1<< // number of clusters track 1
4554 "n0F="<<n0F<< // number of findable
4555 "n1F="<<n1F<< // number of findable
4556 "shared="<<sumShared<< // number of shared clusters
4557 "firstS="<<firstShared<< // first and the last shared row
4558 "lastS="<<lastShared<<
4559 "firstRow="<<firstRow<< // first and the last row with cluster
4560 "lastRow="<<lastRow<< //
4564 // remove track with lower quality
4566 if (ratio0>AliTPCReconstructor::GetRecoParam()->GetCutSharedClusters(0) ||
4567 ratio1>AliTPCReconstructor::GetRecoParam()->GetCutSharedClusters(1)){
4571 MarkSeedFree( array->RemoveAt(index1) );
4576 // 4. Delete temporary array
4586 void AliTPCtrackerMI::FindCurling(const TObjArray * array, AliESDEvent */*esd*/, Int_t iter)
4589 // find Curling tracks
4590 // Use AliTPCReconstructor::StreamLevel()>1 if you want to tune parameters - cuts
4593 // Algorithm done in 2 phases - because of CPU consumption
4594 // it is n^2 algorithm - for lead-lead 20000x20000 combination are investigated
4595 // see detal in MC part what can be used to cut
4599 const Float_t kMaxC = 400; // maximal curvature to of the track
4600 const Float_t kMaxdTheta = 0.15; // maximal distance in theta
4601 const Float_t kMaxdPhi = 0.15; // maximal distance in phi
4602 const Float_t kPtRatio = 0.3; // ratio between pt
4603 const Float_t kMinDCAR = 2.; // distance to the primary vertex in r - see cpipe cut
4606 // Curling tracks cuts
4609 const Float_t kMaxDeltaRMax = 40; // distance in outer radius
4610 const Float_t kMaxDeltaRMin = 5.; // distance in lower radius - see cpipe cut
4611 const Float_t kMinAngle = 2.9; // angle between tracks
4612 const Float_t kMaxDist = 5; // biggest distance
4614 // The cuts can be tuned using the "MC information stored in Multi tree ==> see FindMultiMC
4617 TCut csign("csign","Tr0.fP[4]*Tr1.fP[4]<0"); //opposite sign
4618 TCut cmax("cmax","abs(Tr0.GetC())>1/400");
4619 TCut cda("cda","sqrt(dtheta^2+dfi^2)<0.15");
4620 TCut ccratio("ccratio","abs((Tr0.fP[4]+Tr1.fP[4])/(abs(Tr0.fP[4])+abs(Tr1.fP[4])))<0.3");
4621 TCut cpipe("cpipe", "min(abs(r0-rc0),abs(r1-rc1))>5");
4623 TCut cdrmax("cdrmax","abs(abs(rc0+r0)-abs(rc1+r1))<40")
4624 TCut cdrmin("cdrmin","abs(abs(rc0+r0)-abs(rc1+r1))<10")
4626 Multi->Draw("dfi","iter==0"+csign+cmax+cda+ccratio); ~94% of curling tracks fulfill
4627 Multi->Draw("min(abs(r0-rc0),abs(r1-rc1))","iter==0&&abs(lab1)==abs(lab0)"+csign+cmax+cda+ccratio+cpipe+cdrmin+cdrmax); //80%
4629 Curling2->Draw("dfi","iter==0&&abs(lab0)==abs(lab1)"+csign+cmax+cdtheta+cdfi+ccratio)
4635 Int_t nentries = array->GetEntriesFast();
4636 AliHelix *helixes = new AliHelix[nentries];
4637 for (Int_t i=0;i<nentries;i++){
4638 AliTPCseed* track = (AliTPCseed*)array->At(i);
4639 if (!track) continue;
4640 track->SetCircular(0);
4641 new (&helixes[i]) AliHelix(*track);
4647 Double_t phase[2][2]={{0,0},{0,0}},radius[2]={0,0};
4653 for (Int_t i0=0;i0<nentries;i0++){
4654 AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
4655 if (!track0) continue;
4656 if (TMath::Abs(track0->GetC())<1/kMaxC) continue;
4657 Float_t xc0 = helixes[i0].GetHelix(6);
4658 Float_t yc0 = helixes[i0].GetHelix(7);
4659 Float_t r0 = helixes[i0].GetHelix(8);
4660 Float_t rc0 = TMath::Sqrt(xc0*xc0+yc0*yc0);
4661 Float_t fi0 = TMath::ATan2(yc0,xc0);
4663 for (Int_t i1=i0+1;i1<nentries;i1++){
4664 AliTPCseed * track1 = (AliTPCseed*)array->At(i1);
4665 if (!track1) continue;
4666 if (TMath::Abs(track1->GetC())<1/kMaxC) continue;
4667 Float_t xc1 = helixes[i1].GetHelix(6);
4668 Float_t yc1 = helixes[i1].GetHelix(7);
4669 Float_t r1 = helixes[i1].GetHelix(8);
4670 Float_t rc1 = TMath::Sqrt(xc1*xc1+yc1*yc1);
4671 Float_t fi1 = TMath::ATan2(yc1,xc1);
4673 Float_t dfi = fi0-fi1;
4676 if (dfi>1.5*TMath::Pi()) dfi-=TMath::Pi(); // take care about edge effect
4677 if (dfi<-1.5*TMath::Pi()) dfi+=TMath::Pi(); //
4678 Float_t dtheta = TMath::Abs(track0->GetTgl()-track1->GetTgl())<TMath::Abs(track0->GetTgl()+track1->GetTgl())? track0->GetTgl()-track1->GetTgl():track0->GetTgl()+track1->GetTgl();
4682 if (track0->GetBConstrain()&&track1->GetBConstrain()) continue; // not constrained
4683 if (track1->GetSigned1Pt()*track0->GetSigned1Pt()>0) continue; // not the same sign
4684 if ( TMath::Abs(track1->GetTgl()+track0->GetTgl())>kMaxdTheta) continue; //distance in the Theta
4685 if ( TMath::Abs(dfi)>kMaxdPhi) continue; //distance in phi
4686 if ( TMath::Sqrt(dfi*dfi+dtheta*dtheta)>kMaxdPhi) continue; //common angular offset
4688 Float_t pt0 = track0->GetSignedPt();
4689 Float_t pt1 = track1->GetSignedPt();
4690 if ((TMath::Abs(pt0+pt1)/(TMath::Abs(pt0)+TMath::Abs(pt1)))>kPtRatio) continue;
4691 if ((iter==1) && TMath::Abs(TMath::Abs(rc0+r0)-TMath::Abs(rc1+r1))>kMaxDeltaRMax) continue;
4692 if ((iter!=1) &&TMath::Abs(TMath::Abs(rc0-r0)-TMath::Abs(rc1-r1))>kMaxDeltaRMin) continue;
4693 if (TMath::Min(TMath::Abs(rc0-r0),TMath::Abs(rc1-r1))<kMinDCAR) continue;
4696 // Now find closest approach
4700 Int_t npoints = helixes[i0].GetRPHIintersections(helixes[i1], phase, radius,10);
4701 if (npoints==0) continue;
4702 helixes[i0].GetClosestPhases(helixes[i1], phase);
4706 Double_t hangles[3];
4707 helixes[i0].Evaluate(phase[0][0],xyz0);
4708 helixes[i1].Evaluate(phase[0][1],xyz1);
4710 helixes[i0].GetAngle(phase[0][0],helixes[i1],phase[0][1],hangles);
4711 Double_t deltah[2],deltabest;
4712 if (TMath::Abs(hangles[2])<kMinAngle) continue;
4716 helixes[i0].ParabolicDCA(helixes[i1],phase[0][0],phase[0][1],radius[0],deltah[0],2);
4718 helixes[i0].ParabolicDCA(helixes[i1],phase[1][0],phase[1][1],radius[1],deltah[1],2);
4719 if (deltah[1]<deltah[0]) ibest=1;
4721 deltabest = TMath::Sqrt(deltah[ibest]);
4722 helixes[i0].Evaluate(phase[ibest][0],xyz0);
4723 helixes[i1].Evaluate(phase[ibest][1],xyz1);
4724 helixes[i0].GetAngle(phase[ibest][0],helixes[i1],phase[ibest][1],hangles);
4725 Double_t radiusbest = TMath::Sqrt(radius[ibest]);
4727 if (deltabest>kMaxDist) continue;
4728 // if (mindcar+mindcaz<40 && (TMath::Abs(hangles[2])<kMinAngle ||deltabest>3)) continue;
4729 Bool_t sign =kFALSE;
4730 if (hangles[2]>kMinAngle) sign =kTRUE;
4733 // circular[i0] = kTRUE;
4734 // circular[i1] = kTRUE;
4735 if (track0->OneOverPt()<track1->OneOverPt()){
4736 track0->SetCircular(track0->GetCircular()+1);
4737 track1->SetCircular(track1->GetCircular()+2);
4740 track1->SetCircular(track1->GetCircular()+1);
4741 track0->SetCircular(track0->GetCircular()+2);
4744 if (AliTPCReconstructor::StreamLevel()>2){
4746 //debug stream to tune "fine" cuts
4747 Int_t lab0=track0->GetLabel();
4748 Int_t lab1=track1->GetLabel();
4749 TTreeSRedirector &cstream = *fDebugStreamer;
4750 cstream<<"Curling2"<<
4766 "npoints="<<npoints<<
4767 "hangles0="<<hangles[0]<<
4768 "hangles1="<<hangles[1]<<
4769 "hangles2="<<hangles[2]<<
4772 "radius="<<radiusbest<<
4773 "deltabest="<<deltabest<<
4774 "phase0="<<phase[ibest][0]<<
4775 "phase1="<<phase[ibest][1]<<
4783 if (AliTPCReconstructor::StreamLevel()>1) {
4784 AliInfo("Time for curling tracks removal");
4790 void AliTPCtrackerMI::FindKinks(TObjArray * array, AliESDEvent *esd)
4796 // RS something is wrong in this routine: not all seeds are assigned to daughters and mothers array, but they all are queried
4799 TObjArray *kinks= new TObjArray(10000);
4800 // TObjArray *v0s= new TObjArray(10000);
4801 Int_t nentries = array->GetEntriesFast();
4802 AliHelix *helixes = new AliHelix[nentries];
4803 Int_t *sign = new Int_t[nentries];
4804 Int_t *nclusters = new Int_t[nentries];
4805 Float_t *alpha = new Float_t[nentries];
4806 AliKink *kink = new AliKink();
4807 Int_t * usage = new Int_t[nentries];
4808 Float_t *zm = new Float_t[nentries];
4809 Float_t *z0 = new Float_t[nentries];
4810 Float_t *fim = new Float_t[nentries];
4811 Float_t *shared = new Float_t[nentries];
4812 Bool_t *circular = new Bool_t[nentries];
4813 Float_t *dca = new Float_t[nentries];
4814 //const AliESDVertex * primvertex = esd->GetVertex();
4816 // nentries = array->GetEntriesFast();
4821 for (Int_t i=0;i<nentries;i++){
4824 AliTPCseed* track = (AliTPCseed*)array->At(i);
4825 if (!track) continue;
4826 track->SetCircular(0);
4828 track->UpdatePoints();
4829 if (( track->GetPoints()[2]- track->GetPoints()[0])>5 && track->GetPoints()[3]>0.8){
4831 nclusters[i]=track->GetNumberOfClusters();
4832 alpha[i] = track->GetAlpha();
4833 new (&helixes[i]) AliHelix(*track);
4835 helixes[i].Evaluate(0,xyz);
4836 sign[i] = (track->GetC()>0) ? -1:1;
4839 if (track->GetProlongation(x,y,z)){
4841 fim[i] = alpha[i]+TMath::ATan2(y,x);
4844 zm[i] = track->GetZ();
4848 circular[i]= kFALSE;
4849 if (track->GetProlongation(0,y,z)) z0[i] = z;
4850 dca[i] = track->GetD(0,0);
4856 Int_t ncandidates =0;
4859 Double_t phase[2][2]={{0,0},{0,0}},radius[2]={0,0};
4862 // Find circling track
4864 for (Int_t i0=0;i0<nentries;i0++){
4865 AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
4866 if (!track0) continue;
4867 if (track0->GetNumberOfClusters()<40) continue;
4868 if (TMath::Abs(1./track0->GetC())>200) continue;
4869 for (Int_t i1=i0+1;i1<nentries;i1++){
4870 AliTPCseed * track1 = (AliTPCseed*)array->At(i1);
4871 if (!track1) continue;
4872 if (track1->GetNumberOfClusters()<40) continue;
4873 if ( TMath::Abs(track1->GetTgl()+track0->GetTgl())>0.1) continue;
4874 if (track0->GetBConstrain()&&track1->GetBConstrain()) continue;
4875 if (TMath::Abs(1./track1->GetC())>200) continue;
4876 if (track1->GetSigned1Pt()*track0->GetSigned1Pt()>0) continue;
4877 if (track1->GetTgl()*track0->GetTgl()>0) continue;
4878 if (TMath::Max(TMath::Abs(1./track0->GetC()),TMath::Abs(1./track1->GetC()))>190) continue;
4879 if (track0->GetBConstrain()&&track1->OneOverPt()<track0->OneOverPt()) continue; //returning - lower momenta
4880 if (track1->GetBConstrain()&&track0->OneOverPt()<track1->OneOverPt()) continue; //returning - lower momenta
4882 Float_t mindcar = TMath::Min(TMath::Abs(dca[i0]),TMath::Abs(dca[i1]));
4883 if (mindcar<5) continue;
4884 Float_t mindcaz = TMath::Min(TMath::Abs(z0[i0]-GetZ()),TMath::Abs(z0[i1]-GetZ()));
4885 if (mindcaz<5) continue;
4886 if (mindcar+mindcaz<20) continue;
4889 Float_t xc0 = helixes[i0].GetHelix(6);
4890 Float_t yc0 = helixes[i0].GetHelix(7);
4891 Float_t r0 = helixes[i0].GetHelix(8);
4892 Float_t xc1 = helixes[i1].GetHelix(6);
4893 Float_t yc1 = helixes[i1].GetHelix(7);
4894 Float_t r1 = helixes[i1].GetHelix(8);
4896 Float_t rmean = (r0+r1)*0.5;
4897 Float_t delta =TMath::Sqrt((xc1-xc0)*(xc1-xc0)+(yc1-yc0)*(yc1-yc0));
4898 //if (delta>30) continue;
4899 if (delta>rmean*0.25) continue;
4900 if (TMath::Abs(r0-r1)/rmean>0.3) continue;
4902 Int_t npoints = helixes[i0].GetRPHIintersections(helixes[i1], phase, radius,10);
4903 if (npoints==0) continue;
4904 helixes[i0].GetClosestPhases(helixes[i1], phase);
4908 Double_t hangles[3];
4909 helixes[i0].Evaluate(phase[0][0],xyz0);
4910 helixes[i1].Evaluate(phase[0][1],xyz1);
4912 helixes[i0].GetAngle(phase[0][0],helixes[i1],phase[0][1],hangles);
4913 Double_t deltah[2],deltabest;
4914 if (hangles[2]<2.8) continue;
4917 helixes[i0].ParabolicDCA(helixes[i1],phase[0][0],phase[0][1],radius[0],deltah[0],2);
4919 helixes[i0].ParabolicDCA(helixes[i1],phase[1][0],phase[1][1],radius[1],deltah[1],2);
4920 if (deltah[1]<deltah[0]) ibest=1;
4922 deltabest = TMath::Sqrt(deltah[ibest]);
4923 helixes[i0].Evaluate(phase[ibest][0],xyz0);
4924 helixes[i1].Evaluate(phase[ibest][1],xyz1);
4925 helixes[i0].GetAngle(phase[ibest][0],helixes[i1],phase[ibest][1],hangles);
4926 Double_t radiusbest = TMath::Sqrt(radius[ibest]);
4928 if (deltabest>6) continue;
4929 if (mindcar+mindcaz<40 && (hangles[2]<3.12||deltabest>3)) continue;
4930 Bool_t lsign =kFALSE;
4931 if (hangles[2]>3.06) lsign =kTRUE;
4934 circular[i0] = kTRUE;
4935 circular[i1] = kTRUE;
4936 if (track0->OneOverPt()<track1->OneOverPt()){
4937 track0->SetCircular(track0->GetCircular()+1);
4938 track1->SetCircular(track1->GetCircular()+2);
4941 track1->SetCircular(track1->GetCircular()+1);
4942 track0->SetCircular(track0->GetCircular()+2);
4945 if (lsign&&AliTPCReconstructor::StreamLevel()>1){
4947 Int_t lab0=track0->GetLabel();
4948 Int_t lab1=track1->GetLabel();
4949 TTreeSRedirector &cstream = *fDebugStreamer;
4950 cstream<<"Curling"<<
4957 "mindcar="<<mindcar<<
4958 "mindcaz="<<mindcaz<<
4961 "npoints="<<npoints<<
4962 "hangles0="<<hangles[0]<<
4963 "hangles2="<<hangles[2]<<
4968 "radius="<<radiusbest<<
4969 "deltabest="<<deltabest<<
4970 "phase0="<<phase[ibest][0]<<
4971 "phase1="<<phase[ibest][1]<<
4981 for (Int_t i =0;i<nentries;i++){
4982 if (sign[i]==0) continue;
4983 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
4990 Double_t cradius0 = 40*40;
4991 Double_t cradius1 = 270*270;
4994 Double_t cdist3=0.55;
4995 for (Int_t j =i+1;j<nentries;j++){
4997 if (sign[j]*sign[i]<1) continue;
4998 if ( (nclusters[i]+nclusters[j])>200) continue;
4999 if ( (nclusters[i]+nclusters[j])<80) continue;
5000 if ( TMath::Abs(zm[i]-zm[j])>60.) continue;
5001 if ( TMath::Abs(fim[i]-fim[j])>0.6 && TMath::Abs(fim[i]-fim[j])<5.7 ) continue;
5002 //AliTPCseed * track1 = (AliTPCseed*)array->At(j); Double_t phase[2][2],radius[2];
5003 Int_t npoints = helixes[i].GetRPHIintersections(helixes[j], phase, radius,20);
5004 if (npoints<1) continue;
5007 if (radius[0]<cradius0||radius[0]>cradius1) continue;
5010 if ( (radius[0]<cradius0||radius[0]>cradius1) && (radius[1]<cradius0||radius[1]>cradius1) ) continue;
5013 Double_t delta1=10000,delta2=10000;
5014 // cuts on the intersection radius
5015 helixes[i].LinearDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta1);
5016 if (radius[0]<20&&delta1<1) continue; //intersection at vertex
5017 if (radius[0]<10&&delta1<3) continue; //intersection at vertex
5019 helixes[i].LinearDCA(helixes[j],phase[1][0],phase[1][1],radius[1],delta2);
5020 if (radius[1]<20&&delta2<1) continue; //intersection at vertex
5021 if (radius[1]<10&&delta2<3) continue; //intersection at vertex
5024 Double_t distance1 = TMath::Min(delta1,delta2);
5025 if (distance1>cdist1) continue; // cut on DCA linear approximation
5027 npoints = helixes[i].GetRPHIintersections(helixes[j], phase, radius,20);
5028 helixes[i].ParabolicDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta1);
5029 if (radius[0]<20&&delta1<1) continue; //intersection at vertex
5030 if (radius[0]<10&&delta1<3) continue; //intersection at vertex
5033 helixes[i].ParabolicDCA(helixes[j],phase[1][0],phase[1][1],radius[1],delta2);
5034 if (radius[1]<20&&delta2<1) continue; //intersection at vertex
5035 if (radius[1]<10&&delta2<3) continue; //intersection at vertex
5037 distance1 = TMath::Min(delta1,delta2);
5040 rkink = TMath::Sqrt(radius[0]);
5043 rkink = TMath::Sqrt(radius[1]);
5045 if (distance1>cdist2) continue;
5048 AliTPCseed * track1 = (AliTPCseed*)array->At(j);
5051 Int_t row0 = GetRowNumber(rkink);
5052 if (row0<10) continue;
5053 if (row0>150) continue;
5056 Float_t dens00=-1,dens01=-1;
5057 Float_t dens10=-1,dens11=-1;
5059 Int_t found,foundable,ishared;
5060 track0->GetClusterStatistic(0,row0-5, found, foundable,ishared,kFALSE);
5061 if (foundable>5) dens00 = Float_t(found)/Float_t(foundable);
5062 track0->GetClusterStatistic(row0+5,155, found, foundable,ishared,kFALSE);
5063 if (foundable>5) dens01 = Float_t(found)/Float_t(foundable);
5065 track1->GetClusterStatistic(0,row0-5, found, foundable,ishared,kFALSE);
5066 if (foundable>10) dens10 = Float_t(found)/Float_t(foundable);
5067 track1->GetClusterStatistic(row0+5,155, found, foundable,ishared,kFALSE);
5068 if (foundable>10) dens11 = Float_t(found)/Float_t(foundable);
5070 if (dens00<dens10 && dens01<dens11) continue;
5071 if (dens00>dens10 && dens01>dens11) continue;
5072 if (TMath::Max(dens00,dens10)<0.1) continue;
5073 if (TMath::Max(dens01,dens11)<0.3) continue;
5075 if (TMath::Min(dens00,dens10)>0.6) continue;
5076 if (TMath::Min(dens01,dens11)>0.6) continue;
5079 AliTPCseed * ktrack0, *ktrack1;
5088 if (TMath::Abs(ktrack0->GetC())>5) continue; // cut on the curvature for mother particle
5089 AliExternalTrackParam paramm(*ktrack0);
5090 AliExternalTrackParam paramd(*ktrack1);
5091 if (row0>60&&ktrack1->GetReference().GetX()>90.)new (¶md) AliExternalTrackParam(ktrack1->GetReference());
5094 kink->SetMother(paramm);
5095 kink->SetDaughter(paramd);
5098 Float_t x[3] = { kink->GetPosition()[0],kink->GetPosition()[1],kink->GetPosition()[2]};
5100 fkParam->Transform0to1(x,index);
5101 fkParam->Transform1to2(x,index);
5102 row0 = GetRowNumber(x[0]);
5104 if (kink->GetR()<100) continue;
5105 if (kink->GetR()>240) continue;
5106 if (kink->GetPosition()[2]/kink->GetR()>AliTPCReconstructor::GetCtgRange()) continue; //out of fiducial volume
5107 if (kink->GetDistance()>cdist3) continue;
5108 Float_t dird = kink->GetDaughterP()[0]*kink->GetPosition()[0]+kink->GetDaughterP()[1]*kink->GetPosition()[1]; // rough direction estimate
5109 if (dird<0) continue;
5111 Float_t dirm = kink->GetMotherP()[0]*kink->GetPosition()[0]+kink->GetMotherP()[1]*kink->GetPosition()[1]; // rough direction estimate
5112 if (dirm<0) continue;
5113 Float_t mpt = TMath::Sqrt(kink->GetMotherP()[0]*kink->GetMotherP()[0]+kink->GetMotherP()[1]*kink->GetMotherP()[1]);
5114 if (mpt<0.2) continue;
5117 //for high momenta momentum not defined well in first iteration
5118 Double_t qt = TMath::Sin(kink->GetAngle(2))*ktrack1->GetP();
5119 if (qt>0.35) continue;
5122 kink->SetLabel(CookLabel(ktrack0,0.4,0,row0),0);
5123 kink->SetLabel(CookLabel(ktrack1,0.4,row0,160),1);
5125 kink->SetTPCDensity(dens00,0,0);
5126 kink->SetTPCDensity(dens01,0,1);
5127 kink->SetTPCDensity(dens10,1,0);
5128 kink->SetTPCDensity(dens11,1,1);
5129 kink->SetIndex(i,0);
5130 kink->SetIndex(j,1);
5133 kink->SetTPCDensity(dens10,0,0);
5134 kink->SetTPCDensity(dens11,0,1);
5135 kink->SetTPCDensity(dens00,1,0);
5136 kink->SetTPCDensity(dens01,1,1);
5137 kink->SetIndex(j,0);
5138 kink->SetIndex(i,1);
5141 if (mpt<1||kink->GetAngle(2)>0.1){
5142 // angle and densities not defined yet
5143 if (kink->GetTPCDensityFactor()<0.8) continue;
5144 if ((2-kink->GetTPCDensityFactor())*kink->GetDistance() >0.25) continue;
5145 if (kink->GetAngle(2)*ktrack0->GetP()<0.003) continue; //too small angle
5146 if (kink->GetAngle(2)>0.2&&kink->GetTPCDensityFactor()<1.15) continue;
5147 if (kink->GetAngle(2)>0.2&&kink->GetTPCDensity(0,1)>0.05) continue;
5149 Float_t criticalangle = track0->GetSigmaSnp2()+track0->GetSigmaTgl2();
5150 criticalangle+= track1->GetSigmaSnp2()+track1->GetSigmaTgl2();
5151 criticalangle= 3*TMath::Sqrt(criticalangle);
5152 if (criticalangle>0.02) criticalangle=0.02;
5153 if (kink->GetAngle(2)<criticalangle) continue;
5156 Int_t drow = Int_t(2.+0.5/(0.05+kink->GetAngle(2))); // overlap region defined
5157 Float_t shapesum =0;
5159 for ( Int_t row = row0-drow; row<row0+drow;row++){
5160 if (row<0) continue;
5161 if (row>155) continue;
5162 if (ktrack0->GetClusterPointer(row)){
5163 AliTPCTrackerPoint *point =ktrack0->GetTrackPoint(row);
5164 shapesum+=point->GetSigmaY()+point->GetSigmaZ();
5167 if (ktrack1->GetClusterPointer(row)){
5168 AliTPCTrackerPoint *point =ktrack1->GetTrackPoint(row);
5169 shapesum+=point->GetSigmaY()+point->GetSigmaZ();
5174 kink->SetShapeFactor(-1.);
5177 kink->SetShapeFactor(shapesum/sum);
5179 // esd->AddKink(kink);
5181 // kink->SetMother(paramm);
5182 //kink->SetDaughter(paramd);
5184 Double_t chi2P2 = paramm.GetParameter()[2]-paramd.GetParameter()[2];
5186 chi2P2/=paramm.GetCovariance()[5]+paramd.GetCovariance()[5];
5187 Double_t chi2P3 = paramm.GetParameter()[3]-paramd.GetParameter()[3];
5189 chi2P3/=paramm.GetCovariance()[9]+paramd.GetCovariance()[9];
5191 if (AliTPCReconstructor::StreamLevel()>1) {
5192 (*fDebugStreamer)<<"kinkLpt"<<
5200 if ( chi2P2+chi2P3<AliTPCReconstructor::GetRecoParam()->GetKinkAngleCutChi2(0)){
5204 kinks->AddLast(kink);
5210 // sort the kinks according quality - and refit them towards vertex
5212 Int_t nkinks = kinks->GetEntriesFast();
5213 Float_t *quality = new Float_t[nkinks];
5214 Int_t *indexes = new Int_t[nkinks];
5215 AliTPCseed *mothers = new AliTPCseed[nkinks];
5216 AliTPCseed *daughters = new AliTPCseed[nkinks];
5219 for (Int_t i=0;i<nkinks;i++){
5221 AliKink *kinkl = (AliKink*)kinks->At(i);
5223 // refit kinks towards vertex
5225 Int_t index0 = kinkl->GetIndex(0);
5226 Int_t index1 = kinkl->GetIndex(1);
5227 AliTPCseed * ktrack0 = (AliTPCseed*)array->At(index0);
5228 AliTPCseed * ktrack1 = (AliTPCseed*)array->At(index1);
5230 Int_t sumn=ktrack0->GetNumberOfClusters()+ktrack1->GetNumberOfClusters();
5232 // Refit Kink under if too small angle
5234 if (kinkl->GetAngle(2)<0.05){
5235 kinkl->SetTPCRow0(GetRowNumber(kinkl->GetR()));
5236 Int_t row0 = kinkl->GetTPCRow0();
5237 Int_t drow = Int_t(2.+0.5/(0.05+kinkl->GetAngle(2)));
5240 Int_t last = row0-drow;
5241 if (last<40) last=40;
5242 if (last<ktrack0->GetFirstPoint()+25) last = ktrack0->GetFirstPoint()+25;
5243 AliTPCseed* seed0 = ReSeed(ktrack0,last,kFALSE);
5246 Int_t first = row0+drow;
5247 if (first>130) first=130;
5248 if (first>ktrack1->GetLastPoint()-25) first = TMath::Max(ktrack1->GetLastPoint()-25,30);
5249 AliTPCseed* seed1 = ReSeed(ktrack1,first,kTRUE);
5251 if (seed0 && seed1){
5252 kinkl->SetStatus(1,8);
5253 if (RefitKink(*seed0,*seed1,*kinkl)) kinkl->SetStatus(1,9);
5254 row0 = GetRowNumber(kinkl->GetR());
5255 sumn = seed0->GetNumberOfClusters()+seed1->GetNumberOfClusters();
5256 mothers[i] = *seed0;
5257 daughters[i] = *seed1;
5260 delete kinks->RemoveAt(i);
5261 if (seed0) MarkSeedFree( seed0 );
5262 if (seed1) MarkSeedFree( seed1 );
5265 if (kinkl->GetDistance()>0.5 || kinkl->GetR()<110 || kinkl->GetR()>240) {
5266 delete kinks->RemoveAt(i);
5267 if (seed0) MarkSeedFree( seed0 );
5268 if (seed1) MarkSeedFree( seed1 );
5272 MarkSeedFree( seed0 );
5273 MarkSeedFree( seed1 );
5276 if (kinkl) quality[i] = 160*((0.1+kinkl->GetDistance())*(2.-kinkl->GetTPCDensityFactor()))/(sumn+40.); //the longest -clossest will win
5278 TMath::Sort(nkinks,quality,indexes,kFALSE);
5280 //remove double find kinks
5282 for (Int_t ikink0=1;ikink0<nkinks;ikink0++){
5283 AliKink * kink0 = (AliKink*) kinks->At(indexes[ikink0]);
5284 if (!kink0) continue;
5286 for (Int_t ikink1=0;ikink1<ikink0;ikink1++){
5287 kink0 = (AliKink*) kinks->At(indexes[ikink0]);
5288 if (!kink0) continue;
5289 AliKink * kink1 = (AliKink*) kinks->At(indexes[ikink1]);
5290 if (!kink1) continue;
5291 // if not close kink continue
5292 if (TMath::Abs(kink1->GetPosition()[2]-kink0->GetPosition()[2])>10) continue;
5293 if (TMath::Abs(kink1->GetPosition()[1]-kink0->GetPosition()[1])>10) continue;
5294 if (TMath::Abs(kink1->GetPosition()[0]-kink0->GetPosition()[0])>10) continue;
5296 AliTPCseed &mother0 = mothers[indexes[ikink0]];
5297 AliTPCseed &daughter0 = daughters[indexes[ikink0]];
5298 AliTPCseed &mother1 = mothers[indexes[ikink1]];
5299 AliTPCseed &daughter1 = daughters[indexes[ikink1]];
5300 Int_t row0 = (kink0->GetTPCRow0()+kink1->GetTPCRow0())/2;
5309 for (Int_t i=0;i<row0;i++){
5310 if (mother0.GetClusterIndex(i)>0 && mother1.GetClusterIndex(i)>0){
5313 if (mother0.GetClusterIndex(i)==mother1.GetClusterIndex(i)){
5320 for (Int_t i=row0;i<158;i++){
5321 //if (daughter0.GetClusterIndex(i)>0 && daughter0.GetClusterIndex(i)>0){ // RS: Bug?
5322 if (daughter0.GetClusterIndex(i)>0 && daughter1.GetClusterIndex(i)>0){
5325 if (mother0.GetClusterIndex(i)==mother1.GetClusterIndex(i)){
5331 Float_t ratio = Float_t(same+1)/Float_t(both+1);
5332 Float_t ratiom = Float_t(samem+1)/Float_t(bothm+1);
5333 Float_t ratiod = Float_t(samed+1)/Float_t(bothd+1);
5334 if (ratio>0.3 && ratiom>0.5 &&ratiod>0.5) {
5335 Int_t sum0 = mother0.GetNumberOfClusters()+daughter0.GetNumberOfClusters();
5336 Int_t sum1 = mother1.GetNumberOfClusters()+daughter1.GetNumberOfClusters();
5338 shared[kink0->GetIndex(0)]= kTRUE;
5339 shared[kink0->GetIndex(1)]= kTRUE;
5340 delete kinks->RemoveAt(indexes[ikink0]);
5344 shared[kink1->GetIndex(0)]= kTRUE;
5345 shared[kink1->GetIndex(1)]= kTRUE;
5346 delete kinks->RemoveAt(indexes[ikink1]);
5353 for (Int_t i=0;i<nkinks;i++){
5354 AliKink * kinkl = (AliKink*) kinks->At(indexes[i]);
5355 if (!kinkl) continue;
5356 kinkl->SetTPCRow0(GetRowNumber(kinkl->GetR()));
5357 Int_t index0 = kinkl->GetIndex(0);
5358 Int_t index1 = kinkl->GetIndex(1);
5359 if (circular[index0]||(circular[index1]&&kinkl->GetDistance()>0.2)) continue;
5360 kinkl->SetMultiple(usage[index0],0);
5361 kinkl->SetMultiple(usage[index1],1);
5362 if (kinkl->GetMultiple()[0]+kinkl->GetMultiple()[1]>2) continue;
5363 if (kinkl->GetMultiple()[0]+kinkl->GetMultiple()[1]>0 && quality[indexes[i]]>0.2) continue;
5364 if (kinkl->GetMultiple()[0]+kinkl->GetMultiple()[1]>0 && kinkl->GetDistance()>0.2) continue;
5365 if (circular[index0]||(circular[index1]&&kinkl->GetDistance()>0.1)) continue;
5367 AliTPCseed * ktrack0 = (AliTPCseed*)array->At(index0);
5368 AliTPCseed * ktrack1 = (AliTPCseed*)array->At(index1);
5369 if (!ktrack0 || !ktrack1) continue;
5370 Int_t index = esd->AddKink(kinkl);
5373 if ( ktrack0->GetKinkIndex(0)==0 && ktrack1->GetKinkIndex(0)==0) { //best kink
5374 if (mothers[indexes[i]].GetNumberOfClusters()>20 && daughters[indexes[i]].GetNumberOfClusters()>20 && (mothers[indexes[i]].GetNumberOfClusters()+daughters[indexes[i]].GetNumberOfClusters())>100){
5375 *ktrack0 = mothers[indexes[i]];
5376 *ktrack1 = daughters[indexes[i]];
5380 ktrack0->SetKinkIndex(usage[index0],-(index+1));
5381 ktrack1->SetKinkIndex(usage[index1], (index+1));
5386 // Remove tracks corresponding to shared kink's
5388 for (Int_t i=0;i<nentries;i++){
5389 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5390 if (!track0) continue;
5391 if (track0->GetKinkIndex(0)!=0) continue;
5392 if (shared[i]) MarkSeedFree( array->RemoveAt(i) );
5397 RemoveUsed2(array,0.5,0.4,30);
5399 for (Int_t i=0;i<nentries;i++){
5400 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5401 if (!track0) continue;
5402 track0->CookdEdx(0.02,0.6);
5406 for (Int_t i=0;i<nentries;i++){
5407 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5408 if (!track0) continue;
5409 if (track0->Pt()<1.4) continue;
5410 //remove double high momenta tracks - overlapped with kink candidates
5413 for (Int_t icl=track0->GetFirstPoint();icl<track0->GetLastPoint(); icl++){
5414 if (track0->GetClusterPointer(icl)!=0){
5416 if (track0->GetClusterPointer(icl)->IsUsed(10)) ishared++;
5419 if (Float_t(ishared+1)/Float_t(all+1)>0.5) {
5420 MarkSeedFree( array->RemoveAt(i) );
5424 if (track0->GetKinkIndex(0)!=0) continue;
5425 if (track0->GetNumberOfClusters()<80) continue;
5427 AliTPCseed *pmother = new AliTPCseed();
5428 AliTPCseed *pdaughter = new AliTPCseed();
5429 AliKink *pkink = new AliKink;
5431 AliTPCseed & mother = *pmother;
5432 AliTPCseed & daughter = *pdaughter;
5433 AliKink & kinkl = *pkink;
5434 if (CheckKinkPoint(track0,mother,daughter, kinkl)){
5435 if (mother.GetNumberOfClusters()<30||daughter.GetNumberOfClusters()<20) {
5439 continue; //too short tracks
5441 if (mother.Pt()<1.4) {
5447 Int_t row0= kinkl.GetTPCRow0();
5448 if (kinkl.GetDistance()>0.5 || kinkl.GetR()<110. || kinkl.GetR()>240.) {
5455 Int_t index = esd->AddKink(&kinkl);
5456 mother.SetKinkIndex(0,-(index+1));
5457 daughter.SetKinkIndex(0,index+1);
5458 if (mother.GetNumberOfClusters()>50) {
5459 MarkSeedFree( array->RemoveAt(i) );
5460 AliTPCseed* mtc = new( NextFreeSeed() ) AliTPCseed(mother);
5461 mtc->SetPoolID(fLastSeedID);
5462 array->AddAt(mtc,i);
5465 AliTPCseed* mtc = new( NextFreeSeed() ) AliTPCseed(mother);
5466 mtc->SetPoolID(fLastSeedID);
5467 array->AddLast(mtc);
5469 AliTPCseed* dtc = new( NextFreeSeed() ) AliTPCseed(daughter);
5470 dtc->SetPoolID(fLastSeedID);
5471 array->AddLast(dtc);
5472 for (Int_t icl=0;icl<row0;icl++) {
5473 if (mother.GetClusterPointer(icl)) mother.GetClusterPointer(icl)->Use(20);
5476 for (Int_t icl=row0;icl<158;icl++) {
5477 if (daughter.GetClusterPointer(icl)) daughter.GetClusterPointer(icl)->Use(20);
5486 delete [] daughters;
5508 AliInfo(Form("Ncandidates=\t%d\t%d\t%d\t%d\n",esd->GetNumberOfKinks(),ncandidates,ntracks,nall));
5514 void AliTPCtrackerMI::FindKinks(TObjArray * array, AliESDEvent *esd)
5521 TObjArray *kinks= new TObjArray(10000);
5522 // TObjArray *v0s= new TObjArray(10000);
5523 Int_t nentries = array->GetEntriesFast();
5524 AliHelix *helixes = new AliHelix[nentries];
5525 Int_t *sign = new Int_t[nentries];
5526 Int_t *nclusters = new Int_t[nentries];
5527 Float_t *alpha = new Float_t[nentries];
5528 AliKink *kink = new AliKink();
5529 Int_t * usage = new Int_t[nentries];
5530 Float_t *zm = new Float_t[nentries];
5531 Float_t *z0 = new Float_t[nentries];
5532 Float_t *fim = new Float_t[nentries];
5533 Float_t *shared = new Float_t[nentries];
5534 Bool_t *circular = new Bool_t[nentries];
5535 Float_t *dca = new Float_t[nentries];
5536 //const AliESDVertex * primvertex = esd->GetVertex();
5538 // nentries = array->GetEntriesFast();
5543 for (Int_t i=0;i<nentries;i++){
5546 AliTPCseed* track = (AliTPCseed*)array->At(i);
5547 if (!track) continue;
5548 track->SetCircular(0);
5550 track->UpdatePoints();
5551 if (( track->GetPoints()[2]- track->GetPoints()[0])>5 && track->GetPoints()[3]>0.8){
5553 nclusters[i]=track->GetNumberOfClusters();
5554 alpha[i] = track->GetAlpha();
5555 new (&helixes[i]) AliHelix(*track);
5557 helixes[i].Evaluate(0,xyz);
5558 sign[i] = (track->GetC()>0) ? -1:1;
5561 if (track->GetProlongation(x,y,z)){
5563 fim[i] = alpha[i]+TMath::ATan2(y,x);
5566 zm[i] = track->GetZ();
5570 circular[i]= kFALSE;
5571 if (track->GetProlongation(0,y,z)) z0[i] = z;
5572 dca[i] = track->GetD(0,0);
5578 Int_t ncandidates =0;
5581 Double_t phase[2][2]={{0,0},{0,0}},radius[2]={0,0};
5584 // Find circling track
5586 for (Int_t i0=0;i0<nentries;i0++){
5587 AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
5588 if (!track0) continue;
5589 if (track0->GetNumberOfClusters()<40) continue;
5590 if (TMath::Abs(1./track0->GetC())>200) continue;
5591 for (Int_t i1=i0+1;i1<nentries;i1++){
5592 AliTPCseed * track1 = (AliTPCseed*)array->At(i1);
5593 if (!track1) continue;
5594 if (track1->GetNumberOfClusters()<40) continue;
5595 if ( TMath::Abs(track1->GetTgl()+track0->GetTgl())>0.1) continue;
5596 if (track0->GetBConstrain()&&track1->GetBConstrain()) continue;
5597 if (TMath::Abs(1./track1->GetC())>200) continue;
5598 if (track1->GetSigned1Pt()*track0->GetSigned1Pt()>0) continue;
5599 if (track1->GetTgl()*track0->GetTgl()>0) continue;
5600 if (TMath::Max(TMath::Abs(1./track0->GetC()),TMath::Abs(1./track1->GetC()))>190) continue;
5601 if (track0->GetBConstrain()&&track1->OneOverPt()<track0->OneOverPt()) continue; //returning - lower momenta
5602 if (track1->GetBConstrain()&&track0->OneOverPt()<track1->OneOverPt()) continue; //returning - lower momenta
5604 Float_t mindcar = TMath::Min(TMath::Abs(dca[i0]),TMath::Abs(dca[i1]));
5605 if (mindcar<5) continue;
5606 Float_t mindcaz = TMath::Min(TMath::Abs(z0[i0]-GetZ()),TMath::Abs(z0[i1]-GetZ()));
5607 if (mindcaz<5) continue;
5608 if (mindcar+mindcaz<20) continue;
5611 Float_t xc0 = helixes[i0].GetHelix(6);
5612 Float_t yc0 = helixes[i0].GetHelix(7);
5613 Float_t r0 = helixes[i0].GetHelix(8);
5614 Float_t xc1 = helixes[i1].GetHelix(6);
5615 Float_t yc1 = helixes[i1].GetHelix(7);
5616 Float_t r1 = helixes[i1].GetHelix(8);
5618 Float_t rmean = (r0+r1)*0.5;
5619 Float_t delta =TMath::Sqrt((xc1-xc0)*(xc1-xc0)+(yc1-yc0)*(yc1-yc0));
5620 //if (delta>30) continue;
5621 if (delta>rmean*0.25) continue;
5622 if (TMath::Abs(r0-r1)/rmean>0.3) continue;
5624 Int_t npoints = helixes[i0].GetRPHIintersections(helixes[i1], phase, radius,10);
5625 if (npoints==0) continue;
5626 helixes[i0].GetClosestPhases(helixes[i1], phase);
5630 Double_t hangles[3];
5631 helixes[i0].Evaluate(phase[0][0],xyz0);
5632 helixes[i1].Evaluate(phase[0][1],xyz1);
5634 helixes[i0].GetAngle(phase[0][0],helixes[i1],phase[0][1],hangles);
5635 Double_t deltah[2],deltabest;
5636 if (hangles[2]<2.8) continue;
5639 helixes[i0].ParabolicDCA(helixes[i1],phase[0][0],phase[0][1],radius[0],deltah[0],2);
5641 helixes[i0].ParabolicDCA(helixes[i1],phase[1][0],phase[1][1],radius[1],deltah[1],2);
5642 if (deltah[1]<deltah[0]) ibest=1;
5644 deltabest = TMath::Sqrt(deltah[ibest]);
5645 helixes[i0].Evaluate(phase[ibest][0],xyz0);
5646 helixes[i1].Evaluate(phase[ibest][1],xyz1);
5647 helixes[i0].GetAngle(phase[ibest][0],helixes[i1],phase[ibest][1],hangles);
5648 Double_t radiusbest = TMath::Sqrt(radius[ibest]);
5650 if (deltabest>6) continue;
5651 if (mindcar+mindcaz<40 && (hangles[2]<3.12||deltabest>3)) continue;
5652 Bool_t lsign =kFALSE;
5653 if (hangles[2]>3.06) lsign =kTRUE;
5656 circular[i0] = kTRUE;
5657 circular[i1] = kTRUE;
5658 if (track0->OneOverPt()<track1->OneOverPt()){
5659 track0->SetCircular(track0->GetCircular()+1);
5660 track1->SetCircular(track1->GetCircular()+2);
5663 track1->SetCircular(track1->GetCircular()+1);
5664 track0->SetCircular(track0->GetCircular()+2);
5667 if (lsign&&AliTPCReconstructor::StreamLevel()>1){
5669 Int_t lab0=track0->GetLabel();
5670 Int_t lab1=track1->GetLabel();
5671 TTreeSRedirector &cstream = *fDebugStreamer;
5672 cstream<<"Curling"<<
5679 "mindcar="<<mindcar<<
5680 "mindcaz="<<mindcaz<<
5683 "npoints="<<npoints<<
5684 "hangles0="<<hangles[0]<<
5685 "hangles2="<<hangles[2]<<
5690 "radius="<<radiusbest<<
5691 "deltabest="<<deltabest<<
5692 "phase0="<<phase[ibest][0]<<
5693 "phase1="<<phase[ibest][1]<<
5703 for (Int_t i =0;i<nentries;i++){
5704 if (sign[i]==0) continue;
5705 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5712 Double_t cradius0 = 40*40;
5713 Double_t cradius1 = 270*270;
5716 Double_t cdist3=0.55;
5717 for (Int_t j =i+1;j<nentries;j++){
5719 if (sign[j]*sign[i]<1) continue;
5720 if ( (nclusters[i]+nclusters[j])>200) continue;
5721 if ( (nclusters[i]+nclusters[j])<80) continue;
5722 if ( TMath::Abs(zm[i]-zm[j])>60.) continue;
5723 if ( TMath::Abs(fim[i]-fim[j])>0.6 && TMath::Abs(fim[i]-fim[j])<5.7 ) continue;
5724 //AliTPCseed * track1 = (AliTPCseed*)array->At(j); Double_t phase[2][2],radius[2];
5725 Int_t npoints = helixes[i].GetRPHIintersections(helixes[j], phase, radius,20);
5726 if (npoints<1) continue;
5729 if (radius[0]<cradius0||radius[0]>cradius1) continue;
5732 if ( (radius[0]<cradius0||radius[0]>cradius1) && (radius[1]<cradius0||radius[1]>cradius1) ) continue;
5735 Double_t delta1=10000,delta2=10000;
5736 // cuts on the intersection radius
5737 helixes[i].LinearDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta1);
5738 if (radius[0]<20&&delta1<1) continue; //intersection at vertex
5739 if (radius[0]<10&&delta1<3) continue; //intersection at vertex
5741 helixes[i].LinearDCA(helixes[j],phase[1][0],phase[1][1],radius[1],delta2);
5742 if (radius[1]<20&&delta2<1) continue; //intersection at vertex
5743 if (radius[1]<10&&delta2<3) continue; //intersection at vertex
5746 Double_t distance1 = TMath::Min(delta1,delta2);
5747 if (distance1>cdist1) continue; // cut on DCA linear approximation
5749 npoints = helixes[i].GetRPHIintersections(helixes[j], phase, radius,20);
5750 helixes[i].ParabolicDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta1);
5751 if (radius[0]<20&&delta1<1) continue; //intersection at vertex
5752 if (radius[0]<10&&delta1<3) continue; //intersection at vertex
5755 helixes[i].ParabolicDCA(helixes[j],phase[1][0],phase[1][1],radius[1],delta2);
5756 if (radius[1]<20&&delta2<1) continue; //intersection at vertex
5757 if (radius[1]<10&&delta2<3) continue; //intersection at vertex
5759 distance1 = TMath::Min(delta1,delta2);
5762 rkink = TMath::Sqrt(radius[0]);
5765 rkink = TMath::Sqrt(radius[1]);
5767 if (distance1>cdist2) continue;
5770 AliTPCseed * track1 = (AliTPCseed*)array->At(j);
5773 Int_t row0 = GetRowNumber(rkink);
5774 if (row0<10) continue;
5775 if (row0>150) continue;
5778 Float_t dens00=-1,dens01=-1;
5779 Float_t dens10=-1,dens11=-1;
5781 Int_t found,foundable,ishared;
5782 track0->GetClusterStatistic(0,row0-5, found, foundable,ishared,kFALSE);
5783 if (foundable>5) dens00 = Float_t(found)/Float_t(foundable);
5784 track0->GetClusterStatistic(row0+5,155, found, foundable,ishared,kFALSE);
5785 if (foundable>5) dens01 = Float_t(found)/Float_t(foundable);
5787 track1->GetClusterStatistic(0,row0-5, found, foundable,ishared,kFALSE);
5788 if (foundable>10) dens10 = Float_t(found)/Float_t(foundable);
5789 track1->GetClusterStatistic(row0+5,155, found, foundable,ishared,kFALSE);
5790 if (foundable>10) dens11 = Float_t(found)/Float_t(foundable);
5792 if (dens00<dens10 && dens01<dens11) continue;
5793 if (dens00>dens10 && dens01>dens11) continue;
5794 if (TMath::Max(dens00,dens10)<0.1) continue;
5795 if (TMath::Max(dens01,dens11)<0.3) continue;
5797 if (TMath::Min(dens00,dens10)>0.6) continue;
5798 if (TMath::Min(dens01,dens11)>0.6) continue;
5801 AliTPCseed * ktrack0, *ktrack1;
5810 if (TMath::Abs(ktrack0->GetC())>5) continue; // cut on the curvature for mother particle
5811 AliExternalTrackParam paramm(*ktrack0);
5812 AliExternalTrackParam paramd(*ktrack1);
5813 if (row0>60&&ktrack1->GetReference().GetX()>90.)new (¶md) AliExternalTrackParam(ktrack1->GetReference());
5816 kink->SetMother(paramm);
5817 kink->SetDaughter(paramd);
5820 Float_t x[3] = { kink->GetPosition()[0],kink->GetPosition()[1],kink->GetPosition()[2]};
5822 fkParam->Transform0to1(x,index);
5823 fkParam->Transform1to2(x,index);
5824 row0 = GetRowNumber(x[0]);
5826 if (kink->GetR()<100) continue;
5827 if (kink->GetR()>240) continue;
5828 if (kink->GetPosition()[2]/kink->GetR()>AliTPCReconstructor::GetCtgRange()) continue; //out of fiducial volume
5829 if (kink->GetDistance()>cdist3) continue;
5830 Float_t dird = kink->GetDaughterP()[0]*kink->GetPosition()[0]+kink->GetDaughterP()[1]*kink->GetPosition()[1]; // rough direction estimate
5831 if (dird<0) continue;
5833 Float_t dirm = kink->GetMotherP()[0]*kink->GetPosition()[0]+kink->GetMotherP()[1]*kink->GetPosition()[1]; // rough direction estimate
5834 if (dirm<0) continue;
5835 Float_t mpt = TMath::Sqrt(kink->GetMotherP()[0]*kink->GetMotherP()[0]+kink->GetMotherP()[1]*kink->GetMotherP()[1]);
5836 if (mpt<0.2) continue;
5839 //for high momenta momentum not defined well in first iteration
5840 Double_t qt = TMath::Sin(kink->GetAngle(2))*ktrack1->GetP();
5841 if (qt>0.35) continue;
5844 kink->SetLabel(CookLabel(ktrack0,0.4,0,row0),0);
5845 kink->SetLabel(CookLabel(ktrack1,0.4,row0,160),1);
5847 kink->SetTPCDensity(dens00,0,0);
5848 kink->SetTPCDensity(dens01,0,1);
5849 kink->SetTPCDensity(dens10,1,0);
5850 kink->SetTPCDensity(dens11,1,1);
5851 kink->SetIndex(i,0);
5852 kink->SetIndex(j,1);
5855 kink->SetTPCDensity(dens10,0,0);
5856 kink->SetTPCDensity(dens11,0,1);
5857 kink->SetTPCDensity(dens00,1,0);
5858 kink->SetTPCDensity(dens01,1,1);
5859 kink->SetIndex(j,0);
5860 kink->SetIndex(i,1);
5863 if (mpt<1||kink->GetAngle(2)>0.1){
5864 // angle and densities not defined yet
5865 if (kink->GetTPCDensityFactor()<0.8) continue;
5866 if ((2-kink->GetTPCDensityFactor())*kink->GetDistance() >0.25) continue;
5867 if (kink->GetAngle(2)*ktrack0->GetP()<0.003) continue; //too small angle
5868 if (kink->GetAngle(2)>0.2&&kink->GetTPCDensityFactor()<1.15) continue;
5869 if (kink->GetAngle(2)>0.2&&kink->GetTPCDensity(0,1)>0.05) continue;
5871 Float_t criticalangle = track0->GetSigmaSnp2()+track0->GetSigmaTgl2();
5872 criticalangle+= track1->GetSigmaSnp2()+track1->GetSigmaTgl2();
5873 criticalangle= 3*TMath::Sqrt(criticalangle);
5874 if (criticalangle>0.02) criticalangle=0.02;
5875 if (kink->GetAngle(2)<criticalangle) continue;
5878 Int_t drow = Int_t(2.+0.5/(0.05+kink->GetAngle(2))); // overlap region defined
5879 Float_t shapesum =0;
5881 for ( Int_t row = row0-drow; row<row0+drow;row++){
5882 if (row<0) continue;
5883 if (row>155) continue;
5884 if (ktrack0->GetClusterPointer(row)){
5885 AliTPCTrackerPoint *point =ktrack0->GetTrackPoint(row);
5886 shapesum+=point->GetSigmaY()+point->GetSigmaZ();
5889 if (ktrack1->GetClusterPointer(row)){
5890 AliTPCTrackerPoint *point =ktrack1->GetTrackPoint(row);
5891 shapesum+=point->GetSigmaY()+point->GetSigmaZ();
5896 kink->SetShapeFactor(-1.);
5899 kink->SetShapeFactor(shapesum/sum);
5901 // esd->AddKink(kink);
5903 // kink->SetMother(paramm);
5904 //kink->SetDaughter(paramd);
5906 Double_t chi2P2 = paramm.GetParameter()[2]-paramd.GetParameter()[2];
5908 chi2P2/=paramm.GetCovariance()[5]+paramd.GetCovariance()[5];
5909 Double_t chi2P3 = paramm.GetParameter()[3]-paramd.GetParameter()[3];
5911 chi2P3/=paramm.GetCovariance()[9]+paramd.GetCovariance()[9];
5913 if (AliTPCReconstructor::StreamLevel()>1) {
5914 (*fDebugStreamer)<<"kinkLpt"<<
5922 if ( chi2P2+chi2P3<AliTPCReconstructor::GetRecoParam()->GetKinkAngleCutChi2(0)){
5926 kinks->AddLast(kink);
5932 // sort the kinks according quality - and refit them towards vertex
5934 Int_t nkinks = kinks->GetEntriesFast();
5935 Float_t *quality = new Float_t[nkinks];
5936 Int_t *indexes = new Int_t[nkinks];
5937 AliTPCseed **mothers = new AliTPCseed*[nkinks]; memset(mothers, 0, nkinks*sizeof(AliTPCseed*));
5938 AliTPCseed **daughters = new AliTPCseed*[nkinks]; memset(daughters, 0, nkinks*sizeof(AliTPCseed*));
5941 for (Int_t i=0;i<nkinks;i++){
5943 AliKink *kinkl = (AliKink*)kinks->At(i);
5945 // refit kinks towards vertex
5947 Int_t index0 = kinkl->GetIndex(0);
5948 Int_t index1 = kinkl->GetIndex(1);
5949 AliTPCseed * ktrack0 = (AliTPCseed*)array->At(index0);
5950 AliTPCseed * ktrack1 = (AliTPCseed*)array->At(index1);
5952 Int_t sumn=ktrack0->GetNumberOfClusters()+ktrack1->GetNumberOfClusters();
5954 // Refit Kink under if too small angle
5956 if (kinkl->GetAngle(2)<0.05){
5957 kinkl->SetTPCRow0(GetRowNumber(kinkl->GetR()));
5958 Int_t row0 = kinkl->GetTPCRow0();
5959 Int_t drow = Int_t(2.+0.5/(0.05+kinkl->GetAngle(2)));
5962 Int_t last = row0-drow;
5963 if (last<40) last=40;
5964 if (last<ktrack0->GetFirstPoint()+25) last = ktrack0->GetFirstPoint()+25;
5965 AliTPCseed* seed0 = ReSeed(ktrack0,last,kFALSE);
5968 Int_t first = row0+drow;
5969 if (first>130) first=130;
5970 if (first>ktrack1->GetLastPoint()-25) first = TMath::Max(ktrack1->GetLastPoint()-25,30);
5971 AliTPCseed* seed1 = ReSeed(ktrack1,first,kTRUE);
5973 if (seed0 && seed1){
5974 kinkl->SetStatus(1,8);
5975 if (RefitKink(*seed0,*seed1,*kinkl)) kinkl->SetStatus(1,9);
5976 row0 = GetRowNumber(kinkl->GetR());
5977 sumn = seed0->GetNumberOfClusters()+seed1->GetNumberOfClusters();
5978 mothers[i] = new ( NextFreeSeed() ) AliTPCseed(*seed0);
5979 mothers[i]->SetPoolID(fLastSeedID);
5980 daughters[i] = new (NextFreeSeed() ) AliTPCseed(*seed1);
5981 daughters[i]->SetPoolID(fLastSeedID);
5984 delete kinks->RemoveAt(i);
5985 if (seed0) MarkSeedFree( seed0 );
5986 if (seed1) MarkSeedFree( seed1 );
5989 if (kinkl->GetDistance()>0.5 || kinkl->GetR()<110 || kinkl->GetR()>240) {
5990 delete kinks->RemoveAt(i);
5991 if (seed0) MarkSeedFree( seed0 );
5992 if (seed1) MarkSeedFree( seed1 );
5996 MarkSeedFree( seed0 );
5997 MarkSeedFree( seed1 );
6000 if (kinkl) quality[i] = 160*((0.1+kinkl->GetDistance())*(2.-kinkl->GetTPCDensityFactor()))/(sumn+40.); //the longest -clossest will win
6002 TMath::Sort(nkinks,quality,indexes,kFALSE);
6004 //remove double find kinks
6006 for (Int_t ikink0=1;ikink0<nkinks;ikink0++){
6007 AliKink * kink0 = (AliKink*) kinks->At(indexes[ikink0]);
6008 if (!kink0) continue;
6010 for (Int_t ikink1=0;ikink1<ikink0;ikink1++){
6011 kink0 = (AliKink*) kinks->At(indexes[ikink0]);
6012 if (!kink0) continue;
6013 AliKink * kink1 = (AliKink*) kinks->At(indexes[ikink1]);
6014 if (!kink1) continue;
6015 // if not close kink continue
6016 if (TMath::Abs(kink1->GetPosition()[2]-kink0->GetPosition()[2])>10) continue;
6017 if (TMath::Abs(kink1->GetPosition()[1]-kink0->GetPosition()[1])>10) continue;
6018 if (TMath::Abs(kink1->GetPosition()[0]-kink0->GetPosition()[0])>10) continue;
6020 AliTPCseed &mother0 = *mothers[indexes[ikink0]];
6021 AliTPCseed &daughter0 = *daughters[indexes[ikink0]];
6022 AliTPCseed &mother1 = *mothers[indexes[ikink1]];
6023 AliTPCseed &daughter1 = *daughters[indexes[ikink1]];
6024 Int_t row0 = (kink0->GetTPCRow0()+kink1->GetTPCRow0())/2;
6033 for (Int_t i=0;i<row0;i++){
6034 if (mother0.GetClusterIndex(i)>0 && mother1.GetClusterIndex(i)>0){
6037 if (mother0.GetClusterIndex(i)==mother1.GetClusterIndex(i)){
6044 for (Int_t i=row0;i<158;i++){
6045 //if (daughter0.GetClusterIndex(i)>0 && daughter0.GetClusterIndex(i)>0){ // RS: Bug?
6046 if (daughter0.GetClusterIndex(i)>0 && daughter1.GetClusterIndex(i)>0){
6049 if (mother0.GetClusterIndex(i)==mother1.GetClusterIndex(i)){
6055 Float_t ratio = Float_t(same+1)/Float_t(both+1);
6056 Float_t ratiom = Float_t(samem+1)/Float_t(bothm+1);
6057 Float_t ratiod = Float_t(samed+1)/Float_t(bothd+1);
6058 if (ratio>0.3 && ratiom>0.5 &&ratiod>0.5) {
6059 Int_t sum0 = mother0.GetNumberOfClusters()+daughter0.GetNumberOfClusters();
6060 Int_t sum1 = mother1.GetNumberOfClusters()+daughter1.GetNumberOfClusters();
6062 shared[kink0->GetIndex(0)]= kTRUE;
6063 shared[kink0->GetIndex(1)]= kTRUE;
6064 delete kinks->RemoveAt(indexes[ikink0]);
6068 shared[kink1->GetIndex(0)]= kTRUE;
6069 shared[kink1->GetIndex(1)]= kTRUE;
6070 delete kinks->RemoveAt(indexes[ikink1]);
6077 for (Int_t i=0;i<nkinks;i++){
6078 AliKink * kinkl = (AliKink*) kinks->At(indexes[i]);
6079 if (!kinkl) continue;
6080 kinkl->SetTPCRow0(GetRowNumber(kinkl->GetR()));
6081 Int_t index0 = kinkl->GetIndex(0);
6082 Int_t index1 = kinkl->GetIndex(1);
6083 if (circular[index0]||(circular[index1]&&kinkl->GetDistance()>0.2)) continue;
6084 kinkl->SetMultiple(usage[index0],0);
6085 kinkl->SetMultiple(usage[index1],1);
6086 if (kinkl->GetMultiple()[0]+kinkl->GetMultiple()[1]>2) continue;
6087 if (kinkl->GetMultiple()[0]+kinkl->GetMultiple()[1]>0 && quality[indexes[i]]>0.2) continue;
6088 if (kinkl->GetMultiple()[0]+kinkl->GetMultiple()[1]>0 && kinkl->GetDistance()>0.2) continue;
6089 if (circular[index0]||(circular[index1]&&kinkl->GetDistance()>0.1)) continue;
6091 AliTPCseed * ktrack0 = (AliTPCseed*)array->At(index0);
6092 AliTPCseed * ktrack1 = (AliTPCseed*)array->At(index1);
6093 if (!ktrack0 || !ktrack1) continue;
6094 Int_t index = esd->AddKink(kinkl);
6097 if ( ktrack0->GetKinkIndex(0)==0 && ktrack1->GetKinkIndex(0)==0) { //best kink
6098 if (mothers[indexes[i]]->GetNumberOfClusters()>20 && daughters[indexes[i]]->GetNumberOfClusters()>20 &&
6099 (mothers[indexes[i]]->GetNumberOfClusters()+daughters[indexes[i]]->GetNumberOfClusters())>100){
6100 *ktrack0 = *mothers[indexes[i]];
6101 *ktrack1 = *daughters[indexes[i]];
6105 ktrack0->SetKinkIndex(usage[index0],-(index+1));
6106 ktrack1->SetKinkIndex(usage[index1], (index+1));
6111 // Remove tracks corresponding to shared kink's
6113 for (Int_t i=0;i<nentries;i++){
6114 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
6115 if (!track0) continue;
6116 if (track0->GetKinkIndex(0)!=0) continue;
6117 if (shared[i]) MarkSeedFree( array->RemoveAt(i) );
6122 RemoveUsed2(array,0.5,0.4,30);
6124 for (Int_t i=0;i<nentries;i++){
6125 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
6126 if (!track0) continue;
6127 track0->CookdEdx(0.02,0.6);
6131 for (Int_t i=0;i<nentries;i++){
6132 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
6133 if (!track0) continue;
6134 if (track0->Pt()<1.4) continue;
6135 //remove double high momenta tracks - overlapped with kink candidates
6138 for (Int_t icl=track0->GetFirstPoint();icl<track0->GetLastPoint(); icl++){
6139 if (track0->GetClusterPointer(icl)!=0){
6141 if (track0->GetClusterPointer(icl)->IsUsed(10)) ishared++;
6144 if (Float_t(ishared+1)/Float_t(all+1)>0.5) {
6145 MarkSeedFree( array->RemoveAt(i) );
6149 if (track0->GetKinkIndex(0)!=0) continue;
6150 if (track0->GetNumberOfClusters()<80) continue;
6152 AliTPCseed *pmother = new( NextFreeSeed() ) AliTPCseed();
6153 pmother->SetPoolID(fLastSeedID);
6154 AliTPCseed *pdaughter = new( NextFreeSeed() ) AliTPCseed();
6155 pdaughter->SetPoolID(fLastSeedID);
6156 AliKink *pkink = new AliKink;
6158 AliTPCseed & mother = *pmother;
6159 AliTPCseed & daughter = *pdaughter;
6160 AliKink & kinkl = *pkink;
6161 if (CheckKinkPoint(track0,mother,daughter, kinkl)){
6162 if (mother.GetNumberOfClusters()<30||daughter.GetNumberOfClusters()<20) {
6163 MarkSeedFree( pmother );
6164 MarkSeedFree( pdaughter );
6166 continue; //too short tracks
6168 if (mother.Pt()<1.4) {
6169 MarkSeedFree( pmother );
6170 MarkSeedFree( pdaughter );
6174 Int_t row0= kinkl.GetTPCRow0();
6175 if (kinkl.GetDistance()>0.5 || kinkl.GetR()<110. || kinkl.GetR()>240.) {
6176 MarkSeedFree( pmother );
6177 MarkSeedFree( pdaughter );
6182 Int_t index = esd->AddKink(&kinkl);
6183 mother.SetKinkIndex(0,-(index+1));
6184 daughter.SetKinkIndex(0,index+1);
6185 if (mother.GetNumberOfClusters()>50) {
6186 MarkSeedFree( array->RemoveAt(i) );
6187 AliTPCseed* mtc = new( NextFreeSeed() ) AliTPCseed(mother);
6188 mtc->SetPoolID(fLastSeedID);
6189 array->AddAt(mtc,i);
6192 AliTPCseed* mtc = new( NextFreeSeed() ) AliTPCseed(mother);
6193 mtc->SetPoolID(fLastSeedID);
6194 array->AddLast(mtc);
6196 AliTPCseed* dtc = new( NextFreeSeed() ) AliTPCseed(daughter);
6197 dtc->SetPoolID(fLastSeedID);
6198 array->AddLast(dtc);
6199 for (Int_t icl=0;icl<row0;icl++) {
6200 if (mother.GetClusterPointer(icl)) mother.GetClusterPointer(icl)->Use(20);
6203 for (Int_t icl=row0;icl<158;icl++) {
6204 if (daughter.GetClusterPointer(icl)) daughter.GetClusterPointer(icl)->Use(20);
6208 MarkSeedFree( pmother );
6209 MarkSeedFree( pdaughter );
6213 delete [] daughters;
6235 AliInfo(Form("Ncandidates=\t%d\t%d\t%d\t%d\n",esd->GetNumberOfKinks(),ncandidates,ntracks,nall));
6240 Int_t AliTPCtrackerMI::RefitKink(AliTPCseed &mother, AliTPCseed &daughter, const AliESDkink &knk)
6243 // refit kink towards to the vertex
6246 AliKink &kink=(AliKink &)knk;
6248 Int_t row0 = GetRowNumber(kink.GetR());
6249 FollowProlongation(mother,0);
6250 mother.Reset(kFALSE);
6252 FollowProlongation(daughter,row0);
6253 daughter.Reset(kFALSE);
6254 FollowBackProlongation(daughter,158);
6255 daughter.Reset(kFALSE);
6256 Int_t first = TMath::Max(row0-20,30);
6257 Int_t last = TMath::Min(row0+20,140);
6259 const Int_t kNdiv =5;
6260 AliTPCseed param0[kNdiv]; // parameters along the track
6261 AliTPCseed param1[kNdiv]; // parameters along the track
6262 AliKink kinks[kNdiv]; // corresponding kink parameters
6265 for (Int_t irow=0; irow<kNdiv;irow++){
6266 rows[irow] = first +((last-first)*irow)/(kNdiv-1);
6268 // store parameters along the track
6270 for (Int_t irow=0;irow<kNdiv;irow++){
6271 FollowBackProlongation(mother, rows[irow]);
6272 FollowProlongation(daughter,rows[kNdiv-1-irow]);
6273 param0[irow] = mother;
6274 param1[kNdiv-1-irow] = daughter;
6278 for (Int_t irow=0; irow<kNdiv-1;irow++){
6279 if (param0[irow].GetNumberOfClusters()<kNdiv||param1[irow].GetNumberOfClusters()<kNdiv) continue;
6280 kinks[irow].SetMother(param0[irow]);
6281 kinks[irow].SetDaughter(param1[irow]);
6282 kinks[irow].Update();
6285 // choose kink with best "quality"
6287 Double_t mindist = 10000;
6288 for (Int_t irow=0;irow<kNdiv;irow++){
6289 if (param0[irow].GetNumberOfClusters()<20||param1[irow].GetNumberOfClusters()<20) continue;
6290 if (TMath::Abs(kinks[irow].GetR())>240.) continue;
6291 if (TMath::Abs(kinks[irow].GetR())<100.) continue;
6293 Float_t normdist = TMath::Abs(param0[irow].GetX()-kinks[irow].GetR())*(0.1+kink.GetDistance());
6294 normdist/= (param0[irow].GetNumberOfClusters()+param1[irow].GetNumberOfClusters()+40.);
6295 if (normdist < mindist){
6301 if (index==-1) return 0;
6304 param0[index].Reset(kTRUE);
6305 FollowProlongation(param0[index],0);
6307 mother = param0[index];
6308 daughter = param1[index]; // daughter in vertex
6310 kink.SetMother(mother);
6311 kink.SetDaughter(daughter);
6313 kink.SetTPCRow0(GetRowNumber(kink.GetR()));
6314 kink.SetTPCncls(param0[index].GetNumberOfClusters(),0);
6315 kink.SetTPCncls(param1[index].GetNumberOfClusters(),1);
6316 kink.SetLabel(CookLabel(&mother,0.4, 0,kink.GetTPCRow0()),0);
6317 kink.SetLabel(CookLabel(&daughter,0.4, kink.GetTPCRow0(),160),1);
6318 mother.SetLabel(kink.GetLabel(0));
6319 daughter.SetLabel(kink.GetLabel(1));
6325 void AliTPCtrackerMI::UpdateKinkQualityM(AliTPCseed * seed){
6327 // update Kink quality information for mother after back propagation
6329 if (seed->GetKinkIndex(0)>=0) return;
6330 for (Int_t ikink=0;ikink<3;ikink++){
6331 Int_t index = seed->GetKinkIndex(ikink);
6332 if (index>=0) break;
6333 index = TMath::Abs(index)-1;
6334 AliESDkink * kink = fEvent->GetKink(index);
6335 kink->SetTPCDensity(-1,0,0);
6336 kink->SetTPCDensity(1,0,1);
6338 Int_t row0 = kink->GetTPCRow0() - 2 - Int_t( 0.5/ (0.05+kink->GetAngle(2)));
6339 if (row0<15) row0=15;
6341 Int_t row1 = kink->GetTPCRow0() + 2 + Int_t( 0.5/ (0.05+kink->GetAngle(2)));
6342 if (row1>145) row1=145;
6344 Int_t found,foundable,shared;
6345 seed->GetClusterStatistic(0,row0, found, foundable,shared,kFALSE);
6346 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),0,0);
6347 seed->GetClusterStatistic(row1,155, found, foundable,shared,kFALSE);
6348 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),0,1);
6353 void AliTPCtrackerMI::UpdateKinkQualityD(AliTPCseed * seed){
6355 // update Kink quality information for daughter after refit
6357 if (seed->GetKinkIndex(0)<=0) return;
6358 for (Int_t ikink=0;ikink<3;ikink++){
6359 Int_t index = seed->GetKinkIndex(ikink);
6360 if (index<=0) break;
6361 index = TMath::Abs(index)-1;
6362 AliESDkink * kink = fEvent->GetKink(index);
6363 kink->SetTPCDensity(-1,1,0);
6364 kink->SetTPCDensity(-1,1,1);
6366 Int_t row0 = kink->GetTPCRow0() -2 - Int_t( 0.5/ (0.05+kink->GetAngle(2)));
6367 if (row0<15) row0=15;
6369 Int_t row1 = kink->GetTPCRow0() +2 + Int_t( 0.5/ (0.05+kink->GetAngle(2)));
6370 if (row1>145) row1=145;
6372 Int_t found,foundable,shared;
6373 seed->GetClusterStatistic(0,row0, found, foundable,shared,kFALSE);
6374 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),1,0);
6375 seed->GetClusterStatistic(row1,155, found, foundable,shared,kFALSE);
6376 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),1,1);
6382 Int_t AliTPCtrackerMI::CheckKinkPoint(AliTPCseed*seed,AliTPCseed &mother, AliTPCseed &daughter, const AliESDkink &knk)
6385 // check kink point for given track
6386 // if return value=0 kink point not found
6387 // otherwise seed0 correspond to mother particle
6388 // seed1 correspond to daughter particle
6389 // kink parameter of kink point
6390 AliKink &kink=(AliKink &)knk;
6392 Int_t middlerow = (seed->GetFirstPoint()+seed->GetLastPoint())/2;
6393 Int_t first = seed->GetFirstPoint();
6394 Int_t last = seed->GetLastPoint();
6395 if (last-first<20) return 0; // shortest length - 2*30 = 60 pad-rows
6398 AliTPCseed *seed1 = ReSeed(seed,middlerow+20, kTRUE); //middle of chamber
6399 if (!seed1) return 0;
6400 FollowProlongation(*seed1,seed->GetLastPoint()-20);
6401 seed1->Reset(kTRUE);
6402 FollowProlongation(*seed1,158);
6403 seed1->Reset(kTRUE);
6404 last = seed1->GetLastPoint();
6406 AliTPCseed *seed0 = new( NextFreeSeed() ) AliTPCseed(*seed);
6407 seed0->SetPoolID(fLastSeedID);
6408 seed0->Reset(kFALSE);
6411 AliTPCseed param0[20]; // parameters along the track
6412 AliTPCseed param1[20]; // parameters along the track
6413 AliKink kinks[20]; // corresponding kink parameters
6415 for (Int_t irow=0; irow<20;irow++){
6416 rows[irow] = first +((last-first)*irow)/19;
6418 // store parameters along the track
6420 for (Int_t irow=0;irow<20;irow++){
6421 FollowBackProlongation(*seed0, rows[irow]);
6422 FollowProlongation(*seed1,rows[19-irow]);
6423 param0[irow] = *seed0;
6424 param1[19-irow] = *seed1;
6428 for (Int_t irow=0; irow<19;irow++){
6429 kinks[irow].SetMother(param0[irow]);
6430 kinks[irow].SetDaughter(param1[irow]);
6431 kinks[irow].Update();
6434 // choose kink with biggest change of angle
6436 Double_t maxchange= 0;
6437 for (Int_t irow=1;irow<19;irow++){
6438 if (TMath::Abs(kinks[irow].GetR())>240.) continue;
6439 if (TMath::Abs(kinks[irow].GetR())<110.) continue;
6440 Float_t quality = TMath::Abs(kinks[irow].GetAngle(2))/(3.+TMath::Abs(kinks[irow].GetR()-param0[irow].GetX()));
6441 if ( quality > maxchange){
6442 maxchange = quality;
6447 MarkSeedFree( seed0 );
6448 MarkSeedFree( seed1 );
6449 if (index<0) return 0;
6451 Int_t row0 = GetRowNumber(kinks[index].GetR()); //row 0 estimate
6452 seed0 = new( NextFreeSeed() ) AliTPCseed(param0[index]);
6453 seed0->SetPoolID(fLastSeedID);
6454 seed1 = new( NextFreeSeed() ) AliTPCseed(param1[index]);
6455 seed1->SetPoolID(fLastSeedID);
6456 seed0->Reset(kFALSE);
6457 seed1->Reset(kFALSE);
6458 seed0->ResetCovariance(10.);
6459 seed1->ResetCovariance(10.);
6460 FollowProlongation(*seed0,0);
6461 FollowBackProlongation(*seed1,158);
6462 mother = *seed0; // backup mother at position 0
6463 seed0->Reset(kFALSE);
6464 seed1->Reset(kFALSE);
6465 seed0->ResetCovariance(10.);
6466 seed1->ResetCovariance(10.);
6468 first = TMath::Max(row0-20,0);
6469 last = TMath::Min(row0+20,158);
6471 for (Int_t irow=0; irow<20;irow++){
6472 rows[irow] = first +((last-first)*irow)/19;
6474 // store parameters along the track
6476 for (Int_t irow=0;irow<20;irow++){
6477 FollowBackProlongation(*seed0, rows[irow]);
6478 FollowProlongation(*seed1,rows[19-irow]);
6479 param0[irow] = *seed0;
6480 param1[19-irow] = *seed1;
6484 for (Int_t irow=0; irow<19;irow++){
6485 kinks[irow].SetMother(param0[irow]);
6486 kinks[irow].SetDaughter(param1[irow]);
6487 // param0[irow].Dump();
6488 //param1[irow].Dump();
6489 kinks[irow].Update();
6492 // choose kink with biggest change of angle
6495 for (Int_t irow=0;irow<20;irow++){
6496 if (TMath::Abs(kinks[irow].GetR())>250.) continue;
6497 if (TMath::Abs(kinks[irow].GetR())<90.) continue;
6498 Float_t quality = TMath::Abs(kinks[irow].GetAngle(2))/(3.+TMath::Abs(kinks[irow].GetR()-param0[irow].GetX()));
6499 if ( quality > maxchange){
6500 maxchange = quality;
6507 if (index==-1 || param0[index].GetNumberOfClusters()+param1[index].GetNumberOfClusters()<100){
6508 MarkSeedFree( seed0 );
6509 MarkSeedFree( seed1 );
6513 // Float_t anglesigma = TMath::Sqrt(param0[index].fC22+param0[index].fC33+param1[index].fC22+param1[index].fC33);
6515 kink.SetMother(param0[index]);
6516 kink.SetDaughter(param1[index]);
6519 Double_t chi2P2 = param0[index].GetParameter()[2]-param1[index].GetParameter()[2];
6521 chi2P2/=param0[index].GetCovariance()[5]+param1[index].GetCovariance()[5];
6522 Double_t chi2P3 = param0[index].GetParameter()[3]-param1[index].GetParameter()[3];
6524 chi2P3/=param0[index].GetCovariance()[9]+param1[index].GetCovariance()[9];
6526 if (AliTPCReconstructor::StreamLevel()>1) {
6527 (*fDebugStreamer)<<"kinkHpt"<<
6530 "p0.="<<¶m0[index]<<
6531 "p1.="<<¶m1[index]<<
6535 if ( chi2P2+chi2P3<AliTPCReconstructor::GetRecoParam()->GetKinkAngleCutChi2(0)){
6536 MarkSeedFree( seed0 );
6537 MarkSeedFree( seed1 );
6542 row0 = GetRowNumber(kink.GetR());
6543 kink.SetTPCRow0(row0);
6544 kink.SetLabel(CookLabel(seed0,0.5,0,row0),0);
6545 kink.SetLabel(CookLabel(seed1,0.5,row0,158),1);
6546 kink.SetIndex(-10,0);
6547 kink.SetIndex(int(param0[index].GetNumberOfClusters()+param1[index].GetNumberOfClusters()),1);
6548 kink.SetTPCncls(param0[index].GetNumberOfClusters(),0);
6549 kink.SetTPCncls(param1[index].GetNumberOfClusters(),1);
6552 // new (&mother) AliTPCseed(param0[index]);
6553 daughter = param1[index];
6554 daughter.SetLabel(kink.GetLabel(1));
6555 param0[index].Reset(kTRUE);
6556 FollowProlongation(param0[index],0);
6557 mother = param0[index];
6558 mother.SetLabel(kink.GetLabel(0));
6559 if ( chi2P2+chi2P3<AliTPCReconstructor::GetRecoParam()->GetKinkAngleCutChi2(1)){
6562 MarkSeedFree( seed0 );
6563 MarkSeedFree( seed1 );
6571 AliTPCseed* AliTPCtrackerMI::ReSeed(AliTPCseed *t)
6574 // reseed - refit - track
6577 // Int_t last = fSectors->GetNRows()-1;
6579 if (fSectors == fOuterSec){
6580 first = TMath::Max(first, t->GetFirstPoint()-fInnerSec->GetNRows());
6584 first = t->GetFirstPoint();
6586 AliTPCseed * seed = MakeSeed(t,0.1,0.5,0.9);
6587 FollowBackProlongation(*t,fSectors->GetNRows()-1);
6589 FollowProlongation(*t,first);
6599 //_____________________________________________________________________________
6600 Int_t AliTPCtrackerMI::ReadSeeds(const TFile *inp) {
6601 //-----------------------------------------------------------------
6602 // This function reades track seeds.
6603 //-----------------------------------------------------------------
6604 TDirectory *savedir=gDirectory;
6606 TFile *in=(TFile*)inp;
6607 if (!in->IsOpen()) {
6608 cerr<<"AliTPCtrackerMI::ReadSeeds(): input file is not open !\n";
6613 TTree *seedTree=(TTree*)in->Get("Seeds");
6615 cerr<<"AliTPCtrackerMI::ReadSeeds(): ";
6616 cerr<<"can't get a tree with track seeds !\n";
6619 AliTPCtrack *seed=new AliTPCtrack;
6620 seedTree->SetBranchAddress("tracks",&seed);
6622 if (fSeeds==0) fSeeds=new TObjArray(15000);
6624 Int_t n=(Int_t)seedTree->GetEntries();
6625 for (Int_t i=0; i<n; i++) {
6626 seedTree->GetEvent(i);
6627 AliTPCseed* sdc = new( NextFreeSeed() ) AliTPCseed(*seed/*,seed->GetAlpha()*/);
6628 sdc->SetPoolID(fLastSeedID);
6629 fSeeds->AddLast(sdc);
6632 delete seed; // RS: this seed is not from the pool, delete it !!!
6638 Int_t AliTPCtrackerMI::Clusters2Tracks (AliESDEvent *const esd)
6641 // clusters to tracks
6642 if (fSeeds) DeleteSeeds();
6643 else ResetSeedsPool();
6646 AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ;
6647 transform->SetCurrentTimeStamp( esd->GetTimeStamp());
6648 transform->SetCurrentRun(esd->GetRunNumber());
6651 if (!fSeeds) return 1;
6653 if (AliTPCReconstructor::StreamLevel()>3) DumpClusters(0,fSeeds);
6659 //_____________________________________________________________________________
6660 Int_t AliTPCtrackerMI::Clusters2Tracks() {
6661 //-----------------------------------------------------------------
6662 // This is a track finder.
6663 //-----------------------------------------------------------------
6664 TDirectory *savedir=gDirectory;
6668 fSeeds = Tracking();
6671 Info("Clusters2Tracks","Time for tracking: \t");timer.Print();timer.Start();
6673 //activate again some tracks
6674 for (Int_t i=0; i<fSeeds->GetEntriesFast(); i++) {
6675 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
6677 Int_t nc=t.GetNumberOfClusters();
6679 MarkSeedFree( fSeeds->RemoveAt(i) );
6683 if (pt->GetRemoval()==10) {
6684 if (pt->GetDensityFirst(20)>0.8 || pt->GetDensityFirst(30)>0.8 || pt->GetDensityFirst(40)>0.7)
6685 pt->Desactivate(10); // make track again active // MvL: should be 0 ?
6687 pt->Desactivate(20);
6688 MarkSeedFree( fSeeds->RemoveAt(i) );
6693 RemoveUsed2(fSeeds,0.85,0.85,0);
6694 if (AliTPCReconstructor::GetRecoParam()->GetDoKinks()) FindKinks(fSeeds,fEvent);
6695 //FindCurling(fSeeds, fEvent,0);
6696 if (AliTPCReconstructor::StreamLevel()>5) FindMultiMC(fSeeds, fEvent,-1); // find multi found tracks
6697 RemoveUsed2(fSeeds,0.5,0.4,20);
6698 FindSplitted(fSeeds, fEvent,0); // find multi found tracks
6699 if (AliTPCReconstructor::StreamLevel()>5) FindMultiMC(fSeeds, fEvent,0); // find multi found tracks
6702 // // refit short tracks
6704 Int_t nseed=fSeeds->GetEntriesFast();
6707 for (Int_t i=0; i<nseed; i++) {
6708 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
6710 Int_t nc=t.GetNumberOfClusters();
6712 MarkSeedFree( fSeeds->RemoveAt(i) );
6715 CookLabel(pt,0.1); //For comparison only
6716 //if ((pt->IsActive() || (pt->fRemoval==10) )&& nc>50 &&pt->GetNumberOfClusters()>0.4*pt->fNFoundable){
6717 if ((pt->IsActive() || (pt->GetRemoval()==10) )){
6719 if (fDebug>0) cerr<<found<<'\r';
6723 MarkSeedFree( fSeeds->RemoveAt(i) );
6727 //RemoveOverlap(fSeeds,0.99,7,kTRUE);
6729 //RemoveUsed(fSeeds,0.9,0.9,6);
6731 nseed=fSeeds->GetEntriesFast();
6733 for (Int_t i=0; i<nseed; i++) {
6734 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
6736 Int_t nc=t.GetNumberOfClusters();
6738 MarkSeedFree( fSeeds->RemoveAt(i) );
6742 t.CookdEdx(0.02,0.6);
6743 // CheckKinkPoint(&t,0.05);
6744 //if ((pt->IsActive() || (pt->fRemoval==10) )&& nc>50 &&pt->GetNumberOfClusters()>0.4*pt->fNFoundable){
6745 if ((pt->IsActive() || (pt->GetRemoval()==10) )){
6753 MarkSeedFree( fSeeds->RemoveAt(i) );
6754 //AliTPCseed * seed1 = ReSeed(pt,0.05,0.5,1);
6756 // FollowProlongation(*seed1,0);
6757 // Int_t n = seed1->GetNumberOfClusters();
6758 // printf("fP4\t%f\t%f\n",seed1->GetC(),pt->GetC());
6759 // printf("fN\t%d\t%d\n", seed1->GetNumberOfClusters(),pt->GetNumberOfClusters());
6762 //AliTPCseed * seed2 = ReSeed(pt,0.95,0.5,0.05);
6766 SortTracks(fSeeds, 1);
6770 PrepareForBackProlongation(fSeeds,5.);
6771 PropagateBack(fSeeds);
6772 printf("Time for back propagation: \t");timer.Print();timer.Start();
6776 PrepareForProlongation(fSeeds,5.);
6777 PropagateForard2(fSeeds);
6779 printf("Time for FORWARD propagation: \t");timer.Print();timer.Start();
6780 // RemoveUsed(fSeeds,0.7,0.7,6);
6781 //RemoveOverlap(fSeeds,0.9,7,kTRUE);
6783 nseed=fSeeds->GetEntriesFast();
6785 for (Int_t i=0; i<nseed; i++) {
6786 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
6788 Int_t nc=t.GetNumberOfClusters();
6790 MarkSeedFree( fSeeds->RemoveAt(i) );
6793 t.CookdEdx(0.02,0.6);
6794 // CookLabel(pt,0.1); //For comparison only
6795 //if ((pt->IsActive() || (pt->fRemoval==10) )&& nc>50 &&pt->GetNumberOfClusters()>0.4*pt->fNFoundable){
6796 if ((pt->IsActive() || (pt->fRemoval==10) )){
6797 cerr<<found++<<'\r';
6800 MarkSeedFree( fSeeds->RemoveAt(i) );
6805 // fNTracks = found;
6807 Info("Clusters2Tracks","Time for overlap removal, track writing and dedx cooking: \t"); timer.Print();timer.Start();
6810 // cerr<<"Number of found tracks : "<<"\t"<<found<<endl;
6811 Info("Clusters2Tracks","Number of found tracks %d",found);
6813 // UnloadClusters();
6818 void AliTPCtrackerMI::Tracking(TObjArray * arr)
6821 // tracking of the seeds
6824 fSectors = fOuterSec;
6825 ParallelTracking(arr,150,63);
6826 fSectors = fOuterSec;
6827 ParallelTracking(arr,63,0);
6830 TObjArray * AliTPCtrackerMI::Tracking(Int_t seedtype, Int_t i1, Int_t i2, Float_t cuts[4], Float_t dy, Int_t dsec)
6835 static TObjArray arrTracks;
6836 TObjArray * arr = &arrTracks;
6838 fSectors = fOuterSec;
6841 for (Int_t sec=0;sec<fkNOS;sec++){
6842 if (seedtype==3) MakeSeeds3(arr,sec,i1,i2,cuts,dy, dsec);
6843 if (seedtype==4) MakeSeeds5(arr,sec,i1,i2,cuts,dy);
6844 if (seedtype==2) MakeSeeds2(arr,sec,i1,i2,cuts,dy);
6847 Info("Tracking","\nSeeding - %d\t%d\t%d\t%d\n",seedtype,i1,i2,arr->GetEntriesFast());
6859 TObjArray * AliTPCtrackerMI::Tracking()
6863 if (AliTPCReconstructor::GetRecoParam()->GetSpecialSeeding()) return TrackingSpecial();
6866 Int_t nup=fOuterSec->GetNRows()+fInnerSec->GetNRows();
6868 TObjArray * seeds = new TObjArray;
6870 Int_t fLastSeedRowSec=AliTPCReconstructor::GetRecoParam()->GetLastSeedRowSec();
6871 Int_t gapPrim = AliTPCReconstructor::GetRecoParam()->GetSeedGapPrim();
6872 Int_t gapSec = AliTPCReconstructor::GetRecoParam()->GetSeedGapSec();
6880 Float_t fnumber = 3.0;
6881 Float_t fdensity = 3.0;
6886 for (Int_t delta = 0; delta<18; delta+=gapPrim){
6890 arr = Tracking(3,nup-1-delta,nup-1-delta-gap,cuts,-1,1);
6891 SumTracks(seeds,arr);
6892 SignClusters(seeds,fnumber,fdensity);
6894 for (Int_t i=2;i<6;i+=2){
6895 // seed high pt tracks
6898 arr = Tracking(3,nup-i-delta,nup-i-delta-gap,cuts,-1,0);
6899 SumTracks(seeds,arr);
6900 SignClusters(seeds,fnumber,fdensity);
6905 // RemoveUsed(seeds,0.9,0.9,1);
6906 // UnsignClusters();
6907 // SignClusters(seeds,fnumber,fdensity);
6911 for (Int_t delta = 20; delta<120; delta+=gapPrim){
6913 // seed high pt tracks
6917 arr = Tracking(3,nup-delta,nup-delta-gap,cuts,-1);
6918 SumTracks(seeds,arr);
6919 SignClusters(seeds,fnumber,fdensity);
6924 arr = Tracking(3,nup-delta-5,nup-delta-5-gap,cuts,-1);
6925 SumTracks(seeds,arr);
6926 SignClusters(seeds,fnumber,fdensity);
6937 Info("Tracking()","\n\nPrimary seeding\t%d\n\n",seeds->GetEntriesFast());
6941 // RemoveUsed(seeds,0.75,0.75,1);
6943 //SignClusters(seeds,fnumber,fdensity);
6952 arr = Tracking(4,nup-1,nup-1-gap,cuts,-1);
6953 SumTracks(seeds,arr);
6954 SignClusters(seeds,fnumber,fdensity);
6956 arr = Tracking(4,nup-2,nup-2-gap,cuts,-1);
6957 SumTracks(seeds,arr);
6958 SignClusters(seeds,fnumber,fdensity);
6960 arr = Tracking(4,nup-3,nup-3-gap,cuts,-1);
6961 SumTracks(seeds,arr);
6962 SignClusters(seeds,fnumber,fdensity);
6964 arr = Tracking(4,nup-5,nup-5-gap,cuts,-1);
6965 SumTracks(seeds,arr);
6966 SignClusters(seeds,fnumber,fdensity);
6968 arr = Tracking(4,nup-7,nup-7-gap,cuts,-1);
6969 SumTracks(seeds,arr);
6970 SignClusters(seeds,fnumber,fdensity);
6973 arr = Tracking(4,nup-9,nup-9-gap,cuts,-1);
6974 SumTracks(seeds,arr);
6975 SignClusters(seeds,fnumber,fdensity);
6979 for (Int_t delta = 9; delta<30; delta+=gapSec){
6985 arr = Tracking(4,nup-1-delta,nup-1-delta-gap,cuts,-1);
6986 SumTracks(seeds,arr);
6987 SignClusters(seeds,fnumber,fdensity);
6989 arr = Tracking(4,nup-3-delta,nup-5-delta-gap,cuts,4);
6990 SumTracks(seeds,arr);
6991 SignClusters(seeds,fnumber,fdensity);
7004 for (Int_t delta = 30; delta<fLastSeedRowSec; delta+=gapSec){
7010 arr = Tracking(4,nup-1-delta,nup-1-delta-gap,cuts,-1);
7011 SumTracks(seeds,arr);
7012 SignClusters(seeds,fnumber,fdensity);
7014 arr = Tracking(4,nup-5-delta,nup-5-delta-gap,cuts,5 );
7015 SumTracks(seeds,arr);
7016 SignClusters(seeds,fnumber,fdensity);
7020 Info("Tracking()","\n\nSecondary seeding\t%d\n\n",seeds->GetEntriesFast());
7031 TObjArray * AliTPCtrackerMI::TrackingSpecial()
7034 // seeding adjusted for laser and cosmic tests - short tracks with big inclination angle
7035 // no primary vertex seeding tried
7039 Int_t nup=fOuterSec->GetNRows()+fInnerSec->GetNRows();
7041 TObjArray * seeds = new TObjArray;
7046 Float_t fnumber = 3.0;
7047 Float_t fdensity = 3.0;
7050 cuts[0] = AliTPCReconstructor::GetRecoParam()->GetMaxC(); // max curvature
7051 cuts[1] = 3.5; // max tan(phi) angle for seeding
7052 cuts[2] = 3.; // not used (cut on z primary vertex)
7053 cuts[3] = 3.5; // max tan(theta) angle for seeding
7055 for (Int_t delta = 0; nup-delta-gap-1>0; delta+=3){
7057 arr = Tracking(4,nup-1-delta,nup-1-delta-gap,cuts,-1);
7058 SumTracks(seeds,arr);
7059 SignClusters(seeds,fnumber,fdensity);
7063 Info("Tracking()","\n\nSecondary seeding\t%d\n\n",seeds->GetEntriesFast());
7074 void AliTPCtrackerMI::SumTracks(TObjArray *arr1,TObjArray *&arr2)
7077 //sum tracks to common container
7078 //remove suspicious tracks
7079 // RS: Attention: supplied tracks come in the static array, don't delete them
7080 Int_t nseed = arr2->GetEntriesFast();
7081 for (Int_t i=0;i<nseed;i++){
7082 AliTPCseed *pt=(AliTPCseed*)arr2->UncheckedAt(i);
7085 // remove tracks with too big curvature
7087 if (TMath::Abs(pt->GetC())>AliTPCReconstructor::GetRecoParam()->GetMaxC()){
7088 MarkSeedFree( arr2->RemoveAt(i) );
7091 // REMOVE VERY SHORT TRACKS
7092 if (pt->GetNumberOfClusters()<20){
7093 MarkSeedFree( arr2->RemoveAt(i) );
7096 // NORMAL ACTIVE TRACK
7097 if (pt->IsActive()){
7098 arr1->AddLast(arr2->RemoveAt(i));
7101 //remove not usable tracks
7102 if (pt->GetRemoval()!=10){
7103 MarkSeedFree( arr2->RemoveAt(i) );
7107 // ENABLE ONLY ENOUGH GOOD STOPPED TRACKS
7108 if (pt->GetDensityFirst(20)>0.8 || pt->GetDensityFirst(30)>0.8 || pt->GetDensityFirst(40)>0.7)
7109 arr1->AddLast(arr2->RemoveAt(i));
7111 MarkSeedFree( arr2->RemoveAt(i) );
7115 // delete arr2; arr2 = 0; // RS: this is static array, don't delete it
7120 void AliTPCtrackerMI::ParallelTracking(TObjArray *const arr, Int_t rfirst, Int_t rlast)
7123 // try to track in parralel
7125 Int_t nseed=arr->GetEntriesFast();
7126 //prepare seeds for tracking
7127 for (Int_t i=0; i<nseed; i++) {
7128 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i), &t=*pt;
7130 if (!t.IsActive()) continue;
7131 // follow prolongation to the first layer
7132 if ( (fSectors ==fInnerSec) || (t.GetFirstPoint()-fkParam->GetNRowLow()>rfirst+1) )
7133 FollowProlongation(t, rfirst+1);
7138 for (Int_t nr=rfirst; nr>=rlast; nr--){
7139 if (nr<fInnerSec->GetNRows())
7140 fSectors = fInnerSec;
7142 fSectors = fOuterSec;
7143 // make indexes with the cluster tracks for given
7145 // find nearest cluster
7146 for (Int_t i=0; i<nseed; i++) {
7147 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i), &t=*pt;
7149 if (nr==80) pt->UpdateReference();
7150 if (!pt->IsActive()) continue;
7151 // if ( (fSectors ==fOuterSec) && (pt->fFirstPoint-fkParam->GetNRowLow())<nr) continue;
7152 if (pt->GetRelativeSector()>17) {
7155 UpdateClusters(t,nr);
7157 // prolonagate to the nearest cluster - if founded
7158 for (Int_t i=0; i<nseed; i++) {
7159 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
7161 if (!pt->IsActive()) continue;
7162 // if ((fSectors ==fOuterSec) && (pt->fFirstPoint-fkParam->GetNRowLow())<nr) continue;
7163 if (pt->GetRelativeSector()>17) {
7166 FollowToNextCluster(*pt,nr);
7171 void AliTPCtrackerMI::PrepareForBackProlongation(const TObjArray *const arr,Float_t fac) const
7175 // if we use TPC track itself we have to "update" covariance
7177 Int_t nseed= arr->GetEntriesFast();
7178 for (Int_t i=0;i<nseed;i++){
7179 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
7183 //rotate to current local system at first accepted point
7184 Int_t index = pt->GetClusterIndex2(pt->GetFirstPoint());
7185 Int_t sec = (index&0xff000000)>>24;
7187 Float_t angle1 = fInnerSec->GetAlpha()*sec+fInnerSec->GetAlphaShift();
7188 if (angle1>TMath::Pi())
7189 angle1-=2.*TMath::Pi();
7190 Float_t angle2 = pt->GetAlpha();
7192 if (TMath::Abs(angle1-angle2)>0.001){
7193 if (!pt->Rotate(angle1-angle2)) return;
7194 //angle2 = pt->GetAlpha();
7195 //pt->fRelativeSector = pt->GetAlpha()/fInnerSec->GetAlpha();
7196 //if (pt->GetAlpha()<0)
7197 // pt->fRelativeSector+=18;
7198 //sec = pt->fRelativeSector;
7207 void AliTPCtrackerMI::PrepareForProlongation(TObjArray *const arr, Float_t fac) const
7211 // if we use TPC track itself we have to "update" covariance
7213 Int_t nseed= arr->GetEntriesFast();
7214 for (Int_t i=0;i<nseed;i++){
7215 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
7218 pt->SetFirstPoint(pt->GetLastPoint());
7226 Int_t AliTPCtrackerMI::PropagateBack(const TObjArray *const arr)
7229 // make back propagation
7231 Int_t nseed= arr->GetEntriesFast();
7232 for (Int_t i=0;i<nseed;i++){
7233 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
7234 if (pt&& pt->GetKinkIndex(0)<=0) {
7235 //AliTPCseed *pt2 = new AliTPCseed(*pt);
7236 fSectors = fInnerSec;
7237 //FollowBackProlongation(*pt,fInnerSec->GetNRows()-1);
7238 //fSectors = fOuterSec;
7239 FollowBackProlongation(*pt,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1,1);
7240 //if (pt->GetNumberOfClusters()<(pt->fEsd->GetTPCclusters(0)) ){
7241 // Error("PropagateBack","Not prolonged track %d",pt->GetLabel());
7242 // FollowBackProlongation(*pt2,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1);
7245 if (pt&& pt->GetKinkIndex(0)>0) {
7246 AliESDkink * kink = fEvent->GetKink(pt->GetKinkIndex(0)-1);
7247 pt->SetFirstPoint(kink->GetTPCRow0());
7248 fSectors = fInnerSec;
7249 FollowBackProlongation(*pt,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1,1);
7257 Int_t AliTPCtrackerMI::PropagateForward2(const TObjArray *const arr)
7260 // make forward propagation
7262 Int_t nseed= arr->GetEntriesFast();
7264 for (Int_t i=0;i<nseed;i++){
7265 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
7267 FollowProlongation(*pt,0,1,1);
7276 Int_t AliTPCtrackerMI::PropagateForward()
7279 // propagate track forward
7281 Int_t nseed = fSeeds->GetEntriesFast();
7282 for (Int_t i=0;i<nseed;i++){
7283 AliTPCseed *pt = (AliTPCseed*)fSeeds->UncheckedAt(i);
7285 AliTPCseed &t = *pt;
7286 Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
7287 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
7288 if (alpha < 0. ) alpha += 2.*TMath::Pi();
7289 t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN);
7293 fSectors = fOuterSec;
7294 ParallelTracking(fSeeds,fOuterSec->GetNRows()+fInnerSec->GetNRows()-1,fInnerSec->GetNRows());
7295 fSectors = fInnerSec;
7296 ParallelTracking(fSeeds,fInnerSec->GetNRows()-1,0);
7305 Int_t AliTPCtrackerMI::PropagateBack(AliTPCseed *const pt, Int_t row0, Int_t row1)
7308 // make back propagation, in between row0 and row1
7312 fSectors = fInnerSec;
7315 if (row1<fSectors->GetNRows())
7318 r1 = fSectors->GetNRows()-1;
7320 if (row0<fSectors->GetNRows()&& r1>0 )
7321 FollowBackProlongation(*pt,r1);
7322 if (row1<=fSectors->GetNRows())
7325 r1 = row1 - fSectors->GetNRows();
7326 if (r1<=0) return 0;
7327 if (r1>=fOuterSec->GetNRows()) return 0;
7328 fSectors = fOuterSec;
7329 return FollowBackProlongation(*pt,r1);
7337 void AliTPCtrackerMI::GetShape(AliTPCseed * seed, Int_t row)
7339 // gets cluster shape
7341 AliTPCClusterParam * clparam = AliTPCcalibDB::Instance()->GetClusterParam();
7342 Float_t zdrift = TMath::Abs((fkParam->GetZLength(0)-TMath::Abs(seed->GetZ())));
7343 Int_t type = (seed->GetSector() < fkParam->GetNSector()/2) ? 0: (row>126) ? 1:2;
7344 Double_t angulary = seed->GetSnp();
7346 if (TMath::Abs(angulary)>AliTPCReconstructor::GetMaxSnpTracker()) {
7347 angulary = TMath::Sign(AliTPCReconstructor::GetMaxSnpTracker(),angulary);
7350 angulary = angulary*angulary/((1.-angulary)*(1.+angulary));
7351 Double_t angularz = seed->GetTgl()*seed->GetTgl()*(1.+angulary);
7353 Double_t sigmay = clparam->GetRMS0(0,type,zdrift,TMath::Sqrt(TMath::Abs(angulary)));
7354 Double_t sigmaz = clparam->GetRMS0(1,type,zdrift,TMath::Sqrt(TMath::Abs(angularz)));
7355 seed->SetCurrentSigmaY2(sigmay*sigmay);
7356 seed->SetCurrentSigmaZ2(sigmaz*sigmaz);
7357 // Float_t sd2 = TMath::Abs((fkParam->GetZLength(0)-TMath::Abs(seed->GetZ())))*fkParam->GetDiffL()*fkParam->GetDiffL();
7358 // // Float_t padlength = fkParam->GetPadPitchLength(seed->fSector);
7359 // Float_t padlength = GetPadPitchLength(row);
7361 // Float_t sresy = (seed->GetSector() < fkParam->GetNSector()/2) ? 0.2 :0.3;
7362 // seed->SetCurrentSigmaY2(sd2+padlength*padlength*angulary/12.+sresy*sresy);
7364 // Float_t sresz = fkParam->GetZSigma();
7365 // seed->SetCurrentSigmaZ2(sd2+padlength*padlength*angularz*angularz*(1+angulary)/12.+sresz*sresz);
7367 Float_t wy = GetSigmaY(seed);
7368 Float_t wz = GetSigmaZ(seed);
7371 if (TMath::Abs(wy/seed->fCurrentSigmaY2-1)>0.0001 || TMath::Abs(wz/seed->fCurrentSigmaZ2-1)>0.0001 ){
7372 printf("problem\n");
7379 //__________________________________________________________________________
7380 void AliTPCtrackerMI::CookLabel(AliKalmanTrack *tk, Float_t wrong) const {
7381 //--------------------------------------------------------------------
7382 //This function "cooks" a track label. If label<0, this track is fake.
7383 //--------------------------------------------------------------------
7384 AliTPCseed * t = dynamic_cast<AliTPCseed*>(tk);
7386 printf("%s:%d wrong type \n",(char*)__FILE__,__LINE__);
7390 Int_t noc=t->GetNumberOfClusters();
7392 //printf("\nnot founded prolongation\n\n\n");
7398 AliTPCclusterMI *clusters[160];
7400 for (Int_t i=0;i<160;i++) {
7407 for (i=0; i<160 && current<noc; i++) {
7409 Int_t index=t->GetClusterIndex2(i);
7410 if (index<=0) continue;
7411 if (index&0x8000) continue;
7413 //clusters[current]=GetClusterMI(index);
7414 if (t->GetClusterPointer(i)){
7415 clusters[current]=t->GetClusterPointer(i);
7421 Int_t lab=123456789;
7422 for (i=0; i<noc; i++) {
7423 AliTPCclusterMI *c=clusters[i];
7425 lab=TMath::Abs(c->GetLabel(0));
7427 for (j=0; j<noc; j++) if (lb[j]==lab || mx[j]==0) break;
7433 for (i=0; i<noc; i++) if (mx[i]>max) {max=mx[i]; lab=lb[i];}
7435 for (i=0; i<noc; i++) {
7436 AliTPCclusterMI *c=clusters[i];
7438 if (TMath::Abs(c->GetLabel(1)) == lab ||
7439 TMath::Abs(c->GetLabel(2)) == lab ) max++;
7441 if (noc<=0) { lab=-1; return;}
7442 if ((1.- Float_t(max)/(noc)) > wrong) lab=-lab;
7445 Int_t tail=Int_t(0.10*noc);
7448 for (i=1; i<160&&ind<tail; i++) {
7449 // AliTPCclusterMI *c=clusters[noc-i];
7450 AliTPCclusterMI *c=clusters[i];
7452 if (lab == TMath::Abs(c->GetLabel(0)) ||
7453 lab == TMath::Abs(c->GetLabel(1)) ||
7454 lab == TMath::Abs(c->GetLabel(2))) max++;
7457 if (max < Int_t(0.5*tail)) lab=-lab;
7464 //delete[] clusters;
7468 //__________________________________________________________________________
7469 Int_t AliTPCtrackerMI::CookLabel(AliTPCseed *const t, Float_t wrong,Int_t first, Int_t last) const {
7470 //--------------------------------------------------------------------
7471 //This function "cooks" a track label. If label<0, this track is fake.
7472 //--------------------------------------------------------------------
7473 Int_t noc=t->GetNumberOfClusters();
7475 //printf("\nnot founded prolongation\n\n\n");
7481 AliTPCclusterMI *clusters[160];
7483 for (Int_t i=0;i<160;i++) {
7490 for (i=0; i<160 && current<noc; i++) {
7491 if (i<first) continue;
7492 if (i>last) continue;
7493 Int_t index=t->GetClusterIndex2(i);
7494 if (index<=0) continue;
7495 if (index&0x8000) continue;
7497 //clusters[current]=GetClusterMI(index);
7498 if (t->GetClusterPointer(i)){
7499 clusters[current]=t->GetClusterPointer(i);
7504 //if (noc<5) return -1;
7505 Int_t lab=123456789;
7506 for (i=0; i<noc; i++) {
7507 AliTPCclusterMI *c=clusters[i];
7509 lab=TMath::Abs(c->GetLabel(0));
7511 for (j=0; j<noc; j++) if (lb[j]==lab || mx[j]==0) break;
7517 for (i=0; i<noc; i++) if (mx[i]>max) {max=mx[i]; lab=lb[i];}
7519 for (i=0; i<noc; i++) {
7520 AliTPCclusterMI *c=clusters[i];
7522 if (TMath::Abs(c->GetLabel(1)) == lab ||
7523 TMath::Abs(c->GetLabel(2)) == lab ) max++;
7525 if (noc<=0) { lab=-1; return -1;}
7526 if ((1.- Float_t(max)/(noc)) > wrong) lab=-lab;
7529 Int_t tail=Int_t(0.10*noc);
7532 for (i=1; i<160&&ind<tail; i++) {
7533 // AliTPCclusterMI *c=clusters[noc-i];
7534 AliTPCclusterMI *c=clusters[i];
7536 if (lab == TMath::Abs(c->GetLabel(0)) ||
7537 lab == TMath::Abs(c->GetLabel(1)) ||
7538 lab == TMath::Abs(c->GetLabel(2))) max++;
7541 if (max < Int_t(0.5*tail)) lab=-lab;
7544 // t->SetLabel(lab);
7548 //delete[] clusters;
7552 Int_t AliTPCtrackerMI::GetRowNumber(Double_t x[3]) const
7554 //return pad row number for given x vector
7555 Float_t phi = TMath::ATan2(x[1],x[0]);
7556 if(phi<0) phi=2.*TMath::Pi()+phi;
7557 // Get the local angle in the sector philoc
7558 const Float_t kRaddeg = 180/3.14159265358979312;
7559 Float_t phiangle = (Int_t (phi*kRaddeg/20.) + 0.5)*20./kRaddeg;
7560 Double_t localx = x[0]*TMath::Cos(phiangle)-x[1]*TMath::Sin(phiangle);
7561 return GetRowNumber(localx);
7566 void AliTPCtrackerMI::MakeESDBitmaps(AliTPCseed *t, AliESDtrack *esd)
7568 //-----------------------------------------------------------------------
7569 // Fill the cluster and sharing bitmaps of the track
7570 //-----------------------------------------------------------------------
7572 Int_t firstpoint = 0;
7573 Int_t lastpoint = 159;
7574 AliTPCTrackerPoint *point;
7575 AliTPCclusterMI *cluster;
7578 TBits clusterMap(159);
7579 TBits sharedMap(159);
7581 for (int iter=firstpoint; iter<lastpoint; iter++) {
7582 // Change to cluster pointers to see if we have a cluster at given padrow
7584 cluster = t->GetClusterPointer(iter);
7586 clusterMap.SetBitNumber(iter, kTRUE);
7587 point = t->GetTrackPoint(iter);
7588 if (point->IsShared())
7589 sharedMap.SetBitNumber(iter,kTRUE);
7591 if (t->GetClusterIndex(iter) >= 0 && (t->GetClusterIndex(iter) & 0x8000) == 0) {
7592 fitMap.SetBitNumber(iter, kTRUE);
7596 esd->SetTPCClusterMap(clusterMap);
7597 esd->SetTPCSharedMap(sharedMap);
7598 esd->SetTPCFitMap(fitMap);
7599 if (nclsf != t->GetNumberOfClusters())
7600 AliWarning(Form("Inconsistency between ncls %d and indices %d (found %d)",t->GetNumberOfClusters(),nclsf,esd->GetTPCClusterMap().CountBits()));
7603 Bool_t AliTPCtrackerMI::IsFindable(AliTPCseed & track){
7605 // return flag if there is findable cluster at given position
7608 Float_t z = track.GetZ();
7610 if (TMath::Abs(z)<(AliTPCReconstructor::GetCtgRange()*track.GetX()+kDeltaZ) &&
7611 TMath::Abs(z)<fkParam->GetZLength(0) &&
7612 (TMath::Abs(track.GetSnp())<AliTPCReconstructor::GetMaxSnpTracker()))
7618 void AliTPCtrackerMI::AddCovariance(AliTPCseed * seed){
7620 // Adding systematic error estimate to the covariance matrix
7621 // !!!! the systematic error for element 4 is in 1/cm not in pt
7623 const Double_t *param = AliTPCReconstructor::GetRecoParam()->GetSystematicError();
7625 // use only the diagonal part if not specified otherwise
7626 if (!AliTPCReconstructor::GetRecoParam()->GetUseSystematicCorrelation()) return AddCovarianceAdd(seed);
7628 Double_t *covarS= (Double_t*)seed->GetCovariance();
7629 Double_t factor[5]={1,1,1,1,1};
7630 Double_t facC = AliTracker::GetBz()*kB2C;
7631 factor[0]= TMath::Sqrt(TMath::Abs((covarS[0] + param[0]*param[0])/covarS[0]));
7632 factor[1]= TMath::Sqrt(TMath::Abs((covarS[2] + param[1]*param[1])/covarS[2]));
7633 factor[2]= TMath::Sqrt(TMath::Abs((covarS[5] + param[2]*param[2])/covarS[5]));
7634 factor[3]= TMath::Sqrt(TMath::Abs((covarS[9] + param[3]*param[3])/covarS[9]));
7635 factor[4]= TMath::Sqrt(TMath::Abs((covarS[14] + facC*facC*param[4]*param[4])/covarS[14]));
7641 for (Int_t i=0; i<5; i++){
7642 for (Int_t j=i; j<5; j++){
7643 Int_t index=seed->GetIndex(i,j);
7644 covarS[index]*=factor[i]*factor[j];
7650 void AliTPCtrackerMI::AddCovarianceAdd(AliTPCseed * seed){
7652 // Adding systematic error - as additive factor without correlation
7654 const Double_t *param = AliTPCReconstructor::GetRecoParam()->GetSystematicError();
7655 Double_t *covarIn= (Double_t*)seed->GetCovariance();
7657 for (Int_t i=0;i<15;i++) covar[i]=0;
7663 covar[0] = param[0]*param[0];
7664 covar[2] = param[1]*param[1];
7665 covar[5] = param[2]*param[2];
7666 covar[9] = param[3]*param[3];
7667 Double_t facC = AliTracker::GetBz()*kB2C;
7668 covar[14]= param[4]*param[4]*facC*facC;
7670 covar[1]=TMath::Sqrt((covar[0]*covar[2]))*covarIn[1]/TMath::Sqrt((covarIn[0]*covarIn[2]));
7672 covar[3]=TMath::Sqrt((covar[0]*covar[5]))*covarIn[3]/TMath::Sqrt((covarIn[0]*covarIn[5]));
7673 covar[4]=TMath::Sqrt((covar[2]*covar[5]))*covarIn[4]/TMath::Sqrt((covarIn[2]*covarIn[5]));
7675 covar[6]=TMath::Sqrt((covar[0]*covar[9]))*covarIn[6]/TMath::Sqrt((covarIn[0]*covarIn[9]));
7676 covar[7]=TMath::Sqrt((covar[2]*covar[9]))*covarIn[7]/TMath::Sqrt((covarIn[2]*covarIn[9]));
7677 covar[8]=TMath::Sqrt((covar[5]*covar[9]))*covarIn[8]/TMath::Sqrt((covarIn[5]*covarIn[9]));
7679 covar[10]=TMath::Sqrt((covar[0]*covar[14]))*covarIn[10]/TMath::Sqrt((covarIn[0]*covarIn[14]));
7680 covar[11]=TMath::Sqrt((covar[2]*covar[14]))*covarIn[11]/TMath::Sqrt((covarIn[2]*covarIn[14]));
7681 covar[12]=TMath::Sqrt((covar[5]*covar[14]))*covarIn[12]/TMath::Sqrt((covarIn[5]*covarIn[14]));
7682 covar[13]=TMath::Sqrt((covar[9]*covar[14]))*covarIn[13]/TMath::Sqrt((covarIn[9]*covarIn[14]));
7684 seed->AddCovariance(covar);
7687 //_____________________________________________________________________________
7688 Bool_t AliTPCtrackerMI::IsTPCHVDipEvent(AliESDEvent const *esdEvent) {
7690 // check events affected by TPC HV dip
7692 if(!esdEvent) return kFALSE;
7695 if(!AliTPCcalibDB::Instance()) return kFALSE;
7696 AliTPCcalibDB::Instance()->SetRun(esdEvent->GetRunNumber());
7698 // Get HV TPC chamber sensors and calculate the median
7699 AliDCSSensorArray *voltageArray= AliTPCcalibDB::Instance()->GetVoltageSensors(esdEvent->GetRunNumber());
7700 if(!voltageArray) return kFALSE;
7702 TString sensorName="";
7703 Double_t kTPCHVdip = 2.0; // allow for 2V dip as compared to median from given sensor
7706 for(Int_t sector=0; sector<72; sector++)
7708 Char_t sideName='A';
7709 if ((sector/18)%2==1) sideName='C';
7712 sensorName=Form("TPC_ANODE_I_%c%02d_VMEAS",sideName,sector%18);
7715 sensorName=Form("TPC_ANODE_O_%c%02d_0_VMEAS",sideName,sector%18);
7718 AliDCSSensor* sensor = voltageArray->GetSensor(sensorName.Data());
7719 if(!sensor) continue;
7720 TGraph *graph = sensor->GetGraph();
7721 if(!graph) continue;
7722 Double_t median = TMath::Median(graph->GetN(), graph->GetY());
7723 if(median == 0) continue;
7725 //printf("chamber %d, sensor %s, HV %f, median %f\n", sector, sensorName.Data(), sensor->GetValue(esdEvent->GetTimeStamp()), median);
7727 if(TMath::Abs(sensor->GetValue(esdEvent->GetTimeStamp())-median)>kTPCHVdip) {
7735 //________________________________________
7736 void AliTPCtrackerMI::MarkSeedFree(TObject *sd)
7738 // account that this seed is "deleted"
7739 AliTPCseed* seed = dynamic_cast<AliTPCseed*>(sd);
7740 if (!sd) {AliError(Form("Freeing of non-AliTPCseed %p from the pool is requested",sd)); return;}
7741 int id = seed->GetPoolID();
7742 if (id<0) {AliError(Form("Freeing of seed %p NOT from the pool is requested",sd)); return;}
7743 // AliInfo(Form("%d %p",id, seed));
7744 fSeedsPool->RemoveAt(id);
7745 if (fFreeSeedsID.GetSize()<=fNFreeSeeds) fFreeSeedsID.Set( 2*fNFreeSeeds + 100 );
7746 fFreeSeedsID.GetArray()[fNFreeSeeds++] = id;
7749 //________________________________________
7750 TObject *&AliTPCtrackerMI::NextFreeSeed()
7752 // return next free slot where the seed can be created
7753 fLastSeedID = fNFreeSeeds ? fFreeSeedsID.GetArray()[--fNFreeSeeds] : fSeedsPool->GetEntriesFast();
7754 // AliInfo(Form("%d",fLastSeedID));
7755 return (*fSeedsPool)[ fLastSeedID ];
7759 //________________________________________
7760 void AliTPCtrackerMI::ResetSeedsPool()
7762 // mark all seeds in the pool as unused
7763 AliInfo(Form("CurrentSize: %d, BookedUpTo: %d, free: %d",fSeedsPool->GetSize(),fSeedsPool->GetEntriesFast(),fNFreeSeeds));
7765 fSeedsPool->Clear("C"); // RS: nominally the seeds may allocate memory...