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