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