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