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) {
383 if(!AliTPCReconstructor::GetRecoParam()->GetUseHLTOnePadCluster()) {
384 if (TMath::Abs(cluster->GetSigmaY2()) < kAlmost0)
396 //_____________________________________________________________________________
397 AliTPCtrackerMI::AliTPCtrackerMI(const AliTPCParam *par):
399 fkNIS(par->GetNInnerSector()/2),
401 fkNOS(par->GetNOuterSector()/2),
423 //---------------------------------------------------------------------
424 // The main TPC tracker constructor
425 //---------------------------------------------------------------------
426 fInnerSec=new AliTPCtrackerSector[fkNIS];
427 fOuterSec=new AliTPCtrackerSector[fkNOS];
430 for (i=0; i<fkNIS; i++) fInnerSec[i].Setup(par,0);
431 for (i=0; i<fkNOS; i++) fOuterSec[i].Setup(par,1);
434 Int_t nrowlow = par->GetNRowLow();
435 Int_t nrowup = par->GetNRowUp();
438 for (i=0;i<nrowlow;i++){
439 fXRow[i] = par->GetPadRowRadiiLow(i);
440 fPadLength[i]= par->GetPadPitchLength(0,i);
441 fYMax[i] = fXRow[i]*TMath::Tan(0.5*par->GetInnerAngle());
445 for (i=0;i<nrowup;i++){
446 fXRow[i+nrowlow] = par->GetPadRowRadiiUp(i);
447 fPadLength[i+nrowlow] = par->GetPadPitchLength(60,i);
448 fYMax[i+nrowlow] = fXRow[i+nrowlow]*TMath::Tan(0.5*par->GetOuterAngle());
451 if (AliTPCReconstructor::StreamLevel()>0) {
452 fDebugStreamer = new TTreeSRedirector("TPCdebug.root");
455 fSeedsPool = new TClonesArray("AliTPCseed",1000);
457 //________________________________________________________________________
458 AliTPCtrackerMI::AliTPCtrackerMI(const AliTPCtrackerMI &t):
484 //------------------------------------
485 // dummy copy constructor
486 //------------------------------------------------------------------
488 for (Int_t irow=0; irow<200; irow++){
495 AliTPCtrackerMI & AliTPCtrackerMI::operator=(const AliTPCtrackerMI& /*r*/)
497 //------------------------------
499 //--------------------------------------------------------------
502 //_____________________________________________________________________________
503 AliTPCtrackerMI::~AliTPCtrackerMI() {
504 //------------------------------------------------------------------
505 // TPC tracker destructor
506 //------------------------------------------------------------------
513 if (fDebugStreamer) delete fDebugStreamer;
514 if (fSeedsPool) delete fSeedsPool;
518 void AliTPCtrackerMI::FillESD(const TObjArray* arr)
522 //fill esds using updated tracks
525 // write tracks to the event
526 // store index of the track
527 Int_t nseed=arr->GetEntriesFast();
528 //FindKinks(arr,fEvent);
529 for (Int_t i=0; i<nseed; i++) {
530 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
534 if (AliTPCReconstructor::StreamLevel()>1) {
535 (*fDebugStreamer)<<"Track0"<<
539 // pt->PropagateTo(fkParam->GetInnerRadiusLow());
540 if (pt->GetKinkIndex(0)<=0){ //don't propagate daughter tracks
541 pt->PropagateTo(fkParam->GetInnerRadiusLow());
544 if (( pt->GetPoints()[2]- pt->GetPoints()[0])>5 && pt->GetPoints()[3]>0.8){
546 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
547 iotrack.SetTPCPoints(pt->GetPoints());
548 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
549 iotrack.SetV0Indexes(pt->GetV0Indexes());
550 // iotrack.SetTPCpid(pt->fTPCr);
551 //iotrack.SetTPCindex(i);
552 fEvent->AddTrack(&iotrack);
556 if ( (pt->GetNumberOfClusters()>70)&& (Float_t(pt->GetNumberOfClusters())/Float_t(pt->GetNFoundable()))>0.55) {
558 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
559 iotrack.SetTPCPoints(pt->GetPoints());
560 //iotrack.SetTPCindex(i);
561 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
562 iotrack.SetV0Indexes(pt->GetV0Indexes());
563 // iotrack.SetTPCpid(pt->fTPCr);
564 fEvent->AddTrack(&iotrack);
568 // short tracks - maybe decays
570 if ( (pt->GetNumberOfClusters()>30) && (Float_t(pt->GetNumberOfClusters())/Float_t(pt->GetNFoundable()))>0.70) {
571 Int_t found,foundable,shared;
572 pt->GetClusterStatistic(0,60,found, foundable,shared,kFALSE);
573 if ( (found>20) && (pt->GetNShared()/float(pt->GetNumberOfClusters())<0.2)){
575 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
576 //iotrack.SetTPCindex(i);
577 iotrack.SetTPCPoints(pt->GetPoints());
578 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
579 iotrack.SetV0Indexes(pt->GetV0Indexes());
580 //iotrack.SetTPCpid(pt->fTPCr);
581 fEvent->AddTrack(&iotrack);
586 if ( (pt->GetNumberOfClusters()>20) && (Float_t(pt->GetNumberOfClusters())/Float_t(pt->GetNFoundable()))>0.8) {
587 Int_t found,foundable,shared;
588 pt->GetClusterStatistic(0,60,found, foundable,shared,kFALSE);
589 if (found<20) continue;
590 if (pt->GetNShared()/float(pt->GetNumberOfClusters())>0.2) continue;
593 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
594 iotrack.SetTPCPoints(pt->GetPoints());
595 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
596 iotrack.SetV0Indexes(pt->GetV0Indexes());
597 //iotrack.SetTPCpid(pt->fTPCr);
598 //iotrack.SetTPCindex(i);
599 fEvent->AddTrack(&iotrack);
602 // short tracks - secondaties
604 if ( (pt->GetNumberOfClusters()>30) ) {
605 Int_t found,foundable,shared;
606 pt->GetClusterStatistic(128,158,found, foundable,shared,kFALSE);
607 if ( (found>20) && (pt->GetNShared()/float(pt->GetNumberOfClusters())<0.2) &&float(found)/float(foundable)>0.8){
609 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
610 iotrack.SetTPCPoints(pt->GetPoints());
611 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
612 iotrack.SetV0Indexes(pt->GetV0Indexes());
613 //iotrack.SetTPCpid(pt->fTPCr);
614 //iotrack.SetTPCindex(i);
615 fEvent->AddTrack(&iotrack);
620 if ( (pt->GetNumberOfClusters()>15)) {
621 Int_t found,foundable,shared;
622 pt->GetClusterStatistic(138,158,found, foundable,shared,kFALSE);
623 if (found<15) continue;
624 if (foundable<=0) continue;
625 if (pt->GetNShared()/float(pt->GetNumberOfClusters())>0.2) continue;
626 if (float(found)/float(foundable)<0.8) continue;
629 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
630 iotrack.SetTPCPoints(pt->GetPoints());
631 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
632 iotrack.SetV0Indexes(pt->GetV0Indexes());
633 // iotrack.SetTPCpid(pt->fTPCr);
634 //iotrack.SetTPCindex(i);
635 fEvent->AddTrack(&iotrack);
639 // >> account for suppressed tracks in the kink indices (RS)
640 int nESDtracks = fEvent->GetNumberOfTracks();
641 for (int it=nESDtracks;it--;) {
642 AliESDtrack* esdTr = fEvent->GetTrack(it);
643 if (!esdTr || !esdTr->GetKinkIndex(0)) continue;
644 for (int ik=0;ik<3;ik++) {
646 if (!(knkId=esdTr->GetKinkIndex(ik))) break; // no more kinks for this track
647 AliESDkink* kink = fEvent->GetKink(TMath::Abs(knkId)-1);
649 AliError(Form("ESDTrack%d refers to non-existing kink %d",it,TMath::Abs(knkId)-1));
652 kink->SetIndex(it, knkId<0 ? 0:1); // update track index of the kink: mother at 0, daughter at 1
655 // << account for suppressed tracks in the kink indices (RS)
656 AliInfo(Form("Number of filled ESDs-\t%d\n",fEvent->GetNumberOfTracks()));
664 Double_t AliTPCtrackerMI::ErrY2(AliTPCseed* seed, const AliTPCclusterMI * cl){
667 // Use calibrated cluster error from OCDB
669 AliTPCClusterParam * clparam = AliTPCcalibDB::Instance()->GetClusterParam();
671 Float_t z = TMath::Abs(fkParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
672 Int_t ctype = cl->GetType();
673 Int_t type = (cl->GetRow()<63) ? 0: (cl->GetRow()>126) ? 1:2;
674 Double_t angle = seed->GetSnp()*seed->GetSnp();
675 angle = TMath::Sqrt(TMath::Abs(angle/(1.-angle)));
676 Double_t erry2 = clparam->GetError0Par(0,type, z,angle);
678 erry2+=0.5; // edge cluster
681 seed->SetErrorY2(erry2);
685 //calculate look-up table at the beginning
686 // static Bool_t ginit = kFALSE;
687 // static Float_t gnoise1,gnoise2,gnoise3;
688 // static Float_t ggg1[10000];
689 // static Float_t ggg2[10000];
690 // static Float_t ggg3[10000];
691 // static Float_t glandau1[10000];
692 // static Float_t glandau2[10000];
693 // static Float_t glandau3[10000];
695 // static Float_t gcor01[500];
696 // static Float_t gcor02[500];
697 // static Float_t gcorp[500];
701 // if (ginit==kFALSE){
702 // for (Int_t i=1;i<500;i++){
703 // Float_t rsigma = float(i)/100.;
704 // gcor02[i] = TMath::Max(0.78 +TMath::Exp(7.4*(rsigma-1.2)),0.6);
705 // gcor01[i] = TMath::Max(0.72 +TMath::Exp(3.36*(rsigma-1.2)),0.6);
706 // gcorp[i] = TMath::Max(TMath::Power((rsigma+0.5),1.5),1.2);
710 // for (Int_t i=3;i<10000;i++){
714 // Float_t amp = float(i);
715 // Float_t padlength =0.75;
716 // gnoise1 = 0.0004/padlength;
717 // Float_t nel = 0.268*amp;
718 // Float_t nprim = 0.155*amp;
719 // ggg1[i] = fkParam->GetDiffT()*fkParam->GetDiffT()*(2+0.001*nel/(padlength*padlength))/nel;
720 // glandau1[i] = (2.+0.12*nprim)*0.5* (2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
721 // if (glandau1[i]>1) glandau1[i]=1;
722 // glandau1[i]*=padlength*padlength/12.;
726 // gnoise2 = 0.0004/padlength;
728 // nprim = 0.133*amp;
729 // ggg2[i] = fkParam->GetDiffT()*fkParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
730 // glandau2[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
731 // if (glandau2[i]>1) glandau2[i]=1;
732 // glandau2[i]*=padlength*padlength/12.;
737 // gnoise3 = 0.0004/padlength;
739 // nprim = 0.133*amp;
740 // ggg3[i] = fkParam->GetDiffT()*fkParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
741 // glandau3[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
742 // if (glandau3[i]>1) glandau3[i]=1;
743 // glandau3[i]*=padlength*padlength/12.;
751 // Int_t amp = int(TMath::Abs(cl->GetQ()));
753 // seed->SetErrorY2(1.);
757 // Float_t z = TMath::Abs(fkParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
758 // Int_t ctype = cl->GetType();
759 // Float_t padlength= GetPadPitchLength(seed->GetRow());
760 // Double_t angle2 = seed->GetSnp()*seed->GetSnp();
761 // angle2 = angle2/(1-angle2);
763 // //cluster "quality"
764 // Int_t rsigmay = int(100.*cl->GetSigmaY2()/(seed->GetCurrentSigmaY2()));
767 // if (fSectors==fInnerSec){
768 // snoise2 = gnoise1;
769 // res = ggg1[amp]*z+glandau1[amp]*angle2;
770 // if (ctype==0) res *= gcor01[rsigmay];
773 // res*= gcorp[rsigmay];
777 // if (padlength<1.1){
778 // snoise2 = gnoise2;
779 // res = ggg2[amp]*z+glandau2[amp]*angle2;
780 // if (ctype==0) res *= gcor02[rsigmay];
783 // res*= gcorp[rsigmay];
787 // snoise2 = gnoise3;
788 // res = ggg3[amp]*z+glandau3[amp]*angle2;
789 // if (ctype==0) res *= gcor02[rsigmay];
792 // res*= gcorp[rsigmay];
799 // res*=2.4; // overestimate error 2 times
803 // if (res<2*snoise2)
806 // seed->SetErrorY2(res);
814 Double_t AliTPCtrackerMI::ErrZ2(AliTPCseed* seed, const AliTPCclusterMI * cl){
817 // Use calibrated cluster error from OCDB
819 AliTPCClusterParam * clparam = AliTPCcalibDB::Instance()->GetClusterParam();
821 Float_t z = TMath::Abs(fkParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
822 Int_t ctype = cl->GetType();
823 Int_t type = (cl->GetRow()<63) ? 0: (cl->GetRow()>126) ? 1:2;
825 Double_t angle2 = seed->GetSnp()*seed->GetSnp();
826 angle2 = seed->GetTgl()*seed->GetTgl()*(1+angle2/(1-angle2));
827 Double_t angle = TMath::Sqrt(TMath::Abs(angle2));
828 Double_t errz2 = clparam->GetError0Par(1,type, z,angle);
830 errz2+=0.5; // edge cluster
833 seed->SetErrorZ2(errz2);
839 // //seed->SetErrorY2(0.1);
841 // //calculate look-up table at the beginning
842 // static Bool_t ginit = kFALSE;
843 // static Float_t gnoise1,gnoise2,gnoise3;
844 // static Float_t ggg1[10000];
845 // static Float_t ggg2[10000];
846 // static Float_t ggg3[10000];
847 // static Float_t glandau1[10000];
848 // static Float_t glandau2[10000];
849 // static Float_t glandau3[10000];
851 // static Float_t gcor01[1000];
852 // static Float_t gcor02[1000];
853 // static Float_t gcorp[1000];
857 // if (ginit==kFALSE){
858 // for (Int_t i=1;i<1000;i++){
859 // Float_t rsigma = float(i)/100.;
860 // gcor02[i] = TMath::Max(0.81 +TMath::Exp(6.8*(rsigma-1.2)),0.6);
861 // gcor01[i] = TMath::Max(0.72 +TMath::Exp(2.04*(rsigma-1.2)),0.6);
862 // gcorp[i] = TMath::Max(TMath::Power((rsigma+0.5),1.5),1.2);
866 // for (Int_t i=3;i<10000;i++){
870 // Float_t amp = float(i);
871 // Float_t padlength =0.75;
872 // gnoise1 = 0.0004/padlength;
873 // Float_t nel = 0.268*amp;
874 // Float_t nprim = 0.155*amp;
875 // ggg1[i] = fkParam->GetDiffT()*fkParam->GetDiffT()*(2+0.001*nel/(padlength*padlength))/nel;
876 // glandau1[i] = (2.+0.12*nprim)*0.5* (2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
877 // if (glandau1[i]>1) glandau1[i]=1;
878 // glandau1[i]*=padlength*padlength/12.;
882 // gnoise2 = 0.0004/padlength;
884 // nprim = 0.133*amp;
885 // ggg2[i] = fkParam->GetDiffT()*fkParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
886 // glandau2[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
887 // if (glandau2[i]>1) glandau2[i]=1;
888 // glandau2[i]*=padlength*padlength/12.;
893 // gnoise3 = 0.0004/padlength;
895 // nprim = 0.133*amp;
896 // ggg3[i] = fkParam->GetDiffT()*fkParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
897 // glandau3[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
898 // if (glandau3[i]>1) glandau3[i]=1;
899 // glandau3[i]*=padlength*padlength/12.;
907 // Int_t amp = int(TMath::Abs(cl->GetQ()));
909 // seed->SetErrorY2(1.);
913 // Float_t z = TMath::Abs(fkParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
914 // Int_t ctype = cl->GetType();
915 // Float_t padlength= GetPadPitchLength(seed->GetRow());
917 // Double_t angle2 = seed->GetSnp()*seed->GetSnp();
918 // // if (angle2<0.6) angle2 = 0.6;
919 // angle2 = seed->GetTgl()*seed->GetTgl()*(1+angle2/(1-angle2));
921 // //cluster "quality"
922 // Int_t rsigmaz = int(100.*cl->GetSigmaZ2()/(seed->GetCurrentSigmaZ2()));
925 // if (fSectors==fInnerSec){
926 // snoise2 = gnoise1;
927 // res = ggg1[amp]*z+glandau1[amp]*angle2;
928 // if (ctype==0) res *= gcor01[rsigmaz];
931 // res*= gcorp[rsigmaz];
935 // if (padlength<1.1){
936 // snoise2 = gnoise2;
937 // res = ggg2[amp]*z+glandau2[amp]*angle2;
938 // if (ctype==0) res *= gcor02[rsigmaz];
941 // res*= gcorp[rsigmaz];
945 // snoise2 = gnoise3;
946 // res = ggg3[amp]*z+glandau3[amp]*angle2;
947 // if (ctype==0) res *= gcor02[rsigmaz];
950 // res*= gcorp[rsigmaz];
959 // if ((ctype<0) &&<70){
964 // if (res<2*snoise2)
966 // if (res>3) res =3;
967 // seed->SetErrorZ2(res);
975 void AliTPCtrackerMI::RotateToLocal(AliTPCseed *seed)
977 //rotate to track "local coordinata
978 Float_t x = seed->GetX();
979 Float_t y = seed->GetY();
980 Float_t ymax = x*TMath::Tan(0.5*fSectors->GetAlpha());
983 seed->SetRelativeSector((seed->GetRelativeSector()+1) % fN);
984 if (!seed->Rotate(fSectors->GetAlpha()))
986 } else if (y <-ymax) {
987 seed->SetRelativeSector((seed->GetRelativeSector()-1+fN) % fN);
988 if (!seed->Rotate(-fSectors->GetAlpha()))
996 //_____________________________________________________________________________
997 Double_t AliTPCtrackerMI::F1old(Double_t x1,Double_t y1,
998 Double_t x2,Double_t y2,
999 Double_t x3,Double_t y3) const
1001 //-----------------------------------------------------------------
1002 // Initial approximation of the track curvature
1003 //-----------------------------------------------------------------
1004 Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
1005 Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
1006 (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
1007 Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
1008 (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
1010 Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
1011 if ( xr*xr+yr*yr<=0.00000000000001) return 100;
1012 return -xr*yr/sqrt(xr*xr+yr*yr);
1017 //_____________________________________________________________________________
1018 Double_t AliTPCtrackerMI::F1(Double_t x1,Double_t y1,
1019 Double_t x2,Double_t y2,
1020 Double_t x3,Double_t y3) const
1022 //-----------------------------------------------------------------
1023 // Initial approximation of the track curvature
1024 //-----------------------------------------------------------------
1030 Double_t det = x3*y2-x2*y3;
1031 if (TMath::Abs(det)<1e-10){
1035 Double_t u = 0.5* (x2*(x2-x3)+y2*(y2-y3))/det;
1036 Double_t x0 = x3*0.5-y3*u;
1037 Double_t y0 = y3*0.5+x3*u;
1038 Double_t c2 = 1/TMath::Sqrt(x0*x0+y0*y0);
1044 Double_t AliTPCtrackerMI::F2(Double_t x1,Double_t y1,
1045 Double_t x2,Double_t y2,
1046 Double_t x3,Double_t y3) const
1048 //-----------------------------------------------------------------
1049 // Initial approximation of the track curvature
1050 //-----------------------------------------------------------------
1056 Double_t det = x3*y2-x2*y3;
1057 if (TMath::Abs(det)<1e-10) {
1061 Double_t u = 0.5* (x2*(x2-x3)+y2*(y2-y3))/det;
1062 Double_t x0 = x3*0.5-y3*u;
1063 Double_t y0 = y3*0.5+x3*u;
1064 Double_t c2 = 1/TMath::Sqrt(x0*x0+y0*y0);
1073 //_____________________________________________________________________________
1074 Double_t AliTPCtrackerMI::F2old(Double_t x1,Double_t y1,
1075 Double_t x2,Double_t y2,
1076 Double_t x3,Double_t y3) const
1078 //-----------------------------------------------------------------
1079 // Initial approximation of the track curvature times center of curvature
1080 //-----------------------------------------------------------------
1081 Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
1082 Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
1083 (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
1084 Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
1085 (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
1087 Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
1089 return -a/(d*y1-b)*xr/sqrt(xr*xr+yr*yr);
1092 //_____________________________________________________________________________
1093 Double_t AliTPCtrackerMI::F3(Double_t x1,Double_t y1,
1094 Double_t x2,Double_t y2,
1095 Double_t z1,Double_t z2) const
1097 //-----------------------------------------------------------------
1098 // Initial approximation of the tangent of the track dip angle
1099 //-----------------------------------------------------------------
1100 return (z1 - z2)/sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
1104 Double_t AliTPCtrackerMI::F3n(Double_t x1,Double_t y1,
1105 Double_t x2,Double_t y2,
1106 Double_t z1,Double_t z2, Double_t c) const
1108 //-----------------------------------------------------------------
1109 // Initial approximation of the tangent of the track dip angle
1110 //-----------------------------------------------------------------
1114 //angle1 = (z1-z2)*c/(TMath::ASin(c*x1-ni)-TMath::ASin(c*x2-ni));
1116 Double_t d = TMath::Sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
1117 if (TMath::Abs(d*c*0.5)>1) return 0;
1118 // Double_t angle2 = TMath::ASin(d*c*0.5);
1119 // Double_t angle2 = AliTPCFastMath::FastAsin(d*c*0.5);
1120 Double_t angle2 = (d*c*0.5>0.1)? TMath::ASin(d*c*0.5): AliTPCFastMath::FastAsin(d*c*0.5);
1122 angle2 = (z1-z2)*c/(angle2*2.);
1126 Bool_t AliTPCtrackerMI::GetProlongation(Double_t x1, Double_t x2, Double_t x[5], Double_t &y, Double_t &z) const
1127 {//-----------------------------------------------------------------
1128 // This function find proloncation of a track to a reference plane x=x2.
1129 //-----------------------------------------------------------------
1133 if (TMath::Abs(x[4]*x1 - x[2]) >= 0.999) {
1137 Double_t c1=x[4]*x1 - x[2], r1=TMath::Sqrt((1.-c1)*(1.+c1));
1138 Double_t c2=x[4]*x2 - x[2], r2=TMath::Sqrt((1.-c2)*(1.+c2));
1142 Double_t dy = dx*(c1+c2)/(r1+r2);
1145 Double_t delta = x[4]*dx*(c1+c2)/(c1*r2 + c2*r1);
1147 if (TMath::Abs(delta)>0.01){
1148 dz = x[3]*TMath::ASin(delta)/x[4];
1150 dz = x[3]*AliTPCFastMath::FastAsin(delta)/x[4];
1153 //dz = x[3]*AliTPCFastMath::FastAsin(delta)/x[4];
1161 Int_t AliTPCtrackerMI::LoadClusters (TTree *const tree)
1166 return LoadClusters();
1170 Int_t AliTPCtrackerMI::LoadClusters(const TObjArray *arr)
1173 // load clusters to the memory
1174 AliTPCClustersRow *clrow = new AliTPCClustersRow("AliTPCclusterMI");
1175 Int_t lower = arr->LowerBound();
1176 Int_t entries = arr->GetEntriesFast();
1178 for (Int_t i=lower; i<entries; i++) {
1179 clrow = (AliTPCClustersRow*) arr->At(i);
1180 if(!clrow) continue;
1181 if(!clrow->GetArray()) continue;
1185 fkParam->AdjustSectorRow(clrow->GetID(),sec,row);
1187 for (Int_t icl=0; icl<clrow->GetArray()->GetEntriesFast(); icl++){
1188 Transform((AliTPCclusterMI*)(clrow->GetArray()->At(icl)));
1191 if (clrow->GetArray()->GetEntriesFast()<=0) continue;
1192 AliTPCtrackerRow * tpcrow=0;
1195 tpcrow = &(fInnerSec[sec%fkNIS][row]);
1199 tpcrow = &(fOuterSec[(sec-fkNIS*2)%fkNOS][row]);
1200 left = (sec-fkNIS*2)/fkNOS;
1203 tpcrow->SetN1(clrow->GetArray()->GetEntriesFast());
1204 tpcrow->SetClusters1(new AliTPCclusterMI[tpcrow->GetN1()]);
1205 for (Int_t j=0;j<tpcrow->GetN1();++j)
1206 tpcrow->SetCluster1(j, *(AliTPCclusterMI*)(clrow->GetArray()->At(j)));
1209 tpcrow->SetN2(clrow->GetArray()->GetEntriesFast());
1210 tpcrow->SetClusters2(new AliTPCclusterMI[tpcrow->GetN2()]);
1211 for (Int_t j=0;j<tpcrow->GetN2();++j)
1212 tpcrow->SetCluster2(j,*(AliTPCclusterMI*)(clrow->GetArray()->At(j)));
1214 clrow->GetArray()->Clear("C");
1223 Int_t AliTPCtrackerMI::LoadClusters(const TClonesArray *arr)
1226 // load clusters to the memory from one
1229 AliTPCclusterMI *clust=0;
1230 Int_t count[72][96] = { {0} , {0} };
1232 // loop over clusters
1233 for (Int_t icl=0; icl<arr->GetEntriesFast(); icl++) {
1234 clust = (AliTPCclusterMI*)arr->At(icl);
1235 if(!clust) continue;
1236 //printf("cluster: det %d, row %d \n", clust->GetDetector(),clust->GetRow());
1238 // transform clusters
1241 // count clusters per pad row
1242 count[clust->GetDetector()][clust->GetRow()]++;
1245 // insert clusters to sectors
1246 for (Int_t icl=0; icl<arr->GetEntriesFast(); icl++) {
1247 clust = (AliTPCclusterMI*)arr->At(icl);
1248 if(!clust) continue;
1250 Int_t sec = clust->GetDetector();
1251 Int_t row = clust->GetRow();
1253 // filter overlapping pad rows needed by HLT
1254 if(sec<fkNIS*2) { //IROCs
1255 if(row == 30) continue;
1258 if(row == 27 || row == 76) continue;
1264 fInnerSec[sec%fkNIS].InsertCluster(clust, count[sec][row], fkParam);
1267 left = (sec-fkNIS*2)/fkNOS;
1268 fOuterSec[(sec-fkNIS*2)%fkNOS].InsertCluster(clust, count[sec][row], fkParam);
1272 // Load functions must be called behind LoadCluster(TClonesArray*)
1274 //LoadOuterSectors();
1275 //LoadInnerSectors();
1281 Int_t AliTPCtrackerMI::LoadClusters()
1284 // load clusters to the memory
1285 AliTPCClustersRow *clrow= new AliTPCClustersRow("AliTPCclusterMI");
1287 // TTree * tree = fClustersArray.GetTree();
1289 TTree * tree = fInput;
1290 TBranch * br = tree->GetBranch("Segment");
1291 br->SetAddress(&clrow);
1293 // Conversion of pad, row coordinates in local tracking coords.
1294 // Could be skipped here; is already done in clusterfinder
1296 Int_t j=Int_t(tree->GetEntries());
1297 for (Int_t i=0; i<j; i++) {
1301 fkParam->AdjustSectorRow(clrow->GetID(),sec,row);
1302 for (Int_t icl=0; icl<clrow->GetArray()->GetEntriesFast(); icl++){
1303 Transform((AliTPCclusterMI*)(clrow->GetArray()->At(icl)));
1306 AliTPCtrackerRow * tpcrow=0;
1309 tpcrow = &(fInnerSec[sec%fkNIS][row]);
1313 tpcrow = &(fOuterSec[(sec-fkNIS*2)%fkNOS][row]);
1314 left = (sec-fkNIS*2)/fkNOS;
1317 tpcrow->SetN1(clrow->GetArray()->GetEntriesFast());
1318 tpcrow->SetClusters1(new AliTPCclusterMI[tpcrow->GetN1()]);
1319 for (Int_t k=0;k<tpcrow->GetN1();++k)
1320 tpcrow->SetCluster1(k, *(AliTPCclusterMI*)(clrow->GetArray()->At(k)));
1323 tpcrow->SetN2(clrow->GetArray()->GetEntriesFast());
1324 tpcrow->SetClusters2(new AliTPCclusterMI[tpcrow->GetN2()]);
1325 for (Int_t k=0;k<tpcrow->GetN2();++k)
1326 tpcrow->SetCluster2(k,*(AliTPCclusterMI*)(clrow->GetArray()->At(k)));
1337 void AliTPCtrackerMI::UnloadClusters()
1340 // unload clusters from the memory
1342 Int_t nrows = fOuterSec->GetNRows();
1343 for (Int_t sec = 0;sec<fkNOS;sec++)
1344 for (Int_t row = 0;row<nrows;row++){
1345 AliTPCtrackerRow* tpcrow = &(fOuterSec[sec%fkNOS][row]);
1347 // if (tpcrow->fClusters1) delete []tpcrow->fClusters1;
1348 // if (tpcrow->fClusters2) delete []tpcrow->fClusters2;
1350 tpcrow->ResetClusters();
1353 nrows = fInnerSec->GetNRows();
1354 for (Int_t sec = 0;sec<fkNIS;sec++)
1355 for (Int_t row = 0;row<nrows;row++){
1356 AliTPCtrackerRow* tpcrow = &(fInnerSec[sec%fkNIS][row]);
1358 // if (tpcrow->fClusters1) delete []tpcrow->fClusters1;
1359 //if (tpcrow->fClusters2) delete []tpcrow->fClusters2;
1361 tpcrow->ResetClusters();
1367 void AliTPCtrackerMI::FillClusterArray(TObjArray* array) const{
1369 // Filling cluster to the array - For visualization purposes
1372 nrows = fOuterSec->GetNRows();
1373 for (Int_t sec = 0;sec<fkNOS;sec++)
1374 for (Int_t row = 0;row<nrows;row++){
1375 AliTPCtrackerRow* tpcrow = &(fOuterSec[sec%fkNOS][row]);
1376 if (!tpcrow) continue;
1377 for (Int_t icl = 0;icl<tpcrow->GetN();icl++){
1378 array->AddLast((TObject*)((*tpcrow)[icl]));
1381 nrows = fInnerSec->GetNRows();
1382 for (Int_t sec = 0;sec<fkNIS;sec++)
1383 for (Int_t row = 0;row<nrows;row++){
1384 AliTPCtrackerRow* tpcrow = &(fInnerSec[sec%fkNIS][row]);
1385 if (!tpcrow) continue;
1386 for (Int_t icl = 0;icl<tpcrow->GetN();icl++){
1387 array->AddLast((TObject*)(*tpcrow)[icl]);
1393 void AliTPCtrackerMI::Transform(AliTPCclusterMI * cluster){
1397 AliTPCcalibDB * calibDB = AliTPCcalibDB::Instance();
1398 AliTPCTransform *transform = calibDB->GetTransform() ;
1400 AliFatal("Tranformations not in calibDB");
1403 transform->SetCurrentRecoParam((AliTPCRecoParam*)AliTPCReconstructor::GetRecoParam());
1404 Double_t x[3]={cluster->GetRow(),cluster->GetPad(),cluster->GetTimeBin()};
1405 Int_t i[1]={cluster->GetDetector()};
1406 transform->Transform(x,i,0,1);
1407 // if (cluster->GetDetector()%36>17){
1412 // in debug mode check the transformation
1414 if (AliTPCReconstructor::StreamLevel()>2) {
1416 cluster->GetGlobalXYZ(gx);
1417 Int_t event = (fEvent==NULL)? 0: fEvent->GetEventNumberInFile();
1418 TTreeSRedirector &cstream = *fDebugStreamer;
1419 cstream<<"Transform"<<
1430 cluster->SetX(x[0]);
1431 cluster->SetY(x[1]);
1432 cluster->SetZ(x[2]);
1437 //if (!fkParam->IsGeoRead()) fkParam->ReadGeoMatrices();
1438 if (AliTPCReconstructor::GetRecoParam()->GetUseSectorAlignment() && (!calibDB->HasAlignmentOCDB())){
1439 TGeoHMatrix *mat = fkParam->GetClusterMatrix(cluster->GetDetector());
1441 Double_t pos[3]= {cluster->GetX(),cluster->GetY(),cluster->GetZ()};
1442 Double_t posC[3]={cluster->GetX(),cluster->GetY(),cluster->GetZ()};
1443 if (mat) mat->LocalToMaster(pos,posC);
1445 // chack Loading of Geo matrices from GeoManager - TEMPORARY FIX
1447 cluster->SetX(posC[0]);
1448 cluster->SetY(posC[1]);
1449 cluster->SetZ(posC[2]);
1453 //_____________________________________________________________________________
1454 Int_t AliTPCtrackerMI::LoadOuterSectors() {
1455 //-----------------------------------------------------------------
1456 // This function fills outer TPC sectors with clusters.
1457 //-----------------------------------------------------------------
1458 Int_t nrows = fOuterSec->GetNRows();
1460 for (Int_t sec = 0;sec<fkNOS;sec++)
1461 for (Int_t row = 0;row<nrows;row++){
1462 AliTPCtrackerRow* tpcrow = &(fOuterSec[sec%fkNOS][row]);
1463 Int_t sec2 = sec+2*fkNIS;
1465 Int_t ncl = tpcrow->GetN1();
1467 AliTPCclusterMI *c= (tpcrow->GetCluster1(ncl));
1468 index=(((sec2<<8)+row)<<16)+ncl;
1469 tpcrow->InsertCluster(c,index);
1472 ncl = tpcrow->GetN2();
1474 AliTPCclusterMI *c= (tpcrow->GetCluster2(ncl));
1475 index=((((sec2+fkNOS)<<8)+row)<<16)+ncl;
1476 tpcrow->InsertCluster(c,index);
1479 // write indexes for fast acces
1481 for (Int_t i=0;i<510;i++)
1482 tpcrow->SetFastCluster(i,-1);
1483 for (Int_t i=0;i<tpcrow->GetN();i++){
1484 Int_t zi = Int_t((*tpcrow)[i]->GetZ()+255.);
1485 tpcrow->SetFastCluster(zi,i); // write index
1488 for (Int_t i=0;i<510;i++){
1489 if (tpcrow->GetFastCluster(i)<0)
1490 tpcrow->SetFastCluster(i,last);
1492 last = tpcrow->GetFastCluster(i);
1501 //_____________________________________________________________________________
1502 Int_t AliTPCtrackerMI::LoadInnerSectors() {
1503 //-----------------------------------------------------------------
1504 // This function fills inner TPC sectors with clusters.
1505 //-----------------------------------------------------------------
1506 Int_t nrows = fInnerSec->GetNRows();
1508 for (Int_t sec = 0;sec<fkNIS;sec++)
1509 for (Int_t row = 0;row<nrows;row++){
1510 AliTPCtrackerRow* tpcrow = &(fInnerSec[sec%fkNIS][row]);
1513 Int_t ncl = tpcrow->GetN1();
1515 AliTPCclusterMI *c= (tpcrow->GetCluster1(ncl));
1516 index=(((sec<<8)+row)<<16)+ncl;
1517 tpcrow->InsertCluster(c,index);
1520 ncl = tpcrow->GetN2();
1522 AliTPCclusterMI *c= (tpcrow->GetCluster2(ncl));
1523 index=((((sec+fkNIS)<<8)+row)<<16)+ncl;
1524 tpcrow->InsertCluster(c,index);
1527 // write indexes for fast acces
1529 for (Int_t i=0;i<510;i++)
1530 tpcrow->SetFastCluster(i,-1);
1531 for (Int_t i=0;i<tpcrow->GetN();i++){
1532 Int_t zi = Int_t((*tpcrow)[i]->GetZ()+255.);
1533 tpcrow->SetFastCluster(zi,i); // write index
1536 for (Int_t i=0;i<510;i++){
1537 if (tpcrow->GetFastCluster(i)<0)
1538 tpcrow->SetFastCluster(i,last);
1540 last = tpcrow->GetFastCluster(i);
1552 //_________________________________________________________________________
1553 AliTPCclusterMI *AliTPCtrackerMI::GetClusterMI(Int_t index) const {
1554 //--------------------------------------------------------------------
1555 // Return pointer to a given cluster
1556 //--------------------------------------------------------------------
1557 if (index<0) return 0; // no cluster
1558 Int_t sec=(index&0xff000000)>>24;
1559 Int_t row=(index&0x00ff0000)>>16;
1560 Int_t ncl=(index&0x00007fff)>>00;
1562 const AliTPCtrackerRow * tpcrow=0;
1563 AliTPCclusterMI * clrow =0;
1565 if (sec<0 || sec>=fkNIS*4) {
1566 AliWarning(Form("Wrong sector %d",sec));
1571 AliTPCtrackerSector& tracksec = fInnerSec[sec%fkNIS];
1572 if (tracksec.GetNRows()<=row) return 0;
1573 tpcrow = &(tracksec[row]);
1574 if (tpcrow==0) return 0;
1577 if (tpcrow->GetN1()<=ncl) return 0;
1578 clrow = tpcrow->GetClusters1();
1581 if (tpcrow->GetN2()<=ncl) return 0;
1582 clrow = tpcrow->GetClusters2();
1586 AliTPCtrackerSector& tracksec = fOuterSec[(sec-fkNIS*2)%fkNOS];
1587 if (tracksec.GetNRows()<=row) return 0;
1588 tpcrow = &(tracksec[row]);
1589 if (tpcrow==0) return 0;
1591 if (sec-2*fkNIS<fkNOS) {
1592 if (tpcrow->GetN1()<=ncl) return 0;
1593 clrow = tpcrow->GetClusters1();
1596 if (tpcrow->GetN2()<=ncl) return 0;
1597 clrow = tpcrow->GetClusters2();
1601 return &(clrow[ncl]);
1607 Int_t AliTPCtrackerMI::FollowToNext(AliTPCseed& t, Int_t nr) {
1608 //-----------------------------------------------------------------
1609 // This function tries to find a track prolongation to next pad row
1610 //-----------------------------------------------------------------
1612 Double_t x= GetXrow(nr), ymax=GetMaxY(nr);
1615 AliTPCclusterMI *cl=0;
1616 Int_t tpcindex= t.GetClusterIndex2(nr);
1618 // update current shape info every 5 pad-row
1619 // if ( (nr%5==0) || t.GetNumberOfClusters()<2 || (t.fCurrentSigmaY2<0.0001) ){
1623 if (fIteration>0 && tpcindex>=-1){ //if we have already clusters
1625 if (tpcindex==-1) return 0; //track in dead zone
1626 if (tpcindex >= 0){ //
1627 cl = t.GetClusterPointer(nr);
1628 //if (cl==0) cl = GetClusterMI(tpcindex);
1629 if (!cl) cl = GetClusterMI(tpcindex);
1630 t.SetCurrentClusterIndex1(tpcindex);
1633 Int_t relativesector = ((tpcindex&0xff000000)>>24)%18; // if previously accepted cluster in different sector
1634 Float_t angle = relativesector*fSectors->GetAlpha()+fSectors->GetAlphaShift();
1636 if (angle<-TMath::Pi()) angle += 2*TMath::Pi();
1637 if (angle>=TMath::Pi()) angle -= 2*TMath::Pi();
1639 if (TMath::Abs(angle-t.GetAlpha())>0.001){
1640 Double_t rotation = angle-t.GetAlpha();
1641 t.SetRelativeSector(relativesector);
1642 if (!t.Rotate(rotation)) return 0;
1644 if (!t.PropagateTo(x)) return 0;
1646 t.SetCurrentCluster(cl);
1648 Int_t accept = AcceptCluster(&t,t.GetCurrentCluster());
1649 if ((tpcindex&0x8000)==0) accept =0;
1651 //if founded cluster is acceptible
1652 if (cl->IsUsed(11)) { // id cluster is shared inrease uncertainty
1653 t.SetErrorY2(t.GetErrorY2()+0.03);
1654 t.SetErrorZ2(t.GetErrorZ2()+0.03);
1655 t.SetErrorY2(t.GetErrorY2()*3);
1656 t.SetErrorZ2(t.GetErrorZ2()*3);
1658 t.SetNFoundable(t.GetNFoundable()+1);
1659 UpdateTrack(&t,accept);
1662 else { // Remove old cluster from track
1663 t.SetClusterIndex(nr, -3);
1664 t.SetClusterPointer(nr, 0);
1668 if (TMath::Abs(t.GetSnp())>AliTPCReconstructor::GetMaxSnpTracker()) return 0; // cut on angle
1669 if (fIteration>1 && IsFindable(t)){
1670 // not look for new cluster during refitting
1671 t.SetNFoundable(t.GetNFoundable()+1);
1676 // if (TMath::Abs(t.GetSnp())>0.95 || TMath::Abs(x*t.GetC()-t.GetEta())>0.95) return 0;// patch 28 fev 06
1677 if (!t.PropagateTo(x)) {
1678 if (fIteration==0) t.SetRemoval(10);
1681 Double_t y = t.GetY();
1682 if (TMath::Abs(y)>ymax){
1684 t.SetRelativeSector((t.GetRelativeSector()+1) % fN);
1685 if (!t.Rotate(fSectors->GetAlpha()))
1687 } else if (y <-ymax) {
1688 t.SetRelativeSector((t.GetRelativeSector()-1+fN) % fN);
1689 if (!t.Rotate(-fSectors->GetAlpha()))
1692 if (!t.PropagateTo(x)) {
1693 if (fIteration==0) t.SetRemoval(10);
1699 Double_t z=t.GetZ();
1702 if (!IsActive(t.GetRelativeSector(),nr)) {
1704 t.SetClusterIndex2(nr,-1);
1707 //AliInfo(Form("A - Sector%d phi %f - alpha %f", t.fRelativeSector,y/x, t.GetAlpha()));
1708 Bool_t isActive = IsActive(t.GetRelativeSector(),nr);
1709 Bool_t isActive2 = (nr>=fInnerSec->GetNRows()) ? fOuterSec[t.GetRelativeSector()][nr-fInnerSec->GetNRows()].GetN()>0:fInnerSec[t.GetRelativeSector()][nr].GetN()>0;
1711 if (!isActive || !isActive2) return 0;
1713 const AliTPCtrackerRow &krow=GetRow(t.GetRelativeSector(),nr);
1714 if ( (t.GetSigmaY2()<0) || t.GetSigmaZ2()<0) return 0;
1716 Double_t roadz = 1.;
1718 if (TMath::Abs(TMath::Abs(y)-ymax)<krow.GetDeadZone()){
1720 t.SetClusterIndex2(nr,-1);
1726 // if (TMath::Abs(z)<(AliTPCReconstructor::GetCtgRange()*x+10) && TMath::Abs(z)<fkParam->GetZLength(0) && (TMath::Abs(t.GetSnp())<AliTPCReconstructor::GetMaxSnpTracker()))
1727 t.SetNFoundable(t.GetNFoundable()+1);
1733 // cl = krow.FindNearest2(y+10.,z,roady,roadz,index);
1734 cl = krow.FindNearest2(y,z,roady,roadz,index);
1735 if (cl) t.SetCurrentClusterIndex1(krow.GetIndex(index));
1738 t.SetCurrentCluster(cl);
1740 if (fIteration==2&&cl->IsUsed(10)) return 0;
1741 Int_t accept = AcceptCluster(&t,t.GetCurrentCluster());
1742 if (fIteration==2&&cl->IsUsed(11)) {
1743 t.SetErrorY2(t.GetErrorY2()+0.03);
1744 t.SetErrorZ2(t.GetErrorZ2()+0.03);
1745 t.SetErrorY2(t.GetErrorY2()*3);
1746 t.SetErrorZ2(t.GetErrorZ2()*3);
1749 if (t.fCurrentCluster->IsUsed(10)){
1754 if (t.fNShared>0.7*t.GetNumberOfClusters()) {
1760 if (accept<3) UpdateTrack(&t,accept);
1763 if ( fIteration==0 && t.GetNFoundable()*0.5 > t.GetNumberOfClusters()) t.SetRemoval(10);
1771 //_________________________________________________________________________
1772 Bool_t AliTPCtrackerMI::GetTrackPoint(Int_t index, AliTrackPoint &p ) const
1774 // Get track space point by index
1775 // return false in case the cluster doesn't exist
1776 AliTPCclusterMI *cl = GetClusterMI(index);
1777 if (!cl) return kFALSE;
1778 Int_t sector = (index&0xff000000)>>24;
1779 // Int_t row = (index&0x00ff0000)>>16;
1781 // xyz[0] = fkParam->GetPadRowRadii(sector,row);
1782 xyz[0] = cl->GetX();
1783 xyz[1] = cl->GetY();
1784 xyz[2] = cl->GetZ();
1786 fkParam->AdjustCosSin(sector,cos,sin);
1787 Float_t x = cos*xyz[0]-sin*xyz[1];
1788 Float_t y = cos*xyz[1]+sin*xyz[0];
1790 Float_t sigmaY2 = 0.027*cl->GetSigmaY2();
1791 if (sector < fkParam->GetNInnerSector()) sigmaY2 *= 2.07;
1792 Float_t sigmaZ2 = 0.066*cl->GetSigmaZ2();
1793 if (sector < fkParam->GetNInnerSector()) sigmaZ2 *= 1.77;
1794 cov[0] = sin*sin*sigmaY2;
1795 cov[1] = -sin*cos*sigmaY2;
1797 cov[3] = cos*cos*sigmaY2;
1800 p.SetXYZ(x,y,xyz[2],cov);
1801 AliGeomManager::ELayerID iLayer;
1803 if (sector < fkParam->GetNInnerSector()) {
1804 iLayer = AliGeomManager::kTPC1;
1808 iLayer = AliGeomManager::kTPC2;
1809 idet = sector - fkParam->GetNInnerSector();
1811 UShort_t volid = AliGeomManager::LayerToVolUID(iLayer,idet);
1812 p.SetVolumeID(volid);
1818 Int_t AliTPCtrackerMI::UpdateClusters(AliTPCseed& t, Int_t nr) {
1819 //-----------------------------------------------------------------
1820 // This function tries to find a track prolongation to next pad row
1821 //-----------------------------------------------------------------
1822 t.SetCurrentCluster(0);
1823 t.SetCurrentClusterIndex1(-3);
1825 Double_t xt=t.GetX();
1826 Int_t row = GetRowNumber(xt)-1;
1827 Double_t ymax= GetMaxY(nr);
1829 if (row < nr) return 1; // don't prolongate if not information until now -
1830 // if (TMath::Abs(t.GetSnp())>0.9 && t.GetNumberOfClusters()>40. && fIteration!=2) {
1832 // return 0; // not prolongate strongly inclined tracks
1834 // if (TMath::Abs(t.GetSnp())>0.95) {
1836 // return 0; // not prolongate strongly inclined tracks
1837 // }// patch 28 fev 06
1839 Double_t x= GetXrow(nr);
1841 //t.PropagateTo(x+0.02);
1842 //t.PropagateTo(x+0.01);
1843 if (!t.PropagateTo(x)){
1850 if (TMath::Abs(y)>ymax){
1852 t.SetRelativeSector((t.GetRelativeSector()+1) % fN);
1853 if (!t.Rotate(fSectors->GetAlpha()))
1855 } else if (y <-ymax) {
1856 t.SetRelativeSector((t.GetRelativeSector()-1+fN) % fN);
1857 if (!t.Rotate(-fSectors->GetAlpha()))
1860 // if (!t.PropagateTo(x)){
1867 if (TMath::Abs(t.GetSnp())>AliTPCReconstructor::GetMaxSnpTracker()) return 0;
1869 if (!IsActive(t.GetRelativeSector(),nr)) {
1871 t.SetClusterIndex2(nr,-1);
1874 //AliInfo(Form("A - Sector%d phi %f - alpha %f", t.fRelativeSector,y/x, t.GetAlpha()));
1876 AliTPCtrackerRow &krow=GetRow(t.GetRelativeSector(),nr);
1878 if (TMath::Abs(TMath::Abs(y)-ymax)<krow.GetDeadZone()){
1880 t.SetClusterIndex2(nr,-1);
1886 // if (TMath::Abs(t.GetZ())<(AliTPCReconstructor::GetCtgRange()*t.GetX()+10) && (TMath::Abs(t.GetSnp())<AliTPCReconstructor::GetMaxSnpTracker()))
1887 if (IsFindable(t)) t.SetNFoundable(t.GetNFoundable()+1);
1893 if ( (nr%2==0) || t.GetNumberOfClusters()<2){
1894 // t.fCurrentSigmaY = GetSigmaY(&t);
1895 //t.fCurrentSigmaZ = GetSigmaZ(&t);
1899 AliTPCclusterMI *cl=0;
1902 Double_t roady = 1.;
1903 Double_t roadz = 1.;
1907 index = t.GetClusterIndex2(nr);
1908 if ( (index >= 0) && (index&0x8000)==0){
1909 cl = t.GetClusterPointer(nr);
1910 if ( (cl==0) && (index >= 0)) cl = GetClusterMI(index);
1911 t.SetCurrentClusterIndex1(index);
1913 t.SetCurrentCluster(cl);
1919 // if (index<0) return 0;
1920 UInt_t uindex = TMath::Abs(index);
1923 //cl = krow.FindNearest2(y+10,z,roady,roadz,uindex);
1924 cl = krow.FindNearest2(y,z,roady,roadz,uindex);
1927 if (cl) t.SetCurrentClusterIndex1(krow.GetIndex(uindex));
1928 t.SetCurrentCluster(cl);
1934 Int_t AliTPCtrackerMI::FollowToNextCluster(AliTPCseed & t, Int_t nr) {
1935 //-----------------------------------------------------------------
1936 // This function tries to find a track prolongation to next pad row
1937 //-----------------------------------------------------------------
1939 //update error according neighborhoud
1941 if (t.GetCurrentCluster()) {
1943 Int_t accept = AcceptCluster(&t,t.GetCurrentCluster());
1945 if (t.GetCurrentCluster()->IsUsed(10)){
1950 t.SetNShared(t.GetNShared()+1);
1951 if (t.GetNShared()>0.7*t.GetNumberOfClusters()) {
1956 if (fIteration>0) accept = 0;
1957 if (accept<3) UpdateTrack(&t,accept);
1961 if ( t.GetNumberOfClusters()>18 && ( (t.GetSigmaY2()+t.GetSigmaZ2())>0.16)) t.SetRemoval(10);
1962 if ( t.GetNumberOfClusters()>18 && t.GetChi2()/t.GetNumberOfClusters()>6 ) t.SetRemoval(10);
1964 if (( (t.GetNFoundable()*0.5 > t.GetNumberOfClusters()) || t.GetNoCluster()>15)) t.SetRemoval(10);
1972 //_____________________________________________________________________________
1973 Int_t AliTPCtrackerMI::FollowProlongation(AliTPCseed& t, Int_t rf, Int_t step, Bool_t fromSeeds) {
1974 //-----------------------------------------------------------------
1975 // This function tries to find a track prolongation.
1976 //-----------------------------------------------------------------
1977 Double_t xt=t.GetX();
1979 Double_t alpha=t.GetAlpha();
1980 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
1981 if (alpha < 0. ) alpha += 2.*TMath::Pi();
1983 t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN);
1985 Int_t first = GetRowNumber(xt);
1990 for (Int_t nr= first; nr>=rf; nr-=step) {
1992 if (t.GetKinkIndexes()[0]>0){
1993 for (Int_t i=0;i<3;i++){
1994 Int_t index = t.GetKinkIndexes()[i];
1995 if (index==0) break;
1996 if (index<0) continue;
1998 AliKink * kink = (AliKink*)fEvent->GetKink(index-1);
2000 printf("PROBLEM\n");
2003 Int_t kinkrow = kink->GetTPCRow0()+2+Int_t(0.5/(0.05+kink->GetAngle(2)));
2005 AliExternalTrackParam paramd(t);
2006 kink->SetDaughter(paramd);
2007 kink->SetStatus(2,5);
2014 if (nr==80) t.UpdateReference();
2015 if (nr<fInnerSec->GetNRows())
2016 fSectors = fInnerSec;
2018 fSectors = fOuterSec;
2019 if (FollowToNext(t,nr)==0)
2032 Int_t AliTPCtrackerMI::FollowBackProlongation(AliTPCseed& t, Int_t rf, Bool_t fromSeeds) {
2033 //-----------------------------------------------------------------
2034 // This function tries to find a track prolongation.
2035 //-----------------------------------------------------------------
2037 Double_t xt=t.GetX();
2038 Double_t alpha=t.GetAlpha();
2039 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
2040 if (alpha < 0. ) alpha += 2.*TMath::Pi();
2041 t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN);
2043 Int_t first = t.GetFirstPoint();
2044 Int_t ri = GetRowNumber(xt);
2048 if (first<ri) first = ri;
2050 if (first<0) first=0;
2051 for (Int_t nr=first; nr<=rf; nr++) {
2052 // if ( (TMath::Abs(t.GetSnp())>0.95)) break;//patch 28 fev 06
2053 if (t.GetKinkIndexes()[0]<0){
2054 for (Int_t i=0;i<3;i++){
2055 Int_t index = t.GetKinkIndexes()[i];
2056 if (index==0) break;
2057 if (index>0) continue;
2058 index = TMath::Abs(index);
2059 AliKink * kink = (AliKink*)fEvent->GetKink(index-1);
2061 printf("PROBLEM\n");
2064 Int_t kinkrow = kink->GetTPCRow0()-2-Int_t(0.5/(0.05+kink->GetAngle(2)));
2066 AliExternalTrackParam paramm(t);
2067 kink->SetMother(paramm);
2068 kink->SetStatus(2,1);
2075 if (nr<fInnerSec->GetNRows())
2076 fSectors = fInnerSec;
2078 fSectors = fOuterSec;
2089 Float_t AliTPCtrackerMI::OverlapFactor(AliTPCseed * s1, AliTPCseed * s2, Int_t &sum1, Int_t & sum2)
2091 // overlapping factor
2097 Float_t dz2 =(s1->GetZ() - s2->GetZ());
2100 Float_t dy2 =TMath::Abs((s1->GetY() - s2->GetY()));
2102 Float_t distance = TMath::Sqrt(dz2+dy2);
2103 if (distance>4.) return 0; // if there are far away - not overlap - to reduce combinatorics
2106 Int_t firstpoint = TMath::Min(s1->GetFirstPoint(),s2->GetFirstPoint());
2107 Int_t lastpoint = TMath::Max(s1->GetLastPoint(),s2->GetLastPoint());
2112 if (firstpoint>lastpoint) {
2113 firstpoint =lastpoint;
2118 for (Int_t i=firstpoint-1;i<lastpoint+1;i++){
2119 if (s1->GetClusterIndex2(i)>0) sum1++;
2120 if (s2->GetClusterIndex2(i)>0) sum2++;
2121 if (s1->GetClusterIndex2(i)==s2->GetClusterIndex2(i) && s1->GetClusterIndex2(i)>0) {
2125 if (sum<5) return 0;
2127 Float_t summin = TMath::Min(sum1+1,sum2+1);
2128 Float_t ratio = (sum+1)/Float_t(summin);
2132 void AliTPCtrackerMI::SignShared(AliTPCseed * s1, AliTPCseed * s2)
2136 Float_t thetaCut = 0.2;//+10.*TMath::Sqrt(s1->GetSigmaTglZ()+ s2->GetSigmaTglZ());
2137 if (TMath::Abs(s1->GetTgl()-s2->GetTgl())>thetaCut) return;
2138 Float_t minCl = TMath::Min(s1->GetNumberOfClusters(),s2->GetNumberOfClusters());
2139 Int_t cutN0 = TMath::Max(5,TMath::Nint(0.1*minCl));
2144 //Int_t firstpoint = TMath::Max(s1->GetFirstPoint(),s2->GetFirstPoint());
2145 //Int_t lastpoint = TMath::Min(s1->GetLastPoint(),s2->GetLastPoint());
2146 Int_t firstpoint = 0;
2147 Int_t lastpoint = 160;
2149 // if (firstpoint>=lastpoint-5) return;;
2151 for (Int_t i=firstpoint;i<lastpoint;i++){
2152 // if ( (s1->GetClusterIndex2(i)&0xFFFF8FFF)==(s2->GetClusterIndex2(i)&0xFFFF8FFF) && s1->GetClusterIndex2(i)>0) {
2153 if ( (s1->GetClusterIndex2(i))==(s2->GetClusterIndex2(i)) && s1->GetClusterIndex2(i)>0) {
2157 if (sumshared>cutN0){
2160 for (Int_t i=firstpoint;i<lastpoint;i++){
2161 // if ( (s1->GetClusterIndex2(i)&0xFFFF8FFF)==(s2->GetClusterIndex2(i)&0xFFFF8FFF) && s1->GetClusterIndex2(i)>0) {
2162 if ( (s1->GetClusterIndex2(i))==(s2->GetClusterIndex2(i)) && s1->GetClusterIndex2(i)>0) {
2163 AliTPCTrackerPoint *p1 = s1->GetTrackPoint(i);
2164 AliTPCTrackerPoint *p2 = s2->GetTrackPoint(i);;
2165 if (s1->IsActive()&&s2->IsActive()){
2166 p1->SetShared(kTRUE);
2167 p2->SetShared(kTRUE);
2173 if (sumshared>cutN0){
2174 for (Int_t i=0;i<4;i++){
2175 if (s1->GetOverlapLabel(3*i)==0){
2176 s1->SetOverlapLabel(3*i, s2->GetLabel());
2177 s1->SetOverlapLabel(3*i+1,sumshared);
2178 s1->SetOverlapLabel(3*i+2,s2->GetUniqueID());
2182 for (Int_t i=0;i<4;i++){
2183 if (s2->GetOverlapLabel(3*i)==0){
2184 s2->SetOverlapLabel(3*i, s1->GetLabel());
2185 s2->SetOverlapLabel(3*i+1,sumshared);
2186 s2->SetOverlapLabel(3*i+2,s1->GetUniqueID());
2193 void AliTPCtrackerMI::SignShared(TObjArray * arr)
2196 //sort trackss according sectors
2198 for (Int_t i=0; i<arr->GetEntriesFast(); i++) {
2199 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2201 //if (pt) RotateToLocal(pt);
2205 arr->Sort(); // sorting according relative sectors
2206 arr->Expand(arr->GetEntries());
2209 Int_t nseed=arr->GetEntriesFast();
2210 for (Int_t i=0; i<nseed; i++) {
2211 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2213 for (Int_t j=0;j<12;j++){
2214 pt->SetOverlapLabel(j,0);
2217 for (Int_t i=0; i<nseed; i++) {
2218 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2220 if (pt->GetRemoval()>10) continue;
2221 for (Int_t j=i+1; j<nseed; j++){
2222 AliTPCseed *pt2=(AliTPCseed*)arr->UncheckedAt(j);
2223 if (TMath::Abs(pt->GetRelativeSector()-pt2->GetRelativeSector())>1) continue;
2225 if (pt2->GetRemoval()<=10) {
2226 //if ( TMath::Abs(pt->GetRelativeSector()-pt2->GetRelativeSector())>0) break;
2234 void AliTPCtrackerMI::SortTracks(TObjArray * arr, Int_t mode) const
2237 //sort tracks in array according mode criteria
2238 Int_t nseed = arr->GetEntriesFast();
2239 for (Int_t i=0; i<nseed; i++) {
2240 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2251 void AliTPCtrackerMI::RemoveUsed2(TObjArray * arr, Float_t factor1, Float_t factor2, Int_t minimal)
2254 // Loop over all tracks and remove overlaped tracks (with lower quality)
2256 // 1. Unsign clusters
2257 // 2. Sort tracks according quality
2258 // Quality is defined by the number of cluster between first and last points
2260 // 3. Loop over tracks - decreasing quality order
2261 // a.) remove - If the fraction of shared cluster less than factor (1- n or 2)
2262 // b.) remove - If the minimal number of clusters less than minimal and not ITS
2263 // c.) if track accepted - sign clusters
2265 //Called in - AliTPCtrackerMI::Clusters2Tracks()
2266 // - AliTPCtrackerMI::PropagateBack()
2267 // - AliTPCtrackerMI::RefitInward()
2270 // factor1 - factor for constrained
2271 // factor2 - for non constrained tracks
2272 // if (Float_t(shared+1)/Float_t(found+1)>factor) - DELETE
2276 Int_t nseed = arr->GetEntriesFast();
2277 Float_t * quality = new Float_t[nseed];
2278 Int_t * indexes = new Int_t[nseed];
2282 for (Int_t i=0; i<nseed; i++) {
2283 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2288 pt->UpdatePoints(); //select first last max dens points
2289 Float_t * points = pt->GetPoints();
2290 if (points[3]<0.8) quality[i] =-1;
2291 quality[i] = (points[2]-points[0])+pt->GetNumberOfClusters();
2292 //prefer high momenta tracks if overlaps
2293 quality[i] *= TMath::Sqrt(TMath::Abs(pt->Pt())+0.5);
2295 TMath::Sort(nseed,quality,indexes);
2298 for (Int_t itrack=0; itrack<nseed; itrack++) {
2299 Int_t trackindex = indexes[itrack];
2300 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(trackindex);
2303 if (quality[trackindex]<0){
2304 MarkSeedFree( arr->RemoveAt(trackindex) );
2309 Int_t first = Int_t(pt->GetPoints()[0]);
2310 Int_t last = Int_t(pt->GetPoints()[2]);
2311 Double_t factor = (pt->GetBConstrain()) ? factor1: factor2;
2313 Int_t found,foundable,shared;
2314 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
2315 // pt->GetClusterStatistic(0,160, found, foundable,shared,kFALSE);
2316 Bool_t itsgold =kFALSE;
2319 if (pt->GetESD()->GetITSclusters(dummy)>4) itsgold= kTRUE;
2323 if (Float_t(shared+1)/Float_t(found+1)>factor){
2324 if (pt->GetKinkIndexes()[0]!=0) continue; //don't remove tracks - part of the kinks
2325 if( AliTPCReconstructor::StreamLevel()>3){
2326 TTreeSRedirector &cstream = *fDebugStreamer;
2327 cstream<<"RemoveUsed"<<
2328 "iter="<<fIteration<<
2332 MarkSeedFree( arr->RemoveAt(trackindex) );
2335 if (pt->GetNumberOfClusters()<50&&(found-0.5*shared)<minimal){ //remove short tracks
2336 if (pt->GetKinkIndexes()[0]!=0) continue; //don't remove tracks - part of the kinks
2337 if( AliTPCReconstructor::StreamLevel()>3){
2338 TTreeSRedirector &cstream = *fDebugStreamer;
2339 cstream<<"RemoveShort"<<
2340 "iter="<<fIteration<<
2344 MarkSeedFree( arr->RemoveAt(trackindex) );
2350 //if (sharedfactor>0.4) continue;
2351 if (pt->GetKinkIndexes()[0]>0) continue;
2352 //Remove tracks with undefined properties - seems
2353 if (pt->GetSigmaY2()<kAlmost0) continue; // ? what is the origin ?
2355 for (Int_t i=first; i<last; i++) {
2356 Int_t index=pt->GetClusterIndex2(i);
2357 // if (index<0 || index&0x8000 ) continue;
2358 if (index<0 || index&0x8000 ) continue;
2359 AliTPCclusterMI *c= pt->GetClusterPointer(i);
2366 Info("RemoveUsed","\n*****\nNumber of good tracks after shared removal\t%d\n",fNtracks);
2372 void AliTPCtrackerMI::DumpClusters(Int_t iter, TObjArray *trackArray)
2375 // Dump clusters after reco
2376 // signed and unsigned cluster can be visualized
2377 // 1. Unsign all cluster
2378 // 2. Sign all used clusters
2381 Int_t nseed = trackArray->GetEntries();
2382 for (Int_t i=0; i<nseed; i++){
2383 AliTPCseed *pt=(AliTPCseed*)trackArray->UncheckedAt(i);
2387 Bool_t isKink=pt->GetKinkIndex(0)!=0;
2388 for (Int_t j=0; j<160; ++j) {
2389 Int_t index=pt->GetClusterIndex2(j);
2390 if (index<0) continue;
2391 AliTPCclusterMI *c= pt->GetClusterPointer(j);
2393 if (isKink) c->Use(100); // kink
2394 c->Use(10); // by default usage 10
2399 for (Int_t sec=0;sec<fkNIS;sec++){
2400 for (Int_t row=0;row<fInnerSec->GetNRows();row++){
2401 AliTPCclusterMI *cl = fInnerSec[sec][row].GetClusters1();
2402 for (Int_t icl =0;icl< fInnerSec[sec][row].GetN1();icl++){
2403 Float_t gx[3]; cl[icl].GetGlobalXYZ(gx);
2404 (*fDebugStreamer)<<"clDump"<<
2412 cl = fInnerSec[sec][row].GetClusters2();
2413 for (Int_t icl =0;icl< fInnerSec[sec][row].GetN2();icl++){
2414 Float_t gx[3]; cl[icl].GetGlobalXYZ(gx);
2415 (*fDebugStreamer)<<"clDump"<<
2426 for (Int_t sec=0;sec<fkNOS;sec++){
2427 for (Int_t row=0;row<fOuterSec->GetNRows();row++){
2428 AliTPCclusterMI *cl = fOuterSec[sec][row].GetClusters1();
2429 for (Int_t icl =0;icl< fOuterSec[sec][row].GetN1();icl++){
2430 Float_t gx[3]; cl[icl].GetGlobalXYZ(gx);
2431 (*fDebugStreamer)<<"clDump"<<
2439 cl = fOuterSec[sec][row].GetClusters2();
2440 for (Int_t icl =0;icl< fOuterSec[sec][row].GetN2();icl++){
2441 Float_t gx[3]; cl[icl].GetGlobalXYZ(gx);
2442 (*fDebugStreamer)<<"clDump"<<
2454 void AliTPCtrackerMI::UnsignClusters()
2457 // loop over all clusters and unsign them
2460 for (Int_t sec=0;sec<fkNIS;sec++){
2461 for (Int_t row=0;row<fInnerSec->GetNRows();row++){
2462 AliTPCclusterMI *cl = fInnerSec[sec][row].GetClusters1();
2463 for (Int_t icl =0;icl< fInnerSec[sec][row].GetN1();icl++)
2464 // if (cl[icl].IsUsed(10))
2466 cl = fInnerSec[sec][row].GetClusters2();
2467 for (Int_t icl =0;icl< fInnerSec[sec][row].GetN2();icl++)
2468 //if (cl[icl].IsUsed(10))
2473 for (Int_t sec=0;sec<fkNOS;sec++){
2474 for (Int_t row=0;row<fOuterSec->GetNRows();row++){
2475 AliTPCclusterMI *cl = fOuterSec[sec][row].GetClusters1();
2476 for (Int_t icl =0;icl< fOuterSec[sec][row].GetN1();icl++)
2477 //if (cl[icl].IsUsed(10))
2479 cl = fOuterSec[sec][row].GetClusters2();
2480 for (Int_t icl =0;icl< fOuterSec[sec][row].GetN2();icl++)
2481 //if (cl[icl].IsUsed(10))
2490 void AliTPCtrackerMI::SignClusters(const TObjArray * arr, Float_t fnumber, Float_t fdensity)
2493 //sign clusters to be "used"
2495 // snumber and sdensity sign number of sigmas - bellow mean value to be accepted
2496 // loop over "primaries"
2510 Int_t nseed = arr->GetEntriesFast();
2511 for (Int_t i=0; i<nseed; i++) {
2512 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2516 if (!(pt->IsActive())) continue;
2517 Float_t dens = pt->GetNumberOfClusters()/Float_t(pt->GetNFoundable());
2518 if ( (dens>0.7) && (pt->GetNumberOfClusters()>70)){
2520 sumdens2+= dens*dens;
2521 sumn += pt->GetNumberOfClusters();
2522 sumn2 += pt->GetNumberOfClusters()*pt->GetNumberOfClusters();
2523 Float_t chi2 = pt->GetChi2()/pt->GetNumberOfClusters();
2526 sumchi2 +=chi2*chi2;
2531 Float_t mdensity = 0.9;
2532 Float_t meann = 130;
2533 Float_t meanchi = 1;
2534 Float_t sdensity = 0.1;
2535 Float_t smeann = 10;
2536 Float_t smeanchi =0.4;
2540 mdensity = sumdens/sum;
2542 meanchi = sumchi/sum;
2544 sdensity = sumdens2/sum-mdensity*mdensity;
2546 sdensity = TMath::Sqrt(sdensity);
2550 smeann = sumn2/sum-meann*meann;
2552 smeann = TMath::Sqrt(smeann);
2556 smeanchi = sumchi2/sum - meanchi*meanchi;
2558 smeanchi = TMath::Sqrt(smeanchi);
2564 //REMOVE SHORT DELTAS or tracks going out of sensitive volume of TPC
2566 for (Int_t i=0; i<nseed; i++) {
2567 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2571 if (pt->GetBSigned()) continue;
2572 if (pt->GetBConstrain()) continue;
2573 //if (!(pt->IsActive())) continue;
2575 Int_t found,foundable,shared;
2576 pt->GetClusterStatistic(0,160,found, foundable,shared);
2577 if (shared/float(found)>0.3) {
2578 if (shared/float(found)>0.9 ){
2579 //MarkSeedFree( arr->RemoveAt(i) );
2584 Bool_t isok =kFALSE;
2585 if ( (pt->GetNShared()/pt->GetNumberOfClusters()<0.5) &&pt->GetNumberOfClusters()>60)
2587 if ((TMath::Abs(1/pt->GetC())<100.) && (pt->GetNShared()/pt->GetNumberOfClusters()<0.7))
2589 if (TMath::Abs(pt->GetZ()/pt->GetX())>1.1)
2591 if ( (TMath::Abs(pt->GetSnp()>0.7) && pt->GetD(0,0)>60.))
2595 for (Int_t j=0; j<160; ++j) {
2596 Int_t index=pt->GetClusterIndex2(j);
2597 if (index<0) continue;
2598 AliTPCclusterMI *c= pt->GetClusterPointer(j);
2600 //if (!(c->IsUsed(10))) c->Use();
2607 Double_t maxchi = meanchi+2.*smeanchi;
2609 for (Int_t i=0; i<nseed; i++) {
2610 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2614 //if (!(pt->IsActive())) continue;
2615 if (pt->GetBSigned()) continue;
2616 Double_t chi = pt->GetChi2()/pt->GetNumberOfClusters();
2617 if (chi>maxchi) continue;
2620 Float_t dens = pt->GetNumberOfClusters()/Float_t(pt->GetNFoundable());
2622 //sign only tracks with enoug big density at the beginning
2624 if ((pt->GetDensityFirst(40)<0.75) && pt->GetNumberOfClusters()<meann) continue;
2627 Double_t mindens = TMath::Max(double(mdensity-sdensity*fdensity*bfactor),0.65);
2628 Double_t minn = TMath::Max(Int_t(meann-fnumber*smeann*bfactor),50);
2630 // if (pt->fBConstrain) mindens = TMath::Max(mdensity-sdensity*fdensity*bfactor,0.65);
2631 if ( (pt->GetRemoval()==10) && (pt->GetSnp()>0.8)&&(dens>mindens))
2634 if ((dens>mindens && pt->GetNumberOfClusters()>minn) && chi<maxchi ){
2635 //Int_t noc=pt->GetNumberOfClusters();
2636 pt->SetBSigned(kTRUE);
2637 for (Int_t j=0; j<160; ++j) {
2639 Int_t index=pt->GetClusterIndex2(j);
2640 if (index<0) continue;
2641 AliTPCclusterMI *c= pt->GetClusterPointer(j);
2643 // if (!(c->IsUsed(10))) c->Use();
2648 // gLastCheck = nseed;
2656 void AliTPCtrackerMI::StopNotActive(const TObjArray * arr, Int_t row0, Float_t th0, Float_t th1, Float_t th2) const
2658 // stop not active tracks
2659 // take th1 as threshold for number of founded to number of foundable on last 10 active rows
2660 // take th2 as threshold for number of founded to number of foundable on last 20 active rows
2661 Int_t nseed = arr->GetEntriesFast();
2663 for (Int_t i=0; i<nseed; i++) {
2664 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2668 if (!(pt->IsActive())) continue;
2669 StopNotActive(pt,row0,th0, th1,th2);
2675 void AliTPCtrackerMI::StopNotActive(AliTPCseed * seed, Int_t row0, Float_t th0, Float_t th1,
2678 // stop not active tracks
2679 // take th1 as threshold for number of founded to number of foundable on last 10 active rows
2680 // take th2 as threshold for number of founded to number of foundable on last 20 active rows
2683 Int_t foundable = 0;
2684 Int_t maxindex = seed->GetLastPoint(); //last foundable row
2685 if (seed->GetNFoundable()*th0 > seed->GetNumberOfClusters()) {
2686 seed->Desactivate(10) ;
2690 for (Int_t i=row0; i<maxindex; i++){
2691 Int_t index = seed->GetClusterIndex2(i);
2692 if (index!=-1) foundable++;
2694 if (foundable<=30) sumgood1++;
2695 if (foundable<=50) {
2702 if (foundable>=30.){
2703 if (sumgood1<(th1*30.)) seed->Desactivate(10);
2706 if (sumgood2<(th2*50.)) seed->Desactivate(10);
2710 Int_t AliTPCtrackerMI::RefitInward(AliESDEvent *event)
2713 // back propagation of ESD tracks
2716 if (!event) return 0;
2717 const Int_t kMaxFriendTracks=2000;
2719 // extract correction object for multiplicity dependence of dEdx
2720 TObjArray * gainCalibArray = AliTPCcalibDB::Instance()->GetTimeGainSplinesRun(event->GetRunNumber());
2722 AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ;
2724 AliFatal("Tranformations not in RefitInward");
2727 transform->SetCurrentRecoParam((AliTPCRecoParam*)AliTPCReconstructor::GetRecoParam());
2728 const AliTPCRecoParam * recoParam = AliTPCcalibDB::Instance()->GetTransform()->GetCurrentRecoParam();
2729 Int_t nContribut = event->GetNumberOfTracks();
2730 TGraphErrors * graphMultDependenceDeDx = 0x0;
2731 if (recoParam && recoParam->GetUseMultiplicityCorrectionDedx() && gainCalibArray) {
2732 if (recoParam->GetUseTotCharge()) {
2733 graphMultDependenceDeDx = (TGraphErrors *) gainCalibArray->FindObject("TGRAPHERRORS_MEANQTOT_MULTIPLICITYDEPENDENCE_BEAM_ALL");
2735 graphMultDependenceDeDx = (TGraphErrors *) gainCalibArray->FindObject("TGRAPHERRORS_MEANQMAX_MULTIPLICITYDEPENDENCE_BEAM_ALL");
2741 //PrepareForProlongation(fSeeds,1);
2742 PropagateForward2(fSeeds);
2743 RemoveUsed2(fSeeds,0.4,0.4,20);
2745 TObjArray arraySeed(fSeeds->GetEntries());
2746 for (Int_t i=0;i<fSeeds->GetEntries();i++) {
2747 arraySeed.AddAt(fSeeds->At(i),i);
2749 SignShared(&arraySeed);
2750 // FindCurling(fSeeds, event,2); // find multi found tracks
2751 FindSplitted(fSeeds, event,2); // find multi found tracks
2752 if (AliTPCReconstructor::StreamLevel()>5) FindMultiMC(fSeeds, fEvent,2); // find multi found tracks
2755 Int_t nseed = fSeeds->GetEntriesFast();
2756 for (Int_t i=0;i<nseed;i++){
2757 AliTPCseed * seed = (AliTPCseed*) fSeeds->UncheckedAt(i);
2758 if (!seed) continue;
2759 if (seed->GetKinkIndex(0)>0) UpdateKinkQualityD(seed); // update quality informations for kinks
2760 AliESDtrack *esd=event->GetTrack(i);
2762 if (seed->GetNumberOfClusters()<60 && seed->GetNumberOfClusters()<(esd->GetTPCclusters(0) -5)*0.8){
2763 AliExternalTrackParam paramIn;
2764 AliExternalTrackParam paramOut;
2765 Int_t ncl = seed->RefitTrack(seed,¶mIn,¶mOut);
2766 if (AliTPCReconstructor::StreamLevel()>2) {
2767 (*fDebugStreamer)<<"RecoverIn"<<
2771 "pout.="<<¶mOut<<
2776 seed->Set(paramIn.GetX(),paramIn.GetAlpha(),paramIn.GetParameter(),paramIn.GetCovariance());
2777 seed->SetNumberOfClusters(ncl);
2781 seed->PropagateTo(fkParam->GetInnerRadiusLow());
2782 seed->UpdatePoints();
2783 AddCovariance(seed);
2784 MakeESDBitmaps(seed, esd);
2785 seed->CookdEdx(0.02,0.6);
2786 CookLabel(seed,0.1); //For comparison only
2788 if (AliTPCReconstructor::StreamLevel()>1 && seed!=0) {
2789 TTreeSRedirector &cstream = *fDebugStreamer;
2796 if (seed->GetNumberOfClusters()>15){
2797 esd->UpdateTrackParams(seed,AliESDtrack::kTPCrefit);
2798 esd->SetTPCPoints(seed->GetPoints());
2799 esd->SetTPCPointsF(seed->GetNFoundable());
2800 Int_t ndedx = seed->GetNCDEDX(0);
2801 Float_t sdedx = seed->GetSDEDX(0);
2802 Float_t dedx = seed->GetdEdx();
2803 // apply mutliplicity dependent dEdx correction if available
2804 if (graphMultDependenceDeDx) {
2805 Double_t corrGain = AliTPCcalibDButil::EvalGraphConst(graphMultDependenceDeDx, nContribut);
2806 dedx += (1 - corrGain)*50.; // MIP is normalized to 50
2808 esd->SetTPCsignal(dedx, sdedx, ndedx);
2810 // fill new dEdx information
2812 Double32_t signal[4];
2816 for(Int_t iarr=0;iarr<3;iarr++) {
2817 signal[iarr] = seed->GetDEDXregion(iarr+1);
2818 ncl[iarr] = seed->GetNCDEDX(iarr+1);
2819 nrows[iarr] = seed->GetNCDEDXInclThres(iarr+1);
2821 signal[3] = seed->GetDEDXregion(4);
2823 AliTPCdEdxInfo * infoTpcPid = new AliTPCdEdxInfo();
2824 infoTpcPid->SetTPCSignalRegionInfo(signal, ncl, nrows);
2825 esd->SetTPCdEdxInfo(infoTpcPid);
2827 // add seed to the esd track in Calib level
2829 Bool_t storeFriend = gRandom->Rndm()<(kMaxFriendTracks)/Float_t(nseed);
2830 if (AliTPCReconstructor::StreamLevel()>0 &&storeFriend){
2831 // RS: this is the only place where the seed is created not in the pool,
2832 // since it should belong to ESDevent
2833 AliTPCseed * seedCopy = new AliTPCseed(*seed, kTRUE);
2834 esd->AddCalibObject(seedCopy);
2839 //printf("problem\n");
2842 //FindKinks(fSeeds,event);
2843 if (AliTPCReconstructor::StreamLevel()>3) DumpClusters(2,fSeeds);
2844 Info("RefitInward","Number of refitted tracks %d",ntracks);
2849 Int_t AliTPCtrackerMI::PropagateBack(AliESDEvent *event)
2852 // back propagation of ESD tracks
2854 if (!event) return 0;
2858 PropagateBack(fSeeds);
2859 RemoveUsed2(fSeeds,0.4,0.4,20);
2860 //FindCurling(fSeeds, fEvent,1);
2861 FindSplitted(fSeeds, event,1); // find multi found tracks
2862 if (AliTPCReconstructor::StreamLevel()>5) FindMultiMC(fSeeds, fEvent,1); // find multi found tracks
2865 Int_t nseed = fSeeds->GetEntriesFast();
2867 for (Int_t i=0;i<nseed;i++){
2868 AliTPCseed * seed = (AliTPCseed*) fSeeds->UncheckedAt(i);
2869 if (!seed) continue;
2870 if (seed->GetKinkIndex(0)<0) UpdateKinkQualityM(seed); // update quality informations for kinks
2871 seed->UpdatePoints();
2872 AddCovariance(seed);
2873 AliESDtrack *esd=event->GetTrack(i);
2874 if (!esd) continue; //never happen
2875 if (seed->GetNumberOfClusters()<60 && seed->GetNumberOfClusters()<(esd->GetTPCclusters(0) -5)*0.8){
2876 AliExternalTrackParam paramIn;
2877 AliExternalTrackParam paramOut;
2878 Int_t ncl = seed->RefitTrack(seed,¶mIn,¶mOut);
2879 if (AliTPCReconstructor::StreamLevel()>2) {
2880 (*fDebugStreamer)<<"RecoverBack"<<
2884 "pout.="<<¶mOut<<
2889 seed->Set(paramOut.GetX(),paramOut.GetAlpha(),paramOut.GetParameter(),paramOut.GetCovariance());
2890 seed->SetNumberOfClusters(ncl);
2893 seed->CookdEdx(0.02,0.6);
2894 CookLabel(seed,0.1); //For comparison only
2895 if (seed->GetNumberOfClusters()>15){
2896 esd->UpdateTrackParams(seed,AliESDtrack::kTPCout);
2897 esd->SetTPCPoints(seed->GetPoints());
2898 esd->SetTPCPointsF(seed->GetNFoundable());
2899 Int_t ndedx = seed->GetNCDEDX(0);
2900 Float_t sdedx = seed->GetSDEDX(0);
2901 Float_t dedx = seed->GetdEdx();
2902 esd->SetTPCsignal(dedx, sdedx, ndedx);
2904 Int_t eventnumber = event->GetEventNumberInFile();// patch 28 fev 06
2905 // This is most likely NOT the event number you'd like to use. It has nothing to do with the 'real' event number
2906 if (AliTPCReconstructor::StreamLevel()>1 && esd) {
2907 (*fDebugStreamer)<<"Cback"<<
2910 "EventNrInFile="<<eventnumber<<
2915 if (AliTPCReconstructor::StreamLevel()>3) DumpClusters(1,fSeeds);
2916 //FindKinks(fSeeds,event);
2917 Info("PropagateBack","Number of back propagated tracks %d",ntracks);
2924 void AliTPCtrackerMI::DeleteSeeds()
2933 void AliTPCtrackerMI::ReadSeeds(const AliESDEvent *const event, Int_t direction)
2936 //read seeds from the event
2938 Int_t nentr=event->GetNumberOfTracks();
2940 Info("ReadSeeds", "Number of ESD tracks: %d\n", nentr);
2945 fSeeds = new TObjArray(nentr);
2949 for (Int_t i=0; i<nentr; i++) {
2950 AliESDtrack *esd=event->GetTrack(i);
2951 ULong_t status=esd->GetStatus();
2952 if (!(status&AliESDtrack::kTPCin)) continue;
2953 AliTPCtrack t(*esd);
2954 t.SetNumberOfClusters(0);
2955 // AliTPCseed *seed = new AliTPCseed(t,t.GetAlpha());
2956 AliTPCseed *seed = new( NextFreeSeed() ) AliTPCseed(t/*,t.GetAlpha()*/);
2957 seed->SetPoolID(fLastSeedID);
2958 seed->SetUniqueID(esd->GetID());
2959 AddCovariance(seed); //add systematic ucertainty
2960 for (Int_t ikink=0;ikink<3;ikink++) {
2961 Int_t index = esd->GetKinkIndex(ikink);
2962 seed->GetKinkIndexes()[ikink] = index;
2963 if (index==0) continue;
2964 index = TMath::Abs(index);
2965 AliESDkink * kink = fEvent->GetKink(index-1);
2966 if (kink&&esd->GetKinkIndex(ikink)<0){
2967 if ((status & AliESDtrack::kTRDrefit) != 0) kink->SetStatus(1,2);
2968 if ((status & AliESDtrack::kITSout) != 0) kink->SetStatus(1,0);
2970 if (kink&&esd->GetKinkIndex(ikink)>0){
2971 if ((status & AliESDtrack::kTRDrefit) != 0) kink->SetStatus(1,6);
2972 if ((status & AliESDtrack::kITSout) != 0) kink->SetStatus(1,4);
2976 if (((status&AliESDtrack::kITSout)==0)&&(direction==1)) seed->ResetCovariance(10.);
2977 if ( direction ==2 &&(status & AliESDtrack::kTRDrefit) == 0 ) seed->ResetCovariance(10.);
2978 //if ( direction ==2 && ((status & AliESDtrack::kTPCout) == 0) ) {
2979 // fSeeds->AddAt(0,i);
2980 // MarkSeedFree( seed );
2983 if ( direction ==2 &&(status & AliESDtrack::kTRDrefit) > 0 ) {
2984 Double_t par0[5],par1[5],alpha,x;
2985 esd->GetInnerExternalParameters(alpha,x,par0);
2986 esd->GetExternalParameters(x,par1);
2987 Double_t delta1 = TMath::Abs(par0[4]-par1[4])/(0.000000001+TMath::Abs(par0[4]+par1[4]));
2988 Double_t delta2 = TMath::Abs(par0[3]-par1[3]);
2990 if (esd->GetTRDncls()>0) trdchi2 = esd->GetTRDchi2()/esd->GetTRDncls();
2991 //reset covariance if suspicious
2992 if ( (delta1>0.1) || (delta2>0.006) ||trdchi2>7.)
2993 seed->ResetCovariance(10.);
2998 // rotate to the local coordinate system
3000 fSectors=fInnerSec; fN=fkNIS;
3001 Double_t alpha=seed->GetAlpha();
3002 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
3003 if (alpha < 0. ) alpha += 2.*TMath::Pi();
3004 Int_t ns=Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN;
3005 alpha =ns*fSectors->GetAlpha() + fSectors->GetAlphaShift();
3006 alpha-=seed->GetAlpha();
3007 if (alpha<-TMath::Pi()) alpha += 2*TMath::Pi();
3008 if (alpha>=TMath::Pi()) alpha -= 2*TMath::Pi();
3009 if (TMath::Abs(alpha) > 0.001) { // This should not happen normally
3010 AliWarning(Form("Rotating track over %f",alpha));
3011 if (!seed->Rotate(alpha)) {
3012 MarkSeedFree( seed );
3018 if (esd->GetKinkIndex(0)<=0){
3019 for (Int_t irow=0;irow<160;irow++){
3020 Int_t index = seed->GetClusterIndex2(irow);
3023 AliTPCclusterMI * cl = GetClusterMI(index);
3024 seed->SetClusterPointer(irow,cl);
3026 if ((index & 0x8000)==0){
3027 cl->Use(10); // accepted cluster
3029 cl->Use(6); // close cluster not accepted
3032 Info("ReadSeeds","Not found cluster");
3037 fSeeds->AddAt(seed,i);
3043 //_____________________________________________________________________________
3044 void AliTPCtrackerMI::MakeSeeds3(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, Float_t cuts[4],
3045 Float_t deltay, Int_t ddsec) {
3046 //-----------------------------------------------------------------
3047 // This function creates track seeds.
3048 // SEEDING WITH VERTEX CONSTRAIN
3049 //-----------------------------------------------------------------
3050 // cuts[0] - fP4 cut
3051 // cuts[1] - tan(phi) cut
3052 // cuts[2] - zvertex cut
3053 // cuts[3] - fP3 cut
3061 Double_t x[5], c[15];
3062 // Int_t di = i1-i2;
3064 AliTPCseed * seed = new( NextFreeSeed() ) AliTPCseed();
3065 seed->SetPoolID(fLastSeedID);
3066 Double_t alpha=fSectors->GetAlpha(), shift=fSectors->GetAlphaShift();
3067 Double_t cs=cos(alpha), sn=sin(alpha);
3069 // Double_t x1 =fOuterSec->GetX(i1);
3070 //Double_t xx2=fOuterSec->GetX(i2);
3072 Double_t x1 =GetXrow(i1);
3073 Double_t xx2=GetXrow(i2);
3075 Double_t x3=GetX(), y3=GetY(), z3=GetZ();
3077 Int_t imiddle = (i2+i1)/2; //middle pad row index
3078 Double_t xm = GetXrow(imiddle); // radius of middle pad-row
3079 const AliTPCtrackerRow& krm=GetRow(sec,imiddle); //middle pad -row
3083 const AliTPCtrackerRow& kr1=GetRow(ns,i1);
3084 Double_t ymax = GetMaxY(i1)-kr1.GetDeadZone()-1.5;
3085 Double_t ymaxm = GetMaxY(imiddle)-kr1.GetDeadZone()-1.5;
3088 // change cut on curvature if it can't reach this layer
3089 // maximal curvature set to reach it
3090 Double_t dvertexmax = TMath::Sqrt((x1-x3)*(x1-x3)+(ymax+5-y3)*(ymax+5-y3));
3091 if (dvertexmax*0.5*cuts[0]>0.85){
3092 cuts[0] = 0.85/(dvertexmax*0.5+1.);
3094 Double_t r2min = 1/(cuts[0]*cuts[0]); //minimal square of radius given by cut
3097 if (deltay>0) ddsec = 0;
3098 // loop over clusters
3099 for (Int_t is=0; is < kr1; is++) {
3101 if (kr1[is]->IsUsed(10)) continue;
3102 Double_t y1=kr1[is]->GetY(), z1=kr1[is]->GetZ();
3103 //if (TMath::Abs(y1)>ymax) continue;
3105 if (deltay>0 && TMath::Abs(ymax-TMath::Abs(y1))> deltay ) continue; // seed only at the edge
3107 // find possible directions
3108 Float_t anglez = (z1-z3)/(x1-x3);
3109 Float_t extraz = z1 - anglez*(x1-xx2); // extrapolated z
3112 //find rotation angles relative to line given by vertex and point 1
3113 Double_t dvertex2 = (x1-x3)*(x1-x3)+(y1-y3)*(y1-y3);
3114 Double_t dvertex = TMath::Sqrt(dvertex2);
3115 Double_t angle13 = TMath::ATan((y1-y3)/(x1-x3));
3116 Double_t cs13 = cos(-angle13), sn13 = sin(-angle13);
3119 // loop over 2 sectors
3125 Double_t dddz1=0; // direction of delta inclination in z axis
3132 for (Int_t dsec = dsec1; dsec<=dsec2;dsec++){
3133 Int_t sec2 = sec + dsec;
3135 // AliTPCtrackerRow& kr2 = fOuterSec[(sec2+fkNOS)%fkNOS][i2];
3136 //AliTPCtrackerRow& kr2m = fOuterSec[(sec2+fkNOS)%fkNOS][imiddle];
3137 AliTPCtrackerRow& kr2 = GetRow((sec2+fkNOS)%fkNOS,i2);
3138 AliTPCtrackerRow& kr2m = GetRow((sec2+fkNOS)%fkNOS,imiddle);
3139 Int_t index1 = TMath::Max(kr2.Find(extraz-0.6-dddz1*TMath::Abs(z1)*0.05)-1,0);
3140 Int_t index2 = TMath::Min(kr2.Find(extraz+0.6+dddz2*TMath::Abs(z1)*0.05)+1,kr2);
3142 // rotation angles to p1-p3
3143 Double_t cs13r = cos(-angle13+dsec*alpha)/dvertex, sn13r = sin(-angle13+dsec*alpha)/dvertex;
3144 Double_t x2, y2, z2;
3146 // Double_t dymax = maxangle*TMath::Abs(x1-xx2);
3149 Double_t dxx0 = (xx2-x3)*cs13r;
3150 Double_t dyy0 = (xx2-x3)*sn13r;
3151 for (Int_t js=index1; js < index2; js++) {
3152 const AliTPCclusterMI *kcl = kr2[js];
3153 if (kcl->IsUsed(10)) continue;
3155 //calcutate parameters
3157 Double_t yy0 = dyy0 +(kcl->GetY()-y3)*cs13r;
3159 if (TMath::Abs(yy0)<0.000001) continue;
3160 Double_t xx0 = dxx0 -(kcl->GetY()-y3)*sn13r;
3161 Double_t y0 = 0.5*(xx0*xx0+yy0*yy0-xx0)/yy0;
3162 Double_t r02 = (0.25+y0*y0)*dvertex2;
3163 //curvature (radius) cut
3164 if (r02<r2min) continue;
3168 Double_t c0 = 1/TMath::Sqrt(r02);
3172 //Double_t dfi0 = 2.*TMath::ASin(dvertex*c0*0.5);
3173 //Double_t dfi1 = 2.*TMath::ASin(TMath::Sqrt(yy0*yy0+(1-xx0)*(1-xx0))*dvertex*c0*0.5);
3174 Double_t dfi0 = 2.*AliTPCFastMath::FastAsin(dvertex*c0*0.5);
3175 Double_t dfi1 = 2.*AliTPCFastMath::FastAsin(TMath::Sqrt(yy0*yy0+(1-xx0)*(1-xx0))*dvertex*c0*0.5);
3178 Double_t z0 = kcl->GetZ();
3179 Double_t zzzz2 = z1-(z1-z3)*dfi1/dfi0;
3180 if (TMath::Abs(zzzz2-z0)>0.5) continue;
3183 Double_t dip = (z1-z0)*c0/dfi1;
3184 Double_t x0 = (0.5*cs13+y0*sn13)*dvertex*c0;
3195 x2= xx2*cs-y2*sn*dsec;
3196 y2=+xx2*sn*dsec+y2*cs;
3206 // do we have cluster at the middle ?
3208 GetProlongation(x1,xm,x,ym,zm);
3210 AliTPCclusterMI * cm=0;
3211 if (TMath::Abs(ym)-ymaxm<0){
3212 cm = krm.FindNearest2(ym,zm,1.0,0.6,dummy);
3213 if ((!cm) || (cm->IsUsed(10))) {
3218 // rotate y1 to system 0
3219 // get state vector in rotated system
3220 Double_t yr1 = (-0.5*sn13+y0*cs13)*dvertex*c0;
3221 Double_t xr2 = x0*cs+yr1*sn*dsec;
3222 Double_t xr[5]={kcl->GetY(),kcl->GetZ(), xr2, dip, c0};
3224 GetProlongation(xx2,xm,xr,ym,zm);
3225 if (TMath::Abs(ym)-ymaxm<0){
3226 cm = kr2m.FindNearest2(ym,zm,1.0,0.6,dummy);
3227 if ((!cm) || (cm->IsUsed(10))) {
3237 dym = ym - cm->GetY();
3238 dzm = zm - cm->GetZ();
3245 Double_t sy1=kr1[is]->GetSigmaY2()*2., sz1=kr1[is]->GetSigmaZ2()*2.;
3246 Double_t sy2=kcl->GetSigmaY2()*2., sz2=kcl->GetSigmaZ2()*2.;
3247 //Double_t sy3=400*3./12., sy=0.1, sz=0.1;
3248 Double_t sy3=25000*x[4]*x[4]+0.1, sy=0.1, sz=0.1;
3249 //Double_t sy3=25000*x[4]*x[4]*60+0.5, sy=0.1, sz=0.1;
3251 Double_t f40=(F1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
3252 Double_t f42=(F1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
3253 Double_t f43=(F1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
3254 Double_t f20=(F2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
3255 Double_t f22=(F2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
3256 Double_t f23=(F2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
3258 Double_t f30=(F3(x1,y1+sy,x2,y2,z1,z2)-x[3])/sy;
3259 Double_t f31=(F3(x1,y1,x2,y2,z1+sz,z2)-x[3])/sz;
3260 Double_t f32=(F3(x1,y1,x2,y2+sy,z1,z2)-x[3])/sy;
3261 Double_t f34=(F3(x1,y1,x2,y2,z1,z2+sz)-x[3])/sz;
3265 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3266 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3267 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3268 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3269 c[13]=f30*sy1*f40+f32*sy2*f42;
3270 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3272 // if (!BuildSeed(kr1[is],kcl,0,x1,x2,x3,x,c)) continue;
3274 UInt_t index=kr1.GetIndex(is);
3275 if (seed) {MarkSeedFree(seed); seed = 0;}
3276 AliTPCseed *track = seed = new( NextFreeSeed() ) AliTPCseed(x1, ns*alpha+shift, x, c, index);
3277 seed->SetPoolID(fLastSeedID);
3278 track->SetIsSeeding(kTRUE);
3279 track->SetSeed1(i1);
3280 track->SetSeed2(i2);
3281 track->SetSeedType(3);
3285 FollowProlongation(*track, (i1+i2)/2,1);
3286 Int_t foundable,found,shared;
3287 track->GetClusterStatistic((i1+i2)/2,i1, found, foundable, shared, kTRUE);
3288 if ((found<0.55*foundable) || shared>0.5*found || (track->GetSigmaY2()+track->GetSigmaZ2())>0.5){
3289 MarkSeedFree(seed); seed = 0;
3295 FollowProlongation(*track, i2,1);
3299 track->SetBConstrain(1);
3300 // track->fLastPoint = i1+fInnerSec->GetNRows(); // first cluster in track position
3301 track->SetLastPoint(i1); // first cluster in track position
3302 track->SetFirstPoint(track->GetLastPoint());
3304 if (track->GetNumberOfClusters()<(i1-i2)*0.5 ||
3305 track->GetNumberOfClusters() < track->GetNFoundable()*0.6 ||
3306 track->GetNShared()>0.4*track->GetNumberOfClusters() ) {
3307 MarkSeedFree(seed); seed = 0;
3311 // Z VERTEX CONDITION
3312 Double_t zv, bz=GetBz();
3313 if ( !track->GetZAt(0.,bz,zv) ) continue;
3314 if (TMath::Abs(zv-z3)>cuts[2]) {
3315 FollowProlongation(*track, TMath::Max(i2-20,0));
3316 if ( !track->GetZAt(0.,bz,zv) ) continue;
3317 if (TMath::Abs(zv-z3)>cuts[2]){
3318 FollowProlongation(*track, TMath::Max(i2-40,0));
3319 if ( !track->GetZAt(0.,bz,zv) ) continue;
3320 if (TMath::Abs(zv-z3)>cuts[2] &&(track->GetNumberOfClusters() > track->GetNFoundable()*0.7)){
3321 // make seed without constrain
3322 AliTPCseed * track2 = MakeSeed(track,0.2,0.5,1.);
3323 FollowProlongation(*track2, i2,1);
3324 track2->SetBConstrain(kFALSE);
3325 track2->SetSeedType(1);
3326 arr->AddLast(track2);
3327 MarkSeedFree( seed ); seed = 0;
3331 MarkSeedFree( seed ); seed = 0;
3338 track->SetSeedType(0);
3339 arr->AddLast(track); // note, track is seed, don't free the seed
3340 seed = new( NextFreeSeed() ) AliTPCseed;
3341 seed->SetPoolID(fLastSeedID);
3343 // don't consider other combinations
3344 if (track->GetNumberOfClusters() > track->GetNFoundable()*0.8)
3350 Info("MakeSeeds3","\nSeeding statistic:\t%d\t%d\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin,nout1,nout2);
3352 if (seed) MarkSeedFree( seed );
3356 void AliTPCtrackerMI::MakeSeeds5(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, Float_t cuts[4],
3361 //-----------------------------------------------------------------
3362 // This function creates track seeds.
3363 //-----------------------------------------------------------------
3364 // cuts[0] - fP4 cut
3365 // cuts[1] - tan(phi) cut
3366 // cuts[2] - zvertex cut
3367 // cuts[3] - fP3 cut
3377 Double_t x[5], c[15];
3379 // make temporary seed
3380 AliTPCseed * seed = new( NextFreeSeed() ) AliTPCseed;
3381 seed->SetPoolID(fLastSeedID);
3382 Double_t alpha=fOuterSec->GetAlpha(), shift=fOuterSec->GetAlphaShift();
3383 // Double_t cs=cos(alpha), sn=sin(alpha);
3388 Double_t x1 = GetXrow(i1-1);
3389 const AliTPCtrackerRow& kr1=GetRow(sec,i1-1);
3390 Double_t y1max = GetMaxY(i1-1)-kr1.GetDeadZone()-1.5;
3392 Double_t x1p = GetXrow(i1);
3393 const AliTPCtrackerRow& kr1p=GetRow(sec,i1);
3395 Double_t x1m = GetXrow(i1-2);
3396 const AliTPCtrackerRow& kr1m=GetRow(sec,i1-2);
3399 //last 3 padrow for seeding
3400 AliTPCtrackerRow& kr3 = GetRow((sec+fkNOS)%fkNOS,i1-7);
3401 Double_t x3 = GetXrow(i1-7);
3402 // Double_t y3max= GetMaxY(i1-7)-kr3.fDeadZone-1.5;
3404 AliTPCtrackerRow& kr3p = GetRow((sec+fkNOS)%fkNOS,i1-6);
3405 Double_t x3p = GetXrow(i1-6);
3407 AliTPCtrackerRow& kr3m = GetRow((sec+fkNOS)%fkNOS,i1-8);
3408 Double_t x3m = GetXrow(i1-8);
3413 Int_t im = i1-4; //middle pad row index
3414 Double_t xm = GetXrow(im); // radius of middle pad-row
3415 const AliTPCtrackerRow& krm=GetRow(sec,im); //middle pad -row
3416 // Double_t ymmax = GetMaxY(im)-kr1.fDeadZone-1.5;
3419 Double_t deltax = x1-x3;
3420 Double_t dymax = deltax*cuts[1];
3421 Double_t dzmax = deltax*cuts[3];
3423 // loop over clusters
3424 for (Int_t is=0; is < kr1; is++) {
3426 if (kr1[is]->IsUsed(10)) continue;
3427 Double_t y1=kr1[is]->GetY(), z1=kr1[is]->GetZ();
3429 if (deltay>0 && TMath::Abs(y1max-TMath::Abs(y1))> deltay ) continue; // seed only at the edge
3431 Int_t index1 = TMath::Max(kr3.Find(z1-dzmax)-1,0);
3432 Int_t index2 = TMath::Min(kr3.Find(z1+dzmax)+1,kr3);
3438 for (Int_t js=index1; js < index2; js++) {
3439 const AliTPCclusterMI *kcl = kr3[js];
3440 if (kcl->IsUsed(10)) continue;
3442 // apply angular cuts
3443 if (TMath::Abs(y1-y3)>dymax) continue;
3446 if (TMath::Abs(z1-z3)>dzmax) continue;
3448 Double_t angley = (y1-y3)/(x1-x3);
3449 Double_t anglez = (z1-z3)/(x1-x3);
3451 Double_t erry = TMath::Abs(angley)*(x1-x1m)*0.5+0.5;
3452 Double_t errz = TMath::Abs(anglez)*(x1-x1m)*0.5+0.5;
3454 Double_t yyym = angley*(xm-x1)+y1;
3455 Double_t zzzm = anglez*(xm-x1)+z1;
3457 const AliTPCclusterMI *kcm = krm.FindNearest2(yyym,zzzm,erry,errz,index);
3459 if (kcm->IsUsed(10)) continue;
3461 erry = TMath::Abs(angley)*(x1-x1m)*0.4+0.5;
3462 errz = TMath::Abs(anglez)*(x1-x1m)*0.4+0.5;
3469 // look around first
3470 const AliTPCclusterMI *kc1m = kr1m.FindNearest2(angley*(x1m-x1)+y1,
3476 if (kc1m->IsUsed(10)) used++;
3478 const AliTPCclusterMI *kc1p = kr1p.FindNearest2(angley*(x1p-x1)+y1,
3484 if (kc1p->IsUsed(10)) used++;
3486 if (used>1) continue;
3487 if (found<1) continue;
3491 const AliTPCclusterMI *kc3m = kr3m.FindNearest2(angley*(x3m-x3)+y3,
3497 if (kc3m->IsUsed(10)) used++;
3501 const AliTPCclusterMI *kc3p = kr3p.FindNearest2(angley*(x3p-x3)+y3,
3507 if (kc3p->IsUsed(10)) used++;
3511 if (used>1) continue;
3512 if (found<3) continue;
3522 x[4]=F1(x1,y1,x2,y2,x3,y3);
3523 //if (TMath::Abs(x[4]) >= cuts[0]) continue;
3526 x[2]=F2(x1,y1,x2,y2,x3,y3);
3529 x[3]=F3n(x1,y1,x2,y2,z1,z2,x[4]);
3530 //if (TMath::Abs(x[3]) > cuts[3]) continue;
3534 Double_t sy1=0.1, sz1=0.1;
3535 Double_t sy2=0.1, sz2=0.1;
3536 Double_t sy3=0.1, sy=0.1, sz=0.1;
3538 Double_t f40=(F1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
3539 Double_t f42=(F1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
3540 Double_t f43=(F1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
3541 Double_t f20=(F2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
3542 Double_t f22=(F2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
3543 Double_t f23=(F2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
3545 Double_t f30=(F3(x1,y1+sy,x2,y2,z1,z2)-x[3])/sy;
3546 Double_t f31=(F3(x1,y1,x2,y2,z1+sz,z2)-x[3])/sz;
3547 Double_t f32=(F3(x1,y1,x2,y2+sy,z1,z2)-x[3])/sy;
3548 Double_t f34=(F3(x1,y1,x2,y2,z1,z2+sz)-x[3])/sz;
3552 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3553 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3554 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3555 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3556 c[13]=f30*sy1*f40+f32*sy2*f42;
3557 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3559 // if (!BuildSeed(kr1[is],kcl,0,x1,x2,x3,x,c)) continue;
3561 index=kr1.GetIndex(is);
3562 if (seed) {MarkSeedFree( seed ); seed = 0;}
3563 AliTPCseed *track = seed = new( NextFreeSeed() ) AliTPCseed(x1, sec*alpha+shift, x, c, index);
3564 seed->SetPoolID(fLastSeedID);
3566 track->SetIsSeeding(kTRUE);
3569 FollowProlongation(*track, i1-7,1);
3570 if (track->GetNumberOfClusters() < track->GetNFoundable()*0.75 ||
3571 track->GetNShared()>0.6*track->GetNumberOfClusters() || ( track->GetSigmaY2()+ track->GetSigmaZ2())>0.6){
3572 MarkSeedFree( seed ); seed = 0;
3578 FollowProlongation(*track, i2,1);
3579 track->SetBConstrain(0);
3580 track->SetLastPoint(i1+fInnerSec->GetNRows()); // first cluster in track position
3581 track->SetFirstPoint(track->GetLastPoint());
3583 if (track->GetNumberOfClusters()<(i1-i2)*0.5 ||
3584 track->GetNumberOfClusters()<track->GetNFoundable()*0.7 ||
3585 track->GetNShared()>2. || track->GetChi2()/track->GetNumberOfClusters()>6 || ( track->GetSigmaY2()+ track->GetSigmaZ2())>0.5 ) {
3586 MarkSeedFree( seed ); seed = 0;
3591 FollowProlongation(*track, TMath::Max(i2-10,0),1);
3592 AliTPCseed * track2 = MakeSeed(track,0.2,0.5,0.9);
3593 FollowProlongation(*track2, i2,1);
3594 track2->SetBConstrain(kFALSE);
3595 track2->SetSeedType(4);
3596 arr->AddLast(track2);
3597 MarkSeedFree( seed ); seed = 0;
3601 //arr->AddLast(track);
3602 //seed = new AliTPCseed;
3608 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);
3610 if (seed) MarkSeedFree(seed);
3614 //_____________________________________________________________________________
3615 void AliTPCtrackerMI::MakeSeeds2(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, Float_t */*cuts[4]*/,
3616 Float_t deltay, Bool_t /*bconstrain*/) {
3617 //-----------------------------------------------------------------
3618 // This function creates track seeds - without vertex constraint
3619 //-----------------------------------------------------------------
3620 // cuts[0] - fP4 cut - not applied
3621 // cuts[1] - tan(phi) cut
3622 // cuts[2] - zvertex cut - not applied
3623 // cuts[3] - fP3 cut
3633 Double_t alpha=fOuterSec->GetAlpha(), shift=fOuterSec->GetAlphaShift();
3634 // Double_t cs=cos(alpha), sn=sin(alpha);
3635 Int_t row0 = (i1+i2)/2;
3636 Int_t drow = (i1-i2)/2;
3637 const AliTPCtrackerRow& kr0=fSectors[sec][row0];
3638 AliTPCtrackerRow * kr=0;
3640 AliTPCpolyTrack polytrack;
3641 Int_t nclusters=fSectors[sec][row0];
3642 AliTPCseed * seed = new( NextFreeSeed() ) AliTPCseed;
3643 seed->SetPoolID(fLastSeedID);
3648 for (Int_t is=0; is < nclusters; is++) { //LOOP over clusters
3650 Int_t nfoundable =0;
3651 for (Int_t iter =1; iter<2; iter++){ //iterations
3652 const AliTPCtrackerRow& krm=fSectors[sec][row0-iter];
3653 const AliTPCtrackerRow& krp=fSectors[sec][row0+iter];
3654 const AliTPCclusterMI * cl= kr0[is];
3656 if (cl->IsUsed(10)) {
3662 Double_t x = kr0.GetX();
3663 // Initialization of the polytrack
3668 Double_t y0= cl->GetY();
3669 Double_t z0= cl->GetZ();
3673 Double_t ymax = fSectors->GetMaxY(row0)-kr0.GetDeadZone()-1.5;
3674 if (deltay>0 && TMath::Abs(ymax-TMath::Abs(y0))> deltay ) continue; // seed only at the edge
3676 erry = (0.5)*cl->GetSigmaY2()/TMath::Sqrt(cl->GetQ())*6;
3677 errz = (0.5)*cl->GetSigmaZ2()/TMath::Sqrt(cl->GetQ())*6;
3678 polytrack.AddPoint(x,y0,z0,erry, errz);
3681 if (cl->IsUsed(10)) sumused++;
3684 Float_t roady = (5*TMath::Sqrt(cl->GetSigmaY2()+0.2)+1.)*iter;
3685 Float_t roadz = (5*TMath::Sqrt(cl->GetSigmaZ2()+0.2)+1.)*iter;
3688 AliTPCclusterMI * cl1 = krm.FindNearest(y0,z0,roady,roadz);
3689 if (cl1 && TMath::Abs(ymax-TMath::Abs(y0))) {
3690 erry = (0.5)*cl1->GetSigmaY2()/TMath::Sqrt(cl1->GetQ())*3;
3691 errz = (0.5)*cl1->GetSigmaZ2()/TMath::Sqrt(cl1->GetQ())*3;
3692 if (cl1->IsUsed(10)) sumused++;
3693 polytrack.AddPoint(x,cl1->GetY(),cl1->GetZ(),erry,errz);
3697 AliTPCclusterMI * cl2 = krp.FindNearest(y0,z0,roady,roadz);
3699 erry = (0.5)*cl2->GetSigmaY2()/TMath::Sqrt(cl2->GetQ())*3;
3700 errz = (0.5)*cl2->GetSigmaZ2()/TMath::Sqrt(cl2->GetQ())*3;
3701 if (cl2->IsUsed(10)) sumused++;
3702 polytrack.AddPoint(x,cl2->GetY(),cl2->GetZ(),erry,errz);
3705 if (sumused>0) continue;
3707 polytrack.UpdateParameters();
3713 nfoundable = polytrack.GetN();
3714 nfound = nfoundable;
3716 for (Int_t ddrow = iter+1; ddrow<drow;ddrow++){
3717 Float_t maxdist = 0.8*(1.+3./(ddrow));
3718 for (Int_t delta = -1;delta<=1;delta+=2){
3719 Int_t row = row0+ddrow*delta;
3720 kr = &(fSectors[sec][row]);
3721 Double_t xn = kr->GetX();
3722 Double_t ymax1 = fSectors->GetMaxY(row)-kr->GetDeadZone()-1.5;
3723 polytrack.GetFitPoint(xn,yn,zn);
3724 if (TMath::Abs(yn)>ymax1) continue;
3726 AliTPCclusterMI * cln = kr->FindNearest(yn,zn,roady,roadz);
3728 Float_t dist = TMath::Sqrt( (yn-cln->GetY())*(yn-cln->GetY())+(zn-cln->GetZ())*(zn-cln->GetZ()));
3731 erry = (dist+0.3)*cln->GetSigmaY2()/TMath::Sqrt(cln->GetQ())*(1.+1./(ddrow));
3732 errz = (dist+0.3)*cln->GetSigmaZ2()/TMath::Sqrt(cln->GetQ())*(1.+1./(ddrow));
3733 if (cln->IsUsed(10)) {
3734 // printf("used\n");
3742 polytrack.AddPoint(xn,cln->GetY(),cln->GetZ(),erry, errz);
3747 if ( (sumused>3) || (sumused>0.5*nfound) || (nfound<0.6*nfoundable)) break;
3748 polytrack.UpdateParameters();
3751 if ( (sumused>3) || (sumused>0.5*nfound)) {
3752 //printf("sumused %d\n",sumused);
3757 polytrack.GetFitDerivation(kr0.GetX(),dy,dz);
3758 AliTPCpolyTrack track2;
3760 polytrack.Refit(track2,0.5+TMath::Abs(dy)*0.3,0.4+TMath::Abs(dz)*0.3);
3761 if (track2.GetN()<0.5*nfoundable) continue;
3764 if ((nfound>0.6*nfoundable) &&( nfoundable>0.4*(i1-i2))) {
3766 // test seed with and without constrain
3767 for (Int_t constrain=0; constrain<=0;constrain++){
3768 // add polytrack candidate
3770 Double_t x[5], c[15];
3771 Double_t x1,x2,x3,y1,y2,y3,z1,z2,z3;
3772 track2.GetBoundaries(x3,x1);
3774 track2.GetFitPoint(x1,y1,z1);
3775 track2.GetFitPoint(x2,y2,z2);
3776 track2.GetFitPoint(x3,y3,z3);
3778 //is track pointing to the vertex ?
3781 polytrack.GetFitPoint(x0,y0,z0);
3794 x[4]=F1(x1,y1,x2,y2,x3,y3);
3796 // if (TMath::Abs(x[4]) >= cuts[0]) continue; //
3797 x[2]=F2(x1,y1,x2,y2,x3,y3);
3799 //if (TMath::Abs(x[4]*x1-x[2]) >= cuts[1]) continue;
3800 //x[3]=F3(x1,y1,x2,y2,z1,z2);
3801 x[3]=F3n(x1,y1,x3,y3,z1,z3,x[4]);
3802 //if (TMath::Abs(x[3]) > cuts[3]) continue;
3805 Double_t sy =0.1, sz =0.1;
3806 Double_t sy1=0.02, sz1=0.02;
3807 Double_t sy2=0.02, sz2=0.02;
3811 sy3=25000*x[4]*x[4]+0.1, sy=0.1, sz=0.1;
3814 Double_t f40=(F1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
3815 Double_t f42=(F1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
3816 Double_t f43=(F1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
3817 Double_t f20=(F2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
3818 Double_t f22=(F2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
3819 Double_t f23=(F2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
3821 Double_t f30=(F3(x1,y1+sy,x3,y3,z1,z3)-x[3])/sy;
3822 Double_t f31=(F3(x1,y1,x3,y3,z1+sz,z3)-x[3])/sz;
3823 Double_t f32=(F3(x1,y1,x3,y3+sy,z1,z3)-x[3])/sy;
3824 Double_t f34=(F3(x1,y1,x3,y3,z1,z3+sz)-x[3])/sz;
3829 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3830 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3831 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3832 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3833 c[13]=f30*sy1*f40+f32*sy2*f42;
3834 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3836 //Int_t row1 = fSectors->GetRowNumber(x1);
3837 Int_t row1 = GetRowNumber(x1);
3841 if (seed) {MarkSeedFree( seed ); seed = 0;}
3842 AliTPCseed *track = seed = new( NextFreeSeed() ) AliTPCseed(x1,sec*alpha+shift,x,c,index);
3843 seed->SetPoolID(fLastSeedID);
3844 track->SetIsSeeding(kTRUE);
3845 Int_t rc=FollowProlongation(*track, i2);
3846 if (constrain) track->SetBConstrain(1);
3848 track->SetBConstrain(0);
3849 track->SetLastPoint(row1+fInnerSec->GetNRows()); // first cluster in track position
3850 track->SetFirstPoint(track->GetLastPoint());
3852 if (rc==0 || track->GetNumberOfClusters()<(i1-i2)*0.5 ||
3853 track->GetNumberOfClusters() < track->GetNFoundable()*0.6 ||
3854 track->GetNShared()>0.4*track->GetNumberOfClusters()) {
3855 MarkSeedFree( seed ); seed = 0;
3858 arr->AddLast(track); // track IS seed, don't free seed
3859 seed = new( NextFreeSeed() ) AliTPCseed;
3860 seed->SetPoolID(fLastSeedID);
3864 } // if accepted seed
3867 Info("MakeSeeds2","\nSeeding statiistic:\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin3);
3869 if (seed) MarkSeedFree( seed );
3873 AliTPCseed *AliTPCtrackerMI::MakeSeed(AliTPCseed *const track, Float_t r0, Float_t r1, Float_t r2)
3877 //reseed using track points
3878 Int_t p0 = int(r0*track->GetNumberOfClusters()); // point 0
3879 Int_t p1 = int(r1*track->GetNumberOfClusters());
3880 Int_t p2 = int(r2*track->GetNumberOfClusters()); // last point
3882 Double_t x0[3],x1[3],x2[3];
3883 for (Int_t i=0;i<3;i++){
3889 // find track position at given ratio of the length
3890 Int_t sec0=0, sec1=0, sec2=0;
3893 for (Int_t i=0;i<160;i++){
3894 if (track->GetClusterPointer(i)){
3896 AliTPCTrackerPoint *trpoint =track->GetTrackPoint(i);
3897 if ( (index<p0) || x0[0]<0 ){
3898 if (trpoint->GetX()>1){
3899 clindex = track->GetClusterIndex2(i);
3901 x0[0] = trpoint->GetX();
3902 x0[1] = trpoint->GetY();
3903 x0[2] = trpoint->GetZ();
3904 sec0 = ((clindex&0xff000000)>>24)%18;
3909 if ( (index<p1) &&(trpoint->GetX()>1)){
3910 clindex = track->GetClusterIndex2(i);
3912 x1[0] = trpoint->GetX();
3913 x1[1] = trpoint->GetY();
3914 x1[2] = trpoint->GetZ();
3915 sec1 = ((clindex&0xff000000)>>24)%18;
3918 if ( (index<p2) &&(trpoint->GetX()>1)){
3919 clindex = track->GetClusterIndex2(i);
3921 x2[0] = trpoint->GetX();
3922 x2[1] = trpoint->GetY();
3923 x2[2] = trpoint->GetZ();
3924 sec2 = ((clindex&0xff000000)>>24)%18;
3931 Double_t alpha, cs,sn, xx2,yy2;
3933 alpha = (sec1-sec2)*fSectors->GetAlpha();
3934 cs = TMath::Cos(alpha);
3935 sn = TMath::Sin(alpha);
3936 xx2= x1[0]*cs-x1[1]*sn;
3937 yy2= x1[0]*sn+x1[1]*cs;
3941 alpha = (sec0-sec2)*fSectors->GetAlpha();
3942 cs = TMath::Cos(alpha);
3943 sn = TMath::Sin(alpha);
3944 xx2= x0[0]*cs-x0[1]*sn;
3945 yy2= x0[0]*sn+x0[1]*cs;
3951 Double_t x[5],c[15];
3955 x[4]=F1(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]);
3956 // if (x[4]>1) return 0;
3957 x[2]=F2(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]);
3958 x[3]=F3n(x2[0],x2[1],x0[0],x0[1],x2[2],x0[2],x[4]);
3959 //if (TMath::Abs(x[3]) > 2.2) return 0;
3960 //if (TMath::Abs(x[2]) > 1.99) return 0;
3962 Double_t sy =0.1, sz =0.1;
3964 Double_t sy1=0.02+track->GetSigmaY2(), sz1=0.02+track->GetSigmaZ2();
3965 Double_t sy2=0.01+track->GetSigmaY2(), sz2=0.01+track->GetSigmaZ2();
3966 Double_t sy3=0.01+track->GetSigmaY2();
3968 Double_t f40=(F1(x2[0],x2[1]+sy,x1[0],x1[1],x0[0],x0[1])-x[4])/sy;
3969 Double_t f42=(F1(x2[0],x2[1],x1[0],x1[1]+sy,x0[0],x0[1])-x[4])/sy;
3970 Double_t f43=(F1(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]+sy)-x[4])/sy;
3971 Double_t f20=(F2(x2[0],x2[1]+sy,x1[0],x1[1],x0[0],x0[1])-x[2])/sy;
3972 Double_t f22=(F2(x2[0],x2[1],x1[0],x1[1]+sy,x0[0],x0[1])-x[2])/sy;
3973 Double_t f23=(F2(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]+sy)-x[2])/sy;
3975 Double_t f30=(F3(x2[0],x2[1]+sy,x0[0],x0[1],x2[2],x0[2])-x[3])/sy;
3976 Double_t f31=(F3(x2[0],x2[1],x0[0],x0[1],x2[2]+sz,x0[2])-x[3])/sz;
3977 Double_t f32=(F3(x2[0],x2[1],x0[0],x0[1]+sy,x2[2],x0[2])-x[3])/sy;
3978 Double_t f34=(F3(x2[0],x2[1],x0[0],x0[1],x2[2],x0[2]+sz)-x[3])/sz;
3983 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3984 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3985 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3986 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3987 c[13]=f30*sy1*f40+f32*sy2*f42;
3988 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3990 // Int_t row1 = fSectors->GetRowNumber(x2[0]);
3991 AliTPCseed *seed = new( NextFreeSeed() ) AliTPCseed(x2[0], sec2*fSectors->GetAlpha()+fSectors->GetAlphaShift(), x, c, 0);
3992 seed->SetPoolID(fLastSeedID);
3993 // Double_t y0,z0,y1,z1, y2,z2;
3994 //seed->GetProlongation(x0[0],y0,z0);
3995 // seed->GetProlongation(x1[0],y1,z1);
3996 //seed->GetProlongation(x2[0],y2,z2);
3998 seed->SetLastPoint(pp2);
3999 seed->SetFirstPoint(pp2);
4006 AliTPCseed *AliTPCtrackerMI::ReSeed(const AliTPCseed *track, Float_t r0, Float_t r1, Float_t r2)
4010 //reseed using founded clusters
4012 // Find the number of clusters
4013 Int_t nclusters = 0;
4014 for (Int_t irow=0;irow<160;irow++){
4015 if (track->GetClusterIndex(irow)>0) nclusters++;
4019 ipos[0] = TMath::Max(int(r0*nclusters),0); // point 0 cluster
4020 ipos[1] = TMath::Min(int(r1*nclusters),nclusters-1); //
4021 ipos[2] = TMath::Min(int(r2*nclusters),nclusters-1); // last point
4024 Double_t xyz[3][3]={{0}};
4025 Int_t row[3]={0},sec[3]={0,0,0};
4027 // find track row position at given ratio of the length
4029 for (Int_t irow=0;irow<160;irow++){
4030 if (track->GetClusterIndex2(irow)<0) continue;
4032 for (Int_t ipoint=0;ipoint<3;ipoint++){
4033 if (index<=ipos[ipoint]) row[ipoint] = irow;
4037 //Get cluster and sector position
4038 for (Int_t ipoint=0;ipoint<3;ipoint++){
4039 Int_t clindex = track->GetClusterIndex2(row[ipoint]);
4040 AliTPCclusterMI * cl = GetClusterMI(clindex);
4043 // AliTPCclusterMI * cl = GetClusterMI(clindex);
4046 sec[ipoint] = ((clindex&0xff000000)>>24)%18;
4047 xyz[ipoint][0] = GetXrow(row[ipoint]);
4048 xyz[ipoint][1] = cl->GetY();
4049 xyz[ipoint][2] = cl->GetZ();
4053 // Calculate seed state vector and covariance matrix
4055 Double_t alpha, cs,sn, xx2,yy2;
4057 alpha = (sec[1]-sec[2])*fSectors->GetAlpha();
4058 cs = TMath::Cos(alpha);
4059 sn = TMath::Sin(alpha);
4060 xx2= xyz[1][0]*cs-xyz[1][1]*sn;
4061 yy2= xyz[1][0]*sn+xyz[1][1]*cs;
4065 alpha = (sec[0]-sec[2])*fSectors->GetAlpha();
4066 cs = TMath::Cos(alpha);
4067 sn = TMath::Sin(alpha);
4068 xx2= xyz[0][0]*cs-xyz[0][1]*sn;
4069 yy2= xyz[0][0]*sn+xyz[0][1]*cs;
4075 Double_t x[5],c[15];
4079 x[4]=F1(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
4080 x[2]=F2(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
4081 x[3]=F3n(xyz[2][0],xyz[2][1],xyz[0][0],xyz[0][1],xyz[2][2],xyz[0][2],x[4]);
4083 Double_t sy =0.1, sz =0.1;
4085 Double_t sy1=0.2, sz1=0.2;
4086 Double_t sy2=0.2, sz2=0.2;
4089 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;
4090 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;
4091 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;
4092 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;
4093 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;
4094 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;
4096 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;
4097 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;
4098 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;
4099 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;
4104 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
4105 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
4106 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
4107 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
4108 c[13]=f30*sy1*f40+f32*sy2*f42;
4109 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
4111 // Int_t row1 = fSectors->GetRowNumber(xyz[2][0]);
4112 AliTPCseed *seed=new( NextFreeSeed() ) AliTPCseed(xyz[2][0], sec[2]*fSectors->GetAlpha()+fSectors->GetAlphaShift(), x, c, 0);
4113 seed->SetPoolID(fLastSeedID);
4114 seed->SetLastPoint(row[2]);
4115 seed->SetFirstPoint(row[2]);
4120 AliTPCseed *AliTPCtrackerMI::ReSeed(AliTPCseed *track,Int_t r0, Bool_t forward)
4124 //reseed using founded clusters
4127 Int_t row[3]={0,0,0};
4128 Int_t sec[3]={0,0,0};
4130 // forward direction
4132 for (Int_t irow=r0;irow<160;irow++){
4133 if (track->GetClusterIndex(irow)>0){
4138 for (Int_t irow=160;irow>r0;irow--){
4139 if (track->GetClusterIndex(irow)>0){
4144 for (Int_t irow=row[2]-15;irow>row[0];irow--){
4145 if (track->GetClusterIndex(irow)>0){
4153 for (Int_t irow=0;irow<r0;irow++){
4154 if (track->GetClusterIndex(irow)>0){
4159 for (Int_t irow=r0;irow>0;irow--){
4160 if (track->GetClusterIndex(irow)>0){
4165 for (Int_t irow=row[2]-15;irow>row[0];irow--){
4166 if (track->GetClusterIndex(irow)>0){
4173 if ((row[2]-row[0])<20) return 0;
4174 if (row[1]==0) return 0;
4177 //Get cluster and sector position
4178 for (Int_t ipoint=0;ipoint<3;ipoint++){
4179 Int_t clindex = track->GetClusterIndex2(row[ipoint]);
4180 AliTPCclusterMI * cl = GetClusterMI(clindex);
4183 // AliTPCclusterMI * cl = GetClusterMI(clindex);
4186 sec[ipoint] = ((clindex&0xff000000)>>24)%18;
4187 xyz[ipoint][0] = GetXrow(row[ipoint]);
4188 AliTPCTrackerPoint * point = track->GetTrackPoint(row[ipoint]);
4189 if (point&&ipoint<2){
4191 xyz[ipoint][1] = point->GetY();
4192 xyz[ipoint][2] = point->GetZ();
4195 xyz[ipoint][1] = cl->GetY();
4196 xyz[ipoint][2] = cl->GetZ();
4203 // Calculate seed state vector and covariance matrix
4205 Double_t alpha, cs,sn, xx2,yy2;
4207 alpha = (sec[1]-sec[2])*fSectors->GetAlpha();
4208 cs = TMath::Cos(alpha);
4209 sn = TMath::Sin(alpha);
4210 xx2= xyz[1][0]*cs-xyz[1][1]*sn;
4211 yy2= xyz[1][0]*sn+xyz[1][1]*cs;
4215 alpha = (sec[0]-sec[2])*fSectors->GetAlpha();
4216 cs = TMath::Cos(alpha);
4217 sn = TMath::Sin(alpha);
4218 xx2= xyz[0][0]*cs-xyz[0][1]*sn;
4219 yy2= xyz[0][0]*sn+xyz[0][1]*cs;
4225 Double_t x[5],c[15];
4229 x[4]=F1(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
4230 x[2]=F2(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
4231 x[3]=F3n(xyz[2][0],xyz[2][1],xyz[0][0],xyz[0][1],xyz[2][2],xyz[0][2],x[4]);
4233 Double_t sy =0.1, sz =0.1;
4235 Double_t sy1=0.2, sz1=0.2;
4236 Double_t sy2=0.2, sz2=0.2;
4239 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;
4240 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;
4241 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;
4242 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;
4243 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;
4244 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;
4246 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;
4247 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;
4248 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;
4249 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;
4254 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
4255 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
4256 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
4257 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
4258 c[13]=f30*sy1*f40+f32*sy2*f42;
4259 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
4261 // Int_t row1 = fSectors->GetRowNumber(xyz[2][0]);
4262 AliTPCseed *seed=new( NextFreeSeed() ) AliTPCseed(xyz[2][0], sec[2]*fSectors->GetAlpha()+fSectors->GetAlphaShift(), x, c, 0);
4263 seed->SetPoolID(fLastSeedID);
4264 seed->SetLastPoint(row[2]);
4265 seed->SetFirstPoint(row[2]);
4266 for (Int_t i=row[0];i<row[2];i++){
4267 seed->SetClusterIndex(i, track->GetClusterIndex(i));
4275 void AliTPCtrackerMI::FindMultiMC(const TObjArray * array, AliESDEvent */*esd*/, Int_t iter)
4278 // find multi tracks - THIS FUNCTION IS ONLY FOR DEBUG PURPOSES
4280 // Use AliTPCReconstructor::StreamLevel()>2 if you want to tune parameters - cuts
4282 // Two reasons to have multiple find tracks
4283 // 1. Curling tracks can be find more than once
4284 // 2. Splitted tracks
4285 // a.) Multiple seeding to increase tracking efficiency - (~ 100% reached)
4286 // b.) Edge effect on the sector boundaries
4289 // Algorithm done in 2 phases - because of CPU consumption
4290 // it is n^2 algorithm - for lead-lead 20000x20000 combination are investigated
4292 // Algorihm for curling tracks sign:
4293 // 1 phase -makes a very rough fast cuts to minimize combinatorics
4294 // a.) opposite sign
4295 // b.) one of the tracks - not pointing to the primary vertex -
4296 // c.) delta tan(theta)
4298 // 2 phase - calculates DCA between tracks - time consument
4303 // General cuts - for splitted tracks and for curling tracks
4305 const Float_t kMaxdPhi = 0.2; // maximal distance in phi
4307 // Curling tracks cuts
4312 Int_t nentries = array->GetEntriesFast();
4313 AliHelix *helixes = new AliHelix[nentries];
4314 Float_t *xm = new Float_t[nentries];
4315 Float_t *dz0 = new Float_t[nentries];
4316 Float_t *dz1 = new Float_t[nentries];
4322 // Find track COG in x direction - point with best defined parameters
4324 for (Int_t i=0;i<nentries;i++){
4325 AliTPCseed* track = (AliTPCseed*)array->At(i);
4326 if (!track) continue;
4327 track->SetCircular(0);
4328 new (&helixes[i]) AliHelix(*track);
4332 track->GetDZ(GetX(),GetY(),GetZ(),GetBz(),dz);
4335 for (Int_t icl=0; icl<160; icl++){
4336 AliTPCclusterMI * cl = track->GetClusterPointer(icl);
4342 if (ncl>0) xm[i]/=Float_t(ncl);
4345 for (Int_t i0=0;i0<nentries;i0++){
4346 AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
4347 if (!track0) continue;
4348 Float_t xc0 = helixes[i0].GetHelix(6);
4349 Float_t yc0 = helixes[i0].GetHelix(7);
4350 Float_t r0 = helixes[i0].GetHelix(8);
4351 Float_t rc0 = TMath::Sqrt(xc0*xc0+yc0*yc0);
4352 Float_t fi0 = TMath::ATan2(yc0,xc0);
4354 for (Int_t i1=i0+1;i1<nentries;i1++){
4355 AliTPCseed * track1 = (AliTPCseed*)array->At(i1);
4356 if (!track1) continue;
4357 Int_t lab0=track0->GetLabel();
4358 Int_t lab1=track1->GetLabel();
4359 if (TMath::Abs(lab0)!=TMath::Abs(lab1)) continue;
4361 Float_t xc1 = helixes[i1].GetHelix(6);
4362 Float_t yc1 = helixes[i1].GetHelix(7);
4363 Float_t r1 = helixes[i1].GetHelix(8);
4364 Float_t rc1 = TMath::Sqrt(xc1*xc1+yc1*yc1);
4365 Float_t fi1 = TMath::ATan2(yc1,xc1);
4367 Float_t dfi = fi0-fi1;
4370 if (dfi>1.5*TMath::Pi()) dfi-=TMath::Pi(); // take care about edge effect
4371 if (dfi<-1.5*TMath::Pi()) dfi+=TMath::Pi(); //
4372 if (TMath::Abs(dfi)>kMaxdPhi&&helixes[i0].GetHelix(4)*helixes[i1].GetHelix(4)<0){
4374 // if short tracks with undefined sign
4375 fi1 = -TMath::ATan2(yc1,-xc1);
4378 Float_t dtheta = TMath::Abs(track0->GetTgl()-track1->GetTgl())<TMath::Abs(track0->GetTgl()+track1->GetTgl())? track0->GetTgl()-track1->GetTgl():track0->GetTgl()+track1->GetTgl();
4381 // debug stream to tune "fast cuts"
4383 Double_t dist[3]; // distance at X
4384 Double_t mdist[3]={0,0,0}; // mean distance X+-40cm
4385 track0->GetDistance(track1,0.5*(xm[i0]+xm[i1])-40.,dist,AliTracker::GetBz());
4386 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[i]);
4387 track0->GetDistance(track1,0.5*(xm[i0]+xm[i1])+40.,dist,AliTracker::GetBz());
4388 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[i]);
4389 track0->GetDistance(track1,0.5*(xm[i0]+xm[i1]),dist,AliTracker::GetBz());
4390 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[i]);
4391 for (Int_t i=0;i<3;i++) mdist[i]*=0.33333;
4395 for (Int_t icl=0; icl<160; icl++){
4396 AliTPCclusterMI * cl0 = track0->GetClusterPointer(icl);
4397 AliTPCclusterMI * cl1 = track1->GetClusterPointer(icl);
4400 if (cl0==cl1) sums++;
4404 if (AliTPCReconstructor::StreamLevel()>5) {
4405 TTreeSRedirector &cstream = *fDebugStreamer;
4410 "Tr0.="<<track0<< // seed0
4411 "Tr1.="<<track1<< // seed1
4412 "h0.="<<&helixes[i0]<<
4413 "h1.="<<&helixes[i1]<<
4415 "sum="<<sum<< //the sum of rows with cl in both
4416 "sums="<<sums<< //the sum of shared clusters
4417 "xm0="<<xm[i0]<< // the center of track
4418 "xm1="<<xm[i1]<< // the x center of track
4419 // General cut variables
4420 "dfi="<<dfi<< // distance in fi angle
4421 "dtheta="<<dtheta<< // distance int theta angle
4427 "dist0="<<dist[0]<< //distance x
4428 "dist1="<<dist[1]<< //distance y
4429 "dist2="<<dist[2]<< //distance z
4430 "mdist0="<<mdist[0]<< //distance x
4431 "mdist1="<<mdist[1]<< //distance y
4432 "mdist2="<<mdist[2]<< //distance z
4448 if (AliTPCReconstructor::StreamLevel()>1) {
4449 AliInfo("Time for curling tracks removal DEBUGGING MC");
4456 void AliTPCtrackerMI::FindSplitted(TObjArray * array, AliESDEvent */*esd*/, Int_t /*iter*/){
4458 // Find Splitted tracks and remove the one with worst quality
4459 // Corresponding debug streamer to tune selections - "Splitted2"
4461 // 0. Sort tracks according quility
4462 // 1. Propagate the tracks to the reference radius
4463 // 2. Double_t loop to select close tracks (only to speed up process)
4464 // 3. Calculate cluster overlap ratio - and remove the track if bigger than a threshold
4465 // 4. Delete temporary parameters
4467 const Double_t xref=GetXrow(63); // reference radius -IROC/OROC boundary
4469 const Double_t kCutP1=10; // delta Z cut 10 cm
4470 const Double_t kCutP2=0.15; // delta snp(fi) cut 0.15
4471 const Double_t kCutP3=0.15; // delta tgl(theta) cut 0.15
4472 const Double_t kCutAlpha=0.15; // delta alpha cut
4473 Int_t firstpoint = 0;
4474 Int_t lastpoint = 160;
4476 Int_t nentries = array->GetEntriesFast();
4477 AliExternalTrackParam *params = new AliExternalTrackParam[nentries];
4483 //0. Sort tracks according quality
4484 //1. Propagate the ext. param to reference radius
4485 Int_t nseed = array->GetEntriesFast();
4486 if (nseed<=0) return;
4487 Float_t * quality = new Float_t[nseed];
4488 Int_t * indexes = new Int_t[nseed];
4489 for (Int_t i=0; i<nseed; i++) {
4490 AliTPCseed *pt=(AliTPCseed*)array->UncheckedAt(i);
4495 pt->UpdatePoints(); //select first last max dens points
4496 Float_t * points = pt->GetPoints();
4497 if (points[3]<0.8) quality[i] =-1;
4498 quality[i] = (points[2]-points[0])+pt->GetNumberOfClusters();
4499 //prefer high momenta tracks if overlaps
4500 quality[i] *= TMath::Sqrt(TMath::Abs(pt->Pt())+0.5);
4502 AliTracker::PropagateTrackToBxByBz(&(params[i]),xref,pt->GetMass(),5.,kTRUE);
4503 AliTracker::PropagateTrackToBxByBz(&(params[i]),xref,pt->GetMass(),1.,kTRUE);
4505 TMath::Sort(nseed,quality,indexes);
4507 // 3. Loop over pair of tracks
4509 for (Int_t i0=0; i0<nseed; i0++) {
4510 Int_t index0=indexes[i0];
4511 if (!(array->UncheckedAt(index0))) continue;
4512 AliTPCseed *s1 = (AliTPCseed*)array->UncheckedAt(index0);
4513 if (!s1->IsActive()) continue;
4514 AliExternalTrackParam &par0=params[index0];
4515 for (Int_t i1=i0+1; i1<nseed; i1++) {
4516 Int_t index1=indexes[i1];
4517 if (!(array->UncheckedAt(index1))) continue;
4518 AliTPCseed *s2 = (AliTPCseed*)array->UncheckedAt(index1);
4519 if (!s2->IsActive()) continue;
4520 if (s2->GetKinkIndexes()[0]!=0)
4521 if (s2->GetKinkIndexes()[0] == -s1->GetKinkIndexes()[0]) continue;
4522 AliExternalTrackParam &par1=params[index1];
4523 if (TMath::Abs(par0.GetParameter()[3]-par1.GetParameter()[3])>kCutP3) continue;
4524 if (TMath::Abs(par0.GetParameter()[1]-par1.GetParameter()[1])>kCutP1) continue;
4525 if (TMath::Abs(par0.GetParameter()[2]-par1.GetParameter()[2])>kCutP2) continue;
4526 Double_t dAlpha= TMath::Abs(par0.GetAlpha()-par1.GetAlpha());
4527 if (dAlpha>TMath::Pi()) dAlpha-=TMath::Pi();
4528 if (TMath::Abs(dAlpha)>kCutAlpha) continue;
4533 Int_t firstShared=lastpoint, lastShared=firstpoint;
4534 Int_t firstRow=lastpoint, lastRow=firstpoint;
4536 for (Int_t i=firstpoint;i<lastpoint;i++){
4537 if (s1->GetClusterIndex2(i)>0) nall0++;
4538 if (s2->GetClusterIndex2(i)>0) nall1++;
4539 if (s1->GetClusterIndex2(i)>0 && s2->GetClusterIndex2(i)>0) {
4540 if (i<firstRow) firstRow=i;
4541 if (i>lastRow) lastRow=i;
4543 if ( (s1->GetClusterIndex2(i))==(s2->GetClusterIndex2(i)) && s1->GetClusterIndex2(i)>0) {
4544 if (i<firstShared) firstShared=i;
4545 if (i>lastShared) lastShared=i;
4549 Double_t ratio0 = Float_t(sumShared)/Float_t(TMath::Min(nall0+1,nall1+1));
4550 Double_t ratio1 = Float_t(sumShared)/Float_t(TMath::Max(nall0+1,nall1+1));
4552 if( AliTPCReconstructor::StreamLevel()>1){
4553 TTreeSRedirector &cstream = *fDebugStreamer;
4554 Int_t n0=s1->GetNumberOfClusters();
4555 Int_t n1=s2->GetNumberOfClusters();
4556 Int_t n0F=s1->GetNFoundable();
4557 Int_t n1F=s2->GetNFoundable();
4558 Int_t lab0=s1->GetLabel();
4559 Int_t lab1=s2->GetLabel();
4561 cstream<<"Splitted2"<<
4562 "iter="<<fIteration<<
4563 "lab0="<<lab0<< // MC label if exist
4564 "lab1="<<lab1<< // MC label if exist
4567 "ratio0="<<ratio0<< // shared ratio
4568 "ratio1="<<ratio1<< // shared ratio
4569 "p0.="<<&par0<< // track parameters
4571 "s0.="<<s1<< // full seed
4573 "n0="<<n0<< // number of clusters track 0
4574 "n1="<<n1<< // number of clusters track 1
4575 "nall0="<<nall0<< // number of clusters track 0
4576 "nall1="<<nall1<< // number of clusters track 1
4577 "n0F="<<n0F<< // number of findable
4578 "n1F="<<n1F<< // number of findable
4579 "shared="<<sumShared<< // number of shared clusters
4580 "firstS="<<firstShared<< // first and the last shared row
4581 "lastS="<<lastShared<<
4582 "firstRow="<<firstRow<< // first and the last row with cluster
4583 "lastRow="<<lastRow<< //
4587 // remove track with lower quality
4589 if (ratio0>AliTPCReconstructor::GetRecoParam()->GetCutSharedClusters(0) ||
4590 ratio1>AliTPCReconstructor::GetRecoParam()->GetCutSharedClusters(1)){
4594 MarkSeedFree( array->RemoveAt(index1) );
4599 // 4. Delete temporary array
4609 void AliTPCtrackerMI::FindCurling(const TObjArray * array, AliESDEvent */*esd*/, Int_t iter)
4612 // find Curling tracks
4613 // Use AliTPCReconstructor::StreamLevel()>1 if you want to tune parameters - cuts
4616 // Algorithm done in 2 phases - because of CPU consumption
4617 // it is n^2 algorithm - for lead-lead 20000x20000 combination are investigated
4618 // see detal in MC part what can be used to cut
4622 const Float_t kMaxC = 400; // maximal curvature to of the track
4623 const Float_t kMaxdTheta = 0.15; // maximal distance in theta
4624 const Float_t kMaxdPhi = 0.15; // maximal distance in phi
4625 const Float_t kPtRatio = 0.3; // ratio between pt
4626 const Float_t kMinDCAR = 2.; // distance to the primary vertex in r - see cpipe cut
4629 // Curling tracks cuts
4632 const Float_t kMaxDeltaRMax = 40; // distance in outer radius
4633 const Float_t kMaxDeltaRMin = 5.; // distance in lower radius - see cpipe cut
4634 const Float_t kMinAngle = 2.9; // angle between tracks
4635 const Float_t kMaxDist = 5; // biggest distance
4637 // The cuts can be tuned using the "MC information stored in Multi tree ==> see FindMultiMC
4640 TCut csign("csign","Tr0.fP[4]*Tr1.fP[4]<0"); //opposite sign
4641 TCut cmax("cmax","abs(Tr0.GetC())>1/400");
4642 TCut cda("cda","sqrt(dtheta^2+dfi^2)<0.15");
4643 TCut ccratio("ccratio","abs((Tr0.fP[4]+Tr1.fP[4])/(abs(Tr0.fP[4])+abs(Tr1.fP[4])))<0.3");
4644 TCut cpipe("cpipe", "min(abs(r0-rc0),abs(r1-rc1))>5");
4646 TCut cdrmax("cdrmax","abs(abs(rc0+r0)-abs(rc1+r1))<40")
4647 TCut cdrmin("cdrmin","abs(abs(rc0+r0)-abs(rc1+r1))<10")
4649 Multi->Draw("dfi","iter==0"+csign+cmax+cda+ccratio); ~94% of curling tracks fulfill
4650 Multi->Draw("min(abs(r0-rc0),abs(r1-rc1))","iter==0&&abs(lab1)==abs(lab0)"+csign+cmax+cda+ccratio+cpipe+cdrmin+cdrmax); //80%
4652 Curling2->Draw("dfi","iter==0&&abs(lab0)==abs(lab1)"+csign+cmax+cdtheta+cdfi+ccratio)
4658 Int_t nentries = array->GetEntriesFast();
4659 AliHelix *helixes = new AliHelix[nentries];
4660 for (Int_t i=0;i<nentries;i++){
4661 AliTPCseed* track = (AliTPCseed*)array->At(i);
4662 if (!track) continue;
4663 track->SetCircular(0);
4664 new (&helixes[i]) AliHelix(*track);
4670 Double_t phase[2][2]={{0,0},{0,0}},radius[2]={0,0};
4676 for (Int_t i0=0;i0<nentries;i0++){
4677 AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
4678 if (!track0) continue;
4679 if (TMath::Abs(track0->GetC())<1/kMaxC) continue;
4680 Float_t xc0 = helixes[i0].GetHelix(6);
4681 Float_t yc0 = helixes[i0].GetHelix(7);
4682 Float_t r0 = helixes[i0].GetHelix(8);
4683 Float_t rc0 = TMath::Sqrt(xc0*xc0+yc0*yc0);
4684 Float_t fi0 = TMath::ATan2(yc0,xc0);
4686 for (Int_t i1=i0+1;i1<nentries;i1++){
4687 AliTPCseed * track1 = (AliTPCseed*)array->At(i1);
4688 if (!track1) continue;
4689 if (TMath::Abs(track1->GetC())<1/kMaxC) continue;
4690 Float_t xc1 = helixes[i1].GetHelix(6);
4691 Float_t yc1 = helixes[i1].GetHelix(7);
4692 Float_t r1 = helixes[i1].GetHelix(8);
4693 Float_t rc1 = TMath::Sqrt(xc1*xc1+yc1*yc1);
4694 Float_t fi1 = TMath::ATan2(yc1,xc1);
4696 Float_t dfi = fi0-fi1;
4699 if (dfi>1.5*TMath::Pi()) dfi-=TMath::Pi(); // take care about edge effect
4700 if (dfi<-1.5*TMath::Pi()) dfi+=TMath::Pi(); //
4701 Float_t dtheta = TMath::Abs(track0->GetTgl()-track1->GetTgl())<TMath::Abs(track0->GetTgl()+track1->GetTgl())? track0->GetTgl()-track1->GetTgl():track0->GetTgl()+track1->GetTgl();
4705 if (track0->GetBConstrain()&&track1->GetBConstrain()) continue; // not constrained
4706 if (track1->GetSigned1Pt()*track0->GetSigned1Pt()>0) continue; // not the same sign
4707 if ( TMath::Abs(track1->GetTgl()+track0->GetTgl())>kMaxdTheta) continue; //distance in the Theta
4708 if ( TMath::Abs(dfi)>kMaxdPhi) continue; //distance in phi
4709 if ( TMath::Sqrt(dfi*dfi+dtheta*dtheta)>kMaxdPhi) continue; //common angular offset
4711 Float_t pt0 = track0->GetSignedPt();
4712 Float_t pt1 = track1->GetSignedPt();
4713 if ((TMath::Abs(pt0+pt1)/(TMath::Abs(pt0)+TMath::Abs(pt1)))>kPtRatio) continue;
4714 if ((iter==1) && TMath::Abs(TMath::Abs(rc0+r0)-TMath::Abs(rc1+r1))>kMaxDeltaRMax) continue;
4715 if ((iter!=1) &&TMath::Abs(TMath::Abs(rc0-r0)-TMath::Abs(rc1-r1))>kMaxDeltaRMin) continue;
4716 if (TMath::Min(TMath::Abs(rc0-r0),TMath::Abs(rc1-r1))<kMinDCAR) continue;
4719 // Now find closest approach
4723 Int_t npoints = helixes[i0].GetRPHIintersections(helixes[i1], phase, radius,10);
4724 if (npoints==0) continue;
4725 helixes[i0].GetClosestPhases(helixes[i1], phase);
4729 Double_t hangles[3];
4730 helixes[i0].Evaluate(phase[0][0],xyz0);
4731 helixes[i1].Evaluate(phase[0][1],xyz1);
4733 helixes[i0].GetAngle(phase[0][0],helixes[i1],phase[0][1],hangles);
4734 Double_t deltah[2],deltabest;
4735 if (TMath::Abs(hangles[2])<kMinAngle) continue;
4739 helixes[i0].ParabolicDCA(helixes[i1],phase[0][0],phase[0][1],radius[0],deltah[0],2);
4741 helixes[i0].ParabolicDCA(helixes[i1],phase[1][0],phase[1][1],radius[1],deltah[1],2);
4742 if (deltah[1]<deltah[0]) ibest=1;
4744 deltabest = TMath::Sqrt(deltah[ibest]);
4745 helixes[i0].Evaluate(phase[ibest][0],xyz0);
4746 helixes[i1].Evaluate(phase[ibest][1],xyz1);
4747 helixes[i0].GetAngle(phase[ibest][0],helixes[i1],phase[ibest][1],hangles);
4748 Double_t radiusbest = TMath::Sqrt(radius[ibest]);
4750 if (deltabest>kMaxDist) continue;
4751 // if (mindcar+mindcaz<40 && (TMath::Abs(hangles[2])<kMinAngle ||deltabest>3)) continue;
4752 Bool_t sign =kFALSE;
4753 if (hangles[2]>kMinAngle) sign =kTRUE;
4756 // circular[i0] = kTRUE;
4757 // circular[i1] = kTRUE;
4758 if (track0->OneOverPt()<track1->OneOverPt()){
4759 track0->SetCircular(track0->GetCircular()+1);
4760 track1->SetCircular(track1->GetCircular()+2);
4763 track1->SetCircular(track1->GetCircular()+1);
4764 track0->SetCircular(track0->GetCircular()+2);
4767 if (AliTPCReconstructor::StreamLevel()>2){
4769 //debug stream to tune "fine" cuts
4770 Int_t lab0=track0->GetLabel();
4771 Int_t lab1=track1->GetLabel();
4772 TTreeSRedirector &cstream = *fDebugStreamer;
4773 cstream<<"Curling2"<<
4789 "npoints="<<npoints<<
4790 "hangles0="<<hangles[0]<<
4791 "hangles1="<<hangles[1]<<
4792 "hangles2="<<hangles[2]<<
4795 "radius="<<radiusbest<<
4796 "deltabest="<<deltabest<<
4797 "phase0="<<phase[ibest][0]<<
4798 "phase1="<<phase[ibest][1]<<
4806 if (AliTPCReconstructor::StreamLevel()>1) {
4807 AliInfo("Time for curling tracks removal");
4813 void AliTPCtrackerMI::FindKinks(TObjArray * array, AliESDEvent *esd)
4819 // RS something is wrong in this routine: not all seeds are assigned to daughters and mothers array, but they all are queried
4822 TObjArray *kinks= new TObjArray(10000);
4823 // TObjArray *v0s= new TObjArray(10000);
4824 Int_t nentries = array->GetEntriesFast();
4825 AliHelix *helixes = new AliHelix[nentries];
4826 Int_t *sign = new Int_t[nentries];
4827 Int_t *nclusters = new Int_t[nentries];
4828 Float_t *alpha = new Float_t[nentries];
4829 AliKink *kink = new AliKink();
4830 Int_t * usage = new Int_t[nentries];
4831 Float_t *zm = new Float_t[nentries];
4832 Float_t *z0 = new Float_t[nentries];
4833 Float_t *fim = new Float_t[nentries];
4834 Float_t *shared = new Float_t[nentries];
4835 Bool_t *circular = new Bool_t[nentries];
4836 Float_t *dca = new Float_t[nentries];
4837 //const AliESDVertex * primvertex = esd->GetVertex();
4839 // nentries = array->GetEntriesFast();
4844 for (Int_t i=0;i<nentries;i++){
4847 AliTPCseed* track = (AliTPCseed*)array->At(i);
4848 if (!track) continue;
4849 track->SetCircular(0);
4851 track->UpdatePoints();
4852 if (( track->GetPoints()[2]- track->GetPoints()[0])>5 && track->GetPoints()[3]>0.8){
4854 nclusters[i]=track->GetNumberOfClusters();
4855 alpha[i] = track->GetAlpha();
4856 new (&helixes[i]) AliHelix(*track);
4858 helixes[i].Evaluate(0,xyz);
4859 sign[i] = (track->GetC()>0) ? -1:1;
4862 if (track->GetProlongation(x,y,z)){
4864 fim[i] = alpha[i]+TMath::ATan2(y,x);
4867 zm[i] = track->GetZ();
4871 circular[i]= kFALSE;
4872 if (track->GetProlongation(0,y,z)) z0[i] = z;
4873 dca[i] = track->GetD(0,0);
4879 Int_t ncandidates =0;
4882 Double_t phase[2][2]={{0,0},{0,0}},radius[2]={0,0};
4885 // Find circling track
4887 for (Int_t i0=0;i0<nentries;i0++){
4888 AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
4889 if (!track0) continue;
4890 if (track0->GetNumberOfClusters()<40) continue;
4891 if (TMath::Abs(1./track0->GetC())>200) continue;
4892 for (Int_t i1=i0+1;i1<nentries;i1++){
4893 AliTPCseed * track1 = (AliTPCseed*)array->At(i1);
4894 if (!track1) continue;
4895 if (track1->GetNumberOfClusters()<40) continue;
4896 if ( TMath::Abs(track1->GetTgl()+track0->GetTgl())>0.1) continue;
4897 if (track0->GetBConstrain()&&track1->GetBConstrain()) continue;
4898 if (TMath::Abs(1./track1->GetC())>200) continue;
4899 if (track1->GetSigned1Pt()*track0->GetSigned1Pt()>0) continue;
4900 if (track1->GetTgl()*track0->GetTgl()>0) continue;
4901 if (TMath::Max(TMath::Abs(1./track0->GetC()),TMath::Abs(1./track1->GetC()))>190) continue;
4902 if (track0->GetBConstrain()&&track1->OneOverPt()<track0->OneOverPt()) continue; //returning - lower momenta
4903 if (track1->GetBConstrain()&&track0->OneOverPt()<track1->OneOverPt()) continue; //returning - lower momenta
4905 Float_t mindcar = TMath::Min(TMath::Abs(dca[i0]),TMath::Abs(dca[i1]));
4906 if (mindcar<5) continue;
4907 Float_t mindcaz = TMath::Min(TMath::Abs(z0[i0]-GetZ()),TMath::Abs(z0[i1]-GetZ()));
4908 if (mindcaz<5) continue;
4909 if (mindcar+mindcaz<20) continue;
4912 Float_t xc0 = helixes[i0].GetHelix(6);
4913 Float_t yc0 = helixes[i0].GetHelix(7);
4914 Float_t r0 = helixes[i0].GetHelix(8);
4915 Float_t xc1 = helixes[i1].GetHelix(6);
4916 Float_t yc1 = helixes[i1].GetHelix(7);
4917 Float_t r1 = helixes[i1].GetHelix(8);
4919 Float_t rmean = (r0+r1)*0.5;
4920 Float_t delta =TMath::Sqrt((xc1-xc0)*(xc1-xc0)+(yc1-yc0)*(yc1-yc0));
4921 //if (delta>30) continue;
4922 if (delta>rmean*0.25) continue;
4923 if (TMath::Abs(r0-r1)/rmean>0.3) continue;
4925 Int_t npoints = helixes[i0].GetRPHIintersections(helixes[i1], phase, radius,10);
4926 if (npoints==0) continue;
4927 helixes[i0].GetClosestPhases(helixes[i1], phase);
4931 Double_t hangles[3];
4932 helixes[i0].Evaluate(phase[0][0],xyz0);
4933 helixes[i1].Evaluate(phase[0][1],xyz1);
4935 helixes[i0].GetAngle(phase[0][0],helixes[i1],phase[0][1],hangles);
4936 Double_t deltah[2],deltabest;
4937 if (hangles[2]<2.8) continue;
4940 helixes[i0].ParabolicDCA(helixes[i1],phase[0][0],phase[0][1],radius[0],deltah[0],2);
4942 helixes[i0].ParabolicDCA(helixes[i1],phase[1][0],phase[1][1],radius[1],deltah[1],2);
4943 if (deltah[1]<deltah[0]) ibest=1;
4945 deltabest = TMath::Sqrt(deltah[ibest]);
4946 helixes[i0].Evaluate(phase[ibest][0],xyz0);
4947 helixes[i1].Evaluate(phase[ibest][1],xyz1);
4948 helixes[i0].GetAngle(phase[ibest][0],helixes[i1],phase[ibest][1],hangles);
4949 Double_t radiusbest = TMath::Sqrt(radius[ibest]);
4951 if (deltabest>6) continue;
4952 if (mindcar+mindcaz<40 && (hangles[2]<3.12||deltabest>3)) continue;
4953 Bool_t lsign =kFALSE;
4954 if (hangles[2]>3.06) lsign =kTRUE;
4957 circular[i0] = kTRUE;
4958 circular[i1] = kTRUE;
4959 if (track0->OneOverPt()<track1->OneOverPt()){
4960 track0->SetCircular(track0->GetCircular()+1);
4961 track1->SetCircular(track1->GetCircular()+2);
4964 track1->SetCircular(track1->GetCircular()+1);
4965 track0->SetCircular(track0->GetCircular()+2);
4968 if (lsign&&AliTPCReconstructor::StreamLevel()>1){
4970 Int_t lab0=track0->GetLabel();
4971 Int_t lab1=track1->GetLabel();
4972 TTreeSRedirector &cstream = *fDebugStreamer;
4973 cstream<<"Curling"<<
4980 "mindcar="<<mindcar<<
4981 "mindcaz="<<mindcaz<<
4984 "npoints="<<npoints<<
4985 "hangles0="<<hangles[0]<<
4986 "hangles2="<<hangles[2]<<
4991 "radius="<<radiusbest<<
4992 "deltabest="<<deltabest<<
4993 "phase0="<<phase[ibest][0]<<
4994 "phase1="<<phase[ibest][1]<<
5004 for (Int_t i =0;i<nentries;i++){
5005 if (sign[i]==0) continue;
5006 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5013 Double_t cradius0 = 40*40;
5014 Double_t cradius1 = 270*270;
5017 Double_t cdist3=0.55;
5018 for (Int_t j =i+1;j<nentries;j++){
5020 if (sign[j]*sign[i]<1) continue;
5021 if ( (nclusters[i]+nclusters[j])>200) continue;
5022 if ( (nclusters[i]+nclusters[j])<80) continue;
5023 if ( TMath::Abs(zm[i]-zm[j])>60.) continue;
5024 if ( TMath::Abs(fim[i]-fim[j])>0.6 && TMath::Abs(fim[i]-fim[j])<5.7 ) continue;
5025 //AliTPCseed * track1 = (AliTPCseed*)array->At(j); Double_t phase[2][2],radius[2];
5026 Int_t npoints = helixes[i].GetRPHIintersections(helixes[j], phase, radius,20);
5027 if (npoints<1) continue;
5030 if (radius[0]<cradius0||radius[0]>cradius1) continue;
5033 if ( (radius[0]<cradius0||radius[0]>cradius1) && (radius[1]<cradius0||radius[1]>cradius1) ) continue;
5036 Double_t delta1=10000,delta2=10000;
5037 // cuts on the intersection radius
5038 helixes[i].LinearDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta1);
5039 if (radius[0]<20&&delta1<1) continue; //intersection at vertex
5040 if (radius[0]<10&&delta1<3) continue; //intersection at vertex
5042 helixes[i].LinearDCA(helixes[j],phase[1][0],phase[1][1],radius[1],delta2);
5043 if (radius[1]<20&&delta2<1) continue; //intersection at vertex
5044 if (radius[1]<10&&delta2<3) continue; //intersection at vertex
5047 Double_t distance1 = TMath::Min(delta1,delta2);
5048 if (distance1>cdist1) continue; // cut on DCA linear approximation
5050 npoints = helixes[i].GetRPHIintersections(helixes[j], phase, radius,20);
5051 helixes[i].ParabolicDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta1);
5052 if (radius[0]<20&&delta1<1) continue; //intersection at vertex
5053 if (radius[0]<10&&delta1<3) continue; //intersection at vertex
5056 helixes[i].ParabolicDCA(helixes[j],phase[1][0],phase[1][1],radius[1],delta2);
5057 if (radius[1]<20&&delta2<1) continue; //intersection at vertex
5058 if (radius[1]<10&&delta2<3) continue; //intersection at vertex
5060 distance1 = TMath::Min(delta1,delta2);
5063 rkink = TMath::Sqrt(radius[0]);
5066 rkink = TMath::Sqrt(radius[1]);
5068 if (distance1>cdist2) continue;
5071 AliTPCseed * track1 = (AliTPCseed*)array->At(j);
5074 Int_t row0 = GetRowNumber(rkink);
5075 if (row0<10) continue;
5076 if (row0>150) continue;
5079 Float_t dens00=-1,dens01=-1;
5080 Float_t dens10=-1,dens11=-1;
5082 Int_t found,foundable,ishared;
5083 track0->GetClusterStatistic(0,row0-5, found, foundable,ishared,kFALSE);
5084 if (foundable>5) dens00 = Float_t(found)/Float_t(foundable);
5085 track0->GetClusterStatistic(row0+5,155, found, foundable,ishared,kFALSE);
5086 if (foundable>5) dens01 = Float_t(found)/Float_t(foundable);
5088 track1->GetClusterStatistic(0,row0-5, found, foundable,ishared,kFALSE);
5089 if (foundable>10) dens10 = Float_t(found)/Float_t(foundable);
5090 track1->GetClusterStatistic(row0+5,155, found, foundable,ishared,kFALSE);
5091 if (foundable>10) dens11 = Float_t(found)/Float_t(foundable);
5093 if (dens00<dens10 && dens01<dens11) continue;
5094 if (dens00>dens10 && dens01>dens11) continue;
5095 if (TMath::Max(dens00,dens10)<0.1) continue;
5096 if (TMath::Max(dens01,dens11)<0.3) continue;
5098 if (TMath::Min(dens00,dens10)>0.6) continue;
5099 if (TMath::Min(dens01,dens11)>0.6) continue;
5102 AliTPCseed * ktrack0, *ktrack1;
5111 if (TMath::Abs(ktrack0->GetC())>5) continue; // cut on the curvature for mother particle
5112 AliExternalTrackParam paramm(*ktrack0);
5113 AliExternalTrackParam paramd(*ktrack1);
5114 if (row0>60&&ktrack1->GetReference().GetX()>90.)new (¶md) AliExternalTrackParam(ktrack1->GetReference());
5117 kink->SetMother(paramm);
5118 kink->SetDaughter(paramd);
5121 Float_t x[3] = { kink->GetPosition()[0],kink->GetPosition()[1],kink->GetPosition()[2]};
5123 fkParam->Transform0to1(x,index);
5124 fkParam->Transform1to2(x,index);
5125 row0 = GetRowNumber(x[0]);
5127 if (kink->GetR()<100) continue;
5128 if (kink->GetR()>240) continue;
5129 if (kink->GetPosition()[2]/kink->GetR()>AliTPCReconstructor::GetCtgRange()) continue; //out of fiducial volume
5130 if (kink->GetDistance()>cdist3) continue;
5131 Float_t dird = kink->GetDaughterP()[0]*kink->GetPosition()[0]+kink->GetDaughterP()[1]*kink->GetPosition()[1]; // rough direction estimate
5132 if (dird<0) continue;
5134 Float_t dirm = kink->GetMotherP()[0]*kink->GetPosition()[0]+kink->GetMotherP()[1]*kink->GetPosition()[1]; // rough direction estimate
5135 if (dirm<0) continue;
5136 Float_t mpt = TMath::Sqrt(kink->GetMotherP()[0]*kink->GetMotherP()[0]+kink->GetMotherP()[1]*kink->GetMotherP()[1]);
5137 if (mpt<0.2) continue;
5140 //for high momenta momentum not defined well in first iteration
5141 Double_t qt = TMath::Sin(kink->GetAngle(2))*ktrack1->GetP();
5142 if (qt>0.35) continue;
5145 kink->SetLabel(CookLabel(ktrack0,0.4,0,row0),0);
5146 kink->SetLabel(CookLabel(ktrack1,0.4,row0,160),1);
5148 kink->SetTPCDensity(dens00,0,0);
5149 kink->SetTPCDensity(dens01,0,1);
5150 kink->SetTPCDensity(dens10,1,0);
5151 kink->SetTPCDensity(dens11,1,1);
5152 kink->SetIndex(i,0);
5153 kink->SetIndex(j,1);
5156 kink->SetTPCDensity(dens10,0,0);
5157 kink->SetTPCDensity(dens11,0,1);
5158 kink->SetTPCDensity(dens00,1,0);
5159 kink->SetTPCDensity(dens01,1,1);
5160 kink->SetIndex(j,0);
5161 kink->SetIndex(i,1);
5164 if (mpt<1||kink->GetAngle(2)>0.1){
5165 // angle and densities not defined yet
5166 if (kink->GetTPCDensityFactor()<0.8) continue;
5167 if ((2-kink->GetTPCDensityFactor())*kink->GetDistance() >0.25) continue;
5168 if (kink->GetAngle(2)*ktrack0->GetP()<0.003) continue; //too small angle
5169 if (kink->GetAngle(2)>0.2&&kink->GetTPCDensityFactor()<1.15) continue;
5170 if (kink->GetAngle(2)>0.2&&kink->GetTPCDensity(0,1)>0.05) continue;
5172 Float_t criticalangle = track0->GetSigmaSnp2()+track0->GetSigmaTgl2();
5173 criticalangle+= track1->GetSigmaSnp2()+track1->GetSigmaTgl2();
5174 criticalangle= 3*TMath::Sqrt(criticalangle);
5175 if (criticalangle>0.02) criticalangle=0.02;
5176 if (kink->GetAngle(2)<criticalangle) continue;
5179 Int_t drow = Int_t(2.+0.5/(0.05+kink->GetAngle(2))); // overlap region defined
5180 Float_t shapesum =0;
5182 for ( Int_t row = row0-drow; row<row0+drow;row++){
5183 if (row<0) continue;
5184 if (row>155) continue;
5185 if (ktrack0->GetClusterPointer(row)){
5186 AliTPCTrackerPoint *point =ktrack0->GetTrackPoint(row);
5187 shapesum+=point->GetSigmaY()+point->GetSigmaZ();
5190 if (ktrack1->GetClusterPointer(row)){
5191 AliTPCTrackerPoint *point =ktrack1->GetTrackPoint(row);
5192 shapesum+=point->GetSigmaY()+point->GetSigmaZ();
5197 kink->SetShapeFactor(-1.);
5200 kink->SetShapeFactor(shapesum/sum);
5202 // esd->AddKink(kink);
5204 // kink->SetMother(paramm);
5205 //kink->SetDaughter(paramd);
5207 Double_t chi2P2 = paramm.GetParameter()[2]-paramd.GetParameter()[2];
5209 chi2P2/=paramm.GetCovariance()[5]+paramd.GetCovariance()[5];
5210 Double_t chi2P3 = paramm.GetParameter()[3]-paramd.GetParameter()[3];
5212 chi2P3/=paramm.GetCovariance()[9]+paramd.GetCovariance()[9];
5214 if (AliTPCReconstructor::StreamLevel()>1) {
5215 (*fDebugStreamer)<<"kinkLpt"<<
5223 if ( chi2P2+chi2P3<AliTPCReconstructor::GetRecoParam()->GetKinkAngleCutChi2(0)){
5227 kinks->AddLast(kink);
5233 // sort the kinks according quality - and refit them towards vertex
5235 Int_t nkinks = kinks->GetEntriesFast();
5236 Float_t *quality = new Float_t[nkinks];
5237 Int_t *indexes = new Int_t[nkinks];
5238 AliTPCseed *mothers = new AliTPCseed[nkinks];
5239 AliTPCseed *daughters = new AliTPCseed[nkinks];
5242 for (Int_t i=0;i<nkinks;i++){
5244 AliKink *kinkl = (AliKink*)kinks->At(i);
5246 // refit kinks towards vertex
5248 Int_t index0 = kinkl->GetIndex(0);
5249 Int_t index1 = kinkl->GetIndex(1);
5250 AliTPCseed * ktrack0 = (AliTPCseed*)array->At(index0);
5251 AliTPCseed * ktrack1 = (AliTPCseed*)array->At(index1);
5253 Int_t sumn=ktrack0->GetNumberOfClusters()+ktrack1->GetNumberOfClusters();
5255 // Refit Kink under if too small angle
5257 if (kinkl->GetAngle(2)<0.05){
5258 kinkl->SetTPCRow0(GetRowNumber(kinkl->GetR()));
5259 Int_t row0 = kinkl->GetTPCRow0();
5260 Int_t drow = Int_t(2.+0.5/(0.05+kinkl->GetAngle(2)));
5263 Int_t last = row0-drow;
5264 if (last<40) last=40;
5265 if (last<ktrack0->GetFirstPoint()+25) last = ktrack0->GetFirstPoint()+25;
5266 AliTPCseed* seed0 = ReSeed(ktrack0,last,kFALSE);
5269 Int_t first = row0+drow;
5270 if (first>130) first=130;
5271 if (first>ktrack1->GetLastPoint()-25) first = TMath::Max(ktrack1->GetLastPoint()-25,30);
5272 AliTPCseed* seed1 = ReSeed(ktrack1,first,kTRUE);
5274 if (seed0 && seed1){
5275 kinkl->SetStatus(1,8);
5276 if (RefitKink(*seed0,*seed1,*kinkl)) kinkl->SetStatus(1,9);
5277 row0 = GetRowNumber(kinkl->GetR());
5278 sumn = seed0->GetNumberOfClusters()+seed1->GetNumberOfClusters();
5279 mothers[i] = *seed0;
5280 daughters[i] = *seed1;
5283 delete kinks->RemoveAt(i);
5284 if (seed0) MarkSeedFree( seed0 );
5285 if (seed1) MarkSeedFree( seed1 );
5288 if (kinkl->GetDistance()>0.5 || kinkl->GetR()<110 || kinkl->GetR()>240) {
5289 delete kinks->RemoveAt(i);
5290 if (seed0) MarkSeedFree( seed0 );
5291 if (seed1) MarkSeedFree( seed1 );
5295 MarkSeedFree( seed0 );
5296 MarkSeedFree( seed1 );
5299 if (kinkl) quality[i] = 160*((0.1+kinkl->GetDistance())*(2.-kinkl->GetTPCDensityFactor()))/(sumn+40.); //the longest -clossest will win
5301 TMath::Sort(nkinks,quality,indexes,kFALSE);
5303 //remove double find kinks
5305 for (Int_t ikink0=1;ikink0<nkinks;ikink0++){
5306 AliKink * kink0 = (AliKink*) kinks->At(indexes[ikink0]);
5307 if (!kink0) continue;
5309 for (Int_t ikink1=0;ikink1<ikink0;ikink1++){
5310 kink0 = (AliKink*) kinks->At(indexes[ikink0]);
5311 if (!kink0) continue;
5312 AliKink * kink1 = (AliKink*) kinks->At(indexes[ikink1]);
5313 if (!kink1) continue;
5314 // if not close kink continue
5315 if (TMath::Abs(kink1->GetPosition()[2]-kink0->GetPosition()[2])>10) continue;
5316 if (TMath::Abs(kink1->GetPosition()[1]-kink0->GetPosition()[1])>10) continue;
5317 if (TMath::Abs(kink1->GetPosition()[0]-kink0->GetPosition()[0])>10) continue;
5319 AliTPCseed &mother0 = mothers[indexes[ikink0]];
5320 AliTPCseed &daughter0 = daughters[indexes[ikink0]];
5321 AliTPCseed &mother1 = mothers[indexes[ikink1]];
5322 AliTPCseed &daughter1 = daughters[indexes[ikink1]];
5323 Int_t row0 = (kink0->GetTPCRow0()+kink1->GetTPCRow0())/2;
5332 for (Int_t i=0;i<row0;i++){
5333 if (mother0.GetClusterIndex(i)>0 && mother1.GetClusterIndex(i)>0){
5336 if (mother0.GetClusterIndex(i)==mother1.GetClusterIndex(i)){
5343 for (Int_t i=row0;i<158;i++){
5344 //if (daughter0.GetClusterIndex(i)>0 && daughter0.GetClusterIndex(i)>0){ // RS: Bug?
5345 if (daughter0.GetClusterIndex(i)>0 && daughter1.GetClusterIndex(i)>0){
5348 if (mother0.GetClusterIndex(i)==mother1.GetClusterIndex(i)){
5354 Float_t ratio = Float_t(same+1)/Float_t(both+1);
5355 Float_t ratiom = Float_t(samem+1)/Float_t(bothm+1);
5356 Float_t ratiod = Float_t(samed+1)/Float_t(bothd+1);
5357 if (ratio>0.3 && ratiom>0.5 &&ratiod>0.5) {
5358 Int_t sum0 = mother0.GetNumberOfClusters()+daughter0.GetNumberOfClusters();
5359 Int_t sum1 = mother1.GetNumberOfClusters()+daughter1.GetNumberOfClusters();
5361 shared[kink0->GetIndex(0)]= kTRUE;
5362 shared[kink0->GetIndex(1)]= kTRUE;
5363 delete kinks->RemoveAt(indexes[ikink0]);
5367 shared[kink1->GetIndex(0)]= kTRUE;
5368 shared[kink1->GetIndex(1)]= kTRUE;
5369 delete kinks->RemoveAt(indexes[ikink1]);
5376 for (Int_t i=0;i<nkinks;i++){
5377 AliKink * kinkl = (AliKink*) kinks->At(indexes[i]);
5378 if (!kinkl) continue;
5379 kinkl->SetTPCRow0(GetRowNumber(kinkl->GetR()));
5380 Int_t index0 = kinkl->GetIndex(0);
5381 Int_t index1 = kinkl->GetIndex(1);
5382 if (circular[index0]||(circular[index1]&&kinkl->GetDistance()>0.2)) continue;
5383 kinkl->SetMultiple(usage[index0],0);
5384 kinkl->SetMultiple(usage[index1],1);
5385 if (kinkl->GetMultiple()[0]+kinkl->GetMultiple()[1]>2) continue;
5386 if (kinkl->GetMultiple()[0]+kinkl->GetMultiple()[1]>0 && quality[indexes[i]]>0.2) continue;
5387 if (kinkl->GetMultiple()[0]+kinkl->GetMultiple()[1]>0 && kinkl->GetDistance()>0.2) continue;
5388 if (circular[index0]||(circular[index1]&&kinkl->GetDistance()>0.1)) continue;
5390 AliTPCseed * ktrack0 = (AliTPCseed*)array->At(index0);
5391 AliTPCseed * ktrack1 = (AliTPCseed*)array->At(index1);
5392 if (!ktrack0 || !ktrack1) continue;
5393 Int_t index = esd->AddKink(kinkl);
5396 if ( ktrack0->GetKinkIndex(0)==0 && ktrack1->GetKinkIndex(0)==0) { //best kink
5397 if (mothers[indexes[i]].GetNumberOfClusters()>20 && daughters[indexes[i]].GetNumberOfClusters()>20 && (mothers[indexes[i]].GetNumberOfClusters()+daughters[indexes[i]].GetNumberOfClusters())>100){
5398 *ktrack0 = mothers[indexes[i]];
5399 *ktrack1 = daughters[indexes[i]];
5403 ktrack0->SetKinkIndex(usage[index0],-(index+1));
5404 ktrack1->SetKinkIndex(usage[index1], (index+1));
5409 // Remove tracks corresponding to shared kink's
5411 for (Int_t i=0;i<nentries;i++){
5412 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5413 if (!track0) continue;
5414 if (track0->GetKinkIndex(0)!=0) continue;
5415 if (shared[i]) MarkSeedFree( array->RemoveAt(i) );
5420 RemoveUsed2(array,0.5,0.4,30);
5422 for (Int_t i=0;i<nentries;i++){
5423 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5424 if (!track0) continue;
5425 track0->CookdEdx(0.02,0.6);
5429 for (Int_t i=0;i<nentries;i++){
5430 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5431 if (!track0) continue;
5432 if (track0->Pt()<1.4) continue;
5433 //remove double high momenta tracks - overlapped with kink candidates
5436 for (Int_t icl=track0->GetFirstPoint();icl<track0->GetLastPoint(); icl++){
5437 if (track0->GetClusterPointer(icl)!=0){
5439 if (track0->GetClusterPointer(icl)->IsUsed(10)) ishared++;
5442 if (Float_t(ishared+1)/Float_t(all+1)>0.5) {
5443 MarkSeedFree( array->RemoveAt(i) );
5447 if (track0->GetKinkIndex(0)!=0) continue;
5448 if (track0->GetNumberOfClusters()<80) continue;
5450 AliTPCseed *pmother = new AliTPCseed();
5451 AliTPCseed *pdaughter = new AliTPCseed();
5452 AliKink *pkink = new AliKink;
5454 AliTPCseed & mother = *pmother;
5455 AliTPCseed & daughter = *pdaughter;
5456 AliKink & kinkl = *pkink;
5457 if (CheckKinkPoint(track0,mother,daughter, kinkl)){
5458 if (mother.GetNumberOfClusters()<30||daughter.GetNumberOfClusters()<20) {
5462 continue; //too short tracks
5464 if (mother.Pt()<1.4) {
5470 Int_t row0= kinkl.GetTPCRow0();
5471 if (kinkl.GetDistance()>0.5 || kinkl.GetR()<110. || kinkl.GetR()>240.) {
5478 Int_t index = esd->AddKink(&kinkl);
5479 mother.SetKinkIndex(0,-(index+1));
5480 daughter.SetKinkIndex(0,index+1);
5481 if (mother.GetNumberOfClusters()>50) {
5482 MarkSeedFree( array->RemoveAt(i) );
5483 AliTPCseed* mtc = new( NextFreeSeed() ) AliTPCseed(mother);
5484 mtc->SetPoolID(fLastSeedID);
5485 array->AddAt(mtc,i);
5488 AliTPCseed* mtc = new( NextFreeSeed() ) AliTPCseed(mother);
5489 mtc->SetPoolID(fLastSeedID);
5490 array->AddLast(mtc);
5492 AliTPCseed* dtc = new( NextFreeSeed() ) AliTPCseed(daughter);
5493 dtc->SetPoolID(fLastSeedID);
5494 array->AddLast(dtc);
5495 for (Int_t icl=0;icl<row0;icl++) {
5496 if (mother.GetClusterPointer(icl)) mother.GetClusterPointer(icl)->Use(20);
5499 for (Int_t icl=row0;icl<158;icl++) {
5500 if (daughter.GetClusterPointer(icl)) daughter.GetClusterPointer(icl)->Use(20);
5509 delete [] daughters;
5531 AliInfo(Form("Ncandidates=\t%d\t%d\t%d\t%d\n",esd->GetNumberOfKinks(),ncandidates,ntracks,nall));
5537 void AliTPCtrackerMI::FindKinks(TObjArray * array, AliESDEvent *esd)
5544 TObjArray *kinks= new TObjArray(10000);
5545 // TObjArray *v0s= new TObjArray(10000);
5546 Int_t nentries = array->GetEntriesFast();
5547 AliHelix *helixes = new AliHelix[nentries];
5548 Int_t *sign = new Int_t[nentries];
5549 Int_t *nclusters = new Int_t[nentries];
5550 Float_t *alpha = new Float_t[nentries];
5551 AliKink *kink = new AliKink();
5552 Int_t * usage = new Int_t[nentries];
5553 Float_t *zm = new Float_t[nentries];
5554 Float_t *z0 = new Float_t[nentries];
5555 Float_t *fim = new Float_t[nentries];
5556 Float_t *shared = new Float_t[nentries];
5557 Bool_t *circular = new Bool_t[nentries];
5558 Float_t *dca = new Float_t[nentries];
5559 //const AliESDVertex * primvertex = esd->GetVertex();
5561 // nentries = array->GetEntriesFast();
5566 for (Int_t i=0;i<nentries;i++){
5569 AliTPCseed* track = (AliTPCseed*)array->At(i);
5570 if (!track) continue;
5571 track->SetCircular(0);
5573 track->UpdatePoints();
5574 if (( track->GetPoints()[2]- track->GetPoints()[0])>5 && track->GetPoints()[3]>0.8){
5576 nclusters[i]=track->GetNumberOfClusters();
5577 alpha[i] = track->GetAlpha();
5578 new (&helixes[i]) AliHelix(*track);
5580 helixes[i].Evaluate(0,xyz);
5581 sign[i] = (track->GetC()>0) ? -1:1;
5584 if (track->GetProlongation(x,y,z)){
5586 fim[i] = alpha[i]+TMath::ATan2(y,x);
5589 zm[i] = track->GetZ();
5593 circular[i]= kFALSE;
5594 if (track->GetProlongation(0,y,z)) z0[i] = z;
5595 dca[i] = track->GetD(0,0);
5601 Int_t ncandidates =0;
5604 Double_t phase[2][2]={{0,0},{0,0}},radius[2]={0,0};
5607 // Find circling track
5609 for (Int_t i0=0;i0<nentries;i0++){
5610 AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
5611 if (!track0) continue;
5612 if (track0->GetNumberOfClusters()<40) continue;
5613 if (TMath::Abs(1./track0->GetC())>200) continue;
5614 for (Int_t i1=i0+1;i1<nentries;i1++){
5615 AliTPCseed * track1 = (AliTPCseed*)array->At(i1);
5616 if (!track1) continue;
5617 if (track1->GetNumberOfClusters()<40) continue;
5618 if ( TMath::Abs(track1->GetTgl()+track0->GetTgl())>0.1) continue;
5619 if (track0->GetBConstrain()&&track1->GetBConstrain()) continue;
5620 if (TMath::Abs(1./track1->GetC())>200) continue;
5621 if (track1->GetSigned1Pt()*track0->GetSigned1Pt()>0) continue;
5622 if (track1->GetTgl()*track0->GetTgl()>0) continue;
5623 if (TMath::Max(TMath::Abs(1./track0->GetC()),TMath::Abs(1./track1->GetC()))>190) continue;
5624 if (track0->GetBConstrain()&&track1->OneOverPt()<track0->OneOverPt()) continue; //returning - lower momenta
5625 if (track1->GetBConstrain()&&track0->OneOverPt()<track1->OneOverPt()) continue; //returning - lower momenta
5627 Float_t mindcar = TMath::Min(TMath::Abs(dca[i0]),TMath::Abs(dca[i1]));
5628 if (mindcar<5) continue;
5629 Float_t mindcaz = TMath::Min(TMath::Abs(z0[i0]-GetZ()),TMath::Abs(z0[i1]-GetZ()));
5630 if (mindcaz<5) continue;
5631 if (mindcar+mindcaz<20) continue;
5634 Float_t xc0 = helixes[i0].GetHelix(6);
5635 Float_t yc0 = helixes[i0].GetHelix(7);
5636 Float_t r0 = helixes[i0].GetHelix(8);
5637 Float_t xc1 = helixes[i1].GetHelix(6);
5638 Float_t yc1 = helixes[i1].GetHelix(7);
5639 Float_t r1 = helixes[i1].GetHelix(8);
5641 Float_t rmean = (r0+r1)*0.5;
5642 Float_t delta =TMath::Sqrt((xc1-xc0)*(xc1-xc0)+(yc1-yc0)*(yc1-yc0));
5643 //if (delta>30) continue;
5644 if (delta>rmean*0.25) continue;
5645 if (TMath::Abs(r0-r1)/rmean>0.3) continue;
5647 Int_t npoints = helixes[i0].GetRPHIintersections(helixes[i1], phase, radius,10);
5648 if (npoints==0) continue;
5649 helixes[i0].GetClosestPhases(helixes[i1], phase);
5653 Double_t hangles[3];
5654 helixes[i0].Evaluate(phase[0][0],xyz0);
5655 helixes[i1].Evaluate(phase[0][1],xyz1);
5657 helixes[i0].GetAngle(phase[0][0],helixes[i1],phase[0][1],hangles);
5658 Double_t deltah[2],deltabest;
5659 if (hangles[2]<2.8) continue;
5662 helixes[i0].ParabolicDCA(helixes[i1],phase[0][0],phase[0][1],radius[0],deltah[0],2);
5664 helixes[i0].ParabolicDCA(helixes[i1],phase[1][0],phase[1][1],radius[1],deltah[1],2);
5665 if (deltah[1]<deltah[0]) ibest=1;
5667 deltabest = TMath::Sqrt(deltah[ibest]);
5668 helixes[i0].Evaluate(phase[ibest][0],xyz0);
5669 helixes[i1].Evaluate(phase[ibest][1],xyz1);
5670 helixes[i0].GetAngle(phase[ibest][0],helixes[i1],phase[ibest][1],hangles);
5671 Double_t radiusbest = TMath::Sqrt(radius[ibest]);
5673 if (deltabest>6) continue;
5674 if (mindcar+mindcaz<40 && (hangles[2]<3.12||deltabest>3)) continue;
5675 Bool_t lsign =kFALSE;
5676 if (hangles[2]>3.06) lsign =kTRUE;
5679 circular[i0] = kTRUE;
5680 circular[i1] = kTRUE;
5681 if (track0->OneOverPt()<track1->OneOverPt()){
5682 track0->SetCircular(track0->GetCircular()+1);
5683 track1->SetCircular(track1->GetCircular()+2);
5686 track1->SetCircular(track1->GetCircular()+1);
5687 track0->SetCircular(track0->GetCircular()+2);
5690 if (lsign&&AliTPCReconstructor::StreamLevel()>1){
5692 Int_t lab0=track0->GetLabel();
5693 Int_t lab1=track1->GetLabel();
5694 TTreeSRedirector &cstream = *fDebugStreamer;
5695 cstream<<"Curling"<<
5702 "mindcar="<<mindcar<<
5703 "mindcaz="<<mindcaz<<
5706 "npoints="<<npoints<<
5707 "hangles0="<<hangles[0]<<
5708 "hangles2="<<hangles[2]<<
5713 "radius="<<radiusbest<<
5714 "deltabest="<<deltabest<<
5715 "phase0="<<phase[ibest][0]<<
5716 "phase1="<<phase[ibest][1]<<
5726 for (Int_t i =0;i<nentries;i++){
5727 if (sign[i]==0) continue;
5728 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5735 Double_t cradius0 = 40*40;
5736 Double_t cradius1 = 270*270;
5739 Double_t cdist3=0.55;
5740 for (Int_t j =i+1;j<nentries;j++){
5742 if (sign[j]*sign[i]<1) continue;
5743 if ( (nclusters[i]+nclusters[j])>200) continue;
5744 if ( (nclusters[i]+nclusters[j])<80) continue;
5745 if ( TMath::Abs(zm[i]-zm[j])>60.) continue;
5746 if ( TMath::Abs(fim[i]-fim[j])>0.6 && TMath::Abs(fim[i]-fim[j])<5.7 ) continue;
5747 //AliTPCseed * track1 = (AliTPCseed*)array->At(j); Double_t phase[2][2],radius[2];
5748 Int_t npoints = helixes[i].GetRPHIintersections(helixes[j], phase, radius,20);
5749 if (npoints<1) continue;
5752 if (radius[0]<cradius0||radius[0]>cradius1) continue;
5755 if ( (radius[0]<cradius0||radius[0]>cradius1) && (radius[1]<cradius0||radius[1]>cradius1) ) continue;
5758 Double_t delta1=10000,delta2=10000;
5759 // cuts on the intersection radius
5760 helixes[i].LinearDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta1);
5761 if (radius[0]<20&&delta1<1) continue; //intersection at vertex
5762 if (radius[0]<10&&delta1<3) continue; //intersection at vertex
5764 helixes[i].LinearDCA(helixes[j],phase[1][0],phase[1][1],radius[1],delta2);
5765 if (radius[1]<20&&delta2<1) continue; //intersection at vertex
5766 if (radius[1]<10&&delta2<3) continue; //intersection at vertex
5769 Double_t distance1 = TMath::Min(delta1,delta2);
5770 if (distance1>cdist1) continue; // cut on DCA linear approximation
5772 npoints = helixes[i].GetRPHIintersections(helixes[j], phase, radius,20);
5773 helixes[i].ParabolicDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta1);
5774 if (radius[0]<20&&delta1<1) continue; //intersection at vertex
5775 if (radius[0]<10&&delta1<3) continue; //intersection at vertex
5778 helixes[i].ParabolicDCA(helixes[j],phase[1][0],phase[1][1],radius[1],delta2);
5779 if (radius[1]<20&&delta2<1) continue; //intersection at vertex
5780 if (radius[1]<10&&delta2<3) continue; //intersection at vertex
5782 distance1 = TMath::Min(delta1,delta2);
5785 rkink = TMath::Sqrt(radius[0]);
5788 rkink = TMath::Sqrt(radius[1]);
5790 if (distance1>cdist2) continue;
5793 AliTPCseed * track1 = (AliTPCseed*)array->At(j);
5796 Int_t row0 = GetRowNumber(rkink);
5797 if (row0<10) continue;
5798 if (row0>150) continue;
5801 Float_t dens00=-1,dens01=-1;
5802 Float_t dens10=-1,dens11=-1;
5804 Int_t found,foundable,ishared;
5805 track0->GetClusterStatistic(0,row0-5, found, foundable,ishared,kFALSE);
5806 if (foundable>5) dens00 = Float_t(found)/Float_t(foundable);
5807 track0->GetClusterStatistic(row0+5,155, found, foundable,ishared,kFALSE);
5808 if (foundable>5) dens01 = Float_t(found)/Float_t(foundable);
5810 track1->GetClusterStatistic(0,row0-5, found, foundable,ishared,kFALSE);
5811 if (foundable>10) dens10 = Float_t(found)/Float_t(foundable);
5812 track1->GetClusterStatistic(row0+5,155, found, foundable,ishared,kFALSE);
5813 if (foundable>10) dens11 = Float_t(found)/Float_t(foundable);
5815 if (dens00<dens10 && dens01<dens11) continue;
5816 if (dens00>dens10 && dens01>dens11) continue;
5817 if (TMath::Max(dens00,dens10)<0.1) continue;
5818 if (TMath::Max(dens01,dens11)<0.3) continue;
5820 if (TMath::Min(dens00,dens10)>0.6) continue;
5821 if (TMath::Min(dens01,dens11)>0.6) continue;
5824 AliTPCseed * ktrack0, *ktrack1;
5833 if (TMath::Abs(ktrack0->GetC())>5) continue; // cut on the curvature for mother particle
5834 AliExternalTrackParam paramm(*ktrack0);
5835 AliExternalTrackParam paramd(*ktrack1);
5836 if (row0>60&&ktrack1->GetReference().GetX()>90.)new (¶md) AliExternalTrackParam(ktrack1->GetReference());
5839 kink->SetMother(paramm);
5840 kink->SetDaughter(paramd);
5843 Float_t x[3] = { kink->GetPosition()[0],kink->GetPosition()[1],kink->GetPosition()[2]};
5845 fkParam->Transform0to1(x,index);
5846 fkParam->Transform1to2(x,index);
5847 row0 = GetRowNumber(x[0]);
5849 if (kink->GetR()<100) continue;
5850 if (kink->GetR()>240) continue;
5851 if (kink->GetPosition()[2]/kink->GetR()>AliTPCReconstructor::GetCtgRange()) continue; //out of fiducial volume
5852 if (kink->GetDistance()>cdist3) continue;
5853 Float_t dird = kink->GetDaughterP()[0]*kink->GetPosition()[0]+kink->GetDaughterP()[1]*kink->GetPosition()[1]; // rough direction estimate
5854 if (dird<0) continue;
5856 Float_t dirm = kink->GetMotherP()[0]*kink->GetPosition()[0]+kink->GetMotherP()[1]*kink->GetPosition()[1]; // rough direction estimate
5857 if (dirm<0) continue;
5858 Float_t mpt = TMath::Sqrt(kink->GetMotherP()[0]*kink->GetMotherP()[0]+kink->GetMotherP()[1]*kink->GetMotherP()[1]);
5859 if (mpt<0.2) continue;
5862 //for high momenta momentum not defined well in first iteration
5863 Double_t qt = TMath::Sin(kink->GetAngle(2))*ktrack1->GetP();
5864 if (qt>0.35) continue;
5867 kink->SetLabel(CookLabel(ktrack0,0.4,0,row0),0);
5868 kink->SetLabel(CookLabel(ktrack1,0.4,row0,160),1);
5870 kink->SetTPCDensity(dens00,0,0);
5871 kink->SetTPCDensity(dens01,0,1);
5872 kink->SetTPCDensity(dens10,1,0);
5873 kink->SetTPCDensity(dens11,1,1);
5874 kink->SetIndex(i,0);
5875 kink->SetIndex(j,1);
5878 kink->SetTPCDensity(dens10,0,0);
5879 kink->SetTPCDensity(dens11,0,1);
5880 kink->SetTPCDensity(dens00,1,0);
5881 kink->SetTPCDensity(dens01,1,1);
5882 kink->SetIndex(j,0);
5883 kink->SetIndex(i,1);
5886 if (mpt<1||kink->GetAngle(2)>0.1){
5887 // angle and densities not defined yet
5888 if (kink->GetTPCDensityFactor()<0.8) continue;
5889 if ((2-kink->GetTPCDensityFactor())*kink->GetDistance() >0.25) continue;
5890 if (kink->GetAngle(2)*ktrack0->GetP()<0.003) continue; //too small angle
5891 if (kink->GetAngle(2)>0.2&&kink->GetTPCDensityFactor()<1.15) continue;
5892 if (kink->GetAngle(2)>0.2&&kink->GetTPCDensity(0,1)>0.05) continue;
5894 Float_t criticalangle = track0->GetSigmaSnp2()+track0->GetSigmaTgl2();
5895 criticalangle+= track1->GetSigmaSnp2()+track1->GetSigmaTgl2();
5896 criticalangle= 3*TMath::Sqrt(criticalangle);
5897 if (criticalangle>0.02) criticalangle=0.02;
5898 if (kink->GetAngle(2)<criticalangle) continue;
5901 Int_t drow = Int_t(2.+0.5/(0.05+kink->GetAngle(2))); // overlap region defined
5902 Float_t shapesum =0;
5904 for ( Int_t row = row0-drow; row<row0+drow;row++){
5905 if (row<0) continue;
5906 if (row>155) continue;
5907 if (ktrack0->GetClusterPointer(row)){
5908 AliTPCTrackerPoint *point =ktrack0->GetTrackPoint(row);
5909 shapesum+=point->GetSigmaY()+point->GetSigmaZ();
5912 if (ktrack1->GetClusterPointer(row)){
5913 AliTPCTrackerPoint *point =ktrack1->GetTrackPoint(row);
5914 shapesum+=point->GetSigmaY()+point->GetSigmaZ();
5919 kink->SetShapeFactor(-1.);
5922 kink->SetShapeFactor(shapesum/sum);
5924 // esd->AddKink(kink);
5926 // kink->SetMother(paramm);
5927 //kink->SetDaughter(paramd);
5929 Double_t chi2P2 = paramm.GetParameter()[2]-paramd.GetParameter()[2];
5931 chi2P2/=paramm.GetCovariance()[5]+paramd.GetCovariance()[5];
5932 Double_t chi2P3 = paramm.GetParameter()[3]-paramd.GetParameter()[3];
5934 chi2P3/=paramm.GetCovariance()[9]+paramd.GetCovariance()[9];
5936 if (AliTPCReconstructor::StreamLevel()>1) {
5937 (*fDebugStreamer)<<"kinkLpt"<<
5945 if ( chi2P2+chi2P3<AliTPCReconstructor::GetRecoParam()->GetKinkAngleCutChi2(0)){
5949 kinks->AddLast(kink);
5955 // sort the kinks according quality - and refit them towards vertex
5957 Int_t nkinks = kinks->GetEntriesFast();
5958 Float_t *quality = new Float_t[nkinks];
5959 Int_t *indexes = new Int_t[nkinks];
5960 AliTPCseed **mothers = new AliTPCseed*[nkinks]; memset(mothers, 0, nkinks*sizeof(AliTPCseed*));
5961 AliTPCseed **daughters = new AliTPCseed*[nkinks]; memset(daughters, 0, nkinks*sizeof(AliTPCseed*));
5964 for (Int_t i=0;i<nkinks;i++){
5966 AliKink *kinkl = (AliKink*)kinks->At(i);
5968 // refit kinks towards vertex
5970 Int_t index0 = kinkl->GetIndex(0);
5971 Int_t index1 = kinkl->GetIndex(1);
5972 AliTPCseed * ktrack0 = (AliTPCseed*)array->At(index0);
5973 AliTPCseed * ktrack1 = (AliTPCseed*)array->At(index1);
5975 Int_t sumn=ktrack0->GetNumberOfClusters()+ktrack1->GetNumberOfClusters();
5977 // Refit Kink under if too small angle
5979 if (kinkl->GetAngle(2)<0.05){
5980 kinkl->SetTPCRow0(GetRowNumber(kinkl->GetR()));
5981 Int_t row0 = kinkl->GetTPCRow0();
5982 Int_t drow = Int_t(2.+0.5/(0.05+kinkl->GetAngle(2)));
5985 Int_t last = row0-drow;
5986 if (last<40) last=40;
5987 if (last<ktrack0->GetFirstPoint()+25) last = ktrack0->GetFirstPoint()+25;
5988 AliTPCseed* seed0 = ReSeed(ktrack0,last,kFALSE);
5991 Int_t first = row0+drow;
5992 if (first>130) first=130;
5993 if (first>ktrack1->GetLastPoint()-25) first = TMath::Max(ktrack1->GetLastPoint()-25,30);
5994 AliTPCseed* seed1 = ReSeed(ktrack1,first,kTRUE);
5996 if (seed0 && seed1){
5997 kinkl->SetStatus(1,8);
5998 if (RefitKink(*seed0,*seed1,*kinkl)) kinkl->SetStatus(1,9);
5999 row0 = GetRowNumber(kinkl->GetR());
6000 sumn = seed0->GetNumberOfClusters()+seed1->GetNumberOfClusters();
6001 mothers[i] = new ( NextFreeSeed() ) AliTPCseed(*seed0);
6002 mothers[i]->SetPoolID(fLastSeedID);
6003 daughters[i] = new (NextFreeSeed() ) AliTPCseed(*seed1);
6004 daughters[i]->SetPoolID(fLastSeedID);
6007 delete kinks->RemoveAt(i);
6008 if (seed0) MarkSeedFree( seed0 );
6009 if (seed1) MarkSeedFree( seed1 );
6012 if (kinkl->GetDistance()>0.5 || kinkl->GetR()<110 || kinkl->GetR()>240) {
6013 delete kinks->RemoveAt(i);
6014 if (seed0) MarkSeedFree( seed0 );
6015 if (seed1) MarkSeedFree( seed1 );
6019 MarkSeedFree( seed0 );
6020 MarkSeedFree( seed1 );
6023 if (kinkl) quality[i] = 160*((0.1+kinkl->GetDistance())*(2.-kinkl->GetTPCDensityFactor()))/(sumn+40.); //the longest -clossest will win
6025 TMath::Sort(nkinks,quality,indexes,kFALSE);
6027 //remove double find kinks
6029 for (Int_t ikink0=1;ikink0<nkinks;ikink0++){
6030 AliKink * kink0 = (AliKink*) kinks->At(indexes[ikink0]);
6031 if (!kink0) continue;
6033 for (Int_t ikink1=0;ikink1<ikink0;ikink1++){
6034 kink0 = (AliKink*) kinks->At(indexes[ikink0]);
6035 if (!kink0) continue;
6036 AliKink * kink1 = (AliKink*) kinks->At(indexes[ikink1]);
6037 if (!kink1) continue;
6038 // if not close kink continue
6039 if (TMath::Abs(kink1->GetPosition()[2]-kink0->GetPosition()[2])>10) continue;
6040 if (TMath::Abs(kink1->GetPosition()[1]-kink0->GetPosition()[1])>10) continue;
6041 if (TMath::Abs(kink1->GetPosition()[0]-kink0->GetPosition()[0])>10) continue;
6043 AliTPCseed &mother0 = *mothers[indexes[ikink0]];
6044 AliTPCseed &daughter0 = *daughters[indexes[ikink0]];
6045 AliTPCseed &mother1 = *mothers[indexes[ikink1]];
6046 AliTPCseed &daughter1 = *daughters[indexes[ikink1]];
6047 Int_t row0 = (kink0->GetTPCRow0()+kink1->GetTPCRow0())/2;
6056 for (Int_t i=0;i<row0;i++){
6057 if (mother0.GetClusterIndex(i)>0 && mother1.GetClusterIndex(i)>0){
6060 if (mother0.GetClusterIndex(i)==mother1.GetClusterIndex(i)){
6067 for (Int_t i=row0;i<158;i++){
6068 //if (daughter0.GetClusterIndex(i)>0 && daughter0.GetClusterIndex(i)>0){ // RS: Bug?
6069 if (daughter0.GetClusterIndex(i)>0 && daughter1.GetClusterIndex(i)>0){
6072 if (mother0.GetClusterIndex(i)==mother1.GetClusterIndex(i)){
6078 Float_t ratio = Float_t(same+1)/Float_t(both+1);
6079 Float_t ratiom = Float_t(samem+1)/Float_t(bothm+1);
6080 Float_t ratiod = Float_t(samed+1)/Float_t(bothd+1);
6081 if (ratio>0.3 && ratiom>0.5 &&ratiod>0.5) {
6082 Int_t sum0 = mother0.GetNumberOfClusters()+daughter0.GetNumberOfClusters();
6083 Int_t sum1 = mother1.GetNumberOfClusters()+daughter1.GetNumberOfClusters();
6085 shared[kink0->GetIndex(0)]= kTRUE;
6086 shared[kink0->GetIndex(1)]= kTRUE;
6087 delete kinks->RemoveAt(indexes[ikink0]);
6091 shared[kink1->GetIndex(0)]= kTRUE;
6092 shared[kink1->GetIndex(1)]= kTRUE;
6093 delete kinks->RemoveAt(indexes[ikink1]);
6100 for (Int_t i=0;i<nkinks;i++){
6101 AliKink * kinkl = (AliKink*) kinks->At(indexes[i]);
6102 if (!kinkl) continue;
6103 kinkl->SetTPCRow0(GetRowNumber(kinkl->GetR()));
6104 Int_t index0 = kinkl->GetIndex(0);
6105 Int_t index1 = kinkl->GetIndex(1);
6106 if (circular[index0]||(circular[index1]&&kinkl->GetDistance()>0.2)) continue;
6107 kinkl->SetMultiple(usage[index0],0);
6108 kinkl->SetMultiple(usage[index1],1);
6109 if (kinkl->GetMultiple()[0]+kinkl->GetMultiple()[1]>2) continue;
6110 if (kinkl->GetMultiple()[0]+kinkl->GetMultiple()[1]>0 && quality[indexes[i]]>0.2) continue;
6111 if (kinkl->GetMultiple()[0]+kinkl->GetMultiple()[1]>0 && kinkl->GetDistance()>0.2) continue;
6112 if (circular[index0]||(circular[index1]&&kinkl->GetDistance()>0.1)) continue;
6114 AliTPCseed * ktrack0 = (AliTPCseed*)array->At(index0);
6115 AliTPCseed * ktrack1 = (AliTPCseed*)array->At(index1);
6116 if (!ktrack0 || !ktrack1) continue;
6117 Int_t index = esd->AddKink(kinkl);
6120 if ( ktrack0->GetKinkIndex(0)==0 && ktrack1->GetKinkIndex(0)==0) { //best kink
6121 if (mothers[indexes[i]]->GetNumberOfClusters()>20 && daughters[indexes[i]]->GetNumberOfClusters()>20 &&
6122 (mothers[indexes[i]]->GetNumberOfClusters()+daughters[indexes[i]]->GetNumberOfClusters())>100){
6123 *ktrack0 = *mothers[indexes[i]];
6124 *ktrack1 = *daughters[indexes[i]];
6128 ktrack0->SetKinkIndex(usage[index0],-(index+1));
6129 ktrack1->SetKinkIndex(usage[index1], (index+1));
6134 // Remove tracks corresponding to shared kink's
6136 for (Int_t i=0;i<nentries;i++){
6137 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
6138 if (!track0) continue;
6139 if (track0->GetKinkIndex(0)!=0) continue;
6140 if (shared[i]) MarkSeedFree( array->RemoveAt(i) );
6145 RemoveUsed2(array,0.5,0.4,30);
6147 for (Int_t i=0;i<nentries;i++){
6148 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
6149 if (!track0) continue;
6150 track0->CookdEdx(0.02,0.6);
6154 for (Int_t i=0;i<nentries;i++){
6155 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
6156 if (!track0) continue;
6157 if (track0->Pt()<1.4) continue;
6158 //remove double high momenta tracks - overlapped with kink candidates
6161 for (Int_t icl=track0->GetFirstPoint();icl<track0->GetLastPoint(); icl++){
6162 if (track0->GetClusterPointer(icl)!=0){
6164 if (track0->GetClusterPointer(icl)->IsUsed(10)) ishared++;
6167 if (Float_t(ishared+1)/Float_t(all+1)>0.5) {
6168 MarkSeedFree( array->RemoveAt(i) );
6172 if (track0->GetKinkIndex(0)!=0) continue;
6173 if (track0->GetNumberOfClusters()<80) continue;
6175 AliTPCseed *pmother = new( NextFreeSeed() ) AliTPCseed();
6176 pmother->SetPoolID(fLastSeedID);
6177 AliTPCseed *pdaughter = new( NextFreeSeed() ) AliTPCseed();
6178 pdaughter->SetPoolID(fLastSeedID);
6179 AliKink *pkink = new AliKink;
6181 AliTPCseed & mother = *pmother;
6182 AliTPCseed & daughter = *pdaughter;
6183 AliKink & kinkl = *pkink;
6184 if (CheckKinkPoint(track0,mother,daughter, kinkl)){
6185 if (mother.GetNumberOfClusters()<30||daughter.GetNumberOfClusters()<20) {
6186 MarkSeedFree( pmother );
6187 MarkSeedFree( pdaughter );
6189 continue; //too short tracks
6191 if (mother.Pt()<1.4) {
6192 MarkSeedFree( pmother );
6193 MarkSeedFree( pdaughter );
6197 Int_t row0= kinkl.GetTPCRow0();
6198 if (kinkl.GetDistance()>0.5 || kinkl.GetR()<110. || kinkl.GetR()>240.) {
6199 MarkSeedFree( pmother );
6200 MarkSeedFree( pdaughter );
6205 Int_t index = esd->AddKink(&kinkl);
6206 mother.SetKinkIndex(0,-(index+1));
6207 daughter.SetKinkIndex(0,index+1);
6208 if (mother.GetNumberOfClusters()>50) {
6209 MarkSeedFree( array->RemoveAt(i) );
6210 AliTPCseed* mtc = new( NextFreeSeed() ) AliTPCseed(mother);
6211 mtc->SetPoolID(fLastSeedID);
6212 array->AddAt(mtc,i);
6215 AliTPCseed* mtc = new( NextFreeSeed() ) AliTPCseed(mother);
6216 mtc->SetPoolID(fLastSeedID);
6217 array->AddLast(mtc);
6219 AliTPCseed* dtc = new( NextFreeSeed() ) AliTPCseed(daughter);
6220 dtc->SetPoolID(fLastSeedID);
6221 array->AddLast(dtc);
6222 for (Int_t icl=0;icl<row0;icl++) {
6223 if (mother.GetClusterPointer(icl)) mother.GetClusterPointer(icl)->Use(20);
6226 for (Int_t icl=row0;icl<158;icl++) {
6227 if (daughter.GetClusterPointer(icl)) daughter.GetClusterPointer(icl)->Use(20);
6231 MarkSeedFree( pmother );
6232 MarkSeedFree( pdaughter );
6236 delete [] daughters;
6258 AliInfo(Form("Ncandidates=\t%d\t%d\t%d\t%d\n",esd->GetNumberOfKinks(),ncandidates,ntracks,nall));
6263 Int_t AliTPCtrackerMI::RefitKink(AliTPCseed &mother, AliTPCseed &daughter, const AliESDkink &knk)
6266 // refit kink towards to the vertex
6269 AliKink &kink=(AliKink &)knk;
6271 Int_t row0 = GetRowNumber(kink.GetR());
6272 FollowProlongation(mother,0);
6273 mother.Reset(kFALSE);
6275 FollowProlongation(daughter,row0);
6276 daughter.Reset(kFALSE);
6277 FollowBackProlongation(daughter,158);
6278 daughter.Reset(kFALSE);
6279 Int_t first = TMath::Max(row0-20,30);
6280 Int_t last = TMath::Min(row0+20,140);
6282 const Int_t kNdiv =5;
6283 AliTPCseed param0[kNdiv]; // parameters along the track
6284 AliTPCseed param1[kNdiv]; // parameters along the track
6285 AliKink kinks[kNdiv]; // corresponding kink parameters
6288 for (Int_t irow=0; irow<kNdiv;irow++){
6289 rows[irow] = first +((last-first)*irow)/(kNdiv-1);
6291 // store parameters along the track
6293 for (Int_t irow=0;irow<kNdiv;irow++){
6294 FollowBackProlongation(mother, rows[irow]);
6295 FollowProlongation(daughter,rows[kNdiv-1-irow]);
6296 param0[irow] = mother;
6297 param1[kNdiv-1-irow] = daughter;
6301 for (Int_t irow=0; irow<kNdiv-1;irow++){
6302 if (param0[irow].GetNumberOfClusters()<kNdiv||param1[irow].GetNumberOfClusters()<kNdiv) continue;
6303 kinks[irow].SetMother(param0[irow]);
6304 kinks[irow].SetDaughter(param1[irow]);
6305 kinks[irow].Update();
6308 // choose kink with best "quality"
6310 Double_t mindist = 10000;
6311 for (Int_t irow=0;irow<kNdiv;irow++){
6312 if (param0[irow].GetNumberOfClusters()<20||param1[irow].GetNumberOfClusters()<20) continue;
6313 if (TMath::Abs(kinks[irow].GetR())>240.) continue;
6314 if (TMath::Abs(kinks[irow].GetR())<100.) continue;
6316 Float_t normdist = TMath::Abs(param0[irow].GetX()-kinks[irow].GetR())*(0.1+kink.GetDistance());
6317 normdist/= (param0[irow].GetNumberOfClusters()+param1[irow].GetNumberOfClusters()+40.);
6318 if (normdist < mindist){
6324 if (index==-1) return 0;
6327 param0[index].Reset(kTRUE);
6328 FollowProlongation(param0[index],0);
6330 mother = param0[index];
6331 daughter = param1[index]; // daughter in vertex
6333 kink.SetMother(mother);
6334 kink.SetDaughter(daughter);
6336 kink.SetTPCRow0(GetRowNumber(kink.GetR()));
6337 kink.SetTPCncls(param0[index].GetNumberOfClusters(),0);
6338 kink.SetTPCncls(param1[index].GetNumberOfClusters(),1);
6339 kink.SetLabel(CookLabel(&mother,0.4, 0,kink.GetTPCRow0()),0);
6340 kink.SetLabel(CookLabel(&daughter,0.4, kink.GetTPCRow0(),160),1);
6341 mother.SetLabel(kink.GetLabel(0));
6342 daughter.SetLabel(kink.GetLabel(1));
6348 void AliTPCtrackerMI::UpdateKinkQualityM(AliTPCseed * seed){
6350 // update Kink quality information for mother after back propagation
6352 if (seed->GetKinkIndex(0)>=0) return;
6353 for (Int_t ikink=0;ikink<3;ikink++){
6354 Int_t index = seed->GetKinkIndex(ikink);
6355 if (index>=0) break;
6356 index = TMath::Abs(index)-1;
6357 AliESDkink * kink = fEvent->GetKink(index);
6358 kink->SetTPCDensity(-1,0,0);
6359 kink->SetTPCDensity(1,0,1);
6361 Int_t row0 = kink->GetTPCRow0() - 2 - Int_t( 0.5/ (0.05+kink->GetAngle(2)));
6362 if (row0<15) row0=15;
6364 Int_t row1 = kink->GetTPCRow0() + 2 + Int_t( 0.5/ (0.05+kink->GetAngle(2)));
6365 if (row1>145) row1=145;
6367 Int_t found,foundable,shared;
6368 seed->GetClusterStatistic(0,row0, found, foundable,shared,kFALSE);
6369 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),0,0);
6370 seed->GetClusterStatistic(row1,155, found, foundable,shared,kFALSE);
6371 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),0,1);
6376 void AliTPCtrackerMI::UpdateKinkQualityD(AliTPCseed * seed){
6378 // update Kink quality information for daughter after refit
6380 if (seed->GetKinkIndex(0)<=0) return;
6381 for (Int_t ikink=0;ikink<3;ikink++){
6382 Int_t index = seed->GetKinkIndex(ikink);
6383 if (index<=0) break;
6384 index = TMath::Abs(index)-1;
6385 AliESDkink * kink = fEvent->GetKink(index);
6386 kink->SetTPCDensity(-1,1,0);
6387 kink->SetTPCDensity(-1,1,1);
6389 Int_t row0 = kink->GetTPCRow0() -2 - Int_t( 0.5/ (0.05+kink->GetAngle(2)));
6390 if (row0<15) row0=15;
6392 Int_t row1 = kink->GetTPCRow0() +2 + Int_t( 0.5/ (0.05+kink->GetAngle(2)));
6393 if (row1>145) row1=145;
6395 Int_t found,foundable,shared;
6396 seed->GetClusterStatistic(0,row0, found, foundable,shared,kFALSE);
6397 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),1,0);
6398 seed->GetClusterStatistic(row1,155, found, foundable,shared,kFALSE);
6399 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),1,1);
6405 Int_t AliTPCtrackerMI::CheckKinkPoint(AliTPCseed*seed,AliTPCseed &mother, AliTPCseed &daughter, const AliESDkink &knk)
6408 // check kink point for given track
6409 // if return value=0 kink point not found
6410 // otherwise seed0 correspond to mother particle
6411 // seed1 correspond to daughter particle
6412 // kink parameter of kink point
6413 AliKink &kink=(AliKink &)knk;
6415 Int_t middlerow = (seed->GetFirstPoint()+seed->GetLastPoint())/2;
6416 Int_t first = seed->GetFirstPoint();
6417 Int_t last = seed->GetLastPoint();
6418 if (last-first<20) return 0; // shortest length - 2*30 = 60 pad-rows
6421 AliTPCseed *seed1 = ReSeed(seed,middlerow+20, kTRUE); //middle of chamber
6422 if (!seed1) return 0;
6423 FollowProlongation(*seed1,seed->GetLastPoint()-20);
6424 seed1->Reset(kTRUE);
6425 FollowProlongation(*seed1,158);
6426 seed1->Reset(kTRUE);
6427 last = seed1->GetLastPoint();
6429 AliTPCseed *seed0 = new( NextFreeSeed() ) AliTPCseed(*seed);
6430 seed0->SetPoolID(fLastSeedID);
6431 seed0->Reset(kFALSE);
6434 AliTPCseed param0[20]; // parameters along the track
6435 AliTPCseed param1[20]; // parameters along the track
6436 AliKink kinks[20]; // corresponding kink parameters
6438 for (Int_t irow=0; irow<20;irow++){
6439 rows[irow] = first +((last-first)*irow)/19;
6441 // store parameters along the track
6443 for (Int_t irow=0;irow<20;irow++){
6444 FollowBackProlongation(*seed0, rows[irow]);
6445 FollowProlongation(*seed1,rows[19-irow]);
6446 param0[irow] = *seed0;
6447 param1[19-irow] = *seed1;
6451 for (Int_t irow=0; irow<19;irow++){
6452 kinks[irow].SetMother(param0[irow]);
6453 kinks[irow].SetDaughter(param1[irow]);
6454 kinks[irow].Update();
6457 // choose kink with biggest change of angle
6459 Double_t maxchange= 0;
6460 for (Int_t irow=1;irow<19;irow++){
6461 if (TMath::Abs(kinks[irow].GetR())>240.) continue;
6462 if (TMath::Abs(kinks[irow].GetR())<110.) continue;
6463 Float_t quality = TMath::Abs(kinks[irow].GetAngle(2))/(3.+TMath::Abs(kinks[irow].GetR()-param0[irow].GetX()));
6464 if ( quality > maxchange){
6465 maxchange = quality;
6470 MarkSeedFree( seed0 );
6471 MarkSeedFree( seed1 );
6472 if (index<0) return 0;
6474 Int_t row0 = GetRowNumber(kinks[index].GetR()); //row 0 estimate
6475 seed0 = new( NextFreeSeed() ) AliTPCseed(param0[index]);
6476 seed0->SetPoolID(fLastSeedID);
6477 seed1 = new( NextFreeSeed() ) AliTPCseed(param1[index]);
6478 seed1->SetPoolID(fLastSeedID);
6479 seed0->Reset(kFALSE);
6480 seed1->Reset(kFALSE);
6481 seed0->ResetCovariance(10.);
6482 seed1->ResetCovariance(10.);
6483 FollowProlongation(*seed0,0);
6484 FollowBackProlongation(*seed1,158);
6485 mother = *seed0; // backup mother at position 0
6486 seed0->Reset(kFALSE);
6487 seed1->Reset(kFALSE);
6488 seed0->ResetCovariance(10.);
6489 seed1->ResetCovariance(10.);
6491 first = TMath::Max(row0-20,0);
6492 last = TMath::Min(row0+20,158);
6494 for (Int_t irow=0; irow<20;irow++){
6495 rows[irow] = first +((last-first)*irow)/19;
6497 // store parameters along the track
6499 for (Int_t irow=0;irow<20;irow++){
6500 FollowBackProlongation(*seed0, rows[irow]);
6501 FollowProlongation(*seed1,rows[19-irow]);
6502 param0[irow] = *seed0;
6503 param1[19-irow] = *seed1;
6507 for (Int_t irow=0; irow<19;irow++){
6508 kinks[irow].SetMother(param0[irow]);
6509 kinks[irow].SetDaughter(param1[irow]);
6510 // param0[irow].Dump();
6511 //param1[irow].Dump();
6512 kinks[irow].Update();
6515 // choose kink with biggest change of angle
6518 for (Int_t irow=0;irow<20;irow++){
6519 if (TMath::Abs(kinks[irow].GetR())>250.) continue;
6520 if (TMath::Abs(kinks[irow].GetR())<90.) continue;
6521 Float_t quality = TMath::Abs(kinks[irow].GetAngle(2))/(3.+TMath::Abs(kinks[irow].GetR()-param0[irow].GetX()));
6522 if ( quality > maxchange){
6523 maxchange = quality;
6530 if (index==-1 || param0[index].GetNumberOfClusters()+param1[index].GetNumberOfClusters()<100){
6531 MarkSeedFree( seed0 );
6532 MarkSeedFree( seed1 );
6536 // Float_t anglesigma = TMath::Sqrt(param0[index].fC22+param0[index].fC33+param1[index].fC22+param1[index].fC33);
6538 kink.SetMother(param0[index]);
6539 kink.SetDaughter(param1[index]);
6542 Double_t chi2P2 = param0[index].GetParameter()[2]-param1[index].GetParameter()[2];
6544 chi2P2/=param0[index].GetCovariance()[5]+param1[index].GetCovariance()[5];
6545 Double_t chi2P3 = param0[index].GetParameter()[3]-param1[index].GetParameter()[3];
6547 chi2P3/=param0[index].GetCovariance()[9]+param1[index].GetCovariance()[9];
6549 if (AliTPCReconstructor::StreamLevel()>1) {
6550 (*fDebugStreamer)<<"kinkHpt"<<
6553 "p0.="<<¶m0[index]<<
6554 "p1.="<<¶m1[index]<<
6558 if ( chi2P2+chi2P3<AliTPCReconstructor::GetRecoParam()->GetKinkAngleCutChi2(0)){
6559 MarkSeedFree( seed0 );
6560 MarkSeedFree( seed1 );
6565 row0 = GetRowNumber(kink.GetR());
6566 kink.SetTPCRow0(row0);
6567 kink.SetLabel(CookLabel(seed0,0.5,0,row0),0);
6568 kink.SetLabel(CookLabel(seed1,0.5,row0,158),1);
6569 kink.SetIndex(-10,0);
6570 kink.SetIndex(int(param0[index].GetNumberOfClusters()+param1[index].GetNumberOfClusters()),1);
6571 kink.SetTPCncls(param0[index].GetNumberOfClusters(),0);
6572 kink.SetTPCncls(param1[index].GetNumberOfClusters(),1);
6575 // new (&mother) AliTPCseed(param0[index]);
6576 daughter = param1[index];
6577 daughter.SetLabel(kink.GetLabel(1));
6578 param0[index].Reset(kTRUE);
6579 FollowProlongation(param0[index],0);
6580 mother = param0[index];
6581 mother.SetLabel(kink.GetLabel(0));
6582 if ( chi2P2+chi2P3<AliTPCReconstructor::GetRecoParam()->GetKinkAngleCutChi2(1)){
6585 MarkSeedFree( seed0 );
6586 MarkSeedFree( seed1 );
6594 AliTPCseed* AliTPCtrackerMI::ReSeed(AliTPCseed *t)
6597 // reseed - refit - track
6600 // Int_t last = fSectors->GetNRows()-1;
6602 if (fSectors == fOuterSec){
6603 first = TMath::Max(first, t->GetFirstPoint()-fInnerSec->GetNRows());
6607 first = t->GetFirstPoint();
6609 AliTPCseed * seed = MakeSeed(t,0.1,0.5,0.9);
6610 FollowBackProlongation(*t,fSectors->GetNRows()-1);
6612 FollowProlongation(*t,first);
6622 //_____________________________________________________________________________
6623 Int_t AliTPCtrackerMI::ReadSeeds(const TFile *inp) {
6624 //-----------------------------------------------------------------
6625 // This function reades track seeds.
6626 //-----------------------------------------------------------------
6627 TDirectory *savedir=gDirectory;
6629 TFile *in=(TFile*)inp;
6630 if (!in->IsOpen()) {
6631 cerr<<"AliTPCtrackerMI::ReadSeeds(): input file is not open !\n";
6636 TTree *seedTree=(TTree*)in->Get("Seeds");
6638 cerr<<"AliTPCtrackerMI::ReadSeeds(): ";
6639 cerr<<"can't get a tree with track seeds !\n";
6642 AliTPCtrack *seed=new AliTPCtrack;
6643 seedTree->SetBranchAddress("tracks",&seed);
6645 if (fSeeds==0) fSeeds=new TObjArray(15000);
6647 Int_t n=(Int_t)seedTree->GetEntries();
6648 for (Int_t i=0; i<n; i++) {
6649 seedTree->GetEvent(i);
6650 AliTPCseed* sdc = new( NextFreeSeed() ) AliTPCseed(*seed/*,seed->GetAlpha()*/);
6651 sdc->SetPoolID(fLastSeedID);
6652 fSeeds->AddLast(sdc);
6655 delete seed; // RS: this seed is not from the pool, delete it !!!
6661 Int_t AliTPCtrackerMI::Clusters2Tracks (AliESDEvent *const esd)
6664 // clusters to tracks
6665 if (fSeeds) DeleteSeeds();
6666 else ResetSeedsPool();
6669 AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ;
6670 transform->SetCurrentTimeStamp( esd->GetTimeStamp());
6671 transform->SetCurrentRun(esd->GetRunNumber());
6674 if (!fSeeds) return 1;
6676 if (AliTPCReconstructor::StreamLevel()>3) DumpClusters(0,fSeeds);
6682 //_____________________________________________________________________________
6683 Int_t AliTPCtrackerMI::Clusters2Tracks() {
6684 //-----------------------------------------------------------------
6685 // This is a track finder.
6686 //-----------------------------------------------------------------
6687 TDirectory *savedir=gDirectory;
6691 fSeeds = Tracking();
6694 Info("Clusters2Tracks","Time for tracking: \t");timer.Print();timer.Start();
6696 //activate again some tracks
6697 for (Int_t i=0; i<fSeeds->GetEntriesFast(); i++) {
6698 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
6700 Int_t nc=t.GetNumberOfClusters();
6702 MarkSeedFree( fSeeds->RemoveAt(i) );
6706 if (pt->GetRemoval()==10) {
6707 if (pt->GetDensityFirst(20)>0.8 || pt->GetDensityFirst(30)>0.8 || pt->GetDensityFirst(40)>0.7)
6708 pt->Desactivate(10); // make track again active // MvL: should be 0 ?
6710 pt->Desactivate(20);
6711 MarkSeedFree( fSeeds->RemoveAt(i) );
6716 RemoveUsed2(fSeeds,0.85,0.85,0);
6717 if (AliTPCReconstructor::GetRecoParam()->GetDoKinks()) FindKinks(fSeeds,fEvent);
6718 //FindCurling(fSeeds, fEvent,0);
6719 if (AliTPCReconstructor::StreamLevel()>5) FindMultiMC(fSeeds, fEvent,-1); // find multi found tracks
6720 RemoveUsed2(fSeeds,0.5,0.4,20);
6721 FindSplitted(fSeeds, fEvent,0); // find multi found tracks
6722 if (AliTPCReconstructor::StreamLevel()>5) FindMultiMC(fSeeds, fEvent,0); // find multi found tracks
6725 // // refit short tracks
6727 Int_t nseed=fSeeds->GetEntriesFast();
6730 for (Int_t i=0; i<nseed; i++) {
6731 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
6733 Int_t nc=t.GetNumberOfClusters();
6735 MarkSeedFree( fSeeds->RemoveAt(i) );
6738 CookLabel(pt,0.1); //For comparison only
6739 //if ((pt->IsActive() || (pt->fRemoval==10) )&& nc>50 &&pt->GetNumberOfClusters()>0.4*pt->fNFoundable){
6740 if ((pt->IsActive() || (pt->GetRemoval()==10) )){
6742 if (fDebug>0) cerr<<found<<'\r';
6746 MarkSeedFree( fSeeds->RemoveAt(i) );
6750 //RemoveOverlap(fSeeds,0.99,7,kTRUE);
6752 //RemoveUsed(fSeeds,0.9,0.9,6);
6754 nseed=fSeeds->GetEntriesFast();
6756 for (Int_t i=0; i<nseed; i++) {
6757 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
6759 Int_t nc=t.GetNumberOfClusters();
6761 MarkSeedFree( fSeeds->RemoveAt(i) );
6765 t.CookdEdx(0.02,0.6);
6766 // CheckKinkPoint(&t,0.05);
6767 //if ((pt->IsActive() || (pt->fRemoval==10) )&& nc>50 &&pt->GetNumberOfClusters()>0.4*pt->fNFoundable){
6768 if ((pt->IsActive() || (pt->GetRemoval()==10) )){
6776 MarkSeedFree( fSeeds->RemoveAt(i) );
6777 //AliTPCseed * seed1 = ReSeed(pt,0.05,0.5,1);
6779 // FollowProlongation(*seed1,0);
6780 // Int_t n = seed1->GetNumberOfClusters();
6781 // printf("fP4\t%f\t%f\n",seed1->GetC(),pt->GetC());
6782 // printf("fN\t%d\t%d\n", seed1->GetNumberOfClusters(),pt->GetNumberOfClusters());
6785 //AliTPCseed * seed2 = ReSeed(pt,0.95,0.5,0.05);
6789 SortTracks(fSeeds, 1);
6793 PrepareForBackProlongation(fSeeds,5.);
6794 PropagateBack(fSeeds);
6795 printf("Time for back propagation: \t");timer.Print();timer.Start();
6799 PrepareForProlongation(fSeeds,5.);
6800 PropagateForard2(fSeeds);
6802 printf("Time for FORWARD propagation: \t");timer.Print();timer.Start();
6803 // RemoveUsed(fSeeds,0.7,0.7,6);
6804 //RemoveOverlap(fSeeds,0.9,7,kTRUE);
6806 nseed=fSeeds->GetEntriesFast();
6808 for (Int_t i=0; i<nseed; i++) {
6809 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
6811 Int_t nc=t.GetNumberOfClusters();
6813 MarkSeedFree( fSeeds->RemoveAt(i) );
6816 t.CookdEdx(0.02,0.6);
6817 // CookLabel(pt,0.1); //For comparison only
6818 //if ((pt->IsActive() || (pt->fRemoval==10) )&& nc>50 &&pt->GetNumberOfClusters()>0.4*pt->fNFoundable){
6819 if ((pt->IsActive() || (pt->fRemoval==10) )){
6820 cerr<<found++<<'\r';
6823 MarkSeedFree( fSeeds->RemoveAt(i) );
6828 // fNTracks = found;
6830 Info("Clusters2Tracks","Time for overlap removal, track writing and dedx cooking: \t"); timer.Print();timer.Start();
6833 // cerr<<"Number of found tracks : "<<"\t"<<found<<endl;
6834 Info("Clusters2Tracks","Number of found tracks %d",found);
6836 // UnloadClusters();
6841 void AliTPCtrackerMI::Tracking(TObjArray * arr)
6844 // tracking of the seeds
6847 fSectors = fOuterSec;
6848 ParallelTracking(arr,150,63);
6849 fSectors = fOuterSec;
6850 ParallelTracking(arr,63,0);
6853 TObjArray * AliTPCtrackerMI::Tracking(Int_t seedtype, Int_t i1, Int_t i2, Float_t cuts[4], Float_t dy, Int_t dsec)
6858 static TObjArray arrTracks;
6859 TObjArray * arr = &arrTracks;
6861 fSectors = fOuterSec;
6864 for (Int_t sec=0;sec<fkNOS;sec++){
6865 if (seedtype==3) MakeSeeds3(arr,sec,i1,i2,cuts,dy, dsec);
6866 if (seedtype==4) MakeSeeds5(arr,sec,i1,i2,cuts,dy);
6867 if (seedtype==2) MakeSeeds2(arr,sec,i1,i2,cuts,dy);
6870 Info("Tracking","\nSeeding - %d\t%d\t%d\t%d\n",seedtype,i1,i2,arr->GetEntriesFast());
6882 TObjArray * AliTPCtrackerMI::Tracking()
6886 if (AliTPCReconstructor::GetRecoParam()->GetSpecialSeeding()) return TrackingSpecial();
6889 Int_t nup=fOuterSec->GetNRows()+fInnerSec->GetNRows();
6891 TObjArray * seeds = new TObjArray;
6893 Int_t fLastSeedRowSec=AliTPCReconstructor::GetRecoParam()->GetLastSeedRowSec();
6894 Int_t gapPrim = AliTPCReconstructor::GetRecoParam()->GetSeedGapPrim();
6895 Int_t gapSec = AliTPCReconstructor::GetRecoParam()->GetSeedGapSec();
6903 Float_t fnumber = 3.0;
6904 Float_t fdensity = 3.0;
6909 for (Int_t delta = 0; delta<18; delta+=gapPrim){
6913 arr = Tracking(3,nup-1-delta,nup-1-delta-gap,cuts,-1,1);
6914 SumTracks(seeds,arr);
6915 SignClusters(seeds,fnumber,fdensity);
6917 for (Int_t i=2;i<6;i+=2){
6918 // seed high pt tracks
6921 arr = Tracking(3,nup-i-delta,nup-i-delta-gap,cuts,-1,0);
6922 SumTracks(seeds,arr);
6923 SignClusters(seeds,fnumber,fdensity);
6928 // RemoveUsed(seeds,0.9,0.9,1);
6929 // UnsignClusters();
6930 // SignClusters(seeds,fnumber,fdensity);
6934 for (Int_t delta = 20; delta<120; delta+=gapPrim){
6936 // seed high pt tracks
6940 arr = Tracking(3,nup-delta,nup-delta-gap,cuts,-1);
6941 SumTracks(seeds,arr);
6942 SignClusters(seeds,fnumber,fdensity);
6947 arr = Tracking(3,nup-delta-5,nup-delta-5-gap,cuts,-1);
6948 SumTracks(seeds,arr);
6949 SignClusters(seeds,fnumber,fdensity);
6960 Info("Tracking()","\n\nPrimary seeding\t%d\n\n",seeds->GetEntriesFast());
6964 // RemoveUsed(seeds,0.75,0.75,1);
6966 //SignClusters(seeds,fnumber,fdensity);
6975 arr = Tracking(4,nup-1,nup-1-gap,cuts,-1);
6976 SumTracks(seeds,arr);
6977 SignClusters(seeds,fnumber,fdensity);
6979 arr = Tracking(4,nup-2,nup-2-gap,cuts,-1);
6980 SumTracks(seeds,arr);
6981 SignClusters(seeds,fnumber,fdensity);
6983 arr = Tracking(4,nup-3,nup-3-gap,cuts,-1);
6984 SumTracks(seeds,arr);
6985 SignClusters(seeds,fnumber,fdensity);
6987 arr = Tracking(4,nup-5,nup-5-gap,cuts,-1);
6988 SumTracks(seeds,arr);
6989 SignClusters(seeds,fnumber,fdensity);
6991 arr = Tracking(4,nup-7,nup-7-gap,cuts,-1);
6992 SumTracks(seeds,arr);
6993 SignClusters(seeds,fnumber,fdensity);
6996 arr = Tracking(4,nup-9,nup-9-gap,cuts,-1);
6997 SumTracks(seeds,arr);
6998 SignClusters(seeds,fnumber,fdensity);
7002 for (Int_t delta = 9; delta<30; delta+=gapSec){
7008 arr = Tracking(4,nup-1-delta,nup-1-delta-gap,cuts,-1);
7009 SumTracks(seeds,arr);
7010 SignClusters(seeds,fnumber,fdensity);
7012 arr = Tracking(4,nup-3-delta,nup-5-delta-gap,cuts,4);
7013 SumTracks(seeds,arr);
7014 SignClusters(seeds,fnumber,fdensity);
7027 for (Int_t delta = 30; delta<fLastSeedRowSec; delta+=gapSec){
7033 arr = Tracking(4,nup-1-delta,nup-1-delta-gap,cuts,-1);
7034 SumTracks(seeds,arr);
7035 SignClusters(seeds,fnumber,fdensity);
7037 arr = Tracking(4,nup-5-delta,nup-5-delta-gap,cuts,5 );
7038 SumTracks(seeds,arr);
7039 SignClusters(seeds,fnumber,fdensity);
7043 Info("Tracking()","\n\nSecondary seeding\t%d\n\n",seeds->GetEntriesFast());
7054 TObjArray * AliTPCtrackerMI::TrackingSpecial()
7057 // seeding adjusted for laser and cosmic tests - short tracks with big inclination angle
7058 // no primary vertex seeding tried
7062 Int_t nup=fOuterSec->GetNRows()+fInnerSec->GetNRows();
7064 TObjArray * seeds = new TObjArray;
7069 Float_t fnumber = 3.0;
7070 Float_t fdensity = 3.0;
7073 cuts[0] = AliTPCReconstructor::GetRecoParam()->GetMaxC(); // max curvature
7074 cuts[1] = 3.5; // max tan(phi) angle for seeding
7075 cuts[2] = 3.; // not used (cut on z primary vertex)
7076 cuts[3] = 3.5; // max tan(theta) angle for seeding
7078 for (Int_t delta = 0; nup-delta-gap-1>0; delta+=3){
7080 arr = Tracking(4,nup-1-delta,nup-1-delta-gap,cuts,-1);
7081 SumTracks(seeds,arr);
7082 SignClusters(seeds,fnumber,fdensity);
7086 Info("Tracking()","\n\nSecondary seeding\t%d\n\n",seeds->GetEntriesFast());
7097 void AliTPCtrackerMI::SumTracks(TObjArray *arr1,TObjArray *&arr2)
7100 //sum tracks to common container
7101 //remove suspicious tracks
7102 // RS: Attention: supplied tracks come in the static array, don't delete them
7103 Int_t nseed = arr2->GetEntriesFast();
7104 for (Int_t i=0;i<nseed;i++){
7105 AliTPCseed *pt=(AliTPCseed*)arr2->UncheckedAt(i);
7108 // remove tracks with too big curvature
7110 if (TMath::Abs(pt->GetC())>AliTPCReconstructor::GetRecoParam()->GetMaxC()){
7111 MarkSeedFree( arr2->RemoveAt(i) );
7114 // REMOVE VERY SHORT TRACKS
7115 if (pt->GetNumberOfClusters()<20){
7116 MarkSeedFree( arr2->RemoveAt(i) );
7119 // NORMAL ACTIVE TRACK
7120 if (pt->IsActive()){
7121 arr1->AddLast(arr2->RemoveAt(i));
7124 //remove not usable tracks
7125 if (pt->GetRemoval()!=10){
7126 MarkSeedFree( arr2->RemoveAt(i) );
7130 // ENABLE ONLY ENOUGH GOOD STOPPED TRACKS
7131 if (pt->GetDensityFirst(20)>0.8 || pt->GetDensityFirst(30)>0.8 || pt->GetDensityFirst(40)>0.7)
7132 arr1->AddLast(arr2->RemoveAt(i));
7134 MarkSeedFree( arr2->RemoveAt(i) );
7138 // delete arr2; arr2 = 0; // RS: this is static array, don't delete it
7143 void AliTPCtrackerMI::ParallelTracking(TObjArray *const arr, Int_t rfirst, Int_t rlast)
7146 // try to track in parralel
7148 Int_t nseed=arr->GetEntriesFast();
7149 //prepare seeds for tracking
7150 for (Int_t i=0; i<nseed; i++) {
7151 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i), &t=*pt;
7153 if (!t.IsActive()) continue;
7154 // follow prolongation to the first layer
7155 if ( (fSectors ==fInnerSec) || (t.GetFirstPoint()-fkParam->GetNRowLow()>rfirst+1) )
7156 FollowProlongation(t, rfirst+1);
7161 for (Int_t nr=rfirst; nr>=rlast; nr--){
7162 if (nr<fInnerSec->GetNRows())
7163 fSectors = fInnerSec;
7165 fSectors = fOuterSec;
7166 // make indexes with the cluster tracks for given
7168 // find nearest cluster
7169 for (Int_t i=0; i<nseed; i++) {
7170 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i), &t=*pt;
7172 if (nr==80) pt->UpdateReference();
7173 if (!pt->IsActive()) continue;
7174 // if ( (fSectors ==fOuterSec) && (pt->fFirstPoint-fkParam->GetNRowLow())<nr) continue;
7175 if (pt->GetRelativeSector()>17) {
7178 UpdateClusters(t,nr);
7180 // prolonagate to the nearest cluster - if founded
7181 for (Int_t i=0; i<nseed; i++) {
7182 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
7184 if (!pt->IsActive()) continue;
7185 // if ((fSectors ==fOuterSec) && (pt->fFirstPoint-fkParam->GetNRowLow())<nr) continue;
7186 if (pt->GetRelativeSector()>17) {
7189 FollowToNextCluster(*pt,nr);
7194 void AliTPCtrackerMI::PrepareForBackProlongation(const TObjArray *const arr,Float_t fac) const
7198 // if we use TPC track itself we have to "update" covariance
7200 Int_t nseed= arr->GetEntriesFast();
7201 for (Int_t i=0;i<nseed;i++){
7202 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
7206 //rotate to current local system at first accepted point
7207 Int_t index = pt->GetClusterIndex2(pt->GetFirstPoint());
7208 Int_t sec = (index&0xff000000)>>24;
7210 Float_t angle1 = fInnerSec->GetAlpha()*sec+fInnerSec->GetAlphaShift();
7211 if (angle1>TMath::Pi())
7212 angle1-=2.*TMath::Pi();
7213 Float_t angle2 = pt->GetAlpha();
7215 if (TMath::Abs(angle1-angle2)>0.001){
7216 if (!pt->Rotate(angle1-angle2)) return;
7217 //angle2 = pt->GetAlpha();
7218 //pt->fRelativeSector = pt->GetAlpha()/fInnerSec->GetAlpha();
7219 //if (pt->GetAlpha()<0)
7220 // pt->fRelativeSector+=18;
7221 //sec = pt->fRelativeSector;
7230 void AliTPCtrackerMI::PrepareForProlongation(TObjArray *const arr, Float_t fac) const
7234 // if we use TPC track itself we have to "update" covariance
7236 Int_t nseed= arr->GetEntriesFast();
7237 for (Int_t i=0;i<nseed;i++){
7238 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
7241 pt->SetFirstPoint(pt->GetLastPoint());
7249 Int_t AliTPCtrackerMI::PropagateBack(const TObjArray *const arr)
7252 // make back propagation
7254 Int_t nseed= arr->GetEntriesFast();
7255 for (Int_t i=0;i<nseed;i++){
7256 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
7257 if (pt&& pt->GetKinkIndex(0)<=0) {
7258 //AliTPCseed *pt2 = new AliTPCseed(*pt);
7259 fSectors = fInnerSec;
7260 //FollowBackProlongation(*pt,fInnerSec->GetNRows()-1);
7261 //fSectors = fOuterSec;
7262 FollowBackProlongation(*pt,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1,1);
7263 //if (pt->GetNumberOfClusters()<(pt->fEsd->GetTPCclusters(0)) ){
7264 // Error("PropagateBack","Not prolonged track %d",pt->GetLabel());
7265 // FollowBackProlongation(*pt2,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1);
7268 if (pt&& pt->GetKinkIndex(0)>0) {
7269 AliESDkink * kink = fEvent->GetKink(pt->GetKinkIndex(0)-1);
7270 pt->SetFirstPoint(kink->GetTPCRow0());
7271 fSectors = fInnerSec;
7272 FollowBackProlongation(*pt,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1,1);
7280 Int_t AliTPCtrackerMI::PropagateForward2(const TObjArray *const arr)
7283 // make forward propagation
7285 Int_t nseed= arr->GetEntriesFast();
7287 for (Int_t i=0;i<nseed;i++){
7288 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
7290 FollowProlongation(*pt,0,1,1);
7299 Int_t AliTPCtrackerMI::PropagateForward()
7302 // propagate track forward
7304 Int_t nseed = fSeeds->GetEntriesFast();
7305 for (Int_t i=0;i<nseed;i++){
7306 AliTPCseed *pt = (AliTPCseed*)fSeeds->UncheckedAt(i);
7308 AliTPCseed &t = *pt;
7309 Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
7310 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
7311 if (alpha < 0. ) alpha += 2.*TMath::Pi();
7312 t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN);
7316 fSectors = fOuterSec;
7317 ParallelTracking(fSeeds,fOuterSec->GetNRows()+fInnerSec->GetNRows()-1,fInnerSec->GetNRows());
7318 fSectors = fInnerSec;
7319 ParallelTracking(fSeeds,fInnerSec->GetNRows()-1,0);
7328 Int_t AliTPCtrackerMI::PropagateBack(AliTPCseed *const pt, Int_t row0, Int_t row1)
7331 // make back propagation, in between row0 and row1
7335 fSectors = fInnerSec;
7338 if (row1<fSectors->GetNRows())
7341 r1 = fSectors->GetNRows()-1;
7343 if (row0<fSectors->GetNRows()&& r1>0 )
7344 FollowBackProlongation(*pt,r1);
7345 if (row1<=fSectors->GetNRows())
7348 r1 = row1 - fSectors->GetNRows();
7349 if (r1<=0) return 0;
7350 if (r1>=fOuterSec->GetNRows()) return 0;
7351 fSectors = fOuterSec;
7352 return FollowBackProlongation(*pt,r1);
7360 void AliTPCtrackerMI::GetShape(AliTPCseed * seed, Int_t row)
7362 // gets cluster shape
7364 AliTPCClusterParam * clparam = AliTPCcalibDB::Instance()->GetClusterParam();
7365 Float_t zdrift = TMath::Abs((fkParam->GetZLength(0)-TMath::Abs(seed->GetZ())));
7366 Int_t type = (seed->GetSector() < fkParam->GetNSector()/2) ? 0: (row>126) ? 1:2;
7367 Double_t angulary = seed->GetSnp();
7369 if (TMath::Abs(angulary)>AliTPCReconstructor::GetMaxSnpTracker()) {
7370 angulary = TMath::Sign(AliTPCReconstructor::GetMaxSnpTracker(),angulary);
7373 angulary = angulary*angulary/((1.-angulary)*(1.+angulary));
7374 Double_t angularz = seed->GetTgl()*seed->GetTgl()*(1.+angulary);
7376 Double_t sigmay = clparam->GetRMS0(0,type,zdrift,TMath::Sqrt(TMath::Abs(angulary)));
7377 Double_t sigmaz = clparam->GetRMS0(1,type,zdrift,TMath::Sqrt(TMath::Abs(angularz)));
7378 seed->SetCurrentSigmaY2(sigmay*sigmay);
7379 seed->SetCurrentSigmaZ2(sigmaz*sigmaz);
7380 // Float_t sd2 = TMath::Abs((fkParam->GetZLength(0)-TMath::Abs(seed->GetZ())))*fkParam->GetDiffL()*fkParam->GetDiffL();
7381 // // Float_t padlength = fkParam->GetPadPitchLength(seed->fSector);
7382 // Float_t padlength = GetPadPitchLength(row);
7384 // Float_t sresy = (seed->GetSector() < fkParam->GetNSector()/2) ? 0.2 :0.3;
7385 // seed->SetCurrentSigmaY2(sd2+padlength*padlength*angulary/12.+sresy*sresy);
7387 // Float_t sresz = fkParam->GetZSigma();
7388 // seed->SetCurrentSigmaZ2(sd2+padlength*padlength*angularz*angularz*(1+angulary)/12.+sresz*sresz);
7390 Float_t wy = GetSigmaY(seed);
7391 Float_t wz = GetSigmaZ(seed);
7394 if (TMath::Abs(wy/seed->fCurrentSigmaY2-1)>0.0001 || TMath::Abs(wz/seed->fCurrentSigmaZ2-1)>0.0001 ){
7395 printf("problem\n");
7402 //__________________________________________________________________________
7403 void AliTPCtrackerMI::CookLabel(AliKalmanTrack *tk, Float_t wrong) const {
7404 //--------------------------------------------------------------------
7405 //This function "cooks" a track label. If label<0, this track is fake.
7406 //--------------------------------------------------------------------
7407 AliTPCseed * t = dynamic_cast<AliTPCseed*>(tk);
7409 printf("%s:%d wrong type \n",(char*)__FILE__,__LINE__);
7413 Int_t noc=t->GetNumberOfClusters();
7415 //printf("\nnot founded prolongation\n\n\n");
7421 AliTPCclusterMI *clusters[160];
7423 for (Int_t i=0;i<160;i++) {
7430 for (i=0; i<160 && current<noc; i++) {
7432 Int_t index=t->GetClusterIndex2(i);
7433 if (index<=0) continue;
7434 if (index&0x8000) continue;
7436 //clusters[current]=GetClusterMI(index);
7437 if (t->GetClusterPointer(i)){
7438 clusters[current]=t->GetClusterPointer(i);
7444 Int_t lab=123456789;
7445 for (i=0; i<noc; i++) {
7446 AliTPCclusterMI *c=clusters[i];
7448 lab=TMath::Abs(c->GetLabel(0));
7450 for (j=0; j<noc; j++) if (lb[j]==lab || mx[j]==0) break;
7456 for (i=0; i<noc; i++) if (mx[i]>max) {max=mx[i]; lab=lb[i];}
7458 for (i=0; i<noc; i++) {
7459 AliTPCclusterMI *c=clusters[i];
7461 if (TMath::Abs(c->GetLabel(1)) == lab ||
7462 TMath::Abs(c->GetLabel(2)) == lab ) max++;
7464 if (noc<=0) { lab=-1; return;}
7465 if ((1.- Float_t(max)/(noc)) > wrong) lab=-lab;
7468 Int_t tail=Int_t(0.10*noc);
7471 for (i=1; i<160&&ind<tail; i++) {
7472 // AliTPCclusterMI *c=clusters[noc-i];
7473 AliTPCclusterMI *c=clusters[i];
7475 if (lab == TMath::Abs(c->GetLabel(0)) ||
7476 lab == TMath::Abs(c->GetLabel(1)) ||
7477 lab == TMath::Abs(c->GetLabel(2))) max++;
7480 if (max < Int_t(0.5*tail)) lab=-lab;
7487 //delete[] clusters;
7491 //__________________________________________________________________________
7492 Int_t AliTPCtrackerMI::CookLabel(AliTPCseed *const t, Float_t wrong,Int_t first, Int_t last) const {
7493 //--------------------------------------------------------------------
7494 //This function "cooks" a track label. If label<0, this track is fake.
7495 //--------------------------------------------------------------------
7496 Int_t noc=t->GetNumberOfClusters();
7498 //printf("\nnot founded prolongation\n\n\n");
7504 AliTPCclusterMI *clusters[160];
7506 for (Int_t i=0;i<160;i++) {
7513 for (i=0; i<160 && current<noc; i++) {
7514 if (i<first) continue;
7515 if (i>last) continue;
7516 Int_t index=t->GetClusterIndex2(i);
7517 if (index<=0) continue;
7518 if (index&0x8000) continue;
7520 //clusters[current]=GetClusterMI(index);
7521 if (t->GetClusterPointer(i)){
7522 clusters[current]=t->GetClusterPointer(i);
7527 //if (noc<5) return -1;
7528 Int_t lab=123456789;
7529 for (i=0; i<noc; i++) {
7530 AliTPCclusterMI *c=clusters[i];
7532 lab=TMath::Abs(c->GetLabel(0));
7534 for (j=0; j<noc; j++) if (lb[j]==lab || mx[j]==0) break;
7540 for (i=0; i<noc; i++) if (mx[i]>max) {max=mx[i]; lab=lb[i];}
7542 for (i=0; i<noc; i++) {
7543 AliTPCclusterMI *c=clusters[i];
7545 if (TMath::Abs(c->GetLabel(1)) == lab ||
7546 TMath::Abs(c->GetLabel(2)) == lab ) max++;
7548 if (noc<=0) { lab=-1; return -1;}
7549 if ((1.- Float_t(max)/(noc)) > wrong) lab=-lab;
7552 Int_t tail=Int_t(0.10*noc);
7555 for (i=1; i<160&&ind<tail; i++) {
7556 // AliTPCclusterMI *c=clusters[noc-i];
7557 AliTPCclusterMI *c=clusters[i];
7559 if (lab == TMath::Abs(c->GetLabel(0)) ||
7560 lab == TMath::Abs(c->GetLabel(1)) ||
7561 lab == TMath::Abs(c->GetLabel(2))) max++;
7564 if (max < Int_t(0.5*tail)) lab=-lab;
7567 // t->SetLabel(lab);
7571 //delete[] clusters;
7575 Int_t AliTPCtrackerMI::GetRowNumber(Double_t x[3]) const
7577 //return pad row number for given x vector
7578 Float_t phi = TMath::ATan2(x[1],x[0]);
7579 if(phi<0) phi=2.*TMath::Pi()+phi;
7580 // Get the local angle in the sector philoc
7581 const Float_t kRaddeg = 180/3.14159265358979312;
7582 Float_t phiangle = (Int_t (phi*kRaddeg/20.) + 0.5)*20./kRaddeg;
7583 Double_t localx = x[0]*TMath::Cos(phiangle)-x[1]*TMath::Sin(phiangle);
7584 return GetRowNumber(localx);
7589 void AliTPCtrackerMI::MakeESDBitmaps(AliTPCseed *t, AliESDtrack *esd)
7591 //-----------------------------------------------------------------------
7592 // Fill the cluster and sharing bitmaps of the track
7593 //-----------------------------------------------------------------------
7595 Int_t firstpoint = 0;
7596 Int_t lastpoint = 159;
7597 AliTPCTrackerPoint *point;
7598 AliTPCclusterMI *cluster;
7601 TBits clusterMap(159);
7602 TBits sharedMap(159);
7604 for (int iter=firstpoint; iter<lastpoint; iter++) {
7605 // Change to cluster pointers to see if we have a cluster at given padrow
7607 cluster = t->GetClusterPointer(iter);
7609 clusterMap.SetBitNumber(iter, kTRUE);
7610 point = t->GetTrackPoint(iter);
7611 if (point->IsShared())
7612 sharedMap.SetBitNumber(iter,kTRUE);
7614 if (t->GetClusterIndex(iter) >= 0 && (t->GetClusterIndex(iter) & 0x8000) == 0) {
7615 fitMap.SetBitNumber(iter, kTRUE);
7619 esd->SetTPCClusterMap(clusterMap);
7620 esd->SetTPCSharedMap(sharedMap);
7621 esd->SetTPCFitMap(fitMap);
7622 if (nclsf != t->GetNumberOfClusters())
7623 AliWarning(Form("Inconsistency between ncls %d and indices %d (found %d)",t->GetNumberOfClusters(),nclsf,esd->GetTPCClusterMap().CountBits()));
7626 Bool_t AliTPCtrackerMI::IsFindable(AliTPCseed & track){
7628 // return flag if there is findable cluster at given position
7631 Float_t z = track.GetZ();
7633 if (TMath::Abs(z)<(AliTPCReconstructor::GetCtgRange()*track.GetX()+kDeltaZ) &&
7634 TMath::Abs(z)<fkParam->GetZLength(0) &&
7635 (TMath::Abs(track.GetSnp())<AliTPCReconstructor::GetMaxSnpTracker()))
7641 void AliTPCtrackerMI::AddCovariance(AliTPCseed * seed){
7643 // Adding systematic error estimate to the covariance matrix
7644 // !!!! the systematic error for element 4 is in 1/cm not in pt
7646 const Double_t *param = AliTPCReconstructor::GetRecoParam()->GetSystematicError();
7648 // use only the diagonal part if not specified otherwise
7649 if (!AliTPCReconstructor::GetRecoParam()->GetUseSystematicCorrelation()) return AddCovarianceAdd(seed);
7651 Double_t *covarS= (Double_t*)seed->GetCovariance();
7652 Double_t factor[5]={1,1,1,1,1};
7653 Double_t facC = AliTracker::GetBz()*kB2C;
7654 factor[0]= TMath::Sqrt(TMath::Abs((covarS[0] + param[0]*param[0])/covarS[0]));
7655 factor[1]= TMath::Sqrt(TMath::Abs((covarS[2] + param[1]*param[1])/covarS[2]));
7656 factor[2]= TMath::Sqrt(TMath::Abs((covarS[5] + param[2]*param[2])/covarS[5]));
7657 factor[3]= TMath::Sqrt(TMath::Abs((covarS[9] + param[3]*param[3])/covarS[9]));
7658 factor[4]= TMath::Sqrt(TMath::Abs((covarS[14] + facC*facC*param[4]*param[4])/covarS[14]));
7664 for (Int_t i=0; i<5; i++){
7665 for (Int_t j=i; j<5; j++){
7666 Int_t index=seed->GetIndex(i,j);
7667 covarS[index]*=factor[i]*factor[j];
7673 void AliTPCtrackerMI::AddCovarianceAdd(AliTPCseed * seed){
7675 // Adding systematic error - as additive factor without correlation
7677 const Double_t *param = AliTPCReconstructor::GetRecoParam()->GetSystematicError();
7678 Double_t *covarIn= (Double_t*)seed->GetCovariance();
7680 for (Int_t i=0;i<15;i++) covar[i]=0;
7686 covar[0] = param[0]*param[0];
7687 covar[2] = param[1]*param[1];
7688 covar[5] = param[2]*param[2];
7689 covar[9] = param[3]*param[3];
7690 Double_t facC = AliTracker::GetBz()*kB2C;
7691 covar[14]= param[4]*param[4]*facC*facC;
7693 covar[1]=TMath::Sqrt((covar[0]*covar[2]))*covarIn[1]/TMath::Sqrt((covarIn[0]*covarIn[2]));
7695 covar[3]=TMath::Sqrt((covar[0]*covar[5]))*covarIn[3]/TMath::Sqrt((covarIn[0]*covarIn[5]));
7696 covar[4]=TMath::Sqrt((covar[2]*covar[5]))*covarIn[4]/TMath::Sqrt((covarIn[2]*covarIn[5]));
7698 covar[6]=TMath::Sqrt((covar[0]*covar[9]))*covarIn[6]/TMath::Sqrt((covarIn[0]*covarIn[9]));
7699 covar[7]=TMath::Sqrt((covar[2]*covar[9]))*covarIn[7]/TMath::Sqrt((covarIn[2]*covarIn[9]));
7700 covar[8]=TMath::Sqrt((covar[5]*covar[9]))*covarIn[8]/TMath::Sqrt((covarIn[5]*covarIn[9]));
7702 covar[10]=TMath::Sqrt((covar[0]*covar[14]))*covarIn[10]/TMath::Sqrt((covarIn[0]*covarIn[14]));
7703 covar[11]=TMath::Sqrt((covar[2]*covar[14]))*covarIn[11]/TMath::Sqrt((covarIn[2]*covarIn[14]));
7704 covar[12]=TMath::Sqrt((covar[5]*covar[14]))*covarIn[12]/TMath::Sqrt((covarIn[5]*covarIn[14]));
7705 covar[13]=TMath::Sqrt((covar[9]*covar[14]))*covarIn[13]/TMath::Sqrt((covarIn[9]*covarIn[14]));
7707 seed->AddCovariance(covar);
7710 //________________________________________
7711 void AliTPCtrackerMI::MarkSeedFree(TObject *sd)
7713 // account that this seed is "deleted"
7714 AliTPCseed* seed = dynamic_cast<AliTPCseed*>(sd);
7715 if (!sd) {AliError(Form("Freeing of non-AliTPCseed %p from the pool is requested",sd)); return;}
7716 int id = seed->GetPoolID();
7717 if (id<0) {AliError(Form("Freeing of seed %p NOT from the pool is requested",sd)); return;}
7718 // AliInfo(Form("%d %p",id, seed));
7719 fSeedsPool->RemoveAt(id);
7720 if (fFreeSeedsID.GetSize()<=fNFreeSeeds) fFreeSeedsID.Set( 2*fNFreeSeeds + 100 );
7721 fFreeSeedsID.GetArray()[fNFreeSeeds++] = id;
7724 //________________________________________
7725 TObject *&AliTPCtrackerMI::NextFreeSeed()
7727 // return next free slot where the seed can be created
7728 fLastSeedID = fNFreeSeeds ? fFreeSeedsID.GetArray()[--fNFreeSeeds] : fSeedsPool->GetEntriesFast();
7729 // AliInfo(Form("%d",fLastSeedID));
7730 return (*fSeedsPool)[ fLastSeedID ];
7734 //________________________________________
7735 void AliTPCtrackerMI::ResetSeedsPool()
7737 // mark all seeds in the pool as unused
7738 AliInfo(Form("CurrentSize: %d, BookedUpTo: %d, free: %d",fSeedsPool->GetSize(),fSeedsPool->GetEntriesFast(),fNFreeSeeds));
7740 fSeedsPool->Clear("C"); // RS: nominally the seeds may allocate memory...