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