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"
148 ClassImp(AliTPCtrackerMI)
152 class AliTPCFastMath {
155 static Double_t FastAsin(Double_t x);
157 static Double_t fgFastAsin[20000]; //lookup table for fast asin computation
160 Double_t AliTPCFastMath::fgFastAsin[20000];
161 AliTPCFastMath gAliTPCFastMath; // needed to fill the LUT
163 AliTPCFastMath::AliTPCFastMath(){
165 // initialized lookup table;
166 for (Int_t i=0;i<10000;i++){
167 fgFastAsin[2*i] = TMath::ASin(i/10000.);
168 fgFastAsin[2*i+1] = (TMath::ASin((i+1)/10000.)-fgFastAsin[2*i]);
172 Double_t AliTPCFastMath::FastAsin(Double_t x){
174 // return asin using lookup table
176 Int_t index = int(x*10000);
177 return fgFastAsin[2*index]+(x*10000.-index)*fgFastAsin[2*index+1];
180 Int_t index = int(x*10000);
181 return -(fgFastAsin[2*index]+(x*10000.-index)*fgFastAsin[2*index+1]);
183 //__________________________________________________________________
184 AliTPCtrackerMI::AliTPCtrackerMI()
211 // default constructor
213 for (Int_t irow=0; irow<200; irow++){
220 //_____________________________________________________________________
224 Int_t AliTPCtrackerMI::UpdateTrack(AliTPCseed * track, Int_t accept){
226 //update track information using current cluster - track->fCurrentCluster
229 AliTPCclusterMI* c =track->GetCurrentCluster();
230 if (accept > 0) //sign not accepted clusters
231 track->SetCurrentClusterIndex1(track->GetCurrentClusterIndex1() | 0x8000);
232 else // unsign accpeted clusters
233 track->SetCurrentClusterIndex1(track->GetCurrentClusterIndex1() & 0xffff7fff);
234 UInt_t i = track->GetCurrentClusterIndex1();
236 Int_t sec=(i&0xff000000)>>24;
237 //Int_t row = (i&0x00ff0000)>>16;
238 track->SetRow((i&0x00ff0000)>>16);
239 track->SetSector(sec);
240 // Int_t index = i&0xFFFF;
241 if (sec>=fkParam->GetNInnerSector()) track->SetRow(track->GetRow()+fkParam->GetNRowLow());
242 track->SetClusterIndex2(track->GetRow(), i);
243 //track->fFirstPoint = row;
244 //if ( track->fLastPoint<row) track->fLastPoint =row;
245 // if (track->fRow<0 || track->fRow>160) {
246 // printf("problem\n");
248 if (track->GetFirstPoint()>track->GetRow())
249 track->SetFirstPoint(track->GetRow());
250 if (track->GetLastPoint()<track->GetRow())
251 track->SetLastPoint(track->GetRow());
254 track->SetClusterPointer(track->GetRow(),c);
257 Double_t angle2 = track->GetSnp()*track->GetSnp();
259 //SET NEW Track Point
261 if (angle2<1) //PH sometimes angle2 is very big. To be investigated...
263 angle2 = TMath::Sqrt(angle2/(1-angle2));
264 AliTPCTrackerPoint &point =*(track->GetTrackPoint(track->GetRow()));
266 point.SetSigmaY(c->GetSigmaY2()/track->GetCurrentSigmaY2());
267 point.SetSigmaZ(c->GetSigmaZ2()/track->GetCurrentSigmaZ2());
268 point.SetErrY(sqrt(track->GetErrorY2()));
269 point.SetErrZ(sqrt(track->GetErrorZ2()));
271 point.SetX(track->GetX());
272 point.SetY(track->GetY());
273 point.SetZ(track->GetZ());
274 point.SetAngleY(angle2);
275 point.SetAngleZ(track->GetTgl());
276 if (point.IsShared()){
277 track->SetErrorY2(track->GetErrorY2()*4);
278 track->SetErrorZ2(track->GetErrorZ2()*4);
282 Double_t chi2 = track->GetPredictedChi2(track->GetCurrentCluster());
284 // track->SetErrorY2(track->GetErrorY2()*1.3);
285 // track->SetErrorY2(track->GetErrorY2()+0.01);
286 // track->SetErrorZ2(track->GetErrorZ2()*1.3);
287 // track->SetErrorZ2(track->GetErrorZ2()+0.005);
289 if (accept>0) return 0;
290 if (track->GetNumberOfClusters()%20==0){
291 // if (track->fHelixIn){
292 // TClonesArray & larr = *(track->fHelixIn);
293 // Int_t ihelix = larr.GetEntriesFast();
294 // new(larr[ihelix]) AliHelix(*track) ;
297 track->SetNoCluster(0);
298 return track->Update(c,chi2,i);
303 Int_t AliTPCtrackerMI::AcceptCluster(AliTPCseed * seed, AliTPCclusterMI * cluster)
306 // decide according desired precision to accept given
307 // cluster for tracking
309 seed->GetProlongation(cluster->GetX(),yt,zt);
310 Double_t sy2=ErrY2(seed,cluster);
311 Double_t sz2=ErrZ2(seed,cluster);
313 Double_t sdistancey2 = sy2+seed->GetSigmaY2();
314 Double_t sdistancez2 = sz2+seed->GetSigmaZ2();
315 Double_t dy=seed->GetCurrentCluster()->GetY()-yt;
316 Double_t dz=seed->GetCurrentCluster()->GetZ()-zt;
317 Double_t rdistancey2 = (seed->GetCurrentCluster()->GetY()-yt)*
318 (seed->GetCurrentCluster()->GetY()-yt)/sdistancey2;
319 Double_t rdistancez2 = (seed->GetCurrentCluster()->GetZ()-zt)*
320 (seed->GetCurrentCluster()->GetZ()-zt)/sdistancez2;
322 Double_t rdistance2 = rdistancey2+rdistancez2;
325 if (AliTPCReconstructor::StreamLevel()>2 && seed->GetNumberOfClusters()>20) {
326 Float_t rmsy2 = seed->GetCurrentSigmaY2();
327 Float_t rmsz2 = seed->GetCurrentSigmaZ2();
328 Float_t rmsy2p30 = seed->GetCMeanSigmaY2p30();
329 Float_t rmsz2p30 = seed->GetCMeanSigmaZ2p30();
330 Float_t rmsy2p30R = seed->GetCMeanSigmaY2p30R();
331 Float_t rmsz2p30R = seed->GetCMeanSigmaZ2p30R();
332 AliExternalTrackParam param(*seed);
333 static TVectorD gcl(3),gtr(3);
335 param.GetXYZ(gcl.GetMatrixArray());
336 cluster->GetGlobalXYZ(gclf);
337 gcl[0]=gclf[0]; gcl[1]=gclf[1]; gcl[2]=gclf[2];
340 if (AliTPCReconstructor::StreamLevel()>2) {
341 (*fDebugStreamer)<<"ErrParam"<<
354 "rmsy2p30="<<rmsy2p30<<
355 "rmsz2p30="<<rmsz2p30<<
356 "rmsy2p30R="<<rmsy2p30R<<
357 "rmsz2p30R="<<rmsz2p30R<<
358 // normalize distance -
359 "rdisty="<<rdistancey2<<
360 "rdistz="<<rdistancez2<<
361 "rdist="<<rdistance2<< //
365 //return 0; // temporary
366 if (rdistance2>32) return 3;
369 if ((rdistancey2>9. || rdistancez2>9.) && cluster->GetType()==0)
370 return 2; //suspisiouce - will be changed
372 if ((rdistancey2>6.25 || rdistancez2>6.25) && cluster->GetType()>0)
373 // strict cut on overlaped cluster
374 return 2; //suspisiouce - will be changed
376 if ( (rdistancey2>1. || rdistancez2>6.25 )
377 && cluster->GetType()<0){
378 seed->SetNFoundable(seed->GetNFoundable()-1);
382 if (fUseHLTClusters == 3 || fUseHLTClusters == 4) {
384 if(!AliTPCReconstructor::GetRecoParam()->GetUseHLTOnePadCluster()) {
385 if (TMath::Abs(cluster->GetSigmaY2()) < kAlmost0)
398 //_____________________________________________________________________________
399 AliTPCtrackerMI::AliTPCtrackerMI(const AliTPCParam *par):
401 fkNIS(par->GetNInnerSector()/2),
403 fkNOS(par->GetNOuterSector()/2),
425 //---------------------------------------------------------------------
426 // The main TPC tracker constructor
427 //---------------------------------------------------------------------
428 fInnerSec=new AliTPCtrackerSector[fkNIS];
429 fOuterSec=new AliTPCtrackerSector[fkNOS];
432 for (i=0; i<fkNIS; i++) fInnerSec[i].Setup(par,0);
433 for (i=0; i<fkNOS; i++) fOuterSec[i].Setup(par,1);
436 Int_t nrowlow = par->GetNRowLow();
437 Int_t nrowup = par->GetNRowUp();
440 for (i=0;i<nrowlow;i++){
441 fXRow[i] = par->GetPadRowRadiiLow(i);
442 fPadLength[i]= par->GetPadPitchLength(0,i);
443 fYMax[i] = fXRow[i]*TMath::Tan(0.5*par->GetInnerAngle());
447 for (i=0;i<nrowup;i++){
448 fXRow[i+nrowlow] = par->GetPadRowRadiiUp(i);
449 fPadLength[i+nrowlow] = par->GetPadPitchLength(60,i);
450 fYMax[i+nrowlow] = fXRow[i+nrowlow]*TMath::Tan(0.5*par->GetOuterAngle());
453 if (AliTPCReconstructor::StreamLevel()>0) {
454 fDebugStreamer = new TTreeSRedirector("TPCdebug.root");
457 fSeedsPool = new TClonesArray("AliTPCseed",1000);
459 //________________________________________________________________________
460 AliTPCtrackerMI::AliTPCtrackerMI(const AliTPCtrackerMI &t):
486 //------------------------------------
487 // dummy copy constructor
488 //------------------------------------------------------------------
490 for (Int_t irow=0; irow<200; irow++){
497 AliTPCtrackerMI & AliTPCtrackerMI::operator=(const AliTPCtrackerMI& /*r*/)
499 //------------------------------
501 //--------------------------------------------------------------
504 //_____________________________________________________________________________
505 AliTPCtrackerMI::~AliTPCtrackerMI() {
506 //------------------------------------------------------------------
507 // TPC tracker destructor
508 //------------------------------------------------------------------
515 if (fDebugStreamer) delete fDebugStreamer;
516 if (fSeedsPool) delete fSeedsPool;
520 void AliTPCtrackerMI::FillESD(const TObjArray* arr)
524 //fill esds using updated tracks
527 // write tracks to the event
528 // store index of the track
529 Int_t nseed=arr->GetEntriesFast();
530 //FindKinks(arr,fEvent);
531 for (Int_t i=0; i<nseed; i++) {
532 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
536 if (AliTPCReconstructor::StreamLevel()>1) {
537 (*fDebugStreamer)<<"Track0"<<
541 // pt->PropagateTo(fkParam->GetInnerRadiusLow());
542 if (pt->GetKinkIndex(0)<=0){ //don't propagate daughter tracks
543 pt->PropagateTo(fkParam->GetInnerRadiusLow());
546 if (( pt->GetPoints()[2]- pt->GetPoints()[0])>5 && pt->GetPoints()[3]>0.8){
548 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
549 iotrack.SetTPCPoints(pt->GetPoints());
550 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
551 iotrack.SetV0Indexes(pt->GetV0Indexes());
552 // iotrack.SetTPCpid(pt->fTPCr);
553 //iotrack.SetTPCindex(i);
554 fEvent->AddTrack(&iotrack);
558 if ( (pt->GetNumberOfClusters()>70)&& (Float_t(pt->GetNumberOfClusters())/Float_t(pt->GetNFoundable()))>0.55) {
560 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
561 iotrack.SetTPCPoints(pt->GetPoints());
562 //iotrack.SetTPCindex(i);
563 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
564 iotrack.SetV0Indexes(pt->GetV0Indexes());
565 // iotrack.SetTPCpid(pt->fTPCr);
566 fEvent->AddTrack(&iotrack);
570 // short tracks - maybe decays
572 if ( (pt->GetNumberOfClusters()>30) && (Float_t(pt->GetNumberOfClusters())/Float_t(pt->GetNFoundable()))>0.70) {
573 Int_t found,foundable,shared;
574 pt->GetClusterStatistic(0,60,found, foundable,shared,kFALSE);
575 if ( (found>20) && (pt->GetNShared()/float(pt->GetNumberOfClusters())<0.2)){
577 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
578 //iotrack.SetTPCindex(i);
579 iotrack.SetTPCPoints(pt->GetPoints());
580 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
581 iotrack.SetV0Indexes(pt->GetV0Indexes());
582 //iotrack.SetTPCpid(pt->fTPCr);
583 fEvent->AddTrack(&iotrack);
588 if ( (pt->GetNumberOfClusters()>20) && (Float_t(pt->GetNumberOfClusters())/Float_t(pt->GetNFoundable()))>0.8) {
589 Int_t found,foundable,shared;
590 pt->GetClusterStatistic(0,60,found, foundable,shared,kFALSE);
591 if (found<20) continue;
592 if (pt->GetNShared()/float(pt->GetNumberOfClusters())>0.2) continue;
595 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
596 iotrack.SetTPCPoints(pt->GetPoints());
597 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
598 iotrack.SetV0Indexes(pt->GetV0Indexes());
599 //iotrack.SetTPCpid(pt->fTPCr);
600 //iotrack.SetTPCindex(i);
601 fEvent->AddTrack(&iotrack);
604 // short tracks - secondaties
606 if ( (pt->GetNumberOfClusters()>30) ) {
607 Int_t found,foundable,shared;
608 pt->GetClusterStatistic(128,158,found, foundable,shared,kFALSE);
609 if ( (found>20) && (pt->GetNShared()/float(pt->GetNumberOfClusters())<0.2) &&float(found)/float(foundable)>0.8){
611 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
612 iotrack.SetTPCPoints(pt->GetPoints());
613 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
614 iotrack.SetV0Indexes(pt->GetV0Indexes());
615 //iotrack.SetTPCpid(pt->fTPCr);
616 //iotrack.SetTPCindex(i);
617 fEvent->AddTrack(&iotrack);
622 if ( (pt->GetNumberOfClusters()>15)) {
623 Int_t found,foundable,shared;
624 pt->GetClusterStatistic(138,158,found, foundable,shared,kFALSE);
625 if (found<15) continue;
626 if (foundable<=0) continue;
627 if (pt->GetNShared()/float(pt->GetNumberOfClusters())>0.2) continue;
628 if (float(found)/float(foundable)<0.8) continue;
631 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
632 iotrack.SetTPCPoints(pt->GetPoints());
633 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
634 iotrack.SetV0Indexes(pt->GetV0Indexes());
635 // iotrack.SetTPCpid(pt->fTPCr);
636 //iotrack.SetTPCindex(i);
637 fEvent->AddTrack(&iotrack);
641 // >> account for suppressed tracks in the kink indices (RS)
642 int nESDtracks = fEvent->GetNumberOfTracks();
643 for (int it=nESDtracks;it--;) {
644 AliESDtrack* esdTr = fEvent->GetTrack(it);
645 if (!esdTr || !esdTr->GetKinkIndex(0)) continue;
646 for (int ik=0;ik<3;ik++) {
648 if (!(knkId=esdTr->GetKinkIndex(ik))) break; // no more kinks for this track
649 AliESDkink* kink = fEvent->GetKink(TMath::Abs(knkId)-1);
651 AliError(Form("ESDTrack%d refers to non-existing kink %d",it,TMath::Abs(knkId)-1));
654 kink->SetIndex(it, knkId<0 ? 0:1); // update track index of the kink: mother at 0, daughter at 1
657 // << account for suppressed tracks in the kink indices (RS)
658 AliInfo(Form("Number of filled ESDs-\t%d\n",fEvent->GetNumberOfTracks()));
666 Double_t AliTPCtrackerMI::ErrY2(AliTPCseed* seed, const AliTPCclusterMI * cl){
669 // Use calibrated cluster error from OCDB
671 AliTPCClusterParam * clparam = AliTPCcalibDB::Instance()->GetClusterParam();
673 Float_t z = TMath::Abs(fkParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
674 Int_t ctype = cl->GetType();
675 Int_t type = (cl->GetRow()<63) ? 0: (cl->GetRow()>126) ? 1:2;
676 Double_t angle = seed->GetSnp()*seed->GetSnp();
677 angle = TMath::Sqrt(TMath::Abs(angle/(1.-angle)));
678 Double_t erry2 = clparam->GetError0Par(0,type, z,angle);
680 erry2+=0.5; // edge cluster
683 seed->SetErrorY2(erry2);
687 //calculate look-up table at the beginning
688 // static Bool_t ginit = kFALSE;
689 // static Float_t gnoise1,gnoise2,gnoise3;
690 // static Float_t ggg1[10000];
691 // static Float_t ggg2[10000];
692 // static Float_t ggg3[10000];
693 // static Float_t glandau1[10000];
694 // static Float_t glandau2[10000];
695 // static Float_t glandau3[10000];
697 // static Float_t gcor01[500];
698 // static Float_t gcor02[500];
699 // static Float_t gcorp[500];
703 // if (ginit==kFALSE){
704 // for (Int_t i=1;i<500;i++){
705 // Float_t rsigma = float(i)/100.;
706 // gcor02[i] = TMath::Max(0.78 +TMath::Exp(7.4*(rsigma-1.2)),0.6);
707 // gcor01[i] = TMath::Max(0.72 +TMath::Exp(3.36*(rsigma-1.2)),0.6);
708 // gcorp[i] = TMath::Max(TMath::Power((rsigma+0.5),1.5),1.2);
712 // for (Int_t i=3;i<10000;i++){
716 // Float_t amp = float(i);
717 // Float_t padlength =0.75;
718 // gnoise1 = 0.0004/padlength;
719 // Float_t nel = 0.268*amp;
720 // Float_t nprim = 0.155*amp;
721 // ggg1[i] = fkParam->GetDiffT()*fkParam->GetDiffT()*(2+0.001*nel/(padlength*padlength))/nel;
722 // glandau1[i] = (2.+0.12*nprim)*0.5* (2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
723 // if (glandau1[i]>1) glandau1[i]=1;
724 // glandau1[i]*=padlength*padlength/12.;
728 // gnoise2 = 0.0004/padlength;
730 // nprim = 0.133*amp;
731 // ggg2[i] = fkParam->GetDiffT()*fkParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
732 // glandau2[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
733 // if (glandau2[i]>1) glandau2[i]=1;
734 // glandau2[i]*=padlength*padlength/12.;
739 // gnoise3 = 0.0004/padlength;
741 // nprim = 0.133*amp;
742 // ggg3[i] = fkParam->GetDiffT()*fkParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
743 // glandau3[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
744 // if (glandau3[i]>1) glandau3[i]=1;
745 // glandau3[i]*=padlength*padlength/12.;
753 // Int_t amp = int(TMath::Abs(cl->GetQ()));
755 // seed->SetErrorY2(1.);
759 // Float_t z = TMath::Abs(fkParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
760 // Int_t ctype = cl->GetType();
761 // Float_t padlength= GetPadPitchLength(seed->GetRow());
762 // Double_t angle2 = seed->GetSnp()*seed->GetSnp();
763 // angle2 = angle2/(1-angle2);
765 // //cluster "quality"
766 // Int_t rsigmay = int(100.*cl->GetSigmaY2()/(seed->GetCurrentSigmaY2()));
769 // if (fSectors==fInnerSec){
770 // snoise2 = gnoise1;
771 // res = ggg1[amp]*z+glandau1[amp]*angle2;
772 // if (ctype==0) res *= gcor01[rsigmay];
775 // res*= gcorp[rsigmay];
779 // if (padlength<1.1){
780 // snoise2 = gnoise2;
781 // res = ggg2[amp]*z+glandau2[amp]*angle2;
782 // if (ctype==0) res *= gcor02[rsigmay];
785 // res*= gcorp[rsigmay];
789 // snoise2 = gnoise3;
790 // res = ggg3[amp]*z+glandau3[amp]*angle2;
791 // if (ctype==0) res *= gcor02[rsigmay];
794 // res*= gcorp[rsigmay];
801 // res*=2.4; // overestimate error 2 times
805 // if (res<2*snoise2)
808 // seed->SetErrorY2(res);
816 Double_t AliTPCtrackerMI::ErrZ2(AliTPCseed* seed, const AliTPCclusterMI * cl){
819 // Use calibrated cluster error from OCDB
821 AliTPCClusterParam * clparam = AliTPCcalibDB::Instance()->GetClusterParam();
823 Float_t z = TMath::Abs(fkParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
824 Int_t ctype = cl->GetType();
825 Int_t type = (cl->GetRow()<63) ? 0: (cl->GetRow()>126) ? 1:2;
827 Double_t angle2 = seed->GetSnp()*seed->GetSnp();
828 angle2 = seed->GetTgl()*seed->GetTgl()*(1+angle2/(1-angle2));
829 Double_t angle = TMath::Sqrt(TMath::Abs(angle2));
830 Double_t errz2 = clparam->GetError0Par(1,type, z,angle);
832 errz2+=0.5; // edge cluster
835 seed->SetErrorZ2(errz2);
841 // //seed->SetErrorY2(0.1);
843 // //calculate look-up table at the beginning
844 // static Bool_t ginit = kFALSE;
845 // static Float_t gnoise1,gnoise2,gnoise3;
846 // static Float_t ggg1[10000];
847 // static Float_t ggg2[10000];
848 // static Float_t ggg3[10000];
849 // static Float_t glandau1[10000];
850 // static Float_t glandau2[10000];
851 // static Float_t glandau3[10000];
853 // static Float_t gcor01[1000];
854 // static Float_t gcor02[1000];
855 // static Float_t gcorp[1000];
859 // if (ginit==kFALSE){
860 // for (Int_t i=1;i<1000;i++){
861 // Float_t rsigma = float(i)/100.;
862 // gcor02[i] = TMath::Max(0.81 +TMath::Exp(6.8*(rsigma-1.2)),0.6);
863 // gcor01[i] = TMath::Max(0.72 +TMath::Exp(2.04*(rsigma-1.2)),0.6);
864 // gcorp[i] = TMath::Max(TMath::Power((rsigma+0.5),1.5),1.2);
868 // for (Int_t i=3;i<10000;i++){
872 // Float_t amp = float(i);
873 // Float_t padlength =0.75;
874 // gnoise1 = 0.0004/padlength;
875 // Float_t nel = 0.268*amp;
876 // Float_t nprim = 0.155*amp;
877 // ggg1[i] = fkParam->GetDiffT()*fkParam->GetDiffT()*(2+0.001*nel/(padlength*padlength))/nel;
878 // glandau1[i] = (2.+0.12*nprim)*0.5* (2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
879 // if (glandau1[i]>1) glandau1[i]=1;
880 // glandau1[i]*=padlength*padlength/12.;
884 // gnoise2 = 0.0004/padlength;
886 // nprim = 0.133*amp;
887 // ggg2[i] = fkParam->GetDiffT()*fkParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
888 // glandau2[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
889 // if (glandau2[i]>1) glandau2[i]=1;
890 // glandau2[i]*=padlength*padlength/12.;
895 // gnoise3 = 0.0004/padlength;
897 // nprim = 0.133*amp;
898 // ggg3[i] = fkParam->GetDiffT()*fkParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
899 // glandau3[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
900 // if (glandau3[i]>1) glandau3[i]=1;
901 // glandau3[i]*=padlength*padlength/12.;
909 // Int_t amp = int(TMath::Abs(cl->GetQ()));
911 // seed->SetErrorY2(1.);
915 // Float_t z = TMath::Abs(fkParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
916 // Int_t ctype = cl->GetType();
917 // Float_t padlength= GetPadPitchLength(seed->GetRow());
919 // Double_t angle2 = seed->GetSnp()*seed->GetSnp();
920 // // if (angle2<0.6) angle2 = 0.6;
921 // angle2 = seed->GetTgl()*seed->GetTgl()*(1+angle2/(1-angle2));
923 // //cluster "quality"
924 // Int_t rsigmaz = int(100.*cl->GetSigmaZ2()/(seed->GetCurrentSigmaZ2()));
927 // if (fSectors==fInnerSec){
928 // snoise2 = gnoise1;
929 // res = ggg1[amp]*z+glandau1[amp]*angle2;
930 // if (ctype==0) res *= gcor01[rsigmaz];
933 // res*= gcorp[rsigmaz];
937 // if (padlength<1.1){
938 // snoise2 = gnoise2;
939 // res = ggg2[amp]*z+glandau2[amp]*angle2;
940 // if (ctype==0) res *= gcor02[rsigmaz];
943 // res*= gcorp[rsigmaz];
947 // snoise2 = gnoise3;
948 // res = ggg3[amp]*z+glandau3[amp]*angle2;
949 // if (ctype==0) res *= gcor02[rsigmaz];
952 // res*= gcorp[rsigmaz];
961 // if ((ctype<0) &&<70){
966 // if (res<2*snoise2)
968 // if (res>3) res =3;
969 // seed->SetErrorZ2(res);
977 void AliTPCtrackerMI::RotateToLocal(AliTPCseed *seed)
979 //rotate to track "local coordinata
980 Float_t x = seed->GetX();
981 Float_t y = seed->GetY();
982 Float_t ymax = x*TMath::Tan(0.5*fSectors->GetAlpha());
985 seed->SetRelativeSector((seed->GetRelativeSector()+1) % fN);
986 if (!seed->Rotate(fSectors->GetAlpha()))
988 } else if (y <-ymax) {
989 seed->SetRelativeSector((seed->GetRelativeSector()-1+fN) % fN);
990 if (!seed->Rotate(-fSectors->GetAlpha()))
998 //_____________________________________________________________________________
999 Double_t AliTPCtrackerMI::F1old(Double_t x1,Double_t y1,
1000 Double_t x2,Double_t y2,
1001 Double_t x3,Double_t y3) const
1003 //-----------------------------------------------------------------
1004 // Initial approximation of the track curvature
1005 //-----------------------------------------------------------------
1006 Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
1007 Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
1008 (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
1009 Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
1010 (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
1012 Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
1013 if ( xr*xr+yr*yr<=0.00000000000001) return 100;
1014 return -xr*yr/sqrt(xr*xr+yr*yr);
1019 //_____________________________________________________________________________
1020 Double_t AliTPCtrackerMI::F1(Double_t x1,Double_t y1,
1021 Double_t x2,Double_t y2,
1022 Double_t x3,Double_t y3) const
1024 //-----------------------------------------------------------------
1025 // Initial approximation of the track curvature
1026 //-----------------------------------------------------------------
1032 Double_t det = x3*y2-x2*y3;
1033 if (TMath::Abs(det)<1e-10){
1037 Double_t u = 0.5* (x2*(x2-x3)+y2*(y2-y3))/det;
1038 Double_t x0 = x3*0.5-y3*u;
1039 Double_t y0 = y3*0.5+x3*u;
1040 Double_t c2 = 1/TMath::Sqrt(x0*x0+y0*y0);
1046 Double_t AliTPCtrackerMI::F2(Double_t x1,Double_t y1,
1047 Double_t x2,Double_t y2,
1048 Double_t x3,Double_t y3) const
1050 //-----------------------------------------------------------------
1051 // Initial approximation of the track curvature
1052 //-----------------------------------------------------------------
1058 Double_t det = x3*y2-x2*y3;
1059 if (TMath::Abs(det)<1e-10) {
1063 Double_t u = 0.5* (x2*(x2-x3)+y2*(y2-y3))/det;
1064 Double_t x0 = x3*0.5-y3*u;
1065 Double_t y0 = y3*0.5+x3*u;
1066 Double_t c2 = 1/TMath::Sqrt(x0*x0+y0*y0);
1075 //_____________________________________________________________________________
1076 Double_t AliTPCtrackerMI::F2old(Double_t x1,Double_t y1,
1077 Double_t x2,Double_t y2,
1078 Double_t x3,Double_t y3) const
1080 //-----------------------------------------------------------------
1081 // Initial approximation of the track curvature times center of curvature
1082 //-----------------------------------------------------------------
1083 Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
1084 Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
1085 (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
1086 Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
1087 (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
1089 Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
1091 return -a/(d*y1-b)*xr/sqrt(xr*xr+yr*yr);
1094 //_____________________________________________________________________________
1095 Double_t AliTPCtrackerMI::F3(Double_t x1,Double_t y1,
1096 Double_t x2,Double_t y2,
1097 Double_t z1,Double_t z2) const
1099 //-----------------------------------------------------------------
1100 // Initial approximation of the tangent of the track dip angle
1101 //-----------------------------------------------------------------
1102 return (z1 - z2)/sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
1106 Double_t AliTPCtrackerMI::F3n(Double_t x1,Double_t y1,
1107 Double_t x2,Double_t y2,
1108 Double_t z1,Double_t z2, Double_t c) const
1110 //-----------------------------------------------------------------
1111 // Initial approximation of the tangent of the track dip angle
1112 //-----------------------------------------------------------------
1116 //angle1 = (z1-z2)*c/(TMath::ASin(c*x1-ni)-TMath::ASin(c*x2-ni));
1118 Double_t d = TMath::Sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
1119 if (TMath::Abs(d*c*0.5)>1) return 0;
1120 // Double_t angle2 = TMath::ASin(d*c*0.5);
1121 // Double_t angle2 = AliTPCFastMath::FastAsin(d*c*0.5);
1122 Double_t angle2 = (d*c*0.5>0.1)? TMath::ASin(d*c*0.5): AliTPCFastMath::FastAsin(d*c*0.5);
1124 angle2 = (z1-z2)*c/(angle2*2.);
1128 Bool_t AliTPCtrackerMI::GetProlongation(Double_t x1, Double_t x2, Double_t x[5], Double_t &y, Double_t &z) const
1129 {//-----------------------------------------------------------------
1130 // This function find proloncation of a track to a reference plane x=x2.
1131 //-----------------------------------------------------------------
1135 if (TMath::Abs(x[4]*x1 - x[2]) >= 0.999) {
1139 Double_t c1=x[4]*x1 - x[2], r1=TMath::Sqrt((1.-c1)*(1.+c1));
1140 Double_t c2=x[4]*x2 - x[2], r2=TMath::Sqrt((1.-c2)*(1.+c2));
1144 Double_t dy = dx*(c1+c2)/(r1+r2);
1147 Double_t delta = x[4]*dx*(c1+c2)/(c1*r2 + c2*r1);
1149 if (TMath::Abs(delta)>0.01){
1150 dz = x[3]*TMath::ASin(delta)/x[4];
1152 dz = x[3]*AliTPCFastMath::FastAsin(delta)/x[4];
1155 //dz = x[3]*AliTPCFastMath::FastAsin(delta)/x[4];
1163 Int_t AliTPCtrackerMI::LoadClusters (TTree *const tree)
1168 return LoadClusters();
1172 Int_t AliTPCtrackerMI::LoadClusters(const TObjArray *arr)
1175 // load clusters to the memory
1176 AliTPCClustersRow *clrow = new AliTPCClustersRow("AliTPCclusterMI");
1177 Int_t lower = arr->LowerBound();
1178 Int_t entries = arr->GetEntriesFast();
1180 for (Int_t i=lower; i<entries; i++) {
1181 clrow = (AliTPCClustersRow*) arr->At(i);
1182 if(!clrow) continue;
1183 if(!clrow->GetArray()) continue;
1187 fkParam->AdjustSectorRow(clrow->GetID(),sec,row);
1189 for (Int_t icl=0; icl<clrow->GetArray()->GetEntriesFast(); icl++){
1190 Transform((AliTPCclusterMI*)(clrow->GetArray()->At(icl)));
1193 if (clrow->GetArray()->GetEntriesFast()<=0) continue;
1194 AliTPCtrackerRow * tpcrow=0;
1197 tpcrow = &(fInnerSec[sec%fkNIS][row]);
1201 tpcrow = &(fOuterSec[(sec-fkNIS*2)%fkNOS][row]);
1202 left = (sec-fkNIS*2)/fkNOS;
1205 tpcrow->SetN1(clrow->GetArray()->GetEntriesFast());
1206 tpcrow->SetClusters1(new AliTPCclusterMI[tpcrow->GetN1()]);
1207 for (Int_t j=0;j<tpcrow->GetN1();++j)
1208 tpcrow->SetCluster1(j, *(AliTPCclusterMI*)(clrow->GetArray()->At(j)));
1211 tpcrow->SetN2(clrow->GetArray()->GetEntriesFast());
1212 tpcrow->SetClusters2(new AliTPCclusterMI[tpcrow->GetN2()]);
1213 for (Int_t j=0;j<tpcrow->GetN2();++j)
1214 tpcrow->SetCluster2(j,*(AliTPCclusterMI*)(clrow->GetArray()->At(j)));
1216 clrow->GetArray()->Clear("C");
1225 Int_t AliTPCtrackerMI::LoadClusters(const TClonesArray *arr)
1228 // load clusters to the memory from one
1231 AliTPCclusterMI *clust=0;
1232 Int_t count[72][96] = { {0} , {0} };
1234 // loop over clusters
1235 for (Int_t icl=0; icl<arr->GetEntriesFast(); icl++) {
1236 clust = (AliTPCclusterMI*)arr->At(icl);
1237 if(!clust) continue;
1238 //printf("cluster: det %d, row %d \n", clust->GetDetector(),clust->GetRow());
1240 // transform clusters
1243 // count clusters per pad row
1244 count[clust->GetDetector()][clust->GetRow()]++;
1247 // insert clusters to sectors
1248 for (Int_t icl=0; icl<arr->GetEntriesFast(); icl++) {
1249 clust = (AliTPCclusterMI*)arr->At(icl);
1250 if(!clust) continue;
1252 Int_t sec = clust->GetDetector();
1253 Int_t row = clust->GetRow();
1255 // filter overlapping pad rows needed by HLT
1256 if(sec<fkNIS*2) { //IROCs
1257 if(row == 30) continue;
1260 if(row == 27 || row == 76) continue;
1266 fInnerSec[sec%fkNIS].InsertCluster(clust, count[sec][row], fkParam);
1269 left = (sec-fkNIS*2)/fkNOS;
1270 fOuterSec[(sec-fkNIS*2)%fkNOS].InsertCluster(clust, count[sec][row], fkParam);
1274 // Load functions must be called behind LoadCluster(TClonesArray*)
1276 //LoadOuterSectors();
1277 //LoadInnerSectors();
1283 Int_t AliTPCtrackerMI::LoadClusters()
1286 // load clusters to the memory
1287 AliTPCClustersRow *clrow= new AliTPCClustersRow("AliTPCclusterMI");
1289 // TTree * tree = fClustersArray.GetTree();
1291 TTree * tree = fInput;
1292 TBranch * br = tree->GetBranch("Segment");
1293 br->SetAddress(&clrow);
1295 // Conversion of pad, row coordinates in local tracking coords.
1296 // Could be skipped here; is already done in clusterfinder
1298 Int_t j=Int_t(tree->GetEntries());
1299 for (Int_t i=0; i<j; i++) {
1303 fkParam->AdjustSectorRow(clrow->GetID(),sec,row);
1304 for (Int_t icl=0; icl<clrow->GetArray()->GetEntriesFast(); icl++){
1305 Transform((AliTPCclusterMI*)(clrow->GetArray()->At(icl)));
1308 AliTPCtrackerRow * tpcrow=0;
1311 tpcrow = &(fInnerSec[sec%fkNIS][row]);
1315 tpcrow = &(fOuterSec[(sec-fkNIS*2)%fkNOS][row]);
1316 left = (sec-fkNIS*2)/fkNOS;
1319 tpcrow->SetN1(clrow->GetArray()->GetEntriesFast());
1320 tpcrow->SetClusters1(new AliTPCclusterMI[tpcrow->GetN1()]);
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 tpcrow->SetClusters2(new AliTPCclusterMI[tpcrow->GetN2()]);
1327 for (Int_t k=0;k<tpcrow->GetN2();++k)
1328 tpcrow->SetCluster2(k,*(AliTPCclusterMI*)(clrow->GetArray()->At(k)));
1339 void AliTPCtrackerMI::UnloadClusters()
1342 // unload clusters from the memory
1344 Int_t nrows = fOuterSec->GetNRows();
1345 for (Int_t sec = 0;sec<fkNOS;sec++)
1346 for (Int_t row = 0;row<nrows;row++){
1347 AliTPCtrackerRow* tpcrow = &(fOuterSec[sec%fkNOS][row]);
1349 // if (tpcrow->fClusters1) delete []tpcrow->fClusters1;
1350 // if (tpcrow->fClusters2) delete []tpcrow->fClusters2;
1352 tpcrow->ResetClusters();
1355 nrows = fInnerSec->GetNRows();
1356 for (Int_t sec = 0;sec<fkNIS;sec++)
1357 for (Int_t row = 0;row<nrows;row++){
1358 AliTPCtrackerRow* tpcrow = &(fInnerSec[sec%fkNIS][row]);
1360 // if (tpcrow->fClusters1) delete []tpcrow->fClusters1;
1361 //if (tpcrow->fClusters2) delete []tpcrow->fClusters2;
1363 tpcrow->ResetClusters();
1369 void AliTPCtrackerMI::FillClusterArray(TObjArray* array) const{
1371 // Filling cluster to the array - For visualization purposes
1374 nrows = fOuterSec->GetNRows();
1375 for (Int_t sec = 0;sec<fkNOS;sec++)
1376 for (Int_t row = 0;row<nrows;row++){
1377 AliTPCtrackerRow* tpcrow = &(fOuterSec[sec%fkNOS][row]);
1378 if (!tpcrow) continue;
1379 for (Int_t icl = 0;icl<tpcrow->GetN();icl++){
1380 array->AddLast((TObject*)((*tpcrow)[icl]));
1383 nrows = fInnerSec->GetNRows();
1384 for (Int_t sec = 0;sec<fkNIS;sec++)
1385 for (Int_t row = 0;row<nrows;row++){
1386 AliTPCtrackerRow* tpcrow = &(fInnerSec[sec%fkNIS][row]);
1387 if (!tpcrow) continue;
1388 for (Int_t icl = 0;icl<tpcrow->GetN();icl++){
1389 array->AddLast((TObject*)(*tpcrow)[icl]);
1395 void AliTPCtrackerMI::Transform(AliTPCclusterMI * cluster){
1399 AliTPCcalibDB * calibDB = AliTPCcalibDB::Instance();
1400 AliTPCTransform *transform = calibDB->GetTransform() ;
1402 AliFatal("Tranformations not in calibDB");
1405 transform->SetCurrentRecoParam((AliTPCRecoParam*)AliTPCReconstructor::GetRecoParam());
1406 Double_t x[3]={cluster->GetRow(),cluster->GetPad(),cluster->GetTimeBin()};
1407 Int_t i[1]={cluster->GetDetector()};
1408 transform->Transform(x,i,0,1);
1409 // if (cluster->GetDetector()%36>17){
1414 // in debug mode check the transformation
1416 if (AliTPCReconstructor::StreamLevel()>2) {
1418 cluster->GetGlobalXYZ(gx);
1419 Int_t event = (fEvent==NULL)? 0: fEvent->GetEventNumberInFile();
1420 TTreeSRedirector &cstream = *fDebugStreamer;
1421 cstream<<"Transform"<<
1432 cluster->SetX(x[0]);
1433 cluster->SetY(x[1]);
1434 cluster->SetZ(x[2]);
1439 //if (!fkParam->IsGeoRead()) fkParam->ReadGeoMatrices();
1440 if (AliTPCReconstructor::GetRecoParam()->GetUseSectorAlignment() && (!calibDB->HasAlignmentOCDB())){
1441 TGeoHMatrix *mat = fkParam->GetClusterMatrix(cluster->GetDetector());
1443 Double_t pos[3]= {cluster->GetX(),cluster->GetY(),cluster->GetZ()};
1444 Double_t posC[3]={cluster->GetX(),cluster->GetY(),cluster->GetZ()};
1445 if (mat) mat->LocalToMaster(pos,posC);
1447 // chack Loading of Geo matrices from GeoManager - TEMPORARY FIX
1449 cluster->SetX(posC[0]);
1450 cluster->SetY(posC[1]);
1451 cluster->SetZ(posC[2]);
1455 //_____________________________________________________________________________
1456 Int_t AliTPCtrackerMI::LoadOuterSectors() {
1457 //-----------------------------------------------------------------
1458 // This function fills outer TPC sectors with clusters.
1459 //-----------------------------------------------------------------
1460 Int_t nrows = fOuterSec->GetNRows();
1462 for (Int_t sec = 0;sec<fkNOS;sec++)
1463 for (Int_t row = 0;row<nrows;row++){
1464 AliTPCtrackerRow* tpcrow = &(fOuterSec[sec%fkNOS][row]);
1465 Int_t sec2 = sec+2*fkNIS;
1467 Int_t ncl = tpcrow->GetN1();
1469 AliTPCclusterMI *c= (tpcrow->GetCluster1(ncl));
1470 index=(((sec2<<8)+row)<<16)+ncl;
1471 tpcrow->InsertCluster(c,index);
1474 ncl = tpcrow->GetN2();
1476 AliTPCclusterMI *c= (tpcrow->GetCluster2(ncl));
1477 index=((((sec2+fkNOS)<<8)+row)<<16)+ncl;
1478 tpcrow->InsertCluster(c,index);
1481 // write indexes for fast acces
1483 for (Int_t i=0;i<510;i++)
1484 tpcrow->SetFastCluster(i,-1);
1485 for (Int_t i=0;i<tpcrow->GetN();i++){
1486 Int_t zi = Int_t((*tpcrow)[i]->GetZ()+255.);
1487 tpcrow->SetFastCluster(zi,i); // write index
1490 for (Int_t i=0;i<510;i++){
1491 if (tpcrow->GetFastCluster(i)<0)
1492 tpcrow->SetFastCluster(i,last);
1494 last = tpcrow->GetFastCluster(i);
1503 //_____________________________________________________________________________
1504 Int_t AliTPCtrackerMI::LoadInnerSectors() {
1505 //-----------------------------------------------------------------
1506 // This function fills inner TPC sectors with clusters.
1507 //-----------------------------------------------------------------
1508 Int_t nrows = fInnerSec->GetNRows();
1510 for (Int_t sec = 0;sec<fkNIS;sec++)
1511 for (Int_t row = 0;row<nrows;row++){
1512 AliTPCtrackerRow* tpcrow = &(fInnerSec[sec%fkNIS][row]);
1515 Int_t ncl = tpcrow->GetN1();
1517 AliTPCclusterMI *c= (tpcrow->GetCluster1(ncl));
1518 index=(((sec<<8)+row)<<16)+ncl;
1519 tpcrow->InsertCluster(c,index);
1522 ncl = tpcrow->GetN2();
1524 AliTPCclusterMI *c= (tpcrow->GetCluster2(ncl));
1525 index=((((sec+fkNIS)<<8)+row)<<16)+ncl;
1526 tpcrow->InsertCluster(c,index);
1529 // write indexes for fast acces
1531 for (Int_t i=0;i<510;i++)
1532 tpcrow->SetFastCluster(i,-1);
1533 for (Int_t i=0;i<tpcrow->GetN();i++){
1534 Int_t zi = Int_t((*tpcrow)[i]->GetZ()+255.);
1535 tpcrow->SetFastCluster(zi,i); // write index
1538 for (Int_t i=0;i<510;i++){
1539 if (tpcrow->GetFastCluster(i)<0)
1540 tpcrow->SetFastCluster(i,last);
1542 last = tpcrow->GetFastCluster(i);
1554 //_________________________________________________________________________
1555 AliTPCclusterMI *AliTPCtrackerMI::GetClusterMI(Int_t index) const {
1556 //--------------------------------------------------------------------
1557 // Return pointer to a given cluster
1558 //--------------------------------------------------------------------
1559 if (index<0) return 0; // no cluster
1560 Int_t sec=(index&0xff000000)>>24;
1561 Int_t row=(index&0x00ff0000)>>16;
1562 Int_t ncl=(index&0x00007fff)>>00;
1564 const AliTPCtrackerRow * tpcrow=0;
1565 AliTPCclusterMI * clrow =0;
1567 if (sec<0 || sec>=fkNIS*4) {
1568 AliWarning(Form("Wrong sector %d",sec));
1573 AliTPCtrackerSector& tracksec = fInnerSec[sec%fkNIS];
1574 if (tracksec.GetNRows()<=row) return 0;
1575 tpcrow = &(tracksec[row]);
1576 if (tpcrow==0) return 0;
1579 if (tpcrow->GetN1()<=ncl) return 0;
1580 clrow = tpcrow->GetClusters1();
1583 if (tpcrow->GetN2()<=ncl) return 0;
1584 clrow = tpcrow->GetClusters2();
1588 AliTPCtrackerSector& tracksec = fOuterSec[(sec-fkNIS*2)%fkNOS];
1589 if (tracksec.GetNRows()<=row) return 0;
1590 tpcrow = &(tracksec[row]);
1591 if (tpcrow==0) return 0;
1593 if (sec-2*fkNIS<fkNOS) {
1594 if (tpcrow->GetN1()<=ncl) return 0;
1595 clrow = tpcrow->GetClusters1();
1598 if (tpcrow->GetN2()<=ncl) return 0;
1599 clrow = tpcrow->GetClusters2();
1603 return &(clrow[ncl]);
1609 Int_t AliTPCtrackerMI::FollowToNext(AliTPCseed& t, Int_t nr) {
1610 //-----------------------------------------------------------------
1611 // This function tries to find a track prolongation to next pad row
1612 //-----------------------------------------------------------------
1614 Double_t x= GetXrow(nr), ymax=GetMaxY(nr);
1617 AliTPCclusterMI *cl=0;
1618 Int_t tpcindex= t.GetClusterIndex2(nr);
1620 // update current shape info every 5 pad-row
1621 // if ( (nr%5==0) || t.GetNumberOfClusters()<2 || (t.fCurrentSigmaY2<0.0001) ){
1625 if (fIteration>0 && tpcindex>=-1){ //if we have already clusters
1627 if (tpcindex==-1) return 0; //track in dead zone
1628 if (tpcindex >= 0){ //
1629 cl = t.GetClusterPointer(nr);
1630 //if (cl==0) cl = GetClusterMI(tpcindex);
1631 if (!cl) cl = GetClusterMI(tpcindex);
1632 t.SetCurrentClusterIndex1(tpcindex);
1635 Int_t relativesector = ((tpcindex&0xff000000)>>24)%18; // if previously accepted cluster in different sector
1636 Float_t angle = relativesector*fSectors->GetAlpha()+fSectors->GetAlphaShift();
1638 if (angle<-TMath::Pi()) angle += 2*TMath::Pi();
1639 if (angle>=TMath::Pi()) angle -= 2*TMath::Pi();
1641 if (TMath::Abs(angle-t.GetAlpha())>0.001){
1642 Double_t rotation = angle-t.GetAlpha();
1643 t.SetRelativeSector(relativesector);
1644 if (!t.Rotate(rotation)) return 0;
1646 if (!t.PropagateTo(x)) return 0;
1648 t.SetCurrentCluster(cl);
1650 Int_t accept = AcceptCluster(&t,t.GetCurrentCluster());
1651 if ((tpcindex&0x8000)==0) accept =0;
1653 //if founded cluster is acceptible
1654 if (cl->IsUsed(11)) { // id cluster is shared inrease uncertainty
1655 t.SetErrorY2(t.GetErrorY2()+0.03);
1656 t.SetErrorZ2(t.GetErrorZ2()+0.03);
1657 t.SetErrorY2(t.GetErrorY2()*3);
1658 t.SetErrorZ2(t.GetErrorZ2()*3);
1660 t.SetNFoundable(t.GetNFoundable()+1);
1661 UpdateTrack(&t,accept);
1664 else { // Remove old cluster from track
1665 t.SetClusterIndex(nr, -3);
1666 t.SetClusterPointer(nr, 0);
1670 if (TMath::Abs(t.GetSnp())>AliTPCReconstructor::GetMaxSnpTracker()) return 0; // cut on angle
1671 if (fIteration>1 && IsFindable(t)){
1672 // not look for new cluster during refitting
1673 t.SetNFoundable(t.GetNFoundable()+1);
1678 // if (TMath::Abs(t.GetSnp())>0.95 || TMath::Abs(x*t.GetC()-t.GetEta())>0.95) return 0;// patch 28 fev 06
1679 if (!t.PropagateTo(x)) {
1680 if (fIteration==0) t.SetRemoval(10);
1683 Double_t y = t.GetY();
1684 if (TMath::Abs(y)>ymax){
1686 t.SetRelativeSector((t.GetRelativeSector()+1) % fN);
1687 if (!t.Rotate(fSectors->GetAlpha()))
1689 } else if (y <-ymax) {
1690 t.SetRelativeSector((t.GetRelativeSector()-1+fN) % fN);
1691 if (!t.Rotate(-fSectors->GetAlpha()))
1694 if (!t.PropagateTo(x)) {
1695 if (fIteration==0) t.SetRemoval(10);
1701 Double_t z=t.GetZ();
1704 if (!IsActive(t.GetRelativeSector(),nr)) {
1706 t.SetClusterIndex2(nr,-1);
1709 //AliInfo(Form("A - Sector%d phi %f - alpha %f", t.fRelativeSector,y/x, t.GetAlpha()));
1710 Bool_t isActive = IsActive(t.GetRelativeSector(),nr);
1711 Bool_t isActive2 = (nr>=fInnerSec->GetNRows()) ? fOuterSec[t.GetRelativeSector()][nr-fInnerSec->GetNRows()].GetN()>0:fInnerSec[t.GetRelativeSector()][nr].GetN()>0;
1713 if (!isActive || !isActive2) return 0;
1715 const AliTPCtrackerRow &krow=GetRow(t.GetRelativeSector(),nr);
1716 if ( (t.GetSigmaY2()<0) || t.GetSigmaZ2()<0) return 0;
1718 Double_t roadz = 1.;
1720 if (TMath::Abs(TMath::Abs(y)-ymax)<krow.GetDeadZone()){
1722 t.SetClusterIndex2(nr,-1);
1728 // if (TMath::Abs(z)<(AliTPCReconstructor::GetCtgRange()*x+10) && TMath::Abs(z)<fkParam->GetZLength(0) && (TMath::Abs(t.GetSnp())<AliTPCReconstructor::GetMaxSnpTracker()))
1729 t.SetNFoundable(t.GetNFoundable()+1);
1735 // cl = krow.FindNearest2(y+10.,z,roady,roadz,index);
1736 cl = krow.FindNearest2(y,z,roady,roadz,index);
1737 if (cl) t.SetCurrentClusterIndex1(krow.GetIndex(index));
1740 t.SetCurrentCluster(cl);
1742 if (fIteration==2&&cl->IsUsed(10)) return 0;
1743 Int_t accept = AcceptCluster(&t,t.GetCurrentCluster());
1744 if (fIteration==2&&cl->IsUsed(11)) {
1745 t.SetErrorY2(t.GetErrorY2()+0.03);
1746 t.SetErrorZ2(t.GetErrorZ2()+0.03);
1747 t.SetErrorY2(t.GetErrorY2()*3);
1748 t.SetErrorZ2(t.GetErrorZ2()*3);
1751 if (t.fCurrentCluster->IsUsed(10)){
1756 if (t.fNShared>0.7*t.GetNumberOfClusters()) {
1762 if (accept<3) UpdateTrack(&t,accept);
1765 if ( fIteration==0 && t.GetNFoundable()*0.5 > t.GetNumberOfClusters()) t.SetRemoval(10);
1773 //_________________________________________________________________________
1774 Bool_t AliTPCtrackerMI::GetTrackPoint(Int_t index, AliTrackPoint &p ) const
1776 // Get track space point by index
1777 // return false in case the cluster doesn't exist
1778 AliTPCclusterMI *cl = GetClusterMI(index);
1779 if (!cl) return kFALSE;
1780 Int_t sector = (index&0xff000000)>>24;
1781 // Int_t row = (index&0x00ff0000)>>16;
1783 // xyz[0] = fkParam->GetPadRowRadii(sector,row);
1784 xyz[0] = cl->GetX();
1785 xyz[1] = cl->GetY();
1786 xyz[2] = cl->GetZ();
1788 fkParam->AdjustCosSin(sector,cos,sin);
1789 Float_t x = cos*xyz[0]-sin*xyz[1];
1790 Float_t y = cos*xyz[1]+sin*xyz[0];
1792 Float_t sigmaY2 = 0.027*cl->GetSigmaY2();
1793 if (sector < fkParam->GetNInnerSector()) sigmaY2 *= 2.07;
1794 Float_t sigmaZ2 = 0.066*cl->GetSigmaZ2();
1795 if (sector < fkParam->GetNInnerSector()) sigmaZ2 *= 1.77;
1796 cov[0] = sin*sin*sigmaY2;
1797 cov[1] = -sin*cos*sigmaY2;
1799 cov[3] = cos*cos*sigmaY2;
1802 p.SetXYZ(x,y,xyz[2],cov);
1803 AliGeomManager::ELayerID iLayer;
1805 if (sector < fkParam->GetNInnerSector()) {
1806 iLayer = AliGeomManager::kTPC1;
1810 iLayer = AliGeomManager::kTPC2;
1811 idet = sector - fkParam->GetNInnerSector();
1813 UShort_t volid = AliGeomManager::LayerToVolUID(iLayer,idet);
1814 p.SetVolumeID(volid);
1820 Int_t AliTPCtrackerMI::UpdateClusters(AliTPCseed& t, Int_t nr) {
1821 //-----------------------------------------------------------------
1822 // This function tries to find a track prolongation to next pad row
1823 //-----------------------------------------------------------------
1824 t.SetCurrentCluster(0);
1825 t.SetCurrentClusterIndex1(-3);
1827 Double_t xt=t.GetX();
1828 Int_t row = GetRowNumber(xt)-1;
1829 Double_t ymax= GetMaxY(nr);
1831 if (row < nr) return 1; // don't prolongate if not information until now -
1832 // if (TMath::Abs(t.GetSnp())>0.9 && t.GetNumberOfClusters()>40. && fIteration!=2) {
1834 // return 0; // not prolongate strongly inclined tracks
1836 // if (TMath::Abs(t.GetSnp())>0.95) {
1838 // return 0; // not prolongate strongly inclined tracks
1839 // }// patch 28 fev 06
1841 Double_t x= GetXrow(nr);
1843 //t.PropagateTo(x+0.02);
1844 //t.PropagateTo(x+0.01);
1845 if (!t.PropagateTo(x)){
1852 if (TMath::Abs(y)>ymax){
1854 t.SetRelativeSector((t.GetRelativeSector()+1) % fN);
1855 if (!t.Rotate(fSectors->GetAlpha()))
1857 } else if (y <-ymax) {
1858 t.SetRelativeSector((t.GetRelativeSector()-1+fN) % fN);
1859 if (!t.Rotate(-fSectors->GetAlpha()))
1862 // if (!t.PropagateTo(x)){
1869 if (TMath::Abs(t.GetSnp())>AliTPCReconstructor::GetMaxSnpTracker()) return 0;
1871 if (!IsActive(t.GetRelativeSector(),nr)) {
1873 t.SetClusterIndex2(nr,-1);
1876 //AliInfo(Form("A - Sector%d phi %f - alpha %f", t.fRelativeSector,y/x, t.GetAlpha()));
1878 AliTPCtrackerRow &krow=GetRow(t.GetRelativeSector(),nr);
1880 if (TMath::Abs(TMath::Abs(y)-ymax)<krow.GetDeadZone()){
1882 t.SetClusterIndex2(nr,-1);
1888 // if (TMath::Abs(t.GetZ())<(AliTPCReconstructor::GetCtgRange()*t.GetX()+10) && (TMath::Abs(t.GetSnp())<AliTPCReconstructor::GetMaxSnpTracker()))
1889 if (IsFindable(t)) t.SetNFoundable(t.GetNFoundable()+1);
1895 if ( (nr%2==0) || t.GetNumberOfClusters()<2){
1896 // t.fCurrentSigmaY = GetSigmaY(&t);
1897 //t.fCurrentSigmaZ = GetSigmaZ(&t);
1901 AliTPCclusterMI *cl=0;
1904 Double_t roady = 1.;
1905 Double_t roadz = 1.;
1909 index = t.GetClusterIndex2(nr);
1910 if ( (index >= 0) && (index&0x8000)==0){
1911 cl = t.GetClusterPointer(nr);
1912 if ( (cl==0) && (index >= 0)) cl = GetClusterMI(index);
1913 t.SetCurrentClusterIndex1(index);
1915 t.SetCurrentCluster(cl);
1921 // if (index<0) return 0;
1922 UInt_t uindex = TMath::Abs(index);
1925 //cl = krow.FindNearest2(y+10,z,roady,roadz,uindex);
1926 cl = krow.FindNearest2(y,z,roady,roadz,uindex);
1929 if (cl) t.SetCurrentClusterIndex1(krow.GetIndex(uindex));
1930 t.SetCurrentCluster(cl);
1936 Int_t AliTPCtrackerMI::FollowToNextCluster(AliTPCseed & t, Int_t nr) {
1937 //-----------------------------------------------------------------
1938 // This function tries to find a track prolongation to next pad row
1939 //-----------------------------------------------------------------
1941 //update error according neighborhoud
1943 if (t.GetCurrentCluster()) {
1945 Int_t accept = AcceptCluster(&t,t.GetCurrentCluster());
1947 if (t.GetCurrentCluster()->IsUsed(10)){
1952 t.SetNShared(t.GetNShared()+1);
1953 if (t.GetNShared()>0.7*t.GetNumberOfClusters()) {
1958 if (fIteration>0) accept = 0;
1959 if (accept<3) UpdateTrack(&t,accept);
1963 if ( t.GetNumberOfClusters()>18 && ( (t.GetSigmaY2()+t.GetSigmaZ2())>0.16)) t.SetRemoval(10);
1964 if ( t.GetNumberOfClusters()>18 && t.GetChi2()/t.GetNumberOfClusters()>6 ) t.SetRemoval(10);
1966 if (( (t.GetNFoundable()*0.5 > t.GetNumberOfClusters()) || t.GetNoCluster()>15)) t.SetRemoval(10);
1974 //_____________________________________________________________________________
1975 Int_t AliTPCtrackerMI::FollowProlongation(AliTPCseed& t, Int_t rf, Int_t step, Bool_t fromSeeds) {
1976 //-----------------------------------------------------------------
1977 // This function tries to find a track prolongation.
1978 //-----------------------------------------------------------------
1979 Double_t xt=t.GetX();
1981 Double_t alpha=t.GetAlpha();
1982 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
1983 if (alpha < 0. ) alpha += 2.*TMath::Pi();
1985 t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN);
1987 Int_t first = GetRowNumber(xt);
1992 for (Int_t nr= first; nr>=rf; nr-=step) {
1994 if (t.GetKinkIndexes()[0]>0){
1995 for (Int_t i=0;i<3;i++){
1996 Int_t index = t.GetKinkIndexes()[i];
1997 if (index==0) break;
1998 if (index<0) continue;
2000 AliKink * kink = (AliKink*)fEvent->GetKink(index-1);
2002 printf("PROBLEM\n");
2005 Int_t kinkrow = kink->GetTPCRow0()+2+Int_t(0.5/(0.05+kink->GetAngle(2)));
2007 AliExternalTrackParam paramd(t);
2008 kink->SetDaughter(paramd);
2009 kink->SetStatus(2,5);
2016 if (nr==80) t.UpdateReference();
2017 if (nr<fInnerSec->GetNRows())
2018 fSectors = fInnerSec;
2020 fSectors = fOuterSec;
2021 if (FollowToNext(t,nr)==0)
2034 Int_t AliTPCtrackerMI::FollowBackProlongation(AliTPCseed& t, Int_t rf, Bool_t fromSeeds) {
2035 //-----------------------------------------------------------------
2036 // This function tries to find a track prolongation.
2037 //-----------------------------------------------------------------
2039 Double_t xt=t.GetX();
2040 Double_t alpha=t.GetAlpha();
2041 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
2042 if (alpha < 0. ) alpha += 2.*TMath::Pi();
2043 t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN);
2045 Int_t first = t.GetFirstPoint();
2046 Int_t ri = GetRowNumber(xt);
2050 if (first<ri) first = ri;
2052 if (first<0) first=0;
2053 for (Int_t nr=first; nr<=rf; nr++) {
2054 // if ( (TMath::Abs(t.GetSnp())>0.95)) break;//patch 28 fev 06
2055 if (t.GetKinkIndexes()[0]<0){
2056 for (Int_t i=0;i<3;i++){
2057 Int_t index = t.GetKinkIndexes()[i];
2058 if (index==0) break;
2059 if (index>0) continue;
2060 index = TMath::Abs(index);
2061 AliKink * kink = (AliKink*)fEvent->GetKink(index-1);
2063 printf("PROBLEM\n");
2066 Int_t kinkrow = kink->GetTPCRow0()-2-Int_t(0.5/(0.05+kink->GetAngle(2)));
2068 AliExternalTrackParam paramm(t);
2069 kink->SetMother(paramm);
2070 kink->SetStatus(2,1);
2077 if (nr<fInnerSec->GetNRows())
2078 fSectors = fInnerSec;
2080 fSectors = fOuterSec;
2091 Float_t AliTPCtrackerMI::OverlapFactor(AliTPCseed * s1, AliTPCseed * s2, Int_t &sum1, Int_t & sum2)
2093 // overlapping factor
2099 Float_t dz2 =(s1->GetZ() - s2->GetZ());
2102 Float_t dy2 =TMath::Abs((s1->GetY() - s2->GetY()));
2104 Float_t distance = TMath::Sqrt(dz2+dy2);
2105 if (distance>4.) return 0; // if there are far away - not overlap - to reduce combinatorics
2108 Int_t firstpoint = TMath::Min(s1->GetFirstPoint(),s2->GetFirstPoint());
2109 Int_t lastpoint = TMath::Max(s1->GetLastPoint(),s2->GetLastPoint());
2114 if (firstpoint>lastpoint) {
2115 firstpoint =lastpoint;
2120 for (Int_t i=firstpoint-1;i<lastpoint+1;i++){
2121 if (s1->GetClusterIndex2(i)>0) sum1++;
2122 if (s2->GetClusterIndex2(i)>0) sum2++;
2123 if (s1->GetClusterIndex2(i)==s2->GetClusterIndex2(i) && s1->GetClusterIndex2(i)>0) {
2127 if (sum<5) return 0;
2129 Float_t summin = TMath::Min(sum1+1,sum2+1);
2130 Float_t ratio = (sum+1)/Float_t(summin);
2134 void AliTPCtrackerMI::SignShared(AliTPCseed * s1, AliTPCseed * s2)
2138 Float_t thetaCut = 0.2;//+10.*TMath::Sqrt(s1->GetSigmaTglZ()+ s2->GetSigmaTglZ());
2139 if (TMath::Abs(s1->GetTgl()-s2->GetTgl())>thetaCut) return;
2140 Float_t minCl = TMath::Min(s1->GetNumberOfClusters(),s2->GetNumberOfClusters());
2141 Int_t cutN0 = TMath::Max(5,TMath::Nint(0.1*minCl));
2146 //Int_t firstpoint = TMath::Max(s1->GetFirstPoint(),s2->GetFirstPoint());
2147 //Int_t lastpoint = TMath::Min(s1->GetLastPoint(),s2->GetLastPoint());
2148 Int_t firstpoint = 0;
2149 Int_t lastpoint = 160;
2151 // if (firstpoint>=lastpoint-5) return;;
2153 for (Int_t i=firstpoint;i<lastpoint;i++){
2154 // if ( (s1->GetClusterIndex2(i)&0xFFFF8FFF)==(s2->GetClusterIndex2(i)&0xFFFF8FFF) && s1->GetClusterIndex2(i)>0) {
2155 if ( (s1->GetClusterIndex2(i))==(s2->GetClusterIndex2(i)) && s1->GetClusterIndex2(i)>0) {
2159 if (sumshared>cutN0){
2162 for (Int_t i=firstpoint;i<lastpoint;i++){
2163 // if ( (s1->GetClusterIndex2(i)&0xFFFF8FFF)==(s2->GetClusterIndex2(i)&0xFFFF8FFF) && s1->GetClusterIndex2(i)>0) {
2164 if ( (s1->GetClusterIndex2(i))==(s2->GetClusterIndex2(i)) && s1->GetClusterIndex2(i)>0) {
2165 AliTPCTrackerPoint *p1 = s1->GetTrackPoint(i);
2166 AliTPCTrackerPoint *p2 = s2->GetTrackPoint(i);;
2167 if (s1->IsActive()&&s2->IsActive()){
2168 p1->SetShared(kTRUE);
2169 p2->SetShared(kTRUE);
2175 if (sumshared>cutN0){
2176 for (Int_t i=0;i<4;i++){
2177 if (s1->GetOverlapLabel(3*i)==0){
2178 s1->SetOverlapLabel(3*i, s2->GetLabel());
2179 s1->SetOverlapLabel(3*i+1,sumshared);
2180 s1->SetOverlapLabel(3*i+2,s2->GetUniqueID());
2184 for (Int_t i=0;i<4;i++){
2185 if (s2->GetOverlapLabel(3*i)==0){
2186 s2->SetOverlapLabel(3*i, s1->GetLabel());
2187 s2->SetOverlapLabel(3*i+1,sumshared);
2188 s2->SetOverlapLabel(3*i+2,s1->GetUniqueID());
2195 void AliTPCtrackerMI::SignShared(TObjArray * arr)
2198 //sort trackss according sectors
2200 for (Int_t i=0; i<arr->GetEntriesFast(); i++) {
2201 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2203 //if (pt) RotateToLocal(pt);
2207 arr->Sort(); // sorting according relative sectors
2208 arr->Expand(arr->GetEntries());
2211 Int_t nseed=arr->GetEntriesFast();
2212 for (Int_t i=0; i<nseed; i++) {
2213 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2215 for (Int_t j=0;j<12;j++){
2216 pt->SetOverlapLabel(j,0);
2219 for (Int_t i=0; i<nseed; i++) {
2220 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2222 if (pt->GetRemoval()>10) continue;
2223 for (Int_t j=i+1; j<nseed; j++){
2224 AliTPCseed *pt2=(AliTPCseed*)arr->UncheckedAt(j);
2225 if (TMath::Abs(pt->GetRelativeSector()-pt2->GetRelativeSector())>1) continue;
2227 if (pt2->GetRemoval()<=10) {
2228 //if ( TMath::Abs(pt->GetRelativeSector()-pt2->GetRelativeSector())>0) break;
2236 void AliTPCtrackerMI::SortTracks(TObjArray * arr, Int_t mode) const
2239 //sort tracks in array according mode criteria
2240 Int_t nseed = arr->GetEntriesFast();
2241 for (Int_t i=0; i<nseed; i++) {
2242 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2253 void AliTPCtrackerMI::RemoveUsed2(TObjArray * arr, Float_t factor1, Float_t factor2, Int_t minimal)
2256 // Loop over all tracks and remove overlaped tracks (with lower quality)
2258 // 1. Unsign clusters
2259 // 2. Sort tracks according quality
2260 // Quality is defined by the number of cluster between first and last points
2262 // 3. Loop over tracks - decreasing quality order
2263 // a.) remove - If the fraction of shared cluster less than factor (1- n or 2)
2264 // b.) remove - If the minimal number of clusters less than minimal and not ITS
2265 // c.) if track accepted - sign clusters
2267 //Called in - AliTPCtrackerMI::Clusters2Tracks()
2268 // - AliTPCtrackerMI::PropagateBack()
2269 // - AliTPCtrackerMI::RefitInward()
2272 // factor1 - factor for constrained
2273 // factor2 - for non constrained tracks
2274 // if (Float_t(shared+1)/Float_t(found+1)>factor) - DELETE
2278 Int_t nseed = arr->GetEntriesFast();
2279 Float_t * quality = new Float_t[nseed];
2280 Int_t * indexes = new Int_t[nseed];
2284 for (Int_t i=0; i<nseed; i++) {
2285 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2290 pt->UpdatePoints(); //select first last max dens points
2291 Float_t * points = pt->GetPoints();
2292 if (points[3]<0.8) quality[i] =-1;
2293 quality[i] = (points[2]-points[0])+pt->GetNumberOfClusters();
2294 //prefer high momenta tracks if overlaps
2295 quality[i] *= TMath::Sqrt(TMath::Abs(pt->Pt())+0.5);
2297 TMath::Sort(nseed,quality,indexes);
2300 for (Int_t itrack=0; itrack<nseed; itrack++) {
2301 Int_t trackindex = indexes[itrack];
2302 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(trackindex);
2305 if (quality[trackindex]<0){
2306 MarkSeedFree( arr->RemoveAt(trackindex) );
2311 Int_t first = Int_t(pt->GetPoints()[0]);
2312 Int_t last = Int_t(pt->GetPoints()[2]);
2313 Double_t factor = (pt->GetBConstrain()) ? factor1: factor2;
2315 Int_t found,foundable,shared;
2316 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
2317 // pt->GetClusterStatistic(0,160, found, foundable,shared,kFALSE);
2318 Bool_t itsgold =kFALSE;
2321 if (pt->GetESD()->GetITSclusters(dummy)>4) itsgold= kTRUE;
2325 if (Float_t(shared+1)/Float_t(found+1)>factor){
2326 if (pt->GetKinkIndexes()[0]!=0) continue; //don't remove tracks - part of the kinks
2327 if( AliTPCReconstructor::StreamLevel()>3){
2328 TTreeSRedirector &cstream = *fDebugStreamer;
2329 cstream<<"RemoveUsed"<<
2330 "iter="<<fIteration<<
2334 MarkSeedFree( arr->RemoveAt(trackindex) );
2337 if (pt->GetNumberOfClusters()<50&&(found-0.5*shared)<minimal){ //remove short tracks
2338 if (pt->GetKinkIndexes()[0]!=0) continue; //don't remove tracks - part of the kinks
2339 if( AliTPCReconstructor::StreamLevel()>3){
2340 TTreeSRedirector &cstream = *fDebugStreamer;
2341 cstream<<"RemoveShort"<<
2342 "iter="<<fIteration<<
2346 MarkSeedFree( arr->RemoveAt(trackindex) );
2352 //if (sharedfactor>0.4) continue;
2353 if (pt->GetKinkIndexes()[0]>0) continue;
2354 //Remove tracks with undefined properties - seems
2355 if (pt->GetSigmaY2()<kAlmost0) continue; // ? what is the origin ?
2357 for (Int_t i=first; i<last; i++) {
2358 Int_t index=pt->GetClusterIndex2(i);
2359 // if (index<0 || index&0x8000 ) continue;
2360 if (index<0 || index&0x8000 ) continue;
2361 AliTPCclusterMI *c= pt->GetClusterPointer(i);
2368 Info("RemoveUsed","\n*****\nNumber of good tracks after shared removal\t%d\n",fNtracks);
2374 void AliTPCtrackerMI::DumpClusters(Int_t iter, TObjArray *trackArray)
2377 // Dump clusters after reco
2378 // signed and unsigned cluster can be visualized
2379 // 1. Unsign all cluster
2380 // 2. Sign all used clusters
2383 Int_t nseed = trackArray->GetEntries();
2384 for (Int_t i=0; i<nseed; i++){
2385 AliTPCseed *pt=(AliTPCseed*)trackArray->UncheckedAt(i);
2389 Bool_t isKink=pt->GetKinkIndex(0)!=0;
2390 for (Int_t j=0; j<160; ++j) {
2391 Int_t index=pt->GetClusterIndex2(j);
2392 if (index<0) continue;
2393 AliTPCclusterMI *c= pt->GetClusterPointer(j);
2395 if (isKink) c->Use(100); // kink
2396 c->Use(10); // by default usage 10
2401 for (Int_t sec=0;sec<fkNIS;sec++){
2402 for (Int_t row=0;row<fInnerSec->GetNRows();row++){
2403 AliTPCclusterMI *cl = fInnerSec[sec][row].GetClusters1();
2404 for (Int_t icl =0;icl< fInnerSec[sec][row].GetN1();icl++){
2405 Float_t gx[3]; cl[icl].GetGlobalXYZ(gx);
2406 (*fDebugStreamer)<<"clDump"<<
2414 cl = fInnerSec[sec][row].GetClusters2();
2415 for (Int_t icl =0;icl< fInnerSec[sec][row].GetN2();icl++){
2416 Float_t gx[3]; cl[icl].GetGlobalXYZ(gx);
2417 (*fDebugStreamer)<<"clDump"<<
2428 for (Int_t sec=0;sec<fkNOS;sec++){
2429 for (Int_t row=0;row<fOuterSec->GetNRows();row++){
2430 AliTPCclusterMI *cl = fOuterSec[sec][row].GetClusters1();
2431 for (Int_t icl =0;icl< fOuterSec[sec][row].GetN1();icl++){
2432 Float_t gx[3]; cl[icl].GetGlobalXYZ(gx);
2433 (*fDebugStreamer)<<"clDump"<<
2441 cl = fOuterSec[sec][row].GetClusters2();
2442 for (Int_t icl =0;icl< fOuterSec[sec][row].GetN2();icl++){
2443 Float_t gx[3]; cl[icl].GetGlobalXYZ(gx);
2444 (*fDebugStreamer)<<"clDump"<<
2456 void AliTPCtrackerMI::UnsignClusters()
2459 // loop over all clusters and unsign them
2462 for (Int_t sec=0;sec<fkNIS;sec++){
2463 for (Int_t row=0;row<fInnerSec->GetNRows();row++){
2464 AliTPCclusterMI *cl = fInnerSec[sec][row].GetClusters1();
2465 for (Int_t icl =0;icl< fInnerSec[sec][row].GetN1();icl++)
2466 // if (cl[icl].IsUsed(10))
2468 cl = fInnerSec[sec][row].GetClusters2();
2469 for (Int_t icl =0;icl< fInnerSec[sec][row].GetN2();icl++)
2470 //if (cl[icl].IsUsed(10))
2475 for (Int_t sec=0;sec<fkNOS;sec++){
2476 for (Int_t row=0;row<fOuterSec->GetNRows();row++){
2477 AliTPCclusterMI *cl = fOuterSec[sec][row].GetClusters1();
2478 for (Int_t icl =0;icl< fOuterSec[sec][row].GetN1();icl++)
2479 //if (cl[icl].IsUsed(10))
2481 cl = fOuterSec[sec][row].GetClusters2();
2482 for (Int_t icl =0;icl< fOuterSec[sec][row].GetN2();icl++)
2483 //if (cl[icl].IsUsed(10))
2492 void AliTPCtrackerMI::SignClusters(const TObjArray * arr, Float_t fnumber, Float_t fdensity)
2495 //sign clusters to be "used"
2497 // snumber and sdensity sign number of sigmas - bellow mean value to be accepted
2498 // loop over "primaries"
2512 Int_t nseed = arr->GetEntriesFast();
2513 for (Int_t i=0; i<nseed; i++) {
2514 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2518 if (!(pt->IsActive())) continue;
2519 Float_t dens = pt->GetNumberOfClusters()/Float_t(pt->GetNFoundable());
2520 if ( (dens>0.7) && (pt->GetNumberOfClusters()>70)){
2522 sumdens2+= dens*dens;
2523 sumn += pt->GetNumberOfClusters();
2524 sumn2 += pt->GetNumberOfClusters()*pt->GetNumberOfClusters();
2525 Float_t chi2 = pt->GetChi2()/pt->GetNumberOfClusters();
2528 sumchi2 +=chi2*chi2;
2533 Float_t mdensity = 0.9;
2534 Float_t meann = 130;
2535 Float_t meanchi = 1;
2536 Float_t sdensity = 0.1;
2537 Float_t smeann = 10;
2538 Float_t smeanchi =0.4;
2542 mdensity = sumdens/sum;
2544 meanchi = sumchi/sum;
2546 sdensity = sumdens2/sum-mdensity*mdensity;
2548 sdensity = TMath::Sqrt(sdensity);
2552 smeann = sumn2/sum-meann*meann;
2554 smeann = TMath::Sqrt(smeann);
2558 smeanchi = sumchi2/sum - meanchi*meanchi;
2560 smeanchi = TMath::Sqrt(smeanchi);
2566 //REMOVE SHORT DELTAS or tracks going out of sensitive volume of TPC
2568 for (Int_t i=0; i<nseed; i++) {
2569 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2573 if (pt->GetBSigned()) continue;
2574 if (pt->GetBConstrain()) continue;
2575 //if (!(pt->IsActive())) continue;
2577 Int_t found,foundable,shared;
2578 pt->GetClusterStatistic(0,160,found, foundable,shared);
2579 if (shared/float(found)>0.3) {
2580 if (shared/float(found)>0.9 ){
2581 //MarkSeedFree( arr->RemoveAt(i) );
2586 Bool_t isok =kFALSE;
2587 if ( (pt->GetNShared()/pt->GetNumberOfClusters()<0.5) &&pt->GetNumberOfClusters()>60)
2589 if ((TMath::Abs(1/pt->GetC())<100.) && (pt->GetNShared()/pt->GetNumberOfClusters()<0.7))
2591 if (TMath::Abs(pt->GetZ()/pt->GetX())>1.1)
2593 if ( (TMath::Abs(pt->GetSnp()>0.7) && pt->GetD(0,0)>60.))
2597 for (Int_t j=0; j<160; ++j) {
2598 Int_t index=pt->GetClusterIndex2(j);
2599 if (index<0) continue;
2600 AliTPCclusterMI *c= pt->GetClusterPointer(j);
2602 //if (!(c->IsUsed(10))) c->Use();
2609 Double_t maxchi = meanchi+2.*smeanchi;
2611 for (Int_t i=0; i<nseed; i++) {
2612 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2616 //if (!(pt->IsActive())) continue;
2617 if (pt->GetBSigned()) continue;
2618 Double_t chi = pt->GetChi2()/pt->GetNumberOfClusters();
2619 if (chi>maxchi) continue;
2622 Float_t dens = pt->GetNumberOfClusters()/Float_t(pt->GetNFoundable());
2624 //sign only tracks with enoug big density at the beginning
2626 if ((pt->GetDensityFirst(40)<0.75) && pt->GetNumberOfClusters()<meann) continue;
2629 Double_t mindens = TMath::Max(double(mdensity-sdensity*fdensity*bfactor),0.65);
2630 Double_t minn = TMath::Max(Int_t(meann-fnumber*smeann*bfactor),50);
2632 // if (pt->fBConstrain) mindens = TMath::Max(mdensity-sdensity*fdensity*bfactor,0.65);
2633 if ( (pt->GetRemoval()==10) && (pt->GetSnp()>0.8)&&(dens>mindens))
2636 if ((dens>mindens && pt->GetNumberOfClusters()>minn) && chi<maxchi ){
2637 //Int_t noc=pt->GetNumberOfClusters();
2638 pt->SetBSigned(kTRUE);
2639 for (Int_t j=0; j<160; ++j) {
2641 Int_t index=pt->GetClusterIndex2(j);
2642 if (index<0) continue;
2643 AliTPCclusterMI *c= pt->GetClusterPointer(j);
2645 // if (!(c->IsUsed(10))) c->Use();
2650 // gLastCheck = nseed;
2659 Int_t AliTPCtrackerMI::RefitInward(AliESDEvent *event)
2662 // back propagation of ESD tracks
2665 if (!event) return 0;
2666 const Int_t kMaxFriendTracks=2000;
2668 // extract correction object for multiplicity dependence of dEdx
2669 TObjArray * gainCalibArray = AliTPCcalibDB::Instance()->GetTimeGainSplinesRun(event->GetRunNumber());
2671 AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ;
2673 AliFatal("Tranformations not in RefitInward");
2676 transform->SetCurrentRecoParam((AliTPCRecoParam*)AliTPCReconstructor::GetRecoParam());
2677 const AliTPCRecoParam * recoParam = AliTPCcalibDB::Instance()->GetTransform()->GetCurrentRecoParam();
2678 Int_t nContribut = event->GetNumberOfTracks();
2679 TGraphErrors * graphMultDependenceDeDx = 0x0;
2680 if (recoParam && recoParam->GetUseMultiplicityCorrectionDedx() && gainCalibArray) {
2681 if (recoParam->GetUseTotCharge()) {
2682 graphMultDependenceDeDx = (TGraphErrors *) gainCalibArray->FindObject("TGRAPHERRORS_MEANQTOT_MULTIPLICITYDEPENDENCE_BEAM_ALL");
2684 graphMultDependenceDeDx = (TGraphErrors *) gainCalibArray->FindObject("TGRAPHERRORS_MEANQMAX_MULTIPLICITYDEPENDENCE_BEAM_ALL");
2690 //PrepareForProlongation(fSeeds,1);
2691 PropagateForward2(fSeeds);
2692 RemoveUsed2(fSeeds,0.4,0.4,20);
2694 TObjArray arraySeed(fSeeds->GetEntries());
2695 for (Int_t i=0;i<fSeeds->GetEntries();i++) {
2696 arraySeed.AddAt(fSeeds->At(i),i);
2698 SignShared(&arraySeed);
2699 // FindCurling(fSeeds, event,2); // find multi found tracks
2700 FindSplitted(fSeeds, event,2); // find multi found tracks
2701 if (AliTPCReconstructor::StreamLevel()>5) FindMultiMC(fSeeds, fEvent,2); // find multi found tracks
2704 Int_t nseed = fSeeds->GetEntriesFast();
2705 for (Int_t i=0;i<nseed;i++){
2706 AliTPCseed * seed = (AliTPCseed*) fSeeds->UncheckedAt(i);
2707 if (!seed) continue;
2708 if (seed->GetKinkIndex(0)>0) UpdateKinkQualityD(seed); // update quality informations for kinks
2709 AliESDtrack *esd=event->GetTrack(i);
2711 if (seed->GetNumberOfClusters()<60 && seed->GetNumberOfClusters()<(esd->GetTPCclusters(0) -5)*0.8){
2712 AliExternalTrackParam paramIn;
2713 AliExternalTrackParam paramOut;
2714 Int_t ncl = seed->RefitTrack(seed,¶mIn,¶mOut);
2715 if (AliTPCReconstructor::StreamLevel()>2) {
2716 (*fDebugStreamer)<<"RecoverIn"<<
2720 "pout.="<<¶mOut<<
2725 seed->Set(paramIn.GetX(),paramIn.GetAlpha(),paramIn.GetParameter(),paramIn.GetCovariance());
2726 seed->SetNumberOfClusters(ncl);
2730 seed->PropagateTo(fkParam->GetInnerRadiusLow());
2731 seed->UpdatePoints();
2732 AddCovariance(seed);
2733 MakeESDBitmaps(seed, esd);
2734 seed->CookdEdx(0.02,0.6);
2735 CookLabel(seed,0.1); //For comparison only
2737 if (AliTPCReconstructor::StreamLevel()>1 && seed!=0) {
2738 TTreeSRedirector &cstream = *fDebugStreamer;
2745 if (seed->GetNumberOfClusters()>15){
2746 esd->UpdateTrackParams(seed,AliESDtrack::kTPCrefit);
2747 esd->SetTPCPoints(seed->GetPoints());
2748 esd->SetTPCPointsF(seed->GetNFoundable());
2749 Int_t ndedx = seed->GetNCDEDX(0);
2750 Float_t sdedx = seed->GetSDEDX(0);
2751 Float_t dedx = seed->GetdEdx();
2752 // apply mutliplicity dependent dEdx correction if available
2753 if (graphMultDependenceDeDx) {
2754 Double_t corrGain = AliTPCcalibDButil::EvalGraphConst(graphMultDependenceDeDx, nContribut);
2755 dedx += (1 - corrGain)*50.; // MIP is normalized to 50
2757 esd->SetTPCsignal(dedx, sdedx, ndedx);
2759 // fill new dEdx information
2761 Double32_t signal[4];
2765 for(Int_t iarr=0;iarr<3;iarr++) {
2766 signal[iarr] = seed->GetDEDXregion(iarr+1);
2767 ncl[iarr] = seed->GetNCDEDX(iarr+1);
2768 nrows[iarr] = seed->GetNCDEDXInclThres(iarr+1);
2770 signal[3] = seed->GetDEDXregion(4);
2772 AliTPCdEdxInfo * infoTpcPid = new AliTPCdEdxInfo();
2773 infoTpcPid->SetTPCSignalRegionInfo(signal, ncl, nrows);
2774 esd->SetTPCdEdxInfo(infoTpcPid);
2776 // add seed to the esd track in Calib level
2778 Bool_t storeFriend = gRandom->Rndm()<(kMaxFriendTracks)/Float_t(nseed);
2779 if (AliTPCReconstructor::StreamLevel()>0 &&storeFriend){
2780 // RS: this is the only place where the seed is created not in the pool,
2781 // since it should belong to ESDevent
2782 AliTPCseed * seedCopy = new AliTPCseed(*seed, kTRUE);
2783 esd->AddCalibObject(seedCopy);
2788 //printf("problem\n");
2791 //FindKinks(fSeeds,event);
2792 if (AliTPCReconstructor::StreamLevel()>3) DumpClusters(2,fSeeds);
2793 Info("RefitInward","Number of refitted tracks %d",ntracks);
2798 Int_t AliTPCtrackerMI::PropagateBack(AliESDEvent *event)
2801 // back propagation of ESD tracks
2803 if (!event) return 0;
2807 PropagateBack(fSeeds);
2808 RemoveUsed2(fSeeds,0.4,0.4,20);
2809 //FindCurling(fSeeds, fEvent,1);
2810 FindSplitted(fSeeds, event,1); // find multi found tracks
2811 if (AliTPCReconstructor::StreamLevel()>5) FindMultiMC(fSeeds, fEvent,1); // find multi found tracks
2814 Int_t nseed = fSeeds->GetEntriesFast();
2816 for (Int_t i=0;i<nseed;i++){
2817 AliTPCseed * seed = (AliTPCseed*) fSeeds->UncheckedAt(i);
2818 if (!seed) continue;
2819 if (seed->GetKinkIndex(0)<0) UpdateKinkQualityM(seed); // update quality informations for kinks
2820 seed->UpdatePoints();
2821 AddCovariance(seed);
2822 AliESDtrack *esd=event->GetTrack(i);
2823 if (!esd) continue; //never happen
2824 if (seed->GetNumberOfClusters()<60 && seed->GetNumberOfClusters()<(esd->GetTPCclusters(0) -5)*0.8){
2825 AliExternalTrackParam paramIn;
2826 AliExternalTrackParam paramOut;
2827 Int_t ncl = seed->RefitTrack(seed,¶mIn,¶mOut);
2828 if (AliTPCReconstructor::StreamLevel()>2) {
2829 (*fDebugStreamer)<<"RecoverBack"<<
2833 "pout.="<<¶mOut<<
2838 seed->Set(paramOut.GetX(),paramOut.GetAlpha(),paramOut.GetParameter(),paramOut.GetCovariance());
2839 seed->SetNumberOfClusters(ncl);
2842 seed->CookdEdx(0.02,0.6);
2843 CookLabel(seed,0.1); //For comparison only
2844 if (seed->GetNumberOfClusters()>15){
2845 esd->UpdateTrackParams(seed,AliESDtrack::kTPCout);
2846 esd->SetTPCPoints(seed->GetPoints());
2847 esd->SetTPCPointsF(seed->GetNFoundable());
2848 Int_t ndedx = seed->GetNCDEDX(0);
2849 Float_t sdedx = seed->GetSDEDX(0);
2850 Float_t dedx = seed->GetdEdx();
2851 esd->SetTPCsignal(dedx, sdedx, ndedx);
2853 Int_t eventnumber = event->GetEventNumberInFile();// patch 28 fev 06
2854 // This is most likely NOT the event number you'd like to use. It has nothing to do with the 'real' event number
2855 if (AliTPCReconstructor::StreamLevel()>1 && esd) {
2856 (*fDebugStreamer)<<"Cback"<<
2859 "EventNrInFile="<<eventnumber<<
2864 if (AliTPCReconstructor::StreamLevel()>3) DumpClusters(1,fSeeds);
2865 //FindKinks(fSeeds,event);
2866 Info("PropagateBack","Number of back propagated tracks %d",ntracks);
2873 void AliTPCtrackerMI::DeleteSeeds()
2882 void AliTPCtrackerMI::ReadSeeds(const AliESDEvent *const event, Int_t direction)
2885 //read seeds from the event
2887 Int_t nentr=event->GetNumberOfTracks();
2889 Info("ReadSeeds", "Number of ESD tracks: %d\n", nentr);
2894 fSeeds = new TObjArray(nentr);
2898 for (Int_t i=0; i<nentr; i++) {
2899 AliESDtrack *esd=event->GetTrack(i);
2900 ULong_t status=esd->GetStatus();
2901 if (!(status&AliESDtrack::kTPCin)) continue;
2902 AliTPCtrack t(*esd);
2903 t.SetNumberOfClusters(0);
2904 // AliTPCseed *seed = new AliTPCseed(t,t.GetAlpha());
2905 AliTPCseed *seed = new( NextFreeSeed() ) AliTPCseed(t/*,t.GetAlpha()*/);
2906 seed->SetPoolID(fLastSeedID);
2907 seed->SetUniqueID(esd->GetID());
2908 AddCovariance(seed); //add systematic ucertainty
2909 for (Int_t ikink=0;ikink<3;ikink++) {
2910 Int_t index = esd->GetKinkIndex(ikink);
2911 seed->GetKinkIndexes()[ikink] = index;
2912 if (index==0) continue;
2913 index = TMath::Abs(index);
2914 AliESDkink * kink = fEvent->GetKink(index-1);
2915 if (kink&&esd->GetKinkIndex(ikink)<0){
2916 if ((status & AliESDtrack::kTRDrefit) != 0) kink->SetStatus(1,2);
2917 if ((status & AliESDtrack::kITSout) != 0) kink->SetStatus(1,0);
2919 if (kink&&esd->GetKinkIndex(ikink)>0){
2920 if ((status & AliESDtrack::kTRDrefit) != 0) kink->SetStatus(1,6);
2921 if ((status & AliESDtrack::kITSout) != 0) kink->SetStatus(1,4);
2925 if (((status&AliESDtrack::kITSout)==0)&&(direction==1)) seed->ResetCovariance(10.);
2926 if ( direction ==2 &&(status & AliESDtrack::kTRDrefit) == 0 ) seed->ResetCovariance(10.);
2927 //if ( direction ==2 && ((status & AliESDtrack::kTPCout) == 0) ) {
2928 // fSeeds->AddAt(0,i);
2929 // MarkSeedFree( seed );
2932 if ( direction ==2 &&(status & AliESDtrack::kTRDrefit) > 0 ) {
2933 Double_t par0[5],par1[5],alpha,x;
2934 esd->GetInnerExternalParameters(alpha,x,par0);
2935 esd->GetExternalParameters(x,par1);
2936 Double_t delta1 = TMath::Abs(par0[4]-par1[4])/(0.000000001+TMath::Abs(par0[4]+par1[4]));
2937 Double_t delta2 = TMath::Abs(par0[3]-par1[3]);
2939 if (esd->GetTRDncls()>0) trdchi2 = esd->GetTRDchi2()/esd->GetTRDncls();
2940 //reset covariance if suspicious
2941 if ( (delta1>0.1) || (delta2>0.006) ||trdchi2>7.)
2942 seed->ResetCovariance(10.);
2947 // rotate to the local coordinate system
2949 fSectors=fInnerSec; fN=fkNIS;
2950 Double_t alpha=seed->GetAlpha();
2951 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
2952 if (alpha < 0. ) alpha += 2.*TMath::Pi();
2953 Int_t ns=Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN;
2954 alpha =ns*fSectors->GetAlpha() + fSectors->GetAlphaShift();
2955 alpha-=seed->GetAlpha();
2956 if (alpha<-TMath::Pi()) alpha += 2*TMath::Pi();
2957 if (alpha>=TMath::Pi()) alpha -= 2*TMath::Pi();
2958 if (TMath::Abs(alpha) > 0.001) { // This should not happen normally
2959 AliWarning(Form("Rotating track over %f",alpha));
2960 if (!seed->Rotate(alpha)) {
2961 MarkSeedFree( seed );
2967 if (esd->GetKinkIndex(0)<=0){
2968 for (Int_t irow=0;irow<160;irow++){
2969 Int_t index = seed->GetClusterIndex2(irow);
2972 AliTPCclusterMI * cl = GetClusterMI(index);
2973 seed->SetClusterPointer(irow,cl);
2975 if ((index & 0x8000)==0){
2976 cl->Use(10); // accepted cluster
2978 cl->Use(6); // close cluster not accepted
2981 Info("ReadSeeds","Not found cluster");
2986 fSeeds->AddAt(seed,i);
2992 //_____________________________________________________________________________
2993 void AliTPCtrackerMI::MakeSeeds3(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, Float_t cuts[4],
2994 Float_t deltay, Int_t ddsec) {
2995 //-----------------------------------------------------------------
2996 // This function creates track seeds.
2997 // SEEDING WITH VERTEX CONSTRAIN
2998 //-----------------------------------------------------------------
2999 // cuts[0] - fP4 cut
3000 // cuts[1] - tan(phi) cut
3001 // cuts[2] - zvertex cut
3002 // cuts[3] - fP3 cut
3010 Double_t x[5], c[15];
3011 // Int_t di = i1-i2;
3013 AliTPCseed * seed = new( NextFreeSeed() ) AliTPCseed();
3014 seed->SetPoolID(fLastSeedID);
3015 Double_t alpha=fSectors->GetAlpha(), shift=fSectors->GetAlphaShift();
3016 Double_t cs=cos(alpha), sn=sin(alpha);
3018 // Double_t x1 =fOuterSec->GetX(i1);
3019 //Double_t xx2=fOuterSec->GetX(i2);
3021 Double_t x1 =GetXrow(i1);
3022 Double_t xx2=GetXrow(i2);
3024 Double_t x3=GetX(), y3=GetY(), z3=GetZ();
3026 Int_t imiddle = (i2+i1)/2; //middle pad row index
3027 Double_t xm = GetXrow(imiddle); // radius of middle pad-row
3028 const AliTPCtrackerRow& krm=GetRow(sec,imiddle); //middle pad -row
3032 const AliTPCtrackerRow& kr1=GetRow(ns,i1);
3033 Double_t ymax = GetMaxY(i1)-kr1.GetDeadZone()-1.5;
3034 Double_t ymaxm = GetMaxY(imiddle)-kr1.GetDeadZone()-1.5;
3037 // change cut on curvature if it can't reach this layer
3038 // maximal curvature set to reach it
3039 Double_t dvertexmax = TMath::Sqrt((x1-x3)*(x1-x3)+(ymax+5-y3)*(ymax+5-y3));
3040 if (dvertexmax*0.5*cuts[0]>0.85){
3041 cuts[0] = 0.85/(dvertexmax*0.5+1.);
3043 Double_t r2min = 1/(cuts[0]*cuts[0]); //minimal square of radius given by cut
3046 if (deltay>0) ddsec = 0;
3047 // loop over clusters
3048 for (Int_t is=0; is < kr1; is++) {
3050 if (kr1[is]->IsUsed(10)) continue;
3051 Double_t y1=kr1[is]->GetY(), z1=kr1[is]->GetZ();
3052 //if (TMath::Abs(y1)>ymax) continue;
3054 if (deltay>0 && TMath::Abs(ymax-TMath::Abs(y1))> deltay ) continue; // seed only at the edge
3056 // find possible directions
3057 Float_t anglez = (z1-z3)/(x1-x3);
3058 Float_t extraz = z1 - anglez*(x1-xx2); // extrapolated z
3061 //find rotation angles relative to line given by vertex and point 1
3062 Double_t dvertex2 = (x1-x3)*(x1-x3)+(y1-y3)*(y1-y3);
3063 Double_t dvertex = TMath::Sqrt(dvertex2);
3064 Double_t angle13 = TMath::ATan((y1-y3)/(x1-x3));
3065 Double_t cs13 = cos(-angle13), sn13 = sin(-angle13);
3068 // loop over 2 sectors
3074 Double_t dddz1=0; // direction of delta inclination in z axis
3081 for (Int_t dsec = dsec1; dsec<=dsec2;dsec++){
3082 Int_t sec2 = sec + dsec;
3084 // AliTPCtrackerRow& kr2 = fOuterSec[(sec2+fkNOS)%fkNOS][i2];
3085 //AliTPCtrackerRow& kr2m = fOuterSec[(sec2+fkNOS)%fkNOS][imiddle];
3086 AliTPCtrackerRow& kr2 = GetRow((sec2+fkNOS)%fkNOS,i2);
3087 AliTPCtrackerRow& kr2m = GetRow((sec2+fkNOS)%fkNOS,imiddle);
3088 Int_t index1 = TMath::Max(kr2.Find(extraz-0.6-dddz1*TMath::Abs(z1)*0.05)-1,0);
3089 Int_t index2 = TMath::Min(kr2.Find(extraz+0.6+dddz2*TMath::Abs(z1)*0.05)+1,kr2);
3091 // rotation angles to p1-p3
3092 Double_t cs13r = cos(-angle13+dsec*alpha)/dvertex, sn13r = sin(-angle13+dsec*alpha)/dvertex;
3093 Double_t x2, y2, z2;
3095 // Double_t dymax = maxangle*TMath::Abs(x1-xx2);
3098 Double_t dxx0 = (xx2-x3)*cs13r;
3099 Double_t dyy0 = (xx2-x3)*sn13r;
3100 for (Int_t js=index1; js < index2; js++) {
3101 const AliTPCclusterMI *kcl = kr2[js];
3102 if (kcl->IsUsed(10)) continue;
3104 //calcutate parameters
3106 Double_t yy0 = dyy0 +(kcl->GetY()-y3)*cs13r;
3108 if (TMath::Abs(yy0)<0.000001) continue;
3109 Double_t xx0 = dxx0 -(kcl->GetY()-y3)*sn13r;
3110 Double_t y0 = 0.5*(xx0*xx0+yy0*yy0-xx0)/yy0;
3111 Double_t r02 = (0.25+y0*y0)*dvertex2;
3112 //curvature (radius) cut
3113 if (r02<r2min) continue;
3117 Double_t c0 = 1/TMath::Sqrt(r02);
3121 //Double_t dfi0 = 2.*TMath::ASin(dvertex*c0*0.5);
3122 //Double_t dfi1 = 2.*TMath::ASin(TMath::Sqrt(yy0*yy0+(1-xx0)*(1-xx0))*dvertex*c0*0.5);
3123 Double_t dfi0 = 2.*AliTPCFastMath::FastAsin(dvertex*c0*0.5);
3124 Double_t dfi1 = 2.*AliTPCFastMath::FastAsin(TMath::Sqrt(yy0*yy0+(1-xx0)*(1-xx0))*dvertex*c0*0.5);
3127 Double_t z0 = kcl->GetZ();
3128 Double_t zzzz2 = z1-(z1-z3)*dfi1/dfi0;
3129 if (TMath::Abs(zzzz2-z0)>0.5) continue;
3132 Double_t dip = (z1-z0)*c0/dfi1;
3133 Double_t x0 = (0.5*cs13+y0*sn13)*dvertex*c0;
3144 x2= xx2*cs-y2*sn*dsec;
3145 y2=+xx2*sn*dsec+y2*cs;
3155 // do we have cluster at the middle ?
3157 GetProlongation(x1,xm,x,ym,zm);
3159 AliTPCclusterMI * cm=0;
3160 if (TMath::Abs(ym)-ymaxm<0){
3161 cm = krm.FindNearest2(ym,zm,1.0,0.6,dummy);
3162 if ((!cm) || (cm->IsUsed(10))) {
3167 // rotate y1 to system 0
3168 // get state vector in rotated system
3169 Double_t yr1 = (-0.5*sn13+y0*cs13)*dvertex*c0;
3170 Double_t xr2 = x0*cs+yr1*sn*dsec;
3171 Double_t xr[5]={kcl->GetY(),kcl->GetZ(), xr2, dip, c0};
3173 GetProlongation(xx2,xm,xr,ym,zm);
3174 if (TMath::Abs(ym)-ymaxm<0){
3175 cm = kr2m.FindNearest2(ym,zm,1.0,0.6,dummy);
3176 if ((!cm) || (cm->IsUsed(10))) {
3186 dym = ym - cm->GetY();
3187 dzm = zm - cm->GetZ();
3194 Double_t sy1=kr1[is]->GetSigmaY2()*2., sz1=kr1[is]->GetSigmaZ2()*2.;
3195 Double_t sy2=kcl->GetSigmaY2()*2., sz2=kcl->GetSigmaZ2()*2.;
3196 //Double_t sy3=400*3./12., sy=0.1, sz=0.1;
3197 Double_t sy3=25000*x[4]*x[4]+0.1, sy=0.1, sz=0.1;
3198 //Double_t sy3=25000*x[4]*x[4]*60+0.5, sy=0.1, sz=0.1;
3200 Double_t f40=(F1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
3201 Double_t f42=(F1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
3202 Double_t f43=(F1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
3203 Double_t f20=(F2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
3204 Double_t f22=(F2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
3205 Double_t f23=(F2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
3207 Double_t f30=(F3(x1,y1+sy,x2,y2,z1,z2)-x[3])/sy;
3208 Double_t f31=(F3(x1,y1,x2,y2,z1+sz,z2)-x[3])/sz;
3209 Double_t f32=(F3(x1,y1,x2,y2+sy,z1,z2)-x[3])/sy;
3210 Double_t f34=(F3(x1,y1,x2,y2,z1,z2+sz)-x[3])/sz;
3214 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3215 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3216 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3217 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3218 c[13]=f30*sy1*f40+f32*sy2*f42;
3219 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3221 // if (!BuildSeed(kr1[is],kcl,0,x1,x2,x3,x,c)) continue;
3223 UInt_t index=kr1.GetIndex(is);
3224 if (seed) {MarkSeedFree(seed); seed = 0;}
3225 AliTPCseed *track = seed = new( NextFreeSeed() ) AliTPCseed(x1, ns*alpha+shift, x, c, index);
3226 seed->SetPoolID(fLastSeedID);
3227 track->SetIsSeeding(kTRUE);
3228 track->SetSeed1(i1);
3229 track->SetSeed2(i2);
3230 track->SetSeedType(3);
3234 FollowProlongation(*track, (i1+i2)/2,1);
3235 Int_t foundable,found,shared;
3236 track->GetClusterStatistic((i1+i2)/2,i1, found, foundable, shared, kTRUE);
3237 if ((found<0.55*foundable) || shared>0.5*found || (track->GetSigmaY2()+track->GetSigmaZ2())>0.5){
3238 MarkSeedFree(seed); seed = 0;
3244 FollowProlongation(*track, i2,1);
3248 track->SetBConstrain(1);
3249 // track->fLastPoint = i1+fInnerSec->GetNRows(); // first cluster in track position
3250 track->SetLastPoint(i1); // first cluster in track position
3251 track->SetFirstPoint(track->GetLastPoint());
3253 if (track->GetNumberOfClusters()<(i1-i2)*0.5 ||
3254 track->GetNumberOfClusters() < track->GetNFoundable()*0.6 ||
3255 track->GetNShared()>0.4*track->GetNumberOfClusters() ) {
3256 MarkSeedFree(seed); seed = 0;
3260 // Z VERTEX CONDITION
3261 Double_t zv, bz=GetBz();
3262 if ( !track->GetZAt(0.,bz,zv) ) continue;
3263 if (TMath::Abs(zv-z3)>cuts[2]) {
3264 FollowProlongation(*track, TMath::Max(i2-20,0));
3265 if ( !track->GetZAt(0.,bz,zv) ) continue;
3266 if (TMath::Abs(zv-z3)>cuts[2]){
3267 FollowProlongation(*track, TMath::Max(i2-40,0));
3268 if ( !track->GetZAt(0.,bz,zv) ) continue;
3269 if (TMath::Abs(zv-z3)>cuts[2] &&(track->GetNumberOfClusters() > track->GetNFoundable()*0.7)){
3270 // make seed without constrain
3271 AliTPCseed * track2 = MakeSeed(track,0.2,0.5,1.);
3272 FollowProlongation(*track2, i2,1);
3273 track2->SetBConstrain(kFALSE);
3274 track2->SetSeedType(1);
3275 arr->AddLast(track2);
3276 MarkSeedFree( seed ); seed = 0;
3280 MarkSeedFree( seed ); seed = 0;
3287 track->SetSeedType(0);
3288 arr->AddLast(track); // note, track is seed, don't free the seed
3289 seed = new( NextFreeSeed() ) AliTPCseed;
3290 seed->SetPoolID(fLastSeedID);
3292 // don't consider other combinations
3293 if (track->GetNumberOfClusters() > track->GetNFoundable()*0.8)
3299 Info("MakeSeeds3","\nSeeding statistic:\t%d\t%d\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin,nout1,nout2);
3301 if (seed) MarkSeedFree( seed );
3305 void AliTPCtrackerMI::MakeSeeds5(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, Float_t cuts[4],
3310 //-----------------------------------------------------------------
3311 // This function creates track seeds.
3312 //-----------------------------------------------------------------
3313 // cuts[0] - fP4 cut
3314 // cuts[1] - tan(phi) cut
3315 // cuts[2] - zvertex cut
3316 // cuts[3] - fP3 cut
3326 Double_t x[5], c[15];
3328 // make temporary seed
3329 AliTPCseed * seed = new( NextFreeSeed() ) AliTPCseed;
3330 seed->SetPoolID(fLastSeedID);
3331 Double_t alpha=fOuterSec->GetAlpha(), shift=fOuterSec->GetAlphaShift();
3332 // Double_t cs=cos(alpha), sn=sin(alpha);
3337 Double_t x1 = GetXrow(i1-1);
3338 const AliTPCtrackerRow& kr1=GetRow(sec,i1-1);
3339 Double_t y1max = GetMaxY(i1-1)-kr1.GetDeadZone()-1.5;
3341 Double_t x1p = GetXrow(i1);
3342 const AliTPCtrackerRow& kr1p=GetRow(sec,i1);
3344 Double_t x1m = GetXrow(i1-2);
3345 const AliTPCtrackerRow& kr1m=GetRow(sec,i1-2);
3348 //last 3 padrow for seeding
3349 AliTPCtrackerRow& kr3 = GetRow((sec+fkNOS)%fkNOS,i1-7);
3350 Double_t x3 = GetXrow(i1-7);
3351 // Double_t y3max= GetMaxY(i1-7)-kr3.fDeadZone-1.5;
3353 AliTPCtrackerRow& kr3p = GetRow((sec+fkNOS)%fkNOS,i1-6);
3354 Double_t x3p = GetXrow(i1-6);
3356 AliTPCtrackerRow& kr3m = GetRow((sec+fkNOS)%fkNOS,i1-8);
3357 Double_t x3m = GetXrow(i1-8);
3362 Int_t im = i1-4; //middle pad row index
3363 Double_t xm = GetXrow(im); // radius of middle pad-row
3364 const AliTPCtrackerRow& krm=GetRow(sec,im); //middle pad -row
3365 // Double_t ymmax = GetMaxY(im)-kr1.fDeadZone-1.5;
3368 Double_t deltax = x1-x3;
3369 Double_t dymax = deltax*cuts[1];
3370 Double_t dzmax = deltax*cuts[3];
3372 // loop over clusters
3373 for (Int_t is=0; is < kr1; is++) {
3375 if (kr1[is]->IsUsed(10)) continue;
3376 Double_t y1=kr1[is]->GetY(), z1=kr1[is]->GetZ();
3378 if (deltay>0 && TMath::Abs(y1max-TMath::Abs(y1))> deltay ) continue; // seed only at the edge
3380 Int_t index1 = TMath::Max(kr3.Find(z1-dzmax)-1,0);
3381 Int_t index2 = TMath::Min(kr3.Find(z1+dzmax)+1,kr3);
3387 for (Int_t js=index1; js < index2; js++) {
3388 const AliTPCclusterMI *kcl = kr3[js];
3389 if (kcl->IsUsed(10)) continue;
3391 // apply angular cuts
3392 if (TMath::Abs(y1-y3)>dymax) continue;
3395 if (TMath::Abs(z1-z3)>dzmax) continue;
3397 Double_t angley = (y1-y3)/(x1-x3);
3398 Double_t anglez = (z1-z3)/(x1-x3);
3400 Double_t erry = TMath::Abs(angley)*(x1-x1m)*0.5+0.5;
3401 Double_t errz = TMath::Abs(anglez)*(x1-x1m)*0.5+0.5;
3403 Double_t yyym = angley*(xm-x1)+y1;
3404 Double_t zzzm = anglez*(xm-x1)+z1;
3406 const AliTPCclusterMI *kcm = krm.FindNearest2(yyym,zzzm,erry,errz,index);
3408 if (kcm->IsUsed(10)) continue;
3410 erry = TMath::Abs(angley)*(x1-x1m)*0.4+0.5;
3411 errz = TMath::Abs(anglez)*(x1-x1m)*0.4+0.5;
3418 // look around first
3419 const AliTPCclusterMI *kc1m = kr1m.FindNearest2(angley*(x1m-x1)+y1,
3425 if (kc1m->IsUsed(10)) used++;
3427 const AliTPCclusterMI *kc1p = kr1p.FindNearest2(angley*(x1p-x1)+y1,
3433 if (kc1p->IsUsed(10)) used++;
3435 if (used>1) continue;
3436 if (found<1) continue;
3440 const AliTPCclusterMI *kc3m = kr3m.FindNearest2(angley*(x3m-x3)+y3,
3446 if (kc3m->IsUsed(10)) used++;
3450 const AliTPCclusterMI *kc3p = kr3p.FindNearest2(angley*(x3p-x3)+y3,
3456 if (kc3p->IsUsed(10)) used++;
3460 if (used>1) continue;
3461 if (found<3) continue;
3471 x[4]=F1(x1,y1,x2,y2,x3,y3);
3472 //if (TMath::Abs(x[4]) >= cuts[0]) continue;
3475 x[2]=F2(x1,y1,x2,y2,x3,y3);
3478 x[3]=F3n(x1,y1,x2,y2,z1,z2,x[4]);
3479 //if (TMath::Abs(x[3]) > cuts[3]) continue;
3483 Double_t sy1=0.1, sz1=0.1;
3484 Double_t sy2=0.1, sz2=0.1;
3485 Double_t sy3=0.1, sy=0.1, sz=0.1;
3487 Double_t f40=(F1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
3488 Double_t f42=(F1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
3489 Double_t f43=(F1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
3490 Double_t f20=(F2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
3491 Double_t f22=(F2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
3492 Double_t f23=(F2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
3494 Double_t f30=(F3(x1,y1+sy,x2,y2,z1,z2)-x[3])/sy;
3495 Double_t f31=(F3(x1,y1,x2,y2,z1+sz,z2)-x[3])/sz;
3496 Double_t f32=(F3(x1,y1,x2,y2+sy,z1,z2)-x[3])/sy;
3497 Double_t f34=(F3(x1,y1,x2,y2,z1,z2+sz)-x[3])/sz;
3501 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3502 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3503 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3504 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3505 c[13]=f30*sy1*f40+f32*sy2*f42;
3506 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3508 // if (!BuildSeed(kr1[is],kcl,0,x1,x2,x3,x,c)) continue;
3510 index=kr1.GetIndex(is);
3511 if (seed) {MarkSeedFree( seed ); seed = 0;}
3512 AliTPCseed *track = seed = new( NextFreeSeed() ) AliTPCseed(x1, sec*alpha+shift, x, c, index);
3513 seed->SetPoolID(fLastSeedID);
3515 track->SetIsSeeding(kTRUE);
3518 FollowProlongation(*track, i1-7,1);
3519 if (track->GetNumberOfClusters() < track->GetNFoundable()*0.75 ||
3520 track->GetNShared()>0.6*track->GetNumberOfClusters() || ( track->GetSigmaY2()+ track->GetSigmaZ2())>0.6){
3521 MarkSeedFree( seed ); seed = 0;
3527 FollowProlongation(*track, i2,1);
3528 track->SetBConstrain(0);
3529 track->SetLastPoint(i1+fInnerSec->GetNRows()); // first cluster in track position
3530 track->SetFirstPoint(track->GetLastPoint());
3532 if (track->GetNumberOfClusters()<(i1-i2)*0.5 ||
3533 track->GetNumberOfClusters()<track->GetNFoundable()*0.7 ||
3534 track->GetNShared()>2. || track->GetChi2()/track->GetNumberOfClusters()>6 || ( track->GetSigmaY2()+ track->GetSigmaZ2())>0.5 ) {
3535 MarkSeedFree( seed ); seed = 0;
3540 FollowProlongation(*track, TMath::Max(i2-10,0),1);
3541 AliTPCseed * track2 = MakeSeed(track,0.2,0.5,0.9);
3542 FollowProlongation(*track2, i2,1);
3543 track2->SetBConstrain(kFALSE);
3544 track2->SetSeedType(4);
3545 arr->AddLast(track2);
3546 MarkSeedFree( seed ); seed = 0;
3550 //arr->AddLast(track);
3551 //seed = new AliTPCseed;
3557 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);
3559 if (seed) MarkSeedFree(seed);
3563 //_____________________________________________________________________________
3564 void AliTPCtrackerMI::MakeSeeds2(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, Float_t */*cuts[4]*/,
3565 Float_t deltay, Bool_t /*bconstrain*/) {
3566 //-----------------------------------------------------------------
3567 // This function creates track seeds - without vertex constraint
3568 //-----------------------------------------------------------------
3569 // cuts[0] - fP4 cut - not applied
3570 // cuts[1] - tan(phi) cut
3571 // cuts[2] - zvertex cut - not applied
3572 // cuts[3] - fP3 cut
3582 Double_t alpha=fOuterSec->GetAlpha(), shift=fOuterSec->GetAlphaShift();
3583 // Double_t cs=cos(alpha), sn=sin(alpha);
3584 Int_t row0 = (i1+i2)/2;
3585 Int_t drow = (i1-i2)/2;
3586 const AliTPCtrackerRow& kr0=fSectors[sec][row0];
3587 AliTPCtrackerRow * kr=0;
3589 AliTPCpolyTrack polytrack;
3590 Int_t nclusters=fSectors[sec][row0];
3591 AliTPCseed * seed = new( NextFreeSeed() ) AliTPCseed;
3592 seed->SetPoolID(fLastSeedID);
3597 for (Int_t is=0; is < nclusters; is++) { //LOOP over clusters
3599 Int_t nfoundable =0;
3600 for (Int_t iter =1; iter<2; iter++){ //iterations
3601 const AliTPCtrackerRow& krm=fSectors[sec][row0-iter];
3602 const AliTPCtrackerRow& krp=fSectors[sec][row0+iter];
3603 const AliTPCclusterMI * cl= kr0[is];
3605 if (cl->IsUsed(10)) {
3611 Double_t x = kr0.GetX();
3612 // Initialization of the polytrack
3617 Double_t y0= cl->GetY();
3618 Double_t z0= cl->GetZ();
3622 Double_t ymax = fSectors->GetMaxY(row0)-kr0.GetDeadZone()-1.5;
3623 if (deltay>0 && TMath::Abs(ymax-TMath::Abs(y0))> deltay ) continue; // seed only at the edge
3625 erry = (0.5)*cl->GetSigmaY2()/TMath::Sqrt(cl->GetQ())*6;
3626 errz = (0.5)*cl->GetSigmaZ2()/TMath::Sqrt(cl->GetQ())*6;
3627 polytrack.AddPoint(x,y0,z0,erry, errz);
3630 if (cl->IsUsed(10)) sumused++;
3633 Float_t roady = (5*TMath::Sqrt(cl->GetSigmaY2()+0.2)+1.)*iter;
3634 Float_t roadz = (5*TMath::Sqrt(cl->GetSigmaZ2()+0.2)+1.)*iter;
3637 AliTPCclusterMI * cl1 = krm.FindNearest(y0,z0,roady,roadz);
3638 if (cl1 && TMath::Abs(ymax-TMath::Abs(y0))) {
3639 erry = (0.5)*cl1->GetSigmaY2()/TMath::Sqrt(cl1->GetQ())*3;
3640 errz = (0.5)*cl1->GetSigmaZ2()/TMath::Sqrt(cl1->GetQ())*3;
3641 if (cl1->IsUsed(10)) sumused++;
3642 polytrack.AddPoint(x,cl1->GetY(),cl1->GetZ(),erry,errz);
3646 AliTPCclusterMI * cl2 = krp.FindNearest(y0,z0,roady,roadz);
3648 erry = (0.5)*cl2->GetSigmaY2()/TMath::Sqrt(cl2->GetQ())*3;
3649 errz = (0.5)*cl2->GetSigmaZ2()/TMath::Sqrt(cl2->GetQ())*3;
3650 if (cl2->IsUsed(10)) sumused++;
3651 polytrack.AddPoint(x,cl2->GetY(),cl2->GetZ(),erry,errz);
3654 if (sumused>0) continue;
3656 polytrack.UpdateParameters();
3662 nfoundable = polytrack.GetN();
3663 nfound = nfoundable;
3665 for (Int_t ddrow = iter+1; ddrow<drow;ddrow++){
3666 Float_t maxdist = 0.8*(1.+3./(ddrow));
3667 for (Int_t delta = -1;delta<=1;delta+=2){
3668 Int_t row = row0+ddrow*delta;
3669 kr = &(fSectors[sec][row]);
3670 Double_t xn = kr->GetX();
3671 Double_t ymax1 = fSectors->GetMaxY(row)-kr->GetDeadZone()-1.5;
3672 polytrack.GetFitPoint(xn,yn,zn);
3673 if (TMath::Abs(yn)>ymax1) continue;
3675 AliTPCclusterMI * cln = kr->FindNearest(yn,zn,roady,roadz);
3677 Float_t dist = TMath::Sqrt( (yn-cln->GetY())*(yn-cln->GetY())+(zn-cln->GetZ())*(zn-cln->GetZ()));
3680 erry = (dist+0.3)*cln->GetSigmaY2()/TMath::Sqrt(cln->GetQ())*(1.+1./(ddrow));
3681 errz = (dist+0.3)*cln->GetSigmaZ2()/TMath::Sqrt(cln->GetQ())*(1.+1./(ddrow));
3682 if (cln->IsUsed(10)) {
3683 // printf("used\n");
3691 polytrack.AddPoint(xn,cln->GetY(),cln->GetZ(),erry, errz);
3696 if ( (sumused>3) || (sumused>0.5*nfound) || (nfound<0.6*nfoundable)) break;
3697 polytrack.UpdateParameters();
3700 if ( (sumused>3) || (sumused>0.5*nfound)) {
3701 //printf("sumused %d\n",sumused);
3706 polytrack.GetFitDerivation(kr0.GetX(),dy,dz);
3707 AliTPCpolyTrack track2;
3709 polytrack.Refit(track2,0.5+TMath::Abs(dy)*0.3,0.4+TMath::Abs(dz)*0.3);
3710 if (track2.GetN()<0.5*nfoundable) continue;
3713 if ((nfound>0.6*nfoundable) &&( nfoundable>0.4*(i1-i2))) {
3715 // test seed with and without constrain
3716 for (Int_t constrain=0; constrain<=0;constrain++){
3717 // add polytrack candidate
3719 Double_t x[5], c[15];
3720 Double_t x1,x2,x3,y1,y2,y3,z1,z2,z3;
3721 track2.GetBoundaries(x3,x1);
3723 track2.GetFitPoint(x1,y1,z1);
3724 track2.GetFitPoint(x2,y2,z2);
3725 track2.GetFitPoint(x3,y3,z3);
3727 //is track pointing to the vertex ?
3730 polytrack.GetFitPoint(x0,y0,z0);
3743 x[4]=F1(x1,y1,x2,y2,x3,y3);
3745 // if (TMath::Abs(x[4]) >= cuts[0]) continue; //
3746 x[2]=F2(x1,y1,x2,y2,x3,y3);
3748 //if (TMath::Abs(x[4]*x1-x[2]) >= cuts[1]) continue;
3749 //x[3]=F3(x1,y1,x2,y2,z1,z2);
3750 x[3]=F3n(x1,y1,x3,y3,z1,z3,x[4]);
3751 //if (TMath::Abs(x[3]) > cuts[3]) continue;
3754 Double_t sy =0.1, sz =0.1;
3755 Double_t sy1=0.02, sz1=0.02;
3756 Double_t sy2=0.02, sz2=0.02;
3760 sy3=25000*x[4]*x[4]+0.1, sy=0.1, sz=0.1;
3763 Double_t f40=(F1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
3764 Double_t f42=(F1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
3765 Double_t f43=(F1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
3766 Double_t f20=(F2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
3767 Double_t f22=(F2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
3768 Double_t f23=(F2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
3770 Double_t f30=(F3(x1,y1+sy,x3,y3,z1,z3)-x[3])/sy;
3771 Double_t f31=(F3(x1,y1,x3,y3,z1+sz,z3)-x[3])/sz;
3772 Double_t f32=(F3(x1,y1,x3,y3+sy,z1,z3)-x[3])/sy;
3773 Double_t f34=(F3(x1,y1,x3,y3,z1,z3+sz)-x[3])/sz;
3778 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3779 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3780 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3781 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3782 c[13]=f30*sy1*f40+f32*sy2*f42;
3783 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3785 //Int_t row1 = fSectors->GetRowNumber(x1);
3786 Int_t row1 = GetRowNumber(x1);
3790 if (seed) {MarkSeedFree( seed ); seed = 0;}
3791 AliTPCseed *track = seed = new( NextFreeSeed() ) AliTPCseed(x1,sec*alpha+shift,x,c,index);
3792 seed->SetPoolID(fLastSeedID);
3793 track->SetIsSeeding(kTRUE);
3794 Int_t rc=FollowProlongation(*track, i2);
3795 if (constrain) track->SetBConstrain(1);
3797 track->SetBConstrain(0);
3798 track->SetLastPoint(row1+fInnerSec->GetNRows()); // first cluster in track position
3799 track->SetFirstPoint(track->GetLastPoint());
3801 if (rc==0 || track->GetNumberOfClusters()<(i1-i2)*0.5 ||
3802 track->GetNumberOfClusters() < track->GetNFoundable()*0.6 ||
3803 track->GetNShared()>0.4*track->GetNumberOfClusters()) {
3804 MarkSeedFree( seed ); seed = 0;
3807 arr->AddLast(track); // track IS seed, don't free seed
3808 seed = new( NextFreeSeed() ) AliTPCseed;
3809 seed->SetPoolID(fLastSeedID);
3813 } // if accepted seed
3816 Info("MakeSeeds2","\nSeeding statiistic:\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin3);
3818 if (seed) MarkSeedFree( seed );
3822 AliTPCseed *AliTPCtrackerMI::MakeSeed(AliTPCseed *const track, Float_t r0, Float_t r1, Float_t r2)
3826 //reseed using track points
3827 Int_t p0 = int(r0*track->GetNumberOfClusters()); // point 0
3828 Int_t p1 = int(r1*track->GetNumberOfClusters());
3829 Int_t p2 = int(r2*track->GetNumberOfClusters()); // last point
3831 Double_t x0[3],x1[3],x2[3];
3832 for (Int_t i=0;i<3;i++){
3838 // find track position at given ratio of the length
3839 Int_t sec0=0, sec1=0, sec2=0;
3842 for (Int_t i=0;i<160;i++){
3843 if (track->GetClusterPointer(i)){
3845 AliTPCTrackerPoint *trpoint =track->GetTrackPoint(i);
3846 if ( (index<p0) || x0[0]<0 ){
3847 if (trpoint->GetX()>1){
3848 clindex = track->GetClusterIndex2(i);
3850 x0[0] = trpoint->GetX();
3851 x0[1] = trpoint->GetY();
3852 x0[2] = trpoint->GetZ();
3853 sec0 = ((clindex&0xff000000)>>24)%18;
3858 if ( (index<p1) &&(trpoint->GetX()>1)){
3859 clindex = track->GetClusterIndex2(i);
3861 x1[0] = trpoint->GetX();
3862 x1[1] = trpoint->GetY();
3863 x1[2] = trpoint->GetZ();
3864 sec1 = ((clindex&0xff000000)>>24)%18;
3867 if ( (index<p2) &&(trpoint->GetX()>1)){
3868 clindex = track->GetClusterIndex2(i);
3870 x2[0] = trpoint->GetX();
3871 x2[1] = trpoint->GetY();
3872 x2[2] = trpoint->GetZ();
3873 sec2 = ((clindex&0xff000000)>>24)%18;
3880 Double_t alpha, cs,sn, xx2,yy2;
3882 alpha = (sec1-sec2)*fSectors->GetAlpha();
3883 cs = TMath::Cos(alpha);
3884 sn = TMath::Sin(alpha);
3885 xx2= x1[0]*cs-x1[1]*sn;
3886 yy2= x1[0]*sn+x1[1]*cs;
3890 alpha = (sec0-sec2)*fSectors->GetAlpha();
3891 cs = TMath::Cos(alpha);
3892 sn = TMath::Sin(alpha);
3893 xx2= x0[0]*cs-x0[1]*sn;
3894 yy2= x0[0]*sn+x0[1]*cs;
3900 Double_t x[5],c[15];
3904 x[4]=F1(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]);
3905 // if (x[4]>1) return 0;
3906 x[2]=F2(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]);
3907 x[3]=F3n(x2[0],x2[1],x0[0],x0[1],x2[2],x0[2],x[4]);
3908 //if (TMath::Abs(x[3]) > 2.2) return 0;
3909 //if (TMath::Abs(x[2]) > 1.99) return 0;
3911 Double_t sy =0.1, sz =0.1;
3913 Double_t sy1=0.02+track->GetSigmaY2(), sz1=0.02+track->GetSigmaZ2();
3914 Double_t sy2=0.01+track->GetSigmaY2(), sz2=0.01+track->GetSigmaZ2();
3915 Double_t sy3=0.01+track->GetSigmaY2();
3917 Double_t f40=(F1(x2[0],x2[1]+sy,x1[0],x1[1],x0[0],x0[1])-x[4])/sy;
3918 Double_t f42=(F1(x2[0],x2[1],x1[0],x1[1]+sy,x0[0],x0[1])-x[4])/sy;
3919 Double_t f43=(F1(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]+sy)-x[4])/sy;
3920 Double_t f20=(F2(x2[0],x2[1]+sy,x1[0],x1[1],x0[0],x0[1])-x[2])/sy;
3921 Double_t f22=(F2(x2[0],x2[1],x1[0],x1[1]+sy,x0[0],x0[1])-x[2])/sy;
3922 Double_t f23=(F2(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]+sy)-x[2])/sy;
3924 Double_t f30=(F3(x2[0],x2[1]+sy,x0[0],x0[1],x2[2],x0[2])-x[3])/sy;
3925 Double_t f31=(F3(x2[0],x2[1],x0[0],x0[1],x2[2]+sz,x0[2])-x[3])/sz;
3926 Double_t f32=(F3(x2[0],x2[1],x0[0],x0[1]+sy,x2[2],x0[2])-x[3])/sy;
3927 Double_t f34=(F3(x2[0],x2[1],x0[0],x0[1],x2[2],x0[2]+sz)-x[3])/sz;
3932 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3933 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3934 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3935 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3936 c[13]=f30*sy1*f40+f32*sy2*f42;
3937 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3939 // Int_t row1 = fSectors->GetRowNumber(x2[0]);
3940 AliTPCseed *seed = new( NextFreeSeed() ) AliTPCseed(x2[0], sec2*fSectors->GetAlpha()+fSectors->GetAlphaShift(), x, c, 0);
3941 seed->SetPoolID(fLastSeedID);
3942 // Double_t y0,z0,y1,z1, y2,z2;
3943 //seed->GetProlongation(x0[0],y0,z0);
3944 // seed->GetProlongation(x1[0],y1,z1);
3945 //seed->GetProlongation(x2[0],y2,z2);
3947 seed->SetLastPoint(pp2);
3948 seed->SetFirstPoint(pp2);
3955 AliTPCseed *AliTPCtrackerMI::ReSeed(const AliTPCseed *track, Float_t r0, Float_t r1, Float_t r2)
3959 //reseed using founded clusters
3961 // Find the number of clusters
3962 Int_t nclusters = 0;
3963 for (Int_t irow=0;irow<160;irow++){
3964 if (track->GetClusterIndex(irow)>0) nclusters++;
3968 ipos[0] = TMath::Max(int(r0*nclusters),0); // point 0 cluster
3969 ipos[1] = TMath::Min(int(r1*nclusters),nclusters-1); //
3970 ipos[2] = TMath::Min(int(r2*nclusters),nclusters-1); // last point
3973 Double_t xyz[3][3]={{0}};
3974 Int_t row[3]={0},sec[3]={0,0,0};
3976 // find track row position at given ratio of the length
3978 for (Int_t irow=0;irow<160;irow++){
3979 if (track->GetClusterIndex2(irow)<0) continue;
3981 for (Int_t ipoint=0;ipoint<3;ipoint++){
3982 if (index<=ipos[ipoint]) row[ipoint] = irow;
3986 //Get cluster and sector position
3987 for (Int_t ipoint=0;ipoint<3;ipoint++){
3988 Int_t clindex = track->GetClusterIndex2(row[ipoint]);
3989 AliTPCclusterMI * cl = GetClusterMI(clindex);
3992 // AliTPCclusterMI * cl = GetClusterMI(clindex);
3995 sec[ipoint] = ((clindex&0xff000000)>>24)%18;
3996 xyz[ipoint][0] = GetXrow(row[ipoint]);
3997 xyz[ipoint][1] = cl->GetY();
3998 xyz[ipoint][2] = cl->GetZ();
4002 // Calculate seed state vector and covariance matrix
4004 Double_t alpha, cs,sn, xx2,yy2;
4006 alpha = (sec[1]-sec[2])*fSectors->GetAlpha();
4007 cs = TMath::Cos(alpha);
4008 sn = TMath::Sin(alpha);
4009 xx2= xyz[1][0]*cs-xyz[1][1]*sn;
4010 yy2= xyz[1][0]*sn+xyz[1][1]*cs;
4014 alpha = (sec[0]-sec[2])*fSectors->GetAlpha();
4015 cs = TMath::Cos(alpha);
4016 sn = TMath::Sin(alpha);
4017 xx2= xyz[0][0]*cs-xyz[0][1]*sn;
4018 yy2= xyz[0][0]*sn+xyz[0][1]*cs;
4024 Double_t x[5],c[15];
4028 x[4]=F1(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
4029 x[2]=F2(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
4030 x[3]=F3n(xyz[2][0],xyz[2][1],xyz[0][0],xyz[0][1],xyz[2][2],xyz[0][2],x[4]);
4032 Double_t sy =0.1, sz =0.1;
4034 Double_t sy1=0.2, sz1=0.2;
4035 Double_t sy2=0.2, sz2=0.2;
4038 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;
4039 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;
4040 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;
4041 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;
4042 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;
4043 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;
4045 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;
4046 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;
4047 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;
4048 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;
4053 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
4054 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
4055 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
4056 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
4057 c[13]=f30*sy1*f40+f32*sy2*f42;
4058 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
4060 // Int_t row1 = fSectors->GetRowNumber(xyz[2][0]);
4061 AliTPCseed *seed=new( NextFreeSeed() ) AliTPCseed(xyz[2][0], sec[2]*fSectors->GetAlpha()+fSectors->GetAlphaShift(), x, c, 0);
4062 seed->SetPoolID(fLastSeedID);
4063 seed->SetLastPoint(row[2]);
4064 seed->SetFirstPoint(row[2]);
4069 AliTPCseed *AliTPCtrackerMI::ReSeed(AliTPCseed *track,Int_t r0, Bool_t forward)
4073 //reseed using founded clusters
4076 Int_t row[3]={0,0,0};
4077 Int_t sec[3]={0,0,0};
4079 // forward direction
4081 for (Int_t irow=r0;irow<160;irow++){
4082 if (track->GetClusterIndex(irow)>0){
4087 for (Int_t irow=160;irow>r0;irow--){
4088 if (track->GetClusterIndex(irow)>0){
4093 for (Int_t irow=row[2]-15;irow>row[0];irow--){
4094 if (track->GetClusterIndex(irow)>0){
4102 for (Int_t irow=0;irow<r0;irow++){
4103 if (track->GetClusterIndex(irow)>0){
4108 for (Int_t irow=r0;irow>0;irow--){
4109 if (track->GetClusterIndex(irow)>0){
4114 for (Int_t irow=row[2]-15;irow>row[0];irow--){
4115 if (track->GetClusterIndex(irow)>0){
4122 if ((row[2]-row[0])<20) return 0;
4123 if (row[1]==0) return 0;
4126 //Get cluster and sector position
4127 for (Int_t ipoint=0;ipoint<3;ipoint++){
4128 Int_t clindex = track->GetClusterIndex2(row[ipoint]);
4129 AliTPCclusterMI * cl = GetClusterMI(clindex);
4132 // AliTPCclusterMI * cl = GetClusterMI(clindex);
4135 sec[ipoint] = ((clindex&0xff000000)>>24)%18;
4136 xyz[ipoint][0] = GetXrow(row[ipoint]);
4137 AliTPCTrackerPoint * point = track->GetTrackPoint(row[ipoint]);
4138 if (point&&ipoint<2){
4140 xyz[ipoint][1] = point->GetY();
4141 xyz[ipoint][2] = point->GetZ();
4144 xyz[ipoint][1] = cl->GetY();
4145 xyz[ipoint][2] = cl->GetZ();
4152 // Calculate seed state vector and covariance matrix
4154 Double_t alpha, cs,sn, xx2,yy2;
4156 alpha = (sec[1]-sec[2])*fSectors->GetAlpha();
4157 cs = TMath::Cos(alpha);
4158 sn = TMath::Sin(alpha);
4159 xx2= xyz[1][0]*cs-xyz[1][1]*sn;
4160 yy2= xyz[1][0]*sn+xyz[1][1]*cs;
4164 alpha = (sec[0]-sec[2])*fSectors->GetAlpha();
4165 cs = TMath::Cos(alpha);
4166 sn = TMath::Sin(alpha);
4167 xx2= xyz[0][0]*cs-xyz[0][1]*sn;
4168 yy2= xyz[0][0]*sn+xyz[0][1]*cs;
4174 Double_t x[5],c[15];
4178 x[4]=F1(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
4179 x[2]=F2(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
4180 x[3]=F3n(xyz[2][0],xyz[2][1],xyz[0][0],xyz[0][1],xyz[2][2],xyz[0][2],x[4]);
4182 Double_t sy =0.1, sz =0.1;
4184 Double_t sy1=0.2, sz1=0.2;
4185 Double_t sy2=0.2, sz2=0.2;
4188 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;
4189 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;
4190 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;
4191 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;
4192 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;
4193 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;
4195 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;
4196 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;
4197 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;
4198 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;
4203 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
4204 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
4205 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
4206 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
4207 c[13]=f30*sy1*f40+f32*sy2*f42;
4208 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
4210 // Int_t row1 = fSectors->GetRowNumber(xyz[2][0]);
4211 AliTPCseed *seed=new( NextFreeSeed() ) AliTPCseed(xyz[2][0], sec[2]*fSectors->GetAlpha()+fSectors->GetAlphaShift(), x, c, 0);
4212 seed->SetPoolID(fLastSeedID);
4213 seed->SetLastPoint(row[2]);
4214 seed->SetFirstPoint(row[2]);
4215 for (Int_t i=row[0];i<row[2];i++){
4216 seed->SetClusterIndex(i, track->GetClusterIndex(i));
4224 void AliTPCtrackerMI::FindMultiMC(const TObjArray * array, AliESDEvent */*esd*/, Int_t iter)
4227 // find multi tracks - THIS FUNCTION IS ONLY FOR DEBUG PURPOSES
4229 // Use AliTPCReconstructor::StreamLevel()>2 if you want to tune parameters - cuts
4231 // Two reasons to have multiple find tracks
4232 // 1. Curling tracks can be find more than once
4233 // 2. Splitted tracks
4234 // a.) Multiple seeding to increase tracking efficiency - (~ 100% reached)
4235 // b.) Edge effect on the sector boundaries
4238 // Algorithm done in 2 phases - because of CPU consumption
4239 // it is n^2 algorithm - for lead-lead 20000x20000 combination are investigated
4241 // Algorihm for curling tracks sign:
4242 // 1 phase -makes a very rough fast cuts to minimize combinatorics
4243 // a.) opposite sign
4244 // b.) one of the tracks - not pointing to the primary vertex -
4245 // c.) delta tan(theta)
4247 // 2 phase - calculates DCA between tracks - time consument
4252 // General cuts - for splitted tracks and for curling tracks
4254 const Float_t kMaxdPhi = 0.2; // maximal distance in phi
4256 // Curling tracks cuts
4261 Int_t nentries = array->GetEntriesFast();
4262 AliHelix *helixes = new AliHelix[nentries];
4263 Float_t *xm = new Float_t[nentries];
4264 Float_t *dz0 = new Float_t[nentries];
4265 Float_t *dz1 = new Float_t[nentries];
4271 // Find track COG in x direction - point with best defined parameters
4273 for (Int_t i=0;i<nentries;i++){
4274 AliTPCseed* track = (AliTPCseed*)array->At(i);
4275 if (!track) continue;
4276 track->SetCircular(0);
4277 new (&helixes[i]) AliHelix(*track);
4281 track->GetDZ(GetX(),GetY(),GetZ(),GetBz(),dz);
4284 for (Int_t icl=0; icl<160; icl++){
4285 AliTPCclusterMI * cl = track->GetClusterPointer(icl);
4291 if (ncl>0) xm[i]/=Float_t(ncl);
4294 for (Int_t i0=0;i0<nentries;i0++){
4295 AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
4296 if (!track0) continue;
4297 Float_t xc0 = helixes[i0].GetHelix(6);
4298 Float_t yc0 = helixes[i0].GetHelix(7);
4299 Float_t r0 = helixes[i0].GetHelix(8);
4300 Float_t rc0 = TMath::Sqrt(xc0*xc0+yc0*yc0);
4301 Float_t fi0 = TMath::ATan2(yc0,xc0);
4303 for (Int_t i1=i0+1;i1<nentries;i1++){
4304 AliTPCseed * track1 = (AliTPCseed*)array->At(i1);
4305 if (!track1) continue;
4306 Int_t lab0=track0->GetLabel();
4307 Int_t lab1=track1->GetLabel();
4308 if (TMath::Abs(lab0)!=TMath::Abs(lab1)) continue;
4310 Float_t xc1 = helixes[i1].GetHelix(6);
4311 Float_t yc1 = helixes[i1].GetHelix(7);
4312 Float_t r1 = helixes[i1].GetHelix(8);
4313 Float_t rc1 = TMath::Sqrt(xc1*xc1+yc1*yc1);
4314 Float_t fi1 = TMath::ATan2(yc1,xc1);
4316 Float_t dfi = fi0-fi1;
4319 if (dfi>1.5*TMath::Pi()) dfi-=TMath::Pi(); // take care about edge effect
4320 if (dfi<-1.5*TMath::Pi()) dfi+=TMath::Pi(); //
4321 if (TMath::Abs(dfi)>kMaxdPhi&&helixes[i0].GetHelix(4)*helixes[i1].GetHelix(4)<0){
4323 // if short tracks with undefined sign
4324 fi1 = -TMath::ATan2(yc1,-xc1);
4327 Float_t dtheta = TMath::Abs(track0->GetTgl()-track1->GetTgl())<TMath::Abs(track0->GetTgl()+track1->GetTgl())? track0->GetTgl()-track1->GetTgl():track0->GetTgl()+track1->GetTgl();
4330 // debug stream to tune "fast cuts"
4332 Double_t dist[3]; // distance at X
4333 Double_t mdist[3]={0,0,0}; // mean distance X+-40cm
4334 track0->GetDistance(track1,0.5*(xm[i0]+xm[i1])-40.,dist,AliTracker::GetBz());
4335 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[i]);
4336 track0->GetDistance(track1,0.5*(xm[i0]+xm[i1])+40.,dist,AliTracker::GetBz());
4337 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[i]);
4338 track0->GetDistance(track1,0.5*(xm[i0]+xm[i1]),dist,AliTracker::GetBz());
4339 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[i]);
4340 for (Int_t i=0;i<3;i++) mdist[i]*=0.33333;
4344 for (Int_t icl=0; icl<160; icl++){
4345 AliTPCclusterMI * cl0 = track0->GetClusterPointer(icl);
4346 AliTPCclusterMI * cl1 = track1->GetClusterPointer(icl);
4349 if (cl0==cl1) sums++;
4353 if (AliTPCReconstructor::StreamLevel()>5) {
4354 TTreeSRedirector &cstream = *fDebugStreamer;
4359 "Tr0.="<<track0<< // seed0
4360 "Tr1.="<<track1<< // seed1
4361 "h0.="<<&helixes[i0]<<
4362 "h1.="<<&helixes[i1]<<
4364 "sum="<<sum<< //the sum of rows with cl in both
4365 "sums="<<sums<< //the sum of shared clusters
4366 "xm0="<<xm[i0]<< // the center of track
4367 "xm1="<<xm[i1]<< // the x center of track
4368 // General cut variables
4369 "dfi="<<dfi<< // distance in fi angle
4370 "dtheta="<<dtheta<< // distance int theta angle
4376 "dist0="<<dist[0]<< //distance x
4377 "dist1="<<dist[1]<< //distance y
4378 "dist2="<<dist[2]<< //distance z
4379 "mdist0="<<mdist[0]<< //distance x
4380 "mdist1="<<mdist[1]<< //distance y
4381 "mdist2="<<mdist[2]<< //distance z
4397 if (AliTPCReconstructor::StreamLevel()>1) {
4398 AliInfo("Time for curling tracks removal DEBUGGING MC");
4405 void AliTPCtrackerMI::FindSplitted(TObjArray * array, AliESDEvent */*esd*/, Int_t /*iter*/){
4407 // Find Splitted tracks and remove the one with worst quality
4408 // Corresponding debug streamer to tune selections - "Splitted2"
4410 // 0. Sort tracks according quility
4411 // 1. Propagate the tracks to the reference radius
4412 // 2. Double_t loop to select close tracks (only to speed up process)
4413 // 3. Calculate cluster overlap ratio - and remove the track if bigger than a threshold
4414 // 4. Delete temporary parameters
4416 const Double_t xref=GetXrow(63); // reference radius -IROC/OROC boundary
4418 const Double_t kCutP1=10; // delta Z cut 10 cm
4419 const Double_t kCutP2=0.15; // delta snp(fi) cut 0.15
4420 const Double_t kCutP3=0.15; // delta tgl(theta) cut 0.15
4421 const Double_t kCutAlpha=0.15; // delta alpha cut
4422 Int_t firstpoint = 0;
4423 Int_t lastpoint = 160;
4425 Int_t nentries = array->GetEntriesFast();
4426 AliExternalTrackParam *params = new AliExternalTrackParam[nentries];
4432 //0. Sort tracks according quality
4433 //1. Propagate the ext. param to reference radius
4434 Int_t nseed = array->GetEntriesFast();
4435 if (nseed<=0) return;
4436 Float_t * quality = new Float_t[nseed];
4437 Int_t * indexes = new Int_t[nseed];
4438 for (Int_t i=0; i<nseed; i++) {
4439 AliTPCseed *pt=(AliTPCseed*)array->UncheckedAt(i);
4444 pt->UpdatePoints(); //select first last max dens points
4445 Float_t * points = pt->GetPoints();
4446 if (points[3]<0.8) quality[i] =-1;
4447 quality[i] = (points[2]-points[0])+pt->GetNumberOfClusters();
4448 //prefer high momenta tracks if overlaps
4449 quality[i] *= TMath::Sqrt(TMath::Abs(pt->Pt())+0.5);
4451 AliTracker::PropagateTrackToBxByBz(&(params[i]),xref,pt->GetMass(),5.,kTRUE);
4452 AliTracker::PropagateTrackToBxByBz(&(params[i]),xref,pt->GetMass(),1.,kTRUE);
4454 TMath::Sort(nseed,quality,indexes);
4456 // 3. Loop over pair of tracks
4458 for (Int_t i0=0; i0<nseed; i0++) {
4459 Int_t index0=indexes[i0];
4460 if (!(array->UncheckedAt(index0))) continue;
4461 AliTPCseed *s1 = (AliTPCseed*)array->UncheckedAt(index0);
4462 if (!s1->IsActive()) continue;
4463 AliExternalTrackParam &par0=params[index0];
4464 for (Int_t i1=i0+1; i1<nseed; i1++) {
4465 Int_t index1=indexes[i1];
4466 if (!(array->UncheckedAt(index1))) continue;
4467 AliTPCseed *s2 = (AliTPCseed*)array->UncheckedAt(index1);
4468 if (!s2->IsActive()) continue;
4469 if (s2->GetKinkIndexes()[0]!=0)
4470 if (s2->GetKinkIndexes()[0] == -s1->GetKinkIndexes()[0]) continue;
4471 AliExternalTrackParam &par1=params[index1];
4472 if (TMath::Abs(par0.GetParameter()[3]-par1.GetParameter()[3])>kCutP3) continue;
4473 if (TMath::Abs(par0.GetParameter()[1]-par1.GetParameter()[1])>kCutP1) continue;
4474 if (TMath::Abs(par0.GetParameter()[2]-par1.GetParameter()[2])>kCutP2) continue;
4475 Double_t dAlpha= TMath::Abs(par0.GetAlpha()-par1.GetAlpha());
4476 if (dAlpha>TMath::Pi()) dAlpha-=TMath::Pi();
4477 if (TMath::Abs(dAlpha)>kCutAlpha) continue;
4482 Int_t firstShared=lastpoint, lastShared=firstpoint;
4483 Int_t firstRow=lastpoint, lastRow=firstpoint;
4485 for (Int_t i=firstpoint;i<lastpoint;i++){
4486 if (s1->GetClusterIndex2(i)>0) nall0++;
4487 if (s2->GetClusterIndex2(i)>0) nall1++;
4488 if (s1->GetClusterIndex2(i)>0 && s2->GetClusterIndex2(i)>0) {
4489 if (i<firstRow) firstRow=i;
4490 if (i>lastRow) lastRow=i;
4492 if ( (s1->GetClusterIndex2(i))==(s2->GetClusterIndex2(i)) && s1->GetClusterIndex2(i)>0) {
4493 if (i<firstShared) firstShared=i;
4494 if (i>lastShared) lastShared=i;
4498 Double_t ratio0 = Float_t(sumShared)/Float_t(TMath::Min(nall0+1,nall1+1));
4499 Double_t ratio1 = Float_t(sumShared)/Float_t(TMath::Max(nall0+1,nall1+1));
4501 if( AliTPCReconstructor::StreamLevel()>1){
4502 TTreeSRedirector &cstream = *fDebugStreamer;
4503 Int_t n0=s1->GetNumberOfClusters();
4504 Int_t n1=s2->GetNumberOfClusters();
4505 Int_t n0F=s1->GetNFoundable();
4506 Int_t n1F=s2->GetNFoundable();
4507 Int_t lab0=s1->GetLabel();
4508 Int_t lab1=s2->GetLabel();
4510 cstream<<"Splitted2"<<
4511 "iter="<<fIteration<<
4512 "lab0="<<lab0<< // MC label if exist
4513 "lab1="<<lab1<< // MC label if exist
4516 "ratio0="<<ratio0<< // shared ratio
4517 "ratio1="<<ratio1<< // shared ratio
4518 "p0.="<<&par0<< // track parameters
4520 "s0.="<<s1<< // full seed
4522 "n0="<<n0<< // number of clusters track 0
4523 "n1="<<n1<< // number of clusters track 1
4524 "nall0="<<nall0<< // number of clusters track 0
4525 "nall1="<<nall1<< // number of clusters track 1
4526 "n0F="<<n0F<< // number of findable
4527 "n1F="<<n1F<< // number of findable
4528 "shared="<<sumShared<< // number of shared clusters
4529 "firstS="<<firstShared<< // first and the last shared row
4530 "lastS="<<lastShared<<
4531 "firstRow="<<firstRow<< // first and the last row with cluster
4532 "lastRow="<<lastRow<< //
4536 // remove track with lower quality
4538 if (ratio0>AliTPCReconstructor::GetRecoParam()->GetCutSharedClusters(0) ||
4539 ratio1>AliTPCReconstructor::GetRecoParam()->GetCutSharedClusters(1)){
4543 MarkSeedFree( array->RemoveAt(index1) );
4548 // 4. Delete temporary array
4558 void AliTPCtrackerMI::FindCurling(const TObjArray * array, AliESDEvent */*esd*/, Int_t iter)
4561 // find Curling tracks
4562 // Use AliTPCReconstructor::StreamLevel()>1 if you want to tune parameters - cuts
4565 // Algorithm done in 2 phases - because of CPU consumption
4566 // it is n^2 algorithm - for lead-lead 20000x20000 combination are investigated
4567 // see detal in MC part what can be used to cut
4571 const Float_t kMaxC = 400; // maximal curvature to of the track
4572 const Float_t kMaxdTheta = 0.15; // maximal distance in theta
4573 const Float_t kMaxdPhi = 0.15; // maximal distance in phi
4574 const Float_t kPtRatio = 0.3; // ratio between pt
4575 const Float_t kMinDCAR = 2.; // distance to the primary vertex in r - see cpipe cut
4578 // Curling tracks cuts
4581 const Float_t kMaxDeltaRMax = 40; // distance in outer radius
4582 const Float_t kMaxDeltaRMin = 5.; // distance in lower radius - see cpipe cut
4583 const Float_t kMinAngle = 2.9; // angle between tracks
4584 const Float_t kMaxDist = 5; // biggest distance
4586 // The cuts can be tuned using the "MC information stored in Multi tree ==> see FindMultiMC
4589 TCut csign("csign","Tr0.fP[4]*Tr1.fP[4]<0"); //opposite sign
4590 TCut cmax("cmax","abs(Tr0.GetC())>1/400");
4591 TCut cda("cda","sqrt(dtheta^2+dfi^2)<0.15");
4592 TCut ccratio("ccratio","abs((Tr0.fP[4]+Tr1.fP[4])/(abs(Tr0.fP[4])+abs(Tr1.fP[4])))<0.3");
4593 TCut cpipe("cpipe", "min(abs(r0-rc0),abs(r1-rc1))>5");
4595 TCut cdrmax("cdrmax","abs(abs(rc0+r0)-abs(rc1+r1))<40")
4596 TCut cdrmin("cdrmin","abs(abs(rc0+r0)-abs(rc1+r1))<10")
4598 Multi->Draw("dfi","iter==0"+csign+cmax+cda+ccratio); ~94% of curling tracks fulfill
4599 Multi->Draw("min(abs(r0-rc0),abs(r1-rc1))","iter==0&&abs(lab1)==abs(lab0)"+csign+cmax+cda+ccratio+cpipe+cdrmin+cdrmax); //80%
4601 Curling2->Draw("dfi","iter==0&&abs(lab0)==abs(lab1)"+csign+cmax+cdtheta+cdfi+ccratio)
4607 Int_t nentries = array->GetEntriesFast();
4608 AliHelix *helixes = new AliHelix[nentries];
4609 for (Int_t i=0;i<nentries;i++){
4610 AliTPCseed* track = (AliTPCseed*)array->At(i);
4611 if (!track) continue;
4612 track->SetCircular(0);
4613 new (&helixes[i]) AliHelix(*track);
4619 Double_t phase[2][2]={{0,0},{0,0}},radius[2]={0,0};
4625 for (Int_t i0=0;i0<nentries;i0++){
4626 AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
4627 if (!track0) continue;
4628 if (TMath::Abs(track0->GetC())<1/kMaxC) continue;
4629 Float_t xc0 = helixes[i0].GetHelix(6);
4630 Float_t yc0 = helixes[i0].GetHelix(7);
4631 Float_t r0 = helixes[i0].GetHelix(8);
4632 Float_t rc0 = TMath::Sqrt(xc0*xc0+yc0*yc0);
4633 Float_t fi0 = TMath::ATan2(yc0,xc0);
4635 for (Int_t i1=i0+1;i1<nentries;i1++){
4636 AliTPCseed * track1 = (AliTPCseed*)array->At(i1);
4637 if (!track1) continue;
4638 if (TMath::Abs(track1->GetC())<1/kMaxC) continue;
4639 Float_t xc1 = helixes[i1].GetHelix(6);
4640 Float_t yc1 = helixes[i1].GetHelix(7);
4641 Float_t r1 = helixes[i1].GetHelix(8);
4642 Float_t rc1 = TMath::Sqrt(xc1*xc1+yc1*yc1);
4643 Float_t fi1 = TMath::ATan2(yc1,xc1);
4645 Float_t dfi = fi0-fi1;
4648 if (dfi>1.5*TMath::Pi()) dfi-=TMath::Pi(); // take care about edge effect
4649 if (dfi<-1.5*TMath::Pi()) dfi+=TMath::Pi(); //
4650 Float_t dtheta = TMath::Abs(track0->GetTgl()-track1->GetTgl())<TMath::Abs(track0->GetTgl()+track1->GetTgl())? track0->GetTgl()-track1->GetTgl():track0->GetTgl()+track1->GetTgl();
4654 if (track0->GetBConstrain()&&track1->GetBConstrain()) continue; // not constrained
4655 if (track1->GetSigned1Pt()*track0->GetSigned1Pt()>0) continue; // not the same sign
4656 if ( TMath::Abs(track1->GetTgl()+track0->GetTgl())>kMaxdTheta) continue; //distance in the Theta
4657 if ( TMath::Abs(dfi)>kMaxdPhi) continue; //distance in phi
4658 if ( TMath::Sqrt(dfi*dfi+dtheta*dtheta)>kMaxdPhi) continue; //common angular offset
4660 Float_t pt0 = track0->GetSignedPt();
4661 Float_t pt1 = track1->GetSignedPt();
4662 if ((TMath::Abs(pt0+pt1)/(TMath::Abs(pt0)+TMath::Abs(pt1)))>kPtRatio) continue;
4663 if ((iter==1) && TMath::Abs(TMath::Abs(rc0+r0)-TMath::Abs(rc1+r1))>kMaxDeltaRMax) continue;
4664 if ((iter!=1) &&TMath::Abs(TMath::Abs(rc0-r0)-TMath::Abs(rc1-r1))>kMaxDeltaRMin) continue;
4665 if (TMath::Min(TMath::Abs(rc0-r0),TMath::Abs(rc1-r1))<kMinDCAR) continue;
4668 // Now find closest approach
4672 Int_t npoints = helixes[i0].GetRPHIintersections(helixes[i1], phase, radius,10);
4673 if (npoints==0) continue;
4674 helixes[i0].GetClosestPhases(helixes[i1], phase);
4678 Double_t hangles[3];
4679 helixes[i0].Evaluate(phase[0][0],xyz0);
4680 helixes[i1].Evaluate(phase[0][1],xyz1);
4682 helixes[i0].GetAngle(phase[0][0],helixes[i1],phase[0][1],hangles);
4683 Double_t deltah[2],deltabest;
4684 if (TMath::Abs(hangles[2])<kMinAngle) continue;
4688 helixes[i0].ParabolicDCA(helixes[i1],phase[0][0],phase[0][1],radius[0],deltah[0],2);
4690 helixes[i0].ParabolicDCA(helixes[i1],phase[1][0],phase[1][1],radius[1],deltah[1],2);
4691 if (deltah[1]<deltah[0]) ibest=1;
4693 deltabest = TMath::Sqrt(deltah[ibest]);
4694 helixes[i0].Evaluate(phase[ibest][0],xyz0);
4695 helixes[i1].Evaluate(phase[ibest][1],xyz1);
4696 helixes[i0].GetAngle(phase[ibest][0],helixes[i1],phase[ibest][1],hangles);
4697 Double_t radiusbest = TMath::Sqrt(radius[ibest]);
4699 if (deltabest>kMaxDist) continue;
4700 // if (mindcar+mindcaz<40 && (TMath::Abs(hangles[2])<kMinAngle ||deltabest>3)) continue;
4701 Bool_t sign =kFALSE;
4702 if (hangles[2]>kMinAngle) sign =kTRUE;
4705 // circular[i0] = kTRUE;
4706 // circular[i1] = kTRUE;
4707 if (track0->OneOverPt()<track1->OneOverPt()){
4708 track0->SetCircular(track0->GetCircular()+1);
4709 track1->SetCircular(track1->GetCircular()+2);
4712 track1->SetCircular(track1->GetCircular()+1);
4713 track0->SetCircular(track0->GetCircular()+2);
4716 if (AliTPCReconstructor::StreamLevel()>2){
4718 //debug stream to tune "fine" cuts
4719 Int_t lab0=track0->GetLabel();
4720 Int_t lab1=track1->GetLabel();
4721 TTreeSRedirector &cstream = *fDebugStreamer;
4722 cstream<<"Curling2"<<
4738 "npoints="<<npoints<<
4739 "hangles0="<<hangles[0]<<
4740 "hangles1="<<hangles[1]<<
4741 "hangles2="<<hangles[2]<<
4744 "radius="<<radiusbest<<
4745 "deltabest="<<deltabest<<
4746 "phase0="<<phase[ibest][0]<<
4747 "phase1="<<phase[ibest][1]<<
4755 if (AliTPCReconstructor::StreamLevel()>1) {
4756 AliInfo("Time for curling tracks removal");
4762 void AliTPCtrackerMI::FindKinks(TObjArray * array, AliESDEvent *esd)
4768 // RS something is wrong in this routine: not all seeds are assigned to daughters and mothers array, but they all are queried
4771 TObjArray *kinks= new TObjArray(10000);
4772 // TObjArray *v0s= new TObjArray(10000);
4773 Int_t nentries = array->GetEntriesFast();
4774 AliHelix *helixes = new AliHelix[nentries];
4775 Int_t *sign = new Int_t[nentries];
4776 Int_t *nclusters = new Int_t[nentries];
4777 Float_t *alpha = new Float_t[nentries];
4778 AliKink *kink = new AliKink();
4779 Int_t * usage = new Int_t[nentries];
4780 Float_t *zm = new Float_t[nentries];
4781 Float_t *z0 = new Float_t[nentries];
4782 Float_t *fim = new Float_t[nentries];
4783 Float_t *shared = new Float_t[nentries];
4784 Bool_t *circular = new Bool_t[nentries];
4785 Float_t *dca = new Float_t[nentries];
4786 //const AliESDVertex * primvertex = esd->GetVertex();
4788 // nentries = array->GetEntriesFast();
4793 for (Int_t i=0;i<nentries;i++){
4796 AliTPCseed* track = (AliTPCseed*)array->At(i);
4797 if (!track) continue;
4798 track->SetCircular(0);
4800 track->UpdatePoints();
4801 if (( track->GetPoints()[2]- track->GetPoints()[0])>5 && track->GetPoints()[3]>0.8){
4803 nclusters[i]=track->GetNumberOfClusters();
4804 alpha[i] = track->GetAlpha();
4805 new (&helixes[i]) AliHelix(*track);
4807 helixes[i].Evaluate(0,xyz);
4808 sign[i] = (track->GetC()>0) ? -1:1;
4811 if (track->GetProlongation(x,y,z)){
4813 fim[i] = alpha[i]+TMath::ATan2(y,x);
4816 zm[i] = track->GetZ();
4820 circular[i]= kFALSE;
4821 if (track->GetProlongation(0,y,z)) z0[i] = z;
4822 dca[i] = track->GetD(0,0);
4828 Int_t ncandidates =0;
4831 Double_t phase[2][2]={{0,0},{0,0}},radius[2]={0,0};
4834 // Find circling track
4836 for (Int_t i0=0;i0<nentries;i0++){
4837 AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
4838 if (!track0) continue;
4839 if (track0->GetNumberOfClusters()<40) continue;
4840 if (TMath::Abs(1./track0->GetC())>200) continue;
4841 for (Int_t i1=i0+1;i1<nentries;i1++){
4842 AliTPCseed * track1 = (AliTPCseed*)array->At(i1);
4843 if (!track1) continue;
4844 if (track1->GetNumberOfClusters()<40) continue;
4845 if ( TMath::Abs(track1->GetTgl()+track0->GetTgl())>0.1) continue;
4846 if (track0->GetBConstrain()&&track1->GetBConstrain()) continue;
4847 if (TMath::Abs(1./track1->GetC())>200) continue;
4848 if (track1->GetSigned1Pt()*track0->GetSigned1Pt()>0) continue;
4849 if (track1->GetTgl()*track0->GetTgl()>0) continue;
4850 if (TMath::Max(TMath::Abs(1./track0->GetC()),TMath::Abs(1./track1->GetC()))>190) continue;
4851 if (track0->GetBConstrain()&&track1->OneOverPt()<track0->OneOverPt()) continue; //returning - lower momenta
4852 if (track1->GetBConstrain()&&track0->OneOverPt()<track1->OneOverPt()) continue; //returning - lower momenta
4854 Float_t mindcar = TMath::Min(TMath::Abs(dca[i0]),TMath::Abs(dca[i1]));
4855 if (mindcar<5) continue;
4856 Float_t mindcaz = TMath::Min(TMath::Abs(z0[i0]-GetZ()),TMath::Abs(z0[i1]-GetZ()));
4857 if (mindcaz<5) continue;
4858 if (mindcar+mindcaz<20) continue;
4861 Float_t xc0 = helixes[i0].GetHelix(6);
4862 Float_t yc0 = helixes[i0].GetHelix(7);
4863 Float_t r0 = helixes[i0].GetHelix(8);
4864 Float_t xc1 = helixes[i1].GetHelix(6);
4865 Float_t yc1 = helixes[i1].GetHelix(7);
4866 Float_t r1 = helixes[i1].GetHelix(8);
4868 Float_t rmean = (r0+r1)*0.5;
4869 Float_t delta =TMath::Sqrt((xc1-xc0)*(xc1-xc0)+(yc1-yc0)*(yc1-yc0));
4870 //if (delta>30) continue;
4871 if (delta>rmean*0.25) continue;
4872 if (TMath::Abs(r0-r1)/rmean>0.3) continue;
4874 Int_t npoints = helixes[i0].GetRPHIintersections(helixes[i1], phase, radius,10);
4875 if (npoints==0) continue;
4876 helixes[i0].GetClosestPhases(helixes[i1], phase);
4880 Double_t hangles[3];
4881 helixes[i0].Evaluate(phase[0][0],xyz0);
4882 helixes[i1].Evaluate(phase[0][1],xyz1);
4884 helixes[i0].GetAngle(phase[0][0],helixes[i1],phase[0][1],hangles);
4885 Double_t deltah[2],deltabest;
4886 if (hangles[2]<2.8) continue;
4889 helixes[i0].ParabolicDCA(helixes[i1],phase[0][0],phase[0][1],radius[0],deltah[0],2);
4891 helixes[i0].ParabolicDCA(helixes[i1],phase[1][0],phase[1][1],radius[1],deltah[1],2);
4892 if (deltah[1]<deltah[0]) ibest=1;
4894 deltabest = TMath::Sqrt(deltah[ibest]);
4895 helixes[i0].Evaluate(phase[ibest][0],xyz0);
4896 helixes[i1].Evaluate(phase[ibest][1],xyz1);
4897 helixes[i0].GetAngle(phase[ibest][0],helixes[i1],phase[ibest][1],hangles);
4898 Double_t radiusbest = TMath::Sqrt(radius[ibest]);
4900 if (deltabest>6) continue;
4901 if (mindcar+mindcaz<40 && (hangles[2]<3.12||deltabest>3)) continue;
4902 Bool_t lsign =kFALSE;
4903 if (hangles[2]>3.06) lsign =kTRUE;
4906 circular[i0] = kTRUE;
4907 circular[i1] = kTRUE;
4908 if (track0->OneOverPt()<track1->OneOverPt()){
4909 track0->SetCircular(track0->GetCircular()+1);
4910 track1->SetCircular(track1->GetCircular()+2);
4913 track1->SetCircular(track1->GetCircular()+1);
4914 track0->SetCircular(track0->GetCircular()+2);
4917 if (lsign&&AliTPCReconstructor::StreamLevel()>1){
4919 Int_t lab0=track0->GetLabel();
4920 Int_t lab1=track1->GetLabel();
4921 TTreeSRedirector &cstream = *fDebugStreamer;
4922 cstream<<"Curling"<<
4929 "mindcar="<<mindcar<<
4930 "mindcaz="<<mindcaz<<
4933 "npoints="<<npoints<<
4934 "hangles0="<<hangles[0]<<
4935 "hangles2="<<hangles[2]<<
4940 "radius="<<radiusbest<<
4941 "deltabest="<<deltabest<<
4942 "phase0="<<phase[ibest][0]<<
4943 "phase1="<<phase[ibest][1]<<
4953 for (Int_t i =0;i<nentries;i++){
4954 if (sign[i]==0) continue;
4955 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
4962 Double_t cradius0 = 40*40;
4963 Double_t cradius1 = 270*270;
4966 Double_t cdist3=0.55;
4967 for (Int_t j =i+1;j<nentries;j++){
4969 if (sign[j]*sign[i]<1) continue;
4970 if ( (nclusters[i]+nclusters[j])>200) continue;
4971 if ( (nclusters[i]+nclusters[j])<80) continue;
4972 if ( TMath::Abs(zm[i]-zm[j])>60.) continue;
4973 if ( TMath::Abs(fim[i]-fim[j])>0.6 && TMath::Abs(fim[i]-fim[j])<5.7 ) continue;
4974 //AliTPCseed * track1 = (AliTPCseed*)array->At(j); Double_t phase[2][2],radius[2];
4975 Int_t npoints = helixes[i].GetRPHIintersections(helixes[j], phase, radius,20);
4976 if (npoints<1) continue;
4979 if (radius[0]<cradius0||radius[0]>cradius1) continue;
4982 if ( (radius[0]<cradius0||radius[0]>cradius1) && (radius[1]<cradius0||radius[1]>cradius1) ) continue;
4985 Double_t delta1=10000,delta2=10000;
4986 // cuts on the intersection radius
4987 helixes[i].LinearDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta1);
4988 if (radius[0]<20&&delta1<1) continue; //intersection at vertex
4989 if (radius[0]<10&&delta1<3) continue; //intersection at vertex
4991 helixes[i].LinearDCA(helixes[j],phase[1][0],phase[1][1],radius[1],delta2);
4992 if (radius[1]<20&&delta2<1) continue; //intersection at vertex
4993 if (radius[1]<10&&delta2<3) continue; //intersection at vertex
4996 Double_t distance1 = TMath::Min(delta1,delta2);
4997 if (distance1>cdist1) continue; // cut on DCA linear approximation
4999 npoints = helixes[i].GetRPHIintersections(helixes[j], phase, radius,20);
5000 helixes[i].ParabolicDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta1);
5001 if (radius[0]<20&&delta1<1) continue; //intersection at vertex
5002 if (radius[0]<10&&delta1<3) continue; //intersection at vertex
5005 helixes[i].ParabolicDCA(helixes[j],phase[1][0],phase[1][1],radius[1],delta2);
5006 if (radius[1]<20&&delta2<1) continue; //intersection at vertex
5007 if (radius[1]<10&&delta2<3) continue; //intersection at vertex
5009 distance1 = TMath::Min(delta1,delta2);
5012 rkink = TMath::Sqrt(radius[0]);
5015 rkink = TMath::Sqrt(radius[1]);
5017 if (distance1>cdist2) continue;
5020 AliTPCseed * track1 = (AliTPCseed*)array->At(j);
5023 Int_t row0 = GetRowNumber(rkink);
5024 if (row0<10) continue;
5025 if (row0>150) continue;
5028 Float_t dens00=-1,dens01=-1;
5029 Float_t dens10=-1,dens11=-1;
5031 Int_t found,foundable,ishared;
5032 track0->GetClusterStatistic(0,row0-5, found, foundable,ishared,kFALSE);
5033 if (foundable>5) dens00 = Float_t(found)/Float_t(foundable);
5034 track0->GetClusterStatistic(row0+5,155, found, foundable,ishared,kFALSE);
5035 if (foundable>5) dens01 = Float_t(found)/Float_t(foundable);
5037 track1->GetClusterStatistic(0,row0-5, found, foundable,ishared,kFALSE);
5038 if (foundable>10) dens10 = Float_t(found)/Float_t(foundable);
5039 track1->GetClusterStatistic(row0+5,155, found, foundable,ishared,kFALSE);
5040 if (foundable>10) dens11 = Float_t(found)/Float_t(foundable);
5042 if (dens00<dens10 && dens01<dens11) continue;
5043 if (dens00>dens10 && dens01>dens11) continue;
5044 if (TMath::Max(dens00,dens10)<0.1) continue;
5045 if (TMath::Max(dens01,dens11)<0.3) continue;
5047 if (TMath::Min(dens00,dens10)>0.6) continue;
5048 if (TMath::Min(dens01,dens11)>0.6) continue;
5051 AliTPCseed * ktrack0, *ktrack1;
5060 if (TMath::Abs(ktrack0->GetC())>5) continue; // cut on the curvature for mother particle
5061 AliExternalTrackParam paramm(*ktrack0);
5062 AliExternalTrackParam paramd(*ktrack1);
5063 if (row0>60&&ktrack1->GetReference().GetX()>90.)new (¶md) AliExternalTrackParam(ktrack1->GetReference());
5066 kink->SetMother(paramm);
5067 kink->SetDaughter(paramd);
5070 Float_t x[3] = { kink->GetPosition()[0],kink->GetPosition()[1],kink->GetPosition()[2]};
5072 fkParam->Transform0to1(x,index);
5073 fkParam->Transform1to2(x,index);
5074 row0 = GetRowNumber(x[0]);
5076 if (kink->GetR()<100) continue;
5077 if (kink->GetR()>240) continue;
5078 if (kink->GetPosition()[2]/kink->GetR()>AliTPCReconstructor::GetCtgRange()) continue; //out of fiducial volume
5079 if (kink->GetDistance()>cdist3) continue;
5080 Float_t dird = kink->GetDaughterP()[0]*kink->GetPosition()[0]+kink->GetDaughterP()[1]*kink->GetPosition()[1]; // rough direction estimate
5081 if (dird<0) continue;
5083 Float_t dirm = kink->GetMotherP()[0]*kink->GetPosition()[0]+kink->GetMotherP()[1]*kink->GetPosition()[1]; // rough direction estimate
5084 if (dirm<0) continue;
5085 Float_t mpt = TMath::Sqrt(kink->GetMotherP()[0]*kink->GetMotherP()[0]+kink->GetMotherP()[1]*kink->GetMotherP()[1]);
5086 if (mpt<0.2) continue;
5089 //for high momenta momentum not defined well in first iteration
5090 Double_t qt = TMath::Sin(kink->GetAngle(2))*ktrack1->GetP();
5091 if (qt>0.35) continue;
5094 kink->SetLabel(CookLabel(ktrack0,0.4,0,row0),0);
5095 kink->SetLabel(CookLabel(ktrack1,0.4,row0,160),1);
5097 kink->SetTPCDensity(dens00,0,0);
5098 kink->SetTPCDensity(dens01,0,1);
5099 kink->SetTPCDensity(dens10,1,0);
5100 kink->SetTPCDensity(dens11,1,1);
5101 kink->SetIndex(i,0);
5102 kink->SetIndex(j,1);
5105 kink->SetTPCDensity(dens10,0,0);
5106 kink->SetTPCDensity(dens11,0,1);
5107 kink->SetTPCDensity(dens00,1,0);
5108 kink->SetTPCDensity(dens01,1,1);
5109 kink->SetIndex(j,0);
5110 kink->SetIndex(i,1);
5113 if (mpt<1||kink->GetAngle(2)>0.1){
5114 // angle and densities not defined yet
5115 if (kink->GetTPCDensityFactor()<0.8) continue;
5116 if ((2-kink->GetTPCDensityFactor())*kink->GetDistance() >0.25) continue;
5117 if (kink->GetAngle(2)*ktrack0->GetP()<0.003) continue; //too small angle
5118 if (kink->GetAngle(2)>0.2&&kink->GetTPCDensityFactor()<1.15) continue;
5119 if (kink->GetAngle(2)>0.2&&kink->GetTPCDensity(0,1)>0.05) continue;
5121 Float_t criticalangle = track0->GetSigmaSnp2()+track0->GetSigmaTgl2();
5122 criticalangle+= track1->GetSigmaSnp2()+track1->GetSigmaTgl2();
5123 criticalangle= 3*TMath::Sqrt(criticalangle);
5124 if (criticalangle>0.02) criticalangle=0.02;
5125 if (kink->GetAngle(2)<criticalangle) continue;
5128 Int_t drow = Int_t(2.+0.5/(0.05+kink->GetAngle(2))); // overlap region defined
5129 Float_t shapesum =0;
5131 for ( Int_t row = row0-drow; row<row0+drow;row++){
5132 if (row<0) continue;
5133 if (row>155) continue;
5134 if (ktrack0->GetClusterPointer(row)){
5135 AliTPCTrackerPoint *point =ktrack0->GetTrackPoint(row);
5136 shapesum+=point->GetSigmaY()+point->GetSigmaZ();
5139 if (ktrack1->GetClusterPointer(row)){
5140 AliTPCTrackerPoint *point =ktrack1->GetTrackPoint(row);
5141 shapesum+=point->GetSigmaY()+point->GetSigmaZ();
5146 kink->SetShapeFactor(-1.);
5149 kink->SetShapeFactor(shapesum/sum);
5151 // esd->AddKink(kink);
5153 // kink->SetMother(paramm);
5154 //kink->SetDaughter(paramd);
5156 Double_t chi2P2 = paramm.GetParameter()[2]-paramd.GetParameter()[2];
5158 chi2P2/=paramm.GetCovariance()[5]+paramd.GetCovariance()[5];
5159 Double_t chi2P3 = paramm.GetParameter()[3]-paramd.GetParameter()[3];
5161 chi2P3/=paramm.GetCovariance()[9]+paramd.GetCovariance()[9];
5163 if (AliTPCReconstructor::StreamLevel()>1) {
5164 (*fDebugStreamer)<<"kinkLpt"<<
5172 if ( chi2P2+chi2P3<AliTPCReconstructor::GetRecoParam()->GetKinkAngleCutChi2(0)){
5176 kinks->AddLast(kink);
5182 // sort the kinks according quality - and refit them towards vertex
5184 Int_t nkinks = kinks->GetEntriesFast();
5185 Float_t *quality = new Float_t[nkinks];
5186 Int_t *indexes = new Int_t[nkinks];
5187 AliTPCseed *mothers = new AliTPCseed[nkinks];
5188 AliTPCseed *daughters = new AliTPCseed[nkinks];
5191 for (Int_t i=0;i<nkinks;i++){
5193 AliKink *kinkl = (AliKink*)kinks->At(i);
5195 // refit kinks towards vertex
5197 Int_t index0 = kinkl->GetIndex(0);
5198 Int_t index1 = kinkl->GetIndex(1);
5199 AliTPCseed * ktrack0 = (AliTPCseed*)array->At(index0);
5200 AliTPCseed * ktrack1 = (AliTPCseed*)array->At(index1);
5202 Int_t sumn=ktrack0->GetNumberOfClusters()+ktrack1->GetNumberOfClusters();
5204 // Refit Kink under if too small angle
5206 if (kinkl->GetAngle(2)<0.05){
5207 kinkl->SetTPCRow0(GetRowNumber(kinkl->GetR()));
5208 Int_t row0 = kinkl->GetTPCRow0();
5209 Int_t drow = Int_t(2.+0.5/(0.05+kinkl->GetAngle(2)));
5212 Int_t last = row0-drow;
5213 if (last<40) last=40;
5214 if (last<ktrack0->GetFirstPoint()+25) last = ktrack0->GetFirstPoint()+25;
5215 AliTPCseed* seed0 = ReSeed(ktrack0,last,kFALSE);
5218 Int_t first = row0+drow;
5219 if (first>130) first=130;
5220 if (first>ktrack1->GetLastPoint()-25) first = TMath::Max(ktrack1->GetLastPoint()-25,30);
5221 AliTPCseed* seed1 = ReSeed(ktrack1,first,kTRUE);
5223 if (seed0 && seed1){
5224 kinkl->SetStatus(1,8);
5225 if (RefitKink(*seed0,*seed1,*kinkl)) kinkl->SetStatus(1,9);
5226 row0 = GetRowNumber(kinkl->GetR());
5227 sumn = seed0->GetNumberOfClusters()+seed1->GetNumberOfClusters();
5228 mothers[i] = *seed0;
5229 daughters[i] = *seed1;
5232 delete kinks->RemoveAt(i);
5233 if (seed0) MarkSeedFree( seed0 );
5234 if (seed1) MarkSeedFree( seed1 );
5237 if (kinkl->GetDistance()>0.5 || kinkl->GetR()<110 || kinkl->GetR()>240) {
5238 delete kinks->RemoveAt(i);
5239 if (seed0) MarkSeedFree( seed0 );
5240 if (seed1) MarkSeedFree( seed1 );
5244 MarkSeedFree( seed0 );
5245 MarkSeedFree( seed1 );
5248 if (kinkl) quality[i] = 160*((0.1+kinkl->GetDistance())*(2.-kinkl->GetTPCDensityFactor()))/(sumn+40.); //the longest -clossest will win
5250 TMath::Sort(nkinks,quality,indexes,kFALSE);
5252 //remove double find kinks
5254 for (Int_t ikink0=1;ikink0<nkinks;ikink0++){
5255 AliKink * kink0 = (AliKink*) kinks->At(indexes[ikink0]);
5256 if (!kink0) continue;
5258 for (Int_t ikink1=0;ikink1<ikink0;ikink1++){
5259 kink0 = (AliKink*) kinks->At(indexes[ikink0]);
5260 if (!kink0) continue;
5261 AliKink * kink1 = (AliKink*) kinks->At(indexes[ikink1]);
5262 if (!kink1) continue;
5263 // if not close kink continue
5264 if (TMath::Abs(kink1->GetPosition()[2]-kink0->GetPosition()[2])>10) continue;
5265 if (TMath::Abs(kink1->GetPosition()[1]-kink0->GetPosition()[1])>10) continue;
5266 if (TMath::Abs(kink1->GetPosition()[0]-kink0->GetPosition()[0])>10) continue;
5268 AliTPCseed &mother0 = mothers[indexes[ikink0]];
5269 AliTPCseed &daughter0 = daughters[indexes[ikink0]];
5270 AliTPCseed &mother1 = mothers[indexes[ikink1]];
5271 AliTPCseed &daughter1 = daughters[indexes[ikink1]];
5272 Int_t row0 = (kink0->GetTPCRow0()+kink1->GetTPCRow0())/2;
5281 for (Int_t i=0;i<row0;i++){
5282 if (mother0.GetClusterIndex(i)>0 && mother1.GetClusterIndex(i)>0){
5285 if (mother0.GetClusterIndex(i)==mother1.GetClusterIndex(i)){
5292 for (Int_t i=row0;i<158;i++){
5293 //if (daughter0.GetClusterIndex(i)>0 && daughter0.GetClusterIndex(i)>0){ // RS: Bug?
5294 if (daughter0.GetClusterIndex(i)>0 && daughter1.GetClusterIndex(i)>0){
5297 if (mother0.GetClusterIndex(i)==mother1.GetClusterIndex(i)){
5303 Float_t ratio = Float_t(same+1)/Float_t(both+1);
5304 Float_t ratiom = Float_t(samem+1)/Float_t(bothm+1);
5305 Float_t ratiod = Float_t(samed+1)/Float_t(bothd+1);
5306 if (ratio>0.3 && ratiom>0.5 &&ratiod>0.5) {
5307 Int_t sum0 = mother0.GetNumberOfClusters()+daughter0.GetNumberOfClusters();
5308 Int_t sum1 = mother1.GetNumberOfClusters()+daughter1.GetNumberOfClusters();
5310 shared[kink0->GetIndex(0)]= kTRUE;
5311 shared[kink0->GetIndex(1)]= kTRUE;
5312 delete kinks->RemoveAt(indexes[ikink0]);
5316 shared[kink1->GetIndex(0)]= kTRUE;
5317 shared[kink1->GetIndex(1)]= kTRUE;
5318 delete kinks->RemoveAt(indexes[ikink1]);
5325 for (Int_t i=0;i<nkinks;i++){
5326 AliKink * kinkl = (AliKink*) kinks->At(indexes[i]);
5327 if (!kinkl) continue;
5328 kinkl->SetTPCRow0(GetRowNumber(kinkl->GetR()));
5329 Int_t index0 = kinkl->GetIndex(0);
5330 Int_t index1 = kinkl->GetIndex(1);
5331 if (circular[index0]||(circular[index1]&&kinkl->GetDistance()>0.2)) continue;
5332 kinkl->SetMultiple(usage[index0],0);
5333 kinkl->SetMultiple(usage[index1],1);
5334 if (kinkl->GetMultiple()[0]+kinkl->GetMultiple()[1]>2) continue;
5335 if (kinkl->GetMultiple()[0]+kinkl->GetMultiple()[1]>0 && quality[indexes[i]]>0.2) continue;
5336 if (kinkl->GetMultiple()[0]+kinkl->GetMultiple()[1]>0 && kinkl->GetDistance()>0.2) continue;
5337 if (circular[index0]||(circular[index1]&&kinkl->GetDistance()>0.1)) continue;
5339 AliTPCseed * ktrack0 = (AliTPCseed*)array->At(index0);
5340 AliTPCseed * ktrack1 = (AliTPCseed*)array->At(index1);
5341 if (!ktrack0 || !ktrack1) continue;
5342 Int_t index = esd->AddKink(kinkl);
5345 if ( ktrack0->GetKinkIndex(0)==0 && ktrack1->GetKinkIndex(0)==0) { //best kink
5346 if (mothers[indexes[i]].GetNumberOfClusters()>20 && daughters[indexes[i]].GetNumberOfClusters()>20 && (mothers[indexes[i]].GetNumberOfClusters()+daughters[indexes[i]].GetNumberOfClusters())>100){
5347 *ktrack0 = mothers[indexes[i]];
5348 *ktrack1 = daughters[indexes[i]];
5352 ktrack0->SetKinkIndex(usage[index0],-(index+1));
5353 ktrack1->SetKinkIndex(usage[index1], (index+1));
5358 // Remove tracks corresponding to shared kink's
5360 for (Int_t i=0;i<nentries;i++){
5361 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5362 if (!track0) continue;
5363 if (track0->GetKinkIndex(0)!=0) continue;
5364 if (shared[i]) MarkSeedFree( array->RemoveAt(i) );
5369 RemoveUsed2(array,0.5,0.4,30);
5371 for (Int_t i=0;i<nentries;i++){
5372 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5373 if (!track0) continue;
5374 track0->CookdEdx(0.02,0.6);
5378 for (Int_t i=0;i<nentries;i++){
5379 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5380 if (!track0) continue;
5381 if (track0->Pt()<1.4) continue;
5382 //remove double high momenta tracks - overlapped with kink candidates
5385 for (Int_t icl=track0->GetFirstPoint();icl<track0->GetLastPoint(); icl++){
5386 if (track0->GetClusterPointer(icl)!=0){
5388 if (track0->GetClusterPointer(icl)->IsUsed(10)) ishared++;
5391 if (Float_t(ishared+1)/Float_t(all+1)>0.5) {
5392 MarkSeedFree( array->RemoveAt(i) );
5396 if (track0->GetKinkIndex(0)!=0) continue;
5397 if (track0->GetNumberOfClusters()<80) continue;
5399 AliTPCseed *pmother = new AliTPCseed();
5400 AliTPCseed *pdaughter = new AliTPCseed();
5401 AliKink *pkink = new AliKink;
5403 AliTPCseed & mother = *pmother;
5404 AliTPCseed & daughter = *pdaughter;
5405 AliKink & kinkl = *pkink;
5406 if (CheckKinkPoint(track0,mother,daughter, kinkl)){
5407 if (mother.GetNumberOfClusters()<30||daughter.GetNumberOfClusters()<20) {
5411 continue; //too short tracks
5413 if (mother.Pt()<1.4) {
5419 Int_t row0= kinkl.GetTPCRow0();
5420 if (kinkl.GetDistance()>0.5 || kinkl.GetR()<110. || kinkl.GetR()>240.) {
5427 Int_t index = esd->AddKink(&kinkl);
5428 mother.SetKinkIndex(0,-(index+1));
5429 daughter.SetKinkIndex(0,index+1);
5430 if (mother.GetNumberOfClusters()>50) {
5431 MarkSeedFree( array->RemoveAt(i) );
5432 AliTPCseed* mtc = new( NextFreeSeed() ) AliTPCseed(mother);
5433 mtc->SetPoolID(fLastSeedID);
5434 array->AddAt(mtc,i);
5437 AliTPCseed* mtc = new( NextFreeSeed() ) AliTPCseed(mother);
5438 mtc->SetPoolID(fLastSeedID);
5439 array->AddLast(mtc);
5441 AliTPCseed* dtc = new( NextFreeSeed() ) AliTPCseed(daughter);
5442 dtc->SetPoolID(fLastSeedID);
5443 array->AddLast(dtc);
5444 for (Int_t icl=0;icl<row0;icl++) {
5445 if (mother.GetClusterPointer(icl)) mother.GetClusterPointer(icl)->Use(20);
5448 for (Int_t icl=row0;icl<158;icl++) {
5449 if (daughter.GetClusterPointer(icl)) daughter.GetClusterPointer(icl)->Use(20);
5458 delete [] daughters;
5480 AliInfo(Form("Ncandidates=\t%d\t%d\t%d\t%d\n",esd->GetNumberOfKinks(),ncandidates,ntracks,nall));
5486 void AliTPCtrackerMI::FindKinks(TObjArray * array, AliESDEvent *esd)
5493 TObjArray *kinks= new TObjArray(10000);
5494 // TObjArray *v0s= new TObjArray(10000);
5495 Int_t nentries = array->GetEntriesFast();
5496 AliHelix *helixes = new AliHelix[nentries];
5497 Int_t *sign = new Int_t[nentries];
5498 Int_t *nclusters = new Int_t[nentries];
5499 Float_t *alpha = new Float_t[nentries];
5500 AliKink *kink = new AliKink();
5501 Int_t * usage = new Int_t[nentries];
5502 Float_t *zm = new Float_t[nentries];
5503 Float_t *z0 = new Float_t[nentries];
5504 Float_t *fim = new Float_t[nentries];
5505 Float_t *shared = new Float_t[nentries];
5506 Bool_t *circular = new Bool_t[nentries];
5507 Float_t *dca = new Float_t[nentries];
5508 //const AliESDVertex * primvertex = esd->GetVertex();
5510 // nentries = array->GetEntriesFast();
5515 for (Int_t i=0;i<nentries;i++){
5518 AliTPCseed* track = (AliTPCseed*)array->At(i);
5519 if (!track) continue;
5520 track->SetCircular(0);
5522 track->UpdatePoints();
5523 if (( track->GetPoints()[2]- track->GetPoints()[0])>5 && track->GetPoints()[3]>0.8){
5525 nclusters[i]=track->GetNumberOfClusters();
5526 alpha[i] = track->GetAlpha();
5527 new (&helixes[i]) AliHelix(*track);
5529 helixes[i].Evaluate(0,xyz);
5530 sign[i] = (track->GetC()>0) ? -1:1;
5533 if (track->GetProlongation(x,y,z)){
5535 fim[i] = alpha[i]+TMath::ATan2(y,x);
5538 zm[i] = track->GetZ();
5542 circular[i]= kFALSE;
5543 if (track->GetProlongation(0,y,z)) z0[i] = z;
5544 dca[i] = track->GetD(0,0);
5550 Int_t ncandidates =0;
5553 Double_t phase[2][2]={{0,0},{0,0}},radius[2]={0,0};
5556 // Find circling track
5558 for (Int_t i0=0;i0<nentries;i0++){
5559 AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
5560 if (!track0) continue;
5561 if (track0->GetNumberOfClusters()<40) continue;
5562 if (TMath::Abs(1./track0->GetC())>200) continue;
5563 for (Int_t i1=i0+1;i1<nentries;i1++){
5564 AliTPCseed * track1 = (AliTPCseed*)array->At(i1);
5565 if (!track1) continue;
5566 if (track1->GetNumberOfClusters()<40) continue;
5567 if ( TMath::Abs(track1->GetTgl()+track0->GetTgl())>0.1) continue;
5568 if (track0->GetBConstrain()&&track1->GetBConstrain()) continue;
5569 if (TMath::Abs(1./track1->GetC())>200) continue;
5570 if (track1->GetSigned1Pt()*track0->GetSigned1Pt()>0) continue;
5571 if (track1->GetTgl()*track0->GetTgl()>0) continue;
5572 if (TMath::Max(TMath::Abs(1./track0->GetC()),TMath::Abs(1./track1->GetC()))>190) continue;
5573 if (track0->GetBConstrain()&&track1->OneOverPt()<track0->OneOverPt()) continue; //returning - lower momenta
5574 if (track1->GetBConstrain()&&track0->OneOverPt()<track1->OneOverPt()) continue; //returning - lower momenta
5576 Float_t mindcar = TMath::Min(TMath::Abs(dca[i0]),TMath::Abs(dca[i1]));
5577 if (mindcar<5) continue;
5578 Float_t mindcaz = TMath::Min(TMath::Abs(z0[i0]-GetZ()),TMath::Abs(z0[i1]-GetZ()));
5579 if (mindcaz<5) continue;
5580 if (mindcar+mindcaz<20) continue;
5583 Float_t xc0 = helixes[i0].GetHelix(6);
5584 Float_t yc0 = helixes[i0].GetHelix(7);
5585 Float_t r0 = helixes[i0].GetHelix(8);
5586 Float_t xc1 = helixes[i1].GetHelix(6);
5587 Float_t yc1 = helixes[i1].GetHelix(7);
5588 Float_t r1 = helixes[i1].GetHelix(8);
5590 Float_t rmean = (r0+r1)*0.5;
5591 Float_t delta =TMath::Sqrt((xc1-xc0)*(xc1-xc0)+(yc1-yc0)*(yc1-yc0));
5592 //if (delta>30) continue;
5593 if (delta>rmean*0.25) continue;
5594 if (TMath::Abs(r0-r1)/rmean>0.3) continue;
5596 Int_t npoints = helixes[i0].GetRPHIintersections(helixes[i1], phase, radius,10);
5597 if (npoints==0) continue;
5598 helixes[i0].GetClosestPhases(helixes[i1], phase);
5602 Double_t hangles[3];
5603 helixes[i0].Evaluate(phase[0][0],xyz0);
5604 helixes[i1].Evaluate(phase[0][1],xyz1);
5606 helixes[i0].GetAngle(phase[0][0],helixes[i1],phase[0][1],hangles);
5607 Double_t deltah[2],deltabest;
5608 if (hangles[2]<2.8) continue;
5611 helixes[i0].ParabolicDCA(helixes[i1],phase[0][0],phase[0][1],radius[0],deltah[0],2);
5613 helixes[i0].ParabolicDCA(helixes[i1],phase[1][0],phase[1][1],radius[1],deltah[1],2);
5614 if (deltah[1]<deltah[0]) ibest=1;
5616 deltabest = TMath::Sqrt(deltah[ibest]);
5617 helixes[i0].Evaluate(phase[ibest][0],xyz0);
5618 helixes[i1].Evaluate(phase[ibest][1],xyz1);
5619 helixes[i0].GetAngle(phase[ibest][0],helixes[i1],phase[ibest][1],hangles);
5620 Double_t radiusbest = TMath::Sqrt(radius[ibest]);
5622 if (deltabest>6) continue;
5623 if (mindcar+mindcaz<40 && (hangles[2]<3.12||deltabest>3)) continue;
5624 Bool_t lsign =kFALSE;
5625 if (hangles[2]>3.06) lsign =kTRUE;
5628 circular[i0] = kTRUE;
5629 circular[i1] = kTRUE;
5630 if (track0->OneOverPt()<track1->OneOverPt()){
5631 track0->SetCircular(track0->GetCircular()+1);
5632 track1->SetCircular(track1->GetCircular()+2);
5635 track1->SetCircular(track1->GetCircular()+1);
5636 track0->SetCircular(track0->GetCircular()+2);
5639 if (lsign&&AliTPCReconstructor::StreamLevel()>1){
5641 Int_t lab0=track0->GetLabel();
5642 Int_t lab1=track1->GetLabel();
5643 TTreeSRedirector &cstream = *fDebugStreamer;
5644 cstream<<"Curling"<<
5651 "mindcar="<<mindcar<<
5652 "mindcaz="<<mindcaz<<
5655 "npoints="<<npoints<<
5656 "hangles0="<<hangles[0]<<
5657 "hangles2="<<hangles[2]<<
5662 "radius="<<radiusbest<<
5663 "deltabest="<<deltabest<<
5664 "phase0="<<phase[ibest][0]<<
5665 "phase1="<<phase[ibest][1]<<
5675 for (Int_t i =0;i<nentries;i++){
5676 if (sign[i]==0) continue;
5677 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5684 Double_t cradius0 = 40*40;
5685 Double_t cradius1 = 270*270;
5688 Double_t cdist3=0.55;
5689 for (Int_t j =i+1;j<nentries;j++){
5691 if (sign[j]*sign[i]<1) continue;
5692 if ( (nclusters[i]+nclusters[j])>200) continue;
5693 if ( (nclusters[i]+nclusters[j])<80) continue;
5694 if ( TMath::Abs(zm[i]-zm[j])>60.) continue;
5695 if ( TMath::Abs(fim[i]-fim[j])>0.6 && TMath::Abs(fim[i]-fim[j])<5.7 ) continue;
5696 //AliTPCseed * track1 = (AliTPCseed*)array->At(j); Double_t phase[2][2],radius[2];
5697 Int_t npoints = helixes[i].GetRPHIintersections(helixes[j], phase, radius,20);
5698 if (npoints<1) continue;
5701 if (radius[0]<cradius0||radius[0]>cradius1) continue;
5704 if ( (radius[0]<cradius0||radius[0]>cradius1) && (radius[1]<cradius0||radius[1]>cradius1) ) continue;
5707 Double_t delta1=10000,delta2=10000;
5708 // cuts on the intersection radius
5709 helixes[i].LinearDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta1);
5710 if (radius[0]<20&&delta1<1) continue; //intersection at vertex
5711 if (radius[0]<10&&delta1<3) continue; //intersection at vertex
5713 helixes[i].LinearDCA(helixes[j],phase[1][0],phase[1][1],radius[1],delta2);
5714 if (radius[1]<20&&delta2<1) continue; //intersection at vertex
5715 if (radius[1]<10&&delta2<3) continue; //intersection at vertex
5718 Double_t distance1 = TMath::Min(delta1,delta2);
5719 if (distance1>cdist1) continue; // cut on DCA linear approximation
5721 npoints = helixes[i].GetRPHIintersections(helixes[j], phase, radius,20);
5722 helixes[i].ParabolicDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta1);
5723 if (radius[0]<20&&delta1<1) continue; //intersection at vertex
5724 if (radius[0]<10&&delta1<3) continue; //intersection at vertex
5727 helixes[i].ParabolicDCA(helixes[j],phase[1][0],phase[1][1],radius[1],delta2);
5728 if (radius[1]<20&&delta2<1) continue; //intersection at vertex
5729 if (radius[1]<10&&delta2<3) continue; //intersection at vertex
5731 distance1 = TMath::Min(delta1,delta2);
5734 rkink = TMath::Sqrt(radius[0]);
5737 rkink = TMath::Sqrt(radius[1]);
5739 if (distance1>cdist2) continue;
5742 AliTPCseed * track1 = (AliTPCseed*)array->At(j);
5745 Int_t row0 = GetRowNumber(rkink);
5746 if (row0<10) continue;
5747 if (row0>150) continue;
5750 Float_t dens00=-1,dens01=-1;
5751 Float_t dens10=-1,dens11=-1;
5753 Int_t found,foundable,ishared;
5754 track0->GetClusterStatistic(0,row0-5, found, foundable,ishared,kFALSE);
5755 if (foundable>5) dens00 = Float_t(found)/Float_t(foundable);
5756 track0->GetClusterStatistic(row0+5,155, found, foundable,ishared,kFALSE);
5757 if (foundable>5) dens01 = Float_t(found)/Float_t(foundable);
5759 track1->GetClusterStatistic(0,row0-5, found, foundable,ishared,kFALSE);
5760 if (foundable>10) dens10 = Float_t(found)/Float_t(foundable);
5761 track1->GetClusterStatistic(row0+5,155, found, foundable,ishared,kFALSE);
5762 if (foundable>10) dens11 = Float_t(found)/Float_t(foundable);
5764 if (dens00<dens10 && dens01<dens11) continue;
5765 if (dens00>dens10 && dens01>dens11) continue;
5766 if (TMath::Max(dens00,dens10)<0.1) continue;
5767 if (TMath::Max(dens01,dens11)<0.3) continue;
5769 if (TMath::Min(dens00,dens10)>0.6) continue;
5770 if (TMath::Min(dens01,dens11)>0.6) continue;
5773 AliTPCseed * ktrack0, *ktrack1;
5782 if (TMath::Abs(ktrack0->GetC())>5) continue; // cut on the curvature for mother particle
5783 AliExternalTrackParam paramm(*ktrack0);
5784 AliExternalTrackParam paramd(*ktrack1);
5785 if (row0>60&&ktrack1->GetReference().GetX()>90.)new (¶md) AliExternalTrackParam(ktrack1->GetReference());
5788 kink->SetMother(paramm);
5789 kink->SetDaughter(paramd);
5792 Float_t x[3] = { kink->GetPosition()[0],kink->GetPosition()[1],kink->GetPosition()[2]};
5794 fkParam->Transform0to1(x,index);
5795 fkParam->Transform1to2(x,index);
5796 row0 = GetRowNumber(x[0]);
5798 if (kink->GetR()<100) continue;
5799 if (kink->GetR()>240) continue;
5800 if (kink->GetPosition()[2]/kink->GetR()>AliTPCReconstructor::GetCtgRange()) continue; //out of fiducial volume
5801 if (kink->GetDistance()>cdist3) continue;
5802 Float_t dird = kink->GetDaughterP()[0]*kink->GetPosition()[0]+kink->GetDaughterP()[1]*kink->GetPosition()[1]; // rough direction estimate
5803 if (dird<0) continue;
5805 Float_t dirm = kink->GetMotherP()[0]*kink->GetPosition()[0]+kink->GetMotherP()[1]*kink->GetPosition()[1]; // rough direction estimate
5806 if (dirm<0) continue;
5807 Float_t mpt = TMath::Sqrt(kink->GetMotherP()[0]*kink->GetMotherP()[0]+kink->GetMotherP()[1]*kink->GetMotherP()[1]);
5808 if (mpt<0.2) continue;
5811 //for high momenta momentum not defined well in first iteration
5812 Double_t qt = TMath::Sin(kink->GetAngle(2))*ktrack1->GetP();
5813 if (qt>0.35) continue;
5816 kink->SetLabel(CookLabel(ktrack0,0.4,0,row0),0);
5817 kink->SetLabel(CookLabel(ktrack1,0.4,row0,160),1);
5819 kink->SetTPCDensity(dens00,0,0);
5820 kink->SetTPCDensity(dens01,0,1);
5821 kink->SetTPCDensity(dens10,1,0);
5822 kink->SetTPCDensity(dens11,1,1);
5823 kink->SetIndex(i,0);
5824 kink->SetIndex(j,1);
5827 kink->SetTPCDensity(dens10,0,0);
5828 kink->SetTPCDensity(dens11,0,1);
5829 kink->SetTPCDensity(dens00,1,0);
5830 kink->SetTPCDensity(dens01,1,1);
5831 kink->SetIndex(j,0);
5832 kink->SetIndex(i,1);
5835 if (mpt<1||kink->GetAngle(2)>0.1){
5836 // angle and densities not defined yet
5837 if (kink->GetTPCDensityFactor()<0.8) continue;
5838 if ((2-kink->GetTPCDensityFactor())*kink->GetDistance() >0.25) continue;
5839 if (kink->GetAngle(2)*ktrack0->GetP()<0.003) continue; //too small angle
5840 if (kink->GetAngle(2)>0.2&&kink->GetTPCDensityFactor()<1.15) continue;
5841 if (kink->GetAngle(2)>0.2&&kink->GetTPCDensity(0,1)>0.05) continue;
5843 Float_t criticalangle = track0->GetSigmaSnp2()+track0->GetSigmaTgl2();
5844 criticalangle+= track1->GetSigmaSnp2()+track1->GetSigmaTgl2();
5845 criticalangle= 3*TMath::Sqrt(criticalangle);
5846 if (criticalangle>0.02) criticalangle=0.02;
5847 if (kink->GetAngle(2)<criticalangle) continue;
5850 Int_t drow = Int_t(2.+0.5/(0.05+kink->GetAngle(2))); // overlap region defined
5851 Float_t shapesum =0;
5853 for ( Int_t row = row0-drow; row<row0+drow;row++){
5854 if (row<0) continue;
5855 if (row>155) continue;
5856 if (ktrack0->GetClusterPointer(row)){
5857 AliTPCTrackerPoint *point =ktrack0->GetTrackPoint(row);
5858 shapesum+=point->GetSigmaY()+point->GetSigmaZ();
5861 if (ktrack1->GetClusterPointer(row)){
5862 AliTPCTrackerPoint *point =ktrack1->GetTrackPoint(row);
5863 shapesum+=point->GetSigmaY()+point->GetSigmaZ();
5868 kink->SetShapeFactor(-1.);
5871 kink->SetShapeFactor(shapesum/sum);
5873 // esd->AddKink(kink);
5875 // kink->SetMother(paramm);
5876 //kink->SetDaughter(paramd);
5878 Double_t chi2P2 = paramm.GetParameter()[2]-paramd.GetParameter()[2];
5880 chi2P2/=paramm.GetCovariance()[5]+paramd.GetCovariance()[5];
5881 Double_t chi2P3 = paramm.GetParameter()[3]-paramd.GetParameter()[3];
5883 chi2P3/=paramm.GetCovariance()[9]+paramd.GetCovariance()[9];
5885 if (AliTPCReconstructor::StreamLevel()>1) {
5886 (*fDebugStreamer)<<"kinkLpt"<<
5894 if ( chi2P2+chi2P3<AliTPCReconstructor::GetRecoParam()->GetKinkAngleCutChi2(0)){
5898 kinks->AddLast(kink);
5904 // sort the kinks according quality - and refit them towards vertex
5906 Int_t nkinks = kinks->GetEntriesFast();
5907 Float_t *quality = new Float_t[nkinks];
5908 Int_t *indexes = new Int_t[nkinks];
5909 AliTPCseed **mothers = new AliTPCseed*[nkinks]; memset(mothers, 0, nkinks*sizeof(AliTPCseed*));
5910 AliTPCseed **daughters = new AliTPCseed*[nkinks]; memset(daughters, 0, nkinks*sizeof(AliTPCseed*));
5913 for (Int_t i=0;i<nkinks;i++){
5915 AliKink *kinkl = (AliKink*)kinks->At(i);
5917 // refit kinks towards vertex
5919 Int_t index0 = kinkl->GetIndex(0);
5920 Int_t index1 = kinkl->GetIndex(1);
5921 AliTPCseed * ktrack0 = (AliTPCseed*)array->At(index0);
5922 AliTPCseed * ktrack1 = (AliTPCseed*)array->At(index1);
5924 Int_t sumn=ktrack0->GetNumberOfClusters()+ktrack1->GetNumberOfClusters();
5926 // Refit Kink under if too small angle
5928 if (kinkl->GetAngle(2)<0.05){
5929 kinkl->SetTPCRow0(GetRowNumber(kinkl->GetR()));
5930 Int_t row0 = kinkl->GetTPCRow0();
5931 Int_t drow = Int_t(2.+0.5/(0.05+kinkl->GetAngle(2)));
5934 Int_t last = row0-drow;
5935 if (last<40) last=40;
5936 if (last<ktrack0->GetFirstPoint()+25) last = ktrack0->GetFirstPoint()+25;
5937 AliTPCseed* seed0 = ReSeed(ktrack0,last,kFALSE);
5940 Int_t first = row0+drow;
5941 if (first>130) first=130;
5942 if (first>ktrack1->GetLastPoint()-25) first = TMath::Max(ktrack1->GetLastPoint()-25,30);
5943 AliTPCseed* seed1 = ReSeed(ktrack1,first,kTRUE);
5945 if (seed0 && seed1){
5946 kinkl->SetStatus(1,8);
5947 if (RefitKink(*seed0,*seed1,*kinkl)) kinkl->SetStatus(1,9);
5948 row0 = GetRowNumber(kinkl->GetR());
5949 sumn = seed0->GetNumberOfClusters()+seed1->GetNumberOfClusters();
5950 mothers[i] = new ( NextFreeSeed() ) AliTPCseed(*seed0);
5951 mothers[i]->SetPoolID(fLastSeedID);
5952 daughters[i] = new (NextFreeSeed() ) AliTPCseed(*seed1);
5953 daughters[i]->SetPoolID(fLastSeedID);
5956 delete kinks->RemoveAt(i);
5957 if (seed0) MarkSeedFree( seed0 );
5958 if (seed1) MarkSeedFree( seed1 );
5961 if (kinkl->GetDistance()>0.5 || kinkl->GetR()<110 || kinkl->GetR()>240) {
5962 delete kinks->RemoveAt(i);
5963 if (seed0) MarkSeedFree( seed0 );
5964 if (seed1) MarkSeedFree( seed1 );
5968 MarkSeedFree( seed0 );
5969 MarkSeedFree( seed1 );
5972 if (kinkl) quality[i] = 160*((0.1+kinkl->GetDistance())*(2.-kinkl->GetTPCDensityFactor()))/(sumn+40.); //the longest -clossest will win
5974 TMath::Sort(nkinks,quality,indexes,kFALSE);
5976 //remove double find kinks
5978 for (Int_t ikink0=1;ikink0<nkinks;ikink0++){
5979 AliKink * kink0 = (AliKink*) kinks->At(indexes[ikink0]);
5980 if (!kink0) continue;
5982 for (Int_t ikink1=0;ikink1<ikink0;ikink1++){
5983 kink0 = (AliKink*) kinks->At(indexes[ikink0]);
5984 if (!kink0) continue;
5985 AliKink * kink1 = (AliKink*) kinks->At(indexes[ikink1]);
5986 if (!kink1) continue;
5987 // if not close kink continue
5988 if (TMath::Abs(kink1->GetPosition()[2]-kink0->GetPosition()[2])>10) continue;
5989 if (TMath::Abs(kink1->GetPosition()[1]-kink0->GetPosition()[1])>10) continue;
5990 if (TMath::Abs(kink1->GetPosition()[0]-kink0->GetPosition()[0])>10) continue;
5992 AliTPCseed &mother0 = *mothers[indexes[ikink0]];
5993 AliTPCseed &daughter0 = *daughters[indexes[ikink0]];
5994 AliTPCseed &mother1 = *mothers[indexes[ikink1]];
5995 AliTPCseed &daughter1 = *daughters[indexes[ikink1]];
5996 Int_t row0 = (kink0->GetTPCRow0()+kink1->GetTPCRow0())/2;
6005 for (Int_t i=0;i<row0;i++){
6006 if (mother0.GetClusterIndex(i)>0 && mother1.GetClusterIndex(i)>0){
6009 if (mother0.GetClusterIndex(i)==mother1.GetClusterIndex(i)){
6016 for (Int_t i=row0;i<158;i++){
6017 //if (daughter0.GetClusterIndex(i)>0 && daughter0.GetClusterIndex(i)>0){ // RS: Bug?
6018 if (daughter0.GetClusterIndex(i)>0 && daughter1.GetClusterIndex(i)>0){
6021 if (mother0.GetClusterIndex(i)==mother1.GetClusterIndex(i)){
6027 Float_t ratio = Float_t(same+1)/Float_t(both+1);
6028 Float_t ratiom = Float_t(samem+1)/Float_t(bothm+1);
6029 Float_t ratiod = Float_t(samed+1)/Float_t(bothd+1);
6030 if (ratio>0.3 && ratiom>0.5 &&ratiod>0.5) {
6031 Int_t sum0 = mother0.GetNumberOfClusters()+daughter0.GetNumberOfClusters();
6032 Int_t sum1 = mother1.GetNumberOfClusters()+daughter1.GetNumberOfClusters();
6034 shared[kink0->GetIndex(0)]= kTRUE;
6035 shared[kink0->GetIndex(1)]= kTRUE;
6036 delete kinks->RemoveAt(indexes[ikink0]);
6040 shared[kink1->GetIndex(0)]= kTRUE;
6041 shared[kink1->GetIndex(1)]= kTRUE;
6042 delete kinks->RemoveAt(indexes[ikink1]);
6049 for (Int_t i=0;i<nkinks;i++){
6050 AliKink * kinkl = (AliKink*) kinks->At(indexes[i]);
6051 if (!kinkl) continue;
6052 kinkl->SetTPCRow0(GetRowNumber(kinkl->GetR()));
6053 Int_t index0 = kinkl->GetIndex(0);
6054 Int_t index1 = kinkl->GetIndex(1);
6055 if (circular[index0]||(circular[index1]&&kinkl->GetDistance()>0.2)) continue;
6056 kinkl->SetMultiple(usage[index0],0);
6057 kinkl->SetMultiple(usage[index1],1);
6058 if (kinkl->GetMultiple()[0]+kinkl->GetMultiple()[1]>2) continue;
6059 if (kinkl->GetMultiple()[0]+kinkl->GetMultiple()[1]>0 && quality[indexes[i]]>0.2) continue;
6060 if (kinkl->GetMultiple()[0]+kinkl->GetMultiple()[1]>0 && kinkl->GetDistance()>0.2) continue;
6061 if (circular[index0]||(circular[index1]&&kinkl->GetDistance()>0.1)) continue;
6063 AliTPCseed * ktrack0 = (AliTPCseed*)array->At(index0);
6064 AliTPCseed * ktrack1 = (AliTPCseed*)array->At(index1);
6065 if (!ktrack0 || !ktrack1) continue;
6066 Int_t index = esd->AddKink(kinkl);
6069 if ( ktrack0->GetKinkIndex(0)==0 && ktrack1->GetKinkIndex(0)==0) { //best kink
6070 if (mothers[indexes[i]]->GetNumberOfClusters()>20 && daughters[indexes[i]]->GetNumberOfClusters()>20 &&
6071 (mothers[indexes[i]]->GetNumberOfClusters()+daughters[indexes[i]]->GetNumberOfClusters())>100){
6072 *ktrack0 = *mothers[indexes[i]];
6073 *ktrack1 = *daughters[indexes[i]];
6077 ktrack0->SetKinkIndex(usage[index0],-(index+1));
6078 ktrack1->SetKinkIndex(usage[index1], (index+1));
6083 // Remove tracks corresponding to shared kink's
6085 for (Int_t i=0;i<nentries;i++){
6086 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
6087 if (!track0) continue;
6088 if (track0->GetKinkIndex(0)!=0) continue;
6089 if (shared[i]) MarkSeedFree( array->RemoveAt(i) );
6094 RemoveUsed2(array,0.5,0.4,30);
6096 for (Int_t i=0;i<nentries;i++){
6097 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
6098 if (!track0) continue;
6099 track0->CookdEdx(0.02,0.6);
6103 for (Int_t i=0;i<nentries;i++){
6104 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
6105 if (!track0) continue;
6106 if (track0->Pt()<1.4) continue;
6107 //remove double high momenta tracks - overlapped with kink candidates
6110 for (Int_t icl=track0->GetFirstPoint();icl<track0->GetLastPoint(); icl++){
6111 if (track0->GetClusterPointer(icl)!=0){
6113 if (track0->GetClusterPointer(icl)->IsUsed(10)) ishared++;
6116 if (Float_t(ishared+1)/Float_t(all+1)>0.5) {
6117 MarkSeedFree( array->RemoveAt(i) );
6121 if (track0->GetKinkIndex(0)!=0) continue;
6122 if (track0->GetNumberOfClusters()<80) continue;
6124 AliTPCseed *pmother = new( NextFreeSeed() ) AliTPCseed();
6125 pmother->SetPoolID(fLastSeedID);
6126 AliTPCseed *pdaughter = new( NextFreeSeed() ) AliTPCseed();
6127 pdaughter->SetPoolID(fLastSeedID);
6128 AliKink *pkink = new AliKink;
6130 AliTPCseed & mother = *pmother;
6131 AliTPCseed & daughter = *pdaughter;
6132 AliKink & kinkl = *pkink;
6133 if (CheckKinkPoint(track0,mother,daughter, kinkl)){
6134 if (mother.GetNumberOfClusters()<30||daughter.GetNumberOfClusters()<20) {
6135 MarkSeedFree( pmother );
6136 MarkSeedFree( pdaughter );
6138 continue; //too short tracks
6140 if (mother.Pt()<1.4) {
6141 MarkSeedFree( pmother );
6142 MarkSeedFree( pdaughter );
6146 Int_t row0= kinkl.GetTPCRow0();
6147 if (kinkl.GetDistance()>0.5 || kinkl.GetR()<110. || kinkl.GetR()>240.) {
6148 MarkSeedFree( pmother );
6149 MarkSeedFree( pdaughter );
6154 Int_t index = esd->AddKink(&kinkl);
6155 mother.SetKinkIndex(0,-(index+1));
6156 daughter.SetKinkIndex(0,index+1);
6157 if (mother.GetNumberOfClusters()>50) {
6158 MarkSeedFree( array->RemoveAt(i) );
6159 AliTPCseed* mtc = new( NextFreeSeed() ) AliTPCseed(mother);
6160 mtc->SetPoolID(fLastSeedID);
6161 array->AddAt(mtc,i);
6164 AliTPCseed* mtc = new( NextFreeSeed() ) AliTPCseed(mother);
6165 mtc->SetPoolID(fLastSeedID);
6166 array->AddLast(mtc);
6168 AliTPCseed* dtc = new( NextFreeSeed() ) AliTPCseed(daughter);
6169 dtc->SetPoolID(fLastSeedID);
6170 array->AddLast(dtc);
6171 for (Int_t icl=0;icl<row0;icl++) {
6172 if (mother.GetClusterPointer(icl)) mother.GetClusterPointer(icl)->Use(20);
6175 for (Int_t icl=row0;icl<158;icl++) {
6176 if (daughter.GetClusterPointer(icl)) daughter.GetClusterPointer(icl)->Use(20);
6180 MarkSeedFree( pmother );
6181 MarkSeedFree( pdaughter );
6185 delete [] daughters;
6207 AliInfo(Form("Ncandidates=\t%d\t%d\t%d\t%d\n",esd->GetNumberOfKinks(),ncandidates,ntracks,nall));
6212 Int_t AliTPCtrackerMI::RefitKink(AliTPCseed &mother, AliTPCseed &daughter, const AliESDkink &knk)
6215 // refit kink towards to the vertex
6218 AliKink &kink=(AliKink &)knk;
6220 Int_t row0 = GetRowNumber(kink.GetR());
6221 FollowProlongation(mother,0);
6222 mother.Reset(kFALSE);
6224 FollowProlongation(daughter,row0);
6225 daughter.Reset(kFALSE);
6226 FollowBackProlongation(daughter,158);
6227 daughter.Reset(kFALSE);
6228 Int_t first = TMath::Max(row0-20,30);
6229 Int_t last = TMath::Min(row0+20,140);
6231 const Int_t kNdiv =5;
6232 AliTPCseed param0[kNdiv]; // parameters along the track
6233 AliTPCseed param1[kNdiv]; // parameters along the track
6234 AliKink kinks[kNdiv]; // corresponding kink parameters
6237 for (Int_t irow=0; irow<kNdiv;irow++){
6238 rows[irow] = first +((last-first)*irow)/(kNdiv-1);
6240 // store parameters along the track
6242 for (Int_t irow=0;irow<kNdiv;irow++){
6243 FollowBackProlongation(mother, rows[irow]);
6244 FollowProlongation(daughter,rows[kNdiv-1-irow]);
6245 param0[irow] = mother;
6246 param1[kNdiv-1-irow] = daughter;
6250 for (Int_t irow=0; irow<kNdiv-1;irow++){
6251 if (param0[irow].GetNumberOfClusters()<kNdiv||param1[irow].GetNumberOfClusters()<kNdiv) continue;
6252 kinks[irow].SetMother(param0[irow]);
6253 kinks[irow].SetDaughter(param1[irow]);
6254 kinks[irow].Update();
6257 // choose kink with best "quality"
6259 Double_t mindist = 10000;
6260 for (Int_t irow=0;irow<kNdiv;irow++){
6261 if (param0[irow].GetNumberOfClusters()<20||param1[irow].GetNumberOfClusters()<20) continue;
6262 if (TMath::Abs(kinks[irow].GetR())>240.) continue;
6263 if (TMath::Abs(kinks[irow].GetR())<100.) continue;
6265 Float_t normdist = TMath::Abs(param0[irow].GetX()-kinks[irow].GetR())*(0.1+kink.GetDistance());
6266 normdist/= (param0[irow].GetNumberOfClusters()+param1[irow].GetNumberOfClusters()+40.);
6267 if (normdist < mindist){
6273 if (index==-1) return 0;
6276 param0[index].Reset(kTRUE);
6277 FollowProlongation(param0[index],0);
6279 mother = param0[index];
6280 daughter = param1[index]; // daughter in vertex
6282 kink.SetMother(mother);
6283 kink.SetDaughter(daughter);
6285 kink.SetTPCRow0(GetRowNumber(kink.GetR()));
6286 kink.SetTPCncls(param0[index].GetNumberOfClusters(),0);
6287 kink.SetTPCncls(param1[index].GetNumberOfClusters(),1);
6288 kink.SetLabel(CookLabel(&mother,0.4, 0,kink.GetTPCRow0()),0);
6289 kink.SetLabel(CookLabel(&daughter,0.4, kink.GetTPCRow0(),160),1);
6290 mother.SetLabel(kink.GetLabel(0));
6291 daughter.SetLabel(kink.GetLabel(1));
6297 void AliTPCtrackerMI::UpdateKinkQualityM(AliTPCseed * seed){
6299 // update Kink quality information for mother after back propagation
6301 if (seed->GetKinkIndex(0)>=0) return;
6302 for (Int_t ikink=0;ikink<3;ikink++){
6303 Int_t index = seed->GetKinkIndex(ikink);
6304 if (index>=0) break;
6305 index = TMath::Abs(index)-1;
6306 AliESDkink * kink = fEvent->GetKink(index);
6307 kink->SetTPCDensity(-1,0,0);
6308 kink->SetTPCDensity(1,0,1);
6310 Int_t row0 = kink->GetTPCRow0() - 2 - Int_t( 0.5/ (0.05+kink->GetAngle(2)));
6311 if (row0<15) row0=15;
6313 Int_t row1 = kink->GetTPCRow0() + 2 + Int_t( 0.5/ (0.05+kink->GetAngle(2)));
6314 if (row1>145) row1=145;
6316 Int_t found,foundable,shared;
6317 seed->GetClusterStatistic(0,row0, found, foundable,shared,kFALSE);
6318 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),0,0);
6319 seed->GetClusterStatistic(row1,155, found, foundable,shared,kFALSE);
6320 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),0,1);
6325 void AliTPCtrackerMI::UpdateKinkQualityD(AliTPCseed * seed){
6327 // update Kink quality information for daughter after refit
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,1,0);
6336 kink->SetTPCDensity(-1,1,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),1,0);
6347 seed->GetClusterStatistic(row1,155, found, foundable,shared,kFALSE);
6348 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),1,1);
6354 Int_t AliTPCtrackerMI::CheckKinkPoint(AliTPCseed*seed,AliTPCseed &mother, AliTPCseed &daughter, const AliESDkink &knk)
6357 // check kink point for given track
6358 // if return value=0 kink point not found
6359 // otherwise seed0 correspond to mother particle
6360 // seed1 correspond to daughter particle
6361 // kink parameter of kink point
6362 AliKink &kink=(AliKink &)knk;
6364 Int_t middlerow = (seed->GetFirstPoint()+seed->GetLastPoint())/2;
6365 Int_t first = seed->GetFirstPoint();
6366 Int_t last = seed->GetLastPoint();
6367 if (last-first<20) return 0; // shortest length - 2*30 = 60 pad-rows
6370 AliTPCseed *seed1 = ReSeed(seed,middlerow+20, kTRUE); //middle of chamber
6371 if (!seed1) return 0;
6372 FollowProlongation(*seed1,seed->GetLastPoint()-20);
6373 seed1->Reset(kTRUE);
6374 FollowProlongation(*seed1,158);
6375 seed1->Reset(kTRUE);
6376 last = seed1->GetLastPoint();
6378 AliTPCseed *seed0 = new( NextFreeSeed() ) AliTPCseed(*seed);
6379 seed0->SetPoolID(fLastSeedID);
6380 seed0->Reset(kFALSE);
6383 AliTPCseed param0[20]; // parameters along the track
6384 AliTPCseed param1[20]; // parameters along the track
6385 AliKink kinks[20]; // corresponding kink parameters
6387 for (Int_t irow=0; irow<20;irow++){
6388 rows[irow] = first +((last-first)*irow)/19;
6390 // store parameters along the track
6392 for (Int_t irow=0;irow<20;irow++){
6393 FollowBackProlongation(*seed0, rows[irow]);
6394 FollowProlongation(*seed1,rows[19-irow]);
6395 param0[irow] = *seed0;
6396 param1[19-irow] = *seed1;
6400 for (Int_t irow=0; irow<19;irow++){
6401 kinks[irow].SetMother(param0[irow]);
6402 kinks[irow].SetDaughter(param1[irow]);
6403 kinks[irow].Update();
6406 // choose kink with biggest change of angle
6408 Double_t maxchange= 0;
6409 for (Int_t irow=1;irow<19;irow++){
6410 if (TMath::Abs(kinks[irow].GetR())>240.) continue;
6411 if (TMath::Abs(kinks[irow].GetR())<110.) continue;
6412 Float_t quality = TMath::Abs(kinks[irow].GetAngle(2))/(3.+TMath::Abs(kinks[irow].GetR()-param0[irow].GetX()));
6413 if ( quality > maxchange){
6414 maxchange = quality;
6419 MarkSeedFree( seed0 );
6420 MarkSeedFree( seed1 );
6421 if (index<0) return 0;
6423 Int_t row0 = GetRowNumber(kinks[index].GetR()); //row 0 estimate
6424 seed0 = new( NextFreeSeed() ) AliTPCseed(param0[index]);
6425 seed0->SetPoolID(fLastSeedID);
6426 seed1 = new( NextFreeSeed() ) AliTPCseed(param1[index]);
6427 seed1->SetPoolID(fLastSeedID);
6428 seed0->Reset(kFALSE);
6429 seed1->Reset(kFALSE);
6430 seed0->ResetCovariance(10.);
6431 seed1->ResetCovariance(10.);
6432 FollowProlongation(*seed0,0);
6433 FollowBackProlongation(*seed1,158);
6434 mother = *seed0; // backup mother at position 0
6435 seed0->Reset(kFALSE);
6436 seed1->Reset(kFALSE);
6437 seed0->ResetCovariance(10.);
6438 seed1->ResetCovariance(10.);
6440 first = TMath::Max(row0-20,0);
6441 last = TMath::Min(row0+20,158);
6443 for (Int_t irow=0; irow<20;irow++){
6444 rows[irow] = first +((last-first)*irow)/19;
6446 // store parameters along the track
6448 for (Int_t irow=0;irow<20;irow++){
6449 FollowBackProlongation(*seed0, rows[irow]);
6450 FollowProlongation(*seed1,rows[19-irow]);
6451 param0[irow] = *seed0;
6452 param1[19-irow] = *seed1;
6456 for (Int_t irow=0; irow<19;irow++){
6457 kinks[irow].SetMother(param0[irow]);
6458 kinks[irow].SetDaughter(param1[irow]);
6459 // param0[irow].Dump();
6460 //param1[irow].Dump();
6461 kinks[irow].Update();
6464 // choose kink with biggest change of angle
6467 for (Int_t irow=0;irow<20;irow++){
6468 if (TMath::Abs(kinks[irow].GetR())>250.) continue;
6469 if (TMath::Abs(kinks[irow].GetR())<90.) continue;
6470 Float_t quality = TMath::Abs(kinks[irow].GetAngle(2))/(3.+TMath::Abs(kinks[irow].GetR()-param0[irow].GetX()));
6471 if ( quality > maxchange){
6472 maxchange = quality;
6479 if (index==-1 || param0[index].GetNumberOfClusters()+param1[index].GetNumberOfClusters()<100){
6480 MarkSeedFree( seed0 );
6481 MarkSeedFree( seed1 );
6485 // Float_t anglesigma = TMath::Sqrt(param0[index].fC22+param0[index].fC33+param1[index].fC22+param1[index].fC33);
6487 kink.SetMother(param0[index]);
6488 kink.SetDaughter(param1[index]);
6491 Double_t chi2P2 = param0[index].GetParameter()[2]-param1[index].GetParameter()[2];
6493 chi2P2/=param0[index].GetCovariance()[5]+param1[index].GetCovariance()[5];
6494 Double_t chi2P3 = param0[index].GetParameter()[3]-param1[index].GetParameter()[3];
6496 chi2P3/=param0[index].GetCovariance()[9]+param1[index].GetCovariance()[9];
6498 if (AliTPCReconstructor::StreamLevel()>1) {
6499 (*fDebugStreamer)<<"kinkHpt"<<
6502 "p0.="<<¶m0[index]<<
6503 "p1.="<<¶m1[index]<<
6507 if ( chi2P2+chi2P3<AliTPCReconstructor::GetRecoParam()->GetKinkAngleCutChi2(0)){
6508 MarkSeedFree( seed0 );
6509 MarkSeedFree( seed1 );
6514 row0 = GetRowNumber(kink.GetR());
6515 kink.SetTPCRow0(row0);
6516 kink.SetLabel(CookLabel(seed0,0.5,0,row0),0);
6517 kink.SetLabel(CookLabel(seed1,0.5,row0,158),1);
6518 kink.SetIndex(-10,0);
6519 kink.SetIndex(int(param0[index].GetNumberOfClusters()+param1[index].GetNumberOfClusters()),1);
6520 kink.SetTPCncls(param0[index].GetNumberOfClusters(),0);
6521 kink.SetTPCncls(param1[index].GetNumberOfClusters(),1);
6524 // new (&mother) AliTPCseed(param0[index]);
6525 daughter = param1[index];
6526 daughter.SetLabel(kink.GetLabel(1));
6527 param0[index].Reset(kTRUE);
6528 FollowProlongation(param0[index],0);
6529 mother = param0[index];
6530 mother.SetLabel(kink.GetLabel(0));
6531 if ( chi2P2+chi2P3<AliTPCReconstructor::GetRecoParam()->GetKinkAngleCutChi2(1)){
6534 MarkSeedFree( seed0 );
6535 MarkSeedFree( seed1 );
6543 AliTPCseed* AliTPCtrackerMI::ReSeed(AliTPCseed *t)
6546 // reseed - refit - track
6549 // Int_t last = fSectors->GetNRows()-1;
6551 if (fSectors == fOuterSec){
6552 first = TMath::Max(first, t->GetFirstPoint()-fInnerSec->GetNRows());
6556 first = t->GetFirstPoint();
6558 AliTPCseed * seed = MakeSeed(t,0.1,0.5,0.9);
6559 FollowBackProlongation(*t,fSectors->GetNRows()-1);
6561 FollowProlongation(*t,first);
6571 //_____________________________________________________________________________
6572 Int_t AliTPCtrackerMI::ReadSeeds(const TFile *inp) {
6573 //-----------------------------------------------------------------
6574 // This function reades track seeds.
6575 //-----------------------------------------------------------------
6576 TDirectory *savedir=gDirectory;
6578 TFile *in=(TFile*)inp;
6579 if (!in->IsOpen()) {
6580 cerr<<"AliTPCtrackerMI::ReadSeeds(): input file is not open !\n";
6585 TTree *seedTree=(TTree*)in->Get("Seeds");
6587 cerr<<"AliTPCtrackerMI::ReadSeeds(): ";
6588 cerr<<"can't get a tree with track seeds !\n";
6591 AliTPCtrack *seed=new AliTPCtrack;
6592 seedTree->SetBranchAddress("tracks",&seed);
6594 if (fSeeds==0) fSeeds=new TObjArray(15000);
6596 Int_t n=(Int_t)seedTree->GetEntries();
6597 for (Int_t i=0; i<n; i++) {
6598 seedTree->GetEvent(i);
6599 AliTPCseed* sdc = new( NextFreeSeed() ) AliTPCseed(*seed/*,seed->GetAlpha()*/);
6600 sdc->SetPoolID(fLastSeedID);
6601 fSeeds->AddLast(sdc);
6604 delete seed; // RS: this seed is not from the pool, delete it !!!
6610 Int_t AliTPCtrackerMI::Clusters2Tracks (AliESDEvent *const esd)
6613 // clusters to tracks
6614 if (fSeeds) DeleteSeeds();
6615 else ResetSeedsPool();
6618 AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ;
6619 transform->SetCurrentTimeStamp( esd->GetTimeStamp());
6620 transform->SetCurrentRun(esd->GetRunNumber());
6623 if (!fSeeds) return 1;
6625 if (AliTPCReconstructor::StreamLevel()>3) DumpClusters(0,fSeeds);
6631 //_____________________________________________________________________________
6632 Int_t AliTPCtrackerMI::Clusters2Tracks() {
6633 //-----------------------------------------------------------------
6634 // This is a track finder.
6635 //-----------------------------------------------------------------
6636 TDirectory *savedir=gDirectory;
6640 fSeeds = Tracking();
6643 Info("Clusters2Tracks","Time for tracking: \t");timer.Print();timer.Start();
6645 //activate again some tracks
6646 for (Int_t i=0; i<fSeeds->GetEntriesFast(); i++) {
6647 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
6649 Int_t nc=t.GetNumberOfClusters();
6651 MarkSeedFree( fSeeds->RemoveAt(i) );
6655 if (pt->GetRemoval()==10) {
6656 if (pt->GetDensityFirst(20)>0.8 || pt->GetDensityFirst(30)>0.8 || pt->GetDensityFirst(40)>0.7)
6657 pt->Desactivate(10); // make track again active // MvL: should be 0 ?
6659 pt->Desactivate(20);
6660 MarkSeedFree( fSeeds->RemoveAt(i) );
6665 RemoveUsed2(fSeeds,0.85,0.85,0);
6666 if (AliTPCReconstructor::GetRecoParam()->GetDoKinks()) FindKinks(fSeeds,fEvent);
6667 //FindCurling(fSeeds, fEvent,0);
6668 if (AliTPCReconstructor::StreamLevel()>5) FindMultiMC(fSeeds, fEvent,-1); // find multi found tracks
6669 RemoveUsed2(fSeeds,0.5,0.4,20);
6670 FindSplitted(fSeeds, fEvent,0); // find multi found tracks
6671 if (AliTPCReconstructor::StreamLevel()>5) FindMultiMC(fSeeds, fEvent,0); // find multi found tracks
6674 // // refit short tracks
6676 Int_t nseed=fSeeds->GetEntriesFast();
6679 for (Int_t i=0; i<nseed; i++) {
6680 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
6682 Int_t nc=t.GetNumberOfClusters();
6684 MarkSeedFree( fSeeds->RemoveAt(i) );
6687 CookLabel(pt,0.1); //For comparison only
6688 //if ((pt->IsActive() || (pt->fRemoval==10) )&& nc>50 &&pt->GetNumberOfClusters()>0.4*pt->fNFoundable){
6689 if ((pt->IsActive() || (pt->GetRemoval()==10) )){
6691 if (fDebug>0) cerr<<found<<'\r';
6695 MarkSeedFree( fSeeds->RemoveAt(i) );
6699 //RemoveOverlap(fSeeds,0.99,7,kTRUE);
6701 //RemoveUsed(fSeeds,0.9,0.9,6);
6703 nseed=fSeeds->GetEntriesFast();
6705 for (Int_t i=0; i<nseed; i++) {
6706 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
6708 Int_t nc=t.GetNumberOfClusters();
6710 MarkSeedFree( fSeeds->RemoveAt(i) );
6714 t.CookdEdx(0.02,0.6);
6715 // CheckKinkPoint(&t,0.05);
6716 //if ((pt->IsActive() || (pt->fRemoval==10) )&& nc>50 &&pt->GetNumberOfClusters()>0.4*pt->fNFoundable){
6717 if ((pt->IsActive() || (pt->GetRemoval()==10) )){
6725 MarkSeedFree( fSeeds->RemoveAt(i) );
6726 //AliTPCseed * seed1 = ReSeed(pt,0.05,0.5,1);
6728 // FollowProlongation(*seed1,0);
6729 // Int_t n = seed1->GetNumberOfClusters();
6730 // printf("fP4\t%f\t%f\n",seed1->GetC(),pt->GetC());
6731 // printf("fN\t%d\t%d\n", seed1->GetNumberOfClusters(),pt->GetNumberOfClusters());
6734 //AliTPCseed * seed2 = ReSeed(pt,0.95,0.5,0.05);
6738 SortTracks(fSeeds, 1);
6742 PrepareForBackProlongation(fSeeds,5.);
6743 PropagateBack(fSeeds);
6744 printf("Time for back propagation: \t");timer.Print();timer.Start();
6748 PrepareForProlongation(fSeeds,5.);
6749 PropagateForard2(fSeeds);
6751 printf("Time for FORWARD propagation: \t");timer.Print();timer.Start();
6752 // RemoveUsed(fSeeds,0.7,0.7,6);
6753 //RemoveOverlap(fSeeds,0.9,7,kTRUE);
6755 nseed=fSeeds->GetEntriesFast();
6757 for (Int_t i=0; i<nseed; i++) {
6758 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
6760 Int_t nc=t.GetNumberOfClusters();
6762 MarkSeedFree( fSeeds->RemoveAt(i) );
6765 t.CookdEdx(0.02,0.6);
6766 // CookLabel(pt,0.1); //For comparison only
6767 //if ((pt->IsActive() || (pt->fRemoval==10) )&& nc>50 &&pt->GetNumberOfClusters()>0.4*pt->fNFoundable){
6768 if ((pt->IsActive() || (pt->fRemoval==10) )){
6769 cerr<<found++<<'\r';
6772 MarkSeedFree( fSeeds->RemoveAt(i) );
6777 // fNTracks = found;
6779 Info("Clusters2Tracks","Time for overlap removal, track writing and dedx cooking: \t"); timer.Print();timer.Start();
6782 // cerr<<"Number of found tracks : "<<"\t"<<found<<endl;
6783 Info("Clusters2Tracks","Number of found tracks %d",found);
6785 // UnloadClusters();
6790 void AliTPCtrackerMI::Tracking(TObjArray * arr)
6793 // tracking of the seeds
6796 fSectors = fOuterSec;
6797 ParallelTracking(arr,150,63);
6798 fSectors = fOuterSec;
6799 ParallelTracking(arr,63,0);
6802 TObjArray * AliTPCtrackerMI::Tracking(Int_t seedtype, Int_t i1, Int_t i2, Float_t cuts[4], Float_t dy, Int_t dsec)
6807 static TObjArray arrTracks;
6808 TObjArray * arr = &arrTracks;
6810 fSectors = fOuterSec;
6813 for (Int_t sec=0;sec<fkNOS;sec++){
6814 if (seedtype==3) MakeSeeds3(arr,sec,i1,i2,cuts,dy, dsec);
6815 if (seedtype==4) MakeSeeds5(arr,sec,i1,i2,cuts,dy);
6816 if (seedtype==2) MakeSeeds2(arr,sec,i1,i2,cuts,dy);
6819 Info("Tracking","\nSeeding - %d\t%d\t%d\t%d\n",seedtype,i1,i2,arr->GetEntriesFast());
6831 TObjArray * AliTPCtrackerMI::Tracking()
6835 if (AliTPCReconstructor::GetRecoParam()->GetSpecialSeeding()) return TrackingSpecial();
6838 Int_t nup=fOuterSec->GetNRows()+fInnerSec->GetNRows();
6840 TObjArray * seeds = new TObjArray;
6842 Int_t fLastSeedRowSec=AliTPCReconstructor::GetRecoParam()->GetLastSeedRowSec();
6843 Int_t gapPrim = AliTPCReconstructor::GetRecoParam()->GetSeedGapPrim();
6844 Int_t gapSec = AliTPCReconstructor::GetRecoParam()->GetSeedGapSec();
6852 Float_t fnumber = 3.0;
6853 Float_t fdensity = 3.0;
6858 for (Int_t delta = 0; delta<18; delta+=gapPrim){
6862 arr = Tracking(3,nup-1-delta,nup-1-delta-gap,cuts,-1,1);
6863 SumTracks(seeds,arr);
6864 SignClusters(seeds,fnumber,fdensity);
6866 for (Int_t i=2;i<6;i+=2){
6867 // seed high pt tracks
6870 arr = Tracking(3,nup-i-delta,nup-i-delta-gap,cuts,-1,0);
6871 SumTracks(seeds,arr);
6872 SignClusters(seeds,fnumber,fdensity);
6877 // RemoveUsed(seeds,0.9,0.9,1);
6878 // UnsignClusters();
6879 // SignClusters(seeds,fnumber,fdensity);
6883 for (Int_t delta = 20; delta<120; delta+=gapPrim){
6885 // seed high pt tracks
6889 arr = Tracking(3,nup-delta,nup-delta-gap,cuts,-1);
6890 SumTracks(seeds,arr);
6891 SignClusters(seeds,fnumber,fdensity);
6896 arr = Tracking(3,nup-delta-5,nup-delta-5-gap,cuts,-1);
6897 SumTracks(seeds,arr);
6898 SignClusters(seeds,fnumber,fdensity);
6909 Info("Tracking()","\n\nPrimary seeding\t%d\n\n",seeds->GetEntriesFast());
6913 // RemoveUsed(seeds,0.75,0.75,1);
6915 //SignClusters(seeds,fnumber,fdensity);
6924 arr = Tracking(4,nup-1,nup-1-gap,cuts,-1);
6925 SumTracks(seeds,arr);
6926 SignClusters(seeds,fnumber,fdensity);
6928 arr = Tracking(4,nup-2,nup-2-gap,cuts,-1);
6929 SumTracks(seeds,arr);
6930 SignClusters(seeds,fnumber,fdensity);
6932 arr = Tracking(4,nup-3,nup-3-gap,cuts,-1);
6933 SumTracks(seeds,arr);
6934 SignClusters(seeds,fnumber,fdensity);
6936 arr = Tracking(4,nup-5,nup-5-gap,cuts,-1);
6937 SumTracks(seeds,arr);
6938 SignClusters(seeds,fnumber,fdensity);
6940 arr = Tracking(4,nup-7,nup-7-gap,cuts,-1);
6941 SumTracks(seeds,arr);
6942 SignClusters(seeds,fnumber,fdensity);
6945 arr = Tracking(4,nup-9,nup-9-gap,cuts,-1);
6946 SumTracks(seeds,arr);
6947 SignClusters(seeds,fnumber,fdensity);
6951 for (Int_t delta = 9; delta<30; delta+=gapSec){
6957 arr = Tracking(4,nup-1-delta,nup-1-delta-gap,cuts,-1);
6958 SumTracks(seeds,arr);
6959 SignClusters(seeds,fnumber,fdensity);
6961 arr = Tracking(4,nup-3-delta,nup-5-delta-gap,cuts,4);
6962 SumTracks(seeds,arr);
6963 SignClusters(seeds,fnumber,fdensity);
6976 for (Int_t delta = 30; delta<fLastSeedRowSec; delta+=gapSec){
6982 arr = Tracking(4,nup-1-delta,nup-1-delta-gap,cuts,-1);
6983 SumTracks(seeds,arr);
6984 SignClusters(seeds,fnumber,fdensity);
6986 arr = Tracking(4,nup-5-delta,nup-5-delta-gap,cuts,5 );
6987 SumTracks(seeds,arr);
6988 SignClusters(seeds,fnumber,fdensity);
6992 Info("Tracking()","\n\nSecondary seeding\t%d\n\n",seeds->GetEntriesFast());
7003 TObjArray * AliTPCtrackerMI::TrackingSpecial()
7006 // seeding adjusted for laser and cosmic tests - short tracks with big inclination angle
7007 // no primary vertex seeding tried
7011 Int_t nup=fOuterSec->GetNRows()+fInnerSec->GetNRows();
7013 TObjArray * seeds = new TObjArray;
7018 Float_t fnumber = 3.0;
7019 Float_t fdensity = 3.0;
7022 cuts[0] = AliTPCReconstructor::GetRecoParam()->GetMaxC(); // max curvature
7023 cuts[1] = 3.5; // max tan(phi) angle for seeding
7024 cuts[2] = 3.; // not used (cut on z primary vertex)
7025 cuts[3] = 3.5; // max tan(theta) angle for seeding
7027 for (Int_t delta = 0; nup-delta-gap-1>0; delta+=3){
7029 arr = Tracking(4,nup-1-delta,nup-1-delta-gap,cuts,-1);
7030 SumTracks(seeds,arr);
7031 SignClusters(seeds,fnumber,fdensity);
7035 Info("Tracking()","\n\nSecondary seeding\t%d\n\n",seeds->GetEntriesFast());
7046 void AliTPCtrackerMI::SumTracks(TObjArray *arr1,TObjArray *&arr2)
7049 //sum tracks to common container
7050 //remove suspicious tracks
7051 // RS: Attention: supplied tracks come in the static array, don't delete them
7052 Int_t nseed = arr2->GetEntriesFast();
7053 for (Int_t i=0;i<nseed;i++){
7054 AliTPCseed *pt=(AliTPCseed*)arr2->UncheckedAt(i);
7057 // remove tracks with too big curvature
7059 if (TMath::Abs(pt->GetC())>AliTPCReconstructor::GetRecoParam()->GetMaxC()){
7060 MarkSeedFree( arr2->RemoveAt(i) );
7063 // REMOVE VERY SHORT TRACKS
7064 if (pt->GetNumberOfClusters()<20){
7065 MarkSeedFree( arr2->RemoveAt(i) );
7068 // NORMAL ACTIVE TRACK
7069 if (pt->IsActive()){
7070 arr1->AddLast(arr2->RemoveAt(i));
7073 //remove not usable tracks
7074 if (pt->GetRemoval()!=10){
7075 MarkSeedFree( arr2->RemoveAt(i) );
7079 // ENABLE ONLY ENOUGH GOOD STOPPED TRACKS
7080 if (pt->GetDensityFirst(20)>0.8 || pt->GetDensityFirst(30)>0.8 || pt->GetDensityFirst(40)>0.7)
7081 arr1->AddLast(arr2->RemoveAt(i));
7083 MarkSeedFree( arr2->RemoveAt(i) );
7087 // delete arr2; arr2 = 0; // RS: this is static array, don't delete it
7092 void AliTPCtrackerMI::ParallelTracking(TObjArray *const arr, Int_t rfirst, Int_t rlast)
7095 // try to track in parralel
7097 Int_t nseed=arr->GetEntriesFast();
7098 //prepare seeds for tracking
7099 for (Int_t i=0; i<nseed; i++) {
7100 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i), &t=*pt;
7102 if (!t.IsActive()) continue;
7103 // follow prolongation to the first layer
7104 if ( (fSectors ==fInnerSec) || (t.GetFirstPoint()-fkParam->GetNRowLow()>rfirst+1) )
7105 FollowProlongation(t, rfirst+1);
7110 for (Int_t nr=rfirst; nr>=rlast; nr--){
7111 if (nr<fInnerSec->GetNRows())
7112 fSectors = fInnerSec;
7114 fSectors = fOuterSec;
7115 // make indexes with the cluster tracks for given
7117 // find nearest cluster
7118 for (Int_t i=0; i<nseed; i++) {
7119 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i), &t=*pt;
7121 if (nr==80) pt->UpdateReference();
7122 if (!pt->IsActive()) continue;
7123 // if ( (fSectors ==fOuterSec) && (pt->fFirstPoint-fkParam->GetNRowLow())<nr) continue;
7124 if (pt->GetRelativeSector()>17) {
7127 UpdateClusters(t,nr);
7129 // prolonagate to the nearest cluster - if founded
7130 for (Int_t i=0; i<nseed; i++) {
7131 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
7133 if (!pt->IsActive()) continue;
7134 // if ((fSectors ==fOuterSec) && (pt->fFirstPoint-fkParam->GetNRowLow())<nr) continue;
7135 if (pt->GetRelativeSector()>17) {
7138 FollowToNextCluster(*pt,nr);
7143 void AliTPCtrackerMI::PrepareForBackProlongation(const TObjArray *const arr,Float_t fac) const
7147 // if we use TPC track itself we have to "update" covariance
7149 Int_t nseed= arr->GetEntriesFast();
7150 for (Int_t i=0;i<nseed;i++){
7151 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
7155 //rotate to current local system at first accepted point
7156 Int_t index = pt->GetClusterIndex2(pt->GetFirstPoint());
7157 Int_t sec = (index&0xff000000)>>24;
7159 Float_t angle1 = fInnerSec->GetAlpha()*sec+fInnerSec->GetAlphaShift();
7160 if (angle1>TMath::Pi())
7161 angle1-=2.*TMath::Pi();
7162 Float_t angle2 = pt->GetAlpha();
7164 if (TMath::Abs(angle1-angle2)>0.001){
7165 if (!pt->Rotate(angle1-angle2)) return;
7166 //angle2 = pt->GetAlpha();
7167 //pt->fRelativeSector = pt->GetAlpha()/fInnerSec->GetAlpha();
7168 //if (pt->GetAlpha()<0)
7169 // pt->fRelativeSector+=18;
7170 //sec = pt->fRelativeSector;
7179 void AliTPCtrackerMI::PrepareForProlongation(TObjArray *const arr, Float_t fac) const
7183 // if we use TPC track itself we have to "update" covariance
7185 Int_t nseed= arr->GetEntriesFast();
7186 for (Int_t i=0;i<nseed;i++){
7187 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
7190 pt->SetFirstPoint(pt->GetLastPoint());
7198 Int_t AliTPCtrackerMI::PropagateBack(const TObjArray *const arr)
7201 // make back propagation
7203 Int_t nseed= arr->GetEntriesFast();
7204 for (Int_t i=0;i<nseed;i++){
7205 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
7206 if (pt&& pt->GetKinkIndex(0)<=0) {
7207 //AliTPCseed *pt2 = new AliTPCseed(*pt);
7208 fSectors = fInnerSec;
7209 //FollowBackProlongation(*pt,fInnerSec->GetNRows()-1);
7210 //fSectors = fOuterSec;
7211 FollowBackProlongation(*pt,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1,1);
7212 //if (pt->GetNumberOfClusters()<(pt->fEsd->GetTPCclusters(0)) ){
7213 // Error("PropagateBack","Not prolonged track %d",pt->GetLabel());
7214 // FollowBackProlongation(*pt2,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1);
7217 if (pt&& pt->GetKinkIndex(0)>0) {
7218 AliESDkink * kink = fEvent->GetKink(pt->GetKinkIndex(0)-1);
7219 pt->SetFirstPoint(kink->GetTPCRow0());
7220 fSectors = fInnerSec;
7221 FollowBackProlongation(*pt,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1,1);
7229 Int_t AliTPCtrackerMI::PropagateForward2(const TObjArray *const arr)
7232 // make forward propagation
7234 Int_t nseed= arr->GetEntriesFast();
7236 for (Int_t i=0;i<nseed;i++){
7237 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
7239 FollowProlongation(*pt,0,1,1);
7248 Int_t AliTPCtrackerMI::PropagateForward()
7251 // propagate track forward
7253 Int_t nseed = fSeeds->GetEntriesFast();
7254 for (Int_t i=0;i<nseed;i++){
7255 AliTPCseed *pt = (AliTPCseed*)fSeeds->UncheckedAt(i);
7257 AliTPCseed &t = *pt;
7258 Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
7259 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
7260 if (alpha < 0. ) alpha += 2.*TMath::Pi();
7261 t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN);
7265 fSectors = fOuterSec;
7266 ParallelTracking(fSeeds,fOuterSec->GetNRows()+fInnerSec->GetNRows()-1,fInnerSec->GetNRows());
7267 fSectors = fInnerSec;
7268 ParallelTracking(fSeeds,fInnerSec->GetNRows()-1,0);
7277 Int_t AliTPCtrackerMI::PropagateBack(AliTPCseed *const pt, Int_t row0, Int_t row1)
7280 // make back propagation, in between row0 and row1
7284 fSectors = fInnerSec;
7287 if (row1<fSectors->GetNRows())
7290 r1 = fSectors->GetNRows()-1;
7292 if (row0<fSectors->GetNRows()&& r1>0 )
7293 FollowBackProlongation(*pt,r1);
7294 if (row1<=fSectors->GetNRows())
7297 r1 = row1 - fSectors->GetNRows();
7298 if (r1<=0) return 0;
7299 if (r1>=fOuterSec->GetNRows()) return 0;
7300 fSectors = fOuterSec;
7301 return FollowBackProlongation(*pt,r1);
7309 void AliTPCtrackerMI::GetShape(AliTPCseed * seed, Int_t row)
7311 // gets cluster shape
7313 AliTPCClusterParam * clparam = AliTPCcalibDB::Instance()->GetClusterParam();
7314 Float_t zdrift = TMath::Abs((fkParam->GetZLength(0)-TMath::Abs(seed->GetZ())));
7315 Int_t type = (seed->GetSector() < fkParam->GetNSector()/2) ? 0: (row>126) ? 1:2;
7316 Double_t angulary = seed->GetSnp();
7318 if (TMath::Abs(angulary)>AliTPCReconstructor::GetMaxSnpTracker()) {
7319 angulary = TMath::Sign(AliTPCReconstructor::GetMaxSnpTracker(),angulary);
7322 angulary = angulary*angulary/((1.-angulary)*(1.+angulary));
7323 Double_t angularz = seed->GetTgl()*seed->GetTgl()*(1.+angulary);
7325 Double_t sigmay = clparam->GetRMS0(0,type,zdrift,TMath::Sqrt(TMath::Abs(angulary)));
7326 Double_t sigmaz = clparam->GetRMS0(1,type,zdrift,TMath::Sqrt(TMath::Abs(angularz)));
7327 seed->SetCurrentSigmaY2(sigmay*sigmay);
7328 seed->SetCurrentSigmaZ2(sigmaz*sigmaz);
7329 // Float_t sd2 = TMath::Abs((fkParam->GetZLength(0)-TMath::Abs(seed->GetZ())))*fkParam->GetDiffL()*fkParam->GetDiffL();
7330 // // Float_t padlength = fkParam->GetPadPitchLength(seed->fSector);
7331 // Float_t padlength = GetPadPitchLength(row);
7333 // Float_t sresy = (seed->GetSector() < fkParam->GetNSector()/2) ? 0.2 :0.3;
7334 // seed->SetCurrentSigmaY2(sd2+padlength*padlength*angulary/12.+sresy*sresy);
7336 // Float_t sresz = fkParam->GetZSigma();
7337 // seed->SetCurrentSigmaZ2(sd2+padlength*padlength*angularz*angularz*(1+angulary)/12.+sresz*sresz);
7339 Float_t wy = GetSigmaY(seed);
7340 Float_t wz = GetSigmaZ(seed);
7343 if (TMath::Abs(wy/seed->fCurrentSigmaY2-1)>0.0001 || TMath::Abs(wz/seed->fCurrentSigmaZ2-1)>0.0001 ){
7344 printf("problem\n");
7351 //__________________________________________________________________________
7352 void AliTPCtrackerMI::CookLabel(AliKalmanTrack *tk, Float_t wrong) const {
7353 //--------------------------------------------------------------------
7354 //This function "cooks" a track label. If label<0, this track is fake.
7355 //--------------------------------------------------------------------
7356 AliTPCseed * t = dynamic_cast<AliTPCseed*>(tk);
7358 printf("%s:%d wrong type \n",(char*)__FILE__,__LINE__);
7362 Int_t noc=t->GetNumberOfClusters();
7364 //printf("\nnot founded prolongation\n\n\n");
7370 AliTPCclusterMI *clusters[160];
7372 for (Int_t i=0;i<160;i++) {
7379 for (i=0; i<160 && current<noc; i++) {
7381 Int_t index=t->GetClusterIndex2(i);
7382 if (index<=0) continue;
7383 if (index&0x8000) continue;
7385 //clusters[current]=GetClusterMI(index);
7386 if (t->GetClusterPointer(i)){
7387 clusters[current]=t->GetClusterPointer(i);
7393 Int_t lab=123456789;
7394 for (i=0; i<noc; i++) {
7395 AliTPCclusterMI *c=clusters[i];
7397 lab=TMath::Abs(c->GetLabel(0));
7399 for (j=0; j<noc; j++) if (lb[j]==lab || mx[j]==0) break;
7405 for (i=0; i<noc; i++) if (mx[i]>max) {max=mx[i]; lab=lb[i];}
7407 for (i=0; i<noc; i++) {
7408 AliTPCclusterMI *c=clusters[i];
7410 if (TMath::Abs(c->GetLabel(1)) == lab ||
7411 TMath::Abs(c->GetLabel(2)) == lab ) max++;
7413 if (noc<=0) { lab=-1; return;}
7414 if ((1.- Float_t(max)/(noc)) > wrong) lab=-lab;
7417 Int_t tail=Int_t(0.10*noc);
7420 for (i=1; i<160&&ind<tail; i++) {
7421 // AliTPCclusterMI *c=clusters[noc-i];
7422 AliTPCclusterMI *c=clusters[i];
7424 if (lab == TMath::Abs(c->GetLabel(0)) ||
7425 lab == TMath::Abs(c->GetLabel(1)) ||
7426 lab == TMath::Abs(c->GetLabel(2))) max++;
7429 if (max < Int_t(0.5*tail)) lab=-lab;
7436 //delete[] clusters;
7440 //__________________________________________________________________________
7441 Int_t AliTPCtrackerMI::CookLabel(AliTPCseed *const t, Float_t wrong,Int_t first, Int_t last) const {
7442 //--------------------------------------------------------------------
7443 //This function "cooks" a track label. If label<0, this track is fake.
7444 //--------------------------------------------------------------------
7445 Int_t noc=t->GetNumberOfClusters();
7447 //printf("\nnot founded prolongation\n\n\n");
7453 AliTPCclusterMI *clusters[160];
7455 for (Int_t i=0;i<160;i++) {
7462 for (i=0; i<160 && current<noc; i++) {
7463 if (i<first) continue;
7464 if (i>last) continue;
7465 Int_t index=t->GetClusterIndex2(i);
7466 if (index<=0) continue;
7467 if (index&0x8000) continue;
7469 //clusters[current]=GetClusterMI(index);
7470 if (t->GetClusterPointer(i)){
7471 clusters[current]=t->GetClusterPointer(i);
7476 //if (noc<5) return -1;
7477 Int_t lab=123456789;
7478 for (i=0; i<noc; i++) {
7479 AliTPCclusterMI *c=clusters[i];
7481 lab=TMath::Abs(c->GetLabel(0));
7483 for (j=0; j<noc; j++) if (lb[j]==lab || mx[j]==0) break;
7489 for (i=0; i<noc; i++) if (mx[i]>max) {max=mx[i]; lab=lb[i];}
7491 for (i=0; i<noc; i++) {
7492 AliTPCclusterMI *c=clusters[i];
7494 if (TMath::Abs(c->GetLabel(1)) == lab ||
7495 TMath::Abs(c->GetLabel(2)) == lab ) max++;
7497 if (noc<=0) { lab=-1; return -1;}
7498 if ((1.- Float_t(max)/(noc)) > wrong) lab=-lab;
7501 Int_t tail=Int_t(0.10*noc);
7504 for (i=1; i<160&&ind<tail; i++) {
7505 // AliTPCclusterMI *c=clusters[noc-i];
7506 AliTPCclusterMI *c=clusters[i];
7508 if (lab == TMath::Abs(c->GetLabel(0)) ||
7509 lab == TMath::Abs(c->GetLabel(1)) ||
7510 lab == TMath::Abs(c->GetLabel(2))) max++;
7513 if (max < Int_t(0.5*tail)) lab=-lab;
7516 // t->SetLabel(lab);
7520 //delete[] clusters;
7524 Int_t AliTPCtrackerMI::GetRowNumber(Double_t x[3]) const
7526 //return pad row number for given x vector
7527 Float_t phi = TMath::ATan2(x[1],x[0]);
7528 if(phi<0) phi=2.*TMath::Pi()+phi;
7529 // Get the local angle in the sector philoc
7530 const Float_t kRaddeg = 180/3.14159265358979312;
7531 Float_t phiangle = (Int_t (phi*kRaddeg/20.) + 0.5)*20./kRaddeg;
7532 Double_t localx = x[0]*TMath::Cos(phiangle)-x[1]*TMath::Sin(phiangle);
7533 return GetRowNumber(localx);
7538 void AliTPCtrackerMI::MakeESDBitmaps(AliTPCseed *t, AliESDtrack *esd)
7540 //-----------------------------------------------------------------------
7541 // Fill the cluster and sharing bitmaps of the track
7542 //-----------------------------------------------------------------------
7544 Int_t firstpoint = 0;
7545 Int_t lastpoint = 159;
7546 AliTPCTrackerPoint *point;
7547 AliTPCclusterMI *cluster;
7550 TBits clusterMap(159);
7551 TBits sharedMap(159);
7553 for (int iter=firstpoint; iter<lastpoint; iter++) {
7554 // Change to cluster pointers to see if we have a cluster at given padrow
7556 cluster = t->GetClusterPointer(iter);
7558 clusterMap.SetBitNumber(iter, kTRUE);
7559 point = t->GetTrackPoint(iter);
7560 if (point->IsShared())
7561 sharedMap.SetBitNumber(iter,kTRUE);
7563 if (t->GetClusterIndex(iter) >= 0 && (t->GetClusterIndex(iter) & 0x8000) == 0) {
7564 fitMap.SetBitNumber(iter, kTRUE);
7568 esd->SetTPCClusterMap(clusterMap);
7569 esd->SetTPCSharedMap(sharedMap);
7570 esd->SetTPCFitMap(fitMap);
7571 if (nclsf != t->GetNumberOfClusters())
7572 AliWarning(Form("Inconsistency between ncls %d and indices %d (found %d)",t->GetNumberOfClusters(),nclsf,esd->GetTPCClusterMap().CountBits()));
7575 Bool_t AliTPCtrackerMI::IsFindable(AliTPCseed & track){
7577 // return flag if there is findable cluster at given position
7580 Float_t z = track.GetZ();
7582 if (TMath::Abs(z)<(AliTPCReconstructor::GetCtgRange()*track.GetX()+kDeltaZ) &&
7583 TMath::Abs(z)<fkParam->GetZLength(0) &&
7584 (TMath::Abs(track.GetSnp())<AliTPCReconstructor::GetMaxSnpTracker()))
7590 void AliTPCtrackerMI::AddCovariance(AliTPCseed * seed){
7592 // Adding systematic error estimate to the covariance matrix
7593 // !!!! the systematic error for element 4 is in 1/cm not in pt
7595 const Double_t *param = AliTPCReconstructor::GetRecoParam()->GetSystematicError();
7597 // use only the diagonal part if not specified otherwise
7598 if (!AliTPCReconstructor::GetRecoParam()->GetUseSystematicCorrelation()) return AddCovarianceAdd(seed);
7600 Double_t *covarS= (Double_t*)seed->GetCovariance();
7601 Double_t factor[5]={1,1,1,1,1};
7602 Double_t facC = AliTracker::GetBz()*kB2C;
7603 factor[0]= TMath::Sqrt(TMath::Abs((covarS[0] + param[0]*param[0])/covarS[0]));
7604 factor[1]= TMath::Sqrt(TMath::Abs((covarS[2] + param[1]*param[1])/covarS[2]));
7605 factor[2]= TMath::Sqrt(TMath::Abs((covarS[5] + param[2]*param[2])/covarS[5]));
7606 factor[3]= TMath::Sqrt(TMath::Abs((covarS[9] + param[3]*param[3])/covarS[9]));
7607 factor[4]= TMath::Sqrt(TMath::Abs((covarS[14] + facC*facC*param[4]*param[4])/covarS[14]));
7613 for (Int_t i=0; i<5; i++){
7614 for (Int_t j=i; j<5; j++){
7615 Int_t index=seed->GetIndex(i,j);
7616 covarS[index]*=factor[i]*factor[j];
7622 void AliTPCtrackerMI::AddCovarianceAdd(AliTPCseed * seed){
7624 // Adding systematic error - as additive factor without correlation
7626 const Double_t *param = AliTPCReconstructor::GetRecoParam()->GetSystematicError();
7627 Double_t *covarIn= (Double_t*)seed->GetCovariance();
7629 for (Int_t i=0;i<15;i++) covar[i]=0;
7635 covar[0] = param[0]*param[0];
7636 covar[2] = param[1]*param[1];
7637 covar[5] = param[2]*param[2];
7638 covar[9] = param[3]*param[3];
7639 Double_t facC = AliTracker::GetBz()*kB2C;
7640 covar[14]= param[4]*param[4]*facC*facC;
7642 covar[1]=TMath::Sqrt((covar[0]*covar[2]))*covarIn[1]/TMath::Sqrt((covarIn[0]*covarIn[2]));
7644 covar[3]=TMath::Sqrt((covar[0]*covar[5]))*covarIn[3]/TMath::Sqrt((covarIn[0]*covarIn[5]));
7645 covar[4]=TMath::Sqrt((covar[2]*covar[5]))*covarIn[4]/TMath::Sqrt((covarIn[2]*covarIn[5]));
7647 covar[6]=TMath::Sqrt((covar[0]*covar[9]))*covarIn[6]/TMath::Sqrt((covarIn[0]*covarIn[9]));
7648 covar[7]=TMath::Sqrt((covar[2]*covar[9]))*covarIn[7]/TMath::Sqrt((covarIn[2]*covarIn[9]));
7649 covar[8]=TMath::Sqrt((covar[5]*covar[9]))*covarIn[8]/TMath::Sqrt((covarIn[5]*covarIn[9]));
7651 covar[10]=TMath::Sqrt((covar[0]*covar[14]))*covarIn[10]/TMath::Sqrt((covarIn[0]*covarIn[14]));
7652 covar[11]=TMath::Sqrt((covar[2]*covar[14]))*covarIn[11]/TMath::Sqrt((covarIn[2]*covarIn[14]));
7653 covar[12]=TMath::Sqrt((covar[5]*covar[14]))*covarIn[12]/TMath::Sqrt((covarIn[5]*covarIn[14]));
7654 covar[13]=TMath::Sqrt((covar[9]*covar[14]))*covarIn[13]/TMath::Sqrt((covarIn[9]*covarIn[14]));
7656 seed->AddCovariance(covar);
7659 //________________________________________
7660 void AliTPCtrackerMI::MarkSeedFree(TObject *sd)
7662 // account that this seed is "deleted"
7663 AliTPCseed* seed = dynamic_cast<AliTPCseed*>(sd);
7664 if (!sd) {AliError(Form("Freeing of non-AliTPCseed %p from the pool is requested",sd)); return;}
7665 int id = seed->GetPoolID();
7666 if (id<0) {AliError(Form("Freeing of seed %p NOT from the pool is requested",sd)); return;}
7667 // AliInfo(Form("%d %p",id, seed));
7668 fSeedsPool->RemoveAt(id);
7669 if (fFreeSeedsID.GetSize()<=fNFreeSeeds) fFreeSeedsID.Set( 2*fNFreeSeeds + 100 );
7670 fFreeSeedsID.GetArray()[fNFreeSeeds++] = id;
7673 //________________________________________
7674 TObject *&AliTPCtrackerMI::NextFreeSeed()
7676 // return next free slot where the seed can be created
7677 fLastSeedID = fNFreeSeeds ? fFreeSeedsID.GetArray()[--fNFreeSeeds] : fSeedsPool->GetEntriesFast();
7678 // AliInfo(Form("%d",fLastSeedID));
7679 return (*fSeedsPool)[ fLastSeedID ];
7683 //________________________________________
7684 void AliTPCtrackerMI::ResetSeedsPool()
7686 // mark all seeds in the pool as unused
7687 AliInfo(Form("CurrentSize: %d, BookedUpTo: %d, free: %d",fSeedsPool->GetSize(),fSeedsPool->GetEntriesFast(),fNFreeSeeds));
7689 fSeedsPool->Clear("C"); // RS: nominally the seeds may allocate memory...