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 (fRecoParam->DumpSignal() &&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   if (!fRecoParam->DumpSignal()) {
 580     cl->SetInfo(0);
 581   }
 582
 583   fNcluster++;
 584 }
 585
 586
 587 //_____________________________________________________________________________
 588 void AliTPCclustererMI::Digits2Clusters()
 589 {
 590   //-----------------------------------------------------------------
 591   // This is a simple cluster finder.
 592   //-----------------------------------------------------------------
 593
 594   if (!fInput) {
 595     Error("Digits2Clusters", "input tree not initialised");
 596     return;
 597   }
 598
 599   if (!fOutput) {
 600     Error("Digits2Clusters", "output tree not initialised");
 601     return;
 602   }
 603
 604   AliTPCCalPad * gainTPC = AliTPCcalibDB::Instance()->GetPadGainFactor();
 605   AliTPCCalPad * noiseTPC = AliTPCcalibDB::Instance()->GetPadNoise();
 606
 607   AliSimDigits digarr, *dummy=&digarr;
 608   fRowDig = dummy;
 609   fInput->GetBranch("Segment")->SetAddress(&dummy);
 610   Stat_t nentries = fInput->GetEntries();
 611
 612   fMaxTime=fParam->GetMaxTBin()+6; // add 3 virtual time bins before and 3   after
 613
 614   Int_t nclusters  = 0;
 615
 616   for (Int_t n=0; n<nentries; n++) {
 617     fInput->GetEvent(n);
 618     if (!fParam->AdjustSectorRow(digarr.GetID(),fSector,fRow)) {
 619       cerr<<"AliTPC warning: invalid segment ID ! "<<digarr.GetID()<<endl;
 620       continue;
 621     }
 622     Int_t row = fRow;
 623     AliTPCCalROC * gainROC = gainTPC->GetCalROC(fSector);  // pad gains per given sector
 624     AliTPCCalROC * noiseROC   = noiseTPC->GetCalROC(fSector); // noise per given sector
 625     //
 626     AliTPCClustersRow *clrow= new AliTPCClustersRow();
 627     fRowCl = clrow;
 628     clrow->SetClass("AliTPCclusterMI");
 629     clrow->SetArray(1);
 630
 631     clrow->SetID(digarr.GetID());
 632     fOutput->GetBranch("Segment")->SetAddress(&clrow);
 633     fRx=fParam->GetPadRowRadii(fSector,row);
 634
 635
 636     const Int_t kNIS=fParam->GetNInnerSector(), kNOS=fParam->GetNOuterSector();
 637     fZWidth = fParam->GetZWidth();
 638     if (fSector < kNIS) {
 639       fMaxPad = fParam->GetNPadsLow(row);
 640       fSign = (fSector < kNIS/2) ? 1 : -1;
 641       fPadLength = fParam->GetPadPitchLength(fSector,row);
 642       fPadWidth = fParam->GetPadPitchWidth();
 643     } else {
 644       fMaxPad = fParam->GetNPadsUp(row);
 645       fSign = ((fSector-kNIS) < kNOS/2) ? 1 : -1;
 646       fPadLength = fParam->GetPadPitchLength(fSector,row);
 647       fPadWidth  = fParam->GetPadPitchWidth();
 648     }
 649
 650
 651     fMaxBin=fMaxTime*(fMaxPad+6);  // add 3 virtual pads  before and 3 after
 652     fBins    =new Float_t[fMaxBin];
 653     fSigBins =new Int_t[fMaxBin];
 654     fNSigBins = 0;
 655     memset(fBins,0,sizeof(Float_t)*fMaxBin);
 656
 657     if (digarr.First()) //MI change
 658       do {
 659         Float_t dig=digarr.CurrentDigit();
 660         if (dig<=fParam->GetZeroSup()) continue;
 661         Int_t j=digarr.CurrentRow()+3, i=digarr.CurrentColumn()+3;
 662         Float_t gain = gainROC->GetValue(row,digarr.CurrentColumn());
 663         Int_t bin = i*fMaxTime+j;
 664         fBins[bin]=dig/gain;
 665         fSigBins[fNSigBins++]=bin;
 666       } while (digarr.Next());
 667     digarr.ExpandTrackBuffer();
 668
 669     FindClusters(noiseROC);
 670
 671     fOutput->Fill();
 672     delete clrow;
 673     nclusters+=fNcluster;
 674     delete[] fBins;
 675     delete[] fSigBins;
 676   }
 677
 678   Info("Digits2Clusters", "Number of found clusters : %d", nclusters);
 679 }
 680
 681 void AliTPCclustererMI::Digits2Clusters(AliRawReader* rawReader)
 682 {
 683 //-----------------------------------------------------------------
 684 // This is a cluster finder for the TPC raw data.
 685 // The method assumes NO ordering of the altro channels.
 686 // The pedestal subtraction can be switched on and off
 687 // using an option of the TPC reconstructor
 688 //-----------------------------------------------------------------
 689
 690   if (!fOutput) {
 691     Error("Digits2Clusters", "output tree not initialised");
 692     return;
 693   }
 694
 695   fRowDig = NULL;
 696   AliTPCROC * roc = AliTPCROC::Instance();
 697   AliTPCCalPad * gainTPC = AliTPCcalibDB::Instance()->GetPadGainFactor();
 698   AliTPCCalPad * pedestalTPC = AliTPCcalibDB::Instance()->GetPedestals();
 699   AliTPCCalPad * noiseTPC = AliTPCcalibDB::Instance()->GetPadNoise();
 700   AliTPCRawStream input(rawReader);
 701   fEventHeader = (AliRawEventHeaderBase*)rawReader->GetEventHeader();
 702   if (fEventHeader){
 703     fTimeStamp = fEventHeader->Get("Timestamp");
 704     fEventType = fEventHeader->Get("Type");
 705   }
 706
 707
 708   Int_t nclusters  = 0;
 709
 710   fMaxTime = fParam->GetMaxTBin() + 6; // add 3 virtual time bins before and 3 after
 711   const Int_t kNIS = fParam->GetNInnerSector();
 712   const Int_t kNOS = fParam->GetNOuterSector();
 713   const Int_t kNS = kNIS + kNOS;
 714   fZWidth = fParam->GetZWidth();
 715   Int_t zeroSup = fParam->GetZeroSup();
 716   //
 717   //alocate memory for sector - maximal case
 718   //
 719   Float_t** allBins = NULL;
 720   Int_t** allSigBins = NULL;
 721   Int_t*  allNSigBins = NULL;
 722   Int_t nRowsMax = roc->GetNRows(roc->GetNSector()-1);
 723   Int_t nPadsMax = roc->GetNPads(roc->GetNSector()-1,nRowsMax-1);
 724   allBins = new Float_t*[nRowsMax];
 725   allSigBins = new Int_t*[nRowsMax];
 726   allNSigBins = new Int_t[nRowsMax];
 727   for (Int_t iRow = 0; iRow < nRowsMax; iRow++) {
 728     //
 729     Int_t maxBin = fMaxTime*(nPadsMax+6);  // add 3 virtual pads  before and 3 after
 730     allBins[iRow] = new Float_t[maxBin];
 731     memset(allBins[iRow],0,sizeof(Float_t)*maxBin);
 732     allSigBins[iRow] = new Int_t[maxBin];
 733     allNSigBins[iRow]=0;
 734   }
 735   //
 736   // Loop over sectors
 737   //
 738   for(fSector = 0; fSector < kNS; fSector++) {
 739
 740     AliTPCCalROC * gainROC    = gainTPC->GetCalROC(fSector);  // pad gains per given sector
 741     AliTPCCalROC * pedestalROC = pedestalTPC->GetCalROC(fSector);  // pedestal per given sector
 742     AliTPCCalROC * noiseROC   = noiseTPC->GetCalROC(fSector);  // noise per given sector
 743
 744     Int_t nRows = 0;
 745     Int_t nDDLs = 0, indexDDL = 0;
 746     if (fSector < kNIS) {
 747       nRows = fParam->GetNRowLow();
 748       fSign = (fSector < kNIS/2) ? 1 : -1;
 749       nDDLs = 2;
 750       indexDDL = fSector * 2;
 751     }
 752     else {
 753       nRows = fParam->GetNRowUp();
 754       fSign = ((fSector-kNIS) < kNOS/2) ? 1 : -1;
 755       nDDLs = 4;
 756       indexDDL = (fSector-kNIS) * 4 + kNIS * 2;
 757     }
 758
 759     for (Int_t iRow = 0; iRow < nRows; iRow++) {
 760       Int_t maxPad;
 761       if (fSector < kNIS)
 762         maxPad = fParam->GetNPadsLow(iRow);
 763       else
 764         maxPad = fParam->GetNPadsUp(iRow);
 765
 766       Int_t maxBin = fMaxTime*(maxPad+6);  // add 3 virtual pads  before and 3 after
 767       memset(allBins[iRow],0,sizeof(Float_t)*maxBin);
 768       allNSigBins[iRow] = 0;
 769     }
 770
 771     // Loas the raw data for corresponding DDLs
 772     rawReader->Reset();
 773     input.SetOldRCUFormat(fIsOldRCUFormat);
 774     rawReader->Select("TPC",indexDDL,indexDDL+nDDLs-1);
 775     Int_t digCounter=0;
 776     // Begin loop over altro data
 777     Bool_t calcPedestal = fRecoParam->GetCalcPedestal();
 778     Float_t gain =1;
 779     Int_t lastPad=-1;
 780     while (input.Next()) {
 781       digCounter++;
 782       if (input.GetSector() != fSector)
 783         AliFatal(Form("Sector index mismatch ! Expected (%d), but got (%d) !",fSector,input.GetSector()));
 784
 785
 786       Int_t iRow = input.GetRow();
 787       if (iRow < 0 || iRow >= nRows)
 788         AliFatal(Form("Pad-row index (%d) outside the range (%d -> %d) !",
 789                       iRow, 0, nRows -1));
 790       //pad
 791       Int_t iPad = input.GetPad();
 792       if (iPad < 0 || iPad >= nPadsMax)
 793         AliFatal(Form("Pad index (%d) outside the range (%d -> %d) !",
 794                       iPad, 0, nPadsMax-1));
 795       if (iPad!=lastPad){
 796         gain    = gainROC->GetValue(iRow,iPad);
 797         lastPad = iPad;
 798       }
 799       iPad+=3;
 800       //time
 801       Int_t iTimeBin = input.GetTime();
 802       if ( iTimeBin < 0 || iTimeBin >= fParam->GetMaxTBin())
 803         AliFatal(Form("Timebin index (%d) outside the range (%d -> %d) !",
 804                       iTimeBin, 0, iTimeBin -1));
 805       iTimeBin+=3;
 806       //signal
 807       Float_t signal = input.GetSignal();
 808       if (!calcPedestal && signal <= zeroSup) continue;
 809       if (!calcPedestal) {
 810         Int_t bin = iPad*fMaxTime+iTimeBin;
 811         allBins[iRow][bin] = signal/gain;
 812         allSigBins[iRow][allNSigBins[iRow]++] = bin;
 813       }else{
 814         allBins[iRow][iPad*fMaxTime+iTimeBin] = signal;
 815       }
 816       allBins[iRow][iPad*fMaxTime+0]=1.;  // pad with signal
 817     } // End of the loop over altro data
 818     //
 819     //
 820     // Now loop over rows and perform pedestal subtraction
 821     if (digCounter==0) continue;
 822     //    if (fPedSubtraction) {
 823     if (calcPedestal) {
 824       for (Int_t iRow = 0; iRow < nRows; iRow++) {
 825         Int_t maxPad;
 826         if (fSector < kNIS)
 827           maxPad = fParam->GetNPadsLow(iRow);
 828         else
 829           maxPad = fParam->GetNPadsUp(iRow);
 830
 831         for (Int_t iPad = 3; iPad < maxPad + 3; iPad++) {
 832           if (allBins[iRow][iPad*fMaxTime+0] <1 ) continue;  // no data
 833           Float_t *p = &allBins[iRow][iPad*fMaxTime+3];
 834           //Float_t pedestal = TMath::Median(fMaxTime, p);
 835           Int_t id[3] = {fSector, iRow, iPad-3};
 836           // calib values
 837           Double_t rmsCalib=  noiseROC->GetValue(iRow,iPad-3);
 838           Double_t pedestalCalib = pedestalROC->GetValue(iRow,iPad-3);
 839           Double_t rmsEvent       = rmsCalib;
 840           Double_t pedestalEvent  = pedestalCalib;
 841           ProcesSignal(p, fMaxTime, id, rmsEvent, pedestalEvent);
 842           if (rmsEvent<rmsCalib) rmsEvent = rmsCalib;   // take worst scenario
 843           if (TMath::Abs(pedestalEvent-pedestalCalib)<1.0) pedestalEvent = pedestalCalib;
 844
 845           //
 846           for (Int_t iTimeBin = 0; iTimeBin < fMaxTime; iTimeBin++) {
 847             Int_t bin = iPad*fMaxTime+iTimeBin;
 848             allBins[iRow][bin] -= pedestalEvent;
 849             if (iTimeBin < AliTPCReconstructor::GetRecoParam()->GetFirstBin())
 850               allBins[iRow][bin] = 0;
 851             if (iTimeBin > AliTPCReconstructor::GetRecoParam()->GetLastBin())
 852               allBins[iRow][bin] = 0;
 853             if (allBins[iRow][iPad*fMaxTime+iTimeBin] < zeroSup)
 854               allBins[iRow][bin] = 0;
 855             if (allBins[iRow][bin] < 3.0*rmsEvent)   // 3 sigma cut on RMS
 856               allBins[iRow][bin] = 0;
 857             if (allBins[iRow][bin]) allSigBins[iRow][allNSigBins[iRow]++] = bin;
 858           }
 859         }
 860       }
 861     }
 862     // Now loop over rows and find clusters
 863     for (fRow = 0; fRow < nRows; fRow++) {
 864       fRowCl = new AliTPCClustersRow;
 865       fRowCl->SetClass("AliTPCclusterMI");
 866       fRowCl->SetArray(1);
 867       fRowCl->SetID(fParam->GetIndex(fSector, fRow));
 868       fOutput->GetBranch("Segment")->SetAddress(&fRowCl);
 869
 870       fRx = fParam->GetPadRowRadii(fSector, fRow);
 871       fPadLength = fParam->GetPadPitchLength(fSector, fRow);
 872       fPadWidth  = fParam->GetPadPitchWidth();
 873       if (fSector < kNIS)
 874         fMaxPad = fParam->GetNPadsLow(fRow);
 875       else
 876         fMaxPad = fParam->GetNPadsUp(fRow);
 877       fMaxBin = fMaxTime*(fMaxPad+6);  // add 3 virtual pads  before and 3 after
 878
 879       fBins = allBins[fRow];
 880       fSigBins = allSigBins[fRow];
 881       fNSigBins = allNSigBins[fRow];
 882
 883       FindClusters(noiseROC);
 884
 885       fOutput->Fill();
 886       delete fRowCl;
 887       nclusters += fNcluster;
 888     } // End of loop to find clusters
 889   } // End of loop over sectors
 890
 891   for (Int_t iRow = 0; iRow < nRowsMax; iRow++) {
 892     delete [] allBins[iRow];
 893     delete [] allSigBins[iRow];
 894   }
 895   delete [] allBins;
 896   delete [] allSigBins;
 897   delete [] allNSigBins;
 898
 899   Info("Digits2Clusters", "File  %s Event\t%d\tNumber of found clusters : %d\n", fOutput->GetName(),*(rawReader->GetEventId()), nclusters);
 900
 901 }
 902
 903 void AliTPCclustererMI::FindClusters(AliTPCCalROC * noiseROC)
 904 {
 905
 906   //
 907   // add virtual charge at the edge
 908   //
 909   Double_t kMaxDumpSize = 500000;
 910   if (fRecoParam->GetCalcPedestal() && fOutput->GetZipBytes()< kMaxDumpSize) fBDumpSignal =kTRUE;   //dump signal flag
 911   //
 912   fNcluster=0;
 913   fLoop=1;
 914   Int_t crtime = Int_t((fParam->GetZLength(fSector)-fRecoParam->GetCtgRange()*fRx)/fZWidth-fParam->GetNTBinsL1()-5);
 915   Float_t minMaxCutAbs       = fRecoParam->GetMinMaxCutAbs();
 916   Float_t minLeftRightCutAbs = fRecoParam->GetMinLeftRightCutAbs();
 917   Float_t minUpDownCutAbs    = fRecoParam->GetMinUpDownCutAbs();
 918   Float_t minMaxCutSigma       = fRecoParam->GetMinMaxCutSigma();
 919   Float_t minLeftRightCutSigma = fRecoParam->GetMinLeftRightCutSigma();
 920   Float_t minUpDownCutSigma    = fRecoParam->GetMinUpDownCutSigma();
 921   for (Int_t iSig = 0; iSig < fNSigBins; iSig++) {
 922     Int_t i = fSigBins[iSig];
 923     if (i%fMaxTime<=crtime) continue;
 924     Float_t *b = &fBins[i];
 925     //absolute custs
 926     if (b[0]<minMaxCutAbs) continue;   //threshold for maxima
 927     //
 928     if (b[-1]+b[1]+b[-fMaxTime]+b[fMaxTime]<=0) continue;  // cut on isolated clusters
 929     //    if (b[-1]+b[1]<=0) continue;               // cut on isolated clusters
 930     //if (b[-fMaxTime]+b[fMaxTime]<=0) continue; // cut on isolated clusters
 931     //
 932     if ((b[0]+b[-1]+b[1])<minUpDownCutAbs) continue;   //threshold for up down  (TRF)
 933     if ((b[0]+b[-fMaxTime]+b[fMaxTime])<minLeftRightCutAbs) continue;   //threshold for left right (PRF)
 934     if (!IsMaximum(*b,fMaxTime,b)) continue;
 935     //
 936     Float_t noise = noiseROC->GetValue(fRow, i/fMaxTime);
 937     // sigma cuts
 938     if (b[0]<minMaxCutSigma*noise) continue;   //threshold form maxima
 939     if ((b[0]+b[-1]+b[1])<minUpDownCutSigma*noise) continue;   //threshold for up town TRF
 940     if ((b[0]+b[-fMaxTime]+b[fMaxTime])<minLeftRightCutSigma*noise) continue;   //threshold for left right (PRF)
 941
 942     AliTPCclusterMI c(kFALSE);   // default cosntruction  without info
 943     Int_t dummy=0;
 944     MakeCluster(i, fMaxTime, fBins, dummy,c);
 945
 946     //}
 947   }
 948 }
 949
 950
 951 Double_t AliTPCclustererMI::ProcesSignal(Float_t *signal, Int_t nchannels, Int_t id[3], Double_t &rmsEvent, Double_t &pedestalEvent){
 952   //
 953   // process signal on given pad - + streaming of additional information in special mode
 954   //
 955   // id[0] - sector
 956   // id[1] - row
 957   // id[2] - pad
 958
 959   //
 960   // ESTIMATE pedestal and the noise
 961   //
 962   const Int_t kPedMax = 100;
 963   Double_t kMaxDebugSize = 5000000.;
 964   Float_t  max    =  0;
 965   Float_t  maxPos =  0;
 966   Int_t    median =  -1;
 967   Int_t    count0 =  0;
 968   Int_t    count1 =  0;
 969   Float_t  rmsCalib   = rmsEvent;       // backup initial value ( from calib)
 970   Float_t  pedestalCalib = pedestalEvent;// backup initial value ( from calib)
 971   Int_t    firstBin = AliTPCReconstructor::GetRecoParam()->GetFirstBin();
 972   //
 973   UShort_t histo[kPedMax];
 974   memset(histo,0,kPedMax*sizeof(UShort_t));
 975   for (Int_t i=0; i<fMaxTime; i++){
 976     if (signal[i]<=0) continue;
 977     if (signal[i]>max && i>firstBin) {
 978       max = signal[i];
 979       maxPos = i;
 980     }
 981     if (signal[i]>kPedMax-1) continue;
 982     histo[int(signal[i]+0.5)]++;
 983     count0++;
 984   }
 985   //
 986   for (Int_t i=1; i<kPedMax; i++){
 987     if (count1<count0*0.5) median=i;
 988     count1+=histo[i];
 989   }
 990   // truncated mean
 991   //
 992   Float_t count10=histo[median] ,mean=histo[median]*median,  rms=histo[median]*median*median ;
 993   Float_t count06=histo[median] ,mean06=histo[median]*median,  rms06=histo[median]*median*median ;
 994   Float_t count09=histo[median] ,mean09=histo[median]*median,  rms09=histo[median]*median*median ;
 995   //
 996   for (Int_t idelta=1; idelta<10; idelta++){
 997     if (median-idelta<=0) continue;
 998     if (median+idelta>kPedMax) continue;
 999     if (count06<0.6*count1){
1000       count06+=histo[median-idelta];
1001       mean06 +=histo[median-idelta]*(median-idelta);
1002       rms06  +=histo[median-idelta]*(median-idelta)*(median-idelta);
1003       count06+=histo[median+idelta];
1004       mean06 +=histo[median+idelta]*(median+idelta);
1005       rms06  +=histo[median+idelta]*(median+idelta)*(median+idelta);
1006     }
1007     if (count09<0.9*count1){
1008       count09+=histo[median-idelta];
1009       mean09 +=histo[median-idelta]*(median-idelta);
1010       rms09  +=histo[median-idelta]*(median-idelta)*(median-idelta);
1011       count09+=histo[median+idelta];
1012       mean09 +=histo[median+idelta]*(median+idelta);
1013       rms09  +=histo[median+idelta]*(median+idelta)*(median+idelta);
1014     }
1015     if (count10<0.95*count1){
1016       count10+=histo[median-idelta];
1017       mean +=histo[median-idelta]*(median-idelta);
1018       rms  +=histo[median-idelta]*(median-idelta)*(median-idelta);
1019       count10+=histo[median+idelta];
1020       mean +=histo[median+idelta]*(median+idelta);
1021       rms  +=histo[median+idelta]*(median+idelta)*(median+idelta);
1022     }
1023   }
1024   mean  /=count10;
1025   mean06/=count06;
1026   mean09/=count09;
1027   rms    = TMath::Sqrt(TMath::Abs(rms/count10-mean*mean));
1028   rms06    = TMath::Sqrt(TMath::Abs(rms06/count06-mean06*mean06));
1029  rms09    = TMath::Sqrt(TMath::Abs(rms09/count09-mean09*mean09));
1030   rmsEvent = rms09;
1031   //
1032   pedestalEvent = median;
1033   if (AliLog::GetDebugLevel("","AliTPCclustererMI")==0) return median;
1034   //
1035   UInt_t uid[3] = {UInt_t(id[0]),UInt_t(id[1]),UInt_t(id[2])};
1036   //
1037   // Dump mean signal info
1038   //
1039   (*fDebugStreamer)<<"Signal"<<
1040     "TimeStamp="<<fTimeStamp<<
1041     "EventType="<<fEventType<<
1042     "Sector="<<uid[0]<<
1043     "Row="<<uid[1]<<
1044     "Pad="<<uid[2]<<
1045     "Max="<<max<<
1046     "MaxPos="<<maxPos<<
1047     //
1048     "Median="<<median<<
1049     "Mean="<<mean<<
1050     "RMS="<<rms<<
1051     "Mean06="<<mean06<<
1052     "RMS06="<<rms06<<
1053     "Mean09="<<mean09<<
1054     "RMS09="<<rms09<<
1055     "RMSCalib="<<rmsCalib<<
1056     "PedCalib="<<pedestalCalib<<
1057     "\n";
1058   //
1059   // fill pedestal histogram
1060   //
1061   AliTPCROC * roc = AliTPCROC::Instance();
1062   if (!fAmplitudeHisto){
1063     fAmplitudeHisto = new TObjArray(72);
1064   }
1065   //
1066   if (uid[0]<roc->GetNSectors()
1067       && uid[1]< roc->GetNRows(uid[0])  &&
1068       uid[2] <roc->GetNPads(uid[0], uid[1])){
1069     TObjArray  * sectorArray = (TObjArray*)fAmplitudeHisto->UncheckedAt(uid[0]);
1070     if (!sectorArray){
1071       Int_t npads =roc->GetNChannels(uid[0]);
1072       sectorArray = new TObjArray(npads);
1073       fAmplitudeHisto->AddAt(sectorArray, uid[0]);
1074     }
1075     Int_t position =  uid[2]+roc->GetRowIndexes(uid[0])[uid[1]];
1076     TH1F * histo = (TH1F*)sectorArray->UncheckedAt(position);
1077     if (!histo){
1078       char hname[100];
1079       sprintf(hname,"Amp_%d_%d_%d",uid[0],uid[1],uid[2]);
1080       TFile * backup = gFile;
1081       fDebugStreamer->GetFile()->cd();
1082       histo = new TH1F(hname, hname, 100, 5,100);
1083       //histo->SetDirectory(0);     // histogram not connected to directory -(File)
1084       sectorArray->AddAt(histo, position);
1085       if (backup) backup->cd();
1086     }
1087     for (Int_t i=0; i<nchannels; i++){
1088       histo->Fill(signal[i]);
1089     }
1090   }
1091   //
1092   //
1093   //
1094   Float_t kMin =fRecoParam->GetDumpAmplitudeMin();   // minimal signal to be dumped
1095   Float_t *dsignal = new Float_t[nchannels];
1096   Float_t *dtime   = new Float_t[nchannels];
1097   for (Int_t i=0; i<nchannels; i++){
1098     dtime[i] = i;
1099     dsignal[i] = signal[i];
1100   }
1101   //
1102   // Digital noise
1103   //
1104  //  if (max-median>30.*TMath::Max(1.,Double_t(rms06)) &&  (((*fDebugStreamer)<<"SignalDN").GetSize()<kMaxDebugSize)){
1105 //     //
1106 //     //
1107 //     TGraph * graph =new TGraph(nchannels, dtime, dsignal);
1108 //     //
1109 //     //
1110 //     // jumps left - right
1111 //     Int_t    njumps0=0;
1112 //     Double_t deltaT0[2000];
1113 //     Double_t deltaA0[2000];
1114 //     Int_t    lastJump0 = fRecoParam->GetFirstBin();
1115 //     Int_t    njumps1=0;
1116 //     Double_t deltaT1[2000];
1117 //     Double_t deltaA1[2000];
1118 //     Int_t    lastJump1 = fRecoParam->GetFirstBin();
1119 //     Int_t    njumps2=0;
1120 //     Double_t deltaT2[2000];
1121 //     Double_t deltaA2[2000];
1122 //     Int_t    lastJump2 = fRecoParam->GetFirstBin();
1123
1124 //     for (Int_t itime=fRecoParam->GetFirstBin()+1; itime<fRecoParam->GetLastBin()-1; itime++){
1125 //       if (TMath::Abs(dsignal[itime]-dsignal[itime-1])>30.*TMath::Max(1.,Double_t(rms06))  &&
1126 //        TMath::Abs(dsignal[itime]-dsignal[itime+1])>30.*TMath::Max(1.,Double_t(rms06))  &&
1127 //        (dsignal[itime-1]-median<5.*rms06) &&
1128 //        (dsignal[itime+1]-median<5.*rms06)
1129 //        ){
1130 //      deltaA0[njumps0] = dsignal[itime]-dsignal[itime-1];
1131 //      deltaT0[njumps0] = itime-lastJump0;
1132 //      lastJump0 = itime;
1133 //      njumps0++;
1134 //       }
1135 //       if (TMath::Abs(dsignal[itime]-dsignal[itime-1])>30.*TMath::Max(1.,Double_t(rms06)) &&
1136 //        (dsignal[itime-1]-median<5.*rms06)
1137 //        ) {
1138 //      deltaA1[njumps1] = dsignal[itime]-dsignal[itime-1];
1139 //      deltaT1[njumps1] = itime-lastJump1;
1140 //      lastJump1 = itime;
1141 //      njumps1++;
1142 //       }
1143 //       if (TMath::Abs(dsignal[itime]-dsignal[itime+1])>30.*TMath::Max(1.,Double_t(rms06)) &&
1144 //        (dsignal[itime+1]-median<5.*rms06)
1145 //        ) {
1146 //      deltaA2[njumps2] = dsignal[itime]-dsignal[itime+1];
1147 //      deltaT2[njumps2] = itime-lastJump2;
1148 //      lastJump2 = itime;
1149 //      njumps2++;
1150 //       }
1151 //     }
1152 //     //
1153 //     if (njumps0>0 || njumps1>0 || njumps2>0){
1154 //       TGraph *graphDN0 = new TGraph(njumps0, deltaT0, deltaA0);
1155 //       TGraph *graphDN1 = new TGraph(njumps1, deltaT1, deltaA1);
1156 //       TGraph *graphDN2 = new TGraph(njumps2, deltaT2, deltaA2);
1157 //       (*fDebugStreamer)<<"SignalDN"<<    //digital - noise pads - or random sample of pads
1158 //      "TimeStamp="<<fTimeStamp<<
1159 //      "EventType="<<fEventType<<
1160 //      "Sector="<<uid[0]<<
1161 //      "Row="<<uid[1]<<
1162 //      "Pad="<<uid[2]<<
1163 //      "Graph="<<graph<<
1164 //      "Max="<<max<<
1165 //      "MaxPos="<<maxPos<<
1166 //      "Graph.="<<graph<<
1167 //      "P0GraphDN0.="<<graphDN0<<
1168 //      "P1GraphDN1.="<<graphDN1<<
1169 //      "P2GraphDN2.="<<graphDN2<<
1170 //      //
1171 //      "Median="<<median<<
1172 //      "Mean="<<mean<<
1173 //      "RMS="<<rms<<
1174 //      "Mean06="<<mean06<<
1175 //      "RMS06="<<rms06<<
1176 //      "Mean09="<<mean09<<
1177 //      "RMS09="<<rms09<<
1178 //      "\n";
1179 //       delete graphDN0;
1180 //       delete graphDN1;
1181 //       delete graphDN2;
1182 //     }
1183 //     delete graph;
1184 //   }
1185
1186   //
1187   // NOISE STUDY  Fourier transform
1188   //
1189   TGraph * graph;
1190   Bool_t random = (gRandom->Rndm()<0.0003);
1191   if (((*fDebugStreamer)<<"SignalN").GetSize()<kMaxDebugSize)
1192     if (max-median>kMin || rms06>1.*fParam->GetZeroSup() || random){
1193     graph =new TGraph(nchannels, dtime, dsignal);
1194     if (rms06>1.*fParam->GetZeroSup() || random){
1195       //Double_t *input, Double_t threshold, Bool_t locMax, Double_t *freq, Double_t *re, Double_t *im, Double_t *mag, Double_t *phi);
1196       Float_t * input = &(dsignal[fRecoParam->GetFirstBin()]);
1197       Float_t freq[2000], re[2000], im[2000], mag[2000], phi[2000];
1198       Int_t npoints = TransformFFT(input, -1,kFALSE, freq, re, im, mag, phi);
1199       TGraph *graphMag0 = new TGraph(npoints, freq, mag);
1200       TGraph *graphPhi0 = new TGraph(npoints, freq, phi);
1201       npoints = TransformFFT(input, 0.5,kTRUE, freq, re, im, mag, phi);
1202       TGraph *graphMag1 = new TGraph(npoints, freq, mag);
1203       TGraph *graphPhi1 = new TGraph(npoints, freq, phi);
1204
1205       (*fDebugStreamer)<<"SignalN"<<    //noise pads - or random sample of pads
1206         "TimeStamp="<<fTimeStamp<<
1207         "EventType="<<fEventType<<
1208         "Sector="<<uid[0]<<
1209         "Row="<<uid[1]<<
1210         "Pad="<<uid[2]<<
1211         "Graph.="<<graph<<
1212         "Max="<<max<<
1213         "MaxPos="<<maxPos<<
1214         //
1215         "Median="<<median<<
1216         "Mean="<<mean<<
1217         "RMS="<<rms<<
1218         "Mean06="<<mean06<<
1219         "RMS06="<<rms06<<
1220         "Mean09="<<mean09<<
1221         "RMS09="<<rms09<<
1222         // FFT part
1223         "Mag0.="<<graphMag0<<
1224         "Mag1.="<<graphMag1<<
1225         "Phi0.="<<graphPhi0<<
1226         "Phi1.="<<graphPhi1<<
1227         "\n";
1228       delete graphMag0;
1229       delete graphMag1;
1230       delete graphPhi0;
1231       delete graphPhi1;
1232     }
1233     //
1234     // Big signals dumping
1235     //
1236
1237     if (max-median>kMin &&maxPos>AliTPCReconstructor::GetRecoParam()->GetFirstBin())
1238       (*fDebugStreamer)<<"SignalB"<<     // pads with signal
1239         "TimeStamp="<<fTimeStamp<<
1240         "EventType="<<fEventType<<
1241         "Sector="<<uid[0]<<
1242         "Row="<<uid[1]<<
1243         "Pad="<<uid[2]<<
1244         "Graph="<<graph<<
1245         "Max="<<max<<
1246         "MaxPos="<<maxPos<<
1247         //
1248         "Median="<<median<<
1249         "Mean="<<mean<<
1250         "RMS="<<rms<<
1251         "Mean06="<<mean06<<
1252         "RMS06="<<rms06<<
1253         "Mean09="<<mean09<<
1254         "RMS09="<<rms09<<
1255         "\n";
1256     delete graph;
1257   }
1258
1259   //
1260   //
1261   //  Central Electrode signal analysis
1262   //
1263   Float_t ceQmax  =0, ceQsum=0, ceTime=0;
1264   Float_t cemean  = mean06, cerms=rms06 ;
1265   Int_t    cemaxpos= 0;
1266   Float_t ceThreshold=5.*cerms;
1267   Float_t ceSumThreshold=8.*cerms;
1268   const Int_t    kCemin=5;  // range for the analysis of the ce signal +- channels from the peak
1269   const Int_t    kCemax=5;
1270   for (Int_t i=nchannels-2; i>nchannels/2; i--){
1271     if ( (dsignal[i]-mean06)>ceThreshold && dsignal[i]>=dsignal[i+1] && dsignal[i]>=dsignal[i-1] ){
1272       cemaxpos=i;
1273       break;
1274     }
1275   }
1276   if (cemaxpos!=0){
1277     ceQmax = 0;
1278     Int_t cemaxpos2=0;
1279     for (Int_t i=cemaxpos-20; i<cemaxpos+5; i++){
1280       if (i<0 || i>nchannels-1) continue;
1281       Double_t val=dsignal[i]- cemean;
1282       if (val>ceQmax){
1283         cemaxpos2=i;
1284         ceQmax = val;
1285       }
1286     }
1287     cemaxpos = cemaxpos2;
1288
1289     for (Int_t i=cemaxpos-kCemin; i<cemaxpos+kCemax; i++){
1290       if (i>0 && i<nchannels&&dsignal[i]- cemean>0){
1291         Double_t val=dsignal[i]- cemean;
1292         ceTime+=val*dtime[i];
1293         ceQsum+=val;
1294         if (val>ceQmax) ceQmax=val;
1295       }
1296     }
1297     if (ceQmax&&ceQsum>ceSumThreshold) {
1298       ceTime/=ceQsum;
1299       (*fDebugStreamer)<<"Signalce"<<
1300         "TimeStamp="<<fTimeStamp<<
1301         "EventType="<<fEventType<<
1302         "Sector="<<uid[0]<<
1303         "Row="<<uid[1]<<
1304         "Pad="<<uid[2]<<
1305         "Max="<<ceQmax<<
1306         "Qsum="<<ceQsum<<
1307         "Time="<<ceTime<<
1308         "RMS06="<<rms06<<
1309         //
1310         "\n";
1311     }
1312   }
1313   // end of ce signal analysis
1314   //
1315
1316   //
1317   //  Gating grid signal analysis
1318   //
1319   Double_t ggQmax  =0, ggQsum=0, ggTime=0;
1320   Double_t ggmean  = mean06, ggrms=rms06 ;
1321   Int_t    ggmaxpos= 0;
1322   Double_t ggThreshold=5.*ggrms;
1323   Double_t ggSumThreshold=8.*ggrms;
1324
1325   for (Int_t i=1; i<nchannels/4; i++){
1326     if ( (dsignal[i]-mean06)>ggThreshold && dsignal[i]>=dsignal[i+1] && dsignal[i]>=dsignal[i-1] &&
1327          (dsignal[i]+dsignal[i+1]+dsignal[i-1]-3*mean06)>ggSumThreshold){
1328       ggmaxpos=i;
1329       if (dsignal[i-1]>dsignal[i+1]) ggmaxpos=i-1;
1330       break;
1331     }
1332   }
1333   if (ggmaxpos!=0){
1334       for (Int_t i=ggmaxpos-1; i<ggmaxpos+3; i++){
1335           if (i>0 && i<nchannels && dsignal[i]-ggmean>0){
1336               Double_t val=dsignal[i]- ggmean;
1337               ggTime+=val*dtime[i];
1338               ggQsum+=val;
1339               if (val>ggQmax) ggQmax=val;
1340           }
1341       }
1342       if (ggQmax&&ggQsum>ggSumThreshold) {
1343           ggTime/=ggQsum;
1344           (*fDebugStreamer)<<"Signalgg"<<
1345             "TimeStamp="<<fTimeStamp<<
1346             "EventType="<<fEventType<<
1347               "Sector="<<uid[0]<<
1348               "Row="<<uid[1]<<
1349               "Pad="<<uid[2]<<
1350               "Max="<<ggQmax<<
1351               "Qsum="<<ggQsum<<
1352               "Time="<<ggTime<<
1353               "RMS06="<<rms06<<
1354               //
1355               "\n";
1356       }
1357   }
1358   // end of gg signal analysis
1359
1360
1361   delete [] dsignal;
1362   delete [] dtime;
1363   if (rms06>fRecoParam->GetMaxNoise()) {
1364     pedestalEvent+=1024.;
1365     return 1024+median; // sign noisy channel in debug mode
1366   }
1367   return median;
1368 }
1369
1370
1371
1372 void AliTPCclustererMI::DumpHistos(){
1373   //
1374   // Dump histogram information
1375   //
1376   if (!fAmplitudeHisto) return;
1377   AliTPCROC * roc = AliTPCROC::Instance();
1378   for (UInt_t isector=0; isector<AliTPCROC::Instance()->GetNSectors(); isector++){
1379     TObjArray * array = (TObjArray*)fAmplitudeHisto->UncheckedAt(isector);
1380     if (!array) continue;
1381     for (UInt_t ipad = 0; ipad <(UInt_t)array->GetEntriesFast(); ipad++){
1382       TH1F * histo = (TH1F*) array->UncheckedAt(ipad);
1383       if (!histo) continue;
1384       if (histo->GetEntries()<100) continue;
1385       histo->Fit("gaus","q");
1386       Float_t mean =  histo->GetMean();
1387       Float_t rms  =  histo->GetRMS();
1388       Float_t gmean = histo->GetFunction("gaus")->GetParameter(1);
1389       Float_t gsigma = histo->GetFunction("gaus")->GetParameter(2);
1390       Float_t gmeanErr = histo->GetFunction("gaus")->GetParError(1);
1391       Float_t gsigmaErr = histo->GetFunction("gaus")->GetParError(2);
1392       Float_t max = histo->GetFunction("gaus")->GetParameter(0);
1393
1394       // get pad number
1395       UInt_t row=0, pad =0;
1396       const UInt_t *indexes =roc->GetRowIndexes(isector);
1397       for (UInt_t irow=0; irow<roc->GetNRows(isector); irow++){
1398         if (indexes[irow]<=ipad){
1399           row = irow;
1400           pad = ipad-indexes[irow];
1401         }
1402       }
1403       Int_t rpad = pad - (AliTPCROC::Instance()->GetNPads(isector,row))/2;
1404       //
1405       (*fDebugStreamer)<<"Fit"<<
1406         "TimeStamp="<<fTimeStamp<<
1407         "EventType="<<fEventType<<
1408         "Sector="<<isector<<
1409         "Row="<<row<<
1410         "Pad="<<pad<<
1411         "RPad="<<rpad<<
1412         "Max="<<max<<
1413         "Mean="<<mean<<
1414         "RMS="<<rms<<
1415         "GMean="<<gmean<<
1416         "GSigma="<<gsigma<<
1417         "GMeanErr="<<gmeanErr<<
1418         "GSigmaErr="<<gsigmaErr<<
1419         "\n";
1420       if (array->UncheckedAt(ipad)) fDebugStreamer->StoreObject(array->UncheckedAt(ipad));
1421     }
1422   }
1423 }
1424
1425
1426
1427 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)
1428 {
1429   //
1430   // calculate fourrie transform
1431   // return only frequncies with mag over threshold
1432   // if locMax is spectified only freque with local maxima over theshold is returned
1433
1434   if (! fFFTr2c) return kFALSE;
1435   if (!freq) return kFALSE;
1436
1437   Int_t current=0;
1438   Int_t nPoints = fRecoParam->GetLastBin()-fRecoParam->GetFirstBin();
1439   Double_t *in = new Double_t[nPoints];
1440   Double_t *rfft = new Double_t[nPoints];
1441   Double_t *ifft = new Double_t[nPoints];
1442   for (Int_t i=0; i<nPoints; i++){in[i]=input[i];}
1443   fFFTr2c->SetPoints(in);
1444   fFFTr2c->Transform();
1445   fFFTr2c->GetPointsComplex(rfft, ifft);
1446   for (Int_t i=3; i<nPoints/2-3; i++){
1447     Float_t lmag =  TMath::Sqrt(rfft[i]*rfft[i]+ifft[i]*ifft[i])/nPoints;
1448     if (lmag<threshold) continue;
1449     if (locMax){
1450       if ( TMath::Sqrt(rfft[i-1]*rfft[i-1]+ifft[i-1]*ifft[i-1])/nPoints>lmag) continue;
1451       if ( TMath::Sqrt(rfft[i+1]*rfft[i+1]+ifft[i+1]*ifft[i+1])/nPoints>lmag) continue;
1452       if ( TMath::Sqrt(rfft[i-2]*rfft[i-2]+ifft[i-2]*ifft[i-2])/nPoints>lmag) continue;
1453       if ( TMath::Sqrt(rfft[i+2]*rfft[i+2]+ifft[i+2]*ifft[i+2])/nPoints>lmag) continue;
1454       if ( TMath::Sqrt(rfft[i-3]*rfft[i-3]+ifft[i-3]*ifft[i-3])/nPoints>lmag) continue;
1455       if ( TMath::Sqrt(rfft[i+3]*rfft[i+3]+ifft[i+3]*ifft[i+3])/nPoints>lmag) continue;
1456     }
1457
1458     freq[current] = Float_t(i)/Float_t(nPoints);
1459     //
1460     re[current] = rfft[i];
1461     im[current] = ifft[i];
1462     mag[current]=lmag;
1463     phi[current]=TMath::ATan2(ifft[i],rfft[i]);
1464     current++;
1465   }
1466   delete [] in;
1467   delete [] rfft;
1468   delete [] ifft;
1469   return current;
1470 }
1471