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