TPC/Rec/AliTPCclusterer.cxx

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