Changes for HLT - speedup interface for HLT tracker
[u/mrichter/AliRoot.git] / TPC / AliTPCclustererMI.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 /* $Id$ */
17
18 //-------------------------------------------------------
19 //          Implementation of the TPC clusterer
20 //
21 //  1. The Input data for reconstruction - Options
22 //      1.a Simulated data  - TTree - invoked Digits2Clusters()
23 //      1.b Raw data        - Digits2Clusters(AliRawReader* rawReader); 
24 //
25 //  2. The Output data
26 //      2.a TTree with clusters - if  SetOutput(TTree * tree) invoked
27 //      2.b TObjArray           - Faster option for HLT
28 //      2.c TClonesArray        - Faster option for HLT (smaller memory consumption), activate with fBClonesArray flag
29 //
30 //  3. Reconstruction setup
31 //     see AliTPCRecoParam for list of parameters 
32 //     The reconstruction parameterization taken from the 
33 //     AliTPCReconstructor::GetRecoParam()
34 //     Possible to setup it in reconstruction macro  AliTPCReconstructor::SetRecoParam(...)
35 //     
36 //
37 //
38 //   Origin: Marian Ivanov 
39 //-------------------------------------------------------
40
41 #include "Riostream.h"
42 #include <TF1.h>
43 #include <TFile.h>
44 #include <TGraph.h>
45 #include <TH1F.h>
46 #include <TObjArray.h>
47 #include <TClonesArray.h>
48 #include <TRandom.h>
49 #include <TTree.h>
50 #include <TTreeStream.h>
51
52 #include "AliDigits.h"
53 #include "AliLoader.h"
54 #include "AliLog.h"
55 #include "AliMathBase.h"
56 #include "AliRawEventHeaderBase.h"
57 #include "AliRawReader.h"
58 #include "AliRunLoader.h"
59 #include "AliSimDigits.h"
60 #include "AliTPCCalPad.h"
61 #include "AliTPCCalROC.h"
62 #include "AliTPCClustersArray.h"
63 #include "AliTPCClustersRow.h"
64 #include "AliTPCParam.h"
65 #include "AliTPCRawStream.h"
66 #include "AliTPCRecoParam.h"
67 #include "AliTPCReconstructor.h"
68 #include "AliTPCcalibDB.h"
69 #include "AliTPCclusterInfo.h"
70 #include "AliTPCclusterMI.h"
71 #include "AliTPCTransform.h"
72 #include "AliTPCclustererMI.h"
73
74 ClassImp(AliTPCclustererMI)
75
76
77
78 AliTPCclustererMI::AliTPCclustererMI(const AliTPCParam* par, const AliTPCRecoParam * recoParam):
79   fBins(0),
80   fSigBins(0),
81   fNSigBins(0),
82   fLoop(0),
83   fMaxBin(0),
84   fMaxTime(0),
85   fMaxPad(0),
86   fSector(-1),
87   fRow(-1),
88   fSign(0),
89   fRx(0),
90   fPadWidth(0),
91   fPadLength(0),
92   fZWidth(0),
93   fPedSubtraction(kFALSE),
94   fEventHeader(0),
95   fTimeStamp(0),
96   fEventType(0),
97   fInput(0),
98   fOutput(0),
99   fOutputArray(0),
100   fOutputClonesArray(0),
101   fRowCl(0),
102   fRowDig(0),
103   fParam(0),
104   fNcluster(0),
105   fNclusters(0),
106   fDebugStreamer(0),
107   fRecoParam(0),
108   fBDumpSignal(kFALSE),
109   fBClonesArray(kFALSE)
110 {
111   //
112   // COSNTRUCTOR
113   // param     - tpc parameters for given file
114   // recoparam - reconstruction parameters 
115   //
116   fInput =0;
117   fParam = par;
118   if (recoParam) {
119     fRecoParam = recoParam;
120   }else{
121     //set default parameters if not specified
122     fRecoParam = AliTPCReconstructor::GetRecoParam();
123     if (!fRecoParam)  fRecoParam = AliTPCRecoParam::GetLowFluxParam();
124   }
125   fDebugStreamer = new TTreeSRedirector("TPCsignal.root");
126   //  Int_t nPoints = fRecoParam->GetLastBin()-fRecoParam->GetFirstBin();
127   fRowCl= new AliTPCClustersRow();
128   fRowCl->SetClass("AliTPCclusterMI");
129   fRowCl->SetArray(1);
130
131 }
132 //______________________________________________________________
133 AliTPCclustererMI::AliTPCclustererMI(const AliTPCclustererMI &param)
134               :TObject(param),
135   fBins(0),
136   fSigBins(0),
137   fNSigBins(0),
138   fLoop(0),
139   fMaxBin(0),
140   fMaxTime(0),
141   fMaxPad(0),
142   fSector(-1),
143   fRow(-1),
144   fSign(0),
145   fRx(0),
146   fPadWidth(0),
147   fPadLength(0),
148   fZWidth(0),
149   fPedSubtraction(kFALSE),
150   fEventHeader(0),
151   fTimeStamp(0),
152   fEventType(0),
153   fInput(0),
154   fOutput(0),
155   fOutputArray(0),
156   fOutputClonesArray(0),
157   fRowCl(0),
158   fRowDig(0),
159   fParam(0),
160   fNcluster(0),
161   fNclusters(0),
162   fDebugStreamer(0),
163   fRecoParam(0),
164   fBDumpSignal(kFALSE),
165   fBClonesArray(kFALSE)
166 {
167   //
168   // dummy
169   //
170   fMaxBin = param.fMaxBin;
171 }
172 //______________________________________________________________
173 AliTPCclustererMI & AliTPCclustererMI::operator =(const AliTPCclustererMI & param)
174 {
175   //
176   // assignment operator - dummy
177   //
178   fMaxBin=param.fMaxBin;
179   return (*this);
180 }
181 //______________________________________________________________
182 AliTPCclustererMI::~AliTPCclustererMI(){
183   //
184   //
185   //
186   if (fDebugStreamer) delete fDebugStreamer;
187   if (fOutputArray){
188     //fOutputArray->Delete();
189     delete fOutputArray;
190   }
191   if (fOutputClonesArray){
192     fOutputClonesArray->Delete();
193     delete fOutputClonesArray;
194   }
195 }
196
197 void AliTPCclustererMI::SetInput(TTree * tree)
198 {
199   //
200   // set input tree with digits
201   //
202   fInput = tree;  
203   if  (!fInput->GetBranch("Segment")){
204     cerr<<"AliTPC::Digits2Clusters(): no porper input tree !\n";
205     fInput=0;
206     return;
207   }
208 }
209
210 void AliTPCclustererMI::SetOutput(TTree * tree) 
211 {
212   //
213   // Set the output tree
214   // If not set the ObjArray used - Option for HLT 
215   //
216   if (!tree) return;
217   fOutput= tree;
218   AliTPCClustersRow clrow;
219   AliTPCClustersRow *pclrow=&clrow;  
220   clrow.SetClass("AliTPCclusterMI");
221   clrow.SetArray(1); // to make Clones array
222   fOutput->Branch("Segment","AliTPCClustersRow",&pclrow,32000,200);    
223 }
224
225
226 void AliTPCclustererMI::FillRow(){
227   //
228   // fill the output container - 
229   // 2 Options possible
230   //          Tree       
231   //          TObjArray
232   //
233   if (fOutput) fOutput->Fill();
234   if (!fOutput && !fBClonesArray){
235     //
236     if (!fOutputArray) fOutputArray = new TObjArray(fParam->GetNRowsTotal());
237     if (fRowCl && fRowCl->GetArray()->GetEntriesFast()>0) fOutputArray->AddAt(fRowCl->Clone(), fRowCl->GetID());
238   }
239 }
240
241 Float_t  AliTPCclustererMI::GetSigmaY2(Int_t iz){
242   // sigma y2 = in digits  - we don't know the angle
243   Float_t z = iz*fParam->GetZWidth()+fParam->GetNTBinsL1()*fParam->GetZWidth();
244   Float_t sd2 = (z*fParam->GetDiffL()*fParam->GetDiffL())/
245     (fPadWidth*fPadWidth);
246   Float_t sres = 0.25;
247   Float_t res = sd2+sres;
248   return res;
249 }
250
251
252 Float_t  AliTPCclustererMI::GetSigmaZ2(Int_t iz){
253   //sigma z2 = in digits - angle estimated supposing vertex constraint
254   Float_t z = iz*fZWidth+fParam->GetNTBinsL1()*fParam->GetZWidth();
255   Float_t sd2 = (z*fParam->GetDiffL()*fParam->GetDiffL())/(fZWidth*fZWidth);
256   Float_t angular = fPadLength*(fParam->GetZLength(fSector)-z)/(fRx*fZWidth);
257   angular*=angular;
258   angular/=12.;
259   Float_t sres = fParam->GetZSigma()/fZWidth;
260   sres *=sres;
261   Float_t res = angular +sd2+sres;
262   return res;
263 }
264
265 void AliTPCclustererMI::MakeCluster(Int_t k,Int_t max,Float_t *bins, UInt_t /*m*/,
266 AliTPCclusterMI &c) 
267 {
268   //
269   //  k    - Make cluster at position k  
270   //  bins - 2 D array of signals mapped to 1 dimensional array - 
271   //  max  - the number of time bins er one dimension
272   //  c    - refernce to cluster to be filled
273   //
274   Int_t i0=k/max;  //central pad
275   Int_t j0=k%max;  //central time bin
276
277   // set pointers to data
278   //Int_t dummy[5] ={0,0,0,0,0};
279   Float_t * matrix[5]; //5x5 matrix with digits  - indexing i = 0 ..4  j = -2..2
280   for (Int_t di=-2;di<=2;di++){
281     matrix[di+2]  =  &bins[k+di*max];
282   }
283   //build matrix with virtual charge
284   Float_t sigmay2= GetSigmaY2(j0);
285   Float_t sigmaz2= GetSigmaZ2(j0);
286
287   Float_t vmatrix[5][5];
288   vmatrix[2][2] = matrix[2][0];
289   c.SetType(0);
290   c.SetMax((UShort_t)(vmatrix[2][2])); // write maximal amplitude
291   for (Int_t di =-1;di <=1;di++)
292     for (Int_t dj =-1;dj <=1;dj++){
293       Float_t amp = matrix[di+2][dj];
294       if ( (amp<2) && (fLoop<2)){
295         // if under threshold  - calculate virtual charge
296         Float_t ratio = TMath::Exp(-1.2*TMath::Abs(di)/sigmay2)*TMath::Exp(-1.2*TMath::Abs(dj)/sigmaz2);
297         amp = ((matrix[2][0]-2)*(matrix[2][0]-2)/(matrix[-di+2][-dj]+2))*ratio;
298         if (amp>2) amp = 2;
299         vmatrix[2+di][2+dj]=amp;
300         vmatrix[2+2*di][2+2*dj]=0;
301         if ( (di*dj)!=0){       
302           //DIAGONAL ELEMENTS
303           vmatrix[2+2*di][2+dj] =0;
304           vmatrix[2+di][2+2*dj] =0;
305         }
306         continue;
307       }
308       if (amp<4){
309         //if small  amplitude - below  2 x threshold  - don't consider other one        
310         vmatrix[2+di][2+dj]=amp;
311         vmatrix[2+2*di][2+2*dj]=0;  // don't take to the account next bin
312         if ( (di*dj)!=0){       
313           //DIAGONAL ELEMENTS
314           vmatrix[2+2*di][2+dj] =0;
315           vmatrix[2+di][2+2*dj] =0;
316         }
317         continue;
318       }
319       //if bigger then take everything
320       vmatrix[2+di][2+dj]=amp;
321       vmatrix[2+2*di][2+2*dj]= matrix[2*di+2][2*dj] ;      
322       if ( (di*dj)!=0){       
323           //DIAGONAL ELEMENTS
324           vmatrix[2+2*di][2+dj] = matrix[2*di+2][dj];
325           vmatrix[2+di][2+2*dj] = matrix[2+di][dj*2];
326         }      
327     }
328
329
330   
331   Float_t sumw=0;
332   Float_t sumiw=0;
333   Float_t sumi2w=0;
334   Float_t sumjw=0;
335   Float_t sumj2w=0;
336   //
337   for (Int_t i=-2;i<=2;i++)
338     for (Int_t j=-2;j<=2;j++){
339       Float_t amp = vmatrix[i+2][j+2];
340
341       sumw    += amp;
342       sumiw   += i*amp;
343       sumi2w  += i*i*amp;
344       sumjw   += j*amp;
345       sumj2w  += j*j*amp;
346     }    
347   //   
348   Float_t meani = sumiw/sumw;
349   Float_t mi2   = sumi2w/sumw-meani*meani;
350   Float_t meanj = sumjw/sumw;
351   Float_t mj2   = sumj2w/sumw-meanj*meanj;
352   //
353   Float_t ry = mi2/sigmay2;
354   Float_t rz = mj2/sigmaz2;
355   
356   //
357   if ( ( (ry<0.6) || (rz<0.6) ) && fLoop==2) return;
358   if ( (ry <1.2) && (rz<1.2) || (!fRecoParam->GetDoUnfold())) {
359     //
360     //if cluster looks like expected or Unfolding not switched on
361     //standard COG is used
362     //+1.2 deviation from expected sigma accepted
363     //    c.fMax = FitMax(vmatrix,meani,meanj,TMath::Sqrt(sigmay2),TMath::Sqrt(sigmaz2));
364
365     meani +=i0;
366     meanj +=j0;
367     //set cluster parameters
368     c.SetQ(sumw);
369     c.SetPad(meani-2.5);
370     c.SetTimeBin(meanj-3);
371     c.SetSigmaY2(mi2);
372     c.SetSigmaZ2(mj2);
373     c.SetType(0);
374     AddCluster(c,(Float_t*)vmatrix,k);
375     return;     
376   }
377   //
378   //unfolding when neccessary  
379   //
380   
381   Float_t * matrix2[7]; //7x7 matrix with digits  - indexing i = 0 ..6  j = -3..3
382   Float_t dummy[7]={0,0,0,0,0,0};
383   for (Int_t di=-3;di<=3;di++){
384     matrix2[di+3] =  &bins[k+di*max];
385     if ((k+di*max)<3)  matrix2[di+3] = &dummy[3];
386     if ((k+di*max)>fMaxBin-3)  matrix2[di+3] = &dummy[3];
387   }
388   Float_t vmatrix2[5][5];
389   Float_t sumu;
390   Float_t overlap;
391   UnfoldCluster(matrix2,vmatrix2,meani,meanj,sumu,overlap);
392   //
393   //  c.fMax = FitMax(vmatrix2,meani,meanj,TMath::Sqrt(sigmay2),TMath::Sqrt(sigmaz2));
394   meani +=i0;
395   meanj +=j0;
396   //set cluster parameters
397   c.SetQ(sumu);
398   c.SetPad(meani-2.5);
399   c.SetTimeBin(meanj-3);
400   c.SetSigmaY2(mi2);
401   c.SetSigmaZ2(mj2);
402   c.SetType(Char_t(overlap)+1);
403   AddCluster(c,(Float_t*)vmatrix,k);
404
405   //unfolding 2
406   meani-=i0;
407   meanj-=j0;
408 }
409
410
411
412 void AliTPCclustererMI::UnfoldCluster(Float_t * matrix2[7], Float_t recmatrix[5][5], Float_t & meani, Float_t & meanj, 
413                                       Float_t & sumu, Float_t & overlap )
414 {
415   //
416   //unfold cluster from input matrix
417   //data corresponding to cluster writen in recmatrix
418   //output meani and meanj
419
420   //take separatelly y and z
421
422   Float_t sum3i[7] = {0,0,0,0,0,0,0};
423   Float_t sum3j[7] = {0,0,0,0,0,0,0};
424
425   for (Int_t k =0;k<7;k++)
426     for (Int_t l = -1; l<=1;l++){
427       sum3i[k]+=matrix2[k][l];
428       sum3j[k]+=matrix2[l+3][k-3];
429     }
430   Float_t mratio[3][3]={{1,1,1},{1,1,1},{1,1,1}};
431   //
432   //unfold  y 
433   Float_t sum3wi    = 0;  //charge minus overlap
434   Float_t sum3wio   = 0;  //full charge
435   Float_t sum3iw    = 0;  //sum for mean value
436   for (Int_t dk=-1;dk<=1;dk++){
437     sum3wio+=sum3i[dk+3];
438     if (dk==0){
439       sum3wi+=sum3i[dk+3];     
440     }
441     else{
442       Float_t ratio =1;
443       if (  ( ((sum3i[dk+3]+3)/(sum3i[3]-3))+1 < (sum3i[2*dk+3]-3)/(sum3i[dk+3]+3))||
444            sum3i[dk+3]<=sum3i[2*dk+3] && sum3i[dk+3]>2 ){
445         Float_t xm2 = sum3i[-dk+3];
446         Float_t xm1 = sum3i[+3];
447         Float_t x1  = sum3i[2*dk+3];
448         Float_t x2  = sum3i[3*dk+3];    
449         Float_t w11   = TMath::Max((Float_t)(4.*xm1-xm2),(Float_t)0.000001);      
450         Float_t w12   = TMath::Max((Float_t)(4 *x1 -x2),(Float_t)0.);
451         ratio = w11/(w11+w12);   
452         for (Int_t dl=-1;dl<=1;dl++)
453           mratio[dk+1][dl+1] *= ratio;
454       }
455       Float_t amp = sum3i[dk+3]*ratio;
456       sum3wi+=amp;
457       sum3iw+= dk*amp;      
458     }
459   }
460   meani = sum3iw/sum3wi;
461   Float_t overlapi = (sum3wio-sum3wi)/sum3wio;
462
463
464
465   //unfold  z 
466   Float_t sum3wj    = 0;  //charge minus overlap
467   Float_t sum3wjo   = 0;  //full charge
468   Float_t sum3jw    = 0;  //sum for mean value
469   for (Int_t dk=-1;dk<=1;dk++){
470     sum3wjo+=sum3j[dk+3];
471     if (dk==0){
472       sum3wj+=sum3j[dk+3];     
473     }
474     else{
475       Float_t ratio =1;
476       if ( ( ((sum3j[dk+3]+3)/(sum3j[3]-3))+1 < (sum3j[2*dk+3]-3)/(sum3j[dk+3]+3)) ||
477            (sum3j[dk+3]<=sum3j[2*dk+3] && sum3j[dk+3]>2)){
478         Float_t xm2 = sum3j[-dk+3];
479         Float_t xm1 = sum3j[+3];
480         Float_t x1  = sum3j[2*dk+3];
481         Float_t x2  = sum3j[3*dk+3];    
482         Float_t w11   = TMath::Max((Float_t)(4.*xm1-xm2),(Float_t)0.000001);      
483         Float_t w12   = TMath::Max((Float_t)(4 *x1 -x2),(Float_t)0.);
484         ratio = w11/(w11+w12);   
485         for (Int_t dl=-1;dl<=1;dl++)
486           mratio[dl+1][dk+1] *= ratio;
487       }
488       Float_t amp = sum3j[dk+3]*ratio;
489       sum3wj+=amp;
490       sum3jw+= dk*amp;      
491     }
492   }
493   meanj = sum3jw/sum3wj;
494   Float_t overlapj = (sum3wjo-sum3wj)/sum3wjo;  
495   overlap = Int_t(100*TMath::Max(overlapi,overlapj)+3);  
496   sumu = (sum3wj+sum3wi)/2.;
497   
498   if (overlap ==3) {
499     //if not overlap detected remove everything
500     for (Int_t di =-2; di<=2;di++)
501       for (Int_t dj =-2; dj<=2;dj++){
502         recmatrix[di+2][dj+2] = matrix2[3+di][dj];
503       }
504   }
505   else{
506     for (Int_t di =-1; di<=1;di++)
507       for (Int_t dj =-1; dj<=1;dj++){
508         Float_t ratio =1;
509         if (mratio[di+1][dj+1]==1){
510           recmatrix[di+2][dj+2]     = matrix2[3+di][dj];
511           if (TMath::Abs(di)+TMath::Abs(dj)>1){
512             recmatrix[2*di+2][dj+2] = matrix2[3+2*di][dj];
513             recmatrix[di+2][2*dj+2] = matrix2[3+di][2*dj];
514           }       
515           recmatrix[2*di+2][2*dj+2] = matrix2[3+2*di][2*dj];
516         }
517         else
518           {
519             //if we have overlap in direction
520             recmatrix[di+2][dj+2] = mratio[di+1][dj+1]* matrix2[3+di][dj];    
521             if (TMath::Abs(di)+TMath::Abs(dj)>1){
522               ratio =  TMath::Min((Float_t)(recmatrix[di+2][dj+2]/(matrix2[3+0*di][1*dj]+1)),(Float_t)1.);
523               recmatrix[2*di+2][dj+2] = ratio*recmatrix[di+2][dj+2];
524               //
525               ratio =  TMath::Min((Float_t)(recmatrix[di+2][dj+2]/(matrix2[3+1*di][0*dj]+1)),(Float_t)1.);
526               recmatrix[di+2][2*dj+2] = ratio*recmatrix[di+2][dj+2];
527             }
528             else{
529               ratio =  recmatrix[di+2][dj+2]/matrix2[3][0];
530               recmatrix[2*di+2][2*dj+2] = ratio*recmatrix[di+2][dj+2];
531             }
532           }
533       }
534   }
535   
536 }
537
538 Float_t AliTPCclustererMI::FitMax(Float_t vmatrix[5][5], Float_t y, Float_t z, Float_t sigmay, Float_t sigmaz)
539 {
540   //
541   // estimate max
542   Float_t sumteor= 0;
543   Float_t sumamp = 0;
544
545   for (Int_t di = -1;di<=1;di++)
546     for (Int_t dj = -1;dj<=1;dj++){
547       if (vmatrix[2+di][2+dj]>2){
548         Float_t teor = TMath::Gaus(di,y,sigmay*1.2)*TMath::Gaus(dj,z,sigmaz*1.2);
549         sumteor += teor*vmatrix[2+di][2+dj];
550         sumamp  += vmatrix[2+di][2+dj]*vmatrix[2+di][2+dj];
551       }
552     }    
553   Float_t max = sumamp/sumteor;
554   return max;
555 }
556
557 void AliTPCclustererMI::AddCluster(AliTPCclusterMI &c, Float_t * matrix, Int_t pos){
558   //
559   //
560   // Transform cluster to the rotated global coordinata
561   // Assign labels to the cluster
562   // add the cluster to the array
563   // for more details - See  AliTPCTranform::Transform(x,i,0,1) 
564   Float_t meani = c.GetPad();
565   Float_t meanj = c.GetTimeBin();
566
567   Int_t ki = TMath::Nint(meani);
568   if (ki<0) ki=0;
569   if (ki>=fMaxPad) ki = fMaxPad-1;
570   Int_t kj = TMath::Nint(meanj);
571   if (kj<0) kj=0;
572   if (kj>=fMaxTime-3) kj=fMaxTime-4;
573   // ki and kj shifted as integers coordinata
574   if (fRowDig) {
575     c.SetLabel(fRowDig->GetTrackIDFast(kj,ki,0)-2,0);
576     c.SetLabel(fRowDig->GetTrackIDFast(kj,ki,1)-2,1);
577     c.SetLabel(fRowDig->GetTrackIDFast(kj,ki,2)-2,2);
578   }
579   
580   c.SetRow(fRow);
581   c.SetDetector(fSector);
582   Float_t s2 = c.GetSigmaY2();
583   Float_t w=fParam->GetPadPitchWidth(fSector);
584   c.SetSigmaY2(s2*w*w);
585   s2 = c.GetSigmaZ2(); 
586   c.SetSigmaZ2(s2*fZWidth*fZWidth);
587   //
588   //
589   //
590   AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ;
591   if (!transform) {
592     AliFatal("Tranformations not in calibDB");
593   }
594   Double_t x[3]={c.GetRow(),c.GetPad(),c.GetTimeBin()};
595   Int_t i[1]={fSector};
596   transform->Transform(x,i,0,1);
597   c.SetX(x[0]);
598   c.SetY(x[1]);
599   c.SetZ(x[2]);
600   //
601   //
602   if (fSector%36>17){
603     c.SetY(-c.GetY());
604   }
605
606   if (ki<=1 || ki>=fMaxPad-1 || kj==1 || kj==fMaxTime-2) {
607     c.SetType(-(c.GetType()+3));  //edge clusters
608   }
609   if (fLoop==2) c.SetType(100);
610   if (!AcceptCluster(&c)) return;
611
612   // select output 
613   TClonesArray * arr = 0;
614   AliTPCclusterMI * cl = 0;
615
616   if(fBClonesArray==kFALSE) {
617      arr = fRowCl->GetArray();
618      cl = new ((*arr)[fNcluster]) AliTPCclusterMI(c);
619   } else {
620      cl = new ((*fOutputClonesArray)[fNclusters+fNcluster]) AliTPCclusterMI(c);
621   }
622
623   // if (fRecoParam->DumpSignal() &&matrix ) {
624 //     Int_t nbins=0;
625 //     Float_t *graph =0;
626 //     if (fRecoParam->GetCalcPedestal() && cl->GetMax()>fRecoParam->GetDumpAmplitudeMin() &&fBDumpSignal){
627 //       nbins = fMaxTime;
628 //       graph = &(fBins[fMaxTime*(pos/fMaxTime)]);
629 //     }
630 //     AliTPCclusterInfo * info = new AliTPCclusterInfo(matrix,nbins,graph);
631 //     cl->SetInfo(info);
632 //   }
633   if (!fRecoParam->DumpSignal()) {
634     cl->SetInfo(0);
635   }
636   
637   if (AliTPCReconstructor::StreamLevel()>1) {
638      (*fDebugStreamer)<<"Clusters"<<
639        "Cl.="<<cl<<
640        "\n";
641   }
642
643   fNcluster++;
644 }
645
646
647 //_____________________________________________________________________________
648 void AliTPCclustererMI::Digits2Clusters()
649 {
650   //-----------------------------------------------------------------
651   // This is a simple cluster finder.
652   //-----------------------------------------------------------------
653
654   if (!fInput) { 
655     Error("Digits2Clusters", "input tree not initialised");
656     return;
657   }
658  
659   AliTPCCalPad * gainTPC = AliTPCcalibDB::Instance()->GetPadGainFactor();
660   AliTPCCalPad * noiseTPC = AliTPCcalibDB::Instance()->GetPadNoise();
661   AliSimDigits digarr, *dummy=&digarr;
662   fRowDig = dummy;
663   fInput->GetBranch("Segment")->SetAddress(&dummy);
664   Stat_t nentries = fInput->GetEntries();
665   
666   fMaxTime=fRecoParam->GetLastBin()+6; // add 3 virtual time bins before and 3   after
667     
668   Int_t nclusters  = 0;
669
670   for (Int_t n=0; n<nentries; n++) {
671     fInput->GetEvent(n);
672     if (!fParam->AdjustSectorRow(digarr.GetID(),fSector,fRow)) {
673       cerr<<"AliTPC warning: invalid segment ID ! "<<digarr.GetID()<<endl;
674       continue;
675     }
676     Int_t row = fRow;
677     AliTPCCalROC * gainROC = gainTPC->GetCalROC(fSector);  // pad gains per given sector
678     AliTPCCalROC * noiseROC   = noiseTPC->GetCalROC(fSector); // noise per given sector
679     //
680
681     fRowCl->SetID(digarr.GetID());
682     if (fOutput) fOutput->GetBranch("Segment")->SetAddress(&fRowCl);
683     fRx=fParam->GetPadRowRadii(fSector,row);
684     
685     
686     const Int_t kNIS=fParam->GetNInnerSector(), kNOS=fParam->GetNOuterSector();
687     fZWidth = fParam->GetZWidth();
688     if (fSector < kNIS) {
689       fMaxPad = fParam->GetNPadsLow(row);
690       fSign = (fSector < kNIS/2) ? 1 : -1;
691       fPadLength = fParam->GetPadPitchLength(fSector,row);
692       fPadWidth = fParam->GetPadPitchWidth();
693     } else {
694       fMaxPad = fParam->GetNPadsUp(row);
695       fSign = ((fSector-kNIS) < kNOS/2) ? 1 : -1;
696       fPadLength = fParam->GetPadPitchLength(fSector,row);
697       fPadWidth  = fParam->GetPadPitchWidth();
698     }
699     
700     
701     fMaxBin=fMaxTime*(fMaxPad+6);  // add 3 virtual pads  before and 3 after
702     fBins    =new Float_t[fMaxBin];
703     fSigBins =new Int_t[fMaxBin];
704     fNSigBins = 0;
705     memset(fBins,0,sizeof(Float_t)*fMaxBin);
706     
707     if (digarr.First()) //MI change
708       do {
709         Float_t dig=digarr.CurrentDigit();
710         if (dig<=fParam->GetZeroSup()) continue;
711         Int_t j=digarr.CurrentRow()+3, i=digarr.CurrentColumn()+3;
712         Float_t gain = gainROC->GetValue(row,digarr.CurrentColumn());
713         Int_t bin = i*fMaxTime+j;
714         fBins[bin]=dig/gain;
715         fSigBins[fNSigBins++]=bin;
716       } while (digarr.Next());
717     digarr.ExpandTrackBuffer();
718
719     FindClusters(noiseROC);
720     FillRow();
721     fRowCl->GetArray()->Clear();    
722     nclusters+=fNcluster;    
723
724     delete[] fBins;
725     delete[] fSigBins;
726   }  
727
728   Info("Digits2Clusters", "Number of found clusters : %d", nclusters);
729 }
730
731 void AliTPCclustererMI::Digits2Clusters(AliRawReader* rawReader)
732 {
733 //-----------------------------------------------------------------
734 // This is a cluster finder for the TPC raw data.
735 // The method assumes NO ordering of the altro channels.
736 // The pedestal subtraction can be switched on and off
737 // using an option of the TPC reconstructor
738 //-----------------------------------------------------------------
739
740
741   fRowDig = NULL;
742   AliTPCROC * roc = AliTPCROC::Instance();
743   AliTPCCalPad * gainTPC = AliTPCcalibDB::Instance()->GetPadGainFactor();
744   AliTPCCalPad * pedestalTPC = AliTPCcalibDB::Instance()->GetPedestals();
745   AliTPCCalPad * noiseTPC = AliTPCcalibDB::Instance()->GetPadNoise();
746   AliTPCAltroMapping** mapping =AliTPCcalibDB::Instance()->GetMapping();
747   //
748   AliTPCRawStream input(rawReader,(AliAltroMapping**)mapping);
749   fEventHeader = (AliRawEventHeaderBase*)rawReader->GetEventHeader();
750   if (fEventHeader){
751     fTimeStamp = fEventHeader->Get("Timestamp");  
752     fEventType = fEventHeader->Get("Type");  
753   }
754
755   // creaate one TClonesArray for all clusters
756   if(fBClonesArray) fOutputClonesArray = new TClonesArray("AliTPCclusterMI",1000);
757
758   // reset counter
759   fNclusters  = 0;
760   
761   fMaxTime = fRecoParam->GetLastBin() + 6; // add 3 virtual time bins before and 3 after
762   const Int_t kNIS = fParam->GetNInnerSector();
763   const Int_t kNOS = fParam->GetNOuterSector();
764   const Int_t kNS = kNIS + kNOS;
765   fZWidth = fParam->GetZWidth();
766   Int_t zeroSup = fParam->GetZeroSup();
767   //
768   //alocate memory for sector - maximal case
769   //
770   Float_t** allBins = NULL;
771   Int_t** allSigBins = NULL;
772   Int_t*  allNSigBins = NULL;
773   Int_t nRowsMax = roc->GetNRows(roc->GetNSector()-1);
774   Int_t nPadsMax = roc->GetNPads(roc->GetNSector()-1,nRowsMax-1);
775   allBins = new Float_t*[nRowsMax];
776   allSigBins = new Int_t*[nRowsMax];
777   allNSigBins = new Int_t[nRowsMax];
778   for (Int_t iRow = 0; iRow < nRowsMax; iRow++) {
779     //
780     Int_t maxBin = fMaxTime*(nPadsMax+6);  // add 3 virtual pads  before and 3 after
781     allBins[iRow] = new Float_t[maxBin];
782     memset(allBins[iRow],0,sizeof(Float_t)*maxBin);
783     allSigBins[iRow] = new Int_t[maxBin];
784     allNSigBins[iRow]=0;
785   }
786   //
787   // Loop over sectors
788   //
789   for(fSector = 0; fSector < kNS; fSector++) {
790
791     Int_t nRows = 0;
792     Int_t nDDLs = 0, indexDDL = 0;
793     if (fSector < kNIS) {
794       nRows = fParam->GetNRowLow();
795       fSign = (fSector < kNIS/2) ? 1 : -1;
796       nDDLs = 2;
797       indexDDL = fSector * 2;
798     }
799     else {
800       nRows = fParam->GetNRowUp();
801       fSign = ((fSector-kNIS) < kNOS/2) ? 1 : -1;
802       nDDLs = 4;
803       indexDDL = (fSector-kNIS) * 4 + kNIS * 2;
804     }
805
806     // load the raw data for corresponding DDLs
807     rawReader->Reset();
808     rawReader->Select("TPC",indexDDL,indexDDL+nDDLs-1);
809
810     // select only good sector 
811     input.Next();
812     if(input.GetSector() != fSector) continue;
813
814     AliTPCCalROC * gainROC    = gainTPC->GetCalROC(fSector);  // pad gains per given sector
815     AliTPCCalROC * pedestalROC = pedestalTPC->GetCalROC(fSector);  // pedestal per given sector
816     AliTPCCalROC * noiseROC   = noiseTPC->GetCalROC(fSector);  // noise per given sector
817     //check the presence of the calibration
818     if (!noiseROC ||!pedestalROC ) {
819       AliError(Form("Missing calibration per sector\t%d\n",fSector));
820       continue;
821     }
822     
823     for (Int_t iRow = 0; iRow < nRows; iRow++) {
824       Int_t maxPad;
825       if (fSector < kNIS)
826         maxPad = fParam->GetNPadsLow(iRow);
827       else
828         maxPad = fParam->GetNPadsUp(iRow);
829       
830       Int_t maxBin = fMaxTime*(maxPad+6);  // add 3 virtual pads  before and 3 after
831       memset(allBins[iRow],0,sizeof(Float_t)*maxBin);
832       allNSigBins[iRow] = 0;
833     }
834     
835     Int_t digCounter=0;
836     // Begin loop over altro data
837     Bool_t calcPedestal = fRecoParam->GetCalcPedestal();
838     Float_t gain =1;
839     Int_t lastPad=-1;
840
841     input.Reset();
842     while (input.Next()) {
843       if (input.GetSector() != fSector)
844         AliFatal(Form("Sector index mismatch ! Expected (%d), but got (%d) !",fSector,input.GetSector()));
845
846       
847       Int_t iRow = input.GetRow();
848       if (iRow < 0){
849         continue;
850       }
851
852       if (iRow < 0 || iRow >= nRows){
853         AliError(Form("Pad-row index (%d) outside the range (%d -> %d) !",
854                       iRow, 0, nRows -1));
855         continue;
856       }
857       //pad
858       Int_t iPad = input.GetPad();
859       if (iPad < 0 || iPad >= nPadsMax) {
860         AliError(Form("Pad index (%d) outside the range (%d -> %d) !",
861                       iPad, 0, nPadsMax-1));
862         continue;
863       }
864       if (iPad!=lastPad){
865         gain    = gainROC->GetValue(iRow,iPad);
866         lastPad = iPad;
867       }
868       iPad+=3;
869       //time
870       Int_t iTimeBin = input.GetTime();
871       if ( iTimeBin < fRecoParam->GetFirstBin() || iTimeBin >= fRecoParam->GetLastBin()){
872         continue;
873         AliFatal(Form("Timebin index (%d) outside the range (%d -> %d) !",
874                       iTimeBin, 0, iTimeBin -1));
875       }
876       iTimeBin+=3;
877
878       //signal
879       Float_t signal = input.GetSignal();
880       if (!calcPedestal && signal <= zeroSup) continue;      
881
882       if (!calcPedestal) {
883         Int_t bin = iPad*fMaxTime+iTimeBin;
884         allBins[iRow][bin] = signal/gain;
885         allSigBins[iRow][allNSigBins[iRow]++] = bin;
886       }else{
887         allBins[iRow][iPad*fMaxTime+iTimeBin] = signal;
888       }
889       allBins[iRow][iPad*fMaxTime+0]+=1.;  // pad with signal
890
891       // Temporary
892       digCounter++;
893     } // End of the loop over altro data
894     //
895     //
896     //
897     //
898     // Now loop over rows and perform pedestal subtraction
899     if (digCounter==0) continue;
900     //    if (calcPedestal) {
901     if (kTRUE) {
902       for (Int_t iRow = 0; iRow < nRows; iRow++) {
903         Int_t maxPad;
904         if (fSector < kNIS)
905           maxPad = fParam->GetNPadsLow(iRow);
906         else
907           maxPad = fParam->GetNPadsUp(iRow);
908
909         for (Int_t iPad = 3; iPad < maxPad + 3; iPad++) {
910           //
911           // Temporary fix for data production - !!!! MARIAN
912           // The noise calibration should take mean and RMS - currently the Gaussian fit used
913           // In case of double peak  - the pad should be rejected
914           //
915           // Line mean - if more than given digits over threshold - make a noise calculation
916           // and pedestal substration
917           if (!calcPedestal && allBins[iRow][iPad*fMaxTime+0]<50) continue;
918           //
919           if (allBins[iRow][iPad*fMaxTime+0] <1 ) continue;  // no data
920           Float_t *p = &allBins[iRow][iPad*fMaxTime+3];
921           //Float_t pedestal = TMath::Median(fMaxTime, p);      
922           Int_t id[3] = {fSector, iRow, iPad-3};
923           // calib values
924           Double_t rmsCalib=  noiseROC->GetValue(iRow,iPad-3);
925           Double_t pedestalCalib = pedestalROC->GetValue(iRow,iPad-3);
926           Double_t rmsEvent       = rmsCalib;
927           Double_t pedestalEvent  = pedestalCalib;
928           ProcesSignal(p, fMaxTime, id, rmsEvent, pedestalEvent);
929           if (rmsEvent<rmsCalib) rmsEvent = rmsCalib;   // take worst scenario
930           if (TMath::Abs(pedestalEvent-pedestalCalib)<1.0) pedestalEvent = pedestalCalib;  
931           
932           //
933           for (Int_t iTimeBin = 0; iTimeBin < fMaxTime; iTimeBin++) {
934             Int_t bin = iPad*fMaxTime+iTimeBin;
935             allBins[iRow][bin] -= pedestalEvent;
936             if (iTimeBin < AliTPCReconstructor::GetRecoParam()->GetFirstBin())  
937               allBins[iRow][bin] = 0;
938             if (iTimeBin > AliTPCReconstructor::GetRecoParam()->GetLastBin())  
939               allBins[iRow][bin] = 0;
940             if (allBins[iRow][iPad*fMaxTime+iTimeBin] < zeroSup)
941               allBins[iRow][bin] = 0;
942             if (allBins[iRow][bin] < 3.0*rmsEvent)   // 3 sigma cut on RMS
943               allBins[iRow][bin] = 0;
944             if (allBins[iRow][bin]) allSigBins[iRow][allNSigBins[iRow]++] = bin;
945           }
946         }
947       }
948     }
949
950     if (AliTPCReconstructor::StreamLevel()>3) {
951       for (Int_t iRow = 0; iRow < nRows; iRow++) {
952         Int_t maxPad;
953         if (fSector < kNIS)
954           maxPad = fParam->GetNPadsLow(iRow);
955         else
956           maxPad = fParam->GetNPadsUp(iRow);
957         
958         for (Int_t iPad = 3; iPad < maxPad + 3; iPad++) {
959           for (Int_t iTimeBin = 0; iTimeBin < fMaxTime; iTimeBin++) {
960             Int_t bin = iPad*fMaxTime+iTimeBin;
961             Float_t signal = allBins[iRow][bin];
962             if (AliTPCReconstructor::StreamLevel()>3 && signal>3) {
963               Double_t x[]={iRow,iPad-3,iTimeBin-3};
964               Int_t i[]={fSector};
965               AliTPCTransform trafo;
966               trafo.Transform(x,i,0,1);
967               Double_t gx[3]={x[0],x[1],x[2]};
968               trafo.RotatedGlobal2Global(fSector,gx);
969               
970               (*fDebugStreamer)<<"Digits"<<
971                 "sec="<<fSector<<
972                 "row="<<iRow<<
973                 "pad="<<iPad<<
974                 "time="<<iTimeBin<<
975                 "sig="<<signal<<
976                 "x="<<x[0]<<
977                 "y="<<x[1]<<
978                 "z="<<x[2]<<      
979                 "gx="<<gx[0]<<
980                 "gy="<<gx[1]<<
981                 "gz="<<gx[2]<<
982                 "\n";
983             }
984           }
985         }
986       }
987     }
988
989     // Now loop over rows and find clusters
990     for (fRow = 0; fRow < nRows; fRow++) {
991       fRowCl->SetID(fParam->GetIndex(fSector, fRow));
992       if (fOutput) fOutput->GetBranch("Segment")->SetAddress(&fRowCl);
993
994       fRx = fParam->GetPadRowRadii(fSector, fRow);
995       fPadLength = fParam->GetPadPitchLength(fSector, fRow);
996       fPadWidth  = fParam->GetPadPitchWidth();
997       if (fSector < kNIS)
998         fMaxPad = fParam->GetNPadsLow(fRow);
999       else
1000         fMaxPad = fParam->GetNPadsUp(fRow);
1001       fMaxBin = fMaxTime*(fMaxPad+6);  // add 3 virtual pads  before and 3 after
1002
1003       fBins = allBins[fRow];
1004       fSigBins = allSigBins[fRow];
1005       fNSigBins = allNSigBins[fRow];
1006
1007       FindClusters(noiseROC);
1008
1009       FillRow();
1010       if(fBClonesArray == kFALSE) fRowCl->GetArray()->Clear();    
1011       fNclusters += fNcluster;    
1012
1013     } // End of loop to find clusters
1014   } // End of loop over sectors
1015
1016   for (Int_t iRow = 0; iRow < nRowsMax; iRow++) {
1017     delete [] allBins[iRow];
1018     delete [] allSigBins[iRow];
1019   }  
1020   delete [] allBins;
1021   delete [] allSigBins;
1022   delete [] allNSigBins;
1023   
1024   if (rawReader->GetEventId() && fOutput ){
1025     Info("Digits2Clusters", "File  %s Event\t%d\tNumber of found clusters : %d\n", fOutput->GetName(),*(rawReader->GetEventId()), fNclusters);
1026   }
1027   
1028   if(rawReader->GetEventId()) {
1029     Info("Digits2Clusters", "Event\t%d\tNumber of found clusters : %d\n",*(rawReader->GetEventId()), fNclusters);
1030   }
1031
1032   if(fBClonesArray) {
1033     //Info("Digits2Clusters", "Number of found clusters : %d\n",fOutputClonesArray->GetEntriesFast());
1034   }
1035 }
1036
1037 void AliTPCclustererMI::FindClusters(AliTPCCalROC * noiseROC)
1038 {
1039   
1040   //
1041   // add virtual charge at the edge   
1042   //
1043   Double_t kMaxDumpSize = 500000;
1044   if (!fOutput) {
1045     fBDumpSignal =kFALSE;
1046   }else{
1047     if (fRecoParam->GetCalcPedestal() && fOutput->GetZipBytes()< kMaxDumpSize) fBDumpSignal =kTRUE;   //dump signal flag
1048   }
1049    
1050   fNcluster=0;
1051   fLoop=1;
1052   Int_t crtime = Int_t((fParam->GetZLength(fSector)-fRecoParam->GetCtgRange()*fRx)/fZWidth-fParam->GetNTBinsL1()-5);
1053   Float_t minMaxCutAbs       = fRecoParam->GetMinMaxCutAbs();
1054   Float_t minLeftRightCutAbs = fRecoParam->GetMinLeftRightCutAbs();
1055   Float_t minUpDownCutAbs    = fRecoParam->GetMinUpDownCutAbs();
1056   Float_t minMaxCutSigma       = fRecoParam->GetMinMaxCutSigma();
1057   Float_t minLeftRightCutSigma = fRecoParam->GetMinLeftRightCutSigma();
1058   Float_t minUpDownCutSigma    = fRecoParam->GetMinUpDownCutSigma();
1059   for (Int_t iSig = 0; iSig < fNSigBins; iSig++) {
1060     Int_t i = fSigBins[iSig];
1061     if (i%fMaxTime<=crtime) continue;
1062     Float_t *b = &fBins[i];
1063     //absolute custs
1064     if (b[0]<minMaxCutAbs) continue;   //threshold for maxima  
1065     //
1066     if (b[-1]+b[1]+b[-fMaxTime]+b[fMaxTime]<=0) continue;  // cut on isolated clusters 
1067     if (b[-1]+b[1]<=0) continue;               // cut on isolated clusters
1068     if (b[-fMaxTime]+b[fMaxTime]<=0) continue; // cut on isolated clusters
1069     //
1070     if ((b[0]+b[-1]+b[1])<minUpDownCutAbs) continue;   //threshold for up down  (TRF) 
1071     if ((b[0]+b[-fMaxTime]+b[fMaxTime])<minLeftRightCutAbs) continue;   //threshold for left right (PRF)    
1072     if (!IsMaximum(*b,fMaxTime,b)) continue;
1073     //
1074     Float_t noise = noiseROC->GetValue(fRow, i/fMaxTime);
1075     if (noise>fRecoParam->GetMaxNoise()) continue;
1076     // sigma cuts
1077     if (b[0]<minMaxCutSigma*noise) continue;   //threshold form maxima  
1078     if ((b[0]+b[-1]+b[1])<minUpDownCutSigma*noise) continue;   //threshold for up town TRF 
1079     if ((b[0]+b[-fMaxTime]+b[fMaxTime])<minLeftRightCutSigma*noise) continue;   //threshold for left right (PRF)    
1080   
1081     AliTPCclusterMI c;   // default cosntruction  without info
1082     Int_t dummy=0;
1083     MakeCluster(i, fMaxTime, fBins, dummy,c);
1084     
1085     //}
1086   }
1087 }
1088
1089 Bool_t AliTPCclustererMI::AcceptCluster(AliTPCclusterMI *cl){
1090   //
1091   // Currently hack to filter digital noise (15.06.2008)
1092   // To be parameterized in the AliTPCrecoParam
1093   // More inteligent way  to be used in future
1094   // Acces to the proper pedestal file needed
1095   //
1096   if (cl->GetMax()<400) return kTRUE;
1097   Double_t ratio = cl->GetQ()/cl->GetMax();
1098   if (cl->GetMax()>700){
1099     if ((ratio - int(ratio)>0.8)) return kFALSE;
1100   }
1101   if ((ratio - int(ratio)<0.95)) return kTRUE;
1102   return kFALSE;
1103 }
1104
1105
1106 Double_t AliTPCclustererMI::ProcesSignal(Float_t *signal, Int_t nchannels, Int_t id[3], Double_t &rmsEvent, Double_t &pedestalEvent){
1107   //
1108   // process signal on given pad - + streaming of additional information in special mode
1109   //
1110   // id[0] - sector
1111   // id[1] - row
1112   // id[2] - pad 
1113
1114   //
1115   // ESTIMATE pedestal and the noise
1116   // 
1117   const Int_t kPedMax = 100;
1118   Float_t  max    =  0;
1119   Float_t  maxPos =  0;
1120   Int_t    median =  -1;
1121   Int_t    count0 =  0;
1122   Int_t    count1 =  0;
1123   Float_t  rmsCalib   = rmsEvent;       // backup initial value ( from calib)
1124   Float_t  pedestalCalib = pedestalEvent;// backup initial value ( from calib)
1125   Int_t    firstBin = AliTPCReconstructor::GetRecoParam()->GetFirstBin();
1126   //
1127   UShort_t histo[kPedMax];
1128   //memset(histo,0,kPedMax*sizeof(UShort_t));
1129   for (Int_t i=0; i<kPedMax; i++) histo[i]=0;
1130   for (Int_t i=0; i<fMaxTime; i++){
1131     if (signal[i]<=0) continue;
1132     if (signal[i]>max && i>firstBin) {
1133       max = signal[i];
1134       maxPos = i;
1135     }
1136     if (signal[i]>kPedMax-1) continue;
1137     histo[int(signal[i]+0.5)]++;
1138     count0++;
1139   }
1140   //
1141   for (Int_t i=1; i<kPedMax; i++){
1142     if (count1<count0*0.5) median=i;
1143     count1+=histo[i];
1144   }
1145   // truncated mean  
1146   //
1147   Float_t count10=histo[median] ,mean=histo[median]*median,  rms=histo[median]*median*median ;
1148   Float_t count06=histo[median] ,mean06=histo[median]*median,  rms06=histo[median]*median*median ;
1149   Float_t count09=histo[median] ,mean09=histo[median]*median,  rms09=histo[median]*median*median ;
1150   //
1151   for (Int_t idelta=1; idelta<10; idelta++){
1152     if (median-idelta<=0) continue;
1153     if (median+idelta>kPedMax) continue;
1154     if (count06<0.6*count1){
1155       count06+=histo[median-idelta];
1156       mean06 +=histo[median-idelta]*(median-idelta);
1157       rms06  +=histo[median-idelta]*(median-idelta)*(median-idelta);
1158       count06+=histo[median+idelta];
1159       mean06 +=histo[median+idelta]*(median+idelta);
1160       rms06  +=histo[median+idelta]*(median+idelta)*(median+idelta);
1161     }
1162     if (count09<0.9*count1){
1163       count09+=histo[median-idelta];
1164       mean09 +=histo[median-idelta]*(median-idelta);
1165       rms09  +=histo[median-idelta]*(median-idelta)*(median-idelta);
1166       count09+=histo[median+idelta];
1167       mean09 +=histo[median+idelta]*(median+idelta);
1168       rms09  +=histo[median+idelta]*(median+idelta)*(median+idelta);
1169     }
1170     if (count10<0.95*count1){
1171       count10+=histo[median-idelta];
1172       mean +=histo[median-idelta]*(median-idelta);
1173       rms  +=histo[median-idelta]*(median-idelta)*(median-idelta);
1174       count10+=histo[median+idelta];
1175       mean +=histo[median+idelta]*(median+idelta);
1176       rms  +=histo[median+idelta]*(median+idelta)*(median+idelta);
1177     }
1178   }
1179   if (count10) {
1180     mean  /=count10;
1181     rms    = TMath::Sqrt(TMath::Abs(rms/count10-mean*mean));
1182   }
1183   if (count06) {
1184     mean06/=count06;
1185     rms06    = TMath::Sqrt(TMath::Abs(rms06/count06-mean06*mean06));
1186   }
1187   if (count09) {
1188     mean09/=count09;
1189     rms09    = TMath::Sqrt(TMath::Abs(rms09/count09-mean09*mean09));
1190   }
1191   rmsEvent = rms09;
1192   //
1193   pedestalEvent = median;
1194   if (AliLog::GetDebugLevel("","AliTPCclustererMI")==0) return median;
1195   //
1196   UInt_t uid[3] = {UInt_t(id[0]),UInt_t(id[1]),UInt_t(id[2])};
1197   //
1198   // Dump mean signal info
1199   //
1200   (*fDebugStreamer)<<"Signal"<<
1201     "TimeStamp="<<fTimeStamp<<
1202     "EventType="<<fEventType<<
1203     "Sector="<<uid[0]<<
1204     "Row="<<uid[1]<<
1205     "Pad="<<uid[2]<<
1206     "Max="<<max<<
1207     "MaxPos="<<maxPos<<
1208     //
1209     "Median="<<median<<
1210     "Mean="<<mean<<
1211     "RMS="<<rms<<      
1212     "Mean06="<<mean06<<
1213     "RMS06="<<rms06<<
1214     "Mean09="<<mean09<<
1215     "RMS09="<<rms09<<
1216     "RMSCalib="<<rmsCalib<<
1217     "PedCalib="<<pedestalCalib<<
1218     "\n";
1219   //
1220   // fill pedestal histogram
1221   //
1222   //
1223   //
1224   //
1225   Float_t kMin =fRecoParam->GetDumpAmplitudeMin();   // minimal signal to be dumped
1226   Float_t *dsignal = new Float_t[nchannels];
1227   Float_t *dtime   = new Float_t[nchannels];
1228   for (Int_t i=0; i<nchannels; i++){
1229     dtime[i] = i;
1230     dsignal[i] = signal[i];
1231   }
1232
1233   TGraph * graph=0;
1234   //
1235   // Big signals dumping
1236   //    
1237   if (max-median>kMin &&maxPos>AliTPCReconstructor::GetRecoParam()->GetFirstBin()) 
1238     (*fDebugStreamer)<<"SignalB"<<     // pads with signal
1239       "TimeStamp="<<fTimeStamp<<
1240       "EventType="<<fEventType<<
1241       "Sector="<<uid[0]<<
1242       "Row="<<uid[1]<<
1243       "Pad="<<uid[2]<<
1244       "Graph="<<graph<<
1245       "Max="<<max<<
1246       "MaxPos="<<maxPos<<       
1247       //
1248       "Median="<<median<<
1249       "Mean="<<mean<<
1250       "RMS="<<rms<<      
1251       "Mean06="<<mean06<<
1252       "RMS06="<<rms06<<
1253       "Mean09="<<mean09<<
1254       "RMS09="<<rms09<<
1255       "\n";
1256   delete graph;  
1257
1258   delete [] dsignal;
1259   delete [] dtime;
1260   if (rms06>fRecoParam->GetMaxNoise()) {
1261     pedestalEvent+=1024.;
1262     return 1024+median; // sign noisy channel in debug mode
1263   }
1264   return median;
1265 }
1266
1267
1268
1269
1270