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