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()>5) {
 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()>5) {
 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     AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ;
 914     transform->SetCurrentTimeStamp(fTimeStamp);
 915     transform->SetCurrentRun(rawReader->GetRunNumber());
 916   }
 917
 918   // creaate one TClonesArray for all clusters
 919   if(fBClonesArray && !fOutputClonesArray) fOutputClonesArray = new TClonesArray("AliTPCclusterMI",1000);
 920   // reset counter
 921   fNclusters  = 0;
 922
 923   fMaxTime = fRecoParam->GetLastBin() + 6; // add 3 virtual time bins before and 3 after
 924 //   const Int_t kNIS = fParam->GetNInnerSector();
 925 //   const Int_t kNOS = fParam->GetNOuterSector();
 926 //   const Int_t kNS = kNIS + kNOS;
 927   fZWidth = fParam->GetZWidth();
 928   Int_t zeroSup = fParam->GetZeroSup();
 929   //
 930   //alocate memory for sector - maximal case
 931   //
 932   Float_t** allBins = NULL;
 933   Int_t** allSigBins = NULL;
 934   Int_t*  allNSigBins = NULL;
 935   Int_t nRowsMax = roc->GetNRows(roc->GetNSector()-1);
 936   Int_t nPadsMax = roc->GetNPads(roc->GetNSector()-1,nRowsMax-1);
 937   allBins = new Float_t*[nRowsMax];
 938   allSigBins = new Int_t*[nRowsMax];
 939   allNSigBins = new Int_t[nRowsMax];
 940   for (Int_t iRow = 0; iRow < nRowsMax; iRow++) {
 941     //
 942     Int_t maxBin = fMaxTime*(nPadsMax+6);  // add 3 virtual pads  before and 3 after
 943     allBins[iRow] = new Float_t[maxBin];
 944     memset(allBins[iRow],0,sizeof(Float_t)*maxBin);
 945     allSigBins[iRow] = new Int_t[maxBin];
 946     allNSigBins[iRow]=0;
 947   }
 948
 949   Int_t prevSector=-1;
 950   rawReader->Reset();
 951   Int_t digCounter=0;
 952   //
 953   // Loop over DDLs
 954   //
 955   const Int_t kNIS = fParam->GetNInnerSector();
 956   const Int_t kNOS = fParam->GetNOuterSector();
 957   const Int_t kNS = kNIS + kNOS;
 958
 959   for(fSector = 0; fSector < kNS; fSector++) {
 960
 961     Int_t nRows = 0;
 962     Int_t nDDLs = 0, indexDDL = 0;
 963     if (fSector < kNIS) {
 964       nRows = fParam->GetNRowLow();
 965       fSign = (fSector < kNIS/2) ? 1 : -1;
 966       nDDLs = 2;
 967       indexDDL = fSector * 2;
 968     }
 969     else {
 970       nRows = fParam->GetNRowUp();
 971       fSign = ((fSector-kNIS) < kNOS/2) ? 1 : -1;
 972       nDDLs = 4;
 973       indexDDL = (fSector-kNIS) * 4 + kNIS * 2;
 974     }
 975
 976     // load the raw data for corresponding DDLs
 977     rawReader->Reset();
 978     rawReader->Select("TPC",indexDDL,indexDDL+nDDLs-1);
 979
 980   while (input.NextDDL()){
 981     if (input.GetSector() != fSector)
 982       AliFatal(Form("Sector index mismatch ! Expected (%d), but got (%d) !",fSector,input.GetSector()));
 983
 984     //Int_t nRows = fParam->GetNRow(fSector);
 985
 986     AliTPCCalROC * gainROC    = gainTPC->GetCalROC(fSector);  // pad gains per given sector
 987     // Begin loop over altro data
 988     Bool_t calcPedestal = fRecoParam->GetCalcPedestal();
 989     Float_t gain =1;
 990
 991     //loop over pads
 992     while ( input.NextChannel() ) {
 993       Int_t iRow = input.GetRow();
 994       if (iRow < 0){
 995         continue;
 996       }
 997       if (iRow >= nRows){
 998         AliError(Form("Pad-row index (%d) outside the range (%d -> %d) !",
 999                       iRow, 0, nRows -1));
1000         continue;
1001       }
1002       //pad
1003       Int_t iPad = input.GetPad();
1004       if (iPad < 0 || iPad >= nPadsMax) {
1005         AliError(Form("Pad index (%d) outside the range (%d -> %d) !",
1006                       iPad, 0, nPadsMax-1));
1007         continue;
1008       }
1009       gain    = gainROC->GetValue(iRow,iPad);
1010       iPad+=3;
1011
1012       //loop over bunches
1013       while ( input.NextBunch() ){
1014         Int_t  startTbin    = (Int_t)input.GetStartTimeBin();
1015         Int_t  bunchlength  = (Int_t)input.GetBunchLength();
1016         const UShort_t *sig = input.GetSignals();
1017         for (Int_t iTime = 0; iTime<bunchlength; iTime++){
1018           Int_t iTimeBin=startTbin-iTime;
1019           if ( iTimeBin < fRecoParam->GetFirstBin() || iTimeBin >= fRecoParam->GetLastBin()){
1020             continue;
1021             AliFatal(Form("Timebin index (%d) outside the range (%d -> %d) !",
1022                           iTimeBin, 0, iTimeBin -1));
1023           }
1024           iTimeBin+=3;
1025           //signal
1026           Float_t signal=(Float_t)sig[iTime];
1027           if (!calcPedestal && signal <= zeroSup) continue;
1028
1029           if (!calcPedestal) {
1030             Int_t bin = iPad*fMaxTime+iTimeBin;
1031             if (gain>0){
1032               allBins[iRow][bin] = signal/gain;
1033             }else{
1034               allBins[iRow][bin] =0;
1035             }
1036             allSigBins[iRow][allNSigBins[iRow]++] = bin;
1037           }else{
1038             allBins[iRow][iPad*fMaxTime+iTimeBin] = signal;
1039           }
1040           allBins[iRow][iPad*fMaxTime+0]+=1.;  // pad with signal
1041
1042           // Temporary
1043           digCounter++;
1044         }// end loop signals in bunch
1045       }// end loop bunches
1046     } // end loop pads
1047     //
1048     //
1049     //
1050     //
1051     // Now loop over rows and perform pedestal subtraction
1052     if (digCounter==0) continue;
1053   } // End of loop over sectors
1054   //process last sector
1055   if ( digCounter>0 ){
1056     ProcessSectorData(allBins, allSigBins, allNSigBins);
1057     for (Int_t iRow = 0; iRow < fParam->GetNRow(fSector); iRow++) {
1058       Int_t maxPad = fParam->GetNPads(fSector,iRow);
1059       Int_t maxBin = fMaxTime*(maxPad+6);  // add 3 virtual pads  before and 3 after
1060       memset(allBins[iRow],0,sizeof(Float_t)*maxBin);
1061       allNSigBins[iRow] = 0;
1062     }
1063     prevSector=fSector;
1064     digCounter=0;
1065   }
1066   }
1067
1068
1069
1070   for (Int_t iRow = 0; iRow < nRowsMax; iRow++) {
1071     delete [] allBins[iRow];
1072     delete [] allSigBins[iRow];
1073   }
1074   delete [] allBins;
1075   delete [] allSigBins;
1076   delete [] allNSigBins;
1077
1078   if (rawReader->GetEventId() && fOutput ){
1079     Info("Digits2Clusters", "File  %s Event\t%d\tNumber of found clusters : %d\n", fOutput->GetName(),*(rawReader->GetEventId()), fNclusters);
1080   }
1081
1082   if(rawReader->GetEventId()) {
1083     Info("Digits2Clusters", "Event\t%d\tNumber of found clusters : %d\n",*(rawReader->GetEventId()), fNclusters);
1084   }
1085
1086   if(fBClonesArray) {
1087     //Info("Digits2Clusters", "Number of found clusters : %d\n",fOutputClonesArray->GetEntriesFast());
1088   }
1089 }
1090
1091
1092
1093
1094
1095 void AliTPCclustererMI::Digits2ClustersOld
1096 (AliRawReader* rawReader)
1097 {
1098 //-----------------------------------------------------------------
1099 // This is a cluster finder for the TPC raw data.
1100 // The method assumes NO ordering of the altro channels.
1101 // The pedestal subtraction can be switched on and off
1102 // using an option of the TPC reconstructor
1103 //-----------------------------------------------------------------
1104   fRecoParam = AliTPCReconstructor::GetRecoParam();
1105   if (!fRecoParam){
1106     AliFatal("Can not get the reconstruction parameters");
1107   }
1108   if(AliTPCReconstructor::StreamLevel()>5) {
1109     AliInfo("Parameter Dumps");
1110     fParam->Dump();
1111     fRecoParam->Dump();
1112   }
1113   fRowDig = NULL;
1114   AliTPCROC * roc = AliTPCROC::Instance();
1115   AliTPCCalPad * gainTPC = AliTPCcalibDB::Instance()->GetPadGainFactor();
1116   AliTPCAltroMapping** mapping =AliTPCcalibDB::Instance()->GetMapping();
1117   //
1118   AliTPCRawStream input(rawReader,(AliAltroMapping**)mapping);
1119   fEventHeader = (AliRawEventHeaderBase*)rawReader->GetEventHeader();
1120   if (fEventHeader){
1121     fTimeStamp = fEventHeader->Get("Timestamp");
1122     fEventType = fEventHeader->Get("Type");
1123   }
1124
1125   // creaate one TClonesArray for all clusters
1126   if(fBClonesArray && !fOutputClonesArray) fOutputClonesArray = new TClonesArray("AliTPCclusterMI",1000);
1127   // reset counter
1128   fNclusters  = 0;
1129
1130   fMaxTime = fRecoParam->GetLastBin() + 6; // add 3 virtual time bins before and 3 after
1131   const Int_t kNIS = fParam->GetNInnerSector();
1132   const Int_t kNOS = fParam->GetNOuterSector();
1133   const Int_t kNS = kNIS + kNOS;
1134   fZWidth = fParam->GetZWidth();
1135   Int_t zeroSup = fParam->GetZeroSup();
1136   //
1137   //alocate memory for sector - maximal case
1138   //
1139   Float_t** allBins = NULL;
1140   Int_t** allSigBins = NULL;
1141   Int_t*  allNSigBins = NULL;
1142   Int_t nRowsMax = roc->GetNRows(roc->GetNSector()-1);
1143   Int_t nPadsMax = roc->GetNPads(roc->GetNSector()-1,nRowsMax-1);
1144   allBins = new Float_t*[nRowsMax];
1145   allSigBins = new Int_t*[nRowsMax];
1146   allNSigBins = new Int_t[nRowsMax];
1147   for (Int_t iRow = 0; iRow < nRowsMax; iRow++) {
1148     //
1149     Int_t maxBin = fMaxTime*(nPadsMax+6);  // add 3 virtual pads  before and 3 after
1150     allBins[iRow] = new Float_t[maxBin];
1151     memset(allBins[iRow],0,sizeof(Float_t)*maxBin);
1152     allSigBins[iRow] = new Int_t[maxBin];
1153     allNSigBins[iRow]=0;
1154   }
1155   //
1156   // Loop over sectors
1157   //
1158   for(fSector = 0; fSector < kNS; fSector++) {
1159
1160     Int_t nRows = 0;
1161     Int_t nDDLs = 0, indexDDL = 0;
1162     if (fSector < kNIS) {
1163       nRows = fParam->GetNRowLow();
1164       fSign = (fSector < kNIS/2) ? 1 : -1;
1165       nDDLs = 2;
1166       indexDDL = fSector * 2;
1167     }
1168     else {
1169       nRows = fParam->GetNRowUp();
1170       fSign = ((fSector-kNIS) < kNOS/2) ? 1 : -1;
1171       nDDLs = 4;
1172       indexDDL = (fSector-kNIS) * 4 + kNIS * 2;
1173     }
1174
1175     // load the raw data for corresponding DDLs
1176     rawReader->Reset();
1177     rawReader->Select("TPC",indexDDL,indexDDL+nDDLs-1);
1178
1179     // select only good sector
1180     input.Next();
1181     if(input.GetSector() != fSector) continue;
1182
1183     AliTPCCalROC * gainROC    = gainTPC->GetCalROC(fSector);  // pad gains per given sector
1184
1185     for (Int_t iRow = 0; iRow < nRows; iRow++) {
1186       Int_t maxPad;
1187       if (fSector < kNIS)
1188         maxPad = fParam->GetNPadsLow(iRow);
1189       else
1190         maxPad = fParam->GetNPadsUp(iRow);
1191
1192       Int_t maxBin = fMaxTime*(maxPad+6);  // add 3 virtual pads  before and 3 after
1193       memset(allBins[iRow],0,sizeof(Float_t)*maxBin);
1194       allNSigBins[iRow] = 0;
1195     }
1196
1197     Int_t digCounter=0;
1198     // Begin loop over altro data
1199     Bool_t calcPedestal = fRecoParam->GetCalcPedestal();
1200     Float_t gain =1;
1201     Int_t lastPad=-1;
1202
1203     input.Reset();
1204     while (input.Next()) {
1205       if (input.GetSector() != fSector)
1206         AliFatal(Form("Sector index mismatch ! Expected (%d), but got (%d) !",fSector,input.GetSector()));
1207
1208
1209       Int_t iRow = input.GetRow();
1210       if (iRow < 0){
1211         continue;
1212       }
1213
1214       if (iRow < 0 || iRow >= nRows){
1215         AliError(Form("Pad-row index (%d) outside the range (%d -> %d) !",
1216                       iRow, 0, nRows -1));
1217         continue;
1218       }
1219       //pad
1220       Int_t iPad = input.GetPad();
1221       if (iPad < 0 || iPad >= nPadsMax) {
1222         AliError(Form("Pad index (%d) outside the range (%d -> %d) !",
1223                       iPad, 0, nPadsMax-1));
1224         continue;
1225       }
1226       if (iPad!=lastPad){
1227         gain    = gainROC->GetValue(iRow,iPad);
1228         lastPad = iPad;
1229       }
1230       iPad+=3;
1231       //time
1232       Int_t iTimeBin = input.GetTime();
1233       if ( iTimeBin < fRecoParam->GetFirstBin() || iTimeBin >= fRecoParam->GetLastBin()){
1234         continue;
1235         AliFatal(Form("Timebin index (%d) outside the range (%d -> %d) !",
1236                       iTimeBin, 0, iTimeBin -1));
1237       }
1238       iTimeBin+=3;
1239
1240       //signal
1241       Float_t signal = input.GetSignal();
1242       if (!calcPedestal && signal <= zeroSup) continue;
1243
1244       if (!calcPedestal) {
1245         Int_t bin = iPad*fMaxTime+iTimeBin;
1246         if (gain>0){
1247           allBins[iRow][bin] = signal/gain;
1248         }else{
1249           allBins[iRow][bin] =0;
1250         }
1251         allSigBins[iRow][allNSigBins[iRow]++] = bin;
1252       }else{
1253         allBins[iRow][iPad*fMaxTime+iTimeBin] = signal;
1254       }
1255       allBins[iRow][iPad*fMaxTime+0]+=1.;  // pad with signal
1256
1257       // Temporary
1258       digCounter++;
1259     } // End of the loop over altro data
1260     //
1261     //
1262     //
1263     //
1264     // Now loop over rows and perform pedestal subtraction
1265     if (digCounter==0) continue;
1266     ProcessSectorData(allBins, allSigBins, allNSigBins);
1267   } // End of loop over sectors
1268
1269   for (Int_t iRow = 0; iRow < nRowsMax; iRow++) {
1270     delete [] allBins[iRow];
1271     delete [] allSigBins[iRow];
1272   }
1273   delete [] allBins;
1274   delete [] allSigBins;
1275   delete [] allNSigBins;
1276
1277   if (rawReader->GetEventId() && fOutput ){
1278     Info("Digits2Clusters", "File  %s Event\t%d\tNumber of found clusters : %d\n", fOutput->GetName(),*(rawReader->GetEventId()), fNclusters);
1279   }
1280
1281   if(rawReader->GetEventId()) {
1282     Info("Digits2Clusters", "Event\t%d\tNumber of found clusters : %d\n",*(rawReader->GetEventId()), fNclusters);
1283   }
1284
1285   if(fBClonesArray) {
1286     //Info("Digits2Clusters", "Number of found clusters : %d\n",fOutputClonesArray->GetEntriesFast());
1287   }
1288 }
1289
1290 void AliTPCclustererMI::FindClusters(AliTPCCalROC * noiseROC)
1291 {
1292
1293   //
1294   // add virtual charge at the edge
1295   //
1296   Double_t kMaxDumpSize = 500000;
1297   if (!fOutput) {
1298     fBDumpSignal =kFALSE;
1299   }else{
1300     if (fRecoParam->GetCalcPedestal() && fOutput->GetZipBytes()< kMaxDumpSize) fBDumpSignal =kTRUE;   //dump signal flag
1301   }
1302
1303   fNcluster=0;
1304   fLoop=1;
1305   Int_t crtime = Int_t((fParam->GetZLength(fSector)-fRecoParam->GetCtgRange()*fRx)/fZWidth-fParam->GetNTBinsL1()-5);
1306   Float_t minMaxCutAbs       = fRecoParam->GetMinMaxCutAbs();
1307   Float_t minLeftRightCutAbs = fRecoParam->GetMinLeftRightCutAbs();
1308   Float_t minUpDownCutAbs    = fRecoParam->GetMinUpDownCutAbs();
1309   Float_t minMaxCutSigma       = fRecoParam->GetMinMaxCutSigma();
1310   Float_t minLeftRightCutSigma = fRecoParam->GetMinLeftRightCutSigma();
1311   Float_t minUpDownCutSigma    = fRecoParam->GetMinUpDownCutSigma();
1312   for (Int_t iSig = 0; iSig < fNSigBins; iSig++) {
1313     Int_t i = fSigBins[iSig];
1314     if (i%fMaxTime<=crtime) continue;
1315     Float_t *b = &fBins[i];
1316     //absolute custs
1317     if (b[0]<minMaxCutAbs) continue;   //threshold for maxima
1318     //
1319     if (b[-1]+b[1]+b[-fMaxTime]+b[fMaxTime]<=0) continue;  // cut on isolated clusters
1320     if (b[-1]+b[1]<=0) continue;               // cut on isolated clusters
1321     if (b[-fMaxTime]+b[fMaxTime]<=0) continue; // cut on isolated clusters
1322     //
1323     if ((b[0]+b[-1]+b[1])<minUpDownCutAbs) continue;   //threshold for up down  (TRF)
1324     if ((b[0]+b[-fMaxTime]+b[fMaxTime])<minLeftRightCutAbs) continue;   //threshold for left right (PRF)
1325     if (!IsMaximum(*b,fMaxTime,b)) continue;
1326     //
1327     Float_t noise = noiseROC->GetValue(fRow, i/fMaxTime);
1328     if (noise>fRecoParam->GetMaxNoise()) continue;
1329     // sigma cuts
1330     if (b[0]<minMaxCutSigma*noise) continue;   //threshold form maxima
1331     if ((b[0]+b[-1]+b[1])<minUpDownCutSigma*noise) continue;   //threshold for up town TRF
1332     if ((b[0]+b[-fMaxTime]+b[fMaxTime])<minLeftRightCutSigma*noise) continue;   //threshold for left right (PRF)
1333
1334     AliTPCclusterMI c;   // default cosntruction  without info
1335     Int_t dummy=0;
1336     MakeCluster(i, fMaxTime, fBins, dummy,c);
1337
1338     //}
1339   }
1340 }
1341
1342 Bool_t AliTPCclustererMI::AcceptCluster(AliTPCclusterMI *cl){
1343   //
1344   // Currently hack to filter digital noise (15.06.2008)
1345   // To be parameterized in the AliTPCrecoParam
1346   // More inteligent way  to be used in future
1347   // Acces to the proper pedestal file needed
1348   //
1349   if (cl->GetMax()<400) return kTRUE;
1350   Double_t ratio = cl->GetQ()/cl->GetMax();
1351   if (cl->GetMax()>700){
1352     if ((ratio - int(ratio)>0.8)) return kFALSE;
1353   }
1354   if ((ratio - int(ratio)<0.95)) return kTRUE;
1355   return kFALSE;
1356 }
1357
1358
1359 Double_t AliTPCclustererMI::ProcesSignal(Float_t *signal, Int_t nchannels, Int_t id[3], Double_t &rmsEvent, Double_t &pedestalEvent){
1360   //
1361   // process signal on given pad - + streaming of additional information in special mode
1362   //
1363   // id[0] - sector
1364   // id[1] - row
1365   // id[2] - pad
1366
1367   //
1368   // ESTIMATE pedestal and the noise
1369   //
1370   const Int_t kPedMax = 100;
1371   Float_t  max    =  0;
1372   Float_t  maxPos =  0;
1373   Int_t    median =  -1;
1374   Int_t    count0 =  0;
1375   Int_t    count1 =  0;
1376   Float_t  rmsCalib   = rmsEvent;       // backup initial value ( from calib)
1377   Float_t  pedestalCalib = pedestalEvent;// backup initial value ( from calib)
1378   Int_t    firstBin = fRecoParam->GetFirstBin();
1379   //
1380   UShort_t histo[kPedMax];
1381   //memset(histo,0,kPedMax*sizeof(UShort_t));
1382   for (Int_t i=0; i<kPedMax; i++) histo[i]=0;
1383   for (Int_t i=0; i<fMaxTime; i++){
1384     if (signal[i]<=0) continue;
1385     if (signal[i]>max && i>firstBin) {
1386       max = signal[i];
1387       maxPos = i;
1388     }
1389     if (signal[i]>kPedMax-1) continue;
1390     histo[int(signal[i]+0.5)]++;
1391     count0++;
1392   }
1393   //
1394   for (Int_t i=1; i<kPedMax; i++){
1395     if (count1<count0*0.5) median=i;
1396     count1+=histo[i];
1397   }
1398   // truncated mean
1399   //
1400   Float_t count10=histo[median] ,mean=histo[median]*median,  rms=histo[median]*median*median ;
1401   Float_t count06=histo[median] ,mean06=histo[median]*median,  rms06=histo[median]*median*median ;
1402   Float_t count09=histo[median] ,mean09=histo[median]*median,  rms09=histo[median]*median*median ;
1403   //
1404   for (Int_t idelta=1; idelta<10; idelta++){
1405     if (median-idelta<=0) continue;
1406     if (median+idelta>kPedMax) continue;
1407     if (count06<0.6*count1){
1408       count06+=histo[median-idelta];
1409       mean06 +=histo[median-idelta]*(median-idelta);
1410       rms06  +=histo[median-idelta]*(median-idelta)*(median-idelta);
1411       count06+=histo[median+idelta];
1412       mean06 +=histo[median+idelta]*(median+idelta);
1413       rms06  +=histo[median+idelta]*(median+idelta)*(median+idelta);
1414     }
1415     if (count09<0.9*count1){
1416       count09+=histo[median-idelta];
1417       mean09 +=histo[median-idelta]*(median-idelta);
1418       rms09  +=histo[median-idelta]*(median-idelta)*(median-idelta);
1419       count09+=histo[median+idelta];
1420       mean09 +=histo[median+idelta]*(median+idelta);
1421       rms09  +=histo[median+idelta]*(median+idelta)*(median+idelta);
1422     }
1423     if (count10<0.95*count1){
1424       count10+=histo[median-idelta];
1425       mean +=histo[median-idelta]*(median-idelta);
1426       rms  +=histo[median-idelta]*(median-idelta)*(median-idelta);
1427       count10+=histo[median+idelta];
1428       mean +=histo[median+idelta]*(median+idelta);
1429       rms  +=histo[median+idelta]*(median+idelta)*(median+idelta);
1430     }
1431   }
1432   if (count10) {
1433     mean  /=count10;
1434     rms    = TMath::Sqrt(TMath::Abs(rms/count10-mean*mean));
1435   }
1436   if (count06) {
1437     mean06/=count06;
1438     rms06    = TMath::Sqrt(TMath::Abs(rms06/count06-mean06*mean06));
1439   }
1440   if (count09) {
1441     mean09/=count09;
1442     rms09    = TMath::Sqrt(TMath::Abs(rms09/count09-mean09*mean09));
1443   }
1444   rmsEvent = rms09;
1445   //
1446   pedestalEvent = median;
1447   if (AliLog::GetDebugLevel("","AliTPCclustererMI")==0) return median;
1448   //
1449   UInt_t uid[3] = {UInt_t(id[0]),UInt_t(id[1]),UInt_t(id[2])};
1450   //
1451   // Dump mean signal info
1452   //
1453     if (AliTPCReconstructor::StreamLevel()>0) {
1454     (*fDebugStreamer)<<"Signal"<<
1455     "TimeStamp="<<fTimeStamp<<
1456     "EventType="<<fEventType<<
1457     "Sector="<<uid[0]<<
1458     "Row="<<uid[1]<<
1459     "Pad="<<uid[2]<<
1460     "Max="<<max<<
1461     "MaxPos="<<maxPos<<
1462     //
1463     "Median="<<median<<
1464     "Mean="<<mean<<
1465     "RMS="<<rms<<
1466     "Mean06="<<mean06<<
1467     "RMS06="<<rms06<<
1468     "Mean09="<<mean09<<
1469     "RMS09="<<rms09<<
1470     "RMSCalib="<<rmsCalib<<
1471     "PedCalib="<<pedestalCalib<<
1472     "\n";
1473     }
1474   //
1475   // fill pedestal histogram
1476   //
1477   //
1478   //
1479   //
1480   Float_t kMin =fRecoParam->GetDumpAmplitudeMin();   // minimal signal to be dumped
1481   Float_t *dsignal = new Float_t[nchannels];
1482   Float_t *dtime   = new Float_t[nchannels];
1483   for (Int_t i=0; i<nchannels; i++){
1484     dtime[i] = i;
1485     dsignal[i] = signal[i];
1486   }
1487
1488   TGraph * graph=0;
1489   //
1490   // Big signals dumping
1491   //
1492   if (AliTPCReconstructor::StreamLevel()>0) {
1493   if (max-median>kMin &&maxPos>fRecoParam->GetFirstBin())
1494     (*fDebugStreamer)<<"SignalB"<<     // pads with signal
1495       "TimeStamp="<<fTimeStamp<<
1496       "EventType="<<fEventType<<
1497       "Sector="<<uid[0]<<
1498       "Row="<<uid[1]<<
1499       "Pad="<<uid[2]<<
1500       "Graph="<<graph<<
1501       "Max="<<max<<
1502       "MaxPos="<<maxPos<<
1503       //
1504       "Median="<<median<<
1505       "Mean="<<mean<<
1506       "RMS="<<rms<<
1507       "Mean06="<<mean06<<
1508       "RMS06="<<rms06<<
1509       "Mean09="<<mean09<<
1510       "RMS09="<<rms09<<
1511       "\n";
1512   delete graph;
1513   }
1514
1515   delete [] dsignal;
1516   delete [] dtime;
1517   if (rms06>fRecoParam->GetMaxNoise()) {
1518     pedestalEvent+=1024.;
1519     return 1024+median; // sign noisy channel in debug mode
1520   }
1521   return median;
1522 }
1523
1524
1525
1526
1527