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 //   Origin: Marian Ivanov
  22 //-------------------------------------------------------
  23
  24 #include "AliTPCReconstructor.h"
  25 #include "AliTPCclustererMI.h"
  26 #include "AliTPCclusterMI.h"
  27 #include <TObjArray.h>
  28 #include <TFile.h>
  29 #include "AliTPCClustersArray.h"
  30 #include "AliTPCClustersRow.h"
  31 #include "AliTPCRawStream.h"
  32 #include "AliDigits.h"
  33 #include "AliSimDigits.h"
  34 #include "AliTPCParam.h"
  35 #include "AliRawReader.h"
  36 #include "AliTPCRawStream.h"
  37 #include "AliRunLoader.h"
  38 #include "AliLoader.h"
  39 #include "Riostream.h"
  40 #include <TTree.h>
  41
  42 #include "AliTPCcalibDB.h"
  43 #include "AliTPCCalPad.h"
  44 #include "AliTPCCalROC.h"
  45
  46
  47 ClassImp(AliTPCclustererMI)
  48
  49
  50
  51 AliTPCclustererMI::AliTPCclustererMI(const AliTPCParam* par)
  52 {
  53   fInput =0;
  54   fOutput=0;
  55   fParam = par;
  56 }
  57 void AliTPCclustererMI::SetInput(TTree * tree)
  58 {
  59   //
  60   // set input tree with digits
  61   //
  62   fInput = tree;
  63   if  (!fInput->GetBranch("Segment")){
  64     cerr<<"AliTPC::Digits2Clusters(): no porper input tree !\n";
  65     fInput=0;
  66     return;
  67   }
  68 }
  69
  70 void AliTPCclustererMI::SetOutput(TTree * tree)
  71 {
  72   //
  73   //
  74   fOutput= tree;
  75   AliTPCClustersRow clrow;
  76   AliTPCClustersRow *pclrow=&clrow;
  77   clrow.SetClass("AliTPCclusterMI");
  78   clrow.SetArray(1); // to make Clones array
  79   fOutput->Branch("Segment","AliTPCClustersRow",&pclrow,32000,200);
  80 }
  81
  82
  83 Float_t  AliTPCclustererMI::GetSigmaY2(Int_t iz){
  84   // sigma y2 = in digits  - we don't know the angle
  85   Float_t z = iz*fParam->GetZWidth()+fParam->GetNTBinsL1()*fParam->GetZWidth();
  86   Float_t sd2 = (z*fParam->GetDiffL()*fParam->GetDiffL())/
  87     (fPadWidth*fPadWidth);
  88   Float_t sres = 0.25;
  89   Float_t res = sd2+sres;
  90   return res;
  91 }
  92
  93
  94 Float_t  AliTPCclustererMI::GetSigmaZ2(Int_t iz){
  95   //sigma z2 = in digits - angle estimated supposing vertex constraint
  96   Float_t z = iz*fZWidth+fParam->GetNTBinsL1()*fParam->GetZWidth();
  97   Float_t sd2 = (z*fParam->GetDiffL()*fParam->GetDiffL())/(fZWidth*fZWidth);
  98   Float_t angular = fPadLength*(fParam->GetZLength()-z)/(fRx*fZWidth);
  99   angular*=angular;
 100   angular/=12.;
 101   Float_t sres = fParam->GetZSigma()/fZWidth;
 102   sres *=sres;
 103   Float_t res = angular +sd2+sres;
 104   return res;
 105 }
 106
 107 void AliTPCclustererMI::MakeCluster(Int_t k,Int_t max,Float_t *bins, UInt_t /*m*/,
 108 AliTPCclusterMI &c)
 109 {
 110   Int_t i0=k/max;  //central pad
 111   Int_t j0=k%max;  //central time bin
 112
 113   // set pointers to data
 114   //Int_t dummy[5] ={0,0,0,0,0};
 115   Float_t * matrix[5]; //5x5 matrix with digits  - indexing i = 0 ..4  j = -2..2
 116   Float_t * resmatrix[5];
 117   for (Int_t di=-2;di<=2;di++){
 118     matrix[di+2]  =  &bins[k+di*max];
 119     resmatrix[di+2]  =  &fResBins[k+di*max];
 120   }
 121   //build matrix with virtual charge
 122   Float_t sigmay2= GetSigmaY2(j0);
 123   Float_t sigmaz2= GetSigmaZ2(j0);
 124
 125   Float_t vmatrix[5][5];
 126   vmatrix[2][2] = matrix[2][0];
 127   c.SetType(0);
 128   c.SetMax((UShort_t)(vmatrix[2][2])); // write maximal amplitude
 129   for (Int_t di =-1;di <=1;di++)
 130     for (Int_t dj =-1;dj <=1;dj++){
 131       Float_t amp = matrix[di+2][dj];
 132       if ( (amp<2) && (fLoop<2)){
 133         // if under threshold  - calculate virtual charge
 134         Float_t ratio = TMath::Exp(-1.2*TMath::Abs(di)/sigmay2)*TMath::Exp(-1.2*TMath::Abs(dj)/sigmaz2);
 135         amp = ((matrix[2][0]-2)*(matrix[2][0]-2)/(matrix[-di+2][-dj]+2))*ratio;
 136         if (amp>2) amp = 2;
 137         vmatrix[2+di][2+dj]=amp;
 138         vmatrix[2+2*di][2+2*dj]=0;
 139         if ( (di*dj)!=0){
 140           //DIAGONAL ELEMENTS
 141           vmatrix[2+2*di][2+dj] =0;
 142           vmatrix[2+di][2+2*dj] =0;
 143         }
 144         continue;
 145       }
 146       if (amp<4){
 147         //if small  amplitude - below  2 x threshold  - don't consider other one
 148         vmatrix[2+di][2+dj]=amp;
 149         vmatrix[2+2*di][2+2*dj]=0;  // don't take to the account next bin
 150         if ( (di*dj)!=0){
 151           //DIAGONAL ELEMENTS
 152           vmatrix[2+2*di][2+dj] =0;
 153           vmatrix[2+di][2+2*dj] =0;
 154         }
 155         continue;
 156       }
 157       //if bigger then take everything
 158       vmatrix[2+di][2+dj]=amp;
 159       vmatrix[2+2*di][2+2*dj]= matrix[2*di+2][2*dj] ;
 160       if ( (di*dj)!=0){
 161           //DIAGONAL ELEMENTS
 162           vmatrix[2+2*di][2+dj] = matrix[2*di+2][dj];
 163           vmatrix[2+di][2+2*dj] = matrix[2+di][dj*2];
 164         }
 165     }
 166
 167
 168
 169   Float_t sumw=0;
 170   Float_t sumiw=0;
 171   Float_t sumi2w=0;
 172   Float_t sumjw=0;
 173   Float_t sumj2w=0;
 174   //
 175   for (Int_t i=-2;i<=2;i++)
 176     for (Int_t j=-2;j<=2;j++){
 177       Float_t amp = vmatrix[i+2][j+2];
 178
 179       sumw    += amp;
 180       sumiw   += i*amp;
 181       sumi2w  += i*i*amp;
 182       sumjw   += j*amp;
 183       sumj2w  += j*j*amp;
 184     }
 185   //
 186   Float_t meani = sumiw/sumw;
 187   Float_t mi2   = sumi2w/sumw-meani*meani;
 188   Float_t meanj = sumjw/sumw;
 189   Float_t mj2   = sumj2w/sumw-meanj*meanj;
 190   //
 191   Float_t ry = mi2/sigmay2;
 192   Float_t rz = mj2/sigmaz2;
 193
 194   //
 195   if ( ( (ry<0.6) || (rz<0.6) ) && fLoop==2) return;
 196   if ( (ry <1.2) && (rz<1.2) ) {
 197     //if cluster looks like expected
 198     //+1.2 deviation from expected sigma accepted
 199     //    c.fMax = FitMax(vmatrix,meani,meanj,TMath::Sqrt(sigmay2),TMath::Sqrt(sigmaz2));
 200
 201     meani +=i0;
 202     meanj +=j0;
 203     //set cluster parameters
 204     c.SetQ(sumw);
 205     c.SetY(meani*fPadWidth);
 206     c.SetZ(meanj*fZWidth);
 207     c.SetSigmaY2(mi2);
 208     c.SetSigmaZ2(mj2);
 209     AddCluster(c);
 210     //remove cluster data from data
 211     for (Int_t di=-2;di<=2;di++)
 212       for (Int_t dj=-2;dj<=2;dj++){
 213         resmatrix[di+2][dj] -= vmatrix[di+2][dj+2];
 214         if (resmatrix[di+2][dj]<0) resmatrix[di+2][dj]=0;
 215       }
 216     resmatrix[2][0] =0;
 217     return;
 218   }
 219   //
 220   //unfolding when neccessary
 221   //
 222
 223   Float_t * matrix2[7]; //7x7 matrix with digits  - indexing i = 0 ..6  j = -3..3
 224   Float_t dummy[7]={0,0,0,0,0,0};
 225   for (Int_t di=-3;di<=3;di++){
 226     matrix2[di+3] =  &bins[k+di*max];
 227     if ((k+di*max)<3)  matrix2[di+3] = &dummy[3];
 228     if ((k+di*max)>fMaxBin-3)  matrix2[di+3] = &dummy[3];
 229   }
 230   Float_t vmatrix2[5][5];
 231   Float_t sumu;
 232   Float_t overlap;
 233   UnfoldCluster(matrix2,vmatrix2,meani,meanj,sumu,overlap);
 234   //
 235   //  c.fMax = FitMax(vmatrix2,meani,meanj,TMath::Sqrt(sigmay2),TMath::Sqrt(sigmaz2));
 236   meani +=i0;
 237   meanj +=j0;
 238   //set cluster parameters
 239   c.SetQ(sumu);
 240   c.SetY(meani*fPadWidth);
 241   c.SetZ(meanj*fZWidth);
 242   c.SetSigmaY2(mi2);
 243   c.SetSigmaZ2(mj2);
 244   c.SetType(Char_t(overlap)+1);
 245   AddCluster(c);
 246
 247   //unfolding 2
 248   meani-=i0;
 249   meanj-=j0;
 250   if (gDebug>4)
 251     printf("%f\t%f\n", vmatrix2[2][2], vmatrix[2][2]);
 252 }
 253
 254
 255
 256 void AliTPCclustererMI::UnfoldCluster(Float_t * matrix2[7], Float_t recmatrix[5][5], Float_t & meani, Float_t & meanj,
 257                                       Float_t & sumu, Float_t & overlap )
 258 {
 259   //
 260   //unfold cluster from input matrix
 261   //data corresponding to cluster writen in recmatrix
 262   //output meani and meanj
 263
 264   //take separatelly y and z
 265
 266   Float_t sum3i[7] = {0,0,0,0,0,0,0};
 267   Float_t sum3j[7] = {0,0,0,0,0,0,0};
 268
 269   for (Int_t k =0;k<7;k++)
 270     for (Int_t l = -1; l<=1;l++){
 271       sum3i[k]+=matrix2[k][l];
 272       sum3j[k]+=matrix2[l+3][k-3];
 273     }
 274   Float_t mratio[3][3]={{1,1,1},{1,1,1},{1,1,1}};
 275   //
 276   //unfold  y
 277   Float_t sum3wi    = 0;  //charge minus overlap
 278   Float_t sum3wio   = 0;  //full charge
 279   Float_t sum3iw    = 0;  //sum for mean value
 280   for (Int_t dk=-1;dk<=1;dk++){
 281     sum3wio+=sum3i[dk+3];
 282     if (dk==0){
 283       sum3wi+=sum3i[dk+3];
 284     }
 285     else{
 286       Float_t ratio =1;
 287       if (  ( ((sum3i[dk+3]+3)/(sum3i[3]-3))+1 < (sum3i[2*dk+3]-3)/(sum3i[dk+3]+3))||
 288            sum3i[dk+3]<=sum3i[2*dk+3] && sum3i[dk+3]>2 ){
 289         Float_t xm2 = sum3i[-dk+3];
 290         Float_t xm1 = sum3i[+3];
 291         Float_t x1  = sum3i[2*dk+3];
 292         Float_t x2  = sum3i[3*dk+3];
 293         Float_t w11   = TMath::Max((Float_t)(4.*xm1-xm2),(Float_t)0.000001);
 294         Float_t w12   = TMath::Max((Float_t)(4 *x1 -x2),(Float_t)0.);
 295         ratio = w11/(w11+w12);
 296         for (Int_t dl=-1;dl<=1;dl++)
 297           mratio[dk+1][dl+1] *= ratio;
 298       }
 299       Float_t amp = sum3i[dk+3]*ratio;
 300       sum3wi+=amp;
 301       sum3iw+= dk*amp;
 302     }
 303   }
 304   meani = sum3iw/sum3wi;
 305   Float_t overlapi = (sum3wio-sum3wi)/sum3wio;
 306
 307
 308
 309   //unfold  z
 310   Float_t sum3wj    = 0;  //charge minus overlap
 311   Float_t sum3wjo   = 0;  //full charge
 312   Float_t sum3jw    = 0;  //sum for mean value
 313   for (Int_t dk=-1;dk<=1;dk++){
 314     sum3wjo+=sum3j[dk+3];
 315     if (dk==0){
 316       sum3wj+=sum3j[dk+3];
 317     }
 318     else{
 319       Float_t ratio =1;
 320       if ( ( ((sum3j[dk+3]+3)/(sum3j[3]-3))+1 < (sum3j[2*dk+3]-3)/(sum3j[dk+3]+3)) ||
 321            (sum3j[dk+3]<=sum3j[2*dk+3] && sum3j[dk+3]>2)){
 322         Float_t xm2 = sum3j[-dk+3];
 323         Float_t xm1 = sum3j[+3];
 324         Float_t x1  = sum3j[2*dk+3];
 325         Float_t x2  = sum3j[3*dk+3];
 326         Float_t w11   = TMath::Max((Float_t)(4.*xm1-xm2),(Float_t)0.000001);
 327         Float_t w12   = TMath::Max((Float_t)(4 *x1 -x2),(Float_t)0.);
 328         ratio = w11/(w11+w12);
 329         for (Int_t dl=-1;dl<=1;dl++)
 330           mratio[dl+1][dk+1] *= ratio;
 331       }
 332       Float_t amp = sum3j[dk+3]*ratio;
 333       sum3wj+=amp;
 334       sum3jw+= dk*amp;
 335     }
 336   }
 337   meanj = sum3jw/sum3wj;
 338   Float_t overlapj = (sum3wjo-sum3wj)/sum3wjo;
 339   overlap = Int_t(100*TMath::Max(overlapi,overlapj)+3);
 340   sumu = (sum3wj+sum3wi)/2.;
 341
 342   if (overlap ==3) {
 343     //if not overlap detected remove everything
 344     for (Int_t di =-2; di<=2;di++)
 345       for (Int_t dj =-2; dj<=2;dj++){
 346         recmatrix[di+2][dj+2] = matrix2[3+di][dj];
 347       }
 348   }
 349   else{
 350     for (Int_t di =-1; di<=1;di++)
 351       for (Int_t dj =-1; dj<=1;dj++){
 352         Float_t ratio =1;
 353         if (mratio[di+1][dj+1]==1){
 354           recmatrix[di+2][dj+2]     = matrix2[3+di][dj];
 355           if (TMath::Abs(di)+TMath::Abs(dj)>1){
 356             recmatrix[2*di+2][dj+2] = matrix2[3+2*di][dj];
 357             recmatrix[di+2][2*dj+2] = matrix2[3+di][2*dj];
 358           }
 359           recmatrix[2*di+2][2*dj+2] = matrix2[3+2*di][2*dj];
 360         }
 361         else
 362           {
 363             //if we have overlap in direction
 364             recmatrix[di+2][dj+2] = mratio[di+1][dj+1]* matrix2[3+di][dj];
 365             if (TMath::Abs(di)+TMath::Abs(dj)>1){
 366               ratio =  TMath::Min((Float_t)(recmatrix[di+2][dj+2]/(matrix2[3+0*di][1*dj]+1)),(Float_t)1.);
 367               recmatrix[2*di+2][dj+2] = ratio*recmatrix[di+2][dj+2];
 368               //
 369               ratio =  TMath::Min((Float_t)(recmatrix[di+2][dj+2]/(matrix2[3+1*di][0*dj]+1)),(Float_t)1.);
 370               recmatrix[di+2][2*dj+2] = ratio*recmatrix[di+2][dj+2];
 371             }
 372             else{
 373               ratio =  recmatrix[di+2][dj+2]/matrix2[3][0];
 374               recmatrix[2*di+2][2*dj+2] = ratio*recmatrix[di+2][dj+2];
 375             }
 376           }
 377       }
 378   }
 379   if (gDebug>4)
 380     printf("%f\n", recmatrix[2][2]);
 381
 382 }
 383
 384 Float_t AliTPCclustererMI::FitMax(Float_t vmatrix[5][5], Float_t y, Float_t z, Float_t sigmay, Float_t sigmaz)
 385 {
 386   //
 387   // estimate max
 388   Float_t sumteor= 0;
 389   Float_t sumamp = 0;
 390
 391   for (Int_t di = -1;di<=1;di++)
 392     for (Int_t dj = -1;dj<=1;dj++){
 393       if (vmatrix[2+di][2+dj]>2){
 394         Float_t teor = TMath::Gaus(di,y,sigmay*1.2)*TMath::Gaus(dj,z,sigmaz*1.2);
 395         sumteor += teor*vmatrix[2+di][2+dj];
 396         sumamp  += vmatrix[2+di][2+dj]*vmatrix[2+di][2+dj];
 397       }
 398     }
 399   Float_t max = sumamp/sumteor;
 400   return max;
 401 }
 402
 403 void AliTPCclustererMI::AddCluster(AliTPCclusterMI &c){
 404   //
 405   // transform cluster to the global coordinata
 406   // add the cluster to the array
 407   //
 408   Float_t meani = c.GetY()/fPadWidth;
 409   Float_t meanj = c.GetZ()/fZWidth;
 410
 411   Int_t ki = TMath::Nint(meani-3);
 412   if (ki<0) ki=0;
 413   if (ki>=fMaxPad) ki = fMaxPad-1;
 414   Int_t kj = TMath::Nint(meanj-3);
 415   if (kj<0) kj=0;
 416   if (kj>=fMaxTime-3) kj=fMaxTime-4;
 417   // ki and kj shifted to "real" coordinata
 418   if (fRowDig) {
 419     c.SetLabel(fRowDig->GetTrackIDFast(kj,ki,0)-2,0);
 420     c.SetLabel(fRowDig->GetTrackIDFast(kj,ki,1)-2,1);
 421     c.SetLabel(fRowDig->GetTrackIDFast(kj,ki,2)-2,2);
 422   }
 423
 424
 425   Float_t s2 = c.GetSigmaY2();
 426   Float_t w=fParam->GetPadPitchWidth(fSector);
 427
 428   c.SetSigmaY2(s2*w*w);
 429   s2 = c.GetSigmaZ2();
 430   w=fZWidth;
 431   c.SetSigmaZ2(s2*w*w);
 432   c.SetY((meani - 2.5 - 0.5*fMaxPad)*fParam->GetPadPitchWidth(fSector));
 433   c.SetZ(fZWidth*(meanj-3));
 434   c.SetZ(c.GetZ() - 3.*fParam->GetZSigma() + fParam->GetNTBinsL1()*fParam->GetZWidth()); // PASA delay + L1 delay
 435   c.SetZ(fSign*(fParam->GetZLength() - c.GetZ()));
 436   c.SetX(fRx);
 437   c.SetDetector(fSector);
 438   c.SetRow(fRow);
 439
 440   if (ki<=1 || ki>=fMaxPad-1 || kj==1 || kj==fMaxTime-2) {
 441     //c.SetSigmaY2(c.GetSigmaY2()*25.);
 442     //c.SetSigmaZ2(c.GetSigmaZ2()*4.);
 443     c.SetType(-(c.GetType()+3));  //edge clusters
 444   }
 445   if (fLoop==2) c.SetType(100);
 446
 447   TClonesArray * arr = fRowCl->GetArray();
 448   // AliTPCclusterMI * cl =
 449   new ((*arr)[fNcluster]) AliTPCclusterMI(c);
 450
 451   fNcluster++;
 452 }
 453
 454
 455 //_____________________________________________________________________________
 456 void AliTPCclustererMI::Digits2Clusters()
 457 {
 458   //-----------------------------------------------------------------
 459   // This is a simple cluster finder.
 460   //-----------------------------------------------------------------
 461
 462   if (!fInput) {
 463     Error("Digits2Clusters", "input tree not initialised");
 464     return;
 465   }
 466
 467   if (!fOutput) {
 468     Error("Digits2Clusters", "output tree not initialised");
 469     return;
 470   }
 471
 472   AliTPCCalPad * gainTPC = AliTPCcalibDB::Instance()->GetPadGainFactor();
 473
 474   AliSimDigits digarr, *dummy=&digarr;
 475   fRowDig = dummy;
 476   fInput->GetBranch("Segment")->SetAddress(&dummy);
 477   Stat_t nentries = fInput->GetEntries();
 478
 479   fMaxTime=fParam->GetMaxTBin()+6; // add 3 virtual time bins before and 3   after
 480
 481   Int_t nclusters  = 0;
 482
 483   for (Int_t n=0; n<nentries; n++) {
 484     fInput->GetEvent(n);
 485     if (!fParam->AdjustSectorRow(digarr.GetID(),fSector,fRow)) {
 486       cerr<<"AliTPC warning: invalid segment ID ! "<<digarr.GetID()<<endl;
 487       continue;
 488     }
 489     Int_t row = fRow;
 490     AliTPCCalROC * gainROC = gainTPC->GetCalROC(fSector);  // pad gains per given sector
 491
 492     //
 493     AliTPCClustersRow *clrow= new AliTPCClustersRow();
 494     fRowCl = clrow;
 495     clrow->SetClass("AliTPCclusterMI");
 496     clrow->SetArray(1);
 497
 498     clrow->SetID(digarr.GetID());
 499     fOutput->GetBranch("Segment")->SetAddress(&clrow);
 500     fRx=fParam->GetPadRowRadii(fSector,row);
 501
 502
 503     const Int_t kNIS=fParam->GetNInnerSector(), kNOS=fParam->GetNOuterSector();
 504     fZWidth = fParam->GetZWidth();
 505     if (fSector < kNIS) {
 506       fMaxPad = fParam->GetNPadsLow(row);
 507       fSign = (fSector < kNIS/2) ? 1 : -1;
 508       fPadLength = fParam->GetPadPitchLength(fSector,row);
 509       fPadWidth = fParam->GetPadPitchWidth();
 510     } else {
 511       fMaxPad = fParam->GetNPadsUp(row);
 512       fSign = ((fSector-kNIS) < kNOS/2) ? 1 : -1;
 513       fPadLength = fParam->GetPadPitchLength(fSector,row);
 514       fPadWidth  = fParam->GetPadPitchWidth();
 515     }
 516
 517
 518     fMaxBin=fMaxTime*(fMaxPad+6);  // add 3 virtual pads  before and 3 after
 519     fBins    =new Float_t[fMaxBin];
 520     fResBins =new Float_t[fMaxBin];  //fBins with residuals after 1 finder loop
 521     memset(fBins,0,sizeof(Float_t)*fMaxBin);
 522     memset(fResBins,0,sizeof(Float_t)*fMaxBin);
 523
 524     if (digarr.First()) //MI change
 525       do {
 526         Float_t dig=digarr.CurrentDigit();
 527         if (dig<=fParam->GetZeroSup()) continue;
 528         Int_t j=digarr.CurrentRow()+3, i=digarr.CurrentColumn()+3;
 529         Float_t gain = gainROC->GetValue(row,digarr.CurrentColumn());
 530         fBins[i*fMaxTime+j]=dig/gain;
 531       } while (digarr.Next());
 532     digarr.ExpandTrackBuffer();
 533
 534     FindClusters();
 535
 536     fOutput->Fill();
 537     delete clrow;
 538     nclusters+=fNcluster;
 539     delete[] fBins;
 540     delete[] fResBins;
 541   }
 542
 543   Info("Digits2Clusters", "Number of found clusters : %d", nclusters);
 544 }
 545
 546 void AliTPCclustererMI::Digits2Clusters(AliRawReader* rawReader)
 547 {
 548 //-----------------------------------------------------------------
 549 // This is a cluster finder for raw data.
 550 //-----------------------------------------------------------------
 551
 552   if (!fOutput) {
 553     Error("Digits2Clusters", "output tree not initialised");
 554     return;
 555   }
 556
 557   AliTPCCalPad * gainTPC = AliTPCcalibDB::Instance()->GetPadGainFactor();
 558
 559   rawReader->Reset();
 560   AliTPCRawStream input(rawReader);
 561
 562   fRowDig = NULL;
 563
 564   Int_t nclusters  = 0;
 565
 566   fMaxTime = fParam->GetMaxTBin() + 6; // add 3 virtual time bins before and 3 after
 567   const Int_t kNIS = fParam->GetNInnerSector();
 568   const Int_t kNOS = fParam->GetNOuterSector();
 569   const Int_t kNS = kNIS + kNOS;
 570   fZWidth = fParam->GetZWidth();
 571   Int_t zeroSup = fParam->GetZeroSup();
 572
 573   fBins = NULL;
 574   Float_t** splitRows = new Float_t* [kNS*2];
 575   Float_t** splitRowsRes = new Float_t* [kNS*2];
 576   for (Int_t iSector = 0; iSector < kNS*2; iSector++)
 577     splitRows[iSector] = NULL;
 578   Int_t iSplitRow = -1;
 579
 580   Bool_t next = kTRUE;
 581   while (next) {
 582     next = input.Next();
 583
 584     // when the sector or row number has changed ...
 585     if (input.IsNewRow() || !next) {
 586
 587       // ... find clusters in the previous pad row, and ...
 588       if (fBins) {
 589         if ((iSplitRow < 0) || splitRows[fSector + kNS*iSplitRow]) {
 590           fRowCl = new AliTPCClustersRow;
 591           fRowCl->SetClass("AliTPCclusterMI");
 592           fRowCl->SetArray(1);
 593           fRowCl->SetID(fParam->GetIndex(fSector, input.GetPrevRow()));
 594           fOutput->GetBranch("Segment")->SetAddress(&fRowCl);
 595
 596           FindClusters();
 597
 598           fOutput->Fill();
 599           delete fRowCl;
 600           nclusters += fNcluster;
 601           delete[] fBins;
 602           delete[] fResBins;
 603           if (iSplitRow >= 0) splitRows[fSector + kNS*iSplitRow] = NULL;
 604
 605         } else if (iSplitRow >= 0) {
 606           splitRows[fSector + kNS*iSplitRow] = fBins;
 607           splitRowsRes[fSector + kNS*iSplitRow] = fResBins;
 608         }
 609       }
 610
 611       if (!next) break;
 612
 613       // ... prepare for the next pad row
 614       fSector = input.GetSector();
 615       fRow    = input.GetRow();
 616       Int_t iRow = input.GetRow();
 617       fRx = fParam->GetPadRowRadii(fSector, iRow);
 618
 619       iSplitRow = -1;
 620       if (fSector < kNIS) {
 621         fMaxPad = fParam->GetNPadsLow(iRow);
 622         fSign = (fSector < kNIS/2) ? 1 : -1;
 623         if (iRow == 30) iSplitRow = 0;
 624       } else {
 625         fMaxPad = fParam->GetNPadsUp(iRow);
 626         fSign = ((fSector-kNIS) < kNOS/2) ? 1 : -1;
 627         if (iRow == 27) iSplitRow = 0;
 628         else if (iRow == 76) iSplitRow = 1;
 629       }
 630       fPadLength = fParam->GetPadPitchLength(fSector, iRow);
 631       fPadWidth  = fParam->GetPadPitchWidth();
 632
 633       fMaxBin = fMaxTime*(fMaxPad+6);  // add 3 virtual pads  before and 3 after
 634       if ((iSplitRow < 0) || !splitRows[fSector + kNS*iSplitRow]) {
 635         fBins    = new Float_t[fMaxBin];
 636         fResBins = new Float_t[fMaxBin];  //fBins with residuals after 1 finder loop
 637         memset(fBins, 0, sizeof(Float_t)*fMaxBin);
 638         memset(fResBins, 0, sizeof(Float_t)*fMaxBin);
 639       } else {
 640         fBins    = splitRows[fSector + kNS*iSplitRow];
 641         fResBins = splitRowsRes[fSector + kNS*iSplitRow];
 642       }
 643     }
 644
 645     // fill fBins with digits data
 646     if (input.GetSignal() <= zeroSup) continue;
 647     Int_t i = input.GetPad() + 3;
 648     Int_t j = input.GetTime() + 3;
 649     AliTPCCalROC * gainROC = gainTPC->GetCalROC(fSector);  // pad gains per given sector
 650     Float_t gain = gainROC->GetValue(fRow,input.GetPad());
 651     fBins[i*fMaxTime+j] = input.GetSignal()/gain;
 652   }
 653
 654   // find clusters in split rows that were skipped until now.
 655   // this can happen if the rows were not splitted
 656   for (fSector = 0; fSector < kNS; fSector++)
 657     for (Int_t iSplit = 0; iSplit < 2; iSplit++)
 658       if (splitRows[fSector + kNS*iSplit]) {
 659
 660         Int_t iRow = -1;
 661         if (fSector < kNIS) {
 662           iRow = 30;
 663           fMaxPad = fParam->GetNPadsLow(iRow);
 664           fSign = (fSector < kNIS/2) ? 1 : -1;
 665         } else {
 666           if (iSplit == 0) iRow = 27; else iRow = 76;
 667           fMaxPad = fParam->GetNPadsUp(iRow);
 668           fSign = ((fSector-kNIS) < kNOS/2) ? 1 : -1;
 669         }
 670         fRx = fParam->GetPadRowRadii(fSector, iRow);
 671         fPadLength = fParam->GetPadPitchLength(fSector, iRow);
 672         fPadWidth  = fParam->GetPadPitchWidth();
 673
 674         fMaxBin = fMaxTime*(fMaxPad+6);  // add 3 virtual pads  before and 3 after
 675         fBins    = splitRows[fSector + kNS*iSplit];
 676         fResBins = splitRowsRes[fSector + kNS*iSplit];
 677
 678         fRowCl = new AliTPCClustersRow;
 679         fRowCl->SetClass("AliTPCclusterMI");
 680         fRowCl->SetArray(1);
 681         fRowCl->SetID(fParam->GetIndex(fSector, iRow));
 682         fOutput->GetBranch("Segment")->SetAddress(&fRowCl);
 683
 684         FindClusters();
 685
 686         fOutput->Fill();
 687         delete fRowCl;
 688         nclusters += fNcluster;
 689         delete[] fBins;
 690         delete[] fResBins;
 691       }
 692
 693   delete[] splitRows;
 694   delete[] splitRowsRes;
 695   Info("Digits2Clusters", "Number of found clusters : %d\n", nclusters);
 696 }
 697
 698 void AliTPCclustererMI::FindClusters()
 699 {
 700   //add virtual charge at the edge
 701   for (Int_t i=0; i<fMaxTime; i++){
 702     Float_t amp1 = fBins[i+3*fMaxTime];
 703     Float_t amp0 =0;
 704     if (amp1>0){
 705       Float_t amp2 = fBins[i+4*fMaxTime];
 706       if (amp2==0) amp2=0.5;
 707       Float_t sigma2 = GetSigmaY2(i);
 708       amp0 = (amp1*amp1/amp2)*TMath::Exp(-1./sigma2);
 709       if (gDebug>4) printf("\n%f\n",amp0);
 710     }
 711     fBins[i+2*fMaxTime] = amp0;
 712     amp0 = 0;
 713     amp1 = fBins[(fMaxPad+2)*fMaxTime+i];
 714     if (amp1>0){
 715       Float_t amp2 = fBins[i+(fMaxPad+1)*fMaxTime];
 716       if (amp2==0) amp2=0.5;
 717       Float_t sigma2 = GetSigmaY2(i);
 718       amp0 = (amp1*amp1/amp2)*TMath::Exp(-1./sigma2);
 719       if (gDebug>4) printf("\n%f\n",amp0);
 720     }
 721     fBins[(fMaxPad+3)*fMaxTime+i] = amp0;
 722   }
 723
 724 //  memcpy(fResBins,fBins, fMaxBin*2);
 725   memcpy(fResBins,fBins, fMaxBin);
 726   //
 727   fNcluster=0;
 728   //first loop - for "gold cluster"
 729   fLoop=1;
 730   Float_t *b=&fBins[-1]+2*fMaxTime;
 731   Int_t crtime = Int_t((fParam->GetZLength()-AliTPCReconstructor::GetCtgRange()*fRx)/fZWidth-fParam->GetNTBinsL1()-5);
 732
 733   for (Int_t i=2*fMaxTime; i<fMaxBin-2*fMaxTime; i++) {
 734     b++;
 735     if (*b<8) continue;   //threshold form maxima
 736     if (i%fMaxTime<crtime) {
 737       Int_t delta = -(i%fMaxTime)+crtime;
 738       b+=delta;
 739       i+=delta;
 740       continue;
 741     }
 742
 743     if (!IsMaximum(*b,fMaxTime,b)) continue;
 744     AliTPCclusterMI c;
 745     Int_t dummy=0;
 746     MakeCluster(i, fMaxTime, fBins, dummy,c);
 747     //}
 748   }
 749   //memcpy(fBins,fResBins, fMaxBin*2);
 750   //second  loop - for rest cluster
 751   /*
 752   fLoop=2;
 753   b=&fResBins[-1]+2*fMaxTime;
 754   for (Int_t i=2*fMaxTime; i<fMaxBin-2*fMaxTime; i++) {
 755     b++;
 756     if (*b<25) continue;   // bigger threshold for maxima
 757     if (!IsMaximum(*b,fMaxTime,b)) continue;
 758     AliTPCclusterMI c;
 759     Int_t dummy;
 760     MakeCluster(i, fMaxTime, fResBins, dummy,c);
 761     //}
 762   }
 763   */
 764 }