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