Bug fix (M.Ivanov)
[u/mrichter/AliRoot.git] / TPC / AliTPCtrackerMI.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
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  **************************************************************************/
15
16
17
18
19
20
21 /*
22   AliTPC parallel tracker - 
23   How to use?  - 
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 */
29
30 //-------------------------------------------------------
31 //          Implementation of the TPC tracker
32 //
33 //   Origin: Marian Ivanov   Marian.Ivanov@cern.ch
34 // 
35 //-------------------------------------------------------
36
37 #include <TObjArray.h>
38 #include <TFile.h>
39 #include <TTree.h>
40 #include "Riostream.h"
41
42 #include "AliTPCtrackerMI.h"
43 #include "AliTPCclusterMI.h"
44 #include "AliTPCParam.h"
45 #include "AliTPCClustersRow.h"
46 #include "AliComplexCluster.h"
47 #include "AliTPCpolyTrack.h"
48 #include "TStopwatch.h"
49
50
51 ClassImp(AliTPCseed)
52
53 AliTPCclusterTracks::AliTPCclusterTracks(){
54   // class for storing overlaping info
55   fTrackIndex[0]=-1;
56   fTrackIndex[1]=-1;
57   fTrackIndex[2]=-1;
58   fDistance[0]=1000;
59   fDistance[1]=1000;
60   fDistance[2]=1000;
61 }
62
63
64
65
66
67 Int_t AliTPCtrackerMI::UpdateTrack(AliTPCseed * track, AliTPCclusterMI* c, Double_t chi2, UInt_t i){
68
69   Int_t sec=(i&0xff000000)>>24; 
70   Int_t row = (i&0x00ff0000)>>16;
71   track->fRow=(i&0x00ff0000)>>16;
72   track->fSector = sec;
73   //  Int_t index = i&0xFFFF;
74   if (sec>=fParam->GetNInnerSector()) track->fRow += fParam->GetNRowLow(); 
75   track->fClusterIndex[track->fRow] = i;
76   track->fFirstPoint = row;
77   //
78
79   AliTPCTrackPoint   *trpoint =track->GetTrackPoint(track->fRow);
80   Float_t angle2 = track->GetSnp()*track->GetSnp();
81   angle2 = TMath::Sqrt(angle2/(1-angle2)); 
82   //
83   //SET NEW Track Point
84   //
85   if (c!=0){
86     //if we have a cluster
87     trpoint->GetCPoint().SetY(c->GetY());
88     trpoint->GetCPoint().SetZ(c->GetZ());    
89     //
90     trpoint->GetCPoint().SetSigmaY(c->GetSigmaY2()/(track->fCurrentSigmaY*track->fCurrentSigmaY));
91     trpoint->GetCPoint().SetSigmaZ(c->GetSigmaZ2()/(track->fCurrentSigmaZ*track->fCurrentSigmaZ));
92     //
93     trpoint->GetCPoint().SetType(c->GetType());
94     trpoint->GetCPoint().SetQ(c->GetQ());
95     trpoint->GetCPoint().SetMax(c->GetMax());
96     //  
97     trpoint->GetCPoint().SetErrY(TMath::Sqrt(track->fErrorY2));
98     trpoint->GetCPoint().SetErrZ(TMath::Sqrt(track->fErrorZ2));
99     //
100   }
101   trpoint->GetTPoint().SetX(track->GetX());
102   trpoint->GetTPoint().SetY(track->GetY());
103   trpoint->GetTPoint().SetZ(track->GetZ());
104   //
105   trpoint->GetTPoint().SetAngleY(angle2);
106   trpoint->GetTPoint().SetAngleZ(track->GetTgl());
107   
108
109   
110   if (chi2>10){
111     //    printf("suspicious chi2 %f\n",chi2);
112   }
113   //  if (track->fIsSeeding){ 
114     track->fErrorY2 *= 1.2;
115     track->fErrorY2 += 0.0064;    
116     track->fErrorZ2 *= 1.2;   
117     track->fErrorY2 += 0.005;    
118  
119     //}
120
121   return track->Update(c,chi2,i);
122
123 }
124 //_____________________________________________________________________________
125 AliTPCtrackerMI::AliTPCtrackerMI(const AliTPCParam *par, Int_t eventn): 
126 AliTracker(), fkNIS(par->GetNInnerSector()/2), fkNOS(par->GetNOuterSector()/2)
127 {
128   //---------------------------------------------------------------------
129   // The main TPC tracker constructor
130   //---------------------------------------------------------------------
131   fInnerSec=new AliTPCSector[fkNIS];         
132   fOuterSec=new AliTPCSector[fkNOS];
133
134   Int_t i;
135   for (i=0; i<fkNIS; i++) fInnerSec[i].Setup(par,0);
136   for (i=0; i<fkNOS; i++) fOuterSec[i].Setup(par,1);
137
138   fN=0;  fSectors=0;
139
140   fClustersArray.Setup(par);
141   fClustersArray.SetClusterType("AliTPCclusterMI");
142
143   char   cname[100];
144   if (eventn==-1) {
145     sprintf(cname,"TreeC_TPC");
146   }
147   else {
148     sprintf(cname,"TreeC_TPC_%d",eventn);
149   }
150
151   fClustersArray.ConnectTree(cname);
152
153   fEventN = eventn;
154   fSeeds=0;
155   fNtracks = 0;
156   fParam = par;
157 }
158
159 //_____________________________________________________________________________
160 AliTPCtrackerMI::~AliTPCtrackerMI() {
161   //------------------------------------------------------------------
162   // TPC tracker destructor
163   //------------------------------------------------------------------
164   delete[] fInnerSec;
165   delete[] fOuterSec;
166   if (fSeeds) {
167     fSeeds->Delete(); 
168     delete fSeeds;
169   }
170 }
171
172
173 Double_t AliTPCtrackerMI::ErrY2(AliTPCseed* seed, AliTPCclusterMI * cl){
174   //
175   //
176   Float_t snoise2;
177   Float_t z = fParam->GetZLength()-TMath::Abs(seed->GetZ());
178
179   //cluster "quality"
180   Float_t rsigmay = 1;
181   Int_t ctype = 0;
182
183   //standard if we don't have cluster - take MIP
184   const Float_t chmip = 50.; 
185   Float_t amp = chmip/0.3;  
186   Float_t nel;
187   Float_t nprim;
188   if (cl){
189     amp = cl->GetQ();
190     rsigmay = cl->GetSigmaY2()/(seed->fCurrentSigmaY*seed->fCurrentSigmaY);
191     ctype = cl->GetType();
192   }
193   
194
195   Float_t landau=2 ;    //landau fluctuation part
196   Float_t gg=2;         // gg fluctuation part
197   Float_t padlength= fSectors->GetPadPitchLength(seed->GetX());
198
199
200   if (fSectors==fInnerSec){
201     snoise2 = 0.0004;
202     nel     = 0.268*amp;
203     nprim   = 0.155*amp;
204     gg      = (2+0.0002*amp)/nel;
205     landau  = (2.+0.12*nprim)*0.5*(amp*amp/40000.+2)/nprim;
206     if (landau>1) landau=1;
207   }
208   else {
209     snoise2 = 0.0004;
210     nel     = 0.3*amp;
211     nprim   = 0.133*amp;
212     gg      = (2+0.0002*amp)/nel;
213     landau  = (2.+0.12*nprim)*0.5*(amp*amp/40000.+2)/nprim;
214     if (landau>1) landau=1;
215   }
216
217
218   Float_t sdiff = gg*fParam->GetDiffT()*fParam->GetDiffT()*z;
219   Float_t angle2 = seed->GetSnp()*seed->GetSnp();
220   angle2 = angle2/(1-angle2); 
221   Float_t angular = landau*angle2*padlength*padlength/12.;
222   Float_t res = sdiff + angular;
223   
224   
225   if ((ctype==0) && (fSectors ==fOuterSec))
226     res *= 0.78 +TMath::Exp(7.4*(rsigmay-1.2));
227
228   if ((ctype==0) && (fSectors ==fInnerSec))
229     res *= 0.72 +TMath::Exp(3.36*(rsigmay-1.2));
230
231
232   if ((ctype>0))
233     res*= TMath::Power((rsigmay+0.5),1.5)+0.0064;
234   
235   if (ctype<0)
236     res*=2.4;  // overestimate error 2 times
237   
238   res+= snoise2;
239  
240   if (res<2*snoise2)
241     res = 2*snoise2;
242
243   seed->SetErrorY2(res);
244   return res;
245 }
246
247
248
249
250 Double_t AliTPCtrackerMI::ErrZ2(AliTPCseed* seed, AliTPCclusterMI * cl){
251   //
252   //
253   Float_t snoise2;
254   Float_t z = fParam->GetZLength()-TMath::Abs(seed->GetZ());
255   //signal quality
256   Float_t rsigmaz=1;
257   Int_t ctype =0;
258
259   const Float_t chmip = 50.;
260   Float_t amp = chmip/0.3;  
261   Float_t nel;
262   Float_t nprim;
263   if (cl){
264     amp = cl->GetQ();
265     rsigmaz = cl->GetSigmaZ2()/(seed->fCurrentSigmaZ*seed->fCurrentSigmaZ);
266     ctype = cl->GetType();
267   }
268
269   //
270   Float_t landau=2 ;    //landau fluctuation part
271   Float_t gg=2;         // gg fluctuation part
272   Float_t padlength= fSectors->GetPadPitchLength(seed->GetX());
273
274   if (fSectors==fInnerSec){
275     snoise2 = 0.0004;
276     nel     = 0.268*amp;
277     nprim   = 0.155*amp;
278     gg      = (2+0.0002*amp)/nel;
279     landau  = (2.+0.12*nprim)*0.5*(amp*amp/40000.+2)/nprim;
280     if (landau>1) landau=1;
281   }
282   else {
283     snoise2 = 0.0004;
284     nel     = 0.3*amp;
285     nprim   = 0.133*amp;
286     gg      = (2+0.0002*amp)/nel;
287     landau  = (2.+0.12*nprim)*0.5*(amp*amp/40000.+2)/nprim;
288     if (landau>1) landau=1;
289   }
290   Float_t sdiff = gg*fParam->GetDiffT()*fParam->GetDiffT()*z;
291   
292   Float_t angle = seed->GetTgl();
293   Float_t angular = landau*angle*angle*padlength*padlength/12.;
294   Float_t res = sdiff + angular;
295
296   if ((ctype==0) && (fSectors ==fOuterSec))
297     res *= 0.81 +TMath::Exp(6.8*(rsigmaz-1.2));
298
299   if ((ctype==0) && (fSectors ==fInnerSec))
300     res *= 0.72 +TMath::Exp(2.04*(rsigmaz-1.2));
301   if ((ctype>0))
302     res*= TMath::Power(rsigmaz+0.5,1.5)+0.0064;  //0.31+0.147*ctype;
303   if (ctype<0)
304     res*=1.3;
305   if ((ctype<0) &&amp<70)
306     res*=1.3;  
307
308   res += snoise2;
309   if (res<2*snoise2)
310      res = 2*snoise2;
311
312   seed->SetErrorZ2(res);
313   return res;
314 }
315
316
317 void AliTPCseed::Reset()
318 {
319   //
320   //PH  SetN(0);
321   fNFoundable = 0;
322   ResetCovariance();
323   SetChi2(0);
324   for (Int_t i=0;i<200;i++) fClusterIndex[i]=-1;
325 }
326
327
328 Int_t  AliTPCseed::GetProlongation(Double_t xk, Double_t &y, Double_t & z) const
329 {
330   //-----------------------------------------------------------------
331   // This function find proloncation of a track to a reference plane x=xk.
332   // doesn't change internal state of the track
333   //-----------------------------------------------------------------
334   
335   Double_t x1=fX, x2=x1+(xk-x1), dx=x2-x1;
336   //  Double_t y1=fP0, z1=fP1;
337   Double_t c1=fP4*x1 - fP2, r1=sqrt(1.- c1*c1);
338   Double_t c2=fP4*x2 - fP2, r2=sqrt(1.- c2*c2);
339   
340   y = fP0;
341   z = fP1;
342   y += dx*(c1+c2)/(r1+r2);
343   z += dx*(c1+c2)/(c1*r2 + c2*r1)*fP3;
344   return 0;  
345 }
346
347
348 //_____________________________________________________________________________
349 Double_t AliTPCseed::GetPredictedChi2(const AliTPCclusterMI *c) const 
350 {
351   //-----------------------------------------------------------------
352   // This function calculates a predicted chi2 increment.
353   //-----------------------------------------------------------------
354   //Double_t r00=c->GetSigmaY2(), r01=0., r11=c->GetSigmaZ2();
355   Double_t r00=fErrorY2, r01=0., r11=fErrorZ2;
356   r00+=fC00; r01+=fC10; r11+=fC11;
357
358   Double_t det=r00*r11 - r01*r01;
359   if (TMath::Abs(det) < 1.e-10) {
360     Int_t n=GetNumberOfClusters();
361     if (n>4) cerr<<n<<" AliKalmanTrack warning: Singular matrix !\n";
362     return 1e10;
363   }
364   Double_t tmp=r00; r00=r11; r11=tmp; r01=-r01;
365   
366   Double_t dy=c->GetY() - fP0, dz=c->GetZ() - fP1;
367   
368   return (dy*r00*dy + 2*r01*dy*dz + dz*r11*dz)/det;
369 }
370
371
372 //_________________________________________________________________________________________
373
374
375 Int_t AliTPCseed::Compare(const TObject *o) const {
376   //-----------------------------------------------------------------
377   // This function compares tracks according to the sector - for given sector according z
378   //-----------------------------------------------------------------
379   AliTPCseed *t=(AliTPCseed*)o;
380   if (t->fSector>fSector) return -1;
381   if (t->fSector<fSector) return 1;
382
383   Double_t z2 = t->GetZ();
384   Double_t z1 = GetZ();
385   if (z2>z1) return 1;
386   if (z2<z1) return -1;
387   return 0;
388 }
389
390
391
392
393
394
395 //_____________________________________________________________________________
396 Int_t AliTPCseed::Update(const AliTPCclusterMI *c, Double_t chisq, UInt_t index) {
397   //-----------------------------------------------------------------
398   // This function associates a cluster with this track.
399   //-----------------------------------------------------------------
400   //  Double_t r00=c->GetSigmaY2(), r01=0., r11=c->GetSigmaZ2();
401   //Double_t r00=sigmay2, r01=0., r11=sigmaz2;
402   Double_t r00=fErrorY2, r01=0., r11=fErrorZ2;
403
404   r00+=fC00; r01+=fC10; r11+=fC11;
405   Double_t det=r00*r11 - r01*r01;
406   Double_t tmp=r00; r00=r11/det; r11=tmp/det; r01=-r01/det;
407
408   Double_t k00=fC00*r00+fC10*r01, k01=fC00*r01+fC10*r11;
409   Double_t k10=fC10*r00+fC11*r01, k11=fC10*r01+fC11*r11;
410   Double_t k20=fC20*r00+fC21*r01, k21=fC20*r01+fC21*r11;
411   Double_t k30=fC30*r00+fC31*r01, k31=fC30*r01+fC31*r11;
412   Double_t k40=fC40*r00+fC41*r01, k41=fC40*r01+fC41*r11;
413
414   Double_t dy=c->GetY() - fP0, dz=c->GetZ() - fP1;
415   Double_t cur=fP4 + k40*dy + k41*dz, eta=fP2 + k20*dy + k21*dz;
416   if (TMath::Abs(cur*fX-eta) >= 0.9) {
417     //    Int_t n=GetNumberOfClusters();
418     //if (n>4) cerr<<n<<" AliTPCtrack warning: Filtering failed !\n";
419     return 0;
420   }
421
422   fP0 += k00*dy + k01*dz;
423   fP1 += k10*dy + k11*dz;
424   fP2  = eta;
425   fP3 += k30*dy + k31*dz;
426   fP4  = cur;
427
428   Double_t c01=fC10, c02=fC20, c03=fC30, c04=fC40;
429   Double_t c12=fC21, c13=fC31, c14=fC41;
430
431   fC00-=k00*fC00+k01*fC10; fC10-=k00*c01+k01*fC11;
432   fC20-=k00*c02+k01*c12;   fC30-=k00*c03+k01*c13;
433   fC40-=k00*c04+k01*c14; 
434
435   fC11-=k10*c01+k11*fC11;
436   fC21-=k10*c02+k11*c12;   fC31-=k10*c03+k11*c13;
437   fC41-=k10*c04+k11*c14; 
438
439   fC22-=k20*c02+k21*c12;   fC32-=k20*c03+k21*c13;
440   fC42-=k20*c04+k21*c14; 
441
442   fC33-=k30*c03+k31*c13;
443   fC43-=k40*c03+k41*c13; 
444
445   fC44-=k40*c04+k41*c14; 
446
447   Int_t n=GetNumberOfClusters();
448   fIndex[n]=index;
449   SetNumberOfClusters(n+1);
450   SetChi2(GetChi2()+chisq);
451
452   return 1;
453 }
454
455
456
457 //_____________________________________________________________________________
458 Double_t AliTPCtrackerMI::f1(Double_t x1,Double_t y1,
459                    Double_t x2,Double_t y2,
460                    Double_t x3,Double_t y3) 
461 {
462   //-----------------------------------------------------------------
463   // Initial approximation of the track curvature
464   //-----------------------------------------------------------------
465   Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
466   Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
467                   (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
468   Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
469                   (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
470
471   Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
472
473   return -xr*yr/sqrt(xr*xr+yr*yr); 
474 }
475
476
477 //_____________________________________________________________________________
478 Double_t AliTPCtrackerMI::f2(Double_t x1,Double_t y1,
479                    Double_t x2,Double_t y2,
480                    Double_t x3,Double_t y3) 
481 {
482   //-----------------------------------------------------------------
483   // Initial approximation of the track curvature times center of curvature
484   //-----------------------------------------------------------------
485   Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
486   Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
487                   (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
488   Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
489                   (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
490
491   Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
492   
493   return -a/(d*y1-b)*xr/sqrt(xr*xr+yr*yr);
494 }
495
496 //_____________________________________________________________________________
497 Double_t AliTPCtrackerMI::f3(Double_t x1,Double_t y1, 
498                    Double_t x2,Double_t y2,
499                    Double_t z1,Double_t z2) 
500 {
501   //-----------------------------------------------------------------
502   // Initial approximation of the tangent of the track dip angle
503   //-----------------------------------------------------------------
504   return (z1 - z2)/sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
505 }
506
507
508 void AliTPCtrackerMI::LoadClusters()
509 {
510   //
511   // load clusters to the memory
512   Int_t j=Int_t(fClustersArray.GetTree()->GetEntries());
513   for (Int_t i=0; i<j; i++) {
514     fClustersArray.LoadEntry(i);
515   }
516 }
517
518 void AliTPCtrackerMI::UnloadClusters()
519 {
520   //
521   // load clusters to the memory
522   Int_t j=Int_t(fClustersArray.GetTree()->GetEntries());
523   for (Int_t i=0; i<j; i++) {
524     fClustersArray.ClearSegment(i);
525   }
526 }
527
528
529
530 //_____________________________________________________________________________
531 void AliTPCtrackerMI::LoadOuterSectors() {
532   //-----------------------------------------------------------------
533   // This function fills outer TPC sectors with clusters.
534   //-----------------------------------------------------------------
535   UInt_t index;
536   //Int_t j=Int_t(fClustersArray.GetTree()->GetEntries());
537   Int_t j = ((AliTPCParam*)fParam)->GetNRowsTotal();
538   for (Int_t i=0; i<j; i++) {
539     //  AliSegmentID *s=fClustersArray.LoadEntry(i);
540     AliSegmentID *s= const_cast<AliSegmentID*>(fClustersArray.At(i));
541     if (!s) continue;
542     Int_t sec,row;
543     AliTPCParam *par=(AliTPCParam*)fClustersArray.GetParam();
544     par->AdjustSectorRow(s->GetID(),sec,row);
545     if (sec<fkNIS*2) continue;
546     AliTPCClustersRow *clrow=fClustersArray.GetRow(sec,row);
547     Int_t ncl=clrow->GetArray()->GetEntriesFast();
548     while (ncl--) {
549       AliTPCclusterMI *c=(AliTPCclusterMI*)(*clrow)[ncl];
550       index=(((sec<<8)+row)<<16)+ncl;
551       fOuterSec[(sec-fkNIS*2)%fkNOS][row].InsertCluster(c,index);
552     }
553   }  
554   fN=fkNOS;
555   fSectors=fOuterSec;
556 }
557
558
559 //_____________________________________________________________________________
560 void AliTPCtrackerMI::LoadInnerSectors() {
561   //-----------------------------------------------------------------
562   // This function fills inner TPC sectors with clusters.
563   //-----------------------------------------------------------------
564   UInt_t index;
565   //Int_t j=Int_t(fClustersArray.GetTree()->GetEntries());
566   Int_t j = ((AliTPCParam*)fParam)->GetNRowsTotal();
567   for (Int_t i=0; i<j; i++) {
568     //   AliSegmentID *s=fClustersArray.LoadEntry(i);
569     AliSegmentID *s=const_cast<AliSegmentID*>(fClustersArray.At(i));
570     if (!s) continue;
571     Int_t sec,row;
572     AliTPCParam *par=(AliTPCParam*)fClustersArray.GetParam();
573     par->AdjustSectorRow(s->GetID(),sec,row);
574     if (sec>=fkNIS*2) continue;
575     AliTPCClustersRow *clrow=fClustersArray.GetRow(sec,row);
576     Int_t ncl=clrow->GetArray()->GetEntriesFast();
577     while (ncl--) {
578       AliTPCclusterMI *c=(AliTPCclusterMI*)(*clrow)[ncl];
579       index=(((sec<<8)+row)<<16)+ncl;
580       fInnerSec[sec%fkNIS][row].InsertCluster(c,index);
581     }
582   }
583
584   fN=fkNIS;
585   fSectors=fInnerSec;
586 }
587
588 Int_t AliTPCtrackerMI::FollowToNext(AliTPCseed& t, Int_t nr) {
589   //-----------------------------------------------------------------
590   // This function tries to find a track prolongation to next pad row
591   //-----------------------------------------------------------------
592   //  Double_t xt=t.GetX();
593   //  Int_t row = fSectors->GetRowNumber(xt)-1;
594   //  if (row < nr) return 1; // don't prolongate if not information until now -
595   //
596   Double_t  x=fSectors->GetX(nr), ymax=fSectors->GetMaxY(nr);
597   //  if (t.GetRadius()>x+10 ) return 0;
598
599   if (!t.PropagateTo(x)) {
600     t.fStopped = kTRUE;
601     return 0;
602   }
603   // update current
604   t.fCurrentSigmaY = GetSigmaY(&t);
605   t.fCurrentSigmaZ = GetSigmaZ(&t);
606   //  
607   AliTPCclusterMI *cl=0;
608   UInt_t index=0;
609   const AliTPCRow &krow=fSectors[t.fRelativeSector][nr];
610   Double_t sy2=ErrY2(&t)*2;
611   Double_t sz2=ErrZ2(&t)*2;
612
613
614   Double_t  roady  =3.*sqrt(t.GetSigmaY2() + sy2);
615   Double_t  roadz = 3 *sqrt(t.GetSigmaZ2() + sz2);
616   Double_t  y=t.GetY(), z=t.GetZ();
617
618   if (TMath::Abs(TMath::Abs(y)-ymax)<krow.fDeadZone){
619     t.fInDead = kTRUE;
620     return 0;
621   } 
622   else
623     {
624       if (TMath::Abs(z)<(1.05*x+10)) t.fNFoundable++;
625       else
626         return 0;
627     }   
628   //calculate 
629   Float_t maxdistance = roady*roady + roadz*roadz;
630   if (krow) {
631     for (Int_t i=krow.Find(z-roadz); i<krow; i++) {
632       AliTPCclusterMI *c=(AliTPCclusterMI*)(krow[i]);
633       if (c->GetZ() > z+roadz) break;
634       if ( (c->GetY()-y) >  roady ) continue;
635       Float_t distance = (c->GetZ()-z)*(c->GetZ()-z)+(c->GetY()-y)*(c->GetY()-y);
636       if (maxdistance>distance) {
637         maxdistance = distance;
638         cl=c;
639         index=krow.GetIndex(i);       
640       }
641     }
642   }      
643   if (cl) {
644     //Double_t sy2=
645     ErrY2(&t,cl);
646     //Double_t sz2=
647     ErrZ2(&t,cl);
648     Double_t chi2= t.GetPredictedChi2(cl);    
649     UpdateTrack(&t,cl,chi2,index);   
650   } else {    
651     if (y > ymax) {
652       t.fRelativeSector= (t.fRelativeSector+1) % fN;
653       if (!t.Rotate(fSectors->GetAlpha())) 
654         return 0;
655     } else if (y <-ymax) {
656       t.fRelativeSector= (t.fRelativeSector-1+fN) % fN;
657       if (!t.Rotate(-fSectors->GetAlpha())) 
658         return 0;
659     }   
660   }
661   return 1;
662 }
663
664
665 Int_t AliTPCtrackerMI::UpdateClusters(AliTPCseed& t,Int_t trindex,  Int_t nr) {
666   //-----------------------------------------------------------------
667   // This function tries to find a track prolongation to next pad row
668   //-----------------------------------------------------------------
669   t.fCurrentCluster  = 0;
670   t.fCurrentClusterIndex1 = 0;   
671   t.fCurrentClusterIndex2 = 0;
672    
673   Double_t xt=t.GetX();
674   Int_t row = fSectors->GetRowNumber(xt)-1;
675   if (row < nr) return 1; // don't prolongate if not information until now -
676   Double_t x=fSectors->GetX(nr);
677   //  if (t.fStopped) return 0;
678   //  if (t.GetRadius()>x+10 ) return 0;
679   if (!t.PropagateTo(x)){
680     t.fStopped =kTRUE;
681     return 0;
682   }
683   // update current
684   t.fCurrentSigmaY = GetSigmaY(&t);
685   t.fCurrentSigmaZ = GetSigmaZ(&t);
686     
687   AliTPCclusterMI *cl=0;
688   UInt_t index=0;
689   AliTPCRow &krow=fSectors[t.fRelativeSector][nr];
690   //
691   Double_t  y=t.GetY(), z=t.GetZ();
692   Double_t roady = 3.* TMath::Sqrt(t.GetSigmaY2() + t.fCurrentSigmaY*t.fCurrentSigmaY);
693   Double_t roadz = 3.* TMath::Sqrt(t.GetSigmaZ2() + t.fCurrentSigmaZ*t.fCurrentSigmaZ);
694   //
695
696   Float_t maxdistance = 1000000;
697   if (krow) {    
698     for (Int_t i=krow.Find(z-roadz); i<krow; i++) {
699       AliTPCclusterMI *c=(AliTPCclusterMI*)(krow[i]);
700       if (c->GetZ() > z+roadz) break;
701       if (TMath::Abs(c->GetY()-y)>roady) continue;
702             
703       //krow.UpdateClusterTrack(i,trindex,&t);
704
705       Float_t dy2 = (c->GetY()- t.GetY());
706       dy2*=dy2;
707       Float_t dz2 = (c->GetZ()- t.GetZ());
708       dz2*=dz2;
709       //
710       Float_t distance = dy2+dz2;
711       //
712       if (distance > maxdistance) continue;
713       maxdistance = distance;
714       cl=c;
715       index=i;       
716     }
717   }
718   t.fCurrentCluster  = cl;
719   t.fCurrentClusterIndex1 = krow.GetIndex(index);   
720   t.fCurrentClusterIndex2 = index;   
721   return 1;
722 }
723
724
725 Int_t AliTPCtrackerMI::FollowToNextCluster(Int_t trindex, Int_t nr) {
726   //-----------------------------------------------------------------
727   // This function tries to find a track prolongation to next pad row
728   //-----------------------------------------------------------------
729   AliTPCseed & t  = *((AliTPCseed*)(fSeeds->At(trindex)));
730   AliTPCRow &krow=fSectors[t.fRelativeSector][nr];
731   //  Double_t pt=t.GetConvConst()/(100/0.299792458/0.2)/t.Get1Pt();
732   Double_t y=t.GetY();
733   Double_t ymax=fSectors->GetMaxY(nr);
734
735   if (TMath::Abs(TMath::Abs(y)-ymax)<krow.fDeadZone){
736     t.fInDead = kTRUE;
737     return 0;
738   } 
739   else
740     {
741       if (TMath::Abs(t.GetZ())<(1.05*t.GetX()+10)) t.fNFoundable++;
742       else
743         return 0;      
744     }
745   
746   if (t.fCurrentCluster) {
747     //    Float_t l=fSectors->GetPadPitchLength();
748     //    AliTPCclusterTracks * cltrack = krow.GetClusterTracks(t.fCurrentClusterIndex1);
749
750     Double_t sy2=ErrY2(&t,t.fCurrentCluster);
751     Double_t sz2=ErrZ2(&t,t.fCurrentCluster);
752
753
754     Double_t sdistancey = TMath::Sqrt(sy2+t.GetSigmaY2());
755     Double_t sdistancez = TMath::Sqrt(sz2+t.GetSigmaZ2());
756
757     Double_t rdistancey = TMath::Abs(t.fCurrentCluster->GetY()-t.GetY());
758     Double_t rdistancez = TMath::Abs(t.fCurrentCluster->GetZ()-t.GetZ());
759     
760     Double_t rdistance  = TMath::Sqrt(TMath::Power(rdistancey/sdistancey,2)+TMath::Power(rdistancez/sdistancez,2));
761
762
763     //    printf("\t%f\t%f\t%f\n",rdistancey/sdistancey,rdistancez/sdistancez,rdistance);
764     if ( (rdistancey>1) || (rdistancez>1)) return 0;
765     if (rdistance>4) return 0;
766
767     if ((rdistancey/sdistancey>2.5 || rdistancez/sdistancez>2.5) && t.fCurrentCluster->GetType()==0)  
768         return 0;  //suspisiouce - will be changed
769
770     if ((rdistancey/sdistancey>2. || rdistancez/sdistancez>2.0) && t.fCurrentCluster->GetType()>0)  
771         // strict cut on overlaped cluster
772         return 0;  //suspisiouce - will be changed
773
774     if ( (rdistancey/sdistancey>1. || rdistancez/sdistancez>2.5 ||t.fCurrentCluster->GetQ()<70 ) 
775          && t.fCurrentCluster->GetType()<0)
776       return 0;
777
778     //    t.SetSampledEdx(0.3*t.fCurrentCluster->GetQ()/l,t.GetNumberOfClusters(), GetSigmaY(&t), GetSigmaZ(&t));
779     UpdateTrack(&t,t.fCurrentCluster,t.GetPredictedChi2(t.fCurrentCluster),t.fCurrentClusterIndex1);
780    
781   } else {
782     if (y > ymax) {
783       t.fRelativeSector= (t.fRelativeSector+1) % fN;
784       if (!t.Rotate(fSectors->GetAlpha())) 
785         return 0;
786     } else if (y <-ymax) {
787       t.fRelativeSector= (t.fRelativeSector-1+fN) % fN;
788       if (!t.Rotate(-fSectors->GetAlpha())) 
789         return 0;
790     }
791   }
792   return 1;
793 }
794
795
796
797
798 /*
799 Int_t AliTPCtrackerMI::FollowProlongationFast(AliTPCseed& t, Int_t step)
800 {
801   //-----------------------------------------------------------------
802   // fast prolongation mathod -
803   // don't update track only after step clusters
804   //-----------------------------------------------------------------
805   Double_t xt=t.GetX();
806   //
807   Double_t alpha=t.GetAlpha(); 
808   alpha =- fSectors->GetAlphaShift();
809   if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();  
810   if (alpha < 0.            ) alpha += 2.*TMath::Pi();  
811   t.fRelativeSector = Int_t(alpha/fSectors->GetAlpha())%fN;
812   Int_t row0 = fSectors->GetRowNumber(xt); 
813   Double_t x    = fSectors->GetX(row0);
814   Double_t ymax = fSectors->GetMaxY(row0);
815   //
816   Double_t sy2=ErrY2(&t)*2;
817   Double_t sz2=ErrZ2(&t)*2;
818   Double_t  roady  =3.*sqrt(t.GetSigmaY2() + sy2);
819   Double_t  roadz = 3 *sqrt(t.GetSigmaZ2() + sz2);
820   Float_t maxdistance = roady*roady + roadz*roadz; 
821   t.fCurrentSigmaY = GetSigmaY(&t);
822   t.fCurrentSigmaZ = GetSigmaZ(&t);
823   //
824   Int_t nclusters = 0;
825   Double_t y;
826   Double_t z;
827   Double_t yy[200];   //track prolongation
828   Double_t zz[200];
829   Double_t cy[200];  // founded cluster position
830   Double_t cz[200];
831   Double_t sy[200];  // founded cluster error
832   Double_t sz[200];
833   Bool_t   hitted[200]; // indication of cluster presence
834   //
835   
836   //
837   for (Int_t drow = step; drow>=0; drow--) {
838     Int_t row = row0-drow;
839     if (row<0) break;
840     Double_t x    = fSectors->GetX(row);
841     Double_t ymax = fSectors->GetMaxY(row);
842     t.GetProlongation(x,y,z);
843     yy[drow] =y;
844     zz[drow] =z;    
845     const AliTPCRow &krow=fSectors[t.fRelativeSector][row];
846     if (TMath::Abs(TMath::Abs(y)-ymax)<krow.fDeadZone){
847       t.fInDead = kTRUE;
848       break;
849     } 
850     else
851       {
852         t.fNFoundable++;
853       }  
854     
855     //find nearest  cluster 
856     AliTPCclusterMI *cl= 0;
857     if (krow) {
858       for (Int_t i=krow.Find(z-roadz); i<krow; i++) {
859         AliTPCclusterMI *c=(AliTPCclusterMI*)(krow[i]);
860         if (c->GetZ() > z+roadz) break;
861         if ( (c->GetY()-y) >  roady ) continue;
862         Float_t distance = (c->GetZ()-z)*(c->GetZ()-z)+(c->GetY()-y)*(c->GetY()-y);
863         if (maxdistance>distance) {
864           maxdistance = distance;
865           cl=c;
866           //      index=krow.GetIndex(i);       
867         }       
868       }              
869     }  // end of seearch
870     //update cluster information
871     if (cl){ 
872       cy[drow] = cl->GetY();
873       cz[drow] = cl->GetZ();
874       sy[drow] = ErrY2(&t,cl);
875       sz[drow] = ErrZ2(&t,cl);
876       hitted[drow] = kTRUE;
877       nclusters++;
878     }
879     else
880       hitted[drow] = kFALSE;
881   }
882   //if we have information - update track
883   if (nclusters>0){
884     Float_t sumyw0   = 0;
885     Float_t sumzw0   = 0;
886     Float_t sumyw   = 0;
887     Float_t sumzw   = 0;
888     Float_t sumrow  = 0;
889     for (Int_t i=0;i<step;i++)
890       if (hitted[i]){
891         sumyw0+= 1/sy[i];
892         sumzw0+= 1/sz[i];
893         //
894         sumyw+= (cy[i]-yy[i])/sy[i];
895         sumzw+= (cz[i]-zz[i])/sz[i];
896         sumrow+= i;
897       }
898     Float_t dy = sumyw/sumyw0;
899     Float_t dz = sumzw/sumzw0;
900     Float_t mrow = sumrow/nclusters+row0;
901     Float_t x = fSectors->GetX(mrow);
902     t.PropagateTo(x);
903     AliTPCclusterMI cvirtual;
904     cvirtual.SetZ(dz+t.GetZ());
905     cvirtual.SetY(dy+t.GetY());
906     t.SetErrorY2(1.2*t.fErrorY2/TMath::Sqrt(Float_t(nclusters)));
907     t.SetErrorZ2(1.2*t.fErrorZ2/TMath::Sqrt(Float_t(nclusters)));
908     Float_t chi2 = t.GetPredictedChi2(&cvirtual);
909     t.Update(&cvirtual,chi2,0);
910     Int_t ncl = t.GetNumberOfClusters();
911     ncl = ncl-1+nclusters;
912     t.SetN(ncl);
913   }     
914   return  1;
915 }   
916 */
917
918 //_____________________________________________________________________________
919 Int_t AliTPCtrackerMI::FollowProlongation(AliTPCseed& t, Int_t rf) {
920   //-----------------------------------------------------------------
921   // This function tries to find a track prolongation.
922   //-----------------------------------------------------------------
923   Double_t xt=t.GetX();
924   //
925   Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
926   if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();  
927   if (alpha < 0.            ) alpha += 2.*TMath::Pi();  
928   t.fRelativeSector = Int_t(alpha/fSectors->GetAlpha())%fN;
929     
930   for (Int_t nr=fSectors->GetRowNumber(xt)-1; nr>=rf; nr--) {
931    
932     if (FollowToNext(t,nr)==0) {
933     }
934   }   
935   return 1;
936 }
937
938
939 //_____________________________________________________________________________
940 Int_t AliTPCtrackerMI::FollowBackProlongation(AliTPCseed& t, Int_t rf) {
941   //-----------------------------------------------------------------
942   // This function tries to find a track prolongation.
943   //-----------------------------------------------------------------
944   Double_t xt=t.GetX();
945   //
946   Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
947   if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();  
948   if (alpha < 0.            ) alpha += 2.*TMath::Pi();  
949   t.fRelativeSector = Int_t(alpha/fSectors->GetAlpha())%fN;
950     
951   for (Int_t nr=fSectors->GetRowNumber(xt)+1; nr<=rf; nr++) {
952     if (t.GetSnp()<0.8)
953       FollowToNext(t,nr);                                                             
954   }   
955   return 1;
956 }
957
958
959
960
961    
962 Float_t AliTPCtrackerMI::OverlapFactor(AliTPCseed * s1, AliTPCseed * s2, Int_t &sum1, Int_t & sum2)
963 {
964   //
965   //
966   sum1=0;
967   sum2=0;
968   Int_t sum=0;
969   if (s1->fSector!=s2->fSector) return 0;
970   //
971   Float_t dz2 =(s1->GetZ() - s2->GetZ());
972   dz2*=dz2;
973   Float_t dy2 =(s1->GetY() - s2->GetY());
974   dy2*=dy2;
975   Float_t distance = TMath::Sqrt(dz2+dy2);
976   if (distance>5.) return 0; // if there are far away  - not overlap - to reduce combinatorics
977  
978   for (Int_t i=0;i<160;i++){
979     if (s1->fClusterIndex[i]>0) sum1++;
980     if (s2->fClusterIndex[i]>0) sum2++;
981     if (s1->fClusterIndex[i]==s2->fClusterIndex[i] && s1->fClusterIndex[i]>0) {
982       sum++;
983     }
984   }
985  
986   Float_t summin = TMath::Min(sum1+1,sum2+1);
987   Float_t ratio = (sum+1)/Float_t(summin);
988   return ratio;
989 }
990
991 void  AliTPCtrackerMI::SignShared(AliTPCseed * s1, AliTPCseed * s2)
992 {
993   //
994   //
995   if (s1->fSector!=s2->fSector) return;
996   //
997   Float_t dz2 =(s1->GetZ() - s2->GetZ());
998   dz2*=dz2;
999   Float_t dy2 =(s1->GetY() - s2->GetY());
1000   dy2*=dy2;
1001   Float_t distance = TMath::Sqrt(dz2+dy2);
1002   if (distance>15.) return ; // if there are far away  - not overlap - to reduce combinatorics
1003   //trpoint = new (pointarray[track->fRow]) AliTPCTrackPoint;
1004   //  TClonesArray &pointarray1 = *(s1->fPoints);
1005   //TClonesArray &pointarray2 = *(s2->fPoints);
1006   //
1007   for (Int_t i=0;i<160;i++){
1008     if (s1->fClusterIndex[i]==s2->fClusterIndex[i] && s1->fClusterIndex[i]>0) {
1009       //  AliTPCTrackPoint *p1  = (AliTPCTrackPoint *)(pointarray1.UncheckedAt(i));
1010       //AliTPCTrackPoint *p2  = (AliTPCTrackPoint *)(pointarray2.UncheckedAt(i)); 
1011       AliTPCTrackPoint *p1  = s1->GetTrackPoint(i);
1012       AliTPCTrackPoint *p2  = s2->GetTrackPoint(i);; 
1013       p1->fIsShared = kTRUE;
1014       p2->fIsShared = kTRUE;
1015     }
1016   } 
1017 }
1018
1019
1020
1021
1022 void  AliTPCtrackerMI::RemoveOverlap(TObjArray * arr, Float_t factor, Int_t removalindex , Bool_t shared){
1023
1024   
1025
1026   //
1027   // remove overlap - used removal factor - removal index stored in the track
1028   arr->Sort();  // sorting according z
1029   arr->Expand(arr->GetEntries());
1030   Int_t nseed=arr->GetEntriesFast();
1031   //  printf("seeds \t%p \t%d\n",arr, nseed);
1032   //  arr->Expand(arr->GetEntries());  //remove 0 pointers
1033   nseed = arr->GetEntriesFast();
1034   Int_t removed = 0;
1035   for (Int_t i=0; i<nseed; i++) {
1036     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);    
1037     if (!pt) {
1038       continue;
1039     }
1040     if (!(pt->IsActive())) continue;
1041     for (Int_t j=i+1; j<nseed; j++){
1042       AliTPCseed *pt2=(AliTPCseed*)arr->UncheckedAt(j);
1043       if ((pt2) && pt2->IsActive())
1044         if (pt->fSector == pt2->fSector)
1045           if (TMath::Abs(pt2->GetZ()-pt->GetZ())<2){
1046             Int_t sum1,sum2;
1047             Float_t ratio = OverlapFactor(pt,pt2,sum1,sum2);
1048             //if (sum1==0) {
1049             //  pt->Desactivate(removalindex); // arr->RemoveAt(i); 
1050             //  break;
1051             //}
1052             if (ratio>factor){
1053               //          if (pt->GetChi2()<pt2->GetChi2()) pt2->Desactivate(removalindex);  // arr->RemoveAt(j);             
1054               Float_t ratio2 = (pt->GetChi2()*sum2)/(pt2->GetChi2()*sum1);
1055               Float_t ratio3 = Float_t(sum1-sum2)/Float_t(sum1+sum2);
1056               removed++;
1057               if (TMath::Abs(ratio3)>0.025){  // if much more points  
1058                 if (sum1>sum2) pt2->Desactivate(removalindex);
1059                 else {
1060                   pt->Desactivate(removalindex); // arr->RemoveAt(i); 
1061                   break;
1062                 }
1063               }
1064               else{  //decide on mean chi2
1065                 if (ratio2<1)  
1066                   pt2->Desactivate(removalindex);
1067                 else {
1068                   pt->Desactivate(removalindex); // arr->RemoveAt(i); 
1069                   break;
1070                 }           
1071               }  
1072               
1073             }  // if suspicious ratio
1074           }
1075           else
1076             break;
1077     }
1078   }
1079   //  printf("removed\t%d\n",removed);
1080   Int_t good =0; 
1081   for (Int_t i=0; i<nseed; i++) {
1082     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);    
1083     if (!pt) break;
1084     if (pt->GetNumberOfClusters() < pt->fNFoundable*0.5) {
1085       //desactivate tracks with small number of points
1086       //      printf("%d\t%d\t%f\n", pt->GetNumberOfClusters(), pt->fNFoundable,pt->GetNumberOfClusters()/Float_t(pt->fNFoundable));
1087       pt->Desactivate(10);  //desactivate  - small muber of points
1088     }
1089     if (!(pt->IsActive())) continue;
1090     good++;
1091   }
1092   
1093   
1094   if (shared)
1095     for (Int_t i=0; i<nseed; i++) {
1096       AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);    
1097       if (!pt) continue;
1098       if (!(pt->IsActive())) continue;
1099       for (Int_t j=i+1; j<nseed; j++){
1100         AliTPCseed *pt2=(AliTPCseed*)arr->UncheckedAt(j);
1101         if ((pt2) && pt2->IsActive()) {
1102           if ( TMath::Abs(pt->fSector-pt2->fSector)>1) break;
1103           SignShared(pt,pt2);
1104         }
1105       }
1106     }
1107   fNtracks = good;
1108   printf("\n*****\nNumber of good tracks after overlap removal\t%d\n",fNtracks);
1109
1110
1111 }
1112 void AliTPCtrackerMI::MakeSeedsAll()
1113 {
1114   if (fSeeds == 0) fSeeds = new TObjArray;
1115   TObjArray * arr;
1116   for (Int_t sec=0;sec<fkNOS;sec+=3){
1117      arr = MakeSeedsSectors(sec,sec+3);
1118      Int_t nseed = arr->GetEntriesFast();
1119      for (Int_t i=0;i<nseed;i++)
1120        fSeeds->AddLast(arr->RemoveAt(i));
1121   }  
1122   //  fSeeds = MakeSeedsSectors(0,fkNOS);
1123 }
1124
1125 TObjArray *  AliTPCtrackerMI::MakeSeedsSectors(Int_t sec1, Int_t sec2)
1126 {
1127   //
1128   // loop over all  sectors and make seed
1129   //find track seeds
1130   TStopwatch timer;
1131   Int_t nup=fOuterSec->GetNRows(), nlow=fInnerSec->GetNRows();
1132   Int_t nrows=nlow+nup;  
1133   Int_t gap=Int_t(0.125*nrows), shift=Int_t(0.5*gap);
1134   //  if (fSeeds==0) fSeeds = new TObjArray;
1135   TObjArray * arr = new TObjArray;
1136
1137   for (Int_t sec=sec1; sec<sec2;sec++){
1138     MakeSeeds(arr, sec, nup-1, nup-1-gap);
1139     MakeSeeds(arr, sec, nup-1-shift, nup-1-shift-gap);
1140   }
1141   gap = Int_t(0.3* nrows);
1142   for (Int_t sec=sec1; sec<sec2;sec++){
1143     //find secondaries
1144     MakeSeeds2(arr, sec, nup-1, nup-1-gap);
1145     MakeSeeds2(arr, sec, nup-1-shift, nup-1-shift-gap);
1146     MakeSeeds2(arr, sec, nup-1-2*shift, nup-1-2*shift-gap);
1147     //MakeSeeds2(arr, sec, nup-1-3*shift, nup-1-3*shift-gap);
1148     MakeSeeds2(arr, sec, 30, 0);
1149   }
1150   
1151   Int_t nseed=arr->GetEntriesFast();
1152   gap=Int_t(0.3*nrows);
1153   // continue seeds 
1154   Int_t i;
1155   for (i=0; i<nseed; i++) {
1156     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i), &t=*pt; 
1157     if (!pt) continue;
1158     if (FollowProlongation(t,nup-gap)) {
1159       pt->fIsSeeding =kFALSE;
1160       continue;
1161     }
1162     delete arr->RemoveAt(i);
1163   }
1164      
1165   //
1166   //remove seeds which overlaps  
1167   RemoveOverlap(arr,0.6,1);       
1168   //delete seeds - which were sign  
1169   nseed=arr->GetEntriesFast();
1170   for (i=0; i<nseed; i++) {
1171     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
1172     //, &t=*pt;
1173     if (!pt) continue;
1174     if ((pt->IsActive())  &&  pt->GetNumberOfClusters() > pt->fNFoundable*0.5 ) {
1175       //pt->Reset();
1176       //FollowBackProlongation(*pt,nup-1);
1177       //if ( pt->GetNumberOfClusters() < pt->fNFoundable*0.5 || pt->GetNumberOfClusters()<10 )
1178       //delete arr->RemoveAt(i);
1179       //else
1180       //        pt->Reset();
1181       continue;
1182     }
1183     delete arr->RemoveAt(i);    
1184   } 
1185   //RemoveOverlap(arr,0.6,1);
1186   return arr;
1187 }
1188
1189
1190
1191 //_____________________________________________________________________________
1192 void AliTPCtrackerMI::MakeSeeds(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2) {
1193   //-----------------------------------------------------------------
1194   // This function creates track seeds.
1195   //-----------------------------------------------------------------
1196   //  if (fSeeds==0) fSeeds=new TObjArray(15000);
1197
1198   Double_t x[5], c[15];
1199
1200   Double_t alpha=fOuterSec->GetAlpha(), shift=fOuterSec->GetAlphaShift();
1201   Double_t cs=cos(alpha), sn=sin(alpha);
1202
1203   Double_t x1 =fOuterSec->GetX(i1);
1204   Double_t xx2=fOuterSec->GetX(i2);
1205   //
1206   //  for (Int_t ns=0; ns<fkNOS; ns++) 
1207   Int_t ns =sec;
1208     {
1209     Int_t nl=fOuterSec[(ns-1+fkNOS)%fkNOS][i2];
1210     Int_t nm=fOuterSec[ns][i2];
1211     Int_t nu=fOuterSec[(ns+1)%fkNOS][i2];
1212     const AliTPCRow& kr1=fOuterSec[ns][i1];
1213     AliTPCRow&  kr21 = fOuterSec[(ns-1+fkNOS)%fkNOS][i2];
1214     AliTPCRow&  kr22 = fOuterSec[(ns)%fkNOS][i2];
1215     AliTPCRow&  kr23 = fOuterSec[(ns+1)%fkNOS][i2];
1216
1217     for (Int_t is=0; is < kr1; is++) {
1218       Double_t y1=kr1[is]->GetY(), z1=kr1[is]->GetZ();
1219       Double_t x3=GetX(), y3=GetY(), z3=GetZ();
1220
1221       Float_t anglez = (z1-z3)/(x1-x3); 
1222       Float_t extraz = z1 - anglez*(x1-xx2);  // extrapolated z
1223
1224       for (Int_t js=0; js < nl+nm+nu; js++) {
1225         const AliTPCclusterMI *kcl;
1226         Double_t x2,   y2,   z2;
1227         if (js<nl) {     
1228           if (js==0) {
1229             js = kr21.Find(extraz-15.);
1230             if (js>=nl) continue;
1231           }       
1232           kcl=kr21[js];
1233           z2=kcl->GetZ();
1234           if ((extraz-z2)>10) continue;   
1235           if ((extraz-z2)<-10) {
1236             js = nl-1;
1237             continue;
1238           }
1239           y2=kcl->GetY(); 
1240           x2= xx2*cs+y2*sn;
1241           y2=-xx2*sn+y2*cs;
1242         } else 
1243           if (js<nl+nm) {
1244             if (js==nl) {
1245               js = nl+kr22.Find(extraz-15.);
1246               if (js>=nl+nm) continue;
1247             }             
1248             kcl=kr22[js-nl];
1249             z2=kcl->GetZ();
1250             if ((extraz-z2)>10) continue;
1251             if ((extraz-z2)<-10) {
1252               js = nl+nm-1;
1253               continue;
1254             }
1255             x2=xx2; y2=kcl->GetY(); 
1256           } else {
1257             //const AliTPCRow& kr2=fOuterSec[(ns+1)%fkNOS][i2];   
1258             if (js==nl+nm) {
1259               js = nl+nm+kr23.Find(extraz-15.);
1260               if (js>=nl+nm+nu) break;
1261             }     
1262             kcl=kr23[js-nl-nm];
1263             z2=kcl->GetZ(); 
1264             if ((extraz-z2)>10) continue;
1265             if ((extraz-z2)<-10) {
1266               break;        
1267             }
1268             y2=kcl->GetY();
1269             x2=xx2*cs-y2*sn;
1270             y2=xx2*sn+y2*cs;
1271           }
1272
1273         Double_t zz=z1 - anglez*(x1-x2); 
1274         if (TMath::Abs(zz-z2)>10.) continue;
1275
1276         Double_t d=(x2-x1)*(0.-y2)-(0.-x2)*(y2-y1);
1277         if (d==0.) {cerr<<"MakeSeeds warning: Straight seed !\n"; continue;}
1278
1279         x[0]=y1;
1280         x[1]=z1;
1281         x[4]=f1(x1,y1,x2,y2,x3,y3);
1282         if (TMath::Abs(x[4]) >= 0.0066) continue;
1283         x[2]=f2(x1,y1,x2,y2,x3,y3);
1284         //if (TMath::Abs(x[4]*x1-x[2]) >= 0.99999) continue;
1285         x[3]=f3(x1,y1,x2,y2,z1,z2);
1286         if (TMath::Abs(x[3]) > 1.2) continue;
1287         Double_t a=asin(x[2]);
1288         Double_t zv=z1 - x[3]/x[4]*(a+asin(x[4]*x1-x[2]));
1289         if (TMath::Abs(zv-z3)>10.) continue; 
1290
1291         Double_t sy1=kr1[is]->GetSigmaY2()*2, sz1=kr1[is]->GetSigmaZ2()*4;
1292         Double_t sy2=kcl->GetSigmaY2()*2,     sz2=kcl->GetSigmaZ2()*4;
1293         //Double_t sy3=400*3./12., sy=0.1, sz=0.1;
1294         Double_t sy3=25000*x[4]*x[4]+0.1, sy=0.1, sz=0.1;
1295         //Double_t sy3=25000*x[4]*x[4]*60+0.5, sy=0.1, sz=0.1;
1296
1297         Double_t f40=(f1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
1298         Double_t f42=(f1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
1299         Double_t f43=(f1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
1300         Double_t f20=(f2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
1301         Double_t f22=(f2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
1302         Double_t f23=(f2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
1303         Double_t f30=(f3(x1,y1+sy,x2,y2,z1,z2)-x[3])/sy;
1304         Double_t f31=(f3(x1,y1,x2,y2,z1+sz,z2)-x[3])/sz;
1305         Double_t f32=(f3(x1,y1,x2,y2+sy,z1,z2)-x[3])/sy;
1306         Double_t f34=(f3(x1,y1,x2,y2,z1,z2+sz)-x[3])/sz;
1307
1308         c[0]=sy1;
1309         c[1]=0.;       c[2]=sz1;
1310         c[3]=f20*sy1;  c[4]=0.;       c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
1311         c[6]=f30*sy1;  c[7]=f31*sz1;  c[8]=f30*sy1*f20+f32*sy2*f22;
1312                        c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
1313         c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
1314         c[13]=f30*sy1*f40+f32*sy2*f42;
1315         c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
1316
1317         UInt_t index=kr1.GetIndex(is);
1318         AliTPCseed *track=new AliTPCseed(index, x, c, x1, ns*alpha+shift);
1319         track->fIsSeeding = kTRUE;
1320         Int_t rc=FollowProlongation(*track, i2);
1321         //FollowProlongationFast(*track, 5);
1322         //FollowProlongationFast(*track, 5);
1323         //FollowProlongationFast(*track, 5);
1324         //FollowProlongationFast(*track, 5);
1325         //Int_t rc = 1;
1326         
1327         track->fLastPoint = i1;  // first cluster in track position
1328         if (rc==0 || track->GetNumberOfClusters()<(i1-i2)/4 || track->GetNumberOfClusters() < track->fNFoundable/2. ) delete track;
1329         else arr->AddLast(track); 
1330       }
1331     }
1332   }
1333 }
1334
1335
1336 //_____________________________________________________________________________
1337 void AliTPCtrackerMI::MakeSeeds2(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2) {
1338   //-----------------------------------------------------------------
1339   // This function creates track seeds - without vertex constraint
1340   //-----------------------------------------------------------------
1341
1342   Double_t alpha=fOuterSec->GetAlpha(), shift=fOuterSec->GetAlphaShift();
1343   //  Double_t cs=cos(alpha), sn=sin(alpha);
1344   Int_t row0 = (i1+i2)/2;
1345   Int_t drow = (i1-i2)/2;
1346   const AliTPCRow& kr0=fSectors[sec][row0];
1347   const AliTPCRow& krm=fSectors[sec][row0-1];
1348   const AliTPCRow& krp=fSectors[sec][row0+1];
1349   AliTPCRow * kr=0;
1350
1351   AliTPCpolyTrack polytrack;
1352   Int_t nclusters=fSectors[sec][row0];
1353
1354   for (Int_t is=0; is < nclusters; is++) {
1355     const AliTPCclusterMI * cl= kr0[is];
1356     Double_t x = kr0.GetX();
1357
1358     // Initialization of the polytrack
1359     polytrack.Reset();
1360
1361     Double_t y0= cl->GetY();
1362     Double_t z0= cl->GetZ();
1363     polytrack.AddPoint(x,y0,z0);
1364     Float_t roady = 5*TMath::Sqrt(cl->GetSigmaY2()+0.2);
1365     Float_t roadz = 5*TMath::Sqrt(cl->GetSigmaZ2()+0.2);
1366     //
1367     x = krm.GetX();
1368     cl = krm.FindNearest(y0,z0,roady,roadz);
1369     if (cl) polytrack.AddPoint(x,cl->GetY(),cl->GetZ(),roady,roadz);
1370     //
1371     x = krp.GetX();
1372     cl = krp.FindNearest(y0,z0,roady,roadz);
1373     if (cl) polytrack.AddPoint(x,cl->GetY(),cl->GetZ(),cl->GetSigmaY2()+0.05,cl->GetSigmaZ2()+0.05);
1374     //
1375     polytrack.UpdateParameters();
1376     // follow polytrack
1377     roadz = 0.6;
1378     roady = 0.6;
1379     //
1380     Double_t yn,zn;
1381     Int_t nfoundable = polytrack.GetN();
1382     Int_t nfound     = nfoundable; 
1383     for (Int_t ddrow = 2; ddrow<drow;ddrow++){
1384       for (Int_t delta = -1;delta<=1;delta+=2){
1385         Int_t row = row0+ddrow*delta;
1386         kr = &(fSectors[sec][row]);
1387         Double_t xn = kr->GetX();
1388         Double_t ymax = fSectors->GetMaxY(row)-kr->fDeadZone;
1389         polytrack.GetFitPoint(xn,yn,zn);
1390         if (TMath::Abs(yn)>ymax) continue;
1391         nfoundable++;
1392         AliTPCclusterMI * cln = kr->FindNearest(yn,zn,roady,roadz);
1393         if (cln) {
1394           polytrack.AddPoint(xn,cln->GetY(),cln->GetZ(),cln->GetSigmaY2()+0.05,cln->GetSigmaZ2()+0.05);
1395           nfound++;
1396         }
1397       }
1398       polytrack.UpdateParameters();
1399     }
1400     if ((nfound>0.5*nfoundable) &&( nfoundable>0.4*(i1-i2))) {
1401       // add polytrack candidate
1402       Double_t x[5], c[15];
1403       Double_t x1,x2,x3,y1,y2,y3,z1,z2,z3;
1404       polytrack.GetBoundaries(x3,x1);
1405       x2 = (x1+x3)/2.;
1406       polytrack.GetFitPoint(x1,y1,z1);
1407       polytrack.GetFitPoint(x2,y2,z2);
1408       polytrack.GetFitPoint(x3,y3,z3);
1409       //
1410       //is track pointing to the vertex ?
1411       Double_t x0,y0,z0;
1412       x0=0;
1413       polytrack.GetFitPoint(x0,y0,z0);
1414       if ( (TMath::Abs(z0-GetZ())<10) && (TMath::Abs(y0-GetY())<5)){ //if yes apply vertex constraint
1415         x3 = 0;
1416         y3 = GetY();
1417         z3 = GetZ();
1418       }
1419
1420       x[0]=y1;
1421       x[1]=z1;
1422       x[4]=f1(x1,y1,x2,y2,x3,y3);
1423       if (TMath::Abs(x[4]) >= 0.0066) continue;
1424       x[2]=f2(x1,y1,x2,y2,x3,y3);
1425       //if (TMath::Abs(x[4]*x1-x[2]) >= 0.99999) continue;
1426       x[3]=f3(x1,y1,x2,y2,z1,z2);
1427       if (TMath::Abs(x[3]) > 1.2) continue;
1428       if (TMath::Abs(x[2]) > 0.99) continue;
1429       //      Double_t a=asin(x[2]);
1430       
1431
1432       Double_t sy=1.5, sz=1.5;
1433       Double_t sy1=1.5, sz1=1.5;
1434       Double_t sy2=1.3, sz2=1.3;
1435       Double_t sy3=1.5;
1436       //sz3=1.5;
1437       //Double_t sy3=400*3./12., sy=0.1, sz=0.1;
1438       //      Double_t sy3=25000*x[4]*x[4]+0.1, sy=0.1, sz=0.1;
1439       //Double_t sy3=25000*x[4]*x[4]*60+0.5, sy=0.1, sz=0.1;
1440       
1441       Double_t f40=(f1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
1442       Double_t f42=(f1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
1443       Double_t f43=(f1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
1444       Double_t f20=(f2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
1445       Double_t f22=(f2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
1446       Double_t f23=(f2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
1447       Double_t f30=(f3(x1,y1+sy,x2,y2,z1,z2)-x[3])/sy;
1448       Double_t f31=(f3(x1,y1,x2,y2,z1+sz,z2)-x[3])/sz;
1449       Double_t f32=(f3(x1,y1,x2,y2+sy,z1,z2)-x[3])/sy;
1450       Double_t f34=(f3(x1,y1,x2,y2,z1,z2+sz)-x[3])/sz;
1451       
1452       c[0]=sy1;
1453       c[1]=0.;       c[2]=sz1;
1454       c[3]=f20*sy1;  c[4]=0.;       c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
1455       c[6]=f30*sy1;  c[7]=f31*sz1;  c[8]=f30*sy1*f20+f32*sy2*f22;
1456       c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
1457       c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
1458       c[13]=f30*sy1*f40+f32*sy2*f42;
1459       c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
1460       
1461       UInt_t index=0;
1462       //kr0.GetIndex(is);
1463       AliTPCseed *track=new AliTPCseed(index, x, c, x1, sec*alpha+shift);
1464       track->fStopped =kFALSE;
1465       track->fIsSeeding = kTRUE;
1466       Int_t rc=FollowProlongation(*track, i2);  
1467         track->fLastPoint = i1;  // first cluster in track position
1468         if (rc==0 || track->GetNumberOfClusters()<(i1-i2)/4 || track->GetNumberOfClusters() < track->fNFoundable/2. ) delete track;
1469         else arr->AddLast(track); 
1470     }
1471   }
1472   
1473   
1474   
1475 }
1476
1477
1478
1479
1480
1481 //_____________________________________________________________________________
1482 Int_t AliTPCtrackerMI::ReadSeeds(const TFile *inp) {
1483   //-----------------------------------------------------------------
1484   // This function reades track seeds.
1485   //-----------------------------------------------------------------
1486   TDirectory *savedir=gDirectory; 
1487
1488   TFile *in=(TFile*)inp;
1489   if (!in->IsOpen()) {
1490      cerr<<"AliTPCtrackerMI::ReadSeeds(): input file is not open !\n";
1491      return 1;
1492   }
1493
1494   in->cd();
1495   TTree *seedTree=(TTree*)in->Get("Seeds");
1496   if (!seedTree) {
1497      cerr<<"AliTPCtrackerMI::ReadSeeds(): ";
1498      cerr<<"can't get a tree with track seeds !\n";
1499      return 2;
1500   }
1501   AliTPCtrack *seed=new AliTPCtrack; 
1502   seedTree->SetBranchAddress("tracks",&seed);
1503   
1504   if (fSeeds==0) fSeeds=new TObjArray(15000);
1505
1506   Int_t n=(Int_t)seedTree->GetEntries();
1507   for (Int_t i=0; i<n; i++) {
1508      seedTree->GetEvent(i);
1509      fSeeds->AddLast(new AliTPCseed(*seed,seed->GetAlpha()));
1510   }
1511   
1512   delete seed;
1513   delete seedTree; 
1514   savedir->cd();
1515   return 0;
1516 }
1517
1518 //_____________________________________________________________________________
1519 Int_t AliTPCtrackerMI::Clusters2Tracks(const TFile *inp, TFile *out) {
1520   //-----------------------------------------------------------------
1521   // This is a track finder.
1522   //-----------------------------------------------------------------
1523   TDirectory *savedir=gDirectory; 
1524
1525   if (inp) {
1526      TFile *in=(TFile*)inp;
1527      if (!in->IsOpen()) {
1528         cerr<<"AliTPCtrackerMI::Clusters2Tracks(): input file is not open !\n";
1529         return 1;
1530      }
1531   }
1532
1533   if (!out->IsOpen()) {
1534      cerr<<"AliTPCtrackerMI::Clusters2Tracks(): output file is not open !\n";
1535      return 2;
1536   }
1537
1538   out->cd();
1539
1540   char   tname[100];
1541   sprintf(tname,"TreeT_TPC_%d",fEventN);
1542   TTree tracktree(tname,"Tree with TPC tracks");
1543   TTree seedtree("Seeds","Seeds");
1544   AliTPCtrack *iotrack=0;
1545   AliTPCseed  *ioseed=0;
1546   tracktree.Branch("tracks","AliTPCtrack",&iotrack,32000,0);
1547   TStopwatch timer;
1548   
1549   printf("Loading clusters \n");
1550   LoadClusters();
1551   printf("Time for loading clusters: \t");timer.Print();timer.Start();
1552
1553   printf("Loading outer sectors\n");
1554   LoadOuterSectors();
1555   printf("Time for loading outer sectors: \t");timer.Print();timer.Start();
1556
1557   printf("Loading inner sectors\n");
1558   LoadInnerSectors();
1559   printf("Time for loading inner sectors: \t");timer.Print();timer.Start();
1560   fSectors = fOuterSec;
1561   fN=fkNOS;
1562   
1563   
1564   //find track seeds
1565   MakeSeedsAll(); 
1566   printf("Time for seeding: \t"); timer.Print();timer.Start();
1567   Int_t nup=fOuterSec->GetNRows(), nlow=fInnerSec->GetNRows();
1568   Int_t nrows=nlow+nup;
1569   
1570   Int_t gap=Int_t(0.3*nrows);
1571   Int_t i;
1572   //RemoveOverlap(fSeeds,0.6,2);
1573   Int_t nseed=fSeeds->GetEntriesFast();
1574   // outer sectors parallel tracking
1575   ParallelTracking(fSectors->GetNRows()-gap-1,0); 
1576   //ParallelTracking(fSectors->GetNRows()-1,0); 
1577   //RemoveOverlap(fSeeds, 0.6,3);       
1578   //  ParallelTracking(49,0); 
1579   printf("Time for parralel tracking outer sectors: \t"); timer.Print();timer.Start();
1580
1581   RemoveOverlap(fSeeds, 0.6,3);  
1582   printf("Time for removal overlap- outer sectors: \t");timer.Print();timer.Start();
1583   //parallel tracking 
1584   fSectors = fInnerSec;
1585   fN=fkNIS;
1586   ParallelTracking(fSectors->GetNRows()-1,0);
1587   printf("Number of tracks after  inner tracking  %d\n",fNtracks); 
1588   printf("Time for parralel tracking inner sectors: \t"); timer.Print();timer.Start();
1589   //
1590   RemoveOverlap(fSeeds,0.6,5,kTRUE);  // remove overlap -  shared points signed 
1591   printf("Time for removal overlap- inner sectors: \t"); timer.Print();timer.Start();
1592   //
1593   // 
1594   
1595   ioseed  = (AliTPCseed*)(fSeeds->UncheckedAt(0));
1596   AliTPCseed * vseed = new AliTPCseed;
1597   vseed->fPoints = new TClonesArray("AliTPCTrackPoint",1);
1598   vseed->fEPoints = new TClonesArray("AliTPCExactPoint",1);
1599   vseed->fPoints->ExpandCreateFast(2);
1600   
1601   TBranch * seedbranch =   seedtree.Branch("seeds","AliTPCseed",&vseed,32000,99);
1602   //delete vseed;
1603   nseed=fSeeds->GetEntriesFast();
1604
1605   Int_t found = 0;
1606   for (i=0; i<nseed; i++) {
1607     AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;    
1608     if (!pt) continue;    
1609     Int_t nc=t.GetNumberOfClusters();
1610     if (nc<20) continue;
1611     t.CookdEdx(0.02,0.6);
1612     CookLabel(pt,0.1); //For comparison only
1613     //   if ((pt->IsActive()) && (nc>Int_t(0.4*nrows))){
1614     if ((pt->IsActive()) && (nc>Int_t(0.5*t.fNFoundable) && (t.fNFoundable>Int_t(0.3*nrows)))){
1615       iotrack=pt;
1616       tracktree.Fill();
1617      cerr<<found++<<'\r';      
1618     }   
1619     else 
1620       if ( (pt->IsActive())) fNtracks--;
1621     pt->RebuildSeed();
1622     seedbranch->SetAddress(&pt);
1623
1624     seedtree.Fill();        
1625     for (Int_t j=0;j<160;j++){
1626       delete pt->fPoints->RemoveAt(j);
1627     }
1628     delete pt->fPoints;
1629     pt->fPoints =0;
1630     delete fSeeds->RemoveAt(i);
1631   }
1632   printf("Time for track writing and dedx cooking: \t"); timer.Print();timer.Start();
1633
1634   UnloadClusters();
1635   printf("Time for unloading cluster: \t"); timer.Print();timer.Start();
1636
1637   tracktree.Write();
1638   seedtree.Write();
1639   cerr<<"Number of found tracks : "<<fNtracks<<"\t"<<found<<endl;
1640   
1641   savedir->cd();
1642   
1643   return 0;
1644 }
1645
1646
1647 void  AliTPCtrackerMI::ParallelTracking(Int_t rfirst, Int_t rlast)
1648 {
1649   //
1650   // try to track in parralel
1651
1652   Int_t nseed=fSeeds->GetEntriesFast();
1653   //prepare seeds for tracking
1654   for (Int_t i=0; i<nseed; i++) {
1655     AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt; 
1656     if (!pt) continue;
1657     if (!t.IsActive()) continue;
1658     // follow prolongation to the first layer
1659     FollowProlongation(t, rfirst+1);
1660   }
1661
1662
1663   //
1664   for (Int_t nr=rfirst; nr>=rlast; nr--){      
1665     // make indexes with the cluster tracks for given       
1666     //    for (Int_t i = 0;i<fN;i++)
1667     //  fSectors[i][nr].MakeClusterTracks();
1668
1669     // find nearest cluster
1670     for (Int_t i=0; i<nseed; i++) {
1671       AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt; 
1672       if (!pt) continue;
1673       if (!pt->IsActive()) continue;
1674       if (pt->fRelativeSector>17) {
1675         continue;
1676       }
1677       UpdateClusters(t,i,nr);
1678     }
1679     // prolonagate to the nearest cluster - if founded
1680     for (Int_t i=0; i<nseed; i++) {
1681       AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i); 
1682       if (!pt) continue;
1683       if (!pt->IsActive()) continue; 
1684       if (pt->fRelativeSector>17) {
1685         continue;
1686       }
1687       FollowToNextCluster(i,nr);
1688     }
1689     //    for (Int_t i= 0;i<fN;i++)
1690     //  fSectors[i][nr].ClearClusterTracks();
1691   }  
1692   
1693 }
1694
1695 Float_t  AliTPCtrackerMI::GetSigmaY(AliTPCseed * seed)
1696 {
1697   //
1698   //  
1699   Float_t sd2 = (fParam->GetZLength()-TMath::Abs(seed->GetZ()))*fParam->GetDiffL()*fParam->GetDiffL();
1700   Float_t padlength =  fParam->GetPadPitchLength(seed->fSector);
1701   Float_t sres = (seed->fSector < fParam->GetNSector()/2) ? 0.2 :0.3;
1702   Float_t angular  = seed->GetSnp();
1703   angular = angular*angular/(1-angular*angular);
1704   //  angular*=angular;
1705   //angular  = TMath::Sqrt(angular/(1-angular));
1706   Float_t res = TMath::Sqrt(sd2+padlength*padlength*angular/12.+sres*sres);
1707   return res;
1708 }
1709 Float_t  AliTPCtrackerMI::GetSigmaZ(AliTPCseed * seed)
1710 {
1711   //
1712   //
1713   Float_t sd2 = (fParam->GetZLength()-TMath::Abs(seed->GetZ()))*fParam->GetDiffL()*fParam->GetDiffL();
1714   Float_t padlength =  fParam->GetPadPitchLength(seed->fSector);
1715   Float_t sres = fParam->GetZSigma();
1716   Float_t angular  = seed->GetTgl();
1717   Float_t res = TMath::Sqrt(sd2+padlength*padlength*angular*angular/12.+sres*sres);
1718   return res;
1719 }
1720
1721
1722
1723 //_________________________________________________________________________
1724 AliTPCclusterMI *AliTPCtrackerMI::GetClusterMI(Int_t index) const {
1725   //--------------------------------------------------------------------
1726   //       Return pointer to a given cluster
1727   //--------------------------------------------------------------------
1728   Int_t sec=(index&0xff000000)>>24; 
1729   Int_t row=(index&0x00ff0000)>>16; 
1730   Int_t ncl=(index&0x0000ffff)>>00;
1731
1732   AliTPCClustersRow *clrow=((AliTPCtrackerMI *) this)->fClustersArray.GetRow(sec,row);
1733   if (!clrow) return 0;
1734   return (AliTPCclusterMI*)(*clrow)[ncl];      
1735 }
1736
1737 //__________________________________________________________________________
1738 void AliTPCtrackerMI::CookLabel(AliKalmanTrack *t, Float_t wrong) const {
1739   //--------------------------------------------------------------------
1740   //This function "cooks" a track label. If label<0, this track is fake.
1741   //--------------------------------------------------------------------
1742   Int_t noc=t->GetNumberOfClusters();
1743   Int_t *lb=new Int_t[noc];
1744   Int_t *mx=new Int_t[noc];
1745   AliTPCclusterMI **clusters=new AliTPCclusterMI*[noc];
1746
1747   Int_t i;
1748   for (i=0; i<noc; i++) {
1749      lb[i]=mx[i]=0;
1750      Int_t index=t->GetClusterIndex(i);
1751      clusters[i]=GetClusterMI(index);
1752   }
1753
1754   Int_t lab=123456789;
1755   for (i=0; i<noc; i++) {
1756     AliTPCclusterMI *c=clusters[i];
1757     if (!clusters[i]) continue;
1758     lab=TMath::Abs(c->GetLabel(0));
1759     Int_t j;
1760     for (j=0; j<noc; j++) if (lb[j]==lab || mx[j]==0) break;
1761     lb[j]=lab;
1762     (mx[j])++;
1763   }
1764
1765   Int_t max=0;
1766   for (i=0; i<noc; i++) if (mx[i]>max) {max=mx[i]; lab=lb[i];}
1767     
1768   for (i=0; i<noc; i++) {
1769     AliTPCclusterMI *c=clusters[i]; 
1770     if (!clusters[i]) continue;
1771     if (TMath::Abs(c->GetLabel(1)) == lab ||
1772         TMath::Abs(c->GetLabel(2)) == lab ) max++;
1773   }
1774
1775   if ((1.- Float_t(max)/noc) > wrong) lab=-lab;
1776
1777   else {
1778      Int_t tail=Int_t(0.10*noc);
1779      max=0;
1780      for (i=1; i<=tail; i++) {
1781        AliTPCclusterMI *c=clusters[noc-i];
1782        if (!clusters[i]) continue;
1783        if (lab == TMath::Abs(c->GetLabel(0)) ||
1784            lab == TMath::Abs(c->GetLabel(1)) ||
1785            lab == TMath::Abs(c->GetLabel(2))) max++;
1786      }
1787      if (max < Int_t(0.5*tail)) lab=-lab;
1788   }
1789
1790   t->SetLabel(lab);
1791
1792   delete[] lb;
1793   delete[] mx;
1794   delete[] clusters;
1795 }
1796
1797 //_________________________________________________________________________
1798 void AliTPCtrackerMI::AliTPCSector::Setup(const AliTPCParam *par, Int_t f) {
1799   //-----------------------------------------------------------------------
1800   // Setup inner sector
1801   //-----------------------------------------------------------------------
1802   if (f==0) {
1803      fAlpha=par->GetInnerAngle();
1804      fAlphaShift=par->GetInnerAngleShift();
1805      fPadPitchWidth=par->GetInnerPadPitchWidth();
1806      fPadPitchLength=par->GetInnerPadPitchLength();
1807      fN=par->GetNRowLow();
1808      fRow=new AliTPCRow[fN];
1809      for (Int_t i=0; i<fN; i++) {
1810        fRow[i].SetX(par->GetPadRowRadiiLow(i));
1811        fRow[i].fDeadZone =1.5;  //1.5 cm of dead zone
1812      }
1813   } else {
1814      fAlpha=par->GetOuterAngle();
1815      fAlphaShift=par->GetOuterAngleShift();
1816      fPadPitchWidth  = par->GetOuterPadPitchWidth();
1817      fPadPitchLength = par->GetOuter1PadPitchLength();
1818      f1PadPitchLength = par->GetOuter1PadPitchLength();
1819      f2PadPitchLength = par->GetOuter2PadPitchLength();
1820
1821      fN=par->GetNRowUp();
1822      fRow=new AliTPCRow[fN];
1823      for (Int_t i=0; i<fN; i++) {
1824        fRow[i].SetX(par->GetPadRowRadiiUp(i)); 
1825        fRow[i].fDeadZone =1.5;  // 1.5 cm of dead zone
1826      }
1827   } 
1828 }
1829
1830
1831 AliTPCtrackerMI::AliTPCRow::~AliTPCRow(){
1832   //
1833   if (fClusterTracks) delete [] fClusterTracks;
1834   fClusterTracks = 0;
1835 }
1836
1837 void AliTPCtrackerMI::AliTPCRow::MakeClusterTracks(){
1838   //create cluster tracks
1839   if (fN>0) 
1840     fClusterTracks = new AliTPCclusterTracks[fN];
1841 }
1842
1843 void AliTPCtrackerMI::AliTPCRow::ClearClusterTracks(){
1844   if (fClusterTracks) delete[] fClusterTracks;
1845   fClusterTracks =0;
1846 }
1847
1848
1849
1850 void AliTPCtrackerMI::AliTPCRow::UpdateClusterTrack(Int_t clindex, Int_t trindex, AliTPCseed * seed){
1851   //
1852   //
1853   // update information of the cluster tracks - if track is nearer then other tracks to the 
1854   // given track
1855   const AliTPCclusterMI * cl = (*this)[clindex];
1856   AliTPCclusterTracks * cltracks = GetClusterTracks(clindex);
1857   // find the distance of the cluster to the track
1858   Float_t dy2 = (cl->GetY()- seed->GetY());
1859   dy2*=dy2;
1860   Float_t dz2 = (cl->GetZ()- seed->GetZ());
1861   dz2*=dz2;
1862   //
1863   Float_t distance = TMath::Sqrt(dy2+dz2);
1864   if (distance> 3) 
1865     return;  // MI - to be changed - AliTPCtrackerParam
1866   
1867   if ( distance < cltracks->fDistance[0]){
1868     cltracks->fDistance[2] =cltracks->fDistance[1];
1869     cltracks->fDistance[1] =cltracks->fDistance[0];
1870     cltracks->fDistance[0] =distance;
1871     cltracks->fTrackIndex[2] =cltracks->fTrackIndex[1];
1872     cltracks->fTrackIndex[1] =cltracks->fTrackIndex[0];
1873     cltracks->fTrackIndex[0] =trindex; 
1874   }
1875   else
1876     if ( distance < cltracks->fDistance[1]){
1877       cltracks->fDistance[2] =cltracks->fDistance[1];  
1878       cltracks->fDistance[1] =distance;
1879       cltracks->fTrackIndex[2] =cltracks->fTrackIndex[1];
1880       cltracks->fTrackIndex[1] =trindex; 
1881     } else
1882       if (distance < cltracks->fDistance[2]){
1883         cltracks->fDistance[2] =distance;
1884         cltracks->fTrackIndex[2] =trindex;
1885       }  
1886 }
1887
1888
1889 //_________________________________________________________________________
1890 void 
1891 AliTPCtrackerMI::AliTPCRow::InsertCluster(const AliTPCclusterMI* c, UInt_t index) {
1892   //-----------------------------------------------------------------------
1893   // Insert a cluster into this pad row in accordence with its y-coordinate
1894   //-----------------------------------------------------------------------
1895   if (fN==kMaxClusterPerRow) {
1896     cerr<<"AliTPCRow::InsertCluster(): Too many clusters !\n"; return;
1897   }
1898   if (fN==0) {fIndex[0]=index; fClusters[fN++]=c; return;}
1899   Int_t i=Find(c->GetZ());
1900   memmove(fClusters+i+1 ,fClusters+i,(fN-i)*sizeof(AliTPCclusterMI*));
1901   memmove(fIndex   +i+1 ,fIndex   +i,(fN-i)*sizeof(UInt_t));
1902   fIndex[i]=index; fClusters[i]=c; fN++;
1903 }
1904
1905 //___________________________________________________________________
1906 Int_t AliTPCtrackerMI::AliTPCRow::Find(Double_t z) const {
1907   //-----------------------------------------------------------------------
1908   // Return the index of the nearest cluster 
1909   //-----------------------------------------------------------------------
1910   if (fN==0) return 0;
1911   if (z <= fClusters[0]->GetZ()) return 0;
1912   if (z > fClusters[fN-1]->GetZ()) return fN;
1913   Int_t b=0, e=fN-1, m=(b+e)/2;
1914   for (; b<e; m=(b+e)/2) {
1915     if (z > fClusters[m]->GetZ()) b=m+1;
1916     else e=m; 
1917   }
1918   return m;
1919 }
1920
1921
1922
1923 //___________________________________________________________________
1924 AliTPCclusterMI * AliTPCtrackerMI::AliTPCRow::FindNearest(Double_t y, Double_t z, Double_t roady, Double_t roadz) const {
1925   //-----------------------------------------------------------------------
1926   // Return the index of the nearest cluster in z y 
1927   //-----------------------------------------------------------------------
1928   Float_t maxdistance = roady*roady + roadz*roadz;
1929
1930   AliTPCclusterMI *cl =0;
1931   for (Int_t i=Find(z-roadz); i<fN; i++) {
1932       AliTPCclusterMI *c=(AliTPCclusterMI*)(fClusters[i]);
1933       if (c->GetZ() > z+roadz) break;
1934       if ( (c->GetY()-y) >  roady ) continue;
1935       Float_t distance = (c->GetZ()-z)*(c->GetZ()-z)+(c->GetY()-y)*(c->GetY()-y);
1936       if (maxdistance>distance) {
1937         maxdistance = distance;
1938         cl=c;       
1939       }
1940   }
1941   return cl;      
1942 }
1943
1944
1945
1946
1947
1948 AliTPCseed::AliTPCseed():AliTPCtrack(){
1949   //
1950   fRow=0; 
1951   fRemoval =0; 
1952   memset(fClusterIndex,0,sizeof(Int_t)*200);
1953   fPoints = 0;
1954   fEPoints = 0;
1955   fNFoundable =0;
1956   fNShared  =0;
1957   fTrackPoints =0;
1958   fRemoval = 0;
1959 }
1960
1961 AliTPCseed::AliTPCseed(const AliTPCtrack &t):AliTPCtrack(t){
1962   fPoints = 0;
1963   fEPoints = 0;
1964   fNShared  =0; 
1965   fTrackPoints =0;
1966   fRemoval =0;
1967 }
1968
1969 AliTPCseed::AliTPCseed(const AliKalmanTrack &t, Double_t a):AliTPCtrack(t,a){
1970   fRow=0;
1971   memset(fClusterIndex,0,sizeof(Int_t)*200); 
1972   fPoints = 0;
1973   fEPoints = 0;
1974   fNFoundable =0; 
1975   fNShared  =0; 
1976   fTrackPoints =0;
1977   fRemoval =0;
1978 }
1979
1980 AliTPCseed::AliTPCseed(UInt_t index, const Double_t xx[5], const Double_t cc[15], 
1981                                         Double_t xr, Double_t alpha):      
1982   AliTPCtrack(index, xx, cc, xr, alpha) {
1983   //
1984   //
1985   fRow =0;
1986   memset(fClusterIndex,0,sizeof(Int_t)*200); 
1987   fPoints = 0;
1988   fEPoints = 0;
1989   fNFoundable =0;
1990   fNShared  = 0;
1991   fTrackPoints =0;
1992   fRemoval =0;
1993 }
1994
1995 AliTPCseed::~AliTPCseed(){
1996   if (fPoints) delete fPoints;
1997   fPoints =0;
1998   fEPoints = 0;
1999   if (fTrackPoints){
2000     for (Int_t i=0;i<8;i++){
2001       delete [] fTrackPoints[i];
2002     }
2003     delete fTrackPoints;
2004     fTrackPoints =0;
2005   }
2006
2007 }
2008
2009 AliTPCTrackPoint * AliTPCseed::GetTrackPoint(Int_t i)
2010 {
2011   //
2012   // 
2013   if (!fTrackPoints) {
2014     fTrackPoints = new AliTPCTrackPoint*[8];
2015     for ( Int_t i=0;i<8;i++)
2016       fTrackPoints[i]=0;
2017   }
2018   Int_t index1 = i/20;
2019   if (!fTrackPoints[index1]) fTrackPoints[index1] = new AliTPCTrackPoint[20];
2020   return &(fTrackPoints[index1][i%20]);
2021 }
2022
2023 void AliTPCseed::RebuildSeed()
2024 {
2025   //
2026   // rebuild seed to be ready for storing
2027   fPoints = new TClonesArray("AliTPCTrackPoint",160);
2028   fPoints->ExpandCreateFast(160);
2029   fEPoints = new TClonesArray("AliTPCExactPoint",1);
2030   for (Int_t i=0;i<160;i++){
2031     AliTPCTrackPoint *trpoint = (AliTPCTrackPoint*)fPoints->UncheckedAt(i);
2032     *trpoint = *(GetTrackPoint(i));
2033   }
2034
2035 }
2036
2037 //_____________________________________________________________________________
2038 void AliTPCseed::CookdEdx(Double_t low, Double_t up) {
2039   //-----------------------------------------------------------------
2040   // This funtion calculates dE/dX within the "low" and "up" cuts.
2041   //-----------------------------------------------------------------
2042
2043   Float_t amp[200];
2044   Float_t angular[200];
2045   Float_t weight[200];
2046   Int_t index[200];
2047   //Int_t nc = 0;
2048   //  TClonesArray & arr = *fPoints; 
2049   Float_t meanlog = 100.;
2050   
2051   Float_t mean[4]  = {0,0,0,0};
2052   Float_t sigma[4] = {1000,1000,1000,1000};
2053   Int_t nc[4]      = {0,0,0,0};
2054   Float_t norm[4]    = {1000,1000,1000,1000};
2055   //
2056   //
2057   fNShared =0;
2058
2059   for (Int_t of =0; of<4; of++){    
2060     for (Int_t i=of;i<160;i+=4)
2061       {
2062         //AliTPCTrackPoint * point = (AliTPCTrackPoint *) arr.At(i);
2063         AliTPCTrackPoint * point = GetTrackPoint(i);
2064         if (point==0) continue;
2065         if (point->fIsShared){
2066           fNShared++;
2067           continue;
2068         }
2069         if (point->GetCPoint().GetMax()<5) continue;
2070         Float_t angley = point->GetTPoint().GetAngleY();
2071         Float_t anglez = point->GetTPoint().GetAngleZ();
2072         Int_t   type   = point->GetCPoint().GetType();
2073         Float_t rsigmay =  point->GetCPoint().GetSigmaY();
2074         Float_t rsigmaz =  point->GetCPoint().GetSigmaZ();
2075         Float_t rsigma = TMath::Sqrt(rsigmay*rsigmaz);
2076
2077         Float_t ampc   = 0;     // normalization to the number of electrons
2078         if (i>64){
2079           ampc = 1.*point->GetCPoint().GetMax();
2080           //ampc = 1.*point->GetCPoint().GetQ();          
2081           //      AliTPCClusterPoint & p = point->GetCPoint();
2082           //      Float_t dy = TMath::Abs(Int_t( TMath::Abs(p.GetY()/0.6)) - TMath::Abs(p.GetY()/0.6)+0.5);
2083           // Float_t iz =  (250.0-TMath::Abs(p.GetZ())+0.11)/0.566;
2084           //Float_t dz = 
2085           //  TMath::Abs( Int_t(iz) - iz + 0.5);
2086           //ampc *= 1.15*(1-0.3*dy);
2087           //ampc *= 1.15*(1-0.3*dz);
2088           //      Float_t zfactor = (1.05-0.0004*TMath::Abs(point->GetCPoint().GetZ()));
2089           //ampc               *=zfactor; 
2090         }
2091         else{ 
2092           ampc = 1.0*point->GetCPoint().GetMax(); 
2093           //ampc = 1.0*point->GetCPoint().GetQ(); 
2094           //AliTPCClusterPoint & p = point->GetCPoint();
2095           // Float_t dy = TMath::Abs(Int_t( TMath::Abs(p.GetY()/0.4)) - TMath::Abs(p.GetY()/0.4)+0.5);
2096           //Float_t iz =  (250.0-TMath::Abs(p.GetZ())+0.11)/0.566;
2097           //Float_t dz = 
2098           //  TMath::Abs( Int_t(iz) - iz + 0.5);
2099
2100           //ampc *= 1.15*(1-0.3*dy);
2101           //ampc *= 1.15*(1-0.3*dz);
2102           //    Float_t zfactor = (1.02-0.000*TMath::Abs(point->GetCPoint().GetZ()));
2103           //ampc               *=zfactor; 
2104
2105         }
2106         ampc *= 2.0;     // put mean value to channel 50
2107         //ampc *= 0.58;     // put mean value to channel 50
2108         Float_t w      =  1.;
2109         //      if (type>0)  w =  1./(type/2.-0.5); 
2110         Float_t z = TMath::Abs(point->GetCPoint().GetZ());
2111         if (i<64) {
2112           ampc /= 0.6;
2113           //ampc /= (1+0.0008*z);
2114         } else
2115           if (i>128){
2116             ampc /=1.5;
2117             //ampc /= (1+0.0008*z);
2118           }else{
2119             //ampc /= (1+0.0008*z);
2120           }
2121         
2122         if (type<0) {  //amp at the border - lower weight
2123           // w*= 2.;
2124           
2125           continue;
2126         }
2127         if (rsigma>1.5) ampc/=1.3;  // if big backround
2128         amp[nc[of]]        = ampc;
2129         angular[nc[of]]    = TMath::Sqrt(1.+angley*angley+anglez*anglez);
2130         weight[nc[of]]     = w;
2131         nc[of]++;
2132       }
2133     
2134     TMath::Sort(nc[of],amp,index,kFALSE);
2135     Float_t sumamp=0;
2136     Float_t sumamp2=0;
2137     Float_t sumw=0;
2138     //meanlog = amp[index[Int_t(nc[of]*0.33)]];
2139     meanlog = 200;
2140     for (Int_t i=int(nc[of]*low+0.5);i<int(nc[of]*up+0.5);i++){
2141       Float_t ampl      = amp[index[i]]/angular[index[i]];
2142       ampl              = meanlog*TMath::Log(1.+ampl/meanlog);
2143       //
2144       sumw    += weight[index[i]]; 
2145       sumamp  += weight[index[i]]*ampl;
2146       sumamp2 += weight[index[i]]*ampl*ampl;
2147       norm[of]    += angular[index[i]]*weight[index[i]];
2148     }
2149     if (sumw<1){ 
2150       SetdEdx(0);  
2151     }
2152     else {
2153       norm[of] /= sumw;
2154       mean[of]  = sumamp/sumw;
2155       sigma[of] = sumamp2/sumw-mean[of]*mean[of];
2156       if (sigma[of]>0.1) 
2157         sigma[of] = TMath::Sqrt(sigma[of]);
2158       else
2159         sigma[of] = 1000;
2160       
2161     mean[of] = (TMath::Exp(mean[of]/meanlog)-1)*meanlog;
2162     //mean  *=(1-0.02*(sigma/(mean*0.17)-1.));
2163     //mean *=(1-0.1*(norm-1.));
2164     }
2165   }
2166
2167   Float_t dedx =0;
2168   fSdEdx =0;
2169   fMAngular =0;
2170   //  mean[0]*= (1-0.05*(sigma[0]/(0.01+mean[1]*0.18)-1));
2171   //  mean[1]*= (1-0.05*(sigma[1]/(0.01+mean[0]*0.18)-1));
2172
2173   
2174   //  dedx = (mean[0]* TMath::Sqrt((1.+nc[0]))+ mean[1]* TMath::Sqrt((1.+nc[1])) )/ 
2175   //  (  TMath::Sqrt((1.+nc[0]))+TMath::Sqrt((1.+nc[1])));
2176
2177   Int_t norm2 = 0;
2178   Int_t norm3 = 0;
2179   for (Int_t i =0;i<4;i++){
2180     if (nc[i]>2&&nc[i]<1000){
2181       dedx      += mean[i] *nc[i];
2182       fSdEdx    += sigma[i]*(nc[i]-2);
2183       fMAngular += norm[i] *nc[i];    
2184       norm2     += nc[i];
2185       norm3     += nc[i]-2;
2186     }
2187     fDEDX[i]  = mean[i];             
2188     fSDEDX[i] = sigma[i];            
2189     fNCDEDX[i]= nc[i]; 
2190   }
2191
2192   if (norm3>0){
2193     dedx   /=norm2;
2194     fSdEdx /=norm3;
2195     fMAngular/=norm2;
2196   }
2197   else{
2198     SetdEdx(0);
2199     return;
2200   }
2201   //  Float_t dedx1 =dedx;
2202   
2203   dedx =0;
2204   for (Int_t i =0;i<4;i++){
2205     if (nc[i]>2&&nc[i]<1000){
2206       mean[i]   = mean[i]*(1-0.12*(sigma[i]/(fSdEdx)-1.));
2207       dedx      += mean[i] *nc[i];
2208     }
2209     fDEDX[i]  = mean[i];                
2210   }
2211   dedx /= norm2;
2212   
2213
2214   
2215   SetdEdx(dedx);
2216     
2217   //mi deDX
2218
2219
2220
2221   //Very rough PID
2222   Double_t p=TMath::Sqrt((1.+ GetTgl()*GetTgl())/(Get1Pt()*Get1Pt()));
2223
2224   if (p<0.6) {
2225     if (dedx < 39.+ 12./(p+0.25)/(p+0.25)) { SetMass(0.13957); return;}
2226     if (dedx < 39.+ 12./p/p) { SetMass(0.49368); return;}
2227     SetMass(0.93827); return;
2228   }
2229
2230   if (p<1.2) {
2231     if (dedx < 39.+ 12./(p+0.25)/(p+0.25)) { SetMass(0.13957); return;}
2232     SetMass(0.93827); return;
2233   }
2234
2235   SetMass(0.13957); return;
2236
2237 }
2238
2239