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