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"
55 ClassImp(AliTPCtrackerMI)
58 class AliTPCFastMath {
61 static Double_t FastAsin(Double_t x);
63 static Double_t fgFastAsin[20000]; //lookup table for fast asin computation
66 Double_t AliTPCFastMath::fgFastAsin[20000];
67 AliTPCFastMath gAliTPCFastMath; // needed to fill the LUT
69 AliTPCFastMath::AliTPCFastMath(){
71 // initialized lookup table;
72 for (Int_t i=0;i<10000;i++){
73 fgFastAsin[2*i] = TMath::ASin(i/10000.);
74 fgFastAsin[2*i+1] = (TMath::ASin((i+1)/10000.)-fgFastAsin[2*i]);
78 Double_t AliTPCFastMath::FastAsin(Double_t x){
80 // return asin using lookup table
82 Int_t index = int(x*10000);
83 return fgFastAsin[2*index]+(x*10000.-index)*fgFastAsin[2*index+1];
86 Int_t index = int(x*10000);
87 return -(fgFastAsin[2*index]+(x*10000.-index)*fgFastAsin[2*index+1]);
93 Int_t AliTPCtrackerMI::UpdateTrack(AliTPCseed * track, Int_t accept){
95 //update track information using current cluster - track->fCurrentCluster
98 AliTPCclusterMI* c =track->fCurrentCluster;
99 if (accept>0) track->fCurrentClusterIndex1 |=0x8000; //sign not accepted clusters
101 UInt_t i = track->fCurrentClusterIndex1;
103 Int_t sec=(i&0xff000000)>>24;
104 //Int_t row = (i&0x00ff0000)>>16;
105 track->fRow=(i&0x00ff0000)>>16;
106 track->fSector = sec;
107 // Int_t index = i&0xFFFF;
108 if (sec>=fParam->GetNInnerSector()) track->fRow += fParam->GetNRowLow();
109 track->SetClusterIndex2(track->fRow, i);
110 //track->fFirstPoint = row;
111 //if ( track->fLastPoint<row) track->fLastPoint =row;
112 // if (track->fRow<0 || track->fRow>160) {
113 // printf("problem\n");
115 if (track->fFirstPoint>track->fRow)
116 track->fFirstPoint = track->fRow;
117 if (track->fLastPoint<track->fRow)
118 track->fLastPoint = track->fRow;
121 track->fClusterPointer[track->fRow] = c;
124 Float_t angle2 = track->GetSnp()*track->GetSnp();
125 angle2 = TMath::Sqrt(angle2/(1-angle2));
127 //SET NEW Track Point
131 AliTPCTrackerPoint &point =*(track->GetTrackPoint(track->fRow));
133 point.SetSigmaY(c->GetSigmaY2()/track->fCurrentSigmaY2);
134 point.SetSigmaZ(c->GetSigmaZ2()/track->fCurrentSigmaZ2);
135 point.SetErrY(sqrt(track->fErrorY2));
136 point.SetErrZ(sqrt(track->fErrorZ2));
138 point.SetX(track->GetX());
139 point.SetY(track->GetY());
140 point.SetZ(track->GetZ());
141 point.SetAngleY(angle2);
142 point.SetAngleZ(track->GetTgl());
143 if (point.fIsShared){
144 track->fErrorY2 *= 4;
145 track->fErrorZ2 *= 4;
149 Double_t chi2 = track->GetPredictedChi2(track->fCurrentCluster);
151 track->fErrorY2 *= 1.3;
152 track->fErrorY2 += 0.01;
153 track->fErrorZ2 *= 1.3;
154 track->fErrorZ2 += 0.005;
156 if (accept>0) return 0;
157 if (track->GetNumberOfClusters()%20==0){
158 // if (track->fHelixIn){
159 // TClonesArray & larr = *(track->fHelixIn);
160 // Int_t ihelix = larr.GetEntriesFast();
161 // new(larr[ihelix]) AliHelix(*track) ;
164 track->fNoCluster =0;
165 return track->Update(c,chi2,i);
170 Int_t AliTPCtrackerMI::AcceptCluster(AliTPCseed * seed, AliTPCclusterMI * cluster, Float_t factor,
171 Float_t cory, Float_t corz)
174 // decide according desired precision to accept given
175 // cluster for tracking
176 Double_t sy2=ErrY2(seed,cluster)*cory;
177 Double_t sz2=ErrZ2(seed,cluster)*corz;
178 //sy2=ErrY2(seed,cluster)*cory;
179 //sz2=ErrZ2(seed,cluster)*cory;
181 Double_t sdistancey2 = sy2+seed->GetSigmaY2();
182 Double_t sdistancez2 = sz2+seed->GetSigmaZ2();
184 Double_t rdistancey2 = (seed->fCurrentCluster->GetY()-seed->GetY())*
185 (seed->fCurrentCluster->GetY()-seed->GetY())/sdistancey2;
186 Double_t rdistancez2 = (seed->fCurrentCluster->GetZ()-seed->GetZ())*
187 (seed->fCurrentCluster->GetZ()-seed->GetZ())/sdistancez2;
189 Double_t rdistance2 = rdistancey2+rdistancez2;
192 if (rdistance2>16) return 3;
195 if ((rdistancey2>9.*factor || rdistancez2>9.*factor) && cluster->GetType()==0)
196 return 2; //suspisiouce - will be changed
198 if ((rdistancey2>6.25*factor || rdistancez2>6.25*factor) && cluster->GetType()>0)
199 // strict cut on overlaped cluster
200 return 2; //suspisiouce - will be changed
202 if ( (rdistancey2>1.*factor || rdistancez2>6.25*factor )
203 && cluster->GetType()<0){
213 //_____________________________________________________________________________
214 AliTPCtrackerMI::AliTPCtrackerMI(const AliTPCParam *par):
215 AliTracker(), fkNIS(par->GetNInnerSector()/2), fkNOS(par->GetNOuterSector()/2)
217 //---------------------------------------------------------------------
218 // The main TPC tracker constructor
219 //---------------------------------------------------------------------
220 fInnerSec=new AliTPCSector[fkNIS];
221 fOuterSec=new AliTPCSector[fkNOS];
224 for (i=0; i<fkNIS; i++) fInnerSec[i].Setup(par,0);
225 for (i=0; i<fkNOS; i++) fOuterSec[i].Setup(par,1);
232 Int_t nrowlow = par->GetNRowLow();
233 Int_t nrowup = par->GetNRowUp();
236 for (Int_t i=0;i<nrowlow;i++){
237 fXRow[i] = par->GetPadRowRadiiLow(i);
238 fPadLength[i]= par->GetPadPitchLength(0,i);
239 fYMax[i] = fXRow[i]*TMath::Tan(0.5*par->GetInnerAngle());
243 for (Int_t i=0;i<nrowup;i++){
244 fXRow[i+nrowlow] = par->GetPadRowRadiiUp(i);
245 fPadLength[i+nrowlow] = par->GetPadPitchLength(60,i);
246 fYMax[i+nrowlow] = fXRow[i+nrowlow]*TMath::Tan(0.5*par->GetOuterAngle());
258 //________________________________________________________________________
259 AliTPCtrackerMI::AliTPCtrackerMI(const AliTPCtrackerMI &t):AliTracker(t){
260 //------------------------------------
261 // dummy copy constructor
262 //------------------------------------------------------------------
264 AliTPCtrackerMI & AliTPCtrackerMI::operator=(const AliTPCtrackerMI& /*r*/){
265 //------------------------------
267 //--------------------------------------------------------------
270 //_____________________________________________________________________________
271 AliTPCtrackerMI::~AliTPCtrackerMI() {
272 //------------------------------------------------------------------
273 // TPC tracker destructor
274 //------------------------------------------------------------------
283 void AliTPCtrackerMI::SetIO()
287 fInput = AliRunLoader::GetTreeR("TPC", kFALSE,AliConfig::GetDefaultEventFolderName());
289 fOutput = AliRunLoader::GetTreeT("TPC", kTRUE,AliConfig::GetDefaultEventFolderName());
291 AliTPCtrack *iotrack= new AliTPCtrack;
292 fOutput->Branch("tracks","AliTPCtrack",&iotrack,32000,100);
298 void AliTPCtrackerMI::SetIO(TTree * input, TTree * output, AliESD * event)
314 AliTPCtrack *iotrack= new AliTPCtrack;
315 // iotrack->fHelixIn = new TClonesArray("AliHelix");
316 //iotrack->fHelixOut = new TClonesArray("AliHelix");
317 fOutput->Branch("tracks","AliTPCtrack",&iotrack,32000,100);
320 if (output && (fDebug&2)){
321 //write the full seed information if specified in debug mode
323 fSeedTree = new TTree("Seeds","Seeds");
324 AliTPCseed * vseed = new AliTPCseed;
326 TClonesArray * arrtr = new TClonesArray("AliTPCTrackPoint",160);
327 arrtr->ExpandCreateFast(160);
328 TClonesArray * arre = new TClonesArray("AliTPCExactPoint",160);
330 vseed->fPoints = arrtr;
331 vseed->fEPoints = arre;
332 // vseed->fClusterPoints = arrcl;
333 fSeedTree->Branch("seeds","AliTPCseed",&vseed,32000,99);
336 fTreeDebug = new TTree("trackDebug","trackDebug");
337 TClonesArray * arrd = new TClonesArray("AliTPCTrackPoint2",0);
338 fTreeDebug->Branch("debug",&arrd,32000,99);
346 void AliTPCtrackerMI::FillESD(TObjArray* arr)
350 //fill esds using updated tracks
352 // write tracks to the event
353 // store index of the track
354 Int_t nseed=arr->GetEntriesFast();
355 for (Int_t i=0; i<nseed; i++) {
356 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
358 pt->PropagateTo(fParam->GetInnerRadiusLow());
359 if ( (pt->GetNumberOfClusters()>70)&& (Float_t(pt->GetNumberOfClusters())/Float_t(pt->fNFoundable))>0.55) {
361 iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
362 //iotrack.SetTPCindex(i);
363 fEvent->AddTrack(&iotrack);
369 void AliTPCtrackerMI::WriteTracks(TTree * tree)
372 // write tracks from seed array to selected tree
376 AliTPCtrack *iotrack= new AliTPCtrack;
377 fOutput->Branch("tracks","AliTPCtrack",&iotrack,32000,100);
382 void AliTPCtrackerMI::WriteTracks()
385 // write tracks to the given output tree -
386 // output specified with SetIO routine
393 AliTPCtrack *iotrack= 0;
394 Int_t nseed=fSeeds->GetEntriesFast();
395 //for (Int_t i=0; i<nseed; i++) {
396 // iotrack= (AliTPCtrack*)fSeeds->UncheckedAt(i);
397 // if (iotrack) break;
399 //TBranch * br = fOutput->Branch("tracks","AliTPCtrack",&iotrack,32000,100);
400 TBranch * br = fOutput->GetBranch("tracks");
401 br->SetAddress(&iotrack);
403 for (Int_t i=0; i<nseed; i++) {
404 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i);
406 AliTPCtrack * track = new AliTPCtrack(*pt);
409 // br->SetAddress(&iotrack);
414 //fOutput->GetDirectory()->cd();
420 //write the full seed information if specified in debug mode
422 AliTPCseed * vseed = new AliTPCseed;
424 TClonesArray * arrtr = new TClonesArray("AliTPCTrackPoint",160);
425 arrtr->ExpandCreateFast(160);
426 //TClonesArray * arrcl = new TClonesArray("AliTPCclusterMI",160);
427 //arrcl->ExpandCreateFast(160);
428 TClonesArray * arre = new TClonesArray("AliTPCExactPoint",160);
430 vseed->fPoints = arrtr;
431 vseed->fEPoints = arre;
432 // vseed->fClusterPoints = arrcl;
433 //TBranch * brseed = seedtree->Branch("seeds","AliTPCseed",&vseed,32000,99);
434 TBranch * brseed = fSeedTree->GetBranch("seeds");
436 Int_t nseed=fSeeds->GetEntriesFast();
438 for (Int_t i=0; i<nseed; i++) {
439 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i);
442 // pt->fClusterPoints = arrcl;
446 brseed->SetAddress(&vseed);
450 // pt->fClusterPoints = 0;
453 if (fTreeDebug) fTreeDebug->Write();
461 Double_t AliTPCtrackerMI::ErrY2(AliTPCseed* seed, AliTPCclusterMI * cl){
464 //seed->SetErrorY2(0.1);
466 //calculate look-up table at the beginning
467 static Bool_t ginit = kFALSE;
468 static Float_t gnoise1,gnoise2,gnoise3;
469 static Float_t ggg1[10000];
470 static Float_t ggg2[10000];
471 static Float_t ggg3[10000];
472 static Float_t glandau1[10000];
473 static Float_t glandau2[10000];
474 static Float_t glandau3[10000];
476 static Float_t gcor01[500];
477 static Float_t gcor02[500];
478 static Float_t gcorp[500];
483 for (Int_t i=1;i<500;i++){
484 Float_t rsigma = float(i)/100.;
485 gcor02[i] = TMath::Max(0.78 +TMath::Exp(7.4*(rsigma-1.2)),0.6);
486 gcor01[i] = TMath::Max(0.72 +TMath::Exp(3.36*(rsigma-1.2)),0.6);
487 gcorp[i] = TMath::Max(TMath::Power((rsigma+0.5),1.5),1.2);
491 for (Int_t i=3;i<10000;i++){
495 Float_t amp = float(i);
496 Float_t padlength =0.75;
497 gnoise1 = 0.0004/padlength;
498 Float_t nel = 0.268*amp;
499 Float_t nprim = 0.155*amp;
500 ggg1[i] = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.001*nel/(padlength*padlength))/nel;
501 glandau1[i] = (2.+0.12*nprim)*0.5* (2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
502 if (glandau1[i]>1) glandau1[i]=1;
503 glandau1[i]*=padlength*padlength/12.;
507 gnoise2 = 0.0004/padlength;
510 ggg2[i] = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
511 glandau2[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
512 if (glandau2[i]>1) glandau2[i]=1;
513 glandau2[i]*=padlength*padlength/12.;
518 gnoise3 = 0.0004/padlength;
521 ggg3[i] = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
522 glandau3[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
523 if (glandau3[i]>1) glandau3[i]=1;
524 glandau3[i]*=padlength*padlength/12.;
532 Int_t amp = int(TMath::Abs(cl->GetQ()));
534 seed->SetErrorY2(1.);
538 Float_t z = TMath::Abs(fParam->GetZLength()-TMath::Abs(seed->GetZ()));
539 Int_t ctype = cl->GetType();
540 Float_t padlength= GetPadPitchLength(seed->fRow);
541 Float_t angle2 = seed->GetSnp()*seed->GetSnp();
542 angle2 = angle2/(1-angle2);
545 Int_t rsigmay = int(100.*cl->GetSigmaY2()/(seed->fCurrentSigmaY2));
548 if (fSectors==fInnerSec){
550 res = ggg1[amp]*z+glandau1[amp]*angle2;
551 if (ctype==0) res *= gcor01[rsigmay];
554 res*= gcorp[rsigmay];
560 res = ggg2[amp]*z+glandau2[amp]*angle2;
561 if (ctype==0) res *= gcor02[rsigmay];
564 res*= gcorp[rsigmay];
569 res = ggg3[amp]*z+glandau3[amp]*angle2;
570 if (ctype==0) res *= gcor02[rsigmay];
573 res*= gcorp[rsigmay];
580 res*=2.4; // overestimate error 2 times
587 seed->SetErrorY2(res);
595 Double_t AliTPCtrackerMI::ErrZ2(AliTPCseed* seed, AliTPCclusterMI * cl){
598 //seed->SetErrorY2(0.1);
600 //calculate look-up table at the beginning
601 static Bool_t ginit = kFALSE;
602 static Float_t gnoise1,gnoise2,gnoise3;
603 static Float_t ggg1[10000];
604 static Float_t ggg2[10000];
605 static Float_t ggg3[10000];
606 static Float_t glandau1[10000];
607 static Float_t glandau2[10000];
608 static Float_t glandau3[10000];
610 static Float_t gcor01[1000];
611 static Float_t gcor02[1000];
612 static Float_t gcorp[1000];
617 for (Int_t i=1;i<1000;i++){
618 Float_t rsigma = float(i)/100.;
619 gcor02[i] = TMath::Max(0.81 +TMath::Exp(6.8*(rsigma-1.2)),0.6);
620 gcor01[i] = TMath::Max(0.72 +TMath::Exp(2.04*(rsigma-1.2)),0.6);
621 gcorp[i] = TMath::Max(TMath::Power((rsigma+0.5),1.5),1.2);
625 for (Int_t i=3;i<10000;i++){
629 Float_t amp = float(i);
630 Float_t padlength =0.75;
631 gnoise1 = 0.0004/padlength;
632 Float_t nel = 0.268*amp;
633 Float_t nprim = 0.155*amp;
634 ggg1[i] = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.001*nel/(padlength*padlength))/nel;
635 glandau1[i] = (2.+0.12*nprim)*0.5* (2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
636 if (glandau1[i]>1) glandau1[i]=1;
637 glandau1[i]*=padlength*padlength/12.;
641 gnoise2 = 0.0004/padlength;
644 ggg2[i] = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
645 glandau2[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
646 if (glandau2[i]>1) glandau2[i]=1;
647 glandau2[i]*=padlength*padlength/12.;
652 gnoise3 = 0.0004/padlength;
655 ggg3[i] = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
656 glandau3[i] = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
657 if (glandau3[i]>1) glandau3[i]=1;
658 glandau3[i]*=padlength*padlength/12.;
666 Int_t amp = int(TMath::Abs(cl->GetQ()));
668 seed->SetErrorY2(1.);
672 Float_t z = TMath::Abs(fParam->GetZLength()-TMath::Abs(seed->GetZ()));
673 Int_t ctype = cl->GetType();
674 Float_t padlength= GetPadPitchLength(seed->fRow);
676 Float_t angle2 = seed->GetSnp()*seed->GetSnp();
677 // if (angle2<0.6) angle2 = 0.6;
678 angle2 = seed->GetTgl()*seed->GetTgl()*(1+angle2/(1-angle2));
681 Int_t rsigmaz = int(100.*cl->GetSigmaZ2()/(seed->fCurrentSigmaZ2));
684 if (fSectors==fInnerSec){
686 res = ggg1[amp]*z+glandau1[amp]*angle2;
687 if (ctype==0) res *= gcor01[rsigmaz];
690 res*= gcorp[rsigmaz];
696 res = ggg2[amp]*z+glandau2[amp]*angle2;
697 if (ctype==0) res *= gcor02[rsigmaz];
700 res*= gcorp[rsigmaz];
705 res = ggg3[amp]*z+glandau3[amp]*angle2;
706 if (ctype==0) res *= gcor02[rsigmaz];
709 res*= gcorp[rsigmaz];
718 if ((ctype<0) &&<70){
726 seed->SetErrorZ2(res);
733 Double_t AliTPCtrackerMI::ErrZ2(AliTPCseed* seed, AliTPCclusterMI * cl){
736 //seed->SetErrorZ2(0.1);
740 Float_t z = TMath::Abs(fParam->GetZLength()-TMath::Abs(seed->GetZ()));
742 Float_t rsigmaz = cl->GetSigmaZ2()/(seed->fCurrentSigmaZ2);
743 Int_t ctype = cl->GetType();
744 Float_t amp = TMath::Abs(cl->GetQ());
749 Float_t landau=2 ; //landau fluctuation part
750 Float_t gg=2; // gg fluctuation part
751 Float_t padlength= GetPadPitchLength(seed->GetX());
753 if (fSectors==fInnerSec){
754 snoise2 = 0.0004/padlength;
757 gg = (2+0.001*nel/(padlength*padlength))/nel;
758 landau = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
759 if (landau>1) landau=1;
762 snoise2 = 0.0004/padlength;
765 gg = (2+0.0008*nel/(padlength*padlength))/nel;
766 landau = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
767 if (landau>1) landau=1;
769 Float_t sdiff = gg*fParam->GetDiffT()*fParam->GetDiffT()*z;
772 Float_t angle2 = seed->GetSnp()*seed->GetSnp();
773 angle2 = TMath::Sqrt((1-angle2));
774 if (angle2<0.6) angle2 = 0.6;
777 Float_t angle = seed->GetTgl()/angle2;
778 Float_t angular = landau*angle*angle*padlength*padlength/12.;
779 Float_t res = sdiff + angular;
782 if ((ctype==0) && (fSectors ==fOuterSec))
783 res *= 0.81 +TMath::Exp(6.8*(rsigmaz-1.2));
785 if ((ctype==0) && (fSectors ==fInnerSec))
786 res *= 0.72 +TMath::Exp(2.04*(rsigmaz-1.2));
790 res*= TMath::Power(rsigmaz+0.5,1.5); //0.31+0.147*ctype;
796 if ((ctype<0) &&<70){
804 seed->SetErrorZ2(res);
811 void AliTPCseed::Reset(Bool_t all)
815 SetNumberOfClusters(0);
821 for (Int_t i=0;i<8;i++){
822 delete [] fTrackPoints[i];
830 for (Int_t i=0;i<200;i++) SetClusterIndex2(i,-3);
831 for (Int_t i=0;i<160;i++) fClusterPointer[i]=0;
837 void AliTPCseed::Modify(Double_t factor)
840 //------------------------------------------------------------------
841 //This function makes a track forget its history :)
842 //------------------------------------------------------------------
848 fC10*=0; fC11*=factor;
849 fC20*=0; fC21*=0; fC22*=factor;
850 fC30*=0; fC31*=0; fC32*=0; fC33*=factor;
851 fC40*=0; fC41*=0; fC42*=0; fC43*=0; fC44*=factor;
852 SetNumberOfClusters(0);
856 fCurrentSigmaY2 = 0.000005;
857 fCurrentSigmaZ2 = 0.000005;
866 Int_t AliTPCseed::GetProlongation(Double_t xk, Double_t &y, Double_t & z) const
868 //-----------------------------------------------------------------
869 // This function find proloncation of a track to a reference plane x=xk.
870 // doesn't change internal state of the track
871 //-----------------------------------------------------------------
873 Double_t x1=fX, x2=x1+(xk-x1), dx=x2-x1;
875 if (TMath::Abs(fP4*xk - fP2) >= 0.999) {
879 // Double_t y1=fP0, z1=fP1;
880 Double_t c1=fP4*x1 - fP2, r1=sqrt(1.- c1*c1);
881 Double_t c2=fP4*x2 - fP2, r2=sqrt(1.- c2*c2);
885 //y += dx*(c1+c2)/(r1+r2);
886 //z += dx*(c1+c2)/(c1*r2 + c2*r1)*fP3;
888 Double_t dy = dx*(c1+c2)/(r1+r2);
891 Double_t delta = fP4*dx*(c1+c2)/(c1*r2 + c2*r1);
893 if (TMath::Abs(delta)>0.0001){
894 dz = fP3*TMath::ASin(delta)/fP4;
896 dz = dx*fP3*(c1+c2)/(c1*r2 + c2*r1);
899 dz = fP3*AliTPCFastMath::FastAsin(delta)/fP4;
909 //_____________________________________________________________________________
910 Double_t AliTPCseed::GetPredictedChi2(const AliTPCclusterMI *c) const
912 //-----------------------------------------------------------------
913 // This function calculates a predicted chi2 increment.
914 //-----------------------------------------------------------------
915 //Double_t r00=c->GetSigmaY2(), r01=0., r11=c->GetSigmaZ2();
916 Double_t r00=fErrorY2, r01=0., r11=fErrorZ2;
917 r00+=fC00; r01+=fC10; r11+=fC11;
919 Double_t det=r00*r11 - r01*r01;
920 if (TMath::Abs(det) < 1.e-10) {
921 Int_t n=GetNumberOfClusters();
922 if (n>4) cerr<<n<<" AliKalmanTrack warning: Singular matrix !\n";
925 Double_t tmp=r00; r00=r11; r11=tmp; r01=-r01;
927 Double_t dy=c->GetY() - fP0, dz=c->GetZ() - fP1;
929 return (dy*r00*dy + 2*r01*dy*dz + dz*r11*dz)/det;
933 //_________________________________________________________________________________________
936 Int_t AliTPCseed::Compare(const TObject *o) const {
937 //-----------------------------------------------------------------
938 // This function compares tracks according to the sector - for given sector according z
939 //-----------------------------------------------------------------
940 AliTPCseed *t=(AliTPCseed*)o;
943 if (t->fRelativeSector>fRelativeSector) return -1;
944 if (t->fRelativeSector<fRelativeSector) return 1;
945 Double_t z2 = t->GetZ();
946 Double_t z1 = GetZ();
948 if (z2<z1) return -1;
953 f2 = 1-20*TMath::Sqrt(t->fC44)/(TMath::Abs(t->GetC())+0.0066);
954 if (t->fBConstrain) f2=1.2;
957 f1 = 1-20*TMath::Sqrt(fC44)/(TMath::Abs(GetC())+0.0066);
959 if (fBConstrain) f1=1.2;
961 if (t->GetNumberOfClusters()*f2 <GetNumberOfClusters()*f1) return -1;
966 void AliTPCtrackerMI::RotateToLocal(AliTPCseed *seed)
968 //rotate to track "local coordinata
969 Float_t x = seed->GetX();
970 Float_t y = seed->GetY();
971 Float_t ymax = x*TMath::Tan(0.5*fSectors->GetAlpha());
974 seed->fRelativeSector= (seed->fRelativeSector+1) % fN;
975 if (!seed->Rotate(fSectors->GetAlpha()))
977 } else if (y <-ymax) {
978 seed->fRelativeSector= (seed->fRelativeSector-1+fN) % fN;
979 if (!seed->Rotate(-fSectors->GetAlpha()))
988 //_____________________________________________________________________________
989 Int_t AliTPCseed::Update(const AliTPCclusterMI *c, Double_t chisq, UInt_t /*index*/) {
990 //-----------------------------------------------------------------
991 // This function associates a cluster with this track.
992 //-----------------------------------------------------------------
993 Double_t r00=fErrorY2, r01=0., r11=fErrorZ2;
995 r00+=fC00; r01+=fC10; r11+=fC11;
996 Double_t det=r00*r11 - r01*r01;
997 Double_t tmp=r00; r00=r11/det; r11=tmp/det; r01=-r01/det;
999 Double_t k00=fC00*r00+fC10*r01, k01=fC00*r01+fC10*r11;
1000 Double_t k10=fC10*r00+fC11*r01, k11=fC10*r01+fC11*r11;
1001 Double_t k20=fC20*r00+fC21*r01, k21=fC20*r01+fC21*r11;
1002 Double_t k30=fC30*r00+fC31*r01, k31=fC30*r01+fC31*r11;
1003 Double_t k40=fC40*r00+fC41*r01, k41=fC40*r01+fC41*r11;
1005 Double_t dy=c->GetY() - fP0, dz=c->GetZ() - fP1;
1006 Double_t cur=fP4 + k40*dy + k41*dz, eta=fP2 + k20*dy + k21*dz;
1007 if (TMath::Abs(cur*fX-eta) >= 0.9) {
1011 fP0 += k00*dy + k01*dz;
1012 fP1 += k10*dy + k11*dz;
1014 fP3 += k30*dy + k31*dz;
1017 Double_t c01=fC10, c02=fC20, c03=fC30, c04=fC40;
1018 Double_t c12=fC21, c13=fC31, c14=fC41;
1020 fC00-=k00*fC00+k01*fC10; fC10-=k00*c01+k01*fC11;
1021 fC20-=k00*c02+k01*c12; fC30-=k00*c03+k01*c13;
1022 fC40-=k00*c04+k01*c14;
1024 fC11-=k10*c01+k11*fC11;
1025 fC21-=k10*c02+k11*c12; fC31-=k10*c03+k11*c13;
1026 fC41-=k10*c04+k11*c14;
1028 fC22-=k20*c02+k21*c12; fC32-=k20*c03+k21*c13;
1029 fC42-=k20*c04+k21*c14;
1031 fC33-=k30*c03+k31*c13;
1032 fC43-=k40*c03+k41*c13;
1034 fC44-=k40*c04+k41*c14;
1036 Int_t n=GetNumberOfClusters();
1038 SetNumberOfClusters(n+1);
1039 SetChi2(GetChi2()+chisq);
1046 //_____________________________________________________________________________
1047 Double_t AliTPCtrackerMI::F1old(Double_t x1,Double_t y1,
1048 Double_t x2,Double_t y2,
1049 Double_t x3,Double_t y3)
1051 //-----------------------------------------------------------------
1052 // Initial approximation of the track curvature
1053 //-----------------------------------------------------------------
1054 Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
1055 Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
1056 (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
1057 Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
1058 (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
1060 Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
1061 if ( xr*xr+yr*yr<=0.00000000000001) return 100;
1062 return -xr*yr/sqrt(xr*xr+yr*yr);
1067 //_____________________________________________________________________________
1068 Double_t AliTPCtrackerMI::F1(Double_t x1,Double_t y1,
1069 Double_t x2,Double_t y2,
1070 Double_t x3,Double_t y3)
1072 //-----------------------------------------------------------------
1073 // Initial approximation of the track curvature
1074 //-----------------------------------------------------------------
1080 Double_t det = x3*y2-x2*y3;
1085 Double_t u = 0.5* (x2*(x2-x3)+y2*(y2-y3))/det;
1086 Double_t x0 = x3*0.5-y3*u;
1087 Double_t y0 = y3*0.5+x3*u;
1088 Double_t c2 = 1/TMath::Sqrt(x0*x0+y0*y0);
1094 Double_t AliTPCtrackerMI::F2(Double_t x1,Double_t y1,
1095 Double_t x2,Double_t y2,
1096 Double_t x3,Double_t y3)
1098 //-----------------------------------------------------------------
1099 // Initial approximation of the track curvature
1100 //-----------------------------------------------------------------
1106 Double_t det = x3*y2-x2*y3;
1111 Double_t u = 0.5* (x2*(x2-x3)+y2*(y2-y3))/det;
1112 Double_t x0 = x3*0.5-y3*u;
1113 Double_t y0 = y3*0.5+x3*u;
1114 Double_t c2 = 1/TMath::Sqrt(x0*x0+y0*y0);
1123 //_____________________________________________________________________________
1124 Double_t AliTPCtrackerMI::F2old(Double_t x1,Double_t y1,
1125 Double_t x2,Double_t y2,
1126 Double_t x3,Double_t y3)
1128 //-----------------------------------------------------------------
1129 // Initial approximation of the track curvature times center of curvature
1130 //-----------------------------------------------------------------
1131 Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
1132 Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
1133 (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
1134 Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
1135 (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
1137 Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
1139 return -a/(d*y1-b)*xr/sqrt(xr*xr+yr*yr);
1142 //_____________________________________________________________________________
1143 Double_t AliTPCtrackerMI::F3(Double_t x1,Double_t y1,
1144 Double_t x2,Double_t y2,
1145 Double_t z1,Double_t z2)
1147 //-----------------------------------------------------------------
1148 // Initial approximation of the tangent of the track dip angle
1149 //-----------------------------------------------------------------
1150 return (z1 - z2)/sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
1154 Double_t AliTPCtrackerMI::F3n(Double_t x1,Double_t y1,
1155 Double_t x2,Double_t y2,
1156 Double_t z1,Double_t z2, Double_t c)
1158 //-----------------------------------------------------------------
1159 // Initial approximation of the tangent of the track dip angle
1160 //-----------------------------------------------------------------
1164 //angle1 = (z1-z2)*c/(TMath::ASin(c*x1-ni)-TMath::ASin(c*x2-ni));
1166 Double_t d = TMath::Sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
1167 if (TMath::Abs(d*c*0.5)>1) return 0;
1168 // Double_t angle2 = TMath::ASin(d*c*0.5);
1169 // Double_t angle2 = AliTPCFastMath::FastAsin(d*c*0.5);
1170 Double_t angle2 = (d*c*0.5>0.1)? TMath::ASin(d*c*0.5): AliTPCFastMath::FastAsin(d*c*0.5);
1172 angle2 = (z1-z2)*c/(angle2*2.);
1176 Bool_t AliTPCtrackerMI::GetProlongation(Double_t x1, Double_t x2, Double_t x[5], Double_t &y, Double_t &z)
1177 {//-----------------------------------------------------------------
1178 // This function find proloncation of a track to a reference plane x=x2.
1179 //-----------------------------------------------------------------
1183 if (TMath::Abs(x[4]*x1 - x[2]) >= 0.999) {
1187 Double_t c1=x[4]*x1 - x[2], r1=sqrt(1.- c1*c1);
1188 Double_t c2=x[4]*x2 - x[2], r2=sqrt(1.- c2*c2);
1192 Double_t dy = dx*(c1+c2)/(r1+r2);
1195 Double_t delta = x[4]*dx*(c1+c2)/(c1*r2 + c2*r1);
1197 if (TMath::Abs(delta)>0.01){
1198 dz = x[3]*TMath::ASin(delta)/x[4];
1200 dz = x[3]*AliTPCFastMath::FastAsin(delta)/x[4];
1203 //dz = x[3]*AliTPCFastMath::FastAsin(delta)/x[4];
1211 Int_t AliTPCtrackerMI::LoadClusters (TTree *tree)
1216 return LoadClusters();
1219 Int_t AliTPCtrackerMI::LoadClusters()
1222 // load clusters to the memory
1223 AliTPCClustersRow *clrow= new AliTPCClustersRow;
1224 clrow->SetClass("AliTPCclusterMI");
1226 clrow->GetArray()->ExpandCreateFast(10000);
1228 // TTree * tree = fClustersArray.GetTree();
1230 TTree * tree = fInput;
1231 TBranch * br = tree->GetBranch("Segment");
1232 br->SetAddress(&clrow);
1234 Int_t j=Int_t(tree->GetEntries());
1235 for (Int_t i=0; i<j; i++) {
1239 fParam->AdjustSectorRow(clrow->GetID(),sec,row);
1241 AliTPCRow * tpcrow=0;
1244 tpcrow = &(fInnerSec[sec%fkNIS][row]);
1248 tpcrow = &(fOuterSec[(sec-fkNIS*2)%fkNOS][row]);
1249 left = (sec-fkNIS*2)/fkNOS;
1252 tpcrow->fN1 = clrow->GetArray()->GetEntriesFast();
1253 tpcrow->fClusters1 = new AliTPCclusterMI[tpcrow->fN1];
1254 for (Int_t i=0;i<tpcrow->fN1;i++)
1255 tpcrow->fClusters1[i] = *(AliTPCclusterMI*)(clrow->GetArray()->At(i));
1258 tpcrow->fN2 = clrow->GetArray()->GetEntriesFast();
1259 tpcrow->fClusters2 = new AliTPCclusterMI[tpcrow->fN2];
1260 for (Int_t i=0;i<tpcrow->fN2;i++)
1261 tpcrow->fClusters2[i] = *(AliTPCclusterMI*)(clrow->GetArray()->At(i));
1272 void AliTPCtrackerMI::UnloadClusters()
1275 // unload clusters from the memory
1277 Int_t nrows = fOuterSec->GetNRows();
1278 for (Int_t sec = 0;sec<fkNOS;sec++)
1279 for (Int_t row = 0;row<nrows;row++){
1280 AliTPCRow* tpcrow = &(fOuterSec[sec%fkNOS][row]);
1282 // if (tpcrow->fClusters1) delete []tpcrow->fClusters1;
1283 // if (tpcrow->fClusters2) delete []tpcrow->fClusters2;
1285 tpcrow->ResetClusters();
1288 nrows = fInnerSec->GetNRows();
1289 for (Int_t sec = 0;sec<fkNIS;sec++)
1290 for (Int_t row = 0;row<nrows;row++){
1291 AliTPCRow* tpcrow = &(fInnerSec[sec%fkNIS][row]);
1293 // if (tpcrow->fClusters1) delete []tpcrow->fClusters1;
1294 //if (tpcrow->fClusters2) delete []tpcrow->fClusters2;
1296 tpcrow->ResetClusters();
1303 //_____________________________________________________________________________
1304 Int_t AliTPCtrackerMI::LoadOuterSectors() {
1305 //-----------------------------------------------------------------
1306 // This function fills outer TPC sectors with clusters.
1307 //-----------------------------------------------------------------
1308 Int_t nrows = fOuterSec->GetNRows();
1310 for (Int_t sec = 0;sec<fkNOS;sec++)
1311 for (Int_t row = 0;row<nrows;row++){
1312 AliTPCRow* tpcrow = &(fOuterSec[sec%fkNOS][row]);
1313 Int_t sec2 = sec+2*fkNIS;
1315 Int_t ncl = tpcrow->fN1;
1317 AliTPCclusterMI *c= &(tpcrow->fClusters1[ncl]);
1318 index=(((sec2<<8)+row)<<16)+ncl;
1319 tpcrow->InsertCluster(c,index);
1324 AliTPCclusterMI *c= &(tpcrow->fClusters2[ncl]);
1325 index=((((sec2+fkNOS)<<8)+row)<<16)+ncl;
1326 tpcrow->InsertCluster(c,index);
1329 // write indexes for fast acces
1331 for (Int_t i=0;i<510;i++)
1332 tpcrow->fFastCluster[i]=-1;
1333 for (Int_t i=0;i<tpcrow->GetN();i++){
1334 Int_t zi = Int_t((*tpcrow)[i]->GetZ()+255.);
1335 tpcrow->fFastCluster[zi]=i; // write index
1338 for (Int_t i=0;i<510;i++){
1339 if (tpcrow->fFastCluster[i]<0)
1340 tpcrow->fFastCluster[i] = last;
1342 last = tpcrow->fFastCluster[i];
1351 //_____________________________________________________________________________
1352 Int_t AliTPCtrackerMI::LoadInnerSectors() {
1353 //-----------------------------------------------------------------
1354 // This function fills inner TPC sectors with clusters.
1355 //-----------------------------------------------------------------
1356 Int_t nrows = fInnerSec->GetNRows();
1358 for (Int_t sec = 0;sec<fkNIS;sec++)
1359 for (Int_t row = 0;row<nrows;row++){
1360 AliTPCRow* tpcrow = &(fInnerSec[sec%fkNIS][row]);
1363 Int_t ncl = tpcrow->fN1;
1365 AliTPCclusterMI *c= &(tpcrow->fClusters1[ncl]);
1366 index=(((sec<<8)+row)<<16)+ncl;
1367 tpcrow->InsertCluster(c,index);
1372 AliTPCclusterMI *c= &(tpcrow->fClusters2[ncl]);
1373 index=((((sec+fkNIS)<<8)+row)<<16)+ncl;
1374 tpcrow->InsertCluster(c,index);
1377 // write indexes for fast acces
1379 for (Int_t i=0;i<510;i++)
1380 tpcrow->fFastCluster[i]=-1;
1381 for (Int_t i=0;i<tpcrow->GetN();i++){
1382 Int_t zi = Int_t((*tpcrow)[i]->GetZ()+255.);
1383 tpcrow->fFastCluster[zi]=i; // write index
1386 for (Int_t i=0;i<510;i++){
1387 if (tpcrow->fFastCluster[i]<0)
1388 tpcrow->fFastCluster[i] = last;
1390 last = tpcrow->fFastCluster[i];
1402 //_________________________________________________________________________
1403 AliTPCclusterMI *AliTPCtrackerMI::GetClusterMI(Int_t index) const {
1404 //--------------------------------------------------------------------
1405 // Return pointer to a given cluster
1406 //--------------------------------------------------------------------
1407 Int_t sec=(index&0xff000000)>>24;
1408 Int_t row=(index&0x00ff0000)>>16;
1409 Int_t ncl=(index&0x00007fff)>>00;
1411 const AliTPCRow * tpcrow=0;
1412 AliTPCclusterMI * clrow =0;
1414 tpcrow = &(fInnerSec[sec%fkNIS][row]);
1416 clrow = tpcrow->fClusters1;
1418 clrow = tpcrow->fClusters2;
1421 tpcrow = &(fOuterSec[(sec-fkNIS*2)%fkNOS][row]);
1422 if (sec-2*fkNIS<fkNOS)
1423 clrow = tpcrow->fClusters1;
1425 clrow = tpcrow->fClusters2;
1427 if (tpcrow==0) return 0;
1428 if (tpcrow->GetN()<=ncl) return 0;
1429 // return (AliTPCclusterMI*)(*tpcrow)[ncl];
1430 return &(clrow[ncl]);
1436 Int_t AliTPCtrackerMI::FollowToNext(AliTPCseed& t, Int_t nr) {
1437 //-----------------------------------------------------------------
1438 // This function tries to find a track prolongation to next pad row
1439 //-----------------------------------------------------------------
1441 Double_t x= GetXrow(nr), ymax=GetMaxY(nr);
1442 AliTPCclusterMI *cl=0;
1443 Int_t tpcindex= t.GetClusterIndex2(nr);
1445 // update current shape info every 5 pad-row
1446 // if ( (nr%5==0) || t.GetNumberOfClusters()<2 || (t.fCurrentSigmaY2<0.0001) ){
1450 if (fIteration>0 && tpcindex>=-1){ //if we have already clusters
1452 if (tpcindex==-1) return 0; //track in dead zone
1454 cl = t.fClusterPointer[nr];
1455 if ( (cl==0) ) cl = GetClusterMI(tpcindex);
1456 t.fCurrentClusterIndex1 = tpcindex;
1459 Int_t relativesector = ((tpcindex&0xff000000)>>24)%18; // if previously accepted cluster in different sector
1460 Float_t angle = relativesector*fSectors->GetAlpha()+fSectors->GetAlphaShift();
1462 if (angle<-TMath::Pi()) angle += 2*TMath::Pi();
1463 if (angle>=TMath::Pi()) angle -= 2*TMath::Pi();
1465 if (TMath::Abs(angle-t.GetAlpha())>0.001){
1466 Double_t rotation = angle-t.GetAlpha();
1467 t.fRelativeSector= relativesector;
1472 t.fCurrentCluster = cl;
1474 Int_t accept = AcceptCluster(&t,t.fCurrentCluster,1.);
1475 if ((tpcindex&0x8000)==0) accept =0;
1477 //if founded cluster is acceptible
1478 if (cl->IsUsed(11)) { // id cluster is shared inrease uncertainty
1485 UpdateTrack(&t,accept);
1490 if (fIteration>1) return 0; // not look for new cluster during refitting
1493 if (TMath::Abs(t.GetSnp())>0.95 || TMath::Abs(x*t.GetC()-t.GetEta())>0.95) return 0;
1494 Double_t y=t.GetYat(x);
1495 if (TMath::Abs(y)>ymax){
1497 t.fRelativeSector= (t.fRelativeSector+1) % fN;
1498 if (!t.Rotate(fSectors->GetAlpha()))
1500 } else if (y <-ymax) {
1501 t.fRelativeSector= (t.fRelativeSector-1+fN) % fN;
1502 if (!t.Rotate(-fSectors->GetAlpha()))
1508 if (!t.PropagateTo(x)) {
1509 if (fIteration==0) t.fRemoval = 10;
1513 Double_t z=t.GetZ();
1515 const AliTPCRow &krow=GetRow(t.fRelativeSector,nr);
1516 if ( (t.GetSigmaY2()<0) || t.GetSigmaZ2()<0) return 0;
1518 Double_t roadz = 1.;
1520 if (TMath::Abs(TMath::Abs(y)-ymax)<krow.fDeadZone){
1522 t.SetClusterIndex2(nr,-1);
1527 if (TMath::Abs(z)<(1.05*x+10)) t.fNFoundable++;
1533 // cl = krow.FindNearest2(y+10.,z,roady,roadz,index);
1534 cl = krow.FindNearest2(y,z,roady,roadz,index);
1535 if (cl) t.fCurrentClusterIndex1 = krow.GetIndex(index);
1538 t.fCurrentCluster = cl;
1540 if (fIteration==2&&cl->IsUsed(10)) return 0;
1541 Int_t accept = AcceptCluster(&t,t.fCurrentCluster,1.);
1542 if (fIteration==2&&cl->IsUsed(11)) {
1549 if (t.fCurrentCluster->IsUsed(10)){
1554 if (t.fNShared>0.7*t.GetNumberOfClusters()) {
1560 if (accept<3) UpdateTrack(&t,accept);
1563 if ( fIteration==0 && t.fNFoundable*0.5 > t.GetNumberOfClusters()) t.fRemoval=10;
1569 Int_t AliTPCtrackerMI::FollowToNextFast(AliTPCseed& t, Int_t nr) {
1570 //-----------------------------------------------------------------
1571 // This function tries to find a track prolongation to next pad row
1572 //-----------------------------------------------------------------
1574 Double_t x= GetXrow(nr), ymax=GetMaxY(nr);
1576 if (!t.GetProlongation(x,y,z)) {
1582 if (TMath::Abs(y)>ymax){
1585 t.fRelativeSector= (t.fRelativeSector+1) % fN;
1586 if (!t.Rotate(fSectors->GetAlpha()))
1588 } else if (y <-ymax) {
1589 t.fRelativeSector= (t.fRelativeSector-1+fN) % fN;
1590 if (!t.Rotate(-fSectors->GetAlpha()))
1593 if (!t.PropagateTo(x)) {
1596 t.GetProlongation(x,y,z);
1599 // update current shape info every 3 pad-row
1600 if ( (nr%6==0) || t.GetNumberOfClusters()<2 || (t.fCurrentSigmaY2<0.0001) ){
1601 // t.fCurrentSigmaY = GetSigmaY(&t);
1602 //t.fCurrentSigmaZ = GetSigmaZ(&t);
1606 AliTPCclusterMI *cl=0;
1611 const AliTPCRow &krow=GetRow(t.fRelativeSector,nr);
1612 if ( (t.GetSigmaY2()<0) || t.GetSigmaZ2()<0) return 0;
1614 Double_t roadz = 1.;
1617 if (TMath::Abs(TMath::Abs(y)-ymax)<krow.fDeadZone){
1619 t.SetClusterIndex2(row,-1);
1624 if (TMath::Abs(z)>(1.05*x+10)) t.SetClusterIndex2(row,-1);
1628 if ((cl==0)&&(krow)) {
1629 // cl = krow.FindNearest2(y+10,z,roady,roadz,index);
1630 cl = krow.FindNearest2(y,z,roady,roadz,index);
1632 if (cl) t.fCurrentClusterIndex1 = krow.GetIndex(index);
1636 t.fCurrentCluster = cl;
1637 // Int_t accept = AcceptCluster(&t,t.fCurrentCluster,1.);
1639 t.SetClusterIndex2(row,index);
1640 t.fClusterPointer[row] = cl;
1648 Int_t AliTPCtrackerMI::UpdateClusters(AliTPCseed& t, Int_t nr) {
1649 //-----------------------------------------------------------------
1650 // This function tries to find a track prolongation to next pad row
1651 //-----------------------------------------------------------------
1652 t.fCurrentCluster = 0;
1653 t.fCurrentClusterIndex1 = 0;
1655 Double_t xt=t.GetX();
1656 Int_t row = GetRowNumber(xt)-1;
1657 Double_t ymax= GetMaxY(nr);
1659 if (row < nr) return 1; // don't prolongate if not information until now -
1660 if (TMath::Abs(t.GetSnp())>0.9 && t.GetNumberOfClusters()>40. && fIteration!=2) {
1662 return 0; // not prolongate strongly inclined tracks
1664 if (TMath::Abs(t.GetSnp())>0.95) {
1666 return 0; // not prolongate strongly inclined tracks
1669 Double_t x= GetXrow(nr);
1671 //t.PropagateTo(x+0.02);
1672 //t.PropagateTo(x+0.01);
1673 if (!t.PropagateTo(x)){
1680 if (TMath::Abs(y)>ymax){
1682 t.fRelativeSector= (t.fRelativeSector+1) % fN;
1683 if (!t.Rotate(fSectors->GetAlpha()))
1685 } else if (y <-ymax) {
1686 t.fRelativeSector= (t.fRelativeSector-1+fN) % fN;
1687 if (!t.Rotate(-fSectors->GetAlpha()))
1690 // if (!t.PropagateTo(x)){
1698 AliTPCRow &krow=GetRow(t.fRelativeSector,nr);
1700 if (TMath::Abs(TMath::Abs(y)-ymax)<krow.fDeadZone){
1702 t.SetClusterIndex2(nr,-1);
1707 if (TMath::Abs(t.GetZ())<(1.05*t.GetX()+10)) t.fNFoundable++;
1713 if ( (nr%6==0) || t.GetNumberOfClusters()<2){
1714 // t.fCurrentSigmaY = GetSigmaY(&t);
1715 //t.fCurrentSigmaZ = GetSigmaZ(&t);
1719 AliTPCclusterMI *cl=0;
1722 Double_t roady = 1.;
1723 Double_t roadz = 1.;
1727 index = t.GetClusterIndex2(nr);
1728 if ( (index>0) && (index&0x8000)==0){
1729 cl = t.fClusterPointer[nr];
1730 if ( (cl==0) && (index>0)) cl = GetClusterMI(index);
1731 t.fCurrentClusterIndex1 = index;
1733 t.fCurrentCluster = cl;
1740 //cl = krow.FindNearest2(y+10,z,roady,roadz,index);
1741 cl = krow.FindNearest2(y,z,roady,roadz,index);
1744 if (cl) t.fCurrentClusterIndex1 = krow.GetIndex(index);
1745 t.fCurrentCluster = cl;
1751 Int_t AliTPCtrackerMI::FollowToNextCluster(AliTPCseed & t, Int_t nr) {
1752 //-----------------------------------------------------------------
1753 // This function tries to find a track prolongation to next pad row
1754 //-----------------------------------------------------------------
1756 //update error according neighborhoud
1758 if (t.fCurrentCluster) {
1760 Int_t accept = AcceptCluster(&t,t.fCurrentCluster,1.);
1762 if (t.fCurrentCluster->IsUsed(10)){
1768 if (t.fNShared>0.7*t.GetNumberOfClusters()) {
1773 if (fIteration>0) accept = 0;
1774 if (accept<3) UpdateTrack(&t,accept);
1778 if ( ( (t.GetSigmaY2()+t.GetSigmaZ2())>0.16)&& t.GetNumberOfClusters()>18) t.fRemoval=10;
1779 if ( t.GetChi2()/t.GetNumberOfClusters()>6 &&t.GetNumberOfClusters()>18) t.fRemoval=10;
1781 if (( (t.fNFoundable*0.5 > t.GetNumberOfClusters()) || t.fNoCluster>15)) t.fRemoval=10;
1789 //_____________________________________________________________________________
1790 Int_t AliTPCtrackerMI::FollowProlongation(AliTPCseed& t, Int_t rf, Int_t step) {
1791 //-----------------------------------------------------------------
1792 // This function tries to find a track prolongation.
1793 //-----------------------------------------------------------------
1794 Double_t xt=t.GetX();
1796 Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
1797 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
1798 if (alpha < 0. ) alpha += 2.*TMath::Pi();
1800 t.fRelativeSector = Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN;
1802 Int_t first = GetRowNumber(xt)-1;
1803 for (Int_t nr= first; nr>=rf; nr-=step) {
1804 if (nr<fInnerSec->GetNRows())
1805 fSectors = fInnerSec;
1807 fSectors = fOuterSec;
1808 if (FollowToNext(t,nr)==0)
1817 //_____________________________________________________________________________
1818 Int_t AliTPCtrackerMI::FollowProlongationFast(AliTPCseed& t, Int_t rf, Int_t step) {
1819 //-----------------------------------------------------------------
1820 // This function tries to find a track prolongation.
1821 //-----------------------------------------------------------------
1822 Double_t xt=t.GetX();
1824 Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
1825 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
1826 if (alpha < 0. ) alpha += 2.*TMath::Pi();
1827 t.fRelativeSector = Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN;
1829 for (Int_t nr=GetRowNumber(xt)-1; nr>=rf; nr-=step) {
1831 if (FollowToNextFast(t,nr)==0)
1832 if (!t.IsActive()) return 0;
1842 Int_t AliTPCtrackerMI::FollowBackProlongation(AliTPCseed& t, Int_t rf) {
1843 //-----------------------------------------------------------------
1844 // This function tries to find a track prolongation.
1845 //-----------------------------------------------------------------
1846 // Double_t xt=t.GetX();
1848 Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
1849 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
1850 if (alpha < 0. ) alpha += 2.*TMath::Pi();
1851 t.fRelativeSector = Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN;
1853 Int_t first = t.fFirstPoint;
1855 if (first<0) first=0;
1856 for (Int_t nr=first; nr<=rf; nr++) {
1857 //if ( (t.GetSnp()<0.9))
1858 if (nr<fInnerSec->GetNRows())
1859 fSectors = fInnerSec;
1861 fSectors = fOuterSec;
1871 Float_t AliTPCtrackerMI::OverlapFactor(AliTPCseed * s1, AliTPCseed * s2, Int_t &sum1, Int_t & sum2)
1879 Float_t dz2 =(s1->GetZ() - s2->GetZ());
1882 Float_t dy2 =TMath::Abs((s1->GetY() - s2->GetY()));
1884 Float_t distance = TMath::Sqrt(dz2+dy2);
1885 if (distance>4.) return 0; // if there are far away - not overlap - to reduce combinatorics
1888 Int_t firstpoint = TMath::Min(s1->fFirstPoint,s2->fFirstPoint);
1889 Int_t lastpoint = TMath::Max(s1->fLastPoint,s2->fLastPoint);
1894 if (firstpoint>lastpoint) {
1895 firstpoint =lastpoint;
1900 for (Int_t i=firstpoint-1;i<lastpoint+1;i++){
1901 if (s1->GetClusterIndex2(i)>0) sum1++;
1902 if (s2->GetClusterIndex2(i)>0) sum2++;
1903 if (s1->GetClusterIndex2(i)==s2->GetClusterIndex2(i) && s1->GetClusterIndex2(i)>0) {
1907 if (sum<5) return 0;
1909 Float_t summin = TMath::Min(sum1+1,sum2+1);
1910 Float_t ratio = (sum+1)/Float_t(summin);
1914 void AliTPCtrackerMI::SignShared(AliTPCseed * s1, AliTPCseed * s2)
1918 if (TMath::Abs(s1->GetC()-s2->GetC())>0.004) return;
1919 if (TMath::Abs(s1->GetTgl()-s2->GetTgl())>0.6) return;
1921 Float_t dz2 =(s1->GetZ() - s2->GetZ());
1923 Float_t dy2 =(s1->GetY() - s2->GetY());
1925 Float_t distance = dz2+dy2;
1926 if (distance>325.) return ; // if there are far away - not overlap - to reduce combinatorics
1931 Int_t firstpoint = TMath::Max(s1->fFirstPoint,s2->fFirstPoint);
1932 Int_t lastpoint = TMath::Min(s1->fLastPoint,s2->fLastPoint);
1934 if (firstpoint>=lastpoint-5) return;;
1936 for (Int_t i=firstpoint;i<lastpoint;i++){
1937 // if ( (s1->GetClusterIndex2(i)&0xFFFF8FFF)==(s2->GetClusterIndex2(i)&0xFFFF8FFF) && s1->GetClusterIndex2(i)>0) {
1938 if ( (s1->GetClusterIndex2(i))==(s2->GetClusterIndex2(i)) && s1->GetClusterIndex2(i)>0) {
1945 for (Int_t i=firstpoint;i<lastpoint;i++){
1946 // if ( (s1->GetClusterIndex2(i)&0xFFFF8FFF)==(s2->GetClusterIndex2(i)&0xFFFF8FFF) && s1->GetClusterIndex2(i)>0) {
1947 if ( (s1->GetClusterIndex2(i))==(s2->GetClusterIndex2(i)) && s1->GetClusterIndex2(i)>0) {
1948 AliTPCTrackerPoint *p1 = s1->GetTrackPoint(i);
1949 AliTPCTrackerPoint *p2 = s2->GetTrackPoint(i);;
1950 if (s1->IsActive()&&s2->IsActive()){
1951 p1->fIsShared = kTRUE;
1952 p2->fIsShared = kTRUE;
1959 for (Int_t i=0;i<4;i++){
1960 if (s1->fOverlapLabels[3*i]==0){
1961 s1->fOverlapLabels[3*i] = s2->GetLabel();
1962 s1->fOverlapLabels[3*i+1] = sumshared;
1963 s1->fOverlapLabels[3*i+2] = s2->GetUniqueID();
1967 for (Int_t i=0;i<4;i++){
1968 if (s2->fOverlapLabels[3*i]==0){
1969 s2->fOverlapLabels[3*i] = s1->GetLabel();
1970 s2->fOverlapLabels[3*i+1] = sumshared;
1971 s2->fOverlapLabels[3*i+2] = s1->GetUniqueID();
1979 void AliTPCtrackerMI::SignShared(TObjArray * arr)
1982 //sort trackss according sectors
1984 for (Int_t i=0; i<arr->GetEntriesFast(); i++) {
1985 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
1987 //if (pt) RotateToLocal(pt);
1991 arr->Sort(); // sorting according z
1992 arr->Expand(arr->GetEntries());
1995 Int_t nseed=arr->GetEntriesFast();
1996 for (Int_t i=0; i<nseed; i++) {
1997 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
1999 for (Int_t j=0;j<=12;j++){
2000 pt->fOverlapLabels[j] =0;
2003 for (Int_t i=0; i<nseed; i++) {
2004 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2006 if (pt->fRemoval>10) continue;
2007 for (Int_t j=i+1; j<nseed; j++){
2008 AliTPCseed *pt2=(AliTPCseed*)arr->UncheckedAt(j);
2010 if (pt2->fRemoval<=10) {
2011 if ( TMath::Abs(pt->fRelativeSector-pt2->fRelativeSector)>0) break;
2018 void AliTPCtrackerMI::RemoveDouble(TObjArray * arr, Float_t factor1, Float_t factor2, Int_t removalindex)
2021 //sort trackss according sectors
2024 Info("RemoveDouble","Number of tracks before double removal- %d\n",arr->GetEntries());
2027 for (Int_t i=0; i<arr->GetEntriesFast(); i++) {
2028 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2033 arr->Sort(); // sorting according z
2034 arr->Expand(arr->GetEntries());
2036 //reset overlap labels
2038 Int_t nseed=arr->GetEntriesFast();
2039 for (Int_t i=0; i<nseed; i++) {
2040 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2043 for (Int_t j=0;j<=12;j++){
2044 pt->fOverlapLabels[j] =0;
2048 //sign shared tracks
2049 for (Int_t i=0; i<nseed; i++) {
2050 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2052 if (pt->fRemoval>10) continue;
2053 Float_t deltac = pt->GetC()*0.1;
2054 for (Int_t j=i+1; j<nseed; j++){
2055 AliTPCseed *pt2=(AliTPCseed*)arr->UncheckedAt(j);
2057 if (pt2->fRemoval<=10) {
2058 if ( TMath::Abs(pt->fRelativeSector-pt2->fRelativeSector)>0) break;
2059 if (TMath::Abs(pt->GetC() -pt2->GetC())>deltac) continue;
2060 if (TMath::Abs(pt->GetTgl()-pt2->GetTgl())>0.05) continue;
2067 // remove highly shared tracks
2068 for (Int_t i=0; i<nseed; i++) {
2069 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2071 if (pt->fRemoval>10) continue;
2074 for (Int_t j=0;j<4;j++){
2075 sumshared = pt->fOverlapLabels[j*3+1];
2077 Float_t factor = factor1;
2078 if (pt->fRemoval>0) factor = factor2;
2079 if (sumshared/pt->GetNumberOfClusters()>factor){
2080 for (Int_t j=0;j<4;j++){
2081 if (pt->fOverlapLabels[3*j]==0) continue;
2082 if (pt->fOverlapLabels[3*j+1]<5) continue;
2083 if (pt->fRemoval==removalindex) continue;
2084 AliTPCseed * pt2 = (AliTPCseed*)arr->UncheckedAt(pt->fOverlapLabels[3*j+2]);
2086 if (pt2->GetSigma2C()<pt->GetSigma2C()){
2087 // pt->fRemoval = removalindex;
2088 delete arr->RemoveAt(i);
2096 Info("RemoveDouble","Number of tracks after double removal- %d\n",arr->GetEntries());
2105 void AliTPCtrackerMI::SortTracks(TObjArray * arr, Int_t mode) const
2108 //sort tracks in array according mode criteria
2109 Int_t nseed = arr->GetEntriesFast();
2110 for (Int_t i=0; i<nseed; i++) {
2111 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2121 void AliTPCtrackerMI::RemoveUsed(TObjArray * arr, Float_t factor1, Float_t factor2, Int_t removalindex)
2124 //Loop over all tracks and remove "overlaps"
2127 Int_t nseed = arr->GetEntriesFast();
2130 for (Int_t i=0; i<nseed; i++) {
2131 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2133 delete arr->RemoveAt(i);
2137 pt->fBSigned = kFALSE;
2141 nseed = arr->GetEntriesFast();
2148 for (Int_t i=0; i<nseed; i++) {
2149 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2153 Int_t found,foundable,shared;
2155 pt->GetClusterStatistic(0,160,found, foundable,shared,kFALSE);
2157 pt->GetClusterStatistic(0,160,found, foundable,shared,kTRUE);
2159 Double_t factor = factor2;
2160 if (pt->fBConstrain) factor = factor1;
2162 if ((Float_t(shared)/Float_t(found))>factor){
2163 pt->Desactivate(removalindex);
2168 for (Int_t i=0; i<160; i++) {
2169 Int_t index=pt->GetClusterIndex2(i);
2170 if (index<0 || index&0x8000 ) continue;
2171 AliTPCclusterMI *c= pt->fClusterPointer[i];
2173 // if (!c->IsUsed(10)) c->Use(10);
2174 //if (pt->IsActive())
2183 Info("RemoveUsed","\n*****\nNumber of good tracks after shared removal\t%d\n",fNtracks);
2187 void AliTPCtrackerMI::UnsignClusters()
2190 // loop over all clusters and unsign them
2193 for (Int_t sec=0;sec<fkNIS;sec++){
2194 for (Int_t row=0;row<fInnerSec->GetNRows();row++){
2195 AliTPCclusterMI *cl = fInnerSec[sec][row].fClusters1;
2196 for (Int_t icl =0;icl< fInnerSec[sec][row].fN1;icl++)
2197 // if (cl[icl].IsUsed(10))
2199 cl = fInnerSec[sec][row].fClusters2;
2200 for (Int_t icl =0;icl< fInnerSec[sec][row].fN2;icl++)
2201 //if (cl[icl].IsUsed(10))
2206 for (Int_t sec=0;sec<fkNOS;sec++){
2207 for (Int_t row=0;row<fOuterSec->GetNRows();row++){
2208 AliTPCclusterMI *cl = fOuterSec[sec][row].fClusters1;
2209 for (Int_t icl =0;icl< fOuterSec[sec][row].fN1;icl++)
2210 //if (cl[icl].IsUsed(10))
2212 cl = fOuterSec[sec][row].fClusters2;
2213 for (Int_t icl =0;icl< fOuterSec[sec][row].fN2;icl++)
2214 //if (cl[icl].IsUsed(10))
2223 void AliTPCtrackerMI::SignClusters(TObjArray * arr, Float_t fnumber, Float_t fdensity)
2226 //sign clusters to be "used"
2228 // snumber and sdensity sign number of sigmas - bellow mean value to be accepted
2229 // loop over "primaries"
2243 Int_t nseed = arr->GetEntriesFast();
2244 for (Int_t i=0; i<nseed; i++) {
2245 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2249 if (!(pt->IsActive())) continue;
2250 Float_t dens = pt->GetNumberOfClusters()/Float_t(pt->fNFoundable);
2251 if ( (dens>0.7) && (pt->GetNumberOfClusters()>70)){
2253 sumdens2+= dens*dens;
2254 sumn += pt->GetNumberOfClusters();
2255 sumn2 += pt->GetNumberOfClusters()*pt->GetNumberOfClusters();
2256 Float_t chi2 = pt->GetChi2()/pt->GetNumberOfClusters();
2259 sumchi2 +=chi2*chi2;
2264 Float_t mdensity = 0.9;
2265 Float_t meann = 130;
2266 Float_t meanchi = 1;
2267 Float_t sdensity = 0.1;
2268 Float_t smeann = 10;
2269 Float_t smeanchi =0.4;
2273 mdensity = sumdens/sum;
2275 meanchi = sumchi/sum;
2277 sdensity = sumdens2/sum-mdensity*mdensity;
2278 sdensity = TMath::Sqrt(sdensity);
2280 smeann = sumn2/sum-meann*meann;
2281 smeann = TMath::Sqrt(smeann);
2283 smeanchi = sumchi2/sum - meanchi*meanchi;
2284 smeanchi = TMath::Sqrt(smeanchi);
2288 //REMOVE SHORT DELTAS or tracks going out of sensitive volume of TPC
2290 for (Int_t i=0; i<nseed; i++) {
2291 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2295 if (pt->fBSigned) continue;
2296 if (pt->fBConstrain) continue;
2297 //if (!(pt->IsActive())) continue;
2299 Int_t found,foundable,shared;
2300 pt->GetClusterStatistic(0,160,found, foundable,shared);
2301 if (shared/float(found)>0.3) {
2302 if (shared/float(found)>0.9 ){
2303 //delete arr->RemoveAt(i);
2308 Bool_t isok =kFALSE;
2309 if ( (pt->fNShared/pt->GetNumberOfClusters()<0.5) &&pt->GetNumberOfClusters()>60)
2311 if ((TMath::Abs(1/pt->GetC())<100.) && (pt->fNShared/pt->GetNumberOfClusters()<0.7))
2313 if (TMath::Abs(pt->GetZ()/pt->GetX())>1.1)
2315 if ( (TMath::Abs(pt->GetSnp()>0.7) && pt->GetD(0,0)>60.))
2319 for (Int_t i=0; i<160; i++) {
2320 Int_t index=pt->GetClusterIndex2(i);
2321 if (index<0) continue;
2322 AliTPCclusterMI *c= pt->fClusterPointer[i];
2324 //if (!(c->IsUsed(10))) c->Use();
2331 Double_t maxchi = meanchi+2.*smeanchi;
2333 for (Int_t i=0; i<nseed; i++) {
2334 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2338 //if (!(pt->IsActive())) continue;
2339 if (pt->fBSigned) continue;
2340 Double_t chi = pt->GetChi2()/pt->GetNumberOfClusters();
2341 if (chi>maxchi) continue;
2344 Float_t dens = pt->GetNumberOfClusters()/Float_t(pt->fNFoundable);
2346 //sign only tracks with enoug big density at the beginning
2348 if ((pt->GetDensityFirst(40)<0.75) && pt->GetNumberOfClusters()<meann) continue;
2351 Double_t mindens = TMath::Max(double(mdensity-sdensity*fdensity*bfactor),0.65);
2352 Double_t minn = TMath::Max(Int_t(meann-fnumber*smeann*bfactor),50);
2354 // if (pt->fBConstrain) mindens = TMath::Max(mdensity-sdensity*fdensity*bfactor,0.65);
2355 if ( (pt->fRemoval==10) && (pt->GetSnp()>0.8)&&(dens>mindens))
2358 if ((dens>mindens && pt->GetNumberOfClusters()>minn) && chi<maxchi ){
2359 //Int_t noc=pt->GetNumberOfClusters();
2360 pt->fBSigned = kTRUE;
2361 for (Int_t i=0; i<160; i++) {
2363 Int_t index=pt->GetClusterIndex2(i);
2364 if (index<0) continue;
2365 AliTPCclusterMI *c= pt->fClusterPointer[i];
2367 // if (!(c->IsUsed(10))) c->Use();
2372 // gLastCheck = nseed;
2380 void AliTPCtrackerMI::StopNotActive(TObjArray * arr, Int_t row0, Float_t th0, Float_t th1, Float_t th2) const
2382 // stop not active tracks
2383 // take th1 as threshold for number of founded to number of foundable on last 10 active rows
2384 // take th2 as threshold for number of founded to number of foundable on last 20 active rows
2385 Int_t nseed = arr->GetEntriesFast();
2387 for (Int_t i=0; i<nseed; i++) {
2388 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2392 if (!(pt->IsActive())) continue;
2393 StopNotActive(pt,row0,th0, th1,th2);
2399 void AliTPCtrackerMI::StopNotActive(AliTPCseed * seed, Int_t row0, Float_t th0, Float_t th1,
2402 // stop not active tracks
2403 // take th1 as threshold for number of founded to number of foundable on last 10 active rows
2404 // take th2 as threshold for number of founded to number of foundable on last 20 active rows
2407 Int_t foundable = 0;
2408 Int_t maxindex = seed->fLastPoint; //last foundable row
2409 if (seed->fNFoundable*th0 > seed->GetNumberOfClusters()) {
2410 seed->Desactivate(10) ;
2414 for (Int_t i=row0; i<maxindex; i++){
2415 Int_t index = seed->GetClusterIndex2(i);
2416 if (index!=-1) foundable++;
2418 if (foundable<=30) sumgood1++;
2419 if (foundable<=50) {
2426 if (foundable>=30.){
2427 if (sumgood1<(th1*30.)) seed->Desactivate(10);
2430 if (sumgood2<(th2*50.)) seed->Desactivate(10);
2434 Int_t AliTPCtrackerMI::RefitInward(AliESD *event)
2437 // back propagation of ESD tracks
2443 //PrepareForProlongation(fSeeds,1);
2444 PropagateForward2(fSeeds);
2446 Int_t nseed = fSeeds->GetEntriesFast();
2447 for (Int_t i=0;i<nseed;i++){
2448 AliTPCseed * seed = (AliTPCseed*) fSeeds->UncheckedAt(i);
2449 if (!seed) continue;
2450 seed->PropagateTo(fParam->GetInnerRadiusLow());
2451 AliESDtrack *esd=event->GetTrack(i);
2452 seed->CookdEdx(0.02,0.6);
2453 CookLabel(seed,0.1); //For comparison only
2454 if (seed->GetNumberOfClusters()>60){
2455 esd->UpdateTrackParams(seed,AliESDtrack::kTPCrefit);
2459 //printf("problem\n");
2462 Info("RefitInward","Number of refitted tracks %d",ntracks);
2469 Int_t AliTPCtrackerMI::PropagateBack(AliESD *event)
2472 // back propagation of ESD tracks
2478 PropagateBack(fSeeds);
2479 Int_t nseed = fSeeds->GetEntriesFast();
2481 for (Int_t i=0;i<nseed;i++){
2482 AliTPCseed * seed = (AliTPCseed*) fSeeds->UncheckedAt(i);
2483 if (!seed) continue;
2484 AliESDtrack *esd=event->GetTrack(i);
2485 seed->CookdEdx(0.02,0.6);
2486 CookLabel(seed,0.1); //For comparison only
2487 if (seed->GetNumberOfClusters()>60){
2488 esd->UpdateTrackParams(seed,AliESDtrack::kTPCout);
2492 Info("PropagateBack","Number of back propagated tracks %d",ntracks);
2499 void AliTPCtrackerMI::DeleteSeeds()
2503 Int_t nseed = fSeeds->GetEntriesFast();
2504 for (Int_t i=0;i<nseed;i++){
2505 AliTPCseed * seed = (AliTPCseed*)fSeeds->At(i);
2506 if (seed) delete fSeeds->RemoveAt(i);
2512 void AliTPCtrackerMI::ReadSeeds(AliESD *event, Int_t direction)
2515 //read seeds from the event
2517 Int_t nentr=event->GetNumberOfTracks();
2519 Info("ReadSeeds", "Number of ESD tracks: %d\n", nentr);
2524 fSeeds = new TObjArray(nentr);
2528 for (Int_t i=0; i<nentr; i++) {
2529 AliESDtrack *esd=event->GetTrack(i);
2530 ULong_t status=esd->GetStatus();
2531 AliTPCtrack t(*esd);
2532 AliTPCseed *seed = new AliTPCseed(t,t.GetAlpha());
2533 if ((status==AliESDtrack::kTPCin)&&(direction==1)) seed->ResetCovariance();
2534 if ( direction ==2 &&(status & AliESDtrack::kTRDrefit) == 0 ) seed->ResetCovariance();
2535 if ( direction ==2 && ((status & AliESDtrack::kTPCout) == 0) ) {
2540 if ( direction ==2 &&(status & AliESDtrack::kTRDrefit) > 0 ) {
2541 Double_t par0[5],par1[5],x;
2542 esd->GetInnerExternalParameters(x,par0);
2543 esd->GetExternalParameters(x,par1);
2544 Double_t delta1 = TMath::Abs(par0[4]-par1[4])/(0.000000001+TMath::Abs(par0[4]+par1[4]));
2545 Double_t delta2 = TMath::Abs(par0[3]-par1[3]);
2546 //reset covariance if suspicious
2547 if ( (delta1>0.1) || (delta2>0.006))
2548 seed->ResetCovariance();
2553 // rotate to the local coordinate system
2555 fSectors=fInnerSec; fN=fkNIS;
2557 Double_t alpha=seed->GetAlpha() - fSectors->GetAlphaShift();
2558 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
2559 if (alpha < 0. ) alpha += 2.*TMath::Pi();
2560 Int_t ns=Int_t(alpha/fSectors->GetAlpha())%fN;
2561 alpha =ns*fSectors->GetAlpha() + fSectors->GetAlphaShift();
2562 if (alpha<-TMath::Pi()) alpha += 2*TMath::Pi();
2563 if (alpha>=TMath::Pi()) alpha -= 2*TMath::Pi();
2564 alpha-=seed->GetAlpha();
2565 if (!seed->Rotate(alpha)) {
2571 //seed->PropagateTo(fSectors->GetX(0));
2573 // Int_t index = esd->GetTPCindex();
2574 //AliTPCseed * seed2= (AliTPCseed*)fSeeds->At(index);
2575 //if (direction==2){
2576 // AliTPCseed * seed2 = ReSeed(seed,0.,0.5,1.);
2584 for (Int_t irow=0;irow<160;irow++){
2585 Int_t index = seed->GetClusterIndex2(irow);
2588 AliTPCclusterMI * cl = GetClusterMI(index);
2589 seed->fClusterPointer[irow] = cl;
2591 if ((index & 0x8000)==0){
2592 cl->Use(10); // accepted cluster
2594 cl->Use(6); // close cluster not accepted
2597 Info("ReadSeeds","Not found cluster");
2601 fSeeds->AddAt(seed,i);
2607 //_____________________________________________________________________________
2608 void AliTPCtrackerMI::MakeSeeds3(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, Float_t cuts[4],
2609 Float_t deltay, Int_t ddsec) {
2610 //-----------------------------------------------------------------
2611 // This function creates track seeds.
2612 // SEEDING WITH VERTEX CONSTRAIN
2613 //-----------------------------------------------------------------
2614 // cuts[0] - fP4 cut
2615 // cuts[1] - tan(phi) cut
2616 // cuts[2] - zvertex cut
2617 // cuts[3] - fP3 cut
2625 Double_t x[5], c[15];
2626 // Int_t di = i1-i2;
2628 AliTPCseed * seed = new AliTPCseed;
2629 Double_t alpha=fSectors->GetAlpha(), shift=fSectors->GetAlphaShift();
2630 Double_t cs=cos(alpha), sn=sin(alpha);
2632 // Double_t x1 =fOuterSec->GetX(i1);
2633 //Double_t xx2=fOuterSec->GetX(i2);
2635 Double_t x1 =GetXrow(i1);
2636 Double_t xx2=GetXrow(i2);
2638 Double_t x3=GetX(), y3=GetY(), z3=GetZ();
2640 Int_t imiddle = (i2+i1)/2; //middle pad row index
2641 Double_t xm = GetXrow(imiddle); // radius of middle pad-row
2642 const AliTPCRow& krm=GetRow(sec,imiddle); //middle pad -row
2646 const AliTPCRow& kr1=GetRow(ns,i1);
2647 Double_t ymax = GetMaxY(i1)-kr1.fDeadZone-1.5;
2648 Double_t ymaxm = GetMaxY(imiddle)-kr1.fDeadZone-1.5;
2651 // change cut on curvature if it can't reach this layer
2652 // maximal curvature set to reach it
2653 Double_t dvertexmax = TMath::Sqrt((x1-x3)*(x1-x3)+(ymax+5-y3)*(ymax+5-y3));
2654 if (dvertexmax*0.5*cuts[0]>0.85){
2655 cuts[0] = 0.85/(dvertexmax*0.5+1.);
2657 Double_t r2min = 1/(cuts[0]*cuts[0]); //minimal square of radius given by cut
2660 if (deltay>0) ddsec = 0;
2661 // loop over clusters
2662 for (Int_t is=0; is < kr1; is++) {
2664 if (kr1[is]->IsUsed(10)) continue;
2665 Double_t y1=kr1[is]->GetY(), z1=kr1[is]->GetZ();
2666 //if (TMath::Abs(y1)>ymax) continue;
2668 if (deltay>0 && TMath::Abs(ymax-TMath::Abs(y1))> deltay ) continue; // seed only at the edge
2670 // find possible directions
2671 Float_t anglez = (z1-z3)/(x1-x3);
2672 Float_t extraz = z1 - anglez*(x1-xx2); // extrapolated z
2675 //find rotation angles relative to line given by vertex and point 1
2676 Double_t dvertex2 = (x1-x3)*(x1-x3)+(y1-y3)*(y1-y3);
2677 Double_t dvertex = TMath::Sqrt(dvertex2);
2678 Double_t angle13 = TMath::ATan((y1-y3)/(x1-x3));
2679 Double_t cs13 = cos(-angle13), sn13 = sin(-angle13);
2682 // loop over 2 sectors
2688 Double_t dddz1=0; // direction of delta inclination in z axis
2695 for (Int_t dsec = dsec1; dsec<=dsec2;dsec++){
2696 Int_t sec2 = sec + dsec;
2698 // AliTPCRow& kr2 = fOuterSec[(sec2+fkNOS)%fkNOS][i2];
2699 //AliTPCRow& kr2m = fOuterSec[(sec2+fkNOS)%fkNOS][imiddle];
2700 AliTPCRow& kr2 = GetRow((sec2+fkNOS)%fkNOS,i2);
2701 AliTPCRow& kr2m = GetRow((sec2+fkNOS)%fkNOS,imiddle);
2702 Int_t index1 = TMath::Max(kr2.Find(extraz-0.6-dddz1*TMath::Abs(z1)*0.05)-1,0);
2703 Int_t index2 = TMath::Min(kr2.Find(extraz+0.6+dddz2*TMath::Abs(z1)*0.05)+1,kr2);
2705 // rotation angles to p1-p3
2706 Double_t cs13r = cos(-angle13+dsec*alpha)/dvertex, sn13r = sin(-angle13+dsec*alpha)/dvertex;
2707 Double_t x2, y2, z2;
2709 // Double_t dymax = maxangle*TMath::Abs(x1-xx2);
2712 Double_t dxx0 = (xx2-x3)*cs13r;
2713 Double_t dyy0 = (xx2-x3)*sn13r;
2714 for (Int_t js=index1; js < index2; js++) {
2715 const AliTPCclusterMI *kcl = kr2[js];
2716 if (kcl->IsUsed(10)) continue;
2718 //calcutate parameters
2720 Double_t yy0 = dyy0 +(kcl->GetY()-y3)*cs13r;
2722 if (TMath::Abs(yy0)<0.000001) continue;
2723 Double_t xx0 = dxx0 -(kcl->GetY()-y3)*sn13r;
2724 Double_t y0 = 0.5*(xx0*xx0+yy0*yy0-xx0)/yy0;
2725 Double_t r02 = (0.25+y0*y0)*dvertex2;
2726 //curvature (radius) cut
2727 if (r02<r2min) continue;
2731 Double_t c0 = 1/TMath::Sqrt(r02);
2735 //Double_t dfi0 = 2.*TMath::ASin(dvertex*c0*0.5);
2736 //Double_t dfi1 = 2.*TMath::ASin(TMath::Sqrt(yy0*yy0+(1-xx0)*(1-xx0))*dvertex*c0*0.5);
2737 Double_t dfi0 = 2.*AliTPCFastMath::FastAsin(dvertex*c0*0.5);
2738 Double_t dfi1 = 2.*AliTPCFastMath::FastAsin(TMath::Sqrt(yy0*yy0+(1-xx0)*(1-xx0))*dvertex*c0*0.5);
2741 Double_t z0 = kcl->GetZ();
2742 Double_t zzzz2 = z1-(z1-z3)*dfi1/dfi0;
2743 if (TMath::Abs(zzzz2-z0)>0.5) continue;
2746 Double_t dip = (z1-z0)*c0/dfi1;
2747 Double_t x0 = (0.5*cs13+y0*sn13)*dvertex*c0;
2758 x2= xx2*cs-y2*sn*dsec;
2759 y2=+xx2*sn*dsec+y2*cs;
2769 // do we have cluster at the middle ?
2771 GetProlongation(x1,xm,x,ym,zm);
2773 AliTPCclusterMI * cm=0;
2774 if (TMath::Abs(ym)-ymaxm<0){
2775 cm = krm.FindNearest2(ym,zm,1.0,0.6,dummy);
2776 if ((!cm) || (cm->IsUsed(10))) {
2781 // rotate y1 to system 0
2782 // get state vector in rotated system
2783 Double_t yr1 = (-0.5*sn13+y0*cs13)*dvertex*c0;
2784 Double_t xr2 = x0*cs+yr1*sn*dsec;
2785 Double_t xr[5]={kcl->GetY(),kcl->GetZ(), xr2, dip, c0};
2787 GetProlongation(xx2,xm,xr,ym,zm);
2788 if (TMath::Abs(ym)-ymaxm<0){
2789 cm = kr2m.FindNearest2(ym,zm,1.0,0.6,dummy);
2790 if ((!cm) || (cm->IsUsed(10))) {
2800 dym = ym - cm->GetY();
2801 dzm = zm - cm->GetZ();
2808 Double_t sy1=kr1[is]->GetSigmaY2()*2., sz1=kr1[is]->GetSigmaZ2()*2.;
2809 Double_t sy2=kcl->GetSigmaY2()*2., sz2=kcl->GetSigmaZ2()*2.;
2810 //Double_t sy3=400*3./12., sy=0.1, sz=0.1;
2811 Double_t sy3=25000*x[4]*x[4]+0.1, sy=0.1, sz=0.1;
2812 //Double_t sy3=25000*x[4]*x[4]*60+0.5, sy=0.1, sz=0.1;
2814 Double_t f40=(F1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
2815 Double_t f42=(F1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
2816 Double_t f43=(F1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
2817 Double_t f20=(F2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
2818 Double_t f22=(F2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
2819 Double_t f23=(F2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
2821 Double_t f30=(F3(x1,y1+sy,x2,y2,z1,z2)-x[3])/sy;
2822 Double_t f31=(F3(x1,y1,x2,y2,z1+sz,z2)-x[3])/sz;
2823 Double_t f32=(F3(x1,y1,x2,y2+sy,z1,z2)-x[3])/sy;
2824 Double_t f34=(F3(x1,y1,x2,y2,z1,z2+sz)-x[3])/sz;
2828 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
2829 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
2830 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
2831 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
2832 c[13]=f30*sy1*f40+f32*sy2*f42;
2833 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
2835 // if (!BuildSeed(kr1[is],kcl,0,x1,x2,x3,x,c)) continue;
2837 UInt_t index=kr1.GetIndex(is);
2838 AliTPCseed *track=new(seed) AliTPCseed(index, x, c, x1, ns*alpha+shift);
2840 track->fIsSeeding = kTRUE;
2847 FollowProlongation(*track, (i1+i2)/2,1);
2848 Int_t foundable,found,shared;
2849 track->GetClusterStatistic((i1+i2)/2,i1, found, foundable, shared, kTRUE);
2850 if ((found<0.55*foundable) || shared>0.5*found || (track->GetSigmaY2()+track->GetSigmaZ2())>0.5){
2852 seed->~AliTPCseed();
2858 FollowProlongation(*track, i2,1);
2862 track->fBConstrain =1;
2863 // track->fLastPoint = i1+fInnerSec->GetNRows(); // first cluster in track position
2864 track->fLastPoint = i1; // first cluster in track position
2865 track->fFirstPoint = track->fLastPoint;
2867 if (track->GetNumberOfClusters()<(i1-i2)*0.5 ||
2868 track->GetNumberOfClusters() < track->fNFoundable*0.6 ||
2869 track->fNShared>0.4*track->GetNumberOfClusters() ) {
2871 seed->~AliTPCseed();
2875 // Z VERTEX CONDITION
2877 zv = track->GetZ()+track->GetTgl()/track->GetC()*
2878 ( asin(-track->GetEta()) - asin(track->GetX()*track->GetC()-track->GetEta()));
2879 if (TMath::Abs(zv-z3)>cuts[2]) {
2880 FollowProlongation(*track, TMath::Max(i2-20,0));
2881 zv = track->GetZ()+track->GetTgl()/track->GetC()*
2882 ( asin(-track->GetEta()) - asin(track->GetX()*track->GetC()-track->GetEta()));
2883 if (TMath::Abs(zv-z3)>cuts[2]){
2884 FollowProlongation(*track, TMath::Max(i2-40,0));
2885 zv = track->GetZ()+track->GetTgl()/track->GetC()*
2886 ( asin(-track->GetEta()) - asin(track->GetX()*track->GetC()-track->GetEta()));
2887 if (TMath::Abs(zv-z3)>cuts[2] &&(track->GetNumberOfClusters() > track->fNFoundable*0.7)){
2888 // make seed without constrain
2889 AliTPCseed * track2 = MakeSeed(track,0.2,0.5,1.);
2890 FollowProlongation(*track2, i2,1);
2891 track2->fBConstrain = kFALSE;
2892 track2->fSeedType = 1;
2893 arr->AddLast(track2);
2895 seed->~AliTPCseed();
2900 seed->~AliTPCseed();
2907 track->fSeedType =0;
2908 arr->AddLast(track);
2909 seed = new AliTPCseed;
2911 // don't consider other combinations
2912 if (track->GetNumberOfClusters() > track->fNFoundable*0.8)
2918 Info("MakeSeeds3","\nSeeding statistic:\t%d\t%d\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin,nout1,nout2);
2924 void AliTPCtrackerMI::MakeSeeds5(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, Float_t cuts[4],
2929 //-----------------------------------------------------------------
2930 // This function creates track seeds.
2931 //-----------------------------------------------------------------
2932 // cuts[0] - fP4 cut
2933 // cuts[1] - tan(phi) cut
2934 // cuts[2] - zvertex cut
2935 // cuts[3] - fP3 cut
2945 Double_t x[5], c[15];
2947 // make temporary seed
2948 AliTPCseed * seed = new AliTPCseed;
2949 Double_t alpha=fOuterSec->GetAlpha(), shift=fOuterSec->GetAlphaShift();
2950 // Double_t cs=cos(alpha), sn=sin(alpha);
2955 Double_t x1 = GetXrow(i1-1);
2956 const AliTPCRow& kr1=GetRow(sec,i1-1);
2957 Double_t y1max = GetMaxY(i1-1)-kr1.fDeadZone-1.5;
2959 Double_t x1p = GetXrow(i1);
2960 const AliTPCRow& kr1p=GetRow(sec,i1);
2962 Double_t x1m = GetXrow(i1-2);
2963 const AliTPCRow& kr1m=GetRow(sec,i1-2);
2966 //last 3 padrow for seeding
2967 AliTPCRow& kr3 = GetRow((sec+fkNOS)%fkNOS,i1-7);
2968 Double_t x3 = GetXrow(i1-7);
2969 // Double_t y3max= GetMaxY(i1-7)-kr3.fDeadZone-1.5;
2971 AliTPCRow& kr3p = GetRow((sec+fkNOS)%fkNOS,i1-6);
2972 Double_t x3p = GetXrow(i1-6);
2974 AliTPCRow& kr3m = GetRow((sec+fkNOS)%fkNOS,i1-8);
2975 Double_t x3m = GetXrow(i1-8);
2980 Int_t im = i1-4; //middle pad row index
2981 Double_t xm = GetXrow(im); // radius of middle pad-row
2982 const AliTPCRow& krm=GetRow(sec,im); //middle pad -row
2983 // Double_t ymmax = GetMaxY(im)-kr1.fDeadZone-1.5;
2986 Double_t deltax = x1-x3;
2987 Double_t dymax = deltax*cuts[1];
2988 Double_t dzmax = deltax*cuts[3];
2990 // loop over clusters
2991 for (Int_t is=0; is < kr1; is++) {
2993 if (kr1[is]->IsUsed(10)) continue;
2994 Double_t y1=kr1[is]->GetY(), z1=kr1[is]->GetZ();
2996 if (deltay>0 && TMath::Abs(y1max-TMath::Abs(y1))> deltay ) continue; // seed only at the edge
2998 Int_t index1 = TMath::Max(kr3.Find(z1-dzmax)-1,0);
2999 Int_t index2 = TMath::Min(kr3.Find(z1+dzmax)+1,kr3);
3005 for (Int_t js=index1; js < index2; js++) {
3006 const AliTPCclusterMI *kcl = kr3[js];
3007 if (kcl->IsUsed(10)) continue;
3009 // apply angular cuts
3010 if (TMath::Abs(y1-y3)>dymax) continue;
3013 if (TMath::Abs(z1-z3)>dzmax) continue;
3015 Double_t angley = (y1-y3)/(x1-x3);
3016 Double_t anglez = (z1-z3)/(x1-x3);
3018 Double_t erry = TMath::Abs(angley)*(x1-x1m)*0.5+0.5;
3019 Double_t errz = TMath::Abs(anglez)*(x1-x1m)*0.5+0.5;
3021 Double_t yyym = angley*(xm-x1)+y1;
3022 Double_t zzzm = anglez*(xm-x1)+z1;
3024 const AliTPCclusterMI *kcm = krm.FindNearest2(yyym,zzzm,erry,errz,index);
3026 if (kcm->IsUsed(10)) continue;
3028 erry = TMath::Abs(angley)*(x1-x1m)*0.4+0.5;
3029 errz = TMath::Abs(anglez)*(x1-x1m)*0.4+0.5;
3036 // look around first
3037 const AliTPCclusterMI *kc1m = kr1m.FindNearest2(angley*(x1m-x1)+y1,
3043 if (kc1m->IsUsed(10)) used++;
3045 const AliTPCclusterMI *kc1p = kr1p.FindNearest2(angley*(x1p-x1)+y1,
3051 if (kc1p->IsUsed(10)) used++;
3053 if (used>1) continue;
3054 if (found<1) continue;
3058 const AliTPCclusterMI *kc3m = kr3m.FindNearest2(angley*(x3m-x3)+y3,
3064 if (kc3m->IsUsed(10)) used++;
3068 const AliTPCclusterMI *kc3p = kr3p.FindNearest2(angley*(x3p-x3)+y3,
3074 if (kc3p->IsUsed(10)) used++;
3078 if (used>1) continue;
3079 if (found<3) continue;
3089 x[4]=F1(x1,y1,x2,y2,x3,y3);
3090 //if (TMath::Abs(x[4]) >= cuts[0]) continue;
3093 x[2]=F2(x1,y1,x2,y2,x3,y3);
3096 x[3]=F3n(x1,y1,x2,y2,z1,z2,x[4]);
3097 //if (TMath::Abs(x[3]) > cuts[3]) continue;
3101 Double_t sy1=0.1, sz1=0.1;
3102 Double_t sy2=0.1, sz2=0.1;
3103 Double_t sy3=0.1, sy=0.1, sz=0.1;
3105 Double_t f40=(F1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
3106 Double_t f42=(F1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
3107 Double_t f43=(F1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
3108 Double_t f20=(F2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
3109 Double_t f22=(F2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
3110 Double_t f23=(F2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
3112 Double_t f30=(F3(x1,y1+sy,x2,y2,z1,z2)-x[3])/sy;
3113 Double_t f31=(F3(x1,y1,x2,y2,z1+sz,z2)-x[3])/sz;
3114 Double_t f32=(F3(x1,y1,x2,y2+sy,z1,z2)-x[3])/sy;
3115 Double_t f34=(F3(x1,y1,x2,y2,z1,z2+sz)-x[3])/sz;
3119 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3120 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3121 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3122 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3123 c[13]=f30*sy1*f40+f32*sy2*f42;
3124 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3126 // if (!BuildSeed(kr1[is],kcl,0,x1,x2,x3,x,c)) continue;
3128 UInt_t index=kr1.GetIndex(is);
3129 AliTPCseed *track=new(seed) AliTPCseed(index, x, c, x1, sec*alpha+shift);
3131 track->fIsSeeding = kTRUE;
3134 FollowProlongation(*track, i1-7,1);
3135 if (track->GetNumberOfClusters() < track->fNFoundable*0.75 ||
3136 track->fNShared>0.6*track->GetNumberOfClusters() || ( track->GetSigmaY2()+ track->GetSigmaZ2())>0.6){
3138 seed->~AliTPCseed();
3144 FollowProlongation(*track, i2,1);
3145 track->fBConstrain =0;
3146 track->fLastPoint = i1+fInnerSec->GetNRows(); // first cluster in track position
3147 track->fFirstPoint = track->fLastPoint;
3149 if (track->GetNumberOfClusters()<(i1-i2)*0.5 ||
3150 track->GetNumberOfClusters()<track->fNFoundable*0.7 ||
3151 track->fNShared>2. || track->GetChi2()/track->GetNumberOfClusters()>6 || ( track->GetSigmaY2()+ track->GetSigmaZ2())>0.5 ) {
3153 seed->~AliTPCseed();
3158 FollowProlongation(*track, TMath::Max(i2-10,0),1);
3159 AliTPCseed * track2 = MakeSeed(track,0.2,0.5,0.9);
3160 FollowProlongation(*track2, i2,1);
3161 track2->fBConstrain = kFALSE;
3162 track2->fSeedType = 4;
3163 arr->AddLast(track2);
3165 seed->~AliTPCseed();
3169 //arr->AddLast(track);
3170 //seed = new AliTPCseed;
3176 Info("MakeSeeds5","\nSeeding statiistic:\t%d\t%d\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin,nout1,nout2,nout3);
3182 //_____________________________________________________________________________
3183 void AliTPCtrackerMI::MakeSeeds2(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, Float_t */*cuts[4]*/,
3184 Float_t deltay, Bool_t /*bconstrain*/) {
3185 //-----------------------------------------------------------------
3186 // This function creates track seeds - without vertex constraint
3187 //-----------------------------------------------------------------
3188 // cuts[0] - fP4 cut - not applied
3189 // cuts[1] - tan(phi) cut
3190 // cuts[2] - zvertex cut - not applied
3191 // cuts[3] - fP3 cut
3201 Double_t alpha=fOuterSec->GetAlpha(), shift=fOuterSec->GetAlphaShift();
3202 // Double_t cs=cos(alpha), sn=sin(alpha);
3203 Int_t row0 = (i1+i2)/2;
3204 Int_t drow = (i1-i2)/2;
3205 const AliTPCRow& kr0=fSectors[sec][row0];
3208 AliTPCpolyTrack polytrack;
3209 Int_t nclusters=fSectors[sec][row0];
3210 AliTPCseed * seed = new AliTPCseed;
3215 for (Int_t is=0; is < nclusters; is++) { //LOOP over clusters
3217 Int_t nfoundable =0;
3218 for (Int_t iter =1; iter<2; iter++){ //iterations
3219 const AliTPCRow& krm=fSectors[sec][row0-iter];
3220 const AliTPCRow& krp=fSectors[sec][row0+iter];
3221 const AliTPCclusterMI * cl= kr0[is];
3223 if (cl->IsUsed(10)) {
3229 Double_t x = kr0.GetX();
3230 // Initialization of the polytrack
3235 Double_t y0= cl->GetY();
3236 Double_t z0= cl->GetZ();
3240 Double_t ymax = fSectors->GetMaxY(row0)-kr0.fDeadZone-1.5;
3241 if (deltay>0 && TMath::Abs(ymax-TMath::Abs(y0))> deltay ) continue; // seed only at the edge
3243 erry = (0.5)*cl->GetSigmaY2()/TMath::Sqrt(cl->GetQ())*6;
3244 errz = (0.5)*cl->GetSigmaZ2()/TMath::Sqrt(cl->GetQ())*6;
3245 polytrack.AddPoint(x,y0,z0,erry, errz);
3248 if (cl->IsUsed(10)) sumused++;
3251 Float_t roady = (5*TMath::Sqrt(cl->GetSigmaY2()+0.2)+1.)*iter;
3252 Float_t roadz = (5*TMath::Sqrt(cl->GetSigmaZ2()+0.2)+1.)*iter;
3255 AliTPCclusterMI * cl1 = krm.FindNearest(y0,z0,roady,roadz);
3256 if (cl1 && TMath::Abs(ymax-TMath::Abs(y0))) {
3257 erry = (0.5)*cl1->GetSigmaY2()/TMath::Sqrt(cl1->GetQ())*3;
3258 errz = (0.5)*cl1->GetSigmaZ2()/TMath::Sqrt(cl1->GetQ())*3;
3259 if (cl1->IsUsed(10)) sumused++;
3260 polytrack.AddPoint(x,cl1->GetY(),cl1->GetZ(),erry,errz);
3264 AliTPCclusterMI * cl2 = krp.FindNearest(y0,z0,roady,roadz);
3266 erry = (0.5)*cl2->GetSigmaY2()/TMath::Sqrt(cl2->GetQ())*3;
3267 errz = (0.5)*cl2->GetSigmaZ2()/TMath::Sqrt(cl2->GetQ())*3;
3268 if (cl2->IsUsed(10)) sumused++;
3269 polytrack.AddPoint(x,cl2->GetY(),cl2->GetZ(),erry,errz);
3272 if (sumused>0) continue;
3274 polytrack.UpdateParameters();
3280 nfoundable = polytrack.GetN();
3281 nfound = nfoundable;
3283 for (Int_t ddrow = iter+1; ddrow<drow;ddrow++){
3284 Float_t maxdist = 0.8*(1.+3./(ddrow));
3285 for (Int_t delta = -1;delta<=1;delta+=2){
3286 Int_t row = row0+ddrow*delta;
3287 kr = &(fSectors[sec][row]);
3288 Double_t xn = kr->GetX();
3289 Double_t ymax = fSectors->GetMaxY(row)-kr->fDeadZone-1.5;
3290 polytrack.GetFitPoint(xn,yn,zn);
3291 if (TMath::Abs(yn)>ymax) continue;
3293 AliTPCclusterMI * cln = kr->FindNearest(yn,zn,roady,roadz);
3295 Float_t dist = TMath::Sqrt( (yn-cln->GetY())*(yn-cln->GetY())+(zn-cln->GetZ())*(zn-cln->GetZ()));
3298 erry = (dist+0.3)*cln->GetSigmaY2()/TMath::Sqrt(cln->GetQ())*(1.+1./(ddrow));
3299 errz = (dist+0.3)*cln->GetSigmaZ2()/TMath::Sqrt(cln->GetQ())*(1.+1./(ddrow));
3300 if (cln->IsUsed(10)) {
3301 // printf("used\n");
3309 polytrack.AddPoint(xn,cln->GetY(),cln->GetZ(),erry, errz);
3314 if ( (sumused>3) || (sumused>0.5*nfound) || (nfound<0.6*nfoundable)) break;
3315 polytrack.UpdateParameters();
3318 if ( (sumused>3) || (sumused>0.5*nfound)) {
3319 //printf("sumused %d\n",sumused);
3324 polytrack.GetFitDerivation(kr0.GetX(),dy,dz);
3325 AliTPCpolyTrack track2;
3327 polytrack.Refit(track2,0.5+TMath::Abs(dy)*0.3,0.4+TMath::Abs(dz)*0.3);
3328 if (track2.GetN()<0.5*nfoundable) continue;
3331 if ((nfound>0.6*nfoundable) &&( nfoundable>0.4*(i1-i2))) {
3333 // test seed with and without constrain
3334 for (Int_t constrain=0; constrain<=0;constrain++){
3335 // add polytrack candidate
3337 Double_t x[5], c[15];
3338 Double_t x1,x2,x3,y1,y2,y3,z1,z2,z3;
3339 track2.GetBoundaries(x3,x1);
3341 track2.GetFitPoint(x1,y1,z1);
3342 track2.GetFitPoint(x2,y2,z2);
3343 track2.GetFitPoint(x3,y3,z3);
3345 //is track pointing to the vertex ?
3348 polytrack.GetFitPoint(x0,y0,z0);
3361 x[4]=F1(x1,y1,x2,y2,x3,y3);
3363 // if (TMath::Abs(x[4]) >= cuts[0]) continue; //
3364 x[2]=F2(x1,y1,x2,y2,x3,y3);
3366 //if (TMath::Abs(x[4]*x1-x[2]) >= cuts[1]) continue;
3367 //x[3]=F3(x1,y1,x2,y2,z1,z2);
3368 x[3]=F3n(x1,y1,x3,y3,z1,z3,x[4]);
3369 //if (TMath::Abs(x[3]) > cuts[3]) continue;
3372 Double_t sy =0.1, sz =0.1;
3373 Double_t sy1=0.02, sz1=0.02;
3374 Double_t sy2=0.02, sz2=0.02;
3378 sy3=25000*x[4]*x[4]+0.1, sy=0.1, sz=0.1;
3381 Double_t f40=(F1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
3382 Double_t f42=(F1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
3383 Double_t f43=(F1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
3384 Double_t f20=(F2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
3385 Double_t f22=(F2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
3386 Double_t f23=(F2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
3388 Double_t f30=(F3(x1,y1+sy,x3,y3,z1,z3)-x[3])/sy;
3389 Double_t f31=(F3(x1,y1,x3,y3,z1+sz,z3)-x[3])/sz;
3390 Double_t f32=(F3(x1,y1,x3,y3+sy,z1,z3)-x[3])/sy;
3391 Double_t f34=(F3(x1,y1,x3,y3,z1,z3+sz)-x[3])/sz;
3396 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3397 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3398 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3399 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3400 c[13]=f30*sy1*f40+f32*sy2*f42;
3401 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3403 //Int_t row1 = fSectors->GetRowNumber(x1);
3404 Int_t row1 = GetRowNumber(x1);
3408 AliTPCseed *track=new (seed) AliTPCseed(index, x, c, x1, sec*alpha+shift);
3409 track->fIsSeeding = kTRUE;
3410 Int_t rc=FollowProlongation(*track, i2);
3411 if (constrain) track->fBConstrain =1;
3413 track->fBConstrain =0;
3414 track->fLastPoint = row1+fInnerSec->GetNRows(); // first cluster in track position
3415 track->fFirstPoint = track->fLastPoint;
3417 if (rc==0 || track->GetNumberOfClusters()<(i1-i2)*0.5 ||
3418 track->GetNumberOfClusters() < track->fNFoundable*0.6 ||
3419 track->fNShared>0.4*track->GetNumberOfClusters()) {
3422 seed->~AliTPCseed();
3425 arr->AddLast(track);
3426 seed = new AliTPCseed;
3430 } // if accepted seed
3433 Info("MakeSeeds2","\nSeeding statiistic:\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin3);
3439 AliTPCseed *AliTPCtrackerMI::MakeSeed(AliTPCseed *track, Float_t r0, Float_t r1, Float_t r2)
3443 //reseed using track points
3444 Int_t p0 = int(r0*track->GetNumberOfClusters()); // point 0
3445 Int_t p1 = int(r1*track->GetNumberOfClusters());
3446 Int_t p2 = int(r2*track->GetNumberOfClusters()); // last point
3448 Double_t x0[3],x1[3],x2[3];
3453 // find track position at given ratio of the length
3454 Int_t sec0, sec1, sec2;
3460 for (Int_t i=0;i<160;i++){
3461 if (track->fClusterPointer[i]){
3463 AliTPCTrackerPoint *trpoint =track->GetTrackPoint(i);
3464 if ( (index<p0) || x0[0]<0 ){
3465 if (trpoint->GetX()>1){
3466 clindex = track->GetClusterIndex2(i);
3468 x0[0] = trpoint->GetX();
3469 x0[1] = trpoint->GetY();
3470 x0[2] = trpoint->GetZ();
3471 sec0 = ((clindex&0xff000000)>>24)%18;
3476 if ( (index<p1) &&(trpoint->GetX()>1)){
3477 clindex = track->GetClusterIndex2(i);
3479 x1[0] = trpoint->GetX();
3480 x1[1] = trpoint->GetY();
3481 x1[2] = trpoint->GetZ();
3482 sec1 = ((clindex&0xff000000)>>24)%18;
3485 if ( (index<p2) &&(trpoint->GetX()>1)){
3486 clindex = track->GetClusterIndex2(i);
3488 x2[0] = trpoint->GetX();
3489 x2[1] = trpoint->GetY();
3490 x2[2] = trpoint->GetZ();
3491 sec2 = ((clindex&0xff000000)>>24)%18;
3498 Double_t alpha, cs,sn, xx2,yy2;
3500 alpha = (sec1-sec2)*fSectors->GetAlpha();
3501 cs = TMath::Cos(alpha);
3502 sn = TMath::Sin(alpha);
3503 xx2= x1[0]*cs-x1[1]*sn;
3504 yy2= x1[0]*sn+x1[1]*cs;
3508 alpha = (sec0-sec2)*fSectors->GetAlpha();
3509 cs = TMath::Cos(alpha);
3510 sn = TMath::Sin(alpha);
3511 xx2= x0[0]*cs-x0[1]*sn;
3512 yy2= x0[0]*sn+x0[1]*cs;
3518 Double_t x[5],c[15];
3522 x[4]=F1(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]);
3523 // if (x[4]>1) return 0;
3524 x[2]=F2(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]);
3525 x[3]=F3n(x2[0],x2[1],x0[0],x0[1],x2[2],x0[2],x[4]);
3526 //if (TMath::Abs(x[3]) > 2.2) return 0;
3527 //if (TMath::Abs(x[2]) > 1.99) return 0;
3529 Double_t sy =0.1, sz =0.1;
3531 Double_t sy1=0.02+track->GetSigmaY2(), sz1=0.02+track->GetSigmaZ2();
3532 Double_t sy2=0.01+track->GetSigmaY2(), sz2=0.01+track->GetSigmaZ2();
3533 Double_t sy3=0.01+track->GetSigmaY2();
3535 Double_t f40=(F1(x2[0],x2[1]+sy,x1[0],x1[1],x0[0],x0[1])-x[4])/sy;
3536 Double_t f42=(F1(x2[0],x2[1],x1[0],x1[1]+sy,x0[0],x0[1])-x[4])/sy;
3537 Double_t f43=(F1(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]+sy)-x[4])/sy;
3538 Double_t f20=(F2(x2[0],x2[1]+sy,x1[0],x1[1],x0[0],x0[1])-x[2])/sy;
3539 Double_t f22=(F2(x2[0],x2[1],x1[0],x1[1]+sy,x0[0],x0[1])-x[2])/sy;
3540 Double_t f23=(F2(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]+sy)-x[2])/sy;
3542 Double_t f30=(F3(x2[0],x2[1]+sy,x0[0],x0[1],x2[2],x0[2])-x[3])/sy;
3543 Double_t f31=(F3(x2[0],x2[1],x0[0],x0[1],x2[2]+sz,x0[2])-x[3])/sz;
3544 Double_t f32=(F3(x2[0],x2[1],x0[0],x0[1]+sy,x2[2],x0[2])-x[3])/sy;
3545 Double_t f34=(F3(x2[0],x2[1],x0[0],x0[1],x2[2],x0[2]+sz)-x[3])/sz;
3550 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3551 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3552 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3553 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3554 c[13]=f30*sy1*f40+f32*sy2*f42;
3555 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3557 // Int_t row1 = fSectors->GetRowNumber(x2[0]);
3558 AliTPCseed *seed=new AliTPCseed(0, x, c, x2[0], sec2*fSectors->GetAlpha()+fSectors->GetAlphaShift());
3559 // Double_t y0,z0,y1,z1, y2,z2;
3560 //seed->GetProlongation(x0[0],y0,z0);
3561 // seed->GetProlongation(x1[0],y1,z1);
3562 //seed->GetProlongation(x2[0],y2,z2);
3564 seed->fLastPoint = pp2;
3565 seed->fFirstPoint = pp2;
3572 AliTPCseed *AliTPCtrackerMI::ReSeed(AliTPCseed *track, Float_t r0, Float_t r1, Float_t r2)
3576 //reseed using founded clusters
3578 // Find the number of clusters
3579 Int_t nclusters = 0;
3580 for (Int_t irow=0;irow<160;irow++){
3581 if (track->GetClusterIndex(irow)>0) nclusters++;
3585 ipos[0] = TMath::Max(int(r0*nclusters),0); // point 0 cluster
3586 ipos[1] = TMath::Min(int(r1*nclusters),nclusters-1); //
3587 ipos[2] = TMath::Min(int(r2*nclusters),nclusters-1); // last point
3591 Int_t row[3],sec[3]={0,0,0};
3593 // find track row position at given ratio of the length
3595 for (Int_t irow=0;irow<160;irow++){
3596 if (track->GetClusterIndex2(irow)<0) continue;
3598 for (Int_t ipoint=0;ipoint<3;ipoint++){
3599 if (index<=ipos[ipoint]) row[ipoint] = irow;
3603 //Get cluster and sector position
3604 for (Int_t ipoint=0;ipoint<3;ipoint++){
3605 Int_t clindex = track->GetClusterIndex2(row[ipoint]);
3606 AliTPCclusterMI * cl = GetClusterMI(clindex);
3609 // AliTPCclusterMI * cl = GetClusterMI(clindex);
3612 sec[ipoint] = ((clindex&0xff000000)>>24)%18;
3613 xyz[ipoint][0] = GetXrow(row[ipoint]);
3614 xyz[ipoint][1] = cl->GetY();
3615 xyz[ipoint][2] = cl->GetZ();
3619 // Calculate seed state vector and covariance matrix
3621 Double_t alpha, cs,sn, xx2,yy2;
3623 alpha = (sec[1]-sec[2])*fSectors->GetAlpha();
3624 cs = TMath::Cos(alpha);
3625 sn = TMath::Sin(alpha);
3626 xx2= xyz[1][0]*cs-xyz[1][1]*sn;
3627 yy2= xyz[1][0]*sn+xyz[1][1]*cs;
3631 alpha = (sec[0]-sec[2])*fSectors->GetAlpha();
3632 cs = TMath::Cos(alpha);
3633 sn = TMath::Sin(alpha);
3634 xx2= xyz[0][0]*cs-xyz[0][1]*sn;
3635 yy2= xyz[0][0]*sn+xyz[0][1]*cs;
3641 Double_t x[5],c[15];
3645 x[4]=F1(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
3646 x[2]=F2(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
3647 x[3]=F3n(xyz[2][0],xyz[2][1],xyz[0][0],xyz[0][1],xyz[2][2],xyz[0][2],x[4]);
3649 Double_t sy =0.1, sz =0.1;
3651 Double_t sy1=0.2, sz1=0.2;
3652 Double_t sy2=0.2, sz2=0.2;
3655 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;
3656 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;
3657 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;
3658 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;
3659 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;
3660 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;
3662 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;
3663 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;
3664 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;
3665 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;
3670 c[3]=f20*sy1; c[4]=0.; c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3671 c[6]=f30*sy1; c[7]=f31*sz1; c[8]=f30*sy1*f20+f32*sy2*f22;
3672 c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3673 c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3674 c[13]=f30*sy1*f40+f32*sy2*f42;
3675 c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3677 // Int_t row1 = fSectors->GetRowNumber(xyz[2][0]);
3678 AliTPCseed *seed=new AliTPCseed(0, x, c, xyz[2][0], sec[2]*fSectors->GetAlpha()+fSectors->GetAlphaShift());
3679 seed->fLastPoint = row[2];
3680 seed->fFirstPoint = row[2];
3684 Int_t AliTPCtrackerMI::CheckKinkPoint(AliTPCseed*seed, Float_t th)
3689 for (Int_t i=0;i<12;i++) seed->fKinkPoint[i]=0;
3691 if (TMath::Abs(seed->GetC())>0.01) return 0;
3694 Float_t x[160], y[160], erry[160], z[160], errz[160];
3696 Float_t xt[160], yt[160], zt[160];
3701 Int_t middle = seed->GetNumberOfClusters()/2;
3704 // find central sector, get local cooordinates
3706 for (Int_t i=seed->fFirstPoint;i<=seed->fLastPoint;i++) {
3707 sec[i]= seed->GetClusterSector(i)%18;
3710 AliTPCclusterMI * cl = seed->fClusterPointer[i];
3711 // if (cl==0) cl = GetClusterMI(seed->GetClusterIndex2(i));
3718 if (i>i2) i2 = i; //last point with cluster
3719 if (i2<i1) i1 = i; //first point with cluster
3722 AliTPCTrackerPoint * point = seed->GetTrackPoint(i);
3724 yt[i] = point->GetY();
3725 zt[i] = point->GetZ();
3727 if (point->GetX()>0){
3728 erry[i] = point->GetErrY();
3729 errz[i] = point->GetErrZ();
3734 secm = sec[i]; //central sector
3735 padm = i; //middle point with cluster
3740 // rotate position to global coordinate system connected to sector at last the point
3742 for (Int_t i=i1;i<=i2;i++){
3744 if (sec[i]<0) continue;
3745 Double_t alpha = (sec[i2]-sec[i])*fSectors->GetAlpha();
3746 Double_t cs = TMath::Cos(alpha);
3747 Double_t sn = TMath::Sin(alpha);
3748 Float_t xx2= x[i]*cs+y[i]*sn;
3749 Float_t yy2= -x[i]*sn+y[i]*cs;
3753 xx2= xt[i]*cs+yt[i]*sn;
3754 yy2= -xt[i]*sn+yt[i]*cs;
3759 //get "state" vector
3760 Double_t xh[5],xm = x[padm];
3763 xh[4]=F1(xt[i2],yt[i2],xt[padm],yt[padm],xt[i1],yt[i1]);
3764 xh[2]=F2(xt[i2],yt[i2],xt[padm],yt[padm],xt[i1],yt[i1]);
3765 xh[3]=F3n(xt[i2],yt[i2],xt[i1],yt[i1],zt[i2],zt[i1],xh[4]);
3768 for (Int_t i=i1;i<=i2;i++){
3770 if (sec[i]<0) continue;
3771 GetProlongation(x[i2], x[i],xh,yy,zz);
3772 if (TMath::Abs(y[i]-yy)>4||TMath::Abs(z[i]-zz)>4){
3774 //xxh[4]=F1old(x[i2],y[i2],x[padm],y[padm],x[i1],y[i1]);
3775 //xxh[2]=F2old(x[i2],y[i2],x[padm],y[padm],x[i1],y[i1]);
3776 Error("AliTPCtrackerMI::CheckKinkPoint","problem\n");
3781 Float_t dyup[160],dydown[160], dzup[160], dzdown[160];
3782 Float_t yup[160], ydown[160], zup[160], zdown[160];
3784 AliTPCpolyTrack ptrack1,ptrack2;
3787 for (Int_t i=i1;i<=i2;i++){
3788 AliTPCclusterMI * cl = seed->fClusterPointer[i];
3790 if (cl->GetType()<0) continue;
3791 if (cl->GetType()>10) continue;
3794 ptrack1.AddPoint(x[i]-xm,y[i],z[i],0.1,0.1);
3796 if (ptrack1.GetN()>4.){
3797 ptrack1.UpdateParameters();
3799 ptrack1.GetFitDerivation(x[i]-xm,ddy,ddz);
3801 ptrack1.GetFitPoint(x[i]-xm,yy,zz);
3810 dyup[i]=0.; //not enough points
3815 for (Int_t i=i2;i>=i1;i--){
3816 AliTPCclusterMI * cl = seed->fClusterPointer[i];
3818 if (cl->GetType()<0) continue;
3819 if (cl->GetType()>10) continue;
3821 ptrack2.AddPoint(x[i]-xm,y[i],z[i],0.1,0.1);
3823 if (ptrack2.GetN()>4){
3824 ptrack2.UpdateParameters();
3826 ptrack2.GetFitDerivation(x[i]-xm,ddy,ddz);
3828 ptrack2.GetFitPoint(x[i]-xm,yy,zz);
3836 dydown[i]=0.; //not enough points
3841 // find maximal difference of the derivation
3842 for (Int_t i=0;i<12;i++) seed->fKinkPoint[i]=0;
3845 for (Int_t i=i1+10;i<i2-10;i++){
3846 if ( (TMath::Abs(dydown[i])<0.00000001) || (TMath::Abs(dyup[i])<0.00000001) ||i<30)continue;
3847 // printf("%f\t%f\t%f\t%f\t%f\n",x[i],dydown[i],dyup[i],dzdown[i],dzup[i]);
3849 Float_t ddy = TMath::Abs(dydown[i]-dyup[i]);
3850 Float_t ddz = TMath::Abs(dzdown[i]-dzup[i]);
3851 if ( (ddy+ddz)> th){
3852 seed->fKinkPoint[0] = i;
3853 seed->fKinkPoint[1] = ddy;
3854 seed->fKinkPoint[2] = ddz;
3861 //write information to the debug tree
3862 TBranch * br = fTreeDebug->GetBranch("debug");
3863 TClonesArray * arr = new TClonesArray("AliTPCTrackPoint2");
3864 arr->ExpandCreateFast(i2-i1);
3865 br->SetAddress(&arr);
3867 AliTPCclusterMI cldummy;
3869 AliTPCTrackPoint2 pdummy;
3870 pdummy.GetTPoint().fIsShared = 10;
3872 Double_t alpha = sec[i2]*fSectors->GetAlpha();
3873 Double_t cs = TMath::Cos(alpha);
3874 Double_t sn = TMath::Sin(alpha);
3876 for (Int_t i=i1;i<i2;i++){
3877 AliTPCTrackPoint2 *trpoint = (AliTPCTrackPoint2*)arr->UncheckedAt(i-i1);
3879 AliTPCclusterMI * cl0 = seed->fClusterPointer[i];
3881 AliTPCTrackerPoint * point = seed->GetTrackPoint(i);
3884 Double_t x = GetXrow(i);
3885 trpoint->GetTPoint() = *point;
3886 trpoint->GetCPoint() = *cl0;
3887 trpoint->GetCPoint().SetQ(TMath::Abs(cl0->GetQ()));
3888 trpoint->fID = seed->GetUniqueID();
3889 trpoint->fLab = seed->GetLabel();
3891 trpoint->fGX = cs *x + sn*point->GetY();
3892 trpoint->fGY = -sn *x + cs*point->GetY() ;
3893 trpoint->fGZ = point->GetZ();
3895 trpoint->fDY = y[i];
3896 trpoint->fDZ = z[i];
3898 trpoint->fDYU = dyup[i];
3899 trpoint->fDZU = dzup[i];
3901 trpoint->fDYD = dydown[i];
3902 trpoint->fDZD = dzdown[i];
3904 if (TMath::Abs(dyup[i])>0.00000000001 &&TMath::Abs(dydown[i])>0.00000000001){
3905 trpoint->fDDY = dydown[i]-dyup[i];
3906 trpoint->fDDZ = dzdown[i]-dzup[i];
3914 trpoint->GetCPoint()= cldummy;
3931 AliTPCseed* AliTPCtrackerMI::ReSeed(AliTPCseed *t)
3934 // reseed - refit - track
3937 // Int_t last = fSectors->GetNRows()-1;
3939 if (fSectors == fOuterSec){
3940 first = TMath::Max(first, t->fFirstPoint-fInnerSec->GetNRows());
3944 first = t->fFirstPoint;
3946 AliTPCseed * seed = MakeSeed(t,0.1,0.5,0.9);
3947 FollowBackProlongation(*t,fSectors->GetNRows()-1);
3949 FollowProlongation(*t,first);
3959 //_____________________________________________________________________________
3960 Int_t AliTPCtrackerMI::ReadSeeds(const TFile *inp) {
3961 //-----------------------------------------------------------------
3962 // This function reades track seeds.
3963 //-----------------------------------------------------------------
3964 TDirectory *savedir=gDirectory;
3966 TFile *in=(TFile*)inp;
3967 if (!in->IsOpen()) {
3968 cerr<<"AliTPCtrackerMI::ReadSeeds(): input file is not open !\n";
3973 TTree *seedTree=(TTree*)in->Get("Seeds");
3975 cerr<<"AliTPCtrackerMI::ReadSeeds(): ";
3976 cerr<<"can't get a tree with track seeds !\n";
3979 AliTPCtrack *seed=new AliTPCtrack;
3980 seedTree->SetBranchAddress("tracks",&seed);
3982 if (fSeeds==0) fSeeds=new TObjArray(15000);
3984 Int_t n=(Int_t)seedTree->GetEntries();
3985 for (Int_t i=0; i<n; i++) {
3986 seedTree->GetEvent(i);
3987 fSeeds->AddLast(new AliTPCseed(*seed,seed->GetAlpha()));
3996 Int_t AliTPCtrackerMI::Clusters2Tracks (AliESD *esd)
3999 if (fSeeds) DeleteSeeds();
4002 if (!fSeeds) return 1;
4009 //_____________________________________________________________________________
4010 Int_t AliTPCtrackerMI::Clusters2Tracks() {
4011 //-----------------------------------------------------------------
4012 // This is a track finder.
4013 //-----------------------------------------------------------------
4014 TDirectory *savedir=gDirectory;
4018 fSeeds = Tracking();
4021 Info("Clusters2Tracks","Time for tracking: \t");timer.Print();timer.Start();
4023 //activate again some tracks
4024 for (Int_t i=0; i<fSeeds->GetEntriesFast(); i++) {
4025 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
4027 Int_t nc=t.GetNumberOfClusters();
4029 delete fSeeds->RemoveAt(i);
4032 if (pt->fRemoval==10) {
4033 if (pt->GetDensityFirst(20)>0.8 || pt->GetDensityFirst(30)>0.8 || pt->GetDensityFirst(40)>0.7)
4034 pt->Desactivate(10); // make track again active
4036 pt->Desactivate(20);
4037 delete fSeeds->RemoveAt(i);
4041 RemoveDouble(fSeeds,0.2,0.6,11);
4042 RemoveUsed(fSeeds,0.5,0.5,6);
4045 Int_t nseed=fSeeds->GetEntriesFast();
4047 for (Int_t i=0; i<nseed; i++) {
4048 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
4050 Int_t nc=t.GetNumberOfClusters();
4052 delete fSeeds->RemoveAt(i);
4055 CookLabel(pt,0.1); //For comparison only
4056 //if ((pt->IsActive() || (pt->fRemoval==10) )&& nc>50 &&pt->GetNumberOfClusters()>0.4*pt->fNFoundable){
4057 if ((pt->IsActive() || (pt->fRemoval==10) )){
4059 if (fDebug>0) cerr<<found<<'\r';
4063 delete fSeeds->RemoveAt(i);
4067 //RemoveOverlap(fSeeds,0.99,7,kTRUE);
4069 //RemoveUsed(fSeeds,0.9,0.9,6);
4071 nseed=fSeeds->GetEntriesFast();
4073 for (Int_t i=0; i<nseed; i++) {
4074 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
4076 Int_t nc=t.GetNumberOfClusters();
4078 delete fSeeds->RemoveAt(i);
4082 t.CookdEdx(0.02,0.6);
4083 // CheckKinkPoint(&t,0.05);
4084 //if ((pt->IsActive() || (pt->fRemoval==10) )&& nc>50 &&pt->GetNumberOfClusters()>0.4*pt->fNFoundable){
4085 if ((pt->IsActive() || (pt->fRemoval==10) )){
4093 delete fSeeds->RemoveAt(i);
4094 //AliTPCseed * seed1 = ReSeed(pt,0.05,0.5,1);
4096 // FollowProlongation(*seed1,0);
4097 // Int_t n = seed1->GetNumberOfClusters();
4098 // printf("fP4\t%f\t%f\n",seed1->GetC(),pt->GetC());
4099 // printf("fN\t%d\t%d\n", seed1->GetNumberOfClusters(),pt->GetNumberOfClusters());
4102 //AliTPCseed * seed2 = ReSeed(pt,0.95,0.5,0.05);
4106 SortTracks(fSeeds, 1);
4110 PrepareForBackProlongation(fSeeds,5.);
4111 PropagateBack(fSeeds);
4112 printf("Time for back propagation: \t");timer.Print();timer.Start();
4116 PrepareForProlongation(fSeeds,5.);
4117 PropagateForward2(fSeeds);
4119 printf("Time for FORWARD propagation: \t");timer.Print();timer.Start();
4120 // RemoveUsed(fSeeds,0.7,0.7,6);
4121 //RemoveOverlap(fSeeds,0.9,7,kTRUE);
4123 nseed=fSeeds->GetEntriesFast();
4125 for (Int_t i=0; i<nseed; i++) {
4126 AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
4128 Int_t nc=t.GetNumberOfClusters();
4130 delete fSeeds->RemoveAt(i);
4133 t.CookdEdx(0.02,0.6);
4134 // CookLabel(pt,0.1); //For comparison only
4135 //if ((pt->IsActive() || (pt->fRemoval==10) )&& nc>50 &&pt->GetNumberOfClusters()>0.4*pt->fNFoundable){
4136 if ((pt->IsActive() || (pt->fRemoval==10) )){
4137 cerr<<found++<<'\r';
4140 delete fSeeds->RemoveAt(i);
4145 // fNTracks = found;
4147 Info("Clusters2Tracks","Time for overlap removal, track writing and dedx cooking: \t"); timer.Print();timer.Start();
4150 // cerr<<"Number of found tracks : "<<"\t"<<found<<endl;
4151 Info("Clusters2Tracks","Number of found tracks %d",found);
4153 // UnloadClusters();
4158 void AliTPCtrackerMI::Tracking(TObjArray * arr)
4161 // tracking of the seeds
4164 fSectors = fOuterSec;
4165 ParallelTracking(arr,150,63);
4166 fSectors = fOuterSec;
4167 ParallelTracking(arr,63,0);
4170 TObjArray * AliTPCtrackerMI::Tracking(Int_t seedtype, Int_t i1, Int_t i2, Float_t cuts[4], Float_t dy, Int_t dsec)
4175 TObjArray * arr = new TObjArray;
4177 fSectors = fOuterSec;
4180 for (Int_t sec=0;sec<fkNOS;sec++){
4181 if (seedtype==3) MakeSeeds3(arr,sec,i1,i2,cuts,dy, dsec);
4182 if (seedtype==4) MakeSeeds5(arr,sec,i1,i2,cuts,dy);
4183 if (seedtype==2) MakeSeeds2(arr,sec,i1,i2,cuts,dy);
4186 Info("Tracking","\nSeeding - %d\t%d\t%d\t%d\n",seedtype,i1,i2,arr->GetEntriesFast());
4198 TObjArray * AliTPCtrackerMI::Tracking()
4204 Int_t nup=fOuterSec->GetNRows()+fInnerSec->GetNRows();
4206 TObjArray * seeds = new TObjArray;
4215 Float_t fnumber = 3.0;
4216 Float_t fdensity = 3.0;
4221 for (Int_t delta = 0; delta<18; delta+=6){
4225 arr = Tracking(3,nup-1-delta,nup-1-delta-gap,cuts,-1,1);
4226 SumTracks(seeds,arr);
4227 SignClusters(seeds,fnumber,fdensity);
4229 for (Int_t i=2;i<6;i+=2){
4230 // seed high pt tracks
4233 arr = Tracking(3,nup-i-delta,nup-i-delta-gap,cuts,-1,0);
4234 SumTracks(seeds,arr);
4235 SignClusters(seeds,fnumber,fdensity);
4240 // RemoveUsed(seeds,0.9,0.9,1);
4241 // UnsignClusters();
4242 // SignClusters(seeds,fnumber,fdensity);
4246 for (Int_t delta = 20; delta<120; delta+=10){
4248 // seed high pt tracks
4252 arr = Tracking(3,nup-delta,nup-delta-gap,cuts,-1);
4253 SumTracks(seeds,arr);
4254 SignClusters(seeds,fnumber,fdensity);
4259 arr = Tracking(3,nup-delta-5,nup-delta-5-gap,cuts,-1);
4260 SumTracks(seeds,arr);
4261 SignClusters(seeds,fnumber,fdensity);
4272 Info("Tracking()","\n\nPrimary seeding\t%d\n\n",seeds->GetEntriesFast());
4276 // RemoveUsed(seeds,0.75,0.75,1);
4278 //SignClusters(seeds,fnumber,fdensity);
4287 arr = Tracking(4,nup-1,nup-1-gap,cuts,-1);
4288 SumTracks(seeds,arr);
4289 SignClusters(seeds,fnumber,fdensity);
4291 arr = Tracking(4,nup-2,nup-2-gap,cuts,-1);
4292 SumTracks(seeds,arr);
4293 SignClusters(seeds,fnumber,fdensity);
4295 arr = Tracking(4,nup-3,nup-3-gap,cuts,-1);
4296 SumTracks(seeds,arr);
4297 SignClusters(seeds,fnumber,fdensity);
4301 for (Int_t delta = 3; delta<30; delta+=5){
4307 arr = Tracking(4,nup-1-delta,nup-1-delta-gap,cuts,-1);
4308 SumTracks(seeds,arr);
4309 SignClusters(seeds,fnumber,fdensity);
4311 arr = Tracking(4,nup-3-delta,nup-5-delta-gap,cuts,4);
4312 SumTracks(seeds,arr);
4313 SignClusters(seeds,fnumber,fdensity);
4325 for (Int_t delta = 30; delta<70; delta+=10){
4331 arr = Tracking(4,nup-1-delta,nup-1-delta-gap,cuts,-1);
4332 SumTracks(seeds,arr);
4333 SignClusters(seeds,fnumber,fdensity);
4335 arr = Tracking(4,nup-5-delta,nup-5-delta-gap,cuts,5 );
4336 SumTracks(seeds,arr);
4337 SignClusters(seeds,fnumber,fdensity);
4341 Info("Tracking()","\n\nSecondary seeding\t%d\n\n",seeds->GetEntriesFast());
4352 void AliTPCtrackerMI::SumTracks(TObjArray *arr1,TObjArray *arr2) const
4355 //sum tracks to common container
4356 //remove suspicious tracks
4357 Int_t nseed = arr2->GetEntriesFast();
4358 for (Int_t i=0;i<nseed;i++){
4359 AliTPCseed *pt=(AliTPCseed*)arr2->UncheckedAt(i);
4362 // NORMAL ACTIVE TRACK
4363 if (pt->IsActive()){
4364 arr1->AddLast(arr2->RemoveAt(i));
4367 //remove not usable tracks
4368 if (pt->fRemoval!=10){
4369 delete arr2->RemoveAt(i);
4372 // REMOVE VERY SHORT TRACKS
4373 if (pt->GetNumberOfClusters()<20){
4374 delete arr2->RemoveAt(i);
4377 // ENABLE ONLY ENOUGH GOOD STOPPED TRACKS
4378 if (pt->GetDensityFirst(20)>0.8 || pt->GetDensityFirst(30)>0.8 || pt->GetDensityFirst(40)>0.7)
4379 arr1->AddLast(arr2->RemoveAt(i));
4381 delete arr2->RemoveAt(i);
4390 void AliTPCtrackerMI::ParallelTracking(TObjArray * arr, Int_t rfirst, Int_t rlast)
4393 // try to track in parralel
4395 Int_t nseed=arr->GetEntriesFast();
4396 //prepare seeds for tracking
4397 for (Int_t i=0; i<nseed; i++) {
4398 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i), &t=*pt;
4400 if (!t.IsActive()) continue;
4401 // follow prolongation to the first layer
4402 if ( (fSectors ==fInnerSec) || (t.fFirstPoint-fParam->GetNRowLow()>rfirst+1) )
4403 FollowProlongation(t, rfirst+1);
4408 for (Int_t nr=rfirst; nr>=rlast; nr--){
4409 if (nr<fInnerSec->GetNRows())
4410 fSectors = fInnerSec;
4412 fSectors = fOuterSec;
4413 // make indexes with the cluster tracks for given
4415 // find nearest cluster
4416 for (Int_t i=0; i<nseed; i++) {
4417 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i), &t=*pt;
4419 if (!pt->IsActive()) continue;
4420 // if ( (fSectors ==fOuterSec) && (pt->fFirstPoint-fParam->GetNRowLow())<nr) continue;
4421 if (pt->fRelativeSector>17) {
4424 UpdateClusters(t,nr);
4426 // prolonagate to the nearest cluster - if founded
4427 for (Int_t i=0; i<nseed; i++) {
4428 AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
4430 if (!pt->IsActive()) continue;
4431 // if ((fSectors ==fOuterSec) && (pt->fFirstPoint-fParam->GetNRowLow())<nr) continue;
4432 if (pt->fRelativeSector>17) {
4435 FollowToNextCluster(*pt,nr);
4440 void AliTPCtrackerMI::PrepareForBackProlongation(TObjArray * arr,Float_t fac) const
4444 // if we use TPC track itself we have to "update" covariance
4446 Int_t nseed= arr->GetEntriesFast();
4447 for (Int_t i=0;i<nseed;i++){
4448 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
4452 //rotate to current local system at first accepted point
4453 Int_t index = pt->GetClusterIndex2(pt->fFirstPoint);
4454 Int_t sec = (index&0xff000000)>>24;
4456 Float_t angle1 = fInnerSec->GetAlpha()*sec+fInnerSec->GetAlphaShift();
4457 if (angle1>TMath::Pi())
4458 angle1-=2.*TMath::Pi();
4459 Float_t angle2 = pt->GetAlpha();
4461 if (TMath::Abs(angle1-angle2)>0.001){
4462 pt->Rotate(angle1-angle2);
4463 //angle2 = pt->GetAlpha();
4464 //pt->fRelativeSector = pt->GetAlpha()/fInnerSec->GetAlpha();
4465 //if (pt->GetAlpha()<0)
4466 // pt->fRelativeSector+=18;
4467 //sec = pt->fRelativeSector;
4476 void AliTPCtrackerMI::PrepareForProlongation(TObjArray * arr, Float_t fac) const
4480 // if we use TPC track itself we have to "update" covariance
4482 Int_t nseed= arr->GetEntriesFast();
4483 for (Int_t i=0;i<nseed;i++){
4484 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
4487 pt->fFirstPoint = pt->fLastPoint;
4495 Int_t AliTPCtrackerMI::PropagateBack(TObjArray * arr)
4498 // make back propagation
4500 Int_t nseed= arr->GetEntriesFast();
4501 for (Int_t i=0;i<nseed;i++){
4502 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
4504 //AliTPCseed *pt2 = new AliTPCseed(*pt);
4505 fSectors = fInnerSec;
4506 //FollowBackProlongation(*pt,fInnerSec->GetNRows()-1);
4507 //fSectors = fOuterSec;
4508 FollowBackProlongation(*pt,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1);
4509 //if (pt->GetNumberOfClusters()<(pt->fEsd->GetTPCclusters(0)) ){
4510 // Error("PropagateBack","Not prolonged track %d",pt->GetLabel());
4511 // FollowBackProlongation(*pt2,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1);
4519 Int_t AliTPCtrackerMI::PropagateForward2(TObjArray * arr)
4522 // make forward propagation
4524 Int_t nseed= arr->GetEntriesFast();
4526 for (Int_t i=0;i<nseed;i++){
4527 AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
4529 FollowProlongation(*pt,0);
4536 Int_t AliTPCtrackerMI::PropagateForward()
4539 // propagate track forward
4541 Int_t nseed = fSeeds->GetEntriesFast();
4542 for (Int_t i=0;i<nseed;i++){
4543 AliTPCseed *pt = (AliTPCseed*)fSeeds->UncheckedAt(i);
4545 AliTPCseed &t = *pt;
4546 Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
4547 if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
4548 if (alpha < 0. ) alpha += 2.*TMath::Pi();
4549 t.fRelativeSector = Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN;
4553 fSectors = fOuterSec;
4554 ParallelTracking(fSeeds,fOuterSec->GetNRows()+fInnerSec->GetNRows()-1,fInnerSec->GetNRows());
4555 fSectors = fInnerSec;
4556 ParallelTracking(fSeeds,fInnerSec->GetNRows()-1,0);
4566 Int_t AliTPCtrackerMI::PropagateBack(AliTPCseed * pt, Int_t row0, Int_t row1)
4569 // make back propagation, in between row0 and row1
4573 fSectors = fInnerSec;
4576 if (row1<fSectors->GetNRows())
4579 r1 = fSectors->GetNRows()-1;
4581 if (row0<fSectors->GetNRows()&& r1>0 )
4582 FollowBackProlongation(*pt,r1);
4583 if (row1<=fSectors->GetNRows())
4586 r1 = row1 - fSectors->GetNRows();
4587 if (r1<=0) return 0;
4588 if (r1>=fOuterSec->GetNRows()) return 0;
4589 fSectors = fOuterSec;
4590 return FollowBackProlongation(*pt,r1);
4598 void AliTPCtrackerMI::GetShape(AliTPCseed * seed, Int_t row)
4602 Float_t sd2 = TMath::Abs((fParam->GetZLength()-TMath::Abs(seed->GetZ())))*fParam->GetDiffL()*fParam->GetDiffL();
4603 // Float_t padlength = fParam->GetPadPitchLength(seed->fSector);
4604 Float_t padlength = GetPadPitchLength(row);
4606 Float_t sresy = (seed->fSector < fParam->GetNSector()/2) ? 0.2 :0.3;
4607 Float_t angulary = seed->GetSnp();
4608 angulary = angulary*angulary/(1-angulary*angulary);
4609 seed->fCurrentSigmaY2 = sd2+padlength*padlength*angulary/12.+sresy*sresy;
4611 Float_t sresz = fParam->GetZSigma();
4612 Float_t angularz = seed->GetTgl();
4613 seed->fCurrentSigmaZ2 = sd2+padlength*padlength*angularz*angularz*(1+angulary)/12.+sresz*sresz;
4615 Float_t wy = GetSigmaY(seed);
4616 Float_t wz = GetSigmaZ(seed);
4619 if (TMath::Abs(wy/seed->fCurrentSigmaY2-1)>0.0001 || TMath::Abs(wz/seed->fCurrentSigmaZ2-1)>0.0001 ){
4620 printf("problem\n");
4626 Float_t AliTPCtrackerMI::GetSigmaY(AliTPCseed * seed)
4630 Float_t sd2 = TMath::Abs((fParam->GetZLength()-TMath::Abs(seed->GetZ())))*fParam->GetDiffL()*fParam->GetDiffL();
4631 Float_t padlength = fParam->GetPadPitchLength(seed->fSector);
4632 Float_t sres = (seed->fSector < fParam->GetNSector()/2) ? 0.2 :0.3;
4633 Float_t angular = seed->GetSnp();
4634 angular = angular*angular/(1-angular*angular);
4635 // angular*=angular;
4636 //angular = TMath::Sqrt(angular/(1-angular));
4637 Float_t res = TMath::Sqrt(sd2+padlength*padlength*angular/12.+sres*sres);
4640 Float_t AliTPCtrackerMI::GetSigmaZ(AliTPCseed * seed)
4644 Float_t sd2 = TMath::Abs((fParam->GetZLength()-TMath::Abs(seed->GetZ())))*fParam->GetDiffL()*fParam->GetDiffL();
4645 Float_t padlength = fParam->GetPadPitchLength(seed->fSector);
4646 Float_t sres = fParam->GetZSigma();
4647 Float_t angular = seed->GetTgl();
4648 Float_t res = TMath::Sqrt(sd2+padlength*padlength*angular*angular/12.+sres*sres);
4655 //__________________________________________________________________________
4656 void AliTPCtrackerMI::CookLabel(AliTPCseed *t, Float_t wrong) const {
4657 //--------------------------------------------------------------------
4658 //This function "cooks" a track label. If label<0, this track is fake.
4659 //--------------------------------------------------------------------
4660 Int_t noc=t->GetNumberOfClusters();
4662 //printf("\nnot founded prolongation\n\n\n");
4668 AliTPCclusterMI *clusters[160];
4670 for (Int_t i=0;i<160;i++) {
4677 for (i=0; i<160 && current<noc; i++) {
4679 Int_t index=t->GetClusterIndex2(i);
4680 if (index<=0) continue;
4681 if (index&0x8000) continue;
4683 //clusters[current]=GetClusterMI(index);
4684 if (t->fClusterPointer[i]){
4685 clusters[current]=t->fClusterPointer[i];
4691 Int_t lab=123456789;
4692 for (i=0; i<noc; i++) {
4693 AliTPCclusterMI *c=clusters[i];
4695 lab=TMath::Abs(c->GetLabel(0));
4697 for (j=0; j<noc; j++) if (lb[j]==lab || mx[j]==0) break;
4703 for (i=0; i<noc; i++) if (mx[i]>max) {max=mx[i]; lab=lb[i];}
4705 for (i=0; i<noc; i++) {
4706 AliTPCclusterMI *c=clusters[i];
4708 if (TMath::Abs(c->GetLabel(1)) == lab ||
4709 TMath::Abs(c->GetLabel(2)) == lab ) max++;
4712 if ((1.- Float_t(max)/noc) > wrong) lab=-lab;
4715 Int_t tail=Int_t(0.10*noc);
4718 for (i=1; i<=160&&ind<tail; i++) {
4719 // AliTPCclusterMI *c=clusters[noc-i];
4720 AliTPCclusterMI *c=clusters[i];
4722 if (lab == TMath::Abs(c->GetLabel(0)) ||
4723 lab == TMath::Abs(c->GetLabel(1)) ||
4724 lab == TMath::Abs(c->GetLabel(2))) max++;
4727 if (max < Int_t(0.5*tail)) lab=-lab;
4734 //delete[] clusters;
4738 Int_t AliTPCtrackerMI::AliTPCSector::GetRowNumber(Double_t x) const
4740 //return pad row number for this x
4743 r=fRow[fN-1].GetX();
4744 if (x > r) return fN;
4746 if (x < r) return -1;
4747 return Int_t((x-r)/fPadPitchLength + 0.5);}
4749 r=fRow[fN-1].GetX();
4750 if (x > r) return fN;
4752 if (x < r) return -1;
4753 Double_t r1=fRow[64].GetX();
4755 return Int_t((x-r)/f1PadPitchLength + 0.5);}
4757 return (Int_t((x-r1)/f2PadPitchLength + 0.5)+64);}
4761 //_________________________________________________________________________
4762 void AliTPCtrackerMI::AliTPCSector::Setup(const AliTPCParam *par, Int_t f) {
4763 //-----------------------------------------------------------------------
4764 // Setup inner sector
4765 //-----------------------------------------------------------------------
4767 fAlpha=par->GetInnerAngle();
4768 fAlphaShift=par->GetInnerAngleShift();
4769 fPadPitchWidth=par->GetInnerPadPitchWidth();
4770 fPadPitchLength=par->GetInnerPadPitchLength();
4771 fN=par->GetNRowLow();
4772 fRow=new AliTPCRow[fN];
4773 for (Int_t i=0; i<fN; i++) {
4774 fRow[i].SetX(par->GetPadRowRadiiLow(i));
4775 fRow[i].fDeadZone =1.5; //1.5 cm of dead zone
4778 fAlpha=par->GetOuterAngle();
4779 fAlphaShift=par->GetOuterAngleShift();
4780 fPadPitchWidth = par->GetOuterPadPitchWidth();
4781 fPadPitchLength = par->GetOuter1PadPitchLength();
4782 f1PadPitchLength = par->GetOuter1PadPitchLength();
4783 f2PadPitchLength = par->GetOuter2PadPitchLength();
4785 fN=par->GetNRowUp();
4786 fRow=new AliTPCRow[fN];
4787 for (Int_t i=0; i<fN; i++) {
4788 fRow[i].SetX(par->GetPadRowRadiiUp(i));
4789 fRow[i].fDeadZone =1.5; // 1.5 cm of dead zone
4794 AliTPCtrackerMI::AliTPCRow::AliTPCRow() {
4796 // default constructor
4804 AliTPCtrackerMI::AliTPCRow::~AliTPCRow(){
4811 //_________________________________________________________________________
4813 AliTPCtrackerMI::AliTPCRow::InsertCluster(const AliTPCclusterMI* c, UInt_t index) {
4814 //-----------------------------------------------------------------------
4815 // Insert a cluster into this pad row in accordence with its y-coordinate
4816 //-----------------------------------------------------------------------
4817 if (fN==kMaxClusterPerRow) {
4818 cerr<<"AliTPCRow::InsertCluster(): Too many clusters !\n"; return;
4820 if (fN==0) {fIndex[0]=index; fClusters[fN++]=c; return;}
4821 Int_t i=Find(c->GetZ());
4822 memmove(fClusters+i+1 ,fClusters+i,(fN-i)*sizeof(AliTPCclusterMI*));
4823 memmove(fIndex +i+1 ,fIndex +i,(fN-i)*sizeof(UInt_t));
4824 fIndex[i]=index; fClusters[i]=c; fN++;
4827 void AliTPCtrackerMI::AliTPCRow::ResetClusters() {
4833 //delete[] fClusterArray;
4834 if (fClusters1) delete []fClusters1;
4835 if (fClusters2) delete []fClusters2;
4842 //___________________________________________________________________
4843 Int_t AliTPCtrackerMI::AliTPCRow::Find(Double_t z) const {
4844 //-----------------------------------------------------------------------
4845 // Return the index of the nearest cluster
4846 //-----------------------------------------------------------------------
4847 if (fN==0) return 0;
4848 if (z <= fClusters[0]->GetZ()) return 0;
4849 if (z > fClusters[fN-1]->GetZ()) return fN;
4850 Int_t b=0, e=fN-1, m=(b+e)/2;
4851 for (; b<e; m=(b+e)/2) {
4852 if (z > fClusters[m]->GetZ()) b=m+1;
4860 //___________________________________________________________________
4861 AliTPCclusterMI * AliTPCtrackerMI::AliTPCRow::FindNearest(Double_t y, Double_t z, Double_t roady, Double_t roadz) const {
4862 //-----------------------------------------------------------------------
4863 // Return the index of the nearest cluster in z y
4864 //-----------------------------------------------------------------------
4865 Float_t maxdistance = roady*roady + roadz*roadz;
4867 AliTPCclusterMI *cl =0;
4868 for (Int_t i=Find(z-roadz); i<fN; i++) {
4869 AliTPCclusterMI *c=(AliTPCclusterMI*)(fClusters[i]);
4870 if (c->GetZ() > z+roadz) break;
4871 if ( (c->GetY()-y) > roady ) continue;
4872 Float_t distance = (c->GetZ()-z)*(c->GetZ()-z)+(c->GetY()-y)*(c->GetY()-y);
4873 if (maxdistance>distance) {
4874 maxdistance = distance;
4881 AliTPCclusterMI * AliTPCtrackerMI::AliTPCRow::FindNearest2(Double_t y, Double_t z, Double_t roady, Double_t roadz,UInt_t & index) const
4883 //-----------------------------------------------------------------------
4884 // Return the index of the nearest cluster in z y
4885 //-----------------------------------------------------------------------
4886 Float_t maxdistance = roady*roady + roadz*roadz;
4887 Int_t iz1 = TMath::Max(fFastCluster[Int_t(z-roadz+254.5)]-1,0);
4888 Int_t iz2 = TMath::Min(fFastCluster[Int_t(z+roadz+255.5)]+1,fN);
4890 AliTPCclusterMI *cl =0;
4891 //FindNearest3(y,z,roady,roadz,index);
4892 // for (Int_t i=Find(z-roadz); i<fN; i++) {
4893 for (Int_t i=iz1; i<iz2; i++) {
4894 AliTPCclusterMI *c=(AliTPCclusterMI*)(fClusters[i]);
4895 if (c->GetZ() > z+roadz) break;
4896 if ( c->GetY()-y > roady ) continue;
4897 if ( y-c->GetY() > roady ) continue;
4898 Float_t distance = (c->GetZ()-z)*(c->GetZ()-z)+(c->GetY()-y)*(c->GetY()-y);
4899 if (maxdistance>distance) {
4900 maxdistance = distance;
4903 //roady = TMath::Sqrt(maxdistance);
4911 AliTPCclusterMI * AliTPCtrackerMI::AliTPCRow::FindNearest3(Double_t y, Double_t z, Double_t roady, Double_t roadz,UInt_t & index) const
4913 //-----------------------------------------------------------------------
4914 // Return the index of the nearest cluster in z y
4915 //-----------------------------------------------------------------------
4916 Float_t maxdistance = roady*roady + roadz*roadz;
4917 // Int_t iz = Int_t(z+255.);
4918 AliTPCclusterMI *cl =0;
4919 for (Int_t i=Find(z-roadz); i<fN; i++) {
4920 //for (Int_t i=fFastCluster[iz-2]; i<fFastCluster[iz+2]; i++) {
4921 AliTPCclusterMI *c=(AliTPCclusterMI*)(fClusters[i]);
4922 if (c->GetZ() > z+roadz) break;
4923 if ( c->GetY()-y > roady ) continue;
4924 if ( y-c->GetY() > roady ) continue;
4925 Float_t distance = (c->GetZ()-z)*(c->GetZ()-z)+(c->GetY()-y)*(c->GetY()-y);
4926 if (maxdistance>distance) {
4927 maxdistance = distance;
4930 //roady = TMath::Sqrt(maxdistance);
4939 AliTPCseed::AliTPCseed():AliTPCtrack(){
4943 for (Int_t i=0;i<200;i++) SetClusterIndex2(i,-3);
4944 for (Int_t i=0;i<160;i++) fClusterPointer[i]=0;
4961 AliTPCseed::AliTPCseed(const AliTPCseed &s):AliTPCtrack(s){
4962 //---------------------
4963 // dummy copy constructor
4964 //-------------------------
4966 AliTPCseed::AliTPCseed(const AliTPCtrack &t):AliTPCtrack(t){
4975 for (Int_t i=0;i<160;i++) {
4976 fClusterPointer[i] = 0;
4977 Int_t index = t.GetClusterIndex(i);
4979 SetClusterIndex2(i,index);
4982 SetClusterIndex2(i,-3);
4995 AliTPCseed::AliTPCseed(const AliKalmanTrack &t, Double_t a):AliTPCtrack(t,a){
4999 for (Int_t i=0;i<160;i++) {
5000 fClusterPointer[i] = 0;
5001 Int_t index = t.GetClusterIndex(i);
5002 SetClusterIndex2(i,index);
5023 AliTPCseed::AliTPCseed(UInt_t index, const Double_t xx[5], const Double_t cc[15],
5024 Double_t xr, Double_t alpha):
5025 AliTPCtrack(index, xx, cc, xr, alpha) {
5030 for (Int_t i=0;i<200;i++) SetClusterIndex2(i,-3);
5031 for (Int_t i=0;i<160;i++) fClusterPointer[i]=0;
5040 // fHelixIn = new TClonesArray("AliHelix",0);
5041 //fHelixOut = new TClonesArray("AliHelix",0);
5051 AliTPCseed::~AliTPCseed(){
5054 if (fPoints) delete fPoints;
5056 if (fEPoints) delete fEPoints;
5061 AliTPCTrackerPoint * AliTPCseed::GetTrackPoint(Int_t i)
5065 return &fTrackPoints[i];
5068 void AliTPCseed::RebuildSeed()
5071 // rebuild seed to be ready for storing
5072 AliTPCclusterMI cldummy;
5074 AliTPCTrackPoint pdummy;
5075 pdummy.GetTPoint().fIsShared = 10;
5076 for (Int_t i=0;i<160;i++){
5077 AliTPCclusterMI * cl0 = fClusterPointer[i];
5078 AliTPCTrackPoint *trpoint = (AliTPCTrackPoint*)fPoints->UncheckedAt(i);
5080 trpoint->GetTPoint() = *(GetTrackPoint(i));
5081 trpoint->GetCPoint() = *cl0;
5082 trpoint->GetCPoint().SetQ(TMath::Abs(cl0->GetQ()));
5086 trpoint->GetCPoint()= cldummy;
5094 Double_t AliTPCseed::GetDensityFirst(Int_t n)
5098 // return cluster for n rows bellow first point
5099 Int_t nfoundable = 1;
5101 for (Int_t i=fLastPoint-1;i>0&&nfoundable<n; i--){
5102 Int_t index = GetClusterIndex2(i);
5103 if (index!=-1) nfoundable++;
5104 if (index>0) nfound++;
5106 if (nfoundable<n) return 0;
5107 return Double_t(nfound)/Double_t(nfoundable);
5112 void AliTPCseed::GetClusterStatistic(Int_t first, Int_t last, Int_t &found, Int_t &foundable, Int_t &shared, Bool_t plus2)
5114 // get cluster stat. on given region
5119 for (Int_t i=first;i<last; i++){
5120 Int_t index = GetClusterIndex2(i);
5121 if (index!=-1) foundable++;
5122 if (fClusterPointer[i]) {
5128 if (fClusterPointer[i]->IsUsed(10)) {
5132 if (!plus2) continue; //take also neighborhoud
5134 if ( (i>0) && fClusterPointer[i-1]){
5135 if (fClusterPointer[i-1]->IsUsed(10)) {
5140 if ( fClusterPointer[i+1]){
5141 if (fClusterPointer[i+1]->IsUsed(10)) {
5148 //if (shared>found){
5149 //Error("AliTPCseed::GetClusterStatistic","problem\n");
5153 //_____________________________________________________________________________
5154 void AliTPCseed::CookdEdx(Double_t low, Double_t up,Int_t i1, Int_t i2, Bool_t onlyused) {
5155 //-----------------------------------------------------------------
5156 // This funtion calculates dE/dX within the "low" and "up" cuts.
5157 //-----------------------------------------------------------------
5160 Float_t angular[200];
5161 Float_t weight[200];
5164 // TClonesArray & arr = *fPoints;
5165 Float_t meanlog = 100.;
5167 Float_t mean[4] = {0,0,0,0};
5168 Float_t sigma[4] = {1000,1000,1000,1000};
5169 Int_t nc[4] = {0,0,0,0};
5170 Float_t norm[4] = {1000,1000,1000,1000};
5175 for (Int_t of =0; of<4; of++){
5176 for (Int_t i=of+i1;i<i2;i+=4)
5178 Int_t index = fIndex[i];
5179 if (index<0||index&0x8000) continue;
5181 //AliTPCTrackPoint * point = (AliTPCTrackPoint *) arr.At(i);
5182 AliTPCTrackerPoint * point = GetTrackPoint(i);
5183 //AliTPCTrackerPoint * pointm = GetTrackPoint(i-1);
5184 //AliTPCTrackerPoint * pointp = 0;
5185 //if (i<159) pointp = GetTrackPoint(i+1);
5187 if (point==0) continue;
5188 AliTPCclusterMI * cl = fClusterPointer[i];
5189 if (cl==0) continue;
5190 if (onlyused && (!cl->IsUsed(10))) continue;
5191 if (cl->IsUsed(11)) {
5195 Int_t type = cl->GetType();
5196 //if (point->fIsShared){
5201 // if (pointm->fIsShared) continue;
5203 // if (pointp->fIsShared) continue;
5205 if (type<0) continue;
5206 //if (type>10) continue;
5207 //if (point->GetErrY()==0) continue;
5208 //if (point->GetErrZ()==0) continue;
5210 //Float_t ddy = (point->GetY()-cl->GetY())/point->GetErrY();
5211 //Float_t ddz = (point->GetZ()-cl->GetZ())/point->GetErrZ();
5212 //if ((ddy*ddy+ddz*ddz)>10) continue;
5215 // if (point->GetCPoint().GetMax()<5) continue;
5216 if (cl->GetMax()<5) continue;
5217 Float_t angley = point->GetAngleY();
5218 Float_t anglez = point->GetAngleZ();
5220 Float_t rsigmay2 = point->GetSigmaY();
5221 Float_t rsigmaz2 = point->GetSigmaZ();
5225 rsigmay += pointm->GetTPoint().GetSigmaY();
5226 rsigmaz += pointm->GetTPoint().GetSigmaZ();
5230 rsigmay += pointp->GetTPoint().GetSigmaY();
5231 rsigmaz += pointp->GetTPoint().GetSigmaZ();
5238 Float_t rsigma = TMath::Sqrt(rsigmay2*rsigmaz2);
5240 Float_t ampc = 0; // normalization to the number of electrons
5242 // ampc = 1.*point->GetCPoint().GetMax();
5243 ampc = 1.*cl->GetMax();
5244 //ampc = 1.*point->GetCPoint().GetQ();
5245 // AliTPCClusterPoint & p = point->GetCPoint();
5246 // Float_t dy = TMath::Abs(Int_t( TMath::Abs(p.GetY()/0.6)) - TMath::Abs(p.GetY()/0.6)+0.5);
5247 // Float_t iz = (250.0-TMath::Abs(p.GetZ())+0.11)/0.566;
5249 // TMath::Abs( Int_t(iz) - iz + 0.5);
5250 //ampc *= 1.15*(1-0.3*dy);
5251 //ampc *= 1.15*(1-0.3*dz);
5252 // Float_t zfactor = (1.05-0.0004*TMath::Abs(point->GetCPoint().GetZ()));
5256 //ampc = 1.0*point->GetCPoint().GetMax();
5257 ampc = 1.0*cl->GetMax();
5258 //ampc = 1.0*point->GetCPoint().GetQ();
5259 //AliTPCClusterPoint & p = point->GetCPoint();
5260 // Float_t dy = TMath::Abs(Int_t( TMath::Abs(p.GetY()/0.4)) - TMath::Abs(p.GetY()/0.4)+0.5);
5261 //Float_t iz = (250.0-TMath::Abs(p.GetZ())+0.11)/0.566;
5263 // TMath::Abs( Int_t(iz) - iz + 0.5);
5265 //ampc *= 1.15*(1-0.3*dy);
5266 //ampc *= 1.15*(1-0.3*dz);
5267 // Float_t zfactor = (1.02-0.000*TMath::Abs(point->GetCPoint().GetZ()));
5271 ampc *= 2.0; // put mean value to channel 50
5272 //ampc *= 0.58; // put mean value to channel 50
5274 // if (type>0) w = 1./(type/2.-0.5);
5275 // Float_t z = TMath::Abs(cl->GetZ());
5278 //ampc /= (1+0.0008*z);
5282 //ampc /= (1+0.0008*z);
5284 //ampc /= (1+0.0008*z);
5287 if (type<0) { //amp at the border - lower weight
5292 if (rsigma>1.5) ampc/=1.3; // if big backround
5294 angular[nc[of]] = TMath::Sqrt(1.+angley*angley+anglez*anglez);
5299 TMath::Sort(nc[of],amp,index,kFALSE);
5303 //meanlog = amp[index[Int_t(nc[of]*0.33)]];
5305 for (Int_t i=int(nc[of]*low+0.5);i<int(nc[of]*up+0.5);i++){
5306 Float_t ampl = amp[index[i]]/angular[index[i]];
5307 ampl = meanlog*TMath::Log(1.+ampl/meanlog);
5309 sumw += weight[index[i]];
5310 sumamp += weight[index[i]]*ampl;
5311 sumamp2 += weight[index[i]]*ampl*ampl;
5312 norm[of] += angular[index[i]]*weight[index[i]];
5319 mean[of] = sumamp/sumw;
5320 sigma[of] = sumamp2/sumw-mean[of]*mean[of];
5322 sigma[of] = TMath::Sqrt(sigma[of]);
5326 mean[of] = (TMath::Exp(mean[of]/meanlog)-1)*meanlog;
5327 //mean *=(1-0.02*(sigma/(mean*0.17)-1.));
5328 //mean *=(1-0.1*(norm-1.));
5335 // mean[0]*= (1-0.05*(sigma[0]/(0.01+mean[1]*0.18)-1));
5336 // mean[1]*= (1-0.05*(sigma[1]/(0.01+mean[0]*0.18)-1));
5339 // dedx = (mean[0]* TMath::Sqrt((1.+nc[0]))+ mean[1]* TMath::Sqrt((1.+nc[1])) )/
5340 // ( TMath::Sqrt((1.+nc[0]))+TMath::Sqrt((1.+nc[1])));
5344 for (Int_t i =0;i<4;i++){
5345 if (nc[i]>2&&nc[i]<1000){
5346 dedx += mean[i] *nc[i];
5347 fSdEdx += sigma[i]*(nc[i]-2);
5348 fMAngular += norm[i] *nc[i];
5353 fSDEDX[i] = sigma[i];
5366 // Float_t dedx1 =dedx;
5369 for (Int_t i =0;i<4;i++){
5370 if (nc[i]>2&&nc[i]<1000){
5371 mean[i] = mean[i]*(1-0.12*(sigma[i]/(fSdEdx)-1.));
5372 dedx += mean[i] *nc[i];
5387 Double_t p=TMath::Sqrt((1.+ GetTgl()*GetTgl())/(Get1Pt()*Get1Pt()));
5390 if (dedx < 39.+ 12./(p+0.25)/(p+0.25)) { SetMass(0.13957); return;}
5391 if (dedx < 39.+ 12./p/p) { SetMass(0.49368); return;}
5392 SetMass(0.93827); return;
5396 if (dedx < 39.+ 12./(p+0.25)/(p+0.25)) { SetMass(0.13957); return;}
5397 SetMass(0.93827); return;
5400 SetMass(0.13957); return;
5410 void AliTPCseed::CookdEdx2(Double_t low, Double_t up) {
5411 //-----------------------------------------------------------------
5412 // This funtion calculates dE/dX within the "low" and "up" cuts.
5413 //-----------------------------------------------------------------
5416 Float_t angular[200];
5417 Float_t weight[200];
5419 Bool_t inlimit[200];
5420 for (Int_t i=0;i<200;i++) inlimit[i]=kFALSE;
5421 for (Int_t i=0;i<200;i++) amp[i]=10000;
5422 for (Int_t i=0;i<200;i++) angular[i]= 1;;
5426 Float_t meanlog = 100.;
5427 Int_t indexde[4]={0,64,128,160};
5434 Float_t mean[4] = {0,0,0,0};
5435 Float_t sigma[4] = {1000,1000,1000,1000};
5436 Int_t nc[4] = {0,0,0,0};
5437 Float_t norm[4] = {1000,1000,1000,1000};
5442 // for (Int_t of =0; of<3; of++){
5443 // for (Int_t i=indexde[of];i<indexde[of+1];i++)
5444 for (Int_t i =0; i<160;i++)
5446 AliTPCTrackPoint * point = GetTrackPoint(i);
5447 if (point==0) continue;
5448 if (point->fIsShared){
5452 Int_t type = point->GetCPoint().GetType();
5453 if (type<0) continue;
5454 if (point->GetCPoint().GetMax()<5) continue;
5455 Float_t angley = point->GetTPoint().GetAngleY();
5456 Float_t anglez = point->GetTPoint().GetAngleZ();
5457 Float_t rsigmay = point->GetCPoint().GetSigmaY();
5458 Float_t rsigmaz = point->GetCPoint().GetSigmaZ();
5459 Float_t rsigma = TMath::Sqrt(rsigmay*rsigmaz);
5461 Float_t ampc = 0; // normalization to the number of electrons
5463 ampc = point->GetCPoint().GetMax();
5466 ampc = point->GetCPoint().GetMax();
5468 ampc *= 2.0; // put mean value to channel 50
5469 // ampc *= 0.565; // put mean value to channel 50
5472 Float_t z = TMath::Abs(point->GetCPoint().GetZ());
5479 if (type<0) { //amp at the border - lower weight
5482 if (rsigma>1.5) ampc/=1.3; // if big backround
5483 angular[i] = TMath::Sqrt(1.+angley*angley+anglez*anglez);
5484 amp[i] = ampc/angular[i];
5489 TMath::Sort(159,amp,index,kFALSE);
5490 for (Int_t i=int(anc*low+0.5);i<int(anc*up+0.5);i++){
5491 inlimit[index[i]] = kTRUE; // take all clusters
5494 // meanlog = amp[index[Int_t(anc*0.3)]];
5496 for (Int_t of =0; of<3; of++){
5500 for (Int_t i=indexde[of];i<indexde[of+1];i++)
5502 if (inlimit[i]==kFALSE) continue;
5503 Float_t ampl = amp[i];
5505 ampl = meanlog*TMath::Log(1.+ampl/meanlog);
5508 sumamp += weight[i]*ampl;
5509 sumamp2 += weight[i]*ampl*ampl;
5510 norm[of] += angular[i]*weight[i];
5518 mean[of] = sumamp/sumw;
5519 sigma[of] = sumamp2/sumw-mean[of]*mean[of];
5521 sigma[of] = TMath::Sqrt(sigma[of]);
5524 mean[of] = (TMath::Exp(mean[of]/meanlog)-1)*meanlog;
5534 Float_t www[3] = {12.,14.,17.};
5535 //Float_t www[3] = {1.,1.,1.};
5537 for (Int_t i =0;i<3;i++){
5538 if (nc[i]>2&&nc[i]<1000){
5539 dedx += mean[i] *nc[i]*www[i]/sigma[i];
5540 fSdEdx += sigma[i]*(nc[i]-2)*www[i]/sigma[i];
5541 fMAngular += norm[i] *nc[i];
5542 norm2 += nc[i]*www[i]/sigma[i];
5543 norm3 += (nc[i]-2)*www[i]/sigma[i];
5546 fSDEDX[i] = sigma[i];
5559 // Float_t dedx1 =dedx;
5563 for (Int_t i =0;i<3;i++){
5564 if (nc[i]>2&&nc[i]<1000&&sigma[i]>3){
5565 //mean[i] = mean[i]*(1+0.08*(sigma[i]/(fSdEdx)-1.));
5566 dedx += mean[i] *(nc[i])/(sigma[i]);
5567 norm4 += (nc[i])/(sigma[i]);
5571 if (norm4>0) dedx /= norm4;