TPC/AliTPCclustererMI.cxx

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