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