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