Adding protection for null clusters in assignment operator
[u/mrichter/AliRoot.git] / TPC / Rec / AliTPCseed.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 //           Implementation of the TPC seed class
22 //        This class is used by the AliTPCtracker class
23 //      Origin: Marian Ivanov, CERN, Marian.Ivanov@cern.ch
24 //-----------------------------------------------------------------
25 #include <TVectorF.h>
26 #include "TClonesArray.h"
27 #include "TGraphErrors.h"
28 #include "AliTPCseed.h"
29 #include "AliTPCReconstructor.h"
30 #include "AliTPCClusterParam.h"
31 #include "AliTPCCalPad.h"
32 #include "AliTPCCalROC.h"
33 #include "AliTPCcalibDB.h"
34 #include "AliTPCParam.h"
35 #include "AliMathBase.h"
36 #include "AliTPCTransform.h"
37 #include "AliSplineFit.h"
38 #include "AliCDBManager.h"
39 #include "AliTPCcalibDButil.h"
40 #include <AliCTPTimeParams.h>
41
42
43 ClassImp(AliTPCseed)
44
45
46
47 AliTPCseed::AliTPCseed():
48   AliTPCtrack(),
49   fEsd(0x0),
50   fClusterOwner(kFALSE),
51   fRow(0),
52   fSector(-1),
53   fRelativeSector(-1),
54   fCurrentSigmaY2(1e10),
55   fCurrentSigmaZ2(1e10),
56   fCMeanSigmaY2p30(-1.),   //! current mean sigma Y2 - mean30%
57   fCMeanSigmaZ2p30(-1.),   //! current mean sigma Z2 - mean30%
58   fCMeanSigmaY2p30R(-1.),   //! current mean sigma Y2 - mean2%
59   fCMeanSigmaZ2p30R(-1.),   //! current mean sigma Z2 - mean2%
60   //
61   fErrorY2(1e10),
62   fErrorZ2(1e10),
63   fCurrentCluster(0x0),
64   fCurrentClusterIndex1(-1),
65   fInDead(kFALSE),
66   fIsSeeding(kFALSE),
67   fNoCluster(0),
68   fSort(0),
69   fBSigned(kFALSE),
70   fSeedType(0),
71   fSeed1(-1),
72   fSeed2(-1),
73   fMAngular(0),
74   fCircular(0),
75   fPoolID(-1)
76 {
77   //
78   for (Int_t i=0;i<160;i++) SetClusterIndex2(i,-3);
79   for (Int_t i=0;i<160;i++) fClusterPointer[i]=0;
80   for (Int_t i=0;i<3;i++)   fKinkIndexes[i]=0;
81   for (Int_t i=0;i<AliPID::kSPECIES;i++)   fTPCr[i]=0.2;
82   for (Int_t i=0;i<4;i++) {
83     fDEDX[i] = 0.;
84     fSDEDX[i] = 1e10;
85     fNCDEDX[i] = 0;
86     fNCDEDXInclThres[i] = 0;
87   }
88   for (Int_t i=0;i<9;i++) fDEDX[i] = 0;
89   for (Int_t i=0;i<12;i++) fOverlapLabels[i] = -1;
90 }
91
92 AliTPCseed::AliTPCseed(const AliTPCseed &s, Bool_t clusterOwner):
93   AliTPCtrack(s),
94   fEsd(0x0),
95   fClusterOwner(clusterOwner),
96   fRow(0),
97   fSector(-1),
98   fRelativeSector(-1),
99   fCurrentSigmaY2(-1),
100   fCurrentSigmaZ2(-1),
101   fCMeanSigmaY2p30(-1.),   //! current mean sigma Y2 - mean30%
102   fCMeanSigmaZ2p30(-1.),   //! current mean sigma Z2 - mean30%
103   fCMeanSigmaY2p30R(-1.),   //! current mean sigma Y2 - mean2%
104   fCMeanSigmaZ2p30R(-1.),   //! current mean sigma Z2 - mean2%
105   fErrorY2(1e10),
106   fErrorZ2(1e10),
107   fCurrentCluster(0x0),
108   fCurrentClusterIndex1(-1),
109   fInDead(kFALSE),
110   fIsSeeding(kFALSE),
111   fNoCluster(0),
112   fSort(0),
113   fBSigned(kFALSE),
114   fSeedType(0),
115   fSeed1(-1),
116   fSeed2(-1),
117   fMAngular(0),
118   fCircular(0),
119   fPoolID(-1)
120 {
121   //---------------------
122   // dummy copy constructor
123   //-------------------------
124   for (Int_t i=0;i<160;i++) {
125     fClusterPointer[i]=0;
126     if (fClusterOwner){
127       if (s.fClusterPointer[i])
128         fClusterPointer[i] = new AliTPCclusterMI(*(s.fClusterPointer[i]));
129     }else{
130       fClusterPointer[i] = s.fClusterPointer[i];
131     }
132     fTrackPoints[i] = s.fTrackPoints[i];
133   }
134   for (Int_t i=0;i<160;i++) fIndex[i] = s.fIndex[i];
135   for (Int_t i=0;i<AliPID::kSPECIES;i++)   fTPCr[i]=s.fTPCr[i];
136   for (Int_t i=0;i<4;i++) {
137     fDEDX[i] = s.fDEDX[i];
138     fSDEDX[i] = s.fSDEDX[i];
139     fNCDEDX[i] = s.fNCDEDX[i];
140     fNCDEDXInclThres[i] = s.fNCDEDXInclThres[i];
141   }
142   for (Int_t i=0;i<9;i++) fDEDX[i] = 0;
143
144   for (Int_t i=0;i<12;i++) fOverlapLabels[i] = s.fOverlapLabels[i];
145
146 }
147
148
149 AliTPCseed::AliTPCseed(const AliTPCtrack &t):
150   AliTPCtrack(t),
151   fEsd(0x0),
152   fClusterOwner(kFALSE),
153   fRow(0),
154   fSector(-1),
155   fRelativeSector(-1),
156   fCurrentSigmaY2(-1),
157   fCurrentSigmaZ2(-1),
158   fCMeanSigmaY2p30(-1.),   //! current mean sigma Y2 - mean30%
159   fCMeanSigmaZ2p30(-1.),   //! current mean sigma Z2 - mean30%
160   fCMeanSigmaY2p30R(-1.),   //! current mean sigma Y2 - mean2%
161   fCMeanSigmaZ2p30R(-1.),   //! current mean sigma Z2 - mean2%
162   fErrorY2(1e10),
163   fErrorZ2(1e10),
164   fCurrentCluster(0x0),
165   fCurrentClusterIndex1(-1),
166   fInDead(kFALSE),
167   fIsSeeding(kFALSE),
168   fNoCluster(0),
169   fSort(0),
170   fBSigned(kFALSE),
171   fSeedType(0),
172   fSeed1(-1),
173   fSeed2(-1),
174   fMAngular(0),
175   fCircular(0),
176   fPoolID(-1)
177 {
178   //
179   // Constructor from AliTPCtrack
180   //
181   fFirstPoint =0;
182   for (Int_t i=0;i<5;i++)   fTPCr[i]=0.2;
183   for (Int_t i=0;i<160;i++) {
184     fClusterPointer[i] = 0;
185     Int_t index = t.GetClusterIndex(i);
186     if (index>=-1){ 
187       SetClusterIndex2(i,index);
188     }
189     else{
190       SetClusterIndex2(i,-3); 
191     }    
192   }
193   for (Int_t i=0;i<4;i++) {
194     fDEDX[i] = 0.;
195     fSDEDX[i] = 1e10;
196     fNCDEDX[i] = 0;
197     fNCDEDXInclThres[i] = 0;
198   }
199   for (Int_t i=0;i<9;i++) fDEDX[i] = fDEDX[i];
200
201   for (Int_t i=0;i<12;i++) fOverlapLabels[i] = -1;
202 }
203
204 AliTPCseed::AliTPCseed(Double_t xr, Double_t alpha, const Double_t xx[5],
205                        const Double_t cc[15], Int_t index):      
206   AliTPCtrack(xr, alpha, xx, cc, index),
207   fEsd(0x0),
208   fClusterOwner(kFALSE),
209   fRow(0),
210   fSector(-1),
211   fRelativeSector(-1),
212   fCurrentSigmaY2(-1),
213   fCurrentSigmaZ2(-1),
214   fCMeanSigmaY2p30(-1.),   //! current mean sigma Y2 - mean30%
215   fCMeanSigmaZ2p30(-1.),   //! current mean sigma Z2 - mean30%
216   fCMeanSigmaY2p30R(-1.),   //! current mean sigma Y2 - mean2%
217   fCMeanSigmaZ2p30R(-1.),   //! current mean sigma Z2 - mean2%
218   fErrorY2(1e10),
219   fErrorZ2(1e10),
220   fCurrentCluster(0x0),
221   fCurrentClusterIndex1(-1),
222   fInDead(kFALSE),
223   fIsSeeding(kFALSE),
224   fNoCluster(0),
225   fSort(0),
226   fBSigned(kFALSE),
227   fSeedType(0),
228   fSeed1(-1),
229   fSeed2(-1),
230   fMAngular(0),
231   fCircular(0),
232   fPoolID(-1)
233 {
234   //
235   // Constructor
236   //
237   fFirstPoint =0;
238   for (Int_t i=0;i<160;i++) SetClusterIndex2(i,-3);
239   for (Int_t i=0;i<160;i++) fClusterPointer[i]=0;
240   for (Int_t i=0;i<5;i++)   fTPCr[i]=0.2;
241   for (Int_t i=0;i<4;i++) {
242     fDEDX[i] = 0.;
243     fSDEDX[i] = 1e10;
244     fNCDEDX[i] = 0;
245     fNCDEDXInclThres[i] = 0;
246   }
247   for (Int_t i=0;i<9;i++) fDEDX[i] = 0;
248
249   for (Int_t i=0;i<12;i++) fOverlapLabels[i] = -1;
250 }
251
252 AliTPCseed::~AliTPCseed(){
253   //
254   // destructor
255   fNoCluster =0;
256   if (fClusterOwner){
257     for (Int_t icluster=0; icluster<160; icluster++){
258       delete fClusterPointer[icluster];
259     }
260   }
261
262 }
263 //_________________________________________________
264 AliTPCseed & AliTPCseed::operator=(const AliTPCseed &param)
265 {
266   //
267   // assignment operator 
268   // don't touch pool ID
269   //
270   if(this!=&param){
271     AliTPCtrack::operator=(param);
272     fEsd =param.fEsd; 
273     fClusterOwner = param.fClusterOwner;
274     if (!fClusterOwner) for(Int_t i = 0;i<160;++i)fClusterPointer[i] = param.fClusterPointer[i];
275     else                for(Int_t i = 0;i<160;++i) {
276         delete fClusterPointer[i];
277         if (param.fClusterPointer[i]) { 
278           fClusterPointer[i] = new AliTPCclusterMI(*(param.fClusterPointer[i]));
279         }
280         else {
281           fClusterPointer[i] = 0x0;
282         }
283       }
284     // leave out fPoint, they are also not copied in the copy ctor...
285     // but deleted in the dtor... strange...
286     fRow            = param.fRow;
287     fSector         = param.fSector;
288     fRelativeSector = param.fRelativeSector;
289     fCurrentSigmaY2 = param.fCurrentSigmaY2;
290     fCurrentSigmaZ2 = param.fCurrentSigmaZ2;
291     fErrorY2        = param.fErrorY2;
292     fErrorZ2        = param.fErrorZ2;
293     fCurrentCluster = param.fCurrentCluster; // this is not allocated by AliTPCSeed
294     fCurrentClusterIndex1 = param.fCurrentClusterIndex1; 
295     fInDead         = param.fInDead;
296     fIsSeeding      = param.fIsSeeding;
297     fNoCluster      = param.fNoCluster;
298     fSort           = param.fSort;
299     fBSigned        = param.fBSigned;
300     for(Int_t i = 0;i<4;++i){
301       fDEDX[i]   = param.fDEDX[i];
302       fSDEDX[i]  = param.fSDEDX[i];
303       fNCDEDX[i] = param.fNCDEDX[i];
304       fNCDEDXInclThres[i] = param.fNCDEDXInclThres[i];
305     }
306     for (Int_t i=0;i<9;i++) fDEDX[i] = 0;
307
308     for(Int_t i = 0;i<AliPID::kSPECIES;++i)fTPCr[i] = param.fTPCr[i];
309     
310     fSeedType = param.fSeedType;
311     fSeed1    = param.fSeed1;
312     fSeed2    = param.fSeed2;
313     for(Int_t i = 0;i<12;++i)fOverlapLabels[i] = param.fOverlapLabels[i];
314     fMAngular = param.fMAngular;
315     fCircular = param.fCircular;
316     for(int i = 0;i<160;++i)fTrackPoints[i] =  param.fTrackPoints[i];
317   }
318   return (*this);
319 }
320 //____________________________________________________
321 AliTPCTrackerPoint * AliTPCseed::GetTrackPoint(Int_t i)
322 {
323   //
324   // 
325   return &fTrackPoints[i];
326 }
327
328
329
330 Double_t AliTPCseed::GetDensityFirst(Int_t n)
331 {
332   //
333   //
334   // return cluster for n rows bellow first point
335   Int_t nfoundable = 1;
336   Int_t nfound      = 1;
337   for (Int_t i=fLastPoint-1;i>0&&nfoundable<n; i--){
338     Int_t index = GetClusterIndex2(i);
339     if (index!=-1) nfoundable++;
340     if (index>0) nfound++;
341   }
342   if (nfoundable<n) return 0;
343   return Double_t(nfound)/Double_t(nfoundable);
344
345 }
346
347
348 void AliTPCseed::GetClusterStatistic(Int_t first, Int_t last, Int_t &found, Int_t &foundable, Int_t &shared, Bool_t plus2)
349 {
350   // get cluster stat.  on given region
351   //
352   found       = 0;
353   foundable   = 0;
354   shared      =0;
355   for (Int_t i=first;i<last; i++){
356     Int_t index = GetClusterIndex2(i);
357     if (index!=-1) foundable++;
358     if (index&0x8000) continue;
359     if (fClusterPointer[i]) {
360       found++;
361     }
362     else 
363       continue;
364
365     if (fClusterPointer[i]->IsUsed(10)) {
366       shared++;
367       continue;
368     }
369     if (!plus2) continue; //take also neighborhoud
370     //
371     if ( (i>0) && fClusterPointer[i-1]){
372       if (fClusterPointer[i-1]->IsUsed(10)) {
373         shared++;
374         continue;
375       }
376     }
377     if ( fClusterPointer[i+1]){
378       if (fClusterPointer[i+1]->IsUsed(10)) {
379         shared++;
380         continue;
381       }
382     }
383     
384   }
385   //if (shared>found){
386     //Error("AliTPCseed::GetClusterStatistic","problem\n");
387   //}
388 }
389
390
391
392
393
394 void AliTPCseed::Reset(Bool_t all)
395 {
396   //
397   //
398   SetNumberOfClusters(0);
399   fNFoundable = 0;
400   SetChi2(0);
401   ResetCovariance(10.);
402   /*
403   if (fTrackPoints){
404     for (Int_t i=0;i<8;i++){
405       delete [] fTrackPoints[i];
406     }
407     delete fTrackPoints;
408     fTrackPoints =0;
409   }
410   */
411
412   if (all){   
413     for (Int_t i=200;i--;) SetClusterIndex2(i,-3);
414     if (!fClusterOwner) for (Int_t i=160;i--;) fClusterPointer[i]=0;
415     else                for (Int_t i=160;i--;) {delete fClusterPointer[i]; fClusterPointer[i]=0;}
416   }
417
418 }
419
420
421 void AliTPCseed::Modify(Double_t factor)
422 {
423
424   //------------------------------------------------------------------
425   //This function makes a track forget its history :)  
426   //------------------------------------------------------------------
427   if (factor<=0) {
428     ResetCovariance(10.);
429     return;
430   }
431   ResetCovariance(factor);
432
433   SetNumberOfClusters(0);
434   fNFoundable =0;
435   SetChi2(0);
436   fRemoval = 0;
437   fCurrentSigmaY2 = 0.000005;
438   fCurrentSigmaZ2 = 0.000005;
439   fNoCluster     = 0;
440   //fFirstPoint = 160;
441   //fLastPoint  = 0;
442 }
443
444
445
446
447 Int_t  AliTPCseed::GetProlongation(Double_t xk, Double_t &y, Double_t & z) const
448 {
449   //-----------------------------------------------------------------
450   // This function find proloncation of a track to a reference plane x=xk.
451   // doesn't change internal state of the track
452   //-----------------------------------------------------------------
453   
454   Double_t x1=GetX(), x2=x1+(xk-x1), dx=x2-x1;
455
456   if (TMath::Abs(GetSnp()+GetC()*dx) >= AliTPCReconstructor::GetMaxSnpTrack()) {   
457     return 0;
458   }
459
460   //  Double_t y1=fP0, z1=fP1;
461   Double_t c1=GetSnp(), r1=sqrt((1.-c1)*(1.+c1));
462   Double_t c2=c1 + GetC()*dx, r2=sqrt((1.-c2)*(1.+c2));
463   
464   y = GetY();
465   z = GetZ();
466   //y += dx*(c1+c2)/(r1+r2);
467   //z += dx*(c1+c2)/(c1*r2 + c2*r1)*fP3;
468   
469   Double_t dy = dx*(c1+c2)/(r1+r2);
470   Double_t dz = 0;
471   //
472   Double_t delta = GetC()*dx*(c1+c2)/(c1*r2 + c2*r1);
473   /*
474   if (TMath::Abs(delta)>0.0001){
475     dz = fP3*TMath::ASin(delta)/fP4;
476   }else{
477     dz = dx*fP3*(c1+c2)/(c1*r2 + c2*r1);
478   }
479   */
480   //  dz =  fP3*AliTPCFastMath::FastAsin(delta)/fP4;
481   dz =  GetTgl()*TMath::ASin(delta)/GetC();
482   //
483   y+=dy;
484   z+=dz;
485   
486
487   return 1;  
488 }
489
490
491 //_____________________________________________________________________________
492 Double_t AliTPCseed::GetPredictedChi2(const AliCluster *c) const 
493 {
494   //-----------------------------------------------------------------
495   // This function calculates a predicted chi2 increment.
496   //-----------------------------------------------------------------
497   Double_t p[2]={c->GetY(), c->GetZ()};
498   Double_t cov[3]={fErrorY2, 0., fErrorZ2};
499
500   Float_t dx = ((AliTPCclusterMI*)c)->GetX()-GetX();
501   if (TMath::Abs(dx)>0){
502     Float_t ty = TMath::Tan(TMath::ASin(GetSnp()));
503     Float_t dy = dx*ty;
504     Float_t dz = dx*TMath::Sqrt(1.+ty*ty)*GetTgl();
505     p[0] = c->GetY()-dy;  
506     p[1] = c->GetZ()-dz;  
507   }
508   return AliExternalTrackParam::GetPredictedChi2(p,cov);
509 }
510
511 //_________________________________________________________________________________________
512
513
514 Int_t AliTPCseed::Compare(const TObject *o) const {
515   //-----------------------------------------------------------------
516   // This function compares tracks according to the sector - for given sector according z
517   //-----------------------------------------------------------------
518   AliTPCseed *t=(AliTPCseed*)o;
519
520   if (fSort == 0){
521     if (t->fRelativeSector>fRelativeSector) return -1;
522     if (t->fRelativeSector<fRelativeSector) return 1;
523     Double_t z2 = t->GetZ();
524     Double_t z1 = GetZ();
525     if (z2>z1) return 1;
526     if (z2<z1) return -1;
527     return 0;
528   }
529   else {
530     Float_t f2 =1;
531     f2 = 1-20*TMath::Sqrt(t->GetSigma1Pt2())/(t->OneOverPt()+0.0066);
532     if (t->fBConstrain) f2=1.2;
533
534     Float_t f1 =1;
535     f1 = 1-20*TMath::Sqrt(GetSigma1Pt2())/(OneOverPt()+0.0066);
536
537     if (fBConstrain)   f1=1.2;
538  
539     if (t->GetNumberOfClusters()*f2 <GetNumberOfClusters()*f1) return -1;
540     else return +1;
541   }
542 }
543
544
545
546
547 //_____________________________________________________________________________
548 Bool_t AliTPCseed::Update(const AliCluster *c, Double_t chisq, Int_t index)
549 {
550   //-----------------------------------------------------------------
551   // This function associates a cluster with this track.
552   //-----------------------------------------------------------------
553   Int_t n=GetNumberOfClusters();
554   Int_t idx=GetClusterIndex(n);    // save the current cluster index
555
556   AliTPCclusterMI cl(*(AliTPCclusterMI*)c);  cl.SetSigmaY2(fErrorY2); cl.SetSigmaZ2(fErrorZ2);
557
558   AliTPCClusterParam * parcl = AliTPCcalibDB::Instance()->GetClusterParam();
559   
560   Float_t ty = TMath::Tan(TMath::ASin(GetSnp()));
561   
562   if(  parcl ){
563     Int_t padSize = 0;                          // short pads
564     if (cl.GetDetector() >= 36) {
565       padSize = 1;                              // medium pads 
566       if (cl.GetRow() > 63) padSize = 2; // long pads
567     }
568     Float_t waveCorr = parcl->GetWaveCorrection( padSize, cl.GetZ(), cl.GetMax(),cl.GetPad(), ty );
569     cl.SetY( cl.GetY() - waveCorr ); 
570   }
571
572   Float_t dx = ((AliTPCclusterMI*)c)->GetX()-GetX();
573   if (TMath::Abs(dx)>0){
574     Float_t dy = dx*ty;
575     Float_t dz = dx*TMath::Sqrt(1.+ty*ty)*GetTgl();
576     cl.SetY(cl.GetY()-dy);  
577     cl.SetZ(cl.GetZ()-dz);  
578   }  
579  
580
581   if (!AliTPCtrack::Update(&cl,chisq,index)) return kFALSE;
582   
583   if (fCMeanSigmaY2p30<0){
584     fCMeanSigmaY2p30= c->GetSigmaY2();   //! current mean sigma Y2 - mean30%
585     fCMeanSigmaZ2p30= c->GetSigmaZ2();   //! current mean sigma Z2 - mean30%    
586     fCMeanSigmaY2p30R = 1;   //! current mean sigma Y2 - mean5%
587     fCMeanSigmaZ2p30R = 1;   //! current mean sigma Z2 - mean5%
588   }
589   //
590   fCMeanSigmaY2p30= 0.70*fCMeanSigmaY2p30 +0.30*c->GetSigmaY2();   
591   fCMeanSigmaZ2p30= 0.70*fCMeanSigmaZ2p30 +0.30*c->GetSigmaZ2();  
592   if (fCurrentSigmaY2>0){
593     fCMeanSigmaY2p30R = 0.7*fCMeanSigmaY2p30R  +0.3*c->GetSigmaY2()/fCurrentSigmaY2;  
594     fCMeanSigmaZ2p30R = 0.7*fCMeanSigmaZ2p30R  +0.3*c->GetSigmaZ2()/fCurrentSigmaZ2;   
595   }
596
597
598   SetClusterIndex(n,idx);          // restore the current cluster index
599   return kTRUE;
600 }
601
602
603
604 //_____________________________________________________________________________
605 Float_t AliTPCseed::CookdEdx(Double_t low, Double_t up,Int_t i1, Int_t i2, Bool_t /* onlyused */) {
606   //-----------------------------------------------------------------
607   // This funtion calculates dE/dX within the "low" and "up" cuts.
608   //-----------------------------------------------------------------
609   // CookdEdxAnalytical(Double_t low, Double_t up, Int_t type, Int_t i1, Int_t i2, Int_t returnVal)
610   AliTPCParam *param = AliTPCcalibDB::Instance()->GetParameters();
611   
612   Int_t row0 = param->GetNRowLow();
613   Int_t row1 = row0+param->GetNRowUp1();
614   Int_t row2 = row1+param->GetNRowUp2();
615   const AliTPCRecoParam * recoParam = AliTPCcalibDB::Instance()->GetTransform()->GetCurrentRecoParam();
616   Int_t useTot = 0;
617   if (recoParam) useTot = (recoParam->GetUseTotCharge())? 0:1;
618   //
619   //
620   TVectorF i1i2;
621   TVectorF  irocTot;
622   TVectorF oroc1Tot;
623   TVectorF oroc2Tot;
624   TVectorF forocTot;
625   //
626   TVectorF  irocMax;
627   TVectorF oroc1Max;
628   TVectorF oroc2Max;
629   TVectorF forocMax;
630
631   CookdEdxAnalytical(low,up,useTot ,i1  ,i2,   0, 2, 0, &i1i2);
632   //
633   CookdEdxAnalytical(low,up,kTRUE ,0   ,row0, 0, 2, 0, &irocTot);
634   CookdEdxAnalytical(low,up,kTRUE ,row0,row1, 0, 2, 0, &oroc1Tot);
635   CookdEdxAnalytical(low,up,kTRUE ,row1,row2, 0, 2, 0, &oroc2Tot);
636   CookdEdxAnalytical(low,up,kTRUE ,row0,row2, 0, 2, 0, &forocTot); // full OROC truncated mean
637   //
638   CookdEdxAnalytical(low,up,kFALSE ,0   ,row0, 0, 2, 0, &irocMax);
639   CookdEdxAnalytical(low,up,kFALSE ,row0,row1, 0, 2, 0, &oroc1Max);
640   CookdEdxAnalytical(low,up,kFALSE ,row1,row2, 0, 2, 0, &oroc2Max);
641   CookdEdxAnalytical(low,up,kFALSE ,row0,row2, 0, 2, 0, &forocMax); // full OROC truncated mean
642
643   fDEDX[0]      = i1i2(0);
644   //
645   fDEDX[1]      =  irocTot(0);
646   fDEDX[2]      = oroc1Tot(0);
647   fDEDX[3]      = oroc2Tot(0);
648   fDEDX[4]      = forocTot(0); // full OROC truncated mean
649   fDEDX[5]      =  irocMax(0);
650   fDEDX[6]      = oroc1Max(0);
651   fDEDX[7]      = oroc2Max(0);
652   fDEDX[8]      = forocMax(0); // full OROC truncated mean
653   //
654   fSDEDX[0]     = i1i2(1);
655   fSDEDX[1]     =  irocTot(1);
656   fSDEDX[2]     = oroc1Tot(1);
657   fSDEDX[3]     = oroc2Tot(1);
658   //
659   fNCDEDX[0]    = TMath::Nint(i1i2(2));
660
661   fNCDEDX[1]    = TMath::Nint( irocTot(2));
662   fNCDEDX[2]    = TMath::Nint(oroc1Tot(2));
663   fNCDEDX[3]    = TMath::Nint(oroc2Tot(2));
664   //
665   fNCDEDXInclThres[0]    = TMath::Nint(i1i2(2)+i1i2(9));
666   fNCDEDXInclThres[1]    = TMath::Nint( irocTot(2)+ irocTot(9));
667   fNCDEDXInclThres[2]    = TMath::Nint(oroc1Tot(2)+oroc1Tot(9));
668   fNCDEDXInclThres[3]    = TMath::Nint(oroc2Tot(2)+oroc2Tot(9));
669   //
670   SetdEdx(fDEDX[0]);
671   return fDEDX[0];
672
673 //  return CookdEdxNorm(low,up,0,i1,i2,1,0,2);
674
675
676 //   Float_t amp[200];
677 //   Float_t angular[200];
678 //   Float_t weight[200];
679 //   Int_t index[200];
680 //   //Int_t nc = 0;
681 //   Float_t meanlog = 100.;
682   
683 //   Float_t mean[4]  = {0,0,0,0};
684 //   Float_t sigma[4] = {1000,1000,1000,1000};
685 //   Int_t nc[4]      = {0,0,0,0};
686 //   Float_t norm[4]    = {1000,1000,1000,1000};
687 //   //
688 //   //
689 //   fNShared =0;
690
691 //   Float_t gainGG = 1;
692 //   if (AliTPCcalibDB::Instance()->GetParameters()){
693 //     gainGG= AliTPCcalibDB::Instance()->GetParameters()->GetGasGain()/20000.;  //relative gas gain
694 //   }
695
696
697 //   for (Int_t of =0; of<4; of++){    
698 //     for (Int_t i=of+i1;i<i2;i+=4)
699 //       {
700 //      Int_t clindex = fIndex[i];
701 //      if (clindex<0||clindex&0x8000) continue;
702
703 //      //AliTPCTrackPoint * point = (AliTPCTrackPoint *) arr.At(i);
704 //      AliTPCTrackerPoint * point = GetTrackPoint(i);
705 //      //AliTPCTrackerPoint * pointm = GetTrackPoint(i-1);
706 //      //AliTPCTrackerPoint * pointp = 0;
707 //      //if (i<159) pointp = GetTrackPoint(i+1);
708
709 //      if (point==0) continue;
710 //      AliTPCclusterMI * cl = fClusterPointer[i];
711 //      if (cl==0) continue;    
712 //      if (onlyused && (!cl->IsUsed(10))) continue;
713 //      if (cl->IsUsed(11)) {
714 //        fNShared++;
715 //        continue;
716 //      }
717 //      Int_t   type   = cl->GetType();
718 //      //if (point->fIsShared){
719 //      //  fNShared++;
720 //      //  continue;
721 //      //}
722 //      //if (pointm) 
723 //      //  if (pointm->fIsShared) continue;
724 //      //if (pointp) 
725 //      //  if (pointp->fIsShared) continue;
726
727 //      if (type<0) continue;
728 //      //if (type>10) continue;       
729 //      //if (point->GetErrY()==0) continue;
730 //      //if (point->GetErrZ()==0) continue;
731
732 //      //Float_t ddy = (point->GetY()-cl->GetY())/point->GetErrY();
733 //      //Float_t ddz = (point->GetZ()-cl->GetZ())/point->GetErrZ();
734 //      //if ((ddy*ddy+ddz*ddz)>10) continue; 
735
736
737 //      //      if (point->GetCPoint().GetMax()<5) continue;
738 //      if (cl->GetMax()<5) continue;
739 //      Float_t angley = point->GetAngleY();
740 //      Float_t anglez = point->GetAngleZ();
741
742 //      Float_t rsigmay2 =  point->GetSigmaY();
743 //      Float_t rsigmaz2 =  point->GetSigmaZ();
744 //      /*
745 //      Float_t ns = 1.;
746 //      if (pointm){
747 //        rsigmay +=  pointm->GetTPoint().GetSigmaY();
748 //        rsigmaz +=  pointm->GetTPoint().GetSigmaZ();
749 //        ns+=1.;
750 //      }
751 //      if (pointp){
752 //        rsigmay +=  pointp->GetTPoint().GetSigmaY();
753 //        rsigmaz +=  pointp->GetTPoint().GetSigmaZ();
754 //        ns+=1.;
755 //      }
756 //      rsigmay/=ns;
757 //      rsigmaz/=ns;
758 //      */
759
760 //      Float_t rsigma = TMath::Sqrt(rsigmay2*rsigmaz2);
761
762 //      Float_t ampc   = 0;     // normalization to the number of electrons
763 //      if (i>64){
764 //        //      ampc = 1.*point->GetCPoint().GetMax();
765 //        ampc = 1.*cl->GetMax();
766 //        //ampc = 1.*point->GetCPoint().GetQ();          
767 //        //      AliTPCClusterPoint & p = point->GetCPoint();
768 //        //      Float_t dy = TMath::Abs(Int_t( TMath::Abs(p.GetY()/0.6)) - TMath::Abs(p.GetY()/0.6)+0.5);
769 //        // Float_t iz =  (250.0-TMath::Abs(p.GetZ())+0.11)/0.566;
770 //        //Float_t dz = 
771 //        //  TMath::Abs( Int_t(iz) - iz + 0.5);
772 //        //ampc *= 1.15*(1-0.3*dy);
773 //        //ampc *= 1.15*(1-0.3*dz);
774 //        //      Float_t zfactor = (AliTPCReconstructor::GetCtgRange()-0.0004*TMath::Abs(point->GetCPoint().GetZ()));
775 //        //ampc               *=zfactor; 
776 //      }
777 //      else{ 
778 //        //ampc = 1.0*point->GetCPoint().GetMax(); 
779 //        ampc = 1.0*cl->GetMax(); 
780 //        //ampc = 1.0*point->GetCPoint().GetQ(); 
781 //        //AliTPCClusterPoint & p = point->GetCPoint();
782 //        // Float_t dy = TMath::Abs(Int_t( TMath::Abs(p.GetY()/0.4)) - TMath::Abs(p.GetY()/0.4)+0.5);
783 //        //Float_t iz =  (250.0-TMath::Abs(p.GetZ())+0.11)/0.566;
784 //        //Float_t dz = 
785 //        //  TMath::Abs( Int_t(iz) - iz + 0.5);
786
787 //        //ampc *= 1.15*(1-0.3*dy);
788 //        //ampc *= 1.15*(1-0.3*dz);
789 //        //    Float_t zfactor = (1.02-0.000*TMath::Abs(point->GetCPoint().GetZ()));
790 //        //ampc               *=zfactor; 
791
792 //      }
793 //      ampc *= 2.0;     // put mean value to channel 50
794 //      //ampc *= 0.58;     // put mean value to channel 50
795 //      Float_t w      =  1.;
796 //      //      if (type>0)  w =  1./(type/2.-0.5); 
797 //      //      Float_t z = TMath::Abs(cl->GetZ());
798 //      if (i<64) {
799 //        ampc /= 0.6;
800 //        //ampc /= (1+0.0008*z);
801 //      } else
802 //        if (i>128){
803 //          ampc /=1.5;
804 //          //ampc /= (1+0.0008*z);
805 //        }else{
806 //          //ampc /= (1+0.0008*z);
807 //        }
808         
809 //      if (type<0) {  //amp at the border - lower weight
810 //        // w*= 2.;
811           
812 //        continue;
813 //      }
814 //      if (rsigma>1.5) ampc/=1.3;  // if big backround
815 //      amp[nc[of]]        = ampc;
816 //      amp[nc[of]]       /=gainGG;
817 //      angular[nc[of]]    = TMath::Sqrt(1.+angley*angley+anglez*anglez);
818 //      weight[nc[of]]     = w;
819 //      nc[of]++;
820 //       }
821     
822 //     TMath::Sort(nc[of],amp,index,kFALSE);
823 //     Float_t sumamp=0;
824 //     Float_t sumamp2=0;
825 //     Float_t sumw=0;
826 //     //meanlog = amp[index[Int_t(nc[of]*0.33)]];
827 //     meanlog = 50;
828 //     for (Int_t i=int(nc[of]*low+0.5);i<int(nc[of]*up+0.5);i++){
829 //       Float_t ampl      = amp[index[i]]/angular[index[i]];
830 //       ampl              = meanlog*TMath::Log(1.+ampl/meanlog);
831 //       //
832 //       sumw    += weight[index[i]]; 
833 //       sumamp  += weight[index[i]]*ampl;
834 //       sumamp2 += weight[index[i]]*ampl*ampl;
835 //       norm[of]    += angular[index[i]]*weight[index[i]];
836 //     }
837 //     if (sumw<1){ 
838 //       SetdEdx(0);  
839 //     }
840 //     else {
841 //       norm[of] /= sumw;
842 //       mean[of]  = sumamp/sumw;
843 //       sigma[of] = sumamp2/sumw-mean[of]*mean[of];
844 //       if (sigma[of]>0.1) 
845 //      sigma[of] = TMath::Sqrt(sigma[of]);
846 //       else
847 //      sigma[of] = 1000;
848       
849 //     mean[of] = (TMath::Exp(mean[of]/meanlog)-1)*meanlog;
850 //     //mean  *=(1-0.02*(sigma/(mean*0.17)-1.));
851 //     //mean *=(1-0.1*(norm-1.));
852 //     }
853 //   }
854
855 //   Float_t dedx =0;
856 //   fSdEdx =0;
857 //   fMAngular =0;
858 //   //  mean[0]*= (1-0.05*(sigma[0]/(0.01+mean[1]*0.18)-1));
859 //   //  mean[1]*= (1-0.05*(sigma[1]/(0.01+mean[0]*0.18)-1));
860
861   
862 //   //  dedx = (mean[0]* TMath::Sqrt((1.+nc[0]))+ mean[1]* TMath::Sqrt((1.+nc[1])) )/ 
863 //   //  (  TMath::Sqrt((1.+nc[0]))+TMath::Sqrt((1.+nc[1])));
864
865 //   Int_t norm2 = 0;
866 //   Int_t norm3 = 0;
867 //   for (Int_t i =0;i<4;i++){
868 //     if (nc[i]>2&&nc[i]<1000){
869 //       dedx      += mean[i] *nc[i];
870 //       fSdEdx    += sigma[i]*(nc[i]-2);
871 //       fMAngular += norm[i] *nc[i];    
872 //       norm2     += nc[i];
873 //       norm3     += nc[i]-2;
874 //     }
875 //     fDEDX[i]  = mean[i];             
876 //     fSDEDX[i] = sigma[i];            
877 //     fNCDEDX[i]= nc[i]; 
878 //   }
879
880 //   if (norm3>0){
881 //     dedx   /=norm2;
882 //     fSdEdx /=norm3;
883 //     fMAngular/=norm2;
884 //   }
885 //   else{
886 //     SetdEdx(0);
887 //     return 0;
888 //   }
889 //   //  Float_t dedx1 =dedx;
890 //   /*
891 //   dedx =0;
892 //   for (Int_t i =0;i<4;i++){
893 //     if (nc[i]>2&&nc[i]<1000){
894 //       mean[i]   = mean[i]*(1-0.12*(sigma[i]/(fSdEdx)-1.));
895 //       dedx      += mean[i] *nc[i];
896 //     }
897 //     fDEDX[i]  = mean[i];                
898 //   }
899 //   dedx /= norm2;
900 //   */
901
902   
903 //   SetdEdx(dedx);
904 //   return dedx;
905 }
906
907 void AliTPCseed::CookPID()
908 {
909   //
910   // cook PID information according dEdx
911   //
912   Double_t fRange = 10.;
913   Double_t fRes   = 0.1;
914   Double_t fMIP   = 47.;
915   //
916   Int_t ns=AliPID::kSPECIES;
917   Double_t sumr =0;
918   for (Int_t j=0; j<ns; j++) {
919     Double_t mass=AliPID::ParticleMass(j);
920     Double_t mom=GetP();
921     Double_t dedx=fdEdx/fMIP;
922     Double_t bethe=AliMathBase::BetheBlochAleph(mom/mass); 
923     Double_t sigma=fRes*bethe;
924     if (sigma>0.001){
925       if (TMath::Abs(dedx-bethe) > fRange*sigma) {
926         fTPCr[j]=TMath::Exp(-0.5*fRange*fRange)/sigma;
927         sumr+=fTPCr[j];
928         continue;
929       }
930       fTPCr[j]=TMath::Exp(-0.5*(dedx-bethe)*(dedx-bethe)/(sigma*sigma))/sigma;
931       sumr+=fTPCr[j];
932     }
933     else{
934       fTPCr[j]=1.;
935       sumr+=fTPCr[j];
936     }
937   }
938   for (Int_t j=0; j<ns; j++) {
939     fTPCr[j]/=sumr;           //normalize
940   }
941 }
942
943 Double_t AliTPCseed::GetYat(Double_t xk) const {
944 //-----------------------------------------------------------------
945 // This function calculates the Y-coordinate of a track at the plane x=xk.
946 //-----------------------------------------------------------------
947   if (TMath::Abs(GetSnp())>AliTPCReconstructor::GetMaxSnpTrack()) return 0.; //patch 01 jan 06
948     Double_t c1=GetSnp(), r1=TMath::Sqrt((1.-c1)*(1.+c1));
949     Double_t c2=c1+GetC()*(xk-GetX());
950     if (TMath::Abs(c2)>AliTPCReconstructor::GetMaxSnpTrack()) return 0;
951     Double_t r2=TMath::Sqrt((1.-c2)*(1.+c2));
952     return GetY() + (xk-GetX())*(c1+c2)/(r1+r2);
953 }
954
955
956
957 Float_t  AliTPCseed::CookdEdxNorm(Double_t low, Double_t up, Int_t type, Int_t i1, Int_t i2, Bool_t shapeNorm,Int_t posNorm, Int_t padNorm, Int_t returnVal){
958  
959   //
960   // calculates dedx using the cluster
961   // low    -  up specify trunc mean range  - default form 0-0.7
962   // type   -  1 - max charge  or 0- total charge in cluster 
963   //           //2- max no corr 3- total+ correction
964   // i1-i2  -  the pad-row range used for calculation
965   // shapeNorm - kTRUE  -taken from OCDB
966   //           
967   // posNorm   - usage of pos normalization 
968   // padNorm   - pad type normalization
969   // returnVal - 0 return mean
970   //           - 1 return RMS
971   //           - 2 return number of clusters
972   //           
973   // normalization parametrization taken from AliTPCClusterParam
974   //
975   AliTPCClusterParam * parcl = AliTPCcalibDB::Instance()->GetClusterParam();
976   AliTPCParam * param = AliTPCcalibDB::Instance()->GetParameters();
977   if (!parcl)  return 0;
978   if (!param) return 0;
979   Int_t row0 = param->GetNRowLow();
980   Int_t row1 = row0+param->GetNRowUp1();
981
982   Float_t amp[160];
983   Int_t   indexes[160];
984   Int_t   ncl=0;
985   //
986   //
987   Float_t gainGG      = 1;  // gas gain factor -always enabled
988   Float_t gainPad     = 1;  // gain map  - used always
989   Float_t corrShape   = 1;  // correction due angular effect, diffusion and electron attachment
990   Float_t corrPos     = 1;  // local position correction - if posNorm enabled
991   Float_t corrPadType = 1;  // pad type correction - if padNorm enabled
992   Float_t corrNorm    = 1;  // normalization factor - set Q to channel 50
993   //   
994   //
995   //
996   if (AliTPCcalibDB::Instance()->GetParameters()){
997     gainGG= AliTPCcalibDB::Instance()->GetParameters()->GetGasGain()/20000;  //relative gas gain
998     gainGG *= AliTPCcalibDB::Instance()->GetParameters()->GetNtot()/36.82;//correction for the ionisation
999   }
1000
1001   const Float_t ktany = TMath::Tan(TMath::DegToRad()*10);
1002   const Float_t kedgey =3.;
1003   //
1004   //
1005   for (Int_t irow=i1; irow<i2; irow++){
1006     AliTPCclusterMI* cluster = GetClusterPointer(irow);
1007     if (!cluster) continue;
1008     if (TMath::Abs(cluster->GetY())>cluster->GetX()*ktany-kedgey) continue; // edge cluster
1009     Float_t charge= (type%2)? cluster->GetMax():cluster->GetQ();
1010     Int_t  ipad= 0;
1011     if (irow>=row0) ipad=1;
1012     if (irow>=row1) ipad=2;    
1013     //
1014     //
1015     //
1016     AliTPCCalPad * gainMap =  AliTPCcalibDB::Instance()->GetDedxGainFactor();
1017     if (gainMap) {
1018       //
1019       // Get gainPad - pad by pad calibration
1020       //
1021       Float_t factor = 1;      
1022       AliTPCCalROC * roc = gainMap->GetCalROC(cluster->GetDetector());
1023       if (irow < row0) { // IROC
1024         factor = roc->GetValue(irow, TMath::Nint(cluster->GetPad()));
1025       } else {         // OROC
1026         factor = roc->GetValue(irow - row0, TMath::Nint(cluster->GetPad()));
1027       }
1028       if (factor>0.5) gainPad=factor;
1029     }
1030     //
1031     //do position and angular normalization
1032     //
1033     if (shapeNorm){
1034       if (type<=1){
1035         //      
1036         AliTPCTrackerPoint * point = GetTrackPoint(irow);
1037         Float_t              ty = TMath::Abs(point->GetAngleY());
1038         Float_t              tz = TMath::Abs(point->GetAngleZ()*TMath::Sqrt(1+ty*ty));
1039         
1040         Float_t dr    = (250.-TMath::Abs(cluster->GetZ()))/250.;
1041         corrShape  = parcl->Qnorm(ipad,type,dr,ty,tz);
1042       }
1043     }
1044     
1045     if (posNorm>0){
1046       //
1047       // Do position normalization - relative distance to 
1048       // center of pad- time bin
1049       // Work in progress
1050       //      corrPos = parcl->QnormPos(ipad,type, cluster->GetPad(),
1051       //                                cluster->GetTimeBin(), cluster->GetZ(),
1052       //                                cluster->GetSigmaY2(),cluster->GetSigmaZ2(),
1053       //                                cluster->GetMax(),cluster->GetQ());
1054       // scaled response function
1055       Float_t yres0 = parcl->GetRMS0(0,ipad,0,0)/param->GetPadPitchWidth(cluster->GetDetector());
1056       Float_t zres0 = parcl->GetRMS0(1,ipad,0,0)/param->GetZWidth();
1057       //
1058       
1059       AliTPCTrackerPoint * point = GetTrackPoint(irow);
1060       Float_t              ty = TMath::Abs(point->GetAngleY());
1061       Float_t              tz = TMath::Abs(point->GetAngleZ()*TMath::Sqrt(1+ty*ty));
1062       
1063       if (type==1) corrPos = 
1064         parcl->QmaxCorrection(cluster->GetDetector(), cluster->GetRow(),cluster->GetPad(), 
1065                               cluster->GetTimeBin(),ty,tz,yres0,zres0,0.4);
1066       if (type==0) corrPos = 
1067         parcl->QtotCorrection(cluster->GetDetector(), cluster->GetRow(),cluster->GetPad(), 
1068                               cluster->GetTimeBin(),ty,tz,yres0,zres0,cluster->GetQ(),2.5,0.4);
1069       if (posNorm==3){
1070         Float_t dr    = (250.-TMath::Abs(cluster->GetZ()))/250.;
1071         Double_t signtgl = (cluster->GetZ()*point->GetAngleZ()>0)? 1:-1;
1072         Double_t p2 = TMath::Abs(TMath::Sin(TMath::ATan(ty)));
1073         Float_t corrHis = parcl->QnormHis(ipad,type,dr,p2,TMath::Abs(point->GetAngleZ())*signtgl);
1074         if (corrHis>0) corrPos*=corrHis;
1075       }
1076
1077     }
1078
1079     if (padNorm==1){
1080       //taken from OCDB
1081       if (type==0 && parcl->QpadTnorm()) corrPadType = (*parcl->QpadTnorm())[ipad];
1082       if (type==1 && parcl->QpadMnorm()) corrPadType = (*parcl->QpadMnorm())[ipad];
1083
1084     }
1085     if (padNorm==2){
1086       corrPadType  =param->GetPadPitchLength(cluster->GetDetector(),cluster->GetRow());
1087       //use hardwired - temp fix
1088       if (type==0) corrNorm=3.;
1089       if (type==1) corrNorm=1.;
1090     }
1091     //
1092     //
1093     amp[ncl]=charge;
1094     amp[ncl]/=gainGG;                 // normalized gas gain
1095     amp[ncl]/=gainPad;                // 
1096     amp[ncl]/=corrShape;
1097     amp[ncl]/=corrPadType;
1098     amp[ncl]/=corrPos;
1099     amp[ncl]/=corrNorm; 
1100     //
1101     ncl++;
1102   }
1103
1104   if (type>3) return ncl; 
1105   TMath::Sort(ncl,amp, indexes, kFALSE);
1106
1107   if (ncl<10) return 0;
1108   
1109   Float_t suma=0;
1110   Float_t suma2=0;  
1111   Float_t sumn=0;
1112   Int_t icl0=TMath::Nint(ncl*low);
1113   Int_t icl1=TMath::Nint(ncl*up);
1114   for (Int_t icl=icl0; icl<icl1;icl++){
1115     suma+=amp[indexes[icl]];
1116     suma2+=amp[indexes[icl]]*amp[indexes[icl]];
1117     sumn++;
1118   }
1119   Float_t mean =suma/sumn;
1120   Float_t rms  =TMath::Sqrt(TMath::Abs(suma2/sumn-mean*mean));
1121   //
1122   // do time-dependent correction for pressure and temperature variations
1123   UInt_t runNumber = 1;
1124   Float_t corrTimeGain = 1;
1125   AliTPCTransform * trans = AliTPCcalibDB::Instance()->GetTransform();
1126   const AliTPCRecoParam * recoParam = AliTPCcalibDB::Instance()->GetTransform()->GetCurrentRecoParam();
1127   if (trans && recoParam->GetUseGainCorrectionTime()>0) {
1128     runNumber = trans->GetCurrentRunNumber();
1129     //AliTPCcalibDB::Instance()->SetRun(runNumber);
1130     TObjArray * timeGainSplines = AliTPCcalibDB::Instance()->GetTimeGainSplinesRun(runNumber);
1131     if (timeGainSplines) {
1132       UInt_t time = trans->GetCurrentTimeStamp();
1133       AliSplineFit * fitMIP = (AliSplineFit *) timeGainSplines->At(0);
1134       AliSplineFit * fitFPcosmic = (AliSplineFit *) timeGainSplines->At(1);
1135       if (fitMIP) {
1136         corrTimeGain = AliTPCcalibDButil::EvalGraphConst(fitMIP, time);/*fitMIP->Eval(time);*/
1137       } else {
1138         if (fitFPcosmic) corrTimeGain = AliTPCcalibDButil::EvalGraphConst(fitFPcosmic, time);/*fitFPcosmic->Eval(time);*/ 
1139       }
1140     }
1141   }
1142   mean /= corrTimeGain;
1143   rms /= corrTimeGain;
1144   //
1145   if (returnVal==1) return rms;
1146   if (returnVal==2) return ncl;
1147   return mean;
1148 }
1149
1150 Float_t  AliTPCseed::CookdEdxAnalytical(Double_t low, Double_t up, Int_t type, Int_t i1, Int_t i2, Int_t returnVal, Int_t rowThres, Int_t mode, TVectorT<float> *returnVec){
1151  
1152   //
1153   // calculates dedx using the cluster
1154   // low    -  up specify trunc mean range  - default form 0-0.7
1155   // type   -  1 - max charge  or 0- total charge in cluster 
1156   //           //2- max no corr 3- total+ correction
1157   // i1-i2  -  the pad-row range used for calculation
1158   //           
1159   // posNorm   - usage of pos normalization 
1160   // returnVal - 0  return mean
1161   //           - 1  return RMS
1162   //           - 2  return number of clusters
1163   //           - 3  ratio
1164   //           - 4  mean upper half
1165   //           - 5  mean  - lower half
1166   //           - 6  third moment
1167   // mode      - 0 - linear
1168   //           - 1 - logatithmic
1169   // rowThres  - number of rows before and after given pad row to check for clusters below threshold
1170   //           
1171   // normalization parametrization taken from AliTPCClusterParam
1172   //
1173   if (returnVec) returnVec->ResizeTo(10);
1174
1175   AliTPCClusterParam * parcl = AliTPCcalibDB::Instance()->GetClusterParam();
1176   AliTPCParam * param = AliTPCcalibDB::Instance()->GetParameters();
1177   AliTPCTransform * trans = AliTPCcalibDB::Instance()->GetTransform();
1178   const AliTPCRecoParam * recoParam = AliTPCcalibDB::Instance()->GetTransform()->GetCurrentRecoParam();
1179   if (!parcl)  return 0;
1180   if (!param) return 0;
1181   Int_t row0 = param->GetNRowLow();
1182   Int_t row1 = row0+param->GetNRowUp1();
1183
1184   Float_t amp[160];
1185   Int_t   indexes[160];
1186   Int_t   ncl=0;
1187   Int_t   nclBelowThr = 0; // counts number of clusters below threshold
1188   //
1189   //
1190   Float_t gainGG      = 1;  // gas gain factor -always enabled
1191   Float_t gainPad     = 1;  // gain map  - used always
1192   Float_t corrPos     = 1;  // local position correction - if posNorm enabled
1193   //   
1194   //
1195   //
1196   if (AliTPCcalibDB::Instance()->GetParameters()){
1197     gainGG= AliTPCcalibDB::Instance()->GetParameters()->GetGasGain()/20000;  //relative gas gain
1198     gainGG *= AliTPCcalibDB::Instance()->GetParameters()->GetNtot()/36.82;//correction for the ionisation
1199   }
1200   Double_t timeCut=0;
1201   if (AliTPCcalibDB::Instance()->IsTrgL0()){
1202     // by defualt we assume L1 trigger is used - make a correction in case of  L0
1203     AliCTPTimeParams* ctp = AliTPCcalibDB::Instance()->GetCTPTimeParams();
1204     Double_t delay = ctp->GetDelayL1L0()*0.000000025;
1205     delay/=param->GetTSample();    
1206     timeCut=delay;
1207   }
1208   timeCut += recoParam->GetSkipTimeBins();
1209
1210   //
1211   // extract time-dependent correction for pressure and temperature variations
1212   //
1213   UInt_t runNumber = 1;
1214   Float_t corrTimeGain = 1;
1215   TObjArray * timeGainSplines = 0x0;
1216   TGraphErrors * grPadEqual = 0x0;
1217   TGraphErrors*  grChamberGain[4]={0x0,0x0,0x0,0x0};
1218   TF1*  funDipAngle[4]={0x0,0x0,0x0,0x0};
1219   //
1220   //
1221   if (recoParam->GetNeighborRowsDedx() == 0) rowThres = 0;      
1222   UInt_t time = 1;//
1223   //
1224   if (trans) {
1225       runNumber = trans->GetCurrentRunNumber();
1226       time = trans->GetCurrentTimeStamp();
1227       //AliTPCcalibDB::Instance()->SetRun(runNumber);
1228       timeGainSplines = AliTPCcalibDB::Instance()->GetTimeGainSplinesRun(runNumber);
1229       if (timeGainSplines && recoParam->GetUseGainCorrectionTime()>0) {
1230         AliSplineFit * fitMIP = (AliSplineFit *) timeGainSplines->At(0);
1231         AliSplineFit * fitFPcosmic = (AliSplineFit *) timeGainSplines->At(1);
1232         if (fitMIP) {
1233           corrTimeGain =  AliTPCcalibDButil::EvalGraphConst(fitMIP, time); /*fitMIP->Eval(time);*/
1234         } else {
1235           if (fitFPcosmic) corrTimeGain = AliTPCcalibDButil::EvalGraphConst(fitFPcosmic, time); /*fitFPcosmic->Eval(time); */
1236         }
1237         //
1238         if (type==1) grPadEqual = (TGraphErrors * ) timeGainSplines->FindObject("TGRAPHERRORS_MEANQMAX_PADREGIONGAIN_BEAM_ALL");
1239         if (type==0) grPadEqual = (TGraphErrors * ) timeGainSplines->FindObject("TGRAPHERRORS_MEANQTOT_PADREGIONGAIN_BEAM_ALL");
1240         const char* names[4]={"SHORT","MEDIUM","LONG","ABSOLUTE"};
1241         for (Int_t iPadRegion=0; iPadRegion<4; ++iPadRegion) {
1242           grChamberGain[iPadRegion]=(TGraphErrors*)timeGainSplines->FindObject(Form("TGRAPHERRORS_MEAN_CHAMBERGAIN_%s_BEAM_ALL",names[iPadRegion]));
1243           if (type==1) funDipAngle[iPadRegion]=(TF1*)timeGainSplines->FindObject(Form("TF1_QMAX_DIPANGLE_%s_BEAM_ALL",names[iPadRegion]));
1244           if (type==0) funDipAngle[iPadRegion]=(TF1*)timeGainSplines->FindObject(Form("TF1_QTOT_DIPANGLE_%s_BEAM_ALL",names[iPadRegion]));
1245         }
1246       }
1247   }
1248   
1249   const Float_t kClusterShapeCut = 1.5; // IMPPRTANT TO DO: move value to AliTPCRecoParam
1250   const Float_t ktany = TMath::Tan(TMath::DegToRad()*10);
1251   const Float_t kedgey =3.;
1252   //
1253   //
1254   for (Int_t irow=i1; irow<i2; irow++){
1255     AliTPCclusterMI* cluster = GetClusterPointer(irow);
1256     if (!cluster && irow > 1 && irow < 157) {
1257       Bool_t isClBefore = kFALSE;
1258       Bool_t isClAfter  = kFALSE;
1259       for(Int_t ithres = 1; ithres <= rowThres; ithres++) {
1260         AliTPCclusterMI * clusterBefore = GetClusterPointer(irow - ithres);
1261         if (clusterBefore) isClBefore = kTRUE;
1262         AliTPCclusterMI * clusterAfter  = GetClusterPointer(irow + ithres);
1263         if (clusterAfter) isClAfter = kTRUE;
1264       }
1265       if (isClBefore && isClAfter) nclBelowThr++;
1266     }
1267     if (!cluster) continue;
1268     if (cluster->GetTimeBin()<timeCut) continue; //reject  clusters at the gating grid opening
1269     //
1270     //
1271     if (TMath::Abs(cluster->GetY())>cluster->GetX()*ktany-kedgey) continue; // edge cluster
1272     //
1273     AliTPCTrackerPoint * point = GetTrackPoint(irow);
1274     if (point==0) continue;    
1275     Float_t rsigmay = TMath::Sqrt(point->GetSigmaY());
1276     if (rsigmay > kClusterShapeCut) continue;
1277     //
1278     if (cluster->IsUsed(11)) continue; // remove shared clusters for PbPb
1279     //
1280     Float_t charge= (type%2)? cluster->GetMax():cluster->GetQ();
1281     Int_t  ipad= 0;
1282     if (irow>=row0) ipad=1;
1283     if (irow>=row1) ipad=2;    
1284     //
1285     //
1286     //
1287     AliTPCCalPad * gainMap =  AliTPCcalibDB::Instance()->GetDedxGainFactor();
1288     if (gainMap) {
1289       //
1290       // Get gainPad - pad by pad calibration
1291       //
1292       Float_t factor = 1;      
1293       AliTPCCalROC * roc = gainMap->GetCalROC(cluster->GetDetector());
1294       if (irow < row0) { // IROC
1295         factor = roc->GetValue(irow, TMath::Nint(cluster->GetPad()));
1296       } else {         // OROC
1297         factor = roc->GetValue(irow - row0, TMath::Nint(cluster->GetPad()));
1298       }
1299       if (factor>0.3) gainPad=factor;
1300     }
1301     //
1302     // Do position normalization - relative distance to 
1303     // center of pad- time bin
1304     
1305     Float_t              ty = TMath::Abs(point->GetAngleY());
1306     Float_t              tz = TMath::Abs(point->GetAngleZ()*TMath::Sqrt(1+ty*ty));
1307     Float_t yres0 = parcl->GetRMS0(0,ipad,0,0)/param->GetPadPitchWidth(cluster->GetDetector());
1308     Float_t zres0 = parcl->GetRMS0(1,ipad,0,0)/param->GetZWidth();
1309
1310     yres0 *=parcl->GetQnormCorr(ipad, type,0);
1311     zres0 *=parcl->GetQnormCorr(ipad, type,1);
1312     Float_t effLength=parcl->GetQnormCorr(ipad, type,4)*0.5;
1313     Float_t effDiff  =(parcl->GetQnormCorr(ipad, type,2)+parcl->GetQnormCorr(ipad, type,3))*0.5;
1314     //
1315     if (type==1) {
1316       corrPos = parcl->GetQnormCorr(ipad, type,5)*
1317         parcl->QmaxCorrection(cluster->GetDetector(), cluster->GetRow(),cluster->GetPad(), 
1318                               cluster->GetTimeBin(),ty,tz,yres0,zres0,effLength,effDiff);
1319       Float_t drm   = 0.5-TMath::Abs(cluster->GetZ()/250.);
1320       corrPos*=(1+parcl->GetQnormCorr(ipad, type+2,0)*drm);
1321       corrPos*=(1+parcl->GetQnormCorr(ipad, type+2,1)*ty*ty);
1322       corrPos*=(1+parcl->GetQnormCorr(ipad, type+2,2)*tz*tz);
1323       //
1324     }
1325     if (type==0) {
1326       corrPos = parcl->GetQnormCorr(ipad, type,5)*
1327         parcl->QtotCorrection(cluster->GetDetector(), cluster->GetRow(),cluster->GetPad(), 
1328                               cluster->GetTimeBin(),ty,tz,yres0,zres0,cluster->GetQ(),2.5,effLength,effDiff);
1329       
1330       Float_t drm   = 0.5-TMath::Abs(cluster->GetZ()/250.);
1331       corrPos*=(1+parcl->GetQnormCorr(ipad, type+2,0)*drm);
1332       corrPos*=(1+parcl->GetQnormCorr(ipad, type+2,1)*ty*ty);
1333       corrPos*=(1+parcl->GetQnormCorr(ipad, type+2,2)*tz*tz);
1334       //
1335     }
1336     //
1337     // pad region equalization outside of cluster param
1338     //
1339     Float_t gainEqualPadRegion = 1;
1340     if (grPadEqual && recoParam->GetUseGainCorrectionTime()>0) gainEqualPadRegion = grPadEqual->Eval(ipad);
1341     //
1342     // chamber-by-chamber equalization outside gain map
1343     //
1344     Float_t gainChamber = 1;
1345     if (grChamberGain[ipad] && recoParam->GetUseGainCorrectionTime()>0) {
1346       gainChamber = grChamberGain[ipad]->Eval(cluster->GetDetector());
1347       if (gainChamber==0) gainChamber=1; // in case old calibation was used before use no correction
1348     }
1349     //
1350     // dip angle correction
1351     //
1352     Float_t corrDipAngle = 1;
1353     Float_t corrDipAngleAbs = 1;
1354     //    if (grDipAngle[ipad]) corrDipAngle = grDipAngle[ipad]->Eval(GetTgl());
1355     Double_t tgl=GetTgl();
1356     if (funDipAngle[ipad]) corrDipAngle = funDipAngle[ipad]->Eval(tgl);
1357     if (funDipAngle[3]) corrDipAngleAbs = funDipAngle[3]->Eval(tgl);
1358     //
1359     // pressure temperature and high voltage correction
1360     //
1361     Double_t correctionHVandPT = AliTPCcalibDB::Instance()->GetGainCorrectionHVandPT(time, runNumber,cluster->GetDetector(), 5 , recoParam->GetGainCorrectionHVandPTMode());
1362     //
1363     amp[ncl]=charge;
1364     amp[ncl]/=gainGG;               // nominal gas gain
1365     amp[ncl]/=correctionHVandPT;    // correction for the HV and P/T - time dependent
1366     amp[ncl]/=gainPad;              // 
1367     amp[ncl]/=corrPos;
1368     amp[ncl]/=gainEqualPadRegion;
1369     amp[ncl]/=gainChamber;
1370     amp[ncl]/=corrDipAngle;
1371     amp[ncl]/=corrDipAngleAbs;
1372     //
1373     ncl++;
1374   }
1375
1376   if (type==2) return ncl; 
1377   TMath::Sort(ncl,amp, indexes, kFALSE);
1378   //
1379   if (ncl<10) return 0;
1380   //
1381   Double_t * ampWithBelow = new Double_t[ncl + nclBelowThr];
1382   for(Int_t iCl = 0; iCl < ncl + nclBelowThr; iCl++) {
1383     if (iCl < nclBelowThr) {
1384       ampWithBelow[iCl] = amp[indexes[0]];
1385     } else {
1386       ampWithBelow[iCl] = amp[indexes[iCl - nclBelowThr]];
1387     }
1388   }
1389   //printf("DEBUG: %i shit %f", nclBelowThr, amp[indexes[0]]);
1390   //
1391   Float_t suma=0;
1392   Float_t suma2=0;  
1393   Float_t suma3=0;  
1394   Float_t sumaS=0;  
1395   Float_t sumn=0;
1396   // upper,and lower part statistic
1397   Float_t sumL=0, sumL2=0, sumLN=0;
1398   Float_t sumD=0, sumD2=0, sumDN=0;
1399
1400   Int_t icl0=TMath::Nint((ncl + nclBelowThr)*low);
1401   Int_t icl1=TMath::Nint((ncl + nclBelowThr)*up);
1402   Int_t iclm=TMath::Nint((ncl + nclBelowThr)*(low +(up+low)*0.5));
1403   //
1404   for (Int_t icl=icl0; icl<icl1;icl++){
1405     if (ampWithBelow[icl]<0.1) continue;
1406     Double_t camp=ampWithBelow[icl]/corrTimeGain;
1407     if (mode==1) camp= TMath::Log(camp);
1408     if (icl<icl1){
1409       suma+=camp;
1410       suma2+=camp*camp;
1411       suma3+=camp*camp*camp;
1412       sumaS+=TMath::Power(TMath::Abs(camp),1./3.);
1413       sumn++;
1414     }
1415     if (icl>iclm){
1416       sumL+=camp;
1417       sumL2+=camp*camp;
1418       sumLN++;
1419       }
1420     if (icl<=iclm){
1421       sumD+=camp;
1422       sumD2+=camp*camp;
1423       sumDN++;
1424     }
1425   }
1426   //
1427   Float_t mean = 0;
1428   Float_t meanL = 0;  
1429   Float_t meanD = 0;           // lower half mean
1430   if (sumn > 1e-30)   mean =suma/sumn;
1431   if (sumLN > 1e-30)  meanL =sumL/sumLN;
1432   if (sumDN > 1e-30)  meanD =(sumD/sumDN);
1433   /*
1434   Float_t mean =suma/sumn;
1435   Float_t meanL = sumL/sumLN;  
1436   Float_t meanD =(sumD/sumDN);           // lower half mean
1437   */
1438
1439   Float_t rms = 0;
1440   Float_t mean2=0;
1441   Float_t mean3=0;
1442   Float_t meanS=0;
1443
1444   if(sumn>0){
1445     rms = TMath::Sqrt(TMath::Abs(suma2/sumn-mean*mean));
1446     mean2=suma2/sumn;
1447     mean3=suma3/sumn;
1448     meanS=sumaS/sumn;
1449   }
1450
1451   if (mean2>0) mean2=TMath::Power(TMath::Abs(mean2),1./2.);
1452   if (mean3>0) mean3=TMath::Power(TMath::Abs(mean3),1./3.);
1453   if (meanS>0) meanS=TMath::Power(TMath::Abs(meanS),3.);
1454   //
1455   if (mode==1) mean=TMath::Exp(mean);
1456   if (mode==1) meanL=TMath::Exp(meanL);  // upper truncation
1457   if (mode==1) meanD=TMath::Exp(meanD);  // lower truncation
1458   //
1459   delete [] ampWithBelow; //return?
1460   
1461
1462   //
1463   if(returnVec){
1464       (*returnVec)(0) = mean;
1465       (*returnVec)(1) = rms;
1466       (*returnVec)(2) = ncl;
1467       (*returnVec)(3) = Double_t(nclBelowThr)/Double_t(nclBelowThr+ncl);
1468       (*returnVec)(4) = meanL;
1469       (*returnVec)(5) = meanD;
1470       (*returnVec)(6) = mean2;
1471       (*returnVec)(7) = mean3;
1472       (*returnVec)(8) = meanS;
1473       (*returnVec)(9) = nclBelowThr;
1474   }
1475
1476   if (returnVal==1) return rms;
1477   if (returnVal==2) return ncl;
1478   if (returnVal==3) return Double_t(nclBelowThr)/Double_t(nclBelowThr+ncl);
1479   if (returnVal==4) return meanL;
1480   if (returnVal==5) return meanD;
1481   if (returnVal==6) return mean2;
1482   if (returnVal==7) return mean3;
1483   if (returnVal==8) return meanS;
1484   if (returnVal==9) return nclBelowThr;
1485   return mean;
1486 }
1487
1488
1489
1490
1491 Float_t  AliTPCseed::CookShape(Int_t type){
1492   //
1493   //
1494   //
1495  //-----------------------------------------------------------------
1496   // This funtion calculates dE/dX within the "low" and "up" cuts.
1497   //-----------------------------------------------------------------
1498   Float_t means=0;
1499   Float_t meanc=0;
1500   for (Int_t i =0; i<160;i++)    {
1501     AliTPCTrackerPoint * point = GetTrackPoint(i);
1502     if (point==0) continue;
1503
1504     AliTPCclusterMI * cl = fClusterPointer[i];
1505     if (cl==0) continue;        
1506     
1507     Float_t rsigmay =  TMath::Sqrt(point->GetSigmaY());
1508     Float_t rsigmaz =  TMath::Sqrt(point->GetSigmaZ());
1509     Float_t rsigma =   (rsigmay+rsigmaz)*0.5;
1510     if (type==0) means+=rsigma;
1511     if (type==1) means+=rsigmay;
1512     if (type==2) means+=rsigmaz;
1513     meanc++;
1514   }
1515   Float_t mean = (meanc>0)? means/meanc:0;
1516   return mean;
1517 }
1518
1519
1520
1521 Int_t  AliTPCseed::RefitTrack(AliTPCseed *seed, AliExternalTrackParam * parin, AliExternalTrackParam * parout){
1522   //
1523   // Refit the track
1524   // return value - number of used clusters
1525   // 
1526   //
1527   const Int_t kMinNcl =10;
1528   AliTPCseed *track=new AliTPCseed(*seed);
1529   Int_t sector=-1;
1530   // reset covariance
1531   //
1532   Double_t covar[15];
1533   for (Int_t i=0;i<15;i++) covar[i]=0;
1534   covar[0]=10.*10.;
1535   covar[2]=10.*10.;
1536   covar[5]=10.*10./(64.*64.);
1537   covar[9]=10.*10./(64.*64.);
1538   covar[14]=1*1;
1539   //
1540
1541   Float_t xmin=1000, xmax=-10000;
1542   Int_t imin=158, imax=0;
1543   for (Int_t i=0;i<160;i++) {
1544     AliTPCclusterMI *c=track->GetClusterPointer(i);
1545     if (!c || (track->GetClusterIndex(i) & 0x8000)) continue; 
1546     if (sector<0) sector = c->GetDetector();
1547     if (c->GetX()<xmin) xmin=c->GetX();
1548     if (c->GetX()>xmax) xmax=c->GetX();
1549     if (i<imin) imin=i;
1550     if (i>imax) imax=i;
1551   }
1552   if(imax-imin<kMinNcl) {
1553     delete track;
1554     return 0 ;
1555   }
1556   // Not succes to rotate
1557   if (!track->Rotate(TMath::DegToRad()*(sector%18*20.+10.)-track->GetAlpha())) {
1558     delete track;
1559     return 0;
1560   }
1561   //
1562   //
1563   // fit from inner to outer row
1564   //
1565   AliExternalTrackParam paramIn;
1566   AliExternalTrackParam paramOut;
1567   Bool_t isOK=kTRUE;
1568   Int_t ncl=0;
1569   //
1570   //
1571   //
1572   for (Int_t i=imin; i<=imax; i++){
1573     AliTPCclusterMI *c=track->GetClusterPointer(i);
1574     if (!c || (track->GetClusterIndex(i) & 0x8000)) continue; 
1575     //    if (RejectCluster(c,track)) continue;
1576     sector = (c->GetDetector()%18);
1577     if (!track->Rotate(TMath::DegToRad()*(sector%18*20.+10.)-track->GetAlpha())) {
1578       //continue;
1579     }
1580     Double_t r[3]={c->GetX(),c->GetY(),c->GetZ()};
1581     Double_t cov[3]={0.01,0.,0.01}; //TODO: correct error parametrisation
1582     if (!track->PropagateTo(r[0])) {
1583       isOK=kFALSE;
1584     }
1585     if ( !((static_cast<AliExternalTrackParam*>(track)->Update(&r[1],cov)))) isOK=kFALSE;
1586   }
1587   if (!isOK) { delete track; return 0;}
1588   track->AddCovariance(covar);
1589   //
1590   //
1591   //
1592   for (Int_t i=imax; i>=imin; i--){
1593     AliTPCclusterMI *c=track->GetClusterPointer(i);
1594     if (!c || (track->GetClusterIndex(i) & 0x8000)) continue;
1595     //if (RejectCluster(c,track)) continue;
1596     sector = (c->GetDetector()%18);
1597     if (!track->Rotate(TMath::DegToRad()*(sector%18*20.+10.)-track->GetAlpha())) {
1598       //continue;
1599     }
1600     Double_t r[3]={c->GetX(),c->GetY(),c->GetZ()};
1601     Double_t cov[3]={0.01,0.,0.01}; //TODO: correct error parametrisation
1602     if (!track->PropagateTo(r[0])) {
1603       isOK=kFALSE;
1604     }
1605     if ( !((static_cast<AliExternalTrackParam*>(track)->Update(&r[1],cov)))) isOK=kFALSE;
1606   }
1607   //if (!isOK) { delete track; return 0;}
1608   paramIn = *track;
1609   track->AddCovariance(covar);
1610   //
1611   //
1612   for (Int_t i=imin; i<=imax; i++){
1613     AliTPCclusterMI *c=track->GetClusterPointer(i);
1614     if (!c || (track->GetClusterIndex(i) & 0x8000)) continue; 
1615     sector = (c->GetDetector()%18);
1616     if (!track->Rotate(TMath::DegToRad()*(sector%18*20.+10.)-track->GetAlpha())) {
1617       //continue;
1618     }
1619     ncl++;
1620     //if (RejectCluster(c,track)) continue;
1621     Double_t r[3]={c->GetX(),c->GetY(),c->GetZ()};
1622     Double_t cov[3]={0.01,0.,0.01}; //TODO: correct error parametrisation
1623     if (!track->PropagateTo(r[0])) {
1624       isOK=kFALSE;
1625     }
1626     if ( !((static_cast<AliExternalTrackParam*>(track)->Update(&r[1],cov)))) isOK=kFALSE;
1627   }
1628   //if (!isOK) { delete track; return 0;}
1629   paramOut=*track;
1630   //
1631   //
1632   //
1633   if (parin) (*parin)=paramIn;
1634   if (parout) (*parout)=paramOut;
1635   delete track;
1636   return ncl;
1637 }
1638
1639
1640
1641 Bool_t AliTPCseed::RefitTrack(AliTPCseed* /*seed*/, Bool_t /*out*/){
1642   //
1643   //
1644   //
1645   return kFALSE;
1646 }
1647
1648
1649
1650
1651
1652
1653 void  AliTPCseed::GetError(AliTPCclusterMI* cluster, AliExternalTrackParam * param, 
1654                                   Double_t& erry, Double_t &errz)
1655 {
1656   //
1657   // Get cluster error at given position
1658   //
1659   AliTPCClusterParam *clusterParam = AliTPCcalibDB::Instance()->GetClusterParam();
1660   Double_t tany,tanz;  
1661   Double_t snp1=param->GetSnp();
1662   tany=snp1/TMath::Sqrt((1.-snp1)*(1.+snp1));
1663   //
1664   Double_t tgl1=param->GetTgl();
1665   tanz=tgl1/TMath::Sqrt((1.-snp1)*(1.+snp1));
1666   //
1667   Int_t padSize = 0;                          // short pads
1668   if (cluster->GetDetector() >= 36) {
1669     padSize = 1;                              // medium pads 
1670     if (cluster->GetRow() > 63) padSize = 2; // long pads
1671   }
1672
1673   erry  = clusterParam->GetError0Par( 0, padSize, (250.0 - TMath::Abs(cluster->GetZ())), TMath::Abs(tany) );
1674   errz  = clusterParam->GetError0Par( 1, padSize, (250.0 - TMath::Abs(cluster->GetZ())), TMath::Abs(tanz) );
1675 }
1676
1677
1678 void  AliTPCseed::GetShape(AliTPCclusterMI* cluster, AliExternalTrackParam * param, 
1679                                   Double_t& rmsy, Double_t &rmsz)
1680 {
1681   //
1682   // Get cluster error at given position
1683   //
1684   AliTPCClusterParam *clusterParam = AliTPCcalibDB::Instance()->GetClusterParam();
1685   Double_t tany,tanz;  
1686   Double_t snp1=param->GetSnp();
1687   tany=snp1/TMath::Sqrt((1.-snp1)*(1.+snp1));
1688   //
1689   Double_t tgl1=param->GetTgl();
1690   tanz=tgl1/TMath::Sqrt((1.-snp1)*(1.+snp1));
1691   //
1692   Int_t padSize = 0;                          // short pads
1693   if (cluster->GetDetector() >= 36) {
1694     padSize = 1;                              // medium pads 
1695     if (cluster->GetRow() > 63) padSize = 2; // long pads
1696   }
1697
1698   rmsy  = clusterParam->GetRMSQ( 0, padSize, (250.0 - TMath::Abs(cluster->GetZ())), TMath::Abs(tany), TMath::Abs(cluster->GetMax()) );
1699   rmsz  = clusterParam->GetRMSQ( 1, padSize, (250.0 - TMath::Abs(cluster->GetZ())), TMath::Abs(tanz) ,TMath::Abs(cluster->GetMax()));
1700 }
1701
1702
1703
1704 Double_t AliTPCseed::GetQCorrGeom(Float_t ty, Float_t tz){
1705   //Geoetrical
1706   //ty    - tangent in local y direction
1707   //tz    - tangent 
1708   //
1709   Float_t norm=TMath::Sqrt(1+ty*ty+tz*tz);
1710   return norm;
1711 }
1712
1713 Double_t AliTPCseed::GetQCorrShape(Int_t ipad, Int_t type,Float_t z, Float_t ty, Float_t tz, Float_t /*q*/, Float_t /*thr*/){
1714   //
1715   // Q normalization
1716   //
1717   // return value =  Q Normalization factor
1718   // Normalization - 1 - shape factor part for full drift          
1719   //                 1 - electron attachment for 0 drift
1720
1721   // Input parameters:
1722   //
1723   // ipad - 0 short pad
1724   //        1 medium pad
1725   //        2 long pad
1726   //
1727   // type - 0 qmax
1728   //      - 1 qtot
1729   //
1730   //z     - z position (-250,250 cm)
1731   //ty    - tangent in local y direction
1732   //tz    - tangent 
1733   //
1734
1735   AliTPCClusterParam * paramCl = AliTPCcalibDB::Instance()->GetClusterParam();
1736   AliTPCParam   * paramTPC = AliTPCcalibDB::Instance()->GetParameters();
1737  
1738   if (!paramCl) return 1;
1739   //
1740   Double_t dr =  250.-TMath::Abs(z); 
1741   Double_t sy =  paramCl->GetRMS0( 0,ipad, dr, TMath::Abs(ty));
1742   Double_t sy0=  paramCl->GetRMS0(0,ipad, 250, 0);
1743   Double_t sz =  paramCl->GetRMS0( 1,ipad, dr, TMath::Abs(tz));
1744   Double_t sz0=  paramCl->GetRMS0(1,ipad, 250, 0);
1745
1746   Double_t sfactorMax = TMath::Sqrt(sy0*sz0/(sy*sz));
1747
1748  
1749   Double_t dt = 1000000*(dr/paramTPC->GetDriftV());  //time in microsecond
1750   Double_t attProb = TMath::Exp(-paramTPC->GetAttCoef()*paramTPC->GetOxyCont()*dt);
1751   //
1752   //
1753   if (type==0) return sfactorMax*attProb;
1754
1755   return attProb;
1756
1757
1758 }
1759
1760
1761 /*
1762 //_______________________________________________________________________
1763 Float_t AliTPCseed::GetTPCClustInfo(Int_t nNeighbours, Int_t type, Int_t row0, Int_t row1, TVectorT<float> *returnVec)
1764 {
1765   //
1766   // TPC cluster information
1767   // type 0: get fraction of found/findable clusters with neighbourhood definition
1768   //      1: found clusters
1769   //      2: findable (number of clusters above and below threshold)
1770   //
1771   // definition of findable clusters:
1772   //            a cluster is defined as findable if there is another cluster
1773   //           within +- nNeighbours pad rows. The idea is to overcome threshold
1774   //           effects with a very simple algorithm.
1775   //
1776
1777   const Float_t kClusterShapeCut = 1.5; // IMPPRTANT TO DO: move value to AliTPCRecoParam
1778   const Float_t ktany = TMath::Tan(TMath::DegToRad()*10);
1779   const Float_t kedgey =3.;
1780   
1781   Float_t ncl = 0;
1782   Float_t nclBelowThr = 0; // counts number of clusters below threshold
1783
1784   for (Int_t irow=row0; irow<row1; irow++){
1785     AliTPCclusterMI* cluster = GetClusterPointer(irow);
1786
1787     if (!cluster && irow > 1 && irow < 157) {
1788       Bool_t isClBefore = kFALSE;
1789       Bool_t isClAfter  = kFALSE;
1790       for(Int_t ithres = 1; ithres <= nNeighbours; ithres++) {
1791         AliTPCclusterMI * clusterBefore = GetClusterPointer(irow - ithres);
1792         if (clusterBefore) isClBefore = kTRUE;
1793         AliTPCclusterMI * clusterAfter  = GetClusterPointer(irow + ithres);
1794         if (clusterAfter) isClAfter = kTRUE;
1795       }
1796       if (isClBefore && isClAfter) nclBelowThr++;
1797     }
1798     if (!cluster) continue;
1799     //
1800     //
1801     if (TMath::Abs(cluster->GetY())>cluster->GetX()*ktany-kedgey) continue; // edge cluster
1802     //
1803     AliTPCTrackerPoint * point = GetTrackPoint(irow);
1804     if (point==0) continue;    
1805     Float_t rsigmay = TMath::Sqrt(point->GetSigmaY());
1806     if (rsigmay > kClusterShapeCut) continue;
1807     //
1808     if (cluster->IsUsed(11)) continue; // remove shared clusters for PbPb
1809     ncl++;
1810   }
1811   if(returnVec->GetNoElements != 3){
1812       returnVec->ResizeTo(3);
1813   }
1814   Float_t nclAll = nclBelowThr+ncl;
1815   returnVec(0) = nclAll>0?ncl/nclAll:0;
1816   returnVec(1) = ncl;
1817   returnVec(2) = nclAll;
1818
1819   if(ncl<10)
1820     return 0;
1821   if(type==0) 
1822     if(nclAll>0)
1823       return ncl/nclAll;
1824   if(type==1)
1825     return ncl;
1826   if(type==2)
1827     return nclAll;
1828   return 0;
1829 }
1830 */
1831 //_______________________________________________________________________
1832 Int_t AliTPCseed::GetNumberOfClustersIndices() {
1833   Int_t ncls = 0;
1834   for (int i=0; i < 160; i++) {
1835     if ((fIndex[i] & 0x8000) == 0)
1836       ncls++;
1837   }
1838   return ncls;
1839 }
1840
1841 //_______________________________________________________________________
1842 void AliTPCseed::Clear(Option_t*)
1843 {
1844   // formally seed may allocate memory for clusters (althought this should not happen for 
1845   // the seeds in the pool). Hence we need this method for fwd. compatibility
1846   if (fClusterOwner) for (int i=160;i--;) {delete fClusterPointer[i]; fClusterPointer[i] = 0;}
1847 }
1848
1849 TObject* AliTPCseed::Clone(const char* /*newname*/) const
1850 {
1851   // temporary override TObject::Clone to avoid crashes in reco
1852   AliTPCseed* src = (AliTPCseed*)this;
1853   AliTPCseed* dst = new AliTPCseed(*src,fClusterOwner);
1854   return dst;
1855 }