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/GeV)
83 // The default values are 0.
85 // The sytematic errors are added to the covariance matrix in following places:
87 // 1. During fisrt itteration - AliTPCtrackerMI::FillESD
88 // 2. Second iteration -
89 // 2.a ITS->TPC - AliTPCtrackerMI::ReadSeeds
90 // 2.b TPC->TRD - AliTPCtrackerMI::PropagateBack
91 // 3. Third iteration -
92 // 3.a TRD->TPC - AliTPCtrackerMI::ReadSeeds
93 // 3.b TPC->ITS - AliTPCtrackerMI::RefitInward
95 // There are several places in the code which can be numerically debuged
96 // This code is keeped in order to enable code development and to check the calibration implementtion
98 // 1. ErrParam stream - dump information about
100 // 2.a) cluster error estimate
101 // 3.a) cluster shape estimate
104 // Debug streamer levels:
106 //-------------------------------------------------------
111 #include "Riostream.h"
112 #include <TClonesArray.h>
114 #include <TObjArray.h>
116 #include <TGraphErrors.h>
118 #include "AliComplexCluster.h"
119 #include "AliESDEvent.h"
120 #include "AliESDtrack.h"
121 #include "AliESDVertex.h"
124 #include "AliHelix.h"
125 #include "AliRunLoader.h"
126 #include "AliTPCClustersRow.h"
127 #include "AliTPCParam.h"
128 #include "AliTPCReconstructor.h"
129 #include "AliTPCpolyTrack.h"
130 #include "AliTPCreco.h"
131 #include "AliTPCseed.h"
133 #include "AliTPCtrackerSector.h"
134 #include "AliTPCtrackerMI.h"
135 #include "TStopwatch.h"
136 #include "AliTPCReconstructor.h"
137 #include "AliAlignObj.h"
138 #include "AliTrackPointArray.h"
140 #include "AliTPCcalibDB.h"
141 #include "AliTPCcalibDButil.h"
142 #include "AliTPCTransform.h"
143 #include "AliTPCClusterParam.h"
144 #include "AliTPCdEdxInfo.h"
145 #include "AliDCSSensorArray.h"
146 #include "AliDCSSensor.h"
148 #include "AliCosmicTracker.h"
152 ClassImp(AliTPCtrackerMI)
156 class AliTPCFastMath {
159 static Double_t FastAsin(Double_t x);
161 static Double_t fgFastAsin[20000]; //lookup table for fast asin computation
164 Double_t AliTPCFastMath::fgFastAsin[20000];
165 AliTPCFastMath gAliTPCFastMath; // needed to fill the LUT
167 AliTPCFastMath::AliTPCFastMath(){
169 // initialized lookup table;
170 for (Int_t i=0;i<10000;i++){
171 fgFastAsin[2*i] = TMath::ASin(i/10000.);
172 fgFastAsin[2*i+1] = (TMath::ASin((i+1)/10000.)-fgFastAsin[2*i]);
176 Double_t AliTPCFastMath::FastAsin(Double_t x){
178 // return asin using lookup table
180 Int_t index = int(x*10000);
181 return fgFastAsin[2*index]+(x*10000.-index)*fgFastAsin[2*index+1];
184 Int_t index = int(x*10000);
185 return -(fgFastAsin[2*index]+(x*10000.-index)*fgFastAsin[2*index+1]);
187 //__________________________________________________________________
188 AliTPCtrackerMI::AliTPCtrackerMI()
215 // default constructor
217 for (Int_t irow=0; irow<200; irow++){
224 //_____________________________________________________________________
228 Int_t AliTPCtrackerMI::UpdateTrack(AliTPCseed * track, Int_t accept){
230 //update track information using current cluster - track->fCurrentCluster
233 AliTPCclusterMI* c =track->GetCurrentCluster();
234 if (accept > 0) //sign not accepted clusters
235 track->SetCurrentClusterIndex1(track->GetCurrentClusterIndex1() | 0x8000);
236 else // unsign accpeted clusters
237 track->SetCurrentClusterIndex1(track->GetCurrentClusterIndex1() & 0xffff7fff);
238 UInt_t i = track->GetCurrentClusterIndex1();
240 Int_t sec=(i&0xff000000)>>24;
241 //Int_t row = (i&0x00ff0000)>>16;
242 track->SetRow((i&0x00ff0000)>>16);
243 track->SetSector(sec);
244 // Int_t index = i&0xFFFF;
245 if (sec>=fkParam->GetNInnerSector()) track->SetRow(track->GetRow()+fkParam->GetNRowLow());
246 track->SetClusterIndex2(track->GetRow(), i);
247 //track->fFirstPoint = row;
248 //if ( track->fLastPoint<row) track->fLastPoint =row;
249 // if (track->fRow<0 || track->fRow>160) {
250 // printf("problem\n");
252 if (track->GetFirstPoint()>track->GetRow())
253 track->SetFirstPoint(track->GetRow());
254 if (track->GetLastPoint()<track->GetRow())
255 track->SetLastPoint(track->GetRow());
258 track->SetClusterPointer(track->GetRow(),c);
261 Double_t angle2 = track->GetSnp()*track->GetSnp();
263 //SET NEW Track Point
265 if (angle2<1) //PH sometimes angle2 is very big. To be investigated...
267 angle2 = TMath::Sqrt(angle2/(1-angle2));
268 AliTPCTrackerPoint &point =*(track->GetTrackPoint(track->GetRow()));
270 point.SetSigmaY(c->GetSigmaY2()/track->GetCurrentSigmaY2());
271 point.SetSigmaZ(c->GetSigmaZ2()/track->GetCurrentSigmaZ2());
272 point.SetErrY(sqrt(track->GetErrorY2()));
273 point.SetErrZ(sqrt(track->GetErrorZ2()));
275 point.SetX(track->GetX());
276 point.SetY(track->GetY());
277 point.SetZ(track->GetZ());
278 point.SetAngleY(angle2);
279 point.SetAngleZ(track->GetTgl());
280 if (point.IsShared()){
281 track->SetErrorY2(track->GetErrorY2()*4);
282 track->SetErrorZ2(track->GetErrorZ2()*4);
286 Double_t chi2 = track->GetPredictedChi2(track->GetCurrentCluster());
288 // track->SetErrorY2(track->GetErrorY2()*1.3);
289 // track->SetErrorY2(track->GetErrorY2()+0.01);
290 // track->SetErrorZ2(track->GetErrorZ2()*1.3);
291 // track->SetErrorZ2(track->GetErrorZ2()+0.005);
293 if (accept>0) return 0;
294 if (track->GetNumberOfClusters()%20==0){
295 // if (track->fHelixIn){
296 // TClonesArray & larr = *(track->fHelixIn);
297 // Int_t ihelix = larr.GetEntriesFast();
298 // new(larr[ihelix]) AliHelix(*track) ;
301 track->SetNoCluster(0);
302 return track->Update(c,chi2,i);
307 Int_t AliTPCtrackerMI::AcceptCluster(AliTPCseed * seed, AliTPCclusterMI * cluster)
310 // decide according desired precision to accept given
311 // cluster for tracking
313 seed->GetProlongation(cluster->GetX(),yt,zt);
314 Double_t sy2=ErrY2(seed,cluster);
315 Double_t sz2=ErrZ2(seed,cluster);
317 Double_t sdistancey2 = sy2+seed->GetSigmaY2();
318 Double_t sdistancez2 = sz2+seed->GetSigmaZ2();
319 Double_t dy=seed->GetCurrentCluster()->GetY()-yt;
320 Double_t dz=seed->GetCurrentCluster()->GetZ()-zt;
321 Double_t rdistancey2 = (seed->GetCurrentCluster()->GetY()-yt)*
322 (seed->GetCurrentCluster()->GetY()-yt)/sdistancey2;
323 Double_t rdistancez2 = (seed->GetCurrentCluster()->GetZ()-zt)*
324 (seed->GetCurrentCluster()->GetZ()-zt)/sdistancez2;
326 Double_t rdistance2 = rdistancey2+rdistancez2;
329 if (AliTPCReconstructor::StreamLevel()>2 && seed->GetNumberOfClusters()>20) {
330 Float_t rmsy2 = seed->GetCurrentSigmaY2();
331 Float_t rmsz2 = seed->GetCurrentSigmaZ2();
332 Float_t rmsy2p30 = seed->GetCMeanSigmaY2p30();
333 Float_t rmsz2p30 = seed->GetCMeanSigmaZ2p30();
334 Float_t rmsy2p30R = seed->GetCMeanSigmaY2p30R();
335 Float_t rmsz2p30R = seed->GetCMeanSigmaZ2p30R();
336 AliExternalTrackParam param(*seed);
337 static TVectorD gcl(3),gtr(3);
339 param.GetXYZ(gcl.GetMatrixArray());
340 cluster->GetGlobalXYZ(gclf);
341 gcl[0]=gclf[0]; gcl[1]=gclf[1]; gcl[2]=gclf[2];
344 if (AliTPCReconstructor::StreamLevel()>2) {
345 (*fDebugStreamer)<<"ErrParam"<<
358 "rmsy2p30="<<rmsy2p30<<
359 "rmsz2p30="<<rmsz2p30<<
360 "rmsy2p30R="<<rmsy2p30R<<
361 "rmsz2p30R="<<rmsz2p30R<<
362 // normalize distance -
363 "rdisty="<<rdistancey2<<
364 "rdistz="<<rdistancez2<<
365 "rdist="<<rdistance2<< //
369 //return 0; // temporary
370 if (rdistance2>32) return 3;
373 if ((rdistancey2>9. || rdistancez2>9.) && cluster->GetType()==0)
374 return 2; //suspisiouce - will be changed
376 if ((rdistancey2>6.25 || rdistancez2>6.25) && cluster->GetType()>0)
377 // strict cut on overlaped cluster
378 return 2; //suspisiouce - will be changed
380 if ( (rdistancey2>1. || rdistancez2>6.25 )
381 && cluster->GetType()<0){
382 seed->SetNFoundable(seed->GetNFoundable()-1);
386 if (fUseHLTClusters == 3 || fUseHLTClusters == 4) {
388 if(!AliTPCReconstructor::GetRecoParam()->GetUseHLTOnePadCluster()) {
389 if (TMath::Abs(cluster->GetSigmaY2()) < kAlmost0)
402 //_____________________________________________________________________________
403 AliTPCtrackerMI::AliTPCtrackerMI(const AliTPCParam *par):
405 fkNIS(par->GetNInnerSector()/2),
407 fkNOS(par->GetNOuterSector()/2),
429 //---------------------------------------------------------------------
430 // The main TPC tracker constructor
431 //---------------------------------------------------------------------
432 fInnerSec=new AliTPCtrackerSector[fkNIS];
433 fOuterSec=new AliTPCtrackerSector[fkNOS];
436 for (i=0; i<fkNIS; i++) fInnerSec[i].Setup(par,0);
437 for (i=0; i<fkNOS; i++) fOuterSec[i].Setup(par,1);
440 Int_t nrowlow = par->GetNRowLow();
441 Int_t nrowup = par->GetNRowUp();
444 for (i=0;i<nrowlow;i++){
445 fXRow[i] = par->GetPadRowRadiiLow(i);
446 fPadLength[i]= par->GetPadPitchLength(0,i);
447 fYMax[i] = fXRow[i]*TMath::Tan(0.5*par->GetInnerAngle());
451 for (i=0;i<nrowup;i++){
452 fXRow[i+nrowlow] = par->GetPadRowRadiiUp(i);
453 fPadLength[i+nrowlow] = par->GetPadPitchLength(60,i);
454 fYMax[i+nrowlow] = fXRow[i+nrowlow]*TMath::Tan(0.5*par->GetOuterAngle());
457 if (AliTPCReconstructor::StreamLevel()>0) {
458 fDebugStreamer = new TTreeSRedirector("TPCdebug.root");
461 fSeedsPool = new TClonesArray("AliTPCseed",1000);
463 //________________________________________________________________________
464 AliTPCtrackerMI::AliTPCtrackerMI(const AliTPCtrackerMI &t):
490 //------------------------------------
491 // dummy copy constructor
492 //------------------------------------------------------------------
494 for (Int_t irow=0; irow<200; irow++){
501 AliTPCtrackerMI & AliTPCtrackerMI::operator=(const AliTPCtrackerMI& /*r*/)
503 //------------------------------
505 //--------------------------------------------------------------
508 //_____________________________________________________________________________
509 AliTPCtrackerMI::~AliTPCtrackerMI() {
510 //------------------------------------------------------------------
511 // TPC tracker destructor
512 //------------------------------------------------------------------
519 if (fDebugStreamer) delete fDebugStreamer;
520 if (fSeedsPool) delete fSeedsPool;
524 void AliTPCtrackerMI::FillESD(const TObjArray* arr)
528 //fill esds using updated tracks
531 // write tracks to the event
532 // store index of the track
533 Int_t nseed=arr->GetEntriesFast();
534 //FindKinks(arr,fEvent);
535 for (Int_t i=0; i<nseed; i++) {
536 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
540 if (AliTPCReconstructor::StreamLevel()>1) {
541 (*fDebugStreamer)<<"Track0"<<
545 // pt->PropagateTo(fkParam->GetInnerRadiusLow());
546 if (pt->GetKinkIndex(0)<=0){ //don't propagate daughter tracks
547 pt->PropagateTo(fkParam->GetInnerRadiusLow());
550 if (( pt->GetPoints()[2]- pt->GetPoints()[0])>5 && pt->GetPoints()[3]>0.8){
552 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
553 iotrack.SetTPCPoints(pt->GetPoints());
554 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
555 iotrack.SetV0Indexes(pt->GetV0Indexes());
556 // iotrack.SetTPCpid(pt->fTPCr);
557 //iotrack.SetTPCindex(i);
558 fEvent->AddTrack(&iotrack);
562 if ( (pt->GetNumberOfClusters()>70)&& (Float_t(pt->GetNumberOfClusters())/Float_t(pt->GetNFoundable()))>0.55) {
564 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
565 iotrack.SetTPCPoints(pt->GetPoints());
566 //iotrack.SetTPCindex(i);
567 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
568 iotrack.SetV0Indexes(pt->GetV0Indexes());
569 // iotrack.SetTPCpid(pt->fTPCr);
570 fEvent->AddTrack(&iotrack);
574 // short tracks - maybe decays
576 if ( (pt->GetNumberOfClusters()>30) && (Float_t(pt->GetNumberOfClusters())/Float_t(pt->GetNFoundable()))>0.70) {
577 Int_t found,foundable,shared;
578 pt->GetClusterStatistic(0,60,found, foundable,shared,kFALSE);
579 if ( (found>20) && (pt->GetNShared()/float(pt->GetNumberOfClusters())<0.2)){
581 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
582 //iotrack.SetTPCindex(i);
583 iotrack.SetTPCPoints(pt->GetPoints());
584 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
585 iotrack.SetV0Indexes(pt->GetV0Indexes());
586 //iotrack.SetTPCpid(pt->fTPCr);
587 fEvent->AddTrack(&iotrack);
592 if ( (pt->GetNumberOfClusters()>20) && (Float_t(pt->GetNumberOfClusters())/Float_t(pt->GetNFoundable()))>0.8) {
593 Int_t found,foundable,shared;
594 pt->GetClusterStatistic(0,60,found, foundable,shared,kFALSE);
595 if (found<20) continue;
596 if (pt->GetNShared()/float(pt->GetNumberOfClusters())>0.2) continue;
599 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
600 iotrack.SetTPCPoints(pt->GetPoints());
601 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
602 iotrack.SetV0Indexes(pt->GetV0Indexes());
603 //iotrack.SetTPCpid(pt->fTPCr);
604 //iotrack.SetTPCindex(i);
605 fEvent->AddTrack(&iotrack);
608 // short tracks - secondaties
610 if ( (pt->GetNumberOfClusters()>30) ) {
611 Int_t found,foundable,shared;
612 pt->GetClusterStatistic(128,158,found, foundable,shared,kFALSE);
613 if ( (found>20) && (pt->GetNShared()/float(pt->GetNumberOfClusters())<0.2) &&float(found)/float(foundable)>0.8){
615 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
616 iotrack.SetTPCPoints(pt->GetPoints());
617 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
618 iotrack.SetV0Indexes(pt->GetV0Indexes());
619 //iotrack.SetTPCpid(pt->fTPCr);
620 //iotrack.SetTPCindex(i);
621 fEvent->AddTrack(&iotrack);
626 if ( (pt->GetNumberOfClusters()>15)) {
627 Int_t found,foundable,shared;
628 pt->GetClusterStatistic(138,158,found, foundable,shared,kFALSE);
629 if (found<15) continue;
630 if (foundable<=0) continue;
631 if (pt->GetNShared()/float(pt->GetNumberOfClusters())>0.2) continue;
632 if (float(found)/float(foundable)<0.8) continue;
635 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
636 iotrack.SetTPCPoints(pt->GetPoints());
637 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
638 iotrack.SetV0Indexes(pt->GetV0Indexes());
639 // iotrack.SetTPCpid(pt->fTPCr);
640 //iotrack.SetTPCindex(i);
641 fEvent->AddTrack(&iotrack);
645 // >> account for suppressed tracks in the kink indices (RS)
646 int nESDtracks = fEvent->GetNumberOfTracks();
647 for (int it=nESDtracks;it--;) {
648 AliESDtrack* esdTr = fEvent->GetTrack(it);
649 if (!esdTr || !esdTr->GetKinkIndex(0)) continue;
650 for (int ik=0;ik<3;ik++) {
652 if (!(knkId=esdTr->GetKinkIndex(ik))) break; // no more kinks for this track
653 AliESDkink* kink = fEvent->GetKink(TMath::Abs(knkId)-1);
655 AliError(Form("ESDTrack%d refers to non-existing kink %d",it,TMath::Abs(knkId)-1));
658 kink->SetIndex(it, knkId<0 ? 0:1); // update track index of the kink: mother at 0, daughter at 1
661 // << account for suppressed tracks in the kink indices (RS)
662 AliInfo(Form("Number of filled ESDs-\t%d\n",fEvent->GetNumberOfTracks()));
670 Double_t AliTPCtrackerMI::ErrY2(AliTPCseed* seed, const AliTPCclusterMI * cl){
673 // Use calibrated cluster error from OCDB
675 AliTPCClusterParam * clparam = AliTPCcalibDB::Instance()->GetClusterParam();
677 Float_t z = TMath::Abs(fkParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
678 Int_t ctype = cl->GetType();
679 Int_t type = (cl->GetRow()<63) ? 0: (cl->GetRow()>126) ? 1:2;
680 Double_t angle = seed->GetSnp()*seed->GetSnp();
681 angle = TMath::Sqrt(TMath::Abs(angle/(1.-angle)));
682 Double_t erry2 = clparam->GetError0Par(0,type, z,angle);
684 erry2+=0.5; // edge cluster
687 seed->SetErrorY2(erry2);
691 //calculate look-up table at the beginning
692 // static Bool_t ginit = kFALSE;
693 // static Float_t gnoise1,gnoise2,gnoise3;
694 // static Float_t ggg1[10000];
695 // static Float_t ggg2[10000];
696 // static Float_t ggg3[10000];
697 // static Float_t glandau1[10000];
698 // static Float_t glandau2[10000];
699 // static Float_t glandau3[10000];
701 // static Float_t gcor01[500];
702 // static Float_t gcor02[500];
703 // static Float_t gcorp[500];
707 // if (ginit==kFALSE){
708 // for (Int_t i=1;i<500;i++){
709 // Float_t rsigma = float(i)/100.;
710 // gcor02[i] = TMath::Max(0.78 +TMath::Exp(7.4*(rsigma-1.2)),0.6);
711 // gcor01[i] = TMath::Max(0.72 +TMath::Exp(3.36*(rsigma-1.2)),0.6);
712 // gcorp[i] = TMath::Max(TMath::Power((rsigma+0.5),1.5),1.2);
716 // for (Int_t i=3;i<10000;i++){
720 // Float_t amp = float(i);
721 // Float_t padlength =0.75;
722 // gnoise1 = 0.0004/padlength;
723 // Float_t nel = 0.268*amp;
724 // Float_t nprim = 0.155*amp;
725 // ggg1[i] = fkParam->GetDiffT()*fkParam->GetDiffT()*(2+0.001*nel/(padlength*padlength))/nel;
726 // glandau1[i] = (2.+0.12*nprim)*0.5* (2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
727 // if (glandau1[i]>1) glandau1[i]=1;
728 // glandau1[i]*=padlength*padlength/12.;
732 // gnoise2 = 0.0004/padlength;
734 // nprim = 0.133*amp;
735 // ggg2[i] = fkParam->GetDiffT()*fkParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
736 // glandau2[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
737 // if (glandau2[i]>1) glandau2[i]=1;
738 // glandau2[i]*=padlength*padlength/12.;
743 // gnoise3 = 0.0004/padlength;
745 // nprim = 0.133*amp;
746 // ggg3[i] = fkParam->GetDiffT()*fkParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
747 // glandau3[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
748 // if (glandau3[i]>1) glandau3[i]=1;
749 // glandau3[i]*=padlength*padlength/12.;
757 // Int_t amp = int(TMath::Abs(cl->GetQ()));
759 // seed->SetErrorY2(1.);
763 // Float_t z = TMath::Abs(fkParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
764 // Int_t ctype = cl->GetType();
765 // Float_t padlength= GetPadPitchLength(seed->GetRow());
766 // Double_t angle2 = seed->GetSnp()*seed->GetSnp();
767 // angle2 = angle2/(1-angle2);
769 // //cluster "quality"
770 // Int_t rsigmay = int(100.*cl->GetSigmaY2()/(seed->GetCurrentSigmaY2()));
773 // if (fSectors==fInnerSec){
774 // snoise2 = gnoise1;
775 // res = ggg1[amp]*z+glandau1[amp]*angle2;
776 // if (ctype==0) res *= gcor01[rsigmay];
779 // res*= gcorp[rsigmay];
783 // if (padlength<1.1){
784 // snoise2 = gnoise2;
785 // res = ggg2[amp]*z+glandau2[amp]*angle2;
786 // if (ctype==0) res *= gcor02[rsigmay];
789 // res*= gcorp[rsigmay];
793 // snoise2 = gnoise3;
794 // res = ggg3[amp]*z+glandau3[amp]*angle2;
795 // if (ctype==0) res *= gcor02[rsigmay];
798 // res*= gcorp[rsigmay];
805 // res*=2.4; // overestimate error 2 times
809 // if (res<2*snoise2)
812 // seed->SetErrorY2(res);
820 Double_t AliTPCtrackerMI::ErrZ2(AliTPCseed* seed, const AliTPCclusterMI * cl){
823 // Use calibrated cluster error from OCDB
825 AliTPCClusterParam * clparam = AliTPCcalibDB::Instance()->GetClusterParam();
827 Float_t z = TMath::Abs(fkParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
828 Int_t ctype = cl->GetType();
829 Int_t type = (cl->GetRow()<63) ? 0: (cl->GetRow()>126) ? 1:2;
831 Double_t angle2 = seed->GetSnp()*seed->GetSnp();
832 angle2 = seed->GetTgl()*seed->GetTgl()*(1+angle2/(1-angle2));
833 Double_t angle = TMath::Sqrt(TMath::Abs(angle2));
834 Double_t errz2 = clparam->GetError0Par(1,type, z,angle);
836 errz2+=0.5; // edge cluster
839 seed->SetErrorZ2(errz2);
845 // //seed->SetErrorY2(0.1);
847 // //calculate look-up table at the beginning
848 // static Bool_t ginit = kFALSE;
849 // static Float_t gnoise1,gnoise2,gnoise3;
850 // static Float_t ggg1[10000];
851 // static Float_t ggg2[10000];
852 // static Float_t ggg3[10000];
853 // static Float_t glandau1[10000];
854 // static Float_t glandau2[10000];
855 // static Float_t glandau3[10000];
857 // static Float_t gcor01[1000];
858 // static Float_t gcor02[1000];
859 // static Float_t gcorp[1000];
863 // if (ginit==kFALSE){
864 // for (Int_t i=1;i<1000;i++){
865 // Float_t rsigma = float(i)/100.;
866 // gcor02[i] = TMath::Max(0.81 +TMath::Exp(6.8*(rsigma-1.2)),0.6);
867 // gcor01[i] = TMath::Max(0.72 +TMath::Exp(2.04*(rsigma-1.2)),0.6);
868 // gcorp[i] = TMath::Max(TMath::Power((rsigma+0.5),1.5),1.2);
872 // for (Int_t i=3;i<10000;i++){
876 // Float_t amp = float(i);
877 // Float_t padlength =0.75;
878 // gnoise1 = 0.0004/padlength;
879 // Float_t nel = 0.268*amp;
880 // Float_t nprim = 0.155*amp;
881 // ggg1[i] = fkParam->GetDiffT()*fkParam->GetDiffT()*(2+0.001*nel/(padlength*padlength))/nel;
882 // glandau1[i] = (2.+0.12*nprim)*0.5* (2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
883 // if (glandau1[i]>1) glandau1[i]=1;
884 // glandau1[i]*=padlength*padlength/12.;
888 // gnoise2 = 0.0004/padlength;
890 // nprim = 0.133*amp;
891 // ggg2[i] = fkParam->GetDiffT()*fkParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
892 // glandau2[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
893 // if (glandau2[i]>1) glandau2[i]=1;
894 // glandau2[i]*=padlength*padlength/12.;
899 // gnoise3 = 0.0004/padlength;
901 // nprim = 0.133*amp;
902 // ggg3[i] = fkParam->GetDiffT()*fkParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
903 // glandau3[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
904 // if (glandau3[i]>1) glandau3[i]=1;
905 // glandau3[i]*=padlength*padlength/12.;
913 // Int_t amp = int(TMath::Abs(cl->GetQ()));
915 // seed->SetErrorY2(1.);
919 // Float_t z = TMath::Abs(fkParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
920 // Int_t ctype = cl->GetType();
921 // Float_t padlength= GetPadPitchLength(seed->GetRow());
923 // Double_t angle2 = seed->GetSnp()*seed->GetSnp();
924 // // if (angle2<0.6) angle2 = 0.6;
925 // angle2 = seed->GetTgl()*seed->GetTgl()*(1+angle2/(1-angle2));
927 // //cluster "quality"
928 // Int_t rsigmaz = int(100.*cl->GetSigmaZ2()/(seed->GetCurrentSigmaZ2()));
931 // if (fSectors==fInnerSec){
932 // snoise2 = gnoise1;
933 // res = ggg1[amp]*z+glandau1[amp]*angle2;
934 // if (ctype==0) res *= gcor01[rsigmaz];
937 // res*= gcorp[rsigmaz];
941 // if (padlength<1.1){
942 // snoise2 = gnoise2;
943 // res = ggg2[amp]*z+glandau2[amp]*angle2;
944 // if (ctype==0) res *= gcor02[rsigmaz];
947 // res*= gcorp[rsigmaz];
951 // snoise2 = gnoise3;
952 // res = ggg3[amp]*z+glandau3[amp]*angle2;
953 // if (ctype==0) res *= gcor02[rsigmaz];
956 // res*= gcorp[rsigmaz];
965 // if ((ctype<0) &&<70){
970 // if (res<2*snoise2)
972 // if (res>3) res =3;
973 // seed->SetErrorZ2(res);
981 void AliTPCtrackerMI::RotateToLocal(AliTPCseed *seed)
983 //rotate to track "local coordinata
984 Float_t x = seed->GetX();
985 Float_t y = seed->GetY();
986 Float_t ymax = x*TMath::Tan(0.5*fSectors->GetAlpha());
989 seed->SetRelativeSector((seed->GetRelativeSector()+1) % fN);
990 if (!seed->Rotate(fSectors->GetAlpha()))
992 } else if (y <-ymax) {
993 seed->SetRelativeSector((seed->GetRelativeSector()-1+fN) % fN);
994 if (!seed->Rotate(-fSectors->GetAlpha()))
1002 //_____________________________________________________________________________
1003 Double_t AliTPCtrackerMI::F1old(Double_t x1,Double_t y1,
1004 Double_t x2,Double_t y2,
1005 Double_t x3,Double_t y3) const
1007 //-----------------------------------------------------------------
1008 // Initial approximation of the track curvature
1009 //-----------------------------------------------------------------
1010 Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
1011 Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
1012 (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
1013 Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
1014 (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
1016 Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
1017 if ( xr*xr+yr*yr<=0.00000000000001) return 100;
1018 return -xr*yr/sqrt(xr*xr+yr*yr);
1023 //_____________________________________________________________________________
1024 Double_t AliTPCtrackerMI::F1(Double_t x1,Double_t y1,
1025 Double_t x2,Double_t y2,
1026 Double_t x3,Double_t y3) const
1028 //-----------------------------------------------------------------
1029 // Initial approximation of the track curvature
1030 //-----------------------------------------------------------------
1036 Double_t det = x3*y2-x2*y3;
1037 if (TMath::Abs(det)<1e-10){
1041 Double_t u = 0.5* (x2*(x2-x3)+y2*(y2-y3))/det;
1042 Double_t x0 = x3*0.5-y3*u;
1043 Double_t y0 = y3*0.5+x3*u;
1044 Double_t c2 = 1/TMath::Sqrt(x0*x0+y0*y0);
1050 Double_t AliTPCtrackerMI::F2(Double_t x1,Double_t y1,
1051 Double_t x2,Double_t y2,
1052 Double_t x3,Double_t y3) const
1054 //-----------------------------------------------------------------
1055 // Initial approximation of the track curvature
1056 //-----------------------------------------------------------------
1062 Double_t det = x3*y2-x2*y3;
1063 if (TMath::Abs(det)<1e-10) {
1067 Double_t u = 0.5* (x2*(x2-x3)+y2*(y2-y3))/det;
1068 Double_t x0 = x3*0.5-y3*u;
1069 Double_t y0 = y3*0.5+x3*u;
1070 Double_t c2 = 1/TMath::Sqrt(x0*x0+y0*y0);
1079 //_____________________________________________________________________________
1080 Double_t AliTPCtrackerMI::F2old(Double_t x1,Double_t y1,
1081 Double_t x2,Double_t y2,
1082 Double_t x3,Double_t y3) const
1084 //-----------------------------------------------------------------
1085 // Initial approximation of the track curvature times center of curvature
1086 //-----------------------------------------------------------------
1087 Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
1088 Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
1089 (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
1090 Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
1091 (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
1093 Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
1095 return -a/(d*y1-b)*xr/sqrt(xr*xr+yr*yr);
1098 //_____________________________________________________________________________
1099 Double_t AliTPCtrackerMI::F3(Double_t x1,Double_t y1,
1100 Double_t x2,Double_t y2,
1101 Double_t z1,Double_t z2) const
1103 //-----------------------------------------------------------------
1104 // Initial approximation of the tangent of the track dip angle
1105 //-----------------------------------------------------------------
1106 return (z1 - z2)/sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
1110 Double_t AliTPCtrackerMI::F3n(Double_t x1,Double_t y1,
1111 Double_t x2,Double_t y2,
1112 Double_t z1,Double_t z2, Double_t c) const
1114 //-----------------------------------------------------------------
1115 // Initial approximation of the tangent of the track dip angle
1116 //-----------------------------------------------------------------
1120 //angle1 = (z1-z2)*c/(TMath::ASin(c*x1-ni)-TMath::ASin(c*x2-ni));
1122 Double_t d = TMath::Sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
1123 if (TMath::Abs(d*c*0.5)>1) return 0;
1124 // Double_t angle2 = TMath::ASin(d*c*0.5);
1125 // Double_t angle2 = AliTPCFastMath::FastAsin(d*c*0.5);
1126 Double_t angle2 = (d*c*0.5>0.1)? TMath::ASin(d*c*0.5): AliTPCFastMath::FastAsin(d*c*0.5);
1128 angle2 = (z1-z2)*c/(angle2*2.);
1132 Bool_t AliTPCtrackerMI::GetProlongation(Double_t x1, Double_t x2, Double_t x[5], Double_t &y, Double_t &z) const
1133 {//-----------------------------------------------------------------
1134 // This function find proloncation of a track to a reference plane x=x2.
1135 //-----------------------------------------------------------------
1139 if (TMath::Abs(x[4]*x1 - x[2]) >= 0.999) {
1143 Double_t c1=x[4]*x1 - x[2], r1=TMath::Sqrt((1.-c1)*(1.+c1));
1144 Double_t c2=x[4]*x2 - x[2], r2=TMath::Sqrt((1.-c2)*(1.+c2));
1148 Double_t dy = dx*(c1+c2)/(r1+r2);
1151 Double_t delta = x[4]*dx*(c1+c2)/(c1*r2 + c2*r1);
1153 if (TMath::Abs(delta)>0.01){
1154 dz = x[3]*TMath::ASin(delta)/x[4];
1156 dz = x[3]*AliTPCFastMath::FastAsin(delta)/x[4];
1159 //dz = x[3]*AliTPCFastMath::FastAsin(delta)/x[4];
1167 Int_t AliTPCtrackerMI::LoadClusters (TTree *const tree)
1172 return LoadClusters();
1176 Int_t AliTPCtrackerMI::LoadClusters(const TObjArray *arr)
1179 // load clusters to the memory
1180 AliTPCClustersRow *clrow = new AliTPCClustersRow("AliTPCclusterMI");
1181 Int_t lower = arr->LowerBound();
1182 Int_t entries = arr->GetEntriesFast();
1184 for (Int_t i=lower; i<entries; i++) {
1185 clrow = (AliTPCClustersRow*) arr->At(i);
1186 if(!clrow) continue;
1187 if(!clrow->GetArray()) continue;
1191 fkParam->AdjustSectorRow(clrow->GetID(),sec,row);
1193 for (Int_t icl=0; icl<clrow->GetArray()->GetEntriesFast(); icl++){
1194 Transform((AliTPCclusterMI*)(clrow->GetArray()->At(icl)));
1197 if (clrow->GetArray()->GetEntriesFast()<=0) continue;
1198 AliTPCtrackerRow * tpcrow=0;
1201 tpcrow = &(fInnerSec[sec%fkNIS][row]);
1205 tpcrow = &(fOuterSec[(sec-fkNIS*2)%fkNOS][row]);
1206 left = (sec-fkNIS*2)/fkNOS;
1209 tpcrow->SetN1(clrow->GetArray()->GetEntriesFast());
1210 for (Int_t j=0;j<tpcrow->GetN1();++j)
1211 tpcrow->SetCluster1(j, *(AliTPCclusterMI*)(clrow->GetArray()->At(j)));
1214 tpcrow->SetN2(clrow->GetArray()->GetEntriesFast());
1215 for (Int_t j=0;j<tpcrow->GetN2();++j)
1216 tpcrow->SetCluster2(j, *(AliTPCclusterMI*)(clrow->GetArray()->At(j)));
1218 clrow->GetArray()->Clear("C");
1227 Int_t AliTPCtrackerMI::LoadClusters(const TClonesArray *arr)
1230 // load clusters to the memory from one
1233 AliTPCclusterMI *clust=0;
1234 Int_t count[72][96] = { {0} , {0} };
1236 // loop over clusters
1237 for (Int_t icl=0; icl<arr->GetEntriesFast(); icl++) {
1238 clust = (AliTPCclusterMI*)arr->At(icl);
1239 if(!clust) continue;
1240 //printf("cluster: det %d, row %d \n", clust->GetDetector(),clust->GetRow());
1242 // transform clusters
1245 // count clusters per pad row
1246 count[clust->GetDetector()][clust->GetRow()]++;
1249 // insert clusters to sectors
1250 for (Int_t icl=0; icl<arr->GetEntriesFast(); icl++) {
1251 clust = (AliTPCclusterMI*)arr->At(icl);
1252 if(!clust) continue;
1254 Int_t sec = clust->GetDetector();
1255 Int_t row = clust->GetRow();
1257 // filter overlapping pad rows needed by HLT
1258 if(sec<fkNIS*2) { //IROCs
1259 if(row == 30) continue;
1262 if(row == 27 || row == 76) continue;
1267 // left = sec/fkNIS;
1268 fInnerSec[sec%fkNIS].InsertCluster(clust, count[sec][row], fkParam);
1271 // left = (sec-fkNIS*2)/fkNOS;
1272 fOuterSec[(sec-fkNIS*2)%fkNOS].InsertCluster(clust, count[sec][row], fkParam);
1276 // Load functions must be called behind LoadCluster(TClonesArray*)
1278 //LoadOuterSectors();
1279 //LoadInnerSectors();
1285 Int_t AliTPCtrackerMI::LoadClusters()
1288 // load clusters to the memory
1289 static AliTPCClustersRow *clrow= new AliTPCClustersRow("AliTPCclusterMI");
1291 // TTree * tree = fClustersArray.GetTree();
1293 TTree * tree = fInput;
1294 TBranch * br = tree->GetBranch("Segment");
1295 br->SetAddress(&clrow);
1297 // Conversion of pad, row coordinates in local tracking coords.
1298 // Could be skipped here; is already done in clusterfinder
1300 Int_t j=Int_t(tree->GetEntries());
1301 for (Int_t i=0; i<j; i++) {
1305 fkParam->AdjustSectorRow(clrow->GetID(),sec,row);
1306 for (Int_t icl=0; icl<clrow->GetArray()->GetEntriesFast(); icl++){
1307 Transform((AliTPCclusterMI*)(clrow->GetArray()->At(icl)));
1310 AliTPCtrackerRow * tpcrow=0;
1313 tpcrow = &(fInnerSec[sec%fkNIS][row]);
1317 tpcrow = &(fOuterSec[(sec-fkNIS*2)%fkNOS][row]);
1318 left = (sec-fkNIS*2)/fkNOS;
1321 tpcrow->SetN1(clrow->GetArray()->GetEntriesFast());
1322 for (Int_t k=0;k<tpcrow->GetN1();++k)
1323 tpcrow->SetCluster1(k, *(AliTPCclusterMI*)(clrow->GetArray()->At(k)));
1326 tpcrow->SetN2(clrow->GetArray()->GetEntriesFast());
1327 for (Int_t k=0;k<tpcrow->GetN2();++k)
1328 tpcrow->SetCluster2(k, *(AliTPCclusterMI*)(clrow->GetArray()->At(k)));
1339 void AliTPCtrackerMI::UnloadClusters()
1342 // unload clusters from the memory
1344 Int_t nrows = fOuterSec->GetNRows();
1345 for (Int_t sec = 0;sec<fkNOS;sec++)
1346 for (Int_t row = 0;row<nrows;row++){
1347 AliTPCtrackerRow* tpcrow = &(fOuterSec[sec%fkNOS][row]);
1349 // if (tpcrow->fClusters1) delete []tpcrow->fClusters1;
1350 // if (tpcrow->fClusters2) delete []tpcrow->fClusters2;
1352 tpcrow->ResetClusters();
1355 nrows = fInnerSec->GetNRows();
1356 for (Int_t sec = 0;sec<fkNIS;sec++)
1357 for (Int_t row = 0;row<nrows;row++){
1358 AliTPCtrackerRow* tpcrow = &(fInnerSec[sec%fkNIS][row]);
1360 // if (tpcrow->fClusters1) delete []tpcrow->fClusters1;
1361 //if (tpcrow->fClusters2) delete []tpcrow->fClusters2;
1363 tpcrow->ResetClusters();
1369 void AliTPCtrackerMI::FillClusterArray(TObjArray* array) const{
1371 // Filling cluster to the array - For visualization purposes
1374 nrows = fOuterSec->GetNRows();
1375 for (Int_t sec = 0;sec<fkNOS;sec++)
1376 for (Int_t row = 0;row<nrows;row++){
1377 AliTPCtrackerRow* tpcrow = &(fOuterSec[sec%fkNOS][row]);
1378 if (!tpcrow) continue;
1379 for (Int_t icl = 0;icl<tpcrow->GetN();icl++){
1380 array->AddLast((TObject*)((*tpcrow)[icl]));
1383 nrows = fInnerSec->GetNRows();
1384 for (Int_t sec = 0;sec<fkNIS;sec++)
1385 for (Int_t row = 0;row<nrows;row++){
1386 AliTPCtrackerRow* tpcrow = &(fInnerSec[sec%fkNIS][row]);
1387 if (!tpcrow) continue;
1388 for (Int_t icl = 0;icl<tpcrow->GetN();icl++){
1389 array->AddLast((TObject*)(*tpcrow)[icl]);
1395 void AliTPCtrackerMI::Transform(AliTPCclusterMI * cluster){
1399 AliTPCcalibDB * calibDB = AliTPCcalibDB::Instance();
1400 AliTPCTransform *transform = calibDB->GetTransform() ;
1402 AliFatal("Tranformations not in calibDB");
1405 transform->SetCurrentRecoParam((AliTPCRecoParam*)AliTPCReconstructor::GetRecoParam());
1406 Double_t x[3]={cluster->GetRow(),cluster->GetPad(),cluster->GetTimeBin()};
1407 Int_t i[1]={cluster->GetDetector()};
1408 transform->Transform(x,i,0,1);
1409 // if (cluster->GetDetector()%36>17){
1414 // in debug mode check the transformation
1416 if (AliTPCReconstructor::StreamLevel()>2) {
1418 cluster->GetGlobalXYZ(gx);
1419 Int_t event = (fEvent==NULL)? 0: fEvent->GetEventNumberInFile();
1420 TTreeSRedirector &cstream = *fDebugStreamer;
1421 cstream<<"Transform"<<
1432 cluster->SetX(x[0]);
1433 cluster->SetY(x[1]);
1434 cluster->SetZ(x[2]);
1439 //if (!fkParam->IsGeoRead()) fkParam->ReadGeoMatrices();
1440 if (AliTPCReconstructor::GetRecoParam()->GetUseSectorAlignment() && (!calibDB->HasAlignmentOCDB())){
1441 TGeoHMatrix *mat = fkParam->GetClusterMatrix(cluster->GetDetector());
1443 Double_t pos[3]= {cluster->GetX(),cluster->GetY(),cluster->GetZ()};
1444 Double_t posC[3]={cluster->GetX(),cluster->GetY(),cluster->GetZ()};
1445 if (mat) mat->LocalToMaster(pos,posC);
1447 // chack Loading of Geo matrices from GeoManager - TEMPORARY FIX
1449 cluster->SetX(posC[0]);
1450 cluster->SetY(posC[1]);
1451 cluster->SetZ(posC[2]);
1455 //_____________________________________________________________________________
1456 Int_t AliTPCtrackerMI::LoadOuterSectors() {
1457 //-----------------------------------------------------------------
1458 // This function fills outer TPC sectors with clusters.
1459 //-----------------------------------------------------------------
1460 Int_t nrows = fOuterSec->GetNRows();
1462 for (Int_t sec = 0;sec<fkNOS;sec++)
1463 for (Int_t row = 0;row<nrows;row++){
1464 AliTPCtrackerRow* tpcrow = &(fOuterSec[sec%fkNOS][row]);
1465 Int_t sec2 = sec+2*fkNIS;
1467 Int_t ncl = tpcrow->GetN1();
1469 AliTPCclusterMI *c= (tpcrow->GetCluster1(ncl));
1470 index=(((sec2<<8)+row)<<16)+ncl;
1471 tpcrow->InsertCluster(c,index);
1474 ncl = tpcrow->GetN2();
1476 AliTPCclusterMI *c= (tpcrow->GetCluster2(ncl));
1477 index=((((sec2+fkNOS)<<8)+row)<<16)+ncl;
1478 tpcrow->InsertCluster(c,index);
1481 // write indexes for fast acces
1483 for (Int_t i=0;i<510;i++)
1484 tpcrow->SetFastCluster(i,-1);
1485 for (Int_t i=0;i<tpcrow->GetN();i++){
1486 Int_t zi = Int_t((*tpcrow)[i]->GetZ()+255.);
1487 tpcrow->SetFastCluster(zi,i); // write index
1490 for (Int_t i=0;i<510;i++){
1491 if (tpcrow->GetFastCluster(i)<0)
1492 tpcrow->SetFastCluster(i,last);
1494 last = tpcrow->GetFastCluster(i);
1503 //_____________________________________________________________________________
1504 Int_t AliTPCtrackerMI::LoadInnerSectors() {
1505 //-----------------------------------------------------------------
1506 // This function fills inner TPC sectors with clusters.
1507 //-----------------------------------------------------------------
1508 Int_t nrows = fInnerSec->GetNRows();
1510 for (Int_t sec = 0;sec<fkNIS;sec++)
1511 for (Int_t row = 0;row<nrows;row++){
1512 AliTPCtrackerRow* tpcrow = &(fInnerSec[sec%fkNIS][row]);
1515 Int_t ncl = tpcrow->GetN1();
1517 AliTPCclusterMI *c= (tpcrow->GetCluster1(ncl));
1518 index=(((sec<<8)+row)<<16)+ncl;
1519 tpcrow->InsertCluster(c,index);
1522 ncl = tpcrow->GetN2();
1524 AliTPCclusterMI *c= (tpcrow->GetCluster2(ncl));
1525 index=((((sec+fkNIS)<<8)+row)<<16)+ncl;
1526 tpcrow->InsertCluster(c,index);
1529 // write indexes for fast acces
1531 for (Int_t i=0;i<510;i++)
1532 tpcrow->SetFastCluster(i,-1);
1533 for (Int_t i=0;i<tpcrow->GetN();i++){
1534 Int_t zi = Int_t((*tpcrow)[i]->GetZ()+255.);
1535 tpcrow->SetFastCluster(zi,i); // write index
1538 for (Int_t i=0;i<510;i++){
1539 if (tpcrow->GetFastCluster(i)<0)
1540 tpcrow->SetFastCluster(i,last);
1542 last = tpcrow->GetFastCluster(i);
1554 //_________________________________________________________________________
1555 AliTPCclusterMI *AliTPCtrackerMI::GetClusterMI(Int_t index) const {
1556 //--------------------------------------------------------------------
1557 // Return pointer to a given cluster
1558 //--------------------------------------------------------------------
1559 if (index<0) return 0; // no cluster
1560 Int_t sec=(index&0xff000000)>>24;
1561 Int_t row=(index&0x00ff0000)>>16;
1562 Int_t ncl=(index&0x00007fff)>>00;
1564 const AliTPCtrackerRow * tpcrow=0;
1565 TClonesArray * clrow =0;
1567 if (sec<0 || sec>=fkNIS*4) {
1568 AliWarning(Form("Wrong sector %d",sec));
1573 AliTPCtrackerSector& tracksec = fInnerSec[sec%fkNIS];
1574 if (tracksec.GetNRows()<=row) return 0;
1575 tpcrow = &(tracksec[row]);
1576 if (tpcrow==0) return 0;
1579 if (tpcrow->GetN1()<=ncl) return 0;
1580 clrow = tpcrow->GetClusters1();
1583 if (tpcrow->GetN2()<=ncl) return 0;
1584 clrow = tpcrow->GetClusters2();
1588 AliTPCtrackerSector& tracksec = fOuterSec[(sec-fkNIS*2)%fkNOS];
1589 if (tracksec.GetNRows()<=row) return 0;
1590 tpcrow = &(tracksec[row]);
1591 if (tpcrow==0) return 0;
1593 if (sec-2*fkNIS<fkNOS) {
1594 if (tpcrow->GetN1()<=ncl) return 0;
1595 clrow = tpcrow->GetClusters1();
1598 if (tpcrow->GetN2()<=ncl) return 0;
1599 clrow = tpcrow->GetClusters2();
1603 return (AliTPCclusterMI*)clrow->At(ncl);
1609 Int_t AliTPCtrackerMI::FollowToNext(AliTPCseed& t, Int_t nr) {
1610 //-----------------------------------------------------------------
1611 // This function tries to find a track prolongation to next pad row
1612 //-----------------------------------------------------------------
1614 Double_t x= GetXrow(nr), ymax=GetMaxY(nr);
1617 AliTPCclusterMI *cl=0;
1618 Int_t tpcindex= t.GetClusterIndex2(nr);
1620 // update current shape info every 5 pad-row
1621 // if ( (nr%5==0) || t.GetNumberOfClusters()<2 || (t.fCurrentSigmaY2<0.0001) ){
1625 if (fIteration>0 && tpcindex>=-1){ //if we have already clusters
1627 if (tpcindex==-1) return 0; //track in dead zone
1628 if (tpcindex >= 0){ //
1629 cl = t.GetClusterPointer(nr);
1630 //if (cl==0) cl = GetClusterMI(tpcindex);
1631 if (!cl) cl = GetClusterMI(tpcindex);
1632 t.SetCurrentClusterIndex1(tpcindex);
1635 Int_t relativesector = ((tpcindex&0xff000000)>>24)%18; // if previously accepted cluster in different sector
1636 Float_t angle = relativesector*fSectors->GetAlpha()+fSectors->GetAlphaShift();
1638 if (angle<-TMath::Pi()) angle += 2*TMath::Pi();
1639 if (angle>=TMath::Pi()) angle -= 2*TMath::Pi();
1641 if (TMath::Abs(angle-t.GetAlpha())>0.001){
1642 Double_t rotation = angle-t.GetAlpha();
1643 t.SetRelativeSector(relativesector);
1644 if (!t.Rotate(rotation)) {
1645 t.SetClusterIndex(nr, t.GetClusterIndex(nr) | 0x8000);
1649 if (!t.PropagateTo(x)) {
1650 t.SetClusterIndex(nr, t.GetClusterIndex(nr) | 0x8000);
1654 t.SetCurrentCluster(cl);
1656 Int_t accept = AcceptCluster(&t,t.GetCurrentCluster());
1657 if ((tpcindex&0x8000)==0) accept =0;
1659 //if founded cluster is acceptible
1660 if (cl->IsUsed(11)) { // id cluster is shared inrease uncertainty
1661 t.SetErrorY2(t.GetErrorY2()+0.03);
1662 t.SetErrorZ2(t.GetErrorZ2()+0.03);
1663 t.SetErrorY2(t.GetErrorY2()*3);
1664 t.SetErrorZ2(t.GetErrorZ2()*3);
1666 t.SetNFoundable(t.GetNFoundable()+1);
1667 UpdateTrack(&t,accept);
1670 else { // Remove old cluster from track
1671 t.SetClusterIndex(nr, -3);
1672 t.SetClusterPointer(nr, 0);
1676 if (TMath::Abs(t.GetSnp())>AliTPCReconstructor::GetMaxSnpTracker()) return 0; // cut on angle
1677 if (fIteration>1 && IsFindable(t)){
1678 // not look for new cluster during refitting
1679 t.SetNFoundable(t.GetNFoundable()+1);
1684 // if (TMath::Abs(t.GetSnp())>0.95 || TMath::Abs(x*t.GetC()-t.GetEta())>0.95) return 0;// patch 28 fev 06
1685 if (!t.PropagateTo(x)) {
1686 if (fIteration==0) t.SetRemoval(10);
1689 Double_t y = t.GetY();
1690 if (TMath::Abs(y)>ymax){
1692 t.SetRelativeSector((t.GetRelativeSector()+1) % fN);
1693 if (!t.Rotate(fSectors->GetAlpha()))
1695 } else if (y <-ymax) {
1696 t.SetRelativeSector((t.GetRelativeSector()-1+fN) % fN);
1697 if (!t.Rotate(-fSectors->GetAlpha()))
1700 if (!t.PropagateTo(x)) {
1701 if (fIteration==0) t.SetRemoval(10);
1707 Double_t z=t.GetZ();
1710 if (!IsActive(t.GetRelativeSector(),nr)) {
1712 t.SetClusterIndex2(nr,-1);
1715 //AliInfo(Form("A - Sector%d phi %f - alpha %f", t.fRelativeSector,y/x, t.GetAlpha()));
1716 Bool_t isActive = IsActive(t.GetRelativeSector(),nr);
1717 Bool_t isActive2 = (nr>=fInnerSec->GetNRows()) ? fOuterSec[t.GetRelativeSector()][nr-fInnerSec->GetNRows()].GetN()>0:fInnerSec[t.GetRelativeSector()][nr].GetN()>0;
1719 if (!isActive || !isActive2) return 0;
1721 const AliTPCtrackerRow &krow=GetRow(t.GetRelativeSector(),nr);
1722 if ( (t.GetSigmaY2()<0) || t.GetSigmaZ2()<0) return 0;
1724 Double_t roadz = 1.;
1726 if (TMath::Abs(TMath::Abs(y)-ymax)<krow.GetDeadZone()){
1728 t.SetClusterIndex2(nr,-1);
1734 // if (TMath::Abs(z)<(AliTPCReconstructor::GetCtgRange()*x+10) && TMath::Abs(z)<fkParam->GetZLength(0) && (TMath::Abs(t.GetSnp())<AliTPCReconstructor::GetMaxSnpTracker()))
1735 t.SetNFoundable(t.GetNFoundable()+1);
1741 // cl = krow.FindNearest2(y+10.,z,roady,roadz,index);
1742 cl = krow.FindNearest2(y,z,roady,roadz,index);
1743 if (cl) t.SetCurrentClusterIndex1(krow.GetIndex(index));
1746 t.SetCurrentCluster(cl);
1748 if (fIteration==2&&cl->IsUsed(10)) return 0;
1749 Int_t accept = AcceptCluster(&t,t.GetCurrentCluster());
1750 if (fIteration==2&&cl->IsUsed(11)) {
1751 t.SetErrorY2(t.GetErrorY2()+0.03);
1752 t.SetErrorZ2(t.GetErrorZ2()+0.03);
1753 t.SetErrorY2(t.GetErrorY2()*3);
1754 t.SetErrorZ2(t.GetErrorZ2()*3);
1757 if (t.fCurrentCluster->IsUsed(10)){
1762 if (t.fNShared>0.7*t.GetNumberOfClusters()) {
1768 if (accept<3) UpdateTrack(&t,accept);
1771 if ( fIteration==0 && t.GetNFoundable()*0.5 > t.GetNumberOfClusters()) t.SetRemoval(10);
1779 //_________________________________________________________________________
1780 Bool_t AliTPCtrackerMI::GetTrackPoint(Int_t index, AliTrackPoint &p ) const
1782 // Get track space point by index
1783 // return false in case the cluster doesn't exist
1784 AliTPCclusterMI *cl = GetClusterMI(index);
1785 if (!cl) return kFALSE;
1786 Int_t sector = (index&0xff000000)>>24;
1787 // Int_t row = (index&0x00ff0000)>>16;
1789 // xyz[0] = fkParam->GetPadRowRadii(sector,row);
1790 xyz[0] = cl->GetX();
1791 xyz[1] = cl->GetY();
1792 xyz[2] = cl->GetZ();
1794 fkParam->AdjustCosSin(sector,cos,sin);
1795 Float_t x = cos*xyz[0]-sin*xyz[1];
1796 Float_t y = cos*xyz[1]+sin*xyz[0];
1798 Float_t sigmaY2 = 0.027*cl->GetSigmaY2();
1799 if (sector < fkParam->GetNInnerSector()) sigmaY2 *= 2.07;
1800 Float_t sigmaZ2 = 0.066*cl->GetSigmaZ2();
1801 if (sector < fkParam->GetNInnerSector()) sigmaZ2 *= 1.77;
1802 cov[0] = sin*sin*sigmaY2;
1803 cov[1] = -sin*cos*sigmaY2;
1805 cov[3] = cos*cos*sigmaY2;
1808 p.SetXYZ(x,y,xyz[2],cov);
1809 AliGeomManager::ELayerID iLayer;
1811 if (sector < fkParam->GetNInnerSector()) {
1812 iLayer = AliGeomManager::kTPC1;
1816 iLayer = AliGeomManager::kTPC2;
1817 idet = sector - fkParam->GetNInnerSector();
1819 UShort_t volid = AliGeomManager::LayerToVolUID(iLayer,idet);
1820 p.SetVolumeID(volid);
1826 Int_t AliTPCtrackerMI::UpdateClusters(AliTPCseed& t, Int_t nr) {
1827 //-----------------------------------------------------------------
1828 // This function tries to find a track prolongation to next pad row
1829 //-----------------------------------------------------------------
1830 t.SetCurrentCluster(0);
1831 t.SetCurrentClusterIndex1(-3);
1833 Double_t xt=t.GetX();
1834 Int_t row = GetRowNumber(xt)-1;
1835 Double_t ymax= GetMaxY(nr);
1837 if (row < nr) return 1; // don't prolongate if not information until now -
1838 // if (TMath::Abs(t.GetSnp())>0.9 && t.GetNumberOfClusters()>40. && fIteration!=2) {
1840 // return 0; // not prolongate strongly inclined tracks
1842 // if (TMath::Abs(t.GetSnp())>0.95) {
1844 // return 0; // not prolongate strongly inclined tracks
1845 // }// patch 28 fev 06
1847 Double_t x= GetXrow(nr);
1849 //t.PropagateTo(x+0.02);
1850 //t.PropagateTo(x+0.01);
1851 if (!t.PropagateTo(x)){
1858 if (TMath::Abs(y)>ymax){
1860 t.SetRelativeSector((t.GetRelativeSector()+1) % fN);
1861 if (!t.Rotate(fSectors->GetAlpha()))
1863 } else if (y <-ymax) {
1864 t.SetRelativeSector((t.GetRelativeSector()-1+fN) % fN);
1865 if (!t.Rotate(-fSectors->GetAlpha()))
1868 // if (!t.PropagateTo(x)){
1875 if (TMath::Abs(t.GetSnp())>AliTPCReconstructor::GetMaxSnpTracker()) return 0;
1877 if (!IsActive(t.GetRelativeSector(),nr)) {
1879 t.SetClusterIndex2(nr,-1);
1882 //AliInfo(Form("A - Sector%d phi %f - alpha %f", t.fRelativeSector,y/x, t.GetAlpha()));
1884 AliTPCtrackerRow &krow=GetRow(t.GetRelativeSector(),nr);
1886 if (TMath::Abs(TMath::Abs(y)-ymax)<krow.GetDeadZone()){
1888 t.SetClusterIndex2(nr,-1);
1894 // if (TMath::Abs(t.GetZ())<(AliTPCReconstructor::GetCtgRange()*t.GetX()+10) && (TMath::Abs(t.GetSnp())<AliTPCReconstructor::GetMaxSnpTracker()))
1895 if (IsFindable(t)) t.SetNFoundable(t.GetNFoundable()+1);
1901 if ( (nr%2==0) || t.GetNumberOfClusters()<2){
1902 // t.fCurrentSigmaY = GetSigmaY(&t);
1903 //t.fCurrentSigmaZ = GetSigmaZ(&t);
1907 AliTPCclusterMI *cl=0;
1910 Double_t roady = 1.;
1911 Double_t roadz = 1.;
1915 index = t.GetClusterIndex2(nr);
1916 if ( (index >= 0) && (index&0x8000)==0){
1917 cl = t.GetClusterPointer(nr);
1918 if ( (cl==0) && (index >= 0)) cl = GetClusterMI(index);
1919 t.SetCurrentClusterIndex1(index);
1921 t.SetCurrentCluster(cl);
1927 // if (index<0) return 0;
1928 UInt_t uindex = TMath::Abs(index);
1931 //cl = krow.FindNearest2(y+10,z,roady,roadz,uindex);
1932 cl = krow.FindNearest2(y,z,roady,roadz,uindex);
1935 if (cl) t.SetCurrentClusterIndex1(krow.GetIndex(uindex));
1936 t.SetCurrentCluster(cl);
1942 Int_t AliTPCtrackerMI::FollowToNextCluster(AliTPCseed & t, Int_t nr) {
1943 //-----------------------------------------------------------------
1944 // This function tries to find a track prolongation to next pad row
1945 //-----------------------------------------------------------------
1947 //update error according neighborhoud
1949 if (t.GetCurrentCluster()) {
1951 Int_t accept = AcceptCluster(&t,t.GetCurrentCluster());
1953 if (t.GetCurrentCluster()->IsUsed(10)){
1958 t.SetNShared(t.GetNShared()+1);
1959 if (t.GetNShared()>0.7*t.GetNumberOfClusters()) {
1964 if (fIteration>0) accept = 0;
1965 if (accept<3) UpdateTrack(&t,accept);
1969 if ( t.GetNumberOfClusters()>18 && ( (t.GetSigmaY2()+t.GetSigmaZ2())>0.16)) t.SetRemoval(10);
1970 if ( t.GetNumberOfClusters()>18 && t.GetChi2()/t.GetNumberOfClusters()>6 ) t.SetRemoval(10);
1972 if (( (t.GetNFoundable()*0.5 > t.GetNumberOfClusters()) || t.GetNoCluster()>15)) t.SetRemoval(10);
1980 //_____________________________________________________________________________
1981 Int_t AliTPCtrackerMI::FollowProlongation(AliTPCseed& t, Int_t rf, Int_t step, Bool_t fromSeeds) {
1982 //-----------------------------------------------------------------
1983 // This function tries to find a track prolongation.
1984 //-----------------------------------------------------------------
1985 Double_t xt=t.GetX();
1987 Double_t alpha=t.GetAlpha();
1988 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
1989 if (alpha < 0. ) alpha += 2.*TMath::Pi();
1991 t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN);
1993 Int_t first = GetRowNumber(xt);
1998 for (Int_t nr= first; nr>=rf; nr-=step) {
2000 if (t.GetKinkIndexes()[0]>0){
2001 for (Int_t i=0;i<3;i++){
2002 Int_t index = t.GetKinkIndexes()[i];
2003 if (index==0) break;
2004 if (index<0) continue;
2006 AliKink * kink = (AliKink*)fEvent->GetKink(index-1);
2008 printf("PROBLEM\n");
2011 Int_t kinkrow = kink->GetTPCRow0()+2+Int_t(0.5/(0.05+kink->GetAngle(2)));
2013 AliExternalTrackParam paramd(t);
2014 kink->SetDaughter(paramd);
2015 kink->SetStatus(2,5);
2022 if (nr==80) t.UpdateReference();
2023 if (nr<fInnerSec->GetNRows())
2024 fSectors = fInnerSec;
2026 fSectors = fOuterSec;
2027 if (FollowToNext(t,nr)==0)
2040 Int_t AliTPCtrackerMI::FollowBackProlongation(AliTPCseed& t, Int_t rf, Bool_t fromSeeds) {
2041 //-----------------------------------------------------------------
2042 // This function tries to find a track prolongation.
2043 //-----------------------------------------------------------------
2045 Double_t xt=t.GetX();
2046 Double_t alpha=t.GetAlpha();
2047 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
2048 if (alpha < 0. ) alpha += 2.*TMath::Pi();
2049 t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN);
2051 Int_t first = t.GetFirstPoint();
2052 Int_t ri = GetRowNumber(xt);
2056 if (first<ri) first = ri;
2058 if (first<0) first=0;
2059 for (Int_t nr=first; nr<=rf; nr++) {
2060 // if ( (TMath::Abs(t.GetSnp())>0.95)) break;//patch 28 fev 06
2061 if (t.GetKinkIndexes()[0]<0){
2062 for (Int_t i=0;i<3;i++){
2063 Int_t index = t.GetKinkIndexes()[i];
2064 if (index==0) break;
2065 if (index>0) continue;
2066 index = TMath::Abs(index);
2067 AliKink * kink = (AliKink*)fEvent->GetKink(index-1);
2069 printf("PROBLEM\n");
2072 Int_t kinkrow = kink->GetTPCRow0()-2-Int_t(0.5/(0.05+kink->GetAngle(2)));
2074 AliExternalTrackParam paramm(t);
2075 kink->SetMother(paramm);
2076 kink->SetStatus(2,1);
2083 if (nr<fInnerSec->GetNRows())
2084 fSectors = fInnerSec;
2086 fSectors = fOuterSec;
2097 Float_t AliTPCtrackerMI::OverlapFactor(AliTPCseed * s1, AliTPCseed * s2, Int_t &sum1, Int_t & sum2)
2099 // overlapping factor
2105 Float_t dz2 =(s1->GetZ() - s2->GetZ());
2108 Float_t dy2 =TMath::Abs((s1->GetY() - s2->GetY()));
2110 Float_t distance = TMath::Sqrt(dz2+dy2);
2111 if (distance>4.) return 0; // if there are far away - not overlap - to reduce combinatorics
2114 Int_t firstpoint = TMath::Min(s1->GetFirstPoint(),s2->GetFirstPoint());
2115 Int_t lastpoint = TMath::Max(s1->GetLastPoint(),s2->GetLastPoint());
2120 if (firstpoint>lastpoint) {
2121 firstpoint =lastpoint;
2126 for (Int_t i=firstpoint-1;i<lastpoint+1;i++){
2127 if (s1->GetClusterIndex2(i)>0) sum1++;
2128 if (s2->GetClusterIndex2(i)>0) sum2++;
2129 if (s1->GetClusterIndex2(i)==s2->GetClusterIndex2(i) && s1->GetClusterIndex2(i)>0) {
2133 if (sum<5) return 0;
2135 Float_t summin = TMath::Min(sum1+1,sum2+1);
2136 Float_t ratio = (sum+1)/Float_t(summin);
2140 void AliTPCtrackerMI::SignShared(AliTPCseed * s1, AliTPCseed * s2)
2144 Float_t thetaCut = 0.2;//+10.*TMath::Sqrt(s1->GetSigmaTglZ()+ s2->GetSigmaTglZ());
2145 if (TMath::Abs(s1->GetTgl()-s2->GetTgl())>thetaCut) return;
2146 Float_t minCl = TMath::Min(s1->GetNumberOfClusters(),s2->GetNumberOfClusters());
2147 Int_t cutN0 = TMath::Max(5,TMath::Nint(0.1*minCl));
2152 //Int_t firstpoint = TMath::Max(s1->GetFirstPoint(),s2->GetFirstPoint());
2153 //Int_t lastpoint = TMath::Min(s1->GetLastPoint(),s2->GetLastPoint());
2154 Int_t firstpoint = 0;
2155 Int_t lastpoint = 160;
2157 // if (firstpoint>=lastpoint-5) return;;
2159 for (Int_t i=firstpoint;i<lastpoint;i++){
2160 // if ( (s1->GetClusterIndex2(i)&0xFFFF8FFF)==(s2->GetClusterIndex2(i)&0xFFFF8FFF) && s1->GetClusterIndex2(i)>0) {
2161 if ( (s1->GetClusterIndex2(i))==(s2->GetClusterIndex2(i)) && s1->GetClusterIndex2(i)>0) {
2165 if (sumshared>cutN0){
2168 for (Int_t i=firstpoint;i<lastpoint;i++){
2169 // if ( (s1->GetClusterIndex2(i)&0xFFFF8FFF)==(s2->GetClusterIndex2(i)&0xFFFF8FFF) && s1->GetClusterIndex2(i)>0) {
2170 if ( (s1->GetClusterIndex2(i))==(s2->GetClusterIndex2(i)) && s1->GetClusterIndex2(i)>0) {
2171 AliTPCTrackerPoint *p1 = s1->GetTrackPoint(i);
2172 AliTPCTrackerPoint *p2 = s2->GetTrackPoint(i);;
2173 if (s1->IsActive()&&s2->IsActive()){
2174 p1->SetShared(kTRUE);
2175 p2->SetShared(kTRUE);
2181 if (sumshared>cutN0){
2182 for (Int_t i=0;i<4;i++){
2183 if (s1->GetOverlapLabel(3*i)==0){
2184 s1->SetOverlapLabel(3*i, s2->GetLabel());
2185 s1->SetOverlapLabel(3*i+1,sumshared);
2186 s1->SetOverlapLabel(3*i+2,s2->GetUniqueID());
2190 for (Int_t i=0;i<4;i++){
2191 if (s2->GetOverlapLabel(3*i)==0){
2192 s2->SetOverlapLabel(3*i, s1->GetLabel());
2193 s2->SetOverlapLabel(3*i+1,sumshared);
2194 s2->SetOverlapLabel(3*i+2,s1->GetUniqueID());
2201 void AliTPCtrackerMI::SignShared(TObjArray * arr)
2204 //sort trackss according sectors
2206 for (Int_t i=0; i<arr->GetEntriesFast(); i++) {
2207 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2209 //if (pt) RotateToLocal(pt);
2213 arr->Sort(); // sorting according relative sectors
2214 arr->Expand(arr->GetEntries());
2217 Int_t nseed=arr->GetEntriesFast();
2218 for (Int_t i=0; i<nseed; i++) {
2219 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2221 for (Int_t j=0;j<12;j++){
2222 pt->SetOverlapLabel(j,0);
2225 for (Int_t i=0; i<nseed; i++) {
2226 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2228 if (pt->GetRemoval()>10) continue;
2229 for (Int_t j=i+1; j<nseed; j++){
2230 AliTPCseed *pt2=(AliTPCseed*)arr->UncheckedAt(j);
2231 if (TMath::Abs(pt->GetRelativeSector()-pt2->GetRelativeSector())>1) continue;
2233 if (pt2->GetRemoval()<=10) {
2234 //if ( TMath::Abs(pt->GetRelativeSector()-pt2->GetRelativeSector())>0) break;
2242 void AliTPCtrackerMI::SortTracks(TObjArray * arr, Int_t mode) const
2245 //sort tracks in array according mode criteria
2246 Int_t nseed = arr->GetEntriesFast();
2247 for (Int_t i=0; i<nseed; i++) {
2248 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2259 void AliTPCtrackerMI::RemoveUsed2(TObjArray * arr, Float_t factor1, Float_t factor2, Int_t minimal)
2262 // Loop over all tracks and remove overlaped tracks (with lower quality)
2264 // 1. Unsign clusters
2265 // 2. Sort tracks according quality
2266 // Quality is defined by the number of cluster between first and last points
2268 // 3. Loop over tracks - decreasing quality order
2269 // a.) remove - If the fraction of shared cluster less than factor (1- n or 2)
2270 // b.) remove - If the minimal number of clusters less than minimal and not ITS
2271 // c.) if track accepted - sign clusters
2273 //Called in - AliTPCtrackerMI::Clusters2Tracks()
2274 // - AliTPCtrackerMI::PropagateBack()
2275 // - AliTPCtrackerMI::RefitInward()
2278 // factor1 - factor for constrained
2279 // factor2 - for non constrained tracks
2280 // if (Float_t(shared+1)/Float_t(found+1)>factor) - DELETE
2284 Int_t nseed = arr->GetEntriesFast();
2285 Float_t * quality = new Float_t[nseed];
2286 Int_t * indexes = new Int_t[nseed];
2290 for (Int_t i=0; i<nseed; i++) {
2291 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2296 pt->UpdatePoints(); //select first last max dens points
2297 Float_t * points = pt->GetPoints();
2298 if (points[3]<0.8) quality[i] =-1;
2299 quality[i] = (points[2]-points[0])+pt->GetNumberOfClusters();
2300 //prefer high momenta tracks if overlaps
2301 quality[i] *= TMath::Sqrt(TMath::Abs(pt->Pt())+0.5);
2303 TMath::Sort(nseed,quality,indexes);
2306 for (Int_t itrack=0; itrack<nseed; itrack++) {
2307 Int_t trackindex = indexes[itrack];
2308 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(trackindex);
2311 if (quality[trackindex]<0){
2312 MarkSeedFree( arr->RemoveAt(trackindex) );
2317 Int_t first = Int_t(pt->GetPoints()[0]);
2318 Int_t last = Int_t(pt->GetPoints()[2]);
2319 Double_t factor = (pt->GetBConstrain()) ? factor1: factor2;
2321 Int_t found,foundable,shared;
2322 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
2323 // pt->GetClusterStatistic(0,160, found, foundable,shared,kFALSE);
2324 Bool_t itsgold =kFALSE;
2327 if (pt->GetESD()->GetITSclusters(dummy)>4) itsgold= kTRUE;
2331 if (Float_t(shared+1)/Float_t(found+1)>factor){
2332 if (pt->GetKinkIndexes()[0]!=0) continue; //don't remove tracks - part of the kinks
2333 if( AliTPCReconstructor::StreamLevel()>3){
2334 TTreeSRedirector &cstream = *fDebugStreamer;
2335 cstream<<"RemoveUsed"<<
2336 "iter="<<fIteration<<
2340 MarkSeedFree( arr->RemoveAt(trackindex) );
2343 if (pt->GetNumberOfClusters()<50&&(found-0.5*shared)<minimal){ //remove short tracks
2344 if (pt->GetKinkIndexes()[0]!=0) continue; //don't remove tracks - part of the kinks
2345 if( AliTPCReconstructor::StreamLevel()>3){
2346 TTreeSRedirector &cstream = *fDebugStreamer;
2347 cstream<<"RemoveShort"<<
2348 "iter="<<fIteration<<
2352 MarkSeedFree( arr->RemoveAt(trackindex) );
2358 //if (sharedfactor>0.4) continue;
2359 if (pt->GetKinkIndexes()[0]>0) continue;
2360 //Remove tracks with undefined properties - seems
2361 if (pt->GetSigmaY2()<kAlmost0) continue; // ? what is the origin ?
2363 for (Int_t i=first; i<last; i++) {
2364 Int_t index=pt->GetClusterIndex2(i);
2365 // if (index<0 || index&0x8000 ) continue;
2366 if (index<0 || index&0x8000 ) continue;
2367 AliTPCclusterMI *c= pt->GetClusterPointer(i);
2374 Info("RemoveUsed","\n*****\nNumber of good tracks after shared removal\t%d\n",fNtracks);
2380 void AliTPCtrackerMI::DumpClusters(Int_t iter, TObjArray *trackArray)
2383 // Dump clusters after reco
2384 // signed and unsigned cluster can be visualized
2385 // 1. Unsign all cluster
2386 // 2. Sign all used clusters
2389 Int_t nseed = trackArray->GetEntries();
2390 for (Int_t i=0; i<nseed; i++){
2391 AliTPCseed *pt=(AliTPCseed*)trackArray->UncheckedAt(i);
2395 Bool_t isKink=pt->GetKinkIndex(0)!=0;
2396 for (Int_t j=0; j<160; ++j) {
2397 Int_t index=pt->GetClusterIndex2(j);
2398 if (index<0) continue;
2399 AliTPCclusterMI *c= pt->GetClusterPointer(j);
2401 if (isKink) c->Use(100); // kink
2402 c->Use(10); // by default usage 10
2407 for (Int_t sec=0;sec<fkNIS;sec++){
2408 for (Int_t row=0;row<fInnerSec->GetNRows();row++){
2409 TClonesArray *cla = fInnerSec[sec][row].GetClusters1();
2410 for (Int_t icl =0;icl< fInnerSec[sec][row].GetN1();icl++){
2411 AliTPCclusterMI* cl = (AliTPCclusterMI*)cla->At(icl);
2412 Float_t gx[3]; cl->GetGlobalXYZ(gx);
2413 (*fDebugStreamer)<<"clDump"<<
2421 cla = fInnerSec[sec][row].GetClusters2();
2422 for (Int_t icl =0;icl< fInnerSec[sec][row].GetN2();icl++){
2423 AliTPCclusterMI* cl = (AliTPCclusterMI*)cla->At(icl);
2424 Float_t gx[3]; cl->GetGlobalXYZ(gx);
2425 (*fDebugStreamer)<<"clDump"<<
2436 for (Int_t sec=0;sec<fkNOS;sec++){
2437 for (Int_t row=0;row<fOuterSec->GetNRows();row++){
2438 TClonesArray *cla = fOuterSec[sec][row].GetClusters1();
2439 for (Int_t icl =0;icl< fOuterSec[sec][row].GetN1();icl++){
2441 AliTPCclusterMI* cl = (AliTPCclusterMI*) cla->At(icl);
2442 cl->GetGlobalXYZ(gx);
2443 (*fDebugStreamer)<<"clDump"<<
2451 cla = fOuterSec[sec][row].GetClusters2();
2452 for (Int_t icl =0;icl< fOuterSec[sec][row].GetN2();icl++){
2454 AliTPCclusterMI* cl = (AliTPCclusterMI*) cla->At(icl);
2455 cl->GetGlobalXYZ(gx);
2456 (*fDebugStreamer)<<"clDump"<<
2468 void AliTPCtrackerMI::UnsignClusters()
2471 // loop over all clusters and unsign them
2474 for (Int_t sec=0;sec<fkNIS;sec++){
2475 for (Int_t row=0;row<fInnerSec->GetNRows();row++){
2476 TClonesArray *cla = fInnerSec[sec][row].GetClusters1();
2477 for (Int_t icl =0;icl< fInnerSec[sec][row].GetN1();icl++)
2478 // if (cl[icl].IsUsed(10))
2479 ((AliTPCclusterMI*) cla->At(icl))->Use(-1);
2480 cla = fInnerSec[sec][row].GetClusters2();
2481 for (Int_t icl =0;icl< fInnerSec[sec][row].GetN2();icl++)
2482 //if (cl[icl].IsUsed(10))
2483 ((AliTPCclusterMI*) cla->At(icl))->Use(-1);
2487 for (Int_t sec=0;sec<fkNOS;sec++){
2488 for (Int_t row=0;row<fOuterSec->GetNRows();row++){
2489 TClonesArray *cla = fOuterSec[sec][row].GetClusters1();
2490 for (Int_t icl =0;icl< fOuterSec[sec][row].GetN1();icl++)
2491 //if (cl[icl].IsUsed(10))
2492 ((AliTPCclusterMI*) cla->At(icl))->Use(-1);
2493 cla = fOuterSec[sec][row].GetClusters2();
2494 for (Int_t icl =0;icl< fOuterSec[sec][row].GetN2();icl++)
2495 //if (cl[icl].IsUsed(10))
2496 ((AliTPCclusterMI*) cla->At(icl))->Use(-1);
2504 void AliTPCtrackerMI::SignClusters(const TObjArray * arr, Float_t fnumber, Float_t fdensity)
2507 //sign clusters to be "used"
2509 // snumber and sdensity sign number of sigmas - bellow mean value to be accepted
2510 // loop over "primaries"
2524 Int_t nseed = arr->GetEntriesFast();
2525 for (Int_t i=0; i<nseed; i++) {
2526 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2530 if (!(pt->IsActive())) continue;
2531 Float_t dens = pt->GetNumberOfClusters()/Float_t(pt->GetNFoundable());
2532 if ( (dens>0.7) && (pt->GetNumberOfClusters()>70)){
2534 sumdens2+= dens*dens;
2535 sumn += pt->GetNumberOfClusters();
2536 sumn2 += pt->GetNumberOfClusters()*pt->GetNumberOfClusters();
2537 Float_t chi2 = pt->GetChi2()/pt->GetNumberOfClusters();
2540 sumchi2 +=chi2*chi2;
2545 Float_t mdensity = 0.9;
2546 Float_t meann = 130;
2547 Float_t meanchi = 1;
2548 Float_t sdensity = 0.1;
2549 Float_t smeann = 10;
2550 Float_t smeanchi =0.4;
2554 mdensity = sumdens/sum;
2556 meanchi = sumchi/sum;
2558 sdensity = sumdens2/sum-mdensity*mdensity;
2560 sdensity = TMath::Sqrt(sdensity);
2564 smeann = sumn2/sum-meann*meann;
2566 smeann = TMath::Sqrt(smeann);
2570 smeanchi = sumchi2/sum - meanchi*meanchi;
2572 smeanchi = TMath::Sqrt(smeanchi);
2578 //REMOVE SHORT DELTAS or tracks going out of sensitive volume of TPC
2580 for (Int_t i=0; i<nseed; i++) {
2581 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2585 if (pt->GetBSigned()) continue;
2586 if (pt->GetBConstrain()) continue;
2587 //if (!(pt->IsActive())) continue;
2589 Int_t found,foundable,shared;
2590 pt->GetClusterStatistic(0,160,found, foundable,shared);
2591 if (shared/float(found)>0.3) {
2592 if (shared/float(found)>0.9 ){
2593 //MarkSeedFree( arr->RemoveAt(i) );
2598 Bool_t isok =kFALSE;
2599 if ( (pt->GetNShared()/pt->GetNumberOfClusters()<0.5) &&pt->GetNumberOfClusters()>60)
2601 if ((TMath::Abs(1/pt->GetC())<100.) && (pt->GetNShared()/pt->GetNumberOfClusters()<0.7))
2603 if (TMath::Abs(pt->GetZ()/pt->GetX())>1.1)
2605 if ( (TMath::Abs(pt->GetSnp()>0.7) && pt->GetD(0,0)>60.))
2609 for (Int_t j=0; j<160; ++j) {
2610 Int_t index=pt->GetClusterIndex2(j);
2611 if (index<0) continue;
2612 AliTPCclusterMI *c= pt->GetClusterPointer(j);
2614 //if (!(c->IsUsed(10))) c->Use();
2621 Double_t maxchi = meanchi+2.*smeanchi;
2623 for (Int_t i=0; i<nseed; i++) {
2624 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2628 //if (!(pt->IsActive())) continue;
2629 if (pt->GetBSigned()) continue;
2630 Double_t chi = pt->GetChi2()/pt->GetNumberOfClusters();
2631 if (chi>maxchi) continue;
2634 Float_t dens = pt->GetNumberOfClusters()/Float_t(pt->GetNFoundable());
2636 //sign only tracks with enoug big density at the beginning
2638 if ((pt->GetDensityFirst(40)<0.75) && pt->GetNumberOfClusters()<meann) continue;
2641 Double_t mindens = TMath::Max(double(mdensity-sdensity*fdensity*bfactor),0.65);
2642 Double_t minn = TMath::Max(Int_t(meann-fnumber*smeann*bfactor),50);
2644 // if (pt->fBConstrain) mindens = TMath::Max(mdensity-sdensity*fdensity*bfactor,0.65);
2645 if ( (pt->GetRemoval()==10) && (pt->GetSnp()>0.8)&&(dens>mindens))
2648 if ((dens>mindens && pt->GetNumberOfClusters()>minn) && chi<maxchi ){
2649 //Int_t noc=pt->GetNumberOfClusters();
2650 pt->SetBSigned(kTRUE);
2651 for (Int_t j=0; j<160; ++j) {
2653 Int_t index=pt->GetClusterIndex2(j);
2654 if (index<0) continue;
2655 AliTPCclusterMI *c= pt->GetClusterPointer(j);
2657 // if (!(c->IsUsed(10))) c->Use();
2662 // gLastCheck = nseed;
2671 Int_t AliTPCtrackerMI::RefitInward(AliESDEvent *event)
2674 // back propagation of ESD tracks
2677 if (!event) return 0;
2678 const Int_t kMaxFriendTracks=2000;
2680 // extract correction object for multiplicity dependence of dEdx
2681 TObjArray * gainCalibArray = AliTPCcalibDB::Instance()->GetTimeGainSplinesRun(event->GetRunNumber());
2683 AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ;
2685 AliFatal("Tranformations not in RefitInward");
2688 transform->SetCurrentRecoParam((AliTPCRecoParam*)AliTPCReconstructor::GetRecoParam());
2689 const AliTPCRecoParam * recoParam = AliTPCcalibDB::Instance()->GetTransform()->GetCurrentRecoParam();
2690 Int_t nContribut = event->GetNumberOfTracks();
2691 TGraphErrors * graphMultDependenceDeDx = 0x0;
2692 if (recoParam && recoParam->GetUseMultiplicityCorrectionDedx() && gainCalibArray) {
2693 if (recoParam->GetUseTotCharge()) {
2694 graphMultDependenceDeDx = (TGraphErrors *) gainCalibArray->FindObject("TGRAPHERRORS_MEANQTOT_MULTIPLICITYDEPENDENCE_BEAM_ALL");
2696 graphMultDependenceDeDx = (TGraphErrors *) gainCalibArray->FindObject("TGRAPHERRORS_MEANQMAX_MULTIPLICITYDEPENDENCE_BEAM_ALL");
2702 //PrepareForProlongation(fSeeds,1);
2703 PropagateForward2(fSeeds);
2704 RemoveUsed2(fSeeds,0.4,0.4,20);
2706 TObjArray arraySeed(fSeeds->GetEntries());
2707 for (Int_t i=0;i<fSeeds->GetEntries();i++) {
2708 arraySeed.AddAt(fSeeds->At(i),i);
2710 SignShared(&arraySeed);
2711 // FindCurling(fSeeds, event,2); // find multi found tracks
2712 FindSplitted(fSeeds, event,2); // find multi found tracks
2713 if (AliTPCReconstructor::StreamLevel()>5) FindMultiMC(fSeeds, fEvent,2); // find multi found tracks
2716 Int_t nseed = fSeeds->GetEntriesFast();
2717 for (Int_t i=0;i<nseed;i++){
2718 AliTPCseed * seed = (AliTPCseed*) fSeeds->UncheckedAt(i);
2719 if (!seed) continue;
2720 if (seed->GetKinkIndex(0)>0) UpdateKinkQualityD(seed); // update quality informations for kinks
2721 AliESDtrack *esd=event->GetTrack(i);
2723 if (seed->GetNumberOfClusters()<60 && seed->GetNumberOfClusters()<(esd->GetTPCclusters(0) -5)*0.8){
2724 AliExternalTrackParam paramIn;
2725 AliExternalTrackParam paramOut;
2726 Int_t ncl = seed->RefitTrack(seed,¶mIn,¶mOut);
2727 if (AliTPCReconstructor::StreamLevel()>2) {
2728 (*fDebugStreamer)<<"RecoverIn"<<
2732 "pout.="<<¶mOut<<
2737 seed->Set(paramIn.GetX(),paramIn.GetAlpha(),paramIn.GetParameter(),paramIn.GetCovariance());
2738 seed->SetNumberOfClusters(ncl);
2742 seed->PropagateTo(fkParam->GetInnerRadiusLow());
2743 seed->UpdatePoints();
2744 AddCovariance(seed);
2745 MakeESDBitmaps(seed, esd);
2746 seed->CookdEdx(0.02,0.6);
2747 CookLabel(seed,0.1); //For comparison only
2749 if (AliTPCReconstructor::StreamLevel()>1 && seed!=0) {
2750 TTreeSRedirector &cstream = *fDebugStreamer;
2757 if (seed->GetNumberOfClusters()>15){
2758 esd->UpdateTrackParams(seed,AliESDtrack::kTPCrefit);
2759 esd->SetTPCPoints(seed->GetPoints());
2760 esd->SetTPCPointsF(seed->GetNFoundable());
2761 Int_t ndedx = seed->GetNCDEDX(0);
2762 Float_t sdedx = seed->GetSDEDX(0);
2763 Float_t dedx = seed->GetdEdx();
2764 // apply mutliplicity dependent dEdx correction if available
2765 if (graphMultDependenceDeDx) {
2766 Double_t corrGain = AliTPCcalibDButil::EvalGraphConst(graphMultDependenceDeDx, nContribut);
2767 dedx += (1 - corrGain)*50.; // MIP is normalized to 50
2769 esd->SetTPCsignal(dedx, sdedx, ndedx);
2771 // fill new dEdx information
2773 Double32_t signal[4];
2777 for(Int_t iarr=0;iarr<3;iarr++) {
2778 signal[iarr] = seed->GetDEDXregion(iarr+1);
2779 ncl[iarr] = seed->GetNCDEDX(iarr+1);
2780 nrows[iarr] = seed->GetNCDEDXInclThres(iarr+1);
2782 signal[3] = seed->GetDEDXregion(4);
2784 AliTPCdEdxInfo * infoTpcPid = new AliTPCdEdxInfo();
2785 infoTpcPid->SetTPCSignalRegionInfo(signal, ncl, nrows);
2786 esd->SetTPCdEdxInfo(infoTpcPid);
2788 // add seed to the esd track in Calib level
2790 Bool_t storeFriend = gRandom->Rndm()<(kMaxFriendTracks)/Float_t(nseed);
2791 if (AliTPCReconstructor::StreamLevel()>0 &&storeFriend){
2792 // RS: this is the only place where the seed is created not in the pool,
2793 // since it should belong to ESDevent
2794 AliTPCseed * seedCopy = new AliTPCseed(*seed, kTRUE);
2795 esd->AddCalibObject(seedCopy);
2800 //printf("problem\n");
2803 //FindKinks(fSeeds,event);
2804 if (AliTPCReconstructor::StreamLevel()>3) DumpClusters(2,fSeeds);
2805 Info("RefitInward","Number of refitted tracks %d",ntracks);
2807 AliCosmicTracker::FindCosmic(event, kTRUE);
2813 Int_t AliTPCtrackerMI::PropagateBack(AliESDEvent *event)
2816 // back propagation of ESD tracks
2818 if (!event) return 0;
2822 PropagateBack(fSeeds);
2823 RemoveUsed2(fSeeds,0.4,0.4,20);
2824 //FindCurling(fSeeds, fEvent,1);
2825 FindSplitted(fSeeds, event,1); // find multi found tracks
2826 if (AliTPCReconstructor::StreamLevel()>5) FindMultiMC(fSeeds, fEvent,1); // find multi found tracks
2829 Int_t nseed = fSeeds->GetEntriesFast();
2831 for (Int_t i=0;i<nseed;i++){
2832 AliTPCseed * seed = (AliTPCseed*) fSeeds->UncheckedAt(i);
2833 if (!seed) continue;
2834 if (seed->GetKinkIndex(0)<0) UpdateKinkQualityM(seed); // update quality informations for kinks
2835 seed->UpdatePoints();
2836 AddCovariance(seed);
2837 AliESDtrack *esd=event->GetTrack(i);
2838 if (!esd) continue; //never happen
2839 if (seed->GetNumberOfClusters()<60 && seed->GetNumberOfClusters()<(esd->GetTPCclusters(0) -5)*0.8){
2840 AliExternalTrackParam paramIn;
2841 AliExternalTrackParam paramOut;
2842 Int_t ncl = seed->RefitTrack(seed,¶mIn,¶mOut);
2843 if (AliTPCReconstructor::StreamLevel()>2) {
2844 (*fDebugStreamer)<<"RecoverBack"<<
2848 "pout.="<<¶mOut<<
2853 seed->Set(paramOut.GetX(),paramOut.GetAlpha(),paramOut.GetParameter(),paramOut.GetCovariance());
2854 seed->SetNumberOfClusters(ncl);
2857 seed->CookdEdx(0.02,0.6);
2858 CookLabel(seed,0.1); //For comparison only
2859 if (seed->GetNumberOfClusters()>15){
2860 esd->UpdateTrackParams(seed,AliESDtrack::kTPCout);
2861 esd->SetTPCPoints(seed->GetPoints());
2862 esd->SetTPCPointsF(seed->GetNFoundable());
2863 Int_t ndedx = seed->GetNCDEDX(0);
2864 Float_t sdedx = seed->GetSDEDX(0);
2865 Float_t dedx = seed->GetdEdx();
2866 esd->SetTPCsignal(dedx, sdedx, ndedx);
2868 Int_t eventnumber = event->GetEventNumberInFile();// patch 28 fev 06
2869 // This is most likely NOT the event number you'd like to use. It has nothing to do with the 'real' event number
2870 if (AliTPCReconstructor::StreamLevel()>1 && esd) {
2871 (*fDebugStreamer)<<"Cback"<<
2874 "EventNrInFile="<<eventnumber<<
2879 if (AliTPCReconstructor::StreamLevel()>3) DumpClusters(1,fSeeds);
2880 //FindKinks(fSeeds,event);
2881 Info("PropagateBack","Number of back propagated tracks %d",ntracks);
2888 Int_t AliTPCtrackerMI::PostProcess(AliESDEvent *event)
2891 // Post process events
2893 if (!event) return 0;
2896 // Set TPC event status
2899 // event affected by HV dip
2901 if(IsTPCHVDipEvent(event)) {
2902 event->ResetDetectorStatus(AliDAQ::kTPC);
2905 //printf("Status %d \n", event->IsDetectorOn(AliDAQ::kTPC));
2911 void AliTPCtrackerMI::DeleteSeeds()
2920 void AliTPCtrackerMI::ReadSeeds(const AliESDEvent *const event, Int_t direction)
2923 //read seeds from the event
2925 Int_t nentr=event->GetNumberOfTracks();
2927 Info("ReadSeeds", "Number of ESD tracks: %d\n", nentr);
2932 fSeeds = new TObjArray(nentr);
2936 for (Int_t i=0; i<nentr; i++) {
2937 AliESDtrack *esd=event->GetTrack(i);
2938 ULong_t status=esd->GetStatus();
2939 if (!(status&AliESDtrack::kTPCin)) continue;
2940 AliTPCtrack t(*esd);
2941 t.SetNumberOfClusters(0);
2942 // AliTPCseed *seed = new AliTPCseed(t,t.GetAlpha());
2943 AliTPCseed *seed = new( NextFreeSeed() ) AliTPCseed(t/*,t.GetAlpha()*/);
2944 seed->SetPoolID(fLastSeedID);
2945 seed->SetUniqueID(esd->GetID());
2946 AddCovariance(seed); //add systematic ucertainty
2947 for (Int_t ikink=0;ikink<3;ikink++) {
2948 Int_t index = esd->GetKinkIndex(ikink);
2949 seed->GetKinkIndexes()[ikink] = index;
2950 if (index==0) continue;
2951 index = TMath::Abs(index);
2952 AliESDkink * kink = fEvent->GetKink(index-1);
2953 if (kink&&esd->GetKinkIndex(ikink)<0){
2954 if ((status & AliESDtrack::kTRDrefit) != 0) kink->SetStatus(1,2);
2955 if ((status & AliESDtrack::kITSout) != 0) kink->SetStatus(1,0);
2957 if (kink&&esd->GetKinkIndex(ikink)>0){
2958 if ((status & AliESDtrack::kTRDrefit) != 0) kink->SetStatus(1,6);
2959 if ((status & AliESDtrack::kITSout) != 0) kink->SetStatus(1,4);
2963 if (((status&AliESDtrack::kITSout)==0)&&(direction==1)) seed->ResetCovariance(10.);
2964 if ( direction ==2 &&(status & AliESDtrack::kTRDrefit) == 0 ) seed->ResetCovariance(10.);
2965 //if ( direction ==2 && ((status & AliESDtrack::kTPCout) == 0) ) {
2966 // fSeeds->AddAt(0,i);
2967 // MarkSeedFree( seed );
2970 if ( direction ==2 &&(status & AliESDtrack::kTRDrefit) > 0 ) {
2971 Double_t par0[5],par1[5],alpha,x;
2972 esd->GetInnerExternalParameters(alpha,x,par0);
2973 esd->GetExternalParameters(x,par1);
2974 Double_t delta1 = TMath::Abs(par0[4]-par1[4])/(0.000000001+TMath::Abs(par0[4]+par1[4]));
2975 Double_t delta2 = TMath::Abs(par0[3]-par1[3]);
2977 if (esd->GetTRDncls()>0) trdchi2 = esd->GetTRDchi2()/esd->GetTRDncls();
2978 //reset covariance if suspicious
2979 if ( (delta1>0.1) || (delta2>0.006) ||trdchi2>7.)
2980 seed->ResetCovariance(10.);
2985 // rotate to the local coordinate system
2987 fSectors=fInnerSec; fN=fkNIS;
2988 Double_t alpha=seed->GetAlpha();
2989 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
2990 if (alpha < 0. ) alpha += 2.*TMath::Pi();
2991 Int_t ns=Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN;
2992 alpha =ns*fSectors->GetAlpha() + fSectors->GetAlphaShift();
2993 alpha-=seed->GetAlpha();
2994 if (alpha<-TMath::Pi()) alpha += 2*TMath::Pi();
2995 if (alpha>=TMath::Pi()) alpha -= 2*TMath::Pi();
2996 if (TMath::Abs(alpha) > 0.001) { // This should not happen normally
2997 AliWarning(Form("Rotating track over %f",alpha));
2998 if (!seed->Rotate(alpha)) {
2999 MarkSeedFree( seed );
3005 if (esd->GetKinkIndex(0)<=0){
3006 for (Int_t irow=0;irow<160;irow++){
3007 Int_t index = seed->GetClusterIndex2(irow);
3010 AliTPCclusterMI * cl = GetClusterMI(index);
3011 seed->SetClusterPointer(irow,cl);
3013 if ((index & 0x8000)==0){
3014 cl->Use(10); // accepted cluster
3016 cl->Use(6); // close cluster not accepted
3019 Info("ReadSeeds","Not found cluster");
3024 fSeeds->AddAt(seed,i);
3030 //_____________________________________________________________________________
3031 void AliTPCtrackerMI::MakeSeeds3(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, Float_t cuts[4],
3032 Float_t deltay, Int_t ddsec) {
3033 //-----------------------------------------------------------------
3034 // This function creates track seeds.
3035 // SEEDING WITH VERTEX CONSTRAIN
3036 //-----------------------------------------------------------------
3037 // cuts[0] - fP4 cut
3038 // cuts[1] - tan(phi) cut
3039 // cuts[2] - zvertex cut
3040 // cuts[3] - fP3 cut
3048 Double_t x[5], c[15];
3049 // Int_t di = i1-i2;
3051 AliTPCseed * seed = new( NextFreeSeed() ) AliTPCseed();
3052 seed->SetPoolID(fLastSeedID);
3053 Double_t alpha=fSectors->GetAlpha(), shift=fSectors->GetAlphaShift();
3054 Double_t cs=cos(alpha), sn=sin(alpha);
3056 // Double_t x1 =fOuterSec->GetX(i1);
3057 //Double_t xx2=fOuterSec->GetX(i2);
3059 Double_t x1 =GetXrow(i1);
3060 Double_t xx2=GetXrow(i2);
3062 Double_t x3=GetX(), y3=GetY(), z3=GetZ();
3064 Int_t imiddle = (i2+i1)/2; //middle pad row index
3065 Double_t xm = GetXrow(imiddle); // radius of middle pad-row
3066 const AliTPCtrackerRow& krm=GetRow(sec,imiddle); //middle pad -row
3070 const AliTPCtrackerRow& kr1=GetRow(ns,i1);
3071 Double_t ymax = GetMaxY(i1)-kr1.GetDeadZone()-1.5;
3072 Double_t ymaxm = GetMaxY(imiddle)-kr1.GetDeadZone()-1.5;
3075 // change cut on curvature if it can't reach this layer
3076 // maximal curvature set to reach it
3077 Double_t dvertexmax = TMath::Sqrt((x1-x3)*(x1-x3)+(ymax+5-y3)*(ymax+5-y3));
3078 if (dvertexmax*0.5*cuts[0]>0.85){
3079 cuts[0] = 0.85/(dvertexmax*0.5+1.);
3081 Double_t r2min = 1/(cuts[0]*cuts[0]); //minimal square of radius given by cut
3084 if (deltay>0) ddsec = 0;
3085 // loop over clusters
3086 for (Int_t is=0; is < kr1; is++) {
3088 if (kr1[is]->IsUsed(10)) continue;
3089 Double_t y1=kr1[is]->GetY(), z1=kr1[is]->GetZ();
3090 //if (TMath::Abs(y1)>ymax) continue;
3092 if (deltay>0 && TMath::Abs(ymax-TMath::Abs(y1))> deltay ) continue; // seed only at the edge
3094 // find possible directions
3095 Float_t anglez = (z1-z3)/(x1-x3);
3096 Float_t extraz = z1 - anglez*(x1-xx2); // extrapolated z
3099 //find rotation angles relative to line given by vertex and point 1
3100 Double_t dvertex2 = (x1-x3)*(x1-x3)+(y1-y3)*(y1-y3);
3101 Double_t dvertex = TMath::Sqrt(dvertex2);
3102 Double_t angle13 = TMath::ATan((y1-y3)/(x1-x3));
3103 Double_t cs13 = cos(-angle13), sn13 = sin(-angle13);
3106 // loop over 2 sectors
3112 Double_t dddz1=0; // direction of delta inclination in z axis
3119 for (Int_t dsec = dsec1; dsec<=dsec2;dsec++){
3120 Int_t sec2 = sec + dsec;
3122 // AliTPCtrackerRow& kr2 = fOuterSec[(sec2+fkNOS)%fkNOS][i2];
3123 //AliTPCtrackerRow& kr2m = fOuterSec[(sec2+fkNOS)%fkNOS][imiddle];
3124 AliTPCtrackerRow& kr2 = GetRow((sec2+fkNOS)%fkNOS,i2);
3125 AliTPCtrackerRow& kr2m = GetRow((sec2+fkNOS)%fkNOS,imiddle);
3126 Int_t index1 = TMath::Max(kr2.Find(extraz-0.6-dddz1*TMath::Abs(z1)*0.05)-1,0);
3127 Int_t index2 = TMath::Min(kr2.Find(extraz+0.6+dddz2*TMath::Abs(z1)*0.05)+1,kr2);
3129 // rotation angles to p1-p3
3130 Double_t cs13r = cos(-angle13+dsec*alpha)/dvertex, sn13r = sin(-angle13+dsec*alpha)/dvertex;
3131 Double_t x2, y2, z2;
3133 // Double_t dymax = maxangle*TMath::Abs(x1-xx2);
3136 Double_t dxx0 = (xx2-x3)*cs13r;
3137 Double_t dyy0 = (xx2-x3)*sn13r;
3138 for (Int_t js=index1; js < index2; js++) {
3139 const AliTPCclusterMI *kcl = kr2[js];
3140 if (kcl->IsUsed(10)) continue;
3142 //calcutate parameters
3144 Double_t yy0 = dyy0 +(kcl->GetY()-y3)*cs13r;
3146 if (TMath::Abs(yy0)<0.000001) continue;
3147 Double_t xx0 = dxx0 -(kcl->GetY()-y3)*sn13r;
3148 Double_t y0 = 0.5*(xx0*xx0+yy0*yy0-xx0)/yy0;
3149 Double_t r02 = (0.25+y0*y0)*dvertex2;
3150 //curvature (radius) cut
3151 if (r02<r2min) continue;
3155 Double_t c0 = 1/TMath::Sqrt(r02);
3159 //Double_t dfi0 = 2.*TMath::ASin(dvertex*c0*0.5);
3160 //Double_t dfi1 = 2.*TMath::ASin(TMath::Sqrt(yy0*yy0+(1-xx0)*(1-xx0))*dvertex*c0*0.5);
3161 Double_t dfi0 = 2.*AliTPCFastMath::FastAsin(dvertex*c0*0.5);
3162 Double_t dfi1 = 2.*AliTPCFastMath::FastAsin(TMath::Sqrt(yy0*yy0+(1-xx0)*(1-xx0))*dvertex*c0*0.5);
3165 Double_t z0 = kcl->GetZ();
3166 Double_t zzzz2 = z1-(z1-z3)*dfi1/dfi0;
3167 if (TMath::Abs(zzzz2-z0)>0.5) continue;
3170 Double_t dip = (z1-z0)*c0/dfi1;
3171 Double_t x0 = (0.5*cs13+y0*sn13)*dvertex*c0;
3182 x2= xx2*cs-y2*sn*dsec;
3183 y2=+xx2*sn*dsec+y2*cs;
3193 // do we have cluster at the middle ?
3195 GetProlongation(x1,xm,x,ym,zm);
3197 AliTPCclusterMI * cm=0;
3198 if (TMath::Abs(ym)-ymaxm<0){
3199 cm = krm.FindNearest2(ym,zm,1.0,0.6,dummy);
3200 if ((!cm) || (cm->IsUsed(10))) {
3205 // rotate y1 to system 0
3206 // get state vector in rotated system
3207 Double_t yr1 = (-0.5*sn13+y0*cs13)*dvertex*c0;
3208 Double_t xr2 = x0*cs+yr1*sn*dsec;
3209 Double_t xr[5]={kcl->GetY(),kcl->GetZ(), xr2, dip, c0};
3211 GetProlongation(xx2,xm,xr,ym,zm);
3212 if (TMath::Abs(ym)-ymaxm<0){
3213 cm = kr2m.FindNearest2(ym,zm,1.0,0.6,dummy);
3214 if ((!cm) || (cm->IsUsed(10))) {
3221 // Double_t dym = 0;
3222 // Double_t dzm = 0;
3224 // dym = ym - cm->GetY();
3225 // dzm = zm - cm->GetZ();
3232 Double_t sy1=kr1[is]->GetSigmaY2()*2., sz1=kr1[is]->GetSigmaZ2()*2.;
3233 Double_t sy2=kcl->GetSigmaY2()*2., sz2=kcl->GetSigmaZ2()*2.;
3234 //Double_t sy3=400*3./12., sy=0.1, sz=0.1;
3235 Double_t sy3=25000*x[4]*x[4]+0.1, sy=0.1, sz=0.1;
3236 //Double_t sy3=25000*x[4]*x[4]*60+0.5, sy=0.1, sz=0.1;
3238 Double_t f40=(F1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
3239 Double_t f42=(F1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
3240 Double_t f43=(F1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
3241 Double_t f20=(F2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
3242 Double_t f22=(F2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
3243 Double_t f23=(F2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
3245 Double_t f30=(F3(x1,y1+sy,x2,y2,z1,z2)-x[3])/sy;
3246 Double_t f31=(F3(x1,y1,x2,y2,z1+sz,z2)-x[3])/sz;
3247 Double_t f32=(F3(x1,y1,x2,y2+sy,z1,z2)-x[3])/sy;
3248 Double_t f34=(F3(x1,y1,x2,y2,z1,z2+sz)-x[3])/sz;
3252 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3253 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3254 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3255 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3256 c[13]=f30*sy1*f40+f32*sy2*f42;
3257 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3259 // if (!BuildSeed(kr1[is],kcl,0,x1,x2,x3,x,c)) continue;
3261 UInt_t index=kr1.GetIndex(is);
3262 if (seed) {MarkSeedFree(seed); seed = 0;}
3263 AliTPCseed *track = seed = new( NextFreeSeed() ) AliTPCseed(x1, ns*alpha+shift, x, c, index);
3264 seed->SetPoolID(fLastSeedID);
3265 track->SetIsSeeding(kTRUE);
3266 track->SetSeed1(i1);
3267 track->SetSeed2(i2);
3268 track->SetSeedType(3);
3272 FollowProlongation(*track, (i1+i2)/2,1);
3273 Int_t foundable,found,shared;
3274 track->GetClusterStatistic((i1+i2)/2,i1, found, foundable, shared, kTRUE);
3275 if ((found<0.55*foundable) || shared>0.5*found || (track->GetSigmaY2()+track->GetSigmaZ2())>0.5){
3276 MarkSeedFree(seed); seed = 0;
3282 FollowProlongation(*track, i2,1);
3286 track->SetBConstrain(1);
3287 // track->fLastPoint = i1+fInnerSec->GetNRows(); // first cluster in track position
3288 track->SetLastPoint(i1); // first cluster in track position
3289 track->SetFirstPoint(track->GetLastPoint());
3291 if (track->GetNumberOfClusters()<(i1-i2)*0.5 ||
3292 track->GetNumberOfClusters() < track->GetNFoundable()*0.6 ||
3293 track->GetNShared()>0.4*track->GetNumberOfClusters() ) {
3294 MarkSeedFree(seed); seed = 0;
3298 // Z VERTEX CONDITION
3299 Double_t zv, bz=GetBz();
3300 if ( !track->GetZAt(0.,bz,zv) ) continue;
3301 if (TMath::Abs(zv-z3)>cuts[2]) {
3302 FollowProlongation(*track, TMath::Max(i2-20,0));
3303 if ( !track->GetZAt(0.,bz,zv) ) continue;
3304 if (TMath::Abs(zv-z3)>cuts[2]){
3305 FollowProlongation(*track, TMath::Max(i2-40,0));
3306 if ( !track->GetZAt(0.,bz,zv) ) continue;
3307 if (TMath::Abs(zv-z3)>cuts[2] &&(track->GetNumberOfClusters() > track->GetNFoundable()*0.7)){
3308 // make seed without constrain
3309 AliTPCseed * track2 = MakeSeed(track,0.2,0.5,1.);
3310 FollowProlongation(*track2, i2,1);
3311 track2->SetBConstrain(kFALSE);
3312 track2->SetSeedType(1);
3313 arr->AddLast(track2);
3314 MarkSeedFree( seed ); seed = 0;
3318 MarkSeedFree( seed ); seed = 0;
3325 track->SetSeedType(0);
3326 arr->AddLast(track); // note, track is seed, don't free the seed
3327 seed = new( NextFreeSeed() ) AliTPCseed;
3328 seed->SetPoolID(fLastSeedID);
3330 // don't consider other combinations
3331 if (track->GetNumberOfClusters() > track->GetNFoundable()*0.8)
3337 Info("MakeSeeds3","\nSeeding statistic:\t%d\t%d\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin,nout1,nout2);
3339 if (seed) MarkSeedFree( seed );
3343 void AliTPCtrackerMI::MakeSeeds5(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, Float_t cuts[4],
3348 //-----------------------------------------------------------------
3349 // This function creates track seeds.
3350 //-----------------------------------------------------------------
3351 // cuts[0] - fP4 cut
3352 // cuts[1] - tan(phi) cut
3353 // cuts[2] - zvertex cut
3354 // cuts[3] - fP3 cut
3364 Double_t x[5], c[15];
3366 // make temporary seed
3367 AliTPCseed * seed = new( NextFreeSeed() ) AliTPCseed;
3368 seed->SetPoolID(fLastSeedID);
3369 Double_t alpha=fOuterSec->GetAlpha(), shift=fOuterSec->GetAlphaShift();
3370 // Double_t cs=cos(alpha), sn=sin(alpha);
3375 Double_t x1 = GetXrow(i1-1);
3376 const AliTPCtrackerRow& kr1=GetRow(sec,i1-1);
3377 Double_t y1max = GetMaxY(i1-1)-kr1.GetDeadZone()-1.5;
3379 Double_t x1p = GetXrow(i1);
3380 const AliTPCtrackerRow& kr1p=GetRow(sec,i1);
3382 Double_t x1m = GetXrow(i1-2);
3383 const AliTPCtrackerRow& kr1m=GetRow(sec,i1-2);
3386 //last 3 padrow for seeding
3387 AliTPCtrackerRow& kr3 = GetRow((sec+fkNOS)%fkNOS,i1-7);
3388 Double_t x3 = GetXrow(i1-7);
3389 // Double_t y3max= GetMaxY(i1-7)-kr3.fDeadZone-1.5;
3391 AliTPCtrackerRow& kr3p = GetRow((sec+fkNOS)%fkNOS,i1-6);
3392 Double_t x3p = GetXrow(i1-6);
3394 AliTPCtrackerRow& kr3m = GetRow((sec+fkNOS)%fkNOS,i1-8);
3395 Double_t x3m = GetXrow(i1-8);
3400 Int_t im = i1-4; //middle pad row index
3401 Double_t xm = GetXrow(im); // radius of middle pad-row
3402 const AliTPCtrackerRow& krm=GetRow(sec,im); //middle pad -row
3403 // Double_t ymmax = GetMaxY(im)-kr1.fDeadZone-1.5;
3406 Double_t deltax = x1-x3;
3407 Double_t dymax = deltax*cuts[1];
3408 Double_t dzmax = deltax*cuts[3];
3410 // loop over clusters
3411 for (Int_t is=0; is < kr1; is++) {
3413 if (kr1[is]->IsUsed(10)) continue;
3414 Double_t y1=kr1[is]->GetY(), z1=kr1[is]->GetZ();
3416 if (deltay>0 && TMath::Abs(y1max-TMath::Abs(y1))> deltay ) continue; // seed only at the edge
3418 Int_t index1 = TMath::Max(kr3.Find(z1-dzmax)-1,0);
3419 Int_t index2 = TMath::Min(kr3.Find(z1+dzmax)+1,kr3);
3425 for (Int_t js=index1; js < index2; js++) {
3426 const AliTPCclusterMI *kcl = kr3[js];
3427 if (kcl->IsUsed(10)) continue;
3429 // apply angular cuts
3430 if (TMath::Abs(y1-y3)>dymax) continue;
3433 if (TMath::Abs(z1-z3)>dzmax) continue;
3435 Double_t angley = (y1-y3)/(x1-x3);
3436 Double_t anglez = (z1-z3)/(x1-x3);
3438 Double_t erry = TMath::Abs(angley)*(x1-x1m)*0.5+0.5;
3439 Double_t errz = TMath::Abs(anglez)*(x1-x1m)*0.5+0.5;
3441 Double_t yyym = angley*(xm-x1)+y1;
3442 Double_t zzzm = anglez*(xm-x1)+z1;
3444 const AliTPCclusterMI *kcm = krm.FindNearest2(yyym,zzzm,erry,errz,index);
3446 if (kcm->IsUsed(10)) continue;
3448 erry = TMath::Abs(angley)*(x1-x1m)*0.4+0.5;
3449 errz = TMath::Abs(anglez)*(x1-x1m)*0.4+0.5;
3456 // look around first
3457 const AliTPCclusterMI *kc1m = kr1m.FindNearest2(angley*(x1m-x1)+y1,
3463 if (kc1m->IsUsed(10)) used++;
3465 const AliTPCclusterMI *kc1p = kr1p.FindNearest2(angley*(x1p-x1)+y1,
3471 if (kc1p->IsUsed(10)) used++;
3473 if (used>1) continue;
3474 if (found<1) continue;
3478 const AliTPCclusterMI *kc3m = kr3m.FindNearest2(angley*(x3m-x3)+y3,
3484 if (kc3m->IsUsed(10)) used++;
3488 const AliTPCclusterMI *kc3p = kr3p.FindNearest2(angley*(x3p-x3)+y3,
3494 if (kc3p->IsUsed(10)) used++;
3498 if (used>1) continue;
3499 if (found<3) continue;
3509 x[4]=F1(x1,y1,x2,y2,x3,y3);
3510 //if (TMath::Abs(x[4]) >= cuts[0]) continue;
3513 x[2]=F2(x1,y1,x2,y2,x3,y3);
3516 x[3]=F3n(x1,y1,x2,y2,z1,z2,x[4]);
3517 //if (TMath::Abs(x[3]) > cuts[3]) continue;
3521 Double_t sy1=0.1, sz1=0.1;
3522 Double_t sy2=0.1, sz2=0.1;
3523 Double_t sy3=0.1, sy=0.1, sz=0.1;
3525 Double_t f40=(F1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
3526 Double_t f42=(F1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
3527 Double_t f43=(F1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
3528 Double_t f20=(F2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
3529 Double_t f22=(F2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
3530 Double_t f23=(F2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
3532 Double_t f30=(F3(x1,y1+sy,x2,y2,z1,z2)-x[3])/sy;
3533 Double_t f31=(F3(x1,y1,x2,y2,z1+sz,z2)-x[3])/sz;
3534 Double_t f32=(F3(x1,y1,x2,y2+sy,z1,z2)-x[3])/sy;
3535 Double_t f34=(F3(x1,y1,x2,y2,z1,z2+sz)-x[3])/sz;
3539 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3540 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3541 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3542 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3543 c[13]=f30*sy1*f40+f32*sy2*f42;
3544 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3546 // if (!BuildSeed(kr1[is],kcl,0,x1,x2,x3,x,c)) continue;
3548 index=kr1.GetIndex(is);
3549 if (seed) {MarkSeedFree( seed ); seed = 0;}
3550 AliTPCseed *track = seed = new( NextFreeSeed() ) AliTPCseed(x1, sec*alpha+shift, x, c, index);
3551 seed->SetPoolID(fLastSeedID);
3553 track->SetIsSeeding(kTRUE);
3556 FollowProlongation(*track, i1-7,1);
3557 if (track->GetNumberOfClusters() < track->GetNFoundable()*0.75 ||
3558 track->GetNShared()>0.6*track->GetNumberOfClusters() || ( track->GetSigmaY2()+ track->GetSigmaZ2())>0.6){
3559 MarkSeedFree( seed ); seed = 0;
3565 FollowProlongation(*track, i2,1);
3566 track->SetBConstrain(0);
3567 track->SetLastPoint(i1+fInnerSec->GetNRows()); // first cluster in track position
3568 track->SetFirstPoint(track->GetLastPoint());
3570 if (track->GetNumberOfClusters()<(i1-i2)*0.5 ||
3571 track->GetNumberOfClusters()<track->GetNFoundable()*0.7 ||
3572 track->GetNShared()>2. || track->GetChi2()/track->GetNumberOfClusters()>6 || ( track->GetSigmaY2()+ track->GetSigmaZ2())>0.5 ) {
3573 MarkSeedFree( seed ); seed = 0;
3578 FollowProlongation(*track, TMath::Max(i2-10,0),1);
3579 AliTPCseed * track2 = MakeSeed(track,0.2,0.5,0.9);
3580 FollowProlongation(*track2, i2,1);
3581 track2->SetBConstrain(kFALSE);
3582 track2->SetSeedType(4);
3583 arr->AddLast(track2);
3584 MarkSeedFree( seed ); seed = 0;
3588 //arr->AddLast(track);
3589 //seed = new AliTPCseed;
3595 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);
3597 if (seed) MarkSeedFree(seed);
3601 //_____________________________________________________________________________
3602 void AliTPCtrackerMI::MakeSeeds2(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, Float_t */*cuts[4]*/,
3603 Float_t deltay, Bool_t /*bconstrain*/) {
3604 //-----------------------------------------------------------------
3605 // This function creates track seeds - without vertex constraint
3606 //-----------------------------------------------------------------
3607 // cuts[0] - fP4 cut - not applied
3608 // cuts[1] - tan(phi) cut
3609 // cuts[2] - zvertex cut - not applied
3610 // cuts[3] - fP3 cut
3620 Double_t alpha=fOuterSec->GetAlpha(), shift=fOuterSec->GetAlphaShift();
3621 // Double_t cs=cos(alpha), sn=sin(alpha);
3622 Int_t row0 = (i1+i2)/2;
3623 Int_t drow = (i1-i2)/2;
3624 const AliTPCtrackerRow& kr0=fSectors[sec][row0];
3625 AliTPCtrackerRow * kr=0;
3627 AliTPCpolyTrack polytrack;
3628 Int_t nclusters=fSectors[sec][row0];
3629 AliTPCseed * seed = new( NextFreeSeed() ) AliTPCseed;
3630 seed->SetPoolID(fLastSeedID);
3635 for (Int_t is=0; is < nclusters; is++) { //LOOP over clusters
3637 Int_t nfoundable =0;
3638 for (Int_t iter =1; iter<2; iter++){ //iterations
3639 const AliTPCtrackerRow& krm=fSectors[sec][row0-iter];
3640 const AliTPCtrackerRow& krp=fSectors[sec][row0+iter];
3641 const AliTPCclusterMI * cl= kr0[is];
3643 if (cl->IsUsed(10)) {
3649 Double_t x = kr0.GetX();
3650 // Initialization of the polytrack
3655 Double_t y0= cl->GetY();
3656 Double_t z0= cl->GetZ();
3660 Double_t ymax = fSectors->GetMaxY(row0)-kr0.GetDeadZone()-1.5;
3661 if (deltay>0 && TMath::Abs(ymax-TMath::Abs(y0))> deltay ) continue; // seed only at the edge
3663 erry = (0.5)*cl->GetSigmaY2()/TMath::Sqrt(cl->GetQ())*6;
3664 errz = (0.5)*cl->GetSigmaZ2()/TMath::Sqrt(cl->GetQ())*6;
3665 polytrack.AddPoint(x,y0,z0,erry, errz);
3668 if (cl->IsUsed(10)) sumused++;
3671 Float_t roady = (5*TMath::Sqrt(cl->GetSigmaY2()+0.2)+1.)*iter;
3672 Float_t roadz = (5*TMath::Sqrt(cl->GetSigmaZ2()+0.2)+1.)*iter;
3675 AliTPCclusterMI * cl1 = krm.FindNearest(y0,z0,roady,roadz);
3676 if (cl1 && TMath::Abs(ymax-TMath::Abs(y0))) {
3677 erry = (0.5)*cl1->GetSigmaY2()/TMath::Sqrt(cl1->GetQ())*3;
3678 errz = (0.5)*cl1->GetSigmaZ2()/TMath::Sqrt(cl1->GetQ())*3;
3679 if (cl1->IsUsed(10)) sumused++;
3680 polytrack.AddPoint(x,cl1->GetY(),cl1->GetZ(),erry,errz);
3684 AliTPCclusterMI * cl2 = krp.FindNearest(y0,z0,roady,roadz);
3686 erry = (0.5)*cl2->GetSigmaY2()/TMath::Sqrt(cl2->GetQ())*3;
3687 errz = (0.5)*cl2->GetSigmaZ2()/TMath::Sqrt(cl2->GetQ())*3;
3688 if (cl2->IsUsed(10)) sumused++;
3689 polytrack.AddPoint(x,cl2->GetY(),cl2->GetZ(),erry,errz);
3692 if (sumused>0) continue;
3694 polytrack.UpdateParameters();
3700 nfoundable = polytrack.GetN();
3701 nfound = nfoundable;
3703 for (Int_t ddrow = iter+1; ddrow<drow;ddrow++){
3704 Float_t maxdist = 0.8*(1.+3./(ddrow));
3705 for (Int_t delta = -1;delta<=1;delta+=2){
3706 Int_t row = row0+ddrow*delta;
3707 kr = &(fSectors[sec][row]);
3708 Double_t xn = kr->GetX();
3709 Double_t ymax1 = fSectors->GetMaxY(row)-kr->GetDeadZone()-1.5;
3710 polytrack.GetFitPoint(xn,yn,zn);
3711 if (TMath::Abs(yn)>ymax1) continue;
3713 AliTPCclusterMI * cln = kr->FindNearest(yn,zn,roady,roadz);
3715 Float_t dist = TMath::Sqrt( (yn-cln->GetY())*(yn-cln->GetY())+(zn-cln->GetZ())*(zn-cln->GetZ()));
3718 erry = (dist+0.3)*cln->GetSigmaY2()/TMath::Sqrt(cln->GetQ())*(1.+1./(ddrow));
3719 errz = (dist+0.3)*cln->GetSigmaZ2()/TMath::Sqrt(cln->GetQ())*(1.+1./(ddrow));
3720 if (cln->IsUsed(10)) {
3721 // printf("used\n");
3729 polytrack.AddPoint(xn,cln->GetY(),cln->GetZ(),erry, errz);
3734 if ( (sumused>3) || (sumused>0.5*nfound) || (nfound<0.6*nfoundable)) break;
3735 polytrack.UpdateParameters();
3738 if ( (sumused>3) || (sumused>0.5*nfound)) {
3739 //printf("sumused %d\n",sumused);
3744 polytrack.GetFitDerivation(kr0.GetX(),dy,dz);
3745 AliTPCpolyTrack track2;
3747 polytrack.Refit(track2,0.5+TMath::Abs(dy)*0.3,0.4+TMath::Abs(dz)*0.3);
3748 if (track2.GetN()<0.5*nfoundable) continue;
3751 if ((nfound>0.6*nfoundable) &&( nfoundable>0.4*(i1-i2))) {
3753 // test seed with and without constrain
3754 for (Int_t constrain=0; constrain<=0;constrain++){
3755 // add polytrack candidate
3757 Double_t x[5], c[15];
3758 Double_t x1,x2,x3,y1,y2,y3,z1,z2,z3;
3759 track2.GetBoundaries(x3,x1);
3761 track2.GetFitPoint(x1,y1,z1);
3762 track2.GetFitPoint(x2,y2,z2);
3763 track2.GetFitPoint(x3,y3,z3);
3765 //is track pointing to the vertex ?
3768 polytrack.GetFitPoint(x0,y0,z0);
3781 x[4]=F1(x1,y1,x2,y2,x3,y3);
3783 // if (TMath::Abs(x[4]) >= cuts[0]) continue; //
3784 x[2]=F2(x1,y1,x2,y2,x3,y3);
3786 //if (TMath::Abs(x[4]*x1-x[2]) >= cuts[1]) continue;
3787 //x[3]=F3(x1,y1,x2,y2,z1,z2);
3788 x[3]=F3n(x1,y1,x3,y3,z1,z3,x[4]);
3789 //if (TMath::Abs(x[3]) > cuts[3]) continue;
3792 Double_t sy =0.1, sz =0.1;
3793 Double_t sy1=0.02, sz1=0.02;
3794 Double_t sy2=0.02, sz2=0.02;
3798 sy3=25000*x[4]*x[4]+0.1, sy=0.1, sz=0.1;
3801 Double_t f40=(F1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
3802 Double_t f42=(F1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
3803 Double_t f43=(F1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
3804 Double_t f20=(F2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
3805 Double_t f22=(F2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
3806 Double_t f23=(F2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
3808 Double_t f30=(F3(x1,y1+sy,x3,y3,z1,z3)-x[3])/sy;
3809 Double_t f31=(F3(x1,y1,x3,y3,z1+sz,z3)-x[3])/sz;
3810 Double_t f32=(F3(x1,y1,x3,y3+sy,z1,z3)-x[3])/sy;
3811 Double_t f34=(F3(x1,y1,x3,y3,z1,z3+sz)-x[3])/sz;
3816 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3817 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3818 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3819 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3820 c[13]=f30*sy1*f40+f32*sy2*f42;
3821 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3823 //Int_t row1 = fSectors->GetRowNumber(x1);
3824 Int_t row1 = GetRowNumber(x1);
3828 if (seed) {MarkSeedFree( seed ); seed = 0;}
3829 AliTPCseed *track = seed = new( NextFreeSeed() ) AliTPCseed(x1,sec*alpha+shift,x,c,index);
3830 seed->SetPoolID(fLastSeedID);
3831 track->SetIsSeeding(kTRUE);
3832 Int_t rc=FollowProlongation(*track, i2);
3833 if (constrain) track->SetBConstrain(1);
3835 track->SetBConstrain(0);
3836 track->SetLastPoint(row1+fInnerSec->GetNRows()); // first cluster in track position
3837 track->SetFirstPoint(track->GetLastPoint());
3839 if (rc==0 || track->GetNumberOfClusters()<(i1-i2)*0.5 ||
3840 track->GetNumberOfClusters() < track->GetNFoundable()*0.6 ||
3841 track->GetNShared()>0.4*track->GetNumberOfClusters()) {
3842 MarkSeedFree( seed ); seed = 0;
3845 arr->AddLast(track); // track IS seed, don't free seed
3846 seed = new( NextFreeSeed() ) AliTPCseed;
3847 seed->SetPoolID(fLastSeedID);
3851 } // if accepted seed
3854 Info("MakeSeeds2","\nSeeding statiistic:\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin3);
3856 if (seed) MarkSeedFree( seed );
3860 AliTPCseed *AliTPCtrackerMI::MakeSeed(AliTPCseed *const track, Float_t r0, Float_t r1, Float_t r2)
3864 //reseed using track points
3865 Int_t p0 = int(r0*track->GetNumberOfClusters()); // point 0
3866 Int_t p1 = int(r1*track->GetNumberOfClusters());
3867 Int_t p2 = int(r2*track->GetNumberOfClusters()); // last point
3869 Double_t x0[3],x1[3],x2[3];
3870 for (Int_t i=0;i<3;i++){
3876 // find track position at given ratio of the length
3877 Int_t sec0=0, sec1=0, sec2=0;
3880 for (Int_t i=0;i<160;i++){
3881 if (track->GetClusterPointer(i)){
3883 AliTPCTrackerPoint *trpoint =track->GetTrackPoint(i);
3884 if ( (index<p0) || x0[0]<0 ){
3885 if (trpoint->GetX()>1){
3886 clindex = track->GetClusterIndex2(i);
3888 x0[0] = trpoint->GetX();
3889 x0[1] = trpoint->GetY();
3890 x0[2] = trpoint->GetZ();
3891 sec0 = ((clindex&0xff000000)>>24)%18;
3896 if ( (index<p1) &&(trpoint->GetX()>1)){
3897 clindex = track->GetClusterIndex2(i);
3899 x1[0] = trpoint->GetX();
3900 x1[1] = trpoint->GetY();
3901 x1[2] = trpoint->GetZ();
3902 sec1 = ((clindex&0xff000000)>>24)%18;
3905 if ( (index<p2) &&(trpoint->GetX()>1)){
3906 clindex = track->GetClusterIndex2(i);
3908 x2[0] = trpoint->GetX();
3909 x2[1] = trpoint->GetY();
3910 x2[2] = trpoint->GetZ();
3911 sec2 = ((clindex&0xff000000)>>24)%18;
3918 Double_t alpha, cs,sn, xx2,yy2;
3920 alpha = (sec1-sec2)*fSectors->GetAlpha();
3921 cs = TMath::Cos(alpha);
3922 sn = TMath::Sin(alpha);
3923 xx2= x1[0]*cs-x1[1]*sn;
3924 yy2= x1[0]*sn+x1[1]*cs;
3928 alpha = (sec0-sec2)*fSectors->GetAlpha();
3929 cs = TMath::Cos(alpha);
3930 sn = TMath::Sin(alpha);
3931 xx2= x0[0]*cs-x0[1]*sn;
3932 yy2= x0[0]*sn+x0[1]*cs;
3938 Double_t x[5],c[15];
3942 x[4]=F1(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]);
3943 // if (x[4]>1) return 0;
3944 x[2]=F2(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]);
3945 x[3]=F3n(x2[0],x2[1],x0[0],x0[1],x2[2],x0[2],x[4]);
3946 //if (TMath::Abs(x[3]) > 2.2) return 0;
3947 //if (TMath::Abs(x[2]) > 1.99) return 0;
3949 Double_t sy =0.1, sz =0.1;
3951 Double_t sy1=0.02+track->GetSigmaY2(), sz1=0.02+track->GetSigmaZ2();
3952 Double_t sy2=0.01+track->GetSigmaY2(), sz2=0.01+track->GetSigmaZ2();
3953 Double_t sy3=0.01+track->GetSigmaY2();
3955 Double_t f40=(F1(x2[0],x2[1]+sy,x1[0],x1[1],x0[0],x0[1])-x[4])/sy;
3956 Double_t f42=(F1(x2[0],x2[1],x1[0],x1[1]+sy,x0[0],x0[1])-x[4])/sy;
3957 Double_t f43=(F1(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]+sy)-x[4])/sy;
3958 Double_t f20=(F2(x2[0],x2[1]+sy,x1[0],x1[1],x0[0],x0[1])-x[2])/sy;
3959 Double_t f22=(F2(x2[0],x2[1],x1[0],x1[1]+sy,x0[0],x0[1])-x[2])/sy;
3960 Double_t f23=(F2(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]+sy)-x[2])/sy;
3962 Double_t f30=(F3(x2[0],x2[1]+sy,x0[0],x0[1],x2[2],x0[2])-x[3])/sy;
3963 Double_t f31=(F3(x2[0],x2[1],x0[0],x0[1],x2[2]+sz,x0[2])-x[3])/sz;
3964 Double_t f32=(F3(x2[0],x2[1],x0[0],x0[1]+sy,x2[2],x0[2])-x[3])/sy;
3965 Double_t f34=(F3(x2[0],x2[1],x0[0],x0[1],x2[2],x0[2]+sz)-x[3])/sz;
3970 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3971 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3972 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3973 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3974 c[13]=f30*sy1*f40+f32*sy2*f42;
3975 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3977 // Int_t row1 = fSectors->GetRowNumber(x2[0]);
3978 AliTPCseed *seed = new( NextFreeSeed() ) AliTPCseed(x2[0], sec2*fSectors->GetAlpha()+fSectors->GetAlphaShift(), x, c, 0);
3979 seed->SetPoolID(fLastSeedID);
3980 // Double_t y0,z0,y1,z1, y2,z2;
3981 //seed->GetProlongation(x0[0],y0,z0);
3982 // seed->GetProlongation(x1[0],y1,z1);
3983 //seed->GetProlongation(x2[0],y2,z2);
3985 seed->SetLastPoint(pp2);
3986 seed->SetFirstPoint(pp2);
3993 AliTPCseed *AliTPCtrackerMI::ReSeed(const AliTPCseed *track, Float_t r0, Float_t r1, Float_t r2)
3997 //reseed using founded clusters
3999 // Find the number of clusters
4000 Int_t nclusters = 0;
4001 for (Int_t irow=0;irow<160;irow++){
4002 if (track->GetClusterIndex(irow)>0) nclusters++;
4006 ipos[0] = TMath::Max(int(r0*nclusters),0); // point 0 cluster
4007 ipos[1] = TMath::Min(int(r1*nclusters),nclusters-1); //
4008 ipos[2] = TMath::Min(int(r2*nclusters),nclusters-1); // last point
4011 Double_t xyz[3][3]={{0}};
4012 Int_t row[3]={0},sec[3]={0,0,0};
4014 // find track row position at given ratio of the length
4016 for (Int_t irow=0;irow<160;irow++){
4017 if (track->GetClusterIndex2(irow)<0) continue;
4019 for (Int_t ipoint=0;ipoint<3;ipoint++){
4020 if (index<=ipos[ipoint]) row[ipoint] = irow;
4024 //Get cluster and sector position
4025 for (Int_t ipoint=0;ipoint<3;ipoint++){
4026 Int_t clindex = track->GetClusterIndex2(row[ipoint]);
4027 AliTPCclusterMI * cl = GetClusterMI(clindex);
4030 // AliTPCclusterMI * cl = GetClusterMI(clindex);
4033 sec[ipoint] = ((clindex&0xff000000)>>24)%18;
4034 xyz[ipoint][0] = GetXrow(row[ipoint]);
4035 xyz[ipoint][1] = cl->GetY();
4036 xyz[ipoint][2] = cl->GetZ();
4040 // Calculate seed state vector and covariance matrix
4042 Double_t alpha, cs,sn, xx2,yy2;
4044 alpha = (sec[1]-sec[2])*fSectors->GetAlpha();
4045 cs = TMath::Cos(alpha);
4046 sn = TMath::Sin(alpha);
4047 xx2= xyz[1][0]*cs-xyz[1][1]*sn;
4048 yy2= xyz[1][0]*sn+xyz[1][1]*cs;
4052 alpha = (sec[0]-sec[2])*fSectors->GetAlpha();
4053 cs = TMath::Cos(alpha);
4054 sn = TMath::Sin(alpha);
4055 xx2= xyz[0][0]*cs-xyz[0][1]*sn;
4056 yy2= xyz[0][0]*sn+xyz[0][1]*cs;
4062 Double_t x[5],c[15];
4066 x[4]=F1(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
4067 x[2]=F2(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
4068 x[3]=F3n(xyz[2][0],xyz[2][1],xyz[0][0],xyz[0][1],xyz[2][2],xyz[0][2],x[4]);
4070 Double_t sy =0.1, sz =0.1;
4072 Double_t sy1=0.2, sz1=0.2;
4073 Double_t sy2=0.2, sz2=0.2;
4076 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;
4077 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;
4078 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;
4079 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;
4080 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;
4081 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;
4083 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;
4084 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;
4085 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;
4086 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;
4091 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
4092 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
4093 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
4094 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
4095 c[13]=f30*sy1*f40+f32*sy2*f42;
4096 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
4098 // Int_t row1 = fSectors->GetRowNumber(xyz[2][0]);
4099 AliTPCseed *seed=new( NextFreeSeed() ) AliTPCseed(xyz[2][0], sec[2]*fSectors->GetAlpha()+fSectors->GetAlphaShift(), x, c, 0);
4100 seed->SetPoolID(fLastSeedID);
4101 seed->SetLastPoint(row[2]);
4102 seed->SetFirstPoint(row[2]);
4107 AliTPCseed *AliTPCtrackerMI::ReSeed(AliTPCseed *track,Int_t r0, Bool_t forward)
4111 //reseed using founded clusters
4114 Int_t row[3]={0,0,0};
4115 Int_t sec[3]={0,0,0};
4117 // forward direction
4119 for (Int_t irow=r0;irow<160;irow++){
4120 if (track->GetClusterIndex(irow)>0){
4125 for (Int_t irow=160;irow>r0;irow--){
4126 if (track->GetClusterIndex(irow)>0){
4131 for (Int_t irow=row[2]-15;irow>row[0];irow--){
4132 if (track->GetClusterIndex(irow)>0){
4140 for (Int_t irow=0;irow<r0;irow++){
4141 if (track->GetClusterIndex(irow)>0){
4146 for (Int_t irow=r0;irow>0;irow--){
4147 if (track->GetClusterIndex(irow)>0){
4152 for (Int_t irow=row[2]-15;irow>row[0];irow--){
4153 if (track->GetClusterIndex(irow)>0){
4160 if ((row[2]-row[0])<20) return 0;
4161 if (row[1]==0) return 0;
4164 //Get cluster and sector position
4165 for (Int_t ipoint=0;ipoint<3;ipoint++){
4166 Int_t clindex = track->GetClusterIndex2(row[ipoint]);
4167 AliTPCclusterMI * cl = GetClusterMI(clindex);
4170 // AliTPCclusterMI * cl = GetClusterMI(clindex);
4173 sec[ipoint] = ((clindex&0xff000000)>>24)%18;
4174 xyz[ipoint][0] = GetXrow(row[ipoint]);
4175 AliTPCTrackerPoint * point = track->GetTrackPoint(row[ipoint]);
4176 if (point&&ipoint<2){
4178 xyz[ipoint][1] = point->GetY();
4179 xyz[ipoint][2] = point->GetZ();
4182 xyz[ipoint][1] = cl->GetY();
4183 xyz[ipoint][2] = cl->GetZ();
4190 // Calculate seed state vector and covariance matrix
4192 Double_t alpha, cs,sn, xx2,yy2;
4194 alpha = (sec[1]-sec[2])*fSectors->GetAlpha();
4195 cs = TMath::Cos(alpha);
4196 sn = TMath::Sin(alpha);
4197 xx2= xyz[1][0]*cs-xyz[1][1]*sn;
4198 yy2= xyz[1][0]*sn+xyz[1][1]*cs;
4202 alpha = (sec[0]-sec[2])*fSectors->GetAlpha();
4203 cs = TMath::Cos(alpha);
4204 sn = TMath::Sin(alpha);
4205 xx2= xyz[0][0]*cs-xyz[0][1]*sn;
4206 yy2= xyz[0][0]*sn+xyz[0][1]*cs;
4212 Double_t x[5],c[15];
4216 x[4]=F1(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
4217 x[2]=F2(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
4218 x[3]=F3n(xyz[2][0],xyz[2][1],xyz[0][0],xyz[0][1],xyz[2][2],xyz[0][2],x[4]);
4220 Double_t sy =0.1, sz =0.1;
4222 Double_t sy1=0.2, sz1=0.2;
4223 Double_t sy2=0.2, sz2=0.2;
4226 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;
4227 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;
4228 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;
4229 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;
4230 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;
4231 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;
4233 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;
4234 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;
4235 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;
4236 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;
4241 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
4242 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
4243 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
4244 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
4245 c[13]=f30*sy1*f40+f32*sy2*f42;
4246 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
4248 // Int_t row1 = fSectors->GetRowNumber(xyz[2][0]);
4249 AliTPCseed *seed=new( NextFreeSeed() ) AliTPCseed(xyz[2][0], sec[2]*fSectors->GetAlpha()+fSectors->GetAlphaShift(), x, c, 0);
4250 seed->SetPoolID(fLastSeedID);
4251 seed->SetLastPoint(row[2]);
4252 seed->SetFirstPoint(row[2]);
4253 for (Int_t i=row[0];i<row[2];i++){
4254 seed->SetClusterIndex(i, track->GetClusterIndex(i));
4262 void AliTPCtrackerMI::FindMultiMC(const TObjArray * array, AliESDEvent */*esd*/, Int_t iter)
4265 // find multi tracks - THIS FUNCTION IS ONLY FOR DEBUG PURPOSES
4267 // Use AliTPCReconstructor::StreamLevel()>2 if you want to tune parameters - cuts
4269 // Two reasons to have multiple find tracks
4270 // 1. Curling tracks can be find more than once
4271 // 2. Splitted tracks
4272 // a.) Multiple seeding to increase tracking efficiency - (~ 100% reached)
4273 // b.) Edge effect on the sector boundaries
4276 // Algorithm done in 2 phases - because of CPU consumption
4277 // it is n^2 algorithm - for lead-lead 20000x20000 combination are investigated
4279 // Algorihm for curling tracks sign:
4280 // 1 phase -makes a very rough fast cuts to minimize combinatorics
4281 // a.) opposite sign
4282 // b.) one of the tracks - not pointing to the primary vertex -
4283 // c.) delta tan(theta)
4285 // 2 phase - calculates DCA between tracks - time consument
4290 // General cuts - for splitted tracks and for curling tracks
4292 const Float_t kMaxdPhi = 0.2; // maximal distance in phi
4294 // Curling tracks cuts
4299 Int_t nentries = array->GetEntriesFast();
4300 AliHelix *helixes = new AliHelix[nentries];
4301 Float_t *xm = new Float_t[nentries];
4302 Float_t *dz0 = new Float_t[nentries];
4303 Float_t *dz1 = new Float_t[nentries];
4309 // Find track COG in x direction - point with best defined parameters
4311 for (Int_t i=0;i<nentries;i++){
4312 AliTPCseed* track = (AliTPCseed*)array->At(i);
4313 if (!track) continue;
4314 track->SetCircular(0);
4315 new (&helixes[i]) AliHelix(*track);
4319 track->GetDZ(GetX(),GetY(),GetZ(),GetBz(),dz);
4322 for (Int_t icl=0; icl<160; icl++){
4323 AliTPCclusterMI * cl = track->GetClusterPointer(icl);
4329 if (ncl>0) xm[i]/=Float_t(ncl);
4332 for (Int_t i0=0;i0<nentries;i0++){
4333 AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
4334 if (!track0) continue;
4335 Float_t xc0 = helixes[i0].GetHelix(6);
4336 Float_t yc0 = helixes[i0].GetHelix(7);
4337 Float_t r0 = helixes[i0].GetHelix(8);
4338 Float_t rc0 = TMath::Sqrt(xc0*xc0+yc0*yc0);
4339 Float_t fi0 = TMath::ATan2(yc0,xc0);
4341 for (Int_t i1=i0+1;i1<nentries;i1++){
4342 AliTPCseed * track1 = (AliTPCseed*)array->At(i1);
4343 if (!track1) continue;
4344 Int_t lab0=track0->GetLabel();
4345 Int_t lab1=track1->GetLabel();
4346 if (TMath::Abs(lab0)!=TMath::Abs(lab1)) continue;
4348 Float_t xc1 = helixes[i1].GetHelix(6);
4349 Float_t yc1 = helixes[i1].GetHelix(7);
4350 Float_t r1 = helixes[i1].GetHelix(8);
4351 Float_t rc1 = TMath::Sqrt(xc1*xc1+yc1*yc1);
4352 Float_t fi1 = TMath::ATan2(yc1,xc1);
4354 Float_t dfi = fi0-fi1;
4357 if (dfi>1.5*TMath::Pi()) dfi-=TMath::Pi(); // take care about edge effect
4358 if (dfi<-1.5*TMath::Pi()) dfi+=TMath::Pi(); //
4359 if (TMath::Abs(dfi)>kMaxdPhi&&helixes[i0].GetHelix(4)*helixes[i1].GetHelix(4)<0){
4361 // if short tracks with undefined sign
4362 fi1 = -TMath::ATan2(yc1,-xc1);
4365 Float_t dtheta = TMath::Abs(track0->GetTgl()-track1->GetTgl())<TMath::Abs(track0->GetTgl()+track1->GetTgl())? track0->GetTgl()-track1->GetTgl():track0->GetTgl()+track1->GetTgl();
4368 // debug stream to tune "fast cuts"
4370 Double_t dist[3]; // distance at X
4371 Double_t mdist[3]={0,0,0}; // mean distance X+-40cm
4372 track0->GetDistance(track1,0.5*(xm[i0]+xm[i1])-40.,dist,AliTracker::GetBz());
4373 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[i]);
4374 track0->GetDistance(track1,0.5*(xm[i0]+xm[i1])+40.,dist,AliTracker::GetBz());
4375 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[i]);
4376 track0->GetDistance(track1,0.5*(xm[i0]+xm[i1]),dist,AliTracker::GetBz());
4377 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[i]);
4378 for (Int_t i=0;i<3;i++) mdist[i]*=0.33333;
4382 for (Int_t icl=0; icl<160; icl++){
4383 AliTPCclusterMI * cl0 = track0->GetClusterPointer(icl);
4384 AliTPCclusterMI * cl1 = track1->GetClusterPointer(icl);
4387 if (cl0==cl1) sums++;
4391 if (AliTPCReconstructor::StreamLevel()>5) {
4392 TTreeSRedirector &cstream = *fDebugStreamer;
4397 "Tr0.="<<track0<< // seed0
4398 "Tr1.="<<track1<< // seed1
4399 "h0.="<<&helixes[i0]<<
4400 "h1.="<<&helixes[i1]<<
4402 "sum="<<sum<< //the sum of rows with cl in both
4403 "sums="<<sums<< //the sum of shared clusters
4404 "xm0="<<xm[i0]<< // the center of track
4405 "xm1="<<xm[i1]<< // the x center of track
4406 // General cut variables
4407 "dfi="<<dfi<< // distance in fi angle
4408 "dtheta="<<dtheta<< // distance int theta angle
4414 "dist0="<<dist[0]<< //distance x
4415 "dist1="<<dist[1]<< //distance y
4416 "dist2="<<dist[2]<< //distance z
4417 "mdist0="<<mdist[0]<< //distance x
4418 "mdist1="<<mdist[1]<< //distance y
4419 "mdist2="<<mdist[2]<< //distance z
4435 if (AliTPCReconstructor::StreamLevel()>1) {
4436 AliInfo("Time for curling tracks removal DEBUGGING MC");
4443 void AliTPCtrackerMI::FindSplitted(TObjArray * array, AliESDEvent */*esd*/, Int_t /*iter*/){
4445 // Find Splitted tracks and remove the one with worst quality
4446 // Corresponding debug streamer to tune selections - "Splitted2"
4448 // 0. Sort tracks according quility
4449 // 1. Propagate the tracks to the reference radius
4450 // 2. Double_t loop to select close tracks (only to speed up process)
4451 // 3. Calculate cluster overlap ratio - and remove the track if bigger than a threshold
4452 // 4. Delete temporary parameters
4454 const Double_t xref=GetXrow(63); // reference radius -IROC/OROC boundary
4456 const Double_t kCutP1=10; // delta Z cut 10 cm
4457 const Double_t kCutP2=0.15; // delta snp(fi) cut 0.15
4458 const Double_t kCutP3=0.15; // delta tgl(theta) cut 0.15
4459 const Double_t kCutAlpha=0.15; // delta alpha cut
4460 Int_t firstpoint = 0;
4461 Int_t lastpoint = 160;
4463 Int_t nentries = array->GetEntriesFast();
4464 AliExternalTrackParam *params = new AliExternalTrackParam[nentries];
4470 //0. Sort tracks according quality
4471 //1. Propagate the ext. param to reference radius
4472 Int_t nseed = array->GetEntriesFast();
4473 if (nseed<=0) return;
4474 Float_t * quality = new Float_t[nseed];
4475 Int_t * indexes = new Int_t[nseed];
4476 for (Int_t i=0; i<nseed; i++) {
4477 AliTPCseed *pt=(AliTPCseed*)array->UncheckedAt(i);
4482 pt->UpdatePoints(); //select first last max dens points
4483 Float_t * points = pt->GetPoints();
4484 if (points[3]<0.8) quality[i] =-1;
4485 quality[i] = (points[2]-points[0])+pt->GetNumberOfClusters();
4486 //prefer high momenta tracks if overlaps
4487 quality[i] *= TMath::Sqrt(TMath::Abs(pt->Pt())+0.5);
4489 AliTracker::PropagateTrackToBxByBz(&(params[i]),xref,pt->GetMass(),5.,kTRUE);
4490 AliTracker::PropagateTrackToBxByBz(&(params[i]),xref,pt->GetMass(),1.,kTRUE);
4492 TMath::Sort(nseed,quality,indexes);
4494 // 3. Loop over pair of tracks
4496 for (Int_t i0=0; i0<nseed; i0++) {
4497 Int_t index0=indexes[i0];
4498 if (!(array->UncheckedAt(index0))) continue;
4499 AliTPCseed *s1 = (AliTPCseed*)array->UncheckedAt(index0);
4500 if (!s1->IsActive()) continue;
4501 AliExternalTrackParam &par0=params[index0];
4502 for (Int_t i1=i0+1; i1<nseed; i1++) {
4503 Int_t index1=indexes[i1];
4504 if (!(array->UncheckedAt(index1))) continue;
4505 AliTPCseed *s2 = (AliTPCseed*)array->UncheckedAt(index1);
4506 if (!s2->IsActive()) continue;
4507 if (s2->GetKinkIndexes()[0]!=0)
4508 if (s2->GetKinkIndexes()[0] == -s1->GetKinkIndexes()[0]) continue;
4509 AliExternalTrackParam &par1=params[index1];
4510 if (TMath::Abs(par0.GetParameter()[3]-par1.GetParameter()[3])>kCutP3) continue;
4511 if (TMath::Abs(par0.GetParameter()[1]-par1.GetParameter()[1])>kCutP1) continue;
4512 if (TMath::Abs(par0.GetParameter()[2]-par1.GetParameter()[2])>kCutP2) continue;
4513 Double_t dAlpha= TMath::Abs(par0.GetAlpha()-par1.GetAlpha());
4514 if (dAlpha>TMath::Pi()) dAlpha-=TMath::Pi();
4515 if (TMath::Abs(dAlpha)>kCutAlpha) continue;
4520 Int_t firstShared=lastpoint, lastShared=firstpoint;
4521 Int_t firstRow=lastpoint, lastRow=firstpoint;
4523 for (Int_t i=firstpoint;i<lastpoint;i++){
4524 if (s1->GetClusterIndex2(i)>0) nall0++;
4525 if (s2->GetClusterIndex2(i)>0) nall1++;
4526 if (s1->GetClusterIndex2(i)>0 && s2->GetClusterIndex2(i)>0) {
4527 if (i<firstRow) firstRow=i;
4528 if (i>lastRow) lastRow=i;
4530 if ( (s1->GetClusterIndex2(i))==(s2->GetClusterIndex2(i)) && s1->GetClusterIndex2(i)>0) {
4531 if (i<firstShared) firstShared=i;
4532 if (i>lastShared) lastShared=i;
4536 Double_t ratio0 = Float_t(sumShared)/Float_t(TMath::Min(nall0+1,nall1+1));
4537 Double_t ratio1 = Float_t(sumShared)/Float_t(TMath::Max(nall0+1,nall1+1));
4539 if( AliTPCReconstructor::StreamLevel()>1){
4540 TTreeSRedirector &cstream = *fDebugStreamer;
4541 Int_t n0=s1->GetNumberOfClusters();
4542 Int_t n1=s2->GetNumberOfClusters();
4543 Int_t n0F=s1->GetNFoundable();
4544 Int_t n1F=s2->GetNFoundable();
4545 Int_t lab0=s1->GetLabel();
4546 Int_t lab1=s2->GetLabel();
4548 cstream<<"Splitted2"<<
4549 "iter="<<fIteration<<
4550 "lab0="<<lab0<< // MC label if exist
4551 "lab1="<<lab1<< // MC label if exist
4554 "ratio0="<<ratio0<< // shared ratio
4555 "ratio1="<<ratio1<< // shared ratio
4556 "p0.="<<&par0<< // track parameters
4558 "s0.="<<s1<< // full seed
4560 "n0="<<n0<< // number of clusters track 0
4561 "n1="<<n1<< // number of clusters track 1
4562 "nall0="<<nall0<< // number of clusters track 0
4563 "nall1="<<nall1<< // number of clusters track 1
4564 "n0F="<<n0F<< // number of findable
4565 "n1F="<<n1F<< // number of findable
4566 "shared="<<sumShared<< // number of shared clusters
4567 "firstS="<<firstShared<< // first and the last shared row
4568 "lastS="<<lastShared<<
4569 "firstRow="<<firstRow<< // first and the last row with cluster
4570 "lastRow="<<lastRow<< //
4574 // remove track with lower quality
4576 if (ratio0>AliTPCReconstructor::GetRecoParam()->GetCutSharedClusters(0) ||
4577 ratio1>AliTPCReconstructor::GetRecoParam()->GetCutSharedClusters(1)){
4581 MarkSeedFree( array->RemoveAt(index1) );
4586 // 4. Delete temporary array
4596 void AliTPCtrackerMI::FindCurling(const TObjArray * array, AliESDEvent */*esd*/, Int_t iter)
4599 // find Curling tracks
4600 // Use AliTPCReconstructor::StreamLevel()>1 if you want to tune parameters - cuts
4603 // Algorithm done in 2 phases - because of CPU consumption
4604 // it is n^2 algorithm - for lead-lead 20000x20000 combination are investigated
4605 // see detal in MC part what can be used to cut
4609 const Float_t kMaxC = 400; // maximal curvature to of the track
4610 const Float_t kMaxdTheta = 0.15; // maximal distance in theta
4611 const Float_t kMaxdPhi = 0.15; // maximal distance in phi
4612 const Float_t kPtRatio = 0.3; // ratio between pt
4613 const Float_t kMinDCAR = 2.; // distance to the primary vertex in r - see cpipe cut
4616 // Curling tracks cuts
4619 const Float_t kMaxDeltaRMax = 40; // distance in outer radius
4620 const Float_t kMaxDeltaRMin = 5.; // distance in lower radius - see cpipe cut
4621 const Float_t kMinAngle = 2.9; // angle between tracks
4622 const Float_t kMaxDist = 5; // biggest distance
4624 // The cuts can be tuned using the "MC information stored in Multi tree ==> see FindMultiMC
4627 TCut csign("csign","Tr0.fP[4]*Tr1.fP[4]<0"); //opposite sign
4628 TCut cmax("cmax","abs(Tr0.GetC())>1/400");
4629 TCut cda("cda","sqrt(dtheta^2+dfi^2)<0.15");
4630 TCut ccratio("ccratio","abs((Tr0.fP[4]+Tr1.fP[4])/(abs(Tr0.fP[4])+abs(Tr1.fP[4])))<0.3");
4631 TCut cpipe("cpipe", "min(abs(r0-rc0),abs(r1-rc1))>5");
4633 TCut cdrmax("cdrmax","abs(abs(rc0+r0)-abs(rc1+r1))<40")
4634 TCut cdrmin("cdrmin","abs(abs(rc0+r0)-abs(rc1+r1))<10")
4636 Multi->Draw("dfi","iter==0"+csign+cmax+cda+ccratio); ~94% of curling tracks fulfill
4637 Multi->Draw("min(abs(r0-rc0),abs(r1-rc1))","iter==0&&abs(lab1)==abs(lab0)"+csign+cmax+cda+ccratio+cpipe+cdrmin+cdrmax); //80%
4639 Curling2->Draw("dfi","iter==0&&abs(lab0)==abs(lab1)"+csign+cmax+cdtheta+cdfi+ccratio)
4645 Int_t nentries = array->GetEntriesFast();
4646 AliHelix *helixes = new AliHelix[nentries];
4647 for (Int_t i=0;i<nentries;i++){
4648 AliTPCseed* track = (AliTPCseed*)array->At(i);
4649 if (!track) continue;
4650 track->SetCircular(0);
4651 new (&helixes[i]) AliHelix(*track);
4657 Double_t phase[2][2]={{0,0},{0,0}},radius[2]={0,0};
4663 for (Int_t i0=0;i0<nentries;i0++){
4664 AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
4665 if (!track0) continue;
4666 if (TMath::Abs(track0->GetC())<1/kMaxC) continue;
4667 Float_t xc0 = helixes[i0].GetHelix(6);
4668 Float_t yc0 = helixes[i0].GetHelix(7);
4669 Float_t r0 = helixes[i0].GetHelix(8);
4670 Float_t rc0 = TMath::Sqrt(xc0*xc0+yc0*yc0);
4671 Float_t fi0 = TMath::ATan2(yc0,xc0);
4673 for (Int_t i1=i0+1;i1<nentries;i1++){
4674 AliTPCseed * track1 = (AliTPCseed*)array->At(i1);
4675 if (!track1) continue;
4676 if (TMath::Abs(track1->GetC())<1/kMaxC) continue;
4677 Float_t xc1 = helixes[i1].GetHelix(6);
4678 Float_t yc1 = helixes[i1].GetHelix(7);
4679 Float_t r1 = helixes[i1].GetHelix(8);
4680 Float_t rc1 = TMath::Sqrt(xc1*xc1+yc1*yc1);
4681 Float_t fi1 = TMath::ATan2(yc1,xc1);
4683 Float_t dfi = fi0-fi1;
4686 if (dfi>1.5*TMath::Pi()) dfi-=TMath::Pi(); // take care about edge effect
4687 if (dfi<-1.5*TMath::Pi()) dfi+=TMath::Pi(); //
4688 Float_t dtheta = TMath::Abs(track0->GetTgl()-track1->GetTgl())<TMath::Abs(track0->GetTgl()+track1->GetTgl())? track0->GetTgl()-track1->GetTgl():track0->GetTgl()+track1->GetTgl();
4692 if (track0->GetBConstrain()&&track1->GetBConstrain()) continue; // not constrained
4693 if (track1->GetSigned1Pt()*track0->GetSigned1Pt()>0) continue; // not the same sign
4694 if ( TMath::Abs(track1->GetTgl()+track0->GetTgl())>kMaxdTheta) continue; //distance in the Theta
4695 if ( TMath::Abs(dfi)>kMaxdPhi) continue; //distance in phi
4696 if ( TMath::Sqrt(dfi*dfi+dtheta*dtheta)>kMaxdPhi) continue; //common angular offset
4698 Float_t pt0 = track0->GetSignedPt();
4699 Float_t pt1 = track1->GetSignedPt();
4700 if ((TMath::Abs(pt0+pt1)/(TMath::Abs(pt0)+TMath::Abs(pt1)))>kPtRatio) continue;
4701 if ((iter==1) && TMath::Abs(TMath::Abs(rc0+r0)-TMath::Abs(rc1+r1))>kMaxDeltaRMax) continue;
4702 if ((iter!=1) &&TMath::Abs(TMath::Abs(rc0-r0)-TMath::Abs(rc1-r1))>kMaxDeltaRMin) continue;
4703 if (TMath::Min(TMath::Abs(rc0-r0),TMath::Abs(rc1-r1))<kMinDCAR) continue;
4706 // Now find closest approach
4710 Int_t npoints = helixes[i0].GetRPHIintersections(helixes[i1], phase, radius,10);
4711 if (npoints==0) continue;
4712 helixes[i0].GetClosestPhases(helixes[i1], phase);
4716 Double_t hangles[3];
4717 helixes[i0].Evaluate(phase[0][0],xyz0);
4718 helixes[i1].Evaluate(phase[0][1],xyz1);
4720 helixes[i0].GetAngle(phase[0][0],helixes[i1],phase[0][1],hangles);
4721 Double_t deltah[2],deltabest;
4722 if (TMath::Abs(hangles[2])<kMinAngle) continue;
4726 helixes[i0].ParabolicDCA(helixes[i1],phase[0][0],phase[0][1],radius[0],deltah[0],2);
4728 helixes[i0].ParabolicDCA(helixes[i1],phase[1][0],phase[1][1],radius[1],deltah[1],2);
4729 if (deltah[1]<deltah[0]) ibest=1;
4731 deltabest = TMath::Sqrt(deltah[ibest]);
4732 helixes[i0].Evaluate(phase[ibest][0],xyz0);
4733 helixes[i1].Evaluate(phase[ibest][1],xyz1);
4734 helixes[i0].GetAngle(phase[ibest][0],helixes[i1],phase[ibest][1],hangles);
4735 Double_t radiusbest = TMath::Sqrt(radius[ibest]);
4737 if (deltabest>kMaxDist) continue;
4738 // if (mindcar+mindcaz<40 && (TMath::Abs(hangles[2])<kMinAngle ||deltabest>3)) continue;
4739 Bool_t sign =kFALSE;
4740 if (hangles[2]>kMinAngle) sign =kTRUE;
4743 // circular[i0] = kTRUE;
4744 // circular[i1] = kTRUE;
4745 if (track0->OneOverPt()<track1->OneOverPt()){
4746 track0->SetCircular(track0->GetCircular()+1);
4747 track1->SetCircular(track1->GetCircular()+2);
4750 track1->SetCircular(track1->GetCircular()+1);
4751 track0->SetCircular(track0->GetCircular()+2);
4754 if (AliTPCReconstructor::StreamLevel()>2){
4756 //debug stream to tune "fine" cuts
4757 Int_t lab0=track0->GetLabel();
4758 Int_t lab1=track1->GetLabel();
4759 TTreeSRedirector &cstream = *fDebugStreamer;
4760 cstream<<"Curling2"<<
4776 "npoints="<<npoints<<
4777 "hangles0="<<hangles[0]<<
4778 "hangles1="<<hangles[1]<<
4779 "hangles2="<<hangles[2]<<
4782 "radius="<<radiusbest<<
4783 "deltabest="<<deltabest<<
4784 "phase0="<<phase[ibest][0]<<
4785 "phase1="<<phase[ibest][1]<<
4793 if (AliTPCReconstructor::StreamLevel()>1) {
4794 AliInfo("Time for curling tracks removal");
4800 void AliTPCtrackerMI::FindKinks(TObjArray * array, AliESDEvent *esd)
4806 // RS something is wrong in this routine: not all seeds are assigned to daughters and mothers array, but they all are queried
4809 TObjArray *kinks= new TObjArray(10000);
4810 // TObjArray *v0s= new TObjArray(10000);
4811 Int_t nentries = array->GetEntriesFast();
4812 AliHelix *helixes = new AliHelix[nentries];
4813 Int_t *sign = new Int_t[nentries];
4814 Int_t *nclusters = new Int_t[nentries];
4815 Float_t *alpha = new Float_t[nentries];
4816 AliKink *kink = new AliKink();
4817 Int_t * usage = new Int_t[nentries];
4818 Float_t *zm = new Float_t[nentries];
4819 Float_t *z0 = new Float_t[nentries];
4820 Float_t *fim = new Float_t[nentries];
4821 Float_t *shared = new Float_t[nentries];
4822 Bool_t *circular = new Bool_t[nentries];
4823 Float_t *dca = new Float_t[nentries];
4824 //const AliESDVertex * primvertex = esd->GetVertex();
4826 // nentries = array->GetEntriesFast();
4831 for (Int_t i=0;i<nentries;i++){
4834 AliTPCseed* track = (AliTPCseed*)array->At(i);
4835 if (!track) continue;
4836 track->SetCircular(0);
4838 track->UpdatePoints();
4839 if (( track->GetPoints()[2]- track->GetPoints()[0])>5 && track->GetPoints()[3]>0.8){
4841 nclusters[i]=track->GetNumberOfClusters();
4842 alpha[i] = track->GetAlpha();
4843 new (&helixes[i]) AliHelix(*track);
4845 helixes[i].Evaluate(0,xyz);
4846 sign[i] = (track->GetC()>0) ? -1:1;
4849 if (track->GetProlongation(x,y,z)){
4851 fim[i] = alpha[i]+TMath::ATan2(y,x);
4854 zm[i] = track->GetZ();
4858 circular[i]= kFALSE;
4859 if (track->GetProlongation(0,y,z)) z0[i] = z;
4860 dca[i] = track->GetD(0,0);
4866 Int_t ncandidates =0;
4869 Double_t phase[2][2]={{0,0},{0,0}},radius[2]={0,0};
4872 // Find circling track
4874 for (Int_t i0=0;i0<nentries;i0++){
4875 AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
4876 if (!track0) continue;
4877 if (track0->GetNumberOfClusters()<40) continue;
4878 if (TMath::Abs(1./track0->GetC())>200) continue;
4879 for (Int_t i1=i0+1;i1<nentries;i1++){
4880 AliTPCseed * track1 = (AliTPCseed*)array->At(i1);
4881 if (!track1) continue;
4882 if (track1->GetNumberOfClusters()<40) continue;
4883 if ( TMath::Abs(track1->GetTgl()+track0->GetTgl())>0.1) continue;
4884 if (track0->GetBConstrain()&&track1->GetBConstrain()) continue;
4885 if (TMath::Abs(1./track1->GetC())>200) continue;
4886 if (track1->GetSigned1Pt()*track0->GetSigned1Pt()>0) continue;
4887 if (track1->GetTgl()*track0->GetTgl()>0) continue;
4888 if (TMath::Max(TMath::Abs(1./track0->GetC()),TMath::Abs(1./track1->GetC()))>190) continue;
4889 if (track0->GetBConstrain()&&track1->OneOverPt()<track0->OneOverPt()) continue; //returning - lower momenta
4890 if (track1->GetBConstrain()&&track0->OneOverPt()<track1->OneOverPt()) continue; //returning - lower momenta
4892 Float_t mindcar = TMath::Min(TMath::Abs(dca[i0]),TMath::Abs(dca[i1]));
4893 if (mindcar<5) continue;
4894 Float_t mindcaz = TMath::Min(TMath::Abs(z0[i0]-GetZ()),TMath::Abs(z0[i1]-GetZ()));
4895 if (mindcaz<5) continue;
4896 if (mindcar+mindcaz<20) continue;
4899 Float_t xc0 = helixes[i0].GetHelix(6);
4900 Float_t yc0 = helixes[i0].GetHelix(7);
4901 Float_t r0 = helixes[i0].GetHelix(8);
4902 Float_t xc1 = helixes[i1].GetHelix(6);
4903 Float_t yc1 = helixes[i1].GetHelix(7);
4904 Float_t r1 = helixes[i1].GetHelix(8);
4906 Float_t rmean = (r0+r1)*0.5;
4907 Float_t delta =TMath::Sqrt((xc1-xc0)*(xc1-xc0)+(yc1-yc0)*(yc1-yc0));
4908 //if (delta>30) continue;
4909 if (delta>rmean*0.25) continue;
4910 if (TMath::Abs(r0-r1)/rmean>0.3) continue;
4912 Int_t npoints = helixes[i0].GetRPHIintersections(helixes[i1], phase, radius,10);
4913 if (npoints==0) continue;
4914 helixes[i0].GetClosestPhases(helixes[i1], phase);
4918 Double_t hangles[3];
4919 helixes[i0].Evaluate(phase[0][0],xyz0);
4920 helixes[i1].Evaluate(phase[0][1],xyz1);
4922 helixes[i0].GetAngle(phase[0][0],helixes[i1],phase[0][1],hangles);
4923 Double_t deltah[2],deltabest;
4924 if (hangles[2]<2.8) continue;
4927 helixes[i0].ParabolicDCA(helixes[i1],phase[0][0],phase[0][1],radius[0],deltah[0],2);
4929 helixes[i0].ParabolicDCA(helixes[i1],phase[1][0],phase[1][1],radius[1],deltah[1],2);
4930 if (deltah[1]<deltah[0]) ibest=1;
4932 deltabest = TMath::Sqrt(deltah[ibest]);
4933 helixes[i0].Evaluate(phase[ibest][0],xyz0);
4934 helixes[i1].Evaluate(phase[ibest][1],xyz1);
4935 helixes[i0].GetAngle(phase[ibest][0],helixes[i1],phase[ibest][1],hangles);
4936 Double_t radiusbest = TMath::Sqrt(radius[ibest]);
4938 if (deltabest>6) continue;
4939 if (mindcar+mindcaz<40 && (hangles[2]<3.12||deltabest>3)) continue;
4940 Bool_t lsign =kFALSE;
4941 if (hangles[2]>3.06) lsign =kTRUE;
4944 circular[i0] = kTRUE;
4945 circular[i1] = kTRUE;
4946 if (track0->OneOverPt()<track1->OneOverPt()){
4947 track0->SetCircular(track0->GetCircular()+1);
4948 track1->SetCircular(track1->GetCircular()+2);
4951 track1->SetCircular(track1->GetCircular()+1);
4952 track0->SetCircular(track0->GetCircular()+2);
4955 if (lsign&&AliTPCReconstructor::StreamLevel()>1){
4957 Int_t lab0=track0->GetLabel();
4958 Int_t lab1=track1->GetLabel();
4959 TTreeSRedirector &cstream = *fDebugStreamer;
4960 cstream<<"Curling"<<
4967 "mindcar="<<mindcar<<
4968 "mindcaz="<<mindcaz<<
4971 "npoints="<<npoints<<
4972 "hangles0="<<hangles[0]<<
4973 "hangles2="<<hangles[2]<<
4978 "radius="<<radiusbest<<
4979 "deltabest="<<deltabest<<
4980 "phase0="<<phase[ibest][0]<<
4981 "phase1="<<phase[ibest][1]<<
4991 for (Int_t i =0;i<nentries;i++){
4992 if (sign[i]==0) continue;
4993 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5000 Double_t cradius0 = 40*40;
5001 Double_t cradius1 = 270*270;
5004 Double_t cdist3=0.55;
5005 for (Int_t j =i+1;j<nentries;j++){
5007 if (sign[j]*sign[i]<1) continue;
5008 if ( (nclusters[i]+nclusters[j])>200) continue;
5009 if ( (nclusters[i]+nclusters[j])<80) continue;
5010 if ( TMath::Abs(zm[i]-zm[j])>60.) continue;
5011 if ( TMath::Abs(fim[i]-fim[j])>0.6 && TMath::Abs(fim[i]-fim[j])<5.7 ) continue;
5012 //AliTPCseed * track1 = (AliTPCseed*)array->At(j); Double_t phase[2][2],radius[2];
5013 Int_t npoints = helixes[i].GetRPHIintersections(helixes[j], phase, radius,20);
5014 if (npoints<1) continue;
5017 if (radius[0]<cradius0||radius[0]>cradius1) continue;
5020 if ( (radius[0]<cradius0||radius[0]>cradius1) && (radius[1]<cradius0||radius[1]>cradius1) ) continue;
5023 Double_t delta1=10000,delta2=10000;
5024 // cuts on the intersection radius
5025 helixes[i].LinearDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta1);
5026 if (radius[0]<20&&delta1<1) continue; //intersection at vertex
5027 if (radius[0]<10&&delta1<3) continue; //intersection at vertex
5029 helixes[i].LinearDCA(helixes[j],phase[1][0],phase[1][1],radius[1],delta2);
5030 if (radius[1]<20&&delta2<1) continue; //intersection at vertex
5031 if (radius[1]<10&&delta2<3) continue; //intersection at vertex
5034 Double_t distance1 = TMath::Min(delta1,delta2);
5035 if (distance1>cdist1) continue; // cut on DCA linear approximation
5037 npoints = helixes[i].GetRPHIintersections(helixes[j], phase, radius,20);
5038 helixes[i].ParabolicDCA(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
5043 helixes[i].ParabolicDCA(helixes[j],phase[1][0],phase[1][1],radius[1],delta2);
5044 if (radius[1]<20&&delta2<1) continue; //intersection at vertex
5045 if (radius[1]<10&&delta2<3) continue; //intersection at vertex
5047 distance1 = TMath::Min(delta1,delta2);
5050 rkink = TMath::Sqrt(radius[0]);
5053 rkink = TMath::Sqrt(radius[1]);
5055 if (distance1>cdist2) continue;
5058 AliTPCseed * track1 = (AliTPCseed*)array->At(j);
5061 Int_t row0 = GetRowNumber(rkink);
5062 if (row0<10) continue;
5063 if (row0>150) continue;
5066 Float_t dens00=-1,dens01=-1;
5067 Float_t dens10=-1,dens11=-1;
5069 Int_t found,foundable,ishared;
5070 track0->GetClusterStatistic(0,row0-5, found, foundable,ishared,kFALSE);
5071 if (foundable>5) dens00 = Float_t(found)/Float_t(foundable);
5072 track0->GetClusterStatistic(row0+5,155, found, foundable,ishared,kFALSE);
5073 if (foundable>5) dens01 = Float_t(found)/Float_t(foundable);
5075 track1->GetClusterStatistic(0,row0-5, found, foundable,ishared,kFALSE);
5076 if (foundable>10) dens10 = Float_t(found)/Float_t(foundable);
5077 track1->GetClusterStatistic(row0+5,155, found, foundable,ishared,kFALSE);
5078 if (foundable>10) dens11 = Float_t(found)/Float_t(foundable);
5080 if (dens00<dens10 && dens01<dens11) continue;
5081 if (dens00>dens10 && dens01>dens11) continue;
5082 if (TMath::Max(dens00,dens10)<0.1) continue;
5083 if (TMath::Max(dens01,dens11)<0.3) continue;
5085 if (TMath::Min(dens00,dens10)>0.6) continue;
5086 if (TMath::Min(dens01,dens11)>0.6) continue;
5089 AliTPCseed * ktrack0, *ktrack1;
5098 if (TMath::Abs(ktrack0->GetC())>5) continue; // cut on the curvature for mother particle
5099 AliExternalTrackParam paramm(*ktrack0);
5100 AliExternalTrackParam paramd(*ktrack1);
5101 if (row0>60&&ktrack1->GetReference().GetX()>90.)new (¶md) AliExternalTrackParam(ktrack1->GetReference());
5104 kink->SetMother(paramm);
5105 kink->SetDaughter(paramd);
5108 Float_t x[3] = { kink->GetPosition()[0],kink->GetPosition()[1],kink->GetPosition()[2]};
5110 fkParam->Transform0to1(x,index);
5111 fkParam->Transform1to2(x,index);
5112 row0 = GetRowNumber(x[0]);
5114 if (kink->GetR()<100) continue;
5115 if (kink->GetR()>240) continue;
5116 if (kink->GetPosition()[2]/kink->GetR()>AliTPCReconstructor::GetCtgRange()) continue; //out of fiducial volume
5117 if (kink->GetDistance()>cdist3) continue;
5118 Float_t dird = kink->GetDaughterP()[0]*kink->GetPosition()[0]+kink->GetDaughterP()[1]*kink->GetPosition()[1]; // rough direction estimate
5119 if (dird<0) continue;
5121 Float_t dirm = kink->GetMotherP()[0]*kink->GetPosition()[0]+kink->GetMotherP()[1]*kink->GetPosition()[1]; // rough direction estimate
5122 if (dirm<0) continue;
5123 Float_t mpt = TMath::Sqrt(kink->GetMotherP()[0]*kink->GetMotherP()[0]+kink->GetMotherP()[1]*kink->GetMotherP()[1]);
5124 if (mpt<0.2) continue;
5127 //for high momenta momentum not defined well in first iteration
5128 Double_t qt = TMath::Sin(kink->GetAngle(2))*ktrack1->GetP();
5129 if (qt>0.35) continue;
5132 kink->SetLabel(CookLabel(ktrack0,0.4,0,row0),0);
5133 kink->SetLabel(CookLabel(ktrack1,0.4,row0,160),1);
5135 kink->SetTPCDensity(dens00,0,0);
5136 kink->SetTPCDensity(dens01,0,1);
5137 kink->SetTPCDensity(dens10,1,0);
5138 kink->SetTPCDensity(dens11,1,1);
5139 kink->SetIndex(i,0);
5140 kink->SetIndex(j,1);
5143 kink->SetTPCDensity(dens10,0,0);
5144 kink->SetTPCDensity(dens11,0,1);
5145 kink->SetTPCDensity(dens00,1,0);
5146 kink->SetTPCDensity(dens01,1,1);
5147 kink->SetIndex(j,0);
5148 kink->SetIndex(i,1);
5151 if (mpt<1||kink->GetAngle(2)>0.1){
5152 // angle and densities not defined yet
5153 if (kink->GetTPCDensityFactor()<0.8) continue;
5154 if ((2-kink->GetTPCDensityFactor())*kink->GetDistance() >0.25) continue;
5155 if (kink->GetAngle(2)*ktrack0->GetP()<0.003) continue; //too small angle
5156 if (kink->GetAngle(2)>0.2&&kink->GetTPCDensityFactor()<1.15) continue;
5157 if (kink->GetAngle(2)>0.2&&kink->GetTPCDensity(0,1)>0.05) continue;
5159 Float_t criticalangle = track0->GetSigmaSnp2()+track0->GetSigmaTgl2();
5160 criticalangle+= track1->GetSigmaSnp2()+track1->GetSigmaTgl2();
5161 criticalangle= 3*TMath::Sqrt(criticalangle);
5162 if (criticalangle>0.02) criticalangle=0.02;
5163 if (kink->GetAngle(2)<criticalangle) continue;
5166 Int_t drow = Int_t(2.+0.5/(0.05+kink->GetAngle(2))); // overlap region defined
5167 Float_t shapesum =0;
5169 for ( Int_t row = row0-drow; row<row0+drow;row++){
5170 if (row<0) continue;
5171 if (row>155) continue;
5172 if (ktrack0->GetClusterPointer(row)){
5173 AliTPCTrackerPoint *point =ktrack0->GetTrackPoint(row);
5174 shapesum+=point->GetSigmaY()+point->GetSigmaZ();
5177 if (ktrack1->GetClusterPointer(row)){
5178 AliTPCTrackerPoint *point =ktrack1->GetTrackPoint(row);
5179 shapesum+=point->GetSigmaY()+point->GetSigmaZ();
5184 kink->SetShapeFactor(-1.);
5187 kink->SetShapeFactor(shapesum/sum);
5189 // esd->AddKink(kink);
5191 // kink->SetMother(paramm);
5192 //kink->SetDaughter(paramd);
5194 Double_t chi2P2 = paramm.GetParameter()[2]-paramd.GetParameter()[2];
5196 chi2P2/=paramm.GetCovariance()[5]+paramd.GetCovariance()[5];
5197 Double_t chi2P3 = paramm.GetParameter()[3]-paramd.GetParameter()[3];
5199 chi2P3/=paramm.GetCovariance()[9]+paramd.GetCovariance()[9];
5201 if (AliTPCReconstructor::StreamLevel()>1) {
5202 (*fDebugStreamer)<<"kinkLpt"<<
5210 if ( chi2P2+chi2P3<AliTPCReconstructor::GetRecoParam()->GetKinkAngleCutChi2(0)){
5214 kinks->AddLast(kink);
5220 // sort the kinks according quality - and refit them towards vertex
5222 Int_t nkinks = kinks->GetEntriesFast();
5223 Float_t *quality = new Float_t[nkinks];
5224 Int_t *indexes = new Int_t[nkinks];
5225 AliTPCseed *mothers = new AliTPCseed[nkinks];
5226 AliTPCseed *daughters = new AliTPCseed[nkinks];
5229 for (Int_t i=0;i<nkinks;i++){
5231 AliKink *kinkl = (AliKink*)kinks->At(i);
5233 // refit kinks towards vertex
5235 Int_t index0 = kinkl->GetIndex(0);
5236 Int_t index1 = kinkl->GetIndex(1);
5237 AliTPCseed * ktrack0 = (AliTPCseed*)array->At(index0);
5238 AliTPCseed * ktrack1 = (AliTPCseed*)array->At(index1);
5240 Int_t sumn=ktrack0->GetNumberOfClusters()+ktrack1->GetNumberOfClusters();
5242 // Refit Kink under if too small angle
5244 if (kinkl->GetAngle(2)<0.05){
5245 kinkl->SetTPCRow0(GetRowNumber(kinkl->GetR()));
5246 Int_t row0 = kinkl->GetTPCRow0();
5247 Int_t drow = Int_t(2.+0.5/(0.05+kinkl->GetAngle(2)));
5250 Int_t last = row0-drow;
5251 if (last<40) last=40;
5252 if (last<ktrack0->GetFirstPoint()+25) last = ktrack0->GetFirstPoint()+25;
5253 AliTPCseed* seed0 = ReSeed(ktrack0,last,kFALSE);
5256 Int_t first = row0+drow;
5257 if (first>130) first=130;
5258 if (first>ktrack1->GetLastPoint()-25) first = TMath::Max(ktrack1->GetLastPoint()-25,30);
5259 AliTPCseed* seed1 = ReSeed(ktrack1,first,kTRUE);
5261 if (seed0 && seed1){
5262 kinkl->SetStatus(1,8);
5263 if (RefitKink(*seed0,*seed1,*kinkl)) kinkl->SetStatus(1,9);
5264 row0 = GetRowNumber(kinkl->GetR());
5265 sumn = seed0->GetNumberOfClusters()+seed1->GetNumberOfClusters();
5266 mothers[i] = *seed0;
5267 daughters[i] = *seed1;
5270 delete kinks->RemoveAt(i);
5271 if (seed0) MarkSeedFree( seed0 );
5272 if (seed1) MarkSeedFree( seed1 );
5275 if (kinkl->GetDistance()>0.5 || kinkl->GetR()<110 || kinkl->GetR()>240) {
5276 delete kinks->RemoveAt(i);
5277 if (seed0) MarkSeedFree( seed0 );
5278 if (seed1) MarkSeedFree( seed1 );
5282 MarkSeedFree( seed0 );
5283 MarkSeedFree( seed1 );
5286 if (kinkl) quality[i] = 160*((0.1+kinkl->GetDistance())*(2.-kinkl->GetTPCDensityFactor()))/(sumn+40.); //the longest -clossest will win
5288 TMath::Sort(nkinks,quality,indexes,kFALSE);
5290 //remove double find kinks
5292 for (Int_t ikink0=1;ikink0<nkinks;ikink0++){
5293 AliKink * kink0 = (AliKink*) kinks->At(indexes[ikink0]);
5294 if (!kink0) continue;
5296 for (Int_t ikink1=0;ikink1<ikink0;ikink1++){
5297 kink0 = (AliKink*) kinks->At(indexes[ikink0]);
5298 if (!kink0) continue;
5299 AliKink * kink1 = (AliKink*) kinks->At(indexes[ikink1]);
5300 if (!kink1) continue;
5301 // if not close kink continue
5302 if (TMath::Abs(kink1->GetPosition()[2]-kink0->GetPosition()[2])>10) continue;
5303 if (TMath::Abs(kink1->GetPosition()[1]-kink0->GetPosition()[1])>10) continue;
5304 if (TMath::Abs(kink1->GetPosition()[0]-kink0->GetPosition()[0])>10) continue;
5306 AliTPCseed &mother0 = mothers[indexes[ikink0]];
5307 AliTPCseed &daughter0 = daughters[indexes[ikink0]];
5308 AliTPCseed &mother1 = mothers[indexes[ikink1]];
5309 AliTPCseed &daughter1 = daughters[indexes[ikink1]];
5310 Int_t row0 = (kink0->GetTPCRow0()+kink1->GetTPCRow0())/2;
5319 for (Int_t i=0;i<row0;i++){
5320 if (mother0.GetClusterIndex(i)>0 && mother1.GetClusterIndex(i)>0){
5323 if (mother0.GetClusterIndex(i)==mother1.GetClusterIndex(i)){
5330 for (Int_t i=row0;i<158;i++){
5331 //if (daughter0.GetClusterIndex(i)>0 && daughter0.GetClusterIndex(i)>0){ // RS: Bug?
5332 if (daughter0.GetClusterIndex(i)>0 && daughter1.GetClusterIndex(i)>0){
5335 if (mother0.GetClusterIndex(i)==mother1.GetClusterIndex(i)){
5341 Float_t ratio = Float_t(same+1)/Float_t(both+1);
5342 Float_t ratiom = Float_t(samem+1)/Float_t(bothm+1);
5343 Float_t ratiod = Float_t(samed+1)/Float_t(bothd+1);
5344 if (ratio>0.3 && ratiom>0.5 &&ratiod>0.5) {
5345 Int_t sum0 = mother0.GetNumberOfClusters()+daughter0.GetNumberOfClusters();
5346 Int_t sum1 = mother1.GetNumberOfClusters()+daughter1.GetNumberOfClusters();
5348 shared[kink0->GetIndex(0)]= kTRUE;
5349 shared[kink0->GetIndex(1)]= kTRUE;
5350 delete kinks->RemoveAt(indexes[ikink0]);
5354 shared[kink1->GetIndex(0)]= kTRUE;
5355 shared[kink1->GetIndex(1)]= kTRUE;
5356 delete kinks->RemoveAt(indexes[ikink1]);
5363 for (Int_t i=0;i<nkinks;i++){
5364 AliKink * kinkl = (AliKink*) kinks->At(indexes[i]);
5365 if (!kinkl) continue;
5366 kinkl->SetTPCRow0(GetRowNumber(kinkl->GetR()));
5367 Int_t index0 = kinkl->GetIndex(0);
5368 Int_t index1 = kinkl->GetIndex(1);
5369 if (circular[index0]||(circular[index1]&&kinkl->GetDistance()>0.2)) continue;
5370 kinkl->SetMultiple(usage[index0],0);
5371 kinkl->SetMultiple(usage[index1],1);
5372 if (kinkl->GetMultiple()[0]+kinkl->GetMultiple()[1]>2) continue;
5373 if (kinkl->GetMultiple()[0]+kinkl->GetMultiple()[1]>0 && quality[indexes[i]]>0.2) continue;
5374 if (kinkl->GetMultiple()[0]+kinkl->GetMultiple()[1]>0 && kinkl->GetDistance()>0.2) continue;
5375 if (circular[index0]||(circular[index1]&&kinkl->GetDistance()>0.1)) continue;
5377 AliTPCseed * ktrack0 = (AliTPCseed*)array->At(index0);
5378 AliTPCseed * ktrack1 = (AliTPCseed*)array->At(index1);
5379 if (!ktrack0 || !ktrack1) continue;
5380 Int_t index = esd->AddKink(kinkl);
5383 if ( ktrack0->GetKinkIndex(0)==0 && ktrack1->GetKinkIndex(0)==0) { //best kink
5384 if (mothers[indexes[i]].GetNumberOfClusters()>20 && daughters[indexes[i]].GetNumberOfClusters()>20 && (mothers[indexes[i]].GetNumberOfClusters()+daughters[indexes[i]].GetNumberOfClusters())>100){
5385 *ktrack0 = mothers[indexes[i]];
5386 *ktrack1 = daughters[indexes[i]];
5390 ktrack0->SetKinkIndex(usage[index0],-(index+1));
5391 ktrack1->SetKinkIndex(usage[index1], (index+1));
5396 // Remove tracks corresponding to shared kink's
5398 for (Int_t i=0;i<nentries;i++){
5399 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5400 if (!track0) continue;
5401 if (track0->GetKinkIndex(0)!=0) continue;
5402 if (shared[i]) MarkSeedFree( array->RemoveAt(i) );
5407 RemoveUsed2(array,0.5,0.4,30);
5409 for (Int_t i=0;i<nentries;i++){
5410 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5411 if (!track0) continue;
5412 track0->CookdEdx(0.02,0.6);
5416 for (Int_t i=0;i<nentries;i++){
5417 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5418 if (!track0) continue;
5419 if (track0->Pt()<1.4) continue;
5420 //remove double high momenta tracks - overlapped with kink candidates
5423 for (Int_t icl=track0->GetFirstPoint();icl<track0->GetLastPoint(); icl++){
5424 if (track0->GetClusterPointer(icl)!=0){
5426 if (track0->GetClusterPointer(icl)->IsUsed(10)) ishared++;
5429 if (Float_t(ishared+1)/Float_t(all+1)>0.5) {
5430 MarkSeedFree( array->RemoveAt(i) );
5434 if (track0->GetKinkIndex(0)!=0) continue;
5435 if (track0->GetNumberOfClusters()<80) continue;
5437 AliTPCseed *pmother = new AliTPCseed();
5438 AliTPCseed *pdaughter = new AliTPCseed();
5439 AliKink *pkink = new AliKink;
5441 AliTPCseed & mother = *pmother;
5442 AliTPCseed & daughter = *pdaughter;
5443 AliKink & kinkl = *pkink;
5444 if (CheckKinkPoint(track0,mother,daughter, kinkl)){
5445 if (mother.GetNumberOfClusters()<30||daughter.GetNumberOfClusters()<20) {
5449 continue; //too short tracks
5451 if (mother.Pt()<1.4) {
5457 Int_t row0= kinkl.GetTPCRow0();
5458 if (kinkl.GetDistance()>0.5 || kinkl.GetR()<110. || kinkl.GetR()>240.) {
5465 Int_t index = esd->AddKink(&kinkl);
5466 mother.SetKinkIndex(0,-(index+1));
5467 daughter.SetKinkIndex(0,index+1);
5468 if (mother.GetNumberOfClusters()>50) {
5469 MarkSeedFree( array->RemoveAt(i) );
5470 AliTPCseed* mtc = new( NextFreeSeed() ) AliTPCseed(mother);
5471 mtc->SetPoolID(fLastSeedID);
5472 array->AddAt(mtc,i);
5475 AliTPCseed* mtc = new( NextFreeSeed() ) AliTPCseed(mother);
5476 mtc->SetPoolID(fLastSeedID);
5477 array->AddLast(mtc);
5479 AliTPCseed* dtc = new( NextFreeSeed() ) AliTPCseed(daughter);
5480 dtc->SetPoolID(fLastSeedID);
5481 array->AddLast(dtc);
5482 for (Int_t icl=0;icl<row0;icl++) {
5483 if (mother.GetClusterPointer(icl)) mother.GetClusterPointer(icl)->Use(20);
5486 for (Int_t icl=row0;icl<158;icl++) {
5487 if (daughter.GetClusterPointer(icl)) daughter.GetClusterPointer(icl)->Use(20);
5496 delete [] daughters;
5518 AliInfo(Form("Ncandidates=\t%d\t%d\t%d\t%d\n",esd->GetNumberOfKinks(),ncandidates,ntracks,nall));
5524 void AliTPCtrackerMI::FindKinks(TObjArray * array, AliESDEvent *esd)
5531 TObjArray *kinks= new TObjArray(10000);
5532 // TObjArray *v0s= new TObjArray(10000);
5533 Int_t nentries = array->GetEntriesFast();
5534 AliHelix *helixes = new AliHelix[nentries];
5535 Int_t *sign = new Int_t[nentries];
5536 Int_t *nclusters = new Int_t[nentries];
5537 Float_t *alpha = new Float_t[nentries];
5538 AliKink *kink = new AliKink();
5539 Int_t * usage = new Int_t[nentries];
5540 Float_t *zm = new Float_t[nentries];
5541 Float_t *z0 = new Float_t[nentries];
5542 Float_t *fim = new Float_t[nentries];
5543 Float_t *shared = new Float_t[nentries];
5544 Bool_t *circular = new Bool_t[nentries];
5545 Float_t *dca = new Float_t[nentries];
5546 //const AliESDVertex * primvertex = esd->GetVertex();
5548 // nentries = array->GetEntriesFast();
5553 for (Int_t i=0;i<nentries;i++){
5556 AliTPCseed* track = (AliTPCseed*)array->At(i);
5557 if (!track) continue;
5558 track->SetCircular(0);
5560 track->UpdatePoints();
5561 if (( track->GetPoints()[2]- track->GetPoints()[0])>5 && track->GetPoints()[3]>0.8){
5563 nclusters[i]=track->GetNumberOfClusters();
5564 alpha[i] = track->GetAlpha();
5565 new (&helixes[i]) AliHelix(*track);
5567 helixes[i].Evaluate(0,xyz);
5568 sign[i] = (track->GetC()>0) ? -1:1;
5571 if (track->GetProlongation(x,y,z)){
5573 fim[i] = alpha[i]+TMath::ATan2(y,x);
5576 zm[i] = track->GetZ();
5580 circular[i]= kFALSE;
5581 if (track->GetProlongation(0,y,z)) z0[i] = z;
5582 dca[i] = track->GetD(0,0);
5588 Int_t ncandidates =0;
5591 Double_t phase[2][2]={{0,0},{0,0}},radius[2]={0,0};
5594 // Find circling track
5596 for (Int_t i0=0;i0<nentries;i0++){
5597 AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
5598 if (!track0) continue;
5599 if (track0->GetNumberOfClusters()<40) continue;
5600 if (TMath::Abs(1./track0->GetC())>200) continue;
5601 for (Int_t i1=i0+1;i1<nentries;i1++){
5602 AliTPCseed * track1 = (AliTPCseed*)array->At(i1);
5603 if (!track1) continue;
5604 if (track1->GetNumberOfClusters()<40) continue;
5605 if ( TMath::Abs(track1->GetTgl()+track0->GetTgl())>0.1) continue;
5606 if (track0->GetBConstrain()&&track1->GetBConstrain()) continue;
5607 if (TMath::Abs(1./track1->GetC())>200) continue;
5608 if (track1->GetSigned1Pt()*track0->GetSigned1Pt()>0) continue;
5609 if (track1->GetTgl()*track0->GetTgl()>0) continue;
5610 if (TMath::Max(TMath::Abs(1./track0->GetC()),TMath::Abs(1./track1->GetC()))>190) continue;
5611 if (track0->GetBConstrain()&&track1->OneOverPt()<track0->OneOverPt()) continue; //returning - lower momenta
5612 if (track1->GetBConstrain()&&track0->OneOverPt()<track1->OneOverPt()) continue; //returning - lower momenta
5614 Float_t mindcar = TMath::Min(TMath::Abs(dca[i0]),TMath::Abs(dca[i1]));
5615 if (mindcar<5) continue;
5616 Float_t mindcaz = TMath::Min(TMath::Abs(z0[i0]-GetZ()),TMath::Abs(z0[i1]-GetZ()));
5617 if (mindcaz<5) continue;
5618 if (mindcar+mindcaz<20) continue;
5621 Float_t xc0 = helixes[i0].GetHelix(6);
5622 Float_t yc0 = helixes[i0].GetHelix(7);
5623 Float_t r0 = helixes[i0].GetHelix(8);
5624 Float_t xc1 = helixes[i1].GetHelix(6);
5625 Float_t yc1 = helixes[i1].GetHelix(7);
5626 Float_t r1 = helixes[i1].GetHelix(8);
5628 Float_t rmean = (r0+r1)*0.5;
5629 Float_t delta =TMath::Sqrt((xc1-xc0)*(xc1-xc0)+(yc1-yc0)*(yc1-yc0));
5630 //if (delta>30) continue;
5631 if (delta>rmean*0.25) continue;
5632 if (TMath::Abs(r0-r1)/rmean>0.3) continue;
5634 Int_t npoints = helixes[i0].GetRPHIintersections(helixes[i1], phase, radius,10);
5635 if (npoints==0) continue;
5636 helixes[i0].GetClosestPhases(helixes[i1], phase);
5640 Double_t hangles[3];
5641 helixes[i0].Evaluate(phase[0][0],xyz0);
5642 helixes[i1].Evaluate(phase[0][1],xyz1);
5644 helixes[i0].GetAngle(phase[0][0],helixes[i1],phase[0][1],hangles);
5645 Double_t deltah[2],deltabest;
5646 if (hangles[2]<2.8) continue;
5649 helixes[i0].ParabolicDCA(helixes[i1],phase[0][0],phase[0][1],radius[0],deltah[0],2);
5651 helixes[i0].ParabolicDCA(helixes[i1],phase[1][0],phase[1][1],radius[1],deltah[1],2);
5652 if (deltah[1]<deltah[0]) ibest=1;
5654 deltabest = TMath::Sqrt(deltah[ibest]);
5655 helixes[i0].Evaluate(phase[ibest][0],xyz0);
5656 helixes[i1].Evaluate(phase[ibest][1],xyz1);
5657 helixes[i0].GetAngle(phase[ibest][0],helixes[i1],phase[ibest][1],hangles);
5658 Double_t radiusbest = TMath::Sqrt(radius[ibest]);
5660 if (deltabest>6) continue;
5661 if (mindcar+mindcaz<40 && (hangles[2]<3.12||deltabest>3)) continue;
5662 Bool_t lsign =kFALSE;
5663 if (hangles[2]>3.06) lsign =kTRUE;
5666 circular[i0] = kTRUE;
5667 circular[i1] = kTRUE;
5668 if (track0->OneOverPt()<track1->OneOverPt()){
5669 track0->SetCircular(track0->GetCircular()+1);
5670 track1->SetCircular(track1->GetCircular()+2);
5673 track1->SetCircular(track1->GetCircular()+1);
5674 track0->SetCircular(track0->GetCircular()+2);
5677 if (lsign&&AliTPCReconstructor::StreamLevel()>1){
5679 Int_t lab0=track0->GetLabel();
5680 Int_t lab1=track1->GetLabel();
5681 TTreeSRedirector &cstream = *fDebugStreamer;
5682 cstream<<"Curling"<<
5689 "mindcar="<<mindcar<<
5690 "mindcaz="<<mindcaz<<
5693 "npoints="<<npoints<<
5694 "hangles0="<<hangles[0]<<
5695 "hangles2="<<hangles[2]<<
5700 "radius="<<radiusbest<<
5701 "deltabest="<<deltabest<<
5702 "phase0="<<phase[ibest][0]<<
5703 "phase1="<<phase[ibest][1]<<
5713 for (Int_t i =0;i<nentries;i++){
5714 if (sign[i]==0) continue;
5715 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5722 Double_t cradius0 = 40*40;
5723 Double_t cradius1 = 270*270;
5726 Double_t cdist3=0.55;
5727 for (Int_t j =i+1;j<nentries;j++){
5729 if (sign[j]*sign[i]<1) continue;
5730 if ( (nclusters[i]+nclusters[j])>200) continue;
5731 if ( (nclusters[i]+nclusters[j])<80) continue;
5732 if ( TMath::Abs(zm[i]-zm[j])>60.) continue;
5733 if ( TMath::Abs(fim[i]-fim[j])>0.6 && TMath::Abs(fim[i]-fim[j])<5.7 ) continue;
5734 //AliTPCseed * track1 = (AliTPCseed*)array->At(j); Double_t phase[2][2],radius[2];
5735 Int_t npoints = helixes[i].GetRPHIintersections(helixes[j], phase, radius,20);
5736 if (npoints<1) continue;
5739 if (radius[0]<cradius0||radius[0]>cradius1) continue;
5742 if ( (radius[0]<cradius0||radius[0]>cradius1) && (radius[1]<cradius0||radius[1]>cradius1) ) continue;
5745 Double_t delta1=10000,delta2=10000;
5746 // cuts on the intersection radius
5747 helixes[i].LinearDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta1);
5748 if (radius[0]<20&&delta1<1) continue; //intersection at vertex
5749 if (radius[0]<10&&delta1<3) continue; //intersection at vertex
5751 helixes[i].LinearDCA(helixes[j],phase[1][0],phase[1][1],radius[1],delta2);
5752 if (radius[1]<20&&delta2<1) continue; //intersection at vertex
5753 if (radius[1]<10&&delta2<3) continue; //intersection at vertex
5756 Double_t distance1 = TMath::Min(delta1,delta2);
5757 if (distance1>cdist1) continue; // cut on DCA linear approximation
5759 npoints = helixes[i].GetRPHIintersections(helixes[j], phase, radius,20);
5760 helixes[i].ParabolicDCA(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
5765 helixes[i].ParabolicDCA(helixes[j],phase[1][0],phase[1][1],radius[1],delta2);
5766 if (radius[1]<20&&delta2<1) continue; //intersection at vertex
5767 if (radius[1]<10&&delta2<3) continue; //intersection at vertex
5769 distance1 = TMath::Min(delta1,delta2);
5772 rkink = TMath::Sqrt(radius[0]);
5775 rkink = TMath::Sqrt(radius[1]);
5777 if (distance1>cdist2) continue;
5780 AliTPCseed * track1 = (AliTPCseed*)array->At(j);
5783 Int_t row0 = GetRowNumber(rkink);
5784 if (row0<10) continue;
5785 if (row0>150) continue;
5788 Float_t dens00=-1,dens01=-1;
5789 Float_t dens10=-1,dens11=-1;
5791 Int_t found,foundable,ishared;
5792 track0->GetClusterStatistic(0,row0-5, found, foundable,ishared,kFALSE);
5793 if (foundable>5) dens00 = Float_t(found)/Float_t(foundable);
5794 track0->GetClusterStatistic(row0+5,155, found, foundable,ishared,kFALSE);
5795 if (foundable>5) dens01 = Float_t(found)/Float_t(foundable);
5797 track1->GetClusterStatistic(0,row0-5, found, foundable,ishared,kFALSE);
5798 if (foundable>10) dens10 = Float_t(found)/Float_t(foundable);
5799 track1->GetClusterStatistic(row0+5,155, found, foundable,ishared,kFALSE);
5800 if (foundable>10) dens11 = Float_t(found)/Float_t(foundable);
5802 if (dens00<dens10 && dens01<dens11) continue;
5803 if (dens00>dens10 && dens01>dens11) continue;
5804 if (TMath::Max(dens00,dens10)<0.1) continue;
5805 if (TMath::Max(dens01,dens11)<0.3) continue;
5807 if (TMath::Min(dens00,dens10)>0.6) continue;
5808 if (TMath::Min(dens01,dens11)>0.6) continue;
5811 AliTPCseed * ktrack0, *ktrack1;
5820 if (TMath::Abs(ktrack0->GetC())>5) continue; // cut on the curvature for mother particle
5821 AliExternalTrackParam paramm(*ktrack0);
5822 AliExternalTrackParam paramd(*ktrack1);
5823 if (row0>60&&ktrack1->GetReference().GetX()>90.)new (¶md) AliExternalTrackParam(ktrack1->GetReference());
5826 kink->SetMother(paramm);
5827 kink->SetDaughter(paramd);
5830 Float_t x[3] = { kink->GetPosition()[0],kink->GetPosition()[1],kink->GetPosition()[2]};
5832 fkParam->Transform0to1(x,index);
5833 fkParam->Transform1to2(x,index);
5834 row0 = GetRowNumber(x[0]);
5836 if (kink->GetR()<100) continue;
5837 if (kink->GetR()>240) continue;
5838 if (kink->GetPosition()[2]/kink->GetR()>AliTPCReconstructor::GetCtgRange()) continue; //out of fiducial volume
5839 if (kink->GetDistance()>cdist3) continue;
5840 Float_t dird = kink->GetDaughterP()[0]*kink->GetPosition()[0]+kink->GetDaughterP()[1]*kink->GetPosition()[1]; // rough direction estimate
5841 if (dird<0) continue;
5843 Float_t dirm = kink->GetMotherP()[0]*kink->GetPosition()[0]+kink->GetMotherP()[1]*kink->GetPosition()[1]; // rough direction estimate
5844 if (dirm<0) continue;
5845 Float_t mpt = TMath::Sqrt(kink->GetMotherP()[0]*kink->GetMotherP()[0]+kink->GetMotherP()[1]*kink->GetMotherP()[1]);
5846 if (mpt<0.2) continue;
5849 //for high momenta momentum not defined well in first iteration
5850 Double_t qt = TMath::Sin(kink->GetAngle(2))*ktrack1->GetP();
5851 if (qt>0.35) continue;
5854 kink->SetLabel(CookLabel(ktrack0,0.4,0,row0),0);
5855 kink->SetLabel(CookLabel(ktrack1,0.4,row0,160),1);
5857 kink->SetTPCDensity(dens00,0,0);
5858 kink->SetTPCDensity(dens01,0,1);
5859 kink->SetTPCDensity(dens10,1,0);
5860 kink->SetTPCDensity(dens11,1,1);
5861 kink->SetIndex(i,0);
5862 kink->SetIndex(j,1);
5865 kink->SetTPCDensity(dens10,0,0);
5866 kink->SetTPCDensity(dens11,0,1);
5867 kink->SetTPCDensity(dens00,1,0);
5868 kink->SetTPCDensity(dens01,1,1);
5869 kink->SetIndex(j,0);
5870 kink->SetIndex(i,1);
5873 if (mpt<1||kink->GetAngle(2)>0.1){
5874 // angle and densities not defined yet
5875 if (kink->GetTPCDensityFactor()<0.8) continue;
5876 if ((2-kink->GetTPCDensityFactor())*kink->GetDistance() >0.25) continue;
5877 if (kink->GetAngle(2)*ktrack0->GetP()<0.003) continue; //too small angle
5878 if (kink->GetAngle(2)>0.2&&kink->GetTPCDensityFactor()<1.15) continue;
5879 if (kink->GetAngle(2)>0.2&&kink->GetTPCDensity(0,1)>0.05) continue;
5881 Float_t criticalangle = track0->GetSigmaSnp2()+track0->GetSigmaTgl2();
5882 criticalangle+= track1->GetSigmaSnp2()+track1->GetSigmaTgl2();
5883 criticalangle= 3*TMath::Sqrt(criticalangle);
5884 if (criticalangle>0.02) criticalangle=0.02;
5885 if (kink->GetAngle(2)<criticalangle) continue;
5888 Int_t drow = Int_t(2.+0.5/(0.05+kink->GetAngle(2))); // overlap region defined
5889 Float_t shapesum =0;
5891 for ( Int_t row = row0-drow; row<row0+drow;row++){
5892 if (row<0) continue;
5893 if (row>155) continue;
5894 if (ktrack0->GetClusterPointer(row)){
5895 AliTPCTrackerPoint *point =ktrack0->GetTrackPoint(row);
5896 shapesum+=point->GetSigmaY()+point->GetSigmaZ();
5899 if (ktrack1->GetClusterPointer(row)){
5900 AliTPCTrackerPoint *point =ktrack1->GetTrackPoint(row);
5901 shapesum+=point->GetSigmaY()+point->GetSigmaZ();
5906 kink->SetShapeFactor(-1.);
5909 kink->SetShapeFactor(shapesum/sum);
5911 // esd->AddKink(kink);
5913 // kink->SetMother(paramm);
5914 //kink->SetDaughter(paramd);
5916 Double_t chi2P2 = paramm.GetParameter()[2]-paramd.GetParameter()[2];
5918 chi2P2/=paramm.GetCovariance()[5]+paramd.GetCovariance()[5];
5919 Double_t chi2P3 = paramm.GetParameter()[3]-paramd.GetParameter()[3];
5921 chi2P3/=paramm.GetCovariance()[9]+paramd.GetCovariance()[9];
5923 if (AliTPCReconstructor::StreamLevel()>1) {
5924 (*fDebugStreamer)<<"kinkLpt"<<
5932 if ( chi2P2+chi2P3<AliTPCReconstructor::GetRecoParam()->GetKinkAngleCutChi2(0)){
5936 kinks->AddLast(kink);
5942 // sort the kinks according quality - and refit them towards vertex
5944 Int_t nkinks = kinks->GetEntriesFast();
5945 Float_t *quality = new Float_t[nkinks];
5946 Int_t *indexes = new Int_t[nkinks];
5947 AliTPCseed **mothers = new AliTPCseed*[nkinks]; memset(mothers, 0, nkinks*sizeof(AliTPCseed*));
5948 AliTPCseed **daughters = new AliTPCseed*[nkinks]; memset(daughters, 0, nkinks*sizeof(AliTPCseed*));
5951 for (Int_t i=0;i<nkinks;i++){
5953 AliKink *kinkl = (AliKink*)kinks->At(i);
5955 // refit kinks towards vertex
5957 Int_t index0 = kinkl->GetIndex(0);
5958 Int_t index1 = kinkl->GetIndex(1);
5959 AliTPCseed * ktrack0 = (AliTPCseed*)array->At(index0);
5960 AliTPCseed * ktrack1 = (AliTPCseed*)array->At(index1);
5962 Int_t sumn=ktrack0->GetNumberOfClusters()+ktrack1->GetNumberOfClusters();
5964 // Refit Kink under if too small angle
5966 if (kinkl->GetAngle(2)<0.05){
5967 kinkl->SetTPCRow0(GetRowNumber(kinkl->GetR()));
5968 Int_t row0 = kinkl->GetTPCRow0();
5969 Int_t drow = Int_t(2.+0.5/(0.05+kinkl->GetAngle(2)));
5972 Int_t last = row0-drow;
5973 if (last<40) last=40;
5974 if (last<ktrack0->GetFirstPoint()+25) last = ktrack0->GetFirstPoint()+25;
5975 AliTPCseed* seed0 = ReSeed(ktrack0,last,kFALSE);
5978 Int_t first = row0+drow;
5979 if (first>130) first=130;
5980 if (first>ktrack1->GetLastPoint()-25) first = TMath::Max(ktrack1->GetLastPoint()-25,30);
5981 AliTPCseed* seed1 = ReSeed(ktrack1,first,kTRUE);
5983 if (seed0 && seed1){
5984 kinkl->SetStatus(1,8);
5985 if (RefitKink(*seed0,*seed1,*kinkl)) kinkl->SetStatus(1,9);
5986 row0 = GetRowNumber(kinkl->GetR());
5987 sumn = seed0->GetNumberOfClusters()+seed1->GetNumberOfClusters();
5988 mothers[i] = new ( NextFreeSeed() ) AliTPCseed(*seed0);
5989 mothers[i]->SetPoolID(fLastSeedID);
5990 daughters[i] = new (NextFreeSeed() ) AliTPCseed(*seed1);
5991 daughters[i]->SetPoolID(fLastSeedID);
5994 delete kinks->RemoveAt(i);
5995 if (seed0) MarkSeedFree( seed0 );
5996 if (seed1) MarkSeedFree( seed1 );
5999 if (kinkl->GetDistance()>0.5 || kinkl->GetR()<110 || kinkl->GetR()>240) {
6000 delete kinks->RemoveAt(i);
6001 if (seed0) MarkSeedFree( seed0 );
6002 if (seed1) MarkSeedFree( seed1 );
6006 MarkSeedFree( seed0 );
6007 MarkSeedFree( seed1 );
6010 if (kinkl) quality[i] = 160*((0.1+kinkl->GetDistance())*(2.-kinkl->GetTPCDensityFactor()))/(sumn+40.); //the longest -clossest will win
6012 TMath::Sort(nkinks,quality,indexes,kFALSE);
6014 //remove double find kinks
6016 for (Int_t ikink0=1;ikink0<nkinks;ikink0++){
6017 AliKink * kink0 = (AliKink*) kinks->At(indexes[ikink0]);
6018 if (!kink0) continue;
6020 for (Int_t ikink1=0;ikink1<ikink0;ikink1++){
6021 kink0 = (AliKink*) kinks->At(indexes[ikink0]);
6022 if (!kink0) continue;
6023 AliKink * kink1 = (AliKink*) kinks->At(indexes[ikink1]);
6024 if (!kink1) continue;
6025 // if not close kink continue
6026 if (TMath::Abs(kink1->GetPosition()[2]-kink0->GetPosition()[2])>10) continue;
6027 if (TMath::Abs(kink1->GetPosition()[1]-kink0->GetPosition()[1])>10) continue;
6028 if (TMath::Abs(kink1->GetPosition()[0]-kink0->GetPosition()[0])>10) continue;
6030 AliTPCseed &mother0 = *mothers[indexes[ikink0]];
6031 AliTPCseed &daughter0 = *daughters[indexes[ikink0]];
6032 AliTPCseed &mother1 = *mothers[indexes[ikink1]];
6033 AliTPCseed &daughter1 = *daughters[indexes[ikink1]];
6034 Int_t row0 = (kink0->GetTPCRow0()+kink1->GetTPCRow0())/2;
6043 for (Int_t i=0;i<row0;i++){
6044 if (mother0.GetClusterIndex(i)>0 && mother1.GetClusterIndex(i)>0){
6047 if (mother0.GetClusterIndex(i)==mother1.GetClusterIndex(i)){
6054 for (Int_t i=row0;i<158;i++){
6055 //if (daughter0.GetClusterIndex(i)>0 && daughter0.GetClusterIndex(i)>0){ // RS: Bug?
6056 if (daughter0.GetClusterIndex(i)>0 && daughter1.GetClusterIndex(i)>0){
6059 if (mother0.GetClusterIndex(i)==mother1.GetClusterIndex(i)){
6065 Float_t ratio = Float_t(same+1)/Float_t(both+1);
6066 Float_t ratiom = Float_t(samem+1)/Float_t(bothm+1);
6067 Float_t ratiod = Float_t(samed+1)/Float_t(bothd+1);
6068 if (ratio>0.3 && ratiom>0.5 &&ratiod>0.5) {
6069 Int_t sum0 = mother0.GetNumberOfClusters()+daughter0.GetNumberOfClusters();
6070 Int_t sum1 = mother1.GetNumberOfClusters()+daughter1.GetNumberOfClusters();
6072 shared[kink0->GetIndex(0)]= kTRUE;
6073 shared[kink0->GetIndex(1)]= kTRUE;
6074 delete kinks->RemoveAt(indexes[ikink0]);
6078 shared[kink1->GetIndex(0)]= kTRUE;
6079 shared[kink1->GetIndex(1)]= kTRUE;
6080 delete kinks->RemoveAt(indexes[ikink1]);
6087 for (Int_t i=0;i<nkinks;i++){
6088 AliKink * kinkl = (AliKink*) kinks->At(indexes[i]);
6089 if (!kinkl) continue;
6090 kinkl->SetTPCRow0(GetRowNumber(kinkl->GetR()));
6091 Int_t index0 = kinkl->GetIndex(0);
6092 Int_t index1 = kinkl->GetIndex(1);
6093 if (circular[index0]||(circular[index1]&&kinkl->GetDistance()>0.2)) continue;
6094 kinkl->SetMultiple(usage[index0],0);
6095 kinkl->SetMultiple(usage[index1],1);
6096 if (kinkl->GetMultiple()[0]+kinkl->GetMultiple()[1]>2) continue;
6097 if (kinkl->GetMultiple()[0]+kinkl->GetMultiple()[1]>0 && quality[indexes[i]]>0.2) continue;
6098 if (kinkl->GetMultiple()[0]+kinkl->GetMultiple()[1]>0 && kinkl->GetDistance()>0.2) continue;
6099 if (circular[index0]||(circular[index1]&&kinkl->GetDistance()>0.1)) continue;
6101 AliTPCseed * ktrack0 = (AliTPCseed*)array->At(index0);
6102 AliTPCseed * ktrack1 = (AliTPCseed*)array->At(index1);
6103 if (!ktrack0 || !ktrack1) continue;
6104 Int_t index = esd->AddKink(kinkl);
6107 if ( ktrack0->GetKinkIndex(0)==0 && ktrack1->GetKinkIndex(0)==0) { //best kink
6108 if (mothers[indexes[i]]->GetNumberOfClusters()>20 && daughters[indexes[i]]->GetNumberOfClusters()>20 &&
6109 (mothers[indexes[i]]->GetNumberOfClusters()+daughters[indexes[i]]->GetNumberOfClusters())>100){
6110 *ktrack0 = *mothers[indexes[i]];
6111 *ktrack1 = *daughters[indexes[i]];
6115 ktrack0->SetKinkIndex(usage[index0],-(index+1));
6116 ktrack1->SetKinkIndex(usage[index1], (index+1));
6121 // Remove tracks corresponding to shared kink's
6123 for (Int_t i=0;i<nentries;i++){
6124 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
6125 if (!track0) continue;
6126 if (track0->GetKinkIndex(0)!=0) continue;
6127 if (shared[i]) MarkSeedFree( array->RemoveAt(i) );
6132 RemoveUsed2(array,0.5,0.4,30);
6134 for (Int_t i=0;i<nentries;i++){
6135 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
6136 if (!track0) continue;
6137 track0->CookdEdx(0.02,0.6);
6141 for (Int_t i=0;i<nentries;i++){
6142 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
6143 if (!track0) continue;
6144 if (track0->Pt()<1.4) continue;
6145 //remove double high momenta tracks - overlapped with kink candidates
6148 for (Int_t icl=track0->GetFirstPoint();icl<track0->GetLastPoint(); icl++){
6149 if (track0->GetClusterPointer(icl)!=0){
6151 if (track0->GetClusterPointer(icl)->IsUsed(10)) ishared++;
6154 if (Float_t(ishared+1)/Float_t(all+1)>0.5) {
6155 MarkSeedFree( array->RemoveAt(i) );
6159 if (track0->GetKinkIndex(0)!=0) continue;
6160 if (track0->GetNumberOfClusters()<80) continue;
6162 AliTPCseed *pmother = new( NextFreeSeed() ) AliTPCseed();
6163 pmother->SetPoolID(fLastSeedID);
6164 AliTPCseed *pdaughter = new( NextFreeSeed() ) AliTPCseed();
6165 pdaughter->SetPoolID(fLastSeedID);
6166 AliKink *pkink = new AliKink;
6168 AliTPCseed & mother = *pmother;
6169 AliTPCseed & daughter = *pdaughter;
6170 AliKink & kinkl = *pkink;
6171 if (CheckKinkPoint(track0,mother,daughter, kinkl)){
6172 if (mother.GetNumberOfClusters()<30||daughter.GetNumberOfClusters()<20) {
6173 MarkSeedFree( pmother );
6174 MarkSeedFree( pdaughter );
6176 continue; //too short tracks
6178 if (mother.Pt()<1.4) {
6179 MarkSeedFree( pmother );
6180 MarkSeedFree( pdaughter );
6184 Int_t row0= kinkl.GetTPCRow0();
6185 if (kinkl.GetDistance()>0.5 || kinkl.GetR()<110. || kinkl.GetR()>240.) {
6186 MarkSeedFree( pmother );
6187 MarkSeedFree( pdaughter );
6192 Int_t index = esd->AddKink(&kinkl);
6193 mother.SetKinkIndex(0,-(index+1));
6194 daughter.SetKinkIndex(0,index+1);
6195 if (mother.GetNumberOfClusters()>50) {
6196 MarkSeedFree( array->RemoveAt(i) );
6197 AliTPCseed* mtc = new( NextFreeSeed() ) AliTPCseed(mother);
6198 mtc->SetPoolID(fLastSeedID);
6199 array->AddAt(mtc,i);
6202 AliTPCseed* mtc = new( NextFreeSeed() ) AliTPCseed(mother);
6203 mtc->SetPoolID(fLastSeedID);
6204 array->AddLast(mtc);
6206 AliTPCseed* dtc = new( NextFreeSeed() ) AliTPCseed(daughter);
6207 dtc->SetPoolID(fLastSeedID);
6208 array->AddLast(dtc);
6209 for (Int_t icl=0;icl<row0;icl++) {
6210 if (mother.GetClusterPointer(icl)) mother.GetClusterPointer(icl)->Use(20);
6213 for (Int_t icl=row0;icl<158;icl++) {
6214 if (daughter.GetClusterPointer(icl)) daughter.GetClusterPointer(icl)->Use(20);
6218 MarkSeedFree( pmother );
6219 MarkSeedFree( pdaughter );
6223 delete [] daughters;
6245 AliInfo(Form("Ncandidates=\t%d\t%d\t%d\t%d\n",esd->GetNumberOfKinks(),ncandidates,ntracks,nall));
6250 Int_t AliTPCtrackerMI::RefitKink(AliTPCseed &mother, AliTPCseed &daughter, const AliESDkink &knk)
6253 // refit kink towards to the vertex
6256 AliKink &kink=(AliKink &)knk;
6258 Int_t row0 = GetRowNumber(kink.GetR());
6259 FollowProlongation(mother,0);
6260 mother.Reset(kFALSE);
6262 FollowProlongation(daughter,row0);
6263 daughter.Reset(kFALSE);
6264 FollowBackProlongation(daughter,158);
6265 daughter.Reset(kFALSE);
6266 Int_t first = TMath::Max(row0-20,30);
6267 Int_t last = TMath::Min(row0+20,140);
6269 const Int_t kNdiv =5;
6270 AliTPCseed param0[kNdiv]; // parameters along the track
6271 AliTPCseed param1[kNdiv]; // parameters along the track
6272 AliKink kinks[kNdiv]; // corresponding kink parameters
6275 for (Int_t irow=0; irow<kNdiv;irow++){
6276 rows[irow] = first +((last-first)*irow)/(kNdiv-1);
6278 // store parameters along the track
6280 for (Int_t irow=0;irow<kNdiv;irow++){
6281 FollowBackProlongation(mother, rows[irow]);
6282 FollowProlongation(daughter,rows[kNdiv-1-irow]);
6283 param0[irow] = mother;
6284 param1[kNdiv-1-irow] = daughter;
6288 for (Int_t irow=0; irow<kNdiv-1;irow++){
6289 if (param0[irow].GetNumberOfClusters()<kNdiv||param1[irow].GetNumberOfClusters()<kNdiv) continue;
6290 kinks[irow].SetMother(param0[irow]);
6291 kinks[irow].SetDaughter(param1[irow]);
6292 kinks[irow].Update();
6295 // choose kink with best "quality"
6297 Double_t mindist = 10000;
6298 for (Int_t irow=0;irow<kNdiv;irow++){
6299 if (param0[irow].GetNumberOfClusters()<20||param1[irow].GetNumberOfClusters()<20) continue;
6300 if (TMath::Abs(kinks[irow].GetR())>240.) continue;
6301 if (TMath::Abs(kinks[irow].GetR())<100.) continue;
6303 Float_t normdist = TMath::Abs(param0[irow].GetX()-kinks[irow].GetR())*(0.1+kink.GetDistance());
6304 normdist/= (param0[irow].GetNumberOfClusters()+param1[irow].GetNumberOfClusters()+40.);
6305 if (normdist < mindist){
6311 if (index==-1) return 0;
6314 param0[index].Reset(kTRUE);
6315 FollowProlongation(param0[index],0);
6317 mother = param0[index];
6318 daughter = param1[index]; // daughter in vertex
6320 kink.SetMother(mother);
6321 kink.SetDaughter(daughter);
6323 kink.SetTPCRow0(GetRowNumber(kink.GetR()));
6324 kink.SetTPCncls(param0[index].GetNumberOfClusters(),0);
6325 kink.SetTPCncls(param1[index].GetNumberOfClusters(),1);
6326 kink.SetLabel(CookLabel(&mother,0.4, 0,kink.GetTPCRow0()),0);
6327 kink.SetLabel(CookLabel(&daughter,0.4, kink.GetTPCRow0(),160),1);
6328 mother.SetLabel(kink.GetLabel(0));
6329 daughter.SetLabel(kink.GetLabel(1));
6335 void AliTPCtrackerMI::UpdateKinkQualityM(AliTPCseed * seed){
6337 // update Kink quality information for mother after back propagation
6339 if (seed->GetKinkIndex(0)>=0) return;
6340 for (Int_t ikink=0;ikink<3;ikink++){
6341 Int_t index = seed->GetKinkIndex(ikink);
6342 if (index>=0) break;
6343 index = TMath::Abs(index)-1;
6344 AliESDkink * kink = fEvent->GetKink(index);
6345 kink->SetTPCDensity(-1,0,0);
6346 kink->SetTPCDensity(1,0,1);
6348 Int_t row0 = kink->GetTPCRow0() - 2 - Int_t( 0.5/ (0.05+kink->GetAngle(2)));
6349 if (row0<15) row0=15;
6351 Int_t row1 = kink->GetTPCRow0() + 2 + Int_t( 0.5/ (0.05+kink->GetAngle(2)));
6352 if (row1>145) row1=145;
6354 Int_t found,foundable,shared;
6355 seed->GetClusterStatistic(0,row0, found, foundable,shared,kFALSE);
6356 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),0,0);
6357 seed->GetClusterStatistic(row1,155, found, foundable,shared,kFALSE);
6358 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),0,1);
6363 void AliTPCtrackerMI::UpdateKinkQualityD(AliTPCseed * seed){
6365 // update Kink quality information for daughter after refit
6367 if (seed->GetKinkIndex(0)<=0) return;
6368 for (Int_t ikink=0;ikink<3;ikink++){
6369 Int_t index = seed->GetKinkIndex(ikink);
6370 if (index<=0) break;
6371 index = TMath::Abs(index)-1;
6372 AliESDkink * kink = fEvent->GetKink(index);
6373 kink->SetTPCDensity(-1,1,0);
6374 kink->SetTPCDensity(-1,1,1);
6376 Int_t row0 = kink->GetTPCRow0() -2 - Int_t( 0.5/ (0.05+kink->GetAngle(2)));
6377 if (row0<15) row0=15;
6379 Int_t row1 = kink->GetTPCRow0() +2 + Int_t( 0.5/ (0.05+kink->GetAngle(2)));
6380 if (row1>145) row1=145;
6382 Int_t found,foundable,shared;
6383 seed->GetClusterStatistic(0,row0, found, foundable,shared,kFALSE);
6384 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),1,0);
6385 seed->GetClusterStatistic(row1,155, found, foundable,shared,kFALSE);
6386 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),1,1);
6392 Int_t AliTPCtrackerMI::CheckKinkPoint(AliTPCseed*seed,AliTPCseed &mother, AliTPCseed &daughter, const AliESDkink &knk)
6395 // check kink point for given track
6396 // if return value=0 kink point not found
6397 // otherwise seed0 correspond to mother particle
6398 // seed1 correspond to daughter particle
6399 // kink parameter of kink point
6400 AliKink &kink=(AliKink &)knk;
6402 Int_t middlerow = (seed->GetFirstPoint()+seed->GetLastPoint())/2;
6403 Int_t first = seed->GetFirstPoint();
6404 Int_t last = seed->GetLastPoint();
6405 if (last-first<20) return 0; // shortest length - 2*30 = 60 pad-rows
6408 AliTPCseed *seed1 = ReSeed(seed,middlerow+20, kTRUE); //middle of chamber
6409 if (!seed1) return 0;
6410 FollowProlongation(*seed1,seed->GetLastPoint()-20);
6411 seed1->Reset(kTRUE);
6412 FollowProlongation(*seed1,158);
6413 seed1->Reset(kTRUE);
6414 last = seed1->GetLastPoint();
6416 AliTPCseed *seed0 = new( NextFreeSeed() ) AliTPCseed(*seed);
6417 seed0->SetPoolID(fLastSeedID);
6418 seed0->Reset(kFALSE);
6421 AliTPCseed param0[20]; // parameters along the track
6422 AliTPCseed param1[20]; // parameters along the track
6423 AliKink kinks[20]; // corresponding kink parameters
6425 for (Int_t irow=0; irow<20;irow++){
6426 rows[irow] = first +((last-first)*irow)/19;
6428 // store parameters along the track
6430 for (Int_t irow=0;irow<20;irow++){
6431 FollowBackProlongation(*seed0, rows[irow]);
6432 FollowProlongation(*seed1,rows[19-irow]);
6433 param0[irow] = *seed0;
6434 param1[19-irow] = *seed1;
6438 for (Int_t irow=0; irow<19;irow++){
6439 kinks[irow].SetMother(param0[irow]);
6440 kinks[irow].SetDaughter(param1[irow]);
6441 kinks[irow].Update();
6444 // choose kink with biggest change of angle
6446 Double_t maxchange= 0;
6447 for (Int_t irow=1;irow<19;irow++){
6448 if (TMath::Abs(kinks[irow].GetR())>240.) continue;
6449 if (TMath::Abs(kinks[irow].GetR())<110.) continue;
6450 Float_t quality = TMath::Abs(kinks[irow].GetAngle(2))/(3.+TMath::Abs(kinks[irow].GetR()-param0[irow].GetX()));
6451 if ( quality > maxchange){
6452 maxchange = quality;
6457 MarkSeedFree( seed0 );
6458 MarkSeedFree( seed1 );
6459 if (index<0) return 0;
6461 Int_t row0 = GetRowNumber(kinks[index].GetR()); //row 0 estimate
6462 seed0 = new( NextFreeSeed() ) AliTPCseed(param0[index]);
6463 seed0->SetPoolID(fLastSeedID);
6464 seed1 = new( NextFreeSeed() ) AliTPCseed(param1[index]);
6465 seed1->SetPoolID(fLastSeedID);
6466 seed0->Reset(kFALSE);
6467 seed1->Reset(kFALSE);
6468 seed0->ResetCovariance(10.);
6469 seed1->ResetCovariance(10.);
6470 FollowProlongation(*seed0,0);
6471 FollowBackProlongation(*seed1,158);
6472 mother = *seed0; // backup mother at position 0
6473 seed0->Reset(kFALSE);
6474 seed1->Reset(kFALSE);
6475 seed0->ResetCovariance(10.);
6476 seed1->ResetCovariance(10.);
6478 first = TMath::Max(row0-20,0);
6479 last = TMath::Min(row0+20,158);
6481 for (Int_t irow=0; irow<20;irow++){
6482 rows[irow] = first +((last-first)*irow)/19;
6484 // store parameters along the track
6486 for (Int_t irow=0;irow<20;irow++){
6487 FollowBackProlongation(*seed0, rows[irow]);
6488 FollowProlongation(*seed1,rows[19-irow]);
6489 param0[irow] = *seed0;
6490 param1[19-irow] = *seed1;
6494 for (Int_t irow=0; irow<19;irow++){
6495 kinks[irow].SetMother(param0[irow]);
6496 kinks[irow].SetDaughter(param1[irow]);
6497 // param0[irow].Dump();
6498 //param1[irow].Dump();
6499 kinks[irow].Update();
6502 // choose kink with biggest change of angle
6505 for (Int_t irow=0;irow<20;irow++){
6506 if (TMath::Abs(kinks[irow].GetR())>250.) continue;
6507 if (TMath::Abs(kinks[irow].GetR())<90.) continue;
6508 Float_t quality = TMath::Abs(kinks[irow].GetAngle(2))/(3.+TMath::Abs(kinks[irow].GetR()-param0[irow].GetX()));
6509 if ( quality > maxchange){
6510 maxchange = quality;
6517 if (index==-1 || param0[index].GetNumberOfClusters()+param1[index].GetNumberOfClusters()<100){
6518 MarkSeedFree( seed0 );
6519 MarkSeedFree( seed1 );
6523 // Float_t anglesigma = TMath::Sqrt(param0[index].fC22+param0[index].fC33+param1[index].fC22+param1[index].fC33);
6525 kink.SetMother(param0[index]);
6526 kink.SetDaughter(param1[index]);
6529 Double_t chi2P2 = param0[index].GetParameter()[2]-param1[index].GetParameter()[2];
6531 chi2P2/=param0[index].GetCovariance()[5]+param1[index].GetCovariance()[5];
6532 Double_t chi2P3 = param0[index].GetParameter()[3]-param1[index].GetParameter()[3];
6534 chi2P3/=param0[index].GetCovariance()[9]+param1[index].GetCovariance()[9];
6536 if (AliTPCReconstructor::StreamLevel()>1) {
6537 (*fDebugStreamer)<<"kinkHpt"<<
6540 "p0.="<<¶m0[index]<<
6541 "p1.="<<¶m1[index]<<
6545 if ( chi2P2+chi2P3<AliTPCReconstructor::GetRecoParam()->GetKinkAngleCutChi2(0)){
6546 MarkSeedFree( seed0 );
6547 MarkSeedFree( seed1 );
6552 row0 = GetRowNumber(kink.GetR());
6553 kink.SetTPCRow0(row0);
6554 kink.SetLabel(CookLabel(seed0,0.5,0,row0),0);
6555 kink.SetLabel(CookLabel(seed1,0.5,row0,158),1);
6556 kink.SetIndex(-10,0);
6557 kink.SetIndex(int(param0[index].GetNumberOfClusters()+param1[index].GetNumberOfClusters()),1);
6558 kink.SetTPCncls(param0[index].GetNumberOfClusters(),0);
6559 kink.SetTPCncls(param1[index].GetNumberOfClusters(),1);
6562 // new (&mother) AliTPCseed(param0[index]);
6563 daughter = param1[index];
6564 daughter.SetLabel(kink.GetLabel(1));
6565 param0[index].Reset(kTRUE);
6566 FollowProlongation(param0[index],0);
6567 mother = param0[index];
6568 mother.SetLabel(kink.GetLabel(0));
6569 if ( chi2P2+chi2P3<AliTPCReconstructor::GetRecoParam()->GetKinkAngleCutChi2(1)){
6572 MarkSeedFree( seed0 );
6573 MarkSeedFree( seed1 );
6581 AliTPCseed* AliTPCtrackerMI::ReSeed(AliTPCseed *t)
6584 // reseed - refit - track
6587 // Int_t last = fSectors->GetNRows()-1;
6589 if (fSectors == fOuterSec){
6590 first = TMath::Max(first, t->GetFirstPoint()-fInnerSec->GetNRows());
6594 first = t->GetFirstPoint();
6596 AliTPCseed * seed = MakeSeed(t,0.1,0.5,0.9);
6597 FollowBackProlongation(*t,fSectors->GetNRows()-1);
6599 FollowProlongation(*t,first);
6609 //_____________________________________________________________________________
6610 Int_t AliTPCtrackerMI::ReadSeeds(const TFile *inp) {
6611 //-----------------------------------------------------------------
6612 // This function reades track seeds.
6613 //-----------------------------------------------------------------
6614 TDirectory *savedir=gDirectory;
6616 TFile *in=(TFile*)inp;
6617 if (!in->IsOpen()) {
6618 cerr<<"AliTPCtrackerMI::ReadSeeds(): input file is not open !\n";
6623 TTree *seedTree=(TTree*)in->Get("Seeds");
6625 cerr<<"AliTPCtrackerMI::ReadSeeds(): ";
6626 cerr<<"can't get a tree with track seeds !\n";
6629 AliTPCtrack *seed=new AliTPCtrack;
6630 seedTree->SetBranchAddress("tracks",&seed);
6632 if (fSeeds==0) fSeeds=new TObjArray(15000);
6634 Int_t n=(Int_t)seedTree->GetEntries();
6635 for (Int_t i=0; i<n; i++) {
6636 seedTree->GetEvent(i);
6637 AliTPCseed* sdc = new( NextFreeSeed() ) AliTPCseed(*seed/*,seed->GetAlpha()*/);
6638 sdc->SetPoolID(fLastSeedID);
6639 fSeeds->AddLast(sdc);
6642 delete seed; // RS: this seed is not from the pool, delete it !!!
6648 Int_t AliTPCtrackerMI::Clusters2Tracks (AliESDEvent *const esd)
6651 // clusters to tracks
6652 if (fSeeds) DeleteSeeds();
6653 else ResetSeedsPool();
6656 AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ;
6657 transform->SetCurrentTimeStamp( esd->GetTimeStamp());
6658 transform->SetCurrentRun(esd->GetRunNumber());
6661 if (!fSeeds) return 1;
6663 if (AliTPCReconstructor::StreamLevel()>3) DumpClusters(0,fSeeds);
6669 //_____________________________________________________________________________
6670 Int_t AliTPCtrackerMI::Clusters2Tracks() {
6671 //-----------------------------------------------------------------
6672 // This is a track finder.
6673 //-----------------------------------------------------------------
6674 TDirectory *savedir=gDirectory;
6678 fSeeds = Tracking();
6681 Info("Clusters2Tracks","Time for tracking: \t");timer.Print();timer.Start();
6683 //activate again some tracks
6684 for (Int_t i=0; i<fSeeds->GetEntriesFast(); i++) {
6685 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
6687 Int_t nc=t.GetNumberOfClusters();
6689 MarkSeedFree( fSeeds->RemoveAt(i) );
6693 if (pt->GetRemoval()==10) {
6694 if (pt->GetDensityFirst(20)>0.8 || pt->GetDensityFirst(30)>0.8 || pt->GetDensityFirst(40)>0.7)
6695 pt->Desactivate(10); // make track again active // MvL: should be 0 ?
6697 pt->Desactivate(20);
6698 MarkSeedFree( fSeeds->RemoveAt(i) );
6703 RemoveUsed2(fSeeds,0.85,0.85,0);
6704 if (AliTPCReconstructor::GetRecoParam()->GetDoKinks()) FindKinks(fSeeds,fEvent);
6705 //FindCurling(fSeeds, fEvent,0);
6706 if (AliTPCReconstructor::StreamLevel()>5) FindMultiMC(fSeeds, fEvent,-1); // find multi found tracks
6707 RemoveUsed2(fSeeds,0.5,0.4,20);
6708 FindSplitted(fSeeds, fEvent,0); // find multi found tracks
6709 if (AliTPCReconstructor::StreamLevel()>5) FindMultiMC(fSeeds, fEvent,0); // find multi found tracks
6712 // // refit short tracks
6714 Int_t nseed=fSeeds->GetEntriesFast();
6717 for (Int_t i=0; i<nseed; i++) {
6718 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
6720 Int_t nc=t.GetNumberOfClusters();
6722 MarkSeedFree( fSeeds->RemoveAt(i) );
6725 CookLabel(pt,0.1); //For comparison only
6726 //if ((pt->IsActive() || (pt->fRemoval==10) )&& nc>50 &&pt->GetNumberOfClusters()>0.4*pt->fNFoundable){
6727 if ((pt->IsActive() || (pt->GetRemoval()==10) )){
6729 if (fDebug>0) cerr<<found<<'\r';
6733 MarkSeedFree( fSeeds->RemoveAt(i) );
6737 //RemoveOverlap(fSeeds,0.99,7,kTRUE);
6739 //RemoveUsed(fSeeds,0.9,0.9,6);
6741 nseed=fSeeds->GetEntriesFast();
6743 for (Int_t i=0; i<nseed; i++) {
6744 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
6746 Int_t nc=t.GetNumberOfClusters();
6748 MarkSeedFree( fSeeds->RemoveAt(i) );
6752 t.CookdEdx(0.02,0.6);
6753 // CheckKinkPoint(&t,0.05);
6754 //if ((pt->IsActive() || (pt->fRemoval==10) )&& nc>50 &&pt->GetNumberOfClusters()>0.4*pt->fNFoundable){
6755 if ((pt->IsActive() || (pt->GetRemoval()==10) )){
6763 MarkSeedFree( fSeeds->RemoveAt(i) );
6764 //AliTPCseed * seed1 = ReSeed(pt,0.05,0.5,1);
6766 // FollowProlongation(*seed1,0);
6767 // Int_t n = seed1->GetNumberOfClusters();
6768 // printf("fP4\t%f\t%f\n",seed1->GetC(),pt->GetC());
6769 // printf("fN\t%d\t%d\n", seed1->GetNumberOfClusters(),pt->GetNumberOfClusters());
6772 //AliTPCseed * seed2 = ReSeed(pt,0.95,0.5,0.05);
6776 SortTracks(fSeeds, 1);
6780 PrepareForBackProlongation(fSeeds,5.);
6781 PropagateBack(fSeeds);
6782 printf("Time for back propagation: \t");timer.Print();timer.Start();
6786 PrepareForProlongation(fSeeds,5.);
6787 PropagateForard2(fSeeds);
6789 printf("Time for FORWARD propagation: \t");timer.Print();timer.Start();
6790 // RemoveUsed(fSeeds,0.7,0.7,6);
6791 //RemoveOverlap(fSeeds,0.9,7,kTRUE);
6793 nseed=fSeeds->GetEntriesFast();
6795 for (Int_t i=0; i<nseed; i++) {
6796 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
6798 Int_t nc=t.GetNumberOfClusters();
6800 MarkSeedFree( fSeeds->RemoveAt(i) );
6803 t.CookdEdx(0.02,0.6);
6804 // CookLabel(pt,0.1); //For comparison only
6805 //if ((pt->IsActive() || (pt->fRemoval==10) )&& nc>50 &&pt->GetNumberOfClusters()>0.4*pt->fNFoundable){
6806 if ((pt->IsActive() || (pt->fRemoval==10) )){
6807 cerr<<found++<<'\r';
6810 MarkSeedFree( fSeeds->RemoveAt(i) );
6815 // fNTracks = found;
6817 Info("Clusters2Tracks","Time for overlap removal, track writing and dedx cooking: \t"); timer.Print();timer.Start();
6820 // cerr<<"Number of found tracks : "<<"\t"<<found<<endl;
6821 Info("Clusters2Tracks","Number of found tracks %d",found);
6823 // UnloadClusters();
6828 void AliTPCtrackerMI::Tracking(TObjArray * arr)
6831 // tracking of the seeds
6834 fSectors = fOuterSec;
6835 ParallelTracking(arr,150,63);
6836 fSectors = fOuterSec;
6837 ParallelTracking(arr,63,0);
6840 TObjArray * AliTPCtrackerMI::Tracking(Int_t seedtype, Int_t i1, Int_t i2, Float_t cuts[4], Float_t dy, Int_t dsec)
6845 static TObjArray arrTracks;
6846 TObjArray * arr = &arrTracks;
6848 fSectors = fOuterSec;
6851 for (Int_t sec=0;sec<fkNOS;sec++){
6852 if (seedtype==3) MakeSeeds3(arr,sec,i1,i2,cuts,dy, dsec);
6853 if (seedtype==4) MakeSeeds5(arr,sec,i1,i2,cuts,dy);
6854 if (seedtype==2) MakeSeeds2(arr,sec,i1,i2,cuts,dy);
6857 Info("Tracking","\nSeeding - %d\t%d\t%d\t%d\n",seedtype,i1,i2,arr->GetEntriesFast());
6869 TObjArray * AliTPCtrackerMI::Tracking()
6873 if (AliTPCReconstructor::GetRecoParam()->GetSpecialSeeding()) return TrackingSpecial();
6876 Int_t nup=fOuterSec->GetNRows()+fInnerSec->GetNRows();
6878 TObjArray * seeds = new TObjArray;
6880 Int_t fLastSeedRowSec=AliTPCReconstructor::GetRecoParam()->GetLastSeedRowSec();
6881 Int_t gapPrim = AliTPCReconstructor::GetRecoParam()->GetSeedGapPrim();
6882 Int_t gapSec = AliTPCReconstructor::GetRecoParam()->GetSeedGapSec();
6890 Float_t fnumber = 3.0;
6891 Float_t fdensity = 3.0;
6896 for (Int_t delta = 0; delta<18; delta+=gapPrim){
6900 arr = Tracking(3,nup-1-delta,nup-1-delta-gap,cuts,-1,1);
6901 SumTracks(seeds,arr);
6902 SignClusters(seeds,fnumber,fdensity);
6904 for (Int_t i=2;i<6;i+=2){
6905 // seed high pt tracks
6908 arr = Tracking(3,nup-i-delta,nup-i-delta-gap,cuts,-1,0);
6909 SumTracks(seeds,arr);
6910 SignClusters(seeds,fnumber,fdensity);
6915 // RemoveUsed(seeds,0.9,0.9,1);
6916 // UnsignClusters();
6917 // SignClusters(seeds,fnumber,fdensity);
6921 for (Int_t delta = 20; delta<120; delta+=gapPrim){
6923 // seed high pt tracks
6927 arr = Tracking(3,nup-delta,nup-delta-gap,cuts,-1);
6928 SumTracks(seeds,arr);
6929 SignClusters(seeds,fnumber,fdensity);
6934 arr = Tracking(3,nup-delta-5,nup-delta-5-gap,cuts,-1);
6935 SumTracks(seeds,arr);
6936 SignClusters(seeds,fnumber,fdensity);
6947 Info("Tracking()","\n\nPrimary seeding\t%d\n\n",seeds->GetEntriesFast());
6951 // RemoveUsed(seeds,0.75,0.75,1);
6953 //SignClusters(seeds,fnumber,fdensity);
6962 arr = Tracking(4,nup-1,nup-1-gap,cuts,-1);
6963 SumTracks(seeds,arr);
6964 SignClusters(seeds,fnumber,fdensity);
6966 arr = Tracking(4,nup-2,nup-2-gap,cuts,-1);
6967 SumTracks(seeds,arr);
6968 SignClusters(seeds,fnumber,fdensity);
6970 arr = Tracking(4,nup-3,nup-3-gap,cuts,-1);
6971 SumTracks(seeds,arr);
6972 SignClusters(seeds,fnumber,fdensity);
6974 arr = Tracking(4,nup-5,nup-5-gap,cuts,-1);
6975 SumTracks(seeds,arr);
6976 SignClusters(seeds,fnumber,fdensity);
6978 arr = Tracking(4,nup-7,nup-7-gap,cuts,-1);
6979 SumTracks(seeds,arr);
6980 SignClusters(seeds,fnumber,fdensity);
6983 arr = Tracking(4,nup-9,nup-9-gap,cuts,-1);
6984 SumTracks(seeds,arr);
6985 SignClusters(seeds,fnumber,fdensity);
6989 for (Int_t delta = 9; delta<30; delta+=gapSec){
6995 arr = Tracking(4,nup-1-delta,nup-1-delta-gap,cuts,-1);
6996 SumTracks(seeds,arr);
6997 SignClusters(seeds,fnumber,fdensity);
6999 arr = Tracking(4,nup-3-delta,nup-5-delta-gap,cuts,4);
7000 SumTracks(seeds,arr);
7001 SignClusters(seeds,fnumber,fdensity);
7014 for (Int_t delta = 30; delta<fLastSeedRowSec; delta+=gapSec){
7020 arr = Tracking(4,nup-1-delta,nup-1-delta-gap,cuts,-1);
7021 SumTracks(seeds,arr);
7022 SignClusters(seeds,fnumber,fdensity);
7024 arr = Tracking(4,nup-5-delta,nup-5-delta-gap,cuts,5 );
7025 SumTracks(seeds,arr);
7026 SignClusters(seeds,fnumber,fdensity);
7030 Info("Tracking()","\n\nSecondary seeding\t%d\n\n",seeds->GetEntriesFast());
7041 TObjArray * AliTPCtrackerMI::TrackingSpecial()
7044 // seeding adjusted for laser and cosmic tests - short tracks with big inclination angle
7045 // no primary vertex seeding tried
7049 Int_t nup=fOuterSec->GetNRows()+fInnerSec->GetNRows();
7051 TObjArray * seeds = new TObjArray;
7056 Float_t fnumber = 3.0;
7057 Float_t fdensity = 3.0;
7060 cuts[0] = AliTPCReconstructor::GetRecoParam()->GetMaxC(); // max curvature
7061 cuts[1] = 3.5; // max tan(phi) angle for seeding
7062 cuts[2] = 3.; // not used (cut on z primary vertex)
7063 cuts[3] = 3.5; // max tan(theta) angle for seeding
7065 for (Int_t delta = 0; nup-delta-gap-1>0; delta+=3){
7067 arr = Tracking(4,nup-1-delta,nup-1-delta-gap,cuts,-1);
7068 SumTracks(seeds,arr);
7069 SignClusters(seeds,fnumber,fdensity);
7073 Info("Tracking()","\n\nSecondary seeding\t%d\n\n",seeds->GetEntriesFast());
7084 void AliTPCtrackerMI::SumTracks(TObjArray *arr1,TObjArray *&arr2)
7087 //sum tracks to common container
7088 //remove suspicious tracks
7089 // RS: Attention: supplied tracks come in the static array, don't delete them
7090 Int_t nseed = arr2->GetEntriesFast();
7091 for (Int_t i=0;i<nseed;i++){
7092 AliTPCseed *pt=(AliTPCseed*)arr2->UncheckedAt(i);
7095 // remove tracks with too big curvature
7097 if (TMath::Abs(pt->GetC())>AliTPCReconstructor::GetRecoParam()->GetMaxC()){
7098 MarkSeedFree( arr2->RemoveAt(i) );
7101 // REMOVE VERY SHORT TRACKS
7102 if (pt->GetNumberOfClusters()<20){
7103 MarkSeedFree( arr2->RemoveAt(i) );
7106 // NORMAL ACTIVE TRACK
7107 if (pt->IsActive()){
7108 arr1->AddLast(arr2->RemoveAt(i));
7111 //remove not usable tracks
7112 if (pt->GetRemoval()!=10){
7113 MarkSeedFree( arr2->RemoveAt(i) );
7117 // ENABLE ONLY ENOUGH GOOD STOPPED TRACKS
7118 if (pt->GetDensityFirst(20)>0.8 || pt->GetDensityFirst(30)>0.8 || pt->GetDensityFirst(40)>0.7)
7119 arr1->AddLast(arr2->RemoveAt(i));
7121 MarkSeedFree( arr2->RemoveAt(i) );
7125 // delete arr2; arr2 = 0; // RS: this is static array, don't delete it
7130 void AliTPCtrackerMI::ParallelTracking(TObjArray *const arr, Int_t rfirst, Int_t rlast)
7133 // try to track in parralel
7135 Int_t nseed=arr->GetEntriesFast();
7136 //prepare seeds for tracking
7137 for (Int_t i=0; i<nseed; i++) {
7138 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i), &t=*pt;
7140 if (!t.IsActive()) continue;
7141 // follow prolongation to the first layer
7142 if ( (fSectors ==fInnerSec) || (t.GetFirstPoint()-fkParam->GetNRowLow()>rfirst+1) )
7143 FollowProlongation(t, rfirst+1);
7148 for (Int_t nr=rfirst; nr>=rlast; nr--){
7149 if (nr<fInnerSec->GetNRows())
7150 fSectors = fInnerSec;
7152 fSectors = fOuterSec;
7153 // make indexes with the cluster tracks for given
7155 // find nearest cluster
7156 for (Int_t i=0; i<nseed; i++) {
7157 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i), &t=*pt;
7159 if (nr==80) pt->UpdateReference();
7160 if (!pt->IsActive()) continue;
7161 // if ( (fSectors ==fOuterSec) && (pt->fFirstPoint-fkParam->GetNRowLow())<nr) continue;
7162 if (pt->GetRelativeSector()>17) {
7165 UpdateClusters(t,nr);
7167 // prolonagate to the nearest cluster - if founded
7168 for (Int_t i=0; i<nseed; i++) {
7169 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
7171 if (!pt->IsActive()) continue;
7172 // if ((fSectors ==fOuterSec) && (pt->fFirstPoint-fkParam->GetNRowLow())<nr) continue;
7173 if (pt->GetRelativeSector()>17) {
7176 FollowToNextCluster(*pt,nr);
7181 void AliTPCtrackerMI::PrepareForBackProlongation(const TObjArray *const arr,Float_t fac) const
7185 // if we use TPC track itself we have to "update" covariance
7187 Int_t nseed= arr->GetEntriesFast();
7188 for (Int_t i=0;i<nseed;i++){
7189 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
7193 //rotate to current local system at first accepted point
7194 Int_t index = pt->GetClusterIndex2(pt->GetFirstPoint());
7195 Int_t sec = (index&0xff000000)>>24;
7197 Float_t angle1 = fInnerSec->GetAlpha()*sec+fInnerSec->GetAlphaShift();
7198 if (angle1>TMath::Pi())
7199 angle1-=2.*TMath::Pi();
7200 Float_t angle2 = pt->GetAlpha();
7202 if (TMath::Abs(angle1-angle2)>0.001){
7203 if (!pt->Rotate(angle1-angle2)) return;
7204 //angle2 = pt->GetAlpha();
7205 //pt->fRelativeSector = pt->GetAlpha()/fInnerSec->GetAlpha();
7206 //if (pt->GetAlpha()<0)
7207 // pt->fRelativeSector+=18;
7208 //sec = pt->fRelativeSector;
7217 void AliTPCtrackerMI::PrepareForProlongation(TObjArray *const arr, Float_t fac) const
7221 // if we use TPC track itself we have to "update" covariance
7223 Int_t nseed= arr->GetEntriesFast();
7224 for (Int_t i=0;i<nseed;i++){
7225 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
7228 pt->SetFirstPoint(pt->GetLastPoint());
7236 Int_t AliTPCtrackerMI::PropagateBack(const TObjArray *const arr)
7239 // make back propagation
7241 Int_t nseed= arr->GetEntriesFast();
7242 for (Int_t i=0;i<nseed;i++){
7243 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
7244 if (pt&& pt->GetKinkIndex(0)<=0) {
7245 //AliTPCseed *pt2 = new AliTPCseed(*pt);
7246 fSectors = fInnerSec;
7247 //FollowBackProlongation(*pt,fInnerSec->GetNRows()-1);
7248 //fSectors = fOuterSec;
7249 FollowBackProlongation(*pt,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1,1);
7250 //if (pt->GetNumberOfClusters()<(pt->fEsd->GetTPCclusters(0)) ){
7251 // Error("PropagateBack","Not prolonged track %d",pt->GetLabel());
7252 // FollowBackProlongation(*pt2,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1);
7255 if (pt&& pt->GetKinkIndex(0)>0) {
7256 AliESDkink * kink = fEvent->GetKink(pt->GetKinkIndex(0)-1);
7257 pt->SetFirstPoint(kink->GetTPCRow0());
7258 fSectors = fInnerSec;
7259 FollowBackProlongation(*pt,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1,1);
7267 Int_t AliTPCtrackerMI::PropagateForward2(const TObjArray *const arr)
7270 // make forward propagation
7272 Int_t nseed= arr->GetEntriesFast();
7274 for (Int_t i=0;i<nseed;i++){
7275 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
7277 FollowProlongation(*pt,0,1,1);
7286 Int_t AliTPCtrackerMI::PropagateForward()
7289 // propagate track forward
7291 Int_t nseed = fSeeds->GetEntriesFast();
7292 for (Int_t i=0;i<nseed;i++){
7293 AliTPCseed *pt = (AliTPCseed*)fSeeds->UncheckedAt(i);
7295 AliTPCseed &t = *pt;
7296 Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
7297 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
7298 if (alpha < 0. ) alpha += 2.*TMath::Pi();
7299 t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN);
7303 fSectors = fOuterSec;
7304 ParallelTracking(fSeeds,fOuterSec->GetNRows()+fInnerSec->GetNRows()-1,fInnerSec->GetNRows());
7305 fSectors = fInnerSec;
7306 ParallelTracking(fSeeds,fInnerSec->GetNRows()-1,0);
7315 Int_t AliTPCtrackerMI::PropagateBack(AliTPCseed *const pt, Int_t row0, Int_t row1)
7318 // make back propagation, in between row0 and row1
7322 fSectors = fInnerSec;
7325 if (row1<fSectors->GetNRows())
7328 r1 = fSectors->GetNRows()-1;
7330 if (row0<fSectors->GetNRows()&& r1>0 )
7331 FollowBackProlongation(*pt,r1);
7332 if (row1<=fSectors->GetNRows())
7335 r1 = row1 - fSectors->GetNRows();
7336 if (r1<=0) return 0;
7337 if (r1>=fOuterSec->GetNRows()) return 0;
7338 fSectors = fOuterSec;
7339 return FollowBackProlongation(*pt,r1);
7347 void AliTPCtrackerMI::GetShape(AliTPCseed * seed, Int_t row)
7349 // gets cluster shape
7351 AliTPCClusterParam * clparam = AliTPCcalibDB::Instance()->GetClusterParam();
7352 Float_t zdrift = TMath::Abs((fkParam->GetZLength(0)-TMath::Abs(seed->GetZ())));
7353 Int_t type = (seed->GetSector() < fkParam->GetNSector()/2) ? 0: (row>126) ? 1:2;
7354 Double_t angulary = seed->GetSnp();
7356 if (TMath::Abs(angulary)>AliTPCReconstructor::GetMaxSnpTracker()) {
7357 angulary = TMath::Sign(AliTPCReconstructor::GetMaxSnpTracker(),angulary);
7360 angulary = angulary*angulary/((1.-angulary)*(1.+angulary));
7361 Double_t angularz = seed->GetTgl()*seed->GetTgl()*(1.+angulary);
7363 Double_t sigmay = clparam->GetRMS0(0,type,zdrift,TMath::Sqrt(TMath::Abs(angulary)));
7364 Double_t sigmaz = clparam->GetRMS0(1,type,zdrift,TMath::Sqrt(TMath::Abs(angularz)));
7365 seed->SetCurrentSigmaY2(sigmay*sigmay);
7366 seed->SetCurrentSigmaZ2(sigmaz*sigmaz);
7367 // Float_t sd2 = TMath::Abs((fkParam->GetZLength(0)-TMath::Abs(seed->GetZ())))*fkParam->GetDiffL()*fkParam->GetDiffL();
7368 // // Float_t padlength = fkParam->GetPadPitchLength(seed->fSector);
7369 // Float_t padlength = GetPadPitchLength(row);
7371 // Float_t sresy = (seed->GetSector() < fkParam->GetNSector()/2) ? 0.2 :0.3;
7372 // seed->SetCurrentSigmaY2(sd2+padlength*padlength*angulary/12.+sresy*sresy);
7374 // Float_t sresz = fkParam->GetZSigma();
7375 // seed->SetCurrentSigmaZ2(sd2+padlength*padlength*angularz*angularz*(1+angulary)/12.+sresz*sresz);
7377 Float_t wy = GetSigmaY(seed);
7378 Float_t wz = GetSigmaZ(seed);
7381 if (TMath::Abs(wy/seed->fCurrentSigmaY2-1)>0.0001 || TMath::Abs(wz/seed->fCurrentSigmaZ2-1)>0.0001 ){
7382 printf("problem\n");
7389 //__________________________________________________________________________
7390 void AliTPCtrackerMI::CookLabel(AliKalmanTrack *tk, Float_t wrong) const {
7391 //--------------------------------------------------------------------
7392 //This function "cooks" a track label. If label<0, this track is fake.
7393 //--------------------------------------------------------------------
7394 AliTPCseed * t = dynamic_cast<AliTPCseed*>(tk);
7396 printf("%s:%d wrong type \n",(char*)__FILE__,__LINE__);
7400 Int_t noc=t->GetNumberOfClusters();
7402 //printf("\nnot founded prolongation\n\n\n");
7408 AliTPCclusterMI *clusters[160];
7410 for (Int_t i=0;i<160;i++) {
7417 for (i=0; i<160 && current<noc; i++) {
7419 Int_t index=t->GetClusterIndex2(i);
7420 if (index<=0) continue;
7421 if (index&0x8000) continue;
7423 //clusters[current]=GetClusterMI(index);
7424 if (t->GetClusterPointer(i)){
7425 clusters[current]=t->GetClusterPointer(i);
7431 Int_t lab=123456789;
7432 for (i=0; i<noc; i++) {
7433 AliTPCclusterMI *c=clusters[i];
7435 lab=TMath::Abs(c->GetLabel(0));
7437 for (j=0; j<noc; j++) if (lb[j]==lab || mx[j]==0) break;
7443 for (i=0; i<noc; i++) if (mx[i]>max) {max=mx[i]; lab=lb[i];}
7445 for (i=0; i<noc; i++) {
7446 AliTPCclusterMI *c=clusters[i];
7448 if (TMath::Abs(c->GetLabel(1)) == lab ||
7449 TMath::Abs(c->GetLabel(2)) == lab ) max++;
7451 if (noc<=0) { lab=-1; return;}
7452 if ((1.- Float_t(max)/(noc)) > wrong) lab=-lab;
7455 Int_t tail=Int_t(0.10*noc);
7458 for (i=1; i<160&&ind<tail; i++) {
7459 // AliTPCclusterMI *c=clusters[noc-i];
7460 AliTPCclusterMI *c=clusters[i];
7462 if (lab == TMath::Abs(c->GetLabel(0)) ||
7463 lab == TMath::Abs(c->GetLabel(1)) ||
7464 lab == TMath::Abs(c->GetLabel(2))) max++;
7467 if (max < Int_t(0.5*tail)) lab=-lab;
7474 //delete[] clusters;
7478 //__________________________________________________________________________
7479 Int_t AliTPCtrackerMI::CookLabel(AliTPCseed *const t, Float_t wrong,Int_t first, Int_t last) const {
7480 //--------------------------------------------------------------------
7481 //This function "cooks" a track label. If label<0, this track is fake.
7482 //--------------------------------------------------------------------
7483 Int_t noc=t->GetNumberOfClusters();
7485 //printf("\nnot founded prolongation\n\n\n");
7491 AliTPCclusterMI *clusters[160];
7493 for (Int_t i=0;i<160;i++) {
7500 for (i=0; i<160 && current<noc; i++) {
7501 if (i<first) continue;
7502 if (i>last) continue;
7503 Int_t index=t->GetClusterIndex2(i);
7504 if (index<=0) continue;
7505 if (index&0x8000) continue;
7507 //clusters[current]=GetClusterMI(index);
7508 if (t->GetClusterPointer(i)){
7509 clusters[current]=t->GetClusterPointer(i);
7514 //if (noc<5) return -1;
7515 Int_t lab=123456789;
7516 for (i=0; i<noc; i++) {
7517 AliTPCclusterMI *c=clusters[i];
7519 lab=TMath::Abs(c->GetLabel(0));
7521 for (j=0; j<noc; j++) if (lb[j]==lab || mx[j]==0) break;
7527 for (i=0; i<noc; i++) if (mx[i]>max) {max=mx[i]; lab=lb[i];}
7529 for (i=0; i<noc; i++) {
7530 AliTPCclusterMI *c=clusters[i];
7532 if (TMath::Abs(c->GetLabel(1)) == lab ||
7533 TMath::Abs(c->GetLabel(2)) == lab ) max++;
7535 if (noc<=0) { lab=-1; return -1;}
7536 if ((1.- Float_t(max)/(noc)) > wrong) lab=-lab;
7539 Int_t tail=Int_t(0.10*noc);
7542 for (i=1; i<160&&ind<tail; i++) {
7543 // AliTPCclusterMI *c=clusters[noc-i];
7544 AliTPCclusterMI *c=clusters[i];
7546 if (lab == TMath::Abs(c->GetLabel(0)) ||
7547 lab == TMath::Abs(c->GetLabel(1)) ||
7548 lab == TMath::Abs(c->GetLabel(2))) max++;
7551 if (max < Int_t(0.5*tail)) lab=-lab;
7554 // t->SetLabel(lab);
7558 //delete[] clusters;
7562 Int_t AliTPCtrackerMI::GetRowNumber(Double_t x[3]) const
7564 //return pad row number for given x vector
7565 Float_t phi = TMath::ATan2(x[1],x[0]);
7566 if(phi<0) phi=2.*TMath::Pi()+phi;
7567 // Get the local angle in the sector philoc
7568 const Float_t kRaddeg = 180/3.14159265358979312;
7569 Float_t phiangle = (Int_t (phi*kRaddeg/20.) + 0.5)*20./kRaddeg;
7570 Double_t localx = x[0]*TMath::Cos(phiangle)-x[1]*TMath::Sin(phiangle);
7571 return GetRowNumber(localx);
7576 void AliTPCtrackerMI::MakeESDBitmaps(AliTPCseed *t, AliESDtrack *esd)
7578 //-----------------------------------------------------------------------
7579 // Fill the cluster and sharing bitmaps of the track
7580 //-----------------------------------------------------------------------
7582 Int_t firstpoint = 0;
7583 Int_t lastpoint = 159;
7584 AliTPCTrackerPoint *point;
7585 AliTPCclusterMI *cluster;
7588 TBits clusterMap(159);
7589 TBits sharedMap(159);
7591 for (int iter=firstpoint; iter<lastpoint; iter++) {
7592 // Change to cluster pointers to see if we have a cluster at given padrow
7594 cluster = t->GetClusterPointer(iter);
7596 clusterMap.SetBitNumber(iter, kTRUE);
7597 point = t->GetTrackPoint(iter);
7598 if (point->IsShared())
7599 sharedMap.SetBitNumber(iter,kTRUE);
7601 if (t->GetClusterIndex(iter) >= 0 && (t->GetClusterIndex(iter) & 0x8000) == 0) {
7602 fitMap.SetBitNumber(iter, kTRUE);
7606 esd->SetTPCClusterMap(clusterMap);
7607 esd->SetTPCSharedMap(sharedMap);
7608 esd->SetTPCFitMap(fitMap);
7609 if (nclsf != t->GetNumberOfClusters())
7610 AliWarning(Form("Inconsistency between ncls %d and indices %d (found %d)",t->GetNumberOfClusters(),nclsf,esd->GetTPCClusterMap().CountBits()));
7613 Bool_t AliTPCtrackerMI::IsFindable(AliTPCseed & track){
7615 // return flag if there is findable cluster at given position
7618 Float_t z = track.GetZ();
7620 if (TMath::Abs(z)<(AliTPCReconstructor::GetCtgRange()*track.GetX()+kDeltaZ) &&
7621 TMath::Abs(z)<fkParam->GetZLength(0) &&
7622 (TMath::Abs(track.GetSnp())<AliTPCReconstructor::GetMaxSnpTracker()))
7628 void AliTPCtrackerMI::AddCovariance(AliTPCseed * seed){
7630 // Adding systematic error estimate to the covariance matrix
7631 // !!!! the systematic error for element 4 is in 1/GeV
7632 // 03.03.2012 MI changed in respect to the previous versions
7633 const Double_t *param = AliTPCReconstructor::GetRecoParam()->GetSystematicError();
7635 // use only the diagonal part if not specified otherwise
7636 if (!AliTPCReconstructor::GetRecoParam()->GetUseSystematicCorrelation()) return AddCovarianceAdd(seed);
7638 Double_t *covarS= (Double_t*)seed->GetCovariance();
7639 Double_t factor[5]={1,1,1,1,1};
7640 factor[0]= TMath::Sqrt(TMath::Abs((covarS[0] + param[0]*param[0])/covarS[0]));
7641 factor[1]= TMath::Sqrt(TMath::Abs((covarS[2] + param[1]*param[1])/covarS[2]));
7642 factor[2]= TMath::Sqrt(TMath::Abs((covarS[5] + param[2]*param[2])/covarS[5]));
7643 factor[3]= TMath::Sqrt(TMath::Abs((covarS[9] + param[3]*param[3])/covarS[9]));
7644 factor[4]= TMath::Sqrt(TMath::Abs((covarS[14] +param[4]*param[4])/covarS[14]));
7646 factor[0]=factor[2];
7647 factor[4]=factor[2];
7653 for (Int_t i=0; i<5; i++){
7654 for (Int_t j=i; j<5; j++){
7655 Int_t index=seed->GetIndex(i,j);
7656 covarS[index]*=factor[i]*factor[j];
7662 void AliTPCtrackerMI::AddCovarianceAdd(AliTPCseed * seed){
7664 // Adding systematic error - as additive factor without correlation
7666 // !!!! the systematic error for element 4 is in 1/GeV
7667 // 03.03.2012 MI changed in respect to the previous versions
7669 const Double_t *param = AliTPCReconstructor::GetRecoParam()->GetSystematicError();
7670 Double_t *covarIn= (Double_t*)seed->GetCovariance();
7672 for (Int_t i=0;i<15;i++) covar[i]=0;
7678 covar[0] = param[0]*param[0];
7679 covar[2] = param[1]*param[1];
7680 covar[5] = param[2]*param[2];
7681 covar[9] = param[3]*param[3];
7682 covar[14]= param[4]*param[4];
7684 covar[1]=TMath::Sqrt((covar[0]*covar[2]))*covarIn[1]/TMath::Sqrt((covarIn[0]*covarIn[2]));
7686 covar[3]=TMath::Sqrt((covar[0]*covar[5]))*covarIn[3]/TMath::Sqrt((covarIn[0]*covarIn[5]));
7687 covar[4]=TMath::Sqrt((covar[2]*covar[5]))*covarIn[4]/TMath::Sqrt((covarIn[2]*covarIn[5]));
7689 covar[6]=TMath::Sqrt((covar[0]*covar[9]))*covarIn[6]/TMath::Sqrt((covarIn[0]*covarIn[9]));
7690 covar[7]=TMath::Sqrt((covar[2]*covar[9]))*covarIn[7]/TMath::Sqrt((covarIn[2]*covarIn[9]));
7691 covar[8]=TMath::Sqrt((covar[5]*covar[9]))*covarIn[8]/TMath::Sqrt((covarIn[5]*covarIn[9]));
7693 covar[10]=TMath::Sqrt((covar[0]*covar[14]))*covarIn[10]/TMath::Sqrt((covarIn[0]*covarIn[14]));
7694 covar[11]=TMath::Sqrt((covar[2]*covar[14]))*covarIn[11]/TMath::Sqrt((covarIn[2]*covarIn[14]));
7695 covar[12]=TMath::Sqrt((covar[5]*covar[14]))*covarIn[12]/TMath::Sqrt((covarIn[5]*covarIn[14]));
7696 covar[13]=TMath::Sqrt((covar[9]*covar[14]))*covarIn[13]/TMath::Sqrt((covarIn[9]*covarIn[14]));
7698 seed->AddCovariance(covar);
7701 //_____________________________________________________________________________
7702 Bool_t AliTPCtrackerMI::IsTPCHVDipEvent(AliESDEvent const *esdEvent) {
7704 // check events affected by TPC HV dip
7706 if(!esdEvent) return kFALSE;
7709 if(!AliTPCcalibDB::Instance()) return kFALSE;
7710 AliTPCcalibDB::Instance()->SetRun(esdEvent->GetRunNumber());
7712 // Get HV TPC chamber sensors and calculate the median
7713 AliDCSSensorArray *voltageArray= AliTPCcalibDB::Instance()->GetVoltageSensors(esdEvent->GetRunNumber());
7714 if(!voltageArray) return kFALSE;
7716 TString sensorName="";
7717 Double_t kTPCHVdip = 2.0; // allow for 2V dip as compared to median from given sensor
7720 for(Int_t sector=0; sector<72; sector++)
7722 Char_t sideName='A';
7723 if ((sector/18)%2==1) sideName='C';
7726 sensorName=Form("TPC_ANODE_I_%c%02d_VMEAS",sideName,sector%18);
7729 sensorName=Form("TPC_ANODE_O_%c%02d_0_VMEAS",sideName,sector%18);
7732 AliDCSSensor* sensor = voltageArray->GetSensor(sensorName.Data());
7733 if(!sensor) continue;
7734 TGraph *graph = sensor->GetGraph();
7735 if(!graph) continue;
7736 Double_t median = TMath::Median(graph->GetN(), graph->GetY());
7737 if(median == 0) continue;
7739 //printf("chamber %d, sensor %s, HV %f, median %f\n", sector, sensorName.Data(), sensor->GetValue(esdEvent->GetTimeStamp()), median);
7741 if(TMath::Abs(sensor->GetValue(esdEvent->GetTimeStamp())-median)>kTPCHVdip) {
7749 //________________________________________
7750 void AliTPCtrackerMI::MarkSeedFree(TObject *sd)
7752 // account that this seed is "deleted"
7753 AliTPCseed* seed = dynamic_cast<AliTPCseed*>(sd);
7755 AliError(Form("Freeing of non-AliTPCseed %p from the pool is requested",sd));
7758 int id = seed->GetPoolID();
7760 AliError(Form("Freeing of seed %p NOT from the pool is requested",sd));
7763 // AliInfo(Form("%d %p",id, seed));
7764 fSeedsPool->RemoveAt(id);
7765 if (fFreeSeedsID.GetSize()<=fNFreeSeeds) fFreeSeedsID.Set( 2*fNFreeSeeds + 100 );
7766 fFreeSeedsID.GetArray()[fNFreeSeeds++] = id;
7769 //________________________________________
7770 TObject *&AliTPCtrackerMI::NextFreeSeed()
7772 // return next free slot where the seed can be created
7773 fLastSeedID = fNFreeSeeds ? fFreeSeedsID.GetArray()[--fNFreeSeeds] : fSeedsPool->GetEntriesFast();
7774 // AliInfo(Form("%d",fLastSeedID));
7775 return (*fSeedsPool)[ fLastSeedID ];
7779 //________________________________________
7780 void AliTPCtrackerMI::ResetSeedsPool()
7782 // mark all seeds in the pool as unused
7783 AliInfo(Form("CurrentSize: %d, BookedUpTo: %d, free: %d",fSeedsPool->GetSize(),fSeedsPool->GetEntriesFast(),fNFreeSeeds));
7785 fSeedsPool->Clear("C"); // RS: nominally the seeds may allocate memory...