TPC/AliTPCclustererMI.cxx

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