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