macro dir
[u/mrichter/AliRoot.git] / 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 }