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