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