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 #include <TVirtualFFT.h>
  34
  35 #include "AliDigits.h"
  36 #include "AliLoader.h"
  37 #include "AliLog.h"
  38 #include "AliMathBase.h"
  39 #include "AliRawEventHeaderBase.h"
  40 #include "AliRawReader.h"
  41 #include "AliRunLoader.h"
  42 #include "AliSimDigits.h"
  43 #include "AliTPCCalPad.h"
  44 #include "AliTPCCalROC.h"
  45 #include "AliTPCClustersArray.h"
  46 #include "AliTPCClustersRow.h"
  47 #include "AliTPCParam.h"
  48 #include "AliTPCRawStream.h"
  49 #include "AliTPCRecoParam.h"
  50 #include "AliTPCReconstructor.h"
  51 #include "AliTPCcalibDB.h"
  52 #include "AliTPCclusterInfo.h"
  53 #include "AliTPCclusterMI.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   fAmplitudeHisto(0),
  87   fDebugStreamer(0),
  88   fRecoParam(0),
  89   fBDumpSignal(kFALSE),
  90   fFFTr2c(0)
  91 {
  92   //
  93   // COSNTRUCTOR
  94   // param     - tpc parameters for given file
  95   // recoparam - reconstruction parameters
  96   //
  97   fIsOldRCUFormat = kFALSE;
  98   fInput =0;
  99   fOutput=0;
 100   fParam = par;
 101   if (recoParam) {
 102     fRecoParam = recoParam;
 103   }else{
 104     //set default parameters if not specified
 105     fRecoParam = AliTPCReconstructor::GetRecoParam();
 106     if (!fRecoParam)  fRecoParam = AliTPCRecoParam::GetLowFluxParam();
 107   }
 108   fDebugStreamer = new TTreeSRedirector("TPCsignal.root");
 109   fAmplitudeHisto = 0;
 110   Int_t nPoints = fRecoParam->GetLastBin()-fRecoParam->GetFirstBin();
 111   fFFTr2c = TVirtualFFT::FFT(1, &nPoints, "R2C  K");
 112 }
 113 //______________________________________________________________
 114 AliTPCclustererMI::AliTPCclustererMI(const AliTPCclustererMI &param)
 115               :TObject(param),
 116   fBins(0),
 117   fSigBins(0),
 118   fNSigBins(0),
 119   fLoop(0),
 120   fMaxBin(0),
 121   fMaxTime(0),
 122   fMaxPad(0),
 123   fSector(-1),
 124   fRow(-1),
 125   fSign(0),
 126   fRx(0),
 127   fPadWidth(0),
 128   fPadLength(0),
 129   fZWidth(0),
 130   fPedSubtraction(kFALSE),
 131   fIsOldRCUFormat(kFALSE),
 132   fEventHeader(0),
 133   fTimeStamp(0),
 134   fEventType(0),
 135   fInput(0),
 136   fOutput(0),
 137   fRowCl(0),
 138   fRowDig(0),
 139   fParam(0),
 140   fNcluster(0),
 141   fAmplitudeHisto(0),
 142   fDebugStreamer(0),
 143   fRecoParam(0),
 144   fBDumpSignal(kFALSE),
 145   fFFTr2c(0)
 146 {
 147   //
 148   // dummy
 149   //
 150   fMaxBin = param.fMaxBin;
 151 }
 152 //______________________________________________________________
 153 AliTPCclustererMI & AliTPCclustererMI::operator =(const AliTPCclustererMI & param)
 154 {
 155   //
 156   // assignment operator - dummy
 157   //
 158   fMaxBin=param.fMaxBin;
 159   return (*this);
 160 }
 161 //______________________________________________________________
 162 AliTPCclustererMI::~AliTPCclustererMI(){
 163   DumpHistos();
 164   if (fAmplitudeHisto) delete fAmplitudeHisto;
 165   if (fDebugStreamer) delete fDebugStreamer;
 166 }
 167
 168 void AliTPCclustererMI::SetInput(TTree * tree)
 169 {
 170   //
 171   // set input tree with digits
 172   //
 173   fInput = tree;
 174   if  (!fInput->GetBranch("Segment")){
 175     cerr<<"AliTPC::Digits2Clusters(): no porper input tree !\n";
 176     fInput=0;
 177     return;
 178   }
 179 }
 180
 181 void AliTPCclustererMI::SetOutput(TTree * tree)
 182 {
 183   //
 184   //
 185   fOutput= tree;
 186   AliTPCClustersRow clrow;
 187   AliTPCClustersRow *pclrow=&clrow;
 188   clrow.SetClass("AliTPCclusterMI");
 189   clrow.SetArray(1); // to make Clones array
 190   fOutput->Branch("Segment","AliTPCClustersRow",&pclrow,32000,200);
 191 }
 192
 193
 194 Float_t  AliTPCclustererMI::GetSigmaY2(Int_t iz){
 195   // sigma y2 = in digits  - we don't know the angle
 196   Float_t z = iz*fParam->GetZWidth()+fParam->GetNTBinsL1()*fParam->GetZWidth();
 197   Float_t sd2 = (z*fParam->GetDiffL()*fParam->GetDiffL())/
 198     (fPadWidth*fPadWidth);
 199   Float_t sres = 0.25;
 200   Float_t res = sd2+sres;
 201   return res;
 202 }
 203
 204
 205 Float_t  AliTPCclustererMI::GetSigmaZ2(Int_t iz){
 206   //sigma z2 = in digits - angle estimated supposing vertex constraint
 207   Float_t z = iz*fZWidth+fParam->GetNTBinsL1()*fParam->GetZWidth();
 208   Float_t sd2 = (z*fParam->GetDiffL()*fParam->GetDiffL())/(fZWidth*fZWidth);
 209   Float_t angular = fPadLength*(fParam->GetZLength(fSector)-z)/(fRx*fZWidth);
 210   angular*=angular;
 211   angular/=12.;
 212   Float_t sres = fParam->GetZSigma()/fZWidth;
 213   sres *=sres;
 214   Float_t res = angular +sd2+sres;
 215   return res;
 216 }
 217
 218 void AliTPCclustererMI::MakeCluster(Int_t k,Int_t max,Float_t *bins, UInt_t /*m*/,
 219 AliTPCclusterMI &c)
 220 {
 221   //
 222   //  k    - Make cluster at position k
 223   //  bins - 2 D array of signals mapped to 1 dimensional array -
 224   //  max  - the number of time bins er one dimension
 225   //  c    - refernce to cluster to be filled
 226   //
 227   Int_t i0=k/max;  //central pad
 228   Int_t j0=k%max;  //central time bin
 229
 230   // set pointers to data
 231   //Int_t dummy[5] ={0,0,0,0,0};
 232   Float_t * matrix[5]; //5x5 matrix with digits  - indexing i = 0 ..4  j = -2..2
 233   for (Int_t di=-2;di<=2;di++){
 234     matrix[di+2]  =  &bins[k+di*max];
 235   }
 236   //build matrix with virtual charge
 237   Float_t sigmay2= GetSigmaY2(j0);
 238   Float_t sigmaz2= GetSigmaZ2(j0);
 239
 240   Float_t vmatrix[5][5];
 241   vmatrix[2][2] = matrix[2][0];
 242   c.SetType(0);
 243   c.SetMax((UShort_t)(vmatrix[2][2])); // write maximal amplitude
 244   for (Int_t di =-1;di <=1;di++)
 245     for (Int_t dj =-1;dj <=1;dj++){
 246       Float_t amp = matrix[di+2][dj];
 247       if ( (amp<2) && (fLoop<2)){
 248         // if under threshold  - calculate virtual charge
 249         Float_t ratio = TMath::Exp(-1.2*TMath::Abs(di)/sigmay2)*TMath::Exp(-1.2*TMath::Abs(dj)/sigmaz2);
 250         amp = ((matrix[2][0]-2)*(matrix[2][0]-2)/(matrix[-di+2][-dj]+2))*ratio;
 251         if (amp>2) amp = 2;
 252         vmatrix[2+di][2+dj]=amp;
 253         vmatrix[2+2*di][2+2*dj]=0;
 254         if ( (di*dj)!=0){
 255           //DIAGONAL ELEMENTS
 256           vmatrix[2+2*di][2+dj] =0;
 257           vmatrix[2+di][2+2*dj] =0;
 258         }
 259         continue;
 260       }
 261       if (amp<4){
 262         //if small  amplitude - below  2 x threshold  - don't consider other one
 263         vmatrix[2+di][2+dj]=amp;
 264         vmatrix[2+2*di][2+2*dj]=0;  // don't take to the account next bin
 265         if ( (di*dj)!=0){
 266           //DIAGONAL ELEMENTS
 267           vmatrix[2+2*di][2+dj] =0;
 268           vmatrix[2+di][2+2*dj] =0;
 269         }
 270         continue;
 271       }
 272       //if bigger then take everything
 273       vmatrix[2+di][2+dj]=amp;
 274       vmatrix[2+2*di][2+2*dj]= matrix[2*di+2][2*dj] ;
 275       if ( (di*dj)!=0){
 276           //DIAGONAL ELEMENTS
 277           vmatrix[2+2*di][2+dj] = matrix[2*di+2][dj];
 278           vmatrix[2+di][2+2*dj] = matrix[2+di][dj*2];
 279         }
 280     }
 281
 282
 283
 284   Float_t sumw=0;
 285   Float_t sumiw=0;
 286   Float_t sumi2w=0;
 287   Float_t sumjw=0;
 288   Float_t sumj2w=0;
 289   //
 290   for (Int_t i=-2;i<=2;i++)
 291     for (Int_t j=-2;j<=2;j++){
 292       Float_t amp = vmatrix[i+2][j+2];
 293
 294       sumw    += amp;
 295       sumiw   += i*amp;
 296       sumi2w  += i*i*amp;
 297       sumjw   += j*amp;
 298       sumj2w  += j*j*amp;
 299     }
 300   //
 301   Float_t meani = sumiw/sumw;
 302   Float_t mi2   = sumi2w/sumw-meani*meani;
 303   Float_t meanj = sumjw/sumw;
 304   Float_t mj2   = sumj2w/sumw-meanj*meanj;
 305   //
 306   Float_t ry = mi2/sigmay2;
 307   Float_t rz = mj2/sigmaz2;
 308
 309   //
 310   if ( ( (ry<0.6) || (rz<0.6) ) && fLoop==2) return;
 311   if ( (ry <1.2) && (rz<1.2) || (!fRecoParam->GetDoUnfold())) {
 312     //
 313     //if cluster looks like expected or Unfolding not switched on
 314     //standard COG is used
 315     //+1.2 deviation from expected sigma accepted
 316     //    c.fMax = FitMax(vmatrix,meani,meanj,TMath::Sqrt(sigmay2),TMath::Sqrt(sigmaz2));
 317
 318     meani +=i0;
 319     meanj +=j0;
 320     //set cluster parameters
 321     c.SetQ(sumw);
 322     c.SetY(meani*fPadWidth);
 323     c.SetZ(meanj*fZWidth);
 324     c.SetPad(meani);
 325     c.SetTimeBin(meanj);
 326     c.SetSigmaY2(mi2);
 327     c.SetSigmaZ2(mj2);
 328     AddCluster(c,(Float_t*)vmatrix,k);
 329
 330     return;
 331   }
 332   //
 333   //unfolding when neccessary
 334   //
 335
 336   Float_t * matrix2[7]; //7x7 matrix with digits  - indexing i = 0 ..6  j = -3..3
 337   Float_t dummy[7]={0,0,0,0,0,0};
 338   for (Int_t di=-3;di<=3;di++){
 339     matrix2[di+3] =  &bins[k+di*max];
 340     if ((k+di*max)<3)  matrix2[di+3] = &dummy[3];
 341     if ((k+di*max)>fMaxBin-3)  matrix2[di+3] = &dummy[3];
 342   }
 343   Float_t vmatrix2[5][5];
 344   Float_t sumu;
 345   Float_t overlap;
 346   UnfoldCluster(matrix2,vmatrix2,meani,meanj,sumu,overlap);
 347   //
 348   //  c.fMax = FitMax(vmatrix2,meani,meanj,TMath::Sqrt(sigmay2),TMath::Sqrt(sigmaz2));
 349   meani +=i0;
 350   meanj +=j0;
 351   //set cluster parameters
 352   c.SetQ(sumu);
 353   c.SetY(meani*fPadWidth);
 354   c.SetZ(meanj*fZWidth);
 355   c.SetPad(meani);
 356   c.SetTimeBin(meanj);
 357   c.SetSigmaY2(mi2);
 358   c.SetSigmaZ2(mj2);
 359   c.SetType(Char_t(overlap)+1);
 360   AddCluster(c,(Float_t*)vmatrix,k);
 361
 362   //unfolding 2
 363   meani-=i0;
 364   meanj-=j0;
 365   if (gDebug>4)
 366     printf("%f\t%f\n", vmatrix2[2][2], vmatrix[2][2]);
 367 }
 368
 369
 370
 371 void AliTPCclustererMI::UnfoldCluster(Float_t * matrix2[7], Float_t recmatrix[5][5], Float_t & meani, Float_t & meanj,
 372                                       Float_t & sumu, Float_t & overlap )
 373 {
 374   //
 375   //unfold cluster from input matrix
 376   //data corresponding to cluster writen in recmatrix
 377   //output meani and meanj
 378
 379   //take separatelly y and z
 380
 381   Float_t sum3i[7] = {0,0,0,0,0,0,0};
 382   Float_t sum3j[7] = {0,0,0,0,0,0,0};
 383
 384   for (Int_t k =0;k<7;k++)
 385     for (Int_t l = -1; l<=1;l++){
 386       sum3i[k]+=matrix2[k][l];
 387       sum3j[k]+=matrix2[l+3][k-3];
 388     }
 389   Float_t mratio[3][3]={{1,1,1},{1,1,1},{1,1,1}};
 390   //
 391   //unfold  y
 392   Float_t sum3wi    = 0;  //charge minus overlap
 393   Float_t sum3wio   = 0;  //full charge
 394   Float_t sum3iw    = 0;  //sum for mean value
 395   for (Int_t dk=-1;dk<=1;dk++){
 396     sum3wio+=sum3i[dk+3];
 397     if (dk==0){
 398       sum3wi+=sum3i[dk+3];
 399     }
 400     else{
 401       Float_t ratio =1;
 402       if (  ( ((sum3i[dk+3]+3)/(sum3i[3]-3))+1 < (sum3i[2*dk+3]-3)/(sum3i[dk+3]+3))||
 403            sum3i[dk+3]<=sum3i[2*dk+3] && sum3i[dk+3]>2 ){
 404         Float_t xm2 = sum3i[-dk+3];
 405         Float_t xm1 = sum3i[+3];
 406         Float_t x1  = sum3i[2*dk+3];
 407         Float_t x2  = sum3i[3*dk+3];
 408         Float_t w11   = TMath::Max((Float_t)(4.*xm1-xm2),(Float_t)0.000001);
 409         Float_t w12   = TMath::Max((Float_t)(4 *x1 -x2),(Float_t)0.);
 410         ratio = w11/(w11+w12);
 411         for (Int_t dl=-1;dl<=1;dl++)
 412           mratio[dk+1][dl+1] *= ratio;
 413       }
 414       Float_t amp = sum3i[dk+3]*ratio;
 415       sum3wi+=amp;
 416       sum3iw+= dk*amp;
 417     }
 418   }
 419   meani = sum3iw/sum3wi;
 420   Float_t overlapi = (sum3wio-sum3wi)/sum3wio;
 421
 422
 423
 424   //unfold  z
 425   Float_t sum3wj    = 0;  //charge minus overlap
 426   Float_t sum3wjo   = 0;  //full charge
 427   Float_t sum3jw    = 0;  //sum for mean value
 428   for (Int_t dk=-1;dk<=1;dk++){
 429     sum3wjo+=sum3j[dk+3];
 430     if (dk==0){
 431       sum3wj+=sum3j[dk+3];
 432     }
 433     else{
 434       Float_t ratio =1;
 435       if ( ( ((sum3j[dk+3]+3)/(sum3j[3]-3))+1 < (sum3j[2*dk+3]-3)/(sum3j[dk+3]+3)) ||
 436            (sum3j[dk+3]<=sum3j[2*dk+3] && sum3j[dk+3]>2)){
 437         Float_t xm2 = sum3j[-dk+3];
 438         Float_t xm1 = sum3j[+3];
 439         Float_t x1  = sum3j[2*dk+3];
 440         Float_t x2  = sum3j[3*dk+3];
 441         Float_t w11   = TMath::Max((Float_t)(4.*xm1-xm2),(Float_t)0.000001);
 442         Float_t w12   = TMath::Max((Float_t)(4 *x1 -x2),(Float_t)0.);
 443         ratio = w11/(w11+w12);
 444         for (Int_t dl=-1;dl<=1;dl++)
 445           mratio[dl+1][dk+1] *= ratio;
 446       }
 447       Float_t amp = sum3j[dk+3]*ratio;
 448       sum3wj+=amp;
 449       sum3jw+= dk*amp;
 450     }
 451   }
 452   meanj = sum3jw/sum3wj;
 453   Float_t overlapj = (sum3wjo-sum3wj)/sum3wjo;
 454   overlap = Int_t(100*TMath::Max(overlapi,overlapj)+3);
 455   sumu = (sum3wj+sum3wi)/2.;
 456
 457   if (overlap ==3) {
 458     //if not overlap detected remove everything
 459     for (Int_t di =-2; di<=2;di++)
 460       for (Int_t dj =-2; dj<=2;dj++){
 461         recmatrix[di+2][dj+2] = matrix2[3+di][dj];
 462       }
 463   }
 464   else{
 465     for (Int_t di =-1; di<=1;di++)
 466       for (Int_t dj =-1; dj<=1;dj++){
 467         Float_t ratio =1;
 468         if (mratio[di+1][dj+1]==1){
 469           recmatrix[di+2][dj+2]     = matrix2[3+di][dj];
 470           if (TMath::Abs(di)+TMath::Abs(dj)>1){
 471             recmatrix[2*di+2][dj+2] = matrix2[3+2*di][dj];
 472             recmatrix[di+2][2*dj+2] = matrix2[3+di][2*dj];
 473           }
 474           recmatrix[2*di+2][2*dj+2] = matrix2[3+2*di][2*dj];
 475         }
 476         else
 477           {
 478             //if we have overlap in direction
 479             recmatrix[di+2][dj+2] = mratio[di+1][dj+1]* matrix2[3+di][dj];
 480             if (TMath::Abs(di)+TMath::Abs(dj)>1){
 481               ratio =  TMath::Min((Float_t)(recmatrix[di+2][dj+2]/(matrix2[3+0*di][1*dj]+1)),(Float_t)1.);
 482               recmatrix[2*di+2][dj+2] = ratio*recmatrix[di+2][dj+2];
 483               //
 484               ratio =  TMath::Min((Float_t)(recmatrix[di+2][dj+2]/(matrix2[3+1*di][0*dj]+1)),(Float_t)1.);
 485               recmatrix[di+2][2*dj+2] = ratio*recmatrix[di+2][dj+2];
 486             }
 487             else{
 488               ratio =  recmatrix[di+2][dj+2]/matrix2[3][0];
 489               recmatrix[2*di+2][2*dj+2] = ratio*recmatrix[di+2][dj+2];
 490             }
 491           }
 492       }
 493   }
 494   if (gDebug>4)
 495     printf("%f\n", recmatrix[2][2]);
 496
 497 }
 498
 499 Float_t AliTPCclustererMI::FitMax(Float_t vmatrix[5][5], Float_t y, Float_t z, Float_t sigmay, Float_t sigmaz)
 500 {
 501   //
 502   // estimate max
 503   Float_t sumteor= 0;
 504   Float_t sumamp = 0;
 505
 506   for (Int_t di = -1;di<=1;di++)
 507     for (Int_t dj = -1;dj<=1;dj++){
 508       if (vmatrix[2+di][2+dj]>2){
 509         Float_t teor = TMath::Gaus(di,y,sigmay*1.2)*TMath::Gaus(dj,z,sigmaz*1.2);
 510         sumteor += teor*vmatrix[2+di][2+dj];
 511         sumamp  += vmatrix[2+di][2+dj]*vmatrix[2+di][2+dj];
 512       }
 513     }
 514   Float_t max = sumamp/sumteor;
 515   return max;
 516 }
 517
 518 void AliTPCclustererMI::AddCluster(AliTPCclusterMI &c, Float_t * matrix, Int_t pos){
 519   //
 520   // transform cluster to the global coordinata
 521   // add the cluster to the array
 522   //
 523   Float_t meani = c.GetY()/fPadWidth;
 524   Float_t meanj = c.GetZ()/fZWidth;
 525
 526   Int_t ki = TMath::Nint(meani-3);
 527   if (ki<0) ki=0;
 528   if (ki>=fMaxPad) ki = fMaxPad-1;
 529   Int_t kj = TMath::Nint(meanj-3);
 530   if (kj<0) kj=0;
 531   if (kj>=fMaxTime-3) kj=fMaxTime-4;
 532   // ki and kj shifted to "real" coordinata
 533   if (fRowDig) {
 534     c.SetLabel(fRowDig->GetTrackIDFast(kj,ki,0)-2,0);
 535     c.SetLabel(fRowDig->GetTrackIDFast(kj,ki,1)-2,1);
 536     c.SetLabel(fRowDig->GetTrackIDFast(kj,ki,2)-2,2);
 537   }
 538
 539
 540   Float_t s2 = c.GetSigmaY2();
 541   Float_t w=fParam->GetPadPitchWidth(fSector);
 542
 543   c.SetSigmaY2(s2*w*w);
 544   s2 = c.GetSigmaZ2();
 545   w=fZWidth;
 546   c.SetSigmaZ2(s2*w*w);
 547   c.SetY((meani - 2.5 - 0.5*fMaxPad)*fParam->GetPadPitchWidth(fSector));
 548   if (!fRecoParam->GetBYMirror()){
 549     if (fSector%36>17){
 550       c.SetY(-(meani - 2.5 - 0.5*fMaxPad)*fParam->GetPadPitchWidth(fSector));
 551     }
 552   }
 553   c.SetZ(fZWidth*(meanj-3));
 554   c.SetZ(c.GetZ() - 3.*fParam->GetZSigma() + fParam->GetNTBinsL1()*fParam->GetZWidth()); // PASA delay + L1 delay
 555   c.SetZ(fSign*(fParam->GetZLength(fSector) - c.GetZ()));
 556   c.SetX(fRx);
 557   c.SetDetector(fSector);
 558   c.SetRow(fRow);
 559
 560   if (ki<=1 || ki>=fMaxPad-1 || kj==1 || kj==fMaxTime-2) {
 561     //c.SetSigmaY2(c.GetSigmaY2()*25.);
 562     //c.SetSigmaZ2(c.GetSigmaZ2()*4.);
 563     c.SetType(-(c.GetType()+3));  //edge clusters
 564   }
 565   if (fLoop==2) c.SetType(100);
 566
 567   TClonesArray * arr = fRowCl->GetArray();
 568   AliTPCclusterMI * cl = new ((*arr)[fNcluster]) AliTPCclusterMI(c);
 569   if (matrix ) {
 570     Int_t nbins=0;
 571     Float_t *graph =0;
 572     if (fRecoParam->GetCalcPedestal() && cl->GetMax()>fRecoParam->GetDumpAmplitudeMin() &&fBDumpSignal){
 573       nbins = fMaxTime;
 574       graph = &(fBins[fMaxTime*(pos/fMaxTime)]);
 575     }
 576     AliTPCclusterInfo * info = new AliTPCclusterInfo(matrix,nbins,graph);
 577     cl->SetInfo(info);
 578   }
 579
 580   fNcluster++;
 581 }
 582
 583
 584 //_____________________________________________________________________________
 585 void AliTPCclustererMI::Digits2Clusters()
 586 {
 587   //-----------------------------------------------------------------
 588   // This is a simple cluster finder.
 589   //-----------------------------------------------------------------
 590
 591   if (!fInput) {
 592     Error("Digits2Clusters", "input tree not initialised");
 593     return;
 594   }
 595
 596   if (!fOutput) {
 597     Error("Digits2Clusters", "output tree not initialised");
 598     return;
 599   }
 600
 601   AliTPCCalPad * gainTPC = AliTPCcalibDB::Instance()->GetPadGainFactor();
 602   AliTPCCalPad * noiseTPC = AliTPCcalibDB::Instance()->GetPadNoise();
 603
 604   AliSimDigits digarr, *dummy=&digarr;
 605   fRowDig = dummy;
 606   fInput->GetBranch("Segment")->SetAddress(&dummy);
 607   Stat_t nentries = fInput->GetEntries();
 608
 609   fMaxTime=fParam->GetMaxTBin()+6; // add 3 virtual time bins before and 3   after
 610
 611   Int_t nclusters  = 0;
 612
 613   for (Int_t n=0; n<nentries; n++) {
 614     fInput->GetEvent(n);
 615     if (!fParam->AdjustSectorRow(digarr.GetID(),fSector,fRow)) {
 616       cerr<<"AliTPC warning: invalid segment ID ! "<<digarr.GetID()<<endl;
 617       continue;
 618     }
 619     Int_t row = fRow;
 620     AliTPCCalROC * gainROC = gainTPC->GetCalROC(fSector);  // pad gains per given sector
 621     AliTPCCalROC * noiseROC   = noiseTPC->GetCalROC(fSector); // noise per given sector
 622     //
 623     AliTPCClustersRow *clrow= new AliTPCClustersRow();
 624     fRowCl = clrow;
 625     clrow->SetClass("AliTPCclusterMI");
 626     clrow->SetArray(1);
 627
 628     clrow->SetID(digarr.GetID());
 629     fOutput->GetBranch("Segment")->SetAddress(&clrow);
 630     fRx=fParam->GetPadRowRadii(fSector,row);
 631
 632
 633     const Int_t kNIS=fParam->GetNInnerSector(), kNOS=fParam->GetNOuterSector();
 634     fZWidth = fParam->GetZWidth();
 635     if (fSector < kNIS) {
 636       fMaxPad = fParam->GetNPadsLow(row);
 637       fSign = (fSector < kNIS/2) ? 1 : -1;
 638       fPadLength = fParam->GetPadPitchLength(fSector,row);
 639       fPadWidth = fParam->GetPadPitchWidth();
 640     } else {
 641       fMaxPad = fParam->GetNPadsUp(row);
 642       fSign = ((fSector-kNIS) < kNOS/2) ? 1 : -1;
 643       fPadLength = fParam->GetPadPitchLength(fSector,row);
 644       fPadWidth  = fParam->GetPadPitchWidth();
 645     }
 646
 647
 648     fMaxBin=fMaxTime*(fMaxPad+6);  // add 3 virtual pads  before and 3 after
 649     fBins    =new Float_t[fMaxBin];
 650     fSigBins =new Int_t[fMaxBin];
 651     fNSigBins = 0;
 652     memset(fBins,0,sizeof(Float_t)*fMaxBin);
 653
 654     if (digarr.First()) //MI change
 655       do {
 656         Float_t dig=digarr.CurrentDigit();
 657         if (dig<=fParam->GetZeroSup()) continue;
 658         Int_t j=digarr.CurrentRow()+3, i=digarr.CurrentColumn()+3;
 659         Float_t gain = gainROC->GetValue(row,digarr.CurrentColumn());
 660         Int_t bin = i*fMaxTime+j;
 661         fBins[bin]=dig/gain;
 662         fSigBins[fNSigBins++]=bin;
 663       } while (digarr.Next());
 664     digarr.ExpandTrackBuffer();
 665
 666     FindClusters(noiseROC);
 667
 668     fOutput->Fill();
 669     delete clrow;
 670     nclusters+=fNcluster;
 671     delete[] fBins;
 672     delete[] fSigBins;
 673   }
 674
 675   Info("Digits2Clusters", "Number of found clusters : %d", nclusters);
 676 }
 677
 678 void AliTPCclustererMI::Digits2Clusters(AliRawReader* rawReader)
 679 {
 680 //-----------------------------------------------------------------
 681 // This is a cluster finder for the TPC raw data.
 682 // The method assumes NO ordering of the altro channels.
 683 // The pedestal subtraction can be switched on and off
 684 // using an option of the TPC reconstructor
 685 //-----------------------------------------------------------------
 686
 687   if (!fOutput) {
 688     Error("Digits2Clusters", "output tree not initialised");
 689     return;
 690   }
 691
 692   fRowDig = NULL;
 693   AliTPCROC * roc = AliTPCROC::Instance();
 694   AliTPCCalPad * gainTPC = AliTPCcalibDB::Instance()->GetPadGainFactor();
 695   AliTPCCalPad * pedestalTPC = AliTPCcalibDB::Instance()->GetPedestals();
 696   AliTPCCalPad * noiseTPC = AliTPCcalibDB::Instance()->GetPadNoise();
 697   AliTPCRawStream input(rawReader);
 698   fEventHeader = (AliRawEventHeaderBase*)rawReader->GetEventHeader();
 699   if (fEventHeader){
 700     fTimeStamp = fEventHeader->Get("Timestamp");
 701     fEventType = fEventHeader->Get("Type");
 702   }
 703
 704
 705   Int_t nclusters  = 0;
 706
 707   fMaxTime = fParam->GetMaxTBin() + 6; // add 3 virtual time bins before and 3 after
 708   const Int_t kNIS = fParam->GetNInnerSector();
 709   const Int_t kNOS = fParam->GetNOuterSector();
 710   const Int_t kNS = kNIS + kNOS;
 711   fZWidth = fParam->GetZWidth();
 712   Int_t zeroSup = fParam->GetZeroSup();
 713   //
 714   //alocate memory for sector - maximal case
 715   //
 716   Float_t** allBins = NULL;
 717   Int_t** allSigBins = NULL;
 718   Int_t*  allNSigBins = NULL;
 719   Int_t nRowsMax = roc->GetNRows(roc->GetNSector()-1);
 720   Int_t nPadsMax = roc->GetNPads(roc->GetNSector()-1,nRowsMax-1);
 721   allBins = new Float_t*[nRowsMax];
 722   allSigBins = new Int_t*[nRowsMax];
 723   allNSigBins = new Int_t[nRowsMax];
 724   for (Int_t iRow = 0; iRow < nRowsMax; iRow++) {
 725     //
 726     Int_t maxBin = fMaxTime*(nPadsMax+6);  // add 3 virtual pads  before and 3 after
 727     allBins[iRow] = new Float_t[maxBin];
 728     memset(allBins[iRow],0,sizeof(Float_t)*maxBin);
 729     allSigBins[iRow] = new Int_t[maxBin];
 730     allNSigBins[iRow]=0;
 731   }
 732   //
 733   // Loop over sectors
 734   //
 735   for(fSector = 0; fSector < kNS; fSector++) {
 736
 737     AliTPCCalROC * gainROC    = gainTPC->GetCalROC(fSector);  // pad gains per given sector
 738     AliTPCCalROC * pedestalROC = pedestalTPC->GetCalROC(fSector);  // pedestal per given sector
 739     AliTPCCalROC * noiseROC   = noiseTPC->GetCalROC(fSector);  // noise per given sector
 740
 741     Int_t nRows = 0;
 742     Int_t nDDLs = 0, indexDDL = 0;
 743     if (fSector < kNIS) {
 744       nRows = fParam->GetNRowLow();
 745       fSign = (fSector < kNIS/2) ? 1 : -1;
 746       nDDLs = 2;
 747       indexDDL = fSector * 2;
 748     }
 749     else {
 750       nRows = fParam->GetNRowUp();
 751       fSign = ((fSector-kNIS) < kNOS/2) ? 1 : -1;
 752       nDDLs = 4;
 753       indexDDL = (fSector-kNIS) * 4 + kNIS * 2;
 754     }
 755
 756     for (Int_t iRow = 0; iRow < nRows; iRow++) {
 757       Int_t maxPad;
 758       if (fSector < kNIS)
 759         maxPad = fParam->GetNPadsLow(iRow);
 760       else
 761         maxPad = fParam->GetNPadsUp(iRow);
 762
 763       Int_t maxBin = fMaxTime*(maxPad+6);  // add 3 virtual pads  before and 3 after
 764       memset(allBins[iRow],0,sizeof(Float_t)*maxBin);
 765       allNSigBins[iRow] = 0;
 766     }
 767
 768     // Loas the raw data for corresponding DDLs
 769     rawReader->Reset();
 770     input.SetOldRCUFormat(fIsOldRCUFormat);
 771     rawReader->Select("TPC",indexDDL,indexDDL+nDDLs-1);
 772     Int_t digCounter=0;
 773     // Begin loop over altro data
 774     Bool_t calcPedestal = fRecoParam->GetCalcPedestal();
 775     Float_t gain =1;
 776     Int_t lastPad=-1;
 777     while (input.Next()) {
 778       digCounter++;
 779       if (input.GetSector() != fSector)
 780         AliFatal(Form("Sector index mismatch ! Expected (%d), but got (%d) !",fSector,input.GetSector()));
 781
 782
 783       Int_t iRow = input.GetRow();
 784       if (iRow < 0 || iRow >= nRows)
 785         AliFatal(Form("Pad-row index (%d) outside the range (%d -> %d) !",
 786                       iRow, 0, nRows -1));
 787       //pad
 788       Int_t iPad = input.GetPad();
 789       if (iPad < 0 || iPad >= nPadsMax)
 790         AliFatal(Form("Pad index (%d) outside the range (%d -> %d) !",
 791                       iPad, 0, nPadsMax-1));
 792       if (iPad!=lastPad){
 793         gain    = gainROC->GetValue(iRow,iPad);
 794         lastPad = iPad;
 795       }
 796       iPad+=3;
 797       //time
 798       Int_t iTimeBin = input.GetTime();
 799       if ( iTimeBin < 0 || iTimeBin >= fParam->GetMaxTBin())
 800         AliFatal(Form("Timebin index (%d) outside the range (%d -> %d) !",
 801                       iTimeBin, 0, iTimeBin -1));
 802       iTimeBin+=3;
 803       //signal
 804       Float_t signal = input.GetSignal();
 805       if (!calcPedestal && signal <= zeroSup) continue;
 806       if (!calcPedestal) {
 807         Int_t bin = iPad*fMaxTime+iTimeBin;
 808         allBins[iRow][bin] = signal/gain;
 809         allSigBins[iRow][allNSigBins[iRow]++] = bin;
 810       }else{
 811         allBins[iRow][iPad*fMaxTime+iTimeBin] = signal;
 812       }
 813       allBins[iRow][iPad*fMaxTime+0]=1.;  // pad with signal
 814     } // End of the loop over altro data
 815     //
 816     //
 817     // Now loop over rows and perform pedestal subtraction
 818     if (digCounter==0) continue;
 819     //    if (fPedSubtraction) {
 820     if (calcPedestal) {
 821       for (Int_t iRow = 0; iRow < nRows; iRow++) {
 822         Int_t maxPad;
 823         if (fSector < kNIS)
 824           maxPad = fParam->GetNPadsLow(iRow);
 825         else
 826           maxPad = fParam->GetNPadsUp(iRow);
 827
 828         for (Int_t iPad = 3; iPad < maxPad + 3; iPad++) {
 829           if (allBins[iRow][iPad*fMaxTime+0] <1 ) continue;  // no data
 830           Float_t *p = &allBins[iRow][iPad*fMaxTime+3];
 831           //Float_t pedestal = TMath::Median(fMaxTime, p);
 832           Int_t id[3] = {fSector, iRow, iPad-3};
 833           // calib values
 834           Double_t rmsCalib=  noiseROC->GetValue(iRow,iPad-3);
 835           Double_t pedestalCalib = pedestalROC->GetValue(iRow,iPad-3);
 836           Double_t rmsEvent       = rmsCalib;
 837           Double_t pedestalEvent  = pedestalCalib;
 838           ProcesSignal(p, fMaxTime, id, rmsEvent, pedestalEvent);
 839           if (rmsEvent<rmsCalib) rmsEvent = rmsCalib;   // take worst scenario
 840           if (TMath::Abs(pedestalEvent-pedestalCalib)<1.0) pedestalEvent = pedestalCalib;
 841
 842           //
 843           for (Int_t iTimeBin = 0; iTimeBin < fMaxTime; iTimeBin++) {
 844             Int_t bin = iPad*fMaxTime+iTimeBin;
 845             allBins[iRow][bin] -= pedestalEvent;
 846             if (iTimeBin < AliTPCReconstructor::GetRecoParam()->GetFirstBin())
 847               allBins[iRow][bin] = 0;
 848             if (iTimeBin > AliTPCReconstructor::GetRecoParam()->GetLastBin())
 849               allBins[iRow][bin] = 0;
 850             if (allBins[iRow][iPad*fMaxTime+iTimeBin] < zeroSup)
 851               allBins[iRow][bin] = 0;
 852             if (allBins[iRow][bin] < 3.0*rmsEvent)   // 3 sigma cut on RMS
 853               allBins[iRow][bin] = 0;
 854             if (allBins[iRow][bin]) allSigBins[iRow][allNSigBins[iRow]++] = bin;
 855           }
 856         }
 857       }
 858     }
 859     // Now loop over rows and find clusters
 860     for (fRow = 0; fRow < nRows; fRow++) {
 861       fRowCl = new AliTPCClustersRow;
 862       fRowCl->SetClass("AliTPCclusterMI");
 863       fRowCl->SetArray(1);
 864       fRowCl->SetID(fParam->GetIndex(fSector, fRow));
 865       fOutput->GetBranch("Segment")->SetAddress(&fRowCl);
 866
 867       fRx = fParam->GetPadRowRadii(fSector, fRow);
 868       fPadLength = fParam->GetPadPitchLength(fSector, fRow);
 869       fPadWidth  = fParam->GetPadPitchWidth();
 870       if (fSector < kNIS)
 871         fMaxPad = fParam->GetNPadsLow(fRow);
 872       else
 873         fMaxPad = fParam->GetNPadsUp(fRow);
 874       fMaxBin = fMaxTime*(fMaxPad+6);  // add 3 virtual pads  before and 3 after
 875
 876       fBins = allBins[fRow];
 877       fSigBins = allSigBins[fRow];
 878       fNSigBins = allNSigBins[fRow];
 879
 880       FindClusters(noiseROC);
 881
 882       fOutput->Fill();
 883       delete fRowCl;
 884       nclusters += fNcluster;
 885     } // End of loop to find clusters
 886   } // End of loop over sectors
 887
 888   for (Int_t iRow = 0; iRow < nRowsMax; iRow++) {
 889     delete [] allBins[iRow];
 890     delete [] allSigBins[iRow];
 891   }
 892   delete [] allBins;
 893   delete [] allSigBins;
 894   delete [] allNSigBins;
 895
 896   Info("Digits2Clusters", "File  %s Event\t%d\tNumber of found clusters : %d\n", fOutput->GetName(),*(rawReader->GetEventId()), nclusters);
 897
 898 }
 899
 900 void AliTPCclustererMI::FindClusters(AliTPCCalROC * noiseROC)
 901 {
 902
 903   //
 904   // add virtual charge at the edge
 905   //
 906   Double_t kMaxDumpSize = 500000;
 907   if (fRecoParam->GetCalcPedestal() && fOutput->GetZipBytes()< kMaxDumpSize) fBDumpSignal =kTRUE;   //dump signal flag
 908   //
 909   fNcluster=0;
 910   fLoop=1;
 911   Int_t crtime = Int_t((fParam->GetZLength(fSector)-fRecoParam->GetCtgRange()*fRx)/fZWidth-fParam->GetNTBinsL1()-5);
 912   Float_t minMaxCutAbs       = fRecoParam->GetMinMaxCutAbs();
 913   Float_t minLeftRightCutAbs = fRecoParam->GetMinLeftRightCutAbs();
 914   Float_t minUpDownCutAbs    = fRecoParam->GetMinUpDownCutAbs();
 915   Float_t minMaxCutSigma       = fRecoParam->GetMinMaxCutSigma();
 916   Float_t minLeftRightCutSigma = fRecoParam->GetMinLeftRightCutSigma();
 917   Float_t minUpDownCutSigma    = fRecoParam->GetMinUpDownCutSigma();
 918   for (Int_t iSig = 0; iSig < fNSigBins; iSig++) {
 919     Int_t i = fSigBins[iSig];
 920     if (i%fMaxTime<=crtime) continue;
 921     Float_t *b = &fBins[i];
 922     //absolute custs
 923     if (b[0]<minMaxCutAbs) continue;   //threshold for maxima
 924     //
 925     if (b[-1]+b[1]+b[-fMaxTime]+b[fMaxTime]<=0) continue;  // cut on isolated clusters
 926     //    if (b[-1]+b[1]<=0) continue;               // cut on isolated clusters
 927     //if (b[-fMaxTime]+b[fMaxTime]<=0) continue; // cut on isolated clusters
 928     //
 929     if ((b[0]+b[-1]+b[1])<minUpDownCutAbs) continue;   //threshold for up down  (TRF)
 930     if ((b[0]+b[-fMaxTime]+b[fMaxTime])<minLeftRightCutAbs) continue;   //threshold for left right (PRF)
 931     if (!IsMaximum(*b,fMaxTime,b)) continue;
 932     //
 933     Float_t noise = noiseROC->GetValue(fRow, i/fMaxTime);
 934     // sigma cuts
 935     if (b[0]<minMaxCutSigma*noise) continue;   //threshold form maxima
 936     if ((b[0]+b[-1]+b[1])<minUpDownCutSigma*noise) continue;   //threshold for up town TRF
 937     if ((b[0]+b[-fMaxTime]+b[fMaxTime])<minLeftRightCutSigma*noise) continue;   //threshold for left right (PRF)
 938
 939     AliTPCclusterMI c(kFALSE);   // default cosntruction  without info
 940     Int_t dummy=0;
 941     MakeCluster(i, fMaxTime, fBins, dummy,c);
 942
 943     //}
 944   }
 945 }
 946
 947
 948 Double_t AliTPCclustererMI::ProcesSignal(Float_t *signal, Int_t nchannels, Int_t id[3], Double_t &rmsEvent, Double_t &pedestalEvent){
 949   //
 950   // process signal on given pad - + streaming of additional information in special mode
 951   //
 952   // id[0] - sector
 953   // id[1] - row
 954   // id[2] - pad
 955
 956   //
 957   // ESTIMATE pedestal and the noise
 958   //
 959   const Int_t kPedMax = 100;
 960   Double_t kMaxDebugSize = 5000000.;
 961   Float_t  max    =  0;
 962   Float_t  maxPos =  0;
 963   Int_t    median =  -1;
 964   Int_t    count0 =  0;
 965   Int_t    count1 =  0;
 966   Float_t  rmsCalib   = rmsEvent;       // backup initial value ( from calib)
 967   Float_t  pedestalCalib = pedestalEvent;// backup initial value ( from calib)
 968   Int_t    firstBin = AliTPCReconstructor::GetRecoParam()->GetFirstBin();
 969   //
 970   UShort_t histo[kPedMax];
 971   memset(histo,0,kPedMax*sizeof(UShort_t));
 972   for (Int_t i=0; i<fMaxTime; i++){
 973     if (signal[i]<=0) continue;
 974     if (signal[i]>max && i>firstBin) {
 975       max = signal[i];
 976       maxPos = i;
 977     }
 978     if (signal[i]>kPedMax-1) continue;
 979     histo[int(signal[i]+0.5)]++;
 980     count0++;
 981   }
 982   //
 983   for (Int_t i=1; i<kPedMax; i++){
 984     if (count1<count0*0.5) median=i;
 985     count1+=histo[i];
 986   }
 987   // truncated mean
 988   //
 989   Float_t count10=histo[median] ,mean=histo[median]*median,  rms=histo[median]*median*median ;
 990   Float_t count06=histo[median] ,mean06=histo[median]*median,  rms06=histo[median]*median*median ;
 991   Float_t count09=histo[median] ,mean09=histo[median]*median,  rms09=histo[median]*median*median ;
 992   //
 993   for (Int_t idelta=1; idelta<10; idelta++){
 994     if (median-idelta<=0) continue;
 995     if (median+idelta>kPedMax) continue;
 996     if (count06<0.6*count1){
 997       count06+=histo[median-idelta];
 998       mean06 +=histo[median-idelta]*(median-idelta);
 999       rms06  +=histo[median-idelta]*(median-idelta)*(median-idelta);
1000       count06+=histo[median+idelta];
1001       mean06 +=histo[median+idelta]*(median+idelta);
1002       rms06  +=histo[median+idelta]*(median+idelta)*(median+idelta);
1003     }
1004     if (count09<0.9*count1){
1005       count09+=histo[median-idelta];
1006       mean09 +=histo[median-idelta]*(median-idelta);
1007       rms09  +=histo[median-idelta]*(median-idelta)*(median-idelta);
1008       count09+=histo[median+idelta];
1009       mean09 +=histo[median+idelta]*(median+idelta);
1010       rms09  +=histo[median+idelta]*(median+idelta)*(median+idelta);
1011     }
1012     if (count10<0.95*count1){
1013       count10+=histo[median-idelta];
1014       mean +=histo[median-idelta]*(median-idelta);
1015       rms  +=histo[median-idelta]*(median-idelta)*(median-idelta);
1016       count10+=histo[median+idelta];
1017       mean +=histo[median+idelta]*(median+idelta);
1018       rms  +=histo[median+idelta]*(median+idelta)*(median+idelta);
1019     }
1020   }
1021   mean  /=count10;
1022   mean06/=count06;
1023   mean09/=count09;
1024   rms    = TMath::Sqrt(TMath::Abs(rms/count10-mean*mean));
1025   rms06    = TMath::Sqrt(TMath::Abs(rms06/count06-mean06*mean06));
1026  rms09    = TMath::Sqrt(TMath::Abs(rms09/count09-mean09*mean09));
1027   rmsEvent = rms09;
1028   //
1029   pedestalEvent = median;
1030   if (AliLog::GetDebugLevel("","AliTPCclustererMI")==0) return median;
1031   //
1032   UInt_t uid[3] = {UInt_t(id[0]),UInt_t(id[1]),UInt_t(id[2])};
1033   //
1034   // Dump mean signal info
1035   //
1036   (*fDebugStreamer)<<"Signal"<<
1037     "TimeStamp="<<fTimeStamp<<
1038     "EventType="<<fEventType<<
1039     "Sector="<<uid[0]<<
1040     "Row="<<uid[1]<<
1041     "Pad="<<uid[2]<<
1042     "Max="<<max<<
1043     "MaxPos="<<maxPos<<
1044     //
1045     "Median="<<median<<
1046     "Mean="<<mean<<
1047     "RMS="<<rms<<
1048     "Mean06="<<mean06<<
1049     "RMS06="<<rms06<<
1050     "Mean09="<<mean09<<
1051     "RMS09="<<rms09<<
1052     "RMSCalib="<<rmsCalib<<
1053     "PedCalib="<<pedestalCalib<<
1054     "\n";
1055   //
1056   // fill pedestal histogram
1057   //
1058   AliTPCROC * roc = AliTPCROC::Instance();
1059   if (!fAmplitudeHisto){
1060     fAmplitudeHisto = new TObjArray(72);
1061   }
1062   //
1063   if (uid[0]<roc->GetNSectors()
1064       && uid[1]< roc->GetNRows(uid[0])  &&
1065       uid[2] <roc->GetNPads(uid[0], uid[1])){
1066     TObjArray  * sectorArray = (TObjArray*)fAmplitudeHisto->UncheckedAt(uid[0]);
1067     if (!sectorArray){
1068       Int_t npads =roc->GetNChannels(uid[0]);
1069       sectorArray = new TObjArray(npads);
1070       fAmplitudeHisto->AddAt(sectorArray, uid[0]);
1071     }
1072     Int_t position =  uid[2]+roc->GetRowIndexes(uid[0])[uid[1]];
1073     TH1F * histo = (TH1F*)sectorArray->UncheckedAt(position);
1074     if (!histo){
1075       char hname[100];
1076       sprintf(hname,"Amp_%d_%d_%d",uid[0],uid[1],uid[2]);
1077       TFile * backup = gFile;
1078       fDebugStreamer->GetFile()->cd();
1079       histo = new TH1F(hname, hname, 100, 5,100);
1080       //histo->SetDirectory(0);     // histogram not connected to directory -(File)
1081       sectorArray->AddAt(histo, position);
1082       if (backup) backup->cd();
1083     }
1084     for (Int_t i=0; i<nchannels; i++){
1085       histo->Fill(signal[i]);
1086     }
1087   }
1088   //
1089   //
1090   //
1091   Float_t kMin =fRecoParam->GetDumpAmplitudeMin();   // minimal signal to be dumped
1092   Float_t *dsignal = new Float_t[nchannels];
1093   Float_t *dtime   = new Float_t[nchannels];
1094   for (Int_t i=0; i<nchannels; i++){
1095     dtime[i] = i;
1096     dsignal[i] = signal[i];
1097   }
1098   //
1099   // Digital noise
1100   //
1101  //  if (max-median>30.*TMath::Max(1.,Double_t(rms06)) &&  (((*fDebugStreamer)<<"SignalDN").GetSize()<kMaxDebugSize)){
1102 //     //
1103 //     //
1104 //     TGraph * graph =new TGraph(nchannels, dtime, dsignal);
1105 //     //
1106 //     //
1107 //     // jumps left - right
1108 //     Int_t    njumps0=0;
1109 //     Double_t deltaT0[2000];
1110 //     Double_t deltaA0[2000];
1111 //     Int_t    lastJump0 = fRecoParam->GetFirstBin();
1112 //     Int_t    njumps1=0;
1113 //     Double_t deltaT1[2000];
1114 //     Double_t deltaA1[2000];
1115 //     Int_t    lastJump1 = fRecoParam->GetFirstBin();
1116 //     Int_t    njumps2=0;
1117 //     Double_t deltaT2[2000];
1118 //     Double_t deltaA2[2000];
1119 //     Int_t    lastJump2 = fRecoParam->GetFirstBin();
1120
1121 //     for (Int_t itime=fRecoParam->GetFirstBin()+1; itime<fRecoParam->GetLastBin()-1; itime++){
1122 //       if (TMath::Abs(dsignal[itime]-dsignal[itime-1])>30.*TMath::Max(1.,Double_t(rms06))  &&
1123 //        TMath::Abs(dsignal[itime]-dsignal[itime+1])>30.*TMath::Max(1.,Double_t(rms06))  &&
1124 //        (dsignal[itime-1]-median<5.*rms06) &&
1125 //        (dsignal[itime+1]-median<5.*rms06)
1126 //        ){
1127 //      deltaA0[njumps0] = dsignal[itime]-dsignal[itime-1];
1128 //      deltaT0[njumps0] = itime-lastJump0;
1129 //      lastJump0 = itime;
1130 //      njumps0++;
1131 //       }
1132 //       if (TMath::Abs(dsignal[itime]-dsignal[itime-1])>30.*TMath::Max(1.,Double_t(rms06)) &&
1133 //        (dsignal[itime-1]-median<5.*rms06)
1134 //        ) {
1135 //      deltaA1[njumps1] = dsignal[itime]-dsignal[itime-1];
1136 //      deltaT1[njumps1] = itime-lastJump1;
1137 //      lastJump1 = itime;
1138 //      njumps1++;
1139 //       }
1140 //       if (TMath::Abs(dsignal[itime]-dsignal[itime+1])>30.*TMath::Max(1.,Double_t(rms06)) &&
1141 //        (dsignal[itime+1]-median<5.*rms06)
1142 //        ) {
1143 //      deltaA2[njumps2] = dsignal[itime]-dsignal[itime+1];
1144 //      deltaT2[njumps2] = itime-lastJump2;
1145 //      lastJump2 = itime;
1146 //      njumps2++;
1147 //       }
1148 //     }
1149 //     //
1150 //     if (njumps0>0 || njumps1>0 || njumps2>0){
1151 //       TGraph *graphDN0 = new TGraph(njumps0, deltaT0, deltaA0);
1152 //       TGraph *graphDN1 = new TGraph(njumps1, deltaT1, deltaA1);
1153 //       TGraph *graphDN2 = new TGraph(njumps2, deltaT2, deltaA2);
1154 //       (*fDebugStreamer)<<"SignalDN"<<    //digital - noise pads - or random sample of pads
1155 //      "TimeStamp="<<fTimeStamp<<
1156 //      "EventType="<<fEventType<<
1157 //      "Sector="<<uid[0]<<
1158 //      "Row="<<uid[1]<<
1159 //      "Pad="<<uid[2]<<
1160 //      "Graph="<<graph<<
1161 //      "Max="<<max<<
1162 //      "MaxPos="<<maxPos<<
1163 //      "Graph.="<<graph<<
1164 //      "P0GraphDN0.="<<graphDN0<<
1165 //      "P1GraphDN1.="<<graphDN1<<
1166 //      "P2GraphDN2.="<<graphDN2<<
1167 //      //
1168 //      "Median="<<median<<
1169 //      "Mean="<<mean<<
1170 //      "RMS="<<rms<<
1171 //      "Mean06="<<mean06<<
1172 //      "RMS06="<<rms06<<
1173 //      "Mean09="<<mean09<<
1174 //      "RMS09="<<rms09<<
1175 //      "\n";
1176 //       delete graphDN0;
1177 //       delete graphDN1;
1178 //       delete graphDN2;
1179 //     }
1180 //     delete graph;
1181 //   }
1182
1183   //
1184   // NOISE STUDY  Fourier transform
1185   //
1186   TGraph * graph;
1187   Bool_t random = (gRandom->Rndm()<0.0003);
1188   if (((*fDebugStreamer)<<"SignalN").GetSize()<kMaxDebugSize)
1189     if (max-median>kMin || rms06>1.*fParam->GetZeroSup() || random){
1190     graph =new TGraph(nchannels, dtime, dsignal);
1191     if (rms06>1.*fParam->GetZeroSup() || random){
1192       //Double_t *input, Double_t threshold, Bool_t locMax, Double_t *freq, Double_t *re, Double_t *im, Double_t *mag, Double_t *phi);
1193       Float_t * input = &(dsignal[fRecoParam->GetFirstBin()]);
1194       Float_t freq[2000], re[2000], im[2000], mag[2000], phi[2000];
1195       Int_t npoints = TransformFFT(input, -1,kFALSE, freq, re, im, mag, phi);
1196       TGraph *graphMag0 = new TGraph(npoints, freq, mag);
1197       TGraph *graphPhi0 = new TGraph(npoints, freq, phi);
1198       npoints = TransformFFT(input, 0.5,kTRUE, freq, re, im, mag, phi);
1199       TGraph *graphMag1 = new TGraph(npoints, freq, mag);
1200       TGraph *graphPhi1 = new TGraph(npoints, freq, phi);
1201
1202       (*fDebugStreamer)<<"SignalN"<<    //noise pads - or random sample of pads
1203         "TimeStamp="<<fTimeStamp<<
1204         "EventType="<<fEventType<<
1205         "Sector="<<uid[0]<<
1206         "Row="<<uid[1]<<
1207         "Pad="<<uid[2]<<
1208         "Graph.="<<graph<<
1209         "Max="<<max<<
1210         "MaxPos="<<maxPos<<
1211         //
1212         "Median="<<median<<
1213         "Mean="<<mean<<
1214         "RMS="<<rms<<
1215         "Mean06="<<mean06<<
1216         "RMS06="<<rms06<<
1217         "Mean09="<<mean09<<
1218         "RMS09="<<rms09<<
1219         // FFT part
1220         "Mag0.="<<graphMag0<<
1221         "Mag1.="<<graphMag1<<
1222         "Phi0.="<<graphPhi0<<
1223         "Phi1.="<<graphPhi1<<
1224         "\n";
1225       delete graphMag0;
1226       delete graphMag1;
1227       delete graphPhi0;
1228       delete graphPhi1;
1229     }
1230     //
1231     // Big signals dumping
1232     //
1233
1234     if (max-median>kMin &&maxPos>AliTPCReconstructor::GetRecoParam()->GetFirstBin())
1235       (*fDebugStreamer)<<"SignalB"<<     // pads with signal
1236         "TimeStamp="<<fTimeStamp<<
1237         "EventType="<<fEventType<<
1238         "Sector="<<uid[0]<<
1239         "Row="<<uid[1]<<
1240         "Pad="<<uid[2]<<
1241         "Graph="<<graph<<
1242         "Max="<<max<<
1243         "MaxPos="<<maxPos<<
1244         //
1245         "Median="<<median<<
1246         "Mean="<<mean<<
1247         "RMS="<<rms<<
1248         "Mean06="<<mean06<<
1249         "RMS06="<<rms06<<
1250         "Mean09="<<mean09<<
1251         "RMS09="<<rms09<<
1252         "\n";
1253     delete graph;
1254   }
1255
1256   //
1257   //
1258   //  Central Electrode signal analysis
1259   //
1260   Float_t ceQmax  =0, ceQsum=0, ceTime=0;
1261   Float_t cemean  = mean06, cerms=rms06 ;
1262   Int_t    cemaxpos= 0;
1263   Float_t ceThreshold=5.*cerms;
1264   Float_t ceSumThreshold=8.*cerms;
1265   const Int_t    kCemin=5;  // range for the analysis of the ce signal +- channels from the peak
1266   const Int_t    kCemax=5;
1267   for (Int_t i=nchannels-2; i>nchannels/2; i--){
1268     if ( (dsignal[i]-mean06)>ceThreshold && dsignal[i]>=dsignal[i+1] && dsignal[i]>=dsignal[i-1] ){
1269       cemaxpos=i;
1270       break;
1271     }
1272   }
1273   if (cemaxpos!=0){
1274     ceQmax = 0;
1275     Int_t cemaxpos2=0;
1276     for (Int_t i=cemaxpos-20; i<cemaxpos+5; i++){
1277       if (i<0 || i>nchannels-1) continue;
1278       Double_t val=dsignal[i]- cemean;
1279       if (val>ceQmax){
1280         cemaxpos2=i;
1281         ceQmax = val;
1282       }
1283     }
1284     cemaxpos = cemaxpos2;
1285
1286     for (Int_t i=cemaxpos-kCemin; i<cemaxpos+kCemax; i++){
1287       if (i>0 && i<nchannels&&dsignal[i]- cemean>0){
1288         Double_t val=dsignal[i]- cemean;
1289         ceTime+=val*dtime[i];
1290         ceQsum+=val;
1291         if (val>ceQmax) ceQmax=val;
1292       }
1293     }
1294     if (ceQmax&&ceQsum>ceSumThreshold) {
1295       ceTime/=ceQsum;
1296       (*fDebugStreamer)<<"Signalce"<<
1297         "TimeStamp="<<fTimeStamp<<
1298         "EventType="<<fEventType<<
1299         "Sector="<<uid[0]<<
1300         "Row="<<uid[1]<<
1301         "Pad="<<uid[2]<<
1302         "Max="<<ceQmax<<
1303         "Qsum="<<ceQsum<<
1304         "Time="<<ceTime<<
1305         "RMS06="<<rms06<<
1306         //
1307         "\n";
1308     }
1309   }
1310   // end of ce signal analysis
1311   //
1312
1313   //
1314   //  Gating grid signal analysis
1315   //
1316   Double_t ggQmax  =0, ggQsum=0, ggTime=0;
1317   Double_t ggmean  = mean06, ggrms=rms06 ;
1318   Int_t    ggmaxpos= 0;
1319   Double_t ggThreshold=5.*ggrms;
1320   Double_t ggSumThreshold=8.*ggrms;
1321
1322   for (Int_t i=1; i<nchannels/4; i++){
1323     if ( (dsignal[i]-mean06)>ggThreshold && dsignal[i]>=dsignal[i+1] && dsignal[i]>=dsignal[i-1] &&
1324          (dsignal[i]+dsignal[i+1]+dsignal[i-1]-3*mean06)>ggSumThreshold){
1325       ggmaxpos=i;
1326       if (dsignal[i-1]>dsignal[i+1]) ggmaxpos=i-1;
1327       break;
1328     }
1329   }
1330   if (ggmaxpos!=0){
1331       for (Int_t i=ggmaxpos-1; i<ggmaxpos+3; i++){
1332           if (i>0 && i<nchannels && dsignal[i]-ggmean>0){
1333               Double_t val=dsignal[i]- ggmean;
1334               ggTime+=val*dtime[i];
1335               ggQsum+=val;
1336               if (val>ggQmax) ggQmax=val;
1337           }
1338       }
1339       if (ggQmax&&ggQsum>ggSumThreshold) {
1340           ggTime/=ggQsum;
1341           (*fDebugStreamer)<<"Signalgg"<<
1342             "TimeStamp="<<fTimeStamp<<
1343             "EventType="<<fEventType<<
1344               "Sector="<<uid[0]<<
1345               "Row="<<uid[1]<<
1346               "Pad="<<uid[2]<<
1347               "Max="<<ggQmax<<
1348               "Qsum="<<ggQsum<<
1349               "Time="<<ggTime<<
1350               "RMS06="<<rms06<<
1351               //
1352               "\n";
1353       }
1354   }
1355   // end of gg signal analysis
1356
1357
1358   delete [] dsignal;
1359   delete [] dtime;
1360   if (rms06>fRecoParam->GetMaxNoise()) {
1361     pedestalEvent+=1024.;
1362     return 1024+median; // sign noisy channel in debug mode
1363   }
1364   return median;
1365 }
1366
1367
1368
1369 void AliTPCclustererMI::DumpHistos(){
1370   //
1371   // Dump histogram information
1372   //
1373   if (!fAmplitudeHisto) return;
1374   AliTPCROC * roc = AliTPCROC::Instance();
1375   for (UInt_t isector=0; isector<AliTPCROC::Instance()->GetNSectors(); isector++){
1376     TObjArray * array = (TObjArray*)fAmplitudeHisto->UncheckedAt(isector);
1377     if (!array) continue;
1378     for (UInt_t ipad = 0; ipad <(UInt_t)array->GetEntriesFast(); ipad++){
1379       TH1F * histo = (TH1F*) array->UncheckedAt(ipad);
1380       if (!histo) continue;
1381       if (histo->GetEntries()<100) continue;
1382       histo->Fit("gaus","q");
1383       Float_t mean =  histo->GetMean();
1384       Float_t rms  =  histo->GetRMS();
1385       Float_t gmean = histo->GetFunction("gaus")->GetParameter(1);
1386       Float_t gsigma = histo->GetFunction("gaus")->GetParameter(2);
1387       Float_t gmeanErr = histo->GetFunction("gaus")->GetParError(1);
1388       Float_t gsigmaErr = histo->GetFunction("gaus")->GetParError(2);
1389       Float_t max = histo->GetFunction("gaus")->GetParameter(0);
1390
1391       // get pad number
1392       UInt_t row=0, pad =0;
1393       const UInt_t *indexes =roc->GetRowIndexes(isector);
1394       for (UInt_t irow=0; irow<roc->GetNRows(isector); irow++){
1395         if (indexes[irow]<=ipad){
1396           row = irow;
1397           pad = ipad-indexes[irow];
1398         }
1399       }
1400       Int_t rpad = pad - (AliTPCROC::Instance()->GetNPads(isector,row))/2;
1401       //
1402       (*fDebugStreamer)<<"Fit"<<
1403         "TimeStamp="<<fTimeStamp<<
1404         "EventType="<<fEventType<<
1405         "Sector="<<isector<<
1406         "Row="<<row<<
1407         "Pad="<<pad<<
1408         "RPad="<<rpad<<
1409         "Max="<<max<<
1410         "Mean="<<mean<<
1411         "RMS="<<rms<<
1412         "GMean="<<gmean<<
1413         "GSigma="<<gsigma<<
1414         "GMeanErr="<<gmeanErr<<
1415         "GSigmaErr="<<gsigmaErr<<
1416         "\n";
1417       if (array->UncheckedAt(ipad)) fDebugStreamer->StoreObject(array->UncheckedAt(ipad));
1418     }
1419   }
1420 }
1421
1422
1423
1424 Int_t  AliTPCclustererMI::TransformFFT(Float_t *input, Float_t threshold, Bool_t locMax, Float_t *freq, Float_t *re, Float_t *im, Float_t *mag, Float_t *phi)
1425 {
1426   //
1427   // calculate fourrie transform
1428   // return only frequncies with mag over threshold
1429   // if locMax is spectified only freque with local maxima over theshold is returned
1430
1431   if (! fFFTr2c) return kFALSE;
1432   if (!freq) return kFALSE;
1433
1434   Int_t current=0;
1435   Int_t nPoints = fRecoParam->GetLastBin()-fRecoParam->GetFirstBin();
1436   Double_t *in = new Double_t[nPoints];
1437   Double_t *rfft = new Double_t[nPoints];
1438   Double_t *ifft = new Double_t[nPoints];
1439   for (Int_t i=0; i<nPoints; i++){in[i]=input[i];}
1440   fFFTr2c->SetPoints(in);
1441   fFFTr2c->Transform();
1442   fFFTr2c->GetPointsComplex(rfft, ifft);
1443   for (Int_t i=3; i<nPoints/2-3; i++){
1444     Float_t lmag =  TMath::Sqrt(rfft[i]*rfft[i]+ifft[i]*ifft[i])/nPoints;
1445     if (lmag<threshold) continue;
1446     if (locMax){
1447       if ( TMath::Sqrt(rfft[i-1]*rfft[i-1]+ifft[i-1]*ifft[i-1])/nPoints>lmag) continue;
1448       if ( TMath::Sqrt(rfft[i+1]*rfft[i+1]+ifft[i+1]*ifft[i+1])/nPoints>lmag) continue;
1449       if ( TMath::Sqrt(rfft[i-2]*rfft[i-2]+ifft[i-2]*ifft[i-2])/nPoints>lmag) continue;
1450       if ( TMath::Sqrt(rfft[i+2]*rfft[i+2]+ifft[i+2]*ifft[i+2])/nPoints>lmag) continue;
1451       if ( TMath::Sqrt(rfft[i-3]*rfft[i-3]+ifft[i-3]*ifft[i-3])/nPoints>lmag) continue;
1452       if ( TMath::Sqrt(rfft[i+3]*rfft[i+3]+ifft[i+3]*ifft[i+3])/nPoints>lmag) continue;
1453     }
1454
1455     freq[current] = Float_t(i)/Float_t(nPoints);
1456     //
1457     re[current] = rfft[i];
1458     im[current] = ifft[i];
1459     mag[current]=lmag;
1460     phi[current]=TMath::ATan2(ifft[i],rfft[i]);
1461     current++;
1462   }
1463   delete [] in;
1464   delete [] rfft;
1465   delete [] ifft;
1466   return current;
1467 }
1468