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