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