1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
17 //-------------------------------------------------------
18 // Implementation of the TPC tracker
20 // Origin: Marian Ivanov Marian.Ivanov@cern.ch
22 // AliTPC parallel tracker
24 // The track fitting is based on Kalaman filtering approach
26 // The track finding steps:
27 // 1. Seeding - with and without vertex constraint
28 // - seeding with vertex constain done at first n^2 proble
29 // - seeding without vertex constraint n^3 problem
30 // 2. Tracking - follow prolongation road - find cluster - update kalman track
32 // The seeding and tracking is repeated several times, in different seeding region.
33 // This approach enables to find the track which cannot be seeded in some region of TPC
34 // This can happen because of low momenta (track do not reach outer radius), or track is currently in the ded region between sectors, or the track is for the moment overlapped with other track (seed quality is poor) ...
36 // With this approach we reach almost 100 % efficiency also for high occupancy events.
37 // (If the seeding efficiency in a region is about 90 % than with logical or of several
38 // regions we will reach 100% (in theory - supposing independence)
40 // Repeating several seeding - tracking procedures some of the tracks can be find
43 // The procedures to remove multi find tacks are impremented:
44 // RemoveUsed2 - fast procedure n problem -
45 // Algorithm - Sorting tracks according quality
46 // remove tracks with some shared fraction
47 // Sharing in respect to all tacks
48 // Signing clusters in gold region
49 // FindSplitted - slower algorithm n^2
50 // Sort the tracks according quality
51 // Loop over pair of tracks
52 // If overlap with other track bigger than threshold - remove track
54 // FindCurling - Finds the pair of tracks which are curling
55 // - About 10% of tracks can be find with this procedure
56 // The combinatorial background is too big to be used in High
57 // multiplicity environment
58 // - n^2 problem - Slow procedure - currently it is disabled because of
61 // The number of splitted tracks can be reduced disabling the sharing of the cluster.
62 // tpcRecoParam-> SetClusterSharing(kFALSE);
63 // IT IS HIGHLY non recomended to use it in high flux enviroonment
64 // Even using this switch some tracks can be found more than once
65 // (because of multiple seeding and low quality tracks which will not cross full chamber)
68 // The tracker itself can be debugged - the information about tracks can be stored in several // phases of the reconstruction
69 // To enable storage of the TPC tracks in the ESD friend track
70 // use AliTPCReconstructor::SetStreamLevel(n);
72 // The debug level - different procedure produce tree for numerical debugging
73 // To enable them set AliTPCReconstructor::SetStreamLevel(n); where nis bigger 1
77 // Adding systematic errors to the covariance:
79 // The systematic errors due to the misalignment and miscalibration are added to the covariance matrix
80 // of the tracks (not to the clusters as they are dependent):
81 // The parameters form AliTPCRecoParam are used AliTPCRecoParam::GetSystematicError
82 // The systematic errors are expressed there in RMS - position (cm), angle (rad), curvature (1/GeV)
83 // The default values are 0.
85 // The sytematic errors are added to the covariance matrix in following places:
87 // 1. During fisrt itteration - AliTPCtrackerMI::FillESD
88 // 2. Second iteration -
89 // 2.a ITS->TPC - AliTPCtrackerMI::ReadSeeds
90 // 2.b TPC->TRD - AliTPCtrackerMI::PropagateBack
91 // 3. Third iteration -
92 // 3.a TRD->TPC - AliTPCtrackerMI::ReadSeeds
93 // 3.b TPC->ITS - AliTPCtrackerMI::RefitInward
95 // There are several places in the code which can be numerically debuged
96 // This code is keeped in order to enable code development and to check the calibration implementtion
98 // 1. ErrParam stream - dump information about
100 // 2.a) cluster error estimate
101 // 3.a) cluster shape estimate
104 // Debug streamer levels:
106 //-------------------------------------------------------
111 #include "Riostream.h"
112 #include <TClonesArray.h>
114 #include <TObjArray.h>
116 #include <TGraphErrors.h>
118 #include "AliComplexCluster.h"
119 #include "AliESDEvent.h"
120 #include "AliESDtrack.h"
121 #include "AliESDVertex.h"
124 #include "AliHelix.h"
125 #include "AliRunLoader.h"
126 #include "AliTPCClustersRow.h"
127 #include "AliTPCParam.h"
128 #include "AliTPCReconstructor.h"
129 #include "AliTPCpolyTrack.h"
130 #include "AliTPCreco.h"
131 #include "AliTPCseed.h"
133 #include "AliTPCtrackerSector.h"
134 #include "AliTPCtrackerMI.h"
135 #include "TStopwatch.h"
136 #include "AliTPCReconstructor.h"
137 #include "AliAlignObj.h"
138 #include "AliTrackPointArray.h"
140 #include "AliTPCcalibDB.h"
141 #include "AliTPCcalibDButil.h"
142 #include "AliTPCTransform.h"
143 #include "AliTPCClusterParam.h"
144 #include "AliTPCdEdxInfo.h"
145 #include "AliDCSSensorArray.h"
146 #include "AliDCSSensor.h"
148 #include "AliCosmicTracker.h"
154 ClassImp(AliTPCtrackerMI)
158 class AliTPCFastMath {
161 static Double_t FastAsin(Double_t x);
163 static Double_t fgFastAsin[20000]; //lookup table for fast asin computation
166 Double_t AliTPCFastMath::fgFastAsin[20000];
167 AliTPCFastMath gAliTPCFastMath; // needed to fill the LUT
169 AliTPCFastMath::AliTPCFastMath(){
171 // initialized lookup table;
172 for (Int_t i=0;i<10000;i++){
173 fgFastAsin[2*i] = TMath::ASin(i/10000.);
174 fgFastAsin[2*i+1] = (TMath::ASin((i+1)/10000.)-fgFastAsin[2*i]);
178 Double_t AliTPCFastMath::FastAsin(Double_t x){
180 // return asin using lookup table
182 Int_t index = int(x*10000);
183 return fgFastAsin[2*index]+(x*10000.-index)*fgFastAsin[2*index+1];
186 Int_t index = int(x*10000);
187 return -(fgFastAsin[2*index]+(x*10000.-index)*fgFastAsin[2*index+1]);
189 //__________________________________________________________________
190 AliTPCtrackerMI::AliTPCtrackerMI()
217 // default constructor
219 for (Int_t irow=0; irow<200; irow++){
226 //_____________________________________________________________________
230 Int_t AliTPCtrackerMI::UpdateTrack(AliTPCseed * track, Int_t accept){
232 //update track information using current cluster - track->fCurrentCluster
235 AliTPCclusterMI* c =track->GetCurrentCluster();
236 if (accept > 0) //sign not accepted clusters
237 track->SetCurrentClusterIndex1(track->GetCurrentClusterIndex1() | 0x8000);
238 else // unsign accpeted clusters
239 track->SetCurrentClusterIndex1(track->GetCurrentClusterIndex1() & 0xffff7fff);
240 UInt_t i = track->GetCurrentClusterIndex1();
242 Int_t sec=(i&0xff000000)>>24;
243 //Int_t row = (i&0x00ff0000)>>16;
244 track->SetRow((i&0x00ff0000)>>16);
245 track->SetSector(sec);
246 // Int_t index = i&0xFFFF;
247 if (sec>=fkParam->GetNInnerSector()) track->SetRow(track->GetRow()+fkParam->GetNRowLow());
248 track->SetClusterIndex2(track->GetRow(), i);
249 //track->fFirstPoint = row;
250 //if ( track->fLastPoint<row) track->fLastPoint =row;
251 // if (track->fRow<0 || track->fRow>160) {
252 // printf("problem\n");
254 if (track->GetFirstPoint()>track->GetRow())
255 track->SetFirstPoint(track->GetRow());
256 if (track->GetLastPoint()<track->GetRow())
257 track->SetLastPoint(track->GetRow());
260 track->SetClusterPointer(track->GetRow(),c);
263 Double_t angle2 = track->GetSnp()*track->GetSnp();
265 //SET NEW Track Point
267 if (angle2<1) //PH sometimes angle2 is very big. To be investigated...
269 angle2 = TMath::Sqrt(angle2/(1-angle2));
270 AliTPCTrackerPoint &point =*(track->GetTrackPoint(track->GetRow()));
272 point.SetSigmaY(c->GetSigmaY2()/track->GetCurrentSigmaY2());
273 point.SetSigmaZ(c->GetSigmaZ2()/track->GetCurrentSigmaZ2());
274 point.SetErrY(sqrt(track->GetErrorY2()));
275 point.SetErrZ(sqrt(track->GetErrorZ2()));
277 point.SetX(track->GetX());
278 point.SetY(track->GetY());
279 point.SetZ(track->GetZ());
280 point.SetAngleY(angle2);
281 point.SetAngleZ(track->GetTgl());
282 if (point.IsShared()){
283 track->SetErrorY2(track->GetErrorY2()*4);
284 track->SetErrorZ2(track->GetErrorZ2()*4);
288 Double_t chi2 = track->GetPredictedChi2(track->GetCurrentCluster());
290 // track->SetErrorY2(track->GetErrorY2()*1.3);
291 // track->SetErrorY2(track->GetErrorY2()+0.01);
292 // track->SetErrorZ2(track->GetErrorZ2()*1.3);
293 // track->SetErrorZ2(track->GetErrorZ2()+0.005);
295 if (accept>0) return 0;
296 if (track->GetNumberOfClusters()%20==0){
297 // if (track->fHelixIn){
298 // TClonesArray & larr = *(track->fHelixIn);
299 // Int_t ihelix = larr.GetEntriesFast();
300 // new(larr[ihelix]) AliHelix(*track) ;
303 track->SetNoCluster(0);
304 return track->Update(c,chi2,i);
309 Int_t AliTPCtrackerMI::AcceptCluster(AliTPCseed * seed, AliTPCclusterMI * cluster)
312 // decide according desired precision to accept given
313 // cluster for tracking
315 seed->GetProlongation(cluster->GetX(),yt,zt);
316 Double_t sy2=ErrY2(seed,cluster);
317 Double_t sz2=ErrZ2(seed,cluster);
319 Double_t sdistancey2 = sy2+seed->GetSigmaY2();
320 Double_t sdistancez2 = sz2+seed->GetSigmaZ2();
321 Double_t dy=seed->GetCurrentCluster()->GetY()-yt;
322 Double_t dz=seed->GetCurrentCluster()->GetZ()-zt;
323 Double_t rdistancey2 = (seed->GetCurrentCluster()->GetY()-yt)*
324 (seed->GetCurrentCluster()->GetY()-yt)/sdistancey2;
325 Double_t rdistancez2 = (seed->GetCurrentCluster()->GetZ()-zt)*
326 (seed->GetCurrentCluster()->GetZ()-zt)/sdistancez2;
328 Double_t rdistance2 = rdistancey2+rdistancez2;
331 if (AliTPCReconstructor::StreamLevel()>2 && seed->GetNumberOfClusters()>20) {
332 Float_t rmsy2 = seed->GetCurrentSigmaY2();
333 Float_t rmsz2 = seed->GetCurrentSigmaZ2();
334 Float_t rmsy2p30 = seed->GetCMeanSigmaY2p30();
335 Float_t rmsz2p30 = seed->GetCMeanSigmaZ2p30();
336 Float_t rmsy2p30R = seed->GetCMeanSigmaY2p30R();
337 Float_t rmsz2p30R = seed->GetCMeanSigmaZ2p30R();
338 AliExternalTrackParam param(*seed);
339 static TVectorD gcl(3),gtr(3);
341 param.GetXYZ(gcl.GetMatrixArray());
342 cluster->GetGlobalXYZ(gclf);
343 gcl[0]=gclf[0]; gcl[1]=gclf[1]; gcl[2]=gclf[2];
346 if (AliTPCReconstructor::StreamLevel()>2) {
347 (*fDebugStreamer)<<"ErrParam"<<
360 "rmsy2p30="<<rmsy2p30<<
361 "rmsz2p30="<<rmsz2p30<<
362 "rmsy2p30R="<<rmsy2p30R<<
363 "rmsz2p30R="<<rmsz2p30R<<
364 // normalize distance -
365 "rdisty="<<rdistancey2<<
366 "rdistz="<<rdistancez2<<
367 "rdist="<<rdistance2<< //
371 //return 0; // temporary
372 if (rdistance2>32) return 3;
375 if ((rdistancey2>9. || rdistancez2>9.) && cluster->GetType()==0)
376 return 2; //suspisiouce - will be changed
378 if ((rdistancey2>6.25 || rdistancez2>6.25) && cluster->GetType()>0)
379 // strict cut on overlaped cluster
380 return 2; //suspisiouce - will be changed
382 if ( (rdistancey2>1. || rdistancez2>6.25 )
383 && cluster->GetType()<0){
384 seed->SetNFoundable(seed->GetNFoundable()-1);
388 if (fUseHLTClusters == 3 || fUseHLTClusters == 4) {
390 if(!AliTPCReconstructor::GetRecoParam()->GetUseHLTOnePadCluster()) {
391 if (TMath::Abs(cluster->GetSigmaY2()) < kAlmost0)
404 //_____________________________________________________________________________
405 AliTPCtrackerMI::AliTPCtrackerMI(const AliTPCParam *par):
407 fkNIS(par->GetNInnerSector()/2),
409 fkNOS(par->GetNOuterSector()/2),
431 //---------------------------------------------------------------------
432 // The main TPC tracker constructor
433 //---------------------------------------------------------------------
434 fInnerSec=new AliTPCtrackerSector[fkNIS];
435 fOuterSec=new AliTPCtrackerSector[fkNOS];
438 for (i=0; i<fkNIS; i++) fInnerSec[i].Setup(par,0);
439 for (i=0; i<fkNOS; i++) fOuterSec[i].Setup(par,1);
442 Int_t nrowlow = par->GetNRowLow();
443 Int_t nrowup = par->GetNRowUp();
446 for (i=0;i<nrowlow;i++){
447 fXRow[i] = par->GetPadRowRadiiLow(i);
448 fPadLength[i]= par->GetPadPitchLength(0,i);
449 fYMax[i] = fXRow[i]*TMath::Tan(0.5*par->GetInnerAngle());
453 for (i=0;i<nrowup;i++){
454 fXRow[i+nrowlow] = par->GetPadRowRadiiUp(i);
455 fPadLength[i+nrowlow] = par->GetPadPitchLength(60,i);
456 fYMax[i+nrowlow] = fXRow[i+nrowlow]*TMath::Tan(0.5*par->GetOuterAngle());
459 if (AliTPCReconstructor::StreamLevel()>0) {
460 fDebugStreamer = new TTreeSRedirector("TPCdebug.root");
463 fSeedsPool = new TClonesArray("AliTPCseed",1000);
465 //________________________________________________________________________
466 AliTPCtrackerMI::AliTPCtrackerMI(const AliTPCtrackerMI &t):
492 //------------------------------------
493 // dummy copy constructor
494 //------------------------------------------------------------------
496 for (Int_t irow=0; irow<200; irow++){
503 AliTPCtrackerMI & AliTPCtrackerMI::operator=(const AliTPCtrackerMI& /*r*/)
505 //------------------------------
507 //--------------------------------------------------------------
510 //_____________________________________________________________________________
511 AliTPCtrackerMI::~AliTPCtrackerMI() {
512 //------------------------------------------------------------------
513 // TPC tracker destructor
514 //------------------------------------------------------------------
521 if (fDebugStreamer) delete fDebugStreamer;
522 if (fSeedsPool) delete fSeedsPool;
526 void AliTPCtrackerMI::FillESD(const TObjArray* arr)
530 //fill esds using updated tracks
533 // write tracks to the event
534 // store index of the track
535 Int_t nseed=arr->GetEntriesFast();
536 //FindKinks(arr,fEvent);
537 for (Int_t i=0; i<nseed; i++) {
538 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
542 if (AliTPCReconstructor::StreamLevel()>1) {
543 (*fDebugStreamer)<<"Track0"<<
547 // pt->PropagateTo(fkParam->GetInnerRadiusLow());
548 if (pt->GetKinkIndex(0)<=0){ //don't propagate daughter tracks
549 pt->PropagateTo(fkParam->GetInnerRadiusLow());
552 if (( pt->GetPoints()[2]- pt->GetPoints()[0])>5 && pt->GetPoints()[3]>0.8){
554 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
555 iotrack.SetTPCPoints(pt->GetPoints());
556 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
557 iotrack.SetV0Indexes(pt->GetV0Indexes());
558 // iotrack.SetTPCpid(pt->fTPCr);
559 //iotrack.SetTPCindex(i);
560 fEvent->AddTrack(&iotrack);
564 if ( (pt->GetNumberOfClusters()>70)&& (Float_t(pt->GetNumberOfClusters())/Float_t(pt->GetNFoundable()))>0.55) {
566 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
567 iotrack.SetTPCPoints(pt->GetPoints());
568 //iotrack.SetTPCindex(i);
569 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
570 iotrack.SetV0Indexes(pt->GetV0Indexes());
571 // iotrack.SetTPCpid(pt->fTPCr);
572 fEvent->AddTrack(&iotrack);
576 // short tracks - maybe decays
578 if ( (pt->GetNumberOfClusters()>30) && (Float_t(pt->GetNumberOfClusters())/Float_t(pt->GetNFoundable()))>0.70) {
579 Int_t found,foundable,shared;
580 pt->GetClusterStatistic(0,60,found, foundable,shared,kFALSE);
581 if ( (found>20) && (pt->GetNShared()/float(pt->GetNumberOfClusters())<0.2)){
583 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
584 //iotrack.SetTPCindex(i);
585 iotrack.SetTPCPoints(pt->GetPoints());
586 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
587 iotrack.SetV0Indexes(pt->GetV0Indexes());
588 //iotrack.SetTPCpid(pt->fTPCr);
589 fEvent->AddTrack(&iotrack);
594 if ( (pt->GetNumberOfClusters()>20) && (Float_t(pt->GetNumberOfClusters())/Float_t(pt->GetNFoundable()))>0.8) {
595 Int_t found,foundable,shared;
596 pt->GetClusterStatistic(0,60,found, foundable,shared,kFALSE);
597 if (found<20) continue;
598 if (pt->GetNShared()/float(pt->GetNumberOfClusters())>0.2) continue;
601 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
602 iotrack.SetTPCPoints(pt->GetPoints());
603 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
604 iotrack.SetV0Indexes(pt->GetV0Indexes());
605 //iotrack.SetTPCpid(pt->fTPCr);
606 //iotrack.SetTPCindex(i);
607 fEvent->AddTrack(&iotrack);
610 // short tracks - secondaties
612 if ( (pt->GetNumberOfClusters()>30) ) {
613 Int_t found,foundable,shared;
614 pt->GetClusterStatistic(128,158,found, foundable,shared,kFALSE);
615 if ( (found>20) && (pt->GetNShared()/float(pt->GetNumberOfClusters())<0.2) &&float(found)/float(foundable)>0.8){
617 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
618 iotrack.SetTPCPoints(pt->GetPoints());
619 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
620 iotrack.SetV0Indexes(pt->GetV0Indexes());
621 //iotrack.SetTPCpid(pt->fTPCr);
622 //iotrack.SetTPCindex(i);
623 fEvent->AddTrack(&iotrack);
628 if ( (pt->GetNumberOfClusters()>15)) {
629 Int_t found,foundable,shared;
630 pt->GetClusterStatistic(138,158,found, foundable,shared,kFALSE);
631 if (found<15) continue;
632 if (foundable<=0) continue;
633 if (pt->GetNShared()/float(pt->GetNumberOfClusters())>0.2) continue;
634 if (float(found)/float(foundable)<0.8) continue;
637 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
638 iotrack.SetTPCPoints(pt->GetPoints());
639 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
640 iotrack.SetV0Indexes(pt->GetV0Indexes());
641 // iotrack.SetTPCpid(pt->fTPCr);
642 //iotrack.SetTPCindex(i);
643 fEvent->AddTrack(&iotrack);
647 // >> account for suppressed tracks in the kink indices (RS)
648 int nESDtracks = fEvent->GetNumberOfTracks();
649 for (int it=nESDtracks;it--;) {
650 AliESDtrack* esdTr = fEvent->GetTrack(it);
651 if (!esdTr || !esdTr->GetKinkIndex(0)) continue;
652 for (int ik=0;ik<3;ik++) {
654 if (!(knkId=esdTr->GetKinkIndex(ik))) break; // no more kinks for this track
655 AliESDkink* kink = fEvent->GetKink(TMath::Abs(knkId)-1);
657 AliError(Form("ESDTrack%d refers to non-existing kink %d",it,TMath::Abs(knkId)-1));
660 kink->SetIndex(it, knkId<0 ? 0:1); // update track index of the kink: mother at 0, daughter at 1
663 // << account for suppressed tracks in the kink indices (RS)
664 AliInfo(Form("Number of filled ESDs-\t%d\n",fEvent->GetNumberOfTracks()));
672 Double_t AliTPCtrackerMI::ErrY2(AliTPCseed* seed, const AliTPCclusterMI * cl){
675 // Use calibrated cluster error from OCDB
677 AliTPCClusterParam * clparam = AliTPCcalibDB::Instance()->GetClusterParam();
679 Float_t z = TMath::Abs(fkParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
680 Int_t ctype = cl->GetType();
681 Int_t type = (cl->GetRow()<63) ? 0: (cl->GetRow()>126) ? 1:2;
682 Double_t angle = seed->GetSnp()*seed->GetSnp();
683 angle = TMath::Sqrt(TMath::Abs(angle/(1.-angle)));
684 Double_t erry2 = clparam->GetError0Par(0,type, z,angle);
686 erry2+=0.5; // edge cluster
689 seed->SetErrorY2(erry2);
693 //calculate look-up table at the beginning
694 // static Bool_t ginit = kFALSE;
695 // static Float_t gnoise1,gnoise2,gnoise3;
696 // static Float_t ggg1[10000];
697 // static Float_t ggg2[10000];
698 // static Float_t ggg3[10000];
699 // static Float_t glandau1[10000];
700 // static Float_t glandau2[10000];
701 // static Float_t glandau3[10000];
703 // static Float_t gcor01[500];
704 // static Float_t gcor02[500];
705 // static Float_t gcorp[500];
709 // if (ginit==kFALSE){
710 // for (Int_t i=1;i<500;i++){
711 // Float_t rsigma = float(i)/100.;
712 // gcor02[i] = TMath::Max(0.78 +TMath::Exp(7.4*(rsigma-1.2)),0.6);
713 // gcor01[i] = TMath::Max(0.72 +TMath::Exp(3.36*(rsigma-1.2)),0.6);
714 // gcorp[i] = TMath::Max(TMath::Power((rsigma+0.5),1.5),1.2);
718 // for (Int_t i=3;i<10000;i++){
722 // Float_t amp = float(i);
723 // Float_t padlength =0.75;
724 // gnoise1 = 0.0004/padlength;
725 // Float_t nel = 0.268*amp;
726 // Float_t nprim = 0.155*amp;
727 // ggg1[i] = fkParam->GetDiffT()*fkParam->GetDiffT()*(2+0.001*nel/(padlength*padlength))/nel;
728 // glandau1[i] = (2.+0.12*nprim)*0.5* (2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
729 // if (glandau1[i]>1) glandau1[i]=1;
730 // glandau1[i]*=padlength*padlength/12.;
734 // gnoise2 = 0.0004/padlength;
736 // nprim = 0.133*amp;
737 // ggg2[i] = fkParam->GetDiffT()*fkParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
738 // glandau2[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
739 // if (glandau2[i]>1) glandau2[i]=1;
740 // glandau2[i]*=padlength*padlength/12.;
745 // gnoise3 = 0.0004/padlength;
747 // nprim = 0.133*amp;
748 // ggg3[i] = fkParam->GetDiffT()*fkParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
749 // glandau3[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
750 // if (glandau3[i]>1) glandau3[i]=1;
751 // glandau3[i]*=padlength*padlength/12.;
759 // Int_t amp = int(TMath::Abs(cl->GetQ()));
761 // seed->SetErrorY2(1.);
765 // Float_t z = TMath::Abs(fkParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
766 // Int_t ctype = cl->GetType();
767 // Float_t padlength= GetPadPitchLength(seed->GetRow());
768 // Double_t angle2 = seed->GetSnp()*seed->GetSnp();
769 // angle2 = angle2/(1-angle2);
771 // //cluster "quality"
772 // Int_t rsigmay = int(100.*cl->GetSigmaY2()/(seed->GetCurrentSigmaY2()));
775 // if (fSectors==fInnerSec){
776 // snoise2 = gnoise1;
777 // res = ggg1[amp]*z+glandau1[amp]*angle2;
778 // if (ctype==0) res *= gcor01[rsigmay];
781 // res*= gcorp[rsigmay];
785 // if (padlength<1.1){
786 // snoise2 = gnoise2;
787 // res = ggg2[amp]*z+glandau2[amp]*angle2;
788 // if (ctype==0) res *= gcor02[rsigmay];
791 // res*= gcorp[rsigmay];
795 // snoise2 = gnoise3;
796 // res = ggg3[amp]*z+glandau3[amp]*angle2;
797 // if (ctype==0) res *= gcor02[rsigmay];
800 // res*= gcorp[rsigmay];
807 // res*=2.4; // overestimate error 2 times
811 // if (res<2*snoise2)
814 // seed->SetErrorY2(res);
822 Double_t AliTPCtrackerMI::ErrZ2(AliTPCseed* seed, const AliTPCclusterMI * cl){
825 // Use calibrated cluster error from OCDB
827 AliTPCClusterParam * clparam = AliTPCcalibDB::Instance()->GetClusterParam();
829 Float_t z = TMath::Abs(fkParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
830 Int_t ctype = cl->GetType();
831 Int_t type = (cl->GetRow()<63) ? 0: (cl->GetRow()>126) ? 1:2;
833 Double_t angle2 = seed->GetSnp()*seed->GetSnp();
834 angle2 = seed->GetTgl()*seed->GetTgl()*(1+angle2/(1-angle2));
835 Double_t angle = TMath::Sqrt(TMath::Abs(angle2));
836 Double_t errz2 = clparam->GetError0Par(1,type, z,angle);
838 errz2+=0.5; // edge cluster
841 seed->SetErrorZ2(errz2);
847 // //seed->SetErrorY2(0.1);
849 // //calculate look-up table at the beginning
850 // static Bool_t ginit = kFALSE;
851 // static Float_t gnoise1,gnoise2,gnoise3;
852 // static Float_t ggg1[10000];
853 // static Float_t ggg2[10000];
854 // static Float_t ggg3[10000];
855 // static Float_t glandau1[10000];
856 // static Float_t glandau2[10000];
857 // static Float_t glandau3[10000];
859 // static Float_t gcor01[1000];
860 // static Float_t gcor02[1000];
861 // static Float_t gcorp[1000];
865 // if (ginit==kFALSE){
866 // for (Int_t i=1;i<1000;i++){
867 // Float_t rsigma = float(i)/100.;
868 // gcor02[i] = TMath::Max(0.81 +TMath::Exp(6.8*(rsigma-1.2)),0.6);
869 // gcor01[i] = TMath::Max(0.72 +TMath::Exp(2.04*(rsigma-1.2)),0.6);
870 // gcorp[i] = TMath::Max(TMath::Power((rsigma+0.5),1.5),1.2);
874 // for (Int_t i=3;i<10000;i++){
878 // Float_t amp = float(i);
879 // Float_t padlength =0.75;
880 // gnoise1 = 0.0004/padlength;
881 // Float_t nel = 0.268*amp;
882 // Float_t nprim = 0.155*amp;
883 // ggg1[i] = fkParam->GetDiffT()*fkParam->GetDiffT()*(2+0.001*nel/(padlength*padlength))/nel;
884 // glandau1[i] = (2.+0.12*nprim)*0.5* (2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
885 // if (glandau1[i]>1) glandau1[i]=1;
886 // glandau1[i]*=padlength*padlength/12.;
890 // gnoise2 = 0.0004/padlength;
892 // nprim = 0.133*amp;
893 // ggg2[i] = fkParam->GetDiffT()*fkParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
894 // glandau2[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
895 // if (glandau2[i]>1) glandau2[i]=1;
896 // glandau2[i]*=padlength*padlength/12.;
901 // gnoise3 = 0.0004/padlength;
903 // nprim = 0.133*amp;
904 // ggg3[i] = fkParam->GetDiffT()*fkParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
905 // glandau3[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
906 // if (glandau3[i]>1) glandau3[i]=1;
907 // glandau3[i]*=padlength*padlength/12.;
915 // Int_t amp = int(TMath::Abs(cl->GetQ()));
917 // seed->SetErrorY2(1.);
921 // Float_t z = TMath::Abs(fkParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
922 // Int_t ctype = cl->GetType();
923 // Float_t padlength= GetPadPitchLength(seed->GetRow());
925 // Double_t angle2 = seed->GetSnp()*seed->GetSnp();
926 // // if (angle2<0.6) angle2 = 0.6;
927 // angle2 = seed->GetTgl()*seed->GetTgl()*(1+angle2/(1-angle2));
929 // //cluster "quality"
930 // Int_t rsigmaz = int(100.*cl->GetSigmaZ2()/(seed->GetCurrentSigmaZ2()));
933 // if (fSectors==fInnerSec){
934 // snoise2 = gnoise1;
935 // res = ggg1[amp]*z+glandau1[amp]*angle2;
936 // if (ctype==0) res *= gcor01[rsigmaz];
939 // res*= gcorp[rsigmaz];
943 // if (padlength<1.1){
944 // snoise2 = gnoise2;
945 // res = ggg2[amp]*z+glandau2[amp]*angle2;
946 // if (ctype==0) res *= gcor02[rsigmaz];
949 // res*= gcorp[rsigmaz];
953 // snoise2 = gnoise3;
954 // res = ggg3[amp]*z+glandau3[amp]*angle2;
955 // if (ctype==0) res *= gcor02[rsigmaz];
958 // res*= gcorp[rsigmaz];
967 // if ((ctype<0) &&<70){
972 // if (res<2*snoise2)
974 // if (res>3) res =3;
975 // seed->SetErrorZ2(res);
983 void AliTPCtrackerMI::RotateToLocal(AliTPCseed *seed)
985 //rotate to track "local coordinata
986 Float_t x = seed->GetX();
987 Float_t y = seed->GetY();
988 Float_t ymax = x*TMath::Tan(0.5*fSectors->GetAlpha());
991 seed->SetRelativeSector((seed->GetRelativeSector()+1) % fN);
992 if (!seed->Rotate(fSectors->GetAlpha()))
994 } else if (y <-ymax) {
995 seed->SetRelativeSector((seed->GetRelativeSector()-1+fN) % fN);
996 if (!seed->Rotate(-fSectors->GetAlpha()))
1004 //_____________________________________________________________________________
1005 Double_t AliTPCtrackerMI::F1old(Double_t x1,Double_t y1,
1006 Double_t x2,Double_t y2,
1007 Double_t x3,Double_t y3) const
1009 //-----------------------------------------------------------------
1010 // Initial approximation of the track curvature
1011 //-----------------------------------------------------------------
1012 Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
1013 Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
1014 (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
1015 Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
1016 (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
1018 Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
1019 if ( xr*xr+yr*yr<=0.00000000000001) return 100;
1020 return -xr*yr/sqrt(xr*xr+yr*yr);
1025 //_____________________________________________________________________________
1026 Double_t AliTPCtrackerMI::F1(Double_t x1,Double_t y1,
1027 Double_t x2,Double_t y2,
1028 Double_t x3,Double_t y3) const
1030 //-----------------------------------------------------------------
1031 // Initial approximation of the track curvature
1032 //-----------------------------------------------------------------
1038 Double_t det = x3*y2-x2*y3;
1039 if (TMath::Abs(det)<1e-10){
1043 Double_t u = 0.5* (x2*(x2-x3)+y2*(y2-y3))/det;
1044 Double_t x0 = x3*0.5-y3*u;
1045 Double_t y0 = y3*0.5+x3*u;
1046 Double_t c2 = 1/TMath::Sqrt(x0*x0+y0*y0);
1052 Double_t AliTPCtrackerMI::F2(Double_t x1,Double_t y1,
1053 Double_t x2,Double_t y2,
1054 Double_t x3,Double_t y3) const
1056 //-----------------------------------------------------------------
1057 // Initial approximation of the track curvature
1058 //-----------------------------------------------------------------
1064 Double_t det = x3*y2-x2*y3;
1065 if (TMath::Abs(det)<1e-10) {
1069 Double_t u = 0.5* (x2*(x2-x3)+y2*(y2-y3))/det;
1070 Double_t x0 = x3*0.5-y3*u;
1071 Double_t y0 = y3*0.5+x3*u;
1072 Double_t c2 = 1/TMath::Sqrt(x0*x0+y0*y0);
1081 //_____________________________________________________________________________
1082 Double_t AliTPCtrackerMI::F2old(Double_t x1,Double_t y1,
1083 Double_t x2,Double_t y2,
1084 Double_t x3,Double_t y3) const
1086 //-----------------------------------------------------------------
1087 // Initial approximation of the track curvature times center of curvature
1088 //-----------------------------------------------------------------
1089 Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
1090 Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
1091 (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
1092 Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
1093 (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
1095 Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
1097 return -a/(d*y1-b)*xr/sqrt(xr*xr+yr*yr);
1100 //_____________________________________________________________________________
1101 Double_t AliTPCtrackerMI::F3(Double_t x1,Double_t y1,
1102 Double_t x2,Double_t y2,
1103 Double_t z1,Double_t z2) const
1105 //-----------------------------------------------------------------
1106 // Initial approximation of the tangent of the track dip angle
1107 //-----------------------------------------------------------------
1108 return (z1 - z2)/sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
1112 Double_t AliTPCtrackerMI::F3n(Double_t x1,Double_t y1,
1113 Double_t x2,Double_t y2,
1114 Double_t z1,Double_t z2, Double_t c) const
1116 //-----------------------------------------------------------------
1117 // Initial approximation of the tangent of the track dip angle
1118 //-----------------------------------------------------------------
1122 //angle1 = (z1-z2)*c/(TMath::ASin(c*x1-ni)-TMath::ASin(c*x2-ni));
1124 Double_t d = TMath::Sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
1125 if (TMath::Abs(d*c*0.5)>1) return 0;
1126 // Double_t angle2 = TMath::ASin(d*c*0.5);
1127 // Double_t angle2 = AliTPCFastMath::FastAsin(d*c*0.5);
1128 Double_t angle2 = (d*c*0.5>0.1)? TMath::ASin(d*c*0.5): AliTPCFastMath::FastAsin(d*c*0.5);
1130 angle2 = (z1-z2)*c/(angle2*2.);
1134 Bool_t AliTPCtrackerMI::GetProlongation(Double_t x1, Double_t x2, Double_t x[5], Double_t &y, Double_t &z) const
1135 {//-----------------------------------------------------------------
1136 // This function find proloncation of a track to a reference plane x=x2.
1137 //-----------------------------------------------------------------
1141 if (TMath::Abs(x[4]*x1 - x[2]) >= 0.999) {
1145 Double_t c1=x[4]*x1 - x[2], r1=TMath::Sqrt((1.-c1)*(1.+c1));
1146 Double_t c2=x[4]*x2 - x[2], r2=TMath::Sqrt((1.-c2)*(1.+c2));
1150 Double_t dy = dx*(c1+c2)/(r1+r2);
1153 Double_t delta = x[4]*dx*(c1+c2)/(c1*r2 + c2*r1);
1155 if (TMath::Abs(delta)>0.01){
1156 dz = x[3]*TMath::ASin(delta)/x[4];
1158 dz = x[3]*AliTPCFastMath::FastAsin(delta)/x[4];
1161 //dz = x[3]*AliTPCFastMath::FastAsin(delta)/x[4];
1169 Int_t AliTPCtrackerMI::LoadClusters (TTree *const tree)
1174 return LoadClusters();
1178 Int_t AliTPCtrackerMI::LoadClusters(const TObjArray *arr)
1181 // load clusters to the memory
1182 AliTPCClustersRow *clrow = new AliTPCClustersRow("AliTPCclusterMI");
1183 Int_t lower = arr->LowerBound();
1184 Int_t entries = arr->GetEntriesFast();
1186 for (Int_t i=lower; i<entries; i++) {
1187 clrow = (AliTPCClustersRow*) arr->At(i);
1188 if(!clrow) continue;
1189 if(!clrow->GetArray()) continue;
1193 fkParam->AdjustSectorRow(clrow->GetID(),sec,row);
1195 for (Int_t icl=0; icl<clrow->GetArray()->GetEntriesFast(); icl++){
1196 Transform((AliTPCclusterMI*)(clrow->GetArray()->At(icl)));
1199 if (clrow->GetArray()->GetEntriesFast()<=0) continue;
1200 AliTPCtrackerRow * tpcrow=0;
1203 tpcrow = &(fInnerSec[sec%fkNIS][row]);
1207 tpcrow = &(fOuterSec[(sec-fkNIS*2)%fkNOS][row]);
1208 left = (sec-fkNIS*2)/fkNOS;
1211 tpcrow->SetN1(clrow->GetArray()->GetEntriesFast());
1212 for (Int_t j=0;j<tpcrow->GetN1();++j)
1213 tpcrow->SetCluster1(j, *(AliTPCclusterMI*)(clrow->GetArray()->At(j)));
1216 tpcrow->SetN2(clrow->GetArray()->GetEntriesFast());
1217 for (Int_t j=0;j<tpcrow->GetN2();++j)
1218 tpcrow->SetCluster2(j, *(AliTPCclusterMI*)(clrow->GetArray()->At(j)));
1220 clrow->GetArray()->Clear("C");
1229 Int_t AliTPCtrackerMI::LoadClusters(const TClonesArray *arr)
1232 // load clusters to the memory from one
1235 AliTPCclusterMI *clust=0;
1236 Int_t count[72][96] = { {0} , {0} };
1238 // loop over clusters
1239 for (Int_t icl=0; icl<arr->GetEntriesFast(); icl++) {
1240 clust = (AliTPCclusterMI*)arr->At(icl);
1241 if(!clust) continue;
1242 //printf("cluster: det %d, row %d \n", clust->GetDetector(),clust->GetRow());
1244 // transform clusters
1247 // count clusters per pad row
1248 count[clust->GetDetector()][clust->GetRow()]++;
1251 // insert clusters to sectors
1252 for (Int_t icl=0; icl<arr->GetEntriesFast(); icl++) {
1253 clust = (AliTPCclusterMI*)arr->At(icl);
1254 if(!clust) continue;
1256 Int_t sec = clust->GetDetector();
1257 Int_t row = clust->GetRow();
1259 // filter overlapping pad rows needed by HLT
1260 if(sec<fkNIS*2) { //IROCs
1261 if(row == 30) continue;
1264 if(row == 27 || row == 76) continue;
1269 // left = sec/fkNIS;
1270 fInnerSec[sec%fkNIS].InsertCluster(clust, count[sec][row], fkParam);
1273 // left = (sec-fkNIS*2)/fkNOS;
1274 fOuterSec[(sec-fkNIS*2)%fkNOS].InsertCluster(clust, count[sec][row], fkParam);
1278 // Load functions must be called behind LoadCluster(TClonesArray*)
1280 //LoadOuterSectors();
1281 //LoadInnerSectors();
1287 Int_t AliTPCtrackerMI::LoadClusters()
1290 // load clusters to the memory
1291 static AliTPCClustersRow *clrow= new AliTPCClustersRow("AliTPCclusterMI");
1293 // TTree * tree = fClustersArray.GetTree();
1295 TTree * tree = fInput;
1296 TBranch * br = tree->GetBranch("Segment");
1297 br->SetAddress(&clrow);
1299 // Conversion of pad, row coordinates in local tracking coords.
1300 // Could be skipped here; is already done in clusterfinder
1302 Int_t j=Int_t(tree->GetEntries());
1303 for (Int_t i=0; i<j; i++) {
1307 fkParam->AdjustSectorRow(clrow->GetID(),sec,row);
1308 for (Int_t icl=0; icl<clrow->GetArray()->GetEntriesFast(); icl++){
1309 Transform((AliTPCclusterMI*)(clrow->GetArray()->At(icl)));
1312 AliTPCtrackerRow * tpcrow=0;
1315 tpcrow = &(fInnerSec[sec%fkNIS][row]);
1319 tpcrow = &(fOuterSec[(sec-fkNIS*2)%fkNOS][row]);
1320 left = (sec-fkNIS*2)/fkNOS;
1323 tpcrow->SetN1(clrow->GetArray()->GetEntriesFast());
1324 for (Int_t k=0;k<tpcrow->GetN1();++k)
1325 tpcrow->SetCluster1(k, *(AliTPCclusterMI*)(clrow->GetArray()->At(k)));
1328 tpcrow->SetN2(clrow->GetArray()->GetEntriesFast());
1329 for (Int_t k=0;k<tpcrow->GetN2();++k)
1330 tpcrow->SetCluster2(k, *(AliTPCclusterMI*)(clrow->GetArray()->At(k)));
1341 void AliTPCtrackerMI::UnloadClusters()
1344 // unload clusters from the memory
1346 Int_t nrows = fOuterSec->GetNRows();
1347 for (Int_t sec = 0;sec<fkNOS;sec++)
1348 for (Int_t row = 0;row<nrows;row++){
1349 AliTPCtrackerRow* tpcrow = &(fOuterSec[sec%fkNOS][row]);
1351 // if (tpcrow->fClusters1) delete []tpcrow->fClusters1;
1352 // if (tpcrow->fClusters2) delete []tpcrow->fClusters2;
1354 tpcrow->ResetClusters();
1357 nrows = fInnerSec->GetNRows();
1358 for (Int_t sec = 0;sec<fkNIS;sec++)
1359 for (Int_t row = 0;row<nrows;row++){
1360 AliTPCtrackerRow* tpcrow = &(fInnerSec[sec%fkNIS][row]);
1362 // if (tpcrow->fClusters1) delete []tpcrow->fClusters1;
1363 //if (tpcrow->fClusters2) delete []tpcrow->fClusters2;
1365 tpcrow->ResetClusters();
1371 void AliTPCtrackerMI::FillClusterArray(TObjArray* array) const{
1373 // Filling cluster to the array - For visualization purposes
1376 nrows = fOuterSec->GetNRows();
1377 for (Int_t sec = 0;sec<fkNOS;sec++)
1378 for (Int_t row = 0;row<nrows;row++){
1379 AliTPCtrackerRow* tpcrow = &(fOuterSec[sec%fkNOS][row]);
1380 if (!tpcrow) continue;
1381 for (Int_t icl = 0;icl<tpcrow->GetN();icl++){
1382 array->AddLast((TObject*)((*tpcrow)[icl]));
1385 nrows = fInnerSec->GetNRows();
1386 for (Int_t sec = 0;sec<fkNIS;sec++)
1387 for (Int_t row = 0;row<nrows;row++){
1388 AliTPCtrackerRow* tpcrow = &(fInnerSec[sec%fkNIS][row]);
1389 if (!tpcrow) continue;
1390 for (Int_t icl = 0;icl<tpcrow->GetN();icl++){
1391 array->AddLast((TObject*)(*tpcrow)[icl]);
1397 void AliTPCtrackerMI::Transform(AliTPCclusterMI * cluster){
1401 AliTPCcalibDB * calibDB = AliTPCcalibDB::Instance();
1402 AliTPCTransform *transform = calibDB->GetTransform() ;
1404 AliFatal("Tranformations not in calibDB");
1407 transform->SetCurrentRecoParam((AliTPCRecoParam*)AliTPCReconstructor::GetRecoParam());
1408 Double_t x[3]={cluster->GetRow(),cluster->GetPad(),cluster->GetTimeBin()};
1409 Int_t i[1]={cluster->GetDetector()};
1410 transform->Transform(x,i,0,1);
1411 // if (cluster->GetDetector()%36>17){
1416 // in debug mode check the transformation
1418 if (AliTPCReconstructor::StreamLevel()>2) {
1420 cluster->GetGlobalXYZ(gx);
1421 Int_t event = (fEvent==NULL)? 0: fEvent->GetEventNumberInFile();
1422 TTreeSRedirector &cstream = *fDebugStreamer;
1423 cstream<<"Transform"<<
1434 cluster->SetX(x[0]);
1435 cluster->SetY(x[1]);
1436 cluster->SetZ(x[2]);
1441 //if (!fkParam->IsGeoRead()) fkParam->ReadGeoMatrices();
1442 if (AliTPCReconstructor::GetRecoParam()->GetUseSectorAlignment() && (!calibDB->HasAlignmentOCDB())){
1443 TGeoHMatrix *mat = fkParam->GetClusterMatrix(cluster->GetDetector());
1445 Double_t pos[3]= {cluster->GetX(),cluster->GetY(),cluster->GetZ()};
1446 Double_t posC[3]={cluster->GetX(),cluster->GetY(),cluster->GetZ()};
1447 if (mat) mat->LocalToMaster(pos,posC);
1449 // chack Loading of Geo matrices from GeoManager - TEMPORARY FIX
1451 cluster->SetX(posC[0]);
1452 cluster->SetY(posC[1]);
1453 cluster->SetZ(posC[2]);
1457 //_____________________________________________________________________________
1458 Int_t AliTPCtrackerMI::LoadOuterSectors() {
1459 //-----------------------------------------------------------------
1460 // This function fills outer TPC sectors with clusters.
1461 //-----------------------------------------------------------------
1462 Int_t nrows = fOuterSec->GetNRows();
1464 for (Int_t sec = 0;sec<fkNOS;sec++)
1465 for (Int_t row = 0;row<nrows;row++){
1466 AliTPCtrackerRow* tpcrow = &(fOuterSec[sec%fkNOS][row]);
1467 Int_t sec2 = sec+2*fkNIS;
1469 Int_t ncl = tpcrow->GetN1();
1471 AliTPCclusterMI *c= (tpcrow->GetCluster1(ncl));
1472 index=(((sec2<<8)+row)<<16)+ncl;
1473 tpcrow->InsertCluster(c,index);
1476 ncl = tpcrow->GetN2();
1478 AliTPCclusterMI *c= (tpcrow->GetCluster2(ncl));
1479 index=((((sec2+fkNOS)<<8)+row)<<16)+ncl;
1480 tpcrow->InsertCluster(c,index);
1483 // write indexes for fast acces
1485 for (Int_t i=0;i<510;i++)
1486 tpcrow->SetFastCluster(i,-1);
1487 for (Int_t i=0;i<tpcrow->GetN();i++){
1488 Int_t zi = Int_t((*tpcrow)[i]->GetZ()+255.);
1489 tpcrow->SetFastCluster(zi,i); // write index
1492 for (Int_t i=0;i<510;i++){
1493 if (tpcrow->GetFastCluster(i)<0)
1494 tpcrow->SetFastCluster(i,last);
1496 last = tpcrow->GetFastCluster(i);
1505 //_____________________________________________________________________________
1506 Int_t AliTPCtrackerMI::LoadInnerSectors() {
1507 //-----------------------------------------------------------------
1508 // This function fills inner TPC sectors with clusters.
1509 //-----------------------------------------------------------------
1510 Int_t nrows = fInnerSec->GetNRows();
1512 for (Int_t sec = 0;sec<fkNIS;sec++)
1513 for (Int_t row = 0;row<nrows;row++){
1514 AliTPCtrackerRow* tpcrow = &(fInnerSec[sec%fkNIS][row]);
1517 Int_t ncl = tpcrow->GetN1();
1519 AliTPCclusterMI *c= (tpcrow->GetCluster1(ncl));
1520 index=(((sec<<8)+row)<<16)+ncl;
1521 tpcrow->InsertCluster(c,index);
1524 ncl = tpcrow->GetN2();
1526 AliTPCclusterMI *c= (tpcrow->GetCluster2(ncl));
1527 index=((((sec+fkNIS)<<8)+row)<<16)+ncl;
1528 tpcrow->InsertCluster(c,index);
1531 // write indexes for fast acces
1533 for (Int_t i=0;i<510;i++)
1534 tpcrow->SetFastCluster(i,-1);
1535 for (Int_t i=0;i<tpcrow->GetN();i++){
1536 Int_t zi = Int_t((*tpcrow)[i]->GetZ()+255.);
1537 tpcrow->SetFastCluster(zi,i); // write index
1540 for (Int_t i=0;i<510;i++){
1541 if (tpcrow->GetFastCluster(i)<0)
1542 tpcrow->SetFastCluster(i,last);
1544 last = tpcrow->GetFastCluster(i);
1556 //_________________________________________________________________________
1557 AliTPCclusterMI *AliTPCtrackerMI::GetClusterMI(Int_t index) const {
1558 //--------------------------------------------------------------------
1559 // Return pointer to a given cluster
1560 //--------------------------------------------------------------------
1561 if (index<0) return 0; // no cluster
1562 Int_t sec=(index&0xff000000)>>24;
1563 Int_t row=(index&0x00ff0000)>>16;
1564 Int_t ncl=(index&0x00007fff)>>00;
1566 const AliTPCtrackerRow * tpcrow=0;
1567 TClonesArray * clrow =0;
1569 if (sec<0 || sec>=fkNIS*4) {
1570 AliWarning(Form("Wrong sector %d",sec));
1575 AliTPCtrackerSector& tracksec = fInnerSec[sec%fkNIS];
1576 if (tracksec.GetNRows()<=row) return 0;
1577 tpcrow = &(tracksec[row]);
1578 if (tpcrow==0) return 0;
1581 if (tpcrow->GetN1()<=ncl) return 0;
1582 clrow = tpcrow->GetClusters1();
1585 if (tpcrow->GetN2()<=ncl) return 0;
1586 clrow = tpcrow->GetClusters2();
1590 AliTPCtrackerSector& tracksec = fOuterSec[(sec-fkNIS*2)%fkNOS];
1591 if (tracksec.GetNRows()<=row) return 0;
1592 tpcrow = &(tracksec[row]);
1593 if (tpcrow==0) return 0;
1595 if (sec-2*fkNIS<fkNOS) {
1596 if (tpcrow->GetN1()<=ncl) return 0;
1597 clrow = tpcrow->GetClusters1();
1600 if (tpcrow->GetN2()<=ncl) return 0;
1601 clrow = tpcrow->GetClusters2();
1605 return (AliTPCclusterMI*)clrow->At(ncl);
1611 Int_t AliTPCtrackerMI::FollowToNext(AliTPCseed& t, Int_t nr) {
1612 //-----------------------------------------------------------------
1613 // This function tries to find a track prolongation to next pad row
1614 //-----------------------------------------------------------------
1616 Double_t x= GetXrow(nr), ymax=GetMaxY(nr);
1619 AliTPCclusterMI *cl=0;
1620 Int_t tpcindex= t.GetClusterIndex2(nr);
1622 // update current shape info every 5 pad-row
1623 // if ( (nr%5==0) || t.GetNumberOfClusters()<2 || (t.fCurrentSigmaY2<0.0001) ){
1627 if (fIteration>0 && tpcindex>=-1){ //if we have already clusters
1629 if (tpcindex==-1) return 0; //track in dead zone
1630 if (tpcindex >= 0){ //
1631 cl = t.GetClusterPointer(nr);
1632 //if (cl==0) cl = GetClusterMI(tpcindex);
1633 if (!cl) cl = GetClusterMI(tpcindex);
1634 t.SetCurrentClusterIndex1(tpcindex);
1637 Int_t relativesector = ((tpcindex&0xff000000)>>24)%18; // if previously accepted cluster in different sector
1638 Float_t angle = relativesector*fSectors->GetAlpha()+fSectors->GetAlphaShift();
1640 if (angle<-TMath::Pi()) angle += 2*TMath::Pi();
1641 if (angle>=TMath::Pi()) angle -= 2*TMath::Pi();
1643 if (TMath::Abs(angle-t.GetAlpha())>0.001){
1644 Double_t rotation = angle-t.GetAlpha();
1645 t.SetRelativeSector(relativesector);
1646 if (!t.Rotate(rotation)) {
1647 t.SetClusterIndex(nr, t.GetClusterIndex(nr) | 0x8000);
1651 if (!t.PropagateTo(x)) {
1652 t.SetClusterIndex(nr, t.GetClusterIndex(nr) | 0x8000);
1656 t.SetCurrentCluster(cl);
1658 Int_t accept = AcceptCluster(&t,t.GetCurrentCluster());
1659 if ((tpcindex&0x8000)==0) accept =0;
1661 //if founded cluster is acceptible
1662 if (cl->IsUsed(11)) { // id cluster is shared inrease uncertainty
1663 t.SetErrorY2(t.GetErrorY2()+0.03);
1664 t.SetErrorZ2(t.GetErrorZ2()+0.03);
1665 t.SetErrorY2(t.GetErrorY2()*3);
1666 t.SetErrorZ2(t.GetErrorZ2()*3);
1668 t.SetNFoundable(t.GetNFoundable()+1);
1669 UpdateTrack(&t,accept);
1672 else { // Remove old cluster from track
1673 t.SetClusterIndex(nr, -3);
1674 t.SetClusterPointer(nr, 0);
1678 if (TMath::Abs(t.GetSnp())>AliTPCReconstructor::GetMaxSnpTracker()) return 0; // cut on angle
1679 if (fIteration>1 && IsFindable(t)){
1680 // not look for new cluster during refitting
1681 t.SetNFoundable(t.GetNFoundable()+1);
1686 // if (TMath::Abs(t.GetSnp())>0.95 || TMath::Abs(x*t.GetC()-t.GetEta())>0.95) return 0;// patch 28 fev 06
1687 if (!t.PropagateTo(x)) {
1688 if (fIteration==0) t.SetRemoval(10);
1691 Double_t y = t.GetY();
1692 if (TMath::Abs(y)>ymax){
1694 t.SetRelativeSector((t.GetRelativeSector()+1) % fN);
1695 if (!t.Rotate(fSectors->GetAlpha()))
1697 } else if (y <-ymax) {
1698 t.SetRelativeSector((t.GetRelativeSector()-1+fN) % fN);
1699 if (!t.Rotate(-fSectors->GetAlpha()))
1702 if (!t.PropagateTo(x)) {
1703 if (fIteration==0) t.SetRemoval(10);
1709 Double_t z=t.GetZ();
1712 if (!IsActive(t.GetRelativeSector(),nr)) {
1714 t.SetClusterIndex2(nr,-1);
1717 //AliInfo(Form("A - Sector%d phi %f - alpha %f", t.fRelativeSector,y/x, t.GetAlpha()));
1718 Bool_t isActive = IsActive(t.GetRelativeSector(),nr);
1719 Bool_t isActive2 = (nr>=fInnerSec->GetNRows()) ? fOuterSec[t.GetRelativeSector()][nr-fInnerSec->GetNRows()].GetN()>0:fInnerSec[t.GetRelativeSector()][nr].GetN()>0;
1721 if (!isActive || !isActive2) return 0;
1723 const AliTPCtrackerRow &krow=GetRow(t.GetRelativeSector(),nr);
1724 if ( (t.GetSigmaY2()<0) || t.GetSigmaZ2()<0) return 0;
1726 Double_t roadz = 1.;
1728 if (TMath::Abs(TMath::Abs(y)-ymax)<krow.GetDeadZone()){
1730 t.SetClusterIndex2(nr,-1);
1736 // if (TMath::Abs(z)<(AliTPCReconstructor::GetCtgRange()*x+10) && TMath::Abs(z)<fkParam->GetZLength(0) && (TMath::Abs(t.GetSnp())<AliTPCReconstructor::GetMaxSnpTracker()))
1737 t.SetNFoundable(t.GetNFoundable()+1);
1743 // cl = krow.FindNearest2(y+10.,z,roady,roadz,index);
1744 cl = krow.FindNearest2(y,z,roady,roadz,index);
1745 if (cl) t.SetCurrentClusterIndex1(krow.GetIndex(index));
1748 t.SetCurrentCluster(cl);
1750 if (fIteration==2&&cl->IsUsed(10)) return 0;
1751 Int_t accept = AcceptCluster(&t,t.GetCurrentCluster());
1752 if (fIteration==2&&cl->IsUsed(11)) {
1753 t.SetErrorY2(t.GetErrorY2()+0.03);
1754 t.SetErrorZ2(t.GetErrorZ2()+0.03);
1755 t.SetErrorY2(t.GetErrorY2()*3);
1756 t.SetErrorZ2(t.GetErrorZ2()*3);
1759 if (t.fCurrentCluster->IsUsed(10)){
1764 if (t.fNShared>0.7*t.GetNumberOfClusters()) {
1770 if (accept<3) UpdateTrack(&t,accept);
1773 if ( fIteration==0 && t.GetNFoundable()*0.5 > t.GetNumberOfClusters()) t.SetRemoval(10);
1781 //_________________________________________________________________________
1782 Bool_t AliTPCtrackerMI::GetTrackPoint(Int_t index, AliTrackPoint &p ) const
1784 // Get track space point by index
1785 // return false in case the cluster doesn't exist
1786 AliTPCclusterMI *cl = GetClusterMI(index);
1787 if (!cl) return kFALSE;
1788 Int_t sector = (index&0xff000000)>>24;
1789 // Int_t row = (index&0x00ff0000)>>16;
1791 // xyz[0] = fkParam->GetPadRowRadii(sector,row);
1792 xyz[0] = cl->GetX();
1793 xyz[1] = cl->GetY();
1794 xyz[2] = cl->GetZ();
1796 fkParam->AdjustCosSin(sector,cos,sin);
1797 Float_t x = cos*xyz[0]-sin*xyz[1];
1798 Float_t y = cos*xyz[1]+sin*xyz[0];
1800 Float_t sigmaY2 = 0.027*cl->GetSigmaY2();
1801 if (sector < fkParam->GetNInnerSector()) sigmaY2 *= 2.07;
1802 Float_t sigmaZ2 = 0.066*cl->GetSigmaZ2();
1803 if (sector < fkParam->GetNInnerSector()) sigmaZ2 *= 1.77;
1804 cov[0] = sin*sin*sigmaY2;
1805 cov[1] = -sin*cos*sigmaY2;
1807 cov[3] = cos*cos*sigmaY2;
1810 p.SetXYZ(x,y,xyz[2],cov);
1811 AliGeomManager::ELayerID iLayer;
1813 if (sector < fkParam->GetNInnerSector()) {
1814 iLayer = AliGeomManager::kTPC1;
1818 iLayer = AliGeomManager::kTPC2;
1819 idet = sector - fkParam->GetNInnerSector();
1821 UShort_t volid = AliGeomManager::LayerToVolUID(iLayer,idet);
1822 p.SetVolumeID(volid);
1828 Int_t AliTPCtrackerMI::UpdateClusters(AliTPCseed& t, Int_t nr) {
1829 //-----------------------------------------------------------------
1830 // This function tries to find a track prolongation to next pad row
1831 //-----------------------------------------------------------------
1832 t.SetCurrentCluster(0);
1833 t.SetCurrentClusterIndex1(-3);
1835 Double_t xt=t.GetX();
1836 Int_t row = GetRowNumber(xt)-1;
1837 Double_t ymax= GetMaxY(nr);
1839 if (row < nr) return 1; // don't prolongate if not information until now -
1840 // if (TMath::Abs(t.GetSnp())>0.9 && t.GetNumberOfClusters()>40. && fIteration!=2) {
1842 // return 0; // not prolongate strongly inclined tracks
1844 // if (TMath::Abs(t.GetSnp())>0.95) {
1846 // return 0; // not prolongate strongly inclined tracks
1847 // }// patch 28 fev 06
1849 Double_t x= GetXrow(nr);
1851 //t.PropagateTo(x+0.02);
1852 //t.PropagateTo(x+0.01);
1853 if (!t.PropagateTo(x)){
1860 if (TMath::Abs(y)>ymax){
1862 t.SetRelativeSector((t.GetRelativeSector()+1) % fN);
1863 if (!t.Rotate(fSectors->GetAlpha()))
1865 } else if (y <-ymax) {
1866 t.SetRelativeSector((t.GetRelativeSector()-1+fN) % fN);
1867 if (!t.Rotate(-fSectors->GetAlpha()))
1870 // if (!t.PropagateTo(x)){
1877 if (TMath::Abs(t.GetSnp())>AliTPCReconstructor::GetMaxSnpTracker()) return 0;
1879 if (!IsActive(t.GetRelativeSector(),nr)) {
1881 t.SetClusterIndex2(nr,-1);
1884 //AliInfo(Form("A - Sector%d phi %f - alpha %f", t.fRelativeSector,y/x, t.GetAlpha()));
1886 AliTPCtrackerRow &krow=GetRow(t.GetRelativeSector(),nr);
1888 if (TMath::Abs(TMath::Abs(y)-ymax)<krow.GetDeadZone()){
1890 t.SetClusterIndex2(nr,-1);
1896 // if (TMath::Abs(t.GetZ())<(AliTPCReconstructor::GetCtgRange()*t.GetX()+10) && (TMath::Abs(t.GetSnp())<AliTPCReconstructor::GetMaxSnpTracker()))
1897 if (IsFindable(t)) t.SetNFoundable(t.GetNFoundable()+1);
1903 if ( (nr%2==0) || t.GetNumberOfClusters()<2){
1904 // t.fCurrentSigmaY = GetSigmaY(&t);
1905 //t.fCurrentSigmaZ = GetSigmaZ(&t);
1909 AliTPCclusterMI *cl=0;
1912 Double_t roady = 1.;
1913 Double_t roadz = 1.;
1917 index = t.GetClusterIndex2(nr);
1918 if ( (index >= 0) && (index&0x8000)==0){
1919 cl = t.GetClusterPointer(nr);
1920 if ( (cl==0) && (index >= 0)) cl = GetClusterMI(index);
1921 t.SetCurrentClusterIndex1(index);
1923 t.SetCurrentCluster(cl);
1929 // if (index<0) return 0;
1930 UInt_t uindex = TMath::Abs(index);
1933 //cl = krow.FindNearest2(y+10,z,roady,roadz,uindex);
1934 cl = krow.FindNearest2(y,z,roady,roadz,uindex);
1937 if (cl) t.SetCurrentClusterIndex1(krow.GetIndex(uindex));
1938 t.SetCurrentCluster(cl);
1944 Int_t AliTPCtrackerMI::FollowToNextCluster(AliTPCseed & t, Int_t nr) {
1945 //-----------------------------------------------------------------
1946 // This function tries to find a track prolongation to next pad row
1947 //-----------------------------------------------------------------
1949 //update error according neighborhoud
1951 if (t.GetCurrentCluster()) {
1953 Int_t accept = AcceptCluster(&t,t.GetCurrentCluster());
1955 if (t.GetCurrentCluster()->IsUsed(10)){
1960 t.SetNShared(t.GetNShared()+1);
1961 if (t.GetNShared()>0.7*t.GetNumberOfClusters()) {
1966 if (fIteration>0) accept = 0;
1967 if (accept<3) UpdateTrack(&t,accept);
1971 if ( t.GetNumberOfClusters()>18 && ( (t.GetSigmaY2()+t.GetSigmaZ2())>0.16)) t.SetRemoval(10);
1972 if ( t.GetNumberOfClusters()>18 && t.GetChi2()/t.GetNumberOfClusters()>6 ) t.SetRemoval(10);
1974 if (( (t.GetNFoundable()*0.5 > t.GetNumberOfClusters()) || t.GetNoCluster()>15)) t.SetRemoval(10);
1982 //_____________________________________________________________________________
1983 Int_t AliTPCtrackerMI::FollowProlongation(AliTPCseed& t, Int_t rf, Int_t step, Bool_t fromSeeds) {
1984 //-----------------------------------------------------------------
1985 // This function tries to find a track prolongation.
1986 //-----------------------------------------------------------------
1987 Double_t xt=t.GetX();
1989 Double_t alpha=t.GetAlpha();
1990 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
1991 if (alpha < 0. ) alpha += 2.*TMath::Pi();
1993 t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN);
1995 Int_t first = GetRowNumber(xt);
2000 for (Int_t nr= first; nr>=rf; nr-=step) {
2002 if (t.GetKinkIndexes()[0]>0){
2003 for (Int_t i=0;i<3;i++){
2004 Int_t index = t.GetKinkIndexes()[i];
2005 if (index==0) break;
2006 if (index<0) continue;
2008 AliKink * kink = (AliKink*)fEvent->GetKink(index-1);
2010 printf("PROBLEM\n");
2013 Int_t kinkrow = kink->GetTPCRow0()+2+Int_t(0.5/(0.05+kink->GetAngle(2)));
2015 AliExternalTrackParam paramd(t);
2016 kink->SetDaughter(paramd);
2017 kink->SetStatus(2,5);
2024 if (nr==80) t.UpdateReference();
2025 if (nr<fInnerSec->GetNRows())
2026 fSectors = fInnerSec;
2028 fSectors = fOuterSec;
2029 if (FollowToNext(t,nr)==0)
2042 Int_t AliTPCtrackerMI::FollowBackProlongation(AliTPCseed& t, Int_t rf, Bool_t fromSeeds) {
2043 //-----------------------------------------------------------------
2044 // This function tries to find a track prolongation.
2045 //-----------------------------------------------------------------
2047 Double_t xt=t.GetX();
2048 Double_t alpha=t.GetAlpha();
2049 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
2050 if (alpha < 0. ) alpha += 2.*TMath::Pi();
2051 t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN);
2053 Int_t first = t.GetFirstPoint();
2054 Int_t ri = GetRowNumber(xt);
2058 if (first<ri) first = ri;
2060 if (first<0) first=0;
2061 for (Int_t nr=first; nr<=rf; nr++) {
2062 // if ( (TMath::Abs(t.GetSnp())>0.95)) break;//patch 28 fev 06
2063 if (t.GetKinkIndexes()[0]<0){
2064 for (Int_t i=0;i<3;i++){
2065 Int_t index = t.GetKinkIndexes()[i];
2066 if (index==0) break;
2067 if (index>0) continue;
2068 index = TMath::Abs(index);
2069 AliKink * kink = (AliKink*)fEvent->GetKink(index-1);
2071 printf("PROBLEM\n");
2074 Int_t kinkrow = kink->GetTPCRow0()-2-Int_t(0.5/(0.05+kink->GetAngle(2)));
2076 AliExternalTrackParam paramm(t);
2077 kink->SetMother(paramm);
2078 kink->SetStatus(2,1);
2085 if (nr<fInnerSec->GetNRows())
2086 fSectors = fInnerSec;
2088 fSectors = fOuterSec;
2099 Float_t AliTPCtrackerMI::OverlapFactor(AliTPCseed * s1, AliTPCseed * s2, Int_t &sum1, Int_t & sum2)
2101 // overlapping factor
2107 Float_t dz2 =(s1->GetZ() - s2->GetZ());
2110 Float_t dy2 =TMath::Abs((s1->GetY() - s2->GetY()));
2112 Float_t distance = TMath::Sqrt(dz2+dy2);
2113 if (distance>4.) return 0; // if there are far away - not overlap - to reduce combinatorics
2116 Int_t firstpoint = TMath::Min(s1->GetFirstPoint(),s2->GetFirstPoint());
2117 Int_t lastpoint = TMath::Max(s1->GetLastPoint(),s2->GetLastPoint());
2122 if (firstpoint>lastpoint) {
2123 firstpoint =lastpoint;
2128 for (Int_t i=firstpoint-1;i<lastpoint+1;i++){
2129 if (s1->GetClusterIndex2(i)>0) sum1++;
2130 if (s2->GetClusterIndex2(i)>0) sum2++;
2131 if (s1->GetClusterIndex2(i)==s2->GetClusterIndex2(i) && s1->GetClusterIndex2(i)>0) {
2135 if (sum<5) return 0;
2137 Float_t summin = TMath::Min(sum1+1,sum2+1);
2138 Float_t ratio = (sum+1)/Float_t(summin);
2142 void AliTPCtrackerMI::SignShared(AliTPCseed * s1, AliTPCseed * s2)
2146 Float_t thetaCut = 0.2;//+10.*TMath::Sqrt(s1->GetSigmaTglZ()+ s2->GetSigmaTglZ());
2147 if (TMath::Abs(s1->GetTgl()-s2->GetTgl())>thetaCut) return;
2148 Float_t minCl = TMath::Min(s1->GetNumberOfClusters(),s2->GetNumberOfClusters());
2149 Int_t cutN0 = TMath::Max(5,TMath::Nint(0.1*minCl));
2154 //Int_t firstpoint = TMath::Max(s1->GetFirstPoint(),s2->GetFirstPoint());
2155 //Int_t lastpoint = TMath::Min(s1->GetLastPoint(),s2->GetLastPoint());
2156 Int_t firstpoint = 0;
2157 Int_t lastpoint = 160;
2159 // if (firstpoint>=lastpoint-5) return;;
2161 for (Int_t i=firstpoint;i<lastpoint;i++){
2162 // if ( (s1->GetClusterIndex2(i)&0xFFFF8FFF)==(s2->GetClusterIndex2(i)&0xFFFF8FFF) && s1->GetClusterIndex2(i)>0) {
2163 if ( (s1->GetClusterIndex2(i))==(s2->GetClusterIndex2(i)) && s1->GetClusterIndex2(i)>0) {
2167 if (sumshared>cutN0){
2170 for (Int_t i=firstpoint;i<lastpoint;i++){
2171 // if ( (s1->GetClusterIndex2(i)&0xFFFF8FFF)==(s2->GetClusterIndex2(i)&0xFFFF8FFF) && s1->GetClusterIndex2(i)>0) {
2172 if ( (s1->GetClusterIndex2(i))==(s2->GetClusterIndex2(i)) && s1->GetClusterIndex2(i)>0) {
2173 AliTPCTrackerPoint *p1 = s1->GetTrackPoint(i);
2174 AliTPCTrackerPoint *p2 = s2->GetTrackPoint(i);;
2175 if (s1->IsActive()&&s2->IsActive()){
2176 p1->SetShared(kTRUE);
2177 p2->SetShared(kTRUE);
2183 if (sumshared>cutN0){
2184 for (Int_t i=0;i<4;i++){
2185 if (s1->GetOverlapLabel(3*i)==0){
2186 s1->SetOverlapLabel(3*i, s2->GetLabel());
2187 s1->SetOverlapLabel(3*i+1,sumshared);
2188 s1->SetOverlapLabel(3*i+2,s2->GetUniqueID());
2192 for (Int_t i=0;i<4;i++){
2193 if (s2->GetOverlapLabel(3*i)==0){
2194 s2->SetOverlapLabel(3*i, s1->GetLabel());
2195 s2->SetOverlapLabel(3*i+1,sumshared);
2196 s2->SetOverlapLabel(3*i+2,s1->GetUniqueID());
2203 void AliTPCtrackerMI::SignShared(TObjArray * arr)
2206 //sort trackss according sectors
2208 for (Int_t i=0; i<arr->GetEntriesFast(); i++) {
2209 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2211 //if (pt) RotateToLocal(pt);
2215 arr->Sort(); // sorting according relative sectors
2216 arr->Expand(arr->GetEntries());
2219 Int_t nseed=arr->GetEntriesFast();
2220 for (Int_t i=0; i<nseed; i++) {
2221 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2223 for (Int_t j=0;j<12;j++){
2224 pt->SetOverlapLabel(j,0);
2227 for (Int_t i=0; i<nseed; i++) {
2228 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2230 if (pt->GetRemoval()>10) continue;
2231 for (Int_t j=i+1; j<nseed; j++){
2232 AliTPCseed *pt2=(AliTPCseed*)arr->UncheckedAt(j);
2233 if (TMath::Abs(pt->GetRelativeSector()-pt2->GetRelativeSector())>1) continue;
2235 if (pt2->GetRemoval()<=10) {
2236 //if ( TMath::Abs(pt->GetRelativeSector()-pt2->GetRelativeSector())>0) break;
2244 void AliTPCtrackerMI::SortTracks(TObjArray * arr, Int_t mode) const
2247 //sort tracks in array according mode criteria
2248 Int_t nseed = arr->GetEntriesFast();
2249 for (Int_t i=0; i<nseed; i++) {
2250 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2261 void AliTPCtrackerMI::RemoveUsed2(TObjArray * arr, Float_t factor1, Float_t factor2, Int_t minimal)
2264 // Loop over all tracks and remove overlaped tracks (with lower quality)
2266 // 1. Unsign clusters
2267 // 2. Sort tracks according quality
2268 // Quality is defined by the number of cluster between first and last points
2270 // 3. Loop over tracks - decreasing quality order
2271 // a.) remove - If the fraction of shared cluster less than factor (1- n or 2)
2272 // b.) remove - If the minimal number of clusters less than minimal and not ITS
2273 // c.) if track accepted - sign clusters
2275 //Called in - AliTPCtrackerMI::Clusters2Tracks()
2276 // - AliTPCtrackerMI::PropagateBack()
2277 // - AliTPCtrackerMI::RefitInward()
2280 // factor1 - factor for constrained
2281 // factor2 - for non constrained tracks
2282 // if (Float_t(shared+1)/Float_t(found+1)>factor) - DELETE
2286 Int_t nseed = arr->GetEntriesFast();
2287 Float_t * quality = new Float_t[nseed];
2288 Int_t * indexes = new Int_t[nseed];
2292 for (Int_t i=0; i<nseed; i++) {
2293 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2298 pt->UpdatePoints(); //select first last max dens points
2299 Float_t * points = pt->GetPoints();
2300 if (points[3]<0.8) quality[i] =-1;
2301 quality[i] = (points[2]-points[0])+pt->GetNumberOfClusters();
2302 //prefer high momenta tracks if overlaps
2303 quality[i] *= TMath::Sqrt(TMath::Abs(pt->Pt())+0.5);
2305 TMath::Sort(nseed,quality,indexes);
2308 for (Int_t itrack=0; itrack<nseed; itrack++) {
2309 Int_t trackindex = indexes[itrack];
2310 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(trackindex);
2313 if (quality[trackindex]<0){
2314 MarkSeedFree( arr->RemoveAt(trackindex) );
2319 Int_t first = Int_t(pt->GetPoints()[0]);
2320 Int_t last = Int_t(pt->GetPoints()[2]);
2321 Double_t factor = (pt->GetBConstrain()) ? factor1: factor2;
2323 Int_t found,foundable,shared;
2324 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
2325 // pt->GetClusterStatistic(0,160, found, foundable,shared,kFALSE);
2326 Bool_t itsgold =kFALSE;
2329 if (pt->GetESD()->GetITSclusters(dummy)>4) itsgold= kTRUE;
2333 if (Float_t(shared+1)/Float_t(found+1)>factor){
2334 if (pt->GetKinkIndexes()[0]!=0) continue; //don't remove tracks - part of the kinks
2335 if( AliTPCReconstructor::StreamLevel()>3){
2336 TTreeSRedirector &cstream = *fDebugStreamer;
2337 cstream<<"RemoveUsed"<<
2338 "iter="<<fIteration<<
2342 MarkSeedFree( arr->RemoveAt(trackindex) );
2345 if (pt->GetNumberOfClusters()<50&&(found-0.5*shared)<minimal){ //remove short tracks
2346 if (pt->GetKinkIndexes()[0]!=0) continue; //don't remove tracks - part of the kinks
2347 if( AliTPCReconstructor::StreamLevel()>3){
2348 TTreeSRedirector &cstream = *fDebugStreamer;
2349 cstream<<"RemoveShort"<<
2350 "iter="<<fIteration<<
2354 MarkSeedFree( arr->RemoveAt(trackindex) );
2360 //if (sharedfactor>0.4) continue;
2361 if (pt->GetKinkIndexes()[0]>0) continue;
2362 //Remove tracks with undefined properties - seems
2363 if (pt->GetSigmaY2()<kAlmost0) continue; // ? what is the origin ?
2365 for (Int_t i=first; i<last; i++) {
2366 Int_t index=pt->GetClusterIndex2(i);
2367 // if (index<0 || index&0x8000 ) continue;
2368 if (index<0 || index&0x8000 ) continue;
2369 AliTPCclusterMI *c= pt->GetClusterPointer(i);
2376 Info("RemoveUsed","\n*****\nNumber of good tracks after shared removal\t%d\n",fNtracks);
2382 void AliTPCtrackerMI::DumpClusters(Int_t iter, TObjArray *trackArray)
2385 // Dump clusters after reco
2386 // signed and unsigned cluster can be visualized
2387 // 1. Unsign all cluster
2388 // 2. Sign all used clusters
2391 Int_t nseed = trackArray->GetEntries();
2392 for (Int_t i=0; i<nseed; i++){
2393 AliTPCseed *pt=(AliTPCseed*)trackArray->UncheckedAt(i);
2397 Bool_t isKink=pt->GetKinkIndex(0)!=0;
2398 for (Int_t j=0; j<160; ++j) {
2399 Int_t index=pt->GetClusterIndex2(j);
2400 if (index<0) continue;
2401 AliTPCclusterMI *c= pt->GetClusterPointer(j);
2403 if (isKink) c->Use(100); // kink
2404 c->Use(10); // by default usage 10
2409 for (Int_t sec=0;sec<fkNIS;sec++){
2410 for (Int_t row=0;row<fInnerSec->GetNRows();row++){
2411 TClonesArray *cla = fInnerSec[sec][row].GetClusters1();
2412 for (Int_t icl =0;icl< fInnerSec[sec][row].GetN1();icl++){
2413 AliTPCclusterMI* cl = (AliTPCclusterMI*)cla->At(icl);
2414 Float_t gx[3]; cl->GetGlobalXYZ(gx);
2415 (*fDebugStreamer)<<"clDump"<<
2423 cla = fInnerSec[sec][row].GetClusters2();
2424 for (Int_t icl =0;icl< fInnerSec[sec][row].GetN2();icl++){
2425 AliTPCclusterMI* cl = (AliTPCclusterMI*)cla->At(icl);
2426 Float_t gx[3]; cl->GetGlobalXYZ(gx);
2427 (*fDebugStreamer)<<"clDump"<<
2438 for (Int_t sec=0;sec<fkNOS;sec++){
2439 for (Int_t row=0;row<fOuterSec->GetNRows();row++){
2440 TClonesArray *cla = fOuterSec[sec][row].GetClusters1();
2441 for (Int_t icl =0;icl< fOuterSec[sec][row].GetN1();icl++){
2443 AliTPCclusterMI* cl = (AliTPCclusterMI*) cla->At(icl);
2444 cl->GetGlobalXYZ(gx);
2445 (*fDebugStreamer)<<"clDump"<<
2453 cla = fOuterSec[sec][row].GetClusters2();
2454 for (Int_t icl =0;icl< fOuterSec[sec][row].GetN2();icl++){
2456 AliTPCclusterMI* cl = (AliTPCclusterMI*) cla->At(icl);
2457 cl->GetGlobalXYZ(gx);
2458 (*fDebugStreamer)<<"clDump"<<
2470 void AliTPCtrackerMI::UnsignClusters()
2473 // loop over all clusters and unsign them
2476 for (Int_t sec=0;sec<fkNIS;sec++){
2477 for (Int_t row=0;row<fInnerSec->GetNRows();row++){
2478 TClonesArray *cla = fInnerSec[sec][row].GetClusters1();
2479 for (Int_t icl =0;icl< fInnerSec[sec][row].GetN1();icl++)
2480 // if (cl[icl].IsUsed(10))
2481 ((AliTPCclusterMI*) cla->At(icl))->Use(-1);
2482 cla = fInnerSec[sec][row].GetClusters2();
2483 for (Int_t icl =0;icl< fInnerSec[sec][row].GetN2();icl++)
2484 //if (cl[icl].IsUsed(10))
2485 ((AliTPCclusterMI*) cla->At(icl))->Use(-1);
2489 for (Int_t sec=0;sec<fkNOS;sec++){
2490 for (Int_t row=0;row<fOuterSec->GetNRows();row++){
2491 TClonesArray *cla = fOuterSec[sec][row].GetClusters1();
2492 for (Int_t icl =0;icl< fOuterSec[sec][row].GetN1();icl++)
2493 //if (cl[icl].IsUsed(10))
2494 ((AliTPCclusterMI*) cla->At(icl))->Use(-1);
2495 cla = fOuterSec[sec][row].GetClusters2();
2496 for (Int_t icl =0;icl< fOuterSec[sec][row].GetN2();icl++)
2497 //if (cl[icl].IsUsed(10))
2498 ((AliTPCclusterMI*) cla->At(icl))->Use(-1);
2506 void AliTPCtrackerMI::SignClusters(const TObjArray * arr, Float_t fnumber, Float_t fdensity)
2509 //sign clusters to be "used"
2511 // snumber and sdensity sign number of sigmas - bellow mean value to be accepted
2512 // loop over "primaries"
2526 Int_t nseed = arr->GetEntriesFast();
2527 for (Int_t i=0; i<nseed; i++) {
2528 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2532 if (!(pt->IsActive())) continue;
2533 Float_t dens = pt->GetNumberOfClusters()/Float_t(pt->GetNFoundable());
2534 if ( (dens>0.7) && (pt->GetNumberOfClusters()>70)){
2536 sumdens2+= dens*dens;
2537 sumn += pt->GetNumberOfClusters();
2538 sumn2 += pt->GetNumberOfClusters()*pt->GetNumberOfClusters();
2539 Float_t chi2 = pt->GetChi2()/pt->GetNumberOfClusters();
2542 sumchi2 +=chi2*chi2;
2547 Float_t mdensity = 0.9;
2548 Float_t meann = 130;
2549 Float_t meanchi = 1;
2550 Float_t sdensity = 0.1;
2551 Float_t smeann = 10;
2552 Float_t smeanchi =0.4;
2556 mdensity = sumdens/sum;
2558 meanchi = sumchi/sum;
2560 sdensity = sumdens2/sum-mdensity*mdensity;
2562 sdensity = TMath::Sqrt(sdensity);
2566 smeann = sumn2/sum-meann*meann;
2568 smeann = TMath::Sqrt(smeann);
2572 smeanchi = sumchi2/sum - meanchi*meanchi;
2574 smeanchi = TMath::Sqrt(smeanchi);
2580 //REMOVE SHORT DELTAS or tracks going out of sensitive volume of TPC
2582 for (Int_t i=0; i<nseed; i++) {
2583 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2587 if (pt->GetBSigned()) continue;
2588 if (pt->GetBConstrain()) continue;
2589 //if (!(pt->IsActive())) continue;
2591 Int_t found,foundable,shared;
2592 pt->GetClusterStatistic(0,160,found, foundable,shared);
2593 if (shared/float(found)>0.3) {
2594 if (shared/float(found)>0.9 ){
2595 //MarkSeedFree( arr->RemoveAt(i) );
2600 Bool_t isok =kFALSE;
2601 if ( (pt->GetNShared()/pt->GetNumberOfClusters()<0.5) &&pt->GetNumberOfClusters()>60)
2603 if ((TMath::Abs(1/pt->GetC())<100.) && (pt->GetNShared()/pt->GetNumberOfClusters()<0.7))
2605 if (TMath::Abs(pt->GetZ()/pt->GetX())>1.1)
2607 if ( (TMath::Abs(pt->GetSnp()>0.7) && pt->GetD(0,0)>60.))
2611 for (Int_t j=0; j<160; ++j) {
2612 Int_t index=pt->GetClusterIndex2(j);
2613 if (index<0) continue;
2614 AliTPCclusterMI *c= pt->GetClusterPointer(j);
2616 //if (!(c->IsUsed(10))) c->Use();
2623 Double_t maxchi = meanchi+2.*smeanchi;
2625 for (Int_t i=0; i<nseed; i++) {
2626 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2630 //if (!(pt->IsActive())) continue;
2631 if (pt->GetBSigned()) continue;
2632 Double_t chi = pt->GetChi2()/pt->GetNumberOfClusters();
2633 if (chi>maxchi) continue;
2636 Float_t dens = pt->GetNumberOfClusters()/Float_t(pt->GetNFoundable());
2638 //sign only tracks with enoug big density at the beginning
2640 if ((pt->GetDensityFirst(40)<0.75) && pt->GetNumberOfClusters()<meann) continue;
2643 Double_t mindens = TMath::Max(double(mdensity-sdensity*fdensity*bfactor),0.65);
2644 Double_t minn = TMath::Max(Int_t(meann-fnumber*smeann*bfactor),50);
2646 // if (pt->fBConstrain) mindens = TMath::Max(mdensity-sdensity*fdensity*bfactor,0.65);
2647 if ( (pt->GetRemoval()==10) && (pt->GetSnp()>0.8)&&(dens>mindens))
2650 if ((dens>mindens && pt->GetNumberOfClusters()>minn) && chi<maxchi ){
2651 //Int_t noc=pt->GetNumberOfClusters();
2652 pt->SetBSigned(kTRUE);
2653 for (Int_t j=0; j<160; ++j) {
2655 Int_t index=pt->GetClusterIndex2(j);
2656 if (index<0) continue;
2657 AliTPCclusterMI *c= pt->GetClusterPointer(j);
2659 // if (!(c->IsUsed(10))) c->Use();
2664 // gLastCheck = nseed;
2673 Int_t AliTPCtrackerMI::RefitInward(AliESDEvent *event)
2676 // back propagation of ESD tracks
2679 if (!event) return 0;
2680 const Int_t kMaxFriendTracks=2000;
2682 // extract correction object for multiplicity dependence of dEdx
2683 TObjArray * gainCalibArray = AliTPCcalibDB::Instance()->GetTimeGainSplinesRun(event->GetRunNumber());
2685 AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ;
2687 AliFatal("Tranformations not in RefitInward");
2690 transform->SetCurrentRecoParam((AliTPCRecoParam*)AliTPCReconstructor::GetRecoParam());
2691 const AliTPCRecoParam * recoParam = AliTPCcalibDB::Instance()->GetTransform()->GetCurrentRecoParam();
2692 Int_t nContribut = event->GetNumberOfTracks();
2693 TGraphErrors * graphMultDependenceDeDx = 0x0;
2694 if (recoParam && recoParam->GetUseMultiplicityCorrectionDedx() && gainCalibArray) {
2695 if (recoParam->GetUseTotCharge()) {
2696 graphMultDependenceDeDx = (TGraphErrors *) gainCalibArray->FindObject("TGRAPHERRORS_MEANQTOT_MULTIPLICITYDEPENDENCE_BEAM_ALL");
2698 graphMultDependenceDeDx = (TGraphErrors *) gainCalibArray->FindObject("TGRAPHERRORS_MEANQMAX_MULTIPLICITYDEPENDENCE_BEAM_ALL");
2704 //PrepareForProlongation(fSeeds,1);
2705 PropagateForward2(fSeeds);
2706 RemoveUsed2(fSeeds,0.4,0.4,20);
2708 TObjArray arraySeed(fSeeds->GetEntries());
2709 for (Int_t i=0;i<fSeeds->GetEntries();i++) {
2710 arraySeed.AddAt(fSeeds->At(i),i);
2712 SignShared(&arraySeed);
2713 // FindCurling(fSeeds, event,2); // find multi found tracks
2714 FindSplitted(fSeeds, event,2); // find multi found tracks
2715 if (AliTPCReconstructor::StreamLevel()>5) FindMultiMC(fSeeds, fEvent,2); // find multi found tracks
2718 Int_t nseed = fSeeds->GetEntriesFast();
2719 for (Int_t i=0;i<nseed;i++){
2720 AliTPCseed * seed = (AliTPCseed*) fSeeds->UncheckedAt(i);
2721 if (!seed) continue;
2722 if (seed->GetKinkIndex(0)>0) UpdateKinkQualityD(seed); // update quality informations for kinks
2723 AliESDtrack *esd=event->GetTrack(i);
2725 if (seed->GetNumberOfClusters()<60 && seed->GetNumberOfClusters()<(esd->GetTPCclusters(0) -5)*0.8){
2726 AliExternalTrackParam paramIn;
2727 AliExternalTrackParam paramOut;
2728 Int_t ncl = seed->RefitTrack(seed,¶mIn,¶mOut);
2729 if (AliTPCReconstructor::StreamLevel()>2) {
2730 (*fDebugStreamer)<<"RecoverIn"<<
2734 "pout.="<<¶mOut<<
2739 seed->Set(paramIn.GetX(),paramIn.GetAlpha(),paramIn.GetParameter(),paramIn.GetCovariance());
2740 seed->SetNumberOfClusters(ncl);
2744 seed->PropagateTo(fkParam->GetInnerRadiusLow());
2745 seed->UpdatePoints();
2746 AddCovariance(seed);
2747 MakeESDBitmaps(seed, esd);
2748 seed->CookdEdx(0.02,0.6);
2749 CookLabel(seed,0.1); //For comparison only
2751 if (AliTPCReconstructor::StreamLevel()>1 && seed!=0) {
2752 TTreeSRedirector &cstream = *fDebugStreamer;
2759 if (seed->GetNumberOfClusters()>15){
2760 esd->UpdateTrackParams(seed,AliESDtrack::kTPCrefit);
2761 esd->SetTPCPoints(seed->GetPoints());
2762 esd->SetTPCPointsF(seed->GetNFoundable());
2763 Int_t ndedx = seed->GetNCDEDX(0);
2764 Float_t sdedx = seed->GetSDEDX(0);
2765 Float_t dedx = seed->GetdEdx();
2766 // apply mutliplicity dependent dEdx correction if available
2767 if (graphMultDependenceDeDx) {
2768 Double_t corrGain = AliTPCcalibDButil::EvalGraphConst(graphMultDependenceDeDx, nContribut);
2769 dedx += (1 - corrGain)*50.; // MIP is normalized to 50
2771 esd->SetTPCsignal(dedx, sdedx, ndedx);
2773 // fill new dEdx information
2775 Double32_t signal[4];
2779 for(Int_t iarr=0;iarr<3;iarr++) {
2780 signal[iarr] = seed->GetDEDXregion(iarr+1);
2781 ncl[iarr] = seed->GetNCDEDX(iarr+1);
2782 nrows[iarr] = seed->GetNCDEDXInclThres(iarr+1);
2784 signal[3] = seed->GetDEDXregion(4);
2786 AliTPCdEdxInfo * infoTpcPid = new AliTPCdEdxInfo();
2787 infoTpcPid->SetTPCSignalRegionInfo(signal, ncl, nrows);
2788 esd->SetTPCdEdxInfo(infoTpcPid);
2790 // add seed to the esd track in Calib level
2792 Bool_t storeFriend = gRandom->Rndm()<(kMaxFriendTracks)/Float_t(nseed);
2793 if (AliTPCReconstructor::StreamLevel()>0 &&storeFriend){
2794 // RS: this is the only place where the seed is created not in the pool,
2795 // since it should belong to ESDevent
2796 AliTPCseed * seedCopy = new AliTPCseed(*seed, kTRUE);
2797 esd->AddCalibObject(seedCopy);
2802 //printf("problem\n");
2805 //FindKinks(fSeeds,event);
2806 if (AliTPCReconstructor::StreamLevel()>3) DumpClusters(2,fSeeds);
2807 Info("RefitInward","Number of refitted tracks %d",ntracks);
2809 AliCosmicTracker::FindCosmic(event, kTRUE);
2815 Int_t AliTPCtrackerMI::PropagateBack(AliESDEvent *event)
2818 // back propagation of ESD tracks
2820 if (!event) return 0;
2824 PropagateBack(fSeeds);
2825 RemoveUsed2(fSeeds,0.4,0.4,20);
2826 //FindCurling(fSeeds, fEvent,1);
2827 FindSplitted(fSeeds, event,1); // find multi found tracks
2828 if (AliTPCReconstructor::StreamLevel()>5) FindMultiMC(fSeeds, fEvent,1); // find multi found tracks
2831 Int_t nseed = fSeeds->GetEntriesFast();
2833 for (Int_t i=0;i<nseed;i++){
2834 AliTPCseed * seed = (AliTPCseed*) fSeeds->UncheckedAt(i);
2835 if (!seed) continue;
2836 if (seed->GetKinkIndex(0)<0) UpdateKinkQualityM(seed); // update quality informations for kinks
2837 seed->UpdatePoints();
2838 AddCovariance(seed);
2839 AliESDtrack *esd=event->GetTrack(i);
2840 if (!esd) continue; //never happen
2841 if (seed->GetNumberOfClusters()<60 && seed->GetNumberOfClusters()<(esd->GetTPCclusters(0) -5)*0.8){
2842 AliExternalTrackParam paramIn;
2843 AliExternalTrackParam paramOut;
2844 Int_t ncl = seed->RefitTrack(seed,¶mIn,¶mOut);
2845 if (AliTPCReconstructor::StreamLevel()>2) {
2846 (*fDebugStreamer)<<"RecoverBack"<<
2850 "pout.="<<¶mOut<<
2855 seed->Set(paramOut.GetX(),paramOut.GetAlpha(),paramOut.GetParameter(),paramOut.GetCovariance());
2856 seed->SetNumberOfClusters(ncl);
2859 seed->CookdEdx(0.02,0.6);
2860 CookLabel(seed,0.1); //For comparison only
2861 if (seed->GetNumberOfClusters()>15){
2862 esd->UpdateTrackParams(seed,AliESDtrack::kTPCout);
2863 esd->SetTPCPoints(seed->GetPoints());
2864 esd->SetTPCPointsF(seed->GetNFoundable());
2865 Int_t ndedx = seed->GetNCDEDX(0);
2866 Float_t sdedx = seed->GetSDEDX(0);
2867 Float_t dedx = seed->GetdEdx();
2868 esd->SetTPCsignal(dedx, sdedx, ndedx);
2870 Int_t eventnumber = event->GetEventNumberInFile();// patch 28 fev 06
2871 // This is most likely NOT the event number you'd like to use. It has nothing to do with the 'real' event number
2872 if (AliTPCReconstructor::StreamLevel()>1 && esd) {
2873 (*fDebugStreamer)<<"Cback"<<
2876 "EventNrInFile="<<eventnumber<<
2881 if (AliTPCReconstructor::StreamLevel()>3) DumpClusters(1,fSeeds);
2882 //FindKinks(fSeeds,event);
2883 Info("PropagateBack","Number of back propagated tracks %d",ntracks);
2890 Int_t AliTPCtrackerMI::PostProcess(AliESDEvent *event)
2893 // Post process events
2895 if (!event) return 0;
2898 // Set TPC event status
2901 // event affected by HV dip
2903 if(IsTPCHVDipEvent(event)) {
2904 event->ResetDetectorStatus(AliDAQ::kTPC);
2907 //printf("Status %d \n", event->IsDetectorOn(AliDAQ::kTPC));
2913 void AliTPCtrackerMI::DeleteSeeds()
2922 void AliTPCtrackerMI::ReadSeeds(const AliESDEvent *const event, Int_t direction)
2925 //read seeds from the event
2927 Int_t nentr=event->GetNumberOfTracks();
2929 Info("ReadSeeds", "Number of ESD tracks: %d\n", nentr);
2934 fSeeds = new TObjArray(nentr);
2938 for (Int_t i=0; i<nentr; i++) {
2939 AliESDtrack *esd=event->GetTrack(i);
2940 ULong_t status=esd->GetStatus();
2941 if (!(status&AliESDtrack::kTPCin)) continue;
2942 AliTPCtrack t(*esd);
2943 t.SetNumberOfClusters(0);
2944 // AliTPCseed *seed = new AliTPCseed(t,t.GetAlpha());
2945 AliTPCseed *seed = new( NextFreeSeed() ) AliTPCseed(t/*,t.GetAlpha()*/);
2946 seed->SetPoolID(fLastSeedID);
2947 seed->SetUniqueID(esd->GetID());
2948 AddCovariance(seed); //add systematic ucertainty
2949 for (Int_t ikink=0;ikink<3;ikink++) {
2950 Int_t index = esd->GetKinkIndex(ikink);
2951 seed->GetKinkIndexes()[ikink] = index;
2952 if (index==0) continue;
2953 index = TMath::Abs(index);
2954 AliESDkink * kink = fEvent->GetKink(index-1);
2955 if (kink&&esd->GetKinkIndex(ikink)<0){
2956 if ((status & AliESDtrack::kTRDrefit) != 0) kink->SetStatus(1,2);
2957 if ((status & AliESDtrack::kITSout) != 0) kink->SetStatus(1,0);
2959 if (kink&&esd->GetKinkIndex(ikink)>0){
2960 if ((status & AliESDtrack::kTRDrefit) != 0) kink->SetStatus(1,6);
2961 if ((status & AliESDtrack::kITSout) != 0) kink->SetStatus(1,4);
2965 if (((status&AliESDtrack::kITSout)==0)&&(direction==1)) seed->ResetCovariance(10.);
2966 if ( direction ==2 &&(status & AliESDtrack::kTRDrefit) == 0 ) seed->ResetCovariance(10.);
2967 //if ( direction ==2 && ((status & AliESDtrack::kTPCout) == 0) ) {
2968 // fSeeds->AddAt(0,i);
2969 // MarkSeedFree( seed );
2972 if ( direction ==2 &&(status & AliESDtrack::kTRDrefit) > 0 ) {
2973 Double_t par0[5],par1[5],alpha,x;
2974 esd->GetInnerExternalParameters(alpha,x,par0);
2975 esd->GetExternalParameters(x,par1);
2976 Double_t delta1 = TMath::Abs(par0[4]-par1[4])/(0.000000001+TMath::Abs(par0[4]+par1[4]));
2977 Double_t delta2 = TMath::Abs(par0[3]-par1[3]);
2979 if (esd->GetTRDncls()>0) trdchi2 = esd->GetTRDchi2()/esd->GetTRDncls();
2980 //reset covariance if suspicious
2981 if ( (delta1>0.1) || (delta2>0.006) ||trdchi2>7.)
2982 seed->ResetCovariance(10.);
2987 // rotate to the local coordinate system
2989 fSectors=fInnerSec; fN=fkNIS;
2990 Double_t alpha=seed->GetAlpha();
2991 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
2992 if (alpha < 0. ) alpha += 2.*TMath::Pi();
2993 Int_t ns=Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN;
2994 alpha =ns*fSectors->GetAlpha() + fSectors->GetAlphaShift();
2995 alpha-=seed->GetAlpha();
2996 if (alpha<-TMath::Pi()) alpha += 2*TMath::Pi();
2997 if (alpha>=TMath::Pi()) alpha -= 2*TMath::Pi();
2998 if (TMath::Abs(alpha) > 0.001) { // This should not happen normally
2999 AliWarning(Form("Rotating track over %f",alpha));
3000 if (!seed->Rotate(alpha)) {
3001 MarkSeedFree( seed );
3007 if (esd->GetKinkIndex(0)<=0){
3008 for (Int_t irow=0;irow<160;irow++){
3009 Int_t index = seed->GetClusterIndex2(irow);
3012 AliTPCclusterMI * cl = GetClusterMI(index);
3013 seed->SetClusterPointer(irow,cl);
3015 if ((index & 0x8000)==0){
3016 cl->Use(10); // accepted cluster
3018 cl->Use(6); // close cluster not accepted
3021 Info("ReadSeeds","Not found cluster");
3026 fSeeds->AddAt(seed,i);
3032 //_____________________________________________________________________________
3033 void AliTPCtrackerMI::MakeSeeds3(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, Float_t cuts[4],
3034 Float_t deltay, Int_t ddsec) {
3035 //-----------------------------------------------------------------
3036 // This function creates track seeds.
3037 // SEEDING WITH VERTEX CONSTRAIN
3038 //-----------------------------------------------------------------
3039 // cuts[0] - fP4 cut
3040 // cuts[1] - tan(phi) cut
3041 // cuts[2] - zvertex cut
3042 // cuts[3] - fP3 cut
3050 Double_t x[5], c[15];
3051 // Int_t di = i1-i2;
3053 AliTPCseed * seed = new( NextFreeSeed() ) AliTPCseed();
3054 seed->SetPoolID(fLastSeedID);
3055 Double_t alpha=fSectors->GetAlpha(), shift=fSectors->GetAlphaShift();
3056 Double_t cs=cos(alpha), sn=sin(alpha);
3058 // Double_t x1 =fOuterSec->GetX(i1);
3059 //Double_t xx2=fOuterSec->GetX(i2);
3061 Double_t x1 =GetXrow(i1);
3062 Double_t xx2=GetXrow(i2);
3064 Double_t x3=GetX(), y3=GetY(), z3=GetZ();
3066 Int_t imiddle = (i2+i1)/2; //middle pad row index
3067 Double_t xm = GetXrow(imiddle); // radius of middle pad-row
3068 const AliTPCtrackerRow& krm=GetRow(sec,imiddle); //middle pad -row
3072 const AliTPCtrackerRow& kr1=GetRow(ns,i1);
3073 Double_t ymax = GetMaxY(i1)-kr1.GetDeadZone()-1.5;
3074 Double_t ymaxm = GetMaxY(imiddle)-kr1.GetDeadZone()-1.5;
3077 // change cut on curvature if it can't reach this layer
3078 // maximal curvature set to reach it
3079 Double_t dvertexmax = TMath::Sqrt((x1-x3)*(x1-x3)+(ymax+5-y3)*(ymax+5-y3));
3080 if (dvertexmax*0.5*cuts[0]>0.85){
3081 cuts[0] = 0.85/(dvertexmax*0.5+1.);
3083 Double_t r2min = 1/(cuts[0]*cuts[0]); //minimal square of radius given by cut
3086 if (deltay>0) ddsec = 0;
3087 // loop over clusters
3088 for (Int_t is=0; is < kr1; is++) {
3090 if (kr1[is]->IsUsed(10)) continue;
3091 Double_t y1=kr1[is]->GetY(), z1=kr1[is]->GetZ();
3092 //if (TMath::Abs(y1)>ymax) continue;
3094 if (deltay>0 && TMath::Abs(ymax-TMath::Abs(y1))> deltay ) continue; // seed only at the edge
3096 // find possible directions
3097 Float_t anglez = (z1-z3)/(x1-x3);
3098 Float_t extraz = z1 - anglez*(x1-xx2); // extrapolated z
3101 //find rotation angles relative to line given by vertex and point 1
3102 Double_t dvertex2 = (x1-x3)*(x1-x3)+(y1-y3)*(y1-y3);
3103 Double_t dvertex = TMath::Sqrt(dvertex2);
3104 Double_t angle13 = TMath::ATan((y1-y3)/(x1-x3));
3105 Double_t cs13 = cos(-angle13), sn13 = sin(-angle13);
3108 // loop over 2 sectors
3114 Double_t dddz1=0; // direction of delta inclination in z axis
3121 for (Int_t dsec = dsec1; dsec<=dsec2;dsec++){
3122 Int_t sec2 = sec + dsec;
3124 // AliTPCtrackerRow& kr2 = fOuterSec[(sec2+fkNOS)%fkNOS][i2];
3125 //AliTPCtrackerRow& kr2m = fOuterSec[(sec2+fkNOS)%fkNOS][imiddle];
3126 AliTPCtrackerRow& kr2 = GetRow((sec2+fkNOS)%fkNOS,i2);
3127 AliTPCtrackerRow& kr2m = GetRow((sec2+fkNOS)%fkNOS,imiddle);
3128 Int_t index1 = TMath::Max(kr2.Find(extraz-0.6-dddz1*TMath::Abs(z1)*0.05)-1,0);
3129 Int_t index2 = TMath::Min(kr2.Find(extraz+0.6+dddz2*TMath::Abs(z1)*0.05)+1,kr2);
3131 // rotation angles to p1-p3
3132 Double_t cs13r = cos(-angle13+dsec*alpha)/dvertex, sn13r = sin(-angle13+dsec*alpha)/dvertex;
3133 Double_t x2, y2, z2;
3135 // Double_t dymax = maxangle*TMath::Abs(x1-xx2);
3138 Double_t dxx0 = (xx2-x3)*cs13r;
3139 Double_t dyy0 = (xx2-x3)*sn13r;
3140 for (Int_t js=index1; js < index2; js++) {
3141 const AliTPCclusterMI *kcl = kr2[js];
3142 if (kcl->IsUsed(10)) continue;
3144 //calcutate parameters
3146 Double_t yy0 = dyy0 +(kcl->GetY()-y3)*cs13r;
3148 if (TMath::Abs(yy0)<0.000001) continue;
3149 Double_t xx0 = dxx0 -(kcl->GetY()-y3)*sn13r;
3150 Double_t y0 = 0.5*(xx0*xx0+yy0*yy0-xx0)/yy0;
3151 Double_t r02 = (0.25+y0*y0)*dvertex2;
3152 //curvature (radius) cut
3153 if (r02<r2min) continue;
3157 Double_t c0 = 1/TMath::Sqrt(r02);
3161 //Double_t dfi0 = 2.*TMath::ASin(dvertex*c0*0.5);
3162 //Double_t dfi1 = 2.*TMath::ASin(TMath::Sqrt(yy0*yy0+(1-xx0)*(1-xx0))*dvertex*c0*0.5);
3163 Double_t dfi0 = 2.*AliTPCFastMath::FastAsin(dvertex*c0*0.5);
3164 Double_t dfi1 = 2.*AliTPCFastMath::FastAsin(TMath::Sqrt(yy0*yy0+(1-xx0)*(1-xx0))*dvertex*c0*0.5);
3167 Double_t z0 = kcl->GetZ();
3168 Double_t zzzz2 = z1-(z1-z3)*dfi1/dfi0;
3169 if (TMath::Abs(zzzz2-z0)>0.5) continue;
3172 Double_t dip = (z1-z0)*c0/dfi1;
3173 Double_t x0 = (0.5*cs13+y0*sn13)*dvertex*c0;
3184 x2= xx2*cs-y2*sn*dsec;
3185 y2=+xx2*sn*dsec+y2*cs;
3195 // do we have cluster at the middle ?
3197 GetProlongation(x1,xm,x,ym,zm);
3199 AliTPCclusterMI * cm=0;
3200 if (TMath::Abs(ym)-ymaxm<0){
3201 cm = krm.FindNearest2(ym,zm,1.0,0.6,dummy);
3202 if ((!cm) || (cm->IsUsed(10))) {
3207 // rotate y1 to system 0
3208 // get state vector in rotated system
3209 Double_t yr1 = (-0.5*sn13+y0*cs13)*dvertex*c0;
3210 Double_t xr2 = x0*cs+yr1*sn*dsec;
3211 Double_t xr[5]={kcl->GetY(),kcl->GetZ(), xr2, dip, c0};
3213 GetProlongation(xx2,xm,xr,ym,zm);
3214 if (TMath::Abs(ym)-ymaxm<0){
3215 cm = kr2m.FindNearest2(ym,zm,1.0,0.6,dummy);
3216 if ((!cm) || (cm->IsUsed(10))) {
3223 // Double_t dym = 0;
3224 // Double_t dzm = 0;
3226 // dym = ym - cm->GetY();
3227 // dzm = zm - cm->GetZ();
3234 Double_t sy1=kr1[is]->GetSigmaY2()*2., sz1=kr1[is]->GetSigmaZ2()*2.;
3235 Double_t sy2=kcl->GetSigmaY2()*2., sz2=kcl->GetSigmaZ2()*2.;
3236 //Double_t sy3=400*3./12., sy=0.1, sz=0.1;
3237 Double_t sy3=25000*x[4]*x[4]+0.1, sy=0.1, sz=0.1;
3238 //Double_t sy3=25000*x[4]*x[4]*60+0.5, sy=0.1, sz=0.1;
3240 Double_t f40=(F1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
3241 Double_t f42=(F1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
3242 Double_t f43=(F1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
3243 Double_t f20=(F2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
3244 Double_t f22=(F2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
3245 Double_t f23=(F2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
3247 Double_t f30=(F3(x1,y1+sy,x2,y2,z1,z2)-x[3])/sy;
3248 Double_t f31=(F3(x1,y1,x2,y2,z1+sz,z2)-x[3])/sz;
3249 Double_t f32=(F3(x1,y1,x2,y2+sy,z1,z2)-x[3])/sy;
3250 Double_t f34=(F3(x1,y1,x2,y2,z1,z2+sz)-x[3])/sz;
3254 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3255 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3256 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3257 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3258 c[13]=f30*sy1*f40+f32*sy2*f42;
3259 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3261 // if (!BuildSeed(kr1[is],kcl,0,x1,x2,x3,x,c)) continue;
3263 UInt_t index=kr1.GetIndex(is);
3264 if (seed) {MarkSeedFree(seed); seed = 0;}
3265 AliTPCseed *track = seed = new( NextFreeSeed() ) AliTPCseed(x1, ns*alpha+shift, x, c, index);
3266 seed->SetPoolID(fLastSeedID);
3267 track->SetIsSeeding(kTRUE);
3268 track->SetSeed1(i1);
3269 track->SetSeed2(i2);
3270 track->SetSeedType(3);
3274 FollowProlongation(*track, (i1+i2)/2,1);
3275 Int_t foundable,found,shared;
3276 track->GetClusterStatistic((i1+i2)/2,i1, found, foundable, shared, kTRUE);
3277 if ((found<0.55*foundable) || shared>0.5*found || (track->GetSigmaY2()+track->GetSigmaZ2())>0.5){
3278 MarkSeedFree(seed); seed = 0;
3284 FollowProlongation(*track, i2,1);
3288 track->SetBConstrain(1);
3289 // track->fLastPoint = i1+fInnerSec->GetNRows(); // first cluster in track position
3290 track->SetLastPoint(i1); // first cluster in track position
3291 track->SetFirstPoint(track->GetLastPoint());
3293 if (track->GetNumberOfClusters()<(i1-i2)*0.5 ||
3294 track->GetNumberOfClusters() < track->GetNFoundable()*0.6 ||
3295 track->GetNShared()>0.4*track->GetNumberOfClusters() ) {
3296 MarkSeedFree(seed); seed = 0;
3300 // Z VERTEX CONDITION
3301 Double_t zv, bz=GetBz();
3302 if ( !track->GetZAt(0.,bz,zv) ) continue;
3303 if (TMath::Abs(zv-z3)>cuts[2]) {
3304 FollowProlongation(*track, TMath::Max(i2-20,0));
3305 if ( !track->GetZAt(0.,bz,zv) ) continue;
3306 if (TMath::Abs(zv-z3)>cuts[2]){
3307 FollowProlongation(*track, TMath::Max(i2-40,0));
3308 if ( !track->GetZAt(0.,bz,zv) ) continue;
3309 if (TMath::Abs(zv-z3)>cuts[2] &&(track->GetNumberOfClusters() > track->GetNFoundable()*0.7)){
3310 // make seed without constrain
3311 AliTPCseed * track2 = MakeSeed(track,0.2,0.5,1.);
3312 FollowProlongation(*track2, i2,1);
3313 track2->SetBConstrain(kFALSE);
3314 track2->SetSeedType(1);
3315 arr->AddLast(track2);
3316 MarkSeedFree( seed ); seed = 0;
3320 MarkSeedFree( seed ); seed = 0;
3327 track->SetSeedType(0);
3328 arr->AddLast(track); // note, track is seed, don't free the seed
3329 seed = new( NextFreeSeed() ) AliTPCseed;
3330 seed->SetPoolID(fLastSeedID);
3332 // don't consider other combinations
3333 if (track->GetNumberOfClusters() > track->GetNFoundable()*0.8)
3339 Info("MakeSeeds3","\nSeeding statistic:\t%d\t%d\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin,nout1,nout2);
3341 if (seed) MarkSeedFree( seed );
3345 void AliTPCtrackerMI::MakeSeeds5(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, Float_t cuts[4],
3350 //-----------------------------------------------------------------
3351 // This function creates track seeds.
3352 //-----------------------------------------------------------------
3353 // cuts[0] - fP4 cut
3354 // cuts[1] - tan(phi) cut
3355 // cuts[2] - zvertex cut
3356 // cuts[3] - fP3 cut
3366 Double_t x[5], c[15];
3368 // make temporary seed
3369 AliTPCseed * seed = new( NextFreeSeed() ) AliTPCseed;
3370 seed->SetPoolID(fLastSeedID);
3371 Double_t alpha=fOuterSec->GetAlpha(), shift=fOuterSec->GetAlphaShift();
3372 // Double_t cs=cos(alpha), sn=sin(alpha);
3377 Double_t x1 = GetXrow(i1-1);
3378 const AliTPCtrackerRow& kr1=GetRow(sec,i1-1);
3379 Double_t y1max = GetMaxY(i1-1)-kr1.GetDeadZone()-1.5;
3381 Double_t x1p = GetXrow(i1);
3382 const AliTPCtrackerRow& kr1p=GetRow(sec,i1);
3384 Double_t x1m = GetXrow(i1-2);
3385 const AliTPCtrackerRow& kr1m=GetRow(sec,i1-2);
3388 //last 3 padrow for seeding
3389 AliTPCtrackerRow& kr3 = GetRow((sec+fkNOS)%fkNOS,i1-7);
3390 Double_t x3 = GetXrow(i1-7);
3391 // Double_t y3max= GetMaxY(i1-7)-kr3.fDeadZone-1.5;
3393 AliTPCtrackerRow& kr3p = GetRow((sec+fkNOS)%fkNOS,i1-6);
3394 Double_t x3p = GetXrow(i1-6);
3396 AliTPCtrackerRow& kr3m = GetRow((sec+fkNOS)%fkNOS,i1-8);
3397 Double_t x3m = GetXrow(i1-8);
3402 Int_t im = i1-4; //middle pad row index
3403 Double_t xm = GetXrow(im); // radius of middle pad-row
3404 const AliTPCtrackerRow& krm=GetRow(sec,im); //middle pad -row
3405 // Double_t ymmax = GetMaxY(im)-kr1.fDeadZone-1.5;
3408 Double_t deltax = x1-x3;
3409 Double_t dymax = deltax*cuts[1];
3410 Double_t dzmax = deltax*cuts[3];
3412 // loop over clusters
3413 for (Int_t is=0; is < kr1; is++) {
3415 if (kr1[is]->IsUsed(10)) continue;
3416 Double_t y1=kr1[is]->GetY(), z1=kr1[is]->GetZ();
3418 if (deltay>0 && TMath::Abs(y1max-TMath::Abs(y1))> deltay ) continue; // seed only at the edge
3420 Int_t index1 = TMath::Max(kr3.Find(z1-dzmax)-1,0);
3421 Int_t index2 = TMath::Min(kr3.Find(z1+dzmax)+1,kr3);
3427 for (Int_t js=index1; js < index2; js++) {
3428 const AliTPCclusterMI *kcl = kr3[js];
3429 if (kcl->IsUsed(10)) continue;
3431 // apply angular cuts
3432 if (TMath::Abs(y1-y3)>dymax) continue;
3435 if (TMath::Abs(z1-z3)>dzmax) continue;
3437 Double_t angley = (y1-y3)/(x1-x3);
3438 Double_t anglez = (z1-z3)/(x1-x3);
3440 Double_t erry = TMath::Abs(angley)*(x1-x1m)*0.5+0.5;
3441 Double_t errz = TMath::Abs(anglez)*(x1-x1m)*0.5+0.5;
3443 Double_t yyym = angley*(xm-x1)+y1;
3444 Double_t zzzm = anglez*(xm-x1)+z1;
3446 const AliTPCclusterMI *kcm = krm.FindNearest2(yyym,zzzm,erry,errz,index);
3448 if (kcm->IsUsed(10)) continue;
3450 erry = TMath::Abs(angley)*(x1-x1m)*0.4+0.5;
3451 errz = TMath::Abs(anglez)*(x1-x1m)*0.4+0.5;
3458 // look around first
3459 const AliTPCclusterMI *kc1m = kr1m.FindNearest2(angley*(x1m-x1)+y1,
3465 if (kc1m->IsUsed(10)) used++;
3467 const AliTPCclusterMI *kc1p = kr1p.FindNearest2(angley*(x1p-x1)+y1,
3473 if (kc1p->IsUsed(10)) used++;
3475 if (used>1) continue;
3476 if (found<1) continue;
3480 const AliTPCclusterMI *kc3m = kr3m.FindNearest2(angley*(x3m-x3)+y3,
3486 if (kc3m->IsUsed(10)) used++;
3490 const AliTPCclusterMI *kc3p = kr3p.FindNearest2(angley*(x3p-x3)+y3,
3496 if (kc3p->IsUsed(10)) used++;
3500 if (used>1) continue;
3501 if (found<3) continue;
3511 x[4]=F1(x1,y1,x2,y2,x3,y3);
3512 //if (TMath::Abs(x[4]) >= cuts[0]) continue;
3515 x[2]=F2(x1,y1,x2,y2,x3,y3);
3518 x[3]=F3n(x1,y1,x2,y2,z1,z2,x[4]);
3519 //if (TMath::Abs(x[3]) > cuts[3]) continue;
3523 Double_t sy1=0.1, sz1=0.1;
3524 Double_t sy2=0.1, sz2=0.1;
3525 Double_t sy3=0.1, sy=0.1, sz=0.1;
3527 Double_t f40=(F1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
3528 Double_t f42=(F1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
3529 Double_t f43=(F1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
3530 Double_t f20=(F2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
3531 Double_t f22=(F2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
3532 Double_t f23=(F2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
3534 Double_t f30=(F3(x1,y1+sy,x2,y2,z1,z2)-x[3])/sy;
3535 Double_t f31=(F3(x1,y1,x2,y2,z1+sz,z2)-x[3])/sz;
3536 Double_t f32=(F3(x1,y1,x2,y2+sy,z1,z2)-x[3])/sy;
3537 Double_t f34=(F3(x1,y1,x2,y2,z1,z2+sz)-x[3])/sz;
3541 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3542 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3543 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3544 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3545 c[13]=f30*sy1*f40+f32*sy2*f42;
3546 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3548 // if (!BuildSeed(kr1[is],kcl,0,x1,x2,x3,x,c)) continue;
3550 index=kr1.GetIndex(is);
3551 if (seed) {MarkSeedFree( seed ); seed = 0;}
3552 AliTPCseed *track = seed = new( NextFreeSeed() ) AliTPCseed(x1, sec*alpha+shift, x, c, index);
3553 seed->SetPoolID(fLastSeedID);
3555 track->SetIsSeeding(kTRUE);
3558 FollowProlongation(*track, i1-7,1);
3559 if (track->GetNumberOfClusters() < track->GetNFoundable()*0.75 ||
3560 track->GetNShared()>0.6*track->GetNumberOfClusters() || ( track->GetSigmaY2()+ track->GetSigmaZ2())>0.6){
3561 MarkSeedFree( seed ); seed = 0;
3567 FollowProlongation(*track, i2,1);
3568 track->SetBConstrain(0);
3569 track->SetLastPoint(i1+fInnerSec->GetNRows()); // first cluster in track position
3570 track->SetFirstPoint(track->GetLastPoint());
3572 if (track->GetNumberOfClusters()<(i1-i2)*0.5 ||
3573 track->GetNumberOfClusters()<track->GetNFoundable()*0.7 ||
3574 track->GetNShared()>2. || track->GetChi2()/track->GetNumberOfClusters()>6 || ( track->GetSigmaY2()+ track->GetSigmaZ2())>0.5 ) {
3575 MarkSeedFree( seed ); seed = 0;
3580 FollowProlongation(*track, TMath::Max(i2-10,0),1);
3581 AliTPCseed * track2 = MakeSeed(track,0.2,0.5,0.9);
3582 FollowProlongation(*track2, i2,1);
3583 track2->SetBConstrain(kFALSE);
3584 track2->SetSeedType(4);
3585 arr->AddLast(track2);
3586 MarkSeedFree( seed ); seed = 0;
3590 //arr->AddLast(track);
3591 //seed = new AliTPCseed;
3597 Info("MakeSeeds5","\nSeeding statiistic:\t%d\t%d\t%d\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin,nout1,nout2,nout3);
3599 if (seed) MarkSeedFree(seed);
3603 //_____________________________________________________________________________
3604 void AliTPCtrackerMI::MakeSeeds2(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, Float_t */*cuts[4]*/,
3605 Float_t deltay, Bool_t /*bconstrain*/) {
3606 //-----------------------------------------------------------------
3607 // This function creates track seeds - without vertex constraint
3608 //-----------------------------------------------------------------
3609 // cuts[0] - fP4 cut - not applied
3610 // cuts[1] - tan(phi) cut
3611 // cuts[2] - zvertex cut - not applied
3612 // cuts[3] - fP3 cut
3622 Double_t alpha=fOuterSec->GetAlpha(), shift=fOuterSec->GetAlphaShift();
3623 // Double_t cs=cos(alpha), sn=sin(alpha);
3624 Int_t row0 = (i1+i2)/2;
3625 Int_t drow = (i1-i2)/2;
3626 const AliTPCtrackerRow& kr0=fSectors[sec][row0];
3627 AliTPCtrackerRow * kr=0;
3629 AliTPCpolyTrack polytrack;
3630 Int_t nclusters=fSectors[sec][row0];
3631 AliTPCseed * seed = new( NextFreeSeed() ) AliTPCseed;
3632 seed->SetPoolID(fLastSeedID);
3637 for (Int_t is=0; is < nclusters; is++) { //LOOP over clusters
3639 Int_t nfoundable =0;
3640 for (Int_t iter =1; iter<2; iter++){ //iterations
3641 const AliTPCtrackerRow& krm=fSectors[sec][row0-iter];
3642 const AliTPCtrackerRow& krp=fSectors[sec][row0+iter];
3643 const AliTPCclusterMI * cl= kr0[is];
3645 if (cl->IsUsed(10)) {
3651 Double_t x = kr0.GetX();
3652 // Initialization of the polytrack
3657 Double_t y0= cl->GetY();
3658 Double_t z0= cl->GetZ();
3662 Double_t ymax = fSectors->GetMaxY(row0)-kr0.GetDeadZone()-1.5;
3663 if (deltay>0 && TMath::Abs(ymax-TMath::Abs(y0))> deltay ) continue; // seed only at the edge
3665 erry = (0.5)*cl->GetSigmaY2()/TMath::Sqrt(cl->GetQ())*6;
3666 errz = (0.5)*cl->GetSigmaZ2()/TMath::Sqrt(cl->GetQ())*6;
3667 polytrack.AddPoint(x,y0,z0,erry, errz);
3670 if (cl->IsUsed(10)) sumused++;
3673 Float_t roady = (5*TMath::Sqrt(cl->GetSigmaY2()+0.2)+1.)*iter;
3674 Float_t roadz = (5*TMath::Sqrt(cl->GetSigmaZ2()+0.2)+1.)*iter;
3677 AliTPCclusterMI * cl1 = krm.FindNearest(y0,z0,roady,roadz);
3678 if (cl1 && TMath::Abs(ymax-TMath::Abs(y0))) {
3679 erry = (0.5)*cl1->GetSigmaY2()/TMath::Sqrt(cl1->GetQ())*3;
3680 errz = (0.5)*cl1->GetSigmaZ2()/TMath::Sqrt(cl1->GetQ())*3;
3681 if (cl1->IsUsed(10)) sumused++;
3682 polytrack.AddPoint(x,cl1->GetY(),cl1->GetZ(),erry,errz);
3686 AliTPCclusterMI * cl2 = krp.FindNearest(y0,z0,roady,roadz);
3688 erry = (0.5)*cl2->GetSigmaY2()/TMath::Sqrt(cl2->GetQ())*3;
3689 errz = (0.5)*cl2->GetSigmaZ2()/TMath::Sqrt(cl2->GetQ())*3;
3690 if (cl2->IsUsed(10)) sumused++;
3691 polytrack.AddPoint(x,cl2->GetY(),cl2->GetZ(),erry,errz);
3694 if (sumused>0) continue;
3696 polytrack.UpdateParameters();
3702 nfoundable = polytrack.GetN();
3703 nfound = nfoundable;
3705 for (Int_t ddrow = iter+1; ddrow<drow;ddrow++){
3706 Float_t maxdist = 0.8*(1.+3./(ddrow));
3707 for (Int_t delta = -1;delta<=1;delta+=2){
3708 Int_t row = row0+ddrow*delta;
3709 kr = &(fSectors[sec][row]);
3710 Double_t xn = kr->GetX();
3711 Double_t ymax1 = fSectors->GetMaxY(row)-kr->GetDeadZone()-1.5;
3712 polytrack.GetFitPoint(xn,yn,zn);
3713 if (TMath::Abs(yn)>ymax1) continue;
3715 AliTPCclusterMI * cln = kr->FindNearest(yn,zn,roady,roadz);
3717 Float_t dist = TMath::Sqrt( (yn-cln->GetY())*(yn-cln->GetY())+(zn-cln->GetZ())*(zn-cln->GetZ()));
3720 erry = (dist+0.3)*cln->GetSigmaY2()/TMath::Sqrt(cln->GetQ())*(1.+1./(ddrow));
3721 errz = (dist+0.3)*cln->GetSigmaZ2()/TMath::Sqrt(cln->GetQ())*(1.+1./(ddrow));
3722 if (cln->IsUsed(10)) {
3723 // printf("used\n");
3731 polytrack.AddPoint(xn,cln->GetY(),cln->GetZ(),erry, errz);
3736 if ( (sumused>3) || (sumused>0.5*nfound) || (nfound<0.6*nfoundable)) break;
3737 polytrack.UpdateParameters();
3740 if ( (sumused>3) || (sumused>0.5*nfound)) {
3741 //printf("sumused %d\n",sumused);
3746 polytrack.GetFitDerivation(kr0.GetX(),dy,dz);
3747 AliTPCpolyTrack track2;
3749 polytrack.Refit(track2,0.5+TMath::Abs(dy)*0.3,0.4+TMath::Abs(dz)*0.3);
3750 if (track2.GetN()<0.5*nfoundable) continue;
3753 if ((nfound>0.6*nfoundable) &&( nfoundable>0.4*(i1-i2))) {
3755 // test seed with and without constrain
3756 for (Int_t constrain=0; constrain<=0;constrain++){
3757 // add polytrack candidate
3759 Double_t x[5], c[15];
3760 Double_t x1,x2,x3,y1,y2,y3,z1,z2,z3;
3761 track2.GetBoundaries(x3,x1);
3763 track2.GetFitPoint(x1,y1,z1);
3764 track2.GetFitPoint(x2,y2,z2);
3765 track2.GetFitPoint(x3,y3,z3);
3767 //is track pointing to the vertex ?
3770 polytrack.GetFitPoint(x0,y0,z0);
3783 x[4]=F1(x1,y1,x2,y2,x3,y3);
3785 // if (TMath::Abs(x[4]) >= cuts[0]) continue; //
3786 x[2]=F2(x1,y1,x2,y2,x3,y3);
3788 //if (TMath::Abs(x[4]*x1-x[2]) >= cuts[1]) continue;
3789 //x[3]=F3(x1,y1,x2,y2,z1,z2);
3790 x[3]=F3n(x1,y1,x3,y3,z1,z3,x[4]);
3791 //if (TMath::Abs(x[3]) > cuts[3]) continue;
3794 Double_t sy =0.1, sz =0.1;
3795 Double_t sy1=0.02, sz1=0.02;
3796 Double_t sy2=0.02, sz2=0.02;
3800 sy3=25000*x[4]*x[4]+0.1, sy=0.1, sz=0.1;
3803 Double_t f40=(F1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
3804 Double_t f42=(F1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
3805 Double_t f43=(F1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
3806 Double_t f20=(F2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
3807 Double_t f22=(F2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
3808 Double_t f23=(F2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
3810 Double_t f30=(F3(x1,y1+sy,x3,y3,z1,z3)-x[3])/sy;
3811 Double_t f31=(F3(x1,y1,x3,y3,z1+sz,z3)-x[3])/sz;
3812 Double_t f32=(F3(x1,y1,x3,y3+sy,z1,z3)-x[3])/sy;
3813 Double_t f34=(F3(x1,y1,x3,y3,z1,z3+sz)-x[3])/sz;
3818 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3819 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3820 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3821 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3822 c[13]=f30*sy1*f40+f32*sy2*f42;
3823 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3825 //Int_t row1 = fSectors->GetRowNumber(x1);
3826 Int_t row1 = GetRowNumber(x1);
3830 if (seed) {MarkSeedFree( seed ); seed = 0;}
3831 AliTPCseed *track = seed = new( NextFreeSeed() ) AliTPCseed(x1,sec*alpha+shift,x,c,index);
3832 seed->SetPoolID(fLastSeedID);
3833 track->SetIsSeeding(kTRUE);
3834 Int_t rc=FollowProlongation(*track, i2);
3835 if (constrain) track->SetBConstrain(1);
3837 track->SetBConstrain(0);
3838 track->SetLastPoint(row1+fInnerSec->GetNRows()); // first cluster in track position
3839 track->SetFirstPoint(track->GetLastPoint());
3841 if (rc==0 || track->GetNumberOfClusters()<(i1-i2)*0.5 ||
3842 track->GetNumberOfClusters() < track->GetNFoundable()*0.6 ||
3843 track->GetNShared()>0.4*track->GetNumberOfClusters()) {
3844 MarkSeedFree( seed ); seed = 0;
3847 arr->AddLast(track); // track IS seed, don't free seed
3848 seed = new( NextFreeSeed() ) AliTPCseed;
3849 seed->SetPoolID(fLastSeedID);
3853 } // if accepted seed
3856 Info("MakeSeeds2","\nSeeding statiistic:\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin3);
3858 if (seed) MarkSeedFree( seed );
3862 AliTPCseed *AliTPCtrackerMI::MakeSeed(AliTPCseed *const track, Float_t r0, Float_t r1, Float_t r2)
3866 //reseed using track points
3867 Int_t p0 = int(r0*track->GetNumberOfClusters()); // point 0
3868 Int_t p1 = int(r1*track->GetNumberOfClusters());
3869 Int_t p2 = int(r2*track->GetNumberOfClusters()); // last point
3871 Double_t x0[3],x1[3],x2[3];
3872 for (Int_t i=0;i<3;i++){
3878 // find track position at given ratio of the length
3879 Int_t sec0=0, sec1=0, sec2=0;
3882 for (Int_t i=0;i<160;i++){
3883 if (track->GetClusterPointer(i)){
3885 AliTPCTrackerPoint *trpoint =track->GetTrackPoint(i);
3886 if ( (index<p0) || x0[0]<0 ){
3887 if (trpoint->GetX()>1){
3888 clindex = track->GetClusterIndex2(i);
3890 x0[0] = trpoint->GetX();
3891 x0[1] = trpoint->GetY();
3892 x0[2] = trpoint->GetZ();
3893 sec0 = ((clindex&0xff000000)>>24)%18;
3898 if ( (index<p1) &&(trpoint->GetX()>1)){
3899 clindex = track->GetClusterIndex2(i);
3901 x1[0] = trpoint->GetX();
3902 x1[1] = trpoint->GetY();
3903 x1[2] = trpoint->GetZ();
3904 sec1 = ((clindex&0xff000000)>>24)%18;
3907 if ( (index<p2) &&(trpoint->GetX()>1)){
3908 clindex = track->GetClusterIndex2(i);
3910 x2[0] = trpoint->GetX();
3911 x2[1] = trpoint->GetY();
3912 x2[2] = trpoint->GetZ();
3913 sec2 = ((clindex&0xff000000)>>24)%18;
3920 Double_t alpha, cs,sn, xx2,yy2;
3922 alpha = (sec1-sec2)*fSectors->GetAlpha();
3923 cs = TMath::Cos(alpha);
3924 sn = TMath::Sin(alpha);
3925 xx2= x1[0]*cs-x1[1]*sn;
3926 yy2= x1[0]*sn+x1[1]*cs;
3930 alpha = (sec0-sec2)*fSectors->GetAlpha();
3931 cs = TMath::Cos(alpha);
3932 sn = TMath::Sin(alpha);
3933 xx2= x0[0]*cs-x0[1]*sn;
3934 yy2= x0[0]*sn+x0[1]*cs;
3940 Double_t x[5],c[15];
3944 x[4]=F1(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]);
3945 // if (x[4]>1) return 0;
3946 x[2]=F2(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]);
3947 x[3]=F3n(x2[0],x2[1],x0[0],x0[1],x2[2],x0[2],x[4]);
3948 //if (TMath::Abs(x[3]) > 2.2) return 0;
3949 //if (TMath::Abs(x[2]) > 1.99) return 0;
3951 Double_t sy =0.1, sz =0.1;
3953 Double_t sy1=0.02+track->GetSigmaY2(), sz1=0.02+track->GetSigmaZ2();
3954 Double_t sy2=0.01+track->GetSigmaY2(), sz2=0.01+track->GetSigmaZ2();
3955 Double_t sy3=0.01+track->GetSigmaY2();
3957 Double_t f40=(F1(x2[0],x2[1]+sy,x1[0],x1[1],x0[0],x0[1])-x[4])/sy;
3958 Double_t f42=(F1(x2[0],x2[1],x1[0],x1[1]+sy,x0[0],x0[1])-x[4])/sy;
3959 Double_t f43=(F1(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]+sy)-x[4])/sy;
3960 Double_t f20=(F2(x2[0],x2[1]+sy,x1[0],x1[1],x0[0],x0[1])-x[2])/sy;
3961 Double_t f22=(F2(x2[0],x2[1],x1[0],x1[1]+sy,x0[0],x0[1])-x[2])/sy;
3962 Double_t f23=(F2(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]+sy)-x[2])/sy;
3964 Double_t f30=(F3(x2[0],x2[1]+sy,x0[0],x0[1],x2[2],x0[2])-x[3])/sy;
3965 Double_t f31=(F3(x2[0],x2[1],x0[0],x0[1],x2[2]+sz,x0[2])-x[3])/sz;
3966 Double_t f32=(F3(x2[0],x2[1],x0[0],x0[1]+sy,x2[2],x0[2])-x[3])/sy;
3967 Double_t f34=(F3(x2[0],x2[1],x0[0],x0[1],x2[2],x0[2]+sz)-x[3])/sz;
3972 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3973 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3974 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3975 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3976 c[13]=f30*sy1*f40+f32*sy2*f42;
3977 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3979 // Int_t row1 = fSectors->GetRowNumber(x2[0]);
3980 AliTPCseed *seed = new( NextFreeSeed() ) AliTPCseed(x2[0], sec2*fSectors->GetAlpha()+fSectors->GetAlphaShift(), x, c, 0);
3981 seed->SetPoolID(fLastSeedID);
3982 // Double_t y0,z0,y1,z1, y2,z2;
3983 //seed->GetProlongation(x0[0],y0,z0);
3984 // seed->GetProlongation(x1[0],y1,z1);
3985 //seed->GetProlongation(x2[0],y2,z2);
3987 seed->SetLastPoint(pp2);
3988 seed->SetFirstPoint(pp2);
3995 AliTPCseed *AliTPCtrackerMI::ReSeed(const AliTPCseed *track, Float_t r0, Float_t r1, Float_t r2)
3999 //reseed using founded clusters
4001 // Find the number of clusters
4002 Int_t nclusters = 0;
4003 for (Int_t irow=0;irow<160;irow++){
4004 if (track->GetClusterIndex(irow)>0) nclusters++;
4008 ipos[0] = TMath::Max(int(r0*nclusters),0); // point 0 cluster
4009 ipos[1] = TMath::Min(int(r1*nclusters),nclusters-1); //
4010 ipos[2] = TMath::Min(int(r2*nclusters),nclusters-1); // last point
4013 Double_t xyz[3][3]={{0}};
4014 Int_t row[3]={0},sec[3]={0,0,0};
4016 // find track row position at given ratio of the length
4018 for (Int_t irow=0;irow<160;irow++){
4019 if (track->GetClusterIndex2(irow)<0) continue;
4021 for (Int_t ipoint=0;ipoint<3;ipoint++){
4022 if (index<=ipos[ipoint]) row[ipoint] = irow;
4026 //Get cluster and sector position
4027 for (Int_t ipoint=0;ipoint<3;ipoint++){
4028 Int_t clindex = track->GetClusterIndex2(row[ipoint]);
4029 AliTPCclusterMI * cl = GetClusterMI(clindex);
4032 // AliTPCclusterMI * cl = GetClusterMI(clindex);
4035 sec[ipoint] = ((clindex&0xff000000)>>24)%18;
4036 xyz[ipoint][0] = GetXrow(row[ipoint]);
4037 xyz[ipoint][1] = cl->GetY();
4038 xyz[ipoint][2] = cl->GetZ();
4042 // Calculate seed state vector and covariance matrix
4044 Double_t alpha, cs,sn, xx2,yy2;
4046 alpha = (sec[1]-sec[2])*fSectors->GetAlpha();
4047 cs = TMath::Cos(alpha);
4048 sn = TMath::Sin(alpha);
4049 xx2= xyz[1][0]*cs-xyz[1][1]*sn;
4050 yy2= xyz[1][0]*sn+xyz[1][1]*cs;
4054 alpha = (sec[0]-sec[2])*fSectors->GetAlpha();
4055 cs = TMath::Cos(alpha);
4056 sn = TMath::Sin(alpha);
4057 xx2= xyz[0][0]*cs-xyz[0][1]*sn;
4058 yy2= xyz[0][0]*sn+xyz[0][1]*cs;
4064 Double_t x[5],c[15];
4068 x[4]=F1(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
4069 x[2]=F2(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
4070 x[3]=F3n(xyz[2][0],xyz[2][1],xyz[0][0],xyz[0][1],xyz[2][2],xyz[0][2],x[4]);
4072 Double_t sy =0.1, sz =0.1;
4074 Double_t sy1=0.2, sz1=0.2;
4075 Double_t sy2=0.2, sz2=0.2;
4078 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;
4079 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;
4080 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;
4081 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;
4082 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;
4083 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;
4085 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;
4086 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;
4087 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;
4088 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;
4093 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
4094 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
4095 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
4096 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
4097 c[13]=f30*sy1*f40+f32*sy2*f42;
4098 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
4100 // Int_t row1 = fSectors->GetRowNumber(xyz[2][0]);
4101 AliTPCseed *seed=new( NextFreeSeed() ) AliTPCseed(xyz[2][0], sec[2]*fSectors->GetAlpha()+fSectors->GetAlphaShift(), x, c, 0);
4102 seed->SetPoolID(fLastSeedID);
4103 seed->SetLastPoint(row[2]);
4104 seed->SetFirstPoint(row[2]);
4109 AliTPCseed *AliTPCtrackerMI::ReSeed(AliTPCseed *track,Int_t r0, Bool_t forward)
4113 //reseed using founded clusters
4116 Int_t row[3]={0,0,0};
4117 Int_t sec[3]={0,0,0};
4119 // forward direction
4121 for (Int_t irow=r0;irow<160;irow++){
4122 if (track->GetClusterIndex(irow)>0){
4127 for (Int_t irow=160;irow>r0;irow--){
4128 if (track->GetClusterIndex(irow)>0){
4133 for (Int_t irow=row[2]-15;irow>row[0];irow--){
4134 if (track->GetClusterIndex(irow)>0){
4142 for (Int_t irow=0;irow<r0;irow++){
4143 if (track->GetClusterIndex(irow)>0){
4148 for (Int_t irow=r0;irow>0;irow--){
4149 if (track->GetClusterIndex(irow)>0){
4154 for (Int_t irow=row[2]-15;irow>row[0];irow--){
4155 if (track->GetClusterIndex(irow)>0){
4162 if ((row[2]-row[0])<20) return 0;
4163 if (row[1]==0) return 0;
4166 //Get cluster and sector position
4167 for (Int_t ipoint=0;ipoint<3;ipoint++){
4168 Int_t clindex = track->GetClusterIndex2(row[ipoint]);
4169 AliTPCclusterMI * cl = GetClusterMI(clindex);
4172 // AliTPCclusterMI * cl = GetClusterMI(clindex);
4175 sec[ipoint] = ((clindex&0xff000000)>>24)%18;
4176 xyz[ipoint][0] = GetXrow(row[ipoint]);
4177 AliTPCTrackerPoint * point = track->GetTrackPoint(row[ipoint]);
4178 if (point&&ipoint<2){
4180 xyz[ipoint][1] = point->GetY();
4181 xyz[ipoint][2] = point->GetZ();
4184 xyz[ipoint][1] = cl->GetY();
4185 xyz[ipoint][2] = cl->GetZ();
4192 // Calculate seed state vector and covariance matrix
4194 Double_t alpha, cs,sn, xx2,yy2;
4196 alpha = (sec[1]-sec[2])*fSectors->GetAlpha();
4197 cs = TMath::Cos(alpha);
4198 sn = TMath::Sin(alpha);
4199 xx2= xyz[1][0]*cs-xyz[1][1]*sn;
4200 yy2= xyz[1][0]*sn+xyz[1][1]*cs;
4204 alpha = (sec[0]-sec[2])*fSectors->GetAlpha();
4205 cs = TMath::Cos(alpha);
4206 sn = TMath::Sin(alpha);
4207 xx2= xyz[0][0]*cs-xyz[0][1]*sn;
4208 yy2= xyz[0][0]*sn+xyz[0][1]*cs;
4214 Double_t x[5],c[15];
4218 x[4]=F1(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
4219 x[2]=F2(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
4220 x[3]=F3n(xyz[2][0],xyz[2][1],xyz[0][0],xyz[0][1],xyz[2][2],xyz[0][2],x[4]);
4222 Double_t sy =0.1, sz =0.1;
4224 Double_t sy1=0.2, sz1=0.2;
4225 Double_t sy2=0.2, sz2=0.2;
4228 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;
4229 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;
4230 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;
4231 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;
4232 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;
4233 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;
4235 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;
4236 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;
4237 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;
4238 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;
4243 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
4244 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
4245 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
4246 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
4247 c[13]=f30*sy1*f40+f32*sy2*f42;
4248 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
4250 // Int_t row1 = fSectors->GetRowNumber(xyz[2][0]);
4251 AliTPCseed *seed=new( NextFreeSeed() ) AliTPCseed(xyz[2][0], sec[2]*fSectors->GetAlpha()+fSectors->GetAlphaShift(), x, c, 0);
4252 seed->SetPoolID(fLastSeedID);
4253 seed->SetLastPoint(row[2]);
4254 seed->SetFirstPoint(row[2]);
4255 for (Int_t i=row[0];i<row[2];i++){
4256 seed->SetClusterIndex(i, track->GetClusterIndex(i));
4264 void AliTPCtrackerMI::FindMultiMC(const TObjArray * array, AliESDEvent */*esd*/, Int_t iter)
4267 // find multi tracks - THIS FUNCTION IS ONLY FOR DEBUG PURPOSES
4269 // Use AliTPCReconstructor::StreamLevel()>2 if you want to tune parameters - cuts
4271 // Two reasons to have multiple find tracks
4272 // 1. Curling tracks can be find more than once
4273 // 2. Splitted tracks
4274 // a.) Multiple seeding to increase tracking efficiency - (~ 100% reached)
4275 // b.) Edge effect on the sector boundaries
4278 // Algorithm done in 2 phases - because of CPU consumption
4279 // it is n^2 algorithm - for lead-lead 20000x20000 combination are investigated
4281 // Algorihm for curling tracks sign:
4282 // 1 phase -makes a very rough fast cuts to minimize combinatorics
4283 // a.) opposite sign
4284 // b.) one of the tracks - not pointing to the primary vertex -
4285 // c.) delta tan(theta)
4287 // 2 phase - calculates DCA between tracks - time consument
4292 // General cuts - for splitted tracks and for curling tracks
4294 const Float_t kMaxdPhi = 0.2; // maximal distance in phi
4296 // Curling tracks cuts
4301 Int_t nentries = array->GetEntriesFast();
4302 AliHelix *helixes = new AliHelix[nentries];
4303 Float_t *xm = new Float_t[nentries];
4304 Float_t *dz0 = new Float_t[nentries];
4305 Float_t *dz1 = new Float_t[nentries];
4311 // Find track COG in x direction - point with best defined parameters
4313 for (Int_t i=0;i<nentries;i++){
4314 AliTPCseed* track = (AliTPCseed*)array->At(i);
4315 if (!track) continue;
4316 track->SetCircular(0);
4317 new (&helixes[i]) AliHelix(*track);
4321 track->GetDZ(GetX(),GetY(),GetZ(),GetBz(),dz);
4324 for (Int_t icl=0; icl<160; icl++){
4325 AliTPCclusterMI * cl = track->GetClusterPointer(icl);
4331 if (ncl>0) xm[i]/=Float_t(ncl);
4334 for (Int_t i0=0;i0<nentries;i0++){
4335 AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
4336 if (!track0) continue;
4337 Float_t xc0 = helixes[i0].GetHelix(6);
4338 Float_t yc0 = helixes[i0].GetHelix(7);
4339 Float_t r0 = helixes[i0].GetHelix(8);
4340 Float_t rc0 = TMath::Sqrt(xc0*xc0+yc0*yc0);
4341 Float_t fi0 = TMath::ATan2(yc0,xc0);
4343 for (Int_t i1=i0+1;i1<nentries;i1++){
4344 AliTPCseed * track1 = (AliTPCseed*)array->At(i1);
4345 if (!track1) continue;
4346 Int_t lab0=track0->GetLabel();
4347 Int_t lab1=track1->GetLabel();
4348 if (TMath::Abs(lab0)!=TMath::Abs(lab1)) continue;
4350 Float_t xc1 = helixes[i1].GetHelix(6);
4351 Float_t yc1 = helixes[i1].GetHelix(7);
4352 Float_t r1 = helixes[i1].GetHelix(8);
4353 Float_t rc1 = TMath::Sqrt(xc1*xc1+yc1*yc1);
4354 Float_t fi1 = TMath::ATan2(yc1,xc1);
4356 Float_t dfi = fi0-fi1;
4359 if (dfi>1.5*TMath::Pi()) dfi-=TMath::Pi(); // take care about edge effect
4360 if (dfi<-1.5*TMath::Pi()) dfi+=TMath::Pi(); //
4361 if (TMath::Abs(dfi)>kMaxdPhi&&helixes[i0].GetHelix(4)*helixes[i1].GetHelix(4)<0){
4363 // if short tracks with undefined sign
4364 fi1 = -TMath::ATan2(yc1,-xc1);
4367 Float_t dtheta = TMath::Abs(track0->GetTgl()-track1->GetTgl())<TMath::Abs(track0->GetTgl()+track1->GetTgl())? track0->GetTgl()-track1->GetTgl():track0->GetTgl()+track1->GetTgl();
4370 // debug stream to tune "fast cuts"
4372 Double_t dist[3]; // distance at X
4373 Double_t mdist[3]={0,0,0}; // mean distance X+-40cm
4374 track0->GetDistance(track1,0.5*(xm[i0]+xm[i1])-40.,dist,AliTracker::GetBz());
4375 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[i]);
4376 track0->GetDistance(track1,0.5*(xm[i0]+xm[i1])+40.,dist,AliTracker::GetBz());
4377 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[i]);
4378 track0->GetDistance(track1,0.5*(xm[i0]+xm[i1]),dist,AliTracker::GetBz());
4379 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[i]);
4380 for (Int_t i=0;i<3;i++) mdist[i]*=0.33333;
4384 for (Int_t icl=0; icl<160; icl++){
4385 AliTPCclusterMI * cl0 = track0->GetClusterPointer(icl);
4386 AliTPCclusterMI * cl1 = track1->GetClusterPointer(icl);
4389 if (cl0==cl1) sums++;
4393 if (AliTPCReconstructor::StreamLevel()>5) {
4394 TTreeSRedirector &cstream = *fDebugStreamer;
4399 "Tr0.="<<track0<< // seed0
4400 "Tr1.="<<track1<< // seed1
4401 "h0.="<<&helixes[i0]<<
4402 "h1.="<<&helixes[i1]<<
4404 "sum="<<sum<< //the sum of rows with cl in both
4405 "sums="<<sums<< //the sum of shared clusters
4406 "xm0="<<xm[i0]<< // the center of track
4407 "xm1="<<xm[i1]<< // the x center of track
4408 // General cut variables
4409 "dfi="<<dfi<< // distance in fi angle
4410 "dtheta="<<dtheta<< // distance int theta angle
4416 "dist0="<<dist[0]<< //distance x
4417 "dist1="<<dist[1]<< //distance y
4418 "dist2="<<dist[2]<< //distance z
4419 "mdist0="<<mdist[0]<< //distance x
4420 "mdist1="<<mdist[1]<< //distance y
4421 "mdist2="<<mdist[2]<< //distance z
4437 if (AliTPCReconstructor::StreamLevel()>1) {
4438 AliInfo("Time for curling tracks removal DEBUGGING MC");
4445 void AliTPCtrackerMI::FindSplitted(TObjArray * array, AliESDEvent */*esd*/, Int_t /*iter*/){
4447 // Find Splitted tracks and remove the one with worst quality
4448 // Corresponding debug streamer to tune selections - "Splitted2"
4450 // 0. Sort tracks according quility
4451 // 1. Propagate the tracks to the reference radius
4452 // 2. Double_t loop to select close tracks (only to speed up process)
4453 // 3. Calculate cluster overlap ratio - and remove the track if bigger than a threshold
4454 // 4. Delete temporary parameters
4456 const Double_t xref=GetXrow(63); // reference radius -IROC/OROC boundary
4458 const Double_t kCutP1=10; // delta Z cut 10 cm
4459 const Double_t kCutP2=0.15; // delta snp(fi) cut 0.15
4460 const Double_t kCutP3=0.15; // delta tgl(theta) cut 0.15
4461 const Double_t kCutAlpha=0.15; // delta alpha cut
4462 Int_t firstpoint = 0;
4463 Int_t lastpoint = 160;
4465 Int_t nentries = array->GetEntriesFast();
4466 AliExternalTrackParam *params = new AliExternalTrackParam[nentries];
4472 //0. Sort tracks according quality
4473 //1. Propagate the ext. param to reference radius
4474 Int_t nseed = array->GetEntriesFast();
4475 if (nseed<=0) return;
4476 Float_t * quality = new Float_t[nseed];
4477 Int_t * indexes = new Int_t[nseed];
4478 for (Int_t i=0; i<nseed; i++) {
4479 AliTPCseed *pt=(AliTPCseed*)array->UncheckedAt(i);
4484 pt->UpdatePoints(); //select first last max dens points
4485 Float_t * points = pt->GetPoints();
4486 if (points[3]<0.8) quality[i] =-1;
4487 quality[i] = (points[2]-points[0])+pt->GetNumberOfClusters();
4488 //prefer high momenta tracks if overlaps
4489 quality[i] *= TMath::Sqrt(TMath::Abs(pt->Pt())+0.5);
4491 AliTracker::PropagateTrackToBxByBz(&(params[i]),xref,pt->GetMass(),5.,kTRUE);
4492 AliTracker::PropagateTrackToBxByBz(&(params[i]),xref,pt->GetMass(),1.,kTRUE);
4494 TMath::Sort(nseed,quality,indexes);
4496 // 3. Loop over pair of tracks
4498 for (Int_t i0=0; i0<nseed; i0++) {
4499 Int_t index0=indexes[i0];
4500 if (!(array->UncheckedAt(index0))) continue;
4501 AliTPCseed *s1 = (AliTPCseed*)array->UncheckedAt(index0);
4502 if (!s1->IsActive()) continue;
4503 AliExternalTrackParam &par0=params[index0];
4504 for (Int_t i1=i0+1; i1<nseed; i1++) {
4505 Int_t index1=indexes[i1];
4506 if (!(array->UncheckedAt(index1))) continue;
4507 AliTPCseed *s2 = (AliTPCseed*)array->UncheckedAt(index1);
4508 if (!s2->IsActive()) continue;
4509 if (s2->GetKinkIndexes()[0]!=0)
4510 if (s2->GetKinkIndexes()[0] == -s1->GetKinkIndexes()[0]) continue;
4511 AliExternalTrackParam &par1=params[index1];
4512 if (TMath::Abs(par0.GetParameter()[3]-par1.GetParameter()[3])>kCutP3) continue;
4513 if (TMath::Abs(par0.GetParameter()[1]-par1.GetParameter()[1])>kCutP1) continue;
4514 if (TMath::Abs(par0.GetParameter()[2]-par1.GetParameter()[2])>kCutP2) continue;
4515 Double_t dAlpha= TMath::Abs(par0.GetAlpha()-par1.GetAlpha());
4516 if (dAlpha>TMath::Pi()) dAlpha-=TMath::Pi();
4517 if (TMath::Abs(dAlpha)>kCutAlpha) continue;
4522 Int_t firstShared=lastpoint, lastShared=firstpoint;
4523 Int_t firstRow=lastpoint, lastRow=firstpoint;
4525 for (Int_t i=firstpoint;i<lastpoint;i++){
4526 if (s1->GetClusterIndex2(i)>0) nall0++;
4527 if (s2->GetClusterIndex2(i)>0) nall1++;
4528 if (s1->GetClusterIndex2(i)>0 && s2->GetClusterIndex2(i)>0) {
4529 if (i<firstRow) firstRow=i;
4530 if (i>lastRow) lastRow=i;
4532 if ( (s1->GetClusterIndex2(i))==(s2->GetClusterIndex2(i)) && s1->GetClusterIndex2(i)>0) {
4533 if (i<firstShared) firstShared=i;
4534 if (i>lastShared) lastShared=i;
4538 Double_t ratio0 = Float_t(sumShared)/Float_t(TMath::Min(nall0+1,nall1+1));
4539 Double_t ratio1 = Float_t(sumShared)/Float_t(TMath::Max(nall0+1,nall1+1));
4541 if( AliTPCReconstructor::StreamLevel()>1){
4542 TTreeSRedirector &cstream = *fDebugStreamer;
4543 Int_t n0=s1->GetNumberOfClusters();
4544 Int_t n1=s2->GetNumberOfClusters();
4545 Int_t n0F=s1->GetNFoundable();
4546 Int_t n1F=s2->GetNFoundable();
4547 Int_t lab0=s1->GetLabel();
4548 Int_t lab1=s2->GetLabel();
4550 cstream<<"Splitted2"<<
4551 "iter="<<fIteration<<
4552 "lab0="<<lab0<< // MC label if exist
4553 "lab1="<<lab1<< // MC label if exist
4556 "ratio0="<<ratio0<< // shared ratio
4557 "ratio1="<<ratio1<< // shared ratio
4558 "p0.="<<&par0<< // track parameters
4560 "s0.="<<s1<< // full seed
4562 "n0="<<n0<< // number of clusters track 0
4563 "n1="<<n1<< // number of clusters track 1
4564 "nall0="<<nall0<< // number of clusters track 0
4565 "nall1="<<nall1<< // number of clusters track 1
4566 "n0F="<<n0F<< // number of findable
4567 "n1F="<<n1F<< // number of findable
4568 "shared="<<sumShared<< // number of shared clusters
4569 "firstS="<<firstShared<< // first and the last shared row
4570 "lastS="<<lastShared<<
4571 "firstRow="<<firstRow<< // first and the last row with cluster
4572 "lastRow="<<lastRow<< //
4576 // remove track with lower quality
4578 if (ratio0>AliTPCReconstructor::GetRecoParam()->GetCutSharedClusters(0) ||
4579 ratio1>AliTPCReconstructor::GetRecoParam()->GetCutSharedClusters(1)){
4583 MarkSeedFree( array->RemoveAt(index1) );
4588 // 4. Delete temporary array
4598 void AliTPCtrackerMI::FindCurling(const TObjArray * array, AliESDEvent */*esd*/, Int_t iter)
4601 // find Curling tracks
4602 // Use AliTPCReconstructor::StreamLevel()>1 if you want to tune parameters - cuts
4605 // Algorithm done in 2 phases - because of CPU consumption
4606 // it is n^2 algorithm - for lead-lead 20000x20000 combination are investigated
4607 // see detal in MC part what can be used to cut
4611 const Float_t kMaxC = 400; // maximal curvature to of the track
4612 const Float_t kMaxdTheta = 0.15; // maximal distance in theta
4613 const Float_t kMaxdPhi = 0.15; // maximal distance in phi
4614 const Float_t kPtRatio = 0.3; // ratio between pt
4615 const Float_t kMinDCAR = 2.; // distance to the primary vertex in r - see cpipe cut
4618 // Curling tracks cuts
4621 const Float_t kMaxDeltaRMax = 40; // distance in outer radius
4622 const Float_t kMaxDeltaRMin = 5.; // distance in lower radius - see cpipe cut
4623 const Float_t kMinAngle = 2.9; // angle between tracks
4624 const Float_t kMaxDist = 5; // biggest distance
4626 // The cuts can be tuned using the "MC information stored in Multi tree ==> see FindMultiMC
4629 TCut csign("csign","Tr0.fP[4]*Tr1.fP[4]<0"); //opposite sign
4630 TCut cmax("cmax","abs(Tr0.GetC())>1/400");
4631 TCut cda("cda","sqrt(dtheta^2+dfi^2)<0.15");
4632 TCut ccratio("ccratio","abs((Tr0.fP[4]+Tr1.fP[4])/(abs(Tr0.fP[4])+abs(Tr1.fP[4])))<0.3");
4633 TCut cpipe("cpipe", "min(abs(r0-rc0),abs(r1-rc1))>5");
4635 TCut cdrmax("cdrmax","abs(abs(rc0+r0)-abs(rc1+r1))<40")
4636 TCut cdrmin("cdrmin","abs(abs(rc0+r0)-abs(rc1+r1))<10")
4638 Multi->Draw("dfi","iter==0"+csign+cmax+cda+ccratio); ~94% of curling tracks fulfill
4639 Multi->Draw("min(abs(r0-rc0),abs(r1-rc1))","iter==0&&abs(lab1)==abs(lab0)"+csign+cmax+cda+ccratio+cpipe+cdrmin+cdrmax); //80%
4641 Curling2->Draw("dfi","iter==0&&abs(lab0)==abs(lab1)"+csign+cmax+cdtheta+cdfi+ccratio)
4647 Int_t nentries = array->GetEntriesFast();
4648 AliHelix *helixes = new AliHelix[nentries];
4649 for (Int_t i=0;i<nentries;i++){
4650 AliTPCseed* track = (AliTPCseed*)array->At(i);
4651 if (!track) continue;
4652 track->SetCircular(0);
4653 new (&helixes[i]) AliHelix(*track);
4659 Double_t phase[2][2]={{0,0},{0,0}},radius[2]={0,0};
4665 for (Int_t i0=0;i0<nentries;i0++){
4666 AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
4667 if (!track0) continue;
4668 if (TMath::Abs(track0->GetC())<1/kMaxC) continue;
4669 Float_t xc0 = helixes[i0].GetHelix(6);
4670 Float_t yc0 = helixes[i0].GetHelix(7);
4671 Float_t r0 = helixes[i0].GetHelix(8);
4672 Float_t rc0 = TMath::Sqrt(xc0*xc0+yc0*yc0);
4673 Float_t fi0 = TMath::ATan2(yc0,xc0);
4675 for (Int_t i1=i0+1;i1<nentries;i1++){
4676 AliTPCseed * track1 = (AliTPCseed*)array->At(i1);
4677 if (!track1) continue;
4678 if (TMath::Abs(track1->GetC())<1/kMaxC) continue;
4679 Float_t xc1 = helixes[i1].GetHelix(6);
4680 Float_t yc1 = helixes[i1].GetHelix(7);
4681 Float_t r1 = helixes[i1].GetHelix(8);
4682 Float_t rc1 = TMath::Sqrt(xc1*xc1+yc1*yc1);
4683 Float_t fi1 = TMath::ATan2(yc1,xc1);
4685 Float_t dfi = fi0-fi1;
4688 if (dfi>1.5*TMath::Pi()) dfi-=TMath::Pi(); // take care about edge effect
4689 if (dfi<-1.5*TMath::Pi()) dfi+=TMath::Pi(); //
4690 Float_t dtheta = TMath::Abs(track0->GetTgl()-track1->GetTgl())<TMath::Abs(track0->GetTgl()+track1->GetTgl())? track0->GetTgl()-track1->GetTgl():track0->GetTgl()+track1->GetTgl();
4694 if (track0->GetBConstrain()&&track1->GetBConstrain()) continue; // not constrained
4695 if (track1->GetSigned1Pt()*track0->GetSigned1Pt()>0) continue; // not the same sign
4696 if ( TMath::Abs(track1->GetTgl()+track0->GetTgl())>kMaxdTheta) continue; //distance in the Theta
4697 if ( TMath::Abs(dfi)>kMaxdPhi) continue; //distance in phi
4698 if ( TMath::Sqrt(dfi*dfi+dtheta*dtheta)>kMaxdPhi) continue; //common angular offset
4700 Float_t pt0 = track0->GetSignedPt();
4701 Float_t pt1 = track1->GetSignedPt();
4702 if ((TMath::Abs(pt0+pt1)/(TMath::Abs(pt0)+TMath::Abs(pt1)))>kPtRatio) continue;
4703 if ((iter==1) && TMath::Abs(TMath::Abs(rc0+r0)-TMath::Abs(rc1+r1))>kMaxDeltaRMax) continue;
4704 if ((iter!=1) &&TMath::Abs(TMath::Abs(rc0-r0)-TMath::Abs(rc1-r1))>kMaxDeltaRMin) continue;
4705 if (TMath::Min(TMath::Abs(rc0-r0),TMath::Abs(rc1-r1))<kMinDCAR) continue;
4708 // Now find closest approach
4712 Int_t npoints = helixes[i0].GetRPHIintersections(helixes[i1], phase, radius,10);
4713 if (npoints==0) continue;
4714 helixes[i0].GetClosestPhases(helixes[i1], phase);
4718 Double_t hangles[3];
4719 helixes[i0].Evaluate(phase[0][0],xyz0);
4720 helixes[i1].Evaluate(phase[0][1],xyz1);
4722 helixes[i0].GetAngle(phase[0][0],helixes[i1],phase[0][1],hangles);
4723 Double_t deltah[2],deltabest;
4724 if (TMath::Abs(hangles[2])<kMinAngle) continue;
4728 helixes[i0].ParabolicDCA(helixes[i1],phase[0][0],phase[0][1],radius[0],deltah[0],2);
4730 helixes[i0].ParabolicDCA(helixes[i1],phase[1][0],phase[1][1],radius[1],deltah[1],2);
4731 if (deltah[1]<deltah[0]) ibest=1;
4733 deltabest = TMath::Sqrt(deltah[ibest]);
4734 helixes[i0].Evaluate(phase[ibest][0],xyz0);
4735 helixes[i1].Evaluate(phase[ibest][1],xyz1);
4736 helixes[i0].GetAngle(phase[ibest][0],helixes[i1],phase[ibest][1],hangles);
4737 Double_t radiusbest = TMath::Sqrt(radius[ibest]);
4739 if (deltabest>kMaxDist) continue;
4740 // if (mindcar+mindcaz<40 && (TMath::Abs(hangles[2])<kMinAngle ||deltabest>3)) continue;
4741 Bool_t sign =kFALSE;
4742 if (hangles[2]>kMinAngle) sign =kTRUE;
4745 // circular[i0] = kTRUE;
4746 // circular[i1] = kTRUE;
4747 if (track0->OneOverPt()<track1->OneOverPt()){
4748 track0->SetCircular(track0->GetCircular()+1);
4749 track1->SetCircular(track1->GetCircular()+2);
4752 track1->SetCircular(track1->GetCircular()+1);
4753 track0->SetCircular(track0->GetCircular()+2);
4756 if (AliTPCReconstructor::StreamLevel()>2){
4758 //debug stream to tune "fine" cuts
4759 Int_t lab0=track0->GetLabel();
4760 Int_t lab1=track1->GetLabel();
4761 TTreeSRedirector &cstream = *fDebugStreamer;
4762 cstream<<"Curling2"<<
4778 "npoints="<<npoints<<
4779 "hangles0="<<hangles[0]<<
4780 "hangles1="<<hangles[1]<<
4781 "hangles2="<<hangles[2]<<
4784 "radius="<<radiusbest<<
4785 "deltabest="<<deltabest<<
4786 "phase0="<<phase[ibest][0]<<
4787 "phase1="<<phase[ibest][1]<<
4795 if (AliTPCReconstructor::StreamLevel()>1) {
4796 AliInfo("Time for curling tracks removal");
4802 void AliTPCtrackerMI::FindKinks(TObjArray * array, AliESDEvent *esd)
4808 // RS something is wrong in this routine: not all seeds are assigned to daughters and mothers array, but they all are queried
4811 TObjArray *kinks= new TObjArray(10000);
4812 // TObjArray *v0s= new TObjArray(10000);
4813 Int_t nentries = array->GetEntriesFast();
4814 AliHelix *helixes = new AliHelix[nentries];
4815 Int_t *sign = new Int_t[nentries];
4816 Int_t *nclusters = new Int_t[nentries];
4817 Float_t *alpha = new Float_t[nentries];
4818 AliKink *kink = new AliKink();
4819 Int_t * usage = new Int_t[nentries];
4820 Float_t *zm = new Float_t[nentries];
4821 Float_t *z0 = new Float_t[nentries];
4822 Float_t *fim = new Float_t[nentries];
4823 Float_t *shared = new Float_t[nentries];
4824 Bool_t *circular = new Bool_t[nentries];
4825 Float_t *dca = new Float_t[nentries];
4826 //const AliESDVertex * primvertex = esd->GetVertex();
4828 // nentries = array->GetEntriesFast();
4833 for (Int_t i=0;i<nentries;i++){
4836 AliTPCseed* track = (AliTPCseed*)array->At(i);
4837 if (!track) continue;
4838 track->SetCircular(0);
4840 track->UpdatePoints();
4841 if (( track->GetPoints()[2]- track->GetPoints()[0])>5 && track->GetPoints()[3]>0.8){
4843 nclusters[i]=track->GetNumberOfClusters();
4844 alpha[i] = track->GetAlpha();
4845 new (&helixes[i]) AliHelix(*track);
4847 helixes[i].Evaluate(0,xyz);
4848 sign[i] = (track->GetC()>0) ? -1:1;
4851 if (track->GetProlongation(x,y,z)){
4853 fim[i] = alpha[i]+TMath::ATan2(y,x);
4856 zm[i] = track->GetZ();
4860 circular[i]= kFALSE;
4861 if (track->GetProlongation(0,y,z)) z0[i] = z;
4862 dca[i] = track->GetD(0,0);
4868 Int_t ncandidates =0;
4871 Double_t phase[2][2]={{0,0},{0,0}},radius[2]={0,0};
4874 // Find circling track
4876 for (Int_t i0=0;i0<nentries;i0++){
4877 AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
4878 if (!track0) continue;
4879 if (track0->GetNumberOfClusters()<40) continue;
4880 if (TMath::Abs(1./track0->GetC())>200) continue;
4881 for (Int_t i1=i0+1;i1<nentries;i1++){
4882 AliTPCseed * track1 = (AliTPCseed*)array->At(i1);
4883 if (!track1) continue;
4884 if (track1->GetNumberOfClusters()<40) continue;
4885 if ( TMath::Abs(track1->GetTgl()+track0->GetTgl())>0.1) continue;
4886 if (track0->GetBConstrain()&&track1->GetBConstrain()) continue;
4887 if (TMath::Abs(1./track1->GetC())>200) continue;
4888 if (track1->GetSigned1Pt()*track0->GetSigned1Pt()>0) continue;
4889 if (track1->GetTgl()*track0->GetTgl()>0) continue;
4890 if (TMath::Max(TMath::Abs(1./track0->GetC()),TMath::Abs(1./track1->GetC()))>190) continue;
4891 if (track0->GetBConstrain()&&track1->OneOverPt()<track0->OneOverPt()) continue; //returning - lower momenta
4892 if (track1->GetBConstrain()&&track0->OneOverPt()<track1->OneOverPt()) continue; //returning - lower momenta
4894 Float_t mindcar = TMath::Min(TMath::Abs(dca[i0]),TMath::Abs(dca[i1]));
4895 if (mindcar<5) continue;
4896 Float_t mindcaz = TMath::Min(TMath::Abs(z0[i0]-GetZ()),TMath::Abs(z0[i1]-GetZ()));
4897 if (mindcaz<5) continue;
4898 if (mindcar+mindcaz<20) continue;
4901 Float_t xc0 = helixes[i0].GetHelix(6);
4902 Float_t yc0 = helixes[i0].GetHelix(7);
4903 Float_t r0 = helixes[i0].GetHelix(8);
4904 Float_t xc1 = helixes[i1].GetHelix(6);
4905 Float_t yc1 = helixes[i1].GetHelix(7);
4906 Float_t r1 = helixes[i1].GetHelix(8);
4908 Float_t rmean = (r0+r1)*0.5;
4909 Float_t delta =TMath::Sqrt((xc1-xc0)*(xc1-xc0)+(yc1-yc0)*(yc1-yc0));
4910 //if (delta>30) continue;
4911 if (delta>rmean*0.25) continue;
4912 if (TMath::Abs(r0-r1)/rmean>0.3) continue;
4914 Int_t npoints = helixes[i0].GetRPHIintersections(helixes[i1], phase, radius,10);
4915 if (npoints==0) continue;
4916 helixes[i0].GetClosestPhases(helixes[i1], phase);
4920 Double_t hangles[3];
4921 helixes[i0].Evaluate(phase[0][0],xyz0);
4922 helixes[i1].Evaluate(phase[0][1],xyz1);
4924 helixes[i0].GetAngle(phase[0][0],helixes[i1],phase[0][1],hangles);
4925 Double_t deltah[2],deltabest;
4926 if (hangles[2]<2.8) continue;
4929 helixes[i0].ParabolicDCA(helixes[i1],phase[0][0],phase[0][1],radius[0],deltah[0],2);
4931 helixes[i0].ParabolicDCA(helixes[i1],phase[1][0],phase[1][1],radius[1],deltah[1],2);
4932 if (deltah[1]<deltah[0]) ibest=1;
4934 deltabest = TMath::Sqrt(deltah[ibest]);
4935 helixes[i0].Evaluate(phase[ibest][0],xyz0);
4936 helixes[i1].Evaluate(phase[ibest][1],xyz1);
4937 helixes[i0].GetAngle(phase[ibest][0],helixes[i1],phase[ibest][1],hangles);
4938 Double_t radiusbest = TMath::Sqrt(radius[ibest]);
4940 if (deltabest>6) continue;
4941 if (mindcar+mindcaz<40 && (hangles[2]<3.12||deltabest>3)) continue;
4942 Bool_t lsign =kFALSE;
4943 if (hangles[2]>3.06) lsign =kTRUE;
4946 circular[i0] = kTRUE;
4947 circular[i1] = kTRUE;
4948 if (track0->OneOverPt()<track1->OneOverPt()){
4949 track0->SetCircular(track0->GetCircular()+1);
4950 track1->SetCircular(track1->GetCircular()+2);
4953 track1->SetCircular(track1->GetCircular()+1);
4954 track0->SetCircular(track0->GetCircular()+2);
4957 if (lsign&&AliTPCReconstructor::StreamLevel()>1){
4959 Int_t lab0=track0->GetLabel();
4960 Int_t lab1=track1->GetLabel();
4961 TTreeSRedirector &cstream = *fDebugStreamer;
4962 cstream<<"Curling"<<
4969 "mindcar="<<mindcar<<
4970 "mindcaz="<<mindcaz<<
4973 "npoints="<<npoints<<
4974 "hangles0="<<hangles[0]<<
4975 "hangles2="<<hangles[2]<<
4980 "radius="<<radiusbest<<
4981 "deltabest="<<deltabest<<
4982 "phase0="<<phase[ibest][0]<<
4983 "phase1="<<phase[ibest][1]<<
4993 for (Int_t i =0;i<nentries;i++){
4994 if (sign[i]==0) continue;
4995 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5002 Double_t cradius0 = 40*40;
5003 Double_t cradius1 = 270*270;
5006 Double_t cdist3=0.55;
5007 for (Int_t j =i+1;j<nentries;j++){
5009 if (sign[j]*sign[i]<1) continue;
5010 if ( (nclusters[i]+nclusters[j])>200) continue;
5011 if ( (nclusters[i]+nclusters[j])<80) continue;
5012 if ( TMath::Abs(zm[i]-zm[j])>60.) continue;
5013 if ( TMath::Abs(fim[i]-fim[j])>0.6 && TMath::Abs(fim[i]-fim[j])<5.7 ) continue;
5014 //AliTPCseed * track1 = (AliTPCseed*)array->At(j); Double_t phase[2][2],radius[2];
5015 Int_t npoints = helixes[i].GetRPHIintersections(helixes[j], phase, radius,20);
5016 if (npoints<1) continue;
5019 if (radius[0]<cradius0||radius[0]>cradius1) continue;
5022 if ( (radius[0]<cradius0||radius[0]>cradius1) && (radius[1]<cradius0||radius[1]>cradius1) ) continue;
5025 Double_t delta1=10000,delta2=10000;
5026 // cuts on the intersection radius
5027 helixes[i].LinearDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta1);
5028 if (radius[0]<20&&delta1<1) continue; //intersection at vertex
5029 if (radius[0]<10&&delta1<3) continue; //intersection at vertex
5031 helixes[i].LinearDCA(helixes[j],phase[1][0],phase[1][1],radius[1],delta2);
5032 if (radius[1]<20&&delta2<1) continue; //intersection at vertex
5033 if (radius[1]<10&&delta2<3) continue; //intersection at vertex
5036 Double_t distance1 = TMath::Min(delta1,delta2);
5037 if (distance1>cdist1) continue; // cut on DCA linear approximation
5039 npoints = helixes[i].GetRPHIintersections(helixes[j], phase, radius,20);
5040 helixes[i].ParabolicDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta1);
5041 if (radius[0]<20&&delta1<1) continue; //intersection at vertex
5042 if (radius[0]<10&&delta1<3) continue; //intersection at vertex
5045 helixes[i].ParabolicDCA(helixes[j],phase[1][0],phase[1][1],radius[1],delta2);
5046 if (radius[1]<20&&delta2<1) continue; //intersection at vertex
5047 if (radius[1]<10&&delta2<3) continue; //intersection at vertex
5049 distance1 = TMath::Min(delta1,delta2);
5052 rkink = TMath::Sqrt(radius[0]);
5055 rkink = TMath::Sqrt(radius[1]);
5057 if (distance1>cdist2) continue;
5060 AliTPCseed * track1 = (AliTPCseed*)array->At(j);
5063 Int_t row0 = GetRowNumber(rkink);
5064 if (row0<10) continue;
5065 if (row0>150) continue;
5068 Float_t dens00=-1,dens01=-1;
5069 Float_t dens10=-1,dens11=-1;
5071 Int_t found,foundable,ishared;
5072 track0->GetClusterStatistic(0,row0-5, found, foundable,ishared,kFALSE);
5073 if (foundable>5) dens00 = Float_t(found)/Float_t(foundable);
5074 track0->GetClusterStatistic(row0+5,155, found, foundable,ishared,kFALSE);
5075 if (foundable>5) dens01 = Float_t(found)/Float_t(foundable);
5077 track1->GetClusterStatistic(0,row0-5, found, foundable,ishared,kFALSE);
5078 if (foundable>10) dens10 = Float_t(found)/Float_t(foundable);
5079 track1->GetClusterStatistic(row0+5,155, found, foundable,ishared,kFALSE);
5080 if (foundable>10) dens11 = Float_t(found)/Float_t(foundable);
5082 if (dens00<dens10 && dens01<dens11) continue;
5083 if (dens00>dens10 && dens01>dens11) continue;
5084 if (TMath::Max(dens00,dens10)<0.1) continue;
5085 if (TMath::Max(dens01,dens11)<0.3) continue;
5087 if (TMath::Min(dens00,dens10)>0.6) continue;
5088 if (TMath::Min(dens01,dens11)>0.6) continue;
5091 AliTPCseed * ktrack0, *ktrack1;
5100 if (TMath::Abs(ktrack0->GetC())>5) continue; // cut on the curvature for mother particle
5101 AliExternalTrackParam paramm(*ktrack0);
5102 AliExternalTrackParam paramd(*ktrack1);
5103 if (row0>60&&ktrack1->GetReference().GetX()>90.)new (¶md) AliExternalTrackParam(ktrack1->GetReference());
5106 kink->SetMother(paramm);
5107 kink->SetDaughter(paramd);
5110 Float_t x[3] = { kink->GetPosition()[0],kink->GetPosition()[1],kink->GetPosition()[2]};
5112 fkParam->Transform0to1(x,index);
5113 fkParam->Transform1to2(x,index);
5114 row0 = GetRowNumber(x[0]);
5116 if (kink->GetR()<100) continue;
5117 if (kink->GetR()>240) continue;
5118 if (kink->GetPosition()[2]/kink->GetR()>AliTPCReconstructor::GetCtgRange()) continue; //out of fiducial volume
5119 if (kink->GetDistance()>cdist3) continue;
5120 Float_t dird = kink->GetDaughterP()[0]*kink->GetPosition()[0]+kink->GetDaughterP()[1]*kink->GetPosition()[1]; // rough direction estimate
5121 if (dird<0) continue;
5123 Float_t dirm = kink->GetMotherP()[0]*kink->GetPosition()[0]+kink->GetMotherP()[1]*kink->GetPosition()[1]; // rough direction estimate
5124 if (dirm<0) continue;
5125 Float_t mpt = TMath::Sqrt(kink->GetMotherP()[0]*kink->GetMotherP()[0]+kink->GetMotherP()[1]*kink->GetMotherP()[1]);
5126 if (mpt<0.2) continue;
5129 //for high momenta momentum not defined well in first iteration
5130 Double_t qt = TMath::Sin(kink->GetAngle(2))*ktrack1->GetP();
5131 if (qt>0.35) continue;
5134 kink->SetLabel(CookLabel(ktrack0,0.4,0,row0),0);
5135 kink->SetLabel(CookLabel(ktrack1,0.4,row0,160),1);
5137 kink->SetTPCDensity(dens00,0,0);
5138 kink->SetTPCDensity(dens01,0,1);
5139 kink->SetTPCDensity(dens10,1,0);
5140 kink->SetTPCDensity(dens11,1,1);
5141 kink->SetIndex(i,0);
5142 kink->SetIndex(j,1);
5145 kink->SetTPCDensity(dens10,0,0);
5146 kink->SetTPCDensity(dens11,0,1);
5147 kink->SetTPCDensity(dens00,1,0);
5148 kink->SetTPCDensity(dens01,1,1);
5149 kink->SetIndex(j,0);
5150 kink->SetIndex(i,1);
5153 if (mpt<1||kink->GetAngle(2)>0.1){
5154 // angle and densities not defined yet
5155 if (kink->GetTPCDensityFactor()<0.8) continue;
5156 if ((2-kink->GetTPCDensityFactor())*kink->GetDistance() >0.25) continue;
5157 if (kink->GetAngle(2)*ktrack0->GetP()<0.003) continue; //too small angle
5158 if (kink->GetAngle(2)>0.2&&kink->GetTPCDensityFactor()<1.15) continue;
5159 if (kink->GetAngle(2)>0.2&&kink->GetTPCDensity(0,1)>0.05) continue;
5161 Float_t criticalangle = track0->GetSigmaSnp2()+track0->GetSigmaTgl2();
5162 criticalangle+= track1->GetSigmaSnp2()+track1->GetSigmaTgl2();
5163 criticalangle= 3*TMath::Sqrt(criticalangle);
5164 if (criticalangle>0.02) criticalangle=0.02;
5165 if (kink->GetAngle(2)<criticalangle) continue;
5168 Int_t drow = Int_t(2.+0.5/(0.05+kink->GetAngle(2))); // overlap region defined
5169 Float_t shapesum =0;
5171 for ( Int_t row = row0-drow; row<row0+drow;row++){
5172 if (row<0) continue;
5173 if (row>155) continue;
5174 if (ktrack0->GetClusterPointer(row)){
5175 AliTPCTrackerPoint *point =ktrack0->GetTrackPoint(row);
5176 shapesum+=point->GetSigmaY()+point->GetSigmaZ();
5179 if (ktrack1->GetClusterPointer(row)){
5180 AliTPCTrackerPoint *point =ktrack1->GetTrackPoint(row);
5181 shapesum+=point->GetSigmaY()+point->GetSigmaZ();
5186 kink->SetShapeFactor(-1.);
5189 kink->SetShapeFactor(shapesum/sum);
5191 // esd->AddKink(kink);
5193 // kink->SetMother(paramm);
5194 //kink->SetDaughter(paramd);
5196 Double_t chi2P2 = paramm.GetParameter()[2]-paramd.GetParameter()[2];
5198 chi2P2/=paramm.GetCovariance()[5]+paramd.GetCovariance()[5];
5199 Double_t chi2P3 = paramm.GetParameter()[3]-paramd.GetParameter()[3];
5201 chi2P3/=paramm.GetCovariance()[9]+paramd.GetCovariance()[9];
5203 if (AliTPCReconstructor::StreamLevel()>1) {
5204 (*fDebugStreamer)<<"kinkLpt"<<
5212 if ( chi2P2+chi2P3<AliTPCReconstructor::GetRecoParam()->GetKinkAngleCutChi2(0)){
5216 kinks->AddLast(kink);
5222 // sort the kinks according quality - and refit them towards vertex
5224 Int_t nkinks = kinks->GetEntriesFast();
5225 Float_t *quality = new Float_t[nkinks];
5226 Int_t *indexes = new Int_t[nkinks];
5227 AliTPCseed *mothers = new AliTPCseed[nkinks];
5228 AliTPCseed *daughters = new AliTPCseed[nkinks];
5231 for (Int_t i=0;i<nkinks;i++){
5233 AliKink *kinkl = (AliKink*)kinks->At(i);
5235 // refit kinks towards vertex
5237 Int_t index0 = kinkl->GetIndex(0);
5238 Int_t index1 = kinkl->GetIndex(1);
5239 AliTPCseed * ktrack0 = (AliTPCseed*)array->At(index0);
5240 AliTPCseed * ktrack1 = (AliTPCseed*)array->At(index1);
5242 Int_t sumn=ktrack0->GetNumberOfClusters()+ktrack1->GetNumberOfClusters();
5244 // Refit Kink under if too small angle
5246 if (kinkl->GetAngle(2)<0.05){
5247 kinkl->SetTPCRow0(GetRowNumber(kinkl->GetR()));
5248 Int_t row0 = kinkl->GetTPCRow0();
5249 Int_t drow = Int_t(2.+0.5/(0.05+kinkl->GetAngle(2)));
5252 Int_t last = row0-drow;
5253 if (last<40) last=40;
5254 if (last<ktrack0->GetFirstPoint()+25) last = ktrack0->GetFirstPoint()+25;
5255 AliTPCseed* seed0 = ReSeed(ktrack0,last,kFALSE);
5258 Int_t first = row0+drow;
5259 if (first>130) first=130;
5260 if (first>ktrack1->GetLastPoint()-25) first = TMath::Max(ktrack1->GetLastPoint()-25,30);
5261 AliTPCseed* seed1 = ReSeed(ktrack1,first,kTRUE);
5263 if (seed0 && seed1){
5264 kinkl->SetStatus(1,8);
5265 if (RefitKink(*seed0,*seed1,*kinkl)) kinkl->SetStatus(1,9);
5266 row0 = GetRowNumber(kinkl->GetR());
5267 sumn = seed0->GetNumberOfClusters()+seed1->GetNumberOfClusters();
5268 mothers[i] = *seed0;
5269 daughters[i] = *seed1;
5272 delete kinks->RemoveAt(i);
5273 if (seed0) MarkSeedFree( seed0 );
5274 if (seed1) MarkSeedFree( seed1 );
5277 if (kinkl->GetDistance()>0.5 || kinkl->GetR()<110 || kinkl->GetR()>240) {
5278 delete kinks->RemoveAt(i);
5279 if (seed0) MarkSeedFree( seed0 );
5280 if (seed1) MarkSeedFree( seed1 );
5284 MarkSeedFree( seed0 );
5285 MarkSeedFree( seed1 );
5288 if (kinkl) quality[i] = 160*((0.1+kinkl->GetDistance())*(2.-kinkl->GetTPCDensityFactor()))/(sumn+40.); //the longest -clossest will win
5290 TMath::Sort(nkinks,quality,indexes,kFALSE);
5292 //remove double find kinks
5294 for (Int_t ikink0=1;ikink0<nkinks;ikink0++){
5295 AliKink * kink0 = (AliKink*) kinks->At(indexes[ikink0]);
5296 if (!kink0) continue;
5298 for (Int_t ikink1=0;ikink1<ikink0;ikink1++){
5299 kink0 = (AliKink*) kinks->At(indexes[ikink0]);
5300 if (!kink0) continue;
5301 AliKink * kink1 = (AliKink*) kinks->At(indexes[ikink1]);
5302 if (!kink1) continue;
5303 // if not close kink continue
5304 if (TMath::Abs(kink1->GetPosition()[2]-kink0->GetPosition()[2])>10) continue;
5305 if (TMath::Abs(kink1->GetPosition()[1]-kink0->GetPosition()[1])>10) continue;
5306 if (TMath::Abs(kink1->GetPosition()[0]-kink0->GetPosition()[0])>10) continue;
5308 AliTPCseed &mother0 = mothers[indexes[ikink0]];
5309 AliTPCseed &daughter0 = daughters[indexes[ikink0]];
5310 AliTPCseed &mother1 = mothers[indexes[ikink1]];
5311 AliTPCseed &daughter1 = daughters[indexes[ikink1]];
5312 Int_t row0 = (kink0->GetTPCRow0()+kink1->GetTPCRow0())/2;
5321 for (Int_t i=0;i<row0;i++){
5322 if (mother0.GetClusterIndex(i)>0 && mother1.GetClusterIndex(i)>0){
5325 if (mother0.GetClusterIndex(i)==mother1.GetClusterIndex(i)){
5332 for (Int_t i=row0;i<158;i++){
5333 //if (daughter0.GetClusterIndex(i)>0 && daughter0.GetClusterIndex(i)>0){ // RS: Bug?
5334 if (daughter0.GetClusterIndex(i)>0 && daughter1.GetClusterIndex(i)>0){
5337 if (mother0.GetClusterIndex(i)==mother1.GetClusterIndex(i)){
5343 Float_t ratio = Float_t(same+1)/Float_t(both+1);
5344 Float_t ratiom = Float_t(samem+1)/Float_t(bothm+1);
5345 Float_t ratiod = Float_t(samed+1)/Float_t(bothd+1);
5346 if (ratio>0.3 && ratiom>0.5 &&ratiod>0.5) {
5347 Int_t sum0 = mother0.GetNumberOfClusters()+daughter0.GetNumberOfClusters();
5348 Int_t sum1 = mother1.GetNumberOfClusters()+daughter1.GetNumberOfClusters();
5350 shared[kink0->GetIndex(0)]= kTRUE;
5351 shared[kink0->GetIndex(1)]= kTRUE;
5352 delete kinks->RemoveAt(indexes[ikink0]);
5356 shared[kink1->GetIndex(0)]= kTRUE;
5357 shared[kink1->GetIndex(1)]= kTRUE;
5358 delete kinks->RemoveAt(indexes[ikink1]);
5365 for (Int_t i=0;i<nkinks;i++){
5366 AliKink * kinkl = (AliKink*) kinks->At(indexes[i]);
5367 if (!kinkl) continue;
5368 kinkl->SetTPCRow0(GetRowNumber(kinkl->GetR()));
5369 Int_t index0 = kinkl->GetIndex(0);
5370 Int_t index1 = kinkl->GetIndex(1);
5371 if (circular[index0]||(circular[index1]&&kinkl->GetDistance()>0.2)) continue;
5372 kinkl->SetMultiple(usage[index0],0);
5373 kinkl->SetMultiple(usage[index1],1);
5374 if (kinkl->GetMultiple()[0]+kinkl->GetMultiple()[1]>2) continue;
5375 if (kinkl->GetMultiple()[0]+kinkl->GetMultiple()[1]>0 && quality[indexes[i]]>0.2) continue;
5376 if (kinkl->GetMultiple()[0]+kinkl->GetMultiple()[1]>0 && kinkl->GetDistance()>0.2) continue;
5377 if (circular[index0]||(circular[index1]&&kinkl->GetDistance()>0.1)) continue;
5379 AliTPCseed * ktrack0 = (AliTPCseed*)array->At(index0);
5380 AliTPCseed * ktrack1 = (AliTPCseed*)array->At(index1);
5381 if (!ktrack0 || !ktrack1) continue;
5382 Int_t index = esd->AddKink(kinkl);
5385 if ( ktrack0->GetKinkIndex(0)==0 && ktrack1->GetKinkIndex(0)==0) { //best kink
5386 if (mothers[indexes[i]].GetNumberOfClusters()>20 && daughters[indexes[i]].GetNumberOfClusters()>20 && (mothers[indexes[i]].GetNumberOfClusters()+daughters[indexes[i]].GetNumberOfClusters())>100){
5387 *ktrack0 = mothers[indexes[i]];
5388 *ktrack1 = daughters[indexes[i]];
5392 ktrack0->SetKinkIndex(usage[index0],-(index+1));
5393 ktrack1->SetKinkIndex(usage[index1], (index+1));
5398 // Remove tracks corresponding to shared kink's
5400 for (Int_t i=0;i<nentries;i++){
5401 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5402 if (!track0) continue;
5403 if (track0->GetKinkIndex(0)!=0) continue;
5404 if (shared[i]) MarkSeedFree( array->RemoveAt(i) );
5409 RemoveUsed2(array,0.5,0.4,30);
5411 for (Int_t i=0;i<nentries;i++){
5412 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5413 if (!track0) continue;
5414 track0->CookdEdx(0.02,0.6);
5418 for (Int_t i=0;i<nentries;i++){
5419 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5420 if (!track0) continue;
5421 if (track0->Pt()<1.4) continue;
5422 //remove double high momenta tracks - overlapped with kink candidates
5425 for (Int_t icl=track0->GetFirstPoint();icl<track0->GetLastPoint(); icl++){
5426 if (track0->GetClusterPointer(icl)!=0){
5428 if (track0->GetClusterPointer(icl)->IsUsed(10)) ishared++;
5431 if (Float_t(ishared+1)/Float_t(all+1)>0.5) {
5432 MarkSeedFree( array->RemoveAt(i) );
5436 if (track0->GetKinkIndex(0)!=0) continue;
5437 if (track0->GetNumberOfClusters()<80) continue;
5439 AliTPCseed *pmother = new AliTPCseed();
5440 AliTPCseed *pdaughter = new AliTPCseed();
5441 AliKink *pkink = new AliKink;
5443 AliTPCseed & mother = *pmother;
5444 AliTPCseed & daughter = *pdaughter;
5445 AliKink & kinkl = *pkink;
5446 if (CheckKinkPoint(track0,mother,daughter, kinkl)){
5447 if (mother.GetNumberOfClusters()<30||daughter.GetNumberOfClusters()<20) {
5451 continue; //too short tracks
5453 if (mother.Pt()<1.4) {
5459 Int_t row0= kinkl.GetTPCRow0();
5460 if (kinkl.GetDistance()>0.5 || kinkl.GetR()<110. || kinkl.GetR()>240.) {
5467 Int_t index = esd->AddKink(&kinkl);
5468 mother.SetKinkIndex(0,-(index+1));
5469 daughter.SetKinkIndex(0,index+1);
5470 if (mother.GetNumberOfClusters()>50) {
5471 MarkSeedFree( array->RemoveAt(i) );
5472 AliTPCseed* mtc = new( NextFreeSeed() ) AliTPCseed(mother);
5473 mtc->SetPoolID(fLastSeedID);
5474 array->AddAt(mtc,i);
5477 AliTPCseed* mtc = new( NextFreeSeed() ) AliTPCseed(mother);
5478 mtc->SetPoolID(fLastSeedID);
5479 array->AddLast(mtc);
5481 AliTPCseed* dtc = new( NextFreeSeed() ) AliTPCseed(daughter);
5482 dtc->SetPoolID(fLastSeedID);
5483 array->AddLast(dtc);
5484 for (Int_t icl=0;icl<row0;icl++) {
5485 if (mother.GetClusterPointer(icl)) mother.GetClusterPointer(icl)->Use(20);
5488 for (Int_t icl=row0;icl<158;icl++) {
5489 if (daughter.GetClusterPointer(icl)) daughter.GetClusterPointer(icl)->Use(20);
5498 delete [] daughters;
5520 AliInfo(Form("Ncandidates=\t%d\t%d\t%d\t%d\n",esd->GetNumberOfKinks(),ncandidates,ntracks,nall));
5526 void AliTPCtrackerMI::FindKinks(TObjArray * array, AliESDEvent *esd)
5533 TObjArray *kinks= new TObjArray(10000);
5534 // TObjArray *v0s= new TObjArray(10000);
5535 Int_t nentries = array->GetEntriesFast();
5536 AliHelix *helixes = new AliHelix[nentries];
5537 Int_t *sign = new Int_t[nentries];
5538 Int_t *nclusters = new Int_t[nentries];
5539 Float_t *alpha = new Float_t[nentries];
5540 AliKink *kink = new AliKink();
5541 Int_t * usage = new Int_t[nentries];
5542 Float_t *zm = new Float_t[nentries];
5543 Float_t *z0 = new Float_t[nentries];
5544 Float_t *fim = new Float_t[nentries];
5545 Float_t *shared = new Float_t[nentries];
5546 Bool_t *circular = new Bool_t[nentries];
5547 Float_t *dca = new Float_t[nentries];
5548 //const AliESDVertex * primvertex = esd->GetVertex();
5550 // nentries = array->GetEntriesFast();
5555 for (Int_t i=0;i<nentries;i++){
5558 AliTPCseed* track = (AliTPCseed*)array->At(i);
5559 if (!track) continue;
5560 track->SetCircular(0);
5562 track->UpdatePoints();
5563 if (( track->GetPoints()[2]- track->GetPoints()[0])>5 && track->GetPoints()[3]>0.8){
5565 nclusters[i]=track->GetNumberOfClusters();
5566 alpha[i] = track->GetAlpha();
5567 new (&helixes[i]) AliHelix(*track);
5569 helixes[i].Evaluate(0,xyz);
5570 sign[i] = (track->GetC()>0) ? -1:1;
5573 if (track->GetProlongation(x,y,z)){
5575 fim[i] = alpha[i]+TMath::ATan2(y,x);
5578 zm[i] = track->GetZ();
5582 circular[i]= kFALSE;
5583 if (track->GetProlongation(0,y,z)) z0[i] = z;
5584 dca[i] = track->GetD(0,0);
5590 Int_t ncandidates =0;
5593 Double_t phase[2][2]={{0,0},{0,0}},radius[2]={0,0};
5596 // Find circling track
5598 for (Int_t i0=0;i0<nentries;i0++){
5599 AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
5600 if (!track0) continue;
5601 if (track0->GetNumberOfClusters()<40) continue;
5602 if (TMath::Abs(1./track0->GetC())>200) continue;
5603 for (Int_t i1=i0+1;i1<nentries;i1++){
5604 AliTPCseed * track1 = (AliTPCseed*)array->At(i1);
5605 if (!track1) continue;
5606 if (track1->GetNumberOfClusters()<40) continue;
5607 if ( TMath::Abs(track1->GetTgl()+track0->GetTgl())>0.1) continue;
5608 if (track0->GetBConstrain()&&track1->GetBConstrain()) continue;
5609 if (TMath::Abs(1./track1->GetC())>200) continue;
5610 if (track1->GetSigned1Pt()*track0->GetSigned1Pt()>0) continue;
5611 if (track1->GetTgl()*track0->GetTgl()>0) continue;
5612 if (TMath::Max(TMath::Abs(1./track0->GetC()),TMath::Abs(1./track1->GetC()))>190) continue;
5613 if (track0->GetBConstrain()&&track1->OneOverPt()<track0->OneOverPt()) continue; //returning - lower momenta
5614 if (track1->GetBConstrain()&&track0->OneOverPt()<track1->OneOverPt()) continue; //returning - lower momenta
5616 Float_t mindcar = TMath::Min(TMath::Abs(dca[i0]),TMath::Abs(dca[i1]));
5617 if (mindcar<5) continue;
5618 Float_t mindcaz = TMath::Min(TMath::Abs(z0[i0]-GetZ()),TMath::Abs(z0[i1]-GetZ()));
5619 if (mindcaz<5) continue;
5620 if (mindcar+mindcaz<20) continue;
5623 Float_t xc0 = helixes[i0].GetHelix(6);
5624 Float_t yc0 = helixes[i0].GetHelix(7);
5625 Float_t r0 = helixes[i0].GetHelix(8);
5626 Float_t xc1 = helixes[i1].GetHelix(6);
5627 Float_t yc1 = helixes[i1].GetHelix(7);
5628 Float_t r1 = helixes[i1].GetHelix(8);
5630 Float_t rmean = (r0+r1)*0.5;
5631 Float_t delta =TMath::Sqrt((xc1-xc0)*(xc1-xc0)+(yc1-yc0)*(yc1-yc0));
5632 //if (delta>30) continue;
5633 if (delta>rmean*0.25) continue;
5634 if (TMath::Abs(r0-r1)/rmean>0.3) continue;
5636 Int_t npoints = helixes[i0].GetRPHIintersections(helixes[i1], phase, radius,10);
5637 if (npoints==0) continue;
5638 helixes[i0].GetClosestPhases(helixes[i1], phase);
5642 Double_t hangles[3];
5643 helixes[i0].Evaluate(phase[0][0],xyz0);
5644 helixes[i1].Evaluate(phase[0][1],xyz1);
5646 helixes[i0].GetAngle(phase[0][0],helixes[i1],phase[0][1],hangles);
5647 Double_t deltah[2],deltabest;
5648 if (hangles[2]<2.8) continue;
5651 helixes[i0].ParabolicDCA(helixes[i1],phase[0][0],phase[0][1],radius[0],deltah[0],2);
5653 helixes[i0].ParabolicDCA(helixes[i1],phase[1][0],phase[1][1],radius[1],deltah[1],2);
5654 if (deltah[1]<deltah[0]) ibest=1;
5656 deltabest = TMath::Sqrt(deltah[ibest]);
5657 helixes[i0].Evaluate(phase[ibest][0],xyz0);
5658 helixes[i1].Evaluate(phase[ibest][1],xyz1);
5659 helixes[i0].GetAngle(phase[ibest][0],helixes[i1],phase[ibest][1],hangles);
5660 Double_t radiusbest = TMath::Sqrt(radius[ibest]);
5662 if (deltabest>6) continue;
5663 if (mindcar+mindcaz<40 && (hangles[2]<3.12||deltabest>3)) continue;
5664 Bool_t lsign =kFALSE;
5665 if (hangles[2]>3.06) lsign =kTRUE;
5668 circular[i0] = kTRUE;
5669 circular[i1] = kTRUE;
5670 if (track0->OneOverPt()<track1->OneOverPt()){
5671 track0->SetCircular(track0->GetCircular()+1);
5672 track1->SetCircular(track1->GetCircular()+2);
5675 track1->SetCircular(track1->GetCircular()+1);
5676 track0->SetCircular(track0->GetCircular()+2);
5679 if (lsign&&AliTPCReconstructor::StreamLevel()>1){
5681 Int_t lab0=track0->GetLabel();
5682 Int_t lab1=track1->GetLabel();
5683 TTreeSRedirector &cstream = *fDebugStreamer;
5684 cstream<<"Curling"<<
5691 "mindcar="<<mindcar<<
5692 "mindcaz="<<mindcaz<<
5695 "npoints="<<npoints<<
5696 "hangles0="<<hangles[0]<<
5697 "hangles2="<<hangles[2]<<
5702 "radius="<<radiusbest<<
5703 "deltabest="<<deltabest<<
5704 "phase0="<<phase[ibest][0]<<
5705 "phase1="<<phase[ibest][1]<<
5715 for (Int_t i =0;i<nentries;i++){
5716 if (sign[i]==0) continue;
5717 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5724 Double_t cradius0 = 40*40;
5725 Double_t cradius1 = 270*270;
5728 Double_t cdist3=0.55;
5729 for (Int_t j =i+1;j<nentries;j++){
5731 if (sign[j]*sign[i]<1) continue;
5732 if ( (nclusters[i]+nclusters[j])>200) continue;
5733 if ( (nclusters[i]+nclusters[j])<80) continue;
5734 if ( TMath::Abs(zm[i]-zm[j])>60.) continue;
5735 if ( TMath::Abs(fim[i]-fim[j])>0.6 && TMath::Abs(fim[i]-fim[j])<5.7 ) continue;
5736 //AliTPCseed * track1 = (AliTPCseed*)array->At(j); Double_t phase[2][2],radius[2];
5737 Int_t npoints = helixes[i].GetRPHIintersections(helixes[j], phase, radius,20);
5738 if (npoints<1) continue;
5741 if (radius[0]<cradius0||radius[0]>cradius1) continue;
5744 if ( (radius[0]<cradius0||radius[0]>cradius1) && (radius[1]<cradius0||radius[1]>cradius1) ) continue;
5747 Double_t delta1=10000,delta2=10000;
5748 // cuts on the intersection radius
5749 helixes[i].LinearDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta1);
5750 if (radius[0]<20&&delta1<1) continue; //intersection at vertex
5751 if (radius[0]<10&&delta1<3) continue; //intersection at vertex
5753 helixes[i].LinearDCA(helixes[j],phase[1][0],phase[1][1],radius[1],delta2);
5754 if (radius[1]<20&&delta2<1) continue; //intersection at vertex
5755 if (radius[1]<10&&delta2<3) continue; //intersection at vertex
5758 Double_t distance1 = TMath::Min(delta1,delta2);
5759 if (distance1>cdist1) continue; // cut on DCA linear approximation
5761 npoints = helixes[i].GetRPHIintersections(helixes[j], phase, radius,20);
5762 helixes[i].ParabolicDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta1);
5763 if (radius[0]<20&&delta1<1) continue; //intersection at vertex
5764 if (radius[0]<10&&delta1<3) continue; //intersection at vertex
5767 helixes[i].ParabolicDCA(helixes[j],phase[1][0],phase[1][1],radius[1],delta2);
5768 if (radius[1]<20&&delta2<1) continue; //intersection at vertex
5769 if (radius[1]<10&&delta2<3) continue; //intersection at vertex
5771 distance1 = TMath::Min(delta1,delta2);
5774 rkink = TMath::Sqrt(radius[0]);
5777 rkink = TMath::Sqrt(radius[1]);
5779 if (distance1>cdist2) continue;
5782 AliTPCseed * track1 = (AliTPCseed*)array->At(j);
5785 Int_t row0 = GetRowNumber(rkink);
5786 if (row0<10) continue;
5787 if (row0>150) continue;
5790 Float_t dens00=-1,dens01=-1;
5791 Float_t dens10=-1,dens11=-1;
5793 Int_t found,foundable,ishared;
5794 track0->GetClusterStatistic(0,row0-5, found, foundable,ishared,kFALSE);
5795 if (foundable>5) dens00 = Float_t(found)/Float_t(foundable);
5796 track0->GetClusterStatistic(row0+5,155, found, foundable,ishared,kFALSE);
5797 if (foundable>5) dens01 = Float_t(found)/Float_t(foundable);
5799 track1->GetClusterStatistic(0,row0-5, found, foundable,ishared,kFALSE);
5800 if (foundable>10) dens10 = Float_t(found)/Float_t(foundable);
5801 track1->GetClusterStatistic(row0+5,155, found, foundable,ishared,kFALSE);
5802 if (foundable>10) dens11 = Float_t(found)/Float_t(foundable);
5804 if (dens00<dens10 && dens01<dens11) continue;
5805 if (dens00>dens10 && dens01>dens11) continue;
5806 if (TMath::Max(dens00,dens10)<0.1) continue;
5807 if (TMath::Max(dens01,dens11)<0.3) continue;
5809 if (TMath::Min(dens00,dens10)>0.6) continue;
5810 if (TMath::Min(dens01,dens11)>0.6) continue;
5813 AliTPCseed * ktrack0, *ktrack1;
5822 if (TMath::Abs(ktrack0->GetC())>5) continue; // cut on the curvature for mother particle
5823 AliExternalTrackParam paramm(*ktrack0);
5824 AliExternalTrackParam paramd(*ktrack1);
5825 if (row0>60&&ktrack1->GetReference().GetX()>90.)new (¶md) AliExternalTrackParam(ktrack1->GetReference());
5828 kink->SetMother(paramm);
5829 kink->SetDaughter(paramd);
5832 Float_t x[3] = { kink->GetPosition()[0],kink->GetPosition()[1],kink->GetPosition()[2]};
5834 fkParam->Transform0to1(x,index);
5835 fkParam->Transform1to2(x,index);
5836 row0 = GetRowNumber(x[0]);
5838 if (kink->GetR()<100) continue;
5839 if (kink->GetR()>240) continue;
5840 if (kink->GetPosition()[2]/kink->GetR()>AliTPCReconstructor::GetCtgRange()) continue; //out of fiducial volume
5841 if (kink->GetDistance()>cdist3) continue;
5842 Float_t dird = kink->GetDaughterP()[0]*kink->GetPosition()[0]+kink->GetDaughterP()[1]*kink->GetPosition()[1]; // rough direction estimate
5843 if (dird<0) continue;
5845 Float_t dirm = kink->GetMotherP()[0]*kink->GetPosition()[0]+kink->GetMotherP()[1]*kink->GetPosition()[1]; // rough direction estimate
5846 if (dirm<0) continue;
5847 Float_t mpt = TMath::Sqrt(kink->GetMotherP()[0]*kink->GetMotherP()[0]+kink->GetMotherP()[1]*kink->GetMotherP()[1]);
5848 if (mpt<0.2) continue;
5851 //for high momenta momentum not defined well in first iteration
5852 Double_t qt = TMath::Sin(kink->GetAngle(2))*ktrack1->GetP();
5853 if (qt>0.35) continue;
5856 kink->SetLabel(CookLabel(ktrack0,0.4,0,row0),0);
5857 kink->SetLabel(CookLabel(ktrack1,0.4,row0,160),1);
5859 kink->SetTPCDensity(dens00,0,0);
5860 kink->SetTPCDensity(dens01,0,1);
5861 kink->SetTPCDensity(dens10,1,0);
5862 kink->SetTPCDensity(dens11,1,1);
5863 kink->SetIndex(i,0);
5864 kink->SetIndex(j,1);
5867 kink->SetTPCDensity(dens10,0,0);
5868 kink->SetTPCDensity(dens11,0,1);
5869 kink->SetTPCDensity(dens00,1,0);
5870 kink->SetTPCDensity(dens01,1,1);
5871 kink->SetIndex(j,0);
5872 kink->SetIndex(i,1);
5875 if (mpt<1||kink->GetAngle(2)>0.1){
5876 // angle and densities not defined yet
5877 if (kink->GetTPCDensityFactor()<0.8) continue;
5878 if ((2-kink->GetTPCDensityFactor())*kink->GetDistance() >0.25) continue;
5879 if (kink->GetAngle(2)*ktrack0->GetP()<0.003) continue; //too small angle
5880 if (kink->GetAngle(2)>0.2&&kink->GetTPCDensityFactor()<1.15) continue;
5881 if (kink->GetAngle(2)>0.2&&kink->GetTPCDensity(0,1)>0.05) continue;
5883 Float_t criticalangle = track0->GetSigmaSnp2()+track0->GetSigmaTgl2();
5884 criticalangle+= track1->GetSigmaSnp2()+track1->GetSigmaTgl2();
5885 criticalangle= 3*TMath::Sqrt(criticalangle);
5886 if (criticalangle>0.02) criticalangle=0.02;
5887 if (kink->GetAngle(2)<criticalangle) continue;
5890 Int_t drow = Int_t(2.+0.5/(0.05+kink->GetAngle(2))); // overlap region defined
5891 Float_t shapesum =0;
5893 for ( Int_t row = row0-drow; row<row0+drow;row++){
5894 if (row<0) continue;
5895 if (row>155) continue;
5896 if (ktrack0->GetClusterPointer(row)){
5897 AliTPCTrackerPoint *point =ktrack0->GetTrackPoint(row);
5898 shapesum+=point->GetSigmaY()+point->GetSigmaZ();
5901 if (ktrack1->GetClusterPointer(row)){
5902 AliTPCTrackerPoint *point =ktrack1->GetTrackPoint(row);
5903 shapesum+=point->GetSigmaY()+point->GetSigmaZ();
5908 kink->SetShapeFactor(-1.);
5911 kink->SetShapeFactor(shapesum/sum);
5913 // esd->AddKink(kink);
5915 // kink->SetMother(paramm);
5916 //kink->SetDaughter(paramd);
5918 Double_t chi2P2 = paramm.GetParameter()[2]-paramd.GetParameter()[2];
5920 chi2P2/=paramm.GetCovariance()[5]+paramd.GetCovariance()[5];
5921 Double_t chi2P3 = paramm.GetParameter()[3]-paramd.GetParameter()[3];
5923 chi2P3/=paramm.GetCovariance()[9]+paramd.GetCovariance()[9];
5925 if (AliTPCReconstructor::StreamLevel()>1) {
5926 (*fDebugStreamer)<<"kinkLpt"<<
5934 if ( chi2P2+chi2P3<AliTPCReconstructor::GetRecoParam()->GetKinkAngleCutChi2(0)){
5938 kinks->AddLast(kink);
5944 // sort the kinks according quality - and refit them towards vertex
5946 Int_t nkinks = kinks->GetEntriesFast();
5947 Float_t *quality = new Float_t[nkinks];
5948 Int_t *indexes = new Int_t[nkinks];
5949 AliTPCseed **mothers = new AliTPCseed*[nkinks]; memset(mothers, 0, nkinks*sizeof(AliTPCseed*));
5950 AliTPCseed **daughters = new AliTPCseed*[nkinks]; memset(daughters, 0, nkinks*sizeof(AliTPCseed*));
5953 for (Int_t i=0;i<nkinks;i++){
5955 AliKink *kinkl = (AliKink*)kinks->At(i);
5957 // refit kinks towards vertex
5959 Int_t index0 = kinkl->GetIndex(0);
5960 Int_t index1 = kinkl->GetIndex(1);
5961 AliTPCseed * ktrack0 = (AliTPCseed*)array->At(index0);
5962 AliTPCseed * ktrack1 = (AliTPCseed*)array->At(index1);
5964 Int_t sumn=ktrack0->GetNumberOfClusters()+ktrack1->GetNumberOfClusters();
5966 // Refit Kink under if too small angle
5968 if (kinkl->GetAngle(2)<0.05){
5969 kinkl->SetTPCRow0(GetRowNumber(kinkl->GetR()));
5970 Int_t row0 = kinkl->GetTPCRow0();
5971 Int_t drow = Int_t(2.+0.5/(0.05+kinkl->GetAngle(2)));
5974 Int_t last = row0-drow;
5975 if (last<40) last=40;
5976 if (last<ktrack0->GetFirstPoint()+25) last = ktrack0->GetFirstPoint()+25;
5977 AliTPCseed* seed0 = ReSeed(ktrack0,last,kFALSE);
5980 Int_t first = row0+drow;
5981 if (first>130) first=130;
5982 if (first>ktrack1->GetLastPoint()-25) first = TMath::Max(ktrack1->GetLastPoint()-25,30);
5983 AliTPCseed* seed1 = ReSeed(ktrack1,first,kTRUE);
5985 if (seed0 && seed1){
5986 kinkl->SetStatus(1,8);
5987 if (RefitKink(*seed0,*seed1,*kinkl)) kinkl->SetStatus(1,9);
5988 row0 = GetRowNumber(kinkl->GetR());
5989 sumn = seed0->GetNumberOfClusters()+seed1->GetNumberOfClusters();
5990 mothers[i] = new ( NextFreeSeed() ) AliTPCseed(*seed0);
5991 mothers[i]->SetPoolID(fLastSeedID);
5992 daughters[i] = new (NextFreeSeed() ) AliTPCseed(*seed1);
5993 daughters[i]->SetPoolID(fLastSeedID);
5996 delete kinks->RemoveAt(i);
5997 if (seed0) MarkSeedFree( seed0 );
5998 if (seed1) MarkSeedFree( seed1 );
6001 if (kinkl->GetDistance()>0.5 || kinkl->GetR()<110 || kinkl->GetR()>240) {
6002 delete kinks->RemoveAt(i);
6003 if (seed0) MarkSeedFree( seed0 );
6004 if (seed1) MarkSeedFree( seed1 );
6008 MarkSeedFree( seed0 );
6009 MarkSeedFree( seed1 );
6012 if (kinkl) quality[i] = 160*((0.1+kinkl->GetDistance())*(2.-kinkl->GetTPCDensityFactor()))/(sumn+40.); //the longest -clossest will win
6014 TMath::Sort(nkinks,quality,indexes,kFALSE);
6016 //remove double find kinks
6018 for (Int_t ikink0=1;ikink0<nkinks;ikink0++){
6019 AliKink * kink0 = (AliKink*) kinks->At(indexes[ikink0]);
6020 if (!kink0) continue;
6022 for (Int_t ikink1=0;ikink1<ikink0;ikink1++){
6023 kink0 = (AliKink*) kinks->At(indexes[ikink0]);
6024 if (!kink0) continue;
6025 AliKink * kink1 = (AliKink*) kinks->At(indexes[ikink1]);
6026 if (!kink1) continue;
6027 // if not close kink continue
6028 if (TMath::Abs(kink1->GetPosition()[2]-kink0->GetPosition()[2])>10) continue;
6029 if (TMath::Abs(kink1->GetPosition()[1]-kink0->GetPosition()[1])>10) continue;
6030 if (TMath::Abs(kink1->GetPosition()[0]-kink0->GetPosition()[0])>10) continue;
6032 AliTPCseed &mother0 = *mothers[indexes[ikink0]];
6033 AliTPCseed &daughter0 = *daughters[indexes[ikink0]];
6034 AliTPCseed &mother1 = *mothers[indexes[ikink1]];
6035 AliTPCseed &daughter1 = *daughters[indexes[ikink1]];
6036 Int_t row0 = (kink0->GetTPCRow0()+kink1->GetTPCRow0())/2;
6045 for (Int_t i=0;i<row0;i++){
6046 if (mother0.GetClusterIndex(i)>0 && mother1.GetClusterIndex(i)>0){
6049 if (mother0.GetClusterIndex(i)==mother1.GetClusterIndex(i)){
6056 for (Int_t i=row0;i<158;i++){
6057 //if (daughter0.GetClusterIndex(i)>0 && daughter0.GetClusterIndex(i)>0){ // RS: Bug?
6058 if (daughter0.GetClusterIndex(i)>0 && daughter1.GetClusterIndex(i)>0){
6061 if (mother0.GetClusterIndex(i)==mother1.GetClusterIndex(i)){
6067 Float_t ratio = Float_t(same+1)/Float_t(both+1);
6068 Float_t ratiom = Float_t(samem+1)/Float_t(bothm+1);
6069 Float_t ratiod = Float_t(samed+1)/Float_t(bothd+1);
6070 if (ratio>0.3 && ratiom>0.5 &&ratiod>0.5) {
6071 Int_t sum0 = mother0.GetNumberOfClusters()+daughter0.GetNumberOfClusters();
6072 Int_t sum1 = mother1.GetNumberOfClusters()+daughter1.GetNumberOfClusters();
6074 shared[kink0->GetIndex(0)]= kTRUE;
6075 shared[kink0->GetIndex(1)]= kTRUE;
6076 delete kinks->RemoveAt(indexes[ikink0]);
6080 shared[kink1->GetIndex(0)]= kTRUE;
6081 shared[kink1->GetIndex(1)]= kTRUE;
6082 delete kinks->RemoveAt(indexes[ikink1]);
6089 for (Int_t i=0;i<nkinks;i++){
6090 AliKink * kinkl = (AliKink*) kinks->At(indexes[i]);
6091 if (!kinkl) continue;
6092 kinkl->SetTPCRow0(GetRowNumber(kinkl->GetR()));
6093 Int_t index0 = kinkl->GetIndex(0);
6094 Int_t index1 = kinkl->GetIndex(1);
6095 if (circular[index0]||(circular[index1]&&kinkl->GetDistance()>0.2)) continue;
6096 kinkl->SetMultiple(usage[index0],0);
6097 kinkl->SetMultiple(usage[index1],1);
6098 if (kinkl->GetMultiple()[0]+kinkl->GetMultiple()[1]>2) continue;
6099 if (kinkl->GetMultiple()[0]+kinkl->GetMultiple()[1]>0 && quality[indexes[i]]>0.2) continue;
6100 if (kinkl->GetMultiple()[0]+kinkl->GetMultiple()[1]>0 && kinkl->GetDistance()>0.2) continue;
6101 if (circular[index0]||(circular[index1]&&kinkl->GetDistance()>0.1)) continue;
6103 AliTPCseed * ktrack0 = (AliTPCseed*)array->At(index0);
6104 AliTPCseed * ktrack1 = (AliTPCseed*)array->At(index1);
6105 if (!ktrack0 || !ktrack1) continue;
6106 Int_t index = esd->AddKink(kinkl);
6109 if ( ktrack0->GetKinkIndex(0)==0 && ktrack1->GetKinkIndex(0)==0) { //best kink
6110 if (mothers[indexes[i]]->GetNumberOfClusters()>20 && daughters[indexes[i]]->GetNumberOfClusters()>20 &&
6111 (mothers[indexes[i]]->GetNumberOfClusters()+daughters[indexes[i]]->GetNumberOfClusters())>100){
6112 *ktrack0 = *mothers[indexes[i]];
6113 *ktrack1 = *daughters[indexes[i]];
6117 ktrack0->SetKinkIndex(usage[index0],-(index+1));
6118 ktrack1->SetKinkIndex(usage[index1], (index+1));
6123 // Remove tracks corresponding to shared kink's
6125 for (Int_t i=0;i<nentries;i++){
6126 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
6127 if (!track0) continue;
6128 if (track0->GetKinkIndex(0)!=0) continue;
6129 if (shared[i]) MarkSeedFree( array->RemoveAt(i) );
6134 RemoveUsed2(array,0.5,0.4,30);
6136 for (Int_t i=0;i<nentries;i++){
6137 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
6138 if (!track0) continue;
6139 track0->CookdEdx(0.02,0.6);
6143 for (Int_t i=0;i<nentries;i++){
6144 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
6145 if (!track0) continue;
6146 if (track0->Pt()<1.4) continue;
6147 //remove double high momenta tracks - overlapped with kink candidates
6150 for (Int_t icl=track0->GetFirstPoint();icl<track0->GetLastPoint(); icl++){
6151 if (track0->GetClusterPointer(icl)!=0){
6153 if (track0->GetClusterPointer(icl)->IsUsed(10)) ishared++;
6156 if (Float_t(ishared+1)/Float_t(all+1)>0.5) {
6157 MarkSeedFree( array->RemoveAt(i) );
6161 if (track0->GetKinkIndex(0)!=0) continue;
6162 if (track0->GetNumberOfClusters()<80) continue;
6164 AliTPCseed *pmother = new( NextFreeSeed() ) AliTPCseed();
6165 pmother->SetPoolID(fLastSeedID);
6166 AliTPCseed *pdaughter = new( NextFreeSeed() ) AliTPCseed();
6167 pdaughter->SetPoolID(fLastSeedID);
6168 AliKink *pkink = new AliKink;
6170 AliTPCseed & mother = *pmother;
6171 AliTPCseed & daughter = *pdaughter;
6172 AliKink & kinkl = *pkink;
6173 if (CheckKinkPoint(track0,mother,daughter, kinkl)){
6174 if (mother.GetNumberOfClusters()<30||daughter.GetNumberOfClusters()<20) {
6175 MarkSeedFree( pmother );
6176 MarkSeedFree( pdaughter );
6178 continue; //too short tracks
6180 if (mother.Pt()<1.4) {
6181 MarkSeedFree( pmother );
6182 MarkSeedFree( pdaughter );
6186 Int_t row0= kinkl.GetTPCRow0();
6187 if (kinkl.GetDistance()>0.5 || kinkl.GetR()<110. || kinkl.GetR()>240.) {
6188 MarkSeedFree( pmother );
6189 MarkSeedFree( pdaughter );
6194 Int_t index = esd->AddKink(&kinkl);
6195 mother.SetKinkIndex(0,-(index+1));
6196 daughter.SetKinkIndex(0,index+1);
6197 if (mother.GetNumberOfClusters()>50) {
6198 MarkSeedFree( array->RemoveAt(i) );
6199 AliTPCseed* mtc = new( NextFreeSeed() ) AliTPCseed(mother);
6200 mtc->SetPoolID(fLastSeedID);
6201 array->AddAt(mtc,i);
6204 AliTPCseed* mtc = new( NextFreeSeed() ) AliTPCseed(mother);
6205 mtc->SetPoolID(fLastSeedID);
6206 array->AddLast(mtc);
6208 AliTPCseed* dtc = new( NextFreeSeed() ) AliTPCseed(daughter);
6209 dtc->SetPoolID(fLastSeedID);
6210 array->AddLast(dtc);
6211 for (Int_t icl=0;icl<row0;icl++) {
6212 if (mother.GetClusterPointer(icl)) mother.GetClusterPointer(icl)->Use(20);
6215 for (Int_t icl=row0;icl<158;icl++) {
6216 if (daughter.GetClusterPointer(icl)) daughter.GetClusterPointer(icl)->Use(20);
6220 MarkSeedFree( pmother );
6221 MarkSeedFree( pdaughter );
6225 delete [] daughters;
6247 AliInfo(Form("Ncandidates=\t%d\t%d\t%d\t%d\n",esd->GetNumberOfKinks(),ncandidates,ntracks,nall));
6252 Int_t AliTPCtrackerMI::RefitKink(AliTPCseed &mother, AliTPCseed &daughter, const AliESDkink &knk)
6255 // refit kink towards to the vertex
6258 AliKink &kink=(AliKink &)knk;
6260 Int_t row0 = GetRowNumber(kink.GetR());
6261 FollowProlongation(mother,0);
6262 mother.Reset(kFALSE);
6264 FollowProlongation(daughter,row0);
6265 daughter.Reset(kFALSE);
6266 FollowBackProlongation(daughter,158);
6267 daughter.Reset(kFALSE);
6268 Int_t first = TMath::Max(row0-20,30);
6269 Int_t last = TMath::Min(row0+20,140);
6271 const Int_t kNdiv =5;
6272 AliTPCseed param0[kNdiv]; // parameters along the track
6273 AliTPCseed param1[kNdiv]; // parameters along the track
6274 AliKink kinks[kNdiv]; // corresponding kink parameters
6277 for (Int_t irow=0; irow<kNdiv;irow++){
6278 rows[irow] = first +((last-first)*irow)/(kNdiv-1);
6280 // store parameters along the track
6282 for (Int_t irow=0;irow<kNdiv;irow++){
6283 FollowBackProlongation(mother, rows[irow]);
6284 FollowProlongation(daughter,rows[kNdiv-1-irow]);
6285 param0[irow] = mother;
6286 param1[kNdiv-1-irow] = daughter;
6290 for (Int_t irow=0; irow<kNdiv-1;irow++){
6291 if (param0[irow].GetNumberOfClusters()<kNdiv||param1[irow].GetNumberOfClusters()<kNdiv) continue;
6292 kinks[irow].SetMother(param0[irow]);
6293 kinks[irow].SetDaughter(param1[irow]);
6294 kinks[irow].Update();
6297 // choose kink with best "quality"
6299 Double_t mindist = 10000;
6300 for (Int_t irow=0;irow<kNdiv;irow++){
6301 if (param0[irow].GetNumberOfClusters()<20||param1[irow].GetNumberOfClusters()<20) continue;
6302 if (TMath::Abs(kinks[irow].GetR())>240.) continue;
6303 if (TMath::Abs(kinks[irow].GetR())<100.) continue;
6305 Float_t normdist = TMath::Abs(param0[irow].GetX()-kinks[irow].GetR())*(0.1+kink.GetDistance());
6306 normdist/= (param0[irow].GetNumberOfClusters()+param1[irow].GetNumberOfClusters()+40.);
6307 if (normdist < mindist){
6313 if (index==-1) return 0;
6316 param0[index].Reset(kTRUE);
6317 FollowProlongation(param0[index],0);
6319 mother = param0[index];
6320 daughter = param1[index]; // daughter in vertex
6322 kink.SetMother(mother);
6323 kink.SetDaughter(daughter);
6325 kink.SetTPCRow0(GetRowNumber(kink.GetR()));
6326 kink.SetTPCncls(param0[index].GetNumberOfClusters(),0);
6327 kink.SetTPCncls(param1[index].GetNumberOfClusters(),1);
6328 kink.SetLabel(CookLabel(&mother,0.4, 0,kink.GetTPCRow0()),0);
6329 kink.SetLabel(CookLabel(&daughter,0.4, kink.GetTPCRow0(),160),1);
6330 mother.SetLabel(kink.GetLabel(0));
6331 daughter.SetLabel(kink.GetLabel(1));
6337 void AliTPCtrackerMI::UpdateKinkQualityM(AliTPCseed * seed){
6339 // update Kink quality information for mother after back propagation
6341 if (seed->GetKinkIndex(0)>=0) return;
6342 for (Int_t ikink=0;ikink<3;ikink++){
6343 Int_t index = seed->GetKinkIndex(ikink);
6344 if (index>=0) break;
6345 index = TMath::Abs(index)-1;
6346 AliESDkink * kink = fEvent->GetKink(index);
6347 kink->SetTPCDensity(-1,0,0);
6348 kink->SetTPCDensity(1,0,1);
6350 Int_t row0 = kink->GetTPCRow0() - 2 - Int_t( 0.5/ (0.05+kink->GetAngle(2)));
6351 if (row0<15) row0=15;
6353 Int_t row1 = kink->GetTPCRow0() + 2 + Int_t( 0.5/ (0.05+kink->GetAngle(2)));
6354 if (row1>145) row1=145;
6356 Int_t found,foundable,shared;
6357 seed->GetClusterStatistic(0,row0, found, foundable,shared,kFALSE);
6358 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),0,0);
6359 seed->GetClusterStatistic(row1,155, found, foundable,shared,kFALSE);
6360 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),0,1);
6365 void AliTPCtrackerMI::UpdateKinkQualityD(AliTPCseed * seed){
6367 // update Kink quality information for daughter after refit
6369 if (seed->GetKinkIndex(0)<=0) return;
6370 for (Int_t ikink=0;ikink<3;ikink++){
6371 Int_t index = seed->GetKinkIndex(ikink);
6372 if (index<=0) break;
6373 index = TMath::Abs(index)-1;
6374 AliESDkink * kink = fEvent->GetKink(index);
6375 kink->SetTPCDensity(-1,1,0);
6376 kink->SetTPCDensity(-1,1,1);
6378 Int_t row0 = kink->GetTPCRow0() -2 - Int_t( 0.5/ (0.05+kink->GetAngle(2)));
6379 if (row0<15) row0=15;
6381 Int_t row1 = kink->GetTPCRow0() +2 + Int_t( 0.5/ (0.05+kink->GetAngle(2)));
6382 if (row1>145) row1=145;
6384 Int_t found,foundable,shared;
6385 seed->GetClusterStatistic(0,row0, found, foundable,shared,kFALSE);
6386 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),1,0);
6387 seed->GetClusterStatistic(row1,155, found, foundable,shared,kFALSE);
6388 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),1,1);
6394 Int_t AliTPCtrackerMI::CheckKinkPoint(AliTPCseed*seed,AliTPCseed &mother, AliTPCseed &daughter, const AliESDkink &knk)
6397 // check kink point for given track
6398 // if return value=0 kink point not found
6399 // otherwise seed0 correspond to mother particle
6400 // seed1 correspond to daughter particle
6401 // kink parameter of kink point
6402 AliKink &kink=(AliKink &)knk;
6404 Int_t middlerow = (seed->GetFirstPoint()+seed->GetLastPoint())/2;
6405 Int_t first = seed->GetFirstPoint();
6406 Int_t last = seed->GetLastPoint();
6407 if (last-first<20) return 0; // shortest length - 2*30 = 60 pad-rows
6410 AliTPCseed *seed1 = ReSeed(seed,middlerow+20, kTRUE); //middle of chamber
6411 if (!seed1) return 0;
6412 FollowProlongation(*seed1,seed->GetLastPoint()-20);
6413 seed1->Reset(kTRUE);
6414 FollowProlongation(*seed1,158);
6415 seed1->Reset(kTRUE);
6416 last = seed1->GetLastPoint();
6418 AliTPCseed *seed0 = new( NextFreeSeed() ) AliTPCseed(*seed);
6419 seed0->SetPoolID(fLastSeedID);
6420 seed0->Reset(kFALSE);
6423 AliTPCseed param0[20]; // parameters along the track
6424 AliTPCseed param1[20]; // parameters along the track
6425 AliKink kinks[20]; // corresponding kink parameters
6427 for (Int_t irow=0; irow<20;irow++){
6428 rows[irow] = first +((last-first)*irow)/19;
6430 // store parameters along the track
6432 for (Int_t irow=0;irow<20;irow++){
6433 FollowBackProlongation(*seed0, rows[irow]);
6434 FollowProlongation(*seed1,rows[19-irow]);
6435 param0[irow] = *seed0;
6436 param1[19-irow] = *seed1;
6440 for (Int_t irow=0; irow<19;irow++){
6441 kinks[irow].SetMother(param0[irow]);
6442 kinks[irow].SetDaughter(param1[irow]);
6443 kinks[irow].Update();
6446 // choose kink with biggest change of angle
6448 Double_t maxchange= 0;
6449 for (Int_t irow=1;irow<19;irow++){
6450 if (TMath::Abs(kinks[irow].GetR())>240.) continue;
6451 if (TMath::Abs(kinks[irow].GetR())<110.) continue;
6452 Float_t quality = TMath::Abs(kinks[irow].GetAngle(2))/(3.+TMath::Abs(kinks[irow].GetR()-param0[irow].GetX()));
6453 if ( quality > maxchange){
6454 maxchange = quality;
6459 MarkSeedFree( seed0 );
6460 MarkSeedFree( seed1 );
6461 if (index<0) return 0;
6463 Int_t row0 = GetRowNumber(kinks[index].GetR()); //row 0 estimate
6464 seed0 = new( NextFreeSeed() ) AliTPCseed(param0[index]);
6465 seed0->SetPoolID(fLastSeedID);
6466 seed1 = new( NextFreeSeed() ) AliTPCseed(param1[index]);
6467 seed1->SetPoolID(fLastSeedID);
6468 seed0->Reset(kFALSE);
6469 seed1->Reset(kFALSE);
6470 seed0->ResetCovariance(10.);
6471 seed1->ResetCovariance(10.);
6472 FollowProlongation(*seed0,0);
6473 FollowBackProlongation(*seed1,158);
6474 mother = *seed0; // backup mother at position 0
6475 seed0->Reset(kFALSE);
6476 seed1->Reset(kFALSE);
6477 seed0->ResetCovariance(10.);
6478 seed1->ResetCovariance(10.);
6480 first = TMath::Max(row0-20,0);
6481 last = TMath::Min(row0+20,158);
6483 for (Int_t irow=0; irow<20;irow++){
6484 rows[irow] = first +((last-first)*irow)/19;
6486 // store parameters along the track
6488 for (Int_t irow=0;irow<20;irow++){
6489 FollowBackProlongation(*seed0, rows[irow]);
6490 FollowProlongation(*seed1,rows[19-irow]);
6491 param0[irow] = *seed0;
6492 param1[19-irow] = *seed1;
6496 for (Int_t irow=0; irow<19;irow++){
6497 kinks[irow].SetMother(param0[irow]);
6498 kinks[irow].SetDaughter(param1[irow]);
6499 // param0[irow].Dump();
6500 //param1[irow].Dump();
6501 kinks[irow].Update();
6504 // choose kink with biggest change of angle
6507 for (Int_t irow=0;irow<20;irow++){
6508 if (TMath::Abs(kinks[irow].GetR())>250.) continue;
6509 if (TMath::Abs(kinks[irow].GetR())<90.) continue;
6510 Float_t quality = TMath::Abs(kinks[irow].GetAngle(2))/(3.+TMath::Abs(kinks[irow].GetR()-param0[irow].GetX()));
6511 if ( quality > maxchange){
6512 maxchange = quality;
6519 if (index==-1 || param0[index].GetNumberOfClusters()+param1[index].GetNumberOfClusters()<100){
6520 MarkSeedFree( seed0 );
6521 MarkSeedFree( seed1 );
6525 // Float_t anglesigma = TMath::Sqrt(param0[index].fC22+param0[index].fC33+param1[index].fC22+param1[index].fC33);
6527 kink.SetMother(param0[index]);
6528 kink.SetDaughter(param1[index]);
6531 Double_t chi2P2 = param0[index].GetParameter()[2]-param1[index].GetParameter()[2];
6533 chi2P2/=param0[index].GetCovariance()[5]+param1[index].GetCovariance()[5];
6534 Double_t chi2P3 = param0[index].GetParameter()[3]-param1[index].GetParameter()[3];
6536 chi2P3/=param0[index].GetCovariance()[9]+param1[index].GetCovariance()[9];
6538 if (AliTPCReconstructor::StreamLevel()>1) {
6539 (*fDebugStreamer)<<"kinkHpt"<<
6542 "p0.="<<¶m0[index]<<
6543 "p1.="<<¶m1[index]<<
6547 if ( chi2P2+chi2P3<AliTPCReconstructor::GetRecoParam()->GetKinkAngleCutChi2(0)){
6548 MarkSeedFree( seed0 );
6549 MarkSeedFree( seed1 );
6554 row0 = GetRowNumber(kink.GetR());
6555 kink.SetTPCRow0(row0);
6556 kink.SetLabel(CookLabel(seed0,0.5,0,row0),0);
6557 kink.SetLabel(CookLabel(seed1,0.5,row0,158),1);
6558 kink.SetIndex(-10,0);
6559 kink.SetIndex(int(param0[index].GetNumberOfClusters()+param1[index].GetNumberOfClusters()),1);
6560 kink.SetTPCncls(param0[index].GetNumberOfClusters(),0);
6561 kink.SetTPCncls(param1[index].GetNumberOfClusters(),1);
6564 // new (&mother) AliTPCseed(param0[index]);
6565 daughter = param1[index];
6566 daughter.SetLabel(kink.GetLabel(1));
6567 param0[index].Reset(kTRUE);
6568 FollowProlongation(param0[index],0);
6569 mother = param0[index];
6570 mother.SetLabel(kink.GetLabel(0));
6571 if ( chi2P2+chi2P3<AliTPCReconstructor::GetRecoParam()->GetKinkAngleCutChi2(1)){
6574 MarkSeedFree( seed0 );
6575 MarkSeedFree( seed1 );
6583 AliTPCseed* AliTPCtrackerMI::ReSeed(AliTPCseed *t)
6586 // reseed - refit - track
6589 // Int_t last = fSectors->GetNRows()-1;
6591 if (fSectors == fOuterSec){
6592 first = TMath::Max(first, t->GetFirstPoint()-fInnerSec->GetNRows());
6596 first = t->GetFirstPoint();
6598 AliTPCseed * seed = MakeSeed(t,0.1,0.5,0.9);
6599 FollowBackProlongation(*t,fSectors->GetNRows()-1);
6601 FollowProlongation(*t,first);
6611 //_____________________________________________________________________________
6612 Int_t AliTPCtrackerMI::ReadSeeds(const TFile *inp) {
6613 //-----------------------------------------------------------------
6614 // This function reades track seeds.
6615 //-----------------------------------------------------------------
6616 TDirectory *savedir=gDirectory;
6618 TFile *in=(TFile*)inp;
6619 if (!in->IsOpen()) {
6620 cerr<<"AliTPCtrackerMI::ReadSeeds(): input file is not open !\n";
6625 TTree *seedTree=(TTree*)in->Get("Seeds");
6627 cerr<<"AliTPCtrackerMI::ReadSeeds(): ";
6628 cerr<<"can't get a tree with track seeds !\n";
6631 AliTPCtrack *seed=new AliTPCtrack;
6632 seedTree->SetBranchAddress("tracks",&seed);
6634 if (fSeeds==0) fSeeds=new TObjArray(15000);
6636 Int_t n=(Int_t)seedTree->GetEntries();
6637 for (Int_t i=0; i<n; i++) {
6638 seedTree->GetEvent(i);
6639 AliTPCseed* sdc = new( NextFreeSeed() ) AliTPCseed(*seed/*,seed->GetAlpha()*/);
6640 sdc->SetPoolID(fLastSeedID);
6641 fSeeds->AddLast(sdc);
6644 delete seed; // RS: this seed is not from the pool, delete it !!!
6650 Int_t AliTPCtrackerMI::Clusters2Tracks (AliESDEvent *const esd)
6653 // clusters to tracks
6654 if (fSeeds) DeleteSeeds();
6655 else ResetSeedsPool();
6658 AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ;
6659 transform->SetCurrentTimeStamp( esd->GetTimeStamp());
6660 transform->SetCurrentRun(esd->GetRunNumber());
6663 if (!fSeeds) return 1;
6665 if (AliTPCReconstructor::StreamLevel()>3) DumpClusters(0,fSeeds);
6671 //_____________________________________________________________________________
6672 Int_t AliTPCtrackerMI::Clusters2Tracks() {
6673 //-----------------------------------------------------------------
6674 // This is a track finder.
6675 //-----------------------------------------------------------------
6676 TDirectory *savedir=gDirectory;
6680 fSeeds = Tracking();
6683 Info("Clusters2Tracks","Time for tracking: \t");timer.Print();timer.Start();
6685 //activate again some tracks
6686 for (Int_t i=0; i<fSeeds->GetEntriesFast(); i++) {
6687 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
6689 Int_t nc=t.GetNumberOfClusters();
6691 MarkSeedFree( fSeeds->RemoveAt(i) );
6695 if (pt->GetRemoval()==10) {
6696 if (pt->GetDensityFirst(20)>0.8 || pt->GetDensityFirst(30)>0.8 || pt->GetDensityFirst(40)>0.7)
6697 pt->Desactivate(10); // make track again active // MvL: should be 0 ?
6699 pt->Desactivate(20);
6700 MarkSeedFree( fSeeds->RemoveAt(i) );
6705 RemoveUsed2(fSeeds,0.85,0.85,0);
6706 if (AliTPCReconstructor::GetRecoParam()->GetDoKinks()) FindKinks(fSeeds,fEvent);
6707 //FindCurling(fSeeds, fEvent,0);
6708 if (AliTPCReconstructor::StreamLevel()>5) FindMultiMC(fSeeds, fEvent,-1); // find multi found tracks
6709 RemoveUsed2(fSeeds,0.5,0.4,20);
6710 FindSplitted(fSeeds, fEvent,0); // find multi found tracks
6711 if (AliTPCReconstructor::StreamLevel()>5) FindMultiMC(fSeeds, fEvent,0); // find multi found tracks
6714 // // refit short tracks
6716 Int_t nseed=fSeeds->GetEntriesFast();
6719 for (Int_t i=0; i<nseed; i++) {
6720 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
6722 Int_t nc=t.GetNumberOfClusters();
6724 MarkSeedFree( fSeeds->RemoveAt(i) );
6727 CookLabel(pt,0.1); //For comparison only
6728 //if ((pt->IsActive() || (pt->fRemoval==10) )&& nc>50 &&pt->GetNumberOfClusters()>0.4*pt->fNFoundable){
6729 if ((pt->IsActive() || (pt->GetRemoval()==10) )){
6731 if (fDebug>0) cerr<<found<<'\r';
6735 MarkSeedFree( fSeeds->RemoveAt(i) );
6739 //RemoveOverlap(fSeeds,0.99,7,kTRUE);
6741 //RemoveUsed(fSeeds,0.9,0.9,6);
6743 nseed=fSeeds->GetEntriesFast();
6745 for (Int_t i=0; i<nseed; i++) {
6746 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
6748 Int_t nc=t.GetNumberOfClusters();
6750 MarkSeedFree( fSeeds->RemoveAt(i) );
6754 t.CookdEdx(0.02,0.6);
6755 // CheckKinkPoint(&t,0.05);
6756 //if ((pt->IsActive() || (pt->fRemoval==10) )&& nc>50 &&pt->GetNumberOfClusters()>0.4*pt->fNFoundable){
6757 if ((pt->IsActive() || (pt->GetRemoval()==10) )){
6765 MarkSeedFree( fSeeds->RemoveAt(i) );
6766 //AliTPCseed * seed1 = ReSeed(pt,0.05,0.5,1);
6768 // FollowProlongation(*seed1,0);
6769 // Int_t n = seed1->GetNumberOfClusters();
6770 // printf("fP4\t%f\t%f\n",seed1->GetC(),pt->GetC());
6771 // printf("fN\t%d\t%d\n", seed1->GetNumberOfClusters(),pt->GetNumberOfClusters());
6774 //AliTPCseed * seed2 = ReSeed(pt,0.95,0.5,0.05);
6778 SortTracks(fSeeds, 1);
6782 PrepareForBackProlongation(fSeeds,5.);
6783 PropagateBack(fSeeds);
6784 printf("Time for back propagation: \t");timer.Print();timer.Start();
6788 PrepareForProlongation(fSeeds,5.);
6789 PropagateForard2(fSeeds);
6791 printf("Time for FORWARD propagation: \t");timer.Print();timer.Start();
6792 // RemoveUsed(fSeeds,0.7,0.7,6);
6793 //RemoveOverlap(fSeeds,0.9,7,kTRUE);
6795 nseed=fSeeds->GetEntriesFast();
6797 for (Int_t i=0; i<nseed; i++) {
6798 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
6800 Int_t nc=t.GetNumberOfClusters();
6802 MarkSeedFree( fSeeds->RemoveAt(i) );
6805 t.CookdEdx(0.02,0.6);
6806 // CookLabel(pt,0.1); //For comparison only
6807 //if ((pt->IsActive() || (pt->fRemoval==10) )&& nc>50 &&pt->GetNumberOfClusters()>0.4*pt->fNFoundable){
6808 if ((pt->IsActive() || (pt->fRemoval==10) )){
6809 cerr<<found++<<'\r';
6812 MarkSeedFree( fSeeds->RemoveAt(i) );
6817 // fNTracks = found;
6819 Info("Clusters2Tracks","Time for overlap removal, track writing and dedx cooking: \t"); timer.Print();timer.Start();
6822 // cerr<<"Number of found tracks : "<<"\t"<<found<<endl;
6823 Info("Clusters2Tracks","Number of found tracks %d",found);
6825 // UnloadClusters();
6830 void AliTPCtrackerMI::Tracking(TObjArray * arr)
6833 // tracking of the seeds
6836 fSectors = fOuterSec;
6837 ParallelTracking(arr,150,63);
6838 fSectors = fOuterSec;
6839 ParallelTracking(arr,63,0);
6842 TObjArray * AliTPCtrackerMI::Tracking(Int_t seedtype, Int_t i1, Int_t i2, Float_t cuts[4], Float_t dy, Int_t dsec)
6847 static TObjArray arrTracks;
6848 TObjArray * arr = &arrTracks;
6850 fSectors = fOuterSec;
6853 for (Int_t sec=0;sec<fkNOS;sec++){
6854 if (seedtype==3) MakeSeeds3(arr,sec,i1,i2,cuts,dy, dsec);
6855 if (seedtype==4) MakeSeeds5(arr,sec,i1,i2,cuts,dy);
6856 if (seedtype==2) MakeSeeds2(arr,sec,i1,i2,cuts,dy);
6859 Info("Tracking","\nSeeding - %d\t%d\t%d\t%d\n",seedtype,i1,i2,arr->GetEntriesFast());
6871 TObjArray * AliTPCtrackerMI::Tracking()
6875 if (AliTPCReconstructor::GetRecoParam()->GetSpecialSeeding()) return TrackingSpecial();
6878 Int_t nup=fOuterSec->GetNRows()+fInnerSec->GetNRows();
6880 TObjArray * seeds = new TObjArray;
6882 Int_t fLastSeedRowSec=AliTPCReconstructor::GetRecoParam()->GetLastSeedRowSec();
6883 Int_t gapPrim = AliTPCReconstructor::GetRecoParam()->GetSeedGapPrim();
6884 Int_t gapSec = AliTPCReconstructor::GetRecoParam()->GetSeedGapSec();
6892 Float_t fnumber = 3.0;
6893 Float_t fdensity = 3.0;
6898 for (Int_t delta = 0; delta<18; delta+=gapPrim){
6902 arr = Tracking(3,nup-1-delta,nup-1-delta-gap,cuts,-1,1);
6903 SumTracks(seeds,arr);
6904 SignClusters(seeds,fnumber,fdensity);
6906 for (Int_t i=2;i<6;i+=2){
6907 // seed high pt tracks
6910 arr = Tracking(3,nup-i-delta,nup-i-delta-gap,cuts,-1,0);
6911 SumTracks(seeds,arr);
6912 SignClusters(seeds,fnumber,fdensity);
6917 // RemoveUsed(seeds,0.9,0.9,1);
6918 // UnsignClusters();
6919 // SignClusters(seeds,fnumber,fdensity);
6923 for (Int_t delta = 20; delta<120; delta+=gapPrim){
6925 // seed high pt tracks
6929 arr = Tracking(3,nup-delta,nup-delta-gap,cuts,-1);
6930 SumTracks(seeds,arr);
6931 SignClusters(seeds,fnumber,fdensity);
6936 arr = Tracking(3,nup-delta-5,nup-delta-5-gap,cuts,-1);
6937 SumTracks(seeds,arr);
6938 SignClusters(seeds,fnumber,fdensity);
6949 Info("Tracking()","\n\nPrimary seeding\t%d\n\n",seeds->GetEntriesFast());
6953 // RemoveUsed(seeds,0.75,0.75,1);
6955 //SignClusters(seeds,fnumber,fdensity);
6964 arr = Tracking(4,nup-1,nup-1-gap,cuts,-1);
6965 SumTracks(seeds,arr);
6966 SignClusters(seeds,fnumber,fdensity);
6968 arr = Tracking(4,nup-2,nup-2-gap,cuts,-1);
6969 SumTracks(seeds,arr);
6970 SignClusters(seeds,fnumber,fdensity);
6972 arr = Tracking(4,nup-3,nup-3-gap,cuts,-1);
6973 SumTracks(seeds,arr);
6974 SignClusters(seeds,fnumber,fdensity);
6976 arr = Tracking(4,nup-5,nup-5-gap,cuts,-1);
6977 SumTracks(seeds,arr);
6978 SignClusters(seeds,fnumber,fdensity);
6980 arr = Tracking(4,nup-7,nup-7-gap,cuts,-1);
6981 SumTracks(seeds,arr);
6982 SignClusters(seeds,fnumber,fdensity);
6985 arr = Tracking(4,nup-9,nup-9-gap,cuts,-1);
6986 SumTracks(seeds,arr);
6987 SignClusters(seeds,fnumber,fdensity);
6991 for (Int_t delta = 9; delta<30; delta+=gapSec){
6997 arr = Tracking(4,nup-1-delta,nup-1-delta-gap,cuts,-1);
6998 SumTracks(seeds,arr);
6999 SignClusters(seeds,fnumber,fdensity);
7001 arr = Tracking(4,nup-3-delta,nup-5-delta-gap,cuts,4);
7002 SumTracks(seeds,arr);
7003 SignClusters(seeds,fnumber,fdensity);
7016 for (Int_t delta = 30; delta<fLastSeedRowSec; delta+=gapSec){
7022 arr = Tracking(4,nup-1-delta,nup-1-delta-gap,cuts,-1);
7023 SumTracks(seeds,arr);
7024 SignClusters(seeds,fnumber,fdensity);
7026 arr = Tracking(4,nup-5-delta,nup-5-delta-gap,cuts,5 );
7027 SumTracks(seeds,arr);
7028 SignClusters(seeds,fnumber,fdensity);
7032 Info("Tracking()","\n\nSecondary seeding\t%d\n\n",seeds->GetEntriesFast());
7043 TObjArray * AliTPCtrackerMI::TrackingSpecial()
7046 // seeding adjusted for laser and cosmic tests - short tracks with big inclination angle
7047 // no primary vertex seeding tried
7051 Int_t nup=fOuterSec->GetNRows()+fInnerSec->GetNRows();
7053 TObjArray * seeds = new TObjArray;
7058 Float_t fnumber = 3.0;
7059 Float_t fdensity = 3.0;
7062 cuts[0] = AliTPCReconstructor::GetRecoParam()->GetMaxC(); // max curvature
7063 cuts[1] = 3.5; // max tan(phi) angle for seeding
7064 cuts[2] = 3.; // not used (cut on z primary vertex)
7065 cuts[3] = 3.5; // max tan(theta) angle for seeding
7067 for (Int_t delta = 0; nup-delta-gap-1>0; delta+=3){
7069 arr = Tracking(4,nup-1-delta,nup-1-delta-gap,cuts,-1);
7070 SumTracks(seeds,arr);
7071 SignClusters(seeds,fnumber,fdensity);
7075 Info("Tracking()","\n\nSecondary seeding\t%d\n\n",seeds->GetEntriesFast());
7086 void AliTPCtrackerMI::SumTracks(TObjArray *arr1,TObjArray *&arr2)
7089 //sum tracks to common container
7090 //remove suspicious tracks
7091 // RS: Attention: supplied tracks come in the static array, don't delete them
7092 Int_t nseed = arr2->GetEntriesFast();
7093 for (Int_t i=0;i<nseed;i++){
7094 AliTPCseed *pt=(AliTPCseed*)arr2->UncheckedAt(i);
7097 // remove tracks with too big curvature
7099 if (TMath::Abs(pt->GetC())>AliTPCReconstructor::GetRecoParam()->GetMaxC()){
7100 MarkSeedFree( arr2->RemoveAt(i) );
7103 // REMOVE VERY SHORT TRACKS
7104 if (pt->GetNumberOfClusters()<20){
7105 MarkSeedFree( arr2->RemoveAt(i) );
7108 // NORMAL ACTIVE TRACK
7109 if (pt->IsActive()){
7110 arr1->AddLast(arr2->RemoveAt(i));
7113 //remove not usable tracks
7114 if (pt->GetRemoval()!=10){
7115 MarkSeedFree( arr2->RemoveAt(i) );
7119 // ENABLE ONLY ENOUGH GOOD STOPPED TRACKS
7120 if (pt->GetDensityFirst(20)>0.8 || pt->GetDensityFirst(30)>0.8 || pt->GetDensityFirst(40)>0.7)
7121 arr1->AddLast(arr2->RemoveAt(i));
7123 MarkSeedFree( arr2->RemoveAt(i) );
7127 // delete arr2; arr2 = 0; // RS: this is static array, don't delete it
7132 void AliTPCtrackerMI::ParallelTracking(TObjArray *const arr, Int_t rfirst, Int_t rlast)
7135 // try to track in parralel
7137 Int_t nseed=arr->GetEntriesFast();
7138 //prepare seeds for tracking
7139 for (Int_t i=0; i<nseed; i++) {
7140 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i), &t=*pt;
7142 if (!t.IsActive()) continue;
7143 // follow prolongation to the first layer
7144 if ( (fSectors ==fInnerSec) || (t.GetFirstPoint()-fkParam->GetNRowLow()>rfirst+1) )
7145 FollowProlongation(t, rfirst+1);
7150 for (Int_t nr=rfirst; nr>=rlast; nr--){
7151 if (nr<fInnerSec->GetNRows())
7152 fSectors = fInnerSec;
7154 fSectors = fOuterSec;
7155 // make indexes with the cluster tracks for given
7157 // find nearest cluster
7158 for (Int_t i=0; i<nseed; i++) {
7159 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i), &t=*pt;
7161 if (nr==80) pt->UpdateReference();
7162 if (!pt->IsActive()) continue;
7163 // if ( (fSectors ==fOuterSec) && (pt->fFirstPoint-fkParam->GetNRowLow())<nr) continue;
7164 if (pt->GetRelativeSector()>17) {
7167 UpdateClusters(t,nr);
7169 // prolonagate to the nearest cluster - if founded
7170 for (Int_t i=0; i<nseed; i++) {
7171 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
7173 if (!pt->IsActive()) continue;
7174 // if ((fSectors ==fOuterSec) && (pt->fFirstPoint-fkParam->GetNRowLow())<nr) continue;
7175 if (pt->GetRelativeSector()>17) {
7178 FollowToNextCluster(*pt,nr);
7183 void AliTPCtrackerMI::PrepareForBackProlongation(const TObjArray *const arr,Float_t fac) const
7187 // if we use TPC track itself we have to "update" covariance
7189 Int_t nseed= arr->GetEntriesFast();
7190 for (Int_t i=0;i<nseed;i++){
7191 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
7195 //rotate to current local system at first accepted point
7196 Int_t index = pt->GetClusterIndex2(pt->GetFirstPoint());
7197 Int_t sec = (index&0xff000000)>>24;
7199 Float_t angle1 = fInnerSec->GetAlpha()*sec+fInnerSec->GetAlphaShift();
7200 if (angle1>TMath::Pi())
7201 angle1-=2.*TMath::Pi();
7202 Float_t angle2 = pt->GetAlpha();
7204 if (TMath::Abs(angle1-angle2)>0.001){
7205 if (!pt->Rotate(angle1-angle2)) return;
7206 //angle2 = pt->GetAlpha();
7207 //pt->fRelativeSector = pt->GetAlpha()/fInnerSec->GetAlpha();
7208 //if (pt->GetAlpha()<0)
7209 // pt->fRelativeSector+=18;
7210 //sec = pt->fRelativeSector;
7219 void AliTPCtrackerMI::PrepareForProlongation(TObjArray *const arr, Float_t fac) const
7223 // if we use TPC track itself we have to "update" covariance
7225 Int_t nseed= arr->GetEntriesFast();
7226 for (Int_t i=0;i<nseed;i++){
7227 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
7230 pt->SetFirstPoint(pt->GetLastPoint());
7238 Int_t AliTPCtrackerMI::PropagateBack(const TObjArray *const arr)
7241 // make back propagation
7243 Int_t nseed= arr->GetEntriesFast();
7244 for (Int_t i=0;i<nseed;i++){
7245 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
7246 if (pt&& pt->GetKinkIndex(0)<=0) {
7247 //AliTPCseed *pt2 = new AliTPCseed(*pt);
7248 fSectors = fInnerSec;
7249 //FollowBackProlongation(*pt,fInnerSec->GetNRows()-1);
7250 //fSectors = fOuterSec;
7251 FollowBackProlongation(*pt,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1,1);
7252 //if (pt->GetNumberOfClusters()<(pt->fEsd->GetTPCclusters(0)) ){
7253 // Error("PropagateBack","Not prolonged track %d",pt->GetLabel());
7254 // FollowBackProlongation(*pt2,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1);
7257 if (pt&& pt->GetKinkIndex(0)>0) {
7258 AliESDkink * kink = fEvent->GetKink(pt->GetKinkIndex(0)-1);
7259 pt->SetFirstPoint(kink->GetTPCRow0());
7260 fSectors = fInnerSec;
7261 FollowBackProlongation(*pt,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1,1);
7269 Int_t AliTPCtrackerMI::PropagateForward2(const TObjArray *const arr)
7272 // make forward propagation
7274 Int_t nseed= arr->GetEntriesFast();
7276 for (Int_t i=0;i<nseed;i++){
7277 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
7279 FollowProlongation(*pt,0,1,1);
7288 Int_t AliTPCtrackerMI::PropagateForward()
7291 // propagate track forward
7293 Int_t nseed = fSeeds->GetEntriesFast();
7294 for (Int_t i=0;i<nseed;i++){
7295 AliTPCseed *pt = (AliTPCseed*)fSeeds->UncheckedAt(i);
7297 AliTPCseed &t = *pt;
7298 Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
7299 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
7300 if (alpha < 0. ) alpha += 2.*TMath::Pi();
7301 t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN);
7305 fSectors = fOuterSec;
7306 ParallelTracking(fSeeds,fOuterSec->GetNRows()+fInnerSec->GetNRows()-1,fInnerSec->GetNRows());
7307 fSectors = fInnerSec;
7308 ParallelTracking(fSeeds,fInnerSec->GetNRows()-1,0);
7317 Int_t AliTPCtrackerMI::PropagateBack(AliTPCseed *const pt, Int_t row0, Int_t row1)
7320 // make back propagation, in between row0 and row1
7324 fSectors = fInnerSec;
7327 if (row1<fSectors->GetNRows())
7330 r1 = fSectors->GetNRows()-1;
7332 if (row0<fSectors->GetNRows()&& r1>0 )
7333 FollowBackProlongation(*pt,r1);
7334 if (row1<=fSectors->GetNRows())
7337 r1 = row1 - fSectors->GetNRows();
7338 if (r1<=0) return 0;
7339 if (r1>=fOuterSec->GetNRows()) return 0;
7340 fSectors = fOuterSec;
7341 return FollowBackProlongation(*pt,r1);
7349 void AliTPCtrackerMI::GetShape(AliTPCseed * seed, Int_t row)
7351 // gets cluster shape
7353 AliTPCClusterParam * clparam = AliTPCcalibDB::Instance()->GetClusterParam();
7354 Float_t zdrift = TMath::Abs((fkParam->GetZLength(0)-TMath::Abs(seed->GetZ())));
7355 Int_t type = (seed->GetSector() < fkParam->GetNSector()/2) ? 0: (row>126) ? 1:2;
7356 Double_t angulary = seed->GetSnp();
7358 if (TMath::Abs(angulary)>AliTPCReconstructor::GetMaxSnpTracker()) {
7359 angulary = TMath::Sign(AliTPCReconstructor::GetMaxSnpTracker(),angulary);
7362 angulary = angulary*angulary/((1.-angulary)*(1.+angulary));
7363 Double_t angularz = seed->GetTgl()*seed->GetTgl()*(1.+angulary);
7365 Double_t sigmay = clparam->GetRMS0(0,type,zdrift,TMath::Sqrt(TMath::Abs(angulary)));
7366 Double_t sigmaz = clparam->GetRMS0(1,type,zdrift,TMath::Sqrt(TMath::Abs(angularz)));
7367 seed->SetCurrentSigmaY2(sigmay*sigmay);
7368 seed->SetCurrentSigmaZ2(sigmaz*sigmaz);
7369 // Float_t sd2 = TMath::Abs((fkParam->GetZLength(0)-TMath::Abs(seed->GetZ())))*fkParam->GetDiffL()*fkParam->GetDiffL();
7370 // // Float_t padlength = fkParam->GetPadPitchLength(seed->fSector);
7371 // Float_t padlength = GetPadPitchLength(row);
7373 // Float_t sresy = (seed->GetSector() < fkParam->GetNSector()/2) ? 0.2 :0.3;
7374 // seed->SetCurrentSigmaY2(sd2+padlength*padlength*angulary/12.+sresy*sresy);
7376 // Float_t sresz = fkParam->GetZSigma();
7377 // seed->SetCurrentSigmaZ2(sd2+padlength*padlength*angularz*angularz*(1+angulary)/12.+sresz*sresz);
7379 Float_t wy = GetSigmaY(seed);
7380 Float_t wz = GetSigmaZ(seed);
7383 if (TMath::Abs(wy/seed->fCurrentSigmaY2-1)>0.0001 || TMath::Abs(wz/seed->fCurrentSigmaZ2-1)>0.0001 ){
7384 printf("problem\n");
7391 //__________________________________________________________________________
7392 void AliTPCtrackerMI::CookLabel(AliKalmanTrack *tk, Float_t wrong) const {
7393 //--------------------------------------------------------------------
7394 //This function "cooks" a track label. If label<0, this track is fake.
7395 //--------------------------------------------------------------------
7396 AliTPCseed * t = dynamic_cast<AliTPCseed*>(tk);
7398 printf("%s:%d wrong type \n",(char*)__FILE__,__LINE__);
7402 Int_t noc=t->GetNumberOfClusters();
7404 //printf("\nnot founded prolongation\n\n\n");
7410 AliTPCclusterMI *clusters[160];
7412 for (Int_t i=0;i<160;i++) {
7419 for (i=0; i<160 && current<noc; i++) {
7421 Int_t index=t->GetClusterIndex2(i);
7422 if (index<=0) continue;
7423 if (index&0x8000) continue;
7425 //clusters[current]=GetClusterMI(index);
7426 if (t->GetClusterPointer(i)){
7427 clusters[current]=t->GetClusterPointer(i);
7433 Int_t lab=123456789;
7434 for (i=0; i<noc; i++) {
7435 AliTPCclusterMI *c=clusters[i];
7437 lab=TMath::Abs(c->GetLabel(0));
7439 for (j=0; j<noc; j++) if (lb[j]==lab || mx[j]==0) break;
7445 for (i=0; i<noc; i++) if (mx[i]>max) {max=mx[i]; lab=lb[i];}
7447 for (i=0; i<noc; i++) {
7448 AliTPCclusterMI *c=clusters[i];
7450 if (TMath::Abs(c->GetLabel(1)) == lab ||
7451 TMath::Abs(c->GetLabel(2)) == lab ) max++;
7453 if (noc<=0) { lab=-1; return;}
7454 if ((1.- Float_t(max)/(noc)) > wrong) lab=-lab;
7457 Int_t tail=Int_t(0.10*noc);
7460 for (i=1; i<160&&ind<tail; i++) {
7461 // AliTPCclusterMI *c=clusters[noc-i];
7462 AliTPCclusterMI *c=clusters[i];
7464 if (lab == TMath::Abs(c->GetLabel(0)) ||
7465 lab == TMath::Abs(c->GetLabel(1)) ||
7466 lab == TMath::Abs(c->GetLabel(2))) max++;
7469 if (max < Int_t(0.5*tail)) lab=-lab;
7476 //delete[] clusters;
7480 //__________________________________________________________________________
7481 Int_t AliTPCtrackerMI::CookLabel(AliTPCseed *const t, Float_t wrong,Int_t first, Int_t last) const {
7482 //--------------------------------------------------------------------
7483 //This function "cooks" a track label. If label<0, this track is fake.
7484 //--------------------------------------------------------------------
7485 Int_t noc=t->GetNumberOfClusters();
7487 //printf("\nnot founded prolongation\n\n\n");
7493 AliTPCclusterMI *clusters[160];
7495 for (Int_t i=0;i<160;i++) {
7502 for (i=0; i<160 && current<noc; i++) {
7503 if (i<first) continue;
7504 if (i>last) continue;
7505 Int_t index=t->GetClusterIndex2(i);
7506 if (index<=0) continue;
7507 if (index&0x8000) continue;
7509 //clusters[current]=GetClusterMI(index);
7510 if (t->GetClusterPointer(i)){
7511 clusters[current]=t->GetClusterPointer(i);
7516 //if (noc<5) return -1;
7517 Int_t lab=123456789;
7518 for (i=0; i<noc; i++) {
7519 AliTPCclusterMI *c=clusters[i];
7521 lab=TMath::Abs(c->GetLabel(0));
7523 for (j=0; j<noc; j++) if (lb[j]==lab || mx[j]==0) break;
7529 for (i=0; i<noc; i++) if (mx[i]>max) {max=mx[i]; lab=lb[i];}
7531 for (i=0; i<noc; i++) {
7532 AliTPCclusterMI *c=clusters[i];
7534 if (TMath::Abs(c->GetLabel(1)) == lab ||
7535 TMath::Abs(c->GetLabel(2)) == lab ) max++;
7537 if (noc<=0) { lab=-1; return -1;}
7538 if ((1.- Float_t(max)/(noc)) > wrong) lab=-lab;
7541 Int_t tail=Int_t(0.10*noc);
7544 for (i=1; i<160&&ind<tail; i++) {
7545 // AliTPCclusterMI *c=clusters[noc-i];
7546 AliTPCclusterMI *c=clusters[i];
7548 if (lab == TMath::Abs(c->GetLabel(0)) ||
7549 lab == TMath::Abs(c->GetLabel(1)) ||
7550 lab == TMath::Abs(c->GetLabel(2))) max++;
7553 if (max < Int_t(0.5*tail)) lab=-lab;
7556 // t->SetLabel(lab);
7560 //delete[] clusters;
7564 Int_t AliTPCtrackerMI::GetRowNumber(Double_t x[3]) const
7566 //return pad row number for given x vector
7567 Float_t phi = TMath::ATan2(x[1],x[0]);
7568 if(phi<0) phi=2.*TMath::Pi()+phi;
7569 // Get the local angle in the sector philoc
7570 const Float_t kRaddeg = 180/3.14159265358979312;
7571 Float_t phiangle = (Int_t (phi*kRaddeg/20.) + 0.5)*20./kRaddeg;
7572 Double_t localx = x[0]*TMath::Cos(phiangle)-x[1]*TMath::Sin(phiangle);
7573 return GetRowNumber(localx);
7578 void AliTPCtrackerMI::MakeESDBitmaps(AliTPCseed *t, AliESDtrack *esd)
7580 //-----------------------------------------------------------------------
7581 // Fill the cluster and sharing bitmaps of the track
7582 //-----------------------------------------------------------------------
7584 Int_t firstpoint = 0;
7585 Int_t lastpoint = 159;
7586 AliTPCTrackerPoint *point;
7587 AliTPCclusterMI *cluster;
7590 TBits clusterMap(159);
7591 TBits sharedMap(159);
7593 for (int iter=firstpoint; iter<lastpoint; iter++) {
7594 // Change to cluster pointers to see if we have a cluster at given padrow
7596 cluster = t->GetClusterPointer(iter);
7598 clusterMap.SetBitNumber(iter, kTRUE);
7599 point = t->GetTrackPoint(iter);
7600 if (point->IsShared())
7601 sharedMap.SetBitNumber(iter,kTRUE);
7603 if (t->GetClusterIndex(iter) >= 0 && (t->GetClusterIndex(iter) & 0x8000) == 0) {
7604 fitMap.SetBitNumber(iter, kTRUE);
7608 esd->SetTPCClusterMap(clusterMap);
7609 esd->SetTPCSharedMap(sharedMap);
7610 esd->SetTPCFitMap(fitMap);
7611 if (nclsf != t->GetNumberOfClusters())
7612 AliWarning(Form("Inconsistency between ncls %d and indices %d (found %d)",t->GetNumberOfClusters(),nclsf,esd->GetTPCClusterMap().CountBits()));
7615 Bool_t AliTPCtrackerMI::IsFindable(AliTPCseed & track){
7617 // return flag if there is findable cluster at given position
7620 Float_t z = track.GetZ();
7622 if (TMath::Abs(z)<(AliTPCReconstructor::GetCtgRange()*track.GetX()+kDeltaZ) &&
7623 TMath::Abs(z)<fkParam->GetZLength(0) &&
7624 (TMath::Abs(track.GetSnp())<AliTPCReconstructor::GetMaxSnpTracker()))
7630 void AliTPCtrackerMI::AddCovariance(AliTPCseed * seed){
7632 // Adding systematic error estimate to the covariance matrix
7633 // !!!! the systematic error for element 4 is in 1/GeV
7634 // 03.03.2012 MI changed in respect to the previous versions
7635 const Double_t *param = AliTPCReconstructor::GetRecoParam()->GetSystematicError();
7637 // use only the diagonal part if not specified otherwise
7638 if (!AliTPCReconstructor::GetRecoParam()->GetUseSystematicCorrelation()) return AddCovarianceAdd(seed);
7640 Double_t *covarS= (Double_t*)seed->GetCovariance();
7641 Double_t factor[5]={1,1,1,1,1};
7642 factor[0]= TMath::Sqrt(TMath::Abs((covarS[0] + param[0]*param[0])/covarS[0]));
7643 factor[1]= TMath::Sqrt(TMath::Abs((covarS[2] + param[1]*param[1])/covarS[2]));
7644 factor[2]= TMath::Sqrt(TMath::Abs((covarS[5] + param[2]*param[2])/covarS[5]));
7645 factor[3]= TMath::Sqrt(TMath::Abs((covarS[9] + param[3]*param[3])/covarS[9]));
7646 factor[4]= TMath::Sqrt(TMath::Abs((covarS[14] +param[4]*param[4])/covarS[14]));
7648 factor[0]=factor[2];
7649 factor[4]=factor[2];
7655 for (Int_t i=0; i<5; i++){
7656 for (Int_t j=i; j<5; j++){
7657 Int_t index=seed->GetIndex(i,j);
7658 covarS[index]*=factor[i]*factor[j];
7664 void AliTPCtrackerMI::AddCovarianceAdd(AliTPCseed * seed){
7666 // Adding systematic error - as additive factor without correlation
7668 // !!!! the systematic error for element 4 is in 1/GeV
7669 // 03.03.2012 MI changed in respect to the previous versions
7671 const Double_t *param = AliTPCReconstructor::GetRecoParam()->GetSystematicError();
7672 Double_t *covarIn= (Double_t*)seed->GetCovariance();
7674 for (Int_t i=0;i<15;i++) covar[i]=0;
7680 covar[0] = param[0]*param[0];
7681 covar[2] = param[1]*param[1];
7682 covar[5] = param[2]*param[2];
7683 covar[9] = param[3]*param[3];
7684 covar[14]= param[4]*param[4];
7686 covar[1]=TMath::Sqrt((covar[0]*covar[2]))*covarIn[1]/TMath::Sqrt((covarIn[0]*covarIn[2]));
7688 covar[3]=TMath::Sqrt((covar[0]*covar[5]))*covarIn[3]/TMath::Sqrt((covarIn[0]*covarIn[5]));
7689 covar[4]=TMath::Sqrt((covar[2]*covar[5]))*covarIn[4]/TMath::Sqrt((covarIn[2]*covarIn[5]));
7691 covar[6]=TMath::Sqrt((covar[0]*covar[9]))*covarIn[6]/TMath::Sqrt((covarIn[0]*covarIn[9]));
7692 covar[7]=TMath::Sqrt((covar[2]*covar[9]))*covarIn[7]/TMath::Sqrt((covarIn[2]*covarIn[9]));
7693 covar[8]=TMath::Sqrt((covar[5]*covar[9]))*covarIn[8]/TMath::Sqrt((covarIn[5]*covarIn[9]));
7695 covar[10]=TMath::Sqrt((covar[0]*covar[14]))*covarIn[10]/TMath::Sqrt((covarIn[0]*covarIn[14]));
7696 covar[11]=TMath::Sqrt((covar[2]*covar[14]))*covarIn[11]/TMath::Sqrt((covarIn[2]*covarIn[14]));
7697 covar[12]=TMath::Sqrt((covar[5]*covar[14]))*covarIn[12]/TMath::Sqrt((covarIn[5]*covarIn[14]));
7698 covar[13]=TMath::Sqrt((covar[9]*covar[14]))*covarIn[13]/TMath::Sqrt((covarIn[9]*covarIn[14]));
7700 seed->AddCovariance(covar);
7703 //_____________________________________________________________________________
7704 Bool_t AliTPCtrackerMI::IsTPCHVDipEvent(AliESDEvent const *esdEvent) {
7706 // check events affected by TPC HV dip
7708 if(!esdEvent) return kFALSE;
7711 if(!AliTPCcalibDB::Instance()) return kFALSE;
7712 AliTPCcalibDB::Instance()->SetRun(esdEvent->GetRunNumber());
7714 // Get HV TPC chamber sensors and calculate the median
7715 AliDCSSensorArray *voltageArray= AliTPCcalibDB::Instance()->GetVoltageSensors(esdEvent->GetRunNumber());
7716 if(!voltageArray) return kFALSE;
7718 TString sensorName="";
7719 Double_t kTPCHVdip = 2.0; // allow for 2V dip as compared to median from given sensor
7722 for(Int_t sector=0; sector<72; sector++)
7724 Char_t sideName='A';
7725 if ((sector/18)%2==1) sideName='C';
7728 sensorName=Form("TPC_ANODE_I_%c%02d_VMEAS",sideName,sector%18);
7731 sensorName=Form("TPC_ANODE_O_%c%02d_0_VMEAS",sideName,sector%18);
7734 AliDCSSensor* sensor = voltageArray->GetSensor(sensorName.Data());
7735 if(!sensor) continue;
7736 TGraph *graph = sensor->GetGraph();
7737 if(!graph) continue;
7738 Double_t median = TMath::Median(graph->GetN(), graph->GetY());
7739 if(median == 0) continue;
7741 //printf("chamber %d, sensor %s, HV %f, median %f\n", sector, sensorName.Data(), sensor->GetValue(esdEvent->GetTimeStamp()), median);
7743 if(TMath::Abs(sensor->GetValue(esdEvent->GetTimeStamp())-median)>kTPCHVdip) {
7751 //________________________________________
7752 void AliTPCtrackerMI::MarkSeedFree(TObject *sd)
7754 // account that this seed is "deleted"
7755 AliTPCseed* seed = dynamic_cast<AliTPCseed*>(sd);
7757 AliError(Form("Freeing of non-AliTPCseed %p from the pool is requested",sd));
7760 int id = seed->GetPoolID();
7762 AliError(Form("Freeing of seed %p NOT from the pool is requested",sd));
7765 // AliInfo(Form("%d %p",id, seed));
7766 fSeedsPool->RemoveAt(id);
7767 if (fFreeSeedsID.GetSize()<=fNFreeSeeds) fFreeSeedsID.Set( 2*fNFreeSeeds + 100 );
7768 fFreeSeedsID.GetArray()[fNFreeSeeds++] = id;
7771 //________________________________________
7772 TObject *&AliTPCtrackerMI::NextFreeSeed()
7774 // return next free slot where the seed can be created
7775 fLastSeedID = fNFreeSeeds ? fFreeSeedsID.GetArray()[--fNFreeSeeds] : fSeedsPool->GetEntriesFast();
7776 // AliInfo(Form("%d",fLastSeedID));
7777 return (*fSeedsPool)[ fLastSeedID ];
7781 //________________________________________
7782 void AliTPCtrackerMI::ResetSeedsPool()
7784 // mark all seeds in the pool as unused
7785 AliInfo(Form("CurrentSize: %d, BookedUpTo: %d, free: %d",fSeedsPool->GetSize(),fSeedsPool->GetEntriesFast(),fNFreeSeeds));
7787 fSeedsPool->Clear("C"); // RS: nominally the seeds may allocate memory...