1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
17 //-------------------------------------------------------
18 // Implementation of the TPC tracker
20 // Origin: Marian Ivanov Marian.Ivanov@cern.ch
22 // AliTPC parallel tracker
24 // The track fitting is based on Kalaman filtering approach
26 // The track finding steps:
27 // 1. Seeding - with and without vertex constraint
28 // - seeding with vertex constain done at first n^2 proble
29 // - seeding without vertex constraint n^3 problem
30 // 2. Tracking - follow prolongation road - find cluster - update kalman track
32 // The seeding and tracking is repeated several times, in different seeding region.
33 // This approach enables to find the track which cannot be seeded in some region of TPC
34 // This can happen because of low momenta (track do not reach outer radius), or track is currently in the ded region between sectors, or the track is for the moment overlapped with other track (seed quality is poor) ...
36 // With this approach we reach almost 100 % efficiency also for high occupancy events.
37 // (If the seeding efficiency in a region is about 90 % than with logical or of several
38 // regions we will reach 100% (in theory - supposing independence)
40 // Repeating several seeding - tracking procedures some of the tracks can be find
43 // The procedures to remove multi find tacks are impremented:
44 // RemoveUsed2 - fast procedure n problem -
45 // Algorithm - Sorting tracks according quality
46 // remove tracks with some shared fraction
47 // Sharing in respect to all tacks
48 // Signing clusters in gold region
49 // FindSplitted - slower algorithm n^2
50 // Sort the tracks according quality
51 // Loop over pair of tracks
52 // If overlap with other track bigger than threshold - remove track
54 // FindCurling - Finds the pair of tracks which are curling
55 // - About 10% of tracks can be find with this procedure
56 // The combinatorial background is too big to be used in High
57 // multiplicity environment
58 // - n^2 problem - Slow procedure - currently it is disabled because of
61 // The number of splitted tracks can be reduced disabling the sharing of the cluster.
62 // tpcRecoParam-> SetClusterSharing(kFALSE);
63 // IT IS HIGHLY non recomended to use it in high flux enviroonment
64 // Even using this switch some tracks can be found more than once
65 // (because of multiple seeding and low quality tracks which will not cross full chamber)
68 // The tracker itself can be debugged - the information about tracks can be stored in several // phases of the reconstruction
69 // To enable storage of the TPC tracks in the ESD friend track
70 // use AliTPCReconstructor::SetStreamLevel(n);
72 // The debug level - different procedure produce tree for numerical debugging
73 // To enable them set AliTPCReconstructor::SetStreamLevel(n); where nis bigger 1
77 // Adding systematic errors to the covariance:
79 // The systematic errors due to the misalignment and miscalibration are added to the covariance matrix
80 // of the tracks (not to the clusters as they are dependent):
81 // The parameters form AliTPCRecoParam are used AliTPCRecoParam::GetSystematicError
82 // The systematic errors are expressed there in RMS - position (cm), angle (rad), curvature (1/cm)
83 // The default values are 0.
85 // The sytematic errors are added to the covariance matrix in following places:
87 // 1. During fisrt itteration - AliTPCtrackerMI::FillESD
88 // 2. Second iteration -
89 // 2.a ITS->TPC - AliTPCtrackerMI::ReadSeeds
90 // 2.b TPC->TRD - AliTPCtrackerMI::PropagateBack
91 // 3. Third iteration -
92 // 3.a TRD->TPC - AliTPCtrackerMI::ReadSeeds
93 // 3.b TPC->ITS - AliTPCtrackerMI::RefitInward
95 // There are several places in the code which can be numerically debuged
96 // This code is keeped in order to enable code development and to check the calibration implementtion
98 // 1. ErrParam stream - dump information about
100 // 2.a) cluster error estimate
101 // 3.a) cluster shape estimate
104 // Debug streamer levels:
106 //-------------------------------------------------------
111 #include "Riostream.h"
112 #include <TClonesArray.h>
114 #include <TObjArray.h>
116 #include <TGraphErrors.h>
118 #include "AliComplexCluster.h"
119 #include "AliESDEvent.h"
120 #include "AliESDtrack.h"
121 #include "AliESDVertex.h"
124 #include "AliHelix.h"
125 #include "AliRunLoader.h"
126 #include "AliTPCClustersRow.h"
127 #include "AliTPCParam.h"
128 #include "AliTPCReconstructor.h"
129 #include "AliTPCpolyTrack.h"
130 #include "AliTPCreco.h"
131 #include "AliTPCseed.h"
133 #include "AliTPCtrackerSector.h"
134 #include "AliTPCtrackerMI.h"
135 #include "TStopwatch.h"
136 #include "AliTPCReconstructor.h"
137 #include "AliAlignObj.h"
138 #include "AliTrackPointArray.h"
140 #include "AliTPCcalibDB.h"
141 #include "AliTPCcalibDButil.h"
142 #include "AliTPCTransform.h"
143 #include "AliTPCClusterParam.h"
144 #include "AliTPCdEdxInfo.h"
145 #include "AliClonesPool.h"
146 #include "AliPoolsSet.h"
150 ClassImp(AliTPCtrackerMI)
154 class AliTPCFastMath {
157 static Double_t FastAsin(Double_t x);
159 static Double_t fgFastAsin[20000]; //lookup table for fast asin computation
162 Double_t AliTPCFastMath::fgFastAsin[20000];
163 AliTPCFastMath gAliTPCFastMath; // needed to fill the LUT
165 AliTPCFastMath::AliTPCFastMath(){
167 // initialized lookup table;
168 for (Int_t i=0;i<10000;i++){
169 fgFastAsin[2*i] = TMath::ASin(i/10000.);
170 fgFastAsin[2*i+1] = (TMath::ASin((i+1)/10000.)-fgFastAsin[2*i]);
174 Double_t AliTPCFastMath::FastAsin(Double_t x){
176 // return asin using lookup table
178 Int_t index = int(x*10000);
179 return fgFastAsin[2*index]+(x*10000.-index)*fgFastAsin[2*index+1];
182 Int_t index = int(x*10000);
183 return -(fgFastAsin[2*index]+(x*10000.-index)*fgFastAsin[2*index+1]);
185 //__________________________________________________________________
186 AliTPCtrackerMI::AliTPCtrackerMI()
211 // default constructor
213 for (Int_t irow=0; irow<200; irow++){
220 //_____________________________________________________________________
224 Int_t AliTPCtrackerMI::UpdateTrack(AliTPCseed * track, Int_t accept){
226 //update track information using current cluster - track->fCurrentCluster
229 AliTPCclusterMI* c =track->GetCurrentCluster();
230 if (accept > 0) //sign not accepted clusters
231 track->SetCurrentClusterIndex1(track->GetCurrentClusterIndex1() | 0x8000);
232 else // unsign accpeted clusters
233 track->SetCurrentClusterIndex1(track->GetCurrentClusterIndex1() & 0xffff7fff);
234 UInt_t i = track->GetCurrentClusterIndex1();
236 Int_t sec=(i&0xff000000)>>24;
237 //Int_t row = (i&0x00ff0000)>>16;
238 track->SetRow((i&0x00ff0000)>>16);
239 track->SetSector(sec);
240 // Int_t index = i&0xFFFF;
241 if (sec>=fkParam->GetNInnerSector()) track->SetRow(track->GetRow()+fkParam->GetNRowLow());
242 track->SetClusterIndex2(track->GetRow(), i);
243 //track->fFirstPoint = row;
244 //if ( track->fLastPoint<row) track->fLastPoint =row;
245 // if (track->fRow<0 || track->fRow>160) {
246 // printf("problem\n");
248 if (track->GetFirstPoint()>track->GetRow())
249 track->SetFirstPoint(track->GetRow());
250 if (track->GetLastPoint()<track->GetRow())
251 track->SetLastPoint(track->GetRow());
254 track->SetClusterPointer(track->GetRow(),c);
257 Double_t angle2 = track->GetSnp()*track->GetSnp();
259 //SET NEW Track Point
261 if (angle2<1) //PH sometimes angle2 is very big. To be investigated...
263 angle2 = TMath::Sqrt(angle2/(1-angle2));
264 AliTPCTrackerPoint &point =*(track->GetTrackPoint(track->GetRow()));
266 point.SetSigmaY(c->GetSigmaY2()/track->GetCurrentSigmaY2());
267 point.SetSigmaZ(c->GetSigmaZ2()/track->GetCurrentSigmaZ2());
268 point.SetErrY(sqrt(track->GetErrorY2()));
269 point.SetErrZ(sqrt(track->GetErrorZ2()));
271 point.SetX(track->GetX());
272 point.SetY(track->GetY());
273 point.SetZ(track->GetZ());
274 point.SetAngleY(angle2);
275 point.SetAngleZ(track->GetTgl());
276 if (point.IsShared()){
277 track->SetErrorY2(track->GetErrorY2()*4);
278 track->SetErrorZ2(track->GetErrorZ2()*4);
282 Double_t chi2 = track->GetPredictedChi2(track->GetCurrentCluster());
284 // track->SetErrorY2(track->GetErrorY2()*1.3);
285 // track->SetErrorY2(track->GetErrorY2()+0.01);
286 // track->SetErrorZ2(track->GetErrorZ2()*1.3);
287 // track->SetErrorZ2(track->GetErrorZ2()+0.005);
289 if (accept>0) return 0;
290 if (track->GetNumberOfClusters()%20==0){
291 // if (track->fHelixIn){
292 // TClonesArray & larr = *(track->fHelixIn);
293 // Int_t ihelix = larr.GetEntriesFast();
294 // new(larr[ihelix]) AliHelix(*track) ;
297 track->SetNoCluster(0);
298 return track->Update(c,chi2,i);
303 Int_t AliTPCtrackerMI::AcceptCluster(AliTPCseed * seed, AliTPCclusterMI * cluster)
306 // decide according desired precision to accept given
307 // cluster for tracking
309 seed->GetProlongation(cluster->GetX(),yt,zt);
310 Double_t sy2=ErrY2(seed,cluster);
311 Double_t sz2=ErrZ2(seed,cluster);
313 Double_t sdistancey2 = sy2+seed->GetSigmaY2();
314 Double_t sdistancez2 = sz2+seed->GetSigmaZ2();
315 Double_t dy=seed->GetCurrentCluster()->GetY()-yt;
316 Double_t dz=seed->GetCurrentCluster()->GetZ()-zt;
317 Double_t rdistancey2 = (seed->GetCurrentCluster()->GetY()-yt)*
318 (seed->GetCurrentCluster()->GetY()-yt)/sdistancey2;
319 Double_t rdistancez2 = (seed->GetCurrentCluster()->GetZ()-zt)*
320 (seed->GetCurrentCluster()->GetZ()-zt)/sdistancez2;
322 Double_t rdistance2 = rdistancey2+rdistancez2;
325 if (AliTPCReconstructor::StreamLevel()>2 && seed->GetNumberOfClusters()>20) {
326 Float_t rmsy2 = seed->GetCurrentSigmaY2();
327 Float_t rmsz2 = seed->GetCurrentSigmaZ2();
328 Float_t rmsy2p30 = seed->GetCMeanSigmaY2p30();
329 Float_t rmsz2p30 = seed->GetCMeanSigmaZ2p30();
330 Float_t rmsy2p30R = seed->GetCMeanSigmaY2p30R();
331 Float_t rmsz2p30R = seed->GetCMeanSigmaZ2p30R();
332 AliExternalTrackParam param(*seed);
333 static TVectorD gcl(3),gtr(3);
335 param.GetXYZ(gcl.GetMatrixArray());
336 cluster->GetGlobalXYZ(gclf);
337 gcl[0]=gclf[0]; gcl[1]=gclf[1]; gcl[2]=gclf[2];
340 if (AliTPCReconstructor::StreamLevel()>2) {
341 (*fDebugStreamer)<<"ErrParam"<<
354 "rmsy2p30="<<rmsy2p30<<
355 "rmsz2p30="<<rmsz2p30<<
356 "rmsy2p30R="<<rmsy2p30R<<
357 "rmsz2p30R="<<rmsz2p30R<<
358 // normalize distance -
359 "rdisty="<<rdistancey2<<
360 "rdistz="<<rdistancez2<<
361 "rdist="<<rdistance2<< //
365 //return 0; // temporary
366 if (rdistance2>32) return 3;
369 if ((rdistancey2>9. || rdistancez2>9.) && cluster->GetType()==0)
370 return 2; //suspisiouce - will be changed
372 if ((rdistancey2>6.25 || rdistancez2>6.25) && cluster->GetType()>0)
373 // strict cut on overlaped cluster
374 return 2; //suspisiouce - will be changed
376 if ( (rdistancey2>1. || rdistancez2>6.25 )
377 && cluster->GetType()<0){
378 seed->SetNFoundable(seed->GetNFoundable()-1);
382 if (fUseHLTClusters == 3 || fUseHLTClusters == 4) {
384 if(!AliTPCReconstructor::GetRecoParam()->GetUseHLTOnePadCluster()) {
385 if (TMath::Abs(cluster->GetSigmaY2()) < kAlmost0)
398 //_____________________________________________________________________________
399 AliTPCtrackerMI::AliTPCtrackerMI(const AliTPCParam *par):
401 fkNIS(par->GetNInnerSector()/2),
403 fkNOS(par->GetNOuterSector()/2),
423 //---------------------------------------------------------------------
424 // The main TPC tracker constructor
425 //---------------------------------------------------------------------
426 fInnerSec=new AliTPCtrackerSector[fkNIS];
427 fOuterSec=new AliTPCtrackerSector[fkNOS];
430 for (i=0; i<fkNIS; i++) fInnerSec[i].Setup(par,0);
431 for (i=0; i<fkNOS; i++) fOuterSec[i].Setup(par,1);
434 Int_t nrowlow = par->GetNRowLow();
435 Int_t nrowup = par->GetNRowUp();
438 for (i=0;i<nrowlow;i++){
439 fXRow[i] = par->GetPadRowRadiiLow(i);
440 fPadLength[i]= par->GetPadPitchLength(0,i);
441 fYMax[i] = fXRow[i]*TMath::Tan(0.5*par->GetInnerAngle());
445 for (i=0;i<nrowup;i++){
446 fXRow[i+nrowlow] = par->GetPadRowRadiiUp(i);
447 fPadLength[i+nrowlow] = par->GetPadPitchLength(60,i);
448 fYMax[i+nrowlow] = fXRow[i+nrowlow]*TMath::Tan(0.5*par->GetOuterAngle());
451 if (AliTPCReconstructor::StreamLevel()>0) {
452 fDebugStreamer = new TTreeSRedirector("TPCdebug.root");
455 fSeedsPool = new AliClonesPool("AliTPCseed",1000);
456 fSeedsPool->SetName("TPCSeed");
457 fKinksPool = new AliClonesPool("AliKink",1000);
458 fKinksPool->SetName("TPCKink");
461 //________________________________________________________________________
462 AliTPCtrackerMI::AliTPCtrackerMI(const AliTPCtrackerMI &t):
486 //------------------------------------
487 // dummy copy constructor
488 //------------------------------------------------------------------
490 for (Int_t irow=0; irow<200; irow++){
497 AliTPCtrackerMI & AliTPCtrackerMI::operator=(const AliTPCtrackerMI& /*r*/)
499 //------------------------------
501 //--------------------------------------------------------------
504 //_____________________________________________________________________________
505 AliTPCtrackerMI::~AliTPCtrackerMI() {
506 //------------------------------------------------------------------
507 // TPC tracker destructor
508 //------------------------------------------------------------------
515 if (fDebugStreamer) delete fDebugStreamer;
516 if (fSeedsPool && !fPools) delete fSeedsPool; // RS: owned by reconstruction pools (if exist)
517 if (fKinksPool && !fPools) delete fKinksPool; // RS: owned by reconstruction pools (if exist)
521 void AliTPCtrackerMI::FillESD(const TObjArray* arr)
525 //fill esds using updated tracks
528 // write tracks to the event
529 // store index of the track
530 Int_t nseed=arr->GetEntriesFast();
531 //FindKinks(arr,fEvent);
532 for (Int_t i=0; i<nseed; i++) {
533 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
537 if (AliTPCReconstructor::StreamLevel()>1) {
538 (*fDebugStreamer)<<"Track0"<<
542 // pt->PropagateTo(fkParam->GetInnerRadiusLow());
543 if (pt->GetKinkIndex(0)<=0){ //don't propagate daughter tracks
544 pt->PropagateTo(fkParam->GetInnerRadiusLow());
547 if (( pt->GetPoints()[2]- pt->GetPoints()[0])>5 && pt->GetPoints()[3]>0.8){
549 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
550 iotrack.SetTPCPoints(pt->GetPoints());
551 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
552 iotrack.SetV0Indexes(pt->GetV0Indexes());
553 // iotrack.SetTPCpid(pt->fTPCr);
554 //iotrack.SetTPCindex(i);
555 fEvent->AddTrack(&iotrack);
559 if ( (pt->GetNumberOfClusters()>70)&& (Float_t(pt->GetNumberOfClusters())/Float_t(pt->GetNFoundable()))>0.55) {
561 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
562 iotrack.SetTPCPoints(pt->GetPoints());
563 //iotrack.SetTPCindex(i);
564 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
565 iotrack.SetV0Indexes(pt->GetV0Indexes());
566 // iotrack.SetTPCpid(pt->fTPCr);
567 fEvent->AddTrack(&iotrack);
571 // short tracks - maybe decays
573 if ( (pt->GetNumberOfClusters()>30) && (Float_t(pt->GetNumberOfClusters())/Float_t(pt->GetNFoundable()))>0.70) {
574 Int_t found,foundable,shared;
575 pt->GetClusterStatistic(0,60,found, foundable,shared,kFALSE);
576 if ( (found>20) && (pt->GetNShared()/float(pt->GetNumberOfClusters())<0.2)){
578 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
579 //iotrack.SetTPCindex(i);
580 iotrack.SetTPCPoints(pt->GetPoints());
581 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
582 iotrack.SetV0Indexes(pt->GetV0Indexes());
583 //iotrack.SetTPCpid(pt->fTPCr);
584 fEvent->AddTrack(&iotrack);
589 if ( (pt->GetNumberOfClusters()>20) && (Float_t(pt->GetNumberOfClusters())/Float_t(pt->GetNFoundable()))>0.8) {
590 Int_t found,foundable,shared;
591 pt->GetClusterStatistic(0,60,found, foundable,shared,kFALSE);
592 if (found<20) continue;
593 if (pt->GetNShared()/float(pt->GetNumberOfClusters())>0.2) continue;
596 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
597 iotrack.SetTPCPoints(pt->GetPoints());
598 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
599 iotrack.SetV0Indexes(pt->GetV0Indexes());
600 //iotrack.SetTPCpid(pt->fTPCr);
601 //iotrack.SetTPCindex(i);
602 fEvent->AddTrack(&iotrack);
605 // short tracks - secondaties
607 if ( (pt->GetNumberOfClusters()>30) ) {
608 Int_t found,foundable,shared;
609 pt->GetClusterStatistic(128,158,found, foundable,shared,kFALSE);
610 if ( (found>20) && (pt->GetNShared()/float(pt->GetNumberOfClusters())<0.2) &&float(found)/float(foundable)>0.8){
612 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
613 iotrack.SetTPCPoints(pt->GetPoints());
614 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
615 iotrack.SetV0Indexes(pt->GetV0Indexes());
616 //iotrack.SetTPCpid(pt->fTPCr);
617 //iotrack.SetTPCindex(i);
618 fEvent->AddTrack(&iotrack);
623 if ( (pt->GetNumberOfClusters()>15)) {
624 Int_t found,foundable,shared;
625 pt->GetClusterStatistic(138,158,found, foundable,shared,kFALSE);
626 if (found<15) continue;
627 if (foundable<=0) continue;
628 if (pt->GetNShared()/float(pt->GetNumberOfClusters())>0.2) continue;
629 if (float(found)/float(foundable)<0.8) continue;
632 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
633 iotrack.SetTPCPoints(pt->GetPoints());
634 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
635 iotrack.SetV0Indexes(pt->GetV0Indexes());
636 // iotrack.SetTPCpid(pt->fTPCr);
637 //iotrack.SetTPCindex(i);
638 fEvent->AddTrack(&iotrack);
642 // >> account for suppressed tracks in the kink indices (RS)
643 int nESDtracks = fEvent->GetNumberOfTracks();
644 for (int it=nESDtracks;it--;) {
645 AliESDtrack* esdTr = fEvent->GetTrack(it);
646 if (!esdTr || !esdTr->GetKinkIndex(0)) continue;
647 for (int ik=0;ik<3;ik++) {
649 if (!(knkId=esdTr->GetKinkIndex(ik))) break; // no more kinks for this track
650 AliESDkink* kink = fEvent->GetKink(TMath::Abs(knkId)-1);
652 AliError(Form("ESDTrack%d refers to non-existing kink %d",it,TMath::Abs(knkId)-1));
655 kink->SetIndex(it, knkId<0 ? 0:1); // update track index of the kink: mother at 0, daughter at 1
658 // << account for suppressed tracks in the kink indices (RS)
659 AliInfo(Form("Number of filled ESDs-\t%d\n",fEvent->GetNumberOfTracks()));
667 Double_t AliTPCtrackerMI::ErrY2(AliTPCseed* seed, const AliTPCclusterMI * cl){
670 // Use calibrated cluster error from OCDB
672 AliTPCClusterParam * clparam = AliTPCcalibDB::Instance()->GetClusterParam();
674 Float_t z = TMath::Abs(fkParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
675 Int_t ctype = cl->GetType();
676 Int_t type = (cl->GetRow()<63) ? 0: (cl->GetRow()>126) ? 1:2;
677 Double_t angle = seed->GetSnp()*seed->GetSnp();
678 angle = TMath::Sqrt(TMath::Abs(angle/(1.-angle)));
679 Double_t erry2 = clparam->GetError0Par(0,type, z,angle);
681 erry2+=0.5; // edge cluster
684 seed->SetErrorY2(erry2);
688 //calculate look-up table at the beginning
689 // static Bool_t ginit = kFALSE;
690 // static Float_t gnoise1,gnoise2,gnoise3;
691 // static Float_t ggg1[10000];
692 // static Float_t ggg2[10000];
693 // static Float_t ggg3[10000];
694 // static Float_t glandau1[10000];
695 // static Float_t glandau2[10000];
696 // static Float_t glandau3[10000];
698 // static Float_t gcor01[500];
699 // static Float_t gcor02[500];
700 // static Float_t gcorp[500];
704 // if (ginit==kFALSE){
705 // for (Int_t i=1;i<500;i++){
706 // Float_t rsigma = float(i)/100.;
707 // gcor02[i] = TMath::Max(0.78 +TMath::Exp(7.4*(rsigma-1.2)),0.6);
708 // gcor01[i] = TMath::Max(0.72 +TMath::Exp(3.36*(rsigma-1.2)),0.6);
709 // gcorp[i] = TMath::Max(TMath::Power((rsigma+0.5),1.5),1.2);
713 // for (Int_t i=3;i<10000;i++){
717 // Float_t amp = float(i);
718 // Float_t padlength =0.75;
719 // gnoise1 = 0.0004/padlength;
720 // Float_t nel = 0.268*amp;
721 // Float_t nprim = 0.155*amp;
722 // ggg1[i] = fkParam->GetDiffT()*fkParam->GetDiffT()*(2+0.001*nel/(padlength*padlength))/nel;
723 // glandau1[i] = (2.+0.12*nprim)*0.5* (2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
724 // if (glandau1[i]>1) glandau1[i]=1;
725 // glandau1[i]*=padlength*padlength/12.;
729 // gnoise2 = 0.0004/padlength;
731 // nprim = 0.133*amp;
732 // ggg2[i] = fkParam->GetDiffT()*fkParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
733 // glandau2[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
734 // if (glandau2[i]>1) glandau2[i]=1;
735 // glandau2[i]*=padlength*padlength/12.;
740 // gnoise3 = 0.0004/padlength;
742 // nprim = 0.133*amp;
743 // ggg3[i] = fkParam->GetDiffT()*fkParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
744 // glandau3[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
745 // if (glandau3[i]>1) glandau3[i]=1;
746 // glandau3[i]*=padlength*padlength/12.;
754 // Int_t amp = int(TMath::Abs(cl->GetQ()));
756 // seed->SetErrorY2(1.);
760 // Float_t z = TMath::Abs(fkParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
761 // Int_t ctype = cl->GetType();
762 // Float_t padlength= GetPadPitchLength(seed->GetRow());
763 // Double_t angle2 = seed->GetSnp()*seed->GetSnp();
764 // angle2 = angle2/(1-angle2);
766 // //cluster "quality"
767 // Int_t rsigmay = int(100.*cl->GetSigmaY2()/(seed->GetCurrentSigmaY2()));
770 // if (fSectors==fInnerSec){
771 // snoise2 = gnoise1;
772 // res = ggg1[amp]*z+glandau1[amp]*angle2;
773 // if (ctype==0) res *= gcor01[rsigmay];
776 // res*= gcorp[rsigmay];
780 // if (padlength<1.1){
781 // snoise2 = gnoise2;
782 // res = ggg2[amp]*z+glandau2[amp]*angle2;
783 // if (ctype==0) res *= gcor02[rsigmay];
786 // res*= gcorp[rsigmay];
790 // snoise2 = gnoise3;
791 // res = ggg3[amp]*z+glandau3[amp]*angle2;
792 // if (ctype==0) res *= gcor02[rsigmay];
795 // res*= gcorp[rsigmay];
802 // res*=2.4; // overestimate error 2 times
806 // if (res<2*snoise2)
809 // seed->SetErrorY2(res);
817 Double_t AliTPCtrackerMI::ErrZ2(AliTPCseed* seed, const AliTPCclusterMI * cl){
820 // Use calibrated cluster error from OCDB
822 AliTPCClusterParam * clparam = AliTPCcalibDB::Instance()->GetClusterParam();
824 Float_t z = TMath::Abs(fkParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
825 Int_t ctype = cl->GetType();
826 Int_t type = (cl->GetRow()<63) ? 0: (cl->GetRow()>126) ? 1:2;
828 Double_t angle2 = seed->GetSnp()*seed->GetSnp();
829 angle2 = seed->GetTgl()*seed->GetTgl()*(1+angle2/(1-angle2));
830 Double_t angle = TMath::Sqrt(TMath::Abs(angle2));
831 Double_t errz2 = clparam->GetError0Par(1,type, z,angle);
833 errz2+=0.5; // edge cluster
836 seed->SetErrorZ2(errz2);
842 // //seed->SetErrorY2(0.1);
844 // //calculate look-up table at the beginning
845 // static Bool_t ginit = kFALSE;
846 // static Float_t gnoise1,gnoise2,gnoise3;
847 // static Float_t ggg1[10000];
848 // static Float_t ggg2[10000];
849 // static Float_t ggg3[10000];
850 // static Float_t glandau1[10000];
851 // static Float_t glandau2[10000];
852 // static Float_t glandau3[10000];
854 // static Float_t gcor01[1000];
855 // static Float_t gcor02[1000];
856 // static Float_t gcorp[1000];
860 // if (ginit==kFALSE){
861 // for (Int_t i=1;i<1000;i++){
862 // Float_t rsigma = float(i)/100.;
863 // gcor02[i] = TMath::Max(0.81 +TMath::Exp(6.8*(rsigma-1.2)),0.6);
864 // gcor01[i] = TMath::Max(0.72 +TMath::Exp(2.04*(rsigma-1.2)),0.6);
865 // gcorp[i] = TMath::Max(TMath::Power((rsigma+0.5),1.5),1.2);
869 // for (Int_t i=3;i<10000;i++){
873 // Float_t amp = float(i);
874 // Float_t padlength =0.75;
875 // gnoise1 = 0.0004/padlength;
876 // Float_t nel = 0.268*amp;
877 // Float_t nprim = 0.155*amp;
878 // ggg1[i] = fkParam->GetDiffT()*fkParam->GetDiffT()*(2+0.001*nel/(padlength*padlength))/nel;
879 // glandau1[i] = (2.+0.12*nprim)*0.5* (2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
880 // if (glandau1[i]>1) glandau1[i]=1;
881 // glandau1[i]*=padlength*padlength/12.;
885 // gnoise2 = 0.0004/padlength;
887 // nprim = 0.133*amp;
888 // ggg2[i] = fkParam->GetDiffT()*fkParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
889 // glandau2[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
890 // if (glandau2[i]>1) glandau2[i]=1;
891 // glandau2[i]*=padlength*padlength/12.;
896 // gnoise3 = 0.0004/padlength;
898 // nprim = 0.133*amp;
899 // ggg3[i] = fkParam->GetDiffT()*fkParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
900 // glandau3[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
901 // if (glandau3[i]>1) glandau3[i]=1;
902 // glandau3[i]*=padlength*padlength/12.;
910 // Int_t amp = int(TMath::Abs(cl->GetQ()));
912 // seed->SetErrorY2(1.);
916 // Float_t z = TMath::Abs(fkParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
917 // Int_t ctype = cl->GetType();
918 // Float_t padlength= GetPadPitchLength(seed->GetRow());
920 // Double_t angle2 = seed->GetSnp()*seed->GetSnp();
921 // // if (angle2<0.6) angle2 = 0.6;
922 // angle2 = seed->GetTgl()*seed->GetTgl()*(1+angle2/(1-angle2));
924 // //cluster "quality"
925 // Int_t rsigmaz = int(100.*cl->GetSigmaZ2()/(seed->GetCurrentSigmaZ2()));
928 // if (fSectors==fInnerSec){
929 // snoise2 = gnoise1;
930 // res = ggg1[amp]*z+glandau1[amp]*angle2;
931 // if (ctype==0) res *= gcor01[rsigmaz];
934 // res*= gcorp[rsigmaz];
938 // if (padlength<1.1){
939 // snoise2 = gnoise2;
940 // res = ggg2[amp]*z+glandau2[amp]*angle2;
941 // if (ctype==0) res *= gcor02[rsigmaz];
944 // res*= gcorp[rsigmaz];
948 // snoise2 = gnoise3;
949 // res = ggg3[amp]*z+glandau3[amp]*angle2;
950 // if (ctype==0) res *= gcor02[rsigmaz];
953 // res*= gcorp[rsigmaz];
962 // if ((ctype<0) &&<70){
967 // if (res<2*snoise2)
969 // if (res>3) res =3;
970 // seed->SetErrorZ2(res);
978 void AliTPCtrackerMI::RotateToLocal(AliTPCseed *seed)
980 //rotate to track "local coordinata
981 Float_t x = seed->GetX();
982 Float_t y = seed->GetY();
983 Float_t ymax = x*TMath::Tan(0.5*fSectors->GetAlpha());
986 seed->SetRelativeSector((seed->GetRelativeSector()+1) % fN);
987 if (!seed->Rotate(fSectors->GetAlpha()))
989 } else if (y <-ymax) {
990 seed->SetRelativeSector((seed->GetRelativeSector()-1+fN) % fN);
991 if (!seed->Rotate(-fSectors->GetAlpha()))
999 //_____________________________________________________________________________
1000 Double_t AliTPCtrackerMI::F1old(Double_t x1,Double_t y1,
1001 Double_t x2,Double_t y2,
1002 Double_t x3,Double_t y3) const
1004 //-----------------------------------------------------------------
1005 // Initial approximation of the track curvature
1006 //-----------------------------------------------------------------
1007 Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
1008 Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
1009 (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
1010 Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
1011 (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
1013 Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
1014 if ( xr*xr+yr*yr<=0.00000000000001) return 100;
1015 return -xr*yr/sqrt(xr*xr+yr*yr);
1020 //_____________________________________________________________________________
1021 Double_t AliTPCtrackerMI::F1(Double_t x1,Double_t y1,
1022 Double_t x2,Double_t y2,
1023 Double_t x3,Double_t y3) const
1025 //-----------------------------------------------------------------
1026 // Initial approximation of the track curvature
1027 //-----------------------------------------------------------------
1033 Double_t det = x3*y2-x2*y3;
1034 if (TMath::Abs(det)<1e-10){
1038 Double_t u = 0.5* (x2*(x2-x3)+y2*(y2-y3))/det;
1039 Double_t x0 = x3*0.5-y3*u;
1040 Double_t y0 = y3*0.5+x3*u;
1041 Double_t c2 = 1/TMath::Sqrt(x0*x0+y0*y0);
1047 Double_t AliTPCtrackerMI::F2(Double_t x1,Double_t y1,
1048 Double_t x2,Double_t y2,
1049 Double_t x3,Double_t y3) const
1051 //-----------------------------------------------------------------
1052 // Initial approximation of the track curvature
1053 //-----------------------------------------------------------------
1059 Double_t det = x3*y2-x2*y3;
1060 if (TMath::Abs(det)<1e-10) {
1064 Double_t u = 0.5* (x2*(x2-x3)+y2*(y2-y3))/det;
1065 Double_t x0 = x3*0.5-y3*u;
1066 Double_t y0 = y3*0.5+x3*u;
1067 Double_t c2 = 1/TMath::Sqrt(x0*x0+y0*y0);
1076 //_____________________________________________________________________________
1077 Double_t AliTPCtrackerMI::F2old(Double_t x1,Double_t y1,
1078 Double_t x2,Double_t y2,
1079 Double_t x3,Double_t y3) const
1081 //-----------------------------------------------------------------
1082 // Initial approximation of the track curvature times center of curvature
1083 //-----------------------------------------------------------------
1084 Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
1085 Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
1086 (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
1087 Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
1088 (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
1090 Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
1092 return -a/(d*y1-b)*xr/sqrt(xr*xr+yr*yr);
1095 //_____________________________________________________________________________
1096 Double_t AliTPCtrackerMI::F3(Double_t x1,Double_t y1,
1097 Double_t x2,Double_t y2,
1098 Double_t z1,Double_t z2) const
1100 //-----------------------------------------------------------------
1101 // Initial approximation of the tangent of the track dip angle
1102 //-----------------------------------------------------------------
1103 return (z1 - z2)/sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
1107 Double_t AliTPCtrackerMI::F3n(Double_t x1,Double_t y1,
1108 Double_t x2,Double_t y2,
1109 Double_t z1,Double_t z2, Double_t c) const
1111 //-----------------------------------------------------------------
1112 // Initial approximation of the tangent of the track dip angle
1113 //-----------------------------------------------------------------
1117 //angle1 = (z1-z2)*c/(TMath::ASin(c*x1-ni)-TMath::ASin(c*x2-ni));
1119 Double_t d = TMath::Sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
1120 if (TMath::Abs(d*c*0.5)>1) return 0;
1121 // Double_t angle2 = TMath::ASin(d*c*0.5);
1122 // Double_t angle2 = AliTPCFastMath::FastAsin(d*c*0.5);
1123 Double_t angle2 = (d*c*0.5>0.1)? TMath::ASin(d*c*0.5): AliTPCFastMath::FastAsin(d*c*0.5);
1125 angle2 = (z1-z2)*c/(angle2*2.);
1129 Bool_t AliTPCtrackerMI::GetProlongation(Double_t x1, Double_t x2, Double_t x[5], Double_t &y, Double_t &z) const
1130 {//-----------------------------------------------------------------
1131 // This function find proloncation of a track to a reference plane x=x2.
1132 //-----------------------------------------------------------------
1136 if (TMath::Abs(x[4]*x1 - x[2]) >= 0.999) {
1140 Double_t c1=x[4]*x1 - x[2], r1=TMath::Sqrt((1.-c1)*(1.+c1));
1141 Double_t c2=x[4]*x2 - x[2], r2=TMath::Sqrt((1.-c2)*(1.+c2));
1145 Double_t dy = dx*(c1+c2)/(r1+r2);
1148 Double_t delta = x[4]*dx*(c1+c2)/(c1*r2 + c2*r1);
1150 if (TMath::Abs(delta)>0.01){
1151 dz = x[3]*TMath::ASin(delta)/x[4];
1153 dz = x[3]*AliTPCFastMath::FastAsin(delta)/x[4];
1156 //dz = x[3]*AliTPCFastMath::FastAsin(delta)/x[4];
1164 Int_t AliTPCtrackerMI::LoadClusters (TTree *const tree)
1169 return LoadClusters();
1173 Int_t AliTPCtrackerMI::LoadClusters(const TObjArray *arr)
1176 // load clusters to the memory
1177 AliTPCClustersRow *clrow = new AliTPCClustersRow("AliTPCclusterMI");
1178 Int_t lower = arr->LowerBound();
1179 Int_t entries = arr->GetEntriesFast();
1181 for (Int_t i=lower; i<entries; i++) {
1182 clrow = (AliTPCClustersRow*) arr->At(i);
1183 if(!clrow) continue;
1184 if(!clrow->GetArray()) continue;
1188 fkParam->AdjustSectorRow(clrow->GetID(),sec,row);
1190 for (Int_t icl=0; icl<clrow->GetArray()->GetEntriesFast(); icl++){
1191 Transform((AliTPCclusterMI*)(clrow->GetArray()->At(icl)));
1194 if (clrow->GetArray()->GetEntriesFast()<=0) continue;
1195 AliTPCtrackerRow * tpcrow=0;
1198 tpcrow = &(fInnerSec[sec%fkNIS][row]);
1202 tpcrow = &(fOuterSec[(sec-fkNIS*2)%fkNOS][row]);
1203 left = (sec-fkNIS*2)/fkNOS;
1206 tpcrow->SetN1(clrow->GetArray()->GetEntriesFast());
1207 for (Int_t j=0;j<tpcrow->GetN1();++j)
1208 tpcrow->SetCluster1(j, *(AliTPCclusterMI*)(clrow->GetArray()->At(j)));
1211 tpcrow->SetN2(clrow->GetArray()->GetEntriesFast());
1212 for (Int_t j=0;j<tpcrow->GetN2();++j)
1213 tpcrow->SetCluster2(j, *(AliTPCclusterMI*)(clrow->GetArray()->At(j)));
1215 clrow->GetArray()->Clear("C");
1224 Int_t AliTPCtrackerMI::LoadClusters(const TClonesArray *arr)
1227 // load clusters to the memory from one
1230 AliTPCclusterMI *clust=0;
1231 Int_t count[72][96] = { {0} , {0} };
1233 // loop over clusters
1234 for (Int_t icl=0; icl<arr->GetEntriesFast(); icl++) {
1235 clust = (AliTPCclusterMI*)arr->At(icl);
1236 if(!clust) continue;
1237 //printf("cluster: det %d, row %d \n", clust->GetDetector(),clust->GetRow());
1239 // transform clusters
1242 // count clusters per pad row
1243 count[clust->GetDetector()][clust->GetRow()]++;
1246 // insert clusters to sectors
1247 for (Int_t icl=0; icl<arr->GetEntriesFast(); icl++) {
1248 clust = (AliTPCclusterMI*)arr->At(icl);
1249 if(!clust) continue;
1251 Int_t sec = clust->GetDetector();
1252 Int_t row = clust->GetRow();
1254 // filter overlapping pad rows needed by HLT
1255 if(sec<fkNIS*2) { //IROCs
1256 if(row == 30) continue;
1259 if(row == 27 || row == 76) continue;
1265 fInnerSec[sec%fkNIS].InsertCluster(clust, count[sec][row], fkParam);
1268 left = (sec-fkNIS*2)/fkNOS;
1269 fOuterSec[(sec-fkNIS*2)%fkNOS].InsertCluster(clust, count[sec][row], fkParam);
1273 // Load functions must be called behind LoadCluster(TClonesArray*)
1275 //LoadOuterSectors();
1276 //LoadInnerSectors();
1282 Int_t AliTPCtrackerMI::LoadClusters()
1285 // load clusters to the memory
1286 static AliTPCClustersRow *clrow= new AliTPCClustersRow("AliTPCclusterMI");
1288 // TTree * tree = fClustersArray.GetTree();
1290 TTree * tree = fInput;
1291 TBranch * br = tree->GetBranch("Segment");
1292 br->SetAddress(&clrow);
1294 // Conversion of pad, row coordinates in local tracking coords.
1295 // Could be skipped here; is already done in clusterfinder
1297 Int_t j=Int_t(tree->GetEntries());
1298 for (Int_t i=0; i<j; i++) {
1302 fkParam->AdjustSectorRow(clrow->GetID(),sec,row);
1303 for (Int_t icl=0; icl<clrow->GetArray()->GetEntriesFast(); icl++){
1304 Transform((AliTPCclusterMI*)(clrow->GetArray()->At(icl)));
1307 AliTPCtrackerRow * tpcrow=0;
1310 tpcrow = &(fInnerSec[sec%fkNIS][row]);
1314 tpcrow = &(fOuterSec[(sec-fkNIS*2)%fkNOS][row]);
1315 left = (sec-fkNIS*2)/fkNOS;
1318 tpcrow->SetN1(clrow->GetArray()->GetEntriesFast());
1319 for (Int_t k=0;k<tpcrow->GetN1();++k)
1320 tpcrow->SetCluster1(k, *(AliTPCclusterMI*)(clrow->GetArray()->At(k)));
1323 tpcrow->SetN2(clrow->GetArray()->GetEntriesFast());
1324 for (Int_t k=0;k<tpcrow->GetN2();++k)
1325 tpcrow->SetCluster2(k, *(AliTPCclusterMI*)(clrow->GetArray()->At(k)));
1336 void AliTPCtrackerMI::UnloadClusters()
1339 // unload clusters from the memory
1341 Int_t nrows = fOuterSec->GetNRows();
1342 for (Int_t sec = 0;sec<fkNOS;sec++)
1343 for (Int_t row = 0;row<nrows;row++){
1344 AliTPCtrackerRow* tpcrow = &(fOuterSec[sec%fkNOS][row]);
1346 // if (tpcrow->fClusters1) delete []tpcrow->fClusters1;
1347 // if (tpcrow->fClusters2) delete []tpcrow->fClusters2;
1349 tpcrow->ResetClusters();
1352 nrows = fInnerSec->GetNRows();
1353 for (Int_t sec = 0;sec<fkNIS;sec++)
1354 for (Int_t row = 0;row<nrows;row++){
1355 AliTPCtrackerRow* tpcrow = &(fInnerSec[sec%fkNIS][row]);
1357 // if (tpcrow->fClusters1) delete []tpcrow->fClusters1;
1358 //if (tpcrow->fClusters2) delete []tpcrow->fClusters2;
1360 tpcrow->ResetClusters();
1366 void AliTPCtrackerMI::FillClusterArray(TObjArray* array) const{
1368 // Filling cluster to the array - For visualization purposes
1371 nrows = fOuterSec->GetNRows();
1372 for (Int_t sec = 0;sec<fkNOS;sec++)
1373 for (Int_t row = 0;row<nrows;row++){
1374 AliTPCtrackerRow* tpcrow = &(fOuterSec[sec%fkNOS][row]);
1375 if (!tpcrow) continue;
1376 for (Int_t icl = 0;icl<tpcrow->GetN();icl++){
1377 array->AddLast((TObject*)((*tpcrow)[icl]));
1380 nrows = fInnerSec->GetNRows();
1381 for (Int_t sec = 0;sec<fkNIS;sec++)
1382 for (Int_t row = 0;row<nrows;row++){
1383 AliTPCtrackerRow* tpcrow = &(fInnerSec[sec%fkNIS][row]);
1384 if (!tpcrow) continue;
1385 for (Int_t icl = 0;icl<tpcrow->GetN();icl++){
1386 array->AddLast((TObject*)(*tpcrow)[icl]);
1392 void AliTPCtrackerMI::Transform(AliTPCclusterMI * cluster){
1396 AliTPCcalibDB * calibDB = AliTPCcalibDB::Instance();
1397 AliTPCTransform *transform = calibDB->GetTransform() ;
1399 AliFatal("Tranformations not in calibDB");
1402 transform->SetCurrentRecoParam((AliTPCRecoParam*)AliTPCReconstructor::GetRecoParam());
1403 Double_t x[3]={cluster->GetRow(),cluster->GetPad(),cluster->GetTimeBin()};
1404 Int_t i[1]={cluster->GetDetector()};
1405 transform->Transform(x,i,0,1);
1406 // if (cluster->GetDetector()%36>17){
1411 // in debug mode check the transformation
1413 if (AliTPCReconstructor::StreamLevel()>2) {
1415 cluster->GetGlobalXYZ(gx);
1416 Int_t event = (fEvent==NULL)? 0: fEvent->GetEventNumberInFile();
1417 TTreeSRedirector &cstream = *fDebugStreamer;
1418 cstream<<"Transform"<<
1429 cluster->SetX(x[0]);
1430 cluster->SetY(x[1]);
1431 cluster->SetZ(x[2]);
1436 //if (!fkParam->IsGeoRead()) fkParam->ReadGeoMatrices();
1437 if (AliTPCReconstructor::GetRecoParam()->GetUseSectorAlignment() && (!calibDB->HasAlignmentOCDB())){
1438 TGeoHMatrix *mat = fkParam->GetClusterMatrix(cluster->GetDetector());
1440 Double_t pos[3]= {cluster->GetX(),cluster->GetY(),cluster->GetZ()};
1441 Double_t posC[3]={cluster->GetX(),cluster->GetY(),cluster->GetZ()};
1442 if (mat) mat->LocalToMaster(pos,posC);
1444 // chack Loading of Geo matrices from GeoManager - TEMPORARY FIX
1446 cluster->SetX(posC[0]);
1447 cluster->SetY(posC[1]);
1448 cluster->SetZ(posC[2]);
1452 //_____________________________________________________________________________
1453 Int_t AliTPCtrackerMI::LoadOuterSectors() {
1454 //-----------------------------------------------------------------
1455 // This function fills outer TPC sectors with clusters.
1456 //-----------------------------------------------------------------
1457 Int_t nrows = fOuterSec->GetNRows();
1459 for (Int_t sec = 0;sec<fkNOS;sec++)
1460 for (Int_t row = 0;row<nrows;row++){
1461 AliTPCtrackerRow* tpcrow = &(fOuterSec[sec%fkNOS][row]);
1462 Int_t sec2 = sec+2*fkNIS;
1464 Int_t ncl = tpcrow->GetN1();
1466 AliTPCclusterMI *c= (tpcrow->GetCluster1(ncl));
1467 index=(((sec2<<8)+row)<<16)+ncl;
1468 tpcrow->InsertCluster(c,index);
1471 ncl = tpcrow->GetN2();
1473 AliTPCclusterMI *c= (tpcrow->GetCluster2(ncl));
1474 index=((((sec2+fkNOS)<<8)+row)<<16)+ncl;
1475 tpcrow->InsertCluster(c,index);
1478 // write indexes for fast acces
1480 for (Int_t i=0;i<510;i++)
1481 tpcrow->SetFastCluster(i,-1);
1482 for (Int_t i=0;i<tpcrow->GetN();i++){
1483 Int_t zi = Int_t((*tpcrow)[i]->GetZ()+255.);
1484 tpcrow->SetFastCluster(zi,i); // write index
1487 for (Int_t i=0;i<510;i++){
1488 if (tpcrow->GetFastCluster(i)<0)
1489 tpcrow->SetFastCluster(i,last);
1491 last = tpcrow->GetFastCluster(i);
1500 //_____________________________________________________________________________
1501 Int_t AliTPCtrackerMI::LoadInnerSectors() {
1502 //-----------------------------------------------------------------
1503 // This function fills inner TPC sectors with clusters.
1504 //-----------------------------------------------------------------
1505 Int_t nrows = fInnerSec->GetNRows();
1507 for (Int_t sec = 0;sec<fkNIS;sec++)
1508 for (Int_t row = 0;row<nrows;row++){
1509 AliTPCtrackerRow* tpcrow = &(fInnerSec[sec%fkNIS][row]);
1512 Int_t ncl = tpcrow->GetN1();
1514 AliTPCclusterMI *c= (tpcrow->GetCluster1(ncl));
1515 index=(((sec<<8)+row)<<16)+ncl;
1516 tpcrow->InsertCluster(c,index);
1519 ncl = tpcrow->GetN2();
1521 AliTPCclusterMI *c= (tpcrow->GetCluster2(ncl));
1522 index=((((sec+fkNIS)<<8)+row)<<16)+ncl;
1523 tpcrow->InsertCluster(c,index);
1526 // write indexes for fast acces
1528 for (Int_t i=0;i<510;i++)
1529 tpcrow->SetFastCluster(i,-1);
1530 for (Int_t i=0;i<tpcrow->GetN();i++){
1531 Int_t zi = Int_t((*tpcrow)[i]->GetZ()+255.);
1532 tpcrow->SetFastCluster(zi,i); // write index
1535 for (Int_t i=0;i<510;i++){
1536 if (tpcrow->GetFastCluster(i)<0)
1537 tpcrow->SetFastCluster(i,last);
1539 last = tpcrow->GetFastCluster(i);
1551 //_________________________________________________________________________
1552 AliTPCclusterMI *AliTPCtrackerMI::GetClusterMI(Int_t index) const {
1553 //--------------------------------------------------------------------
1554 // Return pointer to a given cluster
1555 //--------------------------------------------------------------------
1556 if (index<0) return 0; // no cluster
1557 Int_t sec=(index&0xff000000)>>24;
1558 Int_t row=(index&0x00ff0000)>>16;
1559 Int_t ncl=(index&0x00007fff)>>00;
1561 const AliTPCtrackerRow * tpcrow=0;
1562 TClonesArray * clrow =0;
1564 if (sec<0 || sec>=fkNIS*4) {
1565 AliWarning(Form("Wrong sector %d",sec));
1570 AliTPCtrackerSector& tracksec = fInnerSec[sec%fkNIS];
1571 if (tracksec.GetNRows()<=row) return 0;
1572 tpcrow = &(tracksec[row]);
1573 if (tpcrow==0) return 0;
1576 if (tpcrow->GetN1()<=ncl) return 0;
1577 clrow = tpcrow->GetClusters1();
1580 if (tpcrow->GetN2()<=ncl) return 0;
1581 clrow = tpcrow->GetClusters2();
1585 AliTPCtrackerSector& tracksec = fOuterSec[(sec-fkNIS*2)%fkNOS];
1586 if (tracksec.GetNRows()<=row) return 0;
1587 tpcrow = &(tracksec[row]);
1588 if (tpcrow==0) return 0;
1590 if (sec-2*fkNIS<fkNOS) {
1591 if (tpcrow->GetN1()<=ncl) return 0;
1592 clrow = tpcrow->GetClusters1();
1595 if (tpcrow->GetN2()<=ncl) return 0;
1596 clrow = tpcrow->GetClusters2();
1600 return (AliTPCclusterMI*)clrow->At(ncl);
1606 Int_t AliTPCtrackerMI::FollowToNext(AliTPCseed& t, Int_t nr) {
1607 //-----------------------------------------------------------------
1608 // This function tries to find a track prolongation to next pad row
1609 //-----------------------------------------------------------------
1611 Double_t x= GetXrow(nr), ymax=GetMaxY(nr);
1614 AliTPCclusterMI *cl=0;
1615 Int_t tpcindex= t.GetClusterIndex2(nr);
1617 // update current shape info every 5 pad-row
1618 // if ( (nr%5==0) || t.GetNumberOfClusters()<2 || (t.fCurrentSigmaY2<0.0001) ){
1622 if (fIteration>0 && tpcindex>=-1){ //if we have already clusters
1624 if (tpcindex==-1) return 0; //track in dead zone
1625 if (tpcindex >= 0){ //
1626 cl = t.GetClusterPointer(nr);
1627 //if (cl==0) cl = GetClusterMI(tpcindex);
1628 if (!cl) cl = GetClusterMI(tpcindex);
1629 t.SetCurrentClusterIndex1(tpcindex);
1632 Int_t relativesector = ((tpcindex&0xff000000)>>24)%18; // if previously accepted cluster in different sector
1633 Float_t angle = relativesector*fSectors->GetAlpha()+fSectors->GetAlphaShift();
1635 if (angle<-TMath::Pi()) angle += 2*TMath::Pi();
1636 if (angle>=TMath::Pi()) angle -= 2*TMath::Pi();
1638 if (TMath::Abs(angle-t.GetAlpha())>0.001){
1639 Double_t rotation = angle-t.GetAlpha();
1640 t.SetRelativeSector(relativesector);
1641 if (!t.Rotate(rotation)) return 0;
1643 if (!t.PropagateTo(x)) return 0;
1645 t.SetCurrentCluster(cl);
1647 Int_t accept = AcceptCluster(&t,t.GetCurrentCluster());
1648 if ((tpcindex&0x8000)==0) accept =0;
1650 //if founded cluster is acceptible
1651 if (cl->IsUsed(11)) { // id cluster is shared inrease uncertainty
1652 t.SetErrorY2(t.GetErrorY2()+0.03);
1653 t.SetErrorZ2(t.GetErrorZ2()+0.03);
1654 t.SetErrorY2(t.GetErrorY2()*3);
1655 t.SetErrorZ2(t.GetErrorZ2()*3);
1657 t.SetNFoundable(t.GetNFoundable()+1);
1658 UpdateTrack(&t,accept);
1661 else { // Remove old cluster from track
1662 t.SetClusterIndex(nr, -3);
1663 t.SetClusterPointer(nr, 0);
1667 if (TMath::Abs(t.GetSnp())>AliTPCReconstructor::GetMaxSnpTracker()) return 0; // cut on angle
1668 if (fIteration>1 && IsFindable(t)){
1669 // not look for new cluster during refitting
1670 t.SetNFoundable(t.GetNFoundable()+1);
1675 // if (TMath::Abs(t.GetSnp())>0.95 || TMath::Abs(x*t.GetC()-t.GetEta())>0.95) return 0;// patch 28 fev 06
1676 if (!t.PropagateTo(x)) {
1677 if (fIteration==0) t.SetRemoval(10);
1680 Double_t y = t.GetY();
1681 if (TMath::Abs(y)>ymax){
1683 t.SetRelativeSector((t.GetRelativeSector()+1) % fN);
1684 if (!t.Rotate(fSectors->GetAlpha()))
1686 } else if (y <-ymax) {
1687 t.SetRelativeSector((t.GetRelativeSector()-1+fN) % fN);
1688 if (!t.Rotate(-fSectors->GetAlpha()))
1691 if (!t.PropagateTo(x)) {
1692 if (fIteration==0) t.SetRemoval(10);
1698 Double_t z=t.GetZ();
1701 if (!IsActive(t.GetRelativeSector(),nr)) {
1703 t.SetClusterIndex2(nr,-1);
1706 //AliInfo(Form("A - Sector%d phi %f - alpha %f", t.fRelativeSector,y/x, t.GetAlpha()));
1707 Bool_t isActive = IsActive(t.GetRelativeSector(),nr);
1708 Bool_t isActive2 = (nr>=fInnerSec->GetNRows()) ? fOuterSec[t.GetRelativeSector()][nr-fInnerSec->GetNRows()].GetN()>0:fInnerSec[t.GetRelativeSector()][nr].GetN()>0;
1710 if (!isActive || !isActive2) return 0;
1712 const AliTPCtrackerRow &krow=GetRow(t.GetRelativeSector(),nr);
1713 if ( (t.GetSigmaY2()<0) || t.GetSigmaZ2()<0) return 0;
1715 Double_t roadz = 1.;
1717 if (TMath::Abs(TMath::Abs(y)-ymax)<krow.GetDeadZone()){
1719 t.SetClusterIndex2(nr,-1);
1725 // if (TMath::Abs(z)<(AliTPCReconstructor::GetCtgRange()*x+10) && TMath::Abs(z)<fkParam->GetZLength(0) && (TMath::Abs(t.GetSnp())<AliTPCReconstructor::GetMaxSnpTracker()))
1726 t.SetNFoundable(t.GetNFoundable()+1);
1732 // cl = krow.FindNearest2(y+10.,z,roady,roadz,index);
1733 cl = krow.FindNearest2(y,z,roady,roadz,index);
1734 if (cl) t.SetCurrentClusterIndex1(krow.GetIndex(index));
1737 t.SetCurrentCluster(cl);
1739 if (fIteration==2&&cl->IsUsed(10)) return 0;
1740 Int_t accept = AcceptCluster(&t,t.GetCurrentCluster());
1741 if (fIteration==2&&cl->IsUsed(11)) {
1742 t.SetErrorY2(t.GetErrorY2()+0.03);
1743 t.SetErrorZ2(t.GetErrorZ2()+0.03);
1744 t.SetErrorY2(t.GetErrorY2()*3);
1745 t.SetErrorZ2(t.GetErrorZ2()*3);
1748 if (t.fCurrentCluster->IsUsed(10)){
1753 if (t.fNShared>0.7*t.GetNumberOfClusters()) {
1759 if (accept<3) UpdateTrack(&t,accept);
1762 if ( fIteration==0 && t.GetNFoundable()*0.5 > t.GetNumberOfClusters()) t.SetRemoval(10);
1770 //_________________________________________________________________________
1771 Bool_t AliTPCtrackerMI::GetTrackPoint(Int_t index, AliTrackPoint &p ) const
1773 // Get track space point by index
1774 // return false in case the cluster doesn't exist
1775 AliTPCclusterMI *cl = GetClusterMI(index);
1776 if (!cl) return kFALSE;
1777 Int_t sector = (index&0xff000000)>>24;
1778 // Int_t row = (index&0x00ff0000)>>16;
1780 // xyz[0] = fkParam->GetPadRowRadii(sector,row);
1781 xyz[0] = cl->GetX();
1782 xyz[1] = cl->GetY();
1783 xyz[2] = cl->GetZ();
1785 fkParam->AdjustCosSin(sector,cos,sin);
1786 Float_t x = cos*xyz[0]-sin*xyz[1];
1787 Float_t y = cos*xyz[1]+sin*xyz[0];
1789 Float_t sigmaY2 = 0.027*cl->GetSigmaY2();
1790 if (sector < fkParam->GetNInnerSector()) sigmaY2 *= 2.07;
1791 Float_t sigmaZ2 = 0.066*cl->GetSigmaZ2();
1792 if (sector < fkParam->GetNInnerSector()) sigmaZ2 *= 1.77;
1793 cov[0] = sin*sin*sigmaY2;
1794 cov[1] = -sin*cos*sigmaY2;
1796 cov[3] = cos*cos*sigmaY2;
1799 p.SetXYZ(x,y,xyz[2],cov);
1800 AliGeomManager::ELayerID iLayer;
1802 if (sector < fkParam->GetNInnerSector()) {
1803 iLayer = AliGeomManager::kTPC1;
1807 iLayer = AliGeomManager::kTPC2;
1808 idet = sector - fkParam->GetNInnerSector();
1810 UShort_t volid = AliGeomManager::LayerToVolUID(iLayer,idet);
1811 p.SetVolumeID(volid);
1817 Int_t AliTPCtrackerMI::UpdateClusters(AliTPCseed& t, Int_t nr) {
1818 //-----------------------------------------------------------------
1819 // This function tries to find a track prolongation to next pad row
1820 //-----------------------------------------------------------------
1821 t.SetCurrentCluster(0);
1822 t.SetCurrentClusterIndex1(-3);
1824 Double_t xt=t.GetX();
1825 Int_t row = GetRowNumber(xt)-1;
1826 Double_t ymax= GetMaxY(nr);
1828 if (row < nr) return 1; // don't prolongate if not information until now -
1829 // if (TMath::Abs(t.GetSnp())>0.9 && t.GetNumberOfClusters()>40. && fIteration!=2) {
1831 // return 0; // not prolongate strongly inclined tracks
1833 // if (TMath::Abs(t.GetSnp())>0.95) {
1835 // return 0; // not prolongate strongly inclined tracks
1836 // }// patch 28 fev 06
1838 Double_t x= GetXrow(nr);
1840 //t.PropagateTo(x+0.02);
1841 //t.PropagateTo(x+0.01);
1842 if (!t.PropagateTo(x)){
1849 if (TMath::Abs(y)>ymax){
1851 t.SetRelativeSector((t.GetRelativeSector()+1) % fN);
1852 if (!t.Rotate(fSectors->GetAlpha()))
1854 } else if (y <-ymax) {
1855 t.SetRelativeSector((t.GetRelativeSector()-1+fN) % fN);
1856 if (!t.Rotate(-fSectors->GetAlpha()))
1859 // if (!t.PropagateTo(x)){
1866 if (TMath::Abs(t.GetSnp())>AliTPCReconstructor::GetMaxSnpTracker()) return 0;
1868 if (!IsActive(t.GetRelativeSector(),nr)) {
1870 t.SetClusterIndex2(nr,-1);
1873 //AliInfo(Form("A - Sector%d phi %f - alpha %f", t.fRelativeSector,y/x, t.GetAlpha()));
1875 AliTPCtrackerRow &krow=GetRow(t.GetRelativeSector(),nr);
1877 if (TMath::Abs(TMath::Abs(y)-ymax)<krow.GetDeadZone()){
1879 t.SetClusterIndex2(nr,-1);
1885 // if (TMath::Abs(t.GetZ())<(AliTPCReconstructor::GetCtgRange()*t.GetX()+10) && (TMath::Abs(t.GetSnp())<AliTPCReconstructor::GetMaxSnpTracker()))
1886 if (IsFindable(t)) t.SetNFoundable(t.GetNFoundable()+1);
1892 if ( (nr%2==0) || t.GetNumberOfClusters()<2){
1893 // t.fCurrentSigmaY = GetSigmaY(&t);
1894 //t.fCurrentSigmaZ = GetSigmaZ(&t);
1898 AliTPCclusterMI *cl=0;
1901 Double_t roady = 1.;
1902 Double_t roadz = 1.;
1906 index = t.GetClusterIndex2(nr);
1907 if ( (index >= 0) && (index&0x8000)==0){
1908 cl = t.GetClusterPointer(nr);
1909 if ( (cl==0) && (index >= 0)) cl = GetClusterMI(index);
1910 t.SetCurrentClusterIndex1(index);
1912 t.SetCurrentCluster(cl);
1918 // if (index<0) return 0;
1919 UInt_t uindex = TMath::Abs(index);
1922 //cl = krow.FindNearest2(y+10,z,roady,roadz,uindex);
1923 cl = krow.FindNearest2(y,z,roady,roadz,uindex);
1926 if (cl) t.SetCurrentClusterIndex1(krow.GetIndex(uindex));
1927 t.SetCurrentCluster(cl);
1933 Int_t AliTPCtrackerMI::FollowToNextCluster(AliTPCseed & t, Int_t nr) {
1934 //-----------------------------------------------------------------
1935 // This function tries to find a track prolongation to next pad row
1936 //-----------------------------------------------------------------
1938 //update error according neighborhoud
1940 if (t.GetCurrentCluster()) {
1942 Int_t accept = AcceptCluster(&t,t.GetCurrentCluster());
1944 if (t.GetCurrentCluster()->IsUsed(10)){
1949 t.SetNShared(t.GetNShared()+1);
1950 if (t.GetNShared()>0.7*t.GetNumberOfClusters()) {
1955 if (fIteration>0) accept = 0;
1956 if (accept<3) UpdateTrack(&t,accept);
1960 if ( t.GetNumberOfClusters()>18 && ( (t.GetSigmaY2()+t.GetSigmaZ2())>0.16)) t.SetRemoval(10);
1961 if ( t.GetNumberOfClusters()>18 && t.GetChi2()/t.GetNumberOfClusters()>6 ) t.SetRemoval(10);
1963 if (( (t.GetNFoundable()*0.5 > t.GetNumberOfClusters()) || t.GetNoCluster()>15)) t.SetRemoval(10);
1971 //_____________________________________________________________________________
1972 Int_t AliTPCtrackerMI::FollowProlongation(AliTPCseed& t, Int_t rf, Int_t step, Bool_t fromSeeds) {
1973 //-----------------------------------------------------------------
1974 // This function tries to find a track prolongation.
1975 //-----------------------------------------------------------------
1976 Double_t xt=t.GetX();
1978 Double_t alpha=t.GetAlpha();
1979 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
1980 if (alpha < 0. ) alpha += 2.*TMath::Pi();
1982 t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN);
1984 Int_t first = GetRowNumber(xt);
1989 for (Int_t nr= first; nr>=rf; nr-=step) {
1991 if (t.GetKinkIndexes()[0]>0){
1992 for (Int_t i=0;i<3;i++){
1993 Int_t index = t.GetKinkIndexes()[i];
1994 if (index==0) break;
1995 if (index<0) continue;
1997 AliKink * kink = (AliKink*)fEvent->GetKink(index-1);
1999 printf("PROBLEM\n");
2002 Int_t kinkrow = kink->GetTPCRow0()+2+Int_t(0.5/(0.05+kink->GetAngle(2)));
2004 AliExternalTrackParam paramd(t);
2005 kink->SetDaughter(paramd);
2006 kink->SetStatus(2,5);
2013 if (nr==80) t.UpdateReference();
2014 if (nr<fInnerSec->GetNRows())
2015 fSectors = fInnerSec;
2017 fSectors = fOuterSec;
2018 if (FollowToNext(t,nr)==0)
2031 Int_t AliTPCtrackerMI::FollowBackProlongation(AliTPCseed& t, Int_t rf, Bool_t fromSeeds) {
2032 //-----------------------------------------------------------------
2033 // This function tries to find a track prolongation.
2034 //-----------------------------------------------------------------
2036 Double_t xt=t.GetX();
2037 Double_t alpha=t.GetAlpha();
2038 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
2039 if (alpha < 0. ) alpha += 2.*TMath::Pi();
2040 t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN);
2042 Int_t first = t.GetFirstPoint();
2043 Int_t ri = GetRowNumber(xt);
2047 if (first<ri) first = ri;
2049 if (first<0) first=0;
2050 for (Int_t nr=first; nr<=rf; nr++) {
2051 // if ( (TMath::Abs(t.GetSnp())>0.95)) break;//patch 28 fev 06
2052 if (t.GetKinkIndexes()[0]<0){
2053 for (Int_t i=0;i<3;i++){
2054 Int_t index = t.GetKinkIndexes()[i];
2055 if (index==0) break;
2056 if (index>0) continue;
2057 index = TMath::Abs(index);
2058 AliKink * kink = (AliKink*)fEvent->GetKink(index-1);
2060 printf("PROBLEM\n");
2063 Int_t kinkrow = kink->GetTPCRow0()-2-Int_t(0.5/(0.05+kink->GetAngle(2)));
2065 AliExternalTrackParam paramm(t);
2066 kink->SetMother(paramm);
2067 kink->SetStatus(2,1);
2074 if (nr<fInnerSec->GetNRows())
2075 fSectors = fInnerSec;
2077 fSectors = fOuterSec;
2088 Float_t AliTPCtrackerMI::OverlapFactor(AliTPCseed * s1, AliTPCseed * s2, Int_t &sum1, Int_t & sum2)
2090 // overlapping factor
2096 Float_t dz2 =(s1->GetZ() - s2->GetZ());
2099 Float_t dy2 =TMath::Abs((s1->GetY() - s2->GetY()));
2101 Float_t distance = TMath::Sqrt(dz2+dy2);
2102 if (distance>4.) return 0; // if there are far away - not overlap - to reduce combinatorics
2105 Int_t firstpoint = TMath::Min(s1->GetFirstPoint(),s2->GetFirstPoint());
2106 Int_t lastpoint = TMath::Max(s1->GetLastPoint(),s2->GetLastPoint());
2111 if (firstpoint>lastpoint) {
2112 firstpoint =lastpoint;
2117 for (Int_t i=firstpoint-1;i<lastpoint+1;i++){
2118 if (s1->GetClusterIndex2(i)>0) sum1++;
2119 if (s2->GetClusterIndex2(i)>0) sum2++;
2120 if (s1->GetClusterIndex2(i)==s2->GetClusterIndex2(i) && s1->GetClusterIndex2(i)>0) {
2124 if (sum<5) return 0;
2126 Float_t summin = TMath::Min(sum1+1,sum2+1);
2127 Float_t ratio = (sum+1)/Float_t(summin);
2131 void AliTPCtrackerMI::SignShared(AliTPCseed * s1, AliTPCseed * s2)
2135 Float_t thetaCut = 0.2;//+10.*TMath::Sqrt(s1->GetSigmaTglZ()+ s2->GetSigmaTglZ());
2136 if (TMath::Abs(s1->GetTgl()-s2->GetTgl())>thetaCut) return;
2137 Float_t minCl = TMath::Min(s1->GetNumberOfClusters(),s2->GetNumberOfClusters());
2138 Int_t cutN0 = TMath::Max(5,TMath::Nint(0.1*minCl));
2143 //Int_t firstpoint = TMath::Max(s1->GetFirstPoint(),s2->GetFirstPoint());
2144 //Int_t lastpoint = TMath::Min(s1->GetLastPoint(),s2->GetLastPoint());
2145 Int_t firstpoint = 0;
2146 Int_t lastpoint = 160;
2148 // if (firstpoint>=lastpoint-5) return;;
2150 for (Int_t i=firstpoint;i<lastpoint;i++){
2151 // if ( (s1->GetClusterIndex2(i)&0xFFFF8FFF)==(s2->GetClusterIndex2(i)&0xFFFF8FFF) && s1->GetClusterIndex2(i)>0) {
2152 if ( (s1->GetClusterIndex2(i))==(s2->GetClusterIndex2(i)) && s1->GetClusterIndex2(i)>0) {
2156 if (sumshared>cutN0){
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) {
2162 AliTPCTrackerPoint *p1 = s1->GetTrackPoint(i);
2163 AliTPCTrackerPoint *p2 = s2->GetTrackPoint(i);;
2164 if (s1->IsActive()&&s2->IsActive()){
2165 p1->SetShared(kTRUE);
2166 p2->SetShared(kTRUE);
2172 if (sumshared>cutN0){
2173 for (Int_t i=0;i<4;i++){
2174 if (s1->GetOverlapLabel(3*i)==0){
2175 s1->SetOverlapLabel(3*i, s2->GetLabel());
2176 s1->SetOverlapLabel(3*i+1,sumshared);
2177 s1->SetOverlapLabel(3*i+2,s2->GetUniqueID());
2181 for (Int_t i=0;i<4;i++){
2182 if (s2->GetOverlapLabel(3*i)==0){
2183 s2->SetOverlapLabel(3*i, s1->GetLabel());
2184 s2->SetOverlapLabel(3*i+1,sumshared);
2185 s2->SetOverlapLabel(3*i+2,s1->GetUniqueID());
2192 void AliTPCtrackerMI::SignShared(TObjArray * arr)
2195 //sort trackss according sectors
2197 for (Int_t i=0; i<arr->GetEntriesFast(); i++) {
2198 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2200 //if (pt) RotateToLocal(pt);
2204 arr->Sort(); // sorting according relative sectors
2205 arr->Expand(arr->GetEntries());
2208 Int_t nseed=arr->GetEntriesFast();
2209 for (Int_t i=0; i<nseed; i++) {
2210 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2212 for (Int_t j=0;j<12;j++){
2213 pt->SetOverlapLabel(j,0);
2216 for (Int_t i=0; i<nseed; i++) {
2217 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2219 if (pt->GetRemoval()>10) continue;
2220 for (Int_t j=i+1; j<nseed; j++){
2221 AliTPCseed *pt2=(AliTPCseed*)arr->UncheckedAt(j);
2222 if (TMath::Abs(pt->GetRelativeSector()-pt2->GetRelativeSector())>1) continue;
2224 if (pt2->GetRemoval()<=10) {
2225 //if ( TMath::Abs(pt->GetRelativeSector()-pt2->GetRelativeSector())>0) break;
2233 void AliTPCtrackerMI::SortTracks(TObjArray * arr, Int_t mode) const
2236 //sort tracks in array according mode criteria
2237 Int_t nseed = arr->GetEntriesFast();
2238 for (Int_t i=0; i<nseed; i++) {
2239 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2250 void AliTPCtrackerMI::RemoveUsed2(TObjArray * arr, Float_t factor1, Float_t factor2, Int_t minimal)
2253 // Loop over all tracks and remove overlaped tracks (with lower quality)
2255 // 1. Unsign clusters
2256 // 2. Sort tracks according quality
2257 // Quality is defined by the number of cluster between first and last points
2259 // 3. Loop over tracks - decreasing quality order
2260 // a.) remove - If the fraction of shared cluster less than factor (1- n or 2)
2261 // b.) remove - If the minimal number of clusters less than minimal and not ITS
2262 // c.) if track accepted - sign clusters
2264 //Called in - AliTPCtrackerMI::Clusters2Tracks()
2265 // - AliTPCtrackerMI::PropagateBack()
2266 // - AliTPCtrackerMI::RefitInward()
2269 // factor1 - factor for constrained
2270 // factor2 - for non constrained tracks
2271 // if (Float_t(shared+1)/Float_t(found+1)>factor) - DELETE
2275 Int_t nseed = arr->GetEntriesFast();
2276 Float_t * quality = new Float_t[nseed];
2277 Int_t * indexes = new Int_t[nseed];
2281 for (Int_t i=0; i<nseed; i++) {
2282 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2287 pt->UpdatePoints(); //select first last max dens points
2288 Float_t * points = pt->GetPoints();
2289 if (points[3]<0.8) quality[i] =-1;
2290 quality[i] = (points[2]-points[0])+pt->GetNumberOfClusters();
2291 //prefer high momenta tracks if overlaps
2292 quality[i] *= TMath::Sqrt(TMath::Abs(pt->Pt())+0.5);
2294 TMath::Sort(nseed,quality,indexes);
2297 for (Int_t itrack=0; itrack<nseed; itrack++) {
2298 Int_t trackindex = indexes[itrack];
2299 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(trackindex);
2302 if (quality[trackindex]<0){
2303 fSeedsPool->MarkSlotFree( arr->RemoveAt(trackindex) );
2308 Int_t first = Int_t(pt->GetPoints()[0]);
2309 Int_t last = Int_t(pt->GetPoints()[2]);
2310 Double_t factor = (pt->GetBConstrain()) ? factor1: factor2;
2312 Int_t found,foundable,shared;
2313 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
2314 // pt->GetClusterStatistic(0,160, found, foundable,shared,kFALSE);
2315 Bool_t itsgold =kFALSE;
2318 if (pt->GetESD()->GetITSclusters(dummy)>4) itsgold= kTRUE;
2322 if (Float_t(shared+1)/Float_t(found+1)>factor){
2323 if (pt->GetKinkIndexes()[0]!=0) continue; //don't remove tracks - part of the kinks
2324 if( AliTPCReconstructor::StreamLevel()>3){
2325 TTreeSRedirector &cstream = *fDebugStreamer;
2326 cstream<<"RemoveUsed"<<
2327 "iter="<<fIteration<<
2331 fSeedsPool->MarkSlotFree( arr->RemoveAt(trackindex) );
2334 if (pt->GetNumberOfClusters()<50&&(found-0.5*shared)<minimal){ //remove short tracks
2335 if (pt->GetKinkIndexes()[0]!=0) continue; //don't remove tracks - part of the kinks
2336 if( AliTPCReconstructor::StreamLevel()>3){
2337 TTreeSRedirector &cstream = *fDebugStreamer;
2338 cstream<<"RemoveShort"<<
2339 "iter="<<fIteration<<
2343 fSeedsPool->MarkSlotFree( arr->RemoveAt(trackindex) );
2349 //if (sharedfactor>0.4) continue;
2350 if (pt->GetKinkIndexes()[0]>0) continue;
2351 //Remove tracks with undefined properties - seems
2352 if (pt->GetSigmaY2()<kAlmost0) continue; // ? what is the origin ?
2354 for (Int_t i=first; i<last; i++) {
2355 Int_t index=pt->GetClusterIndex2(i);
2356 // if (index<0 || index&0x8000 ) continue;
2357 if (index<0 || index&0x8000 ) continue;
2358 AliTPCclusterMI *c= pt->GetClusterPointer(i);
2365 Info("RemoveUsed","\n*****\nNumber of good tracks after shared removal\t%d\n",fNtracks);
2371 void AliTPCtrackerMI::DumpClusters(Int_t iter, TObjArray *trackArray)
2374 // Dump clusters after reco
2375 // signed and unsigned cluster can be visualized
2376 // 1. Unsign all cluster
2377 // 2. Sign all used clusters
2380 Int_t nseed = trackArray->GetEntries();
2381 for (Int_t i=0; i<nseed; i++){
2382 AliTPCseed *pt=(AliTPCseed*)trackArray->UncheckedAt(i);
2386 Bool_t isKink=pt->GetKinkIndex(0)!=0;
2387 for (Int_t j=0; j<160; ++j) {
2388 Int_t index=pt->GetClusterIndex2(j);
2389 if (index<0) continue;
2390 AliTPCclusterMI *c= pt->GetClusterPointer(j);
2392 if (isKink) c->Use(100); // kink
2393 c->Use(10); // by default usage 10
2398 for (Int_t sec=0;sec<fkNIS;sec++){
2399 for (Int_t row=0;row<fInnerSec->GetNRows();row++){
2400 TClonesArray *cla = fInnerSec[sec][row].GetClusters1();
2401 for (Int_t icl =0;icl< fInnerSec[sec][row].GetN1();icl++){
2402 AliTPCclusterMI* cl = (AliTPCclusterMI*)cla->At(icl);
2403 Float_t gx[3]; cl->GetGlobalXYZ(gx);
2404 (*fDebugStreamer)<<"clDump"<<
2412 cla = fInnerSec[sec][row].GetClusters2();
2413 for (Int_t icl =0;icl< fInnerSec[sec][row].GetN2();icl++){
2414 AliTPCclusterMI* cl = (AliTPCclusterMI*)cla->At(icl);
2415 Float_t gx[3]; cl->GetGlobalXYZ(gx);
2416 (*fDebugStreamer)<<"clDump"<<
2427 for (Int_t sec=0;sec<fkNOS;sec++){
2428 for (Int_t row=0;row<fOuterSec->GetNRows();row++){
2429 TClonesArray *cla = fOuterSec[sec][row].GetClusters1();
2430 for (Int_t icl =0;icl< fOuterSec[sec][row].GetN1();icl++){
2432 AliTPCclusterMI* cl = (AliTPCclusterMI*) cla->At(icl);
2433 cl->GetGlobalXYZ(gx);
2434 (*fDebugStreamer)<<"clDump"<<
2442 cla = fOuterSec[sec][row].GetClusters2();
2443 for (Int_t icl =0;icl< fOuterSec[sec][row].GetN2();icl++){
2445 AliTPCclusterMI* cl = (AliTPCclusterMI*) cla->At(icl);
2446 cl->GetGlobalXYZ(gx);
2447 (*fDebugStreamer)<<"clDump"<<
2459 void AliTPCtrackerMI::UnsignClusters()
2462 // loop over all clusters and unsign them
2465 for (Int_t sec=0;sec<fkNIS;sec++){
2466 for (Int_t row=0;row<fInnerSec->GetNRows();row++){
2467 TClonesArray *cla = fInnerSec[sec][row].GetClusters1();
2468 for (Int_t icl =0;icl< fInnerSec[sec][row].GetN1();icl++)
2469 // if (cl[icl].IsUsed(10))
2470 ((AliTPCclusterMI*) cla->At(icl))->Use(-1);
2471 cla = fInnerSec[sec][row].GetClusters2();
2472 for (Int_t icl =0;icl< fInnerSec[sec][row].GetN2();icl++)
2473 //if (cl[icl].IsUsed(10))
2474 ((AliTPCclusterMI*) cla->At(icl))->Use(-1);
2478 for (Int_t sec=0;sec<fkNOS;sec++){
2479 for (Int_t row=0;row<fOuterSec->GetNRows();row++){
2480 TClonesArray *cla = fOuterSec[sec][row].GetClusters1();
2481 for (Int_t icl =0;icl< fOuterSec[sec][row].GetN1();icl++)
2482 //if (cl[icl].IsUsed(10))
2483 ((AliTPCclusterMI*) cla->At(icl))->Use(-1);
2484 cla = fOuterSec[sec][row].GetClusters2();
2485 for (Int_t icl =0;icl< fOuterSec[sec][row].GetN2();icl++)
2486 //if (cl[icl].IsUsed(10))
2487 ((AliTPCclusterMI*) cla->At(icl))->Use(-1);
2495 void AliTPCtrackerMI::SignClusters(const TObjArray * arr, Float_t fnumber, Float_t fdensity)
2498 //sign clusters to be "used"
2500 // snumber and sdensity sign number of sigmas - bellow mean value to be accepted
2501 // loop over "primaries"
2515 Int_t nseed = arr->GetEntriesFast();
2516 for (Int_t i=0; i<nseed; i++) {
2517 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2521 if (!(pt->IsActive())) continue;
2522 Float_t dens = pt->GetNumberOfClusters()/Float_t(pt->GetNFoundable());
2523 if ( (dens>0.7) && (pt->GetNumberOfClusters()>70)){
2525 sumdens2+= dens*dens;
2526 sumn += pt->GetNumberOfClusters();
2527 sumn2 += pt->GetNumberOfClusters()*pt->GetNumberOfClusters();
2528 Float_t chi2 = pt->GetChi2()/pt->GetNumberOfClusters();
2531 sumchi2 +=chi2*chi2;
2536 Float_t mdensity = 0.9;
2537 Float_t meann = 130;
2538 Float_t meanchi = 1;
2539 Float_t sdensity = 0.1;
2540 Float_t smeann = 10;
2541 Float_t smeanchi =0.4;
2545 mdensity = sumdens/sum;
2547 meanchi = sumchi/sum;
2549 sdensity = sumdens2/sum-mdensity*mdensity;
2551 sdensity = TMath::Sqrt(sdensity);
2555 smeann = sumn2/sum-meann*meann;
2557 smeann = TMath::Sqrt(smeann);
2561 smeanchi = sumchi2/sum - meanchi*meanchi;
2563 smeanchi = TMath::Sqrt(smeanchi);
2569 //REMOVE SHORT DELTAS or tracks going out of sensitive volume of TPC
2571 for (Int_t i=0; i<nseed; i++) {
2572 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2576 if (pt->GetBSigned()) continue;
2577 if (pt->GetBConstrain()) continue;
2578 //if (!(pt->IsActive())) continue;
2580 Int_t found,foundable,shared;
2581 pt->GetClusterStatistic(0,160,found, foundable,shared);
2582 if (shared/float(found)>0.3) {
2583 if (shared/float(found)>0.9 ){
2584 //fSeedsPool->MarkSlotFree( arr->RemoveAt(i) );
2589 Bool_t isok =kFALSE;
2590 if ( (pt->GetNShared()/pt->GetNumberOfClusters()<0.5) &&pt->GetNumberOfClusters()>60)
2592 if ((TMath::Abs(1/pt->GetC())<100.) && (pt->GetNShared()/pt->GetNumberOfClusters()<0.7))
2594 if (TMath::Abs(pt->GetZ()/pt->GetX())>1.1)
2596 if ( (TMath::Abs(pt->GetSnp()>0.7) && pt->GetD(0,0)>60.))
2600 for (Int_t j=0; j<160; ++j) {
2601 Int_t index=pt->GetClusterIndex2(j);
2602 if (index<0) continue;
2603 AliTPCclusterMI *c= pt->GetClusterPointer(j);
2605 //if (!(c->IsUsed(10))) c->Use();
2612 Double_t maxchi = meanchi+2.*smeanchi;
2614 for (Int_t i=0; i<nseed; i++) {
2615 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2619 //if (!(pt->IsActive())) continue;
2620 if (pt->GetBSigned()) continue;
2621 Double_t chi = pt->GetChi2()/pt->GetNumberOfClusters();
2622 if (chi>maxchi) continue;
2625 Float_t dens = pt->GetNumberOfClusters()/Float_t(pt->GetNFoundable());
2627 //sign only tracks with enoug big density at the beginning
2629 if ((pt->GetDensityFirst(40)<0.75) && pt->GetNumberOfClusters()<meann) continue;
2632 Double_t mindens = TMath::Max(double(mdensity-sdensity*fdensity*bfactor),0.65);
2633 Double_t minn = TMath::Max(Int_t(meann-fnumber*smeann*bfactor),50);
2635 // if (pt->fBConstrain) mindens = TMath::Max(mdensity-sdensity*fdensity*bfactor,0.65);
2636 if ( (pt->GetRemoval()==10) && (pt->GetSnp()>0.8)&&(dens>mindens))
2639 if ((dens>mindens && pt->GetNumberOfClusters()>minn) && chi<maxchi ){
2640 //Int_t noc=pt->GetNumberOfClusters();
2641 pt->SetBSigned(kTRUE);
2642 for (Int_t j=0; j<160; ++j) {
2644 Int_t index=pt->GetClusterIndex2(j);
2645 if (index<0) continue;
2646 AliTPCclusterMI *c= pt->GetClusterPointer(j);
2648 // if (!(c->IsUsed(10))) c->Use();
2653 // gLastCheck = nseed;
2662 Int_t AliTPCtrackerMI::RefitInward(AliESDEvent *event)
2665 // back propagation of ESD tracks
2668 if (!event) return 0;
2669 const Int_t kMaxFriendTracks=2000;
2671 // extract correction object for multiplicity dependence of dEdx
2672 TObjArray * gainCalibArray = AliTPCcalibDB::Instance()->GetTimeGainSplinesRun(event->GetRunNumber());
2674 AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ;
2676 AliFatal("Tranformations not in RefitInward");
2679 transform->SetCurrentRecoParam((AliTPCRecoParam*)AliTPCReconstructor::GetRecoParam());
2680 const AliTPCRecoParam * recoParam = AliTPCcalibDB::Instance()->GetTransform()->GetCurrentRecoParam();
2681 Int_t nContribut = event->GetNumberOfTracks();
2682 TGraphErrors * graphMultDependenceDeDx = 0x0;
2683 if (recoParam && recoParam->GetUseMultiplicityCorrectionDedx() && gainCalibArray) {
2684 if (recoParam->GetUseTotCharge()) {
2685 graphMultDependenceDeDx = (TGraphErrors *) gainCalibArray->FindObject("TGRAPHERRORS_MEANQTOT_MULTIPLICITYDEPENDENCE_BEAM_ALL");
2687 graphMultDependenceDeDx = (TGraphErrors *) gainCalibArray->FindObject("TGRAPHERRORS_MEANQMAX_MULTIPLICITYDEPENDENCE_BEAM_ALL");
2693 //PrepareForProlongation(fSeeds,1);
2694 PropagateForward2(fSeeds);
2695 RemoveUsed2(fSeeds,0.4,0.4,20);
2697 TObjArray arraySeed(fSeeds->GetEntries());
2698 for (Int_t i=0;i<fSeeds->GetEntries();i++) {
2699 arraySeed.AddAt(fSeeds->At(i),i);
2701 SignShared(&arraySeed);
2702 // FindCurling(fSeeds, event,2); // find multi found tracks
2703 FindSplitted(fSeeds, event,2); // find multi found tracks
2704 if (AliTPCReconstructor::StreamLevel()>5) FindMultiMC(fSeeds, fEvent,2); // find multi found tracks
2707 Int_t nseed = fSeeds->GetEntriesFast();
2708 for (Int_t i=0;i<nseed;i++){
2709 AliTPCseed * seed = (AliTPCseed*) fSeeds->UncheckedAt(i);
2710 if (!seed) continue;
2711 if (seed->GetKinkIndex(0)>0) UpdateKinkQualityD(seed); // update quality informations for kinks
2712 AliESDtrack *esd=event->GetTrack(i);
2714 if (seed->GetNumberOfClusters()<60 && seed->GetNumberOfClusters()<(esd->GetTPCclusters(0) -5)*0.8){
2715 AliExternalTrackParam paramIn;
2716 AliExternalTrackParam paramOut;
2717 Int_t ncl = seed->RefitTrack(seed,¶mIn,¶mOut);
2718 if (AliTPCReconstructor::StreamLevel()>2) {
2719 (*fDebugStreamer)<<"RecoverIn"<<
2723 "pout.="<<¶mOut<<
2728 seed->Set(paramIn.GetX(),paramIn.GetAlpha(),paramIn.GetParameter(),paramIn.GetCovariance());
2729 seed->SetNumberOfClusters(ncl);
2733 seed->PropagateTo(fkParam->GetInnerRadiusLow());
2734 seed->UpdatePoints();
2735 AddCovariance(seed);
2736 MakeESDBitmaps(seed, esd);
2737 seed->CookdEdx(0.02,0.6);
2738 CookLabel(seed,0.1); //For comparison only
2740 if (AliTPCReconstructor::StreamLevel()>1 && seed!=0) {
2741 TTreeSRedirector &cstream = *fDebugStreamer;
2748 if (seed->GetNumberOfClusters()>15){
2749 esd->UpdateTrackParams(seed,AliESDtrack::kTPCrefit);
2750 esd->SetTPCPoints(seed->GetPoints());
2751 esd->SetTPCPointsF(seed->GetNFoundable());
2752 Int_t ndedx = seed->GetNCDEDX(0);
2753 Float_t sdedx = seed->GetSDEDX(0);
2754 Float_t dedx = seed->GetdEdx();
2755 // apply mutliplicity dependent dEdx correction if available
2756 if (graphMultDependenceDeDx) {
2757 Double_t corrGain = AliTPCcalibDButil::EvalGraphConst(graphMultDependenceDeDx, nContribut);
2758 dedx += (1 - corrGain)*50.; // MIP is normalized to 50
2760 esd->SetTPCsignal(dedx, sdedx, ndedx);
2762 // fill new dEdx information
2764 Double32_t signal[4];
2768 for(Int_t iarr=0;iarr<3;iarr++) {
2769 signal[iarr] = seed->GetDEDXregion(iarr+1);
2770 ncl[iarr] = seed->GetNCDEDX(iarr+1);
2771 nrows[iarr] = seed->GetNCDEDXInclThres(iarr+1);
2773 signal[3] = seed->GetDEDXregion(4);
2775 AliClonesPool* pooldEdx = fPools ? fPools->GetPoolTPCdEdx() : 0;
2776 AliTPCdEdxInfo * infoTpcPid = 0;
2778 infoTpcPid = new(pooldEdx->NextFreeSlot()) AliTPCdEdxInfo();
2779 pooldEdx->RegisterClone(pooldEdx);
2781 else infoTpcPid = new AliTPCdEdxInfo();
2782 infoTpcPid->SetTPCSignalRegionInfo(signal, ncl, nrows);
2783 esd->SetTPCdEdxInfo(infoTpcPid);
2785 // add seed to the esd track in Calib level
2787 Bool_t storeFriend = gRandom->Rndm()<(kMaxFriendTracks)/Float_t(nseed);
2788 if (AliTPCReconstructor::StreamLevel()>0 &&storeFriend){
2789 AliTPCseed * seedCopy = new( fSeedsPool->NextFreeSlot() ) AliTPCseed(*seed, kTRUE);
2790 fSeedsPool->RegisterClone(seedCopy);
2791 esd->AddCalibObject(seedCopy);
2796 //printf("problem\n");
2799 //FindKinks(fSeeds,event);
2800 if (AliTPCReconstructor::StreamLevel()>3) DumpClusters(2,fSeeds);
2801 Info("RefitInward","Number of refitted tracks %d",ntracks);
2806 Int_t AliTPCtrackerMI::PropagateBack(AliESDEvent *event)
2809 // back propagation of ESD tracks
2811 if (!event) return 0;
2815 PropagateBack(fSeeds);
2816 RemoveUsed2(fSeeds,0.4,0.4,20);
2817 //FindCurling(fSeeds, fEvent,1);
2818 FindSplitted(fSeeds, event,1); // find multi found tracks
2819 if (AliTPCReconstructor::StreamLevel()>5) FindMultiMC(fSeeds, fEvent,1); // find multi found tracks
2822 Int_t nseed = fSeeds->GetEntriesFast();
2824 for (Int_t i=0;i<nseed;i++){
2825 AliTPCseed * seed = (AliTPCseed*) fSeeds->UncheckedAt(i);
2826 if (!seed) continue;
2827 if (seed->GetKinkIndex(0)<0) UpdateKinkQualityM(seed); // update quality informations for kinks
2828 seed->UpdatePoints();
2829 AddCovariance(seed);
2830 AliESDtrack *esd=event->GetTrack(i);
2831 if (!esd) continue; //never happen
2832 if (seed->GetNumberOfClusters()<60 && seed->GetNumberOfClusters()<(esd->GetTPCclusters(0) -5)*0.8){
2833 AliExternalTrackParam paramIn;
2834 AliExternalTrackParam paramOut;
2835 Int_t ncl = seed->RefitTrack(seed,¶mIn,¶mOut);
2836 if (AliTPCReconstructor::StreamLevel()>2) {
2837 (*fDebugStreamer)<<"RecoverBack"<<
2841 "pout.="<<¶mOut<<
2846 seed->Set(paramOut.GetX(),paramOut.GetAlpha(),paramOut.GetParameter(),paramOut.GetCovariance());
2847 seed->SetNumberOfClusters(ncl);
2850 seed->CookdEdx(0.02,0.6);
2851 CookLabel(seed,0.1); //For comparison only
2852 if (seed->GetNumberOfClusters()>15){
2853 esd->UpdateTrackParams(seed,AliESDtrack::kTPCout);
2854 esd->SetTPCPoints(seed->GetPoints());
2855 esd->SetTPCPointsF(seed->GetNFoundable());
2856 Int_t ndedx = seed->GetNCDEDX(0);
2857 Float_t sdedx = seed->GetSDEDX(0);
2858 Float_t dedx = seed->GetdEdx();
2859 esd->SetTPCsignal(dedx, sdedx, ndedx);
2861 Int_t eventnumber = event->GetEventNumberInFile();// patch 28 fev 06
2862 // This is most likely NOT the event number you'd like to use. It has nothing to do with the 'real' event number
2863 if (AliTPCReconstructor::StreamLevel()>1 && esd) {
2864 (*fDebugStreamer)<<"Cback"<<
2867 "EventNrInFile="<<eventnumber<<
2872 if (AliTPCReconstructor::StreamLevel()>3) DumpClusters(1,fSeeds);
2873 //FindKinks(fSeeds,event);
2874 Info("PropagateBack","Number of back propagated tracks %d",ntracks);
2881 void AliTPCtrackerMI::DeleteSeeds()
2885 fSeedsPool->Clear("C"); // nominally allocate memory
2890 void AliTPCtrackerMI::ReadSeeds(const AliESDEvent *const event, Int_t direction)
2893 //read seeds from the event
2895 Int_t nentr=event->GetNumberOfTracks();
2897 Info("ReadSeeds", "Number of ESD tracks: %d\n", nentr);
2902 fSeeds = new TObjArray(nentr);
2906 for (Int_t i=0; i<nentr; i++) {
2907 AliESDtrack *esd=event->GetTrack(i);
2908 ULong_t status=esd->GetStatus();
2909 if (!(status&AliESDtrack::kTPCin)) continue;
2910 AliTPCtrack t(*esd);
2911 t.SetNumberOfClusters(0);
2912 // AliTPCseed *seed = new AliTPCseed(t,t.GetAlpha());
2913 AliTPCseed *seed = new( fSeedsPool->NextFreeSlot() ) AliTPCseed(t/*,t.GetAlpha()*/);
2914 fSeedsPool->RegisterClone(seed);
2915 seed->SetUniqueID(esd->GetID());
2916 AddCovariance(seed); //add systematic ucertainty
2917 for (Int_t ikink=0;ikink<3;ikink++) {
2918 Int_t index = esd->GetKinkIndex(ikink);
2919 seed->GetKinkIndexes()[ikink] = index;
2920 if (index==0) continue;
2921 index = TMath::Abs(index);
2922 AliESDkink * kink = fEvent->GetKink(index-1);
2923 if (kink&&esd->GetKinkIndex(ikink)<0){
2924 if ((status & AliESDtrack::kTRDrefit) != 0) kink->SetStatus(1,2);
2925 if ((status & AliESDtrack::kITSout) != 0) kink->SetStatus(1,0);
2927 if (kink&&esd->GetKinkIndex(ikink)>0){
2928 if ((status & AliESDtrack::kTRDrefit) != 0) kink->SetStatus(1,6);
2929 if ((status & AliESDtrack::kITSout) != 0) kink->SetStatus(1,4);
2933 if (((status&AliESDtrack::kITSout)==0)&&(direction==1)) seed->ResetCovariance(10.);
2934 if ( direction ==2 &&(status & AliESDtrack::kTRDrefit) == 0 ) seed->ResetCovariance(10.);
2935 //if ( direction ==2 && ((status & AliESDtrack::kTPCout) == 0) ) {
2936 // fSeeds->AddAt(0,i);
2937 // fSeedsPool->MarkSlotFree( seed );
2940 if ( direction ==2 &&(status & AliESDtrack::kTRDrefit) > 0 ) {
2941 Double_t par0[5],par1[5],alpha,x;
2942 esd->GetInnerExternalParameters(alpha,x,par0);
2943 esd->GetExternalParameters(x,par1);
2944 Double_t delta1 = TMath::Abs(par0[4]-par1[4])/(0.000000001+TMath::Abs(par0[4]+par1[4]));
2945 Double_t delta2 = TMath::Abs(par0[3]-par1[3]);
2947 if (esd->GetTRDncls()>0) trdchi2 = esd->GetTRDchi2()/esd->GetTRDncls();
2948 //reset covariance if suspicious
2949 if ( (delta1>0.1) || (delta2>0.006) ||trdchi2>7.)
2950 seed->ResetCovariance(10.);
2955 // rotate to the local coordinate system
2957 fSectors=fInnerSec; fN=fkNIS;
2958 Double_t alpha=seed->GetAlpha();
2959 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
2960 if (alpha < 0. ) alpha += 2.*TMath::Pi();
2961 Int_t ns=Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN;
2962 alpha =ns*fSectors->GetAlpha() + fSectors->GetAlphaShift();
2963 alpha-=seed->GetAlpha();
2964 if (alpha<-TMath::Pi()) alpha += 2*TMath::Pi();
2965 if (alpha>=TMath::Pi()) alpha -= 2*TMath::Pi();
2966 if (TMath::Abs(alpha) > 0.001) { // This should not happen normally
2967 AliWarning(Form("Rotating track over %f",alpha));
2968 if (!seed->Rotate(alpha)) {
2969 fSeedsPool->MarkSlotFree( seed );
2975 if (esd->GetKinkIndex(0)<=0){
2976 for (Int_t irow=0;irow<160;irow++){
2977 Int_t index = seed->GetClusterIndex2(irow);
2980 AliTPCclusterMI * cl = GetClusterMI(index);
2981 seed->SetClusterPointer(irow,cl);
2983 if ((index & 0x8000)==0){
2984 cl->Use(10); // accepted cluster
2986 cl->Use(6); // close cluster not accepted
2989 Info("ReadSeeds","Not found cluster");
2994 fSeeds->AddAt(seed,i);
3000 //_____________________________________________________________________________
3001 void AliTPCtrackerMI::MakeSeeds3(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, Float_t cuts[4],
3002 Float_t deltay, Int_t ddsec) {
3003 //-----------------------------------------------------------------
3004 // This function creates track seeds.
3005 // SEEDING WITH VERTEX CONSTRAIN
3006 //-----------------------------------------------------------------
3007 // cuts[0] - fP4 cut
3008 // cuts[1] - tan(phi) cut
3009 // cuts[2] - zvertex cut
3010 // cuts[3] - fP3 cut
3018 Double_t x[5], c[15];
3019 // Int_t di = i1-i2;
3021 AliTPCseed * seed = new( fSeedsPool->NextFreeSlot() ) AliTPCseed();
3022 fSeedsPool->RegisterClone(seed);
3023 Double_t alpha=fSectors->GetAlpha(), shift=fSectors->GetAlphaShift();
3024 Double_t cs=cos(alpha), sn=sin(alpha);
3026 // Double_t x1 =fOuterSec->GetX(i1);
3027 //Double_t xx2=fOuterSec->GetX(i2);
3029 Double_t x1 =GetXrow(i1);
3030 Double_t xx2=GetXrow(i2);
3032 Double_t x3=GetX(), y3=GetY(), z3=GetZ();
3034 Int_t imiddle = (i2+i1)/2; //middle pad row index
3035 Double_t xm = GetXrow(imiddle); // radius of middle pad-row
3036 const AliTPCtrackerRow& krm=GetRow(sec,imiddle); //middle pad -row
3040 const AliTPCtrackerRow& kr1=GetRow(ns,i1);
3041 Double_t ymax = GetMaxY(i1)-kr1.GetDeadZone()-1.5;
3042 Double_t ymaxm = GetMaxY(imiddle)-kr1.GetDeadZone()-1.5;
3045 // change cut on curvature if it can't reach this layer
3046 // maximal curvature set to reach it
3047 Double_t dvertexmax = TMath::Sqrt((x1-x3)*(x1-x3)+(ymax+5-y3)*(ymax+5-y3));
3048 if (dvertexmax*0.5*cuts[0]>0.85){
3049 cuts[0] = 0.85/(dvertexmax*0.5+1.);
3051 Double_t r2min = 1/(cuts[0]*cuts[0]); //minimal square of radius given by cut
3054 if (deltay>0) ddsec = 0;
3055 // loop over clusters
3056 for (Int_t is=0; is < kr1; is++) {
3058 if (kr1[is]->IsUsed(10)) continue;
3059 Double_t y1=kr1[is]->GetY(), z1=kr1[is]->GetZ();
3060 //if (TMath::Abs(y1)>ymax) continue;
3062 if (deltay>0 && TMath::Abs(ymax-TMath::Abs(y1))> deltay ) continue; // seed only at the edge
3064 // find possible directions
3065 Float_t anglez = (z1-z3)/(x1-x3);
3066 Float_t extraz = z1 - anglez*(x1-xx2); // extrapolated z
3069 //find rotation angles relative to line given by vertex and point 1
3070 Double_t dvertex2 = (x1-x3)*(x1-x3)+(y1-y3)*(y1-y3);
3071 Double_t dvertex = TMath::Sqrt(dvertex2);
3072 Double_t angle13 = TMath::ATan((y1-y3)/(x1-x3));
3073 Double_t cs13 = cos(-angle13), sn13 = sin(-angle13);
3076 // loop over 2 sectors
3082 Double_t dddz1=0; // direction of delta inclination in z axis
3089 for (Int_t dsec = dsec1; dsec<=dsec2;dsec++){
3090 Int_t sec2 = sec + dsec;
3092 // AliTPCtrackerRow& kr2 = fOuterSec[(sec2+fkNOS)%fkNOS][i2];
3093 //AliTPCtrackerRow& kr2m = fOuterSec[(sec2+fkNOS)%fkNOS][imiddle];
3094 AliTPCtrackerRow& kr2 = GetRow((sec2+fkNOS)%fkNOS,i2);
3095 AliTPCtrackerRow& kr2m = GetRow((sec2+fkNOS)%fkNOS,imiddle);
3096 Int_t index1 = TMath::Max(kr2.Find(extraz-0.6-dddz1*TMath::Abs(z1)*0.05)-1,0);
3097 Int_t index2 = TMath::Min(kr2.Find(extraz+0.6+dddz2*TMath::Abs(z1)*0.05)+1,kr2);
3099 // rotation angles to p1-p3
3100 Double_t cs13r = cos(-angle13+dsec*alpha)/dvertex, sn13r = sin(-angle13+dsec*alpha)/dvertex;
3101 Double_t x2, y2, z2;
3103 // Double_t dymax = maxangle*TMath::Abs(x1-xx2);
3106 Double_t dxx0 = (xx2-x3)*cs13r;
3107 Double_t dyy0 = (xx2-x3)*sn13r;
3108 for (Int_t js=index1; js < index2; js++) {
3109 const AliTPCclusterMI *kcl = kr2[js];
3110 if (kcl->IsUsed(10)) continue;
3112 //calcutate parameters
3114 Double_t yy0 = dyy0 +(kcl->GetY()-y3)*cs13r;
3116 if (TMath::Abs(yy0)<0.000001) continue;
3117 Double_t xx0 = dxx0 -(kcl->GetY()-y3)*sn13r;
3118 Double_t y0 = 0.5*(xx0*xx0+yy0*yy0-xx0)/yy0;
3119 Double_t r02 = (0.25+y0*y0)*dvertex2;
3120 //curvature (radius) cut
3121 if (r02<r2min) continue;
3125 Double_t c0 = 1/TMath::Sqrt(r02);
3129 //Double_t dfi0 = 2.*TMath::ASin(dvertex*c0*0.5);
3130 //Double_t dfi1 = 2.*TMath::ASin(TMath::Sqrt(yy0*yy0+(1-xx0)*(1-xx0))*dvertex*c0*0.5);
3131 Double_t dfi0 = 2.*AliTPCFastMath::FastAsin(dvertex*c0*0.5);
3132 Double_t dfi1 = 2.*AliTPCFastMath::FastAsin(TMath::Sqrt(yy0*yy0+(1-xx0)*(1-xx0))*dvertex*c0*0.5);
3135 Double_t z0 = kcl->GetZ();
3136 Double_t zzzz2 = z1-(z1-z3)*dfi1/dfi0;
3137 if (TMath::Abs(zzzz2-z0)>0.5) continue;
3140 Double_t dip = (z1-z0)*c0/dfi1;
3141 Double_t x0 = (0.5*cs13+y0*sn13)*dvertex*c0;
3152 x2= xx2*cs-y2*sn*dsec;
3153 y2=+xx2*sn*dsec+y2*cs;
3163 // do we have cluster at the middle ?
3165 GetProlongation(x1,xm,x,ym,zm);
3167 AliTPCclusterMI * cm=0;
3168 if (TMath::Abs(ym)-ymaxm<0){
3169 cm = krm.FindNearest2(ym,zm,1.0,0.6,dummy);
3170 if ((!cm) || (cm->IsUsed(10))) {
3175 // rotate y1 to system 0
3176 // get state vector in rotated system
3177 Double_t yr1 = (-0.5*sn13+y0*cs13)*dvertex*c0;
3178 Double_t xr2 = x0*cs+yr1*sn*dsec;
3179 Double_t xr[5]={kcl->GetY(),kcl->GetZ(), xr2, dip, c0};
3181 GetProlongation(xx2,xm,xr,ym,zm);
3182 if (TMath::Abs(ym)-ymaxm<0){
3183 cm = kr2m.FindNearest2(ym,zm,1.0,0.6,dummy);
3184 if ((!cm) || (cm->IsUsed(10))) {
3194 dym = ym - cm->GetY();
3195 dzm = zm - cm->GetZ();
3202 Double_t sy1=kr1[is]->GetSigmaY2()*2., sz1=kr1[is]->GetSigmaZ2()*2.;
3203 Double_t sy2=kcl->GetSigmaY2()*2., sz2=kcl->GetSigmaZ2()*2.;
3204 //Double_t sy3=400*3./12., sy=0.1, sz=0.1;
3205 Double_t sy3=25000*x[4]*x[4]+0.1, sy=0.1, sz=0.1;
3206 //Double_t sy3=25000*x[4]*x[4]*60+0.5, sy=0.1, sz=0.1;
3208 Double_t f40=(F1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
3209 Double_t f42=(F1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
3210 Double_t f43=(F1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
3211 Double_t f20=(F2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
3212 Double_t f22=(F2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
3213 Double_t f23=(F2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
3215 Double_t f30=(F3(x1,y1+sy,x2,y2,z1,z2)-x[3])/sy;
3216 Double_t f31=(F3(x1,y1,x2,y2,z1+sz,z2)-x[3])/sz;
3217 Double_t f32=(F3(x1,y1,x2,y2+sy,z1,z2)-x[3])/sy;
3218 Double_t f34=(F3(x1,y1,x2,y2,z1,z2+sz)-x[3])/sz;
3222 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3223 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3224 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3225 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3226 c[13]=f30*sy1*f40+f32*sy2*f42;
3227 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3229 // if (!BuildSeed(kr1[is],kcl,0,x1,x2,x3,x,c)) continue;
3231 UInt_t index=kr1.GetIndex(is);
3232 if (seed) {fSeedsPool->MarkSlotFree(seed); seed = 0;}
3233 AliTPCseed *track = seed = new( fSeedsPool->NextFreeSlot() ) AliTPCseed(x1, ns*alpha+shift, x, c, index);
3234 fSeedsPool->RegisterClone(seed);
3235 track->SetIsSeeding(kTRUE);
3236 track->SetSeed1(i1);
3237 track->SetSeed2(i2);
3238 track->SetSeedType(3);
3242 FollowProlongation(*track, (i1+i2)/2,1);
3243 Int_t foundable,found,shared;
3244 track->GetClusterStatistic((i1+i2)/2,i1, found, foundable, shared, kTRUE);
3245 if ((found<0.55*foundable) || shared>0.5*found || (track->GetSigmaY2()+track->GetSigmaZ2())>0.5){
3246 fSeedsPool->MarkSlotFree(seed); seed = 0;
3252 FollowProlongation(*track, i2,1);
3256 track->SetBConstrain(1);
3257 // track->fLastPoint = i1+fInnerSec->GetNRows(); // first cluster in track position
3258 track->SetLastPoint(i1); // first cluster in track position
3259 track->SetFirstPoint(track->GetLastPoint());
3261 if (track->GetNumberOfClusters()<(i1-i2)*0.5 ||
3262 track->GetNumberOfClusters() < track->GetNFoundable()*0.6 ||
3263 track->GetNShared()>0.4*track->GetNumberOfClusters() ) {
3264 fSeedsPool->MarkSlotFree(seed); seed = 0;
3268 // Z VERTEX CONDITION
3269 Double_t zv, bz=GetBz();
3270 if ( !track->GetZAt(0.,bz,zv) ) continue;
3271 if (TMath::Abs(zv-z3)>cuts[2]) {
3272 FollowProlongation(*track, TMath::Max(i2-20,0));
3273 if ( !track->GetZAt(0.,bz,zv) ) continue;
3274 if (TMath::Abs(zv-z3)>cuts[2]){
3275 FollowProlongation(*track, TMath::Max(i2-40,0));
3276 if ( !track->GetZAt(0.,bz,zv) ) continue;
3277 if (TMath::Abs(zv-z3)>cuts[2] &&(track->GetNumberOfClusters() > track->GetNFoundable()*0.7)){
3278 // make seed without constrain
3279 AliTPCseed * track2 = MakeSeed(track,0.2,0.5,1.);
3280 FollowProlongation(*track2, i2,1);
3281 track2->SetBConstrain(kFALSE);
3282 track2->SetSeedType(1);
3283 arr->AddLast(track2);
3284 fSeedsPool->MarkSlotFree( seed ); seed = 0;
3288 fSeedsPool->MarkSlotFree( seed ); seed = 0;
3295 track->SetSeedType(0);
3296 arr->AddLast(track); // note, track is seed, don't free the seed
3297 seed = new( fSeedsPool->NextFreeSlot() ) AliTPCseed;
3298 fSeedsPool->RegisterClone(seed);
3300 // don't consider other combinations
3301 if (track->GetNumberOfClusters() > track->GetNFoundable()*0.8)
3307 Info("MakeSeeds3","\nSeeding statistic:\t%d\t%d\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin,nout1,nout2);
3309 if (seed) fSeedsPool->MarkSlotFree( seed );
3313 void AliTPCtrackerMI::MakeSeeds5(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, Float_t cuts[4],
3318 //-----------------------------------------------------------------
3319 // This function creates track seeds.
3320 //-----------------------------------------------------------------
3321 // cuts[0] - fP4 cut
3322 // cuts[1] - tan(phi) cut
3323 // cuts[2] - zvertex cut
3324 // cuts[3] - fP3 cut
3334 Double_t x[5], c[15];
3336 // make temporary seed
3337 AliTPCseed * seed = new( fSeedsPool->NextFreeSlot() ) AliTPCseed;
3338 fSeedsPool->RegisterClone(seed);
3339 Double_t alpha=fOuterSec->GetAlpha(), shift=fOuterSec->GetAlphaShift();
3340 // Double_t cs=cos(alpha), sn=sin(alpha);
3345 Double_t x1 = GetXrow(i1-1);
3346 const AliTPCtrackerRow& kr1=GetRow(sec,i1-1);
3347 Double_t y1max = GetMaxY(i1-1)-kr1.GetDeadZone()-1.5;
3349 Double_t x1p = GetXrow(i1);
3350 const AliTPCtrackerRow& kr1p=GetRow(sec,i1);
3352 Double_t x1m = GetXrow(i1-2);
3353 const AliTPCtrackerRow& kr1m=GetRow(sec,i1-2);
3356 //last 3 padrow for seeding
3357 AliTPCtrackerRow& kr3 = GetRow((sec+fkNOS)%fkNOS,i1-7);
3358 Double_t x3 = GetXrow(i1-7);
3359 // Double_t y3max= GetMaxY(i1-7)-kr3.fDeadZone-1.5;
3361 AliTPCtrackerRow& kr3p = GetRow((sec+fkNOS)%fkNOS,i1-6);
3362 Double_t x3p = GetXrow(i1-6);
3364 AliTPCtrackerRow& kr3m = GetRow((sec+fkNOS)%fkNOS,i1-8);
3365 Double_t x3m = GetXrow(i1-8);
3370 Int_t im = i1-4; //middle pad row index
3371 Double_t xm = GetXrow(im); // radius of middle pad-row
3372 const AliTPCtrackerRow& krm=GetRow(sec,im); //middle pad -row
3373 // Double_t ymmax = GetMaxY(im)-kr1.fDeadZone-1.5;
3376 Double_t deltax = x1-x3;
3377 Double_t dymax = deltax*cuts[1];
3378 Double_t dzmax = deltax*cuts[3];
3380 // loop over clusters
3381 for (Int_t is=0; is < kr1; is++) {
3383 if (kr1[is]->IsUsed(10)) continue;
3384 Double_t y1=kr1[is]->GetY(), z1=kr1[is]->GetZ();
3386 if (deltay>0 && TMath::Abs(y1max-TMath::Abs(y1))> deltay ) continue; // seed only at the edge
3388 Int_t index1 = TMath::Max(kr3.Find(z1-dzmax)-1,0);
3389 Int_t index2 = TMath::Min(kr3.Find(z1+dzmax)+1,kr3);
3395 for (Int_t js=index1; js < index2; js++) {
3396 const AliTPCclusterMI *kcl = kr3[js];
3397 if (kcl->IsUsed(10)) continue;
3399 // apply angular cuts
3400 if (TMath::Abs(y1-y3)>dymax) continue;
3403 if (TMath::Abs(z1-z3)>dzmax) continue;
3405 Double_t angley = (y1-y3)/(x1-x3);
3406 Double_t anglez = (z1-z3)/(x1-x3);
3408 Double_t erry = TMath::Abs(angley)*(x1-x1m)*0.5+0.5;
3409 Double_t errz = TMath::Abs(anglez)*(x1-x1m)*0.5+0.5;
3411 Double_t yyym = angley*(xm-x1)+y1;
3412 Double_t zzzm = anglez*(xm-x1)+z1;
3414 const AliTPCclusterMI *kcm = krm.FindNearest2(yyym,zzzm,erry,errz,index);
3416 if (kcm->IsUsed(10)) continue;
3418 erry = TMath::Abs(angley)*(x1-x1m)*0.4+0.5;
3419 errz = TMath::Abs(anglez)*(x1-x1m)*0.4+0.5;
3426 // look around first
3427 const AliTPCclusterMI *kc1m = kr1m.FindNearest2(angley*(x1m-x1)+y1,
3433 if (kc1m->IsUsed(10)) used++;
3435 const AliTPCclusterMI *kc1p = kr1p.FindNearest2(angley*(x1p-x1)+y1,
3441 if (kc1p->IsUsed(10)) used++;
3443 if (used>1) continue;
3444 if (found<1) continue;
3448 const AliTPCclusterMI *kc3m = kr3m.FindNearest2(angley*(x3m-x3)+y3,
3454 if (kc3m->IsUsed(10)) used++;
3458 const AliTPCclusterMI *kc3p = kr3p.FindNearest2(angley*(x3p-x3)+y3,
3464 if (kc3p->IsUsed(10)) used++;
3468 if (used>1) continue;
3469 if (found<3) continue;
3479 x[4]=F1(x1,y1,x2,y2,x3,y3);
3480 //if (TMath::Abs(x[4]) >= cuts[0]) continue;
3483 x[2]=F2(x1,y1,x2,y2,x3,y3);
3486 x[3]=F3n(x1,y1,x2,y2,z1,z2,x[4]);
3487 //if (TMath::Abs(x[3]) > cuts[3]) continue;
3491 Double_t sy1=0.1, sz1=0.1;
3492 Double_t sy2=0.1, sz2=0.1;
3493 Double_t sy3=0.1, sy=0.1, sz=0.1;
3495 Double_t f40=(F1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
3496 Double_t f42=(F1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
3497 Double_t f43=(F1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
3498 Double_t f20=(F2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
3499 Double_t f22=(F2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
3500 Double_t f23=(F2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
3502 Double_t f30=(F3(x1,y1+sy,x2,y2,z1,z2)-x[3])/sy;
3503 Double_t f31=(F3(x1,y1,x2,y2,z1+sz,z2)-x[3])/sz;
3504 Double_t f32=(F3(x1,y1,x2,y2+sy,z1,z2)-x[3])/sy;
3505 Double_t f34=(F3(x1,y1,x2,y2,z1,z2+sz)-x[3])/sz;
3509 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3510 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3511 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3512 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3513 c[13]=f30*sy1*f40+f32*sy2*f42;
3514 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3516 // if (!BuildSeed(kr1[is],kcl,0,x1,x2,x3,x,c)) continue;
3518 index=kr1.GetIndex(is);
3519 if (seed) {fSeedsPool->MarkSlotFree( seed ); seed = 0;}
3520 AliTPCseed *track = seed = new( fSeedsPool->NextFreeSlot() ) AliTPCseed(x1, sec*alpha+shift, x, c, index);
3521 fSeedsPool->RegisterClone(seed);
3523 track->SetIsSeeding(kTRUE);
3526 FollowProlongation(*track, i1-7,1);
3527 if (track->GetNumberOfClusters() < track->GetNFoundable()*0.75 ||
3528 track->GetNShared()>0.6*track->GetNumberOfClusters() || ( track->GetSigmaY2()+ track->GetSigmaZ2())>0.6){
3529 fSeedsPool->MarkSlotFree( seed ); seed = 0;
3535 FollowProlongation(*track, i2,1);
3536 track->SetBConstrain(0);
3537 track->SetLastPoint(i1+fInnerSec->GetNRows()); // first cluster in track position
3538 track->SetFirstPoint(track->GetLastPoint());
3540 if (track->GetNumberOfClusters()<(i1-i2)*0.5 ||
3541 track->GetNumberOfClusters()<track->GetNFoundable()*0.7 ||
3542 track->GetNShared()>2. || track->GetChi2()/track->GetNumberOfClusters()>6 || ( track->GetSigmaY2()+ track->GetSigmaZ2())>0.5 ) {
3543 fSeedsPool->MarkSlotFree( seed ); seed = 0;
3548 FollowProlongation(*track, TMath::Max(i2-10,0),1);
3549 AliTPCseed * track2 = MakeSeed(track,0.2,0.5,0.9);
3550 FollowProlongation(*track2, i2,1);
3551 track2->SetBConstrain(kFALSE);
3552 track2->SetSeedType(4);
3553 arr->AddLast(track2);
3554 fSeedsPool->MarkSlotFree( seed ); seed = 0;
3558 //arr->AddLast(track);
3559 //seed = new AliTPCseed;
3565 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);
3567 if (seed) fSeedsPool->MarkSlotFree(seed);
3571 //_____________________________________________________________________________
3572 void AliTPCtrackerMI::MakeSeeds2(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, Float_t */*cuts[4]*/,
3573 Float_t deltay, Bool_t /*bconstrain*/) {
3574 //-----------------------------------------------------------------
3575 // This function creates track seeds - without vertex constraint
3576 //-----------------------------------------------------------------
3577 // cuts[0] - fP4 cut - not applied
3578 // cuts[1] - tan(phi) cut
3579 // cuts[2] - zvertex cut - not applied
3580 // cuts[3] - fP3 cut
3590 Double_t alpha=fOuterSec->GetAlpha(), shift=fOuterSec->GetAlphaShift();
3591 // Double_t cs=cos(alpha), sn=sin(alpha);
3592 Int_t row0 = (i1+i2)/2;
3593 Int_t drow = (i1-i2)/2;
3594 const AliTPCtrackerRow& kr0=fSectors[sec][row0];
3595 AliTPCtrackerRow * kr=0;
3597 AliTPCpolyTrack polytrack;
3598 Int_t nclusters=fSectors[sec][row0];
3599 AliTPCseed * seed = new( fSeedsPool->NextFreeSlot() ) AliTPCseed;
3600 fSeedsPool->RegisterClone(seed);
3605 for (Int_t is=0; is < nclusters; is++) { //LOOP over clusters
3607 Int_t nfoundable =0;
3608 for (Int_t iter =1; iter<2; iter++){ //iterations
3609 const AliTPCtrackerRow& krm=fSectors[sec][row0-iter];
3610 const AliTPCtrackerRow& krp=fSectors[sec][row0+iter];
3611 const AliTPCclusterMI * cl= kr0[is];
3613 if (cl->IsUsed(10)) {
3619 Double_t x = kr0.GetX();
3620 // Initialization of the polytrack
3625 Double_t y0= cl->GetY();
3626 Double_t z0= cl->GetZ();
3630 Double_t ymax = fSectors->GetMaxY(row0)-kr0.GetDeadZone()-1.5;
3631 if (deltay>0 && TMath::Abs(ymax-TMath::Abs(y0))> deltay ) continue; // seed only at the edge
3633 erry = (0.5)*cl->GetSigmaY2()/TMath::Sqrt(cl->GetQ())*6;
3634 errz = (0.5)*cl->GetSigmaZ2()/TMath::Sqrt(cl->GetQ())*6;
3635 polytrack.AddPoint(x,y0,z0,erry, errz);
3638 if (cl->IsUsed(10)) sumused++;
3641 Float_t roady = (5*TMath::Sqrt(cl->GetSigmaY2()+0.2)+1.)*iter;
3642 Float_t roadz = (5*TMath::Sqrt(cl->GetSigmaZ2()+0.2)+1.)*iter;
3645 AliTPCclusterMI * cl1 = krm.FindNearest(y0,z0,roady,roadz);
3646 if (cl1 && TMath::Abs(ymax-TMath::Abs(y0))) {
3647 erry = (0.5)*cl1->GetSigmaY2()/TMath::Sqrt(cl1->GetQ())*3;
3648 errz = (0.5)*cl1->GetSigmaZ2()/TMath::Sqrt(cl1->GetQ())*3;
3649 if (cl1->IsUsed(10)) sumused++;
3650 polytrack.AddPoint(x,cl1->GetY(),cl1->GetZ(),erry,errz);
3654 AliTPCclusterMI * cl2 = krp.FindNearest(y0,z0,roady,roadz);
3656 erry = (0.5)*cl2->GetSigmaY2()/TMath::Sqrt(cl2->GetQ())*3;
3657 errz = (0.5)*cl2->GetSigmaZ2()/TMath::Sqrt(cl2->GetQ())*3;
3658 if (cl2->IsUsed(10)) sumused++;
3659 polytrack.AddPoint(x,cl2->GetY(),cl2->GetZ(),erry,errz);
3662 if (sumused>0) continue;
3664 polytrack.UpdateParameters();
3670 nfoundable = polytrack.GetN();
3671 nfound = nfoundable;
3673 for (Int_t ddrow = iter+1; ddrow<drow;ddrow++){
3674 Float_t maxdist = 0.8*(1.+3./(ddrow));
3675 for (Int_t delta = -1;delta<=1;delta+=2){
3676 Int_t row = row0+ddrow*delta;
3677 kr = &(fSectors[sec][row]);
3678 Double_t xn = kr->GetX();
3679 Double_t ymax1 = fSectors->GetMaxY(row)-kr->GetDeadZone()-1.5;
3680 polytrack.GetFitPoint(xn,yn,zn);
3681 if (TMath::Abs(yn)>ymax1) continue;
3683 AliTPCclusterMI * cln = kr->FindNearest(yn,zn,roady,roadz);
3685 Float_t dist = TMath::Sqrt( (yn-cln->GetY())*(yn-cln->GetY())+(zn-cln->GetZ())*(zn-cln->GetZ()));
3688 erry = (dist+0.3)*cln->GetSigmaY2()/TMath::Sqrt(cln->GetQ())*(1.+1./(ddrow));
3689 errz = (dist+0.3)*cln->GetSigmaZ2()/TMath::Sqrt(cln->GetQ())*(1.+1./(ddrow));
3690 if (cln->IsUsed(10)) {
3691 // printf("used\n");
3699 polytrack.AddPoint(xn,cln->GetY(),cln->GetZ(),erry, errz);
3704 if ( (sumused>3) || (sumused>0.5*nfound) || (nfound<0.6*nfoundable)) break;
3705 polytrack.UpdateParameters();
3708 if ( (sumused>3) || (sumused>0.5*nfound)) {
3709 //printf("sumused %d\n",sumused);
3714 polytrack.GetFitDerivation(kr0.GetX(),dy,dz);
3715 AliTPCpolyTrack track2;
3717 polytrack.Refit(track2,0.5+TMath::Abs(dy)*0.3,0.4+TMath::Abs(dz)*0.3);
3718 if (track2.GetN()<0.5*nfoundable) continue;
3721 if ((nfound>0.6*nfoundable) &&( nfoundable>0.4*(i1-i2))) {
3723 // test seed with and without constrain
3724 for (Int_t constrain=0; constrain<=0;constrain++){
3725 // add polytrack candidate
3727 Double_t x[5], c[15];
3728 Double_t x1,x2,x3,y1,y2,y3,z1,z2,z3;
3729 track2.GetBoundaries(x3,x1);
3731 track2.GetFitPoint(x1,y1,z1);
3732 track2.GetFitPoint(x2,y2,z2);
3733 track2.GetFitPoint(x3,y3,z3);
3735 //is track pointing to the vertex ?
3738 polytrack.GetFitPoint(x0,y0,z0);
3751 x[4]=F1(x1,y1,x2,y2,x3,y3);
3753 // if (TMath::Abs(x[4]) >= cuts[0]) continue; //
3754 x[2]=F2(x1,y1,x2,y2,x3,y3);
3756 //if (TMath::Abs(x[4]*x1-x[2]) >= cuts[1]) continue;
3757 //x[3]=F3(x1,y1,x2,y2,z1,z2);
3758 x[3]=F3n(x1,y1,x3,y3,z1,z3,x[4]);
3759 //if (TMath::Abs(x[3]) > cuts[3]) continue;
3762 Double_t sy =0.1, sz =0.1;
3763 Double_t sy1=0.02, sz1=0.02;
3764 Double_t sy2=0.02, sz2=0.02;
3768 sy3=25000*x[4]*x[4]+0.1, sy=0.1, sz=0.1;
3771 Double_t f40=(F1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
3772 Double_t f42=(F1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
3773 Double_t f43=(F1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
3774 Double_t f20=(F2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
3775 Double_t f22=(F2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
3776 Double_t f23=(F2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
3778 Double_t f30=(F3(x1,y1+sy,x3,y3,z1,z3)-x[3])/sy;
3779 Double_t f31=(F3(x1,y1,x3,y3,z1+sz,z3)-x[3])/sz;
3780 Double_t f32=(F3(x1,y1,x3,y3+sy,z1,z3)-x[3])/sy;
3781 Double_t f34=(F3(x1,y1,x3,y3,z1,z3+sz)-x[3])/sz;
3786 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3787 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3788 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3789 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3790 c[13]=f30*sy1*f40+f32*sy2*f42;
3791 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3793 //Int_t row1 = fSectors->GetRowNumber(x1);
3794 Int_t row1 = GetRowNumber(x1);
3798 if (seed) {fSeedsPool->MarkSlotFree( seed ); seed = 0;}
3799 AliTPCseed *track = seed = new( fSeedsPool->NextFreeSlot() ) AliTPCseed(x1,sec*alpha+shift,x,c,index);
3800 fSeedsPool->RegisterClone(seed);
3801 track->SetIsSeeding(kTRUE);
3802 Int_t rc=FollowProlongation(*track, i2);
3803 if (constrain) track->SetBConstrain(1);
3805 track->SetBConstrain(0);
3806 track->SetLastPoint(row1+fInnerSec->GetNRows()); // first cluster in track position
3807 track->SetFirstPoint(track->GetLastPoint());
3809 if (rc==0 || track->GetNumberOfClusters()<(i1-i2)*0.5 ||
3810 track->GetNumberOfClusters() < track->GetNFoundable()*0.6 ||
3811 track->GetNShared()>0.4*track->GetNumberOfClusters()) {
3812 fSeedsPool->MarkSlotFree( seed ); seed = 0;
3815 arr->AddLast(track); // track IS seed, don't free seed
3816 seed = new( fSeedsPool->NextFreeSlot() ) AliTPCseed;
3817 fSeedsPool->RegisterClone(seed);
3821 } // if accepted seed
3824 Info("MakeSeeds2","\nSeeding statiistic:\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin3);
3826 if (seed) fSeedsPool->MarkSlotFree( seed );
3830 AliTPCseed *AliTPCtrackerMI::MakeSeed(AliTPCseed *const track, Float_t r0, Float_t r1, Float_t r2)
3834 //reseed using track points
3835 Int_t p0 = int(r0*track->GetNumberOfClusters()); // point 0
3836 Int_t p1 = int(r1*track->GetNumberOfClusters());
3837 Int_t p2 = int(r2*track->GetNumberOfClusters()); // last point
3839 Double_t x0[3],x1[3],x2[3];
3840 for (Int_t i=0;i<3;i++){
3846 // find track position at given ratio of the length
3847 Int_t sec0=0, sec1=0, sec2=0;
3850 for (Int_t i=0;i<160;i++){
3851 if (track->GetClusterPointer(i)){
3853 AliTPCTrackerPoint *trpoint =track->GetTrackPoint(i);
3854 if ( (index<p0) || x0[0]<0 ){
3855 if (trpoint->GetX()>1){
3856 clindex = track->GetClusterIndex2(i);
3858 x0[0] = trpoint->GetX();
3859 x0[1] = trpoint->GetY();
3860 x0[2] = trpoint->GetZ();
3861 sec0 = ((clindex&0xff000000)>>24)%18;
3866 if ( (index<p1) &&(trpoint->GetX()>1)){
3867 clindex = track->GetClusterIndex2(i);
3869 x1[0] = trpoint->GetX();
3870 x1[1] = trpoint->GetY();
3871 x1[2] = trpoint->GetZ();
3872 sec1 = ((clindex&0xff000000)>>24)%18;
3875 if ( (index<p2) &&(trpoint->GetX()>1)){
3876 clindex = track->GetClusterIndex2(i);
3878 x2[0] = trpoint->GetX();
3879 x2[1] = trpoint->GetY();
3880 x2[2] = trpoint->GetZ();
3881 sec2 = ((clindex&0xff000000)>>24)%18;
3888 Double_t alpha, cs,sn, xx2,yy2;
3890 alpha = (sec1-sec2)*fSectors->GetAlpha();
3891 cs = TMath::Cos(alpha);
3892 sn = TMath::Sin(alpha);
3893 xx2= x1[0]*cs-x1[1]*sn;
3894 yy2= x1[0]*sn+x1[1]*cs;
3898 alpha = (sec0-sec2)*fSectors->GetAlpha();
3899 cs = TMath::Cos(alpha);
3900 sn = TMath::Sin(alpha);
3901 xx2= x0[0]*cs-x0[1]*sn;
3902 yy2= x0[0]*sn+x0[1]*cs;
3908 Double_t x[5],c[15];
3912 x[4]=F1(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]);
3913 // if (x[4]>1) return 0;
3914 x[2]=F2(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]);
3915 x[3]=F3n(x2[0],x2[1],x0[0],x0[1],x2[2],x0[2],x[4]);
3916 //if (TMath::Abs(x[3]) > 2.2) return 0;
3917 //if (TMath::Abs(x[2]) > 1.99) return 0;
3919 Double_t sy =0.1, sz =0.1;
3921 Double_t sy1=0.02+track->GetSigmaY2(), sz1=0.02+track->GetSigmaZ2();
3922 Double_t sy2=0.01+track->GetSigmaY2(), sz2=0.01+track->GetSigmaZ2();
3923 Double_t sy3=0.01+track->GetSigmaY2();
3925 Double_t f40=(F1(x2[0],x2[1]+sy,x1[0],x1[1],x0[0],x0[1])-x[4])/sy;
3926 Double_t f42=(F1(x2[0],x2[1],x1[0],x1[1]+sy,x0[0],x0[1])-x[4])/sy;
3927 Double_t f43=(F1(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]+sy)-x[4])/sy;
3928 Double_t f20=(F2(x2[0],x2[1]+sy,x1[0],x1[1],x0[0],x0[1])-x[2])/sy;
3929 Double_t f22=(F2(x2[0],x2[1],x1[0],x1[1]+sy,x0[0],x0[1])-x[2])/sy;
3930 Double_t f23=(F2(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]+sy)-x[2])/sy;
3932 Double_t f30=(F3(x2[0],x2[1]+sy,x0[0],x0[1],x2[2],x0[2])-x[3])/sy;
3933 Double_t f31=(F3(x2[0],x2[1],x0[0],x0[1],x2[2]+sz,x0[2])-x[3])/sz;
3934 Double_t f32=(F3(x2[0],x2[1],x0[0],x0[1]+sy,x2[2],x0[2])-x[3])/sy;
3935 Double_t f34=(F3(x2[0],x2[1],x0[0],x0[1],x2[2],x0[2]+sz)-x[3])/sz;
3940 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3941 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3942 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3943 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3944 c[13]=f30*sy1*f40+f32*sy2*f42;
3945 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3947 // Int_t row1 = fSectors->GetRowNumber(x2[0]);
3948 AliTPCseed *seed = new( fSeedsPool->NextFreeSlot() ) AliTPCseed(x2[0], sec2*fSectors->GetAlpha()+fSectors->GetAlphaShift(), x, c, 0);
3949 fSeedsPool->RegisterClone(seed);
3950 // Double_t y0,z0,y1,z1, y2,z2;
3951 //seed->GetProlongation(x0[0],y0,z0);
3952 // seed->GetProlongation(x1[0],y1,z1);
3953 //seed->GetProlongation(x2[0],y2,z2);
3955 seed->SetLastPoint(pp2);
3956 seed->SetFirstPoint(pp2);
3963 AliTPCseed *AliTPCtrackerMI::ReSeed(const AliTPCseed *track, Float_t r0, Float_t r1, Float_t r2)
3967 //reseed using founded clusters
3969 // Find the number of clusters
3970 Int_t nclusters = 0;
3971 for (Int_t irow=0;irow<160;irow++){
3972 if (track->GetClusterIndex(irow)>0) nclusters++;
3976 ipos[0] = TMath::Max(int(r0*nclusters),0); // point 0 cluster
3977 ipos[1] = TMath::Min(int(r1*nclusters),nclusters-1); //
3978 ipos[2] = TMath::Min(int(r2*nclusters),nclusters-1); // last point
3981 Double_t xyz[3][3]={{0}};
3982 Int_t row[3]={0},sec[3]={0,0,0};
3984 // find track row position at given ratio of the length
3986 for (Int_t irow=0;irow<160;irow++){
3987 if (track->GetClusterIndex2(irow)<0) continue;
3989 for (Int_t ipoint=0;ipoint<3;ipoint++){
3990 if (index<=ipos[ipoint]) row[ipoint] = irow;
3994 //Get cluster and sector position
3995 for (Int_t ipoint=0;ipoint<3;ipoint++){
3996 Int_t clindex = track->GetClusterIndex2(row[ipoint]);
3997 AliTPCclusterMI * cl = GetClusterMI(clindex);
4000 // AliTPCclusterMI * cl = GetClusterMI(clindex);
4003 sec[ipoint] = ((clindex&0xff000000)>>24)%18;
4004 xyz[ipoint][0] = GetXrow(row[ipoint]);
4005 xyz[ipoint][1] = cl->GetY();
4006 xyz[ipoint][2] = cl->GetZ();
4010 // Calculate seed state vector and covariance matrix
4012 Double_t alpha, cs,sn, xx2,yy2;
4014 alpha = (sec[1]-sec[2])*fSectors->GetAlpha();
4015 cs = TMath::Cos(alpha);
4016 sn = TMath::Sin(alpha);
4017 xx2= xyz[1][0]*cs-xyz[1][1]*sn;
4018 yy2= xyz[1][0]*sn+xyz[1][1]*cs;
4022 alpha = (sec[0]-sec[2])*fSectors->GetAlpha();
4023 cs = TMath::Cos(alpha);
4024 sn = TMath::Sin(alpha);
4025 xx2= xyz[0][0]*cs-xyz[0][1]*sn;
4026 yy2= xyz[0][0]*sn+xyz[0][1]*cs;
4032 Double_t x[5],c[15];
4036 x[4]=F1(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
4037 x[2]=F2(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
4038 x[3]=F3n(xyz[2][0],xyz[2][1],xyz[0][0],xyz[0][1],xyz[2][2],xyz[0][2],x[4]);
4040 Double_t sy =0.1, sz =0.1;
4042 Double_t sy1=0.2, sz1=0.2;
4043 Double_t sy2=0.2, sz2=0.2;
4046 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;
4047 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;
4048 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;
4049 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;
4050 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;
4051 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;
4053 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;
4054 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;
4055 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;
4056 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;
4061 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
4062 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
4063 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
4064 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
4065 c[13]=f30*sy1*f40+f32*sy2*f42;
4066 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
4068 // Int_t row1 = fSectors->GetRowNumber(xyz[2][0]);
4069 AliTPCseed *seed=new( fSeedsPool->NextFreeSlot() ) AliTPCseed(xyz[2][0], sec[2]*fSectors->GetAlpha()+fSectors->GetAlphaShift(), x, c, 0);
4070 fSeedsPool->RegisterClone(seed);
4071 seed->SetLastPoint(row[2]);
4072 seed->SetFirstPoint(row[2]);
4077 AliTPCseed *AliTPCtrackerMI::ReSeed(AliTPCseed *track,Int_t r0, Bool_t forward)
4081 //reseed using founded clusters
4084 Int_t row[3]={0,0,0};
4085 Int_t sec[3]={0,0,0};
4087 // forward direction
4089 for (Int_t irow=r0;irow<160;irow++){
4090 if (track->GetClusterIndex(irow)>0){
4095 for (Int_t irow=160;irow>r0;irow--){
4096 if (track->GetClusterIndex(irow)>0){
4101 for (Int_t irow=row[2]-15;irow>row[0];irow--){
4102 if (track->GetClusterIndex(irow)>0){
4110 for (Int_t irow=0;irow<r0;irow++){
4111 if (track->GetClusterIndex(irow)>0){
4116 for (Int_t irow=r0;irow>0;irow--){
4117 if (track->GetClusterIndex(irow)>0){
4122 for (Int_t irow=row[2]-15;irow>row[0];irow--){
4123 if (track->GetClusterIndex(irow)>0){
4130 if ((row[2]-row[0])<20) return 0;
4131 if (row[1]==0) return 0;
4134 //Get cluster and sector position
4135 for (Int_t ipoint=0;ipoint<3;ipoint++){
4136 Int_t clindex = track->GetClusterIndex2(row[ipoint]);
4137 AliTPCclusterMI * cl = GetClusterMI(clindex);
4140 // AliTPCclusterMI * cl = GetClusterMI(clindex);
4143 sec[ipoint] = ((clindex&0xff000000)>>24)%18;
4144 xyz[ipoint][0] = GetXrow(row[ipoint]);
4145 AliTPCTrackerPoint * point = track->GetTrackPoint(row[ipoint]);
4146 if (point&&ipoint<2){
4148 xyz[ipoint][1] = point->GetY();
4149 xyz[ipoint][2] = point->GetZ();
4152 xyz[ipoint][1] = cl->GetY();
4153 xyz[ipoint][2] = cl->GetZ();
4160 // Calculate seed state vector and covariance matrix
4162 Double_t alpha, cs,sn, xx2,yy2;
4164 alpha = (sec[1]-sec[2])*fSectors->GetAlpha();
4165 cs = TMath::Cos(alpha);
4166 sn = TMath::Sin(alpha);
4167 xx2= xyz[1][0]*cs-xyz[1][1]*sn;
4168 yy2= xyz[1][0]*sn+xyz[1][1]*cs;
4172 alpha = (sec[0]-sec[2])*fSectors->GetAlpha();
4173 cs = TMath::Cos(alpha);
4174 sn = TMath::Sin(alpha);
4175 xx2= xyz[0][0]*cs-xyz[0][1]*sn;
4176 yy2= xyz[0][0]*sn+xyz[0][1]*cs;
4182 Double_t x[5],c[15];
4186 x[4]=F1(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
4187 x[2]=F2(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
4188 x[3]=F3n(xyz[2][0],xyz[2][1],xyz[0][0],xyz[0][1],xyz[2][2],xyz[0][2],x[4]);
4190 Double_t sy =0.1, sz =0.1;
4192 Double_t sy1=0.2, sz1=0.2;
4193 Double_t sy2=0.2, sz2=0.2;
4196 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;
4197 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;
4198 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;
4199 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;
4200 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;
4201 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;
4203 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;
4204 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;
4205 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;
4206 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;
4211 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
4212 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
4213 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
4214 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
4215 c[13]=f30*sy1*f40+f32*sy2*f42;
4216 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
4218 // Int_t row1 = fSectors->GetRowNumber(xyz[2][0]);
4219 AliTPCseed *seed=new( fSeedsPool->NextFreeSlot() ) AliTPCseed(xyz[2][0], sec[2]*fSectors->GetAlpha()+fSectors->GetAlphaShift(), x, c, 0);
4220 fSeedsPool->RegisterClone(seed);
4221 seed->SetLastPoint(row[2]);
4222 seed->SetFirstPoint(row[2]);
4223 for (Int_t i=row[0];i<row[2];i++){
4224 seed->SetClusterIndex(i, track->GetClusterIndex(i));
4232 void AliTPCtrackerMI::FindMultiMC(const TObjArray * array, AliESDEvent */*esd*/, Int_t iter)
4235 // find multi tracks - THIS FUNCTION IS ONLY FOR DEBUG PURPOSES
4237 // Use AliTPCReconstructor::StreamLevel()>2 if you want to tune parameters - cuts
4239 // Two reasons to have multiple find tracks
4240 // 1. Curling tracks can be find more than once
4241 // 2. Splitted tracks
4242 // a.) Multiple seeding to increase tracking efficiency - (~ 100% reached)
4243 // b.) Edge effect on the sector boundaries
4246 // Algorithm done in 2 phases - because of CPU consumption
4247 // it is n^2 algorithm - for lead-lead 20000x20000 combination are investigated
4249 // Algorihm for curling tracks sign:
4250 // 1 phase -makes a very rough fast cuts to minimize combinatorics
4251 // a.) opposite sign
4252 // b.) one of the tracks - not pointing to the primary vertex -
4253 // c.) delta tan(theta)
4255 // 2 phase - calculates DCA between tracks - time consument
4260 // General cuts - for splitted tracks and for curling tracks
4262 const Float_t kMaxdPhi = 0.2; // maximal distance in phi
4264 // Curling tracks cuts
4269 Int_t nentries = array->GetEntriesFast();
4270 AliHelix *helixes = new AliHelix[nentries];
4271 Float_t *xm = new Float_t[nentries];
4272 Float_t *dz0 = new Float_t[nentries];
4273 Float_t *dz1 = new Float_t[nentries];
4279 // Find track COG in x direction - point with best defined parameters
4281 for (Int_t i=0;i<nentries;i++){
4282 AliTPCseed* track = (AliTPCseed*)array->At(i);
4283 if (!track) continue;
4284 track->SetCircular(0);
4285 new (&helixes[i]) AliHelix(*track);
4289 track->GetDZ(GetX(),GetY(),GetZ(),GetBz(),dz);
4292 for (Int_t icl=0; icl<160; icl++){
4293 AliTPCclusterMI * cl = track->GetClusterPointer(icl);
4299 if (ncl>0) xm[i]/=Float_t(ncl);
4302 for (Int_t i0=0;i0<nentries;i0++){
4303 AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
4304 if (!track0) continue;
4305 Float_t xc0 = helixes[i0].GetHelix(6);
4306 Float_t yc0 = helixes[i0].GetHelix(7);
4307 Float_t r0 = helixes[i0].GetHelix(8);
4308 Float_t rc0 = TMath::Sqrt(xc0*xc0+yc0*yc0);
4309 Float_t fi0 = TMath::ATan2(yc0,xc0);
4311 for (Int_t i1=i0+1;i1<nentries;i1++){
4312 AliTPCseed * track1 = (AliTPCseed*)array->At(i1);
4313 if (!track1) continue;
4314 Int_t lab0=track0->GetLabel();
4315 Int_t lab1=track1->GetLabel();
4316 if (TMath::Abs(lab0)!=TMath::Abs(lab1)) continue;
4318 Float_t xc1 = helixes[i1].GetHelix(6);
4319 Float_t yc1 = helixes[i1].GetHelix(7);
4320 Float_t r1 = helixes[i1].GetHelix(8);
4321 Float_t rc1 = TMath::Sqrt(xc1*xc1+yc1*yc1);
4322 Float_t fi1 = TMath::ATan2(yc1,xc1);
4324 Float_t dfi = fi0-fi1;
4327 if (dfi>1.5*TMath::Pi()) dfi-=TMath::Pi(); // take care about edge effect
4328 if (dfi<-1.5*TMath::Pi()) dfi+=TMath::Pi(); //
4329 if (TMath::Abs(dfi)>kMaxdPhi&&helixes[i0].GetHelix(4)*helixes[i1].GetHelix(4)<0){
4331 // if short tracks with undefined sign
4332 fi1 = -TMath::ATan2(yc1,-xc1);
4335 Float_t dtheta = TMath::Abs(track0->GetTgl()-track1->GetTgl())<TMath::Abs(track0->GetTgl()+track1->GetTgl())? track0->GetTgl()-track1->GetTgl():track0->GetTgl()+track1->GetTgl();
4338 // debug stream to tune "fast cuts"
4340 Double_t dist[3]; // distance at X
4341 Double_t mdist[3]={0,0,0}; // mean distance X+-40cm
4342 track0->GetDistance(track1,0.5*(xm[i0]+xm[i1])-40.,dist,AliTracker::GetBz());
4343 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[i]);
4344 track0->GetDistance(track1,0.5*(xm[i0]+xm[i1])+40.,dist,AliTracker::GetBz());
4345 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[i]);
4346 track0->GetDistance(track1,0.5*(xm[i0]+xm[i1]),dist,AliTracker::GetBz());
4347 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[i]);
4348 for (Int_t i=0;i<3;i++) mdist[i]*=0.33333;
4352 for (Int_t icl=0; icl<160; icl++){
4353 AliTPCclusterMI * cl0 = track0->GetClusterPointer(icl);
4354 AliTPCclusterMI * cl1 = track1->GetClusterPointer(icl);
4357 if (cl0==cl1) sums++;
4361 if (AliTPCReconstructor::StreamLevel()>5) {
4362 TTreeSRedirector &cstream = *fDebugStreamer;
4367 "Tr0.="<<track0<< // seed0
4368 "Tr1.="<<track1<< // seed1
4369 "h0.="<<&helixes[i0]<<
4370 "h1.="<<&helixes[i1]<<
4372 "sum="<<sum<< //the sum of rows with cl in both
4373 "sums="<<sums<< //the sum of shared clusters
4374 "xm0="<<xm[i0]<< // the center of track
4375 "xm1="<<xm[i1]<< // the x center of track
4376 // General cut variables
4377 "dfi="<<dfi<< // distance in fi angle
4378 "dtheta="<<dtheta<< // distance int theta angle
4384 "dist0="<<dist[0]<< //distance x
4385 "dist1="<<dist[1]<< //distance y
4386 "dist2="<<dist[2]<< //distance z
4387 "mdist0="<<mdist[0]<< //distance x
4388 "mdist1="<<mdist[1]<< //distance y
4389 "mdist2="<<mdist[2]<< //distance z
4405 if (AliTPCReconstructor::StreamLevel()>1) {
4406 AliInfo("Time for curling tracks removal DEBUGGING MC");
4413 void AliTPCtrackerMI::FindSplitted(TObjArray * array, AliESDEvent */*esd*/, Int_t /*iter*/){
4415 // Find Splitted tracks and remove the one with worst quality
4416 // Corresponding debug streamer to tune selections - "Splitted2"
4418 // 0. Sort tracks according quility
4419 // 1. Propagate the tracks to the reference radius
4420 // 2. Double_t loop to select close tracks (only to speed up process)
4421 // 3. Calculate cluster overlap ratio - and remove the track if bigger than a threshold
4422 // 4. Delete temporary parameters
4424 const Double_t xref=GetXrow(63); // reference radius -IROC/OROC boundary
4426 const Double_t kCutP1=10; // delta Z cut 10 cm
4427 const Double_t kCutP2=0.15; // delta snp(fi) cut 0.15
4428 const Double_t kCutP3=0.15; // delta tgl(theta) cut 0.15
4429 const Double_t kCutAlpha=0.15; // delta alpha cut
4430 Int_t firstpoint = 0;
4431 Int_t lastpoint = 160;
4433 Int_t nentries = array->GetEntriesFast();
4434 AliExternalTrackParam *params = new AliExternalTrackParam[nentries];
4440 //0. Sort tracks according quality
4441 //1. Propagate the ext. param to reference radius
4442 Int_t nseed = array->GetEntriesFast();
4443 if (nseed<=0) return;
4444 Float_t * quality = new Float_t[nseed];
4445 Int_t * indexes = new Int_t[nseed];
4446 for (Int_t i=0; i<nseed; i++) {
4447 AliTPCseed *pt=(AliTPCseed*)array->UncheckedAt(i);
4452 pt->UpdatePoints(); //select first last max dens points
4453 Float_t * points = pt->GetPoints();
4454 if (points[3]<0.8) quality[i] =-1;
4455 quality[i] = (points[2]-points[0])+pt->GetNumberOfClusters();
4456 //prefer high momenta tracks if overlaps
4457 quality[i] *= TMath::Sqrt(TMath::Abs(pt->Pt())+0.5);
4459 AliTracker::PropagateTrackToBxByBz(&(params[i]),xref,pt->GetMass(),5.,kTRUE);
4460 AliTracker::PropagateTrackToBxByBz(&(params[i]),xref,pt->GetMass(),1.,kTRUE);
4462 TMath::Sort(nseed,quality,indexes);
4464 // 3. Loop over pair of tracks
4466 for (Int_t i0=0; i0<nseed; i0++) {
4467 Int_t index0=indexes[i0];
4468 if (!(array->UncheckedAt(index0))) continue;
4469 AliTPCseed *s1 = (AliTPCseed*)array->UncheckedAt(index0);
4470 if (!s1->IsActive()) continue;
4471 AliExternalTrackParam &par0=params[index0];
4472 for (Int_t i1=i0+1; i1<nseed; i1++) {
4473 Int_t index1=indexes[i1];
4474 if (!(array->UncheckedAt(index1))) continue;
4475 AliTPCseed *s2 = (AliTPCseed*)array->UncheckedAt(index1);
4476 if (!s2->IsActive()) continue;
4477 if (s2->GetKinkIndexes()[0]!=0)
4478 if (s2->GetKinkIndexes()[0] == -s1->GetKinkIndexes()[0]) continue;
4479 AliExternalTrackParam &par1=params[index1];
4480 if (TMath::Abs(par0.GetParameter()[3]-par1.GetParameter()[3])>kCutP3) continue;
4481 if (TMath::Abs(par0.GetParameter()[1]-par1.GetParameter()[1])>kCutP1) continue;
4482 if (TMath::Abs(par0.GetParameter()[2]-par1.GetParameter()[2])>kCutP2) continue;
4483 Double_t dAlpha= TMath::Abs(par0.GetAlpha()-par1.GetAlpha());
4484 if (dAlpha>TMath::Pi()) dAlpha-=TMath::Pi();
4485 if (TMath::Abs(dAlpha)>kCutAlpha) continue;
4490 Int_t firstShared=lastpoint, lastShared=firstpoint;
4491 Int_t firstRow=lastpoint, lastRow=firstpoint;
4493 for (Int_t i=firstpoint;i<lastpoint;i++){
4494 if (s1->GetClusterIndex2(i)>0) nall0++;
4495 if (s2->GetClusterIndex2(i)>0) nall1++;
4496 if (s1->GetClusterIndex2(i)>0 && s2->GetClusterIndex2(i)>0) {
4497 if (i<firstRow) firstRow=i;
4498 if (i>lastRow) lastRow=i;
4500 if ( (s1->GetClusterIndex2(i))==(s2->GetClusterIndex2(i)) && s1->GetClusterIndex2(i)>0) {
4501 if (i<firstShared) firstShared=i;
4502 if (i>lastShared) lastShared=i;
4506 Double_t ratio0 = Float_t(sumShared)/Float_t(TMath::Min(nall0+1,nall1+1));
4507 Double_t ratio1 = Float_t(sumShared)/Float_t(TMath::Max(nall0+1,nall1+1));
4509 if( AliTPCReconstructor::StreamLevel()>1){
4510 TTreeSRedirector &cstream = *fDebugStreamer;
4511 Int_t n0=s1->GetNumberOfClusters();
4512 Int_t n1=s2->GetNumberOfClusters();
4513 Int_t n0F=s1->GetNFoundable();
4514 Int_t n1F=s2->GetNFoundable();
4515 Int_t lab0=s1->GetLabel();
4516 Int_t lab1=s2->GetLabel();
4518 cstream<<"Splitted2"<<
4519 "iter="<<fIteration<<
4520 "lab0="<<lab0<< // MC label if exist
4521 "lab1="<<lab1<< // MC label if exist
4524 "ratio0="<<ratio0<< // shared ratio
4525 "ratio1="<<ratio1<< // shared ratio
4526 "p0.="<<&par0<< // track parameters
4528 "s0.="<<s1<< // full seed
4530 "n0="<<n0<< // number of clusters track 0
4531 "n1="<<n1<< // number of clusters track 1
4532 "nall0="<<nall0<< // number of clusters track 0
4533 "nall1="<<nall1<< // number of clusters track 1
4534 "n0F="<<n0F<< // number of findable
4535 "n1F="<<n1F<< // number of findable
4536 "shared="<<sumShared<< // number of shared clusters
4537 "firstS="<<firstShared<< // first and the last shared row
4538 "lastS="<<lastShared<<
4539 "firstRow="<<firstRow<< // first and the last row with cluster
4540 "lastRow="<<lastRow<< //
4544 // remove track with lower quality
4546 if (ratio0>AliTPCReconstructor::GetRecoParam()->GetCutSharedClusters(0) ||
4547 ratio1>AliTPCReconstructor::GetRecoParam()->GetCutSharedClusters(1)){
4551 fSeedsPool->MarkSlotFree( array->RemoveAt(index1) );
4556 // 4. Delete temporary array
4566 void AliTPCtrackerMI::FindCurling(const TObjArray * array, AliESDEvent */*esd*/, Int_t iter)
4569 // find Curling tracks
4570 // Use AliTPCReconstructor::StreamLevel()>1 if you want to tune parameters - cuts
4573 // Algorithm done in 2 phases - because of CPU consumption
4574 // it is n^2 algorithm - for lead-lead 20000x20000 combination are investigated
4575 // see detal in MC part what can be used to cut
4579 const Float_t kMaxC = 400; // maximal curvature to of the track
4580 const Float_t kMaxdTheta = 0.15; // maximal distance in theta
4581 const Float_t kMaxdPhi = 0.15; // maximal distance in phi
4582 const Float_t kPtRatio = 0.3; // ratio between pt
4583 const Float_t kMinDCAR = 2.; // distance to the primary vertex in r - see cpipe cut
4586 // Curling tracks cuts
4589 const Float_t kMaxDeltaRMax = 40; // distance in outer radius
4590 const Float_t kMaxDeltaRMin = 5.; // distance in lower radius - see cpipe cut
4591 const Float_t kMinAngle = 2.9; // angle between tracks
4592 const Float_t kMaxDist = 5; // biggest distance
4594 // The cuts can be tuned using the "MC information stored in Multi tree ==> see FindMultiMC
4597 TCut csign("csign","Tr0.fP[4]*Tr1.fP[4]<0"); //opposite sign
4598 TCut cmax("cmax","abs(Tr0.GetC())>1/400");
4599 TCut cda("cda","sqrt(dtheta^2+dfi^2)<0.15");
4600 TCut ccratio("ccratio","abs((Tr0.fP[4]+Tr1.fP[4])/(abs(Tr0.fP[4])+abs(Tr1.fP[4])))<0.3");
4601 TCut cpipe("cpipe", "min(abs(r0-rc0),abs(r1-rc1))>5");
4603 TCut cdrmax("cdrmax","abs(abs(rc0+r0)-abs(rc1+r1))<40")
4604 TCut cdrmin("cdrmin","abs(abs(rc0+r0)-abs(rc1+r1))<10")
4606 Multi->Draw("dfi","iter==0"+csign+cmax+cda+ccratio); ~94% of curling tracks fulfill
4607 Multi->Draw("min(abs(r0-rc0),abs(r1-rc1))","iter==0&&abs(lab1)==abs(lab0)"+csign+cmax+cda+ccratio+cpipe+cdrmin+cdrmax); //80%
4609 Curling2->Draw("dfi","iter==0&&abs(lab0)==abs(lab1)"+csign+cmax+cdtheta+cdfi+ccratio)
4615 Int_t nentries = array->GetEntriesFast();
4616 AliHelix *helixes = new AliHelix[nentries];
4617 for (Int_t i=0;i<nentries;i++){
4618 AliTPCseed* track = (AliTPCseed*)array->At(i);
4619 if (!track) continue;
4620 track->SetCircular(0);
4621 new (&helixes[i]) AliHelix(*track);
4627 Double_t phase[2][2]={{0,0},{0,0}},radius[2]={0,0};
4633 for (Int_t i0=0;i0<nentries;i0++){
4634 AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
4635 if (!track0) continue;
4636 if (TMath::Abs(track0->GetC())<1/kMaxC) continue;
4637 Float_t xc0 = helixes[i0].GetHelix(6);
4638 Float_t yc0 = helixes[i0].GetHelix(7);
4639 Float_t r0 = helixes[i0].GetHelix(8);
4640 Float_t rc0 = TMath::Sqrt(xc0*xc0+yc0*yc0);
4641 Float_t fi0 = TMath::ATan2(yc0,xc0);
4643 for (Int_t i1=i0+1;i1<nentries;i1++){
4644 AliTPCseed * track1 = (AliTPCseed*)array->At(i1);
4645 if (!track1) continue;
4646 if (TMath::Abs(track1->GetC())<1/kMaxC) continue;
4647 Float_t xc1 = helixes[i1].GetHelix(6);
4648 Float_t yc1 = helixes[i1].GetHelix(7);
4649 Float_t r1 = helixes[i1].GetHelix(8);
4650 Float_t rc1 = TMath::Sqrt(xc1*xc1+yc1*yc1);
4651 Float_t fi1 = TMath::ATan2(yc1,xc1);
4653 Float_t dfi = fi0-fi1;
4656 if (dfi>1.5*TMath::Pi()) dfi-=TMath::Pi(); // take care about edge effect
4657 if (dfi<-1.5*TMath::Pi()) dfi+=TMath::Pi(); //
4658 Float_t dtheta = TMath::Abs(track0->GetTgl()-track1->GetTgl())<TMath::Abs(track0->GetTgl()+track1->GetTgl())? track0->GetTgl()-track1->GetTgl():track0->GetTgl()+track1->GetTgl();
4662 if (track0->GetBConstrain()&&track1->GetBConstrain()) continue; // not constrained
4663 if (track1->GetSigned1Pt()*track0->GetSigned1Pt()>0) continue; // not the same sign
4664 if ( TMath::Abs(track1->GetTgl()+track0->GetTgl())>kMaxdTheta) continue; //distance in the Theta
4665 if ( TMath::Abs(dfi)>kMaxdPhi) continue; //distance in phi
4666 if ( TMath::Sqrt(dfi*dfi+dtheta*dtheta)>kMaxdPhi) continue; //common angular offset
4668 Float_t pt0 = track0->GetSignedPt();
4669 Float_t pt1 = track1->GetSignedPt();
4670 if ((TMath::Abs(pt0+pt1)/(TMath::Abs(pt0)+TMath::Abs(pt1)))>kPtRatio) continue;
4671 if ((iter==1) && TMath::Abs(TMath::Abs(rc0+r0)-TMath::Abs(rc1+r1))>kMaxDeltaRMax) continue;
4672 if ((iter!=1) &&TMath::Abs(TMath::Abs(rc0-r0)-TMath::Abs(rc1-r1))>kMaxDeltaRMin) continue;
4673 if (TMath::Min(TMath::Abs(rc0-r0),TMath::Abs(rc1-r1))<kMinDCAR) continue;
4676 // Now find closest approach
4680 Int_t npoints = helixes[i0].GetRPHIintersections(helixes[i1], phase, radius,10);
4681 if (npoints==0) continue;
4682 helixes[i0].GetClosestPhases(helixes[i1], phase);
4686 Double_t hangles[3];
4687 helixes[i0].Evaluate(phase[0][0],xyz0);
4688 helixes[i1].Evaluate(phase[0][1],xyz1);
4690 helixes[i0].GetAngle(phase[0][0],helixes[i1],phase[0][1],hangles);
4691 Double_t deltah[2],deltabest;
4692 if (TMath::Abs(hangles[2])<kMinAngle) continue;
4696 helixes[i0].ParabolicDCA(helixes[i1],phase[0][0],phase[0][1],radius[0],deltah[0],2);
4698 helixes[i0].ParabolicDCA(helixes[i1],phase[1][0],phase[1][1],radius[1],deltah[1],2);
4699 if (deltah[1]<deltah[0]) ibest=1;
4701 deltabest = TMath::Sqrt(deltah[ibest]);
4702 helixes[i0].Evaluate(phase[ibest][0],xyz0);
4703 helixes[i1].Evaluate(phase[ibest][1],xyz1);
4704 helixes[i0].GetAngle(phase[ibest][0],helixes[i1],phase[ibest][1],hangles);
4705 Double_t radiusbest = TMath::Sqrt(radius[ibest]);
4707 if (deltabest>kMaxDist) continue;
4708 // if (mindcar+mindcaz<40 && (TMath::Abs(hangles[2])<kMinAngle ||deltabest>3)) continue;
4709 Bool_t sign =kFALSE;
4710 if (hangles[2]>kMinAngle) sign =kTRUE;
4713 // circular[i0] = kTRUE;
4714 // circular[i1] = kTRUE;
4715 if (track0->OneOverPt()<track1->OneOverPt()){
4716 track0->SetCircular(track0->GetCircular()+1);
4717 track1->SetCircular(track1->GetCircular()+2);
4720 track1->SetCircular(track1->GetCircular()+1);
4721 track0->SetCircular(track0->GetCircular()+2);
4724 if (AliTPCReconstructor::StreamLevel()>2){
4726 //debug stream to tune "fine" cuts
4727 Int_t lab0=track0->GetLabel();
4728 Int_t lab1=track1->GetLabel();
4729 TTreeSRedirector &cstream = *fDebugStreamer;
4730 cstream<<"Curling2"<<
4746 "npoints="<<npoints<<
4747 "hangles0="<<hangles[0]<<
4748 "hangles1="<<hangles[1]<<
4749 "hangles2="<<hangles[2]<<
4752 "radius="<<radiusbest<<
4753 "deltabest="<<deltabest<<
4754 "phase0="<<phase[ibest][0]<<
4755 "phase1="<<phase[ibest][1]<<
4763 if (AliTPCReconstructor::StreamLevel()>1) {
4764 AliInfo("Time for curling tracks removal");
4770 void AliTPCtrackerMI::FindKinks(TObjArray * array, AliESDEvent *esd)
4776 // RS something is wrong in this routine: not all seeds are assigned to daughters and mothers array, but they all are queried
4778 if (fPools && !fPools->GetPool(AliPoolsSet::kPoolTPCKink)) fPools->SetPool(fKinksPool,AliPoolsSet::kPoolTPCKink);
4779 fKinksPool->Clear();
4780 // TObjArray *v0s= new TObjArray(10000);
4781 Int_t nentries = array->GetEntriesFast();
4782 AliHelix *helixes = new AliHelix[nentries];
4783 Int_t sign[nentries];
4784 Int_t nclusters[nentries];
4785 Float_t alpha[nentries];
4786 Int_t usage[nentries];
4787 Float_t zm[nentries];
4788 Float_t z0[nentries];
4789 Float_t fim[nentries];
4790 Float_t shared[nentries];
4791 Bool_t circular[nentries];
4792 Float_t dca[nentries];
4793 //const AliESDVertex * primvertex = esd->GetVertex();
4795 // nentries = array->GetEntriesFast();
4798 for (Int_t i=0;i<nentries;i++){
4801 AliTPCseed* track = (AliTPCseed*)array->At(i);
4802 if (!track) continue;
4803 track->SetCircular(0);
4805 track->UpdatePoints();
4806 if (( track->GetPoints()[2]- track->GetPoints()[0])>5 && track->GetPoints()[3]>0.8){
4808 nclusters[i]=track->GetNumberOfClusters();
4809 alpha[i] = track->GetAlpha();
4810 new (&helixes[i]) AliHelix(*track);
4812 helixes[i].Evaluate(0,xyz);
4813 sign[i] = (track->GetC()>0) ? -1:1;
4816 if (track->GetProlongation(x,y,z)){
4818 fim[i] = alpha[i]+TMath::ATan2(y,x);
4821 zm[i] = track->GetZ();
4825 circular[i]= kFALSE;
4826 if (track->GetProlongation(0,y,z)) z0[i] = z;
4827 dca[i] = track->GetD(0,0);
4833 Int_t ncandidates =0;
4836 Double_t phase[2][2]={{0,0},{0,0}},radius[2]={0,0};
4839 // Find circling track
4841 for (Int_t i0=0;i0<nentries;i0++){
4842 AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
4843 if (!track0) continue;
4844 if (track0->GetNumberOfClusters()<40) continue;
4845 if (TMath::Abs(1./track0->GetC())>200) continue;
4846 for (Int_t i1=i0+1;i1<nentries;i1++){
4847 AliTPCseed * track1 = (AliTPCseed*)array->At(i1);
4848 if (!track1) continue;
4849 if (track1->GetNumberOfClusters()<40) continue;
4850 if ( TMath::Abs(track1->GetTgl()+track0->GetTgl())>0.1) continue;
4851 if (track0->GetBConstrain()&&track1->GetBConstrain()) continue;
4852 if (TMath::Abs(1./track1->GetC())>200) continue;
4853 if (track1->GetSigned1Pt()*track0->GetSigned1Pt()>0) continue;
4854 if (track1->GetTgl()*track0->GetTgl()>0) continue;
4855 if (TMath::Max(TMath::Abs(1./track0->GetC()),TMath::Abs(1./track1->GetC()))>190) continue;
4856 if (track0->GetBConstrain()&&track1->OneOverPt()<track0->OneOverPt()) continue; //returning - lower momenta
4857 if (track1->GetBConstrain()&&track0->OneOverPt()<track1->OneOverPt()) continue; //returning - lower momenta
4859 Float_t mindcar = TMath::Min(TMath::Abs(dca[i0]),TMath::Abs(dca[i1]));
4860 if (mindcar<5) continue;
4861 Float_t mindcaz = TMath::Min(TMath::Abs(z0[i0]-GetZ()),TMath::Abs(z0[i1]-GetZ()));
4862 if (mindcaz<5) continue;
4863 if (mindcar+mindcaz<20) continue;
4866 Float_t xc0 = helixes[i0].GetHelix(6);
4867 Float_t yc0 = helixes[i0].GetHelix(7);
4868 Float_t r0 = helixes[i0].GetHelix(8);
4869 Float_t xc1 = helixes[i1].GetHelix(6);
4870 Float_t yc1 = helixes[i1].GetHelix(7);
4871 Float_t r1 = helixes[i1].GetHelix(8);
4873 Float_t rmean = (r0+r1)*0.5;
4874 Float_t delta =TMath::Sqrt((xc1-xc0)*(xc1-xc0)+(yc1-yc0)*(yc1-yc0));
4875 //if (delta>30) continue;
4876 if (delta>rmean*0.25) continue;
4877 if (TMath::Abs(r0-r1)/rmean>0.3) continue;
4879 Int_t npoints = helixes[i0].GetRPHIintersections(helixes[i1], phase, radius,10);
4880 if (npoints==0) continue;
4881 helixes[i0].GetClosestPhases(helixes[i1], phase);
4885 Double_t hangles[3];
4886 helixes[i0].Evaluate(phase[0][0],xyz0);
4887 helixes[i1].Evaluate(phase[0][1],xyz1);
4889 helixes[i0].GetAngle(phase[0][0],helixes[i1],phase[0][1],hangles);
4890 Double_t deltah[2],deltabest;
4891 if (hangles[2]<2.8) continue;
4894 helixes[i0].ParabolicDCA(helixes[i1],phase[0][0],phase[0][1],radius[0],deltah[0],2);
4896 helixes[i0].ParabolicDCA(helixes[i1],phase[1][0],phase[1][1],radius[1],deltah[1],2);
4897 if (deltah[1]<deltah[0]) ibest=1;
4899 deltabest = TMath::Sqrt(deltah[ibest]);
4900 helixes[i0].Evaluate(phase[ibest][0],xyz0);
4901 helixes[i1].Evaluate(phase[ibest][1],xyz1);
4902 helixes[i0].GetAngle(phase[ibest][0],helixes[i1],phase[ibest][1],hangles);
4903 Double_t radiusbest = TMath::Sqrt(radius[ibest]);
4905 if (deltabest>6) continue;
4906 if (mindcar+mindcaz<40 && (hangles[2]<3.12||deltabest>3)) continue;
4907 Bool_t lsign =kFALSE;
4908 if (hangles[2]>3.06) lsign =kTRUE;
4911 circular[i0] = kTRUE;
4912 circular[i1] = kTRUE;
4913 if (track0->OneOverPt()<track1->OneOverPt()){
4914 track0->SetCircular(track0->GetCircular()+1);
4915 track1->SetCircular(track1->GetCircular()+2);
4918 track1->SetCircular(track1->GetCircular()+1);
4919 track0->SetCircular(track0->GetCircular()+2);
4922 if (lsign&&AliTPCReconstructor::StreamLevel()>1){
4924 Int_t lab0=track0->GetLabel();
4925 Int_t lab1=track1->GetLabel();
4926 TTreeSRedirector &cstream = *fDebugStreamer;
4927 cstream<<"Curling"<<
4934 "mindcar="<<mindcar<<
4935 "mindcaz="<<mindcaz<<
4938 "npoints="<<npoints<<
4939 "hangles0="<<hangles[0]<<
4940 "hangles2="<<hangles[2]<<
4945 "radius="<<radiusbest<<
4946 "deltabest="<<deltabest<<
4947 "phase0="<<phase[ibest][0]<<
4948 "phase1="<<phase[ibest][1]<<
4957 AliKink* kink = new( fKinksPool->NextFreeSlot() ) AliKink(); fKinksPool->RegisterClone(kink);
4959 for (Int_t i =0;i<nentries;i++){
4960 if (sign[i]==0) continue;
4961 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
4968 Double_t cradius0 = 40*40;
4969 Double_t cradius1 = 270*270;
4972 Double_t cdist3=0.55;
4973 for (Int_t j =i+1;j<nentries;j++){
4975 if (sign[j]*sign[i]<1) continue;
4976 if ( (nclusters[i]+nclusters[j])>200) continue;
4977 if ( (nclusters[i]+nclusters[j])<80) continue;
4978 if ( TMath::Abs(zm[i]-zm[j])>60.) continue;
4979 if ( TMath::Abs(fim[i]-fim[j])>0.6 && TMath::Abs(fim[i]-fim[j])<5.7 ) continue;
4980 //AliTPCseed * track1 = (AliTPCseed*)array->At(j); Double_t phase[2][2],radius[2];
4981 Int_t npoints = helixes[i].GetRPHIintersections(helixes[j], phase, radius,20);
4982 if (npoints<1) continue;
4985 if (radius[0]<cradius0||radius[0]>cradius1) continue;
4988 if ( (radius[0]<cradius0||radius[0]>cradius1) && (radius[1]<cradius0||radius[1]>cradius1) ) continue;
4991 Double_t delta1=10000,delta2=10000;
4992 // cuts on the intersection radius
4993 helixes[i].LinearDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta1);
4994 if (radius[0]<20&&delta1<1) continue; //intersection at vertex
4995 if (radius[0]<10&&delta1<3) continue; //intersection at vertex
4997 helixes[i].LinearDCA(helixes[j],phase[1][0],phase[1][1],radius[1],delta2);
4998 if (radius[1]<20&&delta2<1) continue; //intersection at vertex
4999 if (radius[1]<10&&delta2<3) continue; //intersection at vertex
5002 Double_t distance1 = TMath::Min(delta1,delta2);
5003 if (distance1>cdist1) continue; // cut on DCA linear approximation
5005 npoints = helixes[i].GetRPHIintersections(helixes[j], phase, radius,20);
5006 helixes[i].ParabolicDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta1);
5007 if (radius[0]<20&&delta1<1) continue; //intersection at vertex
5008 if (radius[0]<10&&delta1<3) continue; //intersection at vertex
5011 helixes[i].ParabolicDCA(helixes[j],phase[1][0],phase[1][1],radius[1],delta2);
5012 if (radius[1]<20&&delta2<1) continue; //intersection at vertex
5013 if (radius[1]<10&&delta2<3) continue; //intersection at vertex
5015 distance1 = TMath::Min(delta1,delta2);
5018 rkink = TMath::Sqrt(radius[0]);
5021 rkink = TMath::Sqrt(radius[1]);
5023 if (distance1>cdist2) continue;
5026 AliTPCseed * track1 = (AliTPCseed*)array->At(j);
5029 Int_t row0 = GetRowNumber(rkink);
5030 if (row0<10) continue;
5031 if (row0>150) continue;
5034 Float_t dens00=-1,dens01=-1;
5035 Float_t dens10=-1,dens11=-1;
5037 Int_t found,foundable,ishared;
5038 track0->GetClusterStatistic(0,row0-5, found, foundable,ishared,kFALSE);
5039 if (foundable>5) dens00 = Float_t(found)/Float_t(foundable);
5040 track0->GetClusterStatistic(row0+5,155, found, foundable,ishared,kFALSE);
5041 if (foundable>5) dens01 = Float_t(found)/Float_t(foundable);
5043 track1->GetClusterStatistic(0,row0-5, found, foundable,ishared,kFALSE);
5044 if (foundable>10) dens10 = Float_t(found)/Float_t(foundable);
5045 track1->GetClusterStatistic(row0+5,155, found, foundable,ishared,kFALSE);
5046 if (foundable>10) dens11 = Float_t(found)/Float_t(foundable);
5048 if (dens00<dens10 && dens01<dens11) continue;
5049 if (dens00>dens10 && dens01>dens11) continue;
5050 if (TMath::Max(dens00,dens10)<0.1) continue;
5051 if (TMath::Max(dens01,dens11)<0.3) continue;
5053 if (TMath::Min(dens00,dens10)>0.6) continue;
5054 if (TMath::Min(dens01,dens11)>0.6) continue;
5057 AliTPCseed * ktrack0, *ktrack1;
5066 if (TMath::Abs(ktrack0->GetC())>5) continue; // cut on the curvature for mother particle
5067 AliExternalTrackParam paramm(*ktrack0);
5068 AliExternalTrackParam paramd(*ktrack1);
5069 if (row0>60&&ktrack1->GetReference().GetX()>90.)new (¶md) AliExternalTrackParam(ktrack1->GetReference());
5072 kink->SetMother(paramm);
5073 kink->SetDaughter(paramd);
5076 Float_t x[3] = { kink->GetPosition()[0],kink->GetPosition()[1],kink->GetPosition()[2]};
5078 fkParam->Transform0to1(x,index);
5079 fkParam->Transform1to2(x,index);
5080 row0 = GetRowNumber(x[0]);
5082 if (kink->GetR()<100) continue;
5083 if (kink->GetR()>240) continue;
5084 if (kink->GetPosition()[2]/kink->GetR()>AliTPCReconstructor::GetCtgRange()) continue; //out of fiducial volume
5085 if (kink->GetDistance()>cdist3) continue;
5086 Float_t dird = kink->GetDaughterP()[0]*kink->GetPosition()[0]+kink->GetDaughterP()[1]*kink->GetPosition()[1]; // rough direction estimate
5087 if (dird<0) continue;
5089 Float_t dirm = kink->GetMotherP()[0]*kink->GetPosition()[0]+kink->GetMotherP()[1]*kink->GetPosition()[1]; // rough direction estimate
5090 if (dirm<0) continue;
5091 Float_t mpt = TMath::Sqrt(kink->GetMotherP()[0]*kink->GetMotherP()[0]+kink->GetMotherP()[1]*kink->GetMotherP()[1]);
5092 if (mpt<0.2) continue;
5095 //for high momenta momentum not defined well in first iteration
5096 Double_t qt = TMath::Sin(kink->GetAngle(2))*ktrack1->GetP();
5097 if (qt>0.35) continue;
5100 kink->SetLabel(CookLabel(ktrack0,0.4,0,row0),0);
5101 kink->SetLabel(CookLabel(ktrack1,0.4,row0,160),1);
5103 kink->SetTPCDensity(dens00,0,0);
5104 kink->SetTPCDensity(dens01,0,1);
5105 kink->SetTPCDensity(dens10,1,0);
5106 kink->SetTPCDensity(dens11,1,1);
5107 kink->SetIndex(i,0);
5108 kink->SetIndex(j,1);
5111 kink->SetTPCDensity(dens10,0,0);
5112 kink->SetTPCDensity(dens11,0,1);
5113 kink->SetTPCDensity(dens00,1,0);
5114 kink->SetTPCDensity(dens01,1,1);
5115 kink->SetIndex(j,0);
5116 kink->SetIndex(i,1);
5119 if (mpt<1||kink->GetAngle(2)>0.1){
5120 // angle and densities not defined yet
5121 if (kink->GetTPCDensityFactor()<0.8) continue;
5122 if ((2-kink->GetTPCDensityFactor())*kink->GetDistance() >0.25) continue;
5123 if (kink->GetAngle(2)*ktrack0->GetP()<0.003) continue; //too small angle
5124 if (kink->GetAngle(2)>0.2&&kink->GetTPCDensityFactor()<1.15) continue;
5125 if (kink->GetAngle(2)>0.2&&kink->GetTPCDensity(0,1)>0.05) continue;
5127 Float_t criticalangle = track0->GetSigmaSnp2()+track0->GetSigmaTgl2();
5128 criticalangle+= track1->GetSigmaSnp2()+track1->GetSigmaTgl2();
5129 criticalangle= 3*TMath::Sqrt(criticalangle);
5130 if (criticalangle>0.02) criticalangle=0.02;
5131 if (kink->GetAngle(2)<criticalangle) continue;
5134 Int_t drow = Int_t(2.+0.5/(0.05+kink->GetAngle(2))); // overlap region defined
5135 Float_t shapesum =0;
5137 for ( Int_t row = row0-drow; row<row0+drow;row++){
5138 if (row<0) continue;
5139 if (row>155) continue;
5140 if (ktrack0->GetClusterPointer(row)){
5141 AliTPCTrackerPoint *point =ktrack0->GetTrackPoint(row);
5142 shapesum+=point->GetSigmaY()+point->GetSigmaZ();
5145 if (ktrack1->GetClusterPointer(row)){
5146 AliTPCTrackerPoint *point =ktrack1->GetTrackPoint(row);
5147 shapesum+=point->GetSigmaY()+point->GetSigmaZ();
5152 kink->SetShapeFactor(-1.);
5155 kink->SetShapeFactor(shapesum/sum);
5157 // esd->AddKink(kink);
5159 // kink->SetMother(paramm);
5160 //kink->SetDaughter(paramd);
5162 Double_t chi2P2 = paramm.GetParameter()[2]-paramd.GetParameter()[2];
5164 chi2P2/=paramm.GetCovariance()[5]+paramd.GetCovariance()[5];
5165 Double_t chi2P3 = paramm.GetParameter()[3]-paramd.GetParameter()[3];
5167 chi2P3/=paramm.GetCovariance()[9]+paramd.GetCovariance()[9];
5169 if (AliTPCReconstructor::StreamLevel()>1) {
5170 (*fDebugStreamer)<<"kinkLpt"<<
5178 if ( chi2P2+chi2P3<AliTPCReconstructor::GetRecoParam()->GetKinkAngleCutChi2(0)){
5182 kink = new( fKinksPool->NextFreeSlot() ) AliKink(); fKinksPool->RegisterClone(kink);
5186 // eliminate last dummy kink
5187 fKinksPool->MarkSlotFree(kink);
5189 if (ncandidates!=fKinksPool->GetEntriesFast() ) {
5190 AliError(Form("Mismatch in number of kinks: %d %d",ncandidates,fKinksPool->GetEntriesFast()));
5192 // sort the kinks according quality - and refit them towards vertex
5194 Int_t nkinks = fKinksPool->GetEntriesFast();
5195 Float_t quality[nkinks];
5196 Int_t indexes[nkinks];
5197 AliTPCseed* mothers[nkinks];;
5198 AliTPCseed* daughters[nkinks];
5201 for (Int_t i=0;i<nkinks;i++){
5203 mothers[i] = daughters[i] = 0;
5204 AliKink *kinkl = (AliKink*)fKinksPool->At(i);
5206 // refit kinks towards vertex
5208 Int_t index0 = kinkl->GetIndex(0);
5209 Int_t index1 = kinkl->GetIndex(1);
5210 AliTPCseed * ktrack0 = (AliTPCseed*)array->At(index0);
5211 AliTPCseed * ktrack1 = (AliTPCseed*)array->At(index1);
5213 Int_t sumn=ktrack0->GetNumberOfClusters()+ktrack1->GetNumberOfClusters();
5215 // Refit Kink under if too small angle
5217 if (kinkl->GetAngle(2)<0.05){
5218 kinkl->SetTPCRow0(GetRowNumber(kinkl->GetR()));
5219 Int_t row0 = kinkl->GetTPCRow0();
5220 Int_t drow = Int_t(2.+0.5/(0.05+kinkl->GetAngle(2)));
5223 Int_t last = row0-drow;
5224 if (last<40) last=40;
5225 if (last<ktrack0->GetFirstPoint()+25) last = ktrack0->GetFirstPoint()+25;
5226 AliTPCseed* seed0 = ReSeed(ktrack0,last,kFALSE);
5229 Int_t first = row0+drow;
5230 if (first>130) first=130;
5231 if (first>ktrack1->GetLastPoint()-25) first = TMath::Max(ktrack1->GetLastPoint()-25,30);
5232 AliTPCseed* seed1 = ReSeed(ktrack1,first,kTRUE);
5234 if (seed0 && seed1) {
5235 kinkl->SetStatus(1,8);
5236 if (RefitKink(*seed0,*seed1,*kinkl)) kinkl->SetStatus(1,9);
5237 row0 = GetRowNumber(kinkl->GetR());
5238 sumn = seed0->GetNumberOfClusters()+seed1->GetNumberOfClusters();
5240 daughters[i] = seed1;
5243 fKinksPool->MarkSlotFree(kinkl);
5244 if (seed0) fSeedsPool->MarkSlotFree( seed0 );
5245 if (seed1) fSeedsPool->MarkSlotFree( seed1 );
5248 if (kinkl->GetDistance()>0.5 || kinkl->GetR()<110 || kinkl->GetR()>240) {
5249 fKinksPool->MarkSlotFree(kinkl);
5250 if (seed0) fSeedsPool->MarkSlotFree( seed0 );
5251 if (seed1) fSeedsPool->MarkSlotFree( seed1 );
5252 mothers[i] = daughters[i] = 0;
5257 if (kinkl) quality[i] = 160*((0.1+kinkl->GetDistance())*(2.-kinkl->GetTPCDensityFactor()))/(sumn+40.); //the longest -clossest will win
5259 TMath::Sort(nkinks,quality,indexes,kFALSE);
5261 //remove double find kinks
5263 for (Int_t ikink0=1;ikink0<nkinks;ikink0++){
5264 AliKink * kink0 = (AliKink*) fKinksPool->At(indexes[ikink0]);
5265 if (!kink0) continue;
5266 AliTPCseed *mother0 = mothers[indexes[ikink0]];
5267 AliTPCseed *daughter0 = daughters[indexes[ikink0]];
5268 if (!mother0 || !daughter0) continue;
5270 for (Int_t ikink1=0;ikink1<ikink0;ikink1++){
5271 AliKink * kink1 = (AliKink*) fKinksPool->At(indexes[ikink1]);
5272 if (!kink1) continue;
5273 AliTPCseed *mother1 = mothers[indexes[ikink1]];
5274 AliTPCseed *daughter1 = daughters[indexes[ikink1]];
5275 if (!mother1 || !daughter1) continue;
5276 // if not close kink continue
5277 if (TMath::Abs(kink1->GetPosition()[2]-kink0->GetPosition()[2])>10) continue;
5278 if (TMath::Abs(kink1->GetPosition()[1]-kink0->GetPosition()[1])>10) continue;
5279 if (TMath::Abs(kink1->GetPosition()[0]-kink0->GetPosition()[0])>10) continue;
5281 Int_t row0 = (kink0->GetTPCRow0()+kink1->GetTPCRow0())/2;
5290 for (Int_t i=0;i<row0;i++){
5291 if (mother0->GetClusterIndex(i)>0 && mother1->GetClusterIndex(i)>0){
5294 if (mother0->GetClusterIndex(i)==mother1->GetClusterIndex(i)){
5301 for (Int_t i=row0;i<158;i++){
5302 //if (daughter0.GetClusterIndex(i)>0 && daughter0.GetClusterIndex(i)>0){ // RS: Bug?
5303 if (daughter0->GetClusterIndex(i)>0 && daughter1->GetClusterIndex(i)>0){
5306 if (daughter0->GetClusterIndex(i)==daughter1->GetClusterIndex(i)){
5312 Float_t ratio = Float_t(same+1)/Float_t(both+1);
5313 Float_t ratiom = Float_t(samem+1)/Float_t(bothm+1);
5314 Float_t ratiod = Float_t(samed+1)/Float_t(bothd+1);
5315 if (ratio>0.3 && ratiom>0.5 &&ratiod>0.5) {
5316 Int_t sum0 = mother0->GetNumberOfClusters()+daughter0->GetNumberOfClusters();
5317 Int_t sum1 = mother1->GetNumberOfClusters()+daughter1->GetNumberOfClusters();
5319 shared[kink0->GetIndex(0)]= kTRUE;
5320 shared[kink0->GetIndex(1)]= kTRUE;
5321 fKinksPool->MarkSlotFree(kink0);
5326 shared[kink1->GetIndex(0)]= kTRUE;
5327 shared[kink1->GetIndex(1)]= kTRUE;
5328 fKinksPool->MarkSlotFree(kink1);
5329 fKinksPool->MarkSlotFree(kink1);
5337 for (Int_t i=0;i<nkinks;i++){
5338 AliKink * kinkl = (AliKink*) fKinksPool->At(indexes[i]);
5339 if (!kinkl) continue;
5340 kinkl->SetTPCRow0(GetRowNumber(kinkl->GetR()));
5341 Int_t index0 = kinkl->GetIndex(0);
5342 Int_t index1 = kinkl->GetIndex(1);
5343 if (circular[index0]||(circular[index1]&&kinkl->GetDistance()>0.2)) continue;
5344 kinkl->SetMultiple(usage[index0],0);
5345 kinkl->SetMultiple(usage[index1],1);
5346 if (kinkl->GetMultiple()[0]+kinkl->GetMultiple()[1]>2) continue;
5347 if (kinkl->GetMultiple()[0]+kinkl->GetMultiple()[1]>0 && quality[indexes[i]]>0.2) continue;
5348 if (kinkl->GetMultiple()[0]+kinkl->GetMultiple()[1]>0 && kinkl->GetDistance()>0.2) continue;
5349 if (circular[index0]||(circular[index1]&&kinkl->GetDistance()>0.1)) continue;
5351 AliTPCseed * ktrack0 = (AliTPCseed*)array->At(index0);
5352 AliTPCseed * ktrack1 = (AliTPCseed*)array->At(index1);
5353 if (!ktrack0 || !ktrack1) continue;
5354 Int_t index = esd->AddKink(kinkl);
5357 if ( ktrack0->GetKinkIndex(0)==0 && ktrack1->GetKinkIndex(0)==0) { //best kink
5358 if (mothers[indexes[i]] && daughters[indexes[i]]) {
5359 if (mothers[indexes[i]]->GetNumberOfClusters()>20 && daughters[indexes[i]]->GetNumberOfClusters()>20 &&
5360 (mothers[indexes[i]]->GetNumberOfClusters()+daughters[indexes[i]]->GetNumberOfClusters())>100) {
5361 *ktrack0 = *mothers[indexes[i]];
5362 *ktrack1 = *daughters[indexes[i]];
5367 ktrack0->SetKinkIndex(usage[index0],-(index+1));
5368 ktrack1->SetKinkIndex(usage[index1], (index+1));
5373 // Remove tracks corresponding to shared kink's
5375 for (Int_t i=0;i<nentries;i++){
5376 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5377 if (!track0) continue;
5378 if (track0->GetKinkIndex(0)!=0) continue;
5379 if (shared[i]) fSeedsPool->MarkSlotFree( array->RemoveAt(i) );
5383 RemoveUsed2(array,0.5,0.4,30);
5385 for (Int_t i=0;i<nentries;i++){
5386 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5387 if (!track0) continue;
5388 track0->CookdEdx(0.02,0.6);
5392 for (Int_t i=0;i<nentries;i++){
5393 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5394 if (!track0) continue;
5395 if (track0->Pt()<1.4) continue;
5396 //remove double high momenta tracks - overlapped with kink candidates
5399 for (Int_t icl=track0->GetFirstPoint();icl<track0->GetLastPoint(); icl++){
5400 if (track0->GetClusterPointer(icl)!=0){
5402 if (track0->GetClusterPointer(icl)->IsUsed(10)) ishared++;
5405 if (Float_t(ishared+1)/Float_t(all+1)>0.5) {
5406 fSeedsPool->MarkSlotFree( array->RemoveAt(i) );
5410 if (track0->GetKinkIndex(0)!=0) continue;
5411 if (track0->GetNumberOfClusters()<80) continue;
5413 AliTPCseed *pmother = new( fSeedsPool->NextFreeSlot() ) AliTPCseed(); fSeedsPool->RegisterClone(pmother);
5414 AliTPCseed *pdaughter = new( fSeedsPool->NextFreeSlot() ) AliTPCseed(); fSeedsPool->RegisterClone(pdaughter);
5415 AliKink *pkink = new( fKinksPool->NextFreeSlot() ) AliKink(); fKinksPool->RegisterClone(pkink);
5417 AliTPCseed & mother = *pmother;
5418 AliTPCseed & daughter = *pdaughter;
5419 AliKink & kinkl = *pkink;
5420 if (CheckKinkPoint(track0,mother,daughter, kinkl)){
5421 if (mother.GetNumberOfClusters()<30||daughter.GetNumberOfClusters()<20) {
5422 fSeedsPool->MarkSlotFree(pmother);
5423 fSeedsPool->MarkSlotFree(pdaughter);
5424 fKinksPool->MarkSlotFree(pkink);
5425 continue; //too short tracks
5427 if (mother.Pt()<1.4) {
5428 fSeedsPool->MarkSlotFree(pmother);
5429 fSeedsPool->MarkSlotFree(pdaughter);
5430 fKinksPool->MarkSlotFree(pkink);
5433 Int_t row0= kinkl.GetTPCRow0();
5434 if (kinkl.GetDistance()>0.5 || kinkl.GetR()<110. || kinkl.GetR()>240.) {
5435 fSeedsPool->MarkSlotFree(pmother);
5436 fSeedsPool->MarkSlotFree(pdaughter);
5437 fKinksPool->MarkSlotFree(pkink);
5441 Int_t index = esd->AddKink(&kinkl);
5442 mother.SetKinkIndex(0,-(index+1));
5443 daughter.SetKinkIndex(0,index+1);
5444 if (mother.GetNumberOfClusters()>50) {
5445 fSeedsPool->MarkSlotFree( array->RemoveAt(i) );
5446 AliTPCseed* mtc = new( fSeedsPool->NextFreeSlot() ) AliTPCseed(mother);
5447 fSeedsPool->RegisterClone(mtc);
5448 array->AddAt(mtc,i);
5451 AliTPCseed* mtc = new( fSeedsPool->NextFreeSlot() ) AliTPCseed(mother);
5452 fSeedsPool->RegisterClone(mtc);
5453 array->AddLast(mtc);
5455 AliTPCseed* dtc = new( fSeedsPool->NextFreeSlot() ) AliTPCseed(daughter);
5456 fSeedsPool->RegisterClone(dtc);
5457 array->AddLast(dtc);
5458 for (Int_t icl=0;icl<row0;icl++) {
5459 if (mother.GetClusterPointer(icl)) mother.GetClusterPointer(icl)->Use(20);
5462 for (Int_t icl=row0;icl<158;icl++) {
5463 if (daughter.GetClusterPointer(icl)) daughter.GetClusterPointer(icl)->Use(20);
5467 fSeedsPool->MarkSlotFree(pmother);
5468 fSeedsPool->MarkSlotFree(pdaughter);
5469 fKinksPool->MarkSlotFree(pkink);
5472 for (int i=nkinks;i--;) {
5473 if (daughters[i]) fSeedsPool->MarkSlotFree(daughters[i]);
5474 if (mothers[i]) fSeedsPool->MarkSlotFree(mothers[i]);
5476 fKinksPool->Clear();
5479 AliInfo(Form("Ncandidates=\t%d\t%d\t%d\t%d\n",esd->GetNumberOfKinks(),ncandidates,ntracks,nall));
5484 Int_t AliTPCtrackerMI::RefitKink(AliTPCseed &mother, AliTPCseed &daughter, const AliESDkink &knk)
5487 // refit kink towards to the vertex
5490 AliKink &kink=(AliKink &)knk;
5492 Int_t row0 = GetRowNumber(kink.GetR());
5493 FollowProlongation(mother,0);
5494 mother.Reset(kFALSE);
5496 FollowProlongation(daughter,row0);
5497 daughter.Reset(kFALSE);
5498 FollowBackProlongation(daughter,158);
5499 daughter.Reset(kFALSE);
5500 Int_t first = TMath::Max(row0-20,30);
5501 Int_t last = TMath::Min(row0+20,140);
5503 const Int_t kNdiv =5;
5504 AliTPCseed param0[kNdiv]; // parameters along the track
5505 AliTPCseed param1[kNdiv]; // parameters along the track
5506 AliKink kinks[kNdiv]; // corresponding kink parameters
5509 for (Int_t irow=0; irow<kNdiv;irow++){
5510 rows[irow] = first +((last-first)*irow)/(kNdiv-1);
5512 // store parameters along the track
5514 for (Int_t irow=0;irow<kNdiv;irow++){
5515 FollowBackProlongation(mother, rows[irow]);
5516 FollowProlongation(daughter,rows[kNdiv-1-irow]);
5517 param0[irow] = mother;
5518 param1[kNdiv-1-irow] = daughter;
5522 for (Int_t irow=0; irow<kNdiv-1;irow++){
5523 if (param0[irow].GetNumberOfClusters()<kNdiv||param1[irow].GetNumberOfClusters()<kNdiv) continue;
5524 kinks[irow].SetMother(param0[irow]);
5525 kinks[irow].SetDaughter(param1[irow]);
5526 kinks[irow].Update();
5529 // choose kink with best "quality"
5531 Double_t mindist = 10000;
5532 for (Int_t irow=0;irow<kNdiv;irow++){
5533 if (param0[irow].GetNumberOfClusters()<20||param1[irow].GetNumberOfClusters()<20) continue;
5534 if (TMath::Abs(kinks[irow].GetR())>240.) continue;
5535 if (TMath::Abs(kinks[irow].GetR())<100.) continue;
5537 Float_t normdist = TMath::Abs(param0[irow].GetX()-kinks[irow].GetR())*(0.1+kink.GetDistance());
5538 normdist/= (param0[irow].GetNumberOfClusters()+param1[irow].GetNumberOfClusters()+40.);
5539 if (normdist < mindist){
5545 if (index==-1) return 0;
5548 param0[index].Reset(kTRUE);
5549 FollowProlongation(param0[index],0);
5551 mother = param0[index];
5552 daughter = param1[index]; // daughter in vertex
5554 kink.SetMother(mother);
5555 kink.SetDaughter(daughter);
5557 kink.SetTPCRow0(GetRowNumber(kink.GetR()));
5558 kink.SetTPCncls(param0[index].GetNumberOfClusters(),0);
5559 kink.SetTPCncls(param1[index].GetNumberOfClusters(),1);
5560 kink.SetLabel(CookLabel(&mother,0.4, 0,kink.GetTPCRow0()),0);
5561 kink.SetLabel(CookLabel(&daughter,0.4, kink.GetTPCRow0(),160),1);
5562 mother.SetLabel(kink.GetLabel(0));
5563 daughter.SetLabel(kink.GetLabel(1));
5569 void AliTPCtrackerMI::UpdateKinkQualityM(AliTPCseed * seed){
5571 // update Kink quality information for mother after back propagation
5573 if (seed->GetKinkIndex(0)>=0) return;
5574 for (Int_t ikink=0;ikink<3;ikink++){
5575 Int_t index = seed->GetKinkIndex(ikink);
5576 if (index>=0) break;
5577 index = TMath::Abs(index)-1;
5578 AliESDkink * kink = fEvent->GetKink(index);
5579 kink->SetTPCDensity(-1,0,0);
5580 kink->SetTPCDensity(1,0,1);
5582 Int_t row0 = kink->GetTPCRow0() - 2 - Int_t( 0.5/ (0.05+kink->GetAngle(2)));
5583 if (row0<15) row0=15;
5585 Int_t row1 = kink->GetTPCRow0() + 2 + Int_t( 0.5/ (0.05+kink->GetAngle(2)));
5586 if (row1>145) row1=145;
5588 Int_t found,foundable,shared;
5589 seed->GetClusterStatistic(0,row0, found, foundable,shared,kFALSE);
5590 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),0,0);
5591 seed->GetClusterStatistic(row1,155, found, foundable,shared,kFALSE);
5592 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),0,1);
5597 void AliTPCtrackerMI::UpdateKinkQualityD(AliTPCseed * seed){
5599 // update Kink quality information for daughter after refit
5601 if (seed->GetKinkIndex(0)<=0) return;
5602 for (Int_t ikink=0;ikink<3;ikink++){
5603 Int_t index = seed->GetKinkIndex(ikink);
5604 if (index<=0) break;
5605 index = TMath::Abs(index)-1;
5606 AliESDkink * kink = fEvent->GetKink(index);
5607 kink->SetTPCDensity(-1,1,0);
5608 kink->SetTPCDensity(-1,1,1);
5610 Int_t row0 = kink->GetTPCRow0() -2 - Int_t( 0.5/ (0.05+kink->GetAngle(2)));
5611 if (row0<15) row0=15;
5613 Int_t row1 = kink->GetTPCRow0() +2 + Int_t( 0.5/ (0.05+kink->GetAngle(2)));
5614 if (row1>145) row1=145;
5616 Int_t found,foundable,shared;
5617 seed->GetClusterStatistic(0,row0, found, foundable,shared,kFALSE);
5618 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),1,0);
5619 seed->GetClusterStatistic(row1,155, found, foundable,shared,kFALSE);
5620 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),1,1);
5626 Int_t AliTPCtrackerMI::CheckKinkPoint(AliTPCseed*seed,AliTPCseed &mother, AliTPCseed &daughter, const AliESDkink &knk)
5629 // check kink point for given track
5630 // if return value=0 kink point not found
5631 // otherwise seed0 correspond to mother particle
5632 // seed1 correspond to daughter particle
5633 // kink parameter of kink point
5634 AliKink &kink=(AliKink &)knk;
5636 Int_t middlerow = (seed->GetFirstPoint()+seed->GetLastPoint())/2;
5637 Int_t first = seed->GetFirstPoint();
5638 Int_t last = seed->GetLastPoint();
5639 if (last-first<20) return 0; // shortest length - 2*30 = 60 pad-rows
5642 AliTPCseed *seed1 = ReSeed(seed,middlerow+20, kTRUE); //middle of chamber
5643 if (!seed1) return 0;
5644 FollowProlongation(*seed1,seed->GetLastPoint()-20);
5645 seed1->Reset(kTRUE);
5646 FollowProlongation(*seed1,158);
5647 seed1->Reset(kTRUE);
5648 last = seed1->GetLastPoint();
5650 AliTPCseed *seed0 = new( fSeedsPool->NextFreeSlot() ) AliTPCseed(*seed);
5651 fSeedsPool->RegisterClone(seed0);
5652 seed0->Reset(kFALSE);
5655 AliTPCseed param0[20]; // parameters along the track
5656 AliTPCseed param1[20]; // parameters along the track
5657 AliKink kinks[20]; // corresponding kink parameters
5659 for (Int_t irow=0; irow<20;irow++){
5660 rows[irow] = first +((last-first)*irow)/19;
5662 // store parameters along the track
5664 for (Int_t irow=0;irow<20;irow++){
5665 FollowBackProlongation(*seed0, rows[irow]);
5666 FollowProlongation(*seed1,rows[19-irow]);
5667 param0[irow] = *seed0;
5668 param1[19-irow] = *seed1;
5672 for (Int_t irow=0; irow<19;irow++){
5673 kinks[irow].SetMother(param0[irow]);
5674 kinks[irow].SetDaughter(param1[irow]);
5675 kinks[irow].Update();
5678 // choose kink with biggest change of angle
5680 Double_t maxchange= 0;
5681 for (Int_t irow=1;irow<19;irow++){
5682 if (TMath::Abs(kinks[irow].GetR())>240.) continue;
5683 if (TMath::Abs(kinks[irow].GetR())<110.) continue;
5684 Float_t quality = TMath::Abs(kinks[irow].GetAngle(2))/(3.+TMath::Abs(kinks[irow].GetR()-param0[irow].GetX()));
5685 if ( quality > maxchange){
5686 maxchange = quality;
5691 fSeedsPool->MarkSlotFree( seed0 );
5692 fSeedsPool->MarkSlotFree( seed1 );
5693 if (index<0) return 0;
5695 Int_t row0 = GetRowNumber(kinks[index].GetR()); //row 0 estimate
5696 seed0 = new( fSeedsPool->NextFreeSlot() ) AliTPCseed(param0[index]);
5697 fSeedsPool->RegisterClone(seed0);
5698 seed1 = new( fSeedsPool->NextFreeSlot() ) AliTPCseed(param1[index]);
5699 fSeedsPool->RegisterClone(seed1);
5700 seed0->Reset(kFALSE);
5701 seed1->Reset(kFALSE);
5702 seed0->ResetCovariance(10.);
5703 seed1->ResetCovariance(10.);
5704 FollowProlongation(*seed0,0);
5705 FollowBackProlongation(*seed1,158);
5706 mother = *seed0; // backup mother at position 0
5707 seed0->Reset(kFALSE);
5708 seed1->Reset(kFALSE);
5709 seed0->ResetCovariance(10.);
5710 seed1->ResetCovariance(10.);
5712 first = TMath::Max(row0-20,0);
5713 last = TMath::Min(row0+20,158);
5715 for (Int_t irow=0; irow<20;irow++){
5716 rows[irow] = first +((last-first)*irow)/19;
5718 // store parameters along the track
5720 for (Int_t irow=0;irow<20;irow++){
5721 FollowBackProlongation(*seed0, rows[irow]);
5722 FollowProlongation(*seed1,rows[19-irow]);
5723 param0[irow] = *seed0;
5724 param1[19-irow] = *seed1;
5728 for (Int_t irow=0; irow<19;irow++){
5729 kinks[irow].SetMother(param0[irow]);
5730 kinks[irow].SetDaughter(param1[irow]);
5731 // param0[irow].Dump();
5732 //param1[irow].Dump();
5733 kinks[irow].Update();
5736 // choose kink with biggest change of angle
5739 for (Int_t irow=0;irow<20;irow++){
5740 if (TMath::Abs(kinks[irow].GetR())>250.) continue;
5741 if (TMath::Abs(kinks[irow].GetR())<90.) continue;
5742 Float_t quality = TMath::Abs(kinks[irow].GetAngle(2))/(3.+TMath::Abs(kinks[irow].GetR()-param0[irow].GetX()));
5743 if ( quality > maxchange){
5744 maxchange = quality;
5751 if (index==-1 || param0[index].GetNumberOfClusters()+param1[index].GetNumberOfClusters()<100){
5752 fSeedsPool->MarkSlotFree( seed0 );
5753 fSeedsPool->MarkSlotFree( seed1 );
5757 // Float_t anglesigma = TMath::Sqrt(param0[index].fC22+param0[index].fC33+param1[index].fC22+param1[index].fC33);
5759 kink.SetMother(param0[index]);
5760 kink.SetDaughter(param1[index]);
5763 Double_t chi2P2 = param0[index].GetParameter()[2]-param1[index].GetParameter()[2];
5765 chi2P2/=param0[index].GetCovariance()[5]+param1[index].GetCovariance()[5];
5766 Double_t chi2P3 = param0[index].GetParameter()[3]-param1[index].GetParameter()[3];
5768 chi2P3/=param0[index].GetCovariance()[9]+param1[index].GetCovariance()[9];
5770 if (AliTPCReconstructor::StreamLevel()>1) {
5771 (*fDebugStreamer)<<"kinkHpt"<<
5774 "p0.="<<¶m0[index]<<
5775 "p1.="<<¶m1[index]<<
5779 if ( chi2P2+chi2P3<AliTPCReconstructor::GetRecoParam()->GetKinkAngleCutChi2(0)){
5780 fSeedsPool->MarkSlotFree( seed0 );
5781 fSeedsPool->MarkSlotFree( seed1 );
5786 row0 = GetRowNumber(kink.GetR());
5787 kink.SetTPCRow0(row0);
5788 kink.SetLabel(CookLabel(seed0,0.5,0,row0),0);
5789 kink.SetLabel(CookLabel(seed1,0.5,row0,158),1);
5790 kink.SetIndex(-10,0);
5791 kink.SetIndex(int(param0[index].GetNumberOfClusters()+param1[index].GetNumberOfClusters()),1);
5792 kink.SetTPCncls(param0[index].GetNumberOfClusters(),0);
5793 kink.SetTPCncls(param1[index].GetNumberOfClusters(),1);
5796 // new (&mother) AliTPCseed(param0[index]);
5797 daughter = param1[index];
5798 daughter.SetLabel(kink.GetLabel(1));
5799 param0[index].Reset(kTRUE);
5800 FollowProlongation(param0[index],0);
5801 mother = param0[index];
5802 mother.SetLabel(kink.GetLabel(0));
5803 if ( chi2P2+chi2P3<AliTPCReconstructor::GetRecoParam()->GetKinkAngleCutChi2(1)){
5806 fSeedsPool->MarkSlotFree( seed0 );
5807 fSeedsPool->MarkSlotFree( seed1 );
5815 AliTPCseed* AliTPCtrackerMI::ReSeed(AliTPCseed *t)
5818 // reseed - refit - track
5821 // Int_t last = fSectors->GetNRows()-1;
5823 if (fSectors == fOuterSec){
5824 first = TMath::Max(first, t->GetFirstPoint()-fInnerSec->GetNRows());
5828 first = t->GetFirstPoint();
5830 AliTPCseed * seed = MakeSeed(t,0.1,0.5,0.9);
5831 FollowBackProlongation(*t,fSectors->GetNRows()-1);
5833 FollowProlongation(*t,first);
5843 //_____________________________________________________________________________
5844 Int_t AliTPCtrackerMI::ReadSeeds(const TFile *inp) {
5845 //-----------------------------------------------------------------
5846 // This function reades track seeds.
5847 //-----------------------------------------------------------------
5848 TDirectory *savedir=gDirectory;
5850 TFile *in=(TFile*)inp;
5851 if (!in->IsOpen()) {
5852 cerr<<"AliTPCtrackerMI::ReadSeeds(): input file is not open !\n";
5857 TTree *seedTree=(TTree*)in->Get("Seeds");
5859 cerr<<"AliTPCtrackerMI::ReadSeeds(): ";
5860 cerr<<"can't get a tree with track seeds !\n";
5863 AliTPCtrack *seed=new AliTPCtrack;
5864 seedTree->SetBranchAddress("tracks",&seed);
5866 if (fSeeds==0) fSeeds=new TObjArray(15000);
5868 Int_t n=(Int_t)seedTree->GetEntries();
5869 for (Int_t i=0; i<n; i++) {
5870 seedTree->GetEvent(i);
5871 AliTPCseed* sdc = new( fSeedsPool->NextFreeSlot() ) AliTPCseed(*seed/*,seed->GetAlpha()*/);
5872 fSeedsPool->RegisterClone(sdc);
5873 fSeeds->AddLast(sdc);
5876 delete seed; // RS: this seed is not from the pool, delete it !!!
5882 Int_t AliTPCtrackerMI::Clusters2Tracks (AliESDEvent *const esd)
5885 // clusters to tracks
5886 if (fPools && !fPools->GetPoolTPCSeed()) fPools->SetPool(fSeedsPool,AliPoolsSet::kPoolTPCSeed);
5887 if (fSeeds) DeleteSeeds();
5888 else fSeedsPool->Clear("C");
5891 AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ;
5892 transform->SetCurrentTimeStamp( esd->GetTimeStamp());
5893 transform->SetCurrentRun(esd->GetRunNumber());
5896 if (!fSeeds) return 1;
5898 if (AliTPCReconstructor::StreamLevel()>3) DumpClusters(0,fSeeds);
5904 //_____________________________________________________________________________
5905 Int_t AliTPCtrackerMI::Clusters2Tracks() {
5906 //-----------------------------------------------------------------
5907 // This is a track finder.
5908 //-----------------------------------------------------------------
5909 TDirectory *savedir=gDirectory;
5913 fSeeds = Tracking();
5916 Info("Clusters2Tracks","Time for tracking: \t");timer.Print();timer.Start();
5918 //activate again some tracks
5919 for (Int_t i=0; i<fSeeds->GetEntriesFast(); i++) {
5920 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
5922 Int_t nc=t.GetNumberOfClusters();
5924 fSeedsPool->MarkSlotFree( fSeeds->RemoveAt(i) );
5928 if (pt->GetRemoval()==10) {
5929 if (pt->GetDensityFirst(20)>0.8 || pt->GetDensityFirst(30)>0.8 || pt->GetDensityFirst(40)>0.7)
5930 pt->Desactivate(10); // make track again active // MvL: should be 0 ?
5932 pt->Desactivate(20);
5933 fSeedsPool->MarkSlotFree( fSeeds->RemoveAt(i) );
5938 RemoveUsed2(fSeeds,0.85,0.85,0);
5939 if (AliTPCReconstructor::GetRecoParam()->GetDoKinks()) FindKinks(fSeeds,fEvent);
5940 //FindCurling(fSeeds, fEvent,0);
5941 if (AliTPCReconstructor::StreamLevel()>5) FindMultiMC(fSeeds, fEvent,-1); // find multi found tracks
5942 RemoveUsed2(fSeeds,0.5,0.4,20);
5943 FindSplitted(fSeeds, fEvent,0); // find multi found tracks
5944 if (AliTPCReconstructor::StreamLevel()>5) FindMultiMC(fSeeds, fEvent,0); // find multi found tracks
5947 // // refit short tracks
5949 Int_t nseed=fSeeds->GetEntriesFast();
5952 for (Int_t i=0; i<nseed; i++) {
5953 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
5955 Int_t nc=t.GetNumberOfClusters();
5957 fSeedsPool->MarkSlotFree( fSeeds->RemoveAt(i) );
5960 CookLabel(pt,0.1); //For comparison only
5961 //if ((pt->IsActive() || (pt->fRemoval==10) )&& nc>50 &&pt->GetNumberOfClusters()>0.4*pt->fNFoundable){
5962 if ((pt->IsActive() || (pt->GetRemoval()==10) )){
5964 if (fDebug>0) cerr<<found<<'\r';
5968 fSeedsPool->MarkSlotFree( fSeeds->RemoveAt(i) );
5972 //RemoveOverlap(fSeeds,0.99,7,kTRUE);
5974 //RemoveUsed(fSeeds,0.9,0.9,6);
5976 nseed=fSeeds->GetEntriesFast();
5978 for (Int_t i=0; i<nseed; i++) {
5979 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
5981 Int_t nc=t.GetNumberOfClusters();
5983 fSeedsPool->MarkSlotFree( fSeeds->RemoveAt(i) );
5987 t.CookdEdx(0.02,0.6);
5988 // CheckKinkPoint(&t,0.05);
5989 //if ((pt->IsActive() || (pt->fRemoval==10) )&& nc>50 &&pt->GetNumberOfClusters()>0.4*pt->fNFoundable){
5990 if ((pt->IsActive() || (pt->GetRemoval()==10) )){
5998 fSeedsPool->MarkSlotFree( fSeeds->RemoveAt(i) );
5999 //AliTPCseed * seed1 = ReSeed(pt,0.05,0.5,1);
6001 // FollowProlongation(*seed1,0);
6002 // Int_t n = seed1->GetNumberOfClusters();
6003 // printf("fP4\t%f\t%f\n",seed1->GetC(),pt->GetC());
6004 // printf("fN\t%d\t%d\n", seed1->GetNumberOfClusters(),pt->GetNumberOfClusters());
6007 //AliTPCseed * seed2 = ReSeed(pt,0.95,0.5,0.05);
6011 SortTracks(fSeeds, 1);
6015 PrepareForBackProlongation(fSeeds,5.);
6016 PropagateBack(fSeeds);
6017 printf("Time for back propagation: \t");timer.Print();timer.Start();
6021 PrepareForProlongation(fSeeds,5.);
6022 PropagateForard2(fSeeds);
6024 printf("Time for FORWARD propagation: \t");timer.Print();timer.Start();
6025 // RemoveUsed(fSeeds,0.7,0.7,6);
6026 //RemoveOverlap(fSeeds,0.9,7,kTRUE);
6028 nseed=fSeeds->GetEntriesFast();
6030 for (Int_t i=0; i<nseed; i++) {
6031 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
6033 Int_t nc=t.GetNumberOfClusters();
6035 fSeedsPool->MarkSlotFree( fSeeds->RemoveAt(i) );
6038 t.CookdEdx(0.02,0.6);
6039 // CookLabel(pt,0.1); //For comparison only
6040 //if ((pt->IsActive() || (pt->fRemoval==10) )&& nc>50 &&pt->GetNumberOfClusters()>0.4*pt->fNFoundable){
6041 if ((pt->IsActive() || (pt->fRemoval==10) )){
6042 cerr<<found++<<'\r';
6045 fSeedsPool->MarkSlotFree( fSeeds->RemoveAt(i) );
6050 // fNTracks = found;
6052 Info("Clusters2Tracks","Time for overlap removal, track writing and dedx cooking: \t"); timer.Print();timer.Start();
6055 // cerr<<"Number of found tracks : "<<"\t"<<found<<endl;
6056 Info("Clusters2Tracks","Number of found tracks %d",found);
6058 // UnloadClusters();
6063 void AliTPCtrackerMI::Tracking(TObjArray * arr)
6066 // tracking of the seeds
6069 fSectors = fOuterSec;
6070 ParallelTracking(arr,150,63);
6071 fSectors = fOuterSec;
6072 ParallelTracking(arr,63,0);
6075 TObjArray * AliTPCtrackerMI::Tracking(Int_t seedtype, Int_t i1, Int_t i2, Float_t cuts[4], Float_t dy, Int_t dsec)
6080 static TObjArray arrTracks;
6081 TObjArray * arr = &arrTracks;
6083 fSectors = fOuterSec;
6086 for (Int_t sec=0;sec<fkNOS;sec++){
6087 if (seedtype==3) MakeSeeds3(arr,sec,i1,i2,cuts,dy, dsec);
6088 if (seedtype==4) MakeSeeds5(arr,sec,i1,i2,cuts,dy);
6089 if (seedtype==2) MakeSeeds2(arr,sec,i1,i2,cuts,dy);
6092 Info("Tracking","\nSeeding - %d\t%d\t%d\t%d\n",seedtype,i1,i2,arr->GetEntriesFast());
6104 TObjArray * AliTPCtrackerMI::Tracking()
6108 if (AliTPCReconstructor::GetRecoParam()->GetSpecialSeeding()) return TrackingSpecial();
6111 Int_t nup=fOuterSec->GetNRows()+fInnerSec->GetNRows();
6113 TObjArray * seeds = new TObjArray;
6115 Int_t fLastSeedRowSec=AliTPCReconstructor::GetRecoParam()->GetLastSeedRowSec();
6116 Int_t gapPrim = AliTPCReconstructor::GetRecoParam()->GetSeedGapPrim();
6117 Int_t gapSec = AliTPCReconstructor::GetRecoParam()->GetSeedGapSec();
6125 Float_t fnumber = 3.0;
6126 Float_t fdensity = 3.0;
6131 for (Int_t delta = 0; delta<18; delta+=gapPrim){
6135 arr = Tracking(3,nup-1-delta,nup-1-delta-gap,cuts,-1,1);
6136 SumTracks(seeds,arr);
6137 SignClusters(seeds,fnumber,fdensity);
6139 for (Int_t i=2;i<6;i+=2){
6140 // seed high pt tracks
6143 arr = Tracking(3,nup-i-delta,nup-i-delta-gap,cuts,-1,0);
6144 SumTracks(seeds,arr);
6145 SignClusters(seeds,fnumber,fdensity);
6150 // RemoveUsed(seeds,0.9,0.9,1);
6151 // UnsignClusters();
6152 // SignClusters(seeds,fnumber,fdensity);
6156 for (Int_t delta = 20; delta<120; delta+=gapPrim){
6158 // seed high pt tracks
6162 arr = Tracking(3,nup-delta,nup-delta-gap,cuts,-1);
6163 SumTracks(seeds,arr);
6164 SignClusters(seeds,fnumber,fdensity);
6169 arr = Tracking(3,nup-delta-5,nup-delta-5-gap,cuts,-1);
6170 SumTracks(seeds,arr);
6171 SignClusters(seeds,fnumber,fdensity);
6182 Info("Tracking()","\n\nPrimary seeding\t%d\n\n",seeds->GetEntriesFast());
6186 // RemoveUsed(seeds,0.75,0.75,1);
6188 //SignClusters(seeds,fnumber,fdensity);
6197 arr = Tracking(4,nup-1,nup-1-gap,cuts,-1);
6198 SumTracks(seeds,arr);
6199 SignClusters(seeds,fnumber,fdensity);
6201 arr = Tracking(4,nup-2,nup-2-gap,cuts,-1);
6202 SumTracks(seeds,arr);
6203 SignClusters(seeds,fnumber,fdensity);
6205 arr = Tracking(4,nup-3,nup-3-gap,cuts,-1);
6206 SumTracks(seeds,arr);
6207 SignClusters(seeds,fnumber,fdensity);
6209 arr = Tracking(4,nup-5,nup-5-gap,cuts,-1);
6210 SumTracks(seeds,arr);
6211 SignClusters(seeds,fnumber,fdensity);
6213 arr = Tracking(4,nup-7,nup-7-gap,cuts,-1);
6214 SumTracks(seeds,arr);
6215 SignClusters(seeds,fnumber,fdensity);
6218 arr = Tracking(4,nup-9,nup-9-gap,cuts,-1);
6219 SumTracks(seeds,arr);
6220 SignClusters(seeds,fnumber,fdensity);
6224 for (Int_t delta = 9; delta<30; delta+=gapSec){
6230 arr = Tracking(4,nup-1-delta,nup-1-delta-gap,cuts,-1);
6231 SumTracks(seeds,arr);
6232 SignClusters(seeds,fnumber,fdensity);
6234 arr = Tracking(4,nup-3-delta,nup-5-delta-gap,cuts,4);
6235 SumTracks(seeds,arr);
6236 SignClusters(seeds,fnumber,fdensity);
6249 for (Int_t delta = 30; delta<fLastSeedRowSec; delta+=gapSec){
6255 arr = Tracking(4,nup-1-delta,nup-1-delta-gap,cuts,-1);
6256 SumTracks(seeds,arr);
6257 SignClusters(seeds,fnumber,fdensity);
6259 arr = Tracking(4,nup-5-delta,nup-5-delta-gap,cuts,5 );
6260 SumTracks(seeds,arr);
6261 SignClusters(seeds,fnumber,fdensity);
6265 Info("Tracking()","\n\nSecondary seeding\t%d\n\n",seeds->GetEntriesFast());
6276 TObjArray * AliTPCtrackerMI::TrackingSpecial()
6279 // seeding adjusted for laser and cosmic tests - short tracks with big inclination angle
6280 // no primary vertex seeding tried
6284 Int_t nup=fOuterSec->GetNRows()+fInnerSec->GetNRows();
6286 TObjArray * seeds = new TObjArray;
6291 Float_t fnumber = 3.0;
6292 Float_t fdensity = 3.0;
6295 cuts[0] = AliTPCReconstructor::GetRecoParam()->GetMaxC(); // max curvature
6296 cuts[1] = 3.5; // max tan(phi) angle for seeding
6297 cuts[2] = 3.; // not used (cut on z primary vertex)
6298 cuts[3] = 3.5; // max tan(theta) angle for seeding
6300 for (Int_t delta = 0; nup-delta-gap-1>0; delta+=3){
6302 arr = Tracking(4,nup-1-delta,nup-1-delta-gap,cuts,-1);
6303 SumTracks(seeds,arr);
6304 SignClusters(seeds,fnumber,fdensity);
6308 Info("Tracking()","\n\nSecondary seeding\t%d\n\n",seeds->GetEntriesFast());
6319 void AliTPCtrackerMI::SumTracks(TObjArray *arr1,TObjArray *&arr2)
6322 //sum tracks to common container
6323 //remove suspicious tracks
6324 // RS: Attention: supplied tracks come in the static array, don't delete them
6325 Int_t nseed = arr2->GetEntriesFast();
6326 for (Int_t i=0;i<nseed;i++){
6327 AliTPCseed *pt=(AliTPCseed*)arr2->UncheckedAt(i);
6330 // remove tracks with too big curvature
6332 if (TMath::Abs(pt->GetC())>AliTPCReconstructor::GetRecoParam()->GetMaxC()){
6333 fSeedsPool->MarkSlotFree( arr2->RemoveAt(i) );
6336 // REMOVE VERY SHORT TRACKS
6337 if (pt->GetNumberOfClusters()<20){
6338 fSeedsPool->MarkSlotFree( arr2->RemoveAt(i) );
6341 // NORMAL ACTIVE TRACK
6342 if (pt->IsActive()){
6343 arr1->AddLast(arr2->RemoveAt(i));
6346 //remove not usable tracks
6347 if (pt->GetRemoval()!=10){
6348 fSeedsPool->MarkSlotFree( arr2->RemoveAt(i) );
6352 // ENABLE ONLY ENOUGH GOOD STOPPED TRACKS
6353 if (pt->GetDensityFirst(20)>0.8 || pt->GetDensityFirst(30)>0.8 || pt->GetDensityFirst(40)>0.7)
6354 arr1->AddLast(arr2->RemoveAt(i));
6356 fSeedsPool->MarkSlotFree( arr2->RemoveAt(i) );
6360 // delete arr2; arr2 = 0; // RS: this is static array, don't delete it
6365 void AliTPCtrackerMI::ParallelTracking(TObjArray *const arr, Int_t rfirst, Int_t rlast)
6368 // try to track in parralel
6370 Int_t nseed=arr->GetEntriesFast();
6371 //prepare seeds for tracking
6372 for (Int_t i=0; i<nseed; i++) {
6373 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i), &t=*pt;
6375 if (!t.IsActive()) continue;
6376 // follow prolongation to the first layer
6377 if ( (fSectors ==fInnerSec) || (t.GetFirstPoint()-fkParam->GetNRowLow()>rfirst+1) )
6378 FollowProlongation(t, rfirst+1);
6383 for (Int_t nr=rfirst; nr>=rlast; nr--){
6384 if (nr<fInnerSec->GetNRows())
6385 fSectors = fInnerSec;
6387 fSectors = fOuterSec;
6388 // make indexes with the cluster tracks for given
6390 // find nearest cluster
6391 for (Int_t i=0; i<nseed; i++) {
6392 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i), &t=*pt;
6394 if (nr==80) pt->UpdateReference();
6395 if (!pt->IsActive()) continue;
6396 // if ( (fSectors ==fOuterSec) && (pt->fFirstPoint-fkParam->GetNRowLow())<nr) continue;
6397 if (pt->GetRelativeSector()>17) {
6400 UpdateClusters(t,nr);
6402 // prolonagate to the nearest cluster - if founded
6403 for (Int_t i=0; i<nseed; i++) {
6404 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
6406 if (!pt->IsActive()) continue;
6407 // if ((fSectors ==fOuterSec) && (pt->fFirstPoint-fkParam->GetNRowLow())<nr) continue;
6408 if (pt->GetRelativeSector()>17) {
6411 FollowToNextCluster(*pt,nr);
6416 void AliTPCtrackerMI::PrepareForBackProlongation(const TObjArray *const arr,Float_t fac) const
6420 // if we use TPC track itself we have to "update" covariance
6422 Int_t nseed= arr->GetEntriesFast();
6423 for (Int_t i=0;i<nseed;i++){
6424 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
6428 //rotate to current local system at first accepted point
6429 Int_t index = pt->GetClusterIndex2(pt->GetFirstPoint());
6430 Int_t sec = (index&0xff000000)>>24;
6432 Float_t angle1 = fInnerSec->GetAlpha()*sec+fInnerSec->GetAlphaShift();
6433 if (angle1>TMath::Pi())
6434 angle1-=2.*TMath::Pi();
6435 Float_t angle2 = pt->GetAlpha();
6437 if (TMath::Abs(angle1-angle2)>0.001){
6438 if (!pt->Rotate(angle1-angle2)) return;
6439 //angle2 = pt->GetAlpha();
6440 //pt->fRelativeSector = pt->GetAlpha()/fInnerSec->GetAlpha();
6441 //if (pt->GetAlpha()<0)
6442 // pt->fRelativeSector+=18;
6443 //sec = pt->fRelativeSector;
6452 void AliTPCtrackerMI::PrepareForProlongation(TObjArray *const arr, Float_t fac) const
6456 // if we use TPC track itself we have to "update" covariance
6458 Int_t nseed= arr->GetEntriesFast();
6459 for (Int_t i=0;i<nseed;i++){
6460 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
6463 pt->SetFirstPoint(pt->GetLastPoint());
6471 Int_t AliTPCtrackerMI::PropagateBack(const TObjArray *const arr)
6474 // make back propagation
6476 Int_t nseed= arr->GetEntriesFast();
6477 for (Int_t i=0;i<nseed;i++){
6478 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
6479 if (pt&& pt->GetKinkIndex(0)<=0) {
6480 //AliTPCseed *pt2 = new AliTPCseed(*pt);
6481 fSectors = fInnerSec;
6482 //FollowBackProlongation(*pt,fInnerSec->GetNRows()-1);
6483 //fSectors = fOuterSec;
6484 FollowBackProlongation(*pt,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1,1);
6485 //if (pt->GetNumberOfClusters()<(pt->fEsd->GetTPCclusters(0)) ){
6486 // Error("PropagateBack","Not prolonged track %d",pt->GetLabel());
6487 // FollowBackProlongation(*pt2,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1);
6490 if (pt&& pt->GetKinkIndex(0)>0) {
6491 AliESDkink * kink = fEvent->GetKink(pt->GetKinkIndex(0)-1);
6492 pt->SetFirstPoint(kink->GetTPCRow0());
6493 fSectors = fInnerSec;
6494 FollowBackProlongation(*pt,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1,1);
6502 Int_t AliTPCtrackerMI::PropagateForward2(const TObjArray *const arr)
6505 // make forward propagation
6507 Int_t nseed= arr->GetEntriesFast();
6509 for (Int_t i=0;i<nseed;i++){
6510 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
6512 FollowProlongation(*pt,0,1,1);
6521 Int_t AliTPCtrackerMI::PropagateForward()
6524 // propagate track forward
6526 Int_t nseed = fSeeds->GetEntriesFast();
6527 for (Int_t i=0;i<nseed;i++){
6528 AliTPCseed *pt = (AliTPCseed*)fSeeds->UncheckedAt(i);
6530 AliTPCseed &t = *pt;
6531 Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
6532 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
6533 if (alpha < 0. ) alpha += 2.*TMath::Pi();
6534 t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN);
6538 fSectors = fOuterSec;
6539 ParallelTracking(fSeeds,fOuterSec->GetNRows()+fInnerSec->GetNRows()-1,fInnerSec->GetNRows());
6540 fSectors = fInnerSec;
6541 ParallelTracking(fSeeds,fInnerSec->GetNRows()-1,0);
6550 Int_t AliTPCtrackerMI::PropagateBack(AliTPCseed *const pt, Int_t row0, Int_t row1)
6553 // make back propagation, in between row0 and row1
6557 fSectors = fInnerSec;
6560 if (row1<fSectors->GetNRows())
6563 r1 = fSectors->GetNRows()-1;
6565 if (row0<fSectors->GetNRows()&& r1>0 )
6566 FollowBackProlongation(*pt,r1);
6567 if (row1<=fSectors->GetNRows())
6570 r1 = row1 - fSectors->GetNRows();
6571 if (r1<=0) return 0;
6572 if (r1>=fOuterSec->GetNRows()) return 0;
6573 fSectors = fOuterSec;
6574 return FollowBackProlongation(*pt,r1);
6582 void AliTPCtrackerMI::GetShape(AliTPCseed * seed, Int_t row)
6584 // gets cluster shape
6586 AliTPCClusterParam * clparam = AliTPCcalibDB::Instance()->GetClusterParam();
6587 Float_t zdrift = TMath::Abs((fkParam->GetZLength(0)-TMath::Abs(seed->GetZ())));
6588 Int_t type = (seed->GetSector() < fkParam->GetNSector()/2) ? 0: (row>126) ? 1:2;
6589 Double_t angulary = seed->GetSnp();
6591 if (TMath::Abs(angulary)>AliTPCReconstructor::GetMaxSnpTracker()) {
6592 angulary = TMath::Sign(AliTPCReconstructor::GetMaxSnpTracker(),angulary);
6595 angulary = angulary*angulary/((1.-angulary)*(1.+angulary));
6596 Double_t angularz = seed->GetTgl()*seed->GetTgl()*(1.+angulary);
6598 Double_t sigmay = clparam->GetRMS0(0,type,zdrift,TMath::Sqrt(TMath::Abs(angulary)));
6599 Double_t sigmaz = clparam->GetRMS0(1,type,zdrift,TMath::Sqrt(TMath::Abs(angularz)));
6600 seed->SetCurrentSigmaY2(sigmay*sigmay);
6601 seed->SetCurrentSigmaZ2(sigmaz*sigmaz);
6602 // Float_t sd2 = TMath::Abs((fkParam->GetZLength(0)-TMath::Abs(seed->GetZ())))*fkParam->GetDiffL()*fkParam->GetDiffL();
6603 // // Float_t padlength = fkParam->GetPadPitchLength(seed->fSector);
6604 // Float_t padlength = GetPadPitchLength(row);
6606 // Float_t sresy = (seed->GetSector() < fkParam->GetNSector()/2) ? 0.2 :0.3;
6607 // seed->SetCurrentSigmaY2(sd2+padlength*padlength*angulary/12.+sresy*sresy);
6609 // Float_t sresz = fkParam->GetZSigma();
6610 // seed->SetCurrentSigmaZ2(sd2+padlength*padlength*angularz*angularz*(1+angulary)/12.+sresz*sresz);
6612 Float_t wy = GetSigmaY(seed);
6613 Float_t wz = GetSigmaZ(seed);
6616 if (TMath::Abs(wy/seed->fCurrentSigmaY2-1)>0.0001 || TMath::Abs(wz/seed->fCurrentSigmaZ2-1)>0.0001 ){
6617 printf("problem\n");
6624 //__________________________________________________________________________
6625 void AliTPCtrackerMI::CookLabel(AliKalmanTrack *tk, Float_t wrong) const {
6626 //--------------------------------------------------------------------
6627 //This function "cooks" a track label. If label<0, this track is fake.
6628 //--------------------------------------------------------------------
6629 AliTPCseed * t = dynamic_cast<AliTPCseed*>(tk);
6631 printf("%s:%d wrong type \n",(char*)__FILE__,__LINE__);
6635 Int_t noc=t->GetNumberOfClusters();
6637 //printf("\nnot founded prolongation\n\n\n");
6643 AliTPCclusterMI *clusters[160];
6645 for (Int_t i=0;i<160;i++) {
6652 for (i=0; i<160 && current<noc; i++) {
6654 Int_t index=t->GetClusterIndex2(i);
6655 if (index<=0) continue;
6656 if (index&0x8000) continue;
6658 //clusters[current]=GetClusterMI(index);
6659 if (t->GetClusterPointer(i)){
6660 clusters[current]=t->GetClusterPointer(i);
6666 Int_t lab=123456789;
6667 for (i=0; i<noc; i++) {
6668 AliTPCclusterMI *c=clusters[i];
6670 lab=TMath::Abs(c->GetLabel(0));
6672 for (j=0; j<noc; j++) if (lb[j]==lab || mx[j]==0) break;
6678 for (i=0; i<noc; i++) if (mx[i]>max) {max=mx[i]; lab=lb[i];}
6680 for (i=0; i<noc; i++) {
6681 AliTPCclusterMI *c=clusters[i];
6683 if (TMath::Abs(c->GetLabel(1)) == lab ||
6684 TMath::Abs(c->GetLabel(2)) == lab ) max++;
6686 if (noc<=0) { lab=-1; return;}
6687 if ((1.- Float_t(max)/(noc)) > wrong) lab=-lab;
6690 Int_t tail=Int_t(0.10*noc);
6693 for (i=1; i<160&&ind<tail; i++) {
6694 // AliTPCclusterMI *c=clusters[noc-i];
6695 AliTPCclusterMI *c=clusters[i];
6697 if (lab == TMath::Abs(c->GetLabel(0)) ||
6698 lab == TMath::Abs(c->GetLabel(1)) ||
6699 lab == TMath::Abs(c->GetLabel(2))) max++;
6702 if (max < Int_t(0.5*tail)) lab=-lab;
6709 //delete[] clusters;
6713 //__________________________________________________________________________
6714 Int_t AliTPCtrackerMI::CookLabel(AliTPCseed *const t, Float_t wrong,Int_t first, Int_t last) const {
6715 //--------------------------------------------------------------------
6716 //This function "cooks" a track label. If label<0, this track is fake.
6717 //--------------------------------------------------------------------
6718 Int_t noc=t->GetNumberOfClusters();
6720 //printf("\nnot founded prolongation\n\n\n");
6726 AliTPCclusterMI *clusters[160];
6728 for (Int_t i=0;i<160;i++) {
6735 for (i=0; i<160 && current<noc; i++) {
6736 if (i<first) continue;
6737 if (i>last) continue;
6738 Int_t index=t->GetClusterIndex2(i);
6739 if (index<=0) continue;
6740 if (index&0x8000) continue;
6742 //clusters[current]=GetClusterMI(index);
6743 if (t->GetClusterPointer(i)){
6744 clusters[current]=t->GetClusterPointer(i);
6749 //if (noc<5) return -1;
6750 Int_t lab=123456789;
6751 for (i=0; i<noc; i++) {
6752 AliTPCclusterMI *c=clusters[i];
6754 lab=TMath::Abs(c->GetLabel(0));
6756 for (j=0; j<noc; j++) if (lb[j]==lab || mx[j]==0) break;
6762 for (i=0; i<noc; i++) if (mx[i]>max) {max=mx[i]; lab=lb[i];}
6764 for (i=0; i<noc; i++) {
6765 AliTPCclusterMI *c=clusters[i];
6767 if (TMath::Abs(c->GetLabel(1)) == lab ||
6768 TMath::Abs(c->GetLabel(2)) == lab ) max++;
6770 if (noc<=0) { lab=-1; return -1;}
6771 if ((1.- Float_t(max)/(noc)) > wrong) lab=-lab;
6774 Int_t tail=Int_t(0.10*noc);
6777 for (i=1; i<160&&ind<tail; i++) {
6778 // AliTPCclusterMI *c=clusters[noc-i];
6779 AliTPCclusterMI *c=clusters[i];
6781 if (lab == TMath::Abs(c->GetLabel(0)) ||
6782 lab == TMath::Abs(c->GetLabel(1)) ||
6783 lab == TMath::Abs(c->GetLabel(2))) max++;
6786 if (max < Int_t(0.5*tail)) lab=-lab;
6789 // t->SetLabel(lab);
6793 //delete[] clusters;
6797 Int_t AliTPCtrackerMI::GetRowNumber(Double_t x[3]) const
6799 //return pad row number for given x vector
6800 Float_t phi = TMath::ATan2(x[1],x[0]);
6801 if(phi<0) phi=2.*TMath::Pi()+phi;
6802 // Get the local angle in the sector philoc
6803 const Float_t kRaddeg = 180/3.14159265358979312;
6804 Float_t phiangle = (Int_t (phi*kRaddeg/20.) + 0.5)*20./kRaddeg;
6805 Double_t localx = x[0]*TMath::Cos(phiangle)-x[1]*TMath::Sin(phiangle);
6806 return GetRowNumber(localx);
6811 void AliTPCtrackerMI::MakeESDBitmaps(AliTPCseed *t, AliESDtrack *esd)
6813 //-----------------------------------------------------------------------
6814 // Fill the cluster and sharing bitmaps of the track
6815 //-----------------------------------------------------------------------
6817 Int_t firstpoint = 0;
6818 Int_t lastpoint = 159;
6819 AliTPCTrackerPoint *point;
6820 AliTPCclusterMI *cluster;
6823 TBits clusterMap(159);
6824 TBits sharedMap(159);
6826 for (int iter=firstpoint; iter<lastpoint; iter++) {
6827 // Change to cluster pointers to see if we have a cluster at given padrow
6829 cluster = t->GetClusterPointer(iter);
6831 clusterMap.SetBitNumber(iter, kTRUE);
6832 point = t->GetTrackPoint(iter);
6833 if (point->IsShared())
6834 sharedMap.SetBitNumber(iter,kTRUE);
6836 if (t->GetClusterIndex(iter) >= 0 && (t->GetClusterIndex(iter) & 0x8000) == 0) {
6837 fitMap.SetBitNumber(iter, kTRUE);
6841 esd->SetTPCClusterMap(clusterMap);
6842 esd->SetTPCSharedMap(sharedMap);
6843 esd->SetTPCFitMap(fitMap);
6844 if (nclsf != t->GetNumberOfClusters())
6845 AliWarning(Form("Inconsistency between ncls %d and indices %d (found %d)",t->GetNumberOfClusters(),nclsf,esd->GetTPCClusterMap().CountBits()));
6848 Bool_t AliTPCtrackerMI::IsFindable(AliTPCseed & track){
6850 // return flag if there is findable cluster at given position
6853 Float_t z = track.GetZ();
6855 if (TMath::Abs(z)<(AliTPCReconstructor::GetCtgRange()*track.GetX()+kDeltaZ) &&
6856 TMath::Abs(z)<fkParam->GetZLength(0) &&
6857 (TMath::Abs(track.GetSnp())<AliTPCReconstructor::GetMaxSnpTracker()))
6863 void AliTPCtrackerMI::AddCovariance(AliTPCseed * seed){
6865 // Adding systematic error estimate to the covariance matrix
6866 // !!!! the systematic error for element 4 is in 1/cm not in pt
6868 const Double_t *param = AliTPCReconstructor::GetRecoParam()->GetSystematicError();
6870 // use only the diagonal part if not specified otherwise
6871 if (!AliTPCReconstructor::GetRecoParam()->GetUseSystematicCorrelation()) return AddCovarianceAdd(seed);
6873 Double_t *covarS= (Double_t*)seed->GetCovariance();
6874 Double_t factor[5]={1,1,1,1,1};
6875 Double_t facC = AliTracker::GetBz()*kB2C;
6876 factor[0]= TMath::Sqrt(TMath::Abs((covarS[0] + param[0]*param[0])/covarS[0]));
6877 factor[1]= TMath::Sqrt(TMath::Abs((covarS[2] + param[1]*param[1])/covarS[2]));
6878 factor[2]= TMath::Sqrt(TMath::Abs((covarS[5] + param[2]*param[2])/covarS[5]));
6879 factor[3]= TMath::Sqrt(TMath::Abs((covarS[9] + param[3]*param[3])/covarS[9]));
6880 factor[4]= TMath::Sqrt(TMath::Abs((covarS[14] + facC*facC*param[4]*param[4])/covarS[14]));
6886 for (Int_t i=0; i<5; i++){
6887 for (Int_t j=i; j<5; j++){
6888 Int_t index=seed->GetIndex(i,j);
6889 covarS[index]*=factor[i]*factor[j];
6895 void AliTPCtrackerMI::AddCovarianceAdd(AliTPCseed * seed){
6897 // Adding systematic error - as additive factor without correlation
6899 const Double_t *param = AliTPCReconstructor::GetRecoParam()->GetSystematicError();
6900 Double_t *covarIn= (Double_t*)seed->GetCovariance();
6902 for (Int_t i=0;i<15;i++) covar[i]=0;
6908 covar[0] = param[0]*param[0];
6909 covar[2] = param[1]*param[1];
6910 covar[5] = param[2]*param[2];
6911 covar[9] = param[3]*param[3];
6912 Double_t facC = AliTracker::GetBz()*kB2C;
6913 covar[14]= param[4]*param[4]*facC*facC;
6915 covar[1]=TMath::Sqrt((covar[0]*covar[2]))*covarIn[1]/TMath::Sqrt((covarIn[0]*covarIn[2]));
6917 covar[3]=TMath::Sqrt((covar[0]*covar[5]))*covarIn[3]/TMath::Sqrt((covarIn[0]*covarIn[5]));
6918 covar[4]=TMath::Sqrt((covar[2]*covar[5]))*covarIn[4]/TMath::Sqrt((covarIn[2]*covarIn[5]));
6920 covar[6]=TMath::Sqrt((covar[0]*covar[9]))*covarIn[6]/TMath::Sqrt((covarIn[0]*covarIn[9]));
6921 covar[7]=TMath::Sqrt((covar[2]*covar[9]))*covarIn[7]/TMath::Sqrt((covarIn[2]*covarIn[9]));
6922 covar[8]=TMath::Sqrt((covar[5]*covar[9]))*covarIn[8]/TMath::Sqrt((covarIn[5]*covarIn[9]));
6924 covar[10]=TMath::Sqrt((covar[0]*covar[14]))*covarIn[10]/TMath::Sqrt((covarIn[0]*covarIn[14]));
6925 covar[11]=TMath::Sqrt((covar[2]*covar[14]))*covarIn[11]/TMath::Sqrt((covarIn[2]*covarIn[14]));
6926 covar[12]=TMath::Sqrt((covar[5]*covar[14]))*covarIn[12]/TMath::Sqrt((covarIn[5]*covarIn[14]));
6927 covar[13]=TMath::Sqrt((covar[9]*covar[14]))*covarIn[13]/TMath::Sqrt((covarIn[9]*covarIn[14]));
6929 seed->AddCovariance(covar);