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