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