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