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