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 // run AliTPCFindClusters.C macro - clusters neccessary for tracker are founded
25 // run AliTPCFindTracksMI.C macro - to find tracks
26 // tracks are written to AliTPCtracks.root file
27 // for comparison also seeds are written to the same file - to special branch
28 //-------------------------------------------------------
35 #include "Riostream.h"
36 #include <TClonesArray.h>
38 #include <TObjArray.h>
41 #include "AliComplexCluster.h"
44 #include "AliRunLoader.h"
45 #include "AliTPCClustersRow.h"
46 #include "AliTPCParam.h"
47 #include "AliTPCclusterMI.h"
48 #include "AliTPCpolyTrack.h"
49 #include "AliTPCreco.h"
50 #include "AliTPCtrackerMI.h"
51 #include "TStopwatch.h"
53 #include "AliTPCReconstructor.h"
54 #include "AliESDkink.h"
55 #include "TTreeStream.h"
59 ClassImp(AliTPCtrackerMI)
62 class AliTPCFastMath {
65 static Double_t FastAsin(Double_t x);
67 static Double_t fgFastAsin[20000]; //lookup table for fast asin computation
70 Double_t AliTPCFastMath::fgFastAsin[20000];
71 AliTPCFastMath gAliTPCFastMath; // needed to fill the LUT
73 AliTPCFastMath::AliTPCFastMath(){
75 // initialized lookup table;
76 for (Int_t i=0;i<10000;i++){
77 fgFastAsin[2*i] = TMath::ASin(i/10000.);
78 fgFastAsin[2*i+1] = (TMath::ASin((i+1)/10000.)-fgFastAsin[2*i]);
82 Double_t AliTPCFastMath::FastAsin(Double_t x){
84 // return asin using lookup table
86 Int_t index = int(x*10000);
87 return fgFastAsin[2*index]+(x*10000.-index)*fgFastAsin[2*index+1];
90 Int_t index = int(x*10000);
91 return -(fgFastAsin[2*index]+(x*10000.-index)*fgFastAsin[2*index+1]);
97 Int_t AliTPCtrackerMI::UpdateTrack(AliTPCseed * track, Int_t accept){
99 //update track information using current cluster - track->fCurrentCluster
102 AliTPCclusterMI* c =track->fCurrentCluster;
103 if (accept>0) track->fCurrentClusterIndex1 |=0x8000; //sign not accepted clusters
105 UInt_t i = track->fCurrentClusterIndex1;
107 Int_t sec=(i&0xff000000)>>24;
108 //Int_t row = (i&0x00ff0000)>>16;
109 track->fRow=(i&0x00ff0000)>>16;
110 track->fSector = sec;
111 // Int_t index = i&0xFFFF;
112 if (sec>=fParam->GetNInnerSector()) track->fRow += fParam->GetNRowLow();
113 track->SetClusterIndex2(track->fRow, i);
114 //track->fFirstPoint = row;
115 //if ( track->fLastPoint<row) track->fLastPoint =row;
116 // if (track->fRow<0 || track->fRow>160) {
117 // printf("problem\n");
119 if (track->fFirstPoint>track->fRow)
120 track->fFirstPoint = track->fRow;
121 if (track->fLastPoint<track->fRow)
122 track->fLastPoint = track->fRow;
125 track->fClusterPointer[track->fRow] = c;
128 Float_t angle2 = track->GetSnp()*track->GetSnp();
129 angle2 = TMath::Sqrt(angle2/(1-angle2));
131 //SET NEW Track Point
135 AliTPCTrackerPoint &point =*(track->GetTrackPoint(track->fRow));
137 point.SetSigmaY(c->GetSigmaY2()/track->fCurrentSigmaY2);
138 point.SetSigmaZ(c->GetSigmaZ2()/track->fCurrentSigmaZ2);
139 point.SetErrY(sqrt(track->fErrorY2));
140 point.SetErrZ(sqrt(track->fErrorZ2));
142 point.SetX(track->GetX());
143 point.SetY(track->GetY());
144 point.SetZ(track->GetZ());
145 point.SetAngleY(angle2);
146 point.SetAngleZ(track->GetTgl());
147 if (point.fIsShared){
148 track->fErrorY2 *= 4;
149 track->fErrorZ2 *= 4;
153 Double_t chi2 = track->GetPredictedChi2(track->fCurrentCluster);
155 track->fErrorY2 *= 1.3;
156 track->fErrorY2 += 0.01;
157 track->fErrorZ2 *= 1.3;
158 track->fErrorZ2 += 0.005;
160 if (accept>0) return 0;
161 if (track->GetNumberOfClusters()%20==0){
162 // if (track->fHelixIn){
163 // TClonesArray & larr = *(track->fHelixIn);
164 // Int_t ihelix = larr.GetEntriesFast();
165 // new(larr[ihelix]) AliHelix(*track) ;
168 track->fNoCluster =0;
169 return track->Update(c,chi2,i);
174 Int_t AliTPCtrackerMI::AcceptCluster(AliTPCseed * seed, AliTPCclusterMI * cluster, Float_t factor,
175 Float_t cory, Float_t corz)
178 // decide according desired precision to accept given
179 // cluster for tracking
180 Double_t sy2=ErrY2(seed,cluster)*cory;
181 Double_t sz2=ErrZ2(seed,cluster)*corz;
182 //sy2=ErrY2(seed,cluster)*cory;
183 //sz2=ErrZ2(seed,cluster)*cory;
185 Double_t sdistancey2 = sy2+seed->GetSigmaY2();
186 Double_t sdistancez2 = sz2+seed->GetSigmaZ2();
188 Double_t rdistancey2 = (seed->fCurrentCluster->GetY()-seed->GetY())*
189 (seed->fCurrentCluster->GetY()-seed->GetY())/sdistancey2;
190 Double_t rdistancez2 = (seed->fCurrentCluster->GetZ()-seed->GetZ())*
191 (seed->fCurrentCluster->GetZ()-seed->GetZ())/sdistancez2;
193 Double_t rdistance2 = rdistancey2+rdistancez2;
196 if (rdistance2>16) return 3;
199 if ((rdistancey2>9.*factor || rdistancez2>9.*factor) && cluster->GetType()==0)
200 return 2; //suspisiouce - will be changed
202 if ((rdistancey2>6.25*factor || rdistancez2>6.25*factor) && cluster->GetType()>0)
203 // strict cut on overlaped cluster
204 return 2; //suspisiouce - will be changed
206 if ( (rdistancey2>1.*factor || rdistancez2>6.25*factor )
207 && cluster->GetType()<0){
217 //_____________________________________________________________________________
218 AliTPCtrackerMI::AliTPCtrackerMI(const AliTPCParam *par):
219 AliTracker(), fkNIS(par->GetNInnerSector()/2), fkNOS(par->GetNOuterSector()/2)
221 //---------------------------------------------------------------------
222 // The main TPC tracker constructor
223 //---------------------------------------------------------------------
224 fInnerSec=new AliTPCSector[fkNIS];
225 fOuterSec=new AliTPCSector[fkNOS];
228 for (i=0; i<fkNIS; i++) fInnerSec[i].Setup(par,0);
229 for (i=0; i<fkNOS; i++) fOuterSec[i].Setup(par,1);
236 Int_t nrowlow = par->GetNRowLow();
237 Int_t nrowup = par->GetNRowUp();
240 for (Int_t i=0;i<nrowlow;i++){
241 fXRow[i] = par->GetPadRowRadiiLow(i);
242 fPadLength[i]= par->GetPadPitchLength(0,i);
243 fYMax[i] = fXRow[i]*TMath::Tan(0.5*par->GetInnerAngle());
247 for (Int_t i=0;i<nrowup;i++){
248 fXRow[i+nrowlow] = par->GetPadRowRadiiUp(i);
249 fPadLength[i+nrowlow] = par->GetPadPitchLength(60,i);
250 fYMax[i+nrowlow] = fXRow[i+nrowlow]*TMath::Tan(0.5*par->GetOuterAngle());
262 //________________________________________________________________________
263 AliTPCtrackerMI::AliTPCtrackerMI(const AliTPCtrackerMI &t):
268 //------------------------------------
269 // dummy copy constructor
270 //------------------------------------------------------------------
272 AliTPCtrackerMI & AliTPCtrackerMI::operator=(const AliTPCtrackerMI& /*r*/){
273 //------------------------------
275 //--------------------------------------------------------------
278 //_____________________________________________________________________________
279 AliTPCtrackerMI::~AliTPCtrackerMI() {
280 //------------------------------------------------------------------
281 // TPC tracker destructor
282 //------------------------------------------------------------------
291 void AliTPCtrackerMI::SetIO()
295 fInput = AliRunLoader::GetTreeR("TPC", kFALSE,AliConfig::GetDefaultEventFolderName());
297 fOutput = AliRunLoader::GetTreeT("TPC", kTRUE,AliConfig::GetDefaultEventFolderName());
299 AliTPCtrack *iotrack= new AliTPCtrack;
300 fOutput->Branch("tracks","AliTPCtrack",&iotrack,32000,100);
306 void AliTPCtrackerMI::SetIO(TTree * input, TTree * output, AliESD * event)
322 AliTPCtrack *iotrack= new AliTPCtrack;
323 // iotrack->fHelixIn = new TClonesArray("AliHelix");
324 //iotrack->fHelixOut = new TClonesArray("AliHelix");
325 fOutput->Branch("tracks","AliTPCtrack",&iotrack,32000,100);
328 if (output && (fDebug&2)){
329 //write the full seed information if specified in debug mode
331 fSeedTree = new TTree("Seeds","Seeds");
332 AliTPCseed * vseed = new AliTPCseed;
334 TClonesArray * arrtr = new TClonesArray("AliTPCTrackPoint",160);
335 arrtr->ExpandCreateFast(160);
336 TClonesArray * arre = new TClonesArray("AliTPCExactPoint",160);
338 vseed->fPoints = arrtr;
339 vseed->fEPoints = arre;
340 // vseed->fClusterPoints = arrcl;
341 fSeedTree->Branch("seeds","AliTPCseed",&vseed,32000,99);
344 fTreeDebug = new TTree("trackDebug","trackDebug");
345 TClonesArray * arrd = new TClonesArray("AliTPCTrackPoint2",0);
346 fTreeDebug->Branch("debug",&arrd,32000,99);
354 void AliTPCtrackerMI::FillESD(TObjArray* arr)
358 //fill esds using updated tracks
360 // write tracks to the event
361 // store index of the track
362 Int_t nseed=arr->GetEntriesFast();
363 //FindKinks(arr,fEvent);
364 for (Int_t i=0; i<nseed; i++) {
365 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
368 // pt->PropagateTo(fParam->GetInnerRadiusLow());
369 if (pt->GetKinkIndex(0)<=0){ //don't propagate daughter tracks
370 pt->PropagateTo(fParam->GetInnerRadiusLow());
373 if (( pt->GetPoints()[2]- pt->GetPoints()[0])>5 && pt->GetPoints()[3]>0.8){
375 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
376 iotrack.SetTPCPoints(pt->GetPoints());
377 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
378 //iotrack.SetTPCindex(i);
379 fEvent->AddTrack(&iotrack);
383 if ( (pt->GetNumberOfClusters()>70)&& (Float_t(pt->GetNumberOfClusters())/Float_t(pt->fNFoundable))>0.55) {
385 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
386 iotrack.SetTPCPoints(pt->GetPoints());
387 //iotrack.SetTPCindex(i);
388 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
389 fEvent->AddTrack(&iotrack);
393 // short tracks - maybe decays
395 if ( (pt->GetNumberOfClusters()>30) && (Float_t(pt->GetNumberOfClusters())/Float_t(pt->fNFoundable))>0.70) {
396 Int_t found,foundable,shared;
397 pt->GetClusterStatistic(0,60,found, foundable,shared,kFALSE);
398 if ( (found>20) && (pt->fNShared/float(pt->GetNumberOfClusters())<0.2)){
400 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
401 //iotrack.SetTPCindex(i);
402 iotrack.SetTPCPoints(pt->GetPoints());
403 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
404 fEvent->AddTrack(&iotrack);
409 if ( (pt->GetNumberOfClusters()>20) && (Float_t(pt->GetNumberOfClusters())/Float_t(pt->fNFoundable))>0.8) {
410 Int_t found,foundable,shared;
411 pt->GetClusterStatistic(0,60,found, foundable,shared,kFALSE);
412 if (found<20) continue;
413 if (pt->fNShared/float(pt->GetNumberOfClusters())>0.2) continue;
416 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
417 iotrack.SetTPCPoints(pt->GetPoints());
418 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
419 //iotrack.SetTPCindex(i);
420 fEvent->AddTrack(&iotrack);
423 // short tracks - secondaties
425 if ( (pt->GetNumberOfClusters()>30) ) {
426 Int_t found,foundable,shared;
427 pt->GetClusterStatistic(128,158,found, foundable,shared,kFALSE);
428 if ( (found>20) && (pt->fNShared/float(pt->GetNumberOfClusters())<0.2) &&float(found)/float(foundable)>0.8){
430 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
431 iotrack.SetTPCPoints(pt->GetPoints());
432 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
433 //iotrack.SetTPCindex(i);
434 fEvent->AddTrack(&iotrack);
439 if ( (pt->GetNumberOfClusters()>15)) {
440 Int_t found,foundable,shared;
441 pt->GetClusterStatistic(138,158,found, foundable,shared,kFALSE);
442 if (found<15) continue;
443 if (pt->fNShared/float(pt->GetNumberOfClusters())>0.2) continue;
444 if (float(found)/float(foundable)<0.8) continue;
447 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
448 iotrack.SetTPCPoints(pt->GetPoints());
449 iotrack.SetKinkIndexes(pt->GetKinkIndexes());
450 //iotrack.SetTPCindex(i);
451 fEvent->AddTrack(&iotrack);
456 printf("Number of filled ESDs-\t%d\n",fEvent->GetNumberOfTracks());
459 void AliTPCtrackerMI::WriteTracks(TTree * tree)
462 // write tracks from seed array to selected tree
466 AliTPCtrack *iotrack= new AliTPCtrack;
467 fOutput->Branch("tracks","AliTPCtrack",&iotrack,32000,100);
472 void AliTPCtrackerMI::WriteTracks()
475 // write tracks to the given output tree -
476 // output specified with SetIO routine
483 AliTPCtrack *iotrack= 0;
484 Int_t nseed=fSeeds->GetEntriesFast();
485 //for (Int_t i=0; i<nseed; i++) {
486 // iotrack= (AliTPCtrack*)fSeeds->UncheckedAt(i);
487 // if (iotrack) break;
489 //TBranch * br = fOutput->Branch("tracks","AliTPCtrack",&iotrack,32000,100);
490 TBranch * br = fOutput->GetBranch("tracks");
491 br->SetAddress(&iotrack);
493 for (Int_t i=0; i<nseed; i++) {
494 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i);
496 AliTPCtrack * track = new AliTPCtrack(*pt);
499 // br->SetAddress(&iotrack);
504 //fOutput->GetDirectory()->cd();
510 //write the full seed information if specified in debug mode
512 AliTPCseed * vseed = new AliTPCseed;
514 TClonesArray * arrtr = new TClonesArray("AliTPCTrackPoint",160);
515 arrtr->ExpandCreateFast(160);
516 //TClonesArray * arrcl = new TClonesArray("AliTPCclusterMI",160);
517 //arrcl->ExpandCreateFast(160);
518 TClonesArray * arre = new TClonesArray("AliTPCExactPoint",160);
520 vseed->fPoints = arrtr;
521 vseed->fEPoints = arre;
522 // vseed->fClusterPoints = arrcl;
523 //TBranch * brseed = seedtree->Branch("seeds","AliTPCseed",&vseed,32000,99);
524 TBranch * brseed = fSeedTree->GetBranch("seeds");
526 Int_t nseed=fSeeds->GetEntriesFast();
528 for (Int_t i=0; i<nseed; i++) {
529 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i);
532 // pt->fClusterPoints = arrcl;
536 brseed->SetAddress(&vseed);
540 // pt->fClusterPoints = 0;
543 if (fTreeDebug) fTreeDebug->Write();
551 Double_t AliTPCtrackerMI::ErrY2(AliTPCseed* seed, AliTPCclusterMI * cl){
554 //seed->SetErrorY2(0.1);
556 //calculate look-up table at the beginning
557 static Bool_t ginit = kFALSE;
558 static Float_t gnoise1,gnoise2,gnoise3;
559 static Float_t ggg1[10000];
560 static Float_t ggg2[10000];
561 static Float_t ggg3[10000];
562 static Float_t glandau1[10000];
563 static Float_t glandau2[10000];
564 static Float_t glandau3[10000];
566 static Float_t gcor01[500];
567 static Float_t gcor02[500];
568 static Float_t gcorp[500];
573 for (Int_t i=1;i<500;i++){
574 Float_t rsigma = float(i)/100.;
575 gcor02[i] = TMath::Max(0.78 +TMath::Exp(7.4*(rsigma-1.2)),0.6);
576 gcor01[i] = TMath::Max(0.72 +TMath::Exp(3.36*(rsigma-1.2)),0.6);
577 gcorp[i] = TMath::Max(TMath::Power((rsigma+0.5),1.5),1.2);
581 for (Int_t i=3;i<10000;i++){
585 Float_t amp = float(i);
586 Float_t padlength =0.75;
587 gnoise1 = 0.0004/padlength;
588 Float_t nel = 0.268*amp;
589 Float_t nprim = 0.155*amp;
590 ggg1[i] = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.001*nel/(padlength*padlength))/nel;
591 glandau1[i] = (2.+0.12*nprim)*0.5* (2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
592 if (glandau1[i]>1) glandau1[i]=1;
593 glandau1[i]*=padlength*padlength/12.;
597 gnoise2 = 0.0004/padlength;
600 ggg2[i] = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
601 glandau2[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
602 if (glandau2[i]>1) glandau2[i]=1;
603 glandau2[i]*=padlength*padlength/12.;
608 gnoise3 = 0.0004/padlength;
611 ggg3[i] = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
612 glandau3[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
613 if (glandau3[i]>1) glandau3[i]=1;
614 glandau3[i]*=padlength*padlength/12.;
622 Int_t amp = int(TMath::Abs(cl->GetQ()));
624 seed->SetErrorY2(1.);
628 Float_t z = TMath::Abs(fParam->GetZLength()-TMath::Abs(seed->GetZ()));
629 Int_t ctype = cl->GetType();
630 Float_t padlength= GetPadPitchLength(seed->fRow);
631 Float_t angle2 = seed->GetSnp()*seed->GetSnp();
632 angle2 = angle2/(1-angle2);
635 Int_t rsigmay = int(100.*cl->GetSigmaY2()/(seed->fCurrentSigmaY2));
638 if (fSectors==fInnerSec){
640 res = ggg1[amp]*z+glandau1[amp]*angle2;
641 if (ctype==0) res *= gcor01[rsigmay];
644 res*= gcorp[rsigmay];
650 res = ggg2[amp]*z+glandau2[amp]*angle2;
651 if (ctype==0) res *= gcor02[rsigmay];
654 res*= gcorp[rsigmay];
659 res = ggg3[amp]*z+glandau3[amp]*angle2;
660 if (ctype==0) res *= gcor02[rsigmay];
663 res*= gcorp[rsigmay];
670 res*=2.4; // overestimate error 2 times
677 seed->SetErrorY2(res);
685 Double_t AliTPCtrackerMI::ErrZ2(AliTPCseed* seed, AliTPCclusterMI * cl){
688 //seed->SetErrorY2(0.1);
690 //calculate look-up table at the beginning
691 static Bool_t ginit = kFALSE;
692 static Float_t gnoise1,gnoise2,gnoise3;
693 static Float_t ggg1[10000];
694 static Float_t ggg2[10000];
695 static Float_t ggg3[10000];
696 static Float_t glandau1[10000];
697 static Float_t glandau2[10000];
698 static Float_t glandau3[10000];
700 static Float_t gcor01[1000];
701 static Float_t gcor02[1000];
702 static Float_t gcorp[1000];
707 for (Int_t i=1;i<1000;i++){
708 Float_t rsigma = float(i)/100.;
709 gcor02[i] = TMath::Max(0.81 +TMath::Exp(6.8*(rsigma-1.2)),0.6);
710 gcor01[i] = TMath::Max(0.72 +TMath::Exp(2.04*(rsigma-1.2)),0.6);
711 gcorp[i] = TMath::Max(TMath::Power((rsigma+0.5),1.5),1.2);
715 for (Int_t i=3;i<10000;i++){
719 Float_t amp = float(i);
720 Float_t padlength =0.75;
721 gnoise1 = 0.0004/padlength;
722 Float_t nel = 0.268*amp;
723 Float_t nprim = 0.155*amp;
724 ggg1[i] = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.001*nel/(padlength*padlength))/nel;
725 glandau1[i] = (2.+0.12*nprim)*0.5* (2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
726 if (glandau1[i]>1) glandau1[i]=1;
727 glandau1[i]*=padlength*padlength/12.;
731 gnoise2 = 0.0004/padlength;
734 ggg2[i] = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
735 glandau2[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
736 if (glandau2[i]>1) glandau2[i]=1;
737 glandau2[i]*=padlength*padlength/12.;
742 gnoise3 = 0.0004/padlength;
745 ggg3[i] = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
746 glandau3[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
747 if (glandau3[i]>1) glandau3[i]=1;
748 glandau3[i]*=padlength*padlength/12.;
756 Int_t amp = int(TMath::Abs(cl->GetQ()));
758 seed->SetErrorY2(1.);
762 Float_t z = TMath::Abs(fParam->GetZLength()-TMath::Abs(seed->GetZ()));
763 Int_t ctype = cl->GetType();
764 Float_t padlength= GetPadPitchLength(seed->fRow);
766 Float_t angle2 = seed->GetSnp()*seed->GetSnp();
767 // if (angle2<0.6) angle2 = 0.6;
768 angle2 = seed->GetTgl()*seed->GetTgl()*(1+angle2/(1-angle2));
771 Int_t rsigmaz = int(100.*cl->GetSigmaZ2()/(seed->fCurrentSigmaZ2));
774 if (fSectors==fInnerSec){
776 res = ggg1[amp]*z+glandau1[amp]*angle2;
777 if (ctype==0) res *= gcor01[rsigmaz];
780 res*= gcorp[rsigmaz];
786 res = ggg2[amp]*z+glandau2[amp]*angle2;
787 if (ctype==0) res *= gcor02[rsigmaz];
790 res*= gcorp[rsigmaz];
795 res = ggg3[amp]*z+glandau3[amp]*angle2;
796 if (ctype==0) res *= gcor02[rsigmaz];
799 res*= gcorp[rsigmaz];
808 if ((ctype<0) &&<70){
816 seed->SetErrorZ2(res);
823 Double_t AliTPCtrackerMI::ErrZ2(AliTPCseed* seed, AliTPCclusterMI * cl){
826 //seed->SetErrorZ2(0.1);
830 Float_t z = TMath::Abs(fParam->GetZLength()-TMath::Abs(seed->GetZ()));
832 Float_t rsigmaz = cl->GetSigmaZ2()/(seed->fCurrentSigmaZ2);
833 Int_t ctype = cl->GetType();
834 Float_t amp = TMath::Abs(cl->GetQ());
839 Float_t landau=2 ; //landau fluctuation part
840 Float_t gg=2; // gg fluctuation part
841 Float_t padlength= GetPadPitchLength(seed->GetX());
843 if (fSectors==fInnerSec){
844 snoise2 = 0.0004/padlength;
847 gg = (2+0.001*nel/(padlength*padlength))/nel;
848 landau = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
849 if (landau>1) landau=1;
852 snoise2 = 0.0004/padlength;
855 gg = (2+0.0008*nel/(padlength*padlength))/nel;
856 landau = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
857 if (landau>1) landau=1;
859 Float_t sdiff = gg*fParam->GetDiffT()*fParam->GetDiffT()*z;
862 Float_t angle2 = seed->GetSnp()*seed->GetSnp();
863 angle2 = TMath::Sqrt((1-angle2));
864 if (angle2<0.6) angle2 = 0.6;
867 Float_t angle = seed->GetTgl()/angle2;
868 Float_t angular = landau*angle*angle*padlength*padlength/12.;
869 Float_t res = sdiff + angular;
872 if ((ctype==0) && (fSectors ==fOuterSec))
873 res *= 0.81 +TMath::Exp(6.8*(rsigmaz-1.2));
875 if ((ctype==0) && (fSectors ==fInnerSec))
876 res *= 0.72 +TMath::Exp(2.04*(rsigmaz-1.2));
880 res*= TMath::Power(rsigmaz+0.5,1.5); //0.31+0.147*ctype;
886 if ((ctype<0) &&<70){
894 seed->SetErrorZ2(res);
901 void AliTPCseed::Reset(Bool_t all)
905 SetNumberOfClusters(0);
911 for (Int_t i=0;i<8;i++){
912 delete [] fTrackPoints[i];
920 for (Int_t i=0;i<200;i++) SetClusterIndex2(i,-3);
921 for (Int_t i=0;i<160;i++) fClusterPointer[i]=0;
927 void AliTPCseed::Modify(Double_t factor)
930 //------------------------------------------------------------------
931 //This function makes a track forget its history :)
932 //------------------------------------------------------------------
938 fC10*=0; fC11*=factor;
939 fC20*=0; fC21*=0; fC22*=factor;
940 fC30*=0; fC31*=0; fC32*=0; fC33*=factor;
941 fC40*=0; fC41*=0; fC42*=0; fC43*=0; fC44*=factor;
942 SetNumberOfClusters(0);
946 fCurrentSigmaY2 = 0.000005;
947 fCurrentSigmaZ2 = 0.000005;
956 Int_t AliTPCseed::GetProlongation(Double_t xk, Double_t &y, Double_t & z) const
958 //-----------------------------------------------------------------
959 // This function find proloncation of a track to a reference plane x=xk.
960 // doesn't change internal state of the track
961 //-----------------------------------------------------------------
963 Double_t x1=fX, x2=x1+(xk-x1), dx=x2-x1;
965 if (TMath::Abs(fP4*xk - fP2) >= 0.999) {
969 // Double_t y1=fP0, z1=fP1;
970 Double_t c1=fP4*x1 - fP2, r1=sqrt(1.- c1*c1);
971 Double_t c2=fP4*x2 - fP2, r2=sqrt(1.- c2*c2);
975 //y += dx*(c1+c2)/(r1+r2);
976 //z += dx*(c1+c2)/(c1*r2 + c2*r1)*fP3;
978 Double_t dy = dx*(c1+c2)/(r1+r2);
981 Double_t delta = fP4*dx*(c1+c2)/(c1*r2 + c2*r1);
983 if (TMath::Abs(delta)>0.0001){
984 dz = fP3*TMath::ASin(delta)/fP4;
986 dz = dx*fP3*(c1+c2)/(c1*r2 + c2*r1);
989 dz = fP3*AliTPCFastMath::FastAsin(delta)/fP4;
999 //_____________________________________________________________________________
1000 Double_t AliTPCseed::GetPredictedChi2(const AliTPCclusterMI *c) const
1002 //-----------------------------------------------------------------
1003 // This function calculates a predicted chi2 increment.
1004 //-----------------------------------------------------------------
1005 //Double_t r00=c->GetSigmaY2(), r01=0., r11=c->GetSigmaZ2();
1006 Double_t r00=fErrorY2, r01=0., r11=fErrorZ2;
1007 r00+=fC00; r01+=fC10; r11+=fC11;
1009 Double_t det=r00*r11 - r01*r01;
1010 if (TMath::Abs(det) < 1.e-10) {
1011 Int_t n=GetNumberOfClusters();
1012 if (n>4) cerr<<n<<" AliKalmanTrack warning: Singular matrix !\n";
1015 Double_t tmp=r00; r00=r11; r11=tmp; r01=-r01;
1017 Double_t dy=c->GetY() - fP0, dz=c->GetZ() - fP1;
1019 return (dy*r00*dy + 2*r01*dy*dz + dz*r11*dz)/det;
1023 //_________________________________________________________________________________________
1026 Int_t AliTPCseed::Compare(const TObject *o) const {
1027 //-----------------------------------------------------------------
1028 // This function compares tracks according to the sector - for given sector according z
1029 //-----------------------------------------------------------------
1030 AliTPCseed *t=(AliTPCseed*)o;
1033 if (t->fRelativeSector>fRelativeSector) return -1;
1034 if (t->fRelativeSector<fRelativeSector) return 1;
1035 Double_t z2 = t->GetZ();
1036 Double_t z1 = GetZ();
1037 if (z2>z1) return 1;
1038 if (z2<z1) return -1;
1043 f2 = 1-20*TMath::Sqrt(t->fC44)/(TMath::Abs(t->GetC())+0.0066);
1044 if (t->fBConstrain) f2=1.2;
1047 f1 = 1-20*TMath::Sqrt(fC44)/(TMath::Abs(GetC())+0.0066);
1049 if (fBConstrain) f1=1.2;
1051 if (t->GetNumberOfClusters()*f2 <GetNumberOfClusters()*f1) return -1;
1056 void AliTPCtrackerMI::RotateToLocal(AliTPCseed *seed)
1058 //rotate to track "local coordinata
1059 Float_t x = seed->GetX();
1060 Float_t y = seed->GetY();
1061 Float_t ymax = x*TMath::Tan(0.5*fSectors->GetAlpha());
1064 seed->fRelativeSector= (seed->fRelativeSector+1) % fN;
1065 if (!seed->Rotate(fSectors->GetAlpha()))
1067 } else if (y <-ymax) {
1068 seed->fRelativeSector= (seed->fRelativeSector-1+fN) % fN;
1069 if (!seed->Rotate(-fSectors->GetAlpha()))
1078 //_____________________________________________________________________________
1079 Int_t AliTPCseed::Update(const AliTPCclusterMI *c, Double_t chisq, UInt_t /*index*/) {
1080 //-----------------------------------------------------------------
1081 // This function associates a cluster with this track.
1082 //-----------------------------------------------------------------
1083 Double_t r00=fErrorY2, r01=0., r11=fErrorZ2;
1085 r00+=fC00; r01+=fC10; r11+=fC11;
1086 Double_t det=r00*r11 - r01*r01;
1087 Double_t tmp=r00; r00=r11/det; r11=tmp/det; r01=-r01/det;
1089 Double_t k00=fC00*r00+fC10*r01, k01=fC00*r01+fC10*r11;
1090 Double_t k10=fC10*r00+fC11*r01, k11=fC10*r01+fC11*r11;
1091 Double_t k20=fC20*r00+fC21*r01, k21=fC20*r01+fC21*r11;
1092 Double_t k30=fC30*r00+fC31*r01, k31=fC30*r01+fC31*r11;
1093 Double_t k40=fC40*r00+fC41*r01, k41=fC40*r01+fC41*r11;
1095 Double_t dy=c->GetY() - fP0, dz=c->GetZ() - fP1;
1096 Double_t cur=fP4 + k40*dy + k41*dz, eta=fP2 + k20*dy + k21*dz;
1097 if (TMath::Abs(cur*fX-eta) >= 0.9) {
1101 fP0 += k00*dy + k01*dz;
1102 fP1 += k10*dy + k11*dz;
1104 fP3 += k30*dy + k31*dz;
1107 Double_t c01=fC10, c02=fC20, c03=fC30, c04=fC40;
1108 Double_t c12=fC21, c13=fC31, c14=fC41;
1110 fC00-=k00*fC00+k01*fC10; fC10-=k00*c01+k01*fC11;
1111 fC20-=k00*c02+k01*c12; fC30-=k00*c03+k01*c13;
1112 fC40-=k00*c04+k01*c14;
1114 fC11-=k10*c01+k11*fC11;
1115 fC21-=k10*c02+k11*c12; fC31-=k10*c03+k11*c13;
1116 fC41-=k10*c04+k11*c14;
1118 fC22-=k20*c02+k21*c12; fC32-=k20*c03+k21*c13;
1119 fC42-=k20*c04+k21*c14;
1121 fC33-=k30*c03+k31*c13;
1122 fC43-=k40*c03+k41*c13;
1124 fC44-=k40*c04+k41*c14;
1126 Int_t n=GetNumberOfClusters();
1128 SetNumberOfClusters(n+1);
1129 SetChi2(GetChi2()+chisq);
1136 //_____________________________________________________________________________
1137 Double_t AliTPCtrackerMI::F1old(Double_t x1,Double_t y1,
1138 Double_t x2,Double_t y2,
1139 Double_t x3,Double_t y3)
1141 //-----------------------------------------------------------------
1142 // Initial approximation of the track curvature
1143 //-----------------------------------------------------------------
1144 Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
1145 Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
1146 (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
1147 Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
1148 (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
1150 Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
1151 if ( xr*xr+yr*yr<=0.00000000000001) return 100;
1152 return -xr*yr/sqrt(xr*xr+yr*yr);
1157 //_____________________________________________________________________________
1158 Double_t AliTPCtrackerMI::F1(Double_t x1,Double_t y1,
1159 Double_t x2,Double_t y2,
1160 Double_t x3,Double_t y3)
1162 //-----------------------------------------------------------------
1163 // Initial approximation of the track curvature
1164 //-----------------------------------------------------------------
1170 Double_t det = x3*y2-x2*y3;
1175 Double_t u = 0.5* (x2*(x2-x3)+y2*(y2-y3))/det;
1176 Double_t x0 = x3*0.5-y3*u;
1177 Double_t y0 = y3*0.5+x3*u;
1178 Double_t c2 = 1/TMath::Sqrt(x0*x0+y0*y0);
1184 Double_t AliTPCtrackerMI::F2(Double_t x1,Double_t y1,
1185 Double_t x2,Double_t y2,
1186 Double_t x3,Double_t y3)
1188 //-----------------------------------------------------------------
1189 // Initial approximation of the track curvature
1190 //-----------------------------------------------------------------
1196 Double_t det = x3*y2-x2*y3;
1201 Double_t u = 0.5* (x2*(x2-x3)+y2*(y2-y3))/det;
1202 Double_t x0 = x3*0.5-y3*u;
1203 Double_t y0 = y3*0.5+x3*u;
1204 Double_t c2 = 1/TMath::Sqrt(x0*x0+y0*y0);
1213 //_____________________________________________________________________________
1214 Double_t AliTPCtrackerMI::F2old(Double_t x1,Double_t y1,
1215 Double_t x2,Double_t y2,
1216 Double_t x3,Double_t y3)
1218 //-----------------------------------------------------------------
1219 // Initial approximation of the track curvature times center of curvature
1220 //-----------------------------------------------------------------
1221 Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
1222 Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
1223 (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
1224 Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
1225 (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
1227 Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
1229 return -a/(d*y1-b)*xr/sqrt(xr*xr+yr*yr);
1232 //_____________________________________________________________________________
1233 Double_t AliTPCtrackerMI::F3(Double_t x1,Double_t y1,
1234 Double_t x2,Double_t y2,
1235 Double_t z1,Double_t z2)
1237 //-----------------------------------------------------------------
1238 // Initial approximation of the tangent of the track dip angle
1239 //-----------------------------------------------------------------
1240 return (z1 - z2)/sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
1244 Double_t AliTPCtrackerMI::F3n(Double_t x1,Double_t y1,
1245 Double_t x2,Double_t y2,
1246 Double_t z1,Double_t z2, Double_t c)
1248 //-----------------------------------------------------------------
1249 // Initial approximation of the tangent of the track dip angle
1250 //-----------------------------------------------------------------
1254 //angle1 = (z1-z2)*c/(TMath::ASin(c*x1-ni)-TMath::ASin(c*x2-ni));
1256 Double_t d = TMath::Sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
1257 if (TMath::Abs(d*c*0.5)>1) return 0;
1258 // Double_t angle2 = TMath::ASin(d*c*0.5);
1259 // Double_t angle2 = AliTPCFastMath::FastAsin(d*c*0.5);
1260 Double_t angle2 = (d*c*0.5>0.1)? TMath::ASin(d*c*0.5): AliTPCFastMath::FastAsin(d*c*0.5);
1262 angle2 = (z1-z2)*c/(angle2*2.);
1266 Bool_t AliTPCtrackerMI::GetProlongation(Double_t x1, Double_t x2, Double_t x[5], Double_t &y, Double_t &z)
1267 {//-----------------------------------------------------------------
1268 // This function find proloncation of a track to a reference plane x=x2.
1269 //-----------------------------------------------------------------
1273 if (TMath::Abs(x[4]*x1 - x[2]) >= 0.999) {
1277 Double_t c1=x[4]*x1 - x[2], r1=sqrt(1.- c1*c1);
1278 Double_t c2=x[4]*x2 - x[2], r2=sqrt(1.- c2*c2);
1282 Double_t dy = dx*(c1+c2)/(r1+r2);
1285 Double_t delta = x[4]*dx*(c1+c2)/(c1*r2 + c2*r1);
1287 if (TMath::Abs(delta)>0.01){
1288 dz = x[3]*TMath::ASin(delta)/x[4];
1290 dz = x[3]*AliTPCFastMath::FastAsin(delta)/x[4];
1293 //dz = x[3]*AliTPCFastMath::FastAsin(delta)/x[4];
1301 Int_t AliTPCtrackerMI::LoadClusters (TTree *tree)
1306 return LoadClusters();
1309 Int_t AliTPCtrackerMI::LoadClusters()
1312 // load clusters to the memory
1313 AliTPCClustersRow *clrow= new AliTPCClustersRow;
1314 clrow->SetClass("AliTPCclusterMI");
1316 clrow->GetArray()->ExpandCreateFast(10000);
1318 // TTree * tree = fClustersArray.GetTree();
1320 TTree * tree = fInput;
1321 TBranch * br = tree->GetBranch("Segment");
1322 br->SetAddress(&clrow);
1324 Int_t j=Int_t(tree->GetEntries());
1325 for (Int_t i=0; i<j; i++) {
1329 fParam->AdjustSectorRow(clrow->GetID(),sec,row);
1331 AliTPCRow * tpcrow=0;
1334 tpcrow = &(fInnerSec[sec%fkNIS][row]);
1338 tpcrow = &(fOuterSec[(sec-fkNIS*2)%fkNOS][row]);
1339 left = (sec-fkNIS*2)/fkNOS;
1342 tpcrow->fN1 = clrow->GetArray()->GetEntriesFast();
1343 tpcrow->fClusters1 = new AliTPCclusterMI[tpcrow->fN1];
1344 for (Int_t i=0;i<tpcrow->fN1;i++)
1345 tpcrow->fClusters1[i] = *(AliTPCclusterMI*)(clrow->GetArray()->At(i));
1348 tpcrow->fN2 = clrow->GetArray()->GetEntriesFast();
1349 tpcrow->fClusters2 = new AliTPCclusterMI[tpcrow->fN2];
1350 for (Int_t i=0;i<tpcrow->fN2;i++)
1351 tpcrow->fClusters2[i] = *(AliTPCclusterMI*)(clrow->GetArray()->At(i));
1362 void AliTPCtrackerMI::UnloadClusters()
1365 // unload clusters from the memory
1367 Int_t nrows = fOuterSec->GetNRows();
1368 for (Int_t sec = 0;sec<fkNOS;sec++)
1369 for (Int_t row = 0;row<nrows;row++){
1370 AliTPCRow* tpcrow = &(fOuterSec[sec%fkNOS][row]);
1372 // if (tpcrow->fClusters1) delete []tpcrow->fClusters1;
1373 // if (tpcrow->fClusters2) delete []tpcrow->fClusters2;
1375 tpcrow->ResetClusters();
1378 nrows = fInnerSec->GetNRows();
1379 for (Int_t sec = 0;sec<fkNIS;sec++)
1380 for (Int_t row = 0;row<nrows;row++){
1381 AliTPCRow* tpcrow = &(fInnerSec[sec%fkNIS][row]);
1383 // if (tpcrow->fClusters1) delete []tpcrow->fClusters1;
1384 //if (tpcrow->fClusters2) delete []tpcrow->fClusters2;
1386 tpcrow->ResetClusters();
1393 //_____________________________________________________________________________
1394 Int_t AliTPCtrackerMI::LoadOuterSectors() {
1395 //-----------------------------------------------------------------
1396 // This function fills outer TPC sectors with clusters.
1397 //-----------------------------------------------------------------
1398 Int_t nrows = fOuterSec->GetNRows();
1400 for (Int_t sec = 0;sec<fkNOS;sec++)
1401 for (Int_t row = 0;row<nrows;row++){
1402 AliTPCRow* tpcrow = &(fOuterSec[sec%fkNOS][row]);
1403 Int_t sec2 = sec+2*fkNIS;
1405 Int_t ncl = tpcrow->fN1;
1407 AliTPCclusterMI *c= &(tpcrow->fClusters1[ncl]);
1408 index=(((sec2<<8)+row)<<16)+ncl;
1409 tpcrow->InsertCluster(c,index);
1414 AliTPCclusterMI *c= &(tpcrow->fClusters2[ncl]);
1415 index=((((sec2+fkNOS)<<8)+row)<<16)+ncl;
1416 tpcrow->InsertCluster(c,index);
1419 // write indexes for fast acces
1421 for (Int_t i=0;i<510;i++)
1422 tpcrow->fFastCluster[i]=-1;
1423 for (Int_t i=0;i<tpcrow->GetN();i++){
1424 Int_t zi = Int_t((*tpcrow)[i]->GetZ()+255.);
1425 tpcrow->fFastCluster[zi]=i; // write index
1428 for (Int_t i=0;i<510;i++){
1429 if (tpcrow->fFastCluster[i]<0)
1430 tpcrow->fFastCluster[i] = last;
1432 last = tpcrow->fFastCluster[i];
1441 //_____________________________________________________________________________
1442 Int_t AliTPCtrackerMI::LoadInnerSectors() {
1443 //-----------------------------------------------------------------
1444 // This function fills inner TPC sectors with clusters.
1445 //-----------------------------------------------------------------
1446 Int_t nrows = fInnerSec->GetNRows();
1448 for (Int_t sec = 0;sec<fkNIS;sec++)
1449 for (Int_t row = 0;row<nrows;row++){
1450 AliTPCRow* tpcrow = &(fInnerSec[sec%fkNIS][row]);
1453 Int_t ncl = tpcrow->fN1;
1455 AliTPCclusterMI *c= &(tpcrow->fClusters1[ncl]);
1456 index=(((sec<<8)+row)<<16)+ncl;
1457 tpcrow->InsertCluster(c,index);
1462 AliTPCclusterMI *c= &(tpcrow->fClusters2[ncl]);
1463 index=((((sec+fkNIS)<<8)+row)<<16)+ncl;
1464 tpcrow->InsertCluster(c,index);
1467 // write indexes for fast acces
1469 for (Int_t i=0;i<510;i++)
1470 tpcrow->fFastCluster[i]=-1;
1471 for (Int_t i=0;i<tpcrow->GetN();i++){
1472 Int_t zi = Int_t((*tpcrow)[i]->GetZ()+255.);
1473 tpcrow->fFastCluster[zi]=i; // write index
1476 for (Int_t i=0;i<510;i++){
1477 if (tpcrow->fFastCluster[i]<0)
1478 tpcrow->fFastCluster[i] = last;
1480 last = tpcrow->fFastCluster[i];
1492 //_________________________________________________________________________
1493 AliTPCclusterMI *AliTPCtrackerMI::GetClusterMI(Int_t index) const {
1494 //--------------------------------------------------------------------
1495 // Return pointer to a given cluster
1496 //--------------------------------------------------------------------
1497 Int_t sec=(index&0xff000000)>>24;
1498 Int_t row=(index&0x00ff0000)>>16;
1499 Int_t ncl=(index&0x00007fff)>>00;
1501 const AliTPCRow * tpcrow=0;
1502 AliTPCclusterMI * clrow =0;
1505 tpcrow = &(fInnerSec[sec%fkNIS][row]);
1506 if (tpcrow==0) return 0;
1509 if (tpcrow->fN1<=ncl) return 0;
1510 clrow = tpcrow->fClusters1;
1513 if (tpcrow->fN2<=ncl) return 0;
1514 clrow = tpcrow->fClusters2;
1518 tpcrow = &(fOuterSec[(sec-fkNIS*2)%fkNOS][row]);
1519 if (tpcrow==0) return 0;
1521 if (sec-2*fkNIS<fkNOS) {
1522 if (tpcrow->fN1<=ncl) return 0;
1523 clrow = tpcrow->fClusters1;
1526 if (tpcrow->fN2<=ncl) return 0;
1527 clrow = tpcrow->fClusters2;
1531 return &(clrow[ncl]);
1537 Int_t AliTPCtrackerMI::FollowToNext(AliTPCseed& t, Int_t nr) {
1538 //-----------------------------------------------------------------
1539 // This function tries to find a track prolongation to next pad row
1540 //-----------------------------------------------------------------
1542 Double_t x= GetXrow(nr), ymax=GetMaxY(nr);
1543 AliTPCclusterMI *cl=0;
1544 Int_t tpcindex= t.GetClusterIndex2(nr);
1546 // update current shape info every 5 pad-row
1547 // if ( (nr%5==0) || t.GetNumberOfClusters()<2 || (t.fCurrentSigmaY2<0.0001) ){
1551 if (fIteration>0 && tpcindex>=-1){ //if we have already clusters
1553 if (tpcindex==-1) return 0; //track in dead zone
1555 cl = t.fClusterPointer[nr];
1556 if ( (cl==0) ) cl = GetClusterMI(tpcindex);
1557 t.fCurrentClusterIndex1 = tpcindex;
1560 Int_t relativesector = ((tpcindex&0xff000000)>>24)%18; // if previously accepted cluster in different sector
1561 Float_t angle = relativesector*fSectors->GetAlpha()+fSectors->GetAlphaShift();
1563 if (angle<-TMath::Pi()) angle += 2*TMath::Pi();
1564 if (angle>=TMath::Pi()) angle -= 2*TMath::Pi();
1566 if (TMath::Abs(angle-t.GetAlpha())>0.001){
1567 Double_t rotation = angle-t.GetAlpha();
1568 t.fRelativeSector= relativesector;
1573 t.fCurrentCluster = cl;
1575 Int_t accept = AcceptCluster(&t,t.fCurrentCluster,1.);
1576 if ((tpcindex&0x8000)==0) accept =0;
1578 //if founded cluster is acceptible
1579 if (cl->IsUsed(11)) { // id cluster is shared inrease uncertainty
1586 UpdateTrack(&t,accept);
1591 if (fIteration>1) return 0; // not look for new cluster during refitting
1594 if (TMath::Abs(t.GetSnp())>0.95 || TMath::Abs(x*t.GetC()-t.GetEta())>0.95) return 0;
1595 Double_t y=t.GetYat(x);
1596 if (TMath::Abs(y)>ymax){
1598 t.fRelativeSector= (t.fRelativeSector+1) % fN;
1599 if (!t.Rotate(fSectors->GetAlpha()))
1601 } else if (y <-ymax) {
1602 t.fRelativeSector= (t.fRelativeSector-1+fN) % fN;
1603 if (!t.Rotate(-fSectors->GetAlpha()))
1609 if (!t.PropagateTo(x)) {
1610 if (fIteration==0) t.fRemoval = 10;
1614 Double_t z=t.GetZ();
1616 const AliTPCRow &krow=GetRow(t.fRelativeSector,nr);
1617 if ( (t.GetSigmaY2()<0) || t.GetSigmaZ2()<0) return 0;
1619 Double_t roadz = 1.;
1621 if (TMath::Abs(TMath::Abs(y)-ymax)<krow.fDeadZone){
1623 t.SetClusterIndex2(nr,-1);
1628 if (TMath::Abs(z)<(AliTPCReconstructor::GetCtgRange()*x+10)) t.fNFoundable++;
1634 // cl = krow.FindNearest2(y+10.,z,roady,roadz,index);
1635 cl = krow.FindNearest2(y,z,roady,roadz,index);
1636 if (cl) t.fCurrentClusterIndex1 = krow.GetIndex(index);
1639 t.fCurrentCluster = cl;
1641 if (fIteration==2&&cl->IsUsed(10)) return 0;
1642 Int_t accept = AcceptCluster(&t,t.fCurrentCluster,1.);
1643 if (fIteration==2&&cl->IsUsed(11)) {
1650 if (t.fCurrentCluster->IsUsed(10)){
1655 if (t.fNShared>0.7*t.GetNumberOfClusters()) {
1661 if (accept<3) UpdateTrack(&t,accept);
1664 if ( fIteration==0 && t.fNFoundable*0.5 > t.GetNumberOfClusters()) t.fRemoval=10;
1670 Int_t AliTPCtrackerMI::FollowToNextFast(AliTPCseed& t, Int_t nr) {
1671 //-----------------------------------------------------------------
1672 // This function tries to find a track prolongation to next pad row
1673 //-----------------------------------------------------------------
1675 Double_t x= GetXrow(nr), ymax=GetMaxY(nr);
1677 if (!t.GetProlongation(x,y,z)) {
1683 if (TMath::Abs(y)>ymax){
1686 t.fRelativeSector= (t.fRelativeSector+1) % fN;
1687 if (!t.Rotate(fSectors->GetAlpha()))
1689 } else if (y <-ymax) {
1690 t.fRelativeSector= (t.fRelativeSector-1+fN) % fN;
1691 if (!t.Rotate(-fSectors->GetAlpha()))
1694 if (!t.PropagateTo(x)) {
1697 t.GetProlongation(x,y,z);
1700 // update current shape info every 3 pad-row
1701 if ( (nr%6==0) || t.GetNumberOfClusters()<2 || (t.fCurrentSigmaY2<0.0001) ){
1702 // t.fCurrentSigmaY = GetSigmaY(&t);
1703 //t.fCurrentSigmaZ = GetSigmaZ(&t);
1707 AliTPCclusterMI *cl=0;
1712 const AliTPCRow &krow=GetRow(t.fRelativeSector,nr);
1713 if ( (t.GetSigmaY2()<0) || t.GetSigmaZ2()<0) return 0;
1715 Double_t roadz = 1.;
1718 if (TMath::Abs(TMath::Abs(y)-ymax)<krow.fDeadZone){
1720 t.SetClusterIndex2(row,-1);
1725 if (TMath::Abs(z)>(AliTPCReconstructor::GetCtgRange()*x+10)) t.SetClusterIndex2(row,-1);
1729 if ((cl==0)&&(krow)) {
1730 // cl = krow.FindNearest2(y+10,z,roady,roadz,index);
1731 cl = krow.FindNearest2(y,z,roady,roadz,index);
1733 if (cl) t.fCurrentClusterIndex1 = krow.GetIndex(index);
1737 t.fCurrentCluster = cl;
1738 // Int_t accept = AcceptCluster(&t,t.fCurrentCluster,1.);
1740 t.SetClusterIndex2(row,index);
1741 t.fClusterPointer[row] = cl;
1749 Int_t AliTPCtrackerMI::UpdateClusters(AliTPCseed& t, Int_t nr) {
1750 //-----------------------------------------------------------------
1751 // This function tries to find a track prolongation to next pad row
1752 //-----------------------------------------------------------------
1753 t.fCurrentCluster = 0;
1754 t.fCurrentClusterIndex1 = 0;
1756 Double_t xt=t.GetX();
1757 Int_t row = GetRowNumber(xt)-1;
1758 Double_t ymax= GetMaxY(nr);
1760 if (row < nr) return 1; // don't prolongate if not information until now -
1761 if (TMath::Abs(t.GetSnp())>0.9 && t.GetNumberOfClusters()>40. && fIteration!=2) {
1763 return 0; // not prolongate strongly inclined tracks
1765 if (TMath::Abs(t.GetSnp())>0.95) {
1767 return 0; // not prolongate strongly inclined tracks
1770 Double_t x= GetXrow(nr);
1772 //t.PropagateTo(x+0.02);
1773 //t.PropagateTo(x+0.01);
1774 if (!t.PropagateTo(x)){
1781 if (TMath::Abs(y)>ymax){
1783 t.fRelativeSector= (t.fRelativeSector+1) % fN;
1784 if (!t.Rotate(fSectors->GetAlpha()))
1786 } else if (y <-ymax) {
1787 t.fRelativeSector= (t.fRelativeSector-1+fN) % fN;
1788 if (!t.Rotate(-fSectors->GetAlpha()))
1791 // if (!t.PropagateTo(x)){
1799 AliTPCRow &krow=GetRow(t.fRelativeSector,nr);
1801 if (TMath::Abs(TMath::Abs(y)-ymax)<krow.fDeadZone){
1803 t.SetClusterIndex2(nr,-1);
1808 if (TMath::Abs(t.GetZ())<(AliTPCReconstructor::GetCtgRange()*t.GetX()+10)) t.fNFoundable++;
1814 if ( (nr%6==0) || t.GetNumberOfClusters()<2){
1815 // t.fCurrentSigmaY = GetSigmaY(&t);
1816 //t.fCurrentSigmaZ = GetSigmaZ(&t);
1820 AliTPCclusterMI *cl=0;
1823 Double_t roady = 1.;
1824 Double_t roadz = 1.;
1828 index = t.GetClusterIndex2(nr);
1829 if ( (index>0) && (index&0x8000)==0){
1830 cl = t.fClusterPointer[nr];
1831 if ( (cl==0) && (index>0)) cl = GetClusterMI(index);
1832 t.fCurrentClusterIndex1 = index;
1834 t.fCurrentCluster = cl;
1841 //cl = krow.FindNearest2(y+10,z,roady,roadz,index);
1842 cl = krow.FindNearest2(y,z,roady,roadz,index);
1845 if (cl) t.fCurrentClusterIndex1 = krow.GetIndex(index);
1846 t.fCurrentCluster = cl;
1852 Int_t AliTPCtrackerMI::FollowToNextCluster(AliTPCseed & t, Int_t nr) {
1853 //-----------------------------------------------------------------
1854 // This function tries to find a track prolongation to next pad row
1855 //-----------------------------------------------------------------
1857 //update error according neighborhoud
1859 if (t.fCurrentCluster) {
1861 Int_t accept = AcceptCluster(&t,t.fCurrentCluster,1.);
1863 if (t.fCurrentCluster->IsUsed(10)){
1869 if (t.fNShared>0.7*t.GetNumberOfClusters()) {
1874 if (fIteration>0) accept = 0;
1875 if (accept<3) UpdateTrack(&t,accept);
1879 if ( ( (t.GetSigmaY2()+t.GetSigmaZ2())>0.16)&& t.GetNumberOfClusters()>18) t.fRemoval=10;
1880 if ( t.GetChi2()/t.GetNumberOfClusters()>6 &&t.GetNumberOfClusters()>18) t.fRemoval=10;
1882 if (( (t.fNFoundable*0.5 > t.GetNumberOfClusters()) || t.fNoCluster>15)) t.fRemoval=10;
1890 //_____________________________________________________________________________
1891 Int_t AliTPCtrackerMI::FollowProlongation(AliTPCseed& t, Int_t rf, Int_t step) {
1892 //-----------------------------------------------------------------
1893 // This function tries to find a track prolongation.
1894 //-----------------------------------------------------------------
1895 Double_t xt=t.GetX();
1897 Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
1898 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
1899 if (alpha < 0. ) alpha += 2.*TMath::Pi();
1901 t.fRelativeSector = Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN;
1903 Int_t first = GetRowNumber(xt)-1;
1904 for (Int_t nr= first; nr>=rf; nr-=step) {
1906 if (t.GetKinkIndexes()[0]>0){
1907 for (Int_t i=0;i<3;i++){
1908 Int_t index = t.GetKinkIndexes()[i];
1909 if (index==0) break;
1910 if (index<0) continue;
1912 AliESDkink * kink = fEvent->GetKink(index-1);
1914 printf("PROBLEM\n");
1917 Int_t kinkrow = kink->GetTPCRow0()+2+Int_t(0.5/(0.05+kink->GetAngle(2)));
1919 AliExternalTrackParam paramd(t);
1920 kink->SetDaughter(paramd);
1921 kink->SetStatus(2,5);
1928 if (nr==80) t.UpdateReference();
1929 if (nr<fInnerSec->GetNRows())
1930 fSectors = fInnerSec;
1932 fSectors = fOuterSec;
1933 if (FollowToNext(t,nr)==0)
1942 //_____________________________________________________________________________
1943 Int_t AliTPCtrackerMI::FollowProlongationFast(AliTPCseed& t, Int_t rf, Int_t step) {
1944 //-----------------------------------------------------------------
1945 // This function tries to find a track prolongation.
1946 //-----------------------------------------------------------------
1947 Double_t xt=t.GetX();
1949 Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
1950 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
1951 if (alpha < 0. ) alpha += 2.*TMath::Pi();
1952 t.fRelativeSector = Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN;
1954 for (Int_t nr=GetRowNumber(xt)-1; nr>=rf; nr-=step) {
1956 if (FollowToNextFast(t,nr)==0)
1957 if (!t.IsActive()) return 0;
1967 Int_t AliTPCtrackerMI::FollowBackProlongation(AliTPCseed& t, Int_t rf) {
1968 //-----------------------------------------------------------------
1969 // This function tries to find a track prolongation.
1970 //-----------------------------------------------------------------
1972 Double_t xt=t.GetX();
1973 Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
1974 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
1975 if (alpha < 0. ) alpha += 2.*TMath::Pi();
1976 t.fRelativeSector = Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN;
1978 Int_t first = t.fFirstPoint;
1979 if (first<GetRowNumber(xt)+1) first = GetRowNumber(xt)+1;
1981 if (first<0) first=0;
1982 for (Int_t nr=first; nr<=rf; nr++) {
1983 //if ( (t.GetSnp()<0.9))
1984 if (t.GetKinkIndexes()[0]<0){
1985 for (Int_t i=0;i<3;i++){
1986 Int_t index = t.GetKinkIndexes()[i];
1987 if (index==0) break;
1988 if (index>0) continue;
1989 index = TMath::Abs(index);
1990 AliESDkink * kink = fEvent->GetKink(index-1);
1992 printf("PROBLEM\n");
1995 Int_t kinkrow = kink->GetTPCRow0()-2-Int_t(0.5/(0.05+kink->GetAngle(2)));
1997 AliExternalTrackParam paramm(t);
1998 kink->SetMother(paramm);
1999 kink->SetStatus(2,1);
2006 if (nr<fInnerSec->GetNRows())
2007 fSectors = fInnerSec;
2009 fSectors = fOuterSec;
2019 Float_t AliTPCtrackerMI::OverlapFactor(AliTPCseed * s1, AliTPCseed * s2, Int_t &sum1, Int_t & sum2)
2027 Float_t dz2 =(s1->GetZ() - s2->GetZ());
2030 Float_t dy2 =TMath::Abs((s1->GetY() - s2->GetY()));
2032 Float_t distance = TMath::Sqrt(dz2+dy2);
2033 if (distance>4.) return 0; // if there are far away - not overlap - to reduce combinatorics
2036 Int_t firstpoint = TMath::Min(s1->fFirstPoint,s2->fFirstPoint);
2037 Int_t lastpoint = TMath::Max(s1->fLastPoint,s2->fLastPoint);
2042 if (firstpoint>lastpoint) {
2043 firstpoint =lastpoint;
2048 for (Int_t i=firstpoint-1;i<lastpoint+1;i++){
2049 if (s1->GetClusterIndex2(i)>0) sum1++;
2050 if (s2->GetClusterIndex2(i)>0) sum2++;
2051 if (s1->GetClusterIndex2(i)==s2->GetClusterIndex2(i) && s1->GetClusterIndex2(i)>0) {
2055 if (sum<5) return 0;
2057 Float_t summin = TMath::Min(sum1+1,sum2+1);
2058 Float_t ratio = (sum+1)/Float_t(summin);
2062 void AliTPCtrackerMI::SignShared(AliTPCseed * s1, AliTPCseed * s2)
2066 if (TMath::Abs(s1->GetC()-s2->GetC())>0.004) return;
2067 if (TMath::Abs(s1->GetTgl()-s2->GetTgl())>0.6) return;
2069 Float_t dz2 =(s1->GetZ() - s2->GetZ());
2071 Float_t dy2 =(s1->GetY() - s2->GetY());
2073 Float_t distance = dz2+dy2;
2074 if (distance>325.) return ; // if there are far away - not overlap - to reduce combinatorics
2079 Int_t firstpoint = TMath::Max(s1->fFirstPoint,s2->fFirstPoint);
2080 Int_t lastpoint = TMath::Min(s1->fLastPoint,s2->fLastPoint);
2082 if (firstpoint>=lastpoint-5) return;;
2084 for (Int_t i=firstpoint;i<lastpoint;i++){
2085 // if ( (s1->GetClusterIndex2(i)&0xFFFF8FFF)==(s2->GetClusterIndex2(i)&0xFFFF8FFF) && s1->GetClusterIndex2(i)>0) {
2086 if ( (s1->GetClusterIndex2(i))==(s2->GetClusterIndex2(i)) && s1->GetClusterIndex2(i)>0) {
2093 for (Int_t i=firstpoint;i<lastpoint;i++){
2094 // if ( (s1->GetClusterIndex2(i)&0xFFFF8FFF)==(s2->GetClusterIndex2(i)&0xFFFF8FFF) && s1->GetClusterIndex2(i)>0) {
2095 if ( (s1->GetClusterIndex2(i))==(s2->GetClusterIndex2(i)) && s1->GetClusterIndex2(i)>0) {
2096 AliTPCTrackerPoint *p1 = s1->GetTrackPoint(i);
2097 AliTPCTrackerPoint *p2 = s2->GetTrackPoint(i);;
2098 if (s1->IsActive()&&s2->IsActive()){
2099 p1->fIsShared = kTRUE;
2100 p2->fIsShared = kTRUE;
2107 for (Int_t i=0;i<4;i++){
2108 if (s1->fOverlapLabels[3*i]==0){
2109 s1->fOverlapLabels[3*i] = s2->GetLabel();
2110 s1->fOverlapLabels[3*i+1] = sumshared;
2111 s1->fOverlapLabels[3*i+2] = s2->GetUniqueID();
2115 for (Int_t i=0;i<4;i++){
2116 if (s2->fOverlapLabels[3*i]==0){
2117 s2->fOverlapLabels[3*i] = s1->GetLabel();
2118 s2->fOverlapLabels[3*i+1] = sumshared;
2119 s2->fOverlapLabels[3*i+2] = s1->GetUniqueID();
2127 void AliTPCtrackerMI::SignShared(TObjArray * arr)
2130 //sort trackss according sectors
2132 for (Int_t i=0; i<arr->GetEntriesFast(); i++) {
2133 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2135 //if (pt) RotateToLocal(pt);
2139 arr->Sort(); // sorting according z
2140 arr->Expand(arr->GetEntries());
2143 Int_t nseed=arr->GetEntriesFast();
2144 for (Int_t i=0; i<nseed; i++) {
2145 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2147 for (Int_t j=0;j<=12;j++){
2148 pt->fOverlapLabels[j] =0;
2151 for (Int_t i=0; i<nseed; i++) {
2152 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2154 if (pt->fRemoval>10) continue;
2155 for (Int_t j=i+1; j<nseed; j++){
2156 AliTPCseed *pt2=(AliTPCseed*)arr->UncheckedAt(j);
2158 if (pt2->fRemoval<=10) {
2159 if ( TMath::Abs(pt->fRelativeSector-pt2->fRelativeSector)>0) break;
2166 void AliTPCtrackerMI::RemoveDouble(TObjArray * arr, Float_t factor1, Float_t factor2, Int_t removalindex)
2169 //sort trackss according sectors
2172 Info("RemoveDouble","Number of tracks before double removal- %d\n",arr->GetEntries());
2175 for (Int_t i=0; i<arr->GetEntriesFast(); i++) {
2176 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2181 arr->Sort(); // sorting according z
2182 arr->Expand(arr->GetEntries());
2184 //reset overlap labels
2186 Int_t nseed=arr->GetEntriesFast();
2187 for (Int_t i=0; i<nseed; i++) {
2188 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2191 for (Int_t j=0;j<=12;j++){
2192 pt->fOverlapLabels[j] =0;
2196 //sign shared tracks
2197 for (Int_t i=0; i<nseed; i++) {
2198 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2200 if (pt->fRemoval>10) continue;
2201 Float_t deltac = pt->GetC()*0.1;
2202 for (Int_t j=i+1; j<nseed; j++){
2203 AliTPCseed *pt2=(AliTPCseed*)arr->UncheckedAt(j);
2205 if (pt2->fRemoval<=10) {
2206 if ( TMath::Abs(pt->fRelativeSector-pt2->fRelativeSector)>0) break;
2207 if (TMath::Abs(pt->GetC() -pt2->GetC())>deltac) continue;
2208 if (TMath::Abs(pt->GetTgl()-pt2->GetTgl())>0.05) continue;
2215 // remove highly shared tracks
2216 for (Int_t i=0; i<nseed; i++) {
2217 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2219 if (pt->fRemoval>10) continue;
2222 for (Int_t j=0;j<4;j++){
2223 sumshared = pt->fOverlapLabels[j*3+1];
2225 Float_t factor = factor1;
2226 if (pt->fRemoval>0) factor = factor2;
2227 if (sumshared/pt->GetNumberOfClusters()>factor){
2228 for (Int_t j=0;j<4;j++){
2229 if (pt->fOverlapLabels[3*j]==0) continue;
2230 if (pt->fOverlapLabels[3*j+1]<5) continue;
2231 if (pt->fRemoval==removalindex) continue;
2232 AliTPCseed * pt2 = (AliTPCseed*)arr->UncheckedAt(pt->fOverlapLabels[3*j+2]);
2234 if (pt2->GetSigma2C()<pt->GetSigma2C()){
2235 // pt->fRemoval = removalindex;
2236 delete arr->RemoveAt(i);
2244 Info("RemoveDouble","Number of tracks after double removal- %d\n",arr->GetEntries());
2253 void AliTPCtrackerMI::SortTracks(TObjArray * arr, Int_t mode) const
2256 //sort tracks in array according mode criteria
2257 Int_t nseed = arr->GetEntriesFast();
2258 for (Int_t i=0; i<nseed; i++) {
2259 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2269 void AliTPCtrackerMI::RemoveUsed(TObjArray * arr, Float_t factor1, Float_t factor2, Int_t removalindex)
2272 //Loop over all tracks and remove "overlaps"
2275 Int_t nseed = arr->GetEntriesFast();
2278 for (Int_t i=0; i<nseed; i++) {
2279 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2281 delete arr->RemoveAt(i);
2285 pt->fBSigned = kFALSE;
2289 nseed = arr->GetEntriesFast();
2296 for (Int_t i=0; i<nseed; i++) {
2297 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2301 Int_t found,foundable,shared;
2303 pt->GetClusterStatistic(0,160,found, foundable,shared,kFALSE);
2305 pt->GetClusterStatistic(0,160,found, foundable,shared,kTRUE);
2307 Double_t factor = factor2;
2308 if (pt->fBConstrain) factor = factor1;
2310 if ((Float_t(shared)/Float_t(found))>factor){
2311 pt->Desactivate(removalindex);
2316 for (Int_t i=0; i<160; i++) {
2317 Int_t index=pt->GetClusterIndex2(i);
2318 if (index<0 || index&0x8000 ) continue;
2319 AliTPCclusterMI *c= pt->fClusterPointer[i];
2321 // if (!c->IsUsed(10)) c->Use(10);
2322 //if (pt->IsActive())
2331 Info("RemoveUsed","\n*****\nNumber of good tracks after shared removal\t%d\n",fNtracks);
2336 void AliTPCtrackerMI::RemoveUsed2(TObjArray * arr, Float_t factor1, Float_t factor2, Int_t minimal)
2339 //Loop over all tracks and remove "overlaps"
2344 Int_t nseed = arr->GetEntriesFast();
2345 Float_t * quality = new Float_t[nseed];
2346 Int_t * indexes = new Int_t[nseed];
2350 for (Int_t i=0; i<nseed; i++) {
2351 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2356 pt->UpdatePoints(); //select first last max dens points
2357 Float_t * points = pt->GetPoints();
2358 if (points[3]<0.8) quality[i] =-1;
2360 quality[i] = (points[2]-points[0])+pt->GetNumberOfClusters();
2362 TMath::Sort(nseed,quality,indexes);
2365 for (Int_t itrack=0; itrack<nseed; itrack++) {
2366 Int_t trackindex = indexes[itrack];
2367 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(trackindex);
2368 if (quality[trackindex]<0){
2370 delete arr->RemoveAt(trackindex);
2373 arr->RemoveAt(trackindex);
2378 Int_t first = Int_t(pt->GetPoints()[0]);
2379 Int_t last = Int_t(pt->GetPoints()[2]);
2380 Double_t factor = (pt->fBConstrain) ? factor1: factor2;
2382 Int_t found,foundable,shared;
2383 pt->GetClusterStatistic(first,last, found, foundable,shared,kFALSE);
2384 Float_t sharedfactor = Float_t(shared+1)/Float_t(found+1);
2386 if (Float_t(shared+1)/Float_t(found+1)>factor){
2387 if (pt->GetKinkIndexes()[0]!=0) continue; //don't remove tracks - part of the kinks
2388 delete arr->RemoveAt(trackindex);
2392 if (pt->GetNumberOfClusters()<50&&(found-0.5*shared)<minimal){ //remove short tracks
2393 if (pt->GetKinkIndexes()[0]!=0) continue; //don't remove tracks - part of the kinks
2394 delete arr->RemoveAt(trackindex);
2399 if (sharedfactor>0.4) continue;
2400 for (Int_t i=first; i<last; i++) {
2401 Int_t index=pt->GetClusterIndex2(i);
2402 // if (index<0 || index&0x8000 ) continue;
2403 if (index<0 || index&0x8000 ) continue;
2404 AliTPCclusterMI *c= pt->fClusterPointer[i];
2411 Info("RemoveUsed","\n*****\nNumber of good tracks after shared removal\t%d\n",fNtracks);
2417 void AliTPCtrackerMI::UnsignClusters()
2420 // loop over all clusters and unsign them
2423 for (Int_t sec=0;sec<fkNIS;sec++){
2424 for (Int_t row=0;row<fInnerSec->GetNRows();row++){
2425 AliTPCclusterMI *cl = fInnerSec[sec][row].fClusters1;
2426 for (Int_t icl =0;icl< fInnerSec[sec][row].fN1;icl++)
2427 // if (cl[icl].IsUsed(10))
2429 cl = fInnerSec[sec][row].fClusters2;
2430 for (Int_t icl =0;icl< fInnerSec[sec][row].fN2;icl++)
2431 //if (cl[icl].IsUsed(10))
2436 for (Int_t sec=0;sec<fkNOS;sec++){
2437 for (Int_t row=0;row<fOuterSec->GetNRows();row++){
2438 AliTPCclusterMI *cl = fOuterSec[sec][row].fClusters1;
2439 for (Int_t icl =0;icl< fOuterSec[sec][row].fN1;icl++)
2440 //if (cl[icl].IsUsed(10))
2442 cl = fOuterSec[sec][row].fClusters2;
2443 for (Int_t icl =0;icl< fOuterSec[sec][row].fN2;icl++)
2444 //if (cl[icl].IsUsed(10))
2453 void AliTPCtrackerMI::SignClusters(TObjArray * arr, Float_t fnumber, Float_t fdensity)
2456 //sign clusters to be "used"
2458 // snumber and sdensity sign number of sigmas - bellow mean value to be accepted
2459 // loop over "primaries"
2473 Int_t nseed = arr->GetEntriesFast();
2474 for (Int_t i=0; i<nseed; i++) {
2475 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2479 if (!(pt->IsActive())) continue;
2480 Float_t dens = pt->GetNumberOfClusters()/Float_t(pt->fNFoundable);
2481 if ( (dens>0.7) && (pt->GetNumberOfClusters()>70)){
2483 sumdens2+= dens*dens;
2484 sumn += pt->GetNumberOfClusters();
2485 sumn2 += pt->GetNumberOfClusters()*pt->GetNumberOfClusters();
2486 Float_t chi2 = pt->GetChi2()/pt->GetNumberOfClusters();
2489 sumchi2 +=chi2*chi2;
2494 Float_t mdensity = 0.9;
2495 Float_t meann = 130;
2496 Float_t meanchi = 1;
2497 Float_t sdensity = 0.1;
2498 Float_t smeann = 10;
2499 Float_t smeanchi =0.4;
2503 mdensity = sumdens/sum;
2505 meanchi = sumchi/sum;
2507 sdensity = sumdens2/sum-mdensity*mdensity;
2508 sdensity = TMath::Sqrt(sdensity);
2510 smeann = sumn2/sum-meann*meann;
2511 smeann = TMath::Sqrt(smeann);
2513 smeanchi = sumchi2/sum - meanchi*meanchi;
2514 smeanchi = TMath::Sqrt(smeanchi);
2518 //REMOVE SHORT DELTAS or tracks going out of sensitive volume of TPC
2520 for (Int_t i=0; i<nseed; i++) {
2521 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2525 if (pt->fBSigned) continue;
2526 if (pt->fBConstrain) continue;
2527 //if (!(pt->IsActive())) continue;
2529 Int_t found,foundable,shared;
2530 pt->GetClusterStatistic(0,160,found, foundable,shared);
2531 if (shared/float(found)>0.3) {
2532 if (shared/float(found)>0.9 ){
2533 //delete arr->RemoveAt(i);
2538 Bool_t isok =kFALSE;
2539 if ( (pt->fNShared/pt->GetNumberOfClusters()<0.5) &&pt->GetNumberOfClusters()>60)
2541 if ((TMath::Abs(1/pt->GetC())<100.) && (pt->fNShared/pt->GetNumberOfClusters()<0.7))
2543 if (TMath::Abs(pt->GetZ()/pt->GetX())>1.1)
2545 if ( (TMath::Abs(pt->GetSnp()>0.7) && pt->GetD(0,0)>60.))
2549 for (Int_t i=0; i<160; i++) {
2550 Int_t index=pt->GetClusterIndex2(i);
2551 if (index<0) continue;
2552 AliTPCclusterMI *c= pt->fClusterPointer[i];
2554 //if (!(c->IsUsed(10))) c->Use();
2561 Double_t maxchi = meanchi+2.*smeanchi;
2563 for (Int_t i=0; i<nseed; i++) {
2564 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2568 //if (!(pt->IsActive())) continue;
2569 if (pt->fBSigned) continue;
2570 Double_t chi = pt->GetChi2()/pt->GetNumberOfClusters();
2571 if (chi>maxchi) continue;
2574 Float_t dens = pt->GetNumberOfClusters()/Float_t(pt->fNFoundable);
2576 //sign only tracks with enoug big density at the beginning
2578 if ((pt->GetDensityFirst(40)<0.75) && pt->GetNumberOfClusters()<meann) continue;
2581 Double_t mindens = TMath::Max(double(mdensity-sdensity*fdensity*bfactor),0.65);
2582 Double_t minn = TMath::Max(Int_t(meann-fnumber*smeann*bfactor),50);
2584 // if (pt->fBConstrain) mindens = TMath::Max(mdensity-sdensity*fdensity*bfactor,0.65);
2585 if ( (pt->fRemoval==10) && (pt->GetSnp()>0.8)&&(dens>mindens))
2588 if ((dens>mindens && pt->GetNumberOfClusters()>minn) && chi<maxchi ){
2589 //Int_t noc=pt->GetNumberOfClusters();
2590 pt->fBSigned = kTRUE;
2591 for (Int_t i=0; i<160; i++) {
2593 Int_t index=pt->GetClusterIndex2(i);
2594 if (index<0) continue;
2595 AliTPCclusterMI *c= pt->fClusterPointer[i];
2597 // if (!(c->IsUsed(10))) c->Use();
2602 // gLastCheck = nseed;
2610 void AliTPCtrackerMI::StopNotActive(TObjArray * arr, Int_t row0, Float_t th0, Float_t th1, Float_t th2) const
2612 // stop not active tracks
2613 // take th1 as threshold for number of founded to number of foundable on last 10 active rows
2614 // take th2 as threshold for number of founded to number of foundable on last 20 active rows
2615 Int_t nseed = arr->GetEntriesFast();
2617 for (Int_t i=0; i<nseed; i++) {
2618 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2622 if (!(pt->IsActive())) continue;
2623 StopNotActive(pt,row0,th0, th1,th2);
2629 void AliTPCtrackerMI::StopNotActive(AliTPCseed * seed, Int_t row0, Float_t th0, Float_t th1,
2632 // stop not active tracks
2633 // take th1 as threshold for number of founded to number of foundable on last 10 active rows
2634 // take th2 as threshold for number of founded to number of foundable on last 20 active rows
2637 Int_t foundable = 0;
2638 Int_t maxindex = seed->fLastPoint; //last foundable row
2639 if (seed->fNFoundable*th0 > seed->GetNumberOfClusters()) {
2640 seed->Desactivate(10) ;
2644 for (Int_t i=row0; i<maxindex; i++){
2645 Int_t index = seed->GetClusterIndex2(i);
2646 if (index!=-1) foundable++;
2648 if (foundable<=30) sumgood1++;
2649 if (foundable<=50) {
2656 if (foundable>=30.){
2657 if (sumgood1<(th1*30.)) seed->Desactivate(10);
2660 if (sumgood2<(th2*50.)) seed->Desactivate(10);
2664 Int_t AliTPCtrackerMI::RefitInward(AliESD *event)
2667 // back propagation of ESD tracks
2673 //PrepareForProlongation(fSeeds,1);
2674 PropagateForward2(fSeeds);
2676 Int_t nseed = fSeeds->GetEntriesFast();
2677 for (Int_t i=0;i<nseed;i++){
2678 AliTPCseed * seed = (AliTPCseed*) fSeeds->UncheckedAt(i);
2679 if (!seed) continue;
2680 if (seed->GetKinkIndex(0)>0) UpdateKinkQualityD(seed); // update quality informations for kinks
2682 seed->PropagateTo(fParam->GetInnerRadiusLow());
2683 seed->UpdatePoints();
2684 AliESDtrack *esd=event->GetTrack(i);
2685 seed->CookdEdx(0.02,0.6);
2686 CookLabel(seed,0.1); //For comparison only
2687 if (seed->GetNumberOfClusters()>15){
2688 esd->UpdateTrackParams(seed,AliESDtrack::kTPCrefit);
2689 esd->SetTPCPoints(seed->GetPoints());
2693 //printf("problem\n");
2696 //FindKinks(fSeeds,event);
2697 Info("RefitInward","Number of refitted tracks %d",ntracks);
2704 Int_t AliTPCtrackerMI::PropagateBack(AliESD *event)
2707 // back propagation of ESD tracks
2713 PropagateBack(fSeeds);
2714 Int_t nseed = fSeeds->GetEntriesFast();
2716 for (Int_t i=0;i<nseed;i++){
2717 AliTPCseed * seed = (AliTPCseed*) fSeeds->UncheckedAt(i);
2718 if (!seed) continue;
2719 if (seed->GetKinkIndex(0)<0) UpdateKinkQualityM(seed); // update quality informations for kinks
2720 seed->UpdatePoints();
2721 AliESDtrack *esd=event->GetTrack(i);
2722 seed->CookdEdx(0.02,0.6);
2723 CookLabel(seed,0.1); //For comparison only
2724 if (seed->GetNumberOfClusters()>15){
2725 esd->UpdateTrackParams(seed,AliESDtrack::kTPCout);
2726 esd->SetTPCPoints(seed->GetPoints());
2730 //FindKinks(fSeeds,event);
2731 Info("PropagateBack","Number of back propagated tracks %d",ntracks);
2738 void AliTPCtrackerMI::DeleteSeeds()
2742 Int_t nseed = fSeeds->GetEntriesFast();
2743 for (Int_t i=0;i<nseed;i++){
2744 AliTPCseed * seed = (AliTPCseed*)fSeeds->At(i);
2745 if (seed) delete fSeeds->RemoveAt(i);
2751 void AliTPCtrackerMI::ReadSeeds(AliESD *event, Int_t direction)
2754 //read seeds from the event
2756 Int_t nentr=event->GetNumberOfTracks();
2758 Info("ReadSeeds", "Number of ESD tracks: %d\n", nentr);
2763 fSeeds = new TObjArray(nentr);
2767 for (Int_t i=0; i<nentr; i++) {
2768 AliESDtrack *esd=event->GetTrack(i);
2769 ULong_t status=esd->GetStatus();
2770 if (!(status&AliESDtrack::kTPCin)) continue;
2771 AliTPCtrack t(*esd);
2772 // AliTPCseed *seed = new AliTPCseed(t,t.GetAlpha());
2773 AliTPCseed *seed = new AliTPCseed(t/*,t.GetAlpha()*/);
2774 for (Int_t ikink=0;ikink<3;ikink++) {
2775 Int_t index = esd->GetKinkIndex(ikink);
2776 seed->GetKinkIndexes()[ikink] = index;
2777 if (index==0) continue;
2778 index = TMath::Abs(index);
2779 AliESDkink * kink = fEvent->GetKink(index-1);
2780 if (kink&&esd->GetKinkIndex(ikink)<0){
2781 if ((status & AliESDtrack::kTRDrefit) != 0) kink->SetStatus(1,2);
2782 if ((status & AliESDtrack::kITSout) != 0) kink->SetStatus(1,0);
2784 if (kink&&esd->GetKinkIndex(ikink)>0){
2785 if ((status & AliESDtrack::kTRDrefit) != 0) kink->SetStatus(1,6);
2786 if ((status & AliESDtrack::kITSout) != 0) kink->SetStatus(1,4);
2790 if ((status==AliESDtrack::kTPCin)&&(direction==1)) seed->ResetCovariance();
2791 if ( direction ==2 &&(status & AliESDtrack::kTRDrefit) == 0 ) seed->ResetCovariance();
2792 if ( direction ==2 && ((status & AliESDtrack::kTPCout) == 0) ) {
2797 if ( direction ==2 &&(status & AliESDtrack::kTRDrefit) > 0 ) {
2798 Double_t par0[5],par1[5],x;
2799 esd->GetInnerExternalParameters(x,par0);
2800 esd->GetExternalParameters(x,par1);
2801 Double_t delta1 = TMath::Abs(par0[4]-par1[4])/(0.000000001+TMath::Abs(par0[4]+par1[4]));
2802 Double_t delta2 = TMath::Abs(par0[3]-par1[3]);
2804 if (esd->GetTRDncls()>0) trdchi2 = esd->GetTRDchi2()/esd->GetTRDncls();
2805 //reset covariance if suspicious
2806 if ( (delta1>0.1) || (delta2>0.006) ||trdchi2>7.)
2807 seed->ResetCovariance();
2812 // rotate to the local coordinate system
2814 fSectors=fInnerSec; fN=fkNIS;
2815 Double_t alpha=seed->GetAlpha() - fSectors->GetAlphaShift();
2816 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
2817 if (alpha < 0. ) alpha += 2.*TMath::Pi();
2818 Int_t ns=Int_t(alpha/fSectors->GetAlpha())%fN;
2819 alpha =ns*fSectors->GetAlpha() + fSectors->GetAlphaShift();
2820 if (alpha<-TMath::Pi()) alpha += 2*TMath::Pi();
2821 if (alpha>=TMath::Pi()) alpha -= 2*TMath::Pi();
2822 alpha-=seed->GetAlpha();
2823 if (!seed->Rotate(alpha)) {
2829 for (Int_t irow=0;irow<160;irow++){
2830 Int_t index = seed->GetClusterIndex2(irow);
2833 AliTPCclusterMI * cl = GetClusterMI(index);
2834 seed->fClusterPointer[irow] = cl;
2836 if ((index & 0x8000)==0){
2837 cl->Use(10); // accepted cluster
2839 cl->Use(6); // close cluster not accepted
2842 Info("ReadSeeds","Not found cluster");
2846 fSeeds->AddAt(seed,i);
2852 //_____________________________________________________________________________
2853 void AliTPCtrackerMI::MakeSeeds3(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, Float_t cuts[4],
2854 Float_t deltay, Int_t ddsec) {
2855 //-----------------------------------------------------------------
2856 // This function creates track seeds.
2857 // SEEDING WITH VERTEX CONSTRAIN
2858 //-----------------------------------------------------------------
2859 // cuts[0] - fP4 cut
2860 // cuts[1] - tan(phi) cut
2861 // cuts[2] - zvertex cut
2862 // cuts[3] - fP3 cut
2870 Double_t x[5], c[15];
2871 // Int_t di = i1-i2;
2873 AliTPCseed * seed = new AliTPCseed;
2874 Double_t alpha=fSectors->GetAlpha(), shift=fSectors->GetAlphaShift();
2875 Double_t cs=cos(alpha), sn=sin(alpha);
2877 // Double_t x1 =fOuterSec->GetX(i1);
2878 //Double_t xx2=fOuterSec->GetX(i2);
2880 Double_t x1 =GetXrow(i1);
2881 Double_t xx2=GetXrow(i2);
2883 Double_t x3=GetX(), y3=GetY(), z3=GetZ();
2885 Int_t imiddle = (i2+i1)/2; //middle pad row index
2886 Double_t xm = GetXrow(imiddle); // radius of middle pad-row
2887 const AliTPCRow& krm=GetRow(sec,imiddle); //middle pad -row
2891 const AliTPCRow& kr1=GetRow(ns,i1);
2892 Double_t ymax = GetMaxY(i1)-kr1.fDeadZone-1.5;
2893 Double_t ymaxm = GetMaxY(imiddle)-kr1.fDeadZone-1.5;
2896 // change cut on curvature if it can't reach this layer
2897 // maximal curvature set to reach it
2898 Double_t dvertexmax = TMath::Sqrt((x1-x3)*(x1-x3)+(ymax+5-y3)*(ymax+5-y3));
2899 if (dvertexmax*0.5*cuts[0]>0.85){
2900 cuts[0] = 0.85/(dvertexmax*0.5+1.);
2902 Double_t r2min = 1/(cuts[0]*cuts[0]); //minimal square of radius given by cut
2905 if (deltay>0) ddsec = 0;
2906 // loop over clusters
2907 for (Int_t is=0; is < kr1; is++) {
2909 if (kr1[is]->IsUsed(10)) continue;
2910 Double_t y1=kr1[is]->GetY(), z1=kr1[is]->GetZ();
2911 //if (TMath::Abs(y1)>ymax) continue;
2913 if (deltay>0 && TMath::Abs(ymax-TMath::Abs(y1))> deltay ) continue; // seed only at the edge
2915 // find possible directions
2916 Float_t anglez = (z1-z3)/(x1-x3);
2917 Float_t extraz = z1 - anglez*(x1-xx2); // extrapolated z
2920 //find rotation angles relative to line given by vertex and point 1
2921 Double_t dvertex2 = (x1-x3)*(x1-x3)+(y1-y3)*(y1-y3);
2922 Double_t dvertex = TMath::Sqrt(dvertex2);
2923 Double_t angle13 = TMath::ATan((y1-y3)/(x1-x3));
2924 Double_t cs13 = cos(-angle13), sn13 = sin(-angle13);
2927 // loop over 2 sectors
2933 Double_t dddz1=0; // direction of delta inclination in z axis
2940 for (Int_t dsec = dsec1; dsec<=dsec2;dsec++){
2941 Int_t sec2 = sec + dsec;
2943 // AliTPCRow& kr2 = fOuterSec[(sec2+fkNOS)%fkNOS][i2];
2944 //AliTPCRow& kr2m = fOuterSec[(sec2+fkNOS)%fkNOS][imiddle];
2945 AliTPCRow& kr2 = GetRow((sec2+fkNOS)%fkNOS,i2);
2946 AliTPCRow& kr2m = GetRow((sec2+fkNOS)%fkNOS,imiddle);
2947 Int_t index1 = TMath::Max(kr2.Find(extraz-0.6-dddz1*TMath::Abs(z1)*0.05)-1,0);
2948 Int_t index2 = TMath::Min(kr2.Find(extraz+0.6+dddz2*TMath::Abs(z1)*0.05)+1,kr2);
2950 // rotation angles to p1-p3
2951 Double_t cs13r = cos(-angle13+dsec*alpha)/dvertex, sn13r = sin(-angle13+dsec*alpha)/dvertex;
2952 Double_t x2, y2, z2;
2954 // Double_t dymax = maxangle*TMath::Abs(x1-xx2);
2957 Double_t dxx0 = (xx2-x3)*cs13r;
2958 Double_t dyy0 = (xx2-x3)*sn13r;
2959 for (Int_t js=index1; js < index2; js++) {
2960 const AliTPCclusterMI *kcl = kr2[js];
2961 if (kcl->IsUsed(10)) continue;
2963 //calcutate parameters
2965 Double_t yy0 = dyy0 +(kcl->GetY()-y3)*cs13r;
2967 if (TMath::Abs(yy0)<0.000001) continue;
2968 Double_t xx0 = dxx0 -(kcl->GetY()-y3)*sn13r;
2969 Double_t y0 = 0.5*(xx0*xx0+yy0*yy0-xx0)/yy0;
2970 Double_t r02 = (0.25+y0*y0)*dvertex2;
2971 //curvature (radius) cut
2972 if (r02<r2min) continue;
2976 Double_t c0 = 1/TMath::Sqrt(r02);
2980 //Double_t dfi0 = 2.*TMath::ASin(dvertex*c0*0.5);
2981 //Double_t dfi1 = 2.*TMath::ASin(TMath::Sqrt(yy0*yy0+(1-xx0)*(1-xx0))*dvertex*c0*0.5);
2982 Double_t dfi0 = 2.*AliTPCFastMath::FastAsin(dvertex*c0*0.5);
2983 Double_t dfi1 = 2.*AliTPCFastMath::FastAsin(TMath::Sqrt(yy0*yy0+(1-xx0)*(1-xx0))*dvertex*c0*0.5);
2986 Double_t z0 = kcl->GetZ();
2987 Double_t zzzz2 = z1-(z1-z3)*dfi1/dfi0;
2988 if (TMath::Abs(zzzz2-z0)>0.5) continue;
2991 Double_t dip = (z1-z0)*c0/dfi1;
2992 Double_t x0 = (0.5*cs13+y0*sn13)*dvertex*c0;
3003 x2= xx2*cs-y2*sn*dsec;
3004 y2=+xx2*sn*dsec+y2*cs;
3014 // do we have cluster at the middle ?
3016 GetProlongation(x1,xm,x,ym,zm);
3018 AliTPCclusterMI * cm=0;
3019 if (TMath::Abs(ym)-ymaxm<0){
3020 cm = krm.FindNearest2(ym,zm,1.0,0.6,dummy);
3021 if ((!cm) || (cm->IsUsed(10))) {
3026 // rotate y1 to system 0
3027 // get state vector in rotated system
3028 Double_t yr1 = (-0.5*sn13+y0*cs13)*dvertex*c0;
3029 Double_t xr2 = x0*cs+yr1*sn*dsec;
3030 Double_t xr[5]={kcl->GetY(),kcl->GetZ(), xr2, dip, c0};
3032 GetProlongation(xx2,xm,xr,ym,zm);
3033 if (TMath::Abs(ym)-ymaxm<0){
3034 cm = kr2m.FindNearest2(ym,zm,1.0,0.6,dummy);
3035 if ((!cm) || (cm->IsUsed(10))) {
3045 dym = ym - cm->GetY();
3046 dzm = zm - cm->GetZ();
3053 Double_t sy1=kr1[is]->GetSigmaY2()*2., sz1=kr1[is]->GetSigmaZ2()*2.;
3054 Double_t sy2=kcl->GetSigmaY2()*2., sz2=kcl->GetSigmaZ2()*2.;
3055 //Double_t sy3=400*3./12., sy=0.1, sz=0.1;
3056 Double_t sy3=25000*x[4]*x[4]+0.1, sy=0.1, sz=0.1;
3057 //Double_t sy3=25000*x[4]*x[4]*60+0.5, sy=0.1, sz=0.1;
3059 Double_t f40=(F1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
3060 Double_t f42=(F1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
3061 Double_t f43=(F1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
3062 Double_t f20=(F2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
3063 Double_t f22=(F2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
3064 Double_t f23=(F2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
3066 Double_t f30=(F3(x1,y1+sy,x2,y2,z1,z2)-x[3])/sy;
3067 Double_t f31=(F3(x1,y1,x2,y2,z1+sz,z2)-x[3])/sz;
3068 Double_t f32=(F3(x1,y1,x2,y2+sy,z1,z2)-x[3])/sy;
3069 Double_t f34=(F3(x1,y1,x2,y2,z1,z2+sz)-x[3])/sz;
3073 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3074 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3075 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3076 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3077 c[13]=f30*sy1*f40+f32*sy2*f42;
3078 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3080 // if (!BuildSeed(kr1[is],kcl,0,x1,x2,x3,x,c)) continue;
3082 UInt_t index=kr1.GetIndex(is);
3083 AliTPCseed *track=new(seed) AliTPCseed(index, x, c, x1, ns*alpha+shift);
3085 track->fIsSeeding = kTRUE;
3092 FollowProlongation(*track, (i1+i2)/2,1);
3093 Int_t foundable,found,shared;
3094 track->GetClusterStatistic((i1+i2)/2,i1, found, foundable, shared, kTRUE);
3095 if ((found<0.55*foundable) || shared>0.5*found || (track->GetSigmaY2()+track->GetSigmaZ2())>0.5){
3097 seed->~AliTPCseed();
3103 FollowProlongation(*track, i2,1);
3107 track->fBConstrain =1;
3108 // track->fLastPoint = i1+fInnerSec->GetNRows(); // first cluster in track position
3109 track->fLastPoint = i1; // first cluster in track position
3110 track->fFirstPoint = track->fLastPoint;
3112 if (track->GetNumberOfClusters()<(i1-i2)*0.5 ||
3113 track->GetNumberOfClusters() < track->fNFoundable*0.6 ||
3114 track->fNShared>0.4*track->GetNumberOfClusters() ) {
3116 seed->~AliTPCseed();
3120 // Z VERTEX CONDITION
3122 zv = track->GetZ()+track->GetTgl()/track->GetC()*
3123 ( asin(-track->GetEta()) - asin(track->GetX()*track->GetC()-track->GetEta()));
3124 if (TMath::Abs(zv-z3)>cuts[2]) {
3125 FollowProlongation(*track, TMath::Max(i2-20,0));
3126 zv = track->GetZ()+track->GetTgl()/track->GetC()*
3127 ( asin(-track->GetEta()) - asin(track->GetX()*track->GetC()-track->GetEta()));
3128 if (TMath::Abs(zv-z3)>cuts[2]){
3129 FollowProlongation(*track, TMath::Max(i2-40,0));
3130 zv = track->GetZ()+track->GetTgl()/track->GetC()*
3131 ( asin(-track->GetEta()) - asin(track->GetX()*track->GetC()-track->GetEta()));
3132 if (TMath::Abs(zv-z3)>cuts[2] &&(track->GetNumberOfClusters() > track->fNFoundable*0.7)){
3133 // make seed without constrain
3134 AliTPCseed * track2 = MakeSeed(track,0.2,0.5,1.);
3135 FollowProlongation(*track2, i2,1);
3136 track2->fBConstrain = kFALSE;
3137 track2->fSeedType = 1;
3138 arr->AddLast(track2);
3140 seed->~AliTPCseed();
3145 seed->~AliTPCseed();
3152 track->fSeedType =0;
3153 arr->AddLast(track);
3154 seed = new AliTPCseed;
3156 // don't consider other combinations
3157 if (track->GetNumberOfClusters() > track->fNFoundable*0.8)
3163 Info("MakeSeeds3","\nSeeding statistic:\t%d\t%d\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin,nout1,nout2);
3169 void AliTPCtrackerMI::MakeSeeds5(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, Float_t cuts[4],
3174 //-----------------------------------------------------------------
3175 // This function creates track seeds.
3176 //-----------------------------------------------------------------
3177 // cuts[0] - fP4 cut
3178 // cuts[1] - tan(phi) cut
3179 // cuts[2] - zvertex cut
3180 // cuts[3] - fP3 cut
3190 Double_t x[5], c[15];
3192 // make temporary seed
3193 AliTPCseed * seed = new AliTPCseed;
3194 Double_t alpha=fOuterSec->GetAlpha(), shift=fOuterSec->GetAlphaShift();
3195 // Double_t cs=cos(alpha), sn=sin(alpha);
3200 Double_t x1 = GetXrow(i1-1);
3201 const AliTPCRow& kr1=GetRow(sec,i1-1);
3202 Double_t y1max = GetMaxY(i1-1)-kr1.fDeadZone-1.5;
3204 Double_t x1p = GetXrow(i1);
3205 const AliTPCRow& kr1p=GetRow(sec,i1);
3207 Double_t x1m = GetXrow(i1-2);
3208 const AliTPCRow& kr1m=GetRow(sec,i1-2);
3211 //last 3 padrow for seeding
3212 AliTPCRow& kr3 = GetRow((sec+fkNOS)%fkNOS,i1-7);
3213 Double_t x3 = GetXrow(i1-7);
3214 // Double_t y3max= GetMaxY(i1-7)-kr3.fDeadZone-1.5;
3216 AliTPCRow& kr3p = GetRow((sec+fkNOS)%fkNOS,i1-6);
3217 Double_t x3p = GetXrow(i1-6);
3219 AliTPCRow& kr3m = GetRow((sec+fkNOS)%fkNOS,i1-8);
3220 Double_t x3m = GetXrow(i1-8);
3225 Int_t im = i1-4; //middle pad row index
3226 Double_t xm = GetXrow(im); // radius of middle pad-row
3227 const AliTPCRow& krm=GetRow(sec,im); //middle pad -row
3228 // Double_t ymmax = GetMaxY(im)-kr1.fDeadZone-1.5;
3231 Double_t deltax = x1-x3;
3232 Double_t dymax = deltax*cuts[1];
3233 Double_t dzmax = deltax*cuts[3];
3235 // loop over clusters
3236 for (Int_t is=0; is < kr1; is++) {
3238 if (kr1[is]->IsUsed(10)) continue;
3239 Double_t y1=kr1[is]->GetY(), z1=kr1[is]->GetZ();
3241 if (deltay>0 && TMath::Abs(y1max-TMath::Abs(y1))> deltay ) continue; // seed only at the edge
3243 Int_t index1 = TMath::Max(kr3.Find(z1-dzmax)-1,0);
3244 Int_t index2 = TMath::Min(kr3.Find(z1+dzmax)+1,kr3);
3250 for (Int_t js=index1; js < index2; js++) {
3251 const AliTPCclusterMI *kcl = kr3[js];
3252 if (kcl->IsUsed(10)) continue;
3254 // apply angular cuts
3255 if (TMath::Abs(y1-y3)>dymax) continue;
3258 if (TMath::Abs(z1-z3)>dzmax) continue;
3260 Double_t angley = (y1-y3)/(x1-x3);
3261 Double_t anglez = (z1-z3)/(x1-x3);
3263 Double_t erry = TMath::Abs(angley)*(x1-x1m)*0.5+0.5;
3264 Double_t errz = TMath::Abs(anglez)*(x1-x1m)*0.5+0.5;
3266 Double_t yyym = angley*(xm-x1)+y1;
3267 Double_t zzzm = anglez*(xm-x1)+z1;
3269 const AliTPCclusterMI *kcm = krm.FindNearest2(yyym,zzzm,erry,errz,index);
3271 if (kcm->IsUsed(10)) continue;
3273 erry = TMath::Abs(angley)*(x1-x1m)*0.4+0.5;
3274 errz = TMath::Abs(anglez)*(x1-x1m)*0.4+0.5;
3281 // look around first
3282 const AliTPCclusterMI *kc1m = kr1m.FindNearest2(angley*(x1m-x1)+y1,
3288 if (kc1m->IsUsed(10)) used++;
3290 const AliTPCclusterMI *kc1p = kr1p.FindNearest2(angley*(x1p-x1)+y1,
3296 if (kc1p->IsUsed(10)) used++;
3298 if (used>1) continue;
3299 if (found<1) continue;
3303 const AliTPCclusterMI *kc3m = kr3m.FindNearest2(angley*(x3m-x3)+y3,
3309 if (kc3m->IsUsed(10)) used++;
3313 const AliTPCclusterMI *kc3p = kr3p.FindNearest2(angley*(x3p-x3)+y3,
3319 if (kc3p->IsUsed(10)) used++;
3323 if (used>1) continue;
3324 if (found<3) continue;
3334 x[4]=F1(x1,y1,x2,y2,x3,y3);
3335 //if (TMath::Abs(x[4]) >= cuts[0]) continue;
3338 x[2]=F2(x1,y1,x2,y2,x3,y3);
3341 x[3]=F3n(x1,y1,x2,y2,z1,z2,x[4]);
3342 //if (TMath::Abs(x[3]) > cuts[3]) continue;
3346 Double_t sy1=0.1, sz1=0.1;
3347 Double_t sy2=0.1, sz2=0.1;
3348 Double_t sy3=0.1, sy=0.1, sz=0.1;
3350 Double_t f40=(F1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
3351 Double_t f42=(F1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
3352 Double_t f43=(F1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
3353 Double_t f20=(F2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
3354 Double_t f22=(F2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
3355 Double_t f23=(F2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
3357 Double_t f30=(F3(x1,y1+sy,x2,y2,z1,z2)-x[3])/sy;
3358 Double_t f31=(F3(x1,y1,x2,y2,z1+sz,z2)-x[3])/sz;
3359 Double_t f32=(F3(x1,y1,x2,y2+sy,z1,z2)-x[3])/sy;
3360 Double_t f34=(F3(x1,y1,x2,y2,z1,z2+sz)-x[3])/sz;
3364 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3365 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3366 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3367 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3368 c[13]=f30*sy1*f40+f32*sy2*f42;
3369 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3371 // if (!BuildSeed(kr1[is],kcl,0,x1,x2,x3,x,c)) continue;
3373 UInt_t index=kr1.GetIndex(is);
3374 AliTPCseed *track=new(seed) AliTPCseed(index, x, c, x1, sec*alpha+shift);
3376 track->fIsSeeding = kTRUE;
3379 FollowProlongation(*track, i1-7,1);
3380 if (track->GetNumberOfClusters() < track->fNFoundable*0.75 ||
3381 track->fNShared>0.6*track->GetNumberOfClusters() || ( track->GetSigmaY2()+ track->GetSigmaZ2())>0.6){
3383 seed->~AliTPCseed();
3389 FollowProlongation(*track, i2,1);
3390 track->fBConstrain =0;
3391 track->fLastPoint = i1+fInnerSec->GetNRows(); // first cluster in track position
3392 track->fFirstPoint = track->fLastPoint;
3394 if (track->GetNumberOfClusters()<(i1-i2)*0.5 ||
3395 track->GetNumberOfClusters()<track->fNFoundable*0.7 ||
3396 track->fNShared>2. || track->GetChi2()/track->GetNumberOfClusters()>6 || ( track->GetSigmaY2()+ track->GetSigmaZ2())>0.5 ) {
3398 seed->~AliTPCseed();
3403 FollowProlongation(*track, TMath::Max(i2-10,0),1);
3404 AliTPCseed * track2 = MakeSeed(track,0.2,0.5,0.9);
3405 FollowProlongation(*track2, i2,1);
3406 track2->fBConstrain = kFALSE;
3407 track2->fSeedType = 4;
3408 arr->AddLast(track2);
3410 seed->~AliTPCseed();
3414 //arr->AddLast(track);
3415 //seed = new AliTPCseed;
3421 Info("MakeSeeds5","\nSeeding statiistic:\t%d\t%d\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin,nout1,nout2,nout3);
3427 //_____________________________________________________________________________
3428 void AliTPCtrackerMI::MakeSeeds2(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, Float_t */*cuts[4]*/,
3429 Float_t deltay, Bool_t /*bconstrain*/) {
3430 //-----------------------------------------------------------------
3431 // This function creates track seeds - without vertex constraint
3432 //-----------------------------------------------------------------
3433 // cuts[0] - fP4 cut - not applied
3434 // cuts[1] - tan(phi) cut
3435 // cuts[2] - zvertex cut - not applied
3436 // cuts[3] - fP3 cut
3446 Double_t alpha=fOuterSec->GetAlpha(), shift=fOuterSec->GetAlphaShift();
3447 // Double_t cs=cos(alpha), sn=sin(alpha);
3448 Int_t row0 = (i1+i2)/2;
3449 Int_t drow = (i1-i2)/2;
3450 const AliTPCRow& kr0=fSectors[sec][row0];
3453 AliTPCpolyTrack polytrack;
3454 Int_t nclusters=fSectors[sec][row0];
3455 AliTPCseed * seed = new AliTPCseed;
3460 for (Int_t is=0; is < nclusters; is++) { //LOOP over clusters
3462 Int_t nfoundable =0;
3463 for (Int_t iter =1; iter<2; iter++){ //iterations
3464 const AliTPCRow& krm=fSectors[sec][row0-iter];
3465 const AliTPCRow& krp=fSectors[sec][row0+iter];
3466 const AliTPCclusterMI * cl= kr0[is];
3468 if (cl->IsUsed(10)) {
3474 Double_t x = kr0.GetX();
3475 // Initialization of the polytrack
3480 Double_t y0= cl->GetY();
3481 Double_t z0= cl->GetZ();
3485 Double_t ymax = fSectors->GetMaxY(row0)-kr0.fDeadZone-1.5;
3486 if (deltay>0 && TMath::Abs(ymax-TMath::Abs(y0))> deltay ) continue; // seed only at the edge
3488 erry = (0.5)*cl->GetSigmaY2()/TMath::Sqrt(cl->GetQ())*6;
3489 errz = (0.5)*cl->GetSigmaZ2()/TMath::Sqrt(cl->GetQ())*6;
3490 polytrack.AddPoint(x,y0,z0,erry, errz);
3493 if (cl->IsUsed(10)) sumused++;
3496 Float_t roady = (5*TMath::Sqrt(cl->GetSigmaY2()+0.2)+1.)*iter;
3497 Float_t roadz = (5*TMath::Sqrt(cl->GetSigmaZ2()+0.2)+1.)*iter;
3500 AliTPCclusterMI * cl1 = krm.FindNearest(y0,z0,roady,roadz);
3501 if (cl1 && TMath::Abs(ymax-TMath::Abs(y0))) {
3502 erry = (0.5)*cl1->GetSigmaY2()/TMath::Sqrt(cl1->GetQ())*3;
3503 errz = (0.5)*cl1->GetSigmaZ2()/TMath::Sqrt(cl1->GetQ())*3;
3504 if (cl1->IsUsed(10)) sumused++;
3505 polytrack.AddPoint(x,cl1->GetY(),cl1->GetZ(),erry,errz);
3509 AliTPCclusterMI * cl2 = krp.FindNearest(y0,z0,roady,roadz);
3511 erry = (0.5)*cl2->GetSigmaY2()/TMath::Sqrt(cl2->GetQ())*3;
3512 errz = (0.5)*cl2->GetSigmaZ2()/TMath::Sqrt(cl2->GetQ())*3;
3513 if (cl2->IsUsed(10)) sumused++;
3514 polytrack.AddPoint(x,cl2->GetY(),cl2->GetZ(),erry,errz);
3517 if (sumused>0) continue;
3519 polytrack.UpdateParameters();
3525 nfoundable = polytrack.GetN();
3526 nfound = nfoundable;
3528 for (Int_t ddrow = iter+1; ddrow<drow;ddrow++){
3529 Float_t maxdist = 0.8*(1.+3./(ddrow));
3530 for (Int_t delta = -1;delta<=1;delta+=2){
3531 Int_t row = row0+ddrow*delta;
3532 kr = &(fSectors[sec][row]);
3533 Double_t xn = kr->GetX();
3534 Double_t ymax = fSectors->GetMaxY(row)-kr->fDeadZone-1.5;
3535 polytrack.GetFitPoint(xn,yn,zn);
3536 if (TMath::Abs(yn)>ymax) continue;
3538 AliTPCclusterMI * cln = kr->FindNearest(yn,zn,roady,roadz);
3540 Float_t dist = TMath::Sqrt( (yn-cln->GetY())*(yn-cln->GetY())+(zn-cln->GetZ())*(zn-cln->GetZ()));
3543 erry = (dist+0.3)*cln->GetSigmaY2()/TMath::Sqrt(cln->GetQ())*(1.+1./(ddrow));
3544 errz = (dist+0.3)*cln->GetSigmaZ2()/TMath::Sqrt(cln->GetQ())*(1.+1./(ddrow));
3545 if (cln->IsUsed(10)) {
3546 // printf("used\n");
3554 polytrack.AddPoint(xn,cln->GetY(),cln->GetZ(),erry, errz);
3559 if ( (sumused>3) || (sumused>0.5*nfound) || (nfound<0.6*nfoundable)) break;
3560 polytrack.UpdateParameters();
3563 if ( (sumused>3) || (sumused>0.5*nfound)) {
3564 //printf("sumused %d\n",sumused);
3569 polytrack.GetFitDerivation(kr0.GetX(),dy,dz);
3570 AliTPCpolyTrack track2;
3572 polytrack.Refit(track2,0.5+TMath::Abs(dy)*0.3,0.4+TMath::Abs(dz)*0.3);
3573 if (track2.GetN()<0.5*nfoundable) continue;
3576 if ((nfound>0.6*nfoundable) &&( nfoundable>0.4*(i1-i2))) {
3578 // test seed with and without constrain
3579 for (Int_t constrain=0; constrain<=0;constrain++){
3580 // add polytrack candidate
3582 Double_t x[5], c[15];
3583 Double_t x1,x2,x3,y1,y2,y3,z1,z2,z3;
3584 track2.GetBoundaries(x3,x1);
3586 track2.GetFitPoint(x1,y1,z1);
3587 track2.GetFitPoint(x2,y2,z2);
3588 track2.GetFitPoint(x3,y3,z3);
3590 //is track pointing to the vertex ?
3593 polytrack.GetFitPoint(x0,y0,z0);
3606 x[4]=F1(x1,y1,x2,y2,x3,y3);
3608 // if (TMath::Abs(x[4]) >= cuts[0]) continue; //
3609 x[2]=F2(x1,y1,x2,y2,x3,y3);
3611 //if (TMath::Abs(x[4]*x1-x[2]) >= cuts[1]) continue;
3612 //x[3]=F3(x1,y1,x2,y2,z1,z2);
3613 x[3]=F3n(x1,y1,x3,y3,z1,z3,x[4]);
3614 //if (TMath::Abs(x[3]) > cuts[3]) continue;
3617 Double_t sy =0.1, sz =0.1;
3618 Double_t sy1=0.02, sz1=0.02;
3619 Double_t sy2=0.02, sz2=0.02;
3623 sy3=25000*x[4]*x[4]+0.1, sy=0.1, sz=0.1;
3626 Double_t f40=(F1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
3627 Double_t f42=(F1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
3628 Double_t f43=(F1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
3629 Double_t f20=(F2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
3630 Double_t f22=(F2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
3631 Double_t f23=(F2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
3633 Double_t f30=(F3(x1,y1+sy,x3,y3,z1,z3)-x[3])/sy;
3634 Double_t f31=(F3(x1,y1,x3,y3,z1+sz,z3)-x[3])/sz;
3635 Double_t f32=(F3(x1,y1,x3,y3+sy,z1,z3)-x[3])/sy;
3636 Double_t f34=(F3(x1,y1,x3,y3,z1,z3+sz)-x[3])/sz;
3641 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3642 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3643 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3644 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3645 c[13]=f30*sy1*f40+f32*sy2*f42;
3646 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3648 //Int_t row1 = fSectors->GetRowNumber(x1);
3649 Int_t row1 = GetRowNumber(x1);
3653 AliTPCseed *track=new (seed) AliTPCseed(index, x, c, x1, sec*alpha+shift);
3654 track->fIsSeeding = kTRUE;
3655 Int_t rc=FollowProlongation(*track, i2);
3656 if (constrain) track->fBConstrain =1;
3658 track->fBConstrain =0;
3659 track->fLastPoint = row1+fInnerSec->GetNRows(); // first cluster in track position
3660 track->fFirstPoint = track->fLastPoint;
3662 if (rc==0 || track->GetNumberOfClusters()<(i1-i2)*0.5 ||
3663 track->GetNumberOfClusters() < track->fNFoundable*0.6 ||
3664 track->fNShared>0.4*track->GetNumberOfClusters()) {
3667 seed->~AliTPCseed();
3670 arr->AddLast(track);
3671 seed = new AliTPCseed;
3675 } // if accepted seed
3678 Info("MakeSeeds2","\nSeeding statiistic:\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin3);
3684 AliTPCseed *AliTPCtrackerMI::MakeSeed(AliTPCseed *track, Float_t r0, Float_t r1, Float_t r2)
3688 //reseed using track points
3689 Int_t p0 = int(r0*track->GetNumberOfClusters()); // point 0
3690 Int_t p1 = int(r1*track->GetNumberOfClusters());
3691 Int_t p2 = int(r2*track->GetNumberOfClusters()); // last point
3693 Double_t x0[3],x1[3],x2[3];
3698 // find track position at given ratio of the length
3699 Int_t sec0, sec1, sec2;
3705 for (Int_t i=0;i<160;i++){
3706 if (track->fClusterPointer[i]){
3708 AliTPCTrackerPoint *trpoint =track->GetTrackPoint(i);
3709 if ( (index<p0) || x0[0]<0 ){
3710 if (trpoint->GetX()>1){
3711 clindex = track->GetClusterIndex2(i);
3713 x0[0] = trpoint->GetX();
3714 x0[1] = trpoint->GetY();
3715 x0[2] = trpoint->GetZ();
3716 sec0 = ((clindex&0xff000000)>>24)%18;
3721 if ( (index<p1) &&(trpoint->GetX()>1)){
3722 clindex = track->GetClusterIndex2(i);
3724 x1[0] = trpoint->GetX();
3725 x1[1] = trpoint->GetY();
3726 x1[2] = trpoint->GetZ();
3727 sec1 = ((clindex&0xff000000)>>24)%18;
3730 if ( (index<p2) &&(trpoint->GetX()>1)){
3731 clindex = track->GetClusterIndex2(i);
3733 x2[0] = trpoint->GetX();
3734 x2[1] = trpoint->GetY();
3735 x2[2] = trpoint->GetZ();
3736 sec2 = ((clindex&0xff000000)>>24)%18;
3743 Double_t alpha, cs,sn, xx2,yy2;
3745 alpha = (sec1-sec2)*fSectors->GetAlpha();
3746 cs = TMath::Cos(alpha);
3747 sn = TMath::Sin(alpha);
3748 xx2= x1[0]*cs-x1[1]*sn;
3749 yy2= x1[0]*sn+x1[1]*cs;
3753 alpha = (sec0-sec2)*fSectors->GetAlpha();
3754 cs = TMath::Cos(alpha);
3755 sn = TMath::Sin(alpha);
3756 xx2= x0[0]*cs-x0[1]*sn;
3757 yy2= x0[0]*sn+x0[1]*cs;
3763 Double_t x[5],c[15];
3767 x[4]=F1(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]);
3768 // if (x[4]>1) return 0;
3769 x[2]=F2(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]);
3770 x[3]=F3n(x2[0],x2[1],x0[0],x0[1],x2[2],x0[2],x[4]);
3771 //if (TMath::Abs(x[3]) > 2.2) return 0;
3772 //if (TMath::Abs(x[2]) > 1.99) return 0;
3774 Double_t sy =0.1, sz =0.1;
3776 Double_t sy1=0.02+track->GetSigmaY2(), sz1=0.02+track->GetSigmaZ2();
3777 Double_t sy2=0.01+track->GetSigmaY2(), sz2=0.01+track->GetSigmaZ2();
3778 Double_t sy3=0.01+track->GetSigmaY2();
3780 Double_t f40=(F1(x2[0],x2[1]+sy,x1[0],x1[1],x0[0],x0[1])-x[4])/sy;
3781 Double_t f42=(F1(x2[0],x2[1],x1[0],x1[1]+sy,x0[0],x0[1])-x[4])/sy;
3782 Double_t f43=(F1(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]+sy)-x[4])/sy;
3783 Double_t f20=(F2(x2[0],x2[1]+sy,x1[0],x1[1],x0[0],x0[1])-x[2])/sy;
3784 Double_t f22=(F2(x2[0],x2[1],x1[0],x1[1]+sy,x0[0],x0[1])-x[2])/sy;
3785 Double_t f23=(F2(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]+sy)-x[2])/sy;
3787 Double_t f30=(F3(x2[0],x2[1]+sy,x0[0],x0[1],x2[2],x0[2])-x[3])/sy;
3788 Double_t f31=(F3(x2[0],x2[1],x0[0],x0[1],x2[2]+sz,x0[2])-x[3])/sz;
3789 Double_t f32=(F3(x2[0],x2[1],x0[0],x0[1]+sy,x2[2],x0[2])-x[3])/sy;
3790 Double_t f34=(F3(x2[0],x2[1],x0[0],x0[1],x2[2],x0[2]+sz)-x[3])/sz;
3795 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3796 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3797 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3798 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3799 c[13]=f30*sy1*f40+f32*sy2*f42;
3800 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3802 // Int_t row1 = fSectors->GetRowNumber(x2[0]);
3803 AliTPCseed *seed=new AliTPCseed(0, x, c, x2[0], sec2*fSectors->GetAlpha()+fSectors->GetAlphaShift());
3804 // Double_t y0,z0,y1,z1, y2,z2;
3805 //seed->GetProlongation(x0[0],y0,z0);
3806 // seed->GetProlongation(x1[0],y1,z1);
3807 //seed->GetProlongation(x2[0],y2,z2);
3809 seed->fLastPoint = pp2;
3810 seed->fFirstPoint = pp2;
3817 AliTPCseed *AliTPCtrackerMI::ReSeed(AliTPCseed *track, Float_t r0, Float_t r1, Float_t r2)
3821 //reseed using founded clusters
3823 // Find the number of clusters
3824 Int_t nclusters = 0;
3825 for (Int_t irow=0;irow<160;irow++){
3826 if (track->GetClusterIndex(irow)>0) nclusters++;
3830 ipos[0] = TMath::Max(int(r0*nclusters),0); // point 0 cluster
3831 ipos[1] = TMath::Min(int(r1*nclusters),nclusters-1); //
3832 ipos[2] = TMath::Min(int(r2*nclusters),nclusters-1); // last point
3836 Int_t row[3],sec[3]={0,0,0};
3838 // find track row position at given ratio of the length
3840 for (Int_t irow=0;irow<160;irow++){
3841 if (track->GetClusterIndex2(irow)<0) continue;
3843 for (Int_t ipoint=0;ipoint<3;ipoint++){
3844 if (index<=ipos[ipoint]) row[ipoint] = irow;
3848 //Get cluster and sector position
3849 for (Int_t ipoint=0;ipoint<3;ipoint++){
3850 Int_t clindex = track->GetClusterIndex2(row[ipoint]);
3851 AliTPCclusterMI * cl = GetClusterMI(clindex);
3854 // AliTPCclusterMI * cl = GetClusterMI(clindex);
3857 sec[ipoint] = ((clindex&0xff000000)>>24)%18;
3858 xyz[ipoint][0] = GetXrow(row[ipoint]);
3859 xyz[ipoint][1] = cl->GetY();
3860 xyz[ipoint][2] = cl->GetZ();
3864 // Calculate seed state vector and covariance matrix
3866 Double_t alpha, cs,sn, xx2,yy2;
3868 alpha = (sec[1]-sec[2])*fSectors->GetAlpha();
3869 cs = TMath::Cos(alpha);
3870 sn = TMath::Sin(alpha);
3871 xx2= xyz[1][0]*cs-xyz[1][1]*sn;
3872 yy2= xyz[1][0]*sn+xyz[1][1]*cs;
3876 alpha = (sec[0]-sec[2])*fSectors->GetAlpha();
3877 cs = TMath::Cos(alpha);
3878 sn = TMath::Sin(alpha);
3879 xx2= xyz[0][0]*cs-xyz[0][1]*sn;
3880 yy2= xyz[0][0]*sn+xyz[0][1]*cs;
3886 Double_t x[5],c[15];
3890 x[4]=F1(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
3891 x[2]=F2(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
3892 x[3]=F3n(xyz[2][0],xyz[2][1],xyz[0][0],xyz[0][1],xyz[2][2],xyz[0][2],x[4]);
3894 Double_t sy =0.1, sz =0.1;
3896 Double_t sy1=0.2, sz1=0.2;
3897 Double_t sy2=0.2, sz2=0.2;
3900 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;
3901 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;
3902 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;
3903 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;
3904 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;
3905 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;
3907 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;
3908 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;
3909 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;
3910 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;
3915 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3916 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3917 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3918 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3919 c[13]=f30*sy1*f40+f32*sy2*f42;
3920 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3922 // Int_t row1 = fSectors->GetRowNumber(xyz[2][0]);
3923 AliTPCseed *seed=new AliTPCseed(0, x, c, xyz[2][0], sec[2]*fSectors->GetAlpha()+fSectors->GetAlphaShift());
3924 seed->fLastPoint = row[2];
3925 seed->fFirstPoint = row[2];
3930 AliTPCseed *AliTPCtrackerMI::ReSeed(AliTPCseed *track,Int_t r0, Bool_t forward)
3934 //reseed using founded clusters
3937 Int_t row[3]={0,0,0},sec[3]={0,0,0};
3939 // forward direction
3941 for (Int_t irow=r0;irow<160;irow++){
3942 if (track->GetClusterIndex(irow)>0){
3947 for (Int_t irow=160;irow>r0;irow--){
3948 if (track->GetClusterIndex(irow)>0){
3953 for (Int_t irow=row[2]-15;irow>row[0];irow--){
3954 if (track->GetClusterIndex(irow)>0){
3962 for (Int_t irow=0;irow<r0;irow++){
3963 if (track->GetClusterIndex(irow)>0){
3968 for (Int_t irow=r0;irow>0;irow--){
3969 if (track->GetClusterIndex(irow)>0){
3974 for (Int_t irow=row[2]-15;irow>row[0];irow--){
3975 if (track->GetClusterIndex(irow)>0){
3982 if ((row[2]-row[0])<20) return 0;
3983 if (row[1]==0) return 0;
3986 //Get cluster and sector position
3987 for (Int_t ipoint=0;ipoint<3;ipoint++){
3988 Int_t clindex = track->GetClusterIndex2(row[ipoint]);
3989 AliTPCclusterMI * cl = GetClusterMI(clindex);
3992 // AliTPCclusterMI * cl = GetClusterMI(clindex);
3995 sec[ipoint] = ((clindex&0xff000000)>>24)%18;
3996 xyz[ipoint][0] = GetXrow(row[ipoint]);
3997 AliTPCTrackerPoint * point = track->GetTrackPoint(row[ipoint]);
3998 if (point&&ipoint<2){
4000 xyz[ipoint][1] = point->GetY();
4001 xyz[ipoint][2] = point->GetZ();
4004 xyz[ipoint][1] = cl->GetY();
4005 xyz[ipoint][2] = cl->GetZ();
4012 // Calculate seed state vector and covariance matrix
4014 Double_t alpha, cs,sn, xx2,yy2;
4016 alpha = (sec[1]-sec[2])*fSectors->GetAlpha();
4017 cs = TMath::Cos(alpha);
4018 sn = TMath::Sin(alpha);
4019 xx2= xyz[1][0]*cs-xyz[1][1]*sn;
4020 yy2= xyz[1][0]*sn+xyz[1][1]*cs;
4024 alpha = (sec[0]-sec[2])*fSectors->GetAlpha();
4025 cs = TMath::Cos(alpha);
4026 sn = TMath::Sin(alpha);
4027 xx2= xyz[0][0]*cs-xyz[0][1]*sn;
4028 yy2= xyz[0][0]*sn+xyz[0][1]*cs;
4034 Double_t x[5],c[15];
4038 x[4]=F1(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
4039 x[2]=F2(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
4040 x[3]=F3n(xyz[2][0],xyz[2][1],xyz[0][0],xyz[0][1],xyz[2][2],xyz[0][2],x[4]);
4042 Double_t sy =0.1, sz =0.1;
4044 Double_t sy1=0.2, sz1=0.2;
4045 Double_t sy2=0.2, sz2=0.2;
4048 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;
4049 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;
4050 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;
4051 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;
4052 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;
4053 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;
4055 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;
4056 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;
4057 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;
4058 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;
4063 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
4064 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
4065 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
4066 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
4067 c[13]=f30*sy1*f40+f32*sy2*f42;
4068 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
4070 // Int_t row1 = fSectors->GetRowNumber(xyz[2][0]);
4071 AliTPCseed *seed=new AliTPCseed(0, x, c, xyz[2][0], sec[2]*fSectors->GetAlpha()+fSectors->GetAlphaShift());
4072 seed->fLastPoint = row[2];
4073 seed->fFirstPoint = row[2];
4074 for (Int_t i=row[0];i<row[2];i++){
4075 seed->fIndex[i] = track->fIndex[i];
4081 void AliTPCtrackerMI::FindKinks(TObjArray * array, AliESD *esd)
4088 TObjArray *kinks= new TObjArray(10000);
4089 Int_t nentries = array->GetEntriesFast();
4090 AliHelix *helixes = new AliHelix[nentries];
4091 Int_t *sign = new Int_t[nentries];
4092 Int_t *nclusters = new Int_t[nentries];
4093 Float_t *alpha = new Float_t[nentries];
4094 AliESDkink * kink = new AliESDkink();
4095 Int_t * usage = new Int_t[nentries];
4096 Float_t *zm = new Float_t[nentries];
4097 Float_t *z0 = new Float_t[nentries];
4098 Float_t *fim = new Float_t[nentries];
4099 Float_t *shared = new Float_t[nentries];
4100 Bool_t *circular = new Bool_t[nentries];
4102 // nentries = array->GetEntriesFast();
4107 for (Int_t i=0;i<nentries;i++){
4110 AliTPCseed* track = (AliTPCseed*)array->At(i);
4111 if (!track) continue;
4113 track->UpdatePoints();
4114 if (( track->GetPoints()[2]- track->GetPoints()[0])>5 && track->GetPoints()[3]>0.8){
4116 nclusters[i]=track->GetNumberOfClusters();
4117 alpha[i] = track->GetAlpha();
4118 new (&helixes[i]) AliHelix(*track);
4120 helixes[i].Evaluate(0,xyz);
4121 sign[i] = (track->GetC()>0) ? -1:1;
4124 if (track->GetProlongation(x,y,z)){
4126 fim[i] = alpha[i]+TMath::ATan2(y,x);
4129 zm[i] = track->GetZ();
4133 circular[i]= kFALSE;
4134 if (track->GetProlongation(0,y,z)) z0[i] = TMath::Abs(z);
4140 Int_t ncandidates =0;
4143 Double_t phase[2][2],radius[2];
4146 // Find circling track
4147 // TTreeSRedirector cstream("circling.root");
4149 for (Int_t i0=0;i0<nentries;i0++){
4150 AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
4151 if (!track0) continue;
4152 if (track0->fN<40) continue;
4153 if (TMath::Abs(1./track0->fP4)>200) continue;
4154 for (Int_t i1=i0+1;i1<nentries;i1++){
4155 AliTPCseed * track1 = (AliTPCseed*)array->At(i1);
4156 if (!track1) continue;
4157 if (TMath::Abs(1./track1->fP4)>200) continue;
4158 if (track1->fN<40) continue;
4159 if (z0[i0]<20&&z0[i1]<20) continue;
4160 if (track1->fP4*track0->fP4>0) continue;
4161 if (track1->fP3*track0->fP3>0) continue;
4162 if (track0->fBConstrain&&track1->fBConstrain) continue;
4163 if (max(TMath::Abs(1./track0->fP4),TMath::Abs(1./track1->fP4))>190) continue;
4164 if (track0->fBConstrain&&TMath::Abs(track1->fP4)<TMath::Abs(track0->fP4)) continue; //returning - lower momenta
4165 if (track1->fBConstrain&&TMath::Abs(track0->fP4)<TMath::Abs(track1->fP4)) continue; //returning - lower momenta
4167 if ( TMath::Abs(track1->fP3+track0->fP3)>0.1) continue;
4168 Float_t xc0 = helixes[i0].GetHelix(6);
4169 Float_t yc0 = helixes[i0].GetHelix(7);
4170 Float_t r0 = helixes[i0].GetHelix(8);
4171 Float_t xc1 = helixes[i1].GetHelix(6);
4172 Float_t yc1 = helixes[i1].GetHelix(7);
4173 Float_t r1 = helixes[i1].GetHelix(8);
4175 Float_t rmean = (r0+r1)*0.5;
4176 Float_t delta =TMath::Sqrt((xc1-xc0)*(xc1-xc0)+(yc1-yc0)*(yc1-yc0));
4177 if (delta>30) continue;
4178 if (delta>rmean*0.25) continue;
4179 if (TMath::Abs(r0-r1)/rmean>0.3) continue;
4181 Int_t npoints = helixes[i0].GetRPHIintersections(helixes[i1], phase, radius,10);
4182 if (npoints==0) continue;
4183 helixes[i0].GetClosestPhases(helixes[i1], phase);
4187 Double_t hangles[3];
4188 helixes[i0].Evaluate(phase[0][0],xyz0);
4189 helixes[i1].Evaluate(phase[0][1],xyz1);
4191 helixes[i0].GetAngle(phase[0][0],helixes[i1],phase[0][1],hangles);
4192 Double_t deltah[2],deltabest;
4193 if (hangles[2]<2.8) continue;
4195 cstream<<"C"<<track0->fLab<<track1->fLab<<
4196 track0->fP3<<track1->fP3<<
4197 track0->fP4<<track1->fP4<<
4198 delta<<rmean<<npoints<<
4199 hangles[0]<<hangles[2]<<
4200 xyz0[2]<<xyz1[2]<<radius[0]<<"\n";
4204 helixes[i0].ParabolicDCA(helixes[i1],phase[0][0],phase[0][1],radius[0],deltah[0],5);
4206 helixes[i0].ParabolicDCA(helixes[i1],phase[1][0],phase[1][1],radius[1],deltah[1],5);
4207 if (deltah[1]<deltah[0]) ibest=1;
4209 deltabest = TMath::Sqrt(deltah[ibest]);
4210 helixes[i0].Evaluate(phase[ibest][0],xyz0);
4211 helixes[i1].Evaluate(phase[ibest][1],xyz1);
4212 helixes[i0].GetAngle(phase[ibest][0],helixes[i1],phase[ibest][1],hangles);
4213 //Double_t radiusbest = TMath::Sqrt(radius[ibest]);
4217 cstream<<"CC1"<<track0->fLab<<track1->fLab<< //0-1
4218 track0->fP3<<track1->fP3<< //2-3
4219 track0->fP4<<track1->fP4<< //4-5
4220 delta<<rmean<<npoints<< //6-8
4221 hangles[0]<<hangles[2]<< //9-10
4222 xyz0[2]<<xyz1[2]<<radiusbest<<deltabest<< //11--14
4223 track0->fFirstPoint<<track1->fFirstPoint<< //15-16
4224 track0->fLastPoint<<track1->fLastPoint<< //17-18
4225 track0<<track1<< //19-20
4226 phase[ibest][0]<<phase[ibest][1]<<"\n"; //21-22
4229 Bool_t sign =kFALSE;
4230 if (deltabest<3&&hangles[2]>3.06) sign =kTRUE;
4231 if (TMath::Min(TMath::Abs(track0->GetD()),TMath::Abs(track1->GetD()))>10&&deltabest<6&&hangles[2]>3.) sign= kTRUE;
4233 circular[i0] = kTRUE;
4234 circular[i1] = kTRUE;
4243 for (Int_t i =0;i<nentries;i++){
4244 if (sign[i]==0) continue;
4245 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
4248 Double_t cradius0 = 40*40;
4249 Double_t cradius1 = 270*270;
4252 Double_t cdist3=0.55;
4253 for (Int_t j =i+1;j<nentries;j++){
4255 if (sign[j]*sign[i]<1) continue;
4256 if ( (nclusters[i]+nclusters[j])>200) continue;
4257 if ( (nclusters[i]+nclusters[j])<80) continue;
4258 if ( TMath::Abs(zm[i]-zm[j])>60.) continue;
4259 if ( TMath::Abs(fim[i]-fim[j])>0.6 && TMath::Abs(fim[i]-fim[j])<5.7 ) continue;
4260 //AliTPCseed * track1 = (AliTPCseed*)array->At(j); Double_t phase[2][2],radius[2];
4261 Int_t npoints = helixes[i].GetRPHIintersections(helixes[j], phase, radius,20);
4262 if (npoints<1) continue;
4265 if (radius[0]<cradius0||radius[0]>cradius1) continue;
4268 if ( (radius[0]<cradius0||radius[0]>cradius1) && (radius[1]<cradius0||radius[1]>cradius1) ) continue;
4271 Double_t delta1=10000,delta2=10000;
4272 // cuts on the intersection radius
4273 helixes[i].LinearDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta1);
4274 if (radius[0]<20&&delta1<1) continue; //intersection at vertex
4275 if (radius[0]<10&&delta1<3) continue; //intersection at vertex
4277 helixes[i].LinearDCA(helixes[j],phase[1][0],phase[1][1],radius[1],delta2);
4278 if (radius[1]<20&&delta2<1) continue; //intersection at vertex
4279 if (radius[1]<10&&delta2<3) continue; //intersection at vertex
4282 Double_t distance1 = TMath::Min(delta1,delta2);
4283 if (distance1>cdist1) continue; // cut on DCA linear approximation
4285 npoints = helixes[i].GetRPHIintersections(helixes[j], phase, radius,20);
4286 helixes[i].ParabolicDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta1);
4287 if (radius[0]<20&&delta1<1) continue; //intersection at vertex
4288 if (radius[0]<10&&delta1<3) continue; //intersection at vertex
4291 helixes[i].ParabolicDCA(helixes[j],phase[1][0],phase[1][1],radius[1],delta2);
4292 if (radius[1]<20&&delta2<1) continue; //intersection at vertex
4293 if (radius[1]<10&&delta2<3) continue; //intersection at vertex
4295 distance1 = TMath::Min(delta1,delta2);
4298 rkink = TMath::Sqrt(radius[0]);
4301 rkink = TMath::Sqrt(radius[1]);
4303 if (distance1>cdist2) continue;
4306 AliTPCseed * track1 = (AliTPCseed*)array->At(j);
4309 Int_t row0 = GetRowNumber(rkink);
4310 if (row0<10) continue;
4311 if (row0>150) continue;
4314 Float_t dens00=-1,dens01=-1;
4315 Float_t dens10=-1,dens11=-1;
4317 Int_t found,foundable,shared;
4318 track0->GetClusterStatistic(0,row0-5, found, foundable,shared,kFALSE);
4319 if (foundable>5) dens00 = Float_t(found)/Float_t(foundable);
4320 track0->GetClusterStatistic(row0+5,155, found, foundable,shared,kFALSE);
4321 if (foundable>5) dens01 = Float_t(found)/Float_t(foundable);
4323 track1->GetClusterStatistic(0,row0-5, found, foundable,shared,kFALSE);
4324 if (foundable>10) dens10 = Float_t(found)/Float_t(foundable);
4325 track1->GetClusterStatistic(row0+5,155, found, foundable,shared,kFALSE);
4326 if (foundable>10) dens11 = Float_t(found)/Float_t(foundable);
4328 if (dens00<dens10 && dens01<dens11) continue;
4329 if (dens00>dens10 && dens01>dens11) continue;
4330 if (TMath::Max(dens00,dens10)<0.1) continue;
4331 if (TMath::Max(dens01,dens11)<0.3) continue;
4333 if (TMath::Min(dens00,dens10)>0.6) continue;
4334 if (TMath::Min(dens01,dens11)>0.6) continue;
4337 AliTPCseed * ktrack0, *ktrack1;
4346 if (TMath::Abs(ktrack0->GetC())>5) continue; // cut on the curvature for mother particle
4347 AliExternalTrackParam paramm(*ktrack0);
4348 AliExternalTrackParam paramd(*ktrack1);
4349 if (row0>60&&ktrack1->GetReference().X()>90.) new (¶md) AliExternalTrackParam(ktrack1->GetReference());
4352 kink->SetMother(paramm);
4353 kink->SetDaughter(paramd);
4356 Float_t x[3] = { kink->GetPosition()[0],kink->GetPosition()[1],kink->GetPosition()[2]};
4358 fParam->Transform0to1(x,index);
4359 fParam->Transform1to2(x,index);
4360 row0 = GetRowNumber(x[0]);
4362 if (kink->GetR()<100) continue;
4363 if (kink->GetR()>240) continue;
4364 if (kink->GetPosition()[2]/kink->GetR()>AliTPCReconstructor::GetCtgRange()) continue; //out of fiducial volume
4365 if (kink->GetDistance()>cdist3) continue;
4366 Float_t dird = kink->GetDaughterP()[0]*kink->GetPosition()[0]+kink->GetDaughterP()[1]*kink->GetPosition()[1]; // rough direction estimate
4367 if (dird<0) continue;
4369 Float_t dirm = kink->GetMotherP()[0]*kink->GetPosition()[0]+kink->GetMotherP()[1]*kink->GetPosition()[1]; // rough direction estimate
4370 if (dirm<0) continue;
4371 Float_t mpt = TMath::Sqrt(kink->GetMotherP()[0]*kink->GetMotherP()[0]+kink->GetMotherP()[1]*kink->GetMotherP()[1]);
4372 if (mpt<0.2) continue;
4375 //for high momenta momentum not defined well in first iteration
4376 Double_t qt = TMath::Sin(kink->GetAngle(2))*ktrack1->P();
4377 if (qt>0.35) continue;
4380 kink->SetLabel(CookLabel(ktrack0,0.4,0,row0),0);
4381 kink->SetLabel(CookLabel(ktrack1,0.4,row0,160),1);
4383 kink->SetTPCDensity(dens00,0,0);
4384 kink->SetTPCDensity(dens01,0,1);
4385 kink->SetTPCDensity(dens10,1,0);
4386 kink->SetTPCDensity(dens11,1,1);
4387 kink->SetIndex(i,0);
4388 kink->SetIndex(j,1);
4391 kink->SetTPCDensity(dens10,0,0);
4392 kink->SetTPCDensity(dens11,0,1);
4393 kink->SetTPCDensity(dens00,1,0);
4394 kink->SetTPCDensity(dens01,1,1);
4395 kink->SetIndex(j,0);
4396 kink->SetIndex(i,1);
4399 if (mpt<1||kink->GetAngle(2)>0.1){
4400 // angle and densities not defined yet
4401 if (kink->GetTPCDensityFactor()<0.8) continue;
4402 if ((2-kink->GetTPCDensityFactor())*kink->GetDistance() >0.25) continue;
4403 if (kink->GetAngle(2)*ktrack0->P()<0.003) continue; //too small angle
4404 if (kink->GetAngle(2)>0.2&&kink->GetTPCDensityFactor()<1.15) continue;
4405 if (kink->GetAngle(2)>0.2&&kink->GetTPCDensity(0,1)>0.05) continue;
4407 Float_t criticalangle = track0->fC22+track0->fC33;
4408 criticalangle+= track1->fC22+track1->fC33;
4409 criticalangle= 3*TMath::Sqrt(criticalangle);
4410 if (criticalangle>0.02) criticalangle=0.02;
4411 if (kink->GetAngle(2)<criticalangle) continue;
4414 Int_t drow = Int_t(2.+0.5/(0.05+kink->GetAngle(2))); // overlap region defined
4415 Float_t shapesum =0;
4417 for ( Int_t row = row0-drow; row<row0+drow;row++){
4418 if (row<0) continue;
4419 if (row>155) continue;
4420 if (ktrack0->fClusterPointer[row]){
4421 AliTPCTrackerPoint *point =ktrack0->GetTrackPoint(row);
4422 shapesum+=point->GetSigmaY()+point->GetSigmaZ();
4425 if (ktrack1->fClusterPointer[row]){
4426 AliTPCTrackerPoint *point =ktrack1->GetTrackPoint(row);
4427 shapesum+=point->GetSigmaY()+point->GetSigmaZ();
4432 kink->SetShapeFactor(-1.);
4435 kink->SetShapeFactor(shapesum/sum);
4437 // esd->AddKink(kink);
4438 kinks->AddLast(kink);
4439 kink = new AliESDkink;
4444 // sort the kinks according quality - and refit them towards vertex
4446 Int_t nkinks = kinks->GetEntriesFast();
4447 Float_t *quality = new Float_t[nkinks];
4448 Int_t *indexes = new Int_t[nkinks];
4449 AliTPCseed *mothers = new AliTPCseed[nkinks];
4450 AliTPCseed *daughters = new AliTPCseed[nkinks];
4453 for (Int_t i=0;i<nkinks;i++){
4455 AliESDkink *kink = (AliESDkink*)kinks->At(i);
4457 // refit kinks towards vertex
4459 Int_t index0 = kink->GetIndex(0);
4460 Int_t index1 = kink->GetIndex(1);
4461 AliTPCseed * ktrack0 = (AliTPCseed*)array->At(index0);
4462 AliTPCseed * ktrack1 = (AliTPCseed*)array->At(index1);
4464 Int_t sumn=ktrack0->fN+ktrack1->fN;
4466 // Refit Kink under if too small angle
4468 if (kink->GetAngle(2)<0.05){
4469 kink->SetTPCRow0(GetRowNumber(kink->GetR()));
4470 Int_t row0 = kink->GetTPCRow0();
4471 Int_t drow = Int_t(2.+0.5/(0.05+kink->GetAngle(2)));
4474 Int_t last = row0-drow;
4475 if (last<40) last=40;
4476 if (last<ktrack0->fFirstPoint+25) last = ktrack0->fFirstPoint+25;
4477 AliTPCseed* seed0 = ReSeed(ktrack0,last,kFALSE);
4480 Int_t first = row0+drow;
4481 if (first>130) first=130;
4482 if (first>ktrack1->fLastPoint-25) first = TMath::Max(ktrack1->fLastPoint-25,30);
4483 AliTPCseed* seed1 = ReSeed(ktrack1,first,kTRUE);
4485 if (seed0 && seed1){
4486 kink->SetStatus(1,8);
4487 if (RefitKink(*seed0,*seed1,*kink)) kink->SetStatus(1,9);
4488 row0 = GetRowNumber(kink->GetR());
4489 sumn = seed0->fN+seed1->fN;
4490 new (&mothers[i]) AliTPCseed(*seed0);
4491 new (&daughters[i]) AliTPCseed(*seed1);
4494 delete kinks->RemoveAt(i);
4495 if (seed0) delete seed0;
4496 if (seed1) delete seed1;
4499 if (kink->GetDistance()>0.5 || kink->GetR()<110 || kink->GetR()>240) {
4500 delete kinks->RemoveAt(i);
4501 if (seed0) delete seed0;
4502 if (seed1) delete seed1;
4510 if (kink) quality[i] = 160*((0.1+kink->GetDistance())*(2.-kink->GetTPCDensityFactor()))/(sumn+40.); //the longest -clossest will win
4512 TMath::Sort(nkinks,quality,indexes,kFALSE);
4514 //remove double find kinks
4516 for (Int_t ikink0=1;ikink0<nkinks;ikink0++){
4517 AliESDkink * kink0 = (AliESDkink*) kinks->At(indexes[ikink0]);
4518 if (!kink0) continue;
4520 for (Int_t ikink1=0;ikink1<ikink0;ikink1++){
4521 if (!kink0) continue;
4522 AliESDkink * kink1 = (AliESDkink*) kinks->At(indexes[ikink1]);
4523 if (!kink1) continue;
4524 // if not close kink continue
4525 if (TMath::Abs(kink1->GetPosition()[2]-kink0->GetPosition()[2])>10) continue;
4526 if (TMath::Abs(kink1->GetPosition()[1]-kink0->GetPosition()[1])>10) continue;
4527 if (TMath::Abs(kink1->GetPosition()[0]-kink0->GetPosition()[0])>10) continue;
4529 AliTPCseed &mother0 = mothers[indexes[ikink0]];
4530 AliTPCseed &daughter0 = daughters[indexes[ikink0]];
4531 AliTPCseed &mother1 = mothers[indexes[ikink1]];
4532 AliTPCseed &daughter1 = daughters[indexes[ikink1]];
4533 Int_t row0 = (kink0->GetTPCRow0()+kink1->GetTPCRow0())/2;
4542 for (Int_t i=0;i<row0;i++){
4543 if (mother0.fIndex[i]>0 && mother1.fIndex[i]>0){
4546 if (mother0.fIndex[i]==mother1.fIndex[i]){
4553 for (Int_t i=row0;i<158;i++){
4554 if (daughter0.fIndex[i]>0 && daughter0.fIndex[i]>0){
4557 if (mother0.fIndex[i]==mother1.fIndex[i]){
4563 Float_t ratio = Float_t(same+1)/Float_t(both+1);
4564 Float_t ratiom = Float_t(samem+1)/Float_t(bothm+1);
4565 Float_t ratiod = Float_t(samed+1)/Float_t(bothd+1);
4566 if (ratio>0.3 && ratiom>0.5 &&ratiod>0.5) {
4567 Int_t sum0 = mother0.fN+daughter0.fN;
4568 Int_t sum1 = mother1.fN+daughter1.fN;
4570 shared[kink0->GetIndex(0)]= kTRUE;
4571 shared[kink0->GetIndex(1)]= kTRUE;
4572 delete kinks->RemoveAt(indexes[ikink0]);
4575 shared[kink1->GetIndex(0)]= kTRUE;
4576 shared[kink1->GetIndex(1)]= kTRUE;
4577 delete kinks->RemoveAt(indexes[ikink1]);
4584 for (Int_t i=0;i<nkinks;i++){
4585 AliESDkink * kink = (AliESDkink*) kinks->At(indexes[i]);
4586 if (!kink) continue;
4587 kink->SetTPCRow0(GetRowNumber(kink->GetR()));
4588 Int_t index0 = kink->GetIndex(0);
4589 Int_t index1 = kink->GetIndex(1);
4590 if (circular[index0]||circular[index1]&&kink->GetDistance()>0.2) continue;
4591 kink->SetMultiple(usage[index0],0);
4592 kink->SetMultiple(usage[index1],1);
4593 if (kink->GetMultiple()[0]+kink->GetMultiple()[1]>2) continue;
4594 if (kink->GetMultiple()[0]+kink->GetMultiple()[1]>0 && quality[indexes[i]]>0.2) continue;
4595 if (kink->GetMultiple()[0]+kink->GetMultiple()[1]>0 && kink->GetDistance()>0.2) continue;
4596 if (circular[index0]||circular[index1]&&kink->GetDistance()>0.1) continue;
4598 AliTPCseed * ktrack0 = (AliTPCseed*)array->At(index0);
4599 AliTPCseed * ktrack1 = (AliTPCseed*)array->At(index1);
4600 if (!ktrack0 || !ktrack1) continue;
4601 Int_t index = esd->AddKink(kink);
4604 if ( ktrack0->fKinkIndexes[0]==0 && ktrack1->fKinkIndexes[0]==0) { //best kink
4605 if (mothers[indexes[i]].fN>20 && daughters[indexes[i]].fN>20 && (mothers[indexes[i]].fN+daughters[indexes[i]].fN)>100){
4606 new (ktrack0) AliTPCseed(mothers[indexes[i]]);
4607 new (ktrack1) AliTPCseed(daughters[indexes[i]]);
4611 ktrack0->fKinkIndexes[usage[index0]] = -(index+1);
4612 ktrack1->fKinkIndexes[usage[index1]] = (index+1);
4617 // Remove tracks corresponding to shared kink's
4619 for (Int_t i=0;i<nentries;i++){
4620 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
4621 if (!track0) continue;
4622 if (track0->fKinkIndexes[0]!=0) continue;
4623 if (shared[i]) delete array->RemoveAt(i);
4628 RemoveUsed2(array,0.5,0.4,30);
4631 for (Int_t i=0;i<nentries;i++){
4632 AliTPCseed * track0 = (AliTPCseed*)array->At(i);
4633 if (!track0) continue;
4634 if (track0->Pt()<1.4) continue;
4635 //remove double high momenta tracks - overlapped with kink candidates
4638 for (Int_t icl=track0->fFirstPoint;icl<track0->fLastPoint; icl++){
4639 if (track0->fClusterPointer[icl]!=0){
4641 if (track0->fClusterPointer[icl]->IsUsed(10)) shared++;
4644 if (Float_t(shared+1)/Float_t(nall+1)>0.5) {
4645 delete array->RemoveAt(i);
4648 if (track0->fKinkIndexes[0]!=0) continue;
4649 if (track0->GetNumberOfClusters()<80) continue;
4651 AliTPCseed daughter;
4653 if (CheckKinkPoint(track0,mother,daughter, kink)){
4654 if (mother.fN<30||daughter.fN<20) continue; //too short tracks
4655 if (mother.Pt()<1.4) continue;
4656 Int_t row0= kink.GetTPCRow0();
4657 if (kink.GetDistance()>0.5 || kink.GetR()<110. || kink.GetR()>240.) {
4661 Int_t index = esd->AddKink(&kink);
4662 mother.fKinkIndexes[0] = -(index+1);
4663 daughter.fKinkIndexes[0] = index+1;
4665 delete array->RemoveAt(i);
4666 array->AddAt(new AliTPCseed(mother),i);
4669 array->AddLast(new AliTPCseed(mother));
4671 array->AddLast(new AliTPCseed(daughter));
4672 for (Int_t icl=0;icl<row0;icl++) {
4673 if (mother.fClusterPointer[icl]) mother.fClusterPointer[icl]->Use(20);
4676 for (Int_t icl=row0;icl<158;icl++) {
4677 if (daughter.fClusterPointer[icl]) daughter.fClusterPointer[icl]->Use(20);
4683 delete [] daughters;
4704 printf("Ncandidates=\t%d\t%d\t%d\t%d\n",esd->GetNumberOfKinks(),ncandidates,ntracks,nall);
4708 Int_t AliTPCtrackerMI::RefitKink(AliTPCseed &mother, AliTPCseed &daughter, AliESDkink &kink)
4711 // refit kink towards to the vertex
4714 Int_t row0 = GetRowNumber(kink.GetR());
4715 FollowProlongation(mother,0);
4716 mother.Reset(kFALSE);
4718 FollowProlongation(daughter,row0);
4719 daughter.Reset(kFALSE);
4720 FollowBackProlongation(daughter,158);
4721 daughter.Reset(kFALSE);
4722 Int_t first = TMath::Max(row0-20,30);
4723 Int_t last = TMath::Min(row0+20,140);
4725 const Int_t kNdiv =5;
4726 AliTPCseed param0[kNdiv]; // parameters along the track
4727 AliTPCseed param1[kNdiv]; // parameters along the track
4728 AliESDkink kinks[kNdiv]; // corresponding kink parameters
4731 for (Int_t irow=0; irow<kNdiv;irow++){
4732 rows[irow] = first +((last-first)*irow)/(kNdiv-1);
4734 // store parameters along the track
4736 for (Int_t irow=0;irow<kNdiv;irow++){
4737 FollowBackProlongation(mother, rows[irow]);
4738 FollowProlongation(daughter,rows[kNdiv-1-irow]);
4739 new(¶m0[irow]) AliTPCseed(mother);
4740 new(¶m1[kNdiv-1-irow]) AliTPCseed(daughter);
4744 for (Int_t irow=0; irow<kNdiv-1;irow++){
4745 if (param0[irow].fN<kNdiv||param1[irow].fN<kNdiv) continue;
4746 kinks[irow].SetMother(param0[irow]);
4747 kinks[irow].SetDaughter(param1[irow]);
4748 kinks[irow].Update();
4751 // choose kink with best "quality"
4753 Double_t mindist = 10000;
4754 for (Int_t irow=0;irow<kNdiv;irow++){
4755 if (param0[irow].fN<20||param1[irow].fN<20) continue;
4756 if (TMath::Abs(kinks[irow].GetR())>240.) continue;
4757 if (TMath::Abs(kinks[irow].GetR())<100.) continue;
4759 Float_t normdist = TMath::Abs(param0[irow].fX-kinks[irow].GetR())*(0.1+kink.GetDistance());
4760 normdist/= (param0[irow].fN+param1[irow].fN+40.);
4761 if (normdist < mindist){
4767 if (index==-1) return 0;
4770 param0[index].Reset(kTRUE);
4771 FollowProlongation(param0[index],0);
4773 new (&mother) AliTPCseed(param0[index]);
4774 new (&daughter) AliTPCseed(param1[index]); // daughter in vertex
4776 kink.SetMother(mother);
4777 kink.SetDaughter(daughter);
4779 kink.SetTPCRow0(GetRowNumber(kink.GetR()));
4780 kink.SetTPCncls(param0[index].fN,0);
4781 kink.SetTPCncls(param1[index].fN,1);
4782 kink.SetLabel(CookLabel(&mother,0.4, 0,kink.GetTPCRow0()),0);
4783 kink.SetLabel(CookLabel(&daughter,0.4, kink.GetTPCRow0(),160),1);
4784 mother.SetLabel(kink.GetLabel(0));
4785 daughter.SetLabel(kink.GetLabel(1));
4791 void AliTPCtrackerMI::UpdateKinkQualityM(AliTPCseed * seed){
4793 // update Kink quality information for mother after back propagation
4795 if (seed->GetKinkIndex(0)>=0) return;
4796 for (Int_t ikink=0;ikink<3;ikink++){
4797 Int_t index = seed->GetKinkIndex(ikink);
4798 if (index>=0) break;
4799 index = TMath::Abs(index)-1;
4800 AliESDkink * kink = fEvent->GetKink(index);
4801 //kink->fTPCdensity2[0][0]=-1;
4802 //kink->fTPCdensity2[0][1]=-1;
4803 kink->SetTPCDensity2(-1,0,0);
4804 kink->SetTPCDensity2(1,0,1);
4806 Int_t row0 = kink->GetTPCRow0() - 2 - Int_t( 0.5/ (0.05+kink->GetAngle(2)));
4807 if (row0<15) row0=15;
4809 Int_t row1 = kink->GetTPCRow0() + 2 + Int_t( 0.5/ (0.05+kink->GetAngle(2)));
4810 if (row1>145) row1=145;
4812 Int_t found,foundable,shared;
4813 seed->GetClusterStatistic(0,row0, found, foundable,shared,kFALSE);
4814 if (foundable>5) kink->SetTPCDensity2(Float_t(found)/Float_t(foundable),0,0);
4815 seed->GetClusterStatistic(row1,155, found, foundable,shared,kFALSE);
4816 if (foundable>5) kink->SetTPCDensity2(Float_t(found)/Float_t(foundable),0,1);
4821 void AliTPCtrackerMI::UpdateKinkQualityD(AliTPCseed * seed){
4823 // update Kink quality information for daughter after refit
4825 if (seed->GetKinkIndex(0)<=0) return;
4826 for (Int_t ikink=0;ikink<3;ikink++){
4827 Int_t index = seed->GetKinkIndex(ikink);
4828 if (index<=0) break;
4829 index = TMath::Abs(index)-1;
4830 AliESDkink * kink = fEvent->GetKink(index);
4831 kink->SetTPCDensity2(-1,1,0);
4832 kink->SetTPCDensity2(-1,1,1);
4834 Int_t row0 = kink->GetTPCRow0() -2 - Int_t( 0.5/ (0.05+kink->GetAngle(2)));
4835 if (row0<15) row0=15;
4837 Int_t row1 = kink->GetTPCRow0() +2 + Int_t( 0.5/ (0.05+kink->GetAngle(2)));
4838 if (row1>145) row1=145;
4840 Int_t found,foundable,shared;
4841 seed->GetClusterStatistic(0,row0, found, foundable,shared,kFALSE);
4842 if (foundable>5) kink->SetTPCDensity2(Float_t(found)/Float_t(foundable),1,0);
4843 seed->GetClusterStatistic(row1,155, found, foundable,shared,kFALSE);
4844 if (foundable>5) kink->SetTPCDensity2(Float_t(found)/Float_t(foundable),1,1);
4850 Int_t AliTPCtrackerMI::CheckKinkPoint(AliTPCseed*seed,AliTPCseed &mother, AliTPCseed &daughter, AliESDkink &kink)
4853 // check kink point for given track
4854 // if return value=0 kink point not found
4855 // otherwise seed0 correspond to mother particle
4856 // seed1 correspond to daughter particle
4857 // kink parameter of kink point
4859 Int_t middlerow = (seed->fFirstPoint+seed->fLastPoint)/2;
4860 Int_t first = seed->fFirstPoint;
4861 Int_t last = seed->fLastPoint;
4862 if (last-first<20) return 0; // shortest length - 2*30 = 60 pad-rows
4865 AliTPCseed *seed1 = ReSeed(seed,middlerow+20, kTRUE); //middle of chamber
4866 if (!seed1) return 0;
4867 FollowProlongation(*seed1,seed->fLastPoint-20);
4868 seed1->Reset(kTRUE);
4869 FollowProlongation(*seed1,158);
4870 seed1->Reset(kTRUE);
4871 last = seed1->fLastPoint;
4873 AliTPCseed *seed0 = new AliTPCseed(*seed);
4874 seed0->Reset(kFALSE);
4877 AliTPCseed param0[20]; // parameters along the track
4878 AliTPCseed param1[20]; // parameters along the track
4879 AliESDkink kinks[20]; // corresponding kink parameters
4881 for (Int_t irow=0; irow<20;irow++){
4882 rows[irow] = first +((last-first)*irow)/19;
4884 // store parameters along the track
4886 for (Int_t irow=0;irow<20;irow++){
4887 FollowBackProlongation(*seed0, rows[irow]);
4888 FollowProlongation(*seed1,rows[19-irow]);
4889 new(¶m0[irow]) AliTPCseed(*seed0);
4890 new(¶m1[19-irow]) AliTPCseed(*seed1);
4894 for (Int_t irow=0; irow<19;irow++){
4895 kinks[irow].SetMother(param0[irow]);
4896 kinks[irow].SetDaughter(param1[irow]);
4897 kinks[irow].Update();
4900 // choose kink with biggest change of angle
4902 Double_t maxchange= 0;
4903 for (Int_t irow=1;irow<19;irow++){
4904 if (TMath::Abs(kinks[irow].GetR())>240.) continue;
4905 if (TMath::Abs(kinks[irow].GetR())<110.) continue;
4906 Float_t quality = TMath::Abs(kinks[irow].GetAngle(2))/(3.+TMath::Abs(kinks[irow].GetR()-param0[irow].fX));
4907 if ( quality > maxchange){
4908 maxchange = quality;
4915 if (index<0) return 0;
4917 Int_t row0 = GetRowNumber(kinks[index].GetR()); //row 0 estimate
4918 seed0 = new AliTPCseed(param0[index]);
4919 seed1 = new AliTPCseed(param1[index]);
4920 seed0->Reset(kFALSE);
4921 seed1->Reset(kFALSE);
4922 seed0->ResetCovariance();
4923 seed1->ResetCovariance();
4924 FollowProlongation(*seed0,0);
4925 FollowBackProlongation(*seed1,158);
4926 new (&mother) AliTPCseed(*seed0); // backup mother at position 0
4927 seed0->Reset(kFALSE);
4928 seed1->Reset(kFALSE);
4929 seed0->ResetCovariance();
4930 seed1->ResetCovariance();
4932 first = TMath::Max(row0-20,0);
4933 last = TMath::Min(row0+20,158);
4935 for (Int_t irow=0; irow<20;irow++){
4936 rows[irow] = first +((last-first)*irow)/19;
4938 // store parameters along the track
4940 for (Int_t irow=0;irow<20;irow++){
4941 FollowBackProlongation(*seed0, rows[irow]);
4942 FollowProlongation(*seed1,rows[19-irow]);
4943 new(¶m0[irow]) AliTPCseed(*seed0);
4944 new(¶m1[19-irow]) AliTPCseed(*seed1);
4948 for (Int_t irow=0; irow<19;irow++){
4949 kinks[irow].SetMother(param0[irow]);
4950 kinks[irow].SetDaughter(param1[irow]);
4951 // param0[irow].Dump();
4952 //param1[irow].Dump();
4953 kinks[irow].Update();
4956 // choose kink with biggest change of angle
4959 for (Int_t irow=0;irow<20;irow++){
4960 if (TMath::Abs(kinks[irow].GetR())>250.) continue;
4961 if (TMath::Abs(kinks[irow].GetR())<90.) continue;
4962 Float_t quality = TMath::Abs(kinks[irow].GetAngle(2))/(3.+TMath::Abs(kinks[irow].GetR()-param0[irow].fX));
4963 if ( quality > maxchange){
4964 maxchange = quality;
4971 if (index==-1 || param0[index].fN+param1[index].fN<100){
4976 // Float_t anglesigma = TMath::Sqrt(param0[index].fC22+param0[index].fC33+param1[index].fC22+param1[index].fC33);
4978 kink.SetMother(param0[index]);
4979 kink.SetDaughter(param1[index]);
4981 row0 = GetRowNumber(kink.GetR());
4982 kink.SetTPCRow0(row0);
4983 kink.SetLabel(CookLabel(seed0,0.5,0,row0),0);
4984 kink.SetLabel(CookLabel(seed1,0.5,row0,158),1);
4985 kink.SetIndex(-10,0);
4986 kink.SetIndex(int(param0[index].fN+param1[index].fN),1);
4987 kink.SetTPCncls(param0[index].fN,0);
4988 kink.SetTPCncls(param1[index].fN,1);
4991 // new (&mother) AliTPCseed(param0[index]);
4992 new (&daughter) AliTPCseed(param1[index]);
4993 daughter.SetLabel(kink.GetLabel(1));
4994 param0[index].Reset(kTRUE);
4995 FollowProlongation(param0[index],0);
4996 new (&mother) AliTPCseed(param0[index]);
4997 mother.SetLabel(kink.GetLabel(0));
5007 AliTPCseed* AliTPCtrackerMI::ReSeed(AliTPCseed *t)
5010 // reseed - refit - track
5013 // Int_t last = fSectors->GetNRows()-1;
5015 if (fSectors == fOuterSec){
5016 first = TMath::Max(first, t->fFirstPoint-fInnerSec->GetNRows());
5020 first = t->fFirstPoint;
5022 AliTPCseed * seed = MakeSeed(t,0.1,0.5,0.9);
5023 FollowBackProlongation(*t,fSectors->GetNRows()-1);
5025 FollowProlongation(*t,first);
5035 //_____________________________________________________________________________
5036 Int_t AliTPCtrackerMI::ReadSeeds(const TFile *inp) {
5037 //-----------------------------------------------------------------
5038 // This function reades track seeds.
5039 //-----------------------------------------------------------------
5040 TDirectory *savedir=gDirectory;
5042 TFile *in=(TFile*)inp;
5043 if (!in->IsOpen()) {
5044 cerr<<"AliTPCtrackerMI::ReadSeeds(): input file is not open !\n";
5049 TTree *seedTree=(TTree*)in->Get("Seeds");
5051 cerr<<"AliTPCtrackerMI::ReadSeeds(): ";
5052 cerr<<"can't get a tree with track seeds !\n";
5055 AliTPCtrack *seed=new AliTPCtrack;
5056 seedTree->SetBranchAddress("tracks",&seed);
5058 if (fSeeds==0) fSeeds=new TObjArray(15000);
5060 Int_t n=(Int_t)seedTree->GetEntries();
5061 for (Int_t i=0; i<n; i++) {
5062 seedTree->GetEvent(i);
5063 fSeeds->AddLast(new AliTPCseed(*seed/*,seed->GetAlpha()*/));
5072 Int_t AliTPCtrackerMI::Clusters2Tracks (AliESD *esd)
5075 if (fSeeds) DeleteSeeds();
5078 if (!fSeeds) return 1;
5085 //_____________________________________________________________________________
5086 Int_t AliTPCtrackerMI::Clusters2Tracks() {
5087 //-----------------------------------------------------------------
5088 // This is a track finder.
5089 //-----------------------------------------------------------------
5090 TDirectory *savedir=gDirectory;
5094 fSeeds = Tracking();
5097 Info("Clusters2Tracks","Time for tracking: \t");timer.Print();timer.Start();
5099 //activate again some tracks
5100 for (Int_t i=0; i<fSeeds->GetEntriesFast(); i++) {
5101 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
5103 Int_t nc=t.GetNumberOfClusters();
5105 delete fSeeds->RemoveAt(i);
5109 if (pt->fRemoval==10) {
5110 if (pt->GetDensityFirst(20)>0.8 || pt->GetDensityFirst(30)>0.8 || pt->GetDensityFirst(40)>0.7)
5111 pt->Desactivate(10); // make track again active
5113 pt->Desactivate(20);
5114 delete fSeeds->RemoveAt(i);
5119 RemoveUsed2(fSeeds,0.85,0.85,0);
5120 FindKinks(fSeeds,fEvent);
5121 RemoveUsed2(fSeeds,0.5,0.4,50);
5122 //RemoveDouble(fSeeds,0.2,0.6,11);
5123 // RemoveUsed(fSeeds,0.5,0.5,6);
5126 Int_t nseed=fSeeds->GetEntriesFast();
5128 for (Int_t i=0; i<nseed; i++) {
5129 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
5131 Int_t nc=t.GetNumberOfClusters();
5133 delete fSeeds->RemoveAt(i);
5136 CookLabel(pt,0.1); //For comparison only
5137 //if ((pt->IsActive() || (pt->fRemoval==10) )&& nc>50 &&pt->GetNumberOfClusters()>0.4*pt->fNFoundable){
5138 if ((pt->IsActive() || (pt->fRemoval==10) )){
5140 if (fDebug>0) cerr<<found<<'\r';
5144 delete fSeeds->RemoveAt(i);
5148 //RemoveOverlap(fSeeds,0.99,7,kTRUE);
5150 //RemoveUsed(fSeeds,0.9,0.9,6);
5152 nseed=fSeeds->GetEntriesFast();
5154 for (Int_t i=0; i<nseed; i++) {
5155 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
5157 Int_t nc=t.GetNumberOfClusters();
5159 delete fSeeds->RemoveAt(i);
5163 t.CookdEdx(0.02,0.6);
5164 // CheckKinkPoint(&t,0.05);
5165 //if ((pt->IsActive() || (pt->fRemoval==10) )&& nc>50 &&pt->GetNumberOfClusters()>0.4*pt->fNFoundable){
5166 if ((pt->IsActive() || (pt->fRemoval==10) )){
5174 delete fSeeds->RemoveAt(i);
5175 //AliTPCseed * seed1 = ReSeed(pt,0.05,0.5,1);
5177 // FollowProlongation(*seed1,0);
5178 // Int_t n = seed1->GetNumberOfClusters();
5179 // printf("fP4\t%f\t%f\n",seed1->GetC(),pt->GetC());
5180 // printf("fN\t%d\t%d\n", seed1->GetNumberOfClusters(),pt->GetNumberOfClusters());
5183 //AliTPCseed * seed2 = ReSeed(pt,0.95,0.5,0.05);
5187 SortTracks(fSeeds, 1);
5191 PrepareForBackProlongation(fSeeds,5.);
5192 PropagateBack(fSeeds);
5193 printf("Time for back propagation: \t");timer.Print();timer.Start();
5197 PrepareForProlongation(fSeeds,5.);
5198 PropagateForward2(fSeeds);
5200 printf("Time for FORWARD propagation: \t");timer.Print();timer.Start();
5201 // RemoveUsed(fSeeds,0.7,0.7,6);
5202 //RemoveOverlap(fSeeds,0.9,7,kTRUE);
5204 nseed=fSeeds->GetEntriesFast();
5206 for (Int_t i=0; i<nseed; i++) {
5207 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
5209 Int_t nc=t.GetNumberOfClusters();
5211 delete fSeeds->RemoveAt(i);
5214 t.CookdEdx(0.02,0.6);
5215 // CookLabel(pt,0.1); //For comparison only
5216 //if ((pt->IsActive() || (pt->fRemoval==10) )&& nc>50 &&pt->GetNumberOfClusters()>0.4*pt->fNFoundable){
5217 if ((pt->IsActive() || (pt->fRemoval==10) )){
5218 cerr<<found++<<'\r';
5221 delete fSeeds->RemoveAt(i);
5226 // fNTracks = found;
5228 Info("Clusters2Tracks","Time for overlap removal, track writing and dedx cooking: \t"); timer.Print();timer.Start();
5231 // cerr<<"Number of found tracks : "<<"\t"<<found<<endl;
5232 Info("Clusters2Tracks","Number of found tracks %d",found);
5234 // UnloadClusters();
5239 void AliTPCtrackerMI::Tracking(TObjArray * arr)
5242 // tracking of the seeds
5245 fSectors = fOuterSec;
5246 ParallelTracking(arr,150,63);
5247 fSectors = fOuterSec;
5248 ParallelTracking(arr,63,0);
5251 TObjArray * AliTPCtrackerMI::Tracking(Int_t seedtype, Int_t i1, Int_t i2, Float_t cuts[4], Float_t dy, Int_t dsec)
5256 TObjArray * arr = new TObjArray;
5258 fSectors = fOuterSec;
5261 for (Int_t sec=0;sec<fkNOS;sec++){
5262 if (seedtype==3) MakeSeeds3(arr,sec,i1,i2,cuts,dy, dsec);
5263 if (seedtype==4) MakeSeeds5(arr,sec,i1,i2,cuts,dy);
5264 if (seedtype==2) MakeSeeds2(arr,sec,i1,i2,cuts,dy);
5267 Info("Tracking","\nSeeding - %d\t%d\t%d\t%d\n",seedtype,i1,i2,arr->GetEntriesFast());
5279 TObjArray * AliTPCtrackerMI::Tracking()
5285 Int_t nup=fOuterSec->GetNRows()+fInnerSec->GetNRows();
5287 TObjArray * seeds = new TObjArray;
5296 Float_t fnumber = 3.0;
5297 Float_t fdensity = 3.0;
5302 for (Int_t delta = 0; delta<18; delta+=6){
5306 arr = Tracking(3,nup-1-delta,nup-1-delta-gap,cuts,-1,1);
5307 SumTracks(seeds,arr);
5308 SignClusters(seeds,fnumber,fdensity);
5310 for (Int_t i=2;i<6;i+=2){
5311 // seed high pt tracks
5314 arr = Tracking(3,nup-i-delta,nup-i-delta-gap,cuts,-1,0);
5315 SumTracks(seeds,arr);
5316 SignClusters(seeds,fnumber,fdensity);
5321 // RemoveUsed(seeds,0.9,0.9,1);
5322 // UnsignClusters();
5323 // SignClusters(seeds,fnumber,fdensity);
5327 for (Int_t delta = 20; delta<120; delta+=10){
5329 // seed high pt tracks
5333 arr = Tracking(3,nup-delta,nup-delta-gap,cuts,-1);
5334 SumTracks(seeds,arr);
5335 SignClusters(seeds,fnumber,fdensity);
5340 arr = Tracking(3,nup-delta-5,nup-delta-5-gap,cuts,-1);
5341 SumTracks(seeds,arr);
5342 SignClusters(seeds,fnumber,fdensity);
5353 Info("Tracking()","\n\nPrimary seeding\t%d\n\n",seeds->GetEntriesFast());
5357 // RemoveUsed(seeds,0.75,0.75,1);
5359 //SignClusters(seeds,fnumber,fdensity);
5368 arr = Tracking(4,nup-1,nup-1-gap,cuts,-1);
5369 SumTracks(seeds,arr);
5370 SignClusters(seeds,fnumber,fdensity);
5372 arr = Tracking(4,nup-2,nup-2-gap,cuts,-1);
5373 SumTracks(seeds,arr);
5374 SignClusters(seeds,fnumber,fdensity);
5376 arr = Tracking(4,nup-3,nup-3-gap,cuts,-1);
5377 SumTracks(seeds,arr);
5378 SignClusters(seeds,fnumber,fdensity);
5382 for (Int_t delta = 3; delta<30; delta+=5){
5388 arr = Tracking(4,nup-1-delta,nup-1-delta-gap,cuts,-1);
5389 SumTracks(seeds,arr);
5390 SignClusters(seeds,fnumber,fdensity);
5392 arr = Tracking(4,nup-3-delta,nup-5-delta-gap,cuts,4);
5393 SumTracks(seeds,arr);
5394 SignClusters(seeds,fnumber,fdensity);
5406 for (Int_t delta = 30; delta<70; delta+=10){
5412 arr = Tracking(4,nup-1-delta,nup-1-delta-gap,cuts,-1);
5413 SumTracks(seeds,arr);
5414 SignClusters(seeds,fnumber,fdensity);
5416 arr = Tracking(4,nup-5-delta,nup-5-delta-gap,cuts,5 );
5417 SumTracks(seeds,arr);
5418 SignClusters(seeds,fnumber,fdensity);
5422 Info("Tracking()","\n\nSecondary seeding\t%d\n\n",seeds->GetEntriesFast());
5433 void AliTPCtrackerMI::SumTracks(TObjArray *arr1,TObjArray *arr2) const
5436 //sum tracks to common container
5437 //remove suspicious tracks
5438 Int_t nseed = arr2->GetEntriesFast();
5439 for (Int_t i=0;i<nseed;i++){
5440 AliTPCseed *pt=(AliTPCseed*)arr2->UncheckedAt(i);
5443 // NORMAL ACTIVE TRACK
5444 if (pt->IsActive()){
5445 arr1->AddLast(arr2->RemoveAt(i));
5448 //remove not usable tracks
5449 if (pt->fRemoval!=10){
5450 delete arr2->RemoveAt(i);
5453 // REMOVE VERY SHORT TRACKS
5454 if (pt->GetNumberOfClusters()<20){
5455 delete arr2->RemoveAt(i);
5458 // ENABLE ONLY ENOUGH GOOD STOPPED TRACKS
5459 if (pt->GetDensityFirst(20)>0.8 || pt->GetDensityFirst(30)>0.8 || pt->GetDensityFirst(40)>0.7)
5460 arr1->AddLast(arr2->RemoveAt(i));
5462 delete arr2->RemoveAt(i);
5471 void AliTPCtrackerMI::ParallelTracking(TObjArray * arr, Int_t rfirst, Int_t rlast)
5474 // try to track in parralel
5476 Int_t nseed=arr->GetEntriesFast();
5477 //prepare seeds for tracking
5478 for (Int_t i=0; i<nseed; i++) {
5479 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i), &t=*pt;
5481 if (!t.IsActive()) continue;
5482 // follow prolongation to the first layer
5483 if ( (fSectors ==fInnerSec) || (t.fFirstPoint-fParam->GetNRowLow()>rfirst+1) )
5484 FollowProlongation(t, rfirst+1);
5489 for (Int_t nr=rfirst; nr>=rlast; nr--){
5490 if (nr<fInnerSec->GetNRows())
5491 fSectors = fInnerSec;
5493 fSectors = fOuterSec;
5494 // make indexes with the cluster tracks for given
5496 // find nearest cluster
5497 for (Int_t i=0; i<nseed; i++) {
5498 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i), &t=*pt;
5500 if (nr==80) pt->UpdateReference();
5501 if (!pt->IsActive()) continue;
5502 // if ( (fSectors ==fOuterSec) && (pt->fFirstPoint-fParam->GetNRowLow())<nr) continue;
5503 if (pt->fRelativeSector>17) {
5506 UpdateClusters(t,nr);
5508 // prolonagate to the nearest cluster - if founded
5509 for (Int_t i=0; i<nseed; i++) {
5510 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
5512 if (!pt->IsActive()) continue;
5513 // if ((fSectors ==fOuterSec) && (pt->fFirstPoint-fParam->GetNRowLow())<nr) continue;
5514 if (pt->fRelativeSector>17) {
5517 FollowToNextCluster(*pt,nr);
5522 void AliTPCtrackerMI::PrepareForBackProlongation(TObjArray * arr,Float_t fac) const
5526 // if we use TPC track itself we have to "update" covariance
5528 Int_t nseed= arr->GetEntriesFast();
5529 for (Int_t i=0;i<nseed;i++){
5530 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
5534 //rotate to current local system at first accepted point
5535 Int_t index = pt->GetClusterIndex2(pt->fFirstPoint);
5536 Int_t sec = (index&0xff000000)>>24;
5538 Float_t angle1 = fInnerSec->GetAlpha()*sec+fInnerSec->GetAlphaShift();
5539 if (angle1>TMath::Pi())
5540 angle1-=2.*TMath::Pi();
5541 Float_t angle2 = pt->GetAlpha();
5543 if (TMath::Abs(angle1-angle2)>0.001){
5544 pt->Rotate(angle1-angle2);
5545 //angle2 = pt->GetAlpha();
5546 //pt->fRelativeSector = pt->GetAlpha()/fInnerSec->GetAlpha();
5547 //if (pt->GetAlpha()<0)
5548 // pt->fRelativeSector+=18;
5549 //sec = pt->fRelativeSector;
5558 void AliTPCtrackerMI::PrepareForProlongation(TObjArray * arr, Float_t fac) const
5562 // if we use TPC track itself we have to "update" covariance
5564 Int_t nseed= arr->GetEntriesFast();
5565 for (Int_t i=0;i<nseed;i++){
5566 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
5569 pt->fFirstPoint = pt->fLastPoint;
5577 Int_t AliTPCtrackerMI::PropagateBack(TObjArray * arr)
5580 // make back propagation
5582 Int_t nseed= arr->GetEntriesFast();
5583 for (Int_t i=0;i<nseed;i++){
5584 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
5585 if (pt&& pt->GetKinkIndex(0)<=0) {
5586 //AliTPCseed *pt2 = new AliTPCseed(*pt);
5587 fSectors = fInnerSec;
5588 //FollowBackProlongation(*pt,fInnerSec->GetNRows()-1);
5589 //fSectors = fOuterSec;
5590 FollowBackProlongation(*pt,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1);
5591 //if (pt->GetNumberOfClusters()<(pt->fEsd->GetTPCclusters(0)) ){
5592 // Error("PropagateBack","Not prolonged track %d",pt->GetLabel());
5593 // FollowBackProlongation(*pt2,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1);
5596 if (pt&& pt->GetKinkIndex(0)>0) {
5597 AliESDkink * kink = fEvent->GetKink(pt->GetKinkIndex(0)-1);
5598 pt->fFirstPoint = kink->GetTPCRow0();
5599 fSectors = fInnerSec;
5600 FollowBackProlongation(*pt,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1);
5608 Int_t AliTPCtrackerMI::PropagateForward2(TObjArray * arr)
5611 // make forward propagation
5613 Int_t nseed= arr->GetEntriesFast();
5615 for (Int_t i=0;i<nseed;i++){
5616 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
5618 FollowProlongation(*pt,0);
5625 Int_t AliTPCtrackerMI::PropagateForward()
5628 // propagate track forward
5630 Int_t nseed = fSeeds->GetEntriesFast();
5631 for (Int_t i=0;i<nseed;i++){
5632 AliTPCseed *pt = (AliTPCseed*)fSeeds->UncheckedAt(i);
5634 AliTPCseed &t = *pt;
5635 Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
5636 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
5637 if (alpha < 0. ) alpha += 2.*TMath::Pi();
5638 t.fRelativeSector = Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN;
5642 fSectors = fOuterSec;
5643 ParallelTracking(fSeeds,fOuterSec->GetNRows()+fInnerSec->GetNRows()-1,fInnerSec->GetNRows());
5644 fSectors = fInnerSec;
5645 ParallelTracking(fSeeds,fInnerSec->GetNRows()-1,0);
5655 Int_t AliTPCtrackerMI::PropagateBack(AliTPCseed * pt, Int_t row0, Int_t row1)
5658 // make back propagation, in between row0 and row1
5662 fSectors = fInnerSec;
5665 if (row1<fSectors->GetNRows())
5668 r1 = fSectors->GetNRows()-1;
5670 if (row0<fSectors->GetNRows()&& r1>0 )
5671 FollowBackProlongation(*pt,r1);
5672 if (row1<=fSectors->GetNRows())
5675 r1 = row1 - fSectors->GetNRows();
5676 if (r1<=0) return 0;
5677 if (r1>=fOuterSec->GetNRows()) return 0;
5678 fSectors = fOuterSec;
5679 return FollowBackProlongation(*pt,r1);
5687 void AliTPCtrackerMI::GetShape(AliTPCseed * seed, Int_t row)
5691 Float_t sd2 = TMath::Abs((fParam->GetZLength()-TMath::Abs(seed->GetZ())))*fParam->GetDiffL()*fParam->GetDiffL();
5692 // Float_t padlength = fParam->GetPadPitchLength(seed->fSector);
5693 Float_t padlength = GetPadPitchLength(row);
5695 Float_t sresy = (seed->fSector < fParam->GetNSector()/2) ? 0.2 :0.3;
5696 Float_t angulary = seed->GetSnp();
5697 angulary = angulary*angulary/(1-angulary*angulary);
5698 seed->fCurrentSigmaY2 = sd2+padlength*padlength*angulary/12.+sresy*sresy;
5700 Float_t sresz = fParam->GetZSigma();
5701 Float_t angularz = seed->GetTgl();
5702 seed->fCurrentSigmaZ2 = sd2+padlength*padlength*angularz*angularz*(1+angulary)/12.+sresz*sresz;
5704 Float_t wy = GetSigmaY(seed);
5705 Float_t wz = GetSigmaZ(seed);
5708 if (TMath::Abs(wy/seed->fCurrentSigmaY2-1)>0.0001 || TMath::Abs(wz/seed->fCurrentSigmaZ2-1)>0.0001 ){
5709 printf("problem\n");
5715 Float_t AliTPCtrackerMI::GetSigmaY(AliTPCseed * seed)
5719 Float_t sd2 = TMath::Abs((fParam->GetZLength()-TMath::Abs(seed->GetZ())))*fParam->GetDiffL()*fParam->GetDiffL();
5720 Float_t padlength = fParam->GetPadPitchLength(seed->fSector);
5721 Float_t sres = (seed->fSector < fParam->GetNSector()/2) ? 0.2 :0.3;
5722 Float_t angular = seed->GetSnp();
5723 angular = angular*angular/(1-angular*angular);
5724 // angular*=angular;
5725 //angular = TMath::Sqrt(angular/(1-angular));
5726 Float_t res = TMath::Sqrt(sd2+padlength*padlength*angular/12.+sres*sres);
5729 Float_t AliTPCtrackerMI::GetSigmaZ(AliTPCseed * seed)
5733 Float_t sd2 = TMath::Abs((fParam->GetZLength()-TMath::Abs(seed->GetZ())))*fParam->GetDiffL()*fParam->GetDiffL();
5734 Float_t padlength = fParam->GetPadPitchLength(seed->fSector);
5735 Float_t sres = fParam->GetZSigma();
5736 Float_t angular = seed->GetTgl();
5737 Float_t res = TMath::Sqrt(sd2+padlength*padlength*angular*angular/12.+sres*sres);
5742 //__________________________________________________________________________
5743 void AliTPCtrackerMI::CookLabel(AliTPCseed *t, Float_t wrong) const {
5744 //--------------------------------------------------------------------
5745 //This function "cooks" a track label. If label<0, this track is fake.
5746 //--------------------------------------------------------------------
5747 Int_t noc=t->GetNumberOfClusters();
5749 //printf("\nnot founded prolongation\n\n\n");
5755 AliTPCclusterMI *clusters[160];
5757 for (Int_t i=0;i<160;i++) {
5764 for (i=0; i<160 && current<noc; i++) {
5766 Int_t index=t->GetClusterIndex2(i);
5767 if (index<=0) continue;
5768 if (index&0x8000) continue;
5770 //clusters[current]=GetClusterMI(index);
5771 if (t->fClusterPointer[i]){
5772 clusters[current]=t->fClusterPointer[i];
5778 Int_t lab=123456789;
5779 for (i=0; i<noc; i++) {
5780 AliTPCclusterMI *c=clusters[i];
5782 lab=TMath::Abs(c->GetLabel(0));
5784 for (j=0; j<noc; j++) if (lb[j]==lab || mx[j]==0) break;
5790 for (i=0; i<noc; i++) if (mx[i]>max) {max=mx[i]; lab=lb[i];}
5792 for (i=0; i<noc; i++) {
5793 AliTPCclusterMI *c=clusters[i];
5795 if (TMath::Abs(c->GetLabel(1)) == lab ||
5796 TMath::Abs(c->GetLabel(2)) == lab ) max++;
5799 if ((1.- Float_t(max)/noc) > wrong) lab=-lab;
5802 Int_t tail=Int_t(0.10*noc);
5805 for (i=1; i<=160&&ind<tail; i++) {
5806 // AliTPCclusterMI *c=clusters[noc-i];
5807 AliTPCclusterMI *c=clusters[i];
5809 if (lab == TMath::Abs(c->GetLabel(0)) ||
5810 lab == TMath::Abs(c->GetLabel(1)) ||
5811 lab == TMath::Abs(c->GetLabel(2))) max++;
5814 if (max < Int_t(0.5*tail)) lab=-lab;
5821 //delete[] clusters;
5825 //__________________________________________________________________________
5826 Int_t AliTPCtrackerMI::CookLabel(AliTPCseed *t, Float_t wrong,Int_t first, Int_t last) const {
5827 //--------------------------------------------------------------------
5828 //This function "cooks" a track label. If label<0, this track is fake.
5829 //--------------------------------------------------------------------
5830 Int_t noc=t->GetNumberOfClusters();
5832 //printf("\nnot founded prolongation\n\n\n");
5838 AliTPCclusterMI *clusters[160];
5840 for (Int_t i=0;i<160;i++) {
5847 for (i=0; i<160 && current<noc; i++) {
5848 if (i<first) continue;
5849 if (i>last) continue;
5850 Int_t index=t->GetClusterIndex2(i);
5851 if (index<=0) continue;
5852 if (index&0x8000) continue;
5854 //clusters[current]=GetClusterMI(index);
5855 if (t->fClusterPointer[i]){
5856 clusters[current]=t->fClusterPointer[i];
5861 if (noc<5) return -1;
5862 Int_t lab=123456789;
5863 for (i=0; i<noc; i++) {
5864 AliTPCclusterMI *c=clusters[i];
5866 lab=TMath::Abs(c->GetLabel(0));
5868 for (j=0; j<noc; j++) if (lb[j]==lab || mx[j]==0) break;
5874 for (i=0; i<noc; i++) if (mx[i]>max) {max=mx[i]; lab=lb[i];}
5876 for (i=0; i<noc; i++) {
5877 AliTPCclusterMI *c=clusters[i];
5879 if (TMath::Abs(c->GetLabel(1)) == lab ||
5880 TMath::Abs(c->GetLabel(2)) == lab ) max++;
5883 if ((1.- Float_t(max)/noc) > wrong) lab=-lab;
5886 Int_t tail=Int_t(0.10*noc);
5889 for (i=1; i<=160&&ind<tail; i++) {
5890 // AliTPCclusterMI *c=clusters[noc-i];
5891 AliTPCclusterMI *c=clusters[i];
5893 if (lab == TMath::Abs(c->GetLabel(0)) ||
5894 lab == TMath::Abs(c->GetLabel(1)) ||
5895 lab == TMath::Abs(c->GetLabel(2))) max++;
5898 if (max < Int_t(0.5*tail)) lab=-lab;
5901 // t->SetLabel(lab);
5905 //delete[] clusters;
5909 Int_t AliTPCtrackerMI::AliTPCSector::GetRowNumber(Double_t x) const
5911 //return pad row number for this x
5914 r=fRow[fN-1].GetX();
5915 if (x > r) return fN;
5917 if (x < r) return -1;
5918 return Int_t((x-r)/fPadPitchLength + 0.5);}
5920 r=fRow[fN-1].GetX();
5921 if (x > r) return fN;
5923 if (x < r) return -1;
5924 Double_t r1=fRow[64].GetX();
5926 return Int_t((x-r)/f1PadPitchLength + 0.5);}
5928 return (Int_t((x-r1)/f2PadPitchLength + 0.5)+64);}
5932 Int_t AliTPCtrackerMI::GetRowNumber(Double_t x[3]) const
5934 //return pad row number for given x vector
5935 Float_t phi = TMath::ATan2(x[1],x[0]);
5936 if(phi<0) phi=2.*TMath::Pi()+phi;
5937 // Get the local angle in the sector philoc
5938 const Float_t kRaddeg = 180/3.14159265358979312;
5939 Float_t phiangle = (Int_t (phi*kRaddeg/20.) + 0.5)*20./kRaddeg;
5940 Double_t localx = x[0]*TMath::Cos(phiangle)-x[1]*TMath::Sin(phiangle);
5941 return GetRowNumber(localx);
5944 //_________________________________________________________________________
5945 void AliTPCtrackerMI::AliTPCSector::Setup(const AliTPCParam *par, Int_t f) {
5946 //-----------------------------------------------------------------------
5947 // Setup inner sector
5948 //-----------------------------------------------------------------------
5950 fAlpha=par->GetInnerAngle();
5951 fAlphaShift=par->GetInnerAngleShift();
5952 fPadPitchWidth=par->GetInnerPadPitchWidth();
5953 fPadPitchLength=par->GetInnerPadPitchLength();
5954 fN=par->GetNRowLow();
5955 fRow=new AliTPCRow[fN];
5956 for (Int_t i=0; i<fN; i++) {
5957 fRow[i].SetX(par->GetPadRowRadiiLow(i));
5958 fRow[i].fDeadZone =1.5; //1.5 cm of dead zone
5961 fAlpha=par->GetOuterAngle();
5962 fAlphaShift=par->GetOuterAngleShift();
5963 fPadPitchWidth = par->GetOuterPadPitchWidth();
5964 fPadPitchLength = par->GetOuter1PadPitchLength();
5965 f1PadPitchLength = par->GetOuter1PadPitchLength();
5966 f2PadPitchLength = par->GetOuter2PadPitchLength();
5968 fN=par->GetNRowUp();
5969 fRow=new AliTPCRow[fN];
5970 for (Int_t i=0; i<fN; i++) {
5971 fRow[i].SetX(par->GetPadRowRadiiUp(i));
5972 fRow[i].fDeadZone =1.5; // 1.5 cm of dead zone
5977 AliTPCtrackerMI::AliTPCRow::AliTPCRow() {
5979 // default constructor
5987 AliTPCtrackerMI::AliTPCRow::~AliTPCRow(){
5994 //_________________________________________________________________________
5996 AliTPCtrackerMI::AliTPCRow::InsertCluster(const AliTPCclusterMI* c, UInt_t index) {
5997 //-----------------------------------------------------------------------
5998 // Insert a cluster into this pad row in accordence with its y-coordinate
5999 //-----------------------------------------------------------------------
6000 if (fN==kMaxClusterPerRow) {
6001 cerr<<"AliTPCRow::InsertCluster(): Too many clusters !\n"; return;
6003 if (fN==0) {fIndex[0]=index; fClusters[fN++]=c; return;}
6004 Int_t i=Find(c->GetZ());
6005 memmove(fClusters+i+1 ,fClusters+i,(fN-i)*sizeof(AliTPCclusterMI*));
6006 memmove(fIndex +i+1 ,fIndex +i,(fN-i)*sizeof(UInt_t));
6007 fIndex[i]=index; fClusters[i]=c; fN++;
6010 void AliTPCtrackerMI::AliTPCRow::ResetClusters() {
6016 //delete[] fClusterArray;
6017 if (fClusters1) delete []fClusters1;
6018 if (fClusters2) delete []fClusters2;
6025 //___________________________________________________________________
6026 Int_t AliTPCtrackerMI::AliTPCRow::Find(Double_t z) const {
6027 //-----------------------------------------------------------------------
6028 // Return the index of the nearest cluster
6029 //-----------------------------------------------------------------------
6030 if (fN==0) return 0;
6031 if (z <= fClusters[0]->GetZ()) return 0;
6032 if (z > fClusters[fN-1]->GetZ()) return fN;
6033 Int_t b=0, e=fN-1, m=(b+e)/2;
6034 for (; b<e; m=(b+e)/2) {
6035 if (z > fClusters[m]->GetZ()) b=m+1;
6043 //___________________________________________________________________
6044 AliTPCclusterMI * AliTPCtrackerMI::AliTPCRow::FindNearest(Double_t y, Double_t z, Double_t roady, Double_t roadz) const {
6045 //-----------------------------------------------------------------------
6046 // Return the index of the nearest cluster in z y
6047 //-----------------------------------------------------------------------
6048 Float_t maxdistance = roady*roady + roadz*roadz;
6050 AliTPCclusterMI *cl =0;
6051 for (Int_t i=Find(z-roadz); i<fN; i++) {
6052 AliTPCclusterMI *c=(AliTPCclusterMI*)(fClusters[i]);
6053 if (c->GetZ() > z+roadz) break;
6054 if ( (c->GetY()-y) > roady ) continue;
6055 Float_t distance = (c->GetZ()-z)*(c->GetZ()-z)+(c->GetY()-y)*(c->GetY()-y);
6056 if (maxdistance>distance) {
6057 maxdistance = distance;
6064 AliTPCclusterMI * AliTPCtrackerMI::AliTPCRow::FindNearest2(Double_t y, Double_t z, Double_t roady, Double_t roadz,UInt_t & index) const
6066 //-----------------------------------------------------------------------
6067 // Return the index of the nearest cluster in z y
6068 //-----------------------------------------------------------------------
6069 Float_t maxdistance = roady*roady + roadz*roadz;
6070 Int_t iz1 = TMath::Max(fFastCluster[Int_t(z-roadz+254.5)]-1,0);
6071 Int_t iz2 = TMath::Min(fFastCluster[Int_t(z+roadz+255.5)]+1,fN);
6073 AliTPCclusterMI *cl =0;
6074 //FindNearest3(y,z,roady,roadz,index);
6075 // for (Int_t i=Find(z-roadz); i<fN; i++) {
6076 for (Int_t i=iz1; i<iz2; i++) {
6077 AliTPCclusterMI *c=(AliTPCclusterMI*)(fClusters[i]);
6078 if (c->GetZ() > z+roadz) break;
6079 if ( c->GetY()-y > roady ) continue;
6080 if ( y-c->GetY() > roady ) continue;
6081 Float_t distance = (c->GetZ()-z)*(c->GetZ()-z)+(c->GetY()-y)*(c->GetY()-y);
6082 if (maxdistance>distance) {
6083 maxdistance = distance;
6086 //roady = TMath::Sqrt(maxdistance);
6094 AliTPCclusterMI * AliTPCtrackerMI::AliTPCRow::FindNearest3(Double_t y, Double_t z, Double_t roady, Double_t roadz,UInt_t & index) const
6096 //-----------------------------------------------------------------------
6097 // Return the index of the nearest cluster in z y
6098 //-----------------------------------------------------------------------
6099 Float_t maxdistance = roady*roady + roadz*roadz;
6100 // Int_t iz = Int_t(z+255.);
6101 AliTPCclusterMI *cl =0;
6102 for (Int_t i=Find(z-roadz); i<fN; i++) {
6103 //for (Int_t i=fFastCluster[iz-2]; i<fFastCluster[iz+2]; i++) {
6104 AliTPCclusterMI *c=(AliTPCclusterMI*)(fClusters[i]);
6105 if (c->GetZ() > z+roadz) break;
6106 if ( c->GetY()-y > roady ) continue;
6107 if ( y-c->GetY() > roady ) continue;
6108 Float_t distance = (c->GetZ()-z)*(c->GetZ()-z)+(c->GetY()-y)*(c->GetY()-y);
6109 if (maxdistance>distance) {
6110 maxdistance = distance;
6113 //roady = TMath::Sqrt(maxdistance);
6122 AliTPCseed::AliTPCseed():AliTPCtrack(){
6126 for (Int_t i=0;i<200;i++) SetClusterIndex2(i,-3);
6127 for (Int_t i=0;i<160;i++) fClusterPointer[i]=0;
6128 for (Int_t i=0;i<3;i++) fKinkIndexes[i]=0;
6144 AliTPCseed::AliTPCseed(const AliTPCseed &s):AliTPCtrack(s){
6145 //---------------------
6146 // dummy copy constructor
6147 //-------------------------
6148 for (Int_t i=0;i<160;i++) fClusterPointer[i] = s.fClusterPointer[i];
6149 for (Int_t i=0;i<160;i++) fIndex[i] = s.fIndex[i];
6153 AliTPCseed::AliTPCseed(const AliTPCtrack &t):AliTPCtrack(t){
6162 for (Int_t i=0;i<3;i++) fKinkIndexes[i]=t.GetKinkIndex(i);
6164 for (Int_t i=0;i<160;i++) {
6165 fClusterPointer[i] = 0;
6166 Int_t index = t.GetClusterIndex(i);
6168 SetClusterIndex2(i,index);
6171 SetClusterIndex2(i,-3);
6184 AliTPCseed::AliTPCseed(const AliTPCtrack &t, Double_t a):AliTPCtrack(t,a){
6188 for (Int_t i=0;i<160;i++) {
6189 fClusterPointer[i] = 0;
6190 Int_t index = t.GetClusterIndex(i);
6191 SetClusterIndex2(i,index);
6193 for (Int_t i=0;i<3;i++) fKinkIndexes[i]=0;
6213 AliTPCseed::AliTPCseed(UInt_t index, const Double_t xx[5], const Double_t cc[15],
6214 Double_t xr, Double_t alpha):
6215 AliTPCtrack(index, xx, cc, xr, alpha) {
6220 for (Int_t i=0;i<200;i++) SetClusterIndex2(i,-3);
6221 for (Int_t i=0;i<160;i++) fClusterPointer[i]=0;
6222 for (Int_t i=0;i<3;i++) fKinkIndexes[i]=0;
6231 // fHelixIn = new TClonesArray("AliHelix",0);
6232 //fHelixOut = new TClonesArray("AliHelix",0);
6242 AliTPCseed::~AliTPCseed(){
6245 if (fPoints) delete fPoints;
6247 if (fEPoints) delete fEPoints;
6252 AliTPCTrackerPoint * AliTPCseed::GetTrackPoint(Int_t i)
6256 return &fTrackPoints[i];
6259 void AliTPCseed::RebuildSeed()
6262 // rebuild seed to be ready for storing
6263 AliTPCclusterMI cldummy;
6265 AliTPCTrackPoint pdummy;
6266 pdummy.GetTPoint().fIsShared = 10;
6267 for (Int_t i=0;i<160;i++){
6268 AliTPCclusterMI * cl0 = fClusterPointer[i];
6269 AliTPCTrackPoint *trpoint = (AliTPCTrackPoint*)fPoints->UncheckedAt(i);
6271 trpoint->GetTPoint() = *(GetTrackPoint(i));
6272 trpoint->GetCPoint() = *cl0;
6273 trpoint->GetCPoint().SetQ(TMath::Abs(cl0->GetQ()));
6277 trpoint->GetCPoint()= cldummy;
6285 Double_t AliTPCseed::GetDensityFirst(Int_t n)
6289 // return cluster for n rows bellow first point
6290 Int_t nfoundable = 1;
6292 for (Int_t i=fLastPoint-1;i>0&&nfoundable<n; i--){
6293 Int_t index = GetClusterIndex2(i);
6294 if (index!=-1) nfoundable++;
6295 if (index>0) nfound++;
6297 if (nfoundable<n) return 0;
6298 return Double_t(nfound)/Double_t(nfoundable);
6303 void AliTPCseed::GetClusterStatistic(Int_t first, Int_t last, Int_t &found, Int_t &foundable, Int_t &shared, Bool_t plus2)
6305 // get cluster stat. on given region
6310 for (Int_t i=first;i<last; i++){
6311 Int_t index = GetClusterIndex2(i);
6312 if (index!=-1) foundable++;
6313 if (fClusterPointer[i]) {
6319 if (fClusterPointer[i]->IsUsed(10)) {
6323 if (!plus2) continue; //take also neighborhoud
6325 if ( (i>0) && fClusterPointer[i-1]){
6326 if (fClusterPointer[i-1]->IsUsed(10)) {
6331 if ( fClusterPointer[i+1]){
6332 if (fClusterPointer[i+1]->IsUsed(10)) {
6339 //if (shared>found){
6340 //Error("AliTPCseed::GetClusterStatistic","problem\n");
6344 //_____________________________________________________________________________
6345 void AliTPCseed::CookdEdx(Double_t low, Double_t up,Int_t i1, Int_t i2, Bool_t onlyused) {
6346 //-----------------------------------------------------------------
6347 // This funtion calculates dE/dX within the "low" and "up" cuts.
6348 //-----------------------------------------------------------------
6351 Float_t angular[200];
6352 Float_t weight[200];
6355 // TClonesArray & arr = *fPoints;
6356 Float_t meanlog = 100.;
6358 Float_t mean[4] = {0,0,0,0};
6359 Float_t sigma[4] = {1000,1000,1000,1000};
6360 Int_t nc[4] = {0,0,0,0};
6361 Float_t norm[4] = {1000,1000,1000,1000};
6366 for (Int_t of =0; of<4; of++){
6367 for (Int_t i=of+i1;i<i2;i+=4)
6369 Int_t index = fIndex[i];
6370 if (index<0||index&0x8000) continue;
6372 //AliTPCTrackPoint * point = (AliTPCTrackPoint *) arr.At(i);
6373 AliTPCTrackerPoint * point = GetTrackPoint(i);
6374 //AliTPCTrackerPoint * pointm = GetTrackPoint(i-1);
6375 //AliTPCTrackerPoint * pointp = 0;
6376 //if (i<159) pointp = GetTrackPoint(i+1);
6378 if (point==0) continue;
6379 AliTPCclusterMI * cl = fClusterPointer[i];
6380 if (cl==0) continue;
6381 if (onlyused && (!cl->IsUsed(10))) continue;
6382 if (cl->IsUsed(11)) {
6386 Int_t type = cl->GetType();
6387 //if (point->fIsShared){
6392 // if (pointm->fIsShared) continue;
6394 // if (pointp->fIsShared) continue;
6396 if (type<0) continue;
6397 //if (type>10) continue;
6398 //if (point->GetErrY()==0) continue;
6399 //if (point->GetErrZ()==0) continue;
6401 //Float_t ddy = (point->GetY()-cl->GetY())/point->GetErrY();
6402 //Float_t ddz = (point->GetZ()-cl->GetZ())/point->GetErrZ();
6403 //if ((ddy*ddy+ddz*ddz)>10) continue;
6406 // if (point->GetCPoint().GetMax()<5) continue;
6407 if (cl->GetMax()<5) continue;
6408 Float_t angley = point->GetAngleY();
6409 Float_t anglez = point->GetAngleZ();
6411 Float_t rsigmay2 = point->GetSigmaY();
6412 Float_t rsigmaz2 = point->GetSigmaZ();
6416 rsigmay += pointm->GetTPoint().GetSigmaY();
6417 rsigmaz += pointm->GetTPoint().GetSigmaZ();
6421 rsigmay += pointp->GetTPoint().GetSigmaY();
6422 rsigmaz += pointp->GetTPoint().GetSigmaZ();
6429 Float_t rsigma = TMath::Sqrt(rsigmay2*rsigmaz2);
6431 Float_t ampc = 0; // normalization to the number of electrons
6433 // ampc = 1.*point->GetCPoint().GetMax();
6434 ampc = 1.*cl->GetMax();
6435 //ampc = 1.*point->GetCPoint().GetQ();
6436 // AliTPCClusterPoint & p = point->GetCPoint();
6437 // Float_t dy = TMath::Abs(Int_t( TMath::Abs(p.GetY()/0.6)) - TMath::Abs(p.GetY()/0.6)+0.5);
6438 // Float_t iz = (250.0-TMath::Abs(p.GetZ())+0.11)/0.566;
6440 // TMath::Abs( Int_t(iz) - iz + 0.5);
6441 //ampc *= 1.15*(1-0.3*dy);
6442 //ampc *= 1.15*(1-0.3*dz);
6443 // Float_t zfactor = (AliTPCReconstructor::GetCtgRange()-0.0004*TMath::Abs(point->GetCPoint().GetZ()));
6447 //ampc = 1.0*point->GetCPoint().GetMax();
6448 ampc = 1.0*cl->GetMax();
6449 //ampc = 1.0*point->GetCPoint().GetQ();
6450 //AliTPCClusterPoint & p = point->GetCPoint();
6451 // Float_t dy = TMath::Abs(Int_t( TMath::Abs(p.GetY()/0.4)) - TMath::Abs(p.GetY()/0.4)+0.5);
6452 //Float_t iz = (250.0-TMath::Abs(p.GetZ())+0.11)/0.566;
6454 // TMath::Abs( Int_t(iz) - iz + 0.5);
6456 //ampc *= 1.15*(1-0.3*dy);
6457 //ampc *= 1.15*(1-0.3*dz);
6458 // Float_t zfactor = (1.02-0.000*TMath::Abs(point->GetCPoint().GetZ()));
6462 ampc *= 2.0; // put mean value to channel 50
6463 //ampc *= 0.58; // put mean value to channel 50
6465 // if (type>0) w = 1./(type/2.-0.5);
6466 // Float_t z = TMath::Abs(cl->GetZ());
6469 //ampc /= (1+0.0008*z);
6473 //ampc /= (1+0.0008*z);
6475 //ampc /= (1+0.0008*z);
6478 if (type<0) { //amp at the border - lower weight
6483 if (rsigma>1.5) ampc/=1.3; // if big backround
6485 angular[nc[of]] = TMath::Sqrt(1.+angley*angley+anglez*anglez);
6490 TMath::Sort(nc[of],amp,index,kFALSE);
6494 //meanlog = amp[index[Int_t(nc[of]*0.33)]];
6496 for (Int_t i=int(nc[of]*low+0.5);i<int(nc[of]*up+0.5);i++){
6497 Float_t ampl = amp[index[i]]/angular[index[i]];
6498 ampl = meanlog*TMath::Log(1.+ampl/meanlog);
6500 sumw += weight[index[i]];
6501 sumamp += weight[index[i]]*ampl;
6502 sumamp2 += weight[index[i]]*ampl*ampl;
6503 norm[of] += angular[index[i]]*weight[index[i]];
6510 mean[of] = sumamp/sumw;
6511 sigma[of] = sumamp2/sumw-mean[of]*mean[of];
6513 sigma[of] = TMath::Sqrt(sigma[of]);
6517 mean[of] = (TMath::Exp(mean[of]/meanlog)-1)*meanlog;
6518 //mean *=(1-0.02*(sigma/(mean*0.17)-1.));
6519 //mean *=(1-0.1*(norm-1.));
6526 // mean[0]*= (1-0.05*(sigma[0]/(0.01+mean[1]*0.18)-1));
6527 // mean[1]*= (1-0.05*(sigma[1]/(0.01+mean[0]*0.18)-1));
6530 // dedx = (mean[0]* TMath::Sqrt((1.+nc[0]))+ mean[1]* TMath::Sqrt((1.+nc[1])) )/
6531 // ( TMath::Sqrt((1.+nc[0]))+TMath::Sqrt((1.+nc[1])));
6535 for (Int_t i =0;i<4;i++){
6536 if (nc[i]>2&&nc[i]<1000){
6537 dedx += mean[i] *nc[i];
6538 fSdEdx += sigma[i]*(nc[i]-2);
6539 fMAngular += norm[i] *nc[i];
6544 fSDEDX[i] = sigma[i];
6557 // Float_t dedx1 =dedx;
6560 for (Int_t i =0;i<4;i++){
6561 if (nc[i]>2&&nc[i]<1000){
6562 mean[i] = mean[i]*(1-0.12*(sigma[i]/(fSdEdx)-1.));
6563 dedx += mean[i] *nc[i];
6578 Double_t p=TMath::Sqrt((1.+ GetTgl()*GetTgl())/(Get1Pt()*Get1Pt()));
6581 if (dedx < 39.+ 12./(p+0.25)/(p+0.25)) { SetMass(0.13957); return;}
6582 if (dedx < 39.+ 12./p/p) { SetMass(0.49368); return;}
6583 SetMass(0.93827); return;
6587 if (dedx < 39.+ 12./(p+0.25)/(p+0.25)) { SetMass(0.13957); return;}
6588 SetMass(0.93827); return;
6591 SetMass(0.13957); return;
6601 void AliTPCseed::CookdEdx2(Double_t low, Double_t up) {
6602 //-----------------------------------------------------------------
6603 // This funtion calculates dE/dX within the "low" and "up" cuts.
6604 //-----------------------------------------------------------------
6607 Float_t angular[200];
6608 Float_t weight[200];
6610 Bool_t inlimit[200];
6611 for (Int_t i=0;i<200;i++) inlimit[i]=kFALSE;
6612 for (Int_t i=0;i<200;i++) amp[i]=10000;
6613 for (Int_t i=0;i<200;i++) angular[i]= 1;;
6617 Float_t meanlog = 100.;
6618 Int_t indexde[4]={0,64,128,160};
6625 Float_t mean[4] = {0,0,0,0};
6626 Float_t sigma[4] = {1000,1000,1000,1000};
6627 Int_t nc[4] = {0,0,0,0};
6628 Float_t norm[4] = {1000,1000,1000,1000};
6633 // for (Int_t of =0; of<3; of++){
6634 // for (Int_t i=indexde[of];i<indexde[of+1];i++)
6635 for (Int_t i =0; i<160;i++)
6637 AliTPCTrackPoint * point = GetTrackPoint(i);
6638 if (point==0) continue;
6639 if (point->fIsShared){
6643 Int_t type = point->GetCPoint().GetType();
6644 if (type<0) continue;
6645 if (point->GetCPoint().GetMax()<5) continue;
6646 Float_t angley = point->GetTPoint().GetAngleY();
6647 Float_t anglez = point->GetTPoint().GetAngleZ();
6648 Float_t rsigmay = point->GetCPoint().GetSigmaY();
6649 Float_t rsigmaz = point->GetCPoint().GetSigmaZ();
6650 Float_t rsigma = TMath::Sqrt(rsigmay*rsigmaz);
6652 Float_t ampc = 0; // normalization to the number of electrons
6654 ampc = point->GetCPoint().GetMax();
6657 ampc = point->GetCPoint().GetMax();
6659 ampc *= 2.0; // put mean value to channel 50
6660 // ampc *= 0.565; // put mean value to channel 50
6663 Float_t z = TMath::Abs(point->GetCPoint().GetZ());
6670 if (type<0) { //amp at the border - lower weight
6673 if (rsigma>1.5) ampc/=1.3; // if big backround
6674 angular[i] = TMath::Sqrt(1.+angley*angley+anglez*anglez);
6675 amp[i] = ampc/angular[i];
6680 TMath::Sort(159,amp,index,kFALSE);
6681 for (Int_t i=int(anc*low+0.5);i<int(anc*up+0.5);i++){
6682 inlimit[index[i]] = kTRUE; // take all clusters
6685 // meanlog = amp[index[Int_t(anc*0.3)]];
6687 for (Int_t of =0; of<3; of++){
6691 for (Int_t i=indexde[of];i<indexde[of+1];i++)
6693 if (inlimit[i]==kFALSE) continue;
6694 Float_t ampl = amp[i];
6696 ampl = meanlog*TMath::Log(1.+ampl/meanlog);
6699 sumamp += weight[i]*ampl;
6700 sumamp2 += weight[i]*ampl*ampl;
6701 norm[of] += angular[i]*weight[i];
6709 mean[of] = sumamp/sumw;
6710 sigma[of] = sumamp2/sumw-mean[of]*mean[of];
6712 sigma[of] = TMath::Sqrt(sigma[of]);
6715 mean[of] = (TMath::Exp(mean[of]/meanlog)-1)*meanlog;
6725 Float_t www[3] = {12.,14.,17.};
6726 //Float_t www[3] = {1.,1.,1.};
6728 for (Int_t i =0;i<3;i++){
6729 if (nc[i]>2&&nc[i]<1000){
6730 dedx += mean[i] *nc[i]*www[i]/sigma[i];
6731 fSdEdx += sigma[i]*(nc[i]-2)*www[i]/sigma[i];
6732 fMAngular += norm[i] *nc[i];
6733 norm2 += nc[i]*www[i]/sigma[i];
6734 norm3 += (nc[i]-2)*www[i]/sigma[i];
6737 fSDEDX[i] = sigma[i];
6750 // Float_t dedx1 =dedx;
6754 for (Int_t i =0;i<3;i++){
6755 if (nc[i]>2&&nc[i]<1000&&sigma[i]>3){
6756 //mean[i] = mean[i]*(1+0.08*(sigma[i]/(fSdEdx)-1.));
6757 dedx += mean[i] *(nc[i])/(sigma[i]);
6758 norm4 += (nc[i])/(sigma[i]);
6762 if (norm4>0) dedx /= norm4;