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