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