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>
117 #include "AliComplexCluster.h"
118 #include "AliESDEvent.h"
119 #include "AliESDtrack.h"
120 #include "AliESDVertex.h"
123 #include "AliHelix.h"
124 #include "AliRunLoader.h"
125 #include "AliTPCClustersRow.h"
126 #include "AliTPCParam.h"
127 #include "AliTPCReconstructor.h"
128 #include "AliTPCpolyTrack.h"
129 #include "AliTPCreco.h"
130 #include "AliTPCseed.h"
132 #include "AliTPCtrackerSector.h"
133 #include "AliTPCtrackerMI.h"
134 #include "TStopwatch.h"
135 #include "AliTPCReconstructor.h"
137 #include "TTreeStream.h"
138 #include "AliAlignObj.h"
139 #include "AliTrackPointArray.h"
141 #include "AliTPCcalibDB.h"
142 #include "AliTPCTransform.h"
143 #include "AliTPCClusterParam.h"
147 ClassImp(AliTPCtrackerMI)
151 class AliTPCFastMath {
154 static Double_t FastAsin(Double_t x);
156 static Double_t fgFastAsin[20000]; //lookup table for fast asin computation
159 Double_t AliTPCFastMath::fgFastAsin[20000];
160 AliTPCFastMath gAliTPCFastMath; // needed to fill the LUT
162 AliTPCFastMath::AliTPCFastMath(){
164 // initialized lookup table;
165 for (Int_t i=0;i<10000;i++){
166 fgFastAsin[2*i] = TMath::ASin(i/10000.);
167 fgFastAsin[2*i+1] = (TMath::ASin((i+1)/10000.)-fgFastAsin[2*i]);
171 Double_t AliTPCFastMath::FastAsin(Double_t x){
173 // return asin using lookup table
175 Int_t index = int(x*10000);
176 return fgFastAsin[2*index]+(x*10000.-index)*fgFastAsin[2*index+1];
179 Int_t index = int(x*10000);
180 return -(fgFastAsin[2*index]+(x*10000.-index)*fgFastAsin[2*index+1]);
182 //__________________________________________________________________
183 AliTPCtrackerMI::AliTPCtrackerMI()
205 // default constructor
208 //_____________________________________________________________________
212 Int_t AliTPCtrackerMI::UpdateTrack(AliTPCseed * track, Int_t accept){
214 //update track information using current cluster - track->fCurrentCluster
217 AliTPCclusterMI* c =track->GetCurrentCluster();
218 if (accept>0) track->SetCurrentClusterIndex1(track->GetCurrentClusterIndex1() | 0x8000); //sign not accepted clusters
220 UInt_t i = track->GetCurrentClusterIndex1();
222 Int_t sec=(i&0xff000000)>>24;
223 //Int_t row = (i&0x00ff0000)>>16;
224 track->SetRow((i&0x00ff0000)>>16);
225 track->SetSector(sec);
226 // Int_t index = i&0xFFFF;
227 if (sec>=fkParam->GetNInnerSector()) track->SetRow(track->GetRow()+fkParam->GetNRowLow());
228 track->SetClusterIndex2(track->GetRow(), i);
229 //track->fFirstPoint = row;
230 //if ( track->fLastPoint<row) track->fLastPoint =row;
231 // if (track->fRow<0 || track->fRow>160) {
232 // printf("problem\n");
234 if (track->GetFirstPoint()>track->GetRow())
235 track->SetFirstPoint(track->GetRow());
236 if (track->GetLastPoint()<track->GetRow())
237 track->SetLastPoint(track->GetRow());
240 track->SetClusterPointer(track->GetRow(),c);
243 Double_t angle2 = track->GetSnp()*track->GetSnp();
245 //SET NEW Track Point
247 if (angle2<1) //PH sometimes angle2 is very big. To be investigated...
249 angle2 = TMath::Sqrt(angle2/(1-angle2));
250 AliTPCTrackerPoint &point =*(track->GetTrackPoint(track->GetRow()));
252 point.SetSigmaY(c->GetSigmaY2()/track->GetCurrentSigmaY2());
253 point.SetSigmaZ(c->GetSigmaZ2()/track->GetCurrentSigmaZ2());
254 point.SetErrY(sqrt(track->GetErrorY2()));
255 point.SetErrZ(sqrt(track->GetErrorZ2()));
257 point.SetX(track->GetX());
258 point.SetY(track->GetY());
259 point.SetZ(track->GetZ());
260 point.SetAngleY(angle2);
261 point.SetAngleZ(track->GetTgl());
262 if (point.IsShared()){
263 track->SetErrorY2(track->GetErrorY2()*4);
264 track->SetErrorZ2(track->GetErrorZ2()*4);
268 Double_t chi2 = track->GetPredictedChi2(track->GetCurrentCluster());
270 // track->SetErrorY2(track->GetErrorY2()*1.3);
271 // track->SetErrorY2(track->GetErrorY2()+0.01);
272 // track->SetErrorZ2(track->GetErrorZ2()*1.3);
273 // track->SetErrorZ2(track->GetErrorZ2()+0.005);
275 if (accept>0) return 0;
276 if (track->GetNumberOfClusters()%20==0){
277 // if (track->fHelixIn){
278 // TClonesArray & larr = *(track->fHelixIn);
279 // Int_t ihelix = larr.GetEntriesFast();
280 // new(larr[ihelix]) AliHelix(*track) ;
283 track->SetNoCluster(0);
284 return track->Update(c,chi2,i);
289 Int_t AliTPCtrackerMI::AcceptCluster(AliTPCseed * seed, AliTPCclusterMI * cluster)
292 // decide according desired precision to accept given
293 // cluster for tracking
295 seed->GetProlongation(cluster->GetX(),yt,zt);
296 Double_t sy2=ErrY2(seed,cluster);
297 Double_t sz2=ErrZ2(seed,cluster);
299 Double_t sdistancey2 = sy2+seed->GetSigmaY2();
300 Double_t sdistancez2 = sz2+seed->GetSigmaZ2();
301 Double_t dy=seed->GetCurrentCluster()->GetY()-yt;
302 Double_t dz=seed->GetCurrentCluster()->GetY()-zt;
303 Double_t rdistancey2 = (seed->GetCurrentCluster()->GetY()-yt)*
304 (seed->GetCurrentCluster()->GetY()-yt)/sdistancey2;
305 Double_t rdistancez2 = (seed->GetCurrentCluster()->GetZ()-zt)*
306 (seed->GetCurrentCluster()->GetZ()-zt)/sdistancez2;
308 Double_t rdistance2 = rdistancey2+rdistancez2;
311 if (AliTPCReconstructor::StreamLevel()>2 && seed->GetNumberOfClusters()>20) {
312 Float_t rmsy2 = seed->GetCurrentSigmaY2();
313 Float_t rmsz2 = seed->GetCurrentSigmaZ2();
314 Float_t rmsy2p30 = seed->GetCMeanSigmaY2p30();
315 Float_t rmsz2p30 = seed->GetCMeanSigmaZ2p30();
316 Float_t rmsy2p30R = seed->GetCMeanSigmaY2p30R();
317 Float_t rmsz2p30R = seed->GetCMeanSigmaZ2p30R();
318 AliExternalTrackParam param(*seed);
319 static TVectorD gcl(3),gtr(3);
321 param.GetXYZ(gcl.GetMatrixArray());
322 cluster->GetGlobalXYZ(gclf);
323 gcl[0]=gclf[0]; gcl[1]=gclf[1]; gcl[2]=gclf[2];
326 if (AliTPCReconstructor::StreamLevel()>2) {
327 (*fDebugStreamer)<<"ErrParam"<<
340 "rmsy2p30="<<rmsy2p30<<
341 "rmsz2p30="<<rmsz2p30<<
342 "rmsy2p30R="<<rmsy2p30R<<
343 "rmsz2p30R="<<rmsz2p30R<<
347 //return 0; // temporary
348 if (rdistance2>32) return 3;
351 if ((rdistancey2>9. || rdistancez2>9.) && cluster->GetType()==0)
352 return 2; //suspisiouce - will be changed
354 if ((rdistancey2>6.25 || rdistancez2>6.25) && cluster->GetType()>0)
355 // strict cut on overlaped cluster
356 return 2; //suspisiouce - will be changed
358 if ( (rdistancey2>1. || rdistancez2>6.25 )
359 && cluster->GetType()<0){
360 seed->SetNFoundable(seed->GetNFoundable()-1);
370 //_____________________________________________________________________________
371 AliTPCtrackerMI::AliTPCtrackerMI(const AliTPCParam *par):
373 fkNIS(par->GetNInnerSector()/2),
375 fkNOS(par->GetNOuterSector()/2),
392 //---------------------------------------------------------------------
393 // The main TPC tracker constructor
394 //---------------------------------------------------------------------
395 fInnerSec=new AliTPCtrackerSector[fkNIS];
396 fOuterSec=new AliTPCtrackerSector[fkNOS];
399 for (i=0; i<fkNIS; i++) fInnerSec[i].Setup(par,0);
400 for (i=0; i<fkNOS; i++) fOuterSec[i].Setup(par,1);
403 Int_t nrowlow = par->GetNRowLow();
404 Int_t nrowup = par->GetNRowUp();
407 for (i=0;i<nrowlow;i++){
408 fXRow[i] = par->GetPadRowRadiiLow(i);
409 fPadLength[i]= par->GetPadPitchLength(0,i);
410 fYMax[i] = fXRow[i]*TMath::Tan(0.5*par->GetInnerAngle());
414 for (i=0;i<nrowup;i++){
415 fXRow[i+nrowlow] = par->GetPadRowRadiiUp(i);
416 fPadLength[i+nrowlow] = par->GetPadPitchLength(60,i);
417 fYMax[i+nrowlow] = fXRow[i+nrowlow]*TMath::Tan(0.5*par->GetOuterAngle());
420 if (AliTPCReconstructor::StreamLevel()>0) {
421 fDebugStreamer = new TTreeSRedirector("TPCdebug.root");
424 //________________________________________________________________________
425 AliTPCtrackerMI::AliTPCtrackerMI(const AliTPCtrackerMI &t):
446 //------------------------------------
447 // dummy copy constructor
448 //------------------------------------------------------------------
451 AliTPCtrackerMI & AliTPCtrackerMI::operator=(const AliTPCtrackerMI& /*r*/)
453 //------------------------------
455 //--------------------------------------------------------------
458 //_____________________________________________________________________________
459 AliTPCtrackerMI::~AliTPCtrackerMI() {
460 //------------------------------------------------------------------
461 // TPC tracker destructor
462 //------------------------------------------------------------------
469 if (fDebugStreamer) delete fDebugStreamer;
473 void AliTPCtrackerMI::FillESD(const TObjArray* arr)
477 //fill esds using updated tracks
479 // write tracks to the event
480 // store index of the track
481 Int_t nseed=arr->GetEntriesFast();
482 //FindKinks(arr,fEvent);
483 for (Int_t i=0; i<nseed; i++) {
484 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
488 if (AliTPCReconstructor::StreamLevel()>1) {
489 (*fDebugStreamer)<<"Track0"<<
493 // pt->PropagateTo(fkParam->GetInnerRadiusLow());
494 if (pt->GetKinkIndex(0)<=0){ //don't propagate daughter tracks
495 pt->PropagateTo(fkParam->GetInnerRadiusLow());
498 if (( pt->GetPoints()[2]- pt->GetPoints()[0])>5 && pt->GetPoints()[3]>0.8){
500 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
501 iotrack.SetTPCPoints(pt->GetPoints());
502 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
503 iotrack.SetV0Indexes(pt->GetV0Indexes());
504 // iotrack.SetTPCpid(pt->fTPCr);
505 //iotrack.SetTPCindex(i);
506 fEvent->AddTrack(&iotrack);
510 if ( (pt->GetNumberOfClusters()>70)&& (Float_t(pt->GetNumberOfClusters())/Float_t(pt->GetNFoundable()))>0.55) {
512 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
513 iotrack.SetTPCPoints(pt->GetPoints());
514 //iotrack.SetTPCindex(i);
515 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
516 iotrack.SetV0Indexes(pt->GetV0Indexes());
517 // iotrack.SetTPCpid(pt->fTPCr);
518 fEvent->AddTrack(&iotrack);
522 // short tracks - maybe decays
524 if ( (pt->GetNumberOfClusters()>30) && (Float_t(pt->GetNumberOfClusters())/Float_t(pt->GetNFoundable()))>0.70) {
525 Int_t found,foundable,shared;
526 pt->GetClusterStatistic(0,60,found, foundable,shared,kFALSE);
527 if ( (found>20) && (pt->GetNShared()/float(pt->GetNumberOfClusters())<0.2)){
529 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
530 //iotrack.SetTPCindex(i);
531 iotrack.SetTPCPoints(pt->GetPoints());
532 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
533 iotrack.SetV0Indexes(pt->GetV0Indexes());
534 //iotrack.SetTPCpid(pt->fTPCr);
535 fEvent->AddTrack(&iotrack);
540 if ( (pt->GetNumberOfClusters()>20) && (Float_t(pt->GetNumberOfClusters())/Float_t(pt->GetNFoundable()))>0.8) {
541 Int_t found,foundable,shared;
542 pt->GetClusterStatistic(0,60,found, foundable,shared,kFALSE);
543 if (found<20) continue;
544 if (pt->GetNShared()/float(pt->GetNumberOfClusters())>0.2) continue;
547 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
548 iotrack.SetTPCPoints(pt->GetPoints());
549 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
550 iotrack.SetV0Indexes(pt->GetV0Indexes());
551 //iotrack.SetTPCpid(pt->fTPCr);
552 //iotrack.SetTPCindex(i);
553 fEvent->AddTrack(&iotrack);
556 // short tracks - secondaties
558 if ( (pt->GetNumberOfClusters()>30) ) {
559 Int_t found,foundable,shared;
560 pt->GetClusterStatistic(128,158,found, foundable,shared,kFALSE);
561 if ( (found>20) && (pt->GetNShared()/float(pt->GetNumberOfClusters())<0.2) &&float(found)/float(foundable)>0.8){
563 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
564 iotrack.SetTPCPoints(pt->GetPoints());
565 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
566 iotrack.SetV0Indexes(pt->GetV0Indexes());
567 //iotrack.SetTPCpid(pt->fTPCr);
568 //iotrack.SetTPCindex(i);
569 fEvent->AddTrack(&iotrack);
574 if ( (pt->GetNumberOfClusters()>15)) {
575 Int_t found,foundable,shared;
576 pt->GetClusterStatistic(138,158,found, foundable,shared,kFALSE);
577 if (found<15) continue;
578 if (foundable<=0) continue;
579 if (pt->GetNShared()/float(pt->GetNumberOfClusters())>0.2) continue;
580 if (float(found)/float(foundable)<0.8) continue;
583 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
584 iotrack.SetTPCPoints(pt->GetPoints());
585 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
586 iotrack.SetV0Indexes(pt->GetV0Indexes());
587 // iotrack.SetTPCpid(pt->fTPCr);
588 //iotrack.SetTPCindex(i);
589 fEvent->AddTrack(&iotrack);
593 // >> account for suppressed tracks in the kink indices (RS)
594 int nESDtracks = fEvent->GetNumberOfTracks();
595 for (int it=nESDtracks;it--;) {
596 AliESDtrack* esdTr = fEvent->GetTrack(it);
597 if (!esdTr || !esdTr->GetKinkIndex(0)) continue;
598 for (int ik=0;ik<3;ik++) {
600 if (!(knkId=esdTr->GetKinkIndex(ik))) break; // no more kinks for this track
601 AliESDkink* kink = fEvent->GetKink(TMath::Abs(knkId)-1);
603 AliError(Form("ESDTrack%d refers to non-existing kink %d",it,TMath::Abs(knkId)-1));
606 kink->SetIndex(it, knkId<0 ? 0:1); // update track index of the kink: mother at 0, daughter at 1
609 // << account for suppressed tracks in the kink indices (RS)
611 printf("Number of filled ESDs-\t%d\n",fEvent->GetNumberOfTracks());
618 Double_t AliTPCtrackerMI::ErrY2(AliTPCseed* seed, const AliTPCclusterMI * cl){
621 // Use calibrated cluster error from OCDB
623 AliTPCClusterParam * clparam = AliTPCcalibDB::Instance()->GetClusterParam();
625 Float_t z = TMath::Abs(fkParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
626 Int_t ctype = cl->GetType();
627 Int_t type = (cl->GetRow()<63) ? 0: (cl->GetRow()>126) ? 1:2;
628 Double_t angle = seed->GetSnp()*seed->GetSnp();
629 angle = TMath::Sqrt(TMath::Abs(angle/(1.-angle)));
630 Double_t erry2 = clparam->GetError0Par(0,type, z,angle);
632 erry2+=0.5; // edge cluster
635 seed->SetErrorY2(erry2);
639 //calculate look-up table at the beginning
640 // static Bool_t ginit = kFALSE;
641 // static Float_t gnoise1,gnoise2,gnoise3;
642 // static Float_t ggg1[10000];
643 // static Float_t ggg2[10000];
644 // static Float_t ggg3[10000];
645 // static Float_t glandau1[10000];
646 // static Float_t glandau2[10000];
647 // static Float_t glandau3[10000];
649 // static Float_t gcor01[500];
650 // static Float_t gcor02[500];
651 // static Float_t gcorp[500];
655 // if (ginit==kFALSE){
656 // for (Int_t i=1;i<500;i++){
657 // Float_t rsigma = float(i)/100.;
658 // gcor02[i] = TMath::Max(0.78 +TMath::Exp(7.4*(rsigma-1.2)),0.6);
659 // gcor01[i] = TMath::Max(0.72 +TMath::Exp(3.36*(rsigma-1.2)),0.6);
660 // gcorp[i] = TMath::Max(TMath::Power((rsigma+0.5),1.5),1.2);
664 // for (Int_t i=3;i<10000;i++){
668 // Float_t amp = float(i);
669 // Float_t padlength =0.75;
670 // gnoise1 = 0.0004/padlength;
671 // Float_t nel = 0.268*amp;
672 // Float_t nprim = 0.155*amp;
673 // ggg1[i] = fkParam->GetDiffT()*fkParam->GetDiffT()*(2+0.001*nel/(padlength*padlength))/nel;
674 // glandau1[i] = (2.+0.12*nprim)*0.5* (2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
675 // if (glandau1[i]>1) glandau1[i]=1;
676 // glandau1[i]*=padlength*padlength/12.;
680 // gnoise2 = 0.0004/padlength;
682 // nprim = 0.133*amp;
683 // ggg2[i] = fkParam->GetDiffT()*fkParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
684 // glandau2[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
685 // if (glandau2[i]>1) glandau2[i]=1;
686 // glandau2[i]*=padlength*padlength/12.;
691 // gnoise3 = 0.0004/padlength;
693 // nprim = 0.133*amp;
694 // ggg3[i] = fkParam->GetDiffT()*fkParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
695 // glandau3[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
696 // if (glandau3[i]>1) glandau3[i]=1;
697 // glandau3[i]*=padlength*padlength/12.;
705 // Int_t amp = int(TMath::Abs(cl->GetQ()));
707 // seed->SetErrorY2(1.);
711 // Float_t z = TMath::Abs(fkParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
712 // Int_t ctype = cl->GetType();
713 // Float_t padlength= GetPadPitchLength(seed->GetRow());
714 // Double_t angle2 = seed->GetSnp()*seed->GetSnp();
715 // angle2 = angle2/(1-angle2);
717 // //cluster "quality"
718 // Int_t rsigmay = int(100.*cl->GetSigmaY2()/(seed->GetCurrentSigmaY2()));
721 // if (fSectors==fInnerSec){
722 // snoise2 = gnoise1;
723 // res = ggg1[amp]*z+glandau1[amp]*angle2;
724 // if (ctype==0) res *= gcor01[rsigmay];
727 // res*= gcorp[rsigmay];
731 // if (padlength<1.1){
732 // snoise2 = gnoise2;
733 // res = ggg2[amp]*z+glandau2[amp]*angle2;
734 // if (ctype==0) res *= gcor02[rsigmay];
737 // res*= gcorp[rsigmay];
741 // snoise2 = gnoise3;
742 // res = ggg3[amp]*z+glandau3[amp]*angle2;
743 // if (ctype==0) res *= gcor02[rsigmay];
746 // res*= gcorp[rsigmay];
753 // res*=2.4; // overestimate error 2 times
757 // if (res<2*snoise2)
760 // seed->SetErrorY2(res);
768 Double_t AliTPCtrackerMI::ErrZ2(AliTPCseed* seed, const AliTPCclusterMI * cl){
771 // Use calibrated cluster error from OCDB
773 AliTPCClusterParam * clparam = AliTPCcalibDB::Instance()->GetClusterParam();
775 Float_t z = TMath::Abs(fkParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
776 Int_t ctype = cl->GetType();
777 Int_t type = (cl->GetRow()<63) ? 0: (cl->GetRow()>126) ? 1:2;
779 Double_t angle2 = seed->GetSnp()*seed->GetSnp();
780 angle2 = seed->GetTgl()*seed->GetTgl()*(1+angle2/(1-angle2));
781 Double_t angle = TMath::Sqrt(TMath::Abs(angle2));
782 Double_t errz2 = clparam->GetError0Par(1,type, z,angle);
784 errz2+=0.5; // edge cluster
787 seed->SetErrorZ2(errz2);
793 // //seed->SetErrorY2(0.1);
795 // //calculate look-up table at the beginning
796 // static Bool_t ginit = kFALSE;
797 // static Float_t gnoise1,gnoise2,gnoise3;
798 // static Float_t ggg1[10000];
799 // static Float_t ggg2[10000];
800 // static Float_t ggg3[10000];
801 // static Float_t glandau1[10000];
802 // static Float_t glandau2[10000];
803 // static Float_t glandau3[10000];
805 // static Float_t gcor01[1000];
806 // static Float_t gcor02[1000];
807 // static Float_t gcorp[1000];
811 // if (ginit==kFALSE){
812 // for (Int_t i=1;i<1000;i++){
813 // Float_t rsigma = float(i)/100.;
814 // gcor02[i] = TMath::Max(0.81 +TMath::Exp(6.8*(rsigma-1.2)),0.6);
815 // gcor01[i] = TMath::Max(0.72 +TMath::Exp(2.04*(rsigma-1.2)),0.6);
816 // gcorp[i] = TMath::Max(TMath::Power((rsigma+0.5),1.5),1.2);
820 // for (Int_t i=3;i<10000;i++){
824 // Float_t amp = float(i);
825 // Float_t padlength =0.75;
826 // gnoise1 = 0.0004/padlength;
827 // Float_t nel = 0.268*amp;
828 // Float_t nprim = 0.155*amp;
829 // ggg1[i] = fkParam->GetDiffT()*fkParam->GetDiffT()*(2+0.001*nel/(padlength*padlength))/nel;
830 // glandau1[i] = (2.+0.12*nprim)*0.5* (2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
831 // if (glandau1[i]>1) glandau1[i]=1;
832 // glandau1[i]*=padlength*padlength/12.;
836 // gnoise2 = 0.0004/padlength;
838 // nprim = 0.133*amp;
839 // ggg2[i] = fkParam->GetDiffT()*fkParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
840 // glandau2[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
841 // if (glandau2[i]>1) glandau2[i]=1;
842 // glandau2[i]*=padlength*padlength/12.;
847 // gnoise3 = 0.0004/padlength;
849 // nprim = 0.133*amp;
850 // ggg3[i] = fkParam->GetDiffT()*fkParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
851 // glandau3[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
852 // if (glandau3[i]>1) glandau3[i]=1;
853 // glandau3[i]*=padlength*padlength/12.;
861 // Int_t amp = int(TMath::Abs(cl->GetQ()));
863 // seed->SetErrorY2(1.);
867 // Float_t z = TMath::Abs(fkParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
868 // Int_t ctype = cl->GetType();
869 // Float_t padlength= GetPadPitchLength(seed->GetRow());
871 // Double_t angle2 = seed->GetSnp()*seed->GetSnp();
872 // // if (angle2<0.6) angle2 = 0.6;
873 // angle2 = seed->GetTgl()*seed->GetTgl()*(1+angle2/(1-angle2));
875 // //cluster "quality"
876 // Int_t rsigmaz = int(100.*cl->GetSigmaZ2()/(seed->GetCurrentSigmaZ2()));
879 // if (fSectors==fInnerSec){
880 // snoise2 = gnoise1;
881 // res = ggg1[amp]*z+glandau1[amp]*angle2;
882 // if (ctype==0) res *= gcor01[rsigmaz];
885 // res*= gcorp[rsigmaz];
889 // if (padlength<1.1){
890 // snoise2 = gnoise2;
891 // res = ggg2[amp]*z+glandau2[amp]*angle2;
892 // if (ctype==0) res *= gcor02[rsigmaz];
895 // res*= gcorp[rsigmaz];
899 // snoise2 = gnoise3;
900 // res = ggg3[amp]*z+glandau3[amp]*angle2;
901 // if (ctype==0) res *= gcor02[rsigmaz];
904 // res*= gcorp[rsigmaz];
913 // if ((ctype<0) &&<70){
918 // if (res<2*snoise2)
920 // if (res>3) res =3;
921 // seed->SetErrorZ2(res);
929 void AliTPCtrackerMI::RotateToLocal(AliTPCseed *seed)
931 //rotate to track "local coordinata
932 Float_t x = seed->GetX();
933 Float_t y = seed->GetY();
934 Float_t ymax = x*TMath::Tan(0.5*fSectors->GetAlpha());
937 seed->SetRelativeSector((seed->GetRelativeSector()+1) % fN);
938 if (!seed->Rotate(fSectors->GetAlpha()))
940 } else if (y <-ymax) {
941 seed->SetRelativeSector((seed->GetRelativeSector()-1+fN) % fN);
942 if (!seed->Rotate(-fSectors->GetAlpha()))
950 //_____________________________________________________________________________
951 Double_t AliTPCtrackerMI::F1old(Double_t x1,Double_t y1,
952 Double_t x2,Double_t y2,
953 Double_t x3,Double_t y3) const
955 //-----------------------------------------------------------------
956 // Initial approximation of the track curvature
957 //-----------------------------------------------------------------
958 Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
959 Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
960 (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
961 Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
962 (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
964 Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
965 if ( xr*xr+yr*yr<=0.00000000000001) return 100;
966 return -xr*yr/sqrt(xr*xr+yr*yr);
971 //_____________________________________________________________________________
972 Double_t AliTPCtrackerMI::F1(Double_t x1,Double_t y1,
973 Double_t x2,Double_t y2,
974 Double_t x3,Double_t y3) const
976 //-----------------------------------------------------------------
977 // Initial approximation of the track curvature
978 //-----------------------------------------------------------------
984 Double_t det = x3*y2-x2*y3;
989 Double_t u = 0.5* (x2*(x2-x3)+y2*(y2-y3))/det;
990 Double_t x0 = x3*0.5-y3*u;
991 Double_t y0 = y3*0.5+x3*u;
992 Double_t c2 = 1/TMath::Sqrt(x0*x0+y0*y0);
998 Double_t AliTPCtrackerMI::F2(Double_t x1,Double_t y1,
999 Double_t x2,Double_t y2,
1000 Double_t x3,Double_t y3) const
1002 //-----------------------------------------------------------------
1003 // Initial approximation of the track curvature
1004 //-----------------------------------------------------------------
1010 Double_t det = x3*y2-x2*y3;
1015 Double_t u = 0.5* (x2*(x2-x3)+y2*(y2-y3))/det;
1016 Double_t x0 = x3*0.5-y3*u;
1017 Double_t y0 = y3*0.5+x3*u;
1018 Double_t c2 = 1/TMath::Sqrt(x0*x0+y0*y0);
1027 //_____________________________________________________________________________
1028 Double_t AliTPCtrackerMI::F2old(Double_t x1,Double_t y1,
1029 Double_t x2,Double_t y2,
1030 Double_t x3,Double_t y3) const
1032 //-----------------------------------------------------------------
1033 // Initial approximation of the track curvature times center of curvature
1034 //-----------------------------------------------------------------
1035 Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
1036 Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
1037 (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
1038 Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
1039 (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
1041 Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
1043 return -a/(d*y1-b)*xr/sqrt(xr*xr+yr*yr);
1046 //_____________________________________________________________________________
1047 Double_t AliTPCtrackerMI::F3(Double_t x1,Double_t y1,
1048 Double_t x2,Double_t y2,
1049 Double_t z1,Double_t z2) const
1051 //-----------------------------------------------------------------
1052 // Initial approximation of the tangent of the track dip angle
1053 //-----------------------------------------------------------------
1054 return (z1 - z2)/sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
1058 Double_t AliTPCtrackerMI::F3n(Double_t x1,Double_t y1,
1059 Double_t x2,Double_t y2,
1060 Double_t z1,Double_t z2, Double_t c) const
1062 //-----------------------------------------------------------------
1063 // Initial approximation of the tangent of the track dip angle
1064 //-----------------------------------------------------------------
1068 //angle1 = (z1-z2)*c/(TMath::ASin(c*x1-ni)-TMath::ASin(c*x2-ni));
1070 Double_t d = TMath::Sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
1071 if (TMath::Abs(d*c*0.5)>1) return 0;
1072 // Double_t angle2 = TMath::ASin(d*c*0.5);
1073 // Double_t angle2 = AliTPCFastMath::FastAsin(d*c*0.5);
1074 Double_t angle2 = (d*c*0.5>0.1)? TMath::ASin(d*c*0.5): AliTPCFastMath::FastAsin(d*c*0.5);
1076 angle2 = (z1-z2)*c/(angle2*2.);
1080 Bool_t AliTPCtrackerMI::GetProlongation(Double_t x1, Double_t x2, Double_t x[5], Double_t &y, Double_t &z)
1081 {//-----------------------------------------------------------------
1082 // This function find proloncation of a track to a reference plane x=x2.
1083 //-----------------------------------------------------------------
1087 if (TMath::Abs(x[4]*x1 - x[2]) >= 0.999) {
1091 Double_t c1=x[4]*x1 - x[2], r1=TMath::Sqrt((1.-c1)*(1.+c1));
1092 Double_t c2=x[4]*x2 - x[2], r2=TMath::Sqrt((1.-c2)*(1.+c2));
1096 Double_t dy = dx*(c1+c2)/(r1+r2);
1099 Double_t delta = x[4]*dx*(c1+c2)/(c1*r2 + c2*r1);
1101 if (TMath::Abs(delta)>0.01){
1102 dz = x[3]*TMath::ASin(delta)/x[4];
1104 dz = x[3]*AliTPCFastMath::FastAsin(delta)/x[4];
1107 //dz = x[3]*AliTPCFastMath::FastAsin(delta)/x[4];
1115 Int_t AliTPCtrackerMI::LoadClusters (TTree *tree)
1120 return LoadClusters();
1124 Int_t AliTPCtrackerMI::LoadClusters(const TObjArray *arr)
1127 // load clusters to the memory
1128 AliTPCClustersRow *clrow = 0x0;
1129 Int_t lower = arr->LowerBound();
1130 Int_t entries = arr->GetEntriesFast();
1132 for (Int_t i=lower; i<entries; i++) {
1133 clrow = (AliTPCClustersRow*) arr->At(i);
1134 if(!clrow) continue;
1135 if(!clrow->GetArray()) continue;
1139 fkParam->AdjustSectorRow(clrow->GetID(),sec,row);
1141 for (Int_t icl=0; icl<clrow->GetArray()->GetEntriesFast(); icl++){
1142 Transform((AliTPCclusterMI*)(clrow->GetArray()->At(icl)));
1145 if (clrow->GetArray()->GetEntriesFast()<=0) continue;
1146 AliTPCtrackerRow * tpcrow=0;
1149 tpcrow = &(fInnerSec[sec%fkNIS][row]);
1153 tpcrow = &(fOuterSec[(sec-fkNIS*2)%fkNOS][row]);
1154 left = (sec-fkNIS*2)/fkNOS;
1157 tpcrow->SetN1(clrow->GetArray()->GetEntriesFast());
1158 tpcrow->SetClusters1(new AliTPCclusterMI[tpcrow->GetN1()]);
1159 for (Int_t j=0;j<tpcrow->GetN1();++j)
1160 tpcrow->SetCluster1(j, *(AliTPCclusterMI*)(clrow->GetArray()->At(j)));
1163 tpcrow->SetN2(clrow->GetArray()->GetEntriesFast());
1164 tpcrow->SetClusters2(new AliTPCclusterMI[tpcrow->GetN2()]);
1165 for (Int_t j=0;j<tpcrow->GetN2();++j)
1166 tpcrow->SetCluster2(j,*(AliTPCclusterMI*)(clrow->GetArray()->At(j)));
1176 Int_t AliTPCtrackerMI::LoadClusters(const TClonesArray *arr)
1179 // load clusters to the memory from one
1182 AliTPCclusterMI *clust=0;
1183 Int_t count[72][96] = { {0} , {0} };
1185 // loop over clusters
1186 for (Int_t icl=0; icl<arr->GetEntriesFast(); icl++) {
1187 clust = (AliTPCclusterMI*)arr->At(icl);
1188 if(!clust) continue;
1189 //printf("cluster: det %d, row %d \n", clust->GetDetector(),clust->GetRow());
1191 // transform clusters
1194 // count clusters per pad row
1195 count[clust->GetDetector()][clust->GetRow()]++;
1198 // insert clusters to sectors
1199 for (Int_t icl=0; icl<arr->GetEntriesFast(); icl++) {
1200 clust = (AliTPCclusterMI*)arr->At(icl);
1201 if(!clust) continue;
1203 Int_t sec = clust->GetDetector();
1204 Int_t row = clust->GetRow();
1206 // filter overlapping pad rows needed by HLT
1207 if(sec<fkNIS*2) { //IROCs
1208 if(row == 30) continue;
1211 if(row == 27 || row == 76) continue;
1217 fInnerSec[sec%fkNIS].InsertCluster(clust, count[sec][row], fkParam);
1220 left = (sec-fkNIS*2)/fkNOS;
1221 fOuterSec[(sec-fkNIS*2)%fkNOS].InsertCluster(clust, count[sec][row], fkParam);
1225 // Load functions must be called behind LoadCluster(TClonesArray*)
1227 //LoadOuterSectors();
1228 //LoadInnerSectors();
1234 Int_t AliTPCtrackerMI::LoadClusters()
1237 // load clusters to the memory
1238 AliTPCClustersRow *clrow= new AliTPCClustersRow;
1239 clrow->SetClass("AliTPCclusterMI");
1241 clrow->GetArray()->ExpandCreateFast(10000);
1243 // TTree * tree = fClustersArray.GetTree();
1245 TTree * tree = fInput;
1246 TBranch * br = tree->GetBranch("Segment");
1247 br->SetAddress(&clrow);
1249 Int_t j=Int_t(tree->GetEntries());
1250 for (Int_t i=0; i<j; i++) {
1254 fkParam->AdjustSectorRow(clrow->GetID(),sec,row);
1255 for (Int_t icl=0; icl<clrow->GetArray()->GetEntriesFast(); icl++){
1256 Transform((AliTPCclusterMI*)(clrow->GetArray()->At(icl)));
1259 AliTPCtrackerRow * tpcrow=0;
1262 tpcrow = &(fInnerSec[sec%fkNIS][row]);
1266 tpcrow = &(fOuterSec[(sec-fkNIS*2)%fkNOS][row]);
1267 left = (sec-fkNIS*2)/fkNOS;
1270 tpcrow->SetN1(clrow->GetArray()->GetEntriesFast());
1271 tpcrow->SetClusters1(new AliTPCclusterMI[tpcrow->GetN1()]);
1272 for (Int_t k=0;k<tpcrow->GetN1();++k)
1273 tpcrow->SetCluster1(k, *(AliTPCclusterMI*)(clrow->GetArray()->At(k)));
1276 tpcrow->SetN2(clrow->GetArray()->GetEntriesFast());
1277 tpcrow->SetClusters2(new AliTPCclusterMI[tpcrow->GetN2()]);
1278 for (Int_t k=0;k<tpcrow->GetN2();++k)
1279 tpcrow->SetCluster2(k,*(AliTPCclusterMI*)(clrow->GetArray()->At(k)));
1290 void AliTPCtrackerMI::UnloadClusters()
1293 // unload clusters from the memory
1295 Int_t nrows = fOuterSec->GetNRows();
1296 for (Int_t sec = 0;sec<fkNOS;sec++)
1297 for (Int_t row = 0;row<nrows;row++){
1298 AliTPCtrackerRow* tpcrow = &(fOuterSec[sec%fkNOS][row]);
1300 // if (tpcrow->fClusters1) delete []tpcrow->fClusters1;
1301 // if (tpcrow->fClusters2) delete []tpcrow->fClusters2;
1303 tpcrow->ResetClusters();
1306 nrows = fInnerSec->GetNRows();
1307 for (Int_t sec = 0;sec<fkNIS;sec++)
1308 for (Int_t row = 0;row<nrows;row++){
1309 AliTPCtrackerRow* tpcrow = &(fInnerSec[sec%fkNIS][row]);
1311 // if (tpcrow->fClusters1) delete []tpcrow->fClusters1;
1312 //if (tpcrow->fClusters2) delete []tpcrow->fClusters2;
1314 tpcrow->ResetClusters();
1320 void AliTPCtrackerMI::FillClusterArray(TObjArray* array) const{
1322 // Filling cluster to the array - For visualization purposes
1325 nrows = fOuterSec->GetNRows();
1326 for (Int_t sec = 0;sec<fkNOS;sec++)
1327 for (Int_t row = 0;row<nrows;row++){
1328 AliTPCtrackerRow* tpcrow = &(fOuterSec[sec%fkNOS][row]);
1329 if (!tpcrow) continue;
1330 for (Int_t icl = 0;icl<tpcrow->GetN();icl++){
1331 array->AddLast((TObject*)((*tpcrow)[icl]));
1334 nrows = fInnerSec->GetNRows();
1335 for (Int_t sec = 0;sec<fkNIS;sec++)
1336 for (Int_t row = 0;row<nrows;row++){
1337 AliTPCtrackerRow* tpcrow = &(fInnerSec[sec%fkNIS][row]);
1338 if (!tpcrow) continue;
1339 for (Int_t icl = 0;icl<tpcrow->GetN();icl++){
1340 array->AddLast((TObject*)(*tpcrow)[icl]);
1346 void AliTPCtrackerMI::Transform(AliTPCclusterMI * cluster){
1350 AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ;
1352 AliFatal("Tranformations not in calibDB");
1354 transform->SetCurrentRecoParam((AliTPCRecoParam*)AliTPCReconstructor::GetRecoParam());
1355 Double_t x[3]={cluster->GetRow(),cluster->GetPad(),cluster->GetTimeBin()};
1356 Int_t i[1]={cluster->GetDetector()};
1357 transform->Transform(x,i,0,1);
1358 // if (cluster->GetDetector()%36>17){
1363 // in debug mode check the transformation
1365 if (AliTPCReconstructor::StreamLevel()>2) {
1367 cluster->GetGlobalXYZ(gx);
1368 Int_t event = (fEvent==NULL)? 0: fEvent->GetEventNumberInFile();
1369 TTreeSRedirector &cstream = *fDebugStreamer;
1370 cstream<<"Transform"<<
1381 cluster->SetX(x[0]);
1382 cluster->SetY(x[1]);
1383 cluster->SetZ(x[2]);
1389 //if (!fkParam->IsGeoRead()) fkParam->ReadGeoMatrices();
1390 TGeoHMatrix *mat = fkParam->GetClusterMatrix(cluster->GetDetector());
1392 Double_t pos[3]= {cluster->GetX(),cluster->GetY(),cluster->GetZ()};
1393 Double_t posC[3]={cluster->GetX(),cluster->GetY(),cluster->GetZ()};
1394 if (mat) mat->LocalToMaster(pos,posC);
1396 // chack Loading of Geo matrices from GeoManager - TEMPORARY FIX
1398 cluster->SetX(posC[0]);
1399 cluster->SetY(posC[1]);
1400 cluster->SetZ(posC[2]);
1403 //_____________________________________________________________________________
1404 Int_t AliTPCtrackerMI::LoadOuterSectors() {
1405 //-----------------------------------------------------------------
1406 // This function fills outer TPC sectors with clusters.
1407 //-----------------------------------------------------------------
1408 Int_t nrows = fOuterSec->GetNRows();
1410 for (Int_t sec = 0;sec<fkNOS;sec++)
1411 for (Int_t row = 0;row<nrows;row++){
1412 AliTPCtrackerRow* tpcrow = &(fOuterSec[sec%fkNOS][row]);
1413 Int_t sec2 = sec+2*fkNIS;
1415 Int_t ncl = tpcrow->GetN1();
1417 AliTPCclusterMI *c= (tpcrow->GetCluster1(ncl));
1418 index=(((sec2<<8)+row)<<16)+ncl;
1419 tpcrow->InsertCluster(c,index);
1422 ncl = tpcrow->GetN2();
1424 AliTPCclusterMI *c= (tpcrow->GetCluster2(ncl));
1425 index=((((sec2+fkNOS)<<8)+row)<<16)+ncl;
1426 tpcrow->InsertCluster(c,index);
1429 // write indexes for fast acces
1431 for (Int_t i=0;i<510;i++)
1432 tpcrow->SetFastCluster(i,-1);
1433 for (Int_t i=0;i<tpcrow->GetN();i++){
1434 Int_t zi = Int_t((*tpcrow)[i]->GetZ()+255.);
1435 tpcrow->SetFastCluster(zi,i); // write index
1438 for (Int_t i=0;i<510;i++){
1439 if (tpcrow->GetFastCluster(i)<0)
1440 tpcrow->SetFastCluster(i,last);
1442 last = tpcrow->GetFastCluster(i);
1451 //_____________________________________________________________________________
1452 Int_t AliTPCtrackerMI::LoadInnerSectors() {
1453 //-----------------------------------------------------------------
1454 // This function fills inner TPC sectors with clusters.
1455 //-----------------------------------------------------------------
1456 Int_t nrows = fInnerSec->GetNRows();
1458 for (Int_t sec = 0;sec<fkNIS;sec++)
1459 for (Int_t row = 0;row<nrows;row++){
1460 AliTPCtrackerRow* tpcrow = &(fInnerSec[sec%fkNIS][row]);
1463 Int_t ncl = tpcrow->GetN1();
1465 AliTPCclusterMI *c= (tpcrow->GetCluster1(ncl));
1466 index=(((sec<<8)+row)<<16)+ncl;
1467 tpcrow->InsertCluster(c,index);
1470 ncl = tpcrow->GetN2();
1472 AliTPCclusterMI *c= (tpcrow->GetCluster2(ncl));
1473 index=((((sec+fkNIS)<<8)+row)<<16)+ncl;
1474 tpcrow->InsertCluster(c,index);
1477 // write indexes for fast acces
1479 for (Int_t i=0;i<510;i++)
1480 tpcrow->SetFastCluster(i,-1);
1481 for (Int_t i=0;i<tpcrow->GetN();i++){
1482 Int_t zi = Int_t((*tpcrow)[i]->GetZ()+255.);
1483 tpcrow->SetFastCluster(zi,i); // write index
1486 for (Int_t i=0;i<510;i++){
1487 if (tpcrow->GetFastCluster(i)<0)
1488 tpcrow->SetFastCluster(i,last);
1490 last = tpcrow->GetFastCluster(i);
1502 //_________________________________________________________________________
1503 AliTPCclusterMI *AliTPCtrackerMI::GetClusterMI(Int_t index) const {
1504 //--------------------------------------------------------------------
1505 // Return pointer to a given cluster
1506 //--------------------------------------------------------------------
1507 if (index<0) return 0; // no cluster
1508 Int_t sec=(index&0xff000000)>>24;
1509 Int_t row=(index&0x00ff0000)>>16;
1510 Int_t ncl=(index&0x00007fff)>>00;
1512 const AliTPCtrackerRow * tpcrow=0;
1513 AliTPCclusterMI * clrow =0;
1515 if (sec<0 || sec>=fkNIS*4) {
1516 AliWarning(Form("Wrong sector %d",sec));
1521 AliTPCtrackerSector& tracksec = fInnerSec[sec%fkNIS];
1522 if (tracksec.GetNRows()<=row) return 0;
1523 tpcrow = &(tracksec[row]);
1524 if (tpcrow==0) return 0;
1527 if (tpcrow->GetN1()<=ncl) return 0;
1528 clrow = tpcrow->GetClusters1();
1531 if (tpcrow->GetN2()<=ncl) return 0;
1532 clrow = tpcrow->GetClusters2();
1536 AliTPCtrackerSector& tracksec = fOuterSec[(sec-fkNIS*2)%fkNOS];
1537 if (tracksec.GetNRows()<=row) return 0;
1538 tpcrow = &(tracksec[row]);
1539 if (tpcrow==0) return 0;
1541 if (sec-2*fkNIS<fkNOS) {
1542 if (tpcrow->GetN1()<=ncl) return 0;
1543 clrow = tpcrow->GetClusters1();
1546 if (tpcrow->GetN2()<=ncl) return 0;
1547 clrow = tpcrow->GetClusters2();
1551 return &(clrow[ncl]);
1557 Int_t AliTPCtrackerMI::FollowToNext(AliTPCseed& t, Int_t nr) {
1558 //-----------------------------------------------------------------
1559 // This function tries to find a track prolongation to next pad row
1560 //-----------------------------------------------------------------
1562 Double_t x= GetXrow(nr), ymax=GetMaxY(nr);
1565 AliTPCclusterMI *cl=0;
1566 Int_t tpcindex= t.GetClusterIndex2(nr);
1568 // update current shape info every 5 pad-row
1569 // if ( (nr%5==0) || t.GetNumberOfClusters()<2 || (t.fCurrentSigmaY2<0.0001) ){
1573 if (fIteration>0 && tpcindex>=-1){ //if we have already clusters
1575 if (tpcindex==-1) return 0; //track in dead zone
1577 cl = t.GetClusterPointer(nr);
1578 if ( (cl==0) ) cl = GetClusterMI(tpcindex);
1579 t.SetCurrentClusterIndex1(tpcindex);
1582 Int_t relativesector = ((tpcindex&0xff000000)>>24)%18; // if previously accepted cluster in different sector
1583 Float_t angle = relativesector*fSectors->GetAlpha()+fSectors->GetAlphaShift();
1585 if (angle<-TMath::Pi()) angle += 2*TMath::Pi();
1586 if (angle>=TMath::Pi()) angle -= 2*TMath::Pi();
1588 if (TMath::Abs(angle-t.GetAlpha())>0.001){
1589 Double_t rotation = angle-t.GetAlpha();
1590 t.SetRelativeSector(relativesector);
1591 if (!t.Rotate(rotation)) return 0;
1593 if (!t.PropagateTo(x)) return 0;
1595 t.SetCurrentCluster(cl);
1597 Int_t accept = AcceptCluster(&t,t.GetCurrentCluster());
1598 if ((tpcindex&0x8000)==0) accept =0;
1600 //if founded cluster is acceptible
1601 if (cl->IsUsed(11)) { // id cluster is shared inrease uncertainty
1602 t.SetErrorY2(t.GetErrorY2()+0.03);
1603 t.SetErrorZ2(t.GetErrorZ2()+0.03);
1604 t.SetErrorY2(t.GetErrorY2()*3);
1605 t.SetErrorZ2(t.GetErrorZ2()*3);
1607 t.SetNFoundable(t.GetNFoundable()+1);
1608 UpdateTrack(&t,accept);
1613 if (TMath::Abs(t.GetSnp())>AliTPCReconstructor::GetMaxSnpTracker()) return 0; // cut on angle
1614 if (fIteration>1 && IsFindable(t)){
1615 // not look for new cluster during refitting
1616 t.SetNFoundable(t.GetNFoundable()+1);
1621 // if (TMath::Abs(t.GetSnp())>0.95 || TMath::Abs(x*t.GetC()-t.GetEta())>0.95) return 0;// patch 28 fev 06
1622 Double_t y=t.GetYat(x);
1623 if (TMath::Abs(y)>ymax){
1625 t.SetRelativeSector((t.GetRelativeSector()+1) % fN);
1626 if (!t.Rotate(fSectors->GetAlpha()))
1628 } else if (y <-ymax) {
1629 t.SetRelativeSector((t.GetRelativeSector()-1+fN) % fN);
1630 if (!t.Rotate(-fSectors->GetAlpha()))
1636 if (!t.PropagateTo(x)) {
1637 if (fIteration==0) t.SetRemoval(10);
1641 Double_t z=t.GetZ();
1644 if (!IsActive(t.GetRelativeSector(),nr)) {
1646 t.SetClusterIndex2(nr,-1);
1649 //AliInfo(Form("A - Sector%d phi %f - alpha %f", t.fRelativeSector,y/x, t.GetAlpha()));
1650 Bool_t isActive = IsActive(t.GetRelativeSector(),nr);
1651 Bool_t isActive2 = (nr>=fInnerSec->GetNRows()) ? fOuterSec[t.GetRelativeSector()][nr-fInnerSec->GetNRows()].GetN()>0:fInnerSec[t.GetRelativeSector()][nr].GetN()>0;
1653 if (!isActive || !isActive2) return 0;
1655 const AliTPCtrackerRow &krow=GetRow(t.GetRelativeSector(),nr);
1656 if ( (t.GetSigmaY2()<0) || t.GetSigmaZ2()<0) return 0;
1658 Double_t roadz = 1.;
1660 if (TMath::Abs(TMath::Abs(y)-ymax)<krow.GetDeadZone()){
1662 t.SetClusterIndex2(nr,-1);
1668 // if (TMath::Abs(z)<(AliTPCReconstructor::GetCtgRange()*x+10) && TMath::Abs(z)<fkParam->GetZLength(0) && (TMath::Abs(t.GetSnp())<AliTPCReconstructor::GetMaxSnpTracker()))
1669 t.SetNFoundable(t.GetNFoundable()+1);
1675 // cl = krow.FindNearest2(y+10.,z,roady,roadz,index);
1676 cl = krow.FindNearest2(y,z,roady,roadz,index);
1677 if (cl) t.SetCurrentClusterIndex1(krow.GetIndex(index));
1680 t.SetCurrentCluster(cl);
1682 if (fIteration==2&&cl->IsUsed(10)) return 0;
1683 Int_t accept = AcceptCluster(&t,t.GetCurrentCluster());
1684 if (fIteration==2&&cl->IsUsed(11)) {
1685 t.SetErrorY2(t.GetErrorY2()+0.03);
1686 t.SetErrorZ2(t.GetErrorZ2()+0.03);
1687 t.SetErrorY2(t.GetErrorY2()*3);
1688 t.SetErrorZ2(t.GetErrorZ2()*3);
1691 if (t.fCurrentCluster->IsUsed(10)){
1696 if (t.fNShared>0.7*t.GetNumberOfClusters()) {
1702 if (accept<3) UpdateTrack(&t,accept);
1705 if ( fIteration==0 && t.GetNFoundable()*0.5 > t.GetNumberOfClusters()) t.SetRemoval(10);
1713 //_________________________________________________________________________
1714 Bool_t AliTPCtrackerMI::GetTrackPoint(Int_t index, AliTrackPoint &p ) const
1716 // Get track space point by index
1717 // return false in case the cluster doesn't exist
1718 AliTPCclusterMI *cl = GetClusterMI(index);
1719 if (!cl) return kFALSE;
1720 Int_t sector = (index&0xff000000)>>24;
1721 // Int_t row = (index&0x00ff0000)>>16;
1723 // xyz[0] = fkParam->GetPadRowRadii(sector,row);
1724 xyz[0] = cl->GetX();
1725 xyz[1] = cl->GetY();
1726 xyz[2] = cl->GetZ();
1728 fkParam->AdjustCosSin(sector,cos,sin);
1729 Float_t x = cos*xyz[0]-sin*xyz[1];
1730 Float_t y = cos*xyz[1]+sin*xyz[0];
1732 Float_t sigmaY2 = 0.027*cl->GetSigmaY2();
1733 if (sector < fkParam->GetNInnerSector()) sigmaY2 *= 2.07;
1734 Float_t sigmaZ2 = 0.066*cl->GetSigmaZ2();
1735 if (sector < fkParam->GetNInnerSector()) sigmaZ2 *= 1.77;
1736 cov[0] = sin*sin*sigmaY2;
1737 cov[1] = -sin*cos*sigmaY2;
1739 cov[3] = cos*cos*sigmaY2;
1742 p.SetXYZ(x,y,xyz[2],cov);
1743 AliGeomManager::ELayerID iLayer;
1745 if (sector < fkParam->GetNInnerSector()) {
1746 iLayer = AliGeomManager::kTPC1;
1750 iLayer = AliGeomManager::kTPC2;
1751 idet = sector - fkParam->GetNInnerSector();
1753 UShort_t volid = AliGeomManager::LayerToVolUID(iLayer,idet);
1754 p.SetVolumeID(volid);
1760 Int_t AliTPCtrackerMI::UpdateClusters(AliTPCseed& t, Int_t nr) {
1761 //-----------------------------------------------------------------
1762 // This function tries to find a track prolongation to next pad row
1763 //-----------------------------------------------------------------
1764 t.SetCurrentCluster(0);
1765 t.SetCurrentClusterIndex1(0);
1767 Double_t xt=t.GetX();
1768 Int_t row = GetRowNumber(xt)-1;
1769 Double_t ymax= GetMaxY(nr);
1771 if (row < nr) return 1; // don't prolongate if not information until now -
1772 // if (TMath::Abs(t.GetSnp())>0.9 && t.GetNumberOfClusters()>40. && fIteration!=2) {
1774 // return 0; // not prolongate strongly inclined tracks
1776 // if (TMath::Abs(t.GetSnp())>0.95) {
1778 // return 0; // not prolongate strongly inclined tracks
1779 // }// patch 28 fev 06
1781 Double_t x= GetXrow(nr);
1783 //t.PropagateTo(x+0.02);
1784 //t.PropagateTo(x+0.01);
1785 if (!t.PropagateTo(x)){
1792 if (TMath::Abs(y)>ymax){
1794 t.SetRelativeSector((t.GetRelativeSector()+1) % fN);
1795 if (!t.Rotate(fSectors->GetAlpha()))
1797 } else if (y <-ymax) {
1798 t.SetRelativeSector((t.GetRelativeSector()-1+fN) % fN);
1799 if (!t.Rotate(-fSectors->GetAlpha()))
1802 // if (!t.PropagateTo(x)){
1809 if (TMath::Abs(t.GetSnp())>AliTPCReconstructor::GetMaxSnpTracker()) return 0;
1811 if (!IsActive(t.GetRelativeSector(),nr)) {
1813 t.SetClusterIndex2(nr,-1);
1816 //AliInfo(Form("A - Sector%d phi %f - alpha %f", t.fRelativeSector,y/x, t.GetAlpha()));
1818 AliTPCtrackerRow &krow=GetRow(t.GetRelativeSector(),nr);
1820 if (TMath::Abs(TMath::Abs(y)-ymax)<krow.GetDeadZone()){
1822 t.SetClusterIndex2(nr,-1);
1828 // if (TMath::Abs(t.GetZ())<(AliTPCReconstructor::GetCtgRange()*t.GetX()+10) && (TMath::Abs(t.GetSnp())<AliTPCReconstructor::GetMaxSnpTracker()))
1829 if (IsFindable(t)) t.SetNFoundable(t.GetNFoundable()+1);
1835 if ( (nr%2==0) || t.GetNumberOfClusters()<2){
1836 // t.fCurrentSigmaY = GetSigmaY(&t);
1837 //t.fCurrentSigmaZ = GetSigmaZ(&t);
1841 AliTPCclusterMI *cl=0;
1844 Double_t roady = 1.;
1845 Double_t roadz = 1.;
1849 index = t.GetClusterIndex2(nr);
1850 if ( (index>0) && (index&0x8000)==0){
1851 cl = t.GetClusterPointer(nr);
1852 if ( (cl==0) && (index>0)) cl = GetClusterMI(index);
1853 t.SetCurrentClusterIndex1(index);
1855 t.SetCurrentCluster(cl);
1861 // if (index<0) return 0;
1862 UInt_t uindex = TMath::Abs(index);
1865 //cl = krow.FindNearest2(y+10,z,roady,roadz,uindex);
1866 cl = krow.FindNearest2(y,z,roady,roadz,uindex);
1869 if (cl) t.SetCurrentClusterIndex1(krow.GetIndex(uindex));
1870 t.SetCurrentCluster(cl);
1876 Int_t AliTPCtrackerMI::FollowToNextCluster(AliTPCseed & t, Int_t nr) {
1877 //-----------------------------------------------------------------
1878 // This function tries to find a track prolongation to next pad row
1879 //-----------------------------------------------------------------
1881 //update error according neighborhoud
1883 if (t.GetCurrentCluster()) {
1885 Int_t accept = AcceptCluster(&t,t.GetCurrentCluster());
1887 if (t.GetCurrentCluster()->IsUsed(10)){
1892 t.SetNShared(t.GetNShared()+1);
1893 if (t.GetNShared()>0.7*t.GetNumberOfClusters()) {
1898 if (fIteration>0) accept = 0;
1899 if (accept<3) UpdateTrack(&t,accept);
1903 if ( ( (t.GetSigmaY2()+t.GetSigmaZ2())>0.16)&& t.GetNumberOfClusters()>18) t.SetRemoval(10);
1904 if ( t.GetChi2()/t.GetNumberOfClusters()>6 &&t.GetNumberOfClusters()>18) t.SetRemoval(10);
1906 if (( (t.GetNFoundable()*0.5 > t.GetNumberOfClusters()) || t.GetNoCluster()>15)) t.SetRemoval(10);
1914 //_____________________________________________________________________________
1915 Int_t AliTPCtrackerMI::FollowProlongation(AliTPCseed& t, Int_t rf, Int_t step) {
1916 //-----------------------------------------------------------------
1917 // This function tries to find a track prolongation.
1918 //-----------------------------------------------------------------
1919 Double_t xt=t.GetX();
1921 Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
1922 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
1923 if (alpha < 0. ) alpha += 2.*TMath::Pi();
1925 t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN);
1927 Int_t first = GetRowNumber(xt)-1;
1928 for (Int_t nr= first; nr>=rf; nr-=step) {
1930 if (t.GetKinkIndexes()[0]>0){
1931 for (Int_t i=0;i<3;i++){
1932 Int_t index = t.GetKinkIndexes()[i];
1933 if (index==0) break;
1934 if (index<0) continue;
1936 AliKink * kink = (AliKink*)fEvent->GetKink(index-1);
1938 printf("PROBLEM\n");
1941 Int_t kinkrow = kink->GetTPCRow0()+2+Int_t(0.5/(0.05+kink->GetAngle(2)));
1943 AliExternalTrackParam paramd(t);
1944 kink->SetDaughter(paramd);
1945 kink->SetStatus(2,5);
1952 if (nr==80) t.UpdateReference();
1953 if (nr<fInnerSec->GetNRows())
1954 fSectors = fInnerSec;
1956 fSectors = fOuterSec;
1957 if (FollowToNext(t,nr)==0)
1970 Int_t AliTPCtrackerMI::FollowBackProlongation(AliTPCseed& t, Int_t rf) {
1971 //-----------------------------------------------------------------
1972 // This function tries to find a track prolongation.
1973 //-----------------------------------------------------------------
1975 Double_t xt=t.GetX();
1976 Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
1977 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
1978 if (alpha < 0. ) alpha += 2.*TMath::Pi();
1979 t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN);
1981 Int_t first = t.GetFirstPoint();
1982 if (first<GetRowNumber(xt)+1) first = GetRowNumber(xt)+1;
1984 if (first<0) first=0;
1985 for (Int_t nr=first; nr<=rf; nr++) {
1986 // if ( (TMath::Abs(t.GetSnp())>0.95)) break;//patch 28 fev 06
1987 if (t.GetKinkIndexes()[0]<0){
1988 for (Int_t i=0;i<3;i++){
1989 Int_t index = t.GetKinkIndexes()[i];
1990 if (index==0) break;
1991 if (index>0) continue;
1992 index = TMath::Abs(index);
1993 AliKink * kink = (AliKink*)fEvent->GetKink(index-1);
1995 printf("PROBLEM\n");
1998 Int_t kinkrow = kink->GetTPCRow0()-2-Int_t(0.5/(0.05+kink->GetAngle(2)));
2000 AliExternalTrackParam paramm(t);
2001 kink->SetMother(paramm);
2002 kink->SetStatus(2,1);
2009 if (nr<fInnerSec->GetNRows())
2010 fSectors = fInnerSec;
2012 fSectors = fOuterSec;
2023 Float_t AliTPCtrackerMI::OverlapFactor(AliTPCseed * s1, AliTPCseed * s2, Int_t &sum1, Int_t & sum2)
2031 Float_t dz2 =(s1->GetZ() - s2->GetZ());
2034 Float_t dy2 =TMath::Abs((s1->GetY() - s2->GetY()));
2036 Float_t distance = TMath::Sqrt(dz2+dy2);
2037 if (distance>4.) return 0; // if there are far away - not overlap - to reduce combinatorics
2040 Int_t firstpoint = TMath::Min(s1->GetFirstPoint(),s2->GetFirstPoint());
2041 Int_t lastpoint = TMath::Max(s1->GetLastPoint(),s2->GetLastPoint());
2046 if (firstpoint>lastpoint) {
2047 firstpoint =lastpoint;
2052 for (Int_t i=firstpoint-1;i<lastpoint+1;i++){
2053 if (s1->GetClusterIndex2(i)>0) sum1++;
2054 if (s2->GetClusterIndex2(i)>0) sum2++;
2055 if (s1->GetClusterIndex2(i)==s2->GetClusterIndex2(i) && s1->GetClusterIndex2(i)>0) {
2059 if (sum<5) return 0;
2061 Float_t summin = TMath::Min(sum1+1,sum2+1);
2062 Float_t ratio = (sum+1)/Float_t(summin);
2066 void AliTPCtrackerMI::SignShared(AliTPCseed * s1, AliTPCseed * s2)
2070 Float_t thetaCut = 0.2;//+10.*TMath::Sqrt(s1->GetSigmaTglZ()+ s2->GetSigmaTglZ());
2071 if (TMath::Abs(s1->GetTgl()-s2->GetTgl())>thetaCut) return;
2072 Float_t minCl = TMath::Min(s1->GetNumberOfClusters(),s2->GetNumberOfClusters());
2073 Int_t cutN0 = TMath::Max(5,TMath::Nint(0.1*minCl));
2078 //Int_t firstpoint = TMath::Max(s1->GetFirstPoint(),s2->GetFirstPoint());
2079 //Int_t lastpoint = TMath::Min(s1->GetLastPoint(),s2->GetLastPoint());
2080 Int_t firstpoint = 0;
2081 Int_t lastpoint = 160;
2083 // if (firstpoint>=lastpoint-5) return;;
2085 for (Int_t i=firstpoint;i<lastpoint;i++){
2086 // if ( (s1->GetClusterIndex2(i)&0xFFFF8FFF)==(s2->GetClusterIndex2(i)&0xFFFF8FFF) && s1->GetClusterIndex2(i)>0) {
2087 if ( (s1->GetClusterIndex2(i))==(s2->GetClusterIndex2(i)) && s1->GetClusterIndex2(i)>0) {
2089 s1->SetSharedMapBit(i, kTRUE);
2090 s2->SetSharedMapBit(i, kTRUE);
2092 if (s1->GetClusterIndex2(i)>0)
2093 s1->SetClusterMapBit(i, kTRUE);
2095 if (sumshared>cutN0){
2098 for (Int_t i=firstpoint;i<lastpoint;i++){
2099 // if ( (s1->GetClusterIndex2(i)&0xFFFF8FFF)==(s2->GetClusterIndex2(i)&0xFFFF8FFF) && s1->GetClusterIndex2(i)>0) {
2100 if ( (s1->GetClusterIndex2(i))==(s2->GetClusterIndex2(i)) && s1->GetClusterIndex2(i)>0) {
2101 AliTPCTrackerPoint *p1 = s1->GetTrackPoint(i);
2102 AliTPCTrackerPoint *p2 = s2->GetTrackPoint(i);;
2103 if (s1->IsActive()&&s2->IsActive()){
2104 p1->SetShared(kTRUE);
2105 p2->SetShared(kTRUE);
2111 if (sumshared>cutN0){
2112 for (Int_t i=0;i<4;i++){
2113 if (s1->GetOverlapLabel(3*i)==0){
2114 s1->SetOverlapLabel(3*i, s2->GetLabel());
2115 s1->SetOverlapLabel(3*i+1,sumshared);
2116 s1->SetOverlapLabel(3*i+2,s2->GetUniqueID());
2120 for (Int_t i=0;i<4;i++){
2121 if (s2->GetOverlapLabel(3*i)==0){
2122 s2->SetOverlapLabel(3*i, s1->GetLabel());
2123 s2->SetOverlapLabel(3*i+1,sumshared);
2124 s2->SetOverlapLabel(3*i+2,s1->GetUniqueID());
2131 void AliTPCtrackerMI::SignShared(TObjArray * arr)
2134 //sort trackss according sectors
2136 for (Int_t i=0; i<arr->GetEntriesFast(); i++) {
2137 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2139 //if (pt) RotateToLocal(pt);
2143 arr->Sort(); // sorting according relative sectors
2144 arr->Expand(arr->GetEntries());
2147 Int_t nseed=arr->GetEntriesFast();
2148 for (Int_t i=0; i<nseed; i++) {
2149 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2151 for (Int_t j=0;j<=12;j++){
2152 pt->SetOverlapLabel(j,0);
2155 for (Int_t i=0; i<nseed; i++) {
2156 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2158 if (pt->GetRemoval()>10) continue;
2159 for (Int_t j=i+1; j<nseed; j++){
2160 AliTPCseed *pt2=(AliTPCseed*)arr->UncheckedAt(j);
2161 if (TMath::Abs(pt->GetRelativeSector()-pt2->GetRelativeSector())>1) continue;
2163 if (pt2->GetRemoval()<=10) {
2164 //if ( TMath::Abs(pt->GetRelativeSector()-pt2->GetRelativeSector())>0) break;
2172 void AliTPCtrackerMI::SortTracks(TObjArray * arr, Int_t mode) const
2175 //sort tracks in array according mode criteria
2176 Int_t nseed = arr->GetEntriesFast();
2177 for (Int_t i=0; i<nseed; i++) {
2178 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2189 void AliTPCtrackerMI::RemoveUsed2(TObjArray * arr, Float_t factor1, Float_t factor2, Int_t minimal)
2192 // Loop over all tracks and remove overlaped tracks (with lower quality)
2194 // 1. Unsign clusters
2195 // 2. Sort tracks according quality
2196 // Quality is defined by the number of cluster between first and last points
2198 // 3. Loop over tracks - decreasing quality order
2199 // a.) remove - If the fraction of shared cluster less than factor (1- n or 2)
2200 // b.) remove - If the minimal number of clusters less than minimal and not ITS
2201 // c.) if track accepted - sign clusters
2203 //Called in - AliTPCtrackerMI::Clusters2Tracks()
2204 // - AliTPCtrackerMI::PropagateBack()
2205 // - AliTPCtrackerMI::RefitInward()
2208 // factor1 - factor for constrained
2209 // factor2 - for non constrained tracks
2210 // if (Float_t(shared+1)/Float_t(found+1)>factor) - DELETE
2214 Int_t nseed = arr->GetEntriesFast();
2215 Float_t * quality = new Float_t[nseed];
2216 Int_t * indexes = new Int_t[nseed];
2220 for (Int_t i=0; i<nseed; i++) {
2221 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2226 pt->UpdatePoints(); //select first last max dens points
2227 Float_t * points = pt->GetPoints();
2228 if (points[3]<0.8) quality[i] =-1;
2229 quality[i] = (points[2]-points[0])+pt->GetNumberOfClusters();
2230 //prefer high momenta tracks if overlaps
2231 quality[i] *= TMath::Sqrt(TMath::Abs(pt->Pt())+0.5);
2233 TMath::Sort(nseed,quality,indexes);
2236 for (Int_t itrack=0; itrack<nseed; itrack++) {
2237 Int_t trackindex = indexes[itrack];
2238 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(trackindex);
2241 if (quality[trackindex]<0){
2243 delete arr->RemoveAt(trackindex);
2246 arr->RemoveAt(trackindex);
2252 Int_t first = Int_t(pt->GetPoints()[0]);
2253 Int_t last = Int_t(pt->GetPoints()[2]);
2254 Double_t factor = (pt->GetBConstrain()) ? factor1: factor2;
2256 Int_t found,foundable,shared;
2257 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
2258 // pt->GetClusterStatistic(0,160, found, foundable,shared,kFALSE);
2259 Bool_t itsgold =kFALSE;
2262 if (pt->GetESD()->GetITSclusters(dummy)>4) itsgold= kTRUE;
2266 if (Float_t(shared+1)/Float_t(found+1)>factor){
2267 if (pt->GetKinkIndexes()[0]!=0) continue; //don't remove tracks - part of the kinks
2268 if( AliTPCReconstructor::StreamLevel()>3){
2269 TTreeSRedirector &cstream = *fDebugStreamer;
2270 cstream<<"RemoveUsed"<<
2271 "iter="<<fIteration<<
2275 delete arr->RemoveAt(trackindex);
2278 if (pt->GetNumberOfClusters()<50&&(found-0.5*shared)<minimal){ //remove short tracks
2279 if (pt->GetKinkIndexes()[0]!=0) continue; //don't remove tracks - part of the kinks
2280 if( AliTPCReconstructor::StreamLevel()>3){
2281 TTreeSRedirector &cstream = *fDebugStreamer;
2282 cstream<<"RemoveShort"<<
2283 "iter="<<fIteration<<
2287 delete arr->RemoveAt(trackindex);
2293 //if (sharedfactor>0.4) continue;
2294 if (pt->GetKinkIndexes()[0]>0) continue;
2295 //Remove tracks with undefined properties - seems
2296 if (pt->GetSigmaY2()<kAlmost0) continue; // ? what is the origin ?
2298 for (Int_t i=first; i<last; i++) {
2299 Int_t index=pt->GetClusterIndex2(i);
2300 // if (index<0 || index&0x8000 ) continue;
2301 if (index<0 || index&0x8000 ) continue;
2302 AliTPCclusterMI *c= pt->GetClusterPointer(i);
2309 Info("RemoveUsed","\n*****\nNumber of good tracks after shared removal\t%d\n",fNtracks);
2315 void AliTPCtrackerMI::UnsignClusters()
2318 // loop over all clusters and unsign them
2321 for (Int_t sec=0;sec<fkNIS;sec++){
2322 for (Int_t row=0;row<fInnerSec->GetNRows();row++){
2323 AliTPCclusterMI *cl = fInnerSec[sec][row].GetClusters1();
2324 for (Int_t icl =0;icl< fInnerSec[sec][row].GetN1();icl++)
2325 // if (cl[icl].IsUsed(10))
2327 cl = fInnerSec[sec][row].GetClusters2();
2328 for (Int_t icl =0;icl< fInnerSec[sec][row].GetN2();icl++)
2329 //if (cl[icl].IsUsed(10))
2334 for (Int_t sec=0;sec<fkNOS;sec++){
2335 for (Int_t row=0;row<fOuterSec->GetNRows();row++){
2336 AliTPCclusterMI *cl = fOuterSec[sec][row].GetClusters1();
2337 for (Int_t icl =0;icl< fOuterSec[sec][row].GetN1();icl++)
2338 //if (cl[icl].IsUsed(10))
2340 cl = fOuterSec[sec][row].GetClusters2();
2341 for (Int_t icl =0;icl< fOuterSec[sec][row].GetN2();icl++)
2342 //if (cl[icl].IsUsed(10))
2351 void AliTPCtrackerMI::SignClusters(const TObjArray * arr, Float_t fnumber, Float_t fdensity)
2354 //sign clusters to be "used"
2356 // snumber and sdensity sign number of sigmas - bellow mean value to be accepted
2357 // loop over "primaries"
2371 Int_t nseed = arr->GetEntriesFast();
2372 for (Int_t i=0; i<nseed; i++) {
2373 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2377 if (!(pt->IsActive())) continue;
2378 Float_t dens = pt->GetNumberOfClusters()/Float_t(pt->GetNFoundable());
2379 if ( (dens>0.7) && (pt->GetNumberOfClusters()>70)){
2381 sumdens2+= dens*dens;
2382 sumn += pt->GetNumberOfClusters();
2383 sumn2 += pt->GetNumberOfClusters()*pt->GetNumberOfClusters();
2384 Float_t chi2 = pt->GetChi2()/pt->GetNumberOfClusters();
2387 sumchi2 +=chi2*chi2;
2392 Float_t mdensity = 0.9;
2393 Float_t meann = 130;
2394 Float_t meanchi = 1;
2395 Float_t sdensity = 0.1;
2396 Float_t smeann = 10;
2397 Float_t smeanchi =0.4;
2401 mdensity = sumdens/sum;
2403 meanchi = sumchi/sum;
2405 sdensity = sumdens2/sum-mdensity*mdensity;
2407 sdensity = TMath::Sqrt(sdensity);
2411 smeann = sumn2/sum-meann*meann;
2413 smeann = TMath::Sqrt(smeann);
2417 smeanchi = sumchi2/sum - meanchi*meanchi;
2419 smeanchi = TMath::Sqrt(smeanchi);
2425 //REMOVE SHORT DELTAS or tracks going out of sensitive volume of TPC
2427 for (Int_t i=0; i<nseed; i++) {
2428 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2432 if (pt->GetBSigned()) continue;
2433 if (pt->GetBConstrain()) continue;
2434 //if (!(pt->IsActive())) continue;
2436 Int_t found,foundable,shared;
2437 pt->GetClusterStatistic(0,160,found, foundable,shared);
2438 if (shared/float(found)>0.3) {
2439 if (shared/float(found)>0.9 ){
2440 //delete arr->RemoveAt(i);
2445 Bool_t isok =kFALSE;
2446 if ( (pt->GetNShared()/pt->GetNumberOfClusters()<0.5) &&pt->GetNumberOfClusters()>60)
2448 if ((TMath::Abs(1/pt->GetC())<100.) && (pt->GetNShared()/pt->GetNumberOfClusters()<0.7))
2450 if (TMath::Abs(pt->GetZ()/pt->GetX())>1.1)
2452 if ( (TMath::Abs(pt->GetSnp()>0.7) && pt->GetD(0,0)>60.))
2456 for (Int_t j=0; j<160; ++j) {
2457 Int_t index=pt->GetClusterIndex2(j);
2458 if (index<0) continue;
2459 AliTPCclusterMI *c= pt->GetClusterPointer(j);
2461 //if (!(c->IsUsed(10))) c->Use();
2468 Double_t maxchi = meanchi+2.*smeanchi;
2470 for (Int_t i=0; i<nseed; i++) {
2471 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2475 //if (!(pt->IsActive())) continue;
2476 if (pt->GetBSigned()) continue;
2477 Double_t chi = pt->GetChi2()/pt->GetNumberOfClusters();
2478 if (chi>maxchi) continue;
2481 Float_t dens = pt->GetNumberOfClusters()/Float_t(pt->GetNFoundable());
2483 //sign only tracks with enoug big density at the beginning
2485 if ((pt->GetDensityFirst(40)<0.75) && pt->GetNumberOfClusters()<meann) continue;
2488 Double_t mindens = TMath::Max(double(mdensity-sdensity*fdensity*bfactor),0.65);
2489 Double_t minn = TMath::Max(Int_t(meann-fnumber*smeann*bfactor),50);
2491 // if (pt->fBConstrain) mindens = TMath::Max(mdensity-sdensity*fdensity*bfactor,0.65);
2492 if ( (pt->GetRemoval()==10) && (pt->GetSnp()>0.8)&&(dens>mindens))
2495 if ((dens>mindens && pt->GetNumberOfClusters()>minn) && chi<maxchi ){
2496 //Int_t noc=pt->GetNumberOfClusters();
2497 pt->SetBSigned(kTRUE);
2498 for (Int_t j=0; j<160; ++j) {
2500 Int_t index=pt->GetClusterIndex2(j);
2501 if (index<0) continue;
2502 AliTPCclusterMI *c= pt->GetClusterPointer(j);
2504 // if (!(c->IsUsed(10))) c->Use();
2509 // gLastCheck = nseed;
2517 void AliTPCtrackerMI::StopNotActive(const TObjArray * arr, Int_t row0, Float_t th0, Float_t th1, Float_t th2) const
2519 // stop not active tracks
2520 // take th1 as threshold for number of founded to number of foundable on last 10 active rows
2521 // take th2 as threshold for number of founded to number of foundable on last 20 active rows
2522 Int_t nseed = arr->GetEntriesFast();
2524 for (Int_t i=0; i<nseed; i++) {
2525 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2529 if (!(pt->IsActive())) continue;
2530 StopNotActive(pt,row0,th0, th1,th2);
2536 void AliTPCtrackerMI::StopNotActive(AliTPCseed * seed, Int_t row0, Float_t th0, Float_t th1,
2539 // stop not active tracks
2540 // take th1 as threshold for number of founded to number of foundable on last 10 active rows
2541 // take th2 as threshold for number of founded to number of foundable on last 20 active rows
2544 Int_t foundable = 0;
2545 Int_t maxindex = seed->GetLastPoint(); //last foundable row
2546 if (seed->GetNFoundable()*th0 > seed->GetNumberOfClusters()) {
2547 seed->Desactivate(10) ;
2551 for (Int_t i=row0; i<maxindex; i++){
2552 Int_t index = seed->GetClusterIndex2(i);
2553 if (index!=-1) foundable++;
2555 if (foundable<=30) sumgood1++;
2556 if (foundable<=50) {
2563 if (foundable>=30.){
2564 if (sumgood1<(th1*30.)) seed->Desactivate(10);
2567 if (sumgood2<(th2*50.)) seed->Desactivate(10);
2571 Int_t AliTPCtrackerMI::RefitInward(AliESDEvent *event)
2574 // back propagation of ESD tracks
2577 const Int_t kMaxFriendTracks=2000;
2581 //PrepareForProlongation(fSeeds,1);
2582 PropagateForward2(fSeeds);
2583 RemoveUsed2(fSeeds,0.4,0.4,20);
2585 TObjArray arraySeed(fSeeds->GetEntries());
2586 for (Int_t i=0;i<fSeeds->GetEntries();i++) {
2587 arraySeed.AddAt(fSeeds->At(i),i);
2589 SignShared(&arraySeed);
2590 // FindCurling(fSeeds, event,2); // find multi found tracks
2591 FindSplitted(fSeeds, event,2); // find multi found tracks
2592 if (AliTPCReconstructor::StreamLevel()>5) FindMultiMC(fSeeds, fEvent,2); // find multi found tracks
2595 Int_t nseed = fSeeds->GetEntriesFast();
2596 for (Int_t i=0;i<nseed;i++){
2597 AliTPCseed * seed = (AliTPCseed*) fSeeds->UncheckedAt(i);
2598 if (!seed) continue;
2599 if (seed->GetKinkIndex(0)>0) UpdateKinkQualityD(seed); // update quality informations for kinks
2600 AliESDtrack *esd=event->GetTrack(i);
2602 if (seed->GetNumberOfClusters()<60 && seed->GetNumberOfClusters()<(esd->GetTPCclusters(0) -5)*0.8){
2603 AliExternalTrackParam paramIn;
2604 AliExternalTrackParam paramOut;
2605 Int_t ncl = seed->RefitTrack(seed,¶mIn,¶mOut);
2606 if (AliTPCReconstructor::StreamLevel()>2) {
2607 (*fDebugStreamer)<<"RecoverIn"<<
2611 "pout.="<<¶mOut<<
2616 seed->Set(paramIn.GetX(),paramIn.GetAlpha(),paramIn.GetParameter(),paramIn.GetCovariance());
2617 seed->SetNumberOfClusters(ncl);
2621 seed->PropagateTo(fkParam->GetInnerRadiusLow());
2622 seed->UpdatePoints();
2623 AddCovariance(seed);
2625 seed->CookdEdx(0.02,0.6);
2626 CookLabel(seed,0.1); //For comparison only
2627 esd->SetTPCClusterMap(seed->GetClusterMap());
2628 esd->SetTPCSharedMap(seed->GetSharedMap());
2630 if (AliTPCReconstructor::StreamLevel()>1 && seed!=0&&esd!=0) {
2631 TTreeSRedirector &cstream = *fDebugStreamer;
2638 if (seed->GetNumberOfClusters()>15){
2639 esd->UpdateTrackParams(seed,AliESDtrack::kTPCrefit);
2640 esd->SetTPCPoints(seed->GetPoints());
2641 esd->SetTPCPointsF(seed->GetNFoundable());
2642 Int_t ndedx = seed->GetNCDEDX(0);
2643 Float_t sdedx = seed->GetSDEDX(0);
2644 Float_t dedx = seed->GetdEdx();
2645 esd->SetTPCsignal(dedx, sdedx, ndedx);
2647 // add seed to the esd track in Calib level
2649 Bool_t storeFriend = gRandom->Rndm()<(kMaxFriendTracks)/Float_t(nseed);
2650 if (AliTPCReconstructor::StreamLevel()>0 &&storeFriend){
2651 AliTPCseed * seedCopy = new AliTPCseed(*seed, kTRUE);
2652 esd->AddCalibObject(seedCopy);
2657 //printf("problem\n");
2660 //FindKinks(fSeeds,event);
2661 Info("RefitInward","Number of refitted tracks %d",ntracks);
2666 Int_t AliTPCtrackerMI::PropagateBack(AliESDEvent *event)
2669 // back propagation of ESD tracks
2675 PropagateBack(fSeeds);
2676 RemoveUsed2(fSeeds,0.4,0.4,20);
2677 //FindCurling(fSeeds, fEvent,1);
2678 FindSplitted(fSeeds, event,1); // find multi found tracks
2679 if (AliTPCReconstructor::StreamLevel()>5) FindMultiMC(fSeeds, fEvent,1); // find multi found tracks
2682 Int_t nseed = fSeeds->GetEntriesFast();
2684 for (Int_t i=0;i<nseed;i++){
2685 AliTPCseed * seed = (AliTPCseed*) fSeeds->UncheckedAt(i);
2686 if (!seed) continue;
2687 if (seed->GetKinkIndex(0)<0) UpdateKinkQualityM(seed); // update quality informations for kinks
2688 seed->UpdatePoints();
2689 AddCovariance(seed);
2690 AliESDtrack *esd=event->GetTrack(i);
2691 if (seed->GetNumberOfClusters()<60 && seed->GetNumberOfClusters()<(esd->GetTPCclusters(0) -5)*0.8){
2692 AliExternalTrackParam paramIn;
2693 AliExternalTrackParam paramOut;
2694 Int_t ncl = seed->RefitTrack(seed,¶mIn,¶mOut);
2695 if (AliTPCReconstructor::StreamLevel()>2) {
2696 (*fDebugStreamer)<<"RecoverBack"<<
2700 "pout.="<<¶mOut<<
2705 seed->Set(paramOut.GetX(),paramOut.GetAlpha(),paramOut.GetParameter(),paramOut.GetCovariance());
2706 seed->SetNumberOfClusters(ncl);
2709 seed->CookdEdx(0.02,0.6);
2710 CookLabel(seed,0.1); //For comparison only
2711 if (seed->GetNumberOfClusters()>15){
2712 esd->UpdateTrackParams(seed,AliESDtrack::kTPCout);
2713 esd->SetTPCPoints(seed->GetPoints());
2714 esd->SetTPCPointsF(seed->GetNFoundable());
2715 Int_t ndedx = seed->GetNCDEDX(0);
2716 Float_t sdedx = seed->GetSDEDX(0);
2717 Float_t dedx = seed->GetdEdx();
2718 esd->SetTPCsignal(dedx, sdedx, ndedx);
2720 Int_t eventnumber = event->GetEventNumberInFile();// patch 28 fev 06
2721 // This is most likely NOT the event number you'd like to use. It has nothing to do with the 'real' event number
2722 if (AliTPCReconstructor::StreamLevel()>1 && esd) {
2723 (*fDebugStreamer)<<"Cback"<<
2726 "EventNrInFile="<<eventnumber<<
2731 //FindKinks(fSeeds,event);
2732 Info("PropagateBack","Number of back propagated tracks %d",ntracks);
2739 void AliTPCtrackerMI::DeleteSeeds()
2744 Int_t nseed = fSeeds->GetEntriesFast();
2745 for (Int_t i=0;i<nseed;i++){
2746 AliTPCseed * seed = (AliTPCseed*)fSeeds->At(i);
2747 if (seed) delete fSeeds->RemoveAt(i);
2754 void AliTPCtrackerMI::ReadSeeds(AliESDEvent *event, Int_t direction)
2757 //read seeds from the event
2759 Int_t nentr=event->GetNumberOfTracks();
2761 Info("ReadSeeds", "Number of ESD tracks: %d\n", nentr);
2766 fSeeds = new TObjArray(nentr);
2770 for (Int_t i=0; i<nentr; i++) {
2771 AliESDtrack *esd=event->GetTrack(i);
2772 ULong_t status=esd->GetStatus();
2773 if (!(status&AliESDtrack::kTPCin)) continue;
2774 AliTPCtrack t(*esd);
2775 t.SetNumberOfClusters(0);
2776 // AliTPCseed *seed = new AliTPCseed(t,t.GetAlpha());
2777 AliTPCseed *seed = new AliTPCseed(t/*,t.GetAlpha()*/);
2778 seed->SetUniqueID(esd->GetID());
2779 AddCovariance(seed); //add systematic ucertainty
2780 for (Int_t ikink=0;ikink<3;ikink++) {
2781 Int_t index = esd->GetKinkIndex(ikink);
2782 seed->GetKinkIndexes()[ikink] = index;
2783 if (index==0) continue;
2784 index = TMath::Abs(index);
2785 AliESDkink * kink = fEvent->GetKink(index-1);
2786 if (kink&&esd->GetKinkIndex(ikink)<0){
2787 if ((status & AliESDtrack::kTRDrefit) != 0) kink->SetStatus(1,2);
2788 if ((status & AliESDtrack::kITSout) != 0) kink->SetStatus(1,0);
2790 if (kink&&esd->GetKinkIndex(ikink)>0){
2791 if ((status & AliESDtrack::kTRDrefit) != 0) kink->SetStatus(1,6);
2792 if ((status & AliESDtrack::kITSout) != 0) kink->SetStatus(1,4);
2796 if (((status&AliESDtrack::kITSout)==0)&&(direction==1)) seed->ResetCovariance(10.);
2797 if ( direction ==2 &&(status & AliESDtrack::kTRDrefit) == 0 ) seed->ResetCovariance(10.);
2798 //if ( direction ==2 && ((status & AliESDtrack::kTPCout) == 0) ) {
2799 // fSeeds->AddAt(0,i);
2803 if ( direction ==2 &&(status & AliESDtrack::kTRDrefit) > 0 ) {
2804 Double_t par0[5],par1[5],alpha,x;
2805 esd->GetInnerExternalParameters(alpha,x,par0);
2806 esd->GetExternalParameters(x,par1);
2807 Double_t delta1 = TMath::Abs(par0[4]-par1[4])/(0.000000001+TMath::Abs(par0[4]+par1[4]));
2808 Double_t delta2 = TMath::Abs(par0[3]-par1[3]);
2810 if (esd->GetTRDncls()>0) trdchi2 = esd->GetTRDchi2()/esd->GetTRDncls();
2811 //reset covariance if suspicious
2812 if ( (delta1>0.1) || (delta2>0.006) ||trdchi2>7.)
2813 seed->ResetCovariance(10.);
2818 // rotate to the local coordinate system
2820 fSectors=fInnerSec; fN=fkNIS;
2821 Double_t alpha=seed->GetAlpha() - fSectors->GetAlphaShift();
2822 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
2823 if (alpha < 0. ) alpha += 2.*TMath::Pi();
2824 Int_t ns=Int_t(alpha/fSectors->GetAlpha())%fN;
2825 alpha =ns*fSectors->GetAlpha() + fSectors->GetAlphaShift();
2826 if (alpha<-TMath::Pi()) alpha += 2*TMath::Pi();
2827 if (alpha>=TMath::Pi()) alpha -= 2*TMath::Pi();
2828 alpha-=seed->GetAlpha();
2829 if (!seed->Rotate(alpha)) {
2835 if (esd->GetKinkIndex(0)<=0){
2836 for (Int_t irow=0;irow<160;irow++){
2837 Int_t index = seed->GetClusterIndex2(irow);
2840 AliTPCclusterMI * cl = GetClusterMI(index);
2841 seed->SetClusterPointer(irow,cl);
2843 if ((index & 0x8000)==0){
2844 cl->Use(10); // accepted cluster
2846 cl->Use(6); // close cluster not accepted
2849 Info("ReadSeeds","Not found cluster");
2854 fSeeds->AddAt(seed,i);
2860 //_____________________________________________________________________________
2861 void AliTPCtrackerMI::MakeSeeds3(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, Float_t cuts[4],
2862 Float_t deltay, Int_t ddsec) {
2863 //-----------------------------------------------------------------
2864 // This function creates track seeds.
2865 // SEEDING WITH VERTEX CONSTRAIN
2866 //-----------------------------------------------------------------
2867 // cuts[0] - fP4 cut
2868 // cuts[1] - tan(phi) cut
2869 // cuts[2] - zvertex cut
2870 // cuts[3] - fP3 cut
2878 Double_t x[5], c[15];
2879 // Int_t di = i1-i2;
2881 AliTPCseed * seed = new AliTPCseed();
2882 Double_t alpha=fSectors->GetAlpha(), shift=fSectors->GetAlphaShift();
2883 Double_t cs=cos(alpha), sn=sin(alpha);
2885 // Double_t x1 =fOuterSec->GetX(i1);
2886 //Double_t xx2=fOuterSec->GetX(i2);
2888 Double_t x1 =GetXrow(i1);
2889 Double_t xx2=GetXrow(i2);
2891 Double_t x3=GetX(), y3=GetY(), z3=GetZ();
2893 Int_t imiddle = (i2+i1)/2; //middle pad row index
2894 Double_t xm = GetXrow(imiddle); // radius of middle pad-row
2895 const AliTPCtrackerRow& krm=GetRow(sec,imiddle); //middle pad -row
2899 const AliTPCtrackerRow& kr1=GetRow(ns,i1);
2900 Double_t ymax = GetMaxY(i1)-kr1.GetDeadZone()-1.5;
2901 Double_t ymaxm = GetMaxY(imiddle)-kr1.GetDeadZone()-1.5;
2904 // change cut on curvature if it can't reach this layer
2905 // maximal curvature set to reach it
2906 Double_t dvertexmax = TMath::Sqrt((x1-x3)*(x1-x3)+(ymax+5-y3)*(ymax+5-y3));
2907 if (dvertexmax*0.5*cuts[0]>0.85){
2908 cuts[0] = 0.85/(dvertexmax*0.5+1.);
2910 Double_t r2min = 1/(cuts[0]*cuts[0]); //minimal square of radius given by cut
2913 if (deltay>0) ddsec = 0;
2914 // loop over clusters
2915 for (Int_t is=0; is < kr1; is++) {
2917 if (kr1[is]->IsUsed(10)) continue;
2918 Double_t y1=kr1[is]->GetY(), z1=kr1[is]->GetZ();
2919 //if (TMath::Abs(y1)>ymax) continue;
2921 if (deltay>0 && TMath::Abs(ymax-TMath::Abs(y1))> deltay ) continue; // seed only at the edge
2923 // find possible directions
2924 Float_t anglez = (z1-z3)/(x1-x3);
2925 Float_t extraz = z1 - anglez*(x1-xx2); // extrapolated z
2928 //find rotation angles relative to line given by vertex and point 1
2929 Double_t dvertex2 = (x1-x3)*(x1-x3)+(y1-y3)*(y1-y3);
2930 Double_t dvertex = TMath::Sqrt(dvertex2);
2931 Double_t angle13 = TMath::ATan((y1-y3)/(x1-x3));
2932 Double_t cs13 = cos(-angle13), sn13 = sin(-angle13);
2935 // loop over 2 sectors
2941 Double_t dddz1=0; // direction of delta inclination in z axis
2948 for (Int_t dsec = dsec1; dsec<=dsec2;dsec++){
2949 Int_t sec2 = sec + dsec;
2951 // AliTPCtrackerRow& kr2 = fOuterSec[(sec2+fkNOS)%fkNOS][i2];
2952 //AliTPCtrackerRow& kr2m = fOuterSec[(sec2+fkNOS)%fkNOS][imiddle];
2953 AliTPCtrackerRow& kr2 = GetRow((sec2+fkNOS)%fkNOS,i2);
2954 AliTPCtrackerRow& kr2m = GetRow((sec2+fkNOS)%fkNOS,imiddle);
2955 Int_t index1 = TMath::Max(kr2.Find(extraz-0.6-dddz1*TMath::Abs(z1)*0.05)-1,0);
2956 Int_t index2 = TMath::Min(kr2.Find(extraz+0.6+dddz2*TMath::Abs(z1)*0.05)+1,kr2);
2958 // rotation angles to p1-p3
2959 Double_t cs13r = cos(-angle13+dsec*alpha)/dvertex, sn13r = sin(-angle13+dsec*alpha)/dvertex;
2960 Double_t x2, y2, z2;
2962 // Double_t dymax = maxangle*TMath::Abs(x1-xx2);
2965 Double_t dxx0 = (xx2-x3)*cs13r;
2966 Double_t dyy0 = (xx2-x3)*sn13r;
2967 for (Int_t js=index1; js < index2; js++) {
2968 const AliTPCclusterMI *kcl = kr2[js];
2969 if (kcl->IsUsed(10)) continue;
2971 //calcutate parameters
2973 Double_t yy0 = dyy0 +(kcl->GetY()-y3)*cs13r;
2975 if (TMath::Abs(yy0)<0.000001) continue;
2976 Double_t xx0 = dxx0 -(kcl->GetY()-y3)*sn13r;
2977 Double_t y0 = 0.5*(xx0*xx0+yy0*yy0-xx0)/yy0;
2978 Double_t r02 = (0.25+y0*y0)*dvertex2;
2979 //curvature (radius) cut
2980 if (r02<r2min) continue;
2984 Double_t c0 = 1/TMath::Sqrt(r02);
2988 //Double_t dfi0 = 2.*TMath::ASin(dvertex*c0*0.5);
2989 //Double_t dfi1 = 2.*TMath::ASin(TMath::Sqrt(yy0*yy0+(1-xx0)*(1-xx0))*dvertex*c0*0.5);
2990 Double_t dfi0 = 2.*AliTPCFastMath::FastAsin(dvertex*c0*0.5);
2991 Double_t dfi1 = 2.*AliTPCFastMath::FastAsin(TMath::Sqrt(yy0*yy0+(1-xx0)*(1-xx0))*dvertex*c0*0.5);
2994 Double_t z0 = kcl->GetZ();
2995 Double_t zzzz2 = z1-(z1-z3)*dfi1/dfi0;
2996 if (TMath::Abs(zzzz2-z0)>0.5) continue;
2999 Double_t dip = (z1-z0)*c0/dfi1;
3000 Double_t x0 = (0.5*cs13+y0*sn13)*dvertex*c0;
3011 x2= xx2*cs-y2*sn*dsec;
3012 y2=+xx2*sn*dsec+y2*cs;
3022 // do we have cluster at the middle ?
3024 GetProlongation(x1,xm,x,ym,zm);
3026 AliTPCclusterMI * cm=0;
3027 if (TMath::Abs(ym)-ymaxm<0){
3028 cm = krm.FindNearest2(ym,zm,1.0,0.6,dummy);
3029 if ((!cm) || (cm->IsUsed(10))) {
3034 // rotate y1 to system 0
3035 // get state vector in rotated system
3036 Double_t yr1 = (-0.5*sn13+y0*cs13)*dvertex*c0;
3037 Double_t xr2 = x0*cs+yr1*sn*dsec;
3038 Double_t xr[5]={kcl->GetY(),kcl->GetZ(), xr2, dip, c0};
3040 GetProlongation(xx2,xm,xr,ym,zm);
3041 if (TMath::Abs(ym)-ymaxm<0){
3042 cm = kr2m.FindNearest2(ym,zm,1.0,0.6,dummy);
3043 if ((!cm) || (cm->IsUsed(10))) {
3053 dym = ym - cm->GetY();
3054 dzm = zm - cm->GetZ();
3061 Double_t sy1=kr1[is]->GetSigmaY2()*2., sz1=kr1[is]->GetSigmaZ2()*2.;
3062 Double_t sy2=kcl->GetSigmaY2()*2., sz2=kcl->GetSigmaZ2()*2.;
3063 //Double_t sy3=400*3./12., sy=0.1, sz=0.1;
3064 Double_t sy3=25000*x[4]*x[4]+0.1, sy=0.1, sz=0.1;
3065 //Double_t sy3=25000*x[4]*x[4]*60+0.5, sy=0.1, sz=0.1;
3067 Double_t f40=(F1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
3068 Double_t f42=(F1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
3069 Double_t f43=(F1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
3070 Double_t f20=(F2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
3071 Double_t f22=(F2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
3072 Double_t f23=(F2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
3074 Double_t f30=(F3(x1,y1+sy,x2,y2,z1,z2)-x[3])/sy;
3075 Double_t f31=(F3(x1,y1,x2,y2,z1+sz,z2)-x[3])/sz;
3076 Double_t f32=(F3(x1,y1,x2,y2+sy,z1,z2)-x[3])/sy;
3077 Double_t f34=(F3(x1,y1,x2,y2,z1,z2+sz)-x[3])/sz;
3081 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3082 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3083 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3084 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3085 c[13]=f30*sy1*f40+f32*sy2*f42;
3086 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3088 // if (!BuildSeed(kr1[is],kcl,0,x1,x2,x3,x,c)) continue;
3090 UInt_t index=kr1.GetIndex(is);
3091 seed->~AliTPCseed(); // this does not set the pointer to 0...
3092 AliTPCseed *track=new(seed) AliTPCseed(x1, ns*alpha+shift, x, c, index);
3094 track->SetIsSeeding(kTRUE);
3095 track->SetSeed1(i1);
3096 track->SetSeed2(i2);
3097 track->SetSeedType(3);
3101 FollowProlongation(*track, (i1+i2)/2,1);
3102 Int_t foundable,found,shared;
3103 track->GetClusterStatistic((i1+i2)/2,i1, found, foundable, shared, kTRUE);
3104 if ((found<0.55*foundable) || shared>0.5*found || (track->GetSigmaY2()+track->GetSigmaZ2())>0.5){
3106 seed->~AliTPCseed();
3112 FollowProlongation(*track, i2,1);
3116 track->SetBConstrain(1);
3117 // track->fLastPoint = i1+fInnerSec->GetNRows(); // first cluster in track position
3118 track->SetLastPoint(i1); // first cluster in track position
3119 track->SetFirstPoint(track->GetLastPoint());
3121 if (track->GetNumberOfClusters()<(i1-i2)*0.5 ||
3122 track->GetNumberOfClusters() < track->GetNFoundable()*0.6 ||
3123 track->GetNShared()>0.4*track->GetNumberOfClusters() ) {
3125 seed->~AliTPCseed();
3129 // Z VERTEX CONDITION
3130 Double_t zv, bz=GetBz();
3131 if ( !track->GetZAt(0.,bz,zv) ) continue;
3132 if (TMath::Abs(zv-z3)>cuts[2]) {
3133 FollowProlongation(*track, TMath::Max(i2-20,0));
3134 if ( !track->GetZAt(0.,bz,zv) ) continue;
3135 if (TMath::Abs(zv-z3)>cuts[2]){
3136 FollowProlongation(*track, TMath::Max(i2-40,0));
3137 if ( !track->GetZAt(0.,bz,zv) ) continue;
3138 if (TMath::Abs(zv-z3)>cuts[2] &&(track->GetNumberOfClusters() > track->GetNFoundable()*0.7)){
3139 // make seed without constrain
3140 AliTPCseed * track2 = MakeSeed(track,0.2,0.5,1.);
3141 FollowProlongation(*track2, i2,1);
3142 track2->SetBConstrain(kFALSE);
3143 track2->SetSeedType(1);
3144 arr->AddLast(track2);
3146 seed->~AliTPCseed();
3151 seed->~AliTPCseed();
3158 track->SetSeedType(0);
3159 arr->AddLast(track);
3160 seed = new AliTPCseed;
3162 // don't consider other combinations
3163 if (track->GetNumberOfClusters() > track->GetNFoundable()*0.8)
3169 Info("MakeSeeds3","\nSeeding statistic:\t%d\t%d\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin,nout1,nout2);
3175 void AliTPCtrackerMI::MakeSeeds5(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, Float_t cuts[4],
3180 //-----------------------------------------------------------------
3181 // This function creates track seeds.
3182 //-----------------------------------------------------------------
3183 // cuts[0] - fP4 cut
3184 // cuts[1] - tan(phi) cut
3185 // cuts[2] - zvertex cut
3186 // cuts[3] - fP3 cut
3196 Double_t x[5], c[15];
3198 // make temporary seed
3199 AliTPCseed * seed = new AliTPCseed;
3200 Double_t alpha=fOuterSec->GetAlpha(), shift=fOuterSec->GetAlphaShift();
3201 // Double_t cs=cos(alpha), sn=sin(alpha);
3206 Double_t x1 = GetXrow(i1-1);
3207 const AliTPCtrackerRow& kr1=GetRow(sec,i1-1);
3208 Double_t y1max = GetMaxY(i1-1)-kr1.GetDeadZone()-1.5;
3210 Double_t x1p = GetXrow(i1);
3211 const AliTPCtrackerRow& kr1p=GetRow(sec,i1);
3213 Double_t x1m = GetXrow(i1-2);
3214 const AliTPCtrackerRow& kr1m=GetRow(sec,i1-2);
3217 //last 3 padrow for seeding
3218 AliTPCtrackerRow& kr3 = GetRow((sec+fkNOS)%fkNOS,i1-7);
3219 Double_t x3 = GetXrow(i1-7);
3220 // Double_t y3max= GetMaxY(i1-7)-kr3.fDeadZone-1.5;
3222 AliTPCtrackerRow& kr3p = GetRow((sec+fkNOS)%fkNOS,i1-6);
3223 Double_t x3p = GetXrow(i1-6);
3225 AliTPCtrackerRow& kr3m = GetRow((sec+fkNOS)%fkNOS,i1-8);
3226 Double_t x3m = GetXrow(i1-8);
3231 Int_t im = i1-4; //middle pad row index
3232 Double_t xm = GetXrow(im); // radius of middle pad-row
3233 const AliTPCtrackerRow& krm=GetRow(sec,im); //middle pad -row
3234 // Double_t ymmax = GetMaxY(im)-kr1.fDeadZone-1.5;
3237 Double_t deltax = x1-x3;
3238 Double_t dymax = deltax*cuts[1];
3239 Double_t dzmax = deltax*cuts[3];
3241 // loop over clusters
3242 for (Int_t is=0; is < kr1; is++) {
3244 if (kr1[is]->IsUsed(10)) continue;
3245 Double_t y1=kr1[is]->GetY(), z1=kr1[is]->GetZ();
3247 if (deltay>0 && TMath::Abs(y1max-TMath::Abs(y1))> deltay ) continue; // seed only at the edge
3249 Int_t index1 = TMath::Max(kr3.Find(z1-dzmax)-1,0);
3250 Int_t index2 = TMath::Min(kr3.Find(z1+dzmax)+1,kr3);
3256 for (Int_t js=index1; js < index2; js++) {
3257 const AliTPCclusterMI *kcl = kr3[js];
3258 if (kcl->IsUsed(10)) continue;
3260 // apply angular cuts
3261 if (TMath::Abs(y1-y3)>dymax) continue;
3264 if (TMath::Abs(z1-z3)>dzmax) continue;
3266 Double_t angley = (y1-y3)/(x1-x3);
3267 Double_t anglez = (z1-z3)/(x1-x3);
3269 Double_t erry = TMath::Abs(angley)*(x1-x1m)*0.5+0.5;
3270 Double_t errz = TMath::Abs(anglez)*(x1-x1m)*0.5+0.5;
3272 Double_t yyym = angley*(xm-x1)+y1;
3273 Double_t zzzm = anglez*(xm-x1)+z1;
3275 const AliTPCclusterMI *kcm = krm.FindNearest2(yyym,zzzm,erry,errz,index);
3277 if (kcm->IsUsed(10)) continue;
3279 erry = TMath::Abs(angley)*(x1-x1m)*0.4+0.5;
3280 errz = TMath::Abs(anglez)*(x1-x1m)*0.4+0.5;
3287 // look around first
3288 const AliTPCclusterMI *kc1m = kr1m.FindNearest2(angley*(x1m-x1)+y1,
3294 if (kc1m->IsUsed(10)) used++;
3296 const AliTPCclusterMI *kc1p = kr1p.FindNearest2(angley*(x1p-x1)+y1,
3302 if (kc1p->IsUsed(10)) used++;
3304 if (used>1) continue;
3305 if (found<1) continue;
3309 const AliTPCclusterMI *kc3m = kr3m.FindNearest2(angley*(x3m-x3)+y3,
3315 if (kc3m->IsUsed(10)) used++;
3319 const AliTPCclusterMI *kc3p = kr3p.FindNearest2(angley*(x3p-x3)+y3,
3325 if (kc3p->IsUsed(10)) used++;
3329 if (used>1) continue;
3330 if (found<3) continue;
3340 x[4]=F1(x1,y1,x2,y2,x3,y3);
3341 //if (TMath::Abs(x[4]) >= cuts[0]) continue;
3344 x[2]=F2(x1,y1,x2,y2,x3,y3);
3347 x[3]=F3n(x1,y1,x2,y2,z1,z2,x[4]);
3348 //if (TMath::Abs(x[3]) > cuts[3]) continue;
3352 Double_t sy1=0.1, sz1=0.1;
3353 Double_t sy2=0.1, sz2=0.1;
3354 Double_t sy3=0.1, sy=0.1, sz=0.1;
3356 Double_t f40=(F1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
3357 Double_t f42=(F1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
3358 Double_t f43=(F1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
3359 Double_t f20=(F2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
3360 Double_t f22=(F2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
3361 Double_t f23=(F2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
3363 Double_t f30=(F3(x1,y1+sy,x2,y2,z1,z2)-x[3])/sy;
3364 Double_t f31=(F3(x1,y1,x2,y2,z1+sz,z2)-x[3])/sz;
3365 Double_t f32=(F3(x1,y1,x2,y2+sy,z1,z2)-x[3])/sy;
3366 Double_t f34=(F3(x1,y1,x2,y2,z1,z2+sz)-x[3])/sz;
3370 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3371 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3372 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3373 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3374 c[13]=f30*sy1*f40+f32*sy2*f42;
3375 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3377 // if (!BuildSeed(kr1[is],kcl,0,x1,x2,x3,x,c)) continue;
3379 index=kr1.GetIndex(is);
3380 seed->~AliTPCseed();
3381 AliTPCseed *track=new(seed) AliTPCseed(x1, sec*alpha+shift, x, c, index);
3383 track->SetIsSeeding(kTRUE);
3386 FollowProlongation(*track, i1-7,1);
3387 if (track->GetNumberOfClusters() < track->GetNFoundable()*0.75 ||
3388 track->GetNShared()>0.6*track->GetNumberOfClusters() || ( track->GetSigmaY2()+ track->GetSigmaZ2())>0.6){
3390 seed->~AliTPCseed();
3396 FollowProlongation(*track, i2,1);
3397 track->SetBConstrain(0);
3398 track->SetLastPoint(i1+fInnerSec->GetNRows()); // first cluster in track position
3399 track->SetFirstPoint(track->GetLastPoint());
3401 if (track->GetNumberOfClusters()<(i1-i2)*0.5 ||
3402 track->GetNumberOfClusters()<track->GetNFoundable()*0.7 ||
3403 track->GetNShared()>2. || track->GetChi2()/track->GetNumberOfClusters()>6 || ( track->GetSigmaY2()+ track->GetSigmaZ2())>0.5 ) {
3405 seed->~AliTPCseed();
3410 FollowProlongation(*track, TMath::Max(i2-10,0),1);
3411 AliTPCseed * track2 = MakeSeed(track,0.2,0.5,0.9);
3412 FollowProlongation(*track2, i2,1);
3413 track2->SetBConstrain(kFALSE);
3414 track2->SetSeedType(4);
3415 arr->AddLast(track2);
3417 seed->~AliTPCseed();
3421 //arr->AddLast(track);
3422 //seed = new AliTPCseed;
3428 Info("MakeSeeds5","\nSeeding statiistic:\t%d\t%d\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin,nout1,nout2,nout3);
3434 //_____________________________________________________________________________
3435 void AliTPCtrackerMI::MakeSeeds2(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, Float_t */*cuts[4]*/,
3436 Float_t deltay, Bool_t /*bconstrain*/) {
3437 //-----------------------------------------------------------------
3438 // This function creates track seeds - without vertex constraint
3439 //-----------------------------------------------------------------
3440 // cuts[0] - fP4 cut - not applied
3441 // cuts[1] - tan(phi) cut
3442 // cuts[2] - zvertex cut - not applied
3443 // cuts[3] - fP3 cut
3453 Double_t alpha=fOuterSec->GetAlpha(), shift=fOuterSec->GetAlphaShift();
3454 // Double_t cs=cos(alpha), sn=sin(alpha);
3455 Int_t row0 = (i1+i2)/2;
3456 Int_t drow = (i1-i2)/2;
3457 const AliTPCtrackerRow& kr0=fSectors[sec][row0];
3458 AliTPCtrackerRow * kr=0;
3460 AliTPCpolyTrack polytrack;
3461 Int_t nclusters=fSectors[sec][row0];
3462 AliTPCseed * seed = new AliTPCseed;
3467 for (Int_t is=0; is < nclusters; is++) { //LOOP over clusters
3469 Int_t nfoundable =0;
3470 for (Int_t iter =1; iter<2; iter++){ //iterations
3471 const AliTPCtrackerRow& krm=fSectors[sec][row0-iter];
3472 const AliTPCtrackerRow& krp=fSectors[sec][row0+iter];
3473 const AliTPCclusterMI * cl= kr0[is];
3475 if (cl->IsUsed(10)) {
3481 Double_t x = kr0.GetX();
3482 // Initialization of the polytrack
3487 Double_t y0= cl->GetY();
3488 Double_t z0= cl->GetZ();
3492 Double_t ymax = fSectors->GetMaxY(row0)-kr0.GetDeadZone()-1.5;
3493 if (deltay>0 && TMath::Abs(ymax-TMath::Abs(y0))> deltay ) continue; // seed only at the edge
3495 erry = (0.5)*cl->GetSigmaY2()/TMath::Sqrt(cl->GetQ())*6;
3496 errz = (0.5)*cl->GetSigmaZ2()/TMath::Sqrt(cl->GetQ())*6;
3497 polytrack.AddPoint(x,y0,z0,erry, errz);
3500 if (cl->IsUsed(10)) sumused++;
3503 Float_t roady = (5*TMath::Sqrt(cl->GetSigmaY2()+0.2)+1.)*iter;
3504 Float_t roadz = (5*TMath::Sqrt(cl->GetSigmaZ2()+0.2)+1.)*iter;
3507 AliTPCclusterMI * cl1 = krm.FindNearest(y0,z0,roady,roadz);
3508 if (cl1 && TMath::Abs(ymax-TMath::Abs(y0))) {
3509 erry = (0.5)*cl1->GetSigmaY2()/TMath::Sqrt(cl1->GetQ())*3;
3510 errz = (0.5)*cl1->GetSigmaZ2()/TMath::Sqrt(cl1->GetQ())*3;
3511 if (cl1->IsUsed(10)) sumused++;
3512 polytrack.AddPoint(x,cl1->GetY(),cl1->GetZ(),erry,errz);
3516 AliTPCclusterMI * cl2 = krp.FindNearest(y0,z0,roady,roadz);
3518 erry = (0.5)*cl2->GetSigmaY2()/TMath::Sqrt(cl2->GetQ())*3;
3519 errz = (0.5)*cl2->GetSigmaZ2()/TMath::Sqrt(cl2->GetQ())*3;
3520 if (cl2->IsUsed(10)) sumused++;
3521 polytrack.AddPoint(x,cl2->GetY(),cl2->GetZ(),erry,errz);
3524 if (sumused>0) continue;
3526 polytrack.UpdateParameters();
3532 nfoundable = polytrack.GetN();
3533 nfound = nfoundable;
3535 for (Int_t ddrow = iter+1; ddrow<drow;ddrow++){
3536 Float_t maxdist = 0.8*(1.+3./(ddrow));
3537 for (Int_t delta = -1;delta<=1;delta+=2){
3538 Int_t row = row0+ddrow*delta;
3539 kr = &(fSectors[sec][row]);
3540 Double_t xn = kr->GetX();
3541 Double_t ymax1 = fSectors->GetMaxY(row)-kr->GetDeadZone()-1.5;
3542 polytrack.GetFitPoint(xn,yn,zn);
3543 if (TMath::Abs(yn)>ymax1) continue;
3545 AliTPCclusterMI * cln = kr->FindNearest(yn,zn,roady,roadz);
3547 Float_t dist = TMath::Sqrt( (yn-cln->GetY())*(yn-cln->GetY())+(zn-cln->GetZ())*(zn-cln->GetZ()));
3550 erry = (dist+0.3)*cln->GetSigmaY2()/TMath::Sqrt(cln->GetQ())*(1.+1./(ddrow));
3551 errz = (dist+0.3)*cln->GetSigmaZ2()/TMath::Sqrt(cln->GetQ())*(1.+1./(ddrow));
3552 if (cln->IsUsed(10)) {
3553 // printf("used\n");
3561 polytrack.AddPoint(xn,cln->GetY(),cln->GetZ(),erry, errz);
3566 if ( (sumused>3) || (sumused>0.5*nfound) || (nfound<0.6*nfoundable)) break;
3567 polytrack.UpdateParameters();
3570 if ( (sumused>3) || (sumused>0.5*nfound)) {
3571 //printf("sumused %d\n",sumused);
3576 polytrack.GetFitDerivation(kr0.GetX(),dy,dz);
3577 AliTPCpolyTrack track2;
3579 polytrack.Refit(track2,0.5+TMath::Abs(dy)*0.3,0.4+TMath::Abs(dz)*0.3);
3580 if (track2.GetN()<0.5*nfoundable) continue;
3583 if ((nfound>0.6*nfoundable) &&( nfoundable>0.4*(i1-i2))) {
3585 // test seed with and without constrain
3586 for (Int_t constrain=0; constrain<=0;constrain++){
3587 // add polytrack candidate
3589 Double_t x[5], c[15];
3590 Double_t x1,x2,x3,y1,y2,y3,z1,z2,z3;
3591 track2.GetBoundaries(x3,x1);
3593 track2.GetFitPoint(x1,y1,z1);
3594 track2.GetFitPoint(x2,y2,z2);
3595 track2.GetFitPoint(x3,y3,z3);
3597 //is track pointing to the vertex ?
3600 polytrack.GetFitPoint(x0,y0,z0);
3613 x[4]=F1(x1,y1,x2,y2,x3,y3);
3615 // if (TMath::Abs(x[4]) >= cuts[0]) continue; //
3616 x[2]=F2(x1,y1,x2,y2,x3,y3);
3618 //if (TMath::Abs(x[4]*x1-x[2]) >= cuts[1]) continue;
3619 //x[3]=F3(x1,y1,x2,y2,z1,z2);
3620 x[3]=F3n(x1,y1,x3,y3,z1,z3,x[4]);
3621 //if (TMath::Abs(x[3]) > cuts[3]) continue;
3624 Double_t sy =0.1, sz =0.1;
3625 Double_t sy1=0.02, sz1=0.02;
3626 Double_t sy2=0.02, sz2=0.02;
3630 sy3=25000*x[4]*x[4]+0.1, sy=0.1, sz=0.1;
3633 Double_t f40=(F1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
3634 Double_t f42=(F1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
3635 Double_t f43=(F1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
3636 Double_t f20=(F2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
3637 Double_t f22=(F2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
3638 Double_t f23=(F2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
3640 Double_t f30=(F3(x1,y1+sy,x3,y3,z1,z3)-x[3])/sy;
3641 Double_t f31=(F3(x1,y1,x3,y3,z1+sz,z3)-x[3])/sz;
3642 Double_t f32=(F3(x1,y1,x3,y3+sy,z1,z3)-x[3])/sy;
3643 Double_t f34=(F3(x1,y1,x3,y3,z1,z3+sz)-x[3])/sz;
3648 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3649 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3650 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3651 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3652 c[13]=f30*sy1*f40+f32*sy2*f42;
3653 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3655 //Int_t row1 = fSectors->GetRowNumber(x1);
3656 Int_t row1 = GetRowNumber(x1);
3660 seed->~AliTPCseed();
3661 AliTPCseed *track=new(seed) AliTPCseed(x1,sec*alpha+shift,x,c,index);
3662 track->SetIsSeeding(kTRUE);
3663 Int_t rc=FollowProlongation(*track, i2);
3664 if (constrain) track->SetBConstrain(1);
3666 track->SetBConstrain(0);
3667 track->SetLastPoint(row1+fInnerSec->GetNRows()); // first cluster in track position
3668 track->SetFirstPoint(track->GetLastPoint());
3670 if (rc==0 || track->GetNumberOfClusters()<(i1-i2)*0.5 ||
3671 track->GetNumberOfClusters() < track->GetNFoundable()*0.6 ||
3672 track->GetNShared()>0.4*track->GetNumberOfClusters()) {
3675 seed->~AliTPCseed();
3678 arr->AddLast(track);
3679 seed = new AliTPCseed;
3683 } // if accepted seed
3686 Info("MakeSeeds2","\nSeeding statiistic:\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin3);
3692 AliTPCseed *AliTPCtrackerMI::MakeSeed(AliTPCseed *track, Float_t r0, Float_t r1, Float_t r2)
3696 //reseed using track points
3697 Int_t p0 = int(r0*track->GetNumberOfClusters()); // point 0
3698 Int_t p1 = int(r1*track->GetNumberOfClusters());
3699 Int_t p2 = int(r2*track->GetNumberOfClusters()); // last point
3701 Double_t x0[3],x1[3],x2[3];
3702 for (Int_t i=0;i<3;i++){
3708 // find track position at given ratio of the length
3709 Int_t sec0=0, sec1=0, sec2=0;
3712 for (Int_t i=0;i<160;i++){
3713 if (track->GetClusterPointer(i)){
3715 AliTPCTrackerPoint *trpoint =track->GetTrackPoint(i);
3716 if ( (index<p0) || x0[0]<0 ){
3717 if (trpoint->GetX()>1){
3718 clindex = track->GetClusterIndex2(i);
3720 x0[0] = trpoint->GetX();
3721 x0[1] = trpoint->GetY();
3722 x0[2] = trpoint->GetZ();
3723 sec0 = ((clindex&0xff000000)>>24)%18;
3728 if ( (index<p1) &&(trpoint->GetX()>1)){
3729 clindex = track->GetClusterIndex2(i);
3731 x1[0] = trpoint->GetX();
3732 x1[1] = trpoint->GetY();
3733 x1[2] = trpoint->GetZ();
3734 sec1 = ((clindex&0xff000000)>>24)%18;
3737 if ( (index<p2) &&(trpoint->GetX()>1)){
3738 clindex = track->GetClusterIndex2(i);
3740 x2[0] = trpoint->GetX();
3741 x2[1] = trpoint->GetY();
3742 x2[2] = trpoint->GetZ();
3743 sec2 = ((clindex&0xff000000)>>24)%18;
3750 Double_t alpha, cs,sn, xx2,yy2;
3752 alpha = (sec1-sec2)*fSectors->GetAlpha();
3753 cs = TMath::Cos(alpha);
3754 sn = TMath::Sin(alpha);
3755 xx2= x1[0]*cs-x1[1]*sn;
3756 yy2= x1[0]*sn+x1[1]*cs;
3760 alpha = (sec0-sec2)*fSectors->GetAlpha();
3761 cs = TMath::Cos(alpha);
3762 sn = TMath::Sin(alpha);
3763 xx2= x0[0]*cs-x0[1]*sn;
3764 yy2= x0[0]*sn+x0[1]*cs;
3770 Double_t x[5],c[15];
3774 x[4]=F1(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]);
3775 // if (x[4]>1) return 0;
3776 x[2]=F2(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]);
3777 x[3]=F3n(x2[0],x2[1],x0[0],x0[1],x2[2],x0[2],x[4]);
3778 //if (TMath::Abs(x[3]) > 2.2) return 0;
3779 //if (TMath::Abs(x[2]) > 1.99) return 0;
3781 Double_t sy =0.1, sz =0.1;
3783 Double_t sy1=0.02+track->GetSigmaY2(), sz1=0.02+track->GetSigmaZ2();
3784 Double_t sy2=0.01+track->GetSigmaY2(), sz2=0.01+track->GetSigmaZ2();
3785 Double_t sy3=0.01+track->GetSigmaY2();
3787 Double_t f40=(F1(x2[0],x2[1]+sy,x1[0],x1[1],x0[0],x0[1])-x[4])/sy;
3788 Double_t f42=(F1(x2[0],x2[1],x1[0],x1[1]+sy,x0[0],x0[1])-x[4])/sy;
3789 Double_t f43=(F1(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]+sy)-x[4])/sy;
3790 Double_t f20=(F2(x2[0],x2[1]+sy,x1[0],x1[1],x0[0],x0[1])-x[2])/sy;
3791 Double_t f22=(F2(x2[0],x2[1],x1[0],x1[1]+sy,x0[0],x0[1])-x[2])/sy;
3792 Double_t f23=(F2(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]+sy)-x[2])/sy;
3794 Double_t f30=(F3(x2[0],x2[1]+sy,x0[0],x0[1],x2[2],x0[2])-x[3])/sy;
3795 Double_t f31=(F3(x2[0],x2[1],x0[0],x0[1],x2[2]+sz,x0[2])-x[3])/sz;
3796 Double_t f32=(F3(x2[0],x2[1],x0[0],x0[1]+sy,x2[2],x0[2])-x[3])/sy;
3797 Double_t f34=(F3(x2[0],x2[1],x0[0],x0[1],x2[2],x0[2]+sz)-x[3])/sz;
3802 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3803 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3804 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3805 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3806 c[13]=f30*sy1*f40+f32*sy2*f42;
3807 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3809 // Int_t row1 = fSectors->GetRowNumber(x2[0]);
3810 AliTPCseed *seed=new AliTPCseed(x2[0], sec2*fSectors->GetAlpha()+fSectors->GetAlphaShift(), x, c, 0);
3811 // Double_t y0,z0,y1,z1, y2,z2;
3812 //seed->GetProlongation(x0[0],y0,z0);
3813 // seed->GetProlongation(x1[0],y1,z1);
3814 //seed->GetProlongation(x2[0],y2,z2);
3816 seed->SetLastPoint(pp2);
3817 seed->SetFirstPoint(pp2);
3824 AliTPCseed *AliTPCtrackerMI::ReSeed(const AliTPCseed *track, Float_t r0, Float_t r1, Float_t r2)
3828 //reseed using founded clusters
3830 // Find the number of clusters
3831 Int_t nclusters = 0;
3832 for (Int_t irow=0;irow<160;irow++){
3833 if (track->GetClusterIndex(irow)>0) nclusters++;
3837 ipos[0] = TMath::Max(int(r0*nclusters),0); // point 0 cluster
3838 ipos[1] = TMath::Min(int(r1*nclusters),nclusters-1); //
3839 ipos[2] = TMath::Min(int(r2*nclusters),nclusters-1); // last point
3843 Int_t row[3],sec[3]={0,0,0};
3845 // find track row position at given ratio of the length
3847 for (Int_t irow=0;irow<160;irow++){
3848 if (track->GetClusterIndex2(irow)<0) continue;
3850 for (Int_t ipoint=0;ipoint<3;ipoint++){
3851 if (index<=ipos[ipoint]) row[ipoint] = irow;
3855 //Get cluster and sector position
3856 for (Int_t ipoint=0;ipoint<3;ipoint++){
3857 Int_t clindex = track->GetClusterIndex2(row[ipoint]);
3858 AliTPCclusterMI * cl = GetClusterMI(clindex);
3861 // AliTPCclusterMI * cl = GetClusterMI(clindex);
3864 sec[ipoint] = ((clindex&0xff000000)>>24)%18;
3865 xyz[ipoint][0] = GetXrow(row[ipoint]);
3866 xyz[ipoint][1] = cl->GetY();
3867 xyz[ipoint][2] = cl->GetZ();
3871 // Calculate seed state vector and covariance matrix
3873 Double_t alpha, cs,sn, xx2,yy2;
3875 alpha = (sec[1]-sec[2])*fSectors->GetAlpha();
3876 cs = TMath::Cos(alpha);
3877 sn = TMath::Sin(alpha);
3878 xx2= xyz[1][0]*cs-xyz[1][1]*sn;
3879 yy2= xyz[1][0]*sn+xyz[1][1]*cs;
3883 alpha = (sec[0]-sec[2])*fSectors->GetAlpha();
3884 cs = TMath::Cos(alpha);
3885 sn = TMath::Sin(alpha);
3886 xx2= xyz[0][0]*cs-xyz[0][1]*sn;
3887 yy2= xyz[0][0]*sn+xyz[0][1]*cs;
3893 Double_t x[5],c[15];
3897 x[4]=F1(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
3898 x[2]=F2(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
3899 x[3]=F3n(xyz[2][0],xyz[2][1],xyz[0][0],xyz[0][1],xyz[2][2],xyz[0][2],x[4]);
3901 Double_t sy =0.1, sz =0.1;
3903 Double_t sy1=0.2, sz1=0.2;
3904 Double_t sy2=0.2, sz2=0.2;
3907 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;
3908 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;
3909 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;
3910 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;
3911 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;
3912 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;
3914 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;
3915 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;
3916 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;
3917 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;
3922 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3923 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3924 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3925 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3926 c[13]=f30*sy1*f40+f32*sy2*f42;
3927 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3929 // Int_t row1 = fSectors->GetRowNumber(xyz[2][0]);
3930 AliTPCseed *seed=new AliTPCseed(xyz[2][0], sec[2]*fSectors->GetAlpha()+fSectors->GetAlphaShift(), x, c, 0);
3931 seed->SetLastPoint(row[2]);
3932 seed->SetFirstPoint(row[2]);
3937 AliTPCseed *AliTPCtrackerMI::ReSeed(AliTPCseed *track,Int_t r0, Bool_t forward)
3941 //reseed using founded clusters
3944 Int_t row[3]={0,0,0};
3945 Int_t sec[3]={0,0,0};
3947 // forward direction
3949 for (Int_t irow=r0;irow<160;irow++){
3950 if (track->GetClusterIndex(irow)>0){
3955 for (Int_t irow=160;irow>r0;irow--){
3956 if (track->GetClusterIndex(irow)>0){
3961 for (Int_t irow=row[2]-15;irow>row[0];irow--){
3962 if (track->GetClusterIndex(irow)>0){
3970 for (Int_t irow=0;irow<r0;irow++){
3971 if (track->GetClusterIndex(irow)>0){
3976 for (Int_t irow=r0;irow>0;irow--){
3977 if (track->GetClusterIndex(irow)>0){
3982 for (Int_t irow=row[2]-15;irow>row[0];irow--){
3983 if (track->GetClusterIndex(irow)>0){
3990 if ((row[2]-row[0])<20) return 0;
3991 if (row[1]==0) return 0;
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 AliTPCTrackerPoint * point = track->GetTrackPoint(row[ipoint]);
4006 if (point&&ipoint<2){
4008 xyz[ipoint][1] = point->GetY();
4009 xyz[ipoint][2] = point->GetZ();
4012 xyz[ipoint][1] = cl->GetY();
4013 xyz[ipoint][2] = cl->GetZ();
4020 // Calculate seed state vector and covariance matrix
4022 Double_t alpha, cs,sn, xx2,yy2;
4024 alpha = (sec[1]-sec[2])*fSectors->GetAlpha();
4025 cs = TMath::Cos(alpha);
4026 sn = TMath::Sin(alpha);
4027 xx2= xyz[1][0]*cs-xyz[1][1]*sn;
4028 yy2= xyz[1][0]*sn+xyz[1][1]*cs;
4032 alpha = (sec[0]-sec[2])*fSectors->GetAlpha();
4033 cs = TMath::Cos(alpha);
4034 sn = TMath::Sin(alpha);
4035 xx2= xyz[0][0]*cs-xyz[0][1]*sn;
4036 yy2= xyz[0][0]*sn+xyz[0][1]*cs;
4042 Double_t x[5],c[15];
4046 x[4]=F1(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
4047 x[2]=F2(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
4048 x[3]=F3n(xyz[2][0],xyz[2][1],xyz[0][0],xyz[0][1],xyz[2][2],xyz[0][2],x[4]);
4050 Double_t sy =0.1, sz =0.1;
4052 Double_t sy1=0.2, sz1=0.2;
4053 Double_t sy2=0.2, sz2=0.2;
4056 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;
4057 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;
4058 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;
4059 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;
4060 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;
4061 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;
4063 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;
4064 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;
4065 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;
4066 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;
4071 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
4072 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
4073 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
4074 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
4075 c[13]=f30*sy1*f40+f32*sy2*f42;
4076 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
4078 // Int_t row1 = fSectors->GetRowNumber(xyz[2][0]);
4079 AliTPCseed *seed=new AliTPCseed(xyz[2][0], sec[2]*fSectors->GetAlpha()+fSectors->GetAlphaShift(), x, c, 0);
4080 seed->SetLastPoint(row[2]);
4081 seed->SetFirstPoint(row[2]);
4082 for (Int_t i=row[0];i<row[2];i++){
4083 seed->SetClusterIndex(i, track->GetClusterIndex(i));
4091 void AliTPCtrackerMI::FindMultiMC(const TObjArray * array, AliESDEvent */*esd*/, Int_t iter)
4094 // find multi tracks - THIS FUNCTION IS ONLY FOR DEBUG PURPOSES
4096 // Use AliTPCReconstructor::StreamLevel()>2 if you want to tune parameters - cuts
4098 // Two reasons to have multiple find tracks
4099 // 1. Curling tracks can be find more than once
4100 // 2. Splitted tracks
4101 // a.) Multiple seeding to increase tracking efficiency - (~ 100% reached)
4102 // b.) Edge effect on the sector boundaries
4105 // Algorithm done in 2 phases - because of CPU consumption
4106 // it is n^2 algorithm - for lead-lead 20000x20000 combination are investigated
4108 // Algorihm for curling tracks sign:
4109 // 1 phase -makes a very rough fast cuts to minimize combinatorics
4110 // a.) opposite sign
4111 // b.) one of the tracks - not pointing to the primary vertex -
4112 // c.) delta tan(theta)
4114 // 2 phase - calculates DCA between tracks - time consument
4119 // General cuts - for splitted tracks and for curling tracks
4121 const Float_t kMaxdPhi = 0.2; // maximal distance in phi
4123 // Curling tracks cuts
4128 Int_t nentries = array->GetEntriesFast();
4129 AliHelix *helixes = new AliHelix[nentries];
4130 Float_t *xm = new Float_t[nentries];
4131 Float_t *dz0 = new Float_t[nentries];
4132 Float_t *dz1 = new Float_t[nentries];
4138 // Find track COG in x direction - point with best defined parameters
4140 for (Int_t i=0;i<nentries;i++){
4141 AliTPCseed* track = (AliTPCseed*)array->At(i);
4142 if (!track) continue;
4143 track->SetCircular(0);
4144 new (&helixes[i]) AliHelix(*track);
4148 track->GetDZ(GetX(),GetY(),GetZ(),GetBz(),dz);
4151 for (Int_t icl=0; icl<160; icl++){
4152 AliTPCclusterMI * cl = track->GetClusterPointer(icl);
4158 if (ncl>0) xm[i]/=Float_t(ncl);
4161 for (Int_t i0=0;i0<nentries;i0++){
4162 AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
4163 if (!track0) continue;
4164 Float_t xc0 = helixes[i0].GetHelix(6);
4165 Float_t yc0 = helixes[i0].GetHelix(7);
4166 Float_t r0 = helixes[i0].GetHelix(8);
4167 Float_t rc0 = TMath::Sqrt(xc0*xc0+yc0*yc0);
4168 Float_t fi0 = TMath::ATan2(yc0,xc0);
4170 for (Int_t i1=i0+1;i1<nentries;i1++){
4171 AliTPCseed * track1 = (AliTPCseed*)array->At(i1);
4172 if (!track1) continue;
4173 Int_t lab0=track0->GetLabel();
4174 Int_t lab1=track1->GetLabel();
4175 if (TMath::Abs(lab0)!=TMath::Abs(lab1)) continue;
4177 Float_t xc1 = helixes[i1].GetHelix(6);
4178 Float_t yc1 = helixes[i1].GetHelix(7);
4179 Float_t r1 = helixes[i1].GetHelix(8);
4180 Float_t rc1 = TMath::Sqrt(xc1*xc1+yc1*yc1);
4181 Float_t fi1 = TMath::ATan2(yc1,xc1);
4183 Float_t dfi = fi0-fi1;
4186 if (dfi>1.5*TMath::Pi()) dfi-=TMath::Pi(); // take care about edge effect
4187 if (dfi<-1.5*TMath::Pi()) dfi+=TMath::Pi(); //
4188 if (TMath::Abs(dfi)>kMaxdPhi&&helixes[i0].GetHelix(4)*helixes[i1].GetHelix(4)<0){
4190 // if short tracks with undefined sign
4191 fi1 = -TMath::ATan2(yc1,-xc1);
4194 Float_t dtheta = TMath::Abs(track0->GetTgl()-track1->GetTgl())<TMath::Abs(track0->GetTgl()+track1->GetTgl())? track0->GetTgl()-track1->GetTgl():track0->GetTgl()+track1->GetTgl();
4197 // debug stream to tune "fast cuts"
4199 Double_t dist[3]; // distance at X
4200 Double_t mdist[3]={0,0,0}; // mean distance X+-40cm
4201 track0->GetDistance(track1,0.5*(xm[i0]+xm[i1])-40.,dist,AliTracker::GetBz());
4202 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[i]);
4203 track0->GetDistance(track1,0.5*(xm[i0]+xm[i1])+40.,dist,AliTracker::GetBz());
4204 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[i]);
4205 track0->GetDistance(track1,0.5*(xm[i0]+xm[i1]),dist,AliTracker::GetBz());
4206 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[i]);
4207 for (Int_t i=0;i<3;i++) mdist[i]*=0.33333;
4211 for (Int_t icl=0; icl<160; icl++){
4212 AliTPCclusterMI * cl0 = track0->GetClusterPointer(icl);
4213 AliTPCclusterMI * cl1 = track1->GetClusterPointer(icl);
4216 if (cl0==cl1) sums++;
4220 if (AliTPCReconstructor::StreamLevel()>5) {
4221 TTreeSRedirector &cstream = *fDebugStreamer;
4226 "Tr0.="<<track0<< // seed0
4227 "Tr1.="<<track1<< // seed1
4228 "h0.="<<&helixes[i0]<<
4229 "h1.="<<&helixes[i1]<<
4231 "sum="<<sum<< //the sum of rows with cl in both
4232 "sums="<<sums<< //the sum of shared clusters
4233 "xm0="<<xm[i0]<< // the center of track
4234 "xm1="<<xm[i1]<< // the x center of track
4235 // General cut variables
4236 "dfi="<<dfi<< // distance in fi angle
4237 "dtheta="<<dtheta<< // distance int theta angle
4243 "dist0="<<dist[0]<< //distance x
4244 "dist1="<<dist[1]<< //distance y
4245 "dist2="<<dist[2]<< //distance z
4246 "mdist0="<<mdist[0]<< //distance x
4247 "mdist1="<<mdist[1]<< //distance y
4248 "mdist2="<<mdist[2]<< //distance z
4262 if (AliTPCReconstructor::StreamLevel()>1) {
4263 AliInfo("Time for curling tracks removal DEBUGGING MC");
4270 void AliTPCtrackerMI::FindSplitted(TObjArray * array, AliESDEvent */*esd*/, Int_t /*iter*/){
4272 // Find Splitted tracks and remove the one with worst quality
4273 // Corresponding debug streamer to tune selections - "Splitted2"
4275 // 0. Sort tracks according quility
4276 // 1. Propagate the tracks to the reference radius
4277 // 2. Double_t loop to select close tracks (only to speed up process)
4278 // 3. Calculate cluster overlap ratio - and remove the track if bigger than a threshold
4279 // 4. Delete temporary parameters
4281 const Double_t xref=GetXrow(63); // reference radius -IROC/OROC boundary
4283 const Double_t kCutP1=10; // delta Z cut 10 cm
4284 const Double_t kCutP2=0.15; // delta snp(fi) cut 0.15
4285 const Double_t kCutP3=0.15; // delta tgl(theta) cut 0.15
4286 const Double_t kCutAlpha=0.15; // delta alpha cut
4287 Int_t firstpoint = 0;
4288 Int_t lastpoint = 160;
4290 Int_t nentries = array->GetEntriesFast();
4291 AliExternalTrackParam *params = new AliExternalTrackParam[nentries];
4297 //0. Sort tracks according quality
4298 //1. Propagate the ext. param to reference radius
4299 Int_t nseed = array->GetEntriesFast();
4300 Float_t * quality = new Float_t[nseed];
4301 Int_t * indexes = new Int_t[nseed];
4302 for (Int_t i=0; i<nseed; i++) {
4303 AliTPCseed *pt=(AliTPCseed*)array->UncheckedAt(i);
4308 pt->UpdatePoints(); //select first last max dens points
4309 Float_t * points = pt->GetPoints();
4310 if (points[3]<0.8) quality[i] =-1;
4311 quality[i] = (points[2]-points[0])+pt->GetNumberOfClusters();
4312 //prefer high momenta tracks if overlaps
4313 quality[i] *= TMath::Sqrt(TMath::Abs(pt->Pt())+0.5);
4315 AliTracker::PropagateTrackToBxByBz(&(params[i]),xref,pt->GetMass(),5.,kTRUE);
4316 AliTracker::PropagateTrackToBxByBz(&(params[i]),xref,pt->GetMass(),1.,kTRUE);
4318 TMath::Sort(nseed,quality,indexes);
4320 // 3. Loop over pair of tracks
4322 for (Int_t i0=0; i0<nseed; i0++) {
4323 Int_t index0=indexes[i0];
4324 if (!(array->UncheckedAt(index0))) continue;
4325 AliTPCseed *s1 = (AliTPCseed*)array->UncheckedAt(index0);
4326 if (!s1->IsActive()) continue;
4327 AliExternalTrackParam &par0=params[index0];
4328 for (Int_t i1=i0+1; i1<nseed; i1++) {
4329 Int_t index1=indexes[i1];
4330 if (!(array->UncheckedAt(index1))) continue;
4331 AliTPCseed *s2 = (AliTPCseed*)array->UncheckedAt(index1);
4332 if (!s2->IsActive()) continue;
4333 if (s2->GetKinkIndexes()[0]!=0)
4334 if (s2->GetKinkIndexes()[0] == -s1->GetKinkIndexes()[0]) continue;
4335 AliExternalTrackParam &par1=params[index1];
4336 if (TMath::Abs(par0.GetParameter()[3]-par1.GetParameter()[3])>kCutP3) continue;
4337 if (TMath::Abs(par0.GetParameter()[1]-par1.GetParameter()[1])>kCutP1) continue;
4338 if (TMath::Abs(par0.GetParameter()[2]-par1.GetParameter()[2])>kCutP2) continue;
4339 Double_t dAlpha= TMath::Abs(par0.GetAlpha()-par1.GetAlpha());
4340 if (dAlpha>TMath::Pi()) dAlpha-=TMath::Pi();
4341 if (TMath::Abs(dAlpha)>kCutAlpha) continue;
4346 Int_t firstShared=lastpoint, lastShared=firstpoint;
4347 Int_t firstRow=lastpoint, lastRow=firstpoint;
4349 for (Int_t i=firstpoint;i<lastpoint;i++){
4350 if (s1->GetClusterIndex2(i)>0) nall0++;
4351 if (s2->GetClusterIndex2(i)>0) nall1++;
4352 if (s1->GetClusterIndex2(i)>0 && s2->GetClusterIndex2(i)>0) {
4353 if (i<firstRow) firstRow=i;
4354 if (i>lastRow) lastRow=i;
4356 if ( (s1->GetClusterIndex2(i))==(s2->GetClusterIndex2(i)) && s1->GetClusterIndex2(i)>0) {
4357 if (i<firstShared) firstShared=i;
4358 if (i>lastShared) lastShared=i;
4362 Double_t ratio0 = Float_t(sumShared)/Float_t(TMath::Min(nall0+1,nall1+1));
4363 Double_t ratio1 = Float_t(sumShared)/Float_t(TMath::Max(nall0+1,nall1+1));
4365 if( AliTPCReconstructor::StreamLevel()>1){
4366 TTreeSRedirector &cstream = *fDebugStreamer;
4367 Int_t n0=s1->GetNumberOfClusters();
4368 Int_t n1=s2->GetNumberOfClusters();
4369 Int_t n0F=s1->GetNFoundable();
4370 Int_t n1F=s2->GetNFoundable();
4371 Int_t lab0=s1->GetLabel();
4372 Int_t lab1=s2->GetLabel();
4374 cstream<<"Splitted2"<<
4375 "iter="<<fIteration<<
4376 "lab0="<<lab0<< // MC label if exist
4377 "lab1="<<lab1<< // MC label if exist
4380 "ratio0="<<ratio0<< // shared ratio
4381 "ratio1="<<ratio1<< // shared ratio
4382 "p0.="<<&par0<< // track parameters
4384 "s0.="<<s1<< // full seed
4386 "n0="<<n0<< // number of clusters track 0
4387 "n1="<<n1<< // number of clusters track 1
4388 "nall0="<<nall0<< // number of clusters track 0
4389 "nall1="<<nall1<< // number of clusters track 1
4390 "n0F="<<n0F<< // number of findable
4391 "n1F="<<n1F<< // number of findable
4392 "shared="<<sumShared<< // number of shared clusters
4393 "firstS="<<firstShared<< // first and the last shared row
4394 "lastS="<<lastShared<<
4395 "firstRow="<<firstRow<< // first and the last row with cluster
4396 "lastRow="<<lastRow<< //
4400 // remove track with lower quality
4402 if (ratio0>AliTPCReconstructor::GetRecoParam()->GetCutSharedClusters(0) ||
4403 ratio1>AliTPCReconstructor::GetRecoParam()->GetCutSharedClusters(1)){
4407 delete array->RemoveAt(index1);
4412 // 4. Delete temporary array
4419 void AliTPCtrackerMI::FindCurling(const TObjArray * array, AliESDEvent */*esd*/, Int_t iter)
4422 // find Curling tracks
4423 // Use AliTPCReconstructor::StreamLevel()>1 if you want to tune parameters - cuts
4426 // Algorithm done in 2 phases - because of CPU consumption
4427 // it is n^2 algorithm - for lead-lead 20000x20000 combination are investigated
4428 // see detal in MC part what can be used to cut
4432 const Float_t kMaxC = 400; // maximal curvature to of the track
4433 const Float_t kMaxdTheta = 0.15; // maximal distance in theta
4434 const Float_t kMaxdPhi = 0.15; // maximal distance in phi
4435 const Float_t kPtRatio = 0.3; // ratio between pt
4436 const Float_t kMinDCAR = 2.; // distance to the primary vertex in r - see cpipe cut
4439 // Curling tracks cuts
4442 const Float_t kMaxDeltaRMax = 40; // distance in outer radius
4443 const Float_t kMaxDeltaRMin = 5.; // distance in lower radius - see cpipe cut
4444 const Float_t kMinAngle = 2.9; // angle between tracks
4445 const Float_t kMaxDist = 5; // biggest distance
4447 // The cuts can be tuned using the "MC information stored in Multi tree ==> see FindMultiMC
4450 TCut csign("csign","Tr0.fP[4]*Tr1.fP[4]<0"); //opposite sign
4451 TCut cmax("cmax","abs(Tr0.GetC())>1/400");
4452 TCut cda("cda","sqrt(dtheta^2+dfi^2)<0.15");
4453 TCut ccratio("ccratio","abs((Tr0.fP[4]+Tr1.fP[4])/(abs(Tr0.fP[4])+abs(Tr1.fP[4])))<0.3");
4454 TCut cpipe("cpipe", "min(abs(r0-rc0),abs(r1-rc1))>5");
4456 TCut cdrmax("cdrmax","abs(abs(rc0+r0)-abs(rc1+r1))<40")
4457 TCut cdrmin("cdrmin","abs(abs(rc0+r0)-abs(rc1+r1))<10")
4459 Multi->Draw("dfi","iter==0"+csign+cmax+cda+ccratio); ~94% of curling tracks fulfill
4460 Multi->Draw("min(abs(r0-rc0),abs(r1-rc1))","iter==0&&abs(lab1)==abs(lab0)"+csign+cmax+cda+ccratio+cpipe+cdrmin+cdrmax); //80%
4462 Curling2->Draw("dfi","iter==0&&abs(lab0)==abs(lab1)"+csign+cmax+cdtheta+cdfi+ccratio)
4468 Int_t nentries = array->GetEntriesFast();
4469 AliHelix *helixes = new AliHelix[nentries];
4470 for (Int_t i=0;i<nentries;i++){
4471 AliTPCseed* track = (AliTPCseed*)array->At(i);
4472 if (!track) continue;
4473 track->SetCircular(0);
4474 new (&helixes[i]) AliHelix(*track);
4480 Double_t phase[2][2],radius[2];
4485 for (Int_t i0=0;i0<nentries;i0++){
4486 AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
4487 if (!track0) continue;
4488 if (TMath::Abs(track0->GetC())<1/kMaxC) continue;
4489 Float_t xc0 = helixes[i0].GetHelix(6);
4490 Float_t yc0 = helixes[i0].GetHelix(7);
4491 Float_t r0 = helixes[i0].GetHelix(8);
4492 Float_t rc0 = TMath::Sqrt(xc0*xc0+yc0*yc0);
4493 Float_t fi0 = TMath::ATan2(yc0,xc0);
4495 for (Int_t i1=i0+1;i1<nentries;i1++){
4496 AliTPCseed * track1 = (AliTPCseed*)array->At(i1);
4497 if (!track1) continue;
4498 if (TMath::Abs(track1->GetC())<1/kMaxC) continue;
4499 Float_t xc1 = helixes[i1].GetHelix(6);
4500 Float_t yc1 = helixes[i1].GetHelix(7);
4501 Float_t r1 = helixes[i1].GetHelix(8);
4502 Float_t rc1 = TMath::Sqrt(xc1*xc1+yc1*yc1);
4503 Float_t fi1 = TMath::ATan2(yc1,xc1);
4505 Float_t dfi = fi0-fi1;
4508 if (dfi>1.5*TMath::Pi()) dfi-=TMath::Pi(); // take care about edge effect
4509 if (dfi<-1.5*TMath::Pi()) dfi+=TMath::Pi(); //
4510 Float_t dtheta = TMath::Abs(track0->GetTgl()-track1->GetTgl())<TMath::Abs(track0->GetTgl()+track1->GetTgl())? track0->GetTgl()-track1->GetTgl():track0->GetTgl()+track1->GetTgl();
4514 if (track0->GetBConstrain()&&track1->GetBConstrain()) continue; // not constrained
4515 if (track1->GetSigned1Pt()*track0->GetSigned1Pt()>0) continue; // not the same sign
4516 if ( TMath::Abs(track1->GetTgl()+track0->GetTgl())>kMaxdTheta) continue; //distance in the Theta
4517 if ( TMath::Abs(dfi)>kMaxdPhi) continue; //distance in phi
4518 if ( TMath::Sqrt(dfi*dfi+dtheta*dtheta)>kMaxdPhi) continue; //common angular offset
4520 Float_t pt0 = track0->GetSignedPt();
4521 Float_t pt1 = track1->GetSignedPt();
4522 if ((TMath::Abs(pt0+pt1)/(TMath::Abs(pt0)+TMath::Abs(pt1)))>kPtRatio) continue;
4523 if ((iter==1) && TMath::Abs(TMath::Abs(rc0+r0)-TMath::Abs(rc1+r1))>kMaxDeltaRMax) continue;
4524 if ((iter!=1) &&TMath::Abs(TMath::Abs(rc0-r0)-TMath::Abs(rc1-r1))>kMaxDeltaRMin) continue;
4525 if (TMath::Min(TMath::Abs(rc0-r0),TMath::Abs(rc1-r1))<kMinDCAR) continue;
4528 // Now find closest approach
4532 Int_t npoints = helixes[i0].GetRPHIintersections(helixes[i1], phase, radius,10);
4533 if (npoints==0) continue;
4534 helixes[i0].GetClosestPhases(helixes[i1], phase);
4538 Double_t hangles[3];
4539 helixes[i0].Evaluate(phase[0][0],xyz0);
4540 helixes[i1].Evaluate(phase[0][1],xyz1);
4542 helixes[i0].GetAngle(phase[0][0],helixes[i1],phase[0][1],hangles);
4543 Double_t deltah[2],deltabest;
4544 if (TMath::Abs(hangles[2])<kMinAngle) continue;
4548 helixes[i0].ParabolicDCA(helixes[i1],phase[0][0],phase[0][1],radius[0],deltah[0],2);
4550 helixes[i0].ParabolicDCA(helixes[i1],phase[1][0],phase[1][1],radius[1],deltah[1],2);
4551 if (deltah[1]<deltah[0]) ibest=1;
4553 deltabest = TMath::Sqrt(deltah[ibest]);
4554 helixes[i0].Evaluate(phase[ibest][0],xyz0);
4555 helixes[i1].Evaluate(phase[ibest][1],xyz1);
4556 helixes[i0].GetAngle(phase[ibest][0],helixes[i1],phase[ibest][1],hangles);
4557 Double_t radiusbest = TMath::Sqrt(radius[ibest]);
4559 if (deltabest>kMaxDist) continue;
4560 // if (mindcar+mindcaz<40 && (TMath::Abs(hangles[2])<kMinAngle ||deltabest>3)) continue;
4561 Bool_t sign =kFALSE;
4562 if (hangles[2]>kMinAngle) sign =kTRUE;
4565 // circular[i0] = kTRUE;
4566 // circular[i1] = kTRUE;
4567 if (track0->OneOverPt()<track1->OneOverPt()){
4568 track0->SetCircular(track0->GetCircular()+1);
4569 track1->SetCircular(track1->GetCircular()+2);
4572 track1->SetCircular(track1->GetCircular()+1);
4573 track0->SetCircular(track0->GetCircular()+2);
4576 if (AliTPCReconstructor::StreamLevel()>2){
4578 //debug stream to tune "fine" cuts
4579 Int_t lab0=track0->GetLabel();
4580 Int_t lab1=track1->GetLabel();
4581 TTreeSRedirector &cstream = *fDebugStreamer;
4582 cstream<<"Curling2"<<
4598 "npoints="<<npoints<<
4599 "hangles0="<<hangles[0]<<
4600 "hangles1="<<hangles[1]<<
4601 "hangles2="<<hangles[2]<<
4604 "radius="<<radiusbest<<
4605 "deltabest="<<deltabest<<
4606 "phase0="<<phase[ibest][0]<<
4607 "phase1="<<phase[ibest][1]<<
4615 if (AliTPCReconstructor::StreamLevel()>1) {
4616 AliInfo("Time for curling tracks removal");
4625 void AliTPCtrackerMI::FindKinks(TObjArray * array, AliESDEvent *esd)
4632 TObjArray *kinks= new TObjArray(10000);
4633 // TObjArray *v0s= new TObjArray(10000);
4634 Int_t nentries = array->GetEntriesFast();
4635 AliHelix *helixes = new AliHelix[nentries];
4636 Int_t *sign = new Int_t[nentries];
4637 Int_t *nclusters = new Int_t[nentries];
4638 Float_t *alpha = new Float_t[nentries];
4639 AliKink *kink = new AliKink();
4640 Int_t * usage = new Int_t[nentries];
4641 Float_t *zm = new Float_t[nentries];
4642 Float_t *z0 = new Float_t[nentries];
4643 Float_t *fim = new Float_t[nentries];
4644 Float_t *shared = new Float_t[nentries];
4645 Bool_t *circular = new Bool_t[nentries];
4646 Float_t *dca = new Float_t[nentries];
4647 //const AliESDVertex * primvertex = esd->GetVertex();
4649 // nentries = array->GetEntriesFast();
4654 for (Int_t i=0;i<nentries;i++){
4657 AliTPCseed* track = (AliTPCseed*)array->At(i);
4658 if (!track) continue;
4659 track->SetCircular(0);
4661 track->UpdatePoints();
4662 if (( track->GetPoints()[2]- track->GetPoints()[0])>5 && track->GetPoints()[3]>0.8){
4664 nclusters[i]=track->GetNumberOfClusters();
4665 alpha[i] = track->GetAlpha();
4666 new (&helixes[i]) AliHelix(*track);
4668 helixes[i].Evaluate(0,xyz);
4669 sign[i] = (track->GetC()>0) ? -1:1;
4672 if (track->GetProlongation(x,y,z)){
4674 fim[i] = alpha[i]+TMath::ATan2(y,x);
4677 zm[i] = track->GetZ();
4681 circular[i]= kFALSE;
4682 if (track->GetProlongation(0,y,z)) z0[i] = z;
4683 dca[i] = track->GetD(0,0);
4689 Int_t ncandidates =0;
4692 Double_t phase[2][2],radius[2];
4695 // Find circling track
4697 for (Int_t i0=0;i0<nentries;i0++){
4698 AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
4699 if (!track0) continue;
4700 if (track0->GetNumberOfClusters()<40) continue;
4701 if (TMath::Abs(1./track0->GetC())>200) continue;
4702 for (Int_t i1=i0+1;i1<nentries;i1++){
4703 AliTPCseed * track1 = (AliTPCseed*)array->At(i1);
4704 if (!track1) continue;
4705 if (track1->GetNumberOfClusters()<40) continue;
4706 if ( TMath::Abs(track1->GetTgl()+track0->GetTgl())>0.1) continue;
4707 if (track0->GetBConstrain()&&track1->GetBConstrain()) continue;
4708 if (TMath::Abs(1./track1->GetC())>200) continue;
4709 if (track1->GetSigned1Pt()*track0->GetSigned1Pt()>0) continue;
4710 if (track1->GetTgl()*track0->GetTgl()>0) continue;
4711 if (TMath::Max(TMath::Abs(1./track0->GetC()),TMath::Abs(1./track1->GetC()))>190) continue;
4712 if (track0->GetBConstrain()&&track1->OneOverPt()<track0->OneOverPt()) continue; //returning - lower momenta
4713 if (track1->GetBConstrain()&&track0->OneOverPt()<track1->OneOverPt()) continue; //returning - lower momenta
4715 Float_t mindcar = TMath::Min(TMath::Abs(dca[i0]),TMath::Abs(dca[i1]));
4716 if (mindcar<5) continue;
4717 Float_t mindcaz = TMath::Min(TMath::Abs(z0[i0]-GetZ()),TMath::Abs(z0[i1]-GetZ()));
4718 if (mindcaz<5) continue;
4719 if (mindcar+mindcaz<20) continue;
4722 Float_t xc0 = helixes[i0].GetHelix(6);
4723 Float_t yc0 = helixes[i0].GetHelix(7);
4724 Float_t r0 = helixes[i0].GetHelix(8);
4725 Float_t xc1 = helixes[i1].GetHelix(6);
4726 Float_t yc1 = helixes[i1].GetHelix(7);
4727 Float_t r1 = helixes[i1].GetHelix(8);
4729 Float_t rmean = (r0+r1)*0.5;
4730 Float_t delta =TMath::Sqrt((xc1-xc0)*(xc1-xc0)+(yc1-yc0)*(yc1-yc0));
4731 //if (delta>30) continue;
4732 if (delta>rmean*0.25) continue;
4733 if (TMath::Abs(r0-r1)/rmean>0.3) continue;
4735 Int_t npoints = helixes[i0].GetRPHIintersections(helixes[i1], phase, radius,10);
4736 if (npoints==0) continue;
4737 helixes[i0].GetClosestPhases(helixes[i1], phase);
4741 Double_t hangles[3];
4742 helixes[i0].Evaluate(phase[0][0],xyz0);
4743 helixes[i1].Evaluate(phase[0][1],xyz1);
4745 helixes[i0].GetAngle(phase[0][0],helixes[i1],phase[0][1],hangles);
4746 Double_t deltah[2],deltabest;
4747 if (hangles[2]<2.8) continue;
4750 helixes[i0].ParabolicDCA(helixes[i1],phase[0][0],phase[0][1],radius[0],deltah[0],2);
4752 helixes[i0].ParabolicDCA(helixes[i1],phase[1][0],phase[1][1],radius[1],deltah[1],2);
4753 if (deltah[1]<deltah[0]) ibest=1;
4755 deltabest = TMath::Sqrt(deltah[ibest]);
4756 helixes[i0].Evaluate(phase[ibest][0],xyz0);
4757 helixes[i1].Evaluate(phase[ibest][1],xyz1);
4758 helixes[i0].GetAngle(phase[ibest][0],helixes[i1],phase[ibest][1],hangles);
4759 Double_t radiusbest = TMath::Sqrt(radius[ibest]);
4761 if (deltabest>6) continue;
4762 if (mindcar+mindcaz<40 && (hangles[2]<3.12||deltabest>3)) continue;
4763 Bool_t lsign =kFALSE;
4764 if (hangles[2]>3.06) lsign =kTRUE;
4767 circular[i0] = kTRUE;
4768 circular[i1] = kTRUE;
4769 if (track0->OneOverPt()<track1->OneOverPt()){
4770 track0->SetCircular(track0->GetCircular()+1);
4771 track1->SetCircular(track1->GetCircular()+2);
4774 track1->SetCircular(track1->GetCircular()+1);
4775 track0->SetCircular(track0->GetCircular()+2);
4778 if (lsign&&AliTPCReconstructor::StreamLevel()>1){
4780 Int_t lab0=track0->GetLabel();
4781 Int_t lab1=track1->GetLabel();
4782 TTreeSRedirector &cstream = *fDebugStreamer;
4783 cstream<<"Curling"<<
4790 "mindcar="<<mindcar<<
4791 "mindcaz="<<mindcaz<<
4794 "npoints="<<npoints<<
4795 "hangles0="<<hangles[0]<<
4796 "hangles2="<<hangles[2]<<
4801 "radius="<<radiusbest<<
4802 "deltabest="<<deltabest<<
4803 "phase0="<<phase[ibest][0]<<
4804 "phase1="<<phase[ibest][1]<<
4814 for (Int_t i =0;i<nentries;i++){
4815 if (sign[i]==0) continue;
4816 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
4823 Double_t cradius0 = 40*40;
4824 Double_t cradius1 = 270*270;
4827 Double_t cdist3=0.55;
4828 for (Int_t j =i+1;j<nentries;j++){
4830 if (sign[j]*sign[i]<1) continue;
4831 if ( (nclusters[i]+nclusters[j])>200) continue;
4832 if ( (nclusters[i]+nclusters[j])<80) continue;
4833 if ( TMath::Abs(zm[i]-zm[j])>60.) continue;
4834 if ( TMath::Abs(fim[i]-fim[j])>0.6 && TMath::Abs(fim[i]-fim[j])<5.7 ) continue;
4835 //AliTPCseed * track1 = (AliTPCseed*)array->At(j); Double_t phase[2][2],radius[2];
4836 Int_t npoints = helixes[i].GetRPHIintersections(helixes[j], phase, radius,20);
4837 if (npoints<1) continue;
4840 if (radius[0]<cradius0||radius[0]>cradius1) continue;
4843 if ( (radius[0]<cradius0||radius[0]>cradius1) && (radius[1]<cradius0||radius[1]>cradius1) ) continue;
4846 Double_t delta1=10000,delta2=10000;
4847 // cuts on the intersection radius
4848 helixes[i].LinearDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta1);
4849 if (radius[0]<20&&delta1<1) continue; //intersection at vertex
4850 if (radius[0]<10&&delta1<3) continue; //intersection at vertex
4852 helixes[i].LinearDCA(helixes[j],phase[1][0],phase[1][1],radius[1],delta2);
4853 if (radius[1]<20&&delta2<1) continue; //intersection at vertex
4854 if (radius[1]<10&&delta2<3) continue; //intersection at vertex
4857 Double_t distance1 = TMath::Min(delta1,delta2);
4858 if (distance1>cdist1) continue; // cut on DCA linear approximation
4860 npoints = helixes[i].GetRPHIintersections(helixes[j], phase, radius,20);
4861 helixes[i].ParabolicDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta1);
4862 if (radius[0]<20&&delta1<1) continue; //intersection at vertex
4863 if (radius[0]<10&&delta1<3) continue; //intersection at vertex
4866 helixes[i].ParabolicDCA(helixes[j],phase[1][0],phase[1][1],radius[1],delta2);
4867 if (radius[1]<20&&delta2<1) continue; //intersection at vertex
4868 if (radius[1]<10&&delta2<3) continue; //intersection at vertex
4870 distance1 = TMath::Min(delta1,delta2);
4873 rkink = TMath::Sqrt(radius[0]);
4876 rkink = TMath::Sqrt(radius[1]);
4878 if (distance1>cdist2) continue;
4881 AliTPCseed * track1 = (AliTPCseed*)array->At(j);
4884 Int_t row0 = GetRowNumber(rkink);
4885 if (row0<10) continue;
4886 if (row0>150) continue;
4889 Float_t dens00=-1,dens01=-1;
4890 Float_t dens10=-1,dens11=-1;
4892 Int_t found,foundable,ishared;
4893 track0->GetClusterStatistic(0,row0-5, found, foundable,ishared,kFALSE);
4894 if (foundable>5) dens00 = Float_t(found)/Float_t(foundable);
4895 track0->GetClusterStatistic(row0+5,155, found, foundable,ishared,kFALSE);
4896 if (foundable>5) dens01 = Float_t(found)/Float_t(foundable);
4898 track1->GetClusterStatistic(0,row0-5, found, foundable,ishared,kFALSE);
4899 if (foundable>10) dens10 = Float_t(found)/Float_t(foundable);
4900 track1->GetClusterStatistic(row0+5,155, found, foundable,ishared,kFALSE);
4901 if (foundable>10) dens11 = Float_t(found)/Float_t(foundable);
4903 if (dens00<dens10 && dens01<dens11) continue;
4904 if (dens00>dens10 && dens01>dens11) continue;
4905 if (TMath::Max(dens00,dens10)<0.1) continue;
4906 if (TMath::Max(dens01,dens11)<0.3) continue;
4908 if (TMath::Min(dens00,dens10)>0.6) continue;
4909 if (TMath::Min(dens01,dens11)>0.6) continue;
4912 AliTPCseed * ktrack0, *ktrack1;
4921 if (TMath::Abs(ktrack0->GetC())>5) continue; // cut on the curvature for mother particle
4922 AliExternalTrackParam paramm(*ktrack0);
4923 AliExternalTrackParam paramd(*ktrack1);
4924 if (row0>60&&ktrack1->GetReference().GetX()>90.)new (¶md) AliExternalTrackParam(ktrack1->GetReference());
4927 kink->SetMother(paramm);
4928 kink->SetDaughter(paramd);
4931 Float_t x[3] = { kink->GetPosition()[0],kink->GetPosition()[1],kink->GetPosition()[2]};
4933 fkParam->Transform0to1(x,index);
4934 fkParam->Transform1to2(x,index);
4935 row0 = GetRowNumber(x[0]);
4937 if (kink->GetR()<100) continue;
4938 if (kink->GetR()>240) continue;
4939 if (kink->GetPosition()[2]/kink->GetR()>AliTPCReconstructor::GetCtgRange()) continue; //out of fiducial volume
4940 if (kink->GetDistance()>cdist3) continue;
4941 Float_t dird = kink->GetDaughterP()[0]*kink->GetPosition()[0]+kink->GetDaughterP()[1]*kink->GetPosition()[1]; // rough direction estimate
4942 if (dird<0) continue;
4944 Float_t dirm = kink->GetMotherP()[0]*kink->GetPosition()[0]+kink->GetMotherP()[1]*kink->GetPosition()[1]; // rough direction estimate
4945 if (dirm<0) continue;
4946 Float_t mpt = TMath::Sqrt(kink->GetMotherP()[0]*kink->GetMotherP()[0]+kink->GetMotherP()[1]*kink->GetMotherP()[1]);
4947 if (mpt<0.2) continue;
4950 //for high momenta momentum not defined well in first iteration
4951 Double_t qt = TMath::Sin(kink->GetAngle(2))*ktrack1->GetP();
4952 if (qt>0.35) continue;
4955 kink->SetLabel(CookLabel(ktrack0,0.4,0,row0),0);
4956 kink->SetLabel(CookLabel(ktrack1,0.4,row0,160),1);
4958 kink->SetTPCDensity(dens00,0,0);
4959 kink->SetTPCDensity(dens01,0,1);
4960 kink->SetTPCDensity(dens10,1,0);
4961 kink->SetTPCDensity(dens11,1,1);
4962 kink->SetIndex(i,0);
4963 kink->SetIndex(j,1);
4966 kink->SetTPCDensity(dens10,0,0);
4967 kink->SetTPCDensity(dens11,0,1);
4968 kink->SetTPCDensity(dens00,1,0);
4969 kink->SetTPCDensity(dens01,1,1);
4970 kink->SetIndex(j,0);
4971 kink->SetIndex(i,1);
4974 if (mpt<1||kink->GetAngle(2)>0.1){
4975 // angle and densities not defined yet
4976 if (kink->GetTPCDensityFactor()<0.8) continue;
4977 if ((2-kink->GetTPCDensityFactor())*kink->GetDistance() >0.25) continue;
4978 if (kink->GetAngle(2)*ktrack0->GetP()<0.003) continue; //too small angle
4979 if (kink->GetAngle(2)>0.2&&kink->GetTPCDensityFactor()<1.15) continue;
4980 if (kink->GetAngle(2)>0.2&&kink->GetTPCDensity(0,1)>0.05) continue;
4982 Float_t criticalangle = track0->GetSigmaSnp2()+track0->GetSigmaTgl2();
4983 criticalangle+= track1->GetSigmaSnp2()+track1->GetSigmaTgl2();
4984 criticalangle= 3*TMath::Sqrt(criticalangle);
4985 if (criticalangle>0.02) criticalangle=0.02;
4986 if (kink->GetAngle(2)<criticalangle) continue;
4989 Int_t drow = Int_t(2.+0.5/(0.05+kink->GetAngle(2))); // overlap region defined
4990 Float_t shapesum =0;
4992 for ( Int_t row = row0-drow; row<row0+drow;row++){
4993 if (row<0) continue;
4994 if (row>155) continue;
4995 if (ktrack0->GetClusterPointer(row)){
4996 AliTPCTrackerPoint *point =ktrack0->GetTrackPoint(row);
4997 shapesum+=point->GetSigmaY()+point->GetSigmaZ();
5000 if (ktrack1->GetClusterPointer(row)){
5001 AliTPCTrackerPoint *point =ktrack1->GetTrackPoint(row);
5002 shapesum+=point->GetSigmaY()+point->GetSigmaZ();
5007 kink->SetShapeFactor(-1.);
5010 kink->SetShapeFactor(shapesum/sum);
5012 // esd->AddKink(kink);
5014 // kink->SetMother(paramm);
5015 //kink->SetDaughter(paramd);
5017 Double_t chi2P2 = paramm.GetParameter()[2]-paramd.GetParameter()[2];
5019 chi2P2/=paramm.GetCovariance()[5]+paramd.GetCovariance()[5];
5020 Double_t chi2P3 = paramm.GetParameter()[3]-paramd.GetParameter()[3];
5022 chi2P3/=paramm.GetCovariance()[9]+paramd.GetCovariance()[9];
5024 if (AliTPCReconstructor::StreamLevel()>1) {
5025 (*fDebugStreamer)<<"kinkLpt"<<
5033 if ( chi2P2+chi2P3<AliTPCReconstructor::GetRecoParam()->GetKinkAngleCutChi2(0)){
5037 kinks->AddLast(kink);
5043 // sort the kinks according quality - and refit them towards vertex
5045 Int_t nkinks = kinks->GetEntriesFast();
5046 Float_t *quality = new Float_t[nkinks];
5047 Int_t *indexes = new Int_t[nkinks];
5048 AliTPCseed *mothers = new AliTPCseed[nkinks];
5049 AliTPCseed *daughters = new AliTPCseed[nkinks];
5052 for (Int_t i=0;i<nkinks;i++){
5054 AliKink *kinkl = (AliKink*)kinks->At(i);
5056 // refit kinks towards vertex
5058 Int_t index0 = kinkl->GetIndex(0);
5059 Int_t index1 = kinkl->GetIndex(1);
5060 AliTPCseed * ktrack0 = (AliTPCseed*)array->At(index0);
5061 AliTPCseed * ktrack1 = (AliTPCseed*)array->At(index1);
5063 Int_t sumn=ktrack0->GetNumberOfClusters()+ktrack1->GetNumberOfClusters();
5065 // Refit Kink under if too small angle
5067 if (kinkl->GetAngle(2)<0.05){
5068 kinkl->SetTPCRow0(GetRowNumber(kinkl->GetR()));
5069 Int_t row0 = kinkl->GetTPCRow0();
5070 Int_t drow = Int_t(2.+0.5/(0.05+kinkl->GetAngle(2)));
5073 Int_t last = row0-drow;
5074 if (last<40) last=40;
5075 if (last<ktrack0->GetFirstPoint()+25) last = ktrack0->GetFirstPoint()+25;
5076 AliTPCseed* seed0 = ReSeed(ktrack0,last,kFALSE);
5079 Int_t first = row0+drow;
5080 if (first>130) first=130;
5081 if (first>ktrack1->GetLastPoint()-25) first = TMath::Max(ktrack1->GetLastPoint()-25,30);
5082 AliTPCseed* seed1 = ReSeed(ktrack1,first,kTRUE);
5084 if (seed0 && seed1){
5085 kinkl->SetStatus(1,8);
5086 if (RefitKink(*seed0,*seed1,*kinkl)) kinkl->SetStatus(1,9);
5087 row0 = GetRowNumber(kinkl->GetR());
5088 sumn = seed0->GetNumberOfClusters()+seed1->GetNumberOfClusters();
5089 mothers[i] = *seed0;
5090 daughters[i] = *seed1;
5093 delete kinks->RemoveAt(i);
5094 if (seed0) delete seed0;
5095 if (seed1) delete seed1;
5098 if (kinkl->GetDistance()>0.5 || kinkl->GetR()<110 || kinkl->GetR()>240) {
5099 delete kinks->RemoveAt(i);
5100 if (seed0) delete seed0;
5101 if (seed1) delete seed1;
5109 if (kinkl) quality[i] = 160*((0.1+kinkl->GetDistance())*(2.-kinkl->GetTPCDensityFactor()))/(sumn+40.); //the longest -clossest will win
5111 TMath::Sort(nkinks,quality,indexes,kFALSE);
5113 //remove double find kinks
5115 for (Int_t ikink0=1;ikink0<nkinks;ikink0++){
5116 AliKink * kink0 = (AliKink*) kinks->At(indexes[ikink0]);
5117 if (!kink0) continue;
5119 for (Int_t ikink1=0;ikink1<ikink0;ikink1++){
5120 if (!kink0) continue;
5121 AliKink * kink1 = (AliKink*) kinks->At(indexes[ikink1]);
5122 if (!kink1) continue;
5123 // if not close kink continue
5124 if (TMath::Abs(kink1->GetPosition()[2]-kink0->GetPosition()[2])>10) continue;
5125 if (TMath::Abs(kink1->GetPosition()[1]-kink0->GetPosition()[1])>10) continue;
5126 if (TMath::Abs(kink1->GetPosition()[0]-kink0->GetPosition()[0])>10) continue;
5128 AliTPCseed &mother0 = mothers[indexes[ikink0]];
5129 AliTPCseed &daughter0 = daughters[indexes[ikink0]];
5130 AliTPCseed &mother1 = mothers[indexes[ikink1]];
5131 AliTPCseed &daughter1 = daughters[indexes[ikink1]];
5132 Int_t row0 = (kink0->GetTPCRow0()+kink1->GetTPCRow0())/2;
5141 for (Int_t i=0;i<row0;i++){
5142 if (mother0.GetClusterIndex(i)>0 && mother1.GetClusterIndex(i)>0){
5145 if (mother0.GetClusterIndex(i)==mother1.GetClusterIndex(i)){
5152 for (Int_t i=row0;i<158;i++){
5153 if (daughter0.GetClusterIndex(i)>0 && daughter0.GetClusterIndex(i)>0){
5156 if (mother0.GetClusterIndex(i)==mother1.GetClusterIndex(i)){
5162 Float_t ratio = Float_t(same+1)/Float_t(both+1);
5163 Float_t ratiom = Float_t(samem+1)/Float_t(bothm+1);
5164 Float_t ratiod = Float_t(samed+1)/Float_t(bothd+1);
5165 if (ratio>0.3 && ratiom>0.5 &&ratiod>0.5) {
5166 Int_t sum0 = mother0.GetNumberOfClusters()+daughter0.GetNumberOfClusters();
5167 Int_t sum1 = mother1.GetNumberOfClusters()+daughter1.GetNumberOfClusters();
5169 shared[kink0->GetIndex(0)]= kTRUE;
5170 shared[kink0->GetIndex(1)]= kTRUE;
5171 delete kinks->RemoveAt(indexes[ikink0]);
5174 shared[kink1->GetIndex(0)]= kTRUE;
5175 shared[kink1->GetIndex(1)]= kTRUE;
5176 delete kinks->RemoveAt(indexes[ikink1]);
5183 for (Int_t i=0;i<nkinks;i++){
5184 AliKink * kinkl = (AliKink*) kinks->At(indexes[i]);
5185 if (!kinkl) continue;
5186 kinkl->SetTPCRow0(GetRowNumber(kinkl->GetR()));
5187 Int_t index0 = kinkl->GetIndex(0);
5188 Int_t index1 = kinkl->GetIndex(1);
5189 if (circular[index0]||(circular[index1]&&kinkl->GetDistance()>0.2)) continue;
5190 kinkl->SetMultiple(usage[index0],0);
5191 kinkl->SetMultiple(usage[index1],1);
5192 if (kinkl->GetMultiple()[0]+kinkl->GetMultiple()[1]>2) continue;
5193 if (kinkl->GetMultiple()[0]+kinkl->GetMultiple()[1]>0 && quality[indexes[i]]>0.2) continue;
5194 if (kinkl->GetMultiple()[0]+kinkl->GetMultiple()[1]>0 && kinkl->GetDistance()>0.2) continue;
5195 if (circular[index0]||(circular[index1]&&kinkl->GetDistance()>0.1)) continue;
5197 AliTPCseed * ktrack0 = (AliTPCseed*)array->At(index0);
5198 AliTPCseed * ktrack1 = (AliTPCseed*)array->At(index1);
5199 if (!ktrack0 || !ktrack1) continue;
5200 Int_t index = esd->AddKink(kinkl);
5203 if ( ktrack0->GetKinkIndex(0)==0 && ktrack1->GetKinkIndex(0)==0) { //best kink
5204 if (mothers[indexes[i]].GetNumberOfClusters()>20 && daughters[indexes[i]].GetNumberOfClusters()>20 && (mothers[indexes[i]].GetNumberOfClusters()+daughters[indexes[i]].GetNumberOfClusters())>100){
5205 *ktrack0 = mothers[indexes[i]];
5206 *ktrack1 = daughters[indexes[i]];
5210 ktrack0->SetKinkIndex(usage[index0],-(index+1));
5211 ktrack1->SetKinkIndex(usage[index1], (index+1));
5216 // Remove tracks corresponding to shared kink's
5218 for (Int_t i=0;i<nentries;i++){
5219 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5220 if (!track0) continue;
5221 if (track0->GetKinkIndex(0)!=0) continue;
5222 if (shared[i]) delete array->RemoveAt(i);
5227 RemoveUsed2(array,0.5,0.4,30);
5229 for (Int_t i=0;i<nentries;i++){
5230 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5231 if (!track0) continue;
5232 track0->CookdEdx(0.02,0.6);
5236 for (Int_t i=0;i<nentries;i++){
5237 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5238 if (!track0) continue;
5239 if (track0->Pt()<1.4) continue;
5240 //remove double high momenta tracks - overlapped with kink candidates
5243 for (Int_t icl=track0->GetFirstPoint();icl<track0->GetLastPoint(); icl++){
5244 if (track0->GetClusterPointer(icl)!=0){
5246 if (track0->GetClusterPointer(icl)->IsUsed(10)) ishared++;
5249 if (Float_t(ishared+1)/Float_t(all+1)>0.5) {
5250 delete array->RemoveAt(i);
5254 if (track0->GetKinkIndex(0)!=0) continue;
5255 if (track0->GetNumberOfClusters()<80) continue;
5257 AliTPCseed *pmother = new AliTPCseed();
5258 AliTPCseed *pdaughter = new AliTPCseed();
5259 AliKink *pkink = new AliKink;
5261 AliTPCseed & mother = *pmother;
5262 AliTPCseed & daughter = *pdaughter;
5263 AliKink & kinkl = *pkink;
5264 if (CheckKinkPoint(track0,mother,daughter, kinkl)){
5265 if (mother.GetNumberOfClusters()<30||daughter.GetNumberOfClusters()<20) {
5269 continue; //too short tracks
5271 if (mother.Pt()<1.4) {
5277 Int_t row0= kinkl.GetTPCRow0();
5278 if (kinkl.GetDistance()>0.5 || kinkl.GetR()<110. || kinkl.GetR()>240.) {
5285 Int_t index = esd->AddKink(&kinkl);
5286 mother.SetKinkIndex(0,-(index+1));
5287 daughter.SetKinkIndex(0,index+1);
5288 if (mother.GetNumberOfClusters()>50) {
5289 delete array->RemoveAt(i);
5290 array->AddAt(new AliTPCseed(mother),i);
5293 array->AddLast(new AliTPCseed(mother));
5295 array->AddLast(new AliTPCseed(daughter));
5296 for (Int_t icl=0;icl<row0;icl++) {
5297 if (mother.GetClusterPointer(icl)) mother.GetClusterPointer(icl)->Use(20);
5300 for (Int_t icl=row0;icl<158;icl++) {
5301 if (daughter.GetClusterPointer(icl)) daughter.GetClusterPointer(icl)->Use(20);
5310 delete [] daughters;
5332 printf("Ncandidates=\t%d\t%d\t%d\t%d\n",esd->GetNumberOfKinks(),ncandidates,ntracks,nall);
5336 void AliTPCtrackerMI::FindV0s(const TObjArray * array, AliESDEvent *esd)
5342 TObjArray *tpcv0s = new TObjArray(100000);
5343 Int_t nentries = array->GetEntriesFast();
5344 AliHelix *helixes = new AliHelix[nentries];
5345 Int_t *sign = new Int_t[nentries];
5346 Float_t *alpha = new Float_t[nentries];
5347 Float_t *z0 = new Float_t[nentries];
5348 Float_t *dca = new Float_t[nentries];
5349 Float_t *sdcar = new Float_t[nentries];
5350 Float_t *cdcar = new Float_t[nentries];
5351 Float_t *pulldcar = new Float_t[nentries];
5352 Float_t *pulldcaz = new Float_t[nentries];
5353 Float_t *pulldca = new Float_t[nentries];
5354 Bool_t *isPrim = new Bool_t[nentries];
5355 const AliESDVertex * primvertex = esd->GetVertex();
5356 Double_t zvertex = primvertex->GetZv();
5358 // nentries = array->GetEntriesFast();
5360 for (Int_t i=0;i<nentries;i++){
5363 AliTPCseed* track = (AliTPCseed*)array->At(i);
5364 if (!track) continue;
5365 track->GetV0Indexes()[0] = 0; //rest v0 indexes
5366 track->GetV0Indexes()[1] = 0; //rest v0 indexes
5367 track->GetV0Indexes()[2] = 0; //rest v0 indexes
5369 alpha[i] = track->GetAlpha();
5370 new (&helixes[i]) AliHelix(*track);
5372 helixes[i].Evaluate(0,xyz);
5373 sign[i] = (track->GetC()>0) ? -1:1;
5377 if (track->GetProlongation(0,y,z)) z0[i] = z;
5378 dca[i] = track->GetD(0,0);
5380 // dca error parrameterezation + pulls
5382 sdcar[i] = TMath::Sqrt(0.150*0.150+(100*track->GetC())*(100*track->GetC()));
5383 if (TMath::Abs(track->GetTgl())>1) sdcar[i]*=2.5;
5384 cdcar[i] = TMath::Exp((TMath::Abs(track->GetC())-0.0106)*525.3);
5385 pulldcar[i] = (dca[i]-cdcar[i])/sdcar[i];
5386 pulldcaz[i] = (z0[i]-zvertex)/sdcar[i];
5387 pulldca[i] = TMath::Sqrt(pulldcar[i]*pulldcar[i]+pulldcaz[i]*pulldcaz[i]);
5388 if (track->TPCrPID(1)+track->TPCrPID(2)+track->TPCrPID(3)>0.5) {
5389 if (pulldca[i]<3.) isPrim[i]=kTRUE; //pion, muon and Kaon 3 sigma cut
5391 if (track->TPCrPID(4)>0.5) {
5392 if (pulldca[i]<0.5) isPrim[i]=kTRUE; //proton 0.5 sigma cut
5394 if (track->TPCrPID(0)>0.4) {
5395 isPrim[i]=kFALSE; //electron no sigma cut
5402 Int_t ncandidates =0;
5405 Double_t phase[2][2],radius[2];
5411 Float_t fprimvertex[3]={GetX(),GetY(),GetZ()};
5413 Double_t cradius0 = 10*10;
5414 Double_t cradius1 = 200*200;
5417 Double_t cpointAngle = 0.95;
5419 Double_t delta[2]={10000,10000};
5420 for (Int_t i =0;i<nentries;i++){
5421 if (sign[i]==0) continue;
5422 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5423 if (!track0) continue;
5424 if (AliTPCReconstructor::StreamLevel()>1){
5425 TTreeSRedirector &cstream = *fDebugStreamer;
5430 "zvertex="<<zvertex<<
5431 "sdcar0="<<sdcar[i]<<
5432 "cdcar0="<<cdcar[i]<<
5433 "pulldcar0="<<pulldcar[i]<<
5434 "pulldcaz0="<<pulldcaz[i]<<
5435 "pulldca0="<<pulldca[i]<<
5436 "isPrim="<<isPrim[i]<<
5440 if (track0->GetSigned1Pt()<0) continue;
5441 if (track0->GetKinkIndex(0)>0||isPrim[i]) continue; //daughter kink
5443 if (TMath::Abs(helixes[i].GetHelix(4))<0.000000001) continue;
5448 for (Int_t j =0;j<nentries;j++){
5449 AliTPCseed * track1 = (AliTPCseed*)array->At(j);
5450 if (!track1) continue;
5451 if (track1->GetKinkIndex(0)>0 || isPrim[j]) continue; //daughter kink
5452 if (sign[j]*sign[i]>0) continue;
5453 if (TMath::Abs(helixes[j].GetHelix(4))<0.000001) continue;
5454 if (track0->GetCircular()+track1->GetCircular()>1) continue; //circling -returning track
5457 // DCA to prim vertex cut
5463 Int_t npoints = helixes[i].GetRPHIintersections(helixes[j], phase, radius,cdist2);
5464 if (npoints<1) continue;
5468 if (radius[0]<cradius0||radius[0]>cradius1) continue;
5469 helixes[i].LinearDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta[0]);
5470 if (delta[0]>cdist1) continue;
5473 if (TMath::Max(radius[0],radius[1])<cradius0|| TMath::Min(radius[0],radius[1])>cradius1) continue;
5474 helixes[i].LinearDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta[0]);
5475 helixes[i].LinearDCA(helixes[j],phase[1][0],phase[1][1],radius[1],delta[1]);
5476 if (delta[1]<delta[0]) iclosest=1;
5477 if (delta[iclosest]>cdist1) continue;
5479 helixes[i].ParabolicDCA(helixes[j],phase[iclosest][0],phase[iclosest][1],radius[iclosest],delta[iclosest]);
5480 if (radius[iclosest]<cradius0 || radius[iclosest]>cradius1 || delta[iclosest]>cdist1) continue;
5482 Double_t pointAngle = helixes[i].GetPointAngle(helixes[j],phase[iclosest],fprimvertex);
5483 if (pointAngle<cpointAngle) continue;
5485 Bool_t isGamma = kFALSE;
5486 vertex.SetParamP(*track0); //track0 - plus
5487 vertex.SetParamN(*track1); //track1 - minus
5488 vertex.Update(fprimvertex);
5489 if (track0->TPCrPID(0)>0.3&&track1->TPCrPID(0)>0.3&&vertex.GetAnglep()[2]<0.15) isGamma=kTRUE; // gamma conversion candidate
5490 Double_t pointAngle2 = vertex.GetV0CosineOfPointingAngle();
5492 if (vertex.GetV0CosineOfPointingAngle()<cpointAngle && (!isGamma)) continue;// point angle cut
5493 //Bo: if (vertex.GetDist2()>2&&(!isGamma)) continue; // point angle cut
5494 if (vertex.GetDcaV0Daughters()>2&&(!isGamma)) continue;//Bo: // point angle cut
5495 Float_t sigmae = 0.15*0.15;
5496 if (vertex.GetRr()<80)
5497 sigmae += (sdcar[i]*sdcar[i]+sdcar[j]*sdcar[j])*(1.-vertex.GetRr()/80.)*(1.-vertex.GetRr()/80.);
5498 sigmae+= TMath::Sqrt(sigmae);
5499 //Bo: if (vertex.GetDist2()/sigmae>3.&&(!isGamma)) continue;
5500 if (vertex.GetDcaV0Daughters()/sigmae>3.&&(!isGamma)) continue;
5501 Float_t densb0=0,densb1=0,densa0=0,densa1=0;
5502 Int_t row0 = GetRowNumber(vertex.GetRr());
5504 //Bo: if (vertex.GetDist2()>0.2) continue;
5505 if (vertex.GetDcaV0Daughters()>0.2) continue;
5506 densb0 = track0->Density2(0,row0-5);
5507 densb1 = track1->Density2(0,row0-5);
5508 if (densb0>0.3|| densb1>0.3) continue; //clusters before vertex
5509 densa0 = track0->Density2(row0+5,row0+40);
5510 densa1 = track1->Density2(row0+5,row0+40);
5511 if ((densa0<0.4|| densa1<0.4)&&(!isGamma)) continue; //missing clusters after vertex
5514 densa0 = track0->Density2(0,40); //cluster density
5515 densa1 = track1->Density2(0,40); //cluster density
5516 if ((vertex.GetRr()<80&&densa0+densa1<1.)&&(!isGamma)) continue;
5518 //Bo: vertex.SetLab(0,track0->GetLabel());
5519 //Bo: vertex.SetLab(1,track1->GetLabel());
5520 vertex.SetChi2After((densa0+densa1)*0.5);
5521 vertex.SetChi2Before((densb0+densb1)*0.5);
5522 vertex.SetIndex(0,i);
5523 vertex.SetIndex(1,j);
5524 //Bo: vertex.SetStatus(1); // TPC v0 candidate
5525 vertex.SetOnFlyStatus(2);//Bo: // TPC v0 candidate
5526 //Bo: vertex.SetRp(track0->TPCrPIDs());
5527 //Bo: vertex.SetRm(track1->TPCrPIDs());
5528 tpcv0s->AddLast(new AliESDv0(vertex));
5531 Int_t eventNr = esd->GetEventNumberInFile(); // This is most likely NOT the event number you'd like to use. It has nothing to do with the 'real' event number
5532 Double_t radiusm= (delta[0]<delta[1])? TMath::Sqrt(radius[0]):TMath::Sqrt(radius[1]);
5533 Double_t deltam= (delta[0]<delta[1])? TMath::Sqrt(delta[0]):TMath::Sqrt(delta[1]);
5534 if (AliTPCReconstructor::StreamLevel()>1) {
5535 Int_t lab0=track0->GetLabel();
5536 Int_t lab1=track1->GetLabel();
5537 Char_t c0=track0->GetCircular();
5538 Char_t c1=track1->GetCircular();
5539 TTreeSRedirector &cstream = *fDebugStreamer;
5542 "vertex.="<<&vertex<<
5545 "Helix0.="<<&helixes[i]<<
5548 "Helix1.="<<&helixes[j]<<
5549 "pointAngle="<<pointAngle<<
5550 "pointAngle2="<<pointAngle2<<
5555 "zvertex="<<zvertex<<
5558 "npoints="<<npoints<<
5559 "radius0="<<radius[0]<<
5560 "delta0="<<delta[0]<<
5561 "radius1="<<radius[1]<<
5562 "delta1="<<delta[1]<<
5563 "radiusm="<<radiusm<<
5565 "sdcar0="<<sdcar[i]<<
5566 "sdcar1="<<sdcar[j]<<
5567 "cdcar0="<<cdcar[i]<<
5568 "cdcar1="<<cdcar[j]<<
5569 "pulldcar0="<<pulldcar[i]<<
5570 "pulldcar1="<<pulldcar[j]<<
5571 "pulldcaz0="<<pulldcaz[i]<<
5572 "pulldcaz1="<<pulldcaz[j]<<
5573 "pulldca0="<<pulldca[i]<<
5574 "pulldca1="<<pulldca[j]<<
5584 Float_t *quality = new Float_t[ncandidates];
5585 Int_t *indexes = new Int_t[ncandidates];
5587 for (Int_t i=0;i<ncandidates;i++){
5589 AliESDv0 *v0 = (AliESDv0*)tpcv0s->At(i);
5590 quality[i] = 1./(1.00001-v0->GetV0CosineOfPointingAngle()); //base point angle
5591 // quality[i] /= (0.5+v0->GetDist2());
5592 // quality[i] *= v0->GetChi2After(); //density factor
5594 Int_t index0 = v0->GetIndex(0);
5595 Int_t index1 = v0->GetIndex(1);
5596 //Bo: Double_t minpulldca = TMath::Min(2.+pulldca[v0->GetIndex(0)],(2.+pulldca[v0->GetIndex(1)]) ); //pull
5597 Double_t minpulldca = TMath::Min(2.+pulldca[index0],(2.+pulldca[index1]) );//Bo:
5601 AliTPCseed * track0 = (AliTPCseed*)array->At(index0);
5602 AliTPCseed * track1 = (AliTPCseed*)array->At(index1);
5603 if (track0->TPCrPID(0)>0.3&&track1->TPCrPID(0)>0.3&&v0->GetAnglep()[2]<0.15) quality[i]+=1000000; // gamma conversion candidate
5604 if (track0->TPCrPID(4)>0.9||(track1->TPCrPID(4)>0.9&&minpulldca>4)) quality[i]*=10; // lambda candidate candidate
5607 TMath::Sort(ncandidates,quality,indexes,kTRUE);
5610 for (Int_t i=0;i<ncandidates;i++){
5611 AliESDv0 * v0 = (AliESDv0*)tpcv0s->At(indexes[i]);
5613 Int_t index0 = v0->GetIndex(0);
5614 Int_t index1 = v0->GetIndex(1);
5615 AliTPCseed * track0 = (AliTPCseed*)array->At(index0);
5616 AliTPCseed * track1 = (AliTPCseed*)array->At(index1);
5617 if (!track0||!track1) {
5621 Bool_t accept =kTRUE; //default accept
5622 Int_t *v0indexes0 = track0->GetV0Indexes();
5623 Int_t *v0indexes1 = track1->GetV0Indexes();
5625 Int_t order0 = (v0indexes0[0]!=0) ? 1:0;
5626 Int_t order1 = (v0indexes1[0]!=0) ? 1:0;
5627 if (v0indexes0[1]!=0) order0 =2;
5628 if (v0indexes1[1]!=0) order1 =2;
5630 if (v0indexes0[2]!=0) {order0=3; accept=kFALSE;}
5631 if (v0indexes0[2]!=0) {order1=3; accept=kFALSE;}
5633 AliESDv0 * v02 = v0;
5635 //Bo: v0->SetOrder(0,order0);
5636 //Bo: v0->SetOrder(1,order1);
5637 //Bo: v0->SetOrder(1,order0+order1);
5638 v0->SetOnFlyStatus(kTRUE);
5639 Int_t index = esd->AddV0(v0);
5640 v02 = esd->GetV0(index);
5641 v0indexes0[order0]=index;
5642 v0indexes1[order1]=index;
5646 Int_t eventNr = esd->GetEventNumberInFile(); // This is most likely NOT the event number you'd like to use. It has nothing to do with the 'real' event number
5647 if (AliTPCReconstructor::StreamLevel()>1) {
5648 Int_t lab0=track0->GetLabel();
5649 Int_t lab1=track1->GetLabel();
5650 TTreeSRedirector &cstream = *fDebugStreamer;
5659 "dca0="<<dca[index0]<<
5660 "dca1="<<dca[index1]<<
5664 "quality="<<quality[i]<<
5665 "pulldca0="<<pulldca[index0]<<
5666 "pulldca1="<<pulldca[index1]<<
5690 printf("TPC V0 finder : naccepted\t%d\tncandidates\t%d\tntracks\t%d\tnall\t%d\n",naccepted,ncandidates,ntracks,nall);
5694 Int_t AliTPCtrackerMI::RefitKink(AliTPCseed &mother, AliTPCseed &daughter, AliESDkink &knk)
5697 // refit kink towards to the vertex
5700 AliKink &kink=(AliKink &)knk;
5702 Int_t row0 = GetRowNumber(kink.GetR());
5703 FollowProlongation(mother,0);
5704 mother.Reset(kFALSE);
5706 FollowProlongation(daughter,row0);
5707 daughter.Reset(kFALSE);
5708 FollowBackProlongation(daughter,158);
5709 daughter.Reset(kFALSE);
5710 Int_t first = TMath::Max(row0-20,30);
5711 Int_t last = TMath::Min(row0+20,140);
5713 const Int_t kNdiv =5;
5714 AliTPCseed param0[kNdiv]; // parameters along the track
5715 AliTPCseed param1[kNdiv]; // parameters along the track
5716 AliKink kinks[kNdiv]; // corresponding kink parameters
5719 for (Int_t irow=0; irow<kNdiv;irow++){
5720 rows[irow] = first +((last-first)*irow)/(kNdiv-1);
5722 // store parameters along the track
5724 for (Int_t irow=0;irow<kNdiv;irow++){
5725 FollowBackProlongation(mother, rows[irow]);
5726 FollowProlongation(daughter,rows[kNdiv-1-irow]);
5727 param0[irow] = mother;
5728 param1[kNdiv-1-irow] = daughter;
5732 for (Int_t irow=0; irow<kNdiv-1;irow++){
5733 if (param0[irow].GetNumberOfClusters()<kNdiv||param1[irow].GetNumberOfClusters()<kNdiv) continue;
5734 kinks[irow].SetMother(param0[irow]);
5735 kinks[irow].SetDaughter(param1[irow]);
5736 kinks[irow].Update();
5739 // choose kink with best "quality"
5741 Double_t mindist = 10000;
5742 for (Int_t irow=0;irow<kNdiv;irow++){
5743 if (param0[irow].GetNumberOfClusters()<20||param1[irow].GetNumberOfClusters()<20) continue;
5744 if (TMath::Abs(kinks[irow].GetR())>240.) continue;
5745 if (TMath::Abs(kinks[irow].GetR())<100.) continue;
5747 Float_t normdist = TMath::Abs(param0[irow].GetX()-kinks[irow].GetR())*(0.1+kink.GetDistance());
5748 normdist/= (param0[irow].GetNumberOfClusters()+param1[irow].GetNumberOfClusters()+40.);
5749 if (normdist < mindist){
5755 if (index==-1) return 0;
5758 param0[index].Reset(kTRUE);
5759 FollowProlongation(param0[index],0);
5761 mother = param0[index];
5762 daughter = param1[index]; // daughter in vertex
5764 kink.SetMother(mother);
5765 kink.SetDaughter(daughter);
5767 kink.SetTPCRow0(GetRowNumber(kink.GetR()));
5768 kink.SetTPCncls(param0[index].GetNumberOfClusters(),0);
5769 kink.SetTPCncls(param1[index].GetNumberOfClusters(),1);
5770 kink.SetLabel(CookLabel(&mother,0.4, 0,kink.GetTPCRow0()),0);
5771 kink.SetLabel(CookLabel(&daughter,0.4, kink.GetTPCRow0(),160),1);
5772 mother.SetLabel(kink.GetLabel(0));
5773 daughter.SetLabel(kink.GetLabel(1));
5779 void AliTPCtrackerMI::UpdateKinkQualityM(AliTPCseed * seed){
5781 // update Kink quality information for mother after back propagation
5783 if (seed->GetKinkIndex(0)>=0) return;
5784 for (Int_t ikink=0;ikink<3;ikink++){
5785 Int_t index = seed->GetKinkIndex(ikink);
5786 if (index>=0) break;
5787 index = TMath::Abs(index)-1;
5788 AliESDkink * kink = fEvent->GetKink(index);
5789 kink->SetTPCDensity(-1,0,0);
5790 kink->SetTPCDensity(1,0,1);
5792 Int_t row0 = kink->GetTPCRow0() - 2 - Int_t( 0.5/ (0.05+kink->GetAngle(2)));
5793 if (row0<15) row0=15;
5795 Int_t row1 = kink->GetTPCRow0() + 2 + Int_t( 0.5/ (0.05+kink->GetAngle(2)));
5796 if (row1>145) row1=145;
5798 Int_t found,foundable,shared;
5799 seed->GetClusterStatistic(0,row0, found, foundable,shared,kFALSE);
5800 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),0,0);
5801 seed->GetClusterStatistic(row1,155, found, foundable,shared,kFALSE);
5802 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),0,1);
5807 void AliTPCtrackerMI::UpdateKinkQualityD(AliTPCseed * seed){
5809 // update Kink quality information for daughter after refit
5811 if (seed->GetKinkIndex(0)<=0) return;
5812 for (Int_t ikink=0;ikink<3;ikink++){
5813 Int_t index = seed->GetKinkIndex(ikink);
5814 if (index<=0) break;
5815 index = TMath::Abs(index)-1;
5816 AliESDkink * kink = fEvent->GetKink(index);
5817 kink->SetTPCDensity(-1,1,0);
5818 kink->SetTPCDensity(-1,1,1);
5820 Int_t row0 = kink->GetTPCRow0() -2 - Int_t( 0.5/ (0.05+kink->GetAngle(2)));
5821 if (row0<15) row0=15;
5823 Int_t row1 = kink->GetTPCRow0() +2 + Int_t( 0.5/ (0.05+kink->GetAngle(2)));
5824 if (row1>145) row1=145;
5826 Int_t found,foundable,shared;
5827 seed->GetClusterStatistic(0,row0, found, foundable,shared,kFALSE);
5828 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),1,0);
5829 seed->GetClusterStatistic(row1,155, found, foundable,shared,kFALSE);
5830 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),1,1);
5836 Int_t AliTPCtrackerMI::CheckKinkPoint(AliTPCseed*seed,AliTPCseed &mother, AliTPCseed &daughter, AliESDkink &knk)
5839 // check kink point for given track
5840 // if return value=0 kink point not found
5841 // otherwise seed0 correspond to mother particle
5842 // seed1 correspond to daughter particle
5843 // kink parameter of kink point
5844 AliKink &kink=(AliKink &)knk;
5846 Int_t middlerow = (seed->GetFirstPoint()+seed->GetLastPoint())/2;
5847 Int_t first = seed->GetFirstPoint();
5848 Int_t last = seed->GetLastPoint();
5849 if (last-first<20) return 0; // shortest length - 2*30 = 60 pad-rows
5852 AliTPCseed *seed1 = ReSeed(seed,middlerow+20, kTRUE); //middle of chamber
5853 if (!seed1) return 0;
5854 FollowProlongation(*seed1,seed->GetLastPoint()-20);
5855 seed1->Reset(kTRUE);
5856 FollowProlongation(*seed1,158);
5857 seed1->Reset(kTRUE);
5858 last = seed1->GetLastPoint();
5860 AliTPCseed *seed0 = new AliTPCseed(*seed);
5861 seed0->Reset(kFALSE);
5864 AliTPCseed param0[20]; // parameters along the track
5865 AliTPCseed param1[20]; // parameters along the track
5866 AliKink kinks[20]; // corresponding kink parameters
5868 for (Int_t irow=0; irow<20;irow++){
5869 rows[irow] = first +((last-first)*irow)/19;
5871 // store parameters along the track
5873 for (Int_t irow=0;irow<20;irow++){
5874 FollowBackProlongation(*seed0, rows[irow]);
5875 FollowProlongation(*seed1,rows[19-irow]);
5876 param0[irow] = *seed0;
5877 param1[19-irow] = *seed1;
5881 for (Int_t irow=0; irow<19;irow++){
5882 kinks[irow].SetMother(param0[irow]);
5883 kinks[irow].SetDaughter(param1[irow]);
5884 kinks[irow].Update();
5887 // choose kink with biggest change of angle
5889 Double_t maxchange= 0;
5890 for (Int_t irow=1;irow<19;irow++){
5891 if (TMath::Abs(kinks[irow].GetR())>240.) continue;
5892 if (TMath::Abs(kinks[irow].GetR())<110.) continue;
5893 Float_t quality = TMath::Abs(kinks[irow].GetAngle(2))/(3.+TMath::Abs(kinks[irow].GetR()-param0[irow].GetX()));
5894 if ( quality > maxchange){
5895 maxchange = quality;
5902 if (index<0) return 0;
5904 Int_t row0 = GetRowNumber(kinks[index].GetR()); //row 0 estimate
5905 seed0 = new AliTPCseed(param0[index]);
5906 seed1 = new AliTPCseed(param1[index]);
5907 seed0->Reset(kFALSE);
5908 seed1->Reset(kFALSE);
5909 seed0->ResetCovariance(10.);
5910 seed1->ResetCovariance(10.);
5911 FollowProlongation(*seed0,0);
5912 FollowBackProlongation(*seed1,158);
5913 mother = *seed0; // backup mother at position 0
5914 seed0->Reset(kFALSE);
5915 seed1->Reset(kFALSE);
5916 seed0->ResetCovariance(10.);
5917 seed1->ResetCovariance(10.);
5919 first = TMath::Max(row0-20,0);
5920 last = TMath::Min(row0+20,158);
5922 for (Int_t irow=0; irow<20;irow++){
5923 rows[irow] = first +((last-first)*irow)/19;
5925 // store parameters along the track
5927 for (Int_t irow=0;irow<20;irow++){
5928 FollowBackProlongation(*seed0, rows[irow]);
5929 FollowProlongation(*seed1,rows[19-irow]);
5930 param0[irow] = *seed0;
5931 param1[19-irow] = *seed1;
5935 for (Int_t irow=0; irow<19;irow++){
5936 kinks[irow].SetMother(param0[irow]);
5937 kinks[irow].SetDaughter(param1[irow]);
5938 // param0[irow].Dump();
5939 //param1[irow].Dump();
5940 kinks[irow].Update();
5943 // choose kink with biggest change of angle
5946 for (Int_t irow=0;irow<20;irow++){
5947 if (TMath::Abs(kinks[irow].GetR())>250.) continue;
5948 if (TMath::Abs(kinks[irow].GetR())<90.) continue;
5949 Float_t quality = TMath::Abs(kinks[irow].GetAngle(2))/(3.+TMath::Abs(kinks[irow].GetR()-param0[irow].GetX()));
5950 if ( quality > maxchange){
5951 maxchange = quality;
5958 if (index==-1 || param0[index].GetNumberOfClusters()+param1[index].GetNumberOfClusters()<100){
5964 // Float_t anglesigma = TMath::Sqrt(param0[index].fC22+param0[index].fC33+param1[index].fC22+param1[index].fC33);
5966 kink.SetMother(param0[index]);
5967 kink.SetDaughter(param1[index]);
5970 Double_t chi2P2 = param0[index].GetParameter()[2]-param1[index].GetParameter()[2];
5972 chi2P2/=param0[index].GetCovariance()[5]+param1[index].GetCovariance()[5];
5973 Double_t chi2P3 = param0[index].GetParameter()[3]-param1[index].GetParameter()[3];
5975 chi2P3/=param0[index].GetCovariance()[9]+param1[index].GetCovariance()[9];
5977 if (AliTPCReconstructor::StreamLevel()>1) {
5978 (*fDebugStreamer)<<"kinkHpt"<<
5981 "p0.="<<¶m0[index]<<
5982 "p1.="<<¶m1[index]<<
5986 if ( chi2P2+chi2P3<AliTPCReconstructor::GetRecoParam()->GetKinkAngleCutChi2(0)){
5993 row0 = GetRowNumber(kink.GetR());
5994 kink.SetTPCRow0(row0);
5995 kink.SetLabel(CookLabel(seed0,0.5,0,row0),0);
5996 kink.SetLabel(CookLabel(seed1,0.5,row0,158),1);
5997 kink.SetIndex(-10,0);
5998 kink.SetIndex(int(param0[index].GetNumberOfClusters()+param1[index].GetNumberOfClusters()),1);
5999 kink.SetTPCncls(param0[index].GetNumberOfClusters(),0);
6000 kink.SetTPCncls(param1[index].GetNumberOfClusters(),1);
6003 // new (&mother) AliTPCseed(param0[index]);
6004 daughter = param1[index];
6005 daughter.SetLabel(kink.GetLabel(1));
6006 param0[index].Reset(kTRUE);
6007 FollowProlongation(param0[index],0);
6008 mother = param0[index];
6009 mother.SetLabel(kink.GetLabel(0));
6010 if ( chi2P2+chi2P3<AliTPCReconstructor::GetRecoParam()->GetKinkAngleCutChi2(1)){
6022 AliTPCseed* AliTPCtrackerMI::ReSeed(AliTPCseed *t)
6025 // reseed - refit - track
6028 // Int_t last = fSectors->GetNRows()-1;
6030 if (fSectors == fOuterSec){
6031 first = TMath::Max(first, t->GetFirstPoint()-fInnerSec->GetNRows());
6035 first = t->GetFirstPoint();
6037 AliTPCseed * seed = MakeSeed(t,0.1,0.5,0.9);
6038 FollowBackProlongation(*t,fSectors->GetNRows()-1);
6040 FollowProlongation(*t,first);
6050 //_____________________________________________________________________________
6051 Int_t AliTPCtrackerMI::ReadSeeds(const TFile *inp) {
6052 //-----------------------------------------------------------------
6053 // This function reades track seeds.
6054 //-----------------------------------------------------------------
6055 TDirectory *savedir=gDirectory;
6057 TFile *in=(TFile*)inp;
6058 if (!in->IsOpen()) {
6059 cerr<<"AliTPCtrackerMI::ReadSeeds(): input file is not open !\n";
6064 TTree *seedTree=(TTree*)in->Get("Seeds");
6066 cerr<<"AliTPCtrackerMI::ReadSeeds(): ";
6067 cerr<<"can't get a tree with track seeds !\n";
6070 AliTPCtrack *seed=new AliTPCtrack;
6071 seedTree->SetBranchAddress("tracks",&seed);
6073 if (fSeeds==0) fSeeds=new TObjArray(15000);
6075 Int_t n=(Int_t)seedTree->GetEntries();
6076 for (Int_t i=0; i<n; i++) {
6077 seedTree->GetEvent(i);
6078 fSeeds->AddLast(new AliTPCseed(*seed/*,seed->GetAlpha()*/));
6087 Int_t AliTPCtrackerMI::Clusters2Tracks (AliESDEvent *esd)
6090 if (fSeeds) DeleteSeeds();
6093 if (!fSeeds) return 1;
6100 //_____________________________________________________________________________
6101 Int_t AliTPCtrackerMI::Clusters2Tracks() {
6102 //-----------------------------------------------------------------
6103 // This is a track finder.
6104 //-----------------------------------------------------------------
6105 TDirectory *savedir=gDirectory;
6109 fSeeds = Tracking();
6112 Info("Clusters2Tracks","Time for tracking: \t");timer.Print();timer.Start();
6114 //activate again some tracks
6115 for (Int_t i=0; i<fSeeds->GetEntriesFast(); i++) {
6116 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
6118 Int_t nc=t.GetNumberOfClusters();
6120 delete fSeeds->RemoveAt(i);
6124 if (pt->GetRemoval()==10) {
6125 if (pt->GetDensityFirst(20)>0.8 || pt->GetDensityFirst(30)>0.8 || pt->GetDensityFirst(40)>0.7)
6126 pt->Desactivate(10); // make track again active
6128 pt->Desactivate(20);
6129 delete fSeeds->RemoveAt(i);
6134 RemoveUsed2(fSeeds,0.85,0.85,0);
6135 if (AliTPCReconstructor::GetRecoParam()->GetDoKinks()) FindKinks(fSeeds,fEvent);
6136 //FindCurling(fSeeds, fEvent,0);
6137 if (AliTPCReconstructor::StreamLevel()>5) FindMultiMC(fSeeds, fEvent,-1); // find multi found tracks
6138 RemoveUsed2(fSeeds,0.5,0.4,20);
6139 FindSplitted(fSeeds, fEvent,0); // find multi found tracks
6140 if (AliTPCReconstructor::StreamLevel()>5) FindMultiMC(fSeeds, fEvent,0); // find multi found tracks
6143 // // refit short tracks
6145 Int_t nseed=fSeeds->GetEntriesFast();
6148 for (Int_t i=0; i<nseed; i++) {
6149 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
6151 Int_t nc=t.GetNumberOfClusters();
6153 delete fSeeds->RemoveAt(i);
6156 CookLabel(pt,0.1); //For comparison only
6157 //if ((pt->IsActive() || (pt->fRemoval==10) )&& nc>50 &&pt->GetNumberOfClusters()>0.4*pt->fNFoundable){
6158 if ((pt->IsActive() || (pt->GetRemoval()==10) )){
6160 if (fDebug>0) cerr<<found<<'\r';
6164 delete fSeeds->RemoveAt(i);
6168 //RemoveOverlap(fSeeds,0.99,7,kTRUE);
6170 //RemoveUsed(fSeeds,0.9,0.9,6);
6172 nseed=fSeeds->GetEntriesFast();
6174 for (Int_t i=0; i<nseed; i++) {
6175 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
6177 Int_t nc=t.GetNumberOfClusters();
6179 delete fSeeds->RemoveAt(i);
6183 t.CookdEdx(0.02,0.6);
6184 // CheckKinkPoint(&t,0.05);
6185 //if ((pt->IsActive() || (pt->fRemoval==10) )&& nc>50 &&pt->GetNumberOfClusters()>0.4*pt->fNFoundable){
6186 if ((pt->IsActive() || (pt->GetRemoval()==10) )){
6194 delete fSeeds->RemoveAt(i);
6195 //AliTPCseed * seed1 = ReSeed(pt,0.05,0.5,1);
6197 // FollowProlongation(*seed1,0);
6198 // Int_t n = seed1->GetNumberOfClusters();
6199 // printf("fP4\t%f\t%f\n",seed1->GetC(),pt->GetC());
6200 // printf("fN\t%d\t%d\n", seed1->GetNumberOfClusters(),pt->GetNumberOfClusters());
6203 //AliTPCseed * seed2 = ReSeed(pt,0.95,0.5,0.05);
6207 SortTracks(fSeeds, 1);
6211 PrepareForBackProlongation(fSeeds,5.);
6212 PropagateBack(fSeeds);
6213 printf("Time for back propagation: \t");timer.Print();timer.Start();
6217 PrepareForProlongation(fSeeds,5.);
6218 PropagateForward2(fSeeds);
6220 printf("Time for FORWARD propagation: \t");timer.Print();timer.Start();
6221 // RemoveUsed(fSeeds,0.7,0.7,6);
6222 //RemoveOverlap(fSeeds,0.9,7,kTRUE);
6224 nseed=fSeeds->GetEntriesFast();
6226 for (Int_t i=0; i<nseed; i++) {
6227 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
6229 Int_t nc=t.GetNumberOfClusters();
6231 delete fSeeds->RemoveAt(i);
6234 t.CookdEdx(0.02,0.6);
6235 // CookLabel(pt,0.1); //For comparison only
6236 //if ((pt->IsActive() || (pt->fRemoval==10) )&& nc>50 &&pt->GetNumberOfClusters()>0.4*pt->fNFoundable){
6237 if ((pt->IsActive() || (pt->fRemoval==10) )){
6238 cerr<<found++<<'\r';
6241 delete fSeeds->RemoveAt(i);
6246 // fNTracks = found;
6248 Info("Clusters2Tracks","Time for overlap removal, track writing and dedx cooking: \t"); timer.Print();timer.Start();
6251 // cerr<<"Number of found tracks : "<<"\t"<<found<<endl;
6252 Info("Clusters2Tracks","Number of found tracks %d",found);
6254 // UnloadClusters();
6259 void AliTPCtrackerMI::Tracking(TObjArray * arr)
6262 // tracking of the seeds
6265 fSectors = fOuterSec;
6266 ParallelTracking(arr,150,63);
6267 fSectors = fOuterSec;
6268 ParallelTracking(arr,63,0);
6271 TObjArray * AliTPCtrackerMI::Tracking(Int_t seedtype, Int_t i1, Int_t i2, Float_t cuts[4], Float_t dy, Int_t dsec)
6276 TObjArray * arr = new TObjArray;
6278 fSectors = fOuterSec;
6281 for (Int_t sec=0;sec<fkNOS;sec++){
6282 if (seedtype==3) MakeSeeds3(arr,sec,i1,i2,cuts,dy, dsec);
6283 if (seedtype==4) MakeSeeds5(arr,sec,i1,i2,cuts,dy);
6284 if (seedtype==2) MakeSeeds2(arr,sec,i1,i2,cuts,dy);
6287 Info("Tracking","\nSeeding - %d\t%d\t%d\t%d\n",seedtype,i1,i2,arr->GetEntriesFast());
6299 TObjArray * AliTPCtrackerMI::Tracking()
6303 if (AliTPCReconstructor::GetRecoParam()->GetSpecialSeeding()) return TrackingSpecial();
6306 Int_t nup=fOuterSec->GetNRows()+fInnerSec->GetNRows();
6308 TObjArray * seeds = new TObjArray;
6317 Float_t fnumber = 3.0;
6318 Float_t fdensity = 3.0;
6323 for (Int_t delta = 0; delta<18; delta+=6){
6327 arr = Tracking(3,nup-1-delta,nup-1-delta-gap,cuts,-1,1);
6328 SumTracks(seeds,arr);
6329 SignClusters(seeds,fnumber,fdensity);
6331 for (Int_t i=2;i<6;i+=2){
6332 // seed high pt tracks
6335 arr = Tracking(3,nup-i-delta,nup-i-delta-gap,cuts,-1,0);
6336 SumTracks(seeds,arr);
6337 SignClusters(seeds,fnumber,fdensity);
6342 // RemoveUsed(seeds,0.9,0.9,1);
6343 // UnsignClusters();
6344 // SignClusters(seeds,fnumber,fdensity);
6348 for (Int_t delta = 20; delta<120; delta+=10){
6350 // seed high pt tracks
6354 arr = Tracking(3,nup-delta,nup-delta-gap,cuts,-1);
6355 SumTracks(seeds,arr);
6356 SignClusters(seeds,fnumber,fdensity);
6361 arr = Tracking(3,nup-delta-5,nup-delta-5-gap,cuts,-1);
6362 SumTracks(seeds,arr);
6363 SignClusters(seeds,fnumber,fdensity);
6374 Info("Tracking()","\n\nPrimary seeding\t%d\n\n",seeds->GetEntriesFast());
6378 // RemoveUsed(seeds,0.75,0.75,1);
6380 //SignClusters(seeds,fnumber,fdensity);
6389 arr = Tracking(4,nup-1,nup-1-gap,cuts,-1);
6390 SumTracks(seeds,arr);
6391 SignClusters(seeds,fnumber,fdensity);
6393 arr = Tracking(4,nup-2,nup-2-gap,cuts,-1);
6394 SumTracks(seeds,arr);
6395 SignClusters(seeds,fnumber,fdensity);
6397 arr = Tracking(4,nup-3,nup-3-gap,cuts,-1);
6398 SumTracks(seeds,arr);
6399 SignClusters(seeds,fnumber,fdensity);
6403 for (Int_t delta = 3; delta<30; delta+=5){
6409 arr = Tracking(4,nup-1-delta,nup-1-delta-gap,cuts,-1);
6410 SumTracks(seeds,arr);
6411 SignClusters(seeds,fnumber,fdensity);
6413 arr = Tracking(4,nup-3-delta,nup-5-delta-gap,cuts,4);
6414 SumTracks(seeds,arr);
6415 SignClusters(seeds,fnumber,fdensity);
6426 Int_t fLastSeedRowSec=AliTPCReconstructor::GetRecoParam()->GetLastSeedRowSec();
6429 for (Int_t delta = 30; delta<fLastSeedRowSec; delta+=10){
6435 arr = Tracking(4,nup-1-delta,nup-1-delta-gap,cuts,-1);
6436 SumTracks(seeds,arr);
6437 SignClusters(seeds,fnumber,fdensity);
6439 arr = Tracking(4,nup-5-delta,nup-5-delta-gap,cuts,5 );
6440 SumTracks(seeds,arr);
6441 SignClusters(seeds,fnumber,fdensity);
6445 Info("Tracking()","\n\nSecondary seeding\t%d\n\n",seeds->GetEntriesFast());
6456 TObjArray * AliTPCtrackerMI::TrackingSpecial()
6459 // seeding adjusted for laser and cosmic tests - short tracks with big inclination angle
6460 // no primary vertex seeding tried
6464 Int_t nup=fOuterSec->GetNRows()+fInnerSec->GetNRows();
6466 TObjArray * seeds = new TObjArray;
6471 Float_t fnumber = 3.0;
6472 Float_t fdensity = 3.0;
6475 cuts[0] = AliTPCReconstructor::GetRecoParam()->GetMaxC(); // max curvature
6476 cuts[1] = 3.5; // max tan(phi) angle for seeding
6477 cuts[2] = 3.; // not used (cut on z primary vertex)
6478 cuts[3] = 3.5; // max tan(theta) angle for seeding
6480 for (Int_t delta = 0; nup-delta-gap-1>0; delta+=3){
6482 arr = Tracking(4,nup-1-delta,nup-1-delta-gap,cuts,-1);
6483 SumTracks(seeds,arr);
6484 SignClusters(seeds,fnumber,fdensity);
6488 Info("Tracking()","\n\nSecondary seeding\t%d\n\n",seeds->GetEntriesFast());
6499 void AliTPCtrackerMI::SumTracks(TObjArray *arr1,TObjArray *&arr2) const
6502 //sum tracks to common container
6503 //remove suspicious tracks
6504 Int_t nseed = arr2->GetEntriesFast();
6505 for (Int_t i=0;i<nseed;i++){
6506 AliTPCseed *pt=(AliTPCseed*)arr2->UncheckedAt(i);
6509 // remove tracks with too big curvature
6511 if (TMath::Abs(pt->GetC())>AliTPCReconstructor::GetRecoParam()->GetMaxC()){
6512 delete arr2->RemoveAt(i);
6515 // REMOVE VERY SHORT TRACKS
6516 if (pt->GetNumberOfClusters()<20){
6517 delete arr2->RemoveAt(i);
6520 // NORMAL ACTIVE TRACK
6521 if (pt->IsActive()){
6522 arr1->AddLast(arr2->RemoveAt(i));
6525 //remove not usable tracks
6526 if (pt->GetRemoval()!=10){
6527 delete arr2->RemoveAt(i);
6531 // ENABLE ONLY ENOUGH GOOD STOPPED TRACKS
6532 if (pt->GetDensityFirst(20)>0.8 || pt->GetDensityFirst(30)>0.8 || pt->GetDensityFirst(40)>0.7)
6533 arr1->AddLast(arr2->RemoveAt(i));
6535 delete arr2->RemoveAt(i);
6539 delete arr2; arr2 = 0;
6544 void AliTPCtrackerMI::ParallelTracking(TObjArray * arr, Int_t rfirst, Int_t rlast)
6547 // try to track in parralel
6549 Int_t nseed=arr->GetEntriesFast();
6550 //prepare seeds for tracking
6551 for (Int_t i=0; i<nseed; i++) {
6552 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i), &t=*pt;
6554 if (!t.IsActive()) continue;
6555 // follow prolongation to the first layer
6556 if ( (fSectors ==fInnerSec) || (t.GetFirstPoint()-fkParam->GetNRowLow()>rfirst+1) )
6557 FollowProlongation(t, rfirst+1);
6562 for (Int_t nr=rfirst; nr>=rlast; nr--){
6563 if (nr<fInnerSec->GetNRows())
6564 fSectors = fInnerSec;
6566 fSectors = fOuterSec;
6567 // make indexes with the cluster tracks for given
6569 // find nearest cluster
6570 for (Int_t i=0; i<nseed; i++) {
6571 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i), &t=*pt;
6573 if (nr==80) pt->UpdateReference();
6574 if (!pt->IsActive()) continue;
6575 // if ( (fSectors ==fOuterSec) && (pt->fFirstPoint-fkParam->GetNRowLow())<nr) continue;
6576 if (pt->GetRelativeSector()>17) {
6579 UpdateClusters(t,nr);
6581 // prolonagate to the nearest cluster - if founded
6582 for (Int_t i=0; i<nseed; i++) {
6583 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
6585 if (!pt->IsActive()) continue;
6586 // if ((fSectors ==fOuterSec) && (pt->fFirstPoint-fkParam->GetNRowLow())<nr) continue;
6587 if (pt->GetRelativeSector()>17) {
6590 FollowToNextCluster(*pt,nr);
6595 void AliTPCtrackerMI::PrepareForBackProlongation(TObjArray * arr,Float_t fac) const
6599 // if we use TPC track itself we have to "update" covariance
6601 Int_t nseed= arr->GetEntriesFast();
6602 for (Int_t i=0;i<nseed;i++){
6603 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
6607 //rotate to current local system at first accepted point
6608 Int_t index = pt->GetClusterIndex2(pt->GetFirstPoint());
6609 Int_t sec = (index&0xff000000)>>24;
6611 Float_t angle1 = fInnerSec->GetAlpha()*sec+fInnerSec->GetAlphaShift();
6612 if (angle1>TMath::Pi())
6613 angle1-=2.*TMath::Pi();
6614 Float_t angle2 = pt->GetAlpha();
6616 if (TMath::Abs(angle1-angle2)>0.001){
6617 pt->Rotate(angle1-angle2);
6618 //angle2 = pt->GetAlpha();
6619 //pt->fRelativeSector = pt->GetAlpha()/fInnerSec->GetAlpha();
6620 //if (pt->GetAlpha()<0)
6621 // pt->fRelativeSector+=18;
6622 //sec = pt->fRelativeSector;
6631 void AliTPCtrackerMI::PrepareForProlongation(TObjArray * arr, Float_t fac) const
6635 // if we use TPC track itself we have to "update" covariance
6637 Int_t nseed= arr->GetEntriesFast();
6638 for (Int_t i=0;i<nseed;i++){
6639 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
6642 pt->SetFirstPoint(pt->GetLastPoint());
6650 Int_t AliTPCtrackerMI::PropagateBack(TObjArray * arr)
6653 // make back propagation
6655 Int_t nseed= arr->GetEntriesFast();
6656 for (Int_t i=0;i<nseed;i++){
6657 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
6658 if (pt&& pt->GetKinkIndex(0)<=0) {
6659 //AliTPCseed *pt2 = new AliTPCseed(*pt);
6660 fSectors = fInnerSec;
6661 //FollowBackProlongation(*pt,fInnerSec->GetNRows()-1);
6662 //fSectors = fOuterSec;
6663 FollowBackProlongation(*pt,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1);
6664 //if (pt->GetNumberOfClusters()<(pt->fEsd->GetTPCclusters(0)) ){
6665 // Error("PropagateBack","Not prolonged track %d",pt->GetLabel());
6666 // FollowBackProlongation(*pt2,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1);
6669 if (pt&& pt->GetKinkIndex(0)>0) {
6670 AliESDkink * kink = fEvent->GetKink(pt->GetKinkIndex(0)-1);
6671 pt->SetFirstPoint(kink->GetTPCRow0());
6672 fSectors = fInnerSec;
6673 FollowBackProlongation(*pt,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1);
6681 Int_t AliTPCtrackerMI::PropagateForward2(TObjArray * arr)
6684 // make forward propagation
6686 Int_t nseed= arr->GetEntriesFast();
6688 for (Int_t i=0;i<nseed;i++){
6689 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
6691 FollowProlongation(*pt,0);
6700 Int_t AliTPCtrackerMI::PropagateForward()
6703 // propagate track forward
6705 Int_t nseed = fSeeds->GetEntriesFast();
6706 for (Int_t i=0;i<nseed;i++){
6707 AliTPCseed *pt = (AliTPCseed*)fSeeds->UncheckedAt(i);
6709 AliTPCseed &t = *pt;
6710 Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
6711 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
6712 if (alpha < 0. ) alpha += 2.*TMath::Pi();
6713 t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN);
6717 fSectors = fOuterSec;
6718 ParallelTracking(fSeeds,fOuterSec->GetNRows()+fInnerSec->GetNRows()-1,fInnerSec->GetNRows());
6719 fSectors = fInnerSec;
6720 ParallelTracking(fSeeds,fInnerSec->GetNRows()-1,0);
6729 Int_t AliTPCtrackerMI::PropagateBack(AliTPCseed * pt, Int_t row0, Int_t row1)
6732 // make back propagation, in between row0 and row1
6736 fSectors = fInnerSec;
6739 if (row1<fSectors->GetNRows())
6742 r1 = fSectors->GetNRows()-1;
6744 if (row0<fSectors->GetNRows()&& r1>0 )
6745 FollowBackProlongation(*pt,r1);
6746 if (row1<=fSectors->GetNRows())
6749 r1 = row1 - fSectors->GetNRows();
6750 if (r1<=0) return 0;
6751 if (r1>=fOuterSec->GetNRows()) return 0;
6752 fSectors = fOuterSec;
6753 return FollowBackProlongation(*pt,r1);
6761 void AliTPCtrackerMI::GetShape(AliTPCseed * seed, Int_t row)
6765 AliTPCClusterParam * clparam = AliTPCcalibDB::Instance()->GetClusterParam();
6766 Float_t zdrift = TMath::Abs((fkParam->GetZLength(0)-TMath::Abs(seed->GetZ())));
6767 Int_t type = (seed->GetSector() < fkParam->GetNSector()/2) ? 0: (row>126) ? 1:2;
6768 Double_t angulary = seed->GetSnp();
6769 angulary = angulary*angulary/((1.-angulary)*(1.+angulary));
6770 Double_t angularz = seed->GetTgl()*seed->GetTgl()*(1.+angulary);
6772 Double_t sigmay = clparam->GetRMS0(0,type,zdrift,TMath::Sqrt(TMath::Abs(angulary)));
6773 Double_t sigmaz = clparam->GetRMS0(1,type,zdrift,TMath::Sqrt(TMath::Abs(angularz)));
6774 seed->SetCurrentSigmaY2(sigmay*sigmay);
6775 seed->SetCurrentSigmaZ2(sigmaz*sigmaz);
6776 // Float_t sd2 = TMath::Abs((fkParam->GetZLength(0)-TMath::Abs(seed->GetZ())))*fkParam->GetDiffL()*fkParam->GetDiffL();
6777 // // Float_t padlength = fkParam->GetPadPitchLength(seed->fSector);
6778 // Float_t padlength = GetPadPitchLength(row);
6780 // Float_t sresy = (seed->GetSector() < fkParam->GetNSector()/2) ? 0.2 :0.3;
6781 // seed->SetCurrentSigmaY2(sd2+padlength*padlength*angulary/12.+sresy*sresy);
6783 // Float_t sresz = fkParam->GetZSigma();
6784 // seed->SetCurrentSigmaZ2(sd2+padlength*padlength*angularz*angularz*(1+angulary)/12.+sresz*sresz);
6786 Float_t wy = GetSigmaY(seed);
6787 Float_t wz = GetSigmaZ(seed);
6790 if (TMath::Abs(wy/seed->fCurrentSigmaY2-1)>0.0001 || TMath::Abs(wz/seed->fCurrentSigmaZ2-1)>0.0001 ){
6791 printf("problem\n");
6798 //__________________________________________________________________________
6799 void AliTPCtrackerMI::CookLabel(AliKalmanTrack *tk, Float_t wrong) const {
6800 //--------------------------------------------------------------------
6801 //This function "cooks" a track label. If label<0, this track is fake.
6802 //--------------------------------------------------------------------
6803 AliTPCseed * t = dynamic_cast<AliTPCseed*>(tk);
6805 printf("%s:%d wrong type \n",(char*)__FILE__,__LINE__);
6809 Int_t noc=t->GetNumberOfClusters();
6811 //printf("\nnot founded prolongation\n\n\n");
6817 AliTPCclusterMI *clusters[160];
6819 for (Int_t i=0;i<160;i++) {
6826 for (i=0; i<160 && current<noc; i++) {
6828 Int_t index=t->GetClusterIndex2(i);
6829 if (index<=0) continue;
6830 if (index&0x8000) continue;
6832 //clusters[current]=GetClusterMI(index);
6833 if (t->GetClusterPointer(i)){
6834 clusters[current]=t->GetClusterPointer(i);
6840 Int_t lab=123456789;
6841 for (i=0; i<noc; i++) {
6842 AliTPCclusterMI *c=clusters[i];
6844 lab=TMath::Abs(c->GetLabel(0));
6846 for (j=0; j<noc; j++) if (lb[j]==lab || mx[j]==0) break;
6852 for (i=0; i<noc; i++) if (mx[i]>max) {max=mx[i]; lab=lb[i];}
6854 for (i=0; i<noc; i++) {
6855 AliTPCclusterMI *c=clusters[i];
6857 if (TMath::Abs(c->GetLabel(1)) == lab ||
6858 TMath::Abs(c->GetLabel(2)) == lab ) max++;
6861 if ((1.- Float_t(max)/noc) > wrong) lab=-lab;
6864 Int_t tail=Int_t(0.10*noc);
6867 for (i=1; i<=160&&ind<tail; i++) {
6868 // AliTPCclusterMI *c=clusters[noc-i];
6869 AliTPCclusterMI *c=clusters[i];
6871 if (lab == TMath::Abs(c->GetLabel(0)) ||
6872 lab == TMath::Abs(c->GetLabel(1)) ||
6873 lab == TMath::Abs(c->GetLabel(2))) max++;
6876 if (max < Int_t(0.5*tail)) lab=-lab;
6883 //delete[] clusters;
6887 //__________________________________________________________________________
6888 Int_t AliTPCtrackerMI::CookLabel(AliTPCseed *t, Float_t wrong,Int_t first, Int_t last) const {
6889 //--------------------------------------------------------------------
6890 //This function "cooks" a track label. If label<0, this track is fake.
6891 //--------------------------------------------------------------------
6892 Int_t noc=t->GetNumberOfClusters();
6894 //printf("\nnot founded prolongation\n\n\n");
6900 AliTPCclusterMI *clusters[160];
6902 for (Int_t i=0;i<160;i++) {
6909 for (i=0; i<160 && current<noc; i++) {
6910 if (i<first) continue;
6911 if (i>last) continue;
6912 Int_t index=t->GetClusterIndex2(i);
6913 if (index<=0) continue;
6914 if (index&0x8000) continue;
6916 //clusters[current]=GetClusterMI(index);
6917 if (t->GetClusterPointer(i)){
6918 clusters[current]=t->GetClusterPointer(i);
6923 if (noc<5) return -1;
6924 Int_t lab=123456789;
6925 for (i=0; i<noc; i++) {
6926 AliTPCclusterMI *c=clusters[i];
6928 lab=TMath::Abs(c->GetLabel(0));
6930 for (j=0; j<noc; j++) if (lb[j]==lab || mx[j]==0) break;
6936 for (i=0; i<noc; i++) if (mx[i]>max) {max=mx[i]; lab=lb[i];}
6938 for (i=0; i<noc; i++) {
6939 AliTPCclusterMI *c=clusters[i];
6941 if (TMath::Abs(c->GetLabel(1)) == lab ||
6942 TMath::Abs(c->GetLabel(2)) == lab ) max++;
6945 if ((1.- Float_t(max)/noc) > wrong) lab=-lab;
6948 Int_t tail=Int_t(0.10*noc);
6951 for (i=1; i<=160&&ind<tail; i++) {
6952 // AliTPCclusterMI *c=clusters[noc-i];
6953 AliTPCclusterMI *c=clusters[i];
6955 if (lab == TMath::Abs(c->GetLabel(0)) ||
6956 lab == TMath::Abs(c->GetLabel(1)) ||
6957 lab == TMath::Abs(c->GetLabel(2))) max++;
6960 if (max < Int_t(0.5*tail)) lab=-lab;
6963 // t->SetLabel(lab);
6967 //delete[] clusters;
6971 Int_t AliTPCtrackerMI::GetRowNumber(Double_t x[3]) const
6973 //return pad row number for given x vector
6974 Float_t phi = TMath::ATan2(x[1],x[0]);
6975 if(phi<0) phi=2.*TMath::Pi()+phi;
6976 // Get the local angle in the sector philoc
6977 const Float_t kRaddeg = 180/3.14159265358979312;
6978 Float_t phiangle = (Int_t (phi*kRaddeg/20.) + 0.5)*20./kRaddeg;
6979 Double_t localx = x[0]*TMath::Cos(phiangle)-x[1]*TMath::Sin(phiangle);
6980 return GetRowNumber(localx);
6985 void AliTPCtrackerMI::MakeBitmaps(AliTPCseed *t)
6987 //-----------------------------------------------------------------------
6988 // Fill the cluster and sharing bitmaps of the track
6989 //-----------------------------------------------------------------------
6991 Int_t firstpoint = 0;
6992 Int_t lastpoint = 159;
6993 AliTPCTrackerPoint *point;
6994 AliTPCclusterMI *cluster;
6996 for (int iter=firstpoint; iter<lastpoint; iter++) {
6997 // Change to cluster pointers to see if we have a cluster at given padrow
6998 cluster = t->GetClusterPointer(iter);
7000 t->SetClusterMapBit(iter, kTRUE);
7001 point = t->GetTrackPoint(iter);
7002 if (point->IsShared())
7003 t->SetSharedMapBit(iter,kTRUE);
7005 t->SetSharedMapBit(iter, kFALSE);
7008 t->SetClusterMapBit(iter, kFALSE);
7009 t->SetSharedMapBit(iter, kFALSE);
7014 Bool_t AliTPCtrackerMI::IsFindable(AliTPCseed & track){
7016 // return flag if there is findable cluster at given position
7019 Float_t z = track.GetZ();
7021 if (TMath::Abs(z)<(AliTPCReconstructor::GetCtgRange()*track.GetX()+kDeltaZ) &&
7022 TMath::Abs(z)<fkParam->GetZLength(0) &&
7023 (TMath::Abs(track.GetSnp())<AliTPCReconstructor::GetMaxSnpTracker()))
7029 void AliTPCtrackerMI::AddCovariance(AliTPCseed * seed){
7031 // Adding systematic error
7032 // !!!! the systematic error for element 4 is in 1/cm not in pt
7034 const Double_t *param = AliTPCReconstructor::GetRecoParam()->GetSystematicError();
7036 for (Int_t i=0;i<15;i++) covar[i]=0;
7042 covar[0] = param[0]*param[0];
7043 covar[2] = param[1]*param[1];
7044 covar[5] = param[2]*param[2];
7045 covar[9] = param[3]*param[3];
7046 Double_t facC = AliTracker::GetBz()*kB2C;
7047 covar[14]= param[4]*param[4]*facC*facC;
7048 seed->AddCovariance(covar);