PWG/EMCAL/AliAnalysisTaskEMCALClusterizeFast.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 //_________________________________________________________________________
  17 // This analysis provides a new list of clusters to be used in other analysis
  18 //
  19 // Author: Constantin Loizides, Salvatore Aiola
  20 //         Adapted from analysis class from Deepa Thomas
  21 //
  22 // $Id$
  23 //_________________________________________________________________________
  24
  25 // --- Root ---
  26 #include <TClonesArray.h>
  27 #include <TGeoManager.h>
  28 #include <TObjArray.h>
  29 #include <TString.h>
  30 #include <TTree.h>
  31
  32 // --- AliRoot ---
  33 #include "AliAODCaloCluster.h"
  34 #include "AliAODEvent.h"
  35 #include "AliAnalysisManager.h"
  36 #include "AliCDBEntry.h"
  37 #include "AliCDBManager.h"
  38 #include "AliCaloCalibPedestal.h"
  39 #include "AliEMCALAfterBurnerUF.h"
  40 #include "AliEMCALCalibData.h"
  41 #include "AliEMCALClusterizerNxN.h"
  42 #include "AliEMCALClusterizerv1.h"
  43 #include "AliEMCALClusterizerv2.h"
  44 #include "AliEMCALClusterizerFixedWindow.h"
  45 #include "AliEMCALDigit.h"
  46 #include "AliEMCALGeometry.h"
  47 #include "AliEMCALRecParam.h"
  48 #include "AliEMCALRecPoint.h"
  49 #include "AliEMCALRecoUtils.h"
  50 #include "AliESDEvent.h"
  51 #include "AliInputEventHandler.h"
  52 #include "AliLog.h"
  53
  54 #include "AliAnalysisTaskEMCALClusterizeFast.h"
  55
  56 ClassImp(AliAnalysisTaskEMCALClusterizeFast)
  57
  58 //________________________________________________________________________
  59 AliAnalysisTaskEMCALClusterizeFast::AliAnalysisTaskEMCALClusterizeFast()
  60   : AliAnalysisTaskSE(),
  61     fRun(-1),
  62     fDigitsArr(0),
  63     fClusterArr(0),
  64     fRecParam(0),
  65     fClusterizer(0),
  66     fUnfolder(0),
  67     fJustUnfold(kFALSE),
  68     fGeomName(),
  69     fGeomMatrixSet(kFALSE),
  70     fLoadGeomMatrices(kFALSE),
  71     fOCDBpath(),
  72     fCalibData(0),
  73     fPedestalData(0),
  74     fOutputAODBranch(0),
  75     fOutputAODBrName(),
  76     fRecoUtils(0),
  77     fLoadCalib(0),
  78     fLoadPed(0),
  79     fAttachClusters(1),
  80     fRecalibOnly(0),
  81     fSubBackground(0),
  82     fOverwrite(0),
  83     fNewClusterArrayName("newCaloClusters"),
  84     fNPhi(4),
  85     fNEta(4),
  86     fShiftPhi(2),
  87     fShiftEta(2),
  88     fTRUShift(0),
  89     fInputCellType(kFEEData),
  90     fTrackName()
  91 {
  92   // Constructor
  93
  94   for(Int_t i = 0; i < 12; ++i)
  95     fGeomMatrix[i] = 0;
  96 }
  97
  98 //________________________________________________________________________
  99 AliAnalysisTaskEMCALClusterizeFast::AliAnalysisTaskEMCALClusterizeFast(const char *name)
 100   : AliAnalysisTaskSE(name),
 101     fRun(-1),
 102     fDigitsArr(0),
 103     fClusterArr(0),
 104     fRecParam(new AliEMCALRecParam),
 105     fClusterizer(0),
 106     fUnfolder(0),
 107     fJustUnfold(kFALSE),
 108     fGeomName(""),
 109     fGeomMatrixSet(kFALSE),
 110     fLoadGeomMatrices(kFALSE),
 111     fOCDBpath(),
 112     fCalibData(0),
 113     fPedestalData(0),
 114     fOutputAODBranch(0),
 115     fOutputAODBrName(),
 116     fRecoUtils(0),
 117     fLoadCalib(0),
 118     fLoadPed(0),
 119     fAttachClusters(1),
 120     fRecalibOnly(0),
 121     fSubBackground(0),
 122     fOverwrite(0),
 123     fNewClusterArrayName("newCaloClusters"),
 124     fNPhi(4),
 125     fNEta(4),
 126     fShiftPhi(2),
 127     fShiftEta(2),
 128     fTRUShift(0),
 129     fInputCellType(kFEEData),
 130     fTrackName()
 131 {
 132   // Constructor
 133
 134   fBranchNames     = "ESD:AliESDHeader.,AliESDRun.,EMCALCells.,EMCALTrigger. AOD:header,emcalCells";
 135   for(Int_t i = 0; i < 12; ++i)
 136     fGeomMatrix[i] = 0;
 137 }
 138
 139 //________________________________________________________________________
 140 AliAnalysisTaskEMCALClusterizeFast::~AliAnalysisTaskEMCALClusterizeFast()
 141 {
 142   // Destructor.
 143
 144   delete fClusterizer;
 145   delete fUnfolder;
 146   delete fRecoUtils;
 147 }
 148
 149 //-------------------------------------------------------------------
 150 void AliAnalysisTaskEMCALClusterizeFast::UserCreateOutputObjects()
 151 {
 152   // Create output objects.
 153
 154   if (!fOutputAODBrName.IsNull()) {
 155     fOutputAODBranch = new TClonesArray("AliAODCaloCluster", 0);
 156     fOutputAODBranch->SetName(fOutputAODBrName);
 157     AddAODBranch("TClonesArray", &fOutputAODBranch);
 158     AliInfo(Form("Created Branch: %s",fOutputAODBrName.Data()));
 159   }
 160 }
 161
 162 //________________________________________________________________________
 163 void AliAnalysisTaskEMCALClusterizeFast::UserExec(Option_t *)
 164 {
 165   // Main loop, called for each event
 166
 167   // remove the contents of output list set in the previous event
 168   if (fOutputAODBranch)
 169     fOutputAODBranch->Clear("C");
 170
 171   AliESDEvent *esdevent = dynamic_cast<AliESDEvent*>(InputEvent());
 172   AliAODEvent *aodevent = dynamic_cast<AliAODEvent*>(InputEvent());
 173
 174   if (!esdevent&&!aodevent) {
 175     Error("UserExec","Event not available");
 176     return;
 177   }
 178
 179   LoadBranches();
 180
 181   UInt_t offtrigger = 0;
 182   if (esdevent) {
 183     UInt_t mask1 = esdevent->GetESDRun()->GetDetectorsInDAQ();
 184     UInt_t mask2 = esdevent->GetESDRun()->GetDetectorsInReco();
 185     Bool_t desc1 = (mask1 >> 18) & 0x1;
 186     Bool_t desc2 = (mask2 >> 18) & 0x1;
 187     if (desc1==0 || desc2==0) { //AliDAQ::OfflineModuleName(180=="EMCAL"
 188       AliError(Form("EMCAL not in DAQ/RECO: %u (%u)/%u (%u)",
 189                     mask1, esdevent->GetESDRun()->GetDetectorsInReco(),
 190                     mask2, esdevent->GetESDRun()->GetDetectorsInDAQ()));
 191       return;
 192     }
 193     AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
 194     offtrigger = ((AliInputEventHandler*)(am->GetInputEventHandler()))->IsEventSelected();
 195   } else if (aodevent) {
 196     offtrigger =  aodevent->GetHeader()->GetOfflineTrigger();
 197   }
 198
 199   if (esdevent) {
 200     AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
 201     Bool_t mcmode = 0;
 202     if (am->GetMCtruthEventHandler())
 203       mcmode = 1;
 204     if (!mcmode) {
 205       if (offtrigger & AliVEvent::kFastOnly) {
 206         AliWarning(Form("EMCAL not in fast only partition"));
 207         return;
 208       }
 209     }
 210   }
 211
 212   Init();
 213
 214   if (fJustUnfold) {
 215     AliWarning("Unfolding not implemented");
 216     return;
 217   }
 218
 219   FillDigitsArray();
 220
 221   if (fRecalibOnly) {
 222     UpdateCells();
 223     return; // not requested to run clusterizer
 224   }
 225
 226   Clusterize();
 227   UpdateCells();
 228   UpdateClusters();
 229
 230   if (fOutputAODBranch)
 231     RecPoints2Clusters(fOutputAODBranch);
 232 }
 233
 234 //________________________________________________________________________
 235 void AliAnalysisTaskEMCALClusterizeFast::Clusterize()
 236 {
 237   // Clusterize
 238
 239   if (fSubBackground) {
 240     fClusterizer->SetInputCalibrated(kTRUE);
 241     fClusterizer->SetCalibrationParameters(0);
 242   }
 243
 244   fClusterizer->Digits2Clusters("");
 245   if (fSubBackground) {
 246     if (fCalibData) {
 247       fClusterizer->SetInputCalibrated(kFALSE);
 248       fClusterizer->SetCalibrationParameters(fCalibData);
 249     }
 250   }
 251 }
 252
 253 //________________________________________________________________________
 254 void AliAnalysisTaskEMCALClusterizeFast::FillDigitsArray()
 255 {
 256   // Fill digits array
 257
 258   AliEMCALGeometry *geom = AliEMCALGeometry::GetInstance();
 259
 260   fDigitsArr->Clear("C");
 261
 262   switch (fInputCellType) {
 263
 264   case kFEEData :
 265     {
 266       AliVCaloCells *cells = InputEvent()->GetEMCALCells();
 267       Double_t avgE        = 0; // for background subtraction
 268       const Int_t ncells   = cells->GetNumberOfCells();
 269       for (Int_t icell = 0, idigit = 0; icell < ncells; ++icell) {
 270         Double_t cellAmplitude=0, cellTime=0, cellEFrac = 0;
 271         Short_t  cellNumber=0, cellMCLabel=-1;
 272         if (cells->GetCell(icell, cellNumber, cellAmplitude, cellTime, cellMCLabel, cellEFrac) != kTRUE)
 273         break;
 274         AliEMCALDigit *digit = static_cast<AliEMCALDigit*>(fDigitsArr->New(idigit));
 275         digit->SetId(cellNumber);
 276         digit->SetTime(cellTime);
 277         digit->SetTimeR(cellTime);
 278         digit->SetIndexInList(idigit);
 279         digit->SetType(AliEMCALDigit::kHG);
 280         if (fRecalibOnly||fSubBackground) {
 281           Float_t energy = cellAmplitude;
 282           Float_t time   = cellTime;
 283           fClusterizer->Calibrate(energy,time,cellNumber);
 284           digit->SetAmplitude(energy);
 285         avgE += energy;
 286         } else {
 287           digit->SetAmplitude(cellAmplitude);
 288         }
 289         idigit++;
 290       }
 291
 292       fDigitsArr->Sort();
 293
 294       if (fSubBackground) {
 295         avgE /= geom->GetNumberOfSuperModules()*48*24;
 296         Int_t ndigis = fDigitsArr->GetEntries();
 297         for (Int_t i = 0; i < ndigis; ++i) {
 298         AliEMCALDigit *digit = static_cast<AliEMCALDigit*>(fDigitsArr->At(i));
 299         Double_t energy = digit->GetAmplitude() - avgE;
 300         if (energy<=0.001) {
 301           digit->SetAmplitude(0);
 302         } else {
 303           digit->SetAmplitude(energy);
 304         }
 305         }
 306       }
 307     }
 308     break;
 309
 310   case kPattern :
 311     {
 312       // Fill digits from a pattern
 313       Int_t maxd = geom->GetNCells() / 4;
 314       for (Int_t idigit = 0; idigit < maxd; idigit++){
 315         if (idigit % 24 == 12) idigit += 12;
 316         AliEMCALDigit *digit = static_cast<AliEMCALDigit*>(fDigitsArr->New(idigit));
 317         digit->SetId(idigit * 4);
 318         digit->SetTime(600);
 319         digit->SetTimeR(600);
 320         digit->SetIndexInList(idigit);
 321         digit->SetType(AliEMCALDigit::kHG);
 322         digit->SetAmplitude(0.1);
 323       }
 324     }
 325     break;
 326
 327   case kL0FastORs    :
 328   case kL0FastORsTC  :
 329   case kL1FastORs    :
 330     {
 331       // Fill digits from FastORs
 332
 333       AliVCaloTrigger *triggers = InputEvent()->GetCaloTrigger("EMCAL");
 334
 335       if (!triggers || !(triggers->GetEntries() > 0))
 336         return;
 337
 338       Int_t idigit = 0;
 339       triggers->Reset();
 340
 341       while ((triggers->Next())) {
 342         Float_t L0Amplitude = 0;
 343         triggers->GetAmplitude(L0Amplitude);
 344
 345         if (L0Amplitude <= 0 && fInputCellType != kL1FastORs)
 346           continue;
 347
 348         Int_t L1Amplitude = 0;
 349         triggers->GetL1TimeSum(L1Amplitude);
 350
 351         if (L1Amplitude <= 0 && fInputCellType == kL1FastORs)
 352           continue;
 353
 354         Int_t triggerTime = 0;
 355         Int_t ntimes = 0;
 356         triggers->GetNL0Times(ntimes);
 357
 358         if (ntimes < 1 && fInputCellType == kL0FastORsTC)
 359           continue;
 360
 361         if (ntimes > 0) {
 362           Int_t trgtimes[25];
 363           triggers->GetL0Times(trgtimes);
 364           triggerTime = trgtimes[0];
 365         }
 366
 367         Int_t triggerCol = 0, triggerRow = 0;
 368         triggers->GetPosition(triggerCol, triggerRow);
 369
 370         Int_t find = -1;
 371         geom->GetAbsFastORIndexFromPositionInEMCAL(triggerCol, triggerRow, find);
 372
 373         if (find < 0)
 374           continue;
 375
 376         Int_t cidx[4] = {-1};
 377         Bool_t ret = geom->GetCellIndexFromFastORIndex(find, cidx);
 378
 379         if (!ret)
 380           continue;
 381
 382         Float_t triggerAmplitude = 0;
 383
 384         if (fInputCellType == kL1FastORs) {
 385           triggerAmplitude = 0.25 * L1Amplitude;  // it will add 4 cells for 1 amplitude
 386         }
 387         else {
 388           triggerAmplitude = L0Amplitude;      // 10 bit truncated, so it is already divided by 4
 389         }
 390
 391         for (Int_t idxpos = 0; idxpos < 4; idxpos++) {
 392           Int_t triggerNumber = cidx[idxpos];
 393           AliEMCALDigit *digit = static_cast<AliEMCALDigit*>(fDigitsArr->New(idigit));
 394           digit->SetId(triggerNumber);
 395           digit->SetTime(triggerTime);
 396           digit->SetTimeR(triggerTime);
 397           digit->SetIndexInList(idigit);
 398           digit->SetType(AliEMCALDigit::kHG);
 399           digit->SetAmplitude(triggerAmplitude);
 400           idigit++;
 401         }
 402       }
 403     }
 404     break;
 405   }
 406 }
 407
 408 //________________________________________________________________________________________
 409 void AliAnalysisTaskEMCALClusterizeFast::RecPoints2Clusters(TClonesArray *clus)
 410 {
 411   // Cluster energy, global position, cells and their amplitude fractions are restored.
 412
 413   AliVCaloCells *cells = InputEvent()->GetEMCALCells();
 414   AliEMCALGeometry *geom = AliEMCALGeometry::GetInstance();
 415   AliVEvent *event = InputEvent();
 416   if (!event) {
 417     AliError("Cannot get the event");
 418     return;
 419   }
 420
 421   // tracks array for track/cluster matching
 422   TClonesArray *tarr = 0;
 423   if (!fTrackName.IsNull()) {
 424     tarr = dynamic_cast<TClonesArray*>(event->FindListObject(fTrackName));
 425     if (!tarr) {
 426       AliError(Form("Cannot get tracks named %s", fTrackName.Data()));
 427     }
 428   }
 429
 430   const Int_t Ncls = fClusterArr->GetEntries();
 431   AliDebug(1, Form("total no of clusters %d", Ncls));
 432   for(Int_t i=0, nout=clus->GetEntries(); i < Ncls; ++i) {
 433     AliEMCALRecPoint *recpoint = static_cast<AliEMCALRecPoint*>(fClusterArr->At(i));
 434     Int_t ncells_true = 0;
 435     const Int_t ncells = recpoint->GetMultiplicity();
 436     UShort_t   absIds[ncells];
 437     Double32_t ratios[ncells];
 438     Int_t *dlist = recpoint->GetDigitsList();
 439     Float_t *elist = recpoint->GetEnergiesList();
 440     for (Int_t c = 0; c < ncells; ++c) {
 441       AliEMCALDigit *digit = static_cast<AliEMCALDigit*>(fDigitsArr->At(dlist[c]));
 442       absIds[ncells_true] = digit->GetId();
 443       ratios[ncells_true] = elist[c]/digit->GetAmplitude();
 444       ++ncells_true;
 445     }
 446
 447     if (ncells_true < 1) {
 448       AliWarning("Skipping cluster with no cells");
 449       continue;
 450     }
 451
 452     // calculate new cluster position
 453     TVector3 gpos;
 454     recpoint->GetGlobalPosition(gpos);
 455     Float_t g[3];
 456     gpos.GetXYZ(g);
 457
 458     AliDebug(1, Form("energy %f", recpoint->GetEnergy()));
 459
 460     AliVCluster *c = static_cast<AliVCluster*>(clus->New(nout++));
 461     c->SetType(AliVCluster::kEMCALClusterv1);
 462     c->SetE(recpoint->GetEnergy());
 463     c->SetPosition(g);
 464     c->SetNCells(ncells_true);
 465     c->SetCellsAbsId(absIds);
 466     c->SetCellsAmplitudeFraction(ratios);
 467     c->SetID(recpoint->GetUniqueID());
 468     c->SetDispersion(recpoint->GetDispersion());
 469     c->SetEmcCpvDistance(-1);
 470     c->SetChi2(-1);
 471     c->SetTOF(recpoint->GetTime()) ;     //time-of-flight
 472     c->SetNExMax(recpoint->GetNExMax()); //number of local maxima
 473     Float_t elipAxis[2];
 474     recpoint->GetElipsAxis(elipAxis);
 475     c->SetM02(elipAxis[0]*elipAxis[0]);
 476     c->SetM20(elipAxis[1]*elipAxis[1]);
 477     if (fPedestalData) {
 478       c->SetDistanceToBadChannel(recpoint->GetDistanceToBadTower());
 479     } else {
 480       if (fRecoUtils && fRecoUtils->IsBadChannelsRemovalSwitchedOn()) {
 481         fRecoUtils->RecalculateClusterDistanceToBadChannel(geom, cells, c);
 482       }
 483     }
 484
 485     if (tarr) {
 486       Double_t dEtaMin  = 1e9;
 487       Double_t dPhiMin  = 1e9;
 488       Int_t    imin     = -1;
 489       Double_t ceta     = gpos.Eta();
 490       Double_t cphi     = gpos.Phi();
 491       const Int_t ntrks = tarr->GetEntries();
 492       for(Int_t t = 0; t<ntrks; ++t) {
 493         AliVTrack *track = static_cast<AliVTrack*>(tarr->At(t));
 494         if (!track)
 495           continue;
 496         const AliExternalTrackParam *outp = track->GetOuterParam();
 497         if (!outp)
 498           continue;
 499         Double_t trkPos[3] = {0.,0.,0.};
 500         if (!outp->GetXYZ(trkPos))
 501           continue;
 502         TVector3 vec(trkPos[0],trkPos[1],trkPos[2]);
 503         Double_t veta = vec.Eta();
 504         Double_t vphi = vec.Phi();
 505         if(vphi<0)
 506           vphi += 2*TMath::Pi();
 507         if (TMath::Abs(veta)>0.75 || (vphi<70*TMath::DegToRad()) || (vphi>190*TMath::DegToRad()))
 508           continue;
 509         Double_t dR = vec.DeltaR(gpos);
 510         if(dR > 25)
 511           continue;
 512         Float_t tmpEta=0, tmpPhi=0;
 513         if (0) {
 514           AliExternalTrackParam trkParTemp(*outp); // retrieve the starting point every time before the extrapolation
 515           Bool_t ret = fRecoUtils->ExtrapolateTrackToCluster(&trkParTemp, c, fRecoUtils->GetMass(), fRecoUtils->GetStep(), tmpEta, tmpPhi);
 516           if (!ret)
 517             continue;
 518         } else {
 519           tmpEta = ceta - veta;
 520           tmpPhi = cphi - vphi;
 521         }
 522         if (TMath::Abs(tmpEta)<TMath::Abs(dEtaMin) && TMath::Abs(tmpPhi)<TMath::Abs(dPhiMin)) {
 523           dEtaMin = tmpEta;
 524           dPhiMin = tmpPhi;
 525           imin = t;
 526         }
 527       }
 528       c->SetEmcCpvDistance(imin);
 529       c->SetTrackDistance(dPhiMin, dEtaMin);
 530     }
 531   }
 532 }
 533
 534 //________________________________________________________________________
 535 void AliAnalysisTaskEMCALClusterizeFast::UpdateCells()
 536 {
 537   // Update cells in case re-calibration was done.
 538
 539   if (!fCalibData&&!fSubBackground)
 540     return;
 541
 542   AliVCaloCells *cells = InputEvent()->GetEMCALCells();
 543   const Int_t   ncells = cells->GetNumberOfCells();
 544   const Int_t   ndigis = fDigitsArr->GetEntries();
 545   if (ncells!=ndigis) {
 546     cells->DeleteContainer();
 547     cells->CreateContainer(ndigis);
 548   }
 549   for (Int_t idigit = 0; idigit < ndigis; ++idigit) {
 550     AliEMCALDigit *digit = static_cast<AliEMCALDigit*>(fDigitsArr->At(idigit));
 551     Double_t cellAmplitude = digit->GetCalibAmp();
 552     Short_t cellNumber     = digit->GetId();
 553     Double_t cellTime      = digit->GetTime();
 554     cells->SetCell(idigit, cellNumber, cellAmplitude, cellTime);
 555   }
 556 }
 557
 558 //________________________________________________________________________
 559 void AliAnalysisTaskEMCALClusterizeFast::UpdateClusters()
 560 {
 561   // Update cells in case re-calibration was done.
 562
 563   if (!fAttachClusters)
 564     return;
 565
 566   TClonesArray *clus = 0;
 567
 568   if (fOverwrite) {
 569     clus = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject("caloClusters"));
 570     if (!clus)
 571       clus = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject("CaloClusters"));
 572     if (!clus)
 573       return;
 574
 575     const Int_t nents = clus->GetEntries();
 576     for (Int_t i=0;i<nents;++i) {
 577       AliVCluster *c = static_cast<AliVCluster*>(clus->At(i));
 578       if (!c)
 579         continue;
 580       if (c->IsEMCAL())
 581         delete clus->RemoveAt(i);
 582     }
 583   } else {
 584     clus = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fNewClusterArrayName));
 585     if (!clus) {
 586       clus = new TClonesArray("AliESDCaloCluster");
 587       clus->SetName(fNewClusterArrayName);
 588       InputEvent()->AddObject(clus);
 589     } else {
 590       clus->Delete();
 591     }
 592   }
 593
 594   RecPoints2Clusters(clus);
 595 }
 596
 597 //________________________________________________________________________________________
 598 void AliAnalysisTaskEMCALClusterizeFast::Init()
 599 {
 600   // Select clusterization/unfolding algorithm and set all the needed parameters.
 601
 602   AliVEvent * event = InputEvent();
 603   if (!event) {
 604     AliWarning("Event not available!!!");
 605     return;
 606   }
 607
 608   if (event->GetRunNumber()==fRun)
 609     return;
 610   fRun = event->GetRunNumber();
 611
 612   if (fJustUnfold){
 613     // init the unfolding afterburner
 614     delete fUnfolder;
 615     fUnfolder = new AliEMCALAfterBurnerUF(fRecParam->GetW0(),fRecParam->GetLocMaxCut(),fRecParam->GetMinECut());
 616     return;
 617   }
 618
 619   AliEMCALGeometry *geometry = 0;
 620   if (fGeomName.Length()>0)
 621     geometry = AliEMCALGeometry::GetInstance(fGeomName);
 622   else
 623     geometry = AliEMCALGeometry::GetInstance();
 624   if (!geometry) {
 625     AliFatal("Geometry not available!!!");
 626     return;
 627   }
 628
 629   if (!fGeomMatrixSet) {
 630     if (fLoadGeomMatrices) {
 631       for(Int_t mod=0; mod < geometry->GetNumberOfSuperModules(); ++mod) {
 632         if (fGeomMatrix[mod]){
 633           if (DebugLevel() > 2)
 634             fGeomMatrix[mod]->Print();
 635           geometry->SetMisalMatrix(fGeomMatrix[mod],mod);
 636         }
 637       }
 638     } else { // get matrix from file (work around bug in aliroot)
 639       for(Int_t mod=0; mod < geometry->GetEMCGeometry()->GetNumberOfSuperModules(); ++mod) {
 640         const TGeoHMatrix *gm = 0;
 641         AliESDEvent *esdevent = dynamic_cast<AliESDEvent*>(event);
 642         if (esdevent) {
 643           gm = esdevent->GetEMCALMatrix(mod);
 644         } else {
 645           AliAODHeader *aodheader = dynamic_cast<AliAODHeader*>(event->GetHeader());
 646           if (aodheader) {
 647             gm = aodheader->GetEMCALMatrix(mod);
 648           }
 649         }
 650         if (gm) {
 651           if (DebugLevel() > 2)
 652             gm->Print();
 653           geometry->SetMisalMatrix(gm,mod);
 654         }
 655       }
 656     }
 657     fGeomMatrixSet=kTRUE;
 658   }
 659
 660   // setup digit array if needed
 661   if (!fDigitsArr) {
 662     fDigitsArr = new TClonesArray("AliEMCALDigit", 1000);
 663     fDigitsArr->SetOwner(1);
 664   }
 665
 666   // then setup clusterizer
 667   delete fClusterizer;
 668   if     (fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerv1)
 669     fClusterizer = new AliEMCALClusterizerv1(geometry);
 670   else if (fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerNxN) {
 671    AliEMCALClusterizerNxN *clusterizer = new AliEMCALClusterizerNxN(geometry);
 672    clusterizer->SetNRowDiff(fRecParam->GetNRowDiff());
 673    clusterizer->SetNColDiff(fRecParam->GetNColDiff());
 674     fClusterizer = clusterizer;
 675   } else if (fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerv2)
 676     fClusterizer = new AliEMCALClusterizerv2(geometry);
 677   else if (fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerFW){
 678     AliEMCALClusterizerFixedWindow *clusterizer = new AliEMCALClusterizerFixedWindow(geometry);
 679     clusterizer->SetNphi(fNPhi);
 680     clusterizer->SetNeta(fNEta);
 681     clusterizer->SetShiftPhi(fShiftPhi);
 682     clusterizer->SetShiftEta(fShiftEta);
 683     clusterizer->SetTRUshift(fTRUShift);
 684     fClusterizer = clusterizer;
 685   }
 686   else{
 687     AliFatal(Form("Clusterizer < %d > not available", fRecParam->GetClusterizerFlag()));
 688   }
 689   fClusterizer->InitParameters(fRecParam);
 690
 691   if ((!fCalibData&&fLoadCalib) || (!fPedestalData&&fLoadPed)) {
 692     AliCDBManager *cdb = AliCDBManager::Instance();
 693     if (!cdb->IsDefaultStorageSet() && !fOCDBpath.IsNull())
 694       cdb->SetDefaultStorage(fOCDBpath);
 695     if (fRun!=cdb->GetRun())
 696       cdb->SetRun(fRun);
 697   }
 698   if (!fCalibData&&fLoadCalib&&fRun>0) {
 699     AliCDBEntry *entry = static_cast<AliCDBEntry*>(AliCDBManager::Instance()->Get("EMCAL/Calib/Data"));
 700     if (entry)
 701       fCalibData =  static_cast<AliEMCALCalibData*>(entry->GetObject());
 702     if (!fCalibData)
 703       AliFatal("Calibration parameters not found in CDB!");
 704   }
 705   if (!fPedestalData&&fLoadPed&&fRun>0) {
 706     AliCDBEntry *entry = static_cast<AliCDBEntry*>(AliCDBManager::Instance()->Get("EMCAL/Calib/Pedestals"));
 707     if (entry)
 708       fPedestalData =  static_cast<AliCaloCalibPedestal*>(entry->GetObject());
 709   }
 710   if (fCalibData) {
 711     fClusterizer->SetInputCalibrated(kFALSE);
 712     fClusterizer->SetCalibrationParameters(fCalibData);
 713   } else {
 714     fClusterizer->SetInputCalibrated(kTRUE);
 715   }
 716   fClusterizer->SetCaloCalibPedestal(fPedestalData);
 717   fClusterizer->SetJustClusters(kTRUE);
 718   fClusterizer->SetDigitsArr(fDigitsArr);
 719   fClusterizer->SetOutput(0);
 720   fClusterArr = const_cast<TObjArray *>(fClusterizer->GetRecPoints());
 721 }