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