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