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