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