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