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); where nis bigger 0
72 // The debug level - different procedure produce tree for numerical debugging
73 // To enable them set AliTPCReconstructor::SetStreamLevel(n); where nis bigger 1
76 //-------------------------------------------------------
81 #include "Riostream.h"
82 #include <TClonesArray.h>
84 #include <TObjArray.h>
87 #include "AliComplexCluster.h"
88 #include "AliESDEvent.h"
89 #include "AliESDtrack.h"
90 #include "AliESDVertex.h"
94 #include "AliRunLoader.h"
95 #include "AliTPCClustersRow.h"
96 #include "AliTPCParam.h"
97 #include "AliTPCReconstructor.h"
98 #include "AliTPCpolyTrack.h"
99 #include "AliTPCreco.h"
100 #include "AliTPCseed.h"
101 #include "AliTPCtrackerMI.h"
102 #include "TStopwatch.h"
103 #include "AliTPCReconstructor.h"
105 #include "TTreeStream.h"
106 #include "AliAlignObj.h"
107 #include "AliTrackPointArray.h"
109 //#include "AliTPCTransform.h"
113 ClassImp(AliTPCtrackerMI)
116 class AliTPCFastMath {
119 static Double_t FastAsin(Double_t x);
121 static Double_t fgFastAsin[20000]; //lookup table for fast asin computation
124 Double_t AliTPCFastMath::fgFastAsin[20000];
125 AliTPCFastMath gAliTPCFastMath; // needed to fill the LUT
127 AliTPCFastMath::AliTPCFastMath(){
129 // initialized lookup table;
130 for (Int_t i=0;i<10000;i++){
131 fgFastAsin[2*i] = TMath::ASin(i/10000.);
132 fgFastAsin[2*i+1] = (TMath::ASin((i+1)/10000.)-fgFastAsin[2*i]);
136 Double_t AliTPCFastMath::FastAsin(Double_t x){
138 // return asin using lookup table
140 Int_t index = int(x*10000);
141 return fgFastAsin[2*index]+(x*10000.-index)*fgFastAsin[2*index+1];
144 Int_t index = int(x*10000);
145 return -(fgFastAsin[2*index]+(x*10000.-index)*fgFastAsin[2*index+1]);
147 //__________________________________________________________________
148 AliTPCtrackerMI::AliTPCtrackerMI()
170 // default constructor
173 //_____________________________________________________________________
177 Int_t AliTPCtrackerMI::UpdateTrack(AliTPCseed * track, Int_t accept){
179 //update track information using current cluster - track->fCurrentCluster
182 AliTPCclusterMI* c =track->GetCurrentCluster();
183 if (accept>0) track->SetCurrentClusterIndex1(track->GetCurrentClusterIndex1() | 0x8000); //sign not accepted clusters
185 UInt_t i = track->GetCurrentClusterIndex1();
187 Int_t sec=(i&0xff000000)>>24;
188 //Int_t row = (i&0x00ff0000)>>16;
189 track->SetRow((i&0x00ff0000)>>16);
190 track->SetSector(sec);
191 // Int_t index = i&0xFFFF;
192 if (sec>=fParam->GetNInnerSector()) track->SetRow(track->GetRow()+fParam->GetNRowLow());
193 track->SetClusterIndex2(track->GetRow(), i);
194 //track->fFirstPoint = row;
195 //if ( track->fLastPoint<row) track->fLastPoint =row;
196 // if (track->fRow<0 || track->fRow>160) {
197 // printf("problem\n");
199 if (track->GetFirstPoint()>track->GetRow())
200 track->SetFirstPoint(track->GetRow());
201 if (track->GetLastPoint()<track->GetRow())
202 track->SetLastPoint(track->GetRow());
205 track->SetClusterPointer(track->GetRow(),c);
208 Double_t angle2 = track->GetSnp()*track->GetSnp();
210 //SET NEW Track Point
212 if (angle2<1) //PH sometimes angle2 is very big. To be investigated...
214 angle2 = TMath::Sqrt(angle2/(1-angle2));
215 AliTPCTrackerPoint &point =*(track->GetTrackPoint(track->GetRow()));
217 point.SetSigmaY(c->GetSigmaY2()/track->GetCurrentSigmaY2());
218 point.SetSigmaZ(c->GetSigmaZ2()/track->GetCurrentSigmaZ2());
219 point.SetErrY(sqrt(track->GetErrorY2()));
220 point.SetErrZ(sqrt(track->GetErrorZ2()));
222 point.SetX(track->GetX());
223 point.SetY(track->GetY());
224 point.SetZ(track->GetZ());
225 point.SetAngleY(angle2);
226 point.SetAngleZ(track->GetTgl());
227 if (point.IsShared()){
228 track->SetErrorY2(track->GetErrorY2()*4);
229 track->SetErrorZ2(track->GetErrorZ2()*4);
233 Double_t chi2 = track->GetPredictedChi2(track->GetCurrentCluster());
235 track->SetErrorY2(track->GetErrorY2()*1.3);
236 track->SetErrorY2(track->GetErrorY2()+0.01);
237 track->SetErrorZ2(track->GetErrorZ2()*1.3);
238 track->SetErrorZ2(track->GetErrorZ2()+0.005);
240 if (accept>0) return 0;
241 if (track->GetNumberOfClusters()%20==0){
242 // if (track->fHelixIn){
243 // TClonesArray & larr = *(track->fHelixIn);
244 // Int_t ihelix = larr.GetEntriesFast();
245 // new(larr[ihelix]) AliHelix(*track) ;
248 track->SetNoCluster(0);
249 return track->Update(c,chi2,i);
254 Int_t AliTPCtrackerMI::AcceptCluster(AliTPCseed * seed, AliTPCclusterMI * cluster, Float_t factor,
255 Float_t cory, Float_t corz)
258 // decide according desired precision to accept given
259 // cluster for tracking
260 Double_t sy2=ErrY2(seed,cluster)*cory;
261 Double_t sz2=ErrZ2(seed,cluster)*corz;
262 //sy2=ErrY2(seed,cluster)*cory;
263 //sz2=ErrZ2(seed,cluster)*cory;
265 Double_t sdistancey2 = sy2+seed->GetSigmaY2();
266 Double_t sdistancez2 = sz2+seed->GetSigmaZ2();
268 Double_t rdistancey2 = (seed->GetCurrentCluster()->GetY()-seed->GetY())*
269 (seed->GetCurrentCluster()->GetY()-seed->GetY())/sdistancey2;
270 Double_t rdistancez2 = (seed->GetCurrentCluster()->GetZ()-seed->GetZ())*
271 (seed->GetCurrentCluster()->GetZ()-seed->GetZ())/sdistancez2;
273 Double_t rdistance2 = rdistancey2+rdistancez2;
276 if (rdistance2>16) return 3;
279 if ((rdistancey2>9.*factor || rdistancez2>9.*factor) && cluster->GetType()==0)
280 return 2; //suspisiouce - will be changed
282 if ((rdistancey2>6.25*factor || rdistancez2>6.25*factor) && cluster->GetType()>0)
283 // strict cut on overlaped cluster
284 return 2; //suspisiouce - will be changed
286 if ( (rdistancey2>1.*factor || rdistancez2>6.25*factor )
287 && cluster->GetType()<0){
288 seed->SetNFoundable(seed->GetNFoundable()-1);
297 //_____________________________________________________________________________
298 AliTPCtrackerMI::AliTPCtrackerMI(const AliTPCParam *par):
300 fkNIS(par->GetNInnerSector()/2),
302 fkNOS(par->GetNOuterSector()/2),
319 //---------------------------------------------------------------------
320 // The main TPC tracker constructor
321 //---------------------------------------------------------------------
322 fInnerSec=new AliTPCSector[fkNIS];
323 fOuterSec=new AliTPCSector[fkNOS];
326 for (i=0; i<fkNIS; i++) fInnerSec[i].Setup(par,0);
327 for (i=0; i<fkNOS; i++) fOuterSec[i].Setup(par,1);
330 Int_t nrowlow = par->GetNRowLow();
331 Int_t nrowup = par->GetNRowUp();
334 for (Int_t i=0;i<nrowlow;i++){
335 fXRow[i] = par->GetPadRowRadiiLow(i);
336 fPadLength[i]= par->GetPadPitchLength(0,i);
337 fYMax[i] = fXRow[i]*TMath::Tan(0.5*par->GetInnerAngle());
341 for (Int_t i=0;i<nrowup;i++){
342 fXRow[i+nrowlow] = par->GetPadRowRadiiUp(i);
343 fPadLength[i+nrowlow] = par->GetPadPitchLength(60,i);
344 fYMax[i+nrowlow] = fXRow[i+nrowlow]*TMath::Tan(0.5*par->GetOuterAngle());
347 fDebugStreamer = new TTreeSRedirector("TPCdebug.root");
349 //________________________________________________________________________
350 AliTPCtrackerMI::AliTPCtrackerMI(const AliTPCtrackerMI &t):
371 //------------------------------------
372 // dummy copy constructor
373 //------------------------------------------------------------------
376 AliTPCtrackerMI & AliTPCtrackerMI::operator=(const AliTPCtrackerMI& /*r*/){
377 //------------------------------
379 //--------------------------------------------------------------
382 //_____________________________________________________________________________
383 AliTPCtrackerMI::~AliTPCtrackerMI() {
384 //------------------------------------------------------------------
385 // TPC tracker destructor
386 //------------------------------------------------------------------
393 if (fDebugStreamer) delete fDebugStreamer;
396 void AliTPCtrackerMI::SetIO()
400 fInput = AliRunLoader::GetTreeR("TPC", kFALSE,AliConfig::GetDefaultEventFolderName());
402 fOutput = AliRunLoader::GetTreeT("TPC", kTRUE,AliConfig::GetDefaultEventFolderName());
404 AliTPCtrack *iotrack= new AliTPCtrack;
405 fOutput->Branch("tracks","AliTPCtrack",&iotrack,32000,100);
411 void AliTPCtrackerMI::SetIO(TTree * input, TTree * output, AliESDEvent * event)
427 AliTPCtrack *iotrack= new AliTPCtrack;
428 // iotrack->fHelixIn = new TClonesArray("AliHelix");
429 //iotrack->fHelixOut = new TClonesArray("AliHelix");
430 fOutput->Branch("tracks","AliTPCtrack",&iotrack,32000,100);
433 if (output && (fDebug&2)){
434 //write the full seed information if specified in debug mode
436 fSeedTree = new TTree("Seeds","Seeds");
437 AliTPCseed * vseed = new AliTPCseed;
439 TClonesArray * arrtr = new TClonesArray("AliTPCTrackPoint",160);
440 arrtr->ExpandCreateFast(160);
441 TClonesArray * arre = new TClonesArray("AliTPCExactPoint",160);
443 vseed->SetPoints(arrtr);
444 vseed->SetEPoints(arre);
445 // vseed->fClusterPoints = arrcl;
446 fSeedTree->Branch("seeds","AliTPCseed",&vseed,32000,99);
449 fTreeDebug = new TTree("trackDebug","trackDebug");
450 TClonesArray * arrd = new TClonesArray("AliTPCTrackPoint2",0);
451 fTreeDebug->Branch("debug",&arrd,32000,99);
459 void AliTPCtrackerMI::FillESD(TObjArray* arr)
463 //fill esds using updated tracks
465 // write tracks to the event
466 // store index of the track
467 Int_t nseed=arr->GetEntriesFast();
468 //FindKinks(arr,fEvent);
469 for (Int_t i=0; i<nseed; i++) {
470 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
473 // pt->PropagateTo(fParam->GetInnerRadiusLow());
474 if (pt->GetKinkIndex(0)<=0){ //don't propagate daughter tracks
475 pt->PropagateTo(fParam->GetInnerRadiusLow());
478 if (( pt->GetPoints()[2]- pt->GetPoints()[0])>5 && pt->GetPoints()[3]>0.8){
480 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
481 iotrack.SetTPCPoints(pt->GetPoints());
482 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
483 iotrack.SetV0Indexes(pt->GetV0Indexes());
484 // iotrack.SetTPCpid(pt->fTPCr);
485 //iotrack.SetTPCindex(i);
486 fEvent->AddTrack(&iotrack);
490 if ( (pt->GetNumberOfClusters()>70)&& (Float_t(pt->GetNumberOfClusters())/Float_t(pt->GetNFoundable()))>0.55) {
492 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
493 iotrack.SetTPCPoints(pt->GetPoints());
494 //iotrack.SetTPCindex(i);
495 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
496 iotrack.SetV0Indexes(pt->GetV0Indexes());
497 // iotrack.SetTPCpid(pt->fTPCr);
498 fEvent->AddTrack(&iotrack);
502 // short tracks - maybe decays
504 if ( (pt->GetNumberOfClusters()>30) && (Float_t(pt->GetNumberOfClusters())/Float_t(pt->GetNFoundable()))>0.70) {
505 Int_t found,foundable,shared;
506 pt->GetClusterStatistic(0,60,found, foundable,shared,kFALSE);
507 if ( (found>20) && (pt->GetNShared()/float(pt->GetNumberOfClusters())<0.2)){
509 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
510 //iotrack.SetTPCindex(i);
511 iotrack.SetTPCPoints(pt->GetPoints());
512 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
513 iotrack.SetV0Indexes(pt->GetV0Indexes());
514 //iotrack.SetTPCpid(pt->fTPCr);
515 fEvent->AddTrack(&iotrack);
520 if ( (pt->GetNumberOfClusters()>20) && (Float_t(pt->GetNumberOfClusters())/Float_t(pt->GetNFoundable()))>0.8) {
521 Int_t found,foundable,shared;
522 pt->GetClusterStatistic(0,60,found, foundable,shared,kFALSE);
523 if (found<20) continue;
524 if (pt->GetNShared()/float(pt->GetNumberOfClusters())>0.2) continue;
527 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
528 iotrack.SetTPCPoints(pt->GetPoints());
529 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
530 iotrack.SetV0Indexes(pt->GetV0Indexes());
531 //iotrack.SetTPCpid(pt->fTPCr);
532 //iotrack.SetTPCindex(i);
533 fEvent->AddTrack(&iotrack);
536 // short tracks - secondaties
538 if ( (pt->GetNumberOfClusters()>30) ) {
539 Int_t found,foundable,shared;
540 pt->GetClusterStatistic(128,158,found, foundable,shared,kFALSE);
541 if ( (found>20) && (pt->GetNShared()/float(pt->GetNumberOfClusters())<0.2) &&float(found)/float(foundable)>0.8){
543 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
544 iotrack.SetTPCPoints(pt->GetPoints());
545 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
546 iotrack.SetV0Indexes(pt->GetV0Indexes());
547 //iotrack.SetTPCpid(pt->fTPCr);
548 //iotrack.SetTPCindex(i);
549 fEvent->AddTrack(&iotrack);
554 if ( (pt->GetNumberOfClusters()>15)) {
555 Int_t found,foundable,shared;
556 pt->GetClusterStatistic(138,158,found, foundable,shared,kFALSE);
557 if (found<15) continue;
558 if (foundable<=0) continue;
559 if (pt->GetNShared()/float(pt->GetNumberOfClusters())>0.2) continue;
560 if (float(found)/float(foundable)<0.8) continue;
563 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
564 iotrack.SetTPCPoints(pt->GetPoints());
565 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
566 iotrack.SetV0Indexes(pt->GetV0Indexes());
567 // iotrack.SetTPCpid(pt->fTPCr);
568 //iotrack.SetTPCindex(i);
569 fEvent->AddTrack(&iotrack);
574 printf("Number of filled ESDs-\t%d\n",fEvent->GetNumberOfTracks());
577 void AliTPCtrackerMI::WriteTracks(TTree * tree)
580 // write tracks from seed array to selected tree
584 AliTPCtrack *iotrack= new AliTPCtrack;
585 fOutput->Branch("tracks","AliTPCtrack",&iotrack,32000,100);
590 void AliTPCtrackerMI::WriteTracks()
593 // write tracks to the given output tree -
594 // output specified with SetIO routine
601 AliTPCtrack *iotrack= 0;
602 Int_t nseed=fSeeds->GetEntriesFast();
603 //for (Int_t i=0; i<nseed; i++) {
604 // iotrack= (AliTPCtrack*)fSeeds->UncheckedAt(i);
605 // if (iotrack) break;
607 //TBranch * br = fOutput->Branch("tracks","AliTPCtrack",&iotrack,32000,100);
608 TBranch * br = fOutput->GetBranch("tracks");
609 br->SetAddress(&iotrack);
611 for (Int_t i=0; i<nseed; i++) {
612 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i);
614 AliTPCtrack * track = new AliTPCtrack(*pt);
617 // br->SetAddress(&iotrack);
622 //fOutput->GetDirectory()->cd();
628 //write the full seed information if specified in debug mode
630 AliTPCseed * vseed = new AliTPCseed;
632 TClonesArray * arrtr = new TClonesArray("AliTPCTrackPoint",160);
633 arrtr->ExpandCreateFast(160);
634 //TClonesArray * arrcl = new TClonesArray("AliTPCclusterMI",160);
635 //arrcl->ExpandCreateFast(160);
636 TClonesArray * arre = new TClonesArray("AliTPCExactPoint",160);
638 vseed->SetPoints(arrtr);
639 vseed->SetEPoints(arre);
640 // vseed->fClusterPoints = arrcl;
641 //TBranch * brseed = seedtree->Branch("seeds","AliTPCseed",&vseed,32000,99);
642 TBranch * brseed = fSeedTree->GetBranch("seeds");
644 Int_t nseed=fSeeds->GetEntriesFast();
646 for (Int_t i=0; i<nseed; i++) {
647 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i);
649 pt->SetPoints(arrtr);
650 // pt->fClusterPoints = arrcl;
651 pt->SetEPoints(arre);
654 brseed->SetAddress(&vseed);
658 // pt->fClusterPoints = 0;
661 if (fTreeDebug) fTreeDebug->Write();
669 Double_t AliTPCtrackerMI::ErrY2(AliTPCseed* seed, AliTPCclusterMI * cl){
672 //seed->SetErrorY2(0.1);
674 //calculate look-up table at the beginning
675 static Bool_t ginit = kFALSE;
676 static Float_t gnoise1,gnoise2,gnoise3;
677 static Float_t ggg1[10000];
678 static Float_t ggg2[10000];
679 static Float_t ggg3[10000];
680 static Float_t glandau1[10000];
681 static Float_t glandau2[10000];
682 static Float_t glandau3[10000];
684 static Float_t gcor01[500];
685 static Float_t gcor02[500];
686 static Float_t gcorp[500];
691 for (Int_t i=1;i<500;i++){
692 Float_t rsigma = float(i)/100.;
693 gcor02[i] = TMath::Max(0.78 +TMath::Exp(7.4*(rsigma-1.2)),0.6);
694 gcor01[i] = TMath::Max(0.72 +TMath::Exp(3.36*(rsigma-1.2)),0.6);
695 gcorp[i] = TMath::Max(TMath::Power((rsigma+0.5),1.5),1.2);
699 for (Int_t i=3;i<10000;i++){
703 Float_t amp = float(i);
704 Float_t padlength =0.75;
705 gnoise1 = 0.0004/padlength;
706 Float_t nel = 0.268*amp;
707 Float_t nprim = 0.155*amp;
708 ggg1[i] = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.001*nel/(padlength*padlength))/nel;
709 glandau1[i] = (2.+0.12*nprim)*0.5* (2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
710 if (glandau1[i]>1) glandau1[i]=1;
711 glandau1[i]*=padlength*padlength/12.;
715 gnoise2 = 0.0004/padlength;
718 ggg2[i] = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
719 glandau2[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
720 if (glandau2[i]>1) glandau2[i]=1;
721 glandau2[i]*=padlength*padlength/12.;
726 gnoise3 = 0.0004/padlength;
729 ggg3[i] = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
730 glandau3[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
731 if (glandau3[i]>1) glandau3[i]=1;
732 glandau3[i]*=padlength*padlength/12.;
740 Int_t amp = int(TMath::Abs(cl->GetQ()));
742 seed->SetErrorY2(1.);
746 Float_t z = TMath::Abs(fParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
747 Int_t ctype = cl->GetType();
748 Float_t padlength= GetPadPitchLength(seed->GetRow());
749 Double_t angle2 = seed->GetSnp()*seed->GetSnp();
750 angle2 = angle2/(1-angle2);
753 Int_t rsigmay = int(100.*cl->GetSigmaY2()/(seed->GetCurrentSigmaY2()));
756 if (fSectors==fInnerSec){
758 res = ggg1[amp]*z+glandau1[amp]*angle2;
759 if (ctype==0) res *= gcor01[rsigmay];
762 res*= gcorp[rsigmay];
768 res = ggg2[amp]*z+glandau2[amp]*angle2;
769 if (ctype==0) res *= gcor02[rsigmay];
772 res*= gcorp[rsigmay];
777 res = ggg3[amp]*z+glandau3[amp]*angle2;
778 if (ctype==0) res *= gcor02[rsigmay];
781 res*= gcorp[rsigmay];
788 res*=2.4; // overestimate error 2 times
795 seed->SetErrorY2(res);
803 Double_t AliTPCtrackerMI::ErrZ2(AliTPCseed* seed, AliTPCclusterMI * cl){
806 //seed->SetErrorY2(0.1);
808 //calculate look-up table at the beginning
809 static Bool_t ginit = kFALSE;
810 static Float_t gnoise1,gnoise2,gnoise3;
811 static Float_t ggg1[10000];
812 static Float_t ggg2[10000];
813 static Float_t ggg3[10000];
814 static Float_t glandau1[10000];
815 static Float_t glandau2[10000];
816 static Float_t glandau3[10000];
818 static Float_t gcor01[1000];
819 static Float_t gcor02[1000];
820 static Float_t gcorp[1000];
825 for (Int_t i=1;i<1000;i++){
826 Float_t rsigma = float(i)/100.;
827 gcor02[i] = TMath::Max(0.81 +TMath::Exp(6.8*(rsigma-1.2)),0.6);
828 gcor01[i] = TMath::Max(0.72 +TMath::Exp(2.04*(rsigma-1.2)),0.6);
829 gcorp[i] = TMath::Max(TMath::Power((rsigma+0.5),1.5),1.2);
833 for (Int_t i=3;i<10000;i++){
837 Float_t amp = float(i);
838 Float_t padlength =0.75;
839 gnoise1 = 0.0004/padlength;
840 Float_t nel = 0.268*amp;
841 Float_t nprim = 0.155*amp;
842 ggg1[i] = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.001*nel/(padlength*padlength))/nel;
843 glandau1[i] = (2.+0.12*nprim)*0.5* (2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
844 if (glandau1[i]>1) glandau1[i]=1;
845 glandau1[i]*=padlength*padlength/12.;
849 gnoise2 = 0.0004/padlength;
852 ggg2[i] = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
853 glandau2[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
854 if (glandau2[i]>1) glandau2[i]=1;
855 glandau2[i]*=padlength*padlength/12.;
860 gnoise3 = 0.0004/padlength;
863 ggg3[i] = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
864 glandau3[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
865 if (glandau3[i]>1) glandau3[i]=1;
866 glandau3[i]*=padlength*padlength/12.;
874 Int_t amp = int(TMath::Abs(cl->GetQ()));
876 seed->SetErrorY2(1.);
880 Float_t z = TMath::Abs(fParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
881 Int_t ctype = cl->GetType();
882 Float_t padlength= GetPadPitchLength(seed->GetRow());
884 Double_t angle2 = seed->GetSnp()*seed->GetSnp();
885 // if (angle2<0.6) angle2 = 0.6;
886 angle2 = seed->GetTgl()*seed->GetTgl()*(1+angle2/(1-angle2));
889 Int_t rsigmaz = int(100.*cl->GetSigmaZ2()/(seed->GetCurrentSigmaZ2()));
892 if (fSectors==fInnerSec){
894 res = ggg1[amp]*z+glandau1[amp]*angle2;
895 if (ctype==0) res *= gcor01[rsigmaz];
898 res*= gcorp[rsigmaz];
904 res = ggg2[amp]*z+glandau2[amp]*angle2;
905 if (ctype==0) res *= gcor02[rsigmaz];
908 res*= gcorp[rsigmaz];
913 res = ggg3[amp]*z+glandau3[amp]*angle2;
914 if (ctype==0) res *= gcor02[rsigmaz];
917 res*= gcorp[rsigmaz];
926 if ((ctype<0) &&<70){
934 seed->SetErrorZ2(res);
941 Double_t AliTPCtrackerMI::ErrZ2(AliTPCseed* seed, AliTPCclusterMI * cl){
944 //seed->SetErrorZ2(0.1);
948 Float_t z = TMath::Abs(fParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
950 Float_t rsigmaz = cl->GetSigmaZ2()/(seed->fCurrentSigmaZ2);
951 Int_t ctype = cl->GetType();
952 Float_t amp = TMath::Abs(cl->GetQ());
957 Float_t landau=2 ; //landau fluctuation part
958 Float_t gg=2; // gg fluctuation part
959 Float_t padlength= GetPadPitchLength(seed->GetX());
961 if (fSectors==fInnerSec){
962 snoise2 = 0.0004/padlength;
965 gg = (2+0.001*nel/(padlength*padlength))/nel;
966 landau = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
967 if (landau>1) landau=1;
970 snoise2 = 0.0004/padlength;
973 gg = (2+0.0008*nel/(padlength*padlength))/nel;
974 landau = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
975 if (landau>1) landau=1;
977 Float_t sdiff = gg*fParam->GetDiffT()*fParam->GetDiffT()*z;
980 Float_t angle2 = seed->GetSnp()*seed->GetSnp();
981 angle2 = TMath::Sqrt((1-angle2));
982 if (angle2<0.6) angle2 = 0.6;
985 Float_t angle = seed->GetTgl()/angle2;
986 Float_t angular = landau*angle*angle*padlength*padlength/12.;
987 Float_t res = sdiff + angular;
990 if ((ctype==0) && (fSectors ==fOuterSec))
991 res *= 0.81 +TMath::Exp(6.8*(rsigmaz-1.2));
993 if ((ctype==0) && (fSectors ==fInnerSec))
994 res *= 0.72 +TMath::Exp(2.04*(rsigmaz-1.2));
998 res*= TMath::Power(rsigmaz+0.5,1.5); //0.31+0.147*ctype;
1004 if ((ctype<0) &&<70){
1012 seed->SetErrorZ2(res);
1020 void AliTPCtrackerMI::RotateToLocal(AliTPCseed *seed)
1022 //rotate to track "local coordinata
1023 Float_t x = seed->GetX();
1024 Float_t y = seed->GetY();
1025 Float_t ymax = x*TMath::Tan(0.5*fSectors->GetAlpha());
1028 seed->SetRelativeSector((seed->GetRelativeSector()+1) % fN);
1029 if (!seed->Rotate(fSectors->GetAlpha()))
1031 } else if (y <-ymax) {
1032 seed->SetRelativeSector((seed->GetRelativeSector()-1+fN) % fN);
1033 if (!seed->Rotate(-fSectors->GetAlpha()))
1041 //_____________________________________________________________________________
1042 Double_t AliTPCtrackerMI::F1old(Double_t x1,Double_t y1,
1043 Double_t x2,Double_t y2,
1044 Double_t x3,Double_t y3)
1046 //-----------------------------------------------------------------
1047 // Initial approximation of the track curvature
1048 //-----------------------------------------------------------------
1049 Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
1050 Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
1051 (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
1052 Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
1053 (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
1055 Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
1056 if ( xr*xr+yr*yr<=0.00000000000001) return 100;
1057 return -xr*yr/sqrt(xr*xr+yr*yr);
1062 //_____________________________________________________________________________
1063 Double_t AliTPCtrackerMI::F1(Double_t x1,Double_t y1,
1064 Double_t x2,Double_t y2,
1065 Double_t x3,Double_t y3)
1067 //-----------------------------------------------------------------
1068 // Initial approximation of the track curvature
1069 //-----------------------------------------------------------------
1075 Double_t det = x3*y2-x2*y3;
1080 Double_t u = 0.5* (x2*(x2-x3)+y2*(y2-y3))/det;
1081 Double_t x0 = x3*0.5-y3*u;
1082 Double_t y0 = y3*0.5+x3*u;
1083 Double_t c2 = 1/TMath::Sqrt(x0*x0+y0*y0);
1089 Double_t AliTPCtrackerMI::F2(Double_t x1,Double_t y1,
1090 Double_t x2,Double_t y2,
1091 Double_t x3,Double_t y3)
1093 //-----------------------------------------------------------------
1094 // Initial approximation of the track curvature
1095 //-----------------------------------------------------------------
1101 Double_t det = x3*y2-x2*y3;
1106 Double_t u = 0.5* (x2*(x2-x3)+y2*(y2-y3))/det;
1107 Double_t x0 = x3*0.5-y3*u;
1108 Double_t y0 = y3*0.5+x3*u;
1109 Double_t c2 = 1/TMath::Sqrt(x0*x0+y0*y0);
1118 //_____________________________________________________________________________
1119 Double_t AliTPCtrackerMI::F2old(Double_t x1,Double_t y1,
1120 Double_t x2,Double_t y2,
1121 Double_t x3,Double_t y3)
1123 //-----------------------------------------------------------------
1124 // Initial approximation of the track curvature times center of curvature
1125 //-----------------------------------------------------------------
1126 Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
1127 Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
1128 (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
1129 Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
1130 (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
1132 Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
1134 return -a/(d*y1-b)*xr/sqrt(xr*xr+yr*yr);
1137 //_____________________________________________________________________________
1138 Double_t AliTPCtrackerMI::F3(Double_t x1,Double_t y1,
1139 Double_t x2,Double_t y2,
1140 Double_t z1,Double_t z2)
1142 //-----------------------------------------------------------------
1143 // Initial approximation of the tangent of the track dip angle
1144 //-----------------------------------------------------------------
1145 return (z1 - z2)/sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
1149 Double_t AliTPCtrackerMI::F3n(Double_t x1,Double_t y1,
1150 Double_t x2,Double_t y2,
1151 Double_t z1,Double_t z2, Double_t c)
1153 //-----------------------------------------------------------------
1154 // Initial approximation of the tangent of the track dip angle
1155 //-----------------------------------------------------------------
1159 //angle1 = (z1-z2)*c/(TMath::ASin(c*x1-ni)-TMath::ASin(c*x2-ni));
1161 Double_t d = TMath::Sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
1162 if (TMath::Abs(d*c*0.5)>1) return 0;
1163 // Double_t angle2 = TMath::ASin(d*c*0.5);
1164 // Double_t angle2 = AliTPCFastMath::FastAsin(d*c*0.5);
1165 Double_t angle2 = (d*c*0.5>0.1)? TMath::ASin(d*c*0.5): AliTPCFastMath::FastAsin(d*c*0.5);
1167 angle2 = (z1-z2)*c/(angle2*2.);
1171 Bool_t AliTPCtrackerMI::GetProlongation(Double_t x1, Double_t x2, Double_t x[5], Double_t &y, Double_t &z)
1172 {//-----------------------------------------------------------------
1173 // This function find proloncation of a track to a reference plane x=x2.
1174 //-----------------------------------------------------------------
1178 if (TMath::Abs(x[4]*x1 - x[2]) >= 0.999) {
1182 Double_t c1=x[4]*x1 - x[2], r1=sqrt(1.- c1*c1);
1183 Double_t c2=x[4]*x2 - x[2], r2=sqrt(1.- c2*c2);
1187 Double_t dy = dx*(c1+c2)/(r1+r2);
1190 Double_t delta = x[4]*dx*(c1+c2)/(c1*r2 + c2*r1);
1192 if (TMath::Abs(delta)>0.01){
1193 dz = x[3]*TMath::ASin(delta)/x[4];
1195 dz = x[3]*AliTPCFastMath::FastAsin(delta)/x[4];
1198 //dz = x[3]*AliTPCFastMath::FastAsin(delta)/x[4];
1206 Int_t AliTPCtrackerMI::LoadClusters (TTree *tree)
1211 return LoadClusters();
1214 Int_t AliTPCtrackerMI::LoadClusters()
1217 // load clusters to the memory
1218 AliTPCClustersRow *clrow= new AliTPCClustersRow;
1219 clrow->SetClass("AliTPCclusterMI");
1221 clrow->GetArray()->ExpandCreateFast(10000);
1223 // TTree * tree = fClustersArray.GetTree();
1225 TTree * tree = fInput;
1226 TBranch * br = tree->GetBranch("Segment");
1227 br->SetAddress(&clrow);
1229 Int_t j=Int_t(tree->GetEntries());
1230 for (Int_t i=0; i<j; i++) {
1234 fParam->AdjustSectorRow(clrow->GetID(),sec,row);
1235 for (Int_t icl=0; icl<clrow->GetArray()->GetEntriesFast(); icl++){
1236 Transform((AliTPCclusterMI*)(clrow->GetArray()->At(icl)));
1239 AliTPCRow * tpcrow=0;
1242 tpcrow = &(fInnerSec[sec%fkNIS][row]);
1246 tpcrow = &(fOuterSec[(sec-fkNIS*2)%fkNOS][row]);
1247 left = (sec-fkNIS*2)/fkNOS;
1250 tpcrow->SetN1(clrow->GetArray()->GetEntriesFast());
1251 tpcrow->SetClusters1(new AliTPCclusterMI[tpcrow->GetN1()]);
1252 for (Int_t i=0;i<tpcrow->GetN1();i++)
1253 tpcrow->SetCluster1(i, *(AliTPCclusterMI*)(clrow->GetArray()->At(i)));
1256 tpcrow->SetN2(clrow->GetArray()->GetEntriesFast());
1257 tpcrow->SetClusters2(new AliTPCclusterMI[tpcrow->GetN2()]);
1258 for (Int_t i=0;i<tpcrow->GetN2();i++)
1259 tpcrow->SetCluster2(i,*(AliTPCclusterMI*)(clrow->GetArray()->At(i)));
1270 void AliTPCtrackerMI::UnloadClusters()
1273 // unload clusters from the memory
1275 Int_t nrows = fOuterSec->GetNRows();
1276 for (Int_t sec = 0;sec<fkNOS;sec++)
1277 for (Int_t row = 0;row<nrows;row++){
1278 AliTPCRow* tpcrow = &(fOuterSec[sec%fkNOS][row]);
1280 // if (tpcrow->fClusters1) delete []tpcrow->fClusters1;
1281 // if (tpcrow->fClusters2) delete []tpcrow->fClusters2;
1283 tpcrow->ResetClusters();
1286 nrows = fInnerSec->GetNRows();
1287 for (Int_t sec = 0;sec<fkNIS;sec++)
1288 for (Int_t row = 0;row<nrows;row++){
1289 AliTPCRow* tpcrow = &(fInnerSec[sec%fkNIS][row]);
1291 // if (tpcrow->fClusters1) delete []tpcrow->fClusters1;
1292 //if (tpcrow->fClusters2) delete []tpcrow->fClusters2;
1294 tpcrow->ResetClusters();
1300 void AliTPCtrackerMI::Transform(AliTPCclusterMI * cluster){
1304 // AliTPCTransform trafo;
1306 // Double_t x[3]={cluster->GetX(),cluster->GetY(),cluster->GetZ()};
1307 // Int_t i[1]={cluster->GetDetector()};
1309 // trafo.Transform(x,i,0,1);
1311 // cluster->SetX(x[0]);
1312 // cluster->SetY(x[1]);
1313 // cluster->SetZ(x[2]);
1322 //if (!fParam->IsGeoRead()) fParam->ReadGeoMatrices();
1323 TGeoHMatrix *mat = fParam->GetClusterMatrix(cluster->GetDetector());
1325 Double_t pos[3]= {cluster->GetX(),cluster->GetY(),cluster->GetZ()};
1326 Double_t posC[3]={cluster->GetX(),cluster->GetY(),cluster->GetZ()};
1327 if (mat) mat->LocalToMaster(pos,posC);
1329 // chack Loading of Geo matrices from GeoManager - TEMPORARY FIX
1331 //cluster->SetX(posC[0]);
1332 //cluster->SetY(posC[1]);
1333 //cluster->SetZ(posC[2]);
1336 //_____________________________________________________________________________
1337 Int_t AliTPCtrackerMI::LoadOuterSectors() {
1338 //-----------------------------------------------------------------
1339 // This function fills outer TPC sectors with clusters.
1340 //-----------------------------------------------------------------
1341 Int_t nrows = fOuterSec->GetNRows();
1343 for (Int_t sec = 0;sec<fkNOS;sec++)
1344 for (Int_t row = 0;row<nrows;row++){
1345 AliTPCRow* tpcrow = &(fOuterSec[sec%fkNOS][row]);
1346 Int_t sec2 = sec+2*fkNIS;
1348 Int_t ncl = tpcrow->GetN1();
1350 AliTPCclusterMI *c= (tpcrow->GetCluster1(ncl));
1351 index=(((sec2<<8)+row)<<16)+ncl;
1352 tpcrow->InsertCluster(c,index);
1355 ncl = tpcrow->GetN2();
1357 AliTPCclusterMI *c= (tpcrow->GetCluster2(ncl));
1358 index=((((sec2+fkNOS)<<8)+row)<<16)+ncl;
1359 tpcrow->InsertCluster(c,index);
1362 // write indexes for fast acces
1364 for (Int_t i=0;i<510;i++)
1365 tpcrow->SetFastCluster(i,-1);
1366 for (Int_t i=0;i<tpcrow->GetN();i++){
1367 Int_t zi = Int_t((*tpcrow)[i]->GetZ()+255.);
1368 tpcrow->SetFastCluster(zi,i); // write index
1371 for (Int_t i=0;i<510;i++){
1372 if (tpcrow->GetFastCluster(i)<0)
1373 tpcrow->SetFastCluster(i,last);
1375 last = tpcrow->GetFastCluster(i);
1384 //_____________________________________________________________________________
1385 Int_t AliTPCtrackerMI::LoadInnerSectors() {
1386 //-----------------------------------------------------------------
1387 // This function fills inner TPC sectors with clusters.
1388 //-----------------------------------------------------------------
1389 Int_t nrows = fInnerSec->GetNRows();
1391 for (Int_t sec = 0;sec<fkNIS;sec++)
1392 for (Int_t row = 0;row<nrows;row++){
1393 AliTPCRow* tpcrow = &(fInnerSec[sec%fkNIS][row]);
1396 Int_t ncl = tpcrow->GetN1();
1398 AliTPCclusterMI *c= (tpcrow->GetCluster1(ncl));
1399 index=(((sec<<8)+row)<<16)+ncl;
1400 tpcrow->InsertCluster(c,index);
1403 ncl = tpcrow->GetN2();
1405 AliTPCclusterMI *c= (tpcrow->GetCluster2(ncl));
1406 index=((((sec+fkNIS)<<8)+row)<<16)+ncl;
1407 tpcrow->InsertCluster(c,index);
1410 // write indexes for fast acces
1412 for (Int_t i=0;i<510;i++)
1413 tpcrow->SetFastCluster(i,-1);
1414 for (Int_t i=0;i<tpcrow->GetN();i++){
1415 Int_t zi = Int_t((*tpcrow)[i]->GetZ()+255.);
1416 tpcrow->SetFastCluster(zi,i); // write index
1419 for (Int_t i=0;i<510;i++){
1420 if (tpcrow->GetFastCluster(i)<0)
1421 tpcrow->SetFastCluster(i,last);
1423 last = tpcrow->GetFastCluster(i);
1435 //_________________________________________________________________________
1436 AliTPCclusterMI *AliTPCtrackerMI::GetClusterMI(Int_t index) const {
1437 //--------------------------------------------------------------------
1438 // Return pointer to a given cluster
1439 //--------------------------------------------------------------------
1440 if (index<0) return 0; // no cluster
1441 Int_t sec=(index&0xff000000)>>24;
1442 Int_t row=(index&0x00ff0000)>>16;
1443 Int_t ncl=(index&0x00007fff)>>00;
1445 const AliTPCRow * tpcrow=0;
1446 AliTPCclusterMI * clrow =0;
1448 if (sec<0 || sec>=fkNIS*4) {
1449 AliWarning(Form("Wrong sector %d",sec));
1454 tpcrow = &(fInnerSec[sec%fkNIS][row]);
1455 if (tpcrow==0) return 0;
1458 if (tpcrow->GetN1()<=ncl) return 0;
1459 clrow = tpcrow->GetClusters1();
1462 if (tpcrow->GetN2()<=ncl) return 0;
1463 clrow = tpcrow->GetClusters2();
1467 tpcrow = &(fOuterSec[(sec-fkNIS*2)%fkNOS][row]);
1468 if (tpcrow==0) return 0;
1470 if (sec-2*fkNIS<fkNOS) {
1471 if (tpcrow->GetN1()<=ncl) return 0;
1472 clrow = tpcrow->GetClusters1();
1475 if (tpcrow->GetN2()<=ncl) return 0;
1476 clrow = tpcrow->GetClusters2();
1480 return &(clrow[ncl]);
1486 Int_t AliTPCtrackerMI::FollowToNext(AliTPCseed& t, Int_t nr) {
1487 //-----------------------------------------------------------------
1488 // This function tries to find a track prolongation to next pad row
1489 //-----------------------------------------------------------------
1491 Double_t x= GetXrow(nr), ymax=GetMaxY(nr);
1492 AliTPCclusterMI *cl=0;
1493 Int_t tpcindex= t.GetClusterIndex2(nr);
1495 // update current shape info every 5 pad-row
1496 // if ( (nr%5==0) || t.GetNumberOfClusters()<2 || (t.fCurrentSigmaY2<0.0001) ){
1500 if (fIteration>0 && tpcindex>=-1){ //if we have already clusters
1502 if (tpcindex==-1) return 0; //track in dead zone
1504 cl = t.GetClusterPointer(nr);
1505 if ( (cl==0) ) cl = GetClusterMI(tpcindex);
1506 t.SetCurrentClusterIndex1(tpcindex);
1509 Int_t relativesector = ((tpcindex&0xff000000)>>24)%18; // if previously accepted cluster in different sector
1510 Float_t angle = relativesector*fSectors->GetAlpha()+fSectors->GetAlphaShift();
1512 if (angle<-TMath::Pi()) angle += 2*TMath::Pi();
1513 if (angle>=TMath::Pi()) angle -= 2*TMath::Pi();
1515 if (TMath::Abs(angle-t.GetAlpha())>0.001){
1516 Double_t rotation = angle-t.GetAlpha();
1517 t.SetRelativeSector(relativesector);
1518 if (!t.Rotate(rotation)) return 0;
1520 if (!t.PropagateTo(x)) return 0;
1522 t.SetCurrentCluster(cl);
1524 Int_t accept = AcceptCluster(&t,t.GetCurrentCluster(),1.);
1525 if ((tpcindex&0x8000)==0) accept =0;
1527 //if founded cluster is acceptible
1528 if (cl->IsUsed(11)) { // id cluster is shared inrease uncertainty
1529 t.SetErrorY2(t.GetErrorY2()+0.03);
1530 t.SetErrorZ2(t.GetErrorZ2()+0.03);
1531 t.SetErrorY2(t.GetErrorY2()*3);
1532 t.SetErrorZ2(t.GetErrorZ2()*3);
1534 t.SetNFoundable(t.GetNFoundable()+1);
1535 UpdateTrack(&t,accept);
1540 if (TMath::Abs(t.GetSnp())>AliTPCReconstructor::GetMaxSnpTracker()) return 0; // cut on angle
1542 // not look for new cluster during refitting
1543 t.SetNFoundable(t.GetNFoundable()+1);
1548 // if (TMath::Abs(t.GetSnp())>0.95 || TMath::Abs(x*t.GetC()-t.GetEta())>0.95) return 0;// patch 28 fev 06
1549 Double_t y=t.GetYat(x);
1550 if (TMath::Abs(y)>ymax){
1552 t.SetRelativeSector((t.GetRelativeSector()+1) % fN);
1553 if (!t.Rotate(fSectors->GetAlpha()))
1555 } else if (y <-ymax) {
1556 t.SetRelativeSector((t.GetRelativeSector()-1+fN) % fN);
1557 if (!t.Rotate(-fSectors->GetAlpha()))
1563 if (!t.PropagateTo(x)) {
1564 if (fIteration==0) t.SetRemoval(10);
1568 Double_t z=t.GetZ();
1571 if (!IsActive(t.GetRelativeSector(),nr)) {
1573 t.SetClusterIndex2(nr,-1);
1576 //AliInfo(Form("A - Sector%d phi %f - alpha %f", t.fRelativeSector,y/x, t.GetAlpha()));
1577 Bool_t isActive = IsActive(t.GetRelativeSector(),nr);
1578 Bool_t isActive2 = (nr>=fInnerSec->GetNRows()) ? fOuterSec[t.GetRelativeSector()][nr-fInnerSec->GetNRows()].GetN()>0:fInnerSec[t.GetRelativeSector()][nr].GetN()>0;
1580 if (!isActive || !isActive2) return 0;
1582 const AliTPCRow &krow=GetRow(t.GetRelativeSector(),nr);
1583 if ( (t.GetSigmaY2()<0) || t.GetSigmaZ2()<0) return 0;
1585 Double_t roadz = 1.;
1587 if (TMath::Abs(TMath::Abs(y)-ymax)<krow.GetDeadZone()){
1589 t.SetClusterIndex2(nr,-1);
1594 if (TMath::Abs(z)<(AliTPCReconstructor::GetCtgRange()*x+10) && TMath::Abs(z)<fParam->GetZLength(0) && (TMath::Abs(t.GetSnp())<AliTPCReconstructor::GetMaxSnpTracker()))
1595 t.SetNFoundable(t.GetNFoundable()+1);
1601 // cl = krow.FindNearest2(y+10.,z,roady,roadz,index);
1602 cl = krow.FindNearest2(y,z,roady,roadz,index);
1603 if (cl) t.SetCurrentClusterIndex1(krow.GetIndex(index));
1606 t.SetCurrentCluster(cl);
1608 if (fIteration==2&&cl->IsUsed(10)) return 0;
1609 Int_t accept = AcceptCluster(&t,t.GetCurrentCluster(),1.);
1610 if (fIteration==2&&cl->IsUsed(11)) {
1611 t.SetErrorY2(t.GetErrorY2()+0.03);
1612 t.SetErrorZ2(t.GetErrorZ2()+0.03);
1613 t.SetErrorY2(t.GetErrorY2()*3);
1614 t.SetErrorZ2(t.GetErrorZ2()*3);
1617 if (t.fCurrentCluster->IsUsed(10)){
1622 if (t.fNShared>0.7*t.GetNumberOfClusters()) {
1628 if (accept<3) UpdateTrack(&t,accept);
1631 if ( fIteration==0 && t.GetNFoundable()*0.5 > t.GetNumberOfClusters()) t.SetRemoval(10);
1637 Int_t AliTPCtrackerMI::FollowToNextFast(AliTPCseed& t, Int_t nr) {
1638 //-----------------------------------------------------------------
1639 // This function tries to find a track prolongation to next pad row
1640 //-----------------------------------------------------------------
1642 Double_t x= GetXrow(nr), ymax=GetMaxY(nr);
1644 if (!t.GetProlongation(x,y,z)) {
1650 if (TMath::Abs(y)>ymax){
1653 t.SetRelativeSector((t.GetRelativeSector()+1) % fN);
1654 if (!t.Rotate(fSectors->GetAlpha()))
1656 } else if (y <-ymax) {
1657 t.SetRelativeSector((t.GetRelativeSector()-1+fN) % fN);
1658 if (!t.Rotate(-fSectors->GetAlpha()))
1661 if (!t.PropagateTo(x)) {
1664 t.GetProlongation(x,y,z);
1667 // update current shape info every 3 pad-row
1668 if ( (nr%6==0) || t.GetNumberOfClusters()<2 || (t.GetCurrentSigmaY2()<0.0001) ){
1669 // t.fCurrentSigmaY = GetSigmaY(&t);
1670 //t.fCurrentSigmaZ = GetSigmaZ(&t);
1674 AliTPCclusterMI *cl=0;
1679 const AliTPCRow &krow=GetRow(t.GetRelativeSector(),nr);
1680 if ( (t.GetSigmaY2()<0) || t.GetSigmaZ2()<0) return 0;
1682 Double_t roadz = 1.;
1685 if (TMath::Abs(TMath::Abs(y)-ymax)<krow.GetDeadZone()){
1687 t.SetClusterIndex2(row,-1);
1692 if (TMath::Abs(z)>(AliTPCReconstructor::GetCtgRange()*x+10)) t.SetClusterIndex2(row,-1);
1696 if ((cl==0)&&(krow)) {
1697 // cl = krow.FindNearest2(y+10,z,roady,roadz,index);
1698 cl = krow.FindNearest2(y,z,roady,roadz,index);
1700 if (cl) t.SetCurrentClusterIndex1(krow.GetIndex(index));
1704 t.SetCurrentCluster(cl);
1705 // Int_t accept = AcceptCluster(&t,t.fCurrentCluster,1.);
1707 t.SetClusterIndex2(row,index);
1708 t.SetClusterPointer(row, cl);
1716 //_________________________________________________________________________
1717 Bool_t AliTPCtrackerMI::GetTrackPoint(Int_t index, AliTrackPoint &p ) const
1719 // Get track space point by index
1720 // return false in case the cluster doesn't exist
1721 AliTPCclusterMI *cl = GetClusterMI(index);
1722 if (!cl) return kFALSE;
1723 Int_t sector = (index&0xff000000)>>24;
1724 // Int_t row = (index&0x00ff0000)>>16;
1726 // xyz[0] = fParam->GetPadRowRadii(sector,row);
1727 xyz[0] = cl->GetX();
1728 xyz[1] = cl->GetY();
1729 xyz[2] = cl->GetZ();
1731 fParam->AdjustCosSin(sector,cos,sin);
1732 Float_t x = cos*xyz[0]-sin*xyz[1];
1733 Float_t y = cos*xyz[1]+sin*xyz[0];
1735 Float_t sigmaY2 = 0.027*cl->GetSigmaY2();
1736 if (sector < fParam->GetNInnerSector()) sigmaY2 *= 2.07;
1737 Float_t sigmaZ2 = 0.066*cl->GetSigmaZ2();
1738 if (sector < fParam->GetNInnerSector()) sigmaZ2 *= 1.77;
1739 cov[0] = sin*sin*sigmaY2;
1740 cov[1] = -sin*cos*sigmaY2;
1742 cov[3] = cos*cos*sigmaY2;
1745 p.SetXYZ(x,y,xyz[2],cov);
1746 AliGeomManager::ELayerID iLayer;
1748 if (sector < fParam->GetNInnerSector()) {
1749 iLayer = AliGeomManager::kTPC1;
1753 iLayer = AliGeomManager::kTPC2;
1754 idet = sector - fParam->GetNInnerSector();
1756 UShort_t volid = AliGeomManager::LayerToVolUID(iLayer,idet);
1757 p.SetVolumeID(volid);
1763 Int_t AliTPCtrackerMI::UpdateClusters(AliTPCseed& t, Int_t nr) {
1764 //-----------------------------------------------------------------
1765 // This function tries to find a track prolongation to next pad row
1766 //-----------------------------------------------------------------
1767 t.SetCurrentCluster(0);
1768 t.SetCurrentClusterIndex1(0);
1770 Double_t xt=t.GetX();
1771 Int_t row = GetRowNumber(xt)-1;
1772 Double_t ymax= GetMaxY(nr);
1774 if (row < nr) return 1; // don't prolongate if not information until now -
1775 // if (TMath::Abs(t.GetSnp())>0.9 && t.GetNumberOfClusters()>40. && fIteration!=2) {
1777 // return 0; // not prolongate strongly inclined tracks
1779 // if (TMath::Abs(t.GetSnp())>0.95) {
1781 // return 0; // not prolongate strongly inclined tracks
1782 // }// patch 28 fev 06
1784 Double_t x= GetXrow(nr);
1786 //t.PropagateTo(x+0.02);
1787 //t.PropagateTo(x+0.01);
1788 if (!t.PropagateTo(x)){
1795 if (TMath::Abs(y)>ymax){
1797 t.SetRelativeSector((t.GetRelativeSector()+1) % fN);
1798 if (!t.Rotate(fSectors->GetAlpha()))
1800 } else if (y <-ymax) {
1801 t.SetRelativeSector((t.GetRelativeSector()-1+fN) % fN);
1802 if (!t.Rotate(-fSectors->GetAlpha()))
1805 // if (!t.PropagateTo(x)){
1812 if (TMath::Abs(t.GetSnp())>AliTPCReconstructor::GetMaxSnpTracker()) return 0;
1814 if (!IsActive(t.GetRelativeSector(),nr)) {
1816 t.SetClusterIndex2(nr,-1);
1819 //AliInfo(Form("A - Sector%d phi %f - alpha %f", t.fRelativeSector,y/x, t.GetAlpha()));
1821 AliTPCRow &krow=GetRow(t.GetRelativeSector(),nr);
1823 if (TMath::Abs(TMath::Abs(y)-ymax)<krow.GetDeadZone()){
1825 t.SetClusterIndex2(nr,-1);
1830 if (TMath::Abs(t.GetZ())<(AliTPCReconstructor::GetCtgRange()*t.GetX()+10) && (TMath::Abs(t.GetSnp())<AliTPCReconstructor::GetMaxSnpTracker())) t.SetNFoundable(t.GetNFoundable()+1);
1836 if ( (nr%6==0) || t.GetNumberOfClusters()<2){
1837 // t.fCurrentSigmaY = GetSigmaY(&t);
1838 //t.fCurrentSigmaZ = GetSigmaZ(&t);
1842 AliTPCclusterMI *cl=0;
1845 Double_t roady = 1.;
1846 Double_t roadz = 1.;
1850 index = t.GetClusterIndex2(nr);
1851 if ( (index>0) && (index&0x8000)==0){
1852 cl = t.GetClusterPointer(nr);
1853 if ( (cl==0) && (index>0)) cl = GetClusterMI(index);
1854 t.SetCurrentClusterIndex1(index);
1856 t.SetCurrentCluster(cl);
1862 // if (index<0) return 0;
1863 UInt_t uindex = TMath::Abs(index);
1866 //cl = krow.FindNearest2(y+10,z,roady,roadz,uindex);
1867 cl = krow.FindNearest2(y,z,roady,roadz,uindex);
1870 if (cl) t.SetCurrentClusterIndex1(krow.GetIndex(uindex));
1871 t.SetCurrentCluster(cl);
1877 Int_t AliTPCtrackerMI::FollowToNextCluster(AliTPCseed & t, Int_t nr) {
1878 //-----------------------------------------------------------------
1879 // This function tries to find a track prolongation to next pad row
1880 //-----------------------------------------------------------------
1882 //update error according neighborhoud
1884 if (t.GetCurrentCluster()) {
1886 Int_t accept = AcceptCluster(&t,t.GetCurrentCluster(),1.);
1888 if (t.GetCurrentCluster()->IsUsed(10)){
1893 t.SetNShared(t.GetNShared()+1);
1894 if (t.GetNShared()>0.7*t.GetNumberOfClusters()) {
1899 if (fIteration>0) accept = 0;
1900 if (accept<3) UpdateTrack(&t,accept);
1904 if ( ( (t.GetSigmaY2()+t.GetSigmaZ2())>0.16)&& t.GetNumberOfClusters()>18) t.SetRemoval(10);
1905 if ( t.GetChi2()/t.GetNumberOfClusters()>6 &&t.GetNumberOfClusters()>18) t.SetRemoval(10);
1907 if (( (t.GetNFoundable()*0.5 > t.GetNumberOfClusters()) || t.GetNoCluster()>15)) t.SetRemoval(10);
1915 //_____________________________________________________________________________
1916 Int_t AliTPCtrackerMI::FollowProlongation(AliTPCseed& t, Int_t rf, Int_t step) {
1917 //-----------------------------------------------------------------
1918 // This function tries to find a track prolongation.
1919 //-----------------------------------------------------------------
1920 Double_t xt=t.GetX();
1922 Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
1923 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
1924 if (alpha < 0. ) alpha += 2.*TMath::Pi();
1926 t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN);
1928 Int_t first = GetRowNumber(xt)-1;
1929 for (Int_t nr= first; nr>=rf; nr-=step) {
1931 if (t.GetKinkIndexes()[0]>0){
1932 for (Int_t i=0;i<3;i++){
1933 Int_t index = t.GetKinkIndexes()[i];
1934 if (index==0) break;
1935 if (index<0) continue;
1937 AliKink * kink = (AliKink*)fEvent->GetKink(index-1);
1939 printf("PROBLEM\n");
1942 Int_t kinkrow = kink->GetTPCRow0()+2+Int_t(0.5/(0.05+kink->GetAngle(2)));
1944 AliExternalTrackParam paramd(t);
1945 kink->SetDaughter(paramd);
1946 kink->SetStatus(2,5);
1953 if (nr==80) t.UpdateReference();
1954 if (nr<fInnerSec->GetNRows())
1955 fSectors = fInnerSec;
1957 fSectors = fOuterSec;
1958 if (FollowToNext(t,nr)==0)
1967 //_____________________________________________________________________________
1968 Int_t AliTPCtrackerMI::FollowProlongationFast(AliTPCseed& t, Int_t rf, Int_t step) {
1969 //-----------------------------------------------------------------
1970 // This function tries to find a track prolongation.
1971 //-----------------------------------------------------------------
1972 Double_t xt=t.GetX();
1974 Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
1975 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
1976 if (alpha < 0. ) alpha += 2.*TMath::Pi();
1977 t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN);
1979 for (Int_t nr=GetRowNumber(xt)-1; nr>=rf; nr-=step) {
1981 if (FollowToNextFast(t,nr)==0)
1982 if (!t.IsActive()) return 0;
1992 Int_t AliTPCtrackerMI::FollowBackProlongation(AliTPCseed& t, Int_t rf) {
1993 //-----------------------------------------------------------------
1994 // This function tries to find a track prolongation.
1995 //-----------------------------------------------------------------
1997 Double_t xt=t.GetX();
1998 Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
1999 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
2000 if (alpha < 0. ) alpha += 2.*TMath::Pi();
2001 t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN);
2003 Int_t first = t.GetFirstPoint();
2004 if (first<GetRowNumber(xt)+1) first = GetRowNumber(xt)+1;
2006 if (first<0) first=0;
2007 for (Int_t nr=first; nr<=rf; nr++) {
2008 // if ( (TMath::Abs(t.GetSnp())>0.95)) break;//patch 28 fev 06
2009 if (t.GetKinkIndexes()[0]<0){
2010 for (Int_t i=0;i<3;i++){
2011 Int_t index = t.GetKinkIndexes()[i];
2012 if (index==0) break;
2013 if (index>0) continue;
2014 index = TMath::Abs(index);
2015 AliKink * kink = (AliKink*)fEvent->GetKink(index-1);
2017 printf("PROBLEM\n");
2020 Int_t kinkrow = kink->GetTPCRow0()-2-Int_t(0.5/(0.05+kink->GetAngle(2)));
2022 AliExternalTrackParam paramm(t);
2023 kink->SetMother(paramm);
2024 kink->SetStatus(2,1);
2031 if (nr<fInnerSec->GetNRows())
2032 fSectors = fInnerSec;
2034 fSectors = fOuterSec;
2045 Float_t AliTPCtrackerMI::OverlapFactor(AliTPCseed * s1, AliTPCseed * s2, Int_t &sum1, Int_t & sum2)
2053 Float_t dz2 =(s1->GetZ() - s2->GetZ());
2056 Float_t dy2 =TMath::Abs((s1->GetY() - s2->GetY()));
2058 Float_t distance = TMath::Sqrt(dz2+dy2);
2059 if (distance>4.) return 0; // if there are far away - not overlap - to reduce combinatorics
2062 Int_t firstpoint = TMath::Min(s1->GetFirstPoint(),s2->GetFirstPoint());
2063 Int_t lastpoint = TMath::Max(s1->GetLastPoint(),s2->GetLastPoint());
2068 if (firstpoint>lastpoint) {
2069 firstpoint =lastpoint;
2074 for (Int_t i=firstpoint-1;i<lastpoint+1;i++){
2075 if (s1->GetClusterIndex2(i)>0) sum1++;
2076 if (s2->GetClusterIndex2(i)>0) sum2++;
2077 if (s1->GetClusterIndex2(i)==s2->GetClusterIndex2(i) && s1->GetClusterIndex2(i)>0) {
2081 if (sum<5) return 0;
2083 Float_t summin = TMath::Min(sum1+1,sum2+1);
2084 Float_t ratio = (sum+1)/Float_t(summin);
2088 void AliTPCtrackerMI::SignShared(AliTPCseed * s1, AliTPCseed * s2)
2092 Float_t thetaCut = 0.2;//+10.*TMath::Sqrt(s1->GetSigmaTglZ()+ s2->GetSigmaTglZ());
2093 if (TMath::Abs(s1->GetTgl()-s2->GetTgl())>thetaCut) return;
2094 Float_t minCl = TMath::Min(s1->GetNumberOfClusters(),s2->GetNumberOfClusters());
2095 Int_t cutN0 = TMath::Max(5,TMath::Nint(0.1*minCl));
2100 //Int_t firstpoint = TMath::Max(s1->GetFirstPoint(),s2->GetFirstPoint());
2101 //Int_t lastpoint = TMath::Min(s1->GetLastPoint(),s2->GetLastPoint());
2102 Int_t firstpoint = 0;
2103 Int_t lastpoint = 160;
2105 // if (firstpoint>=lastpoint-5) return;;
2107 for (Int_t i=firstpoint;i<lastpoint;i++){
2108 // if ( (s1->GetClusterIndex2(i)&0xFFFF8FFF)==(s2->GetClusterIndex2(i)&0xFFFF8FFF) && s1->GetClusterIndex2(i)>0) {
2109 if ( (s1->GetClusterIndex2(i))==(s2->GetClusterIndex2(i)) && s1->GetClusterIndex2(i)>0) {
2111 s1->SetSharedMapBit(i, kTRUE);
2112 s2->SetSharedMapBit(i, kTRUE);
2114 if (s1->GetClusterIndex2(i)>0)
2115 s1->SetClusterMapBit(i, kTRUE);
2117 if (sumshared>cutN0){
2120 for (Int_t i=firstpoint;i<lastpoint;i++){
2121 // if ( (s1->GetClusterIndex2(i)&0xFFFF8FFF)==(s2->GetClusterIndex2(i)&0xFFFF8FFF) && s1->GetClusterIndex2(i)>0) {
2122 if ( (s1->GetClusterIndex2(i))==(s2->GetClusterIndex2(i)) && s1->GetClusterIndex2(i)>0) {
2123 AliTPCTrackerPoint *p1 = s1->GetTrackPoint(i);
2124 AliTPCTrackerPoint *p2 = s2->GetTrackPoint(i);;
2125 if (s1->IsActive()&&s2->IsActive()){
2126 p1->SetShared(kTRUE);
2127 p2->SetShared(kTRUE);
2133 if (sumshared>cutN0){
2134 for (Int_t i=0;i<4;i++){
2135 if (s1->GetOverlapLabel(3*i)==0){
2136 s1->SetOverlapLabel(3*i, s2->GetLabel());
2137 s1->SetOverlapLabel(3*i+1,sumshared);
2138 s1->SetOverlapLabel(3*i+2,s2->GetUniqueID());
2142 for (Int_t i=0;i<4;i++){
2143 if (s2->GetOverlapLabel(3*i)==0){
2144 s2->SetOverlapLabel(3*i, s1->GetLabel());
2145 s2->SetOverlapLabel(3*i+1,sumshared);
2146 s2->SetOverlapLabel(3*i+2,s1->GetUniqueID());
2153 void AliTPCtrackerMI::SignShared(TObjArray * arr)
2156 //sort trackss according sectors
2158 for (Int_t i=0; i<arr->GetEntriesFast(); i++) {
2159 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2161 //if (pt) RotateToLocal(pt);
2165 arr->Sort(); // sorting according relative sectors
2166 arr->Expand(arr->GetEntries());
2169 Int_t nseed=arr->GetEntriesFast();
2170 for (Int_t i=0; i<nseed; i++) {
2171 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2173 for (Int_t j=0;j<=12;j++){
2174 pt->SetOverlapLabel(j,0);
2177 for (Int_t i=0; i<nseed; i++) {
2178 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2180 if (pt->GetRemoval()>10) continue;
2181 for (Int_t j=i+1; j<nseed; j++){
2182 AliTPCseed *pt2=(AliTPCseed*)arr->UncheckedAt(j);
2183 if (TMath::Abs(pt->GetRelativeSector()-pt2->GetRelativeSector())>1) continue;
2185 if (pt2->GetRemoval()<=10) {
2186 //if ( TMath::Abs(pt->GetRelativeSector()-pt2->GetRelativeSector())>0) break;
2194 void AliTPCtrackerMI::SortTracks(TObjArray * arr, Int_t mode) const
2197 //sort tracks in array according mode criteria
2198 Int_t nseed = arr->GetEntriesFast();
2199 for (Int_t i=0; i<nseed; i++) {
2200 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2211 void AliTPCtrackerMI::RemoveUsed2(TObjArray * arr, Float_t factor1, Float_t factor2, Int_t minimal)
2214 // Loop over all tracks and remove overlaped tracks (with lower quality)
2216 // 1. Unsign clusters
2217 // 2. Sort tracks according quality
2218 // Quality is defined by the number of cluster between first and last points
2220 // 3. Loop over tracks - decreasing quality order
2221 // a.) remove - If the fraction of shared cluster less than factor (1- n or 2)
2222 // b.) remove - If the minimal number of clusters less than minimal and not ITS
2223 // c.) if track accepted - sign clusters
2225 //Called in - AliTPCtrackerMI::Clusters2Tracks()
2226 // - AliTPCtrackerMI::PropagateBack()
2227 // - AliTPCtrackerMI::RefitInward()
2233 Int_t nseed = arr->GetEntriesFast();
2234 Float_t * quality = new Float_t[nseed];
2235 Int_t * indexes = new Int_t[nseed];
2239 for (Int_t i=0; i<nseed; i++) {
2240 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2245 pt->UpdatePoints(); //select first last max dens points
2246 Float_t * points = pt->GetPoints();
2247 if (points[3]<0.8) quality[i] =-1;
2248 quality[i] = (points[2]-points[0])+pt->GetNumberOfClusters();
2249 //prefer high momenta tracks if overlaps
2250 quality[i] *= TMath::Sqrt(TMath::Abs(pt->Pt())+0.5);
2252 TMath::Sort(nseed,quality,indexes);
2255 for (Int_t itrack=0; itrack<nseed; itrack++) {
2256 Int_t trackindex = indexes[itrack];
2257 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(trackindex);
2260 if (quality[trackindex]<0){
2262 delete arr->RemoveAt(trackindex);
2265 arr->RemoveAt(trackindex);
2271 Int_t first = Int_t(pt->GetPoints()[0]);
2272 Int_t last = Int_t(pt->GetPoints()[2]);
2273 Double_t factor = (pt->GetBConstrain()) ? factor1: factor2;
2275 Int_t found,foundable,shared;
2276 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
2277 // pt->GetClusterStatistic(0,160, found, foundable,shared,kFALSE);
2278 Bool_t itsgold =kFALSE;
2281 if (pt->GetESD()->GetITSclusters(dummy)>4) itsgold= kTRUE;
2285 if (Float_t(shared+1)/Float_t(found+1)>factor){
2286 if (pt->GetKinkIndexes()[0]!=0) continue; //don't remove tracks - part of the kinks
2287 delete arr->RemoveAt(trackindex);
2290 if (pt->GetNumberOfClusters()<50&&(found-0.5*shared)<minimal){ //remove short tracks
2291 if (pt->GetKinkIndexes()[0]!=0) continue; //don't remove tracks - part of the kinks
2292 delete arr->RemoveAt(trackindex);
2298 //if (sharedfactor>0.4) continue;
2299 if (pt->GetKinkIndexes()[0]>0) continue;
2300 //Remove tracks with undefined properties - seems
2301 if (pt->GetSigmaY2()<kAlmost0) continue; // ? what is the origin ?
2303 for (Int_t i=first; i<last; i++) {
2304 Int_t index=pt->GetClusterIndex2(i);
2305 // if (index<0 || index&0x8000 ) continue;
2306 if (index<0 || index&0x8000 ) continue;
2307 AliTPCclusterMI *c= pt->GetClusterPointer(i);
2314 Info("RemoveUsed","\n*****\nNumber of good tracks after shared removal\t%d\n",fNtracks);
2320 void AliTPCtrackerMI::UnsignClusters()
2323 // loop over all clusters and unsign them
2326 for (Int_t sec=0;sec<fkNIS;sec++){
2327 for (Int_t row=0;row<fInnerSec->GetNRows();row++){
2328 AliTPCclusterMI *cl = fInnerSec[sec][row].GetClusters1();
2329 for (Int_t icl =0;icl< fInnerSec[sec][row].GetN1();icl++)
2330 // if (cl[icl].IsUsed(10))
2332 cl = fInnerSec[sec][row].GetClusters2();
2333 for (Int_t icl =0;icl< fInnerSec[sec][row].GetN2();icl++)
2334 //if (cl[icl].IsUsed(10))
2339 for (Int_t sec=0;sec<fkNOS;sec++){
2340 for (Int_t row=0;row<fOuterSec->GetNRows();row++){
2341 AliTPCclusterMI *cl = fOuterSec[sec][row].GetClusters1();
2342 for (Int_t icl =0;icl< fOuterSec[sec][row].GetN1();icl++)
2343 //if (cl[icl].IsUsed(10))
2345 cl = fOuterSec[sec][row].GetClusters2();
2346 for (Int_t icl =0;icl< fOuterSec[sec][row].GetN2();icl++)
2347 //if (cl[icl].IsUsed(10))
2356 void AliTPCtrackerMI::SignClusters(TObjArray * arr, Float_t fnumber, Float_t fdensity)
2359 //sign clusters to be "used"
2361 // snumber and sdensity sign number of sigmas - bellow mean value to be accepted
2362 // loop over "primaries"
2376 Int_t nseed = arr->GetEntriesFast();
2377 for (Int_t i=0; i<nseed; i++) {
2378 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2382 if (!(pt->IsActive())) continue;
2383 Float_t dens = pt->GetNumberOfClusters()/Float_t(pt->GetNFoundable());
2384 if ( (dens>0.7) && (pt->GetNumberOfClusters()>70)){
2386 sumdens2+= dens*dens;
2387 sumn += pt->GetNumberOfClusters();
2388 sumn2 += pt->GetNumberOfClusters()*pt->GetNumberOfClusters();
2389 Float_t chi2 = pt->GetChi2()/pt->GetNumberOfClusters();
2392 sumchi2 +=chi2*chi2;
2397 Float_t mdensity = 0.9;
2398 Float_t meann = 130;
2399 Float_t meanchi = 1;
2400 Float_t sdensity = 0.1;
2401 Float_t smeann = 10;
2402 Float_t smeanchi =0.4;
2406 mdensity = sumdens/sum;
2408 meanchi = sumchi/sum;
2410 sdensity = sumdens2/sum-mdensity*mdensity;
2412 sdensity = TMath::Sqrt(sdensity);
2416 smeann = sumn2/sum-meann*meann;
2418 smeann = TMath::Sqrt(smeann);
2422 smeanchi = sumchi2/sum - meanchi*meanchi;
2424 smeanchi = TMath::Sqrt(smeanchi);
2430 //REMOVE SHORT DELTAS or tracks going out of sensitive volume of TPC
2432 for (Int_t i=0; i<nseed; i++) {
2433 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2437 if (pt->GetBSigned()) continue;
2438 if (pt->GetBConstrain()) continue;
2439 //if (!(pt->IsActive())) continue;
2441 Int_t found,foundable,shared;
2442 pt->GetClusterStatistic(0,160,found, foundable,shared);
2443 if (shared/float(found)>0.3) {
2444 if (shared/float(found)>0.9 ){
2445 //delete arr->RemoveAt(i);
2450 Bool_t isok =kFALSE;
2451 if ( (pt->GetNShared()/pt->GetNumberOfClusters()<0.5) &&pt->GetNumberOfClusters()>60)
2453 if ((TMath::Abs(1/pt->GetC())<100.) && (pt->GetNShared()/pt->GetNumberOfClusters()<0.7))
2455 if (TMath::Abs(pt->GetZ()/pt->GetX())>1.1)
2457 if ( (TMath::Abs(pt->GetSnp()>0.7) && pt->GetD(0,0)>60.))
2461 for (Int_t i=0; i<160; i++) {
2462 Int_t index=pt->GetClusterIndex2(i);
2463 if (index<0) continue;
2464 AliTPCclusterMI *c= pt->GetClusterPointer(i);
2466 //if (!(c->IsUsed(10))) c->Use();
2473 Double_t maxchi = meanchi+2.*smeanchi;
2475 for (Int_t i=0; i<nseed; i++) {
2476 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2480 //if (!(pt->IsActive())) continue;
2481 if (pt->GetBSigned()) continue;
2482 Double_t chi = pt->GetChi2()/pt->GetNumberOfClusters();
2483 if (chi>maxchi) continue;
2486 Float_t dens = pt->GetNumberOfClusters()/Float_t(pt->GetNFoundable());
2488 //sign only tracks with enoug big density at the beginning
2490 if ((pt->GetDensityFirst(40)<0.75) && pt->GetNumberOfClusters()<meann) continue;
2493 Double_t mindens = TMath::Max(double(mdensity-sdensity*fdensity*bfactor),0.65);
2494 Double_t minn = TMath::Max(Int_t(meann-fnumber*smeann*bfactor),50);
2496 // if (pt->fBConstrain) mindens = TMath::Max(mdensity-sdensity*fdensity*bfactor,0.65);
2497 if ( (pt->GetRemoval()==10) && (pt->GetSnp()>0.8)&&(dens>mindens))
2500 if ((dens>mindens && pt->GetNumberOfClusters()>minn) && chi<maxchi ){
2501 //Int_t noc=pt->GetNumberOfClusters();
2502 pt->SetBSigned(kTRUE);
2503 for (Int_t i=0; i<160; i++) {
2505 Int_t index=pt->GetClusterIndex2(i);
2506 if (index<0) continue;
2507 AliTPCclusterMI *c= pt->GetClusterPointer(i);
2509 // if (!(c->IsUsed(10))) c->Use();
2514 // gLastCheck = nseed;
2522 void AliTPCtrackerMI::StopNotActive(TObjArray * arr, Int_t row0, Float_t th0, Float_t th1, Float_t th2) const
2524 // stop not active tracks
2525 // take th1 as threshold for number of founded to number of foundable on last 10 active rows
2526 // take th2 as threshold for number of founded to number of foundable on last 20 active rows
2527 Int_t nseed = arr->GetEntriesFast();
2529 for (Int_t i=0; i<nseed; i++) {
2530 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2534 if (!(pt->IsActive())) continue;
2535 StopNotActive(pt,row0,th0, th1,th2);
2541 void AliTPCtrackerMI::StopNotActive(AliTPCseed * seed, Int_t row0, Float_t th0, Float_t th1,
2544 // stop not active tracks
2545 // take th1 as threshold for number of founded to number of foundable on last 10 active rows
2546 // take th2 as threshold for number of founded to number of foundable on last 20 active rows
2549 Int_t foundable = 0;
2550 Int_t maxindex = seed->GetLastPoint(); //last foundable row
2551 if (seed->GetNFoundable()*th0 > seed->GetNumberOfClusters()) {
2552 seed->Desactivate(10) ;
2556 for (Int_t i=row0; i<maxindex; i++){
2557 Int_t index = seed->GetClusterIndex2(i);
2558 if (index!=-1) foundable++;
2560 if (foundable<=30) sumgood1++;
2561 if (foundable<=50) {
2568 if (foundable>=30.){
2569 if (sumgood1<(th1*30.)) seed->Desactivate(10);
2572 if (sumgood2<(th2*50.)) seed->Desactivate(10);
2576 Int_t AliTPCtrackerMI::RefitInward(AliESDEvent *event)
2579 // back propagation of ESD tracks
2582 const Int_t kMaxFriendTracks=2000;
2586 //PrepareForProlongation(fSeeds,1);
2587 PropagateForward2(fSeeds);
2588 RemoveUsed2(fSeeds,0.4,0.4,20);
2590 TObjArray arraySeed(fSeeds->GetEntries());
2591 for (Int_t i=0;i<fSeeds->GetEntries();i++) {
2592 arraySeed.AddAt(fSeeds->At(i),i);
2594 SignShared(&arraySeed);
2595 // FindCurling(fSeeds, event,2); // find multi found tracks
2596 FindSplitted(fSeeds, event,2); // find multi found tracks
2599 Int_t nseed = fSeeds->GetEntriesFast();
2600 for (Int_t i=0;i<nseed;i++){
2601 AliTPCseed * seed = (AliTPCseed*) fSeeds->UncheckedAt(i);
2602 if (!seed) continue;
2603 if (seed->GetKinkIndex(0)>0) UpdateKinkQualityD(seed); // update quality informations for kinks
2605 seed->PropagateTo(fParam->GetInnerRadiusLow());
2606 seed->UpdatePoints();
2608 AliESDtrack *esd=event->GetTrack(i);
2609 seed->CookdEdx(0.02,0.6);
2610 CookLabel(seed,0.1); //For comparison only
2611 esd->SetTPCClusterMap(seed->GetClusterMap());
2612 esd->SetTPCSharedMap(seed->GetSharedMap());
2614 if (AliTPCReconstructor::StreamLevel()>1 && seed!=0&&esd!=0) {
2615 TTreeSRedirector &cstream = *fDebugStreamer;
2622 if (seed->GetNumberOfClusters()>15){
2623 esd->UpdateTrackParams(seed,AliESDtrack::kTPCrefit);
2624 esd->SetTPCPoints(seed->GetPoints());
2625 esd->SetTPCPointsF(seed->GetNFoundable());
2626 Int_t ndedx = seed->GetNCDEDX(0)+seed->GetNCDEDX(1)+seed->GetNCDEDX(2)+seed->GetNCDEDX(3);
2627 Float_t sdedx = (seed->GetSDEDX(0)+seed->GetSDEDX(1)+seed->GetSDEDX(2)+seed->GetSDEDX(3))*0.25;
2628 Float_t dedx = seed->GetdEdx();
2629 esd->SetTPCsignal(dedx, sdedx, ndedx);
2631 // add seed to the esd track in Calib level
2633 Bool_t storeFriend = gRandom->Rndm()<(kMaxFriendTracks)/Float_t(nseed);
2634 if (AliTPCReconstructor::StreamLevel()>0 &&storeFriend){
2635 AliTPCseed * seedCopy = new AliTPCseed(*seed, kTRUE);
2636 esd->AddCalibObject(seedCopy);
2641 //printf("problem\n");
2644 //FindKinks(fSeeds,event);
2645 Info("RefitInward","Number of refitted tracks %d",ntracks);
2652 Int_t AliTPCtrackerMI::PropagateBack(AliESDEvent *event)
2655 // back propagation of ESD tracks
2661 PropagateBack(fSeeds);
2662 RemoveUsed2(fSeeds,0.4,0.4,20);
2663 //FindCurling(fSeeds, fEvent,1);
2664 FindSplitted(fSeeds, event,1); // find multi found tracks
2667 Int_t nseed = fSeeds->GetEntriesFast();
2669 for (Int_t i=0;i<nseed;i++){
2670 AliTPCseed * seed = (AliTPCseed*) fSeeds->UncheckedAt(i);
2671 if (!seed) continue;
2672 if (seed->GetKinkIndex(0)<0) UpdateKinkQualityM(seed); // update quality informations for kinks
2673 seed->UpdatePoints();
2674 AliESDtrack *esd=event->GetTrack(i);
2675 seed->CookdEdx(0.02,0.6);
2676 CookLabel(seed,0.1); //For comparison only
2677 if (seed->GetNumberOfClusters()>15){
2678 esd->UpdateTrackParams(seed,AliESDtrack::kTPCout);
2679 esd->SetTPCPoints(seed->GetPoints());
2680 esd->SetTPCPointsF(seed->GetNFoundable());
2681 Int_t ndedx = seed->GetNCDEDX(0)+seed->GetNCDEDX(1)+seed->GetNCDEDX(2)+seed->GetNCDEDX(3);
2682 Float_t sdedx = (seed->GetSDEDX(0)+seed->GetSDEDX(1)+seed->GetSDEDX(2)+seed->GetSDEDX(3))*0.25;
2683 Float_t dedx = seed->GetdEdx();
2684 esd->SetTPCsignal(dedx, sdedx, ndedx);
2686 Int_t eventnumber = event->GetEventNumberInFile();// patch 28 fev 06
2687 // This is most likely NOT the event number you'd like to use. It has nothing to do with the 'real' event number
2688 if (AliTPCReconstructor::StreamLevel()>1) {
2689 (*fDebugStreamer)<<"Cback"<<
2691 "EventNrInFile="<<eventnumber<<
2692 "\n"; // patch 28 fev 06
2696 //FindKinks(fSeeds,event);
2697 Info("PropagateBack","Number of back propagated tracks %d",ntracks);
2705 void AliTPCtrackerMI::DeleteSeeds()
2710 Int_t nseed = fSeeds->GetEntriesFast();
2711 for (Int_t i=0;i<nseed;i++){
2712 AliTPCseed * seed = (AliTPCseed*)fSeeds->At(i);
2713 if (seed) delete fSeeds->RemoveAt(i);
2720 void AliTPCtrackerMI::ReadSeeds(AliESDEvent *event, Int_t direction)
2723 //read seeds from the event
2725 Int_t nentr=event->GetNumberOfTracks();
2727 Info("ReadSeeds", "Number of ESD tracks: %d\n", nentr);
2732 fSeeds = new TObjArray(nentr);
2736 for (Int_t i=0; i<nentr; i++) {
2737 AliESDtrack *esd=event->GetTrack(i);
2738 ULong_t status=esd->GetStatus();
2739 if (!(status&AliESDtrack::kTPCin)) continue;
2740 AliTPCtrack t(*esd);
2741 t.SetNumberOfClusters(0);
2742 // AliTPCseed *seed = new AliTPCseed(t,t.GetAlpha());
2743 AliTPCseed *seed = new AliTPCseed(t/*,t.GetAlpha()*/);
2744 seed->SetUniqueID(esd->GetID());
2745 for (Int_t ikink=0;ikink<3;ikink++) {
2746 Int_t index = esd->GetKinkIndex(ikink);
2747 seed->GetKinkIndexes()[ikink] = index;
2748 if (index==0) continue;
2749 index = TMath::Abs(index);
2750 AliESDkink * kink = fEvent->GetKink(index-1);
2751 if (kink&&esd->GetKinkIndex(ikink)<0){
2752 if ((status & AliESDtrack::kTRDrefit) != 0) kink->SetStatus(1,2);
2753 if ((status & AliESDtrack::kITSout) != 0) kink->SetStatus(1,0);
2755 if (kink&&esd->GetKinkIndex(ikink)>0){
2756 if ((status & AliESDtrack::kTRDrefit) != 0) kink->SetStatus(1,6);
2757 if ((status & AliESDtrack::kITSout) != 0) kink->SetStatus(1,4);
2761 if (((status&AliESDtrack::kITSout)==0)&&(direction==1)) seed->ResetCovariance(10.);
2762 if ( direction ==2 &&(status & AliESDtrack::kTRDrefit) == 0 ) seed->ResetCovariance(10.);
2763 if ( direction ==2 && ((status & AliESDtrack::kTPCout) == 0) ) {
2768 if ( direction ==2 &&(status & AliESDtrack::kTRDrefit) > 0 ) {
2769 Double_t par0[5],par1[5],alpha,x;
2770 esd->GetInnerExternalParameters(alpha,x,par0);
2771 esd->GetExternalParameters(x,par1);
2772 Double_t delta1 = TMath::Abs(par0[4]-par1[4])/(0.000000001+TMath::Abs(par0[4]+par1[4]));
2773 Double_t delta2 = TMath::Abs(par0[3]-par1[3]);
2775 if (esd->GetTRDncls()>0) trdchi2 = esd->GetTRDchi2()/esd->GetTRDncls();
2776 //reset covariance if suspicious
2777 if ( (delta1>0.1) || (delta2>0.006) ||trdchi2>7.)
2778 seed->ResetCovariance(10.);
2783 // rotate to the local coordinate system
2785 fSectors=fInnerSec; fN=fkNIS;
2786 Double_t alpha=seed->GetAlpha() - fSectors->GetAlphaShift();
2787 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
2788 if (alpha < 0. ) alpha += 2.*TMath::Pi();
2789 Int_t ns=Int_t(alpha/fSectors->GetAlpha())%fN;
2790 alpha =ns*fSectors->GetAlpha() + fSectors->GetAlphaShift();
2791 if (alpha<-TMath::Pi()) alpha += 2*TMath::Pi();
2792 if (alpha>=TMath::Pi()) alpha -= 2*TMath::Pi();
2793 alpha-=seed->GetAlpha();
2794 if (!seed->Rotate(alpha)) {
2800 if (esd->GetKinkIndex(0)<=0){
2801 for (Int_t irow=0;irow<160;irow++){
2802 Int_t index = seed->GetClusterIndex2(irow);
2805 AliTPCclusterMI * cl = GetClusterMI(index);
2806 seed->SetClusterPointer(irow,cl);
2808 if ((index & 0x8000)==0){
2809 cl->Use(10); // accepted cluster
2811 cl->Use(6); // close cluster not accepted
2814 Info("ReadSeeds","Not found cluster");
2819 fSeeds->AddAt(seed,i);
2825 //_____________________________________________________________________________
2826 void AliTPCtrackerMI::MakeSeeds3(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, Float_t cuts[4],
2827 Float_t deltay, Int_t ddsec) {
2828 //-----------------------------------------------------------------
2829 // This function creates track seeds.
2830 // SEEDING WITH VERTEX CONSTRAIN
2831 //-----------------------------------------------------------------
2832 // cuts[0] - fP4 cut
2833 // cuts[1] - tan(phi) cut
2834 // cuts[2] - zvertex cut
2835 // cuts[3] - fP3 cut
2843 Double_t x[5], c[15];
2844 // Int_t di = i1-i2;
2846 AliTPCseed * seed = new AliTPCseed();
2847 Double_t alpha=fSectors->GetAlpha(), shift=fSectors->GetAlphaShift();
2848 Double_t cs=cos(alpha), sn=sin(alpha);
2850 // Double_t x1 =fOuterSec->GetX(i1);
2851 //Double_t xx2=fOuterSec->GetX(i2);
2853 Double_t x1 =GetXrow(i1);
2854 Double_t xx2=GetXrow(i2);
2856 Double_t x3=GetX(), y3=GetY(), z3=GetZ();
2858 Int_t imiddle = (i2+i1)/2; //middle pad row index
2859 Double_t xm = GetXrow(imiddle); // radius of middle pad-row
2860 const AliTPCRow& krm=GetRow(sec,imiddle); //middle pad -row
2864 const AliTPCRow& kr1=GetRow(ns,i1);
2865 Double_t ymax = GetMaxY(i1)-kr1.GetDeadZone()-1.5;
2866 Double_t ymaxm = GetMaxY(imiddle)-kr1.GetDeadZone()-1.5;
2869 // change cut on curvature if it can't reach this layer
2870 // maximal curvature set to reach it
2871 Double_t dvertexmax = TMath::Sqrt((x1-x3)*(x1-x3)+(ymax+5-y3)*(ymax+5-y3));
2872 if (dvertexmax*0.5*cuts[0]>0.85){
2873 cuts[0] = 0.85/(dvertexmax*0.5+1.);
2875 Double_t r2min = 1/(cuts[0]*cuts[0]); //minimal square of radius given by cut
2878 if (deltay>0) ddsec = 0;
2879 // loop over clusters
2880 for (Int_t is=0; is < kr1; is++) {
2882 if (kr1[is]->IsUsed(10)) continue;
2883 Double_t y1=kr1[is]->GetY(), z1=kr1[is]->GetZ();
2884 //if (TMath::Abs(y1)>ymax) continue;
2886 if (deltay>0 && TMath::Abs(ymax-TMath::Abs(y1))> deltay ) continue; // seed only at the edge
2888 // find possible directions
2889 Float_t anglez = (z1-z3)/(x1-x3);
2890 Float_t extraz = z1 - anglez*(x1-xx2); // extrapolated z
2893 //find rotation angles relative to line given by vertex and point 1
2894 Double_t dvertex2 = (x1-x3)*(x1-x3)+(y1-y3)*(y1-y3);
2895 Double_t dvertex = TMath::Sqrt(dvertex2);
2896 Double_t angle13 = TMath::ATan((y1-y3)/(x1-x3));
2897 Double_t cs13 = cos(-angle13), sn13 = sin(-angle13);
2900 // loop over 2 sectors
2906 Double_t dddz1=0; // direction of delta inclination in z axis
2913 for (Int_t dsec = dsec1; dsec<=dsec2;dsec++){
2914 Int_t sec2 = sec + dsec;
2916 // AliTPCRow& kr2 = fOuterSec[(sec2+fkNOS)%fkNOS][i2];
2917 //AliTPCRow& kr2m = fOuterSec[(sec2+fkNOS)%fkNOS][imiddle];
2918 AliTPCRow& kr2 = GetRow((sec2+fkNOS)%fkNOS,i2);
2919 AliTPCRow& kr2m = GetRow((sec2+fkNOS)%fkNOS,imiddle);
2920 Int_t index1 = TMath::Max(kr2.Find(extraz-0.6-dddz1*TMath::Abs(z1)*0.05)-1,0);
2921 Int_t index2 = TMath::Min(kr2.Find(extraz+0.6+dddz2*TMath::Abs(z1)*0.05)+1,kr2);
2923 // rotation angles to p1-p3
2924 Double_t cs13r = cos(-angle13+dsec*alpha)/dvertex, sn13r = sin(-angle13+dsec*alpha)/dvertex;
2925 Double_t x2, y2, z2;
2927 // Double_t dymax = maxangle*TMath::Abs(x1-xx2);
2930 Double_t dxx0 = (xx2-x3)*cs13r;
2931 Double_t dyy0 = (xx2-x3)*sn13r;
2932 for (Int_t js=index1; js < index2; js++) {
2933 const AliTPCclusterMI *kcl = kr2[js];
2934 if (kcl->IsUsed(10)) continue;
2936 //calcutate parameters
2938 Double_t yy0 = dyy0 +(kcl->GetY()-y3)*cs13r;
2940 if (TMath::Abs(yy0)<0.000001) continue;
2941 Double_t xx0 = dxx0 -(kcl->GetY()-y3)*sn13r;
2942 Double_t y0 = 0.5*(xx0*xx0+yy0*yy0-xx0)/yy0;
2943 Double_t r02 = (0.25+y0*y0)*dvertex2;
2944 //curvature (radius) cut
2945 if (r02<r2min) continue;
2949 Double_t c0 = 1/TMath::Sqrt(r02);
2953 //Double_t dfi0 = 2.*TMath::ASin(dvertex*c0*0.5);
2954 //Double_t dfi1 = 2.*TMath::ASin(TMath::Sqrt(yy0*yy0+(1-xx0)*(1-xx0))*dvertex*c0*0.5);
2955 Double_t dfi0 = 2.*AliTPCFastMath::FastAsin(dvertex*c0*0.5);
2956 Double_t dfi1 = 2.*AliTPCFastMath::FastAsin(TMath::Sqrt(yy0*yy0+(1-xx0)*(1-xx0))*dvertex*c0*0.5);
2959 Double_t z0 = kcl->GetZ();
2960 Double_t zzzz2 = z1-(z1-z3)*dfi1/dfi0;
2961 if (TMath::Abs(zzzz2-z0)>0.5) continue;
2964 Double_t dip = (z1-z0)*c0/dfi1;
2965 Double_t x0 = (0.5*cs13+y0*sn13)*dvertex*c0;
2976 x2= xx2*cs-y2*sn*dsec;
2977 y2=+xx2*sn*dsec+y2*cs;
2987 // do we have cluster at the middle ?
2989 GetProlongation(x1,xm,x,ym,zm);
2991 AliTPCclusterMI * cm=0;
2992 if (TMath::Abs(ym)-ymaxm<0){
2993 cm = krm.FindNearest2(ym,zm,1.0,0.6,dummy);
2994 if ((!cm) || (cm->IsUsed(10))) {
2999 // rotate y1 to system 0
3000 // get state vector in rotated system
3001 Double_t yr1 = (-0.5*sn13+y0*cs13)*dvertex*c0;
3002 Double_t xr2 = x0*cs+yr1*sn*dsec;
3003 Double_t xr[5]={kcl->GetY(),kcl->GetZ(), xr2, dip, c0};
3005 GetProlongation(xx2,xm,xr,ym,zm);
3006 if (TMath::Abs(ym)-ymaxm<0){
3007 cm = kr2m.FindNearest2(ym,zm,1.0,0.6,dummy);
3008 if ((!cm) || (cm->IsUsed(10))) {
3018 dym = ym - cm->GetY();
3019 dzm = zm - cm->GetZ();
3026 Double_t sy1=kr1[is]->GetSigmaY2()*2., sz1=kr1[is]->GetSigmaZ2()*2.;
3027 Double_t sy2=kcl->GetSigmaY2()*2., sz2=kcl->GetSigmaZ2()*2.;
3028 //Double_t sy3=400*3./12., sy=0.1, sz=0.1;
3029 Double_t sy3=25000*x[4]*x[4]+0.1, sy=0.1, sz=0.1;
3030 //Double_t sy3=25000*x[4]*x[4]*60+0.5, sy=0.1, sz=0.1;
3032 Double_t f40=(F1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
3033 Double_t f42=(F1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
3034 Double_t f43=(F1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
3035 Double_t f20=(F2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
3036 Double_t f22=(F2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
3037 Double_t f23=(F2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
3039 Double_t f30=(F3(x1,y1+sy,x2,y2,z1,z2)-x[3])/sy;
3040 Double_t f31=(F3(x1,y1,x2,y2,z1+sz,z2)-x[3])/sz;
3041 Double_t f32=(F3(x1,y1,x2,y2+sy,z1,z2)-x[3])/sy;
3042 Double_t f34=(F3(x1,y1,x2,y2,z1,z2+sz)-x[3])/sz;
3046 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3047 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3048 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3049 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3050 c[13]=f30*sy1*f40+f32*sy2*f42;
3051 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3053 // if (!BuildSeed(kr1[is],kcl,0,x1,x2,x3,x,c)) continue;
3055 UInt_t index=kr1.GetIndex(is);
3056 seed->~AliTPCseed(); // this does not set the pointer to 0...
3057 AliTPCseed *track=new(seed) AliTPCseed(x1, ns*alpha+shift, x, c, index);
3059 track->SetIsSeeding(kTRUE);
3060 track->SetSeed1(i1);
3061 track->SetSeed2(i2);
3062 track->SetSeedType(3);
3066 FollowProlongation(*track, (i1+i2)/2,1);
3067 Int_t foundable,found,shared;
3068 track->GetClusterStatistic((i1+i2)/2,i1, found, foundable, shared, kTRUE);
3069 if ((found<0.55*foundable) || shared>0.5*found || (track->GetSigmaY2()+track->GetSigmaZ2())>0.5){
3071 seed->~AliTPCseed();
3077 FollowProlongation(*track, i2,1);
3081 track->SetBConstrain(1);
3082 // track->fLastPoint = i1+fInnerSec->GetNRows(); // first cluster in track position
3083 track->SetLastPoint(i1); // first cluster in track position
3084 track->SetFirstPoint(track->GetLastPoint());
3086 if (track->GetNumberOfClusters()<(i1-i2)*0.5 ||
3087 track->GetNumberOfClusters() < track->GetNFoundable()*0.6 ||
3088 track->GetNShared()>0.4*track->GetNumberOfClusters() ) {
3090 seed->~AliTPCseed();
3094 // Z VERTEX CONDITION
3095 Double_t zv, bz=GetBz();
3096 if ( !track->GetZAt(0.,bz,zv) ) continue;
3097 if (TMath::Abs(zv-z3)>cuts[2]) {
3098 FollowProlongation(*track, TMath::Max(i2-20,0));
3099 if ( !track->GetZAt(0.,bz,zv) ) continue;
3100 if (TMath::Abs(zv-z3)>cuts[2]){
3101 FollowProlongation(*track, TMath::Max(i2-40,0));
3102 if ( !track->GetZAt(0.,bz,zv) ) continue;
3103 if (TMath::Abs(zv-z3)>cuts[2] &&(track->GetNumberOfClusters() > track->GetNFoundable()*0.7)){
3104 // make seed without constrain
3105 AliTPCseed * track2 = MakeSeed(track,0.2,0.5,1.);
3106 FollowProlongation(*track2, i2,1);
3107 track2->SetBConstrain(kFALSE);
3108 track2->SetSeedType(1);
3109 arr->AddLast(track2);
3111 seed->~AliTPCseed();
3116 seed->~AliTPCseed();
3123 track->SetSeedType(0);
3124 arr->AddLast(track);
3125 seed = new AliTPCseed;
3127 // don't consider other combinations
3128 if (track->GetNumberOfClusters() > track->GetNFoundable()*0.8)
3134 Info("MakeSeeds3","\nSeeding statistic:\t%d\t%d\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin,nout1,nout2);
3140 void AliTPCtrackerMI::MakeSeeds5(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, Float_t cuts[4],
3145 //-----------------------------------------------------------------
3146 // This function creates track seeds.
3147 //-----------------------------------------------------------------
3148 // cuts[0] - fP4 cut
3149 // cuts[1] - tan(phi) cut
3150 // cuts[2] - zvertex cut
3151 // cuts[3] - fP3 cut
3161 Double_t x[5], c[15];
3163 // make temporary seed
3164 AliTPCseed * seed = new AliTPCseed;
3165 Double_t alpha=fOuterSec->GetAlpha(), shift=fOuterSec->GetAlphaShift();
3166 // Double_t cs=cos(alpha), sn=sin(alpha);
3171 Double_t x1 = GetXrow(i1-1);
3172 const AliTPCRow& kr1=GetRow(sec,i1-1);
3173 Double_t y1max = GetMaxY(i1-1)-kr1.GetDeadZone()-1.5;
3175 Double_t x1p = GetXrow(i1);
3176 const AliTPCRow& kr1p=GetRow(sec,i1);
3178 Double_t x1m = GetXrow(i1-2);
3179 const AliTPCRow& kr1m=GetRow(sec,i1-2);
3182 //last 3 padrow for seeding
3183 AliTPCRow& kr3 = GetRow((sec+fkNOS)%fkNOS,i1-7);
3184 Double_t x3 = GetXrow(i1-7);
3185 // Double_t y3max= GetMaxY(i1-7)-kr3.fDeadZone-1.5;
3187 AliTPCRow& kr3p = GetRow((sec+fkNOS)%fkNOS,i1-6);
3188 Double_t x3p = GetXrow(i1-6);
3190 AliTPCRow& kr3m = GetRow((sec+fkNOS)%fkNOS,i1-8);
3191 Double_t x3m = GetXrow(i1-8);
3196 Int_t im = i1-4; //middle pad row index
3197 Double_t xm = GetXrow(im); // radius of middle pad-row
3198 const AliTPCRow& krm=GetRow(sec,im); //middle pad -row
3199 // Double_t ymmax = GetMaxY(im)-kr1.fDeadZone-1.5;
3202 Double_t deltax = x1-x3;
3203 Double_t dymax = deltax*cuts[1];
3204 Double_t dzmax = deltax*cuts[3];
3206 // loop over clusters
3207 for (Int_t is=0; is < kr1; is++) {
3209 if (kr1[is]->IsUsed(10)) continue;
3210 Double_t y1=kr1[is]->GetY(), z1=kr1[is]->GetZ();
3212 if (deltay>0 && TMath::Abs(y1max-TMath::Abs(y1))> deltay ) continue; // seed only at the edge
3214 Int_t index1 = TMath::Max(kr3.Find(z1-dzmax)-1,0);
3215 Int_t index2 = TMath::Min(kr3.Find(z1+dzmax)+1,kr3);
3221 for (Int_t js=index1; js < index2; js++) {
3222 const AliTPCclusterMI *kcl = kr3[js];
3223 if (kcl->IsUsed(10)) continue;
3225 // apply angular cuts
3226 if (TMath::Abs(y1-y3)>dymax) continue;
3229 if (TMath::Abs(z1-z3)>dzmax) continue;
3231 Double_t angley = (y1-y3)/(x1-x3);
3232 Double_t anglez = (z1-z3)/(x1-x3);
3234 Double_t erry = TMath::Abs(angley)*(x1-x1m)*0.5+0.5;
3235 Double_t errz = TMath::Abs(anglez)*(x1-x1m)*0.5+0.5;
3237 Double_t yyym = angley*(xm-x1)+y1;
3238 Double_t zzzm = anglez*(xm-x1)+z1;
3240 const AliTPCclusterMI *kcm = krm.FindNearest2(yyym,zzzm,erry,errz,index);
3242 if (kcm->IsUsed(10)) continue;
3244 erry = TMath::Abs(angley)*(x1-x1m)*0.4+0.5;
3245 errz = TMath::Abs(anglez)*(x1-x1m)*0.4+0.5;
3252 // look around first
3253 const AliTPCclusterMI *kc1m = kr1m.FindNearest2(angley*(x1m-x1)+y1,
3259 if (kc1m->IsUsed(10)) used++;
3261 const AliTPCclusterMI *kc1p = kr1p.FindNearest2(angley*(x1p-x1)+y1,
3267 if (kc1p->IsUsed(10)) used++;
3269 if (used>1) continue;
3270 if (found<1) continue;
3274 const AliTPCclusterMI *kc3m = kr3m.FindNearest2(angley*(x3m-x3)+y3,
3280 if (kc3m->IsUsed(10)) used++;
3284 const AliTPCclusterMI *kc3p = kr3p.FindNearest2(angley*(x3p-x3)+y3,
3290 if (kc3p->IsUsed(10)) used++;
3294 if (used>1) continue;
3295 if (found<3) continue;
3305 x[4]=F1(x1,y1,x2,y2,x3,y3);
3306 //if (TMath::Abs(x[4]) >= cuts[0]) continue;
3309 x[2]=F2(x1,y1,x2,y2,x3,y3);
3312 x[3]=F3n(x1,y1,x2,y2,z1,z2,x[4]);
3313 //if (TMath::Abs(x[3]) > cuts[3]) continue;
3317 Double_t sy1=0.1, sz1=0.1;
3318 Double_t sy2=0.1, sz2=0.1;
3319 Double_t sy3=0.1, sy=0.1, sz=0.1;
3321 Double_t f40=(F1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
3322 Double_t f42=(F1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
3323 Double_t f43=(F1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
3324 Double_t f20=(F2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
3325 Double_t f22=(F2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
3326 Double_t f23=(F2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
3328 Double_t f30=(F3(x1,y1+sy,x2,y2,z1,z2)-x[3])/sy;
3329 Double_t f31=(F3(x1,y1,x2,y2,z1+sz,z2)-x[3])/sz;
3330 Double_t f32=(F3(x1,y1,x2,y2+sy,z1,z2)-x[3])/sy;
3331 Double_t f34=(F3(x1,y1,x2,y2,z1,z2+sz)-x[3])/sz;
3335 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3336 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3337 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3338 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3339 c[13]=f30*sy1*f40+f32*sy2*f42;
3340 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3342 // if (!BuildSeed(kr1[is],kcl,0,x1,x2,x3,x,c)) continue;
3344 UInt_t index=kr1.GetIndex(is);
3345 seed->~AliTPCseed();
3346 AliTPCseed *track=new(seed) AliTPCseed(x1, sec*alpha+shift, x, c, index);
3348 track->SetIsSeeding(kTRUE);
3351 FollowProlongation(*track, i1-7,1);
3352 if (track->GetNumberOfClusters() < track->GetNFoundable()*0.75 ||
3353 track->GetNShared()>0.6*track->GetNumberOfClusters() || ( track->GetSigmaY2()+ track->GetSigmaZ2())>0.6){
3355 seed->~AliTPCseed();
3361 FollowProlongation(*track, i2,1);
3362 track->SetBConstrain(0);
3363 track->SetLastPoint(i1+fInnerSec->GetNRows()); // first cluster in track position
3364 track->SetFirstPoint(track->GetLastPoint());
3366 if (track->GetNumberOfClusters()<(i1-i2)*0.5 ||
3367 track->GetNumberOfClusters()<track->GetNFoundable()*0.7 ||
3368 track->GetNShared()>2. || track->GetChi2()/track->GetNumberOfClusters()>6 || ( track->GetSigmaY2()+ track->GetSigmaZ2())>0.5 ) {
3370 seed->~AliTPCseed();
3375 FollowProlongation(*track, TMath::Max(i2-10,0),1);
3376 AliTPCseed * track2 = MakeSeed(track,0.2,0.5,0.9);
3377 FollowProlongation(*track2, i2,1);
3378 track2->SetBConstrain(kFALSE);
3379 track2->SetSeedType(4);
3380 arr->AddLast(track2);
3382 seed->~AliTPCseed();
3386 //arr->AddLast(track);
3387 //seed = new AliTPCseed;
3393 Info("MakeSeeds5","\nSeeding statiistic:\t%d\t%d\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin,nout1,nout2,nout3);
3399 //_____________________________________________________________________________
3400 void AliTPCtrackerMI::MakeSeeds2(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, Float_t */*cuts[4]*/,
3401 Float_t deltay, Bool_t /*bconstrain*/) {
3402 //-----------------------------------------------------------------
3403 // This function creates track seeds - without vertex constraint
3404 //-----------------------------------------------------------------
3405 // cuts[0] - fP4 cut - not applied
3406 // cuts[1] - tan(phi) cut
3407 // cuts[2] - zvertex cut - not applied
3408 // cuts[3] - fP3 cut
3418 Double_t alpha=fOuterSec->GetAlpha(), shift=fOuterSec->GetAlphaShift();
3419 // Double_t cs=cos(alpha), sn=sin(alpha);
3420 Int_t row0 = (i1+i2)/2;
3421 Int_t drow = (i1-i2)/2;
3422 const AliTPCRow& kr0=fSectors[sec][row0];
3425 AliTPCpolyTrack polytrack;
3426 Int_t nclusters=fSectors[sec][row0];
3427 AliTPCseed * seed = new AliTPCseed;
3432 for (Int_t is=0; is < nclusters; is++) { //LOOP over clusters
3434 Int_t nfoundable =0;
3435 for (Int_t iter =1; iter<2; iter++){ //iterations
3436 const AliTPCRow& krm=fSectors[sec][row0-iter];
3437 const AliTPCRow& krp=fSectors[sec][row0+iter];
3438 const AliTPCclusterMI * cl= kr0[is];
3440 if (cl->IsUsed(10)) {
3446 Double_t x = kr0.GetX();
3447 // Initialization of the polytrack
3452 Double_t y0= cl->GetY();
3453 Double_t z0= cl->GetZ();
3457 Double_t ymax = fSectors->GetMaxY(row0)-kr0.GetDeadZone()-1.5;
3458 if (deltay>0 && TMath::Abs(ymax-TMath::Abs(y0))> deltay ) continue; // seed only at the edge
3460 erry = (0.5)*cl->GetSigmaY2()/TMath::Sqrt(cl->GetQ())*6;
3461 errz = (0.5)*cl->GetSigmaZ2()/TMath::Sqrt(cl->GetQ())*6;
3462 polytrack.AddPoint(x,y0,z0,erry, errz);
3465 if (cl->IsUsed(10)) sumused++;
3468 Float_t roady = (5*TMath::Sqrt(cl->GetSigmaY2()+0.2)+1.)*iter;
3469 Float_t roadz = (5*TMath::Sqrt(cl->GetSigmaZ2()+0.2)+1.)*iter;
3472 AliTPCclusterMI * cl1 = krm.FindNearest(y0,z0,roady,roadz);
3473 if (cl1 && TMath::Abs(ymax-TMath::Abs(y0))) {
3474 erry = (0.5)*cl1->GetSigmaY2()/TMath::Sqrt(cl1->GetQ())*3;
3475 errz = (0.5)*cl1->GetSigmaZ2()/TMath::Sqrt(cl1->GetQ())*3;
3476 if (cl1->IsUsed(10)) sumused++;
3477 polytrack.AddPoint(x,cl1->GetY(),cl1->GetZ(),erry,errz);
3481 AliTPCclusterMI * cl2 = krp.FindNearest(y0,z0,roady,roadz);
3483 erry = (0.5)*cl2->GetSigmaY2()/TMath::Sqrt(cl2->GetQ())*3;
3484 errz = (0.5)*cl2->GetSigmaZ2()/TMath::Sqrt(cl2->GetQ())*3;
3485 if (cl2->IsUsed(10)) sumused++;
3486 polytrack.AddPoint(x,cl2->GetY(),cl2->GetZ(),erry,errz);
3489 if (sumused>0) continue;
3491 polytrack.UpdateParameters();
3497 nfoundable = polytrack.GetN();
3498 nfound = nfoundable;
3500 for (Int_t ddrow = iter+1; ddrow<drow;ddrow++){
3501 Float_t maxdist = 0.8*(1.+3./(ddrow));
3502 for (Int_t delta = -1;delta<=1;delta+=2){
3503 Int_t row = row0+ddrow*delta;
3504 kr = &(fSectors[sec][row]);
3505 Double_t xn = kr->GetX();
3506 Double_t ymax = fSectors->GetMaxY(row)-kr->GetDeadZone()-1.5;
3507 polytrack.GetFitPoint(xn,yn,zn);
3508 if (TMath::Abs(yn)>ymax) continue;
3510 AliTPCclusterMI * cln = kr->FindNearest(yn,zn,roady,roadz);
3512 Float_t dist = TMath::Sqrt( (yn-cln->GetY())*(yn-cln->GetY())+(zn-cln->GetZ())*(zn-cln->GetZ()));
3515 erry = (dist+0.3)*cln->GetSigmaY2()/TMath::Sqrt(cln->GetQ())*(1.+1./(ddrow));
3516 errz = (dist+0.3)*cln->GetSigmaZ2()/TMath::Sqrt(cln->GetQ())*(1.+1./(ddrow));
3517 if (cln->IsUsed(10)) {
3518 // printf("used\n");
3526 polytrack.AddPoint(xn,cln->GetY(),cln->GetZ(),erry, errz);
3531 if ( (sumused>3) || (sumused>0.5*nfound) || (nfound<0.6*nfoundable)) break;
3532 polytrack.UpdateParameters();
3535 if ( (sumused>3) || (sumused>0.5*nfound)) {
3536 //printf("sumused %d\n",sumused);
3541 polytrack.GetFitDerivation(kr0.GetX(),dy,dz);
3542 AliTPCpolyTrack track2;
3544 polytrack.Refit(track2,0.5+TMath::Abs(dy)*0.3,0.4+TMath::Abs(dz)*0.3);
3545 if (track2.GetN()<0.5*nfoundable) continue;
3548 if ((nfound>0.6*nfoundable) &&( nfoundable>0.4*(i1-i2))) {
3550 // test seed with and without constrain
3551 for (Int_t constrain=0; constrain<=0;constrain++){
3552 // add polytrack candidate
3554 Double_t x[5], c[15];
3555 Double_t x1,x2,x3,y1,y2,y3,z1,z2,z3;
3556 track2.GetBoundaries(x3,x1);
3558 track2.GetFitPoint(x1,y1,z1);
3559 track2.GetFitPoint(x2,y2,z2);
3560 track2.GetFitPoint(x3,y3,z3);
3562 //is track pointing to the vertex ?
3565 polytrack.GetFitPoint(x0,y0,z0);
3578 x[4]=F1(x1,y1,x2,y2,x3,y3);
3580 // if (TMath::Abs(x[4]) >= cuts[0]) continue; //
3581 x[2]=F2(x1,y1,x2,y2,x3,y3);
3583 //if (TMath::Abs(x[4]*x1-x[2]) >= cuts[1]) continue;
3584 //x[3]=F3(x1,y1,x2,y2,z1,z2);
3585 x[3]=F3n(x1,y1,x3,y3,z1,z3,x[4]);
3586 //if (TMath::Abs(x[3]) > cuts[3]) continue;
3589 Double_t sy =0.1, sz =0.1;
3590 Double_t sy1=0.02, sz1=0.02;
3591 Double_t sy2=0.02, sz2=0.02;
3595 sy3=25000*x[4]*x[4]+0.1, sy=0.1, sz=0.1;
3598 Double_t f40=(F1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
3599 Double_t f42=(F1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
3600 Double_t f43=(F1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
3601 Double_t f20=(F2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
3602 Double_t f22=(F2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
3603 Double_t f23=(F2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
3605 Double_t f30=(F3(x1,y1+sy,x3,y3,z1,z3)-x[3])/sy;
3606 Double_t f31=(F3(x1,y1,x3,y3,z1+sz,z3)-x[3])/sz;
3607 Double_t f32=(F3(x1,y1,x3,y3+sy,z1,z3)-x[3])/sy;
3608 Double_t f34=(F3(x1,y1,x3,y3,z1,z3+sz)-x[3])/sz;
3613 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3614 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3615 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3616 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3617 c[13]=f30*sy1*f40+f32*sy2*f42;
3618 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3620 //Int_t row1 = fSectors->GetRowNumber(x1);
3621 Int_t row1 = GetRowNumber(x1);
3625 seed->~AliTPCseed();
3626 AliTPCseed *track=new(seed) AliTPCseed(x1,sec*alpha+shift,x,c,index);
3627 track->SetIsSeeding(kTRUE);
3628 Int_t rc=FollowProlongation(*track, i2);
3629 if (constrain) track->SetBConstrain(1);
3631 track->SetBConstrain(0);
3632 track->SetLastPoint(row1+fInnerSec->GetNRows()); // first cluster in track position
3633 track->SetFirstPoint(track->GetLastPoint());
3635 if (rc==0 || track->GetNumberOfClusters()<(i1-i2)*0.5 ||
3636 track->GetNumberOfClusters() < track->GetNFoundable()*0.6 ||
3637 track->GetNShared()>0.4*track->GetNumberOfClusters()) {
3640 seed->~AliTPCseed();
3643 arr->AddLast(track);
3644 seed = new AliTPCseed;
3648 } // if accepted seed
3651 Info("MakeSeeds2","\nSeeding statiistic:\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin3);
3657 AliTPCseed *AliTPCtrackerMI::MakeSeed(AliTPCseed *track, Float_t r0, Float_t r1, Float_t r2)
3661 //reseed using track points
3662 Int_t p0 = int(r0*track->GetNumberOfClusters()); // point 0
3663 Int_t p1 = int(r1*track->GetNumberOfClusters());
3664 Int_t p2 = int(r2*track->GetNumberOfClusters()); // last point
3666 Double_t x0[3],x1[3],x2[3];
3667 for (Int_t i=0;i<3;i++){
3673 // find track position at given ratio of the length
3674 Int_t sec0=0, sec1=0, sec2=0;
3677 for (Int_t i=0;i<160;i++){
3678 if (track->GetClusterPointer(i)){
3680 AliTPCTrackerPoint *trpoint =track->GetTrackPoint(i);
3681 if ( (index<p0) || x0[0]<0 ){
3682 if (trpoint->GetX()>1){
3683 clindex = track->GetClusterIndex2(i);
3685 x0[0] = trpoint->GetX();
3686 x0[1] = trpoint->GetY();
3687 x0[2] = trpoint->GetZ();
3688 sec0 = ((clindex&0xff000000)>>24)%18;
3693 if ( (index<p1) &&(trpoint->GetX()>1)){
3694 clindex = track->GetClusterIndex2(i);
3696 x1[0] = trpoint->GetX();
3697 x1[1] = trpoint->GetY();
3698 x1[2] = trpoint->GetZ();
3699 sec1 = ((clindex&0xff000000)>>24)%18;
3702 if ( (index<p2) &&(trpoint->GetX()>1)){
3703 clindex = track->GetClusterIndex2(i);
3705 x2[0] = trpoint->GetX();
3706 x2[1] = trpoint->GetY();
3707 x2[2] = trpoint->GetZ();
3708 sec2 = ((clindex&0xff000000)>>24)%18;
3715 Double_t alpha, cs,sn, xx2,yy2;
3717 alpha = (sec1-sec2)*fSectors->GetAlpha();
3718 cs = TMath::Cos(alpha);
3719 sn = TMath::Sin(alpha);
3720 xx2= x1[0]*cs-x1[1]*sn;
3721 yy2= x1[0]*sn+x1[1]*cs;
3725 alpha = (sec0-sec2)*fSectors->GetAlpha();
3726 cs = TMath::Cos(alpha);
3727 sn = TMath::Sin(alpha);
3728 xx2= x0[0]*cs-x0[1]*sn;
3729 yy2= x0[0]*sn+x0[1]*cs;
3735 Double_t x[5],c[15];
3739 x[4]=F1(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]);
3740 // if (x[4]>1) return 0;
3741 x[2]=F2(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]);
3742 x[3]=F3n(x2[0],x2[1],x0[0],x0[1],x2[2],x0[2],x[4]);
3743 //if (TMath::Abs(x[3]) > 2.2) return 0;
3744 //if (TMath::Abs(x[2]) > 1.99) return 0;
3746 Double_t sy =0.1, sz =0.1;
3748 Double_t sy1=0.02+track->GetSigmaY2(), sz1=0.02+track->GetSigmaZ2();
3749 Double_t sy2=0.01+track->GetSigmaY2(), sz2=0.01+track->GetSigmaZ2();
3750 Double_t sy3=0.01+track->GetSigmaY2();
3752 Double_t f40=(F1(x2[0],x2[1]+sy,x1[0],x1[1],x0[0],x0[1])-x[4])/sy;
3753 Double_t f42=(F1(x2[0],x2[1],x1[0],x1[1]+sy,x0[0],x0[1])-x[4])/sy;
3754 Double_t f43=(F1(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]+sy)-x[4])/sy;
3755 Double_t f20=(F2(x2[0],x2[1]+sy,x1[0],x1[1],x0[0],x0[1])-x[2])/sy;
3756 Double_t f22=(F2(x2[0],x2[1],x1[0],x1[1]+sy,x0[0],x0[1])-x[2])/sy;
3757 Double_t f23=(F2(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]+sy)-x[2])/sy;
3759 Double_t f30=(F3(x2[0],x2[1]+sy,x0[0],x0[1],x2[2],x0[2])-x[3])/sy;
3760 Double_t f31=(F3(x2[0],x2[1],x0[0],x0[1],x2[2]+sz,x0[2])-x[3])/sz;
3761 Double_t f32=(F3(x2[0],x2[1],x0[0],x0[1]+sy,x2[2],x0[2])-x[3])/sy;
3762 Double_t f34=(F3(x2[0],x2[1],x0[0],x0[1],x2[2],x0[2]+sz)-x[3])/sz;
3767 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3768 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3769 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3770 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3771 c[13]=f30*sy1*f40+f32*sy2*f42;
3772 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3774 // Int_t row1 = fSectors->GetRowNumber(x2[0]);
3775 AliTPCseed *seed=new AliTPCseed(x2[0], sec2*fSectors->GetAlpha()+fSectors->GetAlphaShift(), x, c, 0);
3776 // Double_t y0,z0,y1,z1, y2,z2;
3777 //seed->GetProlongation(x0[0],y0,z0);
3778 // seed->GetProlongation(x1[0],y1,z1);
3779 //seed->GetProlongation(x2[0],y2,z2);
3781 seed->SetLastPoint(pp2);
3782 seed->SetFirstPoint(pp2);
3789 AliTPCseed *AliTPCtrackerMI::ReSeed(AliTPCseed *track, Float_t r0, Float_t r1, Float_t r2)
3793 //reseed using founded clusters
3795 // Find the number of clusters
3796 Int_t nclusters = 0;
3797 for (Int_t irow=0;irow<160;irow++){
3798 if (track->GetClusterIndex(irow)>0) nclusters++;
3802 ipos[0] = TMath::Max(int(r0*nclusters),0); // point 0 cluster
3803 ipos[1] = TMath::Min(int(r1*nclusters),nclusters-1); //
3804 ipos[2] = TMath::Min(int(r2*nclusters),nclusters-1); // last point
3808 Int_t row[3],sec[3]={0,0,0};
3810 // find track row position at given ratio of the length
3812 for (Int_t irow=0;irow<160;irow++){
3813 if (track->GetClusterIndex2(irow)<0) continue;
3815 for (Int_t ipoint=0;ipoint<3;ipoint++){
3816 if (index<=ipos[ipoint]) row[ipoint] = irow;
3820 //Get cluster and sector position
3821 for (Int_t ipoint=0;ipoint<3;ipoint++){
3822 Int_t clindex = track->GetClusterIndex2(row[ipoint]);
3823 AliTPCclusterMI * cl = GetClusterMI(clindex);
3826 // AliTPCclusterMI * cl = GetClusterMI(clindex);
3829 sec[ipoint] = ((clindex&0xff000000)>>24)%18;
3830 xyz[ipoint][0] = GetXrow(row[ipoint]);
3831 xyz[ipoint][1] = cl->GetY();
3832 xyz[ipoint][2] = cl->GetZ();
3836 // Calculate seed state vector and covariance matrix
3838 Double_t alpha, cs,sn, xx2,yy2;
3840 alpha = (sec[1]-sec[2])*fSectors->GetAlpha();
3841 cs = TMath::Cos(alpha);
3842 sn = TMath::Sin(alpha);
3843 xx2= xyz[1][0]*cs-xyz[1][1]*sn;
3844 yy2= xyz[1][0]*sn+xyz[1][1]*cs;
3848 alpha = (sec[0]-sec[2])*fSectors->GetAlpha();
3849 cs = TMath::Cos(alpha);
3850 sn = TMath::Sin(alpha);
3851 xx2= xyz[0][0]*cs-xyz[0][1]*sn;
3852 yy2= xyz[0][0]*sn+xyz[0][1]*cs;
3858 Double_t x[5],c[15];
3862 x[4]=F1(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
3863 x[2]=F2(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
3864 x[3]=F3n(xyz[2][0],xyz[2][1],xyz[0][0],xyz[0][1],xyz[2][2],xyz[0][2],x[4]);
3866 Double_t sy =0.1, sz =0.1;
3868 Double_t sy1=0.2, sz1=0.2;
3869 Double_t sy2=0.2, sz2=0.2;
3872 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;
3873 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;
3874 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;
3875 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;
3876 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;
3877 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;
3879 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;
3880 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;
3881 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;
3882 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;
3887 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3888 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3889 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3890 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3891 c[13]=f30*sy1*f40+f32*sy2*f42;
3892 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3894 // Int_t row1 = fSectors->GetRowNumber(xyz[2][0]);
3895 AliTPCseed *seed=new AliTPCseed(xyz[2][0], sec[2]*fSectors->GetAlpha()+fSectors->GetAlphaShift(), x, c, 0);
3896 seed->SetLastPoint(row[2]);
3897 seed->SetFirstPoint(row[2]);
3902 AliTPCseed *AliTPCtrackerMI::ReSeed(AliTPCseed *track,Int_t r0, Bool_t forward)
3906 //reseed using founded clusters
3909 Int_t row[3]={0,0,0};
3910 Int_t sec[3]={0,0,0};
3912 // forward direction
3914 for (Int_t irow=r0;irow<160;irow++){
3915 if (track->GetClusterIndex(irow)>0){
3920 for (Int_t irow=160;irow>r0;irow--){
3921 if (track->GetClusterIndex(irow)>0){
3926 for (Int_t irow=row[2]-15;irow>row[0];irow--){
3927 if (track->GetClusterIndex(irow)>0){
3935 for (Int_t irow=0;irow<r0;irow++){
3936 if (track->GetClusterIndex(irow)>0){
3941 for (Int_t irow=r0;irow>0;irow--){
3942 if (track->GetClusterIndex(irow)>0){
3947 for (Int_t irow=row[2]-15;irow>row[0];irow--){
3948 if (track->GetClusterIndex(irow)>0){
3955 if ((row[2]-row[0])<20) return 0;
3956 if (row[1]==0) return 0;
3959 //Get cluster and sector position
3960 for (Int_t ipoint=0;ipoint<3;ipoint++){
3961 Int_t clindex = track->GetClusterIndex2(row[ipoint]);
3962 AliTPCclusterMI * cl = GetClusterMI(clindex);
3965 // AliTPCclusterMI * cl = GetClusterMI(clindex);
3968 sec[ipoint] = ((clindex&0xff000000)>>24)%18;
3969 xyz[ipoint][0] = GetXrow(row[ipoint]);
3970 AliTPCTrackerPoint * point = track->GetTrackPoint(row[ipoint]);
3971 if (point&&ipoint<2){
3973 xyz[ipoint][1] = point->GetY();
3974 xyz[ipoint][2] = point->GetZ();
3977 xyz[ipoint][1] = cl->GetY();
3978 xyz[ipoint][2] = cl->GetZ();
3985 // Calculate seed state vector and covariance matrix
3987 Double_t alpha, cs,sn, xx2,yy2;
3989 alpha = (sec[1]-sec[2])*fSectors->GetAlpha();
3990 cs = TMath::Cos(alpha);
3991 sn = TMath::Sin(alpha);
3992 xx2= xyz[1][0]*cs-xyz[1][1]*sn;
3993 yy2= xyz[1][0]*sn+xyz[1][1]*cs;
3997 alpha = (sec[0]-sec[2])*fSectors->GetAlpha();
3998 cs = TMath::Cos(alpha);
3999 sn = TMath::Sin(alpha);
4000 xx2= xyz[0][0]*cs-xyz[0][1]*sn;
4001 yy2= xyz[0][0]*sn+xyz[0][1]*cs;
4007 Double_t x[5],c[15];
4011 x[4]=F1(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
4012 x[2]=F2(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
4013 x[3]=F3n(xyz[2][0],xyz[2][1],xyz[0][0],xyz[0][1],xyz[2][2],xyz[0][2],x[4]);
4015 Double_t sy =0.1, sz =0.1;
4017 Double_t sy1=0.2, sz1=0.2;
4018 Double_t sy2=0.2, sz2=0.2;
4021 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;
4022 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;
4023 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;
4024 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;
4025 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;
4026 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;
4028 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;
4029 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;
4030 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;
4031 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;
4036 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
4037 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
4038 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
4039 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
4040 c[13]=f30*sy1*f40+f32*sy2*f42;
4041 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
4043 // Int_t row1 = fSectors->GetRowNumber(xyz[2][0]);
4044 AliTPCseed *seed=new AliTPCseed(xyz[2][0], sec[2]*fSectors->GetAlpha()+fSectors->GetAlphaShift(), x, c, 0);
4045 seed->SetLastPoint(row[2]);
4046 seed->SetFirstPoint(row[2]);
4047 for (Int_t i=row[0];i<row[2];i++){
4048 seed->SetClusterIndex(i, track->GetClusterIndex(i));
4056 void AliTPCtrackerMI::FindMultiMC(TObjArray * array, AliESDEvent */*esd*/, Int_t iter)
4059 // find multi tracks - THIS FUNCTION IS ONLY FOR DEBUG PURPOSES
4061 // Use AliTPCReconstructor::StreamLevel()>2 if you want to tune parameters - cuts
4063 // Two reasons to have multiple find tracks
4064 // 1. Curling tracks can be find more than once
4065 // 2. Splitted tracks
4066 // a.) Multiple seeding to increase tracking efficiency - (~ 100% reached)
4067 // b.) Edge effect on the sector boundaries
4070 // Algorithm done in 2 phases - because of CPU consumption
4071 // it is n^2 algorithm - for lead-lead 20000x20000 combination are investigated
4073 // Algorihm for curling tracks sign:
4074 // 1 phase -makes a very rough fast cuts to minimize combinatorics
4075 // a.) opposite sign
4076 // b.) one of the tracks - not pointing to the primary vertex -
4077 // c.) delta tan(theta)
4079 // 2 phase - calculates DCA between tracks - time consument
4084 // General cuts - for splitted tracks and for curling tracks
4086 const Float_t kMaxdPhi = 0.2; // maximal distance in phi
4088 // Curling tracks cuts
4093 Int_t nentries = array->GetEntriesFast();
4094 AliHelix *helixes = new AliHelix[nentries];
4095 Float_t *xm = new Float_t[nentries];
4096 Float_t *dz0 = new Float_t[nentries];
4097 Float_t *dz1 = new Float_t[nentries];
4103 // Find track COG in x direction - point with best defined parameters
4105 for (Int_t i=0;i<nentries;i++){
4106 AliTPCseed* track = (AliTPCseed*)array->At(i);
4107 if (!track) continue;
4108 track->SetCircular(0);
4109 new (&helixes[i]) AliHelix(*track);
4113 track->GetDZ(GetX(),GetY(),GetZ(),GetBz(),dz);
4116 for (Int_t icl=0; icl<160; icl++){
4117 AliTPCclusterMI * cl = track->GetClusterPointer(icl);
4123 if (ncl>0) xm[i]/=Float_t(ncl);
4125 TTreeSRedirector &cstream = *fDebugStreamer;
4127 for (Int_t i0=0;i0<nentries;i0++){
4128 AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
4129 if (!track0) continue;
4130 Float_t xc0 = helixes[i0].GetHelix(6);
4131 Float_t yc0 = helixes[i0].GetHelix(7);
4132 Float_t r0 = helixes[i0].GetHelix(8);
4133 Float_t rc0 = TMath::Sqrt(xc0*xc0+yc0*yc0);
4134 Float_t fi0 = TMath::ATan2(yc0,xc0);
4136 for (Int_t i1=i0+1;i1<nentries;i1++){
4137 AliTPCseed * track1 = (AliTPCseed*)array->At(i1);
4138 if (!track1) continue;
4139 Int_t lab0=track0->GetLabel();
4140 Int_t lab1=track1->GetLabel();
4141 if (TMath::Abs(lab0)!=TMath::Abs(lab1)) continue;
4143 Float_t xc1 = helixes[i1].GetHelix(6);
4144 Float_t yc1 = helixes[i1].GetHelix(7);
4145 Float_t r1 = helixes[i1].GetHelix(8);
4146 Float_t rc1 = TMath::Sqrt(xc1*xc1+yc1*yc1);
4147 Float_t fi1 = TMath::ATan2(yc1,xc1);
4149 Float_t dfi = fi0-fi1;
4152 if (dfi>1.5*TMath::Pi()) dfi-=TMath::Pi(); // take care about edge effect
4153 if (dfi<-1.5*TMath::Pi()) dfi+=TMath::Pi(); //
4154 if (TMath::Abs(dfi)>kMaxdPhi&&helixes[i0].GetHelix(4)*helixes[i1].GetHelix(4)<0){
4156 // if short tracks with undefined sign
4157 fi1 = -TMath::ATan2(yc1,-xc1);
4160 Float_t dtheta = TMath::Abs(track0->GetTgl()-track1->GetTgl())<TMath::Abs(track0->GetTgl()+track1->GetTgl())? track0->GetTgl()-track1->GetTgl():track0->GetTgl()+track1->GetTgl();
4163 // debug stream to tune "fast cuts"
4165 Double_t dist[3]; // distance at X
4166 Double_t mdist[3]={0,0,0}; // mean distance X+-40cm
4167 track0->GetDistance(track1,0.5*(xm[i0]+xm[i1])-40.,dist,AliTracker::GetBz());
4168 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[i]);
4169 track0->GetDistance(track1,0.5*(xm[i0]+xm[i1])+40.,dist,AliTracker::GetBz());
4170 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[i]);
4171 track0->GetDistance(track1,0.5*(xm[i0]+xm[i1]),dist,AliTracker::GetBz());
4172 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[i]);
4173 for (Int_t i=0;i<3;i++) mdist[i]*=0.33333;
4177 for (Int_t icl=0; icl<160; icl++){
4178 AliTPCclusterMI * cl0 = track0->GetClusterPointer(icl);
4179 AliTPCclusterMI * cl1 = track1->GetClusterPointer(icl);
4182 if (cl0==cl1) sums++;
4190 "Tr0.="<<track0<< // seed0
4191 "Tr1.="<<track1<< // seed1
4192 "h0.="<<&helixes[i0]<<
4193 "h1.="<<&helixes[i1]<<
4195 "sum="<<sum<< //the sum of rows with cl in both
4196 "sums="<<sums<< //the sum of shared clusters
4197 "xm0="<<xm[i0]<< // the center of track
4198 "xm1="<<xm[i1]<< // the x center of track
4199 // General cut variables
4200 "dfi="<<dfi<< // distance in fi angle
4201 "dtheta="<<dtheta<< // distance int theta angle
4207 "dist0="<<dist[0]<< //distance x
4208 "dist1="<<dist[1]<< //distance y
4209 "dist2="<<dist[2]<< //distance z
4210 "mdist0="<<mdist[0]<< //distance x
4211 "mdist1="<<mdist[1]<< //distance y
4212 "mdist2="<<mdist[2]<< //distance z
4225 if (AliTPCReconstructor::StreamLevel()>1) {
4226 AliInfo("Time for curling tracks removal DEBUGGING MC");
4232 void AliTPCtrackerMI::FindSplitted(TObjArray * array, AliESDEvent */*esd*/, Int_t iter)
4236 // Two reasons to have multiple find tracks
4237 // 1. Curling tracks can be find more than once
4238 // 2. Splitted tracks
4239 // a.) Multiple seeding to increase tracking efficiency - (~ 100% reached)
4240 // b.) Edge effect on the sector boundaries
4242 // This function tries to find tracks closed in the parametric space
4244 // cut logic if distance is bigger than cut continue - Do Nothing
4245 const Float_t kMaxdTheta = 0.05; // maximal distance in theta
4246 const Float_t kMaxdPhi = 0.05; // maximal deistance in phi
4247 const Float_t kdelta = 40.; //delta r to calculate track distance
4249 // const Float_t kMaxDist0 = 1.; // maximal distance 0
4250 //const Float_t kMaxDist1 = 0.3; // maximal distance 1 - cut if track in separate rows
4253 TCut csec("csec","abs(Tr0.fRelativeSector-Tr1.fRelativeSector)<2");
4254 TCut cdtheta("cdtheta","abs(dtheta)<0.05");
4259 Int_t nentries = array->GetEntriesFast();
4260 AliHelix *helixes = new AliHelix[nentries];
4261 Float_t *xm = new Float_t[nentries];
4267 //Sort tracks according quality
4269 Int_t nseed = array->GetEntriesFast();
4270 Float_t * quality = new Float_t[nseed];
4271 Int_t * indexes = new Int_t[nseed];
4272 for (Int_t i=0; i<nseed; i++) {
4273 AliTPCseed *pt=(AliTPCseed*)array->UncheckedAt(i);
4278 pt->UpdatePoints(); //select first last max dens points
4279 Float_t * points = pt->GetPoints();
4280 if (points[3]<0.8) quality[i] =-1;
4281 quality[i] = (points[2]-points[0])+pt->GetNumberOfClusters();
4282 //prefer high momenta tracks if overlaps
4283 quality[i] *= TMath::Sqrt(TMath::Abs(pt->Pt())+0.5);
4285 TMath::Sort(nseed,quality,indexes);
4289 // Find track COG in x direction - point with best defined parameters
4291 for (Int_t i=0;i<nentries;i++){
4292 AliTPCseed* track = (AliTPCseed*)array->At(i);
4293 if (!track) continue;
4294 track->SetCircular(0);
4295 new (&helixes[i]) AliHelix(*track);
4298 for (Int_t icl=0; icl<160; icl++){
4299 AliTPCclusterMI * cl = track->GetClusterPointer(icl);
4305 if (ncl>0) xm[i]/=Float_t(ncl);
4307 TTreeSRedirector &cstream = *fDebugStreamer;
4309 for (Int_t is0=0;is0<nentries;is0++){
4310 Int_t i0 = indexes[is0];
4311 AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
4312 if (!track0) continue;
4313 if (track0->GetKinkIndexes()[0]!=0) continue;
4314 Float_t xc0 = helixes[i0].GetHelix(6);
4315 Float_t yc0 = helixes[i0].GetHelix(7);
4316 Float_t fi0 = TMath::ATan2(yc0,xc0);
4318 for (Int_t is1=is0+1;is1<nentries;is1++){
4319 Int_t i1 = indexes[is1];
4320 AliTPCseed * track1 = (AliTPCseed*)array->At(i1);
4321 if (!track1) continue;
4323 if (TMath::Abs(track0->GetRelativeSector()-track1->GetRelativeSector())>1) continue;
4324 if (track1->GetKinkIndexes()[0]>0 &&track0->GetKinkIndexes()[0]<0) continue;
4325 if (track1->GetKinkIndexes()[0]!=0) continue;
4327 Float_t dtheta = TMath::Abs(track0->GetTgl()-track1->GetTgl())<TMath::Abs(track0->GetTgl()+track1->GetTgl())? track0->GetTgl()-track1->GetTgl():track0->GetTgl()+track1->GetTgl();
4328 if (TMath::Abs(dtheta)>kMaxdTheta) continue;
4330 Float_t xc1 = helixes[i1].GetHelix(6);
4331 Float_t yc1 = helixes[i1].GetHelix(7);
4332 Float_t fi1 = TMath::ATan2(yc1,xc1);
4334 Float_t dfi = fi0-fi1;
4335 if (dfi>1.5*TMath::Pi()) dfi-=TMath::Pi(); // take care about edge effect
4336 if (dfi<-1.5*TMath::Pi()) dfi+=TMath::Pi(); //
4337 if (TMath::Abs(dfi)>kMaxdPhi&&helixes[i0].GetHelix(4)*helixes[i1].GetHelix(4)<0){
4339 // if short tracks with undefined sign
4340 fi1 = -TMath::ATan2(yc1,-xc1);
4343 if (TMath::Abs(dfi)>kMaxdPhi) continue;
4350 for (Int_t icl=0; icl<160; icl++){
4351 Int_t index0=track0->GetClusterIndex2(icl);
4352 Int_t index1=track1->GetClusterIndex2(icl);
4353 Bool_t used0 = (index0>0 && !(index0&0x8000));
4354 Bool_t used1 = (index1>0 && !(index1&0x8000));
4356 if (used0) sum0++; // used cluster0
4357 if (used1) sum1++; // used clusters1
4358 if (used0&&used1) sum++;
4359 if (index0==index1 && used0 && used1) sums++;
4363 if (sums<10) continue;
4364 if (sum<40) continue;
4365 if (sums/Float_t(TMath::Min(sum0,sum1))<0.5) continue;
4367 Double_t dist[5][4]; // distance at X
4368 Double_t mdist[4]={0,0,0,0}; // mean distance on range +-delta
4372 track0->GetDistance(track1,xm[i0],dist[0],AliTracker::GetBz());
4373 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[0][i]);
4374 track0->GetDistance(track1,xm[i1],dist[1],AliTracker::GetBz());
4375 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[1][i]);
4377 track0->GetDistance(track1,TMath::Min(xm[i1],xm[i0])-kdelta,dist[2],AliTracker::GetBz());
4378 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[2][i]);
4379 track0->GetDistance(track1,TMath::Max(xm[i1],xm[i0])+kdelta,dist[3],AliTracker::GetBz());
4380 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[3][i]);
4382 track0->GetDistance(track1,(xm[i1]+xm[i0])*0.5,dist[4],AliTracker::GetBz());
4383 for (Int_t i=0;i<3;i++) mdist[i]+=TMath::Abs(dist[4][i]);
4384 for (Int_t i=0;i<3;i++) mdist[i]*=0.2;
4387 Int_t lab0=track0->GetLabel();
4388 Int_t lab1=track1->GetLabel();
4389 cstream<<"Splitted"<<
4393 "Tr0.="<<track0<< // seed0
4394 "Tr1.="<<track1<< // seed1
4395 "h0.="<<&helixes[i0]<<
4396 "h1.="<<&helixes[i1]<<
4398 "sum="<<sum<< //the sum of rows with cl in both
4399 "sum0="<<sum0<< //the sum of rows with cl in 0 track
4400 "sum1="<<sum1<< //the sum of rows with cl in 1 track
4401 "sums="<<sums<< //the sum of shared clusters
4402 "xm0="<<xm[i0]<< // the center of track
4403 "xm1="<<xm[i1]<< // the x center of track
4404 // General cut variables
4405 "dfi="<<dfi<< // distance in fi angle
4406 "dtheta="<<dtheta<< // distance int theta angle
4409 "dist0="<<dist[4][0]<< //distance x
4410 "dist1="<<dist[4][1]<< //distance y
4411 "dist2="<<dist[4][2]<< //distance z
4412 "mdist0="<<mdist[0]<< //distance x
4413 "mdist1="<<mdist[1]<< //distance y
4414 "mdist2="<<mdist[2]<< //distance z
4417 delete array->RemoveAt(i1);
4422 AliInfo("Time for splitted tracks removal");
4428 void AliTPCtrackerMI::FindCurling(TObjArray * array, AliESDEvent *esd, Int_t iter)
4431 // find Curling tracks
4432 // Use AliTPCReconstructor::StreamLevel()>1 if you want to tune parameters - cuts
4435 // Algorithm done in 2 phases - because of CPU consumption
4436 // it is n^2 algorithm - for lead-lead 20000x20000 combination are investigated
4437 // see detal in MC part what can be used to cut
4441 const Float_t kMaxC = 400; // maximal curvature to of the track
4442 const Float_t kMaxdTheta = 0.15; // maximal distance in theta
4443 const Float_t kMaxdPhi = 0.15; // maximal distance in phi
4444 const Float_t kPtRatio = 0.3; // ratio between pt
4445 const Float_t kMinDCAR = 2.; // distance to the primary vertex in r - see cpipe cut
4448 // Curling tracks cuts
4451 const Float_t kMaxDeltaRMax = 40; // distance in outer radius
4452 const Float_t kMaxDeltaRMin = 5.; // distance in lower radius - see cpipe cut
4453 const Float_t kMinAngle = 2.9; // angle between tracks
4454 const Float_t kMaxDist = 5; // biggest distance
4456 // The cuts can be tuned using the "MC information stored in Multi tree ==> see FindMultiMC
4459 TCut csign("csign","Tr0.fP[4]*Tr1.fP[4]<0"); //opposite sign
4460 TCut cmax("cmax","abs(Tr0.GetC())>1/400");
4461 TCut cda("cda","sqrt(dtheta^2+dfi^2)<0.15");
4462 TCut ccratio("ccratio","abs((Tr0.fP[4]+Tr1.fP[4])/(abs(Tr0.fP[4])+abs(Tr1.fP[4])))<0.3");
4463 TCut cpipe("cpipe", "min(abs(r0-rc0),abs(r1-rc1))>5");
4465 TCut cdrmax("cdrmax","abs(abs(rc0+r0)-abs(rc1+r1))<40")
4466 TCut cdrmin("cdrmin","abs(abs(rc0+r0)-abs(rc1+r1))<10")
4468 Multi->Draw("dfi","iter==0"+csign+cmax+cda+ccratio); ~94% of curling tracks fulfill
4469 Multi->Draw("min(abs(r0-rc0),abs(r1-rc1))","iter==0&&abs(lab1)==abs(lab0)"+csign+cmax+cda+ccratio+cpipe+cdrmin+cdrmax); //80%
4471 Curling2->Draw("dfi","iter==0&&abs(lab0)==abs(lab1)"+csign+cmax+cdtheta+cdfi+ccratio)
4477 Int_t nentries = array->GetEntriesFast();
4478 AliHelix *helixes = new AliHelix[nentries];
4479 for (Int_t i=0;i<nentries;i++){
4480 AliTPCseed* track = (AliTPCseed*)array->At(i);
4481 if (!track) continue;
4482 track->SetCircular(0);
4483 new (&helixes[i]) AliHelix(*track);
4489 Double_t phase[2][2],radius[2];
4493 TTreeSRedirector &cstream = *fDebugStreamer;
4495 for (Int_t i0=0;i0<nentries;i0++){
4496 AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
4497 if (!track0) continue;
4498 if (TMath::Abs(track0->GetC())<1/kMaxC) continue;
4499 Float_t xc0 = helixes[i0].GetHelix(6);
4500 Float_t yc0 = helixes[i0].GetHelix(7);
4501 Float_t r0 = helixes[i0].GetHelix(8);
4502 Float_t rc0 = TMath::Sqrt(xc0*xc0+yc0*yc0);
4503 Float_t fi0 = TMath::ATan2(yc0,xc0);
4505 for (Int_t i1=i0+1;i1<nentries;i1++){
4506 AliTPCseed * track1 = (AliTPCseed*)array->At(i1);
4507 if (!track1) continue;
4508 if (TMath::Abs(track1->GetC())<1/kMaxC) continue;
4509 Float_t xc1 = helixes[i1].GetHelix(6);
4510 Float_t yc1 = helixes[i1].GetHelix(7);
4511 Float_t r1 = helixes[i1].GetHelix(8);
4512 Float_t rc1 = TMath::Sqrt(xc1*xc1+yc1*yc1);
4513 Float_t fi1 = TMath::ATan2(yc1,xc1);
4515 Float_t dfi = fi0-fi1;
4518 if (dfi>1.5*TMath::Pi()) dfi-=TMath::Pi(); // take care about edge effect
4519 if (dfi<-1.5*TMath::Pi()) dfi+=TMath::Pi(); //
4520 Float_t dtheta = TMath::Abs(track0->GetTgl()-track1->GetTgl())<TMath::Abs(track0->GetTgl()+track1->GetTgl())? track0->GetTgl()-track1->GetTgl():track0->GetTgl()+track1->GetTgl();
4524 if (track0->GetBConstrain()&&track1->GetBConstrain()) continue; // not constrained
4525 if (track1->GetSigned1Pt()*track0->GetSigned1Pt()>0) continue; // not the same sign
4526 if ( TMath::Abs(track1->GetTgl()+track0->GetTgl())>kMaxdTheta) continue; //distance in the Theta
4527 if ( TMath::Abs(dfi)>kMaxdPhi) continue; //distance in phi
4528 if ( TMath::Sqrt(dfi*dfi+dtheta*dtheta)>kMaxdPhi) continue; //common angular offset
4530 Float_t pt0 = track0->GetSignedPt();
4531 Float_t pt1 = track1->GetSignedPt();
4532 if ((TMath::Abs(pt0+pt1)/(TMath::Abs(pt0)+TMath::Abs(pt1)))>kPtRatio) continue;
4533 if ((iter==1) && TMath::Abs(TMath::Abs(rc0+r0)-TMath::Abs(rc1+r1))>kMaxDeltaRMax) continue;
4534 if ((iter!=1) &&TMath::Abs(TMath::Abs(rc0-r0)-TMath::Abs(rc1-r1))>kMaxDeltaRMin) continue;
4535 if (TMath::Min(TMath::Abs(rc0-r0),TMath::Abs(rc1-r1))<kMinDCAR) continue;
4538 // Now find closest approach
4542 Int_t npoints = helixes[i0].GetRPHIintersections(helixes[i1], phase, radius,10);
4543 if (npoints==0) continue;
4544 helixes[i0].GetClosestPhases(helixes[i1], phase);
4548 Double_t hangles[3];
4549 helixes[i0].Evaluate(phase[0][0],xyz0);
4550 helixes[i1].Evaluate(phase[0][1],xyz1);
4552 helixes[i0].GetAngle(phase[0][0],helixes[i1],phase[0][1],hangles);
4553 Double_t deltah[2],deltabest;
4554 if (TMath::Abs(hangles[2])<kMinAngle) continue;
4558 helixes[i0].ParabolicDCA(helixes[i1],phase[0][0],phase[0][1],radius[0],deltah[0],2);
4560 helixes[i0].ParabolicDCA(helixes[i1],phase[1][0],phase[1][1],radius[1],deltah[1],2);
4561 if (deltah[1]<deltah[0]) ibest=1;
4563 deltabest = TMath::Sqrt(deltah[ibest]);
4564 helixes[i0].Evaluate(phase[ibest][0],xyz0);
4565 helixes[i1].Evaluate(phase[ibest][1],xyz1);
4566 helixes[i0].GetAngle(phase[ibest][0],helixes[i1],phase[ibest][1],hangles);
4567 Double_t radiusbest = TMath::Sqrt(radius[ibest]);
4569 if (deltabest>kMaxDist) continue;
4570 // if (mindcar+mindcaz<40 && (TMath::Abs(hangles[2])<kMinAngle ||deltabest>3)) continue;
4571 Bool_t sign =kFALSE;
4572 if (hangles[2]>kMinAngle) sign =kTRUE;
4575 // circular[i0] = kTRUE;
4576 // circular[i1] = kTRUE;
4577 if (track0->OneOverPt()<track1->OneOverPt()){
4578 track0->SetCircular(track0->GetCircular()+1);
4579 track1->SetCircular(track1->GetCircular()+2);
4582 track1->SetCircular(track1->GetCircular()+1);
4583 track0->SetCircular(track0->GetCircular()+2);
4586 if (AliTPCReconstructor::StreamLevel()>1){
4588 //debug stream to tune "fine" cuts
4589 Int_t lab0=track0->GetLabel();
4590 Int_t lab1=track1->GetLabel();
4591 cstream<<"Curling2"<<
4607 "npoints="<<npoints<<
4608 "hangles0="<<hangles[0]<<
4609 "hangles1="<<hangles[1]<<
4610 "hangles2="<<hangles[2]<<
4613 "radius="<<radiusbest<<
4614 "deltabest="<<deltabest<<
4615 "phase0="<<phase[ibest][0]<<
4616 "phase1="<<phase[ibest][1]<<
4624 if (AliTPCReconstructor::StreamLevel()>1) {
4625 AliInfo("Time for curling tracks removal");
4634 void AliTPCtrackerMI::FindKinks(TObjArray * array, AliESDEvent *esd)
4641 TObjArray *kinks= new TObjArray(10000);
4642 // TObjArray *v0s= new TObjArray(10000);
4643 Int_t nentries = array->GetEntriesFast();
4644 AliHelix *helixes = new AliHelix[nentries];
4645 Int_t *sign = new Int_t[nentries];
4646 Int_t *nclusters = new Int_t[nentries];
4647 Float_t *alpha = new Float_t[nentries];
4648 AliKink *kink = new AliKink();
4649 Int_t * usage = new Int_t[nentries];
4650 Float_t *zm = new Float_t[nentries];
4651 Float_t *z0 = new Float_t[nentries];
4652 Float_t *fim = new Float_t[nentries];
4653 Float_t *shared = new Float_t[nentries];
4654 Bool_t *circular = new Bool_t[nentries];
4655 Float_t *dca = new Float_t[nentries];
4656 //const AliESDVertex * primvertex = esd->GetVertex();
4658 // nentries = array->GetEntriesFast();
4663 for (Int_t i=0;i<nentries;i++){
4666 AliTPCseed* track = (AliTPCseed*)array->At(i);
4667 if (!track) continue;
4668 track->SetCircular(0);
4670 track->UpdatePoints();
4671 if (( track->GetPoints()[2]- track->GetPoints()[0])>5 && track->GetPoints()[3]>0.8){
4673 nclusters[i]=track->GetNumberOfClusters();
4674 alpha[i] = track->GetAlpha();
4675 new (&helixes[i]) AliHelix(*track);
4677 helixes[i].Evaluate(0,xyz);
4678 sign[i] = (track->GetC()>0) ? -1:1;
4681 if (track->GetProlongation(x,y,z)){
4683 fim[i] = alpha[i]+TMath::ATan2(y,x);
4686 zm[i] = track->GetZ();
4690 circular[i]= kFALSE;
4691 if (track->GetProlongation(0,y,z)) z0[i] = z;
4692 dca[i] = track->GetD(0,0);
4698 Int_t ncandidates =0;
4701 Double_t phase[2][2],radius[2];
4704 // Find circling track
4705 TTreeSRedirector &cstream = *fDebugStreamer;
4707 for (Int_t i0=0;i0<nentries;i0++){
4708 AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
4709 if (!track0) continue;
4710 if (track0->GetNumberOfClusters()<40) continue;
4711 if (TMath::Abs(1./track0->GetC())>200) continue;
4712 for (Int_t i1=i0+1;i1<nentries;i1++){
4713 AliTPCseed * track1 = (AliTPCseed*)array->At(i1);
4714 if (!track1) continue;
4715 if (track1->GetNumberOfClusters()<40) continue;
4716 if ( TMath::Abs(track1->GetTgl()+track0->GetTgl())>0.1) continue;
4717 if (track0->GetBConstrain()&&track1->GetBConstrain()) continue;
4718 if (TMath::Abs(1./track1->GetC())>200) continue;
4719 if (track1->GetSigned1Pt()*track0->GetSigned1Pt()>0) continue;
4720 if (track1->GetTgl()*track0->GetTgl()>0) continue;
4721 if (TMath::Max(TMath::Abs(1./track0->GetC()),TMath::Abs(1./track1->GetC()))>190) continue;
4722 if (track0->GetBConstrain()&&track1->OneOverPt()<track0->OneOverPt()) continue; //returning - lower momenta
4723 if (track1->GetBConstrain()&&track0->OneOverPt()<track1->OneOverPt()) continue; //returning - lower momenta
4725 Float_t mindcar = TMath::Min(TMath::Abs(dca[i0]),TMath::Abs(dca[i1]));
4726 if (mindcar<5) continue;
4727 Float_t mindcaz = TMath::Min(TMath::Abs(z0[i0]-GetZ()),TMath::Abs(z0[i1]-GetZ()));
4728 if (mindcaz<5) continue;
4729 if (mindcar+mindcaz<20) continue;
4732 Float_t xc0 = helixes[i0].GetHelix(6);
4733 Float_t yc0 = helixes[i0].GetHelix(7);
4734 Float_t r0 = helixes[i0].GetHelix(8);
4735 Float_t xc1 = helixes[i1].GetHelix(6);
4736 Float_t yc1 = helixes[i1].GetHelix(7);
4737 Float_t r1 = helixes[i1].GetHelix(8);
4739 Float_t rmean = (r0+r1)*0.5;
4740 Float_t delta =TMath::Sqrt((xc1-xc0)*(xc1-xc0)+(yc1-yc0)*(yc1-yc0));
4741 //if (delta>30) continue;
4742 if (delta>rmean*0.25) continue;
4743 if (TMath::Abs(r0-r1)/rmean>0.3) continue;
4745 Int_t npoints = helixes[i0].GetRPHIintersections(helixes[i1], phase, radius,10);
4746 if (npoints==0) continue;
4747 helixes[i0].GetClosestPhases(helixes[i1], phase);
4751 Double_t hangles[3];
4752 helixes[i0].Evaluate(phase[0][0],xyz0);
4753 helixes[i1].Evaluate(phase[0][1],xyz1);
4755 helixes[i0].GetAngle(phase[0][0],helixes[i1],phase[0][1],hangles);
4756 Double_t deltah[2],deltabest;
4757 if (hangles[2]<2.8) continue;
4760 helixes[i0].ParabolicDCA(helixes[i1],phase[0][0],phase[0][1],radius[0],deltah[0],2);
4762 helixes[i0].ParabolicDCA(helixes[i1],phase[1][0],phase[1][1],radius[1],deltah[1],2);
4763 if (deltah[1]<deltah[0]) ibest=1;
4765 deltabest = TMath::Sqrt(deltah[ibest]);
4766 helixes[i0].Evaluate(phase[ibest][0],xyz0);
4767 helixes[i1].Evaluate(phase[ibest][1],xyz1);
4768 helixes[i0].GetAngle(phase[ibest][0],helixes[i1],phase[ibest][1],hangles);
4769 Double_t radiusbest = TMath::Sqrt(radius[ibest]);
4771 if (deltabest>6) continue;
4772 if (mindcar+mindcaz<40 && (hangles[2]<3.12||deltabest>3)) continue;
4773 Bool_t sign =kFALSE;
4774 if (hangles[2]>3.06) sign =kTRUE;
4777 circular[i0] = kTRUE;
4778 circular[i1] = kTRUE;
4779 if (track0->OneOverPt()<track1->OneOverPt()){
4780 track0->SetCircular(track0->GetCircular()+1);
4781 track1->SetCircular(track1->GetCircular()+2);
4784 track1->SetCircular(track1->GetCircular()+1);
4785 track0->SetCircular(track0->GetCircular()+2);
4788 if (sign&&AliTPCReconstructor::StreamLevel()>1){
4790 Int_t lab0=track0->GetLabel();
4791 Int_t lab1=track1->GetLabel();
4792 cstream<<"Curling"<<
4799 "mindcar="<<mindcar<<
4800 "mindcaz="<<mindcaz<<
4803 "npoints="<<npoints<<
4804 "hangles0="<<hangles[0]<<
4805 "hangles2="<<hangles[2]<<
4810 "radius="<<radiusbest<<
4811 "deltabest="<<deltabest<<
4812 "phase0="<<phase[ibest][0]<<
4813 "phase1="<<phase[ibest][1]<<
4823 for (Int_t i =0;i<nentries;i++){
4824 if (sign[i]==0) continue;
4825 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
4832 Double_t cradius0 = 40*40;
4833 Double_t cradius1 = 270*270;
4836 Double_t cdist3=0.55;
4837 for (Int_t j =i+1;j<nentries;j++){
4839 if (sign[j]*sign[i]<1) continue;
4840 if ( (nclusters[i]+nclusters[j])>200) continue;
4841 if ( (nclusters[i]+nclusters[j])<80) continue;
4842 if ( TMath::Abs(zm[i]-zm[j])>60.) continue;
4843 if ( TMath::Abs(fim[i]-fim[j])>0.6 && TMath::Abs(fim[i]-fim[j])<5.7 ) continue;
4844 //AliTPCseed * track1 = (AliTPCseed*)array->At(j); Double_t phase[2][2],radius[2];
4845 Int_t npoints = helixes[i].GetRPHIintersections(helixes[j], phase, radius,20);
4846 if (npoints<1) continue;
4849 if (radius[0]<cradius0||radius[0]>cradius1) continue;
4852 if ( (radius[0]<cradius0||radius[0]>cradius1) && (radius[1]<cradius0||radius[1]>cradius1) ) continue;
4855 Double_t delta1=10000,delta2=10000;
4856 // cuts on the intersection radius
4857 helixes[i].LinearDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta1);
4858 if (radius[0]<20&&delta1<1) continue; //intersection at vertex
4859 if (radius[0]<10&&delta1<3) continue; //intersection at vertex
4861 helixes[i].LinearDCA(helixes[j],phase[1][0],phase[1][1],radius[1],delta2);
4862 if (radius[1]<20&&delta2<1) continue; //intersection at vertex
4863 if (radius[1]<10&&delta2<3) continue; //intersection at vertex
4866 Double_t distance1 = TMath::Min(delta1,delta2);
4867 if (distance1>cdist1) continue; // cut on DCA linear approximation
4869 npoints = helixes[i].GetRPHIintersections(helixes[j], phase, radius,20);
4870 helixes[i].ParabolicDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta1);
4871 if (radius[0]<20&&delta1<1) continue; //intersection at vertex
4872 if (radius[0]<10&&delta1<3) continue; //intersection at vertex
4875 helixes[i].ParabolicDCA(helixes[j],phase[1][0],phase[1][1],radius[1],delta2);
4876 if (radius[1]<20&&delta2<1) continue; //intersection at vertex
4877 if (radius[1]<10&&delta2<3) continue; //intersection at vertex
4879 distance1 = TMath::Min(delta1,delta2);
4882 rkink = TMath::Sqrt(radius[0]);
4885 rkink = TMath::Sqrt(radius[1]);
4887 if (distance1>cdist2) continue;
4890 AliTPCseed * track1 = (AliTPCseed*)array->At(j);
4893 Int_t row0 = GetRowNumber(rkink);
4894 if (row0<10) continue;
4895 if (row0>150) continue;
4898 Float_t dens00=-1,dens01=-1;
4899 Float_t dens10=-1,dens11=-1;
4901 Int_t found,foundable,shared;
4902 track0->GetClusterStatistic(0,row0-5, found, foundable,shared,kFALSE);
4903 if (foundable>5) dens00 = Float_t(found)/Float_t(foundable);
4904 track0->GetClusterStatistic(row0+5,155, found, foundable,shared,kFALSE);
4905 if (foundable>5) dens01 = Float_t(found)/Float_t(foundable);
4907 track1->GetClusterStatistic(0,row0-5, found, foundable,shared,kFALSE);
4908 if (foundable>10) dens10 = Float_t(found)/Float_t(foundable);
4909 track1->GetClusterStatistic(row0+5,155, found, foundable,shared,kFALSE);
4910 if (foundable>10) dens11 = Float_t(found)/Float_t(foundable);
4912 if (dens00<dens10 && dens01<dens11) continue;
4913 if (dens00>dens10 && dens01>dens11) continue;
4914 if (TMath::Max(dens00,dens10)<0.1) continue;
4915 if (TMath::Max(dens01,dens11)<0.3) continue;
4917 if (TMath::Min(dens00,dens10)>0.6) continue;
4918 if (TMath::Min(dens01,dens11)>0.6) continue;
4921 AliTPCseed * ktrack0, *ktrack1;
4930 if (TMath::Abs(ktrack0->GetC())>5) continue; // cut on the curvature for mother particle
4931 AliExternalTrackParam paramm(*ktrack0);
4932 AliExternalTrackParam paramd(*ktrack1);
4933 if (row0>60&&ktrack1->GetReference().GetX()>90.)new (¶md) AliExternalTrackParam(ktrack1->GetReference());
4936 kink->SetMother(paramm);
4937 kink->SetDaughter(paramd);
4940 Float_t x[3] = { kink->GetPosition()[0],kink->GetPosition()[1],kink->GetPosition()[2]};
4942 fParam->Transform0to1(x,index);
4943 fParam->Transform1to2(x,index);
4944 row0 = GetRowNumber(x[0]);
4946 if (kink->GetR()<100) continue;
4947 if (kink->GetR()>240) continue;
4948 if (kink->GetPosition()[2]/kink->GetR()>AliTPCReconstructor::GetCtgRange()) continue; //out of fiducial volume
4949 if (kink->GetDistance()>cdist3) continue;
4950 Float_t dird = kink->GetDaughterP()[0]*kink->GetPosition()[0]+kink->GetDaughterP()[1]*kink->GetPosition()[1]; // rough direction estimate
4951 if (dird<0) continue;
4953 Float_t dirm = kink->GetMotherP()[0]*kink->GetPosition()[0]+kink->GetMotherP()[1]*kink->GetPosition()[1]; // rough direction estimate
4954 if (dirm<0) continue;
4955 Float_t mpt = TMath::Sqrt(kink->GetMotherP()[0]*kink->GetMotherP()[0]+kink->GetMotherP()[1]*kink->GetMotherP()[1]);
4956 if (mpt<0.2) continue;
4959 //for high momenta momentum not defined well in first iteration
4960 Double_t qt = TMath::Sin(kink->GetAngle(2))*ktrack1->GetP();
4961 if (qt>0.35) continue;
4964 kink->SetLabel(CookLabel(ktrack0,0.4,0,row0),0);
4965 kink->SetLabel(CookLabel(ktrack1,0.4,row0,160),1);
4967 kink->SetTPCDensity(dens00,0,0);
4968 kink->SetTPCDensity(dens01,0,1);
4969 kink->SetTPCDensity(dens10,1,0);
4970 kink->SetTPCDensity(dens11,1,1);
4971 kink->SetIndex(i,0);
4972 kink->SetIndex(j,1);
4975 kink->SetTPCDensity(dens10,0,0);
4976 kink->SetTPCDensity(dens11,0,1);
4977 kink->SetTPCDensity(dens00,1,0);
4978 kink->SetTPCDensity(dens01,1,1);
4979 kink->SetIndex(j,0);
4980 kink->SetIndex(i,1);
4983 if (mpt<1||kink->GetAngle(2)>0.1){
4984 // angle and densities not defined yet
4985 if (kink->GetTPCDensityFactor()<0.8) continue;
4986 if ((2-kink->GetTPCDensityFactor())*kink->GetDistance() >0.25) continue;
4987 if (kink->GetAngle(2)*ktrack0->GetP()<0.003) continue; //too small angle
4988 if (kink->GetAngle(2)>0.2&&kink->GetTPCDensityFactor()<1.15) continue;
4989 if (kink->GetAngle(2)>0.2&&kink->GetTPCDensity(0,1)>0.05) continue;
4991 Float_t criticalangle = track0->GetSigmaSnp2()+track0->GetSigmaTgl2();
4992 criticalangle+= track1->GetSigmaSnp2()+track1->GetSigmaTgl2();
4993 criticalangle= 3*TMath::Sqrt(criticalangle);
4994 if (criticalangle>0.02) criticalangle=0.02;
4995 if (kink->GetAngle(2)<criticalangle) continue;
4998 Int_t drow = Int_t(2.+0.5/(0.05+kink->GetAngle(2))); // overlap region defined
4999 Float_t shapesum =0;
5001 for ( Int_t row = row0-drow; row<row0+drow;row++){
5002 if (row<0) continue;
5003 if (row>155) continue;
5004 if (ktrack0->GetClusterPointer(row)){
5005 AliTPCTrackerPoint *point =ktrack0->GetTrackPoint(row);
5006 shapesum+=point->GetSigmaY()+point->GetSigmaZ();
5009 if (ktrack1->GetClusterPointer(row)){
5010 AliTPCTrackerPoint *point =ktrack1->GetTrackPoint(row);
5011 shapesum+=point->GetSigmaY()+point->GetSigmaZ();
5016 kink->SetShapeFactor(-1.);
5019 kink->SetShapeFactor(shapesum/sum);
5021 // esd->AddKink(kink);
5022 kinks->AddLast(kink);
5028 // sort the kinks according quality - and refit them towards vertex
5030 Int_t nkinks = kinks->GetEntriesFast();
5031 Float_t *quality = new Float_t[nkinks];
5032 Int_t *indexes = new Int_t[nkinks];
5033 AliTPCseed *mothers = new AliTPCseed[nkinks];
5034 AliTPCseed *daughters = new AliTPCseed[nkinks];
5037 for (Int_t i=0;i<nkinks;i++){
5039 AliKink *kink = (AliKink*)kinks->At(i);
5041 // refit kinks towards vertex
5043 Int_t index0 = kink->GetIndex(0);
5044 Int_t index1 = kink->GetIndex(1);
5045 AliTPCseed * ktrack0 = (AliTPCseed*)array->At(index0);
5046 AliTPCseed * ktrack1 = (AliTPCseed*)array->At(index1);
5048 Int_t sumn=ktrack0->GetNumberOfClusters()+ktrack1->GetNumberOfClusters();
5050 // Refit Kink under if too small angle
5052 if (kink->GetAngle(2)<0.05){
5053 kink->SetTPCRow0(GetRowNumber(kink->GetR()));
5054 Int_t row0 = kink->GetTPCRow0();
5055 Int_t drow = Int_t(2.+0.5/(0.05+kink->GetAngle(2)));
5058 Int_t last = row0-drow;
5059 if (last<40) last=40;
5060 if (last<ktrack0->GetFirstPoint()+25) last = ktrack0->GetFirstPoint()+25;
5061 AliTPCseed* seed0 = ReSeed(ktrack0,last,kFALSE);
5064 Int_t first = row0+drow;
5065 if (first>130) first=130;
5066 if (first>ktrack1->GetLastPoint()-25) first = TMath::Max(ktrack1->GetLastPoint()-25,30);
5067 AliTPCseed* seed1 = ReSeed(ktrack1,first,kTRUE);
5069 if (seed0 && seed1){
5070 kink->SetStatus(1,8);
5071 if (RefitKink(*seed0,*seed1,*kink)) kink->SetStatus(1,9);
5072 row0 = GetRowNumber(kink->GetR());
5073 sumn = seed0->GetNumberOfClusters()+seed1->GetNumberOfClusters();
5074 mothers[i] = *seed0;
5075 daughters[i] = *seed1;
5078 delete kinks->RemoveAt(i);
5079 if (seed0) delete seed0;
5080 if (seed1) delete seed1;
5083 if (kink->GetDistance()>0.5 || kink->GetR()<110 || kink->GetR()>240) {
5084 delete kinks->RemoveAt(i);
5085 if (seed0) delete seed0;
5086 if (seed1) delete seed1;
5094 if (kink) quality[i] = 160*((0.1+kink->GetDistance())*(2.-kink->GetTPCDensityFactor()))/(sumn+40.); //the longest -clossest will win
5096 TMath::Sort(nkinks,quality,indexes,kFALSE);
5098 //remove double find kinks
5100 for (Int_t ikink0=1;ikink0<nkinks;ikink0++){
5101 AliKink * kink0 = (AliKink*) kinks->At(indexes[ikink0]);
5102 if (!kink0) continue;
5104 for (Int_t ikink1=0;ikink1<ikink0;ikink1++){
5105 if (!kink0) continue;
5106 AliKink * kink1 = (AliKink*) kinks->At(indexes[ikink1]);
5107 if (!kink1) continue;
5108 // if not close kink continue
5109 if (TMath::Abs(kink1->GetPosition()[2]-kink0->GetPosition()[2])>10) continue;
5110 if (TMath::Abs(kink1->GetPosition()[1]-kink0->GetPosition()[1])>10) continue;
5111 if (TMath::Abs(kink1->GetPosition()[0]-kink0->GetPosition()[0])>10) continue;
5113 AliTPCseed &mother0 = mothers[indexes[ikink0]];
5114 AliTPCseed &daughter0 = daughters[indexes[ikink0]];
5115 AliTPCseed &mother1 = mothers[indexes[ikink1]];
5116 AliTPCseed &daughter1 = daughters[indexes[ikink1]];
5117 Int_t row0 = (kink0->GetTPCRow0()+kink1->GetTPCRow0())/2;
5126 for (Int_t i=0;i<row0;i++){
5127 if (mother0.GetClusterIndex(i)>0 && mother1.GetClusterIndex(i)>0){
5130 if (mother0.GetClusterIndex(i)==mother1.GetClusterIndex(i)){
5137 for (Int_t i=row0;i<158;i++){
5138 if (daughter0.GetClusterIndex(i)>0 && daughter0.GetClusterIndex(i)>0){
5141 if (mother0.GetClusterIndex(i)==mother1.GetClusterIndex(i)){
5147 Float_t ratio = Float_t(same+1)/Float_t(both+1);
5148 Float_t ratiom = Float_t(samem+1)/Float_t(bothm+1);
5149 Float_t ratiod = Float_t(samed+1)/Float_t(bothd+1);
5150 if (ratio>0.3 && ratiom>0.5 &&ratiod>0.5) {
5151 Int_t sum0 = mother0.GetNumberOfClusters()+daughter0.GetNumberOfClusters();
5152 Int_t sum1 = mother1.GetNumberOfClusters()+daughter1.GetNumberOfClusters();
5154 shared[kink0->GetIndex(0)]= kTRUE;
5155 shared[kink0->GetIndex(1)]= kTRUE;
5156 delete kinks->RemoveAt(indexes[ikink0]);
5159 shared[kink1->GetIndex(0)]= kTRUE;
5160 shared[kink1->GetIndex(1)]= kTRUE;
5161 delete kinks->RemoveAt(indexes[ikink1]);
5168 for (Int_t i=0;i<nkinks;i++){
5169 AliKink * kink = (AliKink*) kinks->At(indexes[i]);
5170 if (!kink) continue;
5171 kink->SetTPCRow0(GetRowNumber(kink->GetR()));
5172 Int_t index0 = kink->GetIndex(0);
5173 Int_t index1 = kink->GetIndex(1);
5174 if (circular[index0]||circular[index1]&&kink->GetDistance()>0.2) continue;
5175 kink->SetMultiple(usage[index0],0);
5176 kink->SetMultiple(usage[index1],1);
5177 if (kink->GetMultiple()[0]+kink->GetMultiple()[1]>2) continue;
5178 if (kink->GetMultiple()[0]+kink->GetMultiple()[1]>0 && quality[indexes[i]]>0.2) continue;
5179 if (kink->GetMultiple()[0]+kink->GetMultiple()[1]>0 && kink->GetDistance()>0.2) continue;
5180 if (circular[index0]||circular[index1]&&kink->GetDistance()>0.1) continue;
5182 AliTPCseed * ktrack0 = (AliTPCseed*)array->At(index0);
5183 AliTPCseed * ktrack1 = (AliTPCseed*)array->At(index1);
5184 if (!ktrack0 || !ktrack1) continue;
5185 Int_t index = esd->AddKink(kink);
5188 if ( ktrack0->GetKinkIndex(0)==0 && ktrack1->GetKinkIndex(0)==0) { //best kink
5189 if (mothers[indexes[i]].GetNumberOfClusters()>20 && daughters[indexes[i]].GetNumberOfClusters()>20 && (mothers[indexes[i]].GetNumberOfClusters()+daughters[indexes[i]].GetNumberOfClusters())>100){
5190 *ktrack0 = mothers[indexes[i]];
5191 *ktrack1 = daughters[indexes[i]];
5195 ktrack0->SetKinkIndex(usage[index0],-(index+1));
5196 ktrack1->SetKinkIndex(usage[index1], (index+1));
5201 // Remove tracks corresponding to shared kink's
5203 for (Int_t i=0;i<nentries;i++){
5204 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5205 if (!track0) continue;
5206 if (track0->GetKinkIndex(0)!=0) continue;
5207 if (shared[i]) delete array->RemoveAt(i);
5212 RemoveUsed2(array,0.5,0.4,30);
5214 for (Int_t i=0;i<nentries;i++){
5215 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5216 if (!track0) continue;
5217 track0->CookdEdx(0.02,0.6);
5221 for (Int_t i=0;i<nentries;i++){
5222 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5223 if (!track0) continue;
5224 if (track0->Pt()<1.4) continue;
5225 //remove double high momenta tracks - overlapped with kink candidates
5228 for (Int_t icl=track0->GetFirstPoint();icl<track0->GetLastPoint(); icl++){
5229 if (track0->GetClusterPointer(icl)!=0){
5231 if (track0->GetClusterPointer(icl)->IsUsed(10)) shared++;
5234 if (Float_t(shared+1)/Float_t(all+1)>0.5) {
5235 delete array->RemoveAt(i);
5239 if (track0->GetKinkIndex(0)!=0) continue;
5240 if (track0->GetNumberOfClusters()<80) continue;
5242 AliTPCseed *pmother = new AliTPCseed();
5243 AliTPCseed *pdaughter = new AliTPCseed();
5244 AliKink *pkink = new AliKink;
5246 AliTPCseed & mother = *pmother;
5247 AliTPCseed & daughter = *pdaughter;
5248 AliKink & kink = *pkink;
5249 if (CheckKinkPoint(track0,mother,daughter, kink)){
5250 if (mother.GetNumberOfClusters()<30||daughter.GetNumberOfClusters()<20) {
5254 continue; //too short tracks
5256 if (mother.Pt()<1.4) {
5262 Int_t row0= kink.GetTPCRow0();
5263 if (kink.GetDistance()>0.5 || kink.GetR()<110. || kink.GetR()>240.) {
5270 Int_t index = esd->AddKink(&kink);
5271 mother.SetKinkIndex(0,-(index+1));
5272 daughter.SetKinkIndex(0,index+1);
5273 if (mother.GetNumberOfClusters()>50) {
5274 delete array->RemoveAt(i);
5275 array->AddAt(new AliTPCseed(mother),i);
5278 array->AddLast(new AliTPCseed(mother));
5280 array->AddLast(new AliTPCseed(daughter));
5281 for (Int_t icl=0;icl<row0;icl++) {
5282 if (mother.GetClusterPointer(icl)) mother.GetClusterPointer(icl)->Use(20);
5285 for (Int_t icl=row0;icl<158;icl++) {
5286 if (daughter.GetClusterPointer(icl)) daughter.GetClusterPointer(icl)->Use(20);
5295 delete [] daughters;
5317 printf("Ncandidates=\t%d\t%d\t%d\t%d\n",esd->GetNumberOfKinks(),ncandidates,ntracks,nall);
5321 void AliTPCtrackerMI::FindV0s(TObjArray * array, AliESDEvent *esd)
5327 TObjArray *tpcv0s = new TObjArray(100000);
5328 Int_t nentries = array->GetEntriesFast();
5329 AliHelix *helixes = new AliHelix[nentries];
5330 Int_t *sign = new Int_t[nentries];
5331 Float_t *alpha = new Float_t[nentries];
5332 Float_t *z0 = new Float_t[nentries];
5333 Float_t *dca = new Float_t[nentries];
5334 Float_t *sdcar = new Float_t[nentries];
5335 Float_t *cdcar = new Float_t[nentries];
5336 Float_t *pulldcar = new Float_t[nentries];
5337 Float_t *pulldcaz = new Float_t[nentries];
5338 Float_t *pulldca = new Float_t[nentries];
5339 Bool_t *isPrim = new Bool_t[nentries];
5340 const AliESDVertex * primvertex = esd->GetVertex();
5341 Double_t zvertex = primvertex->GetZv();
5343 // nentries = array->GetEntriesFast();
5345 for (Int_t i=0;i<nentries;i++){
5348 AliTPCseed* track = (AliTPCseed*)array->At(i);
5349 if (!track) continue;
5350 track->GetV0Indexes()[0] = 0; //rest v0 indexes
5351 track->GetV0Indexes()[1] = 0; //rest v0 indexes
5352 track->GetV0Indexes()[2] = 0; //rest v0 indexes
5354 alpha[i] = track->GetAlpha();
5355 new (&helixes[i]) AliHelix(*track);
5357 helixes[i].Evaluate(0,xyz);
5358 sign[i] = (track->GetC()>0) ? -1:1;
5362 if (track->GetProlongation(0,y,z)) z0[i] = z;
5363 dca[i] = track->GetD(0,0);
5365 // dca error parrameterezation + pulls
5367 sdcar[i] = TMath::Sqrt(0.150*0.150+(100*track->GetC())*(100*track->GetC()));
5368 if (TMath::Abs(track->GetTgl())>1) sdcar[i]*=2.5;
5369 cdcar[i] = TMath::Exp((TMath::Abs(track->GetC())-0.0106)*525.3);
5370 pulldcar[i] = (dca[i]-cdcar[i])/sdcar[i];
5371 pulldcaz[i] = (z0[i]-zvertex)/sdcar[i];
5372 pulldca[i] = TMath::Sqrt(pulldcar[i]*pulldcar[i]+pulldcaz[i]*pulldcaz[i]);
5373 if (track->TPCrPID(1)+track->TPCrPID(2)+track->TPCrPID(3)>0.5) {
5374 if (pulldca[i]<3.) isPrim[i]=kTRUE; //pion, muon and Kaon 3 sigma cut
5376 if (track->TPCrPID(4)>0.5) {
5377 if (pulldca[i]<0.5) isPrim[i]=kTRUE; //proton 0.5 sigma cut
5379 if (track->TPCrPID(0)>0.4) {
5380 isPrim[i]=kFALSE; //electron no sigma cut
5387 Int_t ncandidates =0;
5390 Double_t phase[2][2],radius[2];
5396 TTreeSRedirector &cstream = *fDebugStreamer;
5397 Float_t fprimvertex[3]={GetX(),GetY(),GetZ()};
5399 Double_t cradius0 = 10*10;
5400 Double_t cradius1 = 200*200;
5403 Double_t cpointAngle = 0.95;
5405 Double_t delta[2]={10000,10000};
5406 for (Int_t i =0;i<nentries;i++){
5407 if (sign[i]==0) continue;
5408 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
5409 if (!track0) continue;
5410 if (AliTPCReconstructor::StreamLevel()>1){
5415 "zvertex="<<zvertex<<
5416 "sdcar0="<<sdcar[i]<<
5417 "cdcar0="<<cdcar[i]<<
5418 "pulldcar0="<<pulldcar[i]<<
5419 "pulldcaz0="<<pulldcaz[i]<<
5420 "pulldca0="<<pulldca[i]<<
5421 "isPrim="<<isPrim[i]<<
5425 if (track0->GetSigned1Pt()<0) continue;
5426 if (track0->GetKinkIndex(0)>0||isPrim[i]) continue; //daughter kink
5428 if (TMath::Abs(helixes[i].GetHelix(4))<0.000000001) continue;
5433 for (Int_t j =0;j<nentries;j++){
5434 AliTPCseed * track1 = (AliTPCseed*)array->At(j);
5435 if (!track1) continue;
5436 if (track1->GetKinkIndex(0)>0 || isPrim[j]) continue; //daughter kink
5437 if (sign[j]*sign[i]>0) continue;
5438 if (TMath::Abs(helixes[j].GetHelix(4))<0.000001) continue;
5439 if (track0->GetCircular()+track1->GetCircular()>1) continue; //circling -returning track
5442 // DCA to prim vertex cut
5448 Int_t npoints = helixes[i].GetRPHIintersections(helixes[j], phase, radius,cdist2);
5449 if (npoints<1) continue;
5453 if (radius[0]<cradius0||radius[0]>cradius1) continue;
5454 helixes[i].LinearDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta[0]);
5455 if (delta[0]>cdist1) continue;
5458 if (TMath::Max(radius[0],radius[1])<cradius0|| TMath::Min(radius[0],radius[1])>cradius1) continue;
5459 helixes[i].LinearDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta[0]);
5460 helixes[i].LinearDCA(helixes[j],phase[1][0],phase[1][1],radius[1],delta[1]);
5461 if (delta[1]<delta[0]) iclosest=1;
5462 if (delta[iclosest]>cdist1) continue;
5464 helixes[i].ParabolicDCA(helixes[j],phase[iclosest][0],phase[iclosest][1],radius[iclosest],delta[iclosest]);
5465 if (radius[iclosest]<cradius0 || radius[iclosest]>cradius1 || delta[iclosest]>cdist1) continue;
5467 Double_t pointAngle = helixes[i].GetPointAngle(helixes[j],phase[iclosest],fprimvertex);
5468 if (pointAngle<cpointAngle) continue;
5470 Bool_t isGamma = kFALSE;
5471 vertex.SetParamP(*track0); //track0 - plus
5472 vertex.SetParamN(*track1); //track1 - minus
5473 vertex.Update(fprimvertex);
5474 if (track0->TPCrPID(0)>0.3&&track1->TPCrPID(0)>0.3&&vertex.GetAnglep()[2]<0.15) isGamma=kTRUE; // gamma conversion candidate
5475 Double_t pointAngle2 = vertex.GetV0CosineOfPointingAngle();
5477 if (vertex.GetV0CosineOfPointingAngle()<cpointAngle && (!isGamma)) continue;// point angle cut
5478 //Bo: if (vertex.GetDist2()>2&&(!isGamma)) continue; // point angle cut
5479 if (vertex.GetDcaV0Daughters()>2&&(!isGamma)) continue;//Bo: // point angle cut
5480 Float_t sigmae = 0.15*0.15;
5481 if (vertex.GetRr()<80)
5482 sigmae += (sdcar[i]*sdcar[i]+sdcar[j]*sdcar[j])*(1.-vertex.GetRr()/80.)*(1.-vertex.GetRr()/80.);
5483 sigmae+= TMath::Sqrt(sigmae);
5484 //Bo: if (vertex.GetDist2()/sigmae>3.&&(!isGamma)) continue;
5485 if (vertex.GetDcaV0Daughters()/sigmae>3.&&(!isGamma)) continue;
5486 Float_t densb0=0,densb1=0,densa0=0,densa1=0;
5487 Int_t row0 = GetRowNumber(vertex.GetRr());
5489 //Bo: if (vertex.GetDist2()>0.2) continue;
5490 if (vertex.GetDcaV0Daughters()>0.2) continue;
5491 densb0 = track0->Density2(0,row0-5);
5492 densb1 = track1->Density2(0,row0-5);
5493 if (densb0>0.3|| densb1>0.3) continue; //clusters before vertex
5494 densa0 = track0->Density2(row0+5,row0+40);
5495 densa1 = track1->Density2(row0+5,row0+40);
5496 if ((densa0<0.4|| densa1<0.4)&&(!isGamma)) continue; //missing clusters after vertex
5499 densa0 = track0->Density2(0,40); //cluster density
5500 densa1 = track1->Density2(0,40); //cluster density
5501 if ((vertex.GetRr()<80&&densa0+densa1<1.)&&(!isGamma)) continue;
5503 //Bo: vertex.SetLab(0,track0->GetLabel());
5504 //Bo: vertex.SetLab(1,track1->GetLabel());
5505 vertex.SetChi2After((densa0+densa1)*0.5);
5506 vertex.SetChi2Before((densb0+densb1)*0.5);
5507 vertex.SetIndex(0,i);
5508 vertex.SetIndex(1,j);
5509 //Bo: vertex.SetStatus(1); // TPC v0 candidate
5510 vertex.SetOnFlyStatus(2);//Bo: // TPC v0 candidate
5511 //Bo: vertex.SetRp(track0->TPCrPIDs());
5512 //Bo: vertex.SetRm(track1->TPCrPIDs());
5513 tpcv0s->AddLast(new AliESDv0(vertex));
5516 Int_t eventNr = esd->GetEventNumberInFile(); // This is most likely NOT the event number you'd like to use. It has nothing to do with the 'real' event number
5517 Double_t radiusm= (delta[0]<delta[1])? TMath::Sqrt(radius[0]):TMath::Sqrt(radius[1]);
5518 Double_t deltam= (delta[0]<delta[1])? TMath::Sqrt(delta[0]):TMath::Sqrt(delta[1]);
5519 if (AliTPCReconstructor::StreamLevel()>1) {
5520 Int_t lab0=track0->GetLabel();
5521 Int_t lab1=track1->GetLabel();
5522 Char_t c0=track0->GetCircular();
5523 Char_t c1=track1->GetCircular();
5526 "vertex.="<<&vertex<<
5529 "Helix0.="<<&helixes[i]<<
5532 "Helix1.="<<&helixes[j]<<
5533 "pointAngle="<<pointAngle<<
5534 "pointAngle2="<<pointAngle2<<
5539 "zvertex="<<zvertex<<
5542 "npoints="<<npoints<<
5543 "radius0="<<radius[0]<<
5544 "delta0="<<delta[0]<<
5545 "radius1="<<radius[1]<<
5546 "delta1="<<delta[1]<<
5547 "radiusm="<<radiusm<<
5549 "sdcar0="<<sdcar[i]<<
5550 "sdcar1="<<sdcar[j]<<
5551 "cdcar0="<<cdcar[i]<<
5552 "cdcar1="<<cdcar[j]<<
5553 "pulldcar0="<<pulldcar[i]<<
5554 "pulldcar1="<<pulldcar[j]<<
5555 "pulldcaz0="<<pulldcaz[i]<<
5556 "pulldcaz1="<<pulldcaz[j]<<
5557 "pulldca0="<<pulldca[i]<<
5558 "pulldca1="<<pulldca[j]<<
5568 Float_t *quality = new Float_t[ncandidates];
5569 Int_t *indexes = new Int_t[ncandidates];
5571 for (Int_t i=0;i<ncandidates;i++){
5573 AliESDv0 *v0 = (AliESDv0*)tpcv0s->At(i);
5574 quality[i] = 1./(1.00001-v0->GetV0CosineOfPointingAngle()); //base point angle
5575 // quality[i] /= (0.5+v0->GetDist2());
5576 // quality[i] *= v0->GetChi2After(); //density factor
5578 Int_t index0 = v0->GetIndex(0);
5579 Int_t index1 = v0->GetIndex(1);
5580 //Bo: Double_t minpulldca = TMath::Min(2.+pulldca[v0->GetIndex(0)],(2.+pulldca[v0->GetIndex(1)]) ); //pull
5581 Double_t minpulldca = TMath::Min(2.+pulldca[index0],(2.+pulldca[index1]) );//Bo:
5585 AliTPCseed * track0 = (AliTPCseed*)array->At(index0);
5586 AliTPCseed * track1 = (AliTPCseed*)array->At(index1);
5587 if (track0->TPCrPID(0)>0.3&&track1->TPCrPID(0)>0.3&&v0->GetAnglep()[2]<0.15) quality[i]+=1000000; // gamma conversion candidate
5588 if (track0->TPCrPID(4)>0.9||track1->TPCrPID(4)>0.9&&minpulldca>4) quality[i]*=10; // lambda candidate candidate
5591 TMath::Sort(ncandidates,quality,indexes,kTRUE);
5594 for (Int_t i=0;i<ncandidates;i++){
5595 AliESDv0 * v0 = (AliESDv0*)tpcv0s->At(indexes[i]);
5597 Int_t index0 = v0->GetIndex(0);
5598 Int_t index1 = v0->GetIndex(1);
5599 AliTPCseed * track0 = (AliTPCseed*)array->At(index0);
5600 AliTPCseed * track1 = (AliTPCseed*)array->At(index1);
5601 if (!track0||!track1) {
5605 Bool_t accept =kTRUE; //default accept
5606 Int_t *v0indexes0 = track0->GetV0Indexes();
5607 Int_t *v0indexes1 = track1->GetV0Indexes();
5609 Int_t order0 = (v0indexes0[0]!=0) ? 1:0;
5610 Int_t order1 = (v0indexes1[0]!=0) ? 1:0;
5611 if (v0indexes0[1]!=0) order0 =2;
5612 if (v0indexes1[1]!=0) order1 =2;
5614 if (v0indexes0[2]!=0) {order0=3; accept=kFALSE;}
5615 if (v0indexes0[2]!=0) {order1=3; accept=kFALSE;}
5617 AliESDv0 * v02 = v0;
5619 //Bo: v0->SetOrder(0,order0);
5620 //Bo: v0->SetOrder(1,order1);
5621 //Bo: v0->SetOrder(1,order0+order1);
5622 v0->SetOnFlyStatus(kTRUE);
5623 Int_t index = esd->AddV0(v0);
5624 v02 = esd->GetV0(index);
5625 v0indexes0[order0]=index;
5626 v0indexes1[order1]=index;
5630 Int_t eventNr = esd->GetEventNumberInFile(); // This is most likely NOT the event number you'd like to use. It has nothing to do with the 'real' event number
5631 if (AliTPCReconstructor::StreamLevel()>1) {
5632 Int_t lab0=track0->GetLabel();
5633 Int_t lab1=track1->GetLabel();
5642 "dca0="<<dca[index0]<<
5643 "dca1="<<dca[index1]<<
5647 "quality="<<quality[i]<<
5648 "pulldca0="<<pulldca[index0]<<
5649 "pulldca1="<<pulldca[index1]<<
5673 printf("TPC V0 finder : naccepted\t%d\tncandidates\t%d\tntracks\t%d\tnall\t%d\n",naccepted,ncandidates,ntracks,nall);
5677 Int_t AliTPCtrackerMI::RefitKink(AliTPCseed &mother, AliTPCseed &daughter, AliESDkink &knk)
5680 // refit kink towards to the vertex
5683 AliKink &kink=(AliKink &)knk;
5685 Int_t row0 = GetRowNumber(kink.GetR());
5686 FollowProlongation(mother,0);
5687 mother.Reset(kFALSE);
5689 FollowProlongation(daughter,row0);
5690 daughter.Reset(kFALSE);
5691 FollowBackProlongation(daughter,158);
5692 daughter.Reset(kFALSE);
5693 Int_t first = TMath::Max(row0-20,30);
5694 Int_t last = TMath::Min(row0+20,140);
5696 const Int_t kNdiv =5;
5697 AliTPCseed param0[kNdiv]; // parameters along the track
5698 AliTPCseed param1[kNdiv]; // parameters along the track
5699 AliKink kinks[kNdiv]; // corresponding kink parameters
5702 for (Int_t irow=0; irow<kNdiv;irow++){
5703 rows[irow] = first +((last-first)*irow)/(kNdiv-1);
5705 // store parameters along the track
5707 for (Int_t irow=0;irow<kNdiv;irow++){
5708 FollowBackProlongation(mother, rows[irow]);
5709 FollowProlongation(daughter,rows[kNdiv-1-irow]);
5710 param0[irow] = mother;
5711 param1[kNdiv-1-irow] = daughter;
5715 for (Int_t irow=0; irow<kNdiv-1;irow++){
5716 if (param0[irow].GetNumberOfClusters()<kNdiv||param1[irow].GetNumberOfClusters()<kNdiv) continue;
5717 kinks[irow].SetMother(param0[irow]);
5718 kinks[irow].SetDaughter(param1[irow]);
5719 kinks[irow].Update();
5722 // choose kink with best "quality"
5724 Double_t mindist = 10000;
5725 for (Int_t irow=0;irow<kNdiv;irow++){
5726 if (param0[irow].GetNumberOfClusters()<20||param1[irow].GetNumberOfClusters()<20) continue;
5727 if (TMath::Abs(kinks[irow].GetR())>240.) continue;
5728 if (TMath::Abs(kinks[irow].GetR())<100.) continue;
5730 Float_t normdist = TMath::Abs(param0[irow].GetX()-kinks[irow].GetR())*(0.1+kink.GetDistance());
5731 normdist/= (param0[irow].GetNumberOfClusters()+param1[irow].GetNumberOfClusters()+40.);
5732 if (normdist < mindist){
5738 if (index==-1) return 0;
5741 param0[index].Reset(kTRUE);
5742 FollowProlongation(param0[index],0);
5744 mother = param0[index];
5745 daughter = param1[index]; // daughter in vertex
5747 kink.SetMother(mother);
5748 kink.SetDaughter(daughter);
5750 kink.SetTPCRow0(GetRowNumber(kink.GetR()));
5751 kink.SetTPCncls(param0[index].GetNumberOfClusters(),0);
5752 kink.SetTPCncls(param1[index].GetNumberOfClusters(),1);
5753 kink.SetLabel(CookLabel(&mother,0.4, 0,kink.GetTPCRow0()),0);
5754 kink.SetLabel(CookLabel(&daughter,0.4, kink.GetTPCRow0(),160),1);
5755 mother.SetLabel(kink.GetLabel(0));
5756 daughter.SetLabel(kink.GetLabel(1));
5762 void AliTPCtrackerMI::UpdateKinkQualityM(AliTPCseed * seed){
5764 // update Kink quality information for mother after back propagation
5766 if (seed->GetKinkIndex(0)>=0) return;
5767 for (Int_t ikink=0;ikink<3;ikink++){
5768 Int_t index = seed->GetKinkIndex(ikink);
5769 if (index>=0) break;
5770 index = TMath::Abs(index)-1;
5771 AliESDkink * kink = fEvent->GetKink(index);
5772 kink->SetTPCDensity(-1,0,0);
5773 kink->SetTPCDensity(1,0,1);
5775 Int_t row0 = kink->GetTPCRow0() - 2 - Int_t( 0.5/ (0.05+kink->GetAngle(2)));
5776 if (row0<15) row0=15;
5778 Int_t row1 = kink->GetTPCRow0() + 2 + Int_t( 0.5/ (0.05+kink->GetAngle(2)));
5779 if (row1>145) row1=145;
5781 Int_t found,foundable,shared;
5782 seed->GetClusterStatistic(0,row0, found, foundable,shared,kFALSE);
5783 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),0,0);
5784 seed->GetClusterStatistic(row1,155, found, foundable,shared,kFALSE);
5785 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),0,1);
5790 void AliTPCtrackerMI::UpdateKinkQualityD(AliTPCseed * seed){
5792 // update Kink quality information for daughter after refit
5794 if (seed->GetKinkIndex(0)<=0) return;
5795 for (Int_t ikink=0;ikink<3;ikink++){
5796 Int_t index = seed->GetKinkIndex(ikink);
5797 if (index<=0) break;
5798 index = TMath::Abs(index)-1;
5799 AliESDkink * kink = fEvent->GetKink(index);
5800 kink->SetTPCDensity(-1,1,0);
5801 kink->SetTPCDensity(-1,1,1);
5803 Int_t row0 = kink->GetTPCRow0() -2 - Int_t( 0.5/ (0.05+kink->GetAngle(2)));
5804 if (row0<15) row0=15;
5806 Int_t row1 = kink->GetTPCRow0() +2 + Int_t( 0.5/ (0.05+kink->GetAngle(2)));
5807 if (row1>145) row1=145;
5809 Int_t found,foundable,shared;
5810 seed->GetClusterStatistic(0,row0, found, foundable,shared,kFALSE);
5811 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),1,0);
5812 seed->GetClusterStatistic(row1,155, found, foundable,shared,kFALSE);
5813 if (foundable>5) kink->SetTPCDensity(Float_t(found)/Float_t(foundable),1,1);
5819 Int_t AliTPCtrackerMI::CheckKinkPoint(AliTPCseed*seed,AliTPCseed &mother, AliTPCseed &daughter, AliESDkink &knk)
5822 // check kink point for given track
5823 // if return value=0 kink point not found
5824 // otherwise seed0 correspond to mother particle
5825 // seed1 correspond to daughter particle
5826 // kink parameter of kink point
5827 AliKink &kink=(AliKink &)knk;
5829 Int_t middlerow = (seed->GetFirstPoint()+seed->GetLastPoint())/2;
5830 Int_t first = seed->GetFirstPoint();
5831 Int_t last = seed->GetLastPoint();
5832 if (last-first<20) return 0; // shortest length - 2*30 = 60 pad-rows
5835 AliTPCseed *seed1 = ReSeed(seed,middlerow+20, kTRUE); //middle of chamber
5836 if (!seed1) return 0;
5837 FollowProlongation(*seed1,seed->GetLastPoint()-20);
5838 seed1->Reset(kTRUE);
5839 FollowProlongation(*seed1,158);
5840 seed1->Reset(kTRUE);
5841 last = seed1->GetLastPoint();
5843 AliTPCseed *seed0 = new AliTPCseed(*seed);
5844 seed0->Reset(kFALSE);
5847 AliTPCseed param0[20]; // parameters along the track
5848 AliTPCseed param1[20]; // parameters along the track
5849 AliKink kinks[20]; // corresponding kink parameters
5851 for (Int_t irow=0; irow<20;irow++){
5852 rows[irow] = first +((last-first)*irow)/19;
5854 // store parameters along the track
5856 for (Int_t irow=0;irow<20;irow++){
5857 FollowBackProlongation(*seed0, rows[irow]);
5858 FollowProlongation(*seed1,rows[19-irow]);
5859 param0[irow] = *seed0;
5860 param1[19-irow] = *seed1;
5864 for (Int_t irow=0; irow<19;irow++){
5865 kinks[irow].SetMother(param0[irow]);
5866 kinks[irow].SetDaughter(param1[irow]);
5867 kinks[irow].Update();
5870 // choose kink with biggest change of angle
5872 Double_t maxchange= 0;
5873 for (Int_t irow=1;irow<19;irow++){
5874 if (TMath::Abs(kinks[irow].GetR())>240.) continue;
5875 if (TMath::Abs(kinks[irow].GetR())<110.) continue;
5876 Float_t quality = TMath::Abs(kinks[irow].GetAngle(2))/(3.+TMath::Abs(kinks[irow].GetR()-param0[irow].GetX()));
5877 if ( quality > maxchange){
5878 maxchange = quality;
5885 if (index<0) return 0;
5887 Int_t row0 = GetRowNumber(kinks[index].GetR()); //row 0 estimate
5888 seed0 = new AliTPCseed(param0[index]);
5889 seed1 = new AliTPCseed(param1[index]);
5890 seed0->Reset(kFALSE);
5891 seed1->Reset(kFALSE);
5892 seed0->ResetCovariance(10.);
5893 seed1->ResetCovariance(10.);
5894 FollowProlongation(*seed0,0);
5895 FollowBackProlongation(*seed1,158);
5896 mother = *seed0; // backup mother at position 0
5897 seed0->Reset(kFALSE);
5898 seed1->Reset(kFALSE);
5899 seed0->ResetCovariance(10.);
5900 seed1->ResetCovariance(10.);
5902 first = TMath::Max(row0-20,0);
5903 last = TMath::Min(row0+20,158);
5905 for (Int_t irow=0; irow<20;irow++){
5906 rows[irow] = first +((last-first)*irow)/19;
5908 // store parameters along the track
5910 for (Int_t irow=0;irow<20;irow++){
5911 FollowBackProlongation(*seed0, rows[irow]);
5912 FollowProlongation(*seed1,rows[19-irow]);
5913 param0[irow] = *seed0;
5914 param1[19-irow] = *seed1;
5918 for (Int_t irow=0; irow<19;irow++){
5919 kinks[irow].SetMother(param0[irow]);
5920 kinks[irow].SetDaughter(param1[irow]);
5921 // param0[irow].Dump();
5922 //param1[irow].Dump();
5923 kinks[irow].Update();
5926 // choose kink with biggest change of angle
5929 for (Int_t irow=0;irow<20;irow++){
5930 if (TMath::Abs(kinks[irow].GetR())>250.) continue;
5931 if (TMath::Abs(kinks[irow].GetR())<90.) continue;
5932 Float_t quality = TMath::Abs(kinks[irow].GetAngle(2))/(3.+TMath::Abs(kinks[irow].GetR()-param0[irow].GetX()));
5933 if ( quality > maxchange){
5934 maxchange = quality;
5941 if (index==-1 || param0[index].GetNumberOfClusters()+param1[index].GetNumberOfClusters()<100){
5946 // Float_t anglesigma = TMath::Sqrt(param0[index].fC22+param0[index].fC33+param1[index].fC22+param1[index].fC33);
5948 kink.SetMother(param0[index]);
5949 kink.SetDaughter(param1[index]);
5951 row0 = GetRowNumber(kink.GetR());
5952 kink.SetTPCRow0(row0);
5953 kink.SetLabel(CookLabel(seed0,0.5,0,row0),0);
5954 kink.SetLabel(CookLabel(seed1,0.5,row0,158),1);
5955 kink.SetIndex(-10,0);
5956 kink.SetIndex(int(param0[index].GetNumberOfClusters()+param1[index].GetNumberOfClusters()),1);
5957 kink.SetTPCncls(param0[index].GetNumberOfClusters(),0);
5958 kink.SetTPCncls(param1[index].GetNumberOfClusters(),1);
5961 // new (&mother) AliTPCseed(param0[index]);
5962 daughter = param1[index];
5963 daughter.SetLabel(kink.GetLabel(1));
5964 param0[index].Reset(kTRUE);
5965 FollowProlongation(param0[index],0);
5966 mother = param0[index];
5967 mother.SetLabel(kink.GetLabel(0));
5977 AliTPCseed* AliTPCtrackerMI::ReSeed(AliTPCseed *t)
5980 // reseed - refit - track
5983 // Int_t last = fSectors->GetNRows()-1;
5985 if (fSectors == fOuterSec){
5986 first = TMath::Max(first, t->GetFirstPoint()-fInnerSec->GetNRows());
5990 first = t->GetFirstPoint();
5992 AliTPCseed * seed = MakeSeed(t,0.1,0.5,0.9);
5993 FollowBackProlongation(*t,fSectors->GetNRows()-1);
5995 FollowProlongation(*t,first);
6005 //_____________________________________________________________________________
6006 Int_t AliTPCtrackerMI::ReadSeeds(const TFile *inp) {
6007 //-----------------------------------------------------------------
6008 // This function reades track seeds.
6009 //-----------------------------------------------------------------
6010 TDirectory *savedir=gDirectory;
6012 TFile *in=(TFile*)inp;
6013 if (!in->IsOpen()) {
6014 cerr<<"AliTPCtrackerMI::ReadSeeds(): input file is not open !\n";
6019 TTree *seedTree=(TTree*)in->Get("Seeds");
6021 cerr<<"AliTPCtrackerMI::ReadSeeds(): ";
6022 cerr<<"can't get a tree with track seeds !\n";
6025 AliTPCtrack *seed=new AliTPCtrack;
6026 seedTree->SetBranchAddress("tracks",&seed);
6028 if (fSeeds==0) fSeeds=new TObjArray(15000);
6030 Int_t n=(Int_t)seedTree->GetEntries();
6031 for (Int_t i=0; i<n; i++) {
6032 seedTree->GetEvent(i);
6033 fSeeds->AddLast(new AliTPCseed(*seed/*,seed->GetAlpha()*/));
6042 Int_t AliTPCtrackerMI::Clusters2Tracks (AliESDEvent *esd)
6045 if (fSeeds) DeleteSeeds();
6048 if (!fSeeds) return 1;
6055 //_____________________________________________________________________________
6056 Int_t AliTPCtrackerMI::Clusters2Tracks() {
6057 //-----------------------------------------------------------------
6058 // This is a track finder.
6059 //-----------------------------------------------------------------
6060 TDirectory *savedir=gDirectory;
6064 fSeeds = Tracking();
6067 Info("Clusters2Tracks","Time for tracking: \t");timer.Print();timer.Start();
6069 //activate again some tracks
6070 for (Int_t i=0; i<fSeeds->GetEntriesFast(); i++) {
6071 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
6073 Int_t nc=t.GetNumberOfClusters();
6075 delete fSeeds->RemoveAt(i);
6079 if (pt->GetRemoval()==10) {
6080 if (pt->GetDensityFirst(20)>0.8 || pt->GetDensityFirst(30)>0.8 || pt->GetDensityFirst(40)>0.7)
6081 pt->Desactivate(10); // make track again active
6083 pt->Desactivate(20);
6084 delete fSeeds->RemoveAt(i);
6089 RemoveUsed2(fSeeds,0.85,0.85,0);
6090 if (AliTPCReconstructor::GetRecoParam()->GetDoKinks()) FindKinks(fSeeds,fEvent);
6091 //FindCurling(fSeeds, fEvent,0);
6092 if (AliTPCReconstructor::StreamLevel()>2) FindMultiMC(fSeeds, fEvent,0); // find multi found tracks
6093 RemoveUsed2(fSeeds,0.5,0.4,20);
6094 FindSplitted(fSeeds, fEvent,0); // find multi found tracks
6097 // // refit short tracks
6099 Int_t nseed=fSeeds->GetEntriesFast();
6102 for (Int_t i=0; i<nseed; i++) {
6103 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
6105 Int_t nc=t.GetNumberOfClusters();
6107 delete fSeeds->RemoveAt(i);
6110 CookLabel(pt,0.1); //For comparison only
6111 //if ((pt->IsActive() || (pt->fRemoval==10) )&& nc>50 &&pt->GetNumberOfClusters()>0.4*pt->fNFoundable){
6112 if ((pt->IsActive() || (pt->GetRemoval()==10) )){
6114 if (fDebug>0) cerr<<found<<'\r';
6118 delete fSeeds->RemoveAt(i);
6122 //RemoveOverlap(fSeeds,0.99,7,kTRUE);
6124 //RemoveUsed(fSeeds,0.9,0.9,6);
6126 nseed=fSeeds->GetEntriesFast();
6128 for (Int_t i=0; i<nseed; i++) {
6129 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
6131 Int_t nc=t.GetNumberOfClusters();
6133 delete fSeeds->RemoveAt(i);
6137 t.CookdEdx(0.02,0.6);
6138 // CheckKinkPoint(&t,0.05);
6139 //if ((pt->IsActive() || (pt->fRemoval==10) )&& nc>50 &&pt->GetNumberOfClusters()>0.4*pt->fNFoundable){
6140 if ((pt->IsActive() || (pt->GetRemoval()==10) )){
6148 delete fSeeds->RemoveAt(i);
6149 //AliTPCseed * seed1 = ReSeed(pt,0.05,0.5,1);
6151 // FollowProlongation(*seed1,0);
6152 // Int_t n = seed1->GetNumberOfClusters();
6153 // printf("fP4\t%f\t%f\n",seed1->GetC(),pt->GetC());
6154 // printf("fN\t%d\t%d\n", seed1->GetNumberOfClusters(),pt->GetNumberOfClusters());
6157 //AliTPCseed * seed2 = ReSeed(pt,0.95,0.5,0.05);
6161 SortTracks(fSeeds, 1);
6165 PrepareForBackProlongation(fSeeds,5.);
6166 PropagateBack(fSeeds);
6167 printf("Time for back propagation: \t");timer.Print();timer.Start();
6171 PrepareForProlongation(fSeeds,5.);
6172 PropagateForward2(fSeeds);
6174 printf("Time for FORWARD propagation: \t");timer.Print();timer.Start();
6175 // RemoveUsed(fSeeds,0.7,0.7,6);
6176 //RemoveOverlap(fSeeds,0.9,7,kTRUE);
6178 nseed=fSeeds->GetEntriesFast();
6180 for (Int_t i=0; i<nseed; i++) {
6181 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
6183 Int_t nc=t.GetNumberOfClusters();
6185 delete fSeeds->RemoveAt(i);
6188 t.CookdEdx(0.02,0.6);
6189 // CookLabel(pt,0.1); //For comparison only
6190 //if ((pt->IsActive() || (pt->fRemoval==10) )&& nc>50 &&pt->GetNumberOfClusters()>0.4*pt->fNFoundable){
6191 if ((pt->IsActive() || (pt->fRemoval==10) )){
6192 cerr<<found++<<'\r';
6195 delete fSeeds->RemoveAt(i);
6200 // fNTracks = found;
6202 Info("Clusters2Tracks","Time for overlap removal, track writing and dedx cooking: \t"); timer.Print();timer.Start();
6205 // cerr<<"Number of found tracks : "<<"\t"<<found<<endl;
6206 Info("Clusters2Tracks","Number of found tracks %d",found);
6208 // UnloadClusters();
6213 void AliTPCtrackerMI::Tracking(TObjArray * arr)
6216 // tracking of the seeds
6219 fSectors = fOuterSec;
6220 ParallelTracking(arr,150,63);
6221 fSectors = fOuterSec;
6222 ParallelTracking(arr,63,0);
6225 TObjArray * AliTPCtrackerMI::Tracking(Int_t seedtype, Int_t i1, Int_t i2, Float_t cuts[4], Float_t dy, Int_t dsec)
6230 TObjArray * arr = new TObjArray;
6232 fSectors = fOuterSec;
6235 for (Int_t sec=0;sec<fkNOS;sec++){
6236 if (seedtype==3) MakeSeeds3(arr,sec,i1,i2,cuts,dy, dsec);
6237 if (seedtype==4) MakeSeeds5(arr,sec,i1,i2,cuts,dy);
6238 if (seedtype==2) MakeSeeds2(arr,sec,i1,i2,cuts,dy);
6241 Info("Tracking","\nSeeding - %d\t%d\t%d\t%d\n",seedtype,i1,i2,arr->GetEntriesFast());
6253 TObjArray * AliTPCtrackerMI::Tracking()
6257 if (AliTPCReconstructor::GetRecoParam()->GetSpecialSeeding()) return TrackingSpecial();
6260 Int_t nup=fOuterSec->GetNRows()+fInnerSec->GetNRows();
6262 TObjArray * seeds = new TObjArray;
6271 Float_t fnumber = 3.0;
6272 Float_t fdensity = 3.0;
6277 for (Int_t delta = 0; delta<18; delta+=6){
6281 arr = Tracking(3,nup-1-delta,nup-1-delta-gap,cuts,-1,1);
6282 SumTracks(seeds,arr);
6283 SignClusters(seeds,fnumber,fdensity);
6285 for (Int_t i=2;i<6;i+=2){
6286 // seed high pt tracks
6289 arr = Tracking(3,nup-i-delta,nup-i-delta-gap,cuts,-1,0);
6290 SumTracks(seeds,arr);
6291 SignClusters(seeds,fnumber,fdensity);
6296 // RemoveUsed(seeds,0.9,0.9,1);
6297 // UnsignClusters();
6298 // SignClusters(seeds,fnumber,fdensity);
6302 for (Int_t delta = 20; delta<120; delta+=10){
6304 // seed high pt tracks
6308 arr = Tracking(3,nup-delta,nup-delta-gap,cuts,-1);
6309 SumTracks(seeds,arr);
6310 SignClusters(seeds,fnumber,fdensity);
6315 arr = Tracking(3,nup-delta-5,nup-delta-5-gap,cuts,-1);
6316 SumTracks(seeds,arr);
6317 SignClusters(seeds,fnumber,fdensity);
6328 Info("Tracking()","\n\nPrimary seeding\t%d\n\n",seeds->GetEntriesFast());
6332 // RemoveUsed(seeds,0.75,0.75,1);
6334 //SignClusters(seeds,fnumber,fdensity);
6343 arr = Tracking(4,nup-1,nup-1-gap,cuts,-1);
6344 SumTracks(seeds,arr);
6345 SignClusters(seeds,fnumber,fdensity);
6347 arr = Tracking(4,nup-2,nup-2-gap,cuts,-1);
6348 SumTracks(seeds,arr);
6349 SignClusters(seeds,fnumber,fdensity);
6351 arr = Tracking(4,nup-3,nup-3-gap,cuts,-1);
6352 SumTracks(seeds,arr);
6353 SignClusters(seeds,fnumber,fdensity);
6357 for (Int_t delta = 3; delta<30; delta+=5){
6363 arr = Tracking(4,nup-1-delta,nup-1-delta-gap,cuts,-1);
6364 SumTracks(seeds,arr);
6365 SignClusters(seeds,fnumber,fdensity);
6367 arr = Tracking(4,nup-3-delta,nup-5-delta-gap,cuts,4);
6368 SumTracks(seeds,arr);
6369 SignClusters(seeds,fnumber,fdensity);
6380 Int_t fLastSeedRowSec=AliTPCReconstructor::GetRecoParam()->GetLastSeedRowSec();
6383 for (Int_t delta = 30; delta<fLastSeedRowSec; delta+=10){
6389 arr = Tracking(4,nup-1-delta,nup-1-delta-gap,cuts,-1);
6390 SumTracks(seeds,arr);
6391 SignClusters(seeds,fnumber,fdensity);
6393 arr = Tracking(4,nup-5-delta,nup-5-delta-gap,cuts,5 );
6394 SumTracks(seeds,arr);
6395 SignClusters(seeds,fnumber,fdensity);
6399 Info("Tracking()","\n\nSecondary seeding\t%d\n\n",seeds->GetEntriesFast());
6410 TObjArray * AliTPCtrackerMI::TrackingSpecial()
6413 // seeding adjusted for laser and cosmic tests - short tracks with big inclination angle
6414 // no primary vertex seeding tried
6418 Int_t nup=fOuterSec->GetNRows()+fInnerSec->GetNRows();
6420 TObjArray * seeds = new TObjArray;
6425 Float_t fnumber = 3.0;
6426 Float_t fdensity = 3.0;
6429 cuts[0] = AliTPCReconstructor::GetRecoParam()->GetMaxC(); // max curvature
6430 cuts[1] = 3.5; // max tan(phi) angle for seeding
6431 cuts[2] = 3.; // not used (cut on z primary vertex)
6432 cuts[3] = 3.5; // max tan(theta) angle for seeding
6434 for (Int_t delta = 0; nup-delta-gap-1>0; delta+=3){
6436 arr = Tracking(4,nup-1-delta,nup-1-delta-gap,cuts,-1);
6437 SumTracks(seeds,arr);
6438 SignClusters(seeds,fnumber,fdensity);
6442 Info("Tracking()","\n\nSecondary seeding\t%d\n\n",seeds->GetEntriesFast());
6453 void AliTPCtrackerMI::SumTracks(TObjArray *arr1,TObjArray *&arr2) const
6456 //sum tracks to common container
6457 //remove suspicious tracks
6458 Int_t nseed = arr2->GetEntriesFast();
6459 for (Int_t i=0;i<nseed;i++){
6460 AliTPCseed *pt=(AliTPCseed*)arr2->UncheckedAt(i);
6463 // remove tracks with too big curvature
6465 if (TMath::Abs(pt->GetC())>AliTPCReconstructor::GetRecoParam()->GetMaxC()){
6466 delete arr2->RemoveAt(i);
6469 // REMOVE VERY SHORT TRACKS
6470 if (pt->GetNumberOfClusters()<20){
6471 delete arr2->RemoveAt(i);
6474 // NORMAL ACTIVE TRACK
6475 if (pt->IsActive()){
6476 arr1->AddLast(arr2->RemoveAt(i));
6479 //remove not usable tracks
6480 if (pt->GetRemoval()!=10){
6481 delete arr2->RemoveAt(i);
6485 // ENABLE ONLY ENOUGH GOOD STOPPED TRACKS
6486 if (pt->GetDensityFirst(20)>0.8 || pt->GetDensityFirst(30)>0.8 || pt->GetDensityFirst(40)>0.7)
6487 arr1->AddLast(arr2->RemoveAt(i));
6489 delete arr2->RemoveAt(i);
6493 delete arr2; arr2 = 0;
6498 void AliTPCtrackerMI::ParallelTracking(TObjArray * arr, Int_t rfirst, Int_t rlast)
6501 // try to track in parralel
6503 Int_t nseed=arr->GetEntriesFast();
6504 //prepare seeds for tracking
6505 for (Int_t i=0; i<nseed; i++) {
6506 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i), &t=*pt;
6508 if (!t.IsActive()) continue;
6509 // follow prolongation to the first layer
6510 if ( (fSectors ==fInnerSec) || (t.GetFirstPoint()-fParam->GetNRowLow()>rfirst+1) )
6511 FollowProlongation(t, rfirst+1);
6516 for (Int_t nr=rfirst; nr>=rlast; nr--){
6517 if (nr<fInnerSec->GetNRows())
6518 fSectors = fInnerSec;
6520 fSectors = fOuterSec;
6521 // make indexes with the cluster tracks for given
6523 // find nearest cluster
6524 for (Int_t i=0; i<nseed; i++) {
6525 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i), &t=*pt;
6527 if (nr==80) pt->UpdateReference();
6528 if (!pt->IsActive()) continue;
6529 // if ( (fSectors ==fOuterSec) && (pt->fFirstPoint-fParam->GetNRowLow())<nr) continue;
6530 if (pt->GetRelativeSector()>17) {
6533 UpdateClusters(t,nr);
6535 // prolonagate to the nearest cluster - if founded
6536 for (Int_t i=0; i<nseed; i++) {
6537 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
6539 if (!pt->IsActive()) continue;
6540 // if ((fSectors ==fOuterSec) && (pt->fFirstPoint-fParam->GetNRowLow())<nr) continue;
6541 if (pt->GetRelativeSector()>17) {
6544 FollowToNextCluster(*pt,nr);
6549 void AliTPCtrackerMI::PrepareForBackProlongation(TObjArray * arr,Float_t fac) const
6553 // if we use TPC track itself we have to "update" covariance
6555 Int_t nseed= arr->GetEntriesFast();
6556 for (Int_t i=0;i<nseed;i++){
6557 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
6561 //rotate to current local system at first accepted point
6562 Int_t index = pt->GetClusterIndex2(pt->GetFirstPoint());
6563 Int_t sec = (index&0xff000000)>>24;
6565 Float_t angle1 = fInnerSec->GetAlpha()*sec+fInnerSec->GetAlphaShift();
6566 if (angle1>TMath::Pi())
6567 angle1-=2.*TMath::Pi();
6568 Float_t angle2 = pt->GetAlpha();
6570 if (TMath::Abs(angle1-angle2)>0.001){
6571 pt->Rotate(angle1-angle2);
6572 //angle2 = pt->GetAlpha();
6573 //pt->fRelativeSector = pt->GetAlpha()/fInnerSec->GetAlpha();
6574 //if (pt->GetAlpha()<0)
6575 // pt->fRelativeSector+=18;
6576 //sec = pt->fRelativeSector;
6585 void AliTPCtrackerMI::PrepareForProlongation(TObjArray * arr, Float_t fac) const
6589 // if we use TPC track itself we have to "update" covariance
6591 Int_t nseed= arr->GetEntriesFast();
6592 for (Int_t i=0;i<nseed;i++){
6593 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
6596 pt->SetFirstPoint(pt->GetLastPoint());
6604 Int_t AliTPCtrackerMI::PropagateBack(TObjArray * arr)
6607 // make back propagation
6609 Int_t nseed= arr->GetEntriesFast();
6610 for (Int_t i=0;i<nseed;i++){
6611 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
6612 if (pt&& pt->GetKinkIndex(0)<=0) {
6613 //AliTPCseed *pt2 = new AliTPCseed(*pt);
6614 fSectors = fInnerSec;
6615 //FollowBackProlongation(*pt,fInnerSec->GetNRows()-1);
6616 //fSectors = fOuterSec;
6617 FollowBackProlongation(*pt,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1);
6618 //if (pt->GetNumberOfClusters()<(pt->fEsd->GetTPCclusters(0)) ){
6619 // Error("PropagateBack","Not prolonged track %d",pt->GetLabel());
6620 // FollowBackProlongation(*pt2,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1);
6623 if (pt&& pt->GetKinkIndex(0)>0) {
6624 AliESDkink * kink = fEvent->GetKink(pt->GetKinkIndex(0)-1);
6625 pt->SetFirstPoint(kink->GetTPCRow0());
6626 fSectors = fInnerSec;
6627 FollowBackProlongation(*pt,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1);
6635 Int_t AliTPCtrackerMI::PropagateForward2(TObjArray * arr)
6638 // make forward propagation
6640 Int_t nseed= arr->GetEntriesFast();
6642 for (Int_t i=0;i<nseed;i++){
6643 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
6645 FollowProlongation(*pt,0);
6654 Int_t AliTPCtrackerMI::PropagateForward()
6657 // propagate track forward
6659 Int_t nseed = fSeeds->GetEntriesFast();
6660 for (Int_t i=0;i<nseed;i++){
6661 AliTPCseed *pt = (AliTPCseed*)fSeeds->UncheckedAt(i);
6663 AliTPCseed &t = *pt;
6664 Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
6665 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
6666 if (alpha < 0. ) alpha += 2.*TMath::Pi();
6667 t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN);
6671 fSectors = fOuterSec;
6672 ParallelTracking(fSeeds,fOuterSec->GetNRows()+fInnerSec->GetNRows()-1,fInnerSec->GetNRows());
6673 fSectors = fInnerSec;
6674 ParallelTracking(fSeeds,fInnerSec->GetNRows()-1,0);
6684 Int_t AliTPCtrackerMI::PropagateBack(AliTPCseed * pt, Int_t row0, Int_t row1)
6687 // make back propagation, in between row0 and row1
6691 fSectors = fInnerSec;
6694 if (row1<fSectors->GetNRows())
6697 r1 = fSectors->GetNRows()-1;
6699 if (row0<fSectors->GetNRows()&& r1>0 )
6700 FollowBackProlongation(*pt,r1);
6701 if (row1<=fSectors->GetNRows())
6704 r1 = row1 - fSectors->GetNRows();
6705 if (r1<=0) return 0;
6706 if (r1>=fOuterSec->GetNRows()) return 0;
6707 fSectors = fOuterSec;
6708 return FollowBackProlongation(*pt,r1);
6716 void AliTPCtrackerMI::GetShape(AliTPCseed * seed, Int_t row)
6720 Float_t sd2 = TMath::Abs((fParam->GetZLength(0)-TMath::Abs(seed->GetZ())))*fParam->GetDiffL()*fParam->GetDiffL();
6721 // Float_t padlength = fParam->GetPadPitchLength(seed->fSector);
6722 Float_t padlength = GetPadPitchLength(row);
6724 Float_t sresy = (seed->GetSector() < fParam->GetNSector()/2) ? 0.2 :0.3;
6725 Double_t angulary = seed->GetSnp();
6726 angulary = angulary*angulary/(1-angulary*angulary);
6727 seed->SetCurrentSigmaY2(sd2+padlength*padlength*angulary/12.+sresy*sresy);
6729 Float_t sresz = fParam->GetZSigma();
6730 Float_t angularz = seed->GetTgl();
6731 seed->SetCurrentSigmaZ2(sd2+padlength*padlength*angularz*angularz*(1+angulary)/12.+sresz*sresz);
6733 Float_t wy = GetSigmaY(seed);
6734 Float_t wz = GetSigmaZ(seed);
6737 if (TMath::Abs(wy/seed->fCurrentSigmaY2-1)>0.0001 || TMath::Abs(wz/seed->fCurrentSigmaZ2-1)>0.0001 ){
6738 printf("problem\n");
6744 Float_t AliTPCtrackerMI::GetSigmaY(AliTPCseed * seed)
6748 Float_t sd2 = TMath::Abs((fParam->GetZLength(0)-TMath::Abs(seed->GetZ())))*fParam->GetDiffL()*fParam->GetDiffL();
6749 Float_t padlength = fParam->GetPadPitchLength(seed->GetSector());
6750 Float_t sres = (seed->GetSector() < fParam->GetNSector()/2) ? 0.2 :0.3;
6751 Float_t angular = seed->GetSnp();
6752 angular = angular*angular/(1-angular*angular);
6753 // angular*=angular;
6754 //angular = TMath::Sqrt(angular/(1-angular));
6755 Float_t res = TMath::Sqrt(sd2+padlength*padlength*angular/12.+sres*sres);
6758 Float_t AliTPCtrackerMI::GetSigmaZ(AliTPCseed * seed)
6762 Float_t sd2 = TMath::Abs((fParam->GetZLength(0)-TMath::Abs(seed->GetZ())))*fParam->GetDiffL()*fParam->GetDiffL();
6763 Float_t padlength = fParam->GetPadPitchLength(seed->GetSector());
6764 Float_t sres = fParam->GetZSigma();
6765 Float_t angular = seed->GetTgl();
6766 Float_t res = TMath::Sqrt(sd2+padlength*padlength*angular*angular/12.+sres*sres);
6771 //__________________________________________________________________________
6772 void AliTPCtrackerMI::CookLabel(AliKalmanTrack *tk, Float_t wrong) const {
6773 //--------------------------------------------------------------------
6774 //This function "cooks" a track label. If label<0, this track is fake.
6775 //--------------------------------------------------------------------
6776 AliTPCseed * t = dynamic_cast<AliTPCseed*>(tk);
6778 printf("%s:%d wrong type \n",(char*)__FILE__,__LINE__);
6782 Int_t noc=t->GetNumberOfClusters();
6784 //printf("\nnot founded prolongation\n\n\n");
6790 AliTPCclusterMI *clusters[160];
6792 for (Int_t i=0;i<160;i++) {
6799 for (i=0; i<160 && current<noc; i++) {
6801 Int_t index=t->GetClusterIndex2(i);
6802 if (index<=0) continue;
6803 if (index&0x8000) continue;
6805 //clusters[current]=GetClusterMI(index);
6806 if (t->GetClusterPointer(i)){
6807 clusters[current]=t->GetClusterPointer(i);
6813 Int_t lab=123456789;
6814 for (i=0; i<noc; i++) {
6815 AliTPCclusterMI *c=clusters[i];
6817 lab=TMath::Abs(c->GetLabel(0));
6819 for (j=0; j<noc; j++) if (lb[j]==lab || mx[j]==0) break;
6825 for (i=0; i<noc; i++) if (mx[i]>max) {max=mx[i]; lab=lb[i];}
6827 for (i=0; i<noc; i++) {
6828 AliTPCclusterMI *c=clusters[i];
6830 if (TMath::Abs(c->GetLabel(1)) == lab ||
6831 TMath::Abs(c->GetLabel(2)) == lab ) max++;
6834 if ((1.- Float_t(max)/noc) > wrong) lab=-lab;
6837 Int_t tail=Int_t(0.10*noc);
6840 for (i=1; i<=160&&ind<tail; i++) {
6841 // AliTPCclusterMI *c=clusters[noc-i];
6842 AliTPCclusterMI *c=clusters[i];
6844 if (lab == TMath::Abs(c->GetLabel(0)) ||
6845 lab == TMath::Abs(c->GetLabel(1)) ||
6846 lab == TMath::Abs(c->GetLabel(2))) max++;
6849 if (max < Int_t(0.5*tail)) lab=-lab;
6856 //delete[] clusters;
6860 //__________________________________________________________________________
6861 Int_t AliTPCtrackerMI::CookLabel(AliTPCseed *t, Float_t wrong,Int_t first, Int_t last) const {
6862 //--------------------------------------------------------------------
6863 //This function "cooks" a track label. If label<0, this track is fake.
6864 //--------------------------------------------------------------------
6865 Int_t noc=t->GetNumberOfClusters();
6867 //printf("\nnot founded prolongation\n\n\n");
6873 AliTPCclusterMI *clusters[160];
6875 for (Int_t i=0;i<160;i++) {
6882 for (i=0; i<160 && current<noc; i++) {
6883 if (i<first) continue;
6884 if (i>last) continue;
6885 Int_t index=t->GetClusterIndex2(i);
6886 if (index<=0) continue;
6887 if (index&0x8000) continue;
6889 //clusters[current]=GetClusterMI(index);
6890 if (t->GetClusterPointer(i)){
6891 clusters[current]=t->GetClusterPointer(i);
6896 if (noc<5) return -1;
6897 Int_t lab=123456789;
6898 for (i=0; i<noc; i++) {
6899 AliTPCclusterMI *c=clusters[i];
6901 lab=TMath::Abs(c->GetLabel(0));
6903 for (j=0; j<noc; j++) if (lb[j]==lab || mx[j]==0) break;
6909 for (i=0; i<noc; i++) if (mx[i]>max) {max=mx[i]; lab=lb[i];}
6911 for (i=0; i<noc; i++) {
6912 AliTPCclusterMI *c=clusters[i];
6914 if (TMath::Abs(c->GetLabel(1)) == lab ||
6915 TMath::Abs(c->GetLabel(2)) == lab ) max++;
6918 if ((1.- Float_t(max)/noc) > wrong) lab=-lab;
6921 Int_t tail=Int_t(0.10*noc);
6924 for (i=1; i<=160&&ind<tail; i++) {
6925 // AliTPCclusterMI *c=clusters[noc-i];
6926 AliTPCclusterMI *c=clusters[i];
6928 if (lab == TMath::Abs(c->GetLabel(0)) ||
6929 lab == TMath::Abs(c->GetLabel(1)) ||
6930 lab == TMath::Abs(c->GetLabel(2))) max++;
6933 if (max < Int_t(0.5*tail)) lab=-lab;
6936 // t->SetLabel(lab);
6940 //delete[] clusters;
6944 Int_t AliTPCtrackerMI::AliTPCSector::GetRowNumber(Double_t x) const
6946 //return pad row number for this x
6949 r=fRow[fN-1].GetX();
6950 if (x > r) return fN;
6952 if (x < r) return -1;
6953 return Int_t((x-r)/fPadPitchLength + 0.5);}
6955 r=fRow[fN-1].GetX();
6956 if (x > r) return fN;
6958 if (x < r) return -1;
6959 Double_t r1=fRow[64].GetX();
6961 return Int_t((x-r)/f1PadPitchLength + 0.5);}
6963 return (Int_t((x-r1)/f2PadPitchLength + 0.5)+64);}
6967 Int_t AliTPCtrackerMI::GetRowNumber(Double_t x[3]) const
6969 //return pad row number for given x vector
6970 Float_t phi = TMath::ATan2(x[1],x[0]);
6971 if(phi<0) phi=2.*TMath::Pi()+phi;
6972 // Get the local angle in the sector philoc
6973 const Float_t kRaddeg = 180/3.14159265358979312;
6974 Float_t phiangle = (Int_t (phi*kRaddeg/20.) + 0.5)*20./kRaddeg;
6975 Double_t localx = x[0]*TMath::Cos(phiangle)-x[1]*TMath::Sin(phiangle);
6976 return GetRowNumber(localx);
6979 //_________________________________________________________________________
6980 void AliTPCtrackerMI::AliTPCSector::Setup(const AliTPCParam *par, Int_t f) {
6981 //-----------------------------------------------------------------------
6982 // Setup inner sector
6983 //-----------------------------------------------------------------------
6985 fAlpha=par->GetInnerAngle();
6986 fAlphaShift=par->GetInnerAngleShift();
6987 fPadPitchWidth=par->GetInnerPadPitchWidth();
6988 fPadPitchLength=par->GetInnerPadPitchLength();
6989 fN=par->GetNRowLow();
6990 if(fRow)delete [] fRow;fRow = 0;
6991 fRow=new AliTPCRow[fN];
6992 for (Int_t i=0; i<fN; i++) {
6993 fRow[i].SetX(par->GetPadRowRadiiLow(i));
6994 fRow[i].SetDeadZone(1.5); //1.5 cm of dead zone
6997 fAlpha=par->GetOuterAngle();
6998 fAlphaShift=par->GetOuterAngleShift();
6999 fPadPitchWidth = par->GetOuterPadPitchWidth();
7000 fPadPitchLength = par->GetOuter1PadPitchLength();
7001 f1PadPitchLength = par->GetOuter1PadPitchLength();
7002 f2PadPitchLength = par->GetOuter2PadPitchLength();
7003 fN=par->GetNRowUp();
7004 if(fRow)delete [] fRow;fRow = 0;
7005 fRow=new AliTPCRow[fN];
7006 for (Int_t i=0; i<fN; i++) {
7007 fRow[i].SetX(par->GetPadRowRadiiUp(i));
7008 fRow[i].SetDeadZone(1.5); // 1.5 cm of dead zone
7013 AliTPCtrackerMI::AliTPCRow::AliTPCRow():
7023 // default constructor
7027 AliTPCtrackerMI::AliTPCRow::~AliTPCRow(){
7029 for (Int_t i = 0; i < fN1; i++)
7030 fClusters1[i].~AliTPCclusterMI();
7031 delete [] fClusters1;
7032 for (Int_t i = 0; i < fN2; i++)
7033 fClusters2[i].~AliTPCclusterMI();
7034 delete [] fClusters2;
7039 //_________________________________________________________________________
7041 AliTPCtrackerMI::AliTPCRow::InsertCluster(const AliTPCclusterMI* c, UInt_t index) {
7042 //-----------------------------------------------------------------------
7043 // Insert a cluster into this pad row in accordence with its y-coordinate
7044 //-----------------------------------------------------------------------
7045 if (fN==kMaxClusterPerRow) {
7046 cerr<<"AliTPCRow::InsertCluster(): Too many clusters !\n"; return;
7048 if (fN==0) {fIndex[0]=index; fClusters[fN++]=c; return;}
7049 Int_t i=Find(c->GetZ());
7050 memmove(fClusters+i+1 ,fClusters+i,(fN-i)*sizeof(AliTPCclusterMI*));
7051 memmove(fIndex +i+1 ,fIndex +i,(fN-i)*sizeof(UInt_t));
7052 fIndex[i]=index; fClusters[i]=c; fN++;
7055 void AliTPCtrackerMI::AliTPCRow::ResetClusters() {
7058 // MvL: Need to call destructors for AliTPCclusterMI, to delete fInfo
7059 for (Int_t i = 0; i < fN1; i++)
7060 fClusters1[i].~AliTPCclusterMI();
7061 delete [] fClusters1; fClusters1=0;
7062 for (Int_t i = 0; i < fN2; i++)
7063 fClusters2[i].~AliTPCclusterMI();
7064 delete [] fClusters2; fClusters2=0;
7069 //delete[] fClusterArray;
7075 //___________________________________________________________________
7076 Int_t AliTPCtrackerMI::AliTPCRow::Find(Double_t z) const {
7077 //-----------------------------------------------------------------------
7078 // Return the index of the nearest cluster
7079 //-----------------------------------------------------------------------
7080 if (fN==0) return 0;
7081 if (z <= fClusters[0]->GetZ()) return 0;
7082 if (z > fClusters[fN-1]->GetZ()) return fN;
7083 Int_t b=0, e=fN-1, m=(b+e)/2;
7084 for (; b<e; m=(b+e)/2) {
7085 if (z > fClusters[m]->GetZ()) b=m+1;
7093 //___________________________________________________________________
7094 AliTPCclusterMI * AliTPCtrackerMI::AliTPCRow::FindNearest(Double_t y, Double_t z, Double_t roady, Double_t roadz) const {
7095 //-----------------------------------------------------------------------
7096 // Return the index of the nearest cluster in z y
7097 //-----------------------------------------------------------------------
7098 Float_t maxdistance = roady*roady + roadz*roadz;
7100 AliTPCclusterMI *cl =0;
7101 for (Int_t i=Find(z-roadz); i<fN; i++) {
7102 AliTPCclusterMI *c=(AliTPCclusterMI*)(fClusters[i]);
7103 if (c->GetZ() > z+roadz) break;
7104 if ( (c->GetY()-y) > roady ) continue;
7105 Float_t distance = (c->GetZ()-z)*(c->GetZ()-z)+(c->GetY()-y)*(c->GetY()-y);
7106 if (maxdistance>distance) {
7107 maxdistance = distance;
7114 AliTPCclusterMI * AliTPCtrackerMI::AliTPCRow::FindNearest2(Double_t y, Double_t z, Double_t roady, Double_t roadz,UInt_t & index) const
7116 //-----------------------------------------------------------------------
7117 // Return the index of the nearest cluster in z y
7118 //-----------------------------------------------------------------------
7119 Float_t maxdistance = roady*roady + roadz*roadz;
7120 AliTPCclusterMI *cl =0;
7122 //PH Check boundaries. 510 is the size of fFastCluster
7123 Int_t iz1 = Int_t(z-roadz+254.5);
7124 if (iz1<0 || iz1>=510) return cl;
7125 iz1 = TMath::Max(GetFastCluster(iz1)-1,0);
7126 Int_t iz2 = Int_t(z+roadz+255.5);
7127 if (iz2<0 || iz2>=510) return cl;
7128 iz2 = TMath::Min(GetFastCluster(iz2)+1,fN);
7129 Bool_t skipUsed = !(AliTPCReconstructor::GetRecoParam()->GetClusterSharing());
7130 //FindNearest3(y,z,roady,roadz,index);
7131 // for (Int_t i=Find(z-roadz); i<fN; i++) {
7132 for (Int_t i=iz1; i<iz2; i++) {
7133 AliTPCclusterMI *c=(AliTPCclusterMI*)(fClusters[i]);
7134 if (c->GetZ() > z+roadz) break;
7135 if ( c->GetY()-y > roady ) continue;
7136 if ( y-c->GetY() > roady ) continue;
7137 if (skipUsed && c->IsUsed(11)) continue;
7138 Float_t distance = (c->GetZ()-z)*(c->GetZ()-z)+(c->GetY()-y)*(c->GetY()-y);
7139 if (maxdistance>distance) {
7140 maxdistance = distance;
7143 //roady = TMath::Sqrt(maxdistance);
7153 void AliTPCtrackerMI::MakeBitmaps(AliTPCseed *t)
7155 //-----------------------------------------------------------------------
7156 // Fill the cluster and sharing bitmaps of the track
7157 //-----------------------------------------------------------------------
7159 Int_t firstpoint = 0;
7160 Int_t lastpoint = 159;
7161 AliTPCTrackerPoint *point;
7163 for (int iter=firstpoint; iter<lastpoint; iter++) {
7164 point = t->GetTrackPoint(iter);
7166 t->SetClusterMapBit(iter, kTRUE);
7167 if (point->IsShared())
7168 t->SetSharedMapBit(iter,kTRUE);
7170 t->SetSharedMapBit(iter, kFALSE);
7173 t->SetClusterMapBit(iter, kFALSE);
7174 t->SetSharedMapBit(iter, kFALSE);