Removing gcc 4.3 warnings
[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     Float_t xyz[3];
639     cl->GetGlobalXYZ(xyz);
640      (*fDebugStreamer)<<"Clusters"<<
641        "Cl.="<<cl<<
642        "gx="<<xyz[0]<<
643        "gy="<<xyz[1]<<
644        "gz="<<xyz[2]<<
645        "\n";
646   }
647
648   fNcluster++;
649 }
650
651
652 //_____________________________________________________________________________
653 void AliTPCclustererMI::Digits2Clusters()
654 {
655   //-----------------------------------------------------------------
656   // This is a simple cluster finder.
657   //-----------------------------------------------------------------
658
659   if (!fInput) { 
660     Error("Digits2Clusters", "input tree not initialised");
661     return;
662   }
663   fRecoParam = AliTPCReconstructor::GetRecoParam();
664   if (!fRecoParam){
665     AliFatal("Can not get the reconstruction parameters");
666   }
667   if(AliTPCReconstructor::StreamLevel()>5) {
668     AliInfo("Parameter Dumps");
669     fParam->Dump();
670     fRecoParam->Dump();
671   }
672
673   AliTPCCalPad * gainTPC = AliTPCcalibDB::Instance()->GetPadGainFactor();
674   AliTPCCalPad * noiseTPC = AliTPCcalibDB::Instance()->GetPadNoise();
675   AliSimDigits digarr, *dummy=&digarr;
676   fRowDig = dummy;
677   fInput->GetBranch("Segment")->SetAddress(&dummy);
678   Stat_t nentries = fInput->GetEntries();
679   
680   fMaxTime=fRecoParam->GetLastBin()+6; // add 3 virtual time bins before and 3   after
681     
682   Int_t nclusters  = 0;
683
684   for (Int_t n=0; n<nentries; n++) {
685     fInput->GetEvent(n);
686     if (!fParam->AdjustSectorRow(digarr.GetID(),fSector,fRow)) {
687       cerr<<"AliTPC warning: invalid segment ID ! "<<digarr.GetID()<<endl;
688       continue;
689     }
690     Int_t row = fRow;
691     AliTPCCalROC * gainROC = gainTPC->GetCalROC(fSector);  // pad gains per given sector
692     AliTPCCalROC * noiseROC   = noiseTPC->GetCalROC(fSector); // noise per given sector
693     //
694
695     fRowCl->SetID(digarr.GetID());
696     if (fOutput) fOutput->GetBranch("Segment")->SetAddress(&fRowCl);
697     fRx=fParam->GetPadRowRadii(fSector,row);
698     
699     
700     const Int_t kNIS=fParam->GetNInnerSector(), kNOS=fParam->GetNOuterSector();
701     fZWidth = fParam->GetZWidth();
702     if (fSector < kNIS) {
703       fMaxPad = fParam->GetNPadsLow(row);
704       fSign = (fSector < kNIS/2) ? 1 : -1;
705       fPadLength = fParam->GetPadPitchLength(fSector,row);
706       fPadWidth = fParam->GetPadPitchWidth();
707     } else {
708       fMaxPad = fParam->GetNPadsUp(row);
709       fSign = ((fSector-kNIS) < kNOS/2) ? 1 : -1;
710       fPadLength = fParam->GetPadPitchLength(fSector,row);
711       fPadWidth  = fParam->GetPadPitchWidth();
712     }
713     
714     
715     fMaxBin=fMaxTime*(fMaxPad+6);  // add 3 virtual pads  before and 3 after
716     fBins    =new Float_t[fMaxBin];
717     fSigBins =new Int_t[fMaxBin];
718     fNSigBins = 0;
719     memset(fBins,0,sizeof(Float_t)*fMaxBin);
720     
721     if (digarr.First()) //MI change
722       do {
723         Float_t dig=digarr.CurrentDigit();
724         if (dig<=fParam->GetZeroSup()) continue;
725         Int_t j=digarr.CurrentRow()+3, i=digarr.CurrentColumn()+3;
726         Float_t gain = gainROC->GetValue(row,digarr.CurrentColumn());
727         Int_t bin = i*fMaxTime+j;
728         if (gain>0){
729           fBins[bin]=dig/gain;
730         }else{
731           fBins[bin]=0;
732         }
733         fSigBins[fNSigBins++]=bin;
734       } while (digarr.Next());
735     digarr.ExpandTrackBuffer();
736
737     FindClusters(noiseROC);
738     FillRow();
739     fRowCl->GetArray()->Clear();    
740     nclusters+=fNcluster;    
741
742     delete[] fBins;
743     delete[] fSigBins;
744   }  
745  
746   Info("Digits2Clusters", "Number of found clusters : %d", nclusters);
747 }
748
749 void AliTPCclustererMI::Digits2Clusters(AliRawReader* rawReader)
750 {
751 //-----------------------------------------------------------------
752 // This is a cluster finder for the TPC raw data.
753 // The method assumes NO ordering of the altro channels.
754 // The pedestal subtraction can be switched on and off
755 // using an option of the TPC reconstructor
756 //-----------------------------------------------------------------
757   fRecoParam = AliTPCReconstructor::GetRecoParam();
758   if (!fRecoParam){
759     AliFatal("Can not get the reconstruction parameters");
760   }
761   if(AliTPCReconstructor::StreamLevel()>5) {
762     AliInfo("Parameter Dumps");
763     fParam->Dump();
764     fRecoParam->Dump();
765   }
766   fRowDig = NULL;
767   AliTPCROC * roc = AliTPCROC::Instance();
768   AliTPCCalPad * gainTPC = AliTPCcalibDB::Instance()->GetPadGainFactor();
769   AliTPCCalPad * pedestalTPC = AliTPCcalibDB::Instance()->GetPedestals();
770   AliTPCCalPad * noiseTPC = AliTPCcalibDB::Instance()->GetPadNoise();
771   AliTPCAltroMapping** mapping =AliTPCcalibDB::Instance()->GetMapping();
772   //
773   AliTPCRawStream input(rawReader,(AliAltroMapping**)mapping);
774   fEventHeader = (AliRawEventHeaderBase*)rawReader->GetEventHeader();
775   if (fEventHeader){
776     fTimeStamp = fEventHeader->Get("Timestamp");  
777     fEventType = fEventHeader->Get("Type");  
778   }
779
780   // creaate one TClonesArray for all clusters
781   if(fBClonesArray && !fOutputClonesArray) fOutputClonesArray = new TClonesArray("AliTPCclusterMI",1000);
782   // reset counter
783   fNclusters  = 0;
784   
785   fMaxTime = fRecoParam->GetLastBin() + 6; // add 3 virtual time bins before and 3 after
786   const Int_t kNIS = fParam->GetNInnerSector();
787   const Int_t kNOS = fParam->GetNOuterSector();
788   const Int_t kNS = kNIS + kNOS;
789   fZWidth = fParam->GetZWidth();
790   Int_t zeroSup = fParam->GetZeroSup();
791   //
792   //alocate memory for sector - maximal case
793   //
794   Float_t** allBins = NULL;
795   Int_t** allSigBins = NULL;
796   Int_t*  allNSigBins = NULL;
797   Int_t nRowsMax = roc->GetNRows(roc->GetNSector()-1);
798   Int_t nPadsMax = roc->GetNPads(roc->GetNSector()-1,nRowsMax-1);
799   allBins = new Float_t*[nRowsMax];
800   allSigBins = new Int_t*[nRowsMax];
801   allNSigBins = new Int_t[nRowsMax];
802   for (Int_t iRow = 0; iRow < nRowsMax; iRow++) {
803     //
804     Int_t maxBin = fMaxTime*(nPadsMax+6);  // add 3 virtual pads  before and 3 after
805     allBins[iRow] = new Float_t[maxBin];
806     memset(allBins[iRow],0,sizeof(Float_t)*maxBin);
807     allSigBins[iRow] = new Int_t[maxBin];
808     allNSigBins[iRow]=0;
809   }
810   //
811   // Loop over sectors
812   //
813   for(fSector = 0; fSector < kNS; fSector++) {
814
815     Int_t nRows = 0;
816     Int_t nDDLs = 0, indexDDL = 0;
817     if (fSector < kNIS) {
818       nRows = fParam->GetNRowLow();
819       fSign = (fSector < kNIS/2) ? 1 : -1;
820       nDDLs = 2;
821       indexDDL = fSector * 2;
822     }
823     else {
824       nRows = fParam->GetNRowUp();
825       fSign = ((fSector-kNIS) < kNOS/2) ? 1 : -1;
826       nDDLs = 4;
827       indexDDL = (fSector-kNIS) * 4 + kNIS * 2;
828     }
829
830     // load the raw data for corresponding DDLs
831     rawReader->Reset();
832     rawReader->Select("TPC",indexDDL,indexDDL+nDDLs-1);
833
834     // select only good sector 
835     input.Next();
836     if(input.GetSector() != fSector) continue;
837
838     AliTPCCalROC * gainROC    = gainTPC->GetCalROC(fSector);  // pad gains per given sector
839     AliTPCCalROC * pedestalROC = pedestalTPC->GetCalROC(fSector);  // pedestal per given sector
840     AliTPCCalROC * noiseROC   = noiseTPC->GetCalROC(fSector);  // noise per given sector
841     //check the presence of the calibration
842     if (!noiseROC ||!pedestalROC ) {
843       AliError(Form("Missing calibration per sector\t%d\n",fSector));
844       continue;
845     }
846     
847     for (Int_t iRow = 0; iRow < nRows; iRow++) {
848       Int_t maxPad;
849       if (fSector < kNIS)
850         maxPad = fParam->GetNPadsLow(iRow);
851       else
852         maxPad = fParam->GetNPadsUp(iRow);
853       
854       Int_t maxBin = fMaxTime*(maxPad+6);  // add 3 virtual pads  before and 3 after
855       memset(allBins[iRow],0,sizeof(Float_t)*maxBin);
856       allNSigBins[iRow] = 0;
857     }
858     
859     Int_t digCounter=0;
860     // Begin loop over altro data
861     Bool_t calcPedestal = fRecoParam->GetCalcPedestal();
862     Float_t gain =1;
863     Int_t lastPad=-1;
864
865     input.Reset();
866     while (input.Next()) {
867       if (input.GetSector() != fSector)
868         AliFatal(Form("Sector index mismatch ! Expected (%d), but got (%d) !",fSector,input.GetSector()));
869
870       
871       Int_t iRow = input.GetRow();
872       if (iRow < 0){
873         continue;
874       }
875
876       if (iRow < 0 || iRow >= nRows){
877         AliError(Form("Pad-row index (%d) outside the range (%d -> %d) !",
878                       iRow, 0, nRows -1));
879         continue;
880       }
881       //pad
882       Int_t iPad = input.GetPad();
883       if (iPad < 0 || iPad >= nPadsMax) {
884         AliError(Form("Pad index (%d) outside the range (%d -> %d) !",
885                       iPad, 0, nPadsMax-1));
886         continue;
887       }
888       if (iPad!=lastPad){
889         gain    = gainROC->GetValue(iRow,iPad);
890         lastPad = iPad;
891       }
892       iPad+=3;
893       //time
894       Int_t iTimeBin = input.GetTime();
895       if ( iTimeBin < fRecoParam->GetFirstBin() || iTimeBin >= fRecoParam->GetLastBin()){
896         continue;
897         AliFatal(Form("Timebin index (%d) outside the range (%d -> %d) !",
898                       iTimeBin, 0, iTimeBin -1));
899       }
900       iTimeBin+=3;
901
902       //signal
903       Float_t signal = input.GetSignal();
904       if (!calcPedestal && signal <= zeroSup) continue;      
905
906       if (!calcPedestal) {
907         Int_t bin = iPad*fMaxTime+iTimeBin;
908         if (gain>0){
909           allBins[iRow][bin] = signal/gain;
910         }else{
911           allBins[iRow][bin] =0;
912         }
913         allSigBins[iRow][allNSigBins[iRow]++] = bin;
914       }else{
915         allBins[iRow][iPad*fMaxTime+iTimeBin] = signal;
916       }
917       allBins[iRow][iPad*fMaxTime+0]+=1.;  // pad with signal
918
919       // Temporary
920       digCounter++;
921     } // End of the loop over altro data
922     //
923     //
924     //
925     //
926     // Now loop over rows and perform pedestal subtraction
927     if (digCounter==0) continue;
928     //    if (calcPedestal) {
929     if (kFALSE ) {
930       for (Int_t iRow = 0; iRow < nRows; iRow++) {
931         Int_t maxPad;
932         if (fSector < kNIS)
933           maxPad = fParam->GetNPadsLow(iRow);
934         else
935           maxPad = fParam->GetNPadsUp(iRow);
936
937         for (Int_t iPad = 3; iPad < maxPad + 3; iPad++) {
938           //
939           // Temporary fix for data production - !!!! MARIAN
940           // The noise calibration should take mean and RMS - currently the Gaussian fit used
941           // In case of double peak  - the pad should be rejected
942           //
943           // Line mean - if more than given digits over threshold - make a noise calculation
944           // and pedestal substration
945           if (!calcPedestal && allBins[iRow][iPad*fMaxTime+0]<50) continue;
946           //
947           if (allBins[iRow][iPad*fMaxTime+0] <1 ) continue;  // no data
948           Float_t *p = &allBins[iRow][iPad*fMaxTime+3];
949           //Float_t pedestal = TMath::Median(fMaxTime, p);      
950           Int_t id[3] = {fSector, iRow, iPad-3};
951           // calib values
952           Double_t rmsCalib=  noiseROC->GetValue(iRow,iPad-3);
953           Double_t pedestalCalib = pedestalROC->GetValue(iRow,iPad-3);
954           Double_t rmsEvent       = rmsCalib;
955           Double_t pedestalEvent  = pedestalCalib;
956           ProcesSignal(p, fMaxTime, id, rmsEvent, pedestalEvent);
957           if (rmsEvent<rmsCalib) rmsEvent = rmsCalib;   // take worst scenario
958           if (TMath::Abs(pedestalEvent-pedestalCalib)<1.0) pedestalEvent = pedestalCalib;  
959           
960           //
961           for (Int_t iTimeBin = 0; iTimeBin < fMaxTime; iTimeBin++) {
962             Int_t bin = iPad*fMaxTime+iTimeBin;
963             allBins[iRow][bin] -= pedestalEvent;
964             if (iTimeBin < fRecoParam->GetFirstBin())  
965               allBins[iRow][bin] = 0;
966             if (iTimeBin > fRecoParam->GetLastBin())  
967               allBins[iRow][bin] = 0;
968             if (allBins[iRow][iPad*fMaxTime+iTimeBin] < zeroSup)
969               allBins[iRow][bin] = 0;
970             if (allBins[iRow][bin] < 3.0*rmsEvent)   // 3 sigma cut on RMS
971               allBins[iRow][bin] = 0;
972             if (allBins[iRow][bin]) allSigBins[iRow][allNSigBins[iRow]++] = bin;
973           }
974         }
975       }
976     }
977
978     if (AliTPCReconstructor::StreamLevel()>3) {
979       for (Int_t iRow = 0; iRow < nRows; iRow++) {
980         Int_t maxPad;
981         if (fSector < kNIS)
982           maxPad = fParam->GetNPadsLow(iRow);
983         else
984           maxPad = fParam->GetNPadsUp(iRow);
985         
986         for (Int_t iPad = 3; iPad < maxPad + 3; iPad++) {
987           for (Int_t iTimeBin = 0; iTimeBin < fMaxTime; iTimeBin++) {
988             Int_t bin = iPad*fMaxTime+iTimeBin;
989             Float_t signal = allBins[iRow][bin];
990             if (AliTPCReconstructor::StreamLevel()>3 && signal>3) {
991               Double_t x[]={iRow,iPad-3,iTimeBin-3};
992               Int_t i[]={fSector};
993               AliTPCTransform trafo;
994               trafo.Transform(x,i,0,1);
995               Double_t gx[3]={x[0],x[1],x[2]};
996               trafo.RotatedGlobal2Global(fSector,gx);
997               //              allSigBins[iRow][allNSigBins[iRow]++]
998               Int_t rowsigBins = allNSigBins[iRow];
999               Int_t first=allSigBins[iRow][0];
1000               Int_t last= 0;
1001               //              if (rowsigBins>0) allSigBins[iRow][allNSigBins[iRow]-1];
1002
1003               if (AliTPCReconstructor::StreamLevel()>0) {
1004               (*fDebugStreamer)<<"Digits"<<
1005                 "sec="<<fSector<<
1006                 "row="<<iRow<<
1007                 "pad="<<iPad<<
1008                 "time="<<iTimeBin<<
1009                 "sig="<<signal<<
1010                 "x="<<x[0]<<
1011                 "y="<<x[1]<<
1012                 "z="<<x[2]<<      
1013                 "gx="<<gx[0]<<
1014                 "gy="<<gx[1]<<
1015                 "gz="<<gx[2]<<
1016                 //
1017                 "rowsigBins="<<rowsigBins<<
1018                 "first="<<first<<
1019                 "last="<<last<<
1020                 "\n";
1021               }
1022             }
1023           }
1024         }
1025       }
1026     }
1027
1028     // Now loop over rows and find clusters
1029     for (fRow = 0; fRow < nRows; fRow++) {
1030       fRowCl->SetID(fParam->GetIndex(fSector, fRow));
1031       if (fOutput) fOutput->GetBranch("Segment")->SetAddress(&fRowCl);
1032
1033       fRx = fParam->GetPadRowRadii(fSector, fRow);
1034       fPadLength = fParam->GetPadPitchLength(fSector, fRow);
1035       fPadWidth  = fParam->GetPadPitchWidth();
1036       if (fSector < kNIS)
1037         fMaxPad = fParam->GetNPadsLow(fRow);
1038       else
1039         fMaxPad = fParam->GetNPadsUp(fRow);
1040       fMaxBin = fMaxTime*(fMaxPad+6);  // add 3 virtual pads  before and 3 after
1041
1042       fBins = allBins[fRow];
1043       fSigBins = allSigBins[fRow];
1044       fNSigBins = allNSigBins[fRow];
1045
1046       FindClusters(noiseROC);
1047
1048       FillRow();
1049       if(fBClonesArray == kFALSE) fRowCl->GetArray()->Clear();    
1050       fNclusters += fNcluster;    
1051
1052     } // End of loop to find clusters
1053   } // End of loop over sectors
1054
1055   for (Int_t iRow = 0; iRow < nRowsMax; iRow++) {
1056     delete [] allBins[iRow];
1057     delete [] allSigBins[iRow];
1058   }  
1059   delete [] allBins;
1060   delete [] allSigBins;
1061   delete [] allNSigBins;
1062   
1063   if (rawReader->GetEventId() && fOutput ){
1064     Info("Digits2Clusters", "File  %s Event\t%d\tNumber of found clusters : %d\n", fOutput->GetName(),*(rawReader->GetEventId()), fNclusters);
1065   }
1066   
1067   if(rawReader->GetEventId()) {
1068     Info("Digits2Clusters", "Event\t%d\tNumber of found clusters : %d\n",*(rawReader->GetEventId()), fNclusters);
1069   }
1070
1071   if(fBClonesArray) {
1072     //Info("Digits2Clusters", "Number of found clusters : %d\n",fOutputClonesArray->GetEntriesFast());
1073   }
1074 }
1075
1076 void AliTPCclustererMI::FindClusters(AliTPCCalROC * noiseROC)
1077 {
1078   
1079   //
1080   // add virtual charge at the edge   
1081   //
1082   Double_t kMaxDumpSize = 500000;
1083   if (!fOutput) {
1084     fBDumpSignal =kFALSE;
1085   }else{
1086     if (fRecoParam->GetCalcPedestal() && fOutput->GetZipBytes()< kMaxDumpSize) fBDumpSignal =kTRUE;   //dump signal flag
1087   }
1088    
1089   fNcluster=0;
1090   fLoop=1;
1091   Int_t crtime = Int_t((fParam->GetZLength(fSector)-fRecoParam->GetCtgRange()*fRx)/fZWidth-fParam->GetNTBinsL1()-5);
1092   Float_t minMaxCutAbs       = fRecoParam->GetMinMaxCutAbs();
1093   Float_t minLeftRightCutAbs = fRecoParam->GetMinLeftRightCutAbs();
1094   Float_t minUpDownCutAbs    = fRecoParam->GetMinUpDownCutAbs();
1095   Float_t minMaxCutSigma       = fRecoParam->GetMinMaxCutSigma();
1096   Float_t minLeftRightCutSigma = fRecoParam->GetMinLeftRightCutSigma();
1097   Float_t minUpDownCutSigma    = fRecoParam->GetMinUpDownCutSigma();
1098   for (Int_t iSig = 0; iSig < fNSigBins; iSig++) {
1099     Int_t i = fSigBins[iSig];
1100     if (i%fMaxTime<=crtime) continue;
1101     Float_t *b = &fBins[i];
1102     //absolute custs
1103     if (b[0]<minMaxCutAbs) continue;   //threshold for maxima  
1104     //
1105     if (b[-1]+b[1]+b[-fMaxTime]+b[fMaxTime]<=0) continue;  // cut on isolated clusters 
1106     if (b[-1]+b[1]<=0) continue;               // cut on isolated clusters
1107     if (b[-fMaxTime]+b[fMaxTime]<=0) continue; // cut on isolated clusters
1108     //
1109     if ((b[0]+b[-1]+b[1])<minUpDownCutAbs) continue;   //threshold for up down  (TRF) 
1110     if ((b[0]+b[-fMaxTime]+b[fMaxTime])<minLeftRightCutAbs) continue;   //threshold for left right (PRF)    
1111     if (!IsMaximum(*b,fMaxTime,b)) continue;
1112     //
1113     Float_t noise = noiseROC->GetValue(fRow, i/fMaxTime);
1114     if (noise>fRecoParam->GetMaxNoise()) continue;
1115     // sigma cuts
1116     if (b[0]<minMaxCutSigma*noise) continue;   //threshold form maxima  
1117     if ((b[0]+b[-1]+b[1])<minUpDownCutSigma*noise) continue;   //threshold for up town TRF 
1118     if ((b[0]+b[-fMaxTime]+b[fMaxTime])<minLeftRightCutSigma*noise) continue;   //threshold for left right (PRF)    
1119   
1120     AliTPCclusterMI c;   // default cosntruction  without info
1121     Int_t dummy=0;
1122     MakeCluster(i, fMaxTime, fBins, dummy,c);
1123     
1124     //}
1125   }
1126 }
1127
1128 Bool_t AliTPCclustererMI::AcceptCluster(AliTPCclusterMI *cl){
1129   //
1130   // Currently hack to filter digital noise (15.06.2008)
1131   // To be parameterized in the AliTPCrecoParam
1132   // More inteligent way  to be used in future
1133   // Acces to the proper pedestal file needed
1134   //
1135   if (cl->GetMax()<400) return kTRUE;
1136   Double_t ratio = cl->GetQ()/cl->GetMax();
1137   if (cl->GetMax()>700){
1138     if ((ratio - int(ratio)>0.8)) return kFALSE;
1139   }
1140   if ((ratio - int(ratio)<0.95)) return kTRUE;
1141   return kFALSE;
1142 }
1143
1144
1145 Double_t AliTPCclustererMI::ProcesSignal(Float_t *signal, Int_t nchannels, Int_t id[3], Double_t &rmsEvent, Double_t &pedestalEvent){
1146   //
1147   // process signal on given pad - + streaming of additional information in special mode
1148   //
1149   // id[0] - sector
1150   // id[1] - row
1151   // id[2] - pad 
1152
1153   //
1154   // ESTIMATE pedestal and the noise
1155   // 
1156   const Int_t kPedMax = 100;
1157   Float_t  max    =  0;
1158   Float_t  maxPos =  0;
1159   Int_t    median =  -1;
1160   Int_t    count0 =  0;
1161   Int_t    count1 =  0;
1162   Float_t  rmsCalib   = rmsEvent;       // backup initial value ( from calib)
1163   Float_t  pedestalCalib = pedestalEvent;// backup initial value ( from calib)
1164   Int_t    firstBin = fRecoParam->GetFirstBin();
1165   //
1166   UShort_t histo[kPedMax];
1167   //memset(histo,0,kPedMax*sizeof(UShort_t));
1168   for (Int_t i=0; i<kPedMax; i++) histo[i]=0;
1169   for (Int_t i=0; i<fMaxTime; i++){
1170     if (signal[i]<=0) continue;
1171     if (signal[i]>max && i>firstBin) {
1172       max = signal[i];
1173       maxPos = i;
1174     }
1175     if (signal[i]>kPedMax-1) continue;
1176     histo[int(signal[i]+0.5)]++;
1177     count0++;
1178   }
1179   //
1180   for (Int_t i=1; i<kPedMax; i++){
1181     if (count1<count0*0.5) median=i;
1182     count1+=histo[i];
1183   }
1184   // truncated mean  
1185   //
1186   Float_t count10=histo[median] ,mean=histo[median]*median,  rms=histo[median]*median*median ;
1187   Float_t count06=histo[median] ,mean06=histo[median]*median,  rms06=histo[median]*median*median ;
1188   Float_t count09=histo[median] ,mean09=histo[median]*median,  rms09=histo[median]*median*median ;
1189   //
1190   for (Int_t idelta=1; idelta<10; idelta++){
1191     if (median-idelta<=0) continue;
1192     if (median+idelta>kPedMax) continue;
1193     if (count06<0.6*count1){
1194       count06+=histo[median-idelta];
1195       mean06 +=histo[median-idelta]*(median-idelta);
1196       rms06  +=histo[median-idelta]*(median-idelta)*(median-idelta);
1197       count06+=histo[median+idelta];
1198       mean06 +=histo[median+idelta]*(median+idelta);
1199       rms06  +=histo[median+idelta]*(median+idelta)*(median+idelta);
1200     }
1201     if (count09<0.9*count1){
1202       count09+=histo[median-idelta];
1203       mean09 +=histo[median-idelta]*(median-idelta);
1204       rms09  +=histo[median-idelta]*(median-idelta)*(median-idelta);
1205       count09+=histo[median+idelta];
1206       mean09 +=histo[median+idelta]*(median+idelta);
1207       rms09  +=histo[median+idelta]*(median+idelta)*(median+idelta);
1208     }
1209     if (count10<0.95*count1){
1210       count10+=histo[median-idelta];
1211       mean +=histo[median-idelta]*(median-idelta);
1212       rms  +=histo[median-idelta]*(median-idelta)*(median-idelta);
1213       count10+=histo[median+idelta];
1214       mean +=histo[median+idelta]*(median+idelta);
1215       rms  +=histo[median+idelta]*(median+idelta)*(median+idelta);
1216     }
1217   }
1218   if (count10) {
1219     mean  /=count10;
1220     rms    = TMath::Sqrt(TMath::Abs(rms/count10-mean*mean));
1221   }
1222   if (count06) {
1223     mean06/=count06;
1224     rms06    = TMath::Sqrt(TMath::Abs(rms06/count06-mean06*mean06));
1225   }
1226   if (count09) {
1227     mean09/=count09;
1228     rms09    = TMath::Sqrt(TMath::Abs(rms09/count09-mean09*mean09));
1229   }
1230   rmsEvent = rms09;
1231   //
1232   pedestalEvent = median;
1233   if (AliLog::GetDebugLevel("","AliTPCclustererMI")==0) return median;
1234   //
1235   UInt_t uid[3] = {UInt_t(id[0]),UInt_t(id[1]),UInt_t(id[2])};
1236   //
1237   // Dump mean signal info
1238   //
1239     if (AliTPCReconstructor::StreamLevel()>0) {
1240     (*fDebugStreamer)<<"Signal"<<
1241     "TimeStamp="<<fTimeStamp<<
1242     "EventType="<<fEventType<<
1243     "Sector="<<uid[0]<<
1244     "Row="<<uid[1]<<
1245     "Pad="<<uid[2]<<
1246     "Max="<<max<<
1247     "MaxPos="<<maxPos<<
1248     //
1249     "Median="<<median<<
1250     "Mean="<<mean<<
1251     "RMS="<<rms<<      
1252     "Mean06="<<mean06<<
1253     "RMS06="<<rms06<<
1254     "Mean09="<<mean09<<
1255     "RMS09="<<rms09<<
1256     "RMSCalib="<<rmsCalib<<
1257     "PedCalib="<<pedestalCalib<<
1258     "\n";
1259     }
1260   //
1261   // fill pedestal histogram
1262   //
1263   //
1264   //
1265   //
1266   Float_t kMin =fRecoParam->GetDumpAmplitudeMin();   // minimal signal to be dumped
1267   Float_t *dsignal = new Float_t[nchannels];
1268   Float_t *dtime   = new Float_t[nchannels];
1269   for (Int_t i=0; i<nchannels; i++){
1270     dtime[i] = i;
1271     dsignal[i] = signal[i];
1272   }
1273
1274   TGraph * graph=0;
1275   //
1276   // Big signals dumping
1277   //    
1278   if (AliTPCReconstructor::StreamLevel()>0) {
1279   if (max-median>kMin &&maxPos>fRecoParam->GetFirstBin()) 
1280     (*fDebugStreamer)<<"SignalB"<<     // pads with signal
1281       "TimeStamp="<<fTimeStamp<<
1282       "EventType="<<fEventType<<
1283       "Sector="<<uid[0]<<
1284       "Row="<<uid[1]<<
1285       "Pad="<<uid[2]<<
1286       "Graph="<<graph<<
1287       "Max="<<max<<
1288       "MaxPos="<<maxPos<<       
1289       //
1290       "Median="<<median<<
1291       "Mean="<<mean<<
1292       "RMS="<<rms<<      
1293       "Mean06="<<mean06<<
1294       "RMS06="<<rms06<<
1295       "Mean09="<<mean09<<
1296       "RMS09="<<rms09<<
1297       "\n";
1298   delete graph;  
1299   }
1300
1301   delete [] dsignal;
1302   delete [] dtime;
1303   if (rms06>fRecoParam->GetMaxNoise()) {
1304     pedestalEvent+=1024.;
1305     return 1024+median; // sign noisy channel in debug mode
1306   }
1307   return median;
1308 }
1309
1310
1311
1312
1313