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