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