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